Print bond energies

[Stefano Bosisio]

Good morning,

I would like to know if there's a method for printing all the bond energies for every bond I have in a molecole.
For example in methane CH4 I would like to print bond energy of every CH I have.
I found that with (Sire rev 2761 OpenMM 6.1)
system.molecules().first().molecule().property('bond').potentials()
I can see the potential formula. But I'd rather want the numerical value.
What can I do?
Thanks :)

[julien]

Hi Stefano, 

Try to adapt the following code to your problem 

mol = system.molecules().first().molecule()

bonds = mol.property("bond")

bond_potentials = bonds.potentials()

space = system.property("space")

for bond_potential in bond_potentials:

    #print (bond_potential)

    at0 = bond_potential.atom0() # grab reference to at0 and 1

    at1 = bond_potential.atom1()

    at0 = mol.select(at0)               # select at0 and at1 so can pull other attributes

    at1 = mol.select(at1)

    at0coords = at0.property("coordinates")

    at1coords = at1.property("coordinates")

    r01 = space.calcDist(at0coords, at1coords) # use space object in case the coordinates should be wrapped

    #print (r01)

    values = Values({Symbol("r"):r01}) # define a symbol r and associate numerical value of r01

    nrg = bond_potential.function().evaluate(values) * kcal_per_mol # evaluate the potential function with the value of r set to r01

    print ("Bond %s-%s : %s " % (at0.name().value(),at1.name().value(),nrg) )

Best wishes, 

Julien

[Stefano Bosisio]

Great Julien!

It works and I can adapt for my purpose

Thanks a lot!


Best wishes,

Stefano