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 64dd328..47583a4 100644 (file)
@@ -134,7 +134,7 @@ C
 cmc
 cmc Sep-06: egb takes care of dynamic ss bonds too
 cmc
-      if (dyn_ss) call dyn_set_nss
+c      if (dyn_ss) call dyn_set_nss
 
 c      print *,"Processor",myrank," computed USCSC"
 #ifdef TIMING
@@ -331,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
@@ -428,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
@@ -1057,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)'/ 
@@ -1589,6 +1590,8 @@ C
             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=itype(j)
@@ -3020,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)
@@ -3035,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)
@@ -3053,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)
@@ -3065,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
@@ -4424,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,
@@ -4775,6 +4789,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
@@ -4784,7 +4800,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(i-3).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
@@ -4819,7 +4836,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
@@ -4925,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
@@ -5951,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
@@ -5987,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)
@@ -5999,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