imax0 replaced with max0
[unres.git] / source / unres / src_MD / energy_p_new_barrier.F
index b45c847..4bccdf3 100644 (file)
@@ -1679,9 +1679,10 @@ c     &        " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
      &        evdwij
             endif
 
-            if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') 
-     &                        'evdw',i,j,evdwij
-
+            if (energy_dec) then
+              write (iout,'(a6,2i5,0pf7.3)') 'evdw',i,j,evdwij
+              call flush(iout)
+            endif
 C Calculate gradient components.
             e1=e1*eps1*eps2rt**2*eps3rt**2
             fac=-expon*(e1+evdwij)*rij_shift
@@ -3023,6 +3024,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 +3042,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 +3064,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 +3078,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
@@ -4482,6 +4494,8 @@ c
       do i=ibondp_start,ibondp_end
         diff = vbld(i)-vbldp0
 c        write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
+        if (energy_dec)    write (iout,'(a7,i5,4f7.3)') 
+     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
         estr=estr+diff*diff
         do j=1,3
           gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
@@ -4500,6 +4514,12 @@ c
             diff=vbld(i+nres)-vbldsc0(1,iti)
 c            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
 c     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
+            if (energy_dec)  then
+              write (iout,*) 
+     &         "estr sc",i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+     &         AKSC(1,iti),AKSC(1,iti)*diff*diff
+              call flush(iout)
+            endif
             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
             do j=1,3
               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)
@@ -4778,6 +4798,8 @@ C
       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
@@ -4787,7 +4809,8 @@ C
           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(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
           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
@@ -4822,7 +4845,7 @@ C
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
           enddo
         enddo
 10      continue
-        lprn1=.true.
         if (lprn1) write (iout,'(a4,i2,3f8.1,9h ethetai ,f10.5)') 
      &  'ebe', i,theta(i)*rad2deg,phii*rad2deg,
      &   phii1*rad2deg,ethetai
-        lprn1=.false.
         etheta=etheta+ethetai
+        if (energy_dec) write (iout,'(a6,i5,0pf7.3)')
+     &      'ebend',i,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
         gloc(nphi+i-2,icg)=wang*dethetai
@@ -5875,12 +5898,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))
@@ -5899,6 +5924,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
@@ -5914,12 +5941,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)
@@ -5956,18 +5988,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
-        esccor_ii=0.0D0
+C        do i=42,42
+
+        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
+         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
@@ -5989,9 +6025,13 @@ 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        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)
@@ -6004,7 +6044,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