Skip to content

Commit

Permalink
NEB changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@dmejiar
dmejiar committed Jan 17, 2024
1 parent bc14faa commit ed43df1
Show file tree
Hide file tree
Showing 2 changed files with 53 additions and 39 deletions.
68 changes: 37 additions & 31 deletions src/optim/neb/neb_drv.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -450,24 +450,26 @@ logical function neb(rtdb)
finishedstep = .false.
alpha = time_step
jj = 0
do while ((.not.finishedstep).and.(jj.lt.9))
jj = jj + 1
c do while ((.not.finishedstep).and.(jj.lt.9))
c jj = jj + 1
sum = ddot(ng,dbl_mb(g0(1)),1,dbl_mb(s(1)),1)
sum2 = dsqrt(ddot(ng,dbl_mb(s(1)),1,dbl_mb(s(1)),1))
> /dble(nbeads)
if (oprint) write(*,*) "neb: |<s|s>|,<g|s>=",sum2,sum
if (sum2.gt.1.0d0) alpha = 0.3d0

if (oprint) write(*,*) "neb: |<s|s>|,<g|s>=",sum2,sum
c
if ((sum.le.0.0d0).or.
> (sum0.gt.sum0_old).or.
> (sum2.gt.1.0d0)) then
c > (sum0.gt.sum0_old).or.
> ((sum2.gt.1.0d0).and.(itm.gt.2))) then
call dcopy(ng,dbl_mb(g0(1)),1,dbl_mb(s(1)),1)
itm = 0
if (oprint) write(*,*) "neb: s=g and itm reset to 0"
if (oprint) write(*,*) "neb: sum,sum0,sum0_old=",
> sum,sum0,sum0_old
alpha = time_step
sum = ddot(ng,dbl_mb(g0(1)),1,dbl_mb(s(1)),1)
sum2 = dsqrt(ddot(ng,dbl_mb(s(1)),1,dbl_mb(s(1)),1))
> /dble(nbeads)
if (sum2.gt.1.0d0) alpha = 0.3d0
if (oprint) write(*,*) "neb: |<s|s>|,<g|s>=",sum2,sum
end if

Expand Down Expand Up @@ -524,35 +526,35 @@ logical function neb(rtdb)
c > x0,x1,x2,f0,f1,f2
c end if

finishedstep = (sum0.le.sum0_old)
if (.not.finishedstep) then
alpha = 0.50d0*alpha
manyjj = manyjj+1
goodjj = 0
if (manyjj.gt.2) then
time_step = 0.5d0*time_step
manyjj = 0
if (oprint)
> write(*,*) "neb: reducing timestep=",time_step
end if
end if
if (finishedstep.and.(jj.eq.1)) then
goodjj = goodjj + 1
if (goodjj.gt.10) then
time_step = 2.0d0*time_step
if (time_step.gt.time_step0) time_step=time_step0
goodjj = 0
end if
end if

c finishedstep = (sum0.le.sum0_old)
c if (.not.finishedstep) then
c alpha = 0.50d0*alpha
c manyjj = manyjj+1
c goodjj = 0
c if (manyjj.gt.2) then
c time_step = 0.5d0*time_step
c manyjj = 0
c if (oprint)
c > write(*,*) "neb: reducing timestep=",time_step
c end if
c end if
c if (finishedstep.and.(jj.eq.1)) then
c goodjj = goodjj + 1
c if (goodjj.gt.10) then
c time_step = 2.0d0*time_step
c if (time_step.gt.time_step0) time_step=time_step0
c goodjj = 0
c end if
c end if
c
if(oprint)
> write(*,*) "neb: sum0a,sum0b,sum0,sum0_old=",
> sum0a,sum0b,sum0,sum0_old,jj,finishedstep,alpha

c if (finishedstep) sum0_old = sum0

end do
if (sum0b.ge.sum0_old) time_step = 0.5d0*time_step
c end do
if (sum0b.ge.sum0_old) time_step = time_step
sum0_old = sum0


Expand Down Expand Up @@ -641,12 +643,16 @@ logical function neb(rtdb)
> dbl_mb(c1(1)),
> Gmax,Grms,Xmax,Xrms)

