projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Zmiana 21 na ntyp1 w unres SRC_MD oraz SRC_MD-M
[unres.git] / source / unres / src_MD / energy_p_new_barrier.F
diff --git a/source/unres/src_MD/energy_p_new_barrier.F b/source/unres/src_MD/energy_p_new_barrier.F
index 23449da..cc43a72 100644 (file)
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -471,7 +471,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'mpif.h'
 #endif
       double precision gradbufc(3,maxres),gradbufx(3,maxres),
       include 'mpif.h'
 #endif
       double precision gradbufc(3,maxres),gradbufx(3,maxres),
-     &  glocbuf(4*maxres),gradbufc_sum(3,maxres)
+     &  glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.SETUP'
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
@@ -483,6 +483,7 @@ cMS$ATTRIBUTES C ::  proc_proc
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
       include 'COMMON.MAXGRAD'
       include 'COMMON.CONTROL'
       include 'COMMON.TIME1'
       include 'COMMON.MAXGRAD'
+      include 'COMMON.SCCOR'
 #ifdef TIMING
 #ifdef MPI
       time01=MPI_Wtime()
 #ifdef TIMING
 #ifdef MPI
       time01=MPI_Wtime()
@@ -755,7 +756,6 @@ c      enddo
      &   +wturn3*gel_loc_turn3(i)
      &   +wturn6*gel_loc_turn6(i)
      &   +wel_loc*gel_loc_loc(i)
      &   +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"
       enddo
 #ifdef DEBUG
       write (iout,*) "gloc after adding corr"
@@ -774,6 +774,21 @@ c      enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
+#define DEBUG
+#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
+#undef DEBUG
+        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
         time00=MPI_Wtime()
         call MPI_Barrier(FG_COMM,IERR)
         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
@@ -784,7 +799,19 @@ 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)
      &    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
         time_reduce=time_reduce+MPI_Wtime()-time00
+#define DEBUG
+#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
+#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
 #ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
@@ -1038,7 +1065,7 @@ C------------------------------------------------------------------------
      & 'ESC=   ',1pE16.6,' WEIGHT=',1pD16.6,' (SC local)'/
      & 'ETORS= ',1pE16.6,' WEIGHT=',1pD16.6,' (torsional)'/
      & 'ETORSD=',1pE16.6,' WEIGHT=',1pD16.6,' (double torsional)'/
      & '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,
+     & 'EHPB=  ',1pE16.6,' WEIGHT=',1pD16.6,
      & ' (SS bridges & dist. cnstr.)'/
      & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
      & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
      & ' (SS bridges & dist. cnstr.)'/
      & 'ECORR4=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
      & 'ECORR5=',1pE16.6,' WEIGHT=',1pD16.6,' (multi-body)'/
@@ -1110,8 +1137,8 @@ C
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1124,7 +1151,7 @@ C
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
           do j=istart(i,iint),iend(i,iint)
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
           do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -1287,8 +1314,8 @@ C
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1297,7 +1324,7 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -1404,8 +1431,8 @@ c     else
 c     endif
       ind=0
       do i=iatsc_s,iatsc_e
 c     endif
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1420,7 +1447,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=itype(j)
+            itypj=iabs(itype(j))
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             chi1=chi(itypi,itypj)
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             chi1=chi(itypi,itypj)
@@ -1540,8 +1567,8 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c     if (icall.eq.0) lprn=.false.
       ind=0
       do i=iatsc_s,iatsc_e
 c     if (icall.eq.0) lprn=.false.
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1558,7 +1585,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=itype(j)
+            itypj=iabs(itype(j))
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
@@ -1699,8 +1726,8 @@ c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
 c     if (icall.eq.0) lprn=.true.
       ind=0
       do i=iatsc_s,iatsc_e
 c     if (icall.eq.0) lprn=.true.
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1715,7 +1742,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=itype(j)
+            itypj=iabs(itype(j))
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             sig0ij=sigma(itypi,itypj)
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             sig0ij=sigma(itypi,itypj)
@@ -2019,8 +2046,8 @@ C
 cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
 cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
-        itypi=itype(i)
-        itypi1=itype(i+1)
+        itypi=iabs(itype(i))
+        itypi1=iabs(itype(i+1))
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -2031,7 +2058,7 @@ C
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
           do j=istart(i,iint),iend(i,iint)
 cd        write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
 cd   &                  'iend=',iend(i,iint)
           do j=istart(i,iint),iend(i,iint)
