Coverage for cclib/method/stockholder.py : 80%
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"""Stockholder partitioning based on cclib data."""
9import copy
10import random
11import numpy
12import logging
13import math
14import os
15import sys
17from cclib.method.calculationmethod import Method
18from cclib.method.volume import electrondensity_spin
19from cclib.parser.utils import convertor
20from cclib.parser.utils import find_package
22from typing import List
25class MissingInputError(Exception):
26 pass
29class Stockholder(Method):
30 """An abstract base class for stockholder-type methods."""
32 # All of these are required for stockholder-type charges.
33 required_attrs = ("homos", "mocoeffs", "nbasis", "gbasis")
35 def __init__(
36 self,
37 data,
38 volume,
39 proatom_path=None,
40 progress=None,
41 loglevel=logging.INFO,
42 logname="Log",
43 ):
44 """ Initialize Stockholder-type method object.
45 Inputs are:
46 data -- ccData object that describe target molecule.
47 volume -- Volume object that describe target Cartesian grid.
48 proatom_path -- path to proatom densities
49 (directory containing atoms.h5 in horton or c2_001_001_000_400_075.txt in chargemol)
50 """
51 super(Stockholder, self).__init__(data, progress, loglevel, logname)
53 self.volume = volume
54 self.proatom_path = proatom_path
56 # Check whether proatom_path is a valid directory or not.
57 assert os.path.isdir(
58 proatom_path
59 ), "Directory that contains proatom densities should be added as an input."
61 # Read in reference charges.
62 self.proatom_density = []
63 self.radial_grid_r = []
64 for atom_number in self.data.atomnos:
65 density, r = self._read_proatom(self.proatom_path, atom_number, 0)
66 self.proatom_density.append(density)
67 self.radial_grid_r.append(r)
69 def __str__(self):
70 """Return a string representation of the object."""
71 return "Stockholder"
73 def __repr__(self):
74 """Return a representation of the object."""
75 return "Stockholder"
77 def _check_required_attributes(self):
78 super(Stockholder, self)._check_required_attributes()
80 def _read_proatom(
81 self, directory, atom_num, charge # type = str # type = int # type = float
82 ):
83 # type: (...) -> numpy.ndarray, numpy.ndarray
84 """Return a list containing proatom reference densities."""
85 # TODO: Treat calculations with psuedopotentials
86 # TODO: Modify so that proatom densities are read only once for horton
87 # [https://github.com/cclib/cclib/pull/914#discussion_r464039991]
88 # File name format:
89 # ** Chargemol **
90 # c2_[atom number]_[nuclear charge]_[electron count]_[cutoff radius]_[# shells]
91 # ** Horton **
92 # atoms.h5
93 # File format:
94 # Starting from line 13, each line contains the charge densities for each shell
95 # If `charge` is not an integer, proatom densities have to be linearly interpolated between
96 # the densities of the ion/atom with floor(charge) and ceiling(charge)
97 charge_floor = int(math.floor(charge))
98 charge_ceil = int(math.ceil(charge))
100 chargemol_path_floor = os.path.join(
101 directory,
102 "c2_{:03d}_{:03d}_{:03d}_500_100.txt".format(
103 atom_num, atom_num, atom_num - charge_floor
104 ),
105 )
106 chargemol_path_ceil = os.path.join(
107 directory,
108 "c2_{:03d}_{:03d}_{:03d}_500_100.txt".format(
109 atom_num, atom_num, atom_num - charge_ceil
110 ),
111 )
112 horton_path = os.path.join(directory, "atoms.h5")
114 if os.path.isfile(chargemol_path_floor) or os.path.isfile(chargemol_path_ceil):
115 # Use chargemol proatom densities
116 # Each shell is .05 angstroms apart (uniform).
117 # *scalefactor* = 10.58354497764173 bohrs in module_global_parameter.f08
118 if atom_num <= charge_floor:
119 density_floor = numpy.array([0])
120 else:
121 density_floor = numpy.loadtxt(chargemol_path_floor, skiprows=12, dtype=float)
122 if atom_num >= charge_ceil:
123 density_ceil = numpy.array([0])
124 else:
125 density_ceil = numpy.loadtxt(chargemol_path_ceil, skiprows=12, dtype=float)
127 density = (charge_ceil - charge) * density_floor + (
128 charge - charge_floor
129 ) * density_ceil
130 radiusgrid = numpy.arange(1, len(density) + 1) * 0.05
132 elif os.path.isfile(horton_path):
133 # Use horton proatom densities
134 assert find_package("h5py"), "h5py is needed to read in proatom densities from horton."
136 import h5py
138 with h5py.File(horton_path, "r") as proatomdb:
139 if atom_num <= charge_floor:
140 density_floor = numpy.array([0])
141 radiusgrid = numpy.array([0])
142 else:
143 keystring_floor = "Z={}_Q={:+d}".format(atom_num, charge_floor)
144 density_floor = numpy.asanyarray(list(proatomdb[keystring_floor]["rho"]))
146 # gridspec is specification of integration grid for proatom densities in horton.
147 # Example -- ['PowerRTransform', '1.1774580743206259e-07', '20.140888089596444', '41']
148 # is constructed using PowerRTransform grid
149 # with rmin = 1.1774580743206259e-07
150 # rmax = 20.140888089596444
151 # and ngrid = 41
152 # PowerRTransform is default in horton-atomdb.py.
153 gridtype, gridmin, gridmax, gridn = (
154 proatomdb[keystring_floor].attrs["rtransform"].split()
155 )
156 gridmin = convertor(float(gridmin), "bohr", "Angstrom")
157 gridmax = convertor(float(gridmax), "bohr", "Angstrom")
158 gridn = int(gridn)
159 if isinstance(gridtype, bytes):
160 gridtype = gridtype.decode("UTF-8")
162 # First verify that it is one of recognized grids
163 assert gridtype in [
164 "LinearRTransform",
165 "ExpRTransform",
166 "PowerRTransform",
167 ], "Grid type not recognized."
169 if gridtype == "LinearRTransform":
170 # Linear transformation. r(t) = rmin + t*(rmax - rmin)/(npoint - 1)
171 gridcoeff = (gridmax - gridmin) / (gridn - 1)
172 radiusgrid = gridmin + numpy.arange(1, gridn + 1) * gridcoeff
173 elif gridtype == "ExpRTransform":
174 # Exponential transformation. r(t) = rmin*exp(t*log(rmax/rmin)/(npoint - 1))
175 gridcoeff = math.log(gridmax / gridmin) / (gridn - 1)
176 radiusgrid = gridmin * numpy.exp(numpy.arange(1, gridn + 1) * gridcoeff)
177 elif gridtype == "PowerRTransform":
178 # Power transformation. r(t) = rmin*t^power
179 # with power = log(rmax/rmin)/log(npoint)
180 gridcoeff = math.log(gridmax / gridmin) / math.log(gridn)
181 radiusgrid = gridmin * numpy.power(numpy.arange(1, gridn + 1), gridcoeff)
183 if atom_num <= charge_ceil:
184 density_ceil = numpy.array([0])
185 else:
186 keystring_ceil = "Z={}_Q={:+d}".format(atom_num, charge_ceil)
187 density_ceil = numpy.asanyarray(list(proatomdb[keystring_ceil]["rho"]))
189 density = (charge_ceil - charge) * density_floor + (
190 charge - charge_floor
191 ) * density_ceil
193 del h5py
195 else:
196 raise MissingInputError("Pro-atom densities were not found in the specified path.")
198 if charge == charge_floor:
199 density = density_floor
201 return density, radiusgrid
203 def calculate(self, indices=None, fupdate=0.05):
204 """ Charge density on a Cartesian grid is a common routine required for Stockholder-type
205 and related methods. This abstract class prepares the grid if input Volume object
206 is empty.
207 """
208 # Obtain charge densities on the grid if it does not contain one.
209 if not numpy.any(self.volume.data):
210 self.logger.info("Calculating charge densities on the provided empty grid.")
211 if len(self.data.mocoeffs) == 1:
212 self.charge_density = electrondensity_spin(
213 self.data, self.volume, [self.data.mocoeffs[0][: self.data.homos[0] + 1]]
214 )
215 self.charge_density.data *= 2
216 else:
217 self.charge_density = electrondensity_spin(
218 self.data,
219 self.volume,
220 [
221 self.data.mocoeffs[0][: self.data.homos[0] + 1],
222 self.data.mocoeffs[1][: self.data.homos[1] + 1],
223 ],
224 )
225 # If charge densities are provided beforehand, log this information
226 # `Volume` object does not contain (nor rely on) information about the constituent atoms.
227 else:
228 self.logger.info("Using charge densities from the provided Volume object.")
229 self.charge_density = self.volume