Source code for qmmm

"""Module with classes to integrate MM charges into
a QM calculation.

"""
import PsiMod
import re
import os
import math
import physconst
from molutil import *
from driver import *


[docs]class Diffuse(object): def __init__(self, molecule, basisname, ribasisname): self.molecule = molecule self.basisname = basisname self.ribasisname = ribasisname self.basis = None self.ribasis = None self.da = None self.Da = None self.wfn = None def __str__(self): s = ' => Diffuse <=\n\n' s = s + ' ' + str(self.molecule) + '\n' s = s + ' ' + self.basisname + '\n' s = s + ' ' + self.ribasisname + '\n' s = s + '\n' return s
[docs] def fitScf(self): """Function to run scf and fit a system of diffuse charges to resulting density. """ basisChanged = PsiMod.has_option_changed("BASIS") ribasisChanged = PsiMod.has_option_changed("DF_BASIS_SCF") scftypeChanged = PsiMod.has_option_changed("SCF_TYPE") basis = PsiMod.get_option("BASIS") ribasis = PsiMod.get_option("DF_BASIS_SCF") scftype = PsiMod.get_option("SCF_TYPE") PsiMod.print_out(" => Diffuse SCF (Determines Da) <=\n\n") activate(self.molecule) PsiMod.set_global_option("BASIS", self.basisname) PsiMod.set_global_option("DF_BASIS_SCF", self.ribasisname) PsiMod.set_global_option("SCF_TYPE", "DF") energy('scf') PsiMod.print_out("\n") self.fitGeneral() PsiMod.clean() PsiMod.set_global_option("BASIS", basis) PsiMod.set_global_option("DF_BASIS_SCF", ribasis) PsiMod.set_global_option("SCF_TYPE", scftype) if not basisChanged: PsiMod.revoke_option_changed("BASIS") if not ribasisChanged: PsiMod.revoke_option_changed("DF_BASIS_SCF") if not scftypeChanged: PsiMod.revoke_option_changed("SCF_TYPE")
[docs] def fitGeneral(self): """Function to perform a general fit of diffuse charges to wavefunction density. """ PsiMod.print_out(" => Diffuse Charge Fitting (Determines da) <=\n\n") self.wfn = PsiMod.reference_wavefunction() self.Da = self.wfn.Da() self.basis = self.wfn.basisset() parser = PsiMod.Gaussian94BasisSetParser() self.ribasis = PsiMod.BasisSet.construct(parser, self.molecule, "DF_BASIS_SCF") fitter = PsiMod.DFChargeFitter() fitter.setPrimary(self.basis) fitter.setAuxiliary(self.ribasis) fitter.setD(self.Da) self.da = fitter.fit() self.da.scale(2.0)
[docs] def populateExtern(self, extern): # Electronic Part extern.addBasis(self.ribasis, self.da) # Nuclear Part for A in range(0, self.molecule.natom()): extern.addCharge(self.molecule.Z(A), self.molecule.x(A), self.molecule.y(A), self.molecule.z(A))
[docs]class QMMM(object): def __init__(self): self.charges = [] self.diffuses = [] self.extern = PsiMod.ExternalPotential()
[docs] def addDiffuse(self, diffuse): """Function to add a diffuse charge field *diffuse*.""" self.diffuses.append(diffuse)
[docs] def addChargeBohr(self, Q, x, y, z): """Function to add a point charge of magnitude *Q* at position (*x*, *y*, *z*) Bohr. """ self.charges.append([Q, x, y, z])
[docs] def addChargeAngstrom(self, Q, x, y, z): """Function to add a point charge of magnitude *Q* at position (*x*, *y*, *z*) Angstroms. """ self.charges.append([Q, x / physconst.psi_bohr2angstroms, y / physconst.psi_bohr2angstroms, z / physconst.psi_bohr2angstroms])
def __str__(self): s = ' ==> QMMM <==\n\n' s = s + ' => Charges (a.u.) <=\n\n' s = s + ' %11s %11s %11s %11s\n' % ('Z', 'x', 'y', 'z') for k in range(0, len(self.charges)): s = s + ' %11.7f %11.3E %11.3E %11.3E\n' % (self.charges[k][0], self.charges[k][1], self.charges[k][2], self.charges[k][3]) s = s + '\n' s = s + ' => Diffuses <=\n\n' for k in range(0, len(self.diffuses)): s = s + str(self.diffuses[k]) return s
[docs] def populateExtern(self): """Function to define a charge field external to the molecule through point and diffuse charges. """ # Charges for charge in self.charges: self.extern.addCharge(charge[0], charge[1], charge[2], charge[3]) # Diffuses for diffuse in self.diffuses: diffuse.populateExtern(self.extern)