From f038962fbb96e0c2c6f1ccb910373ceddd5b387b Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Fri, 31 Jul 2015 09:36:14 +0200 Subject: [PATCH] First trial to introduce steered MD non successful --- source/unres/src_MD-M/COMMON.CHAIN | 4 +-- source/unres/src_MD-M/COMMON.CONTROL | 5 ++-- source/unres/src_MD-M/MD_A-MTS.F | 19 +++++++++++++- source/unres/src_MD-M/energy_p_new_barrier.F | 36 -------------------------- source/unres/src_MD-M/geomout.F | 23 ++++++++++++++-- source/unres/src_MD-M/int_to_cart.f | 8 ++++++ source/unres/src_MD-M/readrtns_CSA.F | 6 +++-- 7 files changed, 56 insertions(+), 45 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index fe0297c..914c402 100644 --- 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 --- 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 --- 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 --- 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 --- 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 --- 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 --- 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 -- 1.7.9.5