import pandas as pd
import numpy as np
import chemcoord as cc
import sympy
sympy.init_printing()
molecule = cc.Cartesian.read_xyz('MIL53_beta.xyz', start_index=1)
r, theta = sympy.symbols('r, theta', real=True)
Let's build the construction table in order to bend one of the terephtalic acid ligands.
fragment = molecule.get_fragment([(12, 17), (55, 60)])
connection = np.array([[3, 99, 1, 12], [17, 3, 99, 12], [60, 3, 17, 12]])
connection = pd.DataFrame(connection[:, 1:], index=connection[:, 0], columns=['b', 'a', 'd'])
c_table = molecule.get_construction_table([(fragment, connection)])
molecule = molecule.loc[c_table.index]
zmolecule = molecule.get_zmat(c_table)
This gives the following movement:
zmolecule_symb = zmolecule.copy()
zmolecule_symb.safe_loc[3, 'angle'] += theta
cc.xyz_functions.view([zmolecule_symb.subs(theta, a).get_cartesian() for a in [-30, 0, 30]])
For the gradients it is very illustrating to compare: $$ f(x + h) \approx f(x) + f'(x) h $$
$f(x + h)$ will be zmolecule2
and $h$ will be dist_zmol
The boolean chain
argument denotes if the movement should be chained or not.
dist_zmol1 = zmolecule.copy()
r = 3
dist_zmol1.unsafe_loc[:, ['bond', 'angle', 'dihedral']] = 0
dist_zmol1.unsafe_loc[3, 'bond'] = r
cc.xyz_functions.view([molecule,
molecule + zmolecule.get_grad_cartesian(chain=False)(dist_zmol1),
molecule + zmolecule.get_grad_cartesian()(dist_zmol1),
(zmolecule + dist_zmol1).get_cartesian()])
angle = 30
dist_zmol2 = zmolecule.copy()
dist_zmol2.unsafe_loc[:, ['bond', 'angle', 'dihedral']] = 0
dist_zmol2.unsafe_loc[3, 'angle'] = angle
cc.xyz_functions.view([molecule,
molecule + zmolecule.get_grad_cartesian(chain=False)(dist_zmol2),
molecule + zmolecule.get_grad_cartesian()(dist_zmol2),
(zmolecule + dist_zmol2).get_cartesian()])
Note that the deviation between $f(x + h)$ and $f(x) + h f'(x)$ is not an error in the implementation but a visualisation of the small angle approximation.
The smaller the angle the better is the linearisation.
x_dist = 2
dist_mol = molecule.copy()
dist_mol.loc[:, ['x', 'y', 'z']] = 0.
dist_mol.loc[13, 'x'] = x_dist
zmat_dist = molecule.get_grad_zmat(c_table)(dist_mol)
It is immediately obvious, that only the ['bond', 'angle', 'dihedral']
of those atoms change,
which are either moved themselves in cartesian space or use moved references.
zmat_dist[(zmat_dist.loc[:, ['bond', 'angle', 'dihedral']] != 0).any(axis=1)]
atom | b | bond | a | angle | d | dihedral | |
---|---|---|---|---|---|---|---|
13 | O | 12 | 1.521692e+00 | 99 | -1.387881e+01 | 2 | 36.649299 |
11 | C | 13 | -8.852016e-01 | 12 | -6.848733e+01 | 99 | 112.390250 |
23 | C | 11 | 7.771561e-16 | 13 | -6.350987e+01 | 12 | -87.668150 |
21 | O | 11 | -3.330669e-16 | 13 | 6.350985e+01 | 12 | -140.238538 |
25 | C | 23 | 9.487976e-17 | 11 | -2.265493e-15 | 13 | -54.808771 |
27 | C | 23 | -5.332157e-17 | 11 | 1.255238e-15 | 13 | -54.808771 |
50 | H | 21 | 1.283378e-16 | 11 | -2.369105e-15 | 13 | 52.372802 |