Section 2. Temperature and pressure
Video 1. Section 2 summarized into a 1:30 minute video. For a detailed explanation, follow the tutorial below.
Our Argon atoms are now being simulated without the need of energy caps to prevent the solver from diverging, but we still need to set the temperature of the system to \(T = 85 K\). Temperature in our simulation can be estimated using the Equipartition Theorem, relating
kinetic energy (\(E_k\)) with temperature (\(T\))
$$T=\left.\left\langle\sum_{m}\frac{{2E}_k}{{3k}_bN}\ \right.\right\rangle\ ,$$
where \(E_k\) is the kinetic energy, \(k_b\) is the Boltzman constant, and \(N\) is the number of particles.
It is very likely that the temperature of the current system is exceedingly high (or very far from being at 85 K), and we still need to change the energy of the system to reduce the temperature.
This can be achieved by rescaling the velocity of the particles by a factor \(\gamma\) so that the temperature is set to a desired value. The factor \(\gamma\) is computed using the following expression
$$\gamma=\frac{3{Nk}_bT_{set}}{2E_k},$$
where \(E_k\) is the kinetic energy of the particles at equilibrium (before rescaling), \(T_{set}\) is the desired system temperature, \(k_b\) is the Boltzman constant, and \(N\) is the number of particles.
Next, we present two different methods that achieve the same objective using Scymol. The first and faster method is done by using Scymol's ability to automatically run multiple simulations in series while rescaling the velocities between each simulation until the temperature value converges
to the desired value of \(T_{set}\). The second and slower method is done by manually running several simulations in series while applying the velocity changes yourself using the Modify tool in Scymol.
Scymol offers an option to run multiple simulations in series while executing a custom Python code between each simulation, which in our case would be to rescale the velocities of the particles so that we reach the desired temperature \(T_{set} = 85K\). It works the same way as if you were to manually run a simulation, multiply the particles' velocities by a factor \(\gamma\), and then run another simulation. Go into the Run tab and type in '10' in the Number of simulations to run in series input box. Ten simulations is usually enough for the value of \(T\) to reach a value of \(T_{set} = 85K\) for this type of system. Now, go into the Dynamics' Solver tab and copy and paste Code 1 onto the Custom Python code to execute after each simulation input box:
flt_kinetic_energy = 0 flt_kb = 0.836031 flt_Tset = 85 flt_T = 0 for dct_particle_index in dct_every_particle_info[-int_total_nbr_of_particles:]: flt_kinetic_energy += 0.5 * dct_particle_index['Mass'] * ( dct_particle_index['Vel.x'] ** 2 + dct_particle_index['Vel.y'] ** 2 + dct_particle_index['Vel.z'] ** 2) flt_gamma_factor = (3 * flt_kb * flt_Tset) / (2 * flt_kinetic_energy / int_total_nbr_of_particles) flt_T = 2 * flt_kinetic_energy / (3 * flt_kb * int_total_nbr_of_particles) print('Current Temperature: ', flt_T) if flt_T < 83 or flt_T > 87: for dct_particle_index in dct_every_particle_info[-int_total_nbr_of_particles:]: dct_particle_index['Vel.x'] *= flt_gamma_factor dct_particle_index['Vel.y'] *= flt_gamma_factor dct_particle_index['Vel.z'] *= flt_gamma_factor
Code 1. This Python code rescales the velocities of the particles on the last step of a simulation by a factor \(\gamma\). It has to be applied several times (running equilibration steps in between).
In Code 1 we are:
1) Estimating the kinetic energy of the system.
2) Calculating the \(\gamma\) factor and the current system's temperature \(T\).
3) Multiplying the particles' velocities by \(\gamma\) if \(T\) is not close to \(T_{set} = 85K\).
For more information on how to write custom Python code into the program and which variables are intrinsic to Scymol please refer to XX. Go ahead and click on Run. After it finishes, go to the Analyze tab and verify that the temperature \(T\) is around 85 K.
Before rescaling the particles' velocities, we have to estimate the current system's kinetic energy so that we can compute \(\gamma\). Navigate into the Analyze tab, select the kinetic energy option from the Physical drop-down menu, and click on the check to the right to plot and get the values of the kinetic energy. Make sure that you obtain the average of the last ~5 values obtained; this is achieved by selecting a range using the boxes at the top of the Analyze tab: From frame A to Frame B. Now that we have computed a \(\gamma\) factor, we need to go ahead and multiply the particles' velocities by this number. Make sure that you are on the last frame of the simulation. Go into the Edit→Select and select every particle by clicking on the Select All button. Then go to the Modify tab and select the Multiply option from the drop-down menu on the top. Then press three times on the Plus sign on the bottom. Select Vel.x, Vel.y, and Vel.z from each drop-down menu that appeared. Type in the \(\gamma\) factor that you computed into the input box next to the drop-down menus. Press on the Check button to apply the changes. The particles' velocities have been multiplied by \(\gamma\). Go back into the Run tab and run the simulation again.
After it finishes, go into the Analyze tab and check if the temperature of the system is indeed around 85 K. If it is not (which is very unlikely to happen just after one simulation step), go back to the beginning of this paragraph and repeat the steps. After ~10 simulation steps, the temperature should equilibrate around \(T_{set} = 85K\).
Before proceeding to section 3, it is important to make sure that the temperature of our system is stable at around 85 K without any further temperature control. To do this, we run 1 simulation without any velocity rescaling for 50 ps. Go into the Run tab and change the number of simulations to be run in series to just 1 and the simulation time to 50 ps. Click on Run. When the simulation finishes, make sure that the temperature of the system at the last frame is indeed close to 85 K (between 81 and 89 K) it is fine. If it is anything outside this range, apply some additional velocity rescaling steps until the temperature stabilizes around 85 K.