View source for Mod:Hunt Research Group/Chemshell: MM Single Point computation

In this tutorial we are going to start learning about the basic commands of ChemShell and the necessary inputs for QM/MM computations. First of all we need to load the ChemShell module after logging in to cx1:
<pre>
module load chemshell/ . . .
</pre>

After that you need to create a new directory on the hpc home page and upload the input files there. All of the files you will need can be found here: [https://wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:Hunt_Research_Group/Chemshell:Chemshell:_MM_Single_Point_computation_inputs MM Single Point Inputs]


The first computation we are going to carry out with ChemShell is a simple MM (molecular mechanics) single point calculation. In other words we are computing the energy of a given atomic configuration, describing the interactions with the MM forcefield provided as one of the necessary inputs. The only output that will be generated is the energy of the configuration. Thus, this is neither a QM/MM computation nor an optimisation. Still, we will get to know the commands and files we will be using for all our QM/MM optimisations. 

The files that you have just uploaded comprise: 

# A .pun file, which contains the (initial) atomic positions (sorted by molecule type), and the charge carried by each atom.
# A .ff file, specifying the force field functional form and coefficients; the format of the .ff file will change accordingly to the package that we use for the MM computations: in our case we are using DL_POLY, therefore the .ff file will have exactly the same format of the DL_POLY FIELD input file.
# A .py file, this is the python module through which we specify the calculation that we want ChemShell to perform and define the value of all of the parameters involved. 

Let's walk through the .py file to understand the meaning of the various commands it contains. Where necessary, fill in the name you chose for the various inputs. In what follows, is useful to remember that ChemShell is written in Python language. After the commented title line, the first command imports all module contained in ChemShell. Remember always to include this line, otherwise your runscript will be useless. 

Moving down we find: 
<pre>
reline = Fragment(coords="<name of the pun file>.pun", connmode=None)
</pre>
Here the new variable 'reline' is initialised as an object of the class 'Fragment'. Fragment objects are the ones containing the atomic position and charges, together with other info (connectivity, labels, ..). We will always need to define a Fragment object to run a computation in ChemShell, either as the starting configuration for an optimisation, or as the investigated configuration in a single point computation. The keyword argument 'coords'
defines the file from which to extract the atomic coordinates. The 'connmode' argument determines the connectivity mode; for the moment you can ignore this command and always leave 'None' as the argument's value. 

The next line in the .py script we find is: 
<pre>
mm = DL_POLY(frag=reline, ff='<name of the field file>.ff', temperature=0.0)
</pre>
In this line we are defining the MM 'theory' for our calculation, i.e. the package, the forcefield, the initial configuration and all other relevant parameters for the MM part of the computation we want to perform. A number of optional keyword arguments is available for this object; let's focus on the required arguments for the moment. 
The object class used in the definition, determines the package that will be used for the MM computations and consequently the accessible keyword arguments in the following definition. In our case we are going to use DL_POLY. 
The 'frag' argument, determines which fragment object should be used as the starting point (in our case the only point) for our computation. Obviously, you need to define the fragment object of interest BEFORE using it as an argument in a definition/function. The subsequent argument 'ff' determines the file from which the forcefield will be imported for this computation. As already stated above, this file will need to have a suitable format depending on the package that we want to use for the MM computation. Finally, the 'temperature' argument trivially determines the temperature of the system during the computation. In most cases, you will be optimising a given structure neglecting any thermal movement (atomic kinetic energy = 0) and the assigned temperature value will be '0.0'

The last thing left to do before running, is to tell ChemShell the type of computation we want to perform. The next line in the .py script does just that:
<pre>
sp = SP(theory=mm)
</pre>
Here we initialise a new variable of a different object type depending on the required computation. In our case, we want to perform a Single Point Calculation, hence we use the object type SP(). SP object only have one necessary argument: 'theory'; with this argument we tell ChemShell the theory object from which to take all of the necessary information for the computation. We will talk about other optional keywords argument for this command later on. 
Now we have all we need to perform a MM single point computation. The last line in the script:

<pre>
sp.run()
</pre>

is just telling ChemShell to run the computation defined within the 'sp' object. 

Close the .py file and run this single point computation with the following command: 

<pre>
chemsh <module name>.py
</pre>

At the end of the resulting output, you should get something like:

<pre>
>>> Running single point calculation

>>> DL_POLY run successfully.
Peak memory used by procedure chemsh.base.run.run: 90.617 MB
Elapsed time of procedure chemsh.interfaces.mm.run: 0.168 sec

----------------------------------------------------------------------
Final SP energy    =           -0.195569750490 a.u.
----------------------------------------------------------------------

Elapsed time of procedure chemsh.tasks.sp.run: 0.169 sec
</pre>

Now that we have successfully computed the MM energy of our initial configuration, let's run an MM optimisation!

To run an optimisation with the same MM 'theory' (same package, forcefield, initial configuration, temperature, ecc..) as the single point just performed we just need to modify the last two lines in the script, changing the computation type from SP (Single Point) to Opt (Optimisation). To make things even easier for you, the two new lines required are already reported at the end of the script, commented out. Just comment the previous last two lines and uncomment the 'opt' lines. Save the changes to the script and run Chemshell again with the same command as before: 

<pre>
chemsh <module name>.py
</pre>

In the 'opt' variable we use the keyword argument 'maxcycle' to set the maximum number of optimisation cycles.
Don't worry if the optimisation doesn't converge. It is not supposed to!

The outputs that are generated by this computations are:
# _chemsh_run.log: a list of the times at which chemshell has been launched in that directory reporting the specific command used each time. 
#CONFIG, FIELD, CONTROL, OUTPUT, REVCON, REVIVE and STATIS files are inputs and outputs of DL_POLY and will always be generated when DL_POLY is used. 

The following outputs appear only if an optimisation is ran:

#"path" files report information on the optimisation trajectory. In an other tutorial we will show how to visualise optimisations in VMD using the "path.xyz" file. 
#_dl_find.pun is the last configuration accepted during the optimisation. You can use this as an initial configuration to restart a crashed or incomplete optimisation.