00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include <stdio.h>
00036 #include <math.h>
00037 #include "utilities.h"
00038 #include "Measure.h"
00039 #include "MeasureSymmetry.h"
00040 #include "AtomSel.h"
00041 #include "Matrix4.h"
00042 #include "MoleculeList.h"
00043 #include "SpatialSearch.h"
00044 #include "MoleculeGraphics.h"
00045 #include "Inform.h"
00046 #include "WKFUtils.h"
00047
00048 #define POINTGROUP_UNKNOWN 0
00049 #define POINTGROUP_C1 1
00050 #define POINTGROUP_CN 2
00051 #define POINTGROUP_CNV 3
00052 #define POINTGROUP_CNH 4
00053 #define POINTGROUP_CINFV 5
00054 #define POINTGROUP_DN 6
00055 #define POINTGROUP_DND 7
00056 #define POINTGROUP_DNH 8
00057 #define POINTGROUP_DINFH 9
00058 #define POINTGROUP_CI 10
00059 #define POINTGROUP_CS 11
00060 #define POINTGROUP_S2N 12
00061 #define POINTGROUP_T 13
00062 #define POINTGROUP_TD 14
00063 #define POINTGROUP_TH 15
00064 #define POINTGROUP_O 16
00065 #define POINTGROUP_OH 17
00066 #define POINTGROUP_I 18
00067 #define POINTGROUP_IH 19
00068 #define POINTGROUP_KH 20
00069
00070 #define OVERLAPCUTOFF 0.4
00071
00072
00073
00074
00075
00076
00077 #define BONDSUMTOL 1.5
00078
00079
00080
00081 #define MINATOMDIST 0.6
00082
00083
00084 #define VERTICALPLANE 1
00085 #define DIHEDRALPLANE 3
00086 #define HORIZONTALPLANE 4
00087
00088
00089 #define INFINITE_ORDER -1
00090 #define PRELIMINARY_C2 -2
00091
00092
00093 #define TETRAHEDRON -4
00094 #define OCTAHEDRON -8
00095 #define DODECAHEDRON -12
00096
00097
00098
00099 #if defined(_MSC_VER)
00100
00101 static double myerfc(double x) {
00102 double p, a1, a2, a3, a4, a5;
00103 double t, erfcx;
00104
00105 p = 0.3275911;
00106 a1 = 0.254829592;
00107 a2 = -0.284496736;
00108 a3 = 1.421413741;
00109 a4 = -1.453152027;
00110 a5 = 1.061405429;
00111
00112 t = 1.0 / (1.0 + p*x);
00113 erfcx = ((a1 + (a2 + (a3 + (a4 + a5*t)*t)*t)*t)*t) * exp(-pow(x,2.0));
00114 return erfcx;
00115 }
00116 #endif
00117
00118
00119
00120
00121 static inline bool collinear(const float *axis1, const float *axis2, float tol);
00122 static inline bool orthogonal(const float *axis1, const float *axis2, float tol);
00123 static inline bool behind_plane(const float *normal, const float *coor);
00124 static int isprime(int x);
00125 static int numprimefactors(int x);
00126 static void align_plane_with_axis(const float *normal, const float *axis, float *alignedplane);
00127 static void assign_atoms(AtomSel *sel, MoleculeList *mlist, float *(&mycoor), int *(&atomtype));
00128 static float trans_overlap(int *atomtype, float *(&coor), int numcoor,
00129 const Matrix4 *trans, float sigma,
00130 bool skipident, int maxnatoms, float &overlappermatch);
00131 static inline float trans_overlap(int *atomtype, float *(&coor), int numcoor,
00132 const Matrix4 *trans, float sigma,
00133 bool skipident, int maxnatoms);
00134
00135
00136
00137 typedef struct {
00138 int pointgroup;
00139 int pointgrouporder;
00140 float rmsd;
00141 float *idealcoor;
00142 Plane *planes;
00143 Axis *axes;
00144 Axis *rotreflections;
00145 bool linear;
00146 bool planar;
00147 bool inversion;
00148 bool sphericaltop;
00149 bool symmetrictop;
00150 } Best;
00151
00152
00153
00154 typedef struct {
00155 char name[6];
00156 ElementSummary summary;
00157 } PointGroupDefinition;
00158
00159 PointGroupDefinition pgdefinitions[] = {
00160
00161 { "C1", {0, 0, 0, {0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00162 { "Cs", {0, 1, 0, {0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00163 { "Ci", {1, 0, 0, {0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00164 { "C2", {0, 0, 0, {0, 0, 1, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00165 { "C3", {0, 0, 0, {0, 0, 0, 1, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00166 { "C4", {0, 0, 0, {0, 0, 1, 0, 1, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00167 { "C5", {0, 0, 0, {0, 0, 0, 0, 0, 1, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00168 { "C6", {0, 0, 0, {0, 0, 1, 1, 0, 0, 1, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00169 { "C7", {0, 0, 0, {0, 0, 0, 0, 0, 0, 0, 1, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00170 { "C8", {0, 0, 0, {0, 0, 1, 0, 1, 0, 0, 0, 1}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00171 { "D2", {0, 0, 0, {0, 0, 3, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00172 { "D3", {0, 0, 0, {0, 0, 3, 1, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00173 { "D4", {0, 0, 0, {0, 0, 5, 0, 1, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00174 { "D5", {0, 0, 0, {0, 0, 5, 0, 0, 1, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00175 { "D6", {0, 0, 0, {0, 0, 7, 1, 0, 0, 1, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00176 { "D7", {0, 0, 0, {0, 0, 7, 0, 0, 0, 0, 1, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00177 { "D8", {0, 0, 0, {0, 0, 9, 0, 1, 0, 0, 0, 1}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00178 { "C2v", {0, 2, 0, {0, 0, 1, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00179 { "C3v", {0, 3, 0, {0, 0, 0, 1, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00180 { "C4v", {0, 4, 0, {0, 0, 1, 0, 1, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00181 { "C5v", {0, 5, 0, {0, 0, 0, 0, 0, 1, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00182 { "C6v", {0, 6, 0, {0, 0, 1, 1, 0, 0, 1, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00183 { "C7v", {0, 7, 0, {0, 0, 0, 0, 0, 0, 0, 1, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00184 { "C8v", {0, 8, 0, {0, 0, 1, 0, 1, 0, 0, 0, 1}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00185 { "C2h", {1, 1, 0, {0, 0, 1, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00186 { "C3h", {0, 1, 0, {0, 0, 0, 1, 0, 0, 0, 0, 0}, {1,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00187 { "C4h", {1, 1, 0, {0, 0, 1, 0, 1, 0, 0, 0, 0}, {0,1,0,0,0,0,0,0,0,0,0,0,0,0}} },
00188 { "C5h", {0, 1, 0, {0, 0, 0, 0, 0, 1, 0, 0, 0}, {0,0,1,0,0,0,0,0,0,0,0,0,0,0}} },
00189 { "C6h", {1, 1, 0, {0, 0, 1, 1, 0, 0, 1, 0, 0}, {0,0,0,1,0,0,0,0,0,0,0,0,0,0}} },
00190 { "C7h", {0, 1, 0, {0, 0, 0, 0, 0, 0, 0, 1, 0}, {0,0,0,0,1,0,0,0,0,0,0,0,0,0}} },
00191 { "C8h", {1, 1, 0, {0, 0, 1, 0, 1, 0, 0, 0, 1}, {0,1,0,0,0,1,0,0,0,0,0,0,0,0}} },
00192 { "D2h", {1, 3, 0, {0, 0, 3, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00193 { "D3h", {0, 4, 0, {0, 0, 3, 1, 0, 0, 0, 0, 0}, {1,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00194 { "D4h", {1, 5, 0, {0, 0, 5, 0, 1, 0, 0, 0, 0}, {0,1,0,0,0,0,0,0,0,0,0,0,0,0}} },
00195 { "D5h", {0, 6, 0, {0, 0, 5, 0, 0, 1, 0, 0, 0}, {0,0,1,0,0,0,0,0,0,0,0,0,0,0}} },
00196 { "D6h", {1, 7, 0, {0, 0, 7, 1, 0, 0, 1, 0, 0}, {1,0,0,1,0,0,0,0,0,0,0,0,0,0}} },
00197 { "D7h", {0, 8, 0, {0, 0, 7, 0, 0, 0, 0, 1, 0}, {0,0,0,0,1,0,0,0,0,0,0,0,0,0}} },
00198 { "D8h", {1, 9, 0, {0, 0, 9, 0, 1, 0, 0, 0, 1}, {0,1,0,0,0,1,0,0,0,0,0,0,0,0}} },
00199 { "D2d", {0, 2, 0, {0, 0, 3, 0, 0, 0, 0, 0, 0}, {0,1,0,0,0,0,0,0,0,0,0,0,0,0}} },
00200 { "D3d", {1, 3, 0, {0, 0, 3, 1, 0, 0, 0, 0, 0}, {0,0,0,1,0,0,0,0,0,0,0,0,0,0}} },
00201 { "D4d", {0, 4, 0, {0, 0, 5, 0, 1, 0, 0, 0, 0}, {0,0,0,0,0,1,0,0,0,0,0,0,0,0}} },
00202 { "D5d", {1, 5, 0, {0, 0, 5, 0, 0, 1, 0, 0, 0}, {0,0,0,0,0,0,0,1,0,0,0,0,0,0}} },
00203 { "D6d", {0, 6, 0, {0, 0, 7, 1, 0, 0, 1, 0, 0}, {0,1,0,0,0,0,0,0,0,1,0,0,0,0}} },
00204 { "D7d", {1, 7, 0, {0, 0, 7, 0, 0, 0, 0, 1, 0}, {0,0,0,0,0,0,0,0,0,0,0,1,0,0}} },
00205 { "D8d", {0, 8, 0, {0, 0, 9, 0, 1, 0, 0, 0, 1}, {0,0,0,0,0,0,0,0,0,0,0,0,0,1}} },
00206 { "S4", {0, 0, 0, {0, 0, 1, 0, 0, 0, 0, 0, 0}, {0,1,0,0,0,0,0,0,0,0,0,0,0,0}} },
00207 { "S6", {1, 0, 0, {0, 0, 0, 1, 0, 0, 0, 0, 0}, {0,0,0,1,0,0,0,0,0,0,0,0,0,0}} },
00208 { "S8", {0, 0, 0, {0, 0, 1, 0, 1, 0, 0, 0, 0}, {0,0,0,0,0,1,0,0,0,0,0,0,0,0}} },
00209 { "T", {0, 0, 0, {0, 0, 3, 4, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00210 { "Th", {1, 3, 0, {0, 0, 3, 4, 0, 0, 0, 0, 0}, {0,0,0,4,0,0,0,0,0,0,0,0,0,0}} },
00211 { "Td", {0, 6, 0, {0, 0, 3, 4, 0, 0, 0, 0, 0}, {0,3,0,0,0,0,0,0,0,0,0,0,0,0}} },
00212 { "O", {0, 0, 0, {0, 0, 9, 4, 3, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00213 { "Oh", {1, 9, 0, {0, 0, 9, 4, 3, 0, 0, 0, 0}, {0,3,0,4,0,0,0,0,0,0,0,0,0,0}} },
00214 { "I", {0, 0, 0, {0, 0,15,10, 0, 6, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00215 { "Ih", {1,15, 0, {0, 0,15,10, 0, 6, 0, 0, 0}, {0,0,0,10,0,0,0,6,0,0,0,0,0,0}} },
00216 { "Cinfv", {0, 1, 1, {0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00217 { "Dinfh", {1, 2, 1, {0, 0, 1, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00218 { "Kh", {1, 1, 1, {0, 0, 0, 0, 0, 0, 0, 0, 0}, {0,0,0,0,0,0,0,0,0,0,0,0,0,0}} },
00219 };
00220
00221 #define NUMPOINTGROUPS (sizeof(pgdefinitions)/sizeof(PointGroupDefinition))
00222
00223
00224 Symmetry::Symmetry(AtomSel *mysel, MoleculeList *mymlist, int verbosity) :
00225 sel(mysel),
00226 mlist(mymlist),
00227 verbose(verbosity)
00228 {
00229 sigma = 0.1f;
00230 sphericaltop = 0;
00231 symmetrictop = 0;
00232 linear = 0;
00233 planar = 0;
00234 inversion = 0;
00235 maxaxisorder = 0;
00236 maxrotreflorder = 0;
00237 maxweight = 0;
00238 maxoverlap = 0.0;
00239 pointgroup = POINTGROUP_UNKNOWN;
00240 pointgrouporder = 0;
00241 numdihedralplanes = 0;
00242 numverticalplanes = 0;
00243 horizontalplane = -1;
00244 maxnatoms = 200;
00245 rmsd = 0.0;
00246 idealcoor = NULL;
00247 coor = NULL;
00248 checkbonds = 0;
00249 bondsum = NULL;
00250 bondsperatom = NULL;
00251 atomtype = NULL;
00252 atomindex = NULL;
00253 uniqueatoms = NULL;
00254 elementsummarystring = NULL;
00255 missingelementstring = new char[2 + 10L*(3+MAXORDERCN+2*MAXORDERCN)];
00256 additionalelementstring = new char[2 + 10L*(3+MAXORDERCN+2*MAXORDERCN)];
00257 missingelementstring[0] = '\0';
00258 additionalelementstring[0] = '\0';
00259
00260 memset(&elementsummary, 0, sizeof(ElementSummary));
00261
00262 if (sel->selected) {
00263 coor = new float[3L*sel->selected];
00264 idealcoor = new float[3L*sel->selected];
00265 bondsum = new float[3L*sel->selected];
00266 bondsperatom = new Bondlist[sel->selected];
00267 atomtype = new int[sel->selected];
00268 atomindex = new int[sel->selected];
00269 uniqueatoms = new int[sel->selected];
00270
00271 memset(bondsperatom, 0, sel->selected*sizeof(Bondlist));
00272 }
00273
00274 memset(&(inertiaaxes[0][0]), 0, 9L*sizeof(float));
00275 memset(&(inertiaeigenval[0]), 0, 3L*sizeof(float));
00276 memset(&(uniqueprimary[0]), 0, 3L*sizeof(int));
00277 memset(&(rcom[0]), 0, 3L*sizeof(float));
00278
00279
00280
00281 assign_atoms(sel, mlist, coor, atomtype);
00282
00283
00284 assign_bonds();
00285
00286
00287
00288 collintol = cosf(float(DEGTORAD(10.0f)));
00289
00290
00291 orthogtol = cosf(float(DEGTORAD(80.0f)));
00292
00293 if (sel->selected <= 10) {
00294 collintol = cosf(float(DEGTORAD(15.0f)));
00295 orthogtol = cosf(float(DEGTORAD(75.0f)));
00296 }
00297 }
00298
00299
00300
00301 Symmetry::~Symmetry(void) {
00302 if (coor) delete [] coor;
00303 if (idealcoor) delete [] idealcoor;
00304 if (bondsum) delete [] bondsum;
00305 if (atomtype) delete [] atomtype;
00306 if (atomindex) delete [] atomindex;
00307 if (uniqueatoms) delete [] uniqueatoms;
00308 if (elementsummarystring) delete [] elementsummarystring;
00309 delete [] missingelementstring;
00310 delete [] additionalelementstring;
00311 if (bondsperatom) {
00312 int i;
00313 for (i=0; i<sel->selected; i++) {
00314 if (bondsperatom[i].bondto) delete [] bondsperatom[i].bondto;
00315 if (bondsperatom[i].length) delete [] bondsperatom[i].length;
00316 }
00317 delete [] bondsperatom;
00318 }
00319 }
00320
00323 void Symmetry::get_pointgroup(char pg[8], int *order) {
00324 char n[3];
00325 if (order==NULL) sprintf(n, "%i", pointgrouporder);
00326 else {
00327 strcpy(n, "n");
00328 *order = pointgrouporder;
00329 }
00330
00331 switch (pointgroup) {
00332 case POINTGROUP_UNKNOWN:
00333 strcpy(pg, "Unknown");
00334 break;
00335 case POINTGROUP_C1:
00336 strcpy(pg, "C1");
00337 break;
00338 case POINTGROUP_CN:
00339 sprintf(pg, "C%s", n);
00340 break;
00341 case POINTGROUP_CNV:
00342 sprintf(pg, "C%sv", n);
00343 break;
00344 case POINTGROUP_CNH:
00345 sprintf(pg, "C%sh", n);
00346 break;
00347 case POINTGROUP_CINFV:
00348 strcpy(pg, "Cinfv");
00349 break;
00350 case POINTGROUP_DN:
00351 sprintf(pg, "D%s", n);
00352 break;
00353 case POINTGROUP_DND:
00354 sprintf(pg, "D%sd", n);
00355 break;
00356 case POINTGROUP_DNH:
00357 sprintf(pg, "D%sh", n);
00358 break;
00359 case POINTGROUP_DINFH:
00360 strcpy(pg, "Dinfh");
00361 break;
00362 case POINTGROUP_CI:
00363 strcpy(pg, "Ci");
00364 break;
00365 case POINTGROUP_CS:
00366 strcpy(pg, "Cs");
00367 break;
00368 case POINTGROUP_S2N:
00369 if (!strcmp(n, "n"))
00370 strcpy(pg, "S2n");
00371 else
00372 sprintf(pg, "S%i", 2*pointgrouporder);
00373 break;
00374 case POINTGROUP_T:
00375 strcpy(pg, "T");
00376 break;
00377 case POINTGROUP_TD:
00378 strcpy(pg, "Td");
00379 break;
00380 case POINTGROUP_TH:
00381 strcpy(pg, "Th");
00382 break;
00383 case POINTGROUP_O:
00384 strcpy(pg, "O");
00385 break;
00386 case POINTGROUP_OH:
00387 strcpy(pg, "Oh");
00388 break;
00389 case POINTGROUP_I:
00390 strcpy(pg, "I");
00391 break;
00392 case POINTGROUP_IH:
00393 strcpy(pg, "Ih");
00394 break;
00395 case POINTGROUP_KH:
00396 strcpy(pg, "Kh");
00397 break;
00398 }
00399 }
00400
00401
00403 int Symmetry::get_axisorder(int n) {
00404 if (n<numaxes()) return axes[n].order;
00405 return 0;
00406 }
00407
00409 int Symmetry::get_rotreflectorder(int n) {
00410 if (n<numrotreflect()) return rotreflections[n].order;
00411 return 0;
00412 }
00413
00415 int Symmetry::numS2n() {
00416 int i=0, count=0;
00417 for (i=0; i<numrotreflect(); i++) {
00418 if (rotreflections[i].order % 2 == 0) count++;
00419 }
00420 return count;
00421 }
00422
00424 int Symmetry::numprimaryaxes() {
00425 int i;
00426 int numprimary = 0;
00427 for (i=0; i<numaxes(); i++) {
00428 if (axes[i].order==maxaxisorder) numprimary++;
00429 }
00430 return numprimary;
00431 }
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441 void Symmetry::impose(int have_inversion,
00442 int nplanes, const float *planev,
00443 int naxes, const float *axisv, const int *axisorder,
00444 int nrotrefl, const float* rotreflv,
00445 const int *rotreflorder) {
00446 int i;
00447 char buf[256] = { 0 };
00448 if (have_inversion) {
00449 inversion = 1;
00450 if (verbose>1) {
00451 sprintf(buf, "imposing inversion\n");
00452 msgInfo << buf << sendmsg;
00453 }
00454 }
00455 planes.clear();
00456 for (i=0; i<nplanes; i++) {
00457 Plane p;
00458 vec_copy(p.v, &planev[3L*i]);
00459 vec_normalize(p.v);
00460 if (norm(p.v)==0.f) continue;
00461 p.overlap = 1.f;
00462 p.weight = 1;
00463 p.type = 0;
00464 planes.append(p);
00465 if (verbose>1) {
00466 sprintf(buf, "imposing plane {% .2f % .2f % .2f}\n", p.v[0], p.v[1], p.v[2]);
00467 msgInfo << buf << sendmsg;
00468 }
00469 }
00470 axes.clear();
00471 for (i=0; i<naxes; i++) {
00472 if (axisorder[i]<=1) continue;
00473 Axis a;
00474 vec_copy(a.v, &axisv[3L*i]);
00475 vec_normalize(a.v);
00476 if (norm(a.v)==0.f) continue;
00477 a.order = axisorder[i];
00478 a.overlap = 1.f;
00479 a.weight = 1;
00480 a.type = 0;
00481 axes.append(a);
00482 if (verbose>1) {
00483 sprintf(buf, "imposing axis {% .2f % .2f % .2f}\n", a.v[0], a.v[1], a.v[2]);
00484 msgInfo << buf << sendmsg;
00485 }
00486 }
00487 rotreflections.clear();
00488 for (i=0; i<nrotrefl; i++) {
00489 if (rotreflorder[i]<=1) continue;
00490 Axis a;
00491 vec_copy(a.v, &rotreflv[3L*i]);
00492 vec_normalize(a.v);
00493 if (norm(a.v)==0.f) continue;
00494 a.order = rotreflorder[i];
00495 a.overlap = 1.f;
00496 a.weight = 1;
00497 a.type = 0;
00498 rotreflections.append(a);
00499 if (verbose>1) {
00500 sprintf(buf, "imposing rraxis {% .2f % .2f % .2f}\n", a.v[0], a.v[1], a.v[2]);
00501 msgInfo << buf << sendmsg;
00502 }
00503 }
00504
00505
00506 float itensor[4][4];
00507 int ret = measure_inertia(sel, mlist, coor, rcom, inertiaaxes,
00508 itensor, inertiaeigenval);
00509 if (ret < 0) msgErr << "measure inertia failed with code " << ret << sendmsg;
00510
00511
00512
00513
00514 assign_bondvectors();
00515
00516
00517 for (i=0; i<2; i++) {
00518
00519
00520
00521 idealize_coordinates();
00522
00523
00524 memcpy(coor, idealcoor, 3L*sel->selected*sizeof(float));
00525 }
00526 }
00527
00528
00531 int Symmetry::guess(float mysigma) {
00532 if (!sel) return MEASURE_ERR_NOSEL;
00533 if (sel->selected == 0) return MEASURE_ERR_NOATOMS;
00534
00535 if (sel->selected == 1) {
00536 pointgroup = POINTGROUP_KH;
00537 elementsummarystring = new char[1];
00538 elementsummarystring[0] = '\0';
00539 uniqueatoms[0] = 1;
00540 memcpy(idealcoor, coor, 3L*sel->selected*sizeof(float));
00541
00542 return MEASURE_NOERR;
00543 }
00544
00545
00546 float maxsigma = mysigma;
00547
00548 wkf_timerhandle timer = wkf_timer_create();
00549 wkf_timer_start(timer);
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574 int beststep = -1;
00575 Best best;
00576 best.pointgroup = POINTGROUP_UNKNOWN;
00577 best.pointgrouporder = 0;
00578 best.idealcoor = NULL;
00579 int step, nstep;
00580 float stepsize = 0.05f;
00581 nstep = int(0.5f+(float)maxsigma/stepsize);
00582
00583 for (step=0; step<nstep; step++) {
00584 sigma = (1+step)*stepsize;
00585 if (verbose>1) {
00586 msgInfo << sendmsg;
00587 msgInfo << " STEP " << step << " sigma=" << sigma << sendmsg
00588 << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << sendmsg << sendmsg;
00589 }
00590
00591
00592
00593 if (!iterate_element_search()) continue;
00594
00595 wkf_timer_stop(timer);
00596 if (verbose>0) {
00597 char buf[256] = { 0 };
00598 sprintf(buf, "Guess %i completed in %f seconds, RMSD=%.3f", step, wkf_timer_time(timer), rmsd);
00599 msgInfo << buf << sendmsg;
00600 }
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610 if (!ideal_coordinate_sanity()) {
00611 if (verbose>1) {
00612 msgInfo << "INSANE IDEAL COORDINATES" << sendmsg;
00613 }
00614 break;
00615 }
00616
00617 if (verbose>1) {
00618 print_statistics();
00619 }
00620
00621
00622 determine_pointgroup();
00623
00624 char pgname[8];
00625 get_pointgroup(pgname, NULL);
00626 compare_element_summary(pgname);
00627
00628 if (verbose>1) {
00629 msgInfo << "Point group: " << pgname << sendmsg << sendmsg;
00630 print_element_summary(pgname);
00631
00632
00633 }
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 if (pointgroup!=best.pointgroup &&
00644 !strlen(additionalelementstring) &&
00645 !strlen(missingelementstring) &&
00646 pointgroup_rank(pointgroup, pointgrouporder) >
00647 pointgroup_rank(best.pointgroup, best.pointgrouporder)) {
00648 beststep = step;
00649 best.pointgroup = pointgroup;
00650 best.pointgrouporder = pointgrouporder;
00651 if (!best.idealcoor) best.idealcoor = new float[3L*sel->selected];
00652 memcpy(best.idealcoor, idealcoor, 3L*sel->selected*sizeof(float));
00653 best.linear = linear;
00654 best.planar = planar;
00655 best.inversion = inversion;
00656 best.sphericaltop = sphericaltop;
00657 }
00658 }
00659
00660
00661
00662 sigma = (1+beststep)*stepsize;
00663 if (best.idealcoor) {
00664 memcpy(coor, best.idealcoor, 3L*sel->selected*sizeof(float));
00665 delete [] best.idealcoor;
00666 }
00667 iterate_element_search();
00668
00669
00670 determine_pointgroup();
00671
00672 if (verbose>0) {
00673 msgInfo << sendmsg
00674 << "RESULT:" << sendmsg
00675 << "=======" << sendmsg << sendmsg;
00676
00677 msgInfo << "Needed tol = " << sigma << " to detect symmetry."
00678 << sendmsg;
00679
00680 print_statistics();
00681 char pgname[8];
00682 get_pointgroup(pgname, NULL);
00683 compare_element_summary(pgname);
00684 msgInfo << "Point group: " << pgname << sendmsg << sendmsg;
00685 print_element_summary(pgname);
00686 }
00687
00688
00689
00690 sort_planes();
00691
00692
00693 unique_coordinates();
00694
00695
00696
00697 orient_molecule();
00698
00699
00700 wkf_timer_destroy(timer);
00701
00702 return MEASURE_NOERR;
00703 }
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729 int Symmetry::iterate_element_search() {
00730 int iteration;
00731 char oldelementsummarystring[2 + 10*(3+MAXORDERCN+2*MAXORDERCN)];
00732 oldelementsummarystring[0] = '\0';
00733 float oldrmsd = 999999.9f;
00734 int converged = 0;
00735
00736 for (iteration=0; iteration<20; iteration++) {
00737
00738 pointgroup = POINTGROUP_UNKNOWN;
00739 inversion = 0;
00740 linear = 0;
00741 planar = 0;
00742 sphericaltop = 0;
00743 symmetrictop = 0;
00744 axes.clear();
00745 planes.clear();
00746 rotreflections.clear();
00747 bonds.clear();
00748 numdihedralplanes = 0;
00749 numverticalplanes = 0;
00750 horizontalplane = -1;
00751
00752 if (verbose>2) {
00753 msgInfo << sendmsg
00754 << "###################" << sendmsg;
00755 msgInfo << "Iteration " << iteration << sendmsg;
00756 msgInfo << "###################" << sendmsg << sendmsg;
00757 }
00758
00759
00760 assign_bondvectors();
00761
00762
00763 float itensor[4][4];
00764 int ret = measure_inertia(sel, mlist, coor, rcom, inertiaaxes,
00765 itensor, inertiaeigenval);
00766 if (ret < 0) msgErr << "measure inertia failed with code " << ret << sendmsg;
00767
00768
00769 if (verbose>2) {
00770 char buf[256] = { 0 };
00771 sprintf(buf, "eigenvalues: %.2f %.2f %.2f", inertiaeigenval[0], inertiaeigenval[1], inertiaeigenval[2]);
00772 msgInfo << buf << sendmsg;
00773 }
00774
00775
00776
00777 find_symmetry_elements();
00778
00779 if (verbose>2) {
00780 if (linear)
00781 msgInfo << "Linear molecule" << sendmsg;
00782 if (planar)
00783 msgInfo << "Planar molecule" << sendmsg;
00784 if (sphericaltop)
00785 msgInfo << "Molecule is a spherical top." << sendmsg;
00786 if (symmetrictop)
00787 msgInfo << "Molecule is a symmetric top." << sendmsg;
00788
00789 int i, numuniqueprimary = 0;
00790 for (i=0; i<3; i++) {
00791 if (uniqueprimary[i]) numuniqueprimary++;
00792 }
00793 msgInfo << "Number of unique primary axes = " << numuniqueprimary << sendmsg;
00794
00795 }
00796
00797
00798
00799
00800
00801
00802
00803
00804 idealize_elements();
00805
00806
00807
00808
00809 idealize_coordinates();
00810
00811
00812
00813
00814
00815
00816 rmsd = ideal_coordinate_rmsd();
00817
00818
00819 build_element_summary();
00820
00821
00822
00823 build_element_summary_string(elementsummary, elementsummarystring);
00824
00825
00826
00827
00828 converged = 0;
00829 if (fabs(rmsd-oldrmsd)<0.01 &&
00830 !strcmp(elementsummarystring, oldelementsummarystring)) {
00831 if (verbose>1) {
00832 msgInfo << "Symmetry search converged after " << iteration+1
00833 << " iterations" << sendmsg;
00834 }
00835 converged = 1;
00836 }
00837
00838
00839 if (verbose>3 && (converged || verbose>2)) {
00840 print_statistics();
00841 }
00842 if (converged) break;
00843
00844
00845 memcpy(coor, idealcoor, 3L*sel->selected*sizeof(float));
00846
00847 oldrmsd = rmsd;
00848 strcpy(oldelementsummarystring, elementsummarystring);
00849 }
00850
00851 if (!converged) return 0;
00852 return 1;
00853 }
00854
00855
00856
00857 void Symmetry::find_symmetry_elements() {
00858
00859
00860
00861
00862 find_elements_from_inertia();
00863
00864
00865 find_planes();
00866
00867
00868 purge_planes(0.5);
00869
00870
00871 classify_planes();
00872
00873
00874 check_add_inversion();
00875
00876
00877
00878
00879
00880 find_axes_from_plane_intersections();
00881
00882
00883 assign_axis_orders();
00884
00885
00886
00887
00888 find_C2axes();
00889
00890 assign_prelimaxis_orders(2);
00891
00892 if (!axes.num() && !planes.num()) {
00893 int j;
00894 for (j=3; j<=8; j++) {
00895 if (isprime(j) && !axes.num()) {
00896 find_axes(j);
00897 assign_prelimaxis_orders(j);
00898 }
00899 }
00900 }
00901
00902
00903 sort_axes();
00904
00905
00906
00907
00908
00909
00910
00911 int i;
00912 for (i=0; i<numaxes(); i++) {
00913 check_add_rotary_reflection(axes[i].v, 2*axes[i].order);
00914 }
00915
00916
00917 maxrotreflorder = 0;
00918 for (i=0; i<numrotreflect(); i++) {
00919 vec_normalize(rotreflections[i].v);
00920
00921 if (rotreflections[i].order>maxrotreflorder)
00922 maxrotreflorder = rotreflections[i].order;
00923 }
00924
00925
00926
00927 purge_axes(0.7f);
00928 purge_rotreflections(0.75f);
00929
00930
00931
00932 classify_planes();
00933
00934
00935 classify_axes();
00936 }
00937
00938
00939
00940
00941
00942
00943
00944
00945 void Symmetry::find_elements_from_inertia() {
00946 int i;
00947 int numuniqueprimary = 0;
00948 memset(uniqueprimary, 0, 3L*sizeof(int));
00949
00950
00951 float eigenval[3];
00952 float e0 = inertiaeigenval[0];
00953 for (i=0; i<3; i++) {
00954 eigenval[i] = inertiaeigenval[i]/e0;
00955 }
00956
00957
00958 if (fabs(eigenval[0]-eigenval[1])<0.05 && eigenval[2]<0.05) {
00959 linear = 1;
00960 }
00961
00962
00963 float overlap[3];
00964 int order[3], primaryCn[3];
00965 for (i=0; i<3; i++) {
00966 order[i] = axis_order(inertiaaxes[i], overlap+i);
00967
00968 if ((order[i]>1 && overlap[i]>1.5*OVERLAPCUTOFF) ||
00969 (linear && eigenval[i]<0.05)) {
00970 primaryCn[i] = 1;
00971 } else {
00972 primaryCn[i] = 0;
00973 }
00974 }
00975
00976 float tol = 0.25f*sigma + 1.0f/(powf(float(sel->selected+1), 1.5f));
00977
00978
00979
00980 if (fabs(eigenval[0]-eigenval[1])<tol &&
00981 fabs(eigenval[1]-eigenval[2])<tol) {
00982
00983 sphericaltop = 1;
00984
00985 } else {
00986
00987
00988
00989 if (fabs(eigenval[0]-eigenval[1])<tol ||
00990 fabs(eigenval[1]-eigenval[2])<tol ||
00991 fabs(eigenval[0]-eigenval[2])<tol) {
00992
00993
00994
00995
00996
00997 if (fabs(eigenval[1]-eigenval[2])<tol && primaryCn[0]) {
00998 uniqueprimary[0] = 1;
00999 numuniqueprimary++;
01000 }
01001 if (fabs(eigenval[0]-eigenval[2])<tol && primaryCn[1]) {
01002 uniqueprimary[1] = 1;
01003 numuniqueprimary++;
01004 }
01005 if (fabs(eigenval[0]-eigenval[1])<tol && primaryCn[2]) {
01006 uniqueprimary[2] = 1;
01007 numuniqueprimary++;
01008 }
01009
01010 if (numuniqueprimary==1) {
01011 symmetrictop = 1;
01012 }
01013 else {
01014 uniqueprimary[0] = 1;
01015 uniqueprimary[1] = 1;
01016 uniqueprimary[2] = 1;
01017 numuniqueprimary = 3;
01018 }
01019 }
01020 else {
01021 uniqueprimary[0] = 1;
01022 uniqueprimary[1] = 1;
01023 uniqueprimary[2] = 1;
01024 numuniqueprimary = 3;
01025 }
01026 }
01027
01028
01029 for (i=0; i<3; i++) {
01030
01031
01032
01033 int planarmol = is_planar(inertiaaxes[i]);
01034 if (planarmol && !linear) {
01035 if (verbose>3) {
01036 msgInfo << "Planar mol: primary axes " << i << " defines plane" << sendmsg;
01037 }
01038 check_add_plane(inertiaaxes[i]);
01039 planar = 1;
01040 }
01041
01042
01043 if (!uniqueprimary[i]) continue;
01044
01045 if (linear) {
01046
01047
01048
01049
01050 Axis a;
01051 vec_copy(a.v, inertiaaxes[i]);
01052 a.order = INFINITE_ORDER;
01053 float overlap1 = score_axis(a.v, 2);
01054 float overlap2 = score_axis(a.v, 4);
01055 a.overlap = 0.5f*(overlap1+overlap2);
01056 a.weight = sel->num_atoms/2;
01057 a.type = 0;
01058 axes.append(a);
01059
01060
01061 Plane p;
01062 vec_copy(p.v, inertiaaxes[0]);
01063 p.overlap = eigenval[0]*eigenval[1] - eigenval[2];
01064 p.weight = sel->num_atoms/2;
01065 p.type = 0;
01066 planes.append(p);
01067
01068 } else {
01069
01070
01071 if (order[i] > 1 && overlap[i]>OVERLAPCUTOFF) {
01072
01073 Axis a;
01074 vec_copy(a.v, inertiaaxes[i]);
01075 a.order = order[i];
01076 a.overlap = overlap[i];
01077 a.weight = 1;
01078 a.type = 0;
01079 axes.append(a);
01080 }
01081
01082
01083
01084
01085
01086
01087 int weight = 1;
01088 if (planarmol) weight = sel->num_atoms/2;
01089
01090 check_add_plane(inertiaaxes[i], weight);
01091 }
01092 }
01093 }
01094
01095
01096
01097
01098 void Symmetry::find_planes() {
01099 int i,j;
01100 float posA[3], posB[3];
01101
01102
01103 for (i=0; i<sel->selected; i++) {
01104
01105
01106
01107 if (sel->selected > maxnatoms && vmd_random() % sel->selected > maxnatoms)
01108 continue;
01109
01110 vec_sub(posA, coor+3L*i, rcom);
01111
01112 float rA = sqrtf(posA[0]*posA[0] + posA[1]*posA[1] + posA[2]*posA[2]);
01113
01114
01115
01116
01117 if (rA < sigma || rA < sel->num_atoms/7.0*sigma) continue;
01118
01119 for (j=0; j<i; j++) {
01120 vec_sub(posB, coor+3L*j, rcom);
01121
01122 float rB = sqrtf(posB[0]*posB[0] + posB[1]*posB[1] + posB[2]*posB[2]);
01123
01124 if (fabs(rA-rB) > sigma) continue;
01125
01126
01127 if (atomtype[i]==atomtype[j]) {
01128
01129
01130
01131 float dist = distance(posA, posB);
01132 if (dist<0.25) continue;
01133
01134
01135
01136
01137
01138
01139 float normal[3];
01140 vec_sub(normal, posA, posB);
01141 vec_normalize(normal);
01142
01143
01144 check_add_plane(normal);
01145 }
01146 }
01147 }
01148
01149
01150 for (i=0; i<numplanes(); i++) {
01151 vec_normalize(planes[i].v);
01152
01153 }
01154 }
01155
01156
01157
01158
01159 void Symmetry::find_C2axes() {
01160 int i,j;
01161 float posA[3], posB[3];
01162 float rA, rB;
01163
01164 for (i=0; i<sel->selected; i++) {
01165
01166
01167 if (sel->selected > maxnatoms && vmd_random() % sel->selected > maxnatoms)
01168 continue;
01169
01170 vec_sub(posA, coor+3L*i, rcom);
01171
01172 rA = sqrtf(posA[0]*posA[0] + posA[1]*posA[1] + posA[2]*posA[2]);
01173
01174
01175
01176
01177 if (rA < sigma || rA < sel->num_atoms/7.0*sigma) continue;
01178
01179 for (j=i+1; j<sel->selected; j++) {
01180 vec_sub(posB, coor+3L*j, rcom);
01181
01182 rB = sqrtf(posB[0]*posB[0] + posB[1]*posB[1] + posB[2]*posB[2]);
01183
01184 if (fabs(rA-rB) > sigma) continue;
01185
01186
01187 if (atomtype[i]==atomtype[j]) {
01188
01189 float alpha = angle(posA, posB);
01190
01191
01192
01193
01194
01195 if (fabs(fabs(alpha)-180) > 10) {
01196 float testaxis[3];
01197 vec_add(testaxis, posA, posB);
01198 vec_normalize(testaxis);
01199
01200
01201 if (!planes.num() || numverticalplanes) {
01202
01203 check_add_axis(testaxis, 2);
01204 }
01205 }
01206 }
01207 }
01208 }
01209
01210
01211 for (i=0; i<numaxes(); i++) {
01212 vec_normalize(axes[i].v);
01213 }
01214 }
01215
01216
01217 static int fac(int n) {
01218 if (n==0) return 1;
01219 int i, x=1;
01220 for (i=1; i<=n; i++) x*=i;
01221 return x;
01222 }
01223
01224
01225
01226
01227 class Combinations {
01228 public:
01229 int *combolist;
01230 int *combo;
01231 int num;
01232 int n;
01233 int k;
01234
01235 Combinations(int N, int K);
01236 ~Combinations();
01237
01238
01239 void create() { recurse(0, -1); }
01240
01241 void recurse(int begin, int level);
01242 void print();
01243
01244
01245 int* get(int i) { return &combolist[i*k]; }
01246 };
01247
01248 Combinations::Combinations(int N, int K) : n(N), k(K) {
01249 if (n>10) n = 10;
01250 combo = new int[k];
01251 combolist = new int[k*fac(n)/(fac(k)*fac(n-k))];
01252 num = 0;
01253 }
01254
01255 Combinations::~Combinations() {
01256 delete [] combo;
01257 delete [] combolist;
01258 }
01259
01260
01261 void Combinations::recurse(int begin, int level) {
01262 int i;
01263 level++;
01264 if (level>=k) {
01265 for (i=0; i<k; i++) {
01266 combolist[num*k+i] = combo[i];
01267 }
01268 num++;
01269 return;
01270 }
01271
01272 for (i=begin; i<n; i++) {
01273 combo[level] = i;
01274 recurse(i+1, level);
01275 }
01276 }
01277
01278
01279 void Combinations::print() {
01280 int i, j;
01281 for (i=0; i<fac(n)/(fac(k)*fac(n-k)); i++) {
01282 printf("combo %d/%d {", i, num);
01283 for (j=0; j<k; j++) {
01284 printf(" %d", combolist[i*k+j]);
01285 }
01286 printf("}\n");
01287 }
01288
01289 }
01290
01291
01292
01293
01294
01295 void Symmetry::find_axes(int order) {
01296 int i,j;
01297 float posA[3], posB[3];
01298 float rA, rB;
01299 int *atomtuple = new int[sel->selected];
01300
01301
01302 for (i=0; i<sel->selected; i++) {
01303
01304
01305 if (sel->selected > maxnatoms && vmd_random() % sel->selected > maxnatoms)
01306 continue;
01307
01308 vec_sub(posA, coor+3L*i, rcom);
01309
01310 rA = sqrtf(posA[0]*posA[0] + posA[1]*posA[1] + posA[2]*posA[2]);
01311
01312
01313
01314
01315 if (rA < sigma || rA < sel->num_atoms/7.0*sigma) continue;
01316
01317 atomtuple[0] = i;
01318 int ntup = 1;
01319
01320 for (j=i+1; j<sel->selected; j++) {
01321
01322
01323
01324 if (sel->selected > maxnatoms && vmd_random() % sel->selected > maxnatoms)
01325 continue;
01326
01327
01328 if (atomtype[j]!=atomtype[i]) continue;
01329
01330 vec_sub(posB, coor+3L*j, rcom);
01331
01332 rB = sqrtf(posB[0]*posB[0] + posB[1]*posB[1] + posB[2]*posB[2]);
01333
01334
01335 if (fabs(rA-rB) > sigma) continue;
01336
01337 atomtuple[ntup++] = j;
01338
01339
01340
01341 if (ntup>=10) break;
01342 }
01343 if (ntup<order) continue;
01344
01345
01346
01347
01348
01349
01350
01351 Combinations combo(ntup, order);
01352 combo.create();
01353
01354
01355 int m;
01356 for (j=0; j<combo.num; j++) {
01357 float testaxis[3], pos[3], pos2[3], cross[3], normal[3];
01358 vec_zero(testaxis);
01359 vec_zero(normal);
01360
01361
01362
01363
01364 int inplane = 1, anyinplane = 0;
01365 for (m=0; m<order; m++) {
01366 int atom = atomtuple[combo.get(j)[m]];
01367
01368 vec_sub(pos, coor+3L*atom, rcom);
01369 vec_add(testaxis, testaxis, pos);
01370
01371
01372
01373
01374
01375 int n, found = 0;
01376 for (n=0; n<order; n++) {
01377 if (n==m) continue;
01378 int atom2 = atomtuple[combo.get(j)[m]];
01379 vec_sub(pos2, coor+3L*atom2, rcom);
01380 if (angle(pos, pos2)>45.f) {
01381 found = 1;
01382 break;
01383 }
01384 }
01385 if (!found) continue;
01386
01387 cross_prod(cross, pos2, pos);
01388
01389
01390
01391
01392 if (collinear(normal, cross, collintol)) {
01393 vec_add(normal, normal, cross);
01394 anyinplane = 1;
01395 } else {
01396
01397 inplane = 0;
01398 break;
01399 }
01400 }
01401
01402
01403
01404
01405
01406
01407 if (!inplane || !anyinplane) continue;
01408
01409 vec_normalize(normal);
01410
01411
01412 printf("Checking C%d axis\n", order);
01413 check_add_axis(normal, order);
01414 }
01415
01416 }
01417
01418 delete [] atomtuple;
01419
01420
01421 for (i=0; i<numaxes(); i++) {
01422 vec_normalize(axes[i].v);
01423 }
01424 }
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434 void Symmetry::find_axes_from_plane_intersections() {
01435
01436 if (numplanes()<2) return;
01437
01438 int i,j;
01439
01440 for (i=0; i<numplanes(); i++) {
01441 for (j=i+1; j<numplanes(); j++) {
01442 float newaxis[3];
01443 cross_prod(newaxis, plane(i), plane(j));
01444
01445 vec_normalize(newaxis);
01446
01447
01448 if (norm(newaxis)<0.05) continue;
01449
01450
01451
01452
01453 int k;
01454 bool found = 0;
01455 for (k=0; k<numaxes(); k++) {
01456 float avgaxis[3];
01457 vec_copy(avgaxis, axis(k));
01458 vec_normalize(avgaxis);
01459
01460 float dot = dot_prod(avgaxis, newaxis);
01461
01462 if (fabs(dot) > collintol) {
01463
01464
01465 if (dot>0) {
01466 vec_incr(axis(k), newaxis);
01467 } else {
01468 vec_sub(axis(k), axis(k), newaxis);
01469 }
01470 axes[k].weight++;
01471 found = 1;
01472 break;
01473 }
01474 }
01475
01476 if (!found) {
01477
01478 Axis a;
01479 vec_copy(a.v, newaxis);
01480 a.type = 0;
01481 a.order = 0;
01482 a.overlap = 0.0;
01483 if (planes[i].weight<planes[j].weight)
01484 a.weight = planes[i].weight;
01485 else
01486 a.weight = planes[j].weight;
01487 axes.append(a);
01488 }
01489 }
01490 }
01491
01492
01493 for (i=0; i<numaxes(); i++) {
01494 vec_normalize(axis(i));
01495 }
01496 }
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507 void Symmetry::classify_planes() {
01508 if (!numplanes()) return;
01509 if (!numaxes()) return;
01510
01511 numdihedralplanes = 0;
01512 numverticalplanes = 0;
01513 horizontalplane = -1;
01514
01515 int i;
01516
01517
01518 for (i=0; i<numplanes(); i++) {
01519 if (collinear(axis(0), plane(i), collintol)) {
01520 horizontalplane = i;
01521 planes[i].type = HORIZONTALPLANE;
01522 break;
01523 }
01524 }
01525
01526
01527 for (i=0; i<numplanes(); i++) {
01528 if (!orthogonal(axis(0), plane(i), orthogtol)) continue;
01529
01530
01531 numverticalplanes++;
01532 planes[i].type = VERTICALPLANE;
01533
01534
01535 int j, k;
01536 for (j=1; j<numaxes(); j++) {
01537 if (axes[j].order != 2) continue;
01538 if (!orthogonal(axis(j), axis(0), orthogtol)) continue;
01539
01540 for (k=j+1; k<numaxes(); k++) {
01541 if (axes[k].order != 2) continue;
01542 if (!orthogonal(axis(k), axis(0), orthogtol)) continue;
01543
01544
01545 float bisect[3];
01546 vec_add(bisect, axis(j), axis(k));
01547 vec_normalize(bisect);
01548 if (orthogonal(bisect, plane(i), orthogtol)) {
01549 planes[i].type = DIHEDRALPLANE;
01550 numdihedralplanes++;
01551 j=numaxes();
01552 break;
01553 }
01554
01555 vec_sub(bisect, axis(j), axis(k));
01556 vec_normalize(bisect);
01557 if (orthogonal(bisect, plane(i), orthogtol)) {
01558 planes[i].type = DIHEDRALPLANE;
01559 numdihedralplanes++;
01560 j=numaxes();
01561 break;
01562 }
01563
01564 }
01565 }
01566 }
01567 }
01568
01569 void Symmetry::classify_axes() {
01570 if (!numaxes()) return;
01571
01572
01573
01574 int numprimary = 0;
01575 int i;
01576 for (i=0; i<numaxes(); i++) {
01577 if (axes[i].order==maxaxisorder) {
01578 axes[i].type = PRINCIPAL_AXIS;
01579 numprimary++;
01580 }
01581 }
01582 for (i=1; i<numaxes(); i++) {
01583 if (orthogonal(axis(0), axis(i), orthogtol))
01584 axes[i].type |= PERPENDICULAR_AXIS;
01585 }
01586 }
01587
01588
01589 float Symmetry::score_inversion() {
01590
01591 Matrix4 inv;
01592 inv.translate(rcom[0], rcom[1], rcom[2]);
01593 inv.scale(-1.0);
01594 inv.translate(-rcom[0], -rcom[1], -rcom[2]);
01595
01596
01597 return trans_overlap(atomtype, coor, sel->selected, &inv, 1.5f*sigma,
01598 NOSKIP_IDENTICAL, maxnatoms);
01599 }
01600
01601
01602
01603 void Symmetry::check_add_inversion() {
01604
01605 float overlap = score_inversion();
01606
01607
01608
01609
01610
01611
01612
01613 if (overlap>0.5-0.2/(1+(planes.num()*axes.num()))) {
01614 inversion = 1;
01615 }
01616 }
01617
01618
01619 float Symmetry::score_plane(const float *normal) {
01620
01621
01622
01623 Matrix4 mirror;
01624 mirror.translate(rcom[0], rcom[1], rcom[2]);
01625 mirror.scale(-1.0);
01626 mirror.rotate_axis(normal, float(VMD_PI));
01627 mirror.translate(-rcom[0], -rcom[1], -rcom[2]);
01628
01629
01630 return trans_overlap(atomtype, coor, sel->selected, &mirror, 1.5f*sigma, NOSKIP_IDENTICAL, maxnatoms);
01631 }
01632
01633
01634
01635
01636
01637
01638
01639 void Symmetry::check_add_plane(const float *normal, int weight) {
01640
01641
01642 float overlap = score_plane(normal);
01643
01644
01645 if (overlap<OVERLAPCUTOFF) return;
01646
01647
01648
01649
01650
01651
01652
01653
01654 int k;
01655 bool found=0;
01656 for (k=0; k<numplanes(); k++) {
01657 float avgnormal[3];
01658 vec_copy(avgnormal, plane(k));
01659 vec_normalize(avgnormal);
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669 float dot = dot_prod(avgnormal, normal);
01670 if (fabs(dot) > collintol) {
01671
01672
01673
01674 if (dot>0) {
01675 vec_incr(plane(k), normal);
01676 } else {
01677 vec_scaled_add(plane(k), -1, normal);
01678 }
01679 planes[k].overlap = (planes[k].weight*planes[k].overlap + overlap)/(planes[k].weight+1);
01680 (planes[k].weight)++;
01681 found = 1;
01682 break;
01683 }
01684 }
01685 if (!found) {
01686
01687 Plane p;
01688 vec_copy(p.v, normal);
01689 p.overlap = overlap;
01690 p.weight = weight;
01691 p.type = 0;
01692 planes.append(p);
01693 }
01694 }
01695
01696
01697
01698 float Symmetry::score_axis(const float *testaxis, int order) {
01699
01700 Matrix4 rot;
01701 rot.translate(rcom[0], rcom[1], rcom[2]);
01702 rot.rotate_axis(testaxis, float(VMD_TWOPI)/order);
01703 rot.translate(-rcom[0], -rcom[1], -rcom[2]);
01704
01705
01706 return trans_overlap(atomtype, coor, sel->selected, &rot, 2*sigma,
01707 NOSKIP_IDENTICAL, maxnatoms);
01708 }
01709
01710
01711
01712
01713
01714
01715
01716 void Symmetry::check_add_axis(const float *testaxis, int order) {
01717 int k;
01718 bool found = 0;
01719 for (k=0; k<numaxes(); k++) {
01720 float avgaxis[3];
01721 vec_copy(avgaxis, axis(k));
01722 vec_normalize(avgaxis);
01723
01724
01725 float dot = dot_prod(avgaxis, testaxis);
01726 if (fabs(dot) > collintol) {
01727
01728
01729
01730
01731 if (axes[k].order==-order) {
01732 if (dot>0) {
01733 vec_incr(axis(k), testaxis);
01734 } else {
01735 vec_scaled_add(axis(k), -1, testaxis);
01736 }
01737 axes[k].weight++;
01738 }
01739 found = 1;
01740 break;
01741 }
01742 }
01743 if (!found) {
01744
01745 float overlap = score_axis(testaxis, order);
01746 if (overlap>OVERLAPCUTOFF) {
01747 Axis a;
01748 vec_copy(a.v, testaxis);
01749
01750
01751
01752 a.order = -order;
01753 a.overlap = overlap;
01754 a.weight = 1;
01755 a.type = 0;
01756 axes.append(a);
01757 }
01758 }
01759 }
01760
01761
01762 float Symmetry::score_rotary_reflection(const float *testaxis, int order) {
01763
01764
01765
01766 Matrix4 rot;
01767 rot.translate(rcom[0], rcom[1], rcom[2]);
01768 rot.rotate_axis(testaxis, float(VMD_TWOPI)/order);
01769 rot.scale(-1.0);
01770 rot.rotate_axis(testaxis, float(VMD_PI));
01771 rot.translate(-rcom[0], -rcom[1], -rcom[2]);
01772
01773
01774 return trans_overlap(atomtype, coor, sel->selected, &rot, sigma,
01775 NOSKIP_IDENTICAL, maxnatoms);
01776 }
01777
01778
01779
01780 void Symmetry::check_add_rotary_reflection(const float *testaxis, int maxorder) {
01781 if (maxorder<4) return;
01782
01783
01784
01785 int n;
01786 for (n=3; n<=maxorder; n++) {
01787 if (n>=9 && n%2) continue;
01788 if (maxorder%n) continue;
01789
01790
01791 float overlap = score_rotary_reflection(testaxis, n);
01792
01793
01794
01795 if (overlap>OVERLAPCUTOFF) {
01796 int k;
01797 bool found = 0;
01798 for (k=0; k<numrotreflect(); k++) {
01799 float avgaxis[3];
01800 vec_copy(avgaxis, rotreflect(k));
01801 vec_normalize(avgaxis);
01802
01803 if (n!=rotreflections[k].order) continue;
01804
01805 float dot = dot_prod(avgaxis, testaxis);
01806 if (fabs(dot) > collintol) {
01807
01808
01809 if (dot>0) {
01810 vec_incr(rotreflect(k), testaxis);
01811 } else {
01812 vec_scaled_add(rotreflect(k), -1, testaxis);
01813 }
01814 rotreflections[k].weight++;
01815 found = 1;
01816 break;
01817 }
01818 }
01819 if (!found) {
01820
01821 Axis a;
01822 vec_copy(a.v, testaxis);
01823 a.order = n;
01824 a.overlap = overlap;
01825 a.weight = 1;
01826 a.type = 0;
01827 rotreflections.append(a);
01828 }
01829 }
01830 }
01831 }
01832
01833
01834
01835
01836 void Symmetry::purge_planes(float cutoff) {
01837 maxweight = 0;
01838 maxoverlap = 0.0;
01839 if (!numplanes()) return;
01840
01841 maxweight = planes[0].weight;
01842 maxoverlap = planes[0].overlap;
01843 int i;
01844 for (i=1; i<numplanes(); i++) {
01845 if (planes[i].overlap>maxoverlap) maxoverlap = planes[i].overlap;
01846 if (planes[i].weight >maxweight) maxweight = planes[i].weight;
01847 }
01848
01849 float tol = cutoff*maxoverlap;
01850 int *keep = new int[planes.num()];
01851 memset(keep, 0, planes.num()*sizeof(int));
01852
01853 for (i=0; i<numplanes(); i++) {
01854 if (planes[i].overlap>tol) {
01855
01856 keep[i] = 1;
01857 }
01858 }
01859
01860 prune_planes(keep);
01861
01862 delete [] keep;
01863 }
01864
01865
01866
01867
01868 void Symmetry::purge_axes(float cutoff) {
01869 if (!numaxes()) return;
01870
01871 int i;
01872 for (i=0; i<numaxes(); i++) {
01873 if (axes[i].overlap>maxoverlap) maxoverlap = axes[i].overlap;
01874 if (axes[i].weight >maxweight) maxweight = axes[i].weight;
01875 }
01876
01877 float tol = cutoff*maxoverlap;
01878 int *keep = new int[axes.num()];
01879 memset(keep, 0, axes.num()*sizeof(int));
01880
01881 for (i=0; i<numaxes(); i++) {
01882 if (axes[i].overlap>tol) {
01883
01884 keep[i] = 1;
01885 }
01886 }
01887
01888 prune_axes(keep);
01889
01890 delete [] keep;
01891 }
01892
01893
01894
01895
01896 void Symmetry::purge_rotreflections(float cutoff) {
01897 if (!numrotreflect()) return;
01898
01899 int i;
01900 for (i=0; i<numrotreflect(); i++) {
01901 if (rotreflections[i].overlap>maxoverlap) maxoverlap = rotreflections[i].overlap;
01902 if (rotreflections[i].weight >maxweight) maxweight = rotreflections[i].weight;
01903 }
01904
01905 float tol = cutoff*maxoverlap;
01906 int *keep = new int[rotreflections.num()];
01907 memset(keep, 0, rotreflections.num()*sizeof(int));
01908
01909 for (i=0; i<numrotreflect(); i++) {
01910 if (rotreflections[i].overlap>tol) {
01911
01912 keep[i] = 1;
01913 }
01914 else if (verbose>4) {
01915
01916 char buf[512] = { 0 };
01917 sprintf(buf, "purged S%i rotary reflection %i (weight=%i, overlap=%.2f tol=%.2f)",
01918 rotreflections[i].order, i, rotreflections[i].weight, rotreflections[i].overlap, tol);
01919 msgInfo << buf << sendmsg;
01920 }
01921 }
01922
01923 prune_rotreflections(keep);
01924
01925 delete [] keep;
01926 }
01927
01928
01929
01930 void Symmetry::prune_planes(int *keep) {
01931 if (!planes.num()) return;
01932
01933 numverticalplanes = 0;
01934 numdihedralplanes = 0;
01935
01936 int numkept = 0;
01937 Plane *keptplanes = new Plane[numplanes()];
01938 int i;
01939 for (i=0; i<planes.num(); i++) {
01940
01941 if (keep[i]) {
01942
01943 keptplanes[numkept] = planes[i];
01944 if (planes[i].type==HORIZONTALPLANE) horizontalplane=numkept;
01945 else if (planes[i].type==VERTICALPLANE) numverticalplanes++;
01946 else if (planes[i].type==DIHEDRALPLANE) {
01947 numdihedralplanes++;
01948 numverticalplanes++;
01949 }
01950 numkept++;
01951
01952 } else {
01953
01954 if (planes[i].type==HORIZONTALPLANE) horizontalplane=-1;
01955 if (verbose>3) {
01956 char buf[256] = { 0 };
01957 sprintf(buf, "removed %s%s%s plane %i (weight=%i, overlap=%.2f)\n",
01958 (planes[i].type==HORIZONTALPLANE?"horiz":""),
01959 (planes[i].type==VERTICALPLANE?"vert":""),
01960 (planes[i].type==DIHEDRALPLANE?"dihed":""),
01961 i, planes[i].weight, planes[i].overlap);
01962 msgInfo << buf << sendmsg;
01963 }
01964 }
01965 }
01966 memcpy(&(planes[0]), keptplanes, numkept*sizeof(Plane));
01967 planes.truncatelastn(numplanes()-numkept);
01968
01969 delete [] keptplanes;
01970 }
01971
01972
01973
01974 void Symmetry::prune_axes(int *keep) {
01975 if (!axes.num()) return;
01976
01977 int numkept = 0;
01978 Axis *keptaxes = new Axis[numaxes()];
01979
01980 maxaxisorder = 0;
01981 int i;
01982 for (i=0; i<axes.num(); i++) {
01983
01984 if (keep[i]) {
01985
01986 keptaxes[numkept++] = axes[i];
01987
01988 if (axes[i].order>maxaxisorder)
01989 maxaxisorder = axes[i].order;
01990 }
01991 }
01992 memcpy(&(axes[0]), keptaxes, numkept*sizeof(Axis));
01993 axes.truncatelastn(numaxes()-numkept);
01994
01995 delete [] keptaxes;
01996 }
01997
01998
01999
02000 void Symmetry::prune_rotreflections(int *keep) {
02001 if (!rotreflections.num()) return;
02002
02003 int numkept = 0;
02004 Axis *keptrotrefl = new Axis[numrotreflect()];
02005
02006 maxrotreflorder = 0;
02007 int i;
02008 for (i=0; i<numrotreflect(); i++) {
02009 if (keep[i]) {
02010
02011 keptrotrefl[numkept++] = rotreflections[i];
02012
02013 if (rotreflections[i].order>maxrotreflorder)
02014 maxrotreflorder = rotreflections[i].order;
02015 }
02016 }
02017
02018 memcpy(&(rotreflections[0]), keptrotrefl, numkept*sizeof(Axis));
02019 rotreflections.truncatelastn(numrotreflect()-numkept);
02020
02021 delete [] keptrotrefl;
02022 }
02023
02024
02025
02026 void Symmetry::assign_axis_orders() {
02027 if (!numaxes()) return;
02028
02029 maxaxisorder = axes[0].order;
02030
02031
02032 int i;
02033 for (i=0; i<numaxes(); i++) {
02034 if (!axes[i].order) axes[i].order = axis_order(axis(i), &(axes[i].overlap));
02035
02036
02037
02038 if (axes[i].order>maxaxisorder) {
02039 maxaxisorder = axes[i].order;
02040 }
02041 }
02042 }
02043
02044
02045
02046 void Symmetry::assign_prelimaxis_orders(int order) {
02047 int i;
02048 for (i=0; i<numaxes(); i++) {
02049 if (axes[i].order==-order) {
02050 axes[i].order=order;
02051 if (axes[i].weight>1) {
02052
02053
02054
02055 axes[i].overlap = score_axis(axis(i), order);
02056 }
02057 }
02058
02059 if (axes[i].order>maxaxisorder) maxaxisorder = axes[i].order;
02060 }
02061 }
02062
02063
02064
02065
02066 void Symmetry::sort_axes() {
02067 int i;
02068
02069
02070 int numsortedaxes = 0;
02071 for (i=0; i<numaxes(); i++) {
02072 if (axes[i].order>1 || axes[i].order==INFINITE_ORDER) numsortedaxes++;
02073 }
02074
02075 Axis *sortedaxes = new Axis[numsortedaxes*sizeof(Axis)];
02076 int j, k=0, have_inf=0;
02077 for (j=0; j<numaxes(); j++) {
02078 if (axes[j].order == INFINITE_ORDER) {
02079 sortedaxes[k] = axes[j];
02080 k++;
02081 have_inf = 1;
02082 }
02083 }
02084
02085 for (i=maxaxisorder; i>1; i--) {
02086 for (j=0; j<numaxes(); j++) {
02087 if (axes[j].order == i) {
02088 if (have_inf &&
02089 collinear(sortedaxes[0].v, axes[j].v, collintol)) {
02090 continue;
02091 }
02092 sortedaxes[k] = axes[j];
02093 k++;
02094 }
02095 }
02096 }
02097
02098 memcpy(&(axes[0]), sortedaxes, numsortedaxes*sizeof(Axis));
02099 axes.truncatelastn(numaxes()-numsortedaxes);
02100
02101 delete [] sortedaxes;
02102 }
02103
02104
02105
02106
02107 void Symmetry::sort_planes() {
02108 Plane *sortedplanes = new Plane[numplanes()*sizeof(Plane)];
02109
02110 int k = 0;
02111 if (horizontalplane>=0) {
02112 sortedplanes[k] = planes[horizontalplane];
02113 horizontalplane = k++;
02114 }
02115
02116 int j;
02117 for (j=0; j<numplanes(); j++) {
02118 if (planes[j].type == DIHEDRALPLANE) {
02119 sortedplanes[k++] = planes[j];
02120 }
02121 }
02122 for (j=0; j<numplanes(); j++) {
02123 if (planes[j].type == VERTICALPLANE) {
02124 sortedplanes[k++] = planes[j];
02125 }
02126 }
02127 for (j=0; j<numplanes(); j++) {
02128 if (planes[j].type == 0) {
02129 sortedplanes[k++] = planes[j];
02130 }
02131 }
02132
02133 memcpy(&(planes[0]), sortedplanes, numplanes()*sizeof(Plane));
02134
02135 delete [] sortedplanes;
02136 }
02137
02138
02139
02140
02141
02142 int Symmetry::axis_order(const float *axis, float *overlap) {
02143 int i;
02144
02145
02146 float overlaparray[MAXORDERCN+1];
02147 float overlappermarray[MAXORDERCN+1];
02148 overlappermarray[1] = 0.0;
02149 for (i=2; i<=MAXORDERCN; i++) {
02150 Matrix4 rot;
02151
02152 rot.translate(rcom[0], rcom[1], rcom[2]);
02153 rot.rotate_axis(axis, float(DEGTORAD(360.0f/i)));
02154 rot.translate(-rcom[0], -rcom[1], -rcom[2]);
02155 overlaparray[i] = trans_overlap(atomtype, coor, sel->selected, &rot, 1.5f*sigma, SKIP_IDENTICAL, maxnatoms, overlappermarray[i]);
02156 }
02157
02158
02159 float maxover = overlaparray[2];
02160
02161 for (i=3; i<=MAXORDERCN; i++) {
02162
02163 if (overlaparray[i]>maxover) {
02164 maxover = overlaparray[i];
02165 }
02166 }
02167
02168
02169
02170 float maxfalseov = 0.0f;
02171 int maxorder = 1;
02172 short int orders[MAXORDERCN+1];
02173 for (i=2; i<=MAXORDERCN; i++) {
02174 if (overlaparray[i]>0.5f*maxover) {
02175 orders[i] = 1;
02176 maxorder = i;
02177 } else {
02178 orders[i] = 0;
02179 if (overlaparray[i]>maxfalseov) maxfalseov = overlaparray[i];
02180 }
02181 }
02182
02183 #if defined(_MSC_VER)
02184 float tol1 = 0.25f + 0.15f*float(myerfc(0.05f*float(sel->selected-50)));
02185 #elif defined(__irix)
02186 float tol1 = 0.25f + 0.15f*erfc(0.05f*float(sel->selected-50));
02187 #else
02188 float tol1 = 0.25f + 0.15f*erfcf(0.05f*float(sel->selected-50));
02189 #endif
02190
02191
02192
02193
02194 if (maxover<tol1+powf(0.09f*maxorder, 6)) {
02195 *overlap = maxover;
02196 return 1;
02197 }
02198
02199
02200
02201
02202 if (orders[2]) {
02203 if (orders[3] && orders[6]) {
02204 float avgov = (overlaparray[2]+overlaparray[3]+overlaparray[6])/3.0f;
02205 *overlap = 0.5f*(maxover + (avgov-maxfalseov)/avgov);
02206 return 6;
02207 } else if (orders[4]) {
02208 if (MAXORDERCN>=8 && orders[8]) {
02209 float avgov = (overlaparray[2]+overlaparray[4]+overlaparray[8])/3.0f;
02210 *overlap = 0.5f*(maxover + (avgov-maxfalseov)/avgov);
02211 return 8;
02212 }
02213 *overlap = 0.5f*(maxover + (overlaparray[4]-maxfalseov)/overlaparray[4]);
02214 return 4;
02215 }
02216
02217 *overlap = 0.5f*(maxover + (overlaparray[2]-maxfalseov)/overlaparray[2]);
02218 return 2;
02219 } else if (orders[3]) {
02220 *overlap = 0.5f*(maxover + (overlaparray[3]-maxfalseov)/overlaparray[3]);
02221 return 3;
02222 } else if (orders[5]) {
02223 *overlap = 0.5f*(maxover + (overlaparray[5]-maxfalseov)/overlaparray[5]);
02224 return 5;
02225 } else if (orders[7]) {
02226 *overlap = 0.5f*(maxover + (overlaparray[7]-maxfalseov)/overlaparray[7]);
02227 return 7;
02228 }
02229
02230 *overlap = maxover;
02231 return 1;
02232 }
02233
02234
02235
02236
02237
02238 void Symmetry::idealize_coordinates() {
02239 int i;
02240
02241
02242 memcpy(idealcoor, coor, 3L*sel->selected*sizeof(float));
02243
02244
02245 int *keep = new int[planes.num()];
02246 memset(keep, 0, planes.num()*sizeof(int));
02247
02248 for (i=0; i<numplanes(); i++) {
02249 Matrix4 mirror;
02250 mirror.translate(rcom[0], rcom[1], rcom[2]);
02251 mirror.scale(-1.0);
02252 mirror.rotate_axis(planes[i].v, float(VMD_PI));
02253 mirror.translate(-rcom[0], -rcom[1], -rcom[2]);
02254
02255
02256 if (average_coordinates(&mirror)) keep[i] = 1;
02257 }
02258
02259 prune_planes(keep);
02260
02261
02262
02263 int *weight = new int[sel->selected];
02264 float *avgcoor = new float[3L*sel->selected];
02265
02266 delete [] keep;
02267 keep = new int[axes.num()];
02268 memset(keep, 0, axes.num()*sizeof(int));
02269
02270 for (i=0; i<numaxes(); i++) {
02271 memset(weight, 0, sel->selected*sizeof(int));
02272 memcpy(avgcoor, idealcoor, 3L*sel->selected*sizeof(int));
02273 int order = axes[i].order;
02274
02275
02276
02277
02278 if (order==INFINITE_ORDER) order = 4;
02279
02280
02281 int success = 1;
02282 int n, j;
02283 for (n=1; n<order; n++) {
02284 Matrix4 rot;
02285 rot.translate(rcom[0], rcom[1], rcom[2]);
02286 rot.rotate_axis(axis(i), n*float(VMD_TWOPI)/order);
02287 rot.translate(-rcom[0], -rcom[1], -rcom[2]);
02288
02289
02290
02291
02292
02293 int nmatches;
02294 int *matchlist = NULL;
02295 identify_transformed_atoms(&rot, nmatches, matchlist);
02296
02297 for(j=0; j<sel->selected; j++) {
02298 int m = matchlist[j];
02299 if (m<0) continue;
02300
02301 if (checkbonds) {
02302
02303
02304 if (!check_bondsum(j, m, &rot)) {
02305 success = 0;
02306 break;
02307 }
02308 }
02309
02310 float tmpcoor[3];
02311 rot.multpoint3d(idealcoor+3L*m, tmpcoor);
02312 vec_incr(avgcoor+3L*j, tmpcoor);
02313 weight[j]++;
02314 }
02315 if (matchlist) delete [] matchlist;
02316
02317 if (!success) break;
02318
02319 }
02320
02321 if (success) {
02322
02323
02324 for(j=0; j<sel->selected; j++) {
02325 vec_scale(idealcoor+3L*j, 1.0f/(1+weight[j]), avgcoor+3L*j);
02326 }
02327
02328 keep[i] = 1;
02329 }
02330 }
02331
02332 prune_axes(keep);
02333
02334 delete [] weight;
02335 delete [] avgcoor;
02336
02337
02338
02339 delete [] keep;
02340 keep = new int[rotreflections.num()];
02341 memset(keep, 0, rotreflections.num()*sizeof(int));
02342
02343 for (i=0; i<numrotreflect(); i++) {
02344 Matrix4 rot;
02345 rot.translate(rcom[0], rcom[1], rcom[2]);
02346 rot.rotate_axis(rotreflect(i), float(VMD_TWOPI)/rotreflections[i].order);
02347 rot.scale(-1.0f);
02348 rot.rotate_axis(rotreflect(i), float(VMD_PI));
02349 rot.translate(-rcom[0], -rcom[1], -rcom[2]);
02350
02351
02352 if (average_coordinates(&rot)) keep[i] = 1;
02353 }
02354
02355 prune_rotreflections(keep);
02356
02357 delete [] keep;
02358
02359
02360 if (inversion) {
02361
02362 Matrix4 inv;
02363 inv.translate(rcom[0], rcom[1], rcom[2]);
02364 inv.scale(-1.0f);
02365 inv.translate(-rcom[0], -rcom[1], -rcom[2]);
02366
02367 if (!average_coordinates(&inv)) inversion = 0;
02368 }
02369 }
02370
02371
02372
02373 int Symmetry::average_coordinates(Matrix4 *trans) {
02374 int j;
02375 int nmatches;
02376 int *matchlist = NULL;
02377 identify_transformed_atoms(trans, nmatches, matchlist);
02378
02379 if (checkbonds) {
02380 int success = 1;
02381 for(j=0; j<sel->selected; j++) {
02382 int m = matchlist[j];
02383
02384 if (m<0) continue;
02385
02386 if (!check_bondsum(j, m, trans)) {
02387 success = 0;
02388 break;
02389 }
02390 }
02391 if (!success) {
02392 if (verbose>3) {
02393 msgInfo << "Transformation messes up bonds!" << sendmsg;
02394 }
02395 if (matchlist) delete [] matchlist;
02396 return 0;
02397 }
02398 }
02399
02400 for(j=0; j<sel->selected; j++) {
02401 int m = matchlist[j];
02402
02403 if (m<0) continue;
02404
02405
02406 float avgcoor[3];
02407 trans->multpoint3d(idealcoor+3L*m, avgcoor);
02408 vec_incr(avgcoor, idealcoor+3L*j);
02409 vec_scale(idealcoor+3L*j, 0.5, avgcoor);
02410 }
02411 if (matchlist) delete [] matchlist;
02412
02413 return 1;
02414 }
02415
02416
02417
02418 int Symmetry::check_bondsum(int j, int m, Matrix4 *trans) {
02419 float tmp[3];
02420 vec_add(tmp, bondsum+3L*m, rcom);
02421 trans->multpoint3d(tmp, tmp);
02422 vec_sub(tmp, tmp, rcom);
02423
02424 if (distance(bondsum+3L*j, tmp)>BONDSUMTOL) {
02425 if (verbose>4) {
02426 char buf[256] = { 0 };
02427 sprintf(buf, "Bond mismatch %i-->%i, bondsum distance = %.2f",
02428 atomindex[j], atomindex[m],
02429 distance(bondsum+3L*j, tmp));
02430 msgInfo << buf << sendmsg;
02431 }
02432
02433
02434
02435 return 0;
02436 }
02437 return 1;
02438 }
02439
02440
02441
02442
02443
02444
02445
02446 void Symmetry::identify_transformed_atoms(Matrix4 *trans,
02447 int &nmatches,
02448 int *(&matchlist)) {
02449
02450 const float *posA = coor;
02451
02452 float *posB = new float[3L*sel->selected];
02453
02454 if (matchlist) delete [] matchlist;
02455 matchlist = new int[sel->selected];
02456
02457
02458 int i;
02459 for(i=0; i<sel->selected; i++) {
02460 trans->multpoint3d(posA+3L*i, posB+3L*i);
02461 }
02462
02463 nmatches = 0;
02464 for(i=0; i<sel->selected; i++) {
02465 float minr2=999999.f;
02466 int k, kmatch = -1;
02467 for(k=0; k<sel->selected; k++) {
02468
02469 if (atomtype[i]==atomtype[k]) {
02470 #if 0
02471 if (checkbonds) {
02472 float imagebondsum[3];
02473 vec_add(imagebondsum, bondsum+3L*k, rcom);
02474 trans->multpoint3d(imagebondsum, imagebondsum);
02475 vec_sub(imagebondsum, imagebondsum, rcom);
02476
02477 if (distance(bondsum+3L*i, imagebondsum)>BONDSUMTOL) {
02478 continue;
02479 }
02480 }
02481 #endif
02482 float r2 = distance2(posA+3L*i, posB+3L*k);
02483
02484 if (r2<minr2) { minr2 = r2; kmatch = k; }
02485 }
02486 }
02487
02488 if (kmatch>=0) {
02489 matchlist[i] = kmatch;
02490 nmatches++;
02491
02492
02493 }
02494 else {
02495
02496 matchlist[i] = i;
02497
02498 if (verbose>3) {
02499 char buf[256] = { 0 };
02500 sprintf(buf, "No match for atom %i (atomtype %i)\n", atomindex[i], atomtype[i]);
02501 msgInfo << buf << sendmsg;
02502 }
02503 }
02504 }
02505
02506 delete [] posB;
02507 }
02508
02509
02510
02511 float Symmetry::ideal_coordinate_rmsd () {
02512
02513 const float *pos = sel->coordinates(mlist);
02514
02515 float rmsdsum = 0.0;
02516 int i, j=0;
02517 for (i=0; i<sel->num_atoms; i++) {
02518 if (sel->on[i]) {
02519 rmsdsum += distance2(pos+3L*i, idealcoor+3L*j);
02520 j++;
02521 }
02522 }
02523
02524 return sqrtf(rmsdsum / sel->selected);
02525 }
02526
02527
02528 int Symmetry::ideal_coordinate_sanity() {
02529 int i, j;
02530 float mindist2 = float(MINATOMDIST*MINATOMDIST);
02531 for (i=0; i<sel->selected; i++) {
02532 for (j=i+1; j<sel->selected; j++) {
02533 if (distance2(idealcoor+3l*i, idealcoor+3L*j)<mindist2)
02534 return 0;
02535 }
02536 }
02537
02538 return 1;
02539 }
02540
02541
02542
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562 void Symmetry::idealize_elements() {
02563 int i;
02564
02565
02566 int *idealp = NULL;
02567 if (planes.num()) {
02568 idealp = new int[planes.num()];
02569 memset(idealp, 0, planes.num()*sizeof(int));
02570 }
02571
02572
02573
02574
02575
02576
02577 if (axes.num()) {
02578
02579
02580 if (horizontalplane>=0) {
02581
02582
02583 vec_copy(planes[horizontalplane].v, axes[0].v);
02584 idealp[horizontalplane] = 1;
02585 }
02586
02587 if (numverticalplanes) {
02588
02589
02590 int vref = -1;
02591 for (i=0; i<planes.num(); i++) {
02592 if (planes[i].type&VERTICALPLANE) {
02593
02594 vref = i;
02595 float normal[3];
02596 cross_prod(normal, axes[0].v, plane(vref));
02597 cross_prod(planes[vref].v, axes[0].v, normal);
02598 vec_normalize(planes[vref].v);
02599 idealp[vref] = 1;
02600 break;
02601 }
02602 }
02603
02604
02605
02606
02607 for (i=vref+1; i<planes.num(); i++) {
02608 if (planes[i].type&VERTICALPLANE) {
02609 align_plane_with_axis(plane(i), axes[0].v, plane(i));
02610
02611
02612 float idealplane[3];
02613 if (!idealize_angle(planes[vref].v, axes[0].v, plane(i), idealplane, axes[0].order)) {
02614 if (verbose>3)
02615 msgInfo << "Couldn't idealize vertical plane " << i << sendmsg;
02616 continue;
02617 }
02618
02619 int first = plane_exists(idealplane);
02620 if (first>=0) {
02621
02622
02623
02624 if (!idealp[first]) {
02625 vec_copy(planes[first].v, idealplane);
02626 idealp[first] = 1;
02627 }
02628 } else {
02629 vec_copy(planes[i].v, idealplane);
02630 idealp[i] = 1;
02631 }
02632 }
02633 }
02634 }
02635
02636
02637 } else {
02638
02639
02640 if (planes.num()<=1) {
02641 if (idealp) delete [] idealp;
02642 return;
02643 }
02644
02645
02646
02647
02648
02649
02650
02651 int bestplane = 0;
02652 float bestscore = 0.f;
02653 for (i=0; i<planes.num(); i++) {
02654 float score = planes[i].overlap*planes[i].weight;
02655 if (score>bestscore) {
02656 bestplane = i;
02657 bestscore = score;
02658 }
02659 idealp[i] = 0;
02660 }
02661 idealp[bestplane] = 1;
02662
02663 if (verbose>4) {
02664 msgErr << "Found planes without intersection axis!" << sendmsg;
02665 char buf[256] = { 0 };
02666 for (i=0; i<planes.num(); i++) {
02667 sprintf(buf, "plane[%i] {% .3f % .3f % .3f} overlap = %f, weight = %d", i,
02668 planes[i].v[0], planes[i].v[1], planes[i].v[2],
02669 planes[i].overlap, planes[i].weight);
02670 msgErr << buf << sendmsg;
02671 }
02672 }
02673 }
02674
02675 int numidealaxes = 0;
02676 int *ideala = NULL;
02677 if (axes.num()) {
02678 ideala = new int[axes.num()];
02679 memset(ideala, 0, axes.num()*sizeof(int));
02680 ideala[0] = 1;
02681 }
02682
02683 int geometry = 0;
02684 if (maxaxisorder==2) geometry = 4;
02685 else if (maxaxisorder==3) geometry = TETRAHEDRON;
02686 else if (maxaxisorder==4) geometry = OCTAHEDRON;
02687 else if (maxaxisorder==5) geometry = DODECAHEDRON;
02688
02689
02690 for (i=1; i<numaxes(); i++) {
02691 if (axes[i].order<maxaxisorder) continue;
02692 int j, vref=-1;
02693
02694
02695
02696
02697
02698 for (j=0; j<numplanes(); j++) {
02699 if (orthogonal(planes[j].v, axes[i].v, orthogtol) &&
02700 orthogonal(planes[j].v, axes[0].v, orthogtol)) {
02701 if (!idealp[j]) {
02702
02703
02704 continue;
02705 }
02706
02707 float idealaxis[3];
02708 if (!idealize_angle(axes[0].v, planes[j].v, axes[i].v,
02709 idealaxis, geometry)) {
02710 if (verbose>4) {
02711 msgInfo << "Couldn't idealize axis " << i
02712 << " wrt reference axis in plane " << j
02713 << "." << sendmsg;
02714 }
02715 continue;
02716 }
02717 vec_copy(axes[i].v, idealaxis);
02718 vref = j;
02719 ideala[i] = 1;
02720 break;
02721 }
02722 }
02723
02724 if (vref<0) continue;
02725
02726
02727
02728 for (j=0; j<planes.num(); j++) {
02729 if (idealp[j]) continue;
02730
02731
02732 if (orthogonal(planes[j].v, axes[i].v, orthogtol)) {
02733
02734 align_plane_with_axis(planes[j].v, axes[i].v, planes[j].v);
02735
02736
02737
02738 float idealplane[3];
02739 if (!idealize_angle(planes[vref].v, axes[i].v, plane(j), idealplane, axes[i].order)) {
02740 if (verbose>4) {
02741 msgInfo << "Vertical plane " <<j<< " couldn't be idealized!" << sendmsg;
02742 }
02743 continue;
02744 }
02745 int first = plane_exists(idealplane);
02746 if (first>=0) {
02747
02748
02749
02750 if (!idealp[first]) {
02751 vec_copy(planes[first].v, idealplane);
02752 idealp[first] = 1;
02753 }
02754 } else {
02755 vec_copy(planes[j].v, idealplane);
02756 idealp[j] = 1;
02757 }
02758 }
02759 }
02760 }
02761
02762
02763 prune_planes(idealp);
02764
02765
02766
02767
02768 for (i=0; i<numplanes(); i++) {
02769 int j;
02770 for (j=i+1; j<numplanes(); j++) {
02771
02772 float intersect[3];
02773 cross_prod(intersect, plane(i), plane(j));
02774 vec_normalize(intersect);
02775
02776 int k;
02777 for (k=0; k<axes.num(); k++) {
02778 if (ideala[k]) continue;
02779
02780 if (collinear(intersect, axes[k].v, collintol)) {
02781
02782 vec_copy(axes[k].v, intersect);
02783 ideala[k] = 1; numidealaxes++;
02784 break;
02785 }
02786 }
02787 }
02788 }
02789
02790 float halfcollintol = cosf(float(DEGTORAD(5.0f)));
02791
02792
02793
02794 for (i=0; i<numplanes(); i++) {
02795 int j;
02796 for (j=i+1; j<numplanes(); j++) {
02797
02798 float bisect1[3];
02799 vec_add(bisect1, planes[i].v, planes[j].v);
02800 vec_normalize(bisect1);
02801
02802
02803 float bisect2[3], tmp[3];
02804 vec_negate(tmp, planes[i].v);
02805 vec_add(bisect2, tmp, planes[j].v);
02806 vec_normalize(bisect2);
02807
02808 int k;
02809 int foundbi1=0, foundbi2=0;
02810 for (k=0; k<axes.num(); k++) {
02811 if (ideala[k]) continue;
02812
02813 if (!foundbi1 && collinear(bisect1, axes[k].v, halfcollintol)) {
02814
02815 vec_copy(axes[k].v, bisect1);
02816 ideala[k] = 1; numidealaxes++;
02817 foundbi1 = 1;
02818 if (foundbi2) break;
02819 }
02820 if (!foundbi2 && collinear(bisect2, axes[k].v, halfcollintol)) {
02821
02822 vec_copy(axes[k].v, bisect2);
02823 ideala[k] = 1; numidealaxes++;
02824 foundbi2 = 1;
02825 if (foundbi1) break;
02826 }
02827 }
02828 }
02829 }
02830
02831
02832
02833
02834
02835 if (numidealaxes<axes.num()) {
02836 int firstorth = -1;
02837 for (i=0; i<numaxes(); i++) {
02838 if (ideala[i]) continue;
02839
02840 if (!orthogonal(axes[i].v, axes[0].v, orthogtol)) continue;
02841
02842
02843 float tmp[3];
02844 cross_prod(tmp, axes[0].v, axes[i].v);
02845 vec_normalize(tmp);
02846 cross_prod(axes[i].v, tmp, axes[0].v);
02847
02848 if (firstorth<0) {
02849 firstorth = i;
02850 ideala[i] = 1;
02851 continue;
02852 }
02853
02854
02855 if (!idealize_angle(axes[firstorth].v, axes[0].v, axes[i].v, tmp, axes[0].order)) {
02856 if (verbose>4) {
02857 msgInfo << "Couldn't idealize axis "<<i<<" to first orthogonal axis "
02858 <<firstorth<< "!" << sendmsg;
02859 }
02860 continue;
02861 }
02862
02863 vec_copy(axes[i].v, tmp);
02864 ideala[i] = 1;
02865 }
02866 }
02867
02868
02869 prune_axes(ideala);
02870
02871 if (ideala) delete [] ideala;
02872 if (idealp) delete [] idealp;
02873
02874
02875
02876 if (numrotreflect()) {
02877 int *idealr = new int[numrotreflect()];
02878 memset(idealr, 0, numrotreflect()*sizeof(int));
02879
02880 for (i=0; i<numrotreflect(); i++) {
02881 int j, success=0;
02882 for (j=0; j<numaxes(); j++) {
02883 float dot = dot_prod(rotreflections[i].v, axes[j].v);
02884 if (fabs(dot) > collintol) {
02885 float tmp[3];
02886 if (dot>0) vec_negate(tmp, axes[j].v);
02887 else vec_copy(tmp, axes[j].v);
02888
02889 vec_copy(rotreflections[i].v, tmp);
02890 rotreflections[i].type = j;
02891 success = 1;
02892 break;
02893 }
02894 }
02895
02896 if (success) {
02897
02898 idealr[i] = 1;
02899 }
02900 }
02901
02902 prune_rotreflections(idealr);
02903 delete [] idealr;
02904 }
02905
02906 }
02907
02908
02909
02910
02911
02912
02913
02914
02915
02916
02917
02918
02919
02920
02921
02922 int Symmetry::idealize_angle(const float *refaxis, const float *hub,
02923 const float *myaxis, float *idealaxis, int reforder) {
02924 float alpha = angle(refaxis, myaxis);
02925
02926 const float tetrahedron = float(RADTODEG(2.0f*atanf(sqrtf(2.0f))));
02927 const float octahedron = 90.0f;
02928 const float dodecahedron = float(RADTODEG(acosf(-1/sqrtf(5.0f))));
02929 const float tol = 5.0f;
02930
02931 int success = 0;
02932 float idealangle=0.0f;
02933
02934 if (reforder==TETRAHEDRON) idealangle = tetrahedron;
02935 else if (reforder==OCTAHEDRON) idealangle = octahedron;
02936 else if (reforder==DODECAHEDRON) idealangle = dodecahedron;
02937
02938 if (reforder<0) {
02939 if (fabs(idealangle-alpha)<tol) {
02940 alpha = idealangle;
02941 success = 1;
02942
02943 } else if (fabs(180.0-idealangle-alpha)<tol) {
02944 alpha = 180-idealangle;
02945 success = 1;
02946 }
02947 }
02948
02949 int i;
02950 for (i=1; i<reforder; i++) {
02951 idealangle = i*180.0f/reforder;
02952 if (fabs(alpha-idealangle)<tol) {
02953 alpha = idealangle;
02954 success = 1;
02955 break;
02956 }
02957 }
02958
02959
02960 if (fabs(alpha)<tol || fabs(alpha-180)<tol) {
02961
02962 alpha = 0;
02963 success = 1;
02964 }
02965
02966 if (!success) {
02967
02968 return 0;
02969 }
02970
02971 float normal[3];
02972 cross_prod(normal, refaxis, myaxis);
02973 if (dot_prod(hub, normal) < 0) alpha = -alpha;
02974
02975
02976
02977 Matrix4 rot;
02978 rot.rotate_axis(hub, float(DEGTORAD(alpha)));
02979 rot.multpoint3d(refaxis, idealaxis);
02980
02981 return 1;
02982 }
02983
02984
02985
02986
02987
02988 void Symmetry::collapse_axis(const float *axis, int order,
02989 int refatom, const int *matchlist,
02990 int *(&connectedunique)) {
02991 int i;
02992 float refcoor[3];
02993 vec_sub(refcoor, coor+3L*refatom, rcom);
02994
02995
02996 float r0[3];
02997 vec_scale(r0, dot_prod(refcoor, axis), axis);
02998 vec_sub(r0, refcoor, r0);
02999
03000 for (i=0; i<sel->selected; i++) {
03001 if (!uniqueatoms[i] || i==matchlist[i]) continue;
03002
03003
03004
03005
03006
03007
03008
03009
03010
03011 int image, found=0;
03012 int k = i;
03013 for (image=0; image<order; image++) {
03014 float tmp[3];
03015 vec_sub(tmp, idealcoor+3L*k, rcom);
03016
03017
03018 float r[3];
03019 vec_scale(r, dot_prod(tmp, axis), axis);
03020 vec_sub(r, tmp, r);
03021
03022
03023
03024
03025 if (angle(r, r0) <= 180.0/order) {
03026 found = 1;
03027 break;
03028 }
03029 k = matchlist[k];
03030 }
03031 if (found && k!=i) {
03032
03033 uniqueatoms[i] = 0;
03034 uniqueatoms[k] = 1;
03035 }
03036 }
03037
03038
03039 wrap_unconnected_unique_atoms(refatom, matchlist, connectedunique);
03040 }
03041
03042
03043
03044
03045
03046
03047 void Symmetry::wrap_unconnected_unique_atoms(int root,
03048 const int *matchlist,
03049 int *(&connectedunique)) {
03050 int i, k=0;
03051 int numswapped = 0;
03052
03053
03054
03055
03056 if (!connectedunique) {
03057 connectedunique = new int[sel->selected];
03058 }
03059 find_connected_unique_atoms(connectedunique, root);
03060
03061
03062
03063
03064 do {
03065 numswapped = 0;
03066 for (i=0; i<sel->selected; i++) {
03067
03068 if (!uniqueatoms[i] || connectedunique[i]) continue;
03069
03070 int swap = 0;
03071 int image = 0;
03072 int j = matchlist[i];
03073
03074 while (j!=i && image<sel->selected) {
03075
03076 if (connectedunique[j]) {
03077
03078 swap = 1;
03079 } else {
03080
03081
03082 int k;
03083 for (k=0; k<bondsperatom[j].numbonds; k++) {
03084 int bondto = bondsperatom[j].bondto[k];
03085 if (connectedunique[bondto]) swap = 1;
03086 }
03087 }
03088 if (swap) break;
03089
03090 j = matchlist[j];
03091
03092 image++;
03093 }
03094
03095 if (swap) {
03096
03097 uniqueatoms[i] = 0;
03098 uniqueatoms[j] = 1;
03099 connectedunique[j] = 1;
03100 numswapped++;
03101 }
03102 }
03103
03104 if (k>=sel->selected) {
03105 msgErr << "Stop looping unconnected unique atoms" << sendmsg;
03106 break;
03107 }
03108 k++;
03109
03110 } while (numswapped);
03111
03112 }
03113
03114
03115
03116
03117
03118 void Symmetry::find_connected_unique_atoms(int *(&connectedunique),
03119 int root) {
03120 ResizeArray<int> leaves;
03121 ResizeArray<int> newleaves;
03122 leaves.append(root);
03123
03124 memset(connectedunique, 0, sel->selected*sizeof(int));
03125 connectedunique[root] = 1;
03126
03127 int numbonded = 1;
03128 int i;
03129
03130 do {
03131 newleaves.clear();
03132
03133
03134 for (i=0; i<leaves.num(); i++) {
03135 int j = leaves[i];
03136
03137 int k;
03138 for (k=0; k<bondsperatom[j].numbonds; k++) {
03139 int bondto = bondsperatom[j].bondto[k];
03140
03141 if (uniqueatoms[bondto]&& !connectedunique[bondto]) {
03142 connectedunique[bondto] = 1;
03143 newleaves.append(bondto);
03144 numbonded++;
03145 }
03146 }
03147
03148 }
03149
03150 leaves.clear();
03151 for (i=0; i<newleaves.num(); i++) {
03152 leaves.append(newleaves[i]);
03153 }
03154
03155 } while (newleaves.num());
03156
03157
03158
03159
03160
03161 }
03162
03163
03164
03165
03166
03167
03168
03169
03170 void Symmetry::unique_coordinates() {
03171 int i;
03172 for(i=0; i<sel->selected; i++) {
03173 uniqueatoms[i] = 1;
03174 }
03175
03176
03177
03178
03179
03180 float refcoor[3];
03181 int refatom = 0;
03182 vec_sub(refcoor, coor, rcom);
03183 if (norm(refcoor)<0.1) {
03184 refatom = 1;
03185 vec_sub(refcoor, coor+3L*refatom, rcom);
03186 }
03187
03188 int *connectedunique = NULL;
03189
03190
03191 if (inversion) {
03192 Matrix4 inv;
03193 inv.translate(rcom[0], rcom[1], rcom[2]);
03194 inv.scale(-1.0);
03195 inv.translate(-rcom[0], -rcom[1], -rcom[2]);
03196
03197
03198 int *matchlist = unique_coordinates(&inv);
03199
03200 int j;
03201 for(j=0; j<sel->selected; j++) {
03202 if (!uniqueatoms[j] || j==matchlist[j]) continue;
03203 uniqueatoms[j] = 0;
03204 uniqueatoms[matchlist[j]] = 1;
03205 }
03206
03207 wrap_unconnected_unique_atoms(refatom, matchlist, connectedunique);
03208
03209 if (matchlist) delete [] matchlist;
03210 }
03211
03212
03213 for (i=0; i<numaxes(); i++) {
03214 Matrix4 rot;
03215 rot.translate(rcom[0], rcom[1], rcom[2]);
03216 rot.rotate_axis(axes[i].v, float(VMD_TWOPI)/axes[i].order);
03217 rot.translate(-rcom[0], -rcom[1], -rcom[2]);
03218
03219
03220 int *matchlist = unique_coordinates(&rot);
03221
03222
03223
03224
03225 collapse_axis(axes[i].v, axes[i].order, refatom, matchlist,
03226 connectedunique);
03227
03228 if (matchlist) delete [] matchlist;
03229 }
03230
03231
03232 for (i=0; i<numplanes(); i++) {
03233 Matrix4 mirror;
03234 mirror.translate(rcom[0], rcom[1], rcom[2]);
03235 mirror.scale(-1.0);
03236 mirror.rotate_axis(planes[i].v, float(VMD_PI));
03237 mirror.translate(-rcom[0], -rcom[1], -rcom[2]);
03238
03239
03240 int *matchlist = unique_coordinates(&mirror);
03241
03242 int j;
03243 for(j=0; j<sel->selected; j++) {
03244 if (!uniqueatoms[j] || j==matchlist[j]) continue;
03245
03246 float tmp[3];
03247 vec_sub(tmp, coor+3L*j, rcom);
03248 if (behind_plane(planes[i].v, tmp)!=behind_plane(planes[i].v, refcoor)) {
03249 uniqueatoms[j] = 0;
03250 uniqueatoms[matchlist[j]] = 1;
03251 }
03252 }
03253 if (matchlist) delete [] matchlist;
03254 }
03255
03256
03257 for (i=0; i<numrotreflect(); i++) {
03258 Matrix4 rotref;
03259 rotref.translate(rcom[0], rcom[1], rcom[2]);
03260 rotref.rotate_axis(rotreflections[i].v, float(VMD_TWOPI)/rotreflections[i].order);
03261 rotref.scale(-1.0);
03262 rotref.rotate_axis(rotreflections[i].v, float(VMD_PI));
03263 rotref.translate(-rcom[0], -rcom[1], -rcom[2]);
03264
03265
03266 int *matchlist = unique_coordinates(&rotref);
03267
03268 collapse_axis(rotreflections[i].v, rotreflections[i].order, refatom, matchlist, connectedunique);
03269
03270 if (matchlist) delete [] matchlist;
03271 }
03272
03273 if (connectedunique) delete [] connectedunique;
03274 }
03275
03276
03277
03278
03279
03280
03281 int* Symmetry::unique_coordinates(Matrix4 *trans) {
03282 int nmatches;
03283 int *matchlist = NULL;
03284
03285
03286
03287
03288
03289 identify_transformed_atoms(trans, nmatches, matchlist);
03290
03291 int j;
03292 for(j=0; j<sel->selected; j++) {
03293 if (!uniqueatoms[j]) continue;
03294 int m = matchlist[j];
03295
03296 if (m<0) continue;
03297
03298 if (m>j) uniqueatoms[m] = 0;
03299 }
03300
03301 return matchlist;
03302 }
03303
03304
03306 void Symmetry::determine_pointgroup() {
03307 if (linear) {
03308 if (inversion) pointgroup = POINTGROUP_DINFH;
03309 else pointgroup = POINTGROUP_CINFV;
03310
03311 pointgrouporder = -1;
03312 }
03313
03314 else if (sphericaltop && maxaxisorder>=3 && axes[0].order>=2) {
03315 if (maxaxisorder==3) {
03316 if (numplanes()) {
03317 if (inversion) pointgroup = POINTGROUP_TH;
03318 else pointgroup = POINTGROUP_TD;
03319 }
03320 else pointgroup = POINTGROUP_T;
03321 }
03322 else if (maxaxisorder==4) {
03323 if (numplanes() || inversion) pointgroup = POINTGROUP_OH;
03324 else pointgroup = POINTGROUP_O;
03325 }
03326 else if (maxaxisorder==5) {
03327 if (numplanes() || inversion) pointgroup = POINTGROUP_IH;
03328 else pointgroup = POINTGROUP_I;
03329
03330 }
03331 else pointgroup = POINTGROUP_UNKNOWN;
03332
03333 }
03334
03335 else if (numaxes()) {
03336
03337
03338 int i;
03339 int perpC2 = 0;
03340 for (i=0; i<numaxes(); i++) {
03341 if (axes[i].order==2 && (axes[i].type & PERPENDICULAR_AXIS)) {
03342 perpC2++;
03343 }
03344 }
03345
03346 if (perpC2==maxaxisorder) {
03347 if (horizontalplane>=0) pointgroup = POINTGROUP_DNH;
03348 else {
03349
03350
03351 if (numdihedralplanes==maxaxisorder) {
03352 pointgroup = POINTGROUP_DND;
03353 }
03354 else {
03355 pointgroup = POINTGROUP_DN;
03356 }
03357 }
03358
03359 pointgrouporder = maxaxisorder;
03360 }
03361
03362 else {
03363 if (horizontalplane>=0) pointgroup = POINTGROUP_CNH;
03364 else {
03365
03366 if (numverticalplanes==maxaxisorder) {
03367 pointgroup = POINTGROUP_CNV;
03368 }
03369 else {
03370 if (numS2n()) {
03371 pointgroup = POINTGROUP_S2N;
03372 }
03373 else {
03374 pointgroup = POINTGROUP_CN;
03375 }
03376 }
03377 }
03378 pointgrouporder = maxaxisorder;
03379
03380 }
03381 }
03382
03383 else {
03384 if (numplanes()==1) pointgroup = POINTGROUP_CS;
03385 else {
03386 if (inversion) pointgroup = POINTGROUP_CI;
03387 else pointgroup = POINTGROUP_C1;
03388 }
03389 }
03390 }
03391
03392
03393
03394
03395
03396
03397 int Symmetry::pointgroup_rank(int pg, int order) {
03398 if (pg==POINTGROUP_C1) return 1;
03399 if (pg==POINTGROUP_CS || pg==POINTGROUP_CI) return 2;
03400 if (pg==POINTGROUP_CN) {
03401 return 1+numprimefactors(order);
03402 }
03403 if (pg==POINTGROUP_S2N) {
03404 return 1+numprimefactors(order*2);
03405 }
03406 if (pg==POINTGROUP_DN || pg==POINTGROUP_CNV ||
03407 pg==POINTGROUP_CNH) {
03408 return 2+numprimefactors(order);
03409 }
03410 if (pg==POINTGROUP_DND || pg==POINTGROUP_DNH) {
03411 return 3+numprimefactors(order);
03412 }
03413 if (pg==POINTGROUP_CINFV) return 3;
03414 if (pg==POINTGROUP_DINFH || pg==POINTGROUP_T) return 4;
03415 if (pg==POINTGROUP_TD || pg==POINTGROUP_TH ||
03416 pg==POINTGROUP_O) {
03417 return 5;
03418 }
03419 if (pg==POINTGROUP_OH || pg==POINTGROUP_I) return 6;
03420 if (pg==POINTGROUP_IH) return 7;
03421 return 0;
03422 }
03423
03424
03425
03426
03427
03428
03429
03430
03431
03432
03433
03434
03435
03436
03437
03438
03439
03440
03441
03442
03443
03444
03445
03446
03447
03448
03449
03450
03451
03452
03453
03454
03455
03456
03457
03458
03459
03460
03461
03462
03463
03464 void Symmetry::orient_molecule() {
03465 if (pointgroup==POINTGROUP_C1) return;
03466
03467
03468
03469
03470 if (linear) {
03471
03472 orient.transvec(0, 0, 1);
03473
03474 orient.transvecinv(axes[0].v[0], axes[0].v[1], axes[0].v[2]);
03475
03476 orient.translate(-rcom[0], -rcom[1], -rcom[2]);
03477 return;
03478 }
03479
03480 int i;
03481 Matrix4 rot;
03482
03483
03484
03485
03486
03487 if (!sphericaltop && numaxes()) {
03488
03489
03490
03491 rot.transvec(0, 0, 1);
03492
03493 rot.transvecinv(axes[0].v[0], axes[0].v[1], axes[0].v[2]);
03494
03495
03496 int j;
03497 for (j=0; j<=2; j++) {
03498 if (orthogonal(axes[0].v, inertiaaxes[j], orthogtol)) break;
03499 }
03500 float *ortho_inertiaaxis = inertiaaxes[j];
03501
03502
03503
03504 Matrix4 m;
03505 float tmp[3];
03506 rot.multpoint3d(ortho_inertiaaxis, tmp);
03507 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03508
03509
03510 m.multmatrix(rot);
03511 rot.loadmatrix(m);
03512 }
03513
03514
03515
03516 int orthC2 = -1;
03517 for (i=1; i<numaxes(); i++) {
03518 if (axes[i].order==2 &&
03519 (orthogonal(axes[i].v, axes[0].v, orthogtol) ||
03520 (pointgroup>=POINTGROUP_T && pointgroup<=POINTGROUP_IH))) {
03521 orthC2 = i;
03522 break;
03523 }
03524 }
03525 if (orthC2>=0) {
03526
03527
03528
03529 float tmp[3];
03530 rot.multpoint3d(axes[orthC2].v, tmp);
03531 Matrix4 m;
03532 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03533
03534
03535 m.multmatrix(rot);
03536 rot.loadmatrix(m);
03537 }
03538
03539
03540
03541 if (numverticalplanes && orthC2<0) {
03542 for (i=0; i<numplanes(); i++) {
03543 if (planes[i].type==VERTICALPLANE) break;
03544 }
03545
03546
03547 Matrix4 m;
03548
03549
03550 float Y[]={0, 1, 0};
03551 m.transvec(Y[0], Y[1], Y[2]);
03552
03553
03554 float tmp[3];
03555 rot.multpoint3d(planes[i].v, tmp);
03556 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03557
03558
03559 m.multmatrix(rot);
03560 rot.loadmatrix(m);
03561 }
03562
03563
03564
03565
03566 if ((horizontalplane>=0 && !numverticalplanes && orthC2<0) ||
03567 pointgroup==POINTGROUP_CS) {
03568
03569 if (pointgroup==POINTGROUP_CS) i = 0;
03570 else i = horizontalplane;
03571
03572 Matrix4 m;
03573 float Z[]={0, 0, 1};
03574
03575
03576 m.transvec(Z[0], Z[1], Z[2]);
03577
03578
03579 float tmp[3];
03580 rot.multpoint3d(planes[i].v, tmp);
03581 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03582
03583
03584 m.multmatrix(rot);
03585 rot.loadmatrix(m);
03586 }
03587
03588 if (pointgroup>=POINTGROUP_T && pointgroup<=POINTGROUP_IH) {
03589
03590
03591
03592 int found = 0;
03593 float uniqueaxis[3];
03594 for (i=1; i<numaxes(); i++) {
03595 if (axes[i].order<maxaxisorder) break;
03596 int j;
03597 for(j=0; j<sel->selected; j++) {
03598 if (!uniqueatoms[j]) continue;
03599
03600 float tmpcoor[3];
03601 vec_sub(tmpcoor, coor+3L*j, rcom);
03602 if (norm(tmpcoor)<0.1) continue;
03603
03604 if (collinear(axes[i].v, tmpcoor, collintol)) {
03605 if (dot_prod(axes[i].v, tmpcoor)>0)
03606 vec_copy(uniqueaxis, axes[i].v);
03607 else
03608 vec_negate(uniqueaxis, axes[i].v);
03609 found = 1;
03610 }
03611 }
03612 }
03613
03614 if (!found)
03615 msgErr << "orient_molecule(): Couldn't find axis with unique atom!" << sendmsg;
03616
03617 Matrix4 m;
03618 float XYZ[]={1, 1, 1};
03619
03620
03621 m.transvec(XYZ[0], XYZ[1], XYZ[2]);
03622
03623
03624 float tmp[3];
03625 rot.multpoint3d(uniqueaxis, tmp);
03626 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03627
03628
03629 m.multmatrix(rot);
03630 rot.loadmatrix(m);
03631 }
03632
03633 orient.multmatrix(rot);
03634 orient.translate(-rcom[0], -rcom[1], -rcom[2]);
03635 }
03636
03637
03640 void Symmetry::print_statistics() {
03641 int i;
03642
03643 #if 0
03644 for (i=0; i<numplanes(); i++) {
03645 printf("plane[%i]: weight=%3i, overlap=%.2f\n", i, planes[i].weight, planes[i].overlap);
03646 }
03647
03648 for (i=0; i<numaxes(); i++) {
03649 printf("axis[%i]: order=%i, weight=%3i, overlap=%.2f {%.2f %.2f %.2f}\n", i, axes[i].order, axes[i].weight, axes[i].overlap, axes[i].v[0], axes[i].v[1], axes[i].v[2]);
03650 }
03651
03652 for (i=0; i<numrotreflect(); i++) {
03653 printf("rotrefl[%i]: order=%i, weight=%3i, overlap=%.2f {%.2f %.2f %.2f}\n", i, rotreflections[i].order, rotreflections[i].weight, rotreflections[i].overlap, rotreflections[i].v[0], rotreflections[i].v[1], rotreflections[i].v[2]);
03654 }
03655 #endif
03656
03657 char buf [256] = { 0 };
03658 msgInfo << sendmsg;
03659 msgInfo << "Summary of symmetry elements:" << sendmsg;
03660 msgInfo << "+===============================+" << sendmsg;
03661 msgInfo << "| inversion center: " << (inversion ? "yes" : "no ")
03662 << " |" << sendmsg;
03663
03664 if (planes.num()) {
03665 int havehorizplane = 0;
03666 msgInfo << "| |" << sendmsg;
03667 if (horizontalplane>=0) {
03668 msgInfo << "| horizonal planes: 1 |" << sendmsg;
03669 havehorizplane = 1;
03670 }
03671 if (numverticalplanes) {
03672 sprintf(buf, "| vertical planes: %2d |", numverticalplanes);
03673 msgInfo << buf << sendmsg;
03674 }
03675 if (numdihedralplanes) {
03676 sprintf(buf, "| (%-2d of them dihedral) |", numdihedralplanes);
03677 msgInfo << buf << sendmsg;
03678 }
03679 int numregplanes = numplanes()-numverticalplanes-havehorizplane;
03680 if (numregplanes) {
03681 sprintf(buf, "| regular planes: %2d |", numregplanes);
03682 msgInfo << buf << sendmsg;
03683 }
03684 if (numplanes()>1) {
03685 msgInfo << "| ----------------------------- |" << sendmsg;
03686 sprintf(buf, "| total planes: %2d |", numplanes());
03687 msgInfo << buf << sendmsg;
03688 }
03689 }
03690
03691 if (axes.num()) {
03692 msgInfo << "| |" << sendmsg;
03693 if (elementsummary.Cinf) {
03694 sprintf(buf, "| Cinf rotation axes: %2d |", (int)elementsummary.Cinf);
03695 msgInfo << buf << sendmsg;
03696 }
03697 for (i=maxaxisorder; i>=1; i--) {
03698 if (elementsummary.C[i]) {
03699 sprintf(buf, "| C%-4i rotation axes: %2d |", i, elementsummary.C[i]);
03700 msgInfo << buf << sendmsg;
03701 }
03702 }
03703 if (numaxes()>1) {
03704 msgInfo << "| ----------------------------- |" << sendmsg;
03705 sprintf(buf, "| total rotation axes: %2d |", numaxes());
03706 msgInfo << buf << sendmsg;
03707 }
03708 }
03709
03710 if (rotreflections.num()) {
03711 msgInfo << "| |" << sendmsg;
03712 for (i=maxrotreflorder; i>=3; i--) {
03713 if (elementsummary.S[i-3]) {
03714 sprintf(buf, "| S%-4i rotary reflections: %2d |", i, elementsummary.S[i-3]);
03715 msgInfo << buf << sendmsg;
03716 }
03717 }
03718 if (rotreflections.num()>1) {
03719 msgInfo << "| ----------------------------- |" << sendmsg;
03720 sprintf(buf, "| total rotary reflections: %2ld |", long(rotreflections.num()));
03721 msgInfo << buf << sendmsg;
03722 }
03723 }
03724 msgInfo << "+===============================+" << sendmsg;
03725
03726 msgInfo << sendmsg;
03727 msgInfo << "Element summary: " << elementsummarystring << sendmsg;
03728 }
03729
03730
03731
03732 void Symmetry::build_element_summary() {
03733 memset(&elementsummary, 0, sizeof(ElementSummary));
03734
03735 elementsummary.inv = inversion;
03736 elementsummary.sigma = planes.num();
03737
03738 int i;
03739 for (i=0; i<numaxes(); i++) {
03740 if (axes[i].order==INFINITE_ORDER) {
03741 elementsummary.Cinf++;
03742 } else if (axes[i].order<=MAXORDERCN) {
03743 int j;
03744 for (j=2; j<=axes[i].order; j++) {
03745 if (axes[i].order % j == 0) {
03746 (elementsummary.C[j])++;
03747 }
03748 }
03749 }
03750 }
03751
03752 for (i=0; i<numrotreflect(); i++) {
03753 if (rotreflections[i].order<=2*MAXORDERCN) {
03754 (elementsummary.S[rotreflections[i].order-3])++;
03755 }
03756 }
03757 }
03758
03759
03760
03761
03762 void Symmetry::build_element_summary_string(ElementSummary summary, char *(&sestring)) {
03763 int i ;
03764 if (sestring) delete [] sestring;
03765 sestring = new char[2 + 10L*(MAXORDERCN+2*MAXORDERCN+summary.Cinf
03766 +(summary.sigma?1:0)+summary.inv)];
03767 char buf[100] = { 0 };
03768 sestring[0] = '\0';
03769
03770 if (inversion) strcat(sestring, "(inv) ");
03771
03772 if (summary.sigma==1) strcat(sestring, "(sigma) ");
03773 if (summary.sigma>1) {
03774 sprintf(buf, "%d*(sigma) ", summary.sigma);
03775 strcat(sestring, buf);
03776 }
03777
03778 if (summary.Cinf==1)
03779 strcat(sestring, "(Cinf) ");
03780 else if (summary.Cinf>1) {
03781 sprintf(buf, "%d*(Cinf) ", summary.Cinf);
03782 strcat(sestring, buf);
03783 }
03784
03785 for (i=MAXORDERCN; i>=2; i--) {
03786 if (summary.C[i]==1) {
03787 sprintf(buf, "(C%d) ", i);
03788 strcat(sestring, buf);
03789 }
03790 else if (summary.C[i]>1) {
03791 sprintf(buf, "%d*(C%d) ", summary.C[i], i);
03792 strcat(sestring, buf);
03793 }
03794 }
03795
03796 for (i=2*MAXORDERCN; i>=3; i--) {
03797 if (summary.S[i-3]==1) {
03798 sprintf(buf, "(S%d) ", i);
03799 strcat(sestring, buf);
03800 }
03801 else if (summary.S[i-3]>1) {
03802 sprintf(buf, "%d*(S%d) ", summary.S[i-3], i);
03803 strcat(sestring, buf);
03804 }
03805 }
03806 }
03807
03808
03809
03810
03811
03812
03813 void Symmetry::compare_element_summary(const char *pointgroupname) {
03814 missingelementstring[0] = '\0';
03815 additionalelementstring[0] = '\0';
03816
03817 if (!strcmp(pointgroupname, "Unknown")) return;
03818
03819 unsigned int i;
03820 for (i=0; i<NUMPOINTGROUPS; i++) {
03821 if (!strcmp(pointgroupname, pgdefinitions[i].name)) {
03822
03823 if (elementsummary.inv<pgdefinitions[i].summary.inv)
03824 strcat(missingelementstring, "(inv) ");
03825 else if (elementsummary.inv>pgdefinitions[i].summary.inv)
03826 strcat(additionalelementstring, "(inv) ");
03827
03828 char buf[100] = { 0 };
03829 if (elementsummary.sigma<pgdefinitions[i].summary.sigma) {
03830 sprintf(buf, "%i*(sigma) ",
03831 pgdefinitions[i].summary.sigma-elementsummary.sigma);
03832 strcat(missingelementstring, buf);
03833 }
03834 else if (elementsummary.sigma>pgdefinitions[i].summary.sigma) {
03835 sprintf(buf, "%i*(sigma) ",
03836 elementsummary.sigma-pgdefinitions[i].summary.sigma);
03837 strcat(additionalelementstring, buf);
03838 }
03839
03840 int j;
03841 for (j=MAXORDERCN; j>=2; j--) {
03842 if (elementsummary.C[j]<pgdefinitions[i].summary.C[j]) {
03843 sprintf(buf, "%i*(C%i) ", pgdefinitions[i].summary.C[j]-elementsummary.C[j], j);
03844 strcat(missingelementstring, buf);
03845 }
03846 if (elementsummary.C[j]>pgdefinitions[i].summary.C[j]) {
03847 sprintf(buf, "%i*(C%i) ", elementsummary.C[j]-pgdefinitions[i].summary.C[j], j);
03848 strcat(additionalelementstring, buf);
03849 }
03850 }
03851
03852 for (j=2*MAXORDERCN; j>=3; j--) {
03853 if (elementsummary.S[j-3]<pgdefinitions[i].summary.S[j-3]) {
03854 sprintf(buf, "%i*(S%i) ",
03855 pgdefinitions[i].summary.S[j-3]-elementsummary.S[j-3], j);
03856 strcat(missingelementstring, buf);
03857 }
03858 if (elementsummary.S[j-3]>pgdefinitions[i].summary.S[j-3]) {
03859 sprintf(buf, "%i*(S%i) ",
03860 elementsummary.S[j-3]-pgdefinitions[i].summary.S[j-3], j);
03861 strcat(additionalelementstring, buf);
03862 }
03863 }
03864 break;
03865 }
03866 }
03867 }
03868
03869
03870 void Symmetry::print_element_summary(const char *pointgroupname) {
03871 int i, found = 0;
03872 for (i=0; i<(int)NUMPOINTGROUPS; i++) {
03873 if (!strcmp(pointgroupname, pgdefinitions[i].name)) {
03874 found = 1;
03875 break;
03876 }
03877 }
03878 if (!found) return;
03879
03880 char *idealstring=NULL;
03881 build_element_summary_string(pgdefinitions[i].summary, idealstring);
03882
03883
03884
03885 if (strcmp(idealstring, elementsummarystring)) {
03886 char buf[256] = { 0 };
03887 sprintf(buf, "Ideal elements (%5s): %s\n", pgdefinitions[i].name, idealstring);
03888 msgInfo << buf << sendmsg;
03889 sprintf(buf, "Found elements (%5s): %s\n", pointgroupname, elementsummarystring);
03890 msgInfo << buf << sendmsg;
03891 if (strlen(additionalelementstring))
03892 msgInfo << "Additional elements: " << additionalelementstring << sendmsg;
03893 if (strlen(missingelementstring))
03894 msgInfo << "Missing elements: " << missingelementstring << sendmsg;
03895 }
03896 delete [] idealstring;
03897 }
03898
03899
03900
03901
03902
03903 void Symmetry::draw_transformed_mol(Matrix4 rot) {
03904 int i;
03905 Molecule *mol = mlist->mol_from_id(sel->molid());
03906 MoleculeGraphics *gmol = mol->moleculeGraphics();
03907 const float *pos = sel->coordinates(mlist);
03908 for (i=0; i<sel->num_atoms; i++) {
03909 switch (mol->atom(i)->atomicnumber) {
03910 case 1:
03911 gmol->use_color(8);
03912 break;
03913 case 6:
03914 gmol->use_color(10);
03915 break;
03916 case 7:
03917 gmol->use_color(0);
03918 break;
03919 case 8:
03920 gmol->use_color(1);
03921 break;
03922 default:
03923 gmol->use_color(2);
03924 break;
03925 }
03926 float p[3];
03927 rot.multpoint3d(pos+3L*i, p);
03928 gmol->add_sphere(p, 2*sigma, 16);
03929 char tmp[32];
03930 sprintf(tmp, " %i", i);
03931 gmol->add_text(p, tmp, 1, 1.0f);
03932 }
03933 }
03934
03935
03936
03937
03938 #if 0
03939 static inline bool coplanar(const float *normal1, const float *normal2, float tol) {
03940 return collinear(normal1, normal2, tol);
03941 }
03942 #endif
03943
03944 static inline bool collinear(const float *axis1, const float *axis2, float tol) {
03945 if (fabs(dot_prod(axis1, axis2)) > tol) return 1;
03946 return 0;
03947 }
03948
03949 static inline bool orthogonal(const float *axis1, const float *axis2, float tol) {
03950 if (fabs(dot_prod(axis1, axis2)) < tol) return 1;
03951 return 0;
03952 }
03953
03954 static inline bool behind_plane(const float *normal, const float *coor) {
03955 return (dot_prod(normal, coor)>0.01);
03956 }
03957
03958
03959 static void align_plane_with_axis(const float *normal, const float *axis, float *alignedplane) {
03960 float inplane[3];
03961 cross_prod(inplane, axis, normal);
03962 cross_prod(alignedplane, inplane, axis);
03963 vec_normalize(alignedplane);
03964 }
03965
03966
03967
03968 static int isprime(int x) {
03969 int i;
03970 for (i=2; i<x; i++) {
03971 if (!(x%i)) return 0;
03972 }
03973 return 1;
03974 }
03975
03976
03977
03978 static int numprimefactors(int x) {
03979 int i, numfac=0;
03980 for (i=2; i<=x; i++) {
03981 if (!isprime(i)) continue;
03982 if (!(x%i)) {
03983 x /= i;
03984 numfac++;
03985 i--;
03986 }
03987 }
03988 return numfac;
03989 }
03990
03991
03992 inline int Symmetry::find_collinear_axis(const float *myaxis) {
03993 int i;
03994 for (i=0; i<axes.num(); i++) {
03995 if (collinear(myaxis, axes[i].v, collintol)) {
03996 return i;
03997 }
03998 }
03999 return -1;
04000 }
04001
04003 inline int Symmetry::plane_exists(const float *myplane) {
04004 int i;
04005 for (i=0; i<planes.num(); i++) {
04006 if (collinear(myplane, planes[i].v, collintol)) {
04007 return i;
04008 }
04009 }
04010 return -1;
04011 }
04012
04013
04014
04015
04016
04017 int Symmetry::is_planar(const float *normal) {
04018 Matrix4 alignx;
04019 alignx.transvecinv(normal[0], normal[1], normal[2]);
04020 int j, k;
04021 float xmin=0.0, xmax=0.0;
04022 for (j=0; j<sel->selected; j++) {
04023 float tmpcoor[3];
04024 alignx.multpoint3d(coor+3L*j, tmpcoor);
04025 xmin = tmpcoor[0];
04026 xmax = tmpcoor[0];
04027 break;
04028 }
04029
04030 for (k=j+1; k<sel->selected; k++) {
04031 float tmpcoor[3];
04032 alignx.multpoint3d(coor+3L*k, tmpcoor);
04033 if (tmpcoor[0]<xmin) xmin = tmpcoor[0];
04034 else if (tmpcoor[0]>xmax) xmax = tmpcoor[0];
04035 }
04036
04037 if (xmax-xmin < 1.5*sigma) return 1;
04038
04039 return 0;
04040 }
04041
04042
04043
04044 void Symmetry::assign_bonds() {
04045 Molecule *mol = mlist->mol_from_id(sel->molid());
04046
04047
04048
04049 int i, j=0;
04050 for (i=0; i<sel->num_atoms; i++) {
04051 if (sel->on[i]) atomindex[j++] = i;
04052 }
04053
04054 for (i=0; i<sel->selected; i++) {
04055 int j = atomindex[i];
04056
04057 bondsperatom[i].numbonds = 0;
04058 if (bondsperatom[i].bondto) delete [] bondsperatom[i].bondto;
04059 if (bondsperatom[i].length) delete [] bondsperatom[i].length;
04060 bondsperatom[i].bondto = new int[mol->atom(j)->bonds];
04061 bondsperatom[i].length = new float[mol->atom(j)->bonds];
04062
04063 int k;
04064 for (k=0; k<mol->atom(j)->bonds; k++) {
04065 int bondto = mol->atom(j)->bondTo[k];
04066
04067
04068 if (!sel->on[bondto]) continue;
04069
04070 float order = mol->getbondorder(j, k);
04071 if (order<0.f) order = 1.f;
04072
04073 if (bondto > j) {
04074
04075 int m;
04076 int found = 0;
04077 for (m=i+1; m<sel->selected; m++) {
04078 if (atomindex[m]==bondto) { found = 1; break; }
04079 }
04080
04081
04082
04083 if (found) {
04084
04085 Bond b;
04086 b.atom0 = i;
04087 b.atom1 = m;
04088 b.order = order;
04089
04090 b.length = distance(coor+3L*i, coor+3L*m);
04091 bonds.append(b);
04092 }
04093 }
04094 }
04095 }
04096
04097 for (i=0; i<bonds.num(); i++) {
04098 if (verbose>3) {
04099 char buf[256] = { 0 };
04100 sprintf(buf, "Bond %d: %d--%d %.1f L=%.2f", i,
04101 atomindex[bonds[i].atom0],
04102 atomindex[bonds[i].atom1],
04103 bonds[i].order, bonds[i].length);
04104 msgInfo << buf << sendmsg;
04105 }
04106 int numbonds;
04107 numbonds = bondsperatom[bonds[i].atom0].numbonds;
04108 bondsperatom[bonds[i].atom0].bondto[numbonds] = bonds[i].atom1;
04109 bondsperatom[bonds[i].atom0].length[numbonds] = bonds[i].length;
04110 bondsperatom[bonds[i].atom0].numbonds++;
04111
04112 numbonds = bondsperatom[bonds[i].atom1].numbonds;
04113 bondsperatom[bonds[i].atom1].bondto[numbonds] = bonds[i].atom0;
04114 bondsperatom[bonds[i].atom1].length[numbonds] = bonds[i].length;
04115 bondsperatom[bonds[i].atom1].numbonds++;
04116 }
04117 }
04118
04119
04120
04121
04122
04123
04124
04125
04126
04127
04128 void Symmetry::assign_bondvectors() {
04129 Molecule *mol = mlist->mol_from_id(sel->molid());
04130 memset(bondsum, 0, 3L*sel->selected*sizeof(float));
04131 int i;
04132 for (i=0; i<sel->selected; i++) {
04133 int k;
04134 for (k=0; k<bondsperatom[i].numbonds; k++) {
04135 int bondto = bondsperatom[i].bondto[k];
04136
04137
04138 if (!sel->on[bondto]) continue;
04139
04140 int j = atomindex[i];
04141 float order = mol->getbondorder(j, k);
04142 if (order<0.f) order = 1.f;
04143
04144 float bondvec[3];
04145 vec_sub(bondvec, coor+3L*bondto, coor+3L*i);
04146 vec_normalize(bondvec);
04147 vec_scaled_add(bondsum+3L*i, order, bondvec);
04148 }
04149 }
04150 }
04151
04152
04153
04154
04155
04156 static void assign_atoms(AtomSel *sel, MoleculeList *mlist, float *(&mycoor), int *(&atomtype)) {
04157
04158 const float *allcoor = sel->coordinates(mlist);
04159
04160 Molecule *mol = mlist->mol_from_id(sel->molid());
04161
04162
04163 char **typestringptr = new char*[sel->selected];
04164 int numtypes = 0;
04165
04166 int i, j=0;
04167 for (i=0; i<sel->num_atoms; i++) {
04168 if (!sel->on[i]) continue;
04169
04170
04171 vec_copy(mycoor+3L*j, allcoor+3L*i);
04172
04173
04174
04175 int k;
04176 int minatomicnum = 999;
04177 int maxatomicnum = 0;
04178 for (k=0; k<mol->atom(i)->bonds; k++) {
04179 int bondto = mol->atom(i)->bondTo[k];
04180 int atomicnum = mol->atom(bondto)->atomicnumber;
04181 if (atomicnum<minatomicnum) minatomicnum = atomicnum;
04182 if (atomicnum>maxatomicnum) maxatomicnum = atomicnum;
04183 }
04184
04185
04186
04187
04188 char *typestring = new char[8L+12L*mol->atom(i)->bonds];
04189 typestring[0] = '\0';
04190 char buf[100] = { 0 };
04191
04192
04193 sprintf(buf, "%i %i ", mol->atom(i)->atomicnumber, mol->atom(i)->bonds);
04194 strcat(typestring, buf);
04195
04196
04197
04198
04199 int m;
04200 for (m=minatomicnum; m<=maxatomicnum; m++) {
04201 unsigned char bondorder[7];
04202 memset(bondorder, 0, 7L*sizeof(unsigned char));
04203
04204 unsigned char bondedatomicnum = 0;
04205 for (k=0; k<mol->atom(i)->bonds; k++) {
04206 if (m == mol->atom(mol->atom(i)->bondTo[k])->atomicnumber) {
04207 bondedatomicnum++;
04208 float bo = mol->getbondorder(i, k);
04209 if (bo<0.f) bo = 1.f;
04210 (bondorder[(long)(2L*bo)])++;
04211 }
04212 }
04213 for (k=0; k<7; k++) {
04214 if (bondorder[k]) {
04215 sprintf(buf, "%i*(%i)%i ", bondorder[k], k, m);
04216 strcat(typestring, buf);
04217 }
04218 }
04219 }
04220
04221
04222
04223 int found = 0;
04224 for (k=0; k<numtypes; k++) {
04225 if (!strcmp(typestringptr[k], typestring)) {
04226 atomtype[j] = k;
04227 found = 1;
04228 delete [] typestring;
04229 break;
04230 }
04231 }
04232 if (!found) {
04233 atomtype[j] = numtypes;
04234 typestringptr[numtypes++] = typestring;
04235 }
04236
04237
04238 j++;
04239 }
04240
04241 for (i=0; i<numtypes; i++) {
04242 delete [] typestringptr[i];
04243 }
04244 delete [] typestringptr;
04245 }
04246
04247
04248
04249
04250
04251
04252
04253
04254
04255
04256
04257 inline static float trans_overlap(int *atomtype, float *(&coor), int numcoor,
04258 const Matrix4 *trans, float sigma,
04259 bool skipident, int maxnatoms) {
04260 float overlappermatch;
04261 return trans_overlap(atomtype, coor, numcoor, trans, sigma, skipident, maxnatoms, overlappermatch);
04262 }
04263
04264 static float trans_overlap(int *atomtype, float *(&coor), int numcoor,
04265 const Matrix4 *trans, float sigma,
04266 bool skipident, int maxnatoms, float &overlappermatch) {
04267
04268 const float *posA;
04269 posA = coor;
04270
04271 int *flgs = new int[numcoor];
04272 float *posB = new float[3L*numcoor];
04273 memset(flgs, 0, numcoor*sizeof(int));
04274
04275
04276 int i, ncompare=0;
04277 for(i=0; i<numcoor; i++) {
04278 trans->multpoint3d(posA+3L*i, posB+3L*i);
04279
04280
04281
04282 if (!(skipident && distance(posA+3L*i, posB+3L*i) < sigma)) {
04283 flgs[i] = 1;
04284 ncompare++;
04285 }
04286 }
04287
04288 if (ncompare<0.5*numcoor) {
04289
04290 delete [] flgs;
04291 delete [] posB;
04292 return 0.0;
04293 }
04294
04295
04296
04297 float dist;
04298
04299 dist = 3*sigma;
04300 float wrongelementpenalty = 100.0f/ncompare;
04301
04302 float overlap = 0.0;
04303 float antioverlap = 0.0;
04304 int i1, nmatches = 0, noverlap = 0, nwrongelement = 0;
04305 float maxr2=powf(1.0f*dist, 2);
04306 float itwosig2 = 1.0f/(2.0f*sigma*sigma);
04307
04308
04309 for (i1=0; i1<numcoor; i1++) {
04310 if (!flgs[i1]) continue;
04311
04312 float minr2 = maxr2+1.0f;
04313 float r2 = 0.0f;
04314
04315
04316 int j, i2=-1;
04317 for (j=0; j<numcoor; j++) {
04318 if (!flgs[j]) continue;
04319
04320 r2 = distance2(posA+3L*i1, posB+3L*j);
04321
04322 if (r2<minr2) { minr2 = r2; i2 = j; }
04323 }
04324
04325
04326 if (minr2<maxr2) {
04327 noverlap++;
04328
04329
04330 if (atomtype[i1]==atomtype[i2]) {
04331
04332 overlap += expf(-itwosig2*minr2);
04333 nmatches++;
04334 }
04335 else {
04336
04337 antioverlap += wrongelementpenalty*expf(-itwosig2*r2);
04338 nwrongelement++;
04339
04340 }
04341 }
04342 }
04343
04344 float nomatchpenalty = 0.0;
04345 overlappermatch = 0.0;
04346 int numnomatch = ncompare-nmatches;
04347
04348
04349
04350
04351
04352
04353
04354
04355
04356
04357 if (nmatches) overlappermatch = overlap/nmatches;
04358
04359 overlap -= antioverlap;
04360
04361 if (!(numnomatch==0)) {
04362
04363 nomatchpenalty = powf(overlappermatch, 5);
04364
04365 overlap -= 8*numnomatch*nomatchpenalty;
04366 }
04367 if (overlap<0) overlap = 0.0f;
04368
04369 overlap /= ncompare;
04370
04371
04372
04373
04374 delete [] flgs;
04375 delete [] posB;
04376
04377 return overlap;
04378 }
04379
04380
04381
04382
04383
04384
04385
04386
04387
04388
04389
04390
04391
04392
04393
04394
04395 int measure_trans_overlap(AtomSel *sel, MoleculeList *mlist, const Matrix4 *trans,
04396 float sigma, bool skipident, int maxnatoms, float &overlap) {
04397 if (!sel) return MEASURE_ERR_NOSEL;
04398 if (sel->num_atoms == 0) return MEASURE_ERR_NOATOMS;
04399
04400 float *coor = new float[3L*sel->selected];
04401
04402 int *atomtypes = new int[sel->selected];
04403 assign_atoms(sel, mlist, coor, atomtypes);
04404
04405 overlap = trans_overlap(atomtypes, coor, sel->selected, trans, sigma, skipident, maxnatoms);
04406
04407 return MEASURE_NOERR;
04408 }
04409
04410
04411
04412
04413
04414
04415
04416
04417 int measure_pointset_overlap(const float *posA, int natomsA, int *flagsA,
04418 const float *posB, int natomsB, int *flagsB,
04419 float sigma, float pairdist, float &overlap) {
04420
04421 int nsmall = natomsA<natomsB ? natomsA : natomsB;
04422
04423 int maxpairs = -1;
04424 GridSearchPair *pairlist, *p;
04425 pairlist = vmd_gridsearch3(posA, natomsA, flagsA, posB, natomsB, flagsB, pairdist,
04426 1, maxpairs);
04427
04428 overlap = 0.0;
04429 int i1, i2;
04430 float dx, dy, dz, r2, itwosig2 = 1.0f/(2.0f*sigma*sigma);
04431 for (p=pairlist; p; p=p->next) {
04432 i1 = p->ind1;
04433 i2 = p->ind2;
04434 dx = posA[3L*i1 ]-posB[3L*i2];
04435 dy = posA[3L*i1+1]-posB[3L*i2+1];
04436 dz = posA[3L*i1+2]-posB[3L*i2+2];
04437 r2 = dx*dx + dy*dy + dz*dz;
04438 overlap += expf(-itwosig2*r2);
04439 }
04440
04441 overlap /= nsmall;
04442
04443 return MEASURE_NOERR;
04444 }
04445