Skip to content

Commit

Permalink
Moved to using force-converged GRP in RunBreakdown
Browse files Browse the repository at this point in the history
  • Loading branch information
@joegilkes
joegilkes committed Apr 13, 2023
1 parent 52cbf14 commit 2616f37
Show file tree
Hide file tree
Showing 2 changed files with 118 additions and 13 deletions.
2 changes: 1 addition & 1 deletion src/pathfinder.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2638,7 +2638,7 @@ Subroutine GraphsToCoords_BD(wcx, printflag, printfile, &
! species bonding at long distances, leading to broken GRP geometries.
! call OptimizeGRP_DoubleEnded(wcx(2), wcx(1), success, gdsrestspring, nbstrength, nbrange, &
! kradius, ngdsrelax, gdsdtrelax)
call OptimizeGRP(wcx(2), success, gdsrestspring, nbstrength, nbrange, &
call OptimizeGRPForceConv(wcx(2), success, gdsrestspring, nbstrength, nbrange, &
& kradius, ngdsrelax, gdsdtrelax)

if (.not. success) then
Expand Down
129 changes: 117 additions & 12 deletions src/structure.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2769,7 +2769,18 @@ Subroutine OptimizeGRP(cx, success, gdsrestspring, nbstrength, nbrange, kradius,
real(8) :: gdsrestspring, nbstrength, nbrange, kradius, gdsdtrelax,sum,rmax
logical :: success, debug

debug = .false.
debug = .true.
if (debug) open(21, file='ori.xyz', status='unknown', position='append')
if (debug) then
write(21, *) cx%na
write(21, *) 'FORCE? = ', norm2(reshape(cx%dvdr(1:3,1:cx%na), (/cx%na*3/)))
do i = 1, cx%na
write(21, '("'//cx%Atomlabel(i)//'", 3(X, F15.7))') cx%r(1:3, i)*bohr_to_ang
enddo
close(21)
endif
if (debug) open(22,file='banana.xyz',status='unknown',position='append')

! Make a copy of the current structure - cx%graph(:,:) should already contain
! the target graph.
!
Expand All @@ -2778,7 +2789,6 @@ Subroutine OptimizeGRP(cx, success, gdsrestspring, nbstrength, nbrange, kradius,
cc2 = 0 ; cc1 = 0
cx%dvdr(:, :) = 0.0D0
call GraphConstraints(cx, gdsrestspring, nbstrength, nbrange, kradius)
! if (debug) open(22,file='banana.xyz',status='unknown')
do it = 1, ngdsrelax
idof = 0
sum = 0.d0
Expand Down Expand Up @@ -2821,16 +2831,16 @@ Subroutine OptimizeGRP(cx, success, gdsrestspring, nbstrength, nbrange, kradius,
!enddo
!flush(22)
endif
! if (debug) then
! write(22,*) cx%na
! write(22,*) 'FORCE? = ', norm2(reshape(cx%dvdr(1:3,1:cx%na),(/cx%na*3/)))
! do i = 1, cx%na
! write(22,'("'//cxtmp2%atomlabel(i)//'",3(X,F15.7))') cx%r(1:3,i)*bohr_to_ang
! enddo
! flush(22)
!print*, 'FORCE? = ', norm2(reshape(cx%dvdr(1:3,1:cx%na),(/cx%na*3/)))
! endif
call GraphConstraints(cx, gdsrestspring, nbstrength, nbrange, kradius)
if (debug) then
write(22, *) cx%na
write(22, *) 'FORCE? = ', norm2(reshape(cx%dvdr(1:3, 1:cx%na), (/cx%na*3/)))
do i = 1, cx%na
write(22, '("'//cx%Atomlabel(i)//'", 3(X, F15.7))') cx%r(1:3, i)*bohr_to_ang
enddo
flush(22)
!print*, 'FORCE? = ', norm2(reshape(cx%dvdr(1:3,1:cx%na),(/cx%na*3/)))
endif
!print*, 'CCX = ', cc2
if (cc2 > 3 ) then
success = .true.
Expand All @@ -2853,11 +2863,106 @@ Subroutine OptimizeGRP(cx, success, gdsrestspring, nbstrength, nbrange, kradius,
!enddo
!if (isum == 0 ) cc2 = cc2 + 1
!stOP
!close(22)
if (debug) close(22)
return
end Subroutine OptimizeGRP


!
!*************************************************************************
!> OptimizeGRPForceConv
!!
!! Minimises the molecules according to the Graph restraining potential.
!!
!! Unlike the original OptimizeGRP routine, uses the force convergence
!! criteria of OptimizeGRPDoubleEnded to determine if optimisations
!! are successful.
!!
!! - cx: The input chemical structure
!! - success: weather we satisfied the contraints after the monimization or not
!! - rest: parameters for the GRP and minimization
!!
!*************************************************************************
!
Subroutine OptimizeGRPForceConv(cx, success, gdsrestspring, nbstrength, nbrange, kradius, ngdsrelax, gdsdtrelax)
implicit none
type(cxs) :: cx, cxtmp1, cxtmp2
integer :: i, j, k, l, n, m, it, cc1, cc2, idof, isum, ngdsrelax, iact
real(8) :: gdsrestspring, nbstrength, nbrange, kradius, gdsdtrelax,sum,rmax
logical :: success, debug

debug = .false.
if (debug) open(21, file='ori.xyz', status='unknown', position='append')
if (debug) then
write(21, *) cx%na
write(21, *) 'FORCE? = ', norm2(reshape(cx%dvdr(1:3,1:cx%na), (/cx%na*3/)))
do i = 1, cx%na
write(21, '("'//cx%Atomlabel(i)//'", 3(X, F15.7))') cx%r(1:3, i)*bohr_to_ang
enddo
close(21)
endif
if (debug) open(22,file='banana.xyz',status='unknown',position='append')

! SD optimisation
success = .false.
outer: do it = 0, ngdsrelax
idof = 0
sum = 0.d0
rmax = -1d6
iact = 0

if (it > 0) then

do i = 1, cx%na
if (.not. cx%fixedatom(i)) Then
do k = 1, 3
idof = idof + 1
if (.not. cx%fixeddof(idof)) then
iact = iact + 1
cx%r(k, i) = cx%r(k, i) - gdsdtrelax * cx%dvdr(k, i)
sum = sum + cx%dvdr(k, i)**2
if (abs( cx%dvdr(k, i)) > rmax) rmax = abs(cx%dvdr(k, i))
endif
enddo
else
idof = idof + 3
endif
enddo

! Check convergence based on RMS forces and maximum force.
sum = dsqrt(sum / dble(idof) )

if (debug) then
print *, 'SUM = ', sum, ' MAX = ', rmax
endif

if (sum < GRPMINTHRESH .and. rmax < GRPMAXTHRESH) then
success = .true.
exit outer
endif
endif

if (debug) then
write(22, *) cx%na
! write(22, *) 'FORCE? = ', norm2(reshape(cx%dvdr(1:3, 1:cx%na), (/cx%na*3/)))
write(22, '("Properties=species:S:1:pos:R:3:forces:R:3")')
do i = 1, cx%na
write(22, '("'//cx%Atomlabel(i)//'", 6(X, F15.7))') cx%r(1:3, i)*bohr_to_ang, -cx%dvdr(1:3, i)*bohr_to_ang
enddo
flush(22)
endif

! Update derivatives.
cx%dvdr(:, :) = 0.0D0
call GraphConstraints(cx, gdsrestspring, nbstrength, nbrange, kradius)

enddo outer
if (debug) close(22)

return
end Subroutine OptimizeGRPForceConv


!
!*************************************************************************
!> OptimizeGRP2
Expand Down

0 comments on commit 2616f37

Please sign in to comment.