|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +# Copyright (C) 2011 Atsushi Togo |
| 4 | +# All rights reserved. |
| 5 | +# |
| 6 | +# This file is part of phonopy. |
| 7 | +# |
| 8 | +# Redistribution and use in source and binary forms, with or without |
| 9 | +# modification, are permitted provided that the following conditions |
| 10 | +# are met: |
| 11 | +# |
| 12 | +# * Redistributions of source code must retain the above copyright |
| 13 | +# notice, this list of conditions and the following disclaimer. |
| 14 | +# |
| 15 | +# * Redistributions in binary form must reproduce the above copyright |
| 16 | +# notice, this list of conditions and the following disclaimer in |
| 17 | +# the documentation and/or other materials provided with the |
| 18 | +# distribution. |
| 19 | +# |
| 20 | +# * Neither the name of the phonopy project nor the names of its |
| 21 | +# contributors may be used to endorse or promote products derived |
| 22 | +# from this software without specific prior written permission. |
| 23 | +# |
| 24 | +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
| 25 | +# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
| 26 | +# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS |
| 27 | +# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE |
| 28 | +# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, |
| 29 | +# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, |
| 30 | +# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; |
| 31 | +# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER |
| 32 | +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
| 33 | +# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN |
| 34 | +# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| 35 | +# POSSIBILITY OF SUCH DAMAGE. |
| 36 | + |
| 37 | +import numpy as np |
| 38 | +import sys |
| 39 | +from phonopy.units import EVAngstromToGPa |
| 40 | +from phonopy import PhonopyQHA |
| 41 | +from phonopy.file_IO import read_thermal_properties_yaml, read_v_e, read_efe |
| 42 | + |
| 43 | + |
| 44 | +# simplified file reader assuming plain text format for all data |
| 45 | +#def read_v_e(f): |
| 46 | +# return *(np.loadtxt(f).T) |
| 47 | + |
| 48 | +#def read_v_efe(f): |
| 49 | +# dat=np.loadtxt(f) |
| 50 | +# return dat[:, 0], dat[:, 1:] |
| 51 | + |
| 52 | +def read_thermal_properties_txt(fs): |
| 53 | + """ |
| 54 | + sample file |
| 55 | +# QHA Per primitive cell: T (K); E (eV); A=E-TS (eV); S (kB); Cv (kB) |
| 56 | +1.000000000000000000e+01 3.414113867041713990e-02 3.414009916110528897e-02 1.206300633020025714e-03 4.016837188045686326e-03 |
| 57 | +2.000000000000000000e+01 3.415606529685803749e-02 3.413557425809424922e-02 1.188943319227290553e-02 3.942826885954845567e-02 |
| 58 | + """ |
| 59 | + R= 8.31446261815324 |
| 60 | + eV2kJ = 96.4853329 |
| 61 | + dat = np.transpose(np.array([np.loadtxt(f) for f in fs]), (1,0,2)) |
| 62 | + return dat[:,0,0], dat[:,:,4]*R, dat[:,:,3]*R, dat[:,:,2]*eV2kJ, [], [] |
| 63 | + |
| 64 | +def get_options(): |
| 65 | + import argparse |
| 66 | + parser = argparse.ArgumentParser( |
| 67 | + description="Phonopy-QHA command-line-tool") |
| 68 | + parser.set_defaults(pressure=0.0, |
| 69 | + is_graph_plot=False, |
| 70 | + is_graph_save=False, |
| 71 | + is_bulk_modulus_only=False, |
| 72 | + efe_file=None, |
| 73 | + eos="vinet", |
| 74 | + thin_number=10, |
| 75 | + tmax=1000.0) |
| 76 | + parser.add_argument( |
| 77 | + "-b", dest="is_bulk_modulus_only", action="store_true", |
| 78 | + help="Just show Bulk modulus from v-e data") |
| 79 | + parser.add_argument( |
| 80 | + "--eos", dest="eos", |
| 81 | + help="Choise of EOS among vinet, birch_murnaghan, and murnaghan") |
| 82 | + parser.add_argument( |
| 83 | + "--exclude_imaginary", dest="exclude_imaginary", action="store_true", |
| 84 | + help="Exclude volumes that show imaginary modes") |
| 85 | + parser.add_argument( |
| 86 | + "-p", "--plot", dest="is_graph_plot", action="store_true", |
| 87 | + help="Plot data") |
| 88 | + parser.add_argument( |
| 89 | + "--pressure", dest="pressure", type=float, |
| 90 | + help="Pressure in GPa") |
| 91 | + parser.add_argument( |
| 92 | + "--efe", "--electronic-free-energy", |
| 93 | + dest="efe_file", nargs=1, |
| 94 | + help="Read electronic free energies at temperatures and volumes") |
| 95 | + parser.add_argument( |
| 96 | + "-s", "--save", dest="is_graph_save", action="store_true", |
| 97 | + help="Save plot data in pdf") |
| 98 | + parser.add_argument( |
| 99 | + "--sparse", dest="thin_number", type=int, |
| 100 | + help=("Thin out the F-V plots of temperature. The value is " |
| 101 | + "used as deviser of number of temperature points.")) |
| 102 | + parser.add_argument( |
| 103 | + "--tmax", dest="tmax", type=float, |
| 104 | + help="Maximum calculated temperature") |
| 105 | + parser.add_argument( |
| 106 | + "filenames", nargs='*', |
| 107 | + help="Filenames of e-v.dat and thermal_properties.yaml's") |
| 108 | + args = parser.parse_args() |
| 109 | + return args |
| 110 | + |
| 111 | + |
| 112 | +def main(args): |
| 113 | + if args.is_graph_save: |
| 114 | + import matplotlib |
| 115 | + matplotlib.use('Agg') |
| 116 | + matplotlib.rc('pdf', fonttype=42) |
| 117 | + |
| 118 | + # Choose EOS |
| 119 | + if args.eos == "birch_murnaghan": |
| 120 | + print("# Third-order Birch-Murnaghan EOS") |
| 121 | + elif args.eos == "murnaghan": |
| 122 | + print("# Murnaghan EOS") |
| 123 | + else: |
| 124 | + print("# Vinet EOS") |
| 125 | + |
| 126 | + # Show bulk modulus of v-e data |
| 127 | + if args.is_bulk_modulus_only: |
| 128 | + if args.efe_file: |
| 129 | + print("--efe optin can't be used with -b option") |
| 130 | + sys.exit(1) |
| 131 | + |
| 132 | + volumes, electronic_energies = read_v_e(args.filenames[0]) |
| 133 | + electronic_energies += volumes * args.pressure / EVAngstromToGPa |
| 134 | + bulk_modulus = PhonopyQHA(volumes, |
| 135 | + electronic_energies=electronic_energies, |
| 136 | + eos=args.eos) |
| 137 | + parameters = bulk_modulus.get_bulk_modulus_parameters() |
| 138 | + print("Volume: %f" % parameters[3]) |
| 139 | + print("Energy: %f" % parameters[0]) |
| 140 | + print("Bulk modulus: %f" % (parameters[1] * EVAngstromToGPa)) |
| 141 | + print("Parameters: %f %f %f %f" % tuple(parameters)) |
| 142 | + if args.is_graph_plot: |
| 143 | + bulk_modulus.plot_bulk_modulus().show() |
| 144 | + sys.exit(0) |
| 145 | + |
| 146 | + ######################## |
| 147 | + # Read data from files # |
| 148 | + ######################## |
| 149 | + volumes, electronic_energies = read_v_e(args.filenames[0]) |
| 150 | + electronic_energies += volumes * args.pressure / EVAngstromToGPa |
| 151 | + # Check number of files in e-v.dat case |
| 152 | + if len(volumes) != len(args.filenames[1:]): |
| 153 | + print("The number of thermal_properites.yaml files (%d) " |
| 154 | + "is inconsisten with" % len(args.filenames[1:])) |
| 155 | + print("the number of e-v data (%d)." % len(volumes)) |
| 156 | + sys.exit(1) |
| 157 | + |
| 158 | + if args.efe_file: |
| 159 | + _temperatures, electronic_energies = read_efe(args.efe_file[0]) |
| 160 | + if len(volumes) != electronic_energies.shape[1]: |
| 161 | + print("%s and %s are inconsistent for the volume points." |
| 162 | + % (args.filenames[0], args.efe_file[0])) |
| 163 | + sys.exit(1) |
| 164 | + |
| 165 | + if args.filenames[1][-4].lower() == 'yaml': |
| 166 | + dat_reader= read_thermal_properties_yaml |
| 167 | + else: |
| 168 | + dat_reader= read_thermal_properties_txt |
| 169 | + (temperatures, |
| 170 | + cv, |
| 171 | + entropy, |
| 172 | + fe_phonon, |
| 173 | + num_modes, |
| 174 | + num_integrated_modes) = dat_reader(args.filenames[1:]) |
| 175 | + |
| 176 | + if args.efe_file: |
| 177 | + if ((len(temperatures) >= len(_temperatures) and |
| 178 | + (np.abs(temperatures[:len(_temperatures)] - _temperatures) |
| 179 | + > 1e-5).any()) or |
| 180 | + (len(temperatures) < len(_temperatures) and |
| 181 | + (np.abs(temperatures - _temperatures[:len(temperatures)]) |
| 182 | + > 1e-5).any())): |
| 183 | + print("Inconsistency is found in temperatures in %s and " |
| 184 | + "thermal_properties.yaml files." % args.filenames[0]) |
| 185 | + sys.exit(1) |
| 186 | + |
| 187 | + ######################################################## |
| 188 | + # Treatment of thermal properties with imaginary modes # |
| 189 | + ######################################################## |
| 190 | + if args.exclude_imaginary and num_modes: |
| 191 | + indices = [] |
| 192 | + num_imag_modes = np.array(num_modes) - np.array(num_integrated_modes) |
| 193 | + for i, nim in enumerate(num_imag_modes): |
| 194 | + if nim < 4: |
| 195 | + indices.append(i) |
| 196 | + else: |
| 197 | + indices = range(len(volumes)) |
| 198 | + |
| 199 | + if args.efe_file: |
| 200 | + electronic_energies = electronic_energies[:, indices] |
| 201 | + else: |
| 202 | + electronic_energies = electronic_energies[indices] |
| 203 | + |
| 204 | + ########################## |
| 205 | + # Analyzing and plotting # |
| 206 | + ########################## |
| 207 | + phonopy_qha = PhonopyQHA(volumes=volumes[indices], |
| 208 | + electronic_energies=electronic_energies, |
| 209 | + eos=args.eos, |
| 210 | + temperatures=temperatures, |
| 211 | + free_energy=fe_phonon[:, indices], |
| 212 | + cv=cv[:, indices], |
| 213 | + entropy=entropy[:, indices], |
| 214 | + t_max=args.tmax, |
| 215 | + verbose=True) |
| 216 | + |
| 217 | + if num_modes: |
| 218 | + num_imag_modes = np.array(num_modes) - np.array(num_integrated_modes) |
| 219 | + for filename, nim in zip(args.filenames[1:(len(volumes)+1)], |
| 220 | + num_imag_modes): |
| 221 | + if nim > 3: |
| 222 | + if args.exclude_imaginary: |
| 223 | + print("# %s has been excluded." % filename) |
| 224 | + else: |
| 225 | + print("# Warning: %s has imaginary modes." % filename) |
| 226 | + |
| 227 | + if args.is_graph_plot and not args.is_graph_save: |
| 228 | + # Plot on display |
| 229 | + # - Volume vs Helmholtz free energy |
| 230 | + # - Volume vs Temperature |
| 231 | + # - Thermal expansion coefficient |
| 232 | + phonopy_qha.plot_qha(thin_number=args.thin_number).show() |
| 233 | + |
| 234 | + if args.is_graph_save: |
| 235 | + # Volume vs Helmholts free energy |
| 236 | + phonopy_qha.plot_pdf_helmholtz_volume(thin_number=args.thin_number) |
| 237 | + |
| 238 | + # Volume vs Temperature |
| 239 | + phonopy_qha.plot_pdf_volume_temperature() |
| 240 | + |
| 241 | + # Thermal expansion coefficient |
| 242 | + phonopy_qha.plot_pdf_thermal_expansion() |
| 243 | + |
| 244 | + # G vs Temperature |
| 245 | + phonopy_qha.plot_pdf_gibbs_temperature() |
| 246 | + |
| 247 | + # Bulk modulus vs Temperature |
| 248 | + phonopy_qha.plot_pdf_bulk_modulus_temperature() |
| 249 | + |
| 250 | + # C_P vs Temperature |
| 251 | + phonopy_qha.plot_pdf_heat_capacity_P_numerical() |
| 252 | + |
| 253 | + # C_P vs Temperature (poly fit) |
| 254 | + phonopy_qha.plot_pdf_heat_capacity_P_polyfit() |
| 255 | + |
| 256 | + # Gruneisen parameter vs Temperature |
| 257 | + phonopy_qha.plot_pdf_gruneisen_temperature() |
| 258 | + |
| 259 | + phonopy_qha.write_helmholtz_volume() |
| 260 | + phonopy_qha.write_volume_temperature() |
| 261 | + phonopy_qha.write_thermal_expansion() |
| 262 | + phonopy_qha.write_gibbs_temperature() |
| 263 | + phonopy_qha.write_bulk_modulus_temperature() |
| 264 | + phonopy_qha.write_heat_capacity_P_numerical() |
| 265 | + phonopy_qha.write_heat_capacity_P_polyfit() |
| 266 | + phonopy_qha.write_gruneisen_temperature() |
| 267 | + |
| 268 | + |
| 269 | +if __name__ == "__main__": |
| 270 | + main(get_options()) |
0 commit comments