Hide keyboard shortcuts

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

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

"""Contains a calculator which given an arbitrary list of clusters and 

associated force constants can calculate the energy and forces of a displaced 

system 

""" 

import math 

from collections import defaultdict as dd 

from typing import List, Tuple 

 

import numpy as np 

from ase import Atoms 

from ase.calculators.calculator import Calculator, all_changes 

from hiphive.force_constants import SortedForceConstants 

from hiphive.utilities import get_displacements 

from .numba_calc import _get_forces 

 

 

class ForceConstantCalculator(Calculator): 

"""This class provides an ASE calculator that can be used in conjunction 

with integrators and optimizers with the `atomic simulation environment 

(ASE) <https://wiki.fysik.dtu.dk/ase/index.html>`_. To initialize an object 

of this class one must provide the ideal atomic configuration along with a 

compatible force constant model. 

 

Parameters 

----------- 

fcs : ForceConstants 

the force constants instance must include the atomic configuration 

max_disp : float 

maximum allowed displacement before calculator raises ValueError 

""" 

 

implemented_properties = ['energy', 'forces'] 

 

def __init__(self, fcs: SortedForceConstants, max_disp: float = 3.0): 

Calculator.__init__(self) 

 

37 ↛ 38line 37 didn't jump to line 38, because the condition on line 37 was never true if not isinstance(fcs, SortedForceConstants): 

raise TypeError('The FC calculator requires sorted FCs.') 

self.atoms_ideal = fcs.supercell.copy() 

 

# Nearest neighbor distance used as maximum displacement allowed, 

# stops exploding MD simulations. 

self.max_allowed_disp = max_disp 

 

self.clusters = dd(list) 

self.force_constants = dd(list) 

self.atom_indices = dd(list) 

self.atom_positions = dd(list) 

self.atom_counts = dd(list) 

self.prefactors = dd(list) 

# The main idea is to precompute the prefactor and multiplicities of 

# belonging to each cluster 

for cluster, fc in fcs.get_fc_dict().items(): 

argsort = np.argsort(cluster) # TODO: is already True? 

cluster = np.array(sorted(cluster)) 

nbody = len(set(cluster)) 

order = len(cluster) 

key = (nbody, order) 

self.clusters[key].append(cluster) 

assert fc.shape == (3,) * order 

self.force_constants[key].append(fc.transpose(argsort)) 

unique = np.unique(cluster, return_index=True, return_counts=True) 

self.atom_indices[key].append(unique[0]) 

self.atom_positions[key].append(unique[1]) 

self.atom_counts[key].append(unique[2]) 

prefac = [-count/np.prod(list(map(math.factorial, unique[2]))) for count in unique[2]] 

self.prefactors[key].append(prefac) 

for d in [self.clusters, 

self.force_constants, 

self.atom_indices, 

self.atom_positions, 

self.atom_counts, 

self.prefactors, 

]: 

for k, v in d.items(): 

d[k] = np.array(v) 

 

def calculate(self, atoms: Atoms = None, 

properties: List[str] = ['energy'], 

system_changes: List[str] = all_changes) -> None: 

Calculator.calculate(self, atoms, properties, system_changes) 

self._check_atoms() 

self._compute_displacements() 

 

85 ↛ exitline 85 didn't return from function 'calculate', because the condition on line 85 was never false if 'forces' in properties or 'energy' in properties: 

E, forces = self.compute_energy_and_forces() 

self.results['forces'] = forces 

self.results['energy'] = E 

 

def _check_atoms(self) -> None: 

"""Checks that the atomic configuration, with which the calculator is 

associated, is compatible with the ideal configuration provided during 

initialization.""" 

94 ↛ 95line 94 didn't jump to line 95, because the condition on line 94 was never true if len(self.atoms) != len(self.atoms_ideal): 

raise ValueError('Length of atoms does not match reference structure') 

96 ↛ 97line 96 didn't jump to line 97, because the condition on line 96 was never true if not all(self.atoms.numbers == self.atoms_ideal.numbers): 

raise ValueError('Atomic numbers do not match reference structure') 

 

def _compute_displacements(self) -> None: 

"""Evaluates the atomic displacements between the current and the ideal 

(reference) configuration.""" 

self.displacements = get_displacements(self.atoms, self.atoms_ideal) 

 

# sanity check that displacements are not too large 

max_disp = np.max(np.linalg.norm(self.displacements, axis=1)) 

106 ↛ 107line 106 didn't jump to line 107, because the condition on line 106 was never true if max_disp > self.max_allowed_disp: 

msg = 'Displacement {:.5f} larger than maximum allowed displacement {:.5f}' 

raise ValueError(msg.format(max_disp, self.max_allowed_disp)) 

 

def compute_energy_and_forces(self) -> Tuple[float, list]: 

"""Computes energy and forces. 

 

Returns 

------- 

float, list(list(float)) 

energy and forces 

""" 

E = np.zeros(1) 

forces = np.zeros((len(self.atoms), 3)) 

 

for key in self.clusters.keys(): 

nbody = key[0] 

order = key[1] 

_get_forces(self.clusters[key], 

self.force_constants[key], 

self.atom_indices[key], 

self.atom_positions[key], 

self.atom_counts[key], 

self.prefactors[key], 

nbody, order, 

forces, E, self.displacements) 

return float(E), forces 

 

def __repr__(self) -> str: 

fc_dict_str = '{{{}: {}, ...}}'.format(self.clusters[0], self.force_constants[0]) 

fcs_str = 'ForceConstants(fc_dict={}, atoms={!r})'.format(fc_dict_str, self.atoms_ideal) 

return 'ForceConstantCalculator({})'.format(fcs_str)