1C $Id$ 2 Subroutine drdy_path2(iopt,delss,delsh,sold,smax,xsign) 3C 4C Calculate steepest descent path from saddle point 5C iopt - option for following MEP 6C =1, use Euler integrator 7C =2, use Page-McIver LQA algorithm 8C =3, use Page-McIver CLQA algorithm 9C =4, use Page-McIver CUBE algorithm 10C 11 Implicit none 12#include "errquit.fh" 13 Integer iopt 14 Double Precision delss,delsh,sold,smax,xsign 15 Integer i,ierr,j,nbar,nosc 16 Double Precision delsf,dotp,dxmag,shess,swf 17C 18 Double Precision eps 19c 20#include "drdyP.fh" 21C 22 data eps/1.d-6/ 23C 24 delsf = xsign*(s-sold) 25 shess = s 26 nosc = 0 27 nbar = (n3*(n3+1))/2 28 do while (s*xsign.lt.smax*xsign-dels*eps) 29 if (iopt.le.1) then 30 call drdy_euler (n3tm,n3,x,dx,dels) 31 elseif (iopt.eq.2) then 32 ierr = 0 33 call drdy_pmlqa (n3tm,n3,nbar,x,dx,f,hess,amass,dels, 34 * vec0,vec1,u0,scr1,vec2,ierr) 35 if (ierr.ne.0) then 36 write(fu6,*) ' error in pmlqa' 37 call errquit('error in pmlqa',555, UNKNOWN_ERR) 38 endif 39 elseif (iopt.eq.3) then 40 ierr = 0 41 call drdy_pmclqa (n3tm,n3,nbar,x,dx,f,fold,hess,amass, 42 & dels,delsf,vec0,vec1,u0,scr1,vec2,ierr) 43 if (ierr.ne.0) then 44 write(fu6,*) ' error in pmclqa' 45 call errquit('error in pmlqa',555, UNKNOWN_ERR) 46 endif 47 elseif (iopt.ge.4) then 48 call drdy_pmcube (n3tm,n3,x,dx,f,fold,amass,dels,delsf, 49 * vec0,vec1,vec2) 50 endif 51 call drdy_center (natom,n3,amass,x,scr1) 52 s = s + dels*xsign 53 delss = delss + dels 54 delsh = delsh + dels 55 if (delss.lt.delsv*(1.0d00-eps) .and. 56 * delsh.lt.delhss*(1.0d00-eps)) then 57 Call drdy_pot1 (x,v,dx,scr1,amass,natom,n3) 58 else 59 if (iopt.ge.3) then 60C Store previous hessian matrix in fold 61 do i = 1,n3 62 do j = 1,n3 63 fold(j,i) = f(j,i) 64 enddo 65 enddo 66 sold = shess 67 delsf = xsign*(s-sold) 68 shess = s 69 endif 70 Call drdy_pot2 (x,v,dx,f,hess,scr1,amass,natom,n3,n3tm, 71 * nbar) 72 delsh = 0.0d00 73 endif 74 write(fu6,610) s,v-vzero,(x(i),i=1,n3) 75 dxmag = 0.0d00 76 do i = 1,n3 77 dxmag = dxmag + dx(i)*dx(i) 78 enddo 79 dxmag = sqrt(dxmag) 80 write(fu6,611) dxmag,(dx(i)/dxmag,i=1,n3) 81C Check for oscillating gradient 82 dotp = 0.0d00 83 do i = 1,n3 84 dotp = dotp + dxold(i)*dx(i)/dxmag 85 dxold(i) = dx(i)/dxmag 86 enddo 87 if (dotp.lt. 0.0d00) then 88 nosc = nosc + 1 89 if (nosc .ge. 5) then 90 write (fu6,6000) 91 call errquit('The gradient is unstable',555, INPUT_ERR) 92 endif 93 else 94 nosc = 0 95 endif 96 if (delss.ge.delsv*(1.0d00-eps)) then 97 swf = s 98 if(lgs(38).gt.0) call drdy_wfiles(swf,0) 99 if (lgs(39).ne.0) then 100 call drdy_potsp(x,vspc,scr1,amass,natom,n3) 101 write (fu6,612) vspc-vzerosp 102 if(lgs(38).eq.4) call drdy_wfiles(swf,1) 103 endif 104 call drdy_rphwrt(7) 105 if(lgs(2).ne.0) then 106 ierr = 0 107 call drdy_projct(ierr) 108 if (ierr.eq.0) call drdy_fdiag(n3,nf(5),ierr,lgs(2)) 109 endif 110 delss = 0.0d00 111 endif 112 enddo 113 return 114610 Format (t4,0pf12.5,1pe13.5,T33,'x=',(T35,0p9f9.4)) 115611 Format (15x,1pe13.5,T32,'dx=',(T35,0p9f9.4)) 116612 Format(11x,'Vsp=',1pe13.5) 1176000 Format (//,' ********** The direction of the gradient vector', 118 * ' has changed by more than 90 degrees for the last 5 steps.',/, 119 * 12x,'Either the step size should be decreased, or another', 120 * ' integrator should be tried.') 121 end 122