from Sire.IO import *

from Sire.Squire import *

from Sire.System import *

from Sire.MM import *

from Sire.Mol import *

from Sire.CAS 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) )

# New stuff :-)

# we need a forcefield to hold the intermolecular

# energy between the water molecules

waterff = InterCLJFF("waterff")

waterff.setSwitchingFunction(cutoff)

waterff.setSpace(space)

waterff.add(waters)

# Now add all of the forcefields together into a single

# simulation system

system = System()

watergroup = MoleculeGroup("waters")

watergroup.add(waters)

system.add(watergroup)

system.add(mmff)

system.add(qmmmff)

system.add(waterff)

# Define the total MM energy

e_mm = Symbol("E_MM")

e_qmmm = Symbol("E_QMMM")

system.setComponent( e_mm,  mmff.components().total() + \

                            waterff.components().total() )

system.setComponent( e_qmmm, qmmmff.components().total() + \

                             mmff.components().lj() + \

                             waterff.components().total() )

# Now define lambda - set it to 0

lam = Symbol("lambda")

system.setConstant(lam, 0)

# Now define the lambda scaled energy, such that

# it equals e_qmmm at lambda = 0, and e_mm at lambda = 1

e_lam = Symbol("E_lambda")

system.setComponent( e_lam, (1-lam)*e_qmmm + lam * e_mm )

# calculate the energy at the value of lambda supplied by the user

system.setConstant(lam, float(sys.argv[3]))

e_l = system.energy(e_lam)

print("\n%s  %s" % (float(sys.argv[3]), e_l.value()))

← Back