Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Restyle feat: Add entropy output #65

Merged
merged 2 commits into from
Mar 12, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 43 additions & 21 deletions qha/calculator.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
from qha.single_configuration import free_energy
from qha.thermodynamics import *
from qha.type_aliases import Vector
from qha.unit_conversion import gpa_to_ry_b3, ry_b3_to_gpa, b3_to_a3, ry_to_j_mol, ry_to_ev
from qha.unit_conversion import gpa_to_ry_b3, ry_b3_to_gpa, b3_to_a3, ry_to_j_mol, ry_to_ev, ry_to_j
from qha.v2p import v2p

# ===================== What can be exported? =====================
Expand All @@ -48,7 +48,8 @@ def __init__(self, user_settings: Dict[str, Any]):
try:
runtime_settings.update({key: user_settings[key]})
except KeyError:
continue # If a key is not set in user settings, use the default.
# If a key is not set in user settings, use the default.
continue

self._settings = runtime_settings

Expand Down Expand Up @@ -106,12 +107,15 @@ def v_ratio(self) -> Optional[float]:

def read_input(self):
try:
formula_unit_number, volumes, static_energies, frequencies, q_weights = read_input(self.settings['input'])
formula_unit_number, volumes, static_energies, frequencies, q_weights = read_input(
self.settings['input'])
except KeyError:
raise KeyError("The 'input' option must be given in your settings!")
raise KeyError(
"The 'input' option must be given in your settings!")

if not qha.tools.is_monotonic_decreasing(volumes):
raise RuntimeError("Check the input file to make sure the volume decreases!")
raise RuntimeError(
"Check the input file to make sure the volume decreases!")

