#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('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()


# In[4]:


print (murn['output/equilibrium_b_prime'])
print (murn['output/equilibrium_volume']**(1/3))


# In[5]:


murn.plot()


# In[6]:


T_debye_low, T_debye_high = murn.fit.debye_temperature
print (T_debye_low)


# 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])


# 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)


# 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)


# In[ ]: