Skip to content

[pull] master from crest-lab:master #19

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pull merged 1 commit into from
Apr 3, 2025
Merged

[pull] master from crest-lab:master #19

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 12 additions & 4 deletions src/algos/parallel.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,13 +127,15 @@ subroutine crest_sploop(env,nat,nall,at,xyz,eread)
do i = 1,T
do j = 1,env%calc%ncalculations
calculations(i)%calcs(j) = env%calc%calcs(j)
!>--- directories
!>--- directories and io preparation
ex = directory_exist(env%calc%calcs(j)%calcspace)
if (.not.ex) then
io = makedir(trim(env%calc%calcs(j)%calcspace))
end if
write (atmp,'(a,"_",i0)') sep,i
calculations(i)%calcs(j)%calcspace = env%calc%calcs(j)%calcspace//trim(atmp)
if(allocated(calculations(i)%calcs(j)%calcfile)) deallocate(calculations(i)%calcs(j)%calcfile)
if(allocated(calculations(i)%calcs(j)%systemcall)) deallocate(calculations(i)%calcs(j)%systemcall)
call calculations(i)%calcs(j)%printid(i,j)
end do
calculations(i)%pr_energies = .false.
Expand Down Expand Up @@ -307,7 +309,7 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc)
do i = 1,T
do j = 1,mycalc%ncalculations
calculations(i)%calcs(j) = mycalc%calcs(j)
!>--- directories
!>--- directories and io preparation
ex = directory_exist(mycalc%calcs(j)%calcspace)
if (.not.ex) then
io = makedir(trim(mycalc%calcs(j)%calcspace))
Expand All @@ -317,6 +319,8 @@ subroutine crest_oloop(env,nat,nall,at,xyz,eread,dump,customcalc)
endif
write (atmp,'(a,"_",i0)') sep,i
calculations(i)%calcs(j)%calcspace = mycalc%calcs(j)%calcspace//trim(atmp)
if(allocated(calculations(i)%calcs(j)%calcfile)) deallocate(calculations(i)%calcs(j)%calcfile)
if(allocated(calculations(i)%calcs(j)%systemcall)) deallocate(calculations(i)%calcs(j)%systemcall)
call calculations(i)%calcs(j)%printid(i,j)
end do
calculations(i)%pr_energies = .false.
Expand Down Expand Up @@ -564,13 +568,15 @@ subroutine crest_search_multimd(env,mol,mddats,nsim)
moltmps(i)%xyz = mol%xyz
do j = 1,env%calc%ncalculations
calculations(i)%calcs(j) = env%calc%calcs(j)
!>--- directories
!>--- directories and io preparation
ex = directory_exist(env%calc%calcs(j)%calcspace)
if (.not.ex) then
io = makedir(trim(env%calc%calcs(j)%calcspace))
end if
write (atmp,'(a,"_",i0)') sep,i
calculations(i)%calcs(j)%calcspace = env%calc%calcs(j)%calcspace//trim(atmp)
if(allocated(calculations(i)%calcs(j)%calcfile)) deallocate(calculations(i)%calcs(j)%calcfile)
if(allocated(calculations(i)%calcs(j)%systemcall)) deallocate(calculations(i)%calcs(j)%systemcall)
call calculations(i)%calcs(j)%printid(i,j)
end do
calculations(i)%pr_energies = .false.
Expand Down Expand Up @@ -866,13 +872,15 @@ subroutine crest_search_multimd2(env,mols,mddats,nsim)
do i = 1,T
do j = 1,env%calc%ncalculations
calculations(i)%calcs(j) = env%calc%calcs(j)
!>--- directories
!>--- directories and io preparation
ex = directory_exist(env%calc%calcs(j)%calcspace)
if (.not.ex) then
io = makedir(trim(env%calc%calcs(j)%calcspace))
end if
write (atmp,'(a,"_",i0)') sep,i
calculations(i)%calcs(j)%calcspace = env%calc%calcs(j)%calcspace//trim(atmp)
if(allocated(calculations(i)%calcs(j)%calcfile)) deallocate(calculations(i)%calcs(j)%calcfile)
if(allocated(calculations(i)%calcs(j)%systemcall)) deallocate(calculations(i)%calcs(j)%systemcall)
call calculations(i)%calcs(j)%printid(i,j)
end do
calculations(i)%pr_energies = .false.
Expand Down
181 changes: 113 additions & 68 deletions src/dynamics/shake_module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,6 @@ module shake_module
public :: init_shake
public :: do_shake


