NWChem

When in doubt, check the Documentation, and test!

Running NWChem

First, login to a system (e.g. circe) with nwchem installed.

Next, load the NWChem module and set up a working dir. <source lang="bash"> module load apps/nwchem/6.1.1 mkdir nw-test cd nw-test </source>

Then create a molecule (you can download an sdf file from the PDB ligand structure database listed in the course reference material (Search by molecule name -> Download Links (on left panel) -> Ideal Molfile).

Next, paste the molecule into a file, e.g. <source lang="bash"> cat >start.sdf <<. FOR.sdf

 -ISIS-            3D

 4  3  0  0  0  0  0  0  0  0  0 V2000
   0.6070    0.0000    0.0000 C   0  0  0  0  0
  -0.6000    0.0000    0.0000 O   0  0  0  0  0
   1.1470    0.9350    0.0020 H   0  0  0  0  0
   1.1470   -0.9350    0.0020 H   0  0  0  0  0
 1  2  2  0  0  0
 1  3  1  0  0  0
 1  4  1  0  0  0

M END $$$$ . </source>

You'll also need to paste one of the templates below (they work without changes) into an NWChem control file called en_scf.nw. Use the 'cat' command to paste into the file as above.

Finally, convert start.sdf to an xyz format for NWChem and run nwchem. <source lang="bash"> babel start.sdf start.xyz nwchem en_scf.nw >en_scf.log </source>

All your output will be in en_scf.log (scp/download it to your local system). In the future, use the batch queue system with <source lang="bash"> cat >en_job.sh <<.

1. $ -N a-test-job
2. $ -cwd
3. $ -o job-$JOB_ID.out
4. $ -e job-$JOB_ID.err
5. $ -l pcpus=4,h_rt=06:00:00

module load apps/nwchem/6.1.1 cd $HOME/nw-test mpirun `which nwchem` en_scf.nw >en_scf.log .

qsub en_job.sh

qstat -u $USER -t </source>

This runs on 4 cpus in parallel, so it will be faster than running (as before) directly on the login node. Of course, you'll also want to change some options in the run file eventually (especially the basis set). You can also manually edit the xyz file and scan the energy vs some coordinate. See the templates below for more ideas.

Note that you can also edit the files on your system and transfer, or edit locally on circe using your favorite linux editor. Nano is easy to use. <source lang="bash"> nano en_scf.nw </source>

Energy Calculation

A simple SCF energy calculation on the input file start.xyz can be done with:

start en_scf
title "Pople SCF energy"

geometry units angstrom
  load start.xyz
end
basis
#   * library cc-pvdz
#   * library 6-31G*
   * library 3-21G
end

task scf energy

Geometry Optimization

Running a geometry optimization just requires replacing the task directive. Here, we've also added a block of parameters controlling the minimization algorithm.

start opt_scf
title "Pople SCF geometry optimization"

geometry units angstrom
  load start.xyz
end
basis spherical
#   * library cc-pvdz
#   * library 6-31G*
   * library 3-21G
end

driver
   loose
   maxiter 150
   xyz opt
end
task scf optimize

DFT

You can also use DFT rather than HF (SCF) by changing the appropriate keywords and adding a DFT block to specify its functional:

start opt_dft

title "Pople DFT geometry optimization"

geometry units angstrom
  load start.xyz
end
basis spherical
#  * library cc-pvdz
#  * library 6-31G*
  * library 3-21G
end
dft
  xc xpbe96 1.0 \
  pw91lda local 1.0 \
  cpbe96 nonlocal 1.0
  direct
  iterations 150
  grid fine
  disp vdw 2
end

driver
  loose
  maxiter 150
  xyz opt
end
task dft optimize

The "xyz opt" causes NWChem to write out coordinate files, e.g. opt-001.xyz, ... as the optimization is proceeding.

Specifying Charge and Spin

Of course, we need to be able to specify the charge and total spin for most systems as well. Here's an input example appropriate for Fe3+ (5 spin-up electrons making up a half-filled d-shell):

start en_scf
title "Pople SCF energy"

geometry start units angstrom
  load start.xyz
end
basis spherical
#   Fe  library "Ahlrichs pVDZ" # better basis
    * library 3-21G
end

set geometry start
charge 3
scf
  sextet
end

task scf energy

Vibrational Frequencies

At a minimum, the second derivatives of the potential energy surface with respect to the nuclear coordinates make up a 3Nx3N matrix. The eigenvectors represent vibrational 'modes,' and the eigenvalues represent their force constants. Analyzing each mode as a harmonic oscillator leads to a set of vibrational frequencies, which give thermochemical information on the molecule as well as IR spectroscopic information.

From a minimized set of coordinates, all you should require is to add the appropriate task (and optionally a block of parameters)

task scf frequencies

Orbital Analysis

Molecular orbital analysis is relatively easy once the HF equations have been solved, since the linear combinations of atomic orbitals that make up each of the molecular orbitals are stored in a file (jobname).movecs. To plot the total density and the contents of orbitals numbered 54 and 55 (from checking through the HF output), use:

dplot
  title dens
  vectors en_scf.movecs
  LimitXYZ
  -5.0 5.0 30
  -5.0 5.0 30
  -5.0 5.0 30
  spin total
  gaussian
  output dens.cube
end
task dplot

dplot
  title homo
  vectors en_scf.movecs
  LimitXYZ
  -5.0 5.0 30
  -5.0 5.0 30
  -5.0 5.0 30
  spin total
  orbitals view; 1; 54
  gaussian
  output homo.cube
end
task dplot

dplot
  Title lumo
  vectors en_scf.movecs
  LimitXYZ
  -5.0 5.0 30
  -5.0 5.0 30
  -5.0 5.0 30
  spin total
  orbitals view; 1; 55
  gaussian
  output lumo.cube
end
task dplot

The Gaussian-format cube files can be opened with pymol and visualized by typing, e.g.

isomesh mesh1, dens, 0.1, all, carve=1.6
isomesh mesh2, lumo, 0.01, all, carve=1.6
isomesh mesh3, lumo, -0.01, all, carve=1.6

to create isosurfaces at the values 0.1 for the density and <math>\pm 0.01</math> for the LUMO's wavefunction.

Adam Hogan noted that (according to [1]) these orbitals are still associated with HF theory. To compute MP2-level quantities (natural orbitals), you instead use a matrix of derivatives of the orbitals to calculate matrix elements of an operator <math>\langle \psi | O | \psi \rangle</math>. To get those, you have to run task mp2 gradient, which creates (jobname).mp2nos. You then replace vectors en_scf.movecs with vectors en_scf.mp2nos inside the dens blocks above to get the high-level densities. Here's an example of using mo2nos values to get the electrostatic potential:

mp2
  tight 
  freeze atomic 
end

task MP2 gradient

set "esp:input vectors" en_scf.mp2nos

esp
  recalculate 
  probe   0.07
  range   0.3
  factor  1
  spacing 0.02
end
task esp