View source for Third year simulation experiment/Running simulations under specific conditions

'''<big><span style="color:blue; ">This is the fourth section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Equilibration|Equilibration]], or jump ahead to the next section, [[Third_year_simulation_experiment/Structural_properties_and_the_radial_distribution_function| Structural Properties and the Radial Distribution Functions]].</span></big>'''

'''<big>THE FILES THAT YOU NEED FOR THIS SECTION ARE FOUND IN THE "NpT" SUBFOLDER.</big>'''

==Changing Ensemble==

So far, we have been able to do some simulations in which the number of particles and the volume of the simulation cell are held constant. The energy is also constant (within a certain degree of error, which is introduced by the approximations that we make to do the simulation). If the simulation is a working properly, then the pressure and temperature of the system should also reach a constant ''average'' value (although there will again be fluctuations). In the statistical thermodynamics lectures, you met the concept of ensembles, which are used in statistical mechanics to represent different sorts of experimental conditions. The simulations we have done so far are described by the ''microcanonical'', or NVE ensemble (the letters represent those thermodynamic quantities which are constant).

As chemists, we often want to understand what happens under particular experimental conditions &mdash; at 298K under 1 atmosphere of pressure, for example. These sorts of conditions are described by different ensembles in statistical mechanics, such as the NVT (''canonical'') or NpT (''isobaric-isothermal'') ensembles.

In this section, we are going to modify our simulations from the previous section to run under NpT conditions, and sketch an equation of state for our model fluid at atmospheric pressure.

==Temperature and Pressure Control==

The file npt.in can be used to perform a constant temperature/pressure simulation of our model fluid. It starts by melting a simple cubic crystal, just as before, so much of this file will look familiar to you. You will notice a new section near the top, however, called '''### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###'''. It contains three ''variables'' &mdash; these are used by the script later on to define the desired temperature, pressure, and timestep. The ellipses need to be replaced by the actual temperature, pressure and timestep that you want to use, so

<pre>
variable T equal 0.5
variable p equal 1.0
variable timestep equal 0.75
</pre>

would run a simulation at <math>T=0.5,\  p=1.0,\  \delta t=0.75</math>. You should remember from the [[Third_year_simulation_experiment/Equilibration|Equilibration]] section that this is a poor choice of timestep!

'''<big>TASK 8: Choose 5 temperatures (above the critical temperature <math>T^* = 1.5</math>), and two pressures (you can get a good idea of what a reasonable pressure is in Lennard-Jones units by looking at the average pressure of your simulations from the last section). This gives ten phase points &mdash; five temperatures at each pressure. Create 10 copies of npt.in, and modify each to run a simulation at one of your chosen <math>\left(p, T\right)</math> points. You should be able to use the results of the previous section to choose a timestep. Submit these ten jobs to the HPC portal. When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density <math>\left(\frac{N}{V}\right)</math>. </big>'''

=== Plot the density as a function of temperature for both pressures that you simulated. Include a line corresponding to the predictions made by the ideal gas law. [3] ===

=== How do your results compare to the ideal gas law? Do deviations increase/decrease with temperature and pressure? Explain. [7] ===

=== Do you expect your simulation results to be in better or worse agreement with the Van der Waals equation of state? Why? [3] ===

===Thermostats and Barostats - controlling the thermodynamic properties===
The statistical thermodynamics lectures will have introduced you to the ''equipartition theorem'', which states that, on average, every degree of freedom in a system at equilibrium will have <math>\frac{1}{2}k_B T</math> of energy. In our system with <math>N</math> atoms, each with 3 degrees of freedom, we can write
<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

At the end of every timestep, we use the left hand side of this equation to calculate the kinetic energy, then divide by <math>\frac{3}{2}Nk_B</math> to get the ''instantaneous'' temperature <math>T</math>. In general, <math>T</math> will fluctuate, and will be different to our ''target'' temperature, <math>\mathfrak{T}</math> (this is whatever value we specify in the input script). We can change the temperature by multiplying every velocity by a constant factor, <math>\gamma</math>.

* If <math> T > \mathfrak{T} </math>, then the kinetic energy of the system is too high, and we need to reduce it. <math>\gamma < 1</math>
* If <math> T < \mathfrak{T} </math>, then the kinetic energy of the system is too low, and we need to increase it. <math>\gamma > 1</math>

We need to choose a scaling parameter <math>\gamma</math> so that the temperature is correct <math>T = \mathfrak{T}</math> if we multiply every velocity <math>\gamma</math>. We can write two equations:'''
<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}</math>

By combining these equations, one can see that <math> \gamma = \sqrt{\frac{\mathfrak{T}}{T}} </math> (satisfy yourself that this is true!). A target value of <math> \gamma </math> of 1 is required and thus, dependent on whether it's larger or smaller than 1 the simulation can target the desired temperature.

Controlling the pressure is a little more involved, but the principle is largely the same: at each timestep, the pressure of the system is calculated; if the pressure is too high, then the simulation box is made a little larger, while if the pressure is too low the box is made smaller. Simulations in which the pressure is controlled are thus in the NpT ensemble &mdash; the volume of the simulation box is not constant!

===Examining the Input Script===

Open one of your input scripts (it doesn't matter which), and look at the section '''### BRING SYSTEM TO REQUIRED STATE ###'''. The line <pre>fix npt all npt temp ${T} ${T} ${tdamp} iso ${p} ${p} ${pdamp}</pre> is the one responsible for switching on the temperature and pressure control. LAMMPS actually allows us to heat or cool the system over the course of a simulation, if we want to &mdash; this is the reason that the temperature appears twice in this line. The first ${T} is the desired starting temperature, and the second is the desired temperature at the end of the simulation. We want a constant average temperature, so we specify the same value twice. The same goes for the pressure.

Now look at the lines near the end of the file:

<pre>
### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2
run 100000
</pre>

The first command, ''thermo_style'', controls which thermodynamic properties are recorded, as before. The next lines are used to measure ''average'' thermodynamic properties for the system. To draw our equations of state, we need to know the average temperature, pressure, and density, and the statistical errors in those quantities. The six variable lines link those quantities (and their squared values, needed for the errors), to variable names that we can use in the averaging command, which is the line starting ''fix aves...''. This command takes a number of input values and averages them every so many timesteps. Exactly how often this happens depends in the values of the three numbers which follow ''ave/time''.

'''<big><span style="color:blue; ">This is the fourth section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Equilibration|Equilibration]], or jump ahead to the next section, [[Third year simulation experiment/Structural properties and the radial distribution function|Structural properties and the radial distribution function]].</span></big>'''