projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Merge branch 'master' of mmka:unres into multichain
[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 3979f64..33b9ff7 100644 (file)
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -131,6 +131,11 @@ C
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
 C Calculate electrostatic (H-bonding) energy of the main chain.
 C
   107 continue
+cmc
+cmc Sep-06: egb takes care of dynamic ss bonds too
+cmc
+c      if (dyn_ss) call dyn_set_nss
+
 c      print *,"Processor",myrank," computed USCSC"
 #ifdef TIMING
 #ifdef MPI
 c      print *,"Processor",myrank," computed USCSC"
 #ifdef TIMING
 #ifdef MPI
@@ -326,6 +331,7 @@ C
       energia(23)=evdw_m
 c      print *," Processor",myrank," calls SUM_ENERGY"
       call sum_energy(energia,.true.)
       energia(23)=evdw_m
 c      print *," Processor",myrank," calls SUM_ENERGY"
       call sum_energy(energia,.true.)
+      if (dyn_ss) call dyn_set_nss
 c      print *," Processor",myrank," left SUM_ENERGY"
 #ifdef TIMING
 #ifdef MPI
 c      print *," Processor",myrank," left SUM_ENERGY"
 #ifdef TIMING
 #ifdef MPI
@@ -423,14 +429,14 @@ cMS$ATTRIBUTES C ::  proc_proc
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*etors+wscloc*escloc
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*etors+wscloc*escloc
-     & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
+     & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
      & +wbond*estr+Uconst+wsccor*esccor
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
      & +wang*ebe+wtor*etors+wscloc*escloc
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
      & +wbond*estr+Uconst+wsccor*esccor
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
      & +wang*ebe+wtor*etors+wscloc*escloc
-     & +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5
+     & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
      & +wbond*estr+Uconst+wsccor*esccor
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
      & +wbond*estr+Uconst+wsccor*esccor
@@ -471,7 +477,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'
@@ -774,6 +780,19 @@ c      enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         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
         time00=MPI_Wtime()
         call MPI_Barrier(FG_COMM,IERR)
         time_barrier_g=time_barrier_g+MPI_Wtime()-time00
@@ -784,8 +803,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)
      &    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
         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)
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
         write (iout,*) i,gloc(i,icg)
@@ -1029,25 +1058,25 @@ C------------------------------------------------------------------------
      &  edihcnstr,ebr*nss,
      &  Uconst,etot
    10 format (/'Virtual-chain energies:'//
      &  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.)'/
      & ' (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)'/ 
      & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
      & 'ESS=   ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
      & 'UCONST= ',1pE16.6,' (Constraint energy)'/ 
@@ -1420,7 +1449,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)
@@ -1529,6 +1558,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
+      include 'COMMON.SBRIDGE'
       logical lprn
       evdw=0.0D0
 ccccc      energy_dec=.false.
       logical lprn
       evdw=0.0D0
 ccccc      energy_dec=.false.
@@ -1557,6 +1587,27 @@ 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)
+            IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
+              call dyn_ssbond_ene(i,j,evdwij)
+              evdw=evdw+evdwij
+              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') 
+     &                        'evdw',i,j,evdwij,' ss'
+C triple bond artifac removal
+             do k=j+1,iend(i,iint) 
+C search over all next residues
+              if (dyn_ss_mask(k)) then
+C check if they are cysteins
+C              write(iout,*) 'k=',k
+              call triple_ssbond_ene(i,j,k,evdwij)
+C call the energy function that removes the artifical triple disulfide
+C bond the soubroutine is located in ssMD.F
+              evdw=evdw+evdwij             
+              if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
+     &                        'evdw',i,j,evdwij,'tss'
+              endif!dyn_ss_mask(k)
+             enddo! k
+            ELSE
+C            cycle
             ind=ind+1
             itypj=itype(j)
 c            dscj_inv=dsc_inv(itypj)
             ind=ind+1
             itypj=itype(j)
 c            dscj_inv=dsc_inv(itypj)
