Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

MeasureSymmetry.C

Go to the documentation of this file.
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 

Generated on Fri Nov 8 02:44:50 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002