!========================================================================================!
!========================================================================================!
contains !> MODULE PROCEDURES START HERE
Expand All @@ -101,7 +100,7 @@ module shake_module
subroutine init_shake(nat,at,xyz,shk,pr)
!********************************************
!* subroutine init_shake
!* initialize SHAKE algorithm by documenting
!* initialize SHAKE algorithm by documenting
!* which pairs of atoms to constrain
!********************************************
implicit none
Expand All @@ -125,29 +124,29 @@ subroutine init_shake(nat,at,xyz,shk,pr)
integer,allocatable :: cons2(:,:)

!>--- if shake was already set up, return
if(shk%initialized) return
if (shk%initialized) return

!>--- reset counter
nconsu = 0
n2 = 0
checkfreeze = associated(shk%freezeptr)
!> gfortran workaround
if(checkfreeze)then !> gfortran check
if (checkfreeze) then !> gfortran check
nfrz = size(shk%freezeptr,1)
else
nfrz = 0
endif
nfrz = 0
end if

!>--- count user-defined bonds
if ((shk%nusr > 0) .and. allocated(shk%conslistu)) then
if ((shk%nusr > 0).and.allocated(shk%conslistu)) then
do i = 1,shk%nusr
nconsu = nconsu + 1
nconsu = nconsu+1
end do
end if

!>--- constrain all X-H only
if (shk%shake_mode == 1) then
ij = nat * (nat + 1) / 2
ij = nat*(nat+1)/2
allocate (cons2(2,ij),source=0)
do i = 1,nat
if (at(i) .eq. 1) then
Expand All @@ -156,11 +155,11 @@ subroutine init_shake(nat,at,xyz,shk,pr)
if (j .ne. i) then

!> check if BOTH are frozen
if(checkfreeze.and.j<=nfrz.and.i<=nfrz)then
if( shk%freezeptr(i) .and. shk%freezeptr(j) ) cycle
endif
rij = (xyz(1,i) - xyz(1,j))**2 + (xyz(2,i) - xyz(2,j))**2 + (xyz(3,i) - xyz(3,j))**2
if (checkfreeze.and.j <= nfrz.and.i <= nfrz) then
if (shk%freezeptr(i).and.shk%freezeptr(j)) cycle
end if

rij = (xyz(1,i)-xyz(1,j))**2+(xyz(2,i)-xyz(2,j))**2+(xyz(3,i)-xyz(3,j))**2
if (rij .lt. minrij) then
minrij = rij
jmin = j
Expand All @@ -169,14 +168,14 @@ subroutine init_shake(nat,at,xyz,shk,pr)
end do
if (at(jmin) .eq. 1) then
if (jmin .gt. i) then
nconsu = nconsu + 1
n2 = n2 + 1
nconsu = nconsu+1
n2 = n2+1
cons2(1,n2) = i
cons2(2,n2) = jmin
end if
else
nconsu = nconsu + 1
n2 = n2 + 1
nconsu = nconsu+1
n2 = n2+1
cons2(1,n2) = i
cons2(2,n2) = jmin
end if
Expand All @@ -187,42 +186,42 @@ subroutine init_shake(nat,at,xyz,shk,pr)
!>--- SHAKE all bonds
if (shk%shake_mode == 2) then
if (allocated(shk%wbo)) then
ij = nat * (nat + 1) / 2
ij = nat*(nat+1)/2
allocate (cons2(2,ij),source=0)
allocate (list(ij),source=0)
do i = 1,nat
do j = 1,nat
if (i .eq. j) cycle

!>--- if both are frozen, no SHAKE!
if(checkfreeze.and.j<=nfrz.and.i<=nfrz)then
if( shk%freezeptr(i) .and. shk%freezeptr(j) ) cycle
endif
if (checkfreeze.and.j <= nfrz.and.i <= nfrz) then
if (shk%freezeptr(i).and.shk%freezeptr(j)) cycle
end if

