ProSHADE  0.7.6.1 (AUG 2021)
Protein Shape Detection
ProSHADE_internal_mapManip Namespace Reference

This namespace contains the internal functions for manipulating maps already present in the internal structures. More...

Functions

proshade_signed myRound (proshade_double x)
 Calls the appropriate version of round function depending on compiler version. More...
 
proshade_signed myRound (proshade_single x)
 Calls the appropriate version of round function depending on compiler version. More...
 
void determinePDBRanges (gemmi::Structure pdbFile, proshade_single *xFrom, proshade_single *xTo, proshade_single *yFrom, proshade_single *yTo, proshade_single *zFrom, proshade_single *zTo, bool firstModel)
 Function for finding the PDB file ranges. More...
 
void findPDBCOMValues (gemmi::Structure pdbFile, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, bool firstModel)
 This function finds the Centre of Mass for the co-ordinate file. More...
 
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. More...
 
void rotatePDBCoordinates (gemmi::Structure *pdbFile, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double xCom, proshade_double yCom, proshade_double zCom, bool firstModel)
 Function for rotating the PDB file co-ordinates by Euler angles. More...
 
void translatePDBCoordinates (gemmi::Structure *pdbFile, proshade_double transX, proshade_double transY, proshade_double transZ, bool firstModel)
 Function for translating the PDB file co-ordinates by given distances in Angstroms. More...
 
void changePDBBFactors (gemmi::Structure *pdbFile, proshade_double newBFactorValue, bool firstModel)
 Function for changing the PDB B-factor values to a specific single value. More...
 
void removeWaters (gemmi::Structure *pdbFile, bool firstModel)
 This function removed all waters from PDB input file. More...
 
void movePDBForMapCalc (gemmi::Structure *pdbFile, proshade_single xMov, proshade_single yMov, proshade_single zMov, bool firstModel)
 Function for moving co-ordinate atoms to better suit theoretical map computation. More...
 
void generateMapFromPDB (gemmi::Structure pdbFile, proshade_double *&map, proshade_single requestedResolution, proshade_single xCell, proshade_single yCell, proshade_single zCell, proshade_signed *xTo, proshade_signed *yTo, proshade_signed *zTo, bool forceP1, bool firstModel)
 This function generates a theoretical map from co-ordinate input files. More...
 
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. More...
 
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. More...
 
void blurSharpenMap (proshade_double *&map, proshade_double *&maskedMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single blurringFactor)
 Function for blurring/sharpening maps. More...
 
void getMaskFromBlurr (proshade_double *&blurMap, proshade_double *&outMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single noIQRs)
 Function for computing mask from blurred map. More...
 
void getNonZeroBounds (proshade_double *map, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_signed *&ret)
 Function for finding the map boundaries enclosing positive only values. More...
 
void addExtraBoundSpace (proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed *&bounds, proshade_single extraSpace)
 This function takes a set of bounds and adds a given number of Angstroms to them. More...
 
void reSampleMapToResolutionTrilinear (proshade_double *&map, proshade_single resolution, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single *&corrs)
 This function re-samples a map to conform to given resolution using tri-linear interpolation. More...
 
void reSampleMapToResolutionFourier (proshade_double *&map, proshade_single resolution, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single *&corrs)
 This function re-samples a map to conform to given resolution using Fourier. More...
 
void allocateResolutionFourierMemory (fftw_complex *&origMap, fftw_complex *&fCoeffs, fftw_complex *&newFCoeffs, fftw_complex *&newMap, fftw_plan &planForwardFourier, fftw_plan &planBackwardRescaledFourier, proshade_unsign xDimOld, proshade_unsign yDimOld, proshade_unsign zDimOld, proshade_unsign xDimNew, proshade_unsign yDimNew, proshade_unsign zDimNew)
 This function allocates and checks the allocatio of the memory required by the Fourier resampling. More...
 
void releaseResolutionFourierMemory (fftw_complex *&origMap, fftw_complex *&fCoeffs, fftw_complex *&newFCoeffs, fftw_complex *&newMap, fftw_plan &planForwardFourier, fftw_plan &planBackwardRescaledFourier)
 This function releases the memory required by the Fourier resampling. More...
 
void changeFourierOrder (fftw_complex *&fCoeffs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, bool negativeFirst)
 This function changes the order of Fourier coefficients in a 3D array between positive first (default) and negative first (mass centered at xMax/2, yMax/2, zMax/2 instead of 0,0,0) More...
 
void removeMapPhase (fftw_complex *&mapCoeffs, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim)
 This function removes the phase from reciprocal (frequency) map. More...
 
void getFakeHalfMap (proshade_double *&map, proshade_double *&fakeHalfMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_signed fakeMapKernel)
 Function for creating "fake" half-maps. More...
 
void getCorrelationMapMask (proshade_double *&map, proshade_double *&fakeHalfMap, proshade_double *&correlationMask, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_signed corrMaskKernel)
 Function for creating the correlation mask. More...
 
proshade_single getIndicesFromAngstroms (proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single dist)
 This function converts distance in Angstroms to distance in map indices. More...
 
void connectMaskBlobs (proshade_double *&mask, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single maskThres)
 This function connects blobs in mask. More...
 
void beautifyBoundaries (proshade_signed *&bounds, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_signed boundsDiffThres)
 Function for modifying boundaries to a mathematically more pleasant values. More...
 
proshade_signed betterClosePrimeFactors (proshade_signed fromRange, proshade_signed toRange)
 Function for finding close numbers with better prime factors. More...
 
void distributeSpaceToBoundaries (proshade_signed &minBound, proshade_signed &maxBound, proshade_signed oldBoundRange, proshade_signed newBoundRange)
 Function for adding space to boundaries within the map confines. More...
 
void copyMapByBounds (proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo, proshade_signed origXFrom, proshade_signed origYFrom, proshade_signed origZFrom, proshade_unsign yDimIndices, proshade_unsign zDimIndices, proshade_unsign origXDimIndices, proshade_unsign origYDimIndices, proshade_unsign origZDimIndices, proshade_double *&newMap, proshade_double *origMap)
 This function copies an old map to a new map with different boundaries. More...
 

Detailed Description

This namespace contains the internal functions for manipulating maps already present in the internal structures.

The ProSHADE_internal_mapManip namespace contains helper functions for map manipulation and processing. However, these functions do make minimum assumptions about the map internal organisation and variables and simply perform tasks on generic variables. None of these functions should be used directly be the user.

Function Documentation

◆ addExtraBoundSpace()

void ProSHADE_internal_mapManip::addExtraBoundSpace ( proshade_unsign  xDim,
proshade_unsign  yDim,
proshade_unsign  zDim,
proshade_single  xAngs,
proshade_single  yAngs,
proshade_single  zAngs,
proshade_signed *&  bounds,
proshade_single  extraSpace 
)

This function takes a set of bounds and adds a given number of Angstroms to them.

This function adds an extra distance to a set of bounds as given by the number of angstroms to be added. If the resulting bounds are larger than the maximum map bounds, the maxima and minima are returned instead, thus partially ignoring the requested extra space. !!! If the resulting index ranges would be odd, extra one idice is added to make them even. !!!

Parameters
[in]xDimThe number of indices along the x axis of the map.
[in]yDimThe number of indices along the y axis of the map.
[in]zDimThe number of indices along the z axis of the map.
[in]xAngsThe size of the x dimension of the map in angstroms.
[in]yAngsThe size of the y dimension of the map in angstroms.
[in]zAngsThe size of the z dimension of the map in angstroms.
[in]boundsA proshade_unsign pointer reference to array of 6 which has the current bounds (0 = minX; 1 = maxX; 2 = minY; 3 = maxY; 4 - minZ; 5 = maxZ).
[in]extraSpaceThe number of angstroms to be added to each dimension (both to start and to end).

Definition at line 1139 of file ProSHADE_mapManip.cpp.

1140 {
1141  //================================================ Convert angstrom distance to indices
1142  proshade_signed xExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( xAngs / static_cast<proshade_single> ( xDim ) ) );
1143  proshade_signed yExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( yAngs / static_cast<proshade_single> ( yDim ) ) );
1144  proshade_signed zExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( zAngs / static_cast<proshade_single> ( zDim ) ) );
1145 
1146  //=============================================== Changed bounds even if exceeding physical map - this will be deal with in the map creation part
1147  bounds[0] = bounds[0] - xExtraInds;
1148  bounds[1] = bounds[1] + xExtraInds;
1149  bounds[2] = bounds[2] - yExtraInds;
1150  bounds[3] = bounds[3] + yExtraInds;
1151  bounds[4] = bounds[4] - zExtraInds;
1152  bounds[5] = bounds[5] + zExtraInds;
1153 
1154  //================================================ Done
1155  return ;
1156 
1157 }

◆ allocateResolutionFourierMemory()

void ProSHADE_internal_mapManip::allocateResolutionFourierMemory ( fftw_complex *&  origMap,
fftw_complex *&  fCoeffs,
fftw_complex *&  newFCoeffs,
fftw_complex *&  newMap,
fftw_plan &  planForwardFourier,
fftw_plan &  planBackwardRescaledFourier,
proshade_unsign  xDimOld,
proshade_unsign  yDimOld,
proshade_unsign  zDimOld,
proshade_unsign  xDimNew,
proshade_unsign  yDimNew,
proshade_unsign  zDimNew 
)

This function allocates and checks the allocatio of the memory required by the Fourier resampling.

This function allocates the memory required for the Fourier space re-sampling of density maps. It allocates the original map and original map coefficients arrays (these need to be of the fftw_complex type) as well as the re-sampled map and re-sampled map coefficients. It then also proceeds to check the memory allocation and creating the FFTW transform plans for both, the forward and re-sampled backward operations.

Parameters
[in]origMapA Reference pointer to an array where the original map data will be stored.
[in]fCoeffsA Reference pointer to an array where the original map Fourier coefficients data will be stored.
[in]newFCoeffsA Reference pointer to an array where the re-sampled map Fourier coefficients data will be stored.
[in]newMapA Reference pointer to an array where the re-sampled map data will be stored.
[in]planForwardFourierFFTW_plan which will compute the original map to Fourier coefficients transform.
[in]planBackwardRescaledFourierFFTW_plan which will compute the re-sampled Fourier coefficients to re-sampled map transform.
[in]xDimOldThe number of indices along the x-axis of the original map.
[in]yDimOldThe number of indices along the y-axis of the original map.
[in]zDimOldThe number of indices along the z-axis of the original map.
[in]xDimNewThe number of indices along the x-axis of the re-sampled map.
[in]yDimNewThe number of indices along the y-axis of the re-sampled map.
[in]zDimNewThe number of indices along the z-axis of the re-sampled map.

Definition at line 1496 of file ProSHADE_mapManip.cpp.

1497 {
1498  //================================================ Initialise memory
1499  origMap = new fftw_complex [xDimOld * yDimOld * zDimOld];
1500  fCoeffs = new fftw_complex [xDimOld * yDimOld * zDimOld];
1501  newFCoeffs = new fftw_complex [xDimNew * yDimNew * zDimNew];
1502  newMap = new fftw_complex [xDimNew * yDimNew * zDimNew];
1503 
1504  //================================================ Check memory allocation
1505  ProSHADE_internal_misc::checkMemoryAllocation ( origMap, __FILE__, __LINE__, __func__ );
1506  ProSHADE_internal_misc::checkMemoryAllocation ( fCoeffs, __FILE__, __LINE__, __func__ );
1507  ProSHADE_internal_misc::checkMemoryAllocation ( newFCoeffs, __FILE__, __LINE__, __func__ );
1508  ProSHADE_internal_misc::checkMemoryAllocation ( newMap, __FILE__, __LINE__, __func__ );
1509 
1510  //================================================ Create plans
1511  planForwardFourier = fftw_plan_dft_3d ( static_cast< int > ( xDimOld ), static_cast< int > ( yDimOld ), static_cast< int > ( zDimOld ), origMap, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
1512  planBackwardRescaledFourier = fftw_plan_dft_3d ( static_cast< int > ( xDimNew ), static_cast< int > ( yDimNew ), static_cast< int > ( zDimNew ), newFCoeffs, newMap, FFTW_BACKWARD, FFTW_ESTIMATE );
1513 
1514  //================================================ Done
1515  return ;
1516 
1517 }

◆ beautifyBoundaries()

void ProSHADE_internal_mapManip::beautifyBoundaries ( proshade_signed *&  bounds,
proshade_unsign  xDim,
proshade_unsign  yDim,
proshade_unsign  zDim,
proshade_signed  boundsDiffThres 
)

Function for modifying boundaries to a mathematically more pleasant values.

This function does two things. Firstly, it attempts to find even boundaries which have a small sum of the prime factors - i.e. this is good for further processing. Secondly, it attempts to make boundaries with sizes within some threshold the same. This is purely human thing :-).

Parameters
[in]boundsA proshade_signed pointer to array of 6 in which the current bounds are (0 = minX; 1 = maxX; 2 = minY; 3 = maxY; 4 - minZ; 5 = maxZ).
[in]xDimThe number of indices along the x axis of the map.
[in]yDimThe number of indices along the y axis of the map.
[in]zDimThe number of indices along the z axis of the map.
[in]boundsDiffThresNumber of indices by which boudaries in different dimensions can differ and still be made the same.

Definition at line 1916 of file ProSHADE_mapManip.cpp.

1917 {
1918  //================================================ If extra bounds space pushed the bounds over the physical map, freely add up to 10 indices over the current value
1919  while ( bounds[1] >= static_cast<proshade_signed> ( xDim ) ) { xDim += 10; }
1920  while ( bounds[3] >= static_cast<proshade_signed> ( yDim ) ) { yDim += 10; }
1921  while ( bounds[5] >= static_cast<proshade_signed> ( zDim ) ) { zDim += 10; }
1922 
1923  //================================================ Find if better bouds exist in terms of prime numbers
1924  proshade_signed addToX = betterClosePrimeFactors ( bounds[1] - bounds[0] + 1, static_cast< proshade_signed > ( xDim ) );
1925  proshade_signed addToY = betterClosePrimeFactors ( bounds[3] - bounds[2] + 1, static_cast< proshade_signed > ( yDim ) );
1926  proshade_signed addToZ = betterClosePrimeFactors ( bounds[5] - bounds[4] + 1, static_cast< proshade_signed > ( zDim ) );
1927 
1928  //================================================ Figure if similar sizes should not be forced to be the same
1929  proshade_signed XtoY = std::abs ( addToX - addToY );
1930  proshade_signed XtoZ = std::abs ( addToX - addToZ );
1931  proshade_signed YtoZ = std::abs ( addToY - addToZ );
1932 
1933  if ( ( ( XtoY < boundsDiffThres ) && ( XtoZ < boundsDiffThres ) ) ||
1934  ( ( XtoY < boundsDiffThres ) && ( YtoZ < boundsDiffThres ) ) ||
1935  ( ( XtoZ < boundsDiffThres ) && ( YtoZ < boundsDiffThres ) ) )
1936  {
1937  //============================================ Dealing with chain ( a <= b <= c type ) where at least two dista are smaller than threshold
1938  proshade_signed maxSize = std::max ( addToX, std::max ( addToY, addToZ ) );
1939  addToX = maxSize;
1940  addToY = maxSize;
1941  addToZ = maxSize;
1942  }
1943  else
1944  {
1945  //============================================ Only two may be similar enough
1946  if ( XtoY <= boundsDiffThres )
1947  {
1948  proshade_signed maxSize = std::max ( addToX, addToY );
1949  addToX = maxSize;
1950  addToY = maxSize;
1951  }
1952  if ( XtoZ <= boundsDiffThres )
1953  {
1954  proshade_signed maxSize = std::max ( addToX, addToZ );
1955  addToX = maxSize;
1956  addToZ = maxSize;
1957  }
1958  if ( YtoZ <= boundsDiffThres )
1959  {
1960  proshade_signed maxSize = std::max ( addToY, addToZ );
1961  addToY = maxSize;
1962  addToZ = maxSize;
1963  }
1964  }
1965 
1966  //================================================ Add the extra space appropriately
1967  distributeSpaceToBoundaries ( bounds[0], bounds[1], ( bounds[1] - bounds[0] + 1 ), addToX );
1968  distributeSpaceToBoundaries ( bounds[2], bounds[3], ( bounds[3] - bounds[2] + 1 ), addToY );
1969  distributeSpaceToBoundaries ( bounds[4], bounds[5], ( bounds[5] - bounds[4] + 1 ), addToZ );
1970 
1971  //================================================ Done
1972  return ;
1973 
1974 }

