LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] [EXTERNAL] Re: Pressure with an applied E-field
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] [EXTERNAL] Re: Pressure with an applied E-field

From: Matthias Kahk <matthiaskahk@...24...>
Date: Mon, 19 Jun 2017 09:26:38 +0800

Hi Aidan,

Thank you very much for looking into this and coming up with a solution. I guess I will keep an eye on new patches/releases on the LAMMPS website.

All the best,


On Sun, 2017-06-18 at 23:21 +0000, Thompson, Aidan wrote:

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 P. Thompson

      01444 Multiscale Science

      Sandia National Laboratories

      PO Box 5800, MS 1322      Phone: 505-844-9702

      Albuquerque, NM 87185     Fax  : 505-845-7442

      E-mail:athomps@...3... Cell : 505-218-1011



From: Steve Plimpton <sjplimp@...24...>
Date: Friday, June 16, 2017 at 7:50 AM
To: Matthias Kahk <matthiaskahk@...24...>
Cc: LAMMPS <>, "Thompson, Aidan" <athomps@...3...>
Subject: [EXTERNAL] Re: [lammps-users] Pressure with an applied E-field


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.



On Fri, Jun 16, 2017 at 1:31 AM, Matthias Kahk <matthiaskahk@...24...> wrote:



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,





On Fri, 2017-06-16 at 00:37 +0800, Matthias Kahk wrote:

Hi Steve,
Thank you for your response. I understand that this is what is going on in the simulation, but I'm not sure whether it makes physical sense. Consider a simulation of some CO gas in a fixed size container. If a potential difference were applied between two opposite faces of the container, most of the molecules would tend to align with the field. By extrapolating the results from the single-molecule simulation, there would then be a large (negative!) contribution to the total pressure, that has nothing to do with the interactions between the molecules, or the velocities of the gas molecules. 
I wouldn't expect a calculated "pressure" that contains this term to be meaningful in the conventional sense, i.e. in terms of the average force per unit area excerted by the gas molecules on the container walls, but I'm open to being proven wrong. (In fact, I might give this simulation a go and report on the results...)
On Thu, 15 Jun 2017 08:06:26 -0600
Steve Plimpton <sjplimp@...24...> wrote:
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.
On Wed, Jun 14, 2017 at 8:18 PM, Matthias Kahk <matthiaskahk@...24...>
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
Best Regards,
Check out the vibrant tech community on one of the world's most
engaging tech sites,!
lammps-users mailing list