zmiany w galezi multichain

diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN
index bf5f688..777cc43 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CHAIN
+++ b/source/unres/src_MD-M/COMMON.CHAIN
@@ -14,5 +14,5 @@
      & nsup,nstart_sup,nstart_seq,
      & chain_length,iprzes,tabperm(maxperm,maxsym),nperm
       common /from_zscore/ nz_start,nz_end,iz_sc
-      double precision boxxsize,boxysize,boxzsize
-      common /box/  boxxsize,boxysize,boxzsize 
+      double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad
+      common /box/  boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad

diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV
index 524d72a..a5e35e2 100644 (file)
--- a/source/unres/src_MD-M/COMMON.DERIV
+++ b/source/unres/src_MD-M/COMMON.DERIV
@@ -9,7 +9,7 @@
      & gvdwc(3,maxres),gelc(3,maxres),gelc_long(3,maxres),
      & gvdwpp(3,maxres),gvdwc_scpp(3,maxres),
      & gradx_scp(3,maxres),gvdwc_scp(3,maxres),ghpbx(3,maxres),
-     & ghpbc(3,maxres),gloc(maxvar,2),gradcorr(3,maxres),
+     & ghpbc(3,maxres),gloc(0:maxvar,2),gradcorr(3,maxres),
      & gradcorr_long(3,maxres),gradcorr5_long(3,maxres),
      & gradcorr6_long(3,maxres),gcorr6_turn_long(3,maxres),
      & gradxorr(3,maxres),gradcorr5(3,maxres),gradcorr6(3,maxres),

diff --git a/source/unres/src_MD-M/chainbuild.F b/source/unres/src_MD-M/chainbuild.F
index 45a1a53..c4970d3 100644 (file)
--- a/source/unres/src_MD-M/chainbuild.F
+++ b/source/unres/src_MD-M/chainbuild.F
@@ -12,9 +12,121 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.NAMES'
       include 'COMMON.INTERACT'
-      logical lprn
+      double precision e1(3),e2(3),e3(3)
+      logical lprn,perbox,fail
 C Set lprn=.true. for debugging
       lprn = .false.
