[chimera-dev] [Chimera-users] Amber Force Field

[Tom Goddard]

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.py
>>
>> or 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.
>>
>>       Tom
>>
>>
>>> Hi:
>>>
>>> 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 at cgl.ucsf.edu
>>>>>> http://www.rbvi.ucsf.edu/mailman/listinfo/chimera-users
>>>>> _______________________________________________
>>>>> Chimera-users mailing list
>>>>> Chimera-users at cgl.ucsf.edu
>>>>> http://www.rbvi.ucsf.edu/mailman/listinfo/chimera-users
>>> _______________________________________________
>>> Chimera-users mailing list
>>> Chimera-users at cgl.ucsf.edu
>>> http://www.rbvi.ucsf.edu/mailman/listinfo/chimera-users
>>>
>>