XCOperator.C 3.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 22 23 24 25 26 27 28 29 30
#include "XCPotential.h"
#include "ExchangeOperator.h"
using namespace std;

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

  sigma_exc_.resize(6);
Francois Gygi committed
34 35 36 37 38

  string functional_name = s.ctrl.xc;

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

    // create the exchange operator with mixing coeff=0.25
65
    xop_ = new ExchangeOperator(s, s.ctrl.alpha_PBE0);
Francois Gygi committed
66 67 68
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = true;
69
    HFmixCoeff_ = s.ctrl.alpha_PBE0;;
Francois Gygi committed
70 71 72 73
  }
  else if ( functional_name == "B3LYP" )
  {
    // create an exchange potential
74
    xcp_ = new XCPotential(cd, functional_name, s.ctrl);
Francois Gygi committed
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

    // create the exchange operator with mixing coeff=0.20
    xop_ = new ExchangeOperator(s, 0.20);
    hasPotential_ = true;
    hasGGA_ = xcp_->isGGA();
    hasHF_ = true;
    HFmixCoeff_ = 0.20;
  }
  else
  {
    throw XCOperatorException(
      "unknown functional name during exchange operator construction");
  }
}

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

////////////////////////////////////////////////////////////////////////////////
97
void XCOperator::update(std::vector<std::vector<double> >& vr, bool compute_stress)
Francois Gygi committed
98
{
99 100
  // update xc potential and self-energy
  // used whenever the charge density and/or wave functions have changed
Francois Gygi committed
101 102 103 104 105 106 107 108
  // compute vxc potential and energy
  if ( hasPotential_ )
  {
    // update LDA/GGA xc potential
    xcp_->update( vr );

    // LDA/GGA exchange energy
    exc_ = xcp_->exc();
109
    dxc_ = xcp_->dxc();
110 111 112

    if ( compute_stress )
      xcp_->compute_stress(sigma_exc_);
Francois Gygi committed
113 114 115 116
  }
  else
  {
    exc_ = 0.0;
117
    dxc_ = 0.0;
118
    sigma_exc_ = 0.0;
Francois Gygi committed
119 120 121 122
  }

  if ( hasHF() )
  {
123 124 125
    double ex_hf = xop_->update_operator(compute_stress);
    exc_ += ex_hf;
    dxc_ -= ex_hf;
126 127
    if ( compute_stress )
      xop_->add_stress(sigma_exc_);
Francois Gygi committed
128 129 130 131
  }
}

////////////////////////////////////////////////////////////////////////////////
132
void XCOperator::apply_self_energy(Wavefunction &dwf)
Francois Gygi committed
133 134
{
  if ( hasHF() )
135
    xop_->apply_operator(dwf);
Francois Gygi committed
136 137 138
}

////////////////////////////////////////////////////////////////////////////////
139
void XCOperator::compute_stress(std::valarray<double>& sigma)
Francois Gygi committed
140
{
141
  sigma = sigma_exc_;
Francois Gygi committed
142
}
143 144 145 146 147 148 149

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