From 8aafca27bccb6ad8deb5a7c8ead021d301ed4afc Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Thu, 6 Aug 2015 08:59:38 +0200 Subject: [PATCH] intoduction of constant velocity for spring --- source/unres/src_MD-M/COMMON.CHAIN | 7 ++-- source/unres/src_MD-M/MD_A-MTS.F | 36 +++++++++++--------- source/unres/src_MD-M/brown_step.F | 1 + source/unres/src_MD-M/energy_p_new_barrier.F | 46 ++++++++++++++++++++++++++ source/unres/src_MD-M/geomout.F | 22 +++++++----- source/unres/src_MD-M/int_to_cart.f | 12 +++---- source/unres/src_MD-M/readrtns_CSA.F | 1 + 7 files changed, 93 insertions(+), 32 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index 914c402..84ebf1b 100644 --- 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 --- 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 --- 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 --- 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 --- 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 --- 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 --- 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 -- 1.7.9.5