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 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) itypj=itype(j) xj=c(1,nres+j)-c(1,nres+i) yj=c(2,nres+j)-c(2,nres+i) zj=c(3,nres+j)-c(3,nres+i) 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(itypi,itypj) ljA=ljA*aa(itypi,itypj) ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(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(itypi,itypj)/aa(itypi,itypj) 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(itypi,itypj) e2=fac*bb(itypi,itypj) eij=eps1*eps2rt*eps3rt*(e1+e2) 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 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 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(itypi,itypj)/aa(itypi,itypj) d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt) d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt- + alf12/eps3rt) d_ljm(1)=d_ljm(1)*(0.5D0*eps2rt_om1/eps2rt-alf1/eps3rt) f1=(rij-ljxm)/(xm-ljxm) f2=(rij-xm)/(ljxm-xm) h1=h_base(f1,hd1) h2=h_base(f2,hd2) eij=Ht*h1+ljm*h2 delta_inv=1.0d0/(ljxm-xm) deltasq_inv=delta_inv*delta_inv fac=Ht*hd1-ljm*hd2 fac1=deltasq_inv*fac*(ljxm-rij) fac2=deltasq_inv*fac*(rij-xm) ed=delta_inv*(ljm*hd2-Ht*hd1) eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1) eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2) eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3) endif c-------END FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE c-------SECOND METHOD, CONTINUOUS SECOND DERIVATIVE c$$$ ssd=rij-ssXs c$$$ ljd=rij-ljXs c$$$ fac1=rij-ljxm c$$$ fac2=rij-ssxm c$$$ c$$$ d_ljB(1)=ljB*(eps2rt_om1/eps2rt-2.0d0*alf1/eps3rt) c$$$ d_ljB(2)=ljB*(eps2rt_om2/eps2rt+2.0d0*alf2/eps3rt) c$$$ d_ljB(3)=ljB*(eps1_om12+eps2rt_om12/eps2rt-2.0d0*alf12/eps3rt) c$$$ c$$$ ssm=ssC-0.25D0*ssB*ssB/ssA c$$$ d_ssm(1)=0.5D0*akct*ssB/ssA c$$$ d_ssm(2)=2.0D0*akth*deltat2-om1*omega-d_ssm(1) c$$$ d_ssm(1)=-2.0D0*akth*deltat1-om2*omega+d_ssm(1) c$$$ d_ssm(3)=omega c$$$ c$$$ ljm=-0.25D0*bb(itypi,itypj)/aa(itypi,itypj) c$$$ do k=1,3 c$$$ d_ljm(k)=ljm*d_ljB(k) c$$$ enddo c$$$ ljm=ljm*ljB c$$$ c$$$ ss=ssA*ssd*ssd+ssB*ssd+ssC c$$$ d_ss(0)=2.0d0*ssA*ssd+ssB c$$$ d_ss(2)=akct*ssd c$$$ d_ss(1)=-d_ss(2)-2.0d0*akth*deltat1-om2*omega c$$$ d_ss(2)=d_ss(2)+2.0d0*akth*deltat2-om1*omega c$$$ d_ss(3)=omega c$$$ c$$$ ljf=bb(itypi,itypj)/aa(itypi,itypj) c$$$ ljf=9.0d0*ljf*(-0.5d0*ljf)**(1.0d0/3.0d0) c$$$ d_ljf(0)=ljf*2.0d0*ljB*fac1 c$$$ do k=1,3 c$$$ d_ljf(k)=d_ljm(k)+ljf*(d_ljB(k)*fac1*fac1- c$$$ & 2.0d0*ljB*fac1*d_ljxm(k)) c$$$ enddo c$$$ ljf=ljm+ljf*ljB*fac1*fac1 c$$$ c$$$ f1=(rij-ljxm)/(ssxm-ljxm) c$$$ f2=(rij-ssxm)/(ljxm-ssxm) c$$$ h1=h_base(f1,hd1) c$$$ h2=h_base(f2,hd2) c$$$ eij=ss*h1+ljf*h2 c$$$ delta_inv=1.0d0/(ljxm-ssxm) c$$$ deltasq_inv=delta_inv*delta_inv c$$$ fac=ljf*hd2-ss*hd1 c$$$ ed=d_ss(0)*h1+d_ljf(0)*h2+delta_inv*fac c$$$ eom1=d_ss(1)*h1+d_ljf(1)*h2+deltasq_inv*fac* c$$$ & (fac1*d_ssxm(1)-fac2*(d_ljxm(1))) c$$$ eom2=d_ss(2)*h1+d_ljf(2)*h2+deltasq_inv*fac* c$$$ & (fac1*d_ssxm(2)-fac2*(d_ljxm(2))) c$$$ eom12=d_ss(3)*h1+d_ljf(3)*h2+deltasq_inv*fac* c$$$ & (fac1*d_ssxm(3)-fac2*(d_ljxm(3))) c$$$ c$$$ havebond=.false. c$$$ if (ed.gt.0.0d0) havebond=.true. c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE endif if (havebond) then #ifndef CLUST #ifndef WHAM c if (dyn_ssbond_ij(i,j).eq.1.0d300) then c write(iout,'(a15,f12.2,f8.1,2i5)') c & "SSBOND_E_FORM",totT,t_bath,i,j c endif #endif #endif write(iout,*), 'DYN_SS_BOND',i,j,eij 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#ifdef CLUST c#ifdef WHAM integer max_fg_procs,maxprocs parameter (max_fg_procs=128) parameter (maxprocs=128) c#endif c#endif cc include 'COMMON.SETUP' #ifndef CLUST #ifndef WHAM include 'COMMON.MD' include 'COMMON.SETUP' #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(max_fg_procs),displ(0:max_fg_procs) 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 if (.not.found.and.fg_rank.eq.0) & write(iout,'(a15,f12.2,f8.1,2i5)') & "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 if (.not.found.and.fg_rank.eq.0) & write(iout,'(a15,f12.2,f8.1,2i5)') & "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 NIEWIADOM 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-----------------------------------------------------------------------------