# -*- coding: utf-8 -*-
"""Contains a class for using the lowest eigenvalue of the geometric distance matrix as a CV."""
from __future__ import print_function
import numpy as np
import itertools as it
from .. import constants as c
from . import cvfunctions as cvf
from .cv import cv
[docs]class cv_firsteigenvalue(cv):
"""
Lowest eigenvalue of the distance matrix as CV, inherits from `metafalcon.cvs.cv.cv`.
Parameters
----------
idx : int
index of collective variable in the input file
symbols : list of str
atomic symbols of the molecule
**kwargs :
See below
Keyword Arguments
-----------------
exclude_h : bool
do not include distances to hydrogen atoms (default: False)
"""
def __init__(self, idx, symbols, **kwargs):
"""Construct first eigenvalue CV object."""
cv.__init__(self, idx, symbols, **kwargs)
self.x = np.linspace(0., 1., 200)
self.exclude_h = self.get_kwargs("exclude_h", default=False)
if self.exclude_h:
self.no_h = [i for i, s in enumerate(self.symbols) if s is not "H"]
self.au2unit = c.a0
self.reset_width()
[docs] def get_distancematrix(self, coords):
"""
Return the geometrical distance matrix.
Parameters
----------
coords : np.2darray (shape: (N, 3))
atomic xyz-coordinates of the molecule (bohr)
Returns
-------
dm : np.2darray (shape: (N, N))
geometrical distance matrix
"""
nat = len(coords)
dm = np.zeros((nat, nat))
for p in it.combinations(range(nat), 2):
if self.exclude_h == True and ("H" in [self.symbols[p_i] for p_i in p]):
continue
dm[p] = cvf.bond(coords[[p]])
dm += dm.T
return dm[~np.all(dm==0, axis=1)][:, ~np.all(dm==0, axis=0)]
[docs] def get_derivative_eigenvalue(self, e, v, k):
"""
Return the derivative of the k-th eigenvalue with respect to distance matrix.
Parameters
----------
e : np.array
eigenvalues (shape: (N))
v : np.2darray
eigenvectors (shape: (N, N))
k : int
index of the eigenvalue to by differentiated
Returns
-------
de_dd : np.2darray
derivative of eigenvalue k wrt distance matrix (shape: (N, N))
"""
derivative = np.zeros_like(v)
for p in it.combinations(range(len(e)), 2):
derivative[p] = 2*v[p[0], k]*v[p[1], k] # factor 2 because of symmetric matrix
derivative += derivative.T
return derivative # only upper triangular matrix, because diagonal elements not needed
[docs] def set_s_and_ds(self, coords):
"""
Calculate first eigenvalue and its gradient for the current set of coordinates.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
"""
dm = self.get_distancematrix(coords) # shape: (nat, nat) or (nat_noh, nat_noh)
e, v = np.linalg.eigh(dm) # shape e: (nat) or (nat_noh)
de_dd = self.get_derivative_eigenvalue(e, v, -1) # shape: (nat, nat) or (nat_noh, nat_noh)
eigen = e[-1]
deigen = np.zeros_like(coords)
for p in it.combinations(range(len(e)), 2):
deigen[[p]] += de_dd[p]*cvf.dbond(coords[[p]])
if self.exclude_h:
deigen_tmp = deigen
deigen = np.zeros_like(coords)
deigen[self.no_h] = deigen_tmp[:len(self.no_h)]
self.s = eigen
self.ds_dr = deigen
print("eigen: ", eigen*self.au2unit)
print("deigen:\n", deigen)