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