◆ betterClosePrimeFactors()

proshade_signed ProSHADE_internal_mapManip::betterClosePrimeFactors ( proshade_signed  fromRange,
proshade_signed  toRange 
)

Function for finding close numbers with better prime factors.

This function takes a range of integer numbers and proceeds to find if a close number in the range to the range start does not have better (smaller sum of) prime factors, taking into account the distance of this better number to the starting position. This is useful for finding boundaries with nice prime number sizes.

Parameters
[in]fromRangeWhich number the search will be from and possibilities compared to.
[in]toRangeThe maximum number to which the posibilities should be searched.
[out]XThe value in the range which has better prime factors, or the first value if none other has.

Definition at line 1986 of file ProSHADE_mapManip.cpp.

1987 {
1988  //================================================ Initialise variables
1989  proshade_signed ret = fromRange;
1990  std::vector < proshade_signed > posibles, hlp;
1991  proshade_signed sum;
1992 
1993  //================================================ Find the sums of prime factors for each number in the whole range
1994  for ( proshade_signed iter = fromRange; iter < toRange; iter++ )
1995  {
1996  sum = 0;
1998  for ( proshade_unsign i = 0; i < static_cast<proshade_unsign> ( hlp.size() ); i++ ) { sum += hlp.at(i); }
1999  hlp.clear ( );
2000  ProSHADE_internal_misc::addToSignedVector ( &posibles, sum );
2001  }
2002 
2003  //================================================ Find a better number
2004  for ( proshade_signed iter = fromRange; iter < toRange; iter++ )
2005  {
2006  //============================================ Ignore odd numbers
2007  if ( iter %2 != 0 ) { continue; }
2008 
2009  //============================================ Better number needs to have lower prime factor sum, but also needs not to be too far
2010  if ( posibles.at( static_cast< size_t > ( iter - fromRange ) ) < ( posibles.at( static_cast< size_t > ( ret - fromRange ) ) - ( iter - ret ) ) ) { ret = iter; }
2011  }
2012 
2013  //================================================ In the case of large prime number, just add one for even number
2014  if ( ( ret % 2 != 0 ) && ( ret < ( toRange - 1 ) ) ) { ret += 1; }
2015 
2016  //================================================ Done
2017  return ( ret );
2018 
2019 }

◆ blurSharpenMap()

void ProSHADE_internal_mapManip::blurSharpenMap ( proshade_double *&  map,
proshade_double *&  blurredMap,
proshade_unsign  xDimS,
proshade_unsign  yDimS,
proshade_unsign  zDimS,
proshade_single  xAngs,
proshade_single  yAngs,
proshade_single  zAngs,
proshade_single  blurringFactor 
)

Function for blurring/sharpening maps.

This function takes the internal map representation information from the parameters as well as the internal map itself and proceeds to compute blue/sharpen the map by changing its overall B-factors by adding a specific value. If this value is negative, then the map is sharpened, while if the value is positive, the map is blurred.

Parameters
[in]mapA Reference Pointer to the map which should be blurred/sharpened.
[in]blurredMapA Reference Pointer to the variable which will store the modified map.
[in]xDimSThe number of indices along the x axis of the map.
[in]yDimSThe number of indices along the y axis of the map.
[in]zDimSThe number of indices along the z axis of the map.
[in]xAngsThe size of the x dimension of the map in angstroms.
[in]yAngsThe size of the y dimension of the map in angstroms.
[in]zAngsThe size of the z dimension of the map in angstroms.
[in]blurringFactorThe amount of B-factor change that should be applied to the map. (I.e. increasing the map overall B-factors by 20 will be done by supplying 20 as a value for this parameter).

Definition at line 932 of file ProSHADE_mapManip.cpp.

933 {
934  //================================================ Set local variables
935  proshade_signed xDim = static_cast< proshade_signed > ( xDimS );
936  proshade_signed yDim = static_cast< proshade_signed > ( yDimS );
937  proshade_signed zDim = static_cast< proshade_signed > ( zDimS );
938  proshade_double real, imag, S, mag, phase;
939  proshade_signed h, k, l;
940  proshade_unsign arrayPos = 0;
941  proshade_double normFactor = static_cast<proshade_double> ( xDim * yDim * zDim );
942 
943  //================================================ Copy map for processing
944  fftw_complex* mapCoeffs = new fftw_complex[xDim * yDim * zDim];
945  fftw_complex* mapMask = new fftw_complex[xDim * yDim * zDim];
946 
947  //================================================ Check memory allocation
948  ProSHADE_internal_misc::checkMemoryAllocation ( mapCoeffs, __FILE__, __LINE__, __func__ );
949  ProSHADE_internal_misc::checkMemoryAllocation ( mapMask, __FILE__, __LINE__, __func__ );
950 
951  //================================================ Copy data to mask
952  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> (xDim * yDim * zDim); iter++ )
953  {
954  mapMask[iter][0] = map[iter];
955  mapMask[iter][1] = 0.0;
956  }
957 
958  //================================================ Prepare FFTW plans
959  fftw_plan forward = fftw_plan_dft_3d ( static_cast< int > ( xDim ), static_cast< int > ( yDim ), static_cast< int > ( zDim ), mapMask, mapCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
960  fftw_plan inverse = fftw_plan_dft_3d ( static_cast< int > ( xDim ), static_cast< int > ( yDim ), static_cast< int > ( zDim ), mapCoeffs, mapMask, FFTW_BACKWARD, FFTW_ESTIMATE );
961 
962  //================================================ Run forward Fourier
963  fftw_execute ( forward );
964 
965  //================================================ Blur the coeffs
966  for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
967  {
968  for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
969  {
970  for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
971  {
972  //==================================== Var init
973  arrayPos = wIt + static_cast< proshade_unsign > ( zDim ) * ( vIt + static_cast< proshade_unsign > ( yDim ) * uIt );
974  real = mapCoeffs[arrayPos][0];
975  imag = mapCoeffs[arrayPos][1];
976 
977  //==================================== Convert to HKL
978  if ( uIt > static_cast< proshade_unsign > ( (xDim+1) / 2) ) { h = static_cast < proshade_signed > ( uIt ) - xDim; } else { h = static_cast < proshade_signed > ( uIt ); }
979  if ( vIt > static_cast< proshade_unsign > ( (yDim+1) / 2) ) { k = static_cast < proshade_signed > ( vIt ) - yDim; } else { k = static_cast < proshade_signed > ( vIt ); }
980  if ( wIt > static_cast< proshade_unsign > ( (zDim+1) / 2) ) { l = static_cast < proshade_signed > ( wIt ) - zDim; } else { l = static_cast < proshade_signed > ( wIt ); }
981 
982  //====================================Get magnitude and phase with mask parameters
983  S = ( pow( static_cast< proshade_double > ( h ) / static_cast< proshade_double > ( xAngs ), 2.0 ) +
984  pow( static_cast< proshade_double > ( k ) / static_cast< proshade_double > ( yAngs ), 2.0 ) +
985  pow( static_cast< proshade_double > ( l ) / static_cast< proshade_double > ( zAngs ), 2.0 ) );
986  mag = std::sqrt ( (real*real) + (imag*imag) ) * std::exp ( - ( ( static_cast< proshade_double > ( blurringFactor ) * S ) / 4.0 ) );
987  phase = std::atan2 ( imag, real );
988 
989  //==================================== Save the mask data
990  mapCoeffs[arrayPos][0] = ( mag * cos(phase) ) / normFactor;
991  mapCoeffs[arrayPos][1] = ( mag * sin(phase) ) / normFactor;
992  }
993  }
994  }
995 
996  //================================================ Run inverse Fourier
997  fftw_execute ( inverse );
998 
999  //================================================ Save the results
1000  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> (xDim * yDim * zDim); iter++ )
1001  {
1002  blurredMap[iter] = mapMask[iter][0];
1003  }
1004 
1005  //================================================ Release memory
1006  delete[] mapMask;
1007  delete[] mapCoeffs;
1008 
1009  //================================================ Delete FFTW plans
1010  fftw_destroy_plan ( forward );
1011  fftw_destroy_plan ( inverse );
1012 
1013  //================================================ Done
1014  return ;
1015 
1016 }

◆ changeFourierOrder()

void ProSHADE_internal_mapManip::changeFourierOrder ( fftw_complex *&  fCoeffs,
proshade_signed  xDim,
proshade_signed  yDim,
proshade_signed  zDim,
bool  negativeFirst 
)

This function changes the order of Fourier coefficients in a 3D array between positive first (default) and negative first (mass centered at xMax/2, yMax/2, zMax/2 instead of 0,0,0)

This function firstly determines the start and end of the positive and negative Fourier coefficients for all three axes (it assumes 3D array); this works for both odd and even axis sizes. To do this properly, in the case of odd axis sizes, the function needs to know whether we are changing the order from FFTW defaul positive first, or if we are changing the other way around - the last parameter serves the purpose of passing this information. Finally, the function then switches the order of the Fourier coefficients as requested using a temporary array that it allocates and deletes within itself.

Parameters
[in]fCoeffsA Reference pointer to an array where the Fourier coefficients for re-ordering are stored.
[in]xDimThe size of the x-axis dimension of the fCoeffs array.
[in]yDimThe size of the x-axis dimension of the fCoeffs array.
[in]zDimThe size of the x-axis dimension of the fCoeffs array.
[in]negativeFirstShould the coefficients be stored negarive first (TRUE), or are we reversing already re-ordered array (FALSE)?

Definition at line 1559 of file ProSHADE_mapManip.cpp.

1560 {
1561  //================================================ Initialise local variables
1562  proshade_signed h = 0, k = 0, l = 0, origSizeArr = 0, newSizeArr = 0;
1563  proshade_signed xSeq1FreqStart, ySeq1FreqStart, zSeq1FreqStart, xSeq2FreqStart, ySeq2FreqStart, zSeq2FreqStart;
1564 
1565  //================================================ Find the positive and negative indices cot-offs
1566  if ( negativeFirst )
1567  {
1568  if ( ( xDim % 2 ) == 0 ) { xSeq1FreqStart = xDim / 2; xSeq2FreqStart = xDim / 2; } else { xSeq1FreqStart = (xDim / 2) + 1; xSeq2FreqStart = xDim / 2; }
1569  if ( ( yDim % 2 ) == 0 ) { ySeq1FreqStart = yDim / 2; ySeq2FreqStart = yDim / 2; } else { ySeq1FreqStart = (yDim / 2) + 1; ySeq2FreqStart = yDim / 2; }
1570  if ( ( zDim % 2 ) == 0 ) { zSeq1FreqStart = zDim / 2; zSeq2FreqStart = zDim / 2; } else { zSeq1FreqStart = (zDim / 2) + 1; zSeq2FreqStart = zDim / 2; }
1571  }
1572  else
1573  {
1574  if ( ( xDim % 2 ) == 0 ) { xSeq1FreqStart = xDim / 2; xSeq2FreqStart = xDim / 2; } else { xSeq1FreqStart = (xDim / 2); xSeq2FreqStart = xDim / 2 + 1; }
1575  if ( ( yDim % 2 ) == 0 ) { ySeq1FreqStart = yDim / 2; ySeq2FreqStart = yDim / 2; } else { ySeq1FreqStart = (yDim / 2); ySeq2FreqStart = yDim / 2 + 1; }
1576  if ( ( zDim % 2 ) == 0 ) { zSeq1FreqStart = zDim / 2; zSeq2FreqStart = zDim / 2; } else { zSeq1FreqStart = (zDim / 2); zSeq2FreqStart = zDim / 2 + 1; }
1577  }
1578 
1579  //================================================ Allocate helper array memory
1580  fftw_complex *hlpFCoeffs = new fftw_complex [xDim * yDim * zDim];
1581  ProSHADE_internal_misc::checkMemoryAllocation ( hlpFCoeffs, __FILE__, __LINE__, __func__ );
1582 
1583  //================================================ Change the coefficients order
1584  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1585  {
1586  //============================================ Find x frequency
1587  if ( xIt < xSeq1FreqStart ) { h = xIt + xSeq2FreqStart; } else { h = xIt - xSeq1FreqStart; }
1588  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1589  {
1590  //======================================== Find y frequency
1591  if ( yIt < ySeq1FreqStart ) { k = yIt + ySeq2FreqStart; } else { k = yIt - ySeq1FreqStart; }
1592 
1593  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1594  {
1595  //==================================== Find z frequency
1596  if ( zIt < zSeq1FreqStart ) { l = zIt + zSeq2FreqStart; } else { l = zIt - zSeq1FreqStart; }
1597 
1598  //==================================== Find array positions
1599  newSizeArr = l + zDim * ( k + yDim * h );
1600  origSizeArr = zIt + zDim * ( yIt + yDim * xIt );
1601 
1602  //==================================== Copy vals
1603  hlpFCoeffs[newSizeArr][0] = fCoeffs[origSizeArr][0];
1604  hlpFCoeffs[newSizeArr][1] = fCoeffs[origSizeArr][1];
1605  }
1606  }
1607  }
1608 
1609  //================================================ Copy the helper array to the input Fourier coefficients array
1610  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { fCoeffs[iter][0] = hlpFCoeffs[iter][0]; fCoeffs[iter][1] = hlpFCoeffs[iter][1]; }
1611 
1612  //================================================ Release helper array memory
1613  delete[] hlpFCoeffs;
1614 
1615  //================================================ Done
1616  return ;
1617 
1618 }

◆ changePDBBFactors()

void ProSHADE_internal_mapManip::changePDBBFactors ( gemmi::Structure *  pdbFile,
proshade_double  newBFactorValue,
bool  firstModel 
)

Function for changing the PDB B-factor values to a specific single value.

This function does a quick read-through the PDB file and changes all the atom B-factor values to a single, pre-set value. This is important for PDB files which have all atoms with B-factor value 0, as these then produce bad maps, which fail further processing by ProSHADE.

Parameters
[in]pdbFileA pointer to gemmi::Structure object read in from the input file.
[in]newBFactorValueThe value with which all atom B-factor values should be replaced.
[in]firstModelShould only the first, or all models be used?

Definition at line 444 of file ProSHADE_mapManip.cpp.