self._formula_unit_number: int = formula_unit_number
self._volumes = volumes
Expand All @@ -127,7 +131,8 @@ def where_negative_frequencies(self) -> Optional[Vector]:
:return:
"""
if self._frequencies is None:
print("Please invoke ``read_input`` method first!") # ``None`` is returned
# ``None`` is returned
print("Please invoke ``read_input`` method first!")
else:
_ = np.transpose(np.where(self._frequencies < 0))
if _.size == 0:
Expand Down Expand Up @@ -156,20 +161,23 @@ def refine_grid(self):
try:
p_min, p_min_modifier, ntv, order = d['P_MIN'], d['p_min_modifier'], d['NTV'], d['order']
except KeyError:
raise KeyError("All the 'P_MIN', 'p_min_modifier', 'NTV', 'order' options must be given in your settings!")
raise KeyError(
"All the 'P_MIN', 'p_min_modifier', 'NTV', 'order' options must be given in your settings!")

r = FinerGrid(p_min - p_min_modifier, ntv, order=order)

if 'volume_ratio' in d:
self._finer_volumes_bohr3, self._f_tv_ry, self._v_ratio = r.refine_grid(self.volumes, self.vib_ry,
ratio=d['volume_ratio'])
else:
self._finer_volumes_bohr3, self._f_tv_ry, self._v_ratio = r.refine_grid(self.volumes, self.vib_ry)
self._finer_volumes_bohr3, self._f_tv_ry, self._v_ratio = r.refine_grid(
self.volumes, self.vib_ry)

@LazyProperty
def vib_ry(self):
# We grep all the arguments once since they are being invoked for thousands of times, and will be an overhead.
args = self.q_weights, self.static_energies, self.frequencies, self.settings['static_only']
args = self.q_weights, self.static_energies, self.frequencies, self.settings[
'static_only']

mat = np.empty((self.temperature_array.size, self.volumes.size))
for i, t in enumerate(self.temperature_array):
Expand Down Expand Up @@ -218,7 +226,8 @@ def desired_pressure_status(self) -> None:
self.p_tv_gpa[:, 0].max(), self.p_tv_gpa[:, -1].min()))

if self.p_tv_gpa[:, -1].min() < self.desired_pressures_gpa.max():
ntv_max = int((self.p_tv_gpa[:, -1].min() - self.desired_pressures_gpa.min()) / d['DELTA_P'])
ntv_max = int(
(self.p_tv_gpa[:, -1].min() - self.desired_pressures_gpa.min()) / d['DELTA_P'])

save_to_output(d['qha_output'], textwrap.dedent("""\
!!!ATTENTION!!!
Expand All @@ -228,12 +237,17 @@ def desired_pressure_status(self) -> None:
Please reduce the NTV accordingly, for example, try to set NTV < {:4d}.
""".format(ntv_max)))

raise ValueError("DESIRED PRESSURE is too high (NTV is too large), qha results might not be right!")
raise ValueError(
"DESIRED PRESSURE is too high (NTV is too large), qha results might not be right!")

@LazyProperty
def p_tv_au(self):
return pressure(self.finer_volumes_bohr3, self.f_tv_ry)

@LazyProperty
def s_tv_j(self):
return ry_to_j(entropy(self.temperature_array, self.f_tv_ry))

@LazyProperty
def f_tv_ev(self):
return ry_to_ev(self.f_tv_ry)
Expand Down Expand Up @@ -373,12 +387,15 @@ def read_input(self):
q_weights = []

for inp in input_data_files:
nm_tmp, volumes_tmp, static_energies_tmp, freq_tmp, weights_tmp = read_input(inp)
nm_tmp, volumes_tmp, static_energies_tmp, freq_tmp, weights_tmp = read_input(
inp)

if not qha.tools.is_monotonic_decreasing(volumes_tmp):
# TODO: Clean this sentence
save_to_output(self.settings['qha_output'], "Check the input file to make sure the volume decreases")
raise ValueError("Check the input file to make sure the volumes are monotonic decreasing!")
save_to_output(
self.settings['qha_output'], "Check the input file to make sure the volume decreases")
raise ValueError(
"Check the input file to make sure the volumes are monotonic decreasing!")

formula_unit_numbers.append(nm_tmp)
volumes.append(volumes_tmp)
Expand All @@ -393,12 +410,15 @@ def read_input(self):
q_weights = np.array(q_weights)

if not len(set(formula_unit_numbers)) == 1:
raise RuntimeError("All the formula unit number in all inputs should be the same!")
raise RuntimeError(
"All the formula unit number in all inputs should be the same!")

if len(volumes.shape) == 1:
raise RuntimeError("All configurations should have same number of volumes!")
raise RuntimeError(
"All configurations should have same number of volumes!")

self._formula_unit_number = formula_unit_numbers[0] # Choose any of them since they are all the same
# Choose any of them since they are all the same
self._formula_unit_number = formula_unit_numbers[0]
self._volumes = volumes
self._static_energies = static_energies
self._frequencies = frequencies
Expand All @@ -408,11 +428,12 @@ def read_input(self):
def vib_ry(self):
# We grep all the arguments once since they are being invoked for thousands of times, and will be an overhead.
args = self.degeneracies, self.q_weights, self.static_energies, self._volumes, self.frequencies, \
self.settings['static_only']
self.settings['static_only']

mat = np.empty((self.temperature_array.size, self._volumes.shape[1]))
for i, t in enumerate(self.temperature_array):
mat[i] = different_phonon_dos.PartitionFunction(t, *(arg for arg in args)).get_free_energies()
mat[i] = different_phonon_dos.PartitionFunction(
t, *(arg for arg in args)).get_free_energies()

return mat

Expand All @@ -424,9 +445,10 @@ def __init__(self, user_settings):
@LazyProperty
def vib_ry(self):
args = self.degeneracies, self.q_weights[0], self.static_energies, self._volumes, self.frequencies[0], \
self.settings['static_only'], self.settings['order']
self.settings['static_only'], self.settings['order']
mat = np.empty((self.temperature_array.size, self._volumes.shape[1]))

for i, t in enumerate(self.temperature_array):
mat[i] = same_phonon_dos.FreeEnergy(t, *(arg for arg in args)).get_free_energies()
mat[i] = same_phonon_dos.FreeEnergy(
t, *(arg for arg in args)).get_free_energies()
return mat
43 changes: 30 additions & 13 deletions qha/cli/runner.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,8 @@ def run(self, namespace):
if not os.path.exists(user_settings['output_directory']):
os.makedirs(user_settings['output_directory'])

user_settings.update({'qha_output': os.path.join(user_settings['output_directory'], 'output.txt')})
user_settings.update({'qha_output': os.path.join(
user_settings['output_directory'], 'output.txt')})

try:
os.remove(user_settings['qha_output'])
Expand All @@ -64,10 +65,12 @@ def run(self, namespace):
print("You have single-configuration calculation assumed.")
elif calculation_type == 'same phonon dos':
calc = SamePhDOSCalculator(user_settings)
print("You have multi-configuration calculation with the same phonon DOS assumed.")
print(
"You have multi-configuration calculation with the same phonon DOS assumed.")
elif calculation_type == 'different phonon dos':
calc = DifferentPhDOSCalculator(user_settings)
print("You have multi-configuration calculation with different phonon DOS assumed.")
print(
"You have multi-configuration calculation with different phonon DOS assumed.")
else:
raise ValueError("The 'calculation' in your settings in not recognized! It must be one of:"
"'single', 'same phonon dos', 'different phonon dos'!")
Expand All @@ -81,7 +84,8 @@ def run(self, namespace):

print("Caution: If negative frequencies found, they are currently treated as 0!")
tmp = calc.where_negative_frequencies
if tmp is not None and not (tmp.T[-1].max() <= 2): # Don't delete this parenthesis!
# Don't delete this parenthesis!
if tmp is not None and not (tmp.T[-1].max() <= 2):
if calc.frequencies.ndim == 4: # Multiple configuration
for indices in tmp:
print(
Expand Down Expand Up @@ -122,20 +126,29 @@ def run(self, namespace):
}

file_ftv_fitted = results_folder / 'f_tv_fitted_ev_ang3.txt'
save_x_tv(calc.f_tv_ev, temperature_array, calc.finer_volumes_ang3, temperature_sample, file_ftv_fitted)
save_x_tv(calc.f_tv_ev, temperature_array,
calc.finer_volumes_ang3, temperature_sample, file_ftv_fitted)

file_ftv_non_fitted = results_folder / 'f_tv_nonfitted_ev_ang3.txt'
save_x_tv(calc.vib_ev, temperature_array, calc.volumes_ang3, temperature_sample, file_ftv_non_fitted)
save_x_tv(calc.vib_ev, temperature_array, calc.volumes_ang3,
temperature_sample, file_ftv_non_fitted)

file_ptv_gpa = results_folder / 'p_tv_gpa.txt'
save_x_tv(calc.p_tv_gpa, temperature_array, calc.finer_volumes_ang3, temperature_sample, file_ptv_gpa)
save_x_tv(calc.p_tv_gpa, temperature_array,
calc.finer_volumes_ang3, temperature_sample, file_ptv_gpa)

file_stv_j = results_folder / 's_tv_j.txt'
save_x_tv(calc.s_tv_j, temperature_array,
calc.finer_volumes_ang3, temperature_sample, file_stv_j)

for idx in calc.settings['thermodynamic_properties']:
if idx in ['F', 'G', 'H', 'U']:
attr_name = calculation_option[idx] + '_' + calc.settings['energy_unit']
attr_name = calculation_option[idx] + \
'_' + calc.settings['energy_unit']
file_name = attr_name + '.txt'
file_dir = results_folder / file_name
save_x_tp(getattr(calc, attr_name), temperature_array, desired_pressures_gpa, p_sample_gpa, file_dir)
save_x_tp(getattr(calc, attr_name), temperature_array,
desired_pressures_gpa, p_sample_gpa, file_dir)

if idx == 'V':
v_bohr3 = calculation_option[idx] + '_' + 'bohr3'
Expand All @@ -145,15 +158,19 @@ def run(self, namespace):
file_name_ang3 = v_ang3 + '.txt'
file_dir_ang3 = results_folder / file_name_ang3

save_x_tp(getattr(calc, v_bohr3), temperature_array, desired_pressures_gpa, p_sample_gpa, file_dir_au)
save_x_tp(getattr(calc, v_ang3), temperature_array, desired_pressures_gpa, p_sample_gpa, file_dir_ang3)
save_x_tp(getattr(calc, v_bohr3), temperature_array,
desired_pressures_gpa, p_sample_gpa, file_dir_au)
save_x_tp(getattr(calc, v_ang3), temperature_array,
desired_pressures_gpa, p_sample_gpa, file_dir_ang3)

if idx in ['Cv', 'Cp', 'Bt', 'Btp', 'Bs', 'alpha', 'gamma']:
attr_name = calculation_option[idx]
file_name = attr_name + '.txt'
file_dir = results_folder / file_name
save_x_tp(getattr(calc, attr_name), temperature_array, desired_pressures_gpa, p_sample_gpa, file_dir)
save_x_tp(getattr(calc, attr_name), temperature_array,
desired_pressures_gpa, p_sample_gpa, file_dir)

end_time_total = time.time()
time_elapsed = end_time_total - start_time_total
save_to_output(user_settings['qha_output'], make_ending_string(time_elapsed))
save_to_output(user_settings['qha_output'],
make_ending_string(time_elapsed))