---

MeasureSymmetry.h

Go to the documentation of this file.

/***************************************************************************
 *cr                                                                       
 *cr            (C) Copyright 1995-2019 The Board of Trustees of the           
 *cr                        University of Illinois                       
 *cr                         All Rights Reserved                        
 *cr                                                                   
 ***************************************************************************/

/***************************************************************************
 * RCS INFORMATION:
 *
 *      $RCSfile: MeasureSymmetry.h,v $
 *      $Author: johns $        $Locker:  $             $State: Exp $
 *      $Revision: 1.34 $       $Date: 2020/10/28 17:42:35 $
 *
 ***************************************************************************
 * DESCRIPTION:
 * Thy Symmetry class is the work horse behind the "measure symmetry"
 * command which guesses the pointgroup of a selection and returns the
 * according symmetry elements (mirror planes rotary axes, rotary reflection
 * axes). The algorithm is is fairly forgiving about molecules where atoms
 * are perturbed from the ideal position and tries its best to guess the
 * correct point group anyway. 
 * The 'forgivingness' can be controlled with the sigma parameter which is
 * the average allowed deviation from the ideal position.
 * Works nice on my 30 non-patholocical test cases. A pathological case
 * would for instance be a system with more than only a few atoms (say 15)
 * where only one atom distinguishes between two point groups.
 * If your selection contains more than a certain number of atoms then only
 * that much randomly chosen atoms are used to find planes and axes in order
 * to save time.
 *
 ***************************************************************************/

#ifndef MEASURESYMMETRY_H
#define MEASURESYMMETRY_H

#include "ResizeArray.h"

// Flags for trans_overlap() specifying wether atoms that underwent
// an identical transformation should be skipped in the evaluation.
#define SKIP_IDENTICAL   1
#define NOSKIP_IDENTICAL 0

// Rotary axis classification
#define PRINCIPAL_AXIS     1
#define PERPENDICULAR_AXIS 2

// Maximum order for Cn and S2n axes
#define MAXORDERCN 8


typedef struct plane {
  float v[3];
  float overlap;
  int weight;
  int type;
} Plane;

typedef struct axis {
  float v[3];
  float overlap;
  int order;
  int weight;
  int type;
} Axis;

// Summary of existing symmetry elements
typedef struct elementsummary {
  char inv;     // inversion center (0/1)?
  char sigma;   // number of sigma planes
  char Cinf;    // number of infinite order axes (0/1)
  // Numbers of C and S axes with order 2, 3, 4, ... repectively
  char C[MAXORDERCN+1]; // 0,1 unused; accessed by order 2, ... MAXORDERCN;
  char S[MAXORDERCN*2+1];
} ElementSummary;

typedef struct bond {
  int atom0;
  int atom1;
  float order;
  float length;
} Bond;

typedef struct bondlist {
  int numbonds;
  int *bondto;
  float *length;
} Bondlist;


class Symmetry {
private:
  AtomSel *sel;              
  MoleculeList *mlist;

  float *coor;               
  float *bondsum;            
  float rcom[3];             
  float inertiaeigenval[3];  
  float inertiaaxes[3][3];   
  int uniqueprimary[3];      
  int numverticalplanes;     
  int numdihedralplanes;     
  int horizontalplane;       
  int maxaxisorder;          
  int maxrotreflorder;       
  int maxnatoms;             
  int checkbonds;            
  int verbose;               
  int *atomtype;             
  int *atomindex;            
  int *uniqueatoms;          
  float *idealcoor;          

  char *elementsummarystring;    
  char *missingelementstring;    
  char *additionalelementstring; 

  Matrix4 orient;            

  // Summarized numbers of elements in table form.
  ElementSummary elementsummary;

  int   maxweight;           
                             // XXX seems to be set but never used;
  float maxoverlap;          
  float rmsd;                

  int pointgroup;            

 
  int pointgrouporder;       

 
  ResizeArray<Plane> planes;         
  ResizeArray<Axis>  axes;           
  ResizeArray<Axis>  rotreflections; 
  ResizeArray<Bond>  bonds;          

  Bondlist *bondsperatom;  


  bool linear;        
  bool planar;        
  bool inversion;     
  bool sphericaltop;  
  bool symmetrictop;  


  int iterate_element_search();

  void find_symmetry_elements();

  void find_elements_from_inertia();

  void find_planes();

  void find_axes_from_plane_intersections();

  void find_C2axes();
  void find_axes(int order);


  void check_add_inversion();   

  void check_add_plane(const float *normal, int weight=1);

  void check_add_axis(const float *testaxis, int order);

  void check_add_rotary_reflection(const float *testaxis, int maxorder);


  void classify_planes();

  void classify_axes();

  void purge_planes(float cutoff);

