How to solve for the Z-directional potential of the movable electrode due to the use of the cg method?

LAMMPS LAMMPS Beginners
1

Dear all,

I used the cg method in the electrode package of lammps-22Dec2022 to simulate the movable electrodes.
How can I get a Z-directional potential distribution for the system in this case,? My system is 2 graphene flat plate electrodes and ionic liquid electrolyte.

Thanks!
Best regards!

2

Hi! Unfortunately I have been swamped with work. As a last resort you can try the compute potential/atom in my USER-CONP2 package – it’s not very well documented, but it should work. (You will need to compile a separate old LAMMPS version for this – sadly, as an academic, I’m really not paid to develop code beyond the bare essentials.)

3

Dear professor,

Thank you for the solution! I’ll try it and come back to give feedback!

Thanks!
Best regards!

4

Dear professor,
I used lammps-22Dec2022 to run the dump trace then used lammps-27May2021 USER-CONP2’s cmpute potential/atom to rerun trj later for potential calculation. But I found that this path does not seem to work, as long as I use fix to output potential/atom, the program directly crashes when it comes to the output step. For example, when I use the following command, lammps stops at 500 steps.

in:
fix potential all ave/chunk 1 500 500 pos c_potential file z_potential.dat

Error：
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= PID 18410 RUNNING AT mu01
= EXIT CODE: 11
= CLEANING UP REMAINING PROCESSES
= YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
** Intel(R) MPI Library troubleshooting guide:**
** Documentation Library**
==================================================================

Is this route not possible, it is not desirable?

5

Hi,

This is too little information to find out what the error was. I myself haven’t used compute potential/atom in almost two years.

I am glad for your sake that it didn’t work, because I just remembered a big problem with that approach: while it returns correct potentials for the atoms in constant potential electrodes (for obvious reasons), you cannot simply calculate the potential at the location of a ionic liquid’s ions, because ions will usually sit near the local minima of electrostatic potential. You need to define a separate potential probe grid, and that is not something easily taught or troubleshooted purely over the internet.

For more ideas, the definitive reference on this is https://iopscience.iop.org/article/10.1088/0953-8984/28/46/464006/meta . Note most methods listed in the paper are based on a simple 1D charge density profile. If all else fails, you can just try rerunning your simulation trajectory using only a coul/long pair style, at which point the potential at a charge is simply the output of compute pe/atom divided by charge.

6

Dear professor,

Thank you for your valuable advice! I’ll try them all and come back with feedback!

Thanks!
Best regards!

7

Dear professor,

When I used the simplest approach so far: only coul/long and then Pe/q rerun my traces with 2.0V and 0.0V applied to the left and right electrodes respectively, I found the following problems.
in:
kspace_style pppm/electrode 1.0e-6 #
fix cpm EL electrode/conp 2.0 1.805 couple ER 0.0 algo cg 1e-6 symm on ffield on etypes on

in_rerun:
kspace_style pppm/electrode 1.0e-6 #
fix cpm EL electrode/conp 2.0 1.805 couple ER 0.0 algo cg 1e-6 symm on ffield on etypes on

variable q atom q
compute qleft EL reduce sum v_q #
compute qright ER reduce sum v_q
compute qbulk bulk reduce sum v_q

compute pe all pe/atom
compute peleft EL reduce sum c_pe #
compute peleright ER reduce sum c_pe #

variable potential atom c_pe/(v_q+1e-99) # avoid charge is 0

compute pleft EL reduce sum v_potential #
compute pright ER reduce sum v_potential
compute pbulk bulk reduce sum v_potential # compute pbulk bulk reduce sum v_potential

variable pleft1 equal c_peleft/c_qleft
variable pright1 equal c_peleright/c_qright
variable pbulk1 equal c_pleft/c_qbulk

thermo_style custom step time temp c_temphot c_tempcold etotal c_qleft c_qright c_peleft c_peleright c_pleft c_pright c_pbulk v_pleft1 v_pright1 v_ pbulk1#
thermo_modify flush yes
thermo 100 #
###############################

  1. The number obtained after summing by “variable potential atom c_pe/(v_q+1e-99)” (c_pleft c_pright c_pbulk) will be very strange and much worse, probably because I don’t know this command well enough, and “variable pleft1 equal c_peleft/c_qleft” shows a normal behavior.
  2. it is clear that the positive and negative voltage difference is 2.0V, but the voltage difference calculated by using “variable pleft1 equal c_peleft/c_qleft” is less than half, that is, only 1.0V (the same is true for the 1.0V trace, which is only 0.5V).
  3. When rerun 0.0V trajectory, there will be positive potential of both positive and negative poles, which is caused by the total charge and total potential of the positive and negative poles themselves.

The above is the result of my test with coul/long, rerun trajectory, I attach here 0.0V and 2.0V rerun.log, please enlighten me!

Thanks!
Best regards!

rerun0.0v.log (64.3 KB)
rerun2.0v.log (1.0 MB)

8

There are several issues here:

  1. Did you read the 2016 paper I linked to, and understand what it said? Have you tried calculating the charge profile and using that to calculate the potential? I explicitly said that you should try compute pe/atom if all else fails. See this for an example approach: Calculating 1D electrostatic potential profile from MD simulations – Alta Fang

  2. The sum of (atomic potential energy / atomic charge) is not the same as the total potential energy divided by the total charge.

  3. compute pe/atom does not know about the Gaussian charge correction used by ELECTRODE (but that’s a small issue)

  4. More importantly, you did not tell me you are using the finite field approach for constant potential. If you do, you need to both adjust the boundary conditions and include the electric field contribution to potential energy.

  5. Calculating potential at ionic positions using compute pe/atom seems to give very inaccurate results since monatomic ions will prefer to sit in local electrostatic minima. Granted, this effect isn’t very well described in the literature, but I mentioned it in my very last post.

  6. You need to consult either a local supervisor or a local expert in classical electrostatics to make sure you are getting correct results. Constant potential molecular dynamics is relatively untested and I have already seen outrageously wrong papers pass peer review. I cannot possibly check your work line by line, and even if I did, I am just a stranger on the Internet whom you shouldn’t necessarily trust. For example, did you know I’m not even a professor? :wink:

9

Dear professor,

Well I’m still used to calling you Professor.I certainly know that you are not a real professor. But in my opinion in the field of constant potential simulation, you already are. Maybe not now, in the future for sure.
I thank you for quickly pointing out the holes in my calculations with coul/long, and I’ll move on to try to integrate the potential using the charge density distribution (which I’ve actually done, and in a fixed electrode, he’s correct. (After getting to the electrode movable, I don’t know if it’s my code that makes the integration completely wrong, I’ll try again and come back with feedback!)

Thanks!
Best regards!