Python Sire.System类代码示例

def test_central(verbose=False):
    testsys = System(system)

    moves = RigidBodyMC(waters)
    moves.setReflectionSphere(Vector(0), radius)

    if verbose:
        print("Performing 5000 moves 10 times...")

        PDB().write(testsys.molecules(), "test_central0000.pdb")

    for i in range(1,11):
        moves.move(testsys, 5000, False)

        if verbose:
            PDB().write(testsys.molecules(), "test_central%0004d.pdb" % i)

    # the reflection sphere should ensure that none of the 
    # water molecules diffuse outside the sphere
    mols = testsys.molecules()

    for molnum in mols.molNums():
        mol = mols[molnum]

        dist = mol.evaluate().center().length()

        if verbose:
            print("%s : %s A" % (molnum.value(), dist))

        assert( dist <= radius.value() )

def createSystem(molecules):

    moleculeNumbers = molecules.molNums()
    moleculeList = []

    for moleculeNumber in moleculeNumbers:
        molecule = molecules.molecule(moleculeNumber).molecule()
        moleculeList.append(molecule)

    all = MoleculeGroup("all")

    for molecule in moleculeList[0:]:
        all.add(molecule)

    gridwater = MoleculeGroup("gridwater")
    otherwater = MoleculeGroup("otherwater")
    solutes = MoleculeGroup("solutes")
  
    # Add these groups to the System
    system = System()

    system.add(all)
    system.add(gridwater)
    system.add(otherwater)
    system.add(solutes)

    return system

def buildSystem(forcefields):
    system = System()
    for forcefield in forcefields:
        system.add(forcefield)

    system.setProperty("space", space)
    system.setProperty("switchingFunction", HarmonicSwitchingFunction(8*angstrom, 7.5*angstrom))

    system.setProperty("shiftDelta", VariantProperty(2.0))
    system.setProperty("coulombPower", VariantProperty(0))

    return system

    mol = mol.edit().rename("T4P") \
                    .setProperty("charge", charges) \
                    .setProperty("LJ", ljs) \
             .commit()

    cljff.add(mol)

    if i > 3:
        free_mols.add(mol)
    else:
        sync_mols.add(mol)

ms = t.elapsed()
print("Parameterised all of the water molecules (in %d ms)!" % ms)

system = System()

system.add(cljff)

print("Initial energy = %s" % system.energy())

mc = RigidBodyMC(free_mols)
sync_mc = RigidBodyMC(sync_mols)
sync_mc.setSynchronisedTranslation(True)
sync_mc.setSynchronisedRotation(True)

nodes = Cluster.getNode()
this_thread = nodes.borrowThisThread()

moves = WeightedMoves()
moves.add(mc, 2)

def test_local(verbose = False):
    testsys = System(system)

    moves = RigidBodyMC(waters)

    local_radius = 1 * angstrom

    if verbose:
        print("Setting local reflection spheres...")

    for molnum in waters.molNums():
        center = waters[molnum].evaluate().center()
        moves.setReflectionSphere(molnum, center, local_radius)

        if verbose:
            print("%s : %s == %s, %s == %s" % (molnum.value(), \
                center, moves.reflectionSphereCenter(molnum), \
                local_radius.value(), moves.reflectionSphereRadius(molnum).value()))

        assert( moves.reflectionSphereCenter(molnum) == center )
        assert( moves.reflectionSphereRadius(molnum) == local_radius )

    if verbose:                               
        print("Performing 5000 moves 10 times...")
        PDB().write(testsys.molecules(), "test_local0000.pdb")

    for i in range(1,11):
        moves.move(testsys, 5000, False)

        if verbose:
            PDB().write(testsys.molecules(), "test_local%0004d.pdb" % i)

    # the reflection sphere should ensure that none of the
    # water molecules diffuse outside the sphere
    mols = testsys.molecules()

    for molnum in mols.molNums():
        mol = mols[molnum]
        oldmol = waters[molnum]

        dist = Vector.distance( mol.evaluate().center(), oldmol.evaluate().center() )

        if verbose:
            print("%s : %s A" % (molnum.value(), dist))

        assert( dist <= local_radius.value() )

def createSystem(molecules):
    # print("Applying flexibility and zmatrix templates...")
    print("Creating the system...")

    moleculeNumbers = molecules.molNums()
    moleculeList = []

    for moleculeNumber in moleculeNumbers:
        molecule = molecules.molecule(moleculeNumber).molecule()
        moleculeList.append(molecule)

    molecules = MoleculeGroup("molecules")
    ions = MoleculeGroup("ions")

    for molecule in moleculeList:
        natoms = molecule.nAtoms()
        if natoms == 1:
            ions.add(molecule)
        else:
            molecules.add(molecule)

    all = MoleculeGroup("all")
    all.add(molecules)
    all.add(ions)

    # Add these groups to the System
    system = System()

    system.add(all)
    system.add(molecules)
    system.add(ions)

    return system

