1#!/usr/local/bin/python3.8 2 3from __future__ import print_function 4 5import numpy 6 7from molmod import * 8from molmod.io import FCHKFile 9 10 11class InternalCoordinate(object): 12 """Abstract base class for all internal coordinates.""" 13 def __init__(self, indexes, icfn, conversion=1.0): 14 """ 15 Arguments: 16 | ``indexes`` -- The indexes of the atoms in the internal 17 coordinate. The order must be the same as the order 18 of the mandatory arguments of icfn. 19 | ``icfn`` -- a function from molmod.ic that can compute the 20 internal coordinate and its derivatives. 21 | ``conversion`` -- In case the internal coordinate does not have a 22 unit of length, then this conversion factor is 23 used to convert it to a length unit. This way, 24 the Jacobian becomes a dimensionless constant. 25 26 All the Jacobian-logic is implemented in this abstract class. 27 """ 28 self.indexes = indexes 29 self.icfn = icfn 30 self.conversion = conversion 31 32 def fill_jacobian_column(self, jaccol, coordinates): 33 """Fill in a column of the Jacobian. 34 35 Arguments: 36 | ``jaccol`` -- The column of Jacobian to which the result must be 37 added. 38 | ``coordinates`` -- A numpy array with Cartesian coordinates, 39 shape=(N,3) 40 """ 41 q, g = self.icfn(coordinates[list(self.indexes)], 1) 42 for i, j in enumerate(self.indexes): 43 jaccol[3*j:3*j+3] += g[i] 44 return jaccol 45 46 47class BondLength(InternalCoordinate): 48 def __init__(self, i0, i1): 49 InternalCoordinate.__init__(self, (i0, i1), bond_length) 50 51 52class BendingAngle(InternalCoordinate): 53 def __init__(self, i0, i1, i2): 54 InternalCoordinate.__init__(self, (i0, i1, i2), bend_angle, angstrom/(5*deg)) 55 56 57class DihedralAngle(InternalCoordinate): 58 def __init__(self, i0, i1, i2, i3): 59 InternalCoordinate.__init__(self, (i0, i1, i2, i3), dihed_angle, angstrom/(5*deg)) 60 61 62def setup_ics(graph): 63 """Make a list of internal coordinates based on the graph 64 65 Argument: 66 | ``graph`` -- A Graph instance. 67 68 The list of internal coordinates will include all bond lengths, all 69 bending angles, and all dihedral angles. 70 """ 71 ics = [] 72 # A) Collect all bonds. 73 for i0, i1 in graph.edges: 74 ics.append(BondLength(i0, i1)) 75 # B) Collect all bends. (see b_bending_angles.py for the explanation) 76 for i1 in range(graph.num_vertices): 77 n = list(graph.neighbors[i1]) 78 for index, i0 in enumerate(n): 79 for i2 in n[:index]: 80 ics.append(BendingAngle(i0, i1, i2)) 81 # C) Collect all dihedrals. 82 for i1, i2 in graph.edges: 83 for i0 in graph.neighbors[i1]: 84 if i0==i2: 85 # All four indexes must be different. 86 continue 87 for i3 in graph.neighbors[i2]: 88 if i3==i1 or i3==i0: 89 # All four indexes must be different. 90 continue 91 ics.append(DihedralAngle(i0, i1, i2, i3)) 92 return ics 93 94 95def compute_jacobian(ics, coordinates): 96 """Construct a Jacobian for the given internal and Cartesian coordinates 97 98 Arguments: 99 | ``ics`` -- A list of internal coordinate objects. 100 | ``coordinates`` -- A numpy array with Cartesian coordinates, 101 shape=(N,3) 102 103 The return value will be a numpy array with the Jacobian matrix. There 104 will be a column for each internal coordinate, and a row for each 105 Cartesian coordinate (3*N rows). 106 """ 107 N3 = coordinates.size 108 jacobian = numpy.zeros((N3, len(ics)), float) 109 for j, ic in enumerate(ics): 110 # Let the ic object fill in each column of the Jacobian. 111 ic.fill_jacobian_column(jacobian[:,j], coordinates) 112 return jacobian 113 114 115# This if block is only executed when this file is ran as a program, and not 116# when it is loaded as a module. 117if __name__ == "__main__": 118 # Load the formatted checkpoint file with the frequency computation. This 119 # file also contains the atomic numbers and the coordinates of the atoms, 120 # and therefore one can access the dopamine molecule object through 121 # fchk.molecule. 122 fchk = FCHKFile("dopamine.fchk") 123 # Set the default graph for the construction of the internal coordinates: 124 fchk.molecule.set_default_graph() 125 # Setup a list of internal coordinates 126 ics = setup_ics(fchk.molecule.graph) 127 # Compute the Jacobian. 128 J = compute_jacobian(ics, fchk.molecule.coordinates) 129 # Compute the pseudo-inverse, using a loose threshold for the singular 130 # values to filter out equivalent internal coordinates. 131 Jinv = numpy.linalg.pinv(J, 1e-5) 132 # Get the Hessian in Cartesian coordinates. 133 H = fchk.get_hessian() 134 # Transform to internal coordinates. 135 K = numpy.dot(Jinv, numpy.dot(H, Jinv.transpose())) 136 # Make a nice printout of K. 137 print("The Hessian in internal coordinates in kcal/mol/angstrom**2") 138 unit = kcalmol/angstrom**2 139 for row in K: 140 print(" ".join("% 5.0f" % (v/unit) for v in row)) 141