projects / unres.git / commitdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: 4c77976)
zabezpieczenie przedczarkowe
authorAdam Sieradzan <adasko@piasek2.(none)>
Thu, 14 Nov 2013 13:32:27 +0000 (14:32 +0100)
committerAdam Sieradzan <adasko@piasek2.(none)>
Thu, 14 Nov 2013 13:41:53 +0000 (14:41 +0100)
PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new [new file with mode: 0644] patch | blob
source/unres/src_MD-M/CMakeLists.txt patch | blob | history
source/unres/src_MD-M/COMMON.CONTACTS patch | blob | history
source/unres/src_MD-M/COMMON.TORSION patch | blob | history
source/unres/src_MD-M/energy_p_new_barrier.F patch | blob | history
source/unres/src_MD-M/parmread.F patch | blob | history

diff --git a/PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new b/PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new
new file mode 100644 (file)
index 0000000..231618e
--- /dev/null
+++ b/PARAM/fourier_opt_ext.parm.1igd_hc_iter3_3_new
@@ -0,0 +1,71 @@
+    3     # Number of local interaction types
+Gly
+   0.000000000000000
+   0.791965124028570
+   0.206068961118571
+   0.000000000000000
+   0.000000000000000
+   2.373462483972307
+  -0.927962753087420
+   0.000000000000000
+   0.000000000000000
+   1.329421814829764
+  -0.370576187607876
+   0.000000000000000
+   0.000000000000000
+   1.000000000000000
+   0.500000000000000
+   0.500000000000000
+   1.000000000000000
+   1.000000000000000
+   1.500000000000000
+   1.500000000000000
+   2.000000000000000
+   2.000000000000000
+Ala
+   0.000000000000000
+   0.500261572719827
+  -0.233786079150650
+  -0.878020534542259
+   1.501220349138902
+  -2.089734050079038
+   2.302365702770721
+  -0.532502145482045
+  -1.596421505165690
+   1.276301651241011
+   0.399942874780603
+  -0.543778421330412
+   0.400478498916355
+   1.000000000000000
+   0.500000000000000
+   0.500000000000000
+   1.000000000000000
+   1.000000000000000
+   1.500000000000000
+   1.500000000000000
+   2.000000000000000
+   2.000000000000000
+Pro
+   0.000000000000000
+  -1.286480000000000
+   0.031808800000000
+  -0.906628000000000
+   1.015390000000000
+  -1.548870000000000
+   1.914370000000000
+   0.664143000000000
+  -0.454839000000000
+  -0.051291400000000
+   0.103179000000000
+   0.316367000000000
+   0.045277200000000
+   1.000000000000000
+   0.500000000000000
+   0.500000000000000
+   1.000000000000000
+   1.000000000000000
+   1.500000000000000
+   1.500000000000000
+   2.000000000000000
+   2.000000000000000
+
diff --git a/source/unres/src_MD-M/CMakeLists.txt b/source/unres/src_MD-M/CMakeLists.txt
index cba4cc2..cfdc654 100644 (file)
--- a/source/unres/src_MD-M/CMakeLists.txt
+++ b/source/unres/src_MD-M/CMakeLists.txt
@@ -210,7 +210,7 @@ if(UNRES_MD_FF STREQUAL "GAB" )
 #=========================================
 elseif(UNRES_MD_FF STREQUAL "E0LL2Y")
   # set preprocesor flags   
-  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DSCCORPDB" )
+  set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DSCCORPDB -DNEWCORR" )
 endif(UNRES_MD_FF STREQUAL "GAB")
 
 
diff --git a/source/unres/src_MD-M/COMMON.CONTACTS b/source/unres/src_MD-M/COMMON.CONTACTS
index 5b3a90d..4b81714 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CONTACTS
+++ b/source/unres/src_MD-M/COMMON.CONTACTS
@@ -27,7 +27,7 @@ c     &  dipderx(3,5,4,maxconts,maxres)
 C 10/30/99 Added other pre-computed vectors and matrices needed 
 C          to calculate three - six-order el-loc correlation terms
       double precision Ug,Ugder,Ug2,Ug2der,obrot,obrot2,obrot_der,
-     &  obrot2_der,Ub2,Ub2der,mu,muder,EUg,EUgder,CUg,CUgder,
+     &  obrot2_der,Ub2,Ub2der,mu,muder,EUg,EUgder,CUg,CUgder,gmu,gUb2
      &  DUg,DUgder,DtUg2,DtUg2der,Ctobr,Ctobrder,Dtobr2,Dtobr2der
       common /rotat/ Ug(2,2,maxres),Ugder(2,2,maxres),Ug2(2,2,maxres),
      &  Ug2der(2,2,maxres),obrot(2,maxres),obrot2(2,maxres),
@@ -35,6 +35,7 @@ C          to calculate three - six-order el-loc correlation terms
 C This common block contains vectors and matrices dependent on a single
 C amino-acid residue.
       common /precomp1/ mu(2,maxres),muder(2,maxres),Ub2(2,maxres),
+     &  gmu(2,maxres),gUb2(2,maxres),
      &  Ub2der(2,maxres),Ctobr(2,maxres),Ctobrder(2,maxres),
      &  Dtobr2(2,maxres),Dtobr2der(2,maxres),
      &  EUg(2,2,maxres),EUgder(2,2,maxres),CUg(2,2,maxres),
