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