Coverage for gpaw/cdft/cdft.py: 71%
595 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1"""Module for constrained DFT
3for review see Chem. Rev., 2012, 112 (1), pp 321-370
4article on GPAW implementation:
5J. Chem. Theory Comput., 2016, 12 (11), pp 5367-5378
6"""
8import copy
9import functools
10from math import pi
12import numpy as np
13from ase.calculators.calculator import Calculator
14from ase.data import atomic_numbers, covalent_radii
15from ase.units import Bohr, Hartree
16from ase.utils import IOContext
17from scipy.optimize import minimize
19from gpaw.external import ExternalPotential
22class CDFT(Calculator):
23 implemented_properties = ['energy', 'forces']
25 def __init__(self,
26 calc,
27 atoms,
28 charge_regions=[],
29 charges=None,
30 spin_regions=[],
31 spins=None,
32 charge_coefs=None,
33 spin_coefs=None,
34 promolecular_constraint=False,
35 txt='-',
36 minimizer_options={'gtol': 0.01},
37 Rc={},
38 mu={'Li': 0.5,
39 'F': 0.7,
40 'O': 0.7,
41 'V': 0.5},
42 method='BFGS',
43 forces='analytical',
44 use_charge_difference=False,
45 compute_forces=True,
46 maxstep=100,
47 tol=1e-3,
48 bounds=None,
49 restart=False,
50 hess=None):
51 """Constrained DFT calculator.
53 calc: GPAW instance
54 DFT calculator object to be constrained.
55 charge_regions: list of list of int
56 Atom indices of atoms in the different charge_regions.
57 spin_regions: list of list of int
58 Atom indices of atoms in the different spin_regions.
59 charges: list of float
60 constrained charges in the different charge_regions.
61 spins: list of float
62 constrained spins in the different charge_regions.
63 Value of 1 sets net magnetisation of one up/alpha electron
64 charge_coefs: list of float
65 Initial values for charge constraint coefficients (eV).
66 spin_coefs: list of float
67 Initial values for spin constraint coefficients (eV).
68 promolecular_constraint: bool
69 Define charge and/or spin constraints from promolecular
70 densities, see: dx.doi.org/10.1021/cr200148b Eq. 29-31
71 If true, user specified charges/spins are overwritten!
72 The atoms (of Atoms object) specifying the charge/spin regions
73 need to contain have correct charge/spin state!
74 (atoms.set_initial_charges([atomic_charges]) and
75 atoms.set_initial_magnetic_moments([moments]))
76 txt: None or str or file descriptor
77 Log file. Default id '-' meaning standard out. Use None for
78 no output.
79 minimizer_options: dict
80 options for scipy optimizers, see:scipy.optimize.minimize
81 method: str
82 One of scipy optimizers, e.g., BFGS, CG
83 forces: str
84 cDFT weight function contribution to forces
85 'fd' for finite difference or 'analytical'
86 difference: bool
87 If True, constrain charge difference between two regions
88 Then charge_regions needs two regions and charges needs
89 only one item which is the charge difference between
90 the two regions, the first beign donor, the second acceptor
92 If False, each region is treated with the corresponding
93 charge constraint
94 compute_forces: bool
95 Should the forces be computed?
96 restart: bool
97 starting from an old calculation from gpw
98 hess: '2-point', '3-point', 'cs'
99 scipy hessian approximation
100 """
102 Calculator.__init__(self)
104 self.calc = calc
105 self.restart = restart
106 self.iocontext = IOContext()
107 self.log = self.iocontext.openfile(txt, calc.world)
108 self.method = method
109 self.forces = forces
110 self.options = minimizer_options
111 self.difference = use_charge_difference
112 self.compute_forces = compute_forces
113 self.Rc = Rc
114 self.mu = mu
115 # set charge constraints and lagrangians
116 self.v_i = np.empty(shape=(0, 0))
118 self.constraints = np.empty(shape=(0, 0))
120 self.charge_regions = charge_regions
121 self.spin_regions = spin_regions
122 self._n_charge_regions = len(self.charge_regions)
123 self._n_spin_regions = len(self.spin_regions)
124 self._n_regions = self._n_charge_regions + self._n_spin_regions
126 self.max_step = maxstep
127 self.tol = tol
128 self.gtol = minimizer_options['gtol']
129 self.bounds = bounds
131 if self.bounds is not None:
132 self.bounds = np.asarray(self.bounds) / Hartree
133 self.hess = hess
135 if self.difference:
136 # difference calculation only for 2 charge regions
137 if self.n_spin_regions != 0 or self._n_charge_regions != 2:
138 raise ValueError('No spin constraints '
139 'for charge difference calculations and'
140 ' only two regions allowed')
141 assert (self.n_charge_regions == 1)
143 if self.n_charge_regions == 0:
144 self.regions = []
146 else:
147 self.charge_i = np.array(charges, dtype=float)
149 if charge_coefs is None: # to Hartree
150 self.v_i = 0.1 * np.sign(self.charge_i)
151 else:
152 self.v_i = np.array(charge_coefs) / Hartree
154 if not self.difference:
155 self.regions = copy.copy(self.charge_regions)
157 # The objective is to constrain the number of electrons (nel)
158 # in a certain region --> convert charge to nel
159 Zn = np.zeros(self.n_charge_regions)
160 for j in range(len(Zn)):
161 for atom in atoms[self.charge_regions[j]]:
162 Zn[j] += atom.number
164 # combined spin and charge constraints
165 self.constraints = Zn - self.charge_i
167 else: # constrain charge between two regions
168 nD = 0. # neutral donor
169 nA = 0. # neutral acceptor
171 for atom in atoms[charge_regions[0]]:
172 nD += atom.number
173 for atom in atoms[charge_regions[1]]:
174 nA += atom.number
176 self.dn_core = nD - nA # difference of core
177 self.constraints = [self.dn_core - charges[0]]
178 self.regions = charge_regions
180 # set spin constraints
181 if self.n_spin_regions != 0 and not self.difference:
182 spin_i = np.array(spins, dtype=float)
183 self.constraints = np.append(self.constraints, spin_i)
185 if spin_coefs is None: # to Hartree
186 v_is = 0.1 * np.sign(spin_i)
187 else:
188 v_is = np.array(spin_coefs) / Hartree
190 self.v_i = np.append(self.v_i, v_is)
191 # self.regions.append(spin_regions)
192 for i in range(self.n_spin_regions):
193 self.regions.append(self.spin_regions[i])
195 # initialise without v_ext
196 atoms.calc = self.calc
197 if not self.restart:
198 atoms.get_potential_energy()
200 assert atoms.calc.wfs.nspins == 2
202 self.cdft_initialised = False
204 self.atoms = atoms
205 self.gd = self.calc.density.finegd
207 if promolecular_constraint:
208 self.constraints = get_promolecular_constraints(
209 calc=self.calc,
210 atoms=self.atoms,
211 charge_regions=self.charge_regions,
212 spin_regions=self.spin_regions,
213 charges=charges,
214 spins=spins)
216 # get number of core electrons at each constrained region
217 # used for pseudo free energy to neglect core contributions
218 # in coupling calculation
219 self.n_core_electrons = np.zeros(len(self.regions))
220 for a in self.atoms:
221 for r in range(len(self.regions[:self.n_charge_regions])):
222 if a.index in self.regions[r] and not self.difference:
223 n_core = a.number - self.calc.wfs.setups[a.index].Nv
224 self.n_core_electrons[r] += n_core
225 elif a.index in self.regions[r] and self.difference:
226 if r == 0:
227 n_core = a.number - self.calc.wfs.setups[a.index].Nv
228 self.n_core_electrons[r] += n_core
229 else:
230 n_core = a.number - self.calc.wfs.setups[a.index].Nv
231 self.n_core_electrons[r] -= n_core
233 w = WeightFunc(self.gd, self.atoms, None, self.Rc, self.mu)
234 self.Rc, self.mu = w.get_Rc_and_mu()
235 # construct cdft potential
236 self.ext = CDFTPotential(regions=self.regions,
237 gd=self.gd,
238 atoms=self.atoms,
239 constraints=self.constraints,
240 n_charge_regions=self.n_charge_regions,
241 difference=self.difference,
242 log=self.log,
243 vi=self.v_i,
244 Rc=self.Rc,
245 mu=self.mu)
247 self.calc.set(external=self.ext)
249 self.w = self.ext.w_ig
251 def __del__(self):
252 self.iocontext.close()
254 def calculate(self, atoms, properties, system_changes):
255 # check we're dealing with same atoms
256 if atoms != self.atoms:
257 self.atoms = atoms
258 if not self.restart:
259 Calculator.calculate(self, self.atoms)
261 Calculator.calculate(self, self.atoms)
263 # update positions and weight functions
264 if 'positions' in system_changes or not self.cdft_initialised:
265 self.ext.set_positions_and_atomic_numbers(
266 self.atoms.positions / Bohr, self.atoms.numbers)
267 self.cdft_initialised = True
269 self.atoms.calc = self.calc
271 p = functools.partial(print, file=self.log)
272 self.iteration = 0
273 self.old_v_i = self.v_i.copy()
275 def f(v_i):
276 # nonlocal iteration
277 # very simple step size control
278 diff = v_i - self.old_v_i
280 if np.any(np.abs(diff) >= self.max_step / Hartree):
281 v_i = self.old_v_i + np.sign(
282 v_i - self.old_v_i) * self.max_step / Hartree
284 self.ext.set_levels(v_i)
285 self.v_i = v_i
286 # cDFT free energy <A|H^KS + V_a w_a|A> = Edft + <A|Vw|A>
287 self.Ecdft = self.atoms.get_potential_energy() # in eV
289 # cDFT corrections to energy
290 self.get_atomic_density_correction()
291 self.Ecdft += self.get_energy_correction() * Hartree
293 # get the cDFT gradient
294 dn_i = []
295 Delta_n = self.get_energy_correction(return_density=True)
297 if self.calc.density.nt_sg is None:
298 self.density.interpolate_pseudo_density()
300 self.nt_ag = self.calc.density.nt_sg[0]
301 self.nt_bg = self.calc.density.nt_sg[1]
302 total_electrons = []
304 if self.n_charge_regions != 0:
305 # total pseudo electron density
306 n_gt = self.nt_ag + self.nt_bg
308 n_electrons = (self.gd.integrate(
309 self.ext.w_ig[0:self.n_charge_regions] * n_gt,
310 global_integral=True))
312 # corrections
313 n_electrons += Delta_n[0:self.n_charge_regions]
315 # constraint
316 diff = n_electrons - self.constraints[0:self.n_charge_regions]
317 total_electrons.append(n_electrons)
318 dn_i.append(diff)
320 if self.n_spin_regions != 0:
321 # difference of pseudo spin densities
322 Dns_gt = (self.nt_ag - self.nt_bg)
323 n_electrons = self.gd.integrate(
324 self.ext.w_ig[self.n_charge_regions:] * Dns_gt,
325 global_integral=True)
326 # corrections
327 n_electrons += Delta_n[self.n_charge_regions:]
328 # constraint
329 diff = n_electrons - self.constraints[self.n_charge_regions:]
330 total_electrons.append(n_electrons)
331 dn_i.append(diff)
333 self.dn_i = np.asarray(dn_i)
334 self.dn_i = self.dn_i.flatten()
335 total_electrons = np.asarray(total_electrons)
336 self.w = self.ext.w_ig
337 # Do not include external potential
338 E_KS = get_ks_energy_wo_external(self.calc)
340 self.Edft = E_KS
341 # pseudo free energy, neglecting core electrons as done for
342 # coupling constant calculation
343 if not self.difference:
345 self.Ecdft = E_KS + np.dot(
346 self.v_i * Hartree,
347 (total_electrons - self.n_core_electrons)[0])
349 if self.iteration == 0:
350 p(f'Optimizer: {self.method}')
351 p(f'Optimizer setups:{self.options}')
353 header = '----------------------------------------\n'
354 header += 'Iteration: {0}\n'
355 header += 'Energy: {1:10.8f} eV\n'
356 header += 'Coefficients [eV]: {2} \n'
357 header += 'Errors [e]: {3} \n'
359 dn_intermediate = self.dn_i.tolist()
360 p(
361 header.format(self.iteration,
362 self.Edft,
363 ''.join(f'{v:4.3f} '
364 for v in self.v_i * Hartree),
365 ''.join(f'{dn:6.4f} '
366 for dn in dn_intermediate),
367 fill='left',
368 align='left'))
370 self.log.flush()
371 self.iteration += 1
373 self.old_v_i = self.v_i.copy()
374 return np.max(np.abs(self.dn_i))
376 m = minimize(f,
377 self.v_i,
378 jac=self.jacobian,
379 hess=self.hess,
380 bounds=self.bounds,
381 tol=self.tol,
382 method=self.method,
383 options=self.options)
385 p(str(m.message) + '\n')
387 self.density = self.calc.density # TS09-vdw needs this
388 self.results['energy'] = self.Edft
390 # print to log
392 p('Final DFT energy : ' + str(self.Edft) + ' eV \n')
393 p('CDFT free energy <A|H+Vw|A> : ' + str(self.Ecdft) + ' eV \n')
395 if self.compute_forces:
396 f = WeightFunc(self.gd,
397 self.atoms,
398 self.regions,
399 self.Rc,
400 self.mu,
401 new=False)
403 f_cdft = f.get_cdft_forces2(dens=self.calc.density,
404 v_i=self.v_i,
405 n_charge_regions=self.n_charge_regions,
406 n_spin_regions=self.n_spin_regions,
407 w_ig=self.w,
408 method=self.forces,
409 difference=self.difference)
411 self.calc.wfs.world.broadcast(f_cdft, 0)
412 self.ext.set_forces(f_cdft)
413 self.results['forces'] = self.atoms.get_forces()
415 def get_weight(self, save=True, name='weight', pad=False):
416 if not pad:
417 w_g = self.w
418 else:
419 c = CDFTPotential(regions=self.regions,
420 gd=self.gd,
421 atoms=self.atoms,
422 constraints=self.constraints,
423 n_charge_regions=self.n_charge_regions,
424 difference=self.difference,
425 vi=self.v_i,
426 Rc=self.Rc,
427 mu=self.mu,
428 log=self.log)
430 w_g = c.initialize_partitioning(self.gd,
431 construct=True,
432 pad=True,
433 global_array=True)
434 if save:
435 w_s = self.gd.collect(w_g, broadcast=True)
437 if self.gd.comm.rank == 0:
438 np.save('coarse_weight', w_s)
439 return w_g
441 def jacobian(self, v_i):
442 if np.all(np.abs(self.dn_i) < self.gtol):
443 # forces scipy opt to converge when gtol is reached
444 return np.zeros(len(self.v_i))
445 else:
446 return -self.dn_i
448 def cdft_free_energy(self):
449 return self.Ecdft
451 def dft_energy(self):
452 return self.Edft
454 def get_lagrangians(self):
455 return self.v_i * Hartree
457 def get_constraints(self):
458 return self.constraints
460 def get_grid(self):
461 return self.gd
463 def get_all_electron_density(self, gridrefinement=2, spin=None):
464 return self.calc.get_all_electron_density(
465 gridrefinement=gridrefinement, spin=spin)
467 def get_atomic_density_correction(self, return_els=False):
468 # eq. 20 of the paper
469 self.dn_s = np.zeros((2, len(self.atoms)))
471 for i in [0, 1]:
472 for a, D_sp in self.calc.density.D_asp.items():
473 self.dn_s[i, a] += np.sqrt(4 * np.pi) * (
474 np.dot(D_sp[i],
475 self.calc.wfs.setups[a].Delta_pL)[0] +
476 self.calc.wfs.setups[a].Delta0 / 2)
478 self.gd.comm.sum(self.dn_s)
479 for a in range(len(self.atoms)):
480 self.dn_s[:, a] += self.atoms[a].number / 2.
481 if return_els:
482 return self.dn_s
484 def get_energy_correction(self, return_density=False):
485 # Delta n^a part of eq 21
487 # for each region
488 n_a = np.zeros(len(self.regions))
490 # int w_i Dn_i for both spins
491 # in spin constraints w_ib = -w_ia
492 # inside augmentation spheres w_i = 1
494 for c in range(len(self.regions)):
495 # sum all atoms in a region
496 n_sa = self.dn_s[0, self.regions[c]].sum()
497 n_sb = self.dn_s[1, self.regions[c]].sum()
498 # total density correction
499 n_a[c] = n_sa + n_sb
501 for s in range(self.n_spin_regions):
502 n_sa = self.dn_s[0, self.regions[self.n_charge_regions + s]].sum()
503 n_sb = self.dn_s[1, self.regions[self.n_charge_regions + s]].sum()
504 n_a[self.n_charge_regions + s] = n_sa - n_sb
506 if return_density:
507 if not self.difference:
508 # Delta n^a, eq 20
509 return n_a
510 else:
511 # the difference of corrections
512 return [n_a[0] - n_a[1]]
514 else:
515 if not self.difference:
516 return (np.dot(self.v_i, n_a))
517 else:
518 # negative for difference acceptor
519 vi_temp = np.array([self.v_i[0], -self.v_i[0]])
520 return (np.dot(vi_temp, n_a))
522 def get_number_of_electrons_on_atoms(self):
523 """Return the number of electrons with each gaussian."""
525 nelectrons = []
526 ae_dens_correction = self.get_atomic_density_correction(
527 return_els=True)
528 if self.calc.density.nt_sg is None:
529 self.calc.density.interpolate_pseudo_density()
531 dens = self.calc.density.nt_sg[0] + self.calc.density.nt_sg[1]
533 for atom in self.atoms:
534 # weight function with one atom
535 f = WeightFunc(self.gd,
536 self.atoms, [atom.index],
537 self.Rc,
538 self.mu,
539 new=False)
541 w = f.construct_weight_function()
542 n_el = (self.gd.integrate(w * dens, global_integral=True))
543 # corrections
544 n_el += (ae_dens_correction[:, atom.index]).sum()
545 nelectrons.append(n_el)
547 return nelectrons
549 def write(self, name, mode=None):
550 self.calc.write(name, mode=mode)
552 def save_parameters(self, name='initial', save_weight=True, save_wfs=True):
553 if self.gd.comm.rank == 0:
554 file = open(name + '.txt', 'w')
555 file.write('NA = %f ,\n' % (self.constraints))
556 file.write('FA = %f , \n' % (self.Ecdft))
557 file.write('EA = %f , \n' % (self.Edft))
558 file.write('Va = %f , \n' % (self.v_i * Hartree))
559 file.write('N_charge_regions_A = %d ,\n' % self.n_charge_regions)
560 file.close()
562 def get_weight_function_on_coarse_grid(self,
563 regions,
564 gd,
565 atoms,
566 constraints,
567 n_charge_regions,
568 difference,
569 save=True):
571 gd = self.calc.density.gd
572 c = CDFTPotential(regions=self.regions,
573 gd=gd,
574 atoms=self.atoms,
575 constraints=self.constraints,
576 n_charge_regions=self.n_charge_regions,
577 difference=self.difference,
578 vi=self.v_i,
579 log=self.log)
581 w_G = c.initialize_partitioning(gd, construct=True)
583 if save:
584 w_s = gd.collect(w_G, broadcast=True)
585 if gd.comm.rank == 0:
586 np.save('coarse_weight', w_s)
587 return w_G
589 def get_coarse_grid(self, save=True):
591 if save:
592 w_s = self.calc.density.gd.collect(self.gd, broadcast=True)
593 if self.calc.density.gd.comm.rank == 0:
594 np.save('coarse_grid', w_s)
596 return self.calc.density.gd
598 @property
599 def n_spin_regions(self):
600 self._n_spin_regions = len(self.spin_regions)
601 return self._n_spin_regions
603 @property
604 def n_charge_regions(self):
605 if not self.difference:
606 self._n_charge_regions = len(self.charge_regions)
607 else:
608 self._n_charge_regions = 1
609 return self._n_charge_regions
611 @property
612 def n_regions(self):
613 if not self.difference:
614 self._n_regions = self._n_charge_regions + self._n_spin_regions
615 else:
616 self._n_regions = 2
617 return self._n_regions
620def gaussians(gd, positions, numbers):
621 r_Rv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
622 radii = covalent_radii[numbers]
623 cutoffs = radii + 3.0
624 sigmas = radii * min(covalent_radii) + 0.5
625 result_R = gd.zeros()
626 for pos, Z, rc, sigma in zip(positions, numbers, cutoffs, sigmas):
627 d2_R = ((r_Rv - pos)**2).sum(3)
628 a_R = Z / (sigma**3 * (2 * np.pi)**1.5) * np.exp(-d2_R /
629 (2 * sigma**2))
630 a_R[d2_R > rc] = 0.0
631 result_R += a_R
632 return result_R
635class CDFTPotential(ExternalPotential):
636 def __init__(self,
637 regions,
638 gd,
639 atoms,
640 constraints,
641 n_charge_regions,
642 difference,
643 vi,
644 log,
645 Rc={},
646 mu={}):
648 self.indices_i = regions
649 self.gd = gd
650 self.iocontext = IOContext()
651 self.log = log
652 self.atoms = atoms
653 self.pos_av = None
654 self.Z_a = None
655 self.w_ig = None
656 self.n_charge_regions = n_charge_regions
657 self.constraints = constraints
658 self.difference = difference
659 self.v_i = vi
660 self.Rc = Rc
661 self.mu = mu
662 self.name = 'CDFTPotential'
664 def __str__(self):
665 self.name = 'CDFTPotential'
666 return 'CDFTPotential'
668 def get_name(self):
669 return self.name
671 def update_ae_density(self, ae_dens):
672 self.ae_dens = ae_dens
674 def get_atoms(self):
675 return self.atoms
677 def get_ae_density(self):
678 return self.ae_dens
680 def get_vi(self):
681 return self.v_i
683 def get_constraints(self):
684 return self.constraints
686 def set_levels(self, v_i):
687 self.v_i = np.array(v_i, dtype=float)
688 self.vext_g = None
689 return self.v_i
691 def set_forces(self, cdft_forces):
692 self.cdft_forces = cdft_forces
694 def get_cdft_forces(self):
695 return self.cdft_forces
697 def spin_polarized_potential(self):
698 return len(self.constraints) != self.n_charge_regions
700 def get_w(self):
701 return self.w_ig
703 def set_positions_and_atomic_numbers(self, pos_av, Z_a):
704 self.pos_av = pos_av
705 self.Z_a = Z_a
706 self.w_ig = None
707 self.vext_g = None
709 def initialize_partitioning(self,
710 gd,
711 construct=False,
712 pad=False,
713 global_array=False):
714 w = get_all_weight_functions(self.atoms, gd, self.indices_i,
715 self.difference, self.Rc, self.mu)
716 if construct:
717 return np.array(w)
719 self.w_ig = np.array(w)
720 p = functools.partial(print, file=self.log)
721 p('Number of charge constrained regions: {n}'.format(
722 n=self.n_charge_regions))
723 if not self.difference:
724 p('Number of spin constrained regions: {n}'.format(
725 n=len(self.indices_i) - self.n_charge_regions))
726 else:
727 p('Number of spin constrained regions: 0')
728 p(f'Charge difference: {self.difference}')
729 p('Parameters')
730 p('Atom Width[A] Rc[A]')
731 for a in self.mu:
732 p(' {atom} {width:.3f} {Rc:.3f}'.format(
733 atom=a, width=self.mu[a] * Bohr, Rc=self.Rc[a] * Bohr))
734 print(file=self.log)
735 self.log.flush()
737 def calculate_potential(self, gd):
738 # return v_ext^{\sigma} = sum_i V_i*w_i^{\sigma}
739 if self.w_ig is None:
740 self.initialize_partitioning(self.gd)
741 pot = np.empty(self.w_ig.shape)
743 for i in range(len(self.v_i)):
744 pot[i] = self.v_i[i] * self.w_ig[i]
746 # first alpha spin
747 vext_sga = np.sum(np.asarray(pot), axis=0)
748 # then beta
749 vext_sgb = np.asarray(pot)
750 # spin constraints with beta spins
751 vext_sgb[self.n_charge_regions:] *= -1.
752 vext_sgb = np.sum(vext_sgb, axis=0)
753 vext_sg = np.array([vext_sga, vext_sgb])
754 # spin-dependent cdft potential
755 self.vext_g = vext_sg
757 def write(self, writer):
758 writer.write(vext='CDFTPotential')
760 def read(self, reader):
761 pass
763 def todict(self):
764 return {
765 'name': 'CDFTPotential',
766 'regions': self.indices_i,
767 'constraints': self.v_i * Hartree,
768 'n_charge_regions': self.n_charge_regions,
769 'difference': self.difference
770 }
772 def get_atomic_hamiltonians(self, setups, atom):
773 h_cdft_a = np.zeros(setups.shape)
774 h_cdft_b = np.zeros(setups.shape)
775 v_i = self.v_i
776 if self.difference:
777 v_i = [v_i, -v_i]
779 for i in range(len(self.indices_i)):
780 if atom in self.indices_i[i]:
781 h_cdft_a += v_i[i] * 2. * np.sqrt(np.pi) * setups
782 h_cdft_b += v_i[i] * 2. * np.sqrt(np.pi) * setups
783 if i >= self.n_charge_regions and not self.difference:
784 h_cdft_b *= -1.
786 return h_cdft_a, h_cdft_b
788 def get_cdft_external_energy(self, dens, nspins, vext_g, vt_g, vbar_g,
789 vt_sg):
790 # cDFT works with all-electron (spin)density
792 # First the potential
793 # spin dependent external potential, array of Vi*wi
794 # vext_g = self.calculate_potential(self.gd).copy()
795 vt_g += vext_g[0]
796 vext_g[1:nspins] = vext_g[1] + vbar_g
797 vt_sg[nspins:] = 0.0
799 # then energy
800 w = self.get_w() # weight functions
801 Vi = self.get_vi()
802 constraints = self.get_constraints()
804 # pseudo electron density on fine grid
805 if dens.nt_sg is None:
806 dens.interpolate_pseudo_density()
808 diff = []
810 if self.n_charge_regions != 0:
811 # pseudo density
812 nt_g = dens.nt_sg[0] + dens.nt_sg[1]
813 charge_diff = (
814 self.gd.integrate(w[0:self.n_charge_regions],
815 nt_g, global_integral=True) -
816 constraints[0:self.n_charge_regions])
817 diff.append(charge_diff)
819 # constrained spins
820 if len(constraints) - self.n_charge_regions != 0:
821 Delta_nt_g = dens.nt_sg[0] - dens.nt_sg[
822 1] # pseudo spin difference density
823 spin_diff = self.gd.integrate(
824 w[self.n_charge_regions:], Delta_nt_g,
825 global_integral=True) - constraints[self.n_charge_regions:]
826 diff.append(spin_diff)
828 diff = np.asarray(diff)
829 # number of domains
830 size = self.gd.comm.size
831 e = np.multiply(Vi, diff[:]).sum()
832 e /= size
834 return e
837class WeightFunc:
838 """ Class which builds a weight function around atoms or molecules
839 given the atom index - using normalized Gaussian with cut-off!
841 The weight function is constructed on the coarse or fine grid and
842 can be used to do charge constraint DFT.
844 """
845 def __init__(self, gd, atoms, indices, Rc={}, mu={}, new=True):
846 """ Given a grid-descriptor, atoms object and an index list
847 construct a weight function defined by:
848 n_i(r-R_i)
849 w(r) = ---------------
850 sum_a n_a(r-R_a)
852 where a runs over all atoms, and i can index
853 an atom or a list of atoms comprising a molecule, etc.
855 The n_i are construced with atom centered gaussians
856 using a pre-defined cut-off Rc_i.
858 """
859 self.gd = gd
860 self.atoms = atoms
861 self.indices_i = indices # Indices of constrained charge_regions
863 # Weight function parameters in Bohr
864 # Cutoffs
865 if new:
866 new_Rc = {}
867 for a in self.atoms:
868 if a.symbol in Rc:
869 new_Rc[a.symbol] = Rc[a.symbol] / Bohr
870 else:
871 element_number = atomic_numbers[a.symbol]
872 cr = covalent_radii[element_number]
873 # Rc to roughly between 3. and 5.
874 new_Rc[a.symbol] = (cr + 2.5) / Bohr
876 self.Rc = new_Rc
877 else:
878 self.Rc = Rc
880 # Construct mu (width) dict
881 # mu only sets the width and height so it's in angstrom
882 if new:
883 new_mu = {}
884 for a in self.atoms:
885 if a.symbol in mu:
886 new_mu[a.symbol] = mu[a.symbol] / Bohr
887 else:
888 element_number = atomic_numbers[a.symbol]
889 cr = covalent_radii[element_number]
890 # mu to be roughly between 0.5 and 1.0 AA
891 cr = (cr * min(covalent_radii) + 0.5)
892 new_mu[a.symbol] = cr / Bohr
894 # "Larger" atoms may need a bit more width
895 self.mu = new_mu
896 else:
897 self.mu = mu
899 def get_Rc_and_mu(self):
900 return self.Rc, self.mu
902 def normalized_gaussian(self, dis, mu, Rc):
903 # Given mu - width, and Rc
904 # a normalized gaussian is constructed
905 # around some atom. This is
906 # placed on the gd (grid) - and truncated
907 # at a given cut-off value Rc. dis
908 # are the distances from atom to grid points.
909 """ Normalized gaussian is:
910 1
911 g(r) = --------------- e^{-(r-Ra)^2/(2mu^2)}
912 mu^3*(2pi)^1.5
914 for |r-Ra| <= Rc, 0 elsewhere
916 """
917 # Check function
918 check = abs(dis) <= Rc
920 # Make gaussian 3D Guassian
921 gauss = 1 / (mu * (2 * pi)**(1 / 2)) * np.exp(-dis**2 / (2.0 * mu**2))
923 # apply cut-off and return
924 return np.array(gauss * check)
926 def get_distance_vectors(self, pos, distance=True):
927 xyz = self.gd.get_grid_point_distance_vectors(pos)
928 if distance:
929 # returns vector norm
930 return np.sqrt((xyz**2).sum(axis=0))
931 else:
932 # gives raw vector
933 return xyz
935 def construct_total_density(self, atoms):
936 # Add to empty grid
937 dens = self.gd.zeros()
939 for atom in atoms:
940 charge = atom.number
941 symbol = atom.symbol
942 pos = atom.position / Bohr
944 dis = self.get_distance_vectors(pos)
946 dens += charge * self.normalized_gaussian(dis,
947 self.mu[symbol],
948 self.Rc[symbol])
949 return dens
951 def construct_weight_function(self):
952 # Grab atomic / molecular density
953 dens_n = self.construct_total_density(self.atoms[self.indices_i])
954 # Grab total density
955 dens = self.construct_total_density(self.atoms)
956 # Check zero elements
957 check = dens == 0
958 # Add value to zeros ...
959 dens += check * 1.0
960 # make weight function
961 return (dens_n / dens)
963 def get_cdft_forces2(self, dens, v_i, n_charge_regions, n_spin_regions,
964 w_ig, method, difference):
965 """ Calculate cDFT force as a sum
966 dF/dRi = Fi(inside) + Fs(surf)
967 due to cutoff (Rc) in gauss
968 / dw(r)
969 Fi_a = -V | ------ n(r) dr
970 / dR_a
971 dw(r)
972 ---- = sum of Gaussian functions...
973 dR_a
975 this is computed in get_dw_dRa
976 dens = density
977 Vc = cDFT constraint value
978 method = 'fd' or 'analytical' for
979 finite difference or analytical
980 dw/dR
981 """
983 cdft_forces = np.zeros((len(self.atoms), 3))
985 if dens.nt_sg is None:
986 dens.interpolate_pseudo_density()
988 rho_kd = self.construct_total_density(self.atoms) # sum_k n_k
989 # Check zero elements
990 check = rho_kd == 0.
991 # Add value to zeros for denominator...
992 rho_kd += check * 1.0
994 for a, atom in enumerate(self.atoms):
995 wn_sg = self.gd.zeros()
996 prefactor = self.get_derivative_prefactor(n_charge_regions,
997 n_spin_regions, w_ig,
998 v_i,
999 difference,
1000 atom, rho_kd)
1002 # make extended array
1003 for c in range(n_charge_regions):
1004 wn_sg += (dens.nt_sg[0] + dens.nt_sg[1]) * prefactor[0]
1006 for s in range(n_spin_regions):
1007 wn_sg += (dens.nt_sg[0] - dens.nt_sg[1]) * prefactor[1]
1009 for i in [0, 1, 2]:
1010 if method == 'analytical':
1011 dG_dRav = self.get_analytical_derivates(atom,
1012 direction=i)
1013 elif method == 'fd':
1014 dG_dRav = self.get_fd_derivatives(atom,
1015 direction=i)
1017 cdft_forces[a][i] -= self.gd.integrate(wn_sg * dG_dRav,
1018 global_integral=True)
1020 return cdft_forces
1022 def get_fd_derivatives(self, atom, direction, dx=1.e-4):
1024 dirs = [[dx, 0., 0.], [0., dx, 0.], [0., 0., dx]]
1025 charge = atom.number
1026 symbol = atom.symbol
1027 mu = self.mu[symbol]
1028 Rc = self.Rc[symbol]
1030 a_posx = atom.position / Bohr + dirs[direction]
1031 a_dis = self.get_distance_vectors(a_posx)
1032 Ga_posx = charge * self.normalized_gaussian(a_dis, mu, Rc)
1033 # move to -dx
1034 a_negx = atom.position / Bohr - dirs[direction]
1035 a_dis = self.get_distance_vectors(a_negx)
1036 Ga_negx = charge * self.normalized_gaussian(a_dis, mu, Rc)
1037 # dG/dx
1038 dGav = (Ga_posx - Ga_negx) / (2 * dx)
1040 return dGav
1042 def get_derivative_prefactor(self, n_charge_regions, n_spin_regions, w_ig,
1043 v_i, difference, atom, rho_kd):
1044 """Computes the dw/dRa array needed for derivatives/forces
1045 eq 31
1046 needed for lfc-derivative/integrals
1047 """
1048 wc = self.gd.zeros()
1049 ws = self.gd.zeros()
1051 for i in range(n_charge_regions):
1052 # build V_i [sum_k rho_k + sum_{j in i}rho_i]
1053 wi = -w_ig[i]
1055 if not difference:
1056 if atom.index in self.indices_i[i]:
1057 wi += 1.
1058 else:
1059 if atom.index in self.indices_i[0]:
1060 wi += 1.
1061 elif atom.index in self.indices_i[1]:
1062 wi -= 1.
1064 wi *= v_i[i]
1065 wc += wi / rho_kd
1067 for i in range(n_spin_regions):
1068 # build V_i [sum_k rho_k + sum_{j in i}rho_i]
1069 wi = -w_ig[n_charge_regions + i]
1070 if atom.index in self.indices_i[n_charge_regions + i]:
1071 wi += 1.
1072 wi *= v_i[n_charge_regions + i]
1073 ws += wi / rho_kd
1075 return [wc, ws]
1077 def get_analytical_derivates(self, atom, direction):
1078 # equations 32,33,34
1080 a_pos = atom.position / Bohr
1081 a_symbol = atom.symbol
1082 a_charge = atom.number
1084 a_dis = self.get_distance_vectors(a_pos)
1085 rRa = -self.get_distance_vectors(a_pos, distance=False)
1086 dist_rRa = self.get_distance_vectors(a_pos, distance=True)
1087 check = dist_rRa == 0
1089 # Add value to zeros ...
1090 dist_rRa += check * 1.0
1091 # eq 33
1092 drRa_di = rRa[direction] / dist_rRa
1094 # Gaussian derivative eq 34
1096 G_a = a_charge * self.normalized_gaussian(a_dis,
1097 self.mu[a_symbol],
1098 self.Rc[a_symbol])
1100 # within cutoff or at surface ? --> heaviside
1101 # inside
1102 check_i = abs(a_dis) <= self.Rc[a_symbol]
1103 rRc = check_i * a_dis
1104 dGa_drRa = -rRc * G_a / (
1105 self.mu[a_symbol])**2 # (\Theta * (r-R_a) n_A) / \sigma^2
1107 # at surface
1108 surf = abs(abs(a_dis) - self.Rc[a_symbol])
1109 check_s = surf <= max(self.gd.get_grid_spacings())
1110 dGa_drRa += check_s * G_a # \ sigma_{A\in i} n_A
1112 # eq 32
1113 return dGa_drRa * drRa_di
1116def get_ks_energy_wo_external(calc):
1117 return (calc.hamiltonian.e_kinetic + calc.hamiltonian.e_coulomb +
1118 calc.hamiltonian.e_zero + calc.hamiltonian.e_xc -
1119 calc.hamiltonian.e_entropy) * Hartree
1122def get_all_weight_functions(atoms,
1123 gd,
1124 indices_i,
1125 difference,
1126 Rc,
1127 mu,
1128 new=False,
1129 return_Rc_mu=False):
1130 w = []
1131 for i in range(len(indices_i)):
1132 wf = WeightFunc(gd, atoms, indices_i[i], Rc, mu, new)
1133 weig = wf.construct_weight_function()
1135 if not difference:
1136 w.append(weig)
1137 else: # for charge difference constraint
1138 if i == 0:
1139 w.append(weig)
1140 else:
1141 w[0] -= weig # negative for acceptor
1142 if return_Rc_mu:
1143 Rc, mu = wf.get_Rc_and_mu()
1144 return w, Rc, mu
1145 else:
1146 return w
1149def get_promolecular_constraints(calc_a,
1150 atoms_a,
1151 calc_b,
1152 atoms_b,
1153 charge_constraint=True,
1154 spin_constraint=False,
1155 restart=False,
1156 Rc={},
1157 mu={}):
1158 """
1159 - calc_a is for the region you're interested in. Its charge
1160 should correspond to the "free charge" of the promolecule
1161 - atoms_a: atoms object for the region of interest
1162 - calc_b is for the other region
1163 - atoms_b: atoms object for the other region
1164 - restart: bool. If true, atoms_a have used calc_a and atoms_b calc_b
1165 """
1167 constraints = []
1168 atoms = atoms_a + atoms_b
1170 if not restart:
1171 atoms_a.calc = calc_a
1172 atoms_b.calc = calc_b
1173 atoms_a.get_potential_energy()
1174 atoms_b.get_potential_energy()
1176 gd = calc_a.density.finegd
1178 if charge_constraint:
1179 # can charge and spin constraints be treated at the same time?
1180 d_a = calc_a.get_all_electron_density(gridrefinement=2, pad=False)
1181 d_b = calc_b.get_all_electron_density(gridrefinement=2, pad=False)
1183 charge_region = [atom.index for atom in atoms_a]
1184 weight = WeightFunc(gd=gd,
1185 atoms=atoms,
1186 indices=charge_region,
1187 Rc=Rc,
1188 mu=mu)
1189 w = weight.construct_weight_function()
1191 dv = atoms.get_volume() / calc_a.get_number_of_grid_points().prod()
1192 gridrefinement = 2.
1193 Nel = (w * (d_a + d_b)).sum() * dv / gridrefinement**3
1195 Zn = 0
1196 for atom in atoms_a:
1197 Zn += atom.number
1199 constraints.append(Zn - Nel)
1201 if spin_constraint:
1202 # can charge and spin constraints be treated at the same time?
1203 d_a = calc_a.get_all_electron_density(gridrefinement=2,
1204 pad=False,
1205 spin=0)
1206 d_b = calc_b.get_all_electron_density(gridrefinement=2,
1207 pad=False,
1208 spin=1)
1210 charge_region = [atom.index for atom in atoms_a]
1211 weight = WeightFunc(gd=gd,
1212 atoms=atoms,
1213 indices=charge_region,
1214 Rc=Rc,
1215 mu=mu)
1216 w = weight.construct_weight_function()
1218 dv = atoms.get_volume() / calc_a.get_number_of_grid_points().prod()
1219 gridrefinement = 2.
1220 Ns = (w * (d_a - d_b)).sum() * dv / gridrefinement**3
1222 constraints.append(Ns)
1224 return constraints