diff --git a/source/unres/src_MD-M/COMMON.TORSION b/source/unres/src_MD-M/COMMON.TORSION
index 5f4b2b7..95472be 100644 (file)
--- a/source/unres/src_MD-M/COMMON.TORSION
+++ b/source/unres/src_MD-M/COMMON.TORSION
@@ -24,10 +24,17 @@ C 6/23/01 - constants for double torsionals
      &    ntermd_2(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
 C 9/18/99 - added Fourier coeffficients of the expansion of local energy 
 C           surface
-      double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde,b1tilde
+      double precision b1,b2,cc,dd,ee,ctilde,dtilde,b2tilde,b1tilde,
+     &bnew1,bnew2,eenew,ggb1,ggb2
       integer nloctyp
-      common/fourier/ b1(2,-maxtor:maxtor),b2(2,-maxtor:maxtor)
-     &    ,cc(2,2,-maxtor:maxtor),
+      common/fourier/ b1(2,maxres),b2(2,maxres),
+     &    bnew1(3,2,-maxtor:maxtor),bnew2(3,2,-maxtor:maxtor),
+     &    cc(2,2,-maxtor:maxtor),
      &    dd(2,2,-maxtor:maxtor),ee(2,2,-maxtor:maxtor),
+     &    eenew(2,-maxtor:maxtor),
      &    ctilde(2,2,-maxtor:maxtor),
-     &    dtilde(2,2,-maxtor:maxtor),b1tilde(2,-maxtor:maxtor),nloctyp
+     &    dtilde(2,2,-maxtor:maxtor),b1tilde(2,maxres),
+     &    b2tilde(2,maxres),
+     &    gtb1(2,maxres),gtb2(2,maxres),
+     &    nloctyp
+     
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 d4e2092..c13d341 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -2252,11 +2252,70 @@ C
 C Compute the virtual-bond-torsional-angle dependent quantities needed
 C to calculate the el-loc multibody terms of various order.
 C
+      write(iout,*) 'nphi=',nphi,nres
 #ifdef PARMAT
       do i=ivec_start+2,ivec_end+2
 #else
       do i=3,nres+1
 #endif
+#ifdef NEWCORR
+        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+          iti = itortyp(itype(i-2))
+        else
+          iti=ntortyp+1
+        endif
+c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+          iti1 = itortyp(itype(i-1))
+        else
+          iti1=ntortyp+1
+        endif
+        write(iout,*),i
+        b1(1,i-2)=bnew1(1,1,iti)*sin(theta(i-1)/2.0)
+     &           +bnew1(2,1,iti)*sin(theta(i-1))
+     &           +bnew1(3,1,iti)*cos(theta(i-1)/2.0)
+        gtb1(1,i-2)=bnew1(1,1,iti)*cos(theta(i-1)/2.0)/2.0
+     &             +bnew1(2,1,iti)*cos(theta(i-1))
+     &             -bnew1(3,1,iti)*sin(theta(i-1)/2.0)/2.0
+c     &           +bnew1(3,1,iti)*sin(alpha(i))*cos(beta(i))
+c     &*(cos(theta(i)/2.0)
+        b2(1,i-2)=bnew2(1,1,iti)*sin(theta(i-1)/2.0)
+     &           +bnew2(2,1,iti)*sin(theta(i-1))
+     &           +bnew2(3,1,iti)*cos(theta(i-1)/2.0)
+c     &           +bnew2(3,1,iti)*sin(alpha(i))*cos(beta(i))
+c     &*(cos(theta(i)/2.0)
+        gtb2(1,i-2)=bnew2(1,1,iti)*cos(theta(i-1)/2.0)/2.0
+     &             +bnew2(2,1,iti)*cos(theta(i-1))
+     &             -bnew2(3,1,iti)*sin(theta(i-1)/2.0)/2.0
+c        if (ggb1(1,i).eq.0.0d0) then
+c        write(iout,*) 'i=',i,ggb1(1,i),
+c     &bnew1(1,1,iti)*cos(theta(i)/2.0)/2.0,
+c     &bnew1(2,1,iti)*cos(theta(i)),
+c     &bnew1(3,1,iti)*sin(theta(i)/2.0)/2.0
+c        endif
+        b1(2,i-2)=bnew1(1,2,iti)
+        gtb1(2,i-2)=0.0
+        b2(2,i-2)=bnew2(1,2,iti)
+        gtb2(2,i-2)=0.0
+c        EE(1,1,iti)=0.0d0
+c        EE(2,2,iti)=0.0d0
+c        EE(1,2,iti)=0.5d0*eenew(1,iti)
+c        EE(2,1,iti)=0.5d0*eenew(1,iti)
+c        b1(2,iti)=bnew1(1,2,iti)*sin(alpha(i))*sin(beta(i))
+c        b2(2,iti)=bnew2(1,2,iti)*sin(alpha(i))*sin(beta(i))
+       b1tilde(1,i-2)=b1(1,i-2)
+       b1tilde(2,i-2)=-b1(2,i-2)
+       b2tilde(1,i-2)=b2(1,i-2)
+       b2tilde(2,i-2)=-b2(2,i-2)
+       write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2)
+       write (iout,*) 'theta=', theta(i-1)
+       enddo
+#ifdef PARMAT
+      do i=ivec_start+2,ivec_end+2
+#else
+      do i=3,nres+1
+#endif
+#endif
         if (i .lt. nres+1) then
           sin1=dsin(phi(i))
           cos1=dcos(phi(i))
@@ -2324,6 +2383,7 @@ C
           Ug2der(2,2,i-2)=0.0d0
         endif
 c        if (i.gt. iatel_s+2 .and. i.lt.iatel_e+5) then
+#ifndef NEWCORR
         if (i.gt. nnt+2 .and. i.lt.nct+2) then
           iti = itortyp(itype(i-2))
         else
@@ -2335,13 +2395,18 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
         else
           iti1=ntortyp+1
         endif
+#endif
 cd        write (iout,*) '*******i',i,' iti1',iti
 cd        write (iout,*) 'b1',b1(:,iti)
 cd        write (iout,*) 'b2',b2(:,iti)
 cd        write (iout,*) 'Ug',Ug(:,:,i-2)
 c        if (i .gt. iatel_s+2) then
         if (i .gt. nnt+2) then
-          call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2))
+          call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2))
+#ifdef NEWCORR
+          call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2))
+c          write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj"
+#endif
           call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
           if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) 
      &    then
@@ -2364,7 +2429,7 @@ c        if (i .gt. iatel_s+2) then
             enddo
           enddo
         endif
-        call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2))
+        call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2))
         call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2))
         do k=1,2
           muder(k,i-2)=Ub2der(k,i-2)
@@ -2380,7 +2445,7 @@ c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
           iti1=ntortyp+1
         endif
         do k=1,2
-          mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
+          mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1)
         enddo
 cd        write (iout,*) 'mu ',mu(:,i-2)
 cd        write (iout,*) 'mu1',mu1(:,i-2)
@@ -2393,7 +2458,7 @@ cd        write (iout,*) 'mu2',mu2(:,i-2)
         call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2))
         call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2))
 C Vectors and matrices dependent on a single virtual-bond dihedral.
-        call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
+        call matvec2(DD(1,1,iti),b1tilde(1,i-1),auxvec(1))
         call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) 
         call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) 
         call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2))
@@ -2709,7 +2774,7 @@ C
       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),
-     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij(4)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
@@ -2801,7 +2866,8 @@ C
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
         num_conti=0
-        call eelecij(i,i+2,ees,evdw1,eel_loc)
+c TU ZLE
+c        call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
         num_cont_hb(i)=num_conti
       enddo
@@ -2819,7 +2885,8 @@ C
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
         num_conti=num_cont_hb(i)
-        call eelecij(i,i+3,ees,evdw1,eel_loc)
+c TU ZLE
+c        call eelecij(i,i+3,ees,evdw1,eel_loc)
         if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) 
      &   call eturn4(i,eello_turn4)
         num_cont_hb(i)=num_conti
@@ -2827,7 +2894,8 @@ C
 c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