def addRDF(rdf, pdbfile, xscfile=None):
    """Add to the RDF 'rdf' the distances calculated using the coordinates 
       from the PDB file 'pdbfile', using the xscfile 'xscfile' to get the 
       dimensions of the periodic box. If 'xscfile' is None, then 
       no periodic boundary conditions are used."""

    #first get the space in which to calculate intermolecular distances
    space = Cartesian()

    if xscfile:
        lines = open(xscfile,"r").readlines()

        words = lines[0].split()
        mincoords = Vector( float(words[0]), float(words[1]), float(words[2]) )
        maxcoords = Vector( float(words[3]), float(words[4]), float(words[5]) )

        space = PeriodicBox(mincoords, maxcoords)

    #now load all of the molecules
    mols = PDB().read(pdbfile)
                                                
    #create a system to hold the molecules, and add them
    system = System()
    system.add( MoleculeGroup("molecules", mols) )

    #give the space to the system
    system.add( InterCLJFF() )  # bug! need to add InterCLJFF so 
                                # that the system has a space property. This
                                # is fixed in new version of Sire, but is broken
                                # in your version

    system.setProperty("space", space)

    #add the RDF - this calculates the RDF for this PDB file and adds it to 'rdf'
    rdf.monitor(system)

def test_fixed_center(verbose = False):
    ligand = Sire.Stream.load("../io/osel.s3")
    ligand = ligand.edit().setProperty("center", wrap(ligand.evaluate().center())).commit()

    old_center = ligand.property("center")

    intraff = InternalFF("intraff")
    intraff.add(ligand)

    intraclj = IntraCLJFF("intraclj")
    intraclj.add(ligand)

    system = System()
    system.add(intraff)
    system.add(intraclj)

    mols = MoleculeGroup("mols")
    mols.add(ligand)
    system.add(mols)

    intramove = InternalMove(mols)
    rbmove = RigidBodyMC(mols)
    rbmove.setMaximumTranslation(0*angstrom)

    moves = WeightedMoves()
    moves.add(intramove, 1)
    moves.add(rbmove, 1)

    for i in range(0,10):
        system = moves.move(system, 25, False)

        if verbose:
            print("Completed 25 moves...")

        ligand = system[ligand.number()].molecule()

        new_center = ligand.property("center")

        if verbose:
            print("Old center = %s" % old_center)
            print("New center = %s" % new_center)

        assert_almost_equal( old_center.x(), new_center.x(), 1 )
        assert_almost_equal( old_center.y(), new_center.y(), 1 )
        assert_almost_equal( old_center.z(), new_center.z(), 1 )

protoms_dir = os.getenv("HOME") + "/Work/ProtoMS"

protoms = ProtoMS( "%s/protoms2" % protoms_dir )

protoms.addParameterFile("%s/parameter/amber99.ff" % protoms_dir )
protoms.addParameterFile("%s/parameter/gaff.ff" % protoms_dir )
protoms.addParameterFile("test/ff/sb2.ff")

sb2 = protoms.parameterise(sb2, ProtoMS.SOLUTE)

qmff = QMFF("MopacFF")
mopac = Mopac()
qmff.setQuantumProgram( mopac )

qmff.add(sb2)

system = System()

system.add(qmff)

zmat_move = ZMatMove(qmff[MGIdx(0)])

moves = SameMoves(zmat_move)

import Sire.Stream
Sire.Stream.save( (system,moves), "test/Squire/mopacsim.s3" )

print("Energy before == %f kcal mol-1" % system.energy().to(kcal_per_mol))
system = moves.move(system)
print("Energy after == %f kcal mol-1" % system.energy().to(kcal_per_mol))

