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 ( );
537 this->reSampleMap ( settings );
541 this->xDimSizeOriginal = this->xDimSize;
542 this->yDimSizeOriginal = this->yDimSize;
543 this->zDimSizeOriginal = this->zDimSize;
546 this->xDimIndicesOriginal = this->xDimIndices;
547 this->yDimIndicesOriginal = this->yDimIndices;
548 this->zDimIndicesOriginal = this->zDimIndices;
551 this->xAxisOriginOriginal = this->xAxisOrigin;
552 this->yAxisOriginOriginal = this->yAxisOrigin;
553 this->zAxisOriginOriginal = this->zAxisOrigin;
556 this->findMapCOM ( );
557 this->originalMapXCom = this->xCom;
558 this->originalMapYCom = this->yCom;
559 this->originalMapZCom = this->zCom;
583 gemmi::Structure pdbFile = gemmi::read_structure ( gemmi::MaybeGzipped ( this->fileName ) );
598 proshade_double xCOMPdb, yCOMPdb, zCOMPdb;
602 proshade_single xF, xT, yF, yT, zF, zT;
606 proshade_single xMov = 20.0 - xF;
607 proshade_single yMov = 20.0 - yF;
608 proshade_single zMov = 20.0 - zF;
612 this->xDimSize = xT - xF + 40.0;
613 this->yDimSize = yT - yF + 40.0;
614 this->zDimSize = zT - zF + 40.0;
617 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 );
620 this->setPDBMapValues ( );
623 proshade_double xCOMMap, yCOMMap, zCOMMap;
625 this->xDimSize, this->yDimSize, this->zDimSize,
626 this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
628 if ( pdbFile.models.size() > 1 )
630 xMov = xCOMMap - xCOMPdb;
631 yMov = yCOMMap - yCOMPdb;
632 zMov = zCOMMap - zCOMPdb;
637 &this->xFrom, &this->xTo, &this->yFrom, &this->yTo, &this->zFrom, &this->zTo,
638 &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
640 this->xDimIndices, this->yDimIndices, this->zDimIndices );
645 this->reSampleMap ( settings );
649 this->xDimSizeOriginal = this->xDimSize;
650 this->yDimSizeOriginal = this->yDimSize;
651 this->zDimSizeOriginal = this->zDimSize;
654 this->xDimIndicesOriginal = this->xDimIndices;
655 this->yDimIndicesOriginal = this->yDimIndices;
656 this->zDimIndicesOriginal = this->zDimIndices;
659 this->xAxisOriginOriginal = this->xAxisOrigin;
660 this->yAxisOriginOriginal = this->yAxisOrigin;
661 this->zAxisOriginOriginal = this->zAxisOrigin;
664 this->findMapCOM ( );
665 this->originalMapXCom = this->xCom;
666 this->originalMapYCom = this->yCom;
667 this->originalMapZCom = this->zCom;
691 this->xDimIndices = this->xTo;
692 this->yDimIndices = this->yTo;
693 this->zDimIndices = this->zTo;
696 this->xGridIndices = this->xDimIndices;
697 this->yGridIndices = this->yDimIndices;
698 this->zGridIndices = this->zDimIndices;
701 this->xAxisOrder = 1;
702 this->yAxisOrder = 2;
703 this->zAxisOrder = 3;
706 this->xAxisOrigin = this->xFrom;
707 this->yAxisOrigin = this->yFrom;
708 this->zAxisOrigin = this->zFrom;
722 this->xTo = this->xFrom + this->xDimIndices - 1;
723 this->yTo = this->yFrom + this->yDimIndices - 1;
724 this->zTo = this->zFrom + this->zDimIndices - 1;
744 gemmi::Grid<float> mapData;
745 mapData.set_unit_cell ( this->xDimSize, this->yDimSize, this->zDimSize, this->aAngle, this->bAngle, this->cAngle );
746 mapData.set_size_without_checking ( this->xDimIndices, this->yDimIndices, this->zDimIndices );
747 mapData.axis_order = gemmi::AxisOrder::XYZ;
748 mapData.spacegroup = &gemmi::get_spacegroup_p1();
751 gemmi::Ccp4<float> map;
753 map.update_ccp4_header ( mode );
757 this->xDimIndices, this->yDimIndices, this->zDimIndices,
758 this->xDimSize, this->yDimSize, this->zDimSize,
759 this->aAngle, this->bAngle, this->cAngle,
760 this->xFrom, this->yFrom, this->zFrom,
761 this->xAxisOrigin, this->yAxisOrigin, this->zAxisOrigin,
762 this->xAxisOrder, this->yAxisOrder, this->zAxisOrder,
763 this->xGridIndices, this->yGridIndices, this->zGridIndices,
767 proshade_unsign arrPos = 0;
768 for ( proshade_unsign uIt = 0; uIt < this->xDimIndices; uIt++ )
770 for ( proshade_unsign vIt = 0; vIt < this->yDimIndices; vIt++ )
772 for ( proshade_unsign wIt = 0; wIt < this->zDimIndices; wIt++ )
774 arrPos = wIt + this->zDimIndices * ( vIt + this->yDimIndices * uIt );
775 map.grid.set_value ( uIt, vIt, wIt,
static_cast<float> ( this->internalMap[arrPos] ) );
781 map.update_ccp4_header ( mode,
true );
784 map.write_ccp4_map ( fName );
811 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." );
815 gemmi::Structure pdbFile = gemmi::read_structure ( gemmi::MaybeGzipped ( this->fileName ) );
818 if ( ( euA != 0.0 ) || ( euB != 0.0 ) || ( euG != 0.0 ) )
828 std::ofstream outCoOrdFile;
829 outCoOrdFile.open ( fName.c_str() );
831 if ( outCoOrdFile.is_open() )
833 gemmi::PdbWriteOptions opt;
834 write_pdb ( pdbFile, outCoOrdFile, opt );
838 std::stringstream hlpMessage;
839 hlpMessage <<
"Failed to open the PDB file " << fName <<
" for output.";
840 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." );
843 outCoOrdFile.close ( );
860 proshade_double* hlpMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
864 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
866 hlpMap[iter] = this->internalMap[iter];
867 this->internalMap[iter] = mask[iter];
871 this->writeMap ( fName );
874 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
876 this->internalMap[iter] = hlpMap[iter];
900 proshade_signed arrayPos, invPos;
903 proshade_double* hlpMap =
new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
907 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
909 hlpMap[iter] = this->internalMap[iter];
913 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
915 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
917 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
920 arrayPos = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
921 invPos = ( (this->zDimIndices-1) - zIt ) + this->zDimIndices * ( ( (this->yDimIndices-1) - yIt ) + this->yDimIndices * ( (this->xDimIndices-1) - xIt ) );
924 this->internalMap[invPos] = hlpMap[arrayPos];
954 std::vector<proshade_double> mapVals ( this->xDimIndices * this->yDimIndices * this->zDimIndices, 0.0 );
957 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
959 mapVals.at(iter) = this->internalMap[iter];
963 proshade_double* meanSD =
new proshade_double[2];
967 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
969 this->internalMap[iter] = ( this->internalMap[iter] - meanSD[0] ) / meanSD[1];
1001 proshade_double* blurredMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1006 this->xDimSize, this->yDimSize, this->zDimSize, settings->
blurFactor );
1012 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 ); } }
1015 delete[] blurredMap;
1043 for ( proshade_unsign iter = 0; iter < 6; iter++ ) { ret[iter] = settings->
forceBounds[iter]; }
1050 this->xDimSize, this->yDimSize, this->zDimSize, ret );
1054 this->xDimSize, this->yDimSize, this->zDimSize, ret, settings->
boundsExtraSpace );
1060 std::stringstream ssHlp;
1061 ssHlp <<
"New boundaries are: " << ret[1] - ret[0] + 1 <<
" x " << ret[3] - ret[2] + 1 <<
" x " << ret[5] - ret[4] + 1;
1067 for ( proshade_unsign iter = 0; iter < 6; iter++ ) { settings->
forceBounds[iter] = ret[iter]; }
1097 newStr->
fileType = ProSHADE_internal_io::MAP;
1100 newStr->
xDimIndices =
static_cast<proshade_signed
> ( newBounds[1] ) -
static_cast<proshade_signed
> ( newBounds[0] ) + 1;
1101 newStr->
yDimIndices =
static_cast<proshade_signed
> ( newBounds[3] ) -
static_cast<proshade_signed
> ( newBounds[2] ) + 1;
1102 newStr->
zDimIndices =
static_cast<proshade_signed
> ( newBounds[5] ) -
static_cast<proshade_signed
> ( newBounds[4] ) + 1;
1104 newStr->
aAngle = this->aAngle;
1105 newStr->
bAngle = this->aAngle;
1106 newStr->
cAngle = this->aAngle;
1108 newStr->
xDimSize =
static_cast<proshade_single
> ( newStr->
xDimIndices ) * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices ) );
1109 newStr->
yDimSize =
static_cast<proshade_single
> ( newStr->
yDimIndices ) * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices ) );
1110 newStr->
zDimSize =
static_cast<proshade_single
> ( newStr->
zDimIndices ) * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices ) );
1120 newStr->
xAxisOrigin = this->xAxisOrigin + newBounds[0];
1121 newStr->
yAxisOrigin = this->yAxisOrigin + newBounds[2];
1122 newStr->
zAxisOrigin = this->zAxisOrigin + newBounds[4];
1124 newStr->
xFrom = this->xFrom + newBounds[0];
1125 newStr->
yFrom = this->yFrom + newBounds[2];
1126 newStr->
zFrom = this->zFrom + newBounds[4];
1128 newStr->
xTo = this->xTo - ( (this->xDimIndices-1) - newBounds[1] );
1129 newStr->
yTo = this->yTo - ( (this->yDimIndices-1) - newBounds[3] );
1130 newStr->
zTo = this->zTo - ( (this->zDimIndices-1) - newBounds[5] );
1139 this->xDimIndices, this->yDimIndices, this->zDimIndices, newStr->
internalMap, this->internalMap );
1159 proshade_single* changeVals =
new proshade_single[6];
1165 this->xDimSize, this->yDimSize, this->zDimSize, changeVals );
1175 this->xDimSize, this->yDimSize, this->zDimSize, changeVals );
1180 this->xDimIndices +=
static_cast<proshade_unsign
> ( changeVals[0] );
1181 this->yDimIndices +=
static_cast<proshade_unsign
> ( changeVals[1] );
1182 this->zDimIndices +=
static_cast<proshade_unsign
> ( changeVals[2] );
1184 this->xGridIndices = this->xDimIndices;
1185 this->yGridIndices = this->yDimIndices;
1186 this->zGridIndices = this->zDimIndices;
1188 this->xTo +=
static_cast<proshade_unsign
> ( changeVals[0] );
1189 this->yTo +=
static_cast<proshade_unsign
> ( changeVals[1] );
1190 this->zTo +=
static_cast<proshade_unsign
> ( changeVals[2] );
1192 this->xDimSize = changeVals[3];
1193 this->yDimSize = changeVals[4];
1194 this->zDimSize = changeVals[5];
1197 proshade_single xMov = -( ( this->xFrom * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices - changeVals[0] ) ) ) -
1198 ( this->xFrom * ( this->xDimSize /
static_cast<proshade_single
> ( this->xDimIndices ) ) ) );
1199 proshade_single yMov = -( ( this->yFrom * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices - changeVals[1] ) ) ) -
1200 ( this->yFrom * ( this->yDimSize /
static_cast<proshade_single
> ( this->yDimIndices ) ) ) );
1201 proshade_single zMov = -( ( this->zFrom * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices - changeVals[2] ) ) ) -
1202 ( this->zFrom * ( this->zDimSize /
static_cast<proshade_single
> ( this->zDimIndices ) ) ) );
1206 &this->yFrom, &this->yTo, &this->zFrom, &this->zTo, &this->xAxisOrigin, &this->yAxisOrigin, &this->zAxisOrigin );
1209 this->xDimIndices, this->yDimIndices, this->zDimIndices );
1212 delete[] changeVals;
1234 proshade_unsign arrPos = 0;
1235 proshade_single xCOM = 0.0;
1236 proshade_single yCOM = 0.0;
1237 proshade_single zCOM = 0.0;
1238 proshade_single totDens = 0.0;
1241 for ( proshade_unsign xIt = 0; xIt < this->xDimIndices; xIt++ )
1243 for ( proshade_unsign yIt = 0; yIt < this->yDimIndices; yIt++ )
1245 for ( proshade_unsign zIt = 0; zIt < this->zDimIndices; zIt++ )
1248 arrPos = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
1251 if ( this->internalMap[arrPos] > 0.0 )
1253 xCOM +=
static_cast<proshade_single
> ( this->internalMap[arrPos] * xIt );
1254 yCOM +=
static_cast<proshade_single
> ( this->internalMap[arrPos] * yIt );
1255 zCOM +=
static_cast<proshade_single
> ( this->internalMap[arrPos] * zIt );
1256 totDens +=
static_cast<proshade_single
> ( this->internalMap[arrPos] );
1266 proshade_single xDist = (
static_cast<proshade_single
> ( this->xDimIndices / 2.0 ) - xCOM ) *
static_cast<proshade_single
> ( this->xDimSize / this->xDimIndices );
1267 proshade_single yDist = (
static_cast<proshade_single
> ( this->yDimIndices / 2.0 ) - yCOM ) *
static_cast<proshade_single
> ( this->yDimSize / this->yDimIndices );
1268 proshade_single zDist = (
static_cast<proshade_single
> ( this->zDimIndices / 2.0 ) - zCOM ) *
static_cast<proshade_single
> ( this->zDimSize / this->zDimIndices );
1292 std::stringstream hlpSS;
1293 hlpSS <<
"Adding extra " << settings->
addExtraSpace <<
" angstroms.";
1302 this->xDimSize +=
static_cast<proshade_single
> ( 2 * xAddIndices ) *
static_cast<proshade_single
> ( this->xDimSize / this->xDimIndices );
1303 this->yDimSize +=
static_cast<proshade_single
> ( 2 * yAddIndices ) *
static_cast<proshade_single
> ( this->yDimSize / this->yDimIndices );
1304 this->zDimSize +=
static_cast<proshade_single
> ( 2 * zAddIndices ) *
static_cast<proshade_single
> ( this->zDimSize / this->zDimIndices );
1306 this->xDimIndices += 2 * xAddIndices;
1307 this->yDimIndices += 2 * yAddIndices;
1308 this->zDimIndices += 2 * zAddIndices;
1310 this->xGridIndices = this->xDimIndices;
1311 this->yGridIndices = this->yDimIndices;
1312 this->zGridIndices = this->zDimIndices;
1314 this->xAxisOrigin -= xAddIndices;
1315 this->yAxisOrigin -= yAddIndices;
1316 this->zAxisOrigin -= zAddIndices;
1318 this->xFrom -= xAddIndices;
1319 this->yFrom -= yAddIndices;
1320 this->zFrom -= zAddIndices;
1322 this->xTo += xAddIndices;
1323 this->yTo += yAddIndices;
1324 this->zTo += zAddIndices;
1327 proshade_double* newMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1331 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1337 proshade_unsign newMapIndex, oldMapIndex;
1338 for ( proshade_unsign xIt = 0; xIt < (this->xDimIndices - xAddIndices); xIt++ )
1341 if ( xIt < xAddIndices ) {
continue; }
1343 for ( proshade_unsign yIt = 0; yIt < (this->yDimIndices - yAddIndices); yIt++ )
1346 if ( yIt < yAddIndices ) {
continue; }
1348 for ( proshade_unsign zIt = 0; zIt < (this->zDimIndices - zAddIndices); zIt++ )
1351 if ( zIt < zAddIndices ) {
continue; }
1354 newMapIndex = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
1355 oldMapIndex = (zIt - zAddIndices) + (this->zDimIndices - ( 2 * zAddIndices ) ) * ( (yIt - yAddIndices) + (this->yDimIndices - ( 2 * yAddIndices ) ) * (xIt - xAddIndices) );
1357 newMap[newMapIndex] = this->internalMap[oldMapIndex];
1363 delete[] this->internalMap;
1365 this->internalMap =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
1368 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
1370 this->internalMap[iter] = newMap[iter];
1401 if ( settings->
invertMap ) { this->invertMirrorMap ( settings ); }
1405 if ( settings->
normaliseMap ) { this->normaliseMap ( settings ); }
1410 if ( settings->
maskMap ) { this->maskMap ( settings ); }
1414 if ( settings->
moveToCOM ) { this->centreMapOnCOM ( settings ); }
1418 if ( settings->
addExtraSpace != 0.0 ) { this->addExtraSpace ( settings ); }
1445 if ( this->spherePos.size() != 0 )
1447 std::stringstream hlpSS;
1448 hlpSS <<
"The sphere distances were determined as " << this->spherePos.at(0);
1449 for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ ) { hlpSS <<
"; " << this->spherePos.at(iter); }
1450 hlpSS <<
" Angstroms.";
1456 proshade_unsign maxDim = std::max ( this->xDimSize, std::max ( this->yDimSize, this->zDimSize ) );
1457 proshade_unsign minDim = std::min ( this->xDimSize, std::min ( this->yDimSize, this->zDimSize ) );
1458 proshade_unsign midDim = 0;
1459 if ( ( this->xDimSize < maxDim ) && ( this->xDimSize > minDim ) ) { midDim = this->xDimSize; }
1460 else if ( ( this->yDimSize < maxDim ) && ( this->yDimSize > minDim ) ) { midDim = this->yDimSize; }
1461 else { midDim = this->zDimSize; }
1463 proshade_single maxDiag = std::sqrt ( std::pow (
static_cast<proshade_single
> ( maxDim ), 2.0 ) +
1464 std::pow (
static_cast<proshade_single
> ( midDim ), 2.0 ) );
1467 for ( proshade_single iter = 0.5; ( iter * settings->
maxSphereDists ) < ( maxDiag / 2.0 ); iter += 1.0 )
1473 this->noSpheres =
static_cast<proshade_unsign
> ( this->spherePos.size() );
1476 std::stringstream hlpSS;
1477 hlpSS <<
"The sphere distances were determined as " << this->spherePos.at(0);
1478 for ( proshade_unsign iter = 1; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ ) { hlpSS <<
"; " << this->spherePos.at(iter); }
1479 hlpSS <<
" Angstroms.";
1509 this->getSpherePositions ( settings );
1514 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( this->spherePos.size() ); iter++ )
1516 std::stringstream ss;
1517 ss <<
"Now mapping sphere " << iter <<
" .";
1521 this->xDimSize, this->yDimSize, this->zDimSize, iter,
1523 this->internalMap, &this->maxShellBand );
1548 this->sphericalHarmonics =
new proshade_complex* [this->noSpheres];
1550 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
1552 this->sphericalHarmonics[iter] =
new proshade_complex [(this->spheres[iter]->getLocalBandwidth() * 2) * (this->spheres[iter]->getLocalBandwidth() * 2)];
1557 for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
1560 std::stringstream ss;
1561 ss <<
"Now decomposing sphere " << iter <<
". " <<
"( Band is: " << this->spheres[iter]->getLocalBandwidth() <<
").";
1590 std::vector< proshade_double* > CSyms = this->getCyclicSymmetriesList ( settings );
1596 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
1597 std::vector< proshade_double* > ISyms = this->getIcosahedralSymmetriesList ( settings, &CSyms );
1598 std::vector< proshade_double* > OSyms; std::vector< proshade_double* > TSyms;
1599 if ( ISyms.size() < 31 ) { OSyms = this->getOctahedralSymmetriesList ( settings, &CSyms );
if ( OSyms.size() < 13 ) { TSyms = this->getTetrahedralSymmetriesList ( settings, &CSyms ); } }
1602 this->saveRecommendedSymmetry ( settings, &CSyms, &DSyms, &TSyms, &OSyms, &ISyms, axes );
1608 this->saveRequestedSymmetryC ( settings, &CSyms, axes );
1614 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
1615 this->saveRequestedSymmetryD ( settings, &DSyms, axes );
1621 std::vector< proshade_double* > TSyms = this->getTetrahedralSymmetriesList ( settings, &CSyms );
1630 std::vector< proshade_double* > OSyms = this->getOctahedralSymmetriesList ( settings, &CSyms );
1639 std::vector< proshade_double* > ISyms = this->getIcosahedralSymmetriesList ( settings, &CSyms );
1647 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." );
1651 bool isArgSameAsSettings =
true;
1652 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSyms.size() ); cIt++ )
1654 std::vector< proshade_double > nextSym;
1663 if ( ( cIt == 0 ) && ( settings->
allDetectedCAxes.size() == 0 ) ) { isArgSameAsSettings =
false; }
1670 for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( CSyms.size() ); it++ ) {
delete[] CSyms.at(it); }
1710 if ( settings->axisErrToleranceDefault )
1712 settings->
axisErrTolerance = std::max ( 0.01, ( 2.0 * M_PI ) / this->maxShellBand );
1719 std::stringstream hlpSS;
1724 proshade_double symThres = 0.0;
1725 std::vector< proshade_double* > CSyms = this->findRequestedCSymmetryFromAngleAxis ( settings, settings->
requestedSymmetryFold, &symThres );
1729 if ( CSyms.size() > 0 )
1735 this->saveDetectedSymmetries ( settings, &CSyms, allCs );
1751 std::vector< proshade_double* > CSyms = getCyclicSymmetriesListFromAngleAxis ( settings );
1754 if ( this->sphereMappedRotFun.size() < 1 )
1756 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." );
1762 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." );
1769 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
1770 std::vector< proshade_double* > ISyms = this->getPredictedIcosahedralSymmetriesList ( settings, &CSyms );
1771 std::vector< proshade_double* > OSyms; std::vector< proshade_double* > TSyms;
1772 if ( ISyms.size() < 31 ) { OSyms = this->getPredictedOctahedralSymmetriesList ( settings, &CSyms );
1773 if ( OSyms.size() < 13 ) { TSyms = this->getPredictedTetrahedralSymmetriesList ( settings, &CSyms ); } }
1776 this->saveRecommendedSymmetry ( settings, &CSyms, &DSyms, &TSyms, &OSyms, &ISyms, axes );
1782 std::vector< proshade_double* > DSyms = this->getDihedralSymmetriesList ( settings, &CSyms );
1783 this->saveRequestedSymmetryD ( settings, &DSyms, axes );
1789 std::vector< proshade_double* > TSyms = this->getPredictedTetrahedralSymmetriesList ( settings, &CSyms );
1797 std::vector< proshade_double* > OSyms = this->getPredictedOctahedralSymmetriesList ( settings, &CSyms );
1805 std::vector< proshade_double* > ISyms = this->getPredictedIcosahedralSymmetriesList ( settings, &CSyms );
1811 this->saveDetectedSymmetries ( settings, &CSyms, allCs );
1834 bool isArgSameAsSettings =
true;
1837 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( CSyms->size() ); cIt++ )
1840 std::vector< proshade_double > nextSym;
1850 if ( ( cIt == 0 ) && ( settings->
allDetectedCAxes.size() == 0 ) ) { isArgSameAsSettings =
false; }
1855 delete[] CSyms->at(cIt);
1876 if ( CSym->size() == 0 ) { *symInd = 0;
return ( 0.0 ); }
1879 proshade_double ret = CSym->at(0)[5];
1881 proshade_double frac = 0.0;
1885 for ( proshade_unsign ind = 1; ind < static_cast<proshade_unsign>( CSym->size() ); ind++ )
1888 if ( CSym->at(ind)[0] > CSym->at(*symInd)[0] )
1891 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] );
1894 if ( frac >= 1.0 && ( ( CSym->at(*symInd)[5] * 0.85 ) < CSym->at(ind)[5] ) )
1898 ret = CSym->at(ind)[5];
1923 proshade_double ret = 0.0;
1924 proshade_double frac = 0.0;
1925 if ( DSym->size() > 0 )
1927 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] );
1930 else {
return ( ret ); }
1934 for ( proshade_unsign ind = 1; ind < static_cast<proshade_unsign>( DSym->size() ); ind++ )
1937 if ( ( DSym->at(ind)[0] + DSym->at(ind)[6] ) > ( DSym->at(*symInd)[0] + DSym->at(*symInd)[6] ) )
1940 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 );
1943 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] ) )
1947 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] );
1968 proshade_double ret = 0.0;
1969 proshade_double foldSum = 0.0;
1972 if ( TSym->size() == 7 )
1975 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( TSym->size() ); cIt++ )
1977 ret += TSym->at(cIt)[0] * TSym->at(cIt)[5];
1978 foldSum += TSym->at(cIt)[0];
2001 proshade_double ret = 0.0;
2002 proshade_double foldSum = 0.0;
2005 if ( OSym->size() == 13 )
2008 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( OSym->size() ); cIt++ )
2010 ret += OSym->at(cIt)[0] * OSym->at(cIt)[5];
2011 foldSum += OSym->at(cIt)[0];
2034 proshade_double ret = 0.0;
2035 proshade_double foldSum = 0.0;
2038 if ( ISym->size() == 31 )
2041 for ( proshade_unsign cIt = 0; cIt < static_cast<proshade_unsign> ( ISym->size() ); cIt++ )
2043 ret += ISym->at(cIt)[0] * ISym->at(cIt)[5];
2044 foldSum += ISym->at(cIt)[0];
2073 proshade_double cScore = 0.0, dScore = 0.0, tScore = 0.0, oScore = 0.0, iScore = 0.0;
2074 proshade_unsign bestCIndex, bestDIndex;
2077 cScore = this->findBestCScore ( CSym, &bestCIndex );
2078 dScore = this->findBestDScore ( DSym, &bestDIndex );
2079 tScore = this->findTScore ( TSym );
2080 oScore = this->findOScore ( OSym );
2081 iScore = this->findIScore ( ISym );
2084 proshade_double bestWeightedScore = std::max ( cScore, std::max ( dScore * 1.1, std::max ( tScore * 3000.0, std::max ( oScore * 4000.0, iScore * 5000.0 ) ) ) );
2089 if ( bestWeightedScore == cScore )
2097 if ( ( ( 360.0 /
static_cast<double> ( CSym->at(bestCIndex)[0] ) ) - ( 360.0 /
static_cast<double> ( CSym->at(bestCIndex)[0] + 1 ) ) ) <
2098 ( 360.0 /
static_cast<double> ( settings->
maxBandwidth * 4.0 ) ) )
2100 std::stringstream hlpSS;
2101 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).";
2105 if ( bestWeightedScore == dScore * 1.1 )
2108 settings->
setRecommendedFold ( std::max ( DSym->at(bestDIndex)[0], DSym->at(bestDIndex)[6] ) );
2118 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 ) ) ) <
2119 ( 360.0 /
static_cast<double> ( settings->
maxBandwidth * 4.0 ) ) )
2121 std::stringstream hlpSS;
2122 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).";
2126 if ( bestWeightedScore == tScore * 3000.0 )
2131 if ( settings->
detectedSymmetry.size() == 0 ) {
for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( TSym->size() ); it++ ) { settings->
setDetectedSymmetry ( TSym->at(it) ); } }
2133 if ( bestWeightedScore == oScore * 4000.0 )
2138 if ( settings->
detectedSymmetry.size() == 0 ) {
for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( OSym->size() ); it++ ) { settings->
setDetectedSymmetry ( OSym->at(it) ); } }
2140 if ( bestWeightedScore == iScore * 5000.0 )
2145 if ( settings->
detectedSymmetry.size() == 0 ) {
for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( ISym->size() ); it++ ) { settings->
setDetectedSymmetry ( ISym->at(it) ); } }
2167 proshade_unsign bestIndex = 0;
2168 proshade_double highestSym = 0.0;
2171 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( CSym->size() ); iter++ )
2177 if ( CSym->at(iter)[5] > highestSym )
2179 highestSym = CSym->at(iter)[5];
2185 if ( highestSym > 0.0 )
2218 proshade_unsign bestIndex = 0;
2219 proshade_double highestSym = 0.0;
2222 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( DSym->size() ); iter++ )
2225 if ( std::max ( DSym->at(iter)[0], DSym->at(iter)[6] ) != settings->
requestedSymmetryFold ) {
continue; }
2228 if ( ( DSym->at(iter)[5] + DSym->at(iter)[11] ) > highestSym )
2230 highestSym = ( DSym->at(iter)[5] + DSym->at(iter)[11] );
2236 if ( highestSym > 0.0 )
2239 settings->
setRecommendedFold ( std::max ( DSym->at(bestIndex)[0], DSym->at(bestIndex)[6] ) );
2271 std::vector< proshade_double > angList;
2272 std::vector<std::vector< proshade_double > > ret;
2275 proshade_double* rotMat =
new proshade_double[9];
2280 proshade_double normF = std::sqrt( std::pow ( xAx, 2.0 ) + std::pow ( yAx, 2.0 ) + std::pow ( zAx, 2.0 ) );
2286 if ( fold % 2 == 0 )
2289 for ( proshade_double iter =
static_cast < proshade_double
> ( -( ( fold / 2 ) - 1 ) ); iter <= static_cast < proshade_double > ( fold / 2 ); iter++ )
2297 for ( proshade_double iter =
static_cast < proshade_double
> ( -fold / 2 ); iter <= static_cast < proshade_double > ( fold / 2 ); iter++ )
2304 for ( proshade_unsign iter = 0; iter < static_cast < proshade_unsign > ( angList.size() ); iter++ )
2310 std::vector < proshade_double > retEl;
2311 for ( proshade_unsign matIt = 0; matIt < 9; matIt++ )
2335 if ( obtainedAxes != requiredAxes )
2337 std::stringstream hlpSS;
2338 hlpSS <<
"The supplied number of axes for group element\n : detection ( >" << obtainedAxes <<
"< ) does not match the group type ( >" << groupType <<
"< ).";
2339 throw ProSHADE_exception (
"Mismatch between supplied number of axes and\n : symmetry type.",
"ES00059", __FILE__, __LINE__, __func__, hlpSS.str() );
2354 bool checkElementAlreadyExists ( std::vector<std::vector< proshade_double > >* elements, std::vector< proshade_double >* elem, proshade_double matrixTolerance )
2357 bool elementFound =
false;
2360 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( elements->size() ); elIt++ )
2364 elementFound =
true;
2370 return ( elementFound );
2383 bool isGroup =
true;
2386 for ( proshade_unsign gr1 = 0; gr1 < static_cast<proshade_unsign> ( elements->size() ); gr1++ )
2388 for ( proshade_unsign gr2 = 1; gr2 < static_cast<proshade_unsign> ( elements->size() ); gr2++ )
2391 if ( gr1 >= gr2 ) {
continue; }
2405 if ( !isGroup ) {
break; }
2424 std::vector< std::vector< proshade_double > > ret;
2427 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( first->size() ); elIt++ )
2436 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( second->size() ); elIt++ )
2447 for ( proshade_unsign gr1 = 0; gr1 < static_cast<proshade_unsign> ( first->size() ); gr1++ )
2449 for ( proshade_unsign gr2 = 0; gr2 < static_cast<proshade_unsign> ( second->size() ); gr2++ )
2490 std::vector<std::vector< proshade_double > > ret;
2493 if ( groupType ==
"C" )
2508 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." );
2511 else if ( groupType ==
"D" )
2533 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." );
2536 else if ( groupType ==
"T" )
2542 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2559 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2579 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." );
2582 else if ( groupType ==
"O" )
2588 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2605 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2622 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2642 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." );
2645 else if ( groupType ==
"I" )
2651 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2668 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2685 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2705 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." );
2708 else if ( groupType ==
"X" )
2711 for ( proshade_unsign grIt = 0; grIt < static_cast<proshade_unsign> ( axesList.size() ); grIt++ )
2727 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." );
2732 std::stringstream hlpSS;
2733 hlpSS <<
"Unknown symmetry type: >" << groupType <<
"<";
2734 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." );
2752 if ( saveTo != NULL )
2759 saveTo =
new proshade_double[this->xDimIndices * this->yDimIndices * this->zDimIndices];
2765 for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
2767 saveTo[iter] = this->internalMap[iter];
2790 std::stringstream ssHlp;
2796 ssHlp.clear(); ssHlp.str (
"" );
2797 ssHlp <<
" Fold X Y Z Angle Height";
2800 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->
detectedSymmetry.size() ); symIt++ )
2802 ssHlp.clear(); ssHlp.str (
"" );
2807 std::stringstream hlpSS3;
2808 ssHlp.clear(); ssHlp.str (
"" );
2809 hlpSS3 << std::endl <<
"However, since the selection of the recommended symmetry needs improvement, here is a list of all detected C symmetries:";
2814 ssHlp.clear(); ssHlp.str (
"" );
2815 ssHlp <<
" Fold X Y Z Angle Height";
2818 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->
allDetectedCAxes.size() ); symIt++ )
2820 ssHlp.clear(); ssHlp.str (
"" );
2825 hlpSS3.clear(); hlpSS3.str (
"" );
2826 hlpSS3 << std::endl <<
"Also, for the same reason, here is a list of all detected D symmetries:";
2831 ssHlp.clear(); ssHlp.str (
"" );
2832 ssHlp <<
" Fold X Y Z Angle Height";
2835 for ( proshade_unsign symIt = 0; symIt < static_cast<proshade_unsign> ( settings->
allDetectedDAxes.size() ); symIt++ )
2837 ssHlp.clear(); ssHlp.str (
"" );
2841 for ( proshade_unsign axIt = 1; axIt < static_cast<proshade_unsign> ( settings->
allDetectedDAxes.at(symIt).size() ); axIt++ )
2843 ssHlp.clear(); ssHlp.str (
"" );
2848 ssHlp.clear(); ssHlp.str (
"" );
2869 proshade_double totNonZeroPoints = 0.0;
2870 proshade_signed mapIt = 0;
2873 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
2875 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
2877 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
2880 mapIt = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
2883 if ( this->internalMap[mapIt] <= 0.0 ) {
continue; }
2886 this->xCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( xIt + this->xFrom );
2887 this->yCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( yIt + this->yFrom );
2888 this->zCom += this->internalMap[mapIt] *
static_cast<proshade_double
> ( zIt + this->zFrom );
2889 totNonZeroPoints += this->internalMap[mapIt];
2894 this->xCom /= totNonZeroPoints;
2895 this->yCom /= totNonZeroPoints;
2896 this->zCom /= totNonZeroPoints;
2899 this->xCom = (
static_cast<proshade_double
> ( this->xFrom ) * ( this->xDimSizeOriginal /
static_cast<proshade_double
> ( this->xDimIndicesOriginal ) ) ) +
2900 ( ( this->xCom -
static_cast<proshade_double
> ( this->xFrom ) ) *
2901 ( this->xDimSizeOriginal /
static_cast<proshade_double
> ( this->xDimIndicesOriginal ) ) );
2902 this->yCom = (
static_cast<proshade_double
> ( this->yFrom ) * ( this->yDimSizeOriginal /
static_cast<proshade_double
> ( this->yDimIndicesOriginal ) ) ) +
2903 ( ( this->yCom -
static_cast<proshade_double
> ( this->yFrom ) ) *
2904 ( this->yDimSizeOriginal /
static_cast<proshade_double
> ( this->yDimIndicesOriginal ) ) );
2905 this->zCom = (
static_cast<proshade_double
> ( this->zFrom ) * ( this->zDimSizeOriginal /
static_cast<proshade_double
> ( this->zDimIndicesOriginal ) ) ) +
2906 ( ( this->zCom -
static_cast<proshade_double
> ( this->zFrom ) ) *
2907 ( this->zDimSizeOriginal /
static_cast<proshade_double
> ( this->zDimIndicesOriginal ) ) );
2921 return ( this->noSpheres );
2932 return ( this->internalMap[pos] );
2942 return ( this->maxShellBand );
2952 return ( this->rrpMatrices[band][sh1][sh2] );
2967 if ( this->spheres[shell]->getLocalBandwidth( ) >= bandVal )
2991 fftw_complex* mapCoeffs =
new fftw_complex[this->xDimIndices * this->yDimIndices * this->zDimIndices];
2992 fftw_complex* pattersonMap =
new fftw_complex[this->xDimIndices * this->yDimIndices * this->zDimIndices];
2999 for ( proshade_unsign iter = 0; iter < (this->xDimIndices * this->yDimIndices * this->zDimIndices); iter++ )
3001 pattersonMap[iter][0] = this->internalMap[iter];
3002 pattersonMap[iter][1] = 0.0;
3006 fftw_plan forward = fftw_plan_dft_3d ( this->xDimIndices, this->yDimIndices, this->zDimIndices,
3007 pattersonMap, mapCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
3008 fftw_plan inverse = fftw_plan_dft_3d ( this->xDimIndices, this->yDimIndices, this->zDimIndices,
3009 mapCoeffs, pattersonMap, FFTW_BACKWARD, FFTW_ESTIMATE );
3012 fftw_execute ( forward );
3018 fftw_execute ( inverse );
3021 proshade_signed mapIt, patIt, patX, patY, patZ;
3022 for ( proshade_signed xIt = 0; xIt < static_cast<proshade_signed> ( this->xDimIndices ); xIt++ )
3024 for ( proshade_signed yIt = 0; yIt < static_cast<proshade_signed> ( this->yDimIndices ); yIt++ )
3026 for ( proshade_signed zIt = 0; zIt < static_cast<proshade_signed> ( this->zDimIndices ); zIt++ )
3029 patX = xIt - (
static_cast<proshade_signed
> ( this->xDimIndices ) / 2 );
if ( patX < 0 ) { patX += this->xDimIndices; }
3030 patY = yIt - (
static_cast<proshade_signed
> ( this->yDimIndices ) / 2 );
if ( patY < 0 ) { patY += this->yDimIndices; }
3031 patZ = zIt - (
static_cast<proshade_signed
> ( this->zDimIndices ) / 2 );
if ( patZ < 0 ) { patZ += this->zDimIndices; }
3034 mapIt = zIt + this->zDimIndices * ( yIt + this->yDimIndices * xIt );
3035 patIt = patZ + this->zDimIndices * ( patY + this->yDimIndices * patX );
3038 this->internalMap[mapIt] = pattersonMap[patIt][0];
3044 delete[] pattersonMap;
3048 fftw_destroy_plan ( forward );
3049 fftw_destroy_plan ( inverse );
3066 return ( &this->sphericalHarmonics[shell][seanindex (
static_cast<proshade_signed
> ( order ) -
static_cast<proshade_signed
> ( band ),
3068 this->spheres[shell]->getLocalBandwidth() )][0] );
3079 return ( &this->sphericalHarmonics[shell][seanindex (
static_cast<proshade_signed
> ( order ) -
static_cast<proshade_signed
> ( band ),
3081 this->spheres[shell]->getLocalBandwidth() )][1] );
3092 return ( this->spheres[shell]->getShellRadius() );
3103 return ( this->integrationWeight );
3115 return ( this->spheres[shell]->getLocalBandwidth ( ) );
3127 return ( this->spherePos.at(shell) );
3139 return ( this->eMatrices[band] );
3154 *valueReal = this->eMatrices[band][order1][order2][0];
3155 *valueImag = this->eMatrices[band][order1][order2][1];
3169 return ( this->so3CoeffsInverse );
3180 return ( this->so3Coeffs );
3191 return ( this->maxCompBand );
3206 *valueReal = this->wignerMatrices[band][order1][order2][0];
3207 *valueImag = this->wignerMatrices[band][order1][order2][1];
3221 return ( this->xDimSize );
3231 return ( this->yDimSize );
3241 return ( this->zDimSize );
3251 return ( this->xDimIndices );
3261 return ( this->yDimIndices );
3271 return ( this->zDimIndices );
3281 return ( &this->xFrom );
3291 return ( &this->yFrom );
3301 return ( &this->zFrom );
3311 return ( &this->xTo );
3321 return ( &this->yTo );
3331 return ( &this->zTo );
3341 return ( &this->xAxisOrigin );
3351 return ( &this->yAxisOrigin );
3361 return ( &this->zAxisOrigin );
3371 return ( this->internalMap );
3381 return ( this->translationMap );
3391 this->integrationWeight = intW;
3405 this->integrationWeight += intW;
3422 this->eMatrices[band][order1][order2][0] = val[0];
3423 this->eMatrices[band][order1][order2][1] = val[1];
3440 this->eMatrices[band][order1][order2][0] /= normF;
3441 this->eMatrices[band][order1][order2][1] /= normF;
3456 this->so3Coeffs[position][0] = val[0];
3457 this->so3Coeffs[position][1] = val[1];
3474 this->wignerMatrices[band][order1][order2][0] = val[0];
3475 this->wignerMatrices[band][order1][order2][1] = val[1];
3492 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3494 eMatsLMReal[iter] =
static_cast<double> ( this->eMatrices[band][order1][iter][0] );
3512 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3514 eMatsLMImag[iter] =
static_cast<double> ( this->eMatrices[band][order1][iter][1] );
3530 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3532 so3CoefsReal[iter] =
static_cast<double> ( this->so3Coeffs[iter][0] );
3548 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3550 so3CoefsImag[iter] =
static_cast<double> ( this->so3Coeffs[iter][1] );
3570 return (
static_cast<int> ( so3CoefLoc ( order1, order2, band, this->getMaxBand() ) ) );
3581 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3583 rotFunReal[iter] =
static_cast<double> ( this->so3CoeffsInverse[iter][0] );
3599 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3601 rotFunImag[iter] =
static_cast<double> ( this->so3CoeffsInverse[iter][1] );
3617 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3619 trsFunReal[iter] =
static_cast<double> ( this->translationMap[iter][0] );
3635 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3637 trsFunImag[iter] =
static_cast<double> ( this->translationMap[iter][1] );
3656 proshade_double eA, eB, eG;
3660 proshade_double* rMat = NULL;
3661 rMat =
new proshade_double[9];
3668 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( len ); iter++ )
3670 rotMat[iter] =
static_cast<double> ( rMat[iter] );
3711 return (
static_cast<proshade_unsign
> ( settings->
detectedSymmetry.size() ) );
3723 if (
static_cast<proshade_unsign
> ( settings->
detectedSymmetry.size() ) <= axisNo )
3726 return ( std::vector< std::string > ( ) );
3730 std::vector< std::string > ret;
3733 std::stringstream ssHlp;
3778 std::stringstream fNameHlp;
3780 this->writeMap ( fNameHlp.str() );
3787 this->writePdb ( fNameHlp.str(), eulA, eulB, eulG, ultimateTranslation->at(0), ultimateTranslation->at(1), ultimateTranslation->at(2), settings->
firstModelOnly );
3793 ultimateTranslation->at(0), ultimateTranslation->at(1), ultimateTranslation->at(2),
3815 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);
3819 proshade_double* rotMat =
new proshade_double[9];
3823 std::stringstream rotMatSS;
3824 rotMatSS << std::setprecision (3) << std::showpos <<
"The rotation matrix about origin is : " << rotMat[0] <<
" " << rotMat[1] <<
" " << rotMat[2] << std::endl;
3825 rotMatSS << std::setprecision (3) << std::showpos <<
" : " << rotMat[3] <<
" " << rotMat[4] <<
" " << rotMat[5] << std::endl;
3826 rotMatSS << std::setprecision (3) << std::showpos <<
" : " << rotMat[6] <<
" " << rotMat[7] <<
" " << rotMat[8];
3832 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);