ProSHADE  0.7.6.1 (AUG 2021)
Protein Shape Detection
ProSHADE_overlay.cpp
Go to the documentation of this file.
1 
22 //==================================================== ProSHADE
23 #include "ProSHADE_overlay.hpp"
24 
36 {
37  //================================================ Report progress
38  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting rotation function computation." );
39 
40  //================================================ Compute un-weighted E matrices and their weights
41  ProSHADE_internal_distances::computeEMatrices ( obj2, this, settings );
42 
43  //================================================ Normalise E matrices by the magnitudes
44  ProSHADE_internal_distances::normaliseEMatrices ( obj2, this, settings );
45 
46  //================================================ Generate SO(3) coefficients
48 
49  //================================================ Compute the inverse SO(3) Fourier Transform (SOFT) on the newly computed coefficients
51 
52  //================================================ Report completion
53  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Rotation function obtained." );
54 
55  //================================================ Done
56  return ;
57 
58 }
59 
74 void ProSHADE_internal_overlay::getOptimalRotation ( ProSHADE_settings* settings, ProSHADE_internal_data::ProSHADE_data* staticStructure, ProSHADE_internal_data::ProSHADE_data* movingStructure, proshade_double* eulA, proshade_double* eulB, proshade_double* eulG )
75 {
76  //================================================ Read in the structures
77  staticStructure->readInStructure ( settings->inputFiles.at(0), 0, settings );
78  movingStructure->readInStructure ( settings->inputFiles.at(1), 1, settings );
79 
80  //================================================ Internal data processing (COM, norm, mask, extra space)
81  staticStructure->processInternalMap ( settings );
82  movingStructure->processInternalMap ( settings );
83 
84  //================================================ Map to sphere
85  staticStructure->mapToSpheres ( settings );
86  movingStructure->mapToSpheres ( settings );
87 
88  //================================================ Get spherical harmonics
89  staticStructure->computeSphericalHarmonics ( settings );
90  movingStructure->computeSphericalHarmonics ( settings );
91 
92  //================================================ Get the rotation function of the pair
93  movingStructure->getOverlayRotationFunction ( settings, staticStructure );
94 
95  //================================================ Get inverse SO(3) map top peak Euler angle values
97  std::min ( staticStructure->getMaxBand(), movingStructure->getMaxBand() ) * 2,
98  eulA, eulB, eulG, settings );
99 
100  //================================================ Done
101  return ;
102 
103 }
104 
123 void ProSHADE_internal_overlay::getOptimalTranslation ( ProSHADE_settings* settings, ProSHADE_internal_data::ProSHADE_data* staticStructure, ProSHADE_internal_data::ProSHADE_data* movingStructure, proshade_double* trsX, proshade_double* trsY, proshade_double* trsZ, proshade_double eulA, proshade_double eulB, proshade_double eulG )
124 {
125  //================================================ Read in the structures
126  staticStructure->readInStructure ( settings->inputFiles.at(0), 0, settings );
127  movingStructure->readInStructure ( settings->inputFiles.at(1), 1, settings );
128 
129  //================================================ Determine spherical harmonics variables from the first structure (otherwise, they would be determined from the second structure)
130  settings->determineAllSHValues ( staticStructure->xDimIndices, staticStructure->yDimIndices,
131  staticStructure->xDimSize, staticStructure->yDimSize, staticStructure->zDimSize );
132 
133  //================================================ Internal data processing (COM, norm, mask, extra space)
134  staticStructure->processInternalMap ( settings );
135  movingStructure->processInternalMap ( settings );
136 
137  //================================================ Rotate map
138  std::vector< proshade_double > rotCen = movingStructure->rotateMapRealSpaceInPlace ( eulA, eulB, eulG );
139 
140  //================================================ Zero padding for smaller structure
141  staticStructure->zeroPaddToDims ( std::max ( staticStructure->getXDim(), movingStructure->getXDim() ),
142  std::max ( staticStructure->getYDim(), movingStructure->getYDim() ),
143  std::max ( staticStructure->getZDim(), movingStructure->getZDim() ) );
144  movingStructure->zeroPaddToDims ( std::max ( staticStructure->getXDim(), movingStructure->getXDim() ),
145  std::max ( staticStructure->getYDim(), movingStructure->getYDim() ),
146  std::max ( staticStructure->getZDim(), movingStructure->getZDim() ) );
147 
148  //================================================ Report progress
149  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 1, "Starting translation function computation." );
150 
151  //================================================ Compute the translation function map
152  movingStructure->computeTranslationMap ( staticStructure );
153 
154  //================================================ Find highest peak in translation function
155  proshade_double mapPeak = 0.0;
156  proshade_unsign xDimS = staticStructure->getXDim();
157  proshade_unsign yDimS = staticStructure->getYDim();
158  proshade_unsign zDimS = staticStructure->getZDim();
159  ProSHADE_internal_overlay::findHighestValueInMap ( movingStructure->getTranslationFnPointer(), xDimS, yDimS, zDimS, trsX, trsY, trsZ, &mapPeak );
160 
161  //================================================ Dont translate over half
162  if ( *trsX > ( static_cast< proshade_double > ( xDimS ) / 2.0 ) ) { *trsX = *trsX - static_cast< proshade_double > ( xDimS ); }
163  if ( *trsY > ( static_cast< proshade_double > ( yDimS ) / 2.0 ) ) { *trsY = *trsY - static_cast< proshade_double > ( yDimS ); }
164  if ( *trsZ > ( static_cast< proshade_double > ( zDimS ) / 2.0 ) ) { *trsZ = *trsZ - static_cast< proshade_double > ( zDimS ); }
165 
166  //================================================ Convert the map positions onto translation in Angstroms and save the results
167  computeTranslationsFromPeak ( staticStructure, movingStructure, trsX, trsY, trsZ );
168 
169  movingStructure->originalPdbTransX = *trsX;
170  movingStructure->originalPdbTransY = *trsY;
171  movingStructure->originalPdbTransZ = *trsZ;
172 
173  proshade_single xMov = static_cast< proshade_single > ( -(*trsX) );
174  proshade_single yMov = static_cast< proshade_single > ( -(*trsY) );
175  proshade_single zMov = static_cast< proshade_single > ( -(*trsZ) );
176 
177  //================================================ Report progress
178  std::stringstream hlpSS;
179  hlpSS << "Optimal map translation distances are " << *trsX << " ; " << *trsY << " ; " << *trsZ << " Angstroms with peak height " << mapPeak / ( static_cast< proshade_double > ( xDimS ) * static_cast< proshade_double > ( yDimS ) * static_cast< proshade_double > ( zDimS ) );
180 
181  //================================================ Save original from variables for PDB output
182  movingStructure->mapMovFromsChangeX = static_cast< proshade_double > ( movingStructure->xFrom );
183  movingStructure->mapMovFromsChangeY = static_cast< proshade_double > ( movingStructure->yFrom );
184  movingStructure->mapMovFromsChangeZ = static_cast< proshade_double > ( movingStructure->zFrom );
185 
186  //================================================ Move the map
187  ProSHADE_internal_mapManip::moveMapByIndices ( &xMov, &yMov, &zMov, movingStructure->getXDimSize(), movingStructure->getYDimSize(), movingStructure->getZDimSize(),
188  movingStructure->getXFromPtr(), movingStructure->getXToPtr(),
189  movingStructure->getYFromPtr(), movingStructure->getYToPtr(),
190  movingStructure->getZFromPtr(), movingStructure->getZToPtr(),
191  movingStructure->getXAxisOrigin(), movingStructure->getYAxisOrigin(), movingStructure->getZAxisOrigin() );
192 
193  ProSHADE_internal_mapManip::moveMapByFourier ( movingStructure->getInternalMap(), xMov, yMov, zMov,
194  movingStructure->getXDimSize(), movingStructure->getYDimSize(), movingStructure->getZDimSize(),
195  static_cast< proshade_signed > ( movingStructure->getXDim() ), static_cast< proshade_signed > ( movingStructure->getYDim() ),
196  static_cast< proshade_signed > ( movingStructure->getZDim() ) );
197 
198  //================================================ Report progress
199  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 3, hlpSS.str() );
200  ProSHADE_internal_messages::printProgressMessage ( settings->verbose, 2, "Translation function computation complete." );
201 
202  //================================================ Keep only the change in from and to variables
203  movingStructure->mapMovFromsChangeX = static_cast< proshade_double > ( movingStructure->xFrom ) - movingStructure->mapMovFromsChangeX;
204  movingStructure->mapMovFromsChangeY = static_cast< proshade_double > ( movingStructure->yFrom ) - movingStructure->mapMovFromsChangeY;
205  movingStructure->mapMovFromsChangeZ = static_cast< proshade_double > ( movingStructure->zFrom ) - movingStructure->mapMovFromsChangeZ;
206 
207  //================================================ Done
208  return ;
209 
210 }
211 
220 void ProSHADE_internal_overlay::computeTranslationsFromPeak ( ProSHADE_internal_data::ProSHADE_data* staticStructure, ProSHADE_internal_data::ProSHADE_data* movingStructure, proshade_double *trsX, proshade_double *trsY, proshade_double *trsZ )
221 {
222  //================================================ Compute map movement
223  proshade_single xSamplRate = ( static_cast< proshade_single > ( staticStructure->getXDimSize() ) / static_cast< proshade_single > ( staticStructure->getXDim() ) );
224  proshade_single ySamplRate = ( static_cast< proshade_single > ( staticStructure->getYDimSize() ) / static_cast< proshade_single > ( staticStructure->getYDim() ) );
225  proshade_single zSamplRate = ( static_cast< proshade_single > ( staticStructure->getZDimSize() ) / static_cast< proshade_single > ( staticStructure->getZDim() ) );
226 
227  proshade_single xCor = static_cast< proshade_single > ( static_cast< proshade_single > ( staticStructure->xFrom - movingStructure->xFrom ) * xSamplRate );
228  proshade_single yCor = static_cast< proshade_single > ( static_cast< proshade_single > ( staticStructure->yFrom - movingStructure->yFrom ) * ySamplRate );
229  proshade_single zCor = static_cast< proshade_single > ( static_cast< proshade_single > ( staticStructure->zFrom - movingStructure->zFrom ) * zSamplRate );
230 
231  proshade_single xMov = static_cast< proshade_single > ( 0.0 - static_cast< proshade_double > ( xCor ) - ( *trsX * static_cast< proshade_double > ( xSamplRate ) ) );
232  proshade_single yMov = static_cast< proshade_single > ( 0.0 - static_cast< proshade_double > ( yCor ) - ( *trsY * static_cast< proshade_double > ( ySamplRate ) ) );
233  proshade_single zMov = static_cast< proshade_single > ( 0.0 - static_cast< proshade_double > ( zCor ) - ( *trsZ * static_cast< proshade_double > ( zSamplRate ) ) );
234 
235  //================================================ Save translation vector back
236  *trsX = static_cast< proshade_double > ( -xMov );
237  *trsY = static_cast< proshade_double > ( -yMov );
238  *trsZ = static_cast< proshade_double > ( -zMov );
239 
240  //================================================ Done
241  return ;
242 
243 }
244 
251 {
252  //================================================ Initialise local variables
253  std::vector< proshade_double > ret;
254  proshade_double mapPeak = 0.0;
255  proshade_unsign xDimS = staticStructure->getXDim();
256  proshade_unsign yDimS = staticStructure->getYDim();
257  proshade_unsign zDimS = staticStructure->getZDim();
258  proshade_double trsX = 0.0, trsY = 0.0, trsZ = 0.0;
259  ProSHADE_internal_overlay::findHighestValueInMap ( this->getTranslationFnPointer(),
260  xDimS,
261  yDimS,
262  zDimS,
263  &trsX,
264  &trsY,
265  &trsZ,
266  &mapPeak );
267 
268  //================================================ Dont translate over half
269  if ( trsX > ( static_cast< proshade_double > ( xDimS ) / 2.0 ) ) { trsX = trsX - static_cast< proshade_double > ( xDimS ); }
270  if ( trsY > ( static_cast< proshade_double > ( yDimS ) / 2.0 ) ) { trsY = trsY - static_cast< proshade_double > ( yDimS ); }
271  if ( trsZ > ( static_cast< proshade_double > ( zDimS ) / 2.0 ) ) { trsZ = trsZ - static_cast< proshade_double > ( zDimS ); }
272 
273  //================================================ Convert the map positions onto translation in Angstroms and save the results
274  ProSHADE_internal_overlay::computeTranslationsFromPeak ( staticStructure, this, &trsX, &trsY, &trsZ );
275 
276  this->originalPdbTransX = trsX;
277  this->originalPdbTransY = trsY;
278  this->originalPdbTransZ = trsZ;
279 
280  proshade_single xMov = static_cast< proshade_single > ( -trsX );
281  proshade_single yMov = static_cast< proshade_single > ( -trsY );
282  proshade_single zMov = static_cast< proshade_single > ( -trsZ );
283 
284  //================================================ Save results as vector
285  ProSHADE_internal_misc::addToDoubleVector ( &ret, static_cast<proshade_double> ( -xMov ) );
286  ProSHADE_internal_misc::addToDoubleVector ( &ret, static_cast<proshade_double> ( -yMov ) );
287  ProSHADE_internal_misc::addToDoubleVector ( &ret, static_cast<proshade_double> ( -zMov ) );
288 
289  //================================================ Save original from variables for PDB output
290  this->mapMovFromsChangeX = static_cast<proshade_double> ( this->xFrom );
291  this->mapMovFromsChangeY = static_cast<proshade_double> ( this->yFrom );
292  this->mapMovFromsChangeZ = static_cast<proshade_double> ( this->zFrom );
293 
294  //================================================ Move the map
295  ProSHADE_internal_mapManip::moveMapByIndices ( &xMov, &yMov, &zMov, this->getXDimSize(), this->getYDimSize(), this->getZDimSize(),
296  this->getXFromPtr(), this->getXToPtr(),
297  this->getYFromPtr(), this->getYToPtr(),
298  this->getZFromPtr(), this->getZToPtr(),
299  this->getXAxisOrigin(), this->getYAxisOrigin(), this->getZAxisOrigin() );
300 
301  ProSHADE_internal_mapManip::moveMapByFourier ( this->getInternalMap(), xMov, yMov, zMov,
302  this->getXDimSize(), this->getYDimSize(), this->getZDimSize(),
303  static_cast< proshade_signed > ( this->getXDim() ), static_cast< proshade_signed > ( this->getYDim() ), static_cast< proshade_signed > ( this->getZDim() ) );
304 
305  //================================================ Keep only the change in from and to variables
306  this->mapMovFromsChangeX = static_cast<proshade_double> ( this->xFrom ) - this->mapMovFromsChangeX;
307  this->mapMovFromsChangeY = static_cast<proshade_double> ( this->yFrom ) - this->mapMovFromsChangeY;
308  this->mapMovFromsChangeZ = static_cast<proshade_double> ( this->zFrom ) - this->mapMovFromsChangeZ;
309 
310  //================================================ Done
311  return ( ret );
312 
313 }
314 
325 {
326  //================================================ Do this using Fourier!
327  fftw_complex *tmpIn1 = nullptr, *tmpOut1 = nullptr, *tmpIn2 = nullptr, *tmpOut2 = nullptr, *resOut = nullptr;
328  fftw_plan forwardFourierObj1, forwardFourierObj2, inverseFourierCombo;
329  proshade_unsign dimMult = staticStructure->getXDim() * staticStructure->getYDim() * staticStructure->getZDim();
330  ProSHADE_internal_overlay::allocateTranslationFunctionMemory ( tmpIn1, tmpOut1, tmpIn2, tmpOut2, this->translationMap, resOut, forwardFourierObj1, forwardFourierObj2, inverseFourierCombo, staticStructure->getXDim(), staticStructure->getYDim(), staticStructure->getZDim() );
331 
332  //================================================ Fill in input data
333  for ( proshade_unsign iter = 0; iter < dimMult; iter++ ) { tmpIn1[iter][0] = staticStructure->getMapValue ( iter ); tmpIn1[iter][1] = 0.0; }
334  for ( proshade_unsign iter = 0; iter < dimMult; iter++ ) { tmpIn2[iter][0] = this->getMapValue ( iter ); tmpIn2[iter][1] = 0.0; }
335 
336  //================================================ Calculate Fourier
337  fftw_execute ( forwardFourierObj1 );
338  fftw_execute ( forwardFourierObj2 );
339 
340  //================================================ Combine Fourier coeffs and invert
341  ProSHADE_internal_overlay::combineFourierForTranslation ( tmpOut1, tmpOut2, resOut, staticStructure->getXDim(), staticStructure->getYDim(), staticStructure->getZDim() );
342  fftw_execute ( inverseFourierCombo );
343 
344  //================================================ Free memory
345  ProSHADE_internal_overlay::freeTranslationFunctionMemory ( tmpIn1, tmpOut1, tmpIn2, tmpOut2, resOut, forwardFourierObj1, forwardFourierObj2, inverseFourierCombo );
346 
347  //================================================ Done
348  return ;
349 
350 }
351 
367 void ProSHADE_internal_overlay::allocateTranslationFunctionMemory ( fftw_complex*& tmpIn1, fftw_complex*& tmpOut1, fftw_complex*& tmpIn2, fftw_complex*& tmpOut2, fftw_complex*& resIn, fftw_complex*& resOut, fftw_plan& forwardFourierObj1, fftw_plan& forwardFourierObj2, fftw_plan& inverseFourierCombo, proshade_unsign xD, proshade_unsign yD, proshade_unsign zD )
368 {
369  //================================================ Allocate memory
370  tmpIn1 = new fftw_complex[xD * yD * zD];
371  tmpOut1 = new fftw_complex[xD * yD * zD];
372  tmpIn2 = new fftw_complex[xD * yD * zD];
373  tmpOut2 = new fftw_complex[xD * yD * zD];
374  resIn = new fftw_complex[xD * yD * zD];
375  resOut = new fftw_complex[xD * yD * zD];
376 
377  //================================================ Check memory allocation
378  ProSHADE_internal_misc::checkMemoryAllocation ( tmpIn1 , __FILE__, __LINE__, __func__ );
379  ProSHADE_internal_misc::checkMemoryAllocation ( tmpOut1, __FILE__, __LINE__, __func__ );
380  ProSHADE_internal_misc::checkMemoryAllocation ( tmpIn2, __FILE__, __LINE__, __func__ );
381  ProSHADE_internal_misc::checkMemoryAllocation ( tmpOut2, __FILE__, __LINE__, __func__ );
382  ProSHADE_internal_misc::checkMemoryAllocation ( resIn, __FILE__, __LINE__, __func__ );
383  ProSHADE_internal_misc::checkMemoryAllocation ( resOut, __FILE__, __LINE__, __func__ );
384 
385  //================================================ Get Fourier transforms of the maps
386  forwardFourierObj1 = fftw_plan_dft_3d ( static_cast< int > ( xD ), static_cast< int > ( yD ), static_cast< int > ( zD ), tmpIn1, tmpOut1, FFTW_FORWARD , FFTW_ESTIMATE );
387  forwardFourierObj2 = fftw_plan_dft_3d ( static_cast< int > ( xD ), static_cast< int > ( yD ), static_cast< int > ( zD ), tmpIn2, tmpOut2, FFTW_FORWARD , FFTW_ESTIMATE );
388  inverseFourierCombo = fftw_plan_dft_3d ( static_cast< int > ( xD ), static_cast< int > ( yD ), static_cast< int > ( zD ), resOut, resIn , FFTW_BACKWARD, FFTW_ESTIMATE );
389 
390  //================================================ Done
391  return ;
392 
393 }
394 
406 void ProSHADE_internal_overlay::freeTranslationFunctionMemory ( fftw_complex*& tmpIn1, fftw_complex*& tmpOut1, fftw_complex*& tmpIn2, fftw_complex*& tmpOut2, fftw_complex*& resOut, fftw_plan& forwardFourierObj1, fftw_plan& forwardFourierObj2, fftw_plan& inverseFourierCombo )
407 {
408  //================================================ Release memory
409  fftw_destroy_plan ( forwardFourierObj1 );
410  fftw_destroy_plan ( forwardFourierObj2 );
411  fftw_destroy_plan ( inverseFourierCombo );
412  delete[] tmpIn1;
413  delete[] tmpIn2;
414  delete[] tmpOut1;
415  delete[] tmpOut2;
416  delete[] resOut;
417 
418  //======================================== Done
419  return ;
420 
421 }
422 
432 void ProSHADE_internal_overlay::combineFourierForTranslation ( fftw_complex* tmpOut1, fftw_complex* tmpOut2, fftw_complex*& resOut, proshade_unsign xD, proshade_unsign yD, proshade_unsign zD )
433 {
434  //================================================ Initialise local variables
435  double normFactor = static_cast<double> ( xD * yD * zD );
436  proshade_signed arrPos;
437 
438  //================================================ Combine the coefficients
439  for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xD ); xIt++ )
440  {
441  for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yD ); yIt++ )
442  {
443  for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zD ); zIt++ )
444  {
445  //==================================== Find indices
446  arrPos = zIt + static_cast< proshade_signed > ( zD ) * ( yIt + static_cast< proshade_signed > ( yD ) * xIt );
447 
448  //==================================== Combine
450  &tmpOut1[arrPos][1],
451  &tmpOut2[arrPos][0],
452  &tmpOut2[arrPos][1],
453  &resOut[arrPos][0],
454  &resOut[arrPos][1] );
455 
456  //==================================== Save
457  resOut[arrPos][0] /= normFactor;
458  resOut[arrPos][1] /= normFactor;
459  }
460  }
461  }
462 
463  //================================================ Done
464  return ;
465 
466 }
467 
479 void ProSHADE_internal_overlay::findHighestValueInMap ( fftw_complex* resIn, proshade_unsign xD, proshade_unsign yD, proshade_unsign zD, proshade_double* trsX, proshade_double* trsY, proshade_double* trsZ, proshade_double* mapPeak )
480 {
481  //================================================ Initialise variables
482  proshade_signed arrPos;
483  *mapPeak = 0.0;
484 
485  //================================================ Search the map
486  for ( proshade_signed uIt = 0; uIt < static_cast<proshade_signed> ( xD ); uIt++ )
487  {
488  for ( proshade_signed vIt = 0; vIt < static_cast<proshade_signed> ( yD ); vIt++ )
489  {
490  for ( proshade_signed wIt = 0; wIt < static_cast<proshade_signed> ( zD ); wIt++ )
491  {
492  arrPos = wIt + static_cast< proshade_signed > ( zD ) * ( vIt + static_cast< proshade_signed > ( yD ) * uIt );
493  if ( resIn[arrPos][0] > *mapPeak )
494  {
495  *mapPeak = resIn[arrPos][0];
496  *trsX = static_cast< proshade_double > ( uIt );
497  *trsY = static_cast< proshade_double > ( vIt );
498  *trsZ = static_cast< proshade_double > ( wIt );
499  }
500  }
501  }
502  }
503 
504  //================================================ Done
505  return ;
506 
507 }
508 
519 void ProSHADE_internal_data::ProSHADE_data::zeroPaddToDims ( proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim )
520 {
521  //================================================ Sanity check
522  if ( ( this->xDimIndices > xDim ) || ( this->yDimIndices > yDim ) || ( this->zDimIndices > zDim ) )
523  {
524  throw ProSHADE_exception ( "Cannot zero-pad in negative direction.", "EO00034", __FILE__, __LINE__, __func__, "The requested padded size of a structure is smaller than\n : the current size. If the user sees this error, there is\n : likely a considerable bug. Please report this error." );
525  }
526 
527  //================================================ If done, do nothing
528  if ( ( this->xDimIndices == xDim ) && ( this->yDimIndices == yDim ) && ( this->zDimIndices == zDim ) ) { return ; }
529 
530  //================================================ Find out how many zeroes need to be added before and after
531  proshade_unsign addXPre, addYPre, addZPre, addXPost, addYPost, addZPost;
532  ProSHADE_internal_overlay::computeBeforeAfterZeroCounts ( &addXPre, &addYPre, &addZPre, &addXPost, &addYPost, &addZPost, xDim, yDim, zDim, this->xDimIndices, this->yDimIndices, this->zDimIndices );
533 
534  //================================================ Create a new map
535  proshade_double* newMap = new proshade_double [xDim * yDim * zDim];
536 
537  //================================================ Do the hard work
538  ProSHADE_internal_overlay::paddMapWithZeroes ( this->internalMap, newMap, xDim, yDim, zDim, this->xDimIndices, this->yDimIndices, this->zDimIndices, addXPre, addYPre, addZPre );
539 
540  //================================================ Create a new internal map and copy
541  delete[] this->internalMap;
542  this->internalMap = new proshade_double [xDim * yDim * zDim];
543  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { this->internalMap[iter] = newMap[iter]; }
544 
545  //================================================ Release memory
546  delete[] newMap;
547 
548  //================================================ Update map related variables
549  this->xDimSize = static_cast< proshade_single > ( xDim ) * ( this->xDimSize / static_cast< proshade_single > ( this->xDimIndices ) );
550  this->yDimSize = static_cast< proshade_single > ( yDim ) * ( this->yDimSize / static_cast< proshade_single > ( this->yDimIndices ) );
551  this->zDimSize = static_cast< proshade_single > ( zDim ) * ( this->zDimSize / static_cast< proshade_single > ( this->zDimIndices ) );
552  this->xDimIndices = xDim ; this->yDimIndices = yDim ; this->zDimIndices = zDim;
553  this->xGridIndices = xDim ; this->yGridIndices = yDim ; this->zGridIndices = zDim;
554  this->xFrom -= addXPre ; this->yFrom -= addYPre ; this->zFrom -= addZPre;
555  this->xTo += addXPost; this->yTo += addYPost; this->zTo += addZPost;
556  this->xAxisOrigin -= addXPre ; this->yAxisOrigin -= addYPre ; this->zAxisOrigin -= addZPre ;
557 
558  //================================================ Done
559  return ;
560 
561 }
562 
578 void ProSHADE_internal_overlay::computeBeforeAfterZeroCounts ( proshade_unsign* addXPre, proshade_unsign* addYPre, proshade_unsign* addZPre, proshade_unsign* addXPost, proshade_unsign* addYPost, proshade_unsign* addZPost, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_unsign xDimIndices, proshade_unsign yDimIndices, proshade_unsign zDimIndices )
579 {
580  //================================================ Compute
581  *addXPre = ( xDim - xDimIndices ) / 2;
582  *addYPre = ( yDim - yDimIndices ) / 2;
583  *addZPre = ( zDim - zDimIndices ) / 2;
584  *addXPost = *addXPre; if ( ( xDim - xDimIndices ) % 2 == 1 ) { *addXPost += 1; }
585  *addYPost = *addYPre; if ( ( yDim - yDimIndices ) % 2 == 1 ) { *addYPost += 1; }
586  *addZPost = *addZPre; if ( ( zDim - zDimIndices ) % 2 == 1 ) { *addZPost += 1; }
587 
588  //================================================ Done
589  return ;
590 
591 }
592 
607 void ProSHADE_internal_overlay::paddMapWithZeroes ( proshade_double* oldMap, proshade_double*& newMap, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_unsign xDimIndices, proshade_unsign yDimIndices, proshade_unsign zDimIndices, proshade_unsign addXPre, proshade_unsign addYPre, proshade_unsign addZPre )
608 {
609  //================================================ Initialise local variables
610  proshade_unsign newMapIndex = 0;
611  proshade_unsign oldMapIndex = 0;
612 
613  //================================================ Copy and padd as appropriate
614  for ( proshade_unsign xIt = 0; xIt < xDim; xIt++ )
615  {
616  for ( proshade_unsign yIt = 0; yIt < yDim; yIt++ )
617  {
618  for ( proshade_unsign zIt = 0; zIt < zDim; zIt++ )
619  {
620  //==================================== Find map location
621  newMapIndex = zIt + zDim * ( yIt + yDim * xIt );
622 
623  //==================================== If less than needed, add zero
624  if ( xIt < addXPre ) { newMap[newMapIndex] = 0.0; continue; }
625  if ( yIt < addYPre ) { newMap[newMapIndex] = 0.0; continue; }
626  if ( zIt < addZPre ) { newMap[newMapIndex] = 0.0; continue; }
627 
628  //==================================== If more than needed, add zero
629  if ( xIt >= ( addXPre + xDimIndices ) ) { newMap[newMapIndex] = 0.0; continue; }
630  if ( yIt >= ( addYPre + yDimIndices ) ) { newMap[newMapIndex] = 0.0; continue; }
631  if ( zIt >= ( addZPre + zDimIndices ) ) { newMap[newMapIndex] = 0.0; continue; }
632 
633  //==================================== If not padding, copy original values
634  oldMapIndex = (zIt-addZPre) + zDimIndices * ( (yIt-addYPre) + yDimIndices * (xIt-addXPre) );
635  newMap[newMapIndex] = oldMap[oldMapIndex];
636  }
637  }
638  }
639 
640  //======================================== Done
641  return ;
642 
643 }
644 
657 void ProSHADE_internal_data::ProSHADE_data::rotateMapReciprocalSpace ( ProSHADE_settings* settings, proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma )
658 {
659  //================================================ Set maximum comparison bandwidth to maximum object bandwidth
660  this->maxCompBand = this->spheres[this->noSpheres-1]->getLocalBandwidth();
661 
662  //================================================ Save map COM after processing but before rotation
663  this->findMapCOM ( );
664  this->mapCOMProcessChangeX += ( this->xCom - this->originalMapXCom );
665  this->mapCOMProcessChangeY += ( this->yCom - this->originalMapYCom );
666  this->mapCOMProcessChangeZ += ( this->zCom - this->originalMapZCom );
667 
668  //================================================ Compute the Wigner D matrices for the Euler angles
669  ProSHADE_internal_wigner::computeWignerMatricesForRotation ( settings, this, -eulerAlpha, eulerBeta, -eulerGamma );
670 
671  //================================================ Initialise rotated Spherical Harmonics memory
672  this->allocateRotatedSHMemory ( );
673 
674  //================================================ Multiply SH coeffs by Wigner
675  this->computeRotatedSH ( );
676 
677  //================================================ Inverse the SH coeffs to shells
678  this->invertSHCoefficients ( );
679 
680  //================================================ Find spherical cut-offs
681  std::vector<proshade_double> lonCO, latCO;
682  ProSHADE_internal_overlay::computeAngularThreshold ( &lonCO, &latCO, settings->maxBandwidth * 2 );
683 
684  //================================================ Allocate memory for the rotated map
685  proshade_double *densityMapRotated = new proshade_double [this->xDimIndices * this->yDimIndices * this->zDimIndices];
686  ProSHADE_internal_misc::checkMemoryAllocation ( densityMapRotated, __FILE__, __LINE__, __func__ );
687  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { densityMapRotated[iter] = 0.0; }
688 
689  //================================================ Interpolate onto cartesian grid
690  this->interpolateMapFromSpheres ( densityMapRotated );
691 
692  //================================================ Copy map
693  for ( proshade_unsign iter = 0; iter < ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
694  {
695  this->internalMap[iter] = densityMapRotated[iter];
696  }
697 
698  //================================================ Release rotated map (original is now rotated)
699  delete[] densityMapRotated;
700 
701  //================================================ Done
702  return ;
703 
704 }
705 
718 std::vector< proshade_double > ProSHADE_internal_data::ProSHADE_data::rotateMapRealSpace ( proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axAng, proshade_double*& map )
719 {
720  //================================================ Initialise local variables
721  bool withinBounds = true;
722  proshade_double c000, c001, c010, c011, c100, c101, c110, c111, c00, c01, c10, c11, c0, c1;
723  size_t arrPos = 0;
724  proshade_double xCOM, yCOM, zCOM;
725  std::vector< proshade_double > ret;
726 
727  //================================================ Store sampling rates
728  proshade_single xSampRate = this->xDimSize / static_cast< proshade_single > ( this->xTo - this->xFrom );
729  proshade_single ySampRate = this->yDimSize / static_cast< proshade_single > ( this->yTo - this->yFrom );
730  proshade_single zSampRate = this->zDimSize / static_cast< proshade_single > ( this->zTo - this->zFrom );
731 
732  //================================================ Compute map COM
733  ProSHADE_internal_mapManip::findMAPCOMValues ( this->internalMap, &xCOM, &yCOM, &zCOM, this->xDimSize, this->yDimSize, this->zDimSize, this->xFrom, this->xTo, this->yFrom, this->yTo, this->zFrom, this->zTo );
734 
735  //================================================ Allocate local variables
736  proshade_single *mins = new proshade_single[3];
737  proshade_single *maxs = new proshade_single[3];
738  proshade_single *rotMat = new proshade_single[9];
739  proshade_single *rotVec;
740  proshade_single *interpMins = new proshade_single[3];
741  proshade_single *interpMaxs = new proshade_single[3];
742  proshade_single *interpDiff = new proshade_single[3];
743  proshade_single *movs = new proshade_single[3];
744  map = new proshade_double[ this->xDimIndices * this->yDimIndices * this->zDimIndices ];
745 
746  //================================================ Check memory allocation
747  ProSHADE_internal_misc::checkMemoryAllocation ( mins, __FILE__, __LINE__, __func__ );
748  ProSHADE_internal_misc::checkMemoryAllocation ( maxs, __FILE__, __LINE__, __func__ );
749  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
750  ProSHADE_internal_misc::checkMemoryAllocation ( interpMins, __FILE__, __LINE__, __func__ );
751  ProSHADE_internal_misc::checkMemoryAllocation ( interpMaxs, __FILE__, __LINE__, __func__ );
752  ProSHADE_internal_misc::checkMemoryAllocation ( interpDiff, __FILE__, __LINE__, __func__ );
753  ProSHADE_internal_misc::checkMemoryAllocation ( movs, __FILE__, __LINE__, __func__ );
754  ProSHADE_internal_misc::checkMemoryAllocation ( map, __FILE__, __LINE__, __func__ );
755 
756  //================================================ Fill map with zeroes
757  for ( size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ ) { map[iter] = 0.0; }
758 
759  //================================================ Determine map max's and min's in terms of the hkl system
760  mins[0] = std::floor ( static_cast< proshade_single > ( this->xDimIndices ) / -2.0f );
761  mins[1] = std::floor ( static_cast< proshade_single > ( this->yDimIndices ) / -2.0f );
762  mins[2] = std::floor ( static_cast< proshade_single > ( this->zDimIndices ) / -2.0f );
763 
764  maxs[0] = -mins[0];
765  maxs[1] = -mins[1];
766  maxs[2] = -mins[2];
767 
768  if ( this->xDimIndices % 2 == 0 ) { maxs[0] -= 1.0f; }
769  if ( this->yDimIndices % 2 == 0 ) { maxs[1] -= 1.0f; }
770  if ( this->zDimIndices % 2 == 0 ) { maxs[2] -= 1.0f; }
771 
772  //================================================ Rotate about COM instead of map midpoint
773  movs[0] = 0.0; //( mins[0] - static_cast< proshade_single > ( this->xFrom ) );
774  movs[1] = 0.0; //( mins[1] - static_cast< proshade_single > ( this->yFrom ) );
775  movs[2] = 0.0; //( mins[2] - static_cast< proshade_single > ( this->zFrom ) );
776 
777  //================================================ Save rotation centre
778  ProSHADE_internal_misc::addToDoubleVector ( &ret, static_cast< proshade_double > ( ( -mins[0] * xSampRate ) + ( static_cast< proshade_single > ( this->xFrom ) * xSampRate ) ) );
779  ProSHADE_internal_misc::addToDoubleVector ( &ret, static_cast< proshade_double > ( ( -mins[1] * ySampRate ) + ( static_cast< proshade_single > ( this->yFrom ) * ySampRate ) ) );
780  ProSHADE_internal_misc::addToDoubleVector ( &ret, static_cast< proshade_double > ( ( -mins[2] * zSampRate ) + ( static_cast< proshade_single > ( this->zFrom ) * zSampRate ) ) );
781 
782  //================================================ Get rotation matrix from Euler angles
783  ProSHADE_internal_maths::getRotationMatrixFromAngleAxis ( rotMat, axX, axY, axZ, axAng );
784 
785  //================================================ For each point
786  for ( proshade_single xIt = mins[0]; xIt <= maxs[0]; xIt += 1.0f )
787  {
788  for ( proshade_single yIt = mins[1]; yIt <= maxs[1]; yIt += 1.0f )
789  {
790  for ( proshade_single zIt = mins[2]; zIt <= maxs[2]; zIt += 1.0f )
791  {
792  //==================================== Compute new point position
793  rotVec = ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication ( rotMat, xIt - movs[0], yIt - movs[1], zIt - movs[2] );
794 
795  //==================================== Find surrounding grid points indices and check for boundaries
796  withinBounds = true;
797  for ( size_t posIt = 0; posIt < 3; posIt++ )
798  {
799  //================================ Determine surrounding points indices in this dimension
800  interpMins[posIt] = std::floor ( rotVec[posIt] );
801  interpMaxs[posIt] = interpMins[posIt] + 1.0f;
802 
803  //================================ Check for boundaries
804  if ( ( maxs[posIt] < interpMins[posIt] ) || ( interpMins[posIt] < mins[posIt] ) || ( maxs[posIt] < interpMaxs[posIt] ) || ( interpMaxs[posIt] < mins[posIt] ) )
805  {
806  withinBounds = false;
807  break;
808  }
809 
810  //================================ Compute the difference from position to min index along this axis
811  interpDiff[posIt] = rotVec[posIt] - interpMins[posIt];
812  }
813  if ( !withinBounds ) { continue; }
814 
815  //==================================== Make sure interpolation max's are within bounds
816  for ( size_t posIt = 0; posIt < 3; posIt++ )
817  {
818  interpMaxs[posIt] = std::min ( maxs[posIt], std::max ( mins[posIt], interpMaxs[posIt] ) );
819  }
820 
821  //==================================== Release memory
822  delete[] rotVec;
823 
824  //==================================== Find the surrounding points values from their indices
825  arrPos = static_cast< size_t > ( ( interpMins[2] - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( interpMins[1] - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( interpMins[0] - mins[0] ) ) );
826  c000 = this->internalMap[arrPos];
827 
828  arrPos = static_cast< size_t > ( ( interpMaxs[2] - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( interpMins[1] - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( interpMins[0] - mins[0] ) ) );
829  c001 = this->internalMap[arrPos];
830 
831  arrPos = static_cast< size_t > ( ( interpMins[2] - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( interpMaxs[1] - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( interpMins[0] - mins[0] ) ) );
832  c010 = this->internalMap[arrPos];
833 
834  arrPos = static_cast< size_t > ( ( interpMaxs[2] - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( interpMaxs[1] - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( interpMins[0] - mins[0] ) ) );
835  c011 = this->internalMap[arrPos];
836 
837  arrPos = static_cast< size_t > ( ( interpMins[2] - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( interpMins[1] - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( interpMaxs[0] - mins[0] ) ) );
838  c100 = this->internalMap[arrPos];
839 
840  arrPos = static_cast< size_t > ( ( interpMaxs[2] - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( interpMins[1] - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( interpMaxs[0] - mins[0] ) ) );
841  c101 = this->internalMap[arrPos];
842 
843  arrPos = static_cast< size_t > ( ( interpMins[2] - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( interpMaxs[1] - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( interpMaxs[0] - mins[0] ) ) );
844  c110 = this->internalMap[arrPos];
845 
846  arrPos = static_cast< size_t > ( ( interpMaxs[2] - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( interpMaxs[1] - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( interpMaxs[0] - mins[0] ) ) );
847  c111 = this->internalMap[arrPos];
848 
849  //==================================== Interpolate along x-axis
850  c00 = ( c000 * ( 1.0 - static_cast< proshade_double > ( interpDiff[0] ) ) ) + ( c100 * static_cast< proshade_double > ( interpDiff[0] ) );
851  c01 = ( c001 * ( 1.0 - static_cast< proshade_double > ( interpDiff[0] ) ) ) + ( c101 * static_cast< proshade_double > ( interpDiff[0] ) );
852  c10 = ( c010 * ( 1.0 - static_cast< proshade_double > ( interpDiff[0] ) ) ) + ( c110 * static_cast< proshade_double > ( interpDiff[0] ) );
853  c11 = ( c011 * ( 1.0 - static_cast< proshade_double > ( interpDiff[0] ) ) ) + ( c111 * static_cast< proshade_double > ( interpDiff[0] ) );
854 
855  //==================================== Interpolate along y-axis
856  c0 = ( c00 * ( 1.0 - static_cast< proshade_double > ( interpDiff[1] ) ) ) + ( c10 * static_cast< proshade_double > ( interpDiff[1] ) );
857  c1 = ( c01 * ( 1.0 - static_cast< proshade_double > ( interpDiff[1] ) ) ) + ( c11 * static_cast< proshade_double > ( interpDiff[1] ) );
858 
859  //==================================== Interpolate along z-axis
860  arrPos = static_cast< size_t > ( ( zIt - mins[2] ) + static_cast< proshade_single > ( this->zDimIndices ) * ( ( yIt - mins[1] ) + static_cast< proshade_single > ( this->yDimIndices ) * ( xIt - mins[0] ) ) );
861  map[arrPos] = ( c0 * ( 1.0 - static_cast< proshade_double > ( interpDiff[2] ) ) ) + ( c1 * static_cast< proshade_double > ( interpDiff[2] ) );
862  }
863  }
864  }
865 
866  //================================================ Release memory
867  delete[] mins;
868  delete[] maxs;
869  delete[] rotMat;
870  delete[] interpMins;
871  delete[] interpMaxs;
872  delete[] interpDiff;
873  delete[] movs;
874 
875  //================================================ Done
876  return ( ret );
877 
878 }
879 
891 std::vector< proshade_double > ProSHADE_internal_data::ProSHADE_data::rotateMapRealSpaceInPlace ( proshade_double eulA, proshade_double eulB, proshade_double eulG )
892 {
893  //================================================ Initialise local variables
894  proshade_double axX, axY, axZ, axAng, tmp, *rMat, *map;
895 
896  //================================================ Allocate local memory
897  rMat = new proshade_double[9];
898 
899  //================================================ Check local memory allocation
900  ProSHADE_internal_misc::checkMemoryAllocation ( rMat, __FILE__, __LINE__, __func__ );
901 
902  //================================================ Convert Euler angles to rotation matrix
904 
905  //================================================ Transpose the rotation matrix
906  tmp = rMat[1];
907  rMat[1] = rMat[3];
908  rMat[3] = tmp;
909 
910  tmp = rMat[2];
911  rMat[2] = rMat[6];
912  rMat[6] = tmp;
913 
914  tmp = rMat[5];
915  rMat[5] = rMat[7];
916  rMat[7] = tmp;
917 
918  //================================================ Convert rotation matrix to angle-axis
919  ProSHADE_internal_maths::getAxisAngleFromRotationMatrix ( rMat, &axX, &axY, &axZ, &axAng );
920 
921  //================================================ Rotate the internal map
922  std::vector< proshade_double > ret = this->rotateMapRealSpace ( axX, axY, axZ, axAng, map );
923 
924  //================================================ Copy the rotated map in place of the internal map
925  for ( size_t iter = 0; iter < static_cast< size_t > ( this->xDimIndices * this->yDimIndices * this->zDimIndices ); iter++ )
926  {
927  this->internalMap[iter] = map[iter];
928  }
929 
930  //================================================ Release local memory
931  delete[] rMat;
932  delete[] map;
933 
934  //================================================ Save rotation centre for co-ordinates
935  this->originalPdbRotCenX = ret.at(0);
936  this->originalPdbRotCenY = ret.at(1);
937  this->originalPdbRotCenZ = ret.at(2);
938 
939  //================================================ Done
940  return ( ret );
941 
942 }
943 
954 void ProSHADE_internal_data::ProSHADE_data::translateMap ( proshade_double trsX, proshade_double trsY, proshade_double trsZ )
955 {
956  //================================================ Initialise local variables
957  proshade_single xMov = static_cast< proshade_single > ( -trsX );
958  proshade_single yMov = static_cast< proshade_single > ( -trsY );
959  proshade_single zMov = static_cast< proshade_single > ( -trsZ );
960 
961  //================================================ Move the whole map frame to minimise the Fourier movement
962  ProSHADE_internal_mapManip::moveMapByIndices ( &xMov, &yMov, &zMov, this->getXDimSize(), this->getYDimSize(), this->getZDimSize(),
963  this->getXFromPtr(), this->getXToPtr(), this->getYFromPtr(), this->getYToPtr(),
964  this->getZFromPtr(), this->getZToPtr(), this->getXAxisOrigin(), this->getYAxisOrigin(), this->getZAxisOrigin() );
965 
966  //================================================ Finalise the movement by in-frame Fourier movement
967  ProSHADE_internal_mapManip::moveMapByFourier ( this->getInternalMap(), xMov, yMov, zMov, this->getXDimSize(), this->getYDimSize(), this->getZDimSize(),
968  static_cast< proshade_signed > ( this->getXDim() ), static_cast< proshade_signed > ( this->getYDim() ),
969  static_cast< proshade_signed > ( this->getZDim() ) );
970 
971  //================================================ Done
972  return ;
973 
974 }
975 
979 {
980  //================================================ Allocate the main pointer and check
981  this->rotSphericalHarmonics = new proshade_complex* [this->noSpheres];
982  ProSHADE_internal_misc::checkMemoryAllocation ( this->rotSphericalHarmonics, __FILE__, __LINE__, __func__ );
983 
984  //================================================ For each sphere
985  for ( proshade_unsign iter = 0; iter < this->noSpheres; iter++ )
986  {
987  //============================================ Allocate the sphere storage space
988  this->rotSphericalHarmonics[iter] = new proshade_complex [static_cast<proshade_unsign> ( pow ( (this->spheres[iter]->getLocalBandwidth() * 2), 2) )];
989  ProSHADE_internal_misc::checkMemoryAllocation ( this->rotSphericalHarmonics[iter], __FILE__, __LINE__, __func__ );
990 
991  //============================================ Set values to zeroes
992  for ( proshade_unsign it = 0; it < static_cast<proshade_unsign> ( pow ( (this->spheres[iter]->getLocalBandwidth() * 2), 2) ); it++ )
993  {
994  this->rotSphericalHarmonics[iter][it][0] = 0.0;
995  this->rotSphericalHarmonics[iter][it][1] = 0.0;
996  }
997  }
998 
999  //================================================ Done
1000  return ;
1001 
1002 }
1003 
1007 {
1008  //================================================ Initialise variables
1009  proshade_double WigDR, WigDI, *ShR, *ShI, retR, retI;
1010  proshade_unsign arrPos;
1011 
1012  //================================================ Compute
1013  for ( proshade_signed shell = 0; shell < static_cast<proshade_signed> ( this->noSpheres ); shell++ )
1014  {
1015  //============================================ For each band
1016  for ( proshade_signed bandIter = 0; bandIter < static_cast<proshade_signed> ( this->spheres[shell]->getLocalBandwidth() ); bandIter++ )
1017  {
1018  //======================================== For each order1
1019  for ( proshade_signed order1 = 0; order1 < ( ( bandIter * 2 ) + 1 ); order1++ )
1020  {
1021  //==================================== Get Spherical Harmonics value
1022  ShR = getRealSphHarmValue ( static_cast< proshade_unsign > ( bandIter ), static_cast< proshade_unsign > ( order1 ), static_cast< proshade_unsign > ( shell ) );
1023  ShI = getImagSphHarmValue ( static_cast< proshade_unsign > ( bandIter ), static_cast< proshade_unsign > ( order1 ), static_cast< proshade_unsign > ( shell ) );
1024 
1025  //==================================== For each order2
1026  for ( proshade_signed order2 = 0; order2 < ( ( bandIter * 2 ) + 1 ); order2++ )
1027  {
1028  //================================ Get Wigner D values
1029  this->getWignerMatrixValue ( static_cast< proshade_unsign > ( bandIter ), static_cast< proshade_unsign > ( order1 ), static_cast< proshade_unsign > ( order2 ), &WigDR, &WigDI );
1030 
1031  //================================ Multiply SH and Wigner
1032  ProSHADE_internal_maths::complexMultiplication ( ShR, ShI, &WigDR, &WigDI, &retR, &retI );
1033 
1034  //================================ Save
1035  arrPos = static_cast<proshade_unsign> ( seanindex ( static_cast< int > ( order2-bandIter ), static_cast< int > ( bandIter ),
1036  static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ) ) );
1037  this->rotSphericalHarmonics[shell][arrPos][0] += retR;
1038  this->rotSphericalHarmonics[shell][arrPos][1] += retI;
1039  }
1040  }
1041  }
1042  }
1043 
1044  //================================================ Done
1045  return ;
1046 
1047 }
1048 
1061 void ProSHADE_internal_overlay::initialiseInverseSHComputation ( proshade_unsign shBand, double*& sigR, double*& sigI, double*& rcoeffs, double*& icoeffs, double*& weights, double*& workspace, fftw_plan& idctPlan, fftw_plan& ifftPlan )
1062 {
1063  //================================================ Initialise internal variables
1064  proshade_unsign oneDim = shBand * 2;
1065 
1066  //================================================ Allocate memory
1067  sigR = new double [(oneDim * oneDim * 4)];
1068  sigI = sigR + (oneDim * oneDim * 2);
1069  rcoeffs = new double [(oneDim * oneDim * 2)];
1070  icoeffs = rcoeffs + (oneDim * oneDim);
1071  weights = new double [4 * shBand];
1072  workspace = new double [( 20 * shBand * shBand ) + ( 24 * shBand )];
1073 
1074  //================================================ Check memory allocation
1075  ProSHADE_internal_misc::checkMemoryAllocation ( sigR, __FILE__, __LINE__, __func__ );
1076  ProSHADE_internal_misc::checkMemoryAllocation ( rcoeffs, __FILE__, __LINE__, __func__ );
1077  ProSHADE_internal_misc::checkMemoryAllocation ( weights, __FILE__, __LINE__, __func__ );
1078  ProSHADE_internal_misc::checkMemoryAllocation ( workspace, __FILE__, __LINE__, __func__ );
1079 
1080  //================================================ Create the cosine/sine transform plan
1081  idctPlan = fftw_plan_r2r_1d ( static_cast< int > ( oneDim ), weights, workspace, FFTW_REDFT01, FFTW_ESTIMATE );
1082 
1083  //================================================ Set up the discrete Fourier transform
1084  int rank, howmany_rank;
1085  fftw_iodim dims[1], howmany_dims[1];
1086 
1087  rank = 1;
1088  dims[0].n = 2 * static_cast< int > ( shBand );
1089  dims[0].is = 2 * static_cast< int > ( shBand );
1090  dims[0].os = 1;
1091  howmany_rank = 1;
1092  howmany_dims[0].n = 2 * static_cast< int > ( shBand );
1093  howmany_dims[0].is = 1;
1094  howmany_dims[0].os = 2 * static_cast< int > ( shBand );
1095 
1096  //================================================ Create the discrete Fourier transform
1097  ifftPlan = fftw_plan_guru_split_dft ( rank, dims, howmany_rank, howmany_dims, sigR, sigI, rcoeffs, icoeffs, FFTW_ESTIMATE );
1098 
1099  //================================================ Done
1100  return ;
1101 
1102 }
1103 
1108 {
1109  //================================================ Initialise local variables
1110  double *sigR = nullptr, *sigI = nullptr, *rcoeffs = nullptr, *icoeffs = nullptr, *weights = nullptr, *workspace = nullptr;
1111  fftw_plan idctPlan, ifftPlan;
1112 
1113  //================================================ For each shell
1114  for ( int shell = 0; shell < static_cast<int> ( this->noSpheres ); shell++ )
1115  {
1116  //=========================================== Initialise internal variables
1117  proshade_unsign oneDim = this->spheres[shell]->getLocalBandwidth() * 2;
1118 
1119  //=========================================== Allocate memory
1120  ProSHADE_internal_overlay::initialiseInverseSHComputation ( this->spheres[shell]->getLocalBandwidth(), sigR, sigI, rcoeffs, icoeffs, weights, workspace, idctPlan, ifftPlan );
1121 
1122  //=========================================== Compute weights for the transform using the appropriate shell related band
1123  makeweights ( static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ), weights );
1124 
1125  //============================================ Allocate rotated shell mapped data memory
1126  this->spheres[shell]->allocateRotatedMap ( );
1127 
1128  //============================================ Load SH coeffs to arrays
1129  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( oneDim * oneDim ); iter++ )
1130  {
1131  rcoeffs[iter] = this->rotSphericalHarmonics[shell][iter][0];
1132  icoeffs[iter] = this->rotSphericalHarmonics[shell][iter][1];
1133  sigR[iter] = 0.0;
1134  sigI[iter] = 0.0;
1135  }
1136 
1137  //============================================ Get inverse spherical harmonics transform for the shell
1138  InvFST_semi_fly ( rcoeffs,
1139  icoeffs,
1140  sigR,
1141  sigI,
1142  static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ),
1143  workspace,
1144  0,
1145  static_cast< int > ( this->spheres[shell]->getLocalBandwidth() ),
1146  &idctPlan,
1147  &ifftPlan );
1148 
1149  //=========================================== Copy the results to the rotated shells array
1150  for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( oneDim * oneDim ); iter++ )
1151  {
1152  this->spheres[shell]->setRotatedMappedData ( iter, sigR[iter] );
1153  }
1154 
1155  //=========================================== Release the plans
1156  fftw_destroy_plan ( idctPlan );
1157  fftw_destroy_plan ( ifftPlan );
1158 
1159  //=========================================== Release the memory
1160  delete[] sigR;
1161  delete[] rcoeffs;
1162  delete[] weights;
1163  delete[] workspace;
1164  }
1165 
1166  //================================================ Done
1167  return ;
1168 
1169 }
1170 
1177 void ProSHADE_internal_overlay::computeAngularThreshold ( std::vector<proshade_double>* lonCO, std::vector<proshade_double>* latCO, proshade_unsign angRes )
1178 {
1179  //================================================ Longitude angular thresholds
1180  for ( proshade_unsign iter = 0; iter <= angRes; iter++ )
1181  {
1182  ProSHADE_internal_misc::addToDoubleVector ( lonCO, static_cast<proshade_double> ( iter ) * ( ( static_cast<proshade_double> ( M_PI ) * 2.0 ) / static_cast<proshade_double> ( angRes ) ) - ( static_cast<double> ( M_PI ) ) );
1183  }
1184 
1185  //================================================ Lattitude angular thresholds
1186  for ( proshade_unsign iter = 0; iter <= angRes; iter++ )
1187  {
1188  ProSHADE_internal_misc::addToDoubleVector ( latCO, static_cast<proshade_double> ( iter ) * ( static_cast<proshade_double> ( M_PI ) / static_cast<proshade_double> ( angRes ) ) - ( static_cast<proshade_double> ( M_PI ) / 2.0 ) );
1189  }
1190 
1191  //================================================ Done
1192  return ;
1193 
1194 }
1195 
1200 void ProSHADE_internal_data::ProSHADE_data::interpolateMapFromSpheres ( proshade_double*& densityMapRotated )
1201 {
1202  //================================================ Initialise variables
1203  proshade_double rad = 0.0, lon = 0.0, lat = 0.0, newU = 0.0, newV = 0.0, newW = 0.0;
1204  proshade_unsign lowerLonL = 0, upperLonL = 0, lowerLonU = 0, upperLonU = 0, lowerLatL = 0, upperLatL = 0, lowerLatU = 0, upperLatU = 0, lowerShell = 0, upperShell = 0;
1205  proshade_double x00 = 0.0, x01 = 0.0, x10 = 0.0, x11 = 0.0, distLLon = 0.0, distLLat = 0.0, distLRad = 0.0, valLLon = 0.0, valULon = 0.0;
1206  proshade_double lowerShellValue = 0.0, upperShellValue = 0.0;
1207  proshade_double xSamplingRate = static_cast<proshade_double> ( this->xDimSize ) / static_cast<proshade_double> ( this->xDimIndices );
1208  proshade_double ySamplingRate = static_cast<proshade_double> ( this->yDimSize ) / static_cast<proshade_double> ( this->yDimIndices );
1209  proshade_double zSamplingRate = static_cast<proshade_double> ( this->zDimSize ) / static_cast<proshade_double> ( this->zDimIndices );
1210  proshade_signed arrPos;
1211  std::vector<proshade_double> lonCOU, latCOU, lonCOL, latCOL;
1212 
1213  for ( proshade_signed uIt = 0; uIt < static_cast<proshade_signed> (this->xDimIndices); uIt++ )
1214  {
1215  for ( proshade_signed vIt = 0; vIt < static_cast<proshade_signed> (this->yDimIndices); vIt++ )
1216  {
1217  for ( proshade_signed wIt = 0; wIt < static_cast<proshade_signed> (this->zDimIndices); wIt++ )
1218  {
1219  //==================================== Convert to centered coords
1220  newU = static_cast<proshade_double> ( uIt - ( static_cast<proshade_signed> (this->xDimIndices) / 2 ) );
1221  newV = static_cast<proshade_double> ( vIt - ( static_cast<proshade_signed> (this->yDimIndices) / 2 ) );
1222  newW = static_cast<proshade_double> ( wIt - ( static_cast<proshade_signed> (this->zDimIndices) / 2 ) );
1223 
1224  //==================================== Deal with 0 ; 0 ; 0
1225  if ( ( newU == 0.0 ) && ( newV == 0.0 ) && ( newW == 0.0 ) )
1226  {
1227  arrPos = wIt + static_cast< proshade_signed > ( this->zDimIndices ) * ( vIt + static_cast< proshade_signed > ( this->yDimIndices ) * uIt );
1228  densityMapRotated[arrPos] = this->internalMap[arrPos];
1229  continue;
1230  }
1231 
1232  //==================================== Convert to spherical coords
1233  rad = sqrt ( pow( ( newU * xSamplingRate ), 2.0 ) +
1234  pow( ( newV * ySamplingRate ), 2.0 ) +
1235  pow( ( newW * zSamplingRate ), 2.0 ) );
1236  lon = atan2 ( ( newV * ySamplingRate ), ( newU * xSamplingRate ) );
1237  lat = asin ( ( newW * zSamplingRate ) / rad );
1238 
1239  //==================================== Deal with nan's
1240  if ( rad != rad ) { rad = 0.0; }
1241  if ( lon != lon ) { lon = 0.0; }
1242  if ( lat != lat ) { lat = 0.0; }
1243 
1244  //==================================== Find shells above and below
1245  lowerShell = 0;
1246  upperShell = 0;
1247  for ( proshade_unsign iter = 0; iter < (this->noSpheres-1); iter++ )
1248  {
1249  if ( ( static_cast< proshade_double > ( this->spherePos.at(iter) ) <= rad ) && ( static_cast< proshade_double > ( this->spherePos.at(iter+1) ) > rad ) )
1250  {
1251  lowerShell = iter;
1252  upperShell = iter+1;
1253  break;
1254  }
1255  }
1256 
1257  if ( upperShell == 0 )
1258  {
1259  arrPos = wIt + static_cast< proshade_signed > ( this->zDimIndices ) * ( vIt + static_cast< proshade_signed > ( this->yDimIndices ) * uIt );
1260  densityMapRotated[arrPos] = 0.0;
1261  continue;
1262  }
1263 
1264  //==================================== Get the longitude and lattitude cut-offs for this shell resolution
1265  lonCOL.clear(); latCOL.clear(); lonCOU.clear(); latCOU.clear();
1266  ProSHADE_internal_overlay::computeAngularThreshold ( &lonCOL, &latCOL, this->spheres[lowerShell]->getLocalAngRes() );
1267  ProSHADE_internal_overlay::computeAngularThreshold ( &lonCOU, &latCOU, this->spheres[upperShell]->getLocalAngRes() );
1268 
1269  //==================================== Find the angle cutoffs around the point
1270  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( lonCOL.size() ); iter++ )
1271  {
1272  if ( iter == ( static_cast<proshade_unsign> ( lonCOL.size() ) - 1 ) )
1273  {
1274  lowerLonL = 0;
1275  upperLonL = 1;
1276  break;
1277  }
1278  if ( ( std::floor(10000. * lonCOL.at(iter)) <= std::floor(10000. * lon) ) && ( std::floor(10000. * lonCOL.at(iter+1)) > std::floor(10000. * lon) ) )
1279  {
1280  lowerLonL = iter;
1281  upperLonL = iter+1;
1282  break;
1283  }
1284  }
1285  if ( upperLonL == this->spheres[lowerShell]->getLocalAngRes() ) { upperLonL = 0; }
1286 
1287  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( lonCOU.size() ); iter++ )
1288  {
1289  if ( iter == ( static_cast<proshade_unsign> ( lonCOU.size() ) - 1 ) )
1290  {
1291  lowerLonU = 0;
1292  upperLonU = 1;
1293  break;
1294  }
1295  if ( ( std::floor(10000. * lonCOU.at(iter)) <= std::floor(10000. * lon) ) && ( std::floor(10000. * lonCOU.at(iter+1)) > std::floor(10000. * lon) ) )
1296  {
1297  lowerLonU = iter;
1298  upperLonU = iter+1;
1299  break;
1300  }
1301  }
1302  if ( upperLonU == this->spheres[upperShell]->getLocalAngRes() ) { upperLonU = 0; }
1303 
1304  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( latCOL.size() ); iter++ )
1305  {
1306  if ( iter == ( static_cast<proshade_unsign> ( latCOL.size() ) - 1 ) )
1307  {
1308  lowerLatL = 0;
1309  upperLatL = 1;
1310  break;
1311  }
1312  if ( ( std::floor(10000. * latCOL.at(iter)) <= std::floor(10000. * lat) ) && ( std::floor(10000. * latCOL.at(iter+1)) > std::floor(10000. * lat) ) )
1313  {
1314  lowerLatL = iter;
1315  upperLatL = iter+1;
1316  break;
1317  }
1318  }
1319  if ( upperLatL == this->spheres[lowerShell]->getLocalAngRes() ) { upperLatL = 0; }
1320 
1321  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( latCOU.size() ); iter++ )
1322  {
1323  if ( iter == ( static_cast<proshade_unsign> ( latCOU.size() ) - 1 ) )
1324  {
1325  lowerLatU = 0;
1326  upperLatU = 1;
1327  break;
1328  }
1329  if ( ( std::floor(10000. * latCOU.at(iter)) <= std::floor(10000. * lat) ) && ( std::floor(10000. * latCOU.at(iter+1)) > std::floor(10000. * lat) ) )
1330  {
1331  lowerLatU = iter;
1332  upperLatU = iter+1;
1333  break;
1334  }
1335  }
1336  if ( upperLatU == this->spheres[upperShell]->getLocalAngRes() ) { upperLatU = 0; }
1337 
1338  //==================================== Interpolate lower shell
1339  x00 = this->spheres[lowerShell]->getRotatedMappedData ( lowerLatL * this->spheres[lowerShell]->getLocalAngRes() + lowerLonL );
1340  x01 = this->spheres[lowerShell]->getRotatedMappedData ( lowerLatL * this->spheres[lowerShell]->getLocalAngRes() + upperLonL );
1341  x10 = this->spheres[lowerShell]->getRotatedMappedData ( upperLatL * this->spheres[lowerShell]->getLocalAngRes() + lowerLonL );
1342  x11 = this->spheres[lowerShell]->getRotatedMappedData ( upperLatL * this->spheres[lowerShell]->getLocalAngRes() + upperLonL );
1343 
1344  distLLon = std::abs ( lon - lonCOL.at(lowerLonL) ) / ( std::abs( lon - lonCOL.at(lowerLonL) ) + std::abs( lon - lonCOL.at(upperLonL) ) );
1345  valLLon = ( ( 1.0 - distLLon ) * x00 ) + ( distLLon * x01 );
1346  valULon = ( ( 1.0 - distLLon ) * x10 ) + ( distLLon * x11 );
1347 
1348  distLLat = std::abs ( lat - latCOL.at(lowerLatL) ) / ( std::abs( lat - latCOL.at(lowerLatL) ) + std::abs( lat - latCOL.at(upperLatL) ) );
1349  lowerShellValue = ( ( 1.0 - distLLat ) * valLLon ) + ( distLLat * valULon );
1350 
1351  //==================================== Interpolate upper shell
1352  x00 = this->spheres[upperShell]->getRotatedMappedData ( lowerLatU * this->spheres[upperShell]->getLocalAngRes() + lowerLonU );
1353  x01 = this->spheres[upperShell]->getRotatedMappedData ( lowerLatU * this->spheres[upperShell]->getLocalAngRes() + upperLonU );
1354  x10 = this->spheres[upperShell]->getRotatedMappedData ( upperLatU * this->spheres[upperShell]->getLocalAngRes() + lowerLonU );
1355  x11 = this->spheres[upperShell]->getRotatedMappedData ( upperLatU * this->spheres[upperShell]->getLocalAngRes() + upperLonU );
1356 
1357  distLLon = std::abs ( lon - lonCOU.at(lowerLonU) ) / ( std::abs( lon - lonCOU.at(lowerLonU) ) + std::abs( lon - lonCOU.at(upperLonU) ) );
1358  valLLon = ( ( 1.0 - distLLon ) * x00 ) + ( distLLon * x01 );
1359  valULon = ( ( 1.0 - distLLon ) * x10 ) + ( distLLon * x11 );
1360 
1361  distLLat = std::abs ( lat - latCOU.at(lowerLatU) ) / ( std::abs( lat - latCOU.at(lowerLatU) ) + std::abs( lat - latCOU.at(upperLatU) ) );
1362  upperShellValue = ( ( 1.0 - distLLat ) * valLLon ) + ( distLLat * valULon );
1363 
1364  //==================================== Interpolate between shells
1365  distLRad = std::abs ( rad - static_cast< proshade_double > ( this->spherePos.at(lowerShell) ) ) / ( std::abs( rad - static_cast< proshade_double > ( this->spherePos.at(lowerShell) ) ) +
1366  std::abs( rad - static_cast< proshade_double > ( this->spherePos.at(upperShell) ) ) );
1367 
1368  arrPos = wIt + static_cast< proshade_signed > ( this->zDimIndices ) * ( vIt + static_cast< proshade_signed > ( this->yDimIndices ) * uIt );
1369  densityMapRotated[arrPos] = ( ( 1.0 - distLRad ) * lowerShellValue ) + ( distLRad * upperShellValue );
1370  }
1371 
1372  }
1373 
1374  }
1375 
1376  //================================================ Done
1377  return ;
1378 
1379 }
1380 
1387 {
1388  //================================================ Initialise local variables
1389  std::vector < proshade_double > ret;
1390  proshade_double eulA, eulB, eulG;
1391 
1392  //================================================ Get inverse SO(3) map top peak Euler angle values
1394  this->getMaxBand() * 2,
1395  &eulA, &eulB, &eulG, settings );
1396 
1397  //================================================ Re-format to vector
1401 
1402  //================================================ Done
1403  return ( ret );
1404 
1405 }
ProSHADE_internal_distances::normaliseEMatrices
void normaliseEMatrices(ProSHADE_internal_data::ProSHADE_data *obj1, ProSHADE_internal_data::ProSHADE_data *obj2, ProSHADE_settings *settings)
This function normalises the E matrices.
Definition: ProSHADE_distances.cpp:572
ProSHADE_settings::maxBandwidth
proshade_unsign maxBandwidth
The bandwidth of spherical harmonics decomposition for the largest sphere.
Definition: ProSHADE_settings.hpp:58
ProSHADE_internal_data::ProSHADE_data::computeSphericalHarmonics
void computeSphericalHarmonics(ProSHADE_settings *settings)
This function computes the spherical harmonics decomposition for the whole structure.
Definition: ProSHADE_data.cpp:1843
ProSHADE_internal_data::ProSHADE_data::mapMovFromsChangeX
proshade_double mapMovFromsChangeX
When the map is translated, the xFrom and xTo values are changed. This variable holds how much they h...
Definition: ProSHADE_data.hpp:94
ProSHADE_internal_data::ProSHADE_data::allocateRotatedSHMemory
void allocateRotatedSHMemory(void)
This function allocates the memory required for storing the rotated Spherical Harmonics coefficients.
Definition: ProSHADE_overlay.cpp:978
ProSHADE_internal_data::ProSHADE_data::getOverlayRotationFunction
void getOverlayRotationFunction(ProSHADE_settings *settings, ProSHADE_internal_data::ProSHADE_data *obj2)
This function computes the overlay rotation function (i.e. the correlation function in SO(3) space).
Definition: ProSHADE_overlay.cpp:35
ProSHADE_internal_data::ProSHADE_data::zFrom
proshade_signed zFrom
This is the starting index along the z axis.
Definition: ProSHADE_data.hpp:112
ProSHADE_internal_data::ProSHADE_data::xDimIndices
proshade_unsign xDimIndices
This is the size of the map cell x dimension in indices.
Definition: ProSHADE_data.hpp:65
ProSHADE_internal_data::ProSHADE_data::getZDimSize
proshade_single getZDimSize(void)
This function allows access to the map size in angstroms along the Z axis.
Definition: ProSHADE_data.cpp:3898
ProSHADE_internal_data::ProSHADE_data::getBestTranslationMapPeaksAngstrom
std::vector< proshade_double > getBestTranslationMapPeaksAngstrom(ProSHADE_internal_data::ProSHADE_data *staticStructure)
This function gets the optimal translation vector and returns it as a standard library vector....
Definition: ProSHADE_overlay.cpp:250
ProSHADE_settings::determineAllSHValues
void determineAllSHValues(proshade_unsign xDim, proshade_unsign yDim, proshade_single xDimAngs, proshade_single yDimAngs, proshade_single zDimAngs)
This function determines all the required values for spherical harmonics computation.
Definition: ProSHADE.cpp:1615
ProSHADE_internal_data::ProSHADE_data::interpolateMapFromSpheres
void interpolateMapFromSpheres(proshade_double *&densityMapRotated)
This function interpolates the density map from the sphere mapped data.
Definition: ProSHADE_overlay.cpp:1200
ProSHADE_internal_data::ProSHADE_data::getXAxisOrigin
proshade_signed * getXAxisOrigin(void)
This function allows access to the map X axis origin value.
Definition: ProSHADE_data.cpp:3998
ProSHADE_exception
This class is the representation of ProSHADE exception.
Definition: ProSHADE_exceptions.hpp:37
ProSHADE_internal_data::ProSHADE_data::getXDimSize
proshade_single getXDimSize(void)
This function allows access to the map size in angstroms along the X axis.
Definition: ProSHADE_data.cpp:3878
ProSHADE_internal_data::ProSHADE_data::getInternalMap
proshade_double *& getInternalMap(void)
This function allows access to the first map array value address.
Definition: ProSHADE_data.cpp:4028
ProSHADE_internal_data::ProSHADE_data
This class contains all inputed and derived data for a single structure.
Definition: ProSHADE_data.hpp:49
ProSHADE_internal_data::ProSHADE_data::getYDimSize
proshade_single getYDimSize(void)
This function allows access to the map size in angstroms along the Y axis.
Definition: ProSHADE_data.cpp:3888
ProSHADE_internal_data::ProSHADE_data::zDimSize
proshade_single zDimSize
This is the size of the map cell z dimension in Angstroms.
Definition: ProSHADE_data.hpp:61
ProSHADE_internal_mapManip::findMAPCOMValues
void findMAPCOMValues(proshade_double *map, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo)
This function finds the Centre of Mass for a map.
Definition: ProSHADE_mapManip.cpp:240
ProSHADE_internal_data::ProSHADE_data::getBestRotationMapPeaksEulerAngles
std::vector< proshade_double > getBestRotationMapPeaksEulerAngles(ProSHADE_settings *settings)
This function returns a vector of three floats, the three Euler angles of the best peak in the rotati...
Definition: ProSHADE_overlay.cpp:1386
ProSHADE_internal_mapManip::moveMapByFourier
void moveMapByFourier(proshade_double *&map, proshade_single xMov, proshade_single yMov, proshade_single zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim)
Function for moving map back to original PDB location by using Fourier transformation.
Definition: ProSHADE_mapManip.cpp:811
ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles
void getRotationMatrixFromEulerZXZAngles(proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double *matrix)
Function to find the rotation matrix from Euler angles (ZXZ convention).
Definition: ProSHADE_maths.cpp:1007
ProSHADE_internal_overlay::computeBeforeAfterZeroCounts
void computeBeforeAfterZeroCounts(proshade_unsign *addXPre, proshade_unsign *addYPre, proshade_unsign *addZPre, proshade_unsign *addXPost, proshade_unsign *addYPost, proshade_unsign *addZPost, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_unsign xDimIndices, proshade_unsign yDimIndices, proshade_unsign zDimIndices)
This function finds the number of zeroes to be added after and before the structure along each dimens...
Definition: ProSHADE_overlay.cpp:578
ProSHADE_internal_maths::complexMultiplication
void complexMultiplication(proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2, proshade_double *retReal, proshade_double *retImag)
Function to multiply two complex numbers.
Definition: ProSHADE_maths.cpp:38
ProSHADE_internal_overlay::getOptimalRotation
void getOptimalRotation(ProSHADE_settings *settings, ProSHADE_internal_data::ProSHADE_data *staticStructure, ProSHADE_internal_data::ProSHADE_data *movingStructure, proshade_double *eulA, proshade_double *eulB, proshade_double *eulG)
This function finds the optimal rotation between two structures as described by the settings object.
Definition: ProSHADE_overlay.cpp:74
ProSHADE_internal_overlay::freeTranslationFunctionMemory
void freeTranslationFunctionMemory(fftw_complex *&tmpIn1, fftw_complex *&tmpOut1, fftw_complex *&tmpIn2, fftw_complex *&tmpOut2, fftw_complex *&resOut, fftw_plan &forwardFourierObj1, fftw_plan &forwardFourierObj2, fftw_plan &inverseFourierCombo)
This function releases the memory for the Fourier transforms required for translation function comput...
Definition: ProSHADE_overlay.cpp:406
ProSHADE_internal_data::ProSHADE_data::getXFromPtr
proshade_signed * getXFromPtr(void)
This function allows access to the map start along the X axis.
Definition: ProSHADE_data.cpp:3938
ProSHADE_settings::verbose
proshade_signed verbose
Should the software report on the progress, or just be quiet? Value between -1 (nothing) and 4 (loud)
Definition: ProSHADE_settings.hpp:142
ProSHADE_internal_data::ProSHADE_data::yDimIndices
proshade_unsign yDimIndices
This is the size of the map cell y dimension in indices.
Definition: ProSHADE_data.hpp:66
ProSHADE_internal_misc::addToDoubleVector
void addToDoubleVector(std::vector< proshade_double > *vecToAddTo, proshade_double elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:77
ProSHADE_internal_data::ProSHADE_data::getXDim
proshade_unsign getXDim(void)
This function allows access to the map size in indices along the X axis.
Definition: ProSHADE_data.cpp:3908
ProSHADE_internal_data::ProSHADE_data::originalPdbTransX
proshade_double originalPdbTransX
The optimal translation vector as it relates to the original PDB positions (and not the ProSHADE inte...
Definition: ProSHADE_data.hpp:105
ProSHADE_internal_data::ProSHADE_data::rotateMapReciprocalSpace
void rotateMapReciprocalSpace(ProSHADE_settings *settings, proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma)
This function rotates a map based on the given Euler angles.
Definition: ProSHADE_overlay.cpp:657
ProSHADE_internal_data::ProSHADE_data::readInStructure
void readInStructure(std::string fName, proshade_unsign inputO, ProSHADE_settings *settings, proshade_double *maskArr=nullptr, proshade_unsign maskXDim=0, proshade_unsign maskYDim=0, proshade_unsign maskZDim=0, proshade_double *weightsArr=nullptr, proshade_unsign weigXDim=0, proshade_unsign weigYDim=0, proshade_unsign weigZDim=0)
This function initialises the basic ProSHADE_data variables and reads in a single structure.
Definition: ProSHADE_data.cpp:509
ProSHADE_internal_data::ProSHADE_data::mapToSpheres
void mapToSpheres(ProSHADE_settings *settings)
This function converts the internal map onto a set of concentric spheres.
Definition: ProSHADE_data.cpp:1799
ProSHADE_internal_data::ProSHADE_data::originalPdbTransY
proshade_double originalPdbTransY
The optimal translation vector as it relates to the original PDB positions (and not the ProSHADE inte...
Definition: ProSHADE_data.hpp:106
ProSHADE_internal_wigner::computeWignerMatricesForRotation
void computeWignerMatricesForRotation(ProSHADE_settings *settings, ProSHADE_internal_data::ProSHADE_data *obj, proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma)
This function computes the Wigner D matrices for a particular set of Euler angles.
Definition: ProSHADE_wignerMatrices.cpp:258
ProSHADE_overlay.hpp
This header file declares the functions required for the structure overlay computation.
ProSHADE_internal_data::ProSHADE_data::xDimSize
proshade_single xDimSize
This is the size of the map cell x dimension in Angstroms.
Definition: ProSHADE_data.hpp:59
ProSHADE_internal_data::ProSHADE_data::computeRotatedSH
void computeRotatedSH(void)
This function multiplies the objects spherical harmonics with the Wigner D matrices,...
Definition: ProSHADE_overlay.cpp:1006
ProSHADE_internal_data::ProSHADE_data::mapMovFromsChangeY
proshade_double mapMovFromsChangeY
When the map is translated, the yFrom and yTo values are changed. This variable holds how much they h...
Definition: ProSHADE_data.hpp:95
ProSHADE_internal_peakSearch::getBestPeakEulerAngsNaive
void getBestPeakEulerAngsNaive(proshade_complex *map, proshade_unsign dim, proshade_double *eulA, proshade_double *eulB, proshade_double *eulG, ProSHADE_settings *settings)
This function finds the highest peaks optimised Euler angles using the "naive" approach.
Definition: ProSHADE_peakSearch.cpp:352
ProSHADE_internal_data::ProSHADE_data::computeTranslationMap
void computeTranslationMap(ProSHADE_internal_data::ProSHADE_data *obj1)
This function does the computation of the translation map and saves results internally.
Definition: ProSHADE_overlay.cpp:324
ProSHADE_internal_data::ProSHADE_data::xFrom
proshade_signed xFrom
This is the starting index along the x axis.
Definition: ProSHADE_data.hpp:110
ProSHADE_internal_maths::getRotationMatrixFromAngleAxis
void getRotationMatrixFromAngleAxis(proshade_double *rotMat, proshade_double x, proshade_double y, proshade_double z, proshade_double ang)
This function converts the axis-angle representation to the rotation matrix representation.
Definition: ProSHADE_maths.cpp:1444
ProSHADE_internal_data::ProSHADE_data::getMaxBand
proshade_unsign getMaxBand(void)
This function returns the maximum band value for the object.
Definition: ProSHADE_data.cpp:3599
ProSHADE_internal_data::ProSHADE_data::translateMap
void translateMap(proshade_double trsX, proshade_double trsY, proshade_double trsZ)
This function simply translates the map by a given number of Angstroms along the three axes.
Definition: ProSHADE_overlay.cpp:954
ProSHADE_internal_data::ProSHADE_data::invertSHCoefficients
void invertSHCoefficients(void)
This function computes the shell mapped data from inverting the Spherical Harmonics coefficients.
Definition: ProSHADE_overlay.cpp:1107
ProSHADE_internal_overlay::getOptimalTranslation
void getOptimalTranslation(ProSHADE_settings *settings, ProSHADE_internal_data::ProSHADE_data *staticStructure, ProSHADE_internal_data::ProSHADE_data *movingStructure, proshade_double *trsX, proshade_double *trsY, proshade_double *trsZ, proshade_double eulA, proshade_double eulB, proshade_double eulG)
This function finds the optimal translation between two structures as described by the settings objec...
Definition: ProSHADE_overlay.cpp:123
ProSHADE_internal_data::ProSHADE_data::getYDim
proshade_unsign getYDim(void)
This function allows access to the map size in indices along the Y axis.
Definition: ProSHADE_data.cpp:3918
ProSHADE_settings
This class stores all the settings and is passed to the executive classes instead of a multitude of p...
Definition: ProSHADE_settings.hpp:37
ProSHADE_internal_data::ProSHADE_data::originalPdbTransZ
proshade_double originalPdbTransZ
The optimal translation vector as it relates to the original PDB positions (and not the ProSHADE inte...
Definition: ProSHADE_data.hpp:107
ProSHADE_internal_data::ProSHADE_data::processInternalMap
void processInternalMap(ProSHADE_settings *settings)
This function simply clusters several other functions which should be called together.
Definition: ProSHADE_data.cpp:1699
ProSHADE_internal_distances::computeEMatrices
void computeEMatrices(ProSHADE_internal_data::ProSHADE_data *obj1, ProSHADE_internal_data::ProSHADE_data *obj2, ProSHADE_settings *settings)
This function computes the complete E matrices and their weights between any two objects.
Definition: ProSHADE_distances.cpp:510
ProSHADE_internal_data::ProSHADE_data::getTranslationFnPointer
proshade_complex * getTranslationFnPointer(void)
This function allows access to the translation function through a pointer.
Definition: ProSHADE_data.cpp:4038
ProSHADE_internal_data::ProSHADE_data::rotateMapRealSpaceInPlace
std::vector< proshade_double > rotateMapRealSpaceInPlace(proshade_double eulA, proshade_double eulB, proshade_double eulG)
This function rotates a map based on the given Euler angles in place.
Definition: ProSHADE_overlay.cpp:891
ProSHADE_internal_data::ProSHADE_data::getZDim
proshade_unsign getZDim(void)
This function allows access to the map size in indices along the Z axis.
Definition: ProSHADE_data.cpp:3928
ProSHADE_internal_overlay::computeAngularThreshold
void computeAngularThreshold(std::vector< proshade_double > *lonCO, std::vector< proshade_double > *latCO, proshade_unsign angRes)
This function computes the angular thresholds for longitude and lattitude angles.
Definition: ProSHADE_overlay.cpp:1177
ProSHADE_internal_overlay::initialiseInverseSHComputation
void initialiseInverseSHComputation(proshade_unsign shBand, double *&sigR, double *&sigI, double *&rcoeffs, double *&icoeffs, double *&weights, double *&workspace, fftw_plan &idctPlan, fftw_plan &ifftPlan)
This function initialises internal variables for inverse Spherical Harmonics computation.
Definition: ProSHADE_overlay.cpp:1061
ProSHADE_internal_overlay::findHighestValueInMap
void findHighestValueInMap(fftw_complex *resIn, proshade_unsign xD, proshade_unsign yD, proshade_unsign zD, proshade_double *trsX, proshade_double *trsY, proshade_double *trsZ, proshade_double *mapPeak)
This function simply finds the highest value in fftw_complex map and returns its position and value.
Definition: ProSHADE_overlay.cpp:479
ProSHADE_internal_data::ProSHADE_data::mapMovFromsChangeZ
proshade_double mapMovFromsChangeZ
When the map is translated, the zFrom and zTo values are changed. This variable holds how much they h...
Definition: ProSHADE_data.hpp:96
ProSHADE_internal_distances::computeInverseSOFTTransform
void computeInverseSOFTTransform(ProSHADE_internal_data::ProSHADE_data *obj1, ProSHADE_internal_data::ProSHADE_data *obj2, ProSHADE_settings *settings)
This function computes the inverse SO(3) transform.
Definition: ProSHADE_distances.cpp:854
ProSHADE_internal_data::ProSHADE_data::getMapValue
proshade_double getMapValue(proshade_unsign pos)
This function returns the internal map representation value of a particular array position.
Definition: ProSHADE_data.cpp:3589
ProSHADE_internal_maths::getAxisAngleFromRotationMatrix
void getAxisAngleFromRotationMatrix(proshade_double *rotMat, proshade_double *x, proshade_double *y, proshade_double *z, proshade_double *ang)
This function converts rotation matrix to the axis-angle representation.
Definition: ProSHADE_maths.cpp:1124
ProSHADE_internal_data::ProSHADE_data::getYToPtr
proshade_signed * getYToPtr(void)
This function allows access to the map last position along the Y axis.
Definition: ProSHADE_data.cpp:3978
ProSHADE_internal_data::ProSHADE_data::rotateMapRealSpace
std::vector< proshade_double > rotateMapRealSpace(proshade_double axX, proshade_double axY, proshade_double axZ, proshade_double axAng, proshade_double *&map)
This function rotates a map based on the given angle-axis rotation.
Definition: ProSHADE_overlay.cpp:718
ProSHADE_internal_data::ProSHADE_data::getZAxisOrigin
proshade_signed * getZAxisOrigin(void)
This function allows access to the map Z axis origin value.
Definition: ProSHADE_data.cpp:4018
ProSHADE_internal_maths::compute3x3MatrixVectorMultiplication
proshade_double * compute3x3MatrixVectorMultiplication(proshade_double *mat, proshade_double x, proshade_double y, proshade_double z)
Function for computing a 3x3 matrix to 3x1 vector multiplication.
Definition: ProSHADE_maths.cpp:1857
ProSHADE_internal_misc::checkMemoryAllocation
void checkMemoryAllocation(chVar checkVar, std::string fileP, unsigned int lineP, std::string funcP, std::string infoP="This error may occurs when ProSHADE requests memory to be\n : allocated to it and this operation fails. This could\n : happen when not enough memory is available, either due to\n : other processes using a lot of memory, or when the machine\n : does not have sufficient memory available. Re-run to see\n : if this problem persists.")
Checks if memory was allocated properly.
Definition: ProSHADE_misc.hpp:67
ProSHADE_internal_data::ProSHADE_data::getXToPtr
proshade_signed * getXToPtr(void)
This function allows access to the map last position along the X axis.
Definition: ProSHADE_data.cpp:3968
ProSHADE_internal_data::ProSHADE_data::getYFromPtr
proshade_signed * getYFromPtr(void)
This function allows access to the map start along the Y axis.
Definition: ProSHADE_data.cpp:3948
ProSHADE_internal_overlay::allocateTranslationFunctionMemory
void allocateTranslationFunctionMemory(fftw_complex *&tmpIn1, fftw_complex *&tmpOut1, fftw_complex *&tmpIn2, fftw_complex *&tmpOut2, fftw_complex *&resIn, fftw_complex *&resOut, fftw_plan &forwardFourierObj1, fftw_plan &forwardFourierObj2, fftw_plan &inverseFourierCombo, proshade_unsign xD, proshade_unsign yD, proshade_unsign zD)
This function allocates the memory for the Fourier transforms required for translation function compu...
Definition: ProSHADE_overlay.cpp:367
ProSHADE_internal_data::ProSHADE_data::getInvSO3Coeffs
proshade_complex * getInvSO3Coeffs(void)
This function allows access to the inverse SO(3) coefficients array.
Definition: ProSHADE_data.cpp:3826
ProSHADE_internal_data::ProSHADE_data::yFrom
proshade_signed yFrom
This is the starting index along the y axis.
Definition: ProSHADE_data.hpp:111
ProSHADE_settings::inputFiles
std::vector< std::string > inputFiles
This vector contains the filenames of all input structure files.
Definition: ProSHADE_settings.hpp:43
ProSHADE_internal_maths::complexMultiplicationConjug
void complexMultiplicationConjug(proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2, proshade_double *retReal, proshade_double *retImag)
Function to multiply two complex numbers by using the second number's conjugate.
Definition: ProSHADE_maths.cpp:62
ProSHADE_internal_data::ProSHADE_data::getYAxisOrigin
proshade_signed * getYAxisOrigin(void)
This function allows access to the map Y axis origin value.
Definition: ProSHADE_data.cpp:4008
ProSHADE_internal_overlay::computeTranslationsFromPeak
void computeTranslationsFromPeak(ProSHADE_internal_data::ProSHADE_data *staticStructure, ProSHADE_internal_data::ProSHADE_data *movingStructure, proshade_double *trsX, proshade_double *trsY, proshade_double *trsZ)
This function computes the translation in Angstroms that corresponds to the translation function peak...
Definition: ProSHADE_overlay.cpp:220
ProSHADE_internal_overlay::paddMapWithZeroes
void paddMapWithZeroes(proshade_double *oldMap, proshade_double *&newMap, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_unsign xDimIndices, proshade_unsign yDimIndices, proshade_unsign zDimIndices, proshade_unsign addXPre, proshade_unsign addYPre, proshade_unsign addZPre)
This function adds zeroes before and after the central map and copies the central map values into a n...
Definition: ProSHADE_overlay.cpp:607
ProSHADE_internal_overlay::combineFourierForTranslation
void combineFourierForTranslation(fftw_complex *tmpOut1, fftw_complex *tmpOut2, fftw_complex *&resOut, proshade_unsign xD, proshade_unsign yD, proshade_unsign zD)
This function combines Fourier coefficients of two structures in a way, so that inverse Fourier of th...
Definition: ProSHADE_overlay.cpp:432
ProSHADE_internal_distances::generateSO3CoeffsFromEMatrices
void generateSO3CoeffsFromEMatrices(ProSHADE_internal_data::ProSHADE_data *obj1, ProSHADE_internal_data::ProSHADE_data *obj2, ProSHADE_settings *settings)
This function converts the E matrices to SO(3) coefficients.
Definition: ProSHADE_distances.cpp:701
ProSHADE_internal_messages::printProgressMessage
void printProgressMessage(proshade_signed verbose, proshade_signed messageLevel, std::string message)
General stdout message printing.
Definition: ProSHADE_messages.cpp:70
ProSHADE_internal_data::ProSHADE_data::yDimSize
proshade_single yDimSize
This is the size of the map cell y dimension in Angstroms.
Definition: ProSHADE_data.hpp:60
ProSHADE_internal_mapManip::moveMapByIndices
void moveMapByIndices(proshade_single *xMov, proshade_single *yMov, proshade_single *zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed *xFrom, proshade_signed *xTo, proshade_signed *yFrom, proshade_signed *yTo, proshade_signed *zFrom, proshade_signed *zTo, proshade_signed *xOrigin, proshade_signed *yOrigin, proshade_signed *zOrigin)
Function for moving map back to original PDB location by changing the indices.
Definition: ProSHADE_mapManip.cpp:763
ProSHADE_internal_data::ProSHADE_data::zeroPaddToDims
void zeroPaddToDims(proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim)
This function changes the size of a structure to fit the supplied new limits.
Definition: ProSHADE_overlay.cpp:519
ProSHADE_internal_data::ProSHADE_data::getZFromPtr
proshade_signed * getZFromPtr(void)
This function allows access to the map start along the Z axis.
Definition: ProSHADE_data.cpp:3958
ProSHADE_internal_data::ProSHADE_data::getZToPtr
proshade_signed * getZToPtr(void)
This function allows access to the map last position along the Z axis.
Definition: ProSHADE_data.cpp:3988