Python System.totalComponent方法代码示例

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
#system.add( DistanceComponent( r34, solute.atom(AtomName("A3")), solute.atom(AtomName("A4")) ) )
#system.add( "bond34", MonitorComponent(r34, RecordValues()), 5)

#theta123 = Symbol("theta123")
#system.add( AngleComponent( theta123, solute.atom(AtomName("A1")), solute.atom(AtomName("A2")), solute.atom(AtomName("A3")) ) )
#system.add( "theta123", MonitorComponent(theta123, RecordValues()), 5)

#theta234 = Symbol("theta234")
#system.add( AngleComponent( theta234, solute.atom(AtomName("A2")), solute.atom(AtomName("A3")), solute.atom(AtomName("A4")) ) )
#system.add( "theta234", MonitorComponent(theta234, RecordValues()), 5)

phi1234 = Symbol("phi1234")
system.add( DihedralComponent( phi1234, solute.atom(AtomName("A1")), solute.atom(AtomName("A2")), solute.atom(AtomName("A3")), solute.atom(AtomName("A4")) ) )
system.add( "phi1234", MonitorComponent(phi1234, RecordValues()), 5)

e_total = system.totalComponent()
system.add( "total_energy", MonitorComponent(e_total, Average()) )

#system.add( "trajectory", TrajectoryMonitor(system[MGIdx(0)]), 1 )

# Run the simulation
for i in range(0,nblocks):
    print("Running nmove...")
    system = moves.move(system, nmoves, True)

print("Analysis...")
# make an histogram...
nbins=31
#histoAnalysis(nbins=nbins,
#              avgnrg = system.monitor( MonitorName("total_energy") ).accumulator().average(),
#              values=system.monitor( MonitorName("bond12") ).accumulator().values(),

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

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
oscillator = protoms.parameterise(oscillator, ProtoMS.SOLUTE)

print("...parameterisation complete!")

internalff = InternalFF("InternalFF")

internalff.add( oscillator )

system = System()
system.add( internalff )

lam = Symbol("lambda")

system.setComponent(lam, 0.01)

system.setComponent(system.totalComponent(), lam * internalff.components().total())

system.add( "average energy", MonitorComponent(system.totalComponent()) )

lambda_values = [ 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 ]

replicas = Replicas(system, len(lambda_values))

replicas.setLambdaComponent(lam)

for i in range(0,len(lambda_values)):
    replicas.setLambdaValue(i, lambda_values[i])

