Commit 8dc36961 by He Ma

automatic determination of amplitude

git-svn-id: http://qboxcode.org/svn/qb/branches/vext@1906 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 6d0d2a37
......@@ -19,6 +19,9 @@
#include <iostream>
#include <fstream>
#include <cassert>
#include <algorithm>
#include <numeric>
#include <functional>
using namespace std;
#include "Basis.h"
......@@ -59,16 +62,39 @@ void ExternalPotential::update(const ChargeDensity& cd)
getline(vfile,tmp); // skip atom coordinates
const int n012 = n_[0] * n_[1] * n_[2];
const int n12 = n_[1] * n_[2];
vext_read.resize(n012);
// build a min heap to store the largest values of vext
// to compute its magnitude
const int heapsize = n012/1000 + 1; // +1 to avoid getting 0
vector<double> heap(heapsize, 0.0);
double tmp_value;
int counter = 0;
for ( int nx = 0; nx < n_[0]; nx++ )
for ( int ny = 0; ny < n_[1]; ny++ )
for ( int nz = 0; nz < n_[2]; nz++ )
{
vfile >> tmp_value;
if ( counter < heapsize )
heap[counter] = abs(tmp_value);
else
{
if ( counter == heapsize )
make_heap(heap.begin(), heap.end());
if ( abs(tmp_value) > heap[0] )
{
pop_heap(heap.begin(), heap.end(), greater<double>());
heap[heapsize-1] = abs(tmp_value);
push_heap(heap.begin(), heap.end(), greater<double>());
}
}
const int ir = nx + ny * n_[0] + nz * n_[0] * n_[1];
vfile >> vext_read[ir];
vext_read[ir] = tmp_value;
counter ++;
}
vfile.close();
magnitude_ = accumulate(heap.begin(), heap.end(), 0.0) / heapsize;
}
else
{
......@@ -86,8 +112,9 @@ void ExternalPotential::update(const ChargeDensity& cd)
<< n_[0] << " " << n_[1] << " " << n_[2] << endl;
}
// broadcast sizes to lower rows
// broadcast sizes and magnitude to lower rows
MPI_Bcast(&n_[0],3,MPI_INT,0,vcomm);
MPI_Bcast(&magnitude_,1,MPI_DOUBLE,0,vcomm);
// create a Basis compatible with the vext grid read from file
const Basis* vbasis = cd.vbasis();
......@@ -107,6 +134,7 @@ void ExternalPotential::update(const ChargeDensity& cd)
if ( cd_ctxt.onpe0() )
{
cout << " ExternalPotential::update: magnitude: " << magnitude_ << endl;
cout << " ExternalPotential::update: ecut0: " << ecut0 << endl;
cout << " ExternalPotential::update: ecut1: " << ecut1 << endl;
cout << " ExternalPotential::update: ecut2: " << ecut2 << endl;
......@@ -160,6 +188,18 @@ void ExternalPotential::update(const ChargeDensity& cd)
ft2.backward(&vext_g_[0],&tmp_r[0]);
for ( int i = 0; i < vext_r_.size(); i++ )
vext_r_[i] = real(tmp_r[i]);
if ( amplitude_ == 0.0 )
{
// If amplitude_ = 0.0, use following scheme to get an amplitude.
// Empirically, an absolute magnitude of 1.0E-3 ~ 1.0E-5 hartree for Vext
// would be suitable. Here the amplitude is set to scale the
// magnitude of vext to be 1.0E-3 hartree
amplitude_ = 1.0E-3 / magnitude_;
if ( cd_ctxt.onpe0() )
cout << " ExternalPotential::update: amplitude automatically determined to be "
<< amplitude_ << endl;
}
}
double ExternalPotential::compute_eext(const ChargeDensity& cd)
......
......@@ -32,7 +32,9 @@ class ExternalPotential
Sample& s_;
int n_[3]; // real space grid size in 3 dimensions
double ecut_;
double amplitude_;
double magnitude_; // the magnitude of external potential, defined as
// the average of its largest 0.1% (absolute) values
double amplitude_; // overall scaling factor of external potential
vector<double> vext_r_; // vext in real space
std::string filename_; // file name for external potential
......@@ -44,6 +46,7 @@ class ExternalPotential
int n(int i) const { return n_[i]; }
double ecut(void) const { return ecut_; }
double magnitude(void) const { return magnitude_; }
double amplitude(void) const { return amplitude_; }
std::string filename(void) const { return filename_; }
double v(size_t i) const { return amplitude_ * vext_r_[i]; }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment