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

138

139

140

141

142

143

"""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 numpy as np 

import math 

from collections import defaultdict as dd 

 

from hiphive.force_constants import SortedForceConstants 

from ase.calculators.calculator import Calculator, all_changes 

from ase.geometry import find_mic 

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 

""" 

 

implemented_properties = ['energy', 'forces'] 

 

def __init__(self, fcs): 

Calculator.__init__(self) 

 

33 ↛ 34line 33 didn't jump to line 34, because the condition on line 33 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 = 2 * np.min(sorted(np.unique( 

self.atoms_ideal.get_all_distances(mic=True).round(4)))[1]) 

 

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=None, properties=['energy'], 

system_changes=all_changes): 

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

self._check_atoms() 

self._compute_displacements() 

 

82 ↛ exitline 82 didn't return from function 'calculate', because the condition on line 82 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): 

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

associated, is compatible with the ideal configuration provided during 

initialization.""" 

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

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

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

raise ValueError('Atomic numbers does not match reference atoms') 

 

def _compute_displacements(self): 

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

(reference) configuration.""" 

displacements = [] 

for pos, ideal_pos in zip(self.atoms.positions, 

self.atoms_ideal.positions): 

v_ij = np.array([pos - ideal_pos]) 

displacements.append(find_mic(v_ij, self.atoms.cell, 

pbc=True)[0][0]) 

self.displacements = np.array(displacements) 

 

# sanity check that displacements are not too large 

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

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

raise ValueError( 

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

' {:.5f}'.format(max_disp, self.max_allowed_disp)) 

 

def compute_energy_and_forces(self): 

"""Compute 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): 

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)