fermion_nonrel.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_NONREL_FERMION_H
24 #define O2SCL_NONREL_FERMION_H
25 
26 /** \file fermion_nonrel.h
27  \brief File defining \ref o2scl::fermion_nonrel_tl
28 */
29 
30 #include <string>
31 #include <iostream>
32 #include <fstream>
33 #include <cmath>
34 #include <o2scl/constants.h>
35 #include <o2scl/mroot.h>
36 #include <o2scl/inte.h>
37 #include <o2scl/root_cern.h>
38 #include <o2scl/root_brent_gsl.h>
39 #include <o2scl/classical.h>
40 #include <o2scl/inte_qagiu_gsl.h>
41 
42 #include <o2scl/fermion.h>
43 #include <o2scl/polylog.h>
44 
45 #ifndef DOXYGEN_NO_O2NS
46 namespace o2scl {
47 #endif
48 
49  /** \brief Nonrelativistic fermion class
50 
51  The effective mass computed by this class and stored in \ref
52  part_tl::ms is the Landau mass, not the Dirac mass, as computed by
53  \ref o2scl::fermion_rel_tl .
54 
55  This class works with both true and false values for either \ref
56  part_tl::non_interacting or \ref part_tl::inc_rest_mass.
57 
58  Pressure is computed with
59  \f[
60  P = 2 \varepsilon/3
61  \f]
62  and entropy density with
63  \f[
64  s = \frac{5 \varepsilon}{3 T} - \frac{n \mu}{T}
65  \f]
66  These relations can be verified with an integration by
67  parts. See, e.g. \ref Callen85 pg. 403 or \ref Landau80 pg. 164.
68 
69  The functions \ref pair_density() and \ref pair_mu() have not
70  been implemented and just call the error handler.
71 
72  \note The function \ref calc_density() calls the error handler
73  at zero density and finite temperature, because the correct
74  answer implies \f$ \mu = - \infty \f$ . At zero density and zero
75  temperature the function \ref calc_density() calls \ref
76  calc_density_zerot() which gives the proper chemical potential
77  of \f$ mu = m \f$ without calling the error handler.
78 
79  \future Implement \ref o2scl::fermion_nonrel_tl::pair_density() and
80  \ref o2scl::fermion_nonrel_tl::pair_mu(). [AWS, 1/23/19: it is not
81  entirely clear to me that antiparticles will be useful.]
82 
83  \future This could be improved by performing a Chebyshev
84  approximation (for example) to invert the density integral so
85  that we don't need to use a solver.
86  */
87  template<class inte_t=class fermion_nr_integ_gsl, class fp_t=double>
88  class fermion_nonrel_tl : public fermion_thermo_tl<fp_t> {
89 
90  public:
91 
92  /// Create a nonrelativistic fermion with mass 'm' and degeneracy 'g'
95  }
96 
97  virtual ~fermion_nonrel_tl() {
98  }
99 
100  /// Object for Fermi-Dirac integrals
101  inte_t integ;
102 
103  /** \brief Zero temperature fermions
104  */
105  virtual void calc_mu_zerot(fermion &f) {
106  if (f.non_interacting) { f.nu=f.mu; f.ms=f.m; }
107  if (f.inc_rest_mass) {
108  f.kf=sqrt(2.0*f.ms*(f.nu-f.m));
109  } else {
110  f.kf=sqrt(2.0*f.ms*f.nu);
111  }
112  f.n=f.kf*f.kf*f.kf*f.g/6.0/o2scl_const::pi2;
113  f.ed=f.g*pow(f.kf,5.0)/20.0/o2scl_const::pi2/f.ms;
114  if (f.inc_rest_mass) f.ed+=f.n*f.m;
115  f.pr=-f.ed+f.n*f.nu;
116  f.en=0.0;
117  return;
118  }
119 
120 
121  /** \brief Zero temperature fermions
122  */
123  virtual void calc_density_zerot(fermion &f) {
124  if (f.non_interacting) { f.ms=f.m; }
125  this->kf_from_density(f);
126  f.nu=f.kf*f.kf/2.0/f.ms;
127  f.ed=f.g*pow(f.kf,5.0)/20.0/o2scl_const::pi2/f.ms;
128  if (f.inc_rest_mass) {
129  f.ed+=f.n*f.m;
130  f.nu+=f.m;
131  }
132  f.pr=-f.ed+f.n*f.nu;
133  f.en=0.0;
134 
135  if (f.non_interacting) { f.mu=f.nu; }
136  return;
137  }
138 
139  /** \brief Calculate properties as function of chemical potential
140  */
141  virtual void calc_mu(fermion &f, fp_t temper) {
142 
143  if (temper<0.0) {
144  O2SCL_ERR("Temperature less than zero in fermion_nonrel::calc_mu().",
145  exc_einval);
146  }
147  if (temper==0.0) {
148  calc_mu_zerot(f);
149  return;
150  }
151 
152  fp_t y, sy, spi, ey, int1, int2;
153 
154  if (f.non_interacting) { f.nu=f.mu; f.ms=f.m; }
155 
156  if (f.ms<0.0) {
157  O2SCL_ERR2("Mass negative in ",
158  "fermion_nonrel::calc_mu().",exc_einval);
159  }
160 
161  if (temper<=0.0) {
162  calc_mu_zerot(f);
163  return;
164  }
165 
166  if (f.inc_rest_mass) {
167  y=(f.nu-f.m)/temper;
168  } else {
169  y=f.nu/temper;
170  }
171 
172  // Number density
173 
174  //f.n=gsl_sf_fermi_dirac_half(y)*sqrt(o2scl_const::pi)/2.0;
175  f.n=integ.calc_1o2(y);
176  f.n*=f.g*pow(2.0*f.ms*temper,1.5)/4.0/o2scl_const::pi2;
177 
178  // Energy density:
179 
180  //f.ed=gsl_sf_fermi_dirac_3half(y)*0.75*sqrt(o2scl_const::pi);
181  f.ed=integ.calc_3o2(y);
182  f.ed*=f.g*pow(2.0*f.ms*temper,2.5)/8.0/o2scl_const::pi2/f.ms;
183 
184  if (f.inc_rest_mass) {
185 
186  // Finish energy density
187  f.ed+=f.n*f.m;
188 
189  // entropy density
190  f.en=(5.0*(f.ed-f.n*f.m)/3.0-(f.nu-f.m)*f.n)/temper;
191 
192  // pressure
193  f.pr=2.0*(f.ed-f.n*f.m)/3.0;
194 
195  } else {
196 
197  // entropy density
198  f.en=(5.0*f.ed/3.0-f.nu*f.n)/temper;
199 
200  // pressure
201  f.pr=2.0*f.ed/3.0;
202 
203  }
204 
205  if (!std::isfinite(f.nu) || !std::isfinite(f.n)) {
206  O2SCL_ERR2("Chemical potential or density in ",
207  "fermion_nonrel::calc_mu().",exc_efailed);
208  }
209 
210  return;
211  }
212 
213 
214  /** \brief Calculate properties as function of density
215 
216  If the density is zero, this function just sets part::mu,
217  part::nu, part::ed, part::pr, and part::en to zero and returns
218  without calling the error handler (even though at
219  zero density and finite temperature, the chemical potentials
220  formally are equal to \f$ -\infty \f$).
221  */
222  virtual int calc_density(fermion &f, fp_t temper) {
223 
224  if (f.m<0.0 || (f.non_interacting==false && f.ms<0.0)) {
225  O2SCL_ERR2("Mass negative in ",
226  "fermion_nonrel::calc_density().",exc_einval);
227  }
228  if (temper<0.0) {
229  O2SCL_ERR2("Temperature less than zero in ",
230  "fermion_nonrel::calc_density().",exc_einval);
231  }
232  if (temper==0.0) {
234  return 0;
235  }
236 
237  // Note it is important to throw if n=0 because the correct chemical
238  // potential in that case is mu=-infty and we don't want encourage
239  // the user to use this code in that case.
240  if (f.n<=0.0) {
241  std::string str=((std::string)"Density, ")+o2scl::dtos(f.n)+
242  ", less than zero when temperature is positive in "+
243  "fermion_nonrel::calc_density().";
244  O2SCL_ERR(str.c_str(),exc_einval);
245  }
246 
247  if (f.non_interacting==true) { f.nu=f.mu; f.ms=f.m; }
248 
249  nu_from_n(f,temper);
250 
251  if (f.non_interacting) { f.mu=f.nu; }
252 
253  fp_t y, spi, ey, sy;
254  if (f.inc_rest_mass) {
255  y=-(f.nu-f.m)/temper;
256  } else {
257  y=-f.nu/temper;
258  }
259 
260  // energy density
261  //f.ed=gsl_sf_fermi_dirac_3half(-y)*sqrt(o2scl_const::pi)*0.75;
262  f.ed=integ.calc_3o2(-y);
263 
264  if (f.inc_rest_mass) {
265 
266  // Finish energy density
267  f.ed*=f.g*pow(2.0*f.ms*temper,2.5)/8.0/o2scl_const::pi2/f.ms;
268  f.ed+=f.n*f.m;
269 
270  // entropy density
271  f.en=(5.0*(f.ed-f.n*f.m)/3.0-(f.nu-f.m)*f.n)/temper;
272 
273  // pressure
274  f.pr=2.0*(f.ed-f.n*f.m)/3.0;
275 
276  } else {
277 
278  // Finish energy density
279  f.ed*=f.g*pow(2.0*f.ms*temper,2.5)/8.0/o2scl_const::pi2/f.ms;
280 
281  // entropy density
282  f.en=(5.0*f.ed/3.0-f.nu*f.n)/temper;
283 
284  // pressure
285  f.pr=2.0*f.ed/3.0;
286 
287  }
288 
289  return 0;
290  }
291 
292  /** \brief Compute thermodynamics with antiparticles at fixed
293  chemical potential (unimplemented)
294  */
295  virtual void pair_mu(fermion &f, fp_t temper) {
296  O2SCL_ERR2("Function fermion_nonrel::pair_mu() not ",
297  "implemented.",exc_eunimpl);
298  return;
299  }
300 
301  /** \brief Compute thermodynamics with antiparticles at fixed
302  density (unimplemented)
303  */
304  virtual int pair_density(fermion &f, fp_t temper) {
305  O2SCL_ERR2("Function fermion_nonrel::pair_density() not ",
306  "implemented.",exc_eunimpl);
307  return 0;
308  }
309 
310  /// Calculate effective chemical potential from density
311  virtual void nu_from_n(fermion &f, fp_t temper) {
312 
313  fp_t init_n=f.n, init_m=f.m, init_ms=f.ms, init_nu=f.nu;
314 
315  // Use initial value of nu for initial guess
316  fp_t nex;
317  if (f.inc_rest_mass) {
318  nex=-(f.nu-f.m)/temper;
319  } else {
320  nex=-f.nu/temper;
321  }
322 
323  // Make a correction if nex is too small and negative
324  // (Note GSL_LOG_DBL_MIN is about -708)
325  if (nex>-GSL_LOG_DBL_MIN*0.9) nex=-GSL_LOG_DBL_MIN/2.0;
326 
327  funct mf=std::bind(std::mem_fn<fp_t(fp_t,fp_t,fp_t)>
329  this,std::placeholders::_1,f.n/f.g,f.ms*temper);
330 
331  // Turn off convergence errors temporarily, since we'll
332  // try again if it fails
333  bool enc=density_root->err_nonconv;
334  density_root->err_nonconv=false;
335  int ret=density_root->solve(nex,mf);
336 
337  // The root_cern class has a hard time when nex is near zero,
338  // so we try a bracketing solver
339  if (ret!=0) {
340  fp_t blow=fabs(nex), bhigh=-blow;
341  fp_t yhigh=mf(bhigh), ylow=mf(blow);
342  for(size_t j=0;j<10 && yhigh*ylow>0.0;j++) {
343  fp_t delta=fabs(blow);
344  blow+=delta;
345  bhigh-=delta;
346  yhigh=mf(bhigh);
347  ylow=mf(blow);
348  }
349  if ((yhigh<0.0 && ylow>0.0) || (yhigh>0.0 && ylow<0.0)) {
351  rbg.err_nonconv=false;
352  ret=rbg.solve_bkt(blow,bhigh,mf);
353  if (ret==0) nex=blow;
354  }
355  }
356 
357  if (ret!=0) {
358 
359  // If it failed, try to get a guess from classical_thermo particle
360 
361  classical_thermo cl;
362  cl.calc_density(f,temper);
363  if (f.inc_rest_mass) {
364  nex=-(f.nu-f.m)/temper;
365  } else {
366  nex=-f.nu/temper;
367  }
368  ret=density_root->solve(nex,mf);
369 
370  // If it failed again, add error information
371  if (ret!=0) {
372  std::cout << "Function fermion_nonrel::nu_from_n() failed."
373  << std::endl;
374  std::cout.precision(14);
375  std::cout << " n,m,ms,T,nu: " << init_n << " " << init_m << " "
376  << init_ms << " " << temper << " " << init_nu << std::endl;
377  std::cout << " ni,irm: " << f.non_interacting << " "
378  << f.inc_rest_mass << std::endl;
379  O2SCL_ERR("Solver failed in fermion_nonrel::nu_from_n().",ret);
380  }
381  }
382 
383  // Restore the value of density_root->err_nonconv
385 
386  if (f.inc_rest_mass) {
387  f.nu=-nex*temper+f.m;
388  } else {
389  f.nu=-nex*temper;
390  }
391 
392  return;
393  }
394 
395 
396  /** \brief Set the solver for use in calculating the chemical
397  potential from the density
398  */
400  density_root=&rp;
401  return;
402  }
403 
404  /// The default solver for calc_density().
406 
407  /// Return string denoting type ("fermion_nonrel")
408  virtual const char *type() { return "fermion_nonrel"; }
409 
410  /** \brief Function to compute chemical potential from density
411 
412  Variable \c nog is the target baryon density divided by
413  the spin degeneracy, and \c msT is the effective mass
414  times the temperature.
415  */
416  fp_t solve_fun(fp_t x, fp_t nog, fp_t msT) {
417 
418  fp_t nden;
419 
420  // If the argument to gsl_sf_fermi_dirac_half() is less
421  // than GSL_LOG_DBL_MIN (which is about -708), then
422  // an underflow occurs. We just set nden to zero in this
423  // case, as this helps the solver find the right root.
424 
425  if (((-x)<GSL_LOG_DBL_MIN) || !std::isfinite(x)) {
426  nden=0.0;
427  } else {
428  //nden=gsl_sf_fermi_dirac_half(-x)*sqrt(o2scl_const::pi)/2.0;
429  nden=integ.calc_1o2(-x);
430 
431  }
432 
433  nden*=pow(2.0*msT,1.5)/4.0/o2scl_const::pi2;
434  fp_t ret=nden/nog-1.0;
435  return ret;
436  }
437 
438  protected:
439 
440 #ifndef DOXYGEN_NO_O2NS
441 
442  /// Solver to compute chemical potential from density
444 
445  private:
446 
448  fermion_nonrel_tl& operator=(const fermion_nonrel_tl&);
449 
450 #endif
451 
452  };
453 
455 
456 #ifndef DOXYGEN_NO_O2NS
457 }
458 #endif
459 
460 #endif
o2scl::classical_thermo_tl::calc_density
virtual void calc_density(part_tl< fp_t > &p, fp_t temper)
Calculate properties as function of density.
Definition: classical.h:130
o2scl::fermion_nonrel_tl::pair_density
virtual int pair_density(fermion &f, fp_t temper)
Compute thermodynamics with antiparticles at fixed density (unimplemented)
Definition: fermion_nonrel.h:304
o2scl::root::solve
virtual int solve(fp_t &x, func_t &func)=0
o2scl::fermion_nonrel_tl::calc_density_zerot
virtual void calc_density_zerot(fermion &f)
Zero temperature fermions.
Definition: fermion_nonrel.h:123
o2scl_const::pi2
const double pi2
o2scl::part_tl::mu
fp_t mu
Chemical potential.
Definition: part.h:118
o2scl::part_tl::ms
fp_t ms
Effective mass (Dirac unless otherwise specified)
Definition: part.h:122
exc_efailed
exc_efailed
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
o2scl::part_tl::n
fp_t n
Number density.
Definition: part.h:112
o2scl::fermion_nonrel_tl::fermion_nonrel_tl
fermion_nonrel_tl()
Create a nonrelativistic fermion with mass 'm' and degeneracy 'g'.
Definition: fermion_nonrel.h:93
o2scl::fermion_nonrel_tl::solve_fun
fp_t solve_fun(fp_t x, fp_t nog, fp_t msT)
Function to compute chemical potential from density.
Definition: fermion_nonrel.h:416
o2scl::fermion_nonrel_tl::pair_mu
virtual void pair_mu(fermion &f, fp_t temper)
Compute thermodynamics with antiparticles at fixed chemical potential (unimplemented)
Definition: fermion_nonrel.h:295
o2scl::fermion_nonrel_tl::calc_mu
virtual void calc_mu(fermion &f, fp_t temper)
Calculate properties as function of chemical potential.
Definition: fermion_nonrel.h:141
o2scl::root_brent_gsl::solve_bkt
virtual int solve_bkt(fp_t &x1, fp_t x2, func_t &f)
o2scl::part_tl::non_interacting
bool non_interacting
True if the particle is non-interacting (default true)
Definition: part.h:130
exc_eunimpl
exc_eunimpl
o2scl::part_tl::en
fp_t en
Entropy density.
Definition: part.h:120
o2scl::root
o2scl::fermion_nonrel_tl::integ
inte_t integ
Object for Fermi-Dirac integrals.
Definition: fermion_nonrel.h:101
o2scl::part_tl::g
fp_t g
Degeneracy (e.g. spin and color if applicable)
Definition: part.h:108
o2scl::fermion_nonrel_tl::calc_density
virtual int calc_density(fermion &f, fp_t temper)
Calculate properties as function of density.
Definition: fermion_nonrel.h:222
o2scl::part_tl::ed
fp_t ed
Energy density.
Definition: part.h:114
o2scl::fermion_nonrel_tl::set_density_root
void set_density_root(root<> &rp)
Set the solver for use in calculating the chemical potential from the density.
Definition: fermion_nonrel.h:399
o2scl::root_brent_gsl
o2scl::fermion_tl< double >
exc_einval
exc_einval
o2scl::fermion_zerot_tl< double >::kf_from_density
void kf_from_density(fermion_tl< double > &f)
Calculate the Fermi momentum from the density.
Definition: fermion.h:157
o2scl::part_tl::inc_rest_mass
bool inc_rest_mass
If true, include the mass in the energy density and chemical potential (default true)
Definition: part.h:128
o2scl::part_tl::m
fp_t m
Mass.
Definition: part.h:110
o2scl::fermion_nonrel_tl::def_density_root
root_cern def_density_root
The default solver for calc_density().
Definition: fermion_nonrel.h:405
o2scl::part_tl::pr
fp_t pr
Pressure.
Definition: part.h:116
o2scl::fermion_tl::kf
fp_t kf
Fermi momentum.
Definition: fermion.h:60
o2scl::root_cern
o2scl::fermion_nonrel_tl::type
virtual const char * type()
Return string denoting type ("fermion_nonrel")
Definition: fermion_nonrel.h:408
o2scl::classical_thermo_tl< double >
O2SCL_ERR
#define O2SCL_ERR(d, n)
o2scl::funct
std::function< double(double)> funct
o2scl::fermion_nonrel_tl::nu_from_n
virtual void nu_from_n(fermion &f, fp_t temper)
Calculate effective chemical potential from density.
Definition: fermion_nonrel.h:311
o2scl::fermion_nonrel_tl::density_root
root * density_root
Solver to compute chemical potential from density.
Definition: fermion_nonrel.h:443
o2scl::fermion_thermo_tl
Fermion with finite-temperature thermodynamics [abstract base].
Definition: fermion.h:309
o2scl::fermion_nonrel_tl
Nonrelativistic fermion class.
Definition: fermion_nonrel.h:88
o2scl::dtos
std::string dtos(double x, std::ostream &format)
o2scl::fermion_nonrel_tl::calc_mu_zerot
virtual void calc_mu_zerot(fermion &f)
Zero temperature fermions.
Definition: fermion_nonrel.h:105
o2scl::root::err_nonconv
bool err_nonconv
o2scl::part_tl::nu
fp_t nu
Effective chemical potential.
Definition: part.h:124

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).