-      do i=iatel_s,iatel_e
+c      do i=iatel_s,iatel_e
+       do i=4,5
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -2840,8 +2908,9 @@ c
         zmedi=c(3,i)+0.5d0*dzi
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
-        do j=ielstart(i),ielend(i)
-c          write (iout,*) i,j,itype(i),itype(j)
+c        do j=ielstart(i),ielend(i)
+         do j=8,9
+          write (iout,*) 'tu wchodze',i,j,itype(i),itype(j)
           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
         enddo ! j
@@ -2880,7 +2949,8 @@ C-------------------------------------------------------------------------------
       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),
-     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4)
+     &    aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
+     &    gmuij2(4),gmuji2(4)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
@@ -3097,6 +3167,7 @@ C   Fourier series in the angles lambda1 and lambda2 (see Nishikawa et al.
 C   Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
 C   are computed for EVERY pair of non-contiguous peptide groups.
 C
+
           if (j.lt.nres-1) then
             j1=j+1
             j2=j-1
@@ -3105,10 +3176,19 @@ C
             j2=j-2
           endif
           kkk=0
+          lll=0
           do k=1,2
             do l=1,2
               kkk=kkk+1
               muij(kkk)=mu(k,i)*mu(l,j)
+#ifdef NEWCORR
+             gmuij1(kkk)=gtb1(k,i)*mu(l,j)
+             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(k,i),k,i
+             gmuij2(kkk)=gUb2(k,i-1)*mu(l,j)
+             gmuji1(kkk)=mu(k,i)*gtb1(l,j)
+             write(iout,*) 'kkk=', gtb1(k,i)*mu(l,j),gtb1(l,j),l,j
+             gmuji2(kkk)=mu(k,i)*gUb2(l,j-1)
+#endif
             enddo
           enddo  
 cd         write (iout,*) 'EELEC: i',i,' j',j
@@ -3274,6 +3354,40 @@ cgrad            endif
 C Contribution to the local-electrostatic energy coming from the i-j pair
           eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
      &     +a33*muij(4)
+C Calculate patrial derivative for theta angle
+#ifdef NEWCORR
+         geel_loc_ij=a22*gmuij1(1)
+     &     +a23*gmuij1(2)
+     &     +a32*gmuij1(3)
+     &     +a33*gmuij1(4)         
+         write(iout,*) "derivative over thatai"
+         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
+     &   a33*gmuij1(4) 
+         gloc(nphi+i,icg)=gloc(nphi+i,icg)+
+     &      geel_loc_ij*wel_loc
+         write(iout,*) "derivative over thatai-1" 
+         write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
+     &   a33*gmuij2(4)
+         geel_loc_ij=a22*gmuij2(1)+a23*gmuij2(2)+a32*gmuij2(3)
+     &     +a33*gmuij2(4)
+         gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
+     &      geel_loc_ij*wel_loc
+         geel_loc_ji=a22*gmuji1(1)+a23*gmuji1(2)+a32*gmuji1(3)
+     &     +a33*gmuji1(4)
+         write(iout,*) "derivative over thataj" 
+         write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3),
+     &   a33*gmuji1(4)
+
+         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
+     &      geel_loc_ji*wel_loc
+         geel_loc_ji=a22*gmuji2(1)+a23*gmuji2(2)+a32*gmuji2(3)
+     &     +a33*gmuji2(4)
+         write(iout,*) "derivative over thataj-1"
+         write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
+     &   a33*gmuji2(4)
+         gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
+     &      geel_loc_ji*wel_loc
+#endif
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
           if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
@@ -3626,7 +3740,7 @@ C Third- and fourth-order contributions from turns
       dimension ggg(3)
       double precision auxmat(2,2),auxmat1(2,2),auxmat2(2,2),pizda(2,2),
      &  e1t(2,2),e2t(2,2),e3t(2,2),e1tder(2,2),e2tder(2,2),e3tder(2,2),
-     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2)
+     &  e1a(2,2),ae3(2,2),ae3e2(2,2),auxvec(2),auxvec1(2),auxgvec(2)
       double precision agg(3,4),aggi(3,4),aggi1(3,4),
      &    aggj(3,4),aggj1(3,4),a_temp(2,2),auxmat3(2,2)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
@@ -3659,14 +3773,31 @@ c        write(iout,*) "iti1",iti1," iti2",iti2," iti3",iti3
         call transpose2(Eug(1,1,i+3),e3t(1,1))
         call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
         call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
-        s1=scalar2(b1(1,iti2),auxvec(1))
+        call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1)) 
+c        s1=0.0
+c        gs1=0.0    
+        s1=scalar2(b1(1,i+2),auxvec(1))
+c        gs1=scalar2(gtb1(1,i+2),auxgvec(1))
         call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
         call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
-        s2=scalar2(b1(1,iti1),auxvec(1))
+        call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1))
+c        s2=0.0
+c        gs2=0.0
+        s2=scalar2(b1(1,i+1),auxvec(1))
+c        gs2=scalar2(gtb1(1,i+1),auxgvec(1))
+c        write(iout,*) gs1,gs2,'i=',i,auxgvec(1),gUb2(1,i+2),ggb1(1,i+2),
+c     &  ggb1(1,i+1)
         call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
         call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         eello_turn4=eello_turn4-(s1+s2+s3)
+#ifdef NEWCORR
+c        geel_loc_ij=-(gs1+gs2)
+c         gloc(nphi+i,icg)=gloc(nphi+i,icg)-
+c     &   gs1
+c         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
+c     &   gs2
+#endif
         if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
      &      'eturn4',i,j,-(s1+s2+s3)
 cd        write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3),
@@ -3675,14 +3806,14 @@ C Derivatives in gamma(i)
         call transpose2(EUgder(1,1,i+1),e1tder(1,1))
         call matmat2(e1tder(1,1),a_temp(1,1),auxmat(1,1))
         call matvec2(auxmat(1,1),Ub2(1,i+3),auxvec(1))
-        s1=scalar2(b1(1,iti2),auxvec(1))
+        s1=scalar2(b1(1,i+2),auxvec(1))
         call matmat2(ae3e2(1,1),e1tder(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
 C Derivatives in gamma(i+1)
         call transpose2(EUgder(1,1,i+2),e2tder(1,1))
         call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) 
-        s2=scalar2(b1(1,iti1),auxvec(1))
+        s2=scalar2(b1(1,i+1),auxvec(1))
         call matmat2(ae3(1,1),e2tder(1,1),auxmat(1,1))
         call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@ -3690,10 +3821,10 @@ C Derivatives in gamma(i+1)
 C Derivatives in gamma(i+2)
         call transpose2(EUgder(1,1,i+3),e3tder(1,1))
         call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1))
-        s1=scalar2(b1(1,iti2),auxvec(1))
+        s1=scalar2(b1(1,i+2),auxvec(1))
         call matmat2(a_temp(1,1),e3tder(1,1),auxmat(1,1))
         call matvec2(auxmat(1,1),Ub2(1,i+2),auxvec(1)) 
