"""Module to run Tc calculation."""
# Ref:
# https://www.nature.com/articles/s41524-022-00933-1
from jarvis.io.qe.outputs import DataFileSchema
from jarvis.core.atoms import Atoms
from jarvis.core.kpoints import Kpoints3D
from jarvis.tasks.qe.qe import QEjob
from jarvis.tasks.qe.converg import converg_kpoints
import numpy as np
import os
from jarvis.core.utils import get_factors
# from jarvis.core.utils import get_factors
[docs]def calc_Tc(wlog=300, lamb=1.0, mu=0.1):
"""Calculate Tc."""
tc = (wlog / 1.2) * np.exp(
-1.04 * (1 + lamb) / (lamb * (1 - 0.62 * mu) - mu)
)
return tc
[docs]def parse_lambda(filename="lambda"):
"""Parse lambda file."""
f = open(filename, "r")
lines = f.read().splitlines()
f.close()
for i in lines:
if "Broadening" in i:
tmp = i.split()
print(i.split())
wlog = float(tmp[-1])
lamb = float(tmp[3])
Tc = calc_Tc(wlog=wlog, lamb=lamb)
print("Tc", Tc)
print()
[docs]def non_prime_kpoints(kpts=[]):
"""Get non prime kpoints."""
mem = []
for i in kpts:
facts = get_factors(i)
if len(facts) == 1:
val = i + 1
else:
val = i
mem.append(val)
return mem
[docs]def very_clean():
"""Clean files."""
cmd = (
"rm -r elph_dir _ph0 *.dos* *.in *.json *.fc *.freq* *.save "
+ "*.out *.dyn* *wfc* *.xml *save lambda *.modes dyn*"
)
os.system(cmd)
[docs]class SuperCond(object):
"""Module to calculate Tc."""
def __init__(
self,
atoms=None,
kp=None,
qp=None,
qe_cmd="pw.x",
relax_calc="'vc-relax'",
pressure=None,
psp_dir=None,
url=None,
psp_temp_name=None,
clean_files=True,
):
"""Initialize the class."""
self.atoms = atoms
self.kp = kp
self.qp = qp
# self.kp = kp
# self.qp = qp
self.relax_calc = relax_calc
self.qe_cmd = qe_cmd
self.psp_dir = psp_dir
self.url = url
self.pressure = pressure
self.psp_temp_name = psp_temp_name
self.clean_files = clean_files
[docs] def to_dict(self):
"""Get dictionary."""
info = {}
info["atoms"] = self.atoms.to_dict()
info["kp"] = self.kp.to_dict()
info["qp"] = self.qp.to_dict()
info["qe_cmd"] = self.qe_cmd
info["psp_dir"] = self.psp_dir
info["psp_temp_name"] = self.psp_temp_name
info["url"] = self.url
info["relax_calc"] = self.relax_calc
info["pressure"] = self.pressure
return info
[docs] @classmethod
def from_dict(self, info={}):
"""Load from a dictionary."""
return SuperCond(
atoms=Atoms.from_dict(info["atoms"]),
kp=Kpoints3D.from_dict(info["kp"]),
qp=Kpoints3D.from_dict(info["qp"]),
qe_cmd=info["qe_cmd"],
pressure=info["pressure"],
relax_calc=info["relax_calc"],
psp_dir=info["psp_dir"],
url=info["url"],
psp_temp_name=info["psp_temp_name"],
)
[docs] def runjob(self):
"""Calculate Tc using QE."""
# Still under development,
atoms = self.atoms
kp = self.kp
qp = self.qp
if not kp._kpoints:
kp_len = converg_kpoints(
atoms=atoms,
qe_cmd=self.qe_cmd,
psp_dir=self.psp_dir,
url=self.url,
psp_temp_name=self.psp_temp_name,
)
kp = Kpoints3D().automatic_length_mesh(
lattice_mat=atoms.lattice_mat, length=kp_len
)
kpts = kp._kpoints[0]
kpts = non_prime_kpoints(kpts)
kp = Kpoints3D(kpoints=[kpts])
print("kpts", kpts)
else:
kpts = kp._kpoints[0]
kpts = non_prime_kpoints(kpts)
print("kpts xyz", kpts)
nq1 = get_factors(kpts[0])[0]
nq2 = get_factors(kpts[1])[0]
nq3 = get_factors(kpts[2])[0]
qp = Kpoints3D(kpoints=[[nq1, nq2, nq3]])
self.kp = kp
self.qp = qp
relax = {
"control": {
# "calculation": "'scf'",
"calculation": self.relax_calc, # "'vc-relax'",
"restart_mode": "'from_scratch'",
"prefix": "'RELAX'",
"outdir": "'./'",
"tstress": ".true.",
"tprnfor": ".true.",
"disk_io": "'nowf'",
"wf_collect": ".true.",
"pseudo_dir": None,
"verbosity": "'high'",
"nstep": 100,
},
"system": {
"ibrav": 0,
"nat": None,
"ntyp": None,
"ecutwfc": 45,
"ecutrho": 250,
"q2sigma": 1,
"ecfixed": 44.5,
"qcutz": 800,
"occupations": "'smearing'",
"degauss": 0.01,
"lda_plus_u": ".false.",
},
"electrons": {
"diagonalization": "'david'",
"mixing_mode": "'local-TF'",
"mixing_beta": 0.3,
"conv_thr": "1d-9",
},
"ions": {"ion_dynamics": "'bfgs'"},
"cell": {"cell_dynamics": "'bfgs'", "cell_dofree": "'all'"},
}
if self.pressure is not None:
relax["cell"]["press"] = self.pressure
qejob_relax = QEjob(
atoms=atoms,
input_params=relax,
output_file="relax.out",
qe_cmd=self.qe_cmd,
jobname="relax",
kpoints=kp,
input_file="arelax.in",
url=self.url,
psp_dir=self.psp_dir,
psp_temp_name=self.psp_temp_name,
)
info = qejob_relax.runjob()
xml_path = info["xml_path"]
xml_dat = DataFileSchema(xml_path)
final_strt = xml_dat.final_structure
print("final_strt", final_strt)
# scf_init =
scf_init = {
"control": {
"calculation": "'scf'",
"restart_mode": "'from_scratch'",
"prefix": "'QE'",
"outdir": "'./'",
"tstress": ".true.",
"tprnfor": ".true.",
"disk_io": "'nowf'",
"pseudo_dir": None,
"verbosity": "'high'",
"nstep": 100,
"etot_conv_thr": "1.0d-5",
"forc_conv_thr": "1.0d-4",
},
"system": {
"ibrav": 0,
"degauss": 0.01,
"nat": None,
"ntyp": None,
"ecutwfc": 45,
"ecutrho": 250,
"occupations": "'smearing'",
"smearing": "'mp'",
"la2F ": ".true.",
},
"electrons": {
"diagonalization": "'david'",
"mixing_mode": "'plain'",
"mixing_beta": 0.7,
"conv_thr": "1d-9",
},
}
qejob_scf_init = QEjob(
atoms=final_strt,
input_params=scf_init,
output_file="scf_init.out",
qe_cmd=self.qe_cmd,
jobname="scf_init",
kpoints=kp,
input_file="ascf_init.in",
url=self.url,
psp_temp_name=self.psp_temp_name,
psp_dir=self.psp_dir,
)
info_scf = qejob_scf_init.runjob()
print(info_scf)
# kpts = kp._kpoints[0]
qpts = qp._kpoints[0]
nq1 = qpts[0] # get_factors(kpts[0])[0]
nq2 = qpts[1] # get_factors(kpts[1])[0]
nq3 = qpts[2] # get_factors(kpts[2])[0]
ph = {
"inputph": {
"prefix": "'QE'",
"fildyn": "'QE.dyn'",
"outdir": "'./'",
"ldisp": ".true.",
"trans": ".true.",
"fildvscf": "'dvscf'",
"electron_phonon": "'interpolated'",
"el_ph_sigma": 0.005,
"nq1": nq1,
"nq2": nq2,
"nq3": nq3,
"tr2_ph": "1.0d-12",
}
}
qejob_ph = QEjob(
atoms=final_strt,
input_params=ph,
output_file="ph.out",
qe_cmd=self.qe_cmd.replace("pw.x", "ph.x"),
jobname="ph",
kpoints=None,
input_file="aph.in",
url=self.url,
psp_temp_name=self.psp_temp_name,
psp_dir=self.psp_dir,
)
qejob_ph.runjob()
# import sys
# sys.exit()
qr = {
"input": {
"zasr": "'simple'",
"fildyn": "'QE.dyn'",
"flfrc": "'QE333.fc'",
"la2F": ".true.",
}
}
qejob_qr = QEjob(
atoms=final_strt,
input_params=qr,
output_file="q2r.out",
# qe_cmd="/home/knc6/Software/qe/q-e/bin/q2r.x",
qe_cmd=self.qe_cmd.replace("pw.x", "q2r.x"),
jobname="qr",
kpoints=None,
input_file="aqr.in",
url=self.url,
psp_temp_name=self.psp_temp_name,
psp_dir=self.psp_dir,
)
qejob_qr.runjob()
kpts = kp._kpoints[0]
matdyn = {
"input": {
"asr": "'simple'",
"flfrc": "'QE333.fc'",
"flfrq": "'QE333.freq'",
"la2F": ".true.",
"dos": ".true.",
"fldos": "'phonon.dos'",
"nk1": kpts[0],
"nk2": kpts[1],
"nk3": kpts[2],
"ndos": 50,
}
}
qejob_matdyn = QEjob(
atoms=final_strt,
input_params=matdyn,
output_file="matdyn.out",
# qe_cmd="/home/knc6/Software/qe/q-e/bin/matdyn.x",
qe_cmd=self.qe_cmd.replace("pw.x", "matdyn.x"),
jobname="matdyn",
kpoints=None,
input_file="amatdyn.in",
url=self.url,
psp_temp_name=self.psp_temp_name,
psp_dir=self.psp_dir,
)
qejob_matdyn.runjob()
parse_lambda()
cmd = "rm -r _ph0 *.wfc* */wfc* */*UPF */*upf */charge-density.dat"
if self.clean_files:
os.system(cmd)