Charge neutrality of a BaTiO3 block after rotation

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps users,
I am trying to simulate Batio3 crystal rotated at various degrees: (15, 30, 45...). I could do the rotation of the crystal from original 100 crystal and then changing the simulation box size according to the size of the model. But after the rotation and changing the box size , when I tried to run an energy minimization it says : it cant compute MSM because the system is not charge neutral.

How can I make the system charge neutral? Do I need to delete some atoms that contribute to the charge imbalance or there is some other way around?

Here is the script:
# barium titanate using buckingham potential mode 1 fracture
dimension 3
boundary p m p
units metal
atom_style charge

neighbor 0.3 bin
neigh_modify delay 5

# create geometry

lattice custom 3.999 origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 &
basis 0.0 0.0 0.0 &
basis 0.5 0.5 0.5 &
basis 0.0 0.5 0.5 &
basis 0.5 0.0 0.5 &
basis 0.5 0.5 0.0 &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 1.00925231307827

region box block -20 20 -20 20 -0.5 0.5
region simbox block -20 20 -20 20 -0.5 0.5
create_box 3 simbox
create_atoms 1 region box &
basis 1 1 &
basis 2 2 &
basis 3 3 &
basis 4 3 &
basis 5 3

mass 1 137.327
mass 2 47.867
mass 3 15.9994

# rotating the simulation box and atoms
group inibox region box
displace_atoms inibox rotate 0 0 0 0 0 1 45
region outbox block -4 4 -4 4 -0.5 0.5 side out
group outbox region outbox
delete_atoms group outbox
change_box inibox x final -4 4 y final -4 4 z final -0.5 0.5 units lattice

#buckingham and coulombic potential
#pair_style buck/coul/msm 12
pair_style buck 5
pair_coeff 1 1 0.0 10.0 0.0
pair_coeff 1 2 0.0 10.0 0.0 #Ti-Ba
pair_coeff 1 3 1214.4 0.352200 8.0 #Ba-O
pair_coeff 2 2 0.0 10.0 0.0 #Ti-Ti
pair_coeff 2 3 877.2 0.38096 9.0 # Ti-O
pair_coeff 3 3 22764.0 0.1490 43.0 # O-O
#kspace_style msm 1.0e-4
# define groups

#region 1 block INF INF INF -6 INF INF
#group lower region 1
#region 2 block INF INF 4 INF INF INF
#group upper region 2
#group boundary union lower upper
#group mobile subtract all boundary

region crack block -0.5 0.5 -0.5 0.5 -0.5 0.5
group crack region crack
delete_atoms group crack

group Ba type 1
group Ti type 2
group O type 3

#set group O charge -2.0
#set group Ti charge 4.0
#set group Ba charge 2.0

compute stress all stress/atom
compute stress1 all reduce sum c_stress[1]
compute stress2 all reduce sum c_stress[2]
compute stress3 all reduce sum c_stress[3]
compute stress4 all reduce sum c_stress[4]
compute stress5 all reduce sum c_stress[5]
compute stress6 all reduce sum c_stress[6]
#dump eqli all xyz 10 batio3i.xyz
# initial velocities
#velocity lower set 0.0 -0.5 0.0
#velocity upper set 0.0 0.3 0.0

# fixes

fix 1 all nve
#fix 2 boundary setforce NULL 0.0 0.0

# run
timestep 0.003
thermo 10
#dump eql2 all custom 10 id type x y z batio3.file
#dump ID group-ID style N file args
dump eql2 all custom 10 batio3.* id type x y z
dump eql all xyz 10 batio3.xyz
dump 1 all cfg 500 stressimage_*.cfg mass type xs ys zs c_stress[2]
run 2000

Best Regards
Sakib

2

Dear Lammps users,
I am trying to simulate Batio3 crystal rotated at various degrees: (15, 30, 45...). I could do the rotation of the crystal from original 100 crystal and then changing the simulation box size according to the size of the model. But after the rotation and changing the box size , when I tried to run an energy minimization it says : it cant compute MSM because the system is not charge neutral.

How can I make the system charge neutral? Do I need to delete some atoms that contribute to the charge imbalance or there is some other way around?

you first need to understand where the charge imbalance is coming
from. is it because you have a polar surface or because your periodic
continuation is not correct at the box boundary. in both cases, it is
most likely that your geometry is not correct.

it is not so much about whether your system is neutral as recognized
by LAMMPS, but whether your geometry *should* be neutral. that is
something *you* have to know.

axel.

3

BATiO3 is a charged system The unit cell (100 direction) has 3 Oxygen atoms each at the face, 1 Ba atom at the corner(0, 0 , 0) and 1 Ti atom at the middle (0.5, 0.5, 0.5).
Ba has 2+ charge,
Ti has 4+ charge,
O has 2- charge.

After the rotation (110 crystal, 45 degree rotation w.r.t z axis ), the unit cell now has 2 O, 2 Ba and 1 Ti atoms. So the unit cell itself is not charge neutral.

Another thing is , if I want to have a non periodic condition in y direction then how would I implement charge neutrality? Are periodic systems always charge neutral?

Regards
Sakib

4

BATiO3 is a charged system The unit cell (100 direction) has 3 Oxygen atoms each at the face, 1 Ba atom at the corner(0, 0 , 0) and 1 Ti atom at the middle (0.5, 0.5, 0.5).
Ba has 2+ charge,
Ti has 4+ charge,
O has 2- charge.

After the rotation (110 crystal, 45 degree rotation w.r.t z axis ), the unit cell now has 2 O, 2 Ba and 1 Ti atoms. So the unit cell itself is not charge neutral.

a unit cell doesn't change because of the rotation. the problem is
that you define your region boundaries exactly on the spot of lattice
atom positions. this is a common beginner's mistake and has been
discussed on the mailing list many times. please re-read the
corresponding parts of the manual and search through mailing list
archives for those previous discussions. also, there are discussions
on the impact of rotating the lattice on the requirements on periodic
boundaries.

Another thing is , if I want to have a non periodic condition in y direction then how would I implement charge neutrality? Are periodic systems always charge neutral?

as i already stated before, *you* have to know what these surfaces
should look like and choose your region boundaries accordingly.
surfaces can be neutral or polar depending on how you cut them from
the bulk system and what your unit cell looks like. what you should
get is your choice and depends on the science of your problem. that is
not really a LAMMPS problem.

axel.