# Compute rdf - higher temperature

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.

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.

the standard radial distribution functions are not a good measure for such a system. they are meant to be used for bulk systems with a homogeneous distribution of particles (you are effectively comparing the actual relative distribution of all pairs of particles against that of an ideal gas). applying g® to non-homogeneous particle distributions makes the interpretation of the results “tricky”.

that said, your quoted input is very long and convoluted. there is not much of a chance, that somebody will debug it for you.
one thing stands out: why so many different coul/wolf substyles? that will likely lead to bogus results. why not use one sub-style and pair_coeff * * coul/wolf ?

axel.

Thank you for the feedback Axel. Though my system is non-homogeneous, I am calculating rdf for each pair of atoms only, like between 1-1, 1-2 like that. Would it still be not good?

Thank you for the feedback Axel. Though my system is non-homogeneous, I am calculating rdf for each pair of atoms only, like between 1-1, 1-2 like that. Would it still be not good?

this is a redundant question. OF COURSE, you compute a rdf only between pairs of atoms. thus i have already answered your question.
so please re-read my previous reply and think about what i have explained.
…and if you don’t trust me, perhaps read up in a statistical mechanics book about what a g® function means.

axel.