445 {
446  //================================================ Use the first model, if it exists
447  if ( pdbFile->models.size() > 0 )
448  {
449  //============================================ For each model
450  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
451  {
452  //======================================== Check if multiple models are allowed
453  if ( firstModel && ( sIt != 0 ) ) { break; }
454 
455  //======================================== Get model
456  gemmi::Model *model = &pdbFile->models.at(sIt);
457 
458  //======================================== For each chain
459  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
460  {
461  //==================================== Get chain
462  gemmi::Chain *chain = &model->chains.at(mIt);
463 
464  //==================================== For each residue
465  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
466  {
467  //================================ Get residue
468  gemmi::Residue *residue = &chain->residues.at(rIt);
469 
470  //================================ For each atom
471  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
472  {
473  //============================ Get atom
474  gemmi::Atom *atom = &residue->atoms.at(aIt);
475 
476  //============================ Change the B-factors
477  atom->b_iso = static_cast< float > ( newBFactorValue );
478  }
479  }
480  }
481  }
482  }
483  else
484  {
485  std::stringstream hlpSS;
486  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
487  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
488  }
489 
490  //================================================ Done
491  return ;
492 
493 }

◆ connectMaskBlobs()

void ProSHADE_internal_mapManip::connectMaskBlobs ( proshade_double *&  mask,
proshade_signed  xDim,
proshade_signed  yDim,
proshade_signed  zDim,
proshade_single  xAngs,
proshade_single  yAngs,
proshade_single  zAngs,
proshade_single  maskThres 
)

This function connects blobs in mask.

Parameters
[in]maskA pointer reference to mask map in which blobs should be connected.
[in]xDimThe number of map indices along the x-axis.
[in]yDimThe number of map indices along the y-axis.
[in]zDimThe number of map indices along the z-axis.
[in]xAngsThe map size in Angstroms along the x-axis.
[in]yAngsThe map size in Angstroms along the y-axis.
[in]zAngsThe map size in Angstroms along the z-axis.
[in]maskThresThe threshold which will be used for applying mask.

Definition at line 1839 of file ProSHADE_mapManip.cpp.

1840 {
1841  //================================================ Initialise variables
1842  proshade_double* hlpMap = new proshade_double[xDim * yDim * zDim];
1843  proshade_signed addSurroundingPoints = static_cast< proshade_signed > ( std::max ( 3L, static_cast<proshade_signed> ( std::ceil ( getIndicesFromAngstroms( static_cast< proshade_unsign > ( xDim ), static_cast< proshade_unsign > ( yDim ), static_cast< proshade_unsign > ( zDim ), xAngs, yAngs, zAngs, static_cast< proshade_single > ( std::max( xAngs, std::max( yAngs, zAngs ) ) * 0.1f ) ) ) ) ) );
1844  proshade_signed currPos, neighXPos, neighYPos, neighZPos, neighArrPos;
1845 
1846  //================================================ Check memory allocation
1847  ProSHADE_internal_misc::checkMemoryAllocation ( hlpMap, __FILE__, __LINE__, __func__ );
1848 
1849  //================================================ Copy the mask map
1850  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { hlpMap[iter] = mask[iter]; }
1851 
1852  //================================================ Repeat as many times as needed
1853  for ( proshade_signed it = 0; it < addSurroundingPoints; it++ )
1854  {
1855  //============================================ For each x, y and z
1856  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1857  {
1858  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1859  {
1860  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1861  {
1862  //================================ Current position
1863  currPos = zIt + zDim * ( yIt + yDim * xIt );
1864 
1865  //================================ If zero, ignore
1866  if ( hlpMap[currPos] < static_cast< proshade_double > ( maskThres ) ) { continue; }
1867 
1868  //================================ Check neighbours
1869  for ( proshade_signed xCh = -1; xCh <= +1; xCh++ )
1870  {
1871  for ( proshade_signed yCh = -1; yCh <= +1; yCh++ )
1872  {
1873  for ( proshade_signed zCh = -1; zCh <= +1; zCh++ )
1874  {
1875  if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) { continue; }
1876 
1877  //==================== Find the nieghbour indices (without periodicity)
1878  neighXPos = xIt + xCh; if ( neighXPos < 0 ) { continue; } if ( neighXPos >= xDim ) { continue; }
1879  neighYPos = yIt + yCh; if ( neighYPos < 0 ) { continue; } if ( neighYPos >= yDim ) { continue; }
1880  neighZPos = zIt + zCh; if ( neighZPos < 0 ) { continue; } if ( neighZPos >= zDim ) { continue; }
1881  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1882 
1883  //==================== Add to mask if this point is below it (as it is a neighbour to a point which is part of the mask)
1884  if ( hlpMap[neighArrPos] < static_cast< proshade_double > ( maskThres ) ) { mask[neighArrPos] = static_cast< proshade_double > ( maskThres ); }
1885  }
1886  }
1887  }
1888  }
1889  }
1890  }
1891 
1892  //============================================ Now copy the updated mask to the temporary helper, unless last iteration
1893  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { hlpMap[iter] = mask[iter]; }
1894  }
1895 
1896  //================================================ Release memory
1897  delete[] hlpMap;
1898 
1899  //================================================ Done
1900  return ;
1901 
1902 }

◆ copyMapByBounds()

void ProSHADE_internal_mapManip::copyMapByBounds ( proshade_signed  xFrom,
proshade_signed  xTo,
proshade_signed  yFrom,
proshade_signed  yTo,
proshade_signed  zFrom,
proshade_signed  zTo,
proshade_signed  origXFrom,
proshade_signed  origYFrom,
proshade_signed  origZFrom,
proshade_unsign  yDimIndices,
proshade_unsign  zDimIndices,
proshade_unsign  origXDimIndices,
proshade_unsign  origYDimIndices,
proshade_unsign  origZDimIndices,
proshade_double *&  newMap,
proshade_double *  origMap 
)

This function copies an old map to a new map with different boundaries.

This function takes old and new structure details and the original structure map. It then proceed to iterate through the new structure map, checking if this particular point is inside the original structure's map and if it is, it will copy the value of this point from the original to the new map. Otherwise, it will set the new map's point value to 0.0. This is used when creating a new structure from an old structure given new bounds, which can be larger or smaller than the originals.

Parameters
[in]xFromThe starting x-axis index of the new map.
[in]xToThe last x-axis index of the new map.
[in]yFromThe starting y-axis index of the new map.
[in]yToThe last y-axis index of the new map.
[in]zFromThe starting z-axis index of the new map.
[in]zToThe last z-axis index of the new map.
[in]origXFromThe starting x-axis index of the original (copied) map.
[in]origYFromThe starting y-axis index of the original (copied) map.
[in]origZFromThe starting z-axis index of the original (copied) map.
[in]yDimIndicesThe size of the y-axis dimension in the new map.
[in]zDimIndicesThe size of the z-axis dimension in the new map.
[in]origXDimIndicesThe size of the x-axis dimension in the old map.
[in]origYDimIndicesThe size of the y-axis dimension in the old map.
[in]origZDimIndicesThe size of the z-axis dimension in the old map.
[in]newMapPointer reference to the array where new map will be saved to.
[in]origMapPointer to the array where the original map can be accessed.

Definition at line 2082 of file ProSHADE_mapManip.cpp.

2083 {
2084  //================================================ Initialise variables
2085  proshade_signed newMapIndex, oldMapIndex, oldX, oldY, oldZ, newX, newY, newZ;
2086 
2087  //=============================================== For all values in the new map
2088  for ( proshade_signed xIt = xFrom; xIt <= xTo; xIt++ )
2089  {
2090  //============================================ Index position init
2091  newX = ( xIt - xFrom );
2092  oldX = ( newX + ( xFrom - origXFrom ) );
2093 
2094  for ( proshade_signed yIt = yFrom; yIt <= yTo; yIt++ )
2095  {
2096  //======================================== Index position init
2097  newY = ( yIt - yFrom );
2098  oldY = ( newY + ( yFrom - origYFrom ) );
2099 
2100  for ( proshade_signed zIt = zFrom; zIt <= zTo; zIt++ )
2101  {
2102  //==================================== Index position init
2103  newZ = ( zIt - zFrom );
2104  oldZ = ( newZ + ( zFrom - origZFrom ) );
2105 
2106  //==================================== Index arrays
2107  newMapIndex = newZ + static_cast< proshade_signed > ( zDimIndices ) * ( newY + static_cast< proshade_signed > ( yDimIndices ) * newX );
2108  oldMapIndex = oldZ + static_cast< proshade_signed > ( origZDimIndices ) * ( oldY + static_cast< proshade_signed > ( origYDimIndices ) * oldX );
2109 
2110  //============================ Check if we are adding new point
2111  if ( ( ( oldX < 0 ) || ( oldX >= static_cast< proshade_signed > ( origXDimIndices ) ) ) ||
2112  ( ( oldY < 0 ) || ( oldY >= static_cast< proshade_signed > ( origYDimIndices ) ) ) ||
2113  ( ( oldZ < 0 ) || ( oldZ >= static_cast< proshade_signed > ( origZDimIndices ) ) ) )
2114  {
2115  //================================ Yes, this is a new point, no known value for it
2116  newMap[newMapIndex] = 0.0;
2117  continue;
2118  }
2119 
2120  //==================================== If we got here, this should be within the old map, so just copy value
2121  newMap[newMapIndex] = origMap[oldMapIndex];
2122  }
2123  }
2124  }
2125 
2126  //================================================ Done
2127  return ;
2128 
2129 }

◆ determinePDBRanges()

void ProSHADE_internal_mapManip::determinePDBRanges ( gemmi::Structure  pdbFile,
proshade_single *  xFrom,
proshade_single *  xTo,
proshade_single *  yFrom,
proshade_single *  yTo,
proshade_single *  zFrom,
proshade_single *  zTo,
bool  firstModel 
)

Function for finding the PDB file ranges.

This function does a quick read-through the PDB file and reports the x, y and z to and from values. This is used to determine if these need to be changed for proper theoretical map computation operation.

Parameters
[in]pdbFileA gemmi::Structure object read in from the input file.
[in]xFromAddress to a variable to save the x axis minimal atom position.
[in]xToAddress to a variable to save the x axis maximum atom position.
[in]yFromAddress to a variable to save the y axis minimal atom position.
[in]yToAddress to a variable to save the y axis maximum atom position.
[in]zFromAddress to a variable to save the z axis minimal atom position.
[in]zToAddress to a variable to save the z axis maximum atom position.
[in]firstModelShould only the first, or all models be used?
Warning
This function ignores hydrogen atoms, as they will not be used in map computation by Gemmi and their inclusion causes mismatches between the map and the co-ordinates as they now contain different contents.

Definition at line 72 of file ProSHADE_mapManip.cpp.

73 {
74  //================================================ Initialise structure crawl
75  bool firstAtom = true;
76 
77  //================================================ Use the first model, if it exists
78  if ( pdbFile.models.size() > 0 )
79  {
80  //============================================ For each model
81  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); sIt++ )
82  {
83  //======================================== Check if multiple models are allowed
84  if ( firstModel && ( sIt != 0 ) ) { break; }
85 
86  //======================================== Get model
87  gemmi::Model model = pdbFile.models.at(sIt);
88 
89  //======================================== For each chain
90  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model.chains.size() ); mIt++ )
91  {
92  //==================================== Get chain
93  gemmi::Chain chain = model.chains.at(mIt);
94 
95  //==================================== For each residue
96  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain.residues.size() ); rIt++ )
97  {
98  //================================ Get residue
99  gemmi::Residue residue = chain.residues.at(rIt);
100 
101  //================================ For each atom
102  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue.atoms.size() ); aIt++ )
103  {
104  //============================ Get atom
105  gemmi::Atom atom = residue.atoms.at(aIt);
106 
107  //============================ Ignore hydrogens, map computations ignore them anyway and inclusion here causes map - co-ordinate mismatches.
108  if ( atom.is_hydrogen() ) { continue; }
109 
110  //============================ Find the coordinate ranges
111  if ( firstAtom )
112  {
113  *xTo = static_cast<proshade_single> ( atom.pos.x );
114  *xFrom = static_cast<proshade_single> ( atom.pos.x );
115  *yTo = static_cast<proshade_single> ( atom.pos.y );
116  *yFrom = static_cast<proshade_single> ( atom.pos.y );
117  *zTo = static_cast<proshade_single> ( atom.pos.z );
118  *zFrom = static_cast<proshade_single> ( atom.pos.z );
119  firstAtom = false;
120  }
121  else
122  {
123  if ( static_cast<proshade_single> ( atom.pos.x ) > *xTo ) { *xTo = static_cast<proshade_single> ( atom.pos.x ); }
124  if ( static_cast<proshade_single> ( atom.pos.x ) < *xFrom ) { *xFrom = static_cast<proshade_single> ( atom.pos.x ); }
125  if ( static_cast<proshade_single> ( atom.pos.y ) > *yTo ) { *yTo = static_cast<proshade_single> ( atom.pos.y ); }
126  if ( static_cast<proshade_single> ( atom.pos.y ) < *yFrom ) { *yFrom = static_cast<proshade_single> ( atom.pos.y ); }
127  if ( static_cast<proshade_single> ( atom.pos.z ) > *zTo ) { *zTo = static_cast<proshade_single> ( atom.pos.z ); }
128  if ( static_cast<proshade_single> ( atom.pos.z ) < *zFrom ) { *zFrom = static_cast<proshade_single> ( atom.pos.z ); }
129  }
130  }
131  }
132  }
133  }
134  }
135  else
136  {
137  std::stringstream hlpSS;
138  hlpSS << "Found 0 models in input file " << pdbFile.name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
139  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
140  }
141 
142  //================================================ Done
143  return ;
144 
145 }

◆ distributeSpaceToBoundaries()

void ProSHADE_internal_mapManip::distributeSpaceToBoundaries ( proshade_signed &  minBound,
proshade_signed &  maxBound,
proshade_signed  oldBoundRange,
proshade_signed  newBoundRange 
)

Function for adding space to boundaries within the map confines.

This function takes the current boudaries and attempts to add the extra space specified by the parameters 3 and 4 equally to both the top and bottom boundaries, making sure the bottom boundary does not go belowe 0 and the top boundary does not exceed the final parameter.

Parameters
[in]minBoundReference to where the lower boundary value is held.
[in]maxBoundReference to where the upper boundary value is held.
[in]oldBoundRangeThe range between the original boundaries.
[in]newBoundRangeThe range which should be between the new boundaries.

Definition at line 2032 of file ProSHADE_mapManip.cpp.

2033 {
2034  //================================================ Only bother when sometings is to be added
2035  if ( newBoundRange > oldBoundRange )
2036  {
2037  //============================================ How much are we adding?
2038  proshade_signed distributeThis = newBoundRange - oldBoundRange;
2039 
2040  while ( distributeThis != 0 )
2041  {
2042  minBound -= 1;
2043  distributeThis -= 1;
2044 
2045  if ( distributeThis != 0 )
2046  {
2047  maxBound += 1;
2048  distributeThis -= 1;
2049  }
2050  }
2051  }
2052 
2053  //================================================ Done
2054  return ;
2055 
2056 }

◆ findMAPCOMValues()

void ProSHADE_internal_mapManip::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.

This function takes the proshade map object and procceds to compute the Centre of Mass of this object in Angstroms in real space.

Parameters
[in]mapA pointer to prshade density map.
[in]xComA pointer to proshade_double variable where the MAP file COM along the X-axis will be saved.
[in]yComA pointer to proshade_double variable where the MAP file COM along the Y-axis will be saved.
[in]zComA pointer to proshade_double variable where the MAP file COM along the Z-axis will be saved.
[in]xAngsThe map size in Angstroms along the X axis.
[in]yAngsThe map size in Angstroms along the Y axis.
[in]zAngsThe map size in Angstroms along the Z axis.
[in]xFromThe initial index of the x dimension of the map.
[in]xToThe terminal index of the x dimension of the map.
[in]yFromThe initial index of the y dimension of the map.
[in]yToThe terminal index of the y dimension of the map.
[in]zFromThe initial index of the z dimension of the map.
[in]zToThe terminal index of the z dimension of the map.

