The induced charge density obtained by fix polarize is always 4π times less than numerical results

Hi everyone, recently I am trying to use DIELECTRIC package in LAMMPS. I am trying to calculate the induced charges on the interface plane. The interface between two media: one with dielectric constant of 80 (epsilon2), the other of 5 (epsilon1) in my in file. The fixed charge is placed in the continuous medium of epsilon2 with a distance of 6 angstroms from the interface.I have tried different epsilon values, different charges and different distances.

It seems that the induced charge obtained from the simulation is always 4π times less than the numerical results in any case. I have tried different epsilon values, different charges and different distances. I have used lammps-23Jun2022 and lammps-8Feb2023 which is the latest version.And the in file is a little different due to the update of the lammps version.The analytical solution of polarization charge can be found in https://doi.org/10.1016/j.cpc.2019.03.006.

Is it caused by an error in my file or a bug?

Without you providing a suitable (simple!) input deck to reproduce this, there is little that can be done to investigate. Please use the most recent LAMMPS version. There were multiple bugfixes and consistency changes, so the package in the stable version can no longer be recommended for use.

Your question is rather specific to the DIELECTRIC package, hopefully @trungnth will see it and can provide some assistance. If there is no reaction within a few days, I suggest you file a proper bug report (including the reproducer input) as a GitHub issue instead.

Here is my in file.
variable epsilon1 index 5
variable epsilon2 index 80

variable data index suc.data

units real
atom_style dielectric
atom_modify map array
dimension 3
boundary f f f

variable method index icc # gmres = BEM/GMRES
# icc = BEM/ICC*
# dof = Direct optimization of the functional
# none

variable ed equal “v_epsilon2 - v_epsilon1”
variable em equal “(v_epsilon2 + v_epsilon1)/2”
variable epsilon equal “v_epsilon1”
variable area equal 21.65 # patch area, same as in the data file

pair_style lj/cut/coul/cut/dielectric 5.61231 100
read_data ${data}

group interface type 1
group cations type 2

group ions type 2 3

set type 2 x 0 y 0 z 24

variable qcations equal “2”
set group cations charge ${qcations}

pair_coeff * * 1.4988012 5
pair_coeff 1 1 0.0 5
mass * 35.450
pair_modify shift yes
write_data LJfinish.data
neighbor 5.0 bin
neigh_modify delay 0 every 1 check yes page 12000000 one 1000000

compute ef all efield/atom
compute v all stress/atom NULL
timestep 1.5

dump charge interface custom 100 charge.txt x y z q
dump charge1 interface custom 10000 charge1.txt x y z q
dump 2 all custom 1000 all.txt id type x y z q

fix freeze cations setforce 0 0 0
velocity cations set 0 0 0
compute charge interface property/atom q
fix move cations nve

if “{method} == gmres" then & "fix 3 interface polarize/bem/gmres 1 1.0e-6" & "fix_modify 3 itr_max 1000 dielectrics {ed} {em} {epsilon} {area} NULL kspace no" & elif "{method} == icc”&
“fix 3 interface polarize/bem/icc 1 1.0e-5 itr_max 50” &
“fix_modify 3 itr_max 1000 dielectrics {ed} {em} {epsilon} {area} NULL kspace no” &
elif “{method} == dof" & "fix 3 interface polarize/functional 1 0.001" & "fix_modify 3 dielectrics {ed} {em} {epsilon} {area} {qscale}” &
else &
"print ‘Unsupported polarization solver’ "

thermo 1
thermo_style custom step temp press evdwl ecoul elong epair f_3[1] f_3[2]
thermo_modify flush yes
fix atom interface ave/atom 1 8000 10000 c_charge
dump new interface custom 10000 avercharge.txt x y z f_atom
run 10000
write_data LJfinish.data

Sorry I can’t upload files because I am a new user. Is there an email address which I can send my files to it? Thank you for your timely reply!

You can always upload a zip or tar.gz archive with the files to a service like Dropbox, Google Drive, One Drive or equivalent and provide a link. However, even better would be to create a tiny input deck that has the minimum number of commands and a data file with as few atoms as possible to reproduce the issue, and then those can be included inline via cut-n-paste here.

(if you do, please place the quoted block into triple backquotes (```) so that the forum software will render it properly).

Thank you very much for your guidance. I have prepared the link.
https://drive.google.com/file/d/1hTE0B1Hz1Eub431AVujwqCT3AW1K2lLs/view?usp=sharing

Wrong tagging! I’m Trung Nguyen but not the author of DIELECTRIC package. His name is Nguyễn Đắc Trung (Trung Nguyen, email: [email protected] )

@Shibo_Pan Thanks for reporting the issue. I will look into it.