Merge branch 'prerelease-3.2.1'
[unres.git] / source / unres / src_MIN / energy_p_new_barrier.F
index 3f19648..740912a 100644 (file)
@@ -441,12 +441,9 @@ cMS$ATTRIBUTES C ::  proc_proc
 #endif
 #ifdef MPI
       include 'mpif.h'
-      double precision gradbufc(3,maxres),gradbufx(3,maxres),
-     &  glocbuf(4*maxres),gradbufc_sum(3,maxres)
-#else
-      double precision gradbufc(3,maxres),gradbufx(3,maxres),
-     &  glocbuf(4*maxres),gradbufc_sum(3,maxres)
 #endif
+      double precision gradbufc(3,maxres),gradbufx(3,maxres),
+     &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
@@ -458,8 +455,13 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
       include 'COMMON.MAXGRAD'
+      include 'COMMON.SCCOR'
 #ifdef TIMING
+#ifdef MPI
       time01=MPI_Wtime()
+#else
+      time01=tcpu()
+#endif
 #endif
 #ifdef DEBUG
       write (iout,*) "sum_gradient gvdwc, gvdwx"
@@ -554,31 +556,61 @@ c      enddo
       enddo
       call flush(iout)
 #endif
-      call MPI_AllReduce(gradbufc(1,1),gradbufc_sum(1,1),3*nres,
-     &    MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR)
-      time_reduce=time_reduce+MPI_Wtime()-time00
-#ifdef DEBUG
-      write (iout,*) "gradbufc_sum after allreduce"
       do i=1,nres
-        write (iout,'(i3,3f10.5)') i,(gradbufc_sum(j,i),j=1,3)
+        do j=1,3
+          gradbufc_sum(j,i)=gradbufc(j,i)
+        enddo
       enddo
-      call flush(iout)
+c      call MPI_AllReduce(gradbufc(1,1),gradbufc_sum(1,1),3*nres,
+c     &    MPI_DOUBLE_PRECISION,MPI_SUM,FG_COMM,IERR)
+c      time_reduce=time_reduce+MPI_Wtime()-time00
+#ifdef DEBUG
+c      write (iout,*) "gradbufc_sum after allreduce"
+c      do i=1,nres
+c        write (iout,'(i3,3f10.5)') i,(gradbufc_sum(j,i),j=1,3)
+c      enddo
+c      call flush(iout)
 #endif
 #ifdef TIMING
-      time_allreduce=time_allreduce+MPI_Wtime()-time00
+c      time_allreduce=time_allreduce+MPI_Wtime()-time00
 #endif
       do i=nnt,nres
         do k=1,3
           gradbufc(k,i)=0.0d0
         enddo
       enddo
-      do i=igrad_start,igrad_end
-        do j=jgrad_start(i),jgrad_end(i)
-          do k=1,3
-            gradbufc(k,i)=gradbufc(k,i)+gradbufc_sum(k,j)
-          enddo
+#ifdef DEBUG
+      write (iout,*) "igrad_start",igrad_start," igrad_end",igrad_end
+      write (iout,*) (i," jgrad_start",jgrad_start(i),
+     &                  " jgrad_end  ",jgrad_end(i),
+     &                  i=igrad_start,igrad_end)
+#endif
+c
+c Obsolete and inefficient code; we can make the effort O(n) and, therefore,
+c do not parallelize this part.
+c
+c      do i=igrad_start,igrad_end
+c        do j=jgrad_start(i),jgrad_end(i)
+c          do k=1,3
+c            gradbufc(k,i)=gradbufc(k,i)+gradbufc_sum(k,j)
+c          enddo
+c        enddo
+c      enddo
+      do j=1,3
+        gradbufc(j,nres-1)=gradbufc_sum(j,nres)
+      enddo
+      do i=nres-2,nnt,-1
+        do j=1,3
+          gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
         enddo
       enddo
+#ifdef DEBUG
+      write (iout,*) "gradbufc after summing"
+      do i=1,nres
+        write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
+      enddo
+      call flush(iout)
+#endif
       else
 #endif
 #ifdef DEBUG
@@ -588,16 +620,37 @@ c      enddo
       enddo
       call flush(iout)
 #endif
-      do i=nnt,nres-1
-        do k=1,3
-          gradbufc(k,i)=0.0d0
+      do i=1,nres
+        do j=1,3
+          gradbufc_sum(j,i)=gradbufc(j,i)
+          gradbufc(j,i)=0.0d0
         enddo
-        do j=i+1,nres
-          do k=1,3
-            gradbufc(k,i)=gradbufc(k,i)+gradbufc(k,j)
-          enddo
+      enddo
+      do j=1,3
+        gradbufc(j,nres-1)=gradbufc_sum(j,nres)
+      enddo
+      do i=nres-2,nnt,-1
+        do j=1,3
+          gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
         enddo
       enddo
