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

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

From: Matthias Kahk <matthiaskahk@...24...>
Date: Fri, 16 Jun 2017 15:31:34 +0800


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. 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 engaging tech sites,! _______________________________________________ lammps-users mailing list