dynamic disulfides are working again in md

[unres.git] / source / unres / src_MD / ssMD.F

diff --git a/source/unres/src_MD/ssMD.F b/source/unres/src_MD/ssMD.F
index b425056..dd210e5 100644 (file)
--- a/source/unres/src_MD/ssMD.F
+++ b/source/unres/src_MD/ssMD.F
@@ -400,10 +400,10 @@ c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
       if (havebond) then
 #ifndef CLUST
 #ifndef WHAM
-        if (dyn_ssbond_ij(i,j).eq.1.0d300) then
-          write(iout,'(a15,f12.2,f8.1,2i5)')
-     &         "SSBOND_E_FORM",totT,t_bath,i,j
-        endif
+c        if (dyn_ssbond_ij(i,j).eq.1.0d300) then
+c          write(iout,'(a15,f12.2,f8.1,2i5)')
+c     &         "SSBOND_E_FORM",totT,t_bath,i,j
+c        endif
 #endif
 #endif
         dyn_ssbond_ij(i,j)=eij
@@ -411,8 +411,8 @@ c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
         dyn_ssbond_ij(i,j)=1.0d300
 #ifndef CLUST
 #ifndef WHAM
-        write(iout,'(a15,f12.2,f8.1,2i5)')
-     &       "SSBOND_E_BREAK",totT,t_bath,i,j
+c        write(iout,'(a15,f12.2,f8.1,2i5)')
+c     &       "SSBOND_E_BREAK",totT,t_bath,i,j
 #endif
 #endif
       endif
@@ -449,10 +449,15 @@ c-------END TESTING CODE
      &       +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
      &       +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
       enddo
-      do k=i,j-1
-        do l=1,3
-          gvdwc(l,k)=gvdwc(l,k)+gg(l)
-        enddo
+cgrad      do k=i,j-1
+cgrad        do l=1,3
+cgrad          gvdwc(l,k)=gvdwc(l,k)+gg(l)
+cgrad        enddo
+cgrad      enddo
+
+      do l=1,3
+        gvdwc(l,i)=gvdwc(l,i)-gg(l)
+        gvdwc(l,j)=gvdwc(l,j)+gg(l)
       enddo
 
       return
@@ -603,6 +608,9 @@ cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
 #endif
       enddo
 
+c CRC this part of code is not clean, 
+c dyn_ss will not work with contrains
+
       if (diff.gt.0) then
         do i=1,diff
           ihpb(nhpb+i)=ihpb(nss+i)
@@ -612,10 +620,12 @@ cmc      write(iout,*)"NEWNSS ",newnss,(newihpb(i),newjhpb(i),i=1,newnss)
         enddo
       else if (diff.lt.0) then
         do i=diff,-1
+         if(nss+i.gt.0.and.nhpb+i.gt.0) then
           ihpb(nss+i)=ihpb(nhpb+i)
           jhpb(nss+i)=jhpb(nhpb+i)
           forcon(nss+i)=forcon(nhpb+i)
           dhpb(nss+i)=dhpb(nhpb+i)
+         endif
         enddo
       endif