Definition at line 240 of file ProSHADE_mapManip.cpp.

241 {
242  //================================================ Initialise computation
243  proshade_double totDensity = 0.0;
244  *xCom = 0.0;
245  *yCom = 0.0;
246  *zCom = 0.0;
247  proshade_signed arrPos = 0;
248  proshade_single xSampRate = xAngs / static_cast< proshade_single > ( xTo - xFrom );
249  proshade_single ySampRate = yAngs / static_cast< proshade_single > ( yTo - yFrom );
250  proshade_single zSampRate = zAngs / static_cast< proshade_single > ( zTo - zFrom );
251 
252  //================================================ For each map point
253  for ( proshade_signed xIt = xFrom; xIt <= xTo; xIt++ )
254  {
255  for ( proshade_signed yIt = yFrom; yIt <= yTo; yIt++ )
256  {
257  for ( proshade_signed zIt = zFrom; zIt <= zTo; zIt++ )
258  {
259  arrPos = (zIt-zFrom) + ( zTo - zFrom + 1 ) * ( ( yIt - yFrom ) + ( yTo - yFrom + 1 ) * ( xIt - xFrom ) );
260 
261  if ( map[arrPos] > 0.0 )
262  {
263  totDensity += map[arrPos];
264  *xCom += static_cast<proshade_double> ( static_cast< proshade_single > ( xIt ) * xSampRate ) * map[arrPos];
265  *yCom += static_cast<proshade_double> ( static_cast< proshade_single > ( yIt ) * ySampRate ) * map[arrPos];
266  *zCom += static_cast<proshade_double> ( static_cast< proshade_single > ( zIt ) * zSampRate ) * map[arrPos];
267  }
268  }
269  }
270  }
271 
272  //================================================ Normalise sums to COM
273  *xCom /= totDensity;
274  *yCom /= totDensity;
275  *zCom /= totDensity;
276 
277  //================================================ Done
278  return ;
279 
280 }

◆ findPDBCOMValues()

void ProSHADE_internal_mapManip::findPDBCOMValues ( gemmi::Structure  pdbFile,
proshade_double *  xCom,
proshade_double *  yCom,
proshade_double *  zCom,
bool  firstModel 
)

This function finds the Centre of Mass for the co-ordinate file.

This function takes the gemmi::Structure object read from a co-odinate file and procceds to compute the Centre of Mass of this object.

Parameters
[in]pdbFileA gemmi::Structure object read in from input file.
[in]xComA pointer to proshade_double variable where the PDB file COM along the X-axis will be saved.
[in]yComA pointer to proshade_double variable where the PDB file COM along the Y-axis will be saved.
[in]zComA pointer to proshade_double variable where the PDB file COM along the Z-axis will be saved.
[in]firstModelShould only the first, or all models be used?

Definition at line 157 of file ProSHADE_mapManip.cpp.

158 {
159  //================================================ Initialise structure crawl
160  proshade_double totAtoms = 0.0;
161  *xCom = 0.0;
162  *yCom = 0.0;
163  *zCom = 0.0;
164 
165  //================================================ Use the first model, if it exists
166  if ( pdbFile.models.size() > 0 )
167  {
168  //============================================ For each model
169  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); sIt++ )
170  {
171  //======================================== Get model
172  gemmi::Model model = pdbFile.models.at(sIt);
173 
174  //======================================== Check if multiple models are allowed
175  if ( firstModel && ( sIt != 0 ) ) { break; }
176 
177  //======================================== For each chain
178  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model.chains.size() ); mIt++ )
179  {
180  //==================================== Get chain
181  gemmi::Chain chain = model.chains.at(mIt);
182 
183  //==================================== For each residue
184  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain.residues.size() ); rIt++ )
185  {
186  //================================ Get residue
187  gemmi::Residue residue = chain.residues.at(rIt);
188 
189  //================================ For each atom
190  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue.atoms.size() ); aIt++ )
191  {
192  //============================ Get atom
193  gemmi::Atom atom = residue.atoms.at(aIt);
194 
195  //============================ Save the COM sums
196  *xCom += atom.pos.x * atom.element.weight();
197  *yCom += atom.pos.y * atom.element.weight();
198  *zCom += atom.pos.z * atom.element.weight();
199  totAtoms += atom.element.weight();
200  }
201  }
202  }
203  }
204  }
205  else
206  {
207  std::stringstream hlpSS;
208  hlpSS << "Found 0 models in input file " << pdbFile.name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
209  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
210  }
211 
212  //================================================ Normalise sums to COM
213  *xCom /= totAtoms;
214  *yCom /= totAtoms;
215  *zCom /= totAtoms;
216 
217  //================================================ Done
218  return ;
219 
220 }

◆ generateMapFromPDB()

void ProSHADE_internal_mapManip::generateMapFromPDB ( gemmi::Structure  pdbFile,
proshade_double *&  map,
proshade_single  requestedResolution,
proshade_single  xCell,
proshade_single  yCell,
proshade_single  zCell,
proshade_signed *  xTo,
proshade_signed *  yTo,
proshade_signed *  zTo,
bool  forceP1,
bool  firstModel 
)

This function generates a theoretical map from co-ordinate input files.

This function makes use of the Gemmi internal functionality to firstly create a properly sized cell object, then to check the input co-ordinate file for containing known elements as well as elements for which Gemmi knows the form factors. Then, this function proceeds to compute the F's using Cromer & Libermann method (from Gemmi) to finally compute the theoretical map using these. This map is then copied to the variable supplied in the second argument, presumably the ProSHADE internal map variable.

Parameters
[in]pdbFileA gemmi::Structure object read in from the input file.
[in]mapPointer reference to a variable to save the map data.
[in]requestedResolutionMap resolution to which the map should be computed.
[in]xCellThe size of the map cell to be created.
[in]yCellThe size of the map cell to be created.
[in]zCellThe size of the map cell to be created.
[in]xToPointer to variable where the map size along the x-axis in indices will be saved.
[in]yToPointer to variable where the map size along the y-axis in indices will be saved.
[in]zToPointer to variable where the map size along the z-axis in indices will be saved.
[in]forceP1Should the P1 spacegroup be forced?
[in]firstModelShould only the first, or all models be used?
Warning
By default, this function will force the P1 spacegroup!

Definition at line 644 of file ProSHADE_mapManip.cpp.

645 {
646  //================================================ Set cell dimensions from the increased ranges (we need to add some space) and re-calculate cell properties
647  if ( forceP1 ) { pdbFile.cell = gemmi::UnitCell(); }
648  pdbFile.cell.a = static_cast< proshade_double > ( xCell );
649  pdbFile.cell.b = static_cast< proshade_double > ( yCell );
650  pdbFile.cell.c = static_cast< proshade_double > ( zCell );
651  pdbFile.cell.calculate_properties ( );
652 
653  //================================================ Get elements in Gemmi format
654  std::string totElString;
655  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); mIt++ )
656  {
657  //============================================ Check if multiple models are allowed
658  if ( firstModel && ( mIt != 0 ) )
659  {
660  std::stringstream hlpSS;
661  hlpSS << "!!! ProSHADE WARNING !!! Found multiple models (" << pdbFile.models.size() << ") in input file " << pdbFile.name << ", while the settings state that only the first PDB file model should be used. If all models should be used, please supply ProSHADE with the \"-x\" option.";
662  ProSHADE_internal_messages::printWarningMessage ( 0, hlpSS.str(), "WP00055" );
663  break;
664  }
665 
666  std::string hlpStr = pdbFile.models[mIt].present_elements ( ).to_string<char,std::char_traits<char>,std::allocator<char> >();
667  totElString = totElString + hlpStr;
668  }
669  std::bitset< static_cast< size_t > ( gemmi::El::END )> present_elems ( totElString );
670 
671  //================================================ Sanity checks
672  if ( present_elems[static_cast<int> ( gemmi::El::X )] )
673  {
674  throw ProSHADE_exception ( "Found unknown element in input file.", "EP00051", __FILE__, __LINE__, __func__, "Gemmi library does not recognise some of the elements in\n : the co-ordinate file. Please check the file for not being\n : corrupted and containing standard elements." );
675  }
676 
677  for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( present_elems.size() ); elIt++ )
678  {
679  if ( present_elems[elIt] && !gemmi::IT92<double>::has ( static_cast<gemmi::El> ( elIt ) ) )
680  {
681  std::stringstream hlpSS;
682  hlpSS << "Missing form factor for element " << element_name ( static_cast<gemmi::El> ( elIt ) );
683  throw ProSHADE_exception ( hlpSS.str().c_str(), "EP00052", __FILE__, __LINE__, __func__, "Gemmi library does not have a form factor value for this\n : reported element. Please report this to the author." );
684  }
685  }
686 
687  //================================================ Compute the f's
688  double wavelength = 0.1;
689  double energy = gemmi::hc() / wavelength;
690 
691  //================================================ Create the density calculator object and fill it in
692  gemmi::DensityCalculator<gemmi::IT92<double>, float> dencalc;
693 
694  dencalc.d_min = static_cast< double > ( requestedResolution );
695  for ( size_t elIt = 0; elIt < present_elems.size(); elIt++ ) { if ( present_elems[elIt] ) { dencalc.addends.set ( static_cast< gemmi::El > ( elIt ), static_cast< float > ( gemmi::cromer_liberman ( static_cast< int > ( elIt ), energy, nullptr ) ) ); } }
696  dencalc.set_grid_cell_and_spacegroup ( pdbFile );
697 
698  //================================================ Force P1 spacegroup
699  if ( forceP1 ) { dencalc.grid.spacegroup = &gemmi::get_spacegroup_p1(); }
700 
701  //================================================ Compute the theoretical map for each model
702  dencalc.grid.data.clear ( );
703  dencalc.grid.set_size_from_spacing ( dencalc.d_min / ( 2.0 * dencalc.rate), true );
704  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); mIt++ )
705  {
706  if ( firstModel && ( mIt != 0 ) ) { break; }
707  dencalc.add_model_density_to_grid ( pdbFile.models[mIt] );
708  dencalc.grid.symmetrize ( [](float a, float b) { return a + b; } );
709  }
710 
711  //================================================ Get the map
712  const gemmi::Grid<float>& grid = dencalc.grid;
713 
714  //================================================ Save the map dimensions
715  *xTo = grid.nu;
716  *yTo = grid.nv;
717  *zTo = grid.nw;
718 
719  //================================================ Copy the gemmi::Grid to my map format
720  map = new proshade_double [(*xTo) * (*yTo) * (*zTo)];
721  ProSHADE_internal_misc::checkMemoryAllocation ( map, __FILE__, __LINE__, __func__ );
722 
723  proshade_signed arrPos = 0;
724  for ( proshade_signed uIt = 0; uIt < (*xTo); uIt++ )
725  {
726  for ( proshade_signed vIt = 0; vIt < (*yTo); vIt++ )
727  {
728  for ( proshade_signed wIt = 0; wIt < (*zTo); wIt++ )
729  {
730  arrPos = wIt + (*zTo) * ( vIt + (*yTo) * uIt );
731  map[arrPos] = static_cast< proshade_double > ( grid.get_value_q( static_cast< int > ( uIt ), static_cast< int > ( vIt ), static_cast< int > ( wIt ) ) );
732  }
733  }
734  }
735 
736  //================================================ Done
737  return ;
738 
739 }

◆ getCorrelationMapMask()

void ProSHADE_internal_mapManip::getCorrelationMapMask ( proshade_double *&  map,
proshade_double *&  fakeHalfMap,
proshade_double *&  correlationMask,
proshade_unsign  xDimS,
proshade_unsign  yDimS,
proshade_unsign  zDimS,
proshade_signed  corrMaskKernel 
)

Function for creating the correlation mask.

This function is currently not used and should probably be deleted. The proper implementation of this masking approach should be available in the EMDA software.

Parameters
[in]mapA Reference Pointer to the map which should be blurred/sharpened.
[in]blurredMapA Reference Pointer to the variable which stores the fake half-map.
[in]correlationMaskA Reference Pointer to empty map where the mask will be saved.
[in]xDimSThe number of indices along the x axis of the map.
[in]yDimSThe number of indices along the y axis of the map.
[in]zDimSThe number of indices along the z axis of the map.
[in]corrMaskKernelThe amount of neighbours in any direction whose correlation is to be used to get the current points correlation.

Definition at line 1745 of file ProSHADE_mapManip.cpp.

1746 {
1747  //================================================ Set local variables
1748  proshade_signed xDim = static_cast< proshade_signed > ( xDimS ), yDim = static_cast< proshade_signed > ( yDimS ), zDim = static_cast< proshade_signed > ( zDimS ), currentPos, neighArrPos, neighXPos, neighYPos, neighZPos, corrIter;
1749  proshade_unsign noCorrVals = static_cast<proshade_unsign> ( pow ( ( ( corrMaskKernel * 2 ) + 1 ), 3 ) );
1750 
1751  //================================================ Alocate memory
1752  proshade_double *origMap = new proshade_double [noCorrVals];
1753  proshade_double *fakeHM = new proshade_double [noCorrVals];
1754 
1755  //================================================ Check memory allocation
1756  ProSHADE_internal_misc::checkMemoryAllocation ( origMap, __FILE__, __LINE__, __func__ );
1757  ProSHADE_internal_misc::checkMemoryAllocation ( fakeHM, __FILE__, __LINE__, __func__ );
1758 
1759  //================================================ Blur the coeffs
1760  for ( proshade_signed uIt = 0; uIt < xDim; uIt++ )
1761  {
1762  for ( proshade_signed vIt = 0; vIt < yDim; vIt++ )
1763  {
1764  for ( proshade_signed wIt = 0; wIt < zDim; wIt++ )
1765  {
1766  //==================================== Var init
1767  currentPos = wIt + zDim * ( vIt + yDim * uIt );
1768  corrIter = 0;
1769 
1770  //==================================== Average neighbours
1771  for ( proshade_signed xCh = -corrMaskKernel; xCh <= +corrMaskKernel; xCh++ )
1772  {
1773  for ( proshade_signed yCh = -corrMaskKernel; yCh <= +corrMaskKernel; yCh++ )
1774  {
1775  for ( proshade_signed zCh = -corrMaskKernel; zCh <= +corrMaskKernel; zCh++ )
1776  {
1777  //======================== Find the nieghbour peak indices (with periodicity)
1778  neighXPos = uIt + xCh; if ( neighXPos >= xDim ) { neighXPos -= xDim; }; if ( neighXPos < 0 ) { neighXPos += xDim; }
1779  neighYPos = vIt + yCh; if ( neighYPos >= yDim ) { neighYPos -= yDim; }; if ( neighYPos < 0 ) { neighYPos += yDim; }
1780  neighZPos = wIt + zCh; if ( neighZPos >= zDim ) { neighZPos -= zDim; }; if ( neighZPos < 0 ) { neighZPos += zDim; }
1781  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1782 
1783  //======================== Add to correct arrays
1784  origMap[corrIter] = map[neighArrPos];
1785  fakeHM[corrIter] = fakeHalfMap[neighArrPos];
1786  corrIter += 1;
1787  }
1788  }
1789  }
1790 
1791  //==================================== Save the correlation comparison result
1792  correlationMask[currentPos] = ProSHADE_internal_maths::pearsonCorrCoeff ( origMap, fakeHM, noCorrVals );
1793  }
1794  }
1795  }
1796 
1797  //================================================ Done
1798  return ;
1799 
1800 }

◆ getFakeHalfMap()

