I've been going through pmemd source a bit. Temperature is definitely
calculated from kinetic energy, number of degrees of freedom and the
Boltzmann constant, but I am still not sure about the details.
Especially, I am a bit confused by an apparent distinction between
solvent and solute atoms during temperature calculation.
The correct relation should be
E_kin = N_dof * T/2 * k_B
with E_kin being the sum of all atomic kinetic energies and N_dof being
the number of translational degrees of freedom in the system (being 3
times the number of atoms in the most simple case).
The kinetic energy in pmemd is stored to `si(si_solute_kin_ene)` via
do j = 1, atm_cnt
eke = eke + mass(j) * &
(vel(1,j) * vel(1,j) + vel(2,j) * vel(2,j) + vel(3,j) * vel(3,j))
end do
eke = 0.5d0 * eke
! New addition - is this right?
si(si_solute_kin_ene) = eke
Why is this called `...solute...`? How and why is pmemd distinguishing
between solvent and solute? Anyway, this is then copied to `si(si_kin_ene)`:
si(si_kin_ene) = si(si_solute_kin_ene)
The temperature is calculated via:
temp = si(si_kin_ene) / fac(1)
Here, `fac(1)` is
fac(1) = boltz2 * rndf
This is the point I am not sure anymore. For the simple formalism, I've
mentioned above, `fac(1)` must be 3N * k_B / 2.
I.e. `boltz2` must be k_B/2 and `rndf` should be the correct number of
degrees of freedom. From the code, indeed, `boltz2` presumably is half
the Boltzman constant in MD-internal units. `rndf` is the number of
degrees of freedom, which I cannot really follow. It is mainly
calculated here:
degcnt.f90:
[...]
if (ibelly .le. 0) then
rndfp = 3*nsolut - rstssl
rndfs = 3*(nat - nsolut) - rstssv
else
rndfp = 3*ibelsl - rstssl
rndfs = 3*ibelsv - rstssv
end if
runmd.f90:
rndf = rndfp + rndfs ! total degrees of freedom in system
! Correct the totaldegrees of freedom for extra points (lone pairs):
rndf = rndf - 3.d0 * dble(numextra)
! BUGBUG - NOTE that rndfp, rndfs are uncorrected in an extra points
context!
Yeah.. I did not exactly try to follow the DOF calculation. I am still
interested why the solvent/solute atoms in the system are considered
differently. I am wondering, if the temperature in an explicit solvent
simulation is calculated from the solute atoms only or from solute and
solvent atoms. Also, I am wondering about this "BUGBUG" note.
If someone is willing to comment on these things, I would be happy. On
the other hand, it's not crucial for me :)
Thanks,
Jan-Philip
On 01/05/2012 09:50 PM, Brian Radak wrote:
> I can't answer your question for sure, but I do know that the instantaneous
> temperature depends on the number of degrees of freedom of the system (so
> it depends on connectivity, constraints, etc.).
>
> I don't know of any instantaneous temperature formalism other than average
> kinetic energy, so I would boldly assume that that is the one at play.
> Although that's not particularly helpful because of at least two other
> things I don't know:
>
> 1.) which velocities are used? since AMBER uses leap frog verlet it might
> depend on the thermostat?
>
> 2) the averaging could probably be done in a few, possibly nonequivalent,
> ways. I don't think averaging over atomic centers and molecular centers is
> the same (geometric vs center of mass velociites?)
>
> The fastest way to find out might be to check the source code.
>
> Brian
>
> On Thu, Jan 5, 2012 at 12:01 PM, Jan-Philip Gehrcke<jgehrcke.googlemail.com
>> wrote:
>
>> Huhu,
>>
>> sander/pmemd calculated the temperature of the simulated system for a
>> given point in time and prints it to e.g. the mdout file.
>>
>> I just want to make sure.. is this value calculated from the average
>> kinetic energy of all atoms in the system and the Boltzmann constant or
>> is some other formalism applied?
>>
>> Thanks,
>>
>> Jan-Philip
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>
>
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Jan 06 2012 - 07:30:02 PST
