diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py index 0a509415..df47a248 100644 --- a/python/triqs_dft_tools/sumk_dft.py +++ b/python/triqs_dft_tools/sumk_dft.py @@ -1955,22 +1955,23 @@ class SumkDFT(object): """ def find_bounds(function, x_init, delta_x, max_loops=1000): - mpi.report("Finding bounds on chemical potential" ) - x= x_init + mpi.report("Finding bounds on chemical potential") + x = x_init # First find the bounds y1 = function(x) eps = np.sign(y1) - x1= x; - x2= x1;y2 = y1 - - nbre_loop=0 - #abort the loop after maxiter is reached or when y1 and y2 have different sign - while (nbre_loop<= max_loops) and (y2*y1)>0: - nbre_loop +=1 - x1=x2 - y1=y2 + x1 = x + x2 = x1 + y2 = y1 - x2 -= eps*delta_x + nbre_loop = 0 + # abort the loop after maxiter is reached or when y1 and y2 have different sign + while (nbre_loop <= max_loops) and (y2*y1) > 0: + nbre_loop += 1 + x1 = x2 + y1 = y2 + + x2 -= eps*delta_x y2 = function(x2) if nbre_loop > (max_loops): @@ -1978,68 +1979,67 @@ class SumkDFT(object): # Make sure that x2 > x1 if x1 > x2: - x1,x2 = x2,x1 - y1,y2 = y2,y1 - + x1, x2 = x2, x1 + y1, y2 = y2, y1 mpi.report(f"mu_interval: [ {x1:.4f} ; {x2:.4f} ]") mpi.report(f"delta to target density interval: [ {y1:.4f} ; {y2:.4f} ]") return x1, x2 - - # previous implementation + def F_bisection(mu): return self.total_density(mu=mu, broadening=broadening).real density = self.density_required - self.charge_below - #using scipy.optimize + # using scipy.optimize + def F_optimize(mu): mpi.report("Trying out mu = {}".format(str(mu))) calc_dens = self.total_density(mu=mu, broadening=broadening).real - density mpi.report(f"Target density = {density}; Delta to target = {calc_dens}") return calc_dens - - #check for lowercase matching for the method variable - if method.lower()=="dichotomy": - mpi.report("\nSUMK calc_mu: Using dichtomy adjustment to find chemical potential\n") + + # check for lowercase matching for the method variable + if method.lower() == "dichotomy": + mpi.report("\nsumk calc_mu: Using dichtomy adjustment to find chemical potential\n") self.chemical_potential = dichotomy.dichotomy(function=F_bisection, x_init=self.chemical_potential, y_value=density, precision_on_y=precision, delta_x=delta, max_loops=max_loops, x_name="Chemical Potential", y_name="Total Density", verbosity=3)[0] - elif method.lower()=="newton": - mpi.report("\nSUMK calc_mu: Using Newton method to find chemical potential\n") - self.chemical_potential = newton(func=F_optimize, - x0=self.chemical_potential, - tol=precision, maxiter=max_loops, - ) + elif method.lower() == "newton": + mpi.report("\nsumk calc_mu: Using Newton method to find chemical potential\n") + self.chemical_potential = newton(func=F_optimize, + x0=self.chemical_potential, + tol=precision, maxiter=max_loops, + ) - elif method.lower()=="brent": - mpi.report("\nSUMK calc_mu: Using Brent method to find chemical potential") - mpi.report("SUMK calc_mu: Finding bounds \n") - - mu_guess_0, mu_guess_1 = find_bounds(function=F_optimize, - x_init=self.chemical_potential, - delta_x=delta, max_loops=max_loops, - ) - mu_guess_1 += 0.01 #scrambles higher lying interval to avoid getting stuck - mpi.report("\nSUMK calc_mu: Searching root with Brent method\n") + elif method.lower() == "brent": + mpi.report("\nsumk calc_mu: Using Brent method to find chemical potential") + mpi.report("sumk calc_mu: Finding bounds \n") + + mu_guess_0, mu_guess_1 = find_bounds(function=F_optimize, + x_init=self.chemical_potential, + delta_x=delta, max_loops=max_loops, + ) + mu_guess_1 += 0.01 # scrambles higher lying interval to avoid getting stuck + mpi.report("\nsumk calc_mu: Searching root with Brent method\n") self.chemical_potential = brenth(f=F_optimize, - a=mu_guess_0, - b=mu_guess_1, - xtol=precision, maxiter=max_loops, - ) - + a=mu_guess_0, + b=mu_guess_1, + xtol=precision, maxiter=max_loops, + ) + else: raise ValueError( - f"SUMK calc_mu: The selected method: {method}, is not implemented\n", - """ + f"sumk calc_mu: The selected method: {method}, is not implemented\n", + """ Please check for typos or select one of the following: * dichotomy: usual bisection algorithm from the TRIQS library * newton: newton method, fastest convergence but more unstable * brent: finds bounds and proceeds with hyperbolic brent method, a compromise between speed and ensuring convergence """ - ) + ) return self.chemical_potential diff --git a/test/python/calc_mu.py b/test/python/calc_mu.py new file mode 100644 index 00000000..da7f2d3a --- /dev/null +++ b/test/python/calc_mu.py @@ -0,0 +1,56 @@ +########################################################################## +# +# TRIQS: a Toolbox for Research in Interacting Quantum Systems +# +# Copyright (C) 2023 by A. Hampel +# +# TRIQS is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. +# +# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# TRIQS. If not, see . +# +########################################################################## + +import unittest +import numpy as np + +from h5 import * +from triqs.gf import MeshImFreq + +from triqs_dft_tools.sumk_dft import * + + +class test_solver(unittest.TestCase): + + def setUp(self): + self.iw_mesh = MeshImFreq(beta=40, S='Fermion', n_iw=300) + # magic reference value for the Wien2k SVO t2g example + self.ref_mu = 0.281 + + def test_dichotomy(self): + sumk = SumkDFT('SrVO3.ref.h5', mesh=self.iw_mesh) + mu = sumk.calc_mu(method='dichotomy', precision=0.001, delta=0.1) + self.assertTrue(abs(self.ref_mu - mu) < 0.01) + + def test_brent(self): + sumk = SumkDFT('SrVO3.ref.h5', mesh=self.iw_mesh) + mu = sumk.calc_mu(method='brent', precision=0.001, delta=0.1) + self.assertTrue(abs(self.ref_mu - mu) < 0.01) + + # Newton seems to be quite unstable here. Not tested right now! + # def test_newton(self): + # sumk = SumkDFT('SrVO3.ref.h5', mesh = self.iw_mesh) + # mu = sumk.calc_mu(method='newton', precision= 0.1, delta=0.01) + # self.assertTrue(abs(self.ref_mu - mu) < 0.01) + + +if __name__ == '__main__': + unittest.main()