本文整理汇总了C++中Baseoperator类的典型用法代码示例。如果您正苦于以下问题：C++ Baseoperator类的具体用法？C++ Baseoperator怎么用？C++ Baseoperator使用的例子？那么, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了Baseoperator类的9个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: 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);
}
}
}
}
}示例2: u
/*
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();
Wavefunction u;
u.resize(leftBraOpSz*leftKetOpSz, rightKetOpSz);
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))
{
int lindex = lQ*leftKetOpSz+lQPrime;
u.allowed(lindex, rQPrime) = true;
u(lindex,rQPrime).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',
u.operator_element(lindex, rQPrime), factor, 0.);
}
}
}
pout << "after first step in tensormultiply"<<endl;
mcheck("before davidson but after all blocks are built");
{
for (int lQrQ = 0; lQrQ<leftBraOpSz*rightBraOpSz; ++lQrQ)
{
int rQ = lQrQ%rightBraOpSz, lQ=lQrQ/rightBraOpSz;
if (v.allowed(lQ, rQ))
for (int rQPrime = 0; rQPrime < rightKetOpSz; rQPrime++)
if (rightOp.allowed(rQ, rQPrime))
for (int lQPrime = 0; lQPrime < leftKetOpSz; lQPrime++)
if (leftOp.allowed(lQ, lQPrime) && u.allowed(lQ*leftKetOpSz+lQPrime, rQPrime))
{
int lindex = lQ*leftKetOpSz+lQPrime;
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 (u.operator_element(lindex, rQPrime), 'n',
rightOp(rQ, rQPrime), TransposeOf(rightOp.conjugacy()), v.operator_element(lQ, rQ), factor*parity);
}
}
}
}
*/
void SpinAdapted::operatorfunctions::OperatorScaleAdd(double scaleV, const SpinBlock& b, const Baseoperator<Matrix>& op1, Baseoperator<Matrix>& op2)
{
const StateInfo& s = b.get_stateInfo();
for (int lQ = 0; lQ< op2.nrows(); lQ++)
for (int rQ = 0; rQ<op2.ncols(); rQ++)
if (op2.allowed(lQ, rQ) && op1.allowed(lQ,rQ))
{
double factor = op1.get_scaling(s.quanta[lQ], s.quanta[rQ]);
if (op1.conjugacy() == 't')
MatrixScaleAdd(scaleV*factor, op1.operator_element(lQ,rQ).t(), op2.operator_element(lQ,rQ));
else
MatrixScaleAdd(scaleV*factor, op1.operator_element(lQ,rQ), op2.operator_element(lQ,rQ));
}
}示例3: assert
void SpinAdapted::operatorfunctions::TensorTrace (const SpinBlock *ablock, const Baseoperator<Matrix>& a, const SpinBlock* cblock, const StateInfo* cstateinfo, DiagonalMatrix& cDiagonal, Real scale)
{
if (fabs(scale) < TINY) return;
try
{
assert (a.get_initialised());
const char conjC = (ablock == cblock->get_leftBlock()) ? 'n' : 't';
const int aSz = ablock->get_stateInfo().quanta.size ();
const int bSz = (conjC == 'n') ? cblock->get_stateInfo().rightStateInfo->quanta.size () : cblock->get_stateInfo().leftStateInfo->quanta.size ();
const StateInfo& s = cblock->get_stateInfo();
const StateInfo* lS = s.leftStateInfo, *rS = s.rightStateInfo;
for (int aQ = 0; aQ < aSz; ++aQ)
if (a.allowed(aQ, aQ))
for (int bQ = 0; bQ < bSz; ++bQ)
if (s.allowedQuanta (aQ, bQ, conjC))
{
int cQ = s.quantaMap (aQ, bQ, conjC)[0];
for (int cQState = 0; cQState < s.quantaStates [cQ]; ++cQState)
{
Real scaleB = 1.0;
int aQState;
int bQState;
if (conjC == 'n')
{
s.UnMapQuantumState (cQState, s.rightStateInfo->quantaStates [bQ], aQState, bQState);
scaleB *= dmrginp.get_ninej()(lS->quanta[aQ].get_s().getirrep() , rS->quanta[bQ].get_s().getirrep(), cstateinfo->quanta[cQ].get_s().getirrep(),
a.get_spin().getirrep(), 0, 0,
lS->quanta[aQ].get_s().getirrep() , rS->quanta[bQ].get_s().getirrep(), cstateinfo->quanta[cQ].get_s().getirrep());
scaleB *= Symmetry::spatial_ninej(lS->quanta[aQ].get_symm().getirrep() , rS->quanta[bQ].get_symm().getirrep(), cstateinfo->quanta[cQ].get_symm().getirrep(),
a.get_symm().getirrep(), 0, 0,
lS->quanta[aQ].get_symm().getirrep() , rS->quanta[bQ].get_symm().getirrep(), cstateinfo->quanta[cQ].get_symm().getirrep());
}
else
{
scaleB *= dmrginp.get_ninej()(lS->quanta[bQ].get_s().getirrep() , rS->quanta[aQ].get_s().getirrep(), cstateinfo->quanta[cQ].get_s().getirrep(),
0, a.get_spin().getirrep(), 0,
lS->quanta[bQ].get_s().getirrep() , rS->quanta[aQ].get_s().getirrep(), cstateinfo->quanta[cQ].get_s().getirrep());
scaleB *= Symmetry::spatial_ninej(lS->quanta[bQ].get_symm().getirrep() , rS->quanta[aQ].get_symm().getirrep(), cstateinfo->quanta[cQ].get_symm().getirrep(),
0, a.get_symm().getirrep(), 0,
lS->quanta[bQ].get_symm().getirrep() , rS->quanta[aQ].get_symm().getirrep(), cstateinfo->quanta[cQ].get_symm().getirrep());
if (a.get_fermion()&& IsFermion(lS->quanta[bQ])) scaleB *= -1.0;
s.UnMapQuantumState (cQState, s.rightStateInfo->quantaStates [aQ], bQState, aQState);
}
cDiagonal(s.unBlockedIndex [cQ] + cQState + 1) += scale * scaleB * a.operator_element(aQ, aQ)(aQState + 1, aQState + 1);
}
}
}
catch (Exception)
{
pout << Exception::what () << endl;
abort ();
}
}示例4: default
void SpinAdapted::operatorfunctions::TensorProduct (const SpinBlock *ablock, const Baseoperator<Matrix>& a, const Baseoperator<Matrix>& b, const SpinBlock *cblock, const StateInfo *cstateinfo, Baseoperator<Matrix>& c, double scale, int num_thrds)
{
if (fabs(scale) < TINY) return;
int rows = c.nrows();
#ifdef _OPENMP
#pragma omp parallel default(shared) num_threads(num_thrds)
#endif
{
#ifdef _OPENMP
#pragma omp for schedule(dynamic)
#endif
for (int cq = 0; cq < rows; ++cq)
for (int cqprime = 0; cqprime < rows; ++cqprime)
if (c.allowed(cq, cqprime)) {
TensorProductElement(ablock, a, b, cblock, cstateinfo, c, c.operator_element(cq, cqprime), cq, cqprime, scale);
}
}
}示例5: MatrixMultiply
void SpinAdapted::operatorfunctions::TensorMultiply(const Baseoperator<Matrix>& a, const StateInfo *brastateinfo, const StateInfo *ketstateinfo, const Wavefunction& c, Wavefunction& v, const SpinQuantum dQ, bool left, double scale, int num_thrds)
{
//Calculate O_{l or r} |\Psi> without building big block.
const StateInfo* lbraS = brastateinfo->leftStateInfo, *lketS = ketstateinfo->leftStateInfo;
const StateInfo* rbraS = brastateinfo->rightStateInfo, *rketS = ketstateinfo->rightStateInfo;
const int leftBraOpSz = brastateinfo->leftStateInfo->quanta.size ();
const int leftKetOpSz = ketstateinfo->leftStateInfo->quanta.size ();
const int rightBraOpSz = brastateinfo->rightStateInfo->quanta.size ();
const int rightKetOpSz = ketstateinfo->rightStateInfo->quanta.size ();
if (left)
{
//#pragma omp parallel default(shared) num_threads(num_thrds)
{
//#pragma omp for schedule(dynamic)
for (int lQ = 0; lQ < leftBraOpSz; ++lQ) {
for (int lQPrime = 0; lQPrime < leftKetOpSz; ++lQPrime)
{
if (a.allowed(lQ, lQPrime))
{
const Matrix& aop = a.operator_element(lQ, lQPrime);
for (int rQ = 0; rQ < rightKetOpSz; ++rQ)
if (c.allowed(lQPrime, rQ) && v.allowed(lQ, rQ))
{
double fac=scale;
fac *= dmrginp.get_ninej()(lketS->quanta[lQPrime].get_s().getirrep(), rketS->quanta[rQ].get_s().getirrep() , c.get_deltaQuantum(0).get_s().getirrep(),
a.get_spin().getirrep(), 0, a.get_spin().getirrep(),
lbraS->quanta[lQ].get_s().getirrep(), rbraS->quanta[rQ].get_s().getirrep() , v.get_deltaQuantum(0).get_s().getirrep());
fac *= Symmetry::spatial_ninej(lketS->quanta[lQPrime].get_symm().getirrep() , rketS->quanta[rQ].get_symm().getirrep(), c.get_symm().getirrep(),
a.get_symm().getirrep(), 0, a.get_symm().getirrep(),
lbraS->quanta[lQ].get_symm().getirrep() , rbraS->quanta[rQ].get_symm().getirrep(), v.get_symm().getirrep());
fac *= a.get_scaling(lbraS->quanta[lQ], lketS->quanta[lQPrime]);
MatrixMultiply (aop, a.conjugacy(), c.operator_element(lQPrime, rQ), c.conjugacy(),
v.operator_element(lQ, rQ), fac);
}
}
}
}
}
}
else
{
//#pragma omp parallel default(shared) num_threads(num_thrds)
{
//#pragma omp for schedule(dynamic)
for (int rQ = 0; rQ < rightBraOpSz; ++rQ) {
for (int rQPrime = 0; rQPrime < rightKetOpSz; ++rQPrime)
if (a.allowed(rQ, rQPrime))
{
const Matrix& aop = a.operator_element(rQ, rQPrime);
for (int lQPrime = 0; lQPrime < leftKetOpSz; ++lQPrime)
if (v.allowed(lQPrime, rQ) && c.allowed(lQPrime, rQPrime)) {
double fac = scale;
fac *= dmrginp.get_ninej()(lketS->quanta[lQPrime].get_s().getirrep(), rketS->quanta[rQPrime].get_s().getirrep() , c.get_deltaQuantum(0).get_s().getirrep(),
0, a.get_spin().getirrep(), a.get_spin().getirrep(),
lbraS->quanta[lQPrime].get_s().getirrep(), rbraS->quanta[rQ].get_s().getirrep() , v.get_deltaQuantum(0).get_s().getirrep());
fac *= Symmetry::spatial_ninej(lketS->quanta[lQPrime].get_symm().getirrep() , rketS->quanta[rQPrime].get_symm().getirrep(), c.get_symm().getirrep(),
0, a.get_symm().getirrep(), a.get_symm().getirrep(),
lbraS->quanta[lQPrime].get_symm().getirrep() , rbraS->quanta[rQ].get_symm().getirrep(), v.get_symm().getirrep());
fac *= a.get_scaling(rbraS->quanta[rQ], rketS->quanta[rQPrime]);
double parity = a.get_fermion() && IsFermion(lketS->quanta[lQPrime]) ? -1 : 1;
MatrixMultiply (c.operator_element(lQPrime, rQPrime), c.conjugacy(),
aop, TransposeOf(a.conjugacy()), v.operator_element(lQPrime, rQ), fac*parity);
}
}
}
}
}
}示例6: MatrixDiagonalScale
void SpinAdapted::operatorfunctions::TensorProduct (const SpinBlock *ablock, const Baseoperator<Matrix>& a, const Baseoperator<Matrix>& b, const SpinBlock* cblock, const StateInfo* cstateinfo, DiagonalMatrix& cDiagonal, double scale)
{
if (fabs(scale) < TINY) return;
const int aSz = a.nrows();
const int bSz = b.nrows();
const char conjC = (cblock->get_leftBlock() == ablock) ? 'n' : 't';
const SpinBlock* bblock = (cblock->get_leftBlock() == ablock) ? cblock->get_rightBlock() : cblock->get_leftBlock();
const StateInfo& s = cblock->get_stateInfo();
const StateInfo* lS = s.leftStateInfo, *rS = s.rightStateInfo;
for (int aQ = 0; aQ < aSz; ++aQ)
if (a.allowed(aQ, aQ))
for (int bQ = 0; bQ < bSz; ++bQ)
if (b.allowed(bQ, bQ))
if (s.allowedQuanta (aQ, bQ, conjC))
{
int cQ = s.quantaMap (aQ, bQ, conjC)[0];
Real scaleA = scale;
Real scaleB = 1;
if (conjC == 'n')
{
scaleB *= dmrginp.get_ninej()(lS->quanta[aQ].get_s().getirrep() , rS->quanta[bQ].get_s().getirrep(), cstateinfo->quanta[cQ].get_s().getirrep(),
a.get_spin().getirrep(), b.get_spin().getirrep(), 0,
lS->quanta[aQ].get_s().getirrep() , rS->quanta[bQ].get_s().getirrep(), cstateinfo->quanta[cQ].get_s().getirrep());
scaleB *= Symmetry::spatial_ninej(lS->quanta[aQ].get_symm().getirrep() , rS->quanta[bQ].get_symm().getirrep(), cstateinfo->quanta[cQ].get_symm().getirrep(),
a.get_symm().getirrep(), b.get_symm().getirrep(), 0,
lS->quanta[aQ].get_symm().getirrep() , rS->quanta[bQ].get_symm().getirrep(), cstateinfo->quanta[cQ].get_symm().getirrep());
if (b.get_fermion() && IsFermion (lS->quanta [aQ])) scaleB *= -1.0;
for (int aQState = 0; aQState < lS->quantaStates[aQ] ; aQState++)
MatrixDiagonalScale(a.operator_element(aQ, aQ)(aQState+1, aQState+1)*scaleA*scaleB, b.operator_element(bQ, bQ),
cDiagonal.Store()+s.unBlockedIndex[cQ]+aQState*rS->quantaStates[bQ]);
}
else
{
scaleB *= dmrginp.get_ninej()(lS->quanta[bQ].get_s().getirrep() , rS->quanta[aQ].get_s().getirrep(), cstateinfo->quanta[cQ].get_s().getirrep(),
b.get_spin().getirrep(), a.get_spin().getirrep(), 0,
lS->quanta[bQ].get_s().getirrep() , rS->quanta[aQ].get_s().getirrep(), cstateinfo->quanta[cQ].get_s().getirrep());
scaleB *= Symmetry::spatial_ninej(lS->quanta[bQ].get_symm().getirrep() , rS->quanta[aQ].get_symm().getirrep(), cstateinfo->quanta[cQ].get_symm().getirrep(),
b.get_symm().getirrep(), a.get_symm().getirrep(), 0,
lS->quanta[bQ].get_symm().getirrep() , rS->quanta[aQ].get_symm().getirrep(), cstateinfo->quanta[cQ].get_symm().getirrep());
if (a.get_fermion()&& IsFermion(lS->quanta[bQ])) scaleB *= -1.0;
for (int bQState = 0; bQState < lS->quantaStates[bQ] ; bQState++)
MatrixDiagonalScale(b.operator_element(bQ, bQ)(bQState+1, bQState+1)*scaleA*scaleB, a.operator_element(aQ, aQ),
cDiagonal.Store()+s.unBlockedIndex[cQ]+bQState*rS->quantaStates[aQ]);
}
}
}示例7: unitMatrix
void SpinAdapted::operatorfunctions::TensorTraceElement(const SpinBlock *ablock, const Baseoperator<Matrix>& a, const SpinBlock *cblock, const StateInfo *cstateinfo, Baseoperator<Matrix>& c, Matrix& cel, int cq, int cqprime, double scale)
{
if (fabs(scale) < TINY) return;
assert(c.allowed(cq, cqprime));
int aq, aqprime, bq, bqprime, bstates;
const char conjC = (ablock == cblock->get_leftBlock()) ? 'n' : 't';
const std::vector<int> oldToNewI = cstateinfo->oldToNewState.at(cq);
const std::vector<int> oldToNewJ = cstateinfo->oldToNewState.at(cqprime);
const StateInfo* rS = cstateinfo->rightStateInfo, *lS = cstateinfo->leftStateInfo;
int rowstride =0, colstride = 0;
for (int oldi =0; oldi < oldToNewI.size(); oldi++) {
colstride = 0;
for (int oldj = 0; oldj < oldToNewJ.size(); oldj++)
{
if (conjC == 'n')
{
aq = cstateinfo->leftUnMapQuanta[ oldToNewI[oldi] ];
aqprime = cstateinfo->leftUnMapQuanta[ oldToNewJ[oldj] ];
bq = cstateinfo->rightUnMapQuanta[ oldToNewI[oldi] ];
bqprime = cstateinfo->rightUnMapQuanta[ oldToNewJ[oldj] ];
bstates = cstateinfo->rightStateInfo->getquantastates(bq);
}
else
{
aq = cstateinfo->rightUnMapQuanta[ oldToNewI[oldi] ];
aqprime = cstateinfo->rightUnMapQuanta[ oldToNewJ[oldj] ];
bq = cstateinfo->leftUnMapQuanta[ oldToNewI[oldi] ];
bqprime = cstateinfo->leftUnMapQuanta[ oldToNewJ[oldj] ];
bstates = cstateinfo->leftStateInfo->getquantastates(bq);
}
if (a.allowed(aq, aqprime) && (bq == bqprime))
{
DiagonalMatrix unitMatrix(bstates);
unitMatrix = 1.;
Matrix unity(bstates, bstates);
unity = unitMatrix;
if (conjC == 'n')
{
double scaleb = dmrginp.get_ninej()(lS->quanta[aqprime].get_s().getirrep() , rS->quanta[bqprime].get_s().getirrep(), cstateinfo->quanta[cqprime].get_s().getirrep(),
a.get_spin().getirrep(), 0, c.get_spin().getirrep(),
lS->quanta[aq].get_s().getirrep() , rS->quanta[bq].get_s().getirrep(), cstateinfo->quanta[cq].get_s().getirrep());
scaleb *= Symmetry::spatial_ninej(lS->quanta[aqprime].get_symm().getirrep() , rS->quanta[bqprime].get_symm().getirrep(), cstateinfo->quanta[cqprime].get_symm().getirrep(),
a.get_symm().getirrep(), 0, c.get_symm().getirrep(),
lS->quanta[aq].get_symm().getirrep() , rS->quanta[bq].get_symm().getirrep(), cstateinfo->quanta[cq].get_symm().getirrep());
MatrixTensorProduct (a.operator_element(aq, aqprime), a.conjugacy(), scale, unity, 'n', scaleb,
cel, rowstride, colstride);
}
else {
double scaleb = dmrginp.get_ninej()(lS->quanta[bqprime].get_s().getirrep(), rS->quanta[aqprime].get_s().getirrep() , cstateinfo->quanta[cqprime].get_s().getirrep(),
0, a.get_spin().getirrep(), c.get_spin().getirrep(),
lS->quanta[bq].get_s().getirrep(), rS->quanta[aq].get_s().getirrep() , cstateinfo->quanta[cq].get_s().getirrep());
scaleb *= Symmetry::spatial_ninej(lS->quanta[bqprime].get_symm().getirrep() , rS->quanta[aqprime].get_symm().getirrep(), cstateinfo->quanta[cqprime].get_symm().getirrep(),
0, a.get_symm().getirrep(), c.get_symm().getirrep(),
lS->quanta[bq].get_symm().getirrep() , rS->quanta[aq].get_symm().getirrep(), cstateinfo->quanta[cq].get_symm().getirrep());
if (a.get_fermion() && IsFermion (cstateinfo->leftStateInfo->quanta[bqprime]) ) scaleb *= -1.;
MatrixTensorProduct (unity, 'n', scaleb, a.operator_element(aq, aqprime), a.conjugacy(), scale,
cel, rowstride, colstride);
}
}
colstride += cstateinfo->unCollectedStateInfo->quantaStates[ oldToNewJ[oldj] ];
}
rowstride += cstateinfo->unCollectedStateInfo->quantaStates[ oldToNewI[oldi] ];
}
}示例8: racah
void SpinAdapted::operatorfunctions::Product (const SpinBlock *ablock, const Baseoperator<Matrix>& a, const Baseoperator<Matrix>& b, Baseoperator<Matrix>& c, double scale)
{
const StateInfo* astate = &ablock->get_stateInfo();
if (fabs(scale) < TINY) return;
int rows = c.nrows();
for (int cq = 0; cq < rows; ++cq)
for (int cqprime = 0; cqprime < rows; ++cqprime)
if (c.allowed(cq, cqprime))
for (int aprime = 0; aprime < rows; aprime++)
if (a.allowed(cq, aprime) && b.allowed(aprime, cqprime))
{
int apj = astate->quanta[aprime].get_s().getirrep(), cqj = astate->quanta[cq].get_s().getirrep(), cqpj = astate->quanta[cqprime].get_s().getirrep();
double factor = a.get_scaling(astate->quanta[cq], astate->quanta[aprime]);
factor *= b.get_scaling(astate->quanta[aprime], astate->quanta[cqprime]);
if(dmrginp.spinAdapted()){
factor *= racah(cqpj, b.get_spin().getirrep(), cqj, a.get_spin().getirrep(), apj, c.get_spin().getirrep()) * pow( (1.0*c.get_spin().getirrep()+1.0)*(1.0*apj+1.0), 0.5 )
*pow(-1.0, static_cast<int>((b.get_spin().getirrep()+a.get_spin().getirrep()-c.get_spin().getirrep())/2.0));
}
MatrixMultiply(a.operator_element(cq, aprime), a.conjugacy(), b.operator_element(aprime, cqprime), b.conjugacy(),
c.operator_element(cq, cqprime), scale*factor, 1.0);
}
}示例9: assert
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)
{
const int leftOpSz = cblock->get_leftBlock()->get_stateInfo().quanta.size ();
const int rightOpSz = cblock->get_rightBlock()->get_stateInfo().quanta.size ();
const StateInfo* rS = cblock->get_stateInfo().rightStateInfo, *lS = cblock->get_stateInfo().leftStateInfo;
assert (cblock->get_leftBlock() == ablock || cblock->get_rightBlock() == ablock);
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();
Wavefunction u;
u.resize(leftOpSz*leftOpSz, rightOpSz);
int totalmem =0;
{
for (int lQrQPrime = 0; lQrQPrime<leftOpSz*rightOpSz; ++lQrQPrime)
{
int rQPrime = lQrQPrime%rightOpSz, lQ = lQrQPrime/rightOpSz;
for (int lQPrime = 0; lQPrime < leftOpSz; lQPrime++)
if (leftOp.allowed(lQ, lQPrime) && c.allowed(lQPrime, rQPrime))
{
int lindex = lQ*leftOpSz+lQPrime;
u.allowed(lindex, rQPrime) = true;
u(lindex,rQPrime).ReSize(lS->getquantastates(lQ), rS->getquantastates(rQPrime));
double factor = leftOp.get_scaling(lS->quanta[lQ], lS->quanta[lQPrime]);
MatrixMultiply (leftOp.operator_element(lQ, lQPrime), leftConj, c.operator_element(lQPrime, rQPrime), 'n',
u.operator_element(lindex, rQPrime), factor, 0.);
}
}
}
{
for (int lQrQ = 0; lQrQ<leftOpSz*rightOpSz; ++lQrQ)
{
int rQ = lQrQ%rightOpSz, lQ=lQrQ/rightOpSz;
if (v.allowed(lQ, rQ))
for (int rQPrime = 0; rQPrime < rightOpSz; rQPrime++)
if (rightOp.allowed(rQ, rQPrime))
for (int lQPrime = 0; lQPrime < leftOpSz; lQPrime++)
if (leftOp.allowed(lQ, lQPrime) && u.allowed(lQ*leftOpSz+lQPrime, rQPrime))
{
int lindex = lQ*leftOpSz+lQPrime;
double factor = scale;
factor *= dmrginp.get_ninej()(lS->quanta[lQPrime].get_s(), rS->quanta[rQPrime].get_s() , c.get_deltaQuantum().get_s(),
leftOp.get_spin(), rightOp.get_spin(), opQ.get_s(),
lS->quanta[lQ].get_s(), rS->quanta[rQ].get_s() , v.get_deltaQuantum().get_s());
factor *= Symmetry::spatial_ninej(lS->quanta[lQPrime].get_symm().getirrep() , rS->quanta[rQPrime].get_symm().getirrep(), c.get_symm().getirrep(),
leftOp.get_symm().getirrep(), rightOp.get_symm().getirrep(), opQ.get_symm().getirrep(),
lS->quanta[lQ].get_symm().getirrep() , rS->quanta[rQ].get_symm().getirrep(), v.get_symm().getirrep());
int parity = rightOp.get_fermion() && IsFermion(lS->quanta[lQPrime]) ? -1 : 1;
factor *= rightOp.get_scaling(rS->quanta[rQ], rS->quanta[rQPrime]);
MatrixMultiply (u.operator_element(lindex, rQPrime), 'n',
rightOp(rQ, rQPrime), TransposeOf(rightOp.conjugacy()), v.operator_element(lQ, rQ), factor*parity);
}
}
}
}