1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################# 9c ## ## 10c ## subroutine etors2 -- atom-by-atom torsional Hessian ## 11c ## ## 12c ############################################################# 13c 14c 15c "etors2" calculates the second derivatives of the torsional 16c energy for a single atom 17c 18c 19 subroutine etors2 (i) 20 use warp 21 implicit none 22 integer i 23c 24c 25c choose standard or potential energy smoothing version 26c 27 if (use_smooth) then 28 call etors2b (i) 29 else 30 call etors2a (i) 31 end if 32 return 33 end 34c 35c 36c ################################################################ 37c ## ## 38c ## subroutine etors2a -- standard torsional angle Hessian ## 39c ## ## 40c ################################################################ 41c 42c 43c "etors2a" calculates the second derivatives of the torsional 44c energy for a single atom using a standard sum of Fourier terms 45c 46c 47 subroutine etors2a (i) 48 use atoms 49 use bound 50 use group 51 use hessn 52 use torpot 53 use tors 54 implicit none 55 integer i,ktors 56 integer ia,ib,ic,id 57 real*8 eps,fgrp 58 real*8 dedphi,d2edphi2 59 real*8 v1,v2,v3,v4,v5,v6 60 real*8 c1,c2,c3,c4,c5,c6 61 real*8 s1,s2,s3,s4,s5,s6 62 real*8 cosine,cosine2 63 real*8 cosine3,cosine4 64 real*8 cosine5,cosine6 65 real*8 sine,sine2,sine3 66 real*8 sine4,sine5,sine6 67 real*8 xia,yia,zia 68 real*8 xib,yib,zib 69 real*8 xic,yic,zic 70 real*8 xid,yid,zid 71 real*8 xba,yba,zba 72 real*8 xcb,ycb,zcb 73 real*8 xdc,ydc,zdc 74 real*8 xca,yca,zca 75 real*8 xdb,ydb,zdb 76 real*8 xt,yt,zt,xu,yu,zu 77 real*8 xtu,ytu,ztu 78 real*8 rt2,ru2,rtru,rcb 79 real*8 dphi1,dphi2,dphi3 80 real*8 dphi4,dphi5,dphi6 81 real*8 d2phi1,d2phi2,d2phi3 82 real*8 d2phi4,d2phi5,d2phi6 83 real*8 dphidxt,dphidyt,dphidzt 84 real*8 dphidxu,dphidyu,dphidzu 85 real*8 dphidxia,dphidyia,dphidzia 86 real*8 dphidxib,dphidyib,dphidzib 87 real*8 dphidxic,dphidyic,dphidzic 88 real*8 dphidxid,dphidyid,dphidzid 89 real*8 xycb2,xzcb2,yzcb2 90 real*8 rcbxt,rcbyt,rcbzt,rcbt2 91 real*8 rcbxu,rcbyu,rcbzu,rcbu2 92 real*8 dphidxibt,dphidyibt,dphidzibt 93 real*8 dphidxibu,dphidyibu,dphidzibu 94 real*8 dphidxict,dphidyict,dphidzict 95 real*8 dphidxicu,dphidyicu,dphidzicu 96 real*8 dxiaxia,dyiayia,dziazia 97 real*8 dxibxib,dyibyib,dzibzib 98 real*8 dxicxic,dyicyic,dziczic 99 real*8 dxidxid,dyidyid,dzidzid 100 real*8 dxiayia,dxiazia,dyiazia 101 real*8 dxibyib,dxibzib,dyibzib 102 real*8 dxicyic,dxiczic,dyiczic 103 real*8 dxidyid,dxidzid,dyidzid 104 real*8 dxiaxib,dxiayib,dxiazib 105 real*8 dyiaxib,dyiayib,dyiazib 106 real*8 dziaxib,dziayib,dziazib 107 real*8 dxiaxic,dxiayic,dxiazic 108 real*8 dyiaxic,dyiayic,dyiazic 109 real*8 dziaxic,dziayic,dziazic 110 real*8 dxiaxid,dxiayid,dxiazid 111 real*8 dyiaxid,dyiayid,dyiazid 112 real*8 dziaxid,dziayid,dziazid 113 real*8 dxibxic,dxibyic,dxibzic 114 real*8 dyibxic,dyibyic,dyibzic 115 real*8 dzibxic,dzibyic,dzibzic 116 real*8 dxibxid,dxibyid,dxibzid 117 real*8 dyibxid,dyibyid,dyibzid 118 real*8 dzibxid,dzibyid,dzibzid 119 real*8 dxicxid,dxicyid,dxiczid 120 real*8 dyicxid,dyicyid,dyiczid 121 real*8 dzicxid,dzicyid,dziczid 122 logical proceed 123c 124c 125c set tolerance for minimum distance and angle values 126c 127 eps = 0.0001d0 128c 129c calculate the torsional angle interaction Hessian elements 130c 131 do ktors = 1, ntors 132 ia = itors(1,ktors) 133 ib = itors(2,ktors) 134 ic = itors(3,ktors) 135 id = itors(4,ktors) 136c 137c decide whether to compute the current interaction 138c 139 proceed = (i.eq.ia .or. i.eq.ib .or. i.eq.ic .or. i.eq.id) 140 if (proceed .and. use_group) 141 & call groups (proceed,fgrp,ia,ib,ic,id,0,0) 142c 143c compute the value of the torsional angle 144c 145 if (proceed) then 146 xia = x(ia) 147 yia = y(ia) 148 zia = z(ia) 149 xib = x(ib) 150 yib = y(ib) 151 zib = z(ib) 152 xic = x(ic) 153 yic = y(ic) 154 zic = z(ic) 155 xid = x(id) 156 yid = y(id) 157 zid = z(id) 158 xba = xib - xia 159 yba = yib - yia 160 zba = zib - zia 161 xcb = xic - xib 162 ycb = yic - yib 163 zcb = zic - zib 164 xdc = xid - xic 165 ydc = yid - yic 166 zdc = zid - zic 167 if (use_polymer) then 168 call image (xba,yba,zba) 169 call image (xcb,ycb,zcb) 170 call image (xdc,ydc,zdc) 171 end if 172 rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps)) 173 xt = yba*zcb - ycb*zba 174 yt = zba*xcb - zcb*xba 175 zt = xba*ycb - xcb*yba 176 xu = ycb*zdc - ydc*zcb 177 yu = zcb*xdc - zdc*xcb 178 zu = xcb*ydc - xdc*ycb 179 xtu = yt*zu - yu*zt 180 ytu = zt*xu - zu*xt 181 ztu = xt*yu - xu*yt 182 rt2 = max(xt*xt+yt*yt+zt*zt,eps) 183 ru2 = max(xu*xu+yu*yu+zu*zu,eps) 184 rtru = sqrt(rt2 * ru2) 185 cosine = (xt*xu + yt*yu + zt*zu) / rtru 186 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru) 187c 188c set the torsional parameters for this angle 189c 190 v1 = tors1(1,ktors) 191 c1 = tors1(3,ktors) 192 s1 = tors1(4,ktors) 193 v2 = tors2(1,ktors) 194 c2 = tors2(3,ktors) 195 s2 = tors2(4,ktors) 196 v3 = tors3(1,ktors) 197 c3 = tors3(3,ktors) 198 s3 = tors3(4,ktors) 199 v4 = tors4(1,ktors) 200 c4 = tors4(3,ktors) 201 s4 = tors4(4,ktors) 202 v5 = tors5(1,ktors) 203 c5 = tors5(3,ktors) 204 s5 = tors5(4,ktors) 205 v6 = tors6(1,ktors) 206 c6 = tors6(3,ktors) 207 s6 = tors6(4,ktors) 208c 209c compute the multiple angle trigonometry and the phase terms 210c 211 cosine2 = cosine*cosine - sine*sine 212 sine2 = 2.0d0 * cosine * sine 213 cosine3 = cosine*cosine2 - sine*sine2 214 sine3 = cosine*sine2 + sine*cosine2 215 cosine4 = cosine*cosine3 - sine*sine3 216 sine4 = cosine*sine3 + sine*cosine3 217 cosine5 = cosine*cosine4 - sine*sine4 218 sine5 = cosine*sine4 + sine*cosine4 219 cosine6 = cosine*cosine5 - sine*sine5 220 sine6 = cosine*sine5 + sine*cosine5 221 dphi1 = (cosine*s1 - sine*c1) 222 dphi2 = 2.0d0 * (cosine2*s2 - sine2*c2) 223 dphi3 = 3.0d0 * (cosine3*s3 - sine3*c3) 224 dphi4 = 4.0d0 * (cosine4*s4 - sine4*c4) 225 dphi5 = 5.0d0 * (cosine5*s5 - sine5*c5) 226 dphi6 = 6.0d0 * (cosine6*s6 - sine6*c6) 227 d2phi1 = -(cosine*c1 + sine*s1) 228 d2phi2 = -4.0d0 * (cosine2*c2 + sine2*s2) 229 d2phi3 = -9.0d0 * (cosine3*c3 + sine3*s3) 230 d2phi4 = -16.0d0 * (cosine4*c4 + sine4*s4) 231 d2phi5 = -25.0d0 * (cosine5*c5 + sine5*s5) 232 d2phi6 = -36.0d0 * (cosine6*c6 + sine6*s6) 233c 234c calculate the torsional master chain rule terms 235c 236 dedphi = torsunit * (v1*dphi1 + v2*dphi2 + v3*dphi3 237 & + v4*dphi4 + v5*dphi5 + v6*dphi6) 238 d2edphi2 = torsunit * (v1*d2phi1 + v2*d2phi2 + v3*d2phi3 239 & + v4*d2phi4 + v5*d2phi5 + v6*d2phi6) 240c 241c scale the interaction based on its group membership 242c 243 if (use_group) then 244 dedphi = dedphi * fgrp 245 d2edphi2 = d2edphi2 * fgrp 246 end if 247c 248c abbreviations for first derivative chain rule terms 249c 250 xca = xic - xia 251 yca = yic - yia 252 zca = zic - zia 253 xdb = xid - xib 254 ydb = yid - yib 255 zdb = zid - zib 256 if (use_polymer) then 257 call image (xca,yca,zca) 258 call image (xdb,ydb,zdb) 259 end if 260 dphidxt = (yt*zcb - ycb*zt) / (rt2*rcb) 261 dphidyt = (zt*xcb - zcb*xt) / (rt2*rcb) 262 dphidzt = (xt*ycb - xcb*yt) / (rt2*rcb) 263 dphidxu = -(yu*zcb - ycb*zu) / (ru2*rcb) 264 dphidyu = -(zu*xcb - zcb*xu) / (ru2*rcb) 265 dphidzu = -(xu*ycb - xcb*yu) / (ru2*rcb) 266c 267c abbreviations for second derivative chain rule terms 268c 269 xycb2 = xcb*xcb + ycb*ycb 270 xzcb2 = xcb*xcb + zcb*zcb 271 yzcb2 = ycb*ycb + zcb*zcb 272 rcbxt = -2.0d0 * rcb * dphidxt 273 rcbyt = -2.0d0 * rcb * dphidyt 274 rcbzt = -2.0d0 * rcb * dphidzt 275 rcbt2 = rcb * rt2 276 rcbxu = 2.0d0 * rcb * dphidxu 277 rcbyu = 2.0d0 * rcb * dphidyu 278 rcbzu = 2.0d0 * rcb * dphidzu 279 rcbu2 = rcb * ru2 280 dphidxibt = yca*dphidzt - zca*dphidyt 281 dphidxibu = zdc*dphidyu - ydc*dphidzu 282 dphidyibt = zca*dphidxt - xca*dphidzt 283 dphidyibu = xdc*dphidzu - zdc*dphidxu 284 dphidzibt = xca*dphidyt - yca*dphidxt 285 dphidzibu = ydc*dphidxu - xdc*dphidyu 286 dphidxict = zba*dphidyt - yba*dphidzt 287 dphidxicu = ydb*dphidzu - zdb*dphidyu 288 dphidyict = xba*dphidzt - zba*dphidxt 289 dphidyicu = zdb*dphidxu - xdb*dphidzu 290 dphidzict = yba*dphidxt - xba*dphidyt 291 dphidzicu = xdb*dphidyu - ydb*dphidxu 292c 293c chain rule terms for first derivative components 294c 295 dphidxia = zcb*dphidyt - ycb*dphidzt 296 dphidyia = xcb*dphidzt - zcb*dphidxt 297 dphidzia = ycb*dphidxt - xcb*dphidyt 298 dphidxib = dphidxibt + dphidxibu 299 dphidyib = dphidyibt + dphidyibu 300 dphidzib = dphidzibt + dphidzibu 301 dphidxic = dphidxict + dphidxicu 302 dphidyic = dphidyict + dphidyicu 303 dphidzic = dphidzict + dphidzicu 304 dphidxid = zcb*dphidyu - ycb*dphidzu 305 dphidyid = xcb*dphidzu - zcb*dphidxu 306 dphidzid = ycb*dphidxu - xcb*dphidyu 307c 308c chain rule terms for second derivative components 309c 310 dxiaxia = rcbxt*dphidxia 311 dxiayia = rcbxt*dphidyia - zcb*rcb/rt2 312 dxiazia = rcbxt*dphidzia + ycb*rcb/rt2 313 dxiaxic = rcbxt*dphidxict + xcb*xt/rcbt2 314 dxiayic = rcbxt*dphidyict - dphidzt 315 & - (xba*zcb*xcb+zba*yzcb2)/rcbt2 316 dxiazic = rcbxt*dphidzict + dphidyt 317 & + (xba*ycb*xcb+yba*yzcb2)/rcbt2 318 dxiaxid = 0.0d0 319 dxiayid = 0.0d0 320 dxiazid = 0.0d0 321 dyiayia = rcbyt*dphidyia 322 dyiazia = rcbyt*dphidzia - xcb*rcb/rt2 323 dyiaxib = rcbyt*dphidxibt - dphidzt 324 & - (yca*zcb*ycb+zca*xzcb2)/rcbt2 325 dyiaxic = rcbyt*dphidxict + dphidzt 326 & + (yba*zcb*ycb+zba*xzcb2)/rcbt2 327 dyiayic = rcbyt*dphidyict + ycb*yt/rcbt2 328 dyiazic = rcbyt*dphidzict - dphidxt 329 & - (yba*xcb*ycb+xba*xzcb2)/rcbt2 330 dyiaxid = 0.0d0 331 dyiayid = 0.0d0 332 dyiazid = 0.0d0 333 dziazia = rcbzt*dphidzia 334 dziaxib = rcbzt*dphidxibt + dphidyt 335 & + (zca*ycb*zcb+yca*xycb2)/rcbt2 336 dziayib = rcbzt*dphidyibt - dphidxt 337 & - (zca*xcb*zcb+xca*xycb2)/rcbt2 338 dziaxic = rcbzt*dphidxict - dphidyt 339 & - (zba*ycb*zcb+yba*xycb2)/rcbt2 340 dziayic = rcbzt*dphidyict + dphidxt 341 & + (zba*xcb*zcb+xba*xycb2)/rcbt2 342 dziazic = rcbzt*dphidzict + zcb*zt/rcbt2 343 dziaxid = 0.0d0 344 dziayid = 0.0d0 345 dziazid = 0.0d0 346 dxibxic = -xcb*dphidxib/(rcb*rcb) 347 & - (yca*(zba*xcb+yt)-zca*(yba*xcb-zt))/rcbt2 348 & - 2.0d0*(yt*zba-yba*zt)*dphidxibt/rt2 349 & - (zdc*(ydb*xcb+zu)-ydc*(zdb*xcb-yu))/rcbu2 350 & + 2.0d0*(yu*zdb-ydb*zu)*dphidxibu/ru2 351 dxibyic = -ycb*dphidxib/(rcb*rcb) + dphidzt + dphidzu 352 & - (yca*(zba*ycb-xt)+zca*(xba*xcb+zcb*zba))/rcbt2 353 & - 2.0d0*(zt*xba-zba*xt)*dphidxibt/rt2 354 & + (zdc*(xdb*xcb+zcb*zdb)+ydc*(zdb*ycb+xu))/rcbu2 355 & + 2.0d0*(zu*xdb-zdb*xu)*dphidxibu/ru2 356 dxibxid = rcbxu*dphidxibu + xcb*xu/rcbu2 357 dxibyid = rcbyu*dphidxibu - dphidzu 358 & - (ydc*zcb*ycb+zdc*xzcb2)/rcbu2 359 dxibzid = rcbzu*dphidxibu + dphidyu 360 & + (zdc*ycb*zcb+ydc*xycb2)/rcbu2 361 dyibzib = ycb*dphidzib/(rcb*rcb) 362 & - (xca*(xca*xcb+zcb*zca)+yca*(ycb*xca+zt))/rcbt2 363 & - 2.0d0*(xt*zca-xca*zt)*dphidzibt/rt2 364 & + (ydc*(xdc*ycb-zu)+xdc*(xdc*xcb+zcb*zdc))/rcbu2 365 & + 2.0d0*(xu*zdc-xdc*zu)*dphidzibu/ru2 366 dyibxic = -xcb*dphidyib/(rcb*rcb) - dphidzt - dphidzu 367 & + (xca*(zba*xcb+yt)+zca*(zba*zcb+ycb*yba))/rcbt2 368 & - 2.0d0*(yt*zba-yba*zt)*dphidyibt/rt2 369 & - (zdc*(zdb*zcb+ycb*ydb)+xdc*(zdb*xcb-yu))/rcbu2 370 & + 2.0d0*(yu*zdb-ydb*zu)*dphidyibu/ru2 371 dyibyic = -ycb*dphidyib/(rcb*rcb) 372 & - (zca*(xba*ycb+zt)-xca*(zba*ycb-xt))/rcbt2 373 & - 2.0d0*(zt*xba-zba*xt)*dphidyibt/rt2 374 & - (xdc*(zdb*ycb+xu)-zdc*(xdb*ycb-zu))/rcbu2 375 & + 2.0d0*(zu*xdb-zdb*xu)*dphidyibu/ru2 376 dyibxid = rcbxu*dphidyibu + dphidzu 377 & + (xdc*zcb*xcb+zdc*yzcb2)/rcbu2 378 dyibyid = rcbyu*dphidyibu + ycb*yu/rcbu2 379 dyibzid = rcbzu*dphidyibu - dphidxu 380 & - (zdc*xcb*zcb+xdc*xycb2)/rcbu2 381 dzibxic = -xcb*dphidzib/(rcb*rcb) + dphidyt + dphidyu 382 & - (xca*(yba*xcb-zt)+yca*(zba*zcb+ycb*yba))/rcbt2 383 & - 2.0d0*(yt*zba-yba*zt)*dphidzibt/rt2 384 & + (ydc*(zdb*zcb+ycb*ydb)+xdc*(ydb*xcb+zu))/rcbu2 385 & + 2.0d0*(yu*zdb-ydb*zu)*dphidzibu/ru2 386 dzibzic = -zcb*dphidzib/(rcb*rcb) 387 & - (xca*(yba*zcb+xt)-yca*(xba*zcb-yt))/rcbt2 388 & - 2.0d0*(xt*yba-xba*yt)*dphidzibt/rt2 389 & - (ydc*(xdb*zcb+yu)-xdc*(ydb*zcb-xu))/rcbu2 390 & + 2.0d0*(xu*ydb-xdb*yu)*dphidzibu/ru2 391 dzibxid = rcbxu*dphidzibu - dphidyu 392 & - (xdc*ycb*xcb+ydc*yzcb2)/rcbu2 393 dzibyid = rcbyu*dphidzibu + dphidxu 394 & + (ydc*xcb*ycb+xdc*xzcb2)/rcbu2 395 dzibzid = rcbzu*dphidzibu + zcb*zu/rcbu2 396 dxicxid = rcbxu*dphidxicu - xcb*(zdb*ycb-ydb*zcb)/rcbu2 397 dxicyid = rcbyu*dphidxicu + dphidzu 398 & + (ydb*zcb*ycb+zdb*xzcb2)/rcbu2 399 dxiczid = rcbzu*dphidxicu - dphidyu 400 & - (zdb*ycb*zcb+ydb*xycb2)/rcbu2 401 dyicxid = rcbxu*dphidyicu - dphidzu 402 & - (xdb*zcb*xcb+zdb*yzcb2)/rcbu2 403 dyicyid = rcbyu*dphidyicu - ycb*(xdb*zcb-zdb*xcb)/rcbu2 404 dyiczid = rcbzu*dphidyicu + dphidxu 405 & + (zdb*xcb*zcb+xdb*xycb2)/rcbu2 406 dzicxid = rcbxu*dphidzicu + dphidyu 407 & + (xdb*ycb*xcb+ydb*yzcb2)/rcbu2 408 dzicyid = rcbyu*dphidzicu - dphidxu 409 & - (ydb*xcb*ycb+xdb*xzcb2)/rcbu2 410 dziczid = rcbzu*dphidzicu - zcb*(ydb*xcb-xdb*ycb)/rcbu2 411 dxidxid = rcbxu*dphidxid 412 dxidyid = rcbxu*dphidyid + zcb*rcb/ru2 413 dxidzid = rcbxu*dphidzid - ycb*rcb/ru2 414 dyidyid = rcbyu*dphidyid 415 dyidzid = rcbyu*dphidzid + xcb*rcb/ru2 416 dzidzid = rcbzu*dphidzid 417c 418c get some second derivative chain rule terms by difference 419c 420 dxiaxib = -dxiaxia - dxiaxic - dxiaxid 421 dxiayib = -dxiayia - dxiayic - dxiayid 422 dxiazib = -dxiazia - dxiazic - dxiazid 423 dyiayib = -dyiayia - dyiayic - dyiayid 424 dyiazib = -dyiazia - dyiazic - dyiazid 425 dziazib = -dziazia - dziazic - dziazid 426 dxibxib = -dxiaxib - dxibxic - dxibxid 427 dxibyib = -dyiaxib - dxibyic - dxibyid 428 dxibzib = -dxiazib - dzibxic - dzibxid 429 dxibzic = -dziaxib - dxibzib - dxibzid 430 dyibyib = -dyiayib - dyibyic - dyibyid 431 dyibzic = -dziayib - dyibzib - dyibzid 432 dzibzib = -dziazib - dzibzic - dzibzid 433 dzibyic = -dyiazib - dyibzib - dzibyid 434 dxicxic = -dxiaxic - dxibxic - dxicxid 435 dxicyic = -dyiaxic - dyibxic - dxicyid 436 dxiczic = -dziaxic - dzibxic - dxiczid 437 dyicyic = -dyiayic - dyibyic - dyicyid 438 dyiczic = -dziayic - dzibyic - dyiczid 439 dziczic = -dziazic - dzibzic - dziczid 440c 441c increment diagonal and off-diagonal Hessian elements 442c 443 if (i .eq. ia) then 444 hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxia 445 & + d2edphi2*dphidxia*dphidxia 446 hessy(1,ia) = hessy(1,ia) + dedphi*dxiayia 447 & + d2edphi2*dphidxia*dphidyia 448 hessz(1,ia) = hessz(1,ia) + dedphi*dxiazia 449 & + d2edphi2*dphidxia*dphidzia 450 hessx(2,ia) = hessx(2,ia) + dedphi*dxiayia 451 & + d2edphi2*dphidxia*dphidyia 452 hessy(2,ia) = hessy(2,ia) + dedphi*dyiayia 453 & + d2edphi2*dphidyia*dphidyia 454 hessz(2,ia) = hessz(2,ia) + dedphi*dyiazia 455 & + d2edphi2*dphidyia*dphidzia 456 hessx(3,ia) = hessx(3,ia) + dedphi*dxiazia 457 & + d2edphi2*dphidxia*dphidzia 458 hessy(3,ia) = hessy(3,ia) + dedphi*dyiazia 459 & + d2edphi2*dphidyia*dphidzia 460 hessz(3,ia) = hessz(3,ia) + dedphi*dziazia 461 & + d2edphi2*dphidzia*dphidzia 462 hessx(1,ib) = hessx(1,ib) + dedphi*dxiaxib 463 & + d2edphi2*dphidxia*dphidxib 464 hessy(1,ib) = hessy(1,ib) + dedphi*dyiaxib 465 & + d2edphi2*dphidyia*dphidxib 466 hessz(1,ib) = hessz(1,ib) + dedphi*dziaxib 467 & + d2edphi2*dphidzia*dphidxib 468 hessx(2,ib) = hessx(2,ib) + dedphi*dxiayib 469 & + d2edphi2*dphidxia*dphidyib 470 hessy(2,ib) = hessy(2,ib) + dedphi*dyiayib 471 & + d2edphi2*dphidyia*dphidyib 472 hessz(2,ib) = hessz(2,ib) + dedphi*dziayib 473 & + d2edphi2*dphidzia*dphidyib 474 hessx(3,ib) = hessx(3,ib) + dedphi*dxiazib 475 & + d2edphi2*dphidxia*dphidzib 476 hessy(3,ib) = hessy(3,ib) + dedphi*dyiazib 477 & + d2edphi2*dphidyia*dphidzib 478 hessz(3,ib) = hessz(3,ib) + dedphi*dziazib 479 & + d2edphi2*dphidzia*dphidzib 480 hessx(1,ic) = hessx(1,ic) + dedphi*dxiaxic 481 & + d2edphi2*dphidxia*dphidxic 482 hessy(1,ic) = hessy(1,ic) + dedphi*dyiaxic 483 & + d2edphi2*dphidyia*dphidxic 484 hessz(1,ic) = hessz(1,ic) + dedphi*dziaxic 485 & + d2edphi2*dphidzia*dphidxic 486 hessx(2,ic) = hessx(2,ic) + dedphi*dxiayic 487 & + d2edphi2*dphidxia*dphidyic 488 hessy(2,ic) = hessy(2,ic) + dedphi*dyiayic 489 & + d2edphi2*dphidyia*dphidyic 490 hessz(2,ic) = hessz(2,ic) + dedphi*dziayic 491 & + d2edphi2*dphidzia*dphidyic 492 hessx(3,ic) = hessx(3,ic) + dedphi*dxiazic 493 & + d2edphi2*dphidxia*dphidzic 494 hessy(3,ic) = hessy(3,ic) + dedphi*dyiazic 495 & + d2edphi2*dphidyia*dphidzic 496 hessz(3,ic) = hessz(3,ic) + dedphi*dziazic 497 & + d2edphi2*dphidzia*dphidzic 498 hessx(1,id) = hessx(1,id) + dedphi*dxiaxid 499 & + d2edphi2*dphidxia*dphidxid 500 hessy(1,id) = hessy(1,id) + dedphi*dyiaxid 501 & + d2edphi2*dphidyia*dphidxid 502 hessz(1,id) = hessz(1,id) + dedphi*dziaxid 503 & + d2edphi2*dphidzia*dphidxid 504 hessx(2,id) = hessx(2,id) + dedphi*dxiayid 505 & + d2edphi2*dphidxia*dphidyid 506 hessy(2,id) = hessy(2,id) + dedphi*dyiayid 507 & + d2edphi2*dphidyia*dphidyid 508 hessz(2,id) = hessz(2,id) + dedphi*dziayid 509 & + d2edphi2*dphidzia*dphidyid 510 hessx(3,id) = hessx(3,id) + dedphi*dxiazid 511 & + d2edphi2*dphidxia*dphidzid 512 hessy(3,id) = hessy(3,id) + dedphi*dyiazid 513 & + d2edphi2*dphidyia*dphidzid 514 hessz(3,id) = hessz(3,id) + dedphi*dziazid 515 & + d2edphi2*dphidzia*dphidzid 516 else if (i .eq. ib) then 517 hessx(1,ib) = hessx(1,ib) + dedphi*dxibxib 518 & + d2edphi2*dphidxib*dphidxib 519 hessy(1,ib) = hessy(1,ib) + dedphi*dxibyib 520 & + d2edphi2*dphidxib*dphidyib 521 hessz(1,ib) = hessz(1,ib) + dedphi*dxibzib 522 & + d2edphi2*dphidxib*dphidzib 523 hessx(2,ib) = hessx(2,ib) + dedphi*dxibyib 524 & + d2edphi2*dphidxib*dphidyib 525 hessy(2,ib) = hessy(2,ib) + dedphi*dyibyib 526 & + d2edphi2*dphidyib*dphidyib 527 hessz(2,ib) = hessz(2,ib) + dedphi*dyibzib 528 & + d2edphi2*dphidyib*dphidzib 529 hessx(3,ib) = hessx(3,ib) + dedphi*dxibzib 530 & + d2edphi2*dphidxib*dphidzib 531 hessy(3,ib) = hessy(3,ib) + dedphi*dyibzib 532 & + d2edphi2*dphidyib*dphidzib 533 hessz(3,ib) = hessz(3,ib) + dedphi*dzibzib 534 & + d2edphi2*dphidzib*dphidzib 535 hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxib 536 & + d2edphi2*dphidxib*dphidxia 537 hessy(1,ia) = hessy(1,ia) + dedphi*dxiayib 538 & + d2edphi2*dphidyib*dphidxia 539 hessz(1,ia) = hessz(1,ia) + dedphi*dxiazib 540 & + d2edphi2*dphidzib*dphidxia 541 hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxib 542 & + d2edphi2*dphidxib*dphidyia 543 hessy(2,ia) = hessy(2,ia) + dedphi*dyiayib 544 & + d2edphi2*dphidyib*dphidyia 545 hessz(2,ia) = hessz(2,ia) + dedphi*dyiazib 546 & + d2edphi2*dphidzib*dphidyia 547 hessx(3,ia) = hessx(3,ia) + dedphi*dziaxib 548 & + d2edphi2*dphidxib*dphidzia 549 hessy(3,ia) = hessy(3,ia) + dedphi*dziayib 550 & + d2edphi2*dphidyib*dphidzia 551 hessz(3,ia) = hessz(3,ia) + dedphi*dziazib 552 & + d2edphi2*dphidzib*dphidzia 553 hessx(1,ic) = hessx(1,ic) + dedphi*dxibxic 554 & + d2edphi2*dphidxib*dphidxic 555 hessy(1,ic) = hessy(1,ic) + dedphi*dyibxic 556 & + d2edphi2*dphidyib*dphidxic 557 hessz(1,ic) = hessz(1,ic) + dedphi*dzibxic 558 & + d2edphi2*dphidzib*dphidxic 559 hessx(2,ic) = hessx(2,ic) + dedphi*dxibyic 560 & + d2edphi2*dphidxib*dphidyic 561 hessy(2,ic) = hessy(2,ic) + dedphi*dyibyic 562 & + d2edphi2*dphidyib*dphidyic 563 hessz(2,ic) = hessz(2,ic) + dedphi*dzibyic 564 & + d2edphi2*dphidzib*dphidyic 565 hessx(3,ic) = hessx(3,ic) + dedphi*dxibzic 566 & + d2edphi2*dphidxib*dphidzic 567 hessy(3,ic) = hessy(3,ic) + dedphi*dyibzic 568 & + d2edphi2*dphidyib*dphidzic 569 hessz(3,ic) = hessz(3,ic) + dedphi*dzibzic 570 & + d2edphi2*dphidzib*dphidzic 571 hessx(1,id) = hessx(1,id) + dedphi*dxibxid 572 & + d2edphi2*dphidxib*dphidxid 573 hessy(1,id) = hessy(1,id) + dedphi*dyibxid 574 & + d2edphi2*dphidyib*dphidxid 575 hessz(1,id) = hessz(1,id) + dedphi*dzibxid 576 & + d2edphi2*dphidzib*dphidxid 577 hessx(2,id) = hessx(2,id) + dedphi*dxibyid 578 & + d2edphi2*dphidxib*dphidyid 579 hessy(2,id) = hessy(2,id) + dedphi*dyibyid 580 & + d2edphi2*dphidyib*dphidyid 581 hessz(2,id) = hessz(2,id) + dedphi*dzibyid 582 & + d2edphi2*dphidzib*dphidyid 583 hessx(3,id) = hessx(3,id) + dedphi*dxibzid 584 & + d2edphi2*dphidxib*dphidzid 585 hessy(3,id) = hessy(3,id) + dedphi*dyibzid 586 & + d2edphi2*dphidyib*dphidzid 587 hessz(3,id) = hessz(3,id) + dedphi*dzibzid 588 & + d2edphi2*dphidzib*dphidzid 589 else if (i .eq. ic) then 590 hessx(1,ic) = hessx(1,ic) + dedphi*dxicxic 591 & + d2edphi2*dphidxic*dphidxic 592 hessy(1,ic) = hessy(1,ic) + dedphi*dxicyic 593 & + d2edphi2*dphidxic*dphidyic 594 hessz(1,ic) = hessz(1,ic) + dedphi*dxiczic 595 & + d2edphi2*dphidxic*dphidzic 596 hessx(2,ic) = hessx(2,ic) + dedphi*dxicyic 597 & + d2edphi2*dphidxic*dphidyic 598 hessy(2,ic) = hessy(2,ic) + dedphi*dyicyic 599 & + d2edphi2*dphidyic*dphidyic 600 hessz(2,ic) = hessz(2,ic) + dedphi*dyiczic 601 & + d2edphi2*dphidyic*dphidzic 602 hessx(3,ic) = hessx(3,ic) + dedphi*dxiczic 603 & + d2edphi2*dphidxic*dphidzic 604 hessy(3,ic) = hessy(3,ic) + dedphi*dyiczic 605 & + d2edphi2*dphidyic*dphidzic 606 hessz(3,ic) = hessz(3,ic) + dedphi*dziczic 607 & + d2edphi2*dphidzic*dphidzic 608 hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxic 609 & + d2edphi2*dphidxic*dphidxia 610 hessy(1,ia) = hessy(1,ia) + dedphi*dxiayic 611 & + d2edphi2*dphidyic*dphidxia 612 hessz(1,ia) = hessz(1,ia) + dedphi*dxiazic 613 & + d2edphi2*dphidzic*dphidxia 614 hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxic 615 & + d2edphi2*dphidxic*dphidyia 616 hessy(2,ia) = hessy(2,ia) + dedphi*dyiayic 617 & + d2edphi2*dphidyic*dphidyia 618 hessz(2,ia) = hessz(2,ia) + dedphi*dyiazic 619 & + d2edphi2*dphidzic*dphidyia 620 hessx(3,ia) = hessx(3,ia) + dedphi*dziaxic 621 & + d2edphi2*dphidxic*dphidzia 622 hessy(3,ia) = hessy(3,ia) + dedphi*dziayic 623 & + d2edphi2*dphidyic*dphidzia 624 hessz(3,ia) = hessz(3,ia) + dedphi*dziazic 625 & + d2edphi2*dphidzic*dphidzia 626 hessx(1,ib) = hessx(1,ib) + dedphi*dxibxic 627 & + d2edphi2*dphidxic*dphidxib 628 hessy(1,ib) = hessy(1,ib) + dedphi*dxibyic 629 & + d2edphi2*dphidyic*dphidxib 630 hessz(1,ib) = hessz(1,ib) + dedphi*dxibzic 631 & + d2edphi2*dphidzic*dphidxib 632 hessx(2,ib) = hessx(2,ib) + dedphi*dyibxic 633 & + d2edphi2*dphidxic*dphidyib 634 hessy(2,ib) = hessy(2,ib) + dedphi*dyibyic 635 & + d2edphi2*dphidyic*dphidyib 636 hessz(2,ib) = hessz(2,ib) + dedphi*dyibzic 637 & + d2edphi2*dphidzic*dphidyib 638 hessx(3,ib) = hessx(3,ib) + dedphi*dzibxic 639 & + d2edphi2*dphidxic*dphidzib 640 hessy(3,ib) = hessy(3,ib) + dedphi*dzibyic 641 & + d2edphi2*dphidyic*dphidzib 642 hessz(3,ib) = hessz(3,ib) + dedphi*dzibzic 643 & + d2edphi2*dphidzic*dphidzib 644 hessx(1,id) = hessx(1,id) + dedphi*dxicxid 645 & + d2edphi2*dphidxic*dphidxid 646 hessy(1,id) = hessy(1,id) + dedphi*dyicxid 647 & + d2edphi2*dphidyic*dphidxid 648 hessz(1,id) = hessz(1,id) + dedphi*dzicxid 649 & + d2edphi2*dphidzic*dphidxid 650 hessx(2,id) = hessx(2,id) + dedphi*dxicyid 651 & + d2edphi2*dphidxic*dphidyid 652 hessy(2,id) = hessy(2,id) + dedphi*dyicyid 653 & + d2edphi2*dphidyic*dphidyid 654 hessz(2,id) = hessz(2,id) + dedphi*dzicyid 655 & + d2edphi2*dphidzic*dphidyid 656 hessx(3,id) = hessx(3,id) + dedphi*dxiczid 657 & + d2edphi2*dphidxic*dphidzid 658 hessy(3,id) = hessy(3,id) + dedphi*dyiczid 659 & + d2edphi2*dphidyic*dphidzid 660 hessz(3,id) = hessz(3,id) + dedphi*dziczid 661 & + d2edphi2*dphidzic*dphidzid 662 else if (i .eq. id) then 663 hessx(1,id) = hessx(1,id) + dedphi*dxidxid 664 & + d2edphi2*dphidxid*dphidxid 665 hessy(1,id) = hessy(1,id) + dedphi*dxidyid 666 & + d2edphi2*dphidxid*dphidyid 667 hessz(1,id) = hessz(1,id) + dedphi*dxidzid 668 & + d2edphi2*dphidxid*dphidzid 669 hessx(2,id) = hessx(2,id) + dedphi*dxidyid 670 & + d2edphi2*dphidxid*dphidyid 671 hessy(2,id) = hessy(2,id) + dedphi*dyidyid 672 & + d2edphi2*dphidyid*dphidyid 673 hessz(2,id) = hessz(2,id) + dedphi*dyidzid 674 & + d2edphi2*dphidyid*dphidzid 675 hessx(3,id) = hessx(3,id) + dedphi*dxidzid 676 & + d2edphi2*dphidxid*dphidzid 677 hessy(3,id) = hessy(3,id) + dedphi*dyidzid 678 & + d2edphi2*dphidyid*dphidzid 679 hessz(3,id) = hessz(3,id) + dedphi*dzidzid 680 & + d2edphi2*dphidzid*dphidzid 681 hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxid 682 & + d2edphi2*dphidxid*dphidxia 683 hessy(1,ia) = hessy(1,ia) + dedphi*dxiayid 684 & + d2edphi2*dphidyid*dphidxia 685 hessz(1,ia) = hessz(1,ia) + dedphi*dxiazid 686 & + d2edphi2*dphidzid*dphidxia 687 hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxid 688 & + d2edphi2*dphidxid*dphidyia 689 hessy(2,ia) = hessy(2,ia) + dedphi*dyiayid 690 & + d2edphi2*dphidyid*dphidyia 691 hessz(2,ia) = hessz(2,ia) + dedphi*dyiazid 692 & + d2edphi2*dphidzid*dphidyia 693 hessx(3,ia) = hessx(3,ia) + dedphi*dziaxid 694 & + d2edphi2*dphidxid*dphidzia 695 hessy(3,ia) = hessy(3,ia) + dedphi*dziayid 696 & + d2edphi2*dphidyid*dphidzia 697 hessz(3,ia) = hessz(3,ia) + dedphi*dziazid 698 & + d2edphi2*dphidzid*dphidzia 699 hessx(1,ib) = hessx(1,ib) + dedphi*dxibxid 700 & + d2edphi2*dphidxid*dphidxib 701 hessy(1,ib) = hessy(1,ib) + dedphi*dxibyid 702 & + d2edphi2*dphidyid*dphidxib 703 hessz(1,ib) = hessz(1,ib) + dedphi*dxibzid 704 & + d2edphi2*dphidzid*dphidxib 705 hessx(2,ib) = hessx(2,ib) + dedphi*dyibxid 706 & + d2edphi2*dphidxid*dphidyib 707 hessy(2,ib) = hessy(2,ib) + dedphi*dyibyid 708 & + d2edphi2*dphidyid*dphidyib 709 hessz(2,ib) = hessz(2,ib) + dedphi*dyibzid 710 & + d2edphi2*dphidzid*dphidyib 711 hessx(3,ib) = hessx(3,ib) + dedphi*dzibxid 712 & + d2edphi2*dphidxid*dphidzib 713 hessy(3,ib) = hessy(3,ib) + dedphi*dzibyid 714 & + d2edphi2*dphidyid*dphidzib 715 hessz(3,ib) = hessz(3,ib) + dedphi*dzibzid 716 & + d2edphi2*dphidzid*dphidzib 717 hessx(1,ic) = hessx(1,ic) + dedphi*dxicxid 718 & + d2edphi2*dphidxid*dphidxic 719 hessy(1,ic) = hessy(1,ic) + dedphi*dxicyid 720 & + d2edphi2*dphidyid*dphidxic 721 hessz(1,ic) = hessz(1,ic) + dedphi*dxiczid 722 & + d2edphi2*dphidzid*dphidxic 723 hessx(2,ic) = hessx(2,ic) + dedphi*dyicxid 724 & + d2edphi2*dphidxid*dphidyic 725 hessy(2,ic) = hessy(2,ic) + dedphi*dyicyid 726 & + d2edphi2*dphidyid*dphidyic 727 hessz(2,ic) = hessz(2,ic) + dedphi*dyiczid 728 & + d2edphi2*dphidzid*dphidyic 729 hessx(3,ic) = hessx(3,ic) + dedphi*dzicxid 730 & + d2edphi2*dphidxid*dphidzic 731 hessy(3,ic) = hessy(3,ic) + dedphi*dzicyid 732 & + d2edphi2*dphidyid*dphidzic 733 hessz(3,ic) = hessz(3,ic) + dedphi*dziczid 734 & + d2edphi2*dphidzid*dphidzic 735 end if 736 end if 737 end do 738 return 739 end 740c 741c 742c ################################################################ 743c ## ## 744c ## subroutine etors2b -- smoothed torsional angle Hessian ## 745c ## ## 746c ################################################################ 747c 748c 749c "etors2b" calculates the second derivatives of the torsional 750c energy for a single atom for use with potential energy 751c smoothing methods 752c 753c 754 subroutine etors2b (i) 755 use atoms 756 use group 757 use hessn 758 use math 759 use torpot 760 use tors 761 use warp 762 implicit none 763 integer i,ktors 764 integer ia,ib,ic,id 765 real*8 eps,fgrp 766 real*8 width,wterm 767 real*8 dedphi,d2edphi2 768 real*8 v1,v2,v3,v4,v5,v6 769 real*8 c1,c2,c3,c4,c5,c6 770 real*8 s1,s2,s3,s4,s5,s6 771 real*8 cosine,cosine2 772 real*8 cosine3,cosine4 773 real*8 cosine5,cosine6 774 real*8 sine,sine2,sine3 775 real*8 sine4,sine5,sine6 776 real*8 damp1,damp2,damp3 777 real*8 damp4,damp5,damp6 778 real*8 xia,yia,zia 779 real*8 xib,yib,zib 780 real*8 xic,yic,zic 781 real*8 xid,yid,zid 782 real*8 xba,yba,zba 783 real*8 xcb,ycb,zcb 784 real*8 xdc,ydc,zdc 785 real*8 xca,yca,zca 786 real*8 xdb,ydb,zdb 787 real*8 xt,yt,zt,xu,yu,zu 788 real*8 xtu,ytu,ztu 789 real*8 rt2,ru2,rtru,rcb 790 real*8 dphi1,dphi2,dphi3 791 real*8 dphi4,dphi5,dphi6 792 real*8 d2phi1,d2phi2,d2phi3 793 real*8 d2phi4,d2phi5,d2phi6 794 real*8 dphidxt,dphidyt,dphidzt 795 real*8 dphidxu,dphidyu,dphidzu 796 real*8 dphidxia,dphidyia,dphidzia 797 real*8 dphidxib,dphidyib,dphidzib 798 real*8 dphidxic,dphidyic,dphidzic 799 real*8 dphidxid,dphidyid,dphidzid 800 real*8 xycb2,xzcb2,yzcb2 801 real*8 rcbxt,rcbyt,rcbzt,rcbt2 802 real*8 rcbxu,rcbyu,rcbzu,rcbu2 803 real*8 dphidxibt,dphidyibt,dphidzibt 804 real*8 dphidxibu,dphidyibu,dphidzibu 805 real*8 dphidxict,dphidyict,dphidzict 806 real*8 dphidxicu,dphidyicu,dphidzicu 807 real*8 dxiaxia,dyiayia,dziazia 808 real*8 dxibxib,dyibyib,dzibzib 809 real*8 dxicxic,dyicyic,dziczic 810 real*8 dxidxid,dyidyid,dzidzid 811 real*8 dxiayia,dxiazia,dyiazia 812 real*8 dxibyib,dxibzib,dyibzib 813 real*8 dxicyic,dxiczic,dyiczic 814 real*8 dxidyid,dxidzid,dyidzid 815 real*8 dxiaxib,dxiayib,dxiazib 816 real*8 dyiaxib,dyiayib,dyiazib 817 real*8 dziaxib,dziayib,dziazib 818 real*8 dxiaxic,dxiayic,dxiazic 819 real*8 dyiaxic,dyiayic,dyiazic 820 real*8 dziaxic,dziayic,dziazic 821 real*8 dxiaxid,dxiayid,dxiazid 822 real*8 dyiaxid,dyiayid,dyiazid 823 real*8 dziaxid,dziayid,dziazid 824 real*8 dxibxic,dxibyic,dxibzic 825 real*8 dyibxic,dyibyic,dyibzic 826 real*8 dzibxic,dzibyic,dzibzic 827 real*8 dxibxid,dxibyid,dxibzid 828 real*8 dyibxid,dyibyid,dyibzid 829 real*8 dzibxid,dzibyid,dzibzid 830 real*8 dxicxid,dxicyid,dxiczid 831 real*8 dyicxid,dyicyid,dyiczid 832 real*8 dzicxid,dzicyid,dziczid 833 logical proceed 834c 835c 836c set tolerance for minimum distance and angle values 837c 838 eps = 0.0001d0 839c 840c set the extent of smoothing to be performed 841c 842 width = difft * deform 843 if (width .le. 0.0d0) then 844 damp1 = 1.0d0 845 damp2 = 1.0d0 846 damp3 = 1.0d0 847 damp4 = 1.0d0 848 damp5 = 1.0d0 849 damp6 = 1.0d0 850 else if (use_dem) then 851 damp1 = exp(-width) 852 damp2 = exp(-4.0d0*width) 853 damp3 = exp(-9.0d0*width) 854 damp4 = exp(-16.0d0*width) 855 damp5 = exp(-25.0d0*width) 856 damp6 = exp(-36.0d0*width) 857 else if (use_gda) then 858 wterm = difft / 12.0d0 859 else if (use_tophat .or. use_stophat) then 860 damp1 = 0.0d0 861 damp2 = 0.0d0 862 damp3 = 0.0d0 863 damp4 = 0.0d0 864 damp5 = 0.0d0 865 damp6 = 0.0d0 866 if (width .lt. pi) damp1 = sin(width) / width 867 wterm = 2.0d0 * width 868 if (wterm .lt. pi) damp2 = sin(wterm) / wterm 869 wterm = 3.0d0 * width 870 if (wterm .lt. pi) damp3 = sin(wterm) / wterm 871 wterm = 4.0d0 * width 872 if (wterm .lt. pi) damp4 = sin(wterm) / wterm 873 wterm = 5.0d0 * width 874 if (wterm .lt. pi) damp5 = sin(wterm) / wterm 875 wterm = 6.0d0 * width 876 if (wterm .lt. pi) damp6 = sin(wterm) / wterm 877 end if 878c 879c calculate the torsional angle energy term 880c 881 do ktors = 1, ntors 882 ia = itors(1,ktors) 883 ib = itors(2,ktors) 884 ic = itors(3,ktors) 885 id = itors(4,ktors) 886c 887c decide whether to compute the current interaction 888c 889 proceed = (i.eq.ia .or. i.eq.ib .or. i.eq.ic .or. i.eq.id) 890 if (proceed .and. use_group) 891 & call groups (proceed,fgrp,ia,ib,ic,id,0,0) 892c 893c compute the value of the torsional angle 894c 895 if (proceed) then 896 xia = x(ia) 897 yia = y(ia) 898 zia = z(ia) 899 xib = x(ib) 900 yib = y(ib) 901 zib = z(ib) 902 xic = x(ic) 903 yic = y(ic) 904 zic = z(ic) 905 xid = x(id) 906 yid = y(id) 907 zid = z(id) 908 xba = xib - xia 909 yba = yib - yia 910 zba = zib - zia 911 xcb = xic - xib 912 ycb = yic - yib 913 zcb = zic - zib 914 xdc = xid - xic 915 ydc = yid - yic 916 zdc = zid - zic 917 rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps)) 918 xt = yba*zcb - ycb*zba 919 yt = zba*xcb - zcb*xba 920 zt = xba*ycb - xcb*yba 921 xu = ycb*zdc - ydc*zcb 922 yu = zcb*xdc - zdc*xcb 923 zu = xcb*ydc - xdc*ycb 924 xtu = yt*zu - yu*zt 925 ytu = zt*xu - zu*xt 926 ztu = xt*yu - xu*yt 927 rt2 = max(xt*xt+yt*yt+zt*zt,eps) 928 ru2 = max(xu*xu+yu*yu+zu*zu,eps) 929 rtru = sqrt(rt2 * ru2) 930 cosine = (xt*xu + yt*yu + zt*zu) / rtru 931 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru) 932c 933c set the torsional parameters for this angle 934c 935 v1 = tors1(1,ktors) 936 c1 = tors1(3,ktors) 937 s1 = tors1(4,ktors) 938 v2 = tors2(1,ktors) 939 c2 = tors2(3,ktors) 940 s2 = tors2(4,ktors) 941 v3 = tors3(1,ktors) 942 c3 = tors3(3,ktors) 943 s3 = tors3(4,ktors) 944 v4 = tors4(1,ktors) 945 c4 = tors4(3,ktors) 946 s4 = tors4(4,ktors) 947 v5 = tors5(1,ktors) 948 c5 = tors5(3,ktors) 949 s5 = tors5(4,ktors) 950 v6 = tors6(1,ktors) 951 c6 = tors6(3,ktors) 952 s6 = tors6(4,ktors) 953c 954c compute the multiple angle trigonometry and the phase terms 955c 956 cosine2 = cosine*cosine - sine*sine 957 sine2 = 2.0d0 * cosine * sine 958 cosine3 = cosine*cosine2 - sine*sine2 959 sine3 = cosine*sine2 + sine*cosine2 960 cosine4 = cosine*cosine3 - sine*sine3 961 sine4 = cosine*sine3 + sine*cosine3 962 cosine5 = cosine*cosine4 - sine*sine4 963 sine5 = cosine*sine4 + sine*cosine4 964 cosine6 = cosine*cosine5 - sine*sine5 965 sine6 = cosine*sine5 + sine*cosine5 966 dphi1 = (cosine*s1 - sine*c1) 967 dphi2 = 2.0d0 * (cosine2*s2 - sine2*c2) 968 dphi3 = 3.0d0 * (cosine3*s3 - sine3*c3) 969 dphi4 = 4.0d0 * (cosine4*s4 - sine4*c4) 970 dphi5 = 5.0d0 * (cosine5*s5 - sine5*c5) 971 dphi6 = 6.0d0 * (cosine6*s6 - sine6*c6) 972 d2phi1 = -(cosine*c1 + sine*s1) 973 d2phi2 = -4.0d0 * (cosine2*c2 + sine2*s2) 974 d2phi3 = -9.0d0 * (cosine3*c3 + sine3*s3) 975 d2phi4 = -16.0d0 * (cosine4*c4 + sine4*s4) 976 d2phi5 = -25.0d0 * (cosine5*c5 + sine5*s5) 977 d2phi6 = -36.0d0 * (cosine6*c6 + sine6*s6) 978c 979c transform the potential function via smoothing 980c 981 if (use_gda) then 982 width = wterm * (m2(ia)+m2(ib)+m2(ic)+m2(id)) 983 damp1 = exp(-width) 984 damp2 = exp(-4.0d0*width) 985 damp3 = exp(-9.0d0*width) 986 damp4 = exp(-16.0d0*width) 987 damp5 = exp(-25.0d0*width) 988 damp6 = exp(-36.0d0*width) 989 end if 990 dphi1 = dphi1 * damp1 991 dphi2 = dphi2 * damp2 992 dphi3 = dphi3 * damp3 993 dphi4 = dphi4 * damp4 994 dphi5 = dphi5 * damp5 995 dphi6 = dphi6 * damp6 996 d2phi1 = d2phi1 * damp1 997 d2phi2 = d2phi2 * damp2 998 d2phi3 = d2phi3 * damp3 999 d2phi4 = d2phi4 * damp4 1000 d2phi5 = d2phi5 * damp5 1001 d2phi6 = d2phi6 * damp6 1002c 1003c calculate the torsional master chain rule terms 1004c 1005 dedphi = torsunit * (v1*dphi1 + v2*dphi2 + v3*dphi3 1006 & + v4*dphi4 + v5*dphi5 + v6*dphi6) 1007 d2edphi2 = torsunit * (v1*d2phi1 + v2*d2phi2 + v3*d2phi3 1008 & + v4*d2phi4 + v5*d2phi5 + v6*d2phi6) 1009c 1010c scale the interaction based on its group membership 1011c 1012 if (use_group) then 1013 dedphi = dedphi * fgrp 1014 d2edphi2 = d2edphi2 * fgrp 1015 end if 1016c 1017c abbreviations for first derivative chain rule terms 1018c 1019 xca = xic - xia 1020 yca = yic - yia 1021 zca = zic - zia 1022 xdb = xid - xib 1023 ydb = yid - yib 1024 zdb = zid - zib 1025 dphidxt = (yt*zcb - ycb*zt) / (rt2*rcb) 1026 dphidyt = (zt*xcb - zcb*xt) / (rt2*rcb) 1027 dphidzt = (xt*ycb - xcb*yt) / (rt2*rcb) 1028 dphidxu = -(yu*zcb - ycb*zu) / (ru2*rcb) 1029 dphidyu = -(zu*xcb - zcb*xu) / (ru2*rcb) 1030 dphidzu = -(xu*ycb - xcb*yu) / (ru2*rcb) 1031c 1032c abbreviations for second derivative chain rule terms 1033c 1034 xycb2 = xcb*xcb + ycb*ycb 1035 xzcb2 = xcb*xcb + zcb*zcb 1036 yzcb2 = ycb*ycb + zcb*zcb 1037 rcbxt = -2.0d0 * rcb * dphidxt 1038 rcbyt = -2.0d0 * rcb * dphidyt 1039 rcbzt = -2.0d0 * rcb * dphidzt 1040 rcbt2 = rcb * rt2 1041 rcbxu = 2.0d0 * rcb * dphidxu 1042 rcbyu = 2.0d0 * rcb * dphidyu 1043 rcbzu = 2.0d0 * rcb * dphidzu 1044 rcbu2 = rcb * ru2 1045 dphidxibt = yca*dphidzt - zca*dphidyt 1046 dphidxibu = zdc*dphidyu - ydc*dphidzu 1047 dphidyibt = zca*dphidxt - xca*dphidzt 1048 dphidyibu = xdc*dphidzu - zdc*dphidxu 1049 dphidzibt = xca*dphidyt - yca*dphidxt 1050 dphidzibu = ydc*dphidxu - xdc*dphidyu 1051 dphidxict = zba*dphidyt - yba*dphidzt 1052 dphidxicu = ydb*dphidzu - zdb*dphidyu 1053 dphidyict = xba*dphidzt - zba*dphidxt 1054 dphidyicu = zdb*dphidxu - xdb*dphidzu 1055 dphidzict = yba*dphidxt - xba*dphidyt 1056 dphidzicu = xdb*dphidyu - ydb*dphidxu 1057c 1058c chain rule terms for first derivative components 1059c 1060 dphidxia = zcb*dphidyt - ycb*dphidzt 1061 dphidyia = xcb*dphidzt - zcb*dphidxt 1062 dphidzia = ycb*dphidxt - xcb*dphidyt 1063 dphidxib = dphidxibt + dphidxibu 1064 dphidyib = dphidyibt + dphidyibu 1065 dphidzib = dphidzibt + dphidzibu 1066 dphidxic = dphidxict + dphidxicu 1067 dphidyic = dphidyict + dphidyicu 1068 dphidzic = dphidzict + dphidzicu 1069 dphidxid = zcb*dphidyu - ycb*dphidzu 1070 dphidyid = xcb*dphidzu - zcb*dphidxu 1071 dphidzid = ycb*dphidxu - xcb*dphidyu 1072c 1073c chain rule terms for second derivative components 1074c 1075 dxiaxia = rcbxt*dphidxia 1076 dxiayia = rcbxt*dphidyia - zcb*rcb/rt2 1077 dxiazia = rcbxt*dphidzia + ycb*rcb/rt2 1078 dxiaxic = rcbxt*dphidxict + xcb*xt/rcbt2 1079 dxiayic = rcbxt*dphidyict - dphidzt 1080 & - (xba*zcb*xcb+zba*yzcb2)/rcbt2 1081 dxiazic = rcbxt*dphidzict + dphidyt 1082 & + (xba*ycb*xcb+yba*yzcb2)/rcbt2 1083 dxiaxid = 0.0d0 1084 dxiayid = 0.0d0 1085 dxiazid = 0.0d0 1086 dyiayia = rcbyt*dphidyia 1087 dyiazia = rcbyt*dphidzia - xcb*rcb/rt2 1088 dyiaxib = rcbyt*dphidxibt - dphidzt 1089 & - (yca*zcb*ycb+zca*xzcb2)/rcbt2 1090 dyiaxic = rcbyt*dphidxict + dphidzt 1091 & + (yba*zcb*ycb+zba*xzcb2)/rcbt2 1092 dyiayic = rcbyt*dphidyict + ycb*yt/rcbt2 1093 dyiazic = rcbyt*dphidzict - dphidxt 1094 & - (yba*xcb*ycb+xba*xzcb2)/rcbt2 1095 dyiaxid = 0.0d0 1096 dyiayid = 0.0d0 1097 dyiazid = 0.0d0 1098 dziazia = rcbzt*dphidzia 1099 dziaxib = rcbzt*dphidxibt + dphidyt 1100 & + (zca*ycb*zcb+yca*xycb2)/rcbt2 1101 dziayib = rcbzt*dphidyibt - dphidxt 1102 & - (zca*xcb*zcb+xca*xycb2)/rcbt2 1103 dziaxic = rcbzt*dphidxict - dphidyt 1104 & - (zba*ycb*zcb+yba*xycb2)/rcbt2 1105 dziayic = rcbzt*dphidyict + dphidxt 1106 & + (zba*xcb*zcb+xba*xycb2)/rcbt2 1107 dziazic = rcbzt*dphidzict + zcb*zt/rcbt2 1108 dziaxid = 0.0d0 1109 dziayid = 0.0d0 1110 dziazid = 0.0d0 1111 dxibxic = -xcb*dphidxib/(rcb*rcb) 1112 & - (yca*(zba*xcb+yt)-zca*(yba*xcb-zt))/rcbt2 1113 & - 2.0d0*(yt*zba-yba*zt)*dphidxibt/rt2 1114 & - (zdc*(ydb*xcb+zu)-ydc*(zdb*xcb-yu))/rcbu2 1115 & + 2.0d0*(yu*zdb-ydb*zu)*dphidxibu/ru2 1116 dxibyic = -ycb*dphidxib/(rcb*rcb) + dphidzt + dphidzu 1117 & - (yca*(zba*ycb-xt)+zca*(xba*xcb+zcb*zba))/rcbt2 1118 & - 2.0d0*(zt*xba-zba*xt)*dphidxibt/rt2 1119 & + (zdc*(xdb*xcb+zcb*zdb)+ydc*(zdb*ycb+xu))/rcbu2 1120 & + 2.0d0*(zu*xdb-zdb*xu)*dphidxibu/ru2 1121 dxibxid = rcbxu*dphidxibu + xcb*xu/rcbu2 1122 dxibyid = rcbyu*dphidxibu - dphidzu 1123 & - (ydc*zcb*ycb+zdc*xzcb2)/rcbu2 1124 dxibzid = rcbzu*dphidxibu + dphidyu 1125 & + (zdc*ycb*zcb+ydc*xycb2)/rcbu2 1126 dyibzib = ycb*dphidzib/(rcb*rcb) 1127 & - (xca*(xca*xcb+zcb*zca)+yca*(ycb*xca+zt))/rcbt2 1128 & - 2.0d0*(xt*zca-xca*zt)*dphidzibt/rt2 1129 & + (ydc*(xdc*ycb-zu)+xdc*(xdc*xcb+zcb*zdc))/rcbu2 1130 & + 2.0d0*(xu*zdc-xdc*zu)*dphidzibu/ru2 1131 dyibxic = -xcb*dphidyib/(rcb*rcb) - dphidzt - dphidzu 1132 & + (xca*(zba*xcb+yt)+zca*(zba*zcb+ycb*yba))/rcbt2 1133 & - 2.0d0*(yt*zba-yba*zt)*dphidyibt/rt2 1134 & - (zdc*(zdb*zcb+ycb*ydb)+xdc*(zdb*xcb-yu))/rcbu2 1135 & + 2.0d0*(yu*zdb-ydb*zu)*dphidyibu/ru2 1136 dyibyic = -ycb*dphidyib/(rcb*rcb) 1137 & - (zca*(xba*ycb+zt)-xca*(zba*ycb-xt))/rcbt2 1138 & - 2.0d0*(zt*xba-zba*xt)*dphidyibt/rt2 1139 & - (xdc*(zdb*ycb+xu)-zdc*(xdb*ycb-zu))/rcbu2 1140 & + 2.0d0*(zu*xdb-zdb*xu)*dphidyibu/ru2 1141 dyibxid = rcbxu*dphidyibu + dphidzu 1142 & + (xdc*zcb*xcb+zdc*yzcb2)/rcbu2 1143 dyibyid = rcbyu*dphidyibu + ycb*yu/rcbu2 1144 dyibzid = rcbzu*dphidyibu - dphidxu 1145 & - (zdc*xcb*zcb+xdc*xycb2)/rcbu2 1146 dzibxic = -xcb*dphidzib/(rcb*rcb) + dphidyt + dphidyu 1147 & - (xca*(yba*xcb-zt)+yca*(zba*zcb+ycb*yba))/rcbt2 1148 & - 2.0d0*(yt*zba-yba*zt)*dphidzibt/rt2 1149 & + (ydc*(zdb*zcb+ycb*ydb)+xdc*(ydb*xcb+zu))/rcbu2 1150 & + 2.0d0*(yu*zdb-ydb*zu)*dphidzibu/ru2 1151 dzibzic = -zcb*dphidzib/(rcb*rcb) 1152 & - (xca*(yba*zcb+xt)-yca*(xba*zcb-yt))/rcbt2 1153 & - 2.0d0*(xt*yba-xba*yt)*dphidzibt/rt2 1154 & - (ydc*(xdb*zcb+yu)-xdc*(ydb*zcb-xu))/rcbu2 1155 & + 2.0d0*(xu*ydb-xdb*yu)*dphidzibu/ru2 1156 dzibxid = rcbxu*dphidzibu - dphidyu 1157 & - (xdc*ycb*xcb+ydc*yzcb2)/rcbu2 1158 dzibyid = rcbyu*dphidzibu + dphidxu 1159 & + (ydc*xcb*ycb+xdc*xzcb2)/rcbu2 1160 dzibzid = rcbzu*dphidzibu + zcb*zu/rcbu2 1161 dxicxid = rcbxu*dphidxicu - xcb*(zdb*ycb-ydb*zcb)/rcbu2 1162 dxicyid = rcbyu*dphidxicu + dphidzu 1163 & + (ydb*zcb*ycb+zdb*xzcb2)/rcbu2 1164 dxiczid = rcbzu*dphidxicu - dphidyu 1165 & - (zdb*ycb*zcb+ydb*xycb2)/rcbu2 1166 dyicxid = rcbxu*dphidyicu - dphidzu 1167 & - (xdb*zcb*xcb+zdb*yzcb2)/rcbu2 1168 dyicyid = rcbyu*dphidyicu - ycb*(xdb*zcb-zdb*xcb)/rcbu2 1169 dyiczid = rcbzu*dphidyicu + dphidxu 1170 & + (zdb*xcb*zcb+xdb*xycb2)/rcbu2 1171 dzicxid = rcbxu*dphidzicu + dphidyu 1172 & + (xdb*ycb*xcb+ydb*yzcb2)/rcbu2 1173 dzicyid = rcbyu*dphidzicu - dphidxu 1174 & - (ydb*xcb*ycb+xdb*xzcb2)/rcbu2 1175 dziczid = rcbzu*dphidzicu - zcb*(ydb*xcb-xdb*ycb)/rcbu2 1176 dxidxid = rcbxu*dphidxid 1177 dxidyid = rcbxu*dphidyid + zcb*rcb/ru2 1178 dxidzid = rcbxu*dphidzid - ycb*rcb/ru2 1179 dyidyid = rcbyu*dphidyid 1180 dyidzid = rcbyu*dphidzid + xcb*rcb/ru2 1181 dzidzid = rcbzu*dphidzid 1182c 1183c get some second derivative chain rule terms by difference 1184c 1185 dxiaxib = -dxiaxia - dxiaxic - dxiaxid 1186 dxiayib = -dxiayia - dxiayic - dxiayid 1187 dxiazib = -dxiazia - dxiazic - dxiazid 1188 dyiayib = -dyiayia - dyiayic - dyiayid 1189 dyiazib = -dyiazia - dyiazic - dyiazid 1190 dziazib = -dziazia - dziazic - dziazid 1191 dxibxib = -dxiaxib - dxibxic - dxibxid 1192 dxibyib = -dyiaxib - dxibyic - dxibyid 1193 dxibzib = -dxiazib - dzibxic - dzibxid 1194 dxibzic = -dziaxib - dxibzib - dxibzid 1195 dyibyib = -dyiayib - dyibyic - dyibyid 1196 dyibzic = -dziayib - dyibzib - dyibzid 1197 dzibzib = -dziazib - dzibzic - dzibzid 1198 dzibyic = -dyiazib - dyibzib - dzibyid 1199 dxicxic = -dxiaxic - dxibxic - dxicxid 1200 dxicyic = -dyiaxic - dyibxic - dxicyid 1201 dxiczic = -dziaxic - dzibxic - dxiczid 1202 dyicyic = -dyiayic - dyibyic - dyicyid 1203 dyiczic = -dziayic - dzibyic - dyiczid 1204 dziczic = -dziazic - dzibzic - dziczid 1205c 1206c increment diagonal and off-diagonal Hessian elements 1207c 1208 if (i .eq. ia) then 1209 hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxia 1210 & + d2edphi2*dphidxia*dphidxia 1211 hessy(1,ia) = hessy(1,ia) + dedphi*dxiayia 1212 & + d2edphi2*dphidxia*dphidyia 1213 hessz(1,ia) = hessz(1,ia) + dedphi*dxiazia 1214 & + d2edphi2*dphidxia*dphidzia 1215 hessx(2,ia) = hessx(2,ia) + dedphi*dxiayia 1216 & + d2edphi2*dphidxia*dphidyia 1217 hessy(2,ia) = hessy(2,ia) + dedphi*dyiayia 1218 & + d2edphi2*dphidyia*dphidyia 1219 hessz(2,ia) = hessz(2,ia) + dedphi*dyiazia 1220 & + d2edphi2*dphidyia*dphidzia 1221 hessx(3,ia) = hessx(3,ia) + dedphi*dxiazia 1222 & + d2edphi2*dphidxia*dphidzia 1223 hessy(3,ia) = hessy(3,ia) + dedphi*dyiazia 1224 & + d2edphi2*dphidyia*dphidzia 1225 hessz(3,ia) = hessz(3,ia) + dedphi*dziazia 1226 & + d2edphi2*dphidzia*dphidzia 1227 hessx(1,ib) = hessx(1,ib) + dedphi*dxiaxib 1228 & + d2edphi2*dphidxia*dphidxib 1229 hessy(1,ib) = hessy(1,ib) + dedphi*dyiaxib 1230 & + d2edphi2*dphidyia*dphidxib 1231 hessz(1,ib) = hessz(1,ib) + dedphi*dziaxib 1232 & + d2edphi2*dphidzia*dphidxib 1233 hessx(2,ib) = hessx(2,ib) + dedphi*dxiayib 1234 & + d2edphi2*dphidxia*dphidyib 1235 hessy(2,ib) = hessy(2,ib) + dedphi*dyiayib 1236 & + d2edphi2*dphidyia*dphidyib 1237 hessz(2,ib) = hessz(2,ib) + dedphi*dziayib 1238 & + d2edphi2*dphidzia*dphidyib 1239 hessx(3,ib) = hessx(3,ib) + dedphi*dxiazib 1240 & + d2edphi2*dphidxia*dphidzib 1241 hessy(3,ib) = hessy(3,ib) + dedphi*dyiazib 1242 & + d2edphi2*dphidyia*dphidzib 1243 hessz(3,ib) = hessz(3,ib) + dedphi*dziazib 1244 & + d2edphi2*dphidzia*dphidzib 1245 hessx(1,ic) = hessx(1,ic) + dedphi*dxiaxic 1246 & + d2edphi2*dphidxia*dphidxic 1247 hessy(1,ic) = hessy(1,ic) + dedphi*dyiaxic 1248 & + d2edphi2*dphidyia*dphidxic 1249 hessz(1,ic) = hessz(1,ic) + dedphi*dziaxic 1250 & + d2edphi2*dphidzia*dphidxic 1251 hessx(2,ic) = hessx(2,ic) + dedphi*dxiayic 1252 & + d2edphi2*dphidxia*dphidyic 1253 hessy(2,ic) = hessy(2,ic) + dedphi*dyiayic 1254 & + d2edphi2*dphidyia*dphidyic 1255 hessz(2,ic) = hessz(2,ic) + dedphi*dziayic 1256 & + d2edphi2*dphidzia*dphidyic 1257 hessx(3,ic) = hessx(3,ic) + dedphi*dxiazic 1258 & + d2edphi2*dphidxia*dphidzic 1259 hessy(3,ic) = hessy(3,ic) + dedphi*dyiazic 1260 & + d2edphi2*dphidyia*dphidzic 1261 hessz(3,ic) = hessz(3,ic) + dedphi*dziazic 1262 & + d2edphi2*dphidzia*dphidzic 1263 hessx(1,id) = hessx(1,id) + dedphi*dxiaxid 1264 & + d2edphi2*dphidxia*dphidxid 1265 hessy(1,id) = hessy(1,id) + dedphi*dyiaxid 1266 & + d2edphi2*dphidyia*dphidxid 1267 hessz(1,id) = hessz(1,id) + dedphi*dziaxid 1268 & + d2edphi2*dphidzia*dphidxid 1269 hessx(2,id) = hessx(2,id) + dedphi*dxiayid 1270 & + d2edphi2*dphidxia*dphidyid 1271 hessy(2,id) = hessy(2,id) + dedphi*dyiayid 1272 & + d2edphi2*dphidyia*dphidyid 1273 hessz(2,id) = hessz(2,id) + dedphi*dziayid 1274 & + d2edphi2*dphidzia*dphidyid 1275 hessx(3,id) = hessx(3,id) + dedphi*dxiazid 1276 & + d2edphi2*dphidxia*dphidzid 1277 hessy(3,id) = hessy(3,id) + dedphi*dyiazid 1278 & + d2edphi2*dphidyia*dphidzid 1279 hessz(3,id) = hessz(3,id) + dedphi*dziazid 1280 & + d2edphi2*dphidzia*dphidzid 1281 else if (i .eq. ib) then 1282 hessx(1,ib) = hessx(1,ib) + dedphi*dxibxib 1283 & + d2edphi2*dphidxib*dphidxib 1284 hessy(1,ib) = hessy(1,ib) + dedphi*dxibyib 1285 & + d2edphi2*dphidxib*dphidyib 1286 hessz(1,ib) = hessz(1,ib) + dedphi*dxibzib 1287 & + d2edphi2*dphidxib*dphidzib 1288 hessx(2,ib) = hessx(2,ib) + dedphi*dxibyib 1289 & + d2edphi2*dphidxib*dphidyib 1290 hessy(2,ib) = hessy(2,ib) + dedphi*dyibyib 1291 & + d2edphi2*dphidyib*dphidyib 1292 hessz(2,ib) = hessz(2,ib) + dedphi*dyibzib 1293 & + d2edphi2*dphidyib*dphidzib 1294 hessx(3,ib) = hessx(3,ib) + dedphi*dxibzib 1295 & + d2edphi2*dphidxib*dphidzib 1296 hessy(3,ib) = hessy(3,ib) + dedphi*dyibzib 1297 & + d2edphi2*dphidyib*dphidzib 1298 hessz(3,ib) = hessz(3,ib) + dedphi*dzibzib 1299 & + d2edphi2*dphidzib*dphidzib 1300 hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxib 1301 & + d2edphi2*dphidxib*dphidxia 1302 hessy(1,ia) = hessy(1,ia) + dedphi*dxiayib 1303 & + d2edphi2*dphidyib*dphidxia 1304 hessz(1,ia) = hessz(1,ia) + dedphi*dxiazib 1305 & + d2edphi2*dphidzib*dphidxia 1306 hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxib 1307 & + d2edphi2*dphidxib*dphidyia 1308 hessy(2,ia) = hessy(2,ia) + dedphi*dyiayib 1309 & + d2edphi2*dphidyib*dphidyia 1310 hessz(2,ia) = hessz(2,ia) + dedphi*dyiazib 1311 & + d2edphi2*dphidzib*dphidyia 1312 hessx(3,ia) = hessx(3,ia) + dedphi*dziaxib 1313 & + d2edphi2*dphidxib*dphidzia 1314 hessy(3,ia) = hessy(3,ia) + dedphi*dziayib 1315 & + d2edphi2*dphidyib*dphidzia 1316 hessz(3,ia) = hessz(3,ia) + dedphi*dziazib 1317 & + d2edphi2*dphidzib*dphidzia 1318 hessx(1,ic) = hessx(1,ic) + dedphi*dxibxic 1319 & + d2edphi2*dphidxib*dphidxic 1320 hessy(1,ic) = hessy(1,ic) + dedphi*dyibxic 1321 & + d2edphi2*dphidyib*dphidxic 1322 hessz(1,ic) = hessz(1,ic) + dedphi*dzibxic 1323 & + d2edphi2*dphidzib*dphidxic 1324 hessx(2,ic) = hessx(2,ic) + dedphi*dxibyic 1325 & + d2edphi2*dphidxib*dphidyic 1326 hessy(2,ic) = hessy(2,ic) + dedphi*dyibyic 1327 & + d2edphi2*dphidyib*dphidyic 1328 hessz(2,ic) = hessz(2,ic) + dedphi*dzibyic 1329 & + d2edphi2*dphidzib*dphidyic 1330 hessx(3,ic) = hessx(3,ic) + dedphi*dxibzic 1331 & + d2edphi2*dphidxib*dphidzic 1332 hessy(3,ic) = hessy(3,ic) + dedphi*dyibzic 1333 & + d2edphi2*dphidyib*dphidzic 1334 hessz(3,ic) = hessz(3,ic) + dedphi*dzibzic 1335 & + d2edphi2*dphidzib*dphidzic 1336 hessx(1,id) = hessx(1,id) + dedphi*dxibxid 1337 & + d2edphi2*dphidxib*dphidxid 1338 hessy(1,id) = hessy(1,id) + dedphi*dyibxid 1339 & + d2edphi2*dphidyib*dphidxid 1340 hessz(1,id) = hessz(1,id) + dedphi*dzibxid 1341 & + d2edphi2*dphidzib*dphidxid 1342 hessx(2,id) = hessx(2,id) + dedphi*dxibyid 1343 & + d2edphi2*dphidxib*dphidyid 1344 hessy(2,id) = hessy(2,id) + dedphi*dyibyid 1345 & + d2edphi2*dphidyib*dphidyid 1346 hessz(2,id) = hessz(2,id) + dedphi*dzibyid 1347 & + d2edphi2*dphidzib*dphidyid 1348 hessx(3,id) = hessx(3,id) + dedphi*dxibzid 1349 & + d2edphi2*dphidxib*dphidzid 1350 hessy(3,id) = hessy(3,id) + dedphi*dyibzid 1351 & + d2edphi2*dphidyib*dphidzid 1352 hessz(3,id) = hessz(3,id) + dedphi*dzibzid 1353 & + d2edphi2*dphidzib*dphidzid 1354 else if (i .eq. ic) then 1355 hessx(1,ic) = hessx(1,ic) + dedphi*dxicxic 1356 & + d2edphi2*dphidxic*dphidxic 1357 hessy(1,ic) = hessy(1,ic) + dedphi*dxicyic 1358 & + d2edphi2*dphidxic*dphidyic 1359 hessz(1,ic) = hessz(1,ic) + dedphi*dxiczic 1360 & + d2edphi2*dphidxic*dphidzic 1361 hessx(2,ic) = hessx(2,ic) + dedphi*dxicyic 1362 & + d2edphi2*dphidxic*dphidyic 1363 hessy(2,ic) = hessy(2,ic) + dedphi*dyicyic 1364 & + d2edphi2*dphidyic*dphidyic 1365 hessz(2,ic) = hessz(2,ic) + dedphi*dyiczic 1366 & + d2edphi2*dphidyic*dphidzic 1367 hessx(3,ic) = hessx(3,ic) + dedphi*dxiczic 1368 & + d2edphi2*dphidxic*dphidzic 1369 hessy(3,ic) = hessy(3,ic) + dedphi*dyiczic 1370 & + d2edphi2*dphidyic*dphidzic 1371 hessz(3,ic) = hessz(3,ic) + dedphi*dziczic 1372 & + d2edphi2*dphidzic*dphidzic 1373 hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxic 1374 & + d2edphi2*dphidxic*dphidxia 1375 hessy(1,ia) = hessy(1,ia) + dedphi*dxiayic 1376 & + d2edphi2*dphidyic*dphidxia 1377 hessz(1,ia) = hessz(1,ia) + dedphi*dxiazic 1378 & + d2edphi2*dphidzic*dphidxia 1379 hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxic 1380 & + d2edphi2*dphidxic*dphidyia 1381 hessy(2,ia) = hessy(2,ia) + dedphi*dyiayic 1382 & + d2edphi2*dphidyic*dphidyia 1383 hessz(2,ia) = hessz(2,ia) + dedphi*dyiazic 1384 & + d2edphi2*dphidzic*dphidyia 1385 hessx(3,ia) = hessx(3,ia) + dedphi*dziaxic 1386 & + d2edphi2*dphidxic*dphidzia 1387 hessy(3,ia) = hessy(3,ia) + dedphi*dziayic 1388 & + d2edphi2*dphidyic*dphidzia 1389 hessz(3,ia) = hessz(3,ia) + dedphi*dziazic 1390 & + d2edphi2*dphidzic*dphidzia 1391 hessx(1,ib) = hessx(1,ib) + dedphi*dxibxic 1392 & + d2edphi2*dphidxic*dphidxib 1393 hessy(1,ib) = hessy(1,ib) + dedphi*dxibyic 1394 & + d2edphi2*dphidyic*dphidxib 1395 hessz(1,ib) = hessz(1,ib) + dedphi*dxibzic 1396 & + d2edphi2*dphidzic*dphidxib 1397 hessx(2,ib) = hessx(2,ib) + dedphi*dyibxic 1398 & + d2edphi2*dphidxic*dphidyib 1399 hessy(2,ib) = hessy(2,ib) + dedphi*dyibyic 1400 & + d2edphi2*dphidyic*dphidyib 1401 hessz(2,ib) = hessz(2,ib) + dedphi*dyibzic 1402 & + d2edphi2*dphidzic*dphidyib 1403 hessx(3,ib) = hessx(3,ib) + dedphi*dzibxic 1404 & + d2edphi2*dphidxic*dphidzib 1405 hessy(3,ib) = hessy(3,ib) + dedphi*dzibyic 1406 & + d2edphi2*dphidyic*dphidzib 1407 hessz(3,ib) = hessz(3,ib) + dedphi*dzibzic 1408 & + d2edphi2*dphidzic*dphidzib 1409 hessx(1,id) = hessx(1,id) + dedphi*dxicxid 1410 & + d2edphi2*dphidxic*dphidxid 1411 hessy(1,id) = hessy(1,id) + dedphi*dyicxid 1412 & + d2edphi2*dphidyic*dphidxid 1413 hessz(1,id) = hessz(1,id) + dedphi*dzicxid 1414 & + d2edphi2*dphidzic*dphidxid 1415 hessx(2,id) = hessx(2,id) + dedphi*dxicyid 1416 & + d2edphi2*dphidxic*dphidyid 1417 hessy(2,id) = hessy(2,id) + dedphi*dyicyid 1418 & + d2edphi2*dphidyic*dphidyid 1419 hessz(2,id) = hessz(2,id) + dedphi*dzicyid 1420 & + d2edphi2*dphidzic*dphidyid 1421 hessx(3,id) = hessx(3,id) + dedphi*dxiczid 1422 & + d2edphi2*dphidxic*dphidzid 1423 hessy(3,id) = hessy(3,id) + dedphi*dyiczid 1424 & + d2edphi2*dphidyic*dphidzid 1425 hessz(3,id) = hessz(3,id) + dedphi*dziczid 1426 & + d2edphi2*dphidzic*dphidzid 1427 else if (i .eq. id) then 1428 hessx(1,id) = hessx(1,id) + dedphi*dxidxid 1429 & + d2edphi2*dphidxid*dphidxid 1430 hessy(1,id) = hessy(1,id) + dedphi*dxidyid 1431 & + d2edphi2*dphidxid*dphidyid 1432 hessz(1,id) = hessz(1,id) + dedphi*dxidzid 1433 & + d2edphi2*dphidxid*dphidzid 1434 hessx(2,id) = hessx(2,id) + dedphi*dxidyid 1435 & + d2edphi2*dphidxid*dphidyid 1436 hessy(2,id) = hessy(2,id) + dedphi*dyidyid 1437 & + d2edphi2*dphidyid*dphidyid 1438 hessz(2,id) = hessz(2,id) + dedphi*dyidzid 1439 & + d2edphi2*dphidyid*dphidzid 1440 hessx(3,id) = hessx(3,id) + dedphi*dxidzid 1441 & + d2edphi2*dphidxid*dphidzid 1442 hessy(3,id) = hessy(3,id) + dedphi*dyidzid 1443 & + d2edphi2*dphidyid*dphidzid 1444 hessz(3,id) = hessz(3,id) + dedphi*dzidzid 1445 & + d2edphi2*dphidzid*dphidzid 1446 hessx(1,ia) = hessx(1,ia) + dedphi*dxiaxid 1447 & + d2edphi2*dphidxid*dphidxia 1448 hessy(1,ia) = hessy(1,ia) + dedphi*dxiayid 1449 & + d2edphi2*dphidyid*dphidxia 1450 hessz(1,ia) = hessz(1,ia) + dedphi*dxiazid 1451 & + d2edphi2*dphidzid*dphidxia 1452 hessx(2,ia) = hessx(2,ia) + dedphi*dyiaxid 1453 & + d2edphi2*dphidxid*dphidyia 1454 hessy(2,ia) = hessy(2,ia) + dedphi*dyiayid 1455 & + d2edphi2*dphidyid*dphidyia 1456 hessz(2,ia) = hessz(2,ia) + dedphi*dyiazid 1457 & + d2edphi2*dphidzid*dphidyia 1458 hessx(3,ia) = hessx(3,ia) + dedphi*dziaxid 1459 & + d2edphi2*dphidxid*dphidzia 1460 hessy(3,ia) = hessy(3,ia) + dedphi*dziayid 1461 & + d2edphi2*dphidyid*dphidzia 1462 hessz(3,ia) = hessz(3,ia) + dedphi*dziazid 1463 & + d2edphi2*dphidzid*dphidzia 1464 hessx(1,ib) = hessx(1,ib) + dedphi*dxibxid 1465 & + d2edphi2*dphidxid*dphidxib 1466 hessy(1,ib) = hessy(1,ib) + dedphi*dxibyid 1467 & + d2edphi2*dphidyid*dphidxib 1468 hessz(1,ib) = hessz(1,ib) + dedphi*dxibzid 1469 & + d2edphi2*dphidzid*dphidxib 1470 hessx(2,ib) = hessx(2,ib) + dedphi*dyibxid 1471 & + d2edphi2*dphidxid*dphidyib 1472 hessy(2,ib) = hessy(2,ib) + dedphi*dyibyid 1473 & + d2edphi2*dphidyid*dphidyib 1474 hessz(2,ib) = hessz(2,ib) + dedphi*dyibzid 1475 & + d2edphi2*dphidzid*dphidyib 1476 hessx(3,ib) = hessx(3,ib) + dedphi*dzibxid 1477 & + d2edphi2*dphidxid*dphidzib 1478 hessy(3,ib) = hessy(3,ib) + dedphi*dzibyid 1479 & + d2edphi2*dphidyid*dphidzib 1480 hessz(3,ib) = hessz(3,ib) + dedphi*dzibzid 1481 & + d2edphi2*dphidzid*dphidzib 1482 hessx(1,ic) = hessx(1,ic) + dedphi*dxicxid 1483 & + d2edphi2*dphidxid*dphidxic 1484 hessy(1,ic) = hessy(1,ic) + dedphi*dxicyid 1485 & + d2edphi2*dphidyid*dphidxic 1486 hessz(1,ic) = hessz(1,ic) + dedphi*dxiczid 1487 & + d2edphi2*dphidzid*dphidxic 1488 hessx(2,ic) = hessx(2,ic) + dedphi*dyicxid 1489 & + d2edphi2*dphidxid*dphidyic 1490 hessy(2,ic) = hessy(2,ic) + dedphi*dyicyid 1491 & + d2edphi2*dphidyid*dphidyic 1492 hessz(2,ic) = hessz(2,ic) + dedphi*dyiczid 1493 & + d2edphi2*dphidzid*dphidyic 1494 hessx(3,ic) = hessx(3,ic) + dedphi*dzicxid 1495 & + d2edphi2*dphidxid*dphidzic 1496 hessy(3,ic) = hessy(3,ic) + dedphi*dzicyid 1497 & + d2edphi2*dphidyid*dphidzic 1498 hessz(3,ic) = hessz(3,ic) + dedphi*dziczid 1499 & + d2edphi2*dphidzid*dphidzic 1500 end if 1501 end if 1502 end do 1503 return 1504 end 1505