Core-Shell model and boundary setting


I’m doing a relaxation simulation for strontium titanate, and when I choose the boundary condition PPP, the relaxation works. But when I choose PPS as the boundary condition, the relaxation goes wrong and the bond atoms is missing. I tried to add a thick enough vacuum layer up and down the box. The boundary condition was still PPP, but the simulation was still wrong.

I would like to ask if there is any error in my code? or if the core-shell model of ferroelectric materials is only suitable for block-like simulations, and its polarization can cause instability when the boundary conditions change.

Here’s my code:

# ------------------------ INITIALIZATION ----------------------------
units 		metal
dimension	3
boundary	p	p	s
atom_style	full
timestep 0.0001
# ----------------------- ATOM DEFINITION ----------------------------
read_data STOtest3.lmp

group cores type 1 2 3
group shells type  4 5 6

neighbor 2.0 bin
comm_modify vel yes

# ------------------------ FORCE FIELDS ------------------------------
pair_style   born/coul/wolf/cs  0.25 6 8 # A, rho, sigma, C, D 
#                  A           rho       sigma   C           D 
pair_coeff   * *   0.0         1.000      0.00    0.00        0.00
pair_coeff   4 6   776.84    0.35867     0.00    0.0       0.00 #Sr-O
pair_coeff   5 6   877.20      0.38096     0.00    9.00      0.00 #Ti-O
pair_coeff   6 6   22764.3      0.1490     0.00    43.0        0.00 #O-O

bond_style class2
#                  R0     K2/2          K3       K4 
bond_coeff 1       0.0    5.703   0.0     0 #Sr
bond_coeff 2       0.0    32987.0   0.0      0 #Ti
bond_coeff 3       0.0    9.205   0.0      0 #O
special_bonds coul 0.0 0.0 0.0 

# ------------------------ Thermo output -------------------------------
compute CSequ all temp/cs cores shells
thermo 10
thermo_style custom step temp etotal press vol   
thermo_modify temp CSequ  

dump D1 all atom 1 stab.txt 

# ------------------------ Equilibration Run -------------------------------
min_style cg 
minimize 1e-8 1e-8 10000 10000

velocity all create 300 983629 dist gaussian mom yes rot no bias yes temp CSequ 
velocity all scale 300 temp CSequ 
fix 1 all nvt temp 300 300 0.001
fix_modify 1 temp CSequ

run 5000 
unfix 1 
undump D1 

A few things:

  • if you report errors, please always report what version of LAMMPS you are using and what platform you are running on
  • does the “Bond atom missing” error happen immediately or later during the simulation (and when, if later?)
  • is there any unexpected behavior of energies? How different are energies at the beginning between fully periodic and partially periodic?
  • you should check that there are no bonds across the periodic boundary you are turning into non-periodic
  • there should be no atoms with non-zero image flags in z-direction for the partially periodic case
  • are there any warnings that LAMMPS produces even with fully periodic boundaries

Thanks for your reply!

  1. I am using the Lammps version “22Jun 2022”. I am running on linux.
  2. “Bond atoms missing” error happen later during the simulation. It happen after 140 steps.
  3. Energy doesn’t change much at the begining, but it changes a lot when error happens. There is a difference of 200 eV around under different boundary conditions at the energy of PPP:-148274.78; total energy of PPS:-148029.17.
  4. There is no bonds across the periodic boundary.
  5. Sorry. I don’t understand what this means.
  6. there is a warning: WARNING: Using ‘neigh_modify every 1 delay 0 check yes’ setting during minimization (…/min.cpp:187). But I think it doesn’t matter.

Then you should set up a test input, that only runs for 200 steps and outputs energies and coordinates in every step.

This means that either your timestep is too large or that your force field parameters are not suitable. With periodic boundaries, atoms/molecules will have to mostly stay in place for solids, but when you remove a boundary, some reconstruction can happen and atoms that otherwise cannot be close, can get close(r). If then either the timestep is too large or the force field missing some repulsion, atoms can get too close, then forces get too large and atoms may be accelerated by a very large amount, which will then ultimately result in atoms moving so fast that LAMMPS cannot keep track of them, and they will be “lost”. For bonds, it is required that atoms forming a bond remain in the same subdomain or close by (as far as the communication cutoff goes).

With the previous suggestion or recording every step, you should be able to identify which atoms get close and have large forces and then check if this can be addressed by either reducing the time step or check whether there are missing or unsuitable parameters in your pair coefficients.

Thak you so much! I will give it a try based on your suggestion.