-            itypj=itype(j)
+            itypj=iabs(itype(j))
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -4032,7 +4059,7 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
-          itypj=itype(j)
+          itypj=iabs(itype(j))
 C Uncomment following three lines for SC-p interactions
 c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
 C Uncomment following three lines for SC-p interactions
 c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
@@ -4126,7 +4153,7 @@ cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
         do iint=1,nscp_gr(i)
 
         do j=iscpstart(i,iint),iscpend(i,iint)
-          itypj=itype(j)
+          itypj=iabs(itype(j))
 C Uncomment following three lines for SC-p interactions
 c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
 C Uncomment following three lines for SC-p interactions
 c         xj=c(1,nres+j)-xi
 c         yj=c(2,nres+j)-yi
@@ -4239,49 +4266,91 @@ C iii and jjj point to the residues for which the distance is assigned.
           iii=ii
           jjj=jj
         endif
           iii=ii
           jjj=jj
         endif
-cd        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
+c        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
+c     &    dhpb(i),dhpb1(i),forcon(i)
 C 24/11/03 AL: SS bridges handled separately because of introducing a specific
 C    distance and angle dependent SS bond potential.
 C 24/11/03 AL: SS bridges handled separately because of introducing a specific
 C    distance and angle dependent SS bond potential.
-        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. iabs(itype(jjj
+     &)).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
 cd          write (iout,*) "eij",eij
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
 cd          write (iout,*) "eij",eij
+        else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+          dd=dist(ii,jj)
+          if (dhpb1(i).gt.0.0d0) then
+            ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+            fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c            write (iout,*) "beta nmr",
+c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
+            dd=dist(ii,jj)
+            rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+            waga=forcon(i)
+C Calculate the contribution to energy.
+            ehpb=ehpb+waga*rdis*rdis
+c            write (iout,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+            fac=waga*rdis/dd
+          endif  
+          do j=1,3
+            ggg(j)=fac*(c(j,jj)-c(j,ii))
+          enddo
+          do j=1,3
+            ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+            ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+          enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
         else
 C Calculate the distance between the two points and its difference from the
 C target distance.
         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)
+          if (dhpb1(i).gt.0.0d0) then
+            ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+            fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c            write (iout,*) "alph nmr",
+c     &        dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+          else
+            rdis=dd-dhpb(i)
 C Get the force constant corresponding to this distance.
 C Get the force constant corresponding to this distance.
-        waga=forcon(i)
+            waga=forcon(i)
 C Calculate the contribution to energy.
 C Calculate the contribution to energy.
-        ehpb=ehpb+waga*rdis*rdis
+            ehpb=ehpb+waga*rdis*rdis
+c            write (iout,*) "alpha reg",dd,waga*rdis*rdis
 C
 C Evaluate gradient.
 C
 C
 C Evaluate gradient.
 C
-        fac=waga*rdis/dd
+            fac=waga*rdis/dd
+          endif
 cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
 cd   &   ' waga=',waga,' fac=',fac
 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).
 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
           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
 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
         endif
       enddo
       ehpb=0.5D0*ehpb
@@ -4306,7 +4375,7 @@ C
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
       double precision erij(3),dcosom1(3),dcosom2(3),gg(3)
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
       double precision erij(3),dcosom1(3),dcosom2(3),gg(3)
-      itypi=itype(i)
+      itypi=iabs(itype(i))
       xi=c(1,nres+i)
       yi=c(2,nres+i)
       zi=c(3,nres+i)
       xi=c(1,nres+i)
       yi=c(2,nres+i)
       zi=c(3,nres+i)
@@ -4315,7 +4384,7 @@ C
       dzi=dc_norm(3,nres+i)
 c      dsci_inv=dsc_inv(itypi)
       dsci_inv=vbld_inv(nres+i)
       dzi=dc_norm(3,nres+i)
 c      dsci_inv=dsc_inv(itypi)
       dsci_inv=vbld_inv(nres+i)
