c----------------------------------------------------------------------------
      subroutine check_energies
c      implicit none

c     Includes
      include 'DIMENSIONS'
      include 'COMMON.CHAIN'
      include 'COMMON.VAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.SBRIDGE'
      include 'COMMON.LOCAL'
      include 'COMMON.GEO'

c     External functions
      double precision ran_number
      external ran_number

c     Local variables
      integer i,j,k,l,lmax,p,pmax
      double precision rmin,rmax
      double precision eij

      double precision d
      double precision wi,rij,tj,pj


c      return

      i=5
      j=14

      d=dsc(1)
      rmin=2.0D0
      rmax=12.0D0

      lmax=10000
      pmax=1

      do k=1,3
        c(k,i)=0.0D0
        c(k,j)=0.0D0
        c(k,nres+i)=0.0D0
        c(k,nres+j)=0.0D0
      enddo

      do l=1,lmax

ct        wi=ran_number(0.0D0,pi)
c        wi=ran_number(0.0D0,pi/6.0D0)
c        wi=0.0D0
ct        tj=ran_number(0.0D0,pi)
ct        pj=ran_number(0.0D0,pi)
c        pj=ran_number(0.0D0,pi/6.0D0)
c        pj=0.0D0

        do p=1,pmax
ct           rij=ran_number(rmin,rmax)

           c(1,j)=d*sin(pj)*cos(tj)
           c(2,j)=d*sin(pj)*sin(tj)
           c(3,j)=d*cos(pj)

           c(3,nres+i)=-rij

           c(1,i)=d*sin(wi)
           c(3,i)=-rij-d*cos(wi)

           do k=1,3
              dc(k,nres+i)=c(k,nres+i)-c(k,i)
              dc_norm(k,nres+i)=dc(k,nres+i)/d
              dc(k,nres+j)=c(k,nres+j)-c(k,j)
              dc_norm(k,nres+j)=dc(k,nres+j)/d
           enddo

           call dyn_ssbond_ene(i,j,eij)
        enddo
      enddo

      call exit(1)

      return
      end

C-----------------------------------------------------------------------------
      subroutine dyn_ssbond_ene(resi,resj,eij)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'COMMON.SBRIDGE'
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.CALC'
      include 'COMMON.NAMES'
#ifndef CLUST
#ifndef WHAM
      include 'COMMON.MD'
#endif
#endif

c     External functions
      double precision h_base
      external h_base

c     Input arguments
      integer resi,resj

c     Output arguments
      double precision eij

c     Local variables
      logical havebond
c      integer itypi,itypj,k,l
      double precision rrij,ssd,deltat1,deltat2,deltat12,cosphi
      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)
      double precision f1,f2,h1,h2,hd1,hd2
      double precision omega,delta_inv,deltasq_inv,fac1,fac2
c-------FIRST METHOD
      double precision xm,d_xm(1:3)
      double precision sslipi,sslipj,ssgradlipi,ssgradlipj
      integer ici,icj,itypi,itypj
      double precision boxshift,sscale,sscagrad
      double precision aa,bb
c-------END FIRST METHOD
c-------SECOND METHOD
c$$$      double precision ss,d_ss(0:3),ljf,d_ljf(0:3)
c-------END SECOND METHOD

c-------TESTING CODE
      logical checkstop,transgrad
      common /sschecks/ checkstop,transgrad

      integer icheck,nicheck,jcheck,njcheck
      double precision echeck(-1:1),deps,ssx0,ljx0,xi,yi,zi
c-------END TESTING CODE


      i=resi
      j=resj
      ici=icys(i)
      icj=icys(j)
      if (ici.eq.0 .or. icj.eq.0) then
        write (*,'(a,i5,2a,a3,i5,5h and ,a3,i5)') 
     &  "Attempt to create",
     &  " a disulfide link between non-cysteine residues ",restyp(i),i,
     &  restyp(j),j
        stop
      endif
      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)
      call to_box(xi,yi,zi)
