10 May 2012

147. Oniom in nwchem -- with a little bit of help from gromacs and openbabel

Example -- I want to do explicit solvent modelling of methanol in water. This is obviously an articifical approach, but generally applicable.

This is a rough approach to doing oniom calculations using nwchem 6.0 -- this is a technical description, not a how-to when it comes to the science.

1. Pre-optimisation
Draw methanol and set up a simple calc using e.g. ecce to pre-optimise the structure with e.g. an implicit solvent model. Here's nwch.nw:

scratch_dir /scratch
Title "pre-oniom"
Start  pre-oniom
charge 0
geometry autosym units angstrom
 C     0.00000     0.00000     0.00000
 H     -0.675500     -0.675500     0.675500
 H     0.675500     -0.675500     -0.675500
 H     -0.675500     0.675500     -0.675500
 O     0.866025     0.866025     0.866025
 H     1.51843     0.213620     1.51843
ecce_print /home/me/jobs/jobs/testing/old/pre-oniom/ecce.out
basis "ao basis" spherical print
  H library "6-31+G*"
  O library "6-31+G*"
  C library "6-31+G*"
  mult 1
  XC b3lyp
task dft optimize
task dft freq

2. Solvation using gromacs
Take the output, nwch.nwout, and use babel to export the optimised structure

babel -inwo nwch.nwout -oxyz molecule.xyz

The next few steps require gromacs:
editconf -f molecule.xyz -o molecule.gro -box 2 2 2
genbox -cp molecule.gro -cs spc216.gro -o solvated.gro
babel -igro solvated.gro -oxyz solvated.xyz

tail -n +3 solvated.xyz > oniom.nw

3. Putting it all together

The only 'trick' is how to define what part of the input belongs to the high level section and what belongs to the low level section -- for a solvation case like this, where whole molecules are treated by one method or another, use model. The last atom to be part of the high level section is the last of the six atoms in methanol, so the keyword is model 6. Atoms 7 to infinity are thus part of the sto-3g part, and atoms 1-6 part of the 6-31g* part.

memory stack 600 mb heap 200 mb global 800 mb

scratch_dir /scratch
Title "oniom"
Start oniom
charge 0

geometry units angstrom
symmetry c1
C         10.12000       10.35000       10.00000
H          9.60000       10.72000       10.89000
H          9.60000       10.72000        9.11000
H         11.14000       10.74000       10.00000
H          0.46000        2.39000        6.86000
H          1.16000        1.73000        8.18000

basis sto-3g spherical
 * library sto-3g

basis 6-31g spherical
 * library 6-31g

    model 6
    low  dft basis sto-3g input "dft\; xc\; end"
    high dft basis "6-31g" input "dft\; xc\; end"

basis "ao basis" spherical print
  H library "6-31+G*"
  O library "6-31+G*"
  C library "6-31+G*"

task oniom

This takes forever to run, but run it does. The memory statement is important -- if the global memory is too small it will crash. Also, be aware that the amount of memory specified is per instance -- if you launch with mpirun -n 3, multiply by 3 to get the amount of memory that needs to be present (here, 2.4 GB physical RAM).

No comments:

Post a Comment