@@ -1666,6 +1717,7 @@ C Calculate angular part of the gradient.
 #else
             call sc_grad
 #endif
 #else
             call sc_grad
 #endif
+            ENDIF    ! dyn_ss            
           enddo      ! j
         enddo        ! iint
       enddo          ! i
           enddo      ! j
         enddo        ! iint
       enddo          ! i
@@ -4239,49 +4291,96 @@ 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
+cmc        if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
+        if (.not.dyn_ss .and. i.le.nss) then
+C 15/02/13 CC dynamic SSbond - additional check
+         if (ii.gt.nres 
+     &       .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then 
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
+         endif
 cd          write (iout,*) "eij",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
@@ -4343,7 +4442,7 @@ c      dscj_inv=dsc_inv(itypj)
       deltat12=om2-om1+2.0d0
       cosphi=om12-om1*om2
       eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
       deltat12=om2-om1+2.0d0
       cosphi=om12-om1*om2
       eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
-     &  +akct*deltad*deltat12
+     &  +akct*deltad*deltat12+ebr
      &  +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
 c      write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
 c     &  " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
      &  +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
 c      write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
 c     &  " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
@@ -4485,6 +4584,18 @@ c     write (*,'(a,i2)') 'EBEND ICG=',icg
 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)
 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)
+        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 +4629,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)
@@ -4724,6 +4847,9 @@ C
           enddo 
         endif
         if (i.lt.nres) then
           enddo 
         endif
         if (i.lt.nres) then
+
+        if (iabs(itype(i+1)).eq.20) iblock=2
+        if (iabs(itype(i+1)).ne.20) iblock=1
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -4744,7 +4870,7 @@ C
             sinph2(k)=0.0d0
           enddo
         endif  
             sinph2(k)=0.0d0
           enddo
         endif  
-        ethetai=aa0thet(ityp1,ityp2,ityp3)
+         ethetai=aa0thet(ityp1,ityp2,ityp3,iblock)
         do k=1,ndouble
           do l=1,k-1
             ccl=cosph1(l)*cosph2(k-l)
         do k=1,ndouble
           do l=1,k-1
             ccl=cosph1(l)*cosph2(k-l)
@@ -4766,11 +4892,12 @@ C
         enddo
         endif
         do k=1,ntheterm
         enddo
         endif
         do k=1,ntheterm
-          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3)*sinkt(k)
-          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3)
+          ethetai=ethetai+aathet(k,ityp1,ityp2,ityp3,iblock)*sinkt(k)
+          dethetai=dethetai+0.5d0*k*aathet(k,ityp1,ityp2,ityp3,iblock)
      &      *coskt(k)
           if (lprn)
      &      *coskt(k)
           if (lprn)
-     &    write (iout,*) "k",k," aathet",aathet(k,ityp1,ityp2,ityp3),
+     &    write (iout,*) "k",k,
+     &    "aathet",aathet(k,ityp1,ityp2,ityp3,iblock),
      &     " ethetai",ethetai
         enddo
         if (lprn) then
      &     " ethetai",ethetai
         enddo
         if (lprn) then
@@ -4789,24 +4916,24 @@ C
         endif
         do m=1,ntheterm2
           do k=1,nsingle
         endif
         do m=1,ntheterm2
           do k=1,nsingle