zmatmove = ZMatMove( internalff.groups()[0] )
zmatmove.setTemperature( 298 * kelvin )

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
def _pvt_calculateEnergy(lamval, verbose):
    cljff01 = InterGroupCLJFF("cljff01")
    cljff02 = InterGroupCLJFF("cljff02")

    cljff01.add( water0, MGIdx(0) )
    cljff01.add( water1, MGIdx(1) )

    cljff02.add( water0, MGIdx(0) )
    cljff02.add( water2, MGIdx(1) )

    soft_cljff01 = InterGroupSoftCLJFF("soft_cljff01")
    soft_cljff02 = InterGroupSoftCLJFF("soft_cljff02")

    soft_cljff01.add( water0, MGIdx(0) )
    soft_cljff01.add( water1, MGIdx(1) )

    soft_cljff02.add( water0, MGIdx(0) )
    soft_cljff02.add( water2, MGIdx(1) )

    ref = MoleculeGroup("ref", water0)
    group_a = MoleculeGroup("group_a", water1)
    group_b = MoleculeGroup("group_b", water2)

    dlam = 0.001
    lamval_f = lamval + dlam

    lam = Symbol("lambda")
    lam_f = Symbol("lambda_f")

    soft_nrg = (1-lam) * soft_cljff01.components().total(0) + lam * soft_cljff02.components().total(0)
    soft_nrg_f = (1-lam_f) * soft_cljff01.components().total(1) + lam_f * soft_cljff02.components().total(1)

    de_soft = soft_nrg_f - soft_nrg

    nrg = ((1-lam) * cljff01.components().total()) + (lam * cljff02.components().total())
    nrg_f = ((1-lam_f) * cljff01.components().total()) + (lam_f * cljff02.components().total())

    de = nrg_f - nrg

    soft_cljff01.setProperty("alpha0", VariantProperty(lamval))
    soft_cljff02.setProperty("alpha0", VariantProperty(1-lamval))
    soft_cljff01.setProperty("alpha1", VariantProperty(lamval_f))
    soft_cljff02.setProperty("alpha1", VariantProperty(1-lamval_f))

    soft_cljff01.setProperty("coulombPower", VariantProperty(0))
    soft_cljff01.setProperty("shiftDelta", VariantProperty(1.1))
    soft_cljff02.setProperty("coulombPower", VariantProperty(0))
    soft_cljff02.setProperty("shiftDelta", VariantProperty(1.1))

    sys = System()
    sys.add(cljff01)
    sys.add(cljff02)
    sys.add(soft_cljff01)
    sys.add(soft_cljff02)
    sys.add(ref)
    sys.add(group_a)
    sys.add(group_b)
    sys.setComponent(lam, lamval)
    sys.setComponent(lam_f, lamval_f)
    sys.setComponent(sys.totalComponent(), nrg)
    sys.setComponent(Symbol("E_{total_f}"), nrg_f)
    sys.setComponent(Symbol("dE"), de)
    sys.setComponent(Symbol("E_soft_{total}"), soft_nrg)
    sys.setComponent(Symbol("E_soft_{total_f}"), soft_nrg_f)
    sys.setComponent(Symbol("dE_soft"), de_soft)

    nrgmon = FreeEnergyMonitor(ref, group_a, group_b)

    soft_nrgmon = FreeEnergyMonitor(ref, group_a, group_b)
    soft_nrgmon.setCoulombPower(0)
    soft_nrgmon.setShiftDelta(1.1)

    sys.add( "nrgmon", nrgmon )
    sys.add( "soft_nrgmon", soft_nrgmon)

    sys.collectStats()

    nrgmon = sys[ MonitorName("nrgmon") ]
    dg = nrgmon.freeEnergies()[0].average()

    soft_nrgmon = sys[ MonitorName("soft_nrgmon") ]
    soft_dg = soft_nrgmon.freeEnergies()[0].average()

    sys_dg = sys.energy(Symbol("dE")).value()
    sys_soft_dg = sys.energy(Symbol("dE_soft")).value()

    if verbose:
        print("%s : %s versus %s (should be equal)" % (lamval,dg,sys_dg))
        print("%s : %s versus %s (should be equal)" % (lamval,soft_dg,sys_soft_dg))

    assert_almost_equal(dg, sys_dg, 5)
    assert_almost_equal(soft_dg, sys_soft_dg, 5)

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
# Sire performs simulations on simulation System objects.
# A System groups together all of the molecules, forcefields etc.
# into a single unit
system = System()
system.add(cljff)

# Here we create the equation used to calculate the total energy of the
# system. We define a symbol, lambda, which we use here to scale up and
# down the coulomb component of cljff
lam = Symbol("lambda")
total_nrg = cljff.components().lj() + lam * cljff.components().coulomb()

# Here we tell the system to use our equation to calculate the total
# energy, and we set the value of lambda to 0.5
system.setComponent( system.totalComponent(), total_nrg )
system.setConstant( lam, 0.5 )

# Now we create a MoleculeGroup that groups together all of the 
# molecules to be moved
mobile_mols = MoleculeGroup("mobile_molecules")
mobile_mols.add(first_water)
mobile_mols.add(second_water)

# We add the molecule group to the system
system.add(mobile_mols)

# We create a periodic boundaries space, and give it to the system
space = PeriodicBox( Vector(5,5,5) )
system.setProperty( "space", space )

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
oldsys = System()
oldsys.add(ligmol)
oldsys.add(residues)
oldsys.add(oldff)
oldsys.add(old0ff)
oldsys.add(old1ff)

etotal = Symbol("E_{total}")
ecoul = Symbol("E_{coul}")
elj = Symbol("E_{lj}")

oldsys.setComponent( etotal, oldff.components().total() - \
                            old0ff.components().total() - \
                            old1ff.components().total() )

if oldsys.totalComponent() != etotal:
    oldsys.setComponent( oldsys.totalComponent(), etotal )

oldsys.setComponent( ecoul, oldff.components().coulomb() - \
                            old0ff.components().coulomb() - \
                            old1ff.components().coulomb() )
oldsys.setComponent( elj, oldff.components().lj() - \
                          old0ff.components().lj() - \
                          old1ff.components().lj() )

newsys = System()
newsys.add(ligmol)
newsys.add(residues)
newsys.add(newff)
newsys.setComponent( ecoul, newff.components().coulomb() )
newsys.setComponent( elj, newff.components().lj() )

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
print("Initial energy = %s" % system.energy())

data = save(system)

system = load(data)

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

mc = RigidBodyMC(cljff.group(MGIdx(0)))

moves = SameMoves(mc)

#give a lambda coordinate that turns off the coulomb energy
lam = Symbol("lambda")

system.setComponent( system.totalComponent(), cljff.components().lj() + \
                                              lam * cljff.components().coulomb() )

system.setComponent( lam, 0.0 )
print(system.energies())