-        s2=scalar2(b1(1,iti1),auxvec(1))
+        s2=scalar2(b1(1,i+1),auxvec(1))
         call matmat2(auxmat(1,1),e2t(1,1),auxmat3(1,1))
         call matmat2(auxmat3(1,1),e1t(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@ -3708,10 +3839,10 @@ C Derivatives of this turn contributions in DC(i+2)
             a_temp(2,2)=agg(l,4)
             call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
             call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
-            s1=scalar2(b1(1,iti2),auxvec(1))
+            s1=scalar2(b1(1,i+2),auxvec(1))
             call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
             call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
-            s2=scalar2(b1(1,iti1),auxvec(1))
+            s2=scalar2(b1(1,i+1),auxvec(1))
             call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
             call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
             s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@ -3727,10 +3858,10 @@ C Remaining derivatives of this turn contribution
           a_temp(2,2)=aggi(l,4)
           call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
           call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
-          s1=scalar2(b1(1,iti2),auxvec(1))
+          s1=scalar2(b1(1,i+2),auxvec(1))
           call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
           call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
-          s2=scalar2(b1(1,iti1),auxvec(1))
+          s2=scalar2(b1(1,i+1),auxvec(1))
           call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
           call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@ -3741,10 +3872,10 @@ C Remaining derivatives of this turn contribution
           a_temp(2,2)=aggi1(l,4)
           call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
           call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
-          s1=scalar2(b1(1,iti2),auxvec(1))
+          s1=scalar2(b1(1,i+2),auxvec(1))
           call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
           call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
-          s2=scalar2(b1(1,iti1),auxvec(1))
+          s2=scalar2(b1(1,i+1),auxvec(1))
           call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
           call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@ -3755,10 +3886,10 @@ C Remaining derivatives of this turn contribution
           a_temp(2,2)=aggj(l,4)
           call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
           call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
-          s1=scalar2(b1(1,iti2),auxvec(1))
+          s1=scalar2(b1(1,i+2),auxvec(1))
           call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
           call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
-          s2=scalar2(b1(1,iti1),auxvec(1))
+          s2=scalar2(b1(1,i+1),auxvec(1))
           call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
           call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@ -3769,10 +3900,10 @@ C Remaining derivatives of this turn contribution
           a_temp(2,2)=aggj1(l,4)
           call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1))
           call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1))
-          s1=scalar2(b1(1,iti2),auxvec(1))
+          s1=scalar2(b1(1,i+2),auxvec(1))
           call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1))
           call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) 
-          s2=scalar2(b1(1,iti1),auxvec(1))
+          s2=scalar2(b1(1,i+1),auxvec(1))
           call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1))
           call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1))
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@ -6816,10 +6947,10 @@ C---------------------------------------------------------------------------
       do iii=1,2
         dipi(iii,1)=Ub2(iii,i)
         dipderi(iii)=Ub2der(iii,i)
-        dipi(iii,2)=b1(iii,iti1)
+        dipi(iii,2)=b1(iii,i+1)
         dipj(iii,1)=Ub2(iii,j)
         dipderj(iii)=Ub2der(iii,j)
-        dipj(iii,2)=b1(iii,itj1)
+        dipj(iii,2)=b1(iii,j+1)
       enddo
       kkk=0
       do iii=1,2
@@ -6999,26 +7130,26 @@ C They are needed only when the fifth- or the sixth-order cumulants are
 C indluded.
         IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) THEN
         call transpose2(AEA(1,1,1),auxmat(1,1))
-        call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
+        call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1))
         call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
         call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
         call transpose2(AEAderg(1,1,1),auxmat(1,1))
-        call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
+        call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1))
         call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
-        call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
-        call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
+        call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1))
+        call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1))
         call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
         call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
         call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
         call transpose2(AEA(1,1,2),auxmat(1,1))
-        call matvec2(auxmat(1,1),b1(1,itj),AEAb1(1,1,2))
+        call matvec2(auxmat(1,1),b1(1,j),AEAb1(1,1,2))
         call matvec2(auxmat(1,1),Ub2(1,j),AEAb2(1,1,2))
         call matvec2(auxmat(1,1),Ub2der(1,j),AEAb2derg(1,2,1,2))
         call transpose2(AEAderg(1,1,2),auxmat(1,1))
-        call matvec2(auxmat(1,1),b1(1,itj),AEAb1derg(1,1,2))
+        call matvec2(auxmat(1,1),b1(1,j),AEAb1derg(1,1,2))
         call matvec2(auxmat(1,1),Ub2(1,j),AEAb2derg(1,1,1,2))
-        call matvec2(AEA(1,1,2),b1(1,itl1),AEAb1(1,2,2))
-        call matvec2(AEAderg(1,1,2),b1(1,itl1),AEAb1derg(1,2,2))
+        call matvec2(AEA(1,1,2),b1(1,l+1),AEAb1(1,2,2))
+        call matvec2(AEAderg(1,1,2),b1(1,l+1),AEAb1derg(1,2,2))
         call matvec2(AEA(1,1,2),Ub2(1,l+1),AEAb2(1,2,2))
         call matvec2(AEAderg(1,1,2),Ub2(1,l+1),AEAb2derg(1,1,2,2))
         call matvec2(AEA(1,1,2),Ub2der(1,l+1),AEAb2derg(1,2,2,2))
@@ -7027,20 +7158,20 @@ C Calculate the Cartesian derivatives of the vectors.
           do kkk=1,5
             do lll=1,3
               call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
-              call matvec2(auxmat(1,1),b1(1,iti),
+              call matvec2(auxmat(1,1),b1(1,i),
      &          AEAb1derx(1,lll,kkk,iii,1,1))
               call matvec2(auxmat(1,1),Ub2(1,i),
      &          AEAb2derx(1,lll,kkk,iii,1,1))
-              call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+              call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
      &          AEAb1derx(1,lll,kkk,iii,2,1))
               call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
      &          AEAb2derx(1,lll,kkk,iii,2,1))
               call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
-              call matvec2(auxmat(1,1),b1(1,itj),
+              call matvec2(auxmat(1,1),b1(1,j),
      &          AEAb1derx(1,lll,kkk,iii,1,2))
               call matvec2(auxmat(1,1),Ub2(1,j),
      &          AEAb2derx(1,lll,kkk,iii,1,2))
-              call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
+              call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,l+1),
      &          AEAb1derx(1,lll,kkk,iii,2,2))
               call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,l+1),
      &          AEAb2derx(1,lll,kkk,iii,2,2))
@@ -7137,26 +7268,26 @@ C indluded.
         IF (wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0 .or.
      &    (wturn6.gt.0.0d0 .and. j.eq.i+4 .and. l.eq.i+3)) THEN
         call transpose2(AEA(1,1,1),auxmat(1,1))
