Anomalous behaviour of PF6- with ECHEMDID

LAMMPS LAMMPS Beginners
,
1

I’m using LAMMPS version 23 Jun 2022, on Centos 6. I have installed the reaxFF package bundled with lammps and ECHEMDID-23jun2022

Hello, I’m trying to simulate the behavior of 1M of LiPF6 in Ethylene Carbonate under an electric field in order to establish a baseline for my future simulations. In this process I decided to use ReaxFF and ECHEMDID, as the force field and electric field implementations to run my simulations. However, upon equilibration, my PF6 molecule splits apart with all the flourine atoms repelled from the P center.

I would like assistance from people who have conducted similar simulations and whether this kind of behaviour has been observed earlier.

Here is a relevant snippet from my input deck

#################################################

#init

units real
atom_style charge
boundary p p s 
read_data electrolyte_LiPF6_40A.data
pair_style reaxff lmp_control 
pair_coeff * * ffield.reax.CHOSLiPF C C C F F F F F F H H H H Li O O O P 


#minimisation
fix chg all qeq/shielded 1 10 1.0e-10 500 reaxff
thermo 1
minimize 1.0e-12 1.0e-12 10000 10000

#ECHEMDID



variable zmin1 equal bound(all,zmin)
variable zmax1 equal ${zmin1}+0.2
variable zmax2 equal bound(all,zmax)
variable zmin2 equal ${zmax2}-0.2


fix             0 all property/atom d_locpot d_lap
compute         0 all property/atom d_locpot d_lap

region bound1 block INF INF INF INF ${zmin1} ${zmax1} units box
group z_left region bound1
region bound2 block INF INF INF INF ${zmin2} ${zmax2} units box
group z_right region bound2

set region bound1 d_locpot 0.0
set region bound2 d_locpot 0.0

fix 2 all echemdid 1 k 6.0 rc 4.0 norm 0.629967 nelec 25 z_left z_right volt 10.0

timestep 0.1
fix tres all nvt temp 1 300 100
run 30000
unfix tres
fix mynvt all nvt temp 300 300 100
run 50000
write_data electrolyte_min.data

################################################

Before_PF6
Before_PF61366×768 172 KB

After_PF6
After_PF61366×768 187 KB

electrolyte_LiPF6_40A.data (308.4 KB)
ffield.reax.CHOSLiPF (24.9 KB)

2

A straightforward guess is that the methodology has reacted poorly to such a strong electric field / large potential difference, mistakenly polarising the phosphorus from the fluorines in PF6 along the intended applied field, resulting in the topologies you’re seeing.

As to why that is happening, I don’t know. Could be any or all of: reaxff+echemdid simply won’t work for this system; your fix qeq needs to come after echemdid, not before; reaxff won’t work with this system because it doesn’t account for long range electrostatics; etc. etc.

A quick search shows that there are a few groups around the world who specialise in this technique – if you have more advanced questions, you may have to ask them directly.

3

Rather than simulating a big complex system, I suggest you run simulation of a small number of individual components (E.g. \textsf{PF}_6) by themselves and validate whether your force field and geometry is adequate. After that, I would try to set up a more complex geometry but still without the ECHEMDID add-on and see if that gives meaningful results, and only if that is successful, too, turn on the add-on feature. That way it is much easier to determine where the issue is (force field, geometry, combined structure, ECHEMDID add-on).

4

I’m also working on a similar topic, and I probably have even less knowledge about it than you do. But I found some minor errors in your code. The first is that Tdamp should be set to 100 times the size of timestep, i.e. 10, because setting it too big tends to cause crashes. Second, the voltage is set too high. Third, the charge-balance method qeq is not very suitable for the liquid model of applied electric field, and there are relevant literature says so. So some groups use acks2/reaxff, but I’m still exploring how this can be combined with echemdid.