28 #define GEMMI_WRITE_IMPLEMENTATION
29 #include <gemmi/to_pdb.hpp>
44 this->
fileType = ProSHADE_internal_io::UNKNOWN;
109 this->
spherePos = std::vector<proshade_single> ( );
160 ProSHADE_internal_data::ProSHADE_data::ProSHADE_data (
ProSHADE_settings* settings, 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 )
164 this->fileName = strName;
165 this->fileType = ProSHADE_internal_io::MAP;
168 this->internalMap = NULL;
171 this->xDimSize = xDmSz;
172 this->yDimSize = yDmSz;
173 this->zDimSize = zDmSz;
177 this->xDimIndices = xDmInd;
178 this->yDimIndices = yDmInd;
179 this->zDimIndices = zDmInd;
180 this->xGridIndices = xDmInd;
181 this->yGridIndices = yDmInd;
182 this->zGridIndices = zDmInd;
183 this->xAxisOrder = 1;
184 this->yAxisOrder = 2;
185 this->zAxisOrder = 3;
186 this->xAxisOrigin = xFr;
187 this->yAxisOrigin = yFr;
188 this->zAxisOrigin = zFr;
194 this->xDimSizeOriginal = 0.0;
195 this->yDimSizeOriginal = 0.0;
196 this->zDimSizeOriginal = 0.0;
197 this->xDimIndicesOriginal = 0;
198 this->yDimIndicesOriginal = 0;
199 this->zDimIndicesOriginal = 0;
200 this->xAxisOriginOriginal = 0;
201 this->yAxisOriginOriginal = 0;
202 this->zAxisOriginOriginal = 0;
203 this->originalMapXCom = 0.0;
204 this->originalMapYCom = 0.0;
205 this->originalMapZCom = 0.0;
206 this->mapMovFromsChangeX = 0.0;
207 this->mapMovFromsChangeY = 0.0;
208 this->mapMovFromsChangeZ = 0.0;
209 this->mapCOMProcessChangeX = 0.0;
210 this->mapCOMProcessChangeY = 0.0;
211 this->mapCOMProcessChangeZ = 0.0;
214 this->originalPdbRotCenX = 0.0;
215 this->originalPdbRotCenY = 0.0;
216 this->originalPdbRotCenZ = 0.0;
217 this->originalPdbTransX = 0.0;
218 this->originalPdbTransY = 0.0;
219 this->originalPdbTransZ = 0.0;
230 this->spherePos = std::vector<proshade_single> ( );
232 this->spheres = NULL;
233 this->sphericalHarmonics = NULL;
234 this->rotSphericalHarmonics = NULL;
235 this->maxShellBand = 0;
238 this->rrpMatrices = NULL;
239 this->eMatrices = NULL;
240 this->so3Coeffs = NULL;
241 this->so3CoeffsInverse = NULL;
242 this->wignerMatrices = NULL;
243 this->integrationWeight = 0.0;
244 this->maxCompBand = 0;
245 this->translationMap = NULL;
248 this->isEmpty =
false;
249 this->inputOrder = inputO;
252 if (
static_cast<proshade_unsign
> ( len ) != ( xDmInd * yDmInd * zDmInd ) )
254 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." );
257 if ( (
static_cast<proshade_signed
> ( xT - xFr ) !=
static_cast<proshade_signed
> ( xDmInd - 1 ) ) ||
258 (
static_cast<proshade_signed
> ( yT - yFr ) !=
static_cast<proshade_signed
> ( yDmInd - 1 ) ) ||
259 (
static_cast<proshade_signed
> ( zT - zFr ) !=
static_cast<proshade_signed
> ( zDmInd - 1 ) ) )
261 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." );
265 this->internalMap =
new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
269 proshade_unsign arrPos = 0;
270 for ( proshade_unsign xIt = 0; xIt < this->xDimIndices; xIt++ )
272 for ( proshade_unsign yIt = 0; yIt < this->yDimIndices; yIt++ )
274 for ( proshade_unsign zIt = 0; zIt < this->zDimIndices; zIt++ )
276 arrPos = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
277 this->internalMap[arrPos] =
static_cast<proshade_double
> ( mapVals[arrPos] );
298 if ( this->internalMap != NULL )
300 delete[] this->internalMap;
304 if ( this->spheres != NULL )
306 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
308 if ( this->spheres[iter] != NULL )
310 delete this->spheres[iter];
311 this->spheres[iter] = NULL;
314 delete[] this->spheres;
318 if ( this->sphericalHarmonics != NULL )
320 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
322 if ( this->sphericalHarmonics[iter] != NULL )
324 delete[] this->sphericalHarmonics[iter];
325 this->sphericalHarmonics[iter] = NULL;
328 delete[] this->sphericalHarmonics;
332 if ( this->rotSphericalHarmonics != NULL )
334 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
336 if ( this->rotSphericalHarmonics[iter] != NULL )
338 delete[] this->rotSphericalHarmonics[iter];
339 this->rotSphericalHarmonics[iter] = NULL;
342 delete[] this->rotSphericalHarmonics;
346 if ( this->rrpMatrices != NULL )
348 for ( proshade_unsign bwIt = 0; bwIt < this->maxShellBand; bwIt++ )
350 if ( this->rrpMatrices[bwIt] != NULL )
352 for ( proshade_unsign shIt = 0; shIt < this->noSpheres; shIt++ )
354 if ( this->rrpMatrices[bwIt][shIt] != NULL )
356 delete[] this->rrpMatrices[bwIt][shIt];
360 delete[] this->rrpMatrices[bwIt];
364 delete[] this->rrpMatrices;
368 if ( this->eMatrices != NULL )
370 for ( proshade_unsign bandIter = 0; bandIter < this->maxCompBand; bandIter++ )
372 if ( this->eMatrices[bandIter] != NULL )
374 for ( proshade_unsign band2Iter = 0; band2Iter < static_cast<proshade_unsign> ( ( bandIter * 2 ) + 1 ); band2Iter++ )
376 if ( this->eMatrices[bandIter][band2Iter] != NULL )
378 delete[] this->eMatrices[bandIter][band2Iter];
382 delete[] this->eMatrices[bandIter];
386 delete[] this->eMatrices;
390 if ( this->so3Coeffs != NULL )
392 delete[] this->so3Coeffs;
394 if ( this->so3CoeffsInverse != NULL )
396 delete[] this->so3CoeffsInverse;
400 if ( this->wignerMatrices != NULL )
402 for ( proshade_unsign bandIter = 1; bandIter < this->maxCompBand; bandIter++ )
404 if ( this->wignerMatrices[bandIter] != NULL )
406 for ( proshade_unsign order1Iter = 0; order1Iter < ( (bandIter * 2) + 1 ); order1Iter++ )
408 if ( this->wignerMatrices[bandIter][order1Iter] != NULL )
410 delete[] this->wignerMatrices[bandIter][order1Iter];
413 delete[] this->wignerMatrices[bandIter];
416 delete[] wignerMatrices;
420 if ( this->translationMap != NULL )
422 delete[] this->translationMap;
426 if ( this->sphereMappedRotFun.size() > 0 )
428 for ( proshade_unsign spIt = 0; spIt < static_cast<proshade_unsign> ( this->sphereMappedRotFun.size() ); spIt++ )
430 delete this->sphereMappedRotFun.at(spIt);
453 if ( !this->isEmpty )
455 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 );
459 this->fileName = fName;
465 this->inputOrder = inputO;
468 switch ( this->fileType )
470 case ProSHADE_internal_io::UNKNOWN:
471 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." );
474 case ProSHADE_internal_io::PDB:
475 this->readInPDB ( settings );
478 case ProSHADE_internal_io::MAP:
479 this->readInMAP ( settings );
484 this->isEmpty =
false;
504 gemmi::Ccp4<float> map;
505 map.read_ccp4 ( gemmi::MaybeGzipped ( this->fileName.c_str() ) );
508 map.setup ( gemmi::GridSetup::ReorderOnly, NAN );
512 &this->xDimIndices, &this->yDimIndices, &this->zDimIndices,
513 &this->xDimSize, &this->yDimSize, &this->zDimSize,
514 &this->aAngle, &this->bAngle, &this->cAngle,
515 &this->xFrom, &this->yFrom, &this->zFrom,
516 &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin,
517 &this->xAxisOrder, &this->yAxisOrder, &this->zAxisOrder,
518 &this->xGridIndices, &this->yGridIndices, &this->zGridIndices );
521 ProSHADE_internal_io::readInMapData ( &map, this->internalMap, this->xDimIndices, this->yDimIndices, this->zDimIndices, this->xAxisOrder, this->yAxisOrder, this->zAxisOrder );
526 settings->
setResolution ( std::min (
static_cast<proshade_double
> ( this->xDimSize ) /
static_cast<proshade_double
> ( this->xDimIndices ),
527 std::min (
static_cast<proshade_double
> ( this->yDimSize ) /
static_cast<proshade_double
> ( this->yDimIndices ),
528 static_cast<proshade_double
> ( this->zDimSize ) /
static_cast<proshade_double
> ( this->zDimIndices ) ) ) * 2.0 );
532 this->figureIndexStartStop ( );
538 proshade_double xMapCOMPreReSampl = 0.0, yMapCOMPreReSampl = 0.0, zMapCOMPreReSampl = 0.0;
540 this->xDimSize, this->yDimSize, this->zDimSize,
541 this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
544 proshade_double xOrPre = this->xFrom * ( this->xDimSize /
static_cast<proshade_double
> ( this->xDimIndices ) );
545 proshade_double yOrPre = this->yFrom * ( this->yDimSize /
static_cast<proshade_double
> ( this->yDimIndices ) );
546 proshade_double zOrPre = this->zFrom * ( this->zDimSize /
static_cast<proshade_double
> ( this->zDimIndices ) );
549 this->reSampleMap ( settings );
552 proshade_double xMapCOMPostReSampl = 0.0, yMapCOMPostReSampl = 0.0, zMapCOMPostReSampl = 0.0;
554 this->xDimSize, this->yDimSize, this->zDimSize,
555 this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
558 proshade_double xOrPst = this->xFrom * ( this->xDimSize /
static_cast<proshade_double
> ( this->xDimIndices ) );
559 proshade_double yOrPst = this->yFrom * ( this->yDimSize /
static_cast<proshade_double
> ( this->yDimIndices ) );
560 proshade_double zOrPst = this->zFrom * ( this->zDimSize /
static_cast<proshade_double
> ( this->zDimIndices ) );
563 proshade_double xFinMov = ( xMapCOMPreReSampl - xMapCOMPostReSampl ) - ( xOrPre - xOrPst );
564 proshade_double yFinMov = ( yMapCOMPreReSampl - yMapCOMPostReSampl ) - ( yOrPre - yOrPst );
565 proshade_double zFinMov = ( zMapCOMPreReSampl - zMapCOMPostReSampl ) - ( zOrPre - zOrPst );
569 this->xDimSize, this->yDimSize, this->zDimSize,
570 this->xDimIndices, this->yDimIndices, this->zDimIndices );
574 this->xDimSizeOriginal = this->xDimSize;
575 this->yDimSizeOriginal = this->yDimSize;
576 this->zDimSizeOriginal = this->zDimSize;
579 this->xDimIndicesOriginal = this->xDimIndices;
580 this->yDimIndicesOriginal = this->yDimIndices;
581 this->zDimIndicesOriginal = this->zDimIndices;
584 this->xAxisOriginOriginal = this->xAxisOrigin;
585 this->yAxisOriginOriginal = this->yAxisOrigin;
586 this->zAxisOriginOriginal = this->zAxisOrigin;
589 this->findMapCOM ( );
590 this->originalMapXCom = this->xCom;
591 this->originalMapYCom = this->yCom;
592 this->originalMapZCom = this->zCom;
616 gemmi::Structure pdbFile = gemmi::read_structure ( gemmi::MaybeGzipped ( this->fileName ) );
631 proshade_double xCOMPdb, yCOMPdb, zCOMPdb;
635 proshade_single xF, xT, yF, yT, zF, zT;
639 proshade_single xMov = 20.0 - xF;
640 proshade_single yMov = 20.0 - yF;
641 proshade_single zMov = 20.0 - zF;
645 this->xDimSize = xT - xF + 40.0;
646 this->yDimSize = yT - yF + 40.0;
647 this->zDimSize = zT - zF + 40.0;
650 ProSHADE_internal_mapManip::generateMapFromPDB ( pdbFile, this->internalMap, settings->
requestedResolution, this->xDimSize, this->yDimSize, this->zDimSize, &this->xTo, &this->yTo, &this->zTo, settings->
forceP1, settings->
firstModelOnly );
653 this->setPDBMapValues ( );
656 proshade_double xCOMMap, yCOMMap, zCOMMap;
658 this->xDimSize, this->yDimSize, this->zDimSize,
659 this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
661 if ( pdbFile.models.size() > 1 )
663 xMov = xCOMMap - xCOMPdb;
664 yMov = yCOMMap - yCOMPdb;
665 zMov = zCOMMap - zCOMPdb;
670 &this->xFrom, &this->xTo, &this->yFrom, &this->yTo, &this->zFrom, &this->zTo,
671 &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
673 this->xDimIndices, this->yDimIndices, this->zDimIndices );
679 proshade_double xMapCOMPreReSampl = 0.0, yMapCOMPreReSampl = 0.0, zMapCOMPreReSampl = 0.0;
681 this->xDimSize, this->yDimSize, this->zDimSize,
682 this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
685 proshade_double xOrPre = this->xFrom * ( this->xDimSize /
static_cast<proshade_double
> ( this->xDimIndices ) );
686 proshade_double yOrPre = this->yFrom * ( this->yDimSize /
static_cast<proshade_double
> ( this->yDimIndices ) );
687 proshade_double zOrPre = this->zFrom * ( this->zDimSize /
static_cast<proshade_double
> ( this->zDimIndices ) );
690 this->reSampleMap ( settings );
693 proshade_double xMapCOMPostReSampl = 0.0, yMapCOMPostReSampl = 0.0, zMapCOMPostReSampl = 0.0;
695 this->xDimSize, this->yDimSize, this->zDimSize,
696 this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
699 proshade_double xOrPst = this->xFrom * ( this->xDimSize /
static_cast<proshade_double
> ( this->xDimIndices ) );
700 proshade_double yOrPst = this->yFrom * ( this->yDimSize /
static_cast<proshade_double
> ( this->yDimIndices ) );
701 proshade_double zOrPst = this->zFrom * ( this->zDimSize /
static_cast<proshade_double
> ( this->zDimIndices ) );
704 proshade_double xFinMov = ( xMapCOMPreReSampl - xMapCOMPostReSampl ) - ( xOrPre - xOrPst );
705 proshade_double yFinMov = ( yMapCOMPreReSampl - yMapCOMPostReSampl ) - ( yOrPre - yOrPst );
706 proshade_double zFinMov = ( zMapCOMPreReSampl - zMapCOMPostReSampl ) - ( zOrPre - zOrPst );
710 this->xDimSize, this->yDimSize, this->zDimSize,
711 this->xDimIndices, this->yDimIndices, this->zDimIndices );
715 this->xDimSizeOriginal = this->xDimSize;
716 this->yDimSizeOriginal = this->yDimSize;
717 this->zDimSizeOriginal = this->zDimSize;
720 this->xDimIndicesOriginal = this->xDimIndices;
721 this->yDimIndicesOriginal = this->yDimIndices;
722 this->zDimIndicesOriginal = this->zDimIndices;
725 this->xAxisOriginOriginal = this->xAxisOrigin;
726 this->yAxisOriginOriginal = this->yAxisOrigin;
727 this->zAxisOriginOriginal = this->zAxisOrigin;
730 this->findMapCOM ( );
731 this->originalMapXCom = this->xCom;
732 this->originalMapYCom = this->yCom;
733 this->originalMapZCom = this->zCom;
757 this->xDimIndices = this->xTo;
758 this->yDimIndices = this->yTo;
759 this->zDimIndices = this->zTo;
762 this->xGridIndices = this->xDimIndices;
763 this->yGridIndices = this->yDimIndices;
764 this->zGridIndices = this->zDimIndices;
767 this->xAxisOrder = 1;
768 this->yAxisOrder = 2;
769 this->zAxisOrder = 3;
772 this->xAxisOrigin = this->xFrom;
773 this->yAxisOrigin = this->yFrom;
774 this->zAxisOrigin = this->zFrom;
788 this->xTo = this->xFrom + this->xDimIndices - 1;
789 this->yTo = this->yFrom + this->yDimIndices - 1;
790 this->zTo = this->zFrom + this->zDimIndices - 1;
810 gemmi::Grid<float> mapData;
811 mapData.set_unit_cell ( this->xDimSize, this->yDimSize, this->zDimSize, this->aAngle, this->bAngle, this->cAngle );
812 mapData.set_size_without_checking ( this->xDimIndices, this->yDimIndices, this->zDimIndices );
813 mapData.axis_order = gemmi::AxisOrder::XYZ;
814 mapData.spacegroup = &gemmi::get_spacegroup_p1();
817 gemmi::Ccp4<float> map;
819 map.update_ccp4_header ( mode );
823 this->xDimIndices, this->yDimIndices, this->zDimIndices,
824 this->xDimSize, this->yDimSize, this->zDimSize,
825 this->aAngle, this->bAngle, this->cAngle,
826 this->xFrom, this->yFrom, this->zFrom,
827 this->xAxisOrigin, this->yAxisOrigin, this->zAxisOrigin,
828 this->xAxisOrder, this->yAxisOrder, this->zAxisOrder,
829 this->xGridIndices, this->yGridIndices, this->zGridIndices,
833 proshade_unsign arrPos = 0;
834 for ( proshade_unsign uIt = 0; uIt < this->xDimIndices; uIt++ )
836 for ( proshade_unsign vIt = 0; vIt < this->yDimIndices; vIt++ )
838 for ( proshade_unsign wIt = 0; wIt < this->zDimIndices; wIt++ )
840 arrPos = wIt + this->zDimIndices * ( vIt + this->yDimIndices * uIt );
841 map.grid.set_value ( uIt, vIt, wIt,
static_cast<float> ( this->internalMap[arrPos] ) );
847 map.update_ccp4_header ( mode,
true );
850 map.write_ccp4_map ( fName );
877 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." );
881 gemmi::Structure pdbFile = gemmi::read_structure ( gemmi::MaybeGzipped ( this->fileName ) );
884 if ( ( euA != 0.0 ) || ( euB != 0.0 ) || ( euG != 0.0 ) )
894 std::ofstream outCoOrdFile;
895 outCoOrdFile.open ( fName.c_str() );
897 if ( outCoOrdFile.is_open() )
899 gemmi::PdbWriteOptions opt;
900 write_pdb ( pdbFile, outCoOrdFile, opt );
904 std::stringstream hlpMessage;
905 hlpMessage <<
"Failed to open the PDB file " << fName <<
" for output.";
906 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." );
909 outCoOrdFile.close ( );
926 proshade_double* hlpMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
930 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
932 hlpMap[iter] = this->internalMap[iter];
933 this->internalMap[iter] = mask[iter];
937 this->writeMap ( fName );
940 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
942 this->internalMap[iter] = hlpMap[iter];
966 proshade_signed arrayPos, invPos;
969 proshade_double* hlpMap =
new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
973 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
975 hlpMap[iter] = this->internalMap[iter];
979 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
981 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
983 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
986 arrayPos = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
987 invPos = ( (this->zDimIndices-1) - zIt ) + this->zDimIndices * ( ( (this->yDimIndices-1) - yIt ) + this->yDimIndices * ( (this->xDimIndices-1) - xIt ) );
990 this->internalMap[invPos] = hlpMap[arrayPos];
1020 std::vector<proshade_double> mapVals ( this->xDimIndices * this->yDimIndices * this->zDimIndices, 0.0 );
1023 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1025 mapVals.at(iter) = this->internalMap[iter];
1029 proshade_double* meanSD =
new proshade_double[2];
1033 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1035 this->internalMap[iter] = ( this->internalMap[iter] - meanSD[0] ) / meanSD[1];
1067 proshade_double* blurredMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1072 this->xDimSize, this->yDimSize, this->zDimSize, settings->
blurFactor );
1078 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 ); } }
1081 delete[] blurredMap;
1109 for ( proshade_unsign iter = 0; iter < 6; iter++ ) { ret[iter] = settings->
forceBounds[iter]; }
1116 this->xDimSize, this->yDimSize, this->zDimSize, ret );
1120 this->xDimSize, this->yDimSize, this->zDimSize, ret, settings->
boundsExtraSpace );
1126 std::stringstream ssHlp;
1127 ssHlp <<
"New boundaries are: " << ret[1] - ret[0] + 1 <<
" x " << ret[3] - ret[2] + 1 <<
" x " << ret[5] - ret[4] + 1;
1133 for ( proshade_unsign iter = 0; iter < 6; iter++ ) { settings->
forceBounds[iter] = ret[iter]; }
1163 newStr->
fileType = ProSHADE_internal_io::MAP;
1166 newStr->
xDimIndices =
static_cast<proshade_signed
> ( newBounds[1] ) -
static_cast<proshade_signed
> ( newBounds[0] ) + 1;
1167 newStr->
yDimIndices =
static_cast<proshade_signed
> ( newBounds[3] ) -
static_cast<proshade_signed
> ( newBounds[2] ) + 1;
1168 newStr->
zDimIndices =
static_cast<proshade_signed
> ( newBounds[5] ) -
static_cast<proshade_signed
> ( newBounds[4] ) + 1;
1170 newStr->
aAngle = this->aAngle;
1171 newStr->
bAngle = this->aAngle;
1172 newStr->
cAngle = this->aAngle;
1174 newStr->
xDimSize =
static_cast<proshade_single
> ( newStr->
xDimIndices ) * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices ) );
1175 newStr->
yDimSize =
static_cast<proshade_single
> ( newStr->
yDimIndices ) * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices ) );
1176 newStr->
zDimSize =
static_cast<proshade_single
> ( newStr->
zDimIndices ) * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices ) );
1186 newStr->
xAxisOrigin = this->xAxisOrigin + newBounds[0];
1187 newStr->
yAxisOrigin = this->yAxisOrigin + newBounds[2];
1188 newStr->
zAxisOrigin = this->zAxisOrigin + newBounds[4];
1190 newStr->
xFrom = this->xFrom + newBounds[0];
1191 newStr->
yFrom = this->yFrom + newBounds[2];
1192 newStr->
zFrom = this->zFrom + newBounds[4];
1194 newStr->
xTo = this->xTo - ( (this->xDimIndices-1) - newBounds[1] );
1195 newStr->
yTo = this->yTo - ( (this->yDimIndices-1) - newBounds[3] );
1196 newStr->
zTo = this->zTo - ( (this->zDimIndices-1) - newBounds[5] );
1205 this->xDimIndices, this->yDimIndices, this->zDimIndices, newStr->
internalMap, this->internalMap );
1225 proshade_single* changeVals =
new proshade_single[6];
1231 this->xDimSize, this->yDimSize, this->zDimSize, changeVals );
1241 this->xDimSize, this->yDimSize, this->zDimSize, changeVals );
1246 this->xDimIndices +=
static_cast<proshade_unsign
> ( changeVals[0] );
1247 this->yDimIndices +=
static_cast<proshade_unsign
> ( changeVals[1] );
1248 this->zDimIndices +=
static_cast<proshade_unsign
> ( changeVals[2] );
1250 this->xGridIndices = this->xDimIndices;
1251 this->yGridIndices = this->yDimIndices;
1252 this->zGridIndices = this->zDimIndices;
1254 this->xTo +=
static_cast<proshade_unsign
> ( changeVals[0] );
1255 this->yTo +=
static_cast<proshade_unsign
> ( changeVals[1] );
1256 this->zTo +=
static_cast<proshade_unsign
> ( changeVals[2] );
1258 this->xDimSize = changeVals[3];
1259 this->yDimSize = changeVals[4];
1260 this->zDimSize = changeVals[5];
1263 proshade_single xMov = -( ( this->xFrom * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices - changeVals[0] ) ) ) -
1264 ( this->xFrom * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices ) ) ) );
1265 proshade_single yMov = -( ( this->yFrom * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices - changeVals[1] ) ) ) -
1266 ( this->yFrom * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices ) ) ) );
1267 proshade_single zMov = -( ( this->zFrom * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices - changeVals[2] ) ) ) -
1268 ( this->zFrom * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices ) ) ) );
1272 &this->yFrom, &this->yTo, &this->zFrom, &this->zTo, &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
1275 this->xDimIndices, this->yDimIndices, this->zDimIndices );
1278 delete[] changeVals;
1300 proshade_unsign arrPos = 0;
1301 proshade_single xCOM = 0.0;
1302 proshade_single yCOM = 0.0;
1303 proshade_single zCOM = 0.0;
1304 proshade_single totDens = 0.0;
1307 for ( proshade_unsign xIt = 0; xIt < this->xDimIndices; xIt++ )
1309 for ( proshade_unsign yIt = 0; yIt < this->yDimIndices; yIt++ )
1311 for ( proshade_unsign zIt = 0; zIt < this->zDimIndices; zIt++ )
1314 arrPos = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
1317 if ( this->internalMap[arrPos] > 0.0 )
1319 xCOM +=
static_cast<proshade_single
> ( this->internalMap[arrPos] * xIt );
1320 yCOM +=
static_cast<proshade_single
> ( this->internalMap[arrPos] * yIt );
1321 zCOM +=
static_cast<proshade_single
> ( this->internalMap[arrPos] * zIt );
1322 totDens +=
static_cast<proshade_single
> ( this->internalMap[arrPos] );
1332 proshade_single xDist = (
static_cast<proshade_single
> ( this->xDimIndices / 2.0 ) - xCOM ) *
static_cast<proshade_single
> ( this->xDimSize / this->xDimIndices );
1333 proshade_single yDist = (
static_cast<proshade_single
> ( this->yDimIndices / 2.0 ) - yCOM ) *
static_cast<proshade_single
> ( this->yDimSize / this->yDimIndices );
1334 proshade_single zDist = (
static_cast<proshade_single
> ( this->zDimIndices / 2.0 ) - zCOM ) *
static_cast<proshade_single
> ( this->zDimSize / this->zDimIndices );
1358 std::stringstream hlpSS;
1359 hlpSS <<
"Adding extra " << settings->
addExtraSpace <<
" angstroms.";
1368 this->xDimSize +=
static_cast<proshade_single
> ( 2 * xAddIndices ) *
static_cast<proshade_single
> ( this->xDimSize / this->xDimIndices );
1369 this->yDimSize +=
static_cast<proshade_single
> ( 2 * yAddIndices ) *
static_cast<proshade_single
> ( this->yDimSize / this->yDimIndices );
1370 this->zDimSize +=
static_cast<proshade_single
> ( 2 * zAddIndices ) *
static_cast<proshade_single
> ( this->zDimSize / this->zDimIndices );
1372 this->xDimIndices += 2 * xAddIndices;
1373 this->yDimIndices += 2 * yAddIndices;
1374 this->zDimIndices += 2 * zAddIndices;
1376 this->xGridIndices = this->xDimIndices;
1377 this->yGridIndices = this->yDimIndices;
1378 this->zGridIndices = this->zDimIndices;
1380 this->xAxisOrigin -= xAddIndices;
1381 this->yAxisOrigin -= yAddIndices;
1382 this->zAxisOrigin -= zAddIndices;
1384 this->xFrom -= xAddIndices;
1385 this->yFrom -= yAddIndices;
1386 this->zFrom -= zAddIndices;
1388 this->xTo += xAddIndices;
1389 this->yTo += yAddIndices;
1390 this->zTo += zAddIndices;
1393 proshade_double* newMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1397 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1403 proshade_unsign newMapIndex, oldMapIndex;
1404 for ( proshade_unsign xIt = 0; xIt < (this->xDimIndices - xAddIndices); xIt++ )
1407 if ( xIt < xAddIndices ) {
continue; }
1409 for ( proshade_unsign yIt = 0; yIt < (this->yDimIndices - yAddIndices); yIt++ )
1412 if ( yIt < yAddIndices ) {
continue; }
1414 for ( proshade_unsign zIt = 0; zIt < (this->zDimIndices - zAddIndices); zIt++ )
1417 if ( zIt < zAddIndices ) {
continue; }
1420 newMapIndex = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
1421 oldMapIndex = (zIt - zAddIndices) + (this->zDimIndices - ( 2 * zAddIndices ) ) * ( (yIt - yAddIndices) + (this->yDimIndices - ( 2 * yAddIndices ) ) * (xIt - xAddIndices) );
1423 newMap[newMapIndex] = this->internalMap[oldMapIndex];
1429 delete[] this->internalMap;
1431 this->internalMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1434 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1436 this->internalMap[iter] = newMap[iter];
1467 if ( settings->
invertMap ) { this->invertMirrorMap ( settings ); }
1471 if ( settings->
normaliseMap ) { this->normaliseMap ( settings ); }
1475 if ( settings->
maskMap ) { this->maskMap ( settings ); }
1479 if ( settings->
moveToCOM ) { this->centreMapOnCOM ( settings ); }
1483 if ( settings->
addExtraSpace != 0.0 ) { this->addExtraSpace ( settings ); }
1510 if ( this->spherePos.size() != 0 )
1512 std::stringstream hlpSS;
1513 hlpSS <<
"The sphere distances were determined as " << this->spherePos.at(0);
1514 for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ ) { hlpSS <<
"; " << this->spherePos.at(iter); }
1515 hlpSS <<
" Angstroms.";
1521 proshade_unsign maxDim = std::max ( this->xDimSize, std::max ( this->yDimSize, this->zDimSize ) );
1522 proshade_unsign minDim = std::min ( this->xDimSize, std::min ( this->yDimSize, this->zDimSize ) );
1523 proshade_unsign midDim = 0;
1524 if ( ( this->xDimSize < maxDim ) && ( this->xDimSize > minDim ) ) { midDim = this->xDimSize; }
1525 else if ( ( this->yDimSize < maxDim ) && ( this->yDimSize > minDim ) ) { midDim = this->yDimSize; }
1526 else { midDim = this->zDimSize; }
1528 proshade_single maxDiag = std::sqrt ( std::pow (
static_cast<proshade_single
> ( maxDim ), 2.0 ) +
1529 std::pow (
static_cast<proshade_single
> ( midDim ), 2.0 ) );
1532 for ( proshade_single iter = 0.5; ( iter * settings->
maxSphereDists ) < ( maxDiag / 2.0 ); iter += 1.0 )
1538 this->noSpheres =
static_cast<proshade_unsign
> ( this->spherePos.size() );
1541 std::stringstream hlpSS;
1542 hlpSS <<
"The sphere distances were determined as " << this->spherePos.at(0);
1543 for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ ) { hlpSS <<
"; " << this->spherePos.at(iter); }
1544 hlpSS <<
" Angstroms.";
1571 this->xDimSize, this->yDimSize, this->zDimSize );
1575 this->getSpherePositions ( settings );
1580 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ )
1582 std::stringstream ss;
1583 ss <<
"Now mapping sphere " << iter <<
" .";
1587 this->xDimSize, this->yDimSize, this->zDimSize, iter,
1589 this->internalMap, &this->maxShellBand );
1614 this->sphericalHarmonics =
new proshade_complex* [this->noSpheres];
1616 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
1618 this->sphericalHarmonics[iter] =
new proshade_complex [(this->spheres[iter]->getLocalBandwidth() * 2) * (this->spheres[iter]->getLocalBandwidth() * 2)];
1623 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
1626 std::stringstream ss;
1627 ss <<
"Now decomposing sphere " << iter <<
". " <<
"( Band is: " << this->spheres[iter]->getLocalBandwidth() <<
").";
1656 std::vector< proshade_double* > CSyms = this->getCyclicSymmetriesList ( settings );
1662 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
1663 std::vector< proshade_double* > ISyms = this->getIcosahedralSymmetriesList ( settings, &CSyms );
1664 std::vector< proshade_double* > OSyms; std::vector< proshade_double* > TSyms;
1665 if ( ISyms.size() < 31 ) { OSyms = this->getOctahedralSymmetriesList ( settings, &CSyms );
if ( OSyms.size() < 13 ) { TSyms = this->getTetrahedralSymmetriesList ( settings, &CSyms ); } }
1668 this->saveRecommendedSymmetry ( settings, &CSyms, &DSyms, &TSyms, &OSyms, &ISyms, axes );
1674 this->saveRequestedSymmetryC ( settings, &CSyms, axes );
1680 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
1681 this->saveRequestedSymmetryD ( settings, &DSyms, axes );
1687 std::vector< proshade_double* > TSyms = this->getTetrahedralSymmetriesList ( settings, &CSyms );
1696 std::vector< proshade_double* > OSyms = this->getOctahedralSymmetriesList ( settings, &CSyms );
1705 std::vector< proshade_double* > ISyms = this->getIcosahedralSymmetriesList ( settings, &CSyms );
1713 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." );
1717 bool isArgSameAsSettings =
true;
1718 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSyms.size() ); cIt++ )
1720 std::vector< proshade_double > nextSym;
1729 if ( ( cIt == 0 ) && ( settings->
allDetectedCAxes.size() == 0 ) ) { isArgSameAsSettings =
false; }
1736 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( CSyms.size() ); it++ ) {
delete[] CSyms.at(it); }
1776 if ( settings->axisErrToleranceDefault )
1778 settings->
axisErrTolerance = std::max ( 0.01, ( 2.0 * M_PI ) / this->maxShellBand );
1785 std::stringstream hlpSS;
1790 proshade_double symThres = 0.0;
1791 std::vector< proshade_double* > CSyms = this->findRequestedCSymmetryFromAngleAxis ( settings, settings->
requestedSymmetryFold, &symThres );
1795 if ( CSyms.size() > 0 )
1801 this->saveDetectedSymmetries ( settings, &CSyms, allCs );
1817 std::vector< proshade_double* > CSyms = getCyclicSymmetriesListFromAngleAxis ( settings );
1820 if ( this->sphereMappedRotFun.size() < 1 )
1822 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." );
1828 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." );
1835 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
1836 std::vector< proshade_double* > ISyms = this->getPredictedIcosahedralSymmetriesList ( settings, &CSyms );
1837 std::vector< proshade_double* > OSyms = this->getPredictedOctahedralSymmetriesList ( settings, &CSyms );
1838 std::vector< proshade_double* > TSyms = this->getPredictedTetrahedralSymmetriesList ( settings, &CSyms );;
1841 this->saveRecommendedSymmetry ( settings, &CSyms, &DSyms, &TSyms, &OSyms, &ISyms, axes );
1847 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
1848 this->saveRequestedSymmetryD ( settings, &DSyms, axes );
1854 std::vector< proshade_double* > TSyms = this->getPredictedTetrahedralSymmetriesList ( settings, &CSyms );
1862 std::vector< proshade_double* > OSyms = this->getPredictedOctahedralSymmetriesList ( settings, &CSyms );
1870 std::vector< proshade_double* > ISyms = this->getPredictedIcosahedralSymmetriesList ( settings, &CSyms );
1876 this->saveDetectedSymmetries ( settings, &CSyms, allCs );
1899 bool isArgSameAsSettings =
true;
1902 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSyms->size() ); cIt++ )
1905 std::vector< proshade_double > nextSym;
1915 if ( ( cIt == 0 ) && ( settings->
allDetectedCAxes.size() == 0 ) ) { isArgSameAsSettings =
false; }
1920 delete[] CSyms->at(cIt);
1941 if ( CSym->size() == 0 ) { *symInd = 0;
return ( 0.0 ); }
1944 proshade_double ret = CSym->at(0)[5];
1946 proshade_double frac = 0.0;
1950 for ( proshade_unsign ind = 1; ind < static_cast<proshade_unsign>( CSym->size() ); ind++ )
1953 if ( CSym->at(ind)[0] > CSym->at(*symInd)[0] )
1956 frac = ( std::abs( CSym->at(ind)[5]- 0.5 ) / std::abs( CSym->at(*symInd)[5] - 0.5 ) ) / ( CSym->at(*symInd)[0] / CSym->at(ind)[0] );
1959 if ( frac >= 1.0 && ( ( CSym->at(*symInd)[5] * 0.85 ) < CSym->at(ind)[5] ) )
1963 ret = CSym->at(ind)[5];
1988 proshade_double ret = 0.0;
1989 proshade_double frac = 0.0;
1990 if ( DSym->size() > 0 )
1992 ret = ( ( DSym->at(0)[0] * DSym->at(0)[5] ) + ( DSym->at(0)[6] * DSym->at(0)[11] ) ) / ( DSym->at(0)[0] + DSym->at(0)[6] );
1995 else {
return ( ret ); }
1999 for ( proshade_unsign ind = 1; ind < static_cast<proshade_unsign>( DSym->size() ); ind++ )
2002 if ( ( DSym->at(ind)[0] + DSym->at(ind)[6] ) > ( DSym->at(*symInd)[0] + DSym->at(*symInd)[6] ) )
2005 frac = std::max ( std::min ( ( ( DSym->at(*symInd)[0] + DSym->at(*symInd)[6] ) / ( DSym->at(ind)[0] + DSym->at(ind)[6] ) ) * 1.5, 0.9 ), 0.6 );
2008 if ( ( ( ( DSym->at(*symInd)[0] * DSym->at(*symInd)[5] ) + ( DSym->at(*symInd)[6] * DSym->at(*symInd)[11] ) ) / ( DSym->at(*symInd)[0] + DSym->at(*symInd)[6] ) * frac ) < ( ( DSym->at(ind)[0] * DSym->at(ind)[5] ) + ( DSym->at(ind)[6] * DSym->at(ind)[11] ) ) / ( DSym->at(ind)[0] + DSym->at(ind)[6] ) )
2012 ret = ( ( DSym->at(ind)[0] * DSym->at(ind)[5] ) + ( DSym->at(ind)[6] * DSym->at(ind)[11] ) ) / ( DSym->at(ind)[0] + DSym->at(ind)[6] );
2033 proshade_double ret = 0.0;
2034 proshade_double foldSum = 0.0;
2037 if ( TSym->size() == 7 )
2040 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( TSym->size() ); cIt++ )
2042 ret += TSym->at(cIt)[0] * TSym->at(cIt)[5];
2043 foldSum += TSym->at(cIt)[0];
2066 proshade_double ret = 0.0;
2067 proshade_double foldSum = 0.0;
2070 if ( OSym->size() == 13 )
2073 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( OSym->size() ); cIt++ )
2075 ret += OSym->at(cIt)[0] * OSym->at(cIt)[5];
2076 foldSum += OSym->at(cIt)[0];
2099 proshade_double ret = 0.0;
2100 proshade_double foldSum = 0.0;
2103 if ( ISym->size() == 31 )
2106 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( ISym->size() ); cIt++ )
2108 ret += ISym->at(cIt)[0] * ISym->at(cIt)[5];
2109 foldSum += ISym->at(cIt)[0];
2138 proshade_double cScore = 0.0, dScore = 0.0, tScore = 0.0, oScore = 0.0, iScore = 0.0;
2139 proshade_unsign bestCIndex, bestDIndex;
2142 cScore = this->findBestCScore ( CSym, &bestCIndex );
2143 dScore = this->findBestDScore ( DSym, &bestDIndex );
2144 tScore = this->findTScore ( TSym );
2145 oScore = this->findOScore ( OSym );
2146 iScore = this->findIScore ( ISym );
2149 proshade_double bestWeightedScore = std::max ( cScore, std::max ( dScore * 1.1, std::max ( tScore * 3000.0, std::max ( oScore * 4000.0, iScore * 5000.0 ) ) ) );
2154 if ( bestWeightedScore == cScore )
2162 if ( ( ( 360.0 /
static_cast<double> ( CSym->at(bestCIndex)[0] ) ) - ( 360.0 /
static_cast<double> ( CSym->at(bestCIndex)[0] + 1 ) ) ) <
2163 ( 360.0 /
static_cast<double> ( settings->
maxBandwidth * 4.0 ) ) )
2165 std::stringstream hlpSS;
2166 hlpSS <<
"!!! ProSHADE WARNING !!! Reporting symmetry C" << CSym->at(bestCIndex)[0] <<
", however, the grid sampling does not provide reasonable accuracy for symmetry with such high fold and therefore ProSHADE cannot responsibly claim this symmetry to be correct. It is suggested that the grid sampling is increased for more accurate symmetry detection. (Set higher resolution using -r).";
2170 if ( bestWeightedScore == dScore * 1.1 )
2173 settings->
setRecommendedFold ( std::max ( DSym->at(bestDIndex)[0], DSym->at(bestDIndex)[6] ) );
2183 if ( ( ( 360.0 /
static_cast<double> ( std::max ( DSym->at(bestDIndex)[0], DSym->at(bestDIndex)[6] ) ) ) - ( 360.0 /
static_cast<double> ( std::max ( DSym->at(bestDIndex)[0], DSym->at(bestDIndex)[6] ) + 1 ) ) ) <
2184 ( 360.0 /
static_cast<double> ( settings->
maxBandwidth * 4.0 ) ) )
2186 std::stringstream hlpSS;
2187 hlpSS <<
"!!! ProSHADE WARNING !!! Reporting symmetry D" << std::max ( DSym->at(bestDIndex)[0], DSym->at(bestDIndex)[6] ) <<
", however, the grid sampling does not provide reasonable accuracy for symmetry with such high fold and therefore ProSHADE cannot responsibly claim this symmetry to be correct. It is suggested that the grid sampling is increased for more accurate symmetry detection. (Set higher resolution using -r).";
2191 if ( bestWeightedScore == tScore * 3000.0 )
2196 if ( settings->
detectedSymmetry.size() == 0 ) {
for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( TSym->size() ); it++ ) { settings->
setDetectedSymmetry ( TSym->at(it) ); } }
2198 if ( bestWeightedScore == oScore * 4000.0 )
2203 if ( settings->
detectedSymmetry.size() == 0 ) {
for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( OSym->size() ); it++ ) { settings->
setDetectedSymmetry ( OSym->at(it) ); } }
2205 if ( bestWeightedScore == iScore * 5000.0 )
2210 if ( settings->
detectedSymmetry.size() == 0 ) {
for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( ISym->size() ); it++ ) { settings->
setDetectedSymmetry ( ISym->at(it) ); } }
2232 proshade_unsign bestIndex = 0;
2233 proshade_double highestSym = 0.0;
2236 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( CSym->size() ); iter++ )
2242 if ( CSym->at(iter)[5] > highestSym )
2244 highestSym = CSym->at(iter)[5];
2250 if ( highestSym > 0.0 )
2283 proshade_unsign bestIndex = 0;
2284 proshade_double highestSym = 0.0;
2287 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( DSym->size() ); iter++ )
2290 if ( std::max ( DSym->at(iter)[0], DSym->at(iter)[6] ) != settings->
requestedSymmetryFold ) {
continue; }
2293 if ( ( DSym->at(iter)[5] + DSym->at(iter)[11] ) > highestSym )
2295 highestSym = ( DSym->at(iter)[5] + DSym->at(iter)[11] );
2301 if ( highestSym > 0.0 )
2304 settings->
setRecommendedFold ( std::max ( DSym->at(bestIndex)[0], DSym->at(bestIndex)[6] ) );
2336 std::vector< proshade_double > angList;
2337 std::vector<std::vector< proshade_double > > ret;
2340 proshade_double* rotMat =
new proshade_double[9];
2345 proshade_double normF = std::sqrt( std::pow ( xAx, 2.0 ) + std::pow ( yAx, 2.0 ) + std::pow ( zAx, 2.0 ) );
2351 if ( fold % 2 == 0 )
2354 for ( proshade_double iter =
static_cast < proshade_double
> ( -( ( fold / 2 ) - 1 ) ); iter <= static_cast < proshade_double > ( fold / 2 ); iter++ )
2362 for ( proshade_double iter =
static_cast < proshade_double
> ( -fold / 2 ); iter <= static_cast < proshade_double > ( fold / 2 ); iter++ )
2369 for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( angList.size() ); iter++ )
2375 std::vector < proshade_double > retEl;
2376 for ( proshade_unsign matIt = 0; matIt < 9; matIt++ )
2400 if ( obtainedAxes != requiredAxes )
2402 std::stringstream hlpSS;
2403 hlpSS <<
"The supplied number of axes for group element\n : detection ( >" << obtainedAxes <<
"< ) does not match the group type ( >" << groupType <<
"< ).";
2404 throw ProSHADE_exception (
"Mismatch between supplied number of axes and\n : symmetry type.",
"ES00059", __FILE__, __LINE__, __func__, hlpSS.str() );
2419 bool checkElementAlreadyExists ( std::vector<std::vector< proshade_double > >* elements, std::vector< proshade_double >* elem, proshade_double matrixTolerance )
2422 bool elementFound =
false;
2425 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( elements->size() ); elIt++ )
2429 elementFound =
true;
2435 return ( elementFound );
2448 bool isGroup =
true;
2451 for ( proshade_unsign gr1 = 0; gr1 < static_cast<proshade_unsign> ( elements->size() ); gr1++ )
2453 for ( proshade_unsign gr2 = 1; gr2 < static_cast<proshade_unsign> ( elements->size() ); gr2++ )
2456 if ( gr1 >= gr2 ) {
continue; }
2470 if ( !isGroup ) {
break; }
2489 std::vector< std::vector< proshade_double > > ret;
2492 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( first->size() ); elIt++ )
2501 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( second->size() ); elIt++ )
2512 for ( proshade_unsign gr1 = 0; gr1 < static_cast<proshade_unsign> ( first->size() ); gr1++ )
2514 for ( proshade_unsign gr2 = 0; gr2 < static_cast<proshade_unsign> ( second->size() ); gr2++ )
2555 std::vector<std::vector< proshade_double > > ret;
2558 if ( groupType ==
"C" )
2573 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." );
2576 else if ( groupType ==
"D" )
2598 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." );
2601 else if ( groupType ==
"T" )
2607 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2624 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2644 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." );
2647 else if ( groupType ==
"O" )
2653 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2670 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2687 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2707 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." );
2710 else if ( groupType ==
"I" )
2716 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2733 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2750 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2770 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." );
2773 else if ( groupType ==
"X" )
2776 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2792 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." );
2797 std::stringstream hlpSS;
2798 hlpSS <<
"Unknown symmetry type: >" << groupType <<
"<";
2799 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." );
2817 if ( saveTo != NULL )
2824 saveTo =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
2830 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
2832 saveTo[iter] = this->internalMap[iter];
2855 std::stringstream ssHlp;
2861 ssHlp.clear(); ssHlp.str (
"" );
2862 ssHlp <<
" Fold X Y Z Angle Height";
2865 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->
detectedSymmetry.size() ); symIt++ )
2867 ssHlp.clear(); ssHlp.str (
"" );
2872 std::stringstream hlpSS3;
2873 ssHlp.clear(); ssHlp.str (
"" );
2874 hlpSS3 << std::endl <<
"However, since the selection of the recommended symmetry needs improvement, here is a list of all detected C symmetries:";
2879 ssHlp.clear(); ssHlp.str (
"" );
2880 ssHlp <<
" Fold X Y Z Angle Height";
2883 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->
allDetectedCAxes.size() ); symIt++ )
2885 ssHlp.clear(); ssHlp.str (
"" );
2890 hlpSS3.clear(); hlpSS3.str (
"" );
2891 hlpSS3 << std::endl <<
"Also, for the same reason, here is a list of all detected D symmetries:";
2896 ssHlp.clear(); ssHlp.str (
"" );
2897 ssHlp <<
" Fold X Y Z Angle Height";
2900 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->
allDetectedDAxes.size() ); symIt++ )
2902 ssHlp.clear(); ssHlp.str (
"" );
2906 for ( proshade_unsign axIt = 1; axIt < static_cast<proshade_unsign> ( settings->
allDetectedDAxes.at(symIt).size() ); axIt++ )
2908 ssHlp.clear(); ssHlp.str (
"" );
2913 ssHlp.clear(); ssHlp.str (
"" );
2934 proshade_double totNonZeroPoints = 0.0;
2935 proshade_signed mapIt = 0;
2938 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
2940 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
2942 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
2945 mapIt = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
2948 if ( this->internalMap[mapIt] <= 0.0 ) {
continue; }
2951 this->xCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( xIt + this->xFrom );
2952 this->yCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( yIt + this->yFrom );
2953 this->zCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( zIt + this->zFrom );
2954 totNonZeroPoints += this->internalMap[mapIt];
2959 this->xCom /= totNonZeroPoints;
2960 this->yCom /= totNonZeroPoints;
2961 this->zCom /= totNonZeroPoints;
2964 this->xCom = (
static_cast<proshade_double
> ( this->xFrom ) * ( this->xDimSizeOriginal /
static_cast<proshade_double
> ( this->xDimIndicesOriginal ) ) ) +
2965 ( ( this->xCom -
static_cast<proshade_double
> ( this->xFrom ) ) *
2966 ( this->xDimSizeOriginal /
static_cast<proshade_double
> ( this->xDimIndicesOriginal ) ) );
2967 this->yCom = (
static_cast<proshade_double
> ( this->yFrom ) * ( this->yDimSizeOriginal /
static_cast<proshade_double
> ( this->yDimIndicesOriginal ) ) ) +
2968 ( ( this->yCom -
static_cast<proshade_double
> ( this->yFrom ) ) *
2969 ( this->yDimSizeOriginal /
static_cast<proshade_double
> ( this->yDimIndicesOriginal ) ) );
2970 this->zCom = (
static_cast<proshade_double
> ( this->zFrom ) * ( this->zDimSizeOriginal /
static_cast<proshade_double
> ( this->zDimIndicesOriginal ) ) ) +
2971 ( ( this->zCom -
static_cast<proshade_double
> ( this->zFrom ) ) *
2972 ( this->zDimSizeOriginal /
static_cast<proshade_double
> ( this->zDimIndicesOriginal ) ) );
2986 return ( this->noSpheres );
2997 return ( this->internalMap[pos] );
3007 return ( this->maxShellBand );
3017 return ( this->rrpMatrices[band][sh1][sh2] );
3032 if ( this->spheres[shell]->getLocalBandwidth( ) >= bandVal )
3056 fftw_complex* mapCoeffs =
new fftw_complex[this->xDimIndices * this->yDimIndices * this->zDimIndices];
3057 fftw_complex* pattersonMap =
new fftw_complex[this->xDimIndices * this->yDimIndices * this->zDimIndices];
3064 for ( proshade_unsign iter = 0; iter < (this->xDimIndices * this->yDimIndices * this->zDimIndices); iter++ )
3066 pattersonMap[iter][0] = this->internalMap[iter];
3067 pattersonMap[iter][1] = 0.0;
3071 fftw_plan forward = fftw_plan_dft_3d ( this->xDimIndices, this->yDimIndices, this->zDimIndices,
3072 pattersonMap, mapCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
3073 fftw_plan inverse = fftw_plan_dft_3d ( this->xDimIndices, this->yDimIndices, this->zDimIndices,
3074 mapCoeffs, pattersonMap, FFTW_BACKWARD, FFTW_ESTIMATE );
3077 fftw_execute ( forward );
3083 fftw_execute ( inverse );
3086 proshade_signed mapIt, patIt, patX, patY, patZ;
3087 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
3089 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
3091 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
3094 patX = xIt - (
static_cast<proshade_signed
> ( this->xDimIndices ) / 2 );
if ( patX < 0 ) { patX += this->xDimIndices; }
3095 patY = yIt - (
static_cast<proshade_signed
> ( this->yDimIndices ) / 2 );
if ( patY < 0 ) { patY += this->yDimIndices; }
3096 patZ = zIt - (
static_cast<proshade_signed
> ( this->zDimIndices ) / 2 );
if ( patZ < 0 ) { patZ += this->zDimIndices; }
3099 mapIt = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
3100 patIt = patZ + this->zDimIndices * ( patY + this->yDimIndices * patX );
3103 this->internalMap[mapIt] = pattersonMap[patIt][0];
3109 delete[] pattersonMap;
3113 fftw_destroy_plan ( forward );
3114 fftw_destroy_plan ( inverse );
3131 return ( &this->sphericalHarmonics[shell][seanindex (
static_cast<proshade_signed
> ( order ) -
static_cast<proshade_signed
> ( band ),
3133 this->spheres[shell]->getLocalBandwidth() )][0] );
3144 return ( &this->sphericalHarmonics[shell][seanindex (
static_cast<proshade_signed
> ( order ) -
static_cast<proshade_signed
> ( band ),
3146 this->spheres[shell]->getLocalBandwidth() )][1] );
3157 return ( this->spheres[shell]->getShellRadius() );
3168 return ( this->integrationWeight );
3180 return ( this->spheres[shell]->getLocalBandwidth ( ) );
3192 return ( this->spherePos.at(shell) );
3204 return ( this->eMatrices[band] );
3219 *valueReal = this->eMatrices[band][order1][order2][0];
3220 *valueImag = this->eMatrices[band][order1][order2][1];
3234 return ( this->so3CoeffsInverse );
3245 return ( this->so3Coeffs );
3256 return ( this->maxCompBand );
3271 *valueReal = this->wignerMatrices[band][order1][order2][0];
3272 *valueImag = this->wignerMatrices[band][order1][order2][1];
3286 return ( this->xDimSize );
3296 return ( this->yDimSize );
3306 return ( this->zDimSize );
3316 return ( this->xDimIndices );
3326 return ( this->yDimIndices );
3336 return ( this->zDimIndices );
3346 return ( &this->xFrom );
3356 return ( &this->yFrom );
3366 return ( &this->zFrom );
3376 return ( &this->xTo );
3386 return ( &this->yTo );
3396 return ( &this->zTo );
3406 return ( &this->xAxisOrigin );
3416 return ( &this->yAxisOrigin );
3426 return ( &this->zAxisOrigin );
3436 return ( this->internalMap );
3446 return ( this->translationMap );
3456 this->integrationWeight = intW;
3470 this->integrationWeight += intW;
3487 this->eMatrices[band][order1][order2][0] = val[0];
3488 this->eMatrices[band][order1][order2][1] = val[1];
3505 this->eMatrices[band][order1][order2][0] /= normF;
3506 this->eMatrices[band][order1][order2][1] /= normF;
3521 this->so3Coeffs[position][0] = val[0];
3522 this->so3Coeffs[position][1] = val[1];
3539 this->wignerMatrices[band][order1][order2][0] = val[0];
3540 this->wignerMatrices[band][order1][order2][1] = val[1];
3557 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3559 eMatsLMReal[iter] =
static_cast<double> ( this->eMatrices[band][order1][iter][0] );
3577 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3579 eMatsLMImag[iter] =
static_cast<double> ( this->eMatrices[band][order1][iter][1] );
3595 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3597 so3CoefsReal[iter] =
static_cast<double> ( this->so3Coeffs[iter][0] );
3613 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3615 so3CoefsImag[iter] =
static_cast<double> ( this->so3Coeffs[iter][1] );
3635 return (
static_cast<int> ( so3CoefLoc ( order1, order2, band, this->getMaxBand() ) ) );
3646 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3648 rotFunReal[iter] =
static_cast<double> ( this->so3CoeffsInverse[iter][0] );
3664 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3666 rotFunImag[iter] =
static_cast<double> ( this->so3CoeffsInverse[iter][1] );
3682 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3684 trsFunReal[iter] =
static_cast<double> ( this->translationMap[iter][0] );
3700 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3702 trsFunImag[iter] =
static_cast<double> ( this->translationMap[iter][1] );
3721 proshade_double eA, eB, eG;
3725 proshade_double* rMat = NULL;
3726 rMat =
new proshade_double[9];
3733 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3735 rotMat[iter] =
static_cast<double> ( rMat[iter] );
3776 return (
static_cast<proshade_unsign
> ( settings->
detectedSymmetry.size() ) );
3788 if (
static_cast<proshade_unsign
> ( settings->
detectedSymmetry.size() ) <= axisNo )
3791 return ( std::vector< std::string > ( ) );
3795 std::vector< std::string > ret;
3798 std::stringstream ssHlp;
3843 std::stringstream fNameHlp;
3845 this->writeMap ( fNameHlp.str() );
3852 this->writePdb ( fNameHlp.str(), eulA, eulB, eulG, ultimateTranslation->at(0), ultimateTranslation->at(1), ultimateTranslation->at(2), settings->
firstModelOnly );
3858 ultimateTranslation->at(0), ultimateTranslation->at(1), ultimateTranslation->at(2),
3880 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);
3884 proshade_double* rotMat =
new proshade_double[9];
3888 std::stringstream rotMatSS;
3889 rotMatSS << std::setprecision (3) << std::showpos <<
"The rotation matrix about origin is : " << rotMat[0] <<
" " << rotMat[1] <<
" " << rotMat[2] << std::endl;
3890 rotMatSS << std::setprecision (3) << std::showpos <<
" : " << rotMat[3] <<
" " << rotMat[4] <<
" " << rotMat[5] << std::endl;
3891 rotMatSS << std::setprecision (3) << std::showpos <<
" : " << rotMat[6] <<
" " << rotMat[7] <<
" " << rotMat[8];
3897 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);