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)