Python System.setProperty方法代码示例

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

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

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

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

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

printEnergies(system.energies())
 
print("\nEnergy with respect to cutoff length\n")
print("  Distance   Group    Shifted    ReactionField   Atomistic")

for i in range(10,200,5):
    x = i*0.1

    switchfunc = HarmonicSwitchingFunction(x*angstrom, (x-0.5)*angstrom)
    system.setProperty("switchingFunction", switchfunc)

    print("%12.8f  %12.8f  %12.8f  %12.8f  %12.8f" % (x, system.energy(group_coul).value(),
              system.energy(shift_coul).value(), system.energy(field_coul).value(),

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

# Tell the system about the periodic box
system.setProperty("space", space)
system.add( SpaceWrapper(Vector(0), waters) )
system.applyConstraints()

# Save a binary representation of the system
# to the file "waterbox.s3"
import Sire.Stream
Sire.Stream.save( system, "waterbox.s3" )

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
def test_props(verbose=False):
    sys = System()

    box0 = PeriodicBox( Vector(10.0,10.0,10.0) )
    box1 = PeriodicBox( Vector(20.0,20.0,20.0) )

    if verbose:
        print(box0)
        print(box0.volume())
        print(box1.volume())

    assert(not sys.containsProperty("space"))

    sys.add( InterCLJFF("cljff") )

    if verbose:
        print(sys)
        print(sys.property("space"))
        print(sys.userProperties().propertyKeys())
        print(sys.builtinProperties().propertyKeys())

    assert(sys.containsProperty("space"))
    assert_equal( sys.property("space"), Cartesian() )

    sys.setProperty( "space0", LinkToProperty("space", FFIdx(0)) )

    if verbose:
        print(sys.property("space0"))

    assert(sys.containsProperty("space0"))

    sys.setProperty("space0", box0)

    if verbose:
        print(sys.property("space"))

    assert_equal(sys.property("space0"), box0)

    sys.setProperty("space1", box1)

    sys.setProperty("combined_space", CombineSpaces("space0", "space1"))

    assert_equal(sys.property("space1"), box1)

    if verbose:
        print(sys.properties().propertyKeys())

        print(sys.property("combined_space"))
        print(sys.property("combined_space").volume())

    assert_almost_equal( sys.property("combined_space").volume().value(), 
                         sys.property("space0").volume().value() + sys.property("space1").volume().value(), 5 )

    space3 = PeriodicBox( Vector(5,5,5) )
    sys.setProperty("space0", space3)

    assert_equal( sys.property("space0"), space3 )

    if verbose:
        print(sys.property("combined_space"))
        print(sys.property("combined_space").volume())

    assert_almost_equal( sys.property("combined_space").volume().value(), 
                         sys.property("space0").volume().value() + sys.property("space1").volume().value(), 5 )

    sys.removeProperty("space0")

    if verbose:
        print(sys.properties().propertyKeys())

    assert( not sys.containsProperty("space0") )

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
field_coul = field_clj.components().coulomb()

system = System()

for forcefield in forcefields:
    forcefield.add(protein)
    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("switchingFunction", HarmonicSwitchingFunction(10*angstrom, 9.5*angstrom))

printEnergies(system.energies())
 
print("\nEnergy with respect to cutoff length\n")
print("  Distance   Group    Shifted    ReactionField  ")

for i in range(10,501,5):
    x = i*0.1

    switchfunc = HarmonicSwitchingFunction(x*angstrom, (x-0.5)*angstrom)
    system.setProperty("switchingFunction", switchfunc)

    print("%12.8f  %12.8f  %12.8f  %12.8f" % (x, system.energy(group_coul).value(),
              system.energy(shift_coul).value(), system.energy(field_coul).value()))

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
# create the periodic box space for the atoms
box = PeriodicBox( Vector(box_size[0].value(),
                          box_size[1].value(),
                          box_size[2].value()) )

# create a forcefield to calculate the intermolecular coulomb and LJ (CLJ)
# energy between all krypton atoms
interff = InterCLJFF("CLJ")
interff.setProperty("space", box)
interff.add(mols)

# create a simulation system to hold the forcefield and atoms
system = System()
system.add(interff)
system.add(mols)
system.setProperty("space", box)

# add a wrapper that wraps the atoms back into the box
system.add( SpaceWrapper( Vector(0), mols ) ) 

# create rigid body translation moves for the atoms
rb_moves = RigidBodyMC(mols)
rb_moves.setMaximumTranslation(max_translate)
rb_moves.setTemperature(temperature)

# create volume moves to change the box size
vol_moves = VolumeMove(mols)
vol_moves.setMaximumVolumeChange( mols.nMolecules() * 0.1 * angstrom3 )
vol_moves.setTemperature(temperature)
vol_moves.setPressure(pressure)

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
box0 = PeriodicBox( Vector(10.0,10.0,10.0) )
box1 = PeriodicBox( Vector(20.0,20.0,20.0) )

print(box0)
print(box0.volume())
print(box1.volume())

from Sire.MM import *
sys.add( InterCLJFF("cljff") )

print(sys)
print(sys.property("space"))
print(sys.userProperties().propertyKeys())
print(sys.builtinProperties().propertyKeys())

sys.setProperty( "space0", LinkToProperty("space", FFIdx(0)) )

print(sys.property("space0"))

sys.setProperty("space0", box0)

print(sys.property("space"))

sys.setProperty("space1", box1)

sys.setProperty("combined_space", CombineSpaces("space0", "space1"))

print(sys.properties().propertyKeys())

print(sys.property("combined_space"))
print(sys.property("combined_space").volume())

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
salt.add(m1)
salt.add(m2)
salt.add(m3)
#salt.add(c)
#salt.add(c2)

cljff = InterCLJFF("salt-salt")
cljff.add(salt)
internalff.add(salt)

system = System()
system.add(salt)
system.add(cljff)
system.add(internalff)

system.setProperty("space", PeriodicBox( Vector(20,20,20) ) )

system.add( SpaceWrapper(Vector(0,0,0),salt) )

t.start()                                       
print("Initial energy = %s" % system.energy())
print("(took %d ms)" % t.elapsed())

mdmove = MolecularDynamics( salt, VelocityVerlet(), 
                            {"velocity generator":MaxwellBoltzmann(25*celsius)} )

mdmove = MolecularDynamics(salt, DLMRigidBody() )

mdmove.setTimeStep(1*femtosecond)
do_mc = False

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
    solvent.update(mol)

system = System()
system.add(solute)
system.add(solvent)

solvent_cljff = InterCLJFF("solvent_cljff")
solvent_cljff.add(solvent)
solute_solvent_cljff = InterGroupCLJFF("solute_solvent_cljff")
solute_solvent_cljff.add( solute, MGIdx(0) )
solute_solvent_cljff.add( solvent, MGIdx(1) )

system.add(solvent_cljff)
system.add(solute_solvent_cljff)

system.setProperty("space", PeriodicBox( Vector(-18.3854,-18.66855,-18.4445),
                                         Vector( 18.3854, 18.66855, 18.4445) ) )

print("Equilibrating the system...")

solvent_move = RigidBodyMC( PrefSampler(solute.moleculeAt(0), solvent, 200*angstrom2) )
solvent_move.setMaximumTranslation( 0.2 * angstrom )
solvent_move.setMaximumRotation( 5 * degrees )

solute_move = RigidBodyMC( solute )
solute_move.setMaximumTranslation( 0.2 * angstrom )
solute_move.setMaximumRotation( 5 * degrees )

moves = WeightedMoves()
moves.add( solvent_move, 100 )
moves.add( solute_move, 1 )

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
  
# The protein-solvent energy 
protein_solventff = InterGroupCLJFF("protein:solvent")
protein_solventff.add(protein, MGIdx(0))
protein_solventff.add(solvent, MGIdx(1))

# Here is the list of all forcefields
forcefields = [ solute_intraff, solute_intraclj,
                solventff, solute_solventff,
                protein_intraff, protein_intraclj,
                solute_proteinff, protein_solventff ] 
# Add these forcefields to the system
for forcefield in forcefields:
    system.add(forcefield)

system.setProperty( "space", space )
system.setProperty( "switchingFunction", 
                    HarmonicSwitchingFunction(coulomb_cutoff, coulomb_feather,
                                              lj_cutoff, lj_feather) )
system.setProperty( "combiningRules", VariantProperty(combining_rules) )

total_nrg = solute_intraclj.components().total() + solute_intraff.components().total() +\
    solventff.components().total() + solute_solventff.components().total() +\
    protein_intraclj.components().total() + protein_intraff.components().total() + \
    solute_proteinff.components().total() + protein_solventff.components().total() 

e_total = system.totalComponent()
system.setComponent( e_total, total_nrg )

# Add a space wrapper that wraps all molecules into the box centered at (0,0,0)
#system.add( SpaceWrapper(Vector(0,0,0), all) )

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
def makeSim(system, ligand_mol, watersys):
    """Create simulation systems with and without the ligand and return those systems together
       with the moves"""

    stage1 = System("with_ligand")
    stage2 = System("without_ligand")

    if system.containsProperty("reflection center"):
        reflection_center = system.property("reflection center").toVector()[0]
        reflection_radius = float(str(system.property("reflection sphere radius")))

        stage1.setProperty("reflection center", AtomCoords(CoordGroup(1,reflection_center)))
        stage1.setProperty("reflection sphere radius", VariantProperty(reflection_radius))

        stage2.setProperty("reflection center", AtomCoords(CoordGroup(1,reflection_center)))
        stage2.setProperty("reflection sphere radius", VariantProperty(reflection_radius))

    # create a molecule group for fixed atoms (everything except the mobile water)
    fixed_group = MoleculeGroup("fixed")

    if MGName("fixed_molecules") in system.mgNames():
        fixed_group.add( system[ MGName("fixed_molecules") ] )

    if MGName("mobile_solutes") in system.mgNames():
        fixed_group.add( system[MGName("mobile_solutes")] )

    if MGName("protein_sidechains") in system.mgNames() or \
       MGName("protein_backbones") in system.mgNames():

        all_proteins = Molecules()

        try:
            protein_sidechains = system[MGName("protein_sidechains")]
            all_proteins.add(protein_sidechains.molecules())
        except:
            pass

        try:
            protein_backbones = system[MGName("protein_backbones")]
            all_proteins.add(protein_backbones.molecules())
        except:
            pass

        try:
            boundary_molecules = system[MGName("boundary_molecules")]
            all_proteins.add(boundary_molecules.molecules())
        except:
            pass

        for molnum in all_proteins.molNums():
            protein_mol = all_proteins[molnum].join()
            fixed_group.add(protein_mol)

    stage1_fixed_group = MoleculeGroup(fixed_group)
    stage2_fixed_group = MoleculeGroup(fixed_group)

    stage1_fixed_group.add(ligand_mol)
    stage2_fixed_group.remove(ligand_mol)

    mobile_group = MoleculeGroup("mobile_group")
    if MGName("mobile_solvents") in system.mgNames():
        mobile_group.add( system[MGName("mobile_solvents")] )

    stage1_mobile_group = MoleculeGroup(mobile_group)
    stage2_mobile_group = MoleculeGroup(mobile_group)

    # now find water molecules from the water system that can be substituted for the ligand
    watermols = findOverlappingWaters(ligand_mol, watersys)

    stage2_mobile_group.add(watermols)

    print("The number of stage 1 fixed non-solvent molecules is %d." % stage1_fixed_group.nMolecules())
    print("The number of stage 1 mobile solvent molecules is %d." % stage1_mobile_group.nMolecules())

    print("The number of stage 2 fixed non-solvent molecules is %d." % stage2_fixed_group.nMolecules())
    print("The number of stage 2 mobile solvent molecules is %d." % stage2_mobile_group.nMolecules())

    # write a PDB of all of the fixed molecules
    PDB().write(stage1_mobile_group, "stage1_mobile_atoms.pdb")
    PDB().write(stage2_mobile_group, "stage2_mobile_atoms.pdb")
    PDB().write(stage1_fixed_group, "stage1_fixed_atoms.pdb")
    PDB().write(stage2_fixed_group, "stage2_fixed_atoms.pdb")

    # create the forcefields

    if use_fast_ff.val:
        stage1_ff = InterFF("ff")
        stage2_ff = InterFF("ff")
        stage1_ff.setCLJFunction( CLJShiftFunction(Cartesian(), coul_cutoff.val, lj_cutoff.val) )
        stage2_ff.setCLJFunction( CLJShiftFunction(Cartesian(), coul_cutoff.val, lj_cutoff.val) )
        
        if disable_grid.val:
            stage1_ff.disableGrid()
            stage2_ff.disableGrid()
        else:
            stage1_ff.enableGrid()
            stage1_ff.setGridSpacing(grid_spacing.val)
            stage1_ff.setGridBuffer(grid_buffer.val)
            stage2_ff.enableGrid()
            stage2_ff.setGridSpacing(grid_spacing.val)
#.........这里部分代码省略.........

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setProperty [as 别名]
def createSystem():
    protomsdir = "%s/Work/ProtoMS" % os.getenv("HOME")

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

    protoms.addParameterFile( "%s/parameter/amber99.ff" % protomsdir )
    protoms.addParameterFile( "%s/parameter/solvents.ff" % protomsdir )
    protoms.addParameterFile( "%s/parameter/gaff.ff" % protomsdir )
    protoms.addParameterFile( solute_params )

    solute = PDB().readMolecule(solute_file)
    solute = solute.edit().rename(solute_name).commit()
    solute = protoms.parameterise(solute, ProtoMS.SOLUTE)

    perturbation = solute.property("perturbations")

    lam = Symbol("lambda")
    lam_fwd = Symbol("lambda_{fwd}")
    lam_bwd = Symbol("lambda_{bwd}")

    initial = Perturbation.symbols().initial()
    final = Perturbation.symbols().final()

    solute = solute.edit().setProperty("perturbations",
                perturbation.recreate( (1-lam)*initial + lam*final ) ).commit()

    solute_fwd = solute.edit().renumber().setProperty("perturbations",
                perturbation.substitute( lam, lam_fwd ) ).commit()
    solute_bwd = solute.edit().renumber().setProperty("perturbations",
                perturbation.substitute( lam, lam_bwd ) ).commit()

    solvent = PDB().read(solvent_file)

    tip4p = solvent.moleculeAt(0).molecule()
    tip4p = tip4p.edit().rename(solvent_name).commit()
    tip4p = protoms.parameterise(tip4p, ProtoMS.SOLVENT)

    tip4p_chgs = tip4p.property("charge")
    tip4p_ljs = tip4p.property("LJ")

    for i in range(0,solvent.nMolecules()):
        tip4p = solvent.moleculeAt(i).molecule()
        tip4p = tip4p.edit().rename(solvent_name) \
                            .setProperty("charge", tip4p_chgs) \
                            .setProperty("LJ", tip4p_ljs) \
                            .commit()

        solvent.update(tip4p)

    system = System()

    solutes = MoleculeGroup("solutes")
    solutes.add(solute)
    solutes.add(solute_fwd)
    solutes.add(solute_bwd)

    solvent = MoleculeGroup("solvent", solvent)

    all = MoleculeGroup("all")
    all.add(solutes)
    all.add(solvent)

    system.add(solutes)
    system.add(solvent)
    system.add(all)

    solventff = InterCLJFF("solvent:solvent")
    solventff.add(solvent)

    solute_intraff = InternalFF("solute_intraff")
    solute_intraff.add(solute)

    solute_fwd_intraff = InternalFF("solute_fwd_intraff")
    solute_fwd_intraff.add(solute_fwd)

    solute_bwd_intraff = InternalFF("solute_bwd_intraff")
    solute_bwd_intraff.add(solute_bwd)

    solute_intraclj = IntraCLJFF("solute_intraclj")
    solute_intraclj.add(solute)

    solute_fwd_intraclj = IntraCLJFF("solute_fwd_intraclj")
    solute_fwd_intraclj.add(solute_fwd)

    solute_bwd_intraclj = IntraCLJFF("solute_bwd_intraclj")
    solute_bwd_intraclj.add(solute_bwd)

    solute_solventff = InterGroupCLJFF("solute:solvent")
    solute_solventff.add(solute, MGIdx(0))
    solute_solventff.add(solvent, MGIdx(1))

    solute_fwd_solventff = InterGroupCLJFF("solute_fwd:solvent")
    solute_fwd_solventff.add(solute_fwd, MGIdx(0))
    solute_fwd_solventff.add(solvent, MGIdx(1))

    solute_bwd_solventff = InterGroupCLJFF("solute_bwd:solvent")
    solute_bwd_solventff.add(solute_bwd, MGIdx(0))
    solute_bwd_solventff.add(solvent, MGIdx(1))

    forcefields = [ solventff, solute_intraff, solute_intraclj, solute_solventff,
#.........这里部分代码省略.........