X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=bfb0b65d541a4f7af22dcb3d8dbb1490c73c2056;hb=0ef28a0babbfafd06d3977c622ecbe98a5f41e86;hp=cbaadfbe8a46f394de092b39d7dadb240c747c20;hpb=845b24e5c8d81526cc47f5f95982661aa6eddcfa;p=unres.git 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 cbaadfb..bfb0b65 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 @@ -1623,7 +1639,7 @@ C what fraction I am in C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then + elseif (zj.gt.bufliptop) then fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick @@ -1825,12 +1841,12 @@ C the energy transfer exist if (zi.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- - & ((positi-bordlipbot)/lipbufthick) + & ((zi-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipi=sscalelip(fracinbuf) ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick) sslipi=sscalelip(fracinbuf) ssgradlipi=sscagradlip(fracinbuf)/lipbufthick else @@ -1893,12 +1909,12 @@ C the energy transfer exist if (zj.lt.buflipbot) then C what fraction I am in fracinbuf=1.0d0- - & ((positi-bordlipbot)/lipbufthick) + & ((zj-bordlipbot)/lipbufthick) C lipbufthick is thickenes of lipid buffore sslipj=sscalelip(fracinbuf) ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick - elseif (zi.gt.bufliptop) then - fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) sslipj=sscalelip(fracinbuf) ssgradlipj=sscagradlip(fracinbuf)/lipbufthick else @@ -2700,12 +2716,26 @@ c write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) c write(iout,*) 'b1=',b1(1,i-2) c write (iout,*) 'theta=', theta(i-1) enddo +#else + b1(1,i-2)=b(3,iti) + b1(2,i-2)=b(5,iti) + b2(1,i-2)=b(2,iti) + b2(2,i-2)=b(4,iti) + b1tilde(1,i-2)=b1(1,i-2) + b1tilde(2,i-2)=-b1(2,i-2) + b2tilde(1,i-2)=b2(1,i-2) + b2tilde(2,i-2)=-b2(2,i-2) + EE(1,2,i-2)=eeold(1,2,iti) + EE(2,1,i-2)=eeold(2,1,iti) + EE(2,2,i-2)=eeold(2,2,iti) + EE(1,1,i-2)=eeold(1,1,iti) + enddo +#endif #ifdef PARMAT do i=ivec_start+2,ivec_end+2 #else do i=3,nres+1 #endif -#endif if (i .lt. nres+1) then sin1=dsin(phi(i)) cos1=dcos(phi(i)) @@ -10141,7 +10171,8 @@ C lipbufthick is thickenes of lipid buffore gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran -C print *,"doing sccale for lower part" +C print *,"doing sccale for lower part" +C print *,i,sslip,fracinbuf,ssgradlip elseif (positi.gt.bufliptop) then fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) sslip=sscalelip(fracinbuf) @@ -10151,6 +10182,7 @@ C print *,"doing sccale for lower part" gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0 C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran C print *, "doing sscalefor top part" +C print *,i,sslip,fracinbuf,ssgradlip else eliptran=eliptran+pepliptran C print *,"I am in true lipid" @@ -10163,7 +10195,7 @@ C print *, "nic nie bylo w lipidzie?" C now multiply all by the peptide group transfer factor C eliptran=eliptran*pepliptran C now the same for side chains -C do i=1,1 +CV do i=1,1 do i=ilip_start,ilip_end if (itype(i).eq.ntyp1) cycle positi=(mod(c(3,i+nres),boxzsize)) @@ -10183,9 +10215,9 @@ C lipbufthick is thickenes of lipid buffore ssgradlip=-sscagradlip(fracinbuf)/lipbufthick eliptran=eliptran+sslip*liptranene(itype(i)) gliptranx(3,i)=gliptranx(3,i) - &+ssgradlip*liptranene(itype(i))/2.0d0 + &+ssgradlip*liptranene(itype(i)) gliptranc(3,i-1)= gliptranc(3,i-1) - &+ssgradlip*liptranene(itype(i))/2.0d0 + &+ssgradlip*liptranene(itype(i)) C print *,"doing sccale for lower part" elseif (positi.gt.bufliptop) then fracinbuf=1.0d0- @@ -10194,9 +10226,9 @@ C print *,"doing sccale for lower part" ssgradlip=sscagradlip(fracinbuf)/lipbufthick eliptran=eliptran+sslip*liptranene(itype(i)) gliptranx(3,i)=gliptranx(3,i) - &+ssgradlip*liptranene(itype(i))/2.0d0 + &+ssgradlip*liptranene(itype(i)) gliptranc(3,i-1)= gliptranc(3,i-1) - &+ssgradlip*liptranene(itype(i))/2.0d0 + &+ssgradlip*liptranene(itype(i)) C print *, "doing sscalefor top part",sslip,fracinbuf else eliptran=eliptran+liptranene(itype(i)) @@ -10208,3 +10240,72 @@ 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 +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 +