ConfinementPotential.C 2.45 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17 18 19 20 21 22
// ConfinementPotential.C
//
////////////////////////////////////////////////////////////////////////////////

#include "ConfinementPotential.h"
#include "Basis.h"

////////////////////////////////////////////////////////////////////////////////
23
ConfinementPotential::ConfinementPotential(double ecuts, double facs,
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
  double sigmas, const Basis& basis): ecuts_(ecuts), facs_(facs),
  sigmas_(sigmas), basis_(basis)
{
  const int ngwloc = basis_.localsize();
  fstress_.resize(ngwloc);
  dfstress_.resize(ngwloc);
  update();
}

////////////////////////////////////////////////////////////////////////////////
void ConfinementPotential::update(void)
{
  // recompute confinement potential and its derivative
  const double sigmas_inv = 1.0 / sigmas_;
  const int ngwloc = basis_.localsize();
39

40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
  for ( int ig = 0; ig < ngwloc; ig++ )
  {
    const double gsq = basis_.kpg2(ig);
    // Next line: 0.5 from 1/2m
    const double arg = ( 0.5 * gsq - ecuts_ ) * sigmas_inv;
    // Next lines: fp = fermi(arg);
    // fm = fermi(-arg) = 1 - fp;
    double fm,fp;
    if ( arg > 50.0 )
    {
      fm = 1.0;
      fp = 0.0;
    }
    else if ( arg < -50.0 )
    {
      fm = 0.0;
      fp = 1.0;
    }
    else
    {
      fp = 1.0 / ( 1.0 + exp ( arg ) );
      fm = 1.0 - fp;
    }

    // f(G) = facs * ( 1 - fermi( (G^2-ecuts)/sigmas ) )
    // fg = f(G)
    const double fg = facs_ * fm;
    // gfgp = G f'(G)
    const double gfgp = gsq * fg * fp * sigmas_inv;
69

70 71 72 73
    if ( ecuts_ > 0.0 )
    {
      // fstress[ig] = G^2 * f(G)
      fstress_[ig] = gsq * fg;
74

75 76 77 78 79 80 81 82
      // dfstress =  2 f(G) + G * f'(G)
      dfstress_[ig] = 2.0 * fg + gfgp;
    }
    else
    {
      fstress_[ig] = 0.0;
      dfstress_[ig] = 0.0;
    }
83 84 85 86 87 88 89

    // ekin = sum_G |c_G|^2  G^2
    // econf = sum_G |c_G|^2 fstress[G]
    // stress_ekin_ij = (1/Omega) sum_G |c_G|^2 * 2 * G_i * G_j
    // stress_econf_ij = (1/Omega) sum_G |c_G|^2 * dfstress[G] * G_i * G_j
  }
}