+c      do i=nnt,nres-1
+c        do k=1,3
+c          gradbufc(k,i)=0.0d0
+c        enddo
+c        do j=i+1,nres
+c          do k=1,3
+c            gradbufc(k,i)=gradbufc(k,i)+gradbufc(k,j)
+c          enddo
+c        enddo
+c      enddo
+#ifdef DEBUG
+      write (iout,*) "gradbufc after summing"
+      do i=1,nres
+        write (iout,'(i3,3f10.5)') i,(gradbufc(j,i),j=1,3)
+      enddo
+      call flush(iout)
+#endif
 #ifdef MPI
       endif
 #endif
@@ -675,7 +728,6 @@ c      enddo
      &   +wturn3*gel_loc_turn3(i)
      &   +wturn6*gel_loc_turn6(i)
      &   +wel_loc*gel_loc_loc(i)
-     &   +wsccor*gsccor_loc(i)
       enddo
 #ifdef DEBUG
       write (iout,*) "gloc after adding corr"
@@ -694,6 +746,19 @@ c      enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
+#ifdef DEBUG
+      write (iout,*) "gloc_sc before reduce"
+      do i=1,nres
+       do j=1,3
+        write (iout,*) i,j,gloc_sc(j,i,icg)
+       enddo
+      enddo
+#endif
+        do i=1,nres
+         do j=1,3
+          gloc_scbuf(j,i)=gloc_sc(j,i,icg)
+         enddo
+        enddo
         time00=MPI_Wtime()
         call MPI_Barrier(FG_COMM,IERR)
         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
