Tutorial:Molecular alignments

This is an example on how to align whole structures or parts of structures in MSL. This tutorial explains functions in the Transform object.

Complete source of example_molecular_alignment.cpp

In MSL the program alignMolecules utilizes the type of code in this tutorial to structurally align multiple PDB files, given a set of reference atoms.

To compile

% make bin/example_molecular_alignment

To run the program

Go to the main directory and run the command (note, the location of the exampleFiles subdirectory needs to be provided as an argument)

% bin/example_molecular_alignment exampleFiles

Program description

First you just read two molecules into AtomContainers as seen on the AtomContainer Tutorial

 	string refFile = "example0005.pdb";
	refFile = (string)argv[1] + "/" + refFile;
	cout << "Create an AtomContainer and read the atoms from " << refFile << endl;


	string moveFile = "example0006.pdb";
	moveFile = (string)argv[1] + "/" + moveFile;
	cout << "Create an AtomContainer and read the atoms from " << moveFile << endl;

 

	// read the PDB into a new AtomContainer "refAtoms"
	AtomContainer refAtoms;
	if (!refAtoms.readPdb(refFile)) {
		// error checking, the PDB could not be read
		cerr << endl;
		cerr << "File " << refFile << " cannot be found, please speficy the path of the \"exampleFiles\" directory as an argument" << endl;
		exit(1);
	} else {
	  cout << "Reference PDB has "<<refAtoms.size()<<" atoms."<<endl;
	}

	// read the PDB into a new AtomContainer "moveAtoms"
	AtomContainer moveAtoms;
	if (!moveAtoms.readPdb(moveFile)) {
		// error checking, the PDB could not be read
		cerr << endl;
		cerr << "File " << moveFile << " cannot be found, please speficy the path of the \"exampleFiles\" directory as an argument" << endl;
		exit(1);
	} else {
	  cout << "Move PDB has "<<moveAtoms.size()<<" atoms."<<endl;
	}

Next we check before we move any molecules what the RMSD is, by using the AtomPointerVector::rmsd(AtomPointerVector) function.

	// Get RMSD before alignment
	double preRMSD = refAtoms.getAtomPointers().rmsd(moveAtoms.getAtomPointers());

Create a Transforms object and align the two AtomPointerVectors. Note AtomPointerVectors must contain the same number of atoms and it will try to align the atoms in order.

	// Create a Transforms object (includes functions to do molecular alignments)
	Transforms alignAtoms;

	// Do molecular alignment first atoms are moving, second atoms are the reference
	if (!alignAtoms.rmsdAlignment(moveAtoms.getAtomPointers(),refAtoms.getAtomPointers())){
	  cerr << "ERROR alignment of molecules using all the atoms of reference: "<<refFile<<" and move: "<<moveFile<<" files."<<endl;
	  exit(0);
	}

Finally, get the post-alignment RMSD and print out the results... including the Transformation matrix and translation

	// Get RMSD after alignment
	double postRMSD = refAtoms.getAtomPointers().rmsd(moveAtoms.getAtomPointers());

	// Query Transform object for the rotation and translations required for alignment
	Matrix rotMatrix           = alignAtoms.getLastRotationMatrix();
	CartesianPoint translation = alignAtoms.getLastTranslation();

	// Output RMSD, Rotation Matrix and Translation Vector
	cout << endl << endl;
	fprintf(stdout,"Before alignment RMSD: %8.3f\n",preRMSD);
	fprintf(stdout,"After  alignment RMSD: %8.3f\n",postRMSD);
	cout << endl << endl;

	cout << "Rotation Matix: " <<endl<<rotMatrix << endl;
	cout << "Translation Vector: "<<endl<<translation << endl;

---

Back to the tutorial page