E. Hybrid Potentials

Calculations with QC methods rapidly become too expensive as the size of a system increases. Nevertheless, for certain problems, such as the simulation of enzymatic or solution-phase reactions, it is essential to have a QC description of at least a small portion of the system. An appropriate solution to this dilemma, in some circumstances, is to use a hybrid potential in which different parts of a system are treated with different potentials. Although there are a large variety of hybrid potentials, pDynamo's employ two regions, one in which the atoms are treated with a QC model, and the other in which they are represented with a MM method.

pDynamo programs that employ hybrid potentials are very similar to their pure QC and MM counterparts. The major differences are that MM, QC and NB models all have to be defined for the system (in that order!) and that the atoms to be included in the QC region need to be specified in the method call to DefineQCModel. The pDynamo distribution includes two hybrid potential examples, one in which different molecules are in the QC and MM regions and the other in which a single molecule is divided between the two regions. The first is Example8.py which calculates the energy of a water dimer:

# . Define the MM, NB and QC models.

mmModel = MMModelOPLS.WithParameterSet ( "bookSmallExamples" )

nbModel = NBModelFull.WithDefaults ( )

qcModel = QCModelMNDO.WithDefaults ( )

# . Define the molecule.

molecule = ImportSystem ( os.path.join ( molPath, "waterDimer_cs.mol" ) )

# . Define the selection for the first molecule.

firstWater = Selection.FromIterable ( [ 0, 1, 2 ] )

# . Define the energy model.

molecule.DefineMMModel ( mmModel )

molecule.DefineQCModel ( qcModel, qcSelection = firstWater )

molecule.DefineNBModel ( nbModel )

molecule.Summary ( )

# . Calculate an energy.

molecule.Energy ( doGradients = True )

The second example is Example9.py which determines the energy of a bALA molecule for which the side-chain methyl group is in the QC region and the remainder of the molecule in the MM region:

# . Define the MM, NB and QC models.

mmModel = MMModelOPLS.WithParameterSet ( "bookSmallExamples" )

nbModel = NBModelFull.WithDefaults ( )

qcModel = QCModelMNDO.WithDefaults ( )

# . Define the molecule.

molecule = ImportSystem ( os.path.join ( molPath, "bala_c7eq.mol" ) )

# . Define the selection for the first molecule.

methylGroup = Selection.FromIterable ( [ 10, 11, 12, 13 ] )

# . Define the energy model.

molecule.DefineMMModel ( mmModel )

molecule.DefineQCModel ( qcModel, qcSelection = methylGroup )

molecule.DefineNBModel ( nbModel )

molecule.Summary ( )

# . Calculate an energy.

molecule.Energy ( doGradients = True )

Exercises

• 1. The interaction energy of a dimer is defined as the difference between the energy of the dimer and the sum of the individual monomer energies. Take a conformation (or several) of the water dimer and compute its interaction energy with various QC (both ab initio and semi-empirical), MM and hybrid potential models. For each hybrid potential, consider the two partitionings in which different monomers are in the QC region. How different are the energies in each case?