28 #if defined ( __GNUC__ )
29 #pragma GCC diagnostic push
30 #pragma GCC diagnostic ignored "-Wpedantic"
31 #pragma GCC diagnostic ignored "-Wshadow"
32 #pragma GCC diagnostic ignored "-Wall"
33 #pragma GCC diagnostic ignored "-Wextra"
34 #pragma GCC diagnostic ignored "-Wdouble-promotion"
35 #pragma GCC diagnostic ignored "-Wconversion"
39 #if defined ( __clang__ )
40 #pragma clang diagnostic push
41 #pragma clang diagnostic ignored "-Wpedantic"
42 #pragma clang diagnostic ignored "-Wshadow"
43 #pragma clang diagnostic ignored "-Wall"
44 #pragma clang diagnostic ignored "-Wextra"
45 #pragma clang diagnostic ignored "-Wdouble-promotion"
46 #pragma clang diagnostic ignored "-Weverything"
50 #if defined ( _MSC_VER )
51 #pragma warning ( disable:4996 )
55 #define GEMMI_WRITE_IMPLEMENTATION
56 #include <gemmi/to_pdb.hpp>
59 #if defined ( _MSC_VER )
60 #pragma warning ( default:4996 )
64 #if defined ( __GNUC__ )
65 #pragma GCC diagnostic pop
69 #if defined ( __clang__ )
70 #pragma clang diagnostic pop
76 proshade_signed
addAxisUnlessSame ( proshade_unsign fold, proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axHeight, proshade_double averageFSC, std::vector< proshade_double* >* prosp, proshade_double axErr );
81 bool checkElementAlreadyExists ( std::vector<std::vector< proshade_double > >* elements, std::vector< proshade_double >* elem, proshade_double matrixTolerance );
82 bool checkElementsFormGroup ( std::vector<std::vector< proshade_double > >* elements, proshade_double matrixTolerance );
98 this->
fileType = ProSHADE_internal_io::UNKNOWN;
163 this->
spherePos = std::vector<proshade_single> ( );
214 ProSHADE_internal_data::ProSHADE_data::ProSHADE_data ( std::string strName,
double *mapVals,
int len, proshade_single xDmSz, proshade_single yDmSz, proshade_single zDmSz, proshade_unsign xDmInd, proshade_unsign yDmInd, proshade_unsign zDmInd, proshade_signed xFr, proshade_signed yFr, proshade_signed zFr, proshade_signed xT, proshade_signed yT, proshade_signed zT, proshade_unsign inputO )
218 this->fileName = strName;
219 this->fileType = ProSHADE_internal_io::MAP;
222 this->internalMap =
nullptr;
225 this->xDimSize = xDmSz;
226 this->yDimSize = yDmSz;
227 this->zDimSize = zDmSz;
231 this->xDimIndices = xDmInd;
232 this->yDimIndices = yDmInd;
233 this->zDimIndices = zDmInd;
234 this->xGridIndices = xDmInd;
235 this->yGridIndices = yDmInd;
236 this->zGridIndices = zDmInd;
237 this->xAxisOrder = 1;
238 this->yAxisOrder = 2;
239 this->zAxisOrder = 3;
240 this->xAxisOrigin = xFr;
241 this->yAxisOrigin = yFr;
242 this->zAxisOrigin = zFr;
248 this->xDimSizeOriginal = 0.0;
249 this->yDimSizeOriginal = 0.0;
250 this->zDimSizeOriginal = 0.0;
251 this->xDimIndicesOriginal = 0;
252 this->yDimIndicesOriginal = 0;
253 this->zDimIndicesOriginal = 0;
254 this->xAxisOriginOriginal = 0;
255 this->yAxisOriginOriginal = 0;
256 this->zAxisOriginOriginal = 0;
257 this->originalMapXCom = 0.0;
258 this->originalMapYCom = 0.0;
259 this->originalMapZCom = 0.0;
260 this->mapMovFromsChangeX = 0.0;
261 this->mapMovFromsChangeY = 0.0;
262 this->mapMovFromsChangeZ = 0.0;
263 this->mapCOMProcessChangeX = 0.0;
264 this->mapCOMProcessChangeY = 0.0;
265 this->mapCOMProcessChangeZ = 0.0;
268 this->originalPdbRotCenX = 0.0;
269 this->originalPdbRotCenY = 0.0;
270 this->originalPdbRotCenZ = 0.0;
271 this->originalPdbTransX = 0.0;
272 this->originalPdbTransY = 0.0;
273 this->originalPdbTransZ = 0.0;
284 this->spherePos = std::vector<proshade_single> ( );
286 this->spheres =
nullptr;
287 this->sphericalHarmonics =
nullptr;
288 this->rotSphericalHarmonics =
nullptr;
289 this->maxShellBand = 0;
292 this->rrpMatrices =
nullptr;
293 this->eMatrices =
nullptr;
294 this->so3Coeffs =
nullptr;
295 this->so3CoeffsInverse =
nullptr;
296 this->wignerMatrices =
nullptr;
297 this->integrationWeight = 0.0;
298 this->maxCompBand = 0;
299 this->translationMap =
nullptr;
302 this->isEmpty =
false;
303 this->inputOrder = inputO;
306 if (
static_cast<proshade_unsign
> ( len ) != ( xDmInd * yDmInd * zDmInd ) )
308 throw ProSHADE_exception (
"Structure class input map has wrong dimensions.",
"EP00044", __FILE__, __LINE__, __func__,
"The supplied map array size has different dimensions to\n : the required map dimensions." );
311 if ( (
static_cast<proshade_signed
> ( xT - xFr ) !=
static_cast<proshade_signed
> ( xDmInd - 1 ) ) ||
312 (
static_cast<proshade_signed
> ( yT - yFr ) !=
static_cast<proshade_signed
> ( yDmInd - 1 ) ) ||
313 (
static_cast<proshade_signed
> ( zT - zFr ) !=
static_cast<proshade_signed
> ( zDmInd - 1 ) ) )
315 throw ProSHADE_exception (
"Structure class input dimensions not in line with map\n : to/from indices.",
"EP00045", __FILE__, __LINE__, __func__,
"The supplied map information does not add up. The\n : dimensions are not in line with the indexing start/stop\n : position distances and therefore proper map indexing\n : cannot be done. Please check the input values." );
319 this->internalMap =
new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
323 proshade_unsign arrPos = 0;
324 for ( proshade_unsign xIt = 0; xIt < this->xDimIndices; xIt++ )
326 for ( proshade_unsign yIt = 0; yIt < this->yDimIndices; yIt++ )
328 for ( proshade_unsign zIt = 0; zIt < this->zDimIndices; zIt++ )
330 arrPos = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
331 this->internalMap[arrPos] =
static_cast<proshade_double
> ( mapVals[arrPos] );
352 if ( this->internalMap !=
nullptr )
354 delete[] this->internalMap;
358 if ( this->spheres !=
nullptr )
360 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
362 if ( this->spheres[iter] !=
nullptr )
364 delete this->spheres[iter];
365 this->spheres[iter] =
nullptr;
368 delete[] this->spheres;
372 if ( this->sphericalHarmonics !=
nullptr )
374 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
376 if ( this->sphericalHarmonics[iter] !=
nullptr )
378 delete[] this->sphericalHarmonics[iter];
379 this->sphericalHarmonics[iter] =
nullptr;
382 delete[] this->sphericalHarmonics;
386 if ( this->rotSphericalHarmonics !=
nullptr )
388 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
390 if ( this->rotSphericalHarmonics[iter] !=
nullptr )
392 delete[] this->rotSphericalHarmonics[iter];
393 this->rotSphericalHarmonics[iter] =
nullptr;
396 delete[] this->rotSphericalHarmonics;
400 if ( this->rrpMatrices !=
nullptr )
402 for ( proshade_unsign bwIt = 0; bwIt < this->maxShellBand; bwIt++ )
404 if ( this->rrpMatrices[bwIt] !=
nullptr )
406 for ( proshade_unsign shIt = 0; shIt < this->noSpheres; shIt++ )
408 if ( this->rrpMatrices[bwIt][shIt] !=
nullptr )
410 delete[] this->rrpMatrices[bwIt][shIt];
414 delete[] this->rrpMatrices[bwIt];
418 delete[] this->rrpMatrices;
422 if ( this->eMatrices !=
nullptr )
424 for ( proshade_unsign bandIter = 0; bandIter < this->maxCompBand; bandIter++ )
426 if ( this->eMatrices[bandIter] !=
nullptr )
428 for ( proshade_unsign band2Iter = 0; band2Iter < static_cast<proshade_unsign> ( ( bandIter * 2 ) + 1 ); band2Iter++ )
430 if ( this->eMatrices[bandIter][band2Iter] !=
nullptr )
432 delete[] this->eMatrices[bandIter][band2Iter];
436 delete[] this->eMatrices[bandIter];
440 delete[] this->eMatrices;
444 if ( this->so3Coeffs !=
nullptr )
446 delete[] this->so3Coeffs;
448 if ( this->so3CoeffsInverse !=
nullptr )
450 delete[] this->so3CoeffsInverse;
454 if ( this->wignerMatrices !=
nullptr )
456 for ( proshade_unsign bandIter = 1; bandIter < this->maxCompBand; bandIter++ )
458 if ( this->wignerMatrices[bandIter] !=
nullptr )
460 for ( proshade_unsign order1Iter = 0; order1Iter < ( (bandIter * 2) + 1 ); order1Iter++ )
462 if ( this->wignerMatrices[bandIter][order1Iter] !=
nullptr )
464 delete[] this->wignerMatrices[bandIter][order1Iter];
467 delete[] this->wignerMatrices[bandIter];
470 delete[] wignerMatrices;
474 if ( this->translationMap !=
nullptr )
476 delete[] this->translationMap;
480 if ( this->sphereMappedRotFun.size() > 0 )
482 for ( proshade_unsign spIt = 0; spIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); spIt++ )
484 delete this->sphereMappedRotFun.at(spIt);
509 void ProSHADE_internal_data::ProSHADE_data::readInStructure ( std::string fName, proshade_unsign inputO,
ProSHADE_settings* settings, proshade_double* maskArr, proshade_unsign maskXDim, proshade_unsign maskYDim, proshade_unsign maskZDim, proshade_double* weightsArr, proshade_unsign weigXDim, proshade_unsign weigYDim, proshade_unsign weigZDim )
515 if ( !this->isEmpty )
517 throw ProSHADE_exception (
"Structure data class not empty.",
"E000005", __FILE__, __LINE__, __func__,
"Attempted to read in structure into a ProSHADE_data\n : object which already does have structure read in\n : i.e. " + this->fileName );
521 this->fileName = fName;
527 this->inputOrder = inputO;
530 switch ( this->fileType )
532 case ProSHADE_internal_io::UNKNOWN:
533 throw ProSHADE_exception (
"Unknown file type.",
"E000006", __FILE__, __LINE__, __func__,
"When attempting to read the file\n : " + this->fileName +
"\n : the file extension was determined as unknown. This could\n : mean either that the file does not exist, or that it is\n : not one of the supported extensions." );
535 case ProSHADE_internal_io::GEMMI:
536 throw ProSHADE_exception (
"Unknown file type.",
"E000006", __FILE__, __LINE__, __func__,
"When attempting to read the file\n : " + this->fileName +
"\n : the file extension was determined as unknown. This could\n : mean either that the file does not exist, or that it is\n : not one of the supported extensions." );
538 case ProSHADE_internal_io::PDB:
539 this->readInPDB ( settings );
542 case ProSHADE_internal_io::MAP:
543 this->readInMAP ( settings, maskArr, maskXDim, maskYDim, maskZDim, weightsArr, weigXDim, weigYDim, weigZDim );
548 this->isEmpty =
false;
570 std::stringstream ss;
571 ss <<
"Starting to load the structure from Gemmi object " << inputO;
575 if ( !this->isEmpty )
577 throw ProSHADE_exception (
"Structure data class not empty.",
"E000005", __FILE__, __LINE__, __func__,
"Attempted to read in structure into a ProSHADE_data\n : object which already does have structure read in\n : i.e. " + this->fileName );
581 this->fileName = gemmiStruct.name;
584 this->fileType = ProSHADE_internal_io::GEMMI;
587 this->inputOrder = inputO;
590 this->readInGemmi ( gemmiStruct, settings );
593 this->isEmpty =
false;
622 gemmi::Ccp4<float> map;
623 map.read_ccp4 ( gemmi::MaybeGzipped ( this->fileName.c_str() ) );
626 map.setup ( gemmi::GridSetup::ReorderOnly, NAN );
630 &this->xDimIndices, &this->yDimIndices, &this->zDimIndices,
631 &this->xDimSize, &this->yDimSize, &this->zDimSize,
632 &this->aAngle, &this->bAngle, &this->cAngle,
633 &this->xFrom, &this->yFrom, &this->zFrom,
634 &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin,
635 &this->xAxisOrder, &this->yAxisOrder, &this->zAxisOrder,
636 &this->xGridIndices, &this->yGridIndices, &this->zGridIndices );
639 ProSHADE_internal_io::readInMapData ( &map, this->internalMap, this->xDimIndices, this->yDimIndices, this->zDimIndices, this->xAxisOrder, this->yAxisOrder, this->zAxisOrder );
643 maskArr, maskXDim, maskYDim, maskZDim );
647 weightsArr, weigXDim, weigYDim, weigZDim );
650 if ( settings->
removeNegativeDensity ) {
for (
size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) {
if ( this->internalMap[iter] < 0.0 ) { this->internalMap[iter] = 0.0; } } }
655 settings->
setResolution ( std::min (
static_cast<proshade_single
> ( this->xDimSize ) /
static_cast<proshade_single
> ( this->xDimIndices ),
656 std::min (
static_cast<proshade_single
> ( this->yDimSize ) /
static_cast<proshade_single
> ( this->yDimIndices ),
657 static_cast<proshade_single
> ( this->zDimSize ) /
static_cast<proshade_single
> ( this->zDimIndices ) ) ) * 2.0f );
661 this->figureIndexStartStop ( );
667 proshade_single xSampRate = this->xDimSize /
static_cast< proshade_single
> ( this->xTo - this->xFrom );
668 proshade_single ySampRate = this->yDimSize /
static_cast< proshade_single
> ( this->yTo - this->yFrom );
669 proshade_single zSampRate = this->zDimSize /
static_cast< proshade_single
> ( this->zTo - this->zFrom );
672 proshade_single xStartPosBefore =
static_cast< proshade_single
> ( this->xFrom ) * xSampRate;
673 proshade_single yStartPosBefore =
static_cast< proshade_single
> ( this->yFrom ) * ySampRate;
674 proshade_single zStartPosBefore =
static_cast< proshade_single
> ( this->zFrom ) * zSampRate;
677 proshade_double xMapCOMPreReSampl = 0.0, yMapCOMPreReSampl = 0.0, zMapCOMPreReSampl = 0.0;
678 ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xMapCOMPreReSampl, &yMapCOMPreReSampl, &zMapCOMPreReSampl, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
681 this->reSampleMap ( settings );
684 xSampRate = this->xDimSize /
static_cast< proshade_single
> ( this->xTo - this->xFrom );
685 ySampRate = this->yDimSize /
static_cast< proshade_single
> ( this->yTo - this->yFrom );
686 zSampRate = this->zDimSize /
static_cast< proshade_single
> ( this->zTo - this->zFrom );
689 proshade_single xStartPosAfter =
static_cast< proshade_single
> ( this->xFrom ) * xSampRate;
690 proshade_single yStartPosAfter =
static_cast< proshade_single
> ( this->yFrom ) * ySampRate;
691 proshade_single zStartPosAfter =
static_cast< proshade_single
> ( this->zFrom ) * zSampRate;
694 proshade_single xMov =
static_cast< proshade_single
> ( xStartPosAfter - xStartPosBefore );
695 proshade_single yMov =
static_cast< proshade_single
> ( yStartPosAfter - yStartPosBefore );
696 proshade_single zMov =
static_cast< proshade_single
> ( zStartPosAfter - zStartPosBefore );
698 &this->xFrom, &this->xTo, &this->yFrom, &this->yTo, &this->zFrom, &this->zTo,
699 &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
702 proshade_double xMapCOMPostReSampl = 0.0, yMapCOMPostReSampl = 0.0, zMapCOMPostReSampl = 0.0;
703 ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xMapCOMPostReSampl, &yMapCOMPostReSampl, &zMapCOMPostReSampl, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
707 static_cast< proshade_single
> ( xMapCOMPreReSampl - xMapCOMPostReSampl ),
708 static_cast< proshade_single
> ( yMapCOMPreReSampl - yMapCOMPostReSampl ),
709 static_cast< proshade_single
> ( zMapCOMPreReSampl - zMapCOMPostReSampl ),
710 this->xDimSize, this->yDimSize, this->zDimSize,
711 static_cast< proshade_signed
> ( this->xDimIndices ),
712 static_cast< proshade_signed
> ( this->yDimIndices ),
713 static_cast< proshade_signed
> ( this->zDimIndices ) );
717 this->xDimSizeOriginal = this->xDimSize;
718 this->yDimSizeOriginal = this->yDimSize;
719 this->zDimSizeOriginal = this->zDimSize;
722 this->xDimIndicesOriginal = this->xDimIndices;
723 this->yDimIndicesOriginal = this->yDimIndices;
724 this->zDimIndicesOriginal = this->zDimIndices;
727 this->xAxisOriginOriginal = this->xAxisOrigin;
728 this->yAxisOriginOriginal = this->yAxisOrigin;
729 this->zAxisOriginOriginal = this->zAxisOrigin;
732 this->findMapCOM ( );
733 this->originalMapXCom = this->xCom;
734 this->originalMapYCom = this->yCom;
735 this->originalMapZCom = this->zCom;
759 gemmi::Structure pdbFile = gemmi::read_structure ( gemmi::MaybeGzipped ( this->fileName ) );
762 this->readInGemmi ( pdbFile, settings );
806 proshade_double xCOMPdb, yCOMPdb, zCOMPdb;
810 proshade_single xF, xT, yF, yT, zF, zT;
814 proshade_single xMov =
static_cast< proshade_single
> ( 20.0f - xF );
815 proshade_single yMov =
static_cast< proshade_single
> ( 20.0f - yF );
816 proshade_single zMov =
static_cast< proshade_single
> ( 20.0f - zF );
820 this->xDimSize =
static_cast< proshade_single
> ( xT - xF + 40.0f );
821 this->yDimSize =
static_cast< proshade_single
> ( yT - yF + 40.0f );
822 this->zDimSize =
static_cast< proshade_single
> ( zT - zF + 40.0f );
825 ProSHADE_internal_mapManip::generateMapFromPDB ( gemmiStruct, this->internalMap, settings->
requestedResolution, this->xDimSize, this->yDimSize, this->zDimSize, &this->xTo, &this->yTo, &this->zTo, settings->
forceP1, settings->
firstModelOnly );
828 if ( settings->
removeNegativeDensity ) {
for (
size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) {
if ( this->internalMap[iter] < 0.0 ) { this->internalMap[iter] = 0.0; } } }
831 this->setPDBMapValues ( );
834 proshade_double xCOMMap, yCOMMap, zCOMMap;
836 this->xDimSize, this->yDimSize, this->zDimSize,
837 this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
839 if ( gemmiStruct.models.size() > 1 )
841 xMov =
static_cast< proshade_single
> ( xCOMMap - xCOMPdb );
842 yMov =
static_cast< proshade_single
> ( yCOMMap - yCOMPdb );
843 zMov =
static_cast< proshade_single
> ( zCOMMap - zCOMPdb );
848 &this->xFrom, &this->xTo, &this->yFrom, &this->yTo, &this->zFrom, &this->zTo,
849 &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
851 static_cast< proshade_signed
> ( this->xDimIndices ),
static_cast< proshade_signed
> ( this->yDimIndices ),
852 static_cast< proshade_signed
> ( this->zDimIndices ) );
858 proshade_single xSampRate = this->xDimSize /
static_cast< proshade_single
> ( this->xTo - this->xFrom );
859 proshade_single ySampRate = this->yDimSize /
static_cast< proshade_single
> ( this->yTo - this->yFrom );
860 proshade_single zSampRate = this->zDimSize /
static_cast< proshade_single
> ( this->zTo - this->zFrom );
863 proshade_single xStartPosBefore =
static_cast< proshade_single
> ( this->xFrom ) * xSampRate;
864 proshade_single yStartPosBefore =
static_cast< proshade_single
> ( this->yFrom ) * ySampRate;
865 proshade_single zStartPosBefore =
static_cast< proshade_single
> ( this->zFrom ) * zSampRate;
868 proshade_double xMapCOMPreReSampl = 0.0, yMapCOMPreReSampl = 0.0, zMapCOMPreReSampl = 0.0;
869 ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xMapCOMPreReSampl, &yMapCOMPreReSampl, &zMapCOMPreReSampl, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
872 this->reSampleMap ( settings );
875 xSampRate = this->xDimSize /
static_cast< proshade_single
> ( this->xTo - this->xFrom );
876 ySampRate = this->yDimSize /
static_cast< proshade_single
> ( this->yTo - this->yFrom );
877 zSampRate = this->zDimSize /
static_cast< proshade_single
> ( this->zTo - this->zFrom );
880 proshade_single xStartPosAfter =
static_cast< proshade_single
> ( this->xFrom ) * xSampRate;
881 proshade_single yStartPosAfter =
static_cast< proshade_single
> ( this->yFrom ) * ySampRate;
882 proshade_single zStartPosAfter =
static_cast< proshade_single
> ( this->zFrom ) * zSampRate;
885 proshade_single xMovHlp =
static_cast< proshade_single
> ( xStartPosAfter - xStartPosBefore );
886 proshade_single yMovHlp =
static_cast< proshade_single
> ( yStartPosAfter - yStartPosBefore );
887 proshade_single zMovHlp =
static_cast< proshade_single
> ( zStartPosAfter - zStartPosBefore );
889 &this->xFrom, &this->xTo, &this->yFrom, &this->yTo, &this->zFrom, &this->zTo,
890 &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
893 proshade_double xMapCOMPostReSampl = 0.0, yMapCOMPostReSampl = 0.0, zMapCOMPostReSampl = 0.0;
894 ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xMapCOMPostReSampl, &yMapCOMPostReSampl, &zMapCOMPostReSampl, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
898 static_cast< proshade_single
> ( xMapCOMPreReSampl - xMapCOMPostReSampl ),
899 static_cast< proshade_single
> ( yMapCOMPreReSampl - yMapCOMPostReSampl ),
900 static_cast< proshade_single
> ( zMapCOMPreReSampl - zMapCOMPostReSampl ),
901 this->xDimSize, this->yDimSize, this->zDimSize,
902 static_cast< proshade_signed
> ( this->xDimIndices ),
903 static_cast< proshade_signed
> ( this->yDimIndices ),
904 static_cast< proshade_signed
> ( this->zDimIndices ) );
908 this->xDimSizeOriginal = this->xDimSize;
909 this->yDimSizeOriginal = this->yDimSize;
910 this->zDimSizeOriginal = this->zDimSize;
913 this->xDimIndicesOriginal = this->xDimIndices;
914 this->yDimIndicesOriginal = this->yDimIndices;
915 this->zDimIndicesOriginal = this->zDimIndices;
918 this->xAxisOriginOriginal = this->xAxisOrigin;
919 this->yAxisOriginOriginal = this->yAxisOrigin;
920 this->zAxisOriginOriginal = this->zAxisOrigin;
923 this->findMapCOM ( );
924 this->originalMapXCom = this->xCom;
925 this->originalMapYCom = this->yCom;
926 this->originalMapZCom = this->zCom;
950 this->xDimIndices =
static_cast< proshade_unsign
> ( this->xTo );
951 this->yDimIndices =
static_cast< proshade_unsign
> ( this->yTo );
952 this->zDimIndices =
static_cast< proshade_unsign
> ( this->zTo );
960 this->xGridIndices = this->xDimIndices;
961 this->yGridIndices = this->yDimIndices;
962 this->zGridIndices = this->zDimIndices;
965 this->xAxisOrder = 1;
966 this->yAxisOrder = 2;
967 this->zAxisOrder = 3;
970 this->xAxisOrigin = this->xFrom;
971 this->yAxisOrigin = this->yFrom;
972 this->zAxisOrigin = this->zFrom;
986 this->xTo = this->xFrom +
static_cast< proshade_signed
> ( this->xDimIndices ) - 1;
987 this->yTo = this->yFrom +
static_cast< proshade_signed
> ( this->yDimIndices ) - 1;
988 this->zTo = this->zFrom +
static_cast< proshade_signed
> ( this->zDimIndices ) - 1;
1008 gemmi::Grid<float> mapData;
1009 mapData.set_unit_cell (
static_cast< double > ( this->xDimSize ),
static_cast< double > ( this->yDimSize ),
static_cast< double > ( this->zDimSize ),
static_cast< double > ( this->aAngle ),
static_cast< double > ( this->bAngle ),
static_cast< double > ( this->cAngle ) );
1010 mapData.set_size_without_checking (
static_cast< int > ( this->xDimIndices ),
static_cast< int > ( this->yDimIndices ),
static_cast< int > ( this->zDimIndices ) );
1011 mapData.axis_order = gemmi::AxisOrder::XYZ;
1012 mapData.spacegroup = &gemmi::get_spacegroup_p1();
1015 gemmi::Ccp4<float> map;
1017 map.update_ccp4_header ( mode );
1021 this->xDimIndices, this->yDimIndices, this->zDimIndices,
1022 this->xDimSize, this->yDimSize, this->zDimSize,
1023 this->aAngle, this->bAngle, this->cAngle,
1024 this->xFrom, this->yFrom, this->zFrom,
1025 this->xAxisOrigin, this->yAxisOrigin, this->zAxisOrigin,
1026 this->xAxisOrder, this->yAxisOrder, this->zAxisOrder,
1027 this->xGridIndices, this->yGridIndices, this->zGridIndices,
1031 proshade_unsign arrPos = 0;
1032 for ( proshade_unsign uIt = 0; uIt < this->xDimIndices; uIt++ )
1034 for ( proshade_unsign vIt = 0; vIt < this->yDimIndices; vIt++ )
1036 for ( proshade_unsign wIt = 0; wIt < this->zDimIndices; wIt++ )
1038 arrPos = wIt + this->zDimIndices * ( vIt + this->yDimIndices * uIt );
1039 map.grid.set_value (
static_cast< int > ( uIt ),
static_cast< int > ( vIt ),
static_cast< int > ( wIt ),
static_cast<float> ( this->internalMap[arrPos] ) );
1045 map.update_ccp4_header ( mode,
true );
1048 map.write_ccp4_map ( fName );
1072 void ProSHADE_internal_data::ProSHADE_data::writePdb ( std::string fName, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double trsX, proshade_double trsY, proshade_double trsZ, proshade_double rotX, proshade_double rotY, proshade_double rotZ,
bool firstModel )
1077 throw ProSHADE_exception (
"Cannot write co-ordinate file if the input file did not contain co-ordinates.",
"EP00047", __FILE__, __LINE__, __func__,
"You have called the WritePDB function on structure which\n : was created by reading in a map. This is not allowed as\n : ProSHADE cannot create co-ordinates from map file." );
1081 gemmi::Structure pdbFile = gemmi::read_structure ( gemmi::MaybeGzipped ( this->fileName ) );
1084 this->writeGemmi ( fName, pdbFile, euA, euB, euG, trsX, trsY, trsZ, rotX, rotY, rotZ, firstModel );
1109 void ProSHADE_internal_data::ProSHADE_data::writeGemmi ( std::string fName, gemmi::Structure gemmiStruct, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double trsX, proshade_double trsY, proshade_double trsZ, proshade_double rotX, proshade_double rotY, proshade_double rotZ,
bool firstModel )
1112 if ( ( euA != 0.0 ) || ( euB != 0.0 ) || ( euG != 0.0 ) )
1122 std::ofstream outCoOrdFile;
1123 outCoOrdFile.open ( fName.c_str() );
1125 if ( outCoOrdFile.is_open() )
1127 gemmi::PdbWriteOptions opt;
1128 write_pdb ( gemmiStruct, outCoOrdFile, opt );
1132 std::stringstream hlpMessage;
1133 hlpMessage <<
"Failed to open the PDB file " << fName <<
" for output.";
1134 throw ProSHADE_exception ( hlpMessage.str().c_str(),
"EP00048", __FILE__, __LINE__, __func__,
"ProSHADE has failed to open the PDB output file. This is\n : likely caused by either not having the write privileges\n : to the required output path, or by making a mistake in\n : the path." );
1137 outCoOrdFile.close ( );
1154 proshade_double* hlpMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1158 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1160 hlpMap[iter] = this->internalMap[iter];
1161 this->internalMap[iter] = mask[iter];
1165 this->writeMap ( fName );
1168 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1170 this->internalMap[iter] = hlpMap[iter];
1194 proshade_signed arrayPos, invPos;
1197 proshade_double* hlpMap =
new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
1201 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1203 hlpMap[iter] = this->internalMap[iter];
1207 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
1209 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
1211 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
1214 arrayPos = zIt +
static_cast< proshade_signed
> ( this->zDimIndices ) * ( yIt +
static_cast< proshade_signed
> ( this->yDimIndices ) * xIt );
1215 invPos = (
static_cast< proshade_signed
> ( this->zDimIndices - 1 ) - zIt ) +
static_cast< proshade_signed
> ( this->zDimIndices ) * ( (
static_cast< proshade_signed
> ( this->yDimIndices - 1 ) - yIt ) +
static_cast< proshade_signed
> ( this->yDimIndices ) * (
static_cast< proshade_signed
> ( this->xDimIndices - 1 ) - xIt ) );
1218 this->internalMap[invPos] = hlpMap[arrayPos];
1248 std::vector<proshade_double> mapVals ( this->xDimIndices * this->yDimIndices * this->zDimIndices, 0.0 );
1251 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1253 mapVals.at(iter) = this->internalMap[iter];
1257 proshade_double* meanSD =
new proshade_double[2];
1261 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1263 this->internalMap[iter] = ( this->internalMap[iter] - meanSD[0] ) / meanSD[1];
1295 proshade_double* blurredMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1300 this->xDimSize, this->yDimSize, this->zDimSize, settings->
blurFactor );
1306 if ( settings->
saveMask ) {
if ( settings->
maskFileName ==
"" ) { this->writeMask (
"proshade_mask.map", blurredMap ); }
else { std::stringstream ss; ss << settings->
maskFileName <<
"_" << this->inputOrder <<
".map"; this->writeMask ( ss.str(), blurredMap ); } }
1309 delete[] blurredMap;
1337 for ( proshade_unsign iter = 0; iter < 6; iter++ ) { ret[iter] = settings->
forceBounds[iter]; }
1344 static_cast< proshade_signed
> ( this->xDimIndices ),
1345 static_cast< proshade_signed
> ( this->yDimIndices ),
1346 static_cast< proshade_signed
> ( this->zDimIndices ),
1351 this->xDimSize, this->yDimSize, this->zDimSize, ret, settings->
boundsExtraSpace );
1357 std::stringstream ssHlp;
1358 ssHlp <<
"New boundaries are: " << ret[1] - ret[0] + 1 <<
" x " << ret[3] - ret[2] + 1 <<
" x " << ret[5] - ret[4] + 1;
1364 for ( proshade_unsign iter = 0; iter < 6; iter++ ) { settings->
forceBounds[iter] = ret[iter]; }
1394 newStr->
fileType = ProSHADE_internal_io::MAP;
1397 newStr->
xDimIndices =
static_cast< proshade_unsign
> ( newBounds[1] ) -
static_cast< proshade_unsign
> ( newBounds[0] ) + 1;
1398 newStr->
yDimIndices =
static_cast< proshade_unsign
> ( newBounds[3] ) -
static_cast< proshade_unsign
> ( newBounds[2] ) + 1;
1399 newStr->
zDimIndices =
static_cast< proshade_unsign
> ( newBounds[5] ) -
static_cast< proshade_unsign
> ( newBounds[4] ) + 1;
1401 newStr->
aAngle = this->aAngle;
1402 newStr->
bAngle = this->aAngle;
1403 newStr->
cAngle = this->aAngle;
1405 newStr->
xDimSize =
static_cast<proshade_single
> ( newStr->
xDimIndices ) * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices ) );
1406 newStr->
yDimSize =
static_cast<proshade_single
> ( newStr->
yDimIndices ) * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices ) );
1407 newStr->
zDimSize =
static_cast<proshade_single
> ( newStr->
zDimIndices ) * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices ) );
1417 newStr->
xAxisOrigin = this->xAxisOrigin + newBounds[0];
1418 newStr->
yAxisOrigin = this->yAxisOrigin + newBounds[2];
1419 newStr->
zAxisOrigin = this->zAxisOrigin + newBounds[4];
1421 newStr->
xFrom = this->xFrom + newBounds[0];
1422 newStr->
yFrom = this->yFrom + newBounds[2];
1423 newStr->
zFrom = this->zFrom + newBounds[4];
1425 newStr->
xTo = this->xTo - (
static_cast< proshade_signed
> ( this->xDimIndices - 1 ) - newBounds[1] );
1426 newStr->
yTo = this->yTo - (
static_cast< proshade_signed
> ( this->yDimIndices - 1 ) - newBounds[3] );
1427 newStr->
zTo = this->zTo - (
static_cast< proshade_signed
> ( this->zDimIndices - 1 ) - newBounds[5] );
1436 this->xDimIndices, this->yDimIndices, this->zDimIndices, newStr->
internalMap, this->internalMap );
1456 proshade_single* changeVals =
new proshade_single[6];
1462 this->xDimSize, this->yDimSize, this->zDimSize, changeVals );
1472 this->xDimSize, this->yDimSize, this->zDimSize, changeVals );
1477 this->xDimIndices +=
static_cast<proshade_unsign
> ( changeVals[0] );
1478 this->yDimIndices +=
static_cast<proshade_unsign
> ( changeVals[1] );
1479 this->zDimIndices +=
static_cast<proshade_unsign
> ( changeVals[2] );
1481 this->xGridIndices = this->xDimIndices;
1482 this->yGridIndices = this->yDimIndices;
1483 this->zGridIndices = this->zDimIndices;
1485 this->xTo +=
static_cast<proshade_unsign
> ( changeVals[0] );
1486 this->yTo +=
static_cast<proshade_unsign
> ( changeVals[1] );
1487 this->zTo +=
static_cast<proshade_unsign
> ( changeVals[2] );
1489 this->xDimSize = changeVals[3];
1490 this->yDimSize = changeVals[4];
1491 this->zDimSize = changeVals[5];
1494 proshade_single xMov = -( (
static_cast<proshade_single
> ( this->xFrom ) * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices ) - changeVals[0] ) ) -
1495 (
static_cast<proshade_single
> ( this->xFrom ) * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices ) ) ) );
1496 proshade_single yMov = -( (
static_cast<proshade_single
> ( this->yFrom ) * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices ) - changeVals[1] ) ) -
1497 (
static_cast<proshade_single
> ( this->yFrom ) * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices ) ) ) );
1498 proshade_single zMov = -( (
static_cast<proshade_single
> ( this->zFrom ) * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices ) - changeVals[2] ) ) -
1499 (
static_cast<proshade_single
> ( this->zFrom ) * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices ) ) ) );
1503 &this->yFrom, &this->yTo, &this->zFrom, &this->zTo, &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
1506 static_cast< proshade_signed
> ( this->xDimIndices ),
static_cast< proshade_signed
> ( this->yDimIndices ),
static_cast< proshade_signed
> ( this->zDimIndices ) );
1509 delete[] changeVals;
1531 proshade_double xCOM = 0.0;
1532 proshade_double yCOM = 0.0;
1533 proshade_double zCOM = 0.0;
1537 this->xDimSize, this->yDimSize, this->xDimSize, this->xFrom,
1538 this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
1541 proshade_single xSampRate =
static_cast< proshade_single
> ( this->xDimSize ) /
static_cast< proshade_single
> ( this->xTo - this->xFrom );
1542 proshade_single ySampRate =
static_cast< proshade_single
> ( this->yDimSize ) /
static_cast< proshade_single
> ( this->yTo - this->yFrom );
1543 proshade_single zSampRate =
static_cast< proshade_single
> ( this->zDimSize ) /
static_cast< proshade_single
> ( this->zTo - this->zFrom );
1546 xCOM /=
static_cast< proshade_double
> ( xSampRate );
1547 yCOM /=
static_cast< proshade_double
> ( ySampRate );
1548 zCOM /=
static_cast< proshade_double
> ( zSampRate );
1550 xCOM -=
static_cast< proshade_double
> ( this->xFrom );
1551 yCOM -=
static_cast< proshade_double
> ( this->yFrom );
1552 zCOM -=
static_cast< proshade_double
> ( this->zFrom );
1555 proshade_double xDist = ( (
static_cast<proshade_double
> ( this->xDimIndices ) / 2.0 ) - xCOM ) *
static_cast<proshade_double
> ( this->xDimSize ) /
static_cast<proshade_double
> ( this->xDimIndices );
1556 proshade_double yDist = ( (
static_cast<proshade_double
> ( this->yDimIndices ) / 2.0 ) - yCOM ) *
static_cast<proshade_double
> ( this->yDimSize ) /
static_cast<proshade_double
> ( this->yDimIndices );
1557 proshade_double zDist = ( (
static_cast<proshade_double
> ( this->zDimIndices ) / 2.0 ) - zCOM ) *
static_cast<proshade_double
> ( this->zDimSize ) /
static_cast<proshade_double
> ( this->zDimIndices );
1561 static_cast< proshade_single
> ( xDist ),
1562 static_cast< proshade_single
> ( yDist ),
1563 static_cast< proshade_single
> ( zDist ),
1564 this->xDimSize, this->yDimSize, this->zDimSize,
1565 static_cast< proshade_signed
> ( this->xDimIndices ),
1566 static_cast< proshade_signed
> ( this->yDimIndices ),
1567 static_cast< proshade_signed
> ( this->zDimIndices ) );
1570 this->mapCOMProcessChangeX -= xDist;
1571 this->mapCOMProcessChangeY -= yDist;
1572 this->mapCOMProcessChangeZ -= zDist;
1593 std::stringstream hlpSS;
1594 hlpSS <<
"Adding extra " << settings->
addExtraSpace <<
" angstroms.";
1603 this->xDimSize += 2 *
static_cast<proshade_single
> ( xAddIndices ) * this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices );
1604 this->yDimSize += 2 *
static_cast<proshade_single
> ( yAddIndices ) * this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices );
1605 this->zDimSize += 2 *
static_cast<proshade_single
> ( zAddIndices ) * this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices );
1607 this->xDimIndices += 2 * xAddIndices;
1608 this->yDimIndices += 2 * yAddIndices;
1609 this->zDimIndices += 2 * zAddIndices;
1611 this->xGridIndices = this->xDimIndices;
1612 this->yGridIndices = this->yDimIndices;
1613 this->zGridIndices = this->zDimIndices;
1615 this->xAxisOrigin -= xAddIndices;
1616 this->yAxisOrigin -= yAddIndices;
1617 this->zAxisOrigin -= zAddIndices;
1619 this->xFrom -= xAddIndices;
1620 this->yFrom -= yAddIndices;
1621 this->zFrom -= zAddIndices;
1623 this->xTo += xAddIndices;
1624 this->yTo += yAddIndices;
1625 this->zTo += zAddIndices;
1628 proshade_double* newMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1632 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1638 proshade_unsign newMapIndex, oldMapIndex;
1639 for ( proshade_unsign xIt = 0; xIt < (this->xDimIndices - xAddIndices); xIt++ )
1642 if ( xIt < xAddIndices ) {
continue; }
1644 for ( proshade_unsign yIt = 0; yIt < (this->yDimIndices - yAddIndices); yIt++ )
1647 if ( yIt < yAddIndices ) {
continue; }
1649 for ( proshade_unsign zIt = 0; zIt < (this->zDimIndices - zAddIndices); zIt++ )
1652 if ( zIt < zAddIndices ) {
continue; }
1655 newMapIndex = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
1656 oldMapIndex = (zIt - zAddIndices) + (this->zDimIndices - ( 2 * zAddIndices ) ) * ( (yIt - yAddIndices) + (this->yDimIndices - ( 2 * yAddIndices ) ) * (xIt - xAddIndices) );
1658 newMap[newMapIndex] = this->internalMap[oldMapIndex];
1664 delete[] this->internalMap;
1666 this->internalMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1669 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1671 this->internalMap[iter] = newMap[iter];
1702 if ( settings->
invertMap ) { this->invertMirrorMap ( settings ); }
1706 if ( settings->
normaliseMap ) { this->normaliseMap ( settings ); }
1710 if ( settings->
maskMap ) { this->maskMap ( settings ); }
1714 if ( settings->
moveToCOM ) { this->centreMapOnCOM ( settings ); }
1718 if ( settings->
addExtraSpace != 0.0f ) { this->addExtraSpace ( settings ); }
1745 if ( this->spherePos.size() != 0 )
1747 std::stringstream hlpSS;
1748 hlpSS <<
"The sphere distances were determined as " << this->spherePos.at(0);
1749 for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ ) { hlpSS <<
"; " << this->spherePos.at(iter); }
1750 hlpSS <<
" Angstroms.";
1756 proshade_unsign maxDim =
static_cast< proshade_unsign
> ( std::max ( this->xDimSize, std::max ( this->yDimSize, this->zDimSize ) ) );
1757 proshade_unsign minDim =
static_cast< proshade_unsign
> ( std::min ( this->xDimSize, std::min ( this->yDimSize, this->zDimSize ) ) );
1758 proshade_unsign midDim =
static_cast< proshade_unsign
> ( 0 );
1759 if ( ( this->xDimSize <
static_cast< proshade_single
> ( maxDim ) ) && ( this->xDimSize >
static_cast< proshade_single
> ( minDim ) ) ) { midDim =
static_cast< proshade_unsign
> ( this->xDimSize ); }
1760 else if ( ( this->yDimSize <
static_cast< proshade_single
> ( maxDim ) ) && ( this->yDimSize >
static_cast< proshade_single
> ( minDim ) ) ) { midDim =
static_cast< proshade_unsign
> ( this->yDimSize ); }
1761 else { midDim =
static_cast< proshade_unsign
> ( this->zDimSize ); }
1763 proshade_single maxDiag =
static_cast< proshade_single
> ( std::sqrt ( std::pow (
static_cast<proshade_single
> ( maxDim ), 2.0 ) +
1764 std::pow (
static_cast<proshade_single
> ( midDim ), 2.0 ) ) );
1767 for ( proshade_single iter = 0.5f; ( iter * settings->
maxSphereDists ) < ( maxDiag / 2.0f ); iter += 1.0f )
1773 this->noSpheres =
static_cast<proshade_unsign
> ( this->spherePos.size() );
1776 std::stringstream hlpSS;
1777 hlpSS <<
"The sphere distances were determined as " << this->spherePos.at(0);
1778 for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ ) { hlpSS <<
"; " << this->spherePos.at(iter); }
1779 hlpSS <<
" Angstroms.";
1806 this->xDimSize, this->yDimSize, this->zDimSize );
1810 this->getSpherePositions ( settings );
1815 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ )
1817 std::stringstream ss;
1818 ss <<
"Now mapping sphere " << iter <<
" .";
1822 this->xDimSize, this->yDimSize, this->zDimSize, iter,
1824 this->internalMap, &this->maxShellBand );
1849 this->sphericalHarmonics =
new proshade_complex* [this->noSpheres];
1851 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
1853 this->sphericalHarmonics[iter] =
new proshade_complex [(this->spheres[iter]->getLocalBandwidth() * 2) * (this->spheres[iter]->getLocalBandwidth() * 2)];
1858 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
1861 std::stringstream ss;
1862 ss <<
"Now decomposing sphere " << iter <<
". " <<
"( Band is: " << this->spheres[iter]->getLocalBandwidth() <<
").";
1886 return ( a[6] > b[6] );
1906 if ( settings->axisErrToleranceDefault )
1908 settings->
axisErrTolerance = std::min ( std::max ( 0.01, ( 2.0 * M_PI ) /
static_cast< proshade_double
> ( this->maxShellBand ) ), 0.05 );
1912 fftw_complex *mapData, *origCoeffs, *fCoeffs;
1913 fftw_plan planForwardFourier;
1914 proshade_double **bindata;
1915 proshade_signed *binIndexing, *binCounts, noBins;
1916 this->prepareFSCFourierMemory ( mapData, origCoeffs, fCoeffs, binIndexing, &noBins, bindata, binCounts, &planForwardFourier );
1922 std::stringstream hlpSS;
1927 proshade_double symThres = 0.0;
1928 std::vector< proshade_double* > CSyms = this->findRequestedCSymmetryFromAngleAxis ( settings, settings->
requestedSymmetryFold, &symThres );
1931 bool possible010PIIssue =
false;
1933 if ( possible010PIIssue )
1935 proshade_double* addAxis =
new proshade_double[7];
1936 addAxis[0] = 2.0; addAxis[1] = 0.0; addAxis[2] = 1.0; addAxis[3] = 0.0; addAxis[4] = M_PI; addAxis[5] = -999.9; addAxis[6] = -1.0;
1942 for (
size_t cIt = 0; cIt < CSyms.size(); cIt++ ) {
const FloatingPoint< proshade_double > lhs ( CSyms.at(cIt)[5] ), rhs ( -999.9 );
if ( ( CSyms.at(cIt)[5] > settings->
peakThresholdMin ) || ( lhs.AlmostEquals ( rhs ) ) ) { this->
computeFSC ( settings, &CSyms, cIt, mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts ); } }
1949 if ( CSyms.size() > 0 )
1951 bool passedTests =
false;
1952 for (
size_t cIt = 0; cIt < CSyms.size(); cIt++ )
1960 this->saveDetectedSymmetries ( settings, &CSyms, allCs );
1981 delete[] origCoeffs;
1983 fftw_destroy_plan ( planForwardFourier );
1984 delete[] binIndexing;
1985 for (
size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) {
delete[] bindata[binIt]; }
1997 std::vector< proshade_double* > CSyms = getCyclicSymmetriesListFromAngleAxis ( settings );
2000 bool possible010PIIssue =
false;
2002 if ( possible010PIIssue )
2004 proshade_double* addAxis =
new proshade_double[7];
2005 addAxis[0] = 2.0; addAxis[1] = 0.0; addAxis[2] = 1.0; addAxis[3] = 0.0; addAxis[4] = M_PI; addAxis[5] = -999.9; addAxis[6] = -1.0;
2010 for (
size_t cIt = 0; cIt < CSyms.size(); cIt++ ) {
const FloatingPoint< proshade_double > lhs ( CSyms.at(cIt)[5] ), rhs ( -999.9 );
if ( lhs.AlmostEquals ( rhs ) ) { this->
computeFSC ( settings, &CSyms, cIt, mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts ); } }
2013 std::stringstream ss;
2014 ss <<
"Detected " << CSyms.size() <<
" C symmetries.";
2018 if ( this->sphereMappedRotFun.size() < 1 )
2020 throw ProSHADE_exception (
"Rotation function was not converted into angle-axis space.",
"ES00062", __FILE__, __LINE__, __func__,
"It seems that the convertRotationFunction() function was\n : not yet called. Therefore, there are no data to detect the\n : symmetry from; please call the convertRotationFunction()\n : function before the detectSymmetryFromAngleAxisSpace()\n : function." );
2026 throw ProSHADE_exception (
"Requested symmetry supplied, but not recognised.",
"ES00032", __FILE__, __LINE__, __func__,
"There are only the following value allowed for the\n : symmetry type request: \"C\", \"D\", \"T\", \"O\" and \"I\". Any\n : other value will result in this error." );
2033 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
2034 std::vector< proshade_double* > ISyms;
2035 std::vector< proshade_double* > OSyms = this->getPredictedOctahedralSymmetriesList ( settings, &CSyms );
2036 std::vector< proshade_double* > TSyms = this->getPredictedTetrahedralSymmetriesList ( settings, &CSyms );
2039 proshade_double fscMax = 0.0;
2040 size_t fscMaxInd = 0;
2041 std::vector < std::vector< proshade_double* > > ISymsHlp = this->getPredictedIcosahedralSymmetriesList ( settings, &CSyms );
2043 for (
size_t icoIt = 0; icoIt < ISymsHlp.size(); icoIt++ )
2046 if ( ISymsHlp.at(icoIt).size() != 31 ) {
continue; }
2049 proshade_double fscVal = 0.0;
2050 proshade_double fscValAvg = 0.0;
2053 for (
size_t aIt = 0; aIt < ISymsHlp.at(icoIt).size(); aIt++ )
2056 proshade_signed matchedPos =
ProSHADE_internal_symmetry::addAxisUnlessSame (
static_cast< proshade_unsign
> ( ISymsHlp.at(icoIt).at(aIt)[0] ), ISymsHlp.at(icoIt).at(aIt)[1], ISymsHlp.at(icoIt).at(aIt)[2], ISymsHlp.at(icoIt).at(aIt)[3], ISymsHlp.at(icoIt).at(aIt)[5], ISymsHlp.at(icoIt).at(aIt)[6], &CSyms, settings->
axisErrTolerance );
2059 fscVal = this->
computeFSC ( settings, &CSyms,
static_cast< size_t > ( matchedPos ), mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts );
2060 ISymsHlp.at(fscMaxInd).at(aIt)[6] = fscVal;
2061 fscValAvg += fscVal;
2068 if ( fscValAvg > fscMax )
2079 for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ISymsHlp.at(fscMaxInd).size() ); retIt++ )
2082 proshade_signed matchedPos =
ProSHADE_internal_symmetry::addAxisUnlessSame (
static_cast< proshade_unsign
> ( ISymsHlp.at(fscMaxInd).at(retIt)[0] ), ISymsHlp.at(fscMaxInd).at(retIt)[1], ISymsHlp.at(fscMaxInd).at(retIt)[2], ISymsHlp.at(fscMaxInd).at(retIt)[3], ISymsHlp.at(fscMaxInd).at(retIt)[5], ISymsHlp.at(fscMaxInd).at(retIt)[6], &CSyms, settings->
axisErrTolerance );
2091 this->saveRecommendedSymmetry ( settings, &CSyms, &DSyms, &TSyms, &OSyms, &ISyms, axes, mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts );
2097 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
2098 this->saveRequestedSymmetryD ( settings, &DSyms, axes, mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts );
2104 std::vector< proshade_double* > TSyms = this->getPredictedTetrahedralSymmetriesList ( settings, &CSyms );
2106 if ( TSyms.size() == 7 )
2109 proshade_double fscVal = 0.0;
2110 proshade_double fscValAvg = 0.0;
2113 for (
size_t tIt = 0; tIt < 7; tIt++ ) {
if ( CSyms.at(settings->
allDetectedTAxes.at(tIt))[5] > settings->
peakThresholdMin ) { fscVal = this->
computeFSC ( settings, &CSyms, settings->
allDetectedTAxes.at(tIt), mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2131 std::vector< proshade_double* > OSyms = this->getPredictedOctahedralSymmetriesList ( settings, &CSyms );
2133 if ( OSyms.size() == 13 )
2136 proshade_double fscVal = 0.0;
2137 proshade_double fscValAvg = 0.0;
2140 for (
size_t oIt = 0; oIt < 13; oIt++ ) {
if ( CSyms.at(settings->
allDetectedOAxes.at(oIt))[5] > settings->
peakThresholdMin ) { fscVal = this->
computeFSC ( settings, &CSyms, settings->
allDetectedOAxes.at(oIt), mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2158 proshade_double fscMax = 0.0;
2159 size_t fscMaxInd = 0;
2160 std::vector < std::vector< proshade_double* > > ISymsHlp = this->getPredictedIcosahedralSymmetriesList ( settings, &CSyms );
2161 std::vector< proshade_double* > ISyms;
2163 for (
size_t icoIt = 0; icoIt < ISymsHlp.size(); icoIt++ )
2166 if ( ISymsHlp.at(icoIt).size() != 31 ) {
continue; }
2169 proshade_double fscVal = 0.0;
2170 proshade_double fscValAvg = 0.0;
2173 for (
size_t aIt = 0; aIt < ISymsHlp.at(icoIt).size(); aIt++ )
2176 proshade_signed matchedPos =
ProSHADE_internal_symmetry::addAxisUnlessSame (
static_cast< proshade_unsign
> ( ISymsHlp.at(icoIt).at(aIt)[0] ), ISymsHlp.at(icoIt).at(aIt)[1], ISymsHlp.at(icoIt).at(aIt)[2], ISymsHlp.at(icoIt).at(aIt)[3], ISymsHlp.at(icoIt).at(aIt)[5], ISymsHlp.at(icoIt).at(aIt)[6], &CSyms, settings->
axisErrTolerance );
2179 fscVal = this->
computeFSC ( settings, &CSyms,
static_cast< size_t > ( matchedPos ), mapData, fCoeffs, origCoeffs, &planForwardFourier, noBins, binIndexing, bindata, binCounts );
2180 ISymsHlp.at(fscMaxInd).at(aIt)[6] = fscVal;
2181 fscValAvg += fscVal;
2188 if ( fscValAvg > fscMax )
2199 for ( proshade_unsign retIt = 0; retIt < static_cast < proshade_unsign > ( ISymsHlp.at(fscMaxInd).size() ); retIt++ )
2202 proshade_signed matchedPos =
ProSHADE_internal_symmetry::addAxisUnlessSame (
static_cast< proshade_unsign
> ( ISymsHlp.at(fscMaxInd).at(retIt)[0] ), ISymsHlp.at(fscMaxInd).at(retIt)[1], ISymsHlp.at(fscMaxInd).at(retIt)[2], ISymsHlp.at(fscMaxInd).at(retIt)[3], ISymsHlp.at(fscMaxInd).at(retIt)[5], ISymsHlp.at(fscMaxInd).at(retIt)[6], &CSyms, settings->
axisErrTolerance );
2211 for (
size_t gIt = 0; gIt < ISymsHlp.size(); gIt++ ) {
for (
size_t aIt = 0; aIt < ISymsHlp.at(gIt).size(); aIt++ ) {
delete[] ISymsHlp.at(gIt).at(aIt); } }
2225 this->saveDetectedSymmetries ( settings, &CSyms, allCs );
2229 delete[] origCoeffs;
2231 fftw_destroy_plan ( planForwardFourier );
2232 delete[] binIndexing;
2233 for (
size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) {
delete[] bindata[binIt]; }
2258 bool isArgSameAsSettings =
true;
2261 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSyms->size() ); cIt++ )
2264 std::vector< proshade_double > nextSym;
2275 if ( ( cIt == 0 ) && ( settings->
allDetectedCAxes.size() == 0 ) ) { isArgSameAsSettings =
false; }
2280 delete[] CSyms->at(cIt);
2304 proshade_double threshold = 0.0;
2305 proshade_signed totSize =
static_cast< proshade_signed
> ( ( 1.0 / step ) + 1 );
2306 std::vector< std::pair < proshade_double, proshade_unsign > > vals;
2307 std::vector< proshade_double > hist (
static_cast< unsigned long int > ( totSize ), 0.0 );
2308 proshade_unsign histPos = 0;
2311 if ( windowSize % 2 == 0 ) { windowSize += 1; }
2314 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( CSym->size() ); symIt++ ) { vals.emplace_back ( std::pair < proshade_double, proshade_unsign > ( CSym->at(symIt)[peakPos], symIt ) ); }
2317 for ( proshade_double it = 0.0; it <= 1.0; it = it + step )
2319 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( vals.size() ); symIt++ )
2322 if ( ( vals.at(symIt).first > it ) && ( vals.at(symIt).first <= ( it + step ) ) ) { hist.at(histPos) += 1.0; }
2336 proshade_signed bestHistPos;
2337 if ( peaks.size() > 0 ) { bestHistPos = peaks.at(peaks.size()-1) + ( ( windowSize - 1 ) / 2 ); }
2338 else { bestHistPos = 0.0; }
2340 threshold = (
static_cast< proshade_double
> ( bestHistPos ) * step ) - (
static_cast< proshade_double
> ( windowSize ) * step );
2343 return ( threshold );
2363 bindata =
new proshade_double*[*noBins];
2364 binCounts =
new proshade_signed [*noBins];
2367 for (
size_t binIt = 0; binIt < static_cast< size_t > ( *noBins ); binIt++ )
2369 bindata[binIt] =
new proshade_double[12];
2374 mapData =
new fftw_complex [this->xDimIndices * this->yDimIndices * this->zDimIndices];
2375 origCoeffs =
new fftw_complex [this->xDimIndices * this->yDimIndices * this->zDimIndices];
2376 fCoeffs =
new fftw_complex [this->xDimIndices * this->yDimIndices * this->zDimIndices];
2386 *planForwardFourier = fftw_plan_dft_3d (
static_cast< int > ( this->xDimIndices ),
static_cast< int > ( this->yDimIndices ),
static_cast< int > ( this->zDimIndices ), mapData, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
2389 for (
size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { mapData[iter][0] = this->internalMap[iter]; mapData[iter][1] = 0.0; }
2390 fftw_execute ( *planForwardFourier );
2391 for (
size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { origCoeffs[iter][0] = fCoeffs[iter][0]; origCoeffs[iter][1] = fCoeffs[iter][1]; }
2425 proshade_double
ProSHADE_internal_data::ProSHADE_data::computeFSC (
ProSHADE_settings* settings, std::vector< proshade_double* >* CSym,
size_t symIndex, fftw_complex* mapData, fftw_complex* fCoeffs, fftw_complex* origCoeffs, fftw_plan* planForwardFourier, proshade_signed noBins, proshade_signed *binIndexing, proshade_double**& bindata, proshade_signed*& binCounts )
2428 if ( symIndex >= CSym->size() )
2430 std::cerr <<
"The supplied symmetry axes vector does not contain element number " << symIndex <<
". Returning FSC 0.0." << std::endl;
2435 const FloatingPoint< proshade_double > lhs1 ( CSym->at(symIndex)[6] ), rhs1 ( -1.0 );
2436 if ( !lhs1.AlmostEquals ( rhs1 ) ) {
return ( CSym->at(symIndex)[6] ); }
2439 std::stringstream ss2;
2440 ss2 <<
"Computing FSC for symmetry C" << CSym->at(symIndex)[0] <<
" ( " << CSym->at(symIndex)[1] <<
" ; " << CSym->at(symIndex)[2] <<
" ; " << CSym->at(symIndex)[3] <<
" ) with peak height " << CSym->at(symIndex)[5];
2444 proshade_double *rotMap;
2447 proshade_double averageFSC = 0.0;
2448 for ( proshade_double rotIter = 1.0; rotIter < CSym->at(symIndex)[0]; rotIter += 1.0 )
2451 this->rotateMapRealSpace ( CSym->at(symIndex)[1], CSym->at(symIndex)[2], CSym->at(symIndex)[3], ( ( 2.0 * M_PI ) / CSym->at(symIndex)[0] ) * rotIter, rotMap );
2454 for (
size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { mapData[iter][0] = rotMap[iter]; mapData[iter][1] = 0.0; }
2455 fftw_execute ( *planForwardFourier );
2458 for (
size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) {
for (
size_t valIt = 0; valIt < 12; valIt++ ) { bindata[binIt][valIt] = 0.0; } }
2459 for (
size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) { binCounts[binIt] = 0; }
2462 averageFSC +=
ProSHADE_internal_maths::computeFSC ( origCoeffs, fCoeffs, this->xDimIndices, this->yDimIndices, this->zDimIndices, noBins, binIndexing, bindata, binCounts );
2469 averageFSC /= ( CSym->at(symIndex)[0] - 1.0 );
2472 CSym->at(symIndex)[6] = averageFSC;
2475 std::stringstream ss3;
2476 ss3 <<
"FSC value is " << averageFSC <<
" .";
2480 return ( averageFSC );
2513 const FloatingPoint< proshade_double > lhs1 ( sym[6] ), rhs1 ( -1.0 );
2514 if ( !lhs1.AlmostEquals ( rhs1 ) ) {
return ( sym[6] ); }
2517 std::stringstream ss2;
2518 ss2 <<
"Computing FSC for symmetry C" << sym[0] <<
" ( " << sym[1] <<
" ; " << sym[2] <<
" ; " << sym[3] <<
" ) with peak height " << sym[5];
2522 proshade_double *rotMap;
2525 proshade_double averageFSC = 0.0;
2526 for ( proshade_double rotIter = 1.0; rotIter < sym[0]; rotIter += 1.0 )
2529 this->rotateMapRealSpace ( sym[1], sym[2], sym[3], ( ( 2.0 * M_PI ) / sym[0] ) * rotIter, rotMap );
2532 for (
size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { mapData[iter][0] = rotMap[iter]; mapData[iter][1] = 0.0; }
2533 fftw_execute ( *planForwardFourier );
2536 for (
size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) {
for (
size_t valIt = 0; valIt < 12; valIt++ ) { bindata[binIt][valIt] = 0.0; } }
2537 for (
size_t binIt = 0; binIt < static_cast< size_t > ( noBins ); binIt++ ) { binCounts[binIt] = 0; }
2540 averageFSC +=
ProSHADE_internal_maths::computeFSC ( origCoeffs, fCoeffs, this->xDimIndices, this->yDimIndices, this->zDimIndices, noBins, binIndexing, bindata, binCounts );
2547 averageFSC /= ( sym[0] - 1.0 );
2550 sym[6] = averageFSC;
2553 std::stringstream ss3;
2554 ss3 <<
"FSC value is " << averageFSC <<
" .";
2558 return ( averageFSC );
2596 void ProSHADE_internal_data::ProSHADE_data::saveRecommendedSymmetry (
ProSHADE_settings* settings, std::vector< proshade_double* >* CSym, std::vector< proshade_double* >* DSym, std::vector< proshade_double* >* TSym, std::vector< proshade_double* >* OSym, std::vector< proshade_double* >* ISym, std::vector< proshade_double* >* axes, fftw_complex* mapData, fftw_complex* origCoeffs, fftw_complex* fCoeffs, fftw_plan* planForwardFourier, proshade_signed noBins, proshade_signed* binIndexing, proshade_double** bindata, proshade_signed* binCounts )
2602 if ( CSym->size() == 0 )
2610 proshade_double step = 0.01;
2611 proshade_double sigma = 0.03;
2612 proshade_signed windowSize = 9;
2613 proshade_double bestHistPeakStart = this->findTopGroupSmooth ( CSym, 5, step, sigma, windowSize );
2617 proshade_unsign noPassed = 0;
for (
size_t cIt = 0; cIt < CSym->size(); cIt++ ) {
if ( CSym->at(cIt)[5] > bestHistPeakStart ) { noPassed += 1; } }
2618 std::stringstream ss;
2619 ss <<
"Smoothening has resolved in " << noPassed <<
" C symmetries.";
2624 bool alreadyDecided =
false;
2625 if ( ISym->size() == 31 )
2628 proshade_double fscVal = 0.0;
2629 proshade_double fscValAvg = 0.0;
2632 for (
size_t iIt = 0; iIt < 31; iIt++ ) {
if ( CSym->at(settings->
allDetectedIAxes.at(iIt))[5] > bestHistPeakStart ) { fscVal = this->
computeFSC ( settings, CSym, settings->
allDetectedIAxes.at(iIt), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2645 alreadyDecided =
true;
2650 if ( ( OSym->size() == 13 ) && !alreadyDecided )
2653 proshade_double fscVal = 0.0;
2654 proshade_double fscValAvg = 0.0;
2657 for (
size_t oIt = 0; oIt < 13; oIt++ ) {
if ( CSym->at(settings->
allDetectedOAxes.at(oIt))[5] > bestHistPeakStart ) { fscVal = this->
computeFSC ( settings, CSym, settings->
allDetectedOAxes.at(oIt), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2670 alreadyDecided =
true;
2675 if ( ( TSym->size() == 7 ) && !alreadyDecided )
2678 proshade_double fscVal = 0.0;
2679 proshade_double fscValAvg = 0.0;
2682 for (
size_t tIt = 0; tIt < 7; tIt++ ) {
if ( CSym->at(settings->
allDetectedTAxes.at(tIt))[5] > bestHistPeakStart ) { fscVal = this->
computeFSC ( settings, CSym, settings->
allDetectedTAxes.at(tIt), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts ); fscValAvg += fscVal; } }
2695 alreadyDecided =
true;
2700 if ( ( settings->
allDetectedDAxes.size() > 0 ) && ( DSym->size() > 0 ) && !alreadyDecided )
2703 proshade_signed bestD = -1;
2704 proshade_unsign bestFold = 0;
2710 if ( dIt > 20 ) {
continue; }
2713 const FloatingPoint< proshade_double > lhs999a ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(0))[5] ), lhs999b ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(1))[5] ), rhs999 (
static_cast< proshade_double
> ( -999.9 ) );
2714 if ( ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(0))[5] < bestHistPeakStart ) && !( lhs999a.AlmostEquals( rhs999 ) ) ) {
continue; }
2715 if ( ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(1))[5] < bestHistPeakStart ) && !( lhs999b.AlmostEquals( rhs999 ) ) ) {
continue; }
2718 this->
computeFSC ( settings, CSym, settings->
allDetectedDAxes.at(dIt).at(0), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2719 this->
computeFSC ( settings, CSym, settings->
allDetectedDAxes.at(dIt).at(1), mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2723 proshade_double bestHistFSCStart = this->findTopGroupSmooth ( CSym, 5, step, sigma, windowSize );
2729 const FloatingPoint< proshade_double > lhs999a2 ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(0))[5] ), lhs999b2 ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(1))[5] ), rhs999 (
static_cast< proshade_double
> ( -999.9 ) );
2730 if ( ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(0))[5] < bestHistPeakStart ) && !( lhs999a2.AlmostEquals( rhs999 ) ) ) {
continue; }
2731 if ( ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(1))[5] < bestHistPeakStart ) && !( lhs999b2.AlmostEquals( rhs999 ) ) ) {
continue; }
2734 if ( ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(0))[0] >
static_cast< proshade_double
> ( bestFold ) ) || ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(1))[0] >
static_cast< proshade_double
> ( bestFold ) ) )
2739 if ( std::max ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(0))[6], CSym->at(settings->
allDetectedDAxes.at(dIt).at(1))[6] ) < bestHistFSCStart ) {
continue; }
2742 bestFold =
static_cast< proshade_unsign
> ( std::max ( CSym->at(settings->
allDetectedDAxes.at(dIt).at(0))[0], CSym->at(settings->
allDetectedDAxes.at(dIt).at(1))[0] ) );
2743 bestD =
static_cast< proshade_signed
> ( dIt );
2752 settings->
setRecommendedFold (
static_cast< proshade_unsign
> ( std::max ( CSym->at(settings->
allDetectedDAxes.at(
static_cast< size_t > ( bestD ) ).at(0))[0], CSym->at(settings->
allDetectedDAxes.at(
static_cast< size_t > ( bestD ) ).at(1))[0] ) ) );
2762 alreadyDecided =
true;
2767 if ( ( CSym->size() > 0 ) && !alreadyDecided )
2770 proshade_signed bestC = -1;
2771 proshade_unsign bestFold = 0;
2774 for (
size_t cIt = 0; cIt < CSym->size(); cIt++ )
2777 if ( cIt > 20 ) {
continue; }
2780 if ( CSym->at(cIt)[5] < bestHistPeakStart ) {
continue; }
2783 this->
computeFSC ( settings, CSym, cIt, mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2787 proshade_double bestHistFSCStart = this->findTopGroupSmooth ( CSym, 6, step, sigma, windowSize );
2790 for (
size_t cIt = 0; cIt < CSym->size(); cIt++ )
2793 if ( CSym->at(cIt)[0] >
static_cast< proshade_double
> ( bestFold ) )
2796 if ( ( CSym->at(cIt)[6] > settings->
fscThreshold ) && ( CSym->at(cIt)[6] >= bestHistFSCStart ) )
2798 bestFold =
static_cast< proshade_unsign
> ( CSym->at(cIt)[0] );
2799 bestC =
static_cast< proshade_signed
> ( cIt );
2809 settings->
setRecommendedFold (
static_cast< proshade_unsign
> ( CSym->at(
static_cast< size_t > ( bestC ) )[0] ) );
2814 alreadyDecided =
true;
2837 proshade_unsign bestIndex = 0;
2838 proshade_double highestSym = 0.0;
2841 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( CSym->size() ); iter++ )
2844 const FloatingPoint< proshade_double > lhs1 ( CSym->at(iter)[0] ), rhs1 (
static_cast< proshade_double
> ( settings->
requestedSymmetryFold ) );
2845 if ( !lhs1.AlmostEquals ( rhs1 ) ) {
continue; }
2848 if ( CSym->at(iter)[5] > highestSym )
2850 highestSym = CSym->at(iter)[5];
2856 if ( highestSym > 0.0 )
2859 settings->
setRecommendedFold (
static_cast< proshade_unsign
> ( CSym->at(bestIndex)[0] ) );
2886 void ProSHADE_internal_data::ProSHADE_data::saveRequestedSymmetryD (
ProSHADE_settings* settings, std::vector< proshade_double* >* DSym, std::vector< proshade_double* >* axes, fftw_complex* mapData, fftw_complex* origCoeffs, fftw_complex* fCoeffs, fftw_plan* planForwardFourier, proshade_signed noBins, proshade_signed* binIndexing, proshade_double** bindata, proshade_signed* binCounts )
2889 proshade_unsign bestIndex = 0;
2890 proshade_double highestSym = 0.0;
2893 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( DSym->size() ); iter++ )
2896 const FloatingPoint< proshade_double > lhs1 ( std::max ( DSym->at(iter)[0], DSym->at(iter)[7] ) ), rhs1 (
static_cast< proshade_double
> ( settings->
requestedSymmetryFold ) );
2897 if ( !lhs1.AlmostEquals ( rhs1 ) ) {
continue; }
2900 const FloatingPoint< proshade_double > lhs999a ( DSym->at(iter)[5] ), lhs999b ( DSym->at(iter)[12] ), rhs999 (
static_cast< proshade_double
> ( -999.9 ) );
2901 if ( ( DSym->at(iter)[5] < settings->
peakThresholdMin ) && !( lhs999a.AlmostEquals( rhs999 ) ) ) {
continue; }
2902 if ( ( DSym->at(iter)[12] < settings->
peakThresholdMin ) && !( lhs999b.AlmostEquals( rhs999 ) ) ) {
continue; }
2905 this->
computeFSC ( settings, &DSym->at(iter)[0], mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2906 this->
computeFSC ( settings, &DSym->at(iter)[7], mapData, fCoeffs, origCoeffs, planForwardFourier, noBins, binIndexing, bindata, binCounts );
2909 if ( ( DSym->at(iter)[6] + DSym->at(iter)[13] ) > highestSym )
2911 highestSym = ( DSym->at(iter)[6] + DSym->at(iter)[13] );
2917 if ( highestSym > 0.0 )
2920 settings->
setRecommendedFold (
static_cast< proshade_unsign
> ( std::max ( DSym->at(bestIndex)[0], DSym->at(bestIndex)[7] ) ) );
2952 std::vector< proshade_double > angList;
2953 std::vector<std::vector< proshade_double > > ret;
2956 proshade_double* rotMat =
new proshade_double[9];
2961 proshade_double normF = std::sqrt( std::pow ( xAx, 2.0 ) + std::pow ( yAx, 2.0 ) + std::pow ( zAx, 2.0 ) );
2967 if ( fold % 2 == 0 )
2970 for ( proshade_double iter =
static_cast < proshade_double
> ( -( ( fold / 2 ) - 1 ) ); iter <= static_cast < proshade_double > ( fold / 2 ); iter++ )
2978 for ( proshade_double iter =
static_cast < proshade_double
> ( -fold / 2 ); iter <= static_cast < proshade_double > ( fold / 2 ); iter++ )
2985 for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( angList.size() ); iter++ )
2991 std::vector < proshade_double > retEl;
2992 for ( proshade_unsign matIt = 0; matIt < 9; matIt++ )
3016 if ( obtainedAxes != requiredAxes )
3018 std::stringstream hlpSS;
3019 hlpSS <<
"The supplied number of axes for group element\n : detection ( >" << obtainedAxes <<
"< ) does not match the group type ( >" << groupType <<
"< ).";
3020 throw ProSHADE_exception (
"Mismatch between supplied number of axes and\n : symmetry type.",
"ES00059", __FILE__, __LINE__, __func__, hlpSS.str() );
3035 bool checkElementAlreadyExists ( std::vector<std::vector< proshade_double > >* elements, std::vector< proshade_double >* elem, proshade_double matrixTolerance )
3038 bool elementFound =
false;
3041 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( elements->size() ); elIt++ )
3045 elementFound =
true;
3051 return ( elementFound );
3064 bool isGroup =
true;
3067 for ( proshade_unsign gr1 = 0; gr1 < static_cast<proshade_unsign> ( elements->size() ); gr1++ )
3069 for ( proshade_unsign gr2 = 1; gr2 < static_cast<proshade_unsign> ( elements->size() ); gr2++ )
3072 if ( gr1 >= gr2 ) {
continue; }
3086 if ( !isGroup ) {
break; }
3105 std::vector< std::vector< proshade_double > > ret;
3108 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( first->size() ); elIt++ )
3117 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( second->size() ); elIt++ )
3128 for ( proshade_unsign gr1 = 0; gr1 < static_cast<proshade_unsign> ( first->size() ); gr1++ )
3130 for ( proshade_unsign gr2 = 0; gr2 < static_cast<proshade_unsign> ( second->size() ); gr2++ )
3171 std::vector<std::vector< proshade_double > > ret;
3174 if ( groupType ==
"C" )
3183 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(0)).at(0) ) );
3189 throw ProSHADE_exception (
"Computed point group elements do not form a group.",
"ES00060", __FILE__, __LINE__, __func__,
"The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3192 else if ( groupType ==
"D" )
3201 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(0)).at(0) ) );
3205 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(1)).at(0) ) );
3214 throw ProSHADE_exception (
"Computed point group elements do not form a group.",
"ES00060", __FILE__, __LINE__, __func__,
"The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3217 else if ( groupType ==
"T" )
3223 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3226 const FloatingPoint< proshade_double > lhs1 ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 3.0 );
3227 if ( lhs1.AlmostEquals ( rhs1 ) )
3233 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3241 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3244 const FloatingPoint< proshade_double > lhs1 ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 2.0 );
3245 if ( lhs1.AlmostEquals ( rhs1 ) )
3251 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3262 throw ProSHADE_exception (
"Computed point group elements do not form a group.",
"ES00060", __FILE__, __LINE__, __func__,
"The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3265 else if ( groupType ==
"O" )
3271 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3274 const FloatingPoint< proshade_double > lhs1 ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 4.0 );
3275 if ( lhs1.AlmostEquals ( rhs1 ) )
3281 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3289 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3292 const FloatingPoint< proshade_double > lhs1 ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 3.0 );
3293 if ( lhs1.AlmostEquals ( rhs1 ) )
3299 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3307 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3310 const FloatingPoint< proshade_double > lhs1 ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 2.0 );
3311 if ( lhs1.AlmostEquals ( rhs1 ) )
3317 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3328 throw ProSHADE_exception (
"Computed point group elements do not form a group.",
"ES00060", __FILE__, __LINE__, __func__,
"The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3331 else if ( groupType ==
"I" )
3337 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3340 const FloatingPoint< proshade_double > lhs1 ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 5.0 );
3341 if ( lhs1.AlmostEquals ( rhs1 ) )
3347 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3355 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3358 const FloatingPoint< proshade_double > lhs1 ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 3.0 );
3359 if ( lhs1.AlmostEquals ( rhs1 ) )
3365 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3373 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3376 const FloatingPoint< proshade_double > lhs1 ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ), rhs1 ( 2.0 );
3377 if ( lhs1.AlmostEquals ( rhs1 ) )
3383 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3394 throw ProSHADE_exception (
"Computed point group elements do not form a group.",
"ES00060", __FILE__, __LINE__, __func__,
"The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3397 else if ( groupType ==
"X" )
3400 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
3406 static_cast< proshade_signed
> ( settings->
allDetectedCAxes.at(axesList.at(grIt)).at(0) ) );
3416 throw ProSHADE_exception (
"Computed point group elements do not form a group.",
"ES00060", __FILE__, __LINE__, __func__,
"The supplied cyclic groups list does not form a group and\n : therefore such group's elements cannot be obtained. Please\n : check the cyclic groups list supplied to the\n : getAllGroupElements() function." );
3421 std::stringstream hlpSS;
3422 hlpSS <<
"Unknown symmetry type: >" << groupType <<
"<";
3423 throw ProSHADE_exception ( hlpSS.str().c_str(),
"ES00058", __FILE__, __LINE__, __func__,
"Function getAllGroupElements was called with symmetry type\n : value outside of the allowed values C, D, T, O, I\n : or empty for using all supplied axes." );
3438 if ( saveTo !=
nullptr )
3445 saveTo =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
3451 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
3453 saveTo[iter] = this->internalMap[iter];
3476 std::stringstream ssHlp;
3477 std::vector< proshade_double > comMove = this->getMapCOMProcessChange ( );
3478 ssHlp << std::endl <<
"Detected " << settings->
recommendedSymmetryType <<
" symmetry with fold " << settings->
recommendedSymmetryFold <<
" about point [" << comMove.at(0) <<
" , " << comMove.at(1) <<
" , " << comMove.at(2) <<
"] away from centre of mass .";
3483 ssHlp.clear(); ssHlp.str (
"" );
3484 ssHlp <<
" Fold X Y Z Angle Height Average FSC";
3487 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->
detectedSymmetry.size() ); symIt++ )
3489 ssHlp.clear(); ssHlp.str (
"" );
3494 std::stringstream hlpSS3;
3495 ssHlp.clear(); ssHlp.str (
"" );
3496 hlpSS3 << std::endl <<
"To facilitate manual checking for symmetries, the following is a list of all detected C symmetries:";
3501 ssHlp.clear(); ssHlp.str (
"" );
3502 ssHlp <<
" Fold X Y Z Angle Height Average FSC";
3505 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->
allDetectedCAxes.size() ); symIt++ )
3507 ssHlp.clear(); ssHlp.str (
"" );
3529 proshade_double totNonZeroPoints = 0.0;
3530 proshade_signed mapIt = 0;
3533 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
3535 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
3537 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
3540 mapIt = zIt +
static_cast< proshade_signed
> ( this->zDimIndices ) * ( yIt +
static_cast< proshade_signed
> ( this->yDimIndices ) * xIt );
3543 if ( this->internalMap[mapIt] <= 0.0 ) {
continue; }
3546 this->xCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( xIt + this->xFrom );
3547 this->yCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( yIt + this->yFrom );
3548 this->zCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( zIt + this->zFrom );
3549 totNonZeroPoints += this->internalMap[mapIt];
3554 this->xCom /= totNonZeroPoints;
3555 this->yCom /= totNonZeroPoints;
3556 this->zCom /= totNonZeroPoints;
3559 this->xCom =
static_cast< proshade_double
> ( (
static_cast< proshade_single
> ( this->xFrom ) * ( this->xDimSizeOriginal /
static_cast< proshade_single
> ( this->xDimIndicesOriginal ) ) ) +
3560 ( (
static_cast< proshade_single
> ( this->xCom ) -
static_cast< proshade_single
> ( this->xFrom ) ) *
3561 (
static_cast< proshade_single
> ( this->xDimSizeOriginal ) /
static_cast< proshade_single
> ( this->xDimIndicesOriginal ) ) ) );
3562 this->yCom =
static_cast< proshade_double
> ( (
static_cast< proshade_single
> ( this->yFrom ) * ( this->yDimSizeOriginal /
static_cast< proshade_single
> ( this->yDimIndicesOriginal ) ) ) +
3563 ( (
static_cast< proshade_single
> ( this->yCom ) -
static_cast< proshade_single
> ( this->yFrom ) ) *
3564 (
static_cast< proshade_single
> ( this->yDimSizeOriginal ) /
static_cast< proshade_single
> ( this->yDimIndicesOriginal ) ) ) );
3565 this->zCom =
static_cast< proshade_double
> ( (
static_cast< proshade_single
> ( this->zFrom ) * ( this->zDimSizeOriginal /
static_cast< proshade_single
> ( this->zDimIndicesOriginal ) ) ) +
3566 ( (
static_cast< proshade_single
> ( this->zCom ) -
static_cast< proshade_single
> ( this->zFrom ) ) *
3567 (
static_cast< proshade_single
> ( this->zDimSizeOriginal ) /
static_cast< proshade_single
> ( this->zDimIndicesOriginal ) ) ) );
3581 return ( this->noSpheres );
3592 return ( this->internalMap[pos] );
3602 return ( this->maxShellBand );
3612 return ( this->rrpMatrices[band][sh1][sh2] );
3627 if ( this->spheres[shell]->getLocalBandwidth( ) >= bandVal )
3651 fftw_complex* mapCoeffs =
new fftw_complex[this->xDimIndices * this->yDimIndices * this->zDimIndices];
3652 fftw_complex* pattersonMap =
new fftw_complex[this->xDimIndices * this->yDimIndices * this->zDimIndices];
3659 for ( proshade_unsign iter = 0; iter < (this->xDimIndices * this->yDimIndices * this->zDimIndices); iter++ )
3661 pattersonMap[iter][0] = this->internalMap[iter];
3662 pattersonMap[iter][1] = 0.0;
3666 fftw_plan forward = fftw_plan_dft_3d (
static_cast< int > ( this->xDimIndices ),
static_cast< int > ( this->yDimIndices ),
static_cast< int > ( this->zDimIndices ),
3667 pattersonMap, mapCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
3668 fftw_plan inverse = fftw_plan_dft_3d (
static_cast< int > ( this->xDimIndices ),
static_cast< int > ( this->yDimIndices ),
static_cast< int > ( this->zDimIndices ),
3669 mapCoeffs, pattersonMap, FFTW_BACKWARD, FFTW_ESTIMATE );
3672 fftw_execute ( forward );
3678 fftw_execute ( inverse );
3681 proshade_signed mapIt, patIt, patX, patY, patZ;
3682 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
3684 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
3686 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
3689 patX = xIt - (
static_cast<proshade_signed
> ( this->xDimIndices ) / 2 );
if ( patX < 0 ) { patX += this->xDimIndices; }
3690 patY = yIt - (
static_cast<proshade_signed
> ( this->yDimIndices ) / 2 );
if ( patY < 0 ) { patY += this->yDimIndices; }
3691 patZ = zIt - (
static_cast<proshade_signed
> ( this->zDimIndices ) / 2 );
if ( patZ < 0 ) { patZ += this->zDimIndices; }
3694 mapIt = zIt +
static_cast< proshade_signed
> ( this->zDimIndices ) * ( yIt +
static_cast< proshade_signed
> ( this->yDimIndices ) * xIt );
3695 patIt = patZ +
static_cast< proshade_signed
> ( this->zDimIndices ) * ( patY +
static_cast< proshade_signed
> ( this->yDimIndices ) * patX );
3698 this->internalMap[mapIt] = pattersonMap[patIt][0];
3704 delete[] pattersonMap;
3708 fftw_destroy_plan ( forward );
3709 fftw_destroy_plan ( inverse );
3726 return ( &this->sphericalHarmonics[shell][seanindex (
static_cast< int > ( order ) -
static_cast< int > ( band ),
3727 static_cast< int > ( band ),
3728 static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ) )][0] );
3739 return ( &this->sphericalHarmonics[shell][seanindex (
static_cast< int > ( order ) -
static_cast< int > ( band ),
3740 static_cast< int > ( band ),
3741 static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ) )][1] );
3752 return ( this->spheres[shell]->getShellRadius() );
3763 return ( this->integrationWeight );
3775 return ( this->spheres[shell]->getLocalBandwidth ( ) );
3787 return ( this->spherePos.at(shell) );
3799 return ( this->eMatrices[band] );
3814 *valueReal = this->eMatrices[band][order1][order2][0];
3815 *valueImag = this->eMatrices[band][order1][order2][1];
3829 return ( this->so3CoeffsInverse );
3840 return ( this->so3Coeffs );
3851 return ( this->maxCompBand );
3866 *valueReal = this->wignerMatrices[band][order1][order2][0];
3867 *valueImag = this->wignerMatrices[band][order1][order2][1];
3881 return ( this->xDimSize );
3891 return ( this->yDimSize );
3901 return ( this->zDimSize );
3911 return ( this->xDimIndices );
3921 return ( this->yDimIndices );
3931 return ( this->zDimIndices );
3941 return ( &this->xFrom );
3951 return ( &this->yFrom );
3961 return ( &this->zFrom );
3971 return ( &this->xTo );
3981 return ( &this->yTo );
3991 return ( &this->zTo );
4001 return ( &this->xAxisOrigin );
4011 return ( &this->yAxisOrigin );
4021 return ( &this->zAxisOrigin );
4031 return ( this->internalMap );
4041 return ( this->translationMap );
4051 std::vector< proshade_double > ret;
4069 this->integrationWeight = intW;
4083 this->integrationWeight += intW;
4100 this->eMatrices[band][order1][order2][0] = val[0];
4101 this->eMatrices[band][order1][order2][1] = val[1];
4118 this->eMatrices[band][order1][order2][0] /= normF;
4119 this->eMatrices[band][order1][order2][1] /= normF;
4134 this->so3Coeffs[position][0] = val[0];
4135 this->so3Coeffs[position][1] = val[1];
4152 this->wignerMatrices[band][order1][order2][0] = val[0];
4153 this->wignerMatrices[band][order1][order2][1] = val[1];
4170 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4172 eMatsLMReal[iter] =
static_cast<double> ( this->eMatrices[band][order1][iter][0] );
4190 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4192 eMatsLMImag[iter] =
static_cast<double> ( this->eMatrices[band][order1][iter][1] );
4208 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4210 so3CoefsReal[iter] =
static_cast<double> ( this->so3Coeffs[iter][0] );
4226 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4228 so3CoefsImag[iter] =
static_cast<double> ( this->so3Coeffs[iter][1] );
4248 return (
static_cast<int> ( so3CoefLoc (
static_cast< int > ( order1 ),
static_cast< int > ( order2 ),
static_cast< int > ( band ),
static_cast< int > ( this->getMaxBand() ) ) ) );
4259 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4261 rotFunReal[iter] =
static_cast<double> ( this->so3CoeffsInverse[iter][0] );
4277 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4279 rotFunImag[iter] =
static_cast<double> ( this->so3CoeffsInverse[iter][1] );
4295 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4297 trsFunReal[iter] =
static_cast<double> ( this->translationMap[iter][0] );
4313 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4315 trsFunImag[iter] =
static_cast<double> ( this->translationMap[iter][1] );
4334 proshade_double eA, eB, eG;
4338 proshade_double* rMat =
nullptr;
4339 rMat =
new proshade_double[9];
4346 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
4348 rotMat[iter] =
static_cast<double> ( rMat[iter] );
4389 return (
static_cast<proshade_unsign
> ( settings->
detectedSymmetry.size() ) );
4401 if (
static_cast<proshade_unsign
> ( settings->
detectedSymmetry.size() ) <= axisNo )
4404 return ( std::vector< std::string > ( ) );
4408 std::vector< std::string > ret;
4411 std::stringstream ssHlp;
4460 std::stringstream fNameHlp;
4462 this->writeMap ( fNameHlp.str() );
4469 this->writePdb ( fNameHlp.str(), eulA, eulB, eulG, ultimateTranslation->at(0), ultimateTranslation->at(1), ultimateTranslation->at(2), rotCentre->at(0), rotCentre->at(1), rotCentre->at(2), settings->
firstModelOnly );
4475 ultimateTranslation->at(0), ultimateTranslation->at(1), ultimateTranslation->at(2),
4497 std::stringstream rotCen; rotCen << std::setprecision (3) << std::showpos <<
"The rotation centre to origin translation vector is: " << -rotationCentre->at(0) <<
" " << -rotationCentre->at(1) <<
" " << -rotationCentre->at(2);
4501 proshade_double* rotMat =
new proshade_double[9];
4505 std::stringstream rotMatSS;
4506 rotMatSS << std::setprecision (3) << std::showpos <<
"The rotation matrix about origin is : " << rotMat[0] <<
" " << rotMat[1] <<
" " << rotMat[2] << std::endl;
4507 rotMatSS << std::setprecision (3) << std::showpos <<
" : " << rotMat[3] <<
" " << rotMat[4] <<
" " << rotMat[5] << std::endl;
4508 rotMatSS << std::setprecision (3) << std::showpos <<
" : " << rotMat[6] <<
" " << rotMat[7] <<
" " << rotMat[8];
4514 std::stringstream finTrs; finTrs << std::setprecision (3) << std::showpos <<
"The rotation centre to overlay translation vector is: " << finalTranslation->at(0) <<
" " << finalTranslation->at(1) <<
" " << finalTranslation->at(2);