Get the total energy of the input conformation

[diallobakary4]

Hi all Is there a way to get the total energy of the input conformation without doing any sampling? In the findgmec script for example, where can we get the total energy of the structure 1CC8.ss.pdb used as input?

[cuchaz]

Hello,

Strangely enough, that isn't an operation we do very often. I think something like this might work for you though:

pmol = osprey.c.confspace.ParametricMolecule(strand.mol)
inters = osprey.c.energy.ResidueInteractions()
inters.addComplete(pmol.mol.residues)
energy = ecalc.calcEnergy(pmol, inters).energy
print(energy)

That should print the rigid (unminimized) energy of the strand in the strand variable using the energy calculator in the ecalc variable.

[diallobakary4]

Thank you. This works. Does this take into account the reference energy? Is there a way to do that?

[cuchaz]

No, that technique doesn't account for reference energy. For that, you'll have to pick a conformation from your conformation space (probably you want the wild-type conformation in this case) and the calculate its energy using the higher-level ConfEnergyCalculator.

[diallobakary4]

pmol = confEcalc.confSpace.makeMolecule([0])
inters = osprey.c.energy.ResidueInteractions()
inters.addComplete(conf.mol.residues)
print(confEcalc.ecalc.calcEnergy(pmol, inters).energy)

Using the above? Not sure, but I am assuming [0] to correspond to the wild type? Also not sure if the final line is the right one. The end value is sill different from the first approach. print(confEcalc.calcEnergy(pmol, inters).energy) did not work.

[cuchaz]

That looks really close. The tricky part is getting the right numbers for the conformation. Like I said before, this isn't something we do very often (weird, right?) so there aren't any nice shortcuts.

The closest thing I can find is some testing code (in Java) that shows the logic for building a conformation from only the "wild-type rotamers":

https://github.com/donaldlab/OSPREY3/blob/0d2520b09b90dbd051c65b9de002ceb1d792f87e/src/test/java/edu/duke/cs/osprey/confspace/compiled/TestConfSpace.java#L40

That should help you get the conformation you're looking for.

EDIT: oh, right. For this to work, you'll need to make sure you explicitly include the wild-type rotamers in your conformation space.

[cuchaz]

One more thing, that makeMolecule call will generate continuous degrees of freedom for your molecule, if you configured that in your conf space. If you still want rigid energies here, you'll need to tell the energy calculator to turn off minimization. See here in the python API:

https://github.com/donaldlab/OSPREY3/blob/0d2520b09b90dbd051c65b9de002ceb1d792f87e/src/main/python/osprey/__init__.py#L587

[diallobakary4]

Thank you for the link to the java code

(weird, right?) Yeap that the least we can say :)

From the snippet below, can we assume that -1 corresponds to the wild-type? The output looks like this:

pmol = osprey.c.confspace.ParametricMolecule(strand.mol)
inters = osprey.c.energy.ResidueInteractions()
inters.addComplete(pmol.mol.residues)
energy = ecalc.calcEnergy(pmol, inters).energy
print("Input conf", energy)

# is -1 the number for the input conf?
# as it shows same energy as the one above
pmol = confEcalc.confSpace.makeMolecule([-1])
inters = osprey.c.energy.ResidueInteractions()
inters.addComplete(conf.mol.residues)
print(-1, confEcalc.ecalc.calcEnergy(pmol, inters).energy)

for x in range(20):
    pmol = confEcalc.confSpace.makeMolecule([x])
    inters = osprey.c.energy.ResidueInteractions()
    inters.addComplete(conf.mol.residues)
    print(x, confEcalc.ecalc.calcEnergy(pmol, inters).energy)

[cuchaz]

I had to fiddle with this one for a bit to figure it out. You've found an interesting coincidence in Osprey's energy system, but sadly that code snippet won't do quite what you want.

We usually use the -1 index to mean "not assigned" for a design position. So making a conformation of just [-1] will produce a wild-type molecule without any mutations. Which is actually what we want in this case, but it's not really using the conformation system, so we won't be able to calculate conformation energies using reference energies that way.

Wild-type residue conformations don't get assigned a canonical index, like -1 or anything else. They do get assigned an index, but you have to search to find out what it is after you've created the conformation space. Once you have the wild-type conformation indices, then you can use the conformation energy calculator (confEcalc) to compute the conformation energies with all the bells and whistles, like reference energies.

There's one more catch though. Since you're looking for the total molecule energy (instead of just the energy of the design positions, which is what Osprey calculates by default to save time), you'll have to create custom residue interactions to give to the conformation energy calculator. The function that creates those residue interactions allows you to pass in the reference energies.

Clear as mud, right?

Well, that's the high-level overview. I've made a new example script to walk through all the low-level details. Hopefully that can get you the energies you're looking for:

https://github.com/donaldlab/OSPREY3/blob/main/examples/python.GMEC/wildtypeEnergies.py