本文整理汇总了C++中Wavefunction类的典型用法代码示例。如果您正苦于以下问题:C++ Wavefunction类的具体用法?C++ Wavefunction怎么用?C++ Wavefunction使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了Wavefunction类的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: reinitialize_phibetas
void SwappedSystem::reinitialize_phibetas (const Wavefunction<amplitude_t>::Amplitude &phialpha1, const Wavefunction<amplitude_t>::Amplitude &phialpha2)
{
assert(subsystem_particle_counts_match());
#if defined(DEBUG_VMC_SWAPPED_SYSTEM) || defined(DEBUG_VMC_ALL)
for (unsigned int species = 0; species < copy1_subsystem_indices.size(); ++species)
std::cerr << "swapping " << copy1_subsystem_indices[species].size() << " particles of species " << species << std::endl;
std::cerr << std::endl;
#endif
PositionArguments swapped_r1(phialpha1.get_positions()), swapped_r2(phialpha2.get_positions());
swap_positions(swapped_r1, swapped_r2);
phibeta1 = phialpha1.clone();
phibeta1->reset(swapped_r1);
phibeta1_dirty = false;
phibeta2 = phialpha2.clone();
phibeta2->reset(swapped_r2);
phibeta2_dirty = false;
#ifdef VMC_CAREFUL
verify_phibetas(phialpha1, phialpha2);
#endif
}
开发者ID:garrison,项目名称:vmc,代码行数:25,代码来源:SwappedSystem.cpp
示例2: Ctmp
//------------------------------------------------------------------------------
void BasisHarmonicOscillator::createInitalDiscretization()
{
C = zeros<cx_mat>(nGrid, nSpatialOrbitals);
mat Ctmp(nGrid, nSpatialOrbitals);
Wavefunction *wf;
for(int j=0; j<nSpatialOrbitals; j++){
wf = new HarmonicOscillator(cfg, states[2*j]);
for(int i=0; i<nGrid; i++){
Ctmp(i,j) = wf->evaluate(grid->at(i));
}
delete wf;
}
C.set_real(Ctmp);
// Forcing orthogonality by performing a SVD decomposition
SVD(C);
#ifdef DEBUG
cout << "BasisHarmonicOscillator::discretization1d()" << endl;
for(int i=0; i<nSpatialOrbitals; i++){
for(int j=0; j<nSpatialOrbitals; j++){
cout << "|C(" << i << ", " << j << ")| = " << sqrt(cdot(C.col(i),C.col(j))) << endl ;
}
}
#endif
}
开发者ID:sigvebs,项目名称:MCTDHF,代码行数:27,代码来源:basisharmonicoscillator.cpp
示例3: GetWavefunctionParticleExchange
/*
* Maps the given wavefunction to one where the particles are exchanged
* psi(1,2) -> psi(2,1)
*/
Wavefunction<3>::Ptr GetWavefunctionParticleExchange(Wavefunction<3>::Ptr psi, list angularSymmetrizationPairs)
{
typedef Array<cplx, 3> ArrayType;
ArrayType data = psi->GetData();
int countr = data.extent(1);
typedef stl_input_iterator<tuple> Iterator;
Iterator begin(angularSymmetrizationPairs);
Iterator end;
Wavefunction<3>::Ptr exchgPsi = psi->Copy();
ArrayType exchgData = exchgPsi->GetData();
for (Iterator i=begin; i!=end; i++)
{
int a1 = extract<int>((*i)[0]);
int a2 = extract<int>((*i)[1]);
for (int r1=0; r1<countr; r1++)
{
for (int r2=0; r2<countr; r2++)
{
exchgData(a1, r1, r2) = data(a2, r2, r1);
}
}
}
return exchgPsi;
}
开发者ID:AtomAleks,项目名称:pyprop-helium,代码行数:33,代码来源:analysis.cpp
示例4: verify_phibetas
void SwappedSystem::verify_phibetas (const Wavefunction<amplitude_t>::Amplitude &phialpha1, const Wavefunction<amplitude_t>::Amplitude &phialpha2) const
{
#ifdef NDEBUG
(void) phialpha1;
(void) phialpha2;
#else
const PositionArguments &r1 = phialpha1.get_positions();
const PositionArguments &r2 = phialpha2.get_positions();
assert(r1.get_N_species() == r2.get_N_species());
assert(r1.get_N_sites() == r2.get_N_sites());
assert(copy1_subsystem_indices.size() == r1.get_N_species());
assert(copy2_subsystem_indices.size() == r1.get_N_species());
const Lattice &lattice = phialpha1.get_lattice();
for (unsigned int species = 0; species < r1.get_N_species(); ++species) {
const unsigned int N = r1.get_N_filled(species);
assert(N == r2.get_N_filled(species));
// verify that the subsystem index arrays have everything they need (and no duplicates!)
unsigned int c1 = 0, c2 = 0;
for (unsigned int i = 0; i < N; ++i) {
const Particle particle(i, species);
const bool b1 = vector_find(copy1_subsystem_indices[species], i) != -1;
const bool b2 = vector_find(copy2_subsystem_indices[species], i) != -1;
if (b1)
++c1;
if (b2)
++c2;
assert(b1 == subsystem->position_is_within(r1[particle], lattice));
assert(b2 == subsystem->position_is_within(r2[particle], lattice));
}
assert(c1 == c2);
assert(c1 == copy1_subsystem_indices[species].size());
assert(c2 == copy2_subsystem_indices[species].size());
}
assert(phibeta1 != 0);
assert(phibeta2 != 0);
// verify that the positions in the phibeta's are correct
PositionArguments swapped_r1(phialpha1.get_positions()), swapped_r2(phialpha2.get_positions());
swap_positions(swapped_r1, swapped_r2);
for (unsigned int species = 0; species < r1.get_N_species(); ++species) {
for (unsigned int i = 0; i < r1.get_N_filled(species); ++i) {
const Particle particle(i, species);
assert(swapped_r1[particle] == phibeta1->get_positions()[particle]);
assert(swapped_r2[particle] == phibeta2->get_positions()[particle]);
}
}
#endif
}
开发者ID:garrison,项目名称:vmc,代码行数:55,代码来源:SwappedSystem.cpp
示例5: DotProduct
//solves the equation (H-E)|psi_1> = -QV|\Psi_0> where lowerState[0] contains |\Psi_0> to enforce orthogonality (Q), and lowerState[1] contains V|\Psi_0> so we can calculate its projection in the krylov space
//the algorithm is taken from wikipedia page
void SpinAdapted::Linear::ConjugateGradient(Wavefunction& xi, DiagonalMatrix& h_diag, double E0, double normtol, Davidson_functor& h_multiply, std::vector<Wavefunction> &lowerStates)
{
pout.precision (12);
#ifndef SERIAL
mpi::communicator world;
#endif
int iter = 0;
double levelshift = 0.0;
Wavefunction b = lowerStates[1];
//make b[0] orthogonal to lowerStates[0]
double overlap = DotProduct(b, lowerStates[0]);
double overlap2 = DotProduct(lowerStates[0], lowerStates[0]);
ScaleAdd(-overlap/overlap2, lowerStates[0], b);
//normalise all the guess roots
if(mpigetrank() == 0) {
Normalise(xi);
}
Wavefunction pi, ri;
ri=xi;
ri.Clear();
h_multiply(xi, ri);
ScaleAdd(-E0, xi, ri); // (H-E0)|psi>
ScaleAdd(-1.0, lowerStates[1], ri);
Scale(-1.0, ri);
pi = ri;
double oldError = DotProduct(ri, ri);
while(true) {
Wavefunction Hp = pi;
Hp.Clear();
h_multiply(pi, Hp);
ScaleAdd(-E0, pi, Hp); // (H-E0)|psi>
double alpha = oldError/DotProduct(pi, Hp);
ScaleAdd(alpha, pi, xi);
ScaleAdd(-alpha, Hp, ri);
double Error = DotProduct(ri, ri);
if (Error > normtol)
return;
else {
double beta = Error/oldError;
oldError = Error;
ScaleAdd(1.0/beta, ri, pi);
Scale(beta, pi);
}
}
}
开发者ID:junyang4711,项目名称:Block,代码行数:56,代码来源:linear.C
示例6: assert
void DistributedOverlapMatrix<Rank>::SetupRank(Wavefunction<Rank> &srcPsi, int opRank)
{
//Need a copy of srcPsi for future reference (only on first call to this function)
if (!HasPsi)
{
Psi = srcPsi.Copy();
HasPsi = true;
}
//Check that distribution for opRank has not changed since last call.
//Also check that typeID of representation is the same
/* int curDistribOpRank = Psi->GetRepresentation()->GetDistributedModel()->GetDistribution()(opRank);
int srcDistribOpRank = srcPsi.GetRepresentation()->GetDistributedModel()->GetDistribution()(opRank);
if ( (curDistribOpRank != srcDistribOpRank) )
{
Psi = srcPsi.Copy();
//NB: We reset IsSetup flag for _all_ ranks!
IsSetupRank = false;
}
*/
if (!IsSetupRank(opRank))
{
//Sanity check: operation rank should be less than rank of wavefunction (and nonzero, duh)
assert(opRank < Rank);
assert(opRank > -1);
//Assert non-orthogonal rank opRank
assert (!srcPsi.GetRepresentation()->IsOrthogonalBasis(opRank));
//Create Epetra map for this rank
WavefunctionMaps(opRank) = CreateWavefunctionMultiVectorEpetraMap<Rank>(Psi, opRank);
//Setup overlap matrix
SetupOverlapMatrixRank(srcPsi, opRank);
//Setup work multivectors
blitz::Array<cplx, 3> psiData = MapToRank3(srcPsi.GetData(), opRank, 1);
int numVectors = 2 * psiData.extent(0) * psiData.extent(2);
InputVector(opRank) = Epetra_MultiVector_Ptr( new Epetra_MultiVector(*WavefunctionMaps(opRank), numVectors, false) );
OutputVector(opRank) = Epetra_MultiVector_Ptr( new Epetra_MultiVector(*WavefunctionMaps(opRank), numVectors, false) );
//Allocate mem for multivectors
//InputData.resize(srcPsi.GetData().size(), 1);
//OutputData.resize(srcPsi.GetData().size(), 1);
//Setup Amesos solvers
SetupOverlapSolvers(srcPsi, opRank);
//Flag this rank as set up
IsSetupRank(opRank) = true;
}
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:53,代码来源:distributedoverlapmatrix.cpp
示例7: IsFermion
void SpinAdapted::operatorfunctions::TensorMultiply(const SpinBlock *ablock, const Baseoperator<Matrix>& a, const Baseoperator<Matrix>& b, const SpinBlock *cblock, Wavefunction& c, Wavefunction& v, const SpinQuantum opQ, double scale)
{
// can be used for situation with different bra and ket
const int leftBraOpSz = cblock->get_leftBlock()->get_braStateInfo().quanta.size ();
const int leftKetOpSz = cblock->get_leftBlock()->get_ketStateInfo().quanta.size ();
const int rightBraOpSz = cblock->get_rightBlock()->get_braStateInfo().quanta.size ();
const int rightKetOpSz = cblock->get_rightBlock()->get_ketStateInfo().quanta.size ();
const StateInfo* lbraS = cblock->get_braStateInfo().leftStateInfo, *rbraS = cblock->get_braStateInfo().rightStateInfo;
const StateInfo* lketS = cblock->get_ketStateInfo().leftStateInfo, *rketS = cblock->get_ketStateInfo().rightStateInfo;
const char conjC = (cblock->get_leftBlock() == ablock) ? 'n' : 't';
const Baseoperator<Matrix>& leftOp = (conjC == 'n') ? a : b; // an ugly hack to support the release memory optimisation
const Baseoperator<Matrix>& rightOp = (conjC == 'n') ? b : a;
const char leftConj = (conjC == 'n') ? a.conjugacy() : b.conjugacy();
const char rightConj = (conjC == 'n') ? b.conjugacy() : a.conjugacy();
int totalmem =0;
for (int lQrQPrime = 0; lQrQPrime<leftBraOpSz*rightKetOpSz; ++lQrQPrime)
{
int rQPrime = lQrQPrime%rightKetOpSz, lQ = lQrQPrime/rightKetOpSz;
for (int lQPrime = 0; lQPrime < leftKetOpSz; lQPrime++)
if (leftOp.allowed(lQ, lQPrime) && c.allowed(lQPrime, rQPrime))
{
Matrix m; m.ReSize(lbraS->getquantastates(lQ), rketS->getquantastates(rQPrime));
double factor = leftOp.get_scaling(lbraS->quanta[lQ], lketS->quanta[lQPrime]);
MatrixMultiply (leftOp.operator_element(lQ, lQPrime), leftConj, c.operator_element(lQPrime, rQPrime), 'n',
m, factor, 0.);
for (int rQ = 0; rQ<rightBraOpSz; rQ++) {
if (v.allowed(lQ, rQ) && rightOp.allowed(rQ, rQPrime)) {
double factor = scale;
factor *= dmrginp.get_ninej()(lketS->quanta[lQPrime].get_s().getirrep(), rketS->quanta[rQPrime].get_s().getirrep() , c.get_deltaQuantum(0).get_s().getirrep(),
leftOp.get_spin().getirrep(), rightOp.get_spin().getirrep(), opQ.get_s().getirrep(),
lbraS->quanta[lQ].get_s().getirrep(), rbraS->quanta[rQ].get_s().getirrep() , v.get_deltaQuantum(0).get_s().getirrep());
factor *= Symmetry::spatial_ninej(lketS->quanta[lQPrime].get_symm().getirrep() , rketS->quanta[rQPrime].get_symm().getirrep(), c.get_symm().getirrep(),
leftOp.get_symm().getirrep(), rightOp.get_symm().getirrep(), opQ.get_symm().getirrep(),
lbraS->quanta[lQ].get_symm().getirrep() , rbraS->quanta[rQ].get_symm().getirrep(), v.get_symm().getirrep());
int parity = rightOp.get_fermion() && IsFermion(lketS->quanta[lQPrime]) ? -1 : 1;
factor *= rightOp.get_scaling(rbraS->quanta[rQ], rketS->quanta[rQPrime]);
MatrixMultiply (m, 'n', rightOp(rQ, rQPrime), TransposeOf(rightOp.conjugacy()), v.operator_element(lQ, rQ), factor*parity);
}
}
}
}
}
开发者ID:chrinide,项目名称:Block,代码行数:53,代码来源:operatorfunctions.C
示例8: real
void DistributedOverlapMatrix<Rank>::WavefunctionToMultiVector(Wavefunction<Rank> &psi, Epetra_MultiVector_Ptr vec, int opRank)
{
//Map wavefunction to 3D array (compress before- and after-ranks)
blitz::Array<cplx, 3> psiData = MapToRank3(psi.GetData(), opRank, 1);
int beforeSize = psiData.extent(0);
int opSize = psiData.extent(1);
int afterSize = psiData.extent(2);
int otherSize = beforeSize*afterSize;
//Copy real and imag part of wavefunction into multivector
double realVal, imagVal;
for (int i=0; i<beforeSize; i++)
{
for (int j=0; j<opSize; j++)
{
for (int k=0; k<afterSize; k++)
{
//cout << "cur vec = " << i * afterSize + k << ", cur idx = " << j << endl;
//cout << "cur vec = " << i * afterSize + k + otherSize << ", cur idx = " << j << endl;
realVal = real( psiData(i,j,k) );
imagVal = imag( psiData(i,j,k) );
vec->ReplaceMyValue(j, i*afterSize + k, realVal);
vec->ReplaceMyValue(j, i*afterSize + k + otherSize, imagVal);
}
}
}
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:28,代码来源:distributedoverlapmatrix.cpp
示例9: SetupRank
shared_ptr<Epetra_MultiVector> DistributedOverlapMatrix<Rank>::SetupMultivector(Wavefunction<Rank> &srcPsi, int opRank)
{
//Setup Epetra stuff for opRank
SetupRank(srcPsi, opRank);
//Map wavefunction to 3D array (compress before- and after-ranks)
blitz::Array<cplx, 3> psiData = MapToRank3(srcPsi.GetData(), opRank, 1);
int beforeSize = psiData.extent(0);
int opSize = psiData.extent(1);
int afterSize = psiData.extent(2);
int otherSize = beforeSize*afterSize;
//Copy real and imag part of wavefunction into 2D array
blitz::Array<double, 2> data;
data.resize(2 * otherSize, opSize);
for (int i=0; i<beforeSize; i++)
{
for (int j=0; j<opSize; j++)
{
for (int k=0; k<afterSize; k++)
{
data(i*afterSize + k, j) = real( psiData(i,j,k) );
data(i*afterSize + k + otherSize, j) = imag( psiData(i,j,k) );
}
}
}
//Number of vectors (in multivector)
int numVectors = 2 * otherSize;
//Create Epetra multivector (view of data)
Epetra_MultiVector_Ptr srcVec = Epetra_MultiVector_Ptr(new Epetra_MultiVector(Copy, *WavefunctionMaps(opRank), data.data(), opSize, numVectors));
return srcVec;
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:35,代码来源:distributedoverlapmatrix.cpp
示例10: shape
void DistributedOverlapMatrix<Rank>::MultiVectorToWavefunction(Wavefunction<Rank> &psi, Epetra_MultiVector_Ptr vec, int opRank)
{
//Put result into psi
blitz::Array<cplx, 3> data = MapToRank3(psi.GetData(), opRank, 1);
int beforeSize = data.extent(0);
int opSize = data.extent(1);
int afterSize = data.extent(2);
int otherSize = beforeSize*afterSize;
double *destVecView=0;
vec->ExtractView(&destVecView, &opSize);
blitz::Array<double, 2> data2(destVecView, shape(2 * otherSize, opSize), shape(opSize, 1), neverDeleteData);
//blitz::Array<double, 2> data2;
//data2.resize(2 * otherSize, opSize);
//vec->ExtractCopy(data2.data(), opSize);
for (int i=0; i<beforeSize; i++)
{
for (int j=0; j<opSize; j++)
{
for (int k=0; k<afterSize; k++)
{
double realVal = data2(i*afterSize + k, j);
double imagVal = data2(i*afterSize + k + otherSize, j);
data(i,j,k) = realVal + I * imagVal;
}
}
}
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:27,代码来源:distributedoverlapmatrix.cpp
示例11:
void CombinedRepresentation<Rank>::MultiplyOverlap(Wavefunction<Rank> &srcPsi, Wavefunction<Rank> &dstPsi, int rank)
{
using namespace blitz;
if (this->IsOrthogonalBasis(rank))
{
dstPsi.GetData() = srcPsi.GetData();
}
else
{
DistributedOverlap->MultiplyOverlapRank(srcPsi, dstPsi, rank, true);
//Map the data to a 3D array, where the b-spline part is the middle rank
//Array<cplx, 3> srcData = MapToRank3(srcPsi.Data, rank, 1);
//Array<cplx, 3> dstData = MapToRank3(dstPsi.Data, rank, 1);
//this->GetGlobalOverlapMatrix(rank)->MultiplyOverlapTensor(srcData, dstData);
}
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:19,代码来源:combinedrepresentation.cpp
示例12:
void BSplineTransform<Rank>::SetupStep(Wavefunction<Rank> &psi, BSpline::Ptr bsplineObject, int baseRank)
{
using namespace blitz;
SetBaseRank(baseRank);
BSplineObject = bsplineObject;
/*
* Get shape of wavefunction (we are in the bspline repr).
* Then allocate bspline grid representation buffer of wavefunction
* and store buffer names on object.
*/
TinyVector<int, Rank> gridShape = psi.GetData().shape();
gridShape(baseRank) = BSplineObject->GetQuadratureGridGlobal().extent(0);
BSplineGridDataName = psi.AllocateData(gridShape);
BSplineDataName = psi.GetActiveBufferName();
//TempData.resize(BSplineObject->GetQuadratureGridGlobal().extent(0));
TempData.resize(BSplineObject->NumberOfBSplines);
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:20,代码来源:bsplinetransform.cpp
示例13: count_subsystem_particle_counts_for_match
bool count_subsystem_particle_counts_for_match (const Wavefunction<amplitude_t>::Amplitude &wf1, const Wavefunction<amplitude_t>::Amplitude &wf2,
const Subsystem &subsystem)
{
assert(subsystem.lattice_makes_sense(wf1.get_lattice()));
assert(subsystem.lattice_makes_sense(wf2.get_lattice()));
// (we are also assuming that the lattices are in fact equivalent)
const PositionArguments &r1 = wf1.get_positions();
const PositionArguments &r2 = wf2.get_positions();
assert(r1.get_N_species() == r2.get_N_species());
assert(r1.get_N_sites() == r2.get_N_sites());
for (unsigned int species = 0; species < r1.get_N_species(); ++species) {
assert(r1.get_N_filled(species) == r2.get_N_filled(species));
unsigned int count1 = 0, count2 = 0;
for (unsigned int i = 0; i < r1.get_N_filled(species); ++i) {
const Particle particle(i, species);
if (subsystem.position_is_within(r1[particle], wf1.get_lattice()))
++count1;
if (subsystem.position_is_within(r2[particle], wf2.get_lattice()))
++count2;
}
if (count1 != count2)
return false;
}
return true;
}
开发者ID:garrison,项目名称:vmc,代码行数:32,代码来源:SwappedSystem.cpp
示例14: runtime_error
void CombinedRepresentation<Rank>::MultiplyOverlap(cplx sourceScaling, Wavefunction<Rank> &srcPsi, cplx destScaling, Wavefunction<Rank> &dstPsi, int rank)
{
//Not implemented correctly atm.
throw std::runtime_error("Not implemented properly yet!");
using namespace blitz;
if (this->IsOrthogonalBasis(rank))
{
dstPsi.GetData() = srcPsi.GetData();
}
else
{
//Map the data to a 3D array, where the b-spline part is the middle rank
//Array<cplx, 3> srcData = MapToRank3(srcPsi.Data, rank, 1);
//Array<cplx, 3> dstData = MapToRank3(dstPsi.Data, rank, 1);
//this->GetGlobalOverlapMatrix(rank)->MultiplyOverlapTensor(srcData, dstData);
DistributedOverlap->MultiplyOverlapRank(srcPsi, dstPsi, rank, true);
}
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:22,代码来源:combinedrepresentation.cpp
示例15: runtime_error
void BSplineTransform<Rank>::InverseTransform(Wavefunction<Rank> &psi)
{
using namespace blitz;
if(psi.GetActiveBufferName() != BSplineDataName)
{
throw std::runtime_error("Active databuffer is not what is should be...");
}
//Project input and output into 3D arrays with bsplineRank in the middle.
Array<cplx, Rank> inputData(psi.GetData());
Array<cplx, Rank> outputData(psi.GetData(BSplineGridDataName));
Array<cplx, 3> input3d = MapToRank3(inputData, BaseRank, 1);
Array<cplx, 3> output3d = MapToRank3(outputData, BaseRank, 1);
int psiSliceExtent = input3d.extent(1);
Array<cplx, 1> psiSlice(psiSliceExtent);
output3d = 0;
int preCount = input3d.extent(0);
int postCount = input3d.extent(2);
for (int i=0; i<preCount; i++)
{
for (int j=0; j<postCount; j++)
{
/*
* Copy wavefunction slice along bspline rank to temp
* array.
*/
psiSlice = input3d(i, Range::all(), j).copy();
// Call on BSpline function to perform expansion
output3d(i, Range::all() , j) =
BSplineObject->ConstructFunctionFromBSplineExpansion(psiSlice);
}
}
psi.SetActiveBuffer(BSplineGridDataName);
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:38,代码来源:bsplinetransform.cpp
示例16: initialize
void SwappedSystem::initialize (const Wavefunction<amplitude_t>::Amplitude &phialpha1, const Wavefunction<amplitude_t>::Amplitude &phialpha2)
{
assert(current_state == UNINITIALIZED);
const PositionArguments &r1 = phialpha1.get_positions();
const PositionArguments &r2 = phialpha2.get_positions();
// FIXME: we need a way to assert that phialpha1 and phialpha2 represent
// the same wave function, just with different amplitudes. Then again,
// we're only calling this function from two places in the code where this
// can be easily verified ...
#ifndef NDEBUG
assert(r1.get_N_species() == r2.get_N_species());
for (unsigned int i = 0; i < r1.get_N_species(); ++i)
assert(r1.get_N_filled(i) == r2.get_N_filled(i));
#endif
assert(r1.get_N_sites() == r2.get_N_sites());
assert(subsystem->lattice_makes_sense(phialpha1.get_lattice()));
assert(subsystem->lattice_makes_sense(phialpha2.get_lattice()));
assert(&phialpha1.get_lattice() == &phialpha2.get_lattice());
const unsigned int N_species = r1.get_N_species();
copy1_subsystem_indices.resize(N_species);
copy2_subsystem_indices.resize(N_species);
for (unsigned int species = 0; species < N_species; ++species) {
const unsigned int N = r1.get_N_filled(species);
for (unsigned int i = 0; i < N; ++i) {
const Particle particle(i, species);
if (subsystem->position_is_within(r1[particle], phialpha1.get_lattice()))
copy1_subsystem_indices[species].push_back(i);
if (subsystem->position_is_within(r2[particle], phialpha2.get_lattice()))
copy2_subsystem_indices[species].push_back(i);
}
}
assert(subsystem_particle_counts_match());
reinitialize_phibetas(phialpha1, phialpha2);
current_state = READY;
}
开发者ID:garrison,项目名称:vmc,代码行数:41,代码来源:SwappedSystem.cpp
示例17: weights
void CombinedRepresentation<Rank>::MultiplyIntegrationWeights(Wavefunction<Rank> &psi, int rank)
{
blitz::Array<cplx, 3> data = MapToRank3(psi.GetData(), rank, 1);
if (this->IsOrthogonalBasis(rank))
{
blitz::Array<double, 1> weights = this->GetLocalWeights(rank);
data *= weights(blitz::tensor::j) + 0*blitz::tensor::k;
}
else
{
if (this->GetDistributedModel()->IsDistributedRank(rank))
{
DistributedOverlap->MultiplyOverlapRank(psi, rank, true);
}
else
{
this->GetGlobalOverlapMatrix(rank)->MultiplyOverlapTensor(data);
}
}
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:20,代码来源:combinedrepresentation.cpp
示例18: evExp
void Propagator<Rank>::Setup(const Parameter ¶m, const cplx &dt, const Wavefunction<Rank> &psi, int rank)
{
firstIndex i;
secondIndex j;
thirdIndex k;
//Set class parameters
N = psi.GetRepresentation()->GetFullShape()(rank);
PropagateRank = rank;
Param = param;
//create some temporary arrays
Array<double, 2> evInvReal; //inverse eigenvector matrix
//Call setup routines to get the eigenvector decomposition
setup(N, param, Eigenvalues, Eigenvectors, evInvReal);
//We need the eigenvalues and vectors in complex format
Array<cplx, 2> evExp(Eigenvectors.shape()); //eigenvectors scaled by exp(eigenvalues)
Array<cplx, 2> evDiff(Eigenvectors.shape()); //eigenvectors scaled by eigenvalues
Array<cplx, 2> evInv(Eigenvectors.shape()); //inverse eigenvector matrix
evInv = evInvReal(tensor::i, tensor::j);
//scale eigenvectors by complex rotation
//The missing minus sign in the exponent is included in the matrix.
evExp = Eigenvectors(i,j) * exp( I * dt * Eigenvalues(j) / (2.0 * Mass));
evDiff = - Eigenvectors(i,j) * Eigenvalues(j) / (2.0 * Mass);
//Create full matrix to propagate wavefunction
PropagationMatrix.resize(Eigenvectors.shape());
MatrixMatrixMultiply(evExp, evInv, PropagationMatrix);
//Create full differentiation matrix
DiffMatrix.resize(Eigenvectors.shape());
MatrixMatrixMultiply(evDiff, evInv, DiffMatrix);
//Allocate temp data
TempData.resize(Eigenvalues.extent(0));
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:40,代码来源:transformedgridpropagator.cpp
示例19: runtime_error
void DistributedOverlapMatrix<Rank>::SetupOverlapSolvers(Wavefunction<Rank> &psi, int opRank)
{
//Set up Epetra LinearProblem with overlap for this rank and input/output multivectors
EpetraProblems(opRank) = Epetra_LinearProblem_Ptr( new Epetra_LinearProblem(OverlapMatrices(opRank).get(), OutputVector(opRank).get(), InputVector(opRank).get()) );
//Determine solver type. Use SuperLU_dist if opRank is distributed, else KLU
Amesos Factory;
std::string SolverType;
if (psi.GetRepresentation()->GetDistributedModel()->IsDistributedRank(opRank))
{
SolverType = "Amesos_Superludist";
}
else
{
SolverType = "Amesos_Klu";
}
//Create Amesos solver for this rank
Solvers(opRank) = Amesos_BaseSolver_Ptr( Factory.Create(SolverType, *EpetraProblems(opRank)) );
//Check that requested solver exists in Amesos build
if (Solvers(opRank) == 0)
{
throw std::runtime_error("Specified Amesos solver not available");
}
//Setup the parameter list for the solver (TODO: get these from Python config section)
Teuchos::ParameterList List;
List.set("MatrixType", "symmetric");
List.set("PrintTiming", false);
List.set("PrintStatus", false);
List.set("ComputeTrueResidual", false);
Solvers(opRank)->SetParameters(List);
//Perform symbolic factorization
Solvers(opRank)->SymbolicFactorization();
Solvers(opRank)->NumericFactorization();
}
开发者ID:AtomAleks,项目名称:PyProp,代码行数:38,代码来源:distributedoverlapmatrix.cpp
示例20: diagonal
void SpinAdapted::Linear::precondition(Wavefunction& op, double e, DiagonalMatrix& diagonal, double levelshift)
{
if (!mpigetrank())
{
int index = 1;
for (int lQ = 0; lQ < op.nrows (); ++lQ)
for (int rQ = 0; rQ < op.ncols (); ++rQ)
if (op.allowed(lQ, rQ))
for (int lQState = 0; lQState < op.operator_element(lQ, rQ).Nrows (); ++lQState)
for (int rQState = 0; rQState < op.operator_element(lQ, rQ).Ncols (); ++rQState)
{
if (fabs(e - diagonal(index)) > 1.e-12)
op.operator_element(lQ, rQ).element(lQState, rQState) /= (e - diagonal(index)+levelshift);
++index;
}
}
#ifndef SERIAL
//mpi::communicator world;
//mpi::broadcast(world, op, 0);
#endif
}
开发者ID:junyang4711,项目名称:Block,代码行数:21,代码来源:linear.C
注:本文中的Wavefunction类示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论