Source code for hiphive.self_consistent_phonons.self_consistent_harmonic_model

import numpy as np
from hiphive.force_constant_model import ForceConstantModel
from hiphive import StructureContainer
from hiphive.utilities import prepare_structures
from hiphive.structure_generation import generate_rattled_structures, \
                                         generate_phonon_rattled_structures
from trainstation import Optimizer


[docs]def self_consistent_harmonic_model(atoms_ideal, calc, cs, T, alpha, n_iterations, n_structures, QM_statistics=False, parameters_start=None, imag_freq_factor=1.0, 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 QM_statistics : bool if the amplitude of the quantum harmonic oscillator shoule be used instead of the classical amplitude in the phonon rattler parameters_start : numpy.ndarray parameters from which to start the optimization image_freq_factor: float If a squared frequency, w2, is negative then it is set to w2 = imag_freq_factor * np.abs(w2) 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 """ if not 0 < alpha <= 1: raise ValueError('alpha must be between 0.0 and 1.0') 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 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, False) 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, QM_statistics=QM_statistics, imag_freq_factor=imag_freq_factor) phonon_rattled_structures = prepare_structures( phonon_rattled_structures, atoms_ideal, calc, False) # 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 # print iteration summary disps = [atoms.get_array('displacements') for atoms in phonon_rattled_structures] disp_ave = np.mean(np.abs(disps)) disp_max = np.max(np.abs(disps)) 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}, disp_ave = {:.5f}, disp_max = {:.5f}, ' 'rmse = {:.5f}'.format(x_new_norm, delta_x_norm, disp_ave, disp_max, opt.rmse_train)) parameters_traj.append(parameters_new) parameters_old = parameters_new return parameters_traj