View source for Mod:Hunt Research Group/Chemshell:Chemshell: QM/MM Single Point and optimisation

To run a QM/MM computation we need do modify our .py file and define a new theory variable for the QM computations required. The .pun file will be the same we used for the MM computations. However, the molecules that are included in the QM region must not be 'counted' in the .ff file: if there are 3 molecules of water in the system and one is included in the QM region, the number of water molecules reported in the field file will be 2 (for DL_POLY: nummols=2).

The input files needed for this tutorial can be found at the following link: 

<pre> 
https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/Chemshell:Chemshell:_QM/MM_Single_Point_and_optimisation_inputs
</pre>

Most of the .py is the same as the one we used for the MM single point calculation. The first new line we can find is: 

<pre>
qm = NWChem(method="hf", basis="3-21g", charge=-1)
</pre>

Here we have defined a new theory variable 'qm' in which we store information relative to the QM computations we want to perform. The class of the object determines the program that will be used for the QM computations. For each program chosen, we will have access to different arguments and argument values.
The first keyword 'method' defines the QC method we want to use; for the moment we will use Hartree-Fock (value "hf"). The B3LYP method can be called by assigning the value "b3lyp" to the method keyword argument and add the keyword argument "functional" (look at example below).

With the next keyword "basis" we set the basis set for the computation. In case you wanted to use a basis set that is not included in NWChem (or include an atom for which parameters are not available), you can import the basis set from an external file by adding the keyword argument "basis", e.g.:

<pre>
qm = NWChem(basis='(name of basis file).basis',  method="dft", functional='b3lyp', charge=1, memory=380000, grid='ultrafine')
</pre>

In this example we have defined two optional arguments "grid='ultrafine'" and "memory=380000"; the former provides an additional specification for the input of NWChem, while the latter sets the memory usage for the computation. This last argument is required when submitting the job to the HPC. 
The last keyword in this line is "charge" which sets the total charge of the QM region. 

The MM theory line remains almost the same as for the MM single point: we only remove the fragment keyword as the initial configuration will be later defined. 
Next, we find the line that defines the various parameters of the QM/MM computations:

<pre>
qmmm = QMMM(frag=cluster, qm_region=[42,67,68,69], qm=qm, mm=mm, embedding='electrostatic', coupling='covalent')
</pre>

The initial configuration is defined here with the frag keyword. The "qm_region" keyword lets us define the atoms we want included in the QM region (remember that indices start from 0 not 1!). With the "qm" and "mm" keywords we set the variables from which to take the information relative to the QM and MM part of the computations we want performed. "embedding" and "coupling" keywords are less important and for the moment you can use the provided values for your computations.

Moving down, we find a block of print statements. These have been included only to make the output include more information and are not necessary. 

As with the MM computation, the file ends with the command setting the computation type. Try running the single point after having substituted the names of your input files where necessary.

Your output should look something like: 

<pre>

 **********************************************************************
 QM energy          =         -532.960994138396 a.u.
 MM energy          =            0.020349646086 a.u.
 QM/MM total energy =         -532.940644492311 a.u.
 **********************************************************************

 Peak memory used by procedure chemsh.interfaces.qmmm.run: 149.508 MB

 ----------------------------------------------------------------------
 Final SP energy    =         -532.940644492311 a.u.
 ----------------------------------------------------------------------
</pre>

Now uncomment the opt lines at the end of the .py file and run a QM/MM optimisation. If you want to save the output that you get on your terminal, launch chemshell like this:

<pre>
chemsh <name of .py file>.py > out
</pre>

Also in this case the optimisation is not supposed to converge with the maximum number of cycles set. Now, if you look into your directory a number of new files has been generated: the files whose name starts with _nwchem are the inputs and outputs of NWChem and are useful for debugging: when a problem arises in one of the involved programs you can directly check these files to find the error (the same goes for dl_poly inputs and outputs).