Temperature in a NVE simulation

Qijing Zheng

unread,

Sep 6, 2019, 7:49:40 AM9/6/19

to ipi-users

Hi Everyone

I was doing an test NVE molecular dynamics using i-pi-driver as the force engine. Below is the "input.xml",

where I set the initial temperature to be 300 K.

###################################################################################

DEPENDENCIES h5o2.dms4B.coeff.com.dat h5o2.pes4B.coeff.dat init.xyz

COMMAND(4) i-pi-driver -u -h REGTEST_SOCKET -m zundel

ENDREGTEST-->

</ffsocket>

<total_steps>100</total_steps>

<trajectory stride='1' filename='pos' cell_units='angstrom'>positions{angstrom}</trajectory>

[ step, time{femtosecond}, temperature{kelvin}, conserved{electronvolt},

kinetic_md{electronvolt}, potential{electronvolt}]

</properties>

</output>

<prng>

</prng>

</forces>

<cell>

[ 25.29166, 0, 0, 0, 25.29166, 0, 0, 0, 25.29166 ]

</cell>

</initialize>

</dynamics>

</motion>

</system>

</simulation>

###################################################################################

I am a little confused with the temperature written in the output file, as can be seen below, the temperature

is as high as 4.53115807e+04 K. Shouldn't the temperature in the first few steps be around 300K?

===============================================================================

# column 1 --> step : The current simulation time step.

# column 2 --> time{femtosecond} : The elapsed simulation time.

# column 3 --> temperature{kelvin} : The current temperature, as obtained from the MD kinetic energy.

# column 4 --> conserved{electronvolt} : The value of the conserved energy quantity per bead.

# column 5 --> kinetic_md{electronvolt} : The kinetic energy of the (extended) classical system.

# column 6 --> potential{electronvolt} : The physical system potential energy.

0.00000000e+00 0.00000000e+00 4.53115807e+04 3.11658082e-01 1.81791784e-01 1.05496586e-02

1.00000000e+00 2.50000000e-01 4.53066247e+04 3.11692143e-01 1.77307491e-01 1.50680129e-02

2.00000000e+00 5.00000000e-01 4.52995597e+04 3.11736260e-01 1.70914970e-01 2.15046503e-02

3.00000000e+00 7.50000000e-01 4.52911454e+04 3.11785128e-01 1.63301493e-01 2.91669957e-02

===============================================================================

Best Regards

Qijing

Qijing Zheng

unread,

Sep 6, 2019, 10:24:56 PM9/6/19

to ipi-users

In the output properties, e.g. at the first step, the temperature is 4.53115807e+04 while the kinetic_md is 1.81791784e-01eV.

If I am correct, then the kinetic_md corresponds to 281.28 K instead of 4.54E4 K. Here is how I do the conversion

2 * kinetic_md / kB / (natoms * 3) = 2 * 1.82E-1 / 8.6173E-5 / (5 * 3) = 281.28 K

where kB is the boltzmann constants.

Am I doing something wrong in the conversion? Or did I misunderstand the output properties?

Please, can someone help me with this problem? Any help would be appreciated!

Thank you!

Qijing

Michele Ceriotti

unread,

Sep 7, 2019, 3:43:41 AM9/7/19

to ipi-users

OK, this is complicated.

* On a practical level:

if you want to get meaningful temperature numbers, either set a meaningful temperature in the ensemble, or set <fixcom> False </fixcom> in the <motion> block

* On a fundamental level:

temperature is kind of ill-defined in the NVE ensemble, or better the kinetic energy estimator for the temperature assumes you have a canonical ensemble

* On an implementation level:

defining a temperature estimator that gives you the correct ensemble temperature in NVT is kind of a pain when you have ensembles with fixed atoms constraints. if you look at get_temp in properties.py you'll see it's quite involved, at it requires making assumption of what is the ensemble temperature. if you don't set an <ensemble>, it assumes the temperature is 1 atomic unit, which is a lot. it is hard to change the behavior in the current implementation between ensembles, because the ensemble you're sampling depends on the integrator, and due to the fundamental issues above I've no intention of rehauling the code to introduce a definition that would not be any more "right".

Qijing Zheng

unread,

Sep 7, 2019, 4:11:52 AM9/7/19

to ipi-users

Dear Prof. Ceriotti

Thank you very much for your detail and helpful reply!

- I have set the "fixcom" tag as you suggested and get a somewhat meaningful temperature, the order of magnitude is correct.

# [ step, time{femtosecond}, temperature{kelvin}, conserved{electronvolt}, kinetic_md{electronvolt}, potential{electronvolt}]

0.00000000e+00 0.00000000e+00 3.32782334e+02 3.11658082e-01 3.01108423e-01 1.05496586e-02

1.00000000e+00 2.50000000e-01 3.27826427e+02 3.11692228e-01 2.96624215e-01 1.50680127e-02

2.00000000e+00 5.00000000e-01 3.20761642e+02 3.11736524e-01 2.90231850e-01 2.15046739e-02

3.00000000e+00 7.50000000e-01 3.12347254e+02 3.11785434e-01 2.82618335e-01 2.91670989e-02

- However, the kinetic_md and the temperature still do not match. This confuses me since

in the "get_temp" method, you call "get_kinmd" to get the kinetic energy and return the

temperature by dividing the number of degrees of freedom. How can the kinetic_md and the

temperature differ? I think that there must be something that I missed, please enlighten me!

################################################################

kemd, ncount = self.get_kinmd(atom, bead, nm, return_count=True)

if self.motion.fixcom:

# Removes the fake momentum from the centre of mass.

self.beads.p -= self.beads.m3 * vcm

if len(self.motion.fixatoms) > 0:

# re-fixes the fix atoms

for i in self.motion.fixatoms:

self.beads.p[:, 3 * i:3 * i + 3] = 0.0

return 2.0 * kemd / (Constants.kb * 3.0 * float(ncount) * self.beads.nbeads)

################################################################

Thank you!

Qijing

Michele Ceriotti

unread,

Sep 7, 2019, 5:03:53 AM9/7/19

to ipi-users

what do you mean? 3.01108423e-01 eV /(1.5*7*kb) = 332.78279 K

Qijing Zheng

unread,

Sep 7, 2019, 5:16:49 AM9/7/19

to ipi-users

Dear Prof. Ceriotti

- Please forgive me for my stupidness, I just found that I counted the number of atoms wrong, which is 7 instead of 5 :p

- I now sees that how the "fixcom" tag affects the calculation of the kinetic energy:

#######################################################################

if self.motion.fixcom:

# Adds a fake momentum to the centre of mass. This is the easiest way

# of getting meaningful temperatures for subsets of the system when there

# are fixed components

.....

#######################################################################

Thank you again for you kindness and valuable time!

Best

Qijing

Reply all

Reply to author

Forward