-      itypj=itype(j)
+      itypj=iabs(itype(j))
 c      dscj_inv=dsc_inv(itypj)
       dscj_inv=vbld_inv(nres+j)
       xj=c(1,nres+j)-xi
 c      dscj_inv=dsc_inv(itypj)
       dscj_inv=vbld_inv(nres+j)
       xj=c(1,nres+j)-xi
@@ -4409,7 +4478,7 @@ c
 c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
 c
       do i=ibond_start,ibond_end
 c 09/18/07 AL: multimodal bond potential based on AM1 CA-SC PMF's included
 c
       do i=ibond_start,ibond_end
-        iti=itype(i)
+        iti=iabs(itype(i))
         if (iti.ne.10) then
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
         if (iti.ne.10) then
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
@@ -4484,7 +4553,19 @@ c     write (*,'(a,i2)') 'EBEND ICG=',icg
       do i=ithet_start,ithet_end
 C Zero the energy function and its derivative at 0 or pi.
         call splinthet(theta(i),0.5d0*delta,ss,ssd)
       do i=ithet_start,ithet_end
 C Zero the energy function and its derivative at 0 or pi.
         call splinthet(theta(i),0.5d0*delta,ss,ssd)
-        it=itype(i-1)
+        it=(itype(i-1))
+        ichir1=isign(1,itype(i-2))
+        ichir2=isign(1,itype(i))
+        if (itype(i-2).eq.10) ichir1=isign(1,itype(i-1))
+        if (itype(i).eq.10) ichir2=isign(1,itype(i-1))
+        if (itype(i-1).eq.10) then
+         itype1=isign(10,itype(i-2))
+         ichir11=isign(1,itype(i-2))
+         ichir12=isign(1,itype(i-2))
+         itype2=isign(10,itype(i))
+         ichir21=isign(1,itype(i))
+         ichir22=isign(1,itype(i))
+        endif
         if (i.gt.3) then
 #ifdef OSF
          phii=phi(i)
         if (i.gt.3) then
 #ifdef OSF
          phii=phi(i)
@@ -4518,15 +4599,27 @@ C dependent on the adjacent virtual-bond-valence angles (gamma1 & gamma2).
 C In following comments this theta will be referred to as t_c.
         thet_pred_mean=0.0d0
         do k=1,2
 C In following comments this theta will be referred to as t_c.
         thet_pred_mean=0.0d0
         do k=1,2
-          athetk=athet(k,it)
-          bthetk=bthet(k,it)
+          athetk=athet(k,it,ichir1,ichir2)
+          bthetk=bthet(k,it,ichir1,ichir2)
+        if (it.eq.10) then
+           athetk=athet(k,itype1,ichir11,ichir12)
+           bthetk=bthet(k,itype2,ichir21,ichir22)
+        endif
           thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
         enddo
         dthett=thet_pred_mean*ssd
         thet_pred_mean=thet_pred_mean*ss+a0thet(it)
 C Derivatives of the "mean" values in gamma1 and gamma2.
           thet_pred_mean=thet_pred_mean+athetk*y(k)+bthetk*z(k)
         enddo
         dthett=thet_pred_mean*ssd
         thet_pred_mean=thet_pred_mean*ss+a0thet(it)
 C Derivatives of the "mean" values in gamma1 and gamma2.
-        dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
-        dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
+        dthetg1=(-athet(1,it,ichir1,ichir2)*y(2)
+     &+athet(2,it,ichir1,ichir2)*y(1))*ss
+        dthetg2=(-bthet(1,it,ichir1,ichir2)*z(2)
+     &          +bthet(2,it,ichir1,ichir2)*z(1))*ss
+        if (it.eq.10) then
+      dthetg1=(-athet(1,itype1,ichir11,ichir12)*y(2)
+     &+athet(2,itype1,ichir11,ichir12)*y(1))*ss
+        dthetg2=(-bthet(1,itype2,ichir21,ichir22)*z(2)
+     &         +bthet(2,itype2,ichir21,ichir22)*z(1))*ss
+        endif
         if (theta(i).gt.pi-delta) then
           call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
      &         E_tc0)
         if (theta(i).gt.pi-delta) then
           call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
      &         E_tc0)
@@ -4698,7 +4791,7 @@ C
         dephii=0.0d0
         dephii1=0.0d0
         theti2=0.5d0*theta(i)
         dephii=0.0d0
         dephii1=0.0d0
         theti2=0.5d0*theta(i)
