-
Notifications
You must be signed in to change notification settings - Fork 839
/
ewald.py
775 lines (649 loc) · 27.6 KB
/
ewald.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.
"""
This module provides classes for calculating the ewald sum of a structure.
"""
import bisect
from copy import copy, deepcopy
from datetime import datetime
from math import log, pi, sqrt
from typing import Dict
from warnings import warn
import numpy as np
import scipy.constants as constants
from monty.json import MSONable
from scipy.special import comb, erfc
from pymatgen.core.structure import Structure
__author__ = "Shyue Ping Ong, William Davidson Richard"
__copyright__ = "Copyright 2011, The Materials Project"
__credits__ = "Christopher Fischer"
__version__ = "1.0"
__maintainer__ = "Shyue Ping Ong"
__email__ = "shyuep@gmail.com"
__status__ = "Production"
__date__ = "Aug 1 2012"
class EwaldSummation(MSONable):
"""
Calculates the electrostatic energy of a periodic array of charges using
the Ewald technique.
Ref:
Ewald summation techniques in perspective: a survey
Abdulnour Y. Toukmaji and John A. Board Jr.
DOI: 10.1016/0010-4655(96)00016-1
URL: http://www.ee.duke.edu/~ayt/ewaldpaper/ewaldpaper.html
This matrix can be used to do fast calculations of ewald sums after species
removal.
E = E_recip + E_real + E_point
Atomic units used in the code, then converted to eV.
"""
# Converts unit of q*q/r into eV
CONV_FACT = 1e10 * constants.e / (4 * pi * constants.epsilon_0)
def __init__(
self,
structure,
real_space_cut=None,
recip_space_cut=None,
eta=None,
acc_factor=12.0,
w=1 / sqrt(2),
compute_forces=False,
):
"""
Initializes and calculates the Ewald sum. Default convergence
parameters have been specified, but you can override them if you wish.
Args:
structure (Structure): Input structure that must have proper
Species on all sites, i.e. Element with oxidation state. Use
Structure.add_oxidation_state... for example.
real_space_cut (float): Real space cutoff radius dictating how
many terms are used in the real space sum. Defaults to None,
which means determine automagically using the formula given
in gulp 3.1 documentation.
recip_space_cut (float): Reciprocal space cutoff radius.
Defaults to None, which means determine automagically using
the formula given in gulp 3.1 documentation.
eta (float): The screening parameter. Defaults to None, which means
determine automatically.
acc_factor (float): No. of significant figures each sum is
converged to.
w (float): Weight parameter, w, has been included that represents
the relative computational expense of calculating a term in
real and reciprocal space. Default of 0.7 reproduces result
similar to GULP 4.2. This has little effect on the total
energy, but may influence speed of computation in large
systems. Note that this parameter is used only when the
cutoffs are set to None.
compute_forces (bool): Whether to compute forces. False by
default since it is usually not needed.
"""
self._s = structure
self._charged = abs(structure.charge) > 1e-8
self._vol = structure.volume
self._compute_forces = compute_forces
self._acc_factor = acc_factor
# set screening length
self._eta = eta if eta else (len(structure) * w / (self._vol ** 2)) ** (1 / 3) * pi
self._sqrt_eta = sqrt(self._eta)
# acc factor used to automatically determine the optimal real and
# reciprocal space cutoff radii
self._accf = sqrt(log(10 ** acc_factor))
self._rmax = real_space_cut if real_space_cut else self._accf / self._sqrt_eta
self._gmax = recip_space_cut if recip_space_cut else 2 * self._sqrt_eta * self._accf
# The next few lines pre-compute certain quantities and store them.
# Ewald summation is rather expensive, and these shortcuts are
# necessary to obtain several factors of improvement in speedup.
self._oxi_states = [compute_average_oxidation_state(site) for site in structure]
self._coords = np.array(self._s.cart_coords)
# Define the private attributes to lazy compute reciprocal and real
# space terms.
self._initialized = False
self._recip = None
self._real, self._point = None, None
self._forces = None
# Compute the correction for a charged cell
self._charged_cell_energy = (
-EwaldSummation.CONV_FACT / 2 * np.pi / structure.volume / self._eta * structure.charge ** 2
)
def compute_partial_energy(self, removed_indices):
"""
Gives total ewald energy for certain sites being removed, i.e. zeroed
out.
"""
total_energy_matrix = self.total_energy_matrix.copy()
for i in removed_indices:
total_energy_matrix[i, :] = 0
total_energy_matrix[:, i] = 0
return sum(sum(total_energy_matrix))
def compute_sub_structure(self, sub_structure, tol=1e-3):
"""
Gives total ewald energy for an sub structure in the same
lattice. The sub_structure must be a subset of the original
structure, with possible different charges.
Args:
substructure (Structure): Substructure to compute Ewald sum for.
tol (float): Tolerance for site matching in fractional coordinates.
Returns:
Ewald sum of substructure.
"""
total_energy_matrix = self.total_energy_matrix.copy()
def find_match(site):
for test_site in sub_structure:
frac_diff = abs(np.array(site.frac_coords) - np.array(test_site.frac_coords)) % 1
frac_diff = [abs(a) < tol or abs(a) > 1 - tol for a in frac_diff]
if all(frac_diff):
return test_site
return None
matches = []
for i, site in enumerate(self._s):
matching_site = find_match(site)
if matching_site:
new_charge = compute_average_oxidation_state(matching_site)
old_charge = self._oxi_states[i]
scaling_factor = new_charge / old_charge
matches.append(matching_site)
else:
scaling_factor = 0
total_energy_matrix[i, :] *= scaling_factor
total_energy_matrix[:, i] *= scaling_factor
if len(matches) != len(sub_structure):
output = ["Missing sites."]
for site in sub_structure:
if site not in matches:
output.append("unmatched = {}".format(site))
raise ValueError("\n".join(output))
return sum(sum(total_energy_matrix))
@property
def reciprocal_space_energy(self):
"""
The reciprocal space energy.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
return sum(sum(self._recip))
@property
def reciprocal_space_energy_matrix(self):
"""
The reciprocal space energy matrix. Each matrix element (i, j)
corresponds to the interaction energy between site i and site j in
reciprocal space.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
return self._recip
@property
def real_space_energy(self):
"""
The real space space energy.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
return sum(sum(self._real))
@property
def real_space_energy_matrix(self):
"""
The real space energy matrix. Each matrix element (i, j) corresponds to
the interaction energy between site i and site j in real space.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
return self._real
@property
def point_energy(self):
"""
The point energy.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
return sum(self._point)
@property
def point_energy_matrix(self):
"""
The point space matrix. A diagonal matrix with the point terms for each
site in the diagonal elements.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
return self._point
@property
def total_energy(self):
"""
The total energy.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
return sum(sum(self._recip)) + sum(sum(self._real)) + sum(self._point) + self._charged_cell_energy
@property
def total_energy_matrix(self):
"""
The total energy matrix. Each matrix element (i, j) corresponds to the
total interaction energy between site i and site j.
Note that this does not include the charged-cell energy, which is only important
when the simulation cell is not charge balanced.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
totalenergy = self._recip + self._real
for i in range(len(self._point)):
totalenergy[i, i] += self._point[i]
return totalenergy
@property
def forces(self):
"""
The forces on each site as a Nx3 matrix. Each row corresponds to a
site.
"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
if not self._compute_forces:
raise AttributeError("Forces are available only if compute_forces is True!")
return self._forces
def get_site_energy(self, site_index):
"""Compute the energy for a single site in the structure
Args:
site_index (int): Index of site
ReturnS:
(float) - Energy of that site"""
if not self._initialized:
self._calc_ewald_terms()
self._initialized = True
if self._charged:
warn("Per atom energies for charged structures not supported in EwaldSummation")
return np.sum(self._recip[:, site_index]) + np.sum(self._real[:, site_index]) + self._point[site_index]
def _calc_ewald_terms(self):
"""
Calculates and sets all ewald terms (point, real and reciprocal)
"""
self._recip, recip_forces = self._calc_recip()
self._real, self._point, real_point_forces = self._calc_real_and_point()
if self._compute_forces:
self._forces = recip_forces + real_point_forces
def _calc_recip(self):
"""
Perform the reciprocal space summation. Calculates the quantity
E_recip = 1/(2PiV) sum_{G < Gmax} exp(-(G.G/4/eta))/(G.G) S(G)S(-G)
where
S(G) = sum_{k=1,N} q_k exp(-i G.r_k)
S(G)S(-G) = |S(G)|**2
This method is heavily vectorized to utilize numpy's C backend for
speed.
"""
numsites = self._s.num_sites
prefactor = 2 * pi / self._vol
erecip = np.zeros((numsites, numsites), dtype=np.float_)
forces = np.zeros((numsites, 3), dtype=np.float_)
coords = self._coords
rcp_latt = self._s.lattice.reciprocal_lattice
recip_nn = rcp_latt.get_points_in_sphere([[0, 0, 0]], [0, 0, 0], self._gmax)
frac_coords = [fcoords for (fcoords, dist, i, img) in recip_nn if dist != 0]
gs = rcp_latt.get_cartesian_coords(frac_coords)
g2s = np.sum(gs ** 2, 1)
expvals = np.exp(-g2s / (4 * self._eta))
grs = np.sum(gs[:, None] * coords[None, :], 2)
oxistates = np.array(self._oxi_states)
# create array where q_2[i,j] is qi * qj
qiqj = oxistates[None, :] * oxistates[:, None]
# calculate the structure factor
sreals = np.sum(oxistates[None, :] * np.cos(grs), 1)
simags = np.sum(oxistates[None, :] * np.sin(grs), 1)
for g, g2, gr, expval, sreal, simag in zip(gs, g2s, grs, expvals, sreals, simags):
# Uses the identity sin(x)+cos(x) = 2**0.5 sin(x + pi/4)
m = (gr[None, :] + pi / 4) - gr[:, None]
np.sin(m, m)
m *= expval / g2
erecip += m
if self._compute_forces:
pref = 2 * expval / g2 * oxistates
factor = prefactor * pref * (sreal * np.sin(gr) - simag * np.cos(gr))
forces += factor[:, None] * g[None, :]
forces *= EwaldSummation.CONV_FACT
erecip *= prefactor * EwaldSummation.CONV_FACT * qiqj * 2 ** 0.5
return erecip, forces
def _calc_real_and_point(self):
"""
Determines the self energy -(eta/pi)**(1/2) * sum_{i=1}^{N} q_i**2
"""
fcoords = self._s.frac_coords
forcepf = 2.0 * self._sqrt_eta / sqrt(pi)
coords = self._coords
numsites = self._s.num_sites
ereal = np.empty((numsites, numsites), dtype=np.float_)
forces = np.zeros((numsites, 3), dtype=np.float_)
qs = np.array(self._oxi_states)
epoint = -(qs ** 2) * sqrt(self._eta / pi)
for i in range(numsites):
nfcoords, rij, js, _ = self._s.lattice.get_points_in_sphere(
fcoords, coords[i], self._rmax, zip_results=False
)
# remove the rii term
inds = rij > 1e-8
js = js[inds]
rij = rij[inds]
nfcoords = nfcoords[inds]
qi = qs[i]
qj = qs[js]
erfcval = erfc(self._sqrt_eta * rij)
new_ereals = erfcval * qi * qj / rij
# insert new_ereals
for k in range(numsites):
ereal[k, i] = np.sum(new_ereals[js == k])
if self._compute_forces:
nccoords = self._s.lattice.get_cartesian_coords(nfcoords)
fijpf = qj / rij ** 3 * (erfcval + forcepf * rij * np.exp(-self._eta * rij ** 2))
forces[i] += np.sum(
np.expand_dims(fijpf, 1) * (np.array([coords[i]]) - nccoords) * qi * EwaldSummation.CONV_FACT,
axis=0,
)
ereal *= 0.5 * EwaldSummation.CONV_FACT
epoint *= EwaldSummation.CONV_FACT
return ereal, epoint, forces
@property
def eta(self):
"""
Returns: eta value used in Ewald summation.
"""
return self._eta
def __str__(self):
if self._compute_forces:
output = [
"Real = " + str(self.real_space_energy),
"Reciprocal = " + str(self.reciprocal_space_energy),
"Point = " + str(self.point_energy),
"Total = " + str(self.total_energy),
"Forces:\n" + str(self.forces),
]
else:
output = [
"Real = " + str(self.real_space_energy),
"Reciprocal = " + str(self.reciprocal_space_energy),
"Point = " + str(self.point_energy),
"Total = " + str(self.total_energy),
"Forces were not computed",
]
return "\n".join(output)
def as_dict(self, verbosity: int = 0) -> Dict:
"""
Json-serialization dict representation of EwaldSummation.
Args:
verbosity (int): Verbosity level. Default of 0 only includes the
matrix representation. Set to 1 for more details.
"""
d = {
"@module": self.__class__.__module__,
"@class": self.__class__.__name__,
"structure": self._s.as_dict(),
"compute_forces": self._compute_forces,
"eta": self._eta,
"acc_factor": self._acc_factor,
"real_space_cut": self._rmax,
"recip_space_cut": self._gmax,
"_recip": None if self._recip is None else self._recip.tolist(),
"_real": None if self._real is None else self._real.tolist(),
"_point": None if self._point is None else self._point.tolist(),
"_forces": None if self._forces is None else self._forces.tolist(),
}
return d
@classmethod
def from_dict(cls, d: Dict, fmt: str = None, **kwargs):
"""
Create an EwaldSummation instance from json serialized dictionary.
"""
summation = cls(
structure=Structure.from_dict(d["structure"]),
real_space_cut=d["real_space_cut"],
recip_space_cut=d["recip_space_cut"],
eta=d["eta"],
acc_factor=d["acc_factor"],
compute_forces=d["compute_forces"],
)
# set previously computed private attributes
if d["_recip"] is not None:
summation._recip = np.array(d["_recip"])
summation._real = np.array(d["_real"])
summation._point = np.array(d["_point"])
summation._forces = np.array(d["_forces"])
summation._initialized = True
return summation
class EwaldMinimizer:
"""
This class determines the manipulations that will minimize an ewald matrix,
given a list of possible manipulations. This class does not perform the
manipulations on a structure, but will return the list of manipulations
that should be done on one to produce the minimal structure. It returns the
manipulations for the n lowest energy orderings. This class should be used
to perform fractional species substitution or fractional species removal to
produce a new structure. These manipulations create large numbers of
candidate structures, and this class can be used to pick out those with the
lowest ewald sum.
An alternative (possibly more intuitive) interface to this class is the
order disordered structure transformation.
Author - Will Richards
"""
ALGO_FAST = 0
ALGO_COMPLETE = 1
ALGO_BEST_FIRST = 2
"""
ALGO_TIME_LIMIT: Slowly increases the speed (with the cost of decreasing
accuracy) as the minimizer runs. Attempts to limit the run time to
approximately 30 minutes.
"""
ALGO_TIME_LIMIT = 3
def __init__(self, matrix, m_list, num_to_return=1, algo=ALGO_FAST):
"""
Args:
matrix: A matrix of the ewald sum interaction energies. This is stored
in the class as a diagonally symmetric array and so
self._matrix will not be the same as the input matrix.
m_list: list of manipulations. each item is of the form
(multiplication fraction, number_of_indices, indices, species)
These are sorted such that the first manipulation contains the
most permutations. this is actually evaluated last in the
recursion since I'm using pop.
num_to_return: The minimizer will find the number_returned lowest
energy structures. This is likely to return a number of duplicate
structures so it may be necessary to overestimate and then
remove the duplicates later. (duplicate checking in this
process is extremely expensive)
"""
# Setup and checking of inputs
self._matrix = copy(matrix)
# Make the matrix diagonally symmetric (so matrix[i,:] == matrix[:,j])
for i in range(len(self._matrix)):
for j in range(i, len(self._matrix)):
value = (self._matrix[i, j] + self._matrix[j, i]) / 2
self._matrix[i, j] = value
self._matrix[j, i] = value
# sort the m_list based on number of permutations
self._m_list = sorted(m_list, key=lambda x: comb(len(x[2]), x[1]), reverse=True)
for mlist in self._m_list:
if mlist[0] > 1:
raise ValueError("multiplication fractions must be <= 1")
self._current_minimum = float("inf")
self._num_to_return = num_to_return
self._algo = algo
if algo == EwaldMinimizer.ALGO_COMPLETE:
raise NotImplementedError("Complete algo not yet implemented for " "EwaldMinimizer")
self._output_lists = []
# Tag that the recurse function looks at at each level. If a method
# sets this to true it breaks the recursion and stops the search.
self._finished = False
self._start_time = datetime.utcnow()
self.minimize_matrix()
self._best_m_list = self._output_lists[0][1]
self._minimized_sum = self._output_lists[0][0]
def minimize_matrix(self):
"""
This method finds and returns the permutations that produce the lowest
ewald sum calls recursive function to iterate through permutations
"""
if self._algo == EwaldMinimizer.ALGO_FAST or self._algo == EwaldMinimizer.ALGO_BEST_FIRST:
return self._recurse(self._matrix, self._m_list, set(range(len(self._matrix))))
return None
def add_m_list(self, matrix_sum, m_list):
"""
This adds an m_list to the output_lists and updates the current
minimum if the list is full.
"""
if self._output_lists is None:
self._output_lists = [[matrix_sum, m_list]]
else:
bisect.insort(self._output_lists, [matrix_sum, m_list])
if self._algo == EwaldMinimizer.ALGO_BEST_FIRST and len(self._output_lists) == self._num_to_return:
self._finished = True
if len(self._output_lists) > self._num_to_return:
self._output_lists.pop()
if len(self._output_lists) == self._num_to_return:
self._current_minimum = self._output_lists[-1][0]
def best_case(self, matrix, m_list, indices_left):
"""
Computes a best case given a matrix and manipulation list.
Args:
matrix: the current matrix (with some permutations already
performed)
m_list: [(multiplication fraction, number_of_indices, indices,
species)] describing the manipulation
indices: Set of indices which haven't had a permutation
performed on them.
"""
m_indices = []
fraction_list = []
for m in m_list:
m_indices.extend(m[2])
fraction_list.extend([m[0]] * m[1])
indices = list(indices_left.intersection(m_indices))
interaction_matrix = matrix[indices, :][:, indices]
fractions = np.zeros(len(interaction_matrix)) + 1
fractions[: len(fraction_list)] = fraction_list
fractions = np.sort(fractions)
# Sum associated with each index (disregarding interactions between
# indices)
sums = 2 * np.sum(matrix[indices], axis=1)
sums = np.sort(sums)
# Interaction corrections. Can be reduced to (1-x)(1-y) for x,y in
# fractions each element in a column gets multiplied by (1-x), and then
# the sum of the columns gets multiplied by (1-y) since fractions are
# less than 1, there is no effect of one choice on the other
step1 = np.sort(interaction_matrix) * (1 - fractions)
step2 = np.sort(np.sum(step1, axis=1))
step3 = step2 * (1 - fractions)
interaction_correction = np.sum(step3)
if self._algo == self.ALGO_TIME_LIMIT:
elapsed_time = datetime.utcnow() - self._start_time
speedup_parameter = elapsed_time.total_seconds() / 1800
avg_int = np.sum(interaction_matrix, axis=None)
avg_frac = np.average(np.outer(1 - fractions, 1 - fractions))
average_correction = avg_int * avg_frac
interaction_correction = average_correction * speedup_parameter + interaction_correction * (
1 - speedup_parameter
)
best_case = np.sum(matrix) + np.inner(sums[::-1], fractions - 1) + interaction_correction
return best_case
@classmethod
def get_next_index(cls, matrix, manipulation, indices_left):
"""
Returns an index that should have the most negative effect on the
matrix sum
"""
# pylint: disable=E1126
f = manipulation[0]
indices = list(indices_left.intersection(manipulation[2]))
sums = np.sum(matrix[indices], axis=1)
if f < 1:
next_index = indices[sums.argmax(axis=0)]
else:
next_index = indices[sums.argmin(axis=0)]
return next_index
def _recurse(self, matrix, m_list, indices, output_m_list=[]):
"""
This method recursively finds the minimal permutations using a binary
tree search strategy.
Args:
matrix: The current matrix (with some permutations already
performed).
m_list: The list of permutations still to be performed
indices: Set of indices which haven't had a permutation
performed on them.
"""
# check to see if we've found all the solutions that we need
if self._finished:
return
# if we're done with the current manipulation, pop it off.
while m_list[-1][1] == 0:
m_list = copy(m_list)
m_list.pop()
# if there are no more manipulations left to do check the value
if not m_list:
matrix_sum = np.sum(matrix)
if matrix_sum < self._current_minimum:
self.add_m_list(matrix_sum, output_m_list)
return
# if we wont have enough indices left, return
if m_list[-1][1] > len(indices.intersection(m_list[-1][2])):
return
if len(m_list) == 1 or m_list[-1][1] > 1:
if self.best_case(matrix, m_list, indices) > self._current_minimum:
return
index = self.get_next_index(matrix, m_list[-1], indices)
m_list[-1][2].remove(index)
# Make the matrix and new m_list where we do the manipulation to the
# index that we just got
matrix2 = np.copy(matrix)
m_list2 = deepcopy(m_list)
output_m_list2 = copy(output_m_list)
matrix2[index, :] *= m_list[-1][0]
matrix2[:, index] *= m_list[-1][0]
output_m_list2.append([index, m_list[-1][3]])
indices2 = copy(indices)
indices2.remove(index)
m_list2[-1][1] -= 1
# recurse through both the modified and unmodified matrices
self._recurse(matrix2, m_list2, indices2, output_m_list2)
self._recurse(matrix, m_list, indices, output_m_list)
@property
def best_m_list(self):
"""
Returns: Best m_list found.
"""
return self._best_m_list
@property
def minimized_sum(self):
"""
Returns: Minimized sum
"""
return self._minimized_sum
@property
def output_lists(self):
"""
Returns: output lists.
"""
return self._output_lists
def compute_average_oxidation_state(site):
"""
Calculates the average oxidation state of a site
Args:
site: Site to compute average oxidation state
Returns:
Average oxidation state of site.
"""
try:
avg_oxi = sum([sp.oxi_state * occu for sp, occu in site.species.items() if sp is not None])
return avg_oxi
except AttributeError:
pass
try:
return site.charge
except AttributeError:
raise ValueError(
"Ewald summation can only be performed on structures "
"that are either oxidation state decorated or have "
"site charges."
)