void ProSHADE_internal_mapManip::getFakeHalfMap ( proshade_double *&  map,
proshade_double *&  fakeHalfMap,
proshade_unsign  xDimS,
proshade_unsign  yDimS,
proshade_unsign  zDimS,
proshade_signed  fakeMapKernel 
)

Function for creating "fake" half-maps.

This function takes the internal map and an empty map array and proceeds to find all neighbours within the kernel distance to each map points. It then computes the average of these neighbours and saves this average as the value of the "fake" half-map. These are then useul for map masking, as the correlation between the "fake" half-maps and the original maps give nice masks according to Rangana

  • see him about how well these work.
Parameters
[in]mapA Reference Pointer to the map which should be blurred/sharpened.
[in]blurredMapA Reference Pointer to the variable which will store the modified map.
[in]xDimSThe number of indices along the x axis of the map.
[in]yDimSThe number of indices along the y axis of the map.
[in]zDimSThe number of indices along the z axis of the map.
[in]fakeMapKernelThe amount of neighbours in any direction whose average is to be used to get the current point.

Definition at line 1680 of file ProSHADE_mapManip.cpp.

1681 {
1682  //================================================ Set local variables
1683  proshade_signed xDim = static_cast< proshade_signed > ( xDimS );
1684  proshade_signed yDim = static_cast< proshade_signed > ( yDimS );
1685  proshade_signed zDim = static_cast< proshade_signed > ( zDimS );
1686  proshade_signed currentPos, neighArrPos, neighXPos, neighYPos, neighZPos;
1687  proshade_double neighSum;
1688  proshade_double neighCount = pow ( ( ( fakeMapKernel * 2 ) + 1 ), 3.0 ) - 1.0;
1689 
1690  //================================================ Blur the coeffs
1691  for ( proshade_signed uIt = 0; uIt < xDim; uIt++ )
1692  {
1693  for ( proshade_signed vIt = 0; vIt < yDim; vIt++ )
1694  {
1695  for ( proshade_signed wIt = 0; wIt < zDim; wIt++ )
1696  {
1697  //==================================== Var init
1698  currentPos = wIt + zDim * ( vIt + yDim * uIt );
1699  neighSum = 0.0;
1700 
1701  //==================================== Average neighbours
1702  for ( proshade_signed xCh = -fakeMapKernel; xCh <= +fakeMapKernel; xCh++ )
1703  {
1704  for ( proshade_signed yCh = -fakeMapKernel; yCh <= +fakeMapKernel; yCh++ )
1705  {
1706  for ( proshade_signed zCh = -fakeMapKernel; zCh <= +fakeMapKernel; zCh++ )
1707  {
1708  if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) { continue; }
1709 
1710  //======================== Find the nieghbour peak indices (with periodicity)
1711  neighXPos = uIt + xCh; if ( neighXPos >= xDim ) { neighXPos -= xDim; }; if ( neighXPos < 0 ) { neighXPos += xDim; }
1712  neighYPos = vIt + yCh; if ( neighYPos >= yDim ) { neighYPos -= yDim; }; if ( neighYPos < 0 ) { neighYPos += yDim; }
1713  neighZPos = wIt + zCh; if ( neighZPos >= zDim ) { neighZPos -= zDim; }; if ( neighZPos < 0 ) { neighZPos += zDim; }
1714  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1715 
1716  //======================== Add to average
1717  neighSum += map[neighArrPos];
1718  }
1719  }
1720  }
1721 
1722  //==================================== Save the average to "fake" half-map
1723  fakeHalfMap[currentPos] = neighSum / neighCount;
1724  }
1725  }
1726  }
1727 
1728  //================================================ Done
1729  return ;
1730 
1731 }

◆ getIndicesFromAngstroms()

proshade_single ProSHADE_internal_mapManip::getIndicesFromAngstroms ( proshade_unsign  xDim,
proshade_unsign  yDim,
proshade_unsign  zDim,
proshade_single  xAngs,
proshade_single  yAngs,
proshade_single  zAngs,
proshade_single  dist 
)

This function converts distance in Angstroms to distance in map indices.

This function finds out how many indices are required to cover a space of size "dist" in Angstroms. If you need this rounded in any way (ceil, floor, ...), just apply appropriate function to the output of this function.

Parameters
[in]xDimThe number of map indices along the x-axis.
[in]yDimThe number of map indices along the y-axis.
[in]zDimThe number of map indices along the z-axis.
[in]xAngsThe map size in Angstroms along the x-axis.
[in]yAngsThe map size in Angstroms along the y-axis.
[in]zAngsThe map size in Angstroms along the z-axis.
[in]distThe distance in Angstroms to be converted

Definition at line 1816 of file ProSHADE_mapManip.cpp.

1817 {
1818  //================================================ Compute
1819  proshade_single ret = static_cast< proshade_single > ( ProSHADE_internal_mapManip::myRound ( std::max ( dist / ( xAngs / static_cast<proshade_single> ( xDim ) ),
1820  std::max ( dist / ( yAngs / static_cast<proshade_single> ( yDim ) ),
1821  dist / ( zAngs / static_cast<proshade_single> ( zDim ) ) ) ) ) );
1822 
1823  //================================================ Done
1824  return ( ret );
1825 
1826 }

◆ getMaskFromBlurr()

void ProSHADE_internal_mapManip::getMaskFromBlurr ( proshade_double *&  blurMap,
proshade_double *&  outMap,
proshade_unsign  xDim,
proshade_unsign  yDim,
proshade_unsign  zDim,
proshade_single  noIQRs 
)

Function for computing mask from blurred map.

This function takes a blurred map and proceeds to compute the cut-off threshold for masking. It then applies this mask to the blurred map and applies the masking to the second input map. The assumption is that this second map is the map before blurring, as non-zero mask points will not be changed. Careful about this!!! It, however, does not output the mask.

Parameters
[in]blurMapA Reference Pointer to the map which has been blurred for mask computation.
[in]outMapA Reference Pointer to the map which will be masked - this map should be the map before blurring.
[in]xDimThe number of indices along the x axis of the map.
[in]yDimThe number of indices along the y axis of the map.
[in]zDimThe number of indices along the z axis of the map.
[in]noIQRsThe number of inter-quartile ranges from the median which should be used to compute the threshold for masking.

Definition at line 1032 of file ProSHADE_mapManip.cpp.

1033 {
1034  //================================================ Initialise vector for map data
1035  std::vector<proshade_double> mapVals ( xDim * yDim * zDim, 0.0 );
1036 
1037  //================================================ Save map values in vector
1038  for ( proshade_unsign iter = 0; iter < ( xDim * yDim * zDim ); iter++ )
1039  {
1040  mapVals.at(iter) = blurMap[iter];
1041  }
1042 
1043  //================================================ Find median and IQRs
1044  proshade_double* medAndIQR = new proshade_double[2];
1045  ProSHADE_internal_maths::vectorMedianAndIQR ( &mapVals, medAndIQR );
1046 
1047  //================================================ Find the threshold
1048  proshade_double maskThreshold = medAndIQR[0] + ( medAndIQR[1] * static_cast<proshade_double> ( noIQRs ) );
1049 
1050  //================================================ Apply threshold
1051  for ( proshade_unsign iter = 0; iter < ( xDim * yDim * zDim ); iter++ )
1052  {
1053  if ( blurMap[iter] < maskThreshold )
1054  {
1055  outMap[iter] = 0.0;
1056  blurMap[iter] = 0.0;
1057  }
1058  }
1059 
1060  //================================================ Release vector values
1061  mapVals.clear ( );
1062 
1063  //================================================ Release memory
1064  delete[] medAndIQR;
1065 
1066  //================================================ Done
1067  return ;
1068 
1069 }

◆ getNonZeroBounds()

void ProSHADE_internal_mapManip::getNonZeroBounds ( proshade_double *  map,
proshade_signed  xDim,
proshade_signed  yDim,
proshade_signed  zDim,
proshade_signed *&  ret 
)

Function for finding the map boundaries enclosing positive only values.

This function takes a map and searches for minimum and maximum map indices in all three dimensions, which enclose all the non-zero map values.

Parameters
[in]mapA pointer to the map for which the bounds are to be found.
[in]xDimThe number of indices along the x axis of the map.
[in]yDimThe number of indices along the y axis of the map.
[in]zDimThe number of indices along the z axis of the map.
[in]retA proshade_unsign pointer to array of 6 to which the results will be saved (0 = minX; 1 = maxX; 2 = minY; 3 = maxY; 4 - minZ; 5 = maxZ).

Definition at line 1082 of file ProSHADE_mapManip.cpp.

1083 {
1084  //================================================ Initialise local variables
1085  proshade_signed arrayPos = 0;
1086 
1087  //================================================ Initialise result variable
1088  ret[0] = xDim;
1089  ret[1] = 0;
1090  ret[2] = yDim;
1091  ret[3] = 0;
1092  ret[4] = zDim;
1093  ret[5] = 0;
1094 
1095  //================================================ Iterate through map and check bounds
1096  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1097  {
1098  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1099  {
1100  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1101  {
1102  //==================================== Var init
1103  arrayPos = zIt + zDim * ( yIt + yDim * xIt );
1104 
1105  //==================================== Check bounds
1106  if ( map[arrayPos] > 0.001 )
1107  {
1108  if ( xIt < ret[0] ) { ret[0] = xIt; }
1109  if ( xIt > ret[1] ) { ret[1] = xIt; }
1110  if ( yIt < ret[2] ) { ret[2] = yIt; }
1111  if ( yIt > ret[3] ) { ret[3] = yIt; }
1112  if ( zIt < ret[4] ) { ret[4] = zIt; }
1113  if ( zIt > ret[5] ) { ret[5] = zIt; }
1114  }
1115  }
1116  }
1117  }
1118 
1119  //================================================ Done
1120  return ;
1121 
1122 }

◆ moveMapByFourier()

void ProSHADE_internal_mapManip::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.

This function translates the map by changing the phase information of the map Fourier transform and then doing the inverse Fourier. This allows for movement smaller than one index distance, but should not be used over long distances (could move out of boundaries and meet pariodicity problem such as map from different cell now being moved into this cell).

Parameters
[in]mapA Reference Pointer to the map which should be shifted.
[in]xMovThe NEGATIVE value by how many angstroms should the x axis be moved.
[in]yMovThe NEGATIVE value by how many angstroms should the y axis be moved.
[in]zMovThe NEGATIVE value by how many angstroms should the z axis be moved.
[in]xAngsHow many angstroms are there along the x dimension.
[in]yAngsHow many angstroms are there along the y dimensionzAngs
[in]zAngsHow many angstroms are there along the z dimension.
[in]xDimHow many indices are there along the x dimension.
[in]yDimHow many indices are there along the y dimension.
[in]zDimHow many indices are there along the z dimension.

Definition at line 811 of file ProSHADE_mapManip.cpp.