-        ityp2=ithetyp(itype(i-1))
+        ityp2=ithetyp(iabs(itype(i-1)))
         do k=1,nntheterm
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         do k=1,nntheterm
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
@@ -4710,7 +4803,7 @@ C
 #else
           phii=phi(i)
 #endif
 #else
           phii=phi(i)
 #endif
-          ityp1=ithetyp(itype(i-2))
+          ityp1=ithetyp(iabs(itype(i-2)))
           do k=1,nsingle
             cosph1(k)=dcos(k*phii)
             sinph1(k)=dsin(k*phii)
           do k=1,nsingle
             cosph1(k)=dcos(k*phii)
             sinph1(k)=dsin(k*phii)
@@ -4731,7 +4824,7 @@ C
 #else
           phii1=phi(i+1)
 #endif
 #else
           phii1=phi(i+1)
 #endif
-          ityp3=ithetyp(itype(i))
+          ityp3=ithetyp(iabs(itype(i)))
           do k=1,nsingle
             cosph2(k)=dcos(k*phii1)
             sinph2(k)=dsin(k*phii1)
           do k=1,nsingle
             cosph2(k)=dcos(k*phii1)
             sinph2(k)=dsin(k*phii1)
@@ -4882,7 +4975,7 @@ c     write (iout,'(a)') 'ESC'
       do i=loc_start,loc_end
         it=itype(i)
         if (it.eq.10) goto 1
       do i=loc_start,loc_end
         it=itype(i)
         if (it.eq.10) goto 1
-        nlobit=nlob(it)
+        nlobit=nlob(iabs(it))
 c       print *,'i=',i,' it=',it,' nlobit=',nlobit
 c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
         theti=theta(i+1)-pipol
 c       print *,'i=',i,' it=',it,' nlobit=',nlobit
 c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
         theti=theta(i+1)-pipol
@@ -5039,11 +5132,11 @@ C Compute the contribution to SC energy and derivatives
 
           do j=1,nlobit
 #ifdef OSF
 
           do j=1,nlobit
 #ifdef OSF
-            adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin
+            adexp=bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin
             if(adexp.ne.adexp) adexp=1.0
             expfac=dexp(adexp)
 #else
             if(adexp.ne.adexp) adexp=1.0
             expfac=dexp(adexp)
 #else
-            expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
+            expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin)
 #endif
 cd          print *,'j=',j,' expfac=',expfac
             escloc_i=escloc_i+expfac
 #endif
 cd          print *,'j=',j,' expfac=',expfac
             escloc_i=escloc_i+expfac
@@ -5125,7 +5218,7 @@ C Compute the contribution to SC energy and derivatives
 
       dersc12=0.0d0
       do j=1,nlobit
 
       dersc12=0.0d0
       do j=1,nlobit
-        expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin)
+        expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin)
         escloc_i=escloc_i+expfac
         do k=1,2
           dersc(k)=dersc(k)+Ax(k,j)*expfac
         escloc_i=escloc_i+expfac
         do k=1,2
           dersc(k)=dersc(k)+Ax(k,j)*expfac
@@ -5227,7 +5320,7 @@ c
         do j = 1,3
           xx = xx + x_prime(j)*dc_norm(j,i+nres)
           yy = yy + y_prime(j)*dc_norm(j,i+nres)
         do j = 1,3
           xx = xx + x_prime(j)*dc_norm(j,i+nres)
           yy = yy + y_prime(j)*dc_norm(j,i+nres)
-          zz = zz + z_prime(j)*dc_norm(j,i+nres)
+          zz = zz + dsign(1.0,itype(i))*z_prime(j)*dc_norm(j,i+nres)
         enddo
 
         xxtab(i)=xx
         enddo
 
         xxtab(i)=xx
@@ -5237,7 +5330,7 @@ C
 C Compute the energy of the ith side cbain
 C
 c        write (2,*) "xx",xx," yy",yy," zz",zz
 C Compute the energy of the ith side cbain
 C
 c        write (2,*) "xx",xx," yy",yy," zz",zz
-        it=itype(i)
+        it=iabs(itype(i))
         do j = 1,65
           x(j) = sc_parmin(j,it) 
         enddo
         do j = 1,65
           x(j) = sc_parmin(j,it) 
         enddo