def createSystemFrom(molecules, space, system_name, naming_scheme = NamingScheme()):
    """Create a new System from the passed molecules and space,
       sorting the molecules into different molecule groups based on the
       passed naming scheme"""

    system = System(system_name)

    # If requested, change the water model for all water molecules
    if water_model.val == "tip4p":
        molnums = molecules.molNums()
        new_molecules = Molecules()

        print("Forcing all water molecules to use the %s water model..." % water_model.val)
        print("Converting %d molecules..." % len(molnums))
        i = 0
        for molnum in molnums:
            molecule = molecules[molnum].molecule()

            if i % 100 == 0:
                print("%d" % i)                
                sys.stdout.flush()

            elif i % 10 == 0:
                print(".", end=' ')
                sys.stdout.flush()

            i += 1

            if molecule.nAtoms() == 3:
                # this could be a TIP3P water
                resname =str(molecule.residue().name().value()).lower()

                if resname == "wat" or resname == "t3p":
                    new_molecule = convertTip3PtoTip4P(molecule)
                    if new_molecule:
                        molecule = new_molecule

            new_molecules.add(molecule)

        print("%d" % i)

        molecules = new_molecules

    nmols = molecules.nMolecules()

    print("Number of molecules == %s" % nmols)
    print("System space == %s" % space)

    if nmols == 0:
        return system

    print("Assigning molecules to molecule groups...")
    solute_group = MoleculeGroup(naming_scheme.solutesGroupName().value())
    protein_group = MoleculeGroup(naming_scheme.proteinsGroupName().value())
    solvent_group = MoleculeGroup(naming_scheme.solventsGroupName().value())
    water_group = MoleculeGroup(naming_scheme.watersGroupName().value())
    ion_group = MoleculeGroup(naming_scheme.ionsGroupName().value())
    all_group = MoleculeGroup(naming_scheme.allMoleculesGroupName().value())

    # The all molecules group has all of the molecules
    all_group.add(molecules)

    system.add(all_group)

    # Run through each molecule and decide what type it is...
    molnums = molecules.molNums()
    molnums.sort()

    central_molecule = None

    solutes = []
    proteins = []
    solvents = []
    waters = []
    ions = []

    for molnum in molnums:
        molecule = molecules[molnum].molecule()

        resnams = getResidueNames(molecule)

        if naming_scheme.isSolute(resnams):
            solutes.append(molecule)

        elif naming_scheme.isProtein(resnams):
            proteins.append(molecule)

        elif naming_scheme.isWater(resnams):
            waters.append(molecule)

        elif naming_scheme.isIon(resnams):
            ions.append(molecule)

        elif molecule.nResidues() == 1:
            solvents.append(molecule)

        else:
            solutes.append(molecule)

    # Ok - we have now divided everything up into groups
#.........这里部分代码省略.........

def loadQMMMSystem():
    """This function is called to set up the system. It sets everything
       up, then returns a System object that holds the configured system"""

    print("Loading the system...")

    t = QTime()

    if os.path.exists(s3file.val):
        print("Loading existing s3 file %s..." % s3file.val)
        loadsys = Sire.Stream.load(s3file.val)

    else:
        print("Loading from Amber files %s / %s..." % (topfile.val, crdfile.val))
        # Add the name of the ligand to the list of solute molecules
        sys_scheme = NamingScheme()
        sys_scheme.addSoluteResidueName(ligand_name.val)

        # Load up the system. This will automatically find the protein, solute, water, solvent
        # and ion molecules and assign them to different groups
        loadsys = createSystem(topfile.val, crdfile.val, sys_scheme)
        ligand_mol = findMolecule(loadsys, ligand_name.val)

        if ligand_mol is None:
            print("Cannot find the ligand (%s) in the set of loaded molecules!" % ligand_name.val)
            sys.exit(-1)

        # Center the system with the ligand at (0,0,0)
        loadsys = centerSystem(loadsys, ligand_mol)
        ligand_mol = loadsys[ligand_mol.number()].molecule()

        if reflection_radius.val is None:
            loadsys = addFlexibility(loadsys, naming_scheme=sys_scheme )
        else:
            loadsys = addFlexibility(loadsys, Vector(0), reflection_radius.val, naming_scheme=sys_scheme)

        Sire.Stream.save(loadsys, s3file.val)

    ligand_mol = findMolecule(loadsys, ligand_name.val)

    if ligand_mol is None:
        print("Cannot find the ligand (%s) in the set of loaded molecules!" % ligand_name.val)
        sys.exit(-1)

    # Now build the QM/MM system
    system = System("QMMM system")

    if loadsys.containsProperty("reflection center"):
        reflect_center = loadsys.property("reflection center").toVector()[0]
        reflect_radius = float(str(loadsys.property("reflection sphere radius")))

        system.setProperty("reflection center", AtomCoords(CoordGroup(1,reflect_center)))
        system.setProperty("reflection sphere radius", VariantProperty(reflect_radius))
        space = Cartesian()
    else:
        space = loadsys.property("space")

    if loadsys.containsProperty("average solute translation delta"):
        system.setProperty("average solute translation delta", \
                           loadsys.property("average solute translation delta"))

    if loadsys.containsProperty("average solute rotation delta"):
        system.setProperty("average solute rotation delta", \
                           loadsys.property("average solute rotation delta"))

    # create a molecule group to hold all molecules
    all_group = MoleculeGroup("all")

    # create a molecule group for the ligand
    ligand_group = MoleculeGroup("ligand")
    ligand_group.add(ligand_mol)
    all_group.add(ligand_mol)

    groups = []
    groups.append(ligand_group)

    # pull out the groups that we want from the two systems

    # create a group to hold all of the fixed molecules in the bound leg
    fixed_group = MoleculeGroup("fixed_molecules")
    if MGName("fixed_molecules") in loadsys.mgNames():
        fixed_group.add( loadsys[ MGName("fixed_molecules") ] )

    if save_pdb.val:
        # write a PDB of the fixed atoms in the bound and free legs
        if not os.path.exists(outdir.val):
            os.makedirs(outdir.val)

        PDB().write(fixed_group, "%s/fixed.pdb" % outdir.val)

    # create a group to hold all of the mobile solute molecules
    mobile_solutes_group = MoleculeGroup("mobile_solutes")
    if MGName("mobile_solutes") in loadsys.mgNames():
        mobile_solutes_group.add( loadsys[MGName("mobile_solutes")] )
        mobile_solutes_group.remove(ligand_mol)
        if mobile_solutes_group.nMolecules() > 0:
            all_group.add(mobile_solutes_group)
    
    groups.append(mobile_solutes_group)

