Simulations

Introduction

Assorted topics relevant to programming particle-based simulation codes and to using these codes for the modeling of a wide range of systems, notably complex fluids.

Molecular dynamics simulations

By and large the CSML uses LAMMPS for MD, though depending on the application NAMD, GROMACS, or any of a host of other packages may be useful. Most issues can be resolved by consulting the LAMMPS manual, though some common problems are addressed below.

LAMMPS Special Usage Notes

Temperature Normalization

By default LAMMPS normalizes the temperature by an amount $n_{\text{dof}}-d$, where $n_{\text{dof}}$ is the system's total number of degrees of freedom and $d$ is the system's dimensionality. Subtracting $d$ accounts for the center-of-mass motion of the system. This leads to an incorrect reported value if the system has a proper frame of reference, e.g., when using a Langevin thermostat in which all particles interact with a stationary background solvent. In this case it is necessary to ensure $n_{\text{dof}}$ is used instead of $n_{\text{dof}}-d$. To do this, use compute_modify as follows

compute myTemp all temp
compute_modify myTemp extra 0
thermo_modify temp myTemp

As a note, the above only affects the reported temperature. The dynamics are computed correctly regardless.

Compute RDF using 'rerun' command

The 'rerun' command in LAMMPS performs a post-processing simulation by reading the atom information line-by-line from the dump file(s) created from a previous simulation. The command syntax is as follows:

 rerun file1 file2 ... keyword args ... 

A detailed description of the syntax can be found on the LAMMPS website.

Besides the fact that the atoms' positions (and possibly velocities, etc.) are pre-determined from the dump file(s), we use the rerun command as if we are running a normal simulation (with some differences and limitations, explained below). When the rerun command is called, it invokes the read_dump command to read in lines from the dumpfile(s) line-by-line, each time invoking the run command to output computed energy, forces, and any thermo output or diagnostic info the user has defined. Thus, in the input file for this pseudo simulation, we must define a system, units, dimensions, box, etc, and these will typically be identical to the original simulation.

Commands from the original simulation that will not be included are ones such as dump commands and time integration fixes (e.g. fix nve; rerun only looks at single moments in time and cannot perform time integration). Fixes that constrain forces on atoms (such as fix langevin) can be invoked in general, but it does not make sense to do this for computing the RDF (even though the langevin thermostat may be employed in the original simulation).

As an example, let us consider computing the RDF for a typical one-component shifted truncated Lennard-Jones fluid past the interaction cutoff, and let us assume that we have already generated a dumpfile containing information on the atom positions over some set of timesteps. Then we will run a second simulation that reads in the particle positions from the dumpfile(s) over some subset of the original recorded timesteps (see arguments for the rerun command), and will compute and output the RDF with the following commands:

compute    rdfID groupID rdf N #computes rdf with N bins
fix        fixID groupID ave/time Nevery Nrepeat Nfreq c_[rdfID] file rdf.dat mode vector # see note below
rerun      dump.dat dump x y z # 'dump.dat' is the dumpfile to be read 

Note: Since the compute rdf command will only compute the RDF up to the interaction cutoff distance, we must change this parameter in the pair_style and pair_coeff commands so that we can obtain the RDF over the desired domain (i.e. if we want to compute the RDF up to a cutoff of 4.0, we would set the 'cutoff' arguments in those commands to 4.0).