from Sire.IO import *

from Sire.Squire import *

from Sire.MM import *

from Sire.Mol import *

from Sire.Move import *

from Sire.System import *

from Sire.Maths import *

from Sire.CAS import *

from Sire.Units import *

import os

import sys

import math

# 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() )

# also define the difference in energy

e_diff = Symbol("E_diff")

system.setComponent( e_diff, e_mm - e_qmmm )

# Now add a monitor that will record the average of the

# difference in energy between QM/MM and MM

system.add( "difference", MonitorComponent(e_diff, Average()) )

# now print out the QM/MM energy, MM energy and difference

print("Total QM/MM energy = %s" % system.energy(e_qmmm))

print("Total MM energy = %s" % system.energy(e_mm))

print("Difference = %s" % system.energy(e_diff))

# create a random number generator for the simulation

rand = RanGenerator()

# now create a function to run one multiple-timestep

# Monte Carlo move (using 5000 fast MM moves)

def multiMCMove(system, nfast=5000):

    # calculate the old QM/MM and MM energies

    old_qmmm = system.energy(e_qmmm).value()

    old_mm = system.energy(e_mm).value()

    # Save the old conformation of the system

    old_system = System(system)

    # now perform 'nfast' MC moves using the

    # MM potential

    move = RigidBodyMC(watergroup)

    move.setEnergyComponent(e_mm)

    # actually perform the moves

    move.move(system, nfast, False)

    # calculate the new QM/MM and MM energies

    new_qmmm = system.energy(e_qmmm).value()

    new_mm = system.energy(e_mm).value()

    # now use the multiple timestep MC test

    # to accept or reject the move (at 298 K)

    beta = 1.0 / (k_boltz * 298.0)

    test = math.exp( -beta * ( (new_qmmm - new_mm) - \

                               (old_qmmm - old_mm) ) )

    print("E_QMMM(new) - E_MM(new) = %.2f" % (new_qmmm - new_mm))

    print("E_QMMM(old) - E_MM(old) = %.2f" % (old_qmmm - old_mm))

    print("(E_QMMM(new)-E_MM(new)) - (E_QMMM(old)-E_MM(old)) = %.2f" % \

             ((new_qmmm-new_mm) - (old_qmmm - old_mm)) )

    print("exp{-beta * [(E_QMMM(new)-E_MM(new)) - (E_QMMM(old)-E_MM(old)]} = %.2f" % \

                 test )

    randnum = rand.rand()

    print("Random number between 0 and 1 = %.2f" % randnum)

    if randnum < test:

        print("Random number is smaller: Move accepted!")

        # record the average difference

        system.collectStats()

        return system

    else:

        print("Random number is larger: Move rejected!")

        # record the average difference

        system.collectStats()

        return old_system

# perform 20 multiple timestep MC moves

for i in range(1,21):

    print("\nPerforming move %d" % i)

    system = multiMCMove(system)

    PDB().write(system.molecules(), "coords%003d.pdb" % i)

# print the average energy difference

avg_diff = system[ MonitorName("difference") ].accumulator().average()

print("\nTotal QM/MM energy = %s" % system.energy(e_qmmm))

print("Total MM energy = %s" % system.energy(e_mm))

print("Difference = %s" % system.energy(e_diff))

print("\nCoordinates saved to 'coords*.pdb'")

print("\nAverage energy difference = %s kcal mol-1" % avg_diff)

← Back