%matplotlib inline
import numpy as np
import pylab as plt
from pyiron_atomistics import Project
pr=Project('FreeEnergy')
pot = '1995--Angelo-J-E--Ni-Al-H--LAMMPS--ipr1'
# pr.remove_jobs_silently(recursive=True)
alat = 3.5
basis = pr.create_ase_bulk('Ni')
ham = pr.create_job(pr.job_type.Lammps, 'ref')
ham.structure = basis
ham.potential = pot
murn = pr.create_job(job_type=pr.job_type.Murnaghan, job_name='murnaghan')
murn.ref_job = ham
murn.run()
The job murnaghan was saved and received the ID: 207 The job strain_0_9 was saved and received the ID: 208 The job strain_0_92 was saved and received the ID: 209 The job strain_0_94 was saved and received the ID: 210 The job strain_0_96 was saved and received the ID: 211 The job strain_0_98 was saved and received the ID: 212 The job strain_1_0 was saved and received the ID: 213 The job strain_1_02 was saved and received the ID: 214 The job strain_1_04 was saved and received the ID: 215 The job strain_1_06 was saved and received the ID: 216 The job strain_1_08 was saved and received the ID: 217 The job strain_1_1 was saved and received the ID: 218 job_id: 208 finished job_id: 209 finished job_id: 210 finished job_id: 211 finished job_id: 212 finished job_id: 213 finished job_id: 214 finished job_id: 215 finished job_id: 216 finished job_id: 217 finished job_id: 218 finished
print (murn['output/equilibrium_b_prime'])
print (murn['output/equilibrium_volume']**(1/3))
4.204514938124235 2.2174331116927366
murn.plot()
T_debye_low, T_debye_high = murn.fit.debye_temperature
print (T_debye_low)
[445.76166752 442.54156191 439.35906614 436.21357484 433.10449505 430.03124589 426.99325828 423.98997461 421.02084851 418.08534454 415.18293795 412.31311439 409.47536969 406.6692096 403.89414956 401.14971446 398.43543845 395.75086465 393.09554504 390.46904014 387.87091893 385.30075853 382.75814413 380.24266873 377.75393296 375.29154497 372.85512018 370.44428118 368.05865752 365.69788558 363.36160842 361.04947562 358.76114315 356.49627319 354.25453404 352.03559999 349.83915112 347.66487325 345.51245779 343.38160161 341.27200692 339.18338118 337.11543696 335.06789185 333.04046836 331.03289378 329.04490013 327.07622401 325.12660655 323.1957933 ]
plt.plot(murn.fit.volume, T_debye_high, label='high')
plt.plot(murn.fit.volume, T_debye_low, label='low')
plt.title('$T_{Debye}$ vs volume for Moruzzi low and high parameter')
plt.xlabel('Volume [$\AA^3$]')
plt.ylabel('$T_{Debye}$ [K]')
plt.legend();
pes = pr.create_object(object_type=pr.object_type.ThermoBulk)
pes.temperatures = np.linspace(1, 1500, 200)
pes.volumes = murn.fit.volume
pes.energies = murn.fit.polynomial()
pes.energies += murn.fit.energy_vib(T=pes.temperatures, low_T_limit=True).T
plt.plot(pes.volumes, pes.energies[0])
plt.plot(pes.volumes, pes.energies[30])
[<matplotlib.lines.Line2D at 0x7fc03a382460>]
pes.plot_contourf(show_min_erg_path=True);
pes_fine = pes.copy()
pes_fine.set_volumes(volume_min=murn.equilibrium_volume, volume_max=murn.equilibrium_volume*1.1, volume_steps=100)
murn.fit.volume = pes_fine.volumes
murn_fine = pes_fine.copy()
murn_fine.energies = murn.fit.polynomial()
pes_fine.energies = murn_fine.energies + murn.fit.energy_vib(T=pes_fine.temperatures, low_T_limit=True).T
print (pes_fine.temperatures)
[1.00000000e+00 8.53266332e+00 1.60653266e+01 2.35979899e+01 3.11306533e+01 3.86633166e+01 4.61959799e+01 5.37286432e+01 6.12613065e+01 6.87939698e+01 7.63266332e+01 8.38592965e+01 9.13919598e+01 9.89246231e+01 1.06457286e+02 1.13989950e+02 1.21522613e+02 1.29055276e+02 1.36587940e+02 1.44120603e+02 1.51653266e+02 1.59185930e+02 1.66718593e+02 1.74251256e+02 1.81783920e+02 1.89316583e+02 1.96849246e+02 2.04381910e+02 2.11914573e+02 2.19447236e+02 2.26979899e+02 2.34512563e+02 2.42045226e+02 2.49577889e+02 2.57110553e+02 2.64643216e+02 2.72175879e+02 2.79708543e+02 2.87241206e+02 2.94773869e+02 3.02306533e+02 3.09839196e+02 3.17371859e+02 3.24904523e+02 3.32437186e+02 3.39969849e+02 3.47502513e+02 3.55035176e+02 3.62567839e+02 3.70100503e+02 3.77633166e+02 3.85165829e+02 3.92698492e+02 4.00231156e+02 4.07763819e+02 4.15296482e+02 4.22829146e+02 4.30361809e+02 4.37894472e+02 4.45427136e+02 4.52959799e+02 4.60492462e+02 4.68025126e+02 4.75557789e+02 4.83090452e+02 4.90623116e+02 4.98155779e+02 5.05688442e+02 5.13221106e+02 5.20753769e+02 5.28286432e+02 5.35819095e+02 5.43351759e+02 5.50884422e+02 5.58417085e+02 5.65949749e+02 5.73482412e+02 5.81015075e+02 5.88547739e+02 5.96080402e+02 6.03613065e+02 6.11145729e+02 6.18678392e+02 6.26211055e+02 6.33743719e+02 6.41276382e+02 6.48809045e+02 6.56341709e+02 6.63874372e+02 6.71407035e+02 6.78939698e+02 6.86472362e+02 6.94005025e+02 7.01537688e+02 7.09070352e+02 7.16603015e+02 7.24135678e+02 7.31668342e+02 7.39201005e+02 7.46733668e+02 7.54266332e+02 7.61798995e+02 7.69331658e+02 7.76864322e+02 7.84396985e+02 7.91929648e+02 7.99462312e+02 8.06994975e+02 8.14527638e+02 8.22060302e+02 8.29592965e+02 8.37125628e+02 8.44658291e+02 8.52190955e+02 8.59723618e+02 8.67256281e+02 8.74788945e+02 8.82321608e+02 8.89854271e+02 8.97386935e+02 9.04919598e+02 9.12452261e+02 9.19984925e+02 9.27517588e+02 9.35050251e+02 9.42582915e+02 9.50115578e+02 9.57648241e+02 9.65180905e+02 9.72713568e+02 9.80246231e+02 9.87778894e+02 9.95311558e+02 1.00284422e+03 1.01037688e+03 1.01790955e+03 1.02544221e+03 1.03297487e+03 1.04050754e+03 1.04804020e+03 1.05557286e+03 1.06310553e+03 1.07063819e+03 1.07817085e+03 1.08570352e+03 1.09323618e+03 1.10076884e+03 1.10830151e+03 1.11583417e+03 1.12336683e+03 1.13089950e+03 1.13843216e+03 1.14596482e+03 1.15349749e+03 1.16103015e+03 1.16856281e+03 1.17609548e+03 1.18362814e+03 1.19116080e+03 1.19869347e+03 1.20622613e+03 1.21375879e+03 1.22129146e+03 1.22882412e+03 1.23635678e+03 1.24388945e+03 1.25142211e+03 1.25895477e+03 1.26648744e+03 1.27402010e+03 1.28155276e+03 1.28908543e+03 1.29661809e+03 1.30415075e+03 1.31168342e+03 1.31921608e+03 1.32674874e+03 1.33428141e+03 1.34181407e+03 1.34934673e+03 1.35687940e+03 1.36441206e+03 1.37194472e+03 1.37947739e+03 1.38701005e+03 1.39454271e+03 1.40207538e+03 1.40960804e+03 1.41714070e+03 1.42467337e+03 1.43220603e+03 1.43973869e+03 1.44727136e+03 1.45480402e+03 1.46233668e+03 1.46986935e+03 1.47740201e+03 1.48493467e+03 1.49246734e+03 1.50000000e+03]
pes_fine_high = pes_fine.copy()
pes_fine_high.energies = murn_fine.energies + murn.fit.energy_vib(T=pes_fine_high.temperatures,
low_T_limit=False).T
# print (pes_fine_high.temperatures)
ax = pes_fine_high.plot_contourf()
ax = pes_fine_high.plot_min_energy_path('g--', ax=ax, label='high T limit')
ax = pes_fine.plot_min_energy_path('k--', ax=ax, label='low T limit')
ax.legend(loc='lower right');
pes_fine.contour_pressure()
pes_fine.plot_free_energy()
pes_fine.plot_entropy()
pes_fine.num_atoms = 3
pes_fine.plot_heat_capacity(to_kB=True)
pes_QHA = pes_fine.copy()
pes_QHA.set_volumes(volume_min=np.min(pes_QHA.volumes), volume_max=np.max(pes_QHA.volumes), volume_steps=7)
alat_lst = pes_QHA.volumes**(1/3)
print (alat_lst)
[2.21743311 2.22968437 2.24180246 2.25379094 2.26565321 2.27739256 2.28901211]