Hi,
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,
Matthias