c
c
c     ##########################################################
c     ##  COPYRIGHT (C) 2014 by Chao Lu & Jay William Ponder  ##
c     ##                  All Rights Reserved                 ##
c     ##########################################################
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine estrtor1  --  stretch-torsion energy & derivs  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "estrtor1" calculates the stretch-torsion energy and first
c     derivatives with respect to Cartesian coordinates
c
c
      subroutine estrtor1
      use atoms
      use bndstr
      use bound
      use deriv
      use energi
      use group
      use strtor
      use torpot
      use tors
      use usage
      use virial
      implicit none
      integer i,k,istrtor
      integer ia,ib,ic,id
      real*8 e,eps,dr,fgrp
      real*8 rt2,ru2,rtru
      real*8 rba,rcb,rdc
      real*8 e1,e2,e3
      real*8 xt,yt,zt
      real*8 xu,yu,zu
      real*8 xtu,ytu,ztu
      real*8 v1,v2,v3
      real*8 c1,c2,c3
      real*8 s1,s2,s3
      real*8 sine,cosine
      real*8 sine2,cosine2
      real*8 sine3,cosine3
      real*8 phi1,phi2,phi3
      real*8 dphi1,dphi2,dphi3
      real*8 xia,yia,zia
      real*8 xib,yib,zib
      real*8 xic,yic,zic
      real*8 xid,yid,zid
      real*8 xba,yba,zba
      real*8 xcb,ycb,zcb
      real*8 xdc,ydc,zdc
      real*8 xca,yca,zca
      real*8 xdb,ydb,zdb
      real*8 ddr,dedphi
      real*8 ddrdx,ddrdy,ddrdz
      real*8 dedxt,dedyt,dedzt
      real*8 dedxu,dedyu,dedzu
      real*8 dedxia,dedyia,dedzia
      real*8 dedxib,dedyib,dedzib
      real*8 dedxic,dedyic,dedzic
      real*8 dedxid,dedyid,dedzid
      real*8 vxx,vyy,vzz
      real*8 vyx,vzx,vzy
      logical proceed
c
c
c     zero out the stretch-torsion energy and first derivatives
c
      ebt = 0.0d0
      do i = 1, n
         debt(1,i) = 0.0d0
         debt(2,i) = 0.0d0
         debt(3,i) = 0.0d0
      end do
      if (nstrtor .eq. 0)  return
c
c     set tolerance for minimum distance and angle values
c
      eps = 0.0001d0
c
c     OpenMP directives for the major loop structure
c
!$OMP PARALLEL default(private) shared(nstrtor,ist,itors,kst,bl,
!$OMP& tors1,tors2,tors3,use,x,y,z,storunit,eps,use_group,use_polymer)
!$OMP& shared(ebt,debt,vir)
!$OMP DO reduction(+:ebt,debt,vir)
c
c     calculate the stretch-torsion energy and first derivatives
c
      do istrtor = 1, nstrtor
         i = ist(1,istrtor)
         ia = itors(1,i)
         ib = itors(2,i)
         ic = itors(3,i)
         id = itors(4,i)
c
c     decide whether to compute the current interaction
c
         proceed = .true.
         if (use_group)  call groups (proceed,fgrp,ia,ib,ic,id,0,0)
         if (proceed)  proceed = (use(ia) .or. use(ib) .or.
     &                              use(ic) .or. use(id))
