1# Copyright (C) 2018 Atsushi Togo 2# All rights reserved. 3# 4# This file is part of phonopy. 5# 6# Redistribution and use in source and binary forms, with or without 7# modification, are permitted provided that the following conditions 8# are met: 9# 10# * Redistributions of source code must retain the above copyright 11# notice, this list of conditions and the following disclaimer. 12# 13# * Redistributions in binary form must reproduce the above copyright 14# notice, this list of conditions and the following disclaimer in 15# the documentation and/or other materials provided with the 16# distribution. 17# 18# * Neither the name of the phonopy project nor the names of its 19# contributors may be used to endorse or promote products derived 20# from this software without specific prior written permission. 21# 22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 33# POSSIBILITY OF SUCH DAMAGE. 34 35import numpy as np 36from ase.atoms import Atoms 37from hiphive import ClusterSpace, StructureContainer, ForceConstantPotential 38from hiphive.fitting import Optimizer 39from hiphive.cutoffs import estimate_maximum_cutoff 40from hiphive.input_output.logging_tools import set_config 41set_config(level=30) 42 43 44def phonopy_atoms_to_ase(atoms_phonopy): 45 ase_atoms = Atoms(cell=atoms_phonopy.cell, 46 scaled_positions=atoms_phonopy.scaled_positions, 47 numbers=atoms_phonopy.numbers, 48 pbc=True) 49 return ase_atoms 50 51 52def get_fc2(supercell, 53 primitive, 54 displacements, 55 forces, 56 symprec, 57 atom_list=None, 58 options=None, 59 log_level=0): 60 if log_level: 61 msg = [ 62 "-------------------------------" 63 " hiPhive start " 64 "------------------------------", 65 "hiPhive is a non-trivial force constants calculator. " 66 "Please cite the paper:", 67 "\"The Hiphive Package for the Extraction of High‐Order Force " 68 "Constants", 69 " by Machine Learning\"", 70 "by Fredrik Eriksson, Erik Fransson, and Paul Erhart,", 71 "Advanced Theory and Simulations, DOI:10.1002/adts.201800184 " 72 "(2019)", 73 ""] 74 print("\n".join(msg)) 75 76 fc2 = run_hiphive(supercell=supercell, 77 primitive=primitive, 78 displacements=displacements, 79 forces=forces, 80 options=options, 81 symprec=symprec, 82 log_level=log_level) 83 84 p2s_map = primitive.p2s_map 85 is_compact_fc = (atom_list is not None and 86 (atom_list == p2s_map).all()) 87 if is_compact_fc: 88 fc2 = np.array(fc2[p2s_map], dtype='double', order='C') 89 elif atom_list is not None: 90 fc2 = np.array(fc2[atom_list], dtype='double', order='C') 91 92 if log_level: 93 print("--------------------------------" 94 " hiPhive end " 95 "-------------------------------") 96 97 return fc2 98 99 100def run_hiphive(supercell, 101 primitive, 102 displacements, 103 forces, 104 options, 105 symprec, 106 log_level): 107 """Run hiphive 108 109 supercell : Supercell 110 Perfect supercell. 111 primitive : Primitive 112 Primitive cell. 113 displacements : ndarray 114 Displacements of atoms in supercell. 115 shape=(supercells, natom, 3) 116 forces : ndarray 117 Forces on atoms in supercell. 118 shape=(supercells, natom, 3) 119 options : str 120 Force constants calculation options. 121 log_level : int 122 Log control. 0: quiet, 1: normal, 2: verbose 3: debug 123 124 """ 125 126 ase_supercell = phonopy_atoms_to_ase(supercell) 127 ase_prim = phonopy_atoms_to_ase(primitive) 128 129 # setup training structures 130 structures = [] 131 for d, f in zip(displacements, forces): 132 structure = ase_supercell.copy() 133 structure.new_array('displacements', d) 134 structure.new_array('forces', f) 135 structures.append(structure) 136 137 # parse options 138 if options is None: 139 options_dict = {} 140 else: 141 options_dict = _decode_options(options) 142 143 # select cutoff 144 max_cutoff = estimate_maximum_cutoff(ase_supercell) - 1e-5 145 if 'cutoff' in options_dict: 146 cutoff = options_dict['cutoff'] 147 if cutoff > max_cutoff: 148 raise ValueError('Cutoff {:.4f} is larger than maximum allowed ' 149 'cutoff, {:.4f}, for the given supercell.' 150 '\nDecrease cutoff or provide larger supercells.' 151 .format(cutoff, max_cutoff)) 152 else: 153 cutoff = max_cutoff 154 155 # setup ClusterSpace 156 cutoffs = [cutoff] 157 cs = ClusterSpace(ase_prim, cutoffs, symprec=symprec) 158 cs.print_orbits() 159 160 sc = StructureContainer(cs) 161 for structure in structures: 162 sc.add_structure(structure) 163 n_rows, n_cols = sc.data_shape 164 if n_rows < n_cols: 165 raise ValueError('Fitting problem is under-determined.' 166 '\nProvide more structures or decrease cutoff.') 167 168 # Estimate error 169 opt = Optimizer(sc.get_fit_data(), train_size=0.75) 170 opt.train() 171 print(opt) 172 print('RMSE train : {:.4f}'.format(opt.rmse_train)) 173 print('RMSE test : {:.4f}'.format(opt.rmse_test)) 174 175 # Final train 176 opt = Optimizer(sc.get_fit_data(), train_size=1.0) 177 opt.train() 178 179 # get force constants 180 fcp = ForceConstantPotential(cs, opt.parameters) 181 fcs = fcp.get_force_constants(ase_supercell) 182 fc2 = fcs.get_fc_array(order=2) 183 184 return fc2 185 186 187def _decode_options(options): 188 """This is an example to parse options given in str. 189 190 When options = 'cutoff = 4.0', options is converted to {'cutoff': 4.0}. 191 192 In this implementation (can be modified), using phonopy command line 193 options, ``options`` is passed by --fc-calc-opt such as:: 194 195 phonopy --hiphiveph --fc-calc-opt "cutoff = 4" ... 196 197 """ 198 199 option_dict = {} 200 for pair in options.split(','): 201 key, value = [x.strip() for x in pair.split('=')] 202 if key == 'cutoff': 203 option_dict[key] = float(value) 204 return option_dict 205