projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Commit Adam 6/29/2014
[unres.git] / source / wham / src-NEWSC-NEWCORR / energy_p_new.F
diff --git a/source/wham/src-NEWSC-NEWCORR/energy_p_new.F b/source/wham/src-NEWSC-NEWCORR/energy_p_new.F
index 0ee066f..b847313 100644 (file)
--- a/source/wham/src-NEWSC-NEWCORR/energy_p_new.F
+++ b/source/wham/src-NEWSC-NEWCORR/energy_p_new.F
@@ -3406,11 +3406,11 @@ C
           iti1=ntortyp+1
         endif
 cd        write (iout,*) '*******i',i,' iti1',iti
-cd        write (iout,*) 'b1',b1(:,iti)
-cd        write (iout,*) 'b2',b2(:,iti)
+cd        write (iout,*) 'b1',b1(:,i-2)
+cd        write (iout,*) 'b2',b2(:,i-2)
 cd        write (iout,*) 'Ug',Ug(:,:,i-2)
         if (i .gt. iatel_s+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))
           call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2))
           call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2))
           call matmat2(DD(1,1,iti),Ug(1,1,i-2),DUg(1,1,i-2))
@@ -3430,7 +3430,7 @@ cd        write (iout,*) 'Ug',Ug(:,:,i-2)
             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))
         call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2))
         call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2))
@@ -3446,10 +3446,10 @@ cd        write (iout,*) 'Ug',Ug(:,:,i-2)
           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
 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))
@@ -4346,10 +4346,10 @@ cd        call checkint_turn4(i,a_temp,eello_turn4_num)
         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))
+        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))
@@ -4361,14 +4361,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))
@@ -4376,10 +4376,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),auxmat(1,1))
         call matmat2(auxmat(1,1),e1t(1,1),pizda(1,1))
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
@@ -4394,10 +4394,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))
@@ -4413,10 +4413,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))
@@ -4427,10 +4427,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))
@@ -4441,10 +4441,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))
@@ -4455,10 +4455,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))
@@ -7013,10 +7013,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
@@ -7195,26 +7195,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))
@@ -7223,20 +7223,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))
@@ -7333,26 +7333,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))
@@ -7361,20 +7361,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))
@@ -7679,7 +7679,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))
       if (calc_grad) then
 C Explicit gradient in virtual-dihedral angles.
@@ -7690,11 +7690,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
@@ -7706,7 +7706,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
@@ -7764,7 +7764,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))
         if (calc_grad) then
 C Explicit gradient in virtual-dihedral angles.
@@ -7774,7 +7774,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
@@ -7785,7 +7785,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
@@ -7841,7 +7841,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))
         if (calc_grad) then
 C Explicit gradient in virtual-dihedral angles.
@@ -7851,7 +7851,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
@@ -7862,7 +7862,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
@@ -8127,8 +8127,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)
@@ -8142,8 +8142,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))
@@ -8182,10 +8182,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
@@ -8428,10 +8428,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)
@@ -8446,13 +8446,13 @@ cd      write (2,*) 'eello6_graph3:','s1',s1,' s2',s2,' s3',s3,' s4',s4
 c      eello6_graph3=-s4
       if (.not. calc_grad) return
 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)
@@ -8469,12 +8469,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)
@@ -8563,11 +8563,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))
@@ -8592,11 +8592,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
@@ -8625,11 +8625,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,tl),auxvec1(1))
       endif
       call transpose2(EUgder(1,1,k),auxmat1(1,1))
       call matmat2(AECA(1,1,imat),auxmat1(1,1),pizda(1,1))
@@ -8695,12 +8695,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))
@@ -8797,14 +8797,14 @@ 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
 #else
       s1 = 0.0d0
 #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))
@@ -8821,7 +8821,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
 #else
       s13=0.0d0
@@ -8858,14 +8858,14 @@ 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
 #else
       s1d=0.0d0
 #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))
@@ -8919,9 +8919,9 @@ C Derivatives in gamma(i+5)
 #else
       s1d = 0.0d0
 #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))
@@ -8933,7 +8933,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
 #else
       s13d = 0.0d0
@@ -8961,10 +8961,10 @@ C Cartesian derivatives
 #else
             s1d = 0.0d0
 #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))
@@ -9010,7 +9010,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
UNRESPACK 3.x