C define scaling factor for lipids

C        if (positi.le.0) positi=positi+boxzsize
C        print *,i
C first for peptide groups
c for each residue check if it is in lipid or lipid water border area
      call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
      itypj=itype(j)
      xj=c(1,nres+j)
      yj=c(2,nres+j)
      zj=c(3,nres+j)
      call to_box(xj,yj,zj)
      call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
      aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
     &  +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
      bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
     &  +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
      xj=boxshift(xj-xi,boxxsize)
      yj=boxshift(yj-yi,boxysize)
      zj=boxshift(zj-zi,boxzsize)
      dxj=dc_norm(1,nres+j)
      dyj=dc_norm(2,nres+j)
      dzj=dc_norm(3,nres+j)
      dscj_inv=vbld_inv(j+nres)

      chi1=chi(itypi,itypj)
      chi2=chi(itypj,itypi)
      chi12=chi1*chi2
      chip1=chip(itypi)
      chip2=chip(itypj)
      chip12=chip1*chip2
      alf1=alp(itypi)
      alf2=alp(itypj)
      alf12=0.5D0*(alf1+alf2)

      rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
      rij=dsqrt(rrij)  ! sc_angular needs rij to really be the inverse
            sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
            sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
c     The following are set in sc_angular
c      erij(1)=xj*rij
c      erij(2)=yj*rij
c      erij(3)=zj*rij
c      om1=dxi*erij(1)+dyi*erij(2)+dzi*erij(3)
c      om2=dxj*erij(1)+dyj*erij(2)+dzj*erij(3)
c      om12=dxi*dxj+dyi*dyj+dzi*dzj
      call sc_angular
      rij=1.0D0/rij  ! Reset this so it makes sense

      sig0ij=sigma(itypi,itypj)
      sig=sig0ij*dsqrt(1.0D0/sigsq)

      ljXs=sig-sig0ij
      ljA=eps1*eps2rt**2*eps3rt**2
      ljB=ljA*bb
      ljA=ljA*aa
      ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)

      ssXs=d0cm
      deltat1=1.0d0-om1
      deltat2=1.0d0+om2
      deltat12=om2-om1+2.0d0
      cosphi=om12-om1*om2
      ssA=akcm
      ssB=akct*deltat12
      ssC=ss_depth
     &     +akth*(deltat1*deltat1+deltat2*deltat2)
     &     +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
      ssxm=ssXs-0.5D0*ssB/ssA

