LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-22 20:34:38 +01:00
Code Issues Releases Wiki Activity

Symmetric DFT implementation

Browse Source
This commit is contained in:
Dario Fiore Mosca 2024-11-12 15:22:25 +01:00
parent 44d791d90d
commit 7e99a2f291
9 changed files with 105 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.idea/.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
# Default ignored files
/shelf/
/workspace.xml

12
.idea/dft_tools.src.iml Normal file
View File

@ -0,0 +1,12 @@
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="inheritedJdk" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
<component name="PyDocumentationSettings">
<option name="format" value="PLAIN" />
<option name="myDocStringFormat" value="Plain" />
</component>
</module>

18
.idea/inspectionProfiles/Project_Default.xml Normal file
View File

@ -0,0 +1,18 @@
<component name="InspectionProjectProfileManager">
<profile version="1.0">
<option name="myName" value="Project Default" />
<inspection_tool class="PyInterpreterInspection" enabled="false" level="WARNING" enabled_by_default="false" />
<inspection_tool class="PyPackageRequirementsInspection" enabled="true" level="WARNING" enabled_by_default="true">
<option name="ignoredPackages">
<value>
<list size="4">
<item index="0" class="java.lang.String" itemvalue="python3" />
<item index="1" class="java.lang.String" itemvalue="triqs.3.1.x" />
<item index="2" class="java.lang.String" itemvalue="lapack" />
<item index="3" class="java.lang.String" itemvalue="f90wrap" />
</list>
</value>
</option>
</inspection_tool>
</profile>
</component>

6
.idea/inspectionProfiles/profiles_settings.xml Normal file
View File

@ -0,0 +1,6 @@
<component name="InspectionProjectProfileManager">
<settings>
<option name="USE_PROJECT_PROFILE" value="false" />
<version value="1.0" />
</settings>
</component>

8
.idea/modules.xml Normal file
View File

@ -0,0 +1,8 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/dft_tools.src.iml" filepath="$PROJECT_DIR$/.idea/dft_tools.src.iml" />
</modules>
</component>
</project>

6
.idea/vcs.xml Normal file
View File

@ -0,0 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="VcsDirectoryMappings">
<mapping directory="$PROJECT_DIR$" vcs="Git" />
</component>
</project>

42
python/triqs_dft_tools/converters/plovasp/sc_dmft.py
View File

@ -72,22 +72,26 @@ def is_vasp_running(vasp_pid):
pid_exists = mpi.bcast(pid_exists) pid_exists = mpi.bcast(pid_exists)
return pid_exists return pid_exists
def get_dft_energy(): def get_dft_energy():
""" """
Reads energy from the last line of OSZICAR. Reads DFT energy from the last line of Vasp's OSZICAR or from vasptriqs.h5
""" """
with open('OSZICAR', 'r') as f: h5_energy = False
nextline = f.readline() if os.path.isfile('vaspout.h5'):
while nextline.strip(): with HDFArchive('vaspout.h5', 'r') as h5:
line = nextline if 'oszicar' in h5['intermediate/ion_dynamics']:
nextline = f.readline() dft_energy = h5['intermediate/ion_dynamics/oszicar'][-1,1]
# print "OSZICAR: ", line[:-1] h5_energy = True
try: # as backup use OSZICAR file
if not h5_energy:
with open('OSZICAR', 'r') as file:
nextline = file.readline()
while nextline.strip():
line = nextline
nextline = file.readline()
dft_energy = float(line.split()[2]) dft_energy = float(line.split()[2])
except ValueError:
print("Cannot read energy from OSZICAR, setting it to zero")
dft_energy = 0.0
return dft_energy return dft_energy
@ -98,9 +102,8 @@ class bcolors:
YELLOW = '\033[93m' YELLOW = '\033[93m'
RED = '\033[91m' RED = '\033[91m'
ENDC = '\033[0m' ENDC = '\033[0m'
#
# Main self-consistent cycle # Main self-consistent cycle
#
def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version): def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
""" """
""" """
@ -117,14 +120,13 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
mpi.barrier() mpi.barrier()
while is_vasp_lock_present(): while is_vasp_lock_present():
time.sleep(1) time.sleep(1)
# if debug: print bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC if debug: print(bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC)
if not is_vasp_running(vasp_pid): if not is_vasp_running(vasp_pid):
mpi.report(" VASP stopped") mpi.report(" VASP stopped")
vasp_running = False vasp_running = False
break break
# Tell VASP to stop if the maximum number of iterations is reached # Tell VASP to stop if the maximum number of iterations is reached
if debug: print(bcolors.MAGENTA + "rank %s"%(mpi.rank) + bcolors.ENDC) if debug: print(bcolors.MAGENTA + "rank %s"%(mpi.rank) + bcolors.ENDC)
err = 0 err = 0
@ -161,7 +163,7 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
# electron.F around line 644 # electron.F around line 644
iter_dft = 0 iter_dft = 0
if vasp_version == 'standard': if vasp_version == 'standard' or vasp_version == 'ncl':
copyfile(src='GAMMA',dst='GAMMA_recent') copyfile(src='GAMMA',dst='GAMMA_recent')
while iter_dft < n_iter_dft: while iter_dft < n_iter_dft:
# insert recalculation of GAMMA here # insert recalculation of GAMMA here
@ -190,7 +192,7 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
vasp_running = False vasp_running = False
break break
iter_dft += 1 iter_dft += 1
if vasp_version == 'standard': if vasp_version == 'standard' or vasp_version == 'ncl':
copyfile(src='GAMMA_recent',dst='GAMMA') copyfile(src='GAMMA_recent',dst='GAMMA')
iter += 1 iter += 1
if iter == n_iter: if iter == n_iter:
@ -253,8 +255,8 @@ def main():
except KeyError: except KeyError:
vasp_version = 'standard' vasp_version = 'standard'
if vasp_version != 'standard' and vasp_version != 'no_gamma_write': #if vasp_version != 'standard' and vasp_version != 'no_gamma_write':
raise Exception('vasp_version has to be standard or no_gamma_write') # raise Exception('vasp_version has to be standard or no_gamma_write')
# if len(sys.argv) > 1: # if len(sys.argv) > 1:
# vasp_path = sys.argv[1] # vasp_path = sys.argv[1]

