md_lj_tree - Hierarchical (tree force calculation) simulation of Lennard-Jones liquid

DOWNLOAD

Click the link md_lj_tree to download the code.

ANIMATIONS AND TESTS

In the right panel we show some tests and animations.

SYNOPSIS

md_lj_tree [parameter_file_name] [options]

DESCRIPTION

md_lj_tree - Simulates the evolution of a N-body system interacting with a Lennard-Jones potential. The force is computed by walking a hierarchical tree. Code complexity is O(N Log N). Periodic boundary condition is used and temperature is kept constant.

OPTIONS

The options have the structure

option_name = <option_value>

Options and their possible values are:

• paramfile - in this option you give the name file with the values of the input parameters. Overwrites parameters values below. You may also input this filename by only writing:

md_lj_tree parameters_file_name

Parameter input file may be created by hand with the editor of your choice. Comment lines start with an "%". Follow each name option with a blank space and the option value. The order of the option lines does not matter.

Also you may create an example input file by executing

md_lj_tree

This will run the md_lj_tree code with default values and when it finish you will have in your running directory the file "parameters_null-usedvalues". Now you may edit this file to adapt to your own simulation parameters. It may be helpful to change this file name to whatever apropriate.

• forcecalc_method - is the force calculation method to use. We have implemented four methods: direct (direct), nearest neighborh list (nblist), cells method (cells) and tree walk method (normal) which is the default value.

• theta - is a force acuracy parameter. Not in use by now.

• usequad - is the option to include quadrupoles in the force calculation. Not in use by now.

• temperature - is the temperature of the simulation in units such that eps/kB is the unit of temperature (see below section on "UNITS").

• density - is the density of the liquid in code units.

• stepEquil - is the step to begin equilibrium computations. Its default value is 100 time steps.

• stepVel - is the number of steps jumps to save a velocity histogram.

• stepAvgVel - is the number of velocity histograms to average. Its default value is 4.

• sizeHistVel - is the array size for the velocity histogram. Its default value is 50.

• rangeVel - is the range of velocities for the histogram computation.

• stepRdf - is the number of steps jumps to save a radial distribution function (RDF) histogram.

• stepAvgRdf - is the number of RDF histograms to average. Its default value is 20.

• sizeHistRdf - is the array size for the RDF histogram. Its default value is 200.

• rangeRdf - is the range of positions for the RDF histogram computation.

• nbody - it is the number of bodies to simulate. If icfile option is null the code will generate a initial condition (a test run) were all the particles are distributed uniformly in a cubic box with gaussian random velocities. Therefore, nbody will be of the form n^3.

• dtime - is the time step integration in code units. Can also be given in the form of p/q, where q is expressed as powers of 2 so we can think in terms of integration frequency.

• tstop - is the time to stop the simulation.

• seed - it is the Random number seed for the test run.

• icfile - you give here the name of the file with the N-body initial data.

• icfilefmt - is the format of the initial condition file: 'snap-bin' (binary) or 'snap' (ASCII) or 'snap-pv'. See snapoutfmt option below for a description of the file formats.

• snapout - you give here the name structure for the output of N-body snaps. The format follows as the ones used in C-language for integers ("%0#d").

• snapoutfmt - you tell the code the format of the snaps output. There are three options: the standard ASCII snap n-body format (snap-ascii); the binary standard snap format (snap-bin); and the columns format (snap-pv). The n-body standard format binary or ASCII is a file with n-body data written as follows:

nbody

NDIM

tnow

x y z (position for all the particles)

vx vy vz (velocity for all the particles)

Id (for all the particles)

And in the colummns format the particle data is written in column form as

# nbody NDIM time

# nbody value NDIM value time value

Id mass x y z vx vy vz (for all the particles)

You may have optionally written the particle's potential and acceleration vector if in the option you write 'out-phi' and/or 'out-acc'.

