Skip to content

Commit

Permalink
add compatibility with simpson scipy function
Browse files Browse the repository at this point in the history
  • Loading branch information
abelcarreras committed Nov 7, 2024
1 parent 43ec150 commit a4a9d3c
Showing 1 changed file with 7 additions and 7 deletions.
14 changes: 7 additions & 7 deletions dynaphopy/analysis/thermal_properties.py
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ def n(temp, freq):
total_energy = np.nan_to_num([dos[i] * h_bar * freq * (0.5 + n(temperature, freq))
for i, freq in enumerate(frequency)])

total_energy = simps(total_energy, frequency) * N_a / 1000 # KJ/K/mol
total_energy = simps(total_energy, x=frequency) * N_a / 1000 # KJ/K/mol
return total_energy


Expand All @@ -52,7 +52,7 @@ def get_free_energy(temperature, frequency, dos):
for i, freq in enumerate(frequency)])

free_energy[0] = 0
free_energy = simps(free_energy, frequency) * N_a / 1000 # KJ/K/mol
free_energy = simps(free_energy, x=frequency) * N_a / 1000 # KJ/K/mol
return free_energy


Expand All @@ -64,7 +64,7 @@ def n(temp, freq):
free_energy_c = np.nan_to_num([dos[i] * -h_bar/2 *shift*(n(temperature, freq) + 1 / 2.)
for i, freq in enumerate(frequency)])

free_energy_c = simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
free_energy_c = simps(free_energy_c, x=frequency) * N_a / 1000 # KJ/K/mol
return free_energy_c


Expand All @@ -81,7 +81,7 @@ def n(temp, freq):

free_energy_c = free_energy_1 - free_energy_2

free_energy_c = simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
free_energy_c = simps(free_energy_c, x=frequency) * N_a / 1000 # KJ/K/mol
return free_energy_c


Expand All @@ -93,7 +93,7 @@ def coth(x):
entropy = np.nan_to_num([dos[i]*(1.0 / (2. * temperature) * h_bar * freq * coth(h_bar * freq / (2 * k_b * temperature))
- k_b * np.log(2 * np.sinh(h_bar * freq / (2 * k_b * temperature))))
for i, freq in enumerate(frequency)])
entropy = simps(entropy, frequency) * N_a # J/K/mol
entropy = simps(entropy, x=frequency) * N_a # J/K/mol
return entropy

# Alternative way to calculate entropy (not used)
Expand All @@ -105,7 +105,7 @@ def n(temp, freq):
entropy = np.nan_to_num([dos[i] * k_b * ((n(temperature, freq) + 1) * np.log(n(temperature, freq) + 1)
- n(temperature, freq) * np.log(n(temperature, freq)))
for i, freq in enumerate(frequency)])
entropy = simps(entropy, frequency) * N_a # J/K/mol
entropy = simps(entropy, x=frequency) * N_a # J/K/mol
return entropy


Expand All @@ -116,7 +116,7 @@ def z(temp, freq):

c_v = np.nan_to_num([dos[i] * k_b * pow(z(temperature, freq), 2) * np.exp(z(temperature, freq)) / pow(np.exp(z(temperature, freq)) - 1, 2)
for i, freq in enumerate(frequency)])
c_v = simps(c_v, frequency) * N_a # J/K/mol
c_v = simps(c_v, x=frequency) * N_a # J/K/mol

return c_v

Expand Down

0 comments on commit a4a9d3c

Please sign in to comment.