18
python/triqs_dft_tools/converters/plovasp/vaspio.py
View File

@ -82,7 +82,7 @@ class VaspData:
self.plocar = h5Plocar(h5path) self.plocar = h5Plocar(h5path)
self.poscar = h5Poscar(h5path) self.poscar = h5Poscar(h5path)
self.kpoints = h5Kpoints(h5path) self.kpoints = h5Kpoints(h5path)
self.eigenval = h5Eigenval(h5path) self.eigenval = h5Eigenval(h5path, self.kpoints.ksymmap)
self.doscar = h5Doscar(h5path) self.doscar = h5Doscar(h5path)
else: else:
self.plocar = Plocar() self.plocar = Plocar()
@ -716,28 +716,34 @@ class h5Kpoints:
# h5path = './vasptriqs.h5' # h5path = './vasptriqs.h5'
with HDFArchive(h5path, 'a') as archive: with HDFArchive(h5path, 'a') as archive:
kpoints = archive['results/electron_eigenvalues'] kpoints = archive['results/electron_eigenvalues']
self.nktot = kpoints['kpoints'] self.nkred = kpoints['kpoints']
self.kpts = kpoints['kpoint_coords'] self.kpts = kpoints['kpoint_coords_full']
self.kwghts = kpoints['kpoints_symmetry_weight'] self.nktot = len(self.kpts)
self.kwghts = kpoints['kpoints_symmetry_weight_full']
self.ksymmap = kpoints['kpoints_symmetry_mapping']
self.ksymmap -= 1
try: try:
self.ntet = kpoints['num_tetrahedra'] self.ntet = kpoints['num_tetrahedra']
self.vtet = kpoints['volume_weight_tetrahedra'] self.vtet = kpoints['volume_weight_tetrahedra']
self.itet = kpoints['coordinate_id_tetrahedra'] self.itet = kpoints['coordinate_id_tetrahedra']
except KeyError: except KeyError:
print(" No tetrahedron data found in vasptriqs.h5. Skipping...") print(" No tetrahedron data found in vaspout.h5. Skipping...")
self.ntet = 0 self.ntet = 0
print() print()
print(" {0:>26} {1:d}".format("Reduced number of k-points:", self.nkred))
print(" {0:>26} {1:d}".format("Total number of k-points:", self.nktot)) print(" {0:>26} {1:d}".format("Total number of k-points:", self.nktot))
print(" {0:>26} {1:d}".format("Total number of tetrahedra:", self.ntet)) print(" {0:>26} {1:d}".format("Total number of tetrahedra:", self.ntet))
class h5Eigenval: class h5Eigenval:
def __init__(self, h5path): def __init__(self, h5path, symmap):
with HDFArchive(h5path, 'a') as archive: with HDFArchive(h5path, 'a') as archive:
self.eigs = archive['results/electron_eigenvalues']['eigenvalues'] self.eigs = archive['results/electron_eigenvalues']['eigenvalues']
self.eigs = self.eigs[:, symmap, :]
self.ferw = archive['results/electron_eigenvalues']['fermiweights'] self.ferw = archive['results/electron_eigenvalues']['fermiweights']
self.ferw = self.ferw[:, symmap, :]
# TODO Change the format in VASP to have [kpoints, bands, spin] # TODO Change the format in VASP to have [kpoints, bands, spin]
self.eigs = np.transpose(self.eigs, (1, 2, 0)) self.eigs = np.transpose(self.eigs, (1, 2, 0))
self.ferw = np.transpose(self.ferw, (1, 2, 0)) self.ferw = np.transpose(self.ferw, (1, 2, 0))

