diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 00000000..26d33521 --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,3 @@ +# Default ignored files +/shelf/ +/workspace.xml diff --git a/.idea/dft_tools.src.iml b/.idea/dft_tools.src.iml new file mode 100644 index 00000000..8b8c3954 --- /dev/null +++ b/.idea/dft_tools.src.iml @@ -0,0 +1,12 @@ + + + + + + + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/Project_Default.xml b/.idea/inspectionProfiles/Project_Default.xml new file mode 100644 index 00000000..969ecafa --- /dev/null +++ b/.idea/inspectionProfiles/Project_Default.xml @@ -0,0 +1,18 @@ + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 00000000..105ce2da --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 00000000..a8325ac5 --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 00000000..94a25f7f --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/python/triqs_dft_tools/converters/plovasp/sc_dmft.py b/python/triqs_dft_tools/converters/plovasp/sc_dmft.py index 654185d7..f62eac09 100644 --- a/python/triqs_dft_tools/converters/plovasp/sc_dmft.py +++ b/python/triqs_dft_tools/converters/plovasp/sc_dmft.py @@ -72,22 +72,26 @@ def is_vasp_running(vasp_pid): pid_exists = mpi.bcast(pid_exists) return pid_exists + 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: - nextline = f.readline() - while nextline.strip(): - line = nextline - nextline = f.readline() -# print "OSZICAR: ", line[:-1] + h5_energy = False + if os.path.isfile('vaspout.h5'): + with HDFArchive('vaspout.h5', 'r') as h5: + if 'oszicar' in h5['intermediate/ion_dynamics']: + dft_energy = h5['intermediate/ion_dynamics/oszicar'][-1,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]) - except ValueError: - print("Cannot read energy from OSZICAR, setting it to zero") - dft_energy = 0.0 return dft_energy @@ -98,9 +102,8 @@ class bcolors: YELLOW = '\033[93m' RED = '\033[91m' ENDC = '\033[0m' -# + # Main self-consistent cycle -# def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version): """ """ @@ -117,15 +120,14 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version): mpi.barrier() while is_vasp_lock_present(): 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): mpi.report(" VASP stopped") vasp_running = False 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) err = 0 exc = None @@ -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 iter_dft = 0 - if vasp_version == 'standard': + if vasp_version == 'standard' or vasp_version == 'ncl': copyfile(src='GAMMA',dst='GAMMA_recent') while iter_dft < n_iter_dft: # 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 break iter_dft += 1 - if vasp_version == 'standard': + if vasp_version == 'standard' or vasp_version == 'ncl': copyfile(src='GAMMA_recent',dst='GAMMA') iter += 1 if iter == n_iter: @@ -253,8 +255,8 @@ def main(): except KeyError: vasp_version = 'standard' - if vasp_version != 'standard' and vasp_version != 'no_gamma_write': - raise Exception('vasp_version has to be standard or 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') # if len(sys.argv) > 1: # vasp_path = sys.argv[1] diff --git a/python/triqs_dft_tools/converters/plovasp/vaspio.py b/python/triqs_dft_tools/converters/plovasp/vaspio.py index 4afb0692..41a8d813 100644 --- a/python/triqs_dft_tools/converters/plovasp/vaspio.py +++ b/python/triqs_dft_tools/converters/plovasp/vaspio.py @@ -82,7 +82,7 @@ class VaspData: self.plocar = h5Plocar(h5path) self.poscar = h5Poscar(h5path) self.kpoints = h5Kpoints(h5path) - self.eigenval = h5Eigenval(h5path) + self.eigenval = h5Eigenval(h5path, self.kpoints.ksymmap) self.doscar = h5Doscar(h5path) else: self.plocar = Plocar() @@ -716,28 +716,34 @@ class h5Kpoints: # h5path = './vasptriqs.h5' with HDFArchive(h5path, 'a') as archive: kpoints = archive['results/electron_eigenvalues'] - self.nktot = kpoints['kpoints'] - self.kpts = kpoints['kpoint_coords'] - self.kwghts = kpoints['kpoints_symmetry_weight'] + self.nkred = kpoints['kpoints'] + self.kpts = kpoints['kpoint_coords_full'] + self.nktot = len(self.kpts) + self.kwghts = kpoints['kpoints_symmetry_weight_full'] + self.ksymmap = kpoints['kpoints_symmetry_mapping'] + self.ksymmap -= 1 try: self.ntet = kpoints['num_tetrahedra'] self.vtet = kpoints['volume_weight_tetrahedra'] self.itet = kpoints['coordinate_id_tetrahedra'] except KeyError: - print(" No tetrahedron data found in vasptriqs.h5. Skipping...") + print(" No tetrahedron data found in vaspout.h5. Skipping...") self.ntet = 0 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 tetrahedra:", self.ntet)) class h5Eigenval: - def __init__(self, h5path): + def __init__(self, h5path, symmap): with HDFArchive(h5path, 'a') as archive: self.eigs = archive['results/electron_eigenvalues']['eigenvalues'] + self.eigs = self.eigs[:, symmap, :] self.ferw = archive['results/electron_eigenvalues']['fermiweights'] + self.ferw = self.ferw[:, symmap, :] # TODO Change the format in VASP to have [kpoints, bands, spin] self.eigs = np.transpose(self.eigs, (1, 2, 0)) self.ferw = np.transpose(self.ferw, (1, 2, 0)) diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py index 3bf3b5c2..1cf49deb 100644 --- a/python/triqs_dft_tools/sumk_dft.py +++ b/python/triqs_dft_tools/sumk_dft.py @@ -578,7 +578,8 @@ class SumkDFT(object): mesh_values = self.mesh_values 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): mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh)) elif isinstance(mesh, MeshDLRImFreq): @@ -596,7 +597,7 @@ class SumkDFT(object): gf_struct = [(spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO])] 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] G_latt = BlockGf(name_list=block_ind_list, block_list=glist(), make_copies=False) @@ -610,7 +611,7 @@ class SumkDFT(object): ind = ntoi[spn[ibl]] n_orb = self.n_orbitals[ik, ind] 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]) else: gf.data[:, :, :] = (idmat[ibl] * @@ -655,8 +656,8 @@ class SumkDFT(object): if (isinstance(self.mesh, (MeshImFreq, MeshDLRImFreq)) and all(isinstance(gf.mesh, (MeshImFreq, MeshDLRImFreq)) and - isinstance(gf, Gf) and - gf.mesh == self.mesh for bname, gf in Sigma_imp[0])): + isinstance(gf, Gf) and + gf.mesh == self.mesh for bname, gf in Sigma_imp[0])): # Imaginary frequency Sigma: 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)] @@ -692,7 +693,7 @@ class SumkDFT(object): self.max_band_energy - self.chemical_potential): 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' % ( - 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): r""" transform Sigma from solver to sumk space @@ -1843,7 +1844,8 @@ class SumkDFT(object): for bname, gf in sigma_minus_dc[icrsh]: # Transform dc_imp to global coordinate system 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: gf -= self.dc_imp[icrsh][bname] @@ -2205,7 +2207,8 @@ class SumkDFT(object): filename = 'dens_mat.dat' elif dm_type == 'vasp': # 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': filename = 'DMATDMFT.OUT' elif dm_type == 'qe': @@ -2360,7 +2363,7 @@ class SumkDFT(object): vasp_h5['deltaN'] = deltaN else: 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): ib1 = band_window[0][ik, 0] ib2 = band_window[0][ik, 1] @@ -2372,8 +2375,10 @@ class SumkDFT(object): valim = (deltaN['ud'][ik][inu, imu].imag) / 1.0 f.write(" %.14f %.14f" % (valre, valim)) else: - valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0 - valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0 + valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][ + 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("\n") @@ -2396,7 +2401,7 @@ class SumkDFT(object): mu = self.chemical_potential / self.energy_unit # ouput n_k, nspin and max orbitals - a check 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 ispn in range(n_spin_blocks): # Determine the SO density matrix band indices from the spinor band indices