update

[unres.git] / source / unres / src_MD-M / MD.F

diff --git a/source/unres/src_MD-M/MD.F b/source/unres/src_MD-M/MD.F
index 704947a..0f12f9b 100644 (file)
--- a/source/unres/src_MD-M/MD.F
+++ b/source/unres/src_MD-M/MD.F
@@ -176,7 +176,10 @@ c Variable time step algorithm.
           call brown_step(itime)
         endif
         if (ntwe.ne.0) then
-         if (mod(itime,ntwe).eq.0) call statout(itime)
+         if (mod(itime,ntwe).eq.0) then
+           call statout(itime)
+           call returnbox
+          endif
 #ifdef VOUT
         do j=1,3
           v_work(j)=d_t(j,0)
@@ -189,7 +192,7 @@ c Variable time step algorithm.
           enddo
         enddo
         do i=nnt,nct
-          if (itype(i).ne.10 .and. itype(i).ne.21) then
+          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
             do j=1,3
               ind=ind+1
               v_work(ind)=d_t(j,i+nres)
@@ -291,7 +294,7 @@ c------------------------------------------------
       double precision difftol /1.0d-5/
       nbond=nct-nnt
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) nbond=nbond+1
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) nbond=nbond+1
       enddo
 c
       if (lprn1) then
@@ -313,7 +316,7 @@ c
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           ind1=ind1+1
           do j=1,3
             Bmat(ind+j,ind1)=dC_norm(j,i+nres)
@@ -390,7 +393,7 @@ c
           Td(i)=Td(i)+vbl*Tmat(i,ind)
         enddo
         do k=nnt,nct
-          if (itype(k).ne.10 .and. itype(i).ne.21) then
+          if (itype(k).ne.10 .and. itype(i).ne.ntyp1) then
             ind=ind+1
             Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
           endif
@@ -423,7 +426,7 @@ c
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           do j=1,3
             ind=ind+1
             zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
@@ -494,7 +497,7 @@ c
      &        i,(dC(j,i),j=1,3),xx
         enddo
         do i=nnt,nct
-          if (itype(i).ne.10 .and. itype(i).ne.21) then
+          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
             ind=ind+1
             xx=vbld(i+nres)-vbldsc0(1,itype(i))
             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
@@ -522,7 +525,7 @@ c      do iter=1,maxiter
         endif
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           ind=ind+1
           blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
           ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
@@ -565,7 +568,7 @@ c      do iter=1,maxiter
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           do j=1,3
             dc(j,i+nres)=zapas(ind+j)
             dc_work(ind+j)=zapas(ind+j)
@@ -609,7 +612,7 @@ c   Building the chain from the newly calculated coordinates
      &        i,(dC(j,i),j=1,3),xx
         enddo
         do i=nnt,nct
-          if (itype(i).ne.10 .and. itype(i).ne.21) then
+          if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
             ind=ind+1
             xx=vbld(i+nres)-vbldsc0(1,itype(i))
             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
@@ -1178,7 +1181,7 @@ c  forces).
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3
             d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time
@@ -1221,7 +1224,7 @@ c Applying velocity Verlet algorithm - step 1 to coordinates
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3    
             adt=d_a_old(j,inres)*d_time
@@ -1258,7 +1261,7 @@ c  Step 2 of the velocity Verlet algorithm: update velocities
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3
             d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time
@@ -1360,7 +1363,7 @@ c
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3    
             adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time
@@ -1436,7 +1439,7 @@ c
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3
             d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres)
@@ -1482,7 +1485,7 @@ c Side chains
         do j=1,3
           accel(j)=aux(j)
         enddo
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           do j=1,3 
             accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
           enddo
@@ -1527,7 +1530,7 @@ c            write (iout,*) "back",i,j,epdriftij
           enddo
         endif
 c Side chains
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           do j=1,3 
             epdriftij=
      &       dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i))
@@ -1574,7 +1577,7 @@ c      write(iout,*) "fact", fact
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3
             d_t(j,inres)=fact*d_t(j,inres)
@@ -1628,7 +1631,7 @@ c if the friction coefficients do not depend on surface area
           stdforcp(i)=stdfp*dsqrt(gamp)
         enddo
         do i=nnt,nct
-          stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i)))
+          stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(iabs(itype(i))))
         enddo
       endif
 c Open the pdb file for snapshotshots
@@ -1933,7 +1936,7 @@ c Transfer to the d_t vector
       do i=nnt,nct-1
         do j=1,3 
           ind=ind+1
-          if (itype(i).ne.21 .and. itype(i+1).ne.21) then
+          if (itype(i).ne.ntyp1 .and. itype(i+1).ne.ntyp1) then
             d_t(j,i)=d_t_work(ind)
           else
             d_t(j,i)=0.0d0
@@ -1941,7 +1944,7 @@ c Transfer to the d_t vector
         enddo
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           do j=1,3
             ind=ind+1
             d_t(j,i+nres)=d_t_work(ind)
@@ -2174,7 +2177,7 @@ c      enddo
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           do j=1,3
             dc_work(ind+j)=dc_old(j,i+nres)
             d_t_work(ind+j)=d_t_old(j,i+nres)
@@ -2222,7 +2225,7 @@ c      enddo
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3
             dc(j,inres)=dc_work(ind+j)
@@ -2283,7 +2286,7 @@ c
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3
             d_t(j,inres)=d_t_work(ind+j)
@@ -2440,7 +2443,7 @@ c      enddo
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           do j=1,3
             dc_work(ind+j)=dc_old(j,i+nres)
             d_t_work(ind+j)=d_t_old(j,i+nres)
@@ -2489,7 +2492,7 @@ c      enddo
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3
             dc(j,inres)=dc_work(ind+j)
@@ -2550,7 +2553,7 @@ c          ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j)
         ind=ind+3
       enddo
       do i=nnt,nct
-        if (itype(i).ne.10 .and. itype(i).ne.21) then
+        if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3
             d_t(j,inres)=d_t_work(ind+j)