812 {
813  //================================================ Local variables initialisation
814  proshade_unsign arrayPos = 0;
815  proshade_signed h, k, l;
816  proshade_double real = 0.0;
817  proshade_double imag = 0.0;
818  proshade_double trCoeffReal, trCoeffImag;
819  proshade_double normFactor = static_cast< proshade_double > ( xDim * yDim * zDim );
820  proshade_double exponent = 0.0;
821  proshade_double hlpArrReal;
822  proshade_double hlpArrImag;
823 
824  //================================================ Create Fourier map variable
825  fftw_complex *fCoeffs = new fftw_complex [xDim * yDim * zDim];
826  fftw_complex *translatedMap = new fftw_complex [xDim * yDim * zDim];
827 
828  //================================================ Check memory allocation
829  ProSHADE_internal_misc::checkMemoryAllocation ( fCoeffs, __FILE__, __LINE__, __func__ );
830  ProSHADE_internal_misc::checkMemoryAllocation ( translatedMap, __FILE__, __LINE__, __func__ );
831 
832  //================================================ Create plans
833  fftw_plan planForwardFourier = fftw_plan_dft_3d ( static_cast< int > ( xDim ), static_cast< int > ( yDim ), static_cast< int > ( zDim ), translatedMap, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
834  fftw_plan planBackwardFourier = fftw_plan_dft_3d ( static_cast< int > ( xDim ), static_cast< int > ( yDim ), static_cast< int > ( zDim ), fCoeffs, translatedMap, FFTW_BACKWARD, FFTW_ESTIMATE );
835 
836  //================================================ Copy map to complex format
837  for ( proshade_unsign uIt = 0; uIt < static_cast< proshade_unsign > ( xDim ); uIt++ )
838  {
839  for ( proshade_unsign vIt = 0; vIt < static_cast< proshade_unsign > ( yDim ); vIt++ )
840  {
841  for ( proshade_unsign wIt = 0; wIt < static_cast< proshade_unsign > ( zDim ); wIt++ )
842  {
843  arrayPos = wIt + static_cast< proshade_unsign > ( zDim ) * ( vIt + static_cast< proshade_unsign > ( yDim ) * uIt );
844 
845  const FloatingPoint< proshade_double > lhs ( map[arrayPos] ), rhs ( map[arrayPos] );
846  if ( lhs.AlmostEquals ( rhs ) ) { translatedMap[arrayPos][0] = map[arrayPos]; }
847  else { translatedMap[arrayPos][0] = 0.0; }
848  translatedMap[arrayPos][1] = 0.0;
849  }
850  }
851  }
852 
853  //================================================ Compute Forward Fourier
854  fftw_execute ( planForwardFourier );
855 
856  //================================================ Add the required shift
857  for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
858  {
859  for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
860  {
861  for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
862  {
863  //==================================== Var init
864  arrayPos = wIt + static_cast< proshade_unsign > ( zDim ) * ( vIt + static_cast< proshade_unsign > ( yDim ) * uIt );
865  real = fCoeffs[arrayPos][0];
866  imag = fCoeffs[arrayPos][1];
867 
868  //==================================== Change the B-factors, if required
869  if ( uIt > static_cast< proshade_unsign > ( (xDim+1) / 2) ) { h = static_cast < proshade_signed > ( uIt ) - xDim; } else { h = static_cast < proshade_signed > ( uIt ); }
870  if ( vIt > static_cast< proshade_unsign > ( (yDim+1) / 2) ) { k = static_cast < proshade_signed > ( vIt ) - yDim; } else { k = static_cast < proshade_signed > ( vIt ); }
871  if ( wIt > static_cast< proshade_unsign > ( (zDim+1) / 2) ) { l = static_cast < proshade_signed > ( wIt ) - zDim; } else { l = static_cast < proshade_signed > ( wIt ); }
872 
873  //==================================== Get translation coefficient change
874  exponent = ( ( ( static_cast <proshade_double> ( h ) / static_cast <proshade_double> ( xAngs ) ) * static_cast< proshade_double > ( -xMov ) ) +
875  ( ( static_cast <proshade_double> ( k ) / static_cast <proshade_double> ( yAngs ) ) * static_cast< proshade_double > ( -yMov ) ) +
876  ( ( static_cast <proshade_double> ( l ) / static_cast <proshade_double> ( zAngs ) ) * static_cast< proshade_double > ( -zMov ) ) ) * 2.0 * M_PI;
877 
878  trCoeffReal = cos ( exponent );
879  trCoeffImag = sin ( exponent );
880  ProSHADE_internal_maths::complexMultiplication ( &real, &imag, &trCoeffReal, &trCoeffImag, &hlpArrReal, &hlpArrImag );
881 
882  //==================================== Save the mask data
883  fCoeffs[arrayPos][0] = hlpArrReal / normFactor;
884  fCoeffs[arrayPos][1] = hlpArrImag / normFactor;
885  }
886  }
887  }
888 
889  //================================================ Compute inverse Fourier
890  fftw_execute ( planBackwardFourier );
891 
892  //================================================ Copy back to map
893  for ( proshade_unsign uIt = 0; uIt < static_cast< proshade_unsign > ( xDim ); uIt++ )
894  {
895  for ( proshade_unsign vIt = 0; vIt < static_cast< proshade_unsign > ( yDim ); vIt++ )
896  {
897  for ( proshade_unsign wIt = 0; wIt < static_cast< proshade_unsign > ( zDim ); wIt++ )
898  {
899  arrayPos = wIt + static_cast< proshade_unsign > ( zDim ) * ( vIt + static_cast< proshade_unsign > ( yDim ) * uIt );
900  map[arrayPos] = translatedMap[arrayPos][0];
901  }
902  }
903  }
904 
905  //================================================ Release memory
906  fftw_destroy_plan ( planForwardFourier );
907  fftw_destroy_plan ( planBackwardFourier );
908  delete[] fCoeffs;
909  delete[] translatedMap;
910 
911  //================================================ Done
912  return ;
913 
914 }

◆ moveMapByIndices()

void ProSHADE_internal_mapManip::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.

This function translates the map by changing the to and from index values so that the location of the map will be the same as the location of the original PDB file. This most likely cannot be done exactly as it allows only movement by increments of the sampling rate, but it is a decent start.

Parameters
[in]xMovPointer to the NEGATIVE value by how many angstroms should the x axis be moved.
[in]yMovPointer to the NEGATIVE value by how many angstroms should the y axis be moved.
[in]zMovPointer to the NEGATIVE value by how many angstroms should the z axis be moved.
[in]xAngsHow many angstroms are there along the x dimension.
[in]yAngsHow many angstroms are there along the y dimension.
[in]zAngsHow many angstroms are there along the z dimension.
[in]xFromThe initial index of the x dimension of the map.
[in]xToThe terminal index of the x dimension of the map.
[in]yFromThe initial index of the y dimension of the map.
[in]yToThe terminal index of the y dimension of the map.
[in]zFromThe initial index of the z dimension of the map.
[in]zToThe terminal index of the z dimension of the map.
[in]xOriginThe first value of the x axis index.
[in]yOriginThe first value of the y axis index.
[in]zOriginThe first value of the z axis index.

Definition at line 763 of file ProSHADE_mapManip.cpp.

764 {
765  //================================================ Compute movement in indices
766  proshade_single xIndMove = std::floor ( -(*xMov) / ( xAngs / ( static_cast< proshade_single > ( *xTo ) - static_cast< proshade_single > ( *xFrom ) ) ) );
767  proshade_single yIndMove = std::floor ( -(*yMov) / ( yAngs / ( static_cast< proshade_single > ( *yTo ) - static_cast< proshade_single > ( *yFrom ) ) ) );
768  proshade_single zIndMove = std::floor ( -(*zMov) / ( zAngs / ( static_cast< proshade_single > ( *zTo ) - static_cast< proshade_single > ( *zFrom ) ) ) );
769 
770  //================================================ Set the movs to the remainder
771  *xMov = -( *xMov ) - ( xIndMove * ( xAngs / ( static_cast< proshade_single > ( *xTo ) - static_cast< proshade_single > ( *xFrom ) ) ) );
772  *yMov = -( *yMov ) - ( yIndMove * ( yAngs / ( static_cast< proshade_single > ( *yTo ) - static_cast< proshade_single > ( *yFrom ) ) ) );
773  *zMov = -( *zMov ) - ( zIndMove * ( zAngs / ( static_cast< proshade_single > ( *zTo ) - static_cast< proshade_single > ( *zFrom ) ) ) );
774 
775  //================================================ Move indices by as much
776  *xFrom += static_cast< proshade_signed > ( xIndMove );
777  *xTo += static_cast< proshade_signed > ( xIndMove );
778  *yFrom += static_cast< proshade_signed > ( yIndMove );
779  *yTo += static_cast< proshade_signed > ( yIndMove );
780  *zFrom += static_cast< proshade_signed > ( zIndMove );
781  *zTo += static_cast< proshade_signed > ( zIndMove );
782 
783  //================================================ And set origin to reflect the changes
784  *xOrigin = *xFrom;
785  *yOrigin = *yFrom;
786  *zOrigin = *zFrom;
787 
788  //================================================ Done
789  return ;
790 
791 }

◆ movePDBForMapCalc()

void ProSHADE_internal_mapManip::movePDBForMapCalc ( gemmi::Structure *  pdbFile,
proshade_single  xMov,
proshade_single  yMov,
proshade_single  zMov,
bool  firstModel 
)

Function for moving co-ordinate atoms to better suit theoretical map computation.

This function translates all atoms by a given x, y and z distances. This is required as theoretical map computation can only output map cells starting from 0, 0, 0 and therefore to avoid density being in corners for PDB atoms not located in posivite axes, the atoms need to be moved. This effect should be reversed later.

Parameters
[in]pdbFileA pointer to the gemmi::Structure object as read in from the input file.
[in]xMovHow many angstroms should the atoms be moved along the x axis.
[in]yMovHow many angstroms should the atoms be moved along the y axis.
[in]zMovHow many angstroms should the atoms be moved along the z axis.
[in]firstModelShould only the first, or all models be used?

Definition at line 573 of file ProSHADE_mapManip.cpp.

574 {
575  //================================================ Use the first model, if it exists
576  if ( pdbFile->models.size() > 0 )
577  {
578  //============================================ For each model
579  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
580  {
581  //======================================== Check if multiple models are allowed
582  if ( firstModel && ( sIt != 0 ) ) { break; }
583 
584  //======================================== Get model
585  gemmi::Model *model = &pdbFile->models.at(sIt);
586 
587  //======================================== For each chain
588  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
589  {
590  //==================================== Get chain
591  gemmi::Chain *chain = &model->chains.at(mIt);
592 
593  //==================================== For each residue
594  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
595  {
596  //================================ Get residue
597  gemmi::Residue *residue = &chain->residues.at(rIt);
598 
599  //================================ For each atom
600  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
601  {
602  //============================ Get atom
603  gemmi::Atom *atom = &residue->atoms.at(aIt);
604 
605  //============================ Move the atoms
606  atom->pos = gemmi::Position ( atom->pos.x + static_cast< proshade_double > ( xMov ), atom->pos.y + static_cast< proshade_double > ( yMov ), atom->pos.z + static_cast< proshade_double > ( zMov ) );
607  }
608  }
609  }
610 
611  }
612  }
613  else
614  {
615  std::stringstream hlpSS;
616  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
617  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
618  }
619 
620  //================================================ Done
621  return ;
622 }

◆ myRound() [1/2]

proshade_signed ProSHADE_internal_mapManip::myRound ( proshade_double  x)

Calls the appropriate version of round function depending on compiler version.

Parameters
[in]xA decimal point number to be rounded.
[out]XThe rounded number.

Definition at line 31 of file ProSHADE_mapManip.cpp.

32 {
33 #if __cplusplus >= 201103L
34  return ( static_cast< proshade_signed > ( std::round ( x ) ) );
35 #else
36  return ( static_cast< proshade_signed > ( round ( x ) ) );
37 #endif
38 }

◆ myRound() [2/2]

proshade_signed ProSHADE_internal_mapManip::myRound ( proshade_single  x)

Calls the appropriate version of round function depending on compiler version.

Parameters
[in]xA decimal point number to be rounded.
[out]XThe rounded number.

Definition at line 46 of file ProSHADE_mapManip.cpp.

47 {
48 #if __cplusplus >= 201103L
49  return ( static_cast< proshade_signed > ( std::round ( x ) ) );
50 #else
51  return ( static_cast< proshade_signed > ( round ( x ) ) );
52 #endif
53 }

◆ releaseResolutionFourierMemory()

void ProSHADE_internal_mapManip::releaseResolutionFourierMemory ( fftw_complex *&  origMap,
fftw_complex *&  fCoeffs,
fftw_complex *&  newFCoeffs,
fftw_complex *&  newMap,
fftw_plan &  planForwardFourier,
fftw_plan &  planBackwardRescaledFourier 
)

This function releases the memory required by the Fourier resampling.

This function simply deletes all the memory allocated by the allocateResolutionFourierMemory() function.

Parameters
[in]origMapA Reference pointer to an array where the original map data were stored.
[in]fCoeffsA Reference pointer to an array where the original map Fourier coefficients data were stored.
[in]newFCoeffsA Reference pointer to an array where the re-sampled map Fourier coefficients data were stored.
[in]newMapA Reference pointer to an array where the re-sampled map data were stored.
[in]planForwardFourierFFTW_plan which computed the original map to Fourier coefficients transform.
[in]planBackwardRescaledFourierFFTW_plan which computed the re-sampled Fourier coefficients to re-sampled map transform.

Definition at line 1530 of file ProSHADE_mapManip.cpp.

1531 {
1532  //================================================ Delete the FFTW plans
1533  fftw_destroy_plan ( planForwardFourier );
1534  fftw_destroy_plan ( planBackwardRescaledFourier );
1535 
1536  //================================================ Delete the complex arrays
1537  delete[] origMap;
1538  delete[] fCoeffs;
1539  delete[] newFCoeffs;
1540  delete[] newMap;
1541 
1542  //================================================ Done
1543  return ;
1544 
1545 }

◆ removeMapPhase()

void ProSHADE_internal_mapManip::removeMapPhase ( fftw_complex *&  mapCoeffs,
proshade_unsign  xDim,
proshade_unsign  yDim,
proshade_unsign  zDim 
)

This function removes the phase from reciprocal (frequency) map.

This function takes an already FFTW-ed map and its dimensions as the input and proceeds to remove the phase from the map. It writes over the map and does not release any memory - it is the role of the calling function to deal with both these features.

Parameters
[in]mapCoeffsA Reference Pointer to the frequency map, from which phase is to be removed.
[in]xDimThe number of indices along the x-axis of the input map.
[in]yDimThe number of indices along the y-axis of the input map.
[in]zDimThe number of indices along the z-axis of the input map.

Definition at line 1631 of file ProSHADE_mapManip.cpp.

1632 {
1633  //================================================ Set local variables
1634  proshade_double real, imag, mag, phase;
1635  proshade_unsign arrayPos = 0;
1636  proshade_double normFactor = static_cast<proshade_double> ( xDim * yDim * zDim );
1637 
1638  //================================================ Iterate through the map
1639  for ( proshade_unsign uIt = 0; uIt < xDim; uIt++ )
1640  {
1641  for ( proshade_unsign vIt = 0; vIt < yDim; vIt++ )
1642  {
1643  for ( proshade_unsign wIt = 0; wIt < zDim; wIt++ )
1644  {
1645  //==================================== Var init
1646  arrayPos = wIt + zDim * ( vIt + yDim * uIt );
1647  real = mapCoeffs[arrayPos][0];
1648  imag = mapCoeffs[arrayPos][1];
1649 
1650  //==================================== Get magnitude and phase with mask parameters
1651  mag = std::sqrt ( (real*real) + (imag*imag) );;
1652  phase = 0.0; // This would be std::atan2 ( imag, real ); - but here we remove the phase.
1653 
1654  //==================================== Save the phaseless data
1655  mapCoeffs[arrayPos][0] = ( mag * cos(phase) ) / normFactor;
1656  mapCoeffs[arrayPos][1] = ( mag * sin(phase) ) / normFactor;
1657  }
1658  }
1659  }
1660 
1661  //================================================ Done
1662  return ;
1663 
1664 }

◆ removeWaters()

void ProSHADE_internal_mapManip::removeWaters ( gemmi::Structure *  pdbFile,
bool  firstModel 
)

This function removed all waters from PDB input file.

This function does a quick read-through the PDB file and deletes all water molecules from the PDB file. This is important as comparing two similar structures, one with hundreds of waters (with high occupancy in some cases) and one without will lead to the structures being different in ProSHADE's eyes.

Parameters
[in]pdbFileA pointer to gemmi::Structure object read in from the input file.
[in]firstModelShould only the first, or all models be used?

Definition at line 504 of file ProSHADE_mapManip.cpp.

505 {
506  //================================================ Use the first model, if it exists
507  if ( pdbFile->models.size() > 0 )
508  {
509  //============================================ For each model
510  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
511  {
512  //======================================== Check if multiple models are allowed
513  if ( firstModel && ( sIt != 0 ) ) { break; }
514 
515  //======================================== Get model
516  gemmi::Model *model = &pdbFile->models.at(sIt);
517 
518  //======================================== For each chain
519  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
520  {
521  //==================================== Get chain
522  gemmi::Chain *chain = &model->chains.at(mIt);
523 
524  //==================================== Initialise del vector
525  std::vector< proshade_unsign > delVec;
526 
527  //==================================== For each residue
528  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
529  {
530  //================================ Get residue
531  gemmi::Residue *residue = &chain->residues.at(rIt);
532 
533  //================================ If residue is water
534  if ( residue->is_water() )
535  {
537  }
538  }
539 
540  //==================================== Delete from end to avoid indexing issues
541  std::sort ( delVec.begin(), delVec.end(), std::greater<int>() );
542  for ( proshade_unsign vecIt = 0; vecIt < static_cast<proshade_unsign> ( delVec.size() ); vecIt++ )
543  {
544  chain->residues.erase ( chain->residues.begin() + static_cast< long int > ( delVec.at(vecIt) ) );
545  }
546  }
547  }
548  }
549  else
550  {
551  std::stringstream hlpSS;
552  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
553  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
554  }
555 
556  //================================================ Done
557  return ;
558 
559 }

◆ reSampleMapToResolutionFourier()

void ProSHADE_internal_mapManip::reSampleMapToResolutionFourier ( proshade_double *&  map,
proshade_single  resolution,
proshade_unsign  xDimS,
proshade_unsign  yDimS,
proshade_unsign  zDimS,
proshade_single  xAngs,
proshade_single  yAngs,
proshade_single  zAngs,
proshade_single *&  corrs 
)

This function re-samples a map to conform to given resolution using Fourier.

This function re-samples the internal map to a given resolutution by removing or zero-padding the Fourier (reciprocal space) coefficients and computing the inverse Fourier transform. This is the default option for map re-sampling, should it be required by the user.

Parameters
[in]mapA Reference Pointer to the map for which the bounds are to be found.
[in]resolutionThe required resolution value.
[in]xDimSThe number of indices along the x axis of the map.
[in]yDimSThe number of indices along the y axis of the map.
[in]zDimSThe number of indices along the z axis of the map.
[in]xAngsThe size of the x dimension of the map in angstroms.
[in]yAngsThe size of the y dimension of the map in angstroms.
[in]zAngsThe size of the z dimension of the map in angstroms.
[in]corrsPointer reference to proshade_single array of 6 values with the following meaning: 0 = xAdd; 1 = yAdd; 2 = zAdd; 3 = newXAng; 4 = newYAng; 5 = newZAng

Definition at line 1376 of file ProSHADE_mapManip.cpp.

1377 {
1378  //================================================ Sanity check - the resolution needs to be set
1379  if ( resolution <= 0.0f )
1380  {
1381  throw ProSHADE_exception ( "Requested resolution not set for map re-sampling.", "EM00015", __FILE__, __LINE__, __func__, "There is no resolution value set, but map re-sampling to\n : this unset resolution value is required. This error\n : occurs when a task with no resolution requirement is\n : requested on a map data and the map resolution change is\n : set to \'on\'. Either supply a resolution value, or do not\n : re-sample the map." );
1382  }
1383 
1384  //================================================ Initialise variables
1385  proshade_unsign newXDim = static_cast<proshade_unsign> ( ProSHADE_internal_mapManip::myRound ( xAngs / ( resolution / 2.0f ) ) );
1386  proshade_unsign newYDim = static_cast<proshade_unsign> ( ProSHADE_internal_mapManip::myRound ( yAngs / ( resolution / 2.0f ) ) );
1387  proshade_unsign newZDim = static_cast<proshade_unsign> ( ProSHADE_internal_mapManip::myRound ( zAngs / ( resolution / 2.0f ) ) );
1388 
1389  if ( newXDim % 2 != 0 ) { newXDim += 1; }
1390  if ( newYDim % 2 != 0 ) { newYDim += 1; }
1391  if ( newZDim % 2 != 0 ) { newZDim += 1; }
1392 
1393  proshade_signed preXChange, preYChange, preZChange;
1394  if ( ( xDimS % 2 ) == 0 ) { preXChange = static_cast< proshade_signed > ( std::ceil ( ( static_cast<proshade_signed> ( xDimS ) - static_cast<proshade_signed> ( newXDim ) ) / 2 ) ); }
1395  else { preXChange = static_cast< proshade_signed > ( std::floor ( ( static_cast<proshade_signed> ( xDimS ) - static_cast<proshade_signed> ( newXDim ) ) / 2 ) ); }
1396  if ( ( yDimS % 2 ) == 0 ) { preYChange = static_cast< proshade_signed > ( std::ceil ( ( static_cast<proshade_signed> ( yDimS ) - static_cast<proshade_signed> ( newYDim ) ) / 2 ) ); }
1397  else { preYChange = static_cast< proshade_signed > ( std::floor ( ( static_cast<proshade_signed> ( yDimS ) - static_cast<proshade_signed> ( newYDim ) ) / 2 ) ); }
1398  if ( ( zDimS % 2 ) == 0 ) { preZChange = static_cast< proshade_signed > ( std::ceil ( ( static_cast<proshade_signed> ( zDimS ) - static_cast<proshade_signed> ( newZDim ) ) / 2 ) ); }
1399  else { preZChange = static_cast< proshade_signed > ( std::floor ( ( static_cast<proshade_signed> ( zDimS ) - static_cast<proshade_signed> ( newZDim ) ) / 2 ) ); }
1400 
1401  proshade_signed postXChange = static_cast<proshade_signed> ( xDimS ) - ( preXChange + static_cast<proshade_signed> ( newXDim ) );
1402  proshade_signed postYChange = static_cast<proshade_signed> ( yDimS ) - ( preYChange + static_cast<proshade_signed> ( newYDim ) );
1403  proshade_signed postZChange = static_cast<proshade_signed> ( zDimS ) - ( preZChange + static_cast<proshade_signed> ( newZDim ) );
1404 
1405  proshade_unsign origSizeArr = 0, newSizeArr = 0;
1406  proshade_double normFactor = static_cast<proshade_double> ( xDimS * yDimS * zDimS );
1407 
1408  //================================================ Manage memory
1409  fftw_complex *origMap, *fCoeffs, *newFCoeffs, *newMap;
1410  fftw_plan planForwardFourier, planBackwardRescaledFourier;
1411  allocateResolutionFourierMemory ( origMap, fCoeffs, newFCoeffs, newMap, planForwardFourier, planBackwardRescaledFourier,
1412  xDimS, yDimS, zDimS, newXDim, newYDim, newZDim );
1413 
1414  //================================================ Fill maps with data and zeroes
1415  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDimS * yDimS * zDimS ); iter++ ) { origMap[iter][0] = map[iter]; origMap[iter][1] = 0.0; }
1416  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ ) { newFCoeffs[iter][0] = 0.0; newFCoeffs[iter][1] = 0.0; }
1417 
1418  //================================================ Get the Fourier coeffs
1419  fftw_execute ( planForwardFourier );
1420 
1421  //================================================ Change the order of Fourier coefficients
1422  changeFourierOrder ( fCoeffs, static_cast< proshade_signed > ( xDimS ), static_cast< proshade_signed > ( yDimS ), static_cast< proshade_signed > ( zDimS ), true );
1423 
1424  //================================================ Re-sample the coefficients by removing high frequencies or adding these with 0 values
1425  for ( proshade_unsign xIt = 0; xIt < newXDim; xIt++ )
1426  {
1427  for ( proshade_unsign yIt = 0; yIt < newYDim; yIt++ )
1428  {
1429  for ( proshade_unsign zIt = 0; zIt < newZDim; zIt++ )
1430  {
1431  //==================================== Find the array positions
1432  origSizeArr = ( ( zIt + static_cast< proshade_unsign > ( preZChange ) ) + zDimS *
1433  ( ( yIt + static_cast< proshade_unsign > ( preYChange ) ) + yDimS *
1434  ( xIt + static_cast< proshade_unsign > ( preXChange ) ) ) );
1435  newSizeArr = zIt + newZDim * ( yIt + newYDim * xIt );
1436 
1437  //==================================== If original coefficient for this new coefficient position exists, copy
1438  if ( ( ( -1 < static_cast< proshade_signed > ( xIt ) + preXChange ) && ( -1 < static_cast<proshade_signed> ( yIt ) + preYChange ) && ( -1 < static_cast<proshade_signed> ( zIt ) + preZChange ) ) &&
1439  ( ( xIt < newXDim + static_cast<proshade_unsign> ( postXChange ) ) && ( yIt < newYDim + static_cast<proshade_unsign> ( postYChange ) ) && ( zIt < newZDim + static_cast<proshade_unsign> ( postZChange ) ) ) )
1440  {
1441  //================================ Copy the Fourier coeff
1442  newFCoeffs[newSizeArr][0] = fCoeffs[origSizeArr][0] / normFactor;
1443  newFCoeffs[newSizeArr][1] = fCoeffs[origSizeArr][1] / normFactor;
1444  }
1445  }
1446  }
1447  }
1448 
1449  //================================================ Change the order of the re-sampled Fourier coefficients
1450  changeFourierOrder ( newFCoeffs, static_cast< proshade_signed > ( newXDim ), static_cast< proshade_signed > ( newYDim ), static_cast< proshade_signed > ( newZDim ), false );
1451 
1452  //================================================ Get the new map from the re-sized Fourier coefficients
1453  fftw_execute ( planBackwardRescaledFourier );
1454 
1455  //================================================ Delete the old map and create a new, re-sized one. Then copy the new map values into this new map memory.
1456  delete map;
1457  map = new proshade_double [newXDim * newYDim * newZDim];
1458  ProSHADE_internal_misc::checkMemoryAllocation ( map, __FILE__, __LINE__, __func__ );
1459  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ ) { map[iter] = newMap[iter][0]; }
1460 
1461  //================================================ Release memory
1462  releaseResolutionFourierMemory ( origMap, fCoeffs, newFCoeffs, newMap, planForwardFourier, planBackwardRescaledFourier );
1463 
1464  //================================================ Define change in indices and return it
1465  corrs[0] = static_cast< proshade_single > ( newXDim ) - static_cast< proshade_single > ( xDimS );
1466  corrs[1] = static_cast< proshade_single > ( newYDim ) - static_cast< proshade_single > ( yDimS );
1467  corrs[2] = static_cast< proshade_single > ( newZDim ) - static_cast< proshade_single > ( zDimS );
1468  corrs[3] = static_cast< proshade_single > ( newXDim ) * static_cast< proshade_single > ( resolution / 2.0f );
1469  corrs[4] = static_cast< proshade_single > ( newYDim ) * static_cast< proshade_single > ( resolution / 2.0f );
1470  corrs[5] = static_cast< proshade_single > ( newZDim ) * static_cast< proshade_single > ( resolution / 2.0f );
1471 
1472 
1473  //======================================== Done
1474  return ;
1475 
1476 }

