00001 /*************************************************************************** 00002 *cr 00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 00004 *cr University of Illinois 00005 *cr All Rights Reserved 00006 *cr 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * RCS INFORMATION: 00011 * 00012 * $RCSfile: Measure.h,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.76 $ $Date: 2024年05月16日 20:21:10 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * Code to measure atom distances, angles, dihedrals, etc. 00019 ***************************************************************************/ 00020 00021 #include "ResizeArray.h" 00022 #include "Molecule.h" 00023 00024 class AtomSel; 00025 class Matrix4; 00026 class MoleculeList; 00027 00028 #define MEASURE_NOERR 0 00029 #define MEASURE_ERR_NOSEL -1 00030 #define MEASURE_ERR_NOATOMS -2 00031 #define MEASURE_ERR_BADWEIGHTNUM -3 00032 #define MEASURE_ERR_NOWEIGHT -4 00033 #define MEASURE_ERR_NOCOM -4 00034 #define MEASURE_ERR_NOMINMAXCOORDS -4 00035 #define MEASURE_ERR_NORGYR -4 00036 #define MEASURE_ERR_BADWEIGHTSUM -5 00037 #define MEASURE_ERR_NOMOLECULE -6 00038 #define MEASURE_ERR_BADWEIGHTPARM -7 00039 #define MEASURE_ERR_NONNUMBERPARM -8 00040 #define MEASURE_ERR_MISMATCHEDCNT -9 00041 #define MEASURE_ERR_RGYRMISMATCH -10 00042 #define MEASURE_ERR_NOFRAMEPOS -11 00043 #define MEASURE_ERR_NONZEROJACOBI -12 00044 #define MEASURE_ERR_GENERAL -13 00045 #define MEASURE_ERR_NORADII -14 00046 #define MEASURE_ERR_BADORDERINDEX -15 00047 #define MEASURE_ERR_NOTSUP -16 00048 #define MEASURE_ERR_BADFRAMERANGE -17 00049 #define MEASURE_ERR_MISMATCHEDMOLS -18 00050 #define MEASURE_ERR_REPEATEDATOM -19 00051 #define MEASURE_ERR_NOFRAMES -20 00052 #define MEASURE_ERR_BADATOMID -21 00053 #define MEASURE_ERR_BADCUTOFF -22 00054 #define MEASURE_ERR_ZEROGRIDSIZE -23 00055 00056 #define MEASURE_BOND 2 00057 #define MEASURE_ANGLE 3 00058 #define MEASURE_DIHED 4 00059 #define MEASURE_IMPRP 5 00060 #define MEASURE_VDW 6 00061 #define MEASURE_ELECT 7 00062 00063 // symbolic flags for cluster analysis 00064 enum {MEASURE_DIST_RMSD=0, 00065 MEASURE_DIST_FITRMSD, 00066 MEASURE_DIST_RGYRD, 00067 MEASURE_DIST_CUSTOM, 00068 MEASURE_NUM_DIST}; 00069 00070 extern const char *measure_error(int errnum); 00071 00072 // apply a matrix transformation to the coordinates of a selection 00073 extern int measure_move(const AtomSel *sel, float *framepos, 00074 const Matrix4 &mat); 00075 00076 // Calculate average position of selected atoms over selected frames 00077 extern int measure_avpos(const AtomSel *sel, MoleculeList *mlist, 00078 int start, int end, int step, float *avpos); 00079 00080 // find the center of the selected atoms in sel, using the coordinates in 00081 // framepos and the weights in weight. weight has sel->selected elements. 00082 // Place the answer in com. Return 0 on success, or negative on error. 00083 extern int measure_center(const AtomSel *sel, const float *framepos, 00084 const float *weight, float *com); 00085 00086 // find the center of the selected atoms in sel, using the coordinates in 00087 // framepos and the weights in weight. weight has sel->selected elements. 00088 // Place the answer in com. Return number of residues on success, 00089 // or negative on error. 00090 extern int measure_center_perresidue(MoleculeList *mlist, const AtomSel *sel, 00091 const float *framepos, 00092 const float *weight, float *com); 00093 00094 // find the dipole of the selected atoms in sel, using the coordinates in 00095 // framepos, and the atom charges. 00096 // Place the answer in dipole. Return 0 on success, or negative on error. 00097 // Default units are elementary charges/Angstrom. 00098 // Setting unitsdebye to 1 will scale the results by 4.77350732929 (default 0) 00099 // For charged systems the point of reference for computing the dipole is 00100 // selected with usecenter: 00101 // 1 = geometrical center (default), 00102 // -1 = center of mass, 00103 // 0 = origin 00104 extern int measure_dipole(const AtomSel *sel, MoleculeList *mlist, 00105 float *dipole, int unitsdebye, int usecenter); 00106 00107 extern int measure_hbonds(Molecule *mol, AtomSel *sel1, AtomSel *sel2, 00108 double cut, double maxangle, int *donlist, 00109 int *hydlist, int *acclist, int maxsize); 00110 00111 // Find the transformation which aligns the atoms of sel1 and sel2 optimally, 00112 // meaning it minimizes the RMS distance between them, weighted by weight. 00113 // The returned matrix will have positive determinant, even if an optimal 00114 // alignment would produce a matrix with negative determinant; the last 00115 // row in the matrix is flipped to change the sign of the determinant. 00116 // sel1->selected == sel2->selected == len(weight). 00117 extern int measure_fit(const AtomSel *sel1, const AtomSel *sel2, 00118 const float *x, const float *y, const float *weight, 00119 const int *order, Matrix4 *mat); 00120 00121 // compute the axis aligned aligned bounding box for the selected atoms 00122 // returns 0 if success 00123 // returns <0 if not 00124 // If the selection contains no atoms, return {0 0 0} for min and max. 00125 extern int measure_minmax(int num, const int *on, const float *framepos, 00126 const float *radii, 00127 float *min_coord, float *max_coord); 00128 00129 // Calculate the radius of gyration of the given selection, using the 00130 // given weight, placing the result in rgyr. 00131 extern int measure_rgyr(const AtomSel *sel, MoleculeList *mlist, 00132 const float *weight, float *rgyr); 00133 00134 // Calculate the RMS distance between the atoms in the two selections, 00135 // weighted by weight. Same conditions on sel1, sel2, and weight as for 00136 // measure_fit. 00137 extern int measure_rmsd(const AtomSel *sel1, const AtomSel *sel2, 00138 int num, const float *f1, const float *f2, 00139 float *weight, float *rmsd); 00140 00141 // Calculate the RMS distance between the atoms in the two selections, 00142 // weighted by weight. Same conditions on sel1, sel2, and weight as for 00143 // measure_fit. Now done per residue (so the pointer is expected to be an array) 00144 extern int measure_rmsd_perresidue(const AtomSel *sel1, const AtomSel *sel2, 00145 MoleculeList *mlist, int num, 00146 float *weight, float *rmsd); 00147 00148 // Measure RMS distance between two selections as with measure_rmsd(), 00149 // except that it is computed with an implicit best-fit alignment 00150 // by virtue of the QCP algorithm. 00151 extern int measure_rmsd_qcp(VMDApp *app, 00152 const AtomSel *sel1, const AtomSel *sel2, 00153 int num, const float *f1, const float *f2, 00154 float *weight, float *rmsd); 00155 00156 // Measure matrix of RMS distance between all selected trajectory frames, 00157 // computed with an implicit best-fit alignment by virtue of the QCP algorithm. 00158 extern int measure_rmsdmat_qcp(VMDApp *app, 00159 const AtomSel *sel, MoleculeList *mlist, 00160 int num, float *weight, 00161 int start, int end, int step, 00162 float *rmsd); 00163 00164 // Measure matrix of RMS distance between all selected trajectory frames, 00165 // computed with an implicit best-fit alignment by virtue of the QCP algorithm. 00166 // XXX Experimental out-of-core I/O algorithms 00167 extern int measure_rmsdmat_qcp_ooc(VMDApp *app, 00168 const AtomSel *sel, MoleculeList *mlist, 00169 int nfiles, const char **trjfileset, 00170 int num, float *weight, 00171 int start, int end, int step, 00172 int &framecount, float *&rmsd); 00173 00174 // Given the component sums of QCP inner products, uses 00175 // QCP algorithm to solve for best-fit RMSD and rotational alignment 00176 extern int FastCalcRMSDAndRotation(double *rot, double *A, float *rmsd, 00177 double E0, int len, double minScore); 00178 00179 // Calculate RMS fluctuation of selected atoms over selected frames 00180 extern int measure_rmsf(const AtomSel *sel, MoleculeList *mlist, 00181 int start, int end, int step, float *rmsf); 00182 00183 // Calculate RMSF per residue. 00184 extern int measure_rmsf_perresidue(const AtomSel *sel, MoleculeList *mlist, 00185 int start, int end, int step, float *rmsf); 00186 00187 extern int measure_sumweights(const AtomSel *sel, int numweights, 00188 const float *weights, float *weightsum); 00189 00190 00191 // find the solvent-accessible surface area of atoms in the given selection. 00192 // Use the assigned radii for each atom, and extend this radius by the 00193 // parameter srad to find the points on a sphere that are exposed to solvent. 00194 // Optional parameters (pass NULL to ignore) are: 00195 // pts: fills the given array with the location of the points that make 00196 // up the surface. 00197 // restrictsel: Only solvent accessible points near the given selection will 00198 // be considered. 00199 // nsamples: number of points to use around each atom. 00200 extern int measure_sasa(const AtomSel *sel, const float *framepos, 00201 const float *radius, float srad, float *sasa, ResizeArray<float> *pts, 00202 const AtomSel *restrictsel, const int *nsamples); 00203 00204 // XXX experimental version that processes a list of selections at a time 00205 extern int measure_sasalist(MoleculeList *mlist, 00206 const AtomSel **sellist, int numsels, 00207 float srad, float *sasalist, const int *nsamples); 00208 00209 // compute per-residue SASA 00210 extern int measure_sasa_perresidue(const AtomSel *sel, const float *framepos, 00211 const float *radius, float srad, float *sasa, 00212 ResizeArray<float> *sasapts, const AtomSel *restrictsel, 00213 const int *nsamples, int *rescount,MoleculeList *mlist); 00214 00215 00216 // perform cluster analysis 00217 extern int measure_cluster(AtomSel *sel, MoleculeList *mlist, 00218 const int numcluster, const int algorithm, 00219 const int likeness, const double cutoff, 00220 int *clustersize, int **clusterlist, 00221 int first, int last, int step, int selupdate, 00222 float *weights); 00223 00224 // perform cluster size analysis 00225 extern int measure_clustsize(const AtomSel *sel, MoleculeList *mlist, 00226 const double cutoff, int *clustersize, 00227 int *clusternum, int *clusteridx, 00228 int minsize, int numshared, int usepbc); 00229 00230 // calculate g(r) for two selections 00231 extern int measure_gofr(AtomSel *sel1, AtomSel *sel2, 00232 MoleculeList *mlist, 00233 const int count_h, double *gofr, double *numint, 00234 double *histog, const float delta, 00235 int first, int last, int step, int *framecntr, 00236 int usepbc, int selupdate); 00237 00238 // calculate g(r) for two selections 00239 extern int measure_rdf(VMDApp *app, 00240 AtomSel *sel1, AtomSel *sel2, 00241 MoleculeList *mlist, 00242 const int count_h, double *gofr, double *numint, 00243 double *histog, const float delta, 00244 int first, int last, int step, int *framecntr, 00245 int usepbc, int selupdate); 00246 00247 int measure_geom(MoleculeList *mlist, int *molid, int *atmid, ResizeArray<float> *gValues, 00248 int frame, int first, int last, int defmolid, int geomtype); 00249 00250 // calculate the value of this geometry, and return it 00251 int calculate_bond(MoleculeList *mlist, int *molid, int *atmid, float *value); 00252 00253 // calculate the value of this geometry, and return it 00254 int calculate_angle(MoleculeList *mlist, int *molid, int *atmid, float *value); 00255 00256 // calculate the value of this geometry, and return it 00257 int calculate_dihed(MoleculeList *mlist, int *molid, int *atmid, float *value); 00258 00259 // check whether the given molecule & atom index is OK 00260 // if OK, return Molecule pointer; otherwise, return NULL 00261 int check_mol(Molecule *mol, int a); 00262 00263 // for the given Molecule, find the UNTRANSFORMED coords for the given atom 00264 // return Molecule pointer if successful, NULL otherwise. 00265 int normal_atom_coord(Molecule *m, int a, float *pos); 00266 00267 int measure_energy(MoleculeList *mlist, int *molid, int *atmid, int natoms, ResizeArray<float> *gValues, 00268 int frame, int first, int last, int defmolid, double *params, int geomtype); 00269 int compute_bond_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, 00270 float k, float x0); 00271 int compute_angle_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, 00272 float k, float x0, float kub, float s0); 00273 int compute_dihed_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, 00274 float k, int n, float delta); 00275 int compute_imprp_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, 00276 float k, float x0); 00277 int compute_vdw_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, 00278 float rmin1, float eps1, float rmin2, float eps2, float cutoff, float switchdist); 00279 int compute_elect_energy(MoleculeList *mlist, int *molid, int *atmid, float *energy, 00280 float q1, float q2, bool flag1, bool flag2, float cutoff); 00281 00282 // compute matrix that transforms coordinates from an arbitrary PBC cell 00283 // into an orthonormal unitcell. 00284 int measure_pbc2onc(MoleculeList *mlist, int molid, int frame, const float *center, Matrix4 &transform); 00285 00286 // does the low level work for the above 00287 void get_transform_to_orthonormal_cell(const float *cell, const float center[3], Matrix4 &transform); 00288 00289 // get atoms in PBC neighbor cells 00290 int measure_pbc_neighbors(MoleculeList *mlist, AtomSel *sel, int molid, 00291 int frame, const Matrix4 *alignment, 00292 const float *center, const float *cutoff, const float *box, 00293 ResizeArray<float> *extcoord_array, 00294 ResizeArray<int> *indexmap_array); 00295 00296 // compute the orthogonalized bounding box for the PBC cell. 00297 int compute_pbcminmax(MoleculeList *mlist, int molid, int frame, 00298 const float *center, const Matrix4 *transform, 00299 float *min, float *max); 00300 00301 00302 // Return the list of atoms within the specified distance of the surface 00303 // where the surface depth is specified by sel_dist, the grid resolution is 00304 // approximately gridsz, and atoms are assume to have size radius 00305 // If any of a, b, c, alpha, or gamma are zero, assume non-periodic, 00306 // otherwise assume periodic 00307 // returns 0 if success 00308 // returns <0 if not 00309 extern int measure_surface(AtomSel *sel, MoleculeList *mlist, 00310 const float *framepos, 00311 const double gridsz, 00312 const double radius, 00313 const double sel_dist, 00314 int **surface, int *n_surf); 00315 00316 00317 // Calculate center of mass, principle axes and moments of inertia for 00318 // selected atoms. The corresponding eigenvalues are also returned, 00319 // and might tell you if two axes are equivalent. 00320 extern int measure_inertia(AtomSel *sel, MoleculeList *mlist, const float *coor, float rcom[3], 00321 float priaxes[3][3], float itensor[4][4], float evalue[3]); 00322