Gaussian software (commercial license) is used to study several types of molecular problems, namely catalytic mechanisms. However, if we are interested in studying an enzymatic mechanism, the number of atoms of the system is prohibitive in terms of using Quantum Mechanics (QM) for the entire system. Therefore, the ONIOM QM/MM approach is commonly used to treat this type of systems, where the molecules directly involved in the catalysis are considered using a QM approach, and the remaining atoms are calculated using Molecular Mechanics (MM).
The QM/MM approach has the inconvenience of requiring parameters for the MM calculations. Consequently, the generation of parameters is required before starting the calculations. This task can be accomplished using the AmberTools package (free) to generate a set of two files:
- Parameters file (.prmtop)
- Coordinates file (.inpcrd or .rst)
So, these two files are required before starting to build any QM/MM study. The parameters file must be similar to the following one:
%VERSION VERSION_STAMP = V0001.000 DATE = 02/16/18 09:58:14
%FLAG TITLE
%FORMAT(20a4)
default_name
%FLAG POINTERS
%FORMAT(10I8)
13129 18 6649 6620 14541 8972 28436 22349 0 0
71243 927 6620 8972 22349 77 168 58 50 1
0 0 0 0 0 0 0 0 53 0
0
%FLAG ATOM_NAME
%FORMAT(20a4)
N H1 H2 H3 CA HA CB HB2 HB3 CG HG CD1 HD11HD12HD13CD2 HD21HD22HD23C
O N H CA HA CB HB2 HB3 CG HG2 HG3 CD HD2 HD3 CE HE2 HE3 NZ HZ1 HZ2
HZ3 C O N H CA HA CB HB2 HB3 CG HG2 HG3 CD HD2 HD3 NE HE CZ NH1
HH11HH12NH2 HH21HH22C O N H CA HA CB HB2 HB3 CG HG2 HG3 CD OE1 OE2
--- TRUNCATED FILE ---
8.50000000E-01 8.50000000E-01 8.50000000E-01 8.50000000E-01 8.50000000E-01
8.50000000E-01 8.50000000E-01 8.50000000E-01 8.50000000E-01 8.50000000E-01
8.50000000E-01 8.50000000E-01 8.50000000E-01 8.50000000E-01 8.50000000E-01
8.50000000E-01 8.50000000E-01 8.50000000E-01 8.50000000E-01
%FLAG IPOL
%FORMAT(1I8)
0
The coordinates file could be a PDB file or a coordinates file generated by AmberTools, in that case, it should be similar to the following one:
default_name
13129
33.1680000 113.2160000 67.7670000 33.4743499 114.1782406 67.7855420
33.9665432 112.6180098 67.9245970 32.3468423 113.1255172 68.3480443
32.7070000 112.9020000 66.3800000 33.3041921 112.0868937 65.9712672
--- TRUNCATED FILE ---
15.5422000 94.3770000 51.7510000 14.3450120 95.3036270 51.7510000
4.6010000 93.4390000 51.0530000 5.5582000 93.4390000 51.0530000
4.3610120 94.3656270 51.0530000
ATTENTION: Before you start, please ensure that you have the parameters file and the coordinates file. You also need to install VMD (1.9.3 or later) (free and available here) and the molUP (free and available here)
1. Prepare the structure
At this moment, you have the coordinates of all the atoms of your system, and also all the required parameters for the MM calculations. However, sometimes there are some atoms that you do not want to include the QM/MM calculation, e.g., the entire water box or the entire protein. In that case, we need to remove some of those molecules before we start to build the QM/MM model.
1.1. Open the molecular system using VMD
In this step, you need to open the molecular system on VMD using the parameters and coordinates files generated with AmberTools.
- Launch VMD;
- Go to “File” > “New molecule…” > “Browse…” and choose the parameters file (.prmtop);
- Choose the file type “AMBER7 Parm”;
- Click on “Load”
- Click on “Browser…” and choose the coordinates file (.inpcrd or .rst);
- Choose the file type “AMBER7 Restart”;
- Click on “Load”
- Close the “Molecule File Browser” window.
1.2. Select the atoms to be included in the QM/MM model
Now, you have to select which atoms must be included in the QM/MM model.
- Open the “Graphical Representations” window under the “Graphics” > “Representations” menu;
- Change the “all” representation to the desired selection of atoms, such as:
same residue as (all within 5 of protein)
- Select the molecule on the “VMD main” window;
- Go to “File” > “Save coordinates…”;
- Under the “Selected atoms:” option, choose the atom selection that you defined;
- Choose the file type xbgf;
- Click on “Save…”;
- Quit VMD.
After this procedure, you generated a .xbgffile containing only the atoms that must be included in the QM/MM model.
1.3. Prepare the QM/MM model using the molUP extension for VMD
- Launch VMD;
- Start molUP through the “Extensions” > “VMD Store” > “molUP” menu.
- On the “VMD main” window, go to “File” > “New molecule” > “Browse..”, and choose the xbgf file generated in the previous step;
- Click on “Load” and wait a few seconds (it could take a bit longer depending on the number of atoms that molUP have to process);
- Close the “Molecule File Browser” window.
1.3.1. Setting ONIOM layers
One of the most important tasks at the beginning of any QM/MM study is to chose which atoms must be included in the high-layer (HL) which is the layer that will be considered using the QM approach. In this example, a two-layer system, containing one HL and one low-layer (LL), will be used but molUP and Gaussian also allow you to define a medium-layer (ML) between the HL and the LL.
The following steps describe how molUP can help you choosing the atoms that belong to each layer:
- On the molUP window, navigate to the “Model” > “Layer” tab;
- You can use VMD atom selections, selections from the list of atoms, and pick atoms to choose the atoms that must belong to a certain layer;
- A combobox allows you to choose the layer (H, M, or L);
- The “Apply” button applies the currently selected layer to the selected atoms.
IMPORTANT: Make sure that all the atoms have a layer assigned in order to make the ONIOM QM/MM calculation possible.
1.3.2. Setting fixed atoms
In some cases, it is important to fix the position of certain atoms during the optimization process. For example, in the case of this molecular system, it is important to keep fixed the external coat of water molecules to mimic the solvent surrounds and keep them nearby the protein. In some cases, you can also fix all the atoms that are at least 15 Å from the HL. In this last case, you could use the following atom selection to fix the position of those atoms:
all and not (same residue as (all within 15 of altloc H))
The altloc H
selection corresponds to all the atoms included in the HL. You can also use altloc M
or altloc L
to select atoms included in the ML or LL, respectively.
The process to select the atoms under this tab is the same used on the previous step (1.3.1).
1.3.3. Choose the calculation keywords
This is the last step of this procedure to generate an ONIOM QM/MM input file for Gaussian. Now, the calculations keywords have to be chosen, the charge and spin multiplicity of each ONIOM layer must be set, and the parameters need to be loaded.
- On the molUP window, navigate to the “Input” tab;
- (Optional) Choose a title for your job under the “Job Title” field;
- Fill the “Calculations” section with all the Gaussian calculation keywords; molUP has a preset of calculation templates that can be easily chosen from the “Calculations” menu.
You can also set up your own calculation templates and save them to use later in other calculations. - In the “Charge and Multiplicity” section, molUP provides an automatic tool to calculate the charge of each layer. You can also adjust the charge and spin multiplicity of each ONIOM layer;
- Click on the “Other information (Parameters, Modredundant…)” menu to “Load parameters from a PRMTOP file”;
- Choose the parameters file (the same used at the beginning of this tutorial);
- After the parameters being loaded, add two blank lines before the parameters section. (This is mandatory and cannot be automatically predicted because it depends on the type of calculation).
1.3.4. Save the Gaussian input file
Finally, you are able to save the structure as a Gaussian input file that is ready to be used to start a calculation.
- On the molUP window, go to “File” > “Save” > “Save” and choose a name to save the Gaussian input file, e.g. “optimization.com”;
- Quit VMD.
- Launch the calculation, e.g.:
g09 optimization.com
—
A new tutorial will be published on the next week.
Consequently, the generation of parameters is required before starting the calculations. Could you please and please give a brief guidelines on how to generate the Parameters file (.prmtop) and Coordinates file (.inpcrd or .rst) using the AmberTools
Thank you
LOUIS
Hello Biosim,
I want to know how I can analyze covalent bonds or hydrogen bond break/formation in a protein-ligand complex using Molup