+      perbox=.false.
+      fail=.false.
+       if (perbox) then
+      cost=dcos(theta(3))
+      sint=dsin(theta(3))
+      print *,'before refsys'
+      call refsys(2,3,4,e1,e2,e3,fail)
+      print *,'after refsys'
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+      dc(1,0)=c(1,1)
+      dc(2,0)=c(2,1)
+      dc(3,0)=c(3,1)
+      print *,'dc',dc(1,0),dc(2,0),dc(3,0)
+      dc(1,1)=c(1,2)-c(1,1)
+      dc(2,1)=c(2,2)-c(2,1)
+      dc(3,1)=c(3,2)-c(3,1)
+      dc(1,2)=c(1,3)-c(1,2)
+      dc(2,2)=c(2,3)-c(2,2)
+      dc(3,2)=c(3,3)-c(3,2)
+      t(1,1,1)=e1(1)
+      t(1,2,1)=e1(2)
+      t(1,3,1)=e1(3)
+      t(2,1,1)=e2(1)
+      t(2,2,1)=e2(2)
+      t(2,3,1)=e2(3)
+      t(3,1,1)=e3(1)
+      t(3,2,1)=e3(2)
+      t(3,3,1)=e3(3)
+      veclen=0.0d0
+      do i=1,3
+       veclen=veclen+(c(i,2)-c(i,1))**2
+      enddo
+      veclen=sqrt(veclen)
+      r(1,1,1)= 1.0D0
+      r(1,2,1)= 0.0D0
+      r(1,3,1)= 0.0D0
+      r(2,1,1)= 0.0D0
+      r(2,2,1)= 1.0D0
+      r(2,3,1)= 0.0D0
+      r(3,1,1)= 0.0D0
+      r(3,2,1)= 0.0D0
+      r(3,3,1)= 1.0D0
+      do i=1,3
+        do j=1,3
+          rt(i,j,1)=t(i,j,1)
+        enddo
+      enddo
+      do i=1,3
+        do j=1,3
+          prod(i,j,1)=0.0D0
+          prod(i,j,2)=t(i,j,1)
+        enddo
+        prod(i,i,1)=1.0D0
+      enddo
+        call locate_side_chain(2)
+      do i=4,nres
+#ifdef OSF
+      theti=theta(i)
+      if (theti.ne.theti) theti=100.0
+      phii=phi(i)
+      if (phii.ne.phii) phii=180.0
+#else
+      theti=theta(i)
+      phii=phi(i)
+#endif
+      cost=dcos(theti)
+      sint=dsin(theti)
+      cosphi=dcos(phii)
+      sinphi=dsin(phii)
+* Define the matrices of the rotation about the virtual-bond valence angles
+* theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
+* program), R(i,j,k), and, the cumulative matrices of rotation RT
+      t(1,1,i-2)=-cost
+      t(1,2,i-2)=-sint
+      t(1,3,i-2)= 0.0D0
+      t(2,1,i-2)=-sint
+      t(2,2,i-2)= cost
+      t(2,3,i-2)= 0.0D0
+      t(3,1,i-2)= 0.0D0
+      t(3,2,i-2)= 0.0D0
+      t(3,3,i-2)= 1.0D0
+      r(1,1,i-2)= 1.0D0
+      r(1,2,i-2)= 0.0D0
+      r(1,3,i-2)= 0.0D0
+      r(2,1,i-2)= 0.0D0
+      r(2,2,i-2)=-cosphi
+      r(2,3,i-2)= sinphi
+      r(3,1,i-2)= 0.0D0
+      r(3,2,i-2)= sinphi
+      r(3,3,i-2)= cosphi
+      rt(1,1,i-2)=-cost
+      rt(1,2,i-2)=-sint
+      rt(1,3,i-2)=0.0D0
+      rt(2,1,i-2)=sint*cosphi
+      rt(2,2,i-2)=-cost*cosphi
+      rt(2,3,i-2)=sinphi
+      rt(3,1,i-2)=-sint*sinphi
+      rt(3,2,i-2)=cost*sinphi
+      rt(3,3,i-2)=cosphi
+        call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
+      do j=1,3
+        dc_norm(j,i-1)=prod(j,1,i-1)
+        dc(j,i-1)=vbld(i)*prod(j,1,i-1)
+      enddo
+        call locate_side_chain(i-1)
+       enddo
+      else
 C
 C Define the origin and orientation of the coordinate system and locate the
 C first three CA's and SC(2).
@@ -59,7 +171,7 @@ C
  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
 
       endif
-
+      endif
       return
       end
 c-------------------------------------------------------------------------
@@ -126,8 +238,8 @@ C
       dc_norm(3,1)=0.0D0
       do j=1,3
         dc_norm(j,2)=prod(j,1,2)
-       dc(j,2)=vbld(3)*prod(j,1,2)
-       c(j,3)=c(j,2)+dc(j,2)
+        dc(j,2)=vbld(3)*prod(j,1,2)
+        c(j,3)=c(j,2)+dc(j,2)
       enddo
       call locate_side_chain(2)
       return

diff --git a/source/unres/src_MD-M/checkder_p.F b/source/unres/src_MD-M/checkder_p.F
index 0539e48..32d2366 100644 (file)
--- a/source/unres/src_MD-M/checkder_p.F
+++ b/source/unres/src_MD-M/checkder_p.F
@@ -272,8 +272,8 @@ C Check the gradient of the energy in Cartesian coordinates.
       integer uiparm(1)
       double precision urparm(1)
       external fdum
-      r_cut=2.0d0
-      rlambd=0.3d0
+c      r_cut=2.0d0
+c      rlambd=0.3d0
       icg=1
       nf=0
       nfl=0                

diff --git a/source/unres/src_MD-M/contact.f b/source/unres/src_MD-M/contact.f
index 24b11d6..cc4e0b7 100644 (file)
--- a/source/unres/src_MD-M/contact.f
+++ b/source/unres/src_MD-M/contact.f
@@ -175,7 +175,7 @@ c      do i=1,nharp
 c            write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
 c      enddo
       if (lprint) then
-      write (iout,*) "Hairpins:"
+      write (iout,*) "Hairpins:",nharp
       do i=1,nharp
         i1=iharp(1,i)
         j1=iharp(2,i)

