tmpchem/computational_chemistry

Files used in TMP Chem videos on computational chemistry

TMP Chem Computational Chemistry Source Code

This repository contains scripts, programs, and data files used with the "Computational Chemistry" playlist on the TMP Chem YouTube channel. Primary language is Python3, with toy programs to demonstrate modeling and analysis methods. Repository is always subject to change, and no guarantees are made of the correctness of output.

Table of contents

• Getting started
• Prerequisites
• Installation
• Running the tests
• Running the scripts
  • Geometry analysis
  • Molecular mechanics
• Author
• Acknowledgments

Getting started

Most recent version of project is located in TMP Chem GitHub account. Download by following instructions for Git clone from GitHub. Requires a terminal, IDE, etc. to execute Python scripts and ability to write to file system within project directories.

Prerequisites

Requires Python 3.5 or greater for script execution. Requires access to numpy and matplotlib modules. All prerequisites can be met by downloading and using Python from most recent Anaconda distribution.

Installation

No additional installation necessary after cloning the repository.

Running the tests

WARNING: Test suite is incomplete and subject to change without notice.

To run tests for a project, go to the script directory for that project, [top_level_path]/scripts/[project]. If present, execute run_tests.py for the project.

python run_tests.py

This command executes a test suite of unit tests from each test module present in the mmlib directory. Each test module contains a set of unit tests of methods within the module, confirming proper behavior and protecting against regression errors and system misconfigurations.

Once executed, standard output will indicate success with an 'OK' message and the number of executed unit tests as well as total run time. Any other message indicates failure and will include the nature of the failing tests.

Open an issue to give feedback or report bugs.

Running the scripts

As of 16 Feb 2017, repository contains two projects: geometry_analysis, and molecular_mechanics.

Geometry analysis

The geometry_analysis project contains scripts which take an xyz-format molecular geometry file as input, and output to screen associated geometry data, including bond lengths, bond angles, torsion angles, out-of-plane angles, center of mass, and/or moment of inertia, etc. Sample xyz files are located in [top_level_path]/geom/xyz directory.

Molecular mechanics

The molecular_mechanics project contains scripts to compute molecular mechanics energy of a system (mm.py), molecular dynamics trajectories (md.py), Metropolis Monte Carlo ensembles (mc.py), and optimize molecular coordinates to potential energy minima (opt.py).

The energy function and parameters in all cases is based on AMBER FF94.

Cornell et. al, 
J. Am. Chem. Soc. 1995, 117, 5179-5197.
doi.org/10.1021/ja00124a002

Energy function in Equation 1. Atom types in Table 1. Parameter values in Table 14. Download AmberTools15 from here. After unzipping, parameters located in amber14/dat/leap/parm/parm94.dat.

Sample input files for mm are located in [top_level_path]/geom/[file_type] directories, where file_type] is xyzq or prm. Sample input files for md and mc are located in [top_level_path]/geom/sim directory. Output files for each are demonstrated in samples, and may be written to any accessible file name.

Author

The sole author of this package is Trent M. Parker.

• Email: [email protected]
• LinkedIn: https://www.linkedin.com/in/tmpchem
• Youtube: https://www.youtube.com/tmpchem

Acknowledgments

The author wishes to thank the following individuals: Dr. Michael S. Marshall, for encouraging him to learn the Python language; Dr. Lori S. Burns, for encouraging him to conform to style guidelines for readable Python code; Dr. Michael A. Lewis, for providing the inspiration for the author to initiate studies in this field; and Dr. C. David Sherrill, for providing an enormous number of opportunities to continue to learn and grow as a scientist and a person.

computational_chemistry's Issues

Z matrix to Cartesian converter uses average atomic weights

Hello
Thanks for the great youtube courses and utilities. I watched a many :)

I think the weights should be that of the most abundant isotope or relevant isotope in question, not the average weights. This is how other packages do it for ab-initio.

I use this reference:
https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=&ascii=ascii

I noticed while testing my C++ Z matrix code v psi4 which agree, but yours code gives very slightly different centre of mass coordinates when I add the following code (to your utility) to get the centre of mass adjusted Cartesian coordinates. Ultimately, this is caused by the weights you use.

Sorry. I don't know python to write code in normally, but this works okay.

print xyz coordinates to screen

def print_coords(self, comment):
    angstrom_to_bohr = 1.889725989
    total_mass = 0.0
    
    for i in range(self.n_atoms):
        total_mass += at_masses[self.atoms[i].attype]

    x_com = 0.0
    y_com = 0.0
    z_com = 0.0

    for i in range(self.n_atoms):
        x_com += self.atoms[i].coords[0] * at_masses[self.atoms[i].attype]
        y_com += self.atoms[i].coords[1] * at_masses[self.atoms[i].attype]
        z_com += self.atoms[i].coords[2] * at_masses[self.atoms[i].attype]
    
    x_com /= total_mass
    y_com /= total_mass
    z_com /= total_mass

    for i in range(self.n_atoms):
        self.atoms[i].coords[0] -= x_com
        self.atoms[i].coords[1] -= y_com
        self.atoms[i].coords[2] -= z_com
        
    print('\n')
    print('%s\n' % (comment), end='')
    for i in range(self.n_atoms):
        if self.atoms[i].attype == 'X':
                continue
        
        print('%-2s' % (self.atoms[i].attype), end='')
        
        for j in range(3):
            print(' %14.10f' % (self.atoms[i].coords[j] * angstrom_to_bohr), end='')
        
        print('\n', end='')

