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.
Thanks,
Jorge
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