sadas

diff --git a/source/unres/src_MD/COMMON.SCCOR b/source/unres/src_MD/COMMON.SCCOR
index 5217de7..1c725d7 100644 (file)
--- a/source/unres/src_MD/COMMON.SCCOR
+++ b/source/unres/src_MD/COMMON.SCCOR
@@ -1,6 +1,11 @@
-C Parameters of the SCCOR term
-      double precision v1sccor,v2sccor
-      integer nterm_sccor
+ Parameters of the SCCOR term
+      double precision v1sccor,v2sccor,vlor1sccor,
+     &                 vlor2sccor,vlor3sccor
+      integer nterm_sccor,isccortyp,nsccortyp,nlor_sccor
       common/torsion/v1sccor(maxterm_sccor,20,20),
      &    v2sccor(maxterm_sccor,20,20),
-     &    nterm_sccor
+     &    nterm_sccor(ntyp,ntyp),isccortyp(ntyp),nsccortyp,
+     &    nlor_sccor(ntyp,ntyp),vlor1sccor(maxterm_sccor,20,20),
+     &    vlor2sccor(maxterm_sccor,20,20),
+     &    vlor3sccor(maxterm_sccor,20,20)
+

diff --git a/source/unres/src_MD/COMMON.VAR b/source/unres/src_MD/COMMON.VAR
index 71158b8..c2a0a25 100644 (file)
--- a/source/unres/src_MD/COMMON.VAR
+++ b/source/unres/src_MD/COMMON.VAR
@@ -5,6 +5,7 @@ C Store the geometric variables in the following COMMON block.
      &          thetaref,phiref,costtab,sinttab,cost2tab,sint2tab,
      &          xxtab,yytab,zztab,xxref,yyref,zzref
       common /var/ theta(maxres),phi(maxres),alph(maxres),omeg(maxres),
+     &          omicron(2,maxres),tauangle(3,maxres)
      &          vbld(2*maxres),thetaref(maxres),phiref(maxres),
      &          costtab(maxres), sinttab(maxres), cost2tab(maxres),
      &          sint2tab(maxres),xxtab(maxres),yytab(maxres),

diff --git a/source/unres/src_MD/checkder_p.F b/source/unres/src_MD/checkder_p.F
index 26854e6..ec66acd 100644 (file)
--- a/source/unres/src_MD/checkder_p.F
+++ b/source/unres/src_MD/checkder_p.F
@@ -513,7 +513,20 @@ c-------------------------------------------------------------------------
      &     +(c(j,i+1)-c(j,i))/dnorm2)
         enddo
         be=0.0D0
-        if (i.gt.2) phi(i+1)=beta(i-2,i-1,i,i+1)
+        if (i.gt.2) then
+        phi(i+1)=beta(i-2,i-1,i,i+1)
+        if (itype(i).ne.10).and.(itype(i-1).ne.10) then
+         tauangle(3,i+1)=beta(i+nres-1,i-1,i,i+nres)
+        endif
+        if (itype(i-1).ne.10) then
+         tauangle(1,i+1)=beta(i-1+nres,i-1,i,i+1)
+         omicron(1,i)=alpha(i-2,i-1,i-1+nres)
+         omicron(2,i)=alpha(i-1+nres,i-1,i)
+        endif
+        if (itype(i).ne.10) then
+         tauangle(2,i+1)=beta(i-2,i-1,i,i+nres)
+        endif
+        endif
         omeg(i)=beta(nres+i,i,maxres2,i+1)
         alph(i)=alpha(nres+i,i,maxres2)
         theta(i+1)=alpha(i-1,i,i+1)

diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F
index e3f2c36..23449da 100644 (file)
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -5801,7 +5801,6 @@ c     lprn=.true.
         phii1=phi(i+1)
         gloci1=0.0D0
         gloci2=0.0D0
-C Regular cosine and sine terms
         do j=1,ntermd_1(itori,itori1,itori2)
           v1cij=v1c(1,j,itori,itori1,itori2)
           v1sij=v1s(1,j,itori,itori1,itori2)
@@ -5870,24 +5869,43 @@ c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
       do i=iphi_start,iphi_end
         esccor_ii=0.0D0
