MeasureSymmetry.C

00001 /***************************************************************************
00002 *cr 
00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 
00004 *cr University of Illinois 
00005 *cr All Rights Reserved 
00006 *cr 
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010 * RCS INFORMATION:
00011 *
00012 * $RCSfile: MeasureSymmetry.C,v $
00013 * $Author: johns $ $Locker: $ $State: Exp $
00014 * $Revision: 1.67 $ $Date: 2020年07月08日 04:30:17 $
00015 *
00016 ***************************************************************************
00017 * DESCRIPTION:
00018 * The Symmetry class is the work horse behind the "measure symmetry"
00019 * command which guesses the pointgroup of a selection and returns the
00020 * according symmetry elements (mirror planes rotary axes, rotary reflection
00021 * axes). The algorithm is is fairly forgiving about molecules where atoms
00022 * are perturbed from the ideal position and tries its best to guess the
00023 * correct point group anyway. 
00024 * The tolerance can be controlled with the sigma parameter which is
00025 * the average allowed deviation from the ideal position.
00026 * Works nice on my 30 non-patholocical test cases. A pathological case
00027 * would for instance be a system with more than only a few atoms (say 15)
00028 * where only one atom distinguishes between two point groups.
00029 * If your selection contains more than a certain number of atoms then
00030 * only that much randomly chosen atoms are used to find planes and axes
00031 * in order to save time.
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" // for vmd_gridsearch3()
00044 #include "MoleculeGraphics.h" // needed only for debugging
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 // The bondsum is the sum of all normalized bond vectors
00073 // for an atom. If the endpoints of the bondsum vectors
00074 // of two atoms are further apart than BONDSUMTOL, then
00075 // The two atoms are not considered images of each other
00076 // with respect to the transformation that was tested.
00077 #define BONDSUMTOL 1.5
00078 
00079 // Minimum atom distance before the sanity check for
00080 // idealized coordinates fails.
00081 #define MINATOMDIST 0.6
00082 
00083 // Mirror plane classification
00084 #define VERTICALPLANE 1
00085 #define DIHEDRALPLANE 3
00086 #define HORIZONTALPLANE 4
00087 
00088 // Special axis order values
00089 #define INFINITE_ORDER -1
00090 #define PRELIMINARY_C2 -2
00091 
00092 // Flags for requesting special geometries for idealize_angle()
00093 #define TETRAHEDRON -4
00094 #define OCTAHEDRON -8
00095 #define DODECAHEDRON -12
00096 
00097 
00098 
00099 #if defined(_MSC_VER)
00100 // Microsoft's compiler lacks erfc() so we make our own here:
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 // Forward declarations of helper functions
00120 // static inline bool coplanar(const float *normal1, const float *normal2, float tol);
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 // Store the symmetry characteristics for a structure
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 // Symmetry elements defining a pointgroup
00154 typedef struct {
00155 char name[6]; // name string
00156 ElementSummary summary; 
00157 } PointGroupDefinition;
00158 
00159 PointGroupDefinition pgdefinitions[] = {
00160 // inv sig Cinf C2 C3 C4 C5 C6 C7 C8 S3 S8 S16
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 // Copy coordinates of selected atoms into local array and assign
00280 // and atomtypes based on chemial element and topology.
00281 assign_atoms(sel, mlist, coor, atomtype);
00282 
00283 // Determine the bond topology
00284 assign_bonds();
00285 
00286 
00287 // default collinearity tolerance ~0.9848
00288 collintol = cosf(float(DEGTORAD(10.0f)));
00289 
00290 // default coplanarity tolerance ~0.1736
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 // Impose certain symmetry elements on structure
00435 // by wrapping coordinates around and averaging them.
00436 // This creates symmetry, if there was one before or not.
00437 // If mutually incompatible elements are specified the
00438 // result is hard to predict:
00439 // Planes are idealized first, then rotary axes, then
00440 // rotary reflections and last an inversion.
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 // Abuse measure_inertia() to update the center of mass
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 //printf("rcom={%.2f %.2f %.2f}\n", rcom[0], rcom[1], rcom[2]);
00512 
00513 // Assign the bondsum vectors
00514 assign_bondvectors();
00515 
00516 
00517 for (i=0; i<2; i++) {
00518 
00519 // During idealization the coordinates are wrapped and
00520 // averaged.
00521 idealize_coordinates();
00522 
00523 // Use improved coordinates
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 // We don't know beforehand what is the ideal (sigma) tolerance
00552 // value to use in oder to find the right point group. If we
00553 // choose it too low we might find only a subgroup of the real 
00554 // pointgroup (e.g. Cs instead of D2h). However, if we choose
00555 // it too high we will find impossible combinations of symmetry
00556 // elements caused by false positives.
00557 // Thus we slowly ramp up the tolerance in steps of 0.05 until
00558 // we reach the user-specified maximum. This has some advantages:
00559 // After each step where we find a new non-C1 symmetry we'll
00560 // idealize the coordinatesand get a more symmetric structure
00561 // even if we just found a subgroup of the real one.
00562 // At some point we will find the higher symmetry point group.
00563 // But how do we know when to stop?
00564 // After each step the sanity of the idealized coordinates is
00565 // checked. Wrong point groups (i.e. wrong symmetry elements)
00566 // will screw up the idealized coordinates and bondsums which
00567 // can be easily detected. When this occcurs our search has
00568 // finished and we go back to the last sane symmetry.
00569 // The step where this symmetry occured the first time is
00570 // an indicator of the quality of the guess. The earlier we
00571 // found the right point group the less we had to bend the
00572 // coodinates to make them symmetric and the more trustworthy
00573 // is the guess.
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 // Self consistent iteration of symmetry element search
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 // If the idealized coordinates are not sane, i.e.
00604 // there are atoms closer than MINATOMDIST we stop.
00605 // Coordinates typically collapse during idealization
00606 // when false symmetry elements were found.
00607 // We can stop the search because with a higher tolerance
00608 // in the next step we are going to find the same bad
00609 // symmetry element again (and maybe even more).
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 // Determine the point group from symmetry elements
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 //msgInfo << "Level: " << pointgroup_rank(pointgroup, pointgrouporder)
00632 // << sendmsg;
00633 }
00634 
00635 // Try to decide wether the current step yielded a
00636 // better guess:
00637 // Whenever we detect a new pointgroup we have to
00638 // compare the ranks in order to determine which one
00639 // has higher symmetry. Often we find a subgroup
00640 // of the actual point group first. When the tolerance 
00641 // has increased we find the real point group which
00642 // will have a higher rank than the subgroup.
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 // Reset sigma and idealized coordinates to the values of 
00661 // the best step and recompute the symmetry elements.
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 // Determine the point group from symmetry elements
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 // Sort planes in the order horizontal, dihedral, vertical, rest
00690 sort_planes();
00691 
00692 // Determine the unique coordinates of the molecule
00693 unique_coordinates();
00694 
00695 // Try to get a transformation matrix that orients the
00696 // molecule according to the GAMESS standard orientation.
00697 orient_molecule();
00698 
00699 
00700 wkf_timer_destroy(timer);
00701 
00702 return MEASURE_NOERR;
00703 }
00704 
00705 
00706 // Self consistent iteration of symmetry element search
00707 // ----------------------------------------------------
00708 // If we would base our point group guess on elements that were
00709 // found in a single pass, then the biggest problem is to decide
00710 // which symmetry elements to keep and which to discard, based
00711 // on their overlap score.
00712 // While a point group guess, relying only on certain key
00713 // features, can be quite good, the number of found elements would
00714 // often be wrong. In case additional elements were found this
00715 // could totally screw up the idealized coordinates.
00716 // To avoid these problems we are iterate the element guess.
00717 // In each iteration we only keep the elements with a very good
00718 // score. This ensures that we don't get false positives. We then
00719 // idealize the symmetry elements which means we correct the
00720 // angles between then (i.e. an 87deg angle between two planes
00721 // becomes 90deg). Now all atoms will be wrapped around each
00722 // symmetry element and their positions averaged with their images
00723 // which gives us an updated set of idealized coordinates that will
00724 // be "more symmetric" than the coordinates from the previous
00725 // generation. In the next pass we use these improved coordinates
00726 // for guessing the elements and will probably find a few more this
00727 // time. This is repeated until the number of found symmetry elements
00728 // converges.
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 // In each iteration we start with a clean slate
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 // Reassign the bondsum vectors
00760 assign_bondvectors();
00761 
00762 // Determine primary axes of inertia
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 // Find all rotary axes, mirror planes and rotary reflection axes
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 // Now we have a set of symmetry elements but since they were
00798 // constructed from the possibly inacccurate molecular coordinates
00799 // their relationship is not ideal. For example the angle between
00800 // mirror planes parallel to a 3-fold rotation axis will only
00801 // approximately be 120 degree. We have to idealize the symmetry
00802 // elements so that we have a basis to idealize the coordinates
00803 // and determine the symmetry unique atoms.
00804 idealize_elements();
00805 
00806 // Generate symmetricized coordinates by wrapping atoms around
00807 // the symmetry elements and average with these images with the 
00808 // original positions.
00809 idealize_coordinates();
00810 
00811 // Should check the sanity of the idealized coordinates
00812 // here. I.e. check if there are atoms closer than, 0.6 A
00813 // and check if the bondsums match with the original geometry.
00814 
00815 // Compute RMSD between original and idealized coordinates
00816 rmsd = ideal_coordinate_rmsd();
00817 
00818 // Update the symmetry element statistics
00819 build_element_summary();
00820 
00821 // Create a string representation of the summary
00822 // Example: (inv) 2*(sigma) (Cinf) (C2)
00823 build_element_summary_string(elementsummary, elementsummarystring);
00824 
00825 // Iterations have converged when the rmsd difference between this and
00826 // the last structure is below a threshold and the symmetry elements
00827 // haven't changed
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 // Use improved coordinates for the next pass
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 // Find all mirror planes, rotary axes, rotary reflections.
00857 void Symmetry::find_symmetry_elements() {
00858 
00859 // ----- Planes -----------------
00860 
00861 // Check if the unique primary axes of inertia correspond to rotary axes
00862 find_elements_from_inertia();
00863 
00864 // Find all mirror planes
00865 find_planes();
00866 
00867 // Remove planes with too bad overlap
00868 purge_planes(0.5);
00869 
00870 // Are there horizontal/vertical/dihedral mirror planes?
00871 classify_planes();
00872 
00873 // Check for an inversion center
00874 check_add_inversion();
00875 
00876 
00877 // ------ Axes ------------------
00878 
00879 // Find all rotary axes that are intersections of planes.
00880 find_axes_from_plane_intersections();
00881 
00882 // Determine the order of each axis
00883 assign_axis_orders();
00884 
00885 // Some C2 axes cannot be found from principle axes of inertia
00886 // or through plane intersection. We must get these from the
00887 // atomic coordinates.
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 // Sort axes by decreasing order
00903 sort_axes();
00904 
00905 
00906 // ----- Rotary reflections -----
00907 
00908 // Since all rotary reflections are collinear with a rotation axis
00909 // we can simply go through the list of axes and check if which
00910 // axes correspond to improper rotations.
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 // Normalize the summed up rotary reflection vectors
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 // Remove axes with too bad overlap
00927 purge_axes(0.7f);
00928 purge_rotreflections(0.75f);
00929 
00930 // Must classify planes again because of dihedral planes depending
00931 // on the new axes:
00932 classify_planes();
00933 
00934 // Classify perpendicular axes
00935 classify_axes();
00936 }
00937 
00938 
00939 // Determine principle axes of inertia, check them for uniqueness
00940 // by comparing the according eigenvalues. Then, see if the unique
00941 // axes correspond to rotary axes or rotary reflections and, in
00942 // case, add axes to the list and assing the order.
00943 // We also check if there are horizontal mirror planes for the
00944 // principal axes.
00945 void Symmetry::find_elements_from_inertia() {
00946 int i;
00947 int numuniqueprimary = 0;
00948 memset(uniqueprimary, 0, 3L*sizeof(int));
00949 
00950 // Normalize eigenvalues
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 // Check if the molecule is linear
00958 if (fabs(eigenval[0]-eigenval[1])<0.05 && eigenval[2]<0.05) {
00959 linear = 1;
00960 }
00961 
00962 // Get provisional info about the rotational order of the inertia axes
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 //printf("primary axis of inertia %i: order=%i overlap=%.2f\n", i, order[i], overlap[i]);
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)); // empirical;
00977 
00978 // If all three moments of inertia are the same, the molecule is a spherical top,
00979 // i.e. the possible point groups are T, Td, Th, O, Oh, I, Ih.
00980 if (fabs(eigenval[0]-eigenval[1])<tol &&
00981 fabs(eigenval[1]-eigenval[2])<tol) {
00982 
00983 sphericaltop = 1;
00984 
00985 } else {
00986 // If two moments of inertia are the same, the molecule is a
00987 // symmetric top, i.e. the possible points groups are C1, Cn,
00988 // Cnv, Cnh, Sn, Dn, Dnh, Dnd.
00989 if (fabs(eigenval[0]-eigenval[1])<tol ||
00990 fabs(eigenval[1]-eigenval[2])<tol ||
00991 fabs(eigenval[0]-eigenval[2])<tol) {
00992 
00993 // Determine which is the unique primary axis. Also make sure
00994 // that this axis corresponds to a Cn rotation, otherwise we
00995 // have a false positive, i.e. two eigenvalues are identical
00996 // or similar only by incident.
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 // If the molecule is planar (but not linear) with the current
01031 // primary axis of inertia corresaponding to the plane normal
01032 // then we want to add that plane to the list.
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 // Cinfv and Dinfh symmetry:
01047 // Use the unique primary axis of inertia as Cinf
01048 // For Dinfh an arbitrary C2 axis perpendicular to Cinf will be
01049 // automatically generated later by plane intersection.
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 // We have to add an arbitrary vertical plane.
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 //printf("primary axis of inertia %i: order=%i overlap=%.2f\n", i, order[i], overlap[i]);
01070 
01071 if (order[i] > 1 && overlap[i]>OVERLAPCUTOFF) {
01072 //printf("unique primary axes = %i planarmol = %i\n", i, planarmol);
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 // Add a possible horizontal mirror plane:
01083 // In case the molecule is planar the general plane finding algorithm
01084 // won't recognize the horizontal plane because identical atom transformations
01085 // are skipped. Thus we boost the weight here and provide it to 
01086 // check_add_plane().
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 // Find atoms with same distance to the center of mass and check if the 
01097 // plane that is orthogonal to their connecting vector is a mirror plane.
01098 void Symmetry::find_planes() {
01099 int i,j;
01100 float posA[3], posB[3];
01101 
01102 // Loop over all atoms
01103 for (i=0; i<sel->selected; i++) {
01104 
01105 // If we have more than maxnatoms atoms in the selection then pick 
01106 // only approximately maxnatoms random atoms for the comparison.
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 // Exclude atoms that are close to COM:
01115 // Closer than sigma is too close in any case, otherwise exclude
01116 // more the larger the molecule is.
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 // consider only pairs with identical atom types
01127 if (atomtype[i]==atomtype[j]) {
01128 
01129 // If the atoms are too close to the plane, the error gets
01130 // too big so we skip.
01131 float dist = distance(posA, posB);
01132 if (dist<0.25) continue;
01133 
01134 // We have found a hypothetical rotation axis or mirror plane
01135 
01136 // Subtracting the two position vectors yields a vector that
01137 // defines the normal of a plane that bisects the angle between
01138 // the two atoms, i.e. the plane is a potential mirror plane.
01139 float normal[3];
01140 vec_sub(normal, posA, posB);
01141 vec_normalize(normal);
01142 
01143 // Check plane and possibly add it to the list;
01144 check_add_plane(normal);
01145 }
01146 }
01147 }
01148 
01149 // Normalize the summed up normal vectors
01150 for (i=0; i<numplanes(); i++) {
01151 vec_normalize(planes[i].v);
01152 //printf("plane[%i]: weight=%3i, overlap=%.2f\n", i, planes[i].weight, planes[i].overlap);
01153 }
01154 }
01155 
01156 
01157 // Find C2 rotation axes from the vectors bisecting the angle between
01158 // two atoms and the center of mass.
01159 void Symmetry::find_C2axes() {
01160 int i,j;
01161 float posA[3], posB[3];
01162 float rA, rB;
01163 // Loop over all atoms
01164 for (i=0; i<sel->selected; i++) {
01165 // If we have more than maxnatoms atoms in the selection then pick 
01166 // only approximately maxnatoms random atoms for the comparison.
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 // Exclude atoms that are close to COM:
01175 // Closer than sigma is too close in any case, otherwise exclude
01176 // more the larger the molecule is.
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 // consider only pairs with identical atom types
01187 if (atomtype[i]==atomtype[j]) {
01188 // We have found a hypothetical rotation axis or mirror plane
01189 float alpha = angle(posA, posB);
01190 
01191 // See if the vector bisecting the angle between the two atoms
01192 // and the center of mass corresponds to a C2 axis. We must
01193 // exclude alpha==180 because in this case the bisection is
01194 // not uniquely defined.
01195 if (fabs(fabs(alpha)-180) > 10) {
01196 float testaxis[3];
01197 vec_add(testaxis, posA, posB);
01198 vec_normalize(testaxis);
01199 
01200 // Check axis and possibly add it to the list;
01201 if (!planes.num() || numverticalplanes) {
01202 //printf("Checking C2 axis\n"); 
01203 check_add_axis(testaxis, 2);
01204 }
01205 }
01206 }
01207 }
01208 }
01209 
01210 // Normalize the summed up axis vectors
01211 for (i=0; i<numaxes(); i++) {
01212 vec_normalize(axes[i].v);
01213 }
01214 }
01215 
01216 // Computes the factorial of n
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 // Generate all n!/(k!(n-k)!) combinations of k different
01226 // elements drawn from a total of n elements.
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 // Create the combinations (i.e. populate combolist)
01239 void create() { recurse(0, -1); }
01240 
01241 void recurse(int begin, int level);
01242 void print();
01243 
01244 // Get pointer to the i-th combination
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; // prevent overflow of float by fac(n)
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 // Recursive function to generate combinations
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 // Print the combinations
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 // Find Cn rotation axes from the atoms that lie in the
01293 // same plane. The plane normal defines the hypothetical
01294 // axis
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 // Loop over all atoms
01302 for (i=0; i<sel->selected; i++) {
01303 // If we have more than maxnatoms atoms in the selection then pick 
01304 // only approximately maxnatoms random atoms for the comparison.
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 // Exclude atoms that are close to COM:
01313 // Closer than sigma is too close in any case, otherwise exclude
01314 // more the larger the molecule is.
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 // If we have more than maxnatoms atoms in the selection
01322 // then pick only approximately maxnatoms random atoms
01323 // for the comparison.
01324 if (sel->selected > maxnatoms && vmd_random() % sel->selected > maxnatoms)
01325 continue;
01326 
01327 // Consider only pairs with identical atom types
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 // Atoms should have approximately same distance from COM
01335 if (fabs(rA-rB) > sigma) continue;
01336 
01337 atomtuple[ntup++] = j;
01338 
01339 // We don't consider more then 10 equivalent atoms
01340 // to save resources.
01341 if (ntup>=10) break;
01342 }
01343 if (ntup<order) continue;
01344 
01345 //printf("equiv. atoms: ");
01346 //for (j=0; j<ntup; j++) { printf("%d ",atomtuple[j]); }
01347 //printf("\n");
01348 
01349 // Generate all combinations of tuples with order different
01350 // elements drawn from a total of ntup elements.
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 //printf("combi %d: {", j);
01361 
01362 // Test wether all atoms in the tuple lie aproximately
01363 // within one plane
01364 int inplane = 1, anyinplane = 0;
01365 for (m=0; m<order; m++) {
01366 int atom = atomtuple[combo.get(j)[m]];
01367 //printf(" %d", atom); //combo.get(j)[m]);
01368 vec_sub(pos, coor+3L*atom, rcom);
01369 vec_add(testaxis, testaxis, pos);
01370 
01371 // Find a second atom from the tuple that encloses
01372 // an angle of >45 deg with the first one.
01373 // If we don't find one then we ignore the in-plane
01374 // testing for the current atom.
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 // Check if the cross product of the two atom positions
01390 // is collinear with the plane normal (which was summed
01391 // up from the cross products of the previous pairs).
01392 if (collinear(normal, cross, collintol)) {
01393 vec_add(normal, normal, cross);
01394 anyinplane = 1;
01395 } else {
01396 // This atom is not in plane with the others
01397 inplane = 0;
01398 break;
01399 }
01400 }
01401 
01402 //printf("}\n");
01403 
01404 // If the atoms of the current tuple aren't in one
01405 // plane then the positions cannot be used to define
01406 // a Cn axis and we skip it.
01407 if (!inplane || !anyinplane) continue;
01408 
01409 vec_normalize(normal);
01410 
01411 // Check axis and possibly add it to the list;
01412 printf("Checking C%d axis\n", order); 
01413 check_add_axis(normal, order);
01414 }
01415 
01416 }
01417 
01418 delete [] atomtuple;
01419 
01420 // Normalize the summed up axis vectors
01421 for (i=0; i<numaxes(); i++) {
01422 vec_normalize(axes[i].v);
01423 }
01424 }
01425 
01426 
01427 // Get rotation axes.
01428 // Each intersection of two mirror planes is itself a rotation
01429 // axis. Computing the intersection is easy in this case because
01430 // we know that any two mirror planes go through the center of
01431 // mass. The direction of the intersection line is just the 
01432 // cross product of the two plane normals.
01433 
01434 void Symmetry::find_axes_from_plane_intersections() {
01435 // If there aren't at least 2 planes we won't find any axes.
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 // Ignore intersections of parallel planes
01448 if (norm(newaxis)<0.05) continue;
01449 
01450 // Loop over already existing axes and check if the new axis is
01451 // collinear with one of them. In this case we add the two axes
01452 // in order to obtain an average (after normalizing at the end).
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 //printf("dot=% .4f fabs(dot)=%.4f, 1-collintol=%.2f\n", dot, fabs(dot), collintol);
01462 if (fabs(dot) > collintol) {
01463 // We are summing up the collinear axes to get the average
01464 // of the equivalent axes.
01465 if (dot>0) {
01466 vec_incr(axis(k), newaxis); // axes[k] += normal
01467 } else {
01468 vec_sub(axis(k), axis(k), newaxis); // axes[k] -= newaxis
01469 }
01470 axes[k].weight++;
01471 found = 1;
01472 break;
01473 }
01474 }
01475 
01476 if (!found) {
01477 // We found no existing collinear axis, so add a new one to the list.
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 // Normalize the summed up axis vectors
01493 for (i=0; i<numaxes(); i++) {
01494 vec_normalize(axis(i));
01495 }
01496 }
01497 
01498 
01499 // Determine which of the existing mirror planes are horizontal,
01500 // vertical, or dihedral planes with respect to the first axis.
01501 // The axes must have been sorted already so that the axes with
01502 // the highest order (primary axis) is the first one.
01503 // A horizontal plane is perpendicular to the primary axis.
01504 // A vertical plane includes the primary axis.
01505 // A dihedral plane is vertical to the corresponding axis and 
01506 // bisects the angle formed by a pair of C2 axes.
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 // Is there a plane perpendicular to the highest axis?
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 // The current plane is vertical to the first axis.
01531 numverticalplanes++;
01532 planes[i].type = VERTICALPLANE;
01533 
01534 // Now loop over pairs of C2 axes that are in that plane:
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 // Check if the plane bisects the pair of C2 axes
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(); // stop looping axis pairs
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(); // stop looping axis pairs
01561 break;
01562 }
01563 
01564 }
01565 }
01566 }
01567 }
01568 
01569 void Symmetry::classify_axes() {
01570 if (!numaxes()) return;
01571 
01572 // If n is the higest Cn axis order, are there n C2 axes
01573 // perpendicular to Cn?
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 // Check if center of mass is an inversion center.
01589 float Symmetry::score_inversion() {
01590 // Construct inversion matrix
01591 Matrix4 inv;
01592 inv.translate(rcom[0], rcom[1], rcom[2]);
01593 inv.scale(-1.0); // inversion
01594 inv.translate(-rcom[0], -rcom[1], -rcom[2]);
01595 
01596 // Return score 
01597 return trans_overlap(atomtype, coor, sel->selected, &inv, 1.5f*sigma,
01598 NOSKIP_IDENTICAL, maxnatoms);
01599 }
01600 
01601 // Check if center of mass is an inversion center
01602 // and set the inversion flag accordingly.
01603 void Symmetry::check_add_inversion() {
01604 // Verify inversion center
01605 float overlap = score_inversion();
01606 
01607 //printf("inversion overlap = %.2f\n", overlap);
01608 
01609 // We trust planes and axes more than the inversion since 
01610 // they are averages over multiple detected elements.
01611 // Hencewe make the cutoff range from 0.3 to 0.5, depending
01612 // on the number of other found symmetry elements.
01613 if (overlap>0.5-0.2/(1+(planes.num()*axes.num()))) {
01614 inversion = 1;
01615 }
01616 }
01617 
01618 // Check if the given normal represents a mirror plane.
01619 float Symmetry::score_plane(const float *normal) {
01620 // A mirror symmetry operation is an inversion + a rotation
01621 // by 180 deg about the normal of the mirror plane.
01622 // We also need to shift to origin and back.
01623 Matrix4 mirror;
01624 mirror.translate(rcom[0], rcom[1], rcom[2]);
01625 mirror.scale(-1.0); // inversion
01626 mirror.rotate_axis(normal, float(VMD_PI));
01627 mirror.translate(-rcom[0], -rcom[1], -rcom[2]);
01628 
01629 // Return score
01630 return trans_overlap(atomtype, coor, sel->selected, &mirror, 1.5f*sigma, NOSKIP_IDENTICAL, maxnatoms);
01631 }
01632 
01633 // Append the given plane to the list.
01634 // If a approximately coplanar plane exists already then add the new
01635 // plane normal to the existing one thus getting an average direction.
01636 // The input plane normal must have a length of 1.
01637 // Note: After the list of planes is complete you want to normalize all
01638 // plane vectors.
01639 void Symmetry::check_add_plane(const float *normal, int weight) {
01640 
01641 // Verify the mirror plane
01642 float overlap = score_plane(normal);
01643 
01644 // In case of significant overlap between the original and the mirror
01645 if (overlap<OVERLAPCUTOFF) return;
01646 
01647 
01648 // We have to loop over all existing planes and see if the current
01649 // plane is coplanar with one of them. If so, we add to that existing
01650 // one. The resulting plane normal vector can later be normalized so
01651 // that we get the average plane normal of the equivalent planes.
01652 // If no existing coplanar plane is found a new plane is added to the
01653 // array.
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 // If two planes are coplanar we can use their average.
01662 // The planes are parallel if cross(n1,n2)==0.
01663 // However it is faster to use 
01664 // |cross(n1,n2)|^2 = |n1|*|n2| - dot(n1,n2)^2
01665 // = 1 - dot(n1,n2)^2
01666 // The latter is true because the normals have length 1.
01667 // Hence |cross(n1,n2)| == 0 is equivalent to 
01668 // |dot(n1,n2)| == 1.
01669 float dot = dot_prod(avgnormal, normal);
01670 if (fabs(dot) > collintol) {
01671 
01672 // We are summing up the coplanar planes to get the average
01673 // of the equivalent planes.
01674 if (dot>0) { 
01675 vec_incr(plane(k), normal); // planes[k] += normal
01676 } else {
01677 vec_scaled_add(plane(k), -1, normal); // planes[k] -= 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 // We found no existing coplanar plane, so add a new one to the list.
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 // Check if vector testaxis defines a rotary axis of the given order.
01698 float Symmetry::score_axis(const float *testaxis, int order) {
01699 // Construct rotation matrix.
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 // Verify symmetry axis
01706 return trans_overlap(atomtype, coor, sel->selected, &rot, 2*sigma,
01707 NOSKIP_IDENTICAL, maxnatoms);
01708 }
01709 
01710 
01711 // Append a given axis to the list in case it is C2.
01712 // We look for C2 axes only, because these are the only ones
01713 // that cannot always be generated by intersecting mirror planes.
01714 // Examples are the three C2 axes perpendicular to the primary C3 in 
01715 // the D3 pointgroup.
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 // Preliminary axes (order<0) haven't been averaged and
01728 // collapsed yet. For these cases we are summing up the
01729 // collinear axes to get the average of the equivalent
01730 // axes.
01731 if (axes[k].order==-order) {
01732 if (dot>0) { 
01733 vec_incr(axis(k), testaxis); // axes[k] += testaxis
01734 } else {
01735 vec_scaled_add(axis(k), -1, testaxis); // axes[k] -= testaxis
01736 }
01737 axes[k].weight++;
01738 }
01739 found = 1;
01740 break;
01741 }
01742 }
01743 if (!found) {
01744 // We found no existing collinear axis, so add a new one to the list.
01745 float overlap = score_axis(testaxis, order);
01746 if (overlap>OVERLAPCUTOFF) {
01747 Axis a;
01748 vec_copy(a.v, testaxis);
01749 // We tag the order as preliminary C2 in order to distinguish it 
01750 // from C2 axes found by plane intersection. The latter have a higher
01751 // accuracy since they are based on averaged planes.
01752 a.order = -order; //PRELIMINARY_C2;
01753 a.overlap = overlap;
01754 a.weight = 1;
01755 a.type = 0;
01756 axes.append(a);
01757 }
01758 }
01759 }
01760 
01761 // Get the score for given rotary reflection.
01762 float Symmetry::score_rotary_reflection(const float *testaxis, int order) {
01763 // Construct transformation matrix. An n-fold rotary
01764 // reflection is a rotation by 360/n deg followed by
01765 // a reflection about the plane perpendicular to the axis.
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); // inversion
01770 rot.rotate_axis(testaxis, float(VMD_PI));
01771 rot.translate(-rcom[0], -rcom[1], -rcom[2]);
01772 
01773 // Return score
01774 return trans_overlap(atomtype, coor, sel->selected, &rot, sigma,
01775 NOSKIP_IDENTICAL, maxnatoms);
01776 }
01777 
01778 
01779 // testaxis must be normalized.
01780 void Symmetry::check_add_rotary_reflection(const float *testaxis, int maxorder) {
01781 if (maxorder<4) return;
01782 
01783 //msgInfo << "checking improper: maxorder=" << maxorder << sendmsg;
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 // Get the overlap for each rotary reflection:
01791 float overlap = score_rotary_reflection(testaxis, n);
01792 //printf("rotrefl: n=%i, axis angle = %.2f, overlap = %.2f\n", n, 360.0/n, overlap);
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 // We are summing up the collinear axes to get the average
01808 // of the equivalent axes.
01809 if (dot>0) { 
01810 vec_incr(rotreflect(k), testaxis); // axes[k] += testaxis
01811 } else {
01812 vec_scaled_add(rotreflect(k), -1, testaxis); // axes[k] -= testaxis
01813 }
01814 rotreflections[k].weight++;
01815 found = 1;
01816 break;
01817 }
01818 }
01819 if (!found) {
01820 // We found no existing collinear rr-axis, so add a new one to the list.
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 // Purge planes with an overlap below half of the maximum overlap
01835 // that occurs in planes and axes.
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 // keep this plane
01856 keep[i] = 1;
01857 }
01858 }
01859 
01860 prune_planes(keep);
01861 
01862 delete [] keep;
01863 }
01864 
01865 
01866 // Purge axes with an overlap below half of the maximum overlap
01867 // that occurs in planes and axes.
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 // keep this axis
01884 keep[i] = 1;
01885 }
01886 }
01887 
01888 prune_axes(keep);
01889 
01890 delete [] keep;
01891 }
01892 
01893 
01894 // Purge rotary reflections with an overlap below half of the maximum overlap
01895 // that occurs in planes and axes.
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 // keep this rotary reflection
01912 keep[i] = 1;
01913 }
01914 else if (verbose>4) {
01915 // purge this rotary reflection
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 // For each plane i prune from the list if keep[i]==0.
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 //printf("keep plane[%i] = %i\n", i, keep[i]);
01941 if (keep[i]) {
01942 // keep this plane
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 // purge this plane
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 // For each axis i prune from the list if keep[i]==0.
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 //printf("keep axis[%i] = %i\n", i, keep[i]);
01984 if (keep[i]) {
01985 // keep this axis
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 // For each axis i prune from the list if keep[i]==0.
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 // keep this rotary reflection
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 // Assign orders to rotational symmetry axes.
02026 void Symmetry::assign_axis_orders() {
02027 if (!numaxes()) return;
02028 
02029 maxaxisorder = axes[0].order;
02030 
02031 // Compute all axis orders
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 //printf("axis[%i] {%f %f %f} order = %i, weight = %i overlap = %.2f\n", i,
02037 // axes[i].v[0], axes[i].v[1], axes[i].v[2], axes[i].order, axes[i].weight, axes[i].overlap);
02038 if (axes[i].order>maxaxisorder) {
02039 maxaxisorder = axes[i].order;
02040 }
02041 }
02042 }
02043 
02044 
02045 // Assign orders to preliminary C2 rotational symmetry axes.
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 // We have computed the overlap only for
02053 // the first instance not for the average
02054 // so we want to re-score it here.
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 // Sort the axes according to decreasing order
02064 // Also eliminate C1 axes and Cn axes that are parallel
02065 // to a Cinf axes.
02066 void Symmetry::sort_axes() {
02067 int i;
02068 
02069 // count number of sorted axes.
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 // Sort the planes: horizontal plane first, then dihedral,
02106 // then vertical, then the rest
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 // Find the highest rotational symmetry order (up to 8) for the given axis.
02140 // Rotates the selection by 360/i degrees (i=2...8) and checks the structural
02141 // overlap.
02142 int Symmetry::axis_order(const float *axis, float *overlap) {
02143 int i;
02144 
02145 // Get the overlap for each rotation
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 // need to shift to origin and back
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 // Get the maximum overlap
02159 float maxover = overlaparray[2];
02160 //printf("orders: %.2f ", overlaparray[2]);
02161 for (i=3; i<=MAXORDERCN; i++) {
02162 //printf("%.2f ", overlaparray[i]);
02163 if (overlaparray[i]>maxover) {
02164 maxover = overlaparray[i];
02165 }
02166 }
02167 
02168 // If overlap for a certain rotation is greater than half maxover
02169 // we assume that the axis has the according order or a multiple of that.
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 //printf(" tol1=%.2f tol=%.2f\n", tol1, tol1-tol2+pow(0.09*maxorder,6));
02191 
02192 // Make the tolerance dependent on the maximum order.
02193 // For maxorder=2 tol=0.2 for maxorder=8 tol=0.27 (empirical)
02194 if (maxover<tol1+powf(0.09f*maxorder, 6)) {
02195 *overlap = maxover;
02196 return 1;
02197 }
02198 
02199 // We return the arithmetic mean of the relative and the maximum overlap.
02200 // The relative overlap is the quality of the matching axis orders wrt
02201 // the best non-matching one.
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 // Generate symmetricized coordinates by wrapping atoms around
02236 // the symmetry elements and average these images with the 
02237 // original positions.
02238 void Symmetry::idealize_coordinates() {
02239 int i;
02240 
02241 // copy atom coordinates into idealcoor array
02242 memcpy(idealcoor, coor, 3L*sel->selected*sizeof(float));
02243 
02244 // -----------Planes-------------
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); // inversion
02252 mirror.rotate_axis(planes[i].v, float(VMD_PI));
02253 mirror.translate(-rcom[0], -rcom[1], -rcom[2]);
02254 
02255 //printf("PLANE %d\n", i);
02256 if (average_coordinates(&mirror)) keep[i] = 1;
02257 }
02258 
02259 prune_planes(keep);
02260 
02261 
02262 // -----------Axes-----------
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 //printf("AXIS %i\n", i);
02275 
02276 // In case of a Cinf axis just use order 4 to idealize in the two
02277 // dimensions perpendicular to the axis
02278 if (order==INFINITE_ORDER) order = 4;
02279 
02280 // Average the coordinates over all rotary orientations
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 // For each atom find the closest transformed image atom with
02290 // the same chemical element. The array matchlist will contain
02291 // the closest match atom index (into the collapsed list of
02292 // selected atoms) for each selected atom.
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; // no match found;
02300 
02301 if (checkbonds) {
02302 // Rotate the bondsum vector according to the Cn axis
02303 // and compare
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 //printf("average coor for angle=%.2f\n", n*360.f/order);
02319 }
02320 
02321 if (success) {
02322 // Combine the averaged coordinates for this axis with the existing
02323 // ideal coordinates
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 // -----------Rotary reflections-----------
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); // inversion
02348 rot.rotate_axis(rotreflect(i), float(VMD_PI));
02349 rot.translate(-rcom[0], -rcom[1], -rcom[2]);
02350 
02351 //printf("ROTREF %i\n", i);
02352 if (average_coordinates(&rot)) keep[i] = 1;
02353 }
02354 
02355 prune_rotreflections(keep);
02356 
02357 delete [] keep;
02358 
02359 // -----------Inversion-----------
02360 if (inversion) {
02361 // Construct inversion matrix
02362 Matrix4 inv;
02363 inv.translate(rcom[0], rcom[1], rcom[2]);
02364 inv.scale(-1.0f); // inversion
02365 inv.translate(-rcom[0], -rcom[1], -rcom[2]);
02366 
02367 if (!average_coordinates(&inv)) inversion = 0;
02368 }
02369 }
02370 
02371 
02372 // Averages between original and transformed coordinates.
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; // no match found;
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; // no match found;
02404 
02405 // Average between original and image coordinates 
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 // Check the bondsum for atom j and its image m generated by 
02417 // transformation trans.
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 //printf("bondsum 1: {%.2f %.2f %.2f}\n", bondsum[3L*j], bondsum[3L*j+1], bondsum[3L*j+2]);
02433 //printf("bondsum 2: {%.2f %.2f %.2f}\n", bondsum[3L*m], bondsum[3L*m+1], bondsum[3L*m+2]);
02434 //printf("bondsum 1 transf: {%.2f %.2f %.2f}\n", tmp[0], tmp[1], tmp[2]);
02435 return 0;
02436 }
02437 return 1;
02438 }
02439 
02440 
02441 // For each atom find the closest transformed image atom
02442 // with the same chemical element. The result is an array
02443 // containing the closest match atom index for each
02444 // selected atom.
02445 // XXX maybe put bondsum check in here?!
02446 void Symmetry::identify_transformed_atoms(Matrix4 *trans,
02447 int &nmatches,
02448 int *(&matchlist)) {
02449 // get atom coordinates
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 // generate transformed coordinates
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 // consider only pairs with identical atom types 
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 //printf("atom %i matches %i (atomtype %i)\n",
02492 // atomindex[i], atomindex[kmatch], atomtype[i]);
02493 }
02494 else {
02495 // no match found
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 // Compute the RMSD between original and idealized coordinates
02511 float Symmetry::ideal_coordinate_rmsd () {
02512 // get original atom coordinates
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 // Compute the RMSD between original and idealized coordinates
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 // Idealize all symmetry elements so that they have certain geometries
02543 // with respect to each other. This means, for instance, to make sure
02544 // that the planes vertical to a C3 axis are evenly spaced by 60 degree
02545 // angles and that their intersection is exactly the C3 axis.
02546 //
02547 // We start with idealizing horizontal and vertical planes with respect
02548 // to the primary axis. If there is no unique primary axis (e.g. for
02549 // Oh, Th, Ih) we take each of the highest order axes, idealize it wrt
02550 // to the first one and then idealize the planes vertical and horizontal
02551 // to them. Since all plajnes should be idealized by now we next
02552 // idealize axes that are intersections of two planes. Then we adjust
02553 // axes that are ideally bisectors of two planes sch as they occur in
02554 // Dnd point groups. Finally we idealize axes that are not related to
02555 // any planes such as in Dn groups. They should be perpedicular to the
02556 // highest order axis and evenly spaced around it.
02557 // Rotary reflections are simply brought to collinearity with their
02558 // accompanying rotary axis (An improper axis never comes alone, e.g.
02559 // S4 implies a C2 axis).
02560 // Symmetry elements that could not be idealized are removed from the
02561 // list.
02562 void Symmetry::idealize_elements() {
02563 int i;
02564 
02565 // Array of flags indicating if a plane has been idealized
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 // If there are axes present we will use the first axis as reference.
02574 // Since the axes were sorted by axis order before our reference will
02575 // be the axis with the highest order (or one of them).
02576 // Then align any horizontal and vertical planes with that axis.
02577 if (axes.num()) {
02578 
02579 // Make horizontal refplane truly perpendicular to reference axis
02580 if (horizontalplane>=0) {
02581 //msgInfo << "Idealizing horizontal plane " << horizontalplane << ", numplanes="
02582 // << planes.num() << sendmsg;
02583 vec_copy(planes[horizontalplane].v, axes[0].v);
02584 idealp[horizontalplane] = 1;
02585 } 
02586 
02587 if (numverticalplanes) {
02588 // Find the first vertical plane and align it with the reference axis so it
02589 // becomes our vertical reference plane.
02590 int vref = -1;
02591 for (i=0; i<planes.num(); i++) {
02592 if (planes[i].type&VERTICALPLANE) {
02593 //msgInfo << "Idealizing vertical plane " << i << sendmsg;
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 // Find other vertical planes and idealize them wrt the vertical reference plane.
02605 // To do so, we first have to align the found vertical plane with the reference
02606 // axis and then idealize the angle wrt the vertical reference plane.
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 //msgInfo << "Idealizing plane " << i << " wrt vertical plane " << vref << sendmsg;
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 // Equivalent plane exists already.
02622 // Since idealp[i] will not be set 1, this plane will be deleted
02623 // at the end of this function.
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 { // No axis present
02638 
02639 // If we have only one plane, it is automatically ideal.
02640 if (planes.num()<=1) {
02641 if (idealp) delete [] idealp;
02642 return;
02643 }
02644 
02645 // Actually, if there is more than one plane, at least one axis
02646 // should have been found by intersection. If we don't have it
02647 // here, it was probably purged in find_symmetry_elements().
02648 // That means that at least one of the planes is most likely bad.
02649 // We will keep the only best plane and mark it as ideal.
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 //msgInfo << "Idealize axis " << i << " (C" << axes[i].order << ")" << sendmsg;
02695 
02696 // Find plane that contains both axes[0] and axis[i]
02697 // and idealize axis[i] wrt reference axis[0]
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 // The plane couldn't be idealized in the previous step.
02703 // Skip it.
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 // Find planes that are vertical to the current axis and
02727 // idealize them wrt plane vref.
02728 for (j=0; j<planes.num(); j++) {
02729 if (idealp[j]) continue;
02730 
02731 // Check if plane is vertical to axis[i]; 
02732 if (orthogonal(planes[j].v, axes[i].v, orthogtol)) {
02733 // align the current plane with the idealized axis
02734 align_plane_with_axis(planes[j].v, axes[i].v, planes[j].v);
02735 //msgInfo << "Idealizing plane " <<j<< " wrt vertical plane " <<vref<< sendmsg;
02736 
02737 // idealize vertical plane wrt vref
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 // Equivalent plane exists already.
02748 // Since idealp[i] will not be set 1, this plane will be deleted
02749 // below.
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 // Delete all planes that could not be idealized.
02763 prune_planes(idealp);
02764 
02765 
02766 // By now all planes should be idealized, so we can idealize remaining
02767 // axes that are intersections of planes.
02768 for (i=0; i<numplanes(); i++) {
02769 int j;
02770 for (j=i+1; j<numplanes(); j++) {
02771 // Get intersection of the two planes
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 // Idealize axis by intersection of two planes
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 // Idealize axes that are bisectors of two planes.
02792 // We place this search after the one for plane intersections because 
02793 // it showed that otherwise we would get false positive bisectors.
02794 for (i=0; i<numplanes(); i++) {
02795 int j;
02796 for (j=i+1; j<numplanes(); j++) {
02797 // Get first axis bisecting the two planes
02798 float bisect1[3];
02799 vec_add(bisect1, planes[i].v, planes[j].v);
02800 vec_normalize(bisect1);
02801 
02802 // Get second axis bisecting the two planes
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 //printf("idealized bisect1 %i (C%i)\n", k, axes[k].order);
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 //printf("idealized bisect2 %i (C%i)\n", k, axes[k].order);
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 // We might still have axes that are not related to planes or we didn't
02833 // have planes at all (as in Dn). These axes should be orthogonal to
02834 // the reference.
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 // make axis truly orthogonal to reference 
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 // Idealize angle between current and first orthogonal axis
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 // Delete all axes that could not be idealized
02869 prune_axes(ideala);
02870 
02871 if (ideala) delete [] ideala;
02872 if (idealp) delete [] idealp;
02873 
02874 // Idealize rotary reflections with their corresponding rotary axes.
02875 // If no such axis is found we delete the rotary reflection.
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 // keep this rotary reflection
02898 idealr[i] = 1;
02899 }
02900 }
02901 
02902 prune_rotreflections(idealr);
02903 delete [] idealr;
02904 }
02905 
02906 }
02907 
02908 // Assuming that myaxis is close to a symmetry axis the function computes
02909 // the idealized axis, i.e. the closest symmetry axis. The reference axis
02910 // is rotated by the ideal angle around hub and the result is returned in
02911 // idealaxis. All multiples of 180/reforder are considered ideal angles.
02912 // Typically hub would correspond to an already idealized axis and one
02913 // one would just use the order of thar axis for the parameter reforder.
02914 // For instance, if the order of the principle axis is 3 than one would
02915 // expect that other elements like perpendicular axes or vertical planes
02916 // to be spaced at an angle of 180/3=60 degrees. Special values for
02917 // reforder (TETRAHEDRON, OCTAHEDRON, DODECAHEDRON) request ideal angles
02918 // (109, 90, 117) for the accordeing geometries.
02919 // If the measured angle between myaxis and refaxis is within a tolerance
02920 // of an ideal angle then the ideal axis is set and 1 is returned,
02921 // otherwise the return value is 0.
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)))); // tetrahedral angle ~109 deg
02927 const float octahedron = 90.0f; // octahedral angle = 90 deg
02928 const float dodecahedron = float(RADTODEG(acosf(-1/sqrtf(5.0f)))); // dodecahedral angle ~116.565 deg
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 // Determine the ideal angle of the axis wrt the reference axis
02960 if (fabs(alpha)<tol || fabs(alpha-180)<tol) {
02961 // same axis
02962 alpha = 0;
02963 success = 1;
02964 }
02965 
02966 if (!success) {
02967 //printf("alpha = %.2f, tol = %.2f deg, reforder=%i\n", alpha, tol, reforder);
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 // Idealize the axis by rotating the reference axis
02976 // by the ideal angle around the hub.
02977 Matrix4 rot;
02978 rot.rotate_axis(hub, float(DEGTORAD(alpha)));
02979 rot.multpoint3d(refaxis, idealaxis);
02980 
02981 return 1;
02982 }
02983 
02984 
02985 // Determine which atoms are unique for the given rotary axis.
02986 // Result is a modified array of flags 'uniquelist' where 1 indicates
02987 // that the corresponding atom is unique.
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 // Project reference coordinate on the plane defined by the rotary axis
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 // The unique atoms we have at this point will be scattered
03004 // randomly over the molecule. However, we want to find a
03005 // set of unique coordinates that is connected and confined
03006 // to one segment of the rotation.
03007 // Here we loop over all rotary images of the current unique
03008 // atom (regaring the given axis) and use the image that is
03009 // within the samme rotary segment as the reference atom
03010 // as the new unique atom.
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 // Project coordinate on the plane defined by the rotary axis
03018 float r[3];
03019 vec_scale(r, dot_prod(tmp, axis), axis);
03020 vec_sub(r, tmp, r);
03021 
03022 // Measure angle between projected coordinates of current and
03023 // reference atom. If the atom is outside the angle range we
03024 // swap with its image.
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 //printf("atom[%i] --> %i image=%i\n", i, k, image);
03033 uniqueatoms[i] = 0;
03034 uniqueatoms[k] = 1;
03035 }
03036 }
03037 
03038 // Find a connected set of unique atoms
03039 wrap_unconnected_unique_atoms(refatom, matchlist, connectedunique);
03040 }
03041 
03042 
03043 // Find best set of unique atoms, i.e. the set with the most
03044 // atoms connected to atom 'root'. 
03045 // Array matchlist shall contains the closest match of the
03046 // same atom type for each selected atom.
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 // The first time we call this function we need to construct
03054 // an initial set of unique connected atoms connected to the
03055 // first atom.
03056 if (!connectedunique) {
03057 connectedunique = new int[sel->selected];
03058 }
03059 find_connected_unique_atoms(connectedunique, root);
03060 
03061 // Repeatedly loop through the list of unconnected unique atoms
03062 // and try to swap them with images that are connected of that
03063 // are directly bonded to connected ones.
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 //printf("uniqueatoms[%d] = %d matching %d:\n", i, uniqueatoms[i], j);
03074 while (j!=i && image<sel->selected) {
03075 //printf("i=%d, j=%d\n", i, j);
03076 if (connectedunique[j]) {
03077 // If the image is already in the bondtree we can just swap
03078 swap = 1;
03079 } else {
03080 // See if one of the atoms directly bonded to the image
03081 // is in the bondtree. If yes, we can swap.
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 //printf("nonbonded unique atom %i --> %i (image %i), connected=%i\n", i, j, image, connectedunique[j]);
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 // Find all atoms currently marked unique that are within a bond
03116 // tree of unique atoms rooted at 'root'. As a result the array
03117 // 'connectedunique' is populated with flags.
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 // Loop over all endpoints of the tree
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 // for (i=0; i<sel->selected; i++) {
03158 // printf("connectedunique[%i] = %i, unique = %i\n", i, connectedunique[i], uniqueatoms[i]);
03159 // }
03160 
03161 }
03162 
03163 // Determine the unique coordinates for the whole system.
03164 // Result is a modified array of flags 'uniquelist' where 1
03165 // indicates that the corresponding atom is unique.
03166 // We begin assuming all atoms are unique and go through all
03167 // symmetry elements subsequently flagging all atoms that can
03168 // be obtained by the according symmetry operations as 
03169 // not unique.
03170 void Symmetry::unique_coordinates() {
03171 int i;
03172 for(i=0; i<sel->selected; i++) {
03173 uniqueatoms[i] = 1;
03174 }
03175 
03176 // We use the first atom as starting point for searching 
03177 // connected unique atoms. In case this atom is at the
03178 // center of mass it will always have itself as image and
03179 // we better use the second atom.
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 // -- Inversion --
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 // Flag equivalent atoms
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 // -- Rotary axes --
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 // Flag equivalent atoms
03220 int *matchlist = unique_coordinates(&rot);
03221 
03222 // Turn the unique list into a list of unique, connected
03223 // atoms that are within one segment of the rotation
03224 // (e.g within 90 degree for a 4th order axes).
03225 collapse_axis(axes[i].v, axes[i].order, refatom, matchlist,
03226 connectedunique);
03227 
03228 if (matchlist) delete [] matchlist;
03229 }
03230 
03231 // -- Planes --
03232 for (i=0; i<numplanes(); i++) {
03233 Matrix4 mirror;
03234 mirror.translate(rcom[0], rcom[1], rcom[2]);
03235 mirror.scale(-1.0); // inversion
03236 mirror.rotate_axis(planes[i].v, float(VMD_PI));
03237 mirror.translate(-rcom[0], -rcom[1], -rcom[2]);
03238 
03239 // Flag equivalent atoms
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 // -- Rotary reflections --
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); // inversion
03262 rotref.rotate_axis(rotreflections[i].v, float(VMD_PI));
03263 rotref.translate(-rcom[0], -rcom[1], -rcom[2]);
03264 
03265 // Flag equivalent atoms
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 // Goes through the list of atoms matching with original
03278 // ones after a transformation has been applied and assigns
03279 // a zero to the array of unique atom flags if the matching
03280 // index is greater than the original one.
03281 int* Symmetry::unique_coordinates(Matrix4 *trans) {
03282 int nmatches;
03283 int *matchlist = NULL;
03284 
03285 // For each atom find the closest transformed image atom with
03286 // the same chemical element. The array matchlist will contain
03287 // the closest match atom index (into the collapsed list of
03288 // selected atoms) for each selected atom.
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; // no match found;
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 // If n is the higest Cn axis order, are there n C2 axes
03337 // perpendicular to Cn?
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 // Are there n dihedral mirror planes Sd bisecting the angles
03350 // formed by pairs of C2 axes?
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 // Are there n planes vertical to the highest Cn axis?
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 { // numaxes==0
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 // Determine level in the pointgroup hierarchy
03394 // As far as I understand only crystallographic point
03395 // groups can be ranked in hierarchy, but we can use
03396 // heuristics to rank the others (e.g. I, Ih, C5, ...)
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 // Generate transformation matrix that orients the molecule
03426 // according to the GAMESS 'master frame'.
03427 //
03428 // From the GAMESS documentation:
03429 // ------------------------------
03430 // The 'master frame' is just a standard orientation for 
03431 // the molecule. By default, the 'master frame' assumes that 
03432 // 1. z is the principal rotation axis (if any), 
03433 // 2. x is a perpendicular two-fold axis (if any), 
03434 // 3. xz is the sigma-v plane (if any), and 
03435 // 4. xy is the sigma-h plane (if any). 
03436 // Use the lowest number rule that applies to your molecule.
03437 // 
03438 // Some examples of these rules: 
03439 // Ammonia (C3v): the unique H lies in the XZ plane (R1,R3). 
03440 // Ethane (D3d): the unique H lies in the YZ plane (R1,R2). 
03441 // Methane (Td): the H lies in the XYZ direction (R2). Since 
03442 // there is more than one 3-fold, R1 does not apply. 
03443 // HP=O (Cs): the mirror plane is the XY plane (R4). 
03444 // 
03445 // In general, it is a poor idea to try to reorient the 
03446 // molecule. Certain sections of the program, such as the 
03447 // orbital symmetry assignment, do not know how to deal with 
03448 // cases where the 'master frame' has been changed. 
03449 // 
03450 // Linear molecules (C4v or D4h) must lie along the z axis, 
03451 // so do not try to reorient linear molecules. 
03452 
03453 // Note:
03454 // With perpendicular two-fold axis they seem to mean any C2 axis
03455 // that is not collinear with the principal axis. It does not
03456 // have to be orthogonal to the principal axis. Otherwise the methane
03457 // example would not work.
03458 
03459 // We must use the first or the first two rules that apply. 
03460 // This is always sufficient for proper orientation.
03461 // If more that two rules apply we must ignore the additional
03462 // ones, otherwise the orientation will be messed up.
03463 
03464 void Symmetry::orient_molecule() {
03465 if (pointgroup==POINTGROUP_C1) return;
03466 
03467 //msgInfo << "Creating standard orientation:" << sendmsg;
03468 
03469 // Special case: linear molecules
03470 if (linear) {
03471 // Bring X along Z
03472 orient.transvec(0, 0, 1);
03473 // Bring axis along X
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 // GAMESS rule #1:
03484 // z is the principal rotation axis 
03485 // (x is principal axis of inertia with smallest eigenvalue)
03486 
03487 if (!sphericaltop && numaxes()) {
03488 //msgInfo << " Applying rule 1" << sendmsg;
03489 
03490 // Bring X along Z
03491 rot.transvec(0, 0, 1);
03492 // Bring first axis along X
03493 rot.transvecinv(axes[0].v[0], axes[0].v[1], axes[0].v[2]);
03494 
03495 // Find an axis of inertia that is orthogonal to the primary axis
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 // Apply same rotation to selected orthogonal axis of inerta
03503 // and then get transform m to bring it along X.
03504 Matrix4 m;
03505 float tmp[3];
03506 rot.multpoint3d(ortho_inertiaaxis, tmp);
03507 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03508 
03509 // next 2 lines: postmultiply: m*rot 
03510 m.multmatrix(rot);
03511 rot.loadmatrix(m);
03512 }
03513 
03514 // GAMESS rule #2:
03515 // x is a 'perpendicular' (=noncollinear) two-fold axis (if any)
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 //msgInfo << " Applying rule 2\n" << sendmsg;
03527 
03528 // Bring orth C2 axis along X
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 // next 2 lines: postmultiply: m*rot 
03535 m.multmatrix(rot);
03536 rot.loadmatrix(m);
03537 }
03538 
03539 // GAMESS rule #3:
03540 // xz is the sigma-v plane (if any)
03541 if (numverticalplanes && orthC2<0) {
03542 for (i=0; i<numplanes(); i++) {
03543 if (planes[i].type==VERTICALPLANE) break;
03544 }
03545 //msgInfo << "Applying rule 3\n" << sendmsg;
03546 
03547 Matrix4 m;
03548 
03549 // Bring X along Y
03550 float Y[]={0, 1, 0};
03551 m.transvec(Y[0], Y[1], Y[2]);
03552 
03553 // Bring plane normal along X 
03554 float tmp[3];
03555 rot.multpoint3d(planes[i].v, tmp);
03556 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03557 
03558 // next 2 lines: postmultiply: m*rot 
03559 m.multmatrix(rot);
03560 rot.loadmatrix(m);
03561 }
03562 
03563 // GAMESS rule #4:
03564 // xy is the sigma-h plane (if any)
03565 // If the pointgroup is Cs then treat the only plane as horizontal
03566 if ((horizontalplane>=0 && !numverticalplanes && orthC2<0) ||
03567 pointgroup==POINTGROUP_CS) {
03568 //msgInfo << "Applying rule 4\n" << sendmsg;
03569 if (pointgroup==POINTGROUP_CS) i = 0;
03570 else i = horizontalplane;
03571 
03572 Matrix4 m;
03573 float Z[]={0, 0, 1};
03574 
03575 // Bring X along Z
03576 m.transvec(Z[0], Z[1], Z[2]);
03577 
03578 // Bring plane normal along X 
03579 float tmp[3];
03580 rot.multpoint3d(planes[i].v, tmp);
03581 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03582 
03583 // next 2 lines: postmultiply: m*rot 
03584 m.multmatrix(rot);
03585 rot.loadmatrix(m);
03586 }
03587 
03588 if (pointgroup>=POINTGROUP_T && pointgroup<=POINTGROUP_IH) {
03589 //msgInfo << "Applying rule 5\n" << sendmsg;
03590 
03591 // Find a principal axis with a unique atom
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 // Bring X along XYZ
03621 m.transvec(XYZ[0], XYZ[1], XYZ[2]);
03622 
03623 // Bring unique axis along X
03624 float tmp[3];
03625 rot.multpoint3d(uniqueaxis, tmp);
03626 m.transvecinv(tmp[0], tmp[1], tmp[2]);
03627 
03628 // next 2 lines: postmultiply: m*rot 
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 // Build a summary matrix of found symmetry elements
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 // Create human-readable string summarizing symmetry elements
03761 // given in summary.
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 // For the given point group name compare the ideal numbers of 
03810 // symmetry elements with the found ones and determine which
03811 // elements are missing and which were found in addition to the
03812 // ideal ones.
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 // If the found elements differ from the ideal elements
03884 // print a comparison
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 // Draws a atom-colored spheres for each atom at the transformed
03901 // position and labels them according to the atom ID.
03902 // Use for debuggging only!
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 /*********** HELPER FUNCTIONS ************/
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 // currently unused
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 // Returns 1 if x is a prime number
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 // Returns the number of prime factors for x
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 // Checks if the molecule is planar with respect to the given plane normal.
04015 // Rotates the molecule so that he given normal is aligned with the x-axis
04016 // and then tests if all x coordinates are close to zero.
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 // Assign bond topology information to each atom
04044 void Symmetry::assign_bonds() {
04045 Molecule *mol = mlist->mol_from_id(sel->molid());
04046 
04047 // Assign an array of atom indexes associating local indexes
04048 // with the global ones.
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 // Only consider bonds to atoms within the selection
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 // Find the local index of the bonded atom
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 // If the local index is not found then this means that
04082 // the bonded atom is outside the selection.
04083 if (found) {
04084 // Add a new bond to the list
04085 Bond b;
04086 b.atom0 = i; // index into collapsed atomlist
04087 b.atom1 = m; // index into collapsed atomlist
04088 b.order = order;
04089 //printf("i=%i, m=%i, j=%i, bondto=%i\n", i, m, j, bondto);
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 // Build a list of vectors that are the sum of all bond directions
04121 // weighted by bondorder. This serves like a checksum for bonding
04122 // topology and orientation for each atom. Comparing the bondsum
04123 // between atoms and its transformed images can be used to purge
04124 // symmetry elements that would reorient the bonding pattern.
04125 // For example, two sp2 carbons will always have the same atom type
04126 // but the bondsum of will point in direction of the double bond, 
04127 // thus allowing to detect if the bonding pattern has changed.
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 // Only consider bonds to atoms within the selection
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 // Copy coordinates of selected atoms into local array and assign
04153 // and atomtypes based on chemial element and topology.
04154 // Since we also use this in measure_trans_overlap() we don't make it
04155 // a class member of Symmetry:: and pass parameters instead.
04156 static void assign_atoms(AtomSel *sel, MoleculeList *mlist, float *(&mycoor), int *(&atomtype)) {
04157 // get atom coordinates
04158 const float *allcoor = sel->coordinates(mlist);
04159 
04160 Molecule *mol = mlist->mol_from_id(sel->molid());
04161 
04162 // array of strings describing the atom type
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 // copy coordinates of selected atoms into local array
04171 vec_copy(mycoor+3L*j, allcoor+3L*i);
04172 
04173 
04174 // Calculate lightest and heaviest element bonded to this atom
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 // Build up a string describing the atom type. It is not meant
04186 // to be human-readable, so it's not nice but since the contained
04187 // properties are ordered it can be used to compare two atoms.
04188 char *typestring = new char[8L+12L*mol->atom(i)->bonds];
04189 typestring[0] = '0円';
04190 char buf[100] = { 0 };
04191 
04192 // atomic number and number of bonds
04193 sprintf(buf, "%i %i ", mol->atom(i)->atomicnumber, mol->atom(i)->bonds);
04194 strcat(typestring, buf);
04195 
04196 // For each chemical element get the number of bonds for each 
04197 // bond order. We distinguish half and integer bond orders,
04198 // e.g. 2 and 3 mean bond orders 1 and 1.5 respectively.
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 // Try to find this atom's type in the list, if it doesn't exist
04222 // add a new string and numerical type. 
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 //printf("%i: type=%i {%s}\n", j, atomtype[j], typestringptr[atomtype[j]]);
04238 j++;
04239 }
04240 
04241 for (i=0; i<numtypes; i++) {
04242 delete [] typestringptr[i];
04243 }
04244 delete [] typestringptr;
04245 }
04246 
04247 
04248 // Calculate the structural overlap of a set of points with a copy of
04249 // themselves that was transformed according to the given transformation
04250 // matrix.
04251 // Returns the normalized sum over all gaussian function values of the
04252 // pair distances between the atoms in the original and the transformed
04253 // position.
04254 // Two atoms are only considered overlapping if their atom type is
04255 // identical. Atom types must be provided as an array with an integer
04256 // for each atom.
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 // get atom coordinates
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 // generate transformed coordinates
04276 int i, ncompare=0;
04277 for(i=0; i<numcoor; i++) {
04278 trans->multpoint3d(posA+3L*i, posB+3L*i);
04279 
04280 // Depending on the flag skip atoms that underwent an almost
04281 // identical transformation.
04282 if (!(skipident && distance(posA+3L*i, posB+3L*i) < sigma)) {
04283 flgs[i] = 1;
04284 ncompare++; // # of compared atoms with dist<sigma
04285 }
04286 }
04287 
04288 if (ncompare<0.5*numcoor) {
04289 // Not enough atoms to compare
04290 delete [] flgs;
04291 delete [] posB;
04292 return 0.0;
04293 }
04294 
04295 // If the pair distance gets too small the performance is terrible
04296 // even for small systems. So we limit it to at least 4.5A
04297 float dist;
04298 //if (sigma<1.5) dist = 4.5;
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 // Now go through all atoms and find matching pairs
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 // Find the nearest atom
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 // Compute the score for the closest atom
04326 if (minr2<maxr2) {
04327 noverlap++;
04328 
04329 // consider only pairs with identical atom types
04330 if (atomtype[i1]==atomtype[i2]) {
04331 // Gaussian function of the pair distance
04332 overlap += expf(-itwosig2*minr2);
04333 nmatches++;
04334 }
04335 else {
04336 // wrong element matching 
04337 antioverlap += wrongelementpenalty*expf(-itwosig2*r2);
04338 nwrongelement++;
04339 //printf("wrong elements %i-%i (atoms %i-%i)\n", mol->atom(i1)->atomicnumber,mol->atom(i2)->atomicnumber,i1,i2);
04340 }
04341 }
04342 }
04343 
04344 float nomatchpenalty = 0.0;
04345 overlappermatch = 0.0;
04346 int numnomatch = ncompare-nmatches;
04347 
04348 // We make the penalty for unmatched atoms dependent on the
04349 // average overlap for the overlapping atoms with correct element
04350 // matching. In other words, the better the structures match in
04351 // general the higher the penalty for a nonmatching element.
04352 // This helps for instance to prevent mistaking nitrogens for
04353 // carbons in rings.
04354 // On the other hand the algorithm is, at least for larger structures,
04355 // fairly forgiving about not matching an atom at all. This helps
04356 // recognizing symmetry when one or very few atoms are astray.
04357 if (nmatches) overlappermatch = overlap/nmatches;
04358 
04359 overlap -= antioverlap;
04360 
04361 if (!(numnomatch==0)) {
04362 //nomatchpenalty = powf(overlappermatch, 5)*(1-powf(0.2f+float(numnomatch),-2));
04363 nomatchpenalty = powf(overlappermatch, 5);//*(1-powf(0.1f,4.0f)/powf(float(numnomatch)/ncompare,4.0f));
04364 //overlap *= 1- nomatchpenalty;
04365 overlap -= 8*numnomatch*nomatchpenalty;
04366 }
04367 if (overlap<0) overlap = 0.0f;
04368 
04369 overlap /= ncompare;
04370 
04371 //printf("nsel=%i, wrongelementpen=%.2f, nmatches/noverlap/ncompare=%i/%i/%i, overlap/match=%.2f, nomatchpen=%.2f, ov=%.2f\n",
04372 //numcoor, wrongelementpenalty, nmatches, noverlap, ncompare, overlappermatch, nomatchpenalty, overlap);
04373 
04374 delete [] flgs;
04375 delete [] posB;
04376 
04377 return overlap;
04378 }
04379 
04380 
04381 
04382 /******* OTHER FUNCTIONS WITH TCL INTERFACE *********/
04383 
04384 
04385 
04386 // Backend of the TCL interface for measure transoverlap:
04387 // Calculate the structural overlap of a selection with a copy of itself
04388 // that is transformed according to a given transformation matrix.
04389 // Returns the normalized sum over all gaussian function values of the
04390 // pair distances between atoms in the original and the transformed
04391 // selection.
04392 // Two atoms are only considered overlapping if their atom type is
04393 // identical. The atom type is determined based on the chemical element
04394 // and number and order of bonds to different elements.
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 // Calculate the structural overlap between two point sets.
04413 // The overlap of two atoms is defined as the Gaussian of their distance.
04414 // If two atoms have identical positions the overlap is 1. The parameter
04415 // sigma controls the greediness (i.e. the width) of the overlap function.
04416 // Returns the average overlap for all atoms.
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