Dear lammps users,
I need the radial distribution function of the pairs of atoms in my system at various temperature, to study how the atoms distribution changes with temperature. And I am doing this for spherical particle of sizes varying from 2-10 nm.
Strangely, with increasing temperature, the first peak of rdf shifts leftwards up to a distance of 1 Angstrom. And also the peak intensity is around 3000 a.u for high temperature (3000K). Is this a reliable output?
I have given my code below for better clarification of what I am doing.
# ---------- Initialize Simulation ---------------------
#Type 1-Ba, Type 2-Ta, Type 3-O, Type 4-N
clear
log 94.log
units metal
dimension 3
boundary p p p
atom_style charge
lattice custom 3.866 &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0 &
basis 0.5 0.5 0.5 &
basis 0.5 0.5 0.0 &
basis 0.0 0.5 0.5 &
basis 0.5 0.0 0.5
region whole block -20 20 -20 20 -20 20 units lattice
create_box 4 whole
#-------Defining two spheres-----
lattice custom 3.866 &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0 &
basis 0.5 0.5 0.5 &
basis 0.5 0.5 0.0 &
basis 0.0 0.5 0.5 &
basis 0.5 0.0 0.5
region 1 sphere 7 0 0 2.6 units lattice #creating sphere of radius (2.6*3.866)=10 A (around 1nm radius)
region inner1 sphere 7 0 0 1.3 units lattice
create_atoms 4 region 1 basis 1 1 basis 2 2 &
basis 3 3 basis 4 3 basis 5 4
lattice custom 3.866 &
a1 1.0 0.0 0.0 &
a2 0.0 1.0 0.0 &
a3 0.0 0.0 1.0 &
basis 0.0 0.0 0.0 &
basis 0.5 0.5 0.5 &
basis 0.5 0.5 0.0 &
basis 0.0 0.5 0.5 &
basis 0.5 0.0 0.5
region 2 sphere -5 0 0 2.6 units lattice
#— To group the surface atoms alone —
region inner2 sphere -5 0 0 1.3 units lattice
region surf1 intersect 2 1 inner1 side out
region surf2 intersect 2 2 inner2 side out
create_atoms 4 region 2 basis 1 1 basis 2 2 &
basis 3 3 basis 4 3 basis 5 4
mass 1 137.327
set type 1 charge 2.1
mass 2 180.95
set type 2 charge 3.326
mass 3 16.00
set type 3 charge -1.48
mass 4 14.00
set type 4 charge -2.576
group sphere1 region 1
group sphere2 region 2
group Ba type 1
group Ta type 2
group O type 3
group N type 4
# ---------- Define Interatomic Potential ---------------------
pair_style hybrid morse 10 coul/wolf 0.7 10 coul/wolf 0.7 10 coul/wolf 0.7 10 coul/wolf 0.7 10 coul/wolf 0.7 10 coul/wolf 0.7 10
pair_coeff 1 3 morse 0.5799 2.232 2.737 #Ba-O
pair_coeff 1 4 morse 0.858 2.33 2.677 #Ba-N
pair_coeff 2 3 morse 1.797 2.087 1.839 #Ta-O
pair_coeff 2 4 morse 1.7369 2.16 1.987 #Ta-N
pair_coeff 1 1 coul/wolf 1
pair_coeff 1 2 coul/wolf 2
pair_coeff 2 2 coul/wolf 3
pair_coeff 3 3 coul/wolf 4
pair_coeff 4 4 coul/wolf 5
pair_coeff 3 4 coul/wolf 6
neighbor 2.0 bin
neigh_modify delay 10 check yes
#----defining certain variables
variable time equal step
variable temp equal temp
variable pressure equal press
variable potential equal pe
variable kinetic equal ke
variable volume equal vol
variable energy equal etotal
variable H equal enthalpy
variable den equal density
compute 1 all rdf 50 1 1 1 2 1 3 1 4 2 2 2 3 2 4 3 3 3 4 4 4 cutoff 10
compute coord all coord/atom cutoff 10 *
#-----Equlibration—
thermo 1000
thermo_style custom step temp pe press ke vol etotal enthalpy density
velocity all create 300 20000 mom yes rot yes
fix 1 all npt temp 300 300 0.001 iso 1.0 1.0 1000.0
variable temp equal temp
dump 1 all custom 5000 dump.equ1.* id type x y z c_coord
fix rdf all ave/time 1000 5 50000 c_1[*] file gr.dat mode vector
fix 2 all ave/time 1000 50 50000 v_temp v_pressure v_energy v_H v_den v_volume ave one file avgtemp
fix 3 all ave/time 50000 1 50000 v_energy ave one file energy
run 50000
unfix 1
# unfix rdf
displace_atoms sphere2 move 6 0 0
#-----Sintering to 400K
fix 1 all npt temp 400 400 0.01 iso 1.0 1.0 1000.0
variable temp equal temp
run 50000
unfix 1
thermo 1000
thermo_style custom step pe temp press ke
#-----Sintering to 500K
fix 1 all npt temp 500 500 0.01 iso 1.0 1.0 1000.0
variable temp1 equal temp
run 50000
unfix 1
#-----Sintering to 600K
fix 1 all npt temp 600 600 0.01 iso 1.0 1.0 1000.0
variable temp1 equal temp
run 50000
unfix 1
#-----Sintering to 700K
fix 1 all npt temp 700 700 0.01 iso 1.0 1.0 1000.0
variable temp1 equal temp
run 50000
unfix 1
#-----Sintering to 800K
fix 1 all npt temp 800 800 0.01 iso 1.0 1.0 1000.0
variable temp1 equal temp
run 50000
unfix 1
#-----Sintering to 900K
fix 1 all npt temp 900 900 0.01 iso 1.0 1.0 1000.0
variable temp1 equal temp
run 50000
unfix 1
#-----Sintering to 1000K
fix 1 all npt temp 1000 1000 0.01 iso 1.0 1.0 1000.0
variable temp1 equal temp
run 50000
unfix 1
Thanks in advance.