+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)
+c implicit none
+
+c Includes
+ 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
+
+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)
+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
+c-------END TESTING CODE
+
+
+ i=resi
+ j=resj
+
+ 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)
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+ itypj=itype(j)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zj-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ 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=xj-xi
+ yj=yj-yi
+ zj=zj-zi
+ 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
+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)
+C write(iout,*) eij,'TU?1'
+ eps2der=eij*eps3rt
+ eps3der=eij*eps2rt
+ eij=eij*eps2rt*eps3rt
+
+ sigder=-sig/sigsq
+ e1=e1*eps1*eps2rt**2*eps3rt**2
+ ed=-expon*(e1+eij)/ljd
+ sigder=ed*sigder
+ 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
+C write(iout,*) 'TU?2',ssc,ssd
+ ed=2*akcm*ssd+akct*deltat12
+ 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
+C write(iout,*) eij,'TU?3'
+ 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)
+ 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
+C write(iout,*) 'TU?4',ssA
+ 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)
+ 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
+ write(iout,*) 'havebond',havebond
+ 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(i,j)=eij
+ else if (.not.havebond .and. dyn_ssbond_ij(i,j).lt.1.0d300) then
+ dyn_ssbond_ij(i,j)=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
+
+ 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)
+ & +(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)
+ & +(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)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ 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),allnss,
+ & allihpb(maxdim),alljhpb(maxdim),
+ & newnss,newihpb(maxdim),newjhpb(maxdim)
+ logical found
+ integer i_newnss(1024),displ(0:1024)
+ integer g_newihpb(maxdim),g_newjhpb(maxdim),g_newnss
+
+ allnss=0
+ do i=1,nres-1
+ do j=i+1,nres
+ 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 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
+
+c----------------------------------------------------------------------------
+
+
+C-----------------------------------------------------------------------------
+C-----------------------------------------------------------------------------
+C-----------------------------------------------------------------------------
+C-----------------------------------------------------------------------------
+C-----------------------------------------------------------------------------
+C-----------------------------------------------------------------------------
+C-----------------------------------------------------------------------------
+
+c$$$c-----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine ss_relax(i_in,j_in)
+c$$$ implicit none
+c$$$
+c$$$c Includes
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.IOUNITS'
+c$$$ include 'COMMON.INTERACT'
+c$$$
+c$$$c Input arguments
+c$$$ integer i_in,j_in
+c$$$
+c$$$c Local variables
+c$$$ integer i,iretcode,nfun_sc
+c$$$ logical scfail
+c$$$ double precision var(maxvar),e_sc,etot
+c$$$
+c$$$
+c$$$ mask_r=.true.
+c$$$ do i=nnt,nct
+c$$$ mask_side(i)=0
+c$$$ enddo
+c$$$ mask_side(i_in)=1
+c$$$ mask_side(j_in)=1
+c$$$
+c$$$c Minimize the two selected side-chains
+c$$$ call overlap_sc(scfail) ! Better not fail!
+c$$$ call minimize_sc(e_sc,var,iretcode,nfun_sc)
+c$$$
+c$$$ mask_r=.false.
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$c-------------------------------------------------------------
+c$$$
+c$$$ subroutine minimize_sc(etot_sc,iretcode,nfun)
+c$$$c Minimize side-chains only, starting from geom but without modifying
+c$$$c bond lengths.
+c$$$c If mask_r is already set, only the selected side-chains are minimized,
+c$$$c otherwise all side-chains are minimized keeping the backbone frozen.
+c$$$ implicit none
+c$$$
+c$$$c Includes
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.IOUNITS'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.GEO'
+c$$$ include 'COMMON.MINIM'
+c$$$ integer icall
+c$$$ common /srutu/ icall
+c$$$
+c$$$c Output arguments
+c$$$ double precision etot_sc
+c$$$ integer iretcode,nfun
+c$$$
+c$$$c External functions/subroutines
+c$$$ external func_sc,grad_sc,fdum
+c$$$
+c$$$c Local variables
+c$$$ integer liv,lv
+c$$$ parameter (liv=60,lv=(77+maxvar*(maxvar+17)/2))
+c$$$ integer iv(liv)
+c$$$ double precision rdum(1)
+c$$$ double precision d(maxvar),v(1:lv),x(maxvar),xx(maxvar)
+c$$$ integer idum(1)
+c$$$ integer i,nvar_restr
+c$$$
+c$$$
+c$$$cmc start_minim=.true.
+c$$$ call deflt(2,iv,liv,lv,v)
+c$$$* 12 means fresh start, dont call deflt
+c$$$ iv(1)=12
+c$$$* max num of fun calls
+c$$$ if (maxfun.eq.0) maxfun=500
+c$$$ iv(17)=maxfun
+c$$$* max num of iterations
+c$$$ if (maxmin.eq.0) maxmin=1000
+c$$$ iv(18)=maxmin
+c$$$* controls output
+c$$$ iv(19)=1
+c$$$* selects output unit
+c$$$ iv(21)=0
+c$$$c iv(21)=iout ! DEBUG
+c$$$c iv(21)=8 ! DEBUG
+c$$$* 1 means to print out result
+c$$$ iv(22)=0
+c$$$c iv(22)=1 ! DEBUG
+c$$$* 1 means to print out summary stats
+c$$$ iv(23)=0
+c$$$c iv(23)=1 ! DEBUG
+c$$$* 1 means to print initial x and d
+c$$$ iv(24)=0
+c$$$c iv(24)=1 ! DEBUG
+c$$$* min val for v(radfac) default is 0.1
+c$$$ v(24)=0.1D0
+c$$$* max val for v(radfac) default is 4.0
+c$$$ v(25)=2.0D0
+c$$$c v(25)=4.0D0
+c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
+c$$$* the sumsl default is 0.1
+c$$$ v(26)=0.1D0
+c$$$* false conv if (act fnctn decrease) .lt. v(34)
+c$$$* the sumsl default is 100*machep
+c$$$ v(34)=v(34)/100.0D0
+c$$$* absolute convergence
+c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
+c$$$ v(31)=tolf
+c$$$* relative convergence
+c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-1
+c$$$ v(32)=rtolf
+c$$$* controls initial step size
+c$$$ v(35)=1.0D-1
+c$$$* large vals of d correspond to small components of step
+c$$$ do i=1,nphi
+c$$$ d(i)=1.0D-1
+c$$$ enddo
+c$$$ do i=nphi+1,nvar
+c$$$ d(i)=1.0D-1
+c$$$ enddo
+c$$$
+c$$$ call geom_to_var(nvar,x)
+c$$$ IF (mask_r) THEN
+c$$$ do i=1,nres ! Just in case...
+c$$$ mask_phi(i)=0
+c$$$ mask_theta(i)=0
+c$$$ enddo
+c$$$ call x2xx(x,xx,nvar_restr)
+c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
+c$$$ & iv,liv,lv,v,idum,rdum,fdum)
+c$$$ call xx2x(x,xx)
+c$$$ ELSE
+c$$$c When minimizing ALL side-chains, etotal_sc is a little
+c$$$c faster if we don't set mask_r
+c$$$ do i=1,nres
+c$$$ mask_phi(i)=0
+c$$$ mask_theta(i)=0
+c$$$ mask_side(i)=1
+c$$$ enddo
+c$$$ call x2xx(x,xx,nvar_restr)
+c$$$ call sumsl(nvar_restr,d,xx,func_sc,grad_sc,
+c$$$ & iv,liv,lv,v,idum,rdum,fdum)
+c$$$ call xx2x(x,xx)
+c$$$ ENDIF
+c$$$ call var_to_geom(nvar,x)
+c$$$ call chainbuild_sc
+c$$$ etot_sc=v(10)
+c$$$ iretcode=iv(1)
+c$$$ nfun=iv(6)
+c$$$ return
+c$$$ end
+c$$$
+c$$$C--------------------------------------------------------------------------
+c$$$
+c$$$ subroutine chainbuild_sc
+c$$$ implicit none
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.INTERACT'
+c$$$
+c$$$c Local variables
+c$$$ integer i
+c$$$
+c$$$
+c$$$ do i=nnt,nct
+c$$$ if (.not.mask_r .or. mask_side(i).eq.1) then
+c$$$ call locate_side_chain(i)
+c$$$ endif
+c$$$ enddo
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$C--------------------------------------------------------------------------
+c$$$
+c$$$ subroutine func_sc(n,x,nf,f,uiparm,urparm,ufparm)
+c$$$ implicit none
+c$$$
+c$$$c Includes
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.DERIV'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.MINIM'
+c$$$ include 'COMMON.IOUNITS'
+c$$$
+c$$$c Input arguments
+c$$$ integer n
+c$$$ double precision x(maxvar)
+c$$$ double precision ufparm
+c$$$ external ufparm
+c$$$
+c$$$c Input/Output arguments
+c$$$ integer nf
+c$$$ integer uiparm(1)
+c$$$ double precision urparm(1)
+c$$$
+c$$$c Output arguments
+c$$$ double precision f
+c$$$
+c$$$c Local variables
+c$$$ double precision energia(0:n_ene)
+c$$$#ifdef OSF
+c$$$c Variables used to intercept NaNs
+c$$$ double precision x_sum
+c$$$ integer i_NAN
+c$$$#endif
+c$$$
+c$$$
+c$$$ nfl=nf
+c$$$ icg=mod(nf,2)+1
+c$$$
+c$$$#ifdef OSF
+c$$$c Intercept NaNs in the coordinates, before calling etotal_sc
+c$$$ x_sum=0.D0
+c$$$ do i_NAN=1,n
+c$$$ x_sum=x_sum+x(i_NAN)
+c$$$ enddo
+c$$$c Calculate the energy only if the coordinates are ok
+c$$$ if ((.not.(x_sum.lt.0.D0)) .and. (.not.(x_sum.ge.0.D0))) then
+c$$$ write(iout,*)" *** func_restr_sc : Found NaN in coordinates"
+c$$$ f=1.0D+77
+c$$$ nf=0
+c$$$ else
+c$$$#endif
+c$$$
+c$$$ call var_to_geom_restr(n,x)
+c$$$ call zerograd
+c$$$ call chainbuild_sc
+c$$$ call etotal_sc(energia(0))
+c$$$ f=energia(0)
+c$$$ if (energia(1).eq.1.0D20 .or. energia(0).eq.1.0D99) nf=0
+c$$$
+c$$$#ifdef OSF
+c$$$ endif
+c$$$#endif
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$c-------------------------------------------------------
+c$$$
+c$$$ subroutine grad_sc(n,x,nf,g,uiparm,urparm,ufparm)
+c$$$ implicit none
+c$$$
+c$$$c Includes
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.DERIV'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.INTERACT'
+c$$$ include 'COMMON.MINIM'
+c$$$
+c$$$c Input arguments
+c$$$ integer n
+c$$$ double precision x(maxvar)
+c$$$ double precision ufparm
+c$$$ external ufparm
+c$$$
+c$$$c Input/Output arguments
+c$$$ integer nf
+c$$$ integer uiparm(1)
+c$$$ double precision urparm(1)
+c$$$
+c$$$c Output arguments
+c$$$ double precision g(maxvar)
+c$$$
+c$$$c Local variables
+c$$$ double precision f,gphii,gthetai,galphai,gomegai
+c$$$ integer ig,ind,i,j,k,igall,ij
+c$$$
+c$$$
+c$$$ icg=mod(nf,2)+1
+c$$$ if (nf-nfl+1) 20,30,40
+c$$$ 20 call func_sc(n,x,nf,f,uiparm,urparm,ufparm)
+c$$$c write (iout,*) 'grad 20'
+c$$$ if (nf.eq.0) return
+c$$$ goto 40
+c$$$ 30 call var_to_geom_restr(n,x)
+c$$$ call chainbuild_sc
+c$$$C
+c$$$C Evaluate the derivatives of virtual bond lengths and SC vectors in variables.
+c$$$C
+c$$$ 40 call cartder
+c$$$C
+c$$$C Convert the Cartesian gradient into internal-coordinate gradient.
+c$$$C
+c$$$
+c$$$ ig=0
+c$$$ ind=nres-2
+c$$$ do i=2,nres-2
+c$$$ IF (mask_phi(i+2).eq.1) THEN
+c$$$ gphii=0.0D0
+c$$$ do j=i+1,nres-1
+c$$$ ind=ind+1
+c$$$ do k=1,3
+c$$$ gphii=gphii+dcdv(k+3,ind)*gradc(k,j,icg)
+c$$$ gphii=gphii+dxdv(k+3,ind)*gradx(k,j,icg)
+c$$$ enddo
+c$$$ enddo
+c$$$ ig=ig+1
+c$$$ g(ig)=gphii
+c$$$ ELSE
+c$$$ ind=ind+nres-1-i
+c$$$ ENDIF
+c$$$ enddo
+c$$$
+c$$$
+c$$$ ind=0
+c$$$ do i=1,nres-2
+c$$$ IF (mask_theta(i+2).eq.1) THEN
+c$$$ ig=ig+1
+c$$$ gthetai=0.0D0
+c$$$ do j=i+1,nres-1
+c$$$ ind=ind+1
+c$$$ do k=1,3
+c$$$ gthetai=gthetai+dcdv(k,ind)*gradc(k,j,icg)
+c$$$ gthetai=gthetai+dxdv(k,ind)*gradx(k,j,icg)
+c$$$ enddo
+c$$$ enddo
+c$$$ g(ig)=gthetai
+c$$$ ELSE
+c$$$ ind=ind+nres-1-i
+c$$$ ENDIF
+c$$$ enddo
+c$$$
+c$$$ do i=2,nres-1
+c$$$ if (itype(i).ne.10) then
+c$$$ IF (mask_side(i).eq.1) THEN
+c$$$ ig=ig+1
+c$$$ galphai=0.0D0
+c$$$ do k=1,3
+c$$$ galphai=galphai+dxds(k,i)*gradx(k,i,icg)
+c$$$ enddo
+c$$$ g(ig)=galphai
+c$$$ ENDIF
+c$$$ endif
+c$$$ enddo
+c$$$
+c$$$
+c$$$ do i=2,nres-1
+c$$$ if (itype(i).ne.10) then
+c$$$ IF (mask_side(i).eq.1) THEN
+c$$$ ig=ig+1
+c$$$ gomegai=0.0D0
+c$$$ do k=1,3
+c$$$ gomegai=gomegai+dxds(k+3,i)*gradx(k,i,icg)
+c$$$ enddo
+c$$$ g(ig)=gomegai
+c$$$ ENDIF
+c$$$ endif
+c$$$ enddo
+c$$$
+c$$$C
+c$$$C Add the components corresponding to local energy terms.
+c$$$C
+c$$$
+c$$$ ig=0
+c$$$ igall=0
+c$$$ do i=4,nres
+c$$$ igall=igall+1
+c$$$ if (mask_phi(i).eq.1) then
+c$$$ ig=ig+1
+c$$$ g(ig)=g(ig)+gloc(igall,icg)
+c$$$ endif
+c$$$ enddo
+c$$$
+c$$$ do i=3,nres
+c$$$ igall=igall+1
+c$$$ if (mask_theta(i).eq.1) then
+c$$$ ig=ig+1
+c$$$ g(ig)=g(ig)+gloc(igall,icg)
+c$$$ endif
+c$$$ enddo
+c$$$
+c$$$ do ij=1,2
+c$$$ do i=2,nres-1
+c$$$ if (itype(i).ne.10) then
+c$$$ igall=igall+1
+c$$$ if (mask_side(i).eq.1) then
+c$$$ ig=ig+1
+c$$$ g(ig)=g(ig)+gloc(igall,icg)
+c$$$ endif
+c$$$ endif
+c$$$ enddo
+c$$$ enddo
+c$$$
+c$$$cd do i=1,ig
+c$$$cd write (iout,'(a2,i5,a3,f25.8)') 'i=',i,' g=',g(i)
+c$$$cd enddo
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$C-----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine etotal_sc(energy_sc)
+c$$$ implicit none
+c$$$
+c$$$c Includes
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.INTERACT'
+c$$$ include 'COMMON.DERIV'
+c$$$ include 'COMMON.FFIELD'
+c$$$
+c$$$c Output arguments
+c$$$ double precision energy_sc(0:n_ene)
+c$$$
+c$$$c Local variables
+c$$$ double precision evdw,escloc
+c$$$ integer i,j
+c$$$
+c$$$
+c$$$ do i=1,n_ene
+c$$$ energy_sc(i)=0.0D0
+c$$$ enddo
+c$$$
+c$$$ if (mask_r) then
+c$$$ call egb_sc(evdw)
+c$$$ call esc_sc(escloc)
+c$$$ else
+c$$$ call egb(evdw)
+c$$$ call esc(escloc)
+c$$$ endif
+c$$$
+c$$$ if (evdw.eq.1.0D20) then
+c$$$ energy_sc(0)=evdw
+c$$$ else
+c$$$ energy_sc(0)=wsc*evdw+wscloc*escloc
+c$$$ endif
+c$$$ energy_sc(1)=evdw
+c$$$ energy_sc(12)=escloc
+c$$$
+c$$$C
+c$$$C Sum up the components of the Cartesian gradient.
+c$$$C
+c$$$ do i=1,nct
+c$$$ do j=1,3
+c$$$ gradx(j,i,icg)=wsc*gvdwx(j,i)
+c$$$ enddo
+c$$$ enddo
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$C-----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine egb_sc(evdw)
+c$$$C
+c$$$C This subroutine calculates the interaction energy of nonbonded side chains
+c$$$C assuming the Gay-Berne potential of interaction.
+c$$$C
+c$$$ implicit real*8 (a-h,o-z)
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.GEO'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.LOCAL'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.DERIV'
+c$$$ include 'COMMON.NAMES'
+c$$$ include 'COMMON.INTERACT'
+c$$$ include 'COMMON.IOUNITS'
+c$$$ include 'COMMON.CALC'
+c$$$ include 'COMMON.CONTROL'
+c$$$ logical lprn
+c$$$ evdw=0.0D0
+c$$$ energy_dec=.false.
+c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+c$$$ evdw=0.0D0
+c$$$ lprn=.false.
+c$$$c if (icall.eq.0) lprn=.false.
+c$$$ ind=0
+c$$$ do i=iatsc_s,iatsc_e
+c$$$ itypi=itype(i)
+c$$$ itypi1=itype(i+1)
+c$$$ xi=c(1,nres+i)
+c$$$ yi=c(2,nres+i)
+c$$$ zi=c(3,nres+i)
+c$$$ dxi=dc_norm(1,nres+i)
+c$$$ dyi=dc_norm(2,nres+i)
+c$$$ dzi=dc_norm(3,nres+i)
+c$$$c dsci_inv=dsc_inv(itypi)
+c$$$ dsci_inv=vbld_inv(i+nres)
+c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
+c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
+c$$$C
+c$$$C Calculate SC interaction energy.
+c$$$C
+c$$$ do iint=1,nint_gr(i)
+c$$$ do j=istart(i,iint),iend(i,iint)
+c$$$ IF (mask_side(j).eq.1.or.mask_side(i).eq.1) THEN
+c$$$ ind=ind+1
+c$$$ itypj=itype(j)
+c$$$c dscj_inv=dsc_inv(itypj)
+c$$$ dscj_inv=vbld_inv(j+nres)
+c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
+c$$$c & 1.0d0/vbld(j+nres)
+c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+c$$$ sig0ij=sigma(itypi,itypj)
+c$$$ chi1=chi(itypi,itypj)
+c$$$ chi2=chi(itypj,itypi)
+c$$$ chi12=chi1*chi2
+c$$$ chip1=chip(itypi)
+c$$$ chip2=chip(itypj)
+c$$$ chip12=chip1*chip2
+c$$$ alf1=alp(itypi)
+c$$$ alf2=alp(itypj)
+c$$$ alf12=0.5D0*(alf1+alf2)
+c$$$C For diagnostics only!!!
+c$$$c chi1=0.0D0
+c$$$c chi2=0.0D0
+c$$$c chi12=0.0D0
+c$$$c chip1=0.0D0
+c$$$c chip2=0.0D0
+c$$$c chip12=0.0D0
+c$$$c alf1=0.0D0
+c$$$c alf2=0.0D0
+c$$$c alf12=0.0D0
+c$$$ xj=c(1,nres+j)-xi
+c$$$ yj=c(2,nres+j)-yi
+c$$$ zj=c(3,nres+j)-zi
+c$$$ dxj=dc_norm(1,nres+j)
+c$$$ dyj=dc_norm(2,nres+j)
+c$$$ dzj=dc_norm(3,nres+j)
+c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
+c$$$c write (iout,*) "j",j," dc_norm",
+c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
+c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+c$$$ rij=dsqrt(rrij)
+c$$$C Calculate angle-dependent terms of energy and contributions to their
+c$$$C derivatives.
+c$$$ call sc_angular
+c$$$ sigsq=1.0D0/sigsq
+c$$$ sig=sig0ij*dsqrt(sigsq)
+c$$$ rij_shift=1.0D0/rij-sig+sig0ij
+c$$$c for diagnostics; uncomment
+c$$$c rij_shift=1.2*sig0ij
+c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
+c$$$ if (rij_shift.le.0.0D0) then
+c$$$ evdw=1.0D20
+c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+c$$$cd & restyp(itypi),i,restyp(itypj),j,
+c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
+c$$$ return
+c$$$ endif
+c$$$ sigder=-sig*sigsq
+c$$$c---------------------------------------------------------------
+c$$$ rij_shift=1.0D0/rij_shift
+c$$$ fac=rij_shift**expon
+c$$$ e1=fac*fac*aa(itypi,itypj)
+c$$$ e2=fac*bb(itypi,itypj)
+c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+c$$$ eps2der=evdwij*eps3rt
+c$$$ eps3der=evdwij*eps2rt
+c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
+c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+c$$$ evdwij=evdwij*eps2rt*eps3rt
+c$$$ evdw=evdw+evdwij
+c$$$ if (lprn) then
+c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+c$$$ & restyp(itypi),i,restyp(itypj),j,
+c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
+c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+c$$$ & evdwij
+c$$$ endif
+c$$$
+c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
+c$$$ & 'evdw',i,j,evdwij
+c$$$
+c$$$C Calculate gradient components.
+c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
+c$$$ fac=-expon*(e1+evdwij)*rij_shift
+c$$$ sigder=fac*sigder
+c$$$ fac=rij*fac
+c$$$c fac=0.0d0
+c$$$C Calculate the radial part of the gradient
+c$$$ gg(1)=xj*fac
+c$$$ gg(2)=yj*fac
+c$$$ gg(3)=zj*fac
+c$$$C Calculate angular part of the gradient.
+c$$$ call sc_grad
+c$$$ ENDIF
+c$$$ enddo ! j
+c$$$ enddo ! iint
+c$$$ enddo ! i
+c$$$ energy_dec=.false.
+c$$$ return
+c$$$ end
+c$$$
+c$$$c-----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine esc_sc(escloc)
+c$$$C Calculate the local energy of a side chain and its derivatives in the
+c$$$C corresponding virtual-bond valence angles THETA and the spherical angles
+c$$$C ALPHA and OMEGA.
+c$$$ implicit real*8 (a-h,o-z)
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.GEO'
+c$$$ include 'COMMON.LOCAL'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.INTERACT'
+c$$$ include 'COMMON.DERIV'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.IOUNITS'
+c$$$ include 'COMMON.NAMES'
+c$$$ include 'COMMON.FFIELD'
+c$$$ include 'COMMON.CONTROL'
+c$$$ double precision x(3),dersc(3),xemp(3),dersc0(3),dersc1(3),
+c$$$ & ddersc0(3),ddummy(3),xtemp(3),temp(3)
+c$$$ common /sccalc/ time11,time12,time112,theti,it,nlobit
+c$$$ delta=0.02d0*pi
+c$$$ escloc=0.0D0
+c$$$c write (iout,'(a)') 'ESC'
+c$$$ do i=loc_start,loc_end
+c$$$ IF (mask_side(i).eq.1) THEN
+c$$$ it=itype(i)
+c$$$ if (it.eq.10) goto 1
+c$$$ nlobit=nlob(it)
+c$$$c print *,'i=',i,' it=',it,' nlobit=',nlobit
+c$$$c write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
+c$$$ theti=theta(i+1)-pipol
+c$$$ x(1)=dtan(theti)
+c$$$ x(2)=alph(i)
+c$$$ x(3)=omeg(i)
+c$$$
+c$$$ if (x(2).gt.pi-delta) then
+c$$$ xtemp(1)=x(1)
+c$$$ xtemp(2)=pi-delta
+c$$$ xtemp(3)=x(3)
+c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
+c$$$ xtemp(2)=pi
+c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
+c$$$ call spline1(x(2),pi-delta,delta,escloci0,escloci1,dersc0(2),
+c$$$ & escloci,dersc(2))
+c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
+c$$$ & ddersc0(1),dersc(1))
+c$$$ call spline2(x(2),pi-delta,delta,dersc0(3),dersc1(3),
+c$$$ & ddersc0(3),dersc(3))
+c$$$ xtemp(2)=pi-delta
+c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
+c$$$ xtemp(2)=pi
+c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
+c$$$ call spline1(x(2),pi-delta,delta,esclocbi0,esclocbi1,
+c$$$ & dersc0(2),esclocbi,dersc02)
+c$$$ call spline2(x(2),pi-delta,delta,dersc0(1),dersc1(1),
+c$$$ & dersc12,dersc01)
+c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
+c$$$ dersc0(1)=dersc01
+c$$$ dersc0(2)=dersc02
+c$$$ dersc0(3)=0.0d0
+c$$$ do k=1,3
+c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
+c$$$ enddo
+c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
+c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
+c$$$c & esclocbi,ss,ssd
+c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
+c$$$c escloci=esclocbi
+c$$$c write (iout,*) escloci
+c$$$ else if (x(2).lt.delta) then
+c$$$ xtemp(1)=x(1)
+c$$$ xtemp(2)=delta
+c$$$ xtemp(3)=x(3)
+c$$$ call enesc(xtemp,escloci0,dersc0,ddersc0,.true.)
+c$$$ xtemp(2)=0.0d0
+c$$$ call enesc(xtemp,escloci1,dersc1,ddummy,.false.)
+c$$$ call spline1(x(2),delta,-delta,escloci0,escloci1,dersc0(2),
+c$$$ & escloci,dersc(2))
+c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
+c$$$ & ddersc0(1),dersc(1))
+c$$$ call spline2(x(2),delta,-delta,dersc0(3),dersc1(3),
+c$$$ & ddersc0(3),dersc(3))
+c$$$ xtemp(2)=delta
+c$$$ call enesc_bound(xtemp,esclocbi0,dersc0,dersc12,.true.)
+c$$$ xtemp(2)=0.0d0
+c$$$ call enesc_bound(xtemp,esclocbi1,dersc1,chuju,.false.)
+c$$$ call spline1(x(2),delta,-delta,esclocbi0,esclocbi1,
+c$$$ & dersc0(2),esclocbi,dersc02)
+c$$$ call spline2(x(2),delta,-delta,dersc0(1),dersc1(1),
+c$$$ & dersc12,dersc01)
+c$$$ dersc0(1)=dersc01
+c$$$ dersc0(2)=dersc02
+c$$$ dersc0(3)=0.0d0
+c$$$ call splinthet(x(2),0.5d0*delta,ss,ssd)
+c$$$ do k=1,3
+c$$$ dersc(k)=ss*dersc(k)+(1.0d0-ss)*dersc0(k)
+c$$$ enddo
+c$$$ dersc(2)=dersc(2)+ssd*(escloci-esclocbi)
+c$$$c write (iout,*) 'i=',i,x(2)*rad2deg,escloci0,escloci,
+c$$$c & esclocbi,ss,ssd
+c$$$ escloci=ss*escloci+(1.0d0-ss)*esclocbi
+c$$$c write (iout,*) escloci
+c$$$ else
+c$$$ call enesc(x,escloci,dersc,ddummy,.false.)
+c$$$ endif
+c$$$
+c$$$ escloc=escloc+escloci
+c$$$ if (energy_dec) write (iout,'(a6,i,0pf7.3)')
+c$$$ & 'escloc',i,escloci
+c$$$c write (iout,*) 'i=',i,' escloci=',escloci,' dersc=',dersc
+c$$$
+c$$$ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
+c$$$ & wscloc*dersc(1)
+c$$$ gloc(ialph(i,1),icg)=wscloc*dersc(2)
+c$$$ gloc(ialph(i,1)+nside,icg)=wscloc*dersc(3)
+c$$$ 1 continue
+c$$$ ENDIF
+c$$$ enddo
+c$$$ return
+c$$$ end
+c$$$
+c$$$C-----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine egb_ij(i_sc,j_sc,evdw)
+c$$$C
+c$$$C This subroutine calculates the interaction energy of nonbonded side chains
+c$$$C assuming the Gay-Berne potential of interaction.
+c$$$C
+c$$$ implicit real*8 (a-h,o-z)
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.GEO'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.LOCAL'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.DERIV'
+c$$$ include 'COMMON.NAMES'
+c$$$ include 'COMMON.INTERACT'
+c$$$ include 'COMMON.IOUNITS'
+c$$$ include 'COMMON.CALC'
+c$$$ include 'COMMON.CONTROL'
+c$$$ logical lprn
+c$$$ evdw=0.0D0
+c$$$ energy_dec=.false.
+c$$$c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+c$$$ evdw=0.0D0
+c$$$ lprn=.false.
+c$$$ ind=0
+c$$$c$$$ do i=iatsc_s,iatsc_e
+c$$$ i=i_sc
+c$$$ itypi=itype(i)
+c$$$ itypi1=itype(i+1)
+c$$$ xi=c(1,nres+i)
+c$$$ yi=c(2,nres+i)
+c$$$ zi=c(3,nres+i)
+c$$$ dxi=dc_norm(1,nres+i)
+c$$$ dyi=dc_norm(2,nres+i)
+c$$$ dzi=dc_norm(3,nres+i)
+c$$$c dsci_inv=dsc_inv(itypi)
+c$$$ dsci_inv=vbld_inv(i+nres)
+c$$$c write (iout,*) "i",i,dsc_inv(itypi),dsci_inv,1.0d0/vbld(i+nres)
+c$$$c write (iout,*) "dcnori",dxi*dxi+dyi*dyi+dzi*dzi
+c$$$C
+c$$$C Calculate SC interaction energy.
+c$$$C
+c$$$c$$$ do iint=1,nint_gr(i)
+c$$$c$$$ do j=istart(i,iint),iend(i,iint)
+c$$$ j=j_sc
+c$$$ ind=ind+1
+c$$$ itypj=itype(j)
+c$$$c dscj_inv=dsc_inv(itypj)
+c$$$ dscj_inv=vbld_inv(j+nres)
+c$$$c write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
+c$$$c & 1.0d0/vbld(j+nres)
+c$$$c write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
+c$$$ sig0ij=sigma(itypi,itypj)
+c$$$ chi1=chi(itypi,itypj)
+c$$$ chi2=chi(itypj,itypi)
+c$$$ chi12=chi1*chi2
+c$$$ chip1=chip(itypi)
+c$$$ chip2=chip(itypj)
+c$$$ chip12=chip1*chip2
+c$$$ alf1=alp(itypi)
+c$$$ alf2=alp(itypj)
+c$$$ alf12=0.5D0*(alf1+alf2)
+c$$$C For diagnostics only!!!
+c$$$c chi1=0.0D0
+c$$$c chi2=0.0D0
+c$$$c chi12=0.0D0
+c$$$c chip1=0.0D0
+c$$$c chip2=0.0D0
+c$$$c chip12=0.0D0
+c$$$c alf1=0.0D0
+c$$$c alf2=0.0D0
+c$$$c alf12=0.0D0
+c$$$ xj=c(1,nres+j)-xi
+c$$$ yj=c(2,nres+j)-yi
+c$$$ zj=c(3,nres+j)-zi
+c$$$ dxj=dc_norm(1,nres+j)
+c$$$ dyj=dc_norm(2,nres+j)
+c$$$ dzj=dc_norm(3,nres+j)
+c$$$c write (iout,*) "dcnorj",dxi*dxi+dyi*dyi+dzi*dzi
+c$$$c write (iout,*) "j",j," dc_norm",
+c$$$c & dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
+c$$$ rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+c$$$ rij=dsqrt(rrij)
+c$$$C Calculate angle-dependent terms of energy and contributions to their
+c$$$C derivatives.
+c$$$ call sc_angular
+c$$$ sigsq=1.0D0/sigsq
+c$$$ sig=sig0ij*dsqrt(sigsq)
+c$$$ rij_shift=1.0D0/rij-sig+sig0ij
+c$$$c for diagnostics; uncomment
+c$$$c rij_shift=1.2*sig0ij
+c$$$C I hate to put IF's in the loops, but here don't have another choice!!!!
+c$$$ if (rij_shift.le.0.0D0) then
+c$$$ evdw=1.0D20
+c$$$cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+c$$$cd & restyp(itypi),i,restyp(itypj),j,
+c$$$cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
+c$$$ return
+c$$$ endif
+c$$$ sigder=-sig*sigsq
+c$$$c---------------------------------------------------------------
+c$$$ rij_shift=1.0D0/rij_shift
+c$$$ fac=rij_shift**expon
+c$$$ e1=fac*fac*aa(itypi,itypj)
+c$$$ e2=fac*bb(itypi,itypj)
+c$$$ evdwij=eps1*eps2rt*eps3rt*(e1+e2)
+c$$$ eps2der=evdwij*eps3rt
+c$$$ eps3der=evdwij*eps2rt
+c$$$c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
+c$$$c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
+c$$$ evdwij=evdwij*eps2rt*eps3rt
+c$$$ evdw=evdw+evdwij
+c$$$ if (lprn) then
+c$$$ sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+c$$$ epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+c$$$ write (iout,'(2(a3,i3,2x),17(0pf7.3))')
+c$$$ & restyp(itypi),i,restyp(itypj),j,
+c$$$ & epsi,sigm,chi1,chi2,chip1,chip2,
+c$$$ & eps1,eps2rt**2,eps3rt**2,sig,sig0ij,
+c$$$ & om1,om2,om12,1.0D0/rij,1.0D0/rij_shift,
+c$$$ & evdwij
+c$$$ endif
+c$$$
+c$$$ if (energy_dec) write (iout,'(a6,2i,0pf7.3)')
+c$$$ & 'evdw',i,j,evdwij
+c$$$
+c$$$C Calculate gradient components.
+c$$$ e1=e1*eps1*eps2rt**2*eps3rt**2
+c$$$ fac=-expon*(e1+evdwij)*rij_shift
+c$$$ sigder=fac*sigder
+c$$$ fac=rij*fac
+c$$$c fac=0.0d0
+c$$$C Calculate the radial part of the gradient
+c$$$ gg(1)=xj*fac
+c$$$ gg(2)=yj*fac
+c$$$ gg(3)=zj*fac
+c$$$C Calculate angular part of the gradient.
+c$$$ call sc_grad
+c$$$c$$$ enddo ! j
+c$$$c$$$ enddo ! iint
+c$$$c$$$ enddo ! i
+c$$$ energy_dec=.false.
+c$$$ return
+c$$$ end
+c$$$
+c$$$C-----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine perturb_side_chain(i,angle)
+c$$$ implicit none
+c$$$
+c$$$c Includes
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.GEO'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.LOCAL'
+c$$$ include 'COMMON.IOUNITS'
+c$$$
+c$$$c External functions
+c$$$ external ran_number
+c$$$ double precision ran_number
+c$$$
+c$$$c Input arguments
+c$$$ integer i
+c$$$ double precision angle ! In degrees
+c$$$
+c$$$c Local variables
+c$$$ integer i_sc
+c$$$ double precision rad_ang,rand_v(3),length,cost,sint
+c$$$
+c$$$
+c$$$ i_sc=i+nres
+c$$$ rad_ang=angle*deg2rad
+c$$$
+c$$$ length=0.0
+c$$$ do while (length.lt.0.01)
+c$$$ rand_v(1)=ran_number(0.01D0,1.0D0)
+c$$$ rand_v(2)=ran_number(0.01D0,1.0D0)
+c$$$ rand_v(3)=ran_number(0.01D0,1.0D0)
+c$$$ length=rand_v(1)*rand_v(1)+rand_v(2)*rand_v(2)+
+c$$$ + rand_v(3)*rand_v(3)
+c$$$ length=sqrt(length)
+c$$$ rand_v(1)=rand_v(1)/length
+c$$$ rand_v(2)=rand_v(2)/length
+c$$$ rand_v(3)=rand_v(3)/length
+c$$$ cost=rand_v(1)*dc_norm(1,i_sc)+rand_v(2)*dc_norm(2,i_sc)+
+c$$$ + rand_v(3)*dc_norm(3,i_sc)
+c$$$ length=1.0D0-cost*cost
+c$$$ if (length.lt.0.0D0) length=0.0D0
+c$$$ length=sqrt(length)
+c$$$ rand_v(1)=rand_v(1)-cost*dc_norm(1,i_sc)
+c$$$ rand_v(2)=rand_v(2)-cost*dc_norm(2,i_sc)
+c$$$ rand_v(3)=rand_v(3)-cost*dc_norm(3,i_sc)
+c$$$ enddo
+c$$$ rand_v(1)=rand_v(1)/length
+c$$$ rand_v(2)=rand_v(2)/length
+c$$$ rand_v(3)=rand_v(3)/length
+c$$$
+c$$$ cost=dcos(rad_ang)
+c$$$ sint=dsin(rad_ang)
+c$$$ dc(1,i_sc)=vbld(i_sc)*(dc_norm(1,i_sc)*cost+rand_v(1)*sint)
+c$$$ dc(2,i_sc)=vbld(i_sc)*(dc_norm(2,i_sc)*cost+rand_v(2)*sint)
+c$$$ dc(3,i_sc)=vbld(i_sc)*(dc_norm(3,i_sc)*cost+rand_v(3)*sint)
+c$$$ dc_norm(1,i_sc)=dc(1,i_sc)*vbld_inv(i_sc)
+c$$$ dc_norm(2,i_sc)=dc(2,i_sc)*vbld_inv(i_sc)
+c$$$ dc_norm(3,i_sc)=dc(3,i_sc)*vbld_inv(i_sc)
+c$$$ c(1,i_sc)=c(1,i)+dc(1,i_sc)
+c$$$ c(2,i_sc)=c(2,i)+dc(2,i_sc)
+c$$$ c(3,i_sc)=c(3,i)+dc(3,i_sc)
+c$$$
+c$$$ call chainbuild_cart
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$c----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine ss_relax3(i_in,j_in)
+c$$$ implicit none
+c$$$
+c$$$c Includes
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.IOUNITS'
+c$$$ include 'COMMON.INTERACT'
+c$$$
+c$$$c External functions
+c$$$ external ran_number
+c$$$ double precision ran_number
+c$$$
+c$$$c Input arguments
+c$$$ integer i_in,j_in
+c$$$
+c$$$c Local variables
+c$$$ double precision energy_sc(0:n_ene),etot
+c$$$ double precision org_dc(3),org_dc_norm(3),org_c(3)
+c$$$ double precision ang_pert,rand_fact,exp_fact,beta
+c$$$ integer n,i_pert,i
+c$$$ logical notdone
+c$$$
+c$$$
+c$$$ beta=1.0D0
+c$$$
+c$$$ mask_r=.true.
+c$$$ do i=nnt,nct
+c$$$ mask_side(i)=0
+c$$$ enddo
+c$$$ mask_side(i_in)=1
+c$$$ mask_side(j_in)=1
+c$$$
+c$$$ call etotal_sc(energy_sc)
+c$$$ etot=energy_sc(0)
+c$$$c write(iout,'(a,3d15.5)')" SS_MC_START ",energy_sc(0),
+c$$$c + energy_sc(1),energy_sc(12)
+c$$$
+c$$$ notdone=.true.
+c$$$ n=0
+c$$$ do while (notdone)
+c$$$ if (mod(n,2).eq.0) then
+c$$$ i_pert=i_in
+c$$$ else
+c$$$ i_pert=j_in
+c$$$ endif
+c$$$ n=n+1
+c$$$
+c$$$ do i=1,3
+c$$$ org_dc(i)=dc(i,i_pert+nres)
+c$$$ org_dc_norm(i)=dc_norm(i,i_pert+nres)
+c$$$ org_c(i)=c(i,i_pert+nres)
+c$$$ enddo
+c$$$ ang_pert=ran_number(0.0D0,3.0D0)
+c$$$ call perturb_side_chain(i_pert,ang_pert)
+c$$$ call etotal_sc(energy_sc)
+c$$$ exp_fact=exp(beta*(etot-energy_sc(0)))
+c$$$ rand_fact=ran_number(0.0D0,1.0D0)
+c$$$ if (rand_fact.lt.exp_fact) then
+c$$$c write(iout,'(a,3d15.5)')" SS_MC_ACCEPT ",energy_sc(0),
+c$$$c + energy_sc(1),energy_sc(12)
+c$$$ etot=energy_sc(0)
+c$$$ else
+c$$$c write(iout,'(a,3d15.5)')" SS_MC_REJECT ",energy_sc(0),
+c$$$c + energy_sc(1),energy_sc(12)
+c$$$ do i=1,3
+c$$$ dc(i,i_pert+nres)=org_dc(i)
+c$$$ dc_norm(i,i_pert+nres)=org_dc_norm(i)
+c$$$ c(i,i_pert+nres)=org_c(i)
+c$$$ enddo
+c$$$ endif
+c$$$
+c$$$ if (n.eq.10000.or.etot.lt.30.0D0) notdone=.false.
+c$$$ enddo
+c$$$
+c$$$ mask_r=.false.
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$c----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine ss_relax2(etot,iretcode,nfun,i_in,j_in)
+c$$$ implicit none
+c$$$ include 'DIMENSIONS'
+c$$$ integer liv,lv
+c$$$ parameter (liv=60,lv=(77+maxres6*(maxres6+17)/2))
+c$$$*********************************************************************
+c$$$* OPTIMIZE sets up SUMSL or DFP and provides a simple interface for *
+c$$$* the calling subprogram. *
+c$$$* when d(i)=1.0, then v(35) is the length of the initial step, *
+c$$$* calculated in the usual pythagorean way. *
+c$$$* absolute convergence occurs when the function is within v(31) of *
+c$$$* zero. unless you know the minimum value in advance, abs convg *
+c$$$* is probably not useful. *
+c$$$* relative convergence is when the model predicts that the function *
+c$$$* will decrease by less than v(32)*abs(fun). *
+c$$$*********************************************************************
+c$$$ include 'COMMON.IOUNITS'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.GEO'
+c$$$ include 'COMMON.MINIM'
+c$$$ include 'COMMON.CHAIN'
+c$$$
+c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
+c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
+c$$$ + orig_ss_dist(maxres2,maxres2)
+c$$$
+c$$$ double precision etot
+c$$$ integer iretcode,nfun,i_in,j_in
+c$$$
+c$$$ external dist
+c$$$ double precision dist
+c$$$ external ss_func,fdum
+c$$$ double precision ss_func,fdum
+c$$$
+c$$$ integer iv(liv),uiparm(2)
+c$$$ double precision v(lv),x(maxres6),d(maxres6),rdum
+c$$$ integer i,j,k
+c$$$
+c$$$
+c$$$ call deflt(2,iv,liv,lv,v)
+c$$$* 12 means fresh start, dont call deflt
+c$$$ iv(1)=12
+c$$$* max num of fun calls
+c$$$ if (maxfun.eq.0) maxfun=500
+c$$$ iv(17)=maxfun
+c$$$* max num of iterations
+c$$$ if (maxmin.eq.0) maxmin=1000
+c$$$ iv(18)=maxmin
+c$$$* controls output
+c$$$ iv(19)=2
+c$$$* selects output unit
+c$$$c iv(21)=iout
+c$$$ iv(21)=0
+c$$$* 1 means to print out result
+c$$$ iv(22)=0
+c$$$* 1 means to print out summary stats
+c$$$ iv(23)=0
+c$$$* 1 means to print initial x and d
+c$$$ iv(24)=0
+c$$$* min val for v(radfac) default is 0.1
+c$$$ v(24)=0.1D0
+c$$$* max val for v(radfac) default is 4.0
+c$$$ v(25)=2.0D0
+c$$$c v(25)=4.0D0
+c$$$* check false conv if (act fnctn decrease) .lt. v(26)*(exp decrease)
+c$$$* the sumsl default is 0.1
+c$$$ v(26)=0.1D0
+c$$$* false conv if (act fnctn decrease) .lt. v(34)
+c$$$* the sumsl default is 100*machep
+c$$$ v(34)=v(34)/100.0D0
+c$$$* absolute convergence
+c$$$ if (tolf.eq.0.0D0) tolf=1.0D-4
+c$$$ v(31)=tolf
+c$$$ v(31)=1.0D-1
+c$$$* relative convergence
+c$$$ if (rtolf.eq.0.0D0) rtolf=1.0D-4
+c$$$ v(32)=rtolf
+c$$$ v(32)=1.0D-1
+c$$$* controls initial step size
+c$$$ v(35)=1.0D-1
+c$$$* large vals of d correspond to small components of step
+c$$$ do i=1,6*nres
+c$$$ d(i)=1.0D0
+c$$$ enddo
+c$$$
+c$$$ do i=0,2*nres
+c$$$ do j=1,3
+c$$$ orig_ss_dc(j,i)=dc(j,i)
+c$$$ enddo
+c$$$ enddo
+c$$$ call geom_to_var(nvar,orig_ss_var)
+c$$$
+c$$$ do i=1,nres
+c$$$ do j=i,nres
+c$$$ orig_ss_dist(j,i)=dist(j,i)
+c$$$ orig_ss_dist(j+nres,i)=dist(j+nres,i)
+c$$$ orig_ss_dist(j,i+nres)=dist(j,i+nres)
+c$$$ orig_ss_dist(j+nres,i+nres)=dist(j+nres,i+nres)
+c$$$ enddo
+c$$$ enddo
+c$$$
+c$$$ k=0
+c$$$ do i=1,nres-1
+c$$$ do j=1,3
+c$$$ k=k+1
+c$$$ x(k)=dc(j,i)
+c$$$ enddo
+c$$$ enddo
+c$$$ do i=2,nres-1
+c$$$ if (ialph(i,1).gt.0) then
+c$$$ do j=1,3
+c$$$ k=k+1
+c$$$ x(k)=dc(j,i+nres)
+c$$$ enddo
+c$$$ endif
+c$$$ enddo
+c$$$
+c$$$ uiparm(1)=i_in
+c$$$ uiparm(2)=j_in
+c$$$ call smsno(k,d,x,ss_func,iv,liv,lv,v,uiparm,rdum,fdum)
+c$$$ etot=v(10)
+c$$$ iretcode=iv(1)
+c$$$ nfun=iv(6)+iv(30)
+c$$$
+c$$$ k=0
+c$$$ do i=1,nres-1
+c$$$ do j=1,3
+c$$$ k=k+1
+c$$$ dc(j,i)=x(k)
+c$$$ enddo
+c$$$ enddo
+c$$$ do i=2,nres-1
+c$$$ if (ialph(i,1).gt.0) then
+c$$$ do j=1,3
+c$$$ k=k+1
+c$$$ dc(j,i+nres)=x(k)
+c$$$ enddo
+c$$$ endif
+c$$$ enddo
+c$$$ call chainbuild_cart
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$C-----------------------------------------------------------------------------
+c$$$
+c$$$ subroutine ss_func(n,x,nf,f,uiparm,urparm,ufparm)
+c$$$ implicit none
+c$$$ include 'DIMENSIONS'
+c$$$ include 'COMMON.DERIV'
+c$$$ include 'COMMON.IOUNITS'
+c$$$ include 'COMMON.VAR'
+c$$$ include 'COMMON.CHAIN'
+c$$$ include 'COMMON.INTERACT'
+c$$$ include 'COMMON.SBRIDGE'
+c$$$
+c$$$ double precision orig_ss_dc,orig_ss_var,orig_ss_dist
+c$$$ common /orig_ss/ orig_ss_dc(3,0:maxres2),orig_ss_var(maxvar),
+c$$$ + orig_ss_dist(maxres2,maxres2)
+c$$$
+c$$$ integer n
+c$$$ double precision x(maxres6)
+c$$$ integer nf
+c$$$ double precision f
+c$$$ integer uiparm(2)
+c$$$ real*8 urparm(1)
+c$$$ external ufparm
+c$$$ double precision ufparm
+c$$$
+c$$$ external dist
+c$$$ double precision dist
+c$$$
+c$$$ integer i,j,k,ss_i,ss_j
+c$$$ double precision tempf,var(maxvar)
+c$$$
+c$$$
+c$$$ ss_i=uiparm(1)
+c$$$ ss_j=uiparm(2)
+c$$$ f=0.0D0
+c$$$
+c$$$ k=0
+c$$$ do i=1,nres-1
+c$$$ do j=1,3
+c$$$ k=k+1
+c$$$ dc(j,i)=x(k)
+c$$$ enddo
+c$$$ enddo
+c$$$ do i=2,nres-1
+c$$$ if (ialph(i,1).gt.0) then
+c$$$ do j=1,3
+c$$$ k=k+1
+c$$$ dc(j,i+nres)=x(k)
+c$$$ enddo
+c$$$ endif
+c$$$ enddo
+c$$$ call chainbuild_cart
+c$$$
+c$$$ call geom_to_var(nvar,var)
+c$$$
+c$$$c Constraints on all angles
+c$$$ do i=1,nvar
+c$$$ tempf=var(i)-orig_ss_var(i)
+c$$$ f=f+tempf*tempf
+c$$$ enddo
+c$$$
+c$$$c Constraints on all distances
+c$$$ do i=1,nres-1
+c$$$ if (i.gt.1) then
+c$$$ tempf=dist(i+nres,i)-orig_ss_dist(i+nres,i)
+c$$$ f=f+tempf*tempf
+c$$$ endif
+c$$$ do j=i+1,nres
+c$$$ tempf=dist(j,i)-orig_ss_dist(j,i)
+c$$$ if (tempf.lt.0.0D0 .or. j.eq.i+1) f=f+tempf*tempf
+c$$$ tempf=dist(j+nres,i)-orig_ss_dist(j+nres,i)
+c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
+c$$$ tempf=dist(j,i+nres)-orig_ss_dist(j,i+nres)
+c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
+c$$$ tempf=dist(j+nres,i+nres)-orig_ss_dist(j+nres,i+nres)
+c$$$ if (tempf.lt.0.0D0) f=f+tempf*tempf
+c$$$ enddo
+c$$$ enddo
+c$$$
+c$$$c Constraints for the relevant CYS-CYS
+c$$$ tempf=dist(nres+ss_i,nres+ss_j)-8.0D0
+c$$$ f=f+tempf*tempf
+c$$$CCCCCCCCCCCCCCCCC ADD SOME ANGULAR STUFF
+c$$$
+c$$$c$$$ if (nf.ne.nfl) then
+c$$$c$$$ write(iout,'(a,i10,2d15.5)')"IN DIST_FUNC (NF,F,DIST)",nf,
+c$$$c$$$ + f,dist(5+nres,14+nres)
+c$$$c$$$ endif
+c$$$
+c$$$ nfl=nf
+c$$$
+c$$$ return
+c$$$ end
+c$$$
+c$$$C-----------------------------------------------------------------------------
+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
projects / unres.git / commitdiff