1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################### 9c ## ## 10c ## subroutine eangle3 -- angle bending energy & analysis ## 11c ## ## 12c ############################################################### 13c 14c 15c "eangle3" calculates the angle bending potential energy, also 16c partitions the energy among the atoms; projected in-plane 17c angles at trigonal centers, spceial linear or Fourier angle 18c bending terms are optionally used 19c 20c 21 subroutine eangle3 22 use action 23 use analyz 24 use angbnd 25 use angpot 26 use atomid 27 use atoms 28 use bound 29 use energi 30 use group 31 use inform 32 use iounit 33 use math 34 use usage 35 implicit none 36 integer i,ia,ib,ic,id 37 real*8 e,eps 38 real*8 ideal,force 39 real*8 fold,factor 40 real*8 dot,cosine 41 real*8 angle,fgrp 42 real*8 dt,dt2,dt3,dt4 43 real*8 xia,yia,zia 44 real*8 xib,yib,zib 45 real*8 xic,yic,zic 46 real*8 xid,yid,zid 47 real*8 xab,yab,zab 48 real*8 xcb,ycb,zcb 49 real*8 xad,yad,zad 50 real*8 xbd,ybd,zbd 51 real*8 xcd,ycd,zcd 52 real*8 xip,yip,zip 53 real*8 xap,yap,zap 54 real*8 xcp,ycp,zcp 55 real*8 rab2,rcb2 56 real*8 rap2,rcp2 57 real*8 xt,yt,zt 58 real*8 rt2,delta 59 logical proceed 60 logical header,huge 61 character*9 label 62c 63c 64c zero out the angle bending energy and partitioning terms 65c 66 nea = 0 67 ea = 0.0d0 68 do i = 1, n 69 aea(i) = 0.0d0 70 end do 71 if (nangle .eq. 0) return 72c 73c set tolerance for minimum distance and angle values 74c 75 eps = 0.0001d0 76c 77c print header information if debug output was requested 78c 79 header = .true. 80 if (debug .and. nangle.ne.0) then 81 header = .false. 82 write (iout,10) 83 10 format (/,' Individual Angle Bending Interactions :', 84 & //,' Type',18x,'Atom Names',18x, 85 & 'Ideal',4x,'Actual',6x,'Energy',/) 86 end if 87c 88c OpenMP directives for the major loop structure 89c 90!$OMP PARALLEL default(private) shared(nangle,iang,anat,ak,afld, 91!$OMP& use,x,y,z,cang,qang,pang,sang,angtyp,angunit,eps,use_group, 92!$OMP& use_polymer, 93!$OMP& name,verbose,debug,header,iout) 94!$OMP& shared(ea,nea,aea) 95!$OMP DO reduction(+:ea,nea,aea) schedule(guided) 96c 97c calculate the bond angle bending energy term 98c 99 do i = 1, nangle 100 ia = iang(1,i) 101 ib = iang(2,i) 102 ic = iang(3,i) 103 id = iang(4,i) 104 ideal = anat(i) 105 force = ak(i) 106c 107c decide whether to compute the current interaction 108c 109 proceed = .true. 110 if (angtyp(i) .eq. 'IN-PLANE') then 111 if (use_group) call groups (proceed,fgrp,ia,ib,ic,id,0,0) 112 if (proceed) proceed = (use(ia) .or. use(ib) .or. 113 & use(ic) .or. use(id)) 114 else 115 if (use_group) call groups (proceed,fgrp,ia,ib,ic,0,0,0) 116 if (proceed) proceed = (use(ia) .or. use(ib) .or. use(ic)) 117 end if 118c 119c get the coordinates of the atoms in the angle 120c 121 if (proceed) then 122 xia = x(ia) 123 yia = y(ia) 124 zia = z(ia) 125 xib = x(ib) 126 yib = y(ib) 127 zib = z(ib) 128 xic = x(ic) 129 yic = y(ic) 130 zic = z(ic) 131c 132c compute the bond angle bending energy 133c 134 if (angtyp(i) .ne. 'IN-PLANE') then 135 xab = xia - xib 136 yab = yia - yib 137 zab = zia - zib 138 xcb = xic - xib 139 ycb = yic - yib 140 zcb = zic - zib 141 if (use_polymer) then 142 call image (xab,yab,zab) 143 call image (xcb,ycb,zcb) 144 end if 145 rab2 = max(xab*xab+yab*yab+zab*zab,eps) 146 rcb2 = max(xcb*xcb+ycb*ycb+zcb*zcb,eps) 147 dot = xab*xcb + yab*ycb + zab*zcb 148 cosine = dot / sqrt(rab2*rcb2) 149 cosine = min(1.0d0,max(-1.0d0,cosine)) 150 angle = radian * acos(cosine) 151 if (angtyp(i) .eq. 'HARMONIC') then 152 dt = angle - ideal 153 dt2 = dt * dt 154 dt3 = dt2 * dt 155 dt4 = dt2 * dt2 156 e = angunit * force * dt2 157 & * (1.0d0+cang*dt+qang*dt2+pang*dt3+sang*dt4) 158 else if (angtyp(i) .eq. 'LINEAR') then 159 factor = 2.0d0 * angunit * radian**2 160 e = factor * force * (1.0d0+cosine) 161 else if (angtyp(i) .eq. 'FOURIER') then 162 fold = afld(i) 163 factor = 2.0d0 * angunit * (radian/fold)**2 164 cosine = cos((fold*angle-ideal)/radian) 165 e = factor * force * (1.0d0+cosine) 166 end if 167c 168c scale the interaction based on its group membership 169c 170 if (use_group) e = e * fgrp 171c 172c increment the total bond angle bending energy 173c 174 nea = nea + 1 175 ea = ea + e 176 aea(ib) = aea(ib) + e 177c 178c print a message if the energy of this interaction is large 179c 180 huge = (e .gt. 5.0d0) 181 if (debug .or. (verbose.and.huge)) then 182 if (header) then 183 header = .false. 184 write (iout,20) 185 20 format (/,' Individual Angle Bending', 186 & ' Interactions :', 187 & //,' Type',18x,'Atom Names',18x, 188 & 'Ideal',4x,'Actual',6x,'Energy',/) 189 end if 190 label = 'Angle ' 191 if (angtyp(i) .eq. 'LINEAR') then 192 label = 'Angle-Lin' 193 else if (angtyp(i) .eq. 'FOURIER') then 194 label = 'Angle-Cos' 195 ideal = (ideal+180.0d0) / fold 196 if (angle-ideal .gt. 180.0d0/fold) 197 & ideal = ideal + 360.0d0/fold 198 end if 199 write (iout,30) label,ia,name(ia),ib,name(ib), 200 & ic,name(ic),ideal,angle,e 201 30 format (1x,a9,1x,i7,'-',a3,i7,'-',a3,i7, 202 & '-',a3,2x,2f10.4,f12.4) 203 end if 204c 205c compute the projected in-plane angle bend energy 206c 207 else 208 xid = x(id) 209 yid = y(id) 210 zid = z(id) 211 xad = xia - xid 212 yad = yia - yid 213 zad = zia - zid 214 xbd = xib - xid 215 ybd = yib - yid 216 zbd = zib - zid 217 xcd = xic - xid 218 ycd = yic - yid 219 zcd = zic - zid 220 if (use_polymer) then 221 call image (xad,yad,zad) 222 call image (xbd,ybd,zbd) 223 call image (xcd,ycd,zcd) 224 end if 225 xt = yad*zcd - zad*ycd 226 yt = zad*xcd - xad*zcd 227 zt = xad*ycd - yad*xcd 228 rt2 = xt*xt + yt*yt + zt*zt 229 delta = -(xt*xbd + yt*ybd + zt*zbd) / rt2 230 xip = xib + xt*delta 231 yip = yib + yt*delta 232 zip = zib + zt*delta 233 xap = xia - xip 234 yap = yia - yip 235 zap = zia - zip 236 xcp = xic - xip 237 ycp = yic - yip 238 zcp = zic - zip 239 if (use_polymer) then 240 call image (xap,yap,zap) 241 call image (xcp,ycp,zcp) 242 end if 243 rap2 = max(xap*xap+yap*yap+zap*zap,eps) 244 rcp2 = max(xcp*xcp+ycp*ycp+zcp*zcp,eps) 245 dot = xap*xcp + yap*ycp + zap*zcp 246 cosine = dot / sqrt(rap2*rcp2) 247 cosine = min(1.0d0,max(-1.0d0,cosine)) 248 angle = radian * acos(cosine) 249 dt = angle - ideal 250 dt2 = dt * dt 251 dt3 = dt2 * dt 252 dt4 = dt2 * dt2 253 e = angunit * force * dt2 254 & * (1.0d0+cang*dt+qang*dt2+pang*dt3+sang*dt4) 255c 256c scale the interaction based on its group membership 257c 258 if (use_group) e = e * fgrp 259c 260c increment the total bond angle bending energy 261c 262 nea = nea + 1 263 ea = ea + e 264 aea(ib) = aea(ib) + e 265c 266c print a message if the energy of this interaction is large 267c 268 huge = (e .gt. 5.0d0) 269 if (debug .or. (verbose.and.huge)) then 270 if (header) then 271 header = .false. 272 write (iout,40) 273 40 format (/,' Individual Angle Bending', 274 & ' Interactions :', 275 & //,' Type',18x,'Atom Names',18x, 276 & 'Ideal',4x,'Actual',6x,'Energy',/) 277 end if 278 write (iout,50) ia,name(ia),ib,name(ib),ic, 279 & name(ic),ideal,angle,e 280 50 format (' Angle-IP',2x,i7,'-',a3,i7,'-',a3,i7, 281 & '-',a3,2x,2f10.4,f12.4) 282 end if 283 end if 284 end if 285 end do 286c 287c OpenMP directives for the major loop structure 288c 289!$OMP END DO 290!$OMP END PARALLEL 291 return 292 end 293