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