1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1993 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ################################################################ 9c ## ## 10c ## subroutine eurey1 -- bond stretch energy & derivatives ## 11c ## ## 12c ################################################################ 13c 14c 15c "eurey1" calculates the Urey-Bradley interaction energy and 16c its first derivatives with respect to Cartesian coordinates 17c 18c 19 subroutine eurey1 20 use atoms 21 use bound 22 use deriv 23 use energi 24 use group 25 use urey 26 use urypot 27 use usage 28 use virial 29 implicit none 30 integer i,ia,ic 31 real*8 e,de 32 real*8 ideal,force 33 real*8 dt,dt2,deddt,fgrp 34 real*8 dedx,dedy,dedz 35 real*8 xac,yac,zac,rac 36 real*8 vxx,vyy,vzz 37 real*8 vyx,vzx,vzy 38 logical proceed 39c 40c 41c zero out the Urey-Bradley energy and first derivatives 42c 43 eub = 0.0d0 44 do i = 1, n 45 deub(1,i) = 0.0d0 46 deub(2,i) = 0.0d0 47 deub(3,i) = 0.0d0 48 end do 49 if (nurey .eq. 0) return 50c 51c OpenMP directives for the major loop structure 52c 53!$OMP PARALLEL default(private) shared(nurey,iury,ul,uk, 54!$OMP& use,x,y,z,cury,qury,ureyunit,use_group,use_polymer) 55!$OMP& shared(eub,deub,vir) 56!$OMP DO reduction(+:eub,deub,vir) schedule(guided) 57c 58c calculate the Urey-Bradley 1-3 energy and first derivatives 59c 60 do i = 1, nurey 61 ia = iury(1,i) 62 ic = iury(3,i) 63 ideal = ul(i) 64 force = uk(i) 65c 66c decide whether to compute the current interaction 67c 68 proceed = .true. 69 if (use_group) call groups (proceed,fgrp,ia,ic,0,0,0,0) 70 if (proceed) proceed = (use(ia) .or. use(ic)) 71c 72c compute the value of the 1-3 distance deviation 73c 74 if (proceed) then 75 xac = x(ia) - x(ic) 76 yac = y(ia) - y(ic) 77 zac = z(ia) - z(ic) 78 if (use_polymer) call image (xac,yac,zac) 79 rac = sqrt(xac*xac + yac*yac + zac*zac) 80 dt = rac - ideal 81 dt2 = dt * dt 82 e = ureyunit * force * dt2 * (1.0d0+cury*dt+qury*dt2) 83 deddt = 2.0d0 * ureyunit * force * dt 84 & * (1.0d0+1.5d0*cury*dt+2.0d0*qury*dt2) 85c 86c scale the interaction based on its group membership 87c 88 if (use_group) then 89 e = e * fgrp 90 deddt = deddt * fgrp 91 end if 92c 93c compute chain rule terms needed for derivatives 94c 95 de = deddt / rac 96 dedx = de * xac 97 dedy = de * yac 98 dedz = de * zac 99c 100c increment the total Urey-Bradley energy and first derivatives 101c 102 eub = eub + e 103 deub(1,ia) = deub(1,ia) + dedx 104 deub(2,ia) = deub(2,ia) + dedy 105 deub(3,ia) = deub(3,ia) + dedz 106 deub(1,ic) = deub(1,ic) - dedx 107 deub(2,ic) = deub(2,ic) - dedy 108 deub(3,ic) = deub(3,ic) - dedz 109c 110c increment the internal virial tensor components 111c 112 vxx = xac * dedx 113 vyx = yac * dedx 114 vzx = zac * dedx 115 vyy = yac * dedy 116 vzy = zac * dedy 117 vzz = zac * dedz 118 vir(1,1) = vir(1,1) + vxx 119 vir(2,1) = vir(2,1) + vyx 120 vir(3,1) = vir(3,1) + vzx 121 vir(1,2) = vir(1,2) + vyx 122 vir(2,2) = vir(2,2) + vyy 123 vir(3,2) = vir(3,2) + vzy 124 vir(1,3) = vir(1,3) + vzx 125 vir(2,3) = vir(2,3) + vzy 126 vir(3,3) = vir(3,3) + vzz 127 end if 128 end do 129c 130c OpenMP directives for the major loop structure 131c 132!$OMP END DO 133!$OMP END PARALLEL 134 return 135 end 136