c
c     compute the value of the torsional angle
c
         if (proceed) then
            xia = x(ia)
            yia = y(ia)
            zia = z(ia)
            xib = x(ib)
            yib = y(ib)
            zib = z(ib)
            xic = x(ic)
            yic = y(ic)
            zic = z(ic)
            xid = x(id)
            yid = y(id)
            zid = z(id)
            xba = xib - xia
            yba = yib - yia
            zba = zib - zia
            xcb = xic - xib
            ycb = yic - yib
            zcb = zic - zib
            xdc = xid - xic
            ydc = yid - yic
            zdc = zid - zic
            if (use_polymer) then
               call image (xba,yba,zba)
               call image (xcb,ycb,zcb)
               call image (xdc,ydc,zdc)
            end if
            rba = sqrt(xba*xba + yba*yba + zba*zba)
            rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb)
            rcb = max(rcb,eps)
            rdc = sqrt(xdc*xdc + ydc*ydc + zdc*zdc)
            xt = yba*zcb - ycb*zba
            yt = zba*xcb - zcb*xba
            zt = xba*ycb - xcb*yba
            xu = ycb*zdc - ydc*zcb
            yu = zcb*xdc - zdc*xcb
            zu = xcb*ydc - xdc*ycb
            xtu = yt*zu - yu*zt
            ytu = zt*xu - zu*xt
            ztu = xt*yu - xu*yt
            rt2 = xt*xt + yt*yt + zt*zt
            rt2 = max(rt2,eps)
            ru2 = xu*xu + yu*yu + zu*zu
            ru2 = max(ru2,eps)
            rtru = sqrt(rt2*ru2)
            xca = xic - xia
            yca = yic - yia
            zca = zic - zia
            xdb = xid - xib
            ydb = yid - yib
            zdb = zid - zib
            if (use_polymer) then
               call image (xca,yca,zca)
               call image (xdb,ydb,zdb)
            end if
            cosine = (xt*xu + yt*yu + zt*zu) / rtru
            sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru)
c
c     compute multiple angle trigonometry and phase terms
c
            c1 = tors1(3,i)
            s1 = tors1(4,i)
            c2 = tors2(3,i)
            s2 = tors2(4,i)
            c3 = tors3(3,i)
            s3 = tors3(4,i)
            cosine2 = cosine*cosine - sine*sine
            sine2 = 2.0d0 * cosine * sine
            cosine3 = cosine*cosine2 - sine*sine2
            sine3 = cosine*sine2 + sine*cosine2
            phi1 = 1.0d0 + (cosine*c1 + sine*s1)
            phi2 = 1.0d0 + (cosine2*c2 + sine2*s2)
            phi3 = 1.0d0 + (cosine3*c3 + sine3*s3)
            dphi1 = cosine*s1 - sine*c1
            dphi2 = 2.0d0 * (cosine2*s2 - sine2*c2)
            dphi3 = 3.0d0 * (cosine3*s3 - sine3*c3)
c
c     get the stretch-torsion values for the first bond
c
            v1 = kst(1,istrtor)
            v2 = kst(2,istrtor)
            v3 = kst(3,istrtor)
            k = ist(2,istrtor)
            dr = rba - bl(k)
            e1 = storunit * dr * (v1*phi1 + v2*phi2 + v3*phi3)
            dedphi = storunit * dr * (v1*dphi1 + v2*dphi2 + v3*dphi3)
            rba = max(rba,eps)
            ddr = storunit * (v1*phi1 + v2*phi2 + v3*phi3) / rba
c
c     scale the interaction based on its group membership
c
            if (use_group) then
               e1 = e1 * fgrp
               dedphi = dedphi * fgrp
               ddr = ddr * fgrp
            end if
c
c     compute derivative components for this interaction
c
            ddrdx = xba * ddr
            ddrdy = yba * ddr
            ddrdz = zba * ddr
            dedxt = dedphi * (yt*zcb - ycb*zt) / (rt2*rcb)
            dedyt = dedphi * (zt*xcb - zcb*xt) / (rt2*rcb)
            dedzt = dedphi * (xt*ycb - xcb*yt) / (rt2*rcb)
            dedxu = -dedphi * (yu*zcb - ycb*zu) / (ru2*rcb)
            dedyu = -dedphi * (zu*xcb - zcb*xu) / (ru2*rcb)
            dedzu = -dedphi * (xu*ycb - xcb*yu) / (ru2*rcb)
