Coverage for cclib/method/ddec.py : 94%
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 DDEC charges based on data parsed by cclib."""
9import copy
10import random
11import numpy
12import logging
13import math
14import os
15import sys
17from cclib.method.calculationmethod import Method
18from cclib.method.stockholder import Stockholder
19from cclib.method.volume import electrondensity_spin
20from cclib.parser.utils import convertor
21from cclib.parser.utils import find_package
23from typing import List
26class MissingInputError(Exception):
27 pass
30class ConvergenceError(Exception):
31 pass
34class DDEC6(Stockholder):
35 """DDEC6 charges."""
37 # All of these are required for DDEC6 charges.
38 required_attrs = ("homos", "mocoeffs", "nbasis", "gbasis")
40 def __init__(
41 self,
42 data,
43 volume,
44 proatom_path=None,
45 progress=None,
46 convergence_level=1e-10,
47 max_iteration=50,
48 loglevel=logging.INFO,
49 logname="Log",
50 ):
51 """ Initialize DDEC6 object.
52 Inputs are:
53 data -- ccData object that describe target molecule.
54 volume -- Volume object that describe target Cartesian grid.
55 proatom_path -- path to proatom densities
56 (directory containing atoms.h5 in horton or c2_001_001_000_400_075.txt in chargemol)
57 convergence_level -- convergence level to use for conditioning densities in step 3
58 max_iteration -- maximum iteration to optimize phi in step 3-6
60 Note:
61 Proatom densities are used in DDEC6 algorithm in a similar way to other stockholder
62 partitioning methods. They are used as references to appropriately partition the
63 total densities (which is the density stored in cube files). Proatom densities are
64 densities obtained for single atom or ion in a radial grid that originates from
65 the atom or ion.
66 In DDEC6 algorithm, stockholder partitioning is heavily modified to ensure that the
67 total densities that are partitioned resemble the proatom densities and to prevent
68 the numerical algorithm from failing to converge.
69 """
70 super(DDEC6, self).__init__(data, volume, proatom_path, progress, loglevel, logname)
72 self.convergence_level = convergence_level
73 self.max_iteration = max_iteration
75 if numpy.sum(self.data.coreelectrons) != 0:
76 # TODO: Pseudopotentials should be added back
77 pass
79 def __str__(self):
80 """Return a string representation of the object."""
81 return "DDEC6 charges of {}".format(self.data)
83 def __repr__(self):
84 """Return a representation of the object."""
85 return "DDEC6({})".format(self.data)
87 def _check_required_attributes(self):
88 super(DDEC6, self)._check_required_attributes()
90 def _cartesian_dist(self, pt1, pt2):
91 """ Small utility function that calculates Euclidian distance between two points
92 pt1 and pt2 are numpy arrays representing a point in Cartesian coordinates. """
93 return numpy.sqrt(numpy.dot(pt1 - pt2, pt1 - pt2))
95 def _read_proatom(
96 self, directory, atom_num, charge # type = str # type = int # type = float
97 ):
98 return super(DDEC6, self)._read_proatom(directory, atom_num, charge)
100 def calculate(self, indices=None, fupdate=0.05):
101 """
102 Calculate DDEC6 charges based on doi: 10.1039/c6ra04656h paper.
103 Cartesian, uniformly spaced grids are assumed for this function.
104 """
106 super(DDEC6, self).calculate()
108 # Notify user about the total charge in the density grid
109 integrated_density = self.charge_density.integrate()
110 self.logger.info(
111 "Total charge density in the grid is {}. If this does not match what is expected, using a finer grid may help.".format(
112 integrated_density
113 )
114 )
117 # * STEP 1 *
118 # Carry out step 1 of DDEC6 algorithm [Determining reference charge value]
119 # Refer to equations 49-57 in doi: 10.1039/c6ra04656h
120 self.logger.info("Creating first reference charges. (Step 1/7)")
121 (
122 reference_charges,
123 localized_charges,
124 stockholder_charges,
125 ) = self.calculate_reference_charges()
126 self.reference_charges = [reference_charges]
127 self._localized_charges = [localized_charges]
128 self._stockholder_charges = [stockholder_charges]
130 # * STEP 2 *
131 # Load new proatom densities based on the reference charges determined in step 1.
132 self.logger.info("Creating second reference charges. (Step 2/7)")
133 self.proatom_density = []
134 self.radial_grid_r = []
135 for i, atom_number in enumerate(self.data.atomnos):
136 density, r = self._read_proatom(
137 self.proatom_path, atom_number, float(self.reference_charges[0][i])
138 )
139 self.proatom_density.append(density)
140 self.radial_grid_r.append(r)
142 # Carry out step 2 of DDEC6 algorithm [Determining ion charge value again]
143 ref, loc, stock = self.calculate_reference_charges()
144 self.reference_charges.append(ref)
145 self._localized_charges.append(loc)
146 self._stockholder_charges.append(stock)
148 # * STEP 3 *
149 # Load new proatom densities based on the reference charges determined in step 2.
150 self.proatom_density = []
151 self.radial_grid_r = []
152 self._cond_density = []
153 for i, atom_number in enumerate(self.data.atomnos):
154 density, r = self._read_proatom(
155 self.proatom_path, atom_number, float(self.reference_charges[1][i])
156 )
157 self.proatom_density.append(density)
158 self.radial_grid_r.append(r)
160 # Carry out step 3 of DDEC6 algorithm [Determine conditioned charge density and tau]
161 self.logger.info("Conditioning charge densities. (Step 3/7)")
162 self.condition_densities()
164 # Steps 4 through 7 contain similar routines. Comments the precede each step explain the
165 # differences among them.
167 # * STEP 4 *
168 # In step 4, calculate w, u, g, and h but always skip kappa updates.
169 self.logger.info("Optimizing grid weights. (Step 4/7)")
170 self._kappa = [0.0] * self.data.atomnos
171 # N_A is assigned number of electrons determined using equation 72.
172 self.N_A = []
173 # u_A is determined using equation 77.
174 self.u_A = []
176 self.N_A.append(self._calculate_w_and_u())
178 # Calculate G_A and H_A based on S4.3 in doi: 10.1039/c6ra04656h
179 self.reshape_G()
180 self.calculate_H()
182 # Update weights (w_A) using equation 96 in doi: 10.1039/c6ra04656h
183 # self._cond_density is first created in step 3 by conditioning on the total densities
184 # as described in Figure S1. Then, this quantity is updated in every step that follows
185 # until the last step in the algorithm when the weights placed on the grid is iteratively
186 # updated.
187 for atomi in range(self.data.natom):
188 self._cond_density[atomi] = math.exp(self._kappa[atomi]) * self._h[atomi]
190 # Update rho_cond for next iteration
191 self._update_rho_cond()
193 # * STEPS 5 and 6 *
194 # In step 5 and 6, calculate w and u. Then if update_kappa is found to be true, do not
195 # continue to next step. Otherwise, calculate g and h. In both cases, calculate new
196 # rho_cond.
197 steps = 5
198 self._update_kappa = False
199 while steps < 7:
200 self.logger.info("Optimizing grid weights. (Step {}/7)".format(steps))
201 self.N_A.append(self._calculate_w_and_u())
203 # Determine whether kappa needs to be updated or not based on Figure S4.2
204 # of doi: 10.1039/c6ra04656h
205 self._update_kappa = self._check_kappa()
207 if not self._update_kappa:
208 # Increment steps
209 steps = steps + 1
210 # Calculate G_A and H_A based on S4.3 in doi: 10.1039/c6ra04656h
211 self.reshape_G()
212 self.calculate_H()
214 else:
215 # `steps` not incremented in this case
216 # First update kappa based on equation 93
217 kappa_new = self._kappa - self.N_A[-1] / self.u_A[-1]
218 self._kappa = [x if x > 0 else 0.0 for x in kappa_new]
220 # Update weights (w_A) using equation 96 in doi: 10.1039/c6ra04656h
221 # self._cond_density is first created in step 3 by conditioning on the total densities
222 # as described in Figure S1. Then, this quantity is updated in every step that follows
223 # until the last step in the algorithm when the weights placed on the grid is
224 # iteratively updated.
225 for atomi in range(self.data.natom):
226 self._cond_density[atomi] = math.exp(self._kappa[atomi]) * self._h[atomi]
228 # Update rho_cond for next iteration
229 self._update_rho_cond()
231 # * STEP 7 *
232 # In step 7, calculate w and u. Then, stop to calculate reference_charges.
233 self.logger.info("Optimizing grid weights. (Step 7/7)")
234 self.N_A.append(self._calculate_w_and_u())
236 # Finally, store calculated DDEC6 charges in fragcharges
237 self.logger.info("Creating fragcharges: array[1]")
238 self.fragcharges = numpy.array(self.data.atomnos - self.N_A[-1], dtype=float)
240 def _check_kappa(self):
241 """ Return whether kappa needs to be updated or not based on Figure S4.2
242 of doi: 10.1039/c6ra04656h
243 """
244 if numpy.any([x if x < -1e-5 else 0 for x in self.N_A[-1]]):
245 return True
246 elif (
247 self._update_kappa
248 and numpy.any(numpy.diff(self.N_A)[-1] < 1e-5) # change in N_A during last cycle
249 and numpy.any(numpy.diff(self.N_A)[-2] < 1e-5) # change in N_A during last to
250 # second cycle
251 ):
252 self._kappa = [0.0 for x in self.data.atomnos]
253 return False
254 else:
255 return self._update_kappa
257 def calculate_reference_charges(self):
258 """ Calculate reference charges from proatom density and molecular density
259 [STEP 1 and 2]
261 Function returns calculated reference charges, localized charges, and stockholder
262 charges.
263 """
264 # Generator object to iterate over the grid
265 ngridx, ngridy, ngridz = self.charge_density.data.shape
266 grid_shape = (self.data.natom, ngridx, ngridy, ngridz)
268 stockholder_w = numpy.zeros(grid_shape)
269 localized_w = numpy.zeros(grid_shape)
270 self.closest_r_index = numpy.zeros(grid_shape, dtype=int)
272 indices = numpy.asanyarray(
273 tuple(
274 (x, y, z)
275 for x in range(ngridx)
276 for y in range(ngridy)
277 for z in range(ngridz)
278 )
279 )
280 coordinates = self.charge_density.coordinates(indices)
282 for atomi in range(self.data.natom):
283 # Distance of the grid from atom grid
284 self.closest_r_index[atomi] = numpy.argmin(
285 numpy.abs(
286 self.radial_grid_r[atomi][..., numpy.newaxis] - numpy.linalg.norm(self.data.atomcoords[-1][atomi] - coordinates, axis=1)
287 ),
288 axis=0
289 ).reshape((ngridx, ngridy, ngridz))
291 # Equation 54 in doi: 10.1039/c6ra04656h
292 stockholder_w[atomi] = self.proatom_density[atomi][self.closest_r_index[atomi]]
294 # Equation 55 in doi: 10.1039/c6ra04656h
295 localized_w = numpy.power(stockholder_w, 4)
297 # Equation 53 in doi: 10.1039/c6ra04656h
298 stockholder_bigW = numpy.sum(stockholder_w, axis=0)
299 localized_bigW = numpy.sum(localized_w, axis=0)
301 reference_charges = numpy.zeros((self.data.natom))
302 localizedcharges = numpy.zeros((self.data.natom))
303 stockholdercharges = numpy.zeros((self.data.natom))
305 for atomi in range(self.data.natom):
306 # Equation 52 and 51 in doi: 10.1039/c6ra04656h
307 localizedcharges[atomi] = self.data.atomnos[atomi] - self.charge_density.integrate(
308 weights=(localized_w[atomi] / localized_bigW)
309 )
310 stockholdercharges[atomi] = self.data.atomnos[atomi] - self.charge_density.integrate(
311 weights=(stockholder_w[atomi] / stockholder_bigW)
312 )
314 # In DDEC6, weights of 1/3 and 2/3 are assigned for stockholder and localized charges.
315 # (Equation 50 and 58 in doi: 10.1039/c6ra04656h)
316 reference_charges[atomi] = (stockholdercharges[atomi] / 3.0) + (
317 localizedcharges[atomi] * 2.0 / 3.0
318 )
320 return reference_charges, localizedcharges, stockholdercharges
322 def condition_densities(self):
323 """ Calculate conditioned densities
324 [STEP 3]
325 """
326 # Generator object to iterate over the grid
327 ngridx, ngridy, ngridz = self.charge_density.data.shape
329 self._rho_ref = numpy.zeros((ngridx, ngridy, ngridz))
331 for atomi in range(self.data.natom):
332 # rho_ref -- Equation 41 in doi: 10.1039/c6ra04656h
333 self._rho_ref += self.proatom_density[atomi][
334 self.closest_r_index[atomi]
335 ]
337 self._candidates_bigPhi = []
338 self._candidates_phi = []
340 # Initial conditions are detailed in Figure S1 in doi: 10.1039/c6ra04656h
341 phiAI = numpy.zeros_like(self.data.atomnos, dtype=float)
342 bigphiAI = numpy.zeros_like(self.data.atomnos, dtype=float)
343 self._y_a = []
345 for atomi in range(self.data.natom):
346 # y_a -- equation 40 in doi: 10.1039/c6ra04656h
347 self._y_a.append(self._ya(self.proatom_density, atomi))
348 # rho_A^cond -- equation S97 in doi: 10.1039/c6ra04656h
349 self._cond_density.append(
350 self._y_a[atomi] + bigphiAI[atomi] * numpy.sqrt(self._y_a[atomi])
351 )
352 # Monotonic Decrease Condition (as expressed in equation S99)
353 self._cond_density[atomi] = numpy.minimum.accumulate(self._cond_density[atomi])
354 # phi_A^I -- Equation S100 in doi: 10.1039/c6ra04656h
355 phiAI[atomi] = (
356 self._integrate_from_radial([self._cond_density[atomi]], [atomi])
357 - self.data.atomnos[atomi]
358 + self.reference_charges[-1][atomi]
359 )
361 self._candidates_bigPhi.append([bigphiAI[atomi]])
362 self._candidates_phi.append([phiAI[atomi]])
364 # Attempt to find the point where phiAI is zero iteratively
365 # Refer to S101 in doi: 10.1039/c6ra04656h
366 self._candidates_phi[atomi], self._candidates_bigPhi[atomi] = self._converge_phi(
367 phiAI[atomi], 1, atomi
368 )
370 # Perform parabolic fit to find optimized phiAI
371 # Refer to Figure S1 in doi: 10.1039/c6ra04656h
372 bigphiAI[atomi] = self._parabolic_fit(self._y_a[atomi], 1, atomi)
374 # Set final conditioned density using chosen Phi
375 self._cond_density[atomi] = self._update_phiai(
376 self._y_a[atomi], bigphiAI[atomi], atomi
377 )[1]
379 self.logger.info("Calculating tau and combined conditioned densities.")
381 # Calculate tau(r) and rho^cond(r)
382 # Refer to equation 65 and 66 in doi: 10.1039/c6ra04656h
383 # Assign rho^cond on grid using generator object
384 self.rho_cond = copy.deepcopy(self.charge_density)
385 self.rho_cond.data = numpy.zeros_like(self.rho_cond.data, dtype=float)
386 rho_cond_sqrt = numpy.zeros_like(self.rho_cond.data, dtype=float)
388 # Generator object to iterate over the grid
389 ngridx, ngridy, ngridz = self.charge_density.data.shape
391 self._leftterm = numpy.zeros((self.data.natom, ngridx, ngridy, ngridz), dtype=float)
392 # rho_cond_cartesian is rho^cond projected on Cartesian grid
393 # (used for Step 4 calculation)
394 self._rho_cond_cartesian = numpy.zeros(
395 (self.data.natom, ngridx, ngridy, ngridz), dtype=float
396 )
397 self.tau = []
399 # rho_cond -- equation 65 in doi: 10.1039/c6ra04656h
400 for atomi in range(self.data.natom):
401 self.rho_cond.data += self._cond_density[atomi][
402 self.closest_r_index[atomi]
403 ]
405 rho_cond_sqrt = numpy.sqrt(self.rho_cond.data)
407 for atomi in range(self.data.natom):
408 self.tau.append(numpy.zeros_like(self.proatom_density[atomi], dtype=float))
409 # leftterm is the first spherical average term in equation 66.
410 # <rho^cond_A(r_A) / sqrt(rho^cond(r))>
411 self._rho_cond_cartesian[atomi] = self._cond_density[atomi][
412 self.closest_r_index[atomi]
413 ]
414 self._leftterm[atomi] = self._rho_cond_cartesian[atomi] / rho_cond_sqrt
415 for radiusi in range(len(self.tau[atomi])):
416 grid_filter = self.closest_r_index[atomi] == radiusi
417 num_grid_filter = numpy.count_nonzero(grid_filter)
418 if num_grid_filter < 1:
419 self.tau[atomi][radiusi] = 0.0
420 else:
421 leftaverage = numpy.sum(grid_filter * self._leftterm[atomi]) / num_grid_filter
422 rightaverage = numpy.sum(grid_filter * rho_cond_sqrt) / num_grid_filter
423 if leftaverage < 1e-20:
424 self.tau[atomi][radiusi] = 0.0
425 else:
426 self.tau[atomi][radiusi] = numpy.divide(
427 leftaverage,
428 rightaverage,
429 out=numpy.zeros_like(leftaverage),
430 where=rightaverage != 0.0,
431 )
432 # Make tau monotonic decreasing
433 self.tau[atomi] = numpy.maximum.accumulate(self.tau[atomi][::-1])[::-1]
435 def _ya(self, proatom_density, atomi):
436 # Function that calculates Y_a^avg
437 # See Eq. 40-41 in doi: 10.1039/c6ra04656h
438 rho_ref = self._rho_ref
439 # Y_a^avg -- Equation 40 in doi: 10.1039/c6ra04656h
440 ya = numpy.zeros_like(proatom_density[atomi], dtype=float)
441 weights = self.charge_density.data / rho_ref
443 for radiusi in range(len(ya)):
444 grid_filter = self.closest_r_index[atomi] == radiusi
445 num_grid_filter = numpy.count_nonzero(grid_filter)
446 if num_grid_filter < 1:
447 ya[radiusi] = 0.0
448 else:
449 spherical_avg = numpy.sum(grid_filter * weights) / num_grid_filter
450 ya[radiusi] = proatom_density[atomi][radiusi] * spherical_avg
452 # Make y_a monotonic decreasing
453 # Refer to module_reshaping_functions.f08::77-79
454 ya = numpy.maximum.accumulate(ya[::-1])[::-1]
455 ya = numpy.minimum.accumulate(ya)
457 # Normalize y_a (see module_DDEC6_valence_iterator.f08::284)
458 nelec = self._integrate_from_radial([ya], [atomi])
459 ya *= (self.data.atomnos[atomi] - self.reference_charges[-1][atomi]) / nelec
461 return ya
463 def _calculate_w_and_u(self):
464 """ Calculate weights placed on each integration grid point
465 [STEP 4-7]
466 """
467 # From equation 67, w_A(r_A) = self._cond_density
468 # From equation 8, W(r) = self.rho_cond for the rest of this function
469 ngridx, ngridy, ngridz = self.charge_density.data.shape
471 # Evaluate rho_A(r_A) from equation 68, rho_A^avg(r_A) from equation 69,
472 # theta(r_A) from equation 70, _wavg from equation 71,
473 # N_A from equation 72, rho_wavg from equation 73, and u_A from equation 77.
474 self._rho_A = []
475 self._rho_A_avg = []
476 self._theta = []
477 self._wavg = []
478 N_A = []
479 u_A = []
480 self.rho_wavg = []
482 for atomi in range(self.data.natom):
483 self._rho_A.append(copy.deepcopy(self.charge_density))
484 # Equation 68
485 self._rho_A[atomi].data = numpy.divide(
486 self.charge_density.data * self._rho_cond_cartesian[atomi],
487 self.rho_cond.data,
488 out=numpy.zeros_like(self.charge_density.data, dtype=float),
489 where=self.rho_cond.data != 0,
490 )
491 self._rho_A_avg.append(numpy.zeros_like(self.proatom_density[atomi], dtype=float))
492 self._theta.append(numpy.zeros_like(self.proatom_density[atomi], dtype=float))
493 self._wavg.append(numpy.zeros_like(self.proatom_density[atomi], dtype=float))
494 # Equation 69, 70, and 71
495 self._rho_A_avg[atomi] = self._spherical_average_from_cartesian(
496 self._rho_A[atomi].data, atomi, self.radial_grid_r[atomi]
497 )
498 self._rho_A_avg[atomi] = numpy.maximum.accumulate(self._rho_A_avg[atomi][::-1])[::-1]
499 self._theta[atomi] = self._spherical_average_from_cartesian(
500 (
501 1
502 - numpy.divide(
503 self._rho_cond_cartesian[atomi],
504 self.rho_cond.data,
505 out=numpy.zeros_like(self.rho_cond.data, dtype=float),
506 where=self.rho_cond.data != 0,
507 )
508 )
509 * self._rho_A[atomi].data,
510 atomi,
511 self.radial_grid_r[atomi],
512 )
513 self._theta[atomi] = numpy.maximum.accumulate(self._theta[atomi][::-1])[::-1]
514 self._wavg[atomi] = self._spherical_average_from_cartesian(
515 numpy.divide(
516 self._rho_cond_cartesian[atomi],
517 self.rho_cond.data,
518 out=numpy.zeros_like(self.rho_cond.data, dtype=float),
519 where=self.rho_cond.data != 0,
520 ),
521 atomi,
522 self.radial_grid_r[atomi],
523 )
524 self._wavg[atomi] = numpy.maximum.accumulate(self._wavg[atomi][::-1])[::-1]
525 # Equation 72, 73, and 77
526 N_A.append(self._rho_A[atomi].integrate())
527 self.rho_wavg.append(
528 (self._theta[atomi] + self._rho_A_avg[atomi] * self._wavg[atomi] / 5.0)
529 / (1.0 - (4.0 / 5.0) * self._wavg[atomi])
530 )
531 u_A.append(self._integrate_from_radial([self._theta[atomi]], [atomi]))
532 self.u_A.append(u_A)
534 return N_A
536 def reshape_G(self):
537 """ Calculate G_A(r_A) and reshape densities
539 This is a quantity introduced in DDEC6 as a constraint preventing the tails from being
540 too diffuse.
541 [STEP 4-7]
542 """
543 self._candidates_bigPhi = []
544 self._candidates_phi = []
546 # Initial conditions are detailed in Figure S3 in doi: 10.1039/c6ra04656h
547 phiAII = numpy.zeros_like(self.data.atomnos, dtype=float)
548 bigphiAII = numpy.zeros_like(self.data.atomnos, dtype=float)
549 self._g = []
550 self._eta = []
552 for atomi in range(self.data.natom):
553 # G_A -- equation S102 in doi: 10.1039/c6ra04656h
554 self._g.append(numpy.zeros_like(self.proatom_density[atomi], dtype=float))
555 self._g[atomi] = self.rho_wavg[atomi] + bigphiAII[atomi] * numpy.sqrt(
556 self.rho_wavg[atomi]
557 )
558 # Exponential constraint (as expressed in equation S105)
559 self._eta.append((1 - (self.tau[atomi]) ** 2) * 1.75 * convertor(1, "Angstrom", "bohr"))
560 exp_applied = self._g[atomi][:-1] * numpy.exp(
561 -1 * self._eta[atomi][1:] * numpy.diff(self.radial_grid_r[atomi])
562 )
563 for radiusi in range(1, len(self._g[atomi])):
564 self._g[atomi][radiusi] = min(self._g[atomi][radiusi], exp_applied[radiusi - 1])
565 # phi_A^II -- Equation S106 in doi: 10.1039/c6ra04656h
566 phiAII[atomi] = self._integrate_from_radial(
567 [self._g[atomi] - self.rho_wavg[atomi]], [atomi]
568 )
570 self._candidates_bigPhi.append([bigphiAII[atomi]])
571 self._candidates_phi.append([phiAII[atomi]])
573 # Attempt to find the point where phiAI is zero iteratively
574 # Refer to S101 in doi: 10.1039/c6ra04656h
575 self._candidates_phi[atomi], self._candidates_bigPhi[atomi] = self._converge_phi(
576 phiAII[atomi], 1, atomi
577 )
579 # Perform parabolic fit to find optimized phiAI
580 # Refer to Figure S1 in doi: 10.1039/c6ra04656h
581 bigphiAII[atomi] = self._parabolic_fit(self.rho_wavg[atomi], 2, atomi)
583 # Set final G_A value using chosen Phi
584 self._g[atomi] = self._update_phiaii(self.rho_wavg[atomi], bigphiAII[atomi], atomi)[1]
586 def calculate_H(self):
587 """ Calculate H_A(r_A)
589 This is a quantity introduced in DDEC6 as a constraint preventing the tails from being
590 too contracted.
591 [STEP 4-7]
592 """
593 self._h = []
594 for atomi in range(len(self._g)):
595 # First set H_est as G_A
596 self._h.append(self._g[atomi])
598 # Determine eta_upper using equation 86 in doi: 10.1039/c6ra04656h
599 # and apply upper limit using equation 91.
600 temp = (
601 1 - (self.tau[atomi]) ** 2 + self.convergence_level
602 ) # convergence_level is added to avoid divide-by-zero in next line for highly polar molecules.
603 eta = 2.5 * convertor(1, "Angstrom", "bohr") / temp
604 exp_applied = self._h[atomi][:-1] * numpy.exp(
605 -1 * eta[1:] * numpy.diff(self.radial_grid_r[atomi])
606 )
607 for radiusi in range(1, len(self._h[atomi])):
608 self._h[atomi][radiusi] = max(self._h[atomi][radiusi], exp_applied[radiusi - 1])
610 # Normalize using equation 92 in doi: 10.1039/c6ra04656h.
611 self._h[atomi] = (
612 self._h[atomi]
613 * self._integrate_from_radial([self._g[atomi]], [atomi])
614 / self._integrate_from_radial([self._h[atomi]], [atomi])
615 )
617 def _update_phiai(self, ya, bigphiAI, atomi):
618 # Update phi^a_i and quantity that directly follows (cond_density) in each step of
619 # iterative optimization. (Refer to Figure S1)
620 # Re-evaluate cond_density
621 if isinstance(bigphiAI, float) and not numpy.isinf(bigphiAI):
622 cond_density = ya + bigphiAI * numpy.sqrt(ya)
623 else:
624 cond_density = ya
626 # Monotonic Decrease Condition
627 cond_density = numpy.minimum.accumulate(cond_density)
629 # Re-evaluate phi_AI
630 phiAI = (
631 self._integrate_from_radial([cond_density], [atomi])
632 - self.data.atomnos[atomi]
633 + self.reference_charges[-1][atomi]
634 )
636 return phiAI, cond_density
638 def _update_phiaii(self, rhowavg, bigphiAI, atomi):
639 # Update phi^a_ii and quantity that directly follows (G_A) in each step of
640 # iterative optimization. (Refer to Figure S3)
641 # Re-evaluate g_a
642 # Equations can be found in Figure S3.
643 ga = rhowavg + bigphiAI * numpy.sqrt(rhowavg)
645 # Exponential Decrease Condition
646 exp_applied = ga[:-1] * numpy.exp(
647 -1 * self._eta[atomi][1:] * numpy.diff(self.radial_grid_r[atomi])
648 )
649 for radiusi in range(1, len(ga)):
650 ga[radiusi] = min(ga[radiusi], exp_applied[radiusi - 1])
652 # Re-evaluate phi_AII
653 phiAII = self._integrate_from_radial([ga - rhowavg], [atomi])
655 return phiAII, ga
657 def _integrate_from_radial(self, radial_density_list, atom_list):
658 # Function that reads in list of radial densities, projects it on Cartesian grid,
659 # and returns integrated value
660 grid = copy.deepcopy(self.charge_density)
661 grid.data = numpy.zeros_like(grid.data, dtype=float)
663 for density, atomi in zip(radial_density_list, atom_list):
664 grid.data += density[self.closest_r_index[atomi]]
666 return grid.integrate()
668 def _spherical_average_from_cartesian(self, cartesian_grid, atom_index, radius_list):
669 spherical_average = numpy.zeros(len(radius_list))
670 for radiusi in range(len(radius_list)):
671 grid_filter = self.closest_r_index[atom_index] == radiusi
672 num_grid_filter = numpy.count_nonzero(grid_filter)
673 if num_grid_filter < 1:
674 average = 0.0
675 else:
676 average = numpy.sum(grid_filter * cartesian_grid) / num_grid_filter
677 if average < self.convergence_level:
678 average = 0.0
679 spherical_average[radiusi] = average
680 return spherical_average
682 def _update_rho_cond(self):
683 # Update total weights on Cartesian grid using equation 65 in doi: 10.1039/c6ra04656h
684 ngridx, ngridy, ngridz = self.charge_density.data.shape
686 self.rho_cond.data = numpy.zeros_like(self.rho_cond.data, dtype=float)
687 self._rho_cond_cartesian = numpy.zeros(
688 (self.data.natom, ngridx, ngridy, ngridz), dtype=float
689 )
691 for atomi in range(self.data.natom):
692 self.rho_cond.data += self._cond_density[atomi][
693 self.closest_r_index[atomi]
694 ]
695 self._rho_cond_cartesian[atomi] = self._cond_density[atomi][
696 self.closest_r_index[atomi]
697 ]
699 def _converge_phi(self, phiA, superscript, atomi):
700 """ Update phi until it is positive.
701 This is used in step 3 (for phi_A^I) and in steps 4-6 (for phi_A^II).
703 --- Inputs ---
704 phiA Either phi_A^I or phi_A^II
705 superscript 1 when calculating phi_I (STEP 3)
706 2 when calculating phi_II (STEPS 4-6)
707 atomi Index of target atom as in ccData object
709 Refer to Equation S101, Figure S1 and S3 for an overview.
710 """
711 # Initial value of bigphi is zero (Equation S100)
712 bigphiA = 0.0
714 # List to store candidate values for parabolic fitting
715 candidates_phi = [phiA]
716 candidates_bigphi = [bigphiA]
718 while phiA <= 0:
719 # Iterative algorithm until convergence
720 # Refer to S101 in doi: 10.1039/c6ra04656h
721 if superscript == 1:
722 temp = self._integrate_from_radial([numpy.sqrt(self._y_a[atomi])], [atomi])
723 elif superscript == 2:
724 temp = self._integrate_from_radial([numpy.sqrt(self.rho_wavg[atomi])], [atomi])
726 bigphiA = 2 * bigphiA - phiA / temp
728 # When Phi is updated, related quantities are updated as well
729 # Refer to S100 in doi: 10.1039/c6ra04656h [Step 3]
730 # or S102 in doi: 10.1039/c6ra04656h [Steps 4-6]
732 if superscript == 1:
733 phiA, self._cond_density[atomi] = self._update_phiai(
734 self._y_a[atomi], bigphiA, atomi
735 )
736 elif superscript == 2:
737 phiA, self._g[atomi] = self._update_phiaii(self.rho_wavg[atomi], bigphiA, atomi)
739 candidates_phi.append(phiA)
740 candidates_bigphi.append(bigphiA)
742 return candidates_phi, candidates_bigphi
744 def _parabolic_fit(self, pseudodensity, superscript, atomi):
745 """ Optimize phi using parabolic fitting.
746 This is used in step 3 (for phi_A^I) and in steps 4-6 (for phi_A^II).
748 --- Inputs ---
749 phiA Either phi_A^I or phi_A^II
750 superscript 1 when calculating phi_I (STEP 3)
751 2 when calculating phi_II (STEPS 4-6)
752 atomi Index of target atom as in ccData object
754 Refer to Figure S1 and S3 for an overview.
755 """
756 # Set update methods for phi_A^I or phi_A^II
757 if superscript == 1:
759 def update(pdens, bigPhi, atomi):
760 return self._update_phiai(pdens, bigPhi, atomi)
762 elif superscript == 2:
764 def update(pdens, bigPhi, atomi):
765 return self._update_phiaii(pseudodensity, bigPhi, atomi)
767 # lowerbigPhi is bigPhi that yields biggest negative phi.
768 # upperbigPhi is bigPhi that yields smallest positive phi.
769 # The point here is to find two phi values that are closest to zero (from positive side
770 # and negative side respectively).
771 self._candidates_phi[atomi] = numpy.array(self._candidates_phi[atomi], dtype=float)
772 self._candidates_bigPhi[atomi] = numpy.array(self._candidates_bigPhi[atomi], dtype=float)
773 if numpy.count_nonzero(self._candidates_phi[atomi] < 0) > 0:
774 # If there is at least one candidate phi that is negative
775 lower_ind = numpy.where(
776 self._candidates_phi[atomi]
777 == self._candidates_phi[atomi][self._candidates_phi[atomi] < 0].max()
778 )[0][0]
779 lowerbigPhi = self._candidates_bigPhi[atomi][lower_ind]
780 lowerphi = self._candidates_phi[atomi][lower_ind]
781 else: # assign some large negative number otherwise
782 lowerbigPhi = numpy.NINF
783 lowerphi = numpy.NINF
784 if numpy.count_nonzero(self._candidates_phi[atomi] > 0) > 0:
785 # If there is at least one candidate phi that is positive
786 upper_ind = numpy.where(
787 self._candidates_phi[atomi]
788 == self._candidates_phi[atomi][self._candidates_phi[atomi] > 0].min()
789 )[0][0]
790 upperbigPhi = self._candidates_bigPhi[atomi][upper_ind]
791 upperphi = self._candidates_phi[atomi][upper_ind]
792 else: # assign some large positive number otherwise
793 upperbigPhi = numpy.PINF
794 upperphi = numpy.PINF
796 for iteration in range(self.max_iteration):
797 # Flow diagram on Figure S1 in doi: 10.1039/c6ra04656h details the procedure.
798 # Find midpoint between positive bigPhi that yields phi closest to zero and negative
799 # bigPhi closest to zero. Then, evaluate phi.
800 # This can be thought as linear fitting compared to parabolic fitting below.
801 midbigPhi = (lowerbigPhi + upperbigPhi) / 2.0
802 midphi = update(pseudodensity, midbigPhi, atomi)[0]
803 # Exit conditions -- if any of three phi values are within the convergence level.
804 if abs(lowerphi) < self.convergence_level:
805 return lowerbigPhi
806 elif abs(upperphi) < self.convergence_level:
807 return upperbigPhi
808 elif abs(midphi) < self.convergence_level:
809 return midbigPhi
811 # Parabolic fitting as described on Figure S1 in doi: 10.1039/c6ra04656h
812 # Type casting here converts from size 1 numpy.ndarray to float
813 xpts = numpy.array(
814 [float(lowerbigPhi), float(midbigPhi), float(upperbigPhi)], dtype=float
815 )
816 ypts = numpy.array([float(lowerphi), float(midphi), float(upperphi)], dtype=float)
817 fit = numpy.polyfit(xpts, ypts, 2)
818 roots = numpy.roots(fit) # max two roots (bigPhi) from parabolic fitting
820 # Find phi for two bigPhis that were obtained from parabolic fitting.
821 belowphi = update(pseudodensity, roots.min(), atomi)[0]
822 abovephi = update(pseudodensity, roots.max(), atomi)[0]
824 # If phi values from parabolically fitted bigPhis lie within the convergence level,
825 # exit the iterative algorithm.
826 if abs(abovephi) < self.convergence_level:
827 return roots.min()
828 elif abs(belowphi) < self.convergence_level:
829 return roots.max()
830 else:
831 # Otherwise, corrected phi value is obtained in a way that cuts the numerical
832 # search domain in half in each iteration.
833 if 3 * abs(abovephi) < abs(belowphi):
834 corbigPhi = roots.max() - 2.0 * abovephi * (roots.max() - roots.min()) / (
835 abovephi - belowphi
836 )
837 elif 3 * abs(belowphi) < abs(abovephi):
838 corbigPhi = roots.min() - 2.0 * belowphi * (roots.max() - roots.min()) / (
839 abovephi - belowphi
840 )
841 else:
842 corbigPhi = (roots.max() + roots.min()) / 2.0
843 # New candidates of phi and bigPhi are determined as bigPhi yielding largest
844 # negative phi and bigPhi yielding smallest positve phi. This is analogous to how
845 # the first candidiate phi values are evaluated.
846 corphi = update(pseudodensity, corbigPhi, atomi)[0]
847 self._candidates_bigPhi[atomi] = numpy.array(
848 [lowerbigPhi, midbigPhi, upperbigPhi, roots.max(), roots.min(), corbigPhi,],
849 dtype=float,
850 )
851 self._candidates_phi[atomi] = numpy.array(
852 [lowerphi, midphi, upperphi, abovephi, belowphi, corphi], dtype=float
853 )
855 # Set new upperphi and lowerphi
856 lower_ind = numpy.where(
857 self._candidates_phi[atomi]
858 == self._candidates_phi[atomi][self._candidates_phi[atomi] < 0].max()
859 )[0][0]
860 upper_ind = numpy.where(
861 self._candidates_phi[atomi]
862 == self._candidates_phi[atomi][self._candidates_phi[atomi] > 0].min()
863 )[0][0]
865 lowerphi = self._candidates_phi[atomi][lower_ind]
866 upperphi = self._candidates_phi[atomi][upper_ind]
868 # If new lowerphi or upperphi values are within convergence level, exit the
869 # iterative algorithm. Otherwise, start new linear/parabolic fitting.
870 if abs(lowerphi) < self.convergence_level:
871 return self._candidates_bigPhi[atomi][lower_ind]
872 elif abs(upperphi) < self.convergence_level:
873 return self._candidates_bigPhi[atomi][upper_ind]
874 else:
875 # Fitting needs to continue in this case.
876 lowerbigPhi = self._candidates_bigPhi[atomi][lower_ind]
877 lowerphi = self._candidates_phi[atomi][lower_ind]
878 upperbigPhi = self._candidates_bigPhi[atomi][upper_ind]
879 upperphi = self._candidates_phi[atomi][upper_ind]
881 # Raise Exception if convergence is not achieved within max_iteration.
882 raise ConvergenceError("Iterative conditioning failed to converge.")