Python System.applyConstraints方法代码示例

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import applyConstraints [as 别名]

#.........这里部分代码省略.........
    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
    for solute in solutes:
        solute_group.add(solute)

    for protein in proteins:
        protein_group.add(protein)

    for water in waters:
        solvent_group.add(water)
        water_group.add(water)

    for solvent in solvents:
        solvent_group.add(solvent)
    
    for ion in ions:
        solvent_group.add(ion)
        ion_group.add(ion)

    if solute_group.nMolecules() > 0:
        system.add(solute_group)

    if protein_group.nMolecules() > 0:
        system.add(protein_group)

    if solvent_group.nMolecules() > 0:
        system.add(solvent_group)

    if water_group.nMolecules() > 0:
        system.add(water_group)

    if ion_group.nMolecules() > 0:
        system.add(ion_group)    

    print("Number of solute molecules == %s" % solute_group.nMolecules()) 
    print("Number of protein molecules == %s" % protein_group.nMolecules())
    print("Number of ions == %s" % ion_group.nMolecules())
    print("Number of water molecules == %s" % water_group.nMolecules())
    print("Number of solvent molecules == %s" % solvent_group.nMolecules())
    print("(solvent group is waters + ions + unidentified single-residue molecules)")

    system.setProperty("space", space)
    system.add( SpaceWrapper( Vector(0), all_group ) )
    system.applyConstraints()

    print("Returning the constructed system")

    return system

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import applyConstraints [as 别名]
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())

pol_tip4p = system[MGIdx(0)][pol_tip4p.number()].molecule()

print("MM charges\n",tip4p.property("charge"), \
                     tip4p.evaluate().charge({"charge":"charge"}))
print("Fixed charges\n",pol_tip4p.property("fixed_charge"), \
                        pol_tip4p.evaluate().charge({"charge":"fixed_charge"}))
print("Induced charges\n",pol_tip4p.property("induced_charge"), \
                          pol_tip4p.evaluate().charge({"charge":"induced_charge"}))
print("New charges\n",pol_tip4p.property("charge"), \
                      pol_tip4p.evaluate().charge({"charge":"charge"}))

grid = RegularGrid(tip4p.evaluate().center(), 10, 1.0*angstrom)

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import applyConstraints [as 别名]

#.........这里部分代码省略.........
        mobile_fixed.add(fixed_group, MGIdx(1))
        other_nrg += mobile_fixed.components().total()
        forcefields.append(mobile_fixed)
    else:
        mobile_fixed = GridFF("mobile-fixed")
        mobile_fixed = setCLJProperties(mobile_fixed, space)
        mobile_fixed = setGridProperties(mobile_fixed)

        # we use mobile_buffered_group as this group misses out atoms that are bonded
        # to fixed atoms (thus preventing large energies caused by incorrect non-bonded calculations)
        mobile_fixed.add(mobile_buffered_mols, MGIdx(0))
        mobile_fixed.addFixedAtoms(fixed_group)
        other_nrg += mobile_fixed.components().total()
        forcefields.append(mobile_fixed)

    # intramolecular energy of the protein
    if protein_intra_mols.nMolecules() > 0:
        protein_intraclj = IntraCLJFF("protein_intraclj")
        protein_intraclj = setCLJProperties(protein_intraclj, space)

        protein_intraff = InternalFF("protein_intra")

        for molnum in protein_intra_mols.molNums():
            protein_mol = protein_intra_mols[molnum].join()
            protein_intraclj.add(protein_mol)
            protein_intraff.add(protein_mol)

        other_nrg += protein_intraclj.components().total()
        other_nrg += protein_intraff.components().total()
        forcefields.append(protein_intraclj)
        forcefields.append(protein_intraff)

    # intramolecular energy of any other solutes
    if solute_intra_mols.nMolecules() > 0:
        solute_intraclj = IntraCLJFF("solute_intraclj")
        solute_intraclj = setCLJProperties(solute_intraclj, space)

        solute_intraff = InternalFF("solute_intra")

        for molnum in solute_intra_mols.molNums():
            solute_mol = solute_intra_mols[molnum].join()
            solute_intraclj.add(solute_mol)
            solute_intraff.add(solute_mol)

        other_nrg += solute_intraclj.components().total()
        other_nrg += solute_intraff.components().total()
        forcefields.append(solute_intraclj)
        forcefields.append(solute_intraff)

    ###
    ### NOW ADD THE FORCEFIELDS TO THE SYSTEM
    ###
    ###
    ### SETTING THE FORCEFIELD EXPRESSIONS
    ###

    lam = Symbol("lambda")

    e_slow = ((1-lam) * ligand_qm_nrg) + (lam * ligand_mm_nrg) + other_nrg
    e_fast = ligand_mm_nrg + other_nrg

    de_by_dlam = ligand_mm_nrg - ligand_qm_nrg

    for forcefield in forcefields:
        system.add(forcefield)

    system.setConstant(lam, 0.0)

    system.setComponent(Symbol("E_{fast}"), e_fast)
    system.setComponent(Symbol("E_{slow}"), e_slow)
    system.setComponent(Symbol("dE/dlam"), de_by_dlam)
    system.setComponent( system.totalComponent(), e_slow )
 
    system.setProperty("space", space)
    
    if space.isPeriodic():
        # ensure that all molecules are wrapped into the space with the ligand at the center
        print("Adding in a space wrapper constraint %s, %s" % (space, ligand_mol.evaluate().center()))
        system.add( SpaceWrapper( ligand_mol.evaluate().center(), all_group ) )
        system.applyConstraints()

    print("\nHere are the values of all of the initial energy components...")
    t.start()
    printEnergies(system.energies())
    print("(these took %d ms to evaluate)\n" % t.elapsed())

    # Create a monitor to monitor the free energy average
    system.add( "dG/dlam", MonitorComponent(Symbol("dE/dlam"), AverageAndStddev()) )

    if intermolecular_only.val:
        print("\n\n## This simulation uses QM to model *only* the intermolecular energy between")
        print("## the QM and MM atoms. The intramolecular energy of the QM atoms is still")
        print("## modelled using MM.\n")
    else:
        print("\n\n## This simulation uses QM to model both the intermolecular and intramolecular")
        print("## energies of the QM atoms. Because the this, we have to adjust the 'zero' point")
        print("## of the QM potential. You need to add the value %s kcal mol-1 back onto the" % zero_energy)
        print("## QM->MM free energy calculated using this program.\n")

    return system