# Code modified from lammps example of core-shell model of NaCl
# PbTiO3 core-shell model
# All Pb, Ti, and O atoms are treated as core-shells
# Pb-O, Ti-O , Pb-Ti, and O-O interactions are treated with buckingham and coulombic potentials
# the parameters are borrowed from 
# "Atomic-level simulation of ferroelectricity in oxide materials" by 
# Sepliarsky et al, current opinion in solid state and mat. sci, (2005) 9, 107-113

# Initialize Simulation
clear
units metal
boundary p p p
atom_style full

# Create Atoms

#fix csinfo all property/atom i_CSID
read_data PTO.data #fix csinfo NULL CS-Info
group cores type 1 2 3
group shells type 4 5 6

neighbor 2.0 bin
comm_modify vel yes

# Define Interatomic Potential
kspace_style ewald 1.0e-6
pair_style buck/coul/long/cs 10.0
pair_coeff * * 0.0 1.000 0.00 0.00
pair_coeff 4 5 0.09  2.4201   12.57  # Pb-Ti
pair_coeff 4 6 6766.27  0.2738   127.78  # Pb-O
pair_coeff 5 6 1130.01 0.3597   -160.84  # Ti-O
pair_coeff 6 6 3634.86  0.3144  331.60   # O-O
pair_modify table 0

bond_style class2
#            r0  k2     k3  k4
bond_coeff 1 0   119.5  0   17968.0  # Pb core-shell
bond_coeff 2 0   1428.0 0   0	     # Ti core-shell
bond_coeff 3 0   23.3   0   4514.0   # O  core-shell
#bond_coeff 1 0   59.75  0   748.67  # Pb core-shell
#bond_coeff 2 0   714.0  0   0	     # Ti core-shell
#bond_coeff 3 0   11.65  0   188.083   # O  core-shell

# Equilibration Run
reset_timestep 0
timestep 0.0004 # picosecond
thermo 100
thermo_style custom step etotal pe ke ebond temp press vol

compute CSequ all temp/cs cores shells
compute CSpress all pressure CSequ

compute thermo_press_lmp all pressure thermo_temp # press for correct kinetic scalar
variable volume equal vol

thermo_modify temp CSequ press thermo_press_lmp

compute cs_chunk all chunk/atom molecule
compute cstherm all temp/chunk cs_chunk temp internal com yes cdof 3.0
variable temp_int equal ave(c_cstherm)

dump 1 all atom 300 dump.relax
dump_modify 1 scale no sort id

velocity all create 300 13487 dist gaussian mom yes rot no bias yes temp CSequ
velocity all scale 300 temp CSequ
fix atpos all ave/time 1 1000 1000 c_CSequ c_CSpress c_thermo_press_lmp v_volume v_temp_int file avg_PT.txt
fix a all npt temp 300 300 0.04 iso 0 0 0.4
fix_modify a temp CSequ press thermo_press_lmp # press for correct kinetic scalar
timestep 0.0004
run 20000
unfix a
fix b all npt temp 300 300 0.04 aniso 0 0 0.4
fix_modify b temp CSequ press thermo_press_lmp # press for correct kinetic scalar
run 20000
unfix b
fix c all npt temp 300 300 0.04 tri 0 0 0.4
fix_modify c temp CSequ press thermo_press_lmp # press for correct kinetic scalar
run 20000