From 01399a0a9708cfb77c5aaeda0a8b19ba8a9788f7 Mon Sep 17 00:00:00 2001 From: Yuichi Motoyama Date: Thu, 22 Aug 2024 17:01:55 +0900 Subject: [PATCH] update anacont --- src/dcore/anacont/pade.py | 29 +- src/dcore/anacont/spm.py | 50 +- src/dcore/dcore_anacont.py | 4 +- src/dcore/dcore_post.py | 556 +---------------------- src/dcore/program_options.py | 15 +- tests/non-mpi/anacont_spm/anacont_spm.py | 2 +- 6 files changed, 69 insertions(+), 587 deletions(-) diff --git a/src/dcore/anacont/pade.py b/src/dcore/anacont/pade.py index 50939854..88fda347 100644 --- a/src/dcore/anacont/pade.py +++ b/src/dcore/anacont/pade.py @@ -22,17 +22,20 @@ from dcore.program_options import create_parser, parse_parameters -def _set_n_pade(omega_cutoff, beta, n_min, n_max): +def _set_n_matsubara(omega_cutoff, beta, n_min, n_max): """ - Return (int)n_pade: the number of Matsubara frequencies below the cutoff frequency. + Return (int)n_matsubara: the number of Matsubara frequencies below the cutoff frequency. n_pade is bounded between n_min and n_max """ - n_pade = int((beta * omega_cutoff + numpy.pi) / (2.0 * numpy.pi)) - print("n_pade = {} (evaluated from omega_pade)".format(n_pade)) - n_pade = max(n_pade, n_min) - n_pade = min(n_pade, n_max) - print("n_pade = {}".format(n_pade)) - return n_pade + if omega_cutoff > 0.0: + n_matsubara = int((beta * omega_cutoff + numpy.pi) / (2.0 * numpy.pi)) + print("n_matsubara = {} (evaluated from iomega_max)".format(n_matsubara)) + n_matsubara = max(n_matsubara, n_min) + n_matsubara = min(n_matsubara, n_max) + else: + n_matsubara = n_max + print("n_matsubara = {}".format(n_matsubara)) + return n_matsubara def anacont(sigma_iw_npz, beta, mesh_w, params_pade): @@ -40,7 +43,13 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_pade): n_min = params_pade["n_min"] n_max = params_pade["n_max"] eta = params_pade["eta"] - n_pade = _set_n_pade(iw_max, beta, n_min=n_min, n_max=n_max) + + n_matsubara = _set_n_matsubara(iw_max, beta, n_min=n_min, n_max=n_max) + n = sigma_iw_npz['data0'].shape[0] // 2 + if n_matsubara > n: + n_matsubara = n + print("Warning: n_matsubara is larger than the number of calculated Matsubara frequencies, {}".format(n)) + print(" n_matsubara is reset to {}".format(n)) data_w = {} num_data = numpy.sum([key.startswith("data") for key in sigma_iw_npz.keys()]) @@ -50,7 +59,7 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_pade): mesh_iw = MeshImFreq(beta, "Fermion", data.shape[0] // 2) sigma_iw = GfImFreq(data=data, beta=beta, mesh=mesh_iw) sigma_w = GfReFreq(mesh=mesh_w, target_shape=data.shape[1:]) - sigma_w.set_from_pade(sigma_iw, n_points=n_pade, freq_offset=eta) + sigma_w.set_from_pade(sigma_iw, n_points=n_matsubara, freq_offset=eta) data_w[key] = sigma_w.data return data_w diff --git a/src/dcore/anacont/spm.py b/src/dcore/anacont/spm.py index 57f6293f..431ee6ac 100644 --- a/src/dcore/anacont/spm.py +++ b/src/dcore/anacont/spm.py @@ -16,7 +16,9 @@ # along with this program. If not, see . # +import builtins import sys +import re import numpy as np import cvxpy as cp @@ -26,6 +28,13 @@ from dcore._dispatcher import MeshImFreq, GfImFreq, GfReFreq from dcore.program_options import create_parser, parse_parameters + +def __gettype(name): + t = getattr(builtins, name) + if isinstance(t, type): + return t + raise ValueError(name) + def set_default_values(dictionary, default_values_dict): for key, value in default_values_dict.items(): if key not in dictionary: @@ -91,8 +100,8 @@ def get_single_continuation( sum_rule_const, lambd, verbose=True, - max_iters=100, solver="ECOS", + solver_opts={}, ): U, S, Vt, delta_energy, energies_extract = _get_svd_for_continuation( tau_grid, nsv, beta, emin, emax, num_energies @@ -106,8 +115,8 @@ def get_single_continuation( sum_rule_const, lambd, verbose=verbose, - max_iters=max_iters, solver=solver, + solver_opts=solver_opts, ) rho = np.dot(Vt.T, rho_prime) rho_integrated = np.trapz(y=rho, x=energies_extract) @@ -125,8 +134,8 @@ def get_multiple_continuations( sum_rule_const, lambdas, verbose=True, - max_iters=100, solver="ECOS", + solver_opts={}, ): U, S, Vt, delta_energy, energies_extract = _get_svd_for_continuation( tau_grid, nsv, beta, emin, emax, num_energies @@ -144,8 +153,8 @@ def get_multiple_continuations( sum_rule_const, lambd, verbose=verbose, - max_iters=max_iters, solver=solver, + solver_opts=solver_opts, ) rho = np.dot(Vt.T, rho_prime) rho_integrated = np.trapz(y=rho, x=energies_extract) @@ -201,8 +210,8 @@ def _solveProblem( sum_rule_const, lambd, verbose=True, - max_iters=100, solver="ECOS", + solver_opts={}, ): Q = len(S) Smat = np.diag(S) @@ -216,8 +225,9 @@ def _solveProblem( Vt.T @ rho_prime >= 0, cp.sum(delta_energy * V_mod @ rho_prime) == sum_rule_const, ] # uniform real energy grid is assumed here + prob = cp.Problem(objective, constraints) - _ = prob.solve(verbose=verbose, solver=solver, max_iters=max_iters) + _ = prob.solve(verbose=verbose, solver=solver, **solver_opts) gf_tau_fit = np.dot(U, np.dot(Smat, rho_prime.value)) chi2 = 0.5 * np.linalg.norm(gf_tau - gf_tau_fit, ord=2) ** 2 return rho_prime.value, gf_tau_fit, chi2 @@ -282,13 +292,28 @@ def check_parameters(params): pass def parameters_from_ini(inifile): - parser = create_parser(["post.anacont.spm"]) + parser = create_parser(["post.anacont.spm", "post.anacont.spm.solver"]) parser.read(inifile) params = parser.as_dict() parse_parameters(params) - return params["post.anacont.spm"] - -def anacont(sigma_iw_npz, beta, mesh_w, params_spm): + solver_dict = {} + if "post.anacont.spm.solver" in params: + r = re.compile('^(.*)\{(.*)\}$') + for k,v in params["post.anacont.spm.solver"].items(): + try: + m = r.search(k) + if m is None: + continue + + param_name, param_type_str = m.group(1), m.group(2) + param_type = __gettype(param_type_str) + except RuntimeError: + raise RuntimeError("Unknown type or unrecognized format : " + k) + solver_dict[param_name] = param_type(v) + params["post.anacont.spm.solver"] = solver_dict + return params + +def anacont(sigma_iw_npz, beta, mesh_w, params_spm, params_spm_solver): n_tau = params_spm["n_tau"] n_tail = params_spm["n_tail"] @@ -296,8 +321,7 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_spm): n_sv = params_spm["n_sv"] L1_coeff = params_spm["lambda"] - max_iters = params_spm["max_iters_opt"] - solver_opt = params_spm["solver_opt"] + solver = params_spm["solver"] verbose_opt = params_spm["verbose_opt"] data_w = {} @@ -327,7 +351,7 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_spm): ) density, gf_tau_fit, energies_extract, sum_rule_const, chi2 = get_single_continuation( - tau_grid, gf_tau, n_sv, beta, mesh_w.values()[0], mesh_w.values()[-1], mesh_w.size, sum_rule_const=const_imag_tail, lambd=L1_coeff, verbose=verbose_opt, max_iters=max_iters, solver=solver_opt + tau_grid, gf_tau, n_sv, beta, mesh_w.values()[0], mesh_w.values()[-1], mesh_w.size, sum_rule_const=const_imag_tail, lambd=L1_coeff, verbose=verbose_opt, solver=solver, solver_opts=params_spm_solver ) energies, gf_real, gf_imag = get_kramers_kronig_realpart( energies_extract, dos_to_gf_imag(density) diff --git a/src/dcore/dcore_anacont.py b/src/dcore/dcore_anacont.py index 18871a9d..4104ab7b 100644 --- a/src/dcore/dcore_anacont.py +++ b/src/dcore/dcore_anacont.py @@ -57,7 +57,9 @@ def dcore_anacont(inifile): data_w = pade.anacont(sigma_iw_npz, beta, mesh_w, params_ac) elif solver == "spm": params_ac = spm.parameters_from_ini(inifile) - data_w = spm.anacont(sigma_iw_npz, beta, mesh_w, params_ac) + if "post.anacont.spm.solver" not in params_ac: + params_ac["post.anacont.spm.solver"] = {} + data_w = spm.anacont(sigma_iw_npz, beta, mesh_w, params_ac["post.anacont.spm"], params_ac["post.anacont.spm.solver"]) else: assert False, "Unknown solver: " + solver diff --git a/src/dcore/dcore_post.py b/src/dcore/dcore_post.py index ab517e59..ae899e73 100644 --- a/src/dcore/dcore_post.py +++ b/src/dcore/dcore_post.py @@ -16,566 +16,16 @@ # along with this program. If not, see . # -import os import sys -import numpy -import copy -from itertools import product - -from dcore._dispatcher import * -from dcore.dmft_core import DMFTCoreSolver -from dcore.program_options import create_parser, parse_parameters, parse_bvec, _set_nk, print_parameters, delete_parameters -from dcore.tools import save_Sigma_w_sh_txt -from dcore import impurity_solvers -from dcore.lattice_models import create_lattice_model -from .sumkdft_workers.launcher import run_sumkdft - - -def _read_sigma_w(npz_file, nsh, mesh, block_names): - npz = np.load(npz_file) - Sigma_iw = [] - idx = 0 - for _ in range(nsh): - block_list = [] - for _ in range(len(block_names)): - block_list.append(GfReFreq(data=npz[f'data{idx}'], mesh=mesh)) - idx += 1 - G = BlockGf( - name_list = block_names, - block_list = block_list, make_copies = False) - Sigma_iw.append(G) - return Sigma_iw - -class DMFTPostSolver(DMFTCoreSolver): - def __init__(self, seedname, params, output_file='', output_group='dmft_out'): - - super(DMFTPostSolver, self).__init__(seedname, params, output_file, output_group, read_only=True, restart=True) - - - def calc_dos(self, Sigma_w_sh, mesh, broadening): - """ - - Compute dos in real frequency. - - :param Sigma_w_sh: list - List of real-frequency self-energy - - :param broadening: float - Broadening factor - - :return: tuple - Results are 'dos', 'dosproj', 'dosproj_orb'. - - """ - - params = self._make_sumkdft_params() - params['mu'] = self._chemical_potential - params['Sigma_w_sh'] = Sigma_w_sh - params['mesh'] = mesh - params['broadening'] = broadening - r = run_sumkdft( - 'SumkDFTWorkerDOS', - os.path.abspath(self._seedname+'.h5'), './work/sumkdft_dos', self._mpirun_command, params) - return r['dos'], r['dosproj'], r['dosproj_orb'] - - def calc_spaghettis(self, Sigma_w_sh, mesh, broadening, kmesh_type): - """ - - Compute A(k, omega) - - """ - - params = self._make_sumkdft_params() - params['calc_mode'] = 'spaghettis' - params['mu'] = self._chemical_potential - params['Sigma_w_sh'] = Sigma_w_sh - params['mesh'] = mesh - params['broadening'] = broadening - if kmesh_type == 'line': - params['bands_data'] = 'dft_bands_input' - elif kmesh_type == 'mesh': - params['bands_data'] = 'dft_bands_mesh_input' - else: - raise RuntimeError('Invalid kmesh_type: {}'.format(kmesh_type)) - #r = sumkdft.run(os.path.abspath(self._seedname+'.h5'), './work/sumkdft_spaghettis', self._mpirun_command, params) - r = run_sumkdft( - 'SumkDFTWorkerSpaghettis', - os.path.abspath(self._seedname+'.h5'), './work/sumkdft_spaghettis', self._mpirun_command, params) - return r['akw'] - - def calc_momentum_distribution(self): - """ - - Compute momentum distributions and eigen values of H(k) - Data are taken from bands_data. - - """ - - params = self._make_sumkdft_params() - params['calc_mode'] = 'momentum_distribution' - params['mu'] = self._chemical_potential - r = run_sumkdft( - 'SumkDFTWorkerMomDist', - os.path.abspath(self._seedname+'.h5'), './work/sumkdft_mom_dist', self._mpirun_command, params) - return r['den'] - - def calc_Sigma_w(self, mesh): - """ - Compute Sigma_w for computing band structure - For an imaginary-time solver, a list of Nones will be returned. - - :param mesh: (float, float, int) - real-frequency mesh (min, max, num_points) - - :return: list of Sigma_w - - """ - - solver_name = self._params['impurity_solver']['name'] - Solver = impurity_solvers.solver_classes[solver_name] - if Solver.is_gf_realomega_available(): - Gloc_iw_sh, _, _ = self.calc_Gloc() - _, _, sigma_w = self.solve_impurity_models(Gloc_iw_sh, -1, mesh) - return sigma_w - else: - return [None] * self.n_inequiv_shells - - -def _set_n_pade(omega_cutoff, beta, n_min, n_max): - """ - Return (int)n_pade: the number of Matsubara frequencies below the cutoff frequency. - n_pade is bounded between n_min and n_max - """ - n_pade = int((beta * omega_cutoff + numpy.pi) / (2.0 * numpy.pi)) - print("n_pade = {} (evaluated from omega_pade)".format(n_pade)) - n_pade = max(n_pade, n_min) - n_pade = min(n_pade, n_max) - print("n_pade = {}".format(n_pade)) - return n_pade - - -class DMFTCoreTools: - def __init__(self, seedname, params, n_k, xk, nkdiv_mesh, kvec_mesh, prefix): - """ - Class of posting tool for DCore. - - Parameters - ---------- - :param seedname: string - name for hdf5 file - :param params: dictionary - Input parameters - :param n_k: integer - Number of k points - :param xk: integer array - x-position for plotting band - :param nkdiv_mesh: (int, int, int) - Number of k points along each axis for computing A(k, omega) - :param kvec_mesh: float array of dimension (*, 3) - k points in fractional coordinates for computing A(k, omega) on a mesh - """ - - self._params = copy.deepcopy(params) - # Construct a SumKDFT object - self._n_pade = _set_n_pade(omega_cutoff=params['tool']['omega_pade'], - beta=params['system']['beta'], - n_min=params['tool']['n_pade_min'], - n_max=params['tool']['n_pade_max']) - self._omega_min = float(params['tool']['omega_min']) - self._omega_max = float(params['tool']['omega_max']) - self._Nomega = int(params['tool']['Nomega']) - self._broadening = float(params['tool']['broadening']) - self._eta = float(params['tool']['eta']) - self._seedname = seedname - self._n_k = n_k - self._xk = xk - self._kvec_mesh = kvec_mesh - self._nkdiv_mesh = nkdiv_mesh - self._prefix = prefix - - self._solver = DMFTPostSolver(seedname, self._params, output_file=seedname+'.out.h5') - print("iteration :", self._solver.iteration_number) - - def print_dos(self, dos, dosproj_orb, filename): - """ - - Print DOS to file - - """ - nsh = self._solver.n_inequiv_shells - om_mesh = numpy.linspace(self._omega_min, self._omega_max, self._Nomega) - spin_block_names = self._solver.spin_block_names - inequiv_to_corr = self._solver.inequiv_to_corr - corr_shell_info = [self._solver.corr_shell_info(ish) for ish in range(self._solver._n_corr_shells)] - - with open(filename, 'w') as f: - # - # Description of columns - # - print("# [1] Energy", file=f) - ii = 1 - for isp in spin_block_names: - ii += 1 - print("# [%d] Total DOS of spin %s" % (ii, isp), file=f) - for ish in range(nsh): - block_dim = corr_shell_info[inequiv_to_corr[ish]]['block_dim'] - for isp in spin_block_names: - for iorb in range(block_dim): - ii += 1 - print("# [%d] PDOS of shell%d,spin %s,band%d" % (ii, ish, isp, iorb), file=f) - # - for iom in range(self._Nomega): - print("%f" % om_mesh[iom], file=f, end="") - for isp in spin_block_names: - print(" %f" % dos[isp][iom], file=f, end="") - for ish in range(nsh): - block_dim = corr_shell_info[inequiv_to_corr[ish]]['block_dim'] - for isp in spin_block_names: - for iorb in range(block_dim): - print(" %f" % dosproj_orb[ish][isp][iom, iorb, iorb].real, end="", file=f) - print("", file=f) - print("\n Output {0}".format(filename)) - - def print_band(self, akw, filename): - """ - Print A(k,w) into a file - - Parameters - ---------- - akw - filename - - """ - - om_mesh = numpy.linspace(self._omega_min, self._omega_max, self._Nomega) - with open(filename, 'w') as f: - offset = 0.0 - for isp in self._solver.spin_block_names: - for ik, xk in enumerate(self._xk): - for iom, om in enumerate(om_mesh): - print("%f %f %f" % (xk+offset, om, akw[isp][ik, iom]), file=f) - print("", file=f) - offset = self._xk[-1] * 1.1 - print("", file=f) - print("\n Output {0}".format(filename)) - - def post(self): - """ - Calculate DOS (Density Of State) and energy dispersions. - For Hubbard-I solver, self-energy is calculated in this function. - For cthyb (both TRIQS and ALPS), self-energy is read from hdf5 file. - """ - - print("\n############# Compute Green's Function in the Real Frequency ################\n") - - # - # Real-frequency self-energy - # - mesh = [self._omega_min, self._omega_max, self._Nomega] - sigma_w_sh = self._solver.calc_Sigma_w(mesh) - Sigma_iw_sh = self._solver.Sigma_iw_sh(self._solver.iteration_number) - - filename = self._seedname + '_sigma_w.npz' - if os.path.exists(filename): - print(f"Reading sigma_w from {filename}...") - sigma_w_sh = _read_sigma_w(filename, len(sigma_w_sh), - MeshReFreq(*mesh), self._solver.spin_block_names) - else: - # Backward compatibility - print(f"Not found {filename}. Falling back to pade approximant...") - for ish in range(self._solver.n_inequiv_shells): - if not sigma_w_sh[ish] is None: - continue - # set BlockGf sigma_w - Sigma_iw = Sigma_iw_sh[ish] - block_names = self._solver.spin_block_names - def glist(): - return [GfReFreq(indices=sigma.indices, window=(self._omega_min, self._omega_max), - n_points=self._Nomega, name="sig_pade") for block, sigma in Sigma_iw] - sigma_w_sh[ish] = BlockGf(name_list=block_names, block_list=glist(), make_copies=False) - # Analytic continuation - for bname, sig in Sigma_iw: - sigma_w_sh[ish][bname].set_from_pade(sig, n_points=self._n_pade, freq_offset=self._eta) - - print("\n############# Print Self energy in the Real Frequency ################\n") - filename = self._prefix + self._seedname + '_sigmaw.dat' - print("\n Writing real-freqnecy self-energy into ", filename) - save_Sigma_w_sh_txt(filename, sigma_w_sh, self._solver.spin_block_names) - - # - # (Partial) DOS - # - print("\n############# Compute (partial) DOS ################\n") - dos, dosproj, dosproj_orb = self._solver.calc_dos(sigma_w_sh, mesh, self._broadening) - self.print_dos(dos, dosproj_orb, self._prefix + self._seedname+'_dos.dat') - - - # - # Band structure - # - if not self._xk is None: - print("\n############# Compute Band Structure ################\n") - akw = self._solver.calc_spaghettis(sigma_w_sh, mesh, self._broadening, 'line') - self.print_band(akw, self._prefix + self._seedname + '_akw.dat') - - # - # A(k, omega) on a mesh - # - nk_mesh = numpy.prod(self._nkdiv_mesh) - if nk_mesh > 0: - print("\n############# Compute A(k, omega) on a mesh ################\n") - akw = self._solver.calc_spaghettis(sigma_w_sh, mesh, self._broadening, 'mesh') - #print("debug", mesh) - #print("debugB", len(mesh)) - #print("debugC", self._Nomega) - om_mesh = numpy.linspace(mesh[0], mesh[1], mesh[2]) - bvec = parse_bvec(self._params["model"]["bvec"]) - for bname in self._solver.spin_block_names: - filename = self._prefix + self._seedname + '_akw_mesh_{}.dat'.format(bname) - print("\n Output {0}".format(filename)) - with open(filename, 'w') as f: - print('# {} {} {} {}'.format(self._nkdiv_mesh[0], self._nkdiv_mesh[1], self._nkdiv_mesh[2], self._Nomega), file=f) - for i in range(3): - print('# {} {} {}'.format(*bvec[i, :]), file=f) - for ik in range(nk_mesh): - for iom in range(self._Nomega): - print("%f %f %f %f %f" % (self._kvec_mesh[ik,0], - self._kvec_mesh[ik,1], - self._kvec_mesh[ik,2], - om_mesh[iom], akw[bname][ik, iom]), file=f) - - def momentum_distribution(self): - """ - Calculate Momentum distribution - """ - if self._xk is None: - return - - print("\n############# Momentum Distribution ################\n") - - den = self._solver.calc_momentum_distribution() - - spn = self._solver.spin_block_names - - n_k, n_orbitals = den.shape[0], den.shape[2] - - SO = 1 if self._solver.use_spin_orbit else 0 - - # - # Output momentum distribution to file - # - filename = self._prefix + self._seedname + "_momdist.dat" - print("\n Output Momentum distribution : ", filename) - with open(filename, 'w') as fo: - print("# Momentum distribution", file=fo) - # - # Column information - # - print("# [Column] Data", file=fo) - print("# [1] Distance along k-path", file=fo) - icol = 1 - for isp in spn: - for iorb in range(n_orbitals): - for jorb in range(n_orbitals): - icol += 1 - print("# [%d] Re(MomDist_{spin=%s, %d, %d})" % (icol, isp, iorb, jorb), file=fo) - icol += 1 - print("# [%d] Im(MomDist_{spin=%s, %d, %d})" % (icol, isp, iorb, jorb), file=fo) - # - # Write data - # - for ik in range(n_k): - print("%f " % self._xk[ik], end="", file=fo) - for isp in range(2-SO): - for iorb in range(n_orbitals): - for jorb in range(n_orbitals): - print("%f %f " % (den[ik, isp, iorb, jorb].real, - den[ik, isp, iorb, jorb].imag), end="", file=fo) - print("", file=fo) - - -def __print_parameter(p, param_name): - """ - Print parameters. - - Parameters - ---------- - p : dictionary - Dictionary for parameters - param_name : string - key for p - """ - print(param_name + " = " + str(p[param_name])) - - -def gen_script_gnuplot(xnode, seedname, prefix, spin_orbit): - file_akw_gp = prefix + seedname + '_akw.gp' - - def print_klabel(label, x, f, with_comma=True): - print(' "{}" {}'.format(label, x), end="", file=f) - if with_comma: - print(',', end="", file=f) - print(' \\', file=f) - - k_end = len(xnode) - 1 - - with open(file_akw_gp, 'w') as f: - print("set size 0.95, 1.0", file=f) - print("set xtics (\\", file=f) - if spin_orbit: - for i, node in enumerate(xnode): - print_klabel(node.label, node.x, f, i != k_end) - else: - for node in xnode: - print_klabel(node.label, node.x, f) - offset = xnode[-1].x * 1.1 - for i, node in enumerate(xnode): - print_klabel(node.label, node.x + offset, f, i != k_end) - print(" )", file=f) - print("set pm3d map", file=f) - print("#set pm3d interpolate 5, 5", file=f) - print("unset key", file=f) - print("set ylabel \"Energy\"", file=f) - print("set cblabel \"A(k,w)\"", file=f) - print("splot \"{0}_akw.dat\" u 1:2:(abs($3))".format(seedname), file=f) - print("pause -1", file=f) - - print(" Usage:") - print("\n $ gnuplot {0}".format(os.path.basename(file_akw_gp))) def dcore_post(filename, np=1, prefix=None): """ - Main routine for the post-processing tool - - Parameters - ---------- - filename : string - Input-file name + Removed function. Use dcore_anacont and dcore_spectrum instead. """ - print("\n############ Reading Input File #################\n") - print(" Input File Name : ", filename) - # - # Construct a parser with default values - # - pars = create_parser(['model', 'system', 'impurity_solver', 'tool', 'post', 'mpi']) - # - # Parse keywords and store - # - pars.read(filename) - p = pars.as_dict() - parse_parameters(p) - seedname = p["model"]["seedname"] - p["mpi"]["num_processes"] = np - mpirun_command = p['mpi']['command'].replace('#', str(p['mpi']['num_processes'])) - mpirun_command_np1 = p['mpi']['command'].replace('#', '1') - - # for backward compatibility - if prefix is not None: - p["tool"]["post_dir"] = prefix - - # - # Delete unnecessary parameters - # - delete_parameters(p, block='model', delete=['interaction', 'density_density', 'kanamori', 'slater_f', 'slater_uj', 'slater_basis', 'interaction_file', 'local_potential_matrix', 'local_potential_factor']) - - # Summary of input parameters - print_parameters(p) - - # make directory - post_dir = p["tool"]["post_dir"] - if post_dir: - os.makedirs(post_dir, exist_ok=True) - prefix = post_dir + "/" - else: - prefix = "./" - - # - # Generate k-path and compute H(k) on this path - # - print("\n################ Generating k-path ##################\n") - lattice_model = create_lattice_model(p) - xk, xnode = lattice_model.generate_Hk_path(p) - - if xk is None: - n_k = 0 - print(' A(k,w) calc will be skipped') - else: - n_k = len(xk) - print(" Total number of k =", n_k) - print(" k-point x") - for node in xnode: - print(" %6s %f" %(node.label, node.x)) - - # - # Generate k mesh and compute H(k) on the mesh - # - nk_div = _set_nk(p["tool"]["nk_mesh"], p["tool"]["nk0_mesh"], p["tool"]["nk1_mesh"], p["tool"]["nk2_mesh"]) - kvec_mesh = None - if all(div != 0 for div in nk_div): - print("\n################ Constructing H(k) for compute A(k, omega) on a mesh ##################") - k = [numpy.linspace(0, 2*numpy.pi, div+1)[:-1] for div in nk_div] - kvec_mesh = numpy.array([kxyz for kxyz in product(k[0], k[1], k[2])]) - lattice_model.write_dft_band_input_data(p, kvec_mesh, bands_data='dft_bands_mesh_input') - - # - # Compute DOS and A(k,w) - # - print("\n############# Run DMFTCoreTools ########################\n") - dct = DMFTCoreTools(seedname, p, n_k, xk, nk_div, kvec_mesh, prefix) - dct.post() - dct.momentum_distribution() - - # - # Output gnuplot script - # - if xnode is not None: - print("\n############# Generate GnuPlot Script ########################\n") - gen_script_gnuplot(xnode, seedname, prefix, p["model"]["spin_orbit"]) - - print("\n################# Done #####################\n") + sys.exit("dcore_post is removed. Use dcore_anacont and dcore_spectrum instead.") def run(): - import argparse - from dcore.option_tables import generate_all_description - from dcore.version import version, print_header - - print_header() - - parser = argparse.ArgumentParser( - prog='dcore_post.py', - description='Post-processing script in DCore', - # usage='$ dcore_post input --np 4', - add_help=True, - formatter_class=argparse.RawTextHelpFormatter, - epilog=generate_all_description() - ) - parser.add_argument('path_input_files', - action='store', - default=None, - type=str, - nargs='*', - help="Input filename(s)", - ) - parser.add_argument('--np', help='Number of MPI processes', required=True) - parser.add_argument('--version', action='version', version='DCore {}'.format(version)) - parser.add_argument('--prefix', # Deprecated - action='store', - default=None, - type=str, - # help='prefix for output files (default: post/)' - help='[Deprecated] Use post_dir parameter in [tool] block', - ) - - args = parser.parse_args() - - # for backward compatibility - if args.prefix is not None: - print("DeprecationWarning: --prefix option is deprecated. Use post_dir parameter in [tool] block", file=sys.stderr) - - for path_input_file in args.path_input_files: - if os.path.isfile(path_input_file) is False: - sys.exit(f"Input file '{path_input_file}' does not exist.") - dcore_post(args.path_input_files, int(args.np), args.prefix) + sys.exit("dcore_post is removed from DCore v4. Use dcore_anacont and dcore_spectrum instead.") diff --git a/src/dcore/program_options.py b/src/dcore/program_options.py index 08bc3879..36f5880a 100644 --- a/src/dcore/program_options.py +++ b/src/dcore/program_options.py @@ -45,12 +45,13 @@ def create_parser(target_sections=None): Create a parser for all program options of DCore """ if target_sections is None: - parser = TypedParser(['mpi', 'model', 'pre', 'system', 'impurity_solver', 'control', 'post', 'post.anacont', 'post.anacont.pade', 'post.anacont.spm', 'post.spectrum', 'post.check', 'bse', 'vertex', 'sparse_bse']) + parser = TypedParser(['mpi', 'model', 'pre', 'system', 'impurity_solver', 'control', 'post', 'post.anacont', 'post.anacont.pade', 'post.anacont.spm', 'post.anacont.spm.solver' 'post.spectrum', 'post.check', 'bse', 'vertex', 'sparse_bse']) else: parser = TypedParser(target_sections) # tool is removed but read for warning parser.allow_undefined_options("tool") + parser.allow_undefined_options("post.anacont.spm.solver") # [mpi] parser.add_option("mpi", "command", str, default_mpi_command, "Command for executing a MPI job. # will be relaced by the number of processes.") @@ -140,9 +141,9 @@ def create_parser(target_sections=None): parser.add_option("post.anacont", "save_result", bool, False, "plot result of analytic continuation") # [post.anacont.pade] - parser.add_option( "post.anacont.pade", "iomega_max", float, -1.0, "Cut-off frequency of the Matsubara frequency",) - parser.add_option( "post.anacont.pade", "n_min", int, 0, "lower bound of the order of Pade approximation",) - parser.add_option( "post.anacont.pade", "n_max", int, 100000000, "upper bound of the order of Pade approximation",) + parser.add_option( "post.anacont.pade", "iomega_max", float, 0.0, "Cut-off frequency of the Matsubara frequency",) + parser.add_option( "post.anacont.pade", "n_min", int, 0, "lower bound of the number of used Matsubara frequency",) + parser.add_option( "post.anacont.pade", "n_max", int, 100000000, "upper bound of the number of used Matsubara frequency",) parser.add_option( "post.anacont.pade", "eta", float, 0.01, "Imaginary Frequency shift to avoid divergence",) # [post.anacont.spm] @@ -151,9 +152,7 @@ def create_parser(target_sections=None): parser.add_option("post.anacont.spm", "n_tail", int, 10, "number of matsubara points for tail-fitting") parser.add_option("post.anacont.spm", "n_sv", int, 50, "number of singular values to be used") parser.add_option("post.anacont.spm", "lambda", float, 1e-5, "coefficient of L1 regularization") - parser.add_option("post.anacont.spm", "max_iters_opt", int, 100, "maximum number of iterations") - parser.add_option("post.anacont.spm", "solver_opt", str, "ECOS", "solver to be used") - + parser.add_option("post.anacont.spm", "solver", str, "ECOS", "solver to be used") parser.add_option("post.anacont.spm", "verbose_opt", bool, False, "show optimization progress") parser.add_option("post.anacont.spm", "show_fit", bool, False, "plot result of tail-fitting") @@ -322,8 +321,6 @@ def parse_parameters(params): sys.exit(f"ERROR: n_sv={params['post.anacont.spm']['n_sv']} must be a positive integer.") if params['post.anacont.spm']['lambda'] < 0: sys.exit(f"ERROR: lambda={params['post.anacont.spm']['lambda']} must be a non-negative float.") - if params['post.anacont.spm']['max_iters_opt'] <= 0: - sys.exit(f"ERROR: max_iters_opt={params['post.anacont.spm']['max_iters_opt']} must be a positive integer.") if 'bse' in params: two_options_incompatible(params, ('bse', 'skip_Xloc'), ('bse', 'calc_only_chiloc')) diff --git a/tests/non-mpi/anacont_spm/anacont_spm.py b/tests/non-mpi/anacont_spm/anacont_spm.py index 93c4d530..92c7d60a 100644 --- a/tests/non-mpi/anacont_spm/anacont_spm.py +++ b/tests/non-mpi/anacont_spm/anacont_spm.py @@ -193,7 +193,7 @@ def test_solveProblem(): U, S, Vt, delta_energy, energies_extract = _get_svd_for_continuation(tau_grid, nsv, beta, emin, emax, num_energies) - rho_prime, gf_tau_fit, chi2 = _solveProblem(delta_energy, U, S, Vt, gf_tau, b_test, lambd, verbose=False, max_iters=100, solver='ECOS') + rho_prime, gf_tau_fit, chi2 = _solveProblem(delta_energy, U, S, Vt, gf_tau, b_test, lambd, verbose=False, solver='ECOS') rho = np.dot(Vt.T, rho_prime) d_bhatt = np.trapz(y=np.sqrt(np.multiply(rho, dos)), x=energies_extract)