Re: [chimera-dev] [Chimera-users] Amber Force Field
Hi Forbes,
I think Eric and Conrad may no more about your problem than I do. But I can tell you my best guess why your calculation is abysmally slow. I think the basic trouble is that useRotamer() decides to replace the existing residue rather than simply update its coordinates. In other words it deletes the old side chain atoms and bonds and makes new ones. Then the new ones don't have charges assigned for energy calculation. So the minimize "prep" option adds hydrogens and computes charges. It does this for the entire molecule because it doesn't know how to prep just one atom. So for each residue and each rotamer tried you end up adding hydrogens and recomputing charges for the whole molecule. Turning prep off in the minimize() call will just cause an error because the charges are missing. The "cache" option in the minimize() call also won't work because the molecule has changed. In fact the minimize code is not smart enough to realize the molecule has new atoms and will try to use the cached MMTK copy of the former molecule and this will probably lead to errors (as you describe).
So my 30 minute look at the various parts of the Chimera code without any testing suggests the slowness is because hydrogens/charges are recomputed for every rotamer. The new MMTK molecule object is just the one residue you are calculating the energy of. You realize your code does not take account of any interactions of that residue with any other residues. The other residues are completely ignored. That is what exclres means. So you may as well be asking for the energy of a one residue molecule starting in the rotamer configuration with the actual residue backbone atom positions.
I bet this would be very much faster if using a rotamer simply changed existing atom coordinates, then you add hydrogens and charges just one time for all the calculations, and you can cache the MMTK molecule object for each single residue energy calculation which will speed up calculating the different rotamer energies for that residue.
Eric is the one who would know why the useRotamer() code deletes atoms and adds new ones instead of just changing the coordinates. It would probably be pretty easy to make a new routine to use instead of useRotamer() that simply copied atom coordinates from the rotamer molecules returned by getRotamer() to the original molecule. Or what about just calculating the energy of the rotamer molecules directly without trying to put them into the protein (with all its other residues ignored)? I see when I use the Rotamer dialog and hide the original molecule that the returned rotamers include CA and N backbone atoms but not C, O.
Tom
Hi Tom:
Thanks for your email of May 4. After some experimenting I got the code to work but the execution times are deplorable! As I indicated earlier I want to use the Amber energy calculations to compare various confirmations involving different rotamers. My problem seems to be that I cannot get a fast execution because I am unable to cache the results of various computations and I believe the calculations are starting from scratch for each rotamer combination. Here is a simplified version of the code that I am using (in this case working with only one rotamer and essentaily computing its intrinsic energy (I hope)):
import chimera import Rotamers from Rotamers import getRotamers, useRotamer, NoResidueRotamersError from MMMD import base
def evaluateAmberEnergy(res_L, prot):
fixedAtoms = set() molecules = list([prot]) usingResidues = set(res_L) exclres = set(sum((m.residues for m in molecules), [])) exclres.difference_update(usingResidues) def run(minimizer): minimizer.run() minimizer = base.Minimizer(molecules, nsteps=0, stepsize=0.02, cgsteps=0, cgstepsize=0.02, interval = 10, fixedAtoms=fixedAtoms, nogui=True, addhyd=True, callback=run, exclres=exclres, cache=False, prep=True) energy = minimizer._mi.universe.energy() print "energy value: ", energy return energy#==========================================================================
# Global variable: prot = chimera.openModels.open("1CRN", type="PDB")[0]
for rc in prot.residues: rotCur_L = getRotamers(rc)[1] for rotC in rotCur_L: useRotamer(rc, [rotC]) energyAmber = evaluateAmberEnergy([rc], prot) print rc.type, rc.id.position, energyAmber
This is producing a lot of text in the Python Shell and when the program is run with a large protein it will go for hours...
Any suggestions for speeding this up would be greatly appreciated. I tried different settings of cache and prep but got Chimera errors such as missing C++ atom objects.
Depressing....
Best regards, Forbes
On Fri, 4 May 2012, Tom Goddard wrote:
Hi Forbes,
If you want to get the minimized energy within your Python scriptthen you need to run the minimization using a Python function call rather than a Chimera text command. To figure out how to do that I look at the Python code that is called by the Chimera text command in your Chimera distribution
chimera/share/MMMD/cmdline.pyor also available online
http://plato.cgl.ucsf.edu/trac/chimera/browser/trunk/libs/MMMD/cmdline.py
The minimize() Python function is called using the text command arguments. I see it basically does
from MMMD import base minimizer = base.Minimizer(molecules, nsteps, stepsize, cgsteps,cgstepsize, interval, fixedAtoms, memorize, nogui, addhyd, lambda m: m.run(), exclres, cache, prep)
Looking at base.py in the same directory and tracking back from where "Initial energy..." is printed in file MMTKinter.py I see the energy is obtained using
minimizer._mi.universe.energy()I have not tested this out but perhaps it will move you in the direction you want to go.
TomHi:
Thanks for the response. The code that you suggested seems to work. I am getting the initial energy reported in the Python Shell along with other lines that are reporting on the activities of the command.
Since I am embedding this command in a runCommand statement, I would like to get the initial energy value put into a variable within the Python script instead of the Python Shell window. Is there any way to accomplish this?
Cheers, Forbes
On Fri, 13 Apr 2012, Elaine Meng wrote:
Hi Forbes, The minimization tool reports the total energy, not any breakdown into components (electrostatic, VDW, strain) or interactions between sets of atoms. However, the "minimize" command has a "fragment" option for ignoring all but the specified residues, so you could at least limit what is considered the total system. That also makes the calculation faster, although single-point energy calculation (0 steps minimization) should already be pretty fast. Best, Elaine
Example (and see also "nogui true"):
minimize spec :13.a,17.a fragment true cgsteps 0 nsteps 0
http://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/minimize.html
On Apr 13, 2012, at 10:38 AM, Conrad Huang wrote:
Chimera's "Structure Editing -> Minimize Structure" tool (or the "minimize" command) may be used to do this. The tool reports the initial energy of the system in the Reply Log. If you give zero as the number of steps for both steepest descent and conjugate gradient minimization, then Chimera should return immediately.
Conrad
On 4/13/12 7:09 AM, Forbes J. Burkowski wrote:
Hi:
I have noticed that scripts can start with:
import MMTK import MMTK.ForceFields from MMTK.ForceFields import Amber94ForceField
We are doing some side chain packing studies that would benefit from the availablity of "easy" force field calculations. Currently, the algorithm works with useRotamers() to set two neighbouring rotamers and this is followed by a simple energy calculation to evaluate the interaction energy between the two rotamers. Right now, the energy calculation is only a bit more sophisticated than a steric collision detector.
Question: Can we use some functions in the MMTK.ForceFields module to get an improved energy calculation?
Any suggestions would be appreciated.
Cheers, Forbes _______________________________________________ Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Hi:
Thanks for taking the time to look at this. Yes, the code I sent you does not handle rotamer interactions. I have code to do that but I did not want you to spend a lot of time looking at it, because it is much longer. The code I did send only evaluates intrinsic energy for a single rotamer but there is still this long execution time.
The code that handles rotamer pairs seems to work but only for very small proteins like 1CRN. Even this takes several hours. I tried 1c9o. It ran for over 12 hours and then the OS stepped in and said that I ran out of memory! My machine has 12 GB memory with Intel Core i7 CPU 960 @ 3.2GHz (64 bit Windows 7 OS) so I am reasonably well off for computing power.
Hopefully, Conrad or Eric can provide some suggestions before I start making more extensive changes to the program...
Cheers, Forbes
On Fri, 11 May 2012, Tom Goddard wrote:
Hi Forbes,
I think Eric and Conrad may no more about your problem than I do. But I can tell you my best guess why your calculation is abysmally slow. I think the basic trouble is that useRotamer() decides to replace the existing residue rather than simply update its coordinates. In other words it deletes the old side chain atoms and bonds and makes new ones. Then the new ones don't have charges assigned for energy calculation. So the minimize "prep" option adds hydrogens and computes charges. It does this for the entire molecule because it doesn't know how to prep just one atom. So for each residue and each rotamer tried you end up adding hydrogens and recomputing charges for the whole molecule. Turning prep off in the minimize() call will just cause an error because the charges are missing. The "cache" option in the minimize() call also won't work because the molecule has changed. In fact the minimize code is not smart enough to realize the molecule has new atoms and will try to use the cached MMTK copy of the former molecule and this will probably lead to errors (as you describe).
So my 30 minute look at the various parts of the Chimera code without any testing suggests the slowness is because hydrogens/charges are recomputed for every rotamer. The new MMTK molecule object is just the one residue you are calculating the energy of. You realize your code does not take account of any interactions of that residue with any other residues. The other residues are completely ignored. That is what exclres means. So you may as well be asking for the energy of a one residue molecule starting in the rotamer configuration with the actual residue backbone atom positions.
I bet this would be very much faster if using a rotamer simply changed existing atom coordinates, then you add hydrogens and charges just one time for all the calculations, and you can cache the MMTK molecule object for each single residue energy calculation which will speed up calculating the different rotamer energies for that residue.
Eric is the one who would know why the useRotamer() code deletes atoms and adds new ones instead of just changing the coordinates. It would probably be pretty easy to make a new routine to use instead of useRotamer() that simply copied atom coordinates from the rotamer molecules returned by getRotamer() to the original molecule. Or what about just calculating the energy of the rotamer molecules directly without trying to put them into the protein (with all its other residues ignored)? I see when I use the Rotamer dialog and hide the original molecule that the returned rotamers include CA and N backbone atoms but not C, O.
TomHi Tom:
Thanks for your email of May 4. After some experimenting I got the code to work but the execution times are deplorable! As I indicated earlier I want to use the Amber energy calculations to compare various confirmations involving different rotamers. My problem seems to be that I cannot get a fast execution because I am unable to cache the results of various computations and I believe the calculations are starting from scratch for each rotamer combination. Here is a simplified version of the code that I am using (in this case working with only one rotamer and essentaily computing its intrinsic energy (I hope)):
import chimera import Rotamers from Rotamers import getRotamers, useRotamer, NoResidueRotamersError from MMMD import base
def evaluateAmberEnergy(res_L, prot):
fixedAtoms = set() molecules = list([prot]) usingResidues = set(res_L) exclres = set(sum((m.residues for m in molecules), [])) exclres.difference_update(usingResidues) def run(minimizer): minimizer.run() minimizer = base.Minimizer(molecules, nsteps=0, stepsize=0.02, cgsteps=0, cgstepsize=0.02, interval = 10, fixedAtoms=fixedAtoms, nogui=True, addhyd=True, callback=run, exclres=exclres, cache=False, prep=True) energy = minimizer._mi.universe.energy() print "energy value: ", energy return energy#==========================================================================
# Global variable: prot = chimera.openModels.open("1CRN", type="PDB")[0]
for rc in prot.residues: rotCur_L = getRotamers(rc)[1] for rotC in rotCur_L: useRotamer(rc, [rotC]) energyAmber = evaluateAmberEnergy([rc], prot) print rc.type, rc.id.position, energyAmber
This is producing a lot of text in the Python Shell and when the program is run with a large protein it will go for hours...
Any suggestions for speeding this up would be greatly appreciated. I tried different settings of cache and prep but got Chimera errors such as missing C++ atom objects.
Depressing....
Best regards, Forbes
On Fri, 4 May 2012, Tom Goddard wrote:
Hi Forbes,
If you want to get the minimized energy within your Python scriptthen you need to run the minimization using a Python function call rather than a Chimera text command. To figure out how to do that I look at the Python code that is called by the Chimera text command in your Chimera distribution
chimera/share/MMMD/cmdline.pyor also available online
http://plato.cgl.ucsf.edu/trac/chimera/browser/trunk/libs/MMMD/cmdline.py
The minimize() Python function is called using the text command arguments. I see it basically does
from MMMD import base minimizer = base.Minimizer(molecules, nsteps, stepsize, cgsteps,cgstepsize, interval, fixedAtoms, memorize, nogui, addhyd, lambda m: m.run(), exclres, cache, prep)
Looking at base.py in the same directory and tracking back from where "Initial energy..." is printed in file MMTKinter.py I see the energy is obtained using
minimizer._mi.universe.energy()I have not tested this out but perhaps it will move you in the direction you want to go.
TomHi:
Thanks for the response. The code that you suggested seems to work. I am getting the initial energy reported in the Python Shell along with other lines that are reporting on the activities of the command.
Since I am embedding this command in a runCommand statement, I would like to get the initial energy value put into a variable within the Python script instead of the Python Shell window. Is there any way to accomplish this?
Cheers, Forbes
On Fri, 13 Apr 2012, Elaine Meng wrote:
Hi Forbes, The minimization tool reports the total energy, not any breakdown into components (electrostatic, VDW, strain) or interactions between sets of atoms. However, the "minimize" command has a "fragment" option for ignoring all but the specified residues, so you could at least limit what is considered the total system. That also makes the calculation faster, although single-point energy calculation (0 steps minimization) should already be pretty fast. Best, Elaine
Example (and see also "nogui true"):
minimize spec :13.a,17.a fragment true cgsteps 0 nsteps 0
http://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/minimize.html
On Apr 13, 2012, at 10:38 AM, Conrad Huang wrote:
Chimera's "Structure Editing -> Minimize Structure" tool (or the "minimize" command) may be used to do this. The tool reports the initial energy of the system in the Reply Log. If you give zero as the number of steps for both steepest descent and conjugate gradient minimization, then Chimera should return immediately.
Conrad
On 4/13/12 7:09 AM, Forbes J. Burkowski wrote: > Hi: > > I have noticed that scripts can start with: > > import MMTK > import MMTK.ForceFields > from MMTK.ForceFields import Amber94ForceField > > We are doing some side chain packing studies that would benefit from the > availablity of "easy" force field calculations. Currently, the algorithm > works with useRotamers() to set two neighbouring rotamers and this is > followed by a simple energy calculation to evaluate the interaction energy > between the two rotamers. Right now, the energy calculation is only a bit > more sophisticated than a steric collision detector. > > Question: Can we use some functions in the MMTK.ForceFields module to get > an improved energy calculation? > > Any suggestions would be appreciated. > > Cheers, > Forbes > _______________________________________________ > Chimera-users mailing list > Chimera-users@cgl.ucsf.edu > http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users _______________________________________________ Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Hi Forbes/Tom, Keeping in mind that rotamers only specify heavy-atom positions, not hydrogen positions, the only scenario where changing side-chain coordinates would even be feasible is if:
A) the rotamer and original side chain are the same residue type B) the original side chain has the right atom names C) the original side chain has all the heavy atoms and no hydrogens D) no altloc funny business [including adding the rotamer as an altloc while retaining the original side chain]
...and even this wouldn't work for the energy-evaluation purpose because that needs hydrogens. I cannot come up with a fast way to do this evaluation that will produce values that mean anything. The closest thing I've got is to take the rotamers and instead of calling useRotamer(), for each rotamer take the chi angles and apply them to the existing side chain. This is trivial to do, but the hydrogens will still not be oriented in any intelligent way (they will have chemically reasonable distances/angles/dihedrals though, they just won't be smartly positioned for hydrogen bonding and clash avoidance). Perhaps after the chi angle assignments some steps of local steepest-descent minimization could be performed, though again quite possibly not fast enough. I don't know whether or not steepest descent is likely to get hydrogen bonding going, though it would certainly relieve clashes. The "trivial" chi-angle assignment code would look something like:
rotCur_L = getRotamers(rc)[1] for rotC in rotCur_L: for chi in range(1, 5): chiVal = getattr(rotC, "chi%d" % chi) if chiVal is None: break setattr(rc, "chi%d" % chi, chiVal)
--Eric
On May 11, 2012, at 11:58 AM, Tom Goddard wrote:
Hi Forbes,
I think Eric and Conrad may no more about your problem than I do. But I can tell you my best guess why your calculation is abysmally slow. I think the basic trouble is that useRotamer() decides to replace the existing residue rather than simply update its coordinates. In other words it deletes the old side chain atoms and bonds and makes new ones. Then the new ones don't have charges assigned for energy calculation. So the minimize "prep" option adds hydrogens and computes charges. It does this for the entire molecule because it doesn't know how to prep just one atom. So for each residue and each rotamer tried you end up adding hydrogens and recomputing charges for the whole molecule. Turning prep off in the minimize() call will just cause an error because the charges are missing. The "cache" option in the minimize() call also won't work because the molecule has changed. In fact the minimize code is not smart enough to realize the molecule has new atoms and will try to use the cached MMTK copy of the former molecule and this will probably lead to errors (as you describe).
So my 30 minute look at the various parts of the Chimera code without any testing suggests the slowness is because hydrogens/charges are recomputed for every rotamer. The new MMTK molecule object is just the one residue you are calculating the energy of. You realize your code does not take account of any interactions of that residue with any other residues. The other residues are completely ignored. That is what exclres means. So you may as well be asking for the energy of a one residue molecule starting in the rotamer configuration with the actual residue backbone atom positions.
I bet this would be very much faster if using a rotamer simply changed existing atom coordinates, then you add hydrogens and charges just one time for all the calculations, and you can cache the MMTK molecule object for each single residue energy calculation which will speed up calculating the different rotamer energies for that residue.
Eric is the one who would know why the useRotamer() code deletes atoms and adds new ones instead of just changing the coordinates. It would probably be pretty easy to make a new routine to use instead of useRotamer() that simply copied atom coordinates from the rotamer molecules returned by getRotamer() to the original molecule. Or what about just calculating the energy of the rotamer molecules directly without trying to put them into the protein (with all its other residues ignored)? I see when I use the Rotamer dialog and hide the original molecule that the returned rotamers include CA and N backbone atoms but not C, O.
Tom
Hi Tom:
Thanks for your email of May 4. After some experimenting I got the code to work but the execution times are deplorable! As I indicated earlier I want to use the Amber energy calculations to compare various confirmations involving different rotamers. My problem seems to be that I cannot get a fast execution because I am unable to cache the results of various computations and I believe the calculations are starting from scratch for each rotamer combination. Here is a simplified version of the code that I am using (in this case working with only one rotamer and essentaily computing its intrinsic energy (I hope)):
import chimera import Rotamers from Rotamers import getRotamers, useRotamer, NoResidueRotamersError from MMMD import base
def evaluateAmberEnergy(res_L, prot):
fixedAtoms = set() molecules = list([prot]) usingResidues = set(res_L) exclres = set(sum((m.residues for m in molecules), [])) exclres.difference_update(usingResidues) def run(minimizer): minimizer.run() minimizer = base.Minimizer(molecules, nsteps=0, stepsize=0.02, cgsteps=0, cgstepsize=0.02, interval = 10, fixedAtoms=fixedAtoms, nogui=True, addhyd=True, callback=run, exclres=exclres, cache=False, prep=True) energy = minimizer._mi.universe.energy() print "energy value: ", energy return energy#==========================================================================
# Global variable: prot = chimera.openModels.open("1CRN", type="PDB")[0]
for rc in prot.residues: rotCur_L = getRotamers(rc)[1] for rotC in rotCur_L: useRotamer(rc, [rotC]) energyAmber = evaluateAmberEnergy([rc], prot) print rc.type, rc.id.position, energyAmber
This is producing a lot of text in the Python Shell and when the program is run with a large protein it will go for hours...
Any suggestions for speeding this up would be greatly appreciated. I tried different settings of cache and prep but got Chimera errors such as missing C++ atom objects.
Depressing....
Best regards, Forbes
On Fri, 4 May 2012, Tom Goddard wrote:
Hi Forbes,
If you want to get the minimized energy within your Python script then you need to run the minimization using a Python function call rather than a Chimera text command. To figure out how to do that I look at the Python code that is called by the Chimera text command in your Chimera distribution
chimera/share/MMMD/cmdline.pyor also available online
http://plato.cgl.ucsf.edu/trac/chimera/browser/trunk/libs/MMMD/cmdline.py
The minimize() Python function is called using the text command arguments. I see it basically does
from MMMD import base minimizer = base.Minimizer(molecules, nsteps, stepsize, cgsteps,cgstepsize, interval, fixedAtoms, memorize, nogui, addhyd, lambda m: m.run(), exclres, cache, prep)
Looking at base.py in the same directory and tracking back from where "Initial energy..." is printed in file MMTKinter.py I see the energy is obtained using
minimizer._mi.universe.energy()I have not tested this out but perhaps it will move you in the direction you want to go.
TomHi:
Thanks for the response. The code that you suggested seems to work. I am getting the initial energy reported in the Python Shell along with other lines that are reporting on the activities of the command.
Since I am embedding this command in a runCommand statement, I would like to get the initial energy value put into a variable within the Python script instead of the Python Shell window. Is there any way to accomplish this?
Cheers, Forbes
On Fri, 13 Apr 2012, Elaine Meng wrote:
Hi Forbes, The minimization tool reports the total energy, not any breakdown into components (electrostatic, VDW, strain) or interactions between sets of atoms. However, the "minimize" command has a "fragment" option for ignoring all but the specified residues, so you could at least limit what is considered the total system. That also makes the calculation faster, although single-point energy calculation (0 steps minimization) should already be pretty fast. Best, Elaine
Example (and see also "nogui true"):
minimize spec :13.a,17.a fragment true cgsteps 0 nsteps 0
http://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/minimize.html
On Apr 13, 2012, at 10:38 AM, Conrad Huang wrote:
Chimera's "Structure Editing -> Minimize Structure" tool (or the "minimize" command) may be used to do this. The tool reports the initial energy of the system in the Reply Log. If you give zero as the number of steps for both steepest descent and conjugate gradient minimization, then Chimera should return immediately.
Conrad
On 4/13/12 7:09 AM, Forbes J. Burkowski wrote: > Hi: > > I have noticed that scripts can start with: > > import MMTK > import MMTK.ForceFields > from MMTK.ForceFields import Amber94ForceField > > We are doing some side chain packing studies that would benefit from the > availablity of "easy" force field calculations. Currently, the algorithm > works with useRotamers() to set two neighbouring rotamers and this is > followed by a simple energy calculation to evaluate the interaction energy > between the two rotamers. Right now, the energy calculation is only a bit > more sophisticated than a steric collision detector. > > Question: Can we use some functions in the MMTK.ForceFields module to get > an improved energy calculation? > > Any suggestions would be appreciated. > > Cheers, > Forbes > _______________________________________________ > Chimera-users mailing list > Chimera-users@cgl.ucsf.edu > http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users _______________________________________________ Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Chimera-dev mailing list Chimera-dev@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-dev
Hi Eric,
That is useful. The key point is that rotamers assigned by Chimera don't include hydrogen orientations. I believe the Chimera minimize code which is calculating the energy requires hydrogens. The prep stage of each minimize calculation adds hydrogens including refining their positions to some degree. Forbes will have to address the question of how he wants the hydrogens positioned for his rotamer investigation. I think minimizing just one residue or just the hydrogens on one residue with heavy atoms fixed would likely be much faster then running prep for each new rotamer. But testing would be needed to verify that. And maybe optimizing the nearby hydrogens would be sensible too. I don't think prep moves already placed hydrogens.
Tom
Hi Forbes/Tom, Keeping in mind that rotamers only specify heavy-atom positions, not hydrogen positions, the only scenario where changing side-chain coordinates would even be feasible is if:
A) the rotamer and original side chain are the same residue type B) the original side chain has the right atom names C) the original side chain has all the heavy atoms and no hydrogens D) no altloc funny business [including adding the rotamer as an altloc while retaining the original side chain]
...and even this wouldn't work for the energy-evaluation purpose because that needs hydrogens. I cannot come up with a fast way to do this evaluation that will produce values that mean anything. The closest thing I've got is to take the rotamers and instead of calling useRotamer(), for each rotamer take the chi angles and apply them to the existing side chain. This is trivial to do, but the hydrogens will still not be oriented in any intelligent way (they will have chemically reasonable distances/angles/dihedrals though, they just won't be smartly positioned for hydrogen bonding and clash avoidance). Perhaps after the chi angle assignments some steps of local steepest-descent minimization could be performed, though again quite possibly not fast enough. I don't know whether or not steepest descent is likely to get hydrogen bonding going, though it would certainly relieve clashes. The "trivial" chi-angle assignment code would look something like:
rotCur_L = getRotamers(rc)[1] for rotC in rotCur_L: for chi in range(1, 5): chiVal = getattr(rotC, "chi%d" % chi) if chiVal is None: break setattr(rc, "chi%d" % chi, chiVal)
--Eric
On May 11, 2012, at 11:58 AM, Tom Goddard wrote:
Hi Forbes,
I think Eric and Conrad may no more about your problem than I do. But I can tell you my best guess why your calculation is abysmally slow. I think the basic trouble is that useRotamer() decides to replace the existing residue rather than simply update its coordinates. In other words it deletes the old side chain atoms and bonds and makes new ones. Then the new ones don't have charges assigned for energy calculation. So the minimize "prep" option adds hydrogens and computes charges. It does this for the entire molecule because it doesn't know how to prep just one atom. So for each residue and each rotamer tried you end up adding hydrogens and recomputing charges for the whole molecule. Turning prep off in the minimize() call will just cause an error because the charges are missing. The "cache" option in the minimize() call also won't work because the molecule has changed. In fact the minimize code is not smart enough to realize the molecule has new atoms and will try to use the cached MMTK copy of the former molecule and this will probably lead to errors (as you describe).
So my 30 minute look at the various parts of the Chimera code without any testing suggests the slowness is because hydrogens/charges are recomputed for every rotamer. The new MMTK molecule object is just the one residue you are calculating the energy of. You realize your code does not take account of any interactions of that residue with any other residues. The other residues are completely ignored. That is what exclres means. So you may as well be asking for the energy of a one residue molecule starting in the rotamer configuration with the actual residue backbone atom positions.
I bet this would be very much faster if using a rotamer simply changed existing atom coordinates, then you add hydrogens and charges just one time for all the calculations, and you can cache the MMTK molecule object for each single residue energy calculation which will speed up calculating the different rotamer energies for that residue.
Eric is the one who would know why the useRotamer() code deletes atoms and adds new ones instead of just changing the coordinates. It would probably be pretty easy to make a new routine to use instead of useRotamer() that simply copied atom coordinates from the rotamer molecules returned by getRotamer() to the original molecule. Or what about just calculating the energy of the rotamer molecules directly without trying to put them into the protein (with all its other residues ignored)? I see when I use the Rotamer dialog and hide the original molecule that the returned rotamers include CA and N backbone atoms but not C, O.
TomHi Tom:
Thanks for your email of May 4. After some experimenting I got the code to work but the execution times are deplorable! As I indicated earlier I want to use the Amber energy calculations to compare various confirmations involving different rotamers. My problem seems to be that I cannot get a fast execution because I am unable to cache the results of various computations and I believe the calculations are starting from scratch for each rotamer combination. Here is a simplified version of the code that I am using (in this case working with only one rotamer and essentaily computing its intrinsic energy (I hope)):
import chimera import Rotamers from Rotamers import getRotamers, useRotamer, NoResidueRotamersError from MMMD import base
def evaluateAmberEnergy(res_L, prot):
fixedAtoms = set() molecules = list([prot]) usingResidues = set(res_L) exclres = set(sum((m.residues for m in molecules), [])) exclres.difference_update(usingResidues) def run(minimizer): minimizer.run() minimizer = base.Minimizer(molecules, nsteps=0, stepsize=0.02, cgsteps=0, cgstepsize=0.02, interval = 10, fixedAtoms=fixedAtoms, nogui=True, addhyd=True, callback=run, exclres=exclres, cache=False, prep=True) energy = minimizer._mi.universe.energy() print "energy value: ", energy return energy#==========================================================================
# Global variable: prot = chimera.openModels.open("1CRN", type="PDB")[0]
for rc in prot.residues: rotCur_L = getRotamers(rc)[1] for rotC in rotCur_L: useRotamer(rc, [rotC]) energyAmber = evaluateAmberEnergy([rc], prot) print rc.type, rc.id.position, energyAmber
This is producing a lot of text in the Python Shell and when the program is run with a large protein it will go for hours...
Any suggestions for speeding this up would be greatly appreciated. I tried different settings of cache and prep but got Chimera errors such as missing C++ atom objects.
Depressing....
Best regards, Forbes
On Fri, 4 May 2012, Tom Goddard wrote:
Hi Forbes,
If you want to get the minimized energy within your Python scriptthen you need to run the minimization using a Python function call rather than a Chimera text command. To figure out how to do that I look at the Python code that is called by the Chimera text command in your Chimera distribution
chimera/share/MMMD/cmdline.pyor also available online
http://plato.cgl.ucsf.edu/trac/chimera/browser/trunk/libs/MMMD/cmdline.py
The minimize() Python function is called using the text command arguments. I see it basically does
from MMMD import base minimizer = base.Minimizer(molecules, nsteps, stepsize, cgsteps,cgstepsize, interval, fixedAtoms, memorize, nogui, addhyd, lambda m: m.run(), exclres, cache, prep)
Looking at base.py in the same directory and tracking back from where "Initial energy..." is printed in file MMTKinter.py I see the energy is obtained using
minimizer._mi.universe.energy()I have not tested this out but perhaps it will move you in the direction you want to go.
TomHi:
Thanks for the response. The code that you suggested seems to work. I am getting the initial energy reported in the Python Shell along with other lines that are reporting on the activities of the command.
Since I am embedding this command in a runCommand statement, I would like to get the initial energy value put into a variable within the Python script instead of the Python Shell window. Is there any way to accomplish this?
Cheers, Forbes
On Fri, 13 Apr 2012, Elaine Meng wrote:
Hi Forbes, The minimization tool reports the total energy, not any breakdown into components (electrostatic, VDW, strain) or interactions between sets of atoms. However, the "minimize" command has a "fragment" option for ignoring all but the specified residues, so you could at least limit what is considered the total system. That also makes the calculation faster, although single-point energy calculation (0 steps minimization) should already be pretty fast. Best, Elaine
Example (and see also "nogui true"):
minimize spec :13.a,17.a fragment true cgsteps 0 nsteps 0
http://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/minimize.html
On Apr 13, 2012, at 10:38 AM, Conrad Huang wrote:
> Chimera's "Structure Editing -> Minimize Structure" tool (or the "minimize" command) may be used to do this. The tool reports the initial energy of the system in the Reply Log. If you give zero as the number of steps for both steepest descent and conjugate gradient minimization, then Chimera should return immediately. > > Conrad > > On 4/13/12 7:09 AM, Forbes J. Burkowski wrote: >> Hi: >> >> I have noticed that scripts can start with: >> >> import MMTK >> import MMTK.ForceFields >> from MMTK.ForceFields import Amber94ForceField >> >> We are doing some side chain packing studies that would benefit from the >> availablity of "easy" force field calculations. Currently, the algorithm >> works with useRotamers() to set two neighbouring rotamers and this is >> followed by a simple energy calculation to evaluate the interaction energy >> between the two rotamers. Right now, the energy calculation is only a bit >> more sophisticated than a steric collision detector. >> >> Question: Can we use some functions in the MMTK.ForceFields module to get >> an improved energy calculation? >> >> Any suggestions would be appreciated. >> >> Cheers, >> Forbes >> _______________________________________________ >> Chimera-users mailing list >> Chimera-users@cgl.ucsf.edu >> http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users > _______________________________________________ > Chimera-users mailing list > Chimera-users@cgl.ucsf.edu > http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Chimera-users mailing list Chimera-users@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-users
Chimera-dev mailing list Chimera-dev@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-dev
Hi all, I have a problem associating a vrml object an a molecule : Let's say I have an arrow created using the bild format to which I gave 0.0 0.0 0.0 coordinates. This arrow and the molecule translate/rotate together using the mouse. Now if I translate / rotate my molecule using the setcoord command on the atoms in a python script, the arrow won't follow up the molecule. I tried the addAssociatedModel command to link the arrow and the molecule but still the same result. How can I handle this? Thanks for your attention
The arrow you drew has no idea that the atoms have moved, so of course it doesn't move with them. If you just want to rigidly rotate the whole molecule in Python you can change the matrix that places it in space
mol.openState.xform = xf
instead of changing the coordinates of atoms. Then you could rigidly move the VRML model in the same way so that the arrow stays in the same position relative to the molecule.
vrml.openState.xform = xf
If you are not moving the molecule rigidly then you really do need to move individual atoms and then you will have to delete the old arrow and make a new arrow if it is intended to have its ends follow certain atoms.
Tom
-------- Original Message -------- Subject: [chimera-dev] Moving a residue and an VRML object From: ladam To: chimera-dev@cgl.ucsf.edu Date: 10/17/12 4:45 AM
Hi all, I have a problem associating a vrml object an a molecule : Let's say I have an arrow created using the bild format to which I gave 0.0 0.0 0.0 coordinates. This arrow and the molecule translate/rotate together using the mouse. Now if I translate / rotate my molecule using the setcoord command on the atoms in a python script, the arrow won't follow up the molecule. I tried the addAssociatedModel command to link the arrow and the molecule but still the same result. How can I handle this? Thanks for your attention
Salut,
The other trick is to give the VRML model the same model id and submodel id as its associated model. Then the openState is identical, so there would be no need to copy the transformation matrix. Again, this only works for rigid rotations.
-- Greg
On 10/17/2012 01:25 PM, Tom Goddard wrote:
The arrow you drew has no idea that the atoms have moved, so of course it doesn't move with them. If you just want to rigidly rotate the whole molecule in Python you can change the matrix that places it in space
mol.openState.xform = xfinstead of changing the coordinates of atoms. Then you could rigidly move the VRML model in the same way so that the arrow stays in the same position relative to the molecule.
vrml.openState.xform = xfIf you are not moving the molecule rigidly then you really do need to move individual atoms and then you will have to delete the old arrow and make a new arrow if it is intended to have its ends follow certain atoms.
Tom-------- Original Message -------- Subject: [chimera-dev] Moving a residue and an VRML object From: ladam To: chimera-dev@cgl.ucsf.edu Date: 10/17/12 4:45 AM
Hi all, I have a problem associating a vrml object an a molecule : Let's say I have an arrow created using the bild format to which I gave 0.0 0.0 0.0 coordinates. This arrow and the molecule translate/rotate together using the mouse. Now if I translate / rotate my molecule using the setcoord command on the atoms in a python script, the arrow won't follow up the molecule. I tried the addAssociatedModel command to link the arrow and the molecule but still the same result. How can I handle this? Thanks for your attention
Chimera-dev mailing list Chimera-dev@cgl.ucsf.edu http://plato.cgl.ucsf.edu/mailman/listinfo/chimera-dev
On Oct 17, 2012, at 3:21 PM, Greg Couch wrote:
Salut,
The other trick is to give the VRML model the same model id and submodel id as its associated model. Then the openState is identical, so there would be no need to copy the transformation matrix. Again, this only works for rigid rotations.
...and the simplest way to do this is to supply the "sameAs" keyword, with the associated model as the argument, to the openModels.add() or openModels.open() call that you use.
--Eric
Eric Pettersen UCSF Computer Graphics Lab http://www.cgl.ucsf.edu
participants (5)
-
Eric Pettersen -
Forbes J. Burkowski -
Greg Couch -
ladam -
Tom Goddard