1C $Id$
2      Subroutine drdy_path(ns)
3C
4C  Calculate steepest descent path from saddle point
5C     lgs(31)  -  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      Integer ns
13      Integer i,ierr,ii,iopt,j,nbar
14      Double Precision delss,delsh,dxmag,smax,sold,swf,xsign
15
16      Double Precision eps
17c
18#include "drdyP.fh"
19C
20      Data eps/1.d-6/
21C
22      sold = 0.0d0
23      iopt = lgs(31)
24      nbar = (n3*(n3+1))/2
25      write (fu6,600) dels,delsv,delhss,dir,slm,slp
26      if (iopt.le.1) then
27         write (fu6,601)
28      else if (iopt.eq.2) then
29         write (fu6,602)
30      else if (iopt.eq.3) then
31         write (fu6,603)
32      else
33         write (fu6,604)
34      end if
35C
36      if (lgs(1).ne.0) then
37C Set up for calculation on reactant side (if lgs(1)=0 then no saddle
38C    point and reaction path is in product direction only)
39         xsign = -1.0d00
40         smax = slm
41         if (ns.gt.0 .and. sgrid(1).lt.0.0d00) then
42C Restart from grid point with most negative s value
43            ii = 0
44            dxmag = 0.0d00
45            do i = 1,n3
46               x(i) = xgrid(i,1)
47               dx(i) = dxgrid(i,1)
48               dxmag = dxmag + dx(i)*dx(i)
49               do j = 1,i
50                  ii = ii + 1
51                  f(j,i) = hgrid(ii,1)
52                  f(i,j) = hgrid(ii,1)
53               enddo
54            enddo
55            dxmag = sqrt(dxmag)
56            if (iopt.ge.3) then
57C For PG CLQA and CUBE algorithms we need the Hessian at a previous point
58               if (ns.ge.2) then
59                  sold = sgrid(2)
60                  ii = 0
61                  do i = 1,n3
62                     do j = 1,i
63                        ii = ii + 1
64                        fold(j,i) = hgrid(ii,2)
65                        fold(i,j) = hgrid(ii,2)
66                     enddo
67                  enddo
68               else
69C Use saddle point Hessian as previous one
70                  sold = 0.0
71                  do i = 1,n3
72                     do j = 1,n3
73                        fold(j,i) = fspsv(j,i)
74                     enddo
75                  enddo
76               endif
77            endif
78            write(fu6,605)
79            write(fu6,606)
80            write(fu6,607) sgrid(1),vgrid(1),(x(i),i=1,n3)
81            write(fu6,608) dxmag,(dx(i)/dxmag,i=1,n3)
82            s = sgrid(1)
83            delss = 0.0d00
84            delsh = 0.0d00
85            do i = 1,n3
86               dxold(i) = dx(i)/dxmag
87            enddo
88C
89         else
90C Take step off saddle point in reactant direction
91            write(fu6,609)
92            write(fu6,606)
93            write(fu6,607) 0.,vspsv-vzero,(xspsv(i),i=1,n3)
94            write(fu6,608) 0.,(vec0sv(i),i=1,n3)
95            if (lgs(39).ne.0) write (fu6,6001) vspspsv-vzerosp
96            do i = 1,n3
97               x(i) = xspsv(i)-dels*(vec0sv(i)-0.5d00*dels*vec1sv(i))
98               do j = 1,n3
99                  f(j,i) = fspsv(j,i)
100               enddo
101            enddo
102            if (iopt.ge.3) then
103C For CLQA and CUBE, store hessian matrix in fold
104               sold = 0.0d00
105               do i = 1,n3
106                  do j = 1,n3
107                     fold(j,i) = f(j,i)
108                  enddo
109               enddo
110            endif
111            s = -dels
112            delss = dels
113            delsh = dels
114C Calculate g and optionally f at new geometry
115            if (iopt.lt.3 .and. delss.lt.delsv*(1.0d00-eps) .and.
116     *       delsh.lt.delhss*(1.0d00-eps)) then
117               Call drdy_pot1 (x,v,dx,scr1,amass,natom,n3)
118            else
119               Call drdy_pot2 (x,v,dx,f,hess,scr1,amass,natom,n3,
120     *            n3tm,nbar)
121               delsh = 0.0d00
122            endif
123            write(fu6,607) s,v-vzero,(x(i),i=1,n3)
124            dxmag = 0.0d00
125            do i = 1,n3
126               dxmag = dxmag + dx(i)*dx(i)
127            enddo
128            dxmag = sqrt(dxmag)
129            write(fu6,608) dxmag,(dx(i)/dxmag,i=1,n3)
130            do i = 1,n3
131               dxold(i) = dx(i)/dxmag
132            enddo
133            if (delss.ge.delsv*(1.0d00-eps)) then
134               swf = s
135               if(lgs(38).gt.0) call drdy_wfiles(swf,0)
136               if (lgs(39).ne.0) then
137                  call drdy_potsp(x,vspc,scr1,amass,natom,n3)
138                  write (fu6,6001) vspc-vzerosp
139                  if(lgs(38).eq.4) call drdy_wfiles(swf,1)
140               endif
141               call drdy_rphwrt (7)
142               if(lgs(2).ne.0) then
143                  ierr = 0
144                  call drdy_projct(ierr)
145                  if (ierr.eq.0) call drdy_fdiag(n3,nf(5),ierr,lgs(2))
146               endif
147               delss = 0.0d00
148            endif
149         endif
150C
151         Call drdy_path2(iopt,delss,delsh,sold,smax,xsign)
152C
153      endif
154C
155C Calculation on product side
156      xsign = 1.0d00
157      smax = slp
158      if (((lgs(1).ne.0.and.ns.gt.0) .or. ns.gt.1) .and.
159     *   sgrid(ns).gt.0.0d00) then
160C Restart from grid point with most positive s value
161         ii = 0
162         dxmag = 0.0
163         do i = 1,n3
164            x(i) = xgrid(i,ns)
165            dx(i) = dxgrid(i,ns)
166            dxmag = dxmag + dx(i)*dx(i)
167            do j = 1,i
168               ii = ii + 1
169               f(j,i) = hgrid(ii,ns)
170               f(i,j) = hgrid(ii,ns)
171            enddo
172         enddo
173         dxmag = sqrt(dxmag)
174         if (iopt.ge.3) then
175C For PM CLQA and CUBE algorithms we need the Hessian at a previous point
176            if (ns.ge.2) then
177               sold = sgrid(ns-1)
178               ii = 0
179               do i = 1,n3
180                  do j = 1,i
181                     ii = ii + 1
182                     fold(j,i) = hgrid(ii,ns-1)
183                     fold(i,j) = hgrid(ii,ns-1)
184                  enddo
185               enddo
186            else
187C Use saddle point Hessian as previous one (note that for lgs(1)=0,
188C  ns>1 so this section will not be called
189               sold = 0.0d00
190               do i = 1,n3
191                  do j = 1,n3
192                     fold(j,i) = fspsv(j,i)
193                  enddo
194               enddo
195            endif
196         endif
197         write(fu6,610)
198         write(fu6,606)
199         write(fu6,607) sgrid(ns),vgrid(ns),(x(i),i=1,n3)
200         write(fu6,608) dxmag,(dx(i)/dxmag,i=1,n3)
201         do i = 1,n3
202            dxold(i) = dx(i)/dxmag
203         enddo
204         s = sgrid(ns)
205         delss = 0.0d00
206         delsh = 0.0d00
207C
208      elseif (lgs(1).ne.0) then
209C Calculation for step off saddle point
210         write(fu6,611)
211         write(fu6,606)
212         write(fu6,607) 0.,vspsv-vzero,(xspsv(i),i=1,n3)
213         write(fu6,608) 0.,(vec0sv(i),i=1,n3)
214C Take step off saddle point in product direction
215         do i = 1,n3
216            x(i) = xspsv(i) + dels*(vec0sv(i) +
217     &           0.5d0*dels*vec1sv(i))
218            do j = 1,n3
219               f(j,i) = fspsv(j,i)
220            enddo
221         enddo
222         if (iopt.ge.3) then
223C For CLQA and CUBE, store hessian matrix in fold
224            sold = 0.0d00
225            do i = 1,n3
226               do j = 1,n3
227                  fold(j,i) = f(j,i)
228               enddo
229            enddo
230         endif
231         s = dels
232         delss = dels
233         delsh = dels
234C Calculate g and optionally f at new geometry
235         if (iopt.lt.3 .and. delss.lt.delsv*(1.0d00-eps) .and.
236     *    delsh.lt.delhss*(1.0-eps)) then
237            Call drdy_pot1 (x,v,dx,scr1,amass,natom,n3)
238         else
239            Call drdy_pot2 (x,v,dx,f,hess,scr1,amass,natom,n3,
240     *         n3tm,nbar)
241            delsh = 0.0d00
242         endif
243         write(fu6,607) s,v-vzero,(x(i),i=1,n3)
244         dxmag = 0.0
245         do i = 1,n3
246            dxmag = dxmag + dx(i)*dx(i)
247         enddo
248         dxmag = sqrt(dxmag)
249         write(fu6,608) dxmag,(dx(i)/dxmag,i=1,n3)
250         do i = 1,n3
251            dxold(i) = dx(i)/dxmag
252         enddo
253         if (delss.ge.delsv*(1.0d00-eps)) then
254            swf = s
255            if(lgs(38).gt.0) call drdy_wfiles(swf,0)
256            if (lgs(39).ne.0) then
257               call drdy_potsp(x,vspc,scr1,amass,natom,n3)
258               write (fu6,6001) vspc-vzerosp
259               if(lgs(38).eq.4) call drdy_wfiles(swf,1)
260            endif
261            call drdy_rphwrt (7)
262            if(lgs(2).ne.0) then
263               ierr = 0
264               call drdy_projct(ierr)
265               if (ierr.eq.0) call drdy_fdiag(n3,nf(5),ierr,lgs(2))
266            endif
267            delss = 0.0d00
268         endif
269      else
270C No saddle point, take step from starting geometry
271         dxmag = 0.0d00
272         do i = 1,n3
273            x(i) = xspsv(i)
274            dx(i) = dxspsv(i)
275            dxmag = dxmag + dx(i)*dx(i)
276            do j = 1,n3
277               f(j,i) = fspsv(j,i)
278C For CLQA and CUBE, store hessian matrix in fold
279               fold(j,i) = fspsv(j,i)
280            enddo
281         enddo
282         dxmag = sqrt(dxmag)
283         write(fu6,612)
284         write(fu6,606)
285         write(fu6,607) 0.,vspsv-vzero,(x(i),i=1,n3)
286         write(fu6,608) dxmag,(dx(i)/dxmag,i=1,n3)
287         if (iopt.le.1) then
288            call drdy_euler (n3tm,n3,x,dx,dels)
289         else
290            ierr = 0
291            call drdy_pmlqa (n3tm,n3,nbar,x,dx,f,hess,amass,
292     *       dels,vec0,vec1,u0,scr1,vec2,ierr)
293            sold = 0.0
294         endif
295         s = dels
296         delss = dels
297         delsh = dels
298C Calculate g and f at new geometry
299         if (iopt.lt.3 .and. delss.lt.delsv*(1.0d00-eps) .and.
300     *    delsh.lt.delhss*(1.0d00-eps)) then
301            Call drdy_pot1 (x,v,dx,scr1,amass,natom,n3)
302         else
303            Call drdy_pot2 (x,v,dx,f,hess,scr1,amass,natom,n3,
304     *        n3tm,nbar)
305            delsh = 0.0d00
306         endif
307         write(fu6,607) s,v-vzero,(x(i),i=1,n3)
308         dxmag = 0.0d00
309         do i = 1,n3
310            dxmag = dxmag + dx(i)*dx(i)
311         enddo
312         dxmag = sqrt(dxmag)
313         write(fu6,608) dxmag,(dx(i)/dxmag,i=1,n3)
314         do i = 1,n3
315            dxold(i) = dx(i)/dxmag
316         enddo
317         if (delss.ge.delsv*(1.0d00-eps)) then
318            swf = s
319            if(lgs(38).gt.0) call drdy_wfiles(swf,0)
320            if (lgs(39).ne.0) then
321               call drdy_potsp(x,vspc,scr1,amass,natom,n3)
322               write (fu6,6001) vspc-vzerosp
323               if(lgs(38).eq.4) call drdy_wfiles(swf,1)
324            endif
325            call drdy_rphwrt (7)
326            if(lgs(2).ne.0) then
327               ierr = 0
328               call drdy_projct(ierr)
329               if (ierr.eq.0) call drdy_fdiag(n3,nf(5),ierr,lgs(2))
330            endif
331            delss = 0.0
332         endif
333      endif
334C
335      Call drdy_path2(iopt,delss,delsh,sold,smax,xsign)
336C
337      Return
338600   Format(//,' ********** Calculate points along reaction path',//,
339     *   5x,' dels=',T20,1pe13.5,/,5x,' delsv=',T20,1pe13.5,/,
340     *   5x,' delhss=',T20,1pe13.5,/,5x,' dir=',T20,1pe13.5,/,
341     *   5x,' slm=', T20,1pe13.5,/,5x,' slp=',T20,1pe13.5)
342601   Format(/,' Path followed using Euler integrator')
343602   Format(/,' Path followed using Page-McIver LQA algorithm')
344603   Format(/,' Path followed using Page-McIver CLQA algorithm')
345604   Format(/,' Path followed using Page-McIver CUBE algorithm')
346605   Format(/,' Take step from restart geometry in reactant',
347     *   ' direction')
348606   Format(/,'  All quantities in atomic units, x is vector of ',
349     *   'mass-weighted coordinates,',/,'  dx is the normalized mass-',
350     *   'weighted gradient vector, and |dV/x| is the magnitude of',
351     *   ' the gradient',/,9x,'s',7x,'V,|dV/dx|')
352607   Format (t4,0pf12.5,1pe13.5,T33,'x=',(T35,0p9f9.4))
353608   Format (15x,1pe13.5,T32,'dx=',(T35,0p9f9.4))
3546001  Format(11x,'Vsp=',1pe13.5)
355609   Format(/,' Take step from saddle point in reactant direction')
356610   Format(/,' Take step from restart geometry in product',
357     *   ' direction')
358611   Format(/,' Take step from saddle point in product direction')
359612   Format(/,' Take step from initial geometry along gradient')
360      End
361