1c 2c 3c ############################################################# 4c ## COPYRIGHT (C) 2003 by Pengyu Ren & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ############################################################# 7c 8c ################################################################# 9c ## ## 10c ## subroutine etortor -- torsion-torsion cross term energy ## 11c ## ## 12c ################################################################# 13c 14c 15c "etortor" calculates the torsion-torsion potential energy 16c 17c 18 subroutine etortor 19 use atoms 20 use bitor 21 use bound 22 use energi 23 use group 24 use ktrtor 25 use math 26 use torpot 27 use tortor 28 use usage 29 implicit none 30 integer i,k,itortor 31 integer pos1,pos2 32 integer ia,ib,ic,id,ie 33 integer nlo,nhi,nt 34 integer xlo,ylo 35 real*8 e,fgrp,sign 36 real*8 angle1,angle2 37 real*8 value1,value2 38 real*8 cosine1,cosine2 39 real*8 xt,yt,zt,rt2 40 real*8 xu,yu,zu,ru2 41 real*8 xv,yv,zv,rv2 42 real*8 rtru,rurv 43 real*8 x1l,x1u 44 real*8 y1l,y1u 45 real*8 xia,yia,zia 46 real*8 xib,yib,zib 47 real*8 xic,yic,zic 48 real*8 xid,yid,zid 49 real*8 xie,yie,zie 50 real*8 xba,yba,zba 51 real*8 xdc,ydc,zdc 52 real*8 xcb,ycb,zcb 53 real*8 xed,yed,zed 54 real*8 ftt(4),ft12(4) 55 real*8 ft1(4),ft2(4) 56 logical proceed 57c 58c 59c zero out the torsion-torsion energy 60c 61 ett = 0.0d0 62 if (ntortor .eq. 0) return 63c 64c OpenMP directives for the major loop structure 65c 66!$OMP PARALLEL default(private) shared(ntortor,itt,ibitor, 67!$OMP& use,x,y,z,tnx,ttx,tny,tty,tbf,tbx,tby,tbxy,ttorunit, 68!$OMP& use_group,use_polymer) 69!$OMP& shared(ett) 70!$OMP DO reduction(+:ett) schedule(guided) 71c 72c calculate the torsion-torsion interaction energy term 73c 74 do itortor = 1, ntortor 75 i = itt(1,itortor) 76 k = itt(2,itortor) 77 if (itt(3,itortor) .eq. 1) then 78 ia = ibitor(1,i) 79 ib = ibitor(2,i) 80 ic = ibitor(3,i) 81 id = ibitor(4,i) 82 ie = ibitor(5,i) 83 else 84 ia = ibitor(5,i) 85 ib = ibitor(4,i) 86 ic = ibitor(3,i) 87 id = ibitor(2,i) 88 ie = ibitor(1,i) 89 end if 90c 91c decide whether to compute the current interaction 92c 93 proceed = .true. 94 if (use_group) call groups (proceed,fgrp,ia,ib,ic,id,ie,0) 95 if (proceed) proceed = (use(ia) .or. use(ib) .or. use(ic) 96 & .or. use(id) .or. use(ie)) 97c 98c compute the values of the torsional angles 99c 100 if (proceed) then 101 xia = x(ia) 102 yia = y(ia) 103 zia = z(ia) 104 xib = x(ib) 105 yib = y(ib) 106 zib = z(ib) 107 xic = x(ic) 108 yic = y(ic) 109 zic = z(ic) 110 xid = x(id) 111 yid = y(id) 112 zid = z(id) 113 xie = x(ie) 114 yie = y(ie) 115 zie = z(ie) 116 xba = xib - xia 117 yba = yib - yia 118 zba = zib - zia 119 xcb = xic - xib 120 ycb = yic - yib 121 zcb = zic - zib 122 xdc = xid - xic 123 ydc = yid - yic 124 zdc = zid - zic 125 xed = xie - xid 126 yed = yie - yid 127 zed = zie - zid 128 if (use_polymer) then 129 call image (xba,yba,zba) 130 call image (xcb,ycb,zcb) 131 call image (xdc,ydc,zdc) 132 call image (xed,yed,zed) 133 end if 134 xt = yba*zcb - ycb*zba 135 yt = zba*xcb - zcb*xba 136 zt = xba*ycb - xcb*yba 137 xu = ycb*zdc - ydc*zcb 138 yu = zcb*xdc - zdc*xcb 139 zu = xcb*ydc - xdc*ycb 140 rt2 = xt*xt + yt*yt + zt*zt 141 ru2 = xu*xu + yu*yu + zu*zu 142 rtru = sqrt(rt2 * ru2) 143 xv = ydc*zed - yed*zdc 144 yv = zdc*xed - zed*xdc 145 zv = xdc*yed - xed*ydc 146 rv2 = xv*xv + yv*yv + zv*zv 147 rurv = sqrt(ru2 * rv2) 148 if (rtru.ne.0.0d0 .and. rurv.ne.0.0d0) then 149 cosine1 = (xt*xu + yt*yu + zt*zu) / rtru 150 cosine1 = min(1.0d0,max(-1.0d0,cosine1)) 151 angle1 = radian * acos(cosine1) 152 sign = xba*xu + yba*yu + zba*zu 153 if (sign .lt. 0.0d0) angle1 = -angle1 154 value1 = angle1 155 cosine2 = (xu*xv + yu*yv + zu*zv) / rurv 156 cosine2 = min(1.0d0,max(-1.0d0,cosine2)) 157 angle2 = radian * acos(cosine2) 158 sign = xcb*xv + ycb*yv + zcb*zv 159 if (sign .lt. 0.0d0) angle2 = -angle2 160 value2 = angle2 161c 162c check for inverted chirality at the central atom 163c 164 call chkttor (ib,ic,id,sign,value1,value2) 165c 166c use bicubic interpolation to compute spline values 167c 168 nlo = 1 169 nhi = tnx(k) 170 do while (nhi-nlo .gt. 1) 171 nt = (nhi+nlo) / 2 172 if (ttx(nt,k) .gt. value1) then 173 nhi = nt 174 else 175 nlo = nt 176 end if 177 end do 178 xlo = nlo 179 nlo = 1 180 nhi = tny(k) 181 do while (nhi-nlo .gt. 1) 182 nt = (nhi + nlo)/2 183 if (tty(nt,k) .gt. value2) then 184 nhi = nt 185 else 186 nlo = nt 187 end if 188 end do 189 ylo = nlo 190 x1l = ttx(xlo,k) 191 x1u = ttx(xlo+1,k) 192 y1l = tty(ylo,k) 193 y1u = tty(ylo+1,k) 194 pos2 = ylo*tnx(k) + xlo 195 pos1 = pos2 - tnx(k) 196 ftt(1) = tbf(pos1,k) 197 ftt(2) = tbf(pos1+1,k) 198 ftt(3) = tbf(pos2+1,k) 199 ftt(4) = tbf(pos2,k) 200 ft1(1) = tbx(pos1,k) 201 ft1(2) = tbx(pos1+1,k) 202 ft1(3) = tbx(pos2+1,k) 203 ft1(4) = tbx(pos2,k) 204 ft2(1) = tby(pos1,k) 205 ft2(2) = tby(pos1+1,k) 206 ft2(3) = tby(pos2+1,k) 207 ft2(4) = tby(pos2,k) 208 ft12(1) = tbxy(pos1,k) 209 ft12(2) = tbxy(pos1+1,k) 210 ft12(3) = tbxy(pos2+1,k) 211 ft12(4) = tbxy(pos2,k) 212 call bcuint (ftt,ft1,ft2,ft12,x1l,x1u, 213 & y1l,y1u,value1,value2,e) 214 e = ttorunit * e 215c 216c scale the interaction based on its group membership 217c 218 if (use_group) e = e * fgrp 219c 220c increment the total torsion-torsion energy 221c 222 ett = ett + e 223 end if 224 end if 225 end do 226c 227c OpenMP directives for the major loop structure 228c 229!$OMP END DO 230!$OMP END PARALLEL 231 return 232 end 233c 234c 235c ############################################################### 236c ## ## 237c ## subroutine chkttor -- check torsion-torsion chirality ## 238c ## ## 239c ############################################################### 240c 241c 242c "chkttor" tests the attached atoms at a torsion-torsion central 243c site and changes the sign of the torsion values if the site has 244c opposite chirality to that for the original parameter 245c 246c note that the sign convention used in this version is correct 247c for phi-psi torsion-torsion interactions defined for L-amino 248c acids in a protein force field; the code may need to be altered 249c for other chiral torsion-torsion situations, such as the sugar 250c rings in nucleic acids 251c 252c 253 subroutine chkttor (ib,ic,id,sign,value1,value2) 254 use atomid 255 use atoms 256 use couple 257 implicit none 258 integer i,j,k,m 259 integer ia,ib,ic,id 260 real*8 sign 261 real*8 value1 262 real*8 value2 263 real*8 xac,yac,zac 264 real*8 xbc,ybc,zbc 265 real*8 xdc,ydc,zdc 266 real*8 c1,c2,c3,vol 267c 268c 269c test for chirality at the central torsion-torsion site 270c 271 sign = 1.0d0 272 if (n12(ic) .eq. 4) then 273 j = 0 274 do i = 1, 4 275 m = i12(i,ic) 276 if (m.ne.ib .and. m.ne.id) then 277 if (j .eq. 0) then 278 j = m 279 else 280 k = m 281 end if 282 end if 283 end do 284 ia = 0 285 if (type(j) .gt. type(k)) ia = j 286 if (type(k) .gt. type(j)) ia = k 287 if (atomic(j) .gt. atomic(k)) ia = j 288 if (atomic(k) .gt. atomic(j)) ia = k 289c 290c compute the signed parallelpiped volume at central site 291c 292 if (ia .ne. 0) then 293 xac = x(ia) - x(ic) 294 yac = y(ia) - y(ic) 295 zac = z(ia) - z(ic) 296 xbc = x(ib) - x(ic) 297 ybc = y(ib) - y(ic) 298 zbc = z(ib) - z(ic) 299 xdc = x(id) - x(ic) 300 ydc = y(id) - y(ic) 301 zdc = z(id) - z(ic) 302 c1 = ybc*zdc - zbc*ydc 303 c2 = ydc*zac - zdc*yac 304 c3 = yac*zbc - zac*ybc 305 vol = xac*c1 + xbc*c2 + xdc*c3 306c 307c invert the angle values if chirality has an inverted sign 308c 309 if (vol .lt. 0.0d0) then 310 sign = -1.0d0 311 value1 = -value1 312 value2 = -value2 313 end if 314 end if 315 end if 316 return 317 end 318