Skip to content

Commit

Permalink
refactor geometry optimization (part 2) (#982)
Browse files Browse the repository at this point in the history 
* refactor sub ancopt

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* refactor sub relax

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* refactor sub predchange

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* refactor sub axis2

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* refactor ANC type-bound procedures

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* refactor sub detrotra8

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* refactor sub qsort

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

* add comments to optimization parameters

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>

---------

Signed-off-by: albert <92109627+Albkat@users.noreply.github.com>
  • Loading branch information
@Albkat
Albkat authored May 15, 2024
1 parent 96e6a25 commit 67a6b3b
Show file tree
Hide file tree
Showing 9 changed files with 871 additions and 440 deletions.
198 changes: 115 additions & 83 deletions src/axis_trafo.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
123 changes: 66 additions & 57 deletions src/detrotra.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 1 addition & 1 deletion src/main/property.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
Loading

0 comments on commit 67a6b3b

Please sign in to comment.