XCOperator.C 4.85 KB
Newer Older
Francois Gygi committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed 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.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// XCOperator.C
//
////////////////////////////////////////////////////////////////////////////////
#include "XCOperator.h"
19
#include "ChargeDensity.h"
Francois Gygi committed
20 21
#include "XCPotential.h"
#include "ExchangeOperator.h"
22
#include "HSEFunctional.h"
Francois Gygi committed
23
#include "RSHFunctional.h"
Francois Gygi committed
24 25 26 27 28 29 30 31 32
using namespace std;

////////////////////////////////////////////////////////////////////////////////
XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd)
{
  // set initial values
  xcp_ = 0;
  xop_ = 0;
  exc_ = 0.0 ;
33
  dxc_ = 0.0 ;
34 35

  sigma_exc_.resize(6);
Francois Gygi committed
36 37 38 39 40

  string functional_name = s.ctrl.xc;

  // check the name of the functional
  if ( ( functional_name ==  "LDA" ) ||
Francois Gygi committed
41
       ( functional_name ==  "VWN" ) ||
Francois Gygi committed
42 43 44 45
       ( functional_name ==  "PBE" ) ||
       ( functional_name == "BLYP" ) )
  {
    // create only an xc potential
46
    xcp_ = new XCPotential(cd, functional_name, s.ctrl);
Francois Gygi committed
47 48 49 50 51 52 53
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = false;
  }
  else if ( functional_name == "HF" )
  {
    // create exchange operator with mixing coeff=1
Francois Gygi committed
54
    xop_ = new ExchangeOperator(s, 1.0, 1.0, 0.0);
Francois Gygi committed
55 56 57 58 59 60 61
    hasPotential_ = false;
    hasGGA_ = false;
    hasHF_ = true;
  }
  else if ( functional_name == "PBE0" )
  {
    // create an exchange potential
62
    xcp_ = new XCPotential(cd, functional_name, s.ctrl);
63 64

    // create the exchange operator with mixing coeff=0.25
Francois Gygi committed
65
    xop_ = new ExchangeOperator(s, s.ctrl.alpha_PBE0, s.ctrl.alpha_PBE0, 0.0);
66 67 68 69 70 71 72
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = true;
  }
  else if ( functional_name == "HSE" )
  {
    // create an exchange potential
73
    xcp_ = new XCPotential(cd, functional_name, s.ctrl);
Francois Gygi committed
74 75

    // create the exchange operator with mixing coeff=0.25
Francois Gygi committed
76 77 78 79 80 81 82 83 84 85 86 87 88
    xop_ = new ExchangeOperator(s, 0.0, 0.25, 0.11);
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = true;
  }
  else if ( functional_name == "RSH" )
  {
    // create an exchange potential
    xcp_ = new XCPotential(cd, functional_name, s.ctrl);

    // create the exchange operator with mixing coeff=beta_RSH
    xop_ = new ExchangeOperator(s, s.ctrl.alpha_RSH, s.ctrl.beta_RSH,
      s.ctrl.mu_RSH);
Francois Gygi committed
89 90 91 92 93 94 95
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = true;
  }
  else if ( functional_name == "B3LYP" )
  {
    // create an exchange potential
96
    xcp_ = new XCPotential(cd, functional_name, s.ctrl);
Francois Gygi committed
97 98

    // create the exchange operator with mixing coeff=0.20
Francois Gygi committed
99
    xop_ = new ExchangeOperator(s, 0.20, 0.20, 0.0);
Francois Gygi committed
100 101 102 103
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = true;
  }
104 105 106 107 108 109 110 111 112 113 114
  else if ( functional_name == "BHandHLYP" )
  {
    // create an exchange potential
    xcp_ = new XCPotential(cd, functional_name, s.ctrl);

    // create the exchange operator with mixing coeff=0.50
    xop_ = new ExchangeOperator(s, 0.50, 0.50, 0.0);
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = true;
  }
Francois Gygi committed
115 116 117 118 119 120 121 122 123 124 125 126 127 128
  else
  {
    throw XCOperatorException(
      "unknown functional name during exchange operator construction");
  }
}

XCOperator::~XCOperator()
{
  delete xcp_;
  delete xop_;
}

////////////////////////////////////////////////////////////////////////////////
129
void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stress)
Francois Gygi committed
130
{
131 132
  // update xc potential and self-energy
  // used whenever the charge density and/or wave functions have changed
Francois Gygi committed
133 134 135 136 137 138 139 140
  // compute vxc potential and energy
  if ( hasPotential_ )
  {
    // update LDA/GGA xc potential
    xcp_->update( vr );

    // LDA/GGA exchange energy
    exc_ = xcp_->exc();
141
    dxc_ = xcp_->dxc();
142 143 144

    if ( compute_stress )
      xcp_->compute_stress(sigma_exc_);
Francois Gygi committed
145 146 147 148
  }
  else
  {
    exc_ = 0.0;
149
    dxc_ = 0.0;
150
    sigma_exc_ = 0.0;
Francois Gygi committed
151 152 153 154
  }

  if ( hasHF() )
  {
155 156 157
    double ex_hf = xop_->update_operator(compute_stress);
    exc_ += ex_hf;
    dxc_ -= ex_hf;
158 159
    if ( compute_stress )
      xop_->add_stress(sigma_exc_);
Francois Gygi committed
160 161 162 163
  }
}

////////////////////////////////////////////////////////////////////////////////
164
void XCOperator::apply_self_energy(Wavefunction &dwf)
Francois Gygi committed
165 166
{
  if ( hasHF() )
167
    xop_->apply_operator(dwf);
Francois Gygi committed
168 169 170
}

////////////////////////////////////////////////////////////////////////////////
171
void XCOperator::compute_stress(std::valarray<double>& sigma)
Francois Gygi committed
172
{
173
  sigma = sigma_exc_;
Francois Gygi committed
174
}
175 176 177 178 179 180 181

////////////////////////////////////////////////////////////////////////////////
void XCOperator::cell_moved(void)
{
  if ( hasHF() )
    xop_->cell_moved();
}