-            aux=bbthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)
-     &         +ccthet(k,m,ityp1,ityp2,ityp3)*sinph1(k)
-     &         +ddthet(k,m,ityp1,ityp2,ityp3)*cosph2(k)
-     &         +eethet(k,m,ityp1,ityp2,ityp3)*sinph2(k)
+            aux=bbthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)
+     &         +ccthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k)
+     &         +ddthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)
+     &         +eethet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k)
             ethetai=ethetai+sinkt(m)*aux
             dethetai=dethetai+0.5d0*m*aux*coskt(m)
             dephii=dephii+k*sinkt(m)*(
             ethetai=ethetai+sinkt(m)*aux
             dethetai=dethetai+0.5d0*m*aux*coskt(m)
             dephii=dephii+k*sinkt(m)*(
-     &          ccthet(k,m,ityp1,ityp2,ityp3)*cosph1(k)-
-     &          bbthet(k,m,ityp1,ityp2,ityp3)*sinph1(k))
+     &          ccthet(k,m,ityp1,ityp2,ityp3,iblock)*cosph1(k)-
+     &          bbthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph1(k))
             dephii1=dephii1+k*sinkt(m)*(
             dephii1=dephii1+k*sinkt(m)*(
-     &          eethet(k,m,ityp1,ityp2,ityp3)*cosph2(k)-
-     &          ddthet(k,m,ityp1,ityp2,ityp3)*sinph2(k))
+     &          eethet(k,m,ityp1,ityp2,ityp3,iblock)*cosph2(k)-
+     &          ddthet(k,m,ityp1,ityp2,ityp3,iblock)*sinph2(k))
             if (lprn)
      &      write (iout,*) "m",m," k",k," bbthet",
             if (lprn)
      &      write (iout,*) "m",m," k",k," bbthet",
-     &         bbthet(k,m,ityp1,ityp2,ityp3)," ccthet",
-     &         ccthet(k,m,ityp1,ityp2,ityp3)," ddthet",
-     &         ddthet(k,m,ityp1,ityp2,ityp3)," eethet",
-     &         eethet(k,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+     &         bbthet(k,m,ityp1,ityp2,ityp3,iblock)," ccthet",
+     &         ccthet(k,m,ityp1,ityp2,ityp3,iblock)," ddthet",
+     &         ddthet(k,m,ityp1,ityp2,ityp3,iblock)," eethet",
+     &         eethet(k,m,ityp1,ityp2,ityp3,iblock)," ethetai",ethetai
           enddo
         enddo
         if (lprn)
           enddo
         enddo
         if (lprn)
@@ -4814,28 +4941,33 @@ C
         do m=1,ntheterm3
           do k=2,ndouble
             do l=1,k-1
         do m=1,ntheterm3
           do k=2,ndouble
             do l=1,k-1
-              aux=ffthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)
+       aux=ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+     & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l)+
+     & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+     & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)
+
               ethetai=ethetai+sinkt(m)*aux
               dethetai=dethetai+0.5d0*m*coskt(m)*aux
               dephii=dephii+l*sinkt(m)*(
               ethetai=ethetai+sinkt(m)*aux
               dethetai=dethetai+0.5d0*m*coskt(m)*aux
               dephii=dephii+l*sinkt(m)*(
-     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)-
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)+
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+     & -ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)-
+     &  ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+     &  ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)+
+     &  ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
+
               dephii1=dephii1+(k-l)*sinkt(m)*(
               dephii1=dephii1+(k-l)*sinkt(m)*(
-     &           -ffthet(l,k,m,ityp1,ityp2,ityp3)*sinph1ph2(l,k)+
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)*sinph1ph2(k,l)+
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3)*cosph1ph2(l,k)-
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)*cosph1ph2(k,l))
+     &-ffthet(l,k,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(l,k)+
+     & ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)*sinph1ph2(k,l)+
+     & ggthet(l,k,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(l,k)-
+     & ggthet(k,l,m,ityp1,ityp2,ityp3,iblock)*cosph1ph2(k,l))
+
               if (lprn) then
               write (iout,*) "m",m," k",k," l",l," ffthet",
               if (lprn) then
               write (iout,*) "m",m," k",k," l",l," ffthet",
