mirror of https://github.com/QMCPACK/qmcpack.git
*** empty log message ***
git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@10 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
parent
d1dd9b432d
commit
2d929e6d0c
|
@ -19,6 +19,7 @@ LINK_LIBRARIES(${QWT_LIBRARY} ${QT_LIBRARIES} ${OPENGL_LIBRARIES})
|
||||||
ENDIF(QT_FOUND)
|
ENDIF(QT_FOUND)
|
||||||
|
|
||||||
SET(HFSRCS
|
SET(HFSRCS
|
||||||
|
../OhmmsApp/ProjectData.cpp
|
||||||
../Numerics/Clebsch_Gordan.cpp
|
../Numerics/Clebsch_Gordan.cpp
|
||||||
../SQD/SphericalPotential/RadialPotential.cpp
|
../SQD/SphericalPotential/RadialPotential.cpp
|
||||||
../SQD/SphericalPotential/ZOverRPotential.cpp
|
../SQD/SphericalPotential/ZOverRPotential.cpp
|
||||||
|
|
|
@ -16,6 +16,7 @@
|
||||||
// -*- C++ -*-
|
// -*- C++ -*-
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include "SQD/SQDFrame.h"
|
#include "SQD/SQDFrame.h"
|
||||||
|
#include "OhmmsApp/ProjectData.h"
|
||||||
|
|
||||||
#ifdef HAVE_QT
|
#ifdef HAVE_QT
|
||||||
#include <qapplication.h>
|
#include <qapplication.h>
|
||||||
|
@ -76,6 +77,9 @@ SQDFrame::solve(const char* fname) {
|
||||||
// build an XML tree from a the file;
|
// build an XML tree from a the file;
|
||||||
m_doc = xmlParseFile(fname);
|
m_doc = xmlParseFile(fname);
|
||||||
if (m_doc == NULL) return false;
|
if (m_doc == NULL) return false;
|
||||||
|
|
||||||
|
//using XPath instead of recursive search
|
||||||
|
xmlXPathContextPtr m_context = xmlXPathNewContext(m_doc);
|
||||||
|
|
||||||
// Check the document is of the right kind
|
// Check the document is of the right kind
|
||||||
xmlNodePtr cur = xmlDocGetRootElement(m_doc);
|
xmlNodePtr cur = xmlDocGetRootElement(m_doc);
|
||||||
|
@ -91,6 +95,23 @@ SQDFrame::solve(const char* fname) {
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
//assign the project id
|
||||||
|
|
||||||
|
//project description
|
||||||
|
OHMMS::ProjectData myProject;
|
||||||
|
|
||||||
|
xmlXPathObjectPtr result
|
||||||
|
= xmlXPathEvalExpression((const xmlChar*)"//project",m_context);
|
||||||
|
if(xmlXPathNodeSetIsEmpty(result->nodesetval)) {
|
||||||
|
WARNMSG("Project is not defined")
|
||||||
|
myProject.reset();
|
||||||
|
} else {
|
||||||
|
myProject.put(result->nodesetval->nodeTab[0]);
|
||||||
|
}
|
||||||
|
xmlXPathFreeObject(result);
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
using namespace ohmmshf;
|
using namespace ohmmshf;
|
||||||
|
|
||||||
HFSolver = new HartreeFock(Pot,Psi);
|
HFSolver = new HartreeFock(Pot,Psi);
|
||||||
|
@ -98,10 +119,11 @@ SQDFrame::solve(const char* fname) {
|
||||||
|
|
||||||
if(!success) {
|
if(!success) {
|
||||||
ERRORMSG("The input file does not confirm. Exit")
|
ERRORMSG("The input file does not confirm. Exit")
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
|
|
||||||
//HFSolver->setRoot(fname);
|
HFSolver->setRoot(myProject.CurrentRoot());
|
||||||
|
|
||||||
success = HFSolver->solve();
|
success = HFSolver->solve();
|
||||||
xmlFreeDoc(m_doc);
|
xmlFreeDoc(m_doc);
|
||||||
xmlCleanupParser();
|
xmlCleanupParser();
|
||||||
|
|
|
@ -71,7 +71,7 @@
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
|
|
||||||
#if !defined(LOGMSG)
|
#if !defined(LOGMSG)
|
||||||
#define LOGMSG(msg) {std::cout<< "QMC " << msg << std::endl;}
|
#define LOGMSG(msg) {std::cout<< "HF " << msg << std::endl;}
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
#if !defined(WARNMSG)
|
#if !defined(WARNMSG)
|
||||||
|
|
|
@ -47,8 +47,8 @@ namespace ohmmshf {
|
||||||
*/
|
*/
|
||||||
|
|
||||||
void HartreeFock::setRoot(const string& aroot) {
|
void HartreeFock::setRoot(const string& aroot) {
|
||||||
LogFileName = aroot;
|
RootFileName = aroot;
|
||||||
LogFileName.append(".log");
|
LogFileName = RootFileName + ".log";
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
|
@ -73,6 +73,7 @@ namespace ohmmshf {
|
||||||
//using XPath instead of recursive search
|
//using XPath instead of recursive search
|
||||||
xmlXPathContextPtr m_context = xmlXPathNewContext(d_root->doc);
|
xmlXPathContextPtr m_context = xmlXPathNewContext(d_root->doc);
|
||||||
|
|
||||||
|
|
||||||
xmlXPathObjectPtr result = xmlXPathEvalExpression((const xmlChar*)"//atom",m_context);
|
xmlXPathObjectPtr result = xmlXPathEvalExpression((const xmlChar*)"//atom",m_context);
|
||||||
if(xmlXPathNodeSetIsEmpty(result->nodesetval)) {
|
if(xmlXPathNodeSetIsEmpty(result->nodesetval)) {
|
||||||
ERRORMSG("Missing Atom information. Exit")
|
ERRORMSG("Missing Atom information. Exit")
|
||||||
|
@ -93,8 +94,8 @@ namespace ohmmshf {
|
||||||
|
|
||||||
XMLReport("Atom name = " << AtomName)
|
XMLReport("Atom name = " << AtomName)
|
||||||
XMLReport("Number of closed shells = " << num_closed_shells)
|
XMLReport("Number of closed shells = " << num_closed_shells)
|
||||||
//initialize xmlNode pointers for the grid, wavefunction and potential
|
//initialize xmlNode pointers for the grid, wavefunction and potential
|
||||||
xmlNodePtr cur1 = cur->xmlChildrenNode;
|
xmlNodePtr cur1 = cur->xmlChildrenNode;
|
||||||
while(cur1 != NULL) {
|
while(cur1 != NULL) {
|
||||||
string cname1((const char*)(cur1->name));
|
string cname1((const char*)(cur1->name));
|
||||||
if(cname1 == "grid") {
|
if(cname1 == "grid") {
|
||||||
|
@ -183,8 +184,8 @@ namespace ohmmshf {
|
||||||
|
|
||||||
xmlXPathFreeContext(m_context);
|
xmlXPathFreeContext(m_context);
|
||||||
|
|
||||||
LogFileName = AtomName;
|
// LogFileName = AtomName;
|
||||||
LogFileName.append(".log");
|
// LogFileName.append(".log");
|
||||||
|
|
||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
@ -293,16 +294,16 @@ namespace ohmmshf {
|
||||||
Psi.put(cur);
|
Psi.put(cur);
|
||||||
|
|
||||||
LOGMSG("Total number of orbitals = " << Psi.size());
|
LOGMSG("Total number of orbitals = " << Psi.size());
|
||||||
// LOGMSG("Orbital | Orbital ID");
|
// LOGMSG("Orbital | Orbital ID");
|
||||||
// for(int j=0; j < Psi.size(); j++){
|
// for(int j=0; j < Psi.size(); j++){
|
||||||
// LOGMSG(j << " " << Psi.ID[j]);
|
// LOGMSG(j << " " << Psi.ID[j]);
|
||||||
// }
|
// }
|
||||||
// LOGMSG("ID | IDcount");
|
// LOGMSG("ID | IDcount");
|
||||||
LOGMSG("(Orbital index, Number of Orbitals)")
|
LOGMSG("(Orbital index, Number of Orbitals)")
|
||||||
for(int j=0; j < Psi.NumUniqueOrb; j++){
|
for(int j=0; j < Psi.NumUniqueOrb; j++){
|
||||||
int id = Psi.ID[j];
|
int id = Psi.ID[j];
|
||||||
LOGMSG("(" << id << ", " << Psi.IDcount[id] << ")");
|
LOGMSG("(" << id << ", " << Psi.IDcount[id] << ")");
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
//return false if there is no wave functions
|
//return false if there is no wave functions
|
||||||
|
@ -548,8 +549,8 @@ namespace ohmmshf {
|
||||||
}
|
}
|
||||||
fout << endl;
|
fout << endl;
|
||||||
}
|
}
|
||||||
Psi.print_HDF5(AtomName,GridType,eigVal);
|
Psi.print_HDF5(RootFileName,GridType,eigVal);
|
||||||
Psi.print_basis(AtomName,GridType);
|
Psi.print_basis(AtomName,RootFileName,GridType);
|
||||||
return max_rad_all;
|
return max_rad_all;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
|
@ -142,7 +142,7 @@ namespace ohmmshf {
|
||||||
}
|
}
|
||||||
} else {
|
} else {
|
||||||
ERRORMSG("The potential and grid do not have proper transformation for numerov")
|
ERRORMSG("The potential and grid do not have proper transformation for numerov")
|
||||||
return false;
|
return false;
|
||||||
}
|
}
|
||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
|
@ -31,9 +31,12 @@ namespace ohmmshf {
|
||||||
///number of closed shells
|
///number of closed shells
|
||||||
int num_closed_shells;
|
int num_closed_shells;
|
||||||
value_type eig_tol, scf_tol, ratio;
|
value_type eig_tol, scf_tol, ratio;
|
||||||
|
///
|
||||||
|
string RootFileName;
|
||||||
string LogFileName;
|
string LogFileName;
|
||||||
string AtomName, PotType, GridType;
|
string AtomName, PotType, GridType;
|
||||||
|
|
||||||
|
|
||||||
xmlNodePtr grid_ptr;
|
xmlNodePtr grid_ptr;
|
||||||
xmlNodePtr pot_ptr;
|
xmlNodePtr pot_ptr;
|
||||||
xmlNodePtr orb_ptr;
|
xmlNodePtr orb_ptr;
|
||||||
|
|
|
@ -1,471 +0,0 @@
|
||||||
//////////////////////////////////////////////////////////////////
|
|
||||||
// (c) Copyright 2003 by Jeongnim Kim and Jordan Vincent
|
|
||||||
//////////////////////////////////////////////////////////////////
|
|
||||||
//////////////////////////////////////////////////////////////////
|
|
||||||
// National Center for Supercomputing Applications &
|
|
||||||
// Materials Computation Center
|
|
||||||
// University of Illinois, Urbana-Champaign
|
|
||||||
// Urbana, IL 61801
|
|
||||||
// e-mail: jnkim@ncsa.uiuc.edu
|
|
||||||
// Tel: 217-244-6319 (NCSA) 217-333-3324 (MCC)
|
|
||||||
//
|
|
||||||
// Supported by
|
|
||||||
// National Center for Supercomputing Applications, UIUC
|
|
||||||
// Materials Computation Center, UIUC
|
|
||||||
//////////////////////////////////////////////////////////////////
|
|
||||||
// -*- C++ -*-
|
|
||||||
#include "AtomicHF/HFConfiguration.h"
|
|
||||||
#include "Numerics/Clebsch_Gordan.h"
|
|
||||||
#include "AtomicHF/fillshells.h"
|
|
||||||
#include "AtomicHF/SphericalPotential/ZOverRPotential.h"
|
|
||||||
#include "AtomicHF/SphericalPotential/HarmonicPotential.h"
|
|
||||||
#include "AtomicHF/SphericalPotential/StepPotential.h"
|
|
||||||
#include "AtomicHF/SphericalPotential/RegularLinearTransform.h"
|
|
||||||
#include "AtomicHF/SphericalPotential/RegularLogTransform.h"
|
|
||||||
#include "AtomicHF/SphericalPotential/NuclearLinearTransform.h"
|
|
||||||
#include "AtomicHF/SphericalPotential/NuclearLogTransform.h"
|
|
||||||
#include "Numerics/Numerov.h"
|
|
||||||
#include "Numerics/RadialFunctorUtility.h"
|
|
||||||
#include "SQD/HartreeFock.h"
|
|
||||||
|
|
||||||
namespace ohmmshf {
|
|
||||||
|
|
||||||
HartreeFock::HartreeFock(RadialPotentialSet& pot,
|
|
||||||
SphericalOrbitalTraits::BasisSetType& psi,
|
|
||||||
const xmlpp::Node* root):
|
|
||||||
Pot(pot), Psi(psi), maxiter(1000), eig_tol(1e-12),
|
|
||||||
scf_tol(1e-8), ratio(0.35)
|
|
||||||
{
|
|
||||||
LogFileName = "atomiHF.log";
|
|
||||||
put(root);
|
|
||||||
}
|
|
||||||
|
|
||||||
void HartreeFock::setRoot(const string& aroot) {
|
|
||||||
LogFileName = aroot;
|
|
||||||
LogFileName.append(".log");
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
@param fake Transformation object
|
|
||||||
@param norb the number of eigen vectors to be obtained
|
|
||||||
@brief perform self-consistent Hartree-Fock calculations
|
|
||||||
*
|
|
||||||
*The first argument is used to tell the compiler which transform object is used but
|
|
||||||
*not meaningful for the calculations. This is to force the compiler inline
|
|
||||||
*everything possible. A better solution will be implement using traits.
|
|
||||||
*/
|
|
||||||
template<class Transform_t>
|
|
||||||
inline void
|
|
||||||
HartreeFock::run(Transform_t* fake, int norb)
|
|
||||||
{
|
|
||||||
typedef Numerov<Transform_t, RadialOrbital_t> Numerov_t;
|
|
||||||
value_type Vtotal,KEnew, KEold,E;
|
|
||||||
value_type lowerbound, upperbound;
|
|
||||||
vector<value_type> energy(Pot.size());
|
|
||||||
eigVal.resize(norb);
|
|
||||||
|
|
||||||
int iter = 0;
|
|
||||||
Vtotal = Pot.evaluate(Psi,energy,norb);
|
|
||||||
Pot.mix(0.0);
|
|
||||||
KEnew = Pot.calcKE(Psi,0,norb);
|
|
||||||
|
|
||||||
string label("spdf");
|
|
||||||
std::ofstream log_stream(LogFileName.c_str());
|
|
||||||
log_stream.precision(8);
|
|
||||||
|
|
||||||
do {
|
|
||||||
KEold = KEnew;
|
|
||||||
value_type eigsum = 0.0;
|
|
||||||
//loop over the orbitals
|
|
||||||
for(int ob=0; ob < norb; ob++){
|
|
||||||
|
|
||||||
//set the number of nodes of the eigen vector
|
|
||||||
Pot.V[ob].setNumOfNodes(Pot.getNumOfNodes(Psi.N[ob],Psi.L[ob]));
|
|
||||||
|
|
||||||
//calculate the lower and upper bounds for the eigenvalue
|
|
||||||
Pot.EnergyBound(lowerbound,upperbound);
|
|
||||||
|
|
||||||
//set up the transformer
|
|
||||||
Transform_t es(Pot.V[ob], Psi.N[ob], Psi.L[ob],Psi.CuspParam, Pot.getMass());
|
|
||||||
|
|
||||||
//initialize the numerov solver
|
|
||||||
Numerov_t numerov(es,Psi(ob));
|
|
||||||
|
|
||||||
//calculate the eigenvalue and the corresponding orbital
|
|
||||||
eigsum += (eigVal[ob] =
|
|
||||||
numerov.solve(lowerbound, upperbound, eig_tol));
|
|
||||||
|
|
||||||
log_stream << Psi.N[ob]<< label[Psi.L[ob]] << '\t' << eigVal[ob] << endl;
|
|
||||||
}
|
|
||||||
log_stream << endl;
|
|
||||||
|
|
||||||
//normalize the orbitals
|
|
||||||
Psi.normalize(norb);
|
|
||||||
//restrict the orbitals
|
|
||||||
Psi.applyRestriction(norb);
|
|
||||||
//calculate the new kinetic energy
|
|
||||||
KEnew = Pot.calcKE(Psi,eigsum,norb);
|
|
||||||
//the total energy
|
|
||||||
E = KEnew + Vtotal;
|
|
||||||
//for the new orbitals Psi, calculate the new SCF potentials
|
|
||||||
Vtotal = Pot.evaluate(Psi,energy,norb);
|
|
||||||
Pot.applyRestriction(Psi);
|
|
||||||
Pot.mix(ratio);
|
|
||||||
log_stream.precision(10);
|
|
||||||
log_stream << "Iteration #" << iter+1 << endl;
|
|
||||||
log_stream << "KE = " << setw(15) << KEnew
|
|
||||||
<< " PE = " << setw(15) << Vtotal << endl;
|
|
||||||
log_stream << "PE/KE = " << setw(15) << Vtotal/KEnew
|
|
||||||
<< " Energy = " << setw(15) << E << endl;
|
|
||||||
log_stream << endl;
|
|
||||||
iter++;
|
|
||||||
///continue the loop until the kinetic energy converges
|
|
||||||
}while(fabs(KEnew-KEold)>scf_tol && iter<maxiter);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
@param pottype the type of potential
|
|
||||||
@param gridtype the type of the radial grid
|
|
||||||
@param norb the number of target eigen vectors
|
|
||||||
@brief Instantiate a transformation function based on the potential and grid type and call run.
|
|
||||||
*
|
|
||||||
*/
|
|
||||||
void HartreeFock::solve(string pottype, string gridtype, int norb) {
|
|
||||||
if(pottype == "Harmonic" || pottype == "Step"){
|
|
||||||
if(gridtype == "linear"){
|
|
||||||
RegularLinearTransform<RadialOrbital_t> *afake=NULL;
|
|
||||||
run(afake,norb);
|
|
||||||
} else if(gridtype == "log"){
|
|
||||||
RegularLogTransform<RadialOrbital_t> *afake=NULL;
|
|
||||||
run(afake,norb);
|
|
||||||
}
|
|
||||||
} else if(pottype == "Nuclear"){
|
|
||||||
if(gridtype == "linear"){
|
|
||||||
NuclearLinearTransform<RadialOrbital_t> *afake=NULL;
|
|
||||||
run(afake,norb);
|
|
||||||
} else if(gridtype == "log"){
|
|
||||||
NuclearLogTransform<RadialOrbital_t> *afake=NULL;
|
|
||||||
run(afake,norb);
|
|
||||||
}
|
|
||||||
} else return;
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
@param q the current xml node which contains the parameter definitions for the eigen solver
|
|
||||||
@return true if successful
|
|
||||||
@brief Set the parameters for the eigen solver
|
|
||||||
*
|
|
||||||
*Aailable parameters
|
|
||||||
* - max_iter, the maximum self-consistent iterations, default=1000
|
|
||||||
* - eig_tol, the eigen-value tolerance, default = \f$1 e^{-12}\f$
|
|
||||||
* - en_tol, the tolerance of the self-consistent loops, default = \f$1 e^{-8}\f$
|
|
||||||
* - mix_ratio, the mixing ratio of the charge density, default = 0.35
|
|
||||||
*/
|
|
||||||
bool HartreeFock::put(const xmlpp::Node* q){
|
|
||||||
|
|
||||||
using namespace xmlpp;
|
|
||||||
NodeSet eset = q->find("./EigenSolve");
|
|
||||||
if(eset.empty()){
|
|
||||||
WARNMSG("Using default values for eigensolver.")
|
|
||||||
XMLReport("maximum iterations = " << maxiter)
|
|
||||||
XMLReport("eigentolerance = " << eig_tol)
|
|
||||||
XMLReport("scftolerance = " << scf_tol)
|
|
||||||
XMLReport("ratio = " << ratio)
|
|
||||||
} else {
|
|
||||||
NodeSet pset = eset[0]->find("./Parameter");
|
|
||||||
for(int i=0; i<pset.size(); i++){
|
|
||||||
Element* cur = dynamic_cast<Element*>(pset[i]);
|
|
||||||
Attribute* att = cur->get_attribute("name");
|
|
||||||
if(att) {
|
|
||||||
const string& aname = att->get_value();
|
|
||||||
xmlNode* ccur = cur->cobj();
|
|
||||||
if(aname == "max_iter"){
|
|
||||||
putContent(maxiter,ccur);
|
|
||||||
XMLReport("maximum iterations = " << maxiter)
|
|
||||||
} else if(aname == "eig_tol"){
|
|
||||||
putContent(eig_tol,ccur);
|
|
||||||
XMLReport("eigentolerance = " << eig_tol)
|
|
||||||
} else if(aname == "en_tol"){
|
|
||||||
putContent(scf_tol,ccur);
|
|
||||||
XMLReport("scftolerance = " << scf_tol)
|
|
||||||
} else if(aname == "mix_ratio"){
|
|
||||||
putContent(ratio,ccur);
|
|
||||||
XMLReport("ratio = " << ratio)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
bool parseXMLFile(RadialPotentialSet& Pot,
|
|
||||||
SphericalOrbitalTraits::BasisSetType& Psi,
|
|
||||||
string& name,
|
|
||||||
string& pottype,
|
|
||||||
string& gridtype,
|
|
||||||
const xmlpp::Node* root) {
|
|
||||||
using namespace xmlpp;
|
|
||||||
int nshells = 0;
|
|
||||||
OneDimGridBase<double>* grid;
|
|
||||||
NodeSet aset = root->find("//Atom");
|
|
||||||
if(aset.empty()){
|
|
||||||
ERRORMSG("Must include atomic information.")
|
|
||||||
return false;
|
|
||||||
} else {
|
|
||||||
Element* cur = dynamic_cast<Element*>(aset[0]);
|
|
||||||
Attribute* att1 = cur->get_attribute("name");
|
|
||||||
Attribute* att2 = cur->get_attribute("num_closed_shells");
|
|
||||||
if(att1){
|
|
||||||
name = att1->get_value();
|
|
||||||
XMLReport("Name = " << name)
|
|
||||||
} else {
|
|
||||||
ERRORMSG("Must specify name.")
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
if(att2){
|
|
||||||
nshells = atoi(att2->get_value().c_str());
|
|
||||||
XMLReport("Number of Closed Shells = " << nshells)
|
|
||||||
} else {
|
|
||||||
WARNMSG("Number of Closed Shells = " << nshells)
|
|
||||||
}
|
|
||||||
|
|
||||||
NodeSet gset = aset[0]->find("./Grid");
|
|
||||||
if(gset.empty()){
|
|
||||||
ERRORMSG("No grid information")
|
|
||||||
return false;
|
|
||||||
} else {
|
|
||||||
Element* p = dynamic_cast<Element*>(gset[0]);
|
|
||||||
Element::AttributeList atts = p->get_attributes();
|
|
||||||
Element::AttributeList::iterator it = atts.begin();
|
|
||||||
|
|
||||||
double scale = 1.0;
|
|
||||||
double min = 0.001;
|
|
||||||
double max = 1000.0;
|
|
||||||
int npts = 2001;
|
|
||||||
while(it != atts.end()) {
|
|
||||||
const string& aname = (*it)->get_name();
|
|
||||||
if(aname == "type") {
|
|
||||||
gridtype = (*it)->get_value();
|
|
||||||
} else if(aname == "scale") {
|
|
||||||
scale = atof((*it)->get_value().c_str());
|
|
||||||
}
|
|
||||||
it++;
|
|
||||||
}
|
|
||||||
|
|
||||||
XMLReport("Grid type = " << gridtype)
|
|
||||||
|
|
||||||
NodeSet gpset = gset[0]->find("./Parameter");
|
|
||||||
if(gpset.empty()){
|
|
||||||
WARNMSG("Using default grid values")
|
|
||||||
} else {
|
|
||||||
for(int i=0; i<gpset.size(); i++){
|
|
||||||
Element* cur = dynamic_cast<Element*>(gpset[i]);
|
|
||||||
Attribute* att = cur->get_attribute("name");
|
|
||||||
if(att) {
|
|
||||||
const string& aname = att->get_value();
|
|
||||||
xmlNode* curc = cur->cobj();
|
|
||||||
if(aname == "min") putContent(min,curc);
|
|
||||||
else if(aname == "max") putContent(max,curc);
|
|
||||||
else if(aname == "npts") putContent(npts,curc);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if(gridtype == "log"){
|
|
||||||
grid = new LogGrid<double>;
|
|
||||||
min/=scale;
|
|
||||||
max/=sqrt(scale);
|
|
||||||
grid->set(min,max,npts);
|
|
||||||
XMLReport("rmin = " << min << ", rmax = " << max
|
|
||||||
// << ", dh = " << grid->dh() << ", npts = " << npts)
|
|
||||||
<< ", npts = " << npts)
|
|
||||||
} else if(gridtype == "linear"){
|
|
||||||
grid = new LinearGrid<double>;
|
|
||||||
min/=scale;
|
|
||||||
max/=sqrt(scale);
|
|
||||||
grid->set(min,max,npts);
|
|
||||||
XMLReport("rmin = " << min << ", rmax = " << max
|
|
||||||
// << ", dh = " << grid->dh() << ", npts = " << npts)
|
|
||||||
<< ", npts = " << npts)
|
|
||||||
} else {
|
|
||||||
ERRORMSG("Grid Type Options: Log or Linear.")
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
// pass grid to wavefunction
|
|
||||||
Psi.m_grid = grid;
|
|
||||||
|
|
||||||
// class to generate Clebsh Gordan coeff.
|
|
||||||
Clebsch_Gordan* CG_coeff = NULL;
|
|
||||||
|
|
||||||
NodeSet potset = aset[0]->find("./Potential");
|
|
||||||
if(potset.empty()){
|
|
||||||
ERRORMSG("Must provide potential information.")
|
|
||||||
return false;
|
|
||||||
} else {
|
|
||||||
Element* cur = dynamic_cast<Element*>(potset[0]);
|
|
||||||
xmlNodePtr cur_sub = cur->cobj()->xmlChildrenNode;
|
|
||||||
while(cur_sub != NULL) {
|
|
||||||
string pname((const char*)(cur_sub->name));
|
|
||||||
if(pname == "Parameter") {
|
|
||||||
xmlAttrPtr att=cur_sub->properties;
|
|
||||||
while(att != NULL) {
|
|
||||||
string vname((const char*)(att->children->content));
|
|
||||||
if(vname == "mass") {
|
|
||||||
double m=1.0;
|
|
||||||
putContent(m,cur_sub);
|
|
||||||
Pot.setMass(m);
|
|
||||||
XMLReport("The effective mass is set to " << Pot.getMass())
|
|
||||||
}
|
|
||||||
att = att->next;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
cur_sub = cur_sub->next;
|
|
||||||
}
|
|
||||||
Attribute* att = cur->get_attribute("type");
|
|
||||||
if(att){
|
|
||||||
pottype = att->get_value();
|
|
||||||
|
|
||||||
NodeSet oset = aset[0]->find("./OrbitalSet");
|
|
||||||
if(oset.empty()){
|
|
||||||
ERRORMSG("Must specify OrbitalSet")
|
|
||||||
return false;
|
|
||||||
} else {
|
|
||||||
Element* cur = dynamic_cast<Element*>(oset[0]);
|
|
||||||
Attribute* att = cur->get_attribute("rest_type");
|
|
||||||
if(att) Psi.Restriction = att->get_value();
|
|
||||||
XMLReport("Orbital restriction type = " << Psi.Restriction)
|
|
||||||
xmlNode* curc = cur->cobj();
|
|
||||||
if(pottype == "Harmonic" || pottype == "Step"){
|
|
||||||
if(nshells) FillShellsHarmPot(Psi,nshells);
|
|
||||||
} else {
|
|
||||||
if(nshells) FillShellsNucPot(Psi,nshells);
|
|
||||||
}
|
|
||||||
Psi.put(curc);
|
|
||||||
Pot.initialize(Psi);
|
|
||||||
LOGMSG("No of Orbs = " << Psi.size());
|
|
||||||
|
|
||||||
LOGMSG("Orbital | Orbital ID");
|
|
||||||
for(int j=0; j < Psi.size(); j++){
|
|
||||||
LOGMSG(j << " " << Psi.ID[j]);
|
|
||||||
}
|
|
||||||
LOGMSG("ID | IDcount");
|
|
||||||
for(int j=0; j < Psi.NumUniqueOrb; j++){
|
|
||||||
LOGMSG(j << " " << Psi.IDcount[j]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
XMLReport("Potential Type = " << pottype)
|
|
||||||
} else {
|
|
||||||
ERRORMSG("Must provide Potential Type.")
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
if(pottype == "Nuclear"){
|
|
||||||
XMLReport("Creating a Nuclear Potential.")
|
|
||||||
double Z;
|
|
||||||
NodeSet gpset = potset[0]->find("./Parameter");
|
|
||||||
if(gpset.empty()){
|
|
||||||
Z = Psi.NumOrb;
|
|
||||||
WARNMSG("Default: Z = " << Z)
|
|
||||||
} else {
|
|
||||||
for(int i=0; i<gpset.size(); i++){
|
|
||||||
Element* cur = dynamic_cast<Element*>(gpset[i]);
|
|
||||||
Attribute* att = cur->get_attribute("name");
|
|
||||||
if(att) {
|
|
||||||
const string& aname = att->get_value();
|
|
||||||
xmlNode* curc = cur->cobj();
|
|
||||||
if(aname == "Z") putContent(Z,curc);
|
|
||||||
XMLReport("Z = " << Z)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// determine maximum angular momentum
|
|
||||||
int lmax = 0;
|
|
||||||
for(int i=0; i<Psi.L.size(); i++)
|
|
||||||
if(Psi.L[i] > lmax) lmax = Psi.L[i];
|
|
||||||
XMLReport("Maximum Angular Momentum = " << lmax)
|
|
||||||
CG_coeff = new Clebsch_Gordan(lmax);
|
|
||||||
Pot.add(new ZOverRPotential(Z), true);
|
|
||||||
Pot.add(new HartreePotential(CG_coeff));
|
|
||||||
Pot.add(new ExchangePotential(CG_coeff));
|
|
||||||
Psi.CuspParam = Z;
|
|
||||||
} // if Nuclear
|
|
||||||
else if(pottype == "Harmonic"){
|
|
||||||
XMLReport("Creating a Harmonic Potential.")
|
|
||||||
double Omega=0.5;
|
|
||||||
NodeSet gpset = potset[0]->find("./Parameter");
|
|
||||||
if(gpset.empty()){
|
|
||||||
Omega = Psi.NumOrb;
|
|
||||||
WARNMSG("Default: Omega = " << Omega)
|
|
||||||
} else {
|
|
||||||
for(int i=0; i<gpset.size(); i++){
|
|
||||||
Element* cur = dynamic_cast<Element*>(gpset[i]);
|
|
||||||
Attribute* att = cur->get_attribute("name");
|
|
||||||
if(att) {
|
|
||||||
const string& aname = att->get_value();
|
|
||||||
xmlNode* curc = cur->cobj();
|
|
||||||
if(aname == "Omega") putContent(Omega,curc);
|
|
||||||
else if(aname == "Z") putContent(Omega,curc);
|
|
||||||
XMLReport("Omega = " << Omega)
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
// determine maximum angular momentum
|
|
||||||
int lmax = 0;
|
|
||||||
for(int i=0; i<Psi.L.size(); i++)
|
|
||||||
if(Psi.L[i] > lmax) lmax = Psi.L[i];
|
|
||||||
lmax++; // increment by 1
|
|
||||||
XMLReport("Maximum Angular Momentum = " << lmax)
|
|
||||||
CG_coeff = new Clebsch_Gordan(lmax);
|
|
||||||
Pot.add(new HarmonicPotential(Omega),true);
|
|
||||||
Pot.add(new HartreePotential(CG_coeff));
|
|
||||||
Pot.add(new ExchangePotential(CG_coeff));
|
|
||||||
Psi.CuspParam = 0.0;
|
|
||||||
} // if Harmonic
|
|
||||||
else if(pottype == "Step"){
|
|
||||||
XMLReport("Creating a step-function Potential.")
|
|
||||||
int lmax = 0;
|
|
||||||
for(int i=0; i<Psi.L.size(); i++)
|
|
||||||
if(Psi.L[i] > lmax) lmax = Psi.L[i];
|
|
||||||
lmax++; // increment by 1
|
|
||||||
XMLReport("Maximum Angular Momentum = " << lmax)
|
|
||||||
CG_coeff = new Clebsch_Gordan(lmax);
|
|
||||||
StepPotential* apot = new StepPotential;
|
|
||||||
apot->put(dynamic_cast<Element*>(potset[0])->cobj());
|
|
||||||
Pot.add(apot,true);
|
|
||||||
Pot.add(new HartreePotential(CG_coeff));
|
|
||||||
Pot.add(new ExchangePotential(CG_coeff));
|
|
||||||
Psi.CuspParam = 0.0;
|
|
||||||
} // steps
|
|
||||||
}
|
|
||||||
}
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// int main(int argc, char **argv) {
|
|
||||||
// using namespace ohmmshf;
|
|
||||||
// using namespace xmlpp;
|
|
||||||
// DomParser parser;
|
|
||||||
// parser.parse_file(argv[1]);
|
|
||||||
// Node* root = parser.get_document()->get_root_node(); //deleted by DomParser.
|
|
||||||
// string element;
|
|
||||||
// string pottype;
|
|
||||||
// string gridtype;
|
|
||||||
// RadialPotentialSet Pot;
|
|
||||||
// SphericalOrbitalTraits::BasisSetType Psi;
|
|
||||||
// parseXMLFile(Pot,Psi,element,pottype,gridtype,root);
|
|
||||||
// HartreeFock HFSolver(Pot,Psi,root);
|
|
||||||
// HFSolver.solve(pottype,gridtype,Psi.size());
|
|
||||||
// Psi.print(element);
|
|
||||||
// return 1;
|
|
||||||
// }
|
|
||||||
|
|
||||||
|
|
||||||
/***************************************************************************
|
|
||||||
* $RCSfile$ $Author$
|
|
||||||
* $Revision$ $Date$
|
|
||||||
* $Id$
|
|
||||||
***************************************************************************/
|
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -1,47 +0,0 @@
|
||||||
#ifndef OHMMS_ATOMICHF_HARTREEFOCK_H
|
|
||||||
#define OHMMS_ATOMICHF_HARTREEFOCK_H
|
|
||||||
#include "AtomicHF/HFConfiguration.h"
|
|
||||||
#include "AtomicHF/RadialPotentialSet.h"
|
|
||||||
#include "OhmmsData/libxmldefs.h"
|
|
||||||
#include <libxml++/libxml++.h>
|
|
||||||
#include <fstream>
|
|
||||||
|
|
||||||
namespace ohmmshf {
|
|
||||||
|
|
||||||
/**Hartree-Fock solver
|
|
||||||
*/
|
|
||||||
struct HartreeFock {
|
|
||||||
|
|
||||||
typedef SphericalOrbitalTraits::BasisSetType BasisSetType;
|
|
||||||
typedef SphericalOrbitalTraits::value_type value_type;
|
|
||||||
typedef SphericalOrbitalTraits::RadialOrbital_t RadialOrbital_t;
|
|
||||||
|
|
||||||
int maxiter;
|
|
||||||
value_type eig_tol, scf_tol, ratio;
|
|
||||||
string LogFileName;
|
|
||||||
std::vector<value_type> eigVal;
|
|
||||||
|
|
||||||
RadialPotentialSet& Pot;
|
|
||||||
BasisSetType& Psi;
|
|
||||||
|
|
||||||
HartreeFock(RadialPotentialSet& pot,
|
|
||||||
BasisSetType& psi,
|
|
||||||
const xmlpp::Node* root);
|
|
||||||
|
|
||||||
bool put(const xmlpp::Node* q);
|
|
||||||
|
|
||||||
void setRoot(const string& aroot);
|
|
||||||
|
|
||||||
template<class Transform_t> void run(Transform_t* fake, int norb);
|
|
||||||
|
|
||||||
void solve(string pottype, string gridtype, int norb);
|
|
||||||
};
|
|
||||||
|
|
||||||
bool parseXMLFile(RadialPotentialSet& Pot,
|
|
||||||
SphericalOrbitalTraits::BasisSetType& Psi,
|
|
||||||
string& name,string& pottype,
|
|
||||||
string& gridtype,
|
|
||||||
const xmlpp::Node* root);
|
|
||||||
}
|
|
||||||
|
|
||||||
#endif
|
|
|
@ -60,7 +60,7 @@ bool YlmRnlSet<GT>::put(xmlNodePtr cur){
|
||||||
cur = cur->xmlChildrenNode;
|
cur = cur->xmlChildrenNode;
|
||||||
while(cur != NULL) {
|
while(cur != NULL) {
|
||||||
if (!(xmlStrcmp(cur->name, (const xmlChar *) "orbital"))) {
|
if (!(xmlStrcmp(cur->name, (const xmlChar *) "orbital"))) {
|
||||||
LOGMSG("Found Orbital");
|
XMLReport("Found Orbital");
|
||||||
n = atoi((const char*)xmlGetProp(cur, (const xmlChar *) "n"));
|
n = atoi((const char*)xmlGetProp(cur, (const xmlChar *) "n"));
|
||||||
if(n < 0) {
|
if(n < 0) {
|
||||||
ERRORMSG("Invalid value for n");
|
ERRORMSG("Invalid value for n");
|
||||||
|
@ -95,7 +95,7 @@ bool YlmRnlSet<GT>::put(xmlNodePtr cur){
|
||||||
NLMS_Map_t::iterator it = OccNo.find(nlms);
|
NLMS_Map_t::iterator it = OccNo.find(nlms);
|
||||||
if(it == OccNo.end()) {
|
if(it == OccNo.end()) {
|
||||||
add(n,l,m,s,occ);
|
add(n,l,m,s,occ);
|
||||||
LOGMSG("Adding Orbital: n=" << n << ", l=" << l <<
|
XMLReport("Adding Orbital: n=" << n << ", l=" << l <<
|
||||||
", m=" << m << ", s=" << s << ", occ=" << occ);
|
", m=" << m << ", s=" << s << ", occ=" << occ);
|
||||||
//OccNo[nlmn] = Num-1;
|
//OccNo[nlmn] = Num-1;
|
||||||
OccNo[nlms] = 1;
|
OccNo[nlms] = 1;
|
||||||
|
@ -112,21 +112,21 @@ bool YlmRnlSet<GT>::put(xmlNodePtr cur){
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
*\param elementName name of the element
|
*\param RootName name of the element
|
||||||
*\param GridType the type of grid
|
*\param GridType the type of grid
|
||||||
*\param eigVal the eigenvalues
|
*\param eigVal the eigenvalues
|
||||||
*\return true if succeeds
|
*\return true if succeeds
|
||||||
*\brief Prints the grid and orbital information to an HDF5
|
*\brief Prints the grid and orbital information to an HDF5
|
||||||
*file named "elementName.h5".
|
*file named "RootName.h5".
|
||||||
*/
|
*/
|
||||||
|
|
||||||
template<class GT>
|
template<class GT>
|
||||||
bool YlmRnlSet<GT>::print_HDF5(const std::string& elementName,
|
bool YlmRnlSet<GT>::print_HDF5(const std::string& RootName,
|
||||||
const std::string& GridType,
|
const std::string& GridType,
|
||||||
const std::vector<value_type>& eigVal){
|
const std::vector<value_type>& eigVal){
|
||||||
|
|
||||||
//print to file fname
|
//print to file fname
|
||||||
string HDFfname = elementName + ".h5";
|
string HDFfname = RootName + ".h5";
|
||||||
hid_t afile = H5Fcreate(HDFfname.c_str(),H5F_ACC_TRUNC,
|
hid_t afile = H5Fcreate(HDFfname.c_str(),H5F_ACC_TRUNC,
|
||||||
H5P_DEFAULT,H5P_DEFAULT);
|
H5P_DEFAULT,H5P_DEFAULT);
|
||||||
|
|
||||||
|
@ -207,14 +207,16 @@ bool YlmRnlSet<GT>::print_HDF5(const std::string& elementName,
|
||||||
|
|
||||||
/**
|
/**
|
||||||
*\param elementName name of the element
|
*\param elementName name of the element
|
||||||
|
*\param RootName name of the element
|
||||||
*\param GridType the type of grid
|
*\param GridType the type of grid
|
||||||
*\return true if succeeds
|
*\return true if succeeds
|
||||||
*\brief Prints the basis information to an xml file named
|
*\brief Prints the basis information to an xml file named
|
||||||
*"elementName.basis.xml" for use in qmcPlusPlus.
|
*"RootName.basis.xml" for use in qmcPlusPlus.
|
||||||
*/
|
*/
|
||||||
|
|
||||||
template<class GT>
|
template<class GT>
|
||||||
bool YlmRnlSet<GT>::print_basis(const std::string& elementName,
|
bool YlmRnlSet<GT>::print_basis(const std::string& elementName,
|
||||||
|
const std::string& RootName,
|
||||||
const std::string& GridType){
|
const std::string& GridType){
|
||||||
|
|
||||||
int nup = 0;
|
int nup = 0;
|
||||||
|
@ -236,13 +238,13 @@ bool YlmRnlSet<GT>::print_basis(const std::string& elementName,
|
||||||
|
|
||||||
Mup = 0.0; Mdown = 0.0;
|
Mup = 0.0; Mdown = 0.0;
|
||||||
|
|
||||||
string fnameXML = elementName + ".basis.xml";
|
string fnameXML = RootName + ".basis.xml";
|
||||||
string fnameHDF5 = elementName + ".h5";
|
string fnameHDF5 = RootName + ".h5";
|
||||||
|
|
||||||
ofstream osXML(fnameXML.c_str());
|
ofstream osXML(fnameXML.c_str());
|
||||||
osXML << "<DeterminantSet type=\"MolecularOrbital\">" << endl;
|
osXML << "<determinantset type=\"MolecularOrbital\">" << endl;
|
||||||
osXML << "<BasisSet>" << endl;
|
osXML << "<basisset>" << endl;
|
||||||
osXML << "<Basis type=\"HFNG\" species=\"" << elementName
|
osXML << "<basis type=\"HFNG\" species=\"" << elementName
|
||||||
<< "\" file=\"" << fnameHDF5 << "\">" << endl;
|
<< "\" file=\"" << fnameHDF5 << "\">" << endl;
|
||||||
|
|
||||||
if(Restriction == "none"){
|
if(Restriction == "none"){
|
||||||
|
@ -272,23 +274,23 @@ bool YlmRnlSet<GT>::print_basis(const std::string& elementName,
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
osXML << "</Basis>" << endl;
|
osXML << "</basis>" << endl;
|
||||||
osXML << "</BasisSet>" << endl;
|
osXML << "</basisset>" << endl;
|
||||||
osXML << "<SlaterDeterminant>" << endl;
|
osXML << "<slaterdeterminant>" << endl;
|
||||||
osXML << "<Determinant spin=\"1\" orbitals=\"" << Mup.rows()
|
osXML << "<determinant spin=\"1cu\" orbitals=\"" << Mup.rows()
|
||||||
<< "\">" << endl;
|
<< "\">" << endl;
|
||||||
osXML << "<Var type=\"Array\">" << endl;
|
osXML << "<parameter id=\"1cu\" type=\"Array\">" << endl;
|
||||||
osXML << Mup;
|
osXML << Mup;
|
||||||
osXML << "</Var>" << endl;
|
osXML << "</parameter>" << endl;
|
||||||
osXML << "</Determinant>" << endl;
|
osXML << "</determinant>" << endl;
|
||||||
osXML << "<Determinant spin=\"-1\" orbitals=\"" << Mdown.rows()
|
osXML << "<determinant spin=\"-1\" orbitals=\"" << Mdown.rows()
|
||||||
<< "\">" << endl;
|
<< "\">" << endl;
|
||||||
osXML << "<Var type=\"Array\">" << endl;
|
osXML << "<parameter id=\"1cd\" type=\"Array\">" << endl;
|
||||||
osXML << Mdown;
|
osXML << Mdown;
|
||||||
osXML << "</Var>" << endl;
|
osXML << "</parameter>" << endl;
|
||||||
osXML << "</Determinant>" << endl;
|
osXML << "</determinant>" << endl;
|
||||||
osXML << "</SlaterDeterminant>" << endl;
|
osXML << "</slaterdeterminant>" << endl;
|
||||||
osXML << "</DeterminantSet>" << endl;
|
osXML << "</determinantset>" << endl;
|
||||||
|
|
||||||
osXML.close();
|
osXML.close();
|
||||||
|
|
||||||
|
|
|
@ -237,7 +237,8 @@ struct YlmRnlSet {
|
||||||
bool get(std::ostream& os);
|
bool get(std::ostream& os);
|
||||||
bool print_HDF5(const std::string&, const std::string&,
|
bool print_HDF5(const std::string&, const std::string&,
|
||||||
const std::vector<value_type>&);
|
const std::vector<value_type>&);
|
||||||
bool print_basis(const std::string&, const std::string&);
|
bool print_basis(const std::string&, const std::string&,
|
||||||
|
const std::string&);
|
||||||
|
|
||||||
};
|
};
|
||||||
#include "SQD/YlmRnlSet.cpp"
|
#include "SQD/YlmRnlSet.cpp"
|
||||||
|
|
Loading…
Reference in New Issue