Differences

This shows you the differences between two versions of the page.

--- howto:biochem_qmmm [2017/07/25 11:11] – dvanrompaey
+++ howto:biochem_qmmm [2019/08/22 09:23] – sllabres
@@ Line 69: / Line 69: @@
 Visualize your simulation to check for any abnormalities. Take a look at the output file, paying extra attention to any warnings you might receive. Note: CP2K will warn you about missing forcefield terms. You can have CP2K output these using ''[[inp>FORCE_EVAL/MM/PRINT/FF_INFO]]''.In this case, the missing terms are Urey-Bradley interactions. This is normal, as the Amber force field we are employing here does not use Urey-Bradley terms. A second warning states that our CRD file lacks velocities and box information will not be read. This is not a problem either, as we do not have velocities at this moment and box parameters are provided in the input file.
-In the next step we will perform 5 ps of dynamics in the NVT ensemble, to equilibrate the temperature of our simulation using this ''nvt.inp''. We are using the ''[inp>EXT_RESTART]'' section to restart with the last coordinates of the energy minimization procedure. Once the simulation has finished, take a look at the temperature file ''NVT-1.PARTICLES.temp'' to check that the temperature is oscillating around the correct value (in this case, 298K). Note that we are using the CSVR thermostat with a low time constant of 10 fs to efficiently and quickly control the system temperature. For production simulations a higher time constant such as 100 fs should be chosen to avoid interfering with the dynamics of the system.
+In the next step we will perform 5 ps of dynamics in the NVT ensemble, to equilibrate the temperature of our simulation using this ''nvt.inp''. We are using the ''[[inp>EXT_RESTART]]'' section to restart with the last coordinates of the energy minimization procedure. Once the simulation has finished, take a look at the temperature file ''NVT-1.PARTICLES.temp'' to check that the temperature is oscillating around the correct value (in this case, 298K). Note that we are using the CSVR thermostat with a low time constant of 10 fs to efficiently and quickly control the system temperature. For production simulations a higher time constant such as 100 fs should be chosen to avoid interfering with the dynamics of the system.
 The next step is simulating our system in the constant pressure NPT ensemble to get the correct dimensions for our cell. ''npt.inp'' runs 50 ps of NPT, outputting the cell dimensions into a separate file every 100 timesteps. A C-C distance in of the chorismate molecules is constrained, in order to keep it near the reactive conformation to facilitate our later simulations in the QM/MM phase. Pressure control is achieved by adding a barostat with a reference pressure of 1 bar and a coupling time of 100fs, and changing ''[[inp>MOTION/MD#ENSEMBLE]]'' from NVT to NPT_I, indicating that the box will be scaled isotropically. We expect our cell to shrink at first, as the solvation procedure cannot fill up the cell completely due to overlap with the protein and the edges of the box. After a while the cell size should stabilize. Take a look at the cell dimensions in ''NPT.cell''; we will be using these equilibrated cell dimensions for the next tutorial.
@@ Line 93: / Line 93: @@
 </code>
-Now that the QM atoms are defined, CP2K needs to know how to treat these. We will use the semi-empirical AM1 method as our QM method, electrostatically embedded in the MM system: ''[[inp>FORCE_EVAL/QMMM#E_COUPL]]'' COULOMB. This allows for the QM region to be polarized by the MM environment. Mechanical embedding, where the QM region only interacts with the MM region as point charges and no polarization occurs, can be selected through ''E_COUPL NONE''. In DFT calculations, the highly efficient [[doi>10.1021/ct050123f |GEEP]] method can be used as well. For stability reasons, we need to make some slight modifications to our prmtop file. In the TIP3P water model, the hydrogens do not have a separate Lennard-Jones site. The same is true for the amber serine, threonine and tyrosine HO atoms. This can cause unphysical interactions with the QM system: the hydrogen atoms strongly interact with the QM electron clouds, crashing the system. In order to prevent this we will add LJ sites to water hydrogens and serine/threonine/tyrosine hydrogens with the ambertools utility parmed. We will change their LJ parameters to those of a GAFF2 type alcohol using the ''changeLJSingleType'' command. Take care when using this, as it changes the LJ type for all of the atoms that share the LJ type, not just the one you selected. For instance, changing the serine HO also alters all the threonine and tyrosine HO atoms in the system. In this case, this behaviour is useful, but it is strongly advised to use ''printLJTypes'' to get a good idea of what exactly you are modifying. For a more selective manner of doing this, the ''addLJType'' command is useful. This adds an LJ type for your selected atom or group of atoms, rather than modifying all occurrences of a specific LJ type. Confirm that ''@5724'' corresponds to a water hydrogen and that ''@4843'' is a serine HO using ''printDetails @(atomnumber)'' in parmed. If not, use ''printDetails'' to locate the correct atom numbers and modify the changeLJSingleType commands accordingly.
+Now that the QM atoms are defined, CP2K needs to know how to treat these. We will use the semi-empirical AM1 method as our QM method, electrostatically embedded in the MM system: ''[[inp>FORCE_EVAL/QMMM#E_COUPL]]'' COULOMB. This allows for the QM region to be polarized by the MM environment. Mechanical embedding, where the QM region only interacts with the MM region as point charges and no polarization occurs, can be selected through ''E_COUPL NONE''. In DFT calculations, the highly efficient [[doi>10.1021/ct050123f |GEEP]] method can be used as well.
-<code>
+Before running the QM/MM simulation we need to amend our prmtop file, because in the classical forcefield, several hydrogen atom types do not have separate Lennard-Jones parameters or they are set up to 0.0. The AMBER FF atom types are HO (from the TIP3P water model), HG (from serine, threonine and tyrosine residues). This will cause unphysical interactions with the QM region and the simulation will fail due for stability reasons. The hydrogen atoms with these atom types will strongly interact with the QM electron clouds, eventually crashing the system. To prevent this from happening, we have to add the correct LJ parameters to the mentioned atom types using the ambertools utility parmed.
-parmed complex.prmtop
-changeLJSingleType @5724 0.3019 0.047
-changeLJSingleType @4843 0.3019 0.047
+Parmed reads the prmtop file and allows you to modify the parameters defined in it. The parmed command we are going to use is: ''changeLJSingleType'', which changes the LJ type for all of the atoms that share the LJ type, or ''addLJType'', for more selective modifications. We will change their LJ parameters to those of a GAFF2 type alcohol that has ''0.3019 0.047'' values for the LJ type.
+Either if we are using ''changeLJSingleType'' or ''addLJType'', we need to ensure that we are changing the right atom index. To do so, we need to use ''printDetails'' or ''printLJTypes'' that print the information on the selected atom index.
+<code>
+$ parmed complex.prmtop
+changeLJSingleType :WAT@H1 0.3019 0.047
+changeLJSingleType :*@HO 0.3019 0.047
 outparm complex_LJ_mod.prmtop
 quit