end
C-----------------------------------------------------------------------------
-
subroutine dyn_ssbond_ene(resi,resj,eij)
-c implicit none
-
-c Includes
+ implicit none
include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include 'COMMON.NAMES'
#ifndef CLUST
#ifndef WHAM
-C include 'COMMON.MD'
+ include 'COMMON.MD'
#endif
#endif
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)
common /sschecks/ checkstop,transgrad
integer icheck,nicheck,jcheck,njcheck
- double precision echeck(-1:1),deps,ssx0,ljx0
+ 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)
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
+ 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)
- 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
+ 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
+ & +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
+ & +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)
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
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
+ 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
havebond=.true.
ssd=rij-ssXs
eij=ssA*ssd*ssd+ssB*ssd+ssC
-C write(iout,*) 'TU?2',ssc,ssd
+ 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
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)
+ 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)
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)
+ 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)
c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
endif
- write(iout,*) 'havebond',havebond
+
if (havebond) then
#ifndef CLUST
#ifndef WHAM
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
+ 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)')
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
gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
enddo
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ 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)
+ 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 enddo
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ 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)
return
end
-
c----------------------------------------------------------------------------
-
subroutine dyn_set_nss
c Adjust nss and other relevant variables based on dyn_ssbond_ij
c implicit none
c Local variables
double precision emin
integer i,j,imin
- integer diff,allflag(maxdim),allnss,
- & allihpb(maxdim),alljhpb(maxdim),
- & newnss,newihpb(maxdim),newjhpb(maxdim)
+ 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),g_newjhpb(maxdim),g_newnss
+ integer g_newihpb(maxdim_cont),g_newjhpb(maxdim_cont),g_newnss
allnss=0
- do i=1,nres-1
- do j=i+1,nres
+ 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
end
#endif
#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)
+c$$$C----------------------------------------------------------------------------
+ subroutine triple_ssbond_ene(resi,resj,resk,eij)
include 'DIMENSIONS'
include 'COMMON.SBRIDGE'
include 'COMMON.CHAIN'
projects / unres.git / blobdiff