-        call matvec2(auxmat(1,1),b1(1,iti),AEAb1(1,1,1))
+        call matvec2(auxmat(1,1),b1(1,i),AEAb1(1,1,1))
         call matvec2(auxmat(1,1),Ub2(1,i),AEAb2(1,1,1))
         call matvec2(auxmat(1,1),Ub2der(1,i),AEAb2derg(1,2,1,1))
         call transpose2(AEAderg(1,1,1),auxmat(1,1))
-        call matvec2(auxmat(1,1),b1(1,iti),AEAb1derg(1,1,1))
+        call matvec2(auxmat(1,1),b1(1,i),AEAb1derg(1,1,1))
         call matvec2(auxmat(1,1),Ub2(1,i),AEAb2derg(1,1,1,1))
-        call matvec2(AEA(1,1,1),b1(1,itk1),AEAb1(1,2,1))
-        call matvec2(AEAderg(1,1,1),b1(1,itk1),AEAb1derg(1,2,1))
+        call matvec2(AEA(1,1,1),b1(1,k+1),AEAb1(1,2,1))
+        call matvec2(AEAderg(1,1,1),b1(1,k+1),AEAb1derg(1,2,1))
         call matvec2(AEA(1,1,1),Ub2(1,k+1),AEAb2(1,2,1))
         call matvec2(AEAderg(1,1,1),Ub2(1,k+1),AEAb2derg(1,1,2,1))
         call matvec2(AEA(1,1,1),Ub2der(1,k+1),AEAb2derg(1,2,2,1))
         call transpose2(AEA(1,1,2),auxmat(1,1))
-        call matvec2(auxmat(1,1),b1(1,itj1),AEAb1(1,1,2))
+        call matvec2(auxmat(1,1),b1(1,j+1),AEAb1(1,1,2))
         call matvec2(auxmat(1,1),Ub2(1,l),AEAb2(1,1,2))
         call matvec2(auxmat(1,1),Ub2der(1,l),AEAb2derg(1,2,1,2))
         call transpose2(AEAderg(1,1,2),auxmat(1,1))
-        call matvec2(auxmat(1,1),b1(1,itl),AEAb1(1,1,2))
+        call matvec2(auxmat(1,1),b1(1,l),AEAb1(1,1,2))
         call matvec2(auxmat(1,1),Ub2(1,l),AEAb2derg(1,1,1,2))
-        call matvec2(AEA(1,1,2),b1(1,itj1),AEAb1(1,2,2))
-        call matvec2(AEAderg(1,1,2),b1(1,itj1),AEAb1derg(1,2,2))
+        call matvec2(AEA(1,1,2),b1(1,j+1),AEAb1(1,2,2))
+        call matvec2(AEAderg(1,1,2),b1(1,j+1),AEAb1derg(1,2,2))
         call matvec2(AEA(1,1,2),Ub2(1,j),AEAb2(1,2,2))
         call matvec2(AEAderg(1,1,2),Ub2(1,j),AEAb2derg(1,1,2,2))
         call matvec2(AEA(1,1,2),Ub2der(1,j),AEAb2derg(1,2,2,2))
@@ -7165,20 +7296,20 @@ C Calculate the Cartesian derivatives of the vectors.
           do kkk=1,5
             do lll=1,3
               call transpose2(AEAderx(1,1,lll,kkk,iii,1),auxmat(1,1))
-              call matvec2(auxmat(1,1),b1(1,iti),
+              call matvec2(auxmat(1,1),b1(1,i),
      &          AEAb1derx(1,lll,kkk,iii,1,1))
               call matvec2(auxmat(1,1),Ub2(1,i),
      &          AEAb2derx(1,lll,kkk,iii,1,1))
-              call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+              call matvec2(AEAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
      &          AEAb1derx(1,lll,kkk,iii,2,1))
               call matvec2(AEAderx(1,1,lll,kkk,iii,1),Ub2(1,k+1),
      &          AEAb2derx(1,lll,kkk,iii,2,1))
               call transpose2(AEAderx(1,1,lll,kkk,iii,2),auxmat(1,1))
-              call matvec2(auxmat(1,1),b1(1,itl),
+              call matvec2(auxmat(1,1),b1(1,l),
      &          AEAb1derx(1,lll,kkk,iii,1,2))
               call matvec2(auxmat(1,1),Ub2(1,l),
      &          AEAb2derx(1,lll,kkk,iii,1,2))
-              call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,itj1),
+              call matvec2(AEAderx(1,1,lll,kkk,iii,2),b1(1,j+1),
      &          AEAb1derx(1,lll,kkk,iii,2,2))
               call matvec2(AEAderx(1,1,lll,kkk,iii,2),Ub2(1,j),
      &          AEAb2derx(1,lll,kkk,iii,2,2))
@@ -7475,7 +7606,7 @@ C Contribution from graph II
       call matmat2(auxmat(1,1),AEA(1,1,1),pizda(1,1))
       vv(1)=pizda(1,1)+pizda(2,2)
       vv(2)=pizda(2,1)-pizda(1,2)
-      eello5_2=scalar2(AEAb1(1,2,1),b1(1,itk))
+      eello5_2=scalar2(AEAb1(1,2,1),b1(1,k))
      & -0.5d0*scalar2(vv(1),Ctobr(1,k))
 C Explicit gradient in virtual-dihedral angles.
       g_corr5_loc(k-1)=g_corr5_loc(k-1)
@@ -7485,11 +7616,11 @@ C Explicit gradient in virtual-dihedral angles.
       vv(2)=pizda(2,1)-pizda(1,2)
       if (l.eq.j+1) then
         g_corr5_loc(l-1)=g_corr5_loc(l-1)
-     &   +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
+     &   +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k))
      &   -0.5d0*scalar2(vv(1),Ctobr(1,k)))
       else
         g_corr5_loc(j-1)=g_corr5_loc(j-1)
-     &   +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,itk))
+     &   +ekont*(scalar2(AEAb1derg(1,2,1),b1(1,k))
      &   -0.5d0*scalar2(vv(1),Ctobr(1,k)))
       endif
 C Cartesian gradient
@@ -7501,7 +7632,7 @@ C Cartesian gradient
             vv(1)=pizda(1,1)+pizda(2,2)
             vv(2)=pizda(2,1)-pizda(1,2)
             derx(lll,kkk,iii)=derx(lll,kkk,iii)
-     &       +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,itk))
+     &       +scalar2(AEAb1derx(1,lll,kkk,iii,2,1),b1(1,k))
      &       -0.5d0*scalar2(vv(1),Ctobr(1,k))
           enddo
         enddo
@@ -7556,7 +7687,7 @@ cd1110    continue
         call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
         vv(1)=pizda(1,1)+pizda(2,2)
         vv(2)=pizda(2,1)-pizda(1,2)