@@ -5245,7 +5338,7 @@ c        write (2,*) "xx",xx," yy",yy," zz",zz
 Cc diagnostics - remove later
         xx1 = dcos(alph(2))
         yy1 = dsin(alph(2))*dcos(omeg(2))
 Cc diagnostics - remove later
         xx1 = dcos(alph(2))
         yy1 = dsin(alph(2))*dcos(omeg(2))
-        zz1 = -dsin(alph(2))*dsin(omeg(2))
+        zz1 = -dsign(1.0,itype(i))*dsin(alph(2))*dsin(omeg(2))
         write(2,'(3f8.1,3f9.3,1x,3f9.3)') 
      &    alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
      &    xx1,yy1,zz1
         write(2,'(3f8.1,3f9.3,1x,3f9.3)') 
      &    alph(2)*rad2deg,omeg(2)*rad2deg,theta(3)*rad2deg,xx,yy,zz,
      &    xx1,yy1,zz1
@@ -5650,7 +5743,7 @@ C Proline-Proline pair is a special case...
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
      &  (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
         gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
      &  (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
         gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
-c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
+        write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       enddo
 ! 6/20/98 - dihedral angle constraints
       edihcnstr=0.0d0
       enddo
 ! 6/20/98 - dihedral angle constraints
       edihcnstr=0.0d0
@@ -5704,12 +5797,17 @@ c     lprn=.true.
       etors_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
       etors_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
+        if (iabs(itype(i)).eq.20) then
+        iblock=2
+        else
+        iblock=1
+        endif
         phii=phi(i)
         gloci=0.0D0
 C Regular cosine and sine terms
         phii=phi(i)
         gloci=0.0D0
 C Regular cosine and sine terms
-        do j=1,nterm(itori,itori1)
-          v1ij=v1(j,itori,itori1)
-          v2ij=v2(j,itori,itori1)
+        do j=1,nterm(itori,itori1,iblock)
+          v1ij=v1(j,itori,itori1,iblock)
+          v2ij=v2(j,itori,itori1,iblock)
           cosphi=dcos(j*phii)
           sinphi=dsin(j*phii)
           etors=etors+v1ij*cosphi+v2ij*sinphi
           cosphi=dcos(j*phii)
           sinphi=dsin(j*phii)
           etors=etors+v1ij*cosphi+v2ij*sinphi
@@ -5724,7 +5822,7 @@ C          [v2 cos(phi/2)+v3 sin(phi/2)]^2 + 1
 C
         cosphi=dcos(0.5d0*phii)
         sinphi=dsin(0.5d0*phii)
 C
         cosphi=dcos(0.5d0*phii)
         sinphi=dsin(0.5d0*phii)
-        do j=1,nlor(itori,itori1)
+        do j=1,nlor(itori,itori1,iblock)
           vl1ij=vlor1(j,itori,itori1)
           vl2ij=vlor2(j,itori,itori1)
           vl3ij=vlor3(j,itori,itori1)
           vl1ij=vlor1(j,itori,itori1)
           vl2ij=vlor2(j,itori,itori1)
           vl3ij=vlor3(j,itori,itori1)
@@ -5737,13 +5835,14 @@ C
           gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
         enddo
 C Subtract the constant term
           gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
         enddo
 C Subtract the constant term
-        etors=etors-v0(itori,itori1)
+        etors=etors-v0(itori,itori1,iblock)
           if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
           if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
-     &         'etor',i,etors_ii-v0(itori,itori1)
+     &         'etor',i,etors_ii-v0(itori,itori1,iblock)
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
-     &  (v1(j,itori,itori1),j=1,6),(v2(j,itori,itori1),j=1,6)
+     &  (v1(j,itori,itori1,iblock),j=1,6),
+     &  (v2(j,itori,itori1,iblock),j=1,6)
         gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
 c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       enddo
         gloc(i-3,icg)=gloc(i-3,icg)+wtor*gloci
 c       write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg)
       enddo
@@ -5765,6 +5864,7 @@ c      do i=1,ndih_constr
         else
           difi=0.0
         endif
         else
           difi=0.0
         endif
+c        write (iout,*) "gloci", gloc(i-3,icg)
 cd        write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
 cd     &    rad2deg*phi0(i),  rad2deg*drange(i),
 cd     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
 cd        write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii,
 cd     &    rad2deg*phi0(i),  rad2deg*drange(i),
 cd     &    rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg)
