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

import numpy as np 

from hiphive.force_constant_model import ForceConstantModel 

from hiphive import StructureContainer 

from hiphive.fitting import Optimizer 

from hiphive.utilities import prepare_structures 

from hiphive.structure_generation import generate_rattled_structures, \ 

generate_phonon_rattled_structures 

from hiphive.input_output.logging_tools import set_config 

 

 

# Make poor man shut up 

set_config(level=40) 

 

 

def self_consistent_harmonic_model(atoms_ideal, calc, cs, T, alpha, 

n_iterations, n_structures, 

parameters_start=None, fit_kwargs={}): 

""" 

Constructs a set of self-consistent second-order force constants 

that provides the closest match to the potential energy surface at 

a the specified temperature. 

 

Parameters 

---------- 

atoms_ideal : ase.Atoms 

ideal structure 

calc : ASE calculator object 

`calculator 

<https://wiki.fysik.dtu.dk/ase/ase/calculators/calculators.html>`_ 

to be used as reference potential 

cs : ClusterSpace 

clusterspace onto which to project the reference potential 

T : float 

temperature in K 

alpha : float 

stepsize in optimization algorithm 

n_iterations : int 

number of iterations in poor mans 

n_structures : int 

number of structures to use when fitting 

parameters_start : numpy.ndarray 

parameters from which to start the optimization 

fit_kwargs : dict 

kwargs to be used in the fitting process (via Optimizer) 

 

Returns 

------- 

list(numpy.ndarray) 

sequence of parameter vectors generated while iterating to 

self-consistency 

""" 

 

53 ↛ 54line 53 didn't jump to line 54, because the condition on line 53 was never true if not 0 < alpha <= 1: 

raise ValueError('alpha must be between 0.0 and 1.0') 

 

56 ↛ 57line 56 didn't jump to line 57, because the condition on line 56 was never true if max(cs.cutoffs.orders) != 2: 

raise ValueError('ClusterSpace must be second order') 

 

# initialize things 

sc = StructureContainer(cs) 

fcm = ForceConstantModel(atoms_ideal, cs) 

 

# generate initial model 

64 ↛ 76line 64 didn't jump to line 76, because the condition on line 64 was never false if parameters_start is None: 

print('Creating initial model') 

rattled_structures = generate_rattled_structures(atoms_ideal, n_structures, 0.03) 

rattled_structures = prepare_structures(rattled_structures, atoms_ideal, calc) 

for structure in rattled_structures: 

sc.add_structure(structure) 

opt = Optimizer(sc.get_fit_data(), train_size=1.0, **fit_kwargs) 

opt.train() 

parameters_start = opt.parameters 

sc.delete_all_structures() 

 

# run poor mans self consistent 

parameters_old = parameters_start.copy() 

parameters_traj = [parameters_old] 

 

for i in range(n_iterations): 

# generate structures with old model 

print('Iteration {}'.format(i)) 

fcm.parameters = parameters_old 

fc2 = fcm.get_force_constants().get_fc_array(order=2, format='ase') 

phonon_rattled_structures = generate_phonon_rattled_structures( 

atoms_ideal, fc2, n_structures, T) 

phonon_rattled_structures = prepare_structures(phonon_rattled_structures, atoms_ideal, calc) 

 

# fit new model 

for structure in phonon_rattled_structures: 

sc.add_structure(structure) 

opt = Optimizer(sc.get_fit_data(), train_size=1.0, **fit_kwargs) 

opt.train() 

sc.delete_all_structures() 

 

# update parameters 

parameters_new = alpha * opt.parameters + (1-alpha) * parameters_old 

x_new_norm = np.linalg.norm(parameters_new) 

delta_x_norm = np.linalg.norm(parameters_old-parameters_new) 

print(' |x_new| = {:.5f}, |delta x| = {:.8f}, rmse = {:.5f}'.format( 

x_new_norm, delta_x_norm, opt.rmse_train)) 

parameters_traj.append(parameters_new) 

parameters_old = parameters_new 

 

return parameters_traj