X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2FMD_A-MTS.F;h=6650c5a9e61c84013db53ec920023c7597ebc4d8;hb=bde71ffff7527add147c1bee4bad23f02b4d1ab7;hp=56b3ea89b487f9e074a1044bab92074b2bbf405c;hpb=af72f8e89a5d33f0d86ba898d6c5bbbda4b25b84;p=unres.git diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 56b3ea8..6650c5a 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) +C call enerprint(potEcomp) +C print *,itime,'AFM',Eafmforc,etot + endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) @@ -209,7 +213,7 @@ c Variable time step algorithm. enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 ind=ind+1 v_work(ind)=d_t(j,i+nres) @@ -230,6 +234,9 @@ c Variable time step algorithm. #endif endif if (mod(itime,ntwx).eq.0) then + write(iout,*) 'time=',itime +C call check_ecartint + call returnbox write (tytul,'("time",f8.2)') totT if(mdpdb) then call hairpin(.true.,nharp,iharp) @@ -316,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 @@ -359,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 @@ -404,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 @@ -508,11 +530,26 @@ 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 + totTafm=totT +C print *,totTafm,"TU?" if (d_time.ne.d_time0) then d_time=d_time0 #ifndef LANG0 @@ -537,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 @@ -553,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), @@ -566,6 +667,7 @@ c Backup the coordinates, velocities, and accelerations enddo endif endif + endif return end c------------------------------------------------------------------------------- @@ -610,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 @@ -635,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" @@ -652,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 @@ -680,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 @@ -732,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 @@ -762,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 @@ -838,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 @@ -866,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 @@ -902,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) @@ -910,6 +1103,7 @@ c Compute the complete potential energy potE=potEcomp(0)-potEcomp(20) c potE=energia_short(0)+energia_long(0) totT=totT+d_time + totTafm=totT c Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE @@ -927,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 @@ -955,7 +1216,7 @@ c forces). enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t(j,inres)+0.5d0*d_a(j,inres)*d_time @@ -996,6 +1257,8 @@ c Applying velocity Verlet algorithm - step 1 to coordinates d_t(j,0)=d_t_old(j,0)+adt enddo do i=nnt,nct-1 +C SPYTAC ADAMA +C do i=0,nres do j=1,3 adt=d_a_old(j,i)*d_time adt2=0.5d0*adt @@ -1005,7 +1268,8 @@ c Applying velocity Verlet algorithm - step 1 to coordinates enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then +C do i=0,nres + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 adt=d_a_old(j,inres)*d_time @@ -1049,7 +1313,7 @@ c Step 2 of the velocity Verlet algorithm: update velocities enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+0.5d0*d_a(j,inres)*d_time @@ -1161,7 +1425,7 @@ c ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 adt=(d_a_old(j,inres)+d_af_work(ind+j))*d_time @@ -1226,7 +1490,7 @@ c ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_new(j,inres)+(0.5d0*(d_a(j,inres) @@ -1273,6 +1537,7 @@ c accel(j)=aux(j)+0.5d0*(d_a(j,i)-d_a_old(j,i)) c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) if (dabs(accel(j)).gt.dabs(accel_old(j))) then dacc=dabs(accel(j)-accel_old(j)) +c write (iout,*) i,dacc if (dacc.gt.amax) amax=dacc endif enddo @@ -1291,7 +1556,7 @@ c accel(j)=aux(j) enddo endif do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) accel_old(j)=accel_old(j)+d_a_old(j,i+nres) @@ -1302,6 +1567,7 @@ c accel(j)=accel(j)+d_a(j,i+nres)-d_a_old(j,i+nres) c if (dabs(accel(j)).gt.amax) amax=dabs(accel(j)) if (dabs(accel(j)).gt.dabs(accel_old(j))) then dacc=dabs(accel(j)-accel_old(j)) +c write (iout,*) "side-chain",i,dacc if (dacc.gt.amax) amax=dacc endif enddo @@ -1344,7 +1610,7 @@ c write (iout,*) "back",i,j,epdriftij enddo endif c Side chains - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 epdriftij= & dabs((d_a(j,i+nres)-d_a_old(j,i+nres))*gxcart(j,i)) @@ -1391,7 +1657,7 @@ c write(iout,*) "fact", fact enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=fact*d_t(j,inres) @@ -1446,7 +1712,8 @@ c if the friction coefficients do not depend on surface area stdforcp(i)=stdfp*dsqrt(gamp) enddo do i=nnt,nct - stdforcsc(i)=stdfsc(itype(i))*dsqrt(gamsc(itype(i))) + stdforcsc(i)=stdfsc(iabs(itype(i))) + & *dsqrt(gamsc(iabs(itype(i)))) enddo endif c Open the pdb file for snapshotshots @@ -1556,6 +1823,7 @@ c inquire(file=mremd_rst_name,exist=file_exist) endif call random_vel totT=0.0d0 + totTafm=totT endif else c Generate initial velocities @@ -1563,6 +1831,8 @@ c Generate initial velocities & write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 +CtotTafm is the variable for AFM time which eclipsed during + totTafm=totT endif c rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then @@ -1634,6 +1904,7 @@ c Removing the velocity of the center of mass #endif call zerograd call etotal(potEcomp) + call enerprint(potEcomp) if (large) call enerprint(potEcomp) #ifdef TIMING_ENE #ifdef MPI @@ -1643,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 @@ -1652,11 +1937,11 @@ c Removing the velocity of the center of mass & "Time step reduced to",d_time, & " because of too large initial acceleration." endif - if(me.eq.king.or..not.out1file)then - write(iout,*) "Potential energy and its components" - call enerprint(potEcomp) +C if(me.eq.king.or..not.out1file)then +C write(iout,*) "Potential energy and its components" +C call enerprint(potEcomp) c write(iout,*) (potEcomp(i),i=0,n_ene) - endif +C endif potE=potEcomp(0)-potEcomp(20) totE=EK+potE itime=0 @@ -1714,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 @@ -1785,7 +2080,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 @@ -1799,10 +2094,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 +C if (SELFGUIDE.gt.0) then +C distance=0.0 +C do j=1,3 +C vec_afm(j)=c(j,afmend)-c(j,afmbeg) +C distance=distance+vec_afm(j)**2 +C enddo +C distance=dsqrt(distance) +C do j=1,3 +C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance +C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance +C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3), +C & d_t_work_new(j+(afmend-1)*3) +C enddo + +C endif + c diagnostics c Ek1=0.0d0 c ii=0 @@ -1839,7 +2151,7 @@ c Transfer to the d_t vector enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 ind=ind+1 d_t(j,i+nres)=d_t_work(ind) @@ -2068,7 +2380,7 @@ c enddo ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) @@ -2177,7 +2489,7 @@ c ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) @@ -2334,7 +2646,7 @@ c enddo ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 dc_work(ind+j)=dc_old(j,i+nres) d_t_work(ind+j)=d_t_old(j,i+nres) @@ -2383,7 +2695,7 @@ c enddo ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 dc(j,inres)=dc_work(ind+j) @@ -2444,7 +2756,7 @@ c ddt2=ddt2+vrand_mat2(i,j)*stochforcvecV(j) ind=ind+3 enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 d_t(j,inres)=d_t_work(ind+j) @@ -2455,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 +