MinimizeAProtein

How can I energetically optimize the atoms of a protein?

BALL offers a class Minimizer to optimizes atom positions with respect to a ForceField.

C++

#include <BALL/FORMAT/PDBFile.h>
#include <BALL/MOLMEC/AMBER/amber.h>
#include <BALL/MOLMEC/MINIMIZATION/conjugateGradient.h>
#include <BALL/MOLMEC/COMMON/assignTypes.h>
#include <BALL/KERNEL/selector.h>
#include <BALL/STRUCTURE/defaultProcessors.h>
#include <BALL/STRUCTURE/residueChecker.h>
#include <BALL/STRUCTURE/fragmentDB.h>

using namespace BALL;
using namespace std;

int main(int argc, char** argv)
{
  // open a PDB file with the name of the first argument
  PDBFile file(argv[1]);
      
  // create a system and read the contents of the PDB file
  System S;
  file >> S;
  file.close();

  cout << "read " << S.countAtoms() << " atoms." << endl;

  // prepare the protein by adding missing information
  FragmentDB fragment_db("fragments/Fragments.db");
  S.apply(fragment_db.normalize_names);
  S.apply(fragment_db.add_hydrogens);  
  S.apply(fragment_db.build_bonds);
  ResidueChecker checker(fragment_db);
  S.apply(checker);
 
  // create a forcefield
  AmberFF FF(S);

  // just for curiosity: check how many atoms we are going
  // to optimize
  cout << "optimizing " << FF.getNumberOfMovableAtoms() << " out of " << S.countAtoms() << " atoms" << endl;
        
  // now we create a minimizer object that uses a conjugate 
  // gradient algorithm to optimize the atom positions
  ConjugateGradientMinimizer minimizer;

  // calculate the total energy of the system
  float initial_energy = FF.updateEnergy();
  cout << "initial energy: " << initial_energy << " kJ/mol" << endl;

  // initialize the minimizer and perform (up to)
  // 1000 optimization steps
  minimizer.setup(FF);
  minimizer.setEnergyOutputFrequency(1);
  minimizer.minimize(50);

  // calculate the terminal energy and print it
  float terminal_energy = FF.getEnergy();

  cout << "energy before/after minimization: " << initial_energy << "/" << terminal_energy << " kJ/mol" << endl;

  // write the optimized structure to a file whose
  // name is given as the second command line argument
  file.open(argv[2], ios::out);
  file << S;
  file.close();

  // done.
  return 0;
}

Note: The minimizer can be restricted by a selection, e.g. hydrogen atoms by adding the following lines after the instantiation of the forcefield:

 
S.deselect();
FF.setup(S);
Selector selector("element(H)");
S.apply(selector);

Python

#S = System()
#pdbfile.read(S)

S = getSystems()[0]

fdb = FragmentDB("fragments/Fragments.db")
S.apply(fdb.add_hydrogens)
S.apply(fdb.build_bonds)
S.apply(fdb.normalize_names)

# setup a force field
S.deselect()
FF = AmberFF()
FF.setup(S)

# restrict to some atoms, hydrogens or certain residues
selector = Selector("residue(ASP) OR name(H)")
S.apply(selector)
print "Number of selected atoms:", len(selector.getSelectedAtoms())
print "Optimizing ", FF.getNumberOfMovableAtoms(), " out of ", S.countAtoms(), " atoms"

# now we create a minimizer object that uses a conjugate 
# gradient algorithm to optimize the atom positions
minimizer = ConjugateGradientMinimizer()

# calculate the total energy of the system
initial_energy = FF.updateEnergy()

print "Initial energy: ", initial_energy, " kJ/mol"
# initialize the minimizer and perform (up to)
# 50 optimization steps

minimizer.setup(FF)
minimizer.setEnergyOutputFrequency(1)
minimizer.minimize(50)

# calculate the terminal energy and print it
terminal_energy = FF.getEnergy()
print "Energy before/after minimization: ", initial_energy, "/", terminal_energy, " kJ/mol"

# write the optimized structure to a file
outname = "minimized_file.pdb"

print "Write resulting structure as ", outname, " to disk"
outfile = PDBFile()
outfile.open(outname, OpenMode(File.MODE_OUT))
outfile << S
outfile.close()

# initiate an update of the BALLView scene 
getMainControl().update(S)