AmberEnergyEvaluations

How can I compute the Amber energy of a system?

BALL offers an implementation of the Amber forcefield.

C++

#!cpp
#include <BALL/FORMAT/PDBFile.h>
#include <BALL/MOLMEC/AMBER/amber.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("");
  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 amber;
  
  // set some options
  amber.options.set(AmberFF::Option::OVERWRITE_TYPENAMES, "true");
  amber.options.set(AmberFF::Option::ASSIGN_TYPES, "true");
  amber.options.set(AmberFF::Option::ASSIGN_CHARGES, "true");
  amber.options.set(AmberFF::Option::ASSIGN_TYPENAMES, "true");
  
  // add all atoms of system S to the forcefield
  amber.setup(S);

  // just for curiosity: check how many atoms we are going
  // to optimize
  cout << "optimizing " << amber.getNumberOfMovableAtoms() << " out of " << S.countAtoms() << " atoms" << endl;

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

  // manipulate the system, e.g. transformations

  // calculate the new energy and print it
  float new_energy = amber.getEnergy();

  cout << "energy before/after minpulation: " << initial_energy << "/" << new_energy << " kJ/mol" << endl;

  // write the new 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 forcefield evaluations can be restricted by a selection, e.g. hydrogen atoms by adding the following lines after the instantiation of the forcefield:

S.deselect();
// add all atoms to the field
amber.setup(S);
// select all hydrogens
Selector selector("element(H)");
S.apply(selector);
// compute energy of selection in the field of all atoms in S
float new_energy = amber.getEnergy();