LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Pressure with an applied E-field

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

 From: Matthias Kahk Date: Fri, 16 Jun 2017 15:31:34 +0800

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

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...)

Matthias

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.

Steve

On Wed, Jun 14, 2017 at 8:18 PM, Matthias Kahk <matthiaskahk@...24...>
wrote:

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

------------------------------------------------------------
------------------
Check out the vibrant tech community on one of the world's most