1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1993 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ######################################################## 9c ## ## 10c ## subroutine eimprop -- improper dihedral energy ## 11c ## ## 12c ######################################################## 13c 14c 15c "eimprop" calculates the improper dihedral potential energy 16c 17c 18 subroutine eimprop 19 use atoms 20 use bound 21 use energi 22 use group 23 use improp 24 use math 25 use torpot 26 use usage 27 implicit none 28 integer i,ia,ib,ic,id 29 real*8 e,eps,dt,fgrp 30 real*8 ideal,force 31 real*8 cosine,sine 32 real*8 rcb,angle 33 real*8 xt,yt,zt 34 real*8 xu,yu,zu 35 real*8 xtu,ytu,ztu 36 real*8 rt2,ru2,rtru 37 real*8 xia,yia,zia 38 real*8 xib,yib,zib 39 real*8 xic,yic,zic 40 real*8 xid,yid,zid 41 real*8 xba,yba,zba 42 real*8 xcb,ycb,zcb 43 real*8 xdc,ydc,zdc 44 logical proceed 45c 46c 47c zero out improper dihedral energy 48c 49 eid = 0.0d0 50 if (niprop .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(niprop,iiprop,use, 59!$OMP& x,y,z,kprop,vprop,idihunit,eps,use_group,use_polymer) 60!$OMP& shared(eid) 61!$OMP DO reduction(+:eid) schedule(guided) 62c 63c calculate the improper dihedral angle energy term 64c 65 do i = 1, niprop 66 ia = iiprop(1,i) 67 ib = iiprop(2,i) 68 ic = iiprop(3,i) 69 id = iiprop(4,i) 70c 71c decide whether to compute the current interaction 72c 73 proceed = .true. 74 if (use_group) call groups (proceed,fgrp,ia,ib,ic,id,0,0) 75 if (proceed) proceed = (use(ia) .or. use(ib) .or. 76 & use(ic) .or. use(id)) 77c 78c compute the value of the improper dihedral angle 79c 80 if (proceed) then 81 xia = x(ia) 82 yia = y(ia) 83 zia = z(ia) 84 xib = x(ib) 85 yib = y(ib) 86 zib = z(ib) 87 xic = x(ic) 88 yic = y(ic) 89 zic = z(ic) 90 xid = x(id) 91 yid = y(id) 92 zid = z(id) 93 xba = xib - xia 94 yba = yib - yia 95 zba = zib - zia 96 xcb = xic - xib 97 ycb = yic - yib 98 zcb = zic - zib 99 xdc = xid - xic 100 ydc = yid - yic 101 zdc = zid - zic 102 if (use_polymer) then 103 call image (xba,yba,zba) 104 call image (xcb,ycb,zcb) 105 call image (xdc,ydc,zdc) 106 end if 107 rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps)) 108 xt = yba*zcb - ycb*zba 109 yt = zba*xcb - zcb*xba 110 zt = xba*ycb - xcb*yba 111 xu = ycb*zdc - ydc*zcb 112 yu = zcb*xdc - zdc*xcb 113 zu = xcb*ydc - xdc*ycb 114 xtu = yt*zu - yu*zt 115 ytu = zt*xu - zu*xt 116 ztu = xt*yu - xu*yt 117 rt2 = max(xt*xt+yt*yt+zt*zt,eps) 118 ru2 = max(xu*xu+yu*yu+zu*zu,eps) 119 rtru = sqrt(rt2 * ru2) 120 cosine = (xt*xu + yt*yu + zt*zu) / rtru 121 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru) 122 cosine = min(1.0d0,max(-1.0d0,cosine)) 123 angle = radian * acos(cosine) 124 if (sine .lt. 0.0d0) angle = -angle 125c 126c set the improper dihedral parameters for this angle 127c 128 ideal = vprop(i) 129 force = kprop(i) 130 if (abs(angle+ideal) .lt. abs(angle-ideal)) 131 & ideal = -ideal 132 dt = angle - ideal 133 do while (dt .gt. 180.0d0) 134 dt = dt - 360.0d0 135 end do 136 do while (dt .lt. -180.0d0) 137 dt = dt + 360.0d0 138 end do 139c 140c calculate the improper dihedral energy 141c 142 e = idihunit * force * dt**2 143c 144c scale the interaction based on its group membership 145c 146 if (use_group) e = e * fgrp 147c 148c increment the total improper dihedral energy 149c 150 eid = eid + e 151 end if 152 end do 153c 154c OpenMP directives for the major loop structure 155c 156!$OMP END DO 157!$OMP END PARALLEL 158 return 159 end 160