First chunk of changes. test lno_hubbardI_mag runs through

python/solid_dmft/dmft_tools/afm_mapping.py

@@ -120,7 +120,7 @@ def apply(general_params, solver_params, icrsh, gf_struct_solver, solvers):
         solvers[icrsh].density_matrix = solvers[imp_source].density_matrix
         solvers[icrsh].h_loc_diagonalization = solvers[imp_source].h_loc_diagonalization
 
-    if general_params['measure_chi'] != 'none':
+    if general_params['measure_chi'] is not None:
         solvers[icrsh].O_time = solvers[imp_source].O_time
 
     return solvers

python/solid_dmft/dmft_tools/convergence.py

@@ -169,15 +169,14 @@ def max_G_diff(G1, G2, norm_temp = True):
 
     return norm
 
-def prep_conv_obs(h5_archive, sum_k):
+def prep_conv_obs(h5_archive):
     """
     prepares the conv arrays and files for the DMFT calculation
 
     Parameters
     ----------
     h5_archive: hdf archive instance
         hdf archive for calculation
-    sum_k : SumK Object instances
 
     __Returns:__
     conv_obs : dict

python/solid_dmft/dmft_tools/formatter.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 ################################################################################
 #
 # solid_dmft - A versatile python wrapper to perform DFT+DMFT calculations
@@ -90,13 +89,10 @@ def print_block_sym(sum_k, dm, general_params):
     if dm:
         # Prints dft density matrix
         print('\nDFT density matrix')
+        # Double loop to sort by impurities
         for icrsh in range(sum_k.n_inequiv_shells):
             spins = sum_k.spin_block_names[sum_k.SO]
-            # ftps dens_mat only knows inequiv_shells
-            if general_params['solver_type'] in ['ftps']:
-                shlst = [icrsh]
-            else:
-                shlst = [ish for ish, ineq_shell in enumerate(sum_k.corr_to_inequiv) if ineq_shell == icrsh]
+            shlst = [ish for ish, ineq_shell in enumerate(sum_k.corr_to_inequiv) if ineq_shell == icrsh]
             for sh in shlst:
                 n_orb = sum_k.corr_shells[sh]['dim']
                 if sum_k.SO == 0:

python/solid_dmft/dmft_tools/initial_self_energies.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 ################################################################################
 #
 # solid_dmft - A versatile python wrapper to perform DFT+DMFT calculations
@@ -38,7 +37,7 @@
 from triqs.gf import BlockGf, Gf
 
 
-def calculate_double_counting(sum_k, density_matrix, general_params, advanced_params):
+def calculate_double_counting(sum_k, density_matrix, general_params, advanced_params, solver_type_per_imp):
     """
     Calculates the double counting, including all manipulations from advanced_params.
 
