Aligning Structures and Redundant Conformer Elimination¶
import cctk
is assumed.
Calculating RMSD¶
We can calculate the root-mean-square deviation between two conformers.
Calculating the RMSD between two non-conformers is not supported (yet).
No alignment is performed; the two geometries are compared as-is.
# load molecule
path = "test/static/gaussian_file.out"
gaussian_file = cctk.GaussianFile.read_file(path)
ensemble = gaussian_file.ensemble
m1 = ensemble.molecules[0]
m2 = ensemble.molecules[-1]
# measure RMSD in Angstroms
RMSD = cctk.helper_functions.compute_RMSD(m1,m2)
RMSD == 0.0006419131435567976
Aligning Two Molecules¶
We can align all the conformers in an ensemble onto one of it members.
All conformers are translated to the origin. Then, a different rotation is applied to each conformer until the RMSD between the conformer and the template conformer is minimized.
A new
ConformationalEnsemble
is created and the originalMolecule
objects are not modified.The Kabsch algorithm is used.
Alignment of non-conformers is not supported (yet).
# create a ConformationalEnsemble
path = "test/static/phenylpropane*.out"
conformational_ensemble = cctk.ConformationalEnsemble()
for filename in sorted(glob.glob(path)):
gaussian_file = cctk.GaussianFile.read_file(filename)
ensemble = gaussian_file.ensemble
molecule = ensemble.molecules[-1]
properties_dict = ensemble.get_properties_dict(molecule)
conformational_ensemble.add_molecule(molecule,properties_dict)
# define which atoms will be used for alignment
# alternatively, set comparison_atoms to "heavy" or "all"
comparison_atoms = [1,2,3,4,5,6]
# perform the alignment
# to_geometry is the template geometry (0-indexed)
# set compute_RMSD to True to calculate the RMSD before and after the aligment
aligned_ensemble, before_RMSD, after_RMSD = conformational_ensemble.align(to_geometry=np.int64(0), comparison_atoms=comparison_atoms, compute_RMSD=True)
# write out the result to a Gaussian input file
# open in GaussView and select "Single new molecule group for all files" in the "Open Files" dialog box under "Target"
# select geometries of interest of "select all" to compare the aligned molecules
cctk.GaussianFile.write_ensemble_to_file("test/static/phenylpropane_aligned.gjf", aligned_ensemble, "#p")
# alternatively, write all structures to a .mol2 file and open in Maestro or a similar program
# it be necessary to run molecule.assign_connectivity() on each member of the Ensemble first
# if the bond connectivity information is not available
cctk.MOL2File.write_ensemble_to_file("test/static/phenylpropane_aligned.mol2", aligned_ensemble)
Redundant Conformer Elimination¶
Useful for removing redundant conformers after conformational searches.
Every geometry is aligned and then a new
ConformationalEnsemble
is created that only contains unique conformations.Whenever redundancies are found, the lowest energy conformer is kept.
The current
ConformationalEnsemble
is not modified.
# need a ConformationalEnsemble
assert isinstance(conformational_ensemble, ConformationalEnsemble)
# eliminate redundant conformers
# RMSD_cutoff is in Angstroms
# comparison_atoms is a list of atom numbers, "heavy" (default), or "all"
unique_conformers = conformational_ensemble.eliminate_redundant(RMSD_cutoff=0.5, comparison_atoms="heavy")
# how many conformers are left
len(unique_conformers)
# write out the entire ensemble
# alternatively, index the ensemble to retrieve a subset of the structures first
cctk.GaussianFile.write_ensemble_to_file("test/static/phenylpropane_aligned2.gjf", ensemble2, "#p")