-        itori=itype(i-2)
-        itori1=itype(i-1)
+        isccori=isccortyp(itype(i-2))
+        isccori1=isccortyp(itype(i-1))
         phii=phi(i)
+cccc  Added 9 May 2012
+cc Tauangle is torsional engle depending on the value of first digit 
+c(see comment below)
+cc Omicron is flat angle depending on the value of first digit 
+c(see comment below)
+
         gloci=0.0D0
-        do j=1,nterm_sccor
-          v1ij=v1sccor(j,itori,itori1)
-          v2ij=v2sccor(j,itori,itori1)
-          cosphi=dcos(j*phii)
-          sinphi=dsin(j*phii)
+        do intertyp=1,3
+cc Added 09 May 2012 (Adasko)
+cc  Intertyp means interaction type of backbone mainchain correlation: 
+c   1 = SC...Ca...Ca...Ca
+c   2 = Ca...Ca...Ca...SC
+c   3 = SC...Ca...Ca...SC
+        if (((intertyp.eq.3).and.(itype(i-2).eq.10).or.
+     &      (itype(i-1).eq.10))
+     &    .or. ((intertyp.eq.1).and.(itype(i-2).ne.10))
+     &    .or. ((intertyp.eq.2).and.(itype(i-1).ne.10))) cycle  
+        do j=1,nterm_sccor(isccori,isccori1)
+          v1ij=v1sccor(j,intertyp,isccori,isccori1)
+          v2ij=v2sccor(j,intertyp,isccori,isccori1)
+          cosphi=dcos(j*tauangle(intertyp,i))
+          sinphi=dsin(j*tauangle(intertyp,i))
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
+        gloc_sc(intertyp,i-3,icg)=gloc_sc(i-3,icg)+wtor*gloci
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
-     &  (v1sccor(j,itori,itori1),j=1,6),(v2sccor(j,itori,itori1),j=1,6)
+     &  (v1sccor(j,intertyp,itori,itori1),j=1,6)
+     & ,(v2sccor(j,intertyp,itori,itori1),j=1,6)
         gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
       enddo
+      enddo
       return
       end
 c----------------------------------------------------------------------------

diff --git a/source/unres/src_MD/int_to_cart.f b/source/unres/src_MD/int_to_cart.f
index 0528af7..2564372 100644 (file)
--- a/source/unres/src_MD/int_to_cart.f
+++ b/source/unres/src_MD/int_to_cart.f
@@ -113,6 +113,108 @@ c   The side-chain vector derivatives
             enddo
          endif     
        enddo                                                                                                                                                   
