Constant potential method (CPM) Code in ELECTRODE Package

LAMMPS LAMMPS General Discussion
,
1

Dear Dr Tee,

I hope this message finds you well.

I am Lingzhi Cao, a PhD student in University of Melbourne. I am currently working with the ELECTRODE package in LAMMPS to use the constant potential method for controlling the electric potential at a solid–liquid interface. However, I encountered some discrepancies when attempting to extract the electric potential, and I would greatly appreciate your insights on this.

Here is the procedure I followed:

I began by running the au-aq test case provided in the package. I first ran the simulation using the in.ffield input script to get a trajectory file, and then performed a rerun on this trajectory while switching off vdW interactions (i.e., only retaining the coul/long interactions). During this rerun, I extracted the electrostatic potential energy of each atom using compute pe/atom pair kspace. My intention was to compute the local electric potential by dividing this energy by the corresponding atomic charge.

However, I noticed a discrepancy: the electric potential values derived in this way do not match the applied values specified in the original in.ffield input. Specifically, I had set a potential difference of 2 V between the two electrodes using the command:

fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on

If I interpret the syntax correctly, this should enforce -1 V on the bottom electrode and +1 V on the top. But the potential values I obtained from my postprocessing instead showed -0.19 V for the negative electrode and +0.39 V for the positive one. Interestingly, atoms within each metal surface do appear to share a consistent potential, but the absolute values differ from what I had expected.

Could you kindly advise if there is anything I may have misunderstood in either the interpretation of the electrode/conp settings or the methodology I used for extracting the electric potential?

