exercises:2015_cecam_tutorial:mtd1
Differences
This shows you the differences between two versions of the page.
Both sides previous revisionPrevious revisionNext revision | Previous revision | ||
exercise:simple_metadynamics_simulation_using_the_coordination_numbers_as_variables [2015/07/20 08:51] – style corrections and adding links for keywords/parameters/sections tmueller | exercises:2015_cecam_tutorial:mtd1 [2020/08/21 10:15] (current) – external edit 127.0.0.1 | ||
---|---|---|---|
Line 1: | Line 1: | ||
+ | ====== Simple metadynamics simulation using the coordination numbers as variables ====== | ||
+ | |||
Problem: Dissociation reaction of nitric acid on graphene and atomic rearrangements of a < | Problem: Dissociation reaction of nitric acid on graphene and atomic rearrangements of a < | ||
+ | |||
+ | * Original author: Marcella Iannuzzi | ||
+ | * Complete source and output files: [[http:// | ||
+ | |||
+ | ===== Introduction ===== | ||
For this tutorial some input and output files are given in order to present a complete procedure to solve the given problem. Some hints are also given to help in the analysis of the results. In order to be able to run these examples, some paths need to be correctly set in the input files (i.e. set the variables '' | For this tutorial some input and output files are given in order to present a complete procedure to solve the given problem. Some hints are also given to help in the analysis of the results. In order to be able to run these examples, some paths need to be correctly set in the input files (i.e. set the variables '' | ||
Line 8: | Line 15: | ||
* Set up and preliminary simulations to learn about the dynamics of nitric acid on graphene, as obtained at the DFTB level of theory | * Set up and preliminary simulations to learn about the dynamics of nitric acid on graphene, as obtained at the DFTB level of theory | ||
- | * Mmetadynamics | + | * Metadynamics |
* Set up and preliminary simulations to learn about the dynamics of the selected small Si cluster saturated by H atoms | * Set up and preliminary simulations to learn about the dynamics of the selected small Si cluster saturated by H atoms | ||
- | * Metadynamics simulation aimed at observing atomic rearrangements of the cluster bay changing the coordination of both Si and H species. | + | * Metadynamics simulation aimed at observing atomic rearrangements of the cluster bay changing the coordination of both Si and H species. |
| | ||
- | ====== First task: dynamics of two HNO_3 molecules over a graphene sheet ====== | + | ===== First task: dynamics of two HNO3 molecules over a graphene sheet ===== |
| | ||
The examples on this system are in the directory '' | The examples on this system are in the directory '' | ||
However, the reaction should be catalyzed in presence of C particles (sooth). In the proposed example, the molecules are located in the vicinity of graphene, that should mimic the role of sooth. Graphene indeed is not very reactive, better models can be considered using defective or functionalized graphene. | However, the reaction should be catalyzed in presence of C particles (sooth). In the proposed example, the molecules are located in the vicinity of graphene, that should mimic the role of sooth. Graphene indeed is not very reactive, better models can be considered using defective or functionalized graphene. | ||
| | ||
- | The initial coordinates are given in '' | + | The initial coordinates are given in '' |
< | < | ||
Line 54: | Line 61: | ||
&END | &END | ||
</ | </ | ||
- | Otherwise the '' | + | Otherwise the '' |
- | The output of this short simulation (5 ps) is stored in DFTB_NVT. Kinetic energy (3rd col.), temperature (4th col.), potential energy (5th col), and total energy (6th col) can be monitored from '' | + | The output of this short simulation (5 ps) is stored in '' |
< | < | ||
Line 78: | Line 85: | ||
The next step is to set up a few collective variables (CV) that can be later used for metadynamics (MTD) simulations. It is important to select good CV that can describe the relevant configuration along the reaction path. Moreover, it is useful to learn about the typical behavior of the selected CV along an unbiased MD run. Hence, after selecting a set of CV, preliminary runs should be performed in order to monitor the dynamics of these variables. This can be done setting up MTD simulations, | The next step is to set up a few collective variables (CV) that can be later used for metadynamics (MTD) simulations. It is important to select good CV that can describe the relevant configuration along the reaction path. Moreover, it is useful to learn about the typical behavior of the selected CV along an unbiased MD run. Hence, after selecting a set of CV, preliminary runs should be performed in order to monitor the dynamics of these variables. This can be done setting up MTD simulations, | ||
- | The input '' | + | The input '' |
< | < | ||
Line 147: | Line 154: | ||
</ | </ | ||
- | This last CV controls the distance of the molecules from the layer, which is an important factor to determine whether the dissociation is somehow favored by the presence of graphene. The output of the unbiased MD that monitors the behavior of these 4 CV is in DFTB_MTD_4CV_H0. It is obtained by invoking a MTD run where no penalty potential is added. Hence, the FREE_ENERGY subsection is added within the section MOTION. The MTD rum is controlled from the '' | + | This last CV controls the distance of the molecules from the layer, which is an important factor to determine whether the dissociation is somehow favored by the presence of graphene. The output of the unbiased MD that monitors the behavior of these 4 CV is in '' |
< | < | ||
Line 186: | Line 193: | ||
| | ||
</ | </ | ||
- | With '' | + | With '' |
By plotting the CV as recorded along the short MD trajectory (3 ps), the amplitude of the equilibrium fluctuations can be evaluated and then used to set up the size of the Gaussian hills that build up the biasing potential. The first CV fluctuates close to zero, with fluctuations smaller than 0.2. The second is around 2.8. The fluctuations are smaller due to the stiffness of the three NO bonds. The coordination of H to C is also typically zero, but it can change a lot when the molecules approach the layer, even if there is no dissociation of H and no binding to C. This indicates that this variable is difficult to control and might turn out to be tricky to use it to distinguish among different states of the reaction process. The point to plane distance shows quite large fluctuations and it is clearly not suited to distinguish a specific state along the reaction path. Moreover, its minima, when the two molecules are closer to the layer, correspond to the maxima of the third CV, i.e. the CN of H to C. At least before dissociation, | By plotting the CV as recorded along the short MD trajectory (3 ps), the amplitude of the equilibrium fluctuations can be evaluated and then used to set up the size of the Gaussian hills that build up the biasing potential. The first CV fluctuates close to zero, with fluctuations smaller than 0.2. The second is around 2.8. The fluctuations are smaller due to the stiffness of the three NO bonds. The coordination of H to C is also typically zero, but it can change a lot when the molecules approach the layer, even if there is no dissociation of H and no binding to C. This indicates that this variable is difficult to control and might turn out to be tricky to use it to distinguish among different states of the reaction process. The point to plane distance shows quite large fluctuations and it is clearly not suited to distinguish a specific state along the reaction path. Moreover, its minima, when the two molecules are closer to the layer, correspond to the maxima of the third CV, i.e. the CN of H to C. At least before dissociation, | ||
- | ====== Second task: Metadynamics of the dissociation of HNO_3 over a graphene sheet ====== | + | ===== Second task: Metadynamics of the dissociation of HNO3 over a graphene sheet ===== |
- | The presented MTD run employs as CV only the three CN described above. The related input file is '' | + | The presented MTD run employs as CV only the three CN described above. The related input file is '' |
< | < | ||
Line 233: | Line 240: | ||
One Gaussian hill is spawned every NT_HILLS MD steps, while the height of the hill in hartree is given by WW. These parameters together with the width of the Gaussian hills are important to determine the accuracy of the description of the FES through the MTD biasing potential. Since each variable has in principle different dimensions and different dynamics, the shape of the hills filling up the $N_{\text{CV}}$-dimensional configurations space, as defined by the selected CVs, is not isotropic. The parameter SCALE associated to the $i$-th MTD variable determines the amplitude of the Gaussian in the $i$-th space-direction of the $N_{\text{CV}}$-dimensional configuration space. This parameter, as well as the hill’s height and the frequency of collocation, | One Gaussian hill is spawned every NT_HILLS MD steps, while the height of the hill in hartree is given by WW. These parameters together with the width of the Gaussian hills are important to determine the accuracy of the description of the FES through the MTD biasing potential. Since each variable has in principle different dimensions and different dynamics, the shape of the hills filling up the $N_{\text{CV}}$-dimensional configurations space, as defined by the selected CVs, is not isotropic. The parameter SCALE associated to the $i$-th MTD variable determines the amplitude of the Gaussian in the $i$-th space-direction of the $N_{\text{CV}}$-dimensional configuration space. This parameter, as well as the hill’s height and the frequency of collocation, | ||
- | The '' | + | The '' |
The provided trajectory is about 100 ps long and indeed it shows the dissociation of the two molecules into < | The provided trajectory is about 100 ps long and indeed it shows the dissociation of the two molecules into < | ||
- | Other quantities that can be monitored from the '' | + | Other quantities that can be monitored from the '' |
- | ====== Third task: dynamics of Si_6H_8 ====== | + | ===== Third task: dynamics of Si6H8 ===== |
- | The data file for this example are in SI6_CLU. In this case, a small Si cluster of 6 Si atoms saturated by 8 H atoms is studied. Si clusters show different arrangements. The equilibrium structure should be such that Si atoms keep the preferred tetrahedral coordination. In the presence of H saturating the dangling Si bonds, the structure can be open, like the chair structure that is used here as starting conformation. By loosing H atoms, through the formation of molecular hydrogen, the cluster undergoes some rearrangement. The structure should become more compact in order to saturate the Si coordination shell. | + | The data file for this example are in '' |
- | As in the previous example, preliminary MD runs are carried out to study the dynamics of the system. The electronic structure is computed at the DFT level, using the PBE functional. A standard constant temperature MD is simulated running the input si6_clu_nvt.inp. The temperature is set at 300 K and it is controlled by applying the Nose-Hoover thermostat and the temperature rescaling (the small number of degrees of freedom does not allow thermalization within a few ps). The output of this simulation is in NVT. The energy curve over the 6 ps of simulation time and the visualization of the trajectory confirm that the cluster with the initial stoichiometry can be easily equilibrated in the chair structure. | + | As in the previous example, preliminary MD runs are carried out to study the dynamics of the system. The electronic structure is computed at the DFT level, using the PBE functional. A standard constant temperature MD is simulated running the input '' |
- | Possible changes in stoichiometry (loosing H) and structure are going to induce variations in the coordination shell of the Si atoms. Therefore, three CN are selected as CV : Si to Si, Si to H, and H to H. These variables are monitored over the equilibrium trajectory by running the input si6_clu_mtd_h0_p1.inp. As in the previous example, this is a MTD input where the the variable DO_HILLS is set to false, so that no biasing potential is added. The output obtained running this test is in MTD_H0. By monitoring the three CN along the 10 ps log trajectory (2nd, 3rd and 4th col. of si6_clu_mtd_h0-COLVAR.metadynLog), | + | Possible changes in stoichiometry (loosing H) and structure are going to induce variations in the coordination shell of the Si atoms. Therefore, three CN are selected as CV : Si to Si, Si to H, and H to H. These variables are monitored over the equilibrium trajectory by running the input '' |
< | < | ||
Line 260: | Line 267: | ||
the Si-Si CN oscillates around 2. Actually, in the chair configuration 4 Si are three fold coordinated with neighboring Si atoms and 2 are two fold coordinated. The bond length is about 2.4 Å. By changing the curvature of the function, the average value of the CN can be easily moved towards 3 and it can become more sensitive to fluctuations of the Si-Si bond length. The fluctuations of the two other CN are even smaller. Si-H fluctuates around 1.4 and H-H is very close to zero, since in the initial configuration the H atoms do not see each other. | the Si-Si CN oscillates around 2. Actually, in the chair configuration 4 Si are three fold coordinated with neighboring Si atoms and 2 are two fold coordinated. The bond length is about 2.4 Å. By changing the curvature of the function, the average value of the CN can be easily moved towards 3 and it can become more sensitive to fluctuations of the Si-Si bond length. The fluctuations of the two other CN are even smaller. Si-H fluctuates around 1.4 and H-H is very close to zero, since in the initial configuration the H atoms do not see each other. | ||
- | The output of a second simulation that monitors the same three CVs is in MTD_L_H0, and the corresponding input is si6_clu_mtd_l_h0_p1.inp. Nothing has been changed in the definition of the CVs, but the Lagrangian MTD formalism has been used. With this scheme, an auxiliary variable is associated to each CV, and when the biasing potential is added, it is defined as function of the auxiliary variables rather than of the CVs. The auxiliary variable behaves as additional degree of freedom. Therefore, an inertial mass is associated to it and its dynamics is determined by integrating the same type of equations of motion as for all the other degrees of freedom. The variable is coupled to the corresponding CV through a harmonic potential, and the forces acting on it are those derived from the harmonic potential and from the MTD biasing potential, when it is present. Hence, in the METAVAR input section, two additional parameters are needed, which are the mass of the auxiliary variable and the coupling constant for the harmonic potential: | + | The output of a second simulation that monitors the same three CVs is in MTD_L_H0, and the corresponding input is '' |
< | < | ||
Line 273: | Line 280: | ||
A temperature is associated to the auxiliary variables and can be controlled by temperature rescaling. The use of thermostats for such few degrees of freedom is questionable. The Lagrangian MTD formalism is used in order to better control the kinetics of the CV. This control is obtained through the coupled to the auxiliary variables, whose dynamics depends on the mass and the temperature, | A temperature is associated to the auxiliary variables and can be controlled by temperature rescaling. The use of thermostats for such few degrees of freedom is questionable. The Lagrangian MTD formalism is used in order to better control the kinetics of the CV. This control is obtained through the coupled to the auxiliary variables, whose dynamics depends on the mass and the temperature, | ||
- | Whenever a METAVAR is defined (also in a not Lagrangian scheme), it is possible to limit the range of values that are going to be explored by setting a so called WALL potential. This is typically done to avoid the time-consuming exploration of regions of the configurations space that are not relevant for the process that is under investigation. When auxiliary variables are employed, the WALL potential can also be activated in order to avoid the exploration of unphysical values, that can happen by too large fluctuations away from the corresponding CV. In the case of the CN, negative values are unphysical and must be avoided. For this reason the subsection WALL is added METAVAR coupled to the H-H CN, which is known to oscillate close to zero in the initial state: | + | Whenever a '' |
< | < | ||
Line 286: | Line 293: | ||
</ | </ | ||
- | In this case the potential function is $f(s)=K (s-s_0) ^4$, and it is activated whenever the variable becomes lower (DIRECTION WALL_MINUS) than the given $s_0$ (POSITION). | + | In this case the potential function is $f(s)=K (s-s_0) ^4$, and it is activated whenever the variable becomes lower ('' |
When such Lagrangian scheme is used more columns appear in the COLVAR, which contain all the relevant information. The 1st column is always the time in fs, the next $N_{\text{CV}}$ columns are the instantaneous values of the auxiliary variables, followed by $N_{\text{CV}}$ columns where the instantaneous values of the CVs are reported. The next columns report the values of the potential gradients: $N_{\text{CV}}$ columns for the gradients of the harmonic potential, $N_{\text{CV}}$ for the gradients of the MTD biasing potential, and $N_{\text{CV}}$ for the gradient of the WALL potential (these are zeros when the corresponding potential is not activated). The following $N_{\text{CV}}$ are the velocities of the auxiliary variables. Then there are the instantaneous values of the harmonic potential, of the MTD potential, and of the WALL potential. The last column is the temperature of the auxiliary variables. | When such Lagrangian scheme is used more columns appear in the COLVAR, which contain all the relevant information. The 1st column is always the time in fs, the next $N_{\text{CV}}$ columns are the instantaneous values of the auxiliary variables, followed by $N_{\text{CV}}$ columns where the instantaneous values of the CVs are reported. The next columns report the values of the potential gradients: $N_{\text{CV}}$ columns for the gradients of the harmonic potential, $N_{\text{CV}}$ for the gradients of the MTD biasing potential, and $N_{\text{CV}}$ for the gradient of the WALL potential (these are zeros when the corresponding potential is not activated). The following $N_{\text{CV}}$ are the velocities of the auxiliary variables. Then there are the instantaneous values of the harmonic potential, of the MTD potential, and of the WALL potential. The last column is the temperature of the auxiliary variables. | ||
- | The dynamics of the CVs along the two simulations, | + | The dynamics of the CVs along the two simulations, |
- | ====== Fourth task: Lagrangian MTD of the atomic rearrangement of Si6H8 ====== | + | ===== Fourth task: Lagrangian MTD of the atomic rearrangement of Si6H8 ===== |
- | '' | + | '' |
exercises/2015_cecam_tutorial/mtd1.1437382262.txt.gz · Last modified: 2020/08/21 10:14 (external edit)