Source code for contactpp.utp

import numpy as np
from scipy.optimize import fmin
from numpy.polynomial.polynomial import Polynomial

from troullier import make_troullier_potential
from solve import ShootScattering

class UTPPotential(object):
[docs] """ Object representing a scattering potential. """ def __init__(self,potential_polynomial,cutoff,scattering_length): self._potential_polynomial = potential_polynomial self.cutoff = cutoff self._V = lambda r: np.where(abs(r)<=self.cutoff, self._potential_polynomial(r/self.cutoff),0.) @property def V(self): return self._V @property def coefficients(self): return self._potential_polynomial.coef class UTPGenerator(object):
def __init__(self,a,kf,c=None,init_coeffs=None, integrator="vode",npoints=1000,nks=10,objective_function="rms",dos=None, verbose=True): self.a = a self.kf = float(kf) self.Ef = self.kf**2 self.c = c if self.a < 0.: if self.c is None: raise ValueError("Cutoff cannot be undefined if the \ branch is 'attractive'.") if init_coeffs is not None: if c is None: raise ValueError( "If initial coefficients are specified, the cutoff\ must also be specified explicitly.") self.init_coeffs = init_coeffs else: self.init_coeffs, self.c = \ self._get_init_coeffs_from_troullier() self.integrator = integrator self.npoints = npoints self.objective_function = objective_function if self.objective_function not in ("rms", "max") : raise ValueError( "objective_function argument should be either 'rms'\ or 'max'.") if self.objective_function == "max" and dos is not None: raise ValueError( "Cannot specify DOS when objective function is 'max'. "\ "Specifying a DOS is only valid when the objective "\ "function is 'rms'.") if self.objective_function == "rms" : self.dos = dos if dos is not None else lambda k: k**2 self.OBJECTIVE_FUNCTIONS = { "rms" : self._calc_rms_error, "max" : self._calc_max_error } self.ks = np.linspace(0.,kf,nks) self.verbose = bool(verbose) def _get_init_coeffs_from_troullier(self): Ef = self.kf**2 calibration_energy = 0.6*Ef troullier_pp = make_troullier_potential( "repulsive" if self.a > 0 else "attractive", self.a,calibration_energy,self.c) assert(troullier_pp.coefficients[1] == troullier_pp.coefficients[3] == 0.) troullier_coeffs = troullier_pp.coefficients # convert from polynomial in r to polynomial in r/rc troullier_poly = Polynomial( [ coeff*troullier_pp.cutoff**icoeff for icoeff, coeff in enumerate(troullier_coeffs) ]) # factor out (1-r/r_c)^2 init_poly = troullier_poly / (Polynomial([1.,-1.])**2) init_poly_coeffs = init_poly.coef init_coeffs = np.zeros((len(init_poly_coeffs)-2,)) init_coeffs[0] = init_poly_coeffs[1] # v1 init_coeffs[1] = init_poly_coeffs[2] # v2 init_coeffs[2:] = init_poly_coeffs[4:] return init_coeffs,troullier_pp.cutoff def _calc_tandelta(self,V,k): s = ShootScattering(V,integrator=self.integrator,mass=0.5) E = k**2 delta = s.get_phase_shift(E,1.1*self.c,npts=self.npoints) return float(np.tan(delta)) def _calc_true_tandelta(self,k): return -k*self.a def _calc_tandelta_error(self,V): tandeltas = np.array([ self._calc_tandelta(V,k) for k in self.ks ]) true_tandeltas = np.array([ self._calc_true_tandelta(k) for k in self.ks ]) return true_tandeltas - tandeltas def _calc_rms_error(self,V): """ Objective function. RMS error over all k values, weighted by a density of states. """ tandelta_errors = self._calc_tandelta_error(V) return np.sqrt(sum( (tandelta_errors)**2 * self.dos(self.ks)))/len(self.ks) def _calc_max_error(self,V): """ Objective function. Maximal (L-infinity norm) error over all k-values. """ tandelta_errors = self._calc_tandelta_error(V) return max(abs(tandelta_errors)) def _coeffs2poly(self,coeffs): poly_coeffs = np.zeros((len(self.init_coeffs)+2,)) poly_coeffs[0] = 0.5*coeffs[0] # v1 poly_coeffs[1] = coeffs[0] # v1 poly_coeffs[2] = coeffs[1] # v2 poly_coeffs[3] = 2*coeffs[1]-coeffs[0] poly_coeffs[4:] = coeffs[2:] p = Polynomial(poly_coeffs)*Polynomial([1,-1])**2 assert abs(p.coef[3]) < 1e-10 assert abs(p.coef[1]) < 1e-10 assert abs(p(1.))<1e-10 return p def _make_pseudo(self,coeffs): p = self._coeffs2poly(coeffs) return lambda r: np.where(r<self.c,p(r/self.c),0.) def _make_objective_func(self): f = self.OBJECTIVE_FUNCTIONS[self.objective_function] def objective(coeffs): V = self._make_pseudo(coeffs) error = f(V) return error return objective def _make_callback_func(self): self.niter = 0 if self.verbose: def callback(coeffs): if (self.niter % 10 == 0): V = self._make_pseudo(coeffs) error = self.OBJECTIVE_FUNCTIONS[self.objective_function](V) print "# iterations: {:<4d} | {} error: {}".format( self.niter,self.objective_function, error) self.niter += 1 else: callback = None return callback def _run_optimizer(self): new_coeffs = fmin(self._make_objective_func(), self.init_coeffs, callback=self._make_callback_func(), disp=self.verbose) return new_coeffs def make_pseudopotential(self): new_coeffs = self._run_optimizer() p = self._coeffs2poly(new_coeffs) return UTPPotential(p,self.c,self.a) def make_utp_potential(scattering_length, fermi_energy, cutoff=None, init_coeffs=None,integrator="vode",npoints=1000,nks=10,objective_function="rms", dos=None,verbose=True): """ Create a UTPPotential object. Parameters ---------- scattering_length : float Scattering length. If positive, the pseudopotential is built for the repulsive branch. If negative, it is built for the attractive branch. Note that UTP potentials do not exist for the bound state. fermi_energy : float The UTP constructor will minimize the errors in phase shifts over the range ``k = 0`` to ``kf``, where ``kf = sqrt(Ef)`` is the Fermi wavevector. cutoff : float Choice of cutoff radius. Compulsory for negative scattering lengths. If the branch is "repulsive", the cutoff is chosen as the first antinode of the true wavefunction at ``3/5 * fermi_energy``. Specifying the cutoff explicitly will override this. init_coeffs : iterable, optional Initial coefficients for the optimizer. If not specified, a Troullier pseudopotential is constructed and those coefficients are used. If specified, the coefficients will be interpreted as: ``[ v1, v2, v4, v5 ... v11 ]``, with the potential defined as, ``V(x)/Ef = (1-x^2) * (v1*(0.5+x) + v2*x^2 + (2v2-v1)*x^3 + v4*x^4 + sum_{i=5}^13 v_i x^i`` where ``x = r/cutoff``. integrator : string, optional Which integrator to use when computing the phase shift. See the documentation for ``scipy.integrate.ode`` for details. If the cutoff chosen is very small, it might be worthwhile changing this to ``"dopri5"``. ``"vode"`` by default. npoints : int, optional Number of points to use in the integration. If the cutoff chosen is small, convergence of the results with respect to ``npoints`` should be checked. 1000 by default. nks : int, optional Number of k-points to use when integrating the error in phase shifts. Increasing this may lead to somewhat better results. 10 by default. objective_function : string, optional Type of objective function to minimize: either "rms" to minimize ``int_0^kF { [delta_PP(k) - delta_true(k)]^2 dos(k) dk``, where dos is another optional argument, or "max" to minimize ``max(|delta_PP(k) - delta_true(k)|)``. "rms" by default. dos : function, dos(k) -> density of states, optional When using the RMS objective function, this argument sets the cost function that is optimized when constructing the potential. The objective function is: ``int_0^kF { [delta_PP(k) - delta_true(k)]^2 dos(k) dk``. By default, ``dos(k) = lambda k: k**2``. Must accept a numpy array (of k-values) as input. If the objective_function is not RMS, specifying this argument raises a ValueError. verbose : boolean, optional Display convergence information. True by default. """ pgen = UTPGenerator(scattering_length,np.sqrt(fermi_energy), cutoff,init_coeffs,integrator,npoints,nks,objective_function, dos,verbose) return pgen.make_pseudopotential()