molUP Tutorial 1 – Start a QM/MM study

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.

1. Launch VMD;
2. Go to “File” > “New molecule…” > “Browse…” and choose the parameters file (.prmtop);
3. Choose the file type “AMBER7 Parm”;
4. Click on “Load”
5. Click on “Browser…” and choose the coordinates file (.inpcrd or .rst);
6. Choose the file type “AMBER7 Restart”;
7. Click on “Load”
8. 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.

1. Open the “Graphical Representations” window under the “Graphics” > “Representations” menu;
2. Change the “all” representation to the desired selection of atoms, such as:
   same residue as (all within 5 of protein)
3. Select the molecule on the “VMD main” window;
4. Go to “File” > “Save coordinates…”;
5. Under the “Selected atoms:” option, choose the atom selection that you defined;
6. Choose the file type xbgf;
7. Click on “Save…”;
8. 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

1. Launch VMD;
2. Start molUP through the “Extensions” > “VMD Store” > “molUP” menu.
3. On the “VMD main” window, go to “File” > “New molecule” > “Browse..”, and choose the xbgf file generated in the previous step;
4. 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);
5. 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:

1. On the molUP window, navigate to the “Model” > “Layer” tab;
2. 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;
3. A combobox allows you to choose the layer (H, M, or L);
4. 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.

1. On the molUP window, navigate to the “Input” tab;
2. (Optional) Choose a title for your job under the “Job Title” field;
3. 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.
4. 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;
5. Click on the “Other information (Parameters, Modredundant…)” menu to “Load parameters from a PRMTOP file”;
6. Choose the parameters file (the same used at the beginning of this tutorial);
7. 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.

1. On the molUP window, go to “File” > “Save” > “Save” and choose a name to save the Gaussian input file, e.g. “optimization.com”;
2. Quit VMD.
3. Launch the calculation, e.g.:
   g09 optimization.com

—

A new tutorial will be published on the next week.