-        eello5_4=scalar2(AEAb1(1,2,2),b1(1,itl))
+        eello5_4=scalar2(AEAb1(1,2,2),b1(1,l))
      &   -0.5d0*scalar2(vv(1),Ctobr(1,l))
 C Explicit gradient in virtual-dihedral angles.
         g_corr5_loc(l-1)=g_corr5_loc(l-1)
@@ -7565,7 +7696,7 @@ C Explicit gradient in virtual-dihedral angles.
         vv(1)=pizda(1,1)+pizda(2,2)
         vv(2)=pizda(2,1)-pizda(1,2)
         g_corr5_loc(k-1)=g_corr5_loc(k-1)
-     &   +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itl))
+     &   +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,l))
      &   -0.5d0*scalar2(vv(1),Ctobr(1,l)))
 C Cartesian gradient
         do iii=1,2
@@ -7576,7 +7707,7 @@ C Cartesian gradient
               vv(1)=pizda(1,1)+pizda(2,2)
               vv(2)=pizda(2,1)-pizda(1,2)
               derx(lll,kkk,iii)=derx(lll,kkk,iii)
-     &         +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itl))
+     &         +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,l))
      &         -0.5d0*scalar2(vv(1),Ctobr(1,l))
             enddo
           enddo
@@ -7629,7 +7760,7 @@ C Contribution from graph IV
         call matmat2(auxmat(1,1),AEA(1,1,2),pizda(1,1))
         vv(1)=pizda(1,1)+pizda(2,2)
         vv(2)=pizda(2,1)-pizda(1,2)
-        eello5_4=scalar2(AEAb1(1,2,2),b1(1,itj))
+        eello5_4=scalar2(AEAb1(1,2,2),b1(1,j))
      &   -0.5d0*scalar2(vv(1),Ctobr(1,j))
 C Explicit gradient in virtual-dihedral angles.
         g_corr5_loc(j-1)=g_corr5_loc(j-1)
@@ -7638,7 +7769,7 @@ C Explicit gradient in virtual-dihedral angles.
         vv(1)=pizda(1,1)+pizda(2,2)
         vv(2)=pizda(2,1)-pizda(1,2)
         g_corr5_loc(k-1)=g_corr5_loc(k-1)
-     &   +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,itj))
+     &   +ekont*(scalar2(AEAb1derg(1,2,2),b1(1,j))
      &   -0.5d0*scalar2(vv(1),Ctobr(1,j)))
 C Cartesian gradient
         do iii=1,2
@@ -7649,7 +7780,7 @@ C Cartesian gradient
               vv(1)=pizda(1,1)+pizda(2,2)
               vv(2)=pizda(2,1)-pizda(1,2)
               derx(lll,kkk,3-iii)=derx(lll,kkk,3-iii)
-     &         +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,itj))
+     &         +scalar2(AEAb1derx(1,lll,kkk,iii,2,2),b1(1,j))
      &         -0.5d0*scalar2(vv(1),Ctobr(1,j))
             enddo
           enddo
@@ -7931,8 +8062,8 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       vv1(1)=pizda1(1,1)-pizda1(2,2)
       vv1(2)=pizda1(1,2)+pizda1(2,1)
       s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
-      vv(1)=AEAb1(1,2,imat)*b1(1,itk)-AEAb1(2,2,imat)*b1(2,itk)
-      vv(2)=AEAb1(1,2,imat)*b1(2,itk)+AEAb1(2,2,imat)*b1(1,itk)
+      vv(1)=AEAb1(1,2,imat)*b1(1,k)-AEAb1(2,2,imat)*b1(2,k)
+      vv(2)=AEAb1(1,2,imat)*b1(2,k)+AEAb1(2,2,imat)*b1(1,k)
       s5=scalar2(vv(1),Dtobr2(1,i))
 cd      write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5
       eello6_graph1=-0.5d0*(s1+s2+s3+s4+s5)
@@ -7945,8 +8076,8 @@ cd      write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5
       call matmat2(AEAderg(1,1,imat),auxmat(1,1),pizda1(1,1))
       vv1(1)=pizda1(1,1)-pizda1(2,2)
       vv1(2)=pizda1(1,2)+pizda1(2,1)
-      vv(1)=AEAb1derg(1,2,imat)*b1(1,itk)-AEAb1derg(2,2,imat)*b1(2,itk)
-      vv(2)=AEAb1derg(1,2,imat)*b1(2,itk)+AEAb1derg(2,2,imat)*b1(1,itk)
+      vv(1)=AEAb1derg(1,2,imat)*b1(1,k)-AEAb1derg(2,2,imat)*b1(2,k)
+      vv(2)=AEAb1derg(1,2,imat)*b1(2,k)+AEAb1derg(2,2,imat)*b1(1,k)
       if (l.eq.j+1) then
         g_corr6_loc(l-1)=g_corr6_loc(l-1)
      & +ekont*(-0.5d0*(scalar2(AEAb1derg(1,2,imat),CUgb2(1,i))
@@ -7985,10 +8116,10 @@ cd      write (2,*) 's1',s1,' s2',s2,' s3',s3,' s4', s4,' s5',s5
             vv1(1)=pizda1(1,1)-pizda1(2,2)
             vv1(2)=pizda1(1,2)+pizda1(2,1)
             s4=0.5d0*scalar2(vv1(1),Dtobr2(1,i))
-            vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,itk)
-     &       -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,itk)
-            vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,itk)
-     &       +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,itk)
+            vv(1)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(1,k)
+     &       -AEAb1derx(2,lll,kkk,iii,2,imat)*b1(2,k)
+            vv(2)=AEAb1derx(1,lll,kkk,iii,2,imat)*b1(2,k)
+     &       +AEAb1derx(2,lll,kkk,iii,2,imat)*b1(1,k)
             s5=scalar2(vv(1),Dtobr2(1,i))
             derx(lll,kkk,ind)=derx(lll,kkk,ind)-0.5d0*(s1+s2+s3+s4+s5)
           enddo
@@ -8228,10 +8359,10 @@ C           energy moment and not to the cluster cumulant.
 #ifdef MOMENT
       s1=dip(4,jj,i)*dip(4,kk,k)
 #endif
-      call matvec2(AECA(1,1,1),b1(1,itk1),auxvec(1))
-      s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
-      call matvec2(AECA(1,1,2),b1(1,itl1),auxvec(1))
-      s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+      call matvec2(AECA(1,1,1),b1(1,k+1),auxvec(1))
+      s2=0.5d0*scalar2(b1(1,k),auxvec(1))
+      call matvec2(AECA(1,1,2),b1(1,l+1),auxvec(1))
+      s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
       call transpose2(EE(1,1,itk),auxmat(1,1))
       call matmat2(auxmat(1,1),AECA(1,1,1),pizda(1,1))
       vv(1)=pizda(1,1)+pizda(2,2)
