39 for ( proshade_unsign bwIt = 0; bwIt < this->
maxShellBand; bwIt++ )
45 for ( proshade_unsign shIt = 0; shIt < this->
noSpheres; shIt++ )
67 this->allocateRRPMemory ( settings );
70 proshade_double descValR = 0.0;
71 proshade_unsign arrPos1, arrPos2;
72 for ( proshade_unsign band = 0; band < this->maxShellBand; band++ )
75 for ( proshade_unsign shell1 = 0; shell1 < this->noSpheres; shell1++ )
80 for ( proshade_unsign shell2 = 0; shell2 < this->noSpheres; shell2++ )
82 this->rrpMatrices[band][shell1][shell2] = 0.0;
83 this->rrpMatrices[band][shell2][shell1] = 0.0;
88 for ( proshade_unsign shell2 = 0; shell2 < this->noSpheres; shell2++ )
91 if ( shell1 > shell2 ) {
continue; }
96 this->rrpMatrices[band][shell1][shell2] = 0.0;
97 this->rrpMatrices[band][shell2][shell1] = 0.0;
105 for ( proshade_unsign order = 0; order < static_cast<proshade_unsign> ( ( 2 * band ) + 1 ); order++ )
107 arrPos1 =
static_cast<proshade_unsign
> ( seanindex (
static_cast<proshade_signed
> ( order ) -
108 static_cast<proshade_signed
> ( band ),
109 band, this->spheres[shell1]->getLocalBandwidth() ) );
110 arrPos2 =
static_cast<proshade_unsign
> ( seanindex (
static_cast<proshade_signed
> ( order ) -
111 static_cast<proshade_signed
> ( band ),
112 band, this->spheres[shell2]->getLocalBandwidth() ) );
114 &this->sphericalHarmonics[shell1][arrPos1][1],
115 &this->sphericalHarmonics[shell2][arrPos2][0],
116 &this->sphericalHarmonics[shell2][arrPos2][1] );
120 this->rrpMatrices[band][shell1][shell2] = descValR;
121 this->rrpMatrices[band][shell2][shell1] = descValR;
145 if ( bandInQuestion < spheres[shellInQuestion]->getLocalBandwidth() )
171 proshade_double ret = 0.0;
172 std::vector<proshade_double> bandDists;
177 throw ProSHADE_exception (
"Attempted computing energy levels descriptors when it was not required.",
"ED00017", __FILE__, __LINE__, __func__,
"Attempted to pre-compute the RRP matrices, when the user\n : has specifically stated that these should not be computed.\n : Unless you manipulated the code, this error should never\n : occur; if you see this, I made a large blunder. Please let\n : me know!" );
192 ret =
static_cast<proshade_double
> ( std::accumulate ( bandDists.begin(), bandDists.end(), 0.0 ) ) /
193 static_cast<proshade_double
> ( bandDists.size() );
221 proshade_double *str1Vals =
new proshade_double[minCommonShells * minCommonShells];
222 proshade_double *str2Vals =
new proshade_double[minCommonShells * minCommonShells];
225 proshade_unsign arrIter = 0;
228 for ( proshade_unsign band = 0; band < minCommonBands; band++ )
234 for ( proshade_unsign shell1 = 0; shell1 < minCommonShells; shell1++ )
239 for ( proshade_unsign shell2 = 0; shell2 < minCommonShells; shell2++ )
245 str1Vals[arrIter] = obj1->
getRRPValue ( band, shell1, shell2 ) *
248 str2Vals[arrIter] = obj2->
getRRPValue ( band, shell1, shell2 ) *
283 this->maxCompBand = band;
286 this->eMatrices =
new proshade_complex** [this->maxCompBand];
289 for ( proshade_unsign bandIter = 0; bandIter < this->maxCompBand; bandIter++ )
292 this->eMatrices[bandIter] =
new proshade_complex* [
static_cast<proshade_unsign
> ( ( bandIter * 2 ) + 1 )];
295 for ( proshade_unsign band2Iter = 0; band2Iter < static_cast<proshade_unsign> ( ( bandIter * 2 ) + 1 ); band2Iter++ )
297 this->eMatrices[bandIter][band2Iter] =
new proshade_complex [
static_cast<proshade_unsign
> ( ( bandIter * 2 ) + 1 )];
320 obj1Vals =
new proshade_double [minSpheres];
321 obj2Vals =
new proshade_double [minSpheres];
322 radiiVals =
new proshade_complex[minSpheres];
323 GLabscissas =
new proshade_double [intOrder];
324 GLweights =
new proshade_double [intOrder];
355 *result *= pow (
static_cast<proshade_double
> ( obj->
getAnySphereRadius( radius ) ), 2.0 );
378 proshade_unsign objCombValsIter = 0;
379 proshade_double hlpReal, hlpImag;
380 proshade_complex arrVal;
383 for ( proshade_unsign order2Iter = 0; order2Iter < ( ( bandIter * 2 ) + 1 ); order2Iter++ )
400 &hlpReal, &hlpImag );
403 radiiVals[objCombValsIter][0] = hlpReal * pow ( (
static_cast<proshade_double
> ( obj1->
getAnySphereRadius( radiusIter ) ) ), 2.0 );
404 radiiVals[objCombValsIter][1] = hlpImag * pow ( (
static_cast<proshade_double
> ( obj1->
getAnySphereRadius( radiusIter ) ) ), 2.0 );
406 objCombValsIter += 1;
441 proshade_unsign obj1ValsIter = 0;
442 proshade_unsign obj2ValsIter = 0;
446 proshade_unsign maxSphere = 0;
453 minSphere = std::min ( radiusIter, minSphere );
454 maxSphere = std::max ( radiusIter, maxSphere );
464 proshade_double minSphereRad = obj1->
getSpherePosValue ( minSphere ) - ( sphereDist * 0.5 );
465 proshade_double maxSphereRad = obj1->
getSpherePosValue ( maxSphere ) + ( sphereDist * 0.5 );
471 return ( maxSphereRad - minSphereRad );
489 delete[] GLabscissas;
524 proshade_double *obj1Vals, *obj2Vals, *GLAbscissas, *GLWeights;
525 proshade_complex* radiiVals;
526 proshade_double integRange;
535 for ( proshade_unsign bandIter = 0; bandIter < std::min ( obj1->
getMaxBand(), obj2->
getMaxBand() ); bandIter++ )
538 for ( proshade_unsign orderIter = 0; orderIter < ( ( bandIter * 2 ) + 1 ); orderIter++ )
550 std::stringstream hlpSS;
551 hlpSS <<
"E matrices computed for band " << bandIter;
586 if ( settings->
task == Symmetry ) { eMatNormFactor /= 2.0; }
588 for ( proshade_unsign bandIter = 0; bandIter < std::min ( obj1->
getMaxBand(), obj2->
getMaxBand() ); bandIter++ )
591 for ( proshade_unsign orderIter = 0; orderIter < ( ( bandIter * 2 ) + 1 ); orderIter++ )
593 for ( proshade_unsign order2Iter = 0; order2Iter < ( ( bandIter * 2 ) + 1 ); order2Iter++ )
627 proshade_double ret = 0.0;
632 throw ProSHADE_exception (
"Attempted computing trace sigma descriptors when it was\n : not required.",
"ED00018", __FILE__, __LINE__, __func__,
"Attempted to pre-compute the E matrices, when the user\n : has specifically stated that these should not be computed.\n : Unless you manipulated the code, this error should never\n : occur; if you see this, I made a large blunder. Please let\n : me know!" );
646 double* singularValues =
new double[( ( std::min ( obj1->
getMaxBand(), obj2->
getMaxBand() ) * 2 ) + 1 )];
650 for ( proshade_unsign lIter = 0; lIter < std::min ( obj1->
getMaxBand(), obj2->
getMaxBand() ); lIter++ )
656 for ( proshade_unsign iter = 0; iter < ( ( lIter * 2 ) + 1 ); iter++ )
658 ret += singularValues[iter];
666 delete[] singularValues;
683 this->so3Coeffs =
new fftw_complex [
static_cast<proshade_unsign
>( ( 4 * pow(
static_cast<proshade_double
> ( band ), 3.0 ) -
static_cast<proshade_double
> ( band ) ) / 3.0 )];
684 this->so3CoeffsInverse =
new fftw_complex [
static_cast<proshade_unsign
>( pow(
static_cast<proshade_double
> ( band ) * 2.0, 3.0 ) )];
715 proshade_double wigNorm, hlpValReal, hlpValImag;
716 proshade_double signValue = 1.0;
717 proshade_unsign indexO;
718 proshade_complex hlpVal;
721 for ( proshade_signed bandIter = 0; bandIter < static_cast<proshade_signed> ( std::min ( obj1->
getMaxBand(), obj2->
getMaxBand() ) ); bandIter++ )
724 wigNorm = 2.0 * M_PI * sqrt ( 2.0 / (2.0 * bandIter + 1.0 ) ) ;
727 for ( proshade_signed orderIter = 0; orderIter < ( ( bandIter * 2 ) + 1 ); orderIter++ )
730 if ( ( orderIter - bandIter + bandIter ) % 2 ) { signValue = -1.0 ; }
731 else { signValue = 1.0 ; }
734 for ( proshade_signed order2Iter = 0; order2Iter < ( ( bandIter * 2 ) + 1 ); order2Iter++ )
737 indexO =
static_cast<proshade_unsign
> ( so3CoefLoc ( orderIter - bandIter, order2Iter - bandIter, bandIter, std::min ( obj1->
getMaxBand(), obj2->
getMaxBand() ) ) );
740 obj2->
getEMatrixValue ( bandIter, orderIter, order2Iter, &hlpValReal, &hlpValImag );
741 hlpVal[0] = hlpValReal * wigNorm * signValue;
742 hlpVal[1] = hlpValImag * wigNorm * signValue;
769 work1 =
new proshade_complex[8 *
static_cast<proshade_unsign
> ( pow(
static_cast<double> ( band ), 3.0 ) )];
770 work2 =
new proshade_complex[14 *
static_cast<proshade_unsign
> ( pow(
static_cast<double> ( band ), 2.0 ) ) + (48 * band)];
771 work3 =
new proshade_double [2 *
static_cast<proshade_unsign
> ( pow(
static_cast<double> ( band ), 2.0 ) ) + (24 * band)];
793 int howmany = 4 * band * band;
794 int idist = 2 * band;
795 int odist = 2 * band;
798 int inembed[2], onembed[2];
799 inembed[0] = 2 * band;
800 inembed[1] = 4 * band * band;
801 onembed[0] = 2 * band;
802 onembed[1] = 4 * band * band;
812 *inverseSO3 = fftw_plan_many_dft ( rank,
865 proshade_complex *workspace1, *workspace2;
866 proshade_double *workspace3;
867 fftw_plan inverseSO3;
887 fftw_destroy_plan ( inverseSO3 );
917 proshade_double ret = 0.0;
918 proshade_double eulA, eulB, eulG, EMatR, EMatI, WigDR, WigDI;
923 throw ProSHADE_exception (
"Attempted computing rotation function descriptors when it\n : was not required.",
"ED00023", __FILE__, __LINE__, __func__,
"Attempted to compute the SO(3) transform and the rotation \n : function descriptor when the user did not request this. \n : Unless you manipulated the code, this error should never \n : occur; if you see this, I made a large blunder. \n : Please let me know!" );
942 &eulA, &eulB, &eulG, settings );
948 for ( proshade_unsign bandIter = 0; bandIter < obj2->
getComparisonBand(); bandIter++ )
951 for ( proshade_unsign order1 = 0; order1 < ( ( bandIter * 2 ) + 1 ); order1++ )
954 for ( proshade_unsign order2 = 0; order2 < ( ( bandIter * 2 ) + 1 ); order2++ )