22 May 2012

159. PES scanning of methanol bonds, angles, torsion using nwchem, nwgeom and python

NOTE: there's something dodgy with the potential/bond length plots -- they optimal bond lengths are way too long. I'll leave this post up here anyway, but be WARNED.

This is more of an overview of an idea of how to do it together with some explicit examples. This is more of a sketch than a step-by-step account.

Today's molecule is Methanol.

0. ecce and nwgeom/python
You need to set PYTHONPATH to /opt/nwchem/nwchem-6.1/contrib/python in order for nwchem to find the nwgeom module. That's easy enough on a local system since ~/.bashrc is read -- but it won't read ~/.bashrc on remote systems. For this you need to edit your CONFIG.<machine> files (in ~/.ECCE on your main node)  -- add
setup {
     setenv PYTHONPATH /opt/nwchem/nwchem-6.1/contrib/python

NWChemEnvironment {
          PYTHONPATH /opt/nwchem/nwchem-6.1/contrib/python

1. Draw the molecule,
Draw the Carbon, then the oxygen, then the protons on the carbon, then the protons on the oxygen. Basically, draw the backbone first, then add protons by hand.  Turn it into a residue-based system (under 'build') and optimise the structure using e.g. RHF/6-31G* (this was written with MM/FF parametrisation in mind -- for simple scanning, just do whatever you want )

2. Calculate the partial charges (rhf/6-31g*).  Can skip this
You can e.g. constrain the methyl groups, or force all the methyl protons to be equal. It's a bit of a soft science, really. After the calc has finished, assign charges. I can't claim to understand which method is better (RESP, CRESP, CRESP2 etc.).

 (this was written with MM/FF parametrisation in mind -- for simple scanning, skip this step)

Here's CRESP2 (some variability...):
C -0.721
O -0.368
H 0.240
H 0.240
H 0.240
H 0.368

Also, assign (atom) Types (CT, OH, HC, HC, HC, HO) -- this is done by hand. Pick atom table, select residue view (or something similar), and fill in the Type column.

Then click on Tools, check Residue table, click on the menu-looking icon in the residue view, and select write fragment. Make sure you put the fragment file in a place where ecce and nwchem will find it (e.g. amber_u). For some reason I can't get ecce to actually change the name of my residues, so edit the fragment file by hand and change all instance of UNK to the same name as the fragment file, e.g. TST if you called it TST.frg.

3. Write down the bonds, angles and dihedral angles.
H-C 1.087
C-O 1.400
O-H 0.946

H-O-C 109.467
O-C-H 112.039
H-C-H 108.682

H-C-O-H -61.229

When you scan the parameters using python you want to be able to
1) see if the lowest energy conf make sense and
2) not deviate too much from the ideal angle/bond/torsion.

Things get weird if you do.

4. Try to determine bond strength
This is best done outside of ecce, and you really should have compiled nwchem with python-support to make this easier.

Copy the input file you used for the ESP calc. Call it bonds.nw. Remove anything about esp and all task directives, then add:

Technically I think the bond strengths only really need two data points unless you want to fit the Lennard-Jones equation to them, but it certainly looks neater getting the full behaviour.

Then run using
mpirun -n 2 nwchem bonds.nw | tee bond12.nw

You can also do it in ecce -- but the plotting will have to be done by hand (I open the out file with vim, select the energy/atom-distance columns, :w angle.dat and then plot with gnuplot)

Do the same for atom pair 1-6 and 2-3. Make sure to pull the atoms far enough apart that the energy tails out (the 1-6 pair in the figure needs to be separated more)

5. Angles
As you'll discover, it's not just a matter of throwing in random numbers and scanning -- if you don't collect enough points, or if your first point is far away from the optimal angle, the data will look very odd.

6. Torsion


  1. Lindqvist,
    thank you for the informative info above...I would like to run dihedral relaxed PES scan (similar to what is run in Gaussian:

    A relaxed PES scan (with geometry optimization at each point) is requested with the Opt keyword.

    does one only need to subsitute

    'scf optimize' in place of where scan is used in the following route (input line)



    1. change task_energy to task_optimize
      See here for more information:

  2. In this line :
    " for i in range(0,len(results)):
    print results[i][0][0],results[i][1] "
    How do you know that results[i][0][0] will give our required bond length and results[i][1] will give energy?
    I am new to Python too. Can you suggest some literature to read about or send the links here?

    1. See https://svn.pnl.gov/svn/nwchem/trunk/contrib/python/nwgeom.py -- it's a list of tuples.
      Best way to learn is to create a project, and use google whenever you get stuck.

  3. How to do the same for a mcscf calculation, i.e. specify the active space and other stuff?

    1. The same? It's a completely different type of computation.
      See http://www.nwchem-sw.org/index.php/Release66:Multiconfiguration_SCF