@@ -8246,13 +8377,13 @@ cd     & "sum",-(s2+s3+s4)
 #endif
 c      eello6_graph3=-s4
 C Derivatives in gamma(k-1)
-      call matvec2(AECAderg(1,1,2),b1(1,itl1),auxvec(1))
-      s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+      call matvec2(AECAderg(1,1,2),b1(1,l+1),auxvec(1))
+      s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
       s4=-0.25d0*scalar2(vv(1),Ctobrder(1,k))
       g_corr6_loc(k-1)=g_corr6_loc(k-1)-ekont*(s3+s4)
 C Derivatives in gamma(l-1)
-      call matvec2(AECAderg(1,1,1),b1(1,itk1),auxvec(1))
-      s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
+      call matvec2(AECAderg(1,1,1),b1(1,k+1),auxvec(1))
+      s2=0.5d0*scalar2(b1(1,k),auxvec(1))
       call matmat2(auxmat(1,1),AECAderg(1,1,1),pizda(1,1))
       vv(1)=pizda(1,1)+pizda(2,2)
       vv(2)=pizda(2,1)-pizda(1,2)
@@ -8269,12 +8400,12 @@ C Cartesian derivatives.
               s1=dip(4,jj,i)*dipderx(lll,kkk,4,kk,k)
             endif
 #endif
-            call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,itk1),
+            call matvec2(AECAderx(1,1,lll,kkk,iii,1),b1(1,k+1),
      &        auxvec(1))
-            s2=0.5d0*scalar2(b1(1,itk),auxvec(1))
-            call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,itl1),
+            s2=0.5d0*scalar2(b1(1,k),auxvec(1))
+            call matvec2(AECAderx(1,1,lll,kkk,iii,2),b1(1,l+1),
      &        auxvec(1))
-            s3=0.5d0*scalar2(b1(1,itj1),auxvec(1))
+            s3=0.5d0*scalar2(b1(1,j+1),auxvec(1))
             call matmat2(auxmat(1,1),AECAderx(1,1,lll,kkk,iii,1),
      &        pizda(1,1))
             vv(1)=pizda(1,1)+pizda(2,2)
@@ -8362,11 +8493,11 @@ cd     & ' itl',itl,' itl1',itl1
       call matvec2(AECA(1,1,imat),Ub2(1,k),auxvec(1))
       s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
       if (j.eq.l+1) then
-        call matvec2(ADtEA1(1,1,3-imat),b1(1,itj1),auxvec1(1))
-        s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+        call matvec2(ADtEA1(1,1,3-imat),b1(1,j+1),auxvec1(1))
+        s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
       else
-        call matvec2(ADtEA1(1,1,3-imat),b1(1,itl1),auxvec1(1))
-        s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+        call matvec2(ADtEA1(1,1,3-imat),b1(1,l+1),auxvec1(1))
+        s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
       endif
       call transpose2(EUg(1,1,k),auxmat(1,1))
       call matmat2(AECA(1,1,imat),auxmat(1,1),pizda(1,1))
@@ -8390,11 +8521,11 @@ C Derivatives in gamma(i-1)
 #endif
         s2=0.5d0*scalar2(Ub2der(1,i),auxvec(1))
         if (j.eq.l+1) then
-          call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itj1),auxvec1(1))
-          s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+          call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,j+1),auxvec1(1))
+          s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
         else
-          call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,itl1),auxvec1(1))
-          s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+          call matvec2(ADtEA1derg(1,1,1,3-imat),b1(1,l+1),auxvec1(1))
+          s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
         endif
         s4=0.25d0*scalar2(vv(1),Dtobr2der(1,i))
         if (wturn6.gt.0.0d0 .and. k.eq.l+4 .and. i.eq.j+2) then
@@ -8423,11 +8554,11 @@ C Derivatives in gamma(k-1)
       call matvec2(AECA(1,1,imat),Ub2der(1,k),auxvec1(1))
       s2=0.5d0*scalar2(Ub2(1,i),auxvec1(1))
       if (j.eq.l+1) then
-        call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itj1),auxvec1(1))
-        s3=-0.5d0*scalar2(b1(1,itj),auxvec1(1))
+        call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,j+1),auxvec1(1))
+        s3=-0.5d0*scalar2(b1(1,j),auxvec1(1))
       else
-        call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,itl1),auxvec1(1))
-        s3=-0.5d0*scalar2(b1(1,itl),auxvec1(1))
+        call matvec2(ADtEA1derg(1,1,2,3-imat),b1(1,l+1),auxvec1(1))
+        s3=-0.5d0*scalar2(b1(1,l),auxvec1(1))
       endif
       call transpose2(EUgder(1,1,k),auxmat1(1,1))
       call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1))
@@ -8493,12 +8624,12 @@ C Cartesian derivatives.
             s2=0.5d0*scalar2(Ub2(1,i),auxvec(1))
             if (j.eq.l+1) then
               call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
-     &          b1(1,itj1),auxvec(1))
-              s3=-0.5d0*scalar2(b1(1,itj),auxvec(1))
+     &          b1(1,j+1),auxvec(1))
+              s3=-0.5d0*scalar2(b1(1,j),auxvec(1))
             else
               call matvec2(ADtEA1derx(1,1,lll,kkk,iii,3-imat),
-     &          b1(1,itl1),auxvec(1))
-              s3=-0.5d0*scalar2(b1(1,itl),auxvec(1))
+     &          b1(1,l+1),auxvec(1))
+              s3=-0.5d0*scalar2(b1(1,l),auxvec(1))
             endif
             call matmat2(AECAderx(1,1,lll,kkk,iii,imat),auxmat(1,1),
      &        pizda(1,1))
@@ -8598,12 +8729,12 @@ cd      write (2,*) 'eello6_5',eello6_5
 #ifdef MOMENT
       call transpose2(AEA(1,1,1),auxmat(1,1))
       call matmat2(EUg(1,1,i+1),auxmat(1,1),auxmat(1,1))
-      ss1=scalar2(Ub2(1,i+2),b1(1,itl))
+      ss1=scalar2(Ub2(1,i+2),b1(1,l))
       s1 = (auxmat(1,1)+auxmat(2,2))*ss1
 #endif
-      call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
+      call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1))
       call matvec2(AEA(1,1,1),vtemp1(1),vtemp1(1))
-      s2 = scalar2(b1(1,itk),vtemp1(1))
+      s2 = scalar2(b1(1,k),vtemp1(1))
 #ifdef MOMENT
       call transpose2(AEA(1,1,2),atemp(1,1))
       call matmat2(atemp(1,1),EUg(1,1,i+4),atemp(1,1))
