1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1993 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################## 9c ## ## 10c ## subroutine estrbnd -- stretch-bend cross term energy ## 11c ## ## 12c ############################################################## 13c 14c 15c "estrbnd" calculates the stretch-bend potential energy 16c 17c 18 subroutine estrbnd 19 use angbnd 20 use angpot 21 use atoms 22 use bndstr 23 use bound 24 use energi 25 use group 26 use math 27 use strbnd 28 use usage 29 implicit none 30 integer i,j,k 31 integer istrbnd 32 integer ia,ib,ic 33 real*8 e,eps,dt 34 real*8 dr1,dr2 35 real*8 fgrp,angle 36 real*8 force1,force2 37 real*8 dot,cosine 38 real*8 xia,yia,zia 39 real*8 xib,yib,zib 40 real*8 xic,yic,zic 41 real*8 xab,yab,zab 42 real*8 xcb,ycb,zcb 43 real*8 rab,rcb 44 logical proceed 45c 46c 47c zero out the stretch-bend cross term energy 48c 49 eba = 0.0d0 50 if (nstrbnd .eq. 0) return 51c 52c set tolerance for minimum distance and angle values 53c 54 eps = 0.0001d0 55c 56c OpenMP directives for the major loop structure 57c 58!$OMP PARALLEL default(private) shared(nstrbnd,isb,iang,sbk, 59!$OMP& anat,bl,bk,use,x,y,z,stbnunit,eps,use_group,use_polymer) 60!$OMP& shared(eba) 61!$OMP DO reduction(+:eba) schedule(guided) 62c 63c calculate the stretch-bend interaction energy term 64c 65 do istrbnd = 1, nstrbnd 66 i = isb(1,istrbnd) 67 ia = iang(1,i) 68 ib = iang(2,i) 69 ic = iang(3,i) 70 force1 = sbk(1,istrbnd) 71 force2 = sbk(2,istrbnd) 72c 73c decide whether to compute the current interaction 74c 75 proceed = .true. 76 if (use_group) call groups (proceed,fgrp,ia,ib,ic,0,0,0) 77 if (proceed) proceed = (use(ia) .or. use(ib) .or. use(ic)) 78c 79c get the coordinates of the atoms in the angle 80c 81 if (proceed) then 82 xia = x(ia) 83 yia = y(ia) 84 zia = z(ia) 85 xib = x(ib) 86 yib = y(ib) 87 zib = z(ib) 88 xic = x(ic) 89 yic = y(ic) 90 zic = z(ic) 91c 92c compute the value of the bond angle 93c 94 xab = xia - xib 95 yab = yia - yib 96 zab = zia - zib 97 xcb = xic - xib 98 ycb = yic - yib 99 zcb = zic - zib 100 if (use_polymer) then 101 call image (xab,yab,zab) 102 call image (xcb,ycb,zcb) 103 end if 104 rab = sqrt(max(xab*xab+yab*yab+zab*zab,eps)) 105 rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps)) 106 dot = xab*xcb + yab*ycb + zab*zcb 107 cosine = dot / (rab*rcb) 108 cosine = min(1.0d0,max(-1.0d0,cosine)) 109 angle = radian * acos(cosine) 110 dt = angle - anat(i) 111c 112c get the stretch-bend interaction energy 113c 114 j = isb(2,istrbnd) 115 k = isb(3,istrbnd) 116 dr1 = rab - bl(j) 117 dr2 = rcb - bl(k) 118 e = stbnunit * (force1*dr1+force2*dr2) * dt 119c 120c scale the interaction based on its group membership 121c 122 if (use_group) e = e * fgrp 123c 124c increment the total stretch-bend energy 125c 126 eba = eba + e 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