diff --git a/source/unres/src_MD-M/elecont.f b/source/unres/src_MD-M/elecont.f
index a962630..4faf3e5 100644 (file)
--- a/source/unres/src_MD-M/elecont.f
+++ b/source/unres/src_MD-M/elecont.f
@@ -214,7 +214,7 @@ c--------------------------------------------
       double precision p1,p2
       external freeres
 
-      if(.not.dccart) call chainbuild
+cc????      if(.not.dccart) call chainbuild
 cd      call write_pdb(99,'sec structure',0d0)
       ncont=0
       nbfrag=0

diff --git a/source/unres/src_MD-M/energy_p_new-sep.F b/source/unres/src_MD-M/energy_p_new-sep.F
index 0b8f27b..85f3e1f 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep.F
+++ b/source/unres/src_MD-M/energy_p_new-sep.F
@@ -2,6 +2,7 @@ C-----------------------------------------------------------------------
       double precision function sscale(r)
       double precision r,gamm
       include "COMMON.SPLITELE"
+      include "COMMON.CHAIN"
       if(r.lt.r_cut-rlamb) then
         sscale=1.0d0
       else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
@@ -13,6 +14,23 @@ C-----------------------------------------------------------------------
       return
       end
 C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+      double precision function sscagrad(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+      include "COMMON.CHAIN"
+      if(r.lt.r_cut-rlamb) then
+        sscagrad=0.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscagrad=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscagrad=0.0d0
+      endif
+      return
+      end
+C-----------------------------------------------------------------------
+
       subroutine elj_long(evdw)
 C
 C This subroutine calculates the interaction energy of nonbonded side chains

diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
index 6592ace..574e19b 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
@@ -13,6 +13,21 @@ C-----------------------------------------------------------------------
       return
       end
 C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+      double precision function sscagrad(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+      if(r.lt.r_cut-rlamb) then
+        sscagrad=0.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscagrad=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscagrad=0.0d0
+      endif
+      return
+      end
+C-----------------------------------------------------------------------
       subroutine elj_long(evdw)
 C
 C This subroutine calculates the interaction energy of nonbonded side chains

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 635a41e..8f89ae9 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -24,6 +24,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.MD'
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
+      include 'COMMON.SPLITELE'
 #ifdef MPI      
 c      print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
 c     & " nfgtasks",nfgtasks
@@ -1412,6 +1413,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
       logical lprn
       integer xshift,yshift,zshift
       evdw=0.0D0
@@ -1451,8 +1453,8 @@ C Condition for being inside the proper box
         go to 135
         endif
   136   continue
-        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
-        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
 C Condition for being inside the proper box
         if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
      &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
@@ -1520,8 +1522,8 @@ C Condition for being inside the proper box
         go to 138
         endif
   139   continue
-        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
-        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
 C Condition for being inside the proper box
         if ((zj.gt.((0.5d0)*boxzsize)).or.
      &       (zj.lt.((-0.5d0)*boxzsize))) then
@@ -1539,6 +1541,12 @@ c            write (iout,*) "j",j," dc_norm",
 c     &       dc_norm(1,nres+j),dc_norm(2,nres+j),dc_norm(3,nres+j)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             rij=dsqrt(rrij)
+            sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
+            sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+             
+c            write (iout,'(a7,4f8.3)') 
+c    &      "ssscale",sss,((1.0d0/rij)/sigma(itypi,itypj)),r_cut,rlamb
+            if (sss.gt.0.0d0) then
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -1567,7 +1575,7 @@ c---------------------------------------------------------------
 c            write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
 c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
             evdwij=evdwij*eps2rt*eps3rt
-            evdw=evdw+evdwij
+            evdw=evdw+evdwij*sss
             if (lprn) then
             sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
             epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
@@ -1587,6 +1595,9 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac
+c            print '(2i4,6f8.4)',i,j,sss,sssgrad*
+c     &      evdwij,fac,sigma(itypi,itypj),expon
+            fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
 c            fac=0.0d0
 C Calculate the radial part of the gradient
             gg(1)=xj*fac
@@ -1594,6 +1605,7 @@ C Calculate the radial part of the gradient
             gg(3)=zj*fac
 C Calculate angular part of the gradient.
             call sc_grad
+            endif
           enddo      ! j
         enddo        ! iint
       enddo          ! i
@@ -1807,6 +1819,7 @@ C----------------------------------------------------------------------------
       include 'COMMON.CALC'
       include 'COMMON.IOUNITS'
       double precision dcosom1(3),dcosom2(3)
+cc      print *,'sss=',sss
       eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
       eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
       eom12=evdwij*eps1_om12+eps2der*eps2rt_om12
@@ -1825,16 +1838,16 @@ c      write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
         dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
       enddo
       do k=1,3
-        gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+        gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss
       enddo 
 c      write (iout,*) "gg",(gg(k),k=1,3)
       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
+     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss
         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
+     &            +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss
 c        write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
 c     &            +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
 c        write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
@@ -2770,6 +2783,7 @@ C
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       include 'COMMON.TIME1'
+      include 'COMMON.SPLITELE'
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
@@ -2883,8 +2897,8 @@ C Condition for being inside the proper box
         go to 185
         endif
   186   continue
-        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
-        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
+        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
 C Condition for being inside the proper box
         if ((zmedi.gt.((0.5d0)*boxzsize)).or.
      &       (zmedi.lt.((-0.5d0)*boxzsize))) then
@@ -2926,8 +2940,8 @@ C Condition for being inside the proper box
         go to 195
         endif
   196   continue
-        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
-        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+        if (zmedi.gt.((0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
+        if (zmedi.lt.((-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
 C Condition for being inside the proper box
         if ((zmedi.gt.((0.5d0)*boxzsize)).or.
      &       (zmedi.lt.((-0.5d0)*boxzsize))) then
@@ -2976,8 +2990,8 @@ C Condition for being inside the proper box
         go to 165
         endif
   166   continue
-        if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxxsize
-        if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxxsize
+        if (zmedi.gt.((zshift+0.5d0)*boxzsize)) zmedi=zmedi-boxzsize
+        if (zmedi.lt.((zshift-0.5d0)*boxzsize)) zmedi=zmedi+boxzsize
 C Condition for being inside the proper box
         if ((zmedi.gt.((zshift+0.5d0)*boxzsize)).or.
      &       (zmedi.lt.((zshift-0.5d0)*boxzsize))) then
@@ -3027,6 +3041,7 @@ C-------------------------------------------------------------------------------
       include 'COMMON.VECTORS'
       include 'COMMON.FFIELD'
       include 'COMMON.TIME1'
+      include 'COMMON.SPLITELE'
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
@@ -3085,8 +3100,8 @@ C Condition for being inside the proper box
         go to 175
         endif
   176   continue
-        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
-        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
 C Condition for being inside the proper box
         if ((zj.gt.((0.5d0)*boxzsize)).or.
      &       (zj.lt.((-0.5d0)*boxzsize))) then
@@ -3097,6 +3112,10 @@ C        endif !endPBC condintion
         yj=yj-ymedi
         zj=zj-zmedi
           rij=xj*xj+yj*yj+zj*zj
+
+            sss=sscale(sqrt(rij))
+            sssgrad=sscagrad(sqrt(rij))
+c            if (sss.gt.0.0d0) then  
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
           rmij=1.0D0/rij
@@ -3112,14 +3131,15 @@ c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
           ev2=bbb*r6ij
           fac3=ael6i*r6ij
           fac4=ael3i*r3ij
-          evdwij=ev1+ev2
+          evdwij=(ev1+ev2)
           el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
           el2=fac4*fac       
-          eesij=el1+el2
+C MARYSIA
+          eesij=(el1+el2)
 C 12/26/95 - for the evaluation of multi-body H-bonding interactions
           ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
           ees=ees+eesij
-          evdw1=evdw1+evdwij
+          evdw1=evdw1+evdwij*sss
 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
@@ -3136,7 +3156,7 @@ C
 C Calculate contributions to the Cartesian gradient.
 C
 #ifdef SPLITELE
-          facvdw=-6*rrmij*(ev1+evdwij)
+          facvdw=-6*rrmij*(ev1+evdwij)*sss
           facel=-3*rrmij*(el1+eesij)
           fac1=fac
           erij(1)=xj*rmij
@@ -3188,10 +3208,11 @@ cgrad              gvdwpp(l,k)=gvdwpp(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
 #else
-          facvdw=ev1+evdwij 
-          facel=el1+eesij  
+C MARYSIA
+          facvdw=(ev1+evdwij)*sss
+          facel=(el1+eesij)
           fac1=fac
-          fac=-3*rrmij*(facvdw+facvdw+facel)
+          fac=-3*rrmij*(facvdw+facvdw+facel)+sssgrad*rmij*evdwij
           erij(1)=xj*rmij
           erij(2)=yj*rmij
           erij(3)=zj*rmij
@@ -3220,9 +3241,9 @@ cgrad              gelc(l,k)=gelc(l,k)+ggg(l)
 cgrad            enddo
 cgrad          enddo
 c 9/28/08 AL Gradient compotents will be summed only at the end
-          ggg(1)=facvdw*xj
-          ggg(2)=facvdw*yj
-          ggg(3)=facvdw*zj
+          ggg(1)=facvdw*xj*sss
+          ggg(2)=facvdw*yj*sss
+          ggg(3)=facvdw*zj*sss
           do k=1,3
             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
@@ -3261,14 +3282,16 @@ cgrad            enddo
 cgrad          enddo
           do k=1,3
             gelc(k,i)=gelc(k,i)
-     &               +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
-     &               + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+     &           +(ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
+     &           + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
             gelc(k,j)=gelc(k,j)
-     &               +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
-     &               + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+     &           +(ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
+     &           + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
           enddo
+C MARYSIA
+c          endif !sscale
           IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
      &        .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0 
      &        .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
@@ -3460,9 +3483,9 @@ cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
      &            'eelloc',i,j,eel_loc_ij
-           if (eel_loc_ij.ne.0)
-     &      write (iout,'(a4,2i4,8f9.5)')'chuj',
-     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
+c           if (eel_loc_ij.ne.0)
+c     &      write (iout,'(a4,2i4,8f9.5)')'chuj',
+c     &     i,j,a22,muij(1),a23,muij(2),a32,muij(3),a33,muij(4)
 
           eel_loc=eel_loc+eel_loc_ij
 C Partial derivatives in virtual-bond dihedral angles gamma
@@ -3490,14 +3513,14 @@ cgrad            enddo
 cgrad          enddo
 C Remaining derivatives of eello
           do l=1,3
-            gel_loc(l,i)=gel_loc(l,i)+aggi(l,1)*muij(1)+
-     &          aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)
-            gel_loc(l,i+1)=gel_loc(l,i+1)+aggi1(l,1)*muij(1)+
-     &          aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)
-            gel_loc(l,j)=gel_loc(l,j)+aggj(l,1)*muij(1)+
-     &          aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)
-            gel_loc(l,j1)=gel_loc(l,j1)+aggj1(l,1)*muij(1)+
-     &          aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)
+            gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
+     &        aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
+            gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
+     &     aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
+            gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
+     &       aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
+            gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
+     &     aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
           enddo
           ENDIF
 C Change 12/26/95 to calculate four-body contributions to H-bonding energy
@@ -4050,8 +4073,8 @@ C Condition for being inside the proper box
         go to 135
         endif
   136   continue
-        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
-        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
 C Condition for being inside the proper box
         if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
      &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
@@ -4087,8 +4110,8 @@ C Condition for being inside the proper box
         go to 175
         endif
   176   continue
-        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
-        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
 C Condition for being inside the proper box
         if ((zj.gt.((0.5d0)*boxzsize)).or.
      &       (zj.lt.((-0.5d0)*boxzsize))) then
@@ -4098,6 +4121,7 @@ C Condition for being inside the proper box
           yj=yj-yi
           zj=zj-zi
           rij=xj*xj+yj*yj+zj*zj
+
           r0ij=r0_scp
           r0ijsq=r0ij*r0ij
           if (rij.lt.r0ijsq) then
@@ -4171,6 +4195,7 @@ C
       include 'COMMON.FFIELD'
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
       dimension ggg(3)
       evdw2=0.0D0
       evdw2_14=0.0d0
@@ -4203,8 +4228,8 @@ C Condition for being inside the proper box
         go to 135
         endif
   136   continue
-        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxxsize
-        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxxsize
+        if (zi.gt.((zshift+0.5d0)*boxzsize)) zi=zi-boxzsize
+        if (zi.lt.((zshift-0.5d0)*boxzsize)) zi=zi+boxzsize
 C Condition for being inside the proper box
         if ((zi.gt.((zshift+0.5d0)*boxzsize)).or.
      &       (zi.lt.((zshift-0.5d0)*boxzsize))) then
@@ -4240,8 +4265,8 @@ C Condition for being inside the proper box
         go to 175
         endif
   176   continue
-        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxxsize
-        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxxsize
+        if (zj.gt.((0.5d0)*boxzsize)) zj=zj-boxzsize
+        if (zj.lt.((-0.5d0)*boxzsize)) zj=zj+boxzsize
 C Condition for being inside the proper box
         if ((zj.gt.((0.5d0)*boxzsize)).or.
      &       (zj.lt.((-0.5d0)*boxzsize))) then
@@ -4251,23 +4276,27 @@ C Condition for being inside the proper box
           yj=yj-yi
           zj=zj-zi
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+          sss=sscale(1.0d0/(dsqrt(rrij)))
+          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
+          if (sss.gt.0.0d0) then
           fac=rrij**expon2
           e1=fac*fac*aad(itypj,iteli)
           e2=fac*bad(itypj,iteli)
           if (iabs(j-i) .le. 2) then
             e1=scal14*e1
             e2=scal14*e2
-            evdw2_14=evdw2_14+e1+e2
+            evdw2_14=evdw2_14+(e1+e2)*sss
           endif
           evdwij=e1+e2
-          evdw2=evdw2+evdwij
+          evdw2=evdw2+evdwij*sss
           if (energy_dec) write (iout,'(a6,2i5,0pf7.3,2i3,3e11.3)')
      &        'evdw2',i,j,evdwij,iteli,itypj,fac,aad(itypj,iteli),
      &       bad(itypj,iteli)
 C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
-          fac=-(evdwij+e1)*rrij
+          fac=-(evdwij+e1)*rrij*sss
+          fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
           ggg(1)=xj*fac
           ggg(2)=yj*fac
           ggg(3)=zj*fac
@@ -4302,7 +4331,8 @@ cgrad          enddo
             gvdwc_scpp(k,i)=gvdwc_scpp(k,i)-ggg(k)
             gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
           enddo
-        enddo
+        endif !endif for sscale cutoff
+        enddo ! j
 
         enddo ! iint
       enddo ! i
@@ -4622,7 +4652,9 @@ c      time12=1.0d0
       etheta=0.0D0
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
       do i=ithet_start,ithet_end
+        print *,i,itype(i-1),itype(i),itype(i-2)
         if (itype(i-1).eq.ntyp1) cycle
+        print *,'wchodze',itype(i-1)
 C Zero the energy function and its derivative at 0 or pi.
         call splinthet(theta(i),0.5d0*delta,ss,ssd)
         it=itype(i-1)
@@ -4720,7 +4752,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2.
      &      'ebend',i,ethetai,theta(i),itype(i)
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
-        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
+        gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)+gloc(nphi+i-2,icg)
       enddo
 C Ufff.... We've done all this!!! 
       return
@@ -4739,7 +4771,7 @@ C Calculate the contributions to both Gaussian lobes.
 C 6/6/97 - Deform the Gaussians using the factor of 1/(1+time)
 C The "polynomial part" of the "standard deviation" of this part of 
 C the distributioni.
-        write (iout,*) thetai,thet_pred_mean
+ccc        write (iout,*) thetai,thet_pred_mean
         sig=polthet(3,it)
         do j=2,0,-1
           sig=sig*thet_pred_mean+polthet(j,it)
@@ -4864,6 +4896,7 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
+c        print *,i,itype(i-1),itype(i),itype(i-2)
         if ((itype(i-1).eq.ntyp1)) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
@@ -5029,7 +5062,7 @@ c        lprn1=.false.
         etheta=etheta+ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
-        gloc(nphi+i-2,icg)=wang*dethetai
+        gloc(nphi+i-2,icg)=wang*dethetai+ gloc(nphi+i-2,icg)
       enddo
       return
       end
@@ -5892,9 +5925,11 @@ c     lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
 C ANY TWO ARE DUMMY ATOMS in row CYCLE
-        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
-     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1))  .or.
-     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
+c        if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or.
+c     &      ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1))  .or.
+c     &      ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle
+         if ((itype(i-3).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &  (itype(i-1).eq.ntyp1).or.(itype(i).eq.ntyp1)) cycle
 C For introducing the NH3+ and COO- group please check the etor_d for reference
 C and guidance
         etors_ii=0.0D0

diff --git a/source/unres/src_MD-M/geomout.F b/source/unres/src_MD-M/geomout.F
index c6a3944..4673c4e 100644 (file)
--- a/source/unres/src_MD-M/geomout.F
+++ b/source/unres/src_MD-M/geomout.F
@@ -91,7 +91,7 @@ cmodel      write (iunit,'(a5,i6)') 'MODEL',1
       ires=0
       do i=nnt,nct
         iti=itype(i)
-        if (iti.eq.ntyp1) then
+        if ((iti.eq.ntyp1).and.((itype(i+1)).eq.ntyp1)) then
           ichain=ichain+1
           ires=0
           write (iunit,'(a)') 'TER'
@@ -99,6 +99,7 @@ cmodel      write (iunit,'(a5,i6)') 'MODEL',1
         ires=ires+1
         iatom=iatom+1
         ica(i)=iatom
+        if (iti.ne.ntyp1) then
         write (iunit,10) iatom,restyp(iti),chainid(ichain),
      &     ires,(c(j,i),j=1,3),vtot(i)
         if (iti.ne.10) then
@@ -106,6 +107,7 @@ cmodel      write (iunit,'(a5,i6)') 'MODEL',1
           write (iunit,20) iatom,restyp(iti),chainid(ichain),
      &      ires,(c(j,nres+i),j=1,3),
      &      vtot(i+nres)
+         endif
         endif
         endif
       enddo

diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F
index a650241..70e35ab 100644 (file)
--- a/source/unres/src_MD-M/initialize_p.F
+++ b/source/unres/src_MD-M/initialize_p.F
@@ -267,8 +267,8 @@ C
 C Initialize constants used to split the energy into long- and short-range
 C components
 C
-      r_cut=2.0d0
-      rlamb=0.3d0
+C      r_cut=2.0d0
+C      rlamb=0.3d0
 #ifndef SPLITELE
       nprint_ene=nprint_ene-1
 #endif

diff --git a/source/unres/src_MD-M/intcartderiv.F b/source/unres/src_MD-M/intcartderiv.F
index 369a4f0..d715c19 100644 (file)
--- a/source/unres/src_MD-M/intcartderiv.F
+++ b/source/unres/src_MD-M/intcartderiv.F
@@ -692,7 +692,7 @@ c     Check omega gradient
       enddo
       return
       end
-
+c------------------------------------------------------------
       subroutine chainbuild_cart
       implicit real*8 (a-h,o-z)
       include 'DIMENSIONS'
@@ -736,6 +736,7 @@ c        call flush(iout)
 #endif
       do j=1,3
         c(j,1)=dc(j,0)
+c        c(j,1)=c(j,1)
       enddo
       do i=2,nres
         do j=1,3

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index bd2165b..22236bc 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -591,6 +591,7 @@ C
       read (itorp,*,end=113,err=113) ntortyp,nterm_old
       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
+      itortyp(ntyp1)=0
       do i=1,ntortyp
        do j=1,ntortyp
          read (itorp,'(a)')

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index de4b92e..8e7b152 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -8,6 +8,7 @@
       include 'COMMON.CONTROL'
       include 'COMMON.SBRIDGE'
       include 'COMMON.IOUNITS'
+      include 'COMMON.SPLITELE'
       logical file_exist
 C Read force-field parameters except weights
       call parmread
@@ -79,6 +80,7 @@ C
       include 'COMMON.FFIELD'
       include 'COMMON.INTERACT'
       include 'COMMON.SETUP'
+      include 'COMMON.SPLITELE'
       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
       character*80 ucase
@@ -219,7 +221,9 @@ C Reading the dimensions of box in x,y,z coordinates
       call reada(controlcard,'BOXX',boxxsize,100.0d0)
       call reada(controlcard,'BOXY',boxysize,100.0d0)
       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
-
+c Cutoff range for interactions
+      call reada(controlcard,"R_CUT",r_cut,15.0d0)
+      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
       if (me.eq.king .or. .not.out1file ) 
      &  write (iout,*) "DISTCHAINMAX",distchainmax
       
@@ -346,8 +350,8 @@ C
       ntime_split0=ntime_split
       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
       ntime_split0=ntime_split
-      call reada(controlcard,"R_CUT",r_cut,2.0d0)
-      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
+c      call reada(controlcard,"R_CUT",r_cut,2.0d0)
+c      call reada(controlcard,"LAMBDA",rlamb,0.3d0)
       rest = index(controlcard,"REST").gt.0
       tbf = index(controlcard,"TBF").gt.0
       usampl = index(controlcard,"USAMPL").gt.0

diff --git a/source/unres/src_MD-M/refsys.f b/source/unres/src_MD-M/refsys.f
index 86b0524..dff4aa7 100644 (file)
--- a/source/unres/src_MD-M/refsys.f
+++ b/source/unres/src_MD-M/refsys.f
@@ -14,12 +14,16 @@ c form a linear fragment. Returns vectors e1, e2, and e3.
       include 'COMMON.IOUNITS'
       include 'COMMON.CHAIN'
       double precision coinc/1.0D-13/,align /1.0D-13/
+c      print *,'just initialize'
       fail=.false.
+c      print *,fail
       s1=0.0
       s2=0.0
+      print *,s1,s2
       do 1 i=1,3
       zi=c(i,i2)-c(i,i3)
       ui=c(i,i4)-c(i,i3)
+      print *,zi,ui
       s1=s1+zi*zi
       s2=s2+ui*ui
       z(i)=zi
@@ -38,6 +42,7 @@ c   3 c(i,i1)=0.0D0
       do 5 i=1,3
     5 c(i,i1)=0.0D0
       return
+      print *,'two if pass'
     4 s1=1.0/s1
       s2=1.0/s2
       v1=z(2)*u(3)-z(3)*u(2)
@@ -60,6 +65,7 @@ c   7 c(i,i1)=0.0D0
       e2(1)=e1(3)*e3(2)-e1(2)*e3(3)
       e2(2)=e1(1)*e3(3)-e1(3)*e3(1)
       e2(3)=e1(2)*e3(1)-e1(1)*e3(2)
+      print *,'just before leave'
  1000 format (/1x,' * * * Error - atoms',i4,' and',i4,' coincide.',
      1 'coordinates of atom',i4,' are set to zero.')
  1010 format (/1x,' * * * Error - atoms',2(i4,2h, ),i4,' form a linear',

diff --git a/source/unres/src_MD-M/unres.F b/source/unres/src_MD-M/unres.F
index 0039fcc..83afd9f 100644 (file)
--- a/source/unres/src_MD-M/unres.F
+++ b/source/unres/src_MD-M/unres.F
@@ -189,9 +189,16 @@ c---------------------------------------------------------------------------
       double precision energy(0:n_ene)
       double precision energy_long(0:n_ene),energy_short(0:n_ene)
       double precision varia(maxvar)
-      if (indpdb.eq.0) call chainbuild
+      if (indpdb.eq.0)     call chainbuild
       time00=MPI_Wtime()
+      print *,'dc',c(1,1)
+      if (indpdb.ne.0) then
+      dc(1,0)=c(1,1)
+      dc(2,0)=c(2,1)
+      dc(3,0)=c(3,1)
+      endif
       call chainbuild_cart
+      print *,'dc',dc(1,0),dc(2,0),dc(3,0)
       if (split_ene) then
        print *,"Processor",myrank," after chainbuild"
        icall=1
@@ -216,7 +223,9 @@ c---------------------------------------------------------------------------
       etot =etota
       call enerprint(energy(0))
       call hairpin(.true.,nharp,iharp)
+        print *,'after hairpin'
       call secondary2(.true.)
+        print *,'after secondary'
       if (minim) then
 crc overlap test
         if (overlapsc) then 
@@ -248,8 +257,10 @@ crc overlap test
         evals=nfun/(MPI_WTIME()-time1)
         print *,'# eval/s',evals
         print *,'refstr=',refstr
-        call hairpin(.true.,nharp,iharp)
+        call hairpin(.false.,nharp,iharp)
+        print *,'after hairpin'
         call secondary2(.true.)
+        print *,'after secondary'
         call etotal(energy(0))
         etot = energy(0)
         call enerprint(energy(0))