Xmipp  v3.23.11-Nereus
Functions
lib_pwk.cpp File Reference
#include "situs.h"
#include "lib_pio.h"
#include "lib_pwk.h"
#include "lib_vwk.h"
#include "lib_vec.h"
#include "lib_err.h"
#include "lib_eul.h"
Include dependency graph for lib_pwk.cpp:

Go to the source code of this file.

Functions

void copy_atoms (PDB *pdb_original, PDB *pdb_duplicate, int io, int id, int ktot)
 
void rot_axis (PDB *pdb_original, PDB *pdb_rotate, unsigned num_atoms, char axis, double angle)
 
void rot_euler (PDB *pdb_original, PDB *pdb_rotate, unsigned num_atoms, double psi, double theta, double phi)
 
void translate (PDB *pdb_original, PDB *pdb_move, unsigned num_atoms, double x0, double y0, double z0)
 
void calc_center (PDB *pdb0, unsigned num_atoms, double *cx, double *cy, double *cz)
 
double calc_mass (PDB *pdb0, unsigned num_atoms)
 
double calc_center_mass (PDB *pdb0, unsigned num_atoms, double *cx, double *cy, double *cz)
 
double calc_sphere (PDB *pdb0, unsigned num_atoms, double cx, double cy, double cz)
 
void calc_box (PDB *pdb0, unsigned num_atoms, double *minx, double *miny, double *minz, double *maxx, double *maxy, double *maxz)
 
void project_mass (double **outmap, unsigned long nvox, double widthx, double widthy, double widthz, unsigned extx, unsigned exty, unsigned extz, PDB *inpdb, unsigned num_atoms, double shift[3], unsigned ignored[3])
 
void project_mass_convolve_kernel_corr (double widthx, double widthy, double widthz, unsigned extx, unsigned exty, unsigned extz, PDB *inpdb, unsigned num_atoms, double shift[3], double *kernel, unsigned kernel_size, double normfac, unsigned ignored[3], double *lowmap, double *corr_hi_low)
 
int check_if_inside (double fraction, double widthx, double widthy, double widthz, unsigned extx, unsigned exty, unsigned extz, PDB *inpdb, unsigned num_atoms, double shift[3])
 

Function Documentation

◆ calc_box()

void calc_box ( PDB pdb0,
unsigned  num_atoms,
double *  minx,
double *  miny,
double *  minz,
double *  maxx,
double *  maxy,
double *  maxz 
)

Definition at line 237 of file lib_pwk.cpp.

