[ase-users] temperature setting for NPT calculations

[Ask Hjorth Larsen]

Dear Riping,

Den tir. 4. sep. 2018 kl. 11.59 skrev Wang Riping via ase-users
<ase-users at listserv.fysik.dtu.dk>:
>
> Dear ASE,
>
> I have error in setting temperature for NPT calculations.
>
> When I use the following lines to setting the temperature to be 300 K for NPT calculations, it gives the temperture in md.log file 0 K.
>
> from ase.md.npt import NPT
> dyn = NPT(atoms, 5*units.fs, 300, numpy.array([-100,-100,-100,0,0,0]), None, None)
>
>
> When I use the fillowing lines for NPT, the md.log file gives average tempearture about 150 K.
>
>
>
> from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
> MaxwellBoltzmannDistribution(atoms, 300 * units.kB)

This gives you a distribution with an instantaneous temperature of 300
K.  See below though.

> from ase.md.npt import NPT
>
> dyn = NPT(atoms, 5*units.fs, 300, numpy.array([-100,-100,-100,0,0,0]), None, None)

The NPT takes a temperature in eV.  You should probably put T = 300 *
units.kB here as well, else the temperature will drift to the wrong
value as the NPT compensates.

Now for the problem: You probably have a crystal.  It should have both
potential and kinetic energy when it moves 'naturally', but most
likely you gave it only kinetic energy.  The system will probably have
about equal parts kinetic and potential energy when oscillating
harmonically about its equilibrium positions.  When you only give it
kinetic energy, some of it will be converted to potential energy as
the system oscillates, and your actual temperature becomes half of the
'instantaneous' one you started with.

>
> Therefore, how to set the temperature for NPT?

In ase-3.16.2 you can set the starting temperature to 600 * units.kB
so it quickly goes down to 300 Kelvin.

In ase-3.17 (not released yet) you can calculate the phonon modes and
use the PhononHarmonics velocity distribution which displaces the
atoms /and/ sets the velocities according to thermally excited phonon
modes.  This instantaneously gives you the correct state.

Best regards
Ask

>
> Thank you very much.
>
>
> Riping Wang
>
> 2-18.9.4
>
>
>
>
>
> _______________________________________________
> ase-users mailing list
> ase-users at listserv.fysik.dtu.dk
> https://listserv.fysik.dtu.dk/mailman/listinfo/ase-users