Creating a potential profile from trajectory generated with ELECTRODE package

LAMMPS LAMMPS Beginners
1

Hi, I have run some simulations of mixtures of water and DMSO between two gold electrodes using fix electrode electrode/conp command (using the finite field algorithm). For example:

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

I wonder if it possible to extract potential (or electric field) profiles from the trajectories obtained. I have seen a previous discussion in the forum (How to solve for the Z-directional potential of the movable electrode due to the use of the cg method? - #7 by musong_mm) that seems to suggest it could be possible to extract them using the finite field method if taking care of the boundary conditions during integrations of the charge distributions. Is that correct? The variable charges of the electrode atoms are not needed?

Thanks for any feedback

Alfredo

2

In general, the electrostatic potential is obtained by solving the Poisson equation for \psi(z):

\partial_z^2 \psi(z) = \rho(z)

where \rho(z) is the 1D charge density (and some factor involving \epsilon_0 is added for units), obtained (for example) with fix ave/chunk (note you need ave none to average over time but not per particle).

Since this is a second order differential equation, it has a unique solution when combined with two boundary conditions. One such BC is always determined by the potential difference (which, for finite field, is applied at the simulation box boundaries and not the electrode positions) and the second BC corresponds to a potential offset for the whole system (which has no physical consequence).

In either case, the charge density includes the electrode charges. (How could it not?)

3

Thanks for the reply. I noticed in the user manual for ave/chunk there is not an “ave none” option but there is a “norm none”. It shouldn’t be the latter that needs to be added so the average is not divided by the number of atoms in the bin?

4

Yes, norm none – I was working off memory for that one.

5

I would to check if the following is correct. When integrating the Poisson Eq. once:

dPsi(z)/dz = - integral_0^z (rho(z’)/eps) dz’ + C1

The constant C1 should be C1 = deltaV/L_z (in the finite field method) where deltaV = V_right - V_left is 1.2 V in my case, and L_z the length of the simulation cell in the z direction. Is that correct?

For the second integration (to get Psi(z)) we get C1*z and a second constant C2 (that only produces a vertical shift of the potential, and therefore not physically relevant as you said before).
Thanks

6

That is the correct equation, yes