@@ -5797,15 +5897,17 @@ c     lprn=.true.
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
+        iblock=1
+        if (iabs(itype(i+1)).eq.20) iblock=2
         phii=phi(i)
         phii1=phi(i+1)
         gloci1=0.0D0
         gloci2=0.0D0
         phii=phi(i)
         phii1=phi(i+1)
         gloci1=0.0D0
         gloci2=0.0D0
-        do j=1,ntermd_1(itori,itori1,itori2)
-          v1cij=v1c(1,j,itori,itori1,itori2)
-          v1sij=v1s(1,j,itori,itori1,itori2)
-          v2cij=v1c(2,j,itori,itori1,itori2)
-          v2sij=v1s(2,j,itori,itori1,itori2)
+        do j=1,ntermd_1(itori,itori1,itori2,iblock)
+          v1cij=v1c(1,j,itori,itori1,itori2,iblock)
+          v1sij=v1s(1,j,itori,itori1,itori2,iblock)
+          v2cij=v1c(2,j,itori,itori1,itori2,iblock)
+          v2sij=v1s(2,j,itori,itori1,itori2,iblock)
           cosphi1=dcos(j*phii)
           sinphi1=dsin(j*phii)
           cosphi2=dcos(j*phii1)
           cosphi1=dcos(j*phii)
           sinphi1=dsin(j*phii)
           cosphi2=dcos(j*phii1)
@@ -5815,12 +5917,12 @@ c     lprn=.true.
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
-        do k=2,ntermd_2(itori,itori1,itori2)
+        do k=2,ntermd_2(itori,itori1,itori2,iblock)
           do l=1,k-1
           do l=1,k-1
-            v1cdij = v2c(k,l,itori,itori1,itori2)
-            v2cdij = v2c(l,k,itori,itori1,itori2)
-            v1sdij = v2s(k,l,itori,itori1,itori2)
-            v2sdij = v2s(l,k,itori,itori1,itori2)
+            v1cdij = v2c(k,l,itori,itori1,itori2,iblock)
+            v2cdij = v2c(l,k,itori,itori1,itori2,iblock)
+            v1sdij = v2s(k,l,itori,itori1,itori2,iblock)
+            v2sdij = v2s(l,k,itori,itori1,itori2,iblock)
             cosphi1p2=dcos(l*phii+(k-l)*phii1)
             cosphi1m2=dcos(l*phii-(k-l)*phii1)
             sinphi1p2=dsin(l*phii+(k-l)*phii1)
             cosphi1p2=dcos(l*phii+(k-l)*phii1)
             cosphi1m2=dcos(l*phii-(k-l)*phii1)
             sinphi1p2=dsin(l*phii+(k-l)*phii1)
@@ -5835,6 +5937,7 @@ c     lprn=.true.
         enddo
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
         enddo
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
+c        write (iout,*) "gloci", gloc(i-3,icg)
       enddo
       return
       end
       enddo
       return
       end
@@ -5867,7 +5970,7 @@ C Set lprn=.true. for debugging
 c      lprn=.true.
 c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
 c      lprn=.true.
 c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
-      do i=iphi_start,iphi_end
+      do i=itau_start,itau_end
         esccor_ii=0.0D0
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
         esccor_ii=0.0D0
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
@@ -5878,17 +5981,24 @@ c(see comment below)
 cc Omicron is flat angle depending on the value of first digit 
 c(see comment below)
 
 cc Omicron is flat angle depending on the value of first digit 
 c(see comment below)
 
-        gloci=0.0D0
-        do intertyp=1,3
+        
+        do intertyp=1,3 !intertyp
 cc Added 09 May 2012 (Adasko)
 cc  Intertyp means interaction type of backbone mainchain correlation: 
 c   1 = SC...Ca...Ca...Ca
 c   2 = Ca...Ca...Ca...SC
 cc Added 09 May 2012 (Adasko)
 cc  Intertyp means interaction type of backbone mainchain correlation: 
 c   1 = SC...Ca...Ca...Ca
 c   2 = Ca...Ca...Ca...SC
