intoduction of constant velocity for spring

diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN
index 914c402..84ebf1b 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CHAIN
+++ b/source/unres/src_MD-M/COMMON.CHAIN
@@ -2,7 +2,8 @@
      &  nres0,nstart_seq,chain_length,iprzes,tabperm,nperm,afmend,
      &  afmbeg
       double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r,
-     & prod,rt,dc_work,cref,crefjlee,chain_rep,dc_norm2,velAFMconst
+     & prod,rt,dc_work,cref,crefjlee,chain_rep,dc_norm2,velAFMconst,
+     & totTafm
       common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
      & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
      & dc_norm2(3,0:maxres2),
@@ -20,5 +21,7 @@
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
       common /box/  boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
      & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
-      common /afm/  forceAFMconst, distafminit,afmend,afmbeg,velAFMconst
+      common /afm/  forceAFMconst, distafminit,afmend,afmbeg,
+     & velAFMconst,
+     & totTafm
 

diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F
index 42b3af4..dc58df8 100644 (file)
--- a/source/unres/src_MD-M/MD_A-MTS.F
+++ b/source/unres/src_MD-M/MD_A-MTS.F
@@ -198,7 +198,7 @@ c Variable time step algorithm.
         if (ntwe.ne.0) then
          if (mod(itime,ntwe).eq.0) then
            call statout(itime)
-           call enerprint(potEcomp)
+C           call enerprint(potEcomp)
 C           print *,itime,'AFM',Eafmforc,etot
          endif
 #ifdef VOUT
@@ -520,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
@@ -917,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
@@ -1569,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
@@ -1576,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
@@ -1817,21 +1823,21 @@ 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
+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
 
-       endif
+C       endif
 
 c diagnostics
 c      Ek1=0.0d0

diff --git a/source/unres/src_MD-M/brown_step.F b/source/unres/src_MD-M/brown_step.F
index 0be97f5..49652b8 100644 (file)
--- a/source/unres/src_MD-M/brown_step.F
+++ b/source/unres/src_MD-M/brown_step.F
@@ -381,6 +381,7 @@ c Calculate energy and forces
       potE=potEcomp(0)-potEcomp(20)
       call cartgrad
       totT=totT+d_time
+      totTafm=totT
 c  Calculate the kinetic and total energy and the kinetic temperature
       call kinetic(EK)
 #ifdef MPI

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 2bd6f22..2208b17 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -275,6 +275,8 @@ C      print *,"przed lipidami"
 C      print *,"za lipidami"
       if (AFMlog.gt.0) then
         call AFMforce(Eafmforce)
+      else if (selfguide.gt.0) then
+        call AFMvel(Eafmforce)
       endif
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
@@ -10273,3 +10275,47 @@ C AFM soubroutine for constant force
 C      print *,'AFM',Eafmforce
       return
       end
+C---------------------------------------------------------
+C AFM subroutine with pseudoconstant velocity
+       subroutine AFMvel(Eafmforce)
+       implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CALC'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
+      include 'COMMON.SBRIDGE'
+      real*8 diffafm(3)
+C Only for check grad COMMENT if not used for checkgrad
+C      totT=3.0d0
+C--------------------------------------------------------
+C      print *,"wchodze"
+      dist=0.0d0
+      Eafmforce=0.0d0
+      do i=1,3
+      diffafm(i)=c(i,afmend)-c(i,afmbeg)
+      dist=dist+diffafm(i)**2
+      enddo
+      dist=dsqrt(dist)
+      Eafmforce=0.5d0*forceAFMconst
+     & *(distafminit+totTafm*velAFMconst-dist)**2
+C      Eafmforce=-forceAFMconst*(dist-distafminit)
+      do i=1,3
+      gradafm(i,afmend-1)=-forceAFMconst*
+     &(distafminit+totTafm*velAFMconst-dist)
+     &*diffafm(i)/dist
+      gradafm(i,afmbeg-1)=forceAFMconst*
+     &(distafminit+totTafm*velAFMconst-dist)
+     &*diffafm(i)/dist
+      enddo
+C      print *,'AFM',Eafmforce,totTafm*velAFMconst,dist
+      return
+      end
+

diff --git a/source/unres/src_MD-M/geomout.F b/source/unres/src_MD-M/geomout.F
index f8e8f41..0b711bd 100644 (file)
--- a/source/unres/src_MD-M/geomout.F
+++ b/source/unres/src_MD-M/geomout.F
@@ -453,34 +453,38 @@ c-----------------------------------------------------------------
      &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
      &          potEcomp(23),me
           format1="a133"
-          print *,'A CHUJ',potEcomp(23)
+         else
+C          print *,'A CHUJ',potEcomp(23)
           write (line1,'(i10,f15.2,7f12.3,i5,$)')
      &           itime,totT,EK,potE,totE,
      &           kinetic_T,t_bath,gyrate(),
      &           potEcomp(23),me
           format1="a114"
         endif
-       else if (SELFGUIDE.gt.0) then
+       else if (selfguide.gt.0) then
        distance=0.0
        do j=1,3
-       distance=distance+c(j,afmend)-c(j,afmbeg)
+       distance=distance+(c(j,afmend)-c(j,afmbeg))**2
        enddo
        distance=dsqrt(distance)
        if (refstr) then
          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
-          write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')
+          write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2,
+     &    f9.2,i5,$)')
      &          itime,totT,EK,potE,totE,
      &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
-     &          distance,me
+     &          distance,potEcomp(23),me
           format1="a133"
-          print *,'A CHUJ',potEcomp(23)
-          write (line1,'(i10,f15.2,7f12.3,i5,$)')
+C          print *,"CHUJOWO"
+         else
+C          print *,'A CHUJ',potEcomp(23)
+          write (line1,'(i10,f15.2,8f12.3,i5,$)')
      &           itime,totT,EK,potE,totE,
      &           kinetic_T,t_bath,gyrate(),
-     &           distance,me
+     &           distance,potEcomp(23),me
           format1="a114"
         endif
-
+       else
        if (refstr) then
          call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
           write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)')

diff --git a/source/unres/src_MD-M/int_to_cart.f b/source/unres/src_MD-M/int_to_cart.f
index 106458e..d3a8a92 100644 (file)
--- a/source/unres/src_MD-M/int_to_cart.f
+++ b/source/unres/src_MD-M/int_to_cart.f
@@ -276,13 +276,13 @@ c  Settind dE/ddnres
         enddo
        endif
 c   The side-chain vector derivatives
-      if (SELFGUIDE.gt.0) then
-      do j=1,3
+C      if (SELFGUIDE.gt.0) then
+C      do j=1,3
 C       gcart(j,afmbeg)=gcart(j,afmbeg)+gcart(j,afmend)
-       gcart(j,afmbeg)=0.0d0
-       gcart(j,afmend)=0.0d0
-      enddo
-      endif
+C       gcart(j,afmbeg)=0.0d0
+C       gcart(j,afmend)=0.0d0
+C      enddo
+C      endif
       return
       end      
        

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 6a6efd9..6ab2a4e 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -2202,6 +2202,7 @@ c-------------------------------------------------------------------------------
       include 'COMMON.MD'
       open(irest2,file=rest2name,status='unknown')
       read(irest2,*) totT,EK,potE,totE,t_bath
+      totTafm=totT
       do i=1,2*nres
          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
       enddo