+++ /dev/null
-c-------------------------------------------------------------------------------
- subroutine RESPA_step(itime)
-c-------------------------------------------------------------------------------
-c Perform a single RESPA step.
-c-------------------------------------------------------------------------------
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
-#ifdef MPI
- include 'mpif.h'
- integer IERROR,ERRCODE
-#endif
- include 'COMMON.SETUP'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
-#else
- include 'COMMON.LANGEVIN.lang0'
-#endif
- 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_short(0:n_ene),
- & energia_long(0:n_ene)
- double precision cm(3),L(3),vcm(3),incr(3)
- double precision dc_old0(3,0:maxres2),d_t_old0(3,0:maxres2),
- & d_a_old0(3,0:maxres2)
- integer ilen,count,rstcount
- external ilen
- character*50 tytul
- integer maxcount_scale /10/
- common /gucio/ cm
- double precision stochforcvec(MAXRES6)
- common /stochcalc/ stochforcvec
- integer itime
- logical scale
- common /cipiszcze/ itt
- itt=itime
- if (ntwe.ne.0) then
- if (large.and. mod(itime,ntwe).eq.0) then
- write (iout,*) "***************** RESPA itime",itime
- write (iout,*) "Cartesian and internal coordinates: step 0"
-c call cartprint
- call pdbout(0.0d0,"cipiszcze",iout)
- call intout
- write (iout,*) "Accelerations from long-range forces"
- 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
- endif
-c
-c Perform the initial RESPA step (increment velocities)
-c write (iout,*) "*********************** RESPA ini"
- call RESPA_vel
- if (ntwe.ne.0) then
- 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
- endif
-c Compute the short-range forces
-#ifdef MPI
- tt0 =MPI_Wtime()
-#else
- tt0 = tcpu()
-#endif
-C 7/2/2009 commented out
-c call zerograd
-c call etotal_short(energia_short)
-c call cartgrad
-c call lagrangian
-C 7/2/2009 Copy accelerations due to short-lange forces from previous MD step
- do i=0,2*nres
- do j=1,3
- d_a(j,i)=d_a_short(j,i)
- enddo
- enddo
- if (ntwe.ne.0) then
- if (large.and. mod(itime,ntwe).eq.0) then
- write (iout,*) "energia_short",energia_short(0)
- write (iout,*) "Accelerations from short-range forces"
- 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
- endif
- endif
-#ifdef MPI
- t_enegrad=t_enegrad+MPI_Wtime()-tt0
-#else
- t_enegrad=t_enegrad+tcpu()-tt0
-#endif
- 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)
- d_a_old(j,i)=d_a(j,i)
- enddo
- enddo
-c 6/30/08 A-MTS: attempt at increasing the split number
- do i=0,2*nres
- do j=1,3
- dc_old0(j,i)=dc_old(j,i)
- d_t_old0(j,i)=d_t_old(j,i)
- d_a_old0(j,i)=d_a_old(j,i)
- enddo
- enddo
- if (ntime_split.gt.ntime_split0) ntime_split=ntime_split/2
- if (ntime_split.lt.ntime_split0) ntime_split=ntime_split0
-c
- scale=.true.
- d_time0=d_time
- do while (scale)
-
- scale=.false.
-c write (iout,*) "itime",itime," ntime_split",ntime_split
-c Split the time step
- d_time=d_time0/ntime_split
-c Perform the short-range RESPA steps (velocity Verlet increments of
-c positions and velocities using short-range forces)
-c write (iout,*) "*********************** RESPA split"
- do itsplit=1,ntime_split
- if (lang.eq.1) then
- call sddir_precalc
- else if (lang.eq.2 .or. lang.eq.3) then
-#ifndef LANG0
- call stochastic_force(stochforcvec)
-#else
- write (iout,*)
- & "LANG=2 or 3 NOT SUPPORTED. Recompile without -DLANG0"
-#ifdef MPI
- call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
-#endif
- stop
-#endif
- endif
-c First step of the velocity Verlet algorithm
- if (lang.eq.2) then
-#ifndef LANG0
- call sd_verlet1
-#endif
- else if (lang.eq.3) then
-#ifndef LANG0
- call sd_verlet1_ciccotti
-#endif
- else if (lang.eq.1) then
- call sddir_verlet1
- else
- call verlet1
- endif
-c Build the chain from the newly calculated coordinates
- call chainbuild_cart
- if (rattle) call rattle1
- if (ntwe.ne.0) then
- if (large.and. mod(itime,ntwe).eq.0) then
- write (iout,*) "***** ITSPLIT",itsplit
- write (iout,*) "Cartesian and internal coordinates: step 1"
- call pdbout(0.0d0,"cipiszcze",iout)
-c call cartprint
- call intout
- 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
- endif
-#ifdef MPI
- tt0 = MPI_Wtime()
-#else
- tt0 = tcpu()
-#endif
-c Calculate energy and forces
- call zerograd
- call etotal_short(energia_short)
-#ifdef TIMING_ENE
-#ifdef MPI
- t_eshort=t_eshort+MPI_Wtime()-tt0
-#else
- t_eshort=t_eshort+tcpu()-tt0
-#endif
-#endif
- call cartgrad
-c Get the new accelerations
- call lagrangian
-C 7/2/2009 Copy accelerations due to short-lange forces to an auxiliary array
- do i=0,2*nres
- do j=1,3
- d_a_short(j,i)=d_a(j,i)
- enddo
- enddo
- if (ntwe.ne.0) then
- if (large.and. mod(itime,ntwe).eq.0) then
- write (iout,*)"energia_short",energia_short(0)
- write (iout,*) "Accelerations from short-range forces"
- 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
- endif
- endif
-c 6/30/08 A-MTS
-c Determine maximum acceleration and scale down the timestep if needed
- call max_accel
- amax=amax/ntime_split**2
- call predict_edrift(epdrift)
- if (ntwe.gt.0 .and. large .and. mod(itime,ntwe).eq.0)
- & write (iout,*) "amax",amax," damax",damax,
- & " epdrift",epdrift," epdriftmax",epdriftmax
-c Exit loop and try with increased split number if the change of
-c acceleration is too big
- if (amax.gt.damax .or. epdrift.gt.edriftmax) then
- if (ntime_split.lt.maxtime_split) then
- scale=.true.
- ntime_split=ntime_split*2
- do i=0,2*nres
- do j=1,3
- dc_old(j,i)=dc_old0(j,i)
- d_t_old(j,i)=d_t_old0(j,i)
- d_a_old(j,i)=d_a_old0(j,i)
- enddo
- enddo
- write (iout,*) "acceleration/energy drift too large",amax,
- & epdrift," split increased to ",ntime_split," itime",itime,
- & " itsplit",itsplit
- exit
- else
- write (iout,*)
- & "Uh-hu. Bumpy landscape. Maximum splitting number",
- & maxtime_split,
- & " already reached!!! Trying to carry on!"
- endif
- endif
-#ifdef MPI
- t_enegrad=t_enegrad+MPI_Wtime()-tt0
-#else
- t_enegrad=t_enegrad+tcpu()-tt0
-#endif
-c Second step of the velocity Verlet algorithm
- if (lang.eq.2) then
-#ifndef LANG0
- call sd_verlet2
-#endif
- else if (lang.eq.3) then
-#ifndef LANG0
- call sd_verlet2_ciccotti
-#endif
- else if (lang.eq.1) then
- call sddir_verlet2
- else
- call verlet2
- endif
- if (rattle) call rattle2
-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)
- d_a_old(j,i)=d_a(j,i)
- enddo
- enddo
- enddo
-
- enddo ! while scale
-
-c Restore the time step
- d_time=d_time0
-c Compute long-range forces
-#ifdef MPI
- tt0 =MPI_Wtime()
-#else
- tt0 = tcpu()
-#endif
- call zerograd
- call etotal_long(energia_long)
-#ifdef TIMING_ENE
-#ifdef MPI
- t_elong=t_elong+MPI_Wtime()-tt0
-#else
- t_elong=t_elong+tcpu()-tt0
-#endif
-#endif
- call cartgrad
- call lagrangian
-#ifdef MPI
- t_enegrad=t_enegrad+MPI_Wtime()-tt0
-#else
- t_enegrad=t_enegrad+tcpu()-tt0
-#endif
-c Compute accelerations from long-range forces
- if (ntwe.ne.0) then
- if (large.and. mod(itime,ntwe).eq.0) then
- write (iout,*) "energia_long",energia_long(0)
- write (iout,*) "Cartesian and internal coordinates: step 2"
-c call cartprint
- call pdbout(0.0d0,"cipiszcze",iout)
- call intout
- write (iout,*) "Accelerations from long-range forces"
- 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
- endif
-c Compute the final RESPA step (increment velocities)
-c write (iout,*) "*********************** RESPA fin"
- call RESPA_vel
-c Compute the complete potential energy
- do i=0,n_ene
- potEcomp(i)=energia_short(i)+energia_long(i)
- enddo
- potE=potEcomp(0)-potEcomp(20)
-c 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/(dimen3*Rb)*EK
-c Backup the coordinates, velocities, and accelerations
- if (ntwe.ne.0) then
- 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
- endif
- return
- end
projects / unres.git / blobdiff