From bde71ffff7527add147c1bee4bad23f02b4d1ab7 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Tue, 10 Jan 2017 14:37:10 +0100 Subject: [PATCH] dodanie parametrow do nanotube --- PARAM/tube.parm | 50 +- PARAM/tube_carbonano.parm | 25 + PARAM/tube_sigma.parm | 6 + source/unres/src_MD-M/COMMON.INTERACT | 6 +- source/unres/src_MD-M/COMMON.MD | 14 +- source/unres/src_MD-M/DIMENSIONS | 3 + source/unres/src_MD-M/MD_A-MTS.F | 949 +++++++++++++++++++++++++- source/unres/src_MD-M/energy_p_new_barrier.F | 72 +- source/unres/src_MD-M/readrtns_CSA.F | 45 ++ source/wham/src-M/energy_p_new.F | 146 ++++ 10 files changed, 1263 insertions(+), 53 deletions(-) create mode 100644 PARAM/tube_carbonano.parm create mode 100644 PARAM/tube_sigma.parm diff --git a/PARAM/tube.parm b/PARAM/tube.parm index 182c68f..d3f9616 100644 --- a/PARAM/tube.parm +++ b/PARAM/tube.parm @@ -1,26 +1,24 @@ - 200.00 2.20 3.0 - 1.00 0.50 3.0 - 1.00 0.20 3.0 - 1.00 0.30 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 200.00 2.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - 1.00 1.00 3.0 - + 6.2159898496 + 6.6717260508 + 6.6424306340 + 6.9715947250 + 6.9241225787 + 6.5291325353 + 6.1394821460 + 5.4415971840 + 5.2914044780 + 4.7881880526 + 4.1302408718 + 3.9405275117 + 3.9022731880 + 3.6078928145 + 3.0809713206 + 3.0182595342 + 4.7935507272 + 3.8711928134 + 3.6826168780 + 4.4872661030 + 6.6717260508 + 6.6424306340 + 5.2914044780 + 5.2914044780 diff --git a/PARAM/tube_carbonano.parm b/PARAM/tube_carbonano.parm new file mode 100644 index 0000000..f282790 --- /dev/null +++ b/PARAM/tube_carbonano.parm @@ -0,0 +1,25 @@ + 1.1970470 5.3667307 3.0000000 + 1.5539975 5.6438808 3.0000000 + 1.6679316 5.6689787 3.0000000 + 1.6606077 5.9381499 3.0000000 + 1.7428987 5.8625088 3.0000000 + 1.7310307 5.9950466 3.0000000 + 1.6322831 5.8318806 3.0000000 + 1.5348705 5.4955850 3.0000000 + 1.3603992 5.3937664 3.0000000 + 1.3228511 5.4371481 3.0000000 + 1.1970470 5.3667307 3.0000000 + 1.0325602 5.5439558 3.0000000 + 0.98513186 5.3780737 3.0000000 + 0.97556829 5.3995867 3.0000000 + 0.90197319 5.4184709 3.0000000 + 0.77024281 5.4679136 3.0000000 + 0.75456488 5.4686551 3.0000000 + 1.1983876 5.3466215 3.0000000 + 0.96779823 5.2968884 3.0000000 + 0.92065424 5.3752089 3.0000000 + 1.1218165 5.6721835 3.0000000 + 1.6679316 5.7029562 3.0000000 + 1.6606077 5.9355397 3.0000000 + 1.3228511 5.4343948 3.0000000 + 1.3228511 5.4343948 3.0000000 diff --git a/PARAM/tube_sigma.parm b/PARAM/tube_sigma.parm new file mode 100644 index 0000000..acefcf4 --- /dev/null +++ b/PARAM/tube_sigma.parm @@ -0,0 +1,6 @@ + 2.674806001700000 2.699903447923231 2.969074945285474 2.893433539204819 + 3.025971784817564 2.862805741160594 2.526510086781303 2.424691603147259 + 2.468072853823478 2.397655578798363 2.574880627216491 2.408998508198888 + 2.430511884949466 2.449395893015385 2.498838549850558 2.499580191788341 + 2.377546242648636 2.327813161765345 2.406133938779852 2.703108539059061 + 2.7338810145 2.9664647229 2.4653201213 2.4653201213 diff --git a/source/unres/src_MD-M/COMMON.INTERACT b/source/unres/src_MD-M/COMMON.INTERACT index 14d92ef..be18ffe 100644 --- a/source/unres/src_MD-M/COMMON.INTERACT +++ b/source/unres/src_MD-M/COMMON.INTERACT @@ -1,7 +1,11 @@ double precision aa,bb,augm,aad,bad,app,bpp,ale6,ael3,ael6, &aa_lip,bb_lip,aa_aq,bb_aq,sc_aa_tube_par,sc_bb_tube_par, & pep_aa_tube,pep_bb_tube - + double precision wdti,wdti2,wdti4,wdti8, + & wdtii,wdtii2,wdtii4,wdtii8 + common /nosehoover_dt/ + & wdti(maxyosh),wdti2(maxyosh),wdti4(maxyosh),wdti8(maxyosh), + & wdtii(maxyosh),wdtii2(maxyosh),wdtii4(maxyosh),wdtii8(maxyosh) integer expon,expon2 integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro, & ielstart,ielend,ielstart_vdw,ielend_vdw,nscp_gr,iscpstart, diff --git a/source/unres/src_MD-M/COMMON.MD b/source/unres/src_MD-M/COMMON.MD index b17c722..66cb890 100644 --- a/source/unres/src_MD-M/COMMON.MD +++ b/source/unres/src_MD-M/COMMON.MD @@ -30,7 +30,9 @@ & ifrag_back(3,maxfrag_back,maxprocs/20),ntime_split,ntime_split0, & maxtime_split logical large,print_compon,tbf,rest,reset_moment,reset_vel, - & surfarea,rattle,usampl,mdpdb,RESPA + & surfarea,rattle,usampl,mdpdb,RESPA,tnp,tnp1,tnh,xiresp + integer nresn,nyosh,nnos + double precision glogs,qmass,vlogs,xlogs integer igmult_start,igmult_end,my_ng_count,ng_start,ng_counts, & nginv_start,nginv_counts,myginv_ng_count common /back_constr/ uconst_back,utheta,ugamma,uscdiff, @@ -42,7 +44,7 @@ common /mdpar/ v_ini,d_time,d_time0,scal_fric, & t_bath,tau_bath,dvmax,damax,n_timestep,mdpdb, & ntime_split,ntime_split0,maxtime_split, - & ntwx,ntwe,large,print_compon,tbf,rest + & ntwx,ntwe,large,print_compon,tbf,rest,tnp,tnp1,tnh common /MDcalc/ totT,totE,potE,potEcomp,EK,amax,edriftmax, & kinetic_T common /lagrange/ d_t,d_t_old,d_t_new,d_t_work, @@ -63,3 +65,11 @@ & myginv_ng_count, & ng_start(0:MaxProcs-1),ng_counts(0:MaxProcs-1), & nginv_start(0:MaxProcs),nginv_counts(0:MaxProcs-1) + double precision pi_np,pistar,s_np,s12_np,Q_np,E_old,H0,E_long, + & sold_np,d_t_half,Csplit,hhh + common /nosepoincare/ pi_np,pistar,s_np,s12_np,Q_np,E_old,H0, + & E_long,sold_np,d_t_half(3,0:MAXRES2),Csplit,hhh + common /nosehoover/ glogs(maxmnh),qmass(maxmnh), + & vlogs(maxmnh),xlogs(maxmnh), + & nresn,nyosh,nnos,xiresp + diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index d6c5e46..b14b323 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -140,3 +140,6 @@ C Maximum number of conformation stored in cache on each CPU before sending C to master; depends on nstex / ntwx ratio integer max_cache_traj parameter (max_cache_traj=10) +C NOSE-HOOVER + integer maxmnh,maxyosh + parameter(maxmnh=10,maxyosh=5) diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 3a45231..6650c5a 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -323,6 +323,8 @@ c------------------------------------------------------------------------------- common /stochcalc/ stochforcvec integer itime logical scale + double precision HNose1,HNose,HNose_nh,H,vtnp(maxres6) + double precision vtnp_(maxres6),vtnp_a(maxres6) c scale=.true. icount_scale=0 @@ -366,7 +368,19 @@ c First step of the velocity Verlet algorithm #endif else if (lang.eq.1) then call sddir_verlet1 +C else if (tnp1) then +C call tnp1_step1 +C else if (tnp) then +C call tnp_step1 else + if (tnh) then + call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8) + do i=0,2*nres + do j=1,3 + d_t_old(j,i)=d_t_old(j,i)*scale_nh + enddo + enddo + endif call verlet1 endif c Build the chain from the newly calculated coordinates @@ -411,6 +425,7 @@ c Calculate energy and forces t_etotal=t_etotal+tcpu()-tt0 #endif #endif + E_old=potE potE=potEcomp(0)-potEcomp(20) call cartgrad c Get the new accelerations @@ -515,8 +530,21 @@ c Second step of the velocity Verlet algorithm #endif else if (lang.eq.1) then call sddir_verlet2 +C> else if (tnp1) then +C> call tnp1_step2 +C> else if (tnp) then +C> call tnp_step2 else call verlet2 + if (tnh) then + call kinetic(EK) + call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8) + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t(j,i)*scale_nh + enddo + enddo + endif endif if (rattle) call rattle2 totT=totT+d_time @@ -546,6 +574,15 @@ c Restore the matrices of tinker integrator if the time step has been restored scale=.false. endif enddo + if (tnp .or. tnp1) then + do i=0,2*nres + do j=1,3 + d_t_old(j,i)=d_t(j,i) + d_t(j,i)=d_t(j,i)/s_np + enddo + enddo + endif + c Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE @@ -562,12 +599,67 @@ c Backup the coordinates, velocities, and accelerations do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) - d_t_old(j,i)=d_t(j,i) + if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0 .and. large) then + if(tnp .or. tnp1) then + HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3) + H=(HNose1-H0)*s_np +cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0 +cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np) +cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0 + hhh=h + endif + + if(tnh) then + HNose1=Hnose_nh(EK,potE) + H=HNose1-H0 + hhh=h +cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0 + endif + + if (large) then + itnp=0 + do j=1,3 + itnp=itnp+1 + vtnp(itnp)=d_t(j,0) + enddo + do i=nnt,nct-1 + do j=1,3 + itnp=itnp+1 + vtnp(itnp)=d_t(j,i) + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + itnp=itnp+1 + vtnp(itnp)=d_t(j,inres) + enddo + endif + enddo + +c Transform velocities from UNRES coordinate space to cartesian and Gvec +c eigenvector space + + do i=1,dimen3 + vtnp_(i)=0.0d0 + vtnp_a(i)=0.0d0 + do j=1,dimen3 + vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j) + vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j) + enddo + vtnp_(i)=vtnp_(i)*dsqrt(geigen(i)) + enddo + + do i=1,dimen3 + write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i) + enddo + write (iout,*) "Velocities, step 2" do i=0,nres write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3), @@ -575,6 +667,7 @@ c Backup the coordinates, velocities, and accelerations enddo endif endif + endif return end c------------------------------------------------------------------------------- @@ -619,6 +712,7 @@ c------------------------------------------------------------------------------- double precision stochforcvec(MAXRES6) common /stochcalc/ stochforcvec integer itime + double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6) logical scale common /cipiszcze/ itt itt=itime @@ -644,7 +738,26 @@ c call cartprint c c Perform the initial RESPA step (increment velocities) c write (iout,*) "*********************** RESPA ini" - call RESPA_vel + if (tnp1) then + call tnp_respa_step1 + else if (tnp) then + call tnp_respa_step1 + else + if (tnh.and..not.xiresp) then + call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8) + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t(j,i)*scale_nh + enddo + enddo + endif + call RESPA_vel + endif + +cd if(tnp .or. tnp1) then +cd write (iout,'(a,3f)') "EE1 NP S, pi",totT, s_np, pi_np +cd endif + if (ntwe.ne.0) then if (mod(itime,ntwe).eq.0 .and. large) then write (iout,*) "Velocities, end" @@ -661,6 +774,7 @@ c Compute the short-range forces tt0 = tcpu() #endif C 7/2/2009 commented out + if (tnp.or.tnp1) potE=energia_short(0) c call zerograd c call etotal_short(energia_short) c call cartgrad @@ -689,7 +803,8 @@ C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) - d_t_old(j,i)=d_t(j,i) +C d_t_old(j,i)=d_t(j,i) + if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo @@ -741,7 +856,21 @@ c First step of the velocity Verlet algorithm #endif else if (lang.eq.1) then call sddir_verlet1 + else if (tnp1) then + call tnp1_respa_i_step1 + else if (tnp) then + call tnp_respa_i_step1 else + if (tnh.and.xiresp) then + call kinetic(EK) + call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8) + do i=0,2*nres + do j=1,3 + d_t_old(j,i)=d_t_old(j,i)*scale_nh + enddo + enddo +cd write(iout,*) "SSS1",itsplit,EK,scale_nh + endif call verlet1 endif c Build the chain from the newly calculated coordinates @@ -771,6 +900,9 @@ c Calculate energy and forces call etotal_short(energia_short) if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_short) + E_old=potE + potE=energia_short(0) + #ifdef TIMING_ENE #ifdef MPI t_eshort=t_eshort+MPI_Wtime()-tt0 @@ -847,15 +979,40 @@ c Second step of the velocity Verlet algorithm #endif else if (lang.eq.1) then call sddir_verlet2 + else if (tnp1) then + call tnp1_respa_i_step2 + else if (tnp) then + call tnp_respa_i_step2 else call verlet2 + if (tnh.and.xiresp) then + call kinetic(EK) + call nhcint(EK,scale_nh,wdtii,wdtii2,wdtii4,wdtii8) + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t(j,i)*scale_nh + enddo + enddo +cd write(iout,*) "SSS2",itsplit,EK,scale_nh + endif endif if (rattle) call rattle2 + if (tnp .or. tnp1) then + do i=0,2*nres + do j=1,3 + d_t_old(j,i)=d_t(j,i) + if (tnp) d_t(j,i)=d_t(j,i)/s_np + if (tnp1) d_t(j,i)=d_t(j,i)/s_np + enddo + enddo + endif + + c Backup the coordinates, velocities, and accelerations do i=0,2*nres do j=1,3 dc_old(j,i)=dc(j,i) - d_t_old(j,i)=d_t(j,i) + if(.not.(tnp .or. tnp1)) d_t_old(j,i)=d_t(j,i) d_a_old(j,i)=d_a(j,i) enddo enddo @@ -875,6 +1032,8 @@ c Compute long-range forces call etotal_long(energia_long) if (large.and. mod(itime,ntwe).eq.0) & call enerprint(energia_long) + E_long=energia_long(0) + potE=energia_short(0)+energia_long(0) #ifdef TIMING_ENE #ifdef MPI t_elong=t_elong+MPI_Wtime()-tt0 @@ -911,7 +1070,32 @@ c call cartprint endif c Compute the final RESPA step (increment velocities) c write (iout,*) "*********************** RESPA fin" - call RESPA_vel +C call RESPA_vel + if (tnp1) then + call tnp_respa_step2 + else if (tnp) then + call tnp_respa_step2 + else + call RESPA_vel + if (tnh.and..not.xiresp) then + call kinetic(EK) + call nhcint(EK,scale_nh,wdti,wdti2,wdti4,wdti8) + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t(j,i)*scale_nh + enddo + enddo + endif + endif + + if (tnp .or. tnp1) then + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t_old(j,i)/s_np + enddo + enddo + endif + c Compute the complete potential energy do i=0,n_ene potEcomp(i)=energia_short(i)+energia_long(i) @@ -937,6 +1121,73 @@ c Backup the coordinates, velocities, and accelerations & (d_t(j,i+nres),j=1,3) enddo endif + if (mod(itime,ntwe).eq.0) then + + if(tnp .or. tnp1) then +#ifndef G77 + write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit, + & E_long,energia_short(0) +#else + write (iout,'(a3,7f20.10)') "TTT",EK,s_np,potE,pi_np,Csplit, + & E_long,energia_short(0) +#endif + HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3) + H=(HNose1-H0)*s_np +cd write (iout,'(a,10f)') "hhh",EK,s_np,potE,pi_np,H0 +cd & ,EK+potE+pi_np**2/(2*Q_np)+dimen3*0.001986d0*t_bath*log(s_np) +cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0 + hhh=h +cd write (iout,'(a,3f)') "EE2 NP S, pi",totT, s_np, pi_np + endif + + if(tnh) then + HNose1=Hnose_nh(EK,potE) + H=HNose1-H0 +cd write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0 + hhh=h + endif + + + if (large) then + itnp=0 + do j=1,3 + itnp=itnp+1 + vtnp(itnp)=d_t(j,0) + enddo + do i=nnt,nct-1 + do j=1,3 + itnp=itnp+1 + vtnp(itnp)=d_t(j,i) + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + itnp=itnp+1 + vtnp(itnp)=d_t(j,inres) + enddo + endif + enddo +c Transform velocities from UNRES coordinate space to cartesian and Gvec +c eigenvector space + + do i=1,dimen3 + vtnp_(i)=0.0d0 + vtnp_a(i)=0.0d0 + do j=1,dimen3 + vtnp_(i)=vtnp_(i)+Gvec(j,i)*vtnp(j) + vtnp_a(i)=vtnp_a(i)+A(i,j)*vtnp(j) + enddo + vtnp_(i)=vtnp_(i)*dsqrt(geigen(i)) + enddo + + do i=1,dimen3 + write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i) + enddo + + endif + endif endif return end @@ -1663,6 +1914,20 @@ c Removing the velocity of the center of mass #endif #endif potE=potEcomp(0) + if(tnp .or. tnp1) then + s_np=1.0 + pi_np=0.0 + HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen3) + H0=Hnose1 + write(iout,*) 'H0= ',H0 + endif + + if(tnh) then + HNose1=Hnose_nh(EK,potE) + H0=HNose1 + write (iout,*) 'H0= ',H0 + endif + call cartgrad call lagrangian call max_accel @@ -1734,6 +1999,16 @@ c write (iout,*) "dc_old",i,(dc_old(j,i),j=1,3) t_eshort=t_eshort+tcpu()-tt0 #endif #endif + if(tnp .or. tnp1) then + E_short=energia_short(0) + HNose1=Hnose(EK,s_np,E_short,pi_np,Q_np,t_bath,dimen3) + Csplit=Hnose1 +c Csplit =110 +c_new_var_csplit Csplit=H0-E_long +c Csplit = H0-energia_short(0) + write(iout,*) 'Csplit= ',Csplit + endif + call cartgrad call lagrangian if(.not.out1file .and. large) then @@ -2492,3 +2767,667 @@ c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j) return end #endif + double precision function HNose(ek,s,e,pi,Q,t_bath,dimenl) + implicit none + double precision ek,s,e,pi,Q,t_bath,Rb + integer dimenl + Rb=0.001986d0 + HNose=ek+e+pi**2/(2*Q)+dimenl*Rb*t_bath*log(s) +c print '(6f15.5,i5,a2,2f15.5)',ek,s,e,pi,Q,t_bath,dimenl,"--", +c & pi**2/(2*Q),dimenl*Rb*t_bath*log(s) + return + end +c----------------------------------------------------------------- + double precision function HNose_nh(eki,e) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.MD' + HNose_nh=eki+e+dimen3*Rb*t_bath*xlogs(1)+qmass(1)*vlogs(1)**2/2 + do i=2,nnos + HNose_nh=HNose_nh+qmass(i)*vlogs(i)**2/2+Rb*t_bath*xlogs(i) + enddo +c write(4,'(5e15.5)') +c & vlogs(1),xlogs(1),HNose,eki,e + return + end +c----------------------------------------------------------------- + SUBROUTINE NHCINT(akin,scale,wdti,wdti2,wdti4,wdti8) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.MD' + double precision akin,gnkt,dt,aa,gkt,scale + double precision wdti(maxyosh),wdti2(maxyosh), + & wdti4(maxyosh),wdti8(maxyosh) + integer i,iresn,iyosh,inos,nnos1 + + dt=d_time + nnos1=nnos+1 + GKT = Rb*t_bath + GNKT = dimen3*GKT + akin=akin*2 + + +C THIS ROUTINE DOES THE NOSE-HOOVER PART OF THE +C INTEGRATION FROM t=0 TO t=DT/2 +C GET THE TOTAL KINETIC ENERGY + SCALE = 1.D0 +c CALL GETKINP(MASS,VX,VY,VZ,AKIN) +C UPDATE THE FORCES + GLOGS(1) = (AKIN - GNKT)/QMASS(1) +C START THE MULTIPLE TIME STEP PROCEDURE + DO IRESN = 1,NRESN + DO IYOSH = 1,NYOSH +C UPDATE THE THERMOSTAT VELOCITIES + VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH) + DO INOS = 1,NNOS-1 + AA = EXP(-WDTI8(IYOSH)*VLOGS(NNOS1-INOS) ) + VLOGS(NNOS-INOS) = VLOGS(NNOS-INOS)*AA*AA + & + WDTI4(IYOSH)*GLOGS(NNOS-INOS)*AA + ENDDO +C UPDATE THE PARTICLE VELOCITIES + AA = EXP(-WDTI2(IYOSH)*VLOGS(1) ) + SCALE = SCALE*AA +C UPDATE THE FORCES + GLOGS(1) = (SCALE*SCALE*AKIN - GNKT)/QMASS(1) +C UPDATE THE THERMOSTAT POSITIONS + DO INOS = 1,NNOS + XLOGS(INOS) = XLOGS(INOS) + VLOGS(INOS)*WDTI2(IYOSH) + ENDDO +C UPDATE THE THERMOSTAT VELOCITIES + DO INOS = 1,NNOS-1 + AA = EXP(-WDTI8(IYOSH)*VLOGS(INOS+1) ) + VLOGS(INOS) = VLOGS(INOS)*AA*AA + & + WDTI4(IYOSH)*GLOGS(INOS)*AA + GLOGS(INOS+1) = (QMASS(INOS)*VLOGS(INOS)*VLOGS(INOS) + & -GKT)/QMASS(INOS+1) + ENDDO + VLOGS(NNOS) = VLOGS(NNOS) + GLOGS(NNOS)*WDTI4(IYOSH) + ENDDO + ENDDO +C UPDATE THE PARTICLE VELOCITIES +c outside of this subroutine +c DO I = 1,N +c VX(I) = VX(I)*SCALE +c VY(I) = VY(I)*SCALE +c VZ(I) = VZ(I)*SCALE +c ENDDO + RETURN + END +c----------------------------------------------------------------- + subroutine tnp1_respa_i_step1 +c Applying Nose-Poincare algorithm - step 1 to coordinates +c JPSJ 70 75 (2001) S. Nose +c +c d_t is not updated here +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + double precision adt,adt2,tmp + + tmp=1+pi_np/(2*Q_np)*0.5*d_time + s12_np=s_np*tmp**2 + pistar=pi_np/tmp + s12_dt=d_time/s12_np + d_time_s12=d_time*0.5*s12_np + + do j=1,3 + d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12 + dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt + enddo + do i=nnt,nct-1 + do j=1,3 + d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12 + dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12 + dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt + enddo + endif + enddo + return + end +c--------------------------------------------------------------------- + subroutine tnp1_respa_i_step2 +c Step 2 of the velocity Verlet algorithm: update velocities + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + + double precision d_time_s12 + + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t_new(j,i) + enddo + enddo + + call kinetic(EK) + EK=EK/s12_np**2 + + d_time_s12=0.5d0*s12_np*d_time + + do j=1,3 + d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12 + enddo + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12 + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12 + enddo + endif + enddo + + pistar=pistar+(EK-0.5*(E_old+potE) + & -dimen3*Rb*t_bath*log(s12_np)+Csplit-dimen3*Rb*t_bath)*d_time + tmp=1+pistar/(2*Q_np)*0.5*d_time + s_np=s12_np*tmp**2 + pi_np=pistar/tmp + + return + end +c------------------------------------------------------- + + subroutine tnp1_step1 +c Applying Nose-Poincare algorithm - step 1 to coordinates +c JPSJ 70 75 (2001) S. Nose +c +c d_t is not updated here +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + double precision adt,adt2,tmp + + tmp=1+pi_np/(2*Q_np)*0.5*d_time + s12_np=s_np*tmp**2 + pistar=pi_np/tmp + s12_dt=d_time/s12_np + d_time_s12=d_time*0.5*s12_np + + do j=1,3 + d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12 + dc(j,0)=dc_old(j,0)+d_t_new(j,0)*s12_dt + enddo + do i=nnt,nct-1 + do j=1,3 + d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12 + dc(j,i)=dc_old(j,i)+d_t_new(j,i)*s12_dt + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s12 + dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*s12_dt + enddo + endif + enddo + return + end +c--------------------------------------------------------------------- + subroutine tnp1_step2 +c Step 2 of the velocity Verlet algorithm: update velocities + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + + double precision d_time_s12 + + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t_new(j,i) + enddo + enddo + + call kinetic(EK) + EK=EK/s12_np**2 + + d_time_s12=0.5d0*s12_np*d_time + + do j=1,3 + d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s12 + enddo + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s12 + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s12 + enddo + endif + enddo + +cd write(iout,*) 'pistar',pistar,EK,E_old,potE,s12_np + pistar=pistar+(EK-0.5*(E_old+potE) + & -dimen3*Rb*t_bath*log(s12_np)+H0-dimen3*Rb*t_bath)*d_time + tmp=1+pistar/(2*Q_np)*0.5*d_time + s_np=s12_np*tmp**2 + pi_np=pistar/tmp + + return + end + +c----------------------------------------------------------------- + subroutine tnp_respa_i_step1 +c Applying Nose-Poincare algorithm - step 1 to coordinates +c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird +c +c d_t is not updated here, it is destroyed +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + double precision C_np,d_time_s,tmp,d_time_ss + + d_time_s=d_time*0.5*s_np +ct2 d_time_s=d_time*0.5*s12_np + + do j=1,3 + d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s + enddo + do i=nnt,nct-1 + do j=1,3 + d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s + enddo + endif + enddo + + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t_new(j,i) + enddo + enddo + + call kinetic(EK) + EK=EK/s_np**2 + + C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-Csplit) + & -pi_np + + pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np)) + tmp=0.5*d_time*pistar/Q_np + s12_np=s_np*(1.0+tmp)/(1.0-tmp) + + d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np) +ct2 d_time_ss=d_time/s12_np +c d_time_ss=0.5*d_time*(1.0/sold_np+1.0/s_np) + + do j=1,3 + dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss + enddo + do i=nnt,nct-1 + do j=1,3 + dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss + enddo + endif + enddo + + return + end +c--------------------------------------------------------------------- + + subroutine tnp_respa_i_step2 +c Step 2 of the velocity Verlet algorithm: update velocities + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + + double precision d_time_s + + EK=EK*(s_np/s12_np)**2 + HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3) + pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath + & -HNose1+Csplit) + +cr print '(a,5f)','i_step2',EK,potE,HNose1,pi_np,E_long + d_time_s=d_time*0.5*s12_np +c d_time_s=d_time*0.5*s_np + + do j=1,3 + d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s + enddo + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s + enddo + endif + enddo + + s_np=s12_np + + return + end +c----------------------------------------------------------------- + subroutine tnp_respa_step1 +c Applying Nose-Poincare algorithm - step 1 to vel for RESPA +c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird +c +c d_t is not updated here, it is destroyed +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + double precision C_np,d_time_s,tmp,d_time_ss + double precision energia(0:n_ene) + + d_time_s=d_time*0.5*s_np + + do j=1,3 + d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s + enddo + do i=nnt,nct-1 + do j=1,3 + d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s + enddo + endif + enddo + + +c C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0) +c & -pi_np +c +c pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np)) +c tmp=0.5*d_time*pistar/Q_np +c s12_np=s_np*(1.0+tmp)/(1.0-tmp) +c write(iout,*) 'tnp_respa_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp + +ct1 pi_np=pistar +c sold_np=s_np +c s_np=s12_np + +c------------------------------------- +c test of reviewer's comment + pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0) +cr print '(a,3f)','1 pi_np,s_np',pi_np,s_np,E_long +c------------------------------------- + + return + end +c--------------------------------------------------------------------- + subroutine tnp_respa_step2 +c Step 2 of the velocity Verlet algorithm: update velocities for RESPA + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + + double precision d_time_s + +ct1 s12_np=s_np +ct2 pistar=pi_np + +ct call kinetic(EK) +ct HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3) +ct pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath) +ct & -0.5*d_time*(HNose1-H0) + +c------------------------------------- +c test of reviewer's comment + pi_np=pi_np-0.5*d_time*(E_long+Csplit-H0) +cr print '(a,3f)','2 pi_np,s_np',pi_np,s_np,E_long +c------------------------------------- + d_time_s=d_time*0.5*s_np + + do j=1,3 + d_t_old(j,0)=d_t_old(j,0)+d_a(j,0)*d_time_s + enddo + do i=nnt,nct-1 + do j=1,3 + d_t_old(j,i)=d_t_old(j,i)+d_a(j,i)*d_time_s + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t_old(j,inres)=d_t_old(j,inres)+d_a(j,inres)*d_time_s + enddo + endif + enddo + +cd s_np=s12_np + + return + end +c--------------------------------------------------------------------- + subroutine tnp_step1 +c Applying Nose-Poincare algorithm - step 1 to coordinates +c J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird +c +c d_t is not updated here, it is destroyed +c + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + double precision C_np,d_time_s,tmp,d_time_ss + + d_time_s=d_time*0.5*s_np + + do j=1,3 + d_t_new(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s + enddo + do i=nnt,nct-1 + do j=1,3 + d_t_new(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t_new(j,inres)=d_t_old(j,inres)+d_a_old(j,inres)*d_time_s + enddo + endif + enddo + + do i=0,2*nres + do j=1,3 + d_t(j,i)=d_t_new(j,i) + enddo + enddo + + call kinetic(EK) + EK=EK/s_np**2 + + C_np=0.5*d_time*(dimen3*Rb*t_bath*(1.0+log(s_np))-EK+potE-H0) + & -pi_np + + pistar=-2.0*C_np/(1.0+sqrt(1.0-C_np*d_time/Q_np)) + tmp=0.5*d_time*pistar/Q_np + s12_np=s_np*(1.0+tmp)/(1.0-tmp) +c write(iout,*) 'tnp_step1',s_np,s12_np,EK,potE,C_np,pistar,tmp + + d_time_ss=0.5*d_time*(1.0/s12_np+1.0/s_np) + + do j=1,3 + dc(j,0)=dc_old(j,0)+d_t_new(j,0)*d_time_ss + enddo + do i=nnt,nct-1 + do j=1,3 + dc(j,i)=dc_old(j,i)+d_t_new(j,i)*d_time_ss + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + dc(j,inres)=dc_old(j,inres)+d_t_new(j,inres)*d_time_ss + enddo + endif + enddo + + return + end +c----------------------------------------------------------------- + subroutine tnp_step2 +c Step 2 of the velocity Verlet algorithm: update velocities + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.VAR' + include 'COMMON.MD' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.GEO' + include 'COMMON.LOCAL' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.NAMES' + + double precision d_time_s + + EK=EK*(s_np/s12_np)**2 + HNose1=Hnose(EK,s12_np,potE,pistar,Q_np,t_bath,dimen3) + pi_np=pistar+0.5*d_time*(2*EK-dimen3*Rb*t_bath) + & -0.5*d_time*(HNose1-H0) + +cd write(iout,'(a,4f)') 'mmm',EK,potE,HNose1,pi_np + d_time_s=d_time*0.5*s12_np + + do j=1,3 + d_t(j,0)=d_t_new(j,0)+d_a(j,0)*d_time_s + enddo + do i=nnt,nct-1 + do j=1,3 + d_t(j,i)=d_t_new(j,i)+d_a(j,i)*d_time_s + enddo + enddo + do i=nnt,nct + if (itype(i).ne.10) then + inres=i+nres + do j=1,3 + d_t(j,inres)=d_t_new(j,inres)+d_a(j,inres)*d_time_s + enddo + endif + enddo + + s_np=s12_np + + return + end + 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 f7ff34c..51867aa 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -12106,10 +12106,27 @@ C for UNRES C lets ommit dummy atoms for now if ((itype(i).eq.ntyp1).or.(itype(i+1).eq.ntyp1)) cycle C now calculate distance from center of tube and direction vectors - vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) - if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize - vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) - if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp vectube(1)=vectube(1)-tubecenter(1) vectube(2)=vectube(2)-tubecenter(2) @@ -12129,12 +12146,12 @@ C calculte rdiffrence between r and r0 C and its 6 power rdiff6=rdiff**6.0d0 C for vectorization reasons we will sumup at the end to avoid depenence of previous - enetube(i)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6 + enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6 C write(iout,*) "TU13",i,rdiff6,enetube(i) C print *,rdiff,rdiff6,pep_aa_tube C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 C now we calculate gradient - fac=(-12.0d0*pep_aa_tube/rdiff6+ + fac=(-12.0d0*pep_aa_tube/rdiff6- & 6.0d0*pep_bb_tube)/rdiff6/rdiff C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), C &rdiff,fac @@ -12154,13 +12171,29 @@ C lets ommit dummy atoms for now C in UNRES uncomment the line below as GLY has no side-chain... C .or.(iti.eq.10) & ) cycle - vectube(1)=c(1,i+nres) - vectube(1)=mod(vectube(1),boxxsize) - if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize - vectube(2)=c(2,i+nres) - vectube(2)=mod(vectube(2),boxysize) - if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize - + xmin=boxxsize + ymin=boxysize + do j=-1,1 + vectube(1)=mod((c(1,i+nres)),boxxsize) + vectube(1)=vectube(1)+boxxsize*j + vectube(2)=mod((c(2,i+nres)),boxysize) + vectube(2)=vectube(2)+boxysize*j + + xminact=abs(vectube(1)-tubecenter(1)) + yminact=abs(vectube(2)-tubecenter(2)) + if (xmin.gt.xminact) then + xmin=xminact + xtemp=vectube(1) + endif + if (ymin.gt.yminact) then + ymin=yminact + ytemp=vectube(2) + endif + enddo + vectube(1)=xtemp + vectube(2)=ytemp +C write(iout,*), "tututu", vectube(1),tubecenter(1),vectube(2), +C & tubecenter(2) vectube(1)=vectube(1)-tubecenter(1) vectube(2)=vectube(2)-tubecenter(2) @@ -12172,6 +12205,7 @@ C now calculte the distance C now normalize vector vectube(1)=vectube(1)/tub_r vectube(2)=vectube(2)/tub_r + C calculte rdiffrence between r and r0 rdiff=tub_r-tubeR0 C and its 6 power @@ -12179,10 +12213,10 @@ C and its 6 power C for vectorization reasons we will sumup at the end to avoid depenence of previous sc_aa_tube=sc_aa_tube_par(iti) sc_bb_tube=sc_bb_tube_par(iti) - enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0-sc_bb_tube/rdiff6 + enetube(i+nres)=sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6 C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 C now we calculate gradient - fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff+ + fac=-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- & 6.0d0*sc_bb_tube/rdiff6/rdiff C now direction of gg_tube vector do j=1,3 @@ -12265,12 +12299,12 @@ C calculte rdiffrence between r and r0 C and its 6 power rdiff6=rdiff**6.0d0 C for vectorization reasons we will sumup at the end to avoid depenence of previous - enetube(i)=pep_aa_tube/rdiff6**2.0d0-pep_bb_tube/rdiff6 + enetube(i)=pep_aa_tube/rdiff6**2.0d0+pep_bb_tube/rdiff6 C write(iout,*) "TU13",i,rdiff6,enetube(i) C print *,rdiff,rdiff6,pep_aa_tube C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 C now we calculate gradient - fac=(-12.0d0*pep_aa_tube/rdiff6+ + fac=(-12.0d0*pep_aa_tube/rdiff6- & 6.0d0*pep_bb_tube)/rdiff6/rdiff C write(iout,'(a5,i4,f12.1,3f12.5)') "TU13",i,rdiff6,enetube(i), C &rdiff,fac @@ -12360,11 +12394,11 @@ C and its 6 power C for vectorization reasons we will sumup at the end to avoid depenence of previous sc_aa_tube=sc_aa_tube_par(iti) sc_bb_tube=sc_bb_tube_par(iti) - enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0-sc_bb_tube/rdiff6) + enetube(i+nres)=(sc_aa_tube/rdiff6**2.0d0+sc_bb_tube/rdiff6) & *sstube+enetube(i+nres) C pep_aa_tube and pep_bb_tube are precomputed values A=4eps*sigma^12 B=4eps*sigma^6 C now we calculate gradient - fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff+ + fac=(-12.0d0*sc_aa_tube/rdiff6**2.0d0/rdiff- & 6.0d0*sc_bb_tube/rdiff6/rdiff)*sstube C now direction of gg_tube vector do j=1,3 diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index fa4c75a..577c1ee 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -424,6 +424,14 @@ c call reada(controlcard,"R_CUT",r_cut,2.0d0) c call reada(controlcard,"LAMBDA",rlamb,0.3d0) rest = index(controlcard,"REST").gt.0 tbf = index(controlcard,"TBF").gt.0 + tnp = index(controlcard,"NOSEPOINCARE99").gt.0 + tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0 + tnh = index(controlcard,"NOSEHOOVER96").gt.0 + if (RESPA.and.tnh)then + xiresp = index(controlcard,"XIRESP").gt.0 + endif + write(iout,*) "tnh",tnh + call reada(controlcard,"Q_NP",Q_np,0.1d0) usampl = index(controlcard,"USAMPL").gt.0 mdpdb = index(controlcard,"MDPDB").gt.0 @@ -546,6 +554,43 @@ c Calculate friction coefficients and bounds of stochastic forces & "Velocities will be reset at random every",count_reset_vel, & " steps" endif + else if (tnp .or. tnp1 .or. tnh) then + if (tnp .or. tnp1) then + write (iout,'(a)') "Nose-Poincare bath calculation" + if (tnp) write (iout,'(a)') + & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird" + if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose" + else + write (iout,'(a)') "Nose-Hoover bath calculation" + write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al." + nresn=1 + nyosh=1 + nnos=1 + do i=1,nnos + qmass(i)=Q_np + xlogs(i)=1.0 + vlogs(i)=0.0 + enddo + do i=1,nyosh + WDTI(i) = 1.0*d_time/nresn + WDTI2(i)=WDTI(i)/2 + WDTI4(i)=WDTI(i)/4 + WDTI8(i)=WDTI(i)/8 + enddo + if (RESPA) then + if(xiresp) then + write (iout,'(a)') "NVT-XI-RESPA algorithm" + else + write (iout,'(a)') "NVT-XO-RESPA algorithm" + endif + do i=1,nyosh + WDTIi(i) = 1.0*d_time/nresn/ntime_split + WDTIi2(i)=WDTIi(i)/2 + WDTIi4(i)=WDTIi(i)/4 + WDTIi8(i)=WDTIi(i)/8 + enddo + endif + endif else if(me.eq.king.or..not.out1file) & write (iout,'(a31)') "Microcanonical mode calculation" diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index a57f604..b5f5705 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -250,11 +250,28 @@ C & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(2)*gsccorx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + else gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i) & +fact(1)*wscp*gvdwc_scp(j,i)+ @@ -312,11 +329,28 @@ C & wturn6*fact(5)*gcorr6_turn(j,i)+ & wsccor*fact(2)*gsccorc(j,i) & +wliptran*gliptranc(j,i) + & +welec*gshieldc(j,i) + & +welec*gshieldc_loc(j,i) + & +wcorr*gshieldc_ec(j,i) + & +wcorr*gshieldc_loc_ec(j,i) + & +wturn3*gshieldc_t3(j,i) + & +wturn3*gshieldc_loc_t3(j,i) + & +wturn4*gshieldc_t4(j,i) + & +wturn4*gshieldc_loc_t4(j,i) + & +wel_loc*gshieldc_ll(j,i) + & +wel_loc*gshieldc_loc_ll(j,i) + gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*fact(1)*gsccorx(j,i) & +wliptran*gliptranx(j,i) + & +welec*gshieldx(j,i) + & +wcorr*gshieldx_ec(j,i) + & +wturn3*gshieldx_t3(j,i) + & +wturn4*gshieldx_t4(j,i) + & +wel_loc*gshieldx_ll(j,i) + else gradc(j,i,icg)=fact(1)*wsc*gvdwc(j,i)+ & fact(1)*wscp*gvdwc_scp(j,i)+ @@ -2171,6 +2205,31 @@ C endif if (ymedi.lt.0) ymedi=ymedi+boxysize zmedi=mod(zmedi,boxzsize) if (zmedi.lt.0) zmedi=zmedi+boxzsize + zmedi2=mod(zmedi,boxzsize) + if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize + if ((zmedi2.gt.bordlipbot) + &.and.(zmedi2.lt.bordliptop)) then +C the energy transfer exist + if (zmedi2.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zmedi2-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipi=sscalelip(fracinbuf) + ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick + elseif (zmedi2.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zmedi2)/lipbufthick) + sslipi=sscalelip(fracinbuf) + ssgradlipi=sscagradlip(fracinbuf)/lipbufthick + else + sslipi=1.0d0 + ssgradlipi=0.0d0 + endif + else + sslipi=0.0d0 + ssgradlipi=0.0d0 + endif + num_conti=0 C write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) do j=ielstart(i),ielend(i) @@ -2216,6 +2275,28 @@ C End diagnostics if (yj.lt.0) yj=yj+boxysize zj=mod(zj,boxzsize) if (zj.lt.0) zj=zj+boxzsize + if ((zj.gt.bordlipbot) + &.and.(zj.lt.bordliptop)) then +C the energy transfer exist + if (zj.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zj-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2 xj_safe=xj yj_safe=yj @@ -2722,6 +2803,7 @@ C fac_shield(j)=0.6 endif eel_loc_ij=eel_loc_ij & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij C write (iout,'(a6,2i5,0pf7.3)') C & 'eelloc',i,j,eel_loc_ij @@ -2779,11 +2861,13 @@ C Partial derivatives in virtual-bond dihedral angles gamma & (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j) & +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc_loc(j-1)=gel_loc_loc(j-1)+ & (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j) & +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) cd call checkint3(i,j,mu1,mu2,a22,a23,a32,a33,acipa,eel_loc_ij) cd write(iout,*) 'agg ',agg @@ -2797,6 +2881,7 @@ C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2) ggg(l)=(agg(l,1)*muij(1)+ & agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo do k=i+2,j2 @@ -2809,18 +2894,22 @@ C Remaining derivatives of eello gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+ & aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+ & aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+ & aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+ & aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo endif @@ -3083,6 +3172,37 @@ C Third- and fourth-order contributions from turns double precision agg(3,4),aggi(3,4),aggi1(3,4), & aggj(3,4),aggj1(3,4),a_temp(2,2) common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2 + zj=(c(3,j)+c(3,j+1))/2.0d0 +C xj=mod(xj,boxxsize) +C if (xj.lt.0) xj=xj+boxxsize +C yj=mod(yj,boxysize) +C if (yj.lt.0) yj=yj+boxysize + zj=mod(zj,boxzsize) + if (zj.lt.0) zj=zj+boxzsize +C if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ" + if ((zj.gt.bordlipbot) + &.and.(zj.lt.bordliptop)) then +C the energy transfer exist + if (zj.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((zj-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick) + sslipj=sscalelip(fracinbuf) + ssgradlipj=sscagradlip(fracinbuf)/lipbufthick + else + sslipj=1.0d0 + ssgradlipj=0.0 + endif + else + sslipj=0.0d0 + ssgradlipj=0.0 + endif + if (j.eq.i+2) then if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds @@ -3120,8 +3240,11 @@ C fac_shield(j)=0.6 eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + eello_t3=0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) cd write (2,*) 'i,',i,' j',j,'eello_turn3', cd & 0.5d0*(pizda(1,1)+pizda(2,2)), @@ -3176,6 +3299,8 @@ C Derivatives in gamma(i) call matmat2(a_temp(1,1),pizda(1,1),pizda(1,1)) gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) + C Derivatives in gamma(i+1) call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1)) call transpose2(auxmat2(1,1),pizda(1,1)) @@ -3183,6 +3308,7 @@ C Derivatives in gamma(i+1) gel_loc_turn3(i+1)=gel_loc_turn3(i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Cartesian derivatives do l=1,3 @@ -3194,6 +3320,7 @@ C Cartesian derivatives gcorr3_turn(l,i)=gcorr3_turn(l,i) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) @@ -3203,6 +3330,7 @@ C Cartesian derivatives gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) @@ -3212,6 +3340,7 @@ C Cartesian derivatives gcorr3_turn(l,j)=gcorr3_turn(l,j) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) @@ -3221,6 +3350,7 @@ C Cartesian derivatives gcorr3_turn(l,j1)=gcorr3_turn(l,j1) & +0.5d0*(pizda(1,1)+pizda(2,2)) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo endif @@ -3330,6 +3460,7 @@ C & *2.0 s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Derivatives in gamma(i+1) call transpose2(EUgder(1,1,i+2),e2tder(1,1)) @@ -3340,6 +3471,7 @@ C Derivatives in gamma(i+1) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Derivatives in gamma(i+2) call transpose2(EUgder(1,1,i+3),e3tder(1,1)) @@ -3353,6 +3485,7 @@ C Derivatives in gamma(i+2) s3=0.5d0*(pizda(1,1)+pizda(2,2)) gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) C Cartesian derivatives @@ -3375,6 +3508,7 @@ C Derivatives of this turn contributions in DC(i+2) ggg(l)=-(s1+s2+s3) gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo endif @@ -3395,6 +3529,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggi1(l,1) a_temp(1,2)=aggi1(l,2) @@ -3411,6 +3546,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj(l,1) a_temp(1,2)=aggj(l,2) @@ -3427,6 +3563,7 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) a_temp(1,1)=aggj1(l,1) a_temp(1,2)=aggj1(l,2) @@ -3443,8 +3580,17 @@ C Remaining derivatives of this turn contribution s3=0.5d0*(pizda(1,1)+pizda(2,2)) gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3) & *fac_shield(i)*fac_shield(j) + &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0) enddo + gshieldc_t4(3,i)=gshieldc_t4(3,i)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j)=gshieldc_t4(3,j)+ + & ssgradlipj*eello_t4/4.0d0*lipscale + gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+ + & ssgradlipi*eello_t4/4.0d0*lipscale + gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+ + & ssgradlipj*eello_t4/4.0d0*lipscale endif 178 continue endif -- 1.7.9.5