Coverage for cclib/method/bader.py : 93%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
2#
3# Copyright (c) 2020, the cclib development team
4#
5# This file is part of cclib (http://cclib.github.io) and is distributed under
6# the terms of the BSD 3-Clause License.
8"""Calculation of Bader's QTAIM charges based on data parsed by cclib."""
9import copy
10import random
11import numpy
12import logging
13import math
15from cclib.method.calculationmethod import Method
16from cclib.method.volume import electrondensity_spin
17from cclib.parser.utils import convertor
19# Distance between two adjacent grids (sqrt[2] or sqrt[3] for uniform Cartesian grid).
20_griddist = numpy.array(
21 [
22 [[1.73205, 1.41421, 1.73205], [1.41421, 1, 1.41421], [1.73205, 1.41421, 1.73205],],
23 [[1.41421, 1, 1.41421], [1, 1, 1], [1.41421, 1, 1.41421]],
24 [[1.73205, 1.41421, 1.73205], [1.41421, 1, 1.41421], [1.73205, 1.41421, 1.73205],],
25 ],
26 dtype=float,
27)
30class MissingInputError(Exception):
31 pass
34def __cartesian_dist(pt1, pt2):
35 """ Small utility function that calculates Euclidian distance between two points
36 pt1 and pt2 are numpy arrays representing a point in Cartesian coordinates. """
37 return numpy.sqrt(numpy.einsum("ij,ij->j", pt1 - pt2, pt1 - pt2))
40class Bader(Method):
41 """Bader's QTAIM charges."""
43 # All of these are required for QTAIM charges.
44 required_attrs = ("homos", "mocoeffs", "nbasis", "gbasis")
46 def __init__(self, data, volume, progress=None, loglevel=logging.INFO, logname="Log"):
47 super(Bader, self).__init__(data, progress, loglevel, logname)
49 self.volume = volume
50 self.fragresults = None
52 if numpy.sum(self.data.coreelectrons) != 0:
53 # Pseudopotentials can cause Bader spaces to be inaccurate, as suggested by the
54 # original paper.
55 self.logger.info(
56 "It looks like pseudopotentials were used to generate this output. Please note that the Bader charges may not be accurate and may report unexpected results. Consult the original paper (doi:10.1016/j.commatsci.2005.04.010) for more information."
57 )
59 def __str__(self):
60 """Return a string representation of the object."""
61 return "Bader's QTAIM charges of {}".format(self.data)
63 def __repr__(self):
64 """Return a representation of the object."""
65 return "Bader({})".format(self.data)
67 def _check_required_attributes(self):
68 super(Bader, self)._check_required_attributes()
70 def calculate(self, indices=None, fupdate=0.05):
71 """Calculate Bader's QTAIM charges using on-grid algorithm proposed by Henkelman group
72 in doi:10.1016/j.commatsci.2005.04.010
74 Cartesian, uniformly spaced grids are assumed for this function.
75 """
77 # Obtain charge densities on the grid if it does not contain one.
78 if not numpy.any(self.volume.data):
79 self.logger.info("Calculating charge densities on the provided empty grid.")
80 if len(self.data.mocoeffs) == 1:
81 self.chgdensity = electrondensity_spin(
82 self.data, self.volume, [self.data.mocoeffs[0][: self.data.homos[0]]]
83 )
84 self.chgdensity.data *= 2
85 else:
86 self.chgdensity = electrondensity_spin(
87 self.data,
88 self.volume,
89 [
90 self.data.mocoeffs[0][: self.data.homos[0]],
91 self.data.mocoeffs[1][: self.data.homos[1]],
92 ],
93 )
95 # If charge densities are provided beforehand, log this information
96 # `Volume` object does not contain (nor rely on) information about the constituent atoms.
97 else:
98 self.logger.info("Using charge densities from the provided Volume object.")
99 self.chgdensity = self.volume
101 # Assign each grid point to Bader areas
102 self.fragresults = numpy.zeros(self.chgdensity.data.shape, "d")
103 next_index = 1
105 self.logger.info("Partitioning space into Bader areas.")
107 # Generator to iterate over the elements excluding the outermost positions
108 xshape, yshape, zshape = self.chgdensity.data.shape
109 indices = (
110 (x, y, z)
111 for x in range(1, xshape - 1)
112 for y in range(1, yshape - 1)
113 for z in range(1, zshape - 1)
114 )
116 for xindex, yindex, zindex in indices:
117 if self.fragresults[xindex, yindex, zindex] != 0:
118 # index has already been assigned for this grid point
119 continue
120 else:
121 listcoord = []
122 local_max_reached = False
124 while not local_max_reached:
125 # Here, `delta_rho` corresponds to equation 2,
126 # and `grad_rho_dot_r` corresponds to equation 1 in the aforementioned
127 # paper (doi:10.1016/j.commatsci.2005.04.010)
128 delta_rho = (
129 self.chgdensity.data[
130 xindex - 1 : xindex + 2,
131 yindex - 1 : yindex + 2,
132 zindex - 1 : zindex + 2,
133 ]
134 - self.chgdensity.data[xindex, yindex, zindex]
135 )
136 grad_rho_dot_r = delta_rho / _griddist
137 maxat = numpy.where(grad_rho_dot_r == numpy.amax(grad_rho_dot_r))
139 directions = list(zip(maxat[0], maxat[1], maxat[2]))
140 next_direction = [ind - 1 for ind in directions[0]]
142 if len(directions) > 1:
143 # when one or more directions indicate max grad (of 0), prioritize
144 # to include all points in the Bader space
145 if directions[0] == [1, 1, 1]:
146 next_direction = [ind - 1 for ind in directions[1]]
148 listcoord.append((xindex, yindex, zindex))
149 bader_candidate_index = self.fragresults[
150 xindex + next_direction[0],
151 yindex + next_direction[1],
152 zindex + next_direction[2],
153 ]
155 if bader_candidate_index != 0:
156 # Path arrived at a point that has already been assigned with an index
157 bader_index = bader_candidate_index
158 listcoord = tuple(numpy.array(listcoord).T)
159 self.fragresults[listcoord] = bader_index
161 local_max_reached = True
163 elif (
164 next_direction == [0, 0, 0]
165 or xindex + next_direction[0] == 0
166 or xindex + next_direction[0] == (len(self.chgdensity.data) - 1)
167 or yindex + next_direction[1] == 0
168 or yindex + next_direction[1] == (len(self.chgdensity.data[0]) - 1)
169 or zindex + next_direction[2] == 0
170 or zindex + next_direction[2] == (len(self.chgdensity.data[0][0]) - 1)
171 ):
172 # When next_direction is [0, 0, 0] -- local maximum
173 # Other conditions indicate that the path is heading out to edge of
174 # the grid. Here, assign new Bader space to avoid exiting the grid.
175 bader_index = next_index
176 next_index += 1
178 listcoord = tuple(numpy.array(listcoord).T)
179 self.fragresults[listcoord] = bader_index
181 local_max_reached = True
183 else:
184 # Advance to the next point according to the direction of
185 # maximum gradient
186 xindex += next_direction[0]
187 yindex += next_direction[1]
188 zindex += next_direction[2]
190 # Now try to identify each Bader region to individual atom.
191 # Try to find an area that captures enough representation
192 self.matches = numpy.zeros_like(self.data.atomnos)
193 for pos in range(len(self.data.atomcoords[-1])):
194 gridpt = numpy.round(
195 (self.data.atomcoords[-1][pos] - self.volume.origin) / self.volume.spacing
196 )
197 xgrid = int(gridpt[0])
198 ygrid = int(gridpt[1])
199 zgrid = int(gridpt[2])
200 self.matches[pos] = self.fragresults[xgrid, ygrid, zgrid]
202 assert (
203 0 not in self.matches
204 ), "Failed to assign Bader regions to atoms. Try with a finer grid. Content of Bader area matches: {}".format(
205 self.matches
206 )
207 assert len(
208 numpy.unique(self.matches) != len(self.data.atomnos)
209 ), "Failed to assign unique Bader regions to each atom. Try with a finer grid."
211 # Finally integrate the assigned Bader areas
212 self.logger.info("Creating fragcharges: array[1]")
213 self.fragcharges = numpy.zeros(len(self.data.atomcoords[-1]), "d")
215 for atom_index, baderarea_index in enumerate(self.matches):
216 self.fragcharges[atom_index] = self.chgdensity.integrate(
217 weights=(self.fragresults == baderarea_index)
218 )
220 return True