1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1993 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################# 9c ## ## 10c ## subroutine egeom1 -- restraint energy & derivatives ## 11c ## ## 12c ############################################################# 13c 14c 15c "egeom1" calculates the energy and first derivatives 16c with respect to Cartesian coordinates due to restraints 17c on positions, distances, angles and torsions as well as 18c Gaussian basin and spherical droplet restraints 19c 20c 21 subroutine egeom1 22 use atomid 23 use atoms 24 use bound 25 use boxes 26 use cell 27 use deriv 28 use energi 29 use group 30 use molcul 31 use math 32 use restrn 33 use usage 34 use virial 35 implicit none 36 integer i,j,k 37 integer ia,ib,ic,id 38 real*8 e,eps,fgrp 39 real*8 de,dt,dt2,deddt 40 real*8 xr,yr,zr 41 real*8 r,r2,r6,r12 42 real*8 dedx,dedy,dedz 43 real*8 angle,target 44 real*8 dot,force 45 real*8 cosine,sine 46 real*8 terma,termc 47 real*8 rab2,rcb2 48 real*8 xia,yia,zia 49 real*8 xib,yib,zib 50 real*8 xic,yic,zic 51 real*8 xid,yid,zid 52 real*8 xab,yab,zab 53 real*8 xba,yba,zba 54 real*8 xcb,ycb,zcb 55 real*8 xdc,ydc,zdc 56 real*8 xca,yca,zca 57 real*8 xdb,ydb,zdb 58 real*8 xad,yad,zad 59 real*8 xbd,ybd,zbd 60 real*8 xcd,ycd,zcd 61 real*8 xp,yp,zp,rp 62 real*8 xt,yt,zt 63 real*8 xu,yu,zu 64 real*8 xtu,ytu,ztu 65 real*8 rt2,ru2,rtru 66 real*8 rcb,dedphi 67 real*8 dedxt,dedyt,dedzt 68 real*8 dedxu,dedyu,dedzu 69 real*8 dedxia,dedyia,dedzia 70 real*8 dedxib,dedyib,dedzib 71 real*8 dedxic,dedyic,dedzic 72 real*8 dedxid,dedyid,dedzid 73 real*8 df1,df2 74 real*8 af1,af2 75 real*8 tf1,tf2,t1,t2 76 real*8 gf1,gf2 77 real*8 weigh,ratio 78 real*8 weigha,weighb 79 real*8 xcm,ycm,zcm 80 real*8 cf1,cf2,vol 81 real*8 c1,c2,c3 82 real*8 xi,yi,zi,ri 83 real*8 a,b,buffer,term 84 real*8 vxx,vyy,vzz 85 real*8 vyx,vzx,vzy 86 real*8 xorig,xorig2 87 real*8 yorig,yorig2 88 real*8 zorig,zorig2 89 logical proceed,intermol 90c 91c 92c zero out the restraint energy term and first derivatives 93c 94 eg = 0.0d0 95 do i = 1, n 96 deg(1,i) = 0.0d0 97 deg(2,i) = 0.0d0 98 deg(3,i) = 0.0d0 99 end do 100c 101c set tolerance for minimum distance and angle values 102c 103 eps = 0.0001d0 104c 105c disable replica mechanism when computing restraint terms 106c 107 if (use_replica) then 108 xorig = xcell 109 yorig = ycell 110 zorig = zcell 111 xorig2 = xcell2 112 yorig2 = ycell2 113 zorig2 = zcell2 114 xcell = xbox 115 ycell = ybox 116 zcell = zbox 117 xcell2 = xbox2 118 ycell2 = ybox2 119 zcell2 = zbox2 120 end if 121c 122c get energy and derivatives for position restraint terms 123c 124 do i = 1, npfix 125 ia = ipfix(i) 126 proceed = .true. 127 if (use_group) call groups (proceed,fgrp,ia,0,0,0,0,0) 128 if (proceed) proceed = (use(ia)) 129 if (proceed) then 130 xr = 0.0d0 131 yr = 0.0d0 132 zr = 0.0d0 133 if (kpfix(1,i) .ne. 0) xr = x(ia) - xpfix(i) 134 if (kpfix(2,i) .ne. 0) yr = y(ia) - ypfix(i) 135 if (kpfix(3,i) .ne. 0) zr = z(ia) - zpfix(i) 136 if (use_bounds) call image (xr,yr,zr) 137 r = sqrt(xr*xr + yr*yr + zr*zr) 138 force = pfix(1,i) 139 dt = max(0.0d0,r-pfix(2,i)) 140 dt2 = dt * dt 141 e = force * dt2 142 de = 2.0d0 * force * dt / max(r,eps) 143c 144c scale the interaction based on its group membership 145c 146 if (use_group) then 147 e = e * fgrp 148 de = de * fgrp 149 end if 150c 151c compute chain rule terms needed for derivatives 152c 153 dedx = de * xr 154 dedy = de * yr 155 dedz = de * zr 156c 157c increment the total energy and first derivatives 158c 159 eg = eg + e 160 deg(1,ia) = deg(1,ia) + dedx 161 deg(2,ia) = deg(2,ia) + dedy 162 deg(3,ia) = deg(3,ia) + dedz 163c 164c increment the internal virial tensor components 165c 166 vxx = xr * dedx 167 vyx = yr * dedx 168 vzx = zr * dedx 169 vyy = yr * dedy 170 vzy = zr * dedy 171 vzz = zr * dedz 172 vir(1,1) = vir(1,1) + vxx 173 vir(2,1) = vir(2,1) + vyx 174 vir(3,1) = vir(3,1) + vzx 175 vir(1,2) = vir(1,2) + vyx 176 vir(2,2) = vir(2,2) + vyy 177 vir(3,2) = vir(3,2) + vzy 178 vir(1,3) = vir(1,3) + vzx 179 vir(2,3) = vir(2,3) + vzy 180 vir(3,3) = vir(3,3) + vzz 181 end if 182 end do 183c 184c get energy and derivatives for distance restraint terms 185c 186 do i = 1, ndfix 187 ia = idfix(1,i) 188 ib = idfix(2,i) 189 proceed = .true. 190 if (use_group) call groups (proceed,fgrp,ia,ib,0,0,0,0) 191 if (proceed) proceed = (use(ia) .or. use(ib)) 192 if (proceed) then 193 xr = x(ia) - x(ib) 194 yr = y(ia) - y(ib) 195 zr = z(ia) - z(ib) 196 intermol = (molcule(ia) .ne. molcule(ib)) 197 if (use_bounds .and. intermol) call image (xr,yr,zr) 198 r = sqrt(xr*xr + yr*yr + zr*zr) 199 force = dfix(1,i) 200 df1 = dfix(2,i) 201 df2 = dfix(3,i) 202 target = r 203 if (r .lt. df1) target = df1 204 if (r .gt. df2) target = df2 205 dt = r - target 206 dt2 = dt * dt 207 e = force * dt2 208 de = 2.0d0 * force * dt / max(r,eps) 209c 210c scale the interaction based on its group membership 211c 212 if (use_group) then 213 e = e * fgrp 214 de = de * fgrp 215 end if 216c 217c compute chain rule terms needed for derivatives 218c 219 dedx = de * xr 220 dedy = de * yr 221 dedz = de * zr 222c 223c increment the total energy and first derivatives 224c 225 eg = eg + e 226 deg(1,ia) = deg(1,ia) + dedx 227 deg(2,ia) = deg(2,ia) + dedy 228 deg(3,ia) = deg(3,ia) + dedz 229 deg(1,ib) = deg(1,ib) - dedx 230 deg(2,ib) = deg(2,ib) - dedy 231 deg(3,ib) = deg(3,ib) - dedz 232c 233c increment the internal virial tensor components 234c 235 vxx = xr * dedx 236 vyx = yr * dedx 237 vzx = zr * dedx 238 vyy = yr * dedy 239 vzy = zr * dedy 240 vzz = zr * dedz 241 vir(1,1) = vir(1,1) + vxx 242 vir(2,1) = vir(2,1) + vyx 243 vir(3,1) = vir(3,1) + vzx 244 vir(1,2) = vir(1,2) + vyx 245 vir(2,2) = vir(2,2) + vyy 246 vir(3,2) = vir(3,2) + vzy 247 vir(1,3) = vir(1,3) + vzx 248 vir(2,3) = vir(2,3) + vzy 249 vir(3,3) = vir(3,3) + vzz 250 end if 251 end do 252c 253c get energy and derivatives for angle restraint terms 254c 255 do i = 1, nafix 256 ia = iafix(1,i) 257 ib = iafix(2,i) 258 ic = iafix(3,i) 259 proceed = .true. 260 if (use_group) call groups (proceed,fgrp,ia,ib,ic,0,0,0) 261 if (proceed) proceed = (use(ia) .or. use(ib) .or. use(ic)) 262 if (proceed) then 263 xia = x(ia) 264 yia = y(ia) 265 zia = z(ia) 266 xib = x(ib) 267 yib = y(ib) 268 zib = z(ib) 269 xic = x(ic) 270 yic = y(ic) 271 zic = z(ic) 272 xab = xia - xib 273 yab = yia - yib 274 zab = zia - zib 275 xcb = xic - xib 276 ycb = yic - yib 277 zcb = zic - zib 278 rab2 = max(xab*xab+yab*yab+zab*zab,eps) 279 rcb2 = max(xcb*xcb+ycb*ycb+zcb*zcb,eps) 280 xp = ycb*zab - zcb*yab 281 yp = zcb*xab - xcb*zab 282 zp = xcb*yab - ycb*xab 283 rp = sqrt(max(xp*xp+yp*yp+zp*zp,eps)) 284 dot = xab*xcb + yab*ycb + zab*zcb 285 cosine = dot / sqrt(rab2*rcb2) 286 cosine = min(1.0d0,max(-1.0d0,cosine)) 287 angle = radian * acos(cosine) 288 force = afix(1,i) 289 af1 = afix(2,i) 290 af2 = afix(3,i) 291 target = angle 292 if (angle .lt. af1) target = af1 293 if (angle .gt. af2) target = af2 294 dt = angle - target 295 dt = dt / radian 296 dt2 = dt * dt 297 e = force * dt2 298 deddt = 2.0d0 * force * dt 299c 300c scale the interaction based on its group membership 301c 302 if (use_group) then 303 e = e * fgrp 304 deddt = deddt * fgrp 305 end if 306c 307c compute derivative components for this interaction 308c 309 terma = -deddt / (rab2*rp) 310 termc = deddt / (rcb2*rp) 311 dedxia = terma * (yab*zp-zab*yp) 312 dedyia = terma * (zab*xp-xab*zp) 313 dedzia = terma * (xab*yp-yab*xp) 314 dedxic = termc * (ycb*zp-zcb*yp) 315 dedyic = termc * (zcb*xp-xcb*zp) 316 dedzic = termc * (xcb*yp-ycb*xp) 317 dedxib = -dedxia - dedxic 318 dedyib = -dedyia - dedyic 319 dedzib = -dedzia - dedzic 320c 321c increment the overall energy term and derivatives 322c 323 eg = eg + e 324 deg(1,ia) = deg(1,ia) + dedxia 325 deg(2,ia) = deg(2,ia) + dedyia 326 deg(3,ia) = deg(3,ia) + dedzia 327 deg(1,ib) = deg(1,ib) + dedxib 328 deg(2,ib) = deg(2,ib) + dedyib 329 deg(3,ib) = deg(3,ib) + dedzib 330 deg(1,ic) = deg(1,ic) + dedxic 331 deg(2,ic) = deg(2,ic) + dedyic 332 deg(3,ic) = deg(3,ic) + dedzic 333c 334c increment the internal virial tensor components 335c 336 vxx = xab*dedxia + xcb*dedxic 337 vyx = yab*dedxia + ycb*dedxic 338 vzx = zab*dedxia + zcb*dedxic 339 vyy = yab*dedyia + ycb*dedyic 340 vzy = zab*dedyia + zcb*dedyic 341 vzz = zab*dedzia + zcb*dedzic 342 vir(1,1) = vir(1,1) + vxx 343 vir(2,1) = vir(2,1) + vyx 344 vir(3,1) = vir(3,1) + vzx 345 vir(1,2) = vir(1,2) + vyx 346 vir(2,2) = vir(2,2) + vyy 347 vir(3,2) = vir(3,2) + vzy 348 vir(1,3) = vir(1,3) + vzx 349 vir(2,3) = vir(2,3) + vzy 350 vir(3,3) = vir(3,3) + vzz 351 end if 352 end do 353c 354c get energy and derivatives for torsion restraint terms 355c 356 do i = 1, ntfix 357 ia = itfix(1,i) 358 ib = itfix(2,i) 359 ic = itfix(3,i) 360 id = itfix(4,i) 361 proceed = .true. 362 if (use_group) call groups (proceed,fgrp,ia,ib,ic,id,0,0) 363 if (proceed) proceed = (use(ia) .or. use(ib) .or. 364 & use(ic) .or. use(id)) 365 if (proceed) then 366 xia = x(ia) 367 yia = y(ia) 368 zia = z(ia) 369 xib = x(ib) 370 yib = y(ib) 371 zib = z(ib) 372 xic = x(ic) 373 yic = y(ic) 374 zic = z(ic) 375 xid = x(id) 376 yid = y(id) 377 zid = z(id) 378 xba = xib - xia 379 yba = yib - yia 380 zba = zib - zia 381 xcb = xic - xib 382 ycb = yic - yib 383 zcb = zic - zib 384 rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps)) 385 xdc = xid - xic 386 ydc = yid - yic 387 zdc = zid - zic 388 xt = yba*zcb - ycb*zba 389 yt = zba*xcb - zcb*xba 390 zt = xba*ycb - xcb*yba 391 xu = ycb*zdc - ydc*zcb 392 yu = zcb*xdc - zdc*xcb 393 zu = xcb*ydc - xdc*ycb 394 xtu = yt*zu - yu*zt 395 ytu = zt*xu - zu*xt 396 ztu = xt*yu - xu*yt 397 rt2 = max(xt*xt+yt*yt+zt*zt,eps) 398 ru2 = max(xu*xu+yu*yu+zu*zu,eps) 399 rtru = sqrt(rt2 * ru2) 400 cosine = (xt*xu + yt*yu + zt*zu) / rtru 401 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru) 402 cosine = min(1.0d0,max(-1.0d0,cosine)) 403 angle = radian * acos(cosine) 404 if (sine .lt. 0.0d0) angle = -angle 405 force = tfix(1,i) 406 tf1 = tfix(2,i) 407 tf2 = tfix(3,i) 408 if (angle.gt.tf1 .and. angle.lt.tf2) then 409 target = angle 410 else if (angle.gt.tf1 .and. tf1.gt.tf2) then 411 target = angle 412 else if (angle.lt.tf2 .and. tf1.gt.tf2) then 413 target = angle 414 else 415 t1 = angle - tf1 416 t2 = angle - tf2 417 if (t1 .gt. 180.0d0) then 418 t1 = t1 - 360.0d0 419 else if (t1 .lt. -180.0d0) then 420 t1 = t1 + 360.0d0 421 end if 422 if (t2 .gt. 180.0d0) then 423 t2 = t2 - 360.0d0 424 else if (t2 .lt. -180.0d0) then 425 t2 = t2 + 360.0d0 426 end if 427 if (abs(t1) .lt. abs(t2)) then 428 target = tf1 429 else 430 target = tf2 431 end if 432 end if 433 dt = angle - target 434 if (dt .gt. 180.0d0) then 435 dt = dt - 360.0d0 436 else if (dt .lt. -180.0d0) then 437 dt = dt + 360.0d0 438 end if 439 dt = dt / radian 440 dt2 = dt * dt 441 e = force * dt2 442 dedphi = 2.0d0 * force * dt 443c 444c scale the interaction based on its group membership 445c 446 if (use_group) then 447 e = e * fgrp 448 dedphi = dedphi * fgrp 449 end if 450c 451c chain rule terms for first derivative components 452c 453 xca = xic - xia 454 yca = yic - yia 455 zca = zic - zia 456 xdb = xid - xib 457 ydb = yid - yib 458 zdb = zid - zib 459 dedxt = dedphi * (yt*zcb - ycb*zt) / (rt2*rcb) 460 dedyt = dedphi * (zt*xcb - zcb*xt) / (rt2*rcb) 461 dedzt = dedphi * (xt*ycb - xcb*yt) / (rt2*rcb) 462 dedxu = -dedphi * (yu*zcb - ycb*zu) / (ru2*rcb) 463 dedyu = -dedphi * (zu*xcb - zcb*xu) / (ru2*rcb) 464 dedzu = -dedphi * (xu*ycb - xcb*yu) / (ru2*rcb) 465c 466c compute derivative components for this interaction 467c 468 dedxia = zcb*dedyt - ycb*dedzt 469 dedyia = xcb*dedzt - zcb*dedxt 470 dedzia = ycb*dedxt - xcb*dedyt 471 dedxib = yca*dedzt - zca*dedyt + zdc*dedyu - ydc*dedzu 472 dedyib = zca*dedxt - xca*dedzt + xdc*dedzu - zdc*dedxu 473 dedzib = xca*dedyt - yca*dedxt + ydc*dedxu - xdc*dedyu 474 dedxic = zba*dedyt - yba*dedzt + ydb*dedzu - zdb*dedyu 475 dedyic = xba*dedzt - zba*dedxt + zdb*dedxu - xdb*dedzu 476 dedzic = yba*dedxt - xba*dedyt + xdb*dedyu - ydb*dedxu 477 dedxid = zcb*dedyu - ycb*dedzu 478 dedyid = xcb*dedzu - zcb*dedxu 479 dedzid = ycb*dedxu - xcb*dedyu 480c 481c increment the overall energy term and derivatives 482c 483 eg = eg + e 484 deg(1,ia) = deg(1,ia) + dedxia 485 deg(2,ia) = deg(2,ia) + dedyia 486 deg(3,ia) = deg(3,ia) + dedzia 487 deg(1,ib) = deg(1,ib) + dedxib 488 deg(2,ib) = deg(2,ib) + dedyib 489 deg(3,ib) = deg(3,ib) + dedzib 490 deg(1,ic) = deg(1,ic) + dedxic 491 deg(2,ic) = deg(2,ic) + dedyic 492 deg(3,ic) = deg(3,ic) + dedzic 493 deg(1,id) = deg(1,id) + dedxid 494 deg(2,id) = deg(2,id) + dedyid 495 deg(3,id) = deg(3,id) + dedzid 496c 497c increment the internal virial tensor components 498c 499 vxx = xcb*(dedxic+dedxid) - xba*dedxia + xdc*dedxid 500 vyx = ycb*(dedxic+dedxid) - yba*dedxia + ydc*dedxid 501 vzx = zcb*(dedxic+dedxid) - zba*dedxia + zdc*dedxid 502 vyy = ycb*(dedyic+dedyid) - yba*dedyia + ydc*dedyid 503 vzy = zcb*(dedyic+dedyid) - zba*dedyia + zdc*dedyid 504 vzz = zcb*(dedzic+dedzid) - zba*dedzia + zdc*dedzid 505 vir(1,1) = vir(1,1) + vxx 506 vir(2,1) = vir(2,1) + vyx 507 vir(3,1) = vir(3,1) + vzx 508 vir(1,2) = vir(1,2) + vyx 509 vir(2,2) = vir(2,2) + vyy 510 vir(3,2) = vir(3,2) + vzy 511 vir(1,3) = vir(1,3) + vzx 512 vir(2,3) = vir(2,3) + vzy 513 vir(3,3) = vir(3,3) + vzz 514 end if 515 end do 516c 517c get energy and derivatives for group distance restraint terms 518c 519 do i = 1, ngfix 520 ia = igfix(1,i) 521 ib = igfix(2,i) 522 xcm = 0.0d0 523 ycm = 0.0d0 524 zcm = 0.0d0 525 do j = igrp(1,ia), igrp(2,ia) 526 k = kgrp(j) 527 weigh = mass(k) 528 xcm = xcm + x(k)*weigh 529 ycm = ycm + y(k)*weigh 530 zcm = zcm + z(k)*weigh 531 end do 532 weigha = max(1.0d0,grpmass(ia)) 533 xr = xcm / weigha 534 yr = ycm / weigha 535 zr = zcm / weigha 536 xcm = 0.0d0 537 ycm = 0.0d0 538 zcm = 0.0d0 539 do j = igrp(1,ib), igrp(2,ib) 540 k = kgrp(j) 541 weigh = mass(k) 542 xcm = xcm + x(k)*weigh 543 ycm = ycm + y(k)*weigh 544 zcm = zcm + z(k)*weigh 545 end do 546 weighb = max(1.0d0,grpmass(ib)) 547 xr = xr - xcm/weighb 548 yr = yr - ycm/weighb 549 zr = zr - zcm/weighb 550 intermol = (molcule(kgrp(igrp(1,ia))) .ne. 551 & molcule(kgrp(igrp(1,ib)))) 552 if (use_bounds .and. intermol) call image (xr,yr,zr) 553 r = sqrt(xr*xr + yr*yr + zr*zr) 554 force = gfix(1,i) 555 gf1 = gfix(2,i) 556 gf2 = gfix(3,i) 557 target = r 558 if (r .lt. gf1) target = gf1 559 if (r .gt. gf2) target = gf2 560 dt = r - target 561 dt2 = dt * dt 562 e = force * dt2 563 de = 2.0d0 * force * dt / max(r,eps) 564c 565c compute chain rule terms needed for derivatives 566c 567 dedx = de * xr 568 dedy = de * yr 569 dedz = de * zr 570c 571c increment the total energy and first derivatives 572c 573 eg = eg + e 574 do j = igrp(1,ia), igrp(2,ia) 575 k = kgrp(j) 576 ratio = mass(k) / weigha 577 deg(1,k) = deg(1,k) + dedx*ratio 578 deg(2,k) = deg(2,k) + dedy*ratio 579 deg(3,k) = deg(3,k) + dedz*ratio 580 end do 581 do j = igrp(1,ib), igrp(2,ib) 582 k = kgrp(j) 583 ratio = mass(k) / weighb 584 deg(1,k) = deg(1,k) - dedx*ratio 585 deg(2,k) = deg(2,k) - dedy*ratio 586 deg(3,k) = deg(3,k) - dedz*ratio 587 end do 588c 589c increment the internal virial tensor components 590c 591 vxx = xr * dedx 592 vyx = yr * dedx 593 vzx = zr * dedx 594 vyy = yr * dedy 595 vzy = zr * dedy 596 vzz = zr * dedz 597 vir(1,1) = vir(1,1) + vxx 598 vir(2,1) = vir(2,1) + vyx 599 vir(3,1) = vir(3,1) + vzx 600 vir(1,2) = vir(1,2) + vyx 601 vir(2,2) = vir(2,2) + vyy 602 vir(3,2) = vir(3,2) + vzy 603 vir(1,3) = vir(1,3) + vzx 604 vir(2,3) = vir(2,3) + vzy 605 vir(3,3) = vir(3,3) + vzz 606 end do 607c 608c get energy and derivatives for chirality restraint terms 609c 610 do i = 1, nchir 611 ia = ichir(1,i) 612 ib = ichir(2,i) 613 ic = ichir(3,i) 614 id = ichir(4,i) 615 proceed = .true. 616 if (use_group) call groups (proceed,fgrp,ia,ib,ic,id,0,0) 617 if (proceed) proceed = (use(ia) .or. use(ib) .or. 618 & use(ic) .or. use(id)) 619 if (proceed) then 620 xad = x(ia) - x(id) 621 yad = y(ia) - y(id) 622 zad = z(ia) - z(id) 623 xbd = x(ib) - x(id) 624 ybd = y(ib) - y(id) 625 zbd = z(ib) - z(id) 626 xcd = x(ic) - x(id) 627 ycd = y(ic) - y(id) 628 zcd = z(ic) - z(id) 629 c1 = ybd*zcd - zbd*ycd 630 c2 = ycd*zad - zcd*yad 631 c3 = yad*zbd - zad*ybd 632 vol = xad*c1 + xbd*c2 + xcd*c3 633 force = chir(1,i) 634 cf1 = chir(2,i) 635 cf2 = chir(3,i) 636 target = vol 637 if (vol .lt. min(cf1,cf2)) target = min(cf1,cf2) 638 if (vol .gt. max(cf1,cf2)) target = max(cf1,cf2) 639 dt = vol - target 640 dt2 = dt * dt 641 e = force * dt2 642 deddt = 2.0d0 * force * dt 643c 644c scale the interaction based on its group membership 645c 646 if (use_group) then 647 e = e * fgrp 648 deddt = deddt * fgrp 649 end if 650c 651c compute derivative components for this interaction 652c 653 dedxia = deddt * (ybd*zcd - zbd*ycd) 654 dedyia = deddt * (zbd*xcd - xbd*zcd) 655 dedzia = deddt * (xbd*ycd - ybd*xcd) 656 dedxib = deddt * (zad*ycd - yad*zcd) 657 dedyib = deddt * (xad*zcd - zad*xcd) 658 dedzib = deddt * (yad*xcd - xad*ycd) 659 dedxic = deddt * (yad*zbd - zad*ybd) 660 dedyic = deddt * (zad*xbd - xad*zbd) 661 dedzic = deddt * (xad*ybd - yad*xbd) 662 dedxid = -dedxia - dedxib - dedxic 663 dedyid = -dedyia - dedyib - dedyic 664 dedzid = -dedzia - dedzib - dedzic 665c 666c increment the overall energy term and derivatives 667c 668 eg = eg + e 669 deg(1,ia) = deg(1,ia) + dedxia 670 deg(2,ia) = deg(2,ia) + dedyia 671 deg(3,ia) = deg(3,ia) + dedzia 672 deg(1,ib) = deg(1,ib) + dedxib 673 deg(2,ib) = deg(2,ib) + dedyib 674 deg(3,ib) = deg(3,ib) + dedzib 675 deg(1,ic) = deg(1,ic) + dedxic 676 deg(2,ic) = deg(2,ic) + dedyic 677 deg(3,ic) = deg(3,ic) + dedzic 678 deg(1,id) = deg(1,id) + dedxid 679 deg(2,id) = deg(2,id) + dedyid 680 deg(3,id) = deg(3,id) + dedzid 681c 682c increment the internal virial tensor components 683c 684 vxx = xad*dedxia + xbd*dedxib + xcd*dedxic 685 vyx = yad*dedxia + ybd*dedxib + ycd*dedxic 686 vzx = zad*dedxia + zbd*dedxib + zcd*dedxic 687 vyy = yad*dedyia + ybd*dedyib + ycd*dedyic 688 vzy = zad*dedyia + zbd*dedyib + zcd*dedyic 689 vzz = zad*dedzia + zbd*dedzib + zcd*dedzic 690 vir(1,1) = vir(1,1) + vxx 691 vir(2,1) = vir(2,1) + vyx 692 vir(3,1) = vir(3,1) + vzx 693 vir(1,2) = vir(1,2) + vyx 694 vir(2,2) = vir(2,2) + vyy 695 vir(3,2) = vir(3,2) + vzy 696 vir(1,3) = vir(1,3) + vzx 697 vir(2,3) = vir(2,3) + vzy 698 vir(3,3) = vir(3,3) + vzz 699 end if 700 end do 701c 702c get energy and derivatives for a Gaussian basin restraint 703c 704 if (use_basin) then 705 do i = 1, n-1 706 xi = x(i) 707 yi = y(i) 708 zi = z(i) 709 do k = i+1, n 710 proceed = .true. 711 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 712 if (proceed) proceed = (use(i) .or. use(k)) 713 if (proceed) then 714 xr = xi - x(k) 715 yr = yi - y(k) 716 zr = zi - z(k) 717 r2 = xr*xr + yr*yr + zr*zr 718 term = -width * r2 719 e = 0.0d0 720 if (term .gt. -50.0d0) e = depth * exp(term) 721 de = -2.0d0 * width * e 722 e = e - depth 723c 724c scale the interaction based on its group membership 725c 726 if (use_group) then 727 e = e * fgrp 728 de = de * fgrp 729 end if 730c 731c compute chain rule terms needed for derivatives 732c 733 dedx = de * xr 734 dedy = de * yr 735 dedz = de * zr 736c 737c increment the overall energy term and derivatives 738c 739 eg = eg + e 740 deg(1,i) = deg(1,i) + dedx 741 deg(2,i) = deg(2,i) + dedy 742 deg(3,i) = deg(3,i) + dedz 743 deg(1,k) = deg(1,k) - dedx 744 deg(2,k) = deg(2,k) - dedy 745 deg(3,k) = deg(3,k) - dedz 746c 747c increment the internal virial tensor components 748c 749 vxx = xr * dedx 750 vyx = yr * dedx 751 vzx = zr * dedx 752 vyy = yr * dedy 753 vzy = zr * dedy 754 vzz = zr * dedz 755 vir(1,1) = vir(1,1) + vxx 756 vir(2,1) = vir(2,1) + vyx 757 vir(3,1) = vir(3,1) + vzx 758 vir(1,2) = vir(1,2) + vyx 759 vir(2,2) = vir(2,2) + vyy 760 vir(3,2) = vir(3,2) + vzy 761 vir(1,3) = vir(1,3) + vzx 762 vir(2,3) = vir(2,3) + vzy 763 vir(3,3) = vir(3,3) + vzz 764 end if 765 end do 766 end do 767 end if 768c 769c get energy and derivatives for a spherical droplet restraint 770c 771 if (use_wall) then 772 buffer = 2.5d0 773 a = 2048.0d0 774 b = 64.0d0 775 do i = 1, n 776 proceed = .true. 777 if (use_group) call groups (proceed,fgrp,i,0,0,0,0,0) 778 if (proceed) proceed = (use(i)) 779 if (proceed) then 780 xi = x(i) 781 yi = y(i) 782 zi = z(i) 783 ri = sqrt(xi**2 + yi**2 + zi**2) 784 r = rwall + buffer - ri 785 r2 = r * r 786 r6 = r2 * r2 * r2 787 r12 = r6 * r6 788 e = a/r12 - b/r6 789 if (ri .eq. 0.0d0) ri = 1.0d0 790 de = (12.0d0*a/r12 - 6.0d0*b/r6) / (r*ri) 791c 792c scale the interaction based on its group membership 793c 794 if (use_group) then 795 e = e * fgrp 796 de = de * fgrp 797 end if 798c 799c compute chain rule terms needed for derivatives 800c 801 dedx = de * xi 802 dedy = de * yi 803 dedz = de * zi 804c 805c increment the overall energy term and derivatives 806c 807 eg = eg + e 808 deg(1,i) = deg(1,i) + dedx 809 deg(2,i) = deg(2,i) + dedy 810 deg(3,i) = deg(3,i) + dedz 811c 812c increment the internal virial tensor components 813c 814 xr = r * xi/ri 815 yr = r * yi/ri 816 zr = r * zi/ri 817 vxx = xr * dedx 818 vyx = yr * dedx 819 vzx = zr * dedx 820 vyy = yr * dedy 821 vzy = zr * dedy 822 vzz = zr * dedz 823 vir(1,1) = vir(1,1) + vxx 824 vir(2,1) = vir(2,1) + vyx 825 vir(3,1) = vir(3,1) + vzx 826 vir(1,2) = vir(1,2) + vyx 827 vir(2,2) = vir(2,2) + vyy 828 vir(3,2) = vir(3,2) + vzy 829 vir(1,3) = vir(1,3) + vzx 830 vir(2,3) = vir(2,3) + vzy 831 vir(3,3) = vir(3,3) + vzz 832 end if 833 end do 834 end if 835c 836c reinstate the replica mechanism if it is being used 837c 838 if (use_replica) then 839 xcell = xorig 840 ycell = yorig 841 zcell = zorig 842 xcell2 = xorig2 843 ycell2 = yorig2 844 zcell2 = zorig2 845 end if 846 return 847 end 848