Merge pull request #82 from TRIQS/solver_interface_refactoring

Simplified solver interface for adding new solvers

python/solid_dmft/dmft_cycle.py

@@ -48,7 +48,7 @@
 from solid_dmft.version import version as solid_dmft_version
 from solid_dmft.dmft_tools.observables import (calc_dft_kin_en, add_dmft_observables, calc_bandcorr_man, write_obs,
                                          add_dft_values_as_zeroth_iteration, write_header_to_file, prep_observables)
-from solid_dmft.dmft_tools.solver import SolverStructure
+from solid_dmft.dmft_tools.solver import create_solver
 from solid_dmft.dmft_tools import convergence
 from solid_dmft.dmft_tools import formatter
 from solid_dmft.dmft_tools import interaction_hamiltonian
@@ -514,9 +514,17 @@ def dmft_cycle(general_params, solver_params, advanced_params, dft_params,
     solvers = [None] * sum_k.n_inequiv_shells
     for icrsh in range(sum_k.n_inequiv_shells):
         # Construct the Solver instances
-        solvers[icrsh] = SolverStructure(general_params, solver_params[map_imp_solver[icrsh]],
-                                         gw_params, advanced_params, sum_k, icrsh, h_int[icrsh],
-                                         iteration_offset, deg_orbs_ftps)
+        mpi.report("Creating solver with new interface")
+        solvers[icrsh] = create_solver(general_params = general_params,
+                                       solver_params = solver_params[map_imp_solver[icrsh]],
+                                       gw_params = gw_params,
+                                       advanced_params = advanced_params,
+                                       sum_k= sum_k,
+                                       h_int = h_int[icrsh],
+                                       icrsh= icrsh,
+                                       iteration_offset=iteration_offset,
+                                       deg_orbs_ftps=deg_orbs_ftps)
+
 
     # store solver hash to archive
     if mpi.is_master_node():

python/solid_dmft/dmft_tools/common.py

@@ -0,0 +1,26 @@
+# miscellanous functions for the DMFT calculation used in other modules
+
+def get_n_orbitals(sum_k):
+    """
+    determines the number of orbitals within the
+    solver block structure.
+
+    Parameters
+    ----------
+    sum_k : dft_tools sumk object
+
+    Returns
+    -------
+    n_orb : dict of int
+        number of orbitals for up / down as dict for SOC calculation
+        without up / down block up holds the number of orbitals
+    """
+    n_orbitals = [{'up': 0, 'down': 0} for i in range(sum_k.n_inequiv_shells)]
+    for icrsh in range(sum_k.n_inequiv_shells):
+        for block, n_orb in sum_k.gf_struct_solver[icrsh].items():
+            if 'down' in block:
+                n_orbitals[icrsh]['down'] += sum_k.gf_struct_solver[icrsh][block]
+            else:
+                n_orbitals[icrsh]['up'] += sum_k.gf_struct_solver[icrsh][block]
+
+    return n_orbitals

python/solid_dmft/dmft_tools/convergence.py

@@ -31,15 +31,15 @@
 # triqs
 from h5 import HDFArchive
 from triqs.gf import MeshImFreq, MeshImTime, MeshReFreq, BlockGf
-from solid_dmft.dmft_tools import solver
+from solid_dmft.dmft_tools import common
 
 def _generate_header(general_params, sum_k):
     """
     Generates the headers that are used in write_header_to_file.
     Returns a dict with {file_name: header_string}
     """
 
-    n_orb = solver.get_n_orbitals(sum_k)
+    n_orb = common.get_n_orbitals(sum_k)
 
     header_energy_mask = '| {:>11} '
     header_energy = header_energy_mask.format('δE_tot')
@@ -77,7 +77,7 @@ def write_conv(conv_obs, sum_k, general_params):
 
     """
 
-    n_orb = solver.get_n_orbitals(sum_k)
+    n_orb = common.get_n_orbitals(sum_k)
 
     for icrsh in range(sum_k.n_inequiv_shells):
         line = '{:3d} | '.format(conv_obs['iteration'][-1])

python/solid_dmft/dmft_tools/initial_self_energies.py

@@ -520,11 +520,13 @@ def determine_dc_and_initial_sigma(general_params, gw_params, advanced_params, s
                                                            mesh=sum_k.mesh)
                            for iineq in range(sum_k.n_inequiv_shells)]
             for icrsh in range(sum_k.n_inequiv_shells):
+                n_orb = sum_k.corr_shells[icrsh]['dim']
                 dc_pot = sum_k.block_structure.convert_matrix(sum_k.dc_imp[sum_k.inequiv_to_corr[icrsh]],
                                                                ish_from=icrsh,
                                                                space_from='sumk', space_to='solver')
 
                 if (general_params['magnetic'] and general_params['magmom'] and sum_k.SO == 0):
+                    mpi.report(f'\n*** Adding magnetic bias to initial sigma and DC for impurity {icrsh} ***')
                     # if we are doing a magnetic calculation and initial magnetic moments
                     # are set, manipulate the initial sigma accordingly
                     fac = general_params['magmom'][icrsh]
@@ -533,9 +535,10 @@ def determine_dc_and_initial_sigma(general_params, gw_params, advanced_params, s
                     # if magmom positive the up channel will be favored
                     for spin_channel in sum_k.gf_struct_solver[icrsh].keys():
                         if 'up' in spin_channel:
-                            start_sigma[icrsh][spin_channel] << -np.eye(dc_pot[spin_channel].shape[0])*fac + dc_pot[spin_channel]
+                            start_sigma[icrsh][spin_channel] <<  dc_pot[spin_channel] - fac*np.eye(n_orb)
                         else:
-                            start_sigma[icrsh][spin_channel] << np.eye(dc_pot[spin_channel].shape[0])*fac + dc_pot[spin_channel]
+                            start_sigma[icrsh][spin_channel] << dc_pot[spin_channel] + fac*np.eye(n_orb)
+
                 else:
                     for spin_channel in sum_k.gf_struct_solver[icrsh].keys():
                         start_sigma[icrsh][spin_channel] << dc_pot[spin_channel]
@@ -545,6 +548,7 @@ def determine_dc_and_initial_sigma(general_params, gw_params, advanced_params, s
             start_sigma = [sum_k.block_structure.create_gf(ish=iineq, gf_function=Gf, space='solver', mesh=sum_k.mesh)
                                 for iineq in range(sum_k.n_inequiv_shells)]
             for icrsh in range(sum_k.n_inequiv_shells):
+                n_orb = sum_k.corr_shells[icrsh]['dim']
                 if (general_params['magnetic'] and general_params['magmom'] and sum_k.SO == 0):
                     mpi.report(f'\n*** Adding magnetic bias to initial sigma for impurity {icrsh} ***')
                     # if we are doing a magnetic calculation and initial magnetic moments
@@ -554,9 +558,9 @@ def determine_dc_and_initial_sigma(general_params, gw_params, advanced_params, s
                     # if magmom positive the up channel will be favored
                     for spin_channel in sum_k.gf_struct_solver[icrsh].keys():
                         if 'up' in spin_channel:
-                            start_sigma[icrsh][spin_channel] << -fac
+                            start_sigma[icrsh][spin_channel] << -fac*np.eye(n_orb)
                         else:
-                            start_sigma[icrsh][spin_channel] << fac
+                            start_sigma[icrsh][spin_channel] << fac*np.eye(n_orb)
         else:
             start_sigma = [sum_k.block_structure.create_gf(ish=iineq, gf_function=Gf, space='solver', mesh=sum_k.mesh)
                            for iineq in range(sum_k.n_inequiv_shells)]

python/solid_dmft/dmft_tools/interaction_hamiltonian.py

@@ -37,7 +37,7 @@
 import triqs.utility.mpi as mpi
 from triqs.gf import make_gf_imfreq
 from triqs.operators import util, n, c, c_dag, Operator
-from solid_dmft.dmft_tools import solver
+from solid_dmft.dmft_tools import common
 
 
 try:
@@ -67,7 +67,7 @@ def _round_to_int(data):
         first_index_shell = 0
         for ish in range(sum_k.n_corr_shells):
             icrsh = sum_k.corr_to_inequiv[ish]
-            n_orb = solver.get_n_orbitals(sum_k)[icrsh]['up']
+            n_orb = common.get_n_orbitals(sum_k)[icrsh]['up']
             u_matrix_temp = u_matrix_four_indices[first_index_shell:first_index_shell+n_orb,
                                                   first_index_shell:first_index_shell+n_orb,
                                                   first_index_shell:first_index_shell+n_orb,
@@ -186,7 +186,7 @@ def _construct_kanamori(sum_k, general_params, solver_type_per_imp, icrsh, den_d
     from the parameters U and J.
     """
 
-    n_orb = solver.get_n_orbitals(sum_k)[icrsh]['up']
+    n_orb = common.get_n_orbitals(sum_k)[icrsh]['up']
     if sum_k.SO == 1:
         assert n_orb % 2 == 0
         n_orb = n_orb // 2
@@ -331,7 +331,7 @@ def _construct_dynamic(sum_k, general_params, icrsh):
             U_onsite = ar['dynamic_U']['U_scr']
     U_onsite = mpi.bcast(U_onsite)
 
-    n_orb = solver.get_n_orbitals(sum_k)[icrsh]['up']
+    n_orb = common.get_n_orbitals(sum_k)[icrsh]['up']
     if sum_k.SO == 1:
         raise ValueError('dynamic U not implemented for SO!=0')
     if n_orb > 1:
@@ -353,7 +353,7 @@ def _generate_four_index_u_matrix(sum_k, general_params, icrsh):
 
     # ish points to the shell representative of the current group
     ish = sum_k.inequiv_to_corr[icrsh]
-    n_orb = solver.get_n_orbitals(sum_k)[icrsh]['up']
+    n_orb = common.get_n_orbitals(sum_k)[icrsh]['up']
     if sum_k.SO == 1:
         assert n_orb % 2 == 0
         n_orb = n_orb // 2
@@ -426,7 +426,7 @@ def _construct_density_density(sum_k, general_params, Umat_full_rotated, icrsh):
     """
 
     # Constructs Hamiltonian from Umat_full_rotated
-    n_orb = solver.get_n_orbitals(sum_k)[icrsh]['up']
+    n_orb = common.get_n_orbitals(sum_k)[icrsh]['up']
 
     Umat, Upmat = util.reduce_4index_to_2index(Umat_full_rotated)
     h_int = util.h_int_density(sum_k.spin_block_names[sum_k.SO], n_orb,
@@ -442,7 +442,7 @@ def _construct_slater(sum_k, general_params, Umat_full_rotated, icrsh):
     matrix.
     """
 
-    n_orb = solver.get_n_orbitals(sum_k)[icrsh]['up']
+    n_orb = common.get_n_orbitals(sum_k)[icrsh]['up']
 
     h_int = util.h_int_slater(sum_k.spin_block_names[sum_k.SO], n_orb,
                               map_operator_structure=sum_k.sumk_to_solver[icrsh],
@@ -589,7 +589,7 @@ def construct(sum_k, general_params, solver_type_per_imp,  gw_params=None):
 
         if general_params['h_int_type'][icrsh] == 'simple_intra':
             h_int[icrsh] = h_int_simple_intra(sum_k.spin_block_names[sum_k.SO],
-                                              solver.get_n_orbitals(sum_k)[icrsh]['up'],
+                                              common.get_n_orbitals(sum_k)[icrsh]['up'],
                                               map_operator_structure=sum_k.sumk_to_solver[icrsh],
                                               U=general_params['U'][icrsh],
                                               H_dump=os.path.join(general_params['jobname'], f'H_imp{icrsh}.txt'))

python/solid_dmft/dmft_tools/observables.py

@@ -35,7 +35,7 @@
 from triqs.gf import Gf, MeshImTime
 from triqs.atom_diag import trace_rho_op
 from triqs.gf.descriptors import Fourier
-from solid_dmft.dmft_tools import solver
+from solid_dmft.dmft_tools import common
 
 def prep_observables(h5_archive, sum_k):
     """
@@ -93,7 +93,7 @@ def _generate_header(general_params, sum_k):
     Generates the headers that are used in write_header_to_file.
     Returns a dict with {file_name: header_string}
     """
-    n_orb = solver.get_n_orbitals(sum_k)
+    n_orb = common.get_n_orbitals(sum_k)
 
     header_energy_mask = ' | {:>10} | {:>10}   {:>10}   {:>10}   {:>10}'
     header_energy = header_energy_mask.format('E_tot', 'E_DFT', 'E_bandcorr', 'E_int_imp', 'E_DC')
@@ -416,7 +416,7 @@ def write_obs(observables, sum_k, general_params):
 
     """
 
-    n_orb = solver.get_n_orbitals(sum_k)
+    n_orb = common.get_n_orbitals(sum_k)
 
     for icrsh in range(sum_k.n_inequiv_shells):
         if general_params['magnetic'] and sum_k.SO == 0:

python/solid_dmft/dmft_tools/results_to_archive.py

@@ -89,7 +89,7 @@ def _compile_information(sum_k, general_params, solver_params, solvers, map_imp_
                 write_to_h5['G_time_dlr_{}'.format(icrsh)] = solvers[icrsh].G_time_dlr
                 write_to_h5['Sigma_dlr_{}'.format(icrsh)] = solvers[icrsh].Sigma_dlr
 
-        if solver_type_per_imp[icrsh] == 'ctint' and solver_params[isolvsec]['measure_histogram']:
+        if solver_type_per_imp[icrsh] == 'ctint' and solver_params[isolvsec]['measure_pert_order']:
             write_to_h5['pert_order_imp_{}'.format(icrsh)] = solvers[icrsh].perturbation_order
 
         if solver_type_per_imp[icrsh] == 'hubbardI':