from aiida import load_profile, orm
load_profile('FermiEnergyPaper')
from aiida.orm import load_node, load_code, StructureData, load_group
from aiida.plugins import WorkflowFactory, DataFactory, CalculationFactory
from aiida.engine import submit
from aiida.tools import delete_nodes
from pymatgen.core.structure import Structure
import numpy as np
import matplotlib.pyplot as plt
group_label = 'force_vs_degauss_gauss_FD2.565'
code = load_code(label='QE6.8-pw@prnmarvelsrv4')
PwBaseWorkChain = WorkflowFactory('quantumespresso.pw.base')
#High resolution
params3 = {'font.size': 18,
'font.family': 'serif',
'legend.fontsize': 12,
'mathtext.fontset': 'dejavuserif',
'figure.figsize': (7, 4.5) ,
# 'figure.dpi' : 200,
'figure.dpi' : 100,
'legend.handlelength': 2}
marker_list = [ 'o', 's', '^', 'v', '*', 'x', 'P', 'D', 'h', 'p', '', '', '', '', '', '', '', '', '', '', '' ]
def query(group_label, smearing, degauss) :
from aiida import orm
PwCalculation = CalculationFactory('quantumespresso.pw')
qb = orm.QueryBuilder()
qb.append(orm.Group , tag='group', filters={'label': group_label})
qb.append(orm.WorkChainNode, tag='base_workchains' , with_group='group', filters={'attributes.process_label': 'PwBaseWorkChain', 'attributes.exit_status': 0}, project='*')
qb.append(orm.Dict, with_outgoing='base_workchains', edge_filters = {'label': 'pw__parameters'}, filters={ 'and' : [
{'attributes.SYSTEM.degauss' : {'==' : degauss } },
{'attributes.SYSTEM.smearing' : {'==' : smearing} } ] } )
print(' Total number of entries:', qb.count() )
return qb
corrected_energy = -2149.849420861203
def plot_kpoints_convergence(smearing, degauss_list, property='energy') :
plt.rcParams.update(params3)
marker_counter = 0
print('smearing=', smearing)
xmin= 3; xmax= 20
if property == 'energy' :
plt.plot( [xmin,xmax], [corrected_energy,corrected_energy] , linestyle=':', color='gray')
elif property in ['stress']:
plt.plot( [xmin,xmax], [0,0] , linestyle=':', color='gray')
for degauss in degauss_list :
print(' degauss=', degauss)
qb = query(group_label, smearing, degauss)
ys = []
xs = []
for item in qb.all() :
chain = item[0]
calculation_smearing = chain.inputs.pw__parameters['SYSTEM']['smearing']
calculation_degauss = chain.inputs.pw__parameters['SYSTEM']['degauss']
kpoint_mesh = chain.inputs.kpoints.get_kpoints_mesh()[0][0]
if kpoint_mesh >= xmin :
if property == 'energy' :
if smearing == 'gaussian' :
total_energy = chain.outputs.output_parameters.get_dict()['energy'] - 0.5*chain.outputs.output_parameters.get_dict()['energy_smearing']
else :
total_energy = chain.outputs.output_parameters.get_dict()['energy']
ys.append(total_energy)
elif property == 'stress':
stress_tensor = chain.outputs.output_trajectory.get_array('stress')[-1]
stress = np.trace( stress_tensor )
ys.append(stress)
elif property == 'force':
# total force: kind of the sum of the absolute values of all forces components of all atoms
force = chain.outputs.output_parameters.attributes['total_force']
ys.append(force)
elif property == 'force2':
# x component of force in atom 0
force = chain.outputs.output_trajectory.get_array('forces')[0][0][0]
ys.append(force)
else :
raise( 'Property not recognized. The available properties are: energy, stress, force.' )
xs.append(kpoint_mesh)
# print(calculation.id)
if len(xs) > 0:
xs, ys = zip(*sorted(zip(xs, ys)))
plt.plot( xs, ys, marker=marker_list[marker_counter], linewidth=2, label=smearing+' σ='+str(degauss) + ' Ry')
marker_counter += 1
# print(ys[-1])
plt.xticks( xs )
plt.xticks( [i for i in range(1,21,2)] )
deltay = 0.10
if property == 'energy' :
plt.ylim(corrected_energy-deltay, corrected_energy+deltay)
plt.ylabel("Total energy (eV)")
elif property == 'stress':
plt.ylabel("Tr[ stress tensor ] (GPascal)")
# plt.yscale('log')
# plt.ylim(-0.05, 0.5)
elif property == 'force':
plt.ylabel("Total force (eV/Å)")
elif property == 'force2':
plt.ylabel("F_x of atom 1 (eV/Å)")
plt.xlim(xmin, xmax)
plt.xlabel("kpoint mesh [#, #, #]")
plt.legend()
plt.tight_layout()
plt.savefig('force_kpoint.png')
plt.show()
def plot_kpoints_convergence_iteration(smearing, degauss_list) :
plt.rcParams.update(params3)
marker_counter = 0
print('smearing=', smearing)
xmin= 1; xmax= 20
# plt.plot( [xmin,xmax], [corrected_energy,corrected_energy] , linestyle=':', color='gray')
for degauss in degauss_list :
print(' degauss=', degauss)
qb = query(group_label, smearing, degauss)
ys = []
xs = []
for item in qb.all() :
chain = item[0]
calculation_smearing = chain.inputs.pw__parameters['SYSTEM']['smearing']
calculation_degauss = chain.inputs.pw__parameters['SYSTEM']['degauss']
kpoint_mesh = chain.inputs.kpoints.get_kpoints_mesh()[0][0]
if kpoint_mesh >= xmin :
try :
num_iteration = chain.outputs.output_trajectory.get_array('scf_iterations')
# print(num_iteration)
if len(chain.outputs.output_trajectory.get_array('scf_iterations')) > 1 :
print('array scf_iterations has more than one value:', chain.outputs.output_trajectory.get_array('scf_iterations') )
ys.append(num_iteration)
xs.append(kpoint_mesh)
# print(calculation.id)
except :
pass
if len(xs) > 0:
xs, ys = zip(*sorted(zip(xs, ys)))
plt.plot( xs, ys, marker=marker_list[marker_counter], linewidth=2, label=smearing+' σ='+str(degauss) + ' Ry')
marker_counter += 1
# print(ys[-1])
# plt.xticks( [ i for i in range( xs[0], xs[-1]+1 ) ] )
# plt.xticks( [ i for i in range( xs[0], 11 ) ] + [ i for i in range( 12, xs[-1]+1, 2 ) ] )
plt.xticks( [ i for i in range(1, 21, 2 ) ] )
deltay = 0.10
# plt.ylim(corrected_energy-deltay, corrected_energy+deltay)
plt.xlim(xmin, xmax)
plt.xlabel("kpoint mesh [#, #, #]")
plt.yscale('log')
plt.ylabel("Num. iterations")
plt.legend()
plt.tight_layout()
# plt.savefig('iteration_kpoint.png')
plt.show()
qb = query('kpoint_convergence', 'gaussian', 0.01)
Total number of entries: 20
group_label = 'kpoint_convergence'
plot_kpoints_convergence_iteration('cold', [0.000002, 0.01, 0.02])
smearing= cold degauss= 2e-06 Total number of entries: 18 degauss= 0.01 Total number of entries: 20 degauss= 0.02 Total number of entries: 20
def query_energy_smearing(group_label, smearing) :
from aiida import orm
PwCalculation = CalculationFactory('quantumespresso.pw')
qb = orm.QueryBuilder()
qb.append(orm.Group , tag='group', filters={'label': group_label})
qb.append(orm.WorkChainNode, tag='base_workchains' , with_group='group', filters={'attributes.process_label': 'PwBaseWorkChain', 'attributes.exit_status': 0 }, project='*')
qb.append(orm.Dict, with_outgoing='base_workchains', edge_filters = {'label': 'pw__parameters'}, filters={ 'attributes.SYSTEM.smearing' : smearing } )
qb.append(orm.KpointsData, with_outgoing='base_workchains', edge_filters = {'label': 'kpoints'}, project='*' )
print(' Total number of entries:', qb.count() )
return qb
corrected_energy = -2149.849420861203
#High resolution for the paper
params4 = {'font.size': 18,
'font.family': 'serif',
'legend.fontsize': 12,
'mathtext.fontset': 'dejavuserif',
'figure.figsize': (8, 4.5) ,
# 'figure.dpi' : 200,
'figure.dpi' : 100,
'legend.handlelength': 2}
def plot_energy_vs_smearing2(smearing_list, kpoint_mesh, property, scale=2.25) :
plt.rcParams.update(params4)
marker_counter = 0
xmin= 0; xmax= 0.1
if property == 'energy' :
# plt.plot( [xmin,xmax], [corrected_energy,corrected_energy] , linestyle=':', color='gray')
pass
elif property in ['stressSSSS']:
plt.plot( [xmin,xmax], [0,0] , linestyle=':', color='gray')
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
for smearing in smearing_list :
print(' smearing=', smearing)
qb = query_energy_smearing(group_label, smearing)
ys = []
ys2 = []
xs = []
ys_28 = []
xs_28 = []
ys_54 = []
num_bands_list = []
for item in qb.all() :
chain = item[0]
kpoint = item[1]
if kpoint.attributes['mesh'][0] == kpoint_mesh :
mesh = kpoint.attributes['mesh']
calculation_degauss = chain.inputs.pw__parameters['SYSTEM']['degauss']
number_of_bands = chain.outputs.output_parameters.get_dict()['number_of_bands']
num_bands_list.append(number_of_bands)
if property == 'energy' :
if smearing in ['gaussian', 'fermi-dirac' ] :
total_energy = chain.outputs.output_parameters.get_dict()['energy']
total_energy2 = chain.outputs.output_parameters.get_dict()['energy'] - 0.5*chain.outputs.output_parameters.get_dict()['energy_smearing']
ys2.append(total_energy2)
force = chain.outputs.output_parameters.attributes['total_force']
ys_54.append(force)
else :
total_energy = chain.outputs.output_parameters.get_dict()['energy']
ys.append(total_energy)
else :
raise( 'Property not recognized. The available properties are: energy, stress, force.' )
if(smearing == 'fermi-dirac') :
xs.append(calculation_degauss*scale)
# xs.append(calculation_degauss*2.25)
else :
xs.append(calculation_degauss)
if kpoint.attributes['mesh'][0] == 28 :
mesh = kpoint.attributes['mesh']
calculation_degauss = chain.inputs.pw__parameters['SYSTEM']['degauss']
if property == 'stress':
stress_tensor = chain.outputs.output_trajectory.get_array('stress')[-1]
# print(stress_tensor)
# print(chain.id)
stress = np.trace( stress_tensor ) / 3.0
# print('stress', stress)
ys_28.append(stress)
xs_28.append(calculation_degauss)
if len(xs) > 0:
if smearing in ['gaussian', 'fermi-dirac'] and property == 'energy' :
# xs, ys, ys2 = zip(*sorted(zip(xs, ys, ys2)))
xs, ys, ys2, ys_54 = zip(*sorted(zip(xs, ys, ys2, ys_54)))
else :
xs, ys = zip(*sorted(zip(xs, ys)))
if smearing in ['cold', 'methfessel-paxton'] :
xs_28, ys_28 = zip(*sorted(zip(xs_28, ys_28)))
if smearing == 'methfessel-paxton' :
label = 'Methfessel-Paxton'
elif smearing == 'fermi-dirac' :
label = 'Fermi-Dirac'
elif smearing == 'cold' :
label = 'Cold'
elif smearing == 'gaussian' :
label = 'Gaussian'
else :
label = smearing
style_list = ['-', '--']
style_list = ['', '']
color_list1 = ['tab:blue', 'tab:pink']
color_list2 = ['tab:red', 'tab:orange']
size_list = [8, 6]
marker_list = ['s','o']
ax1.plot( xs, np.array(ys) + 2.1498E3 + 0.026, color=color_list1[marker_counter], marker=marker_list[marker_counter], linestyle=style_list[marker_counter], label='Energy - '+label, linewidth=2*(marker_counter+1), markersize=size_list[marker_counter] )
if smearing in ['gaussian', 'fermi-dirac'] and property == 'energy':
ax2.plot( xs, ys_54, color=color_list2[marker_counter], marker=marker_list[marker_counter], linestyle=style_list[marker_counter], label= 'Force - '+label, linewidth=2*(marker_counter+1), markersize=size_list[marker_counter])
# # if smearing in ['cold', 'methfessel-paxton'] :
# # plt.plot( xs_28, ys_28, marker=marker_list[marker_counter], label=label+' 28x28x28', linewidth=2 )
marker_counter += 1
print( 'kpoint mesh=', mesh )
print( 'num_bands_list=',num_bands_list)
deltay = 0.80
if property == 'energy' :
# # plt.ylim(corrected_energy-deltay, corrected_energy+deltay*0.1)
# # plt.autoscale
# plt.xlim(xs[0], 0.02)
# # plt.ylim(-2149.8475, -2149.8490)
# # if not isinstance(kpoint_mesh, list) :
# # # plt.ylim(-2150.1, -2149.8)
# # plt.ylim(-2149.92, -2149.84) # new range
# plt.ylabel("Energy (eV)")
# ax1.xlim(xs[0], 0.02)
ax1.set_xlim(-0.001, 0.021)
color = 'tab:blue'
ax1.set_ylabel("Energy (eV)", color=color)
ax1.tick_params(axis='y', labelcolor=color)
color = 'tab:red'
ax2.set_ylabel("Total force (eV/Å)", color=color)
ax2.tick_params(axis='y', labelcolor=color)
elif property == 'stress':
# plt.ylabel("Tr[ stress tensor ] (GPascal)")
plt.ylabel("Pressure (GPascal)")
plt.xlim(xs[0], 0.02)
# plt.ylim(-0.15, 0.2)
# plt.yscale('log')
# plt.ylim(-0.05, 0.5)
elif property == 'force':
plt.ylabel("Total force (eV/Å)")
plt.xlim(xs[0], xs[-1])
if len(kpoint_mesh) > 1 :
plt.ylim(0.5677, 0.5749)
if smearing == 'fermi-dirac':
plt.xlim(0, 0.0045)
# plt.ylim(0.567, 0.575)
# plt.ylim(0.563, 0.579)
# plt.ylim(0.53, 0.58)
plt.ylim(0.564, 0.578)
if smearing == 'gaussian':
plt.xlim(0, 0.01)
elif property == 'force2':
plt.ylabel("F_x of atom 1 (eV/Å)")
plt.xlim(xs[0], xs[-1])
# plt.semilogx()
# plt.xlabel("smearing temperature (Ry)")
ax1.set_xlabel("smearing temperature (Ry)")
# plt.legend(bbox_to_anchor=(1.21,0.8))
# plt.legend(mode="expand", borderaxespad=0.)
# plt.legend(ncol=2)
# fig.legend(ncol=1)
fig.legend(loc="lower left", bbox_to_anchor=(0.3,0.2))
# plt.legend()
# plt.tight_layout()
# ax1.ticklabel_format(axis="y", style="sci", scilimits=(0,0))
# ax2.ticklabel_format(axis="y", style="", scilimits=(0,0))
fig.tight_layout()
# if property == 'force' :
# plt.savefig('force_smearing_cold.png')
# if property == 'energy' :
# # plt.savefig('energy_smearing.png')
# plt.savefig('energy_smearing_new_range.png')
# elif property == 'stress' :
# plt.savefig('stress_smearing_new.png')
# plt.show()
plt.savefig('force_vs_degauss_FD_GAUS_scale2.565_kpoint30_newzero.png')
params4 = {'font.size': 18,
'font.family': 'serif',
'legend.fontsize': 11,
'mathtext.fontset': 'dejavuserif',
'figure.figsize': (7.5, 4.5) ,
# 'figure.dpi' : 200,
'figure.dpi' : 100,
'legend.handlelength': 2}
smearing_list = [
'gaussian',
'fermi-dirac'
]
group_label = 'force_vs_degauss_gauss_FD2.565'
plot_energy_vs_smearing2(smearing_list, kpoint_mesh= 30, property='energy', scale=2.565)
smearing= gaussian Total number of entries: 13
/home/fdossantos/venv/aiida_git/lib/python3.8/site-packages/aiida/orm/utils/managers.py:103: AiidaDeprecationWarning: dereferencing nodes with links containing double underscores is deprecated, simply replace the double underscores with a single dot instead. For example: `self.inputs.some__label` can be written as `self.inputs.some.label` instead. Support for double underscores will be removed in the future. warnings.warn(
kpoint mesh= [30, 30, 30] num_bands_list= [10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10] smearing= fermi-dirac Total number of entries: 13
/home/fdossantos/venv/aiida_git/lib/python3.8/site-packages/aiida/orm/utils/managers.py:103: AiidaDeprecationWarning: dereferencing nodes with links containing double underscores is deprecated, simply replace the double underscores with a single dot instead. For example: `self.inputs.some__label` can be written as `self.inputs.some.label` instead. Support for double underscores will be removed in the future. warnings.warn(
kpoint mesh= [30, 30, 30] num_bands_list= [10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10]
corrected_energy = -2149.849420861203
#High resolution for the paper
params3 = {'font.size': 18,
'font.family': 'serif',
'legend.fontsize': 12,
'mathtext.fontset': 'dejavuserif',
'figure.figsize': (7, 4.5) ,
'figure.dpi' : 200,
# 'figure.dpi' : 100,
'legend.handlelength': 2}
def plot_energy_vs_smearing(smearing_list, kpoint_mesh, property) :
plt.rcParams.update(params3)
marker_counter = 0
xmin= 0; xmax= 0.1
if property == 'energy' :
# plt.plot( [xmin,xmax], [corrected_energy,corrected_energy] , linestyle=':', color='gray')
pass
elif property in ['stressSSSS']:
plt.plot( [xmin,xmax], [0,0] , linestyle=':', color='gray')
for smearing in smearing_list :
print(' smearing=', smearing)
qb = query_energy_smearing(group_label, smearing)
ys = []
ys2 = []
xs = []
ys_28 = []
xs_28 = []
if isinstance(kpoint_mesh, list) :
for kpoint_mesh_item in kpoint_mesh :
ys = []
ys2 = []
xs = []
ys_28 = []
xs_28 = []
for item in qb.all() :
# print(item)
# calculation = item[1]
chain = item[0]
kpoint = item[1]
if kpoint.attributes['mesh'][0] == kpoint_mesh_item :
mesh = kpoint.attributes['mesh']
calculation_degauss = chain.inputs.pw__parameters['SYSTEM']['degauss']
if property == 'energy' :
if smearing in ['gaussian', 'fermi-dirac' ] :
total_energy = chain.outputs.output_parameters.get_dict()['energy']
total_energy2 = chain.outputs.output_parameters.get_dict()['energy'] - 0.5*chain.outputs.output_parameters.get_dict()['energy_smearing']
ys2.append(total_energy2)
# print( chain.inputs.pw__parameters.attributes )
else :
total_energy = chain.outputs.output_parameters.get_dict()['energy']
ys.append(total_energy)
elif property == 'stress':
stress_tensor = chain.outputs.output_trajectory.get_array('stress')[-1]
# print(stress_tensor)
# print(chain.id)
stress = np.trace( stress_tensor ) / 3.0
# print('stress', stress)
ys.append(stress)
elif property == 'force':
# total force: kind of the sum of the absolute values of all forces components of all atoms
# force = chain.outputs.output_parameters.attributes['total_force'][-1]
force = chain.outputs.output_parameters.attributes['total_force']
ys.append(force)
elif property == 'force2':
# x component of force in atom 0
force = chain.outputs.output_trajectory.get_array('forces')[0][0][0]
ys.append(force)
else :
raise( 'Property not recognized. The available properties are: energy, stress, force.' )
xs.append(calculation_degauss)
if len(xs) > 0:
if smearing in ['gaussian', 'fermi-dirac'] and property == 'energy' :
xs, ys, ys2 = zip(*sorted(zip(xs, ys, ys2)))
else :
xs, ys = zip(*sorted(zip(xs, ys)))
# if smearing in ['cold', 'methfessel-paxton'] :
# xs_28, ys_28 = zip(*sorted(zip(xs_28, ys_28)))
linestyle = '-'
if smearing == 'methfessel-paxton' :
label = ''
label=label+f'${kpoint_mesh_item}^3$'
if len(smearing_list) == 1 and property == 'force':
plt.text(0.013, 0.5741, 'Methfessel-Paxton', fontsize=12)
elif smearing == 'fermi-dirac' :
label = ''
if len(smearing_list) == 1 and property == 'force':
plt.text(0.002925, 0.5741, 'Fermi-Dirac', fontsize=12)
if len(smearing_list) == 1 :
label=label+f'${kpoint_mesh_item}^3$'
elif smearing == 'cold' :
label = ''
if len(smearing_list) == 1 and property == 'force':
plt.text(0.013, 0.5741, 'Cold smearing', fontsize=12)
if len(smearing_list) == 1 :
label=label+f'${kpoint_mesh_item}^3$'
# linestyle = '--'
elif smearing == 'gaussian' :
label = ''
if len(smearing_list) == 1 and property == 'force':
plt.text(0.0065, 0.5741, 'Gaussian', fontsize=12)
if len(smearing_list) == 1 :
label=label+f'${kpoint_mesh_item}^3$'
else :
label = smearing
# label=label+str(kpoint_mesh_item)
plt.plot( xs, ys, marker=marker_list[marker_counter], label=label, linewidth=2, linestyle=linestyle )
if smearing in ['gaussian', 'fermi-dirac'] and property == 'energy':
plt.plot( xs, ys2, color='gray', marker=marker_list[marker_counter], label= 'corrected '+label, linewidth=2)
# if smearing in ['cold', 'methfessel-paxton'] :
# plt.plot( xs_28, ys_28, marker=marker_list[marker_counter], label=label+' 28x28x28', linewidth=2 )
marker_counter += 1
# print( 'kpoint mesh=', mesh )
marker_counter = 0
plt.gca().set_prop_cycle(None)
else : # kpoint_mesh is not a list
for item in qb.all() :
# print(item)
# calculation = item[1]
chain = item[0]
kpoint = item[1]
if kpoint.attributes['mesh'][0] == kpoint_mesh :
mesh = kpoint.attributes['mesh']
calculation_degauss = chain.inputs.pw__parameters['SYSTEM']['degauss']
if property == 'energy' :
if smearing in ['gaussian', 'fermi-dirac' ] :
total_energy = chain.outputs.output_parameters.get_dict()['energy']
total_energy2 = chain.outputs.output_parameters.get_dict()['energy'] - 0.5*chain.outputs.output_parameters.get_dict()['energy_smearing']
ys2.append(total_energy2)
# print( chain.inputs.pw__parameters.attributes )
else :
total_energy = chain.outputs.output_parameters.get_dict()['energy']
ys.append(total_energy)
elif property == 'stress':
stress_tensor = chain.outputs.output_trajectory.get_array('stress')[-1]
# print(stress_tensor)
# print(chain.id)
stress = np.trace( stress_tensor ) / 3.0
# print('stress', stress)
ys.append(stress)
elif property == 'force':
# total force: kind of the sum of the absolute values of all forces components of all atoms
force = chain.outputs.output_parameters.attributes['total_force'][-1]
ys.append(force)
elif property == 'force2':
# x component of force in atom 0
force = chain.outputs.output_trajectory.get_array('forces')[0][0][0]
ys.append(force)
else :
raise( 'Property not recognized. The available properties are: energy, stress, force.' )
xs.append(calculation_degauss)
if kpoint.attributes['mesh'][0] == 28 :
mesh = kpoint.attributes['mesh']
calculation_degauss = chain.inputs.pw__parameters['SYSTEM']['degauss']
if property == 'stress':
stress_tensor = chain.outputs.output_trajectory.get_array('stress')[-1]
# print(stress_tensor)
# print(chain.id)
stress = np.trace( stress_tensor ) / 3.0
# print('stress', stress)
ys_28.append(stress)
xs_28.append(calculation_degauss)
if len(xs) > 0:
if smearing in ['gaussian', 'fermi-dirac'] and property == 'energy' :
xs, ys, ys2 = zip(*sorted(zip(xs, ys, ys2)))
else :
xs, ys = zip(*sorted(zip(xs, ys)))
if smearing in ['cold', 'methfessel-paxton'] :
# print('lens', len(xs_28), len(ys_28))
if kpoint.attributes['mesh'][0] == 28 :
xs_28, ys_28 = zip(*sorted(zip(xs_28, ys_28)))
if smearing == 'methfessel-paxton' :
label = 'Methfessel-Paxton'
elif smearing == 'fermi-dirac' :
label = 'Fermi-Dirac'
elif smearing == 'cold' :
label = 'cold smearing'
elif smearing == 'gaussian' :
label = 'Gaussian'
else :
label = smearing
plt.plot( xs, np.array(ys)+2.149E3+0.85, marker=marker_list[marker_counter], label=label, linewidth=2 )
if smearing in ['gaussian', 'fermi-dirac'] and property == 'energy':
plt.plot( xs, np.array(ys2)+2.149E3+0.85, color='gray', marker=marker_list[marker_counter], label= 'corrected '+label, linewidth=2)
# if smearing in ['cold', 'methfessel-paxton'] :
# plt.plot( xs_28, ys_28, marker=marker_list[marker_counter], label=label+' 28x28x28', linewidth=2 )
marker_counter += 1
print( 'kpoint mesh=', mesh )
deltay = 0.80
if property == 'energy' :
# plt.ylim(corrected_energy-deltay, corrected_energy+deltay*0.1)
# plt.autoscale
plt.xlim(xs[0], 0.02)
# plt.ylim(-2149.8475, -2149.8490)
if not isinstance(kpoint_mesh, list) :
# plt.ylim(-2150.1, -2149.8)
plt.ylim(-2149.92, -2149.84) # new range
plt.ylim(-2149.92+2.149E3+0.85, -2149.84+2.149E3+0.85) # new zero
plt.ylabel("Energy (eV)")
elif property == 'stress':
# plt.ylabel("Tr[ stress tensor ] (GPascal)")
plt.ylabel("Pressure (GPascal)")
plt.xlim(xs[0], 0.02)
plt.ylim(-0.15, 0.2)
# plt.yscale('log')
# plt.ylim(-0.05, 0.5)
elif property == 'force':
plt.ylabel("Total force (eV/Å)")
plt.xlim(xs[0], xs[-1])
if len(kpoint_mesh) > 1 :
plt.ylim(0.5677, 0.5749)
if smearing == 'fermi-dirac':
plt.xlim(0, 0.0045)
# plt.ylim(0.567, 0.575)
# plt.ylim(0.563, 0.579)
# plt.ylim(0.53, 0.58)
plt.ylim(0.564, 0.578)
if smearing == 'gaussian':
plt.xlim(0, 0.01)
elif property == 'force2':
plt.ylabel("F_x of atom 1 (eV/Å)")
plt.xlim(xs[0], xs[-1])
# plt.semilogx()
plt.xlabel("smearing temperature (Ry)")
# plt.legend(bbox_to_anchor=(1.21,0.8))
# plt.legend(mode="expand", borderaxespad=0.)
plt.legend(ncol=1)
# plt.legend()
plt.tight_layout()
if property == 'force' :
plt.savefig('force_smearing_cold.png')
if property == 'energy' :
# plt.savefig('energy_smearing.png')
plt.savefig('energy_smearing_new_range_newzero.png')
elif property == 'stress' :
plt.savefig('stress_smearing_new.png')
plt.show()
smearing_list = [
'gaussian',
'methfessel-paxton',
'cold',
'fermi-dirac'
]
group_label = 'energy_vs_degauss'
plot_energy_vs_smearing(smearing_list, kpoint_mesh= 14, property='energy')
smearing= gaussian Total number of entries: 30
/home/fdossantos/venv/aiida_git/lib/python3.8/site-packages/aiida/orm/utils/managers.py:103: AiidaDeprecationWarning: dereferencing nodes with links containing double underscores is deprecated, simply replace the double underscores with a single dot instead. For example: `self.inputs.some__label` can be written as `self.inputs.some.label` instead. Support for double underscores will be removed in the future. warnings.warn(
kpoint mesh= [14, 14, 14] smearing= methfessel-paxton Total number of entries: 40
/home/fdossantos/venv/aiida_git/lib/python3.8/site-packages/aiida/orm/utils/managers.py:103: AiidaDeprecationWarning: dereferencing nodes with links containing double underscores is deprecated, simply replace the double underscores with a single dot instead. For example: `self.inputs.some__label` can be written as `self.inputs.some.label` instead. Support for double underscores will be removed in the future. warnings.warn(
kpoint mesh= [14, 14, 14] smearing= cold Total number of entries: 40
/home/fdossantos/venv/aiida_git/lib/python3.8/site-packages/aiida/orm/utils/managers.py:103: AiidaDeprecationWarning: dereferencing nodes with links containing double underscores is deprecated, simply replace the double underscores with a single dot instead. For example: `self.inputs.some__label` can be written as `self.inputs.some.label` instead. Support for double underscores will be removed in the future. warnings.warn(
kpoint mesh= [14, 14, 14] smearing= fermi-dirac Total number of entries: 30
/home/fdossantos/venv/aiida_git/lib/python3.8/site-packages/aiida/orm/utils/managers.py:103: AiidaDeprecationWarning: dereferencing nodes with links containing double underscores is deprecated, simply replace the double underscores with a single dot instead. For example: `self.inputs.some__label` can be written as `self.inputs.some.label` instead. Support for double underscores will be removed in the future. warnings.warn(
kpoint mesh= [14, 14, 14]
smearing_list = [
# 'gaussian',
# 'methfessel-paxton',
'cold',
# 'fermi-dirac'
]
group_label = 'force_vs_degauss_nicola'
plot_energy_vs_smearing(smearing_list, kpoint_mesh= [30,36,42,48,54,60,70,80,90,100], property='force')
# plot_energy_vs_smearing(smearing_list, kpoint_mesh= [12,18,24], property='force')
smearing= cold Total number of entries: 169
plot_energy_vs_smearing(smearing_list, kpoint_mesh= 14, property='stress')
smearing= gaussian Total number of entries: 30 kpoint mesh= [14, 14, 14] smearing= methfessel-paxton Total number of entries: 40 kpoint mesh= [14, 14, 14] smearing= cold Total number of entries: 40 kpoint mesh= [14, 14, 14] smearing= fermi-dirac Total number of entries: 30 kpoint mesh= [14, 14, 14]
