# Pressure with an applied E-field

Hi,

I’ve run another set of test simulations to illustrate the problem.

The system consists of a gas of rigid, polar, diatomic molecules, in a fixed size box. The simulations were run in the NVE ensemble, either without an applied field, or with an applied field of 1.0 V/Angstrom along the z-axis. The simulation box contained 27 molecules (54 atoms) and the simulations were run over 20 million timesteps.

I calculated the pressure in three ways in both cases: (i) from the ideal gas law (ii) the value that LAMMPS prints under the thermo command (iii) by applying wall boundaries to all faces of the simulation cell, and calculating the average force per unit area from all of the walls. The results are given below.

Case 1: no external fields.

Temperature: 612 K
Pressure (ideal gas law): 2.85 x 10^4 Pa
Pressure (LAMMPS thermo): 2.88 x 10^4 Pa
Pressure (forces on the walls): 2.91 x 10^4 Pa

Case 2: 1.0 V/Angstrom external field.

Temperature: 627 K
Pressure (ideal gas law): 2.92 x 10^4 Pa
Pressure (LAMMPS thermo): -1.42 x 10^5 Pa
Pressure (forces on the walls): 2.96 * 10^4 Pa

Does anybody know of a way to get the thermo keyword to return physically meaningful pressures in simulations with an applied electric field?

Best Regards,

Matthias

A couple things you can try.

a) use flexible molecules instead of SHAKE, see if the results are different

b) define a compute pressure command that includes keywords that

do not include the contribution from fixes (fix shake in this case),

to see if a pressure w/out that term is what you are looking for

note that you can make the pressure from that compute you define

be output with thermo, or even replace the default value that “press”

reports (see the thermo_modify press command)

I’m CCing Aidan, as he may have better comments. He recently worked

on options to enable the forces from walls to be included in the virial/pressure,

so that is another option.

Steve

I agree with Matthias’ analysis. I was able to modify fix efield to recover P=0 by adding code to optionally add in the virial contribution from fix efield which exactly cancels the contribution from fix shake. This should be released shortly. The reason this will not be the default is that the decision to add or not add the virial contribution from a fix must be handled on a case by case basis, and so is best left up to the user.

Aidan

I agree with Matthias’ analysis. I was able to modify fix efield to recover P=0 by adding code to optionally add in the virial contribution from fix efield which exactly cancels the contribution from fix shake. This should be released shortly. The reason this will not be the default is that the decision to add or not add the virial contribution from a fix must be handled on a case by case basis, and so is best left up to the user.

Aidan