29
python/triqs_dft_tools/sumk_dft.py
View File

@ -578,7 +578,8 @@ class SumkDFT(object):
mesh_values = self.mesh_values mesh_values = self.mesh_values
elif not mesh is None: elif not mesh is None:
assert isinstance(mesh, (MeshReFreq, MeshDLRImFreq, MeshImFreq)), "mesh must be a triqs MeshReFreq or MeshImFreq" assert isinstance(mesh,
(MeshReFreq, MeshDLRImFreq, MeshImFreq)), "mesh must be a triqs MeshReFreq or MeshImFreq"
if isinstance(mesh, MeshImFreq): if isinstance(mesh, MeshImFreq):
mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh)) mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
elif isinstance(mesh, MeshDLRImFreq): elif isinstance(mesh, MeshDLRImFreq):
@ -596,7 +597,7 @@ class SumkDFT(object):
gf_struct = [(spn[isp], block_structure[isp]) gf_struct = [(spn[isp], block_structure[isp])
for isp in range(self.n_spin_blocks[self.SO])] for isp in range(self.n_spin_blocks[self.SO])]
block_ind_list = [block for block, inner in gf_struct] block_ind_list = [block for block, inner in gf_struct]
glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)]) glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner), len(inner)])
for block, inner in gf_struct] for block, inner in gf_struct]
G_latt = BlockGf(name_list=block_ind_list, G_latt = BlockGf(name_list=block_ind_list,
block_list=glist(), make_copies=False) block_list=glist(), make_copies=False)
@ -610,7 +611,7 @@ class SumkDFT(object):
ind = ntoi[spn[ibl]] ind = ntoi[spn[ibl]]
n_orb = self.n_orbitals[ik, ind] n_orb = self.n_orbitals[ik, ind]
if isinstance(mesh, (MeshImFreq, MeshDLRImFreq)): if isinstance(mesh, (MeshImFreq, MeshDLRImFreq)):
gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl)) gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field * (1 - 2 * ibl))
- self.hopping[ik, ind, 0:n_orb, 0:n_orb]) - self.hopping[ik, ind, 0:n_orb, 0:n_orb])
else: else:
gf.data[:, :, :] = (idmat[ibl] * gf.data[:, :, :] = (idmat[ibl] *
@ -655,8 +656,8 @@ class SumkDFT(object):
if (isinstance(self.mesh, (MeshImFreq, MeshDLRImFreq)) and if (isinstance(self.mesh, (MeshImFreq, MeshDLRImFreq)) and
all(isinstance(gf.mesh, (MeshImFreq, MeshDLRImFreq)) and all(isinstance(gf.mesh, (MeshImFreq, MeshDLRImFreq)) and
isinstance(gf, Gf) and isinstance(gf, Gf) and
gf.mesh == self.mesh for bname, gf in Sigma_imp[0])): gf.mesh == self.mesh for bname, gf in Sigma_imp[0])):
# Imaginary frequency Sigma: # Imaginary frequency Sigma:
self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk') self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk')
for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells)]
@ -692,7 +693,7 @@ class SumkDFT(object):
self.max_band_energy - self.chemical_potential): self.max_band_energy - self.chemical_potential):
warn( warn(
'The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f' % ( 'The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f' % (
mesh[0], mesh[-1], self.min_band_energy, self.max_band_energy)) mesh[0], mesh[-1], self.min_band_energy, self.max_band_energy))
def transform_to_sumk_blocks(self, Sigma_imp, Sigma_out=None): def transform_to_sumk_blocks(self, Sigma_imp, Sigma_out=None):
r""" transform Sigma from solver to sumk space r""" transform Sigma from solver to sumk space
@ -1843,7 +1844,8 @@ class SumkDFT(object):
for bname, gf in sigma_minus_dc[icrsh]: for bname, gf in sigma_minus_dc[icrsh]:
# Transform dc_imp to global coordinate system # Transform dc_imp to global coordinate system
if self.use_rotations: if self.use_rotations:
gf -= np.dot(self.rot_mat[icrsh], np.dot(self.dc_imp[icrsh][bname], self.rot_mat[icrsh].conjugate().transpose())) gf -= np.dot(self.rot_mat[icrsh],
np.dot(self.dc_imp[icrsh][bname], self.rot_mat[icrsh].conjugate().transpose()))
else: else:
gf -= self.dc_imp[icrsh][bname] gf -= self.dc_imp[icrsh][bname]
@ -2205,7 +2207,8 @@ class SumkDFT(object):
filename = 'dens_mat.dat' filename = 'dens_mat.dat'
elif dm_type == 'vasp': elif dm_type == 'vasp':
# use new h5 interface to vasp by default, if not wanted specify dm_type='vasp' + filename='GAMMA' # use new h5 interface to vasp by default, if not wanted specify dm_type='vasp' + filename='GAMMA'
filename = 'vaspgamma.h5' # filename = 'vaspgamma.h5'
filename = 'GAMMA'
elif dm_type == 'elk': elif dm_type == 'elk':
filename = 'DMATDMFT.OUT' filename = 'DMATDMFT.OUT'
elif dm_type == 'qe': elif dm_type == 'qe':
@ -2360,7 +2363,7 @@ class SumkDFT(object):
vasp_h5['deltaN'] = deltaN vasp_h5['deltaN'] = deltaN
else: else:
with open(filename, 'w') as f: with open(filename, 'w') as f:
f.write(" %i -1 ! Number of k-points, default number of bands\n" % len(kpts_to_write)) f.write(" -1 -1 ! Number of k-points, default number of bands\n") # % len(kpts_to_write))
for index, ik in enumerate(kpts_to_write): for index, ik in enumerate(kpts_to_write):
ib1 = band_window[0][ik, 0] ib1 = band_window[0][ik, 0]
ib2 = band_window[0][ik, 1] ib2 = band_window[0][ik, 1]
@ -2372,8 +2375,10 @@ class SumkDFT(object):
valim = (deltaN['ud'][ik][inu, imu].imag) / 1.0 valim = (deltaN['ud'][ik][inu, imu].imag) / 1.0
f.write(" %.14f %.14f" % (valre, valim)) f.write(" %.14f %.14f" % (valre, valim))
else: else:
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0 valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0 inu, imu].real) / 2.0
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][
inu, imu].imag) / 2.0
f.write(" %.14f %.14f" % (valre, valim)) f.write(" %.14f %.14f" % (valre, valim))
f.write("\n") f.write("\n")
@ -2396,7 +2401,7 @@ class SumkDFT(object):
mu = self.chemical_potential / self.energy_unit mu = self.chemical_potential / self.energy_unit
# ouput n_k, nspin and max orbitals - a check # ouput n_k, nspin and max orbitals - a check
f.write(" %d %d %d %.14f %.14f ! nkpt, nspin, nstmax, beta, mu\n" % ( f.write(" %d %d %d %.14f %.14f ! nkpt, nspin, nstmax, beta, mu\n" % (
self.n_k, n_spin_blocks, nbmax, beta, mu)) self.n_k, n_spin_blocks, nbmax, beta, mu))
for ik in range(self.n_k): for ik in range(self.n_k):
for ispn in range(n_spin_blocks): for ispn in range(n_spin_blocks):
# Determine the SO density matrix band indices from the spinor band indices # Determine the SO density matrix band indices from the spinor band indices