First trial to introduce steered MD non successful

[unres.git] / source / unres / src_MD-M / MD_A-MTS.F

diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F
index 56b3ea8..42b3af4 100644 (file)
--- a/source/unres/src_MD-M/MD_A-MTS.F
+++ b/source/unres/src_MD-M/MD_A-MTS.F
@@ -196,7 +196,11 @@ c Variable time step algorithm.
 #endif
         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 enerprint(potEcomp)
+C           print *,itime,'AFM',Eafmforc,etot
+         endif
 #ifdef VOUT
         do j=1,3
           v_work(j)=d_t(j,0)
@@ -209,7 +213,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)
@@ -230,6 +234,9 @@ c Variable time step algorithm.
 #endif
         endif
         if (mod(itime,ntwx).eq.0) then
+          write(iout,*) 'time=',itime
+C          call check_ecartint
+          call returnbox
           write (tytul,'("time",f8.2)') totT
           if(mdpdb) then
              call hairpin(.true.,nharp,iharp)
@@ -955,7 +962,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
@@ -996,6 +1003,8 @@ c Applying velocity Verlet algorithm - step 1 to coordinates
         d_t(j,0)=d_t_old(j,0)+adt
       enddo
       do i=nnt,nct-1   
+C      SPYTAC ADAMA
+C       do i=0,nres
         do j=1,3    
           adt=d_a_old(j,i)*d_time
           adt2=0.5d0*adt
@@ -1005,7 +1014,8 @@ 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
+C        do i=0,nres
+        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
@@ -1049,7 +1059,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
@@ -1161,7 +1171,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
@@ -1226,7 +1236,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)
@@ -1273,6 +1283,7 @@ c            accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i))
 c            if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
             if (dabs(accel(j)).gt.dabs(accel_old(j))) then
               dacc=dabs(accel(j)-accel_old(j))
+c              write (iout,*) i,dacc
               if (dacc.gt.amax) amax=dacc
             endif
           enddo
@@ -1291,7 +1302,7 @@ c        accel(j)=aux(j)
         enddo
       endif
       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 
 c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
             accel_old(j)=accel_old(j)+d_a_old(j,i+nres)
@@ -1302,6 +1313,7 @@ c            accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres)
 c          if (dabs(accel(j)).gt.amax) amax=dabs(accel(j))
           if (dabs(accel(j)).gt.dabs(accel_old(j))) then
             dacc=dabs(accel(j)-accel_old(j))
+c            write (iout,*) "side-chain",i,dacc
             if (dacc.gt.amax) amax=dacc
           endif
         enddo
@@ -1344,7 +1356,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))
@@ -1391,7 +1403,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)
@@ -1446,7 +1458,8 @@ 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(iabs(itype(i)))
+     &                *dsqrt(gamsc(iabs(itype(i))))
         enddo
       endif
 c Open the pdb file for snapshotshots
@@ -1785,7 +1798,7 @@ c-----------------------------------------------------------
       include 'COMMON.IOUNITS'
       include 'COMMON.NAMES'
       include 'COMMON.TIME1'
-      double precision xv,sigv,lowb,highb
+      double precision xv,sigv,lowb,highb,vec_afm(3)
 c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m 
 c First generate velocities in the eigenspace of the G matrix
 c      write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3
@@ -1799,10 +1812,27 @@ c      call flush(iout)
           lowb=-5*sigv
           highb=5*sigv
           d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb)
+
 c          write (iout,*) "i",i," ii",ii," geigen",geigen(i),
 c     &      " d_t_work_new",d_t_work_new(ii)
         enddo
       enddo
+       if (SELFGUIDE.gt.0) then
+       distance=0.0
+       do j=1,3
+       vec_afm(j)=c(j,afmend)-c(j,afmbeg)  
+       distance=distance+vec_afm(j)**2
+       enddo
+       distance=dsqrt(distance)
+       do j=1,3
+         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
+         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
+         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
+     &    d_t_work_new(j+(afmend-1)*3)
+       enddo
+
+       endif
+
 c diagnostics
 c      Ek1=0.0d0
 c      ii=0
@@ -1839,7 +1869,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)
@@ -2068,7 +2098,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)
@@ -2177,7 +2207,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)
@@ -2334,7 +2364,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)
@@ -2383,7 +2413,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)
@@ -2444,7 +2474,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)