Python System.energy方法代码示例

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energy [as 别名]
    mol = mol.edit().rename("T4P") \
                    .setProperty("charge", charges) \
                    .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)

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() + \

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energy [as 别名]
for i in range(0,len(lambda_values)):
    replicas.setLambdaValue(i, lambda_values[i])

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

nsubmoves = 1000

replicas.setSubMoves( SameMoves(zmatmove) )
replicas.setNSubMoves(nsubmoves)

# Average energy should be 1/2 kT
theo_nrg = 0.5 * gasr * 298

print("Running a simulation - initial energy = %f kcal mol-1" % system.energy().to(kcal_per_mol))

repexmove = RepExMove()

lambda_trajectory = []

i = -1

def printInfo(replicas):

    lamtraj = replicas.lambdaTrajectory()
    lambda_trajectory.append(lamtraj)

    ids = replicas.replicaIDs()

    for j in range(0,replicas.count()):

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

moves = SameMoves(rbmc)

PDB().write(grid_system2.molecules(), "test0000.pdb")

t = QTime()

for i in range(1,11):
    print("Moving the system...")

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

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

# Add the CLJ forcefield to the system
system.add(cljff)

# Add the molecule group containing the waters to the system - we will give
# this molecule group the name 'waters' so that we can extract it from
# the system at a later point in the script
waters.setName("waters")
system.add(waters)

# Add a constraint that waters are mapped back into the periodic box
system.add( SpaceWrapper(Vector(0,0,0), waters) )

print "Calculating the starting energy... (should be -16364.5 kcal mol-1)"
print "...Initial energy = %s" % system.energy()

# Now create the rigid body move (RigidBodyMC, from Sire.Move) that act on the waters
rbmc = RigidBodyMC(waters)

# Set the maximum amount by which to translate and rotate the water molecules
# to 0.25 A and 15 degrees (angstrom and degrees are from Sire.Units)
rbmc.setMaximumTranslation( 0.25*angstrom )
rbmc.setMaximumRotation( 15*degrees )

# Set the temperature of the moves to 25 Celsius (Celsius is from Sire.Units)
rbmc.setTemperature( 25*celsius )

# Now create the volume move (VolumeMove, from Sire.Move) that maintains a constant pressure
volmc = VolumeMove(waters)

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

# group these two moves together
moves = WeightedMoves()
moves.add( rb_moves, mols.nMolecules() )
moves.add( vol_moves, 1 )

# print the initial energy and coordinates
print("0:   %s" % system.energy())
PDB().write(system.molecules(), "output000000.pdb") 

# now run the simulation in blocks of 1000 moves
nmoves = 0
while nmoves < num_moves:
    system = moves.move(system, 1000, False)
    nmoves += 1000

    # print out the energy and coordinates every 1000 moves
    print("%s  %s" % (nmoves, system.energy()))
    print(moves)

    PDB().write(system.molecules(), "output%000006d.pdb" % nmoves)

print("Complete!")

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energy [as 别名]
    mol = mols.moleculeAt(i).molecule()

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

    cljff.add(mol)
    mols.update(mol)

system = System()

system.add(cljff)

print("System energy equals...")
print(system.energy())

group0 = MoleculeGroup("group0")
group1 = MoleculeGroup("group1")

group0.add( mols.moleculeAt(100) )
group0.add( mols.moleculeAt(101) )

group1.add( mols.moleculeAt(102) )
group1.add( mols.moleculeAt(103) )
group1.add( mols.moleculeAt(104) )