  void purge_axes(float cutoff);

  void purge_rotreflections(float cutoff);

  void prune_planes(int *keep);

  void prune_axes(int *keep);

  void prune_rotreflections(int *keep);

  void assign_axis_orders();

  void assign_prelimaxis_orders(int order);

  int axis_order(const float *axis, float *overlap);

  void sort_axes();

  void sort_planes();

  void print_statistics();

  void build_element_summary();

  void build_element_summary_string(ElementSummary summary, char *(&sestring));

  void compare_element_summary(const char *pointgroupname);

  void print_element_summary(const char *pointgroupname);

  void idealize_elements();

  int idealize_angle(const float *refaxis, const float *hub,
                    const float *myaxis, float *idealaxis, int reforder);

  void idealize_coordinates();

  int average_coordinates(Matrix4 *trans);

  // transformation trans.
  int check_bondsum(int j, int m, Matrix4 *trans);

  void identify_transformed_atoms(Matrix4 *trans, int &nmatches, int *(&matchlist));

  float ideal_coordinate_rmsd ();

  int ideal_coordinate_sanity();

  void unique_coordinates();

  int* unique_coordinates(Matrix4 *trans);

  void collapse_axis(const float *axis, int order, int refatom, const int *matchlist, int *(&connectedunique));

  int pointgroup_rank(int pg, int order);

  void orient_molecule();

  void determine_pointgroup();

  inline int find_collinear_axis(const float *myaxis);

  inline int plane_exists(const float *myplane);

  int is_planar(const float *axis);

  void assign_bonds();

  void assign_bondvectors();

  void wrap_unconnected_unique_atoms(int root, const int *matchlist, int *(&connectedunique));

  // Find all unique atoms that are connnected to root
  void find_connected_unique_atoms(int *(&connectedunique), int root);

  void draw_transformed_mol(Matrix4 rot);

public:
  Symmetry(AtomSel *sel, MoleculeList *mlist, int verbosity);
  ~Symmetry(void);

  float sigma;     // atomic overlap tolerance
  float collintol; // collinearity tolerance in radians
  float orthogtol; // coplanarity tolerance in radians

  // Check if center of mass is an inversion center.
  float score_inversion();

  float score_axis(const float *testaxis, int order);

  float score_plane(const float *normal);

  float score_rotary_reflection(const float *testaxis, int order);

  void impose(int have_inversion, 
              int nplanes, const float *planev,
              int naxes, const float *axisv, const int *axisorder,
              int nrotrefl, const float* rotreflv,
              const int *rotreflorder);

  int guess(float mysigma);

  void get_pointgroup(char pg[8], int *order);

  float get_rmsd() { return rmsd; }

  float* get_inertia_axes()      { return &(inertiaaxes[0][0]); }

  float* get_inertia_eigenvals() { return &(inertiaeigenval[0]); }

  int*   get_inertia_unique()    { return &(uniqueprimary[0]); }

  int  get_axisorder(int n);

  int  get_rotreflectorder(int n);

  Matrix4* get_orientation() { return &orient; }

  const char *get_element_summary()     { return elementsummarystring; }

  const char *get_missing_elements()    { return missingelementstring; }

  const char *get_additional_elements() { return additionalelementstring; }

  void set_maxnumatoms(int n) { maxnatoms = n; }

  void set_overlaptol(float tol) { sigma = tol; }

  void set_checkbonds(int flag) { checkbonds = flag; }

  int numplanes()     { return int(planes.num()); }

  int numaxes()       { return int(axes.num()); }

  int numprimaryaxes();

  int numrotreflect() { return int(rotreflections.num()); }

  int numS2n();

  int has_inversion() { return inversion; }

  int is_spherical_top() { return sphericaltop; }

  int is_symmetric_top() { return symmetrictop; }

  float* center_of_mass() { return rcom; }

  float* plane(int i)      { return (float *) &(planes[i].v); }

  float* axis(int i)       { return (float *) &(axes[i].v); }

  float* rotreflect(int i) { return (float *) &(rotreflections[i].v); }

  int get_planetype(int i)     { return planes[i].type; }

  int get_axistype(int i)      { return axes[i].type; }

  int get_rotrefltype(int i)   { return rotreflections[i].type; }

  float *idealpos(int i)   { return idealcoor+3*i; }

  int get_unique_atom(int i)   { return uniqueatoms[i]; }
};



int measure_trans_overlap(AtomSel *sel, MoleculeList *mlist, const Matrix4 *trans,
                          float sigma, bool skipident, int maxnatoms, float &overlap);


int measure_pointset_overlap(const float *posA, int natomsA, int *flagsA,
                             const float *posB, int natomsB, int *flagsB,
                             float sigma, float pairdist, float &similarity);




#endif

---