Merge branch 'prerelease-3.2.1' of mmka.chem.univ.gda.pl:unres into prerelease-3.2.1
[unres.git] / source / unres / src_MD / energy_p_new_barrier.F
index e943ce5..47583a4 100644 (file)
@@ -131,6 +131,11 @@ C
 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
@@ -326,6 +331,7 @@ C
       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
@@ -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
-     & +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
-     & +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
@@ -774,7 +780,6 @@ c      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
@@ -783,7 +788,6 @@ c      enddo
        enddo
       enddo
 #endif
-#undef DEBUG
         do i=1,nres
          do j=1,3
           gloc_scbuf(j,i)=gloc_sc(j,i,icg)
@@ -802,7 +806,6 @@ c      enddo
         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
-#define DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc_sc after reduce"
       do i=1,nres
@@ -811,7 +814,6 @@ c      enddo
        enddo
       enddo
 #endif
-#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
@@ -1056,25 +1058,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)'/
-     & 'EHPB=  ',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)'/ 
@@ -1137,8 +1139,8 @@ C
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
-        itypi1=iabs(itype(i+1))
+        itypi=itype(i)
+        itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1151,7 +1153,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)
-            itypj=iabs(itype(j))
+            itypj=itype(j)
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -1314,8 +1316,8 @@ C
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
-        itypi1=iabs(itype(i+1))
+        itypi=itype(i)
+        itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1324,7 +1326,7 @@ C Calculate SC interaction energy.
 C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
-            itypj=iabs(itype(j))
+            itypj=itype(j)
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -1431,8 +1433,8 @@ c     else
 c     endif
       ind=0
       do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
-        itypi1=iabs(itype(i+1))
+        itypi=itype(i)
+        itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1556,6 +1558,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
+      include 'COMMON.SBRIDGE'
       logical lprn
       evdw=0.0D0
 ccccc      energy_dec=.false.
@@ -1567,8 +1570,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
-        itypi=iabs(itype(i))
-        itypi1=iabs(itype(i+1))
+        itypi=itype(i)
+        itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1584,8 +1587,14 @@ C Calculate SC interaction energy.
 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'
+            ELSE
             ind=ind+1
-            itypj=iabs(itype(j))
+            itypj=itype(j)
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
 c            write (iout,*) "j",j,dsc_inv(itypj),dscj_inv,
@@ -1693,6 +1702,7 @@ C Calculate angular part of the gradient.
 #else
             call sc_grad
 #endif
+            ENDIF    ! dyn_ss            
           enddo      ! j
         enddo        ! iint
       enddo          ! i
@@ -1726,8 +1736,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
-        itypi=iabs(itype(i))
-        itypi1=iabs(itype(i+1))
+        itypi=itype(i)
+        itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -1742,7 +1752,7 @@ C
         do iint=1,nint_gr(i)
           do j=istart(i,iint),iend(i,iint)
             ind=ind+1
-            itypj=iabs(itype(j))
+            itypj=itype(j)
 c            dscj_inv=dsc_inv(itypj)
             dscj_inv=vbld_inv(j+nres)
             sig0ij=sigma(itypi,itypj)
@@ -2046,8 +2056,8 @@ C
 cd    print *,'Entering Esoft_sphere nnt=',nnt,' nct=',nct
       evdw=0.0D0
       do i=iatsc_s,iatsc_e
-        itypi=iabs(itype(i))
-        itypi1=iabs(itype(i+1))
+        itypi=itype(i)
+        itypi1=itype(i+1)
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
@@ -2058,7 +2068,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)
-            itypj=iabs(itype(j))
+            itypj=itype(j)
             xj=c(1,nres+j)-xi
             yj=c(2,nres+j)-yi
             zj=c(3,nres+j)-zi
@@ -3013,6 +3023,9 @@ C
 C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
       do i=iturn3_start,iturn3_end
+C        if (itype(i).eq.21 .or. itype(i+1).eq.21
+C     &  .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21)
+C     &  cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -3028,6 +3041,10 @@ C
         num_cont_hb(i)=num_conti
       enddo
       do i=iturn4_start,iturn4_end