-c   3 = SC...Ca...Ca...SC
-        if (((intertyp.eq.3).and.(itype(i-2).eq.10).or.
-     &      (itype(i-1).eq.10))
-     &    .or. ((intertyp.eq.1).and.(itype(i-2).ne.10))
-     &    .or. ((intertyp.eq.2).and.(itype(i-1).ne.10))) cycle  
+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)))
+     &    .or. ((intertyp.eq.1).and.((itype(i-2).eq.10)
+     &     .or.(itype(i-2).eq.ntyp1)))
+     &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
+     &      (itype(i-1).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))
+     & cycle
         do j=1,nterm_sccor(isccori,isccori1)
           v1ij=v1sccor(j,intertyp,isccori,isccori1)
           v2ij=v2sccor(j,intertyp,isccori,isccori1)
         do j=1,nterm_sccor(isccori,isccori1)
           v1ij=v1sccor(j,intertyp,isccori,isccori1)
           v2ij=v2sccor(j,intertyp,isccori,isccori1)
@@ -5897,15 +6007,20 @@ c   3 = SC...Ca...Ca...SC
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
-        gloc_sc(intertyp,i-3,icg)=gloc_sc(i-3,icg)+wtor*gloci
+        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,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
         if (lprn)
      &  write (iout,'(2(a3,2x,i3,2x),2i3,6f8.3/26x,6f8.3/)')
      &  restyp(itype(i-2)),i-2,restyp(itype(i-1)),i-1,itori,itori1,
      &  (v1sccor(j,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
       enddo
-      enddo
+c        do i=1,nres
+c        write (iout,*) "W@T@F",  gloc_sc(1,i,icg),gloc(i,icg)
+c        enddo
       return
       end
 c----------------------------------------------------------------------------
       return
       end
 c----------------------------------------------------------------------------
@@ -8012,9 +8127,9 @@ C
 C      Parallel       Antiparallel
 C                                             
 C          o             o         
 C      Parallel       Antiparallel
 C                                             
 C          o             o         
-C         /l\           /j\       
-C        /   \         /   \      
-C       /| o |         | o |\     
+C         /l\           /j\
+C        /   \         /   \
+C       /| o |         | o |\
 C     \ j|/k\|  /   \  |/k\|l /   
 C      \ /   \ /     \ /   \ /    
 C       o     o       o     o                
 C     \ j|/k\|  /   \  |/k\|l /   
 C      \ /   \ /     \ /   \ /    
 C       o     o       o     o                
@@ -8113,18 +8228,18 @@ c----------------------------------------------------------------------------
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       logical lprn
       common /kutas/ lprn
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C     \   /l\           /j\   /   
-C      \ /   \         /   \ /    
-C       o| o |         | o |o     
-C     \ j|/k\|      \  |/k\|l     
-C      \ /   \       \ /   \      
-C       o             o                      
-C       i             i                     
-C
+C                                                                              C
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C
+C     \   /l\           /j\   /                                                C
+C      \ /   \         /   \ /                                                 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           
 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, 
 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, 
@@ -8295,18 +8410,18 @@ c----------------------------------------------------------------------------
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C         /l\   /   \   /j\       
-C        /   \ /     \ /   \      
-C       /| o |o       o| o |\     
-C       j|/k\|  /      |/k\|l /   
-C        /   \ /       /   \ /    
-C       /     o       /     o                
-C       i             i                     
-C
+C                                                                              C 
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C 
+C         /l\   /   \   /j\                                                    C 
+C        /   \ /     \ /   \                                                   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
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
@@ -8412,18 +8527,18 @@ c----------------------------------------------------------------------------
      & auxvec1(2),auxmat1(2,2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      & auxvec1(2),auxmat1(2,2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                              
-C      Parallel       Antiparallel
-C                                             
-C          o             o         
-C         /l\   /   \   /j\       
-C        /   \ /     \ /   \      
-C       /| o |o       o| o |\     
-C     \ j|/k\|      \  |/k\|l     
-C      \ /   \       \ /   \      
-C       o     \       o     \                
-C       i             i                     
-C
+C                                                                              C                       
+C      Parallel       Antiparallel                                             C
+C                                                                              C
+C          o             o                                                     C
+C         /l\   /   \   /j\                                                    C
+C        /   \ /     \ /   \                                                   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 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective 
UNRESPACK 3.x