energy of unres src_MD-M match prerelease

[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 5eca704..fc95e25 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -3324,21 +3324,21 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups
 C
 C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
       do i=iturn3_start,iturn3_end
-        if (i.le.1) cycle
+CAna        if (i.le.1) cycle
 C        write(iout,*) "tu jest i",i
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((i+4).gt.nres)
-     & .or.((i-1).le.0)
+CAna     & .or.((i+4).gt.nres)
+CAna     & .or.((i-1).le.0)
 C end of changes by Ana
      &  .or. itype(i+2).eq.ntyp1
      &  .or. itype(i+3).eq.ntyp1) cycle
-        if(i.gt.1)then
-          if(itype(i-1).eq.ntyp1)cycle
-        end if
-        if(i.LT.nres-3)then
-          if (itype(i+4).eq.ntyp1) cycle
-        end if
+CAna        if(i.gt.1)then
+CAna          if(itype(i-1).eq.ntyp1)cycle
+CAna        end if
+CAna        if(i.LT.nres-3)then
+CAna          if (itype(i+4).eq.ntyp1) cycle
+CAna        end if
         dxi=dc(1,i)
         dyi=dc(2,i)
         dzi=dc(3,i)
@@ -3360,17 +3360,17 @@ C end of changes by Ana
         num_cont_hb(i)=num_conti
       enddo
       do i=iturn4_start,iturn4_end
-        if (i.le.1) cycle
+cAna        if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((i+5).gt.nres)
-     & .or.((i-1).le.0)
+cAna     & .or.((i+5).gt.nres)
+cAna     & .or.((i-1).le.0)
 C end of changes suggested by Ana
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1
-     &    .or. itype(i+5).eq.ntyp1
-     &    .or. itype(i).eq.ntyp1
-     &    .or. itype(i-1).eq.ntyp1
+cAna     &    .or. itype(i+5).eq.ntyp1
+cAna     &    .or. itype(i).eq.ntyp1
+cAna     &    .or. itype(i-1).eq.ntyp1
      &                             ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -3428,14 +3428,14 @@ 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
-        if (i.le.1) cycle
+cAna        if (i.le.1) cycle
         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((i+2).gt.nres)
-     & .or.((i-1).le.0)
+cAna     & .or.((i+2).gt.nres)
+cAna     & .or.((i-1).le.0)
 C end of changes by Ana
-     &  .or. itype(i+2).eq.ntyp1
-     &  .or. itype(i-1).eq.ntyp1
+cAna     &  .or. itype(i+2).eq.ntyp1
+cAna     &  .or. itype(i-1).eq.ntyp1
      &                ) cycle
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -3486,14 +3486,14 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         num_conti=num_cont_hb(i)
         do j=ielstart(i),ielend(i)
 C          write (iout,*) i,j
-         if (j.le.1) cycle
+cAna         if (j.le.1) cycle
           if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((j+2).gt.nres)
-     & .or.((j-1).le.0)
+cAna     & .or.((j+2).gt.nres)
+cAna     & .or.((j-1).le.0)
 C end of changes by Ana
-     & .or.itype(j+2).eq.ntyp1
-     & .or.itype(j-1).eq.ntyp1
+cAna     & .or.itype(j+2).eq.ntyp1
+cAna     & .or.itype(j-1).eq.ntyp1
      &) cycle
           call eelecij(i,j,ees,evdw1,eel_loc)
         enddo ! j
@@ -5684,7 +5684,7 @@ C
       logical lprn /.false./, lprn1 /.false./
       etheta=0.0D0
       do i=ithet_start,ithet_end
-        if (i.eq.2) cycle
+c        if (i.eq.2) cycle
 c        print *,i,itype(i-1),itype(i),itype(i-2)
         if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1)
      &  .or.(itype(i).eq.ntyp1)) cycle
@@ -5701,7 +5701,7 @@ C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -5716,7 +5716,7 @@ C propagation of chirality for glycine type
           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
@@ -5737,7 +5737,7 @@ C propagation of chirality for glycine type
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -7550,12 +7550,12 @@ c      write (iout,*) "EBACK_SC_COR",itau_start,itau_end
       esccor=0.0D0
       do i=itau_start,itau_end
         if ((itype(i-2).eq.ntyp1).or.(itype(i-1).eq.ntyp1)) cycle
-        esccor_ii=0.0D0
         isccori=isccortyp(itype(i-2))
         isccori1=isccortyp(itype(i-1))
 c      write (iout,*) "EBACK_SC_COR",i,nterm_sccor(isccori,isccori1)
         phii=phi(i)
         do intertyp=1,3 !intertyp
+         esccor_ii=0.0D0
 cc Added 09 May 2012 (Adasko)
 cc  Intertyp means interaction type of backbone mainchain correlation: 
 c   1 = SC...Ca...Ca...Ca
@@ -7579,9 +7579,12 @@ c   3 = SC...Ca...Ca...SCi
           v2ij=v2sccor(j,intertyp,isccori,isccori1)
           cosphi=dcos(j*tauangle(intertyp,i))
           sinphi=dsin(j*tauangle(intertyp,i))
+          if (energy_dec) esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi
           esccor=esccor+v1ij*cosphi+v2ij*sinphi
           gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
         enddo
+         if (energy_dec) write (iout,'(a6,i5,i2,0pf7.3)')
+     &         'esccor',i,intertyp,esccor_ii
 c      write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp
         gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
         if (lprn)