c-------TESTING CODE
c$$$c     Some extra output
c$$$      ssm=ssC-0.25D0*ssB*ssB/ssA
c$$$      ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
c$$$      ssx0=ssB*ssB-4.0d0*ssA*ssC
c$$$      if (ssx0.gt.0.0d0) then
c$$$        ssx0=ssXs+0.5d0*(-ssB+sqrt(ssx0))/ssA
c$$$      else
c$$$        ssx0=ssxm
c$$$      endif
c$$$      ljx0=ljXs+(-aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
c$$$      write(iout,'(a,4f8.2,2f15.2,3f6.2)')"SSENERGIES ",
c$$$     &     ssxm,ljxm,ssx0,ljx0,ssm,ljm,om1,om2,om12
c$$$      return
c-------END TESTING CODE

c-------TESTING CODE
c     Stop and plot energy and derivative as a function of distance
      if (checkstop) then
        ssm=ssC-0.25D0*ssB*ssB/ssA
        ljm=-0.25D0*ljB*bb/aa
        if (ssm.lt.ljm .and.
     &       dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
          nicheck=1000
          njcheck=1
          deps=0.5d-7
        else
          checkstop=.false.
        endif
      endif
      if (.not.checkstop) then
        nicheck=0
        njcheck=-1
      endif

      do icheck=0,nicheck
      do jcheck=-1,njcheck
      if (checkstop) rij=(ssxm-1.0d0)+
     &       ((ljxm-ssxm+2.0d0)*icheck)/nicheck+jcheck*deps
c-------END TESTING CODE

      if (rij.gt.ljxm) then
        havebond=.false.
        ljd=rij-ljXs
        fac=(1.0D0/ljd)**expon
        e1=fac*fac*aa
        e2=fac*bb
        eij=eps1*eps2rt*eps3rt*(e1+e2)
        eps2der=eij*eps3rt
        eps3der=eij*eps2rt
        eij=eij*eps2rt*eps3rt*sss

        sigder=-sig/sigsq
        e1=e1*eps1*eps2rt**2*eps3rt**2
        ed=-expon*(e1+eij)/ljd
        sigder=ed*sigder
        ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
        eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
        eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
        eom12=eij*eps1_om12+eps2der*eps2rt_om12
     &       -2.0D0*alf12*eps3der+sigder*sigsq_om12
      else if (rij.lt.ssxm) then
        havebond=.true.
        ssd=rij-ssXs
        eij=ssA*ssd*ssd+ssB*ssd+ssC
        eij=eij*sss        
        ed=2*akcm*ssd+akct*deltat12
        ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
        pom1=akct*ssd
        pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
        eom1=-2*akth*deltat1-pom1-om2*pom2
        eom2= 2*akth*deltat2+pom1-om1*pom2
        eom12=pom2
      else
        omega=v1ss+2.0d0*v2ss*cosphi+3.0d0*v3ss*cosphi*cosphi

        d_ssxm(1)=0.5D0*akct/ssA
        d_ssxm(2)=-d_ssxm(1)
        d_ssxm(3)=0.0D0

        d_ljxm(1)=sig0ij/sqrt(sigsq**3)
        d_ljxm(2)=d_ljxm(1)*sigsq_om2
        d_ljxm(3)=d_ljxm(1)*sigsq_om12
        d_ljxm(1)=d_ljxm(1)*sigsq_om1

c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
        xm=0.5d0*(ssxm+ljxm)
        do k=1,3
          d_xm(k)=0.5d0*(d_ssxm(k)+d_ljxm(k))
        enddo
        if (rij.lt.xm) then
          havebond=.true.
          ssm=ssC-0.25D0*ssB*ssB/ssA
          d_ssm(1)=0.5D0*akct*ssB/ssA
          d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
          d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
          d_ssm(3)=omega
          f1=(rij-xm)/(ssxm-xm)
          f2=(rij-ssxm)/(xm-ssxm)
          h1=h_base(f1,hd1)
          h2=h_base(f2,hd2)
          eij=ssm*h1+Ht*h2
          delta_inv=1.0d0/(xm-ssxm)
          deltasq_inv=delta_inv*delta_inv
          fac=ssm*hd1-Ht*hd2
          fac1=deltasq_inv*fac*(xm-rij)
          fac2=deltasq_inv*fac*(rij-ssxm)
          ed=delta_inv*(Ht*hd2-ssm*hd1)
          eij=eij*sss
          ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
          eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
          eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
          eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
        else
          havebond=.false.
          ljm=-0.25D0*ljB*bb/aa
          d_ljm(1)=-0.5D0*bb/aa*ljB
          d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
          d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
     +         alf12/eps3rt)
          d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt)
          f1=(rij-ljxm)/(xm-ljxm)
          f2=(rij-xm)/(ljxm-xm)
          h1=h_base(f1,hd1)
          h2=h_base(f2,hd2)
          eij=Ht*h1+ljm*h2
          delta_inv=1.0d0/(ljxm-xm)
          deltasq_inv=delta_inv*delta_inv
          fac=Ht*hd1-ljm*hd2
          fac1=deltasq_inv*fac*(ljxm-rij)
          fac2=deltasq_inv*fac*(rij-xm)
          ed=delta_inv*(ljm*hd2-Ht*hd1)
          eij=eij*sss
          ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
          eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
          eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
          eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
        endif