-     &            ffthet(l,k,m,ityp1,ityp2,ityp3),
-     &            ffthet(k,l,m,ityp1,ityp2,ityp3)," ggthet",
-     &            ggthet(l,k,m,ityp1,ityp2,ityp3),
-     &            ggthet(k,l,m,ityp1,ityp2,ityp3)," ethetai",ethetai
+     &            ffthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+     &            ffthet(k,l,m,ityp1,ityp2,ityp3,iblock)," ggthet",
+     &            ggthet(l,k,m,ityp1,ityp2,ityp3,iblock),
+     &            ggthet(k,l,m,ityp1,ityp2,ityp3,iblock),
+     &            " ethetai",ethetai
+
               write (iout,*) cosph1ph2(l,k)*sinkt(m),
      &            cosph1ph2(k,l)*sinkt(m),
      &            sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
               write (iout,*) cosph1ph2(l,k)*sinkt(m),
      &            cosph1ph2(k,l)*sinkt(m),
      &            sinph1ph2(l,k)*sinkt(m),sinph1ph2(k,l)*sinkt(m)
@@ -4844,9 +4976,11 @@ C
           enddo
         enddo
 10      continue
           enddo
         enddo
 10      continue
-        if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
-     &   i,theta(i)*rad2deg,phii*rad2deg,
+c        lprn1=.true.
+        if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') 
+     &  'ebe', i,theta(i)*rad2deg,phii*rad2deg,
      &   phii1*rad2deg,ethetai
      &   phii1*rad2deg,ethetai
+c        lprn1=.false.
         etheta=etheta+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
         etheta=etheta+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
@@ -5187,7 +5321,7 @@ C
         cosfac=dsqrt(cosfac2)
         sinfac2=0.5d0/(1.0d0-costtab(i+1))
         sinfac=dsqrt(sinfac2)
         cosfac=dsqrt(cosfac2)
         sinfac2=0.5d0/(1.0d0-costtab(i+1))
         sinfac=dsqrt(sinfac2)
-        it=itype(i)
+        it=iabs(itype(i))
         if (it.eq.10) goto 1
 c
 C  Compute the axes of tghe local cartesian coordinates system; store in
         if (it.eq.10) goto 1
 c
 C  Compute the axes of tghe local cartesian coordinates system; store in
@@ -5205,7 +5339,7 @@ C     &   dc_norm(3,i+nres)
           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
         enddo
         do j = 1,3
           y_prime(j) = (dc_norm(j,i) + dc_norm(j,i-1))*sinfac
         enddo
         do j = 1,3
-          z_prime(j) = -uz(j,i-1)
+          z_prime(j) = -uz(j,i-1)*dsign(1.0d0,dfloat(itype(i)))
         enddo     
 c       write (2,*) "i",i
 c       write (2,*) "x_prime",(x_prime(j),j=1,3)
         enddo     
 c       write (2,*) "i",i
 c       write (2,*) "x_prime",(x_prime(j),j=1,3)
@@ -5237,7 +5371,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 +5379,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, dfloat(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
@@ -5415,8 +5549,10 @@ c     &   (dC_norm(j,i-1),j=1,3)," vbld_inv",vbld_inv(i+1),vbld_inv(i)
          dZZ_Ci1(k)=0.0d0
          dZZ_Ci(k)=0.0d0
          do j=1,3
          dZZ_Ci1(k)=0.0d0
          dZZ_Ci(k)=0.0d0
          do j=1,3
-           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)*dC_norm(j,i+nres)
-           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)*dC_norm(j,i+nres)
+           dZZ_Ci(k)=dZZ_Ci(k)-uzgrad(j,k,2,i-1)
+     &     *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
+           dZZ_Ci1(k)=dZZ_Ci1(k)-uzgrad(j,k,1,i-1)
+     &     *dsign(1.0d0,dfloat(itype(i)))*dC_norm(j,i+nres)
          enddo
           
          dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
          enddo
           
          dXX_XYZ(k)=vbld_inv(i+nres)*(x_prime(k)-xx*dC_norm(k,i+nres))
