Problem calculating Madelung constant

Hi all,

I am following one of the user-contributed scripts to calculate the Madelung constant of a NaCl crystal. Just to check my understanding of the lattice, region, create_box and create_atoms commands, I am trying to replicate exactly what is done in the mentioned script but, instead of reading the coordinates of the atoms from a file, I am generating the crystal using Lammps commands. The strange thing is that I am able to get exactly the same coordinates as in the file, but the the energy and value of the Madelung constant that I obtain are wrong. Can anybody see why? I am attaching my script below. Any help will be greatly appreciated.



PS: I’m using the keyword “origin” in the lattice command so that the atoms are exactly in the same position as those in the file “data.Madelung”

units lj
atom_style full
dimension 3
neighbor 0.5 bin
neigh_modify every 1 delay 3 check yes
boundary p p p

dielectric 1.0

pair_style lj/cut/coul/long 2.5 10.0

Position of ions

lattice fcc 6.25e-2 origin 0.25 0.25 0.25
region region1 block -2.5 2.5 -2.5 2.5 -2.5 2.5 units lattice
create_box 2 region1
create_atoms 1 box
mass 1 1
lattice fcc 6.25e-2 origin 0.75 0.25 0.25
create_atoms 2 box
mass 2 1
set atom 1 charge -1
set atom 2 charge 1

velocity all set 0.0 0.0 0.0 sum no units box

pair_coeff * * 0.0 1.0 2.5 # zero pairwise interactions

variable el equal ecoul+elong
variable Mc equal (ecoul+elong)*-4 # Madelung constant

thermo 1
thermo_style custom step pe v_el v_Mc
thermo_modify flush yes norm yes

kspace_style pppm 0.00001
run 0

Hi Jorge,

Looks like you need to set the charges for all the atoms,

set type 1 charge 1.0
set type 2 charge -1.0

instead of setting for only two atoms

set atom 1 charge 1.0
set atom 2 charge -1.0