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.

This tutorial is in progress, there may be missing example files, source files or bugs in the code.

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 [[]]

 	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 me 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.

	// 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