Python System.energies方法代码示例

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

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

    ligand_intraff = InternalFF("ligand:intra")
    ligand_intraff.add(ligand_mols)

    forcefields.append(ligand_intraclj)
    forcefields.append(ligand_intraff)

    ligand_mm_nrg = ligand_intraclj.components().total() + ligand_intraff.components().total()

    ###
    ### FORCEFIELDS INVOLVING THE LIGAND/CLUSTER AND OTHER ATOMS
    ###

    # forcefield holding the energy between the ligand and the mobile atoms in the
    # bound leg
    ligand_mobile = InterGroupCLJFF("system:ligand-mobile")
    ligand_mobile = setCLJProperties(ligand_mobile, space)

    ligand_mobile.add(ligand_mols, MGIdx(0))
    ligand_mobile.add(mobile_mols, MGIdx(1))

    qm_ligand = QMMMFF("system:ligand-QM")    
    qm_ligand = setQMProperties(qm_ligand, space)

    qm_ligand.add(ligand_mols, MGIdx(0))

    zero_energy = 0

    if not intermolecular_only.val:
        if qm_zero_energy.val is None:
            # calculate the delta value for the system - this is the difference between
            # the MM and QM intramolecular energy of the ligand
            t.start()
            print("\nComparing the MM and QM energies of the ligand...")
            mm_intra = ligand_intraclj.energy().value() + ligand_intraff.energy().value()
            print("MM energy = %s kcal mol-1 (took %s ms)" % (mm_intra, t.elapsed()))

            t.start()
            qm_intra = qm_ligand.energy().value()
            print("QM energy = %s kcal mol-1 (took %s ms)" % (qm_intra, t.elapsed()))

            print("\nSetting the QM zero energy to %s kcal mol-1" % (qm_intra - mm_intra))
            qm_ligand.setZeroEnergy( (qm_intra-mm_intra) * kcal_per_mol )
            zero_energy = qm_intra - mm_intra
        else:
            print("\nManually setting the QM zero energy to %s" % qm_zero_energy.val)
            qm_ligand.setZeroEnergy( qm_zero_energy.val )
            zero_energy = qm_zero_energy.val

    qm_ligand.add(mobile_mols, MGIdx(1))

    ligand_mm_nrg += ligand_mobile.components().total()
    ligand_qm_nrg = qm_ligand.components().total() + ligand_mobile.components().lj()

    if intermolecular_only.val:
        # the QM model still uses the MM intramolecular energy of the ligand
        ligand_qm_nrg += ligand_intraclj.components().total() + ligand_intraff.components().total()

    forcefields.append(ligand_mobile)
    forcefields.append(qm_ligand)

    if fixed_group.nMolecules() > 0:
        # there are fixed molecules

        # Whether or not to disable the grid and calculate all energies atomisticly
        if disable_grid:

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energies [as 别名]
    waters.update(water)

print("Constructing the forcefields...")

pol_tip4p = waters.moleculeAt(0).molecule()
waters.remove(pol_tip4p)

cljff = InterGroupCLJFF("pol_tip4p-water")
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"), \

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energies [as 别名]
move.setMaximumRotation(5 * degrees)

# We now group all of the moves to be performed into a single Moves
# object. In this case, a WeightedMoves will draw moves randomly
# according to their weight
moves = WeightedMoves()
moves.add( move, 1 )

# Lets perform 100 moves. The moves are performed on a copy of 'system',
# with the updated version of 'system' after the moves returned by this
# function
print("Running 100 moves...")
new_system = moves.move(system, 100, True)

# Now lets run a simulation, writing out a PDB trajectory
PDB().write(system.molecules(), "output000.pdb")
print(system.energies())

# Here we run 10 blocks of 1000 moves, printing out the 
# energies and a PDB of the coordinates after each block
for i in range(1,11):
    system = moves.move(system, 1000, True)
    print("%d: %s" % (i, system.energies()))
    PDB().write(system.molecules(), "output%003d.pdb" % i)

# Finally, we print out information about how many moves
# were accepted and rejected.
print("Move information:")
print(moves)

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

    oldsys = System()
    newsys = System()

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

    oldsys.add(oldff)
    newsys.add(newff)

    t = QElapsedTimer()

    oldsys.mustNowRecalculateFromScratch()
    newsys.mustNowRecalculateFromScratch()

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

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

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

    newcnrg = newsys.energy( newff.components().coulomb() ).value()
    newljnrg = newsys.energy( newff.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))

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

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

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

    moves.setGenerator( RanGenerator( 42 ) )
    moves.clearStatistics()

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

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

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

    oldcnrg = oldsys.energy( oldff.components().coulomb() ).value()
    oldljnrg = oldsys.energy( oldff.components().lj() ).value()
    
    newcnrg = newsys.energy( newff.components().coulomb() ).value()
    newljnrg = newsys.energy( newff.components().lj() ).value()

    if verbose:
        print("\nMoves: old %s ms vs. new %s ms" % (0.000001*move_oldns, 0.000001*move_newns))
        print("OLD SYS:  %s  %s  %s  : %s ms" % (oldcnrg+oldljnrg,oldcnrg,oldljnrg,
                                                 0.000001*oldns))
        print("nAccepted() = %s, nRejected() = %s" % (old_naccepted, old_nrejected))
        print("NEW SYS:  %s  %s  %s  : %s ms" % (newcnrg+newljnrg,newcnrg,newljnrg,
                                                 0.000001*newns))
        print("nAccepted() = %s, nRejected() = %s" % (new_naccepted, new_nrejected))
    
    oldsys.mustNowRecalculateFromScratch()
    newsys.mustNowRecalculateFromScratch()
    
    t.start()
    nrgs = oldsys.energies()
    oldns = t.nsecsElapsed()
    
    t.start()
    nrgs = newsys.energies()
    newns = t.nsecsElapsed()

    r_oldcnrg = oldsys.energy( oldff.components().coulomb() ).value()
    r_oldljnrg = oldsys.energy( oldff.components().lj() ).value()

    r_newcnrg = newsys.energy( newff.components().coulomb() ).value()
    r_newljnrg = newsys.energy( newff.components().lj() ).value()

    if verbose:
        print("\nRecalculated energy")
#.........这里部分代码省略.........

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

moves = SameMoves(rbmc)

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

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

moves = SameMoves(rbmc)

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

# 需要导入模块: from Sire import System [as 别名]
# 或者: from Sire.System import energies [as 别名]
system.add(solvent)
system.setProperty("space", space)

print("Initialising the ions...")
titrator.applyTo(system)
print("Randomising the location of ions...")
titrator.randomiseCharge(3)
titrator.applyTo(system)
print("System is ready for simulation :-)")

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

move = TitrationMove()
move.setTemperature( 25*celsius )

moves = WeightedMoves()
moves.add(move, 1)

move = RigidBodyMC(solvent)
moves.add(move, 1)

print("Start: %s" % system.energies())

for i in range(1,11):
    system = moves.move(system, 1000, False)
    print("%5d: %s" % (i, system.energies()))
    PDB().write(system.molecules(), "test%0004d.pdb" % i)

print(moves)