Comparing the output from MD and MC simulations
Introduction
In this module you will compare the output from molecular dynamics simulations (MD) and Monte Carlo simulations (MC) on liquid argon. The molecular dynamics simulations will be carried out using the DISCOVER MD program. The Monte Carlo simulations will be carried out using the program NVT_LJ that applies specifically to spherical molecules (or atoms) written by Dr. Pikunic. We will start with a box of argon atoms in a face-centered cubic lattice. The first issue that arises from this initial configuration is to find the appropriate randomized configuration of argon atoms. The second issue that arises is to equilibrate the argon atoms at the temperature of the simulation. Then we will propagate the argon in a number of steps. In MD simulations these are explicit time steps. In MC simulations these are steps in configuration space. Once the simulations are completed we will compare the results of the two methods in terms of their accuracy and advantages they each possess for determination of relevant physical properties. Please note that the DISCOVER force field assumes distance in units of Å and energy in units of kcal/mol. The MC program (NVT_LJ) is in reduced units.
Reduced units
The code that you will use in this exercise performs Monte Carlo simulations of a Lennard-Jones fluid in the Canonical ensemble (constant NVT). All the variables (in the simulation code, input and output files) are reduced using the Lennard-Jones parameters:
Lengths
Density
Temperature
Energy
Pressure
To compare the output from such an MC simulation you will need to dimension the output accordingly. This requires an analysis of the units in the MD simulations. In the CVFF force field used in DISCOVER and other related force fields used for biological applications (e.g. CHARMM and AMBER), the Lennard-Jones potential is expressed as:
In the MD simulations the potential energy u(r) has units of kcal/mol and distances are in units of Å. Thus, the A and B parameters have units of Å12kcal/mol and Å6kcal/mol, respectively. An alternative expression is:
Thus, we can see that A = 4es12 and B = 4es6. The often used values for argon are e/k = 120 K and s = 3.4 Å. This implies that e = 1.65 x 10-21 J per molecule or e = 0.997 kJ/mol. Given that 4.184 J = 1 cal, we have e = 0.238 kcal/mol. For Argon the parameters for the Lennard-Jones potential are A = 298302 and B = 193. These values are fixed in the CVFF force field that you will use for the MD simulation.
Since the MC simulation is in reduced units, you will have to convert the output energy, pressure, fluctuations, and radial distribution into units of kcal/mol and Å in order to compare the output from the two simulations.
Converting input files
The program DISCOVER uses two files ([name].car and [name].mdf) for the coordinates and atom type. As DISCOVER runs it outputs the configurations to a history file ([name].his) that is later used for analysis of the molecular dynamics trajectory. Since we will follow the calculation originally carried out by Rahman using 864 argon atoms in a cubic box with a dimension of 10s (where s is the LJ parameter), the input files you will need are:
Analysis of the results in DECIPHER is performed after the trajectory is complete.
The program NVT_LJ uses a file called initial.conf to input the initial configuration. The input parameters are read from inputfile. The output energy (U) and pressure (P) are written to a file called blkaverages. The output radial distribution function is written to a file called gr
The computer program MC_to_MD takes a file of the format of initial.conf and writes a DISCOVER input coordinate file.
The computer program MD_to_MC takes a file of the format of [name].car and writes a NVT_LJ input coordinate file.
Note that you can read a file with the format of initial.conf directly into InsightII provided you use the free_format option under the Get command and specify a format file mc_file.frc. You must specify the path to this file as well.