system.setComponent( lam, 1.0 )
print(system.energies())

#create 5 replicas that map from lambda=0 to lambda=1
replicas = Replicas(system, 5)
replicas.setSubMoves(moves)
replicas.setNSubMoves(1000)

replicas.setLambdaComponent(lam)

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

#.........这里部分代码省略.........
    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)
            stage2_ff.setGridBuffer(grid_buffer.val)

        stage1_ff.add(stage1_mobile_group)
        stage1_ff.setFixedAtoms(stage1_fixed_group.molecules())
        stage2_ff.add(stage2_mobile_group)
        stage2_ff.setFixedAtoms(stage2_fixed_group.molecules())

        stage1.add(stage1_ff)
        stage1.setComponent(stage1.totalComponent(), stage1_ff.components().total())
        stage2.add(stage1_ff)
        stage2.setComponent(stage2.totalComponent(), stage2_ff.components().total())

    else:
        # forcefield holding the energy between the mobile atoms and  
        # the fixed atoms
        if disable_grid.val:
            stage1_mobile_fixed = InterGroupCLJFF("mobile-fixed")
            stage1_mobile_fixed = setCLJProperties(stage1_mobile_fixed)
            stage1_mobile_fixed = setFakeGridProperties(stage1_mobile_fixed)
            stage1_mobile_fixed.add(stage1_mobile_group, MGIdx(0))
            stage1_mobile_fixed.add(stage1_fixed_group, MGIdx(1))

            stage2_mobile_fixed = InterGroupCLJFF("mobile-fixed")
            stage2_mobile_fixed = setCLJProperties(stage2_mobile_fixed)
            stage2_mobile_fixed = setFakeGridProperties(stage2_mobile_fixed)
            stage2_mobile_fixed.add(stage2_mobile_group, MGIdx(0))
            stage2_mobile_fixed.add(stage2_fixed_group, MGIdx(1))
        else:
            stage1_mobile_fixed = GridFF("mobile-fixed")
            stage1_mobile_fixed = setCLJProperties(stage1_mobile_fixed)
            stage1_mobile_fixed = setGridProperties(stage1_mobile_fixed)

            stage1_mobile_fixed.add(stage1_mobile_group, MGIdx(0))
            stage1_mobile_fixed.addFixedAtoms(stage1_fixed_group)

            stage2_mobile_fixed = GridFF("mobile-fixed")
            stage2_mobile_fixed = setCLJProperties(stage2_mobile_fixed)
            stage2_mobile_fixed = setGridProperties(stage2_mobile_fixed)

            stage2_mobile_fixed.add(stage2_mobile_group, MGIdx(0))
            stage2_mobile_fixed.addFixedAtoms(stage2_fixed_group)

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
                    .setProperty("LJ", ljs) \
             .commit()

    cljff.add(mol)

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

system = System()

system.add(cljff)

lam = Symbol("lambda")

system.setComponent( lam, 0.2 )
system.setComponent( system.totalComponent(), lam * cljff.components().total() )

mc = RigidBodyMC(cljff.group(MGIdx(0)))

moves = SameMoves(mc)

def testStream(c):
    t.start()

    data = Sire.Stream.save(c)

    ms = t.elapsed()

    print(("Streaming %s took %d ms" % (c.what(), ms)))
    print(("%s takes up %d bytes" % (c.what(),data.size())))

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
cljff.add(swapwaters, MGIdx(0))
cljff.add(waters, MGIdx(1))
cljff.setSpace( Cartesian() )
cljff.setShiftElectrostatics(True)

swap_swapff = InterCLJFF("swap-swap")
swap_swapff.setSpace(Cartesian())
swap_swapff.setSwitchingFunction( HarmonicSwitchingFunction(25*angstrom, 25*angstrom, 10*angstrom, 10*angstrom) )
swap_swapff.add(swapwaters)
grid_system.add(swap_swapff)
exp_system.add(swap_swapff)

grid_system.add(gridff)
exp_system.add(cljff)

grid_system.setComponent( grid_system.totalComponent(),  \
                          gridff.components().total() + swap_swapff.components().total() )
exp_system.setComponent( exp_system.totalComponent(), \
                              cljff.components().total() + swap_swapff.components().total() )

print((grid_system.energies()))
print((exp_system.energies()))

print(("\nGrid energy equals: %s. Explicit energy equals: %s." % \
          (grid_system.energy(), exp_system.energy())))

diff = grid_system.energy() - exp_system.energy()
print(("The difference is %s\n" % diff))

rbmc = RigidBodyMC(swapwaters)
rbmc.setReflectionSphere(center_point, 7.5*angstrom)

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import totalComponent [as 别名]
gridff2.setSpace( Cartesian() )
gridff2.setShiftElectrostatics(True)
#gridff2.setUseReactionField(True)
#gridff2.setUseAtomisticCutoff(True)

