+++ /dev/null
-c-------------------------------------------------------------------------------
- subroutine RESPA_step(itime)
-c-------------------------------------------------------------------------------
-c Perform a single RESPA step.
-c-------------------------------------------------------------------------------
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
- include 'COMMON.LANGEVIN'
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- include 'COMMON.TIME1'
- double precision energia(0:n_ene),energia_short(0:n_ene),
- & energia_long(0:n_ene)
- double precision cm(3),L(3),vcm(3),incr(3)
- integer ilen,count,rstcount
- external ilen
- character*50 tytul
- integer maxcount_scale /10/
- common /gucio/ cm
- double precision stochforcvec(MAXRES6)
- double precision grad_tmp(3,0:maxres2)
- common /stochcalc/ stochforcvec
- integer itime
- logical scale
- double precision vtnp(maxres6), vtnp_(maxres6), vtnp_a(maxres6)
- if (large.and. mod(itime,ntwe).eq.0) then
- write (iout,*) "***************** RESPA itime",itime
- write (iout,*) "Cartesian and internal coordinates: step 0"
- call cartprint
- call intout
- write (iout,*) "Accelerations"
- do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
- & (d_a(j,i+nres),j=1,3)
- enddo
- write (iout,*) "Velocities, step 0"
- do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
- & (d_t(j,i+nres),j=1,3)
- enddo
- endif
-c
-c Perform the initial RESPA step (increment velocities)
-c write (iout,*) "*********************** RESPA ini"
-
- if (tnp1) then
-creview call tnp1_respa_step1
- 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 (mod(itime,ntwe).eq.0 .and. large) then
- write (iout,*) "Velocities, end"
- do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
- & (d_t(j,i+nres),j=1,3)
- enddo
- endif
-c
-cr if (tnp) then
-cr do i=0,nres
-cr do j=1,3
-cr grad_tmp(j,i)=gcart(j,i)
-cr grad_tmp(j,i+nres)=gxcart(j,i)
-cr enddo
-cr enddo
-cr endif
-c
-c Compute the short-range forces
- call zerograd
- call etotal_short(energia_short)
- if (tnp.or.tnp1) potE=energia_short(0)
- call cartgrad
- call lagrangian
- do i=0,2*nres
- do j=1,3
- dc_old(j,i)=dc(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
- d_time0=d_time
-c Split the time step
- d_time=d_time/ntime_split
-ctest E_long2=E_long
-c Perform the short-range RESPSA steps (velocity Verlet increments of
-c positions and velocities using short-range forces)
-c write (iout,*) "*********************** RESPA split"
-creview E_old=potE
- do itsplit=1,ntime_split
- if (lang.eq.1) then
- call sddir_precalc
- else if (lang.eq.2 .or. lang.eq.3) then
- call stochastic_force(stochforcvec)
- endif
-c First step of the velocity Verlet algorithm
- if (lang.eq.2) then
- call sd_verlet1
- else if (lang.eq.3) then
- call sd_verlet1_ciccotti
- else if (lang.eq.1) then
- call sddir_verlet1
- else if (tnp1) then
- call tnp1_respa_i_step1
-cr if(itsplit.eq.1)then
-cr d_time_s12=d_time0*0.5*s12_np
-cr do j=1,3
-cr d_t_half(j,0)=d_t_old(j,0)+d_a_old(j,0)*d_time_s12
-cr enddo
-cr do i=nnt,nct-1
-cr do j=1,3
-cr d_t_half(j,i)=d_t_old(j,i)+d_a_old(j,i)*d_time_s12
-cr enddo
-cr enddo
-cr do i=nnt,nct
-cr if (itype(i).ne.10) then
-cr inres=i+nres
-cr do j=1,3
-cr d_t_half(j,inres)=d_t_old(j,inres)
-cr & +d_a_old(j,inres)*d_time_s12
-cr enddo
-cr endif
-cr enddo
-cr endif
- 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
- call chainbuild_cart
- if (rattle) call rattle1
- if (large.and. mod(itime,ntwe).eq.0) then
- write (iout,*) "Cartesian and internal coordinates: step 1"
- call cartprint
- call intout
- write (iout,*) "Accelerations"
- do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
- & (d_a(j,i+nres),j=1,3)
- enddo
- write (iout,*) "Velocities, step 1"
- do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
- & (d_t(j,i+nres),j=1,3)
- enddo
- endif
- tt0 = tcpu()
-c Calculate energy and forces
-c
-c E_long aproximation
-cr if (tnp .or. tnp1) then
-cr dtmp=0.5*d_time*(1.0/s12_np+1.0/s_np)
-cr do i=0,2*nres
-cr do j=1,3
-cr E_long=E_long+d_t_new(j,i)*dtmp*grad_tmp(j,i)
-cr enddo
-cr enddo
-cr endif
-c-------------------------------------
-c test of reviewer's comment
-cr E_long=0
-c-------------------------------------
-
-
-c
-ctest call etotal_long(energia_long)
-ctest E_long=energia_long(0)
-ctest
- call zerograd
-
- call etotal_short(energia_short)
- E_old=potE
- potE=energia_short(0)
-
-c if(tnp .or. tnp1) then
-c write (iout,*) "kkk",E_long2,E_long
-c endif
-
-
- call cartgrad
-c Get the new accelerations
- call lagrangian
- t_enegrad=t_enegrad+tcpu()-tt0
-c Second step of the velocity Verlet algorithm
- if (lang.eq.2) then
- call sd_verlet2
- else if (lang.eq.3) then
- call sd_verlet2_ciccotti
- 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
-c Backup the coordinates, velocities, and accelerations
- 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
- do i=0,2*nres
- do j=1,3
- dc_old(j,i)=dc(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
-
-cd if(tnp .or. tnp1) then
-cd call kinetic(EK)
-cd HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
-cd H=(HNose1-H0)*s_np
-cd write (iout,*) "jjj",EK,potE
-cd write (iout,*) "iii H=",H,abs(HNose1-H0)/H0
-cd write (iout,'(a,3f)')
-cd & "III NP S, pi",totT+itsplit*d_time, s_np, pi_np
-cd endif
-
-
- enddo
-c Restore the time step
- d_time=d_time0
-c Compute long-range forces
- call zerograd
- call etotal_long(energia_long)
- E_long=energia_long(0)
- potE=energia_short(0)+energia_long(0)
- call cartgrad
-c Compute accelerations from long-range forces
- call lagrangian
- if (large.and. mod(itime,ntwe).eq.0) then
- write (iout,*) "Cartesian and internal coordinates: step 2"
- call cartprint
- call intout
- write (iout,*) "Accelerations"
- do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_a(j,i),j=1,3),
- & (d_a(j,i+nres),j=1,3)
- 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),
- & (d_t(j,i+nres),j=1,3)
- enddo
- endif
-c Compute the final RESPA step (increment velocities)
-c write (iout,*) "*********************** RESPA fin"
- if (tnp1) then
-creview call tnp1_respa_step2
- 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
-cc potE=energia_short(0)+energia_long(0)
- totT=totT+d_time
-c Calculate the kinetic and the total energy and the kinetic temperature
- call kinetic(EK)
- totE=EK+potE
-
-c Couple the system to Berendsen bath if needed
- if (tbf .and. lang.eq.0) then
- call verlet_bath
- endif
- kinetic_T=2.0d0/(dimen*Rb)*EK
-c Backup the coordinates, velocities, and accelerations
- if (mod(itime,ntwe).eq.0 .and. large) then
- write (iout,*) "Velocities, end"
- do i=0,nres
- write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
- & (d_t(j,i+nres),j=1,3)
- enddo
- endif
-
-c-----review
-c if(tnp .or. tnp1) then
-c HNose1=Hnose(EK,s_np,energia_short(0),pi_np,Q_np,t_bath,dimen)
-c_new_var_csplit Csplit=H0-E_long
-c Csplit=H0-energia_short(0)
-c endif
-c----------
-
-
- if (mod(itime,ntwe).eq.0) then
-
- if(tnp .or. tnp1) then
- write (iout,'(a3,7f)') "TTT",EK,s_np,potE,pi_np,Csplit,
- & E_long,energia_short(0)
- HNose1=Hnose(EK,s_np,potE,pi_np,Q_np,t_bath,dimen)
- 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)+dimen*0.001986d0*t_bath*log(s_np)
- write (iout,*) "HHH H=",H,abs(HNose1-H0)/H0
-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
- 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,dimen
- vtnp_(i)=0.0d0
- vtnp_a(i)=0.0d0
- do j=1,dimen
- 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,dimen
- write (iout,'("WWW",i3,3f10.5)') i,vtnp(i),vtnp_(i),vtnp_a(i)
- enddo
-
- endif
- endif
- return
- end
projects / unres.git / blobdiff