rij = (xyz(1,i) - xyz(1,j))**2 + (xyz(2,i) - xyz(2,j))**2 + (xyz(3,i) - xyz(3,j))**2
rco = (atomicRad(at(j)) + atomicRad(at(i))) * autoaa
rij = (xyz(1,i)-xyz(1,j))**2+(xyz(2,i)-xyz(2,j))**2+(xyz(3,i)-xyz(3,j))**2
rco = (atomicRad(at(j))+atomicRad(at(i)))*autoaa
ij = shake_lin(i,j)
!>--- to consider?
rcut = (0.52917726_wp * sqrt(rij) .lt. 1.2_wp * rco) .and. (list(ij) .eq. 0)
rcut = (0.52917726_wp*sqrt(rij) .lt. 1.2_wp*rco).and.(list(ij) .eq. 0)
!>--- metal bond?
metalbond = metal(at(i)) .eq. 1 .or. metal(at(j)) .eq. 1
metalbond = metal(at(i)) .eq. 1.or.metal(at(j)) .eq. 1
!>--- WBO threshold. if WBO > thr, the bond is constrained
wthr = 0.5_wp
!>--- modify wthr, e.g. K...O in phosphates have WBO around 0.15
if (metalbond) wthr = 0.1_wp
!>--- check relevant
if (rcut .and. shk%wbo(i,j) .gt. wthr) then
if (rcut.and.shk%wbo(i,j) .gt. wthr) then
!>--- do not constrain M bonds except Li/Be
metalbond = metal(at(i)) .eq. 1 .or. metal(at(j)) .eq. 1
if (at(i) .eq. 3 .or. at(i) .eq. 4 .or. at(j) .eq. 3 .or. at(j) .eq. 4) &
& metalbond = .false.
metalbond = metal(at(i)) .eq. 1.or.metal(at(j)) .eq. 1
if (at(i) .eq. 3.or.at(i) .eq. 4.or.at(j) .eq. 3.or.at(j) .eq. 4) &
& metalbond = .false.
if (metalbond) then
if ((i .lt. j) .and. pr) &
& write (*,*) 'init_shake: metal bond ',i,j,'not constrained'
if ((i .lt. j).and.pr) &
& write (*,*) 'init_shake: metal bond ',i,j,'not constrained'
else
list(ij) = 1
nconsu = nconsu + 1
n2 = n2 + 1
nconsu = nconsu+1
n2 = n2+1
cons2(1,n2) = i
cons2(2,n2) = j
end if
Expand Down Expand Up @@ -250,7 +249,7 @@ subroutine init_shake(nat,at,xyz,shk,pr)

if (shk%nusr > 0) then
shk%conslist(1:2,1:shk%nusr) = shk%conslistu(1:2,1:shk%nusr)
j = shk%nusr + 1
j = shk%nusr+1
else
j = 1
end if
Expand All @@ -260,8 +259,8 @@ subroutine init_shake(nat,at,xyz,shk,pr)
do i = 1,shk%ncons
iat = shk%conslist(1,i)
jat = shk%conslist(2,i)
drij(1:3) = xyz(1:3,iat) - xyz(1:3,jat)
shk%distcons(i) = drij(1)**2 + drij(2)**2 + drij(3)**2
drij(1:3) = xyz(1:3,iat)-xyz(1:3,jat)
shk%distcons(i) = drij(1)**2+drij(2)**2+drij(3)**2
end do

if (allocated(cons2)) deallocate (cons2)
Expand Down Expand Up @@ -291,22 +290,31 @@ subroutine do_shake(nat,xyzo,xyz,velo,acc,mass,tstep,shk,pr,iostat)
real(wp) :: r,dev,dist,denom
real(wp) :: gcons,rmi,rmj
real(wp) :: tau1,tau2
integer :: jmaxdev
integer :: jmaxdev,nfrz
logical :: checkfreeze

conv = .false. !> assume not converged
io = 1 !> and therefore unsuccessfull
icyc = 0

if(.not.allocated(shk%xyzt)) allocate(shk%xyzt(3,nat))
if (.not.allocated(shk%xyzt)) allocate (shk%xyzt(3,nat))
shk%xyzt = xyz

tau1 = 1.0_wp / tstep
tau2 = tau1 * tau1
checkfreeze = associated(shk%freezeptr)
!> gfortran workaround
if (checkfreeze) then !> gfortran check
nfrz = size(shk%freezeptr,1)
else
nfrz = 0
end if

tau1 = 1.0_wp/tstep
tau2 = tau1*tau1

do i = 1,shk%ncons
iat = shk%conslist(1,i)
jat = shk%conslist(2,i)
shk%dro(1:3,i) = xyzo(1:3,iat) - xyzo(1:3,jat)
shk%dro(1:3,i) = xyzo(1:3,iat)-xyzo(1:3,jat)
end do

!>--- iterative SHAKE loop
Expand All @@ -316,10 +324,10 @@ subroutine do_shake(nat,xyzo,xyz,velo,acc,mass,tstep,shk,pr,iostat)
do i = 1,shk%ncons
iat = shk%conslist(1,i)
jat = shk%conslist(2,i)
shk%dr(1:3,i) = shk%xyzt(1:3,iat) - shk%xyzt(1:3,jat)
shk%dr(4,i) = shk%dr(1,i)**2 + shk%dr(2,i)**2 + shk%dr(3,i)**2
shk%dr(1:3,i) = shk%xyzt(1:3,iat)-shk%xyzt(1:3,jat)
shk%dr(4,i) = shk%dr(1,i)**2+shk%dr(2,i)**2+shk%dr(3,i)**2
dist = shk%distcons(i)
dev = abs(shk%dr(4,i) - dist) / dist
dev = abs(shk%dr(4,i)-dist)/dist
if (dev .gt. maxdev) then
maxdev = dev
jmaxdev = i
Expand All @@ -328,38 +336,75 @@ subroutine do_shake(nat,xyzo,xyz,velo,acc,mass,tstep,shk,pr,iostat)

