diff --git a/src/optim/neb/neb_drv.F b/src/optim/neb/neb_drv.F index e50c1f5272..67a3e3eb3e 100644 --- a/src/optim/neb/neb_drv.F +++ b/src/optim/neb/neb_drv.F @@ -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: ||,=",sum2,sum + if (sum2.gt.1.0d0) alpha = 0.3d0 + if (oprint) write(*,*) "neb: ||,=",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: ||,=",sum2,sum end if @@ -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 @@ -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. diff --git a/src/optim/neb/neb_utils.F b/src/optim/neb/neb_utils.F index 8024aa622a..443e54d8a2 100644 --- a/src/optim/neb/neb_utils.F +++ b/src/optim/neb/neb_utils.F @@ -341,16 +341,16 @@ 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 @@ -358,7 +358,7 @@ subroutine neb_tangent(nbeads,nion,c,e,t,amat) 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 @@ -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 @@ -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 @@ -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))