#.........这里部分代码省略.........

field_cljff.setUseReactionField(True)
field_cljff.setReactionFieldDielectric(78.0)

atom_cljff = InterCLJFF("atomistic_cutoff")
atom_cljff.setUseAtomisticCutoff(True)

group_coul = group_cljff.components().coulomb()
shift_coul = shift_cljff.components().coulomb()
field_coul = field_cljff.components().coulomb()
atom_coul = atom_cljff.components().coulomb()

forcefields = [ group_cljff, shift_cljff, field_cljff, atom_cljff ]

switchfunc = HarmonicSwitchingFunction(10*angstrom, 9.5*angstrom)

system = System()

for forcefield in forcefields:
    forcefield.add(waters)
    system.add(forcefield)

def printEnergies(nrgs):
    keys = list(nrgs.keys())
    keys.sort()

    for key in keys:
        print("%25s : %12.8f" % (key, nrgs[key]))

system.setProperty("space", space)
system.setProperty("switchingFunction", switchfunc)

from Sire.Vol import *
from Sire.Maths import *
from Sire.Units import *

## Simple script to load up a water box from waterbox.top/.crd
## and to create a System that can be used for Monte Carlo
## simulation of these waters

# Load the waters and the periodic box from the crd/top files
(waters, space) = Amber().readCrdTop("waterbox.crd", "waterbox.top")

# Rename the 'waters' group to 'water'
waters.setName("water")

# Create a system to hold the system to be simulated
system = System()

# Add the waters to the system
system.add(waters)

# create a forcefield to calculate the intermolecular
# Coulomb and Lennard-Jones energies of the waters (CLJ)
cljff = InterFF("cljff")

# Use a shift-electrostatics with a 10 angstrom cutoff
cljff.setCLJFunction( CLJShiftFunction(10*angstrom) )

# Add the waters and set the cutoff to 10 angstroms
cljff.add(waters)

# Add the forcefield to the system

words = open("test/io/water.xsc", "r").readline().split()
space = PeriodicBox( Vector( float(words[0]), float(words[1]), float(words[2]) ),
                     Vector( float(words[3]), float(words[4]), float(words[5]) ) )

qm_water = waters.moleculeAt(0).molecule()
waters.remove(qm_water)

molpro = Molpro()

qmff = QMMMFF("qmmm")
qmff.setQuantumProgram(molpro)

qmff.add(qm_water, MGIdx(0))
qmff.add(waters, MGIdx(1))

system = System()
system.add(qmff)
system.setProperty("switchingFunction", switchfunc)
system.setProperty("space", space)

print("Calculating the system energies...")
print(system.energies())

polchgs = PolariseCharges(waters, qmff.components().total(),
                          CoulombProbe(1*mod_electron))

system.add(polchgs)
system.add(polchgs.selfEnergyFF())

print("Applying the polarisation constraint...")
system.applyConstraints()

                        .setProperty("polarisability", tip4p_pol) \
                        .setProperty("connectivity", connectivity) \
                 .commit()

    waters.update(water)

print("Constructing the forcefields...")

pol_tip4p = waters.moleculeAt(0).molecule()
waters.remove(pol_tip4p)

cljff = InterGroupCLJFF("pol_tip4p-water")
cljff.add(pol_tip4p, MGIdx(0))
cljff.add(waters, MGIdx(1))

system = System()
system.add(cljff)

print(system.energies())

polchgs = PolariseCharges(cljff[MGIdx(0)], cljff.components().coulomb(),
                          CoulombProbe(1*mod_electron))

system.add(polchgs)
system.add(polchgs.selfEnergyFF())

print("Applying the polarisation constraint...")
system.applyConstraints()

print(system.energies())