intoduction of constant velocity for spring

[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 0fe17b0..dc58df8 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)
+C           call enerprint(potEcomp)
+C           print *,itime,'AFM',Eafmforc,etot
+         endif
 #ifdef VOUT
         do j=1,3
           v_work(j)=d_t(j,0)
@@ -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)
@@ -513,6 +520,8 @@ c Second step of the velocity Verlet algorithm
           endif                    
           if (rattle) call rattle2
           totT=totT+d_time
+          totTafm=totT
+C          print *,totTafm,"TU?"
           if (d_time.ne.d_time0) then
             d_time=d_time0
 #ifndef   LANG0
@@ -910,6 +919,7 @@ c Compute the complete potential energy
       potE=potEcomp(0)-potEcomp(20)
 c      potE=energia_short(0)+energia_long(0)
       totT=totT+d_time
+      totTafm=totT
 c Calculate the kinetic and the total energy and the kinetic temperature
       call kinetic(EK)
       totE=EK+potE
@@ -996,6 +1006,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,6 +1017,7 @@ c Applying velocity Verlet algorithm - step 1 to coordinates
         enddo
       enddo
       do i=nnt,nct
+C        do i=0,nres
         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
           inres=i+nres
           do j=1,3    
@@ -1087,12 +1100,25 @@ c Applying velocity Verlet algorithm - step 1 to coordinates
 c
 c Compute friction and stochastic forces
 c
+#ifdef MPI
       time00=MPI_Wtime()
+#else
+      time00=tcpu()
+#endif
       call friction_force
+#ifdef MPI
       time_fric=time_fric+MPI_Wtime()-time00
       time00=MPI_Wtime()
+#else
+      time_fric=time_fric+tcpu()-time00
+      time00=tcpu()
+#endif
       call stochastic_force(stochforcvec) 
+#ifdef MPI
       time_stoch=time_stoch+MPI_Wtime()-time00
+#else
+      time_stoch=time_stoch+tcpu()-time00
+#endif
 c
 c Compute the acceleration due to friction forces (d_af_work) and stochastic
 c forces (d_as_work)
@@ -1510,11 +1536,13 @@ c        statname=statname(:ilen(statname)-5)//qstr(:ipos-1)//'.stat'
        if (restart1file) then
          if (me.eq.king)
      &     inquire(file=mremd_rst_name,exist=file_exist)
+#ifdef MPI
            write (*,*) me," Before broadcast: file_exist",file_exist
          call MPI_Bcast(file_exist,1,MPI_LOGICAL,king,CG_COMM,
      &          IERR)
          write (*,*) me," After broadcast: file_exist",file_exist
 c        inquire(file=mremd_rst_name,exist=file_exist)
+#endif
         if(me.eq.king.or..not.out1file)
      &   write(iout,*) "Initial state read by master and distributed"
        else
@@ -1544,6 +1572,7 @@ c        inquire(file=mremd_rst_name,exist=file_exist)
         endif
         call random_vel
         totT=0.0d0
+        totTafm=totT
        endif
       else
 c Generate initial velocities
@@ -1551,6 +1580,8 @@ c Generate initial velocities
      &   write(iout,*) "Initial velocities randomly generated"
         call random_vel
         totT=0.0d0
+CtotTafm is the variable for AFM time which eclipsed during  
+        totTafm=totT
       endif
 c      rest2name = prefix(:ilen(prefix))//'.rst'
       if(me.eq.king.or..not.out1file)then
@@ -1608,7 +1639,7 @@ c Removing the velocity of the center of mass
       call chainbuild_cart
       call kinetic(EK)
       if (tbf) then
-        call verlet_bath(EK)
+        call verlet_bath
       endif      
       kinetic_T=2.0d0/(dimen3*Rb)*EK
       if(me.eq.king.or..not.out1file)then
@@ -1773,7 +1804,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
@@ -1787,10 +1818,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
+C       if (SELFGUIDE.gt.0) then
+C       distance=0.0
+C       do j=1,3
+C       vec_afm(j)=c(j,afmend)-c(j,afmbeg)  
+C       distance=distance+vec_afm(j)**2
+C       enddo
+C       distance=dsqrt(distance)
+C       do j=1,3
+C         d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance
+C         d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance
+C         write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3),
+C     &    d_t_work_new(j+(afmend-1)*3)
+C       enddo
+
+C       endif
+
 c diagnostics
 c      Ek1=0.0d0
 c      ii=0