cljff2 = InterGroupCLJFF("group_energy")
cljff2.add( group0, MGIdx(0) )
cljff2.add( group1, MGIdx(1) )
cljff2.setSpace(vol)

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

    oldsys = System()
    newsys = System()
    parsys = System()

    oldsys.add(mols)
    newsys.add(mols)
    parsys.add(mols)

    oldsys.add(newff)
    newsys.add(newff)
    parsys.add(parff)

    t = QElapsedTimer()

    oldsys.mustNowRecalculateFromScratch()
    newsys.mustNowRecalculateFromScratch()
    parsys.mustNowRecalculateFromScratch()

    t.start()
    nrgs = oldsys.energies()
    oldns = t.nsecsElapsed()

    t.start()
    nrgs = newsys.energies()
    newns = t.nsecsElapsed()

    t.start()
    nrgs = parsys.energies()
    parns = t.nsecsElapsed()

    oldcnrg = oldsys.energy( newff.components().coulomb() ).value()
    oldljnrg = oldsys.energy( newff.components().lj() ).value()

    newcnrg = newsys.energy( newff.components().coulomb() ).value()
    newljnrg = newsys.energy( newff.components().lj() ).value()

    parcnrg = parsys.energy( parff.components().coulomb() ).value()
    parljnrg = parsys.energy( parff.components().lj() ).value()

    if verbose:
        print("\nStarting energy")
        print("OLD SYS:  %s  %s  %s  : %s ms" % (oldcnrg+oldljnrg,oldcnrg,oldljnrg,
                                                 0.000001*oldns))
        print("NEW SYS:  %s  %s  %s  : %s ms" % (newcnrg+newljnrg,newcnrg,newljnrg,
                                                 0.000001*newns))
        print("PAR SYS:  %s  %s  %s  : %s ms" % (parcnrg+parljnrg,parcnrg,parljnrg,
                                                 0.000001*parns))

    assert_almost_equal( oldcnrg, newcnrg )
    assert_almost_equal( oldljnrg, newljnrg )

    assert_almost_equal( parcnrg, newcnrg )
    assert_almost_equal( parljnrg, newljnrg )

    moves = RigidBodyMC(mols)
    moves.disableOptimisedMoves()
    moves.setGenerator( RanGenerator( 42 ) )

    optmoves = RigidBodyMC(mols)
    optmoves.enableOptimisedMoves()
    optmoves.setGenerator( RanGenerator( 42 ) )

    parmoves = RigidBodyMC(mols)
    parmoves.enableOptimisedMoves()
    parmoves.setGenerator( RanGenerator( 42 ) )

    t.start()
    moves.move(oldsys, nmoves, False)
    move_oldns = t.nsecsElapsed()

    old_naccepted = moves.nAccepted()
    old_nrejected = moves.nRejected()

    t.start()
    optmoves.move(newsys, nmoves, False)
    move_newns = t.nsecsElapsed()

    new_naccepted = optmoves.nAccepted()
    new_nrejected = optmoves.nRejected()

    t.start()
    parmoves.move(parsys, nmoves, False)
    move_parns = t.nsecsElapsed()

    par_naccepted = parmoves.nAccepted()
    par_nrejected = parmoves.nRejected()

    t.start()
    nrgs = oldsys.energies()
    oldns = t.nsecsElapsed()
    
    t.start()
    nrgs = newsys.energies()
    newns = t.nsecsElapsed()
    
    t.start()
    nrgs = parsys.energies()
    parns = t.nsecsElapsed()
#.........这里部分代码省略.........

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

print("\nTotal energy ")
print(system.energy())

print("Components energies ")
for component in list(system.energyComponents().keys()):
    print(component, system.energyComponents().value(component) * kcal_per_mol)

# Note that tip3p water are likely to have bonds between hydrogen atoms.
PDB().write(all, "out.pdb")

