Ewald summation vs Wolf summation for thin films using core-shell interatomic potential

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS users and developers,

I am trying to run barium titanate thin film using this core-shell interatomic potential: https://doi.org/10.1063/1.4827475

The problem I ran into is that the thin film simulation crashes within a couple hundreds steps (Bond loss at the surface) using the ewald summation method, while it does not with the wolf summation method.

I have benchmarked the phase transition of the core-shell model with the bulk material and the phase transition was well reproduced. In this benchmark, the ewald summation was used to compute Coulombic interactions. I have tested the other bulk systems with the ewald summation and all showed stable energy profile.

The thin film system was created by adding vacuum along the z direction (following this paper: https://doi.org/10.1016/j.surfin.2024.105589). When testing with the ewald summation, the simulations crash within a couple hundreds steps with the error message (ERROR on proc 38: Bond atoms 7497 7498 missing on proc 38 at step 107 (…/ntopo_bond_all.cpp:63)). Visualization of the trajectories show that the core-shell bonds at the surface are broken off and ions are penetrating to the other ions’ sites.

When switching from the Ewald summation to the Wolf summation using the Wolf parameters from the papers (https://doi.org/10.1016/j.surfin.2024.105589, https://doi.org/10.1038/s41598-017-01002-0), the simulations don’t crash and the energy profiles look very stable.

Ultimately, I would like to use the ewald summation or pppm method, so I would love to understand what I am doing wrong with the ewald summation.

If anyone could help me understand why this is the case, it would be greatly appreciated.

Best regards,

Sue

PS. I am using the LAMMPS 2019 version.

Here’s the input file:

------------------------ INITIALIZATION ----------------------------

units metal
dimension 3
boundary p p p
atom_style full

----------------------- ATOM DEFINITION ----------------------------

read_data data.BTO_10_10_10

group cores type 1 2 3 #Ba Ti O
group shells type 4 5 6 #Ba Ti O

group Ba type 1
group Ti type 2
group O type 3
group Ba_s type 4
group Ti_s type 5
group O_s type 6

neigh_modify delay 0 every 1 check yes page 500000 one 50000
comm_modify vel yes # yes to communicate velocity info with ghost atoms

------------------------ FORCE FIELDS ------------------------------

K space: long-range Coulomb is computed in recipe. space

#kspace_style pppm 1.0e-6
kspace_style ewald 1.0e-6

Short-range interaction (Buckingham potential) only between shells of different atoms

Coulmb interaction between all cores and shells excpet the core and the shell in the same atom

#pair_style born/coul/wolf/cs 0.25 16.0 14.0
pair_style born/coul/long/cs 16.0 14.0

Born para:i j A_ij rho_ij sigma C_ij D_ij

pair_coeff * * 0.00 0.1000 0.00 0.00 0.00
pair_coeff 4 6 7149.81 0.3019 0.00 0.00 0.00 #Ba-O
pair_coeff 5 6 7200.27 0.2303 0.00 0.00 0.00 #Ti-O
pair_coeff 6 6 3719.60 0.3408 0.00 597.17 0.00 #O-O
pair_modify table 0

special_bonds coul 0.0 1.0 1.0

Core-shell bonds: k2r^2/2 + k4r^4/24

NOTE: values from reference paper must be rescaled to match

“class2” bond style, i.e. k2/2 and k4/24

bond_style class2 #r0, K2, K3, K4 for anharmonic spring

para: i r0 k2 k3 k4

bond_coeff 1 0.0 149.255 0.00 0.000 # Ba
bond_coeff 2 0.0 153.070 0.00 20.833 # Ti
bond_coeff 3 0.0 18.465 0.00 208.333 # O

------------------------ Equilibration Run -------------------------------

##COMPUTES

compute pe all pe/atom # potential energy
compute s all stress/atom NULL
compute u all displace/atom
compute CSequ all temp/cs cores shells
compute thermo_press_lmp all pressure thermo_temp

thermo_style custom step time temp press etotal pe vol lx ly lz press # print these very N step
thermo 1 # outptu thermodynamics every N steps
thermo_modify temp CSequ

#dump
dump dump1 all dcd 250 BTO.dcd
dump dump2 all atom 1 BTO.dump

velocity bias option

velocity all create 300 200 bias yes temp CSequ
velocity all scale 300 temp CSequ

timestep 0.0004 # 0.4 fs timestep
fix 1 all nvt temp 300 300 0.1
fix_modify 1 temp CSequ

------------------------ Dynamic Run -------------------------------

run 20000
undump dump1

2

I haven’t looked at this in depth, but you mention thin film, vacuum along the z direction, so it seems like perhaps you should be using the quasi-2D Ewald sum instead of 3D, see kspace_modify command — LAMMPS documentation, e.g. kspace_modify slab 3.0.

Also highly recommend updating your LAMMPS version to at least the latest stable in case relevant bugs were fixed in the last 5 years.

3

Thank you so much for the quick response. I will try the 2d ewald summation.

Besides, I was not able to attach the input files as I am a new user. I can send a data file via email if needed.

Best,

Sue

4

More fundamentally, Wolf charges do not have long range interactions. Any attempt to perform Ewald summation with charge values from a “Wolf force field” is technically using a completely different force field. The results should be somewhat similar for a suitable bulk homogeneous system, but any system with long range inhomogeneity or anisotropy should show different results.

Note that long-ranged electrostatics are not always automatically “more correct”. For example, the oxDNA2 model (pair_style oxdna2/excv command — LAMMPS documentation) represents coarse-grained DNA components with charges interacting through a truncated Debye-Huckel model, with an adjustable parameter for the solvent ionic strength. This interaction is not long-ranged. But it does correctly reproduce DNA structure and dynamics very well because on the target length scale, DNA charges are not long-ranged – they are screened by dissolved ions. Using Ewald summation or PPPM with oxDNA2 charge values would make the model worse, not better, while requiring more running time.

5

Hi,

Thank you so much for your insights. I have read the related literatures on the Ewald summation and Wolf summation approaches and understand that they are completely different.

I believe that the original force field (https://doi.org/10.1063/1.4827475) of this material is not built just for the Wolf summation. Unfortunately, the original paper doesn’t state which method is used to compute Coulombic interactions. However, this paper (https://doi.org/10.1038/s41598-017-01002-0) has benchmarked the Wolf summation parameters to reproduce the average Coulombic energy and pressure per ion when using the Ewald summation. Besides, I have seen other papers (e.g. https://doi.org/10.1063/5.0090231) using the Coulombic interaction for this shell model, although it wasn’t thin film.

Therefore, I believe that this shell model should be compatible with the ewald summation.

Best,
Sue