1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 2003 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ################################################################# 9c ## ## 10c ## subroutine epitors2 -- pi-system torsion Hessian; numer ## 11c ## ## 12c ################################################################# 13c 14c 15c "epitors2" calculates the second derivatives of the pi-system 16c torsion energy for a single atom using finite difference methods 17c 18c 19 subroutine epitors2 (i) 20 use angbnd 21 use atoms 22 use group 23 use hessn 24 use pitors 25 use usage 26 implicit none 27 integer i,j,ipitors 28 integer ia,ib,ic 29 integer id,ie,ig 30 real*8 eps,fgrp 31 real*8 old,term 32 real*8, allocatable :: de(:,:) 33 real*8, allocatable :: d0(:,:) 34 logical proceed 35 logical twosided 36c 37c 38c set stepsize for derivatives and default group weight 39c 40 eps = 1.0d-5 41 fgrp = 1.0d0 42 twosided = .false. 43 if (n .le. 50) twosided = .true. 44c 45c perform dynamic allocation of some local arrays 46c 47 allocate (de(3,n)) 48 allocate (d0(3,n)) 49c 50c calculate numerical pi-system torsion Hessian for current atom 51c 52 do ipitors = 1, npitors 53 ia = ipit(1,ipitors) 54 ib = ipit(2,ipitors) 55 ic = ipit(3,ipitors) 56 id = ipit(4,ipitors) 57 ie = ipit(5,ipitors) 58 ig = ipit(6,ipitors) 59c 60c decide whether to compute the current interaction 61c 62 proceed = .true. 63 if (use_group) call groups (proceed,fgrp,ia,ib,ic,id,ie,ig) 64 if (proceed) proceed = (use(ia) .or. use(ib) .or. use(ic) .or. 65 & use(id) .or. use(ie) .or. use(ig)) 66 if (proceed) then 67 term = fgrp / eps 68c 69c find first derivatives for the base structure 70c 71 if (.not. twosided) then 72 call epitors2a (ipitors,de) 73 do j = 1, 3 74 d0(j,ia) = de(j,ia) 75 d0(j,ib) = de(j,ib) 76 d0(j,ic) = de(j,ic) 77 d0(j,id) = de(j,id) 78 d0(j,ie) = de(j,ie) 79 d0(j,ig) = de(j,ig) 80 end do 81 end if 82c 83c find numerical x-components via perturbed structures 84c 85 old = x(i) 86 if (twosided) then 87 x(i) = x(i) - 0.5d0*eps 88 call epitors2a (ipitors,de) 89 do j = 1, 3 90 d0(j,ia) = de(j,ia) 91 d0(j,ib) = de(j,ib) 92 d0(j,ic) = de(j,ic) 93 d0(j,id) = de(j,id) 94 d0(j,ie) = de(j,ie) 95 d0(j,ig) = de(j,ig) 96 end do 97 end if 98 x(i) = x(i) + eps 99 call epitors2a (ipitors,de) 100 x(i) = old 101 do j = 1, 3 102 hessx(j,ia) = hessx(j,ia) + term*(de(j,ia)-d0(j,ia)) 103 hessx(j,ib) = hessx(j,ib) + term*(de(j,ib)-d0(j,ib)) 104 hessx(j,ic) = hessx(j,ic) + term*(de(j,ic)-d0(j,ic)) 105 hessx(j,id) = hessx(j,id) + term*(de(j,id)-d0(j,id)) 106 hessx(j,ie) = hessx(j,ie) + term*(de(j,ie)-d0(j,ie)) 107 hessx(j,ig) = hessx(j,ig) + term*(de(j,ig)-d0(j,ig)) 108 end do 109c 110c find numerical y-components via perturbed structures 111c 112 old = y(i) 113 if (twosided) then 114 y(i) = y(i) - 0.5d0*eps 115 call epitors2a (ipitors,de) 116 do j = 1, 3 117 d0(j,ia) = de(j,ia) 118 d0(j,ib) = de(j,ib) 119 d0(j,ic) = de(j,ic) 120 d0(j,id) = de(j,id) 121 d0(j,ie) = de(j,ie) 122 d0(j,ig) = de(j,ig) 123 end do 124 end if 125 y(i) = y(i) + eps 126 call epitors2a (ipitors,de) 127 y(i) = old 128 do j = 1, 3 129 hessy(j,ia) = hessy(j,ia) + term*(de(j,ia)-d0(j,ia)) 130 hessy(j,ib) = hessy(j,ib) + term*(de(j,ib)-d0(j,ib)) 131 hessy(j,ic) = hessy(j,ic) + term*(de(j,ic)-d0(j,ic)) 132 hessy(j,id) = hessy(j,id) + term*(de(j,id)-d0(j,id)) 133 hessy(j,ie) = hessy(j,ie) + term*(de(j,ie)-d0(j,ie)) 134 hessy(j,ig) = hessy(j,ig) + term*(de(j,ig)-d0(j,ig)) 135 end do 136c 137c find numerical z-components via perturbed structures 138c 139 old = z(i) 140 if (twosided) then 141 z(i) = z(i) - 0.5d0*eps 142 call epitors2a (ipitors,de) 143 do j = 1, 3 144 d0(j,ia) = de(j,ia) 145 d0(j,ib) = de(j,ib) 146 d0(j,ic) = de(j,ic) 147 d0(j,id) = de(j,id) 148 d0(j,ie) = de(j,ie) 149 d0(j,ig) = de(j,ig) 150 end do 151 end if 152 z(i) = z(i) + eps 153 call epitors2a (ipitors,de) 154 z(i) = old 155 do j = 1, 3 156 hessz(j,ia) = hessz(j,ia) + term*(de(j,ia)-d0(j,ia)) 157 hessz(j,ib) = hessz(j,ib) + term*(de(j,ib)-d0(j,ib)) 158 hessz(j,ic) = hessz(j,ic) + term*(de(j,ic)-d0(j,ic)) 159 hessz(j,id) = hessz(j,id) + term*(de(j,id)-d0(j,id)) 160 hessz(j,ie) = hessz(j,ie) + term*(de(j,ie)-d0(j,ie)) 161 hessz(j,ig) = hessz(j,ig) + term*(de(j,ig)-d0(j,ig)) 162 end do 163 end if 164 end do 165c 166c perform deallocation of some local arrays 167c 168 deallocate (d0) 169 return 170 end 171c 172c 173c ############################################################### 174c ## ## 175c ## subroutine epitors2a -- pi-system torsion derivatives ## 176c ## ## 177c ############################################################### 178c 179c 180c "epitors2a" calculates the pi-system torsion first derivatives; 181c used in computation of finite difference second derivatives 182c 183c 184 subroutine epitors2a (i,de) 185 use atoms 186 use bound 187 use deriv 188 use pitors 189 use torpot 190 implicit none 191 integer i,ia,ib,ic 192 integer id,ie,ig 193 real*8 dedphi,dphi2 194 real*8 xt,yt,zt,rt2 195 real*8 xu,yu,zu,ru2 196 real*8 xtu,ytu,ztu 197 real*8 rdc,rtru 198 real*8 v2,c2,s2 199 real*8 sine,cosine 200 real*8 sine2,cosine2 201 real*8 xia,yia,zia 202 real*8 xib,yib,zib 203 real*8 xic,yic,zic 204 real*8 xid,yid,zid 205 real*8 xie,yie,zie 206 real*8 xig,yig,zig 207 real*8 xip,yip,zip 208 real*8 xiq,yiq,ziq 209 real*8 xad,yad,zad 210 real*8 xbd,ybd,zbd 211 real*8 xec,yec,zec 212 real*8 xgc,ygc,zgc 213 real*8 xcp,ycp,zcp 214 real*8 xdc,ydc,zdc 215 real*8 xqd,yqd,zqd 216 real*8 xdp,ydp,zdp 217 real*8 xqc,yqc,zqc 218 real*8 dedxt,dedyt,dedzt 219 real*8 dedxu,dedyu,dedzu 220 real*8 dedxia,dedyia,dedzia 221 real*8 dedxib,dedyib,dedzib 222 real*8 dedxic,dedyic,dedzic 223 real*8 dedxid,dedyid,dedzid 224 real*8 dedxie,dedyie,dedzie 225 real*8 dedxig,dedyig,dedzig 226 real*8 dedxip,dedyip,dedzip 227 real*8 dedxiq,dedyiq,dedziq 228 real*8 de(3,*) 229c 230c 231c set the atom numbers for this pi-system torsion 232c 233 ia = ipit(1,i) 234 ib = ipit(2,i) 235 ic = ipit(3,i) 236 id = ipit(4,i) 237 ie = ipit(5,i) 238 ig = ipit(6,i) 239c 240c compute the value of the pi-system torsion angle 241c 242 xia = x(ia) 243 yia = y(ia) 244 zia = z(ia) 245 xib = x(ib) 246 yib = y(ib) 247 zib = z(ib) 248 xic = x(ic) 249 yic = y(ic) 250 zic = z(ic) 251 xid = x(id) 252 yid = y(id) 253 zid = z(id) 254 xie = x(ie) 255 yie = y(ie) 256 zie = z(ie) 257 xig = x(ig) 258 yig = y(ig) 259 zig = z(ig) 260 xad = xia - xid 261 yad = yia - yid 262 zad = zia - zid 263 xbd = xib - xid 264 ybd = yib - yid 265 zbd = zib - zid 266 xec = xie - xic 267 yec = yie - yic 268 zec = zie - zic 269 xgc = xig - xic 270 ygc = yig - yic 271 zgc = zig - zic 272 if (use_polymer) then 273 call image (xad,yad,zad) 274 call image (xbd,ybd,zbd) 275 call image (xec,yec,zec) 276 call image (xgc,ygc,zgc) 277 end if 278 xip = yad*zbd - ybd*zad + xic 279 yip = zad*xbd - zbd*xad + yic 280 zip = xad*ybd - xbd*yad + zic 281 xiq = yec*zgc - ygc*zec + xid 282 yiq = zec*xgc - zgc*xec + yid 283 ziq = xec*ygc - xgc*yec + zid 284 xcp = xic - xip 285 ycp = yic - yip 286 zcp = zic - zip 287 xdc = xid - xic 288 ydc = yid - yic 289 zdc = zid - zic 290 xqd = xiq - xid 291 yqd = yiq - yid 292 zqd = ziq - zid 293 if (use_polymer) then 294 call image (xcp,ycp,zcp) 295 call image (xdc,ydc,zdc) 296 call image (xqd,yqd,zqd) 297 end if 298 xt = ycp*zdc - ydc*zcp 299 yt = zcp*xdc - zdc*xcp 300 zt = xcp*ydc - xdc*ycp 301 xu = ydc*zqd - yqd*zdc 302 yu = zdc*xqd - zqd*xdc 303 zu = xdc*yqd - xqd*ydc 304 xtu = yt*zu - yu*zt 305 ytu = zt*xu - zu*xt 306 ztu = xt*yu - xu*yt 307 rt2 = xt*xt + yt*yt + zt*zt 308 ru2 = xu*xu + yu*yu + zu*zu 309 rtru = sqrt(rt2 * ru2) 310 if (rtru .ne. 0.0d0) then 311 rdc = sqrt(xdc*xdc + ydc*ydc + zdc*zdc) 312 cosine = (xt*xu + yt*yu + zt*zu) / rtru 313 sine = (xdc*xtu + ydc*ytu + zdc*ztu) / (rdc*rtru) 314c 315c set the pi-system torsion parameters for this angle 316c 317 v2 = kpit(i) 318 c2 = -1.0d0 319 s2 = 0.0d0 320c 321c compute the multiple angle trigonometry and the phase terms 322c 323 cosine2 = cosine*cosine - sine*sine 324 sine2 = 2.0d0 * cosine * sine 325 dphi2 = 2.0d0 * (cosine2*s2 - sine2*c2) 326c 327c calculate pi-system torsion energy and master chain rule term 328c 329 dedphi = ptorunit * v2 * dphi2 330c 331c chain rule terms for first derivative components 332c 333 xdp = xid - xip 334 ydp = yid - yip 335 zdp = zid - zip 336 xqc = xiq - xic 337 yqc = yiq - yic 338 zqc = ziq - zic 339 dedxt = dedphi * (yt*zdc - ydc*zt) / (rt2*rdc) 340 dedyt = dedphi * (zt*xdc - zdc*xt) / (rt2*rdc) 341 dedzt = dedphi * (xt*ydc - xdc*yt) / (rt2*rdc) 342 dedxu = -dedphi * (yu*zdc - ydc*zu) / (ru2*rdc) 343 dedyu = -dedphi * (zu*xdc - zdc*xu) / (ru2*rdc) 344 dedzu = -dedphi * (xu*ydc - xdc*yu) / (ru2*rdc) 345c 346c compute first derivative components for pi-system angle 347c 348 dedxip = zdc*dedyt - ydc*dedzt 349 dedyip = xdc*dedzt - zdc*dedxt 350 dedzip = ydc*dedxt - xdc*dedyt 351 dedxic = ydp*dedzt - zdp*dedyt + zqd*dedyu - yqd*dedzu 352 dedyic = zdp*dedxt - xdp*dedzt + xqd*dedzu - zqd*dedxu 353 dedzic = xdp*dedyt - ydp*dedxt + yqd*dedxu - xqd*dedyu 354 dedxid = zcp*dedyt - ycp*dedzt + yqc*dedzu - zqc*dedyu 355 dedyid = xcp*dedzt - zcp*dedxt + zqc*dedxu - xqc*dedzu 356 dedzid = ycp*dedxt - xcp*dedyt + xqc*dedyu - yqc*dedxu 357 dedxiq = zdc*dedyu - ydc*dedzu 358 dedyiq = xdc*dedzu - zdc*dedxu 359 dedziq = ydc*dedxu - xdc*dedyu 360c 361c compute first derivative components for individual atoms 362c 363 dedxia = ybd*dedzip - zbd*dedyip 364 dedyia = zbd*dedxip - xbd*dedzip 365 dedzia = xbd*dedyip - ybd*dedxip 366 dedxib = zad*dedyip - yad*dedzip 367 dedyib = xad*dedzip - zad*dedxip 368 dedzib = yad*dedxip - xad*dedyip 369 dedxie = ygc*dedziq - zgc*dedyiq 370 dedyie = zgc*dedxiq - xgc*dedziq 371 dedzie = xgc*dedyiq - ygc*dedxiq 372 dedxig = zec*dedyiq - yec*dedziq 373 dedyig = xec*dedziq - zec*dedxiq 374 dedzig = yec*dedxiq - xec*dedyiq 375 dedxic = dedxic + dedxip - dedxie - dedxig 376 dedyic = dedyic + dedyip - dedyie - dedyig 377 dedzic = dedzic + dedzip - dedzie - dedzig 378 dedxid = dedxid + dedxiq - dedxia - dedxib 379 dedyid = dedyid + dedyiq - dedyia - dedyib 380 dedzid = dedzid + dedziq - dedzia - dedzib 381c 382c increment the total pi-system torsion energy and gradient 383c 384 de(1,ia) = dedxia 385 de(2,ia) = dedyia 386 de(3,ia) = dedzia 387 de(1,ib) = dedxib 388 de(2,ib) = dedyib 389 de(3,ib) = dedzib 390 de(1,ic) = dedxic 391 de(2,ic) = dedyic 392 de(3,ic) = dedzic 393 de(1,id) = dedxid 394 de(2,id) = dedyid 395 de(3,id) = dedzid 396 de(1,ie) = dedxie 397 de(2,ie) = dedyie 398 de(3,ie) = dedzie 399 de(1,ig) = dedxig 400 de(2,ig) = dedyig 401 de(3,ig) = dedzig 402 end if 403 return 404 end 405