# -*- coding: utf-8 -*-
"""Contains a class for using the Mulliken charge as a CV."""
from __future__ import print_function
import numpy as np
import os
import multiprocessing as mp
from . import mulliken as mul
from .cv import cv
[docs]class cv_mullcharge(cv):
"""
Mulliken charge 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
-----------------
atoms : list of int
list of three atom indices for the calculation of the angle
iface : str
interface to QC code
method : str
'analytic' or 'numeric' derivatives
nproc : int
number of processors to be used for numeric differentiation
"""
def __init__(self, idx, symbols, **kwargs):
"""Construct Mulliken charge CV object."""
cv.__init__(self, idx, symbols, **kwargs)
if type(kwargs["atoms"]) == type([]):
self.atoms = kwargs["atoms"]
else:
self.atoms = [kwargs["atoms"]]
self.iface = self.get_kwargs("iface", default="turbomole")
self.x = np.linspace(-1.0, 1.0, 200)
self.ticks = None
if self.iface == "dftbaby":
self.method = self.get_kwargs("method")
self.tmp = "tmp_charges.npy"
if self.iface == "turbomole":
self.check_ridft()
self.nproc = self.get_kwargs("nproc", default=1)
[docs] def check_ridft(self):
"""Check if a ridft-calculation is performed."""
self.ridft = False
with open("tm/control") as f:
for line in f:
if "$rij" in line:
self.ridft = True
break
[docs] def set_dftb_object(self, dftb):
"""Set the `DFTB.PES.PotentialEnergySurfaces` object."""
self.dftb = dftb
[docs] def set_s_and_ds(self, coords):
"""
Calculate Mulliken charge and its gradient for the current set of coordinates.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
"""
nat = len(self.symbols)
if self.iface == "dftbaby":
charge = np.sum(self.dftb.tddftb.dftb2.dq[self.atoms])
if self.method == "analytic":
dcharges = self.dftb.grads.getChargeGradients()
dcharge = np.sum(dcharges[:, self.atoms], axis=1).reshape(nat, 3)
print("charge: ", charge)
print("dcharge:\n", dcharge)
return charge, dcharge
elif self.iface == "turbomole":
charge = np.sum(mul.get_charges(self.symbols)[self.atoms])
# numerical charge gradients:
h = 1e-3
dcharge = np.zeros(coords.shape)
if os.path.exists(self.tmp):
charges_h = np.load(self.tmp)
for i in range(3*nat):
dcharge[i/3, i % 3] = (np.sum(charges_h[i/3, i % 3, self.atoms])-charge)/h
else:
pool = mp.Pool(processes=self.nproc)
coords_h = []
for i in range(3*nat):
coords_h.append(np.copy(coords))
coords_h[i][i/3, i % 3] += h
i_sym_coords_h = [(i, self.symbols, coord) for i, coord in enumerate(coords_h)]
if self.iface == "dftbaby":
charges_h = np.array(pool.map(mul.get_new_charges_dftb, i_sym_coords_h))
elif self.iface == "turbomole":
charges_h = np.array(pool.map(mul.get_new_charges_tm, i_sym_coords_h))
pool.close()
pool.join()
dcharge = (np.sum(charges_h[:, self.atoms], axis=1).reshape(nat, 3)-charge)/h
np.save(self.tmp, charges_h)
print("charge: ", charge)
print("dcharge:\n", dcharge)
self.s = charge
self.ds_dr = dcharge