• dtout - this is the time for output a snap file. The out files will be written every dtout time. Can also be given in the form of p/q, where q is expressed as powers of 2.

• dtoutinfo - this is the time for output in the stadar output (stdout). The output to the stdout will be written every dtoutinfo time. Can also be given in the form of p/q, where q is expressed as powers of 2. This option is useful for controlling the cpu time consumed by output processing.

• statefile - you give here the name of a file where the run state will be saved. If it is null no run state will be saved.

• restorefile - if it is not null a run will be restarted from the data stored in this file.

• options - you may give here various code behavior options. They are, "reset-time" (inputdata); "out-phi" (outputdata); "out-acc" (outputdata). If you save a state file, this parameter is saved.

UNITS

Units are such that eps=sigma=m=1, where eps is the potential depth of the Lennard-Jones potential, sigma the cut radius, and m is the mass of each particle. We also have that the Boltzmann constant kB=1 such that eps/kB is the unit of temperature.

OUTPUT AND THERMODYNAMICS

md_lj_tree code produces several output files by default: md_lj_tree.log, thermo.dat, snap.dat, vel.dat, rdf.dat, and if you instruct other additional files are produced, such as those instructed by snapout options. The file md_lj_tree.log is the log file were some simulation parameters are save such as size of the simulation box, number of particle created if a initial condition file were not used, and so on. The files snap.dat, vel.dat, and rdf.dat are intended to be used in combination with gnuplot. In files vel.dat and rdf.dat are saved the velocity and radial distribution function histograms, they are produced in steps greater than stepEquil, and its frequency is controled by the options stepVel, stepAvgVel, stepRdf, and stepAvgRdf. In snap.dat file a snapshot of the system is saved and this file is produced every dtout/dtime steps. Thermodynamics parameters such as pressure and its fluctuation are saved in file thermo.dat every dtout/dtime steps.

SPATIAL DIMENSIONS

md_lj_tree code may be run in two or three spatial dimensions. To choose two dimensions edit the file "vectdefs.h" in directory "NagBody_pkg/General_libs" and choose TWODIM. Recompile again the code.

STOPPING THE CODE

Once md_lj_tree is running you may allways stop it by executing the command

echo > stop

You must be in the same directory were the process were lunched. A file named 'stop-state' is saved containing the state of the simulation run. This option permits us to stop the simulation and change numerical and/or physical parameters such as temperature or density.

EXAMPLES

Executing

md_lj_tree

will run the default simulation which is consistent with the experimental data reported in: J.L. Yarnell, et al, "Structure factor and radial distribution function for liquid Argon at 85 K", Phys. Rev. A7 (1972) 2130. Parameter eps in the Lennard-Jones potential that fit the experimental data is eps/kB = 119 K.

With the command

md_lj_tree nbody=4096 out=snap%03d dtout=2/256

we run a simulation with 4096 bodies and save a snap data every other time step.

With the command

md_lj_tree statefile=state

we run a simulation were a state file is saved every dtout/dtime steps. If for some reasons the simulation run is interrupted you may always restart the run ejecuting the command

md_lj_tree restorefile=state

Or if you stop it the run using the command 'echo > stop' you may restart the run with the command

md_lj_tree restorefile=stop-state

ANIMATIONS

You may run gnuplot to see animation plots. Make the following script:

a=a+1

plot "rdf.dat" w l

pause 1

if(a<50000) reread

Save the script in a file named "rdf.gnu". At the same time that the simulation is running start gnuplot and execute the commands:

a=0

load "rdf.gnu"

Also you may see animation of the particles positions. Make the following script:

a=a+1

plot "snap.dat" w p pointtype 6

pause 1

if(a<50000) reread

Save the script in a file named "snap.gnu". Again, at the same time that the simulation is running start gnuplot and execute the command:

a=0

set xrange [-5:5]

set yrange [-5:5]

load "snap.gnu"

The values of "xrange" and "yrange" options may change according the particle positions ranges.