c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE

c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
c$$$        ssd=rij-ssXs
c$$$        ljd=rij-ljXs
c$$$        fac1=rij-ljxm
c$$$        fac2=rij-ssxm
c$$$
c$$$        d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt)
c$$$        d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt)
c$$$        d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt)
c$$$
c$$$        ssm=ssC-0.25D0*ssB*ssB/ssA
c$$$        d_ssm(1)=0.5D0*akct*ssB/ssA
c$$$        d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1)
c$$$        d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1)
c$$$        d_ssm(3)=omega
c$$$
c$$$        ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj)
c$$$        do k=1,3
c$$$          d_ljm(k)=ljm*d_ljB(k)
c$$$        enddo
c$$$        ljm=ljm*ljB
c$$$
c$$$        ss=ssA*ssd*ssd+ssB*ssd+ssC
c$$$        d_ss(0)=2.0d0*ssA*ssd+ssB
c$$$        d_ss(2)=akct*ssd
c$$$        d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega
c$$$        d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega
c$$$        d_ss(3)=omega
c$$$
c$$$        ljf=bb(itypi,itypj)/aa(itypi,itypj)
c$$$        ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0)
c$$$        d_ljf(0)=ljf*2.0d0*ljB*fac1
c$$$        do k=1,3
c$$$          d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1-
c$$$     &         2.0d0*ljB*fac1*d_ljxm(k))
c$$$        enddo
c$$$        ljf=ljm+ljf*ljB*fac1*fac1
c$$$
c$$$        f1=(rij-ljxm)/(ssxm-ljxm)
c$$$        f2=(rij-ssxm)/(ljxm-ssxm)
c$$$        h1=h_base(f1,hd1)
c$$$        h2=h_base(f2,hd2)
c$$$        eij=ss*h1+ljf*h2
c$$$        delta_inv=1.0d0/(ljxm-ssxm)
c$$$        deltasq_inv=delta_inv*delta_inv
c$$$        fac=ljf*hd2-ss*hd1
c$$$        ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac
c$$$        eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac*
c$$$     &       (fac1*d_ssxm(1)-fac2*(d_ljxm(1)))
c$$$        eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac*
c$$$     &       (fac1*d_ssxm(2)-fac2*(d_ljxm(2)))
c$$$        eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac*
c$$$     &       (fac1*d_ssxm(3)-fac2*(d_ljxm(3)))
c$$$
c$$$        havebond=.false.
c$$$        if (ed.gt.0.0d0) havebond=.true.
c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE

      endif

      if (havebond) then
#ifndef CLUST
#ifndef WHAM
c        if (dyn_ssbond_ij(i,j).eq.1.0d300) then
c          write(iout,'(a15,f12.2,f8.1,2i5)')
c     &         "SSBOND_E_FORM",totT,t_bath,i,j
c        endif
#endif
#endif
        dyn_ssbond_ij(ici,icj)=eij
      else if (.not.havebond .and. dyn_ssbond_ij(ici,icj).lt.1.0d300) 
     &then
        dyn_ssbond_ij(ici,icj)=1.0d300
#ifndef CLUST
#ifndef WHAM
c        write(iout,'(a15,f12.2,f8.1,2i5)')
c     &       "SSBOND_E_BREAK",totT,t_bath,i,j
#endif
#endif
      endif

c-------TESTING CODE
      if (checkstop) then
        if (jcheck.eq.0) write(iout,'(a,3f15.8,$)')
     &       "CHECKSTOP",rij,eij,ed
        echeck(jcheck)=eij
      endif
      enddo
      if (checkstop) then
        write(iout,'(f15.8)')(echeck(1)-echeck(-1))*0.5d0/deps
      endif
      enddo
      if (checkstop) then
        transgrad=.true.
        checkstop=.false.
      endif
