Hi,
What version of cpptraj are you using? Check 'cpptraj
--internal-version' ; if it is 4.23.1 or earlier, there was a bug in
the temperature calculation where removal of global degrees of freedom
was not properly taken into account. If I do a test with the current
cpptraj (from GitHub or AmberTools 20) things seem ok, see below.
-Dan
```
CPPTRAJ: Trajectory Analysis. V4.29.0 (GitHub)
...
[setvelocity ntc 2]
Random_Number: seed is <= 0, using wallclock time as seed (6139)
SETVELOCITY: Assigning velocities for atoms in mask '*'
Temperature= 300.00, using Maxwellian distribution.
Constraints on bonds to H
Time step= 0.002 ps, epsilon = 1e-07
[temperature MyTemp ntc 2]
TEMPERATURE: Calculate temperature for atoms in mask [*]
Constraints: bonds to H
[run]
...
ACTION SETUP FOR PARM 'tz2.parm7' (2 actions):
0: [setvelocity ntc 2]
Mask [*] corresponds to 223 atoms.
Constraints on 106 bonds to hydrogen.
# of degrees of freedom = 563
1: [temperature MyTemp ntc 2]
Mask [*] corresponds to 223 atoms.
Constraints on 106 bonds to hydrogen.
# of degrees of freedom = 563
...
[printdata MyTemp]
#Frame MyTemp
1 303.2060
```
On Tue, Aug 11, 2020 at 7:46 PM German P. Barletta <pbarletta.gmail.com> wrote:
>
> Thank you for your answers. I tried setting the temperature with ptraj but
> ran into a new problem.
>
> Taking the 7th frame and assigning new velocities to the whole system and
> letting cpptraj know I intend to use SHAKE for all hydrogens:
> ----
> parm lb3f.prmtop
> trajin lb3f.nc 7 7
> setvelocity ntc 2
> trajout a.rst7
> go
> ----
>
> Checking everything went ok:
>
> ----
> parm lb3f.prmtop
> trajin a.rst7
> temperature ntc 2 out temp
> go
> ----
>
> But the bad news is temp file looks like this:
>
> ----
> #Frame Tdata_00001
> 1 213.6664
> ----from
>
> I checked the source code and *setvelocity* looks fine but *temperature* does
> a weird thing:
> ------------
> double Frame::CalcTemperature(AtomMask const& mask, int deg_of_freedom)
> const {
> ...
> double fac = (Constants::GASK_KCAL/2.0) * deg_of_freedom; // Estimate for
> DoF
> ...
> total_KE *= 0.5;
> return total_KE / fac;
> ------------
>
> If I'm not mistaken, the 2.0 and the 0.5 cancel each other.
>
> Being suspicious of the *temperature* function, I launched a replica from
> my *a.rst7* frame without restarting simulation and reading the velocities.
> I also set *ntpr = 1* so I could get the temperature after the 1st step.
> Sure enough, *pmemd.cuda* agrees with *setvelocity*:
> ---
> NSTEP = 1 TIME(PS) = 0.002 TEMP(K) = 303.02 PRESS =
> 0.0
> ---
>
> So, I think there's a bug in the *temperature* function.
> _______________________________________________
> 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 Wed Aug 12 2020 - 06:30:03 PDT