c
c     determine chain rule components for the first bond
c
            dedxia = zcb*dedyt - ycb*dedzt - ddrdx
            dedyia = xcb*dedzt - zcb*dedxt - ddrdy
            dedzia = ycb*dedxt - xcb*dedyt - ddrdz
            dedxib = yca*dedzt - zca*dedyt + zdc*dedyu
     &                  - ydc*dedzu + ddrdx
            dedyib = zca*dedxt - xca*dedzt + xdc*dedzu
     &                  - zdc*dedxu + ddrdy
            dedzib = xca*dedyt - yca*dedxt + ydc*dedxu
     &                  - xdc*dedyu + ddrdz
            dedxic = zba*dedyt - yba*dedzt + ydb*dedzu
     &                  - zdb*dedyu
            dedyic = xba*dedzt - zba*dedxt + zdb*dedxu
     &                  - xdb*dedzu
            dedzic = yba*dedxt - xba*dedyt + xdb*dedyu
     &                  - ydb*dedxu
            dedxid = zcb*dedyu - ycb*dedzu
            dedyid = xcb*dedzu - zcb*dedxu
            dedzid = ycb*dedxu - xcb*dedyu
c
c     get the stretch-torsion values for the second bond
c
            v1 = kst(4,istrtor)
            v2 = kst(5,istrtor)
            v3 = kst(6,istrtor)
            k = ist(3,istrtor)
            dr = rcb - bl(k)
            e2 = storunit * dr * (v1*phi1 + v2*phi2 + v3*phi3)
            dedphi = storunit * dr * (v1*dphi1 + v2*dphi2 + v3*dphi3)
            ddr = storunit * (v1*phi1 + v2*phi2 + v3*phi3) / rcb
c
c     scale the interaction based on its group membership
c
            if (use_group) then
               e2 = e2 * fgrp
               dedphi = dedphi * fgrp
               ddr = ddr * fgrp
            end if
c
c     compute derivative components for this interaction
c
            ddrdx = xcb * ddr
            ddrdy = ycb * ddr
            ddrdz = zcb * ddr
            dedxt = dedphi * (yt*zcb - ycb*zt) / (rt2*rcb)
            dedyt = dedphi * (zt*xcb - zcb*xt) / (rt2*rcb)
            dedzt = dedphi * (xt*ycb - xcb*yt) / (rt2*rcb)
            dedxu = -dedphi * (yu*zcb - ycb*zu) / (ru2*rcb)
            dedyu = -dedphi * (zu*xcb - zcb*xu) / (ru2*rcb)
            dedzu = -dedphi * (xu*ycb - xcb*yu) / (ru2*rcb)
c
c     increment chain rule components for the second bond
c
            dedxia = dedxia + zcb*dedyt - ycb*dedzt
            dedyia = dedyia + xcb*dedzt - zcb*dedxt
            dedzia = dedzia + ycb*dedxt - xcb*dedyt
            dedxib = dedxib + yca*dedzt - zca*dedyt + zdc*dedyu
     &                  - ydc*dedzu - ddrdx
            dedyib = dedyib + zca*dedxt - xca*dedzt + xdc*dedzu
     &                  - zdc*dedxu - ddrdy
            dedzib = dedzib + xca*dedyt - yca*dedxt + ydc*dedxu
     &                  - xdc*dedyu - ddrdz
            dedxic = dedxic + zba*dedyt - yba*dedzt + ydb*dedzu
     &                  - zdb*dedyu + ddrdx
            dedyic = dedyic + xba*dedzt - zba*dedxt + zdb*dedxu
     &                  - xdb*dedzu + ddrdy
            dedzic = dedzic + yba*dedxt - xba*dedyt + xdb*dedyu
     &                  - ydb*dedxu + ddrdz
            dedxid = dedxid + zcb*dedyu - ycb*dedzu
            dedyid = dedyid + xcb*dedzu - zcb*dedxu
            dedzid = dedzid + ycb*dedxu - xcb*dedyu
c
c     get the stretch-torsion values for the third bond
c
            v1 = kst(7,istrtor)
            v2 = kst(8,istrtor)
            v3 = kst(9,istrtor)
            k = ist(4,istrtor)
            dr = rdc - bl(k)
            e3 = storunit * dr * (v1*phi1 + v2*phi2 + v3*phi3)
            dedphi = storunit * dr * (v1*dphi1 + v2*dphi2 + v3*dphi3)
            rdc = max(rdc,eps)
            ddr = storunit * (v1*phi1 + v2*phi2 + v3*phi3) / rdc
c
c     scale the interaction based on its group membership
c
            if (use_group) then
               e3 = e3 * fgrp
               dedphi = dedphi * fgrp
               ddr = ddr * fgrp
            end if
