MD Model Potential
We are going to run all the simulations in the HPC facility of IIT Kanpur. For this you need to login into your HPC account by typing
ssh -X username@hpc2010.hpc.iitk.ac.in
Then enter your password, which is already given to you.
The files required for running the simulation are already present in the home folder of your account.
For today's tutorial go to the folder T1_MD by following command
cd T1_MD/
Introduction
The aim of this tutorial is to run a molecular dynamics (NVT) simulation on a model system and to generate the free energy surface from the probability distribution.
The model system:
The model system that we are going to use for our simulation is a one-dimensional double well potential with barrier equal to 0.002 a.u. and having two minima at -1 and 1 (shown below). The potential has form
This function can be plotted by the gnuplot software.
For opening gnuplot, type gnuplot in the terminal. Once gnuplot terminal appears, then you need to type the following:
load "01-plot.gpi"
The following plot will appear on your screen.
To close gnuplot, type q in the gnuplot terminal and press ENTER.
Now, we will perform MD simulation (NVT) on this model potential energy surface at a temperature of 65 K (kBT<<barrier height).
Enter the directory MD_65K by typing
cd MD_65K/
In this folder, we have an input file 01-input_md that contains following input parameters:
01-input_md
1.d0 #initial position
10.d0 #mass(in a.u.)
6.d0 #timestep(in a.u)
65.d0 #temperature at which the simulation should be performed
.T. #.T. or .F. depending on the case if Langevin thermostat should be on or off
0.000001 #friction coefficient to be used for Langevin thermostat
99000000 #number of steps
10 #frequency at which the information should be printed in the output files
.F. #.T. or .F. if the restart is required from the previous simulation
Now, we are going to submit the run for MD simulation by the following command:
qsub sge_submit.sh
You can check the status of your run by typing:
qstat -u username
After the run is over, the files TRAJ, AVGTEMP and ENERGIES will be generated as output containing the following information:
- TRAJ: step number, position, velocity, potential energy
- AVGTEMP: step number, instantaneous temperature
- ENERGIES: step number, potential energy, kinetic energy, total energy
- free_energy: coordinate, probability, free energy
Before analysing the results, we are going to perform the same MD simulation at a higher temperature (650 K; kBT>barrier height).
Now go to the directory MD_650K
cd ../MD_650K
Now, we are going to submit the run for MD simulation by the following command:
qsub sge_submit.sh
Once the runs are over, we will transfer the data to your system.
For this we need to open a new terminal and make a new directory Workshop_Tutorial
mkdir Workshop_Tutorial
cd Workshop_Tutorial
Now copy the data from HPC to our system by the following command
scp -rp username@hpc2010.hpc.iitk.ac.in:/home/username/T1_MD .
Analysis of Results:
Now, we can see the trajectory of the particle moving on the potential energy surface. This can be visualized from the TRAJ file, by plotting the coordinates versus potential energy together with the function f(x).
In case of simulation at 65 K,
- we need to go to the directory MD_65K.
- Open gnuplot
- Load the script 02-plot.gpi
Here we observe the particle is stuck in one of the well and the barrier crossing did not take place.
In order to plot the probability as a function of x, load the script 01-prob.gpi in gnuplot.
To plot free energy as a function of x, load the script 02-prob.gpi in gnuplot
In the case of low temperature, it can be seen that the free energy surface is mimicking the bottom of one of the wells, but not the full potential energy landscape.
In case of simulation at 650 K,
- we need to go to the directory MD_650K.
- Open gnuplot
- Load the script 02-plot.gpi
In order to plot the probability as a function of x, load the script 01-prob.gpi in gnuplot.
The free energy as a function of x can be reconstructed in the similar way.
To plot free energy as a function of x, load the script 02-prob.gpi in gnuplot