https://github.com/bdallapi/gpvmc
Tip revision: 19575b3e6fe0b27b464a092d92072e74b9ea59dd authored by Bastien Dalla Piazza on 29 October 2014, 21:50:40 UTC
tuning
tuning
Tip revision: 19575b3
LatticeStepper.cpp
#include "LatticeStepper.h"
#include "SlaterDeterminant.h"
#include "LatticeState.h"
#include "Lattice.h"
#include "WaveFunction.h"
#include "linalg.h"
#include "RanGen.h"
#include "Timer.h"
#include <numeric>
using namespace std;
LatticeStepper::LatticeStepper(LatticeState* latstate, WaveFunction* wav, SlaterDeterminant* amp)
:Stepper(latstate,wav,amp),
m_Nfl(latstate->GetNfl()),
m_Nifs(latstate->GetNifs()),
m_prev(latstate->GetNfl()),
m_flavor_flip(latstate->GetNfl(),false),
m_khop(-1), m_weight(-1)
{}
void LatticeStepper::Reset()
{
m_prev=vector<hop_path_t>(m_Nfl);
m_khop=-1;
m_prev_weight=-1;
m_weight=-1;
}
void LatticeStepper::SetFlavorFlip(const vector<bool>& flavor_flip)
{
m_flavor_flip=flavor_flip;
}
BigDouble LatticeStepper::trystep()
{
#ifdef PROFILE
Timer::tic("LatticeStepper::trystep");
#endif
const Lattice* lat=m_latstate->GetLattice();
const Vertex* v1,*v2;
// it is assumed all vertices hold exactly one fermion
// pick up a random vertex
v1=lat->GetVertices()[size_t(lat->GetNv()*RanGen::uniform())];
// pick up a neighbour at random
v2=v1->edges[size_t(v1->edges.size()*RanGen::uniform())]->GetOther(v1);
vector<uint_vec_t> s1,s2;
m_latstate->GetLatOc(v1->idx,s1);
m_latstate->GetLatOc(v2->idx,s2);
// get s1 occupied flavour
size_t f1=max_element(s1.begin(),s1.end(),
[](const uint_vec_t& a, const uint_vec_t& b)
{return a.size()<b.size();})-s1.begin();
// get s2 occupied flavour
size_t f2=max_element(s2.begin(),s2.end(),
[](const uint_vec_t& a, const uint_vec_t& b)
{return a.size()<b.size();})-s2.begin();
// from here it is assumed there is only one fermion per site.
size_t c= m_Nifs[f1]*RanGen::uniform();
if(!m_flavor_flip[f1] || c==s1[f1][0]){
// exchange neighbouring fermions
m_prev=vector<hop_path_t>(m_Nfl);
if(f1==f2 && s1[f1][0]==s2[f2][0]){
// do nothing if neighbouring fermions are in same state
m_khop=-1;
m_prev_weight=m_weight;
return m_weight;
} else {
m_prev[f1].push_back(hop_t(v1->idx*m_Nifs[f1]+s1[f1][0],
v2->idx*m_Nifs[f1]+s1[f1][0]));
m_prev[f2].push_back(hop_t(v2->idx*m_Nifs[f2]+s2[f2][0],
v1->idx*m_Nifs[f2]+s2[f2][0]));
}
} else {
// flip fermion on v1 to any of the
// other internal degrees of freedom
m_prev=vector<hop_path_t>(m_Nfl);
m_prev[f1].push_back(hop_t(v1->idx*m_Nifs[f1]+s1[f1][0],v1->idx*m_Nifs[f1]+c));
}
const vector<vector<hop_path_t> >& khops=m_wav->GetHops();
vector<BigComplex> qs(khops.size());
#ifdef PROFILE
Timer::tic("LatticeStepper::trystep/VirtUpdate");
#endif
m_amp->VirtUpdate(vector<vector<hop_path_t> >(1,m_prev),khops,qs);
#ifdef PROFILE
Timer::toc("LatticeStepper::trystep/VirtUpdate");
#endif
BigDouble out(0);
for_each(qs.begin(),qs.end(),[&](const BigComplex& aq){out+=norm(aq);});
m_khop=max_element(qs.begin(),qs.end(),
[](const BigComplex& a,const BigComplex& b)
{return norm(a)<norm(b);}) - qs.begin();
m_prev_weight=out;
#ifdef PROFILE
Timer::toc("LatticeStepper::trystep");
#endif
return out;
}
void LatticeStepper::step()
{
static int count=0;
m_latstate->Hop(m_prev);
try{
if(m_khop>=0){
const vector<hop_path_t>& khop_p=m_wav->GetHop(m_khop);
m_wav->Hop(m_khop);
m_amp->Update(m_prev,khop_p);
} else {
m_amp->Update(m_prev,vector<hop_path_t>(m_Nfl));
}
} catch(std::logic_error e){
cout<<*m_latstate<<endl<<*m_wav<<endl;
#ifndef USEMPI
m_latstate->Save();
#endif
throw;
}
if(m_amp->Amp()==0){
cerr<<"LatticeStepper::step(): Stepped to a zero overlap!"<<endl;
abort();
}
BigDouble ow=m_weight;
#ifdef DNDEBUG
if(m_prev_weight>0.0)
m_weight=m_prev_weight;
else {
m_weight=-1;
weight();
}
#else
m_weight=-1;
weight();
if(m_prev_weight>=0.0 && abs(double(m_prev_weight/m_weight)-1)>1e-6){
cerr<<"LatticeStepper::step: error: abs(double(m_prev_weight/m_weight))="<<abs(double(m_prev_weight/m_weight))-1<<" m_prev_weight="<<m_prev_weight<<" and m_weight="<<m_weight<<endl;
abort();
}
#endif
if(abs(double(m_weight/ow))<1e-6){
cerr<<"LatticeStepper::step: warning, abs(double(m_weight/ow)) is VERY small! : "<<abs(double(m_weight/ow))<<endl;
}
m_prev=vector<hop_path_t>(m_Nfl);
m_khop=-1;
count++;
}
BigDouble LatticeStepper::weight()
{
if(m_weight<0.0){
BigDouble out(0);
const vector<vector<hop_path_t> >& khop=m_wav->GetHops();
vector<vector<hop_path_t> > rhop(1,vector<hop_path_t>(m_Nfl));
vector<BigComplex> qs(khop.size(),0);
m_amp->VirtUpdate(rhop,khop,qs);
for(size_t k=0;k<khop.size();++k){
out+=norm(qs[k]);
}
m_weight=out;
if(m_weight==0.0){
std::cout<<"Warning: LatticeStepper::weight: "
"state has no overlap."<<std::endl;
}
return m_weight;
} else {
return m_weight;
}
}
double LatticeStepper::transprob()
{
return 1.0;
}