mirror of https://github.com/QMCPACK/qmcpack.git
adding dynamic structure factor (density(0)density(tau)) hamiltonian element
git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@5600 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Jeremy McMinis
parent
12b5a6541f
commit
d94114a7f4
|
|
@ -14,6 +14,7 @@ SET(HAMSRCS
|
|||
DensityEstimator.cpp
|
||||
SkPot.cpp
|
||||
SkEstimator.cpp
|
||||
DynSkEstimator.cpp
|
||||
MomentumEstimator.cpp
|
||||
ForceBase.cpp
|
||||
)
|
||||
|
|
|
|||
|
|
@ -0,0 +1,219 @@
|
|||
//////////////////////////////////////////////////////////////////
|
||||
// (c) Copyright 2008- by Jeongnim Kim
|
||||
//////////////////////////////////////////////////////////////////
|
||||
//////////////////////////////////////////////////////////////////
|
||||
// National Center for Supercomputing Applications &
|
||||
// Materials Computation Center
|
||||
// University of Illinois, Urbana-Champaign
|
||||
// Urbana, IL 61801
|
||||
// e-mail: jnkim@ncsa.uiuc.edu
|
||||
//
|
||||
// Supported by
|
||||
// National Center for Supercomputing Applications, UIUC
|
||||
// Materials Computation Center, UIUC
|
||||
//////////////////////////////////////////////////////////////////
|
||||
// -*- C++ -*-
|
||||
#include <QMCHamiltonians/DynSkEstimator.h>
|
||||
#include <LongRange/StructFact.h>
|
||||
#include <Utilities/IteratorUtility.h>
|
||||
#include <OhmmsData/AttributeSet.h>
|
||||
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
DynSkEstimator::DynSkEstimator(ParticleSet& source)
|
||||
{
|
||||
sourcePtcl = &source;
|
||||
UpdateMode.set(COLLECTABLE,1);
|
||||
NumSpecies=source.getSpeciesSet().getTotalNum();
|
||||
NumK=source.SK->KLists.numk;
|
||||
OneOverN=1.0/static_cast<RealType>(source.getTotalNum());
|
||||
Kshell=source.SK->KLists.kshell;
|
||||
MaxKshell=Kshell.size()-1;
|
||||
RhokTot.resize(NumK);
|
||||
values.resize(NumK);
|
||||
|
||||
Kmag.resize(MaxKshell);
|
||||
OneOverDnk.resize(MaxKshell);
|
||||
for(int ks=0, k=0; ks<MaxKshell; ks++)
|
||||
{
|
||||
Kmag[ks]=std::sqrt(source.SK->KLists.ksq[Kshell[ks]]);
|
||||
OneOverDnk[ks]=1.0/static_cast<RealType>(Kshell[ks+1]-Kshell[ks]);
|
||||
}
|
||||
hdf5_out=false;
|
||||
}
|
||||
|
||||
void DynSkEstimator::resetTargetParticleSet(ParticleSet& P)
|
||||
{
|
||||
sourcePtcl = &P;
|
||||
}
|
||||
|
||||
DynSkEstimator::Return_t DynSkEstimator::calculate(ParticleSet& P)
|
||||
{
|
||||
for(int j(0);j<NumT-MinT;j++)
|
||||
{
|
||||
// j steps ago, all pindex synched
|
||||
int oindex=tWalker->PHindex[pindx]-1-MinT;
|
||||
if(oindex<0)
|
||||
oindex+=NumT;
|
||||
cindx = oindex-j;
|
||||
if(cindx<0)
|
||||
cindx+=NumT;
|
||||
|
||||
for(int i=0; i<NumK; ++i)
|
||||
{
|
||||
RealType then_real=tWalker->PropertyHistory[pindx+2*i][cindx];
|
||||
RealType now_real=tWalker->PropertyHistory[pindx+2*i][oindex];
|
||||
|
||||
RealType then_imag=tWalker->PropertyHistory[pindx+2*i+1][cindx];
|
||||
RealType now_imag=tWalker->PropertyHistory[pindx+2*i+1][oindex];
|
||||
|
||||
values[i+NumK*j]=OneOverN*(now_real*then_real+now_imag*then_imag );
|
||||
}
|
||||
}
|
||||
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
|
||||
DynSkEstimator::Return_t DynSkEstimator::evaluate(ParticleSet& P)
|
||||
{
|
||||
//sum over species for t=0 (now)
|
||||
std::copy(P.SK->rhok[0],P.SK->rhok[0]+NumK,RhokTot.begin());
|
||||
for(int i=1; i<NumSpecies; ++i)
|
||||
accumulate_elements(P.SK->rhok[i],P.SK->rhok[i]+NumK,RhokTot.begin());
|
||||
|
||||
for(int i=0; i<NumK; ++i)
|
||||
{
|
||||
tWalker->addPropertyHistoryPoint(pindx+2*i , RhokTot[i].real());
|
||||
tWalker->addPropertyHistoryPoint(pindx+2*i+1, RhokTot[i].imag());
|
||||
}
|
||||
|
||||
calculate(P);
|
||||
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
void DynSkEstimator::addObservables(PropertySetType& plist, BufferType& collectables)
|
||||
{
|
||||
// if(hdf5_out)
|
||||
// {
|
||||
// myIndex=collectables.size();
|
||||
// vector<RealType> tmp(NumK);
|
||||
// collectables.add(tmp.begin(),tmp.end());
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
myIndex=plist.size();
|
||||
for (int i=MinT;i<NumT;i++)
|
||||
{
|
||||
for(int j(0);j<NumK;j++)
|
||||
{
|
||||
std::stringstream sstr;
|
||||
sstr << "dsk_" <<j<<"_"<<i;
|
||||
int id=plist.add(sstr.str());
|
||||
}
|
||||
}
|
||||
// }
|
||||
}
|
||||
|
||||
void DynSkEstimator::addObservables(PropertySetType& plist )
|
||||
{
|
||||
myIndex=plist.size();
|
||||
for (int i=MinT;i<NumT;i++)
|
||||
{
|
||||
for(int j(0);j<NumK;j++)
|
||||
{
|
||||
std::stringstream sstr;
|
||||
sstr << "dsk_" <<j<<"_"<<i;
|
||||
int id=plist.add(sstr.str());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void DynSkEstimator::setObservables(PropertySetType& plist)
|
||||
{
|
||||
if (!hdf5_out)
|
||||
std::copy(values.begin(),values.end(),plist.begin()+myIndex);
|
||||
}
|
||||
|
||||
void DynSkEstimator::setParticlePropertyList(PropertySetType& plist
|
||||
, int offset)
|
||||
{
|
||||
if (!hdf5_out)
|
||||
std::copy(values.begin(),values.end(),plist.begin()+myIndex+offset);
|
||||
}
|
||||
|
||||
|
||||
void DynSkEstimator::registerCollectables(vector<observable_helper*>& h5desc
|
||||
, hid_t gid) const
|
||||
{
|
||||
// if (hdf5_out)
|
||||
// {
|
||||
// vector<int> ndim(1,NumK);
|
||||
// observable_helper* h5o=new observable_helper(myName);
|
||||
// h5o->set_dimensions(ndim,myIndex);
|
||||
// h5o->open(gid);
|
||||
// h5desc.push_back(h5o);
|
||||
//
|
||||
// hsize_t kdims[2];
|
||||
// kdims[0] = NumK;
|
||||
// kdims[1] = OHMMS_DIM;
|
||||
// string kpath = myName + "/kpoints";
|
||||
// hid_t k_space = H5Screate_simple(2,kdims, NULL);
|
||||
// hid_t k_set = H5Dcreate (gid, kpath.c_str(), H5T_NATIVE_DOUBLE, k_space, H5P_DEFAULT);
|
||||
// hid_t mem_space = H5Screate_simple (2, kdims, NULL);
|
||||
// double *ptr = &(sourcePtcl->SK->KLists.kpts_cart[0][0]);
|
||||
// herr_t ret = H5Dwrite(k_set, H5T_NATIVE_DOUBLE, mem_space, k_space, H5P_DEFAULT, ptr);
|
||||
// H5Dclose (k_set);
|
||||
// H5Sclose (mem_space);
|
||||
// H5Sclose (k_space);
|
||||
// H5Fflush(gid, H5F_SCOPE_GLOBAL);
|
||||
// }
|
||||
}
|
||||
|
||||
bool DynSkEstimator::putSpecial(xmlNodePtr cur, ParticleSet& P)
|
||||
{
|
||||
string tagName("dSK");
|
||||
NumT=10;
|
||||
MinT=0;
|
||||
|
||||
OhmmsAttributeSet Tattrib;
|
||||
Tattrib.add(tagName,"name");
|
||||
Tattrib.add(NumT,"max");
|
||||
Tattrib.add(MinT,"min");
|
||||
Tattrib.put(cur);
|
||||
|
||||
app_log()<<" "<<NumT<<" dynamic SK values will be calculated"<<endl;
|
||||
|
||||
// property history is real, so we must add two (complex)
|
||||
pindx = P.addPropertyHistory(NumT);
|
||||
for(int j(1);j<NumK*2;j++)
|
||||
P.addPropertyHistory(NumT);
|
||||
|
||||
// we need NumK*NumT values for all observables
|
||||
values.resize(NumK*(NumT-MinT));
|
||||
cindx=0;
|
||||
return true;
|
||||
}
|
||||
|
||||
bool DynSkEstimator::get(std::ostream& os) const
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
QMCHamiltonianBase* DynSkEstimator::makeClone(ParticleSet& qp
|
||||
, TrialWaveFunction& psi)
|
||||
{
|
||||
DynSkEstimator* myclone = new DynSkEstimator(*this);
|
||||
// myclone->hdf5_out=hdf5_out;
|
||||
// myclone->myIndex=myIndex;
|
||||
return myclone;
|
||||
}
|
||||
}
|
||||
|
||||
/***************************************************************************
|
||||
* $RCSfile$ $Author: jnkim $
|
||||
* $Revision: 2945 $ $Date: 2008-08-05 10:21:33 -0500 (Tue, 05 Aug 2008) $
|
||||
* $Id: ForceBase.h 2945 2008-08-05 15:21:33Z jnkim $
|
||||
***************************************************************************/
|
||||
|
|
@ -0,0 +1,110 @@
|
|||
//////////////////////////////////////////////////////////////////
|
||||
// (c) Copyright 2008- by Jeongnim Kim
|
||||
//////////////////////////////////////////////////////////////////
|
||||
//////////////////////////////////////////////////////////////////
|
||||
// National Center for Supercomputing Applications &
|
||||
// Materials Computation Center
|
||||
// University of Illinois, Urbana-Champaign
|
||||
// Urbana, IL 61801
|
||||
// e-mail: jnkim@ncsa.uiuc.edu
|
||||
//
|
||||
// Supported by
|
||||
// National Center for Supercomputing Applications, UIUC
|
||||
// Materials Computation Center, UIUC
|
||||
//////////////////////////////////////////////////////////////////
|
||||
// -*- C++ -*-
|
||||
/** @file DynSkEstimator.h
|
||||
* @brief Declare DynSkEstimator
|
||||
*/
|
||||
#ifndef QMCPLUSPLUS_DYN_SK_ESTIMATOR_H
|
||||
#define QMCPLUSPLUS_DYN_SK_ESTIMATOR_H
|
||||
#include <QMCHamiltonians/QMCHamiltonianBase.h>
|
||||
namespace qmcplusplus
|
||||
{
|
||||
|
||||
/** DynSkEstimator evaluate the dynamic structure factor of the target particleset
|
||||
*
|
||||
* <estimator name="dsk" type="dynsk" debug="no"/>
|
||||
*/
|
||||
class DynSkEstimator: public QMCHamiltonianBase
|
||||
{
|
||||
public:
|
||||
|
||||
DynSkEstimator(ParticleSet& elns);
|
||||
|
||||
void resetTargetParticleSet(ParticleSet& P);
|
||||
|
||||
Return_t evaluate(ParticleSet& P);
|
||||
|
||||
Return_t calculate(ParticleSet& P);
|
||||
|
||||
inline Return_t evaluate(ParticleSet& P, vector<NonLocalData>& Txy)
|
||||
{
|
||||
return evaluate(P);
|
||||
}
|
||||
|
||||
inline Return_t rejectedMove(ParticleSet& P)
|
||||
{
|
||||
int lastindex = tWalker->PHindex[pindx]-1;
|
||||
if (lastindex<0) lastindex += NumT;
|
||||
|
||||
for (int i=0;i<NumK*2;i++)
|
||||
tWalker->addPropertyHistoryPoint(pindx+i, tWalker->PropertyHistory[pindx+i][lastindex]);
|
||||
|
||||
calculate(P);
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
void addObservables(PropertySetType& plist);
|
||||
void addObservables(PropertySetType& plist, BufferType& collectables);
|
||||
void registerCollectables(vector<observable_helper*>& h5desc, hid_t gid) const ;
|
||||
void setObservables(PropertySetType& plist);
|
||||
void setParticlePropertyList(PropertySetType& plist, int offset);
|
||||
bool putSpecial(xmlNodePtr cur, ParticleSet& P);
|
||||
bool put(xmlNodePtr cur) {};
|
||||
bool get(std::ostream& os) const;
|
||||
QMCHamiltonianBase* makeClone(ParticleSet& qp, TrialWaveFunction& psi);
|
||||
|
||||
protected:
|
||||
ParticleSet *sourcePtcl;
|
||||
/** number of species */
|
||||
int NumSpecies;
|
||||
/** number of kpoints */
|
||||
int NumK;
|
||||
/** number of kshells */
|
||||
int MaxKshell;
|
||||
/** normalization factor */
|
||||
RealType OneOverN;
|
||||
/** kshell counters */
|
||||
vector<int> Kshell;
|
||||
/** instantaneous structure factor */
|
||||
vector<RealType> Kmag;
|
||||
/** 1.0/degenracy for a ksell */
|
||||
vector<RealType> OneOverDnk;
|
||||
/** \f$rho_k = \sum_{\alpha} \rho_k^{\alpha} \f$ for species index \f$\alpha\f$ */
|
||||
Vector<ComplexType> RhokTot;
|
||||
Vector<RealType> values;
|
||||
/** resize the internal data
|
||||
*
|
||||
* The argument list is not completed
|
||||
*/
|
||||
void resize();
|
||||
|
||||
bool hdf5_out;
|
||||
|
||||
// storing old rhoK
|
||||
/** starting index in walker of stored rho values, last index(circular queue **/
|
||||
int pindx,cindx;
|
||||
/** length of stored rhoK **/
|
||||
int NumT,MinT;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
/***************************************************************************
|
||||
* $RCSfile$ $Author: jnkim $
|
||||
* $Revision: 2945 $ $Date: 2008-08-05 10:21:33 -0500 (Tue, 05 Aug 2008) $
|
||||
* $Id: ForceBase.h 2945 2008-08-05 15:21:33Z jnkim $
|
||||
***************************************************************************/
|
||||
|
|
@ -38,6 +38,7 @@
|
|||
#include "QMCHamiltonians/LocalMomentEstimator.h"
|
||||
#include "QMCHamiltonians/DensityEstimator.h"
|
||||
#include "QMCHamiltonians/SkEstimator.h"
|
||||
#include "QMCHamiltonians/DynSkEstimator.h"
|
||||
#include "QMCHamiltonians/MomentumEstimator.h"
|
||||
#if OHMMS_DIM == 3
|
||||
#include "QMCHamiltonians/LocalCorePolPotential.h"
|
||||
|
|
@ -349,7 +350,7 @@ namespace qmcplusplus {
|
|||
apot->put(cur);
|
||||
targetH->addOperator(apot,potName,false);
|
||||
}
|
||||
else if(potType == "localmoment")
|
||||
else if(potType == "localmoment")
|
||||
{
|
||||
string SourceName = "ion0";
|
||||
OhmmsAttributeSet hAttrib;
|
||||
|
|
@ -383,7 +384,17 @@ namespace qmcplusplus {
|
|||
targetH->addOperator(apot,potName,false);
|
||||
}
|
||||
}
|
||||
else if(potType == "sk")
|
||||
else if(potType == "dynsk")
|
||||
{
|
||||
if(PBCType)//only if perioidic
|
||||
{
|
||||
DynSkEstimator* apot=new DynSkEstimator(*targetPtcl);
|
||||
apot->putSpecial(cur,*targetPtcl);
|
||||
targetH->addOperator(apot,potName,false);
|
||||
app_log()<<"Adding dynamic S(k) estimator"<<endl;
|
||||
}
|
||||
}
|
||||
else if(potType == "sk")
|
||||
{
|
||||
if(PBCType)//only if perioidic
|
||||
{
|
||||
|
|
|
|||
Loading…
Reference in New Issue