To help clarify my issue, I’ve uploaded my input scripts to dropbox (https://www.dropbox.com/scl/fi/b14qbtpkdhwgkymcivwqu/input_scripts_lc.tar.xz?rlkey=6tgdmfkp75m1kwmvx6rnmlts7&st=w6t6c4d9&dl=0). I am using LAMMPS version 29Aug2024, compiled with Intel OneAPI and the Intel acceleration package.

Thanks @srtee very much for your excellent work on the ELECTRODE package. I’m truly grateful for your time and any guidance you can provide.

Kind regards,

2

Hi Lingzhi,

Thank you for your interest in our code. Please double check that (1) the same electric field has been applied in your rerun as in the original run; (2) the electric field contribution is included in the per-atom coulombic potential energy contribution. (For example you could artificially impose different fix efield magnitudes in the reruns and see if that influences the readout of compute pe/atom).

3

Dear Dr Tee,

Thank you very much for your prompt reply.

I just wanted to clarify one point in case I didn’t express myself clearly earlier. In both the initial run and the rerun, I did not use fix efield to apply any explicit external electric field. Instead, I relied solely on the fix electrode/conp command to impose the constant potential boundary condition across the two electrodes.

My main goal was to verify whether the local electric potential at each electrode surface could be recovered from the per-atom electrostatic potential energy (from compute pe/atom pair kspace) divided by the atomic charge. I performed a rerun with van der Waals interactions removed, keeping only coul/long, and extracted the potential energy this way.

Given this setup, I had expected to recover values close to –1 V and +1 V (as specified by fix electrode/conp). However, the results I obtained instead showed approximately –0.19 V for the negative electrode and +0.39 V for the positive one. This leads me to wonder whether:

  1. The effect of the imposed boundary condition from electrode/conp is retained during rerun, or
  2. The compute pe/atom command includes this contribution when evaluating the electrostatic energy.

If you happen to have a moment, I would really appreciate it if you could take a quick look at my input script. I’ve included the relevant files below — the system is small and should run fairly quickly, and the discrepancy I mentioned appears quite clearly in the output.

Here is the input script for your reference (can also find them on dropbox):

# Input script
boundary p p p # ffield uses periodic z-boundary and no slab
units real
atom_style full
pair_style lj/cut/coul/long 15
bond_style harmonic
angle_style harmonic
kspace_style pppm/electrode 1e-7

read_data "data.au-aq"

group bot type 6
group top type 7

group SPC type 1 2 3
group electrolyte type 1 2 3 4 5

fix nvt electrolyte nvt temp 298.0 298.0 241
fix shake SPC shake 1e-4 20 0 b 1 2 a 1

variable q atom q
variable qz atom q*(z-lz/2)
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qztop top reduce sum v_qz
compute qzbot bot reduce sum v_qz
compute ctemp electrolyte temp

fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on

thermo 50
thermo_style custom step temp c_ctemp epair etotal c_qtop c_qbot c_qztop c_qzbot
dump            1 all  custom  1000 flow.lammpstrj id type x y z q
dump_modify     1 sort id
run 5000


# Rerun script

boundary p p p # ffield uses periodic z-boundary and no slab
units real
atom_style full
pair_style coul/long 15
bond_style harmonic
angle_style harmonic
kspace_style pppm/electrode 1e-7

read_data "data_no_lj.au-aq"
pair_coeff * *

group bot type 6
group top type 7

group SPC type 1 2 3
group electrolyte type 1 2 3 4 5

fix nvt electrolyte nvt temp 298.0 298.0 241
fix shake SPC shake 1e-4 20 0 b 1 2 a 1

compute potential_per_atom all pe/atom pair kspace

fix conp bot electrode/conp -1.0 1.805132 couple top 1.0 symm on ffield yes etypes on

thermo 50
thermo_style custom step temp
dump            1 all  custom  1000 rerun.lammpstrj id type x y z q c_potential_per_atom
dump_modify     1 sort id

rerun           ./flow.lammpstrj every 1000 dump x y z q
run 0

Thank you again for your time and your excellent work on the ELECTRODE package.

Kind regards,
Lingzhi

4

Hi Lingzhi,

The finite field method imposes an electric field across the box to achieve the specified potential difference – see the documentation under “ffield on” for details and references. For your particular simulation the name of that fix would be conp_efield (in general the fix is named by adding _efield to the CONP fix’s name).

Now fix efield does not add energy contributions by default (you need to fix_modify ... energy yes) and your compute pe/atom does not currently count fix contributions (you need to add the fix keyword to pair and kspace). Once you have done that you should get the expected results.

5

Hi Dr Tee,

Thanks again for your clarifications—your explanation about how ffield yes works and how the conp_efield fix contributes to the potential energy was very helpful.

To isolate the issue more clearly, I decided to test this in a simpler setup—the planar example provided with the ELECTRODE package—since it avoids complications from van der Waals interactions. As per your suggestion, I added:

fix_modify conp_efield energy yes

and used:

compute pe/atom all pe/atom pair kspace fix

to include the fix contribution in the electrostatic energy. I then calculated the electric potential per atom by dividing the per-atom energy by the atomic charge.

However, the result still seems to deviate from the expected value. For example, in the 5th step (which I believe corresponds to a surface potential of 1 V), the computed potential on the Au surface atoms is around 0.1244 V—well below the expected value.

I’m attaching below the input script I used for this planar system test. The simulation is quite lightweight and reproduces the issue consistently:

boundary p p p # finite field, fully periodic
kspace_style pppm/electrode 1.0e-7

# set boundary in main script because ffield is periodic
units real
# distribute electrode atoms among all processors:
if "$(extract_setting(world_size) % 2) == 0" then "processors * * 2"

atom_style full
pair_style lj/cut/coul/long 14
# kspace_style and _modify in main script to test different approaches

read_data "data.au-vac"

group bot molecule 1
group top molecule 2

# ramping potential difference
variable v equal ramp(0,2)

# get electrode charges
variable q atom q
compute qbot bot reduce sum v_q
compute qtop top reduce sum v_q

# get theoretical charges:
# calculate distance dz between electrodes
compute zbot bot reduce max z
compute ztop top reduce min z
variable dz equal c_ztop-c_zbot

# calculate theoretical capacitance as eps0 * area / dz
variable eps0 equal 55.26349406/10000 # epsilon zero
variable capac equal "v_eps0 * lx * ly / v_dz"

# calculate theoretical charges and deviation of constant potential charges from theory
variable qtheory equal "v_capac * v_v"
variable percdev equal "100 * (c_qtop - v_qtheory) / (v_qtheory + 1e-9)" # avoid divide-by-zero

# constant potential electrodes with ramping potential difference
fix conp bot electrode/conp 0 1.979 couple top v_v symm on ffield yes
fix_modify conp_efield energy yes # As suggested

thermo 1
# thermo: step, imposed potential, bottom charge, top charge, theory charge, percent deviation
thermo_style custom step v_v c_qbot c_qtop v_qtheory v_percdev

compute potential_per_atom all pe/atom pair kspace fix
variable calculated_V atom c_potential_per_atom*0.043/(v_q+1e-8) # unit conversion and avoid divide-by-zero error

dump            1 all  custom  1 ffield_on.lammpstrj id type x y z q c_potential_per_atom v_calculated_V
dump_modify     1 sort id
run 10

I would greatly appreciate if you could have a look and let me know if I’ve missed anything. Thanks again for your time and continued support.

Kind regards,
Lingzhi

6

Hi Lingzhi,

Thank you for taking time to continue testing. This is a very helpful issue to consider as my collaborators in Hamburg (who honestly do more of the coding work than me) work on improving the ELECTRODE package, namely that the electrostatic potential is not easy to extract as a compute right now.

The main source of discrepancy is, surprisingly, that fix efield does not calculate per-atom energy contributions, only an overall contribution. I can think of good technical reasons for this (to do with periodic boundaries) but it still seems like an oversight to be fixed some day. In the meantime you can still add it manually to your calculations with the variable

variable efield_pot atom -v_conp_ffield_zfield*z

which returns the field-associated potential directly in volts (or generally in units of electric field * distance based on the units setting).

This will make up most of the discrepancy for the aqueous (au-aq) case. It still leaves a significant difference for the planar case because of the next reason: the electric potential of the Gaussian charges is not properly adjusted for their Gaussian shapes.

Historically, the first versions of fix conp focused on charge equilibration and the subsequent dynamics and did not have bells and whistles like per-atom energy correctly added in, which has persisted until today. This was compounded by most tests (including mine) focusing on planar electrodes, for which the most efficient way to calculate the potential profile is still calculating the transverse charge density and double-integrating (by Poisson’s equation), rather than per-particle “potential probes”. In realistic simulations the electrolyte charges are never greatly affected by the Gaussian shape of the electrode charges, and therefore the effect of Gaussianity on such potential profiles is minimal.

In the aqueous au-aq case, the majority of potential on the electrode charges is due to nearby electrolyte point charges, so after accounting for electric field you should recover the stated potential difference (with minor deviations due to Gaussianity). In the planar case, with only Gaussian charges, it is not currently possible to directly recover the correct electric potential from the standard LAMMPS pe/atom calculation.

Note that this is only an incorrect representation of the in-plane, per-atom potential energy – Gaussian forces are correctly imposed and the overall system energy correctly adjusted. Also, a quick direct calculation shows that the constant potential code charges the planar capacitor with the correct charge density of two plates separated by vacuum (indeed, historically, \eta was chosen precisely to minimize the discrepancy between a plane capacitor and the atomistic representation, rather than any sophisticated quantum chemistry parametrization!).

I hope this makes sense to you. In the long term, the ELECTRODE team’s plan is to directly account for Gaussian charges in pair styles (as they should) and make pair styles “responsible” for consistently calculating electric potential and electric potential energy, at which point you should get consistent results with your method (which is theoretically sound).