Coverage for hiphive/self_consistent_phonons/self_consistent_harmonic_model.py: 92%

1import numpy as np 

2from hiphive.force_constant_model import ForceConstantModel 

3from hiphive import StructureContainer 

4from hiphive.utilities import prepare_structures 

5from hiphive.structure_generation import generate_rattled_structures, \ 

6 generate_phonon_rattled_structures 

7from trainstation import Optimizer 

8 

9 

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

11 n_iterations, n_structures, 

12 QM_statistics=False, parameters_start=None, 

13 imag_freq_factor=1.0, fit_kwargs={}): 

14 """ 

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

16 that provides the closest match to the potential energy surface at 

17 a the specified temperature. 

18 

19 Parameters 

20 ---------- 

21 atoms_ideal : ase.Atoms 

22 ideal structure 

23 calc : ASE calculator object 

24 `calculator 

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

26 to be used as reference potential 

27 cs : ClusterSpace 

28 clusterspace onto which to project the reference potential 

29 T : float 

30 temperature in K 

31 alpha : float 

32 stepsize in optimization algorithm 

33 n_iterations : int 

34 number of iterations in poor mans 

35 n_structures : int 

36 number of structures to use when fitting 

37 QM_statistics : bool 

38 if the amplitude of the quantum harmonic oscillator shoule be used 

39 instead of the classical amplitude in the phonon rattler 

40 parameters_start : numpy.ndarray 

41 parameters from which to start the optimization 

42 image_freq_factor: float 

43 If a squared frequency, w2, is negative then it is set to 

44 w2 = imag_freq_factor * np.abs(w2) 

45 fit_kwargs : dict 

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

47 

48 Returns 

49 ------- 

50 list(numpy.ndarray) 

51 sequence of parameter vectors generated while iterating to 

52 self-consistency 

53 """ 

54 

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

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

57 

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

59 raise ValueError('ClusterSpace must be second order') 

60 

61 # initialize things 

62 sc = StructureContainer(cs) 

63 fcm = ForceConstantModel(atoms_ideal, cs) 

64 

65 # generate initial model 

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

67 print('Creating initial model') 

68 rattled_structures = generate_rattled_structures(atoms_ideal, n_structures, 0.03) 

69 rattled_structures = prepare_structures(rattled_structures, atoms_ideal, calc, False) 

70 for structure in rattled_structures: 

71 sc.add_structure(structure) 

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

73 opt.train() 

74 parameters_start = opt.parameters 

75 sc.delete_all_structures() 

76 

77 # run poor mans self consistent 

78 parameters_old = parameters_start.copy() 

79 parameters_traj = [parameters_old] 

80 

81 for i in range(n_iterations): 

82 # generate structures with old model 

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

84 fcm.parameters = parameters_old 

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

86 phonon_rattled_structures = generate_phonon_rattled_structures( 

87 atoms_ideal, fc2, n_structures, T, QM_statistics=QM_statistics, 

88 imag_freq_factor=imag_freq_factor) 

89 phonon_rattled_structures = prepare_structures( 

90 phonon_rattled_structures, atoms_ideal, calc, False) 

91 

92 # fit new model 

93 for structure in phonon_rattled_structures: 

94 sc.add_structure(structure) 

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

96 opt.train() 

97 sc.delete_all_structures() 

98 

99 # update parameters 

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

101 

102 # print iteration summary 

103 disps = [atoms.get_array('displacements') for atoms in phonon_rattled_structures] 

104 disp_ave = np.mean(np.abs(disps)) 

105 disp_max = np.max(np.abs(disps)) 

106 x_new_norm = np.linalg.norm(parameters_new) 

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

108 print(' |x_new| = {:.5f}, |delta x| = {:.8f}, disp_ave = {:.5f}, disp_max = {:.5f}, ' 

109 'rmse = {:.5f}'.format(x_new_norm, delta_x_norm, disp_ave, disp_max, opt.rmse_train)) 

110 parameters_traj.append(parameters_new) 

111 parameters_old = parameters_new 

112 

113 return parameters_traj