On Fri, 2015-07-31 at 15:48 +0800, whb wrote:
> Dear all
>
> When I restarted the pH-REMD simulation, I’ve got a problem about
> the temperature, it fell soon from 300K to 4K.
>
> At first it runs well when I restarted from the equilibrium step(runs in
> different pH, but did not exchange pH protonation with each other).And the
> temperature jumps between 260K to 300K, and I’ve accepted Jason’s
> explanation it’s the reason that the protein was frozen sometimes, and the
> temp seems a bit lower.
>
> I collected 10ns in the first simulation with Expilict solvent pH-REMD, it
> seems well in VMD.
>
> But when I want to restart this simulation using the same input file, the
> temp fell soon from 300K to 4K.
>
> I don’t know what the problem is and I need your help!
This sounds like exactly the same problem as before. Note that these
simulations have two different types of MC exchange attempts -- attempts
*within a replica* to change the protonation state of a titratable
residue and attempts to swap two different replicas.
The solute is frozen *only* following the protonation state change
attempts (not the replica exchange attempts). So even if replicas never
exchange with each other, any protonation state changes will cause the
temperature to be miscalculated (you would have to look at the cpout
file to determine if a protonation state change succeeded).
What's more is that I can't see anything that I'm doing wrong with the
velocities in the code itself. From
$AMBERHOME/src/pmemd/src/constantph.F90, lines 1437 to 1469 read like
this:
#ifdef MPI
call mpi_allgathervec(atm_cnt, vel)
#endif
do i = 1, cphfirst_sol - 1
vtemp(:, i) = vel(:, i)
vel(:, i) = 0.d0
end do
#ifdef CUDA
call gpu_upload_vel(vel)
#endif
! Call routine to do relaxation dynamics
#ifdef MPI
call relaxmd(atm_cnt, crd, mass, frc, vel, last_vel, my_atm_lst, &
mobile_atoms, ntrelax)
call mpi_allgathervec(atm_cnt, crd)
call mpi_allgathervec(atm_cnt, vel)
#else
call relaxmd(atm_cnt, crd, mass, frc, vel, last_vel, &
mobile_atoms, ntrelax)
#endif
! Restore the original velocities
#ifdef CUDA
call gpu_download_vel(vel)
#endif
do i = 1, cphfirst_sol - 1
vel(:,i) = vtemp(:,i)
end do
The first thing that happens in parallel is that all of the velocities
are distributed to all threads so everybody knows how fast every atom is
moving. Then, the X, Y, and Z-values of the velocity array is backed up
for the whole solute in a temporary array, then they are zeroed. The
relaxation dynamics routine is called, and both the coordinates and
velocities are gathered once again to all threads. Then, the velocities
for the solute are restored to the original velocity array.
That said, I think I've found the problem. Look at lines 1432 to 1438
of $AMBERHOME/src/pmemd/src/runmd.F90 (shown below):
#ifdef MPI
do atm_lst_idx = 1, my_atm_cnt
j = my_atm_lst(atm_lst_idx)
#else
do j = 1, atm_cnt
#endif
eke = eke + mass(j) * 0.25d0 * c_ave *&
((vel(1,j) + last_vel(1,j))**2 + &
(vel(2,j) + last_vel(2,j))**2 + &
(vel(3,j) + last_vel(3,j))**2)
end do
This approach to computing the velocities is done because with the
leapfrog integrator, the velocities are a half timestep ahead of the
positions, so the velocities at the same timestep as the positions are
taken as the average of the velocities at the half timestep before and
after the current positions (0.5^2 = 0.25, which is where the 0.25 comes
from). However, following the relaxation dynamics, the "last_vel" of
all of the solute atoms is 0, which means that the velocities at the
current step are estimated incorrectly -- they are half as large as they
*should* be, which will clearly reduce the temperature.
Fortunately, this quantity is computed for two purposes -- to compute
the temperature and to compute the temperature scaling factor that needs
to be used for the Berendsen thermostat. But you are using Langevin
dynamics, so the second one does not apply to you. As a result, the
temperature will be calculated wrong, but everything else will be fine.
I will work on a fix for this. In the meantime, just don't use the
Berendsen thermostat (but that is ALWAYS good advice :)).
All the best,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Jul 31 2015 - 05:30:02 PDT