@@ -8618,7 +8749,7 @@ cd      write (2,*) 'eello6_5',eello6_5
       call matmat2(achuj_temp(1,1),EUg(1,1,i+2),gtemp(1,1))
       call matmat2(gtemp(1,1),EUg(1,1,i+3),gtemp(1,1)) 
       call matvec2(a_chuj(1,1,jj,i),Ub2(1,i+4),vtemp4(1)) 
-      ss13 = scalar2(b1(1,itk),vtemp4(1))
+      ss13 = scalar2(b1(1,k),vtemp4(1))
       s13 = (gtemp(1,1)+gtemp(2,2))*ss13
 #endif
 c      write (2,*) 's1,s2,s8,s12,s13',s1,s2,s8,s12,s13
@@ -8652,12 +8783,12 @@ C Derivatives in gamma(i+3)
 #ifdef MOMENT
       call transpose2(AEA(1,1,1),auxmatd(1,1))
       call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
-      ss1d=scalar2(Ub2der(1,i+2),b1(1,itl))
+      ss1d=scalar2(Ub2der(1,i+2),b1(1,l))
       s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1d
 #endif
-      call matvec2(EUgder(1,1,i+2),b1(1,itl),vtemp1d(1))
+      call matvec2(EUgder(1,1,i+2),b1(1,l),vtemp1d(1))
       call matvec2(AEA(1,1,1),vtemp1d(1),vtemp1d(1))
-      s2d = scalar2(b1(1,itk),vtemp1d(1))
+      s2d = scalar2(b1(1,k),vtemp1d(1))
 #ifdef MOMENT
       call matvec2(Ug2der(1,1,i+2),dd(1,1,itk1),vtemp2d(1))
       s8d = -(atemp(1,1)+atemp(2,2))*scalar2(cc(1,1,itl),vtemp2d(1))
@@ -8705,9 +8836,9 @@ C Derivatives in gamma(i+5)
       call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
       s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
 #endif
-      call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1d(1))
+      call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1d(1))
       call matvec2(AEAderg(1,1,1),vtemp1d(1),vtemp1d(1))
-      s2d = scalar2(b1(1,itk),vtemp1d(1))
+      s2d = scalar2(b1(1,k),vtemp1d(1))
 #ifdef MOMENT
       call transpose2(AEA(1,1,2),atempd(1,1))
       call matmat2(atempd(1,1),EUgder(1,1,i+4),atempd(1,1))
@@ -8717,7 +8848,7 @@ C Derivatives in gamma(i+5)
       s12d = scalar2(Ub2(1,i+2),vtemp3d(1))
 #ifdef MOMENT
       call matvec2(a_chuj(1,1,jj,i),Ub2der(1,i+4),vtemp4d(1)) 
-      ss13d = scalar2(b1(1,itk),vtemp4d(1))
+      ss13d = scalar2(b1(1,k),vtemp4d(1))
       s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
 #endif
 c      s1d=0.0d0
@@ -8741,10 +8872,10 @@ C Cartesian derivatives
             call matmat2(EUg(1,1,i+1),auxmatd(1,1),auxmatd(1,1))
             s1d = (auxmatd(1,1)+auxmatd(2,2))*ss1
 #endif
-            call matvec2(EUg(1,1,i+2),b1(1,itl),vtemp1(1))
+            call matvec2(EUg(1,1,i+2),b1(1,l),vtemp1(1))
             call matvec2(AEAderx(1,1,lll,kkk,iii,1),vtemp1(1),
      &          vtemp1d(1))
-            s2d = scalar2(b1(1,itk),vtemp1d(1))
+            s2d = scalar2(b1(1,k),vtemp1d(1))
 #ifdef MOMENT
             call transpose2(AEAderx(1,1,lll,kkk,iii,2),atempd(1,1))
             call matmat2(atempd(1,1),EUg(1,1,i+4),atempd(1,1))
@@ -8788,7 +8919,7 @@ c      s13d=0.0d0
           derx_turn(lll,kkk,2) = derx_turn(lll,kkk,2)-0.5d0*s13d
           call matvec2(a_chuj_der(1,1,lll,kkk,jj,i),Ub2(1,i+4),
      &      vtemp4d(1)) 
-          ss13d = scalar2(b1(1,itk),vtemp4d(1))
+          ss13d = scalar2(b1(1,k),vtemp4d(1))
           s13d = (gtemp(1,1)+gtemp(2,2))*ss13d
           derx_turn(lll,kkk,1) = derx_turn(lll,kkk,1)-0.5d0*s13d
         enddo
diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 6ea0840..9db9b60 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -907,29 +907,37 @@ C
         write (iout,*) "Coefficients of the cumulants"
       endif
       read (ifourier,*) nloctyp
+
       do i=0,nloctyp-1
         read (ifourier,*,end=115,err=115)
         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
+#ifdef NEWCORR
+        read (ifourier,*,end=115,err=115) (bnew1(ii,1,i),ii=1,3)
+        read (ifourier,*,end=115,err=115) (bnew2(ii,1,i),ii=1,3)
+        read (ifourier,*,end=115,err=115) (bnew1(ii,2,i),ii=1,1)
+        read (ifourier,*,end=115,err=115) (bnew2(ii,2,i),ii=1,1)
+        read (ifourier,*,end=115,err=115) (eenew(ii,i),ii=1,1)
+#endif 
         if (lprint) then
         write (iout,*) 'Type',i
         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
         endif
-        B1(1,i)  = b(3)
-        B1(2,i)  = b(5)
-        B1(1,-i) = b(3)
-        B1(2,-i) = -b(5)
+c        B1(1,i)  = b(3)
+c        B1(2,i)  = b(5)
+c        B1(1,-i) = b(3)
+c        B1(2,-i) = -b(5)
 c        b1(1,i)=0.0d0
 c        b1(2,i)=0.0d0
-        B1tilde(1,i) = b(3)
-        B1tilde(2,i) =-b(5)
-        B1tilde(1,-i) =-b(3)
-        B1tilde(2,-i) =b(5)
+c        B1tilde(1,i) = b(3)
+c        B1tilde(2,i) =-b(5)
+c        B1tilde(1,-i) =-b(3)
+c        B1tilde(2,-i) =b(5)
 c        b1tilde(1,i)=0.0d0
 c        b1tilde(2,i)=0.0d0
-        B2(1,i)  = b(2)
-        B2(2,i)  = b(4)
-        B2(1,-i)  =b(2)
-        B2(2,-i)  =-b(4)
+c        B2(1,i)  = b(2)
+c        B2(2,i)  = b(4)
+c        B2(1,-i)  =b(2)
+c        B2(2,-i)  =-b(4)
 
 c        b2(1,i)=0.0d0
 c        b2(2,i)=0.0d0
UNRESPACK 3.x