swap_swapff = InterCLJFF("swap-swap")
swap_swapff.setSpace(Cartesian())
swap_swapff.setSwitchingFunction( HarmonicSwitchingFunction(25*angstrom, 25*angstrom, 10*angstrom, 10*angstrom) )
swap_swapff.add(swapwaters)
grid_system.add(swap_swapff)
grid_system2.add(swap_swapff)

grid_system.add(gridff)
grid_system2.add(gridff2)

grid_system.setComponent( grid_system.totalComponent(),  \
                          gridff.components().total() + swap_swapff.components().total() )
grid_system2.setComponent( grid_system2.totalComponent(), \
                           gridff2.components().total() + swap_swapff.components().total() )

print(grid_system.energies())
print(grid_system2.energies())

print("\nOld Grid energy equals: %s. New GridFF energy equals: %s." % \
          (grid_system.energy(), grid_system2.energy()))

diff = grid_system.energy() - grid_system2.energy()
print("The difference is %s\n" % diff)

rbmc = RigidBodyMC(swapwaters)
rbmc.setReflectionSphere(center_point, 7.5*angstrom)

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

#.........这里部分代码省略.........
    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,
                               solute_fwd_intraff, solute_fwd_intraclj, solute_fwd_solventff,
                               solute_bwd_intraff, solute_bwd_intraclj, solute_bwd_solventff ]

    for forcefield in forcefields:
        system.add(forcefield)

    xsc_line = open(solvent_file, "r").readlines()[0]
    words = xsc_line.split()

    #HEADER box  -12.5  -12.5  -12.5   12.5   12.5   12.5
    space = PeriodicBox( Vector( float(words[2]), float(words[3]), float(words[4]) ),
                         Vector( float(words[5]), float(words[6]), float(words[7]) ) )

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

    e_total = system.totalComponent()
    e_fwd = Symbol("E_{fwd}")
    e_bwd = Symbol("E_{bwd}")

    total_nrg = solventff.components().total() + \
                solute_intraclj.components().total() + solute_intraff.components().total() + \
                solute_solventff.components().total()

    fwd_nrg = solventff.components().total() + \
              solute_fwd_intraclj.components().total() + solute_fwd_intraff.components().total() + \
              solute_fwd_solventff.components().total()

    bwd_nrg = solventff.components().total() + \
              solute_bwd_intraclj.components().total() + solute_bwd_intraff.components().total() + \
              solute_bwd_solventff.components().total()

    system.setComponent( e_total, total_nrg )
    system.setComponent( e_fwd, fwd_nrg )
    system.setComponent( e_bwd, bwd_nrg )

    system.setConstant(lam, 0.0)
    system.setConstant(lam_fwd, 0.0)
    system.setConstant(lam_bwd, 0.0)

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

    system.add( PerturbationConstraint(solutes) )

    system.add( ComponentConstraint( lam_fwd, Min( lam + delta_lambda, 1 ) ) )
    system.add( ComponentConstraint( lam_bwd, Max( lam - delta_lambda, 0 ) ) )

    de_fwd = Symbol("de_fwd")
    de_bwd = Symbol("de_bwd")

    system.setComponent( de_fwd, fwd_nrg - total_nrg )
    system.setComponent( de_bwd, total_nrg - bwd_nrg )

    system.add( "total_energy", MonitorComponent(e_total, Average()) )
    system.add( "de_fwd", MonitorComponent(de_fwd, FreeEnergyAverage(temperature)) )
    system.add( "de_bwd", MonitorComponent(de_bwd, FreeEnergyAverage(temperature)) )

    system.setComponent(lam, 0.0)
    print "LAMBDA=0   : Energy = %f kcal mol-1" % system.energy().to(kcal_per_mol)
    print "             (%f, %f)" % (system.energy(e_fwd).to(kcal_per_mol),
                                   system.energy(e_bwd).to(kcal_per_mol))

    system.setComponent(lam, 0.5)
    print "LAMBDA=0.5 : Energy = %f kcal mol-1" % system.energy().to(kcal_per_mol)
    print "             (%f, %f)" % (system.energy(e_fwd).to(kcal_per_mol),
                                   system.energy(e_bwd).to(kcal_per_mol))

    system.setComponent(lam, 1.0)
    print "LAMBDA=1.0 : Energy = %f kcal mol-1" % system.energy().to(kcal_per_mol)
    print "             (%f, %f)" % (system.energy(e_fwd).to(kcal_per_mol),
                                   system.energy(e_bwd).to(kcal_per_mol))
    return system