Issue with the calculation of a radial distribution function

LAMMPS
1

#############################

P R E L I M I N A R I E S

#############################

variable Tmp equal 0.68

units lj # Use reduced LJ style units
atom_style atomic # Uncharged point particles
atom_modify map hash # Needed for indexing of clusters etc

lattice fcc 1.0 spacing 1 1 1
region box block 0 4 0 4 0 4 units lattice
create_box 1 box
create_atoms 1 box

variable ra uloop 3

variable rand equal (25373+${ra})

log cryst_t${Tmp}_r${rand}.log

Set the mass of the first (and only) atom type.

mass 1 1.0

Lennard-Jones interactions between particles, 3.5 sigma cut-off. Apply

long range tail corrections to energy and pressure

pair_style lj/cut 3.5
pair_modify tail yes
pair_coeff 1 1 1.0 1.0

#############################

M E L T S Y S T E M

#############################
velocity all create 2.4 ${rand} mom yes dist gaussian # Assign velocities

timestep 0.002 # simulation timestep
thermo 100 # output thermodynamic data every 100 steps

trajectory first 13000 steps

dump traj_first all custom 100 traj_first_${ra}.dat id type x y z vx vy vz

###############
#first fix
fix 1 all npt temp 2.4 2.4 0.1 iso 5.0 5.0 0.5 mtk yes tchain 5 pchain 5
run 10000

###############
###############
#second fix
fix 1 all npt temp 1.25 1.25 0.1 iso 5.68 5.68 0.5 mtk yes tchain 5 pchain 5
run 1000

################

################
#third fix
fix 1 all npt temp 1.25 ${Tmp} 0.1 iso 5.68 5.68 0.5 mtk yes tchain 5 pchain 5
run 2000

#g(r) third fix
compute myRDF1 all rdf 150
fix rdf_calc1 all ave/time 50 10 500 c_myRDF1[*] file rdf_first_${ra}.dat mode vector
run 0
###############

undump traj_first

#############################

F R E E Z E S Y S T E M

#############################

Define solid atoms. This closely follows ten Wolde Ruiz-Montero and Frenkel Faraday Discuss 1996 104 93-110

Compute components of the per-atom q6 vector

compute q6 all orientorder/atom degrees 1 6 components 6 nnn NULL cutoff 1.3

get number of connections

compute coord_number all coord/atom orientorder q6 0.5

An atom is solid if it has 8 or more connections

variable is_solid atom c_coord_number>=8
group solid dynamic all var is_solid # Must be dynamic to update

do clustering

compute cluster solid cluster/atom 1.3

define chunks one chunk per cluster

compute clus_chunks solid chunk/atom c_cluster

count the size of each chunk

compute size_chunks solid property/chunk clus_chunks count

Find the maximum entry in the vector of chunk sizes

variable max_n equal max(c_size_chunks)

Thermo style which includes this variable

thermo_style custom step temp pe vol v_max_n

thermo 100 # Print the thermo information every 100 steps

Reset the npt fix at a lower temperature (below freezing)

fix 1 all npt temp {Tmp} {Tmp} 0.1 iso 5.0 5.0 0.5 mtk yes tchain 5 pchain 5

print positions and velocities

dump traj all custom 100 traj_${ra}.dat id type x y z vx vy vz

run 100000 # Run for this many steps

g(r) final state

compute myRDF2 all rdf 1024

fix rdf_calc2 all ave/time 1 1 1 c_myRDF2[*] file rdf_final_${ra}.dat mode vector

run 0

clear
next ra
jump SELF

Hello everyone, could someone kindly help me?

This is my code for studying nucleation. In the code, the RDF is calculated twice. However, I can’t understand why the second calculation produces an empty file (with no data) when I use numbers other than 1 1 1 in the command

fix rdf_calc2 all ave/time 1 1 1 c_myRDF2[*] file rdf_final_${ra}.dat mode vector

For example, using 50 10 500, as I did for the first RDF calculation, doesn’t work. Does anyone know the reason for this behavior?

P.S. The code is used to study the nucleation of the system, which starts at T=2.4 in LJ units and is gradually cooled down to T=0.68. The goal is to calculate the q6 parameter and observe how many atoms crystallize during the evolution.

Thank you!

2

This should be because you have a run 0 after the command, so basically there is just one frame. Is run 0 really what you intended to write after that second compute rdf you are defining?

PS: also, note that summoning compute rdf commands after some run xxx command has been summoned can be troubling if you dont handle the Nfreq keyword in the fix ave/time correctly: the fix ave/time will output things only at timesteps that are a multiple of Nfreq, so using reset_timestep (see reset_timestep command — LAMMPS documentation) might be a good idea in general.

3

I used run 0 because I want to calculate the RDF immediately after 13,000 steps and again after the final 100,000 steps. My goal is to observe a liquid state in the first case and a solid state in the second. I have conducted several tests, and only with run 0 did I get any results.

4

Hi @Gian_Maria_Lombardi,

As @ceciliaalvares hints, your issue is because of how fix ave/time works and not compute rdf.

Data are written to the file only when N_{repeat} values have been accumulated, this includes the initial computation. When using 1 1 1 as an input, calling the fix during the system setup trivially computes 1 value, thus writes the file. However, no averaging was performed. I suggest you go through the fix ave/time manual page where this is detailed.

Also have a look at the forum guidelines concerning the formatting of your posts. You’ll also find general advice to troubleshoot common issues and what to provide for better help.

5

Well, then you are achieving your goals successfully in the second compute rdf command as the 1 1 1 combination is the way to go. It should also be this combination for the first compute rdf command (additionally, in the former case, you may want to unfix the first fix ave/time you define (unfixing to be done after a run 0 command) in order not to have a huge output in the rdf file associated with the command). Note that the combination you used in your first compute rdf command is not giving you what you want given the combination of Nevery, Nfreq, Nrepeat you chose in the input file shown above.

As Germain said, check the fix ave/time page of the manual carefully, specifically the part where Nevery, Nfreq and Nrepeat are defined, and you will solve your problem :slight_smile: