Coverage for gpaw/analyse/hirshfeld.py: 98%
107 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import numpy as np
3from ase.units import Bohr
4from ase.utils.timing import Timer
6from gpaw.density import RealSpaceDensity
7from gpaw.lfc import BasisFunctions
8from gpaw.setup import Setups
9from gpaw.xc import XC
10from gpaw.utilities.tools import coordinates
11from gpaw.utilities.partition import AtomPartition
12from gpaw.mpi import world
13from gpaw.io.logger import GPAWLogger
16class HirshfeldDensity(RealSpaceDensity):
17 """Density as sum of atomic densities."""
19 def __init__(self, calculator, log=None):
20 self.calculator = calculator
21 dens = calculator.density
22 super().__init__(dens.gd, dens.finegd,
23 dens.nspins, collinear=True, charge=0.0,
24 stencil=getattr(dens, 'stencil', 3),
25 redistributor=dens.redistributor)
26 self.log = GPAWLogger(world=world)
27 if log is None:
28 self.log.fd = None
29 else:
30 self.log.fd = log
32 def set_positions(self, spos_ac, atom_partition):
33 """HirshfeldDensity builds a hack density object to calculate
34 all electron density
35 of atoms. This methods overrides the parallel distribution of
36 atomic density matrices
37 in density.py"""
38 self.atom_partition = atom_partition
39 self.atomdist = self.redistributor.get_atom_distributions(spos_ac)
40 self.nct.set_positions(spos_ac)
41 self.ghat.set_positions(spos_ac)
42 self.mixer.reset()
43 # self.nt_sG = None
44 self.nt_sg = None
45 self.nt_g = None
46 self.rhot_g = None
47 self.Q_aL = None
48 self.nct_G = self.gd.zeros()
49 self.nct.add(self.nct_G, 1.0 / self.nspins)
51 def get_density(self, atom_indices=None, gridrefinement=2):
52 """Get sum of atomic densities from the given atom list.
54 Parameters
55 ----------
56 atom_indices : list_like
57 All atoms are taken if the list is not given.
58 gridrefinement : 1, 2, 4
59 Gridrefinement given to get_all_electron_density
61 Returns
62 -------
63 type
64 spin summed density, grid_descriptor
65 """
67 all_atoms = self.calculator.get_atoms()
68 if atom_indices is None:
69 atom_indices = range(len(all_atoms))
71 # select atoms
72 atoms = self.calculator.get_atoms()[atom_indices]
73 spos_ac = atoms.get_scaled_positions()
74 Z_a = atoms.get_atomic_numbers()
76 par = self.calculator.parameters
77 setups = Setups(Z_a, par.setups, par.basis,
78 XC(par.xc),
79 world=self.calculator.wfs.world)
81 # initialize
82 self.initialize(setups,
83 self.calculator.timer,
84 np.zeros((len(atoms), 3)), False)
85 self.set_mixer(None)
86 rank_a = self.gd.get_ranks_from_positions(spos_ac)
87 self.set_positions(spos_ac, AtomPartition(self.gd.comm, rank_a))
88 basis_functions = BasisFunctions(self.gd,
89 [setup.basis_functions_J
90 for setup in self.setups],
91 cut=True)
92 basis_functions.set_positions(spos_ac)
93 self.initialize_from_atomic_densities(basis_functions)
95 aed_sg, gd = self.get_all_electron_density(atoms,
96 gridrefinement)
97 return aed_sg.sum(axis=0), gd
100class HirshfeldPartitioning:
101 """Partion space according to the Hirshfeld method.
103 After: F. L. Hirshfeld Theoret. Chim.Acta 44 (1977) 129-138
104 """
106 def __init__(self, calculator, density_cutoff=1.e-12):
107 self.calculator = calculator
108 self.density_cutoff = density_cutoff
110 if hasattr(self.calculator, 'timer'):
111 self.timer = self.calculator.timer
112 else:
113 self.timer = Timer()
115 def initialize(self):
116 with self.timer('HirshfeldPartitioning initialize'):
117 self.atoms = self.calculator.get_atoms()
118 self.hdensity = HirshfeldDensity(self.calculator)
119 density_g, gd = self.hdensity.get_density()
120 self.invweight_g = 0. * density_g
121 density_ok = np.where(density_g > self.density_cutoff)
122 self.invweight_g[density_ok] = 1.0 / density_g[density_ok]
124 den_sg, gd = self.calculator.density.get_all_electron_density(
125 self.atoms)
126 assert gd == self.calculator.density.finegd
127 self.den_g = den_sg.sum(axis=0)
129 def get_calculator(self):
130 return self.calculator
132 def get_effective_volume_ratios(self):
133 """Return the list of effective volume to free volume ratios."""
134 self.initialize()
135 with self.timer('HirshfeldPartitioning ratios'):
136 kptband_comm = self.calculator.comms['D']
137 ratios = []
138 for a, atom in enumerate(self.atoms):
139 ratios.append(self.get_effective_volume_ratio(a))
141 ratios = np.array(ratios)
142 kptband_comm.broadcast(ratios, 0)
143 return ratios
145 def get_effective_volume_ratio(self, atom_index):
146 """Effective volume to free volume ratio.
148 After: Tkatchenko and Scheffler PRL 102 (2009) 073005, eq. (7)
149 """
150 finegd = self.calculator.density.finegd
151 denfree_g, gd = self.hdensity.get_density([atom_index])
152 assert gd == finegd
154 # the atoms r^3 grid
155 position = self.atoms[atom_index].position / Bohr
156 r_vg, r2_g = coordinates(finegd, origin=position)
157 r3_g = r2_g * np.sqrt(r2_g)
159 weight_g = denfree_g * self.invweight_g
161 nom = finegd.integrate(r3_g * self.den_g * weight_g)
162 denom = finegd.integrate(r3_g * denfree_g)
164 return nom / denom
166 def get_weight(self, atom_index):
167 denfree_g, gd = self.hdensity.get_density([atom_index])
168 weight_g = denfree_g * self.invweight_g
169 return weight_g
171 def get_charges(self, den_g=None):
172 """Charge on the atom according to the Hirshfeld partitioning
174 Can be applied to any density den_g.
175 """
176 self.initialize()
177 finegd = self.calculator.density.finegd
179 if den_g is None:
180 den_sg, gd = self.calculator.density.get_all_electron_density(
181 self.atoms)
182 den_g = den_sg.sum(axis=0)
183 assert den_g.shape == tuple(finegd.n_c)
185 charges = []
186 for ia, atom in enumerate(self.atoms):
187 weight_g = self.get_weight(ia)
188# charge = atom.number - finegd.integrate(weight_g * den_g)
189 charges.append(atom.number - finegd.integrate(weight_g * den_g))
190 return charges