if (maxdev .lt. shk%tolshake) conv = .true.

if (.not. conv) then
do i = 1,shk%ncons
iat = shk%conslist(1,i)
jat = shk%conslist(2,i)
dist = shk%distcons(i)
rmi = 1.0_wp / mass(iat)
rmj = 1.0_wp / mass(jat)
denom = 2.0_wp * (rmi + rmj) * (shk%dr(1,i) * shk%dro(1,i) + &
& shk%dr(2,i) * shk%dro(2,i) + &
& shk%dr(3,i) * shk%dro(3,i))
gcons = (dist - shk%dr(4,i)) / denom
shk%xyzt(1:3,iat) = shk%xyzt(1:3,iat) + rmi * gcons * shk%dro(1:3,i)
shk%xyzt(1:3,jat) = shk%xyzt(1:3,jat) - rmj * gcons * shk%dro(1:3,i)
end do
if (.not.conv) then
if (.not.checkfreeze.or.nfrz > 0) then !> Non-frozen case

do i = 1,shk%ncons
iat = shk%conslist(1,i)
jat = shk%conslist(2,i)
dist = shk%distcons(i)
rmi = 1.0_wp/mass(iat)
rmj = 1.0_wp/mass(jat)
denom = 2.0_wp*(rmi+rmj)*(shk%dr(1,i)*shk%dro(1,i)+ &
& shk%dr(2,i)*shk%dro(2,i)+ &
& shk%dr(3,i)*shk%dro(3,i))
gcons = (dist-shk%dr(4,i))/denom
shk%xyzt(1:3,iat) = shk%xyzt(1:3,iat)+rmi*gcons*shk%dro(1:3,i)
shk%xyzt(1:3,jat) = shk%xyzt(1:3,jat)-rmj*gcons*shk%dro(1:3,i)
end do

else !> Frozen case

do i = 1,shk%ncons
iat = shk%conslist(1,i)
jat = shk%conslist(2,i)
dist = shk%distcons(i)
rmi = 1.0_wp/mass(iat)
rmj = 1.0_wp/mass(jat)
denom = 2.0_wp*(rmi+rmj)*(shk%dr(1,i)*shk%dro(1,i)+ &
& shk%dr(2,i)*shk%dro(2,i)+ &
& shk%dr(3,i)*shk%dro(3,i))
gcons = (dist-shk%dr(4,i))/denom
!> update only non-frozen
if (.not.shk%freezeptr(iat)) then
shk%xyzt(1:3,iat) = shk%xyzt(1:3,iat)+rmi*gcons*shk%dro(1:3,i)
end if
if (.not.shk%freezeptr(jat)) then
shk%xyzt(1:3,jat) = shk%xyzt(1:3,jat)-rmj*gcons*shk%dro(1:3,i)
end if
end do
end if

end if
icyc = icyc + 1
if (.not. conv .and. icyc .le. shk%maxcyc) cycle
icyc = icyc+1
if (.not.conv.and.icyc .le. shk%maxcyc) cycle
exit
end do

if (conv) then
velo = velo + (shk%xyzt - xyz) * tau1
acc = acc + (shk%xyzt - xyz) * tau2
xyz = shk%xyzt
io = 0 !> successful termination

if (nfrz .eq. 0) then
velo = velo+(shk%xyzt-xyz)*tau1
acc = acc+(shk%xyzt-xyz)*tau2
xyz = shk%xyzt
io = 0 !> successful termination
else
do i = 1,nat
if (.not.shk%freezeptr(i)) then
velo(:,i) = velo(:,i)+(shk%xyzt(:,i)-xyz(:,i))*tau1
acc(:,i) = acc(:,i)+(shk%xyzt(:,i)-xyz(:,i))*tau2
xyz(:,i) = shk%xyzt(:,i)
end if
end do
end if
io = 0
else if (pr) then
write (*,*) 'SHAKE did not converge! maxdev=',maxdev
end if

if(present(iostat))then
if (present(iostat)) then
iostat = io
endif
end if

return
end subroutine do_shake
Expand All @@ -372,7 +417,7 @@ integer function shake_lin(i1,i2)
integer :: i1,i2,idum1,idum2
idum1 = max(i1,i2)
idum2 = min(i1,i2)
shake_lin = idum2 + idum1 * (idum1 - 1) / 2
shake_lin = idum2+idum1*(idum1-1)/2
return
end function shake_lin

Expand Down
Loading