Coordinate systems and the CoordSys class¶
In NITROGEN, a coordinate system (CS) refers to any mapping between a set of curvilinear
coordinates and a set of Cartesian-like or rectangular coordinates (i.e.,
those with a constant, diagonal metric tensor). The most common situation for
molecular nuclear motion problems is a set of internal coordinates (like bond
lengths and angles) that defines the values of 3N atomic Cartesian coordinates,
but other CSs that do not explicitly refer to 3-D space are equally
valid. A CS is implemented with the CoordSys
class and its sub-classes, all of which are themselves
sub-classes of DFun.
Getting started¶
Let’s explore CoordSys objects with a simple
built-in CS for triatomic systems, Valence3,
a sub-class of CoordSys:
>>> import nitrogen as n2
>>> cs = n2.coordsys.Valence3()
>>> isinstance(cs, n2.coordsys.CoordSys)
True
CoordSys objects have some simple descriptor
attributes such as a name, the number of input (curvilinear) coordinates Q
and output (rectangular) coordinates X, and coordinate labels.
>>> cs.name
'Triatomic valence'
>>> cs.nQ, cs.nX
(3, 9)
>>> cs.Qstr
['r1', 'r2', 'theta']
>>> cs.Xstr
['x0', 'y0', 'z0', 'x1', 'y1', 'z1', 'x2', 'y2', 'z2']
Valence3 has three input coordinates: two bond
lengths (r1 and r2) and a bond angle (theta, in radians). The
output coordinates are the Cartesian positions of the three atoms placed at
\((0,0,-r_1)\), the origin, and
\((0, r_2 \sin\theta, -r_2 \cos\theta )\), respectively.
Their values
and derivatives are calculated with the Q2X()
class method.
>>> cs.Q2X(np.array([1.0, 1.3, 2*n2.pi/3])) # deriv = 0 by default
array([[ 0. , 0. , -1. , 0. , 0. ,
0. , 0. , 1.12583302, 0.65 ]])
>>> cs.Q2X(np.array([1.0, 1.3, 2*n2.pi/3]), deriv = 1)
array([[ 0. , 0. , -1. , 0. , 0. ,
0. , 0. , 1.12583302, 0.65 ],
[ 0. , 0. , -1. , 0. , 0. ,
0. , 0. , 0. , 0. ],
[ 0. , 0. , -0. , 0. , 0. ,
0. , 0. , 0.8660254 , 0.5 ],
[ 0. , 0. , -0. , 0. , 0. ,
0. , 0. , -0.65 , 1.12583302]])
CoordSys objects are also instances of
DFun with the expected attributes.
>>> cs.nQ == cs.nx # DFun number of inputs
True
>>> cs.nX == cs.nf # DFun number of outputs
True
Q2X() is just a wrapper for
DFun.f() and returns the same standard
derivative array. (See this section on derivative arrays.)
Z-matrix coordinates¶
The standard Z-matrix provides a useful general purpose
internal coordinate system. It is implemented with the
built-in CoordSys sub-class ZMAT. Building
a Z-matrix CS requires passing a Z-matrix definition string to the
ZMAT initializer. This example creates a Z-matrix
CS that is functionally equivalent to Valence3:
>>> zmat = """
... H
... O 1 r1
... H 2 r2 1 theta
... """
>>> cs = n2.coordsys.ZMAT(zmat)
(Note that the default unit for ZMAT
angles is degrees, not radians.)
The definition string uses standard Z-matrix conventions with a few
additional features. Any Z-matrix value with a valid variable label
(i.e. beginning with a letter, not a number) becomes a coordinate.
The same coordinate label can be used twice, which will constrain
multiple Z-matrix values to be equal:
>>> zmat = """
... H
... O 1 r
... H 2 r 1 theta
... """
The order of coordinates in the CS is the order of their first appearance in the Z-matrix definition string. A particular Z-matrix value can be held fixed and excluded from the CS by giving it a literal numeric value.
>>> zmat = """
... H
... O 1 r1
... H 2 r2 1 105.0
... """
Coordinate labels can also be prefixed with a sign (+ or -) or
a constant numerical coefficient:
>>> zmat = """
... H
... O 1 +r
... H 2 -r 1 2.3*a
... """
The atom label at the start of each line is required, but has
no meaning other than to indicate a dummy atom, for which the label is X
or x.
>>> zmat = """
... H
... O 1 r1
... H 2 r2 1 a1
... X 3 r3 2 a2 1 tau
... """
The above Z-matrix will have six input coordinates, but generate only nine Cartesian output coordinates for the first three atoms. The Cartesian coordinates for dummy atoms are not included as output coordinates.
Coordinate transformations¶
Because CoordSys objects are also instances of
DFun, general non-linear coordinate transformations
can be implemented using DFun composition. (See this
section of the DFun tutorial.) A CoordSys
object XQ that implements \(X(Q)\) and a DFun object
QQp that implements \(Q(Q')\), where \(Q'\) are the new, transformed coordinates
can be combined to generate the transformed coordinate system \(X(Q(Q'))\) via
>>> XQp = QQp ** XQ
>>> XQp = XQ @ QQp # functionally equivalent
In this case, because QQp is only an instance of DFun,
the most generic parent class, the resulting composition is also an instance of
DFun, not CoordSys.
To formalize the special role of coordinate transformations, NITROGEN provides
yet another DFun sub-class, called
CoordTrans, for implementing coordinate transformation
functions. A simple and commonly used transformation is a linear one,
\(Q_i = T_{ij} Q'_j\),
where \(T_{ij}\) is a constant matrix. This is implemented by the built-in
LinearTrans, a sub-class of
CoordTrans. As an example, let’s use a linear
transformation to construct symmetrized coordinates for water. We begin again
with the Z-matrix CS:
>>> zmat = """
... H
... O 1 r1
... H 2 r2 1 theta
... """
>>> cs = n2.coordsys.ZMAT(zmat)
We define \(r_s = (r_1 + r_2)/2\) and \(r_a = (r_1 - r_2)/2\), so the transformation matrix is \(T = ((1, 1, 0), (1, -1, 0), (0, 0, 1))\)
>>> T = np.array([[1.0, 1.0, 0], [1.0, -1.0, 0], [0, 0, 1.0]])
>>> ct = n2.coordsys.LinearTrans(T, Qpstr = ['rs', 'ra', 'theta'], name = 'sym')
>>> cs_sym = ct ** cs
>>> isinstance(cs_sym, n2.coordsys.CoordSys)
True
Coordinate system diagrams¶
Coordinate system and transformation objects provide a
diagram() method that generates a
string-based pictorial representation of the CS. This can help visualize
a chain of possibly many transformations to keep track of its sequence and order.
Note that the diagram strings use unicode characters, and the appearance may
be affected by your font settings.
>>> print(cs_sym.diagram())
│↓ ↑│
│Q'[3] │
╔═╧══════╗ │
║ ║ │
║ sym ║ │
║ ║ │
╚═╤══════╝ │
│Q [3] │
│↓ ↑│
│Q [3] [9] X│
╔═╪════════════════╪═╗
║ │ ┌────────────┐ │ ║
║ ╰─┤ Z-matrix ├─╯ ║
║ └────────────┘ ║
╚════════════════════╝
This diagram shows that the inputs are now the new \(Q'\) coordinates, which are transformed to the original Z-matrix \(Q\) coordinates, which are finally used to compute the output Cartesian coordinates \(X\). The numbers in brackets indicate the number of coordinates in each segment.
Atomic coordinate systems¶
Atomic coordinate systems refer specifically to CoordSys
objects whose output coordinates are the \(3N\) Cartesian positions of
\(N\) particles in 3-D space. The coordinates must be ordered
\(x_0, y_0, z_0, x_1, y_1, z_1, \ldots\). The object attribute
isatomic is True for atomic coordinate systems.
Certain class methods, such as Q2g(), have
different behavior or options for atomic vs. non-atomic coordinate systems.
CoordSys objects for molecular problems,
like the built-in class ZMAT, should
generally be atomic. The isatomic value of
transformed coordinate systems is inherited from that of the
untransformed coordinate sytem:
>>> cs.isatomic # instance of ZMAT
True
>>> cs_sym.isatomic # symmetry-transformed ZMAT coordinates
True