print("The AMBER11/sander energies for this system are ") 
print("""
# NSTEP =        0   TIME(PS) =       0.000  TEMP(K) =     0.00  PRESS =     0.0
# Etot   =    -47010.2216  EKtot   =         0.0000  EPtot      =    -47010.2216
# BOND   =       898.1982  ANGLE   =      5310.2620  DIHED      =      2922.5644
# 1-4 NB =       790.8755  1-4 EEL =      7702.0145  VDWAALS    =      7345.0484
# EELEC  =    -71979.1846  EHBOND  =         0.0000  RESTRAINT  =         0.0000

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energy [as 别名]
print(bonds.potentials())
print(angles.potentials())
print(dihedrals.potentials())

#intraff = InternalFF("intraff")
intraclj = IntraCLJFF("intraclj")

#intraff.add(ethane)
intraclj.add(ethane)

solute = MoleculeGroup("solute", ethane)
solute.add(ethane)

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

md = MolecularDynamics(solute, VelocityVerlet()) 
  #                     {"velocity generator" : MaxwellBoltzmann(25*celsius)})

md.setTimeStep(1*femtosecond)

PDB().write(system.molecules(), "test0000.pdb")

for i in range(1,250):
    md.move(system, 1)
    print(system.energy(), md.kineticEnergy(), (system.energy()+md.kineticEnergy()))
    PDB().write(system.molecules(), "test%0004d.pdb" % i)

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energy [as 别名]
    print((header.systemInfo()))

    c2 = Sire.Stream.load(data)

    ms = t.elapsed()
  
    print(("Reading the data took %d ms" % ms))
    print(c)
    print(c2)

testStream(system)

data = Sire.Stream.save(system)

print("Probing the system...")
print((system.energy()))
print((system.energies()))

system = Sire.Stream.load(data)

print((system.energy()))
print((system.energies()))

print("\nGetting data info...")

t.start()
Sire.Stream.save( system, "test/SireStream/tmp_testdata.sire" )
ms = t.elapsed()

print(("Saving a system to a file took %d ms" % ms))

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energy [as 别名]
def test_moves(verbose = False):
    newsys = System()
    oldsys = System()
    optsys = System()
    parsys = System()

    res = MoleculeGroup("residues")
    protein = mols[ MolWithResID("ALA") ].molecule()

    for residue in protein.residues():
        res.add(residue)

    newsys.add(newff)
    oldsys.add(oldff)
    optsys.add(newff)
    parsys.add(parff)

    newsys.add(res)
    oldsys.add(res)
    optsys.add(res)
    parsys.add(res)

    moves = RigidBodyMC(res)
    moves.disableOptimisedMoves()

    moves.setMaximumTranslation( 0.2 * angstrom )
    moves.setMaximumRotation( 0.5 * degrees )

    opt_moves = RigidBodyMC(res)
    opt_moves.enableOptimisedMoves()

    opt_moves.setMaximumTranslation( 0.2 * angstrom )
    opt_moves.setMaximumRotation( 0.5 * degrees )

    t = QElapsedTimer()

    if verbose:
        print("Calculating initial energy...")

    t.start()
    oldnrg = oldsys.energy().value()
    oldns = t.nsecsElapsed()

    t.start()
    newnrg = newsys.energy().value()
    newns = t.nsecsElapsed()

    t.start()
    optnrg = optsys.energy().value()
    optns = t.nsecsElapsed()

    t.start()
    parnrg = parsys.energy().value()
    parns = t.nsecsElapsed()

    if verbose:
        print("\nTIMES: old = %s ms, new = %s ms, opt = %s ms, par = %s ms" % \
                        (oldns*0.000001,newns*0.000001,optns*0.000001,parns*0.000001))
        print("\nENERGIES: old = %s, new = %s, opt = %s, par = %s" % \
                        (oldnrg,newnrg,optnrg,parnrg))

    assert_almost_equal( oldnrg, newnrg, 0.5 )
    assert_almost_equal( optnrg, newnrg, 0.1 )
    assert_almost_equal( parnrg, newnrg, 0.1 )

    if verbose:
        print("\nPerforming simulation...")

    moves.clearStatistics()
    moves.setGenerator( RanGenerator(42) )
    t.start()
    moves.move(oldsys, nmoves, False)
    oldns = t.nsecsElapsed()
    old_naccept = moves.nAccepted()
    old_nreject = moves.nRejected()

    oldcnrg = oldsys.energy( oldff.components().coulomb() ).value()
    oldljnrg = oldsys.energy( oldff.components().lj() ).value()

    oldsys.mustNowRecalculateFromScratch()    

    check_oldcnrg = oldsys.energy( oldff.components().coulomb() ).value()
    check_oldljnrg = oldsys.energy( oldff.components().lj() ).value()

    moves.clearStatistics()
    moves.setGenerator( RanGenerator(42) )
    t.start()
    moves.move(newsys, nmoves, False)
    newns = t.nsecsElapsed()
    new_naccept = moves.nAccepted()
    new_nreject = moves.nRejected()

    newcnrg = newsys.energy( newff.components().coulomb() ).value()
    newljnrg = newsys.energy( newff.components().lj() ).value()

    newsys.mustNowRecalculateFromScratch()    

    check_newcnrg = newsys.energy( newff.components().coulomb() ).value()
    check_newljnrg = newsys.energy( newff.components().lj() ).value()

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

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energy [as 别名]
        cljff_a.add(mol)
        cljff_a_b.add( mol, MGIdx(0) )
    else:
        cljff_b.add(mol)
        cljff_a_b.add( mol, MGIdx(1) )

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

system = System()
system.add(solvent)

system.add(cljff)

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

system2 = System()
system2.add(solvent)

system2.add(cljff_a)
system2.add(cljff_b)
system2.add(cljff_a_b)

print(system.groupNumbers())
print(system.groupNames())
print(system.energies())
print(system2.groupNumbers())

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

moves = SameMoves(rbmc)

PDB().write(grid_system.molecules(), "test0000.pdb")

t = QTime()

for i in range(1,11):
    print("Moving the system...")

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