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),
& 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
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,
& 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
& 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),
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)
#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)
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
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"
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
& +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
& wturn6*gcorr6_turn_long(j,i)+
& wstrain*ghpbc(j,i)
& +wliptran*gliptranc(j,i)
+ & +gradafm(j,i)
enddo
enddo
& wturn6*gcorr6_turn_long(j,i)+
& wstrain*ghpbc(j,i)
& +wliptran*gliptranc(j,i)
+ & +gradafm(j,i)
+
enddo
enddo
#endif
& 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)+
& 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)+
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,
& 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)'/
& '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,
& 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)'/
& '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
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
+
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
hfrag(i,j)=hfrag(i,j)-ishift
enddo
enddo
-
return
end
c---------------------------------------------------------------------------
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)
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
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)
projects / unres.git / commitdiff