c
c     compute derivative components for this interaction
c
            ddrdx = xdc * ddr
            ddrdy = ydc * ddr
            ddrdz = zdc * ddr
            dedxt = dedphi * (yt*zcb - ycb*zt) / (rt2*rcb)
            dedyt = dedphi * (zt*xcb - zcb*xt) / (rt2*rcb)
            dedzt = dedphi * (xt*ycb - xcb*yt) / (rt2*rcb)
            dedxu = -dedphi * (yu*zcb - ycb*zu) / (ru2*rcb)
            dedyu = -dedphi * (zu*xcb - zcb*xu) / (ru2*rcb)
            dedzu = -dedphi * (xu*ycb - xcb*yu) / (ru2*rcb)
c
c     increment chain rule components for the third bond
c
            dedxia = dedxia + zcb*dedyt - ycb*dedzt
            dedyia = dedyia + xcb*dedzt - zcb*dedxt
            dedzia = dedzia + ycb*dedxt - xcb*dedyt
            dedxib = dedxib + yca*dedzt - zca*dedyt + zdc*dedyu
     &                  - ydc*dedzu
            dedyib = dedyib + zca*dedxt - xca*dedzt + xdc*dedzu
     &                  - zdc*dedxu
            dedzib = dedzib + xca*dedyt - yca*dedxt + ydc*dedxu
     &                  - xdc*dedyu
            dedxic = dedxic + zba*dedyt - yba*dedzt + ydb*dedzu
     &                  - zdb*dedyu - ddrdx
            dedyic = dedyic + xba*dedzt - zba*dedxt + zdb*dedxu
     &                  - xdb*dedzu - ddrdy
            dedzic = dedzic + yba*dedxt - xba*dedyt + xdb*dedyu
     &                  - ydb*dedxu - ddrdz
            dedxid = dedxid + zcb*dedyu - ycb*dedzu + ddrdx
            dedyid = dedyid + xcb*dedzu - zcb*dedxu + ddrdy
            dedzid = dedzid + ycb*dedxu - xcb*dedyu + ddrdz
c
c     increment the stretch-torsion energy and gradient
c
            e = e1 + e2 + e3
            ebt = ebt + e
            debt(1,ia) = debt(1,ia) + dedxia
            debt(2,ia) = debt(2,ia) + dedyia
            debt(3,ia) = debt(3,ia) + dedzia
            debt(1,ib) = debt(1,ib) + dedxib
            debt(2,ib) = debt(2,ib) + dedyib
            debt(3,ib) = debt(3,ib) + dedzib
            debt(1,ic) = debt(1,ic) + dedxic
            debt(2,ic) = debt(2,ic) + dedyic
            debt(3,ic) = debt(3,ic) + dedzic
            debt(1,id) = debt(1,id) + dedxid
            debt(2,id) = debt(2,id) + dedyid
            debt(3,id) = debt(3,id) + dedzid
c
c     increment the internal virial tensor components
c
            vxx = xcb*(dedxic+dedxid) - xba*dedxia + xdc*dedxid
            vyx = ycb*(dedxic+dedxid) - yba*dedxia + ydc*dedxid
            vzx = zcb*(dedxic+dedxid) - zba*dedxia + zdc*dedxid
            vyy = ycb*(dedyic+dedyid) - yba*dedyia + ydc*dedyid
            vzy = zcb*(dedyic+dedyid) - zba*dedyia + zdc*dedyid
            vzz = zcb*(dedzic+dedzid) - zba*dedzia + zdc*dedzid
            vir(1,1) = vir(1,1) + vxx
            vir(2,1) = vir(2,1) + vyx
            vir(3,1) = vir(3,1) + vzx
            vir(1,2) = vir(1,2) + vyx
            vir(2,2) = vir(2,2) + vyy
            vir(3,2) = vir(3,2) + vzy
            vir(1,3) = vir(1,3) + vzx
            vir(2,3) = vir(2,3) + vzy
            vir(3,3) = vir(3,3) + vzz
         end if
      end do
c
c     OpenMP directives for the major loop structure
c
!$OMP END DO
!$OMP END PARALLEL
      return
      end