scnc
/
abipy
mirror of https://github.com/abinit/abipy.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #6 from abinit/gp_qha 

Use kptopt 2 for strain inputs
This commit is contained in:
Guido Petretto 2018-07-31 15:17:44 +02:00 committed by GitHub
parent 548d8470ef dab7ccf7a4
commit b96f42faca
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 194 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
abipy/abio/inputs.py
View File

@ -1132,7 +1132,7 @@ class AbinitInput(six.with_metaclass(abc.ABCMeta, AbiAbstractInput, MSONable, Ha
#def pop_relax_vars(self):
# """
# Remove all the `ird*` variables present in self.
# Remove all the relax variables present in self.
# Return dictionary with the variables that have been removed.
# """
# return scf_input.pop_vars(["ionmov", "optcell", "ntime", "dilatmx"])
@ -1425,7 +1425,7 @@ class AbinitInput(six.with_metaclass(abc.ABCMeta, AbiAbstractInput, MSONable, Ha
return multi
def make_strain_perts_inputs(self, tolerance=None, manager=None, phonon_pert=True, kptopt=3):
def make_strain_perts_inputs(self, tolerance=None, manager=None, phonon_pert=True, kptopt=2):
"""
Return inputs for the strain perturbation calculation.
This functions should be called with an input that represents a GS run.
@ -1435,14 +1435,16 @@ class AbinitInput(six.with_metaclass(abc.ABCMeta, AbiAbstractInput, MSONable, Ha
Defaults to {"tolvrs": 1.0e-12}.
manager: |TaskManager| of the task. If None, the manager is initialized from the config file.
phonon_pert: is True also the phonon perturbations will be considered. Default False.
kptopt: 2 to take into account time-reversal symmetry. Default 3 for backward compatibility.
kptopt: 2 to take into account time-reversal symmetry.
"""
if tolerance is None:
tolerance = {"tolvrs": 1.0e-12}
if len(tolerance) != 1 or any(k not in _TOLVARS for k in tolerance):
raise self.Error("Invalid tolerance: {}".format(str(tolerance)))
perts = self.abiget_irred_strainperts(kptopt=2, manager=manager, phonon_pert=phonon_pert)
perts = self.abiget_irred_strainperts(kptopt=kptopt, manager=manager, phonon_pert=phonon_pert)
#print("Stress perts:", perts)
# Build list of datasets (one input per perturbation)
multi = MultiDataset.replicate_input(input=self, ndtset=len(perts))
@ -1450,6 +1452,7 @@ class AbinitInput(six.with_metaclass(abc.ABCMeta, AbiAbstractInput, MSONable, Ha
for pert, inp in zip(perts, multi):
rfdir = 3 * [0]
rfdir[pert.idir -1] = 1
if pert.ipert <= len(self.structure):
inp.set_vars(rfphon=1, # Activate the calculation of the atomic dispacement perturbations
rfatpol=[pert.ipert, pert.ipert],
@ -1457,34 +1460,32 @@ class AbinitInput(six.with_metaclass(abc.ABCMeta, AbiAbstractInput, MSONable, Ha
nqpt=1, # One wavevector is to be considered
qpt=(0, 0, 0), # q-wavevector.
kptopt=kptopt, # No symmetries
iscf=7,
#iscf=7,
paral_kgb=0
)
elif pert.ipert == len(self.structure) + 3:
inp.set_vars(rfstrs=1, # Activate the calculation of the strain perturbations (uniaxial)
rfdir=rfdir,
nqpt=1, # One wavevector is to be considered
qpt=(0, 0, 0), # q-wavevector.
kptopt=kptopt, # No symmetries
iscf=7,
#iscf=7,
paral_kgb=0
)
elif pert.ipert == len(self.structure) + 4:
inp.set_vars(rfstrs=2, # Activate the calculation of the strain perturbations (shear)
rfdir=rfdir,
nqpt=1, # One wavevector is to be considered
qpt=(0, 0, 0), # q-wavevector.
kptopt=kptopt, # No symmetries
iscf=7,
#iscf=7,
paral_kgb=0
)
inp.pop_tolerances()
inp.set_vars(tolerance)
# Adding buffer to help convergence ...
if 'nbdbuf' not in inp:
nbdbuf = max(int(0.1*inp['nband']), 4)
inp.set_vars(nband=inp['nband']+nbdbuf, nbdbuf=nbdbuf)
return multi

2
abipy/core/structure.py
View File

@ -1274,6 +1274,8 @@ class Structure(pymatgen.Structure, NotebookWriter):
"""
if fmt in ("abivars", "abinit"):
return self.abi_string
elif fmt == "abipython":
return pformat(self.to_abivars(), indent=4)
elif fmt == "qe":
from pymatgen.io.pwscf import PWInput
return str(PWInput(self, pseudo={s: s + ".pseudo" for s in self.symbol_set}))

119
abipy/dfpt/ddb.py
View File

@ -152,8 +152,12 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
app("")
app("Number of q-points in DDB: %d" % len(self.qpoints))
app("guessed_ngqpt: %s (guess for the q-mesh divisions made by AbiPy)" % self.guessed_ngqpt)
app("Has electric-field perturbation: %s" % self.has_emacro_terms())
app("Has Born effective charges: %s" % self.has_bec_terms())
app("Has (at least one) atomic pertubation: %s" % self.has_at_least_one_atomic_perturbation())
app("Has (at least one) electric-field perturbation: %s" % self.has_emacro_terms(select="at_least_one"))
app("Has (at least one) Born effective charge: %s" % self.has_bec_terms(select="at_least_one"))
app("Has (all) strain terms: %s" % self.has_strain_terms(select="all"))
app("Has (all) internal strain terms: %s" % self.has_internalstrain_terms(select="all"))
app("Has (all) piezoelectric terms: %s" % self.has_piezoelectric_terms(select="all"))
if verbose:
app(self.qpoints.to_string(verbose=verbose, title="Q-points in DDB"))
@ -571,6 +575,23 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
"""
return self.has_emacro_terms(select=select) and self.has_bec_terms(select=select)
@lru_cache(typed=True)
def has_at_least_one_atomic_perturbation(self):
"""
True if the DDB file contains info on (at least one) atomic perturbation.
"""
natom = len(self.structure)
ap_list = list(itertools.product(range(1, 4), range(1, natom + 1)))
for qpt, df in self.computed_dynmat.items():
index_set = set(df.index)
for p1 in ap_list:
for p2 in ap_list:
p12 = p1 + p2
if p12 in index_set: return True
return False
@lru_cache(typed=True)
def has_emacro_terms(self, select="at_least_one"):
"""
@ -600,7 +621,7 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
else:
raise ValueError("Wrong select %s" % str(select))
return False
return False if select == "at_least_one" else True
@lru_cache(typed=True)
def has_bec_terms(self, select="at_least_one"):
@ -631,18 +652,24 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
else:
raise ValueError("Wrong select %s" % str(select))
return False
return False if select == "at_least_one" else True
@lru_cache(typed=True)
def has_strain_terms(self, select="at_least_one"):
def has_strain_terms(self, select="all"):
"""
True if the DDB file contains info on the strain perturbation.
True if the DDB file contains info on the (clamped-ion) strain perturbation
(i.e. 2nd order derivatives wrt strain)
Args:
select: Possible values in ["at_least_one", "all"]
If select == "at_least_one", we check if there's at least one entry associated to the strain.
and we assume that anaddb will be able to reconstruct the full tensor by symmetry.
If select == "all", all tensor components must be present in the DDB file.
.. note::
As anaddb is not yet able to reconstruct the strain terms by symmetry,
the default value for select is "all"
"""
gamma = Kpoint.gamma(self.structure.reciprocal_lattice)
if gamma not in self.computed_dynmat:
@ -654,6 +681,46 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
sp_list = list(itertools.product(range(1, 4), [natom + 3, natom + 4]))
for p1 in sp_list:
for p2 in sp_list:
p12 = p1 + p2
if select == "at_least_one":
if p12 in index_set: return True
elif select == "all":
if p12 not in index_set:
#print("p12", p12, "not in index_set")
return False
else:
raise ValueError("Wrong select %s" % str(select))
return False if select == "at_least_one" else True
@lru_cache(typed=True)
def has_internalstrain_terms(self, select="all"):
"""
True if the DDB file contains internal strain terms
i.e "off-diagonal" 2nd order derivatives wrt (strain, atomic displacement)
Args:
select: Possible values in ["at_least_one", "all"]
If select == "at_least_one", we check if there's at least one entry associated to the strain.
and we assume that anaddb will be able to reconstruct the full tensor by symmetry.
If select == "all", all tensor components must be present in the DDB file.
.. note::
As anaddb is not yet able to reconstruct the strain terms by symmetry,
the default value for select is "all"
"""
gamma = Kpoint.gamma(self.structure.reciprocal_lattice)
if gamma not in self.computed_dynmat:
return False
index_set = set(self.computed_dynmat[gamma].index)
natom = len(self.structure)
sp_list = list(itertools.product(range(1, 4), [natom + 3, natom + 4]))
ap_list = list(itertools.product(range(1, 4), range(1, natom + 1)))
for p1 in sp_list:
for p2 in ap_list:
p12 = p1 + p2
if select == "at_least_one":
if p12 in index_set: return True
@ -662,7 +729,45 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
else:
raise ValueError("Wrong select %s" % str(select))
return False
return False if select == "at_least_one" else True
@lru_cache(typed=True)
def has_piezoelectric_terms(self, select="all"):
"""
True if the DDB file contains piezoelectric terms
i.e "off-diagonal" 2nd order derivatives wrt (electric_field, strain)
Args:
select: Possible values in ["at_least_one", "all"]
If select == "at_least_one", we check if there's at least one entry associated to the strain.
and we assume that anaddb will be able to reconstruct the full tensor by symmetry.
If select == "all", all tensor components must be present in the DDB file.
.. note::
As anaddb is not yet able to reconstruct the (strain, electric) terms by symmetry,
the default value for select is "all"
"""
gamma = Kpoint.gamma(self.structure.reciprocal_lattice)
if gamma not in self.computed_dynmat:
return False
index_set = set(self.computed_dynmat[gamma].index)
natom = len(self.structure)
sp_list = list(itertools.product(range(1, 4), [natom + 3, natom + 4]))
ep_list = list(itertools.product(range(1, 4), [natom + 2]))
for p1 in sp_list:
for p2 in ep_list:
p12 = p1 + p2
if select == "at_least_one":
if p12 in index_set: return True
elif select == "all":
if p12 not in index_set: return False
else:
raise ValueError("Wrong select %s" % str(select))
return False if select == "at_least_one" else True
def view_phononwebsite(self, browser=None, verbose=0, dryrun=False, **kwargs):
"""

15
abipy/dfpt/tests/test_ddb.py
View File

@ -129,6 +129,10 @@ class DdbTest(AbipyTest):
assert not ddb.has_bec_terms(select="all")
assert not ddb.has_emacro_terms()
assert not ddb.has_lo_to_data()
assert not ddb.has_internalstrain_terms()
assert not ddb.has_piezoelectric_terms()
assert not ddb.has_strain_terms()
assert ddb.has_at_least_one_atomic_perturbation()
ref_qpoints = np.reshape([
0.00000000E+00, 0.00000000E+00, 0.00000000E+00,
@ -222,8 +226,8 @@ class DdbTest(AbipyTest):
# Test Lru_cache as well
assert ddb.has_bec_terms(select="at_least_one")
assert ddb.has_bec_terms(select="at_least_one")
assert not ddb.has_bec_terms(select="all")
assert not ddb.has_bec_terms(select="all")
assert ddb.has_bec_terms(select="all")
assert ddb.has_bec_terms(select="all")
assert ddb.has_emacro_terms()
assert ddb.has_lo_to_data()
@ -354,8 +358,11 @@ class DdbTest(AbipyTest):
self.skip_if_abinit_not_ge("8.9.3")
with abilab.abiopen(abidata.ref_file("refs/alas_elastic_dfpt/AlAs_elastic_DDB")) as ddb:
assert ddb.has_strain_terms()
assert not ddb.has_strain_terms("all")
assert ddb.has_strain_terms(select="at_least_one")
assert ddb.has_strain_terms(select="all")
#assert ddb.has_internalstrain_terms(select="all")
#assert ddb.has_piezoelectric_terms_terms(select="all")
#assert ddb.has_at_least_one_atomic_perturbation()
e = ddb.anaget_elastic(has_dde=True, has_gamma_ph=True, verbose=2)
self.assert_almost_equal(e.elastic_relaxed[0,0,0,0], 120.41874336082199)
self.assert_almost_equal(e.piezo_relaxed[0,1,2], -0.030391022487094244)

58
abipy/flowtk/works.py
View File

@ -3,4 +3,60 @@ from __future__ import print_function, division, unicode_literals, absolute_impo
import numpy as np
from pymatgen.io.abinit.works import Work
from pymatgen.io.abinit.works import Work, MergeDdb
from abipy.abio import input_tags as tags
class ElasticWork(Work, MergeDdb):
"""
Work for the calculation of elastic constants and (optionally piezoelectric tensor
"""
@classmethod
def from_scf_input(cls, scf_input, tolerance=None, with_internal_strain=True, with_piezoelectric=False, manager=None):
"""
Similar to `from_scf_task`, the difference is that this method requires
an input for SCF calculation instead of a ScfTask. All the tasks (Scf + Phonon)
are packed in a single Work whereas in the previous case we usually have multiple works.
"""
new = cls(manager=manager)
# Register task for SCF calculation.
scf_task = new.register_scf_task(scf_input)
multi = scf_task.input.make_strain_perts_inputs(tolerance=tolerance, manager=manager, phonon_pert=False, kptopt=2)
ddk_tasks = []
if with_piezoelectric:
#sfc_task.input.make_ddk_inputs(tolerance=None, manager=None):
for inp in multi.filter_by_tags(tags=tags.DDK):
ddk_task = new.register_ddk_task(inp, deps={scf_task: "WFK"})
ddk_tasks.append(ddk_task)
assert len(ddk_tasks) == 3
if with_internal_strain:
ph_deps = {scf_task: "WFK"}
#if with_piezoelectric: ph_deps.update()
for inp in multi.filter_by_tags(tags=tags.PHONON):
new.register_phonon_task(inp, deps=ph_deps)
#bec_deps = {ddk_task: "DDK" for ddk_task in ddk_tasks}
elast_deps = {scf_task: "WFK"}
#if with_piezoelectric: ph_deps.update()
#for inp in multi.filter_by_tags(tags=tags.STRAIN)
for inp in multi:
new.register_elastic_task(inp, deps=elast_deps)
return new
def on_all_ok(self):
"""
This method is called when all the tasks reach S_OK.
Ir runs `mrgddb` in sequential on the local machine to produce
the final DDB file in the outdir of the `Work`.
"""
# Merge DDB files.
out_ddb = self.merge_ddb_files(delete_source_ddbs=False)
results = self.Results(node=self, returncode=0, message="DDB merge done")
return results