c-------END TESTING CODE
            gg_lipi(3)=ssgradlipi*eij
            gg_lipj(3)=ssgradlipj*eij

      do k=1,3
        dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
        dcosom2(k)=(dc_norm(k,nres+j)-om2*erij(k))/rij
      enddo
      do k=1,3
        gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
      enddo
      do k=1,3
        gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
     &       +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
     &       +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
        gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
     &       +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
     &       +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
      enddo
cgrad      do k=i,j-1
cgrad        do l=1,3
cgrad          gvdwc(l,k)=gvdwc(l,k)+gg(l)
cgrad        enddo
cgrad      enddo

      do l=1,3
        gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
        gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
      enddo

      return
      end
C-----------------------------------------------------------------------------

      double precision function h_base(x,deriv)
c     A smooth function going 0->1 in range [0,1]
c     It should NOT be called outside range [0,1], it will not work there.
      implicit none

c     Input arguments
      double precision x

c     Output arguments
      double precision deriv

c     Local variables
      double precision xsq


c     Two parabolas put together.  First derivative zero at extrema
c$$$      if (x.lt.0.5D0) then
c$$$        h_base=2.0D0*x*x
c$$$        deriv=4.0D0*x
c$$$      else
c$$$        deriv=1.0D0-x
c$$$        h_base=1.0D0-2.0D0*deriv*deriv
c$$$        deriv=4.0D0*deriv
c$$$      endif

c     Third degree polynomial.  First derivative zero at extrema
      h_base=x*x*(3.0d0-2.0d0*x)
      deriv=6.0d0*x*(1.0d0-x)

c     Fifth degree polynomial.  First and second derivatives zero at extrema
c$$$      xsq=x*x
c$$$      h_base=x*xsq*(6.0d0*xsq-15.0d0*x+10.0d0)
c$$$      deriv=x-1.0d0
c$$$      deriv=deriv*deriv
c$$$      deriv=30.0d0*xsq*deriv

      return
      end
c----------------------------------------------------------------------------
      subroutine dyn_set_nss
c     Adjust nss and other relevant variables based on dyn_ssbond_ij
c      implicit none

c     Includes
      include 'DIMENSIONS'
#ifdef MPI
      include "mpif.h"
#endif
      include 'COMMON.SBRIDGE'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
C      include 'COMMON.SETUP'
#ifndef CLUST
#ifndef WHAM
C      include 'COMMON.MD'
#endif
#endif

c     Local variables
      double precision emin
      integer i,j,imin
      integer diff,allflag(maxdim_cont),allnss,
     &     allihpb(maxdim_cont),alljhpb(maxdim_cont),
     &     newnss,newihpb(maxdim_cont),newjhpb(maxdim_cont)
      logical found
      integer i_newnss(1024),displ(0:1024)
      integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss

      allnss=0
      do i=1,ns-1
        do j=i+1,ns
          if (dyn_ssbond_ij(i,j).lt.1.0d300) then
            allnss=allnss+1
            allflag(allnss)=0
            allihpb(allnss)=i
            alljhpb(allnss)=j
          endif
        enddo
      enddo

cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)

 1    emin=1.0d300
      do i=1,allnss
        if (allflag(i).eq.0 .and.
     &       dyn_ssbond_ij(allihpb(i),alljhpb(i)).lt.emin) then
          emin=dyn_ssbond_ij(allihpb(i),alljhpb(i))
          imin=i
        endif
      enddo
      if (emin.lt.1.0d300) then
        allflag(imin)=1
        do i=1,allnss
          if (allflag(i).eq.0 .and.
     &         (allihpb(i).eq.allihpb(imin) .or.
     &         alljhpb(i).eq.allihpb(imin) .or.
     &         allihpb(i).eq.alljhpb(imin) .or.
     &         alljhpb(i).eq.alljhpb(imin))) then
            allflag(i)=-1
          endif
        enddo
        goto 1
      endif