◆ reSampleMapToResolutionTrilinear()

void ProSHADE_internal_mapManip::reSampleMapToResolutionTrilinear ( proshade_double *&  map,
proshade_single  resolution,
proshade_unsign  xDimS,
proshade_unsign  yDimS,
proshade_unsign  zDimS,
proshade_single  xAngs,
proshade_single  yAngs,
proshade_single  zAngs,
proshade_single *&  corrs 
)

This function re-samples a map to conform to given resolution using tri-linear interpolation.

This function takes a map and resolution value and it proceeds to create a new map, which has sampling resolution/2 and is large enough to contain the original map. It then proceeds to interpolate the new map values from the old map values, re-writing the old map once interpolation is done.

Parameters
[in]mapA Reference Pointer to the map for which the bounds are to be found.
[in]resolutionThe required resolution value.
[in]xDimSThe number of indices along the x axis of the map.
[in]yDimSThe number of indices along the y axis of the map.
[in]zDimSThe number of indices along the z axis of the map.
[in]xAngsThe size of the x dimension of the map in angstroms.
[in]yAngsThe size of the y dimension of the map in angstroms.
[in]zAngsThe size of the z dimension of the map in angstroms.
[in]corrsPointer reference to proshade_single array of 6 values with the following meaning: 0 = xAdd; 1 = yAdd; 2 = zAdd; 3 = newXAng; 4 = newYAng; 5 = newZAng

Definition at line 1175 of file ProSHADE_mapManip.cpp.

1176 {
1177  //================================================ Sanity check - the resolution needs to be set
1178  if ( resolution <= 0.0f )
1179  {
1180  throw ProSHADE_exception ( "Requested resolution not set for map re-sampling.", "EM00015", __FILE__, __LINE__, __func__, "There is no resolution value set, but map re-sampling to\n : this unset resolution value is required. This error\n : occurs when a task with no resolution requirement is\n : requested on a map data and the map resolution change is\n : set to \'on\'. Either supply a resolution value, or do not\n : re-sample the map." );
1181  }
1182 
1183  //================================================ Initialise local variables
1184  proshade_signed xDim = static_cast<proshade_signed> ( xDimS );
1185  proshade_signed yDim = static_cast<proshade_signed> ( yDimS );
1186  proshade_signed zDim = static_cast<proshade_signed> ( zDimS );
1187  proshade_single oldXSample = ( xAngs / static_cast<proshade_single> ( xDim ) );
1188  proshade_single oldYSample = ( yAngs / static_cast<proshade_single> ( yDim ) );
1189  proshade_single oldZSample = ( zAngs / static_cast<proshade_single> ( zDim ) );
1190  proshade_single newXSample = static_cast< proshade_single > ( resolution / 2.0f );
1191  proshade_single newYSample = static_cast< proshade_single > ( resolution / 2.0f );
1192  proshade_single newZSample = static_cast< proshade_single > ( resolution / 2.0f );
1193 
1194  //================================================ Compute required grid size
1195  proshade_signed newXDim = static_cast<proshade_signed> ( std::ceil ( xAngs / newXSample ) );
1196  proshade_signed newYDim = static_cast<proshade_signed> ( std::ceil ( yAngs / newYSample ) );
1197  proshade_signed newZDim = static_cast<proshade_signed> ( std::ceil ( zAngs / newZSample ) );
1198 
1199  //================================================ Create a new map variable
1200  proshade_double* newMap = new proshade_double [newXDim * newYDim * newZDim];
1201 
1202  //================================================ For each new map point
1203  proshade_signed xBottom = 0, xTop, yBottom = 0, yTop, zBottom = 0, zTop, oldMapIndex, newMapIndex;
1204  std::vector<proshade_double> c000 = std::vector<proshade_double> ( 4, 0.0 );
1205  std::vector<proshade_double> c001 = std::vector<proshade_double> ( 4, 0.0 );
1206  std::vector<proshade_double> c010 = std::vector<proshade_double> ( 4, 0.0 );
1207  std::vector<proshade_double> c011 = std::vector<proshade_double> ( 4, 0.0 );
1208  std::vector<proshade_double> c100 = std::vector<proshade_double> ( 4, 0.0 );
1209  std::vector<proshade_double> c101 = std::vector<proshade_double> ( 4, 0.0 );
1210  std::vector<proshade_double> c110 = std::vector<proshade_double> ( 4, 0.0 );
1211  std::vector<proshade_double> c111 = std::vector<proshade_double> ( 4, 0.0 );
1212  std::vector<proshade_double> c00 = std::vector<proshade_double> ( 4, 0.0 );
1213  std::vector<proshade_double> c01 = std::vector<proshade_double> ( 4, 0.0 );
1214  std::vector<proshade_double> c10 = std::vector<proshade_double> ( 4, 0.0 );
1215  std::vector<proshade_double> c11 = std::vector<proshade_double> ( 4, 0.0 );
1216  std::vector<proshade_double> c0 = std::vector<proshade_double> ( 4, 0.0 );
1217  std::vector<proshade_double> c1 = std::vector<proshade_double> ( 4, 0.0 );
1218  proshade_double xRelative, yRelative, zRelative;
1219 
1220  for ( proshade_signed xIt = 0; xIt < newXDim; xIt++ )
1221  {
1222  for ( proshade_signed yIt = 0; yIt < newYDim; yIt++ )
1223  {
1224  for ( proshade_signed zIt = 0; zIt < newZDim; zIt++ )
1225  {
1226  //==================================== Get this point's index
1227  newMapIndex = zIt + newZDim * ( yIt + newYDim * xIt );
1228 
1229  //==================================== Find this points bottom and top positions in the old map (including periodicity)
1230  for ( proshade_signed ox = 0; ox < ( static_cast< proshade_signed > ( xDimS ) - 1 ); ox++ ) { if ( ( ( static_cast< proshade_single > ( xIt ) * newXSample ) >= ( static_cast< proshade_single > ( ox ) * oldXSample ) ) && ( ( static_cast< proshade_single > ( xIt ) * newXSample ) <= ( ( static_cast< proshade_single > ( ox ) + 1 ) * oldXSample ) ) ) { xBottom = ox; break; } }
1231  for ( proshade_signed oy = 0; oy < ( static_cast< proshade_signed > ( yDimS ) - 1 ); oy++ ) { if ( ( ( static_cast< proshade_single > ( yIt ) * newYSample ) >= ( static_cast< proshade_single > ( oy ) * oldYSample ) ) && ( ( static_cast< proshade_single > ( yIt ) * newYSample ) <= ( ( static_cast< proshade_single > ( oy ) + 1 ) * oldYSample ) ) ) { yBottom = oy; break; } }
1232  for ( proshade_signed oz = 0; oz < ( static_cast< proshade_signed > ( zDimS ) - 1 ); oz++ ) { if ( ( ( static_cast< proshade_single > ( zIt ) * newZSample ) >= ( static_cast< proshade_single > ( oz ) * oldZSample ) ) && ( ( static_cast< proshade_single > ( zIt ) * newZSample ) <= ( ( static_cast< proshade_single > ( oz ) + 1 ) * oldZSample ) ) ) { zBottom = oz; break; } }
1233  xTop = xBottom + 1;
1234  yTop = yBottom + 1;
1235  zTop = zBottom + 1;
1236 
1237  //==================================== Find the surrounding point's values from the original map
1238  oldMapIndex = zBottom + static_cast< proshade_signed > ( zDimS ) * ( yBottom + static_cast< proshade_signed > ( yDimS ) * xBottom );
1239  c000.at(0) = static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample );
1240  c000.at(1) = static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample );
1241  c000.at(2) = static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample );
1242  c000.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1243 
1244  oldMapIndex = zTop + static_cast< proshade_signed > ( zDimS ) * ( yBottom + static_cast< proshade_signed > ( yDimS ) * xBottom );
1245  c001.at(0) = static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample );
1246  c001.at(1) = static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample );
1247  c001.at(2) = static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample );
1248  c001.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1249 
1250  oldMapIndex = zBottom + static_cast< proshade_signed > ( zDimS ) * ( yTop + static_cast< proshade_signed > ( yDimS ) * xBottom );
1251  c010.at(0) = static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample );
1252  c010.at(1) = static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample );
1253  c010.at(2) = static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample );
1254  c010.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1255 
1256  oldMapIndex = zTop + static_cast< proshade_signed > ( zDimS ) * ( yTop + static_cast< proshade_signed > ( yDimS ) * xBottom );
1257  c011.at(0) = static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample );
1258  c011.at(1) = static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample );
1259  c011.at(2) = static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample );
1260  c011.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1261 
1262  oldMapIndex = zBottom + static_cast< proshade_signed > ( zDimS ) * ( yBottom + static_cast< proshade_signed > ( yDimS ) * xTop );
1263  c100.at(0) = static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample );
1264  c100.at(1) = static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample );
1265  c100.at(2) = static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample );
1266  c100.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1267 
1268  oldMapIndex = zTop + static_cast< proshade_signed > ( zDimS ) * ( yBottom + static_cast< proshade_signed > ( yDimS ) * xTop );
1269  c101.at(0) = static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample );
1270  c101.at(1) = static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample );
1271  c101.at(2) = static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample );
1272  c101.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1273 
1274  oldMapIndex = zBottom + static_cast< proshade_signed > ( zDimS ) * ( yTop + static_cast< proshade_signed > ( yDimS ) * xTop );
1275  c110.at(0) = static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample );
1276  c110.at(1) = static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample );
1277  c110.at(2) = static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample );
1278  c110.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1279 
1280  oldMapIndex = zTop + static_cast< proshade_signed > ( zDimS ) * ( yTop + static_cast< proshade_signed > ( yDimS ) * xTop );
1281  c111.at(0) = static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample );
1282  c111.at(1) = static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample );
1283  c111.at(2) = static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample );
1284  c111.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1285 
1286  //==================================== Interpolate to the new grid along X
1287  xRelative = ( ( static_cast<proshade_double> ( xIt ) * static_cast<proshade_double> ( newXSample ) ) - ( static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample ) ) ) / ( ( static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample ) ) - ( static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample ) ) );
1288 
1289  //==================================== Interpolate for the less less point
1290  c00.at(0) = ( static_cast< proshade_double > ( newXSample ) * xRelative ) + c000.at(0);
1291  c00.at(1) = c000.at(1);
1292  c00.at(2) = c000.at(2);
1293  c00.at(3) = ( c000.at(3) * ( 1.0 - xRelative ) ) + ( c100.at(3) * xRelative );
1294 
1295  //==================================== Interpolate for the less more point
1296  c01.at(0) = ( static_cast< proshade_double > ( newXSample ) * xRelative ) + c001.at(0);
1297  c01.at(1) = c001.at(1);
1298  c01.at(2) = c001.at(2);
1299  c01.at(3) = ( c001.at(3) * ( 1.0 - xRelative ) ) + ( c101.at(3) * xRelative );
1300 
1301  //==================================== Interpolate for the more less point
1302  c10.at(0) = ( static_cast< proshade_double > ( newXSample ) * xRelative ) + c010.at(0);
1303  c10.at(1) = c010.at(1);
1304  c10.at(2) = c010.at(2);
1305  c10.at(3) = ( c010.at(3) * ( 1.0 - xRelative ) ) + ( c110.at(3) * xRelative );
1306 
1307  //==================================== Interpolate for the more more point
1308  c11.at(0) = ( static_cast< proshade_double > ( newXSample ) * xRelative ) + c011.at(0);
1309  c11.at(1) = c011.at(1);
1310  c11.at(2) = c011.at(2);
1311  c11.at(3) = ( c011.at(3) * ( 1.0 - xRelative ) ) + ( c111.at(3) * xRelative );
1312 
1313  //==================================== Interpolate to the new grid along Y
1314  yRelative = ( ( static_cast<proshade_double> ( yIt ) * static_cast<proshade_double> ( newYSample ) ) - ( static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample ) ) ) / ( ( static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample ) ) - ( static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample ) ) );
1315 
1316  //==================================== Interpolate for the less point
1317  c0.at(0) = c00.at(0);
1318  c0.at(1) = ( static_cast< proshade_double > ( newYSample ) * yRelative ) + c00.at(1);
1319  c0.at(2) = c00.at(2);
1320  c0.at(3) = ( c00.at(3) * ( 1.0 - yRelative ) ) + ( c10.at(3) * yRelative );
1321 
1322  //==================================== Interpolate for the more point
1323  c1.at(0) = c01.at(0);
1324  c1.at(1) = ( static_cast< proshade_double > ( newYSample ) * yRelative ) + c01.at(1);
1325  c1.at(2) = c01.at(2);
1326  c1.at(3) = ( c01.at(3) * ( 1.0 - yRelative ) ) + ( c11.at(3) * yRelative );
1327 
1328  //==================================== Interpolate to the new grid along Z
1329  zRelative = ( ( static_cast<proshade_double> ( zIt ) * static_cast< proshade_double > ( newZSample ) ) - ( static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample ) ) ) / static_cast< proshade_double > ( ( static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample ) ) - ( static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample ) ) );
1330  newMap[newMapIndex] = ( c0.at(3) * ( 1.0 - zRelative ) ) + ( c1.at(3) * zRelative );
1331  }
1332  }
1333  }
1334 
1335  //================================================ Delete old map and allocate new memory
1336  delete[] map;
1337  map = new proshade_double [newXDim * newYDim * newZDim];
1338 
1339  //================================================ Copy map
1340  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ )
1341  {
1342  map[iter] = newMap[iter];
1343  }
1344 
1345  //================================================ Release memory
1346  delete[] newMap;
1347 
1348  //================================================ Define change in indices and return it
1349  corrs[0] = static_cast< proshade_single > ( newXDim - xDim );
1350  corrs[1] = static_cast< proshade_single > ( newYDim - yDim );
1351  corrs[2] = static_cast< proshade_single > ( newZDim - zDim );
1352  corrs[3] = static_cast< proshade_single > ( newXDim ) * static_cast< proshade_single > ( newXSample );
1353  corrs[4] = static_cast< proshade_single > ( newYDim ) * static_cast< proshade_single > ( newYSample );
1354  corrs[5] = static_cast< proshade_single > ( newZDim ) * static_cast< proshade_single > ( newZSample );
1355 
1356  //======================================== Done
1357  return ;
1358 
1359 }