@@ -63,34 +62,39 @@ def calculate_double_counting(sum_k, density_matrix, general_params, advanced_pa
     # copy the density matrix to not change it
     density_matrix_DC = deepcopy(density_matrix)
 
-    if general_params['solver_type'] == 'hartree':
-        mpi.report('\nSOLID_DMFT: Hartree solver detected, zeroing out the DC correction. This gets computed at the solver level')
-        for icrsh in range(sum_k.n_inequiv_shells):
-            sum_k.calc_dc(density_matrix_DC[icrsh], orb=icrsh,
-                          use_dc_value=0.0)
-        return sum_k
+    def iterate_except_hartree():
+        for iineq, type in enumerate(solver_type_per_imp):
+            if type == 'hartree':
+                mpi.report('\nSOLID_DMFT: Hartree solver for impurity {iineq} detected. '
+                           'Zeroing out the DC correction. This gets computed at the solver level')
+            else:
+                yield iineq
+
+    # TODO: suppress print when reseting DC to zero
+    for icrsh in range(sum_k.n_inequiv_shells):
+        sum_k.calc_dc(density_matrix_DC[icrsh], orb=icrsh,
+                      use_dc_value=0.0)
 
     # Sets the DC and exits the function if advanced_params['dc_fixed_value'] is specified
-    if advanced_params['dc_fixed_value'] != 'none':
-        for icrsh in range(sum_k.n_inequiv_shells):
+    if advanced_params['dc_fixed_value'] is not None:
+        for icrsh in iterate_except_hartree():
             sum_k.calc_dc(density_matrix_DC[icrsh], orb=icrsh,
                           use_dc_value=advanced_params['dc_fixed_value'])
         return sum_k
 
     # use DC for CPA
     if general_params['dc'] and general_params['dc_type'] == 4:
         zetas = general_params['cpa_zeta']
-        for icrsh in range(sum_k.n_inequiv_shells):
+        for icrsh in iterate_except_hartree():
             sum_k.calc_dc(density_matrix_DC[icrsh], orb=icrsh, use_dc_value=zetas[icrsh])
 
         return sum_k
 
-
-    if advanced_params['dc_fixed_occ'] != 'none':
+    if advanced_params['dc_fixed_occ'] is not None:
         mpi.report('Fixing occupation for DC potential to provided value')
 
         assert sum_k.n_inequiv_shells == len(advanced_params['dc_fixed_occ']), "give exactly one occupation per correlated shell"
-        for icrsh in range(sum_k.n_inequiv_shells):
+        for icrsh in iterate_except_hartree():
             mpi.report('fixing occupation for impurity '+str(icrsh)+' to n='+str(advanced_params['dc_fixed_occ'][icrsh]))
             n_orb = sum_k.corr_shells[icrsh]['dim']
             # we need to handover a matrix to calc_dc so calc occ per orb per spin channel
@@ -100,7 +104,7 @@ def calculate_double_counting(sum_k, density_matrix, general_params, advanced_pa
                 np.fill_diagonal(inner, orb_occ+0.0j)
 
     # The regular way: calculates the DC based on U, J and the dc_type
-    for icrsh in range(sum_k.n_inequiv_shells):
+    for icrsh in iterate_except_hartree():
         if general_params['dc_type'] == 3:
             # this is FLL for eg orbitals only as done in Seth PRB 96 205139 2017 eq 10
             # this setting for U and J is reasonable as it is in the spirit of F0 and Javg
@@ -118,6 +122,8 @@ def calculate_double_counting(sum_k, density_matrix, general_params, advanced_pa
     # for the fixed DC according to https://doi.org/10.1103/PhysRevB.90.075136
     # dc_imp is calculated with fixed occ but dc_energ is calculated with given n
     if advanced_params['dc_nominal']:
+        if 'Hartree' in solver_type_per_imp:
+            raise NotImplementedError('dc_nominal not implemented in presence of Hartree solver')
         mpi.report('\ncalculating DC energy with fixed DC potential from above\n'
                    + ' for the original density matrix doi.org/10.1103/PhysRevB.90.075136\n'
                    + ' aka nominal DC')
@@ -147,7 +153,9 @@ def calculate_double_counting(sum_k, density_matrix, general_params, advanced_pa
             mpi.report('DC energy for shell {} = {}'.format(icrsh, energy_per_shell))
 
     # Rescales DC if advanced_params['dc_factor'] is given
-    if advanced_params['dc_factor'] != 'none':
+    if advanced_params['dc_factor'] is not None:
+        if 'Hartree' in solver_type_per_imp:
+            raise NotImplementedError('dc_factor not implemented in presence of Hartree solver')
         rescaled_dc_imp = [{spin: advanced_params['dc_factor'] * dc_per_spin
                             for spin, dc_per_spin in dc_per_shell.items()}
                           for dc_per_shell in sum_k.dc_imp]
@@ -162,7 +170,9 @@ def calculate_double_counting(sum_k, density_matrix, general_params, advanced_pa
                 mpi.report('DC for shell {} and block {} = {}'.format(icrsh, spin, dc_per_spin[0][0]))
             mpi.report('DC energy for shell {} = {}'.format(icrsh, energy_per_shell))
 
-    if advanced_params['dc_orb_shift'] != 'none':
+    if advanced_params['dc_orb_shift'] is not None:
+        if 'Hartree' in solver_type_per_imp:
+            raise NotImplementedError('dc_orb_shift not implemented in presence of Hartree solver')
         mpi.report('adding an extra orbital dependent shift per impurity')
         tot_norb = 0
         dc_orb_shift = []
@@ -350,7 +360,8 @@ def _set_loaded_sigma(sum_k, loaded_sigma, loaded_dc_imp, general_params):
 
 
 def determine_dc_and_initial_sigma(general_params, advanced_params, sum_k,
-                                   archive, iteration_offset, density_mat_dft, solvers):
+                                   archive, iteration_offset, density_mat_dft, solvers,
+                                   solver_type_per_imp):
     """
     Determines the double counting (DC) and the initial Sigma. This can happen
     in five different ways:
@@ -399,7 +410,7 @@ def determine_dc_and_initial_sigma(general_params, advanced_params, sum_k,
             print('\nFrom previous calculation:', end=' ')
             start_sigma, sum_k.dc_imp, sum_k.dc_energ, last_g0,  _ = _load_sigma_from_h5(archive, -1)
             if general_params['csc'] and not general_params['dc_dmft']:
-                sum_k = calculate_double_counting(sum_k, density_mat_dft, general_params, advanced_params)
+                sum_k = calculate_double_counting(sum_k, density_mat_dft, general_params, advanced_params, solver_type_per_imp)
         # Loads Sigma from different calculation
         elif general_params['load_sigma']:
             print('\nFrom {}:'.format(general_params['path_to_sigma']), end=' ')
@@ -411,25 +422,24 @@ def determine_dc_and_initial_sigma(general_params, advanced_params, sum_k,
             if general_params['dc']:
                 if general_params['dc_dmft']:
                     sum_k = calculate_double_counting(sum_k, loaded_density_matrix,
-                                                       general_params, advanced_params)
+                                                      general_params, advanced_params,
+                                                      solver_type_per_imp)
                 else:
                     sum_k = calculate_double_counting(sum_k, density_mat_dft,
-                                                       general_params, advanced_params)
+                                                      general_params, advanced_params,
+                                                      solver_type_per_imp)
 
             start_sigma = _set_loaded_sigma(sum_k, loaded_sigma, loaded_dc_imp, general_params)
 
         # Sets DC as Sigma because no initial Sigma given
         elif general_params['dc']:
-            sum_k = calculate_double_counting(sum_k, density_mat_dft, general_params, advanced_params)
+            sum_k = calculate_double_counting(sum_k, density_mat_dft, general_params, advanced_params,
+                                              solver_type_per_imp)
 
             # initialize Sigma from sum_k
-            if general_params['solver_type'] in ['ftps']:
-                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)]
-            else:
-                start_sigma = [sum_k.block_structure.create_gf(ish=iineq, space='solver', mesh=sum_k.mesh)
-                               for iineq in range(sum_k.n_inequiv_shells)]
+            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):
                 dc_value = sum_k.dc_imp[sum_k.inequiv_to_corr[icrsh]][sum_k.spin_block_names[sum_k.SO][0]][0, 0]
 
@@ -466,13 +476,8 @@ def determine_dc_and_initial_sigma(general_params, advanced_params, sum_k,
                         else:
                             start_sigma[icrsh][spin_channel] << fac
         else:
-            if general_params['solver_type'] in ['ftps']:
-                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)]
-            else:
-                start_sigma = [sum_k.block_structure.create_gf(ish=iineq, space='solver', mesh=sum_k.mesh)
-                               for iineq in range(sum_k.n_inequiv_shells)]
+            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)]
 
     # Adds random, frequency-independent noise in zeroth iteration to break symmetries
     if not np.isclose(general_params['noise_level_initial_sigma'], 0) and iteration_offset == 0:
@@ -500,11 +505,11 @@ def determine_dc_and_initial_sigma(general_params, advanced_params, sum_k,
     sum_k.put_Sigma([solvers[icrsh].Sigma_freq for icrsh in range(sum_k.n_inequiv_shells)])
 
     # load sigma as first guess in the hartree solver if applicable
-    if general_params['solver_type'] == 'hartree':
+    for icrsh in range(sum_k.n_inequiv_shells):
         # TODO:
         # should this be moved to before the solve() call? Having it only here means there is a mismatch
         # between the mixing at the level of the solver and the sumk (solver mixes always 100%)
-        for icrsh in range(sum_k.n_inequiv_shells):
+        if solver_type_per_imp[icrsh] == 'hartree':
             mpi.report(f"SOLID_DMFT: setting first guess hartree solver for impurity {icrsh}")
             solvers[icrsh].triqs_solver.reinitialize_sigma(start_sigma[icrsh])
 

python/solid_dmft/dmft_tools/interaction_hamiltonian.py

@@ -42,33 +42,6 @@
 except ImportError:
     pass
 
-
-def _extract_U_J_list(param_name, n_inequiv_shells, general_params):
-    """
-    Checks if param_name ('U', 'U_prime', or 'J') are a single value or
-    different per inequivalent shell. If just a single value is given,
-    this value is applied to each shell.
-    """
-
-    if not isinstance(param_name, str):
-        formatted_param = ['none' if p == 'none' else '{:.2f}'.format(p)
-                           for p in general_params[param_name]]
-    else:
-        formatted_param = general_params[param_name]
-
-    if len(general_params[param_name]) == 1:
-        mpi.report('Assuming {} = '.format(param_name)
-                   + '{} for all correlated shells'.format(formatted_param[0]))
-        general_params[param_name] *= n_inequiv_shells
-    elif len(general_params[param_name]) == n_inequiv_shells:
-        mpi.report('{} list for correlated shells: {}'.format(param_name, formatted_param))
-    else:
-        raise IndexError('Property list {} '.format(general_params[param_name])
-                         + 'must have length 1 or n_inequiv_shells')
-
-    return general_params
-
-
 def _load_crpa_interaction_matrix(sum_k, filename='UIJKL'):
     """
     Loads VASP cRPA data to use as an interaction Hamiltonian.
@@ -173,7 +146,7 @@ def _adapt_U_4index_for_SO(Umat_full):
     return Umat_full_SO
 
 
-def _construct_kanamori(sum_k, general_params, icrsh):
+def _construct_kanamori(sum_k, general_params, solver_type_per_imp, icrsh):
     """
     Constructs the Kanamori interaction Hamiltonian. Only Kanamori does not
     need the full four-index matrix. Therefore, we can construct it directly
@@ -195,7 +168,7 @@ def _construct_kanamori(sum_k, general_params, icrsh):
     else:
         U_prime = general_params['U_prime'][icrsh]
 
-    if general_params['solver_type'] == 'ftps':
+    if solver_type_per_imp[icrsh] == 'ftps':
         # 1-band modell requires J and U' equals zero
         if n_orb == 1:
             up, j = 0.0, 0.0
@@ -355,7 +328,7 @@ def _generate_four_index_u_matrix(sum_k, general_params, icrsh):
         U = general_params['U'][icrsh]
         J = general_params['J'][icrsh]
         R = general_params['ratio_F4_F2'][icrsh]
-        R = 0.63 if R == 'none' else R
+        R = 0.63 if R is None else R
         slater_integrals = np.array([U, 14*J/(1+R), 14*J*R/(1+R)])
 
     mpi.report('\nImpurity {}: The corresponding slater integrals are'.format(icrsh))
@@ -497,7 +470,7 @@ def h_int_simple_intra(spin_names,n_orb,U,off_diag=None,map_operator_structure=N
     return H
 
 
-def construct(sum_k, general_params, advanced_params):
+def construct(sum_k, general_params, solver_type_per_imp):
     """
     Constructs the interaction Hamiltonian. Currently implemented are the
     Kanamori Hamiltonian (usually for 2 or 3 orbitals), the density-density and
@@ -522,18 +495,6 @@ def construct(sum_k, general_params, advanced_params):
     # Extracts U and J
     mpi.report('*** interaction parameters ***')
 
-    general_params = _extract_U_J_list('h_int_type', sum_k.n_inequiv_shells, general_params)
-    for param_name in ('U', 'J'):
-        general_params = _extract_U_J_list(param_name, sum_k.n_inequiv_shells, general_params)
-    if 'kanamori' in general_params['h_int_type']:
-        general_params = _extract_U_J_list('U_prime', sum_k.n_inequiv_shells, general_params)
-    for param_name in ('dc_U', 'dc_J'):
-        advanced_params = _extract_U_J_list(param_name, sum_k.n_inequiv_shells, advanced_params)
-
-    # Extracts ratio_F4_F2 if any shell uses a solver supporting it
-    if 'density_density' in general_params['h_int_type'] or 'full_slater' in general_params['h_int_type']:
-        general_params = _extract_U_J_list('ratio_F4_F2', sum_k.n_inequiv_shells, general_params)
-
     # Constructs the interaction Hamiltonian. Needs to come after setting sum_k.rot_mat
     mpi.report('\nConstructing the interaction Hamiltonians')
     h_int = [None] * sum_k.n_inequiv_shells
@@ -542,22 +503,22 @@ def construct(sum_k, general_params, advanced_params):
 
         # Kanamori
         if general_params['h_int_type'][icrsh] == 'kanamori':
-            h_int[icrsh] = _construct_kanamori(sum_k, general_params, icrsh)
+            h_int[icrsh] = _construct_kanamori(sum_k, general_params, solver_type_per_imp, icrsh)
             continue
 
         # for density density or full slater get full four-index U matrix
         if general_params['h_int_type'][icrsh] in ('density_density', 'full_slater'):
             mpi.report('\nNote: The input parameters U and J here are orbital-averaged parameters.')
             mpi.report('Note: The order of the orbitals is important. See also the doc string of this method.')
-            # Checks that all entries are l == 2 or R == 'none'
+            # Checks that all entries are l == 2 or R is None
             if (sum_k.corr_shells[sum_k.inequiv_to_corr[icrsh]]['l'] != 2
-                    and general_params['ratio_F4_F2'][icrsh] != 'none'):
+                    and general_params['ratio_F4_F2'][icrsh] is not None):
                 raise ValueError('Ratio F4/F2 only implemented for d-shells '
                                  + 'but set in impurity {}'.format(icrsh))
 
-            if general_params['h_int_type'][icrsh] == 'density_density' and general_params['solver_type'] == 'ftps':
+            if general_params['h_int_type'][icrsh] == 'density_density' and solver_type_per_imp[icrsh] == 'ftps':
                 # TODO: implement
-                raise NotImplementedError('\nNote: Density-density not implemented for ftps.')
+                raise NotImplementedError('Note: Density-density not implemented for ftps.')
 
             Umat_full = _generate_four_index_u_matrix(sum_k, general_params, icrsh)