energy_dec works with etors_d

[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 0bcd6f3..0c6dd2e 100644 (file)
--- a/source/unres/src_MD/energy_p_new_barrier.F
+++ b/source/unres/src_MD/energy_p_new_barrier.F
@@ -3023,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)
@@ -3038,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)
@@ -3056,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)
@@ -3068,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
@@ -4931,11 +4942,9 @@ C        if (i.gt.3) then
           enddo
         enddo
 10      continue
-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
-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
@@ -5878,12 +5887,14 @@ C 6/23/01 Compute double torsional energy
       include 'COMMON.IOUNITS'
       include 'COMMON.FFIELD'
       include 'COMMON.TORCNSTR'
+      include 'COMMON.CONTROL'
       logical lprn
 C Set lprn=.true. for debugging
       lprn=.false.
 c     lprn=.true.
       etors_d=0.0D0
       do i=iphid_start,iphid_end
+        etors_d_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         itori2=itortyp(itype(i))
@@ -5902,6 +5913,8 @@ c     lprn=.true.
           sinphi2=dsin(j*phii1)
           etors_d=etors_d+v1cij*cosphi1+v1sij*sinphi1+
      &     v2cij*cosphi2+v2sij*sinphi2
+          if (energy_dec) etors_d_ii=etors_d_ii+
+     &     v1cij*cosphi1+v1sij*sinphi1+v2cij*cosphi2+v2sij*sinphi2
           gloci1=gloci1+j*(v1sij*cosphi1-v1cij*sinphi1)
           gloci2=gloci2+j*(v2sij*cosphi2-v2cij*sinphi2)
         enddo
@@ -5917,12 +5930,17 @@ c     lprn=.true.
             sinphi1m2=dsin(l*phii-(k-l)*phii1)
             etors_d=etors_d+v1cdij*cosphi1p2+v2cdij*cosphi1m2+
      &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
+            if (energy_dec) etors_d_ii=etors_d_ii+
+     &        v1cdij*cosphi1p2+v2cdij*cosphi1m2+
+     &        v1sdij*sinphi1p2+v2sdij*sinphi1m2
             gloci1=gloci1+l*(v1sdij*cosphi1p2+v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2-v2cdij*sinphi1m2)
             gloci2=gloci2+(k-l)*(v1sdij*cosphi1p2-v2sdij*cosphi1m2
      &        -v1cdij*sinphi1p2+v2cdij*sinphi1m2) 
           enddo
         enddo
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &        'etor_d',i,etors_d_ii
         gloc(i-3,icg)=gloc(i-3,icg)+wtor_d*gloci1
         gloc(i-2,icg)=gloc(i-2,icg)+wtor_d*gloci2
 c        write (iout,*) "gloci", gloc(i-3,icg)
@@ -5959,19 +5977,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
@@ -5996,6 +6017,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)
@@ -6008,7 +6030,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