1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1993 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ################################################################ 9c ## ## 10c ## subroutine eurey2 -- atom-by-atom Urey-Bradley Hessian ## 11c ## ## 12c ################################################################ 13c 14c 15c "eurey2" calculates second derivatives of the Urey-Bradley 16c interaction energy for a single atom at a time 17c 18c 19 subroutine eurey2 (i) 20 use atoms 21 use bound 22 use couple 23 use group 24 use hessn 25 use urey 26 use urypot 27 implicit none 28 integer i,j,ia,ic,iurey 29 real*8 ideal,force,fgrp 30 real*8 xac,yac,zac 31 real*8 rac,rac2 32 real*8 dt,dt2,term 33 real*8 termx,termy,termz 34 real*8 de,deddt,d2eddt2 35 real*8 d2e(3,3) 36 logical proceed 37c 38c 39c calculate the Urey-Bradley interaction Hessian elements 40c 41 do iurey = 1, nurey 42 ia = iury(1,iurey) 43 ic = iury(3,iurey) 44 ideal = ul(iurey) 45 force = uk(iurey) 46c 47c decide whether to compute the current interaction 48c 49 proceed = (i.eq.ia .or. i.eq.ic) 50 if (proceed .and. use_group) 51 & call groups (proceed,fgrp,ia,ic,0,0,0,0) 52c 53c compute the value of the 1-3 distance deviation 54c 55 if (proceed) then 56 if (i .eq. ic) then 57 ic = ia 58 ia = i 59 end if 60 xac = x(ia) - x(ic) 61 yac = y(ia) - y(ic) 62 zac = z(ia) - z(ic) 63 if (use_polymer) call image (xac,yac,zac) 64 rac2 = xac*xac + yac*yac + zac*zac 65 rac = sqrt(rac2) 66 dt = rac - ideal 67 dt2 = dt * dt 68 deddt = 2.0d0 * ureyunit * force * dt 69 & * (1.0d0+1.5d0*cury*dt+2.0d0*qury*dt2) 70 d2eddt2 = 2.0d0 * ureyunit * force 71 & * (1.0d0+3.0d0*cury*dt+6.0d0*qury*dt2) 72c 73c scale the interaction based on its group membership 74c 75 if (use_group) then 76 deddt = deddt * fgrp 77 d2eddt2 = d2eddt2 * fgrp 78 end if 79c 80c set the chain rule terms for the Hessian elements 81c 82 de = deddt / rac 83 term = (d2eddt2-de) / rac2 84 termx = term * xac 85 termy = term * yac 86 termz = term * zac 87 d2e(1,1) = termx*xac + de 88 d2e(1,2) = termx*yac 89 d2e(1,3) = termx*zac 90 d2e(2,1) = d2e(1,2) 91 d2e(2,2) = termy*yac + de 92 d2e(2,3) = termy*zac 93 d2e(3,1) = d2e(1,3) 94 d2e(3,2) = d2e(2,3) 95 d2e(3,3) = termz*zac + de 96c 97c increment diagonal and non-diagonal Hessian elements 98c 99 do j = 1, 3 100 hessx(j,ia) = hessx(j,ia) + d2e(1,j) 101 hessy(j,ia) = hessy(j,ia) + d2e(2,j) 102 hessz(j,ia) = hessz(j,ia) + d2e(3,j) 103 hessx(j,ic) = hessx(j,ic) - d2e(1,j) 104 hessy(j,ic) = hessy(j,ic) - d2e(2,j) 105 hessz(j,ic) = hessz(j,ic) - d2e(3,j) 106 end do 107 end if 108 end do 109 return 110 end 111