stalled = ((Gmax.ge.gmax0).or.(Grms.ge.grms0)).and.(it.gt.1)
stalled=((Gmax.ge.1.3d0*gmax0).or.(Grms.ge.1.3d0**grms0)).and.
> (it.gt.1)

!*** neb algorithm 3 - switch between fixed point and damped Verlet ***
if (neb_algorithm.ge.3) then

!*** switch to fixed point if Gmax less than 0.5 ***
if (ga_nodeid().eq.0)write(*,*)gmax.lt.0.5d0,gmax.lt.gmax0,
> grms.lt.grms0

if ((Gmax.lt.0.5d0).and.
> (Gmax.lt.Gmax0).and.
> (Grms.lt.Grms0).and.
Expand Down
24 changes: 16 additions & 8 deletions src/optim/neb/neb_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -341,24 +341,24 @@ subroutine neb_tangent(nbeads,nion,c,e,t,amat)
c > c(index),1,
c > t(rp), 1)

if ( ((e(i+1)-e(i)).gt.1.0e-2).and.
> ((e(i)-e(i-1)).gt.1.0d-2) ) then
if ( ((e(i+1)-e(i)).gt.0d0).and.
> ((e(i)-e(i-1)).gt.0d0) ) then
call dcopy(3*nion,t(rp),1,t(index),1)
else if ( ((e(i-1)-e(i)).gt.1.0d-2) .and.
> ((e(i)-e(i+1)).gt.1.0d-2) ) then
else if ( ((e(i-1)-e(i)).gt.0d0) .and.
> ((e(i)-e(i+1)).gt.0d0) ) then
call dcopy(3*nion,t(rm),1,t(index),1)
else

if ( (dabs(e(i+1)-e(i))-dabs(e(i-1)-e(i)))
> .gt.1.0d-3) then
> .gt.0d0) then
dVmax = dabs(e(i+1)-e(i))
dVmin = dabs(e(i-1)-e(i))
else
dVmax = dabs(e(i-1)-e(i))
dVmin = dabs(e(i+1)-e(i))
end if

if ((e(i+1)-e(i-1)).gt.1.0d-2) then
if ((e(i+1)-e(i-1)).gt.0d0) then
call dscal(3*nion,dVmax,t(rp),1)
call dscal(3*nion,dVmin,t(rm),1)
else
Expand Down Expand Up @@ -1145,7 +1145,7 @@ subroutine neb_verlet_update(ng,c0,c1,v1,dti,g1)
* *** NEED: CONSTRAINED DYNAMICS (FROZEN ATOMS)
do i=1,ng
v1(i) = c1(i)-c0(i)
if (v1(i)*g1(i).gt.0.0d0) v1(i) = 0.0d0
c if (v1(i)*g1(i).gt.0.0d0) v1(i) = 0.0d0
end do
call dcopy(ng,c1(1),1,c0(1),1)
do i=1,ng
Expand Down Expand Up @@ -1278,7 +1278,7 @@ subroutine neb_lmbfgs(n,m,x,g,hg)
tmp = tmp - ddot(n,x(1,k+1),1,g(1,k), 1)
tmp = tmp - ddot(n,x(1,k), 1,g(1,k+1),1)
tmp = tmp + ddot(n,x(1,k), 1,g(1,k), 1)
if (dabs(tmp).gt.1.0d-9) then
if (dabs(tmp).gt.1.0d-12) then
rho(k) = 1.0d0/tmp
else
rho(k) = 0.0d0
Expand All @@ -1295,6 +1295,14 @@ subroutine neb_lmbfgs(n,m,x,g,hg)
call daxpy(n,( alpha(k)),g(1,k), 1,hg,1)
end do

if (m.gt.2) then
tmp = ddot(n,g(1,m-1),1,g(1,m-1),1) -
$ 2.0d0*ddot(n,g(1,m-1),1,g(1,m-2),1) +
$ ddot(n,g(1,m-2),1,g(1,m-2),1)
tmp = 1.0d0/rho(m-2)/tmp
call dscal(n,tmp,hg,1)
endif

do k = 1,(m-1)
beta(k) = rho(k)
> *(ddot(n,g(1,k+1),1,hg,1) - ddot(n,g(1,k),1,hg,1))
Expand Down

0 comments on commit ed43df1

Please sign in to comment.