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