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. 

7 

8"""Calculation of Hirshfeld charges based on data parsed by cclib.""" 

9import copy 

10import random 

11import numpy 

12import logging 

13import math 

14import os 

15import sys 

16 

17from cclib.method.calculationmethod import Method 

18from cclib.method.stockholder import Stockholder 

19from cclib.method.volume import electrondensity_spin 

20from cclib.parser.utils import find_package 

21 

22from typing import List 

23 

24 

25class MissingInputError(Exception): 

26 pass 

27 

28 

29class ConvergenceError(Exception): 

30 pass 

31 

32 

33class Hirshfeld(Stockholder): 

34 """Hirshfeld charges.""" 

35 

36 # All of these are required for DDEC6 charges. 

37 required_attrs = ("homos", "mocoeffs", "nbasis", "gbasis") 

38 

39 def __init__( 

40 self, 

41 data, 

42 volume, 

43 proatom_path=None, 

44 progress=None, 

45 loglevel=logging.INFO, 

46 logname="Log", 

47 ): 

48 """ Initialize Hirshfeld object. 

49 Inputs are: 

50 data -- ccData object that describe target molecule. 

51 volume -- Volume object that describe target Cartesian grid. 

52 proatom_path -- path to proatom densities 

53 (directory containing atoms.h5 in horton or c2_001_001_000_400_075.txt in chargemol) 

54 """ 

55 super(Hirshfeld, self).__init__(data, volume, proatom_path, progress, loglevel, logname) 

56 

57 def __str__(self): 

58 """Return a string representation of the object.""" 

59 return "Hirshfeld charges of {}".format(self.data) 

60 

61 def __repr__(self): 

62 """Return a representation of the object.""" 

63 return "Hirshfeld({})".format(self.data) 

64 

65 def _check_required_attributes(self): 

66 super(Hirshfeld, self)._check_required_attributes() 

67 

68 def _cartesian_dist(self, pt1, pt2): 

69 """Small utility function that calculates Euclidian distance between two points. 

70  

71 Arguments pt1 and pt2 are NumPy arrays representing points in Cartesian coordinates. 

72 """ 

73 return numpy.sqrt(numpy.dot(pt1 - pt2, pt1 - pt2)) 

74 

75 def _read_proatom( 

76 self, directory, atom_num, charge # type = str # type = int # type = float 

77 ): 

78 return super(Hirshfeld, self)._read_proatom(directory, atom_num, charge) 

79 

80 def calculate(self): 

81 """Calculate Hirshfeld charges.""" 

82 super(Hirshfeld, self).calculate() 

83 

84 # Generator object to iterate over the grid. 

85 ngridx, ngridy, ngridz = self.charge_density.data.shape 

86 indices = ( 

87 (i, x, y, z) 

88 for i in range(self.data.natom) 

89 for x in range(ngridx) 

90 for y in range(ngridy) 

91 for z in range(ngridz) 

92 ) 

93 grid_shape = (self.data.natom, ngridx, ngridy, ngridz) 

94 

95 stockholder_w = numpy.zeros(grid_shape) 

96 self.closest_r_index = numpy.zeros(grid_shape, dtype=int) 

97 

98 for atomi, xindex, yindex, zindex in indices: 

99 # Distance of the grid from atom grid. 

100 dist_r = self._cartesian_dist( 

101 self.data.atomcoords[-1][atomi], 

102 self.charge_density.coordinates([xindex, yindex, zindex]), 

103 ) 

104 self.closest_r_index[atomi][xindex][yindex][zindex] = numpy.abs( 

105 self.radial_grid_r[atomi] - dist_r 

106 ).argmin() 

107 

108 # stockholder_w is radial proatom density projected on Cartesian grid 

109 stockholder_w[atomi][xindex][yindex][zindex] = self.proatom_density[atomi][ 

110 self.closest_r_index[atomi][xindex][yindex][zindex] 

111 ] 

112 

113 # Equation 3 in doi:10.1007/BF01113058 

114 # stockholder_bigW represents the denominator in this equation 

115 # \sum{eta_chi(r_chi)}_chi 

116 # where eta is proatom density, chi is index of atoms in the molecule, and r_chi is 

117 # radial distance from center of atom chi. 

118 stockholder_bigW = numpy.sum(stockholder_w, axis=0) 

119 

120 self.fragcharges = numpy.zeros((self.data.natom)) 

121 self.logger.info("Creating fragcharges: array[1]") 

122 

123 for atomi in range(self.data.natom): 

124 # Equation 1 in doi:10.1007/BF01113058 

125 # The weights supplied for integrate function correspond to W_A in the equation. 

126 # Q_A = Z_A - \integral(W_A(r) * q_A(r) dr) 

127 # where Q_A is Hirshfeld charges, W_A is weights determined on each grid and on each 

128 # atom, and q_A is total charge densities on the Cartesian grid. 

129 self.fragcharges[atomi] = self.data.atomnos[atomi] - self.charge_density.integrate( 

130 weights=(stockholder_w[atomi] / stockholder_bigW) 

131 ) 

132 

133 return True