Merge branch 'devel' of mmka:unres into devel

[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 b7e216c..5ad53f1 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -121,6 +121,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
       time01=MPI_Wtime() 
@@ -300,6 +305,7 @@ c    Here are the energies showed per procesor if the are more processors
 c    per molecule then we sum it up in sum_energy subroutine 
 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
       time_sumene=time_sumene+MPI_Wtime()-time00
@@ -389,14 +395,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
@@ -710,7 +716,7 @@ c      enddo
         do i=1,4*nres
           glocbuf(i)=gloc(i,icg)
         enddo
-#define DEBUG
+c#define DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc_sc before reduce"
       do i=1,nres
@@ -719,7 +725,7 @@ c      enddo
        enddo
       enddo
 #endif
-#undef DEBUG
+c#undef DEBUG
         do i=1,nres
          do j=1,3
           gloc_scbuf(j,i)=gloc_sc(j,i,icg)
@@ -739,7 +745,7 @@ 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
+c#define DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc_sc after reduce"
       do i=1,nres
@@ -748,7 +754,7 @@ c      enddo
        enddo
       enddo
 #endif
-#undef DEBUG
+c#undef DEBUG
 #ifdef DEBUG
       write (iout,*) "gloc after reduce"
       do i=1,4*nres
@@ -1412,6 +1418,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
+      include 'COMMON.SBRIDGE'
       logical lprn
       evdw=0.0D0
 ccccc      energy_dec=.false.
@@ -1439,6 +1446,12 @@ 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))
             if (itypj.eq.ntyp1) cycle
@@ -1533,6 +1546,7 @@ C Calculate the radial part of the gradient
             gg(3)=zj*fac
 C Calculate angular part of the gradient.
             call sc_grad
+            ENDIF    ! dyn_ss            
           enddo      ! j
         enddo        ! iint
       enddo          ! i
@@ -4059,14 +4073,20 @@ C iii and jjj point to the residues for which the distance is assigned.
           iii=ii
           jjj=jj
         endif
-cd        write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
+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.
+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. iabs(itype(iii)).eq.1 .and.
      & iabs(itype(jjj)).eq.1) then
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
 cd          write (iout,*) "eij",eij
+         endif
         else
 C Calculate the distance between the two points and its difference from the
 C target distance.
@@ -4165,7 +4185,7 @@ c      dscj_inv=dsc_inv(itypj)
       cosphi=om12-om1*om2
       eij=akcm*deltad*deltad+akth*(deltat1*deltat1+deltat2*deltat2)
      &  +akct*deltad*deltat12
-     &  +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi
+     &  +v1ss*cosphi+v2ss*cosphi*cosphi+v3ss*cosphi*cosphi*cosphi+ebr
 c      write(iout,*) i,j,"rij",rij,"d0cm",d0cm," akcm",akcm," akth",akth,
 c     &  " akct",akct," deltad",deltad," deltat",deltat1,deltat2,
 c     &  " deltat12",deltat12," eij",eij