# -*- coding: utf-8 -*-
"""Contains a class for using the 3D Balaban index 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_balaban(cv):
"""
Balaban index as collective variable, 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)
update_adjacency : bool
update the adjacency in every dynamics step (default: False)
"""
def __init__(self, idx, symbols, **kwargs):
"""Construct Balaban index 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)
self.update_am = self.get_kwargs("update_adjacency", default=False)
self.am = np.zeros((len(self.symbols), len(self.symbols)))
self.au2unit = (c.a0*c.a0)**(-0.5)
self.reset_width()
[docs] def get_adjacencymatrix(self, coords):
"""
Calculate the adjacency matrix for the given coordinates.
Parameters
----------
coords : np.2darray (shape: (N, 3))
atomic xyz-coordinates of the molecule (bohr)
Returns
-------
am : np.2darray (shape: (N, N))
adjacency matrix
"""
nat = len(coords)
rcov = {"C": 0.77/c.a0, "H": 0.37/c.a0, "O": 0.66/c.a0}
threshold = 0.2
am = 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
rcov_sum = rcov[self.symbols[p[0]]] + rcov[self.symbols[p[1]]]
if np.abs(cvf.bond(coords[[p]]) - rcov_sum) < threshold * rcov_sum:
am[p] = 1.0 # cvf.bond(coords[[p]])
am += am.T
return am[~np.all(am==0, axis=1)][:, ~np.all(am==0, axis=0)]
[docs] def get_edges(self, am):
"""
Return edges from the adjacency matrix.
Parameters
----------
am : np.2darray (shape: (N, N))
adjacency matrix
Returns
-------
edges : np.2darray (shape: (Nedges, 2))
pairs of connected atoms
"""
return np.transpose(np.nonzero(np.triu(am)))
[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_ddistancematrix(self, coords):
"""
Return the gradient of the geometrical distance matrix.
Parameters
----------
coords : np.2darray (shape: (N, 3))
atomic xyz-coordinates of the molecule (bohr)
Returns
-------
ddm : np.2darray (shape: (N, N, N, 3))
gradient of the geometrical distance matrix
"""
nat = len(coords)
ddm = np.zeros((nat, nat, nat, 3))
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
ddm[p[0], p[1]][[p]] = cvf.dbond(coords[[p]])
ddm[p[1], p[0]][[p]] = ddm[p[0], p[1]][[p]]
return ddm[~np.all(ddm==0, axis=(1, 2, 3))][:, ~np.all(ddm==0, axis=(0, 2, 3))]
[docs] def get_balaban(self, dm, edges):
"""
Return Balaban index.
Parameters
----------
dm : np.2darray (shape: (N, N))
geometrical distance matrix
edges : np.2darray (shape: (Nedges, 2))
pairs of connected atoms
Returns
-------
balaban : float
3D Balaban index
"""
distsum = np.sum(dm, axis=0)
n = np.size(distsum)
m = np.size(edges, axis=0)
balaban = 0.0
for p in edges:
balaban += np.prod(distsum[[p]])**(-0.5)
return (m/(m-n+2))*balaban
[docs] def get_dbalaban(self, dm, ddm, edges):
"""
Return the gradient of the Balaban index.
Parameters
----------
dm : np.2darray (shape: (N, N))
geometrical distance matrix
ddm : np.2darray (shape: (N, N, N, 3))
gradient of the geometrical distance matrix
edges : np.2darray (shape: (Nedges, 2))
pairs of connected atoms
Returns
-------
dbalaban : np.2darray (shape: (N, 3))
gradient of the 3D Balaban index
"""
distsum = np.sum(dm, axis=0)
ddistsum = np.sum(ddm, axis=0)
n = np.size(dm, axis=0)
m = np.size(edges, axis=0)
dbalaban = np.zeros_like(ddm[0, 0])
for p in edges:
dbalaban += -0.5*np.prod(distsum[[p]])**(-1.5)* \
(distsum[p[0]]*ddistsum[p[1]]+distsum[p[1]]*ddistsum[p[0]])
return (m/(m-n+2))*dbalaban
[docs] def set_s_and_ds(self, coords):
"""
Calculate 3D Balaban index and its gradient for the current set of coordinates.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
"""
if not self.am.any() or self.update_am:
self.am = self.get_adjacencymatrix(coords)
edges = self.get_edges(self.am)
dm = self.get_distancematrix(coords)
ddm = self.get_ddistancematrix(coords)
balaban = self.get_balaban(dm, edges)
dbalaban = self.get_dbalaban(dm, ddm, edges)
self.s = balaban
self.ds_dr = dbalaban
print("balaban: ", balaban*self.au2unit)
print("dbalaban:\n", dbalaban)