[lammps-users] Questions about temperature in an gas simulation

Dear Sir/Madam,

I am trying to simulate a simple gas with, made of let’s say A and B particles. The forces are derived from a Lennard-Jones potential between the pair A-B and a zero potential between similar pairs A-A and B-B.

The positions and the velocities of the atoms are generated for the input files. The velocities follow the Maxwell-Boltzammn distribution. The simulation runs fine but the temperature doesn’t correspond to the ensemble I have generated. Indeed it is systematically higher. I have check the temperature of my ensemble before running LAMMPS and it is correct, about T=296K while it is about T=320K in LAMMPS.

For the moment II suspect a problem with the units. In the documentation, it says that in the units electron, which is the one I use, the velocities are in Bohr/atomic time unit. So to build my input file, after having generated the velocities in m/s I convert them in Bohr /atomic time unit by simply multiplying them by (1.03275/5.29)*1e-4. Maybe I am doing something wrong here ?

I copy hereafter the code I am using so maybe it can help you to understand my problem.

units electron
atom_style atomic

read_data input/positions_initial.lammps

group mol1 type 1
group mol2 type 2

pair_style hybrid lj/cut 30.0 zero 30.0
pair_coeff 1 2 lj/cut 0.00064 7.058
pair_coeff 1 1 zero
pair_coeff 2 2 zero

timestep 1

fix 1 all nve

compute myRDF all rdf 100 1 2
fix 3 all ave/time 1 1 1 c_myRDF[1] c_myRDF[2] file rdf.dat mode vector

thermo 1000
run 10000

I hope that you can help to understand the problem.

Best regards

Albert,

There is not enough information to make any specific recommendations. For example, it would be extremely important to see a log file of your simulation.

I have two general points I would like to make and for you to consider:

  1. You have to look at the evolution of potential energy versus, kinetic energy and the total energy. Your initial configuration will most likely not be in equilibrium, which means that either some potential energy will be converted into kinetic energy (i.e. your temperature will rise) or vice versa
  2. You have to look at the stability of the time integration. I.e. is a timestep of 1 time unit providing a sufficient conservation of the total energy? Doing the simulation for “only” 10,000 steps may be too little to tell.
  3. Also, you have to look at the system size. Is your system large enough with respect to your conditions (temperature, density, pressure) to be a proper representation and to allow for a suitable average time between collisions/interactions between particles. If the system is too small, you may have unexpected behavior due to finite size effects.

Based on the limited information, point 1) is the most likely explanation, but please make some (careful) tests to confirm that this is the case before deciding on a remedy. You also may have a combination of the 3) options happening, or something else that I am currently not thinking of.

The suggestion of a remedy for point 1) would be to run with a thermostat added (i.e. replace fix nve with fix nvt) for a while until the expected temperature has reached an (time window averaged) value that you desired and then switch to fix nve.
The suggestion for point 2) would be to reduce the timestep.
The suggestion for point 3) would be to use a “replicate 2 2 2” to double the system in each direction and observe if results would improve in the desired direction.

Axel.

If the initial temperature does not match your computation, then most likely your computation is not using the same conversion factors as LAMMPS.

You can easily look up the details in the source code of “compute temp” and Update::set_units().

LAMMPS computes temperature from mv^2 multiplied by: mvv2e / (N_dof * k_B)
For simple atoms (point particles) in 3d N_dof is 3
natom - 3
For units electron:
mvv2e is 1.06657236
k_B is 3.16681534e-6