本文整理汇总了C++中SweepParams::get_keep_states方法的典型用法代码示例。如果您正苦于以下问题:C++ SweepParams::get_keep_states方法的具体用法?C++ SweepParams::get_keep_states怎么用?C++ SweepParams::get_keep_states使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类SweepParams
的用法示例。
在下文中一共展示了SweepParams::get_keep_states方法的8个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: do_one
double SweepTwopdm::do_one(SweepParams &sweepParams, const bool &warmUp, const bool &forward, const bool &restart, const int &restartSize, int state)
{
Timer sweeptimer;
int integralIndex = 0;
if (dmrginp.hamiltonian() == BCS) {
pout << "Two PDM with BCS calculations is not implemented" << endl;
exit(0);
}
pout.precision(12);
SpinBlock system;
const int nroots = dmrginp.nroots();
std::vector<double> finalEnergy(nroots,0.);
std::vector<double> finalEnergy_spins(nroots,0.);
double finalError = 0.;
sweepParams.set_sweep_parameters();
// a new renormalisation sweep routine
pout << ((forward) ? "\t\t\t Starting renormalisation sweep in forwards direction" : "\t\t\t Starting renormalisation sweep in backwards direction") << endl;
pout << "\t\t\t ============================================================================ " << endl;
InitBlocks::InitStartingBlock (system,forward, sweepParams.current_root(), sweepParams.current_root(), sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), restartSize, restart, warmUp, integralIndex);
if(!restart)
sweepParams.set_block_iter() = 0;
pout << "\t\t\t Starting block is :: " << endl << system << endl;
if (!restart)
SpinBlock::store (forward, system.get_sites(), system, sweepParams.current_root(), sweepParams.current_root()); // if restart, just restoring an existing block --
sweepParams.savestate(forward, system.get_sites().size());
bool dot_with_sys = true;
array_4d<double> twopdm(2*dmrginp.last_site(), 2*dmrginp.last_site(), 2*dmrginp.last_site(), 2*dmrginp.last_site());
twopdm.Clear();
save_twopdm_binary(twopdm, state, state);
for (; sweepParams.get_block_iter() < sweepParams.get_n_iters(); )
{
pout << "\n\t\t\t Block Iteration :: " << sweepParams.get_block_iter() << endl;
pout << "\t\t\t ----------------------------" << endl;
if (forward)
p1out << "\t\t\t Current direction is :: Forwards " << endl;
else
p1out << "\t\t\t Current direction is :: Backwards " << endl;
//if (SHOW_MORE) pout << "system block" << endl << system << endl;
if (dmrginp.no_transform())
sweepParams.set_guesstype() = BASIC;
else if (!warmUp && sweepParams.get_block_iter() != 0)
sweepParams.set_guesstype() = TRANSFORM;
else if (!warmUp && sweepParams.get_block_iter() == 0 &&
((dmrginp.algorithm_method() == TWODOT_TO_ONEDOT && dmrginp.twodot_to_onedot_iter() != sweepParams.get_sweep_iter()) ||
dmrginp.algorithm_method() != TWODOT_TO_ONEDOT))
sweepParams.set_guesstype() = TRANSPOSE;
else
sweepParams.set_guesstype() = BASIC;
p1out << "\t\t\t Blocking and Decimating " << endl;
SpinBlock newSystem;
BlockAndDecimate (sweepParams, system, newSystem, warmUp, dot_with_sys, state);
for(int j=0;j<nroots;++j)
pout << "\t\t\t Total block energy for State [ " << j <<
" ] with " << sweepParams.get_keep_states()<<" :: " << sweepParams.get_lowest_energy()[j] <<endl;
finalEnergy_spins = ((sweepParams.get_lowest_energy()[0] < finalEnergy[0]) ? sweepParams.get_lowest_energy_spins() : finalEnergy_spins);
finalEnergy = ((sweepParams.get_lowest_energy()[0] < finalEnergy[0]) ? sweepParams.get_lowest_energy() : finalEnergy);
finalError = max(sweepParams.get_lowest_error(),finalError);
system = newSystem;
pout << system<<endl;
SpinBlock::store (forward, system.get_sites(), system, sweepParams.current_root(), sweepParams.current_root());
p1out << "\t\t\t saving state " << system.get_sites().size() << endl;
++sweepParams.set_block_iter();
//sweepParams.savestate(forward, system.get_sites().size());
}
//for(int j=0;j<nroots;++j)
{int j = state;
pout << "\t\t\t Finished Sweep with " << sweepParams.get_keep_states() << " states and sweep energy for State [ " << j
<< " ] with Spin [ " << dmrginp.molecule_quantum().get_s() << " ] :: " << finalEnergy[j] << endl;
}
pout << "\t\t\t Largest Error for Sweep with " << sweepParams.get_keep_states() << " states is " << finalError << endl;
pout << "\t\t\t ============================================================================ " << endl;
int i = state, j = state;
//for (int j=0; j<=i; j++) {
load_twopdm_binary(twopdm, i, j);
//calcenergy(twopdm, i);
save_twopdm_text(twopdm, i, j);
save_spatial_twopdm_text(twopdm, i, j);
save_spatial_twopdm_binary(twopdm, i, j);
// update the static number of iterations
//.........这里部分代码省略.........
示例2: BlockAndDecimate
void SweepTwopdm::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, int state)
{
//mcheck("at the start of block and decimate");
// figure out if we are going forward or backwards
dmrginp.guessgenT -> start();
bool forward = (system.get_sites() [0] == 0);
SpinBlock systemDot;
SpinBlock envDot;
int systemDotStart, systemDotEnd;
int systemDotSize = sweepParams.get_sys_add() - 1;
if (forward)
{
systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
systemDotEnd = systemDotStart + systemDotSize;
}
else
{
systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
systemDotEnd = systemDotStart - systemDotSize;
}
vector<int> spindotsites(2);
spindotsites[0] = systemDotStart;
spindotsites[1] = systemDotEnd;
//if (useSlater) {
systemDot = SpinBlock(systemDotStart, systemDotEnd, system.get_integralIndex(), true);
//SpinBlock::store(true, systemDot.get_sites(), systemDot);
//}
//else
//SpinBlock::restore(true, spindotsites, systemDot);
SpinBlock environment, environmentDot, newEnvironment;
int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;
const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();
system.addAdditionalCompOps();
InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, sweepParams.current_root(), sweepParams.current_root(), sweepParams.get_sys_add(), dmrginp.direct(), system.get_integralIndex(), DISTRIBUTED_STORAGE, true, true);
InitBlocks::InitNewEnvironmentBlock(environment, systemDot, newEnvironment, system, systemDot, sweepParams.current_root(), sweepParams.current_root(),
sweepParams.get_sys_add(), sweepParams.get_env_add(), forward, dmrginp.direct(),
sweepParams.get_onedot(), nexact, useSlater, system.get_integralIndex(), true, true, true);
SpinBlock big;
newSystem.set_loopblock(true);
system.set_loopblock(false);
newEnvironment.set_loopblock(false);
InitBlocks::InitBigBlock(newSystem, newEnvironment, big);
const int nroots = dmrginp.nroots();
std::vector<Wavefunction> solution(1);
DiagonalMatrix e;
GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, true, 0.0);
#ifndef SERIAL
mpi::communicator world;
mpi::broadcast(world, solution, 0);
#endif
std::vector<Matrix> rotateMatrix;
DensityMatrix tracedMatrix(newSystem.get_stateInfo());
tracedMatrix.allocate(newSystem.get_stateInfo());
tracedMatrix.makedensitymatrix(solution, big, std::vector<double>(1,1.0), 0.0, 0.0, false);
rotateMatrix.clear();
if (!mpigetrank())
double error = makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
#ifndef SERIAL
mpi::broadcast(world,rotateMatrix,0);
#endif
#ifdef SERIAL
const int numprocs = 1;
#endif
#ifndef SERIAL
const int numprocs = world.size();
#endif
if (sweepParams.get_block_iter() == 0)
compute_twopdm_initial(solution, system, systemDot, newSystem, newEnvironment, big, numprocs, state);
compute_twopdm_sweep(solution, system, systemDot, newSystem, newEnvironment, big, numprocs, state);
if (sweepParams.get_block_iter() == sweepParams.get_n_iters() - 1)
compute_twopdm_final(solution, system, systemDot, newSystem, newEnvironment, big, numprocs, state);
SaveRotationMatrix (newSystem.get_sites(), rotateMatrix, state);
//for(int i=0;i<dmrginp.nroots();++i)
solution[0].SaveWavefunctionInfo (big.get_stateInfo(), big.get_leftBlock()->get_sites(), state);
newSystem.transform_operators(rotateMatrix);
}
示例3: dmrg
void dmrg(double sweep_tol)
{
double last_fe = 10.e6;
double last_be = 10.e6;
double old_fe = 0.;
double old_be = 0.;
SweepParams sweepParams;
int old_states=sweepParams.get_keep_states();
int new_states;
double old_error=0.0;
double old_energy=0.0;
// warm up sweep ...
bool dodiis = false;
int domoreIter = 0;
bool direction;
//this is regular dmrg calculation
if(!dmrginp.setStateSpecific()) {
sweepParams.current_root() = -1;
last_fe = Sweep::do_one(sweepParams, true, true, false, 0);
direction = false;
while ((fabs(last_fe - old_fe) > sweep_tol) || (fabs(last_be - old_be) > sweep_tol) ||
(dmrginp.algorithm_method() == TWODOT_TO_ONEDOT && dmrginp.twodot_to_onedot_iter()+1 >= sweepParams.get_sweep_iter()) )
{
old_fe = last_fe;
old_be = last_be;
if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
break;
last_be = Sweep::do_one(sweepParams, false, false, false, 0);
direction = true;
pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
break;
//For obtaining the extrapolated energy
old_states=sweepParams.get_keep_states();
new_states=sweepParams.get_keep_states_ls();
last_fe = Sweep::do_one(sweepParams, false, true, false, 0);
direction = false;
new_states=sweepParams.get_keep_states();
pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
if (domoreIter == 2) {
dodiis = true;
break;
}
}
}
else { //this is state specific calculation
const int nroots = dmrginp.nroots();
bool direction=true;
int restartsize;
//sweepParams.restorestate(direction, restartsize);
//sweepParams.set_sweep_iter() = 0;
//sweepParams.set_restart_iter() = 0;
algorithmTypes atype;
pout << "STARTING STATE SPECIFIC CALCULATION "<<endl;
for (int i=0; i<nroots; i++) {
atype = dmrginp.algorithm_method();
dmrginp.set_algorithm_method() = ONEDOT;
sweepParams.current_root() = i;
p1out << "RUNNING GENERATE BLOCKS FOR STATE "<<i<<endl;
if (mpigetrank()==0) {
Sweep::InitializeStateInfo(sweepParams, direction, i);
Sweep::InitializeStateInfo(sweepParams, !direction, i);
Sweep::CanonicalizeWavefunction(sweepParams, direction, i);
Sweep::CanonicalizeWavefunction(sweepParams, !direction, i);
Sweep::CanonicalizeWavefunction(sweepParams, direction, i);
Sweep::InitializeStateInfo(sweepParams, direction, i);
Sweep::InitializeStateInfo(sweepParams, !direction, i);
}
for (int j=0; j<i ; j++) {
int integralIndex = 0;
Sweep::InitializeOverlapSpinBlocks(sweepParams, direction, i, j, integralIndex);
Sweep::InitializeOverlapSpinBlocks(sweepParams, !direction, i, j, integralIndex);
}
dmrginp.set_algorithm_method() = atype;
p1out << "RUNNING GENERATE BLOCKS FOR STATE "<<i<<endl;
SweepGenblock::do_one(sweepParams, false, !direction, false, 0, i, i);
sweepParams.set_sweep_iter() = 0;
sweepParams.set_restart_iter() = 0;
sweepParams.savestate(!direction, restartsize);
pout << "STATE SPECIFIC CALCULATION FOR STATE: "<<i<<endl;
//.........这里部分代码省略.........
示例4: dmrg_stateSpecific
void dmrg_stateSpecific(double sweep_tol, int targetState)
{
double last_fe = 10.e6;
double last_be = 10.e6;
double old_fe = 0.;
double old_be = 0.;
int ls_count=0;
SweepParams sweepParams;
int old_states=sweepParams.get_keep_states();
int new_states;
double old_error=0.0;
double old_energy=0.0;
// warm up sweep ...
bool direction;
int restartsize;
sweepParams.restorestate(direction, restartsize);
//initialize array of size m_maxiter or dmrginp.max_iter() for dw and energy
sweepParams.current_root() = targetState;
last_fe = Sweep::do_one(sweepParams, false, direction, true, restartsize);
while ((fabs(last_fe - old_fe) > sweep_tol) || (fabs(last_be - old_be) > sweep_tol) )
{
old_fe = last_fe;
old_be = last_be;
if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
break;
last_be = Sweep::do_one(sweepParams, false, !direction, false, 0);
pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
break;
last_fe = Sweep::do_one(sweepParams, false, direction, false, 0);
new_states=sweepParams.get_keep_states();
pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
}
pout << "Converged Energy " << sweepParams.get_lowest_energy()[0]<< std::endl;
if(dmrginp.max_iter() <= sweepParams.get_sweep_iter()) {
pout << "Maximum sweep iterations achieved " << std::endl;
}
//one has to canonicalize the wavefunction with atleast 3 sweeps, this is a quirk of the way
//we transform wavefunction
if (mpigetrank()==0) {
Sweep::InitializeStateInfo(sweepParams, !direction, targetState);
Sweep::InitializeStateInfo(sweepParams, direction, targetState);
Sweep::CanonicalizeWavefunction(sweepParams, !direction, targetState);
Sweep::CanonicalizeWavefunction(sweepParams, direction, targetState);
Sweep::CanonicalizeWavefunction(sweepParams, !direction, targetState);
Sweep::InitializeStateInfo(sweepParams, !direction, targetState);
Sweep::InitializeStateInfo(sweepParams, direction, targetState);
}
}
示例5: dotSystem
void SpinAdapted::mps_nevpt::type1::Startup(const SweepParams &sweepParams, const bool &forward, perturber& pb, int baseState) {
#ifndef SERIAL
mpi::communicator world;
#endif
assert(forward);
SpinBlock system;
system.nonactive_orb() =pb.orb();
bool restart=false, warmUp = false;
int forward_starting_size=1, backward_starting_size=0, restartSize =0;
InitBlocks::InitStartingBlock(system, forward, pb.wavenumber(), baseState, forward_starting_size, backward_starting_size, restartSize, restart, warmUp, 0,pb.braquanta, pb.ketquanta);
SpinBlock::store (forward, system.get_sites(), system, pb.wavenumber(), baseState); // if restart, just restoring an existing block --
for (int i=0; i<mps_nevpt::sweepIters; i++) {
SpinBlock newSystem;
SpinBlock dotSystem(i+1,i+1,pb.orb(),false);
system.addAdditionalCompOps();
//newSystem.default_op_components(true, system, dotSystem, true, true, false);
newSystem.perturb_op_components(false, system, dotSystem, pb);
newSystem.setstoragetype(DISTRIBUTED_STORAGE);
newSystem.BuildSumBlock(LessThanQ, system, dotSystem, pb.braquanta, pb.ketquanta);
newSystem.printOperatorSummary();
//SpinBlock Environment, big;
//SpinBlock::restore (!forward, newSystem.get_complementary_sites() , Environment, baseState, baseState);
//TODO
//SpinBlock::restore (!forward, newSystem.get_complementary_sites() , Environment,sweepParams.current_root(),sweepParams.current_root());
//big.BuildSumBlock(PARTICLE_SPIN_NUMBER_CONSTRAINT, newSystem, Environment, pb.braquanta, pb.ketquanta);
//StateInfo envStateInfo;
StateInfo ketStateInfo;
StateInfo braStateInfo;
StateInfo halfbraStateInfo;// It has the same left and right StateInfo as braStateInfo. However, its total quanta is pb.ketquanta.
// It is used to project solution into to braStateInfo.
std::vector<Wavefunction> solution; solution.resize(1);
std::vector<Wavefunction> outputState; outputState.resize(1);
std::vector<Wavefunction> solutionprojector; solutionprojector.resize(1);
solution[0].LoadWavefunctionInfo(ketStateInfo, newSystem.get_sites(), baseState);
#ifndef SERIAL
broadcast(world, ketStateInfo, 0);
broadcast(world, solution, 0);
#endif
outputState[0].AllowQuantaFor(newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta);
outputState[0].set_onedot(solution[0].get_onedot());
outputState[0].Clear();
solutionprojector[0].AllowQuantaFor(newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.ketquanta);
solutionprojector[0].set_onedot(solution[0].get_onedot());
solutionprojector[0].Clear();
//TensorProduct (newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta[0], EqualQ, braStateInfo);
//TODO
//TensorProduct do not support const StateInfo&
TensorProduct (newSystem.set_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta[0], EqualQ, braStateInfo);
TensorProduct (newSystem.set_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.ketquanta[0], EqualQ, halfbraStateInfo);
//StateInfo::restore(forward, environmentsites, envStateInfo, baseState);
//DiagonalMatrix e;
//if(i == 0)
// GuessWave::guess_wavefunctions(solution, e, big, TRANSPOSE, true, true, 0.0, baseState);
//else
// GuessWave::guess_wavefunctions(solution, e, big, TRANSFORM, true, true, 0.0, baseState);
//SpinAdapted::operatorfunctions::Product(&newSystem, ccd, solution[0], &ketStateInfo, stateb.getw(), temp, SpinQuantum(0, SpinSpace(0), IrrepSpace(0)), true, 1.0);
boost::shared_ptr<SparseMatrix> O;
if (pb.type() == TwoPerturbType::Va)
O = newSystem.get_op_array(CDD_SUM).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
if (pb.type() == TwoPerturbType::Vi)
O = newSystem.get_op_array(CCD_SUM).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
boost::shared_ptr<SparseMatrix> overlap = newSystem.get_op_array(OVERLAP).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
SpinAdapted::operatorfunctions::TensorMultiply(*O, &braStateInfo, &ketStateInfo , solution[0], outputState[0], pb.delta, true, 1.0);
SpinAdapted::operatorfunctions::TensorMultiply(*overlap, &halfbraStateInfo, &ketStateInfo , solution[0], solutionprojector[0], overlap->get_deltaQuantum(0), true, 1.0);
DensityMatrix bratracedMatrix(newSystem.get_braStateInfo());
bratracedMatrix.allocate(newSystem.get_braStateInfo());
double norm = DotProduct(outputState[0], outputState[0]);
if(norm > NUMERICAL_ZERO)
SpinAdapted::operatorfunctions::MultiplyProduct(outputState[0], Transpose(const_cast<Wavefunction&> (outputState[0])), bratracedMatrix, 0.5/norm);
SpinAdapted::operatorfunctions::MultiplyProduct(solutionprojector[0], Transpose(const_cast<Wavefunction&> (solutionprojector[0])), bratracedMatrix, 0.5);
std::vector<Matrix> brarotateMatrix, ketrotateMatrix;
LoadRotationMatrix (newSystem.get_sites(), ketrotateMatrix, baseState);
double error;
if (!mpigetrank())
error = makeRotateMatrix(bratracedMatrix, brarotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
#ifndef SERIAL
broadcast(world, ketrotateMatrix, 0);
broadcast(world, brarotateMatrix, 0);
#endif
SaveRotationMatrix (newSystem.get_sites(), brarotateMatrix, pb.wavenumber());
newSystem.transform_operators(brarotateMatrix,ketrotateMatrix);
SpinBlock::store (forward, newSystem.get_sites(), newSystem, pb.wavenumber(), baseState); // if restart, just restoring an existing block --
system=newSystem;
}
//TODO
//.........这里部分代码省略.........
示例6: finalEnergy
//.........这里部分代码省略.........
if (dmrginp.outputlevel() > 0) {
if (forward)
{
pout << "\t\t\t Current direction is :: Forwards " << endl;
}
else
{
pout << "\t\t\t Current direction is :: Backwards " << endl;
}
}
if (sweepParams.get_block_iter() != 0)
sweepParams.set_guesstype() = TRANSFORM;
else
sweepParams.set_guesstype() = TRANSPOSE;
if (dmrginp.outputlevel() > 0)
pout << "\t\t\t Blocking and Decimating " << endl;
SpinBlock newSystem; // new system after blocking and decimating
newSystem.nonactive_orb() = pb.orb();
//Need to substitute by:
// if (warmUp )
// Startup(sweepParams, system, newSystem, dot_with_sys, pb.wavenumber(), baseState);
// else {
// BlockDecimateAndCompress (sweepParams, system, newSystem, false, dot_with_sys, pb.wavenumber(), baseState);
// }
BlockDecimateAndCompress (sweepParams, system, newSystem, warmUp, dot_with_sys,pb, baseState);
//Need to substitute by?
system = newSystem;
if (dmrginp.outputlevel() > 0){
pout << system<<endl;
pout << system.get_braStateInfo()<<endl;
system.printOperatorSummary();
}
//system size is going to be less than environment size
if (forward && system.get_complementary_sites()[0] >= dmrginp.last_site()/2)
dot_with_sys = false;
if (!forward && system.get_sites()[0]-1 < dmrginp.last_site()/2)
dot_with_sys = false;
SpinBlock::store (forward, system.get_sites(), system, pb.wavenumber(), baseState);
syssites = system.get_sites();
if (dmrginp.outputlevel() > 0)
pout << "\t\t\t saving state " << syssites.size() << endl;
++sweepParams.set_block_iter();
#ifndef SERIAL
mpi::communicator world;
world.barrier();
#endif
sweepParams.savestate(forward, syssites.size());
if (dmrginp.outputlevel() > 0)
mcheck("at the end of sweep iteration");
}
//FIXME
//It does not seem necessary.
//when we are doing twodot, we still need to do the last sweep to make sure that the
//correctionVector and base wavefunction are propogated correctly across sweeps
// //especially when we switch from twodot to onedot algorithm
// if (!sweepParams.get_onedot() && !warmUp) {
// pout << "\t\t\t Block Iteration :: " << sweepParams.get_block_iter() << endl;
// pout << "\t\t\t ----------------------------" << endl;
// if (dmrginp.outputlevel() > 0) {
// if (forward)
// pout << "\t\t\t Current direction is :: Forwards " << endl;
// else
// pout << "\t\t\t Current direction is :: Backwards " << endl;
// }
// sweepParams.set_onedot() = true;
// sweepParams.set_env_add() = 0;
// bool dot_with_sys = true;
// WavefunctionCanonicalize(sweepParams, system, warmUp, dot_with_sys, targetState, baseState);
// sweepParams.set_onedot() = false;
// sweepParams.set_env_add() = 1;
// }
//
pout << "\t\t\t Largest Error for Sweep with " << sweepParams.get_keep_states() << " states is " << finalError << endl;
pout << "\t\t\t Largest overlap for Sweep with " << sweepParams.get_keep_states() << " states is " << finalEnergy[0] << endl;
sweepParams.set_largest_dw() = finalError;
pout << "\t\t\t ============================================================================ " << endl;
// update the static number of iterations
++sweepParams.set_sweep_iter();
return finalError;
}
示例7: bratracedMatrix
//.........这里部分代码省略.........
mpi::communicator world;
broadcast(world, solution, 0);
#endif
outputState[0].AllowQuantaFor(big.get_leftBlock()->get_braStateInfo(), big.get_rightBlock()->get_braStateInfo(),pb.braquanta);
outputState[0].set_onedot(sweepParams.get_onedot());
outputState[0].Clear();
if (pb.type() == TwoPerturbType::Va)
big.multiplyCDD_sum(solution[0],&(outputState[0]),MAX_THRD);
if (pb.type() == TwoPerturbType::Vi)
big.multiplyCCD_sum(solution[0],&(outputState[0]),MAX_THRD);
//davidson_f(solution[0], outputState[0]);
SpinBlock newbig;
if (sweepParams.get_onedot() && !dot_with_sys)
{
InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, baseState, pb.wavenumber(), systemDot.size(), dmrginp.direct(), system.get_integralIndex(), DISTRIBUTED_STORAGE, false, true,NO_PARTICLE_SPIN_NUMBER_CONSTRAINT,pb.braquanta,pb.ketquanta);
InitBlocks::InitBigBlock(newSystem, environment, newbig,pb.braquanta,pb.ketquanta);
Wavefunction tempwave = outputState[0];
GuessWave::onedot_shufflesysdot(big.get_braStateInfo(), newbig.get_braStateInfo(), outputState[0], tempwave);
outputState[0] = tempwave;
tempwave = solution[0];
GuessWave::onedot_shufflesysdot(big.get_ketStateInfo(), newbig.get_ketStateInfo(), solution[0], tempwave);
solution[0] = tempwave;
big.get_rightBlock()->clear();
big.clear();
}
else
newbig = big;
DensityMatrix bratracedMatrix(newSystem.get_braStateInfo());
bratracedMatrix.allocate(newSystem.get_braStateInfo());
//bratracedMatrix.makedensitymatrix(outputState, newbig, dmrginp.weights(sweepiter), 0.0, 0.0, true);
bratracedMatrix.makedensitymatrix(outputState, newbig, std::vector<double>(1,1.0), 0.0, 0.0, true);
if (sweepParams.get_noise() > NUMERICAL_ZERO) {
pout << "adding noise "<<trace(bratracedMatrix)<<" "<<sweepiter<<" "<<dmrginp.weights(sweepiter)[0]<<endl;
bratracedMatrix.add_onedot_noise_forCompression(solution[0], newbig, sweepParams.get_noise()*max(1.0,trace(bratracedMatrix)));
if (trace(bratracedMatrix) <1e-14)
bratracedMatrix.SymmetricRandomise();
pout << "after noise "<<trace(bratracedMatrix)<<" "<<sweepParams.get_noise()<<endl;
}
environment.clear();
newEnvironment.clear();
std::vector<Matrix> brarotateMatrix, ketrotateMatrix;
LoadRotationMatrix (newSystem.get_sites(), ketrotateMatrix, baseState);
double braerror;
if (!mpigetrank()) {
braerror = makeRotateMatrix(bratracedMatrix, brarotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
}
#ifndef SERIAL
broadcast(world, ketrotateMatrix, 0);
broadcast(world, brarotateMatrix, 0);
#endif
if (dmrginp.outputlevel() > 0)
pout << "\t\t\t Total bra discarded weight "<<braerror<<endl<<endl;
sweepParams.set_lowest_error() = braerror;
SaveRotationMatrix (newbig.get_leftBlock()->get_sites(), brarotateMatrix, pb.wavenumber());
//FIXME
//It is neccessary for twodot algorithm to save baseState wavefuntion.
//I do not know why.
solution[0].SaveWavefunctionInfo (newbig.get_ketStateInfo(), newbig.get_leftBlock()->get_sites(), baseState);
outputState[0].SaveWavefunctionInfo (newbig.get_braStateInfo(), newbig.get_leftBlock()->get_sites(), pb.wavenumber());
//TODO
//Why do I need this?
//They should have been consistent.
// solution[0].SaveWavefunctionInfo (newbig.get_ketStateInfo(), newbig.get_leftBlock()->get_sites(), baseState);
// SaveRotationMatrix (newbig.get_leftBlock()->get_sites(), ketrotateMatrix, baseState);
if (dmrginp.outputlevel() > 0)
pout <<"\t\t\t Performing Renormalization "<<endl;
newSystem.transform_operators(brarotateMatrix, ketrotateMatrix);
if (dmrginp.outputlevel() > 0)
mcheck("after rotation and transformation of block");
if (dmrginp.outputlevel() > 0){
pout << *dmrginp.guessgenT<<" "<<*dmrginp.multiplierT<<" "<<*dmrginp.operrotT<< " "<<globaltimer.totalwalltime()<<" timer "<<endl;
pout << *dmrginp.makeopsT<<" makeops "<<endl;
pout << *dmrginp.datatransfer<<" datatransfer "<<endl;
pout <<"oneindexopmult twoindexopmult Hc couplingcoeff"<<endl;
pout << *dmrginp.oneelecT<<" "<<*dmrginp.twoelecT<<" "<<*dmrginp.hmultiply<<" "<<*dmrginp.couplingcoeff<<" hmult"<<endl;
pout << *dmrginp.buildsumblock<<" "<<*dmrginp.buildblockops<<" build block"<<endl;
pout << "addnoise S_0_opxop S_1_opxop S_2_opxop"<<endl;
pout << *dmrginp.addnoise<<" "<<*dmrginp.s0time<<" "<<*dmrginp.s1time<<" "<<*dmrginp.s2time<<endl;
}
}
示例8: BlockAndDecimate
void SweepGenblock::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, int state)
{
if (dmrginp.outputlevel() > 0)
mcheck("at the start of block and decimate");
// figure out if we are going forward or backwards
pout << "\t\t\t Performing Blocking"<<endl;
dmrginp.guessgenT -> start();
bool forward = (system.get_sites() [0] == 0);
SpinBlock systemDot;
int systemDotStart, systemDotEnd;
int systemDotSize = sweepParams.get_sys_add() - 1;
if (forward)
{
systemDotStart = *system.get_sites().rbegin () + 1;
systemDotEnd = systemDotStart + systemDotSize;
}
else
{
systemDotStart = system.get_sites() [0] - 1;
systemDotEnd = systemDotStart - systemDotSize;
}
vector<int> spindotsites(2);
spindotsites[0] = systemDotStart;
spindotsites[1] = systemDotEnd;
systemDot = SpinBlock(systemDotStart, systemDotEnd);
const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();
system.addAdditionalCompOps();
InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, sweepParams.get_sys_add(), dmrginp.direct(), DISTRIBUTED_STORAGE, dot_with_sys, true);
pout << "\t\t\t System Block"<<newSystem;
if (dmrginp.outputlevel() > 0)
newSystem.printOperatorSummary();
std::vector<Matrix> rotateMatrix;
if (!dmrginp.get_fullrestart()) {
//this should be done when we actually have wavefunctions stored, otherwise not!!
SpinBlock environment, environmentDot, newEnvironment;
int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;
InitBlocks::InitNewEnvironmentBlock(environment, systemDot, newEnvironment, system, systemDot,
sweepParams.get_sys_add(), sweepParams.get_env_add(), forward, dmrginp.direct(),
sweepParams.get_onedot(), nexact, useSlater, true, true, true);
SpinBlock big;
InitBlocks::InitBigBlock(newSystem, newEnvironment, big);
DiagonalMatrix e;
std::vector<Wavefunction> solution(1);
GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, true, 0.0);
solution[0].SaveWavefunctionInfo (big.get_stateInfo(), big.get_leftBlock()->get_sites(), state);
DensityMatrix tracedMatrix;
tracedMatrix.allocate(newSystem.get_stateInfo());
tracedMatrix.makedensitymatrix(solution, big, std::vector<double>(1, 1.0), 0.0, 0.0, false);
rotateMatrix.clear();
if (!mpigetrank())
double error = newSystem.makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
}
else
LoadRotationMatrix (newSystem.get_sites(), rotateMatrix, state);
#ifndef SERIAL
mpi::communicator world;
broadcast(world, rotateMatrix, 0);
#endif
if (!dmrginp.get_fullrestart())
SaveRotationMatrix (newSystem.get_sites(), rotateMatrix, state);
pout <<"\t\t\t Performing Renormalization "<<endl<<endl;
newSystem.transform_operators(rotateMatrix);
if (dmrginp.outputlevel() > 0)
mcheck("after rotation and transformation of block");
if (dmrginp.outputlevel() > 0)
pout <<newSystem<<endl;
if (dmrginp.outputlevel() > 0)
newSystem.printOperatorSummary();
//mcheck("After renorm transform");
}