+C        if (itype(i).eq.21 .or. itype(i+1).eq.21
+C     &  .or. itype(i+2).eq.21 .or. itype(i+3).eq.21.or.itype(i+4).eq.21
+C     &  .or. itype(i+5).eq.21)
+C     & cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -3046,6 +3063,8 @@ c
 c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
 c
       do i=iatel_s,iatel_e
+C          if (itype(i).eq.21 .or. itype(i+1).eq.21
+C     &.or.itype(i+2)) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -3058,6 +3077,8 @@ c
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
+C          if (itype(j).eq.21 .or. itype(j+1).eq.21
+C     &.or.itype(j+2)) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
         enddo ! j
         num_cont_hb(i)=num_conti
@@ -4059,7 +4080,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)
-          itypj=iabs(itype(j))
+          itypj=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
@@ -4153,7 +4174,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)
-          itypj=iabs(itype(j))
+          itypj=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
@@ -4270,10 +4291,15 @@ 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.
-        if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. iabs(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
+         endif
 cd          write (iout,*) "eij",eij
         else if (ii.gt.nres .and. jj.gt.nres) then
 c Restraints from contact prediction
@@ -4375,7 +4401,7 @@ C
       include 'COMMON.VAR'
       include 'COMMON.IOUNITS'
       double precision erij(3),dcosom1(3),dcosom2(3),gg(3)
-      itypi=iabs(itype(i))
+      itypi=itype(i)
       xi=c(1,nres+i)
       yi=c(2,nres+i)
       zi=c(3,nres+i)
@@ -4384,7 +4410,7 @@ C
       dzi=dc_norm(3,nres+i)
 c      dsci_inv=dsc_inv(itypi)
       dsci_inv=vbld_inv(nres+i)
-      itypj=iabs(itype(j))
+      itypj=itype(j)
 c      dscj_inv=dsc_inv(itypj)
       dscj_inv=vbld_inv(nres+j)
       xj=c(1,nres+j)-xi
@@ -4412,7 +4438,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)
-     &  +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,
@@ -4478,7 +4504,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
-        iti=iabs(itype(i))
+        iti=itype(i)
         if (iti.ne.10) then
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
@@ -4553,19 +4579,7 @@ 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)
-        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
+        it=itype(i-1)
         if (i.gt.3) then
 #ifdef OSF
          phii=phi(i)
@@ -4599,27 +4613,15 @@ 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
-          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
+          athetk=athet(k,it)
+          bthetk=bthet(k,it)
           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,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
+        dthetg1=(-athet(1,it)*y(2)+athet(2,it)*y(1))*ss
+        dthetg2=(-bthet(1,it)*z(2)+bthet(2,it)*z(1))*ss
         if (theta(i).gt.pi-delta) then
           call theteng(pi-delta,thet_pred_mean,theta0(it),f0,fprim0,
      &         E_tc0)
       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
         theti2=0.5d0*theta(i)
-        ityp2=ithetyp(iabs(itype(i-1)))
+        ityp2=ithetyp(itype(i-1))
         do k=1,nntheterm
           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(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
 #else
           phii=phi(i)
 #endif
-          ityp1=ithetyp(iabs(itype(i-2)))
+          ityp1=ithetyp(itype(i-2))
           do k=1,nsingle
             cosph1(k)=dcos(k*phii)
             sinph1(k)=dsin(k*phii)
           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
 #else
           phii1=phi(i+1)
 #endif
-          ityp3=ithetyp(iabs(itype(i)))
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=dcos(k*phii1)
             sinph2(k)=dsin(k*phii1)
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -4937,8 +4942,8 @@ C
           enddo
         enddo
 10      continue
-        if (lprn1) write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') 
-     &   i,theta(i)*rad2deg,phii*rad2deg,
+        if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') 
+     &  'ebe', i,theta(i)*rad2deg,phii*rad2deg,
      &   phii1*rad2deg,ethetai
         etheta=etheta+ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
@@ -4975,7 +4980,7 @@ c     write (iout,'(a)') 'ESC'
       do i=loc_start,loc_end
         it=itype(i)
         if (it.eq.10) goto 1
-        nlobit=nlob(iabs(it))
+        nlobit=nlob(it)
 c       print *,'i=',i,' it=',it,' nlobit=',nlobit
 c       write (iout,*) 'i=',i,' ssa=',ssa,' ssad=',ssad
         theti=theta(i+1)-pipol
@@ -5132,11 +5137,11 @@ C Compute the contribution to SC energy and derivatives
 
           do j=1,nlobit
 #ifdef OSF
-            adexp=bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin
+            adexp=bsc(j,it)-0.5D0*contr(j,iii)+emin
             if(adexp.ne.adexp) adexp=1.0
             expfac=dexp(adexp)
 #else
-            expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j,iii)+emin)
+            expfac=dexp(bsc(j,it)-0.5D0*contr(j,iii)+emin)
 #endif
 cd          print *,'j=',j,' expfac=',expfac
             escloc_i=escloc_i+expfac
