refactor geometry optimization (part 2)

src/axis_trafo.f90

@@ -112,96 +112,128 @@ subroutine axis(numat,nat,xyz,aa,bb,cc)
     return
   end subroutine axis
 
-  subroutine axis2(numat,nat,xyz,aa,bb,cc,avmom,sumw)
-    use xtb_splitparam
-    implicit double precision (a-h,o-z)
-    dimension xyz(3,numat)
-    integer nat(numat)
 
-    PARAMETER (BOHR=0.52917726)
-    dimension t(6), rot(3), xyzmom(3), eig(3), evec(3,3)
-    dimension x(numat),y(numat),z(numat),coord(3,numat)
-    data t /6*0.d0/
-    !************************************************************************
-    !*     const1 =  10**40/(n*a*a)
-    !*               n = avergadro's number
-    !*               a = cm in an angstrom
-    !*               10**40 is to allow units to be 10**(-40)gram-cm**2
-    !*
-    !************************************************************************
-    const1 = 1.66053d0
-    !************************************************************************
-    !*
-    !*     const2 = conversion factor from angstrom-amu to cm**(-1)
-    !*
-    !*            = (planck's constant*n*10**16)/(8*pi*pi*c)
-    !*            = 6.62618*10**(-27)[erg-sec]*6.02205*10**23*10**16/
-    !*              (8*(3.1415926535)**2*2.997925*10**10[cm/sec])
-    !*
-    !************************************************************************
-    const2=16.8576522d0
-
-    sumw=1.d-20
-    sumwx=0.d0
-    sumwy=0.d0
-    sumwz=0.d0
 
-    coord(1:3,1:numat)=xyz(1:3,1:numat)*bohr
+subroutine axis2(n,xyz,aa,bb,cc,avmom,sumw)
+
+   use xtb_splitparam, only : atmass
+   use xtb_mctc_convert, only : autoaa
+
+   implicit double precision (a-h,o-z)
+
+   !> number of atoms 
+   integer, intent(in) :: n
+
+   !> cartesian coordinates
+   real(wp), intent(in) :: xyz(:,:)
+
+   !> rotational constants
+   real(wp), intent(out) :: aa,bb,cc
+
+   !> average moment of inertia
+   real(wp), intent(out) :: avmom
+
+   !> sum of atomic masses
+   real(wp), intent(out) :: sumw
+
+   !> const1 =  10**40/(n*a*a)
+   !>            n = avergadro's number
+   !>            a = cm in an angstrom
+   !>            10**40 is to allow units to be 10**(-40)gram-cm**2
+   real(wp), parameter :: const1 = 1.66053_wp
+
+   !> const2 = conversion factor from angstrom-amu to cm**(-1)
+   !>        = (planck's constant*n*10**16)/(8*pi*pi*c)
+   !>        = 6.62618*10**(-27)[erg-sec]*6.02205*10**23*10**16/
+   !>         (8*(3.1415926535)**2*2.997925*10**10[cm/sec])
+   real(wp), parameter :: const2 = 16.8576522_wp
+
+   !> temporary variables
+   real(wp) :: xsum,eps,ams
+
+   !> weighted sum of coordinates
+   real(wp) :: sumwx,sumwy,sumwz
+
+   !> loop variables
+   integer :: i
+
+   !> temporary arrays
+   real(wp) :: t(6), rot(3), xyzmom(3), eig(3), evec(3,3)
+   real(wp) :: x(n),y(n),z(n)
+
+   !> Cartesian coordinates in Bohr
+   real(wp) :: coord(3,n)  
+
+   t(:)  = 0.0_wp 
+   sumw  = 0.0_wp
+   sumwx = 0.0_wp
+   sumwy = 0.0_wp
+   sumwz = 0.0_wp
+
+   ! convert Bohr to Angstrom ! 
+   coord(:,:) = xyz(:,:) * autoaa
+
+   ! cma !
+   do i=1,n
+      sumw=sumw+atmass(i)
+      sumwx=sumwx+atmass(i)*coord(1,i)
+      sumwy=sumwy+atmass(i)*coord(2,i)
+      sumwz=sumwz+atmass(i)*coord(3,i)
+   enddo
+
+   sumwx=sumwx/sumw
+   sumwy=sumwy/sumw
+   sumwz=sumwz/sumw
+
+   do i=1,n
+      x(i)=coord(1,i)-sumwx
+      y(i)=coord(2,i)-sumwy
+      z(i)=coord(3,i)-sumwz
+   enddo
+
+   !************************************************************************
+   !*    matrix for moments of inertia is of form
+   !*
+   !*           |   y**2+z**2                         |
+   !*           |    -y*x       z**2+x**2             | -i =0
+   !*           |    -z*x        -z*y       x**2+y**2 |
+   !*
+   !************************************************************************
+   do i=1,6
+      t(i)=dble(i)*1.0d-10
+   enddo
+
+   do i=1,n
+      t(1)=t(1)+atmass(i)*(y(i)**2+z(i)**2)
+      t(2)=t(2)-atmass(i)*x(i)*y(i)
+      t(3)=t(3)+atmass(i)*(z(i)**2+x(i)**2)
+      t(4)=t(4)-atmass(i)*z(i)*x(i)
+      t(5)=t(5)-atmass(i)*y(i)*z(i)
+      t(6)=t(6)+atmass(i)*(x(i)**2+y(i)**2)
+   enddo
+
+   call rsp(t,3,3,eig,evec)
+
+   do i=1,3
+
+      if(eig(i).lt.3.d-4) then
+         eig(i)=0.d0
+         rot(i)=0.d0
+      else
+         rot(i)=2.9979245d+4*const2/eig(i)
+      endif
 
-    do i=1,numat
-       sumw=sumw+atmass(i)
-       sumwx=sumwx+atmass(i)*coord(1,i)
-       sumwy=sumwy+atmass(i)*coord(2,i)
-       sumwz=sumwz+atmass(i)*coord(3,i)
-    enddo
+      xyzmom(i)=eig(i)*const1
 
-    sumwx=sumwx/sumw
-    sumwy=sumwy/sumw
-    sumwz=sumwz/sumw
-    f=1.0d0/bohr
-    do i=1,numat
-       x(i)=coord(1,i)-sumwx
-       y(i)=coord(2,i)-sumwy
-       z(i)=coord(3,i)-sumwz
-    enddo
+   enddo
 
-    !************************************************************************
-    !*    matrix for moments of inertia is of form
-    !*
-    !*           |   y**2+z**2                         |
-    !*           |    -y*x       z**2+x**2             | -i =0
-    !*           |    -z*x        -z*y       x**2+y**2 |
-    !*
-    !************************************************************************
-    do i=1,6
-       t(i)=dble(i)*1.0d-10
-    enddo
-    do i=1,numat
-       t(1)=t(1)+atmass(i)*(y(i)**2+z(i)**2)
-       t(2)=t(2)-atmass(i)*x(i)*y(i)
-       t(3)=t(3)+atmass(i)*(z(i)**2+x(i)**2)
-       t(4)=t(4)-atmass(i)*z(i)*x(i)
-       t(5)=t(5)-atmass(i)*y(i)*z(i)
-       t(6)=t(6)+atmass(i)*(x(i)**2+y(i)**2)
-    enddo
-    call rsp(t,3,3,eig,evec)
-    do i=1,3
-       if(eig(i).lt.3.d-4) then
-          eig(i)=0.d0
-          rot(i)=0.d0
-       else
-          rot(i)=2.9979245d+4*const2/eig(i)
-       endif
-       xyzmom(i)=eig(i)*const1
-    enddo
+   aa=rot(3)/2.9979245d+4
+   bb=rot(2)/2.9979245d+4
+   cc=rot(1)/2.9979245d+4
 
-    aa=rot(3)/2.9979245d+4
-    bb=rot(2)/2.9979245d+4
-    cc=rot(1)/2.9979245d+4
-    avmom=1.d-47*(xyzmom(1)+xyzmom(2)+xyzmom(3))/3.
+   avmom=1.d-47*(xyzmom(1)+xyzmom(2)+xyzmom(3))/3.
 
-    return
-  end subroutine axis2
+end subroutine axis2
 
   subroutine axisvec(numat,nat,xyz,aa,bb,cc,evec)
     use xtb_splitparam

src/detrotra.f90

@@ -77,62 +77,71 @@ subroutine detrotra4(linear,mol,h,eig)
     enddo
 
   end subroutine detrotra4
-
-  subroutine detrotra8(linear,n,xyz,h,eig)
-    use xtb_mctc_accuracy, only : wp
-    use xtb_type_molecule
-    implicit none
-    integer, intent(in)     :: n
-    real(wp), intent(in)    :: xyz(3,n)    ! values from projected Lindh diag
-    real(wp), intent(in)    :: h(3*n,3*n)    ! values from projected Lindh diag
-    real(wp), intent(inout) :: eig(3*n)          ! eigenvectors from projected Lindh diag
-    logical, intent(in)     :: linear
-
-    integer               :: i,j,k,kk,ii,nn,n3,nend
-    integer,allocatable   :: ind(:)
-    real(wp), allocatable :: tmp(:,:)
-    real(wp), allocatable :: e(:)
-    real(wp)              :: a0,b0,c0
-
-    n3 = 3*n
-
-    allocate(tmp(3,n),e(n3),ind(n3))
-
-    nn = 0
-    do ii=1, n3
-       if(eig(ii).gt.0.05) cycle  ! only lowest checked
-
-       kk=0
-       do j=1,n
-          do k=1,3
-             kk=kk+1
-             tmp(k,j)=xyz(k,j)+h(kk,ii) ! distort along mode ii
-          enddo
-       enddo
-
-       c0=0            ! compared all interatomic distances of original and distortet geom.
-       do i=2,n
-          do j=1,i-1
-             a0 = sqrt((xyz(1,i)-xyz(1,j))**2+(xyz(2,i)-xyz(2,j))**2+(xyz(3,i)-xyz(3,j))**2)
-             b0 = sqrt((tmp(1,i)-tmp(1,j))**2+(tmp(2,i)-tmp(2,j))**2+(tmp(3,i)-tmp(3,j))**2)
-             c0 = c0+(a0-b0)**2   
-          enddo
-       enddo
-       nn = nn + 1
-       e(nn)=sqrt(c0/n)*abs(eig(ii))  ! weight by Lindh eigenvalue
-       ind(nn)=nn
-
-    enddo
-
-    call qsort(e, 1, nn, ind)   ! sort
-
-    nend = 6
-    if (linear) nend = 5
-
-    do i=1,nend
-       eig(ind(i)) = 0.0  ! identifier for rot/tra
-    enddo
-
-  end subroutine detrotra8
+
+!> determine rotational and translational modes 
+subroutine detrotra8(linear,n,xyz,h,eig)
+
+   use xtb_mctc_accuracy, only : wp
+   use xtb_type_molecule
+   implicit none
+
+   integer, intent(in)     :: n
+   real(wp), intent(in)    :: xyz(3,n)    ! values from projected Lindh diag
+   real(wp), intent(in)    :: h(3*n,3*n)    ! values from projected Lindh diag
+   real(wp), intent(inout) :: eig(3*n)          ! eigenvalues from projected Lindh diag
+   logical, intent(in)     :: linear
+
+   integer               :: i,j,k,kk,ii,nn,n3,nend
+   integer,allocatable   :: ind(:)
+   real(wp), allocatable :: tmpxyz(:,:)
+   real(wp), allocatable :: e(:)
+   real(wp)              :: a0,b0,c0
+
+   n3 = 3*n
+
+   allocate(tmpxyz(3,n),e(n3),ind(n3))
+
+   nn = 0
+   do ii=1, n3
+
+      if(eig(ii).gt.0.05) cycle  ! check only low-lying modes 
+
+      ! distort initial geometry along ii-th mode !
+      kk=0
+      do j=1,n
+         do k=1,3
+            kk=kk+1
+            tmpxyz(k,j)=xyz(k,j)+h(kk,ii)
+         enddo
+      enddo
+
+      ! compare all interatomic distances of original and distortet geom. !
+      c0=0           
+      do i=2,n
+         do j=1,i-1
+            a0 = sqrt((xyz(1,i)-xyz(1,j))**2+(xyz(2,i)-xyz(2,j))**2+(xyz(3,i)-xyz(3,j))**2)
+            b0 = sqrt((tmpxyz(1,i)-tmpxyz(1,j))**2+(tmpxyz(2,i)-tmpxyz(2,j))**2+(tmpxyz(3,i)-tmpxyz(3,j))**2)
+            c0 = c0+(a0-b0)**2 ! sum of squared differences  
+         enddo
+      enddo
+
+      nn = nn + 1 
+      e(nn)=sqrt(c0/n)*abs(eig(ii))  ! weight by Lindh eigenvalue
+      ind(nn)=nn
+
+   enddo
+
+   ! sort in ascending order !
+   call qsort(e,1,nn,ind)   
+
+   nend = 6
+   if (linear) nend = 5
+
+   ! set lowest eigenvalues to zero !
+   do i=1,nend
+      eig(ind(i)) = 0.0  
+   enddo
+
+end subroutine detrotra8
 
 end module xtb_detrotra

src/main/property.F90

@@ -1239,7 +1239,7 @@ subroutine print_thermo(iunit, nat, nvib_in, at, xyz, freq, etot, htot, gtot, ni
       nvib = 0
       nimag = 0
 
-      call axis2(nat, at, xyz, aa, bb, cc, avmom, wt)
+   call axis2(nat,xyz,aa,bb,cc,avmom,wt)
 
       nvib_theo = 3 * nat - 6
       if (cc < 1.d-10) linear = .true.