@@ -704,8 +769,18 @@ c      enddo
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         call MPI_Reduce(glocbuf(1),gloc(1,icg),4*nres,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
+        call MPI_Reduce(gloc_scbuf(1,1),gloc_sc(1,1,icg),3*nres,
+     &    MPI_DOUBLE_PRECISION,MPI_SUM,king,FG_COMM,IERR)
         time_reduce=time_reduce+MPI_Wtime()-time00
 #ifdef DEBUG
+      write (iout,*) "gloc_sc after reduce"
+      do i=1,nres
+       do j=1,3
+        write (iout,*) i,j,gloc_sc(j,i,icg)
+       enddo
+      enddo
+#endif
+#ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
         write (iout,*) i,gloc(i,icg)
@@ -831,7 +906,11 @@ c
       enddo 
 #endif
 #ifdef TIMING
+#ifdef MPI
       time_sumgradient=time_sumgradient+MPI_Wtime()-time01
+#else
+      time_sumgradient=time_sumgradient+tcpu()-time01
+#endif
 #endif
       return
       end
@@ -945,25 +1024,25 @@ C------------------------------------------------------------------------
      &  edihcnstr,ebr*nss,
      &  Uconst,etot
    10 format (/'Virtual-chain energies:'//
-     & 'EVDW=  ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
-     & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
-     & 'EES=   ',1pE16.6,' WEIGHT=',1pD16.6,' (p-p)'/
-     & 'EVDWPP=',1pE16.6,' WEIGHT=',1pD16.6,' (p-p VDW)'/
-     & 'ESTR=  ',1pE16.6,' WEIGHT=',1pD16.6,' (stretching)'/
-     & 'EBE=   ',1pE16.6,' WEIGHT=',1pD16.6,' (bending)'/
-     & 'ESC=   ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
-     & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
-     & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
-     & 'EHBP=  ',1pE16.6,' WEIGHT=',1pD16.6,
+     & 'EVDW=  ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-SC)'/
+     & 'EVDW2= ',1pE16.6,' WEIGHT=',1pE16.6,' (SC-p)'/
+     & 'EES=   ',1pE16.6,' WEIGHT=',1pE16.6,' (p-p)'/
+     & 'EVDWPP=',1pE16.6,' WEIGHT=',1pE16.6,' (p-p VDW)'/
+     & 'ESTR=  ',1pE16.6,' WEIGHT=',1pE16.6,' (stretching)'/
+     & 'EBE=   ',1pE16.6,' WEIGHT=',1pE16.6,' (bending)'/
+     & 'ESC=   ',1pE16.6,' WEIGHT=',1pE16.6,' (SC local)'/
+     & 'ETORS= ',1pE16.6,' WEIGHT=',1pE16.6,' (torsional)'/
+     & 'ETORSD=',1pE16.6,' WEIGHT=',1pE16.6,' (double torsional)'/
+     & 'EHPB=  ',1pE16.6,' WEIGHT=',1pE16.6,
      & ' (SS bridges & dist. cnstr.)'/
-     & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-     & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-     & 'ECORR6=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
-     & 'EELLO= ',1pE16.6,' WEIGHT=',1pD16.6,' (electrostatic-local)'/
-     & 'ETURN3=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 3rd order)'/
-     & 'ETURN4=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 4th order)'/
-     & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/
-     & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/
+     & 'ECORR4=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+     & 'ECORR5=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+     & 'ECORR6=',1pE16.6,' WEIGHT=',1pE16.6,' (multi-body)'/
+     & 'EELLO= ',1pE16.6,' WEIGHT=',1pE16.6,' (electrostatic-local)'/
+     & 'ETURN3=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 3rd order)'/
+     & 'ETURN4=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 4th order)'/
+     & 'ETURN6=',1pE16.6,' WEIGHT=',1pE16.6,' (turns, 6th order)'/
+     & 'ESCCOR=',1pE16.6,' WEIGHT=',1pE16.6,' (backbone-rotamer corr)'/
      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
@@ -4165,39 +4244,39 @@ cd          write (iout,*) "eij",eij
         else
 C Calculate the distance between the two points and its difference from the
 C target distance.
-        dd=dist(ii,jj)
-        rdis=dd-dhpb(i)
+          dd=dist(ii,jj)
+            rdis=dd-dhpb(i)
 C Get the force constant corresponding to this distance.
-        waga=forcon(i)
+            waga=forcon(i)
 C Calculate the contribution to energy.
-        ehpb=ehpb+waga*rdis*rdis
+            ehpb=ehpb+waga*rdis*rdis
 C
 C Evaluate gradient.
 C
-        fac=waga*rdis/dd
+            fac=waga*rdis/dd
 cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
 cd   &   ' waga=',waga,' fac=',fac
-        do j=1,3
-          ggg(j)=fac*(c(j,jj)-c(j,ii))
-        enddo
+            do j=1,3
+              ggg(j)=fac*(c(j,jj)-c(j,ii))
+            enddo
 cd      print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
 C If this is a SC-SC distance, we need to calculate the contributions to the
 C Cartesian gradient in the SC vectors (ghpbx).
-        if (iii.lt.ii) then
+          if (iii.lt.ii) then
           do j=1,3
             ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
             ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
           enddo
-        endif
+          endif
 cgrad        do j=iii,jjj-1
 cgrad          do k=1,3
 cgrad            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
 cgrad          enddo
 cgrad        enddo
-        do k=1,3
-          ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
-          ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
-        enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
         endif
       enddo
       ehpb=0.5D0*ehpb
@@ -4314,6 +4393,8 @@ c
       do i=ibondp_start,ibondp_end
         diff = vbld(i)-vbldp0
 c        write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+        if (energy_dec)    write (iout,'(a7,i5,4f7.3)') 
+     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
         do j=1,3
           gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
@@ -4332,6 +4413,12 @@ c
             diff=vbld(i+nres)-vbldsc0(1,iti)
 c            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
 c     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
+            if (energy_dec)  then
+              write (iout,*) 
+     &         "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+     &         AKSC(1,iti),AKSC(1,iti)*diff*diff
+              call flush(iout)
+            endif
             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
             do j=1,3
               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
@@ -4610,6 +4697,8 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
+        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &(itype(i).eq.ntyp1)) cycle
         dethetai=0.0d0
         dephii=0.0d0
         dephii1=0.0d0
@@ -4619,7 +4708,8 @@ C
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3) then
+C        if (i.gt.3) then
+        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+          ityp1=ithetyp(itype(i-2))
           do k=1,nsingle
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres) then
+        if ((i.lt.nres).and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -4654,7 +4744,7 @@ C
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -4764,6 +4854,8 @@ C
      &   i,theta(i)*rad2deg,phii*rad2deg,
      &   phii1*rad2deg,ethetai
         etheta=etheta+ethetai
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &      'ebend',i,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
         gloc(nphi+i-2,icg)=wang*dethetai
@@ -5704,12 +5796,14 @@ C 6/23/01 Compute double torsional energy
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
       logical lprn
 C Set lprn=.true. for debugging
       lprn=.false.
 c     lprn=.true.
       etors_d=0.0D0
       do i=iphid_start,iphid_end
+        etors_d_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
@@ -5729,6 +5823,8 @@ C Regular cosine and sine terms
           sinphi2=dsin(j*phii1)
           etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
      &     v2cij*cosphi2+v2sij*sinphi2
+          if (energy_dec) etors_d_ii=etors_d_ii+
+     &     v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
@@ -5744,12 +5840,17 @@ C Regular cosine and sine terms
             sinphi1m2=dsin(l*phii-(k-l)*phii1)
             etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
      &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
+            if (energy_dec) etors_d_ii=etors_d_ii+
+     &        v1cdij*cosphi1p2+v2cdij*cosphi1m2+
+     &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
             gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
             gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2+v2cdij*sinphi1m2) 
           enddo
         enddo
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &        'etor_d',i,etors_d_ii
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
       enddo