+C INTERTYP=1 SC...Ca...Ca...Ca
+c   calculating dE/ddc1      
+       intertyp=1
+       if (nres.lt.2) return
+       if ((nres.lt.3).and.(itype(1).eq.1)) return
+       if (itype(1).ne.10) then
+        do j=1,3
+         gxcart(j,1)=gxcart(j,1)+gloc_sc(intertyp,1,icg)*
+     &     dtauangle(j,1,1,4)
+c As potetnial DO NOT 
+c     &     +gloc_sc(intertyp,nres-2,icg)*dtheta(j,1,3)         
+        enddo
+       endif
+c     Calculating the remainder of dE/ddc2
+       do j=1,3
+         gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+
+     &  gloc(nres-2,icg)*dtheta(j,2,3)+gloc(nres-1,icg)*dtheta(j,1,4)
+        if(itype(2).ne.10) then
+          gcart(j,2)=gcart(j,2)+gloc(ialph(2,1),icg)*dalpha(j,2,2)+
+     &    gloc(ialph(2,1)+nside,icg)*domega(j,2,2)
+        endif
+               if(itype(3).ne.10) then
+         gcart(j,2)=gcart(j,2)+gloc(ialph(3,1),icg)*dalpha(j,1,3)+
+     &    gloc(ialph(3,1)+nside,icg)*domega(j,1,3)
+        endif
+        if(nres.gt.4) then
+          gcart(j,2)=gcart(j,2)+gloc(2,icg)*dphi(j,1,5)
+        endif                  
+       enddo
+c  If there are only five residues       
+       if(nres.eq.5) then
+         do j=1,3
+           gcart(j,3)=gcart(j,3)+gloc(1,icg)*dphi(j,3,4)+gloc(2,icg)*
+     &     dphi(j,2,5)+gloc(nres-1,icg)*dtheta(j,2,4)+gloc(nres,icg)*
+     &     dtheta(j,1,5)
+         if(itype(3).ne.10) then
+          gcart(j,3)=gcart(j,3)+gloc(ialph(3,1),icg)*
+     &    dalpha(j,2,3)+gloc(ialph(3,1)+nside,icg)*domega(j,2,3)
+         endif
+        if(itype(4).ne.10) then
+          gcart(j,3)=gcart(j,3)+gloc(ialph(4,1),icg)*
+     &    dalpha(j,1,4)+gloc(ialph(4,1)+nside,icg)*domega(j,1,4)
+         endif
+       enddo
+       endif
+c    If there are more than five residues
+      if(nres.gt.5) then                          
+        do i=3,nres-3
+         do j=1,3
+          gcart(j,i)=gcart(j,i)+gloc(i-2,icg)*dphi(j,3,i+1)
+     &    +gloc(i-1,icg)*dphi(j,2,i+2)+
+     &    gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+
+     &    gloc(nres+i-3,icg)*dtheta(j,1,i+2)
+          if(itype(i).ne.10) then
+           gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+
+     &     gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
+          endif
+          if(itype(i+1).ne.10) then
+           gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1)
+     &     +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
+          endif
+         enddo
+        enddo
+      endif    
+c  Setting dE/ddnres-2       
+      if(nres.gt.5) then
+         do j=1,3
+           gcart(j,nres-2)=gcart(j,nres-2)+gloc(nres-4,icg)*
+     &    dphi(j,3,nres-1)+gloc(nres-3,icg)*dphi(j,2,nres)
+     &     +gloc(2*nres-6,icg)*
+     &     dtheta(j,2,nres-1)+gloc(2*nres-5,icg)*dtheta(j,1,nres)
+          if(itype(nres-2).ne.10) then
+              gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-2,1),icg)*
+     &       dalpha(j,2,nres-2)+gloc(ialph(nres-2,1)+nside,icg)*
+     &        domega(j,2,nres-2)
+          endif
+          if(itype(nres-1).ne.10) then
+             gcart(j,nres-2)=gcart(j,nres-2)+gloc(ialph(nres-1,1),icg)*
+     &      dalpha(j,1,nres-1)+gloc(ialph(nres-1,1)+nside,icg)*
+     &       domega(j,1,nres-1)
+          endif
+         enddo
+      endif 
+c  Settind dE/ddnres-1       
+       do j=1,3
+        gcart(j,nres-1)=gcart(j,nres-1)+gloc(nres-3,icg)*dphi(j,3,nres)+
+     & gloc(2*nres-5,icg)*dtheta(j,2,nres)
+        if(itype(nres-1).ne.10) then
+          gcart(j,nres-1)=gcart(j,nres-1)+gloc(ialph(nres-1,1),icg)*
+     &   dalpha(j,2,nres-1)+gloc(ialph(nres-1,1)+nside,icg)*
+     &    domega(j,2,nres-1)
+        endif
+        enddo
+c   The side-chain vector derivatives
+        do i=2,nres-1
+         if(itype(i).ne.10) then       
+            do j=1,3   
+              gxcart(j,i)=gxcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,3,i)
+     &        +gloc(ialph(i,1)+nside,icg)*domega(j,3,i)
+            enddo
+         endif     
+       enddo                                                                                                                                                   
       return
       end      
        

diff --git a/source/unres/src_MD/intcartderiv.F b/source/unres/src_MD/intcartderiv.F
index 5fea875..c7cf95c 100644 (file)
--- a/source/unres/src_MD/intcartderiv.F
+++ b/source/unres/src_MD/intcartderiv.F
@@ -44,6 +44,43 @@ c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
           dtheta(j,2,i)=-1/sint*dcostheta(j,2,i)     
         enddo
       enddo
