Merge branch 'prerelease-3.2.1' into czarek
[unres.git] / source / unres / src_MD-M / energy_p_new_barrier.F
index dee7c2f..91046c3 100644 (file)
@@ -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() 
@@ -298,6 +303,7 @@ C
       energia(21)=esccor
 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
@@ -1382,6 +1388,7 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.CALC'
       include 'COMMON.CONTROL'
+      include 'COMMON.SBRIDGE'
       logical lprn
       evdw=0.0D0
 ccccc      energy_dec=.false.
@@ -1409,6 +1416,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=itype(j)
             if (itypj.eq.21) cycle
@@ -1503,6 +1516,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
@@ -4021,49 +4035,56 @@ 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. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+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. itype(iii).eq.1 .and. itype(jjj).eq.1) then 
           call ssbond_ene(iii,jjj,eij)
           ehpb=ehpb+2*eij
+         endif
 cd          write (iout,*) "eij",eij
         else
 C Calculate the distance between the two points and its difference from the
 C target distance.
-        dd=dist(ii,jj)
-        rdis=dd-dhpb(i)
+          dd=dist(ii,jj)
+            rdis=dd-dhpb(i)
 C Get the force constant corresponding to this distance.
-        waga=forcon(i)
+            waga=forcon(i)
 C Calculate the contribution to energy.
-        ehpb=ehpb+waga*rdis*rdis
+            ehpb=ehpb+waga*rdis*rdis
 C
 C Evaluate gradient.
 C
-        fac=waga*rdis/dd
+            fac=waga*rdis/dd
 cd      print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
 cd   &   ' waga=',waga,' fac=',fac
-        do j=1,3
-          ggg(j)=fac*(c(j,jj)-c(j,ii))
-        enddo
+            do j=1,3
+              ggg(j)=fac*(c(j,jj)-c(j,ii))
+            enddo
 cd      print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
 C If this is a SC-SC distance, we need to calculate the contributions to the
 C Cartesian gradient in the SC vectors (ghpbx).
-        if (iii.lt.ii) then
+          if (iii.lt.ii) then
           do j=1,3
             ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
             ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
           enddo
-        endif
+          endif
 cgrad        do j=iii,jjj-1
 cgrad          do k=1,3
 cgrad            ghpbc(k,j)=ghpbc(k,j)+ggg(k)
 cgrad          enddo
 cgrad        enddo
-        do k=1,3
-          ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
-          ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
-        enddo
+          do k=1,3
+            ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+            ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+          enddo
         endif
       enddo
       ehpb=0.5D0*ehpb
@@ -5504,7 +5525,7 @@ c     lprn=.true.
       do i=iphi_start,iphi_end
         if (itype(i-2).eq.21 .or. itype(i-1).eq.21 
      &       .or. itype(i).eq.21) cycle
-        etors_ii=0.0D0
+      etors_ii=0.0D0
         itori=itortyp(itype(i-2))
         itori1=itortyp(itype(i-1))
         phii=phi(i)
@@ -7907,12 +7928,12 @@ C                                                                              C
 C          o             o                                                     C
 C     \   /l\           /j\   /                                                C
 C      \ /   \         /   \ /                                                 C
-C       o| o |         | o |o                                                  C
+C       o| o |         | o |o                                                  C                
 C     \ j|/k\|      \  |/k\|l                                                  C
 C      \ /   \       \ /   \                                                   C
 C       o             o                                                        C
-C       i             i                                                        C
-C                                                                              C
+C       i             i                                                        C 
+C                                                                              C           
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 cd      write (2,*) 'eello6_graph2: i,',i,' j',j,' k',k,' l',l
 C AL 7/4/01 s1 would occur in the sixth-order moment, 
@@ -8083,10 +8104,10 @@ c----------------------------------------------------------------------------
       double precision vv(2),pizda(2,2),auxmat(2,2),auxvec(2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                                                              C
+C                                                                              C 
 C      Parallel       Antiparallel                                             C
 C                                                                              C
-C          o             o                                                     C
+C          o             o                                                     C 
 C         /l\   /   \   /j\                                                    C 
 C        /   \ /     \ /   \                                                   C
 C       /| o |o       o| o |\                                                  C
@@ -8200,7 +8221,7 @@ c----------------------------------------------------------------------------
      & auxvec1(2),auxmat1(2,2)
       logical swap
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-C                                                                              C
+C                                                                              C                       
 C      Parallel       Antiparallel                                             C
 C                                                                              C
 C          o             o                                                     C
@@ -8208,10 +8229,10 @@ C         /l\   /   \   /j\                                                    C
 C        /   \ /     \ /   \                                                   C
 C       /| o |o       o| o |\                                                  C
 C     \ j|/k\|      \  |/k\|l                                                  C
-C      \ /   \       \ /   \                                                   C
+C      \ /   \       \ /   \                                                   C 
 C       o     \       o     \                                                  C
 C       i             i                                                        C
-C                                                                              C
+C                                                                              C 
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C 4/7/01 AL Component s1 was removed, because it pertains to the respective