X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fcluster%2Fwham%2Fsrc-HCD%2FssMD.F;fp=source%2Fcluster%2Fwham%2Fsrc-HCD%2FssMD.F;h=9b2908fe24bb3252e9e934267ad378b2f3ef5db8;hb=57038e4bdff4cc9534106b25bfbd4b9a844d47fd;hp=9c23fe0ae64922c40964af01c25a67c62b4f6fb9;hpb=32caa3b64eb94b90fa9fd402b77263ea89efffa1;p=unres.git diff --git a/source/cluster/wham/src-HCD/ssMD.F b/source/cluster/wham/src-HCD/ssMD.F index 9c23fe0..9b2908f 100644 --- a/source/cluster/wham/src-HCD/ssMD.F +++ b/source/cluster/wham/src-HCD/ssMD.F @@ -3,7 +3,6 @@ c---------------------------------------------------------------------------- c implicit none c Includes - implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.VAR' @@ -83,12 +82,8 @@ ct rij=ran_number(rmin,rmax) end C----------------------------------------------------------------------------- - subroutine dyn_ssbond_ene(resi,resj,eij) -c implicit none - -c Includes - implicit real*8 (a-h,o-z) + implicit none include 'DIMENSIONS' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN' @@ -98,9 +93,10 @@ c Includes 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 @@ -130,7 +126,10 @@ c integer itypi,itypj,k,l double precision omega,delta_inv,deltasq_inv,fac1,fac2 c-------FIRST METHOD double precision xm,d_xm(1:3) - integer xshift,yshift,zshift + 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) @@ -141,129 +140,50 @@ c-------TESTING CODE 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) 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=dmod(xi,boxxsize) - if (xi.lt.0) xi=xi+boxxsize - yi=dmod(yi,boxysize) - if (yi.lt.0) yi=yi+boxysize - zi=dmod(zi,boxzsize) - if (zi.lt.0) zi=zi+boxzsize + xi=c(1,nres+i) + yi=c(2,nres+i) + zi=c(3,nres+i) + call to_box(xi,yi,zi) C define scaling factor for lipids C if (positi.le.0) positi=positi+boxzsize C print *,i C first for peptide groups c for each residue check if it is in lipid or lipid water border area - 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- - & ((positi-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-positi)/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 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=dmod(xj,boxxsize) - if (xj.lt.0) xj=xj+boxxsize - yj=dmod(yj,boxysize) - if (yj.lt.0) yj=yj+boxysize - zj=dmod(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- - & ((positi-bordlipbot)/lipbufthick) -C lipbufthick is thickenes of lipid buffore - sslipj=sscalelip(fracinbuf) - ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) - sslipj=sscalelip(fracinbuf) - ssgradlipj=sscagradlip(fracinbuf)/lipbufthick - else - sslipj=1.0d0 - ssgradlipj=0.0 - endif - else - sslipj=0.0d0 - ssgradlipj=0.0 - endif + xj=c(1,nres+j) + yj=c(2,nres+j) + zj=c(3,nres+j) + call to_box(xj,yj,zj) + call lipid_layer(xj,yj,zj,sslipj,ssgradlipj) aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0 & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0 - - dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - xj_safe=xj - yj_safe=yj - zj_safe=zj - subchap=0 - do xshift=-1,1 - do yshift=-1,1 - do zshift=-1,1 - xj=xj_safe+xshift*boxxsize - yj=yj_safe+yshift*boxysize - zj=zj_safe+zshift*boxzsize - dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 - if(dist_temp.lt.dist_init) then - dist_init=dist_temp - xj_temp=xj - yj_temp=yj - zj_temp=zj - subchap=1 - endif - enddo - enddo - enddo - if (subchap.eq.1) then - xj=xj_temp-xi - yj=yj_temp-yi - zj=zj_temp-zi - else - xj=xj_safe-xi - yj=yj_safe-yi - zj=zj_safe-zi - endif - + 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) @@ -281,6 +201,8 @@ C lipbufthick is thickenes of lipid buffore 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 @@ -360,15 +282,15 @@ c-------END TESTING CODE 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 @@ -377,8 +299,9 @@ C write(iout,*) eij,'TU?1' 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 @@ -413,13 +336,14 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE 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) @@ -436,13 +360,14 @@ C write(iout,*) eij,'TU?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) @@ -508,7 +433,7 @@ c$$$ if (ed.gt.0.0d0) havebond=.true. c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE endif -C write(iout,*) 'havebond',havebond + if (havebond) then #ifndef CLUST #ifndef WHAM @@ -518,9 +443,10 @@ 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 + 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)') @@ -545,6 +471,8 @@ c-------TESTING CODE 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 @@ -554,10 +482,10 @@ c-------END TESTING CODE 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 @@ -568,13 +496,12 @@ cgrad 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) @@ -615,15 +542,12 @@ 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 - implicit real*8 (a-h,o-z) include 'DIMENSIONS' #ifdef MPI include "mpif.h" @@ -641,16 +565,16 @@ C include 'COMMON.MD' 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 @@ -768,1268 +692,34 @@ c & "SSBOND_FORM",totT,t_bath,newihpb(i),newjhpb(i) return end +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) - implicit real*8 (a-h,o-z) +#ifdef SSREAD +#ifdef WHAM + subroutine read_ssHist + implicit none + +c Includes + include 'DIMENSIONS' + include "DIMENSIONS.FREE" + include 'COMMON.FREE' + +c Local variables + integer i,j + character*80 controlcard + + do i=1,dyn_nssHist + call card_concat(controlcard,.true.) + read(controlcard,*) + & dyn_ssHist(i,0),(dyn_ssHist(i,j),j=1,2*dyn_ssHist(i,0)) + enddo + + return + end +#endif +#endif +c$$$C---------------------------------------------------------------------------- + subroutine triple_ssbond_ene(resi,resj,resk,eij) include 'DIMENSIONS' include 'COMMON.SBRIDGE' include 'COMMON.CHAIN'