I think a bit higher precision in the printed output and definitions of weights would also be useful due to getting rounding errors you get from the cartesians you produce when used in calculations when comparing with others.

o' course weight are not used anyway in the utility to convert, so it is minor. The coordinates without COM adjustment are correct up to that point to within rounding errors.

Best regards
...

Kinetic pressure calculation

I was looking at the calculations for ensemble properties of the system, and I found something confusing in regarding the calculation for kinetic pressure, defined as follows:

def get_virial(mol):
    """Clausius virial function for all atoms, force,s and coordinates.
    
    Args:
        mol (mmlib.molecule.Molecule): Molecule object with coordinate
            and force data.
    """
    mol.virial = 0.0
    for i in range(mol.n_atoms):
        for j in range(3):
            mol.virial += -mol.atoms[i].coords[j] * mol.g_total[i][j]

def get_pressure(mol):
    """Update total pressure of a system of molecules with boundary.
    
    Args:
        mol (mmlib.molecule.Molecule): Molecule object with temperature,
            volume, and virial data.
    """
    get_virial(mol)
    pv = mol.n_atoms * energy.kb() * mol.temp
    pv += mol.virial / (3.0 * mol.n_atoms)
    mol.press = kcalamol2pa() * pv / mol.vol

In equation form:

After spending an afternoon going through as much MD literature I could find, I've yet to see any paper that lists that N dividing the virial when calculating the pressure. Some examples:

Page 2: https://arxiv.org/pdf/1111.2705.pdf
Page 1: https://spiral.imperial.ac.uk/bitstream/10044/1/262/1/edm2006jcp.pdf
Page 5: http://phys.ubbcluj.ro/~tbeu/MD/C2_for.pdf

Not only that, but the literature is consistent as expected, and apparently in disagreement with this implementation. I am probably missing something, since I'm no expert in MD and the notation could be the culprit here. In any case, I'd like to know what's the problem.

Sorting convention for principal axis

Hi tmpchem

I'm wondering what is the sorting convention you are using to align molecules with respect to the principal axis (ascending or descending)?. I'm trying to compute the gyration tensor of a macromolecule knowing beforehand its relationship to the inertia tensor (see appendix of dx.doi.org/10.1021/jp2065612).
In my case, all atoms have the same mass, and I would expect to get the same coordinates for both, the gyration tensor aligned molecule and inertia tensor aligned molecule.

computational_chemistry/scripts/geometry_analysis/geometry_analysis.py

Lines 368 to 387 in 329559f

	# rotate molecule to inertial frame of principal moments
	def get_inertial_coords(geom, moi):
	moi_eigvals, moi_eigvecs = numpy.linalg.eig(moi)
	coords = geom[1]
	coords = numpy.array(numpy.dot(coords, moi_eigvecs))
	geom[1] = coords
	moi = get_moi(geom)
	order = [0, 1, 2]
	for p in range(3):
	for q in range(p+1, 3):
	if (moi.item(p, p) < moi.item(q, q)):
	temp = order[p]
	order[p] = order[q]
	order[q] = temp
	moveaxes = numpy.zeros((3, 3))
	for p in range(3):
	moveaxes[p][order[p]] = 1.0
	coords = numpy.dot(coords, moveaxes)
	geom[1] = coords
	return geom

Many thanks in advance for your time

encountered error

when i run below line in jupyter notebook
run ./scripts/molecular_mechanics/mm.py ./geom/prm/co.prm

i get the following errors
ValueError Traceback (most recent call last)
~\computational_chemistry-master\scripts\molecular_mechanics\mm.py in
24
25 # Read in molecular geometry and topology.
---> 26 molecule = molecule.Molecule(input_file_name)
27
28 # Calculate potential energy.

~\computational_chemistry-master\scripts\molecular_mechanics\mmlib\molecule.py in init(self, infile_name)
480 self.GetTopology()
481 elif (self.filetype == 'prm'):
--> 482 self.ReadInPrm()
483 self.UpdateInternals()
484

~\computational_chemistry-master\scripts\molecular_mechanics\mmlib\molecule.py in ReadInPrm(self)
504 input_rows = fileio.GetFileStringArray(self.infile)
505
--> 506 self.atoms = fileio.GetAtomsFromPrm(input_rows)
507 self.bonds = fileio.GetBondsFromPrm(input_rows)
508 self.angles = fileio.GetAnglesFromPrm(input_rows)

~\computational_chemistry-master\scripts\molecular_mechanics\mmlib\fileio.py in GetAtomsFromPrm(rows)
401 for row in rows:
402 if row and row[0].upper() == 'ATOM':
--> 403 atoms.append(GetAtomFromPrm(row))
404 return atoms
405

~\computational_chemistry-master\scripts\molecular_mechanics\mmlib\fileio.py in GetAtomFromPrm(row)
231 if not _IsType(float, ro) or not float(ro) > 0.0:
232 raise ValueError(
--> 233 'Atomic vdw radius must be positive numeric value: %s' % ro)
234 if not _IsType(float, eps) or not float(eps) >= 0.0:
235 raise ValueError(

ValueError: Atomic vdw radius must be positive numeric value: 0.0

R_ij divisor terms

computational_chemistry/scripts/molecular_mechanics/mmlib/gradient.py

Line 162 in 608b092

gdir1 = geomcalc.GetUcp(u_21, cp) / r_21

Why is the radius showing up in the fraction when the vectors are supposed to be normalized already? It's skewing up the direction of max increasing angle for the middle atom in the direction of the shortest side. What's the interpretation behind it?

PS: This also appears to be present in torsion calculations and oop angle gradients.