---

more/phys/bcs.h

Go to the documentation of this file.

//  Copyright (C) 2000--2002  Petter Urkedal (petter.urkedal@matfys.lth.se)
//
//  This file 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 2 of the License, or
//  (at your option) any later version.
//
//  This file 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 this program; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//  As a special exception, you may use this file as part of a free
//  software library without restriction.  Specifically, if other files
//  instantiate templates or use macros or inline functions from this
//  file, or you compile this file and link it with other files to
//  produce an executable, this file does not by itself cause the
//  resulting executable to be covered by the GNU General Public
//  License.  This exception does not however invalidate any other
//  reasons why the executable file might be covered by the GNU General
//  Public License.
//
//  $Id: bcs.h,v 1.1 2002/05/30 18:01:40 petter_urkedal Exp $

/** \file more/phys/bcs.h
 ** Provides a BCS solver. */

#ifndef PHYSICS_BCS_H
#define PHYSICS_BCS_H

#include <algorithm>
#include <vector>
#include <more/gen/iterator.h>
#include <more/math/math.h>


namespace more {
namespace phys {

  /** \class BCS_solver bcs.h more/phys/bcs.h
   ** A solver of the Bardeen–Cooper–Schrieffer (BCS) equations.
   ** Energy */
  template<typename T>
    struct BCS_solver
    {
        typedef T real_type;
      private:
        typedef std::vector<real_type> energy_container;
        typedef typename energy_container::iterator energy_iterator;
      public:
        typedef typename energy_container::size_type size_type;

        BCS_solver(int N, real_type G,
                   real_type cut_low, real_type cut_high,
                   real_type cut_width)
            : m_is_solved(false),
              m_n_occ(N), m_G(G),
              m_do_smooth_cut(true),
              m_cut_low(cut_low), m_cut_high(cut_high),
              m_cut_width(cut_width),
              m_prec(s_prec_default),
              m_n_iter_max(200),
              m_is_verbose(false) {}

        BCS_solver(int N, real_type G)
            : m_is_solved(false),
              m_n_occ(N), m_G(G),
              m_do_smooth_cut(false),
              m_prec(s_prec_default),
              m_n_iter_max(200),
              m_is_verbose(false) {}

        BCS_solver()
            : m_is_solved(false),
              m_n_occ(0), m_G(0),
              m_do_smooth_cut(false),
              m_prec(s_prec_default),
              m_n_iter_max(200),
              m_is_verbose(false) {}

        /** Set the coupling strength. */
        void set_coupling_strength(real_type G) { m_G = G; }

        /** Return the coupling strength. */
        real_type coupling_strength() const { return m_G; }

        /** Set the number of occupied states = the number of
         ** particles.  */
        void set_occupied_count(size_type n) { m_n_occ = n; }

        /** Return the number of occupied states. */
        real_type occupied_count() const { return m_n_occ; }

        /** Set the wanted precision for Δ and λ. */
        void set_precision(real_type prec) { m_prec = prec; }

        /** Return the precision. */
        real_type precision() const { return m_prec; }

        /** Enable smooth cut to window with energies from \a cut_low
         ** to \a cut_high, bondaries being smeared over a range \a
         ** cut_width. */
        void set_smooth_cutoff(real_type cut_low, real_type cut_high,
                               real_type cut_width);

        /** If doing smooth cut, return true and return the details in
         ** the arguments, else return false. */
        bool smooth_cut(real_type& cut_low, real_type& cut_high,
                        real_type& cut_width);

        void solve();

        /** Take the energy levels and solve the BCS equations. */
        template<typename Iterator>
          void solve(Iterator first, Iterator last)
          {
              m_energies.clear();
              m_energies.resize(std::distance(first, last));
              std::copy(first, last, m_energies.begin());
              solve();
          }

        /** Return $v$ for the given energy.
         ** \pre The equations are solved. */
        real_type sqr_v(real_type epsilon)
        {
            req_solved();
            real_type eml = epsilon - m_lambda;
            return .5 - .5*eml/sqrt(eml*eml + m_Delta2);
        }

        /** Return $u$ for the given energy.
         ** \pre The equations are solved. */
        real_type sqr_u(real_type epsilon)
        {
            req_solved();
            real_type eml = epsilon - m_lambda;
            return .5 + .5*eml/sqrt(eml*eml + m_Delta2);
        }

        void write_to(std::ostream& out) const;

        /** Return the pairing energy.
         ** \pre The equations are solved. */
        real_type pairing_energy() { req_solved(); return m_E_pair; }

        /** Return the pairing gap Δ.
         ** \pre The equations are solved. */
        real_type pairing_Delta() { req_solved(); return m_Delta; }

        /** Return the Fermi energy λ.
         ** \pre The equations are solved. */
        real_type Fermi_energy() { req_solved(); return m_lambda; }

        real_type cut_function(real_type epsilon)
        {
            return internal_cut_function(epsilon-m_lambda);
        }

        /** Show each step of the iteration iff \a value is \c true. */
        void be_verbose(bool value = true) { m_is_verbose = value; }

        /** Return the value of the verbosity option. */
        bool is_verbose() const { return m_is_verbose; }

      private:
        bool m_is_solved;
        size_type m_n_occ;
        real_type m_G;

        bool m_do_smooth_cut;
        real_type m_cut_low, m_cut_high, m_cut_width;

        real_type m_E_pair, m_E_pair_nocorr;
        real_type m_Delta, m_Delta2;
        real_type m_lambda;

        real_type m_prec;
        int m_n_iter_max;
        std::vector<real_type> m_energies;

        bool m_is_verbose;

        static real_type s_prec_default;

        real_type
        internal_cut_function(real_type eml) const
        {
            if(m_do_smooth_cut)
                return 1/(1 + exp((eml-m_cut_high)/m_cut_width)
                            + exp(-(eml-m_cut_low)/m_cut_width));
            else
                return 1.0;
        }

        void
        req_solved() const
        {
            if (!m_is_solved)
                throw std::logic_error("more::phys::BCS_solver: "
                                       "The equations are not solved yet.");
        }
    };

  template<typename T>
    inline std::ostream&
    operator<<(std::ostream& os, const BCS_solver<T>& bcs)
    {
        bcs.write_to(os); return os;
    }

}} // more::phys

#endif

---