Python System.setConstant方法代码示例

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import setConstant [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 setConstant [as 别名]
mol = mol.edit().rename("SB2").commit()

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

perturbations = mol.property("perturbations")

print(perturbations)

print(perturbations.requiredSymbols())
print(perturbations.requiredProperties())

lam = perturbations.symbols().Lambda()

system = System()

solute = MoleculeGroup("solute", mol)

system.add(solute)

system.setConstant(lam, 0.0)
system.add( PerturbationConstraint(solute) )

print(system.constraintsSatisfied())

for i in range(0,101,10):
    system.setConstant(lam, 0.01 * i)

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

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

# We also add a SpaceWrapper that ensures any moves on mobile_mols

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