cmc      write(iout,*)"ALLNSS ",allnss,(allihpb(i),alljhpb(i),i=1,allnss)

      newnss=0
      do i=1,allnss
        if (allflag(i).eq.1) then
          newnss=newnss+1
          newihpb(newnss)=allihpb(i)
          newjhpb(newnss)=alljhpb(i)
        endif
      enddo

#ifdef MPI
      if (nfgtasks.gt.1)then

        call MPI_Reduce(newnss,g_newnss,1,
     &    MPI_INTEGER,MPI_SUM,king,FG_COMM,IERR)
        call MPI_Gather(newnss,1,MPI_INTEGER,
     &                  i_newnss,1,MPI_INTEGER,king,FG_COMM,IERR)
        displ(0)=0
        do i=1,nfgtasks-1,1
          displ(i)=i_newnss(i-1)+displ(i-1)
        enddo
        call MPI_Gatherv(newihpb,newnss,MPI_INTEGER,
     &                   g_newihpb,i_newnss,displ,MPI_INTEGER,
     &                   king,FG_COMM,IERR)     
        call MPI_Gatherv(newjhpb,newnss,MPI_INTEGER,
     &                   g_newjhpb,i_newnss,displ,MPI_INTEGER,
     &                   king,FG_COMM,IERR)     
        if(fg_rank.eq.0) then
c         print *,'g_newnss',g_newnss
c         print *,'g_newihpb',(g_newihpb(i),i=1,g_newnss)
c         print *,'g_newjhpb',(g_newjhpb(i),i=1,g_newnss)
         newnss=g_newnss  
         do i=1,newnss
          newihpb(i)=g_newihpb(i)
          newjhpb(i)=g_newjhpb(i)
         enddo
        endif
      endif
#endif

      diff=newnss-nss

cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)

      do i=1,nss
        found=.false.
        do j=1,newnss
          if (idssb(i).eq.newihpb(j) .and.
     &         jdssb(i).eq.newjhpb(j)) found=.true.
        enddo
#ifndef CLUST
#ifndef WHAM
c        if (.not.found.and.fg_rank.eq.0) 
c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
c     &       "SSBOND_BREAK",totT,t_bath,idssb(i),jdssb(i)
#endif
#endif
      enddo

      do i=1,newnss
        found=.false.
        do j=1,nss
          if (newihpb(i).eq.idssb(j) .and.
     &         newjhpb(i).eq.jdssb(j)) found=.true.
        enddo
#ifndef CLUST
#ifndef WHAM
c        if (.not.found.and.fg_rank.eq.0) 
c     &      write(iout,'(a15,f12.2,f8.1,2i5)')
c     &       "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i)
#endif
#endif
      enddo

      nss=newnss
      do i=1,nss
        idssb(i)=newihpb(i)
        jdssb(i)=newjhpb(i)
      enddo

      return
      end

c----------------------------------------------------------------------------

#ifdef SSREAD
#ifdef WHAM
      subroutine read_ssHist
      implicit none

c     Includes
      include 'DIMENSIONS'
      include "DIMENSIONS.FREE"
      include 'COMMON.FREE'

c     Local variables
      integer i,j
      character*80 controlcard

      do i=1,dyn_nssHist
        call card_concat(controlcard,.true.)
        read(controlcard,*)
     &       dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0))
      enddo

      return
      end
#endif
#endif
c$$$C----------------------------------------------------------------------------
      subroutine triple_ssbond_ene(resi,resj,resk,eij)
      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
C      include 'COMMON.MD'
#endif
#endif

c     External functions
      double precision h_base
      external h_base

c     Input arguments
      integer resi,resj,resk

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
      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)

      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.
      eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
C second case jth atom is center
      eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
C the third case kth atom is the center
      eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
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)+2.0*btriss*(rij+rik))
     &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))  
            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)+2.0*btriss*(rij+rik))
     &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
            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)+2.0*btriss*(rij+rjk))-
     &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
            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