◆ rotatePDBCoordinates()

void ProSHADE_internal_mapManip::rotatePDBCoordinates ( gemmi::Structure *  pdbFile,
proshade_double  euA,
proshade_double  euB,
proshade_double  euG,
proshade_double  xCom,
proshade_double  yCom,
proshade_double  zCom,
bool  firstModel 
)

Function for rotating the PDB file co-ordinates by Euler angles.

This function takes the three Euler angles and a pointer to a gemmi::Structure and it then proceeds to compute the rotation matrix from the Euler angles. This matrix is then applied to the co-ordinates in the gemmi::Structure in a way so that the rotation is done over the centre of the co-ordinates, but the co-ordinate positions stay unchanged.

Parameters
[in]pdbFilePointer to a gemmi::Structure object which will have its co-ordinates rotated.
[in]euAThe Euler angle alpha by which the co-ordinates should be rotated.
[in]euBThe Euler angle beta by which the co-ordinates should be rotated.
[in]euGThe Euler angle gamma by which the co-ordinates should be rotated.
[in]xComThe x-axis position around which the rotation should be applied.
[in]yComThe y-axis position around which the rotation should be applied.
[in]zComThe z-axis position around which the rotation should be applied.
[in]firstModelShould only the first, or all models be used?

Definition at line 296 of file ProSHADE_mapManip.cpp.

298 {
299  //================================================ Convert Euler angles to rotation matrix
300  proshade_double *rotMat = new proshade_double[9];
301  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
303 
304  //================================================ Initialise internal variables
305  proshade_double xTmp, yTmp, zTmp;
306 
307  //================================================ Use the first model, if it exists
308  if ( pdbFile->models.size() > 0 )
309  {
310  //============================================ For each model
311  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
312  {
313  //======================================== Get model
314  gemmi::Model *model = &pdbFile->models.at(sIt);
315 
316  //======================================== Check if multiple models are allowed
317  if ( firstModel && ( sIt != 0 ) ) { break; }
318 
319  //======================================== For each chain
320  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
321  {
322  //==================================== Get chain
323  gemmi::Chain *chain = &model->chains.at(mIt);
324 
325  //==================================== For each residue
326  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
327  {
328  //================================ Get residue
329  gemmi::Residue *residue = &chain->residues.at(rIt);
330 
331  //================================ For each atom
332  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
333  {
334  //============================ Get atom
335  gemmi::Atom *atom = &residue->atoms.at(aIt);
336 
337  //============================ Move to mid-point
338  xTmp = static_cast< proshade_double > ( atom->pos.x - xCom );
339  yTmp = static_cast< proshade_double > ( atom->pos.y - yCom );
340  zTmp = static_cast< proshade_double > ( atom->pos.z - zCom );
341 
342  //============================ Rotate the atom position
343  atom->pos.x = ( xTmp * rotMat[0] ) + ( yTmp * rotMat[1] ) + ( zTmp * rotMat[2] );
344  atom->pos.y = ( xTmp * rotMat[3] ) + ( yTmp * rotMat[4] ) + ( zTmp * rotMat[5] );
345  atom->pos.z = ( xTmp * rotMat[6] ) + ( yTmp * rotMat[7] ) + ( zTmp * rotMat[8] );
346 
347  //============================ Move back
348  atom->pos.x = atom->pos.x + xCom;
349  atom->pos.y = atom->pos.y + yCom;
350  atom->pos.z = atom->pos.z + zCom;
351  }
352  }
353  }
354  }
355  }
356  else
357  {
358  std::stringstream hlpSS;
359  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
360  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
361  }
362 
363  //================================================ Release memory
364  delete[] rotMat;
365 
366  //================================================ Done
367  return ;
368 
369 }

◆ translatePDBCoordinates()

void ProSHADE_internal_mapManip::translatePDBCoordinates ( gemmi::Structure *  pdbFile,
proshade_double  transX,
proshade_double  transY,
proshade_double  transZ,
bool  firstModel 
)

Function for translating the PDB file co-ordinates by given distances in Angstroms.

This function simply iterates through the given structure object and adds the required shift to all atoms in the first model along the three axes.

Parameters
[in]pdbFilePointer to a gemmi::Structure object whose co-ordinates are to be translated.
[in]transXThe
[in]transYThe
[in]transZThe
[in]firstModelShould only the first, or all models be used?

Definition at line 381 of file ProSHADE_mapManip.cpp.

382 {
383  //================================================ Use the first model, if it exists
384  if ( pdbFile->models.size() > 0 )
385  {
386  //============================================ For each model
387  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
388  {
389  //======================================== Check if multiple models are allowed
390  if ( firstModel && ( sIt != 0 ) ) { break; }
391 
392  //======================================== Get model
393  gemmi::Model *model = &pdbFile->models.at(sIt);
394 
395  //======================================== For each chain
396  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
397  {
398  //==================================== Get chain
399  gemmi::Chain *chain = &model->chains.at(mIt);
400 
401  //==================================== For each residue
402  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
403  {
404  //================================ Get residue
405  gemmi::Residue *residue = &chain->residues.at(rIt);
406 
407  //================================ For each atom
408  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
409  {
410  //============================ Get atom
411  gemmi::Atom *atom = &residue->atoms.at(aIt);
412 
413  //============================ Translate
414  atom->pos.x += transX;
415  atom->pos.y += transY;
416  atom->pos.z += transZ;
417  }
418  }
419  }
420  }
421  }
422  else
423  {
424  std::stringstream hlpSS;
425  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
426  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
427  }
428 
429  //================================================ Done
430  return ;
431 
432 }
ProSHADE_exception
This class is the representation of ProSHADE exception.
Definition: ProSHADE_exceptions.hpp:37
ProSHADE_internal_mapManip::distributeSpaceToBoundaries
void distributeSpaceToBoundaries(proshade_signed &minBound, proshade_signed &maxBound, proshade_signed oldBoundRange, proshade_signed newBoundRange)
Function for adding space to boundaries within the map confines.
Definition: ProSHADE_mapManip.cpp:2032
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_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_messages::printWarningMessage
void printWarningMessage(proshade_signed verbose, std::string message, std::string warnCode)
General stderr message printing (used for warnings).
Definition: ProSHADE_messages.cpp:101
ProSHADE_internal_mapManip::changeFourierOrder
void changeFourierOrder(fftw_complex *&fCoeffs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, bool negativeFirst)
This function changes the order of Fourier coefficients in a 3D array between positive first (default...
Definition: ProSHADE_mapManip.cpp:1559
ProSHADE_internal_mapManip::releaseResolutionFourierMemory
void releaseResolutionFourierMemory(fftw_complex *&origMap, fftw_complex *&fCoeffs, fftw_complex *&newFCoeffs, fftw_complex *&newMap, fftw_plan &planForwardFourier, fftw_plan &planBackwardRescaledFourier)
This function releases the memory required by the Fourier resampling.
Definition: ProSHADE_mapManip.cpp:1530
ProSHADE_internal_maths::pearsonCorrCoeff
proshade_double pearsonCorrCoeff(proshade_double *valSet1, proshade_double *valSet2, proshade_unsign length)
Function for computing the Pearson's correlation coefficient.
Definition: ProSHADE_maths.cpp:246
ProSHADE_internal_misc::addToSignedVector
void addToSignedVector(std::vector< proshade_signed > *vecToAddTo, proshade_signed elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:121
ProSHADE_internal_mapManip::getIndicesFromAngstroms
proshade_single getIndicesFromAngstroms(proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single dist)
This function converts distance in Angstroms to distance in map indices.
Definition: ProSHADE_mapManip.cpp:1816
ProSHADE_internal_mapManip::allocateResolutionFourierMemory
void allocateResolutionFourierMemory(fftw_complex *&origMap, fftw_complex *&fCoeffs, fftw_complex *&newFCoeffs, fftw_complex *&newMap, fftw_plan &planForwardFourier, fftw_plan &planBackwardRescaledFourier, proshade_unsign xDimOld, proshade_unsign yDimOld, proshade_unsign zDimOld, proshade_unsign xDimNew, proshade_unsign yDimNew, proshade_unsign zDimNew)
This function allocates and checks the allocatio of the memory required by the Fourier resampling.
Definition: ProSHADE_mapManip.cpp:1496
ProSHADE_internal_maths::primeFactorsDecomp
std::vector< proshade_signed > primeFactorsDecomp(proshade_signed number)
Function to find prime factors of an integer.
Definition: ProSHADE_maths.cpp:1707
ProSHADE_internal_mapManip::betterClosePrimeFactors
proshade_signed betterClosePrimeFactors(proshade_signed fromRange, proshade_signed toRange)
Function for finding close numbers with better prime factors.
Definition: ProSHADE_mapManip.cpp:1986
ProSHADE_internal_mapManip::myRound
proshade_signed myRound(proshade_double x)
Calls the appropriate version of round function depending on compiler version.
Definition: ProSHADE_mapManip.cpp:31
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_maths::vectorMedianAndIQR
void vectorMedianAndIQR(std::vector< proshade_double > *vec, proshade_double *&ret)
Function to get vector median and inter-quartile range.
Definition: ProSHADE_maths.cpp:149
ProSHADE_internal_misc::addToUnsignVector
void addToUnsignVector(std::vector< proshade_unsign > *vecToAddTo, proshade_unsign elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:99