+!--------------------------------------------------------------------------
+ subroutine triple_ssbond_ene(resi,resj,resk,eij)
+! implicit none
+! Includes
+ use calc_data
+ use comm_sschecks
+! include 'DIMENSIONS'
+! include 'COMMON.SBRIDGE'
+! include 'COMMON.CHAIN'
+! include 'COMMON.DERIV'
+! include 'COMMON.LOCAL'
+! include 'COMMON.INTERACT'
+! include 'COMMON.VAR'
+! include 'COMMON.IOUNITS'
+! include 'COMMON.CALC'
+#ifndef CLUST
+#ifndef WHAM
+ use MD_data
+! include 'COMMON.MD'
+! use MD, only: totT,t_bath
+#endif
+#endif
+ double precision h_base
+ external h_base
+
+!c Input arguments
+ integer resi,resj,resk,m,itypi,itypj,itypk
+
+!c Output arguments
+ double precision eij,eij1,eij2,eij3
+
+!c Local variables
+ logical havebond
+!c integer itypi,itypj,k,l
+ double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
+ double precision rrik,rrjk,rik,rjk,xi,xk,yi,yk,zi,zk,xij,yij,zij
+ double precision xik,yik,zik,xjk,yjk,zjk,dxk,dyk,dzk
+ double precision sig0ij,ljd,sig,fac,e1,e2
+ double precision dcosom1(3),dcosom2(3),ed
+ double precision pom1,pom2
+ double precision ljA,ljB,ljXs
+ double precision d_ljB(1:3)
+ double precision ssA,ssB,ssC,ssXs
+ double precision ssxm,ljxm,ssm,ljm
+ double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+ eij=0.0
+ if (dtriss.eq.0) return
+ i=resi
+ j=resj
+ k=resk
+!C write(iout,*) resi,resj,resk
+ itypi=itype(i)
+ dxi=dc_norm(1,nres+i)
+ dyi=dc_norm(2,nres+i)
+ dzi=dc_norm(3,nres+i)
+ dsci_inv=vbld_inv(i+nres)
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ itypj=itype(j)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ dzj=dc_norm(3,nres+j)
+ dscj_inv=vbld_inv(j+nres)
+ itypk=itype(k)
+ xk=c(1,nres+k)
+ yk=c(2,nres+k)
+ zk=c(3,nres+k)
+
+ dxk=dc_norm(1,nres+k)
+ dyk=dc_norm(2,nres+k)
+ dzk=dc_norm(3,nres+k)
+ dscj_inv=vbld_inv(k+nres)
+ xij=xj-xi
+ xik=xk-xi
+ xjk=xk-xj
+ yij=yj-yi
+ yik=yk-yi
+ yjk=yk-yj
+ zij=zj-zi
+ zik=zk-zi
+ zjk=zk-zj
+ rrij=(xij*xij+yij*yij+zij*zij)
+ rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+ rrik=(xik*xik+yik*yik+zik*zik)
+ rik=dsqrt(rrik)
+ rrjk=(xjk*xjk+yjk*yjk+zjk*zjk)
+ rjk=dsqrt(rrjk)
+!C there are three combination of distances for each trisulfide bonds
+!C The first case the ith atom is the center
+!C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
+!C distance y is second distance the a,b,c,d are parameters derived for
+!C this problem d parameter was set as a penalty currenlty set to 1.
+ if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
+ eij1=0.0d0
+ else
+ eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
+ endif
+!C second case jth atom is center
+ if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
+ eij2=0.0d0
+ else
+ eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
+ endif
+!C the third case kth atom is the center
+ if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
+ eij3=0.0d0
+ else
+ eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
+ endif
+!C eij2=0.0
+!C eij3=0.0
+!C eij1=0.0
+ eij=eij1+eij2+eij3
+!C write(iout,*)i,j,k,eij
+!C The energy penalty calculated now time for the gradient part
+!C derivative over rij
+ fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)
+ gg(1)=xij*fac/rij
+ gg(2)=yij*fac/rij
+ gg(3)=zij*fac/rij
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,j)=gvdwx(m,j)+gg(m)
+ enddo
+
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ enddo
+!C now derivative over rik
+ fac=-eij1**2/dtriss* &
+ (-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5) &
+ -eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xik*fac/rik
+ gg(2)=yik*fac/rik
+ gg(3)=zik*fac/rik
+ do m=1,3
+ gvdwx(m,i)=gvdwx(m,i)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+!C now derivative over rjk
+ fac=-eij2**2/dtriss* &
+ (-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)- &
+ eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
+ gg(1)=xjk*fac/rjk
+ gg(2)=yjk*fac/rjk
+ gg(3)=zjk*fac/rjk
+ do m=1,3
+ gvdwx(m,j)=gvdwx(m,j)-gg(m)
+ gvdwx(m,k)=gvdwx(m,k)+gg(m)
+ enddo
+ do l=1,3
+ gvdwc(l,j)=gvdwc(l,j)-gg(l)
+ gvdwc(l,k)=gvdwc(l,k)+gg(l)
+ enddo
+ return
+ end subroutine triple_ssbond_ene
+
+
+