In [1]:
%matplotlib inline
import numpy as np
import pylab as plt
from pyiron import Project
In [2]:
pr=Project('FreeEnergy')
pot = '1995--Angelo-J-E--Ni-Al-H--LAMMPS--ipr1'
# pr.remove_jobs_silently(recursive=True)
In [3]:
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
In [4]:
print (murn['output/equilibrium_b_prime'])
print (murn['output/equilibrium_volume']**(1/3))
4.204514938124235
2.2174331116927366
In [5]:
murn.plot()
In [6]:
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 ]
In [7]:
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();
In [8]:
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
In [9]:
plt.plot(pes.volumes, pes.energies[0])
plt.plot(pes.volumes, pes.energies[30])
Out[9]:
[<matplotlib.lines.Line2D at 0x7fc03a382460>]
In [10]:
pes.plot_contourf(show_min_erg_path=True);
In [11]:
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
In [12]:
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]
In [13]:
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
In [14]:
# print (pes_fine_high.temperatures)
In [15]:
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');
In [16]:
pes_fine.contour_pressure()
In [17]:
pes_fine.plot_free_energy()
In [18]:
pes_fine.plot_entropy()
In [19]:
pes_fine.num_atoms = 3
pes_fine.plot_heat_capacity(to_kB=True)
In [20]:
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]
In [ ]: