Talk:Mod:Hunt Research Group/QMMM SP Aqeuous Cu(II)

QM/MM Single Point Energy Calculation of Cu(II) and 2 Solvation Shells

In this example, the first solvation shell is modelled at the B3LYP 6-311G level (via Gaussian) while the second solvation shell is modelled using a classical pairwise potential (via DL_POLY). ChemShell employs an additive scheme to calculate the total QM/MM energy i.e. the total energy is the sum of the MM energy of the MM region and the QM energy of the QM region. In addition, link atoms are placed along the axis of bonds at the QM-MM boundary. More details can be found here.

Create an input file called energy.chm with the following content:

energy energy=e coords=cu_18water.pun \
  theory= hybrid : { coupling=shift
                     qm_region = {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22}
                     qm_theory=gaussian : { basis=6-311g g98_mem=80000000 charge=2 mult=2 hamiltonian=b3lyp}
                     mm_theory=dl_poly : mm_defs=ff.dat
                   }

If you are modifying a opt.chm file generated in MolCluster, these are the only lines that are needed, so remove the other lines in the file such as "atom types", "groups", and "active atoms".

The co-ordinates of the system are read from cu_18water.pun (created here). The first 22 atoms (Cu(II) + 7 waters in the first solvation shell) are modelled at the B3LYP 6-311G level. The rest of the system is modelled using the force field defined in ff.dat (see here for details). The QM and MM regions are coupled via the charge shift scheme.

The calculation should be submitted to the PBS queue to avoid overloading the login shell and inconveniencing other cx1 users (details here). Don't forget to load the DL_POLY and Gaussian modules in your submit script!

The energy can be found at the end of the output (energy.out):

-----------------------------------------------------------------------------------
Energy (     hybrid):                                           -2175.199009 (a.u.)
-----------------------------------------------------------------------------------
ChemShell exiting code 0