From: Adam Sieradzan Date: Fri, 17 Apr 2015 17:44:02 +0000 (+0200) Subject: Zaimplementowany AFM X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=84bd5165b1f28ca095851e740cbcc1ff2a6be29c Zaimplementowany AFM --- diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index 6fb7b6f..fe0297c 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -1,5 +1,6 @@ integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc, - & nres0,nstart_seq,chain_length,iprzes,tabperm,nperm + & 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 common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2), @@ -19,4 +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 diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index 40346c7..f9184f8 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 + & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec, & sideadd,lsecondary,read_cart,unres_pdb, & vdisulf,searchsc,lmuca,dccart,extconf,out1file, @@ -8,6 +8,6 @@ & 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 + & constr_dist,gnorm_check,gradout,split_ene,symetr,AFMlog C... minim = .true. means DO minimization. C... energy_dec = .true. means print energy decomposition matrix diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV index 6a5e74d..8082b0a 100644 --- a/source/unres/src_MD-M/COMMON.DERIV +++ b/source/unres/src_MD-M/COMMON.DERIV @@ -10,6 +10,7 @@ & gvdwpp(3,-1:maxres),gvdwc_scpp(3,-1:maxres), & gliptranc(3,-1:maxres), & gliptranx(3,-1:maxres), + & gradafm(3,-1:maxres), & gradx_scp(3,-1:maxres),gvdwc_scp(3,-1:maxres), & ghpbx(3,-1:maxres), & ghpbc(3,-1:maxres),gloc(maxvar,2),gradcorr(3,-1:maxres), diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index fa5cf66..da975d3 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -95,7 +95,7 @@ C Max. number of conformations in the pool parameter (max_pool=10) C Number of energy components integer n_ene,n_ene2 - parameter (n_ene=22,n_ene2=2*n_ene) + parameter (n_ene=23,n_ene2=2*n_ene) C Number of threads in deformation integer max_thread,max_thread2 parameter (max_thread=4,max_thread2=2*max_thread) diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 2dc00c1..bb87090 100644 --- 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) + call enerprint(potEcomp) +C print *,itime,'AFM',Eafmforc,etot + endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,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 0044ec5..d00102c 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -273,6 +273,9 @@ C print *,"przed lipidami" call Eliptransfer(eliptran) endif C print *,"za lipidami" + if (AFMlog.gt.0) then + call AFMforce(Eafmforce) + endif #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -315,6 +318,7 @@ C energia(20)=Uconst+Uconst_back energia(21)=esccor energia(22)=eliptran + energia(23)=Eafmforce c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" @@ -407,13 +411,14 @@ cMS$ATTRIBUTES C :: proc_proc Uconst=energia(20) esccor=energia(21) eliptran=energia(22) + Eafmforce=energia(23) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran + & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc @@ -421,6 +426,7 @@ cMS$ATTRIBUTES C :: proc_proc & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran + & +Eafmforce #endif energia(0)=etot c detecting NaNQ @@ -524,6 +530,7 @@ c enddo & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) + & +gradafm(j,i) enddo enddo @@ -541,6 +548,8 @@ c enddo & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) & +wliptran*gliptranc(j,i) + & +gradafm(j,i) + enddo enddo #endif @@ -677,6 +686,7 @@ c enddo & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) & +wliptran*gliptranc(j,i) + & +gradafm(j,i) #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ @@ -697,6 +707,8 @@ c enddo & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) & +wliptran*gliptranc(j,i) + & +gradafm(j,i) + #endif gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ @@ -993,6 +1005,7 @@ C------------------------------------------------------------------------ Uconst=energia(20) esccor=energia(21) eliptran=energia(22) + Eafmforce=energia(23) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, @@ -1001,7 +1014,7 @@ C------------------------------------------------------------------------ & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, & edihcnstr,ebr*nss, - & Uconst,eliptran,wliptran,etot + & Uconst,eliptran,wliptran,Eafmforce,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1026,7 +1039,9 @@ C------------------------------------------------------------------------ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ + & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETOT= ',1pE16.6,' (total)') + #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec, & estr,wbond,ebe,wang, @@ -1034,7 +1049,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, - & ebr*nss,Uconst,eliptran,wliptran,etot + & ebr*nss,Uconst,eliptran,wliptran,Eafmforc,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1058,6 +1073,7 @@ C------------------------------------------------------------------------ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ + & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -10224,3 +10240,37 @@ C eliptran=elpitran+0.0 ! I am in water enddo return end +C--------------------------------------------------------- +C AFM soubroutine for constant force + subroutine AFMforce(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) + 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=-forceAFMconst*(dist-distafminit) + do i=1,3 + gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist + gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist + enddo +C print *,'AFM',Eafmforce + return + end + diff --git a/source/unres/src_MD-M/gradient_p.F b/source/unres/src_MD-M/gradient_p.F index 3d90077..acd4472 100644 --- a/source/unres/src_MD-M/gradient_p.F +++ b/source/unres/src_MD-M/gradient_p.F @@ -382,7 +382,8 @@ C gscloc(j,i)=0.0d0 gsclocx(j,i)=0.0d0 gliptranc(j,i)=0.0d0 - gliptranx(j,i)=0.0d0 + gliptranx(j,i)=0.0d0 + gradafm(j,i)=0.0d0 do intertyp=1,3 gloc_sc(intertyp,i,icg)=0.0d0 enddo diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F index 5da38bd..03e928d 100644 --- a/source/unres/src_MD-M/readpdb.F +++ b/source/unres/src_MD-M/readpdb.F @@ -369,7 +369,6 @@ cc enddiag hfrag(i,j)=hfrag(i,j)-ishift enddo enddo - return end c--------------------------------------------------------------------------- diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index d02b83f..85e9ee3 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -137,6 +137,8 @@ C Set up the time limit (caution! The time must be input in minutes!) refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0) indpdb=index(controlcard,'PDBSTART') extconf=(index(controlcard,'EXTCONF').gt.0) + AFMlog=(index(controlcard,'AFM')) + print *,'AFMlog',AFMlog call readi(controlcard,'IPRINT',iprint,0) call readi(controlcard,'MAXGEN',maxgen,10000) call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) @@ -955,6 +957,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 call hpb_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co @@ -2253,6 +2256,35 @@ c------------------------------------------------------------------------------- enddo return end +C--------------------------------------------------------------------------- + subroutine read_afminp + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.SBRIDGE' + character*320 afmcard + print *, "wchodze" + call card_concat(afmcard) + call readi(afmcard,"BEG",afmbeg,0) + call readi(afmcard,"END",afmend,0) + call reada(afmcard,"FORCE",forceAFMconst,0.0d0) + print *,'FORCE=' ,forceAFMconst +CCCC NOW PROPERTIES FOR AFM + distafminit=0.0d0 + do i=1,3 + distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit + enddo + distafminit=dsqrt(distafminit) + print *,'initdist',distafminit + return + end + c------------------------------------------------------------------------------- subroutine read_dist_constr implicit real*8 (a-h,o-z)