XCOperator.C 5.16 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"
Francois Gygi committed
19
#include "Sample.h"
20
#include "ChargeDensity.h"
Francois Gygi committed
21 22
#include "XCPotential.h"
#include "ExchangeOperator.h"
23
#include "HSEFunctional.h"
Francois Gygi committed
24
#include "RSHFunctional.h"
Francois Gygi committed
25 26 27
using namespace std;

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
28
XCOperator::XCOperator(Sample& s, const ChargeDensity& cd) :cd_(cd), s_(s)
Francois Gygi committed
29 30 31 32 33
{
  // set initial values
  xcp_ = 0;
  xop_ = 0;
  exc_ = 0.0 ;
34
  dxc_ = 0.0 ;
Francois Gygi committed
35 36
  hasHF_ = false;
  hasMeta_ = false;
37 38

  sigma_exc_.resize(6);
Francois Gygi committed
39 40 41 42 43

  string functional_name = s.ctrl.xc;

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

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

    // create the exchange operator with mixing coeff=0.25
Francois Gygi committed
79 80 81 82 83 84 85 86
    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
Francois Gygi committed
87
    xcp_ = new XCPotential(cd, functional_name, s_);
Francois Gygi committed
88 89 90 91

    // 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
92 93 94 95 96 97 98
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = true;
  }
  else if ( functional_name == "B3LYP" )
  {
    // create an exchange potential
Francois Gygi committed
99
    xcp_ = new XCPotential(cd, functional_name, s_);
Francois Gygi committed
100 101

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

    // 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
118 119 120 121 122 123 124 125
  else if ( functional_name == "SCAN" )
  {
    // create an exchange potential
    xcp_ = new XCPotential(cd, functional_name, s_);
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasMeta_ = true;
  }
Francois Gygi committed
126 127 128 129 130 131 132 133 134 135 136 137 138 139
  else
  {
    throw XCOperatorException(
      "unknown functional name during exchange operator construction");
  }
}

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

////////////////////////////////////////////////////////////////////////////////
140
void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stress)
Francois Gygi committed
141
{
142 143
  // update xc potential and self-energy
  // used whenever the charge density and/or wave functions have changed
Francois Gygi committed
144 145 146
  // compute vxc potential and energy
  if ( hasPotential_ )
  {
Francois Gygi committed
147
    // update LDA/GGA/MetaGGA xc potential
Francois Gygi committed
148 149 150 151
    xcp_->update( vr );

    // LDA/GGA exchange energy
    exc_ = xcp_->exc();
152
    dxc_ = xcp_->dxc();
153 154 155

    if ( compute_stress )
      xcp_->compute_stress(sigma_exc_);
Francois Gygi committed
156 157 158 159
  }
  else
  {
    exc_ = 0.0;
160
    dxc_ = 0.0;
161
    sigma_exc_ = 0.0;
Francois Gygi committed
162 163 164 165
  }

  if ( hasHF() )
  {
166 167 168
    double ex_hf = xop_->update_operator(compute_stress);
    exc_ += ex_hf;
    dxc_ -= ex_hf;
169 170
    if ( compute_stress )
      xop_->add_stress(sigma_exc_);
Francois Gygi committed
171 172 173 174
  }
}

////////////////////////////////////////////////////////////////////////////////
175
void XCOperator::apply_self_energy(Wavefunction &dwf)
Francois Gygi committed
176 177
{
  if ( hasHF() )
178
    xop_->apply_operator(dwf);
Francois Gygi committed
179 180
  if ( hasMeta() )
    xcp_->apply_meta_operator(dwf);
Francois Gygi committed
181 182 183
}

////////////////////////////////////////////////////////////////////////////////
184
void XCOperator::compute_stress(std::valarray<double>& sigma)
Francois Gygi committed
185
{
186
  sigma = sigma_exc_;
Francois Gygi committed
187
}
188 189 190 191 192 193 194

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