@@ -5218,7 +5223,7 @@ C Compute the contribution to SC energy and derivatives
 
       dersc12=0.0d0
       do j=1,nlobit
-        expfac=dexp(bsc(j,iabs(it))-0.5D0*contr(j)+emin)
+        expfac=dexp(bsc(j,it)-0.5D0*contr(j)+emin)
         escloc_i=escloc_i+expfac
         do k=1,2
           dersc(k)=dersc(k)+Ax(k,j)*expfac
@@ -5797,17 +5802,12 @@ c     lprn=.true.
       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
-        do j=1,nterm(itori,itori1,iblock)
-          v1ij=v1(j,itori,itori1,iblock)
-          v2ij=v2(j,itori,itori1,iblock)
+        do j=1,nterm(itori,itori1)
+          v1ij=v1(j,itori,itori1)
+          v2ij=v2(j,itori,itori1)
           cosphi=dcos(j*phii)
           sinphi=dsin(j*phii)
           etors=etors+v1ij*cosphi+v2ij*sinphi
@@ -5822,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)
-        do j=1,nlor(itori,itori1,iblock)
+        do j=1,nlor(itori,itori1)
           vl1ij=vlor1(j,itori,itori1)
           vl2ij=vlor2(j,itori,itori1)
           vl3ij=vlor3(j,itori,itori1)
           gloci=gloci+vl1ij*(vl3ij*cosphi-vl2ij*sinphi)*pom
         enddo
 C Subtract the constant term
-        etors=etors-v0(itori,itori1,iblock)
+        etors=etors-v0(itori,itori1)
           if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
-     &         'etor',i,etors_ii-v0(itori,itori1,iblock)
+     &         'etor',i,etors_ii-v0(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,iblock),j=1,6),
-     &  (v2(j,itori,itori1,iblock),j=1,6)
+     &  (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)
       enddo
@@ -5897,17 +5896,15 @@ c     lprn=.true.
         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
-        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)
+        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)
           cosphi1=dcos(j*phii)
           sinphi1=dsin(j*phii)
           cosphi2=dcos(j*phii1)
@@ -5917,12 +5914,12 @@ c     lprn=.true.
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
-        do k=2,ntermd_2(itori,itori1,itori2,iblock)
+        do k=2,ntermd_2(itori,itori1,itori2)
           do l=1,k-1
-            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)
+            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)
             cosphi1p2=dcos(l*phii+(k-l)*phii1)
             cosphi1m2=dcos(l*phii-(k-l)*phii1)
             sinphi1p2=dsin(l*phii+(k-l)*phii1)
@@ -5971,18 +5968,22 @@ c      lprn=.true.
 c      write (iout,*) "EBACK_SC_COR",iphi_start,iphi_end,nterm_sccor
       esccor=0.0D0
       do i=itau_start,itau_end
+C        do i=42,42
         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)
+
 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
+c        do intertyp=1,3 !intertyp
+        do intertyp=2,2 !intertyp
 cc Added 09 May 2012 (Adasko)
 cc  Intertyp means interaction type of backbone mainchain correlation: 
 c   1 = SC...Ca...Ca...Ca
@@ -6007,6 +6008,7 @@ c   3 = SC...Ca...Ca...SCi
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
+C        print *,i,tauangle(1,i),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)
@@ -6019,7 +6021,8 @@ c     &gloc_sc(intertyp,i-3,icg)
        enddo !intertyp
       enddo
 c        do i=1,nres
-c        write (iout,*) "W@T@F",  gloc_sc(1,i,icg),gloc(i,icg)
+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
@@ -8224,7 +8227,7 @@ c----------------------------------------------------------------------------
       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