@@ -5787,10 +5888,19 @@ c      write (iout,*) "EBACK_SC_COR",itau_start,itau_end
       esccor=0.0D0
       do i=itau_start,itau_end
         esccor_ii=0.0D0
+        if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
 c      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
         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)
+C        print *,i,tauangle(1,i)
+        
         do intertyp=1,3 !intertyp
 cc Added 09 May 2012 (Adasko)
 cc  Intertyp means interaction type of backbone mainchain correlation: 
@@ -5799,18 +5909,16 @@ c   2 = Ca...Ca...Ca...SC
 c   3 = SC...Ca...Ca...SCi
         gloci=0.0D0
         if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
-     &      (itype(i-1).eq.10).or.(itype(i-2).eq.ntyp1).or.
-     &      (itype(i-1).eq.ntyp1)))
+     &      (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
+     &      (itype(i-1).eq.21)))
      &    .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
-     &     .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)
-     &     .or.(itype(i).eq.ntyp1)))
+     &     .or.(itype(i-2).eq.21)))
      &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
-     &      (itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
-     &      (itype(i-3).eq.ntyp1)))) cycle
-        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.ntyp1)) cycle
-        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.ntyp1))
+     &      (itype(i-1).eq.21)))) cycle  
+        if ((intertyp.eq.2).and.(i.eq.4).and.(itype(1).eq.21)) cycle
+        if ((intertyp.eq.1).and.(i.eq.nres).and.(itype(nres).eq.21))
      & cycle
-       do j=1,nterm_sccor(isccori,isccori1)
+        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))
@@ -5820,14 +5928,20 @@ c   3 = SC...Ca...Ca...SCi
         enddo
 c      write (iout,*) "EBACK_SC_COR",i,esccor,intertyp
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
+c        write (iout,*) "WTF",intertyp,i,itype(i),v1ij*cosphi+v2ij*sinphi
+c     &gloc_sc(intertyp,i-3,icg)
         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,isccori,isccori1,
-     &  (v1sccor(j,intertyp,isccori,isccori1),j=1,6)
-     & ,(v2sccor(j,intertyp,isccori,isccori1),j=1,6)
-C        gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
+     &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
+     &  (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 !intertyp
       enddo
+c        do i=1,nres
+c        write (iout,*) "W@T@F",  gloc_sc(1,i,icg),gloc_sc(2,i,icg),
+c     &   gloc_sc(3,i,icg)
+c        enddo
       return
       end
 c----------------------------------------------------------------------------
@@ -8041,12 +8155,12 @@ C                                                                              C
 C          o             o                                                     C
 C     \   /l\           /j\   /                                                C
 C      \ /   \         /   \ /                                                 C
-C       o| o |         | o |o                                                  C
+C       o| o |         | o |o                                                  C                
 C     \ j|/k\|      \  |/k\|l                                                  C
 C      \ /   \       \ /   \                                                   C
 C       o             o                                                        C
-C       i             i                                                        C
-C                                                                              C
+C       i             i                                                        C 
+C                                                                              C           
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 cd      write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
 C AL 7/4/01 s1 would occur in the sixth-order moment, 
@@ -8217,11 +8331,11 @@ c----------------------------------------------------------------------------
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                                                              C
+C                                                                              C 
 C      Parallel       Antiparallel                                             C
 C                                                                              C
 C          o             o                                                     C 
-C         /l\   /   \   /j\                                                    C
+C         /l\   /   \   /j\                                                    C 
 C        /   \ /     \ /   \                                                   C
 C       /| o |o       o| o |\                                                  C
 C       j|/k\|  /      |/k\|l /                                                C
@@ -8334,18 +8448,18 @@ c----------------------------------------------------------------------------
      & auxvec1(2),auxmat1(2,2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                                                              C
+C                                                                              C                       
 C      Parallel       Antiparallel                                             C
 C                                                                              C
-C          o             o                                                     C 
+C          o             o                                                     C
 C         /l\   /   \   /j\                                                    C
 C        /   \ /     \ /   \                                                   C
-C       /| o |o       o| o |\                                                  C 
+C       /| o |o       o| o |\                                                  C
 C     \ j|/k\|      \  |/k\|l                                                  C
-C      \ /   \       \ /   \                                                   C
+C      \ /   \       \ /   \                                                   C 
 C       o     \       o     \                                                  C
 C       i             i                                                        C
-C                                                                              C
+C                                                                              C 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective