Pressure with an applied E-field


I have been running some MD simulations using fix efield, and noticed that “something strange” was going on with the pressures that I was getting out… so I did some tests.

Consider a single, rigid, diatomic molecule, aligned with the z-axis, in a non-periodic simulation cell, at zero temperature. (I.e. all of the atoms are at rest.)

The molecule carries no net charge, but the individual atoms have charges of +1e and -1e respectively. There are no forces defined between the atoms, and the geometry of the molecule is kept constant using fix shake.

What do I do?

I run the simulation, from this starting point, using the NVE time integrator, and monitor the pressure.

What do I expect?

That the pressure should be zero, because the atoms are at rest, and there are no forces between them. Even under an applied (constant, uniform) field, a rigid dipole could experience a torque, but no net force.

What do I observe?

That the pressure returned by the “thermo” keyword scales linearly with the applied field.

Why might this be?

From reading the LAMMPS manual and previous emails on this list, I have gathered that (i) External forces are not included in the virial of the pressure. (ii) Forces due to fixes that apply constraints are. Whilst I am not questioning the reasons for which these choices where originally made, I believe that in the current scenario this results in unphysical behaviour. (Pressure in a box with one rigid molecule at rest.)

Is pressure a well-defined quantity in simulations under an applied electric field?

Maybe not in general, however… in a system that satisfies the following criteria: (i) no “free” charge carriers (ii) homogeneous (iii) electric field is constant (iv) electric field is uniform; I cannot see why it shouldn’t be.

What would be a solution?

It might be that including the forces due to the electric field in the virial of the pressure would solve the problem. Alternatively, a more advanced treatment of the shake forces in which the response to particle-particle interactions can be distinguished from the response to external forces could probably achieve the same result.

Am I thinking along the right lines here, and do you think that a solution would be feasible? Any help would be greatly appreciated.

I can post the input and output files of my test job if they are of interest.

Best Regards,


Not sure I undestand. But if you apply an electric field

to a diatomic with + and - charges, then the larger the

field the more the 2 particles will want to violate the

SHAKE constraint, i.e. the farther they would move apart

or together on a timestep. So the SHAKE constraint

has to apply larger forces to keep them at the proper

distance. Hence the internal virial due to SHAKE increases.

Hence the pressure goes up.