+
+#if defined(MPI) && defined(PARINTDER)
+c We need dtheta(:,:,i-1) to compute dphi(:,:,i)
+      do i=max0(ithet_start-1,3),ithet_end
+#else
+      do i=3,nres
+#endif
+      if (itype(i-1).ne.10) then
+        cost1=dcos(omicron(1,i))
+       sint1=sqrt(1-cost1*cost1)
+        cost2=dcos(omicron(2,i))
+        sint2=sqrt(1-cost2*cost2)
+        do j=1,3
+CC Calculate derivative over first omicron (Cai-2,Cai-1,SCi-1) 
+          dcosomicron(j,1,1,i)=-(-dc_norm(j,i-1+nres)+
+     &    cost1*dc_norm(j,i-2))/
+     &   vbld(i-1)
+          domicron(j,1,1,i)=-1/sint1*dcosomicron(j,1,1,i)     
+          dcosomicron(j,1,2,i)=-(dc_norm(j,i-2)
+     &    +cost1*(-dc_norm(j,i-1+nres)))/
+     &   vbld(i+nres)
+          domicron(j,1,2,i)=-1/sint1*dcosomicron(j,1,2,i)  
+CC Calculate derivative over second omicron Sci-1,Cai-1 Cai
+CC Looks messy but better than if in loop
+          dcosomicron(j,2,1,i)=-(-dc_norm(j,i-1+nres)
+     &    +cost2*dc_norm(j,i-1))/
+     &    vbld(i)
+          domicron(j,2,1,i)=-1/sint2*dcosomicron(j,2,1,i)
+          dcosomicron(j,2,2,i)=-(dc_norm(j,i-1)
+     &     +cost1*(-dc_norm(j,i-1+nres)))/
+     &    vbld(i+nres)
+          domicron(j,2,2,i)=-1/sint2*dcosomicron(j,2,2,i)   
+        enddo
+       endif
+      enddo
+
+
       
 c Derivatives of phi:
 c If phi is 0 or 180 degrees, then the formulas 
@@ -109,6 +146,210 @@ c   Obtaining the gamma derivatives from cosine derivative
          enddo
         endif                                                                                           
       enddo