@@ -5701,15 +5837,22 @@ C Set lprn=.true. for debugging
 c     lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
 c     lprn=.true.
       etors=0.0D0
       do i=iphi_start,iphi_end
-      etors_ii=0.0D0
+        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1 
+     &       .or. itype(i).eq.ntyp1) cycle
+        etors_ii=0.0D0
+         if (iabs(itype(i)).eq.20) then
+         iblock=2
+         else
+         iblock=1
+         endif
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         phii=phi(i)
         gloci=0.0D0
 C Regular cosine and sine terms
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         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 +5867,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 +5880,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
@@ -5794,7 +5938,10 @@ C Set lprn=.true. for debugging
       lprn=.false.
 c     lprn=.true.
       etors_d=0.0D0
       lprn=.false.
 c     lprn=.true.
       etors_d=0.0D0
+c      write(iout,*) "a tu??"
       do i=iphid_start,iphid_end
       do i=iphid_start,iphid_end
+        if (itype(i-2).eq.ntyp1 .or. itype(i-1).eq.ntyp1
+     &      .or. itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
         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))
@@ -5802,11 +5949,15 @@ c     lprn=.true.
         phii1=phi(i+1)
         gloci1=0.0D0
         gloci2=0.0D0
         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)
+        iblock=1
+        if (iabs(itype(i+1)).eq.20) iblock=2
+
+C Regular cosine and sine terms
+        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)
@@ -5816,12 +5967,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)
@@ -5836,7 +5987,6 @@ 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,10 +6017,11 @@ c        amino-acid residues.
 C Set lprn=.true. for debugging
       lprn=.false.
 c      lprn=.true.
 C Set lprn=.true. for debugging
       lprn=.false.
 c      lprn=.true.
-c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
+c     write (iout,*) "EBACK_SC_COR",itau_start,itau_end
       esccor=0.0D0
       esccor=0.0D0
-      do i=iphi_start-1,iphi_end+1
+      do i=itau_start,itau_end
         esccor_ii=0.0D0
         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))
         phii=phi(i)
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
         phii=phi(i)
@@ -5880,22 +6031,26 @@ 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,1 !intertyp
+        
+        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
+c   3 = SC...Ca...Ca...SCi
+        gloci=0.0D0
         if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
         if (((intertyp.eq.3).and.((itype(i-2).eq.10).or.
-     &      (itype(i-1).eq.10).or.(itype(i-2).eq.21).or.
-     &      (itype(i-1).eq.21)))
+     &      (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. ((intertyp.eq.1).and.((itype(i-2).eq.10)
-     &     .or.(itype(i-2).eq.21)))
+     &     .or.(itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)
+     &     .or.(itype(i).eq.ntyp1)))
      &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
      &    .or.((intertyp.eq.2).and.((itype(i-1).eq.10).or.
-     &      (itype(i-1).eq.21)))) cycle  
-        if ((intertyp.eq.2).and.(i.le.iphi_start-1)) cycle
-        if ((intertyp.eq.1).and.(i.ge.iphi_end+1)) cycle
+     &      (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))
+     & 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)
@@ -5905,13 +6060,13 @@ c   3 = SC...Ca...Ca...SC
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
-c        write (iout,*) "WTF",intertyp,i,itype(i),
+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/)')
 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)
+     &  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)
         gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
        enddo !intertyp
       enddo
         gsccor_loc(i-3)=gsccor_loc(i-3)+gloci
        enddo !intertyp
       enddo
@@ -8024,9 +8179,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                
@@ -8121,22 +8276,22 @@ c----------------------------------------------------------------------------
       include 'COMMON.GEO'
       logical swap
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2),
       include 'COMMON.GEO'
       logical swap
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2),
-     & auxvec1(2),auxvec2(1),auxmat1(2,2)
+     & auxvec1(2),auxvec2(2),auxmat1(2,2)
       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, 
@@ -8307,18 +8462,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 
@@ -8424,18 +8579,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