1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1993 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################### 9c ## ## 10c ## subroutine estrbnd3 -- stretch-bend energy & analysis ## 11c ## ## 12c ############################################################### 13c 14c 15c "estrbnd3" calculates the stretch-bend potential energy; 16c also partitions the energy among the atoms 17c 18c 19 subroutine estrbnd3 20 use action 21 use analyz 22 use angbnd 23 use angpot 24 use atomid 25 use atoms 26 use bndstr 27 use bound 28 use energi 29 use group 30 use inform 31 use iounit 32 use math 33 use strbnd 34 use usage 35 implicit none 36 integer i,j,k 37 integer istrbnd 38 integer ia,ib,ic 39 real*8 e,eps,dt 40 real*8 dr1,dr2 41 real*8 fgrp,angle 42 real*8 force1,force2 43 real*8 dot,cosine 44 real*8 xia,yia,zia 45 real*8 xib,yib,zib 46 real*8 xic,yic,zic 47 real*8 xab,yab,zab 48 real*8 xcb,ycb,zcb 49 real*8 rab,rcb 50 logical proceed 51 logical header,huge 52c 53c 54c zero out the energy component and partitioning terms 55c 56 neba = 0 57 eba = 0.0d0 58 do i = 1, n 59 aeba(i) = 0.0d0 60 end do 61 if (nstrbnd .eq. 0) return 62c 63c set tolerance for minimum distance and angle values 64c 65 eps = 0.0001d0 66c 67c print header information if debug output was requested 68c 69 header = .true. 70 if (debug .and. nstrbnd.ne.0) then 71 header = .false. 72 write (iout,10) 73 10 format (/,' Individual Stretch-Bend Interactions :', 74 & //,' Type',18x,'Atom Names',18x,'dSB 1', 75 & 5x,'dSB 2',6x,'Energy',/) 76 end if 77c 78c OpenMP directives for the major loop structure 79c 80!$OMP PARALLEL default(private) shared(nstrbnd,isb,iang,sbk, 81!$OMP& anat,bl,bk,use,x,y,z,stbnunit,eps,use_group,use_polymer, 82!$OMP& name,verbose,debug,header,iout) 83!$OMP& shared(eba,neba,aeba) 84!$OMP DO reduction(+:eba,neba,aeba) schedule(guided) 85c 86c calculate the stretch-bend interaction energy term 87c 88 do istrbnd = 1, nstrbnd 89 i = isb(1,istrbnd) 90 ia = iang(1,i) 91 ib = iang(2,i) 92 ic = iang(3,i) 93 force1 = sbk(1,istrbnd) 94 force2 = sbk(2,istrbnd) 95c 96c decide whether to compute the current interaction 97c 98 proceed = .true. 99 if (use_group) call groups (proceed,fgrp,ia,ib,ic,0,0,0) 100 if (proceed) proceed = (use(ia) .or. use(ib) .or. use(ic)) 101c 102c get the coordinates of the atoms in the angle 103c 104 if (proceed) then 105 xia = x(ia) 106 yia = y(ia) 107 zia = z(ia) 108 xib = x(ib) 109 yib = y(ib) 110 zib = z(ib) 111 xic = x(ic) 112 yic = y(ic) 113 zic = z(ic) 114c 115c compute the value of the bond angle 116c 117 xab = xia - xib 118 yab = yia - yib 119 zab = zia - zib 120 xcb = xic - xib 121 ycb = yic - yib 122 zcb = zic - zib 123 if (use_polymer) then 124 call image (xab,yab,zab) 125 call image (xcb,ycb,zcb) 126 end if 127 rab = sqrt(max(xab*xab+yab*yab+zab*zab,eps)) 128 rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps)) 129 dot = xab*xcb + yab*ycb + zab*zcb 130 cosine = dot / (rab*rcb) 131 cosine = min(1.0d0,max(-1.0d0,cosine)) 132 angle = radian * acos(cosine) 133 dt = angle - anat(i) 134c 135c get the stretch-bend interaction energy 136c 137 j = isb(2,istrbnd) 138 k = isb(3,istrbnd) 139 dr1 = rab - bl(j) 140 dr2 = rcb - bl(k) 141 e = stbnunit * (force1*dr1+force2*dr2) * dt 142c 143c scale the interaction based on its group membership 144c 145 if (use_group) e = e * fgrp 146c 147c increment the total stretch-bend energy 148c 149 neba = neba + 1 150 eba = eba + e 151 aeba(ib) = aeba(ib) + e 152c 153c print a message if the energy of this interaction is large 154c 155 huge = (abs(e) .gt. 5.0d0) 156 if (debug .or. (verbose.and.huge)) then 157 if (header) then 158 header = .false. 159 write (iout,20) 160 20 format (/,' Individual Stretch-Bend', 161 & ' Interactions :', 162 & //,' Type',18x,'Atom Names',18x,'dSB 1', 163 & 5x,'dSB 2',6x,'Energy',/) 164 end if 165 write (iout,30) ia,name(ia),ib,name(ib), 166 & ic,name(ic),dr1*dt,dr2*dt,e 167 30 format (' StrBend',3x,3(i7,'-',a3),2x,2f10.4,f12.4) 168 end if 169 end if 170 end do 171c 172c OpenMP directives for the major loop structure 173c 174!$OMP END DO 175!$OMP END PARALLEL 176 return 177 end 178