+
+Calculate derivative of Tauangle
+#ifdef PARINTDER
+      do i=iphi1_start,iphi1_end
+#else
+      do i=4,nres
+#endif
+cc INTERTYP=1 SC...Ca...Ca..Ca
+c the conventional case
+        sint=dsin(theta(i))
+        sint1=dsin(omicron(2,i-1))
+        sing=dsin(tauangle(1,i))
+        cost=dcos(theta(i))
+        cost1=dcos(omicron(2,i-1))
+        cosg=dcos(tauangle(1,i))
+        do j=1,3
+        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+        enddo
+        scalp=scalar(dc_norm2(1,i-2+nres),dc_norm(1,i-1))
+        fac0=1.0d0/(sint1*sint)
+        fac1=cost*fac0
+        fac2=cost1*fac0
+        fac3=cosg*cost1/(sint1*sint1)
+        fac4=cosg*cost/(sint*sint)
+c    Obtaining the gamma derivatives from sine derivative                                
+       if (tauangle(1,i).gt.-pi4.and.tauangle(1,i).le.pi4.or.
+     &     tauangle(1,i).gt.pi34.and.tauangle(1,i).le.pi.or.
+     &     tauangle(1,i).gt.-pi.and.tauangle(1,i).le.-pi34) then
+         do j=1,3
+         dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+         enddo
+         call vecpr(dc_norm(1,i-1),dc_norm(1,i-2),vp1)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-1),vp2)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
+        do j=1,3
+            ctgt=cost/sint
+            ctgt1=cost1/sint1
+            cosg_inv=1.0d0/cosg
+            dsintau(j,1,1,i)=-sing*ctgt1*domicron(j,2,1,i-1)
+     &-(fac0*vp1(j)+sing*(-dc_norm(j,i-2+nres)))
+     & *vbld_inv(i-2+nres)
+            dtauangle(j,1,1,i)=cosg_inv*dsintau(j,1,1,i)
+            dsintau(j,1,2,i)=
+     &        -sing*(ctgt1*domicron(j,2,2,i-1)+ctgt*dtheta(j,1,i))
+     &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+            dtauangle(j,1,2,i)=cosg_inv*dsintau(j,1,2,i)
+c Bug fixed 3/24/05 (AL)
+            dsintau(j,1,3,i)=-sing*ctgt*dtheta(j,2,i)
+     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i)
+c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+            dtauangle(j,1,3,i)=cosg_inv*dsinphi(j,1,3,i)
+         enddo                                          
+c   Obtaining the gamma derivatives from cosine derivative
+        else
+           do j=1,3
+           dcostau(j,1,1,i)=fac1*dcosomicron(j,2,1,i-1)+fac3*
+     &     dcosomicron(j,2,1,i-1)-fac0*(dc_norm(j,i-1)-scalp*
+     &     (-dc_norm(j,i-2+nres)))/vbld(i-2+nres)
+           dtauangle(j,1,1,i)=-1/sing*dcostau(j,1,1,i)
+           dcostau(j,1,2,i)=fac1*dcosomicron(j,2,2,i-1)+fac2*
+     &     dcostheta(j,1,i)+fac3*dcosomicron(j,2,2,i-1)+fac4*
+     &     dcostheta(j,1,i)
+           dtauangle(j,1,2,i)=-1/sing*dcostau(j,1,2,i)
+           dcostau(j,1,3,i)=fac2*dcostheta(j,2,i)+fac4*
+     &     dcostheta(j,2,i)-fac0*(-dc_norm(j,i-2+nres)-scalp*
+     &     dc_norm(j,i-1))/vbld(i)
+           dtauangle(j,1,3,i)=-1/sing*dcostau(j,1,3,i)
+         enddo
+        endif                                                                                            
+      enddo
+CC Second case Ca...Ca...Ca...SC
+#ifdef PARINTDER
+      do i=iphi1_start,iphi1_end
+#else
+      do i=4,nres
+#endif
+c the conventional case
+        sint=dsin(omicron(1,i))
+        sint1=dsin(theta(i-1))
+        sing=dsin(tauangle(2,i)
+        cost=dcos(omicron(1,i)
+        cost1=dcos(theta(i-1))
+        cosg=dcos(tauangle(2,i))
+        do j=1,3
+        dc_norm2(j,i-1+nres)=-dc_norm(j,i-1+nres)
+        enddo
+        scalp=scalar(dc_norm(1,i-3),dc_norm2(1,i-1+nres))
+        fac0=1.0d0/(sint1*sint)
+        fac1=cost*fac0
+        fac2=cost1*fac0
+        fac3=cosg*cost1/(sint1*sint1)
+        fac4=cosg*cost/(sint*sint)
+c    Obtaining the gamma derivatives from sine derivative                                
+       if (tauangle(2,i).gt.-pi4.and.tauangle(2,i).le.pi4.or.
+     &     tauangle(2,i).gt.pi34.and.tauangle(2,i).le.pi.or.
+     &     tauangle(2,i).gt.-pi.and.tauangle(2,i).le.-pi34) then
+         call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
+         call vecpr(dc_norm(1,i-3),dc_norm2(1,i-1+nres),vp2)
+         call vecpr(dc_norm(1,i-3),dc_norm(1,i-2),vp3)
+        do j=1,3
+            ctgt=cost/sint
+            ctgt1=cost1/sint1
+            cosg_inv=1.0d0/cosg
+            dsintau(j,2,1,i)=-sing*ctgt1*dtheta(j,1,i-1)
+     &        -(fac0*vp1(j)+sing*dc_norm(j,i-3))*vbld_inv(i-2)
+            dtauangle(j,2,1,i)=cosg_inv*dsintau(j,2,1,i)
+            dsintau(j,2,2,i)=
+     &        -sing*(ctgt1*dtheta(j,2,i-1)+ctgt*domicron(j,1,1,i))
+     &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1)
+            dphi(j,2,i)=cosg_inv*dsintau(j,2,2,i)
+c Bug fixed 3/24/05 (AL)
+            dsintau(j,2,3,i)=-sing*ctgt*domicron(j,1,2,i)
+     &        +(fac0*vp3(j)-sing*dc_norm2(j,i-1+nres))*vbld_inv(i)
+c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+            dtauangle(j,2,3,i)=cosg_inv*dsinphi(j,2,3,i)
+         enddo                                          
+c   Obtaining the gamma derivatives from cosine derivative
+        else
+           do j=1,3
+           dcostau(j,2,1,i)=fac1*dcostheta(j,1,i-1)+fac3*
+     &     dcostheta(j,1,i-1)-fac0*(dc_norm2(j,i-1+nres)-scalp*
+     &     dc_norm(j,i-3))/vbld(i-2)
+           dtauangle(j,2,1,i)=-1/sing*dcostau(j,2,1,i)
+           dcostau(j,2,2,i)=fac1*dcostheta(j,2,i-1)+fac2*
+     &     dcosomicron(j,1,1,i)+fac3*dcostheta(j,2,i-1)+fac4*
+     &     dcosomicron(j,1,1,i)
+           dtauanlge(j,2,2,i)=-1/sing*dcostau(j,2,2,i)
+           dcostau(j,2,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
+     &     dcostheta(j,1,2,i)-fac0*(dc_norm(j,i-3)-scalp*
+     &     dc_norm2(j,i-1+nres))/vbld(i-1+nres)
+           dtauanlge(j,2,3,i)=-1/sing*dcosphi(j,3,i)
+         enddo
+        endif                                                                                            
+      enddo
+
+
+CCC third case SC...Ca...Ca...SC
+#ifdef PARINTDER
+
+      do i=iphi1_start,iphi1_end
+#else
+      do i=4,nres
+#endif
+c the conventional case
+        sint=dsin(omicron(1,i))
+        sint1=dsin(omicron(2,i-1))
+        sing=dsin(tauangle(3,i))
+        cost=dcos(omicron(1,i))
+        cost1=dcos(omicron(2,i-1))
+        cosg=dcos(tauangle(3,i))
+        do j=1,3
+        dc_norm2(j,i-2+nres)=-dc_norm(j,i-2+nres)
+        dc_norm2(j,i-1+nres)=-dc_norm(j,i-2+nres)
+        enddo
+        scalp=scalar(dc_norm2(1,i-2+nres),dc_norm2(1,i-1+nres))
+        fac0=1.0d0/(sint1*sint)
+        fac1=cost*fac0
+        fac2=cost1*fac0
+        fac3=cosg*cost1/(sint1*sint1)
+        fac4=cosg*cost/(sint*sint)
+c    Obtaining the gamma derivatives from sine derivative                                
+       if (tauangle(3,i).gt.-pi4.and.tauangle(3,i).le.pi4.or.
+     &     tauangle(3,i).gt.pi34.and.tauangle(3,i).le.pi.or.
+     &     tauangle(3,i).gt.-pi.and.tauangle(3,i).le.-pi34) then
+         call vecpr(dc_norm2(1,i-1+nres),dc_norm(1,i-2),vp1)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm2(1,i-1+nres),vp2)
+         call vecpr(dc_norm2(1,i-2+nres),dc_norm(1,i-2),vp3)
+        do j=1,3
+            ctgt=cost/sint
+            ctgt1=cost1/sint1
+            cosg_inv=1.0d0/cosg
+            dsintau(j,3,1,i)=-sing*ctgt1*domicron(j,2,1,i-1)
+     &        -(fac0*vp1(j)+sing*dc_norm2(j,i-2+nres))
+     &        *vbld_inv(i-2+nres)
+            dtauangle(j,3,1,i)=cosg_inv*dsintau(j,3,1,i)
+            dsintau(j,3,2,i)=
+     &        -sing*(ctgt1*domicron(j,2,2,i-1)+ctgt*domicron(j,1,1,i))
+     &        -(fac0*vp2(j)+sing*dc_norm(j,i-2))*vbld_inv(i-1+nres)
+            dtauangle(j,3,2,i)=cosg_inv*dsintau(j,3,2,i)
+c Bug fixed 3/24/05 (AL)
+            dsintau(j,3,3,i)=-sing*ctgt*domicron(j,1,2,i)
+     &        +(fac0*vp3(j)-sing*dc_norm2(j,i-1+nres))
+     &        *vbld_inv(i-1+nres)
+c     &        +(fac0*vp3(j)-sing*dc_norm(j,i-1))*vbld_inv(i-1)
+            dphi(j,3,3,i)=cosg_inv*dsintau(j,3,3,i)
+         enddo                                          
+c   Obtaining the gamma derivatives from cosine derivative
+        else
+           do j=1,3
+           dcostau(j,3,1,i)=fac1*dcosomicron(j,2,1,i-1)+fac3*
+     &     dcostheta(j,1,i-1)-fac0*(dc_norm2(j,i-1+nres)-scalp*
+     &     dc_norm2(j,i-2+nres))/vbld(i-2+nres)
+           dtauangle(j,3,1,i)=-1/sing*dcostau(j,3,1,i)
+           dcostau(j,3,2,i)=fac1*dcosomicron(j,2,2,i-1)+fac2*
+     &     dcosomicron(j,1,1,i)+fac3*dcosomicron(j,2,2,i-1)+fac4*
+     &     dcosomicron(j,1,1,i)
+           dtauangle(j,3,2,i)=-1/sing*dcostau(j,3,2,i)
+           dcostau(j,3,3,i)=fac2*dcosomicron(j,1,2,i)+fac4*
+     &     dcostau(j,3,2,i)-fac0*(dc_norm2(j,i-2+nres)-scalp*
+     &     dc_norm2(j,i-1+nres))/vbld(i-1+nres)
+           dtauangle(j,3,3,i)=-1/sing*dcostau(j,3,3,i)
+         enddo
+        endif                                                                                            
+      enddo
 #ifdef CRYST_SC
 c   Derivatives of side-chain angles alpha and omega
 #if defined(MPI) && defined(PARINTDER)

diff --git a/source/unres/src_MD/parmread.F b/source/unres/src_MD/parmread.F
index be5f8b8..be1a193 100644 (file)
--- a/source/unres/src_MD/parmread.F
+++ b/source/unres/src_MD/parmread.F
@@ -553,33 +553,57 @@ C
       enddo
       endif
 #endif
+C Read of Side-chain backbone correlation parameters
+C Modified 11 May 2012 by Adasko
+CCC
 C
-C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
-C         correlation energies.
-C
-      read (isccor,*,end=119,err=119) nterm_sccor
-      do i=1,20
-       do j=1,20
-          read (isccor,'(a)')
-         do k=1,nterm_sccor
-           read (isccor,*,end=119,err=119) kk,v1sccor(k,i,j),
-     &        v2sccor(k,i,j) 
+      read (isccor,*,end=113,err=113) nsccortyp
+      read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
+c      write (iout,*) 'ntortyp',ntortyp
+      maxinter=3
+cc maxinter is maximum interaction sites
+      do l=1,maxinter
+      do i=1,nsccortyp
+       do j=1,nsccortyp
+         read (isccor,*,end=113,err=113) nterm_sccor(i,j),nlor_sccor(i,j)
+          v0ijsccor=0.0d0
+          si=-1.0d0
+         do k=1,nterm_sccor(i,j)
+           read (isccor,*,end=113,err=113) kk,v1sccor(k,i,j),v2sccor(k,i,j) 
+            v0ijsccor=v0ijsccor+si*v1sccor(k,i,j)
+            si=-si
+          enddo
+         do k=1,nlor_sccor(i,j)
+            read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
+     &        vlor2sccor(k,i,j),vlor3sccor(k,i,j) 
+            v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
+     &(1+vlor3sccor(k,i,j)**2)
           enddo
+          v0sccor(i,j)=v0ijsccor
         enddo
       enddo
+      enddo
       close (isccor)
+      
       if (lprint) then
-       write (iout,'(/a/)') 'Torsional constants of SCCORR:'
-       do i=1,20
-         do j=1,20
+       write (iout,'(/a/)') 'Torsional constants:'
+       do i=1,nsccortyp
+         do j=1,nsccortyp
             write (iout,*) 'ityp',i,' jtyp',j
-            do k=1,nterm_sccor
+            write (iout,*) 'Fourier constants'
+            do k=1,nterm_sccor(i,j)
              write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
             enddo
+            write (iout,*) 'Lorenz constants'
+            do k=1,nlor_sccor(i,j)
+             write (iout,'(3(1pe15.5))') 
+     &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
+            enddo
           enddo
         enddo
       endif
 C
+C
 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
 C         interaction energy of the Gly, Ala, and Pro prototypes.
 C