Mod:Hunt Research Group/MoltenSaltSimulation
Setting up and running Molten Salt Molecular Dynamic Simulations
There are several pieces of information needed before staring any MD simulation
- The density of the system at the final temperature as it is used to determine the size of the simulation box
- Parameters for the Field File
- Potential to use in the Field File
Once the information is found, the input structure can be generated.
MD simulations on Molten Salts should be ran from a crystal structure. I don't know where to find crystal structures
Using Packmol to generate starting configurations
The density of the salt at the final temperature is needed to determine the box size for the simulation. The size of the simulation box can be determined by employing the following equation.
Volume = molar mass / 0.6022 (density / number of molecules)
Taking the cube root of the volume gives the length of the box size.
For each molecule considered in the MD simulation an input structure needs to be generated. These files need to have a the name ‘FILENAME.xyz’. For chlorine the .xyz file looks like this
1 cl CL 0.00000 0.00000 0.00000
Once the input structure files are generated (one for each ion/molecule in the system) another file needs to be created which is the generation file that packmol uses. The generation file needs the size of the box, so it should have already been determined. The packmol generation file has to be named ‘FILENAME.inp’.
# # Generation of NACL # tolerance 3.0 filetype xyz output nacl_200ips.xyz structure na.xyz number 200 inside box -10. -10. -10. 10. 10. 10. end structure structure cl.xyz number 200 inside box -10. -10. -10. 10. 10. 10. end structure
The tolerance is....
Output nacl_200ips.xyz is the name of the file that is generated by packmol.
Structure na.xyz reads the molecular structure in the file na.xyz
Number 200 packs 200 Na atoms (can be moleculaes) into the box.
Inside box -10. -10. -10. 10. 10. 10. is the size of the box that is generated. On a grid the axes are drawn from -10 to 10, meaning the length of each side of the box is 20 Å.
End structure tells packmol that is everything that needs to be known for that particular molecule. The same process the occurs for the Cl ions.
Running Packmol from the Terminal
This is a relatively simple step! Start in the same directory as the .inp and the .xyz files. From here navigate to the directory where packmol is stored. On my Mac it is in Applications.
../../../../Applications/packmol/packmol <FILENAME.inp
doing this will run packmol and generate a file with the name specified in the .inp file (nacl_200ips.xyz).
Its now necessary to move the newly generated .xyz file to the hpc, this can be done using the scope command in the terminal. Each MD simulation that is run needs to have its own directry, as all the files are called the same thing.
Crystal Structures
Making the CONFIG file
Converting Packmol to CONFIG
Once the file is on the HPC it is necessary to get the file into the correct format for dlpoly to read. To change the file format a perl conversion script, conv_xyz2conf.pl, needs to be used.
Conversion script is shown below for an NaCl system
#!/usr/bin/perl
$name=$ARGV[0];
$i=1;
open(PDB,"${name}.xyz");
while($line=<PDB>){
if($line=~/ NA /){
chomp($line);
@data=split(/ +/,$line);
$atomType[$i]=$data[1];$atomX[$i]=$data[2];$atomY[$i]=$data[3];$atomZ[$i]=$data[4];
$i++;
}
if($line=~/ CL /){
chomp($line);
@data=split(/ +/,$line);
$atomType[$i]=$data[1];$atomX[$i]=$data[2];$atomY[$i]=$data[3];$atomZ[$i]=$data[4];
$i++;
}
}
open(CFG,">${name}.cfg");
$header="Converted from Packmol XYZ\n";
print CFG $header;
printf CFG "%10d%10d%10d%20.6f\n",0,0,0,0;
printf CFG "%20.12f%20.12f%20.12f\n",0,0,0;
printf CFG "%20.12f%20.12f%20.12f\n",0,0,0;
printf CFG "%20.12f%20.12f%20.12f\n",0,0,0;
for($i=1;$i<=$#atomType;$i++){
printf CFG "%-8s%10d\n",$atomType[$i],$i;
printf CFG "%20.6f%20.6f%20.6f\n",$atomX[$i],$atomY[$i],$atomZ[$i];
}
close(CFG); close(PDB);
For each new system being considered the conversion script needs to be modified. Below is the section that needs to be modified.
if($line=~/ NA /){
chomp($line);
@data=split(/ +/,$line);
$atomType[$i]=$data[1];$atomX[$i]=$data[2];$atomY[$i]=$data[3];$atomZ[$i]=$data[4];
$i++;
}
For each ion/molecule in the system this section needs to be in the conversion scirpt. So for NaCl there will be this section twice, once with NA on the first line, the second time with CL on the first line.
To run the conversion script on the HPC simply type ‘perl conv_xyz2conf.pl FILENAME’, without .xyz.
Once the conversion script has run a new file is generate, named FILENAME.cfg. Change this file name to CONFIG.
CONFIG file
The CONFIG file will look something like this
Converted from Packmol XYZ
0 0 0 0.000000
0.000000000000 0.000000000000 0.000000000000
0.000000000000 0.000000000000 0.000000000000
0.000000000000 0.000000000000 0.000000000000
NA 1
-2.187765 3.553628 0.510867
NA 2
0.656191 -5.190594 -3.045550
NA 3
-8.111144 -7.639522 5.169185
…
…
…
It is necessary to modify the top section of the CONFIG file. The CONFIG file is hugely format driven, changing a 0 to 20 means that its necessary to remove both a 0 and a space.
The number 0 indicated that for each atom in the simulation there is 1 piece of information, its location. This number does change if additional information is known for each atom.
The number 1 idicates cubic boundary conditions.
The 400 indicates the number of ions in the simulation, 200 Na and 200 Cl.
The final number 0.0000 is the energy of the system. It is 0 because it hasn't yet been determined.
The next 3 lines of the file indicate the box size. The box size is 10% larger than that generated in packmol. This is done to ensure that all atoms are housed within the box, packmol can generate a system with the end of a molecule (although not considered here) hanging out the edge of the box. The size increase is continued as it is good practice to allow the system to expand if needed in the initial steps.
After that are the atomic coordinates.
Once the top few lines of the CONFIG file is modified it should look something like this.
Converted from Packmol XYZ
0 1 400 0.000000
22.000000000000 0.000000000000 0.000000000000
0.000000000000 22.000000000000 0.000000000000
0.000000000000 0.000000000000 22.000000000000
NA 1
-2.187765 3.553628 0.510867
NA 2
0.656191 -5.190594 -3.045550
NA 3
-8.111144 -7.639522 5.169185
…
The next steps involve setting up the files that are used by DL_poly. The CONFIG file has already been generated, the other 2 files that ned to be generated at the FIELD file and the CONTROL file. The FIELD file is normally made from values found in the literature. The CONTROL file tells DL_poly what to run.
FIELD File
The first line in the FIELD file is the title, although I havent done it here, I’d recommend that the reference for the parameters used later in the file are inculded in the title line. This will allow you to look back at the work years later and easily find the paper.
The second line is simply the unit of the calculation.
molecules shows the number of molecules that are considered in the system, in NaCl there are 2 types, Na and Cl.
The next section is about each molecule/atom in the simulation. The name of the atom (NA) The number of NA atoms in the simulation (200) The number of atoms in the molecule (1) Atomic mass and charge (22.98... 1.00... 1) Finish – end of section on Na Same for Cl
vdW3 – the number of non-bonded interactions
The next section contains the parameters and the potential (BHM) being employed
900 NACL at 1200.00 K units kJ molecules 2 NA nummols 475 atoms 1 NA 22.9898000000 1.00000 1 finish CL nummols 475 atoms 1 CL 35.4530000000 -1.00000 1 finish vdW 3 CL CL BHM 0.26370 0.317 2.340 1.05 0.499 NA NA BHM 0.21096 0.317 2.755 6.99 8.68 CL NA BHM 0.15822 0.317 3.170 72.4 145.4 close
Submitting a simulation
To submit jobs the submission script (Dlrun.sh) is employed. To run the script you need to be in the directory with the CONTROL, CONFIG, FIELD and Dlrun.sh files. Using the terminal simply type ‘qsub Dlrun.sh’ and enter.
Submission script for DL_POLY
Within the submission script it is possible to edit the number of processors and the amount of memory that is required.
Visualising the Output
Right click on the XQuartz icon (back X with an orange ring around its centre) and open a terminal, via the applications tab. In this new terminal window sign into the HPC and navigate to the directory containing the simulation.
Once in the directory type module load vmd vmd
This will launch VMD via the HPC, so it isnt necessary to instal VMD on your local machine. Although it is always useful to have!
Once VMD has opened select ‘New Molecule’ from the ‘File’ tab at the top left and browse to find the HISTORY file. Change the file type to ‘DLPOLY V2 History’. If you don’t select the correct file type the program will freeze and you’ll need to start the whole procedure again.
Once the HISTORY file and correct file type are selected click load. On the black display screen there should be lots of small dots moving around. To change the size of the atoms to something that is more visible click on the ‘Graphics’ tab and select representations. Using the ‘Drawing Method’ drop down list select CPK.
Procedure
EACH SIMULATION NEEDS TO BE IN IT'S OWN DIRECTORY
The initial CONFIG file can be run at a temperature of about 300K. In the 'Getting started: generating a solvated structure and "relaxing" it' section of the Wiki, it is advised that the first run should be at 5K. As we are dealing with molten salts that don't melt until several hundred degrees, we can start at a higher temperature. The first run employs the NVE ensemble and is relatively short 500 steps and ensures that the parameters being employed are sensible, if they are not then the simulation will run into problems.
The next simulation run has been carried out above the melting temperature, 1200K, again employing the NVE ensemble. This run is longer at 60000 steps meaning 500ps with a tilmestep of 0.005. The starting CONFIG file has been generated from the REVCON file from the short run at 300K. The REVCON file needs to be renamed as CONFIG and the temperature in the CONTROL file changed.
Next a simulation employing the NVT ensemble has been carried out, again at 1200K using the REVCON file from the 1200K NVE simulation. This enables the size of the box to change accounting for the different temperatures that are being used.
Following the NPT simulation, the temperature is increased by 100K, to 1300K where an NVE simulation is carried out. The REVCON file from the 1200K NPT simulation is the new CONFIG file.
Starting at 1200K an NVE simulation followed by an NPT simulation has been completed. The temperature is then increased by 100K and an NVE followed by NPT simulation is completed. The CONTROL file is not changed, other than the temperature and the type of ensemble. This is repeated until the temperature is well above the melting point of the salt.
Checking the system is equilibrated
There are 3 sections that need to be commonly checked for basic analysis that should be carried out at the end of each simulation.
First thing to check is the energy of the system. Is The energy getting to a stable value before the ensemble freezes it in an NVE simulation? The second thing is the temperature. This will always vary, but should remain within 10% (ideally 5%) of the set value in an NPT simulation. The final thing to check in an NPT is the volume of the box. Has the volume changed and remains constant during the NPT simulation.
To get the data points to plot the graphs needed for the above checks a very useful script has been written by Vincent. [1] To get the script simply copy it from the link and save it as a new file call statis2xmgr.
to run the script simply type './statis2xmgr 1' this will give the energy data points in a file called stat2dl. It is necessary to change the file name to something else as each time you run the script it calls the file the same thing, overwriting what was there before. Typing no number at all after the script name gives the options available, temperature is 2 and volume is 19. To plot the points generated using statis2xmgr you can use gun plot. gnu plot is very powerful and I only ever use it to plot a scatter graph, ensuring the temp / energy / volume is not misbehaving. To ue gnu plot type 'gnu plot' in the HPC terminal, this will run gnuplot. then type "plot 'filename' " (you need the apostrophes) and it will plot that data and a popup will appear.
Calculating RDF
The method I am currently using is the in-built RDF from DL_poly.
In the CONTROL file simply type 'rdf' and 'print rdf' this will calculate and write a new file called RDFDAT. To bring the RDFDAT file back the submission script might need to be modified.
Modifying Dlrun.sh
simply include this line towards the end of the submission script (with all the other that look the same)
cp $TMPDIR/RDFDAT $PBS_O_WORKDIR/
This will copy the RDFDAT file back to your directory from the HPC where the calculation was completed.
Points to consider
Sometimes when you have modified a file and resubmit the calculations it appears that the modifications that you have made have not occurred. This is because the HPC only reads the files every minute or so, it is worth modifying the files, waiting a minute or two then re-submitting the modified file.
For the work I have carried out finding reliable parameters for the FIELD file caused some major issues. I expect most of the parameters being used come from literature. If this is the case it is worth trying to find multiple papers with the same parameters, or at east very similar. Many papers say they use parameters from paper X, but is worth trying to find two papers that state the parameters and not just reference them. I found that there are commonly problems in reporting the parameters in the correct units!
With a KCl system voids kept forming in the simultation box. I thought this would be due to problems with the parameters being used in the FIELD file. I tried several different sets of parameters, all to no avail. I tried changing the density (within what I could find in literature) which didn’t work. Eventually I tried taking a NaCl simulation that was working and changing the Na ions to K ions. It seems that the packmol structure generated, even for differeny box sizes (accounting for the different densities), didn’t change enough to allow the simulation to be succesful.
To change the Na ions to K the following command was used. ‘sed –i ‘s/NA/K /’ FILENAME’ you should notice that there is a space after the K, this is because DL_poly is very format driven. Replacing ‘Na’ with ‘K’ only replaces 1 of the 2 characters being removed.