from Sire.IO import *

from Sire.Squire import *

from Sire.MM import *

from Sire.Mol import *

from Sire.Units import *

import os

import sys

# Get the name of the solute from the first

# command line argument (should be 'PHE' or 'TRP')

molname = sys.argv[1]

# Load the solute in a box of water, together with the

# information about the periodic box

(mols, space) = Amber().readCrdTop("%s.crd" % molname, "%s.top" % molname)

# pick out the solute from the box of molecules

mol = mols[ MolWithResID(molname) ].molecule()

# now get all of the waters (every molecule except the solute)

waters = mols.molecules()

waters.remove(mol)

# Set the cutoff distance as 10 angstroms

cutoff = HarmonicSwitchingFunction( 10 * angstrom )

# Create an MM forcefield to calculate the

# MM electrostatic energy between the solute and the waters

mmff = InterGroupCLJFF("mmff")

mmff.setSwitchingFunction(cutoff)

mmff.setSpace(space)

# Create a QM/MM forcefield to calculate the QM/MM

# electrostatic energy between the solute and the waters

qmmmff = QMMMFF("qmmmff")

qmmmff.setSwitchingFunction(cutoff)

qmmmff.setSpace(space)

qmmmff.setIntermolecularOnly(True)

# We will use the SQM program to calculate QM energies.

# Create the program and set the semi-empirical method

# from the passed string

qmprog = SQM()

qmprog.setMethod(sys.argv[2])

# Tell qmprog where to find the SQM executable

amberhome = os.getenv("AMBERHOME")

if amberhome is None:

    print("Please set the AMBERHOME environment variable to point")

    print("to the Amber14 installation.")

    sys.exit(-1)

qmprog.setExecutable( "%s/bin/sqm" % amberhome )

# Tell the QM/MM forcefield to use this program

qmmmff.setQuantumProgram(qmprog)

# Now add the solute to group 0 of both forcefields

mmff.add( mol, MGIdx(0) )

qmmmff.add( mol, MGIdx(0) )

# Now add the waters to group 1 of both forcefields

mmff.add( waters, MGIdx(1) )

qmmmff.add( waters, MGIdx(1) )

# Now calculate the QM and MM energies

mm_energy = mmff.energy(mmff.components().coulomb())

qmmm_energy = qmmmff.energy()

# Now print the QM/MM energy

print("\nQM/MM energy = %s" % qmmm_energy )

# Now print the MM coulomb energy

print("\nMM electrostatic energy = %s" % mm_energy )

# Calculate the difference between the QM/MM and MM energies

print("\nDifference = %s" % (mm_energy - qmmm_energy))

← Back