First trial to introduce steered MD non successful

diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN
index fe0297c..914c402 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CHAIN
+++ b/source/unres/src_MD-M/COMMON.CHAIN
@@ -2,7 +2,7 @@
      &  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
+     & prod,rt,dc_work,cref,crefjlee,chain_rep,dc_norm2,velAFMconst
       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 +20,5 @@
      & 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
+      common /afm/  forceAFMconst, distafminit,afmend,afmbeg,velAFMconst
 

diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL
index f9184f8..b8a775e 100644 (file)
--- a/source/unres/src_MD-M/COMMON.CONTROL
+++ b/source/unres/src_MD-M/COMMON.CONTROL
@@ -1,5 +1,5 @@
       integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad,
-     & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog
+     & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide
       logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec,
      &                 sideadd,lsecondary,read_cart,unres_pdb,
      &                 vdisulf,searchsc,lmuca,dccart,extconf,out1file,
@@ -8,6 +8,7 @@
      & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint,
      & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb
      & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file,
-     & constr_dist,gnorm_check,gradout,split_ene,symetr,AFMlog
+     & constr_dist,gnorm_check,gradout,split_ene,symetr,AFMlog,
+     & selfguide
 C... minim = .true. means DO minimization.
 C... energy_dec = .true. means print energy decomposition matrix

diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F
index bb87090..42b3af4 100644 (file)
--- a/source/unres/src_MD-M/MD_A-MTS.F
+++ b/source/unres/src_MD-M/MD_A-MTS.F
@@ -1798,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
@@ -1812,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

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 bfb0b65..2bd6f22 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -10273,39 +10273,3 @@ C AFM soubroutine for constant force
 C      print *,'AFM',Eafmforce
       return
       end
-C AFM soubroutine for constant 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'
-      include 'COMMON.MD'
-      real*8 diffafm(3)
-      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=-(dist-distafminit)
-      do i=1,3
-      gradafm(i,afmend-1)=-(velconst*diffafm(i)/dist-d_t(i,afmend-1))
-     & /d_time
-      gradafm(i,afmbeg-1)=(velconst*diffafm(i)/dist-d_t(i,afmbeg-1))
-     & /d_time
-      enddo
-C      print *,'AFM',Eafmforce
-      return
-      end
-      

diff --git a/source/unres/src_MD-M/geomout.F b/source/unres/src_MD-M/geomout.F
index 416122b..f8e8f41 100644 (file)
--- a/source/unres/src_MD-M/geomout.F
+++ b/source/unres/src_MD-M/geomout.F
@@ -454,14 +454,33 @@ c-----------------------------------------------------------------
      &          potEcomp(23),me
           format1="a133"
           print *,'A CHUJ',potEcomp(23)
-        else
           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
+       else if (SELFGUIDE.gt.0) then
+       distance=0.0
+       do j=1,3
+       distance=distance+c(j,afmend)-c(j,afmbeg)
+       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,$)')
+     &          itime,totT,EK,potE,totE,
+     &          rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),
+     &          distance,me
+          format1="a133"
+          print *,'A CHUJ',potEcomp(23)
+          write (line1,'(i10,f15.2,7f12.3,i5,$)')
+     &           itime,totT,EK,potE,totE,
+     &           kinetic_T,t_bath,gyrate(),
+     &           distance,me
+          format1="a114"
+        endif
+
        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 491c676..106458e 100644 (file)
--- a/source/unres/src_MD-M/int_to_cart.f
+++ b/source/unres/src_MD-M/int_to_cart.f
@@ -14,6 +14,7 @@ c-------------------------------------------------------------
       include 'COMMON.MD'
       include 'COMMON.IOUNITS'
       include 'COMMON.SCCOR' 
+      include 'COMMON.CONTROL'
 c   calculating dE/ddc1     
 C       print *,"wchodze22",ialph(2,1) 
        if (nres.lt.3) go to 18
@@ -275,6 +276,13 @@ c  Settind dE/ddnres
         enddo
        endif
 c   The side-chain vector derivatives
+      if (SELFGUIDE.gt.0) then
+      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
       return
       end      
        

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index 85e9ee3..6a6efd9 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -138,7 +138,8 @@ C Set up the time limit (caution! The time must be input in minutes!)
       indpdb=index(controlcard,'PDBSTART')
       extconf=(index(controlcard,'EXTCONF').gt.0)
       AFMlog=(index(controlcard,'AFM'))
-      print *,'AFMlog',AFMlog
+      selfguide=(index(controlcard,'SELFGUIDE'))
+      print *,'AFMlog',AFMlog,selfguide,"KUPA"
       call readi(controlcard,'IPRINT',iprint,0)
       call readi(controlcard,'MAXGEN',maxgen,10000)
       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
@@ -957,7 +958,7 @@ c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
         call flush(iout)
         if (constr_dist.gt.0) call read_dist_constr
         write (iout,*) "After read_dist_constr nhpb",nhpb
-        if (AFMlog.gt.0) call read_afminp
+        if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
         call hpb_partition
         if(me.eq.king.or..not.out1file)
      &   write (iout,*) 'Contact order:',co
@@ -2274,6 +2275,7 @@ C---------------------------------------------------------------------------
       call readi(afmcard,"BEG",afmbeg,0)
       call readi(afmcard,"END",afmend,0)
       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
+      call reada(afmcard,"VEL",velAFMconst,0.0d0)
       print *,'FORCE=' ,forceAFMconst
 CCCC NOW PROPERTIES FOR AFM
        distafminit=0.0d0