240 {
241  unsigned id;
242  const char *program = "lib_pwk";
243  const char *shape = "box";
244 
245  *minx = 1e20;
246  *miny = 1e20;
247  *minz = 1e20;
248  *maxx = -1e20;
249  *maxy = -1e20;
250  *maxz = -1e20;
251  for (id = 0; id < num_atoms; ++id) {
252  if (*minx > pdb0[id].x) *minx = pdb0[id].x;
253  if (*maxx < pdb0[id].x) *maxx = pdb0[id].x;
254  if (*miny > pdb0[id].y) *miny = pdb0[id].y;
255  if (*maxy < pdb0[id].y) *maxy = pdb0[id].y;
256  if (*minz > pdb0[id].z) *minz = pdb0[id].z;
257  if (*maxz < pdb0[id].z) *maxz = pdb0[id].z;
258  }
259  if (*minx < 1e20 && *miny < 1e20 && *minz < 1e20 && *maxx > -1e20 && *maxy > -1e20 && *maxz > -1e20) return;
260  else {
261  error_no_bounding(16040, program, shape);
262  }
263 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
static double * y
void error_no_bounding(int error_number, const char *program, const char *shape)
Definition: lib_err.cpp:258
doublereal * x
float x
Definition: situs.h:46
double z

◆ calc_center()

void calc_center ( PDB pdb0,
unsigned  num_atoms,
double *  cx,
double *  cy,
double *  cz 
)

Definition at line 140 of file lib_pwk.cpp.

142 {
143  unsigned id;
144  const char *program = "lib_pwk";
145  *cx = 0.0;
146  *cy = 0.0;
147  *cz = 0.0;
148  for (id = 0; id < num_atoms; ++id) {
149  *cx += pdb0[id].x;
150  *cy += pdb0[id].y;
151  *cz += pdb0[id].z;
152  }
153  if (num_atoms > 0) {
154  *cx /= (num_atoms * 1.0);
155  *cy /= (num_atoms * 1.0);
156  *cz /= (num_atoms * 1.0);
157  } else {
158  error_divide_zero(16010, program);
159  }
160 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
float x
Definition: situs.h:46
void error_divide_zero(int error_number, const char *program)
Definition: lib_err.cpp:217

◆ calc_center_mass()

double calc_center_mass ( PDB pdb0,
unsigned  num_atoms,
double *  cx,
double *  cy,
double *  cz 
)

Definition at line 179 of file lib_pwk.cpp.

181 {
182  unsigned id;
183  double mtot;
184  const char *program = "lib_pwk";
185 
186  *cx = 0.0;
187  *cy = 0.0;
188  *cz = 0.0;
189  mtot = 0.0;
190  for (id = 0; id < num_atoms; ++id) {
191  mtot += pdb0[id].weight;
192  *cx += pdb0[id].x * pdb0[id].weight;
193  *cy += pdb0[id].y * pdb0[id].weight;
194  *cz += pdb0[id].z * pdb0[id].weight;
195  }
196  if (mtot > 0.0) {
197  *cx /= mtot;
198  *cy /= mtot;
199  *cz /= mtot;
200  } else {
201  error_divide_zero(16020, program);
202  }
203  return mtot;
204 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
float weight
Definition: situs.h:52
float x
Definition: situs.h:46
void error_divide_zero(int error_number, const char *program)
Definition: lib_err.cpp:217

◆ calc_mass()

double calc_mass ( PDB pdb0,
unsigned  num_atoms 
)

Definition at line 165 of file lib_pwk.cpp.

166 {
167  unsigned id;
168  double mtot;
169 
170  mtot = 0.0;
171  for (id = 0; id < num_atoms; ++id)
172  mtot += pdb0[id].weight;
173  return mtot;
174 }

◆ calc_sphere()

double calc_sphere ( PDB pdb0,
unsigned  num_atoms,
double  cx,
double  cy,
double  cz 
)

Definition at line 209 of file lib_pwk.cpp.

211 {
212  unsigned id;
213  double maxradius, currradius;
214  const char *program = "lib_pwk";
215  const char *shape = "sphere";
216 
217  maxradius = -1e20;
218  for (id = 0; id < num_atoms; ++id) {
219  currradius = ((pdb0[id].x - cx) * (pdb0[id].x - cx) +
220  (pdb0[id].y - cy) * (pdb0[id].y - cy) +
221  (pdb0[id].z - cz) * (pdb0[id].z - cz));
222  if (currradius > maxradius) maxradius = currradius;
223  }
224 
225  if (maxradius >= 0.0) {
226  maxradius = sqrt(maxradius);
227  return maxradius;
228  } else {
229  error_no_bounding(16030, program, shape);
230  return 0;
231  }
232 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
void sqrt(Image< double > &op)
static double * y
void error_no_bounding(int error_number, const char *program, const char *shape)
Definition: lib_err.cpp:258
doublereal * x
float x
Definition: situs.h:46
double z

◆ check_if_inside()

int check_if_inside ( double  fraction,
double  widthx,
double  widthy,
double  widthz,
unsigned  extx,
unsigned  exty,
unsigned  extz,
PDB inpdb,
unsigned  num_atoms,
double  shift[3] 
)

Definition at line 441 of file lib_pwk.cpp.

444 {
445 
446  int x0, y0, z0, x1, y1, z1, i;
447  double gx, gy, gz;
448  int outside_error = 0;
449 
450  for (i = 0; i < num_atoms; ++i) {
451  /* compute position within grid */
452  gx = extx / 2.0 + (inpdb[i].x + shift[0]) / widthx;
453  gy = exty / 2.0 + (inpdb[i].y + shift[1]) / widthy;
454  gz = extz / 2.0 + (inpdb[i].z + shift[2]) / widthz;
455  x0 = floor(gx);
456  y0 = floor(gy);
457  z0 = floor(gz);
458  x1 = x0 + 1;
459  y1 = y0 + 1;
460  z1 = z0 + 1;
461  if (x0 < 0) ++outside_error;
462  if (y0 < 0) ++outside_error;
463  if (z0 < 0) ++outside_error;
464  if (x1 >= extx) ++outside_error;
465  if (y1 >= exty) ++outside_error;
466  if (z1 >= extz) ++outside_error;
467  }
468  if (outside_error >= (1 - fraction)*num_atoms) return outside_error;
469  else return 0;
470 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
__host__ __device__ float2 floor(const float2 v)
#define z0
#define i
#define y0
#define x0
float x
Definition: situs.h:46

◆ copy_atoms()

void copy_atoms ( PDB pdb_original,
PDB pdb_duplicate,
int  io,
int  id,
int  ktot 
)

Definition at line 25 of file lib_pwk.cpp.

27 {
28  int j, k;
29  for (k = 0; k < ktot; ++k) {
30  pdb_duplicate[id + k].weight = pdb_original[io].weight;
31  pdb_duplicate[id + k].x = pdb_original[io].x;
32  pdb_duplicate[id + k].y = pdb_original[io].y;
33  pdb_duplicate[id + k].z = pdb_original[io].z;
34  for (j = 0; j < 4; ++j) pdb_duplicate[id + k].segid[j] = pdb_original[io].segid[j];
35  pdb_duplicate[id + k].serial = id + k + 1;
36  for (j = 0; j < 7; ++j) pdb_duplicate[id + k].recd[j] = pdb_original[io].recd[j];
37  for (j = 0; j < 3; ++j) pdb_duplicate[id + k].type[j] = pdb_original[io].type[j];
38  for (j = 0; j < 3; ++j) pdb_duplicate[id + k].loc[j] = pdb_original[io].loc[j];
39  for (j = 0; j < 2; ++j) pdb_duplicate[id + k].alt[j] = pdb_original[io].alt[j];
40  for (j = 0; j < 5; ++j) pdb_duplicate[id + k].res[j] = pdb_original[io].res[j];
41  for (j = 0; j < 2; ++j) pdb_duplicate[id + k].chain[j] = pdb_original[io].chain[j];
42  pdb_duplicate[id + k].seq = pdb_original[io].seq;
43  for (j = 0; j < 2; ++j) pdb_duplicate[id + k].icode[j] = pdb_original[io].icode[j];
44  pdb_duplicate[id + k].occupancy = pdb_original[io].occupancy;
45  pdb_duplicate[id + k].beta = pdb_original[io].beta;
46  pdb_duplicate[id + k].footnote = pdb_original[io].footnote;
47  for (j = 0; j < 3; ++j) pdb_duplicate[id + k].element[j] = pdb_original[io].element[j];
48  for (j = 0; j < 3; ++j) pdb_duplicate[id + k].charge[j] = pdb_original[io].charge[j];
49  }
50  return;
51 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
float weight
Definition: situs.h:52
int footnote
Definition: situs.h:51
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
int seq
Definition: situs.h:45
viol type
float occupancy
Definition: situs.h:49
float x
Definition: situs.h:46
glob_prnt io
#define j
int serial
Definition: situs.h:44
float beta
Definition: situs.h:50

◆ project_mass()

void project_mass ( double **  outmap,
unsigned long  nvox,
double  widthx,
double  widthy,
double  widthz,
unsigned  extx,
unsigned  exty,
unsigned  extz,
PDB inpdb,
unsigned  num_atoms,
double  shift[3],
unsigned  ignored[3] 
)

Definition at line 269 of file lib_pwk.cpp.

271 {
272  int x0, y0, z0, x1, y1, z1, i;
273  double gx, gy, gz;
274  double a, b, c;
275  double ab, ab1, a1b, a1b1;
276  int marglx, margly, marglz, margux, marguy, marguz;
277 
278  /* lower margins */
279  marglx = ignored[0];
280  margly = ignored[1];
281  marglz = ignored[2];
282 
283  /* one voxel safety buffer for boundary test below */
284  if (marglx == 0) marglx = 1;
285  if (margly == 0) margly = 1;
286  if (marglz == 0) marglz = 1;
287 
288  /* upper margins */
289  margux = extx - marglx;
290  marguy = exty - margly;
291  marguz = extz - marglz;
292 
293  zero_vect(*outmap, nvox);
294  for (i = 0; i < num_atoms; ++i) {
295 
296  /* compute position within grid in voxel units */
297 
298  gx = extx / 2.0 + (inpdb[i].x + shift[0]) / widthx;
299  gy = exty / 2.0 + (inpdb[i].y + shift[1]) / widthy;
300  gz = extz / 2.0 + (inpdb[i].z + shift[2]) / widthz;
301  x0 = floor(gx);
302  y0 = floor(gy);
303  z0 = floor(gz);
304  x1 = x0 + 1;
305  y1 = y0 + 1;
306  z1 = z0 + 1;
307 
308  /* boundary check considers if probe cube overlaps at least partially with allowed region (stripped of ignored) */
309  if (x1 >= marglx && x0 < margux && y1 >= margly && y0 < marguy && z1 >= marglz && z0 < marguz) {
310 
311  /* interpolate */
312  a = x1 - gx;
313  b = y1 - gy;
314  c = z1 - gz;
315  ab = a * b;
316  ab1 = a * (1 - b);
317  a1b = (1 - a) * b;
318  a1b1 = (1 - a) * (1 - b);
319  a = (1 - c);
320  *(*outmap + gidz_general(z0, y0, x0, exty, extx)) += ab * c * inpdb[i].weight;
321  *(*outmap + gidz_general(z1, y0, x0, exty, extx)) += ab * a * inpdb[i].weight;
322  *(*outmap + gidz_general(z0, y1, x0, exty, extx)) += ab1 * c * inpdb[i].weight;
323  *(*outmap + gidz_general(z1, y1, x0, exty, extx)) += ab1 * a * inpdb[i].weight;
324  *(*outmap + gidz_general(z0, y0, x1, exty, extx)) += a1b * c * inpdb[i].weight;
325  *(*outmap + gidz_general(z1, y0, x1, exty, extx)) += a1b * a * inpdb[i].weight;
326  *(*outmap + gidz_general(z0, y1, x1, exty, extx)) += a1b1 * c * inpdb[i].weight;
327  *(*outmap + gidz_general(z1, y1, x1, exty, extx)) += a1b1 * a * inpdb[i].weight;
328  }
329  }
330 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
__host__ __device__ float2 floor(const float2 v)
doublereal * c
unsigned long gidz_general(int k, int j, int i, unsigned ey, unsigned ex)
Definition: lib_vwk.cpp:32
#define z0
#define i
doublereal * b
#define y0
#define x0
float x
Definition: situs.h:46
void zero_vect(double *vect, unsigned long len)
Definition: lib_vec.cpp:26
doublereal * a

◆ project_mass_convolve_kernel_corr()

void project_mass_convolve_kernel_corr ( double  widthx,
double  widthy,
double  widthz,
unsigned  extx,
unsigned  exty,
unsigned  extz,
PDB inpdb,
unsigned  num_atoms,
double  shift[3],
double *  kernel,
unsigned  kernel_size,
double  normfac,
unsigned  ignored[3],
double *  lowmap,
double *  corr_hi_low 
)

Definition at line 338 of file lib_pwk.cpp.

343 {
344  int x0, y0, z0, x1, y1, z1, i;
345  double gx, gy, gz;
346  double a, b, c;
347  double ab, ab1, a1b, a1b1;
348  double dval1, dval2, dval3, dval4, dval5, dval6, dval7, dval8;
349  int indx2, indy2, indz2, margin;
350  double inv_normfac;
351  double *idx_kern;
352  double *idx_low1, *idx_low2, *idx_low3, *idx_low4;
353  int marglx, margly, marglz, margux, marguy, marguz;
354 
355  inv_normfac = 1.0 / normfac;
356  margin = (kernel_size - 1) / 2;
357 
358  /* lower margins */
359  marglx = margin + ignored[0];
360  margly = margin + ignored[1];
361  marglz = margin + ignored[2];
362 
363  /* one voxel safety buffer for boundary test below */
364  if (marglx == margin) marglx += 1;
365  if (margly == margin) margly += 1;
366  if (marglz == margin) marglz += 1;
367 
368  /* upper margins */
369  margux = extx - marglx;
370  marguy = exty - margly;
371  marguz = extz - marglz;
372 
373  for (i = 0; i < num_atoms; ++i) {
374 
375  /* compute position of probe cube within grid in voxel units */
376 
377  gx = extx / 2.0 + (inpdb[i].x + shift[0]) / widthx;
378  gy = exty / 2.0 + (inpdb[i].y + shift[1]) / widthy;
379  gz = extz / 2.0 + (inpdb[i].z + shift[2]) / widthz;
380  x0 = floor(gx);
381  y0 = floor(gy);
382  z0 = floor(gz);
383  x1 = x0 + 1;
384  y1 = y0 + 1;
385  z1 = z0 + 1;
386 
387  /* boundary check considers if probe cube overlaps at least partially with original map (stripped of margin + ignored) */
388  if (x1 >= marglx && x0 < margux && y1 >= margly && y0 < marguy && z1 >= marglz && z0 < marguz) {
389 
390  /* interpolate */
391  a = x1 - gx;
392  b = y1 - gy;
393  c = z1 - gz;
394  ab = a * b;
395  ab1 = a * (1 - b);
396  a1b = (1 - a) * b;
397  a1b1 = (1 - a) * (1 - b);
398  a = (1 - c);
399 
400  dval1 = inv_normfac * ab * c * inpdb[i].weight;
401  dval2 = inv_normfac * a1b * c * inpdb[i].weight;
402  dval3 = inv_normfac * ab1 * c * inpdb[i].weight;
403  dval4 = inv_normfac * a1b1 * c * inpdb[i].weight;
404  dval5 = inv_normfac * ab * a * inpdb[i].weight;
405  dval6 = inv_normfac * a1b * a * inpdb[i].weight;
406  dval7 = inv_normfac * ab1 * a * inpdb[i].weight;
407  dval8 = inv_normfac * a1b1 * a * inpdb[i].weight;
408 
409  idx_kern = kernel;
410  for (indz2 = -margin; indz2 <= margin; indz2++) {
411  for (indy2 = -margin; indy2 <= margin; indy2++) {
412  idx_low1 = lowmap + (exty * extx * (z0 + indz2) + extx * (y0 + indy2) + x0);
413  idx_low2 = lowmap + (exty * extx * (z0 + indz2) + extx * (y1 + indy2) + x0);
414  idx_low3 = lowmap + (exty * extx * (z1 + indz2) + extx * (y0 + indy2) + x0);
415  idx_low4 = lowmap + (exty * extx * (z1 + indz2) + extx * (y1 + indy2) + x0);
416  for (indx2 = -margin; indx2 <= margin; indx2++) {
417  *corr_hi_low += *(idx_kern) *
418  (dval1 * *(idx_low1 + indx2)
419  + dval2 * *(idx_low1 + indx2 + 1)
420  + dval3 * *(idx_low2 + indx2)
421  + dval4 * *(idx_low2 + indx2 + 1)
422  + dval5 * *(idx_low3 + indx2)
423  + dval6 * *(idx_low3 + indx2 + 1)
424  + dval7 * *(idx_low4 + indx2)
425  + dval8 * *(idx_low4 + indx2 + 1));
426  idx_kern++;
427  }
428  }
429  } /* end of loop over kernel */
430 
431  } /* end of boundary check */
432 
433  } /* end of loop over atoms */
434 
435 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
__host__ __device__ float2 floor(const float2 v)
doublereal * c
float weight
Definition: situs.h:52
#define z0
#define i
doublereal * b
#define y0
#define x0
float x
Definition: situs.h:46
doublereal * a

◆ rot_axis()

void rot_axis ( PDB pdb_original,
PDB pdb_rotate,
unsigned  num_atoms,
char  axis,
double  angle 
)

Definition at line 56 of file lib_pwk.cpp.

58 {
59  unsigned id;
60  double x, y, z;
61  double sint = sin(angle);
62  double cost = cos(angle);
63 
64  switch (axis) {
65  case ('X'):
66  for (id = 0; id < num_atoms; ++id) {
67  y = pdb_original[id].y;
68  z = pdb_original[id].z;
69  pdb_rotate[id].x = pdb_original[id].x;
70  pdb_rotate[id].y = (cost * y + sint * z);
71  pdb_rotate[id].z = (cost * z - sint * y);
72  }
73  break;
74  case ('Y'):
75  for (id = 0; id < num_atoms; ++id) {
76  x = pdb_original[id].x;
77  z = pdb_original[id].z;
78  pdb_rotate[id].y = pdb_original[id].y;
79  pdb_rotate[id].x = (cost * x + sint * z);
80  pdb_rotate[id].z = (cost * z - sint * x);
81  }
82  break;
83  case ('Z'):
84  for (id = 0; id < num_atoms; ++id) {
85  x = pdb_original[id].x;
86  y = pdb_original[id].y;
87  pdb_rotate[id].z = pdb_original[id].z;
88  pdb_rotate[id].x = (cost * x + sint * y);
89  pdb_rotate[id].y = (cost * y - sint * x);
90  }
91  break;
92  }
93 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
static double * y
doublereal * x
char axis
float x
Definition: situs.h:46
double z

◆ rot_euler()

void rot_euler ( PDB pdb_original,
PDB pdb_rotate,
unsigned  num_atoms,
double  psi,
double  theta,
double  phi 
)

Definition at line 98 of file lib_pwk.cpp.

100 {
101  unsigned id;
102  double rot_matrix[3][3], currx, curry, currz;
103 
104  get_rot_matrix(rot_matrix, psi, theta, phi);
105 
106  for (id = 0; id < num_atoms; ++id) {
107  currx = pdb_original[id].x;
108  curry = pdb_original[id].y;
109  currz = pdb_original[id].z;
110  pdb_rotate[id].x = currx * rot_matrix[0][0] +
111  curry * rot_matrix[0][1] +
112  currz * rot_matrix[0][2];
113  pdb_rotate[id].y = currx * rot_matrix[1][0] +
114  curry * rot_matrix[1][1] +
115  currz * rot_matrix[1][2];
116  pdb_rotate[id].z = currx * rot_matrix[2][0] +
117  curry * rot_matrix[2][1] +
118  currz * rot_matrix[2][2];
119  }
120 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
double theta
float x
Definition: situs.h:46
void get_rot_matrix(double dump[3][3], double psi, double theta, double phi)
Definition: lib_eul.cpp:184
double psi(const double x)

◆ translate()

void translate ( PDB pdb_original,
PDB pdb_move,
unsigned  num_atoms,
double  x0,
double  y0,
double  z0 
)

Definition at line 125 of file lib_pwk.cpp.

127 {
128  unsigned id;
129 
130  for (id = 0; id < num_atoms; ++id) {
131  pdb_move[id].x = pdb_original[id].x + x0;
132  pdb_move[id].y = pdb_original[id].y + y0;
133  pdb_move[id].z = pdb_original[id].z + z0;
134  }
135 }
float y
Definition: situs.h:47
float z
Definition: situs.h:48
#define z0
#define y0
#define x0
float x
Definition: situs.h:46