Subversion Repositories wimsdev

Rev

Rev 3662 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. /*
  2.     Sketch Elements: Chemistry molecular diagram drawing tool.
  3.    
  4.     (c) 2005 Dr. Alex M. Clark
  5.    
  6.     Released as GNUware, under the Gnu Public License (GPL)
  7.    
  8.     See www.gnu.org for details.
  9. */
  10.  
  11. package WIMSchem;
  12.  
  13. import java.util.*;
  14.  
  15. /*
  16.     Internal molecule datastructure, which captures the basic information about a molecule as envisaged by the traditional diagram
  17.     notation. The core data is stored as a list of labelled nodes (atoms) and labelled edges (bonds). Besides very basic editing and
  18.     convenient object manipulation, a few computed properties are available, some of which are cached.
  19. */
  20.  
  21. public class Molecule
  22. {
  23.     class Atom
  24.     {
  25.         String Element;
  26.         double X,Y;
  27.         int Charge,Unpaired;
  28.         int HExplicit;
  29.         int MapNum;
  30.     }
  31.     class Bond
  32.     {
  33.         int From,To;
  34.         int Order,Type;
  35.     }
  36.    
  37.     public static final int HEXPLICIT_UNKNOWN=-1;
  38.    
  39.     public static final int BONDTYPE_NORMAL=0;
  40.     public static final int BONDTYPE_INCLINED=1;
  41.     public static final int BONDTYPE_DECLINED=2;
  42.     public static final int BONDTYPE_UNKNOWN=3;
  43.    
  44.     public static final String[] ELEMENTS=
  45.         {null,"H","He","Li","Be","B","C","N","O","F","Ne","Na","Mg","Al","Si","P","S","Cl","Ar","K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co",
  46.         "Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I",
  47.         "Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W","Re","Os","Ir","Pt","Au",
  48.         "Hg","Tl","Pb","Bi","Po","At","Rn","Fr","Ra","Ac","Th","Pa","U","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr"};
  49.  
  50.     private boolean invalMinMax=false;
  51.     private double minX=0,minY=0,maxX=0,maxY=0;
  52.     private ArrayList<Atom> atoms=new ArrayList<Atom>();
  53.     private ArrayList<Bond> bonds=new ArrayList<Bond>();
  54.    
  55.     // computed graph topology properties
  56.     private int graph[][]=null;
  57.     private int ringID[]=null,compID[]=null,priority[]=null;
  58.     private ArrayList<Object> ring3=null,ring4=null,ring5=null,ring6=null,ring7=null; // common ring sizes
  59.    
  60.     // computed stereochemistry properties
  61.     public static final int STEREO_NONE=0; // topology does not allow any stereoisomers
  62.     public static final int STEREO_POS=1; // R or Z, depending on type
  63.     public static final int STEREO_NEG=2; // S or E, depending on type
  64.     public static final int STEREO_UNKNOWN=3; // topology allows stereoisomers, but wedges/geometry do not resolve which
  65.     private int chiral[]=null; // (per atom)
  66.     private int cistrans[]=null; // (per bond)
  67.  
  68.     // data for calculation of implicit hydrogens
  69.     static final String[] HYVALENCE_EL={"C","N","O","S","P"};
  70.     static final int[]   HYVALENCE_VAL={ 4,  3,  2,  2,  3 };
  71.    
  72.     // ------------------ public functions --------------------
  73.  
  74.     public Molecule()
  75.     {
  76.     }
  77.  
  78.     // returns an equivalent replica of the molecule
  79.     public Molecule clone()
  80.     {
  81.         Molecule mol=new Molecule();
  82.         for (int n=1;n<=numAtoms();n++)
  83.         {
  84.             Atom a=atoms.get(n-1);
  85.             int num=mol.addAtom(a.Element,a.X,a.Y,a.Charge,a.Unpaired);
  86.             mol.setAtomHExplicit(num,atomHExplicit(n));
  87.             mol.setAtomMapNum(num,atomMapNum(n));
  88.         }
  89.         for (int n=1;n<=numBonds();n++)
  90.         {
  91.             Bond b=bonds.get(n-1);
  92.             mol.addBond(b.From,b.To,b.Order,b.Type);
  93.         }
  94.         return mol;
  95.     }
  96.    
  97.     // ----- fetching information directly from the datastructures ----
  98.        
  99.     public int numAtoms() {return atoms.size();}
  100.     public int numBonds() {return bonds.size();}
  101.    
  102.     public String atomElement(int N) {return atoms.get(N-1).Element;}
  103.     public double atomX(int N) {return atoms.get(N-1).X;}
  104.     public double atomY(int N) {return atoms.get(N-1).Y;}
  105.     public int atomCharge(int N) {return atoms.get(N-1).Charge;}
  106.     public int atomUnpaired(int N) {return atoms.get(N-1).Unpaired;}
  107.     public int atomHExplicit(int N) {return atoms.get(N-1).HExplicit;}
  108.     public int atomMapNum(int N) {return atoms.get(N-1).MapNum;}
  109.    
  110.     public int bondFrom(int N) {return bonds.get(N-1).From;}
  111.     public int bondTo(int N) {return bonds.get(N-1).To;}
  112.     public int bondOrder(int N) {return bonds.get(N-1).Order;}
  113.     public int bondType(int N) {return bonds.get(N-1).Type;}
  114.    
  115.     // ----- modifying information in the datastructures ----
  116.  
  117.     // adds an atom to the end of the molecule; it is fairly reasonable to assume that the new one will be at position==numAtoms(),
  118.     // and that all other atoms will be in the same places, but the position is returned for convenience
  119.     public int addAtom(String Element,double X,double Y) {return addAtom(Element,X,Y,0,0);}
  120.     public int addAtom(String Element,double X,double Y,int Charge,int Unpaired)
  121.     {
  122.         if (X<minX || numAtoms()==0) minX=X;
  123.         if (X>maxX || numAtoms()==0) maxX=X;
  124.         if (Y<minY || numAtoms()==0) minY=Y;
  125.         if (Y>maxY || numAtoms()==0) maxY=Y;
  126.        
  127.         Atom a=new Atom();
  128.         a.Element=Element;
  129.         a.X=X;
  130.         a.Y=Y;
  131.         a.Charge=Charge;
  132.         a.Unpaired=Unpaired;
  133.         a.HExplicit=HEXPLICIT_UNKNOWN;
  134.         atoms.add(a);
  135.         trashGraph();
  136.         return atoms.size();
  137.     }
  138.     public void setAtomElement(int N,String Element) {atoms.get(N-1).Element=Element; trashPriority();}
  139.     public void setAtomPos(int N,double X,double Y)
  140.     {
  141.         atoms.get(N-1).X=X;
  142.         atoms.get(N-1).Y=Y;
  143.         invalMinMax=true;
  144.         trashStereo();
  145.     }
  146.     public void setAtomCharge(int N,int Charge) {atoms.get(N-1).Charge=Charge; trashPriority();}
  147.     public void setAtomUnpaired(int N,int Unpaired) {atoms.get(N-1).Unpaired=Unpaired; trashPriority();}
  148.     public void setAtomHExplicit(int N,int HExplicit) {atoms.get(N-1).HExplicit=HExplicit; trashPriority();}
  149.     public void setAtomMapNum(int N,int MapNum) {atoms.get(N-1).MapNum=MapNum;}
  150.    
  151.     // adds an atom the end of the molecule; the position of the new bond is returned for convenience
  152.     public int addBond(int From,int To,int Order) {return addBond(From,To,Order,BONDTYPE_NORMAL);}
  153.     public int addBond(int From,int To,int Order,int Type)
  154.     {
  155.         Bond b=new Bond();
  156.         b.From=From;
  157.         b.To=To;
  158.         b.Order=Order;
  159.         b.Type=Type;
  160.         bonds.add(b);
  161.         trashGraph();
  162.         return bonds.size();
  163.     }
  164.    
  165.     public void setBondFromTo(int N,int From,int To)
  166.     {
  167.         bonds.get(N-1).From=From;
  168.         bonds.get(N-1).To=To;
  169.         trashGraph();
  170.     }
  171.     public void setBondOrder(int N,int Order) {bonds.get(N-1).Order=Order; trashPriority();}
  172.     public void setBondType(int N,int Type) {bonds.get(N-1).Type=Type; trashStereo();}
  173.    
  174.     // deletes the indicated atom, and also any bonds which were connected to it
  175.     public void deleteAtomAndBonds(int N)
  176.     {
  177.         for (int n=0;n<numBonds();)
  178.         {
  179.             Bond b=bonds.get(n);
  180.             if (b.From==N || b.To==N) bonds.remove(n);
  181.             else
  182.             {
  183.                 if (b.From>N) b.From--;
  184.                 if (b.To>N) b.To--;
  185.                 n++;
  186.             }
  187.         }
  188.         atoms.remove(N-1);
  189.         invalMinMax=true;
  190.         trashGraph();
  191.     }
  192.     public void deleteBond(int N)
  193.     {
  194.         bonds.remove(N-1);
  195.         trashGraph();
  196.     }
  197.      
  198.     // ----- fetching information which is computed from the underlying data, usually cached -----
  199.      
  200.     public double minX() {if (invalMinMax) determineMinMax(); return minX;}
  201.     public double maxX() {if (invalMinMax) determineMinMax(); return maxX;}
  202.     public double minY() {if (invalMinMax) determineMinMax(); return minY;}
  203.     public double maxY() {if (invalMinMax) determineMinMax(); return maxY;}
  204.     public double rangeX() {return maxX-minX;}
  205.     public double rangeY() {return maxY-minY;}
  206.    
  207.     // return the index of the bond, if any, which connects the two indicated atoms
  208.     public int findBond(int A1,int A2)
  209.     {
  210.         for (int n=1;n<=numBonds();n++)
  211.         {
  212.             if (bondFrom(n)==A1 && bondTo(n)==A2) return n;
  213.             if (bondTo(n)==A1 && bondFrom(n)==A2) return n;
  214.         }
  215.         return 0;
  216.     }
  217.  
  218.     // for a bond, returns the end which is not==Ref; return value will be From,To or 0    
  219.     public int bondOther(int N,int Ref)
  220.     {
  221.         if (bondFrom(N)==Ref) return bondTo(N);
  222.         if (bondTo(N)==Ref) return bondFrom(N);
  223.         return 0;
  224.     }
  225.  
  226.     // returns whether an atom logically ought to be drawn "explicitly";if false, it is a carbon which should be just a dot
  227.     public boolean atomExplicit(int N)
  228.     {
  229.         if (atoms.get(N-1).Element.compareTo("C")!=0 || atoms.get(N-1).Charge!=0 || atoms.get(N-1).Unpaired!=0) return true;
  230.         for (int n=0;n<bonds.size();n++) if (bonds.get(n).From==N || bonds.get(n).To==N) return false;
  231.         return true;
  232.     }
  233.    
  234.     // uses either explicit or computed number to determine how many hydrogens the atom has; the field for explicit hydrogens takes
  235.     // absolute preference, if it has its default value of 'unknown', the number is computed by looking up the hydrogen capacity for
  236.     // the element (most of which are zero), subtracting the total of bond orders, then returning the difference, or zero; the calculation
  237.     // tends to err on the side of too few, since the concept is just an aid to drawing organic structures, not a cheminformatic attempt
  238.     // to compensate for 2 1/2 decades of bad file formats
  239.     // (note: returns "implicit"+"explicit", but does NOT count "actual" hydrogens, i.e. those which have their own atom nodes)
  240.     public int atomHydrogens(int N)
  241.     {
  242.         int hy=atomHExplicit(N);
  243.         if (hy!=HEXPLICIT_UNKNOWN) return hy;
  244.        
  245.         for (int n=0;n<HYVALENCE_EL.length;n++) if (HYVALENCE_EL[n].compareTo(atomElement(N))==0) {hy=HYVALENCE_VAL[n]; break;}
  246.         if (hy==HEXPLICIT_UNKNOWN) return 0;
  247.         int ch=atomCharge(N);
  248.         if (atomElement(N).compareTo("C")==0) ch=-Math.abs(ch);
  249.         hy+=ch-atomUnpaired(N);
  250.         for (int n=1;n<=numBonds();n++) if (bondFrom(n)==N || bondTo(n)==N) hy-=bondOrder(n);
  251.         return hy<0 ? 0 : hy;
  252.     }
  253.      
  254.     // returns the numerical ID of the ring block in which the atom resides, or 0 if it is not in a ring  
  255.     public int atomRingBlock(int N)
  256.     {
  257.         if (graph==null) buildGraph();
  258.         if (ringID==null) buildRingID();
  259.         return ringID[N-1];
  260.     }
  261.     // returns whether or not a bond resides in a ring (i.e. ring block of each end is the same and nonzero)
  262.     public boolean bondInRing(int N)
  263.     {
  264.         int r1=atomRingBlock(bondFrom(N)),r2=atomRingBlock(bondTo(N));
  265.         return r1>0 && r1==r2;
  266.     }
  267.    
  268.     // returns the connected component ID of the atom, always 1 or more
  269.     public int atomConnComp(int N)
  270.     {
  271.         if (graph==null) buildGraph();
  272.         if (compID==null) buildConnComp();
  273.         return compID[N-1];
  274.     }
  275.    
  276.     // returns the number of neighbours of an atom
  277.     public int atomAdjCount(int N)
  278.     {
  279.         if (graph==null) buildGraph();
  280.         return graph[N-1].length;
  281.     }
  282.    
  283.     // returns the actual adjacency list of an atom; unfortunately this has to make a clone of the array, since the numbers are converted
  284.     // to 1-based indices, and we would not want callers modifying it anyway
  285.     public int[] atomAdjList(int N)
  286.     {
  287.         if (graph==null) buildGraph();
  288.         int[] adj=graph[N-1].clone();
  289.         for (int n=0;n<adj.length;n++) adj[n]++;
  290.         return adj;
  291.     }
  292.    
  293.     // returns atom-based priority according to the Cahn-Ingold-Prelog rules
  294.     public int atomPriority(int N)
  295.     {
  296.         if (graph==null) buildGraph();
  297.         if (compID==null) buildConnComp();
  298.         if (priority==null) buildPriority();
  299.         return priority[N-1];
  300.     }
  301.    
  302.     // return the tetrahedral chirality of the given atom
  303.     public int atomChirality(int N)
  304.     {
  305.         if (graph==null) buildGraph();
  306.         if (compID==null) buildConnComp();
  307.         if (priority==null) buildPriority();
  308.         if (chiral==null) buildChirality();
  309.         return chiral[N-1];
  310.     }
  311.    
  312.     // return the cis/trans style stereochemistry of the given bond
  313.     public int bondStereo(int N)
  314.     {
  315.         if (graph==null) buildGraph();
  316.         if (compID==null) buildConnComp();
  317.         if (priority==null) buildPriority();
  318.         if (cistrans==null) buildCisTrans();
  319.         return cistrans[N-1];
  320.     }
  321.    
  322.     // returns _all_ rings of indicated size; each item in the array list is an array of int[Size], a consecutively ordered array of atom
  323.     // numbers; uses a recursive depth first search, which must be bounded above by Size being small in order to avoid exponential blowup
  324.     public int[][] findRingSize(int Size)
  325.     {
  326.         ArrayList<Object> rings=null;
  327.         if (Size==3 && ring3!=null) rings=ring3;
  328.         if (Size==4 && ring4!=null) rings=ring4;
  329.         if (Size==5 && ring5!=null) rings=ring5;
  330.         if (Size==6 && ring6!=null) rings=ring6;
  331.         if (Size==7 && ring7!=null) rings=ring7;
  332.    
  333.         if (rings==null)
  334.         {
  335.             if (graph==null) buildGraph();
  336.             if (ringID==null) buildRingID();
  337.  
  338.             rings=new ArrayList<Object>();
  339.             for (int n=1;n<=numAtoms();n++) if (ringID[n-1]>0)
  340.             {
  341.                 int path[]=new int[Size];
  342.                 path[0]=n;
  343.                 recursiveRingFind(path,1,Size,ringID[n-1],rings);
  344.             }
  345.  
  346.             if (Size==3) ring3=rings;
  347.             if (Size==4) ring4=rings;
  348.             if (Size==5) ring5=rings;
  349.             if (Size==6) ring6=rings;
  350.             if (Size==7) ring7=rings;
  351.         }
  352.        
  353.         int ret[][]=new int[rings.size()][];
  354.         for (int n=0;n<rings.size();n++) ret[n]=((int[])rings.get(n)).clone();
  355.         return ret;
  356.     }
  357.    
  358.     // compare to another molecule instance; null is equivalent to empty; return values:
  359.     //     -1 if Mol < Other
  360.     //      0 if Mol == Other
  361.     //      1 if Mol > Other
  362.     //  can be used to sort molecules, albeit crudely; does not do any kind of canonical preprocessing
  363.     public int compareTo(Molecule Other)
  364.     {
  365.         if (Other==null) return numAtoms()==0 ? 0 : 1; // null is equivalent to empty
  366.         if (numAtoms()<Other.numAtoms()) return -1;
  367.         if (numAtoms()>Other.numAtoms()) return 1;
  368.         if (numBonds()<Other.numBonds()) return -1;
  369.         if (numBonds()>Other.numBonds()) return 1;
  370.        
  371.         for (int n=1;n<=numAtoms();n++)
  372.         {
  373.             int c=atomElement(n).compareTo(Other.atomElement(n));
  374.             if (c!=0) return c;
  375.             if (atomX(n)<Other.atomX(n)) return -1; if (atomX(n)>Other.atomX(n)) return 1;
  376.             if (atomY(n)<Other.atomY(n)) return -1; if (atomY(n)>Other.atomY(n)) return 1;
  377.             if (atomCharge(n)<Other.atomCharge(n)) return -1; if (atomCharge(n)>Other.atomCharge(n)) return 1;
  378.             if (atomUnpaired(n)<Other.atomUnpaired(n)) return -1; if (atomUnpaired(n)>Other.atomUnpaired(n)) return 1;
  379.             if (atomHExplicit(n)<Other.atomHExplicit(n)) return -1; if (atomHExplicit(n)>Other.atomHExplicit(n)) return 1;
  380.         }
  381.         for (int n=1;n<=numBonds();n++)
  382.         {
  383.             if (bondFrom(n)<Other.bondFrom(n)) return -1; if (bondFrom(n)>Other.bondFrom(n)) return 1;
  384.             if (bondTo(n)<Other.bondTo(n)) return -1; if (bondTo(n)>Other.bondTo(n)) return 1;
  385.             if (bondOrder(n)<Other.bondOrder(n)) return -1; if (bondOrder(n)>Other.bondOrder(n)) return 1;
  386.             if (bondType(n)<Other.bondType(n)) return -1; if (bondType(n)>Other.bondType(n)) return 1;
  387.         }
  388.        
  389.         return 0;
  390.     }
  391.    
  392.     // lookup the atomic number for the element, or return 0 if not in the periodic table
  393.     int atomicNumber(int N) {return atomicNumber(atomElement(N));}
  394.     int atomicNumber(String El)
  395.     {
  396.         for (int n=1;n<ELEMENTS.length;n++) if (ELEMENTS[n].compareTo(El)==0) return n;
  397.         return 0;
  398.     }
  399.  
  400.     // takes the indicated molecular fragment and appends it to the end of the current molecule, with order preserved
  401.     public void append(Molecule Frag)
  402.     {
  403.         int base=numAtoms();
  404.         for (int n=1;n<=Frag.numAtoms();n++)
  405.         {
  406.             int num=addAtom(Frag.atomElement(n),Frag.atomX(n),Frag.atomY(n),Frag.atomCharge(n),Frag.atomUnpaired(n));
  407.             setAtomHExplicit(num,Frag.atomHExplicit(n));
  408.             setAtomMapNum(num,Frag.atomMapNum(n));
  409.         }
  410.         for (int n=1;n<=Frag.numBonds();n++)
  411.         {
  412.             addBond(Frag.bondFrom(n)+base,Frag.bondTo(n)+base,Frag.bondOrder(n),Frag.bondType(n));
  413.         }
  414.     }
  415.  
  416.     // returns a molecule which is smaller than the current one, according to the mask (which must be of size #atoms); boundary cases
  417.     // are a null molecule or cloned copy
  418.     public Molecule subgraph(boolean[] Mask)
  419.     {
  420.         int[] invidx=new int[numAtoms()];
  421.         int sum=0;
  422.         for (int n=0;n<numAtoms();n++) if (Mask[n]) invidx[n]=++sum; else invidx[n]=0;
  423.         if (sum==0) return new Molecule();
  424.         if (sum==numAtoms()) return clone();
  425.        
  426.         Molecule frag=new Molecule();
  427.         for (int n=1;n<=numAtoms();n++) if (Mask[n-1])
  428.         {
  429.             int num=frag.addAtom(atomElement(n),atomX(n),atomY(n),atomCharge(n),atomUnpaired(n));
  430.             frag.setAtomHExplicit(num,atomHExplicit(n));
  431.             frag.setAtomMapNum(num,atomMapNum(n));
  432.         }
  433.         for (int n=1;n<=numBonds();n++)
  434.         {
  435.             int from=invidx[bondFrom(n)-1],to=invidx[bondTo(n)-1];
  436.             if (from>0 && to>0) frag.addBond(from,to,bondOrder(n),bondType(n));
  437.         }
  438.        
  439.         return frag;
  440.     }
  441.  
  442.    
  443.     // convenient debugging mechanism
  444.     public void dump()
  445.     {
  446.         System.out.println("#Atoms="+numAtoms()+" #Bonds="+numBonds());
  447.         for (int n=0;n<numAtoms();n++)
  448.         {
  449.             Atom a=atoms.get(n);
  450.             System.out.println(" A"+(n+1)+": "+a.Element+"="+a.X+","+a.Y+";"+a.Charge+","+a.Unpaired);
  451.         }
  452.         for (int n=0;n<numBonds();n++)
  453.         {
  454.             Bond b=bonds.get(n);
  455.             System.out.println(" B"+(n+1)+": "+b.From+"-"+b.To+"="+b.Order+","+b.Type);
  456.         }
  457.     }
  458.  
  459.     // ------------------ private functions --------------------
  460.    
  461.     // for when the connectivity changes somehow; abandon graph properties
  462.     private void trashGraph()
  463.     {
  464.         graph=null;
  465.         ringID=null;
  466.         ring3=null;
  467.         ring4=null;
  468.         ring5=null;
  469.         ring6=null;
  470.         ring7=null;
  471.         compID=null;
  472.         trashPriority();
  473.     }
  474.     private void trashPriority()
  475.     {
  476.         priority=null;
  477.         trashStereo();
  478.     }
  479.     private void trashStereo()
  480.     {
  481.         chiral=null;
  482.         cistrans=null;
  483.     }
  484.    
  485.     // update the cache of atom range
  486.     private void determineMinMax()
  487.     {
  488.         invalMinMax=false;
  489.         if (numAtoms()==0)
  490.         {
  491.             minX=maxX=minY=maxY=0;
  492.             return;
  493.         }
  494.         minX=maxX=atomX(1);
  495.         minY=maxY=atomY(1);
  496.         for (int n=2;n<=numAtoms();n++)
  497.         {
  498.             double x=atomX(n),y=atomY(n);
  499.             minX=Math.min(minX,x);
  500.             maxX=Math.max(maxX,x);
  501.             minY=Math.min(minY,y);
  502.             maxY=Math.max(maxY,y);
  503.         }
  504.     }
  505.  
  506.     // update the cache of the molecular graph, in adjacency format, rather than the edge list format; the former is much better for
  507.     // graph algorithms, the latter much more suited to capturing sketcher-ish information content; note that the graph indices are
  508.     // zero-based
  509.     private void buildGraph()
  510.     {
  511.         graph=new int[numAtoms()][];
  512.         for (int n=1;n<=numBonds();n++)
  513.         {
  514.             int bf=bondFrom(n)-1,bt=bondTo(n)-1;
  515.             if (bf==bt) continue; // ouch!
  516.             int lf=graph[bf]==null ? 0 : graph[bf].length,lt=graph[bt]==null ? 0 : graph[bt].length;
  517.             int bl[]=new int[lf+1];
  518.             for (int i=0;i<lf;i++) bl[i]=graph[bf][i];
  519.             bl[lf]=bt;
  520.             graph[bf]=bl;
  521.             bl=new int[lt+1];
  522.             for (int i=0;i<lt;i++) bl[i]=graph[bt][i];
  523.             bl[lt]=bf;
  524.             graph[bt]=bl;
  525.         }
  526.         for (int n=0;n<numAtoms();n++) if (graph[n]==null) graph[n]=new int[0];
  527.        
  528.         /*for (int n=0;n<numAtoms();n++)
  529.         {
  530.             System.out.print("#"+n+":");
  531.             for (int i=0;i<graph[n].length;i++) System.out.print(" "+graph[n][i]);
  532.             System.out.println();
  533.         }*/
  534.     }
  535.    
  536.     // passes over the graph establishing which component each atom belongs to
  537.     private void buildConnComp()
  538.     {
  539.         compID=new int[numAtoms()];
  540.         for (int n=0;n<numAtoms();n++) compID[n]=0;
  541.         int comp=1;
  542.         compID[0]=comp;
  543.  
  544.         // (not very efficient, should use a stack-based depth first search)
  545.         while (true)
  546.         {
  547.             boolean anything=false;
  548.             for (int n=0;n<numAtoms();n++) if (compID[n]>0)
  549.             {
  550.                 for (int i=0;i<graph[n].length;i++) if (compID[graph[n][i]]==0)
  551.                 {
  552.                     compID[graph[n][i]]=comp;
  553.                     anything=true;
  554.                 }
  555.             }
  556.            
  557.             if (!anything)
  558.             {
  559.                 for (int n=0;n<numAtoms();n++) if (compID[n]==0) {compID[n]=++comp; anything=true; break;}
  560.                 if (!anything) break;
  561.             }
  562.         }
  563.     }
  564.  
  565.     // generates Cahn-Ingold-Prelog priority for each atom, where degeneracies are indicated by having the same number
  566.     private void buildPriority()
  567.     {
  568.         // build a graph representation which has entries replicated according to bond order
  569.         int cipgr[][]=new int[numAtoms()][];
  570.         for (int n=0;n<numAtoms();n++) if (cipgr[n]==null)
  571.         {
  572.             cipgr[n]=new int[atomHydrogens(n+1)];
  573.             for (int i=0;i<cipgr[n].length;i++) cipgr[n][i]=-1;
  574.         }
  575.         for (int n=1;n<=numBonds();n++)
  576.         {
  577.             int bf=bondFrom(n)-1,bt=bondTo(n)-1,bo=bondOrder(n);
  578.             if (bf==bt || bo==0) continue;
  579.             int lf=cipgr[bf].length,lt=cipgr[bt].length;
  580.             int bl[]=new int[lf+bo];
  581.             for (int i=0;i<lf;i++) bl[i]=cipgr[bf][i];
  582.             for (int i=0;i<bo;i++) bl[lf++]=bt;
  583.             cipgr[bf]=bl;
  584.             bl=new int[lt+bo];
  585.             for (int i=0;i<lt;i++) bl[i]=cipgr[bt][i];
  586.             for (int i=0;i<bo;i++) bl[lt++]=bf;
  587.             cipgr[bt]=bl;
  588.         }
  589.  
  590.         // seed the priorities with atomic number
  591.         priority=new int[numAtoms()];
  592.         boolean anyActualH=false;
  593.         for (int n=0;n<numAtoms();n++) {priority[n]=atomicNumber(n+1); if (priority[n]==1) anyActualH=true;}
  594.        
  595.         // pass through and reassign priorities as many times as necessary, until no change
  596.         int prigr[][]=new int[numAtoms()][];
  597.         while (true)
  598.         {
  599.             // make an equivalent to cipgr which has priorities instead of indices
  600.             for (int n=0;n<numAtoms();n++)
  601.             {
  602.                 prigr[n]=new int[cipgr[n].length];
  603.                 for (int i=0;i<prigr[n].length;i++) if (cipgr[n][i]<0) prigr[n][i]=1; else prigr[n][i]=priority[cipgr[n][i]];
  604.                 int p=0;
  605.                 while (p<prigr[n].length-1)
  606.                 {
  607.                     if (prigr[n][p]<prigr[n][p+1]) {
  608.                         int i=prigr[n][p]; prigr[n][p]=prigr[n][p+1]; prigr[n][p+1]=i;
  609.                         if (p>0) p--;
  610.                     } else p++;
  611.                 }
  612.             }
  613.            
  614.             // divide each priority category into groups, then for each of these groups, split the contents out and reassign
  615.             int groups[][]=sortAndGroup(priority);
  616.             int nextpri=anyActualH ? 0 : 1;
  617.             boolean repartitioned=false;
  618.            
  619.             for (int n=0;n<groups.length;n++)
  620.             {
  621.                 // sort the groups according to their cipgr contents
  622.                 int p=0;
  623.                 while (p<groups[n].length-1)
  624.                 {
  625.                     int i1=groups[n][p],i2=groups[n][p+1];
  626.                     int cmp=0,sz=Math.max(prigr[i1].length,prigr[i2].length);
  627.                     for (int i=0;i<sz;i++)
  628.                     {
  629.                         int v1=i<prigr[i1].length ? prigr[i1][i] : 0,v2=i<prigr[i2].length ? prigr[i2][i] : 0;
  630.                         if (v1<v2) {cmp=-1; break;}
  631.                         if (v1>v2) {cmp=1; break;}
  632.                     }
  633.                     if (cmp>0) {groups[n][p]=i2; groups[n][p+1]=i1; if (p>0) p--;}
  634.                     else p++;
  635.                 }
  636.                
  637.                 for (int i=0;i<groups[n].length;i++)
  638.                 {
  639.                     if (i==0) nextpri++;
  640.                     else if (prigr[groups[n][i]].length!=prigr[groups[n][i-1]].length) {nextpri++; repartitioned=true;}
  641.                     else
  642.                     {
  643.                         for (int j=0;j<prigr[groups[n][i]].length;j++)
  644.                             if (prigr[groups[n][i]][j]!=prigr[groups[n][i-1]][j]) {nextpri++; repartitioned=true; break;}
  645.                     }
  646.                    
  647.                     priority[groups[n][i]]=nextpri;
  648.                 }
  649.             }
  650.  
  651.             if (!repartitioned) break;
  652.         }
  653.     }
  654.  
  655.     // compute the chirality values for each atom centre
  656.     private void buildChirality()
  657.     {
  658.         chiral=new int[numAtoms()];
  659.        
  660.         boolean haswedge[]=new boolean[numAtoms()];
  661.         for (int n=0;n<numAtoms();n++) haswedge[n]=false;
  662.         for (int n=1;n<=numBonds();n++) if (bondType(n)==BONDTYPE_INCLINED || bondType(n)==BONDTYPE_DECLINED) haswedge[bondFrom(n)-1]=true;
  663.        
  664.         int pri[]=new int[4];
  665.         double x[]=new double[4],y[]=new double[4],z[]=new double[4];
  666.        
  667.         for (int n=0;n<numAtoms();n++)
  668.         {
  669.             chiral[n]=STEREO_NONE;
  670.             if (!(graph[n].length==4 || (graph[n].length==3 && atomHydrogens(n+1)==1))) continue;
  671.  
  672.             for (int i=0;i<graph[n].length;i++)
  673.             {
  674.                 pri[i]=priority[graph[n][i]];
  675.                 x[i]=atomX(graph[n][i]+1)-atomX(n+1);
  676.                 y[i]=atomY(graph[n][i]+1)-atomY(n+1);
  677.                 z[i]=0;
  678.                 if (haswedge[n])
  679.                     for (int j=1;j<=numBonds();j++) if (bondFrom(j)==n+1 && bondTo(j)==graph[n][i]+1)
  680.                 {
  681.                     if (bondType(j)==BONDTYPE_INCLINED) z[i]=1;
  682.                     if (bondType(j)==BONDTYPE_DECLINED) z[i]=-1;
  683.                     break;
  684.                 }
  685.             }
  686.             if (graph[n].length==3) {pri[3]=0; x[3]=0; y[3]=0; z[3]=0;}
  687.            
  688.             int p=0;
  689.             while (p<3)
  690.             {
  691.                 if (pri[p]>pri[p+1])
  692.                 {
  693.                     int i; double d;
  694.                     i=pri[p]; pri[p]=pri[p+1]; pri[p+1]=i;
  695.                     d=x[p]; x[p]=x[p+1]; x[p+1]=d;
  696.                     d=y[p]; y[p]=y[p+1]; y[p+1]=d;
  697.                     d=z[p]; z[p]=z[p+1]; z[p+1]=d;
  698.                     if (p>0) p--;
  699.                 }
  700.                 else p++;
  701.             }
  702.             if ((pri[0]==0 && pri[1]==1) || pri[0]==pri[1] || pri[1]==pri[2] || pri[2]==pri[3]) continue; // no topological chirality
  703.            
  704.             chiral[n]=STEREO_UNKNOWN;
  705.             if (z[0]==0 && z[1]==0 && z[2]==0 && z[3]==0) continue; // all atoms are in-plane
  706.            
  707.             boolean sane=true;
  708.             for (int i=0;i<4;i++) if (pri[0]!=0)
  709.             {
  710.                 double r=x[i]*x[i]+y[i]*y[i]+z[i]*z[i];
  711.                 if (r<0.01*0.01) {sane=false; break;}
  712.                 r=1/Math.sqrt(r);
  713.                 x[i]=x[i]*r; y[i]=y[i]*r; z[i]=z[i]*r;
  714.             }
  715.             if (!sane) continue;
  716.            
  717.             if (pri[0]==0) // build a position for the implicit H
  718.             {
  719.                 x[0]=-(x[1]+x[2]+x[3]);
  720.                 y[0]=-(y[1]+y[2]+y[3]);
  721.                 z[0]=-(z[1]+z[2]+z[3]);
  722.                 double r=x[0]*x[0]+y[0]*y[0]+z[0]*z[0];
  723.                 if (r<0.01*0.01) {sane=false; break;}
  724.                 r=1/Math.sqrt(r);
  725.                 x[0]=x[0]*r; y[0]=y[0]*r; z[0]=z[0]*r;
  726.             }
  727.             if (!sane) continue;
  728.            
  729.             double R=0,S=0;
  730.             for (int i=1;i<=6;i++)
  731.             {
  732.                 int a=0,b=0;
  733.                 if      (i==1) {a=1; b=2;} else if (i==2) {a=2; b=3;} else if (i==3) {a=3; b=1;}
  734.                 else if (i==4) {a=2; b=1;} else if (i==5) {a=3; b=2;} else if (i==6) {a=1; b=3;}
  735.                 double xx=y[a]*z[b]-y[b]*z[a] - x[0];
  736.                 double yy=z[a]*x[b]-z[b]*x[a] - y[0];
  737.                 double zz=x[a]*y[b]-x[b]*y[a] - z[0];
  738.                 if (i<=3) R+=xx*xx+yy*yy+zz*zz; else S+=xx*xx+yy*yy+zz*zz;
  739.             }
  740.             chiral[n]=R>S ? STEREO_POS : STEREO_NEG;
  741.         }
  742.     }
  743.    
  744.     // computer the cis/trans stereochemistry for each bond
  745.     private void buildCisTrans()
  746.     {
  747.         cistrans=new int[numBonds()];
  748.         int sf[]=new int[2],st[]=new int[2];
  749.        
  750.         for (int n=0;n<numBonds();n++)
  751.         {
  752.             cistrans[n]=STEREO_NONE;
  753.             int bf=bondFrom(n+1)-1,bt=bondTo(n+1)-1;
  754.             if (bondOrder(n+1)!=2 || graph[bf].length<=1 || graph[bt].length<= 1 || graph[bf].length>3 || graph[bt].length>3) continue;
  755.  
  756.             int nf=0,nt=0;
  757.             for (int i=0;i<graph[bf].length;i++) if (graph[bf][i]!=bt) sf[nf++]=graph[bf][i];
  758.             for (int i=0;i<graph[bt].length;i++) if (graph[bt][i]!=bf) st[nt++]=graph[bt][i];
  759.  
  760.             if (nf==1) {if (atomHydrogens(bf+1)!=1 || priority[sf[0]]==1) continue;}
  761.             else
  762.             {
  763.                 if (priority[sf[0]]==priority[sf[1]]) continue;
  764.                 if (priority[sf[0]]<priority[sf[1]]) {int i=sf[0]; sf[0]=sf[1]; sf[1]=i;}
  765.             }
  766.            
  767.             if (nt==1) {if (atomHydrogens(bt+1)!=1 || priority[st[0]]==1) continue;}
  768.             else
  769.             {
  770.                 if (priority[st[0]]==priority[st[1]]) continue;
  771.                 if (priority[st[0]]<priority[st[1]]) {int i=st[0]; st[0]=st[1]; st[1]=i;}
  772.             }
  773.        
  774.             cistrans[n]=STEREO_UNKNOWN;
  775.  
  776.             double xa=atomX(bf+1),ya=atomY(bf+1),xb=atomX(bt+1),yb=atomY(bt+1);
  777.             double tha0=Math.atan2(yb-ya,xb-xa),thb0=Math.PI+tha0;
  778.             double tha1=Math.atan2(atomY(sf[0]+1)-ya,atomX(sf[0]+1)-xa);
  779.             double tha2=nf==2 ? Math.atan2(atomY(sf[1]+1)-ya,atomX(sf[1]+1)-xa) : thetaObtuse(tha0,tha1);
  780.             double thb1=Math.atan2(atomY(st[0]+1)-yb,atomX(st[0]+1)-xb);
  781.             double thb2=nt==2 ? Math.atan2(atomY(st[1]+1)-yb,atomX(st[1]+1)-xb) : thetaObtuse(thb0,thb1);
  782.             tha1-=tha0; tha2-=tha0; thb1-=thb0; thb2-=thb0;
  783.             tha1+=(tha1<-Math.PI ? 2*Math.PI : 0)+(tha1>Math.PI ? -2*Math.PI : 0);
  784.             tha2+=(tha2<-Math.PI ? 2*Math.PI : 0)+(tha2>Math.PI ? -2*Math.PI : 0);
  785.             thb1+=(thb1<-Math.PI ? 2*Math.PI : 0)+(thb1>Math.PI ? -2*Math.PI : 0);
  786.             thb2+=(thb2<-Math.PI ? 2*Math.PI : 0)+(thb2>Math.PI ? -2*Math.PI : 0);
  787.  
  788.             final double SMALL=5*Math.PI/180; // basically colinear
  789.             if (Math.abs(tha1)<SMALL || Math.abs(tha2)<SMALL || Math.abs(thb1)<SMALL || Math.abs(thb2)<SMALL) continue;
  790.             if (Math.abs(tha1)>Math.PI-SMALL || Math.abs(tha2)>Math.PI-SMALL) continue;
  791.             if (Math.abs(thb1)>Math.PI-SMALL || Math.abs(thb2)>Math.PI-SMALL) continue;
  792.             tha1=Math.signum(tha1); tha2=Math.signum(tha2); thb1=Math.signum(thb1); thb2=Math.signum(thb2);
  793.             if (tha1==tha2 || thb1==thb2) continue;
  794.             if (tha1*thb1<0) cistrans[n]=STEREO_POS;
  795.             if (tha1*thb1>0) cistrans[n]=STEREO_NEG;
  796.         }
  797.     }
  798.  
  799.     // generally useful function which takes a list of numbers and sorts them, then bins the unique values into sub-arrays
  800.     // (note: inefficient implementation, but could be improved easily enough)
  801.     public static int[][] sortAndGroup(int[] Val)
  802.     {
  803.         ArrayList<Integer> unique=new ArrayList<Integer>();
  804.         for (int n=0;n<Val.length;n++) unique.add(Val[n]);
  805.         Collections.sort(unique);
  806.         int p=0;
  807.         while (p<unique.size()-1)
  808.         {
  809.             if (unique.get(p)==unique.get(p+1)) unique.remove(p); else p++;
  810.         }
  811.  
  812.         int ret[][]=new int[unique.size()][];
  813.        
  814.         for (int n=0;n<Val.length;n++)
  815.         {
  816.             int grp=unique.indexOf(Val[n]);
  817.             int[] cnt=new int[ret[grp]==null ? 1 : ret[grp].length+1];
  818.             for (int i=0;i<cnt.length-1;i++) cnt[i]=ret[grp][i];
  819.             cnt[cnt.length-1]=n;
  820.             ret[grp]=cnt;
  821.         }
  822.        
  823.         return ret;
  824.     }
  825.  
  826.     // update the ring-block-identifier for each atom
  827.     private void buildRingID()
  828.     {
  829.         ringID=new int[numAtoms()];
  830.         if (numAtoms()==0) return;
  831.         boolean visited[]=new boolean[numAtoms()];
  832.         for (int n=0;n<numAtoms();n++)
  833.         {
  834.             ringID[n]=0;
  835.             visited[n]=false;
  836.         }
  837.  
  838.         int path[]=new int[numAtoms()+1],plen=0,numVisited=0;
  839.         boolean rewound=false;
  840.         while (true)
  841.         {
  842.             int last,current;
  843.  
  844.             if (plen==0) // find an element of a new component to visit
  845.             {
  846.                 last=-1;
  847.                 for (current=0;visited[current];current++) {}
  848.             }
  849.             else
  850.             {
  851.                 last=path[plen-1];
  852.                 current=-1;
  853.                 for (int n=0;n<graph[last].length;n++) if (!visited[graph[last][n]]) {current=graph[last][n]; break;}
  854.             }
  855.            
  856.             /*System.out.print("numVisited="+numVisited+" last="+last+" current="+current+" path=");
  857.             for (int n=0;n<plen;n++) System.out.print(path[n]+" ");
  858.             System.out.println();*/
  859.            
  860.             if (current>=0 && plen>=2) // path is at least 2 items long, so look for any not-previous visited neighbours
  861.             {
  862.                 int back=path[plen-1];
  863.                 //System.out.println(" back="+back);
  864.                 for (int n=0;n<graph[current].length;n++)
  865.                 {
  866.                     int join=graph[current][n];
  867.                     //System.out.println(" join="+join+" visited[join]="+visited[join]);
  868.                     if (join!=back && visited[join])
  869.                     {
  870.                         //System.out.print(" path:");
  871.  
  872.                         path[plen]=current;
  873.                         for (int i=plen;i==plen || path[i+1]!=join;i--)
  874.                         {
  875.                             //System.out.print(" "+path[i]);
  876.  
  877.                             int id=ringID[path[i]];
  878.                             if (id==0) ringID[path[i]]=last;
  879.                             else if (id!=last)
  880.                             {
  881.                                 for (int j=0;j<numAtoms();j++) if (ringID[j]==id) ringID[j]=last;
  882.                             }
  883.                         }
  884.                         //System.out.println();
  885.                     }
  886.                 }
  887.             }
  888.             if (current>=0)  // can mark the new one as visited
  889.             {
  890.                 visited[current]=true;
  891.                 path[plen++]=current;
  892.                 numVisited++;
  893.                 rewound=false;
  894.             }
  895.             else // otherwise, found nothing and must rewind the path
  896.             {
  897.                 plen--;
  898.                 rewound=true;
  899.             }
  900.            
  901.             if (numVisited==numAtoms()) break;
  902.         }
  903.        
  904.         // the ring ID's are not necessarily consecutive, so reassign them to 0=none, 1..NBlocks
  905.         int nextID=0;
  906.         for (int i=0;i<numAtoms();i++) if (ringID[i]>0)
  907.         {
  908.             nextID--;
  909.             for (int j=numAtoms()-1;j>=i;j--) if (ringID[j]==ringID[i]) ringID[j]=nextID;
  910.         }
  911.         for (int i=0;i<numAtoms();i++) ringID[i]=-ringID[i];
  912.     }
  913.    
  914.     // ring hunter: recursive step; finds, compares and collects
  915.     private void recursiveRingFind(int[] Path,int PSize,int Capacity,int RBlk,ArrayList<Object> Rings)
  916.     {
  917.         // not enough atoms yet, so look for new possibilities
  918.         if (PSize<Capacity)
  919.         {
  920.             int last=Path[PSize-1];
  921.             for (int n=0;n<graph[last-1].length;n++)
  922.             {
  923.                 int adj=graph[last-1][n]+1;
  924.                 if (ringID[adj-1]!=RBlk) continue;
  925.                 boolean fnd=false;
  926.                 for (int i=0;i<PSize;i++) if (Path[i]==adj) {fnd=true; break;}
  927.                 if (!fnd)
  928.                 {
  929.                     int newPath[]=Path.clone();
  930.                     newPath[PSize]=adj;
  931.                     recursiveRingFind(newPath,PSize+1,Capacity,RBlk,Rings);
  932.                 }
  933.             }
  934.             return;
  935.         }
  936.        
  937.         // path is full, so make sure it eats its tail
  938.         int last=Path[PSize-1];
  939.         boolean fnd=false;
  940.         for (int n=0;n<graph[last-1].length;n++) if (graph[last-1][n]+1==Path[0]) {fnd=true; break;}
  941.         if (!fnd) return;
  942.        
  943.         // reorder the array then look for duplicates
  944.         int first=0;
  945.         for (int n=1;n<PSize;n++) if (Path[n]<Path[first]) first=n;
  946.         int fm=(first-1+PSize)%PSize,fp=(first+1)%PSize;
  947.         boolean flip=Path[fm]<Path[fp];
  948.         if (first!=0 || flip)
  949.         {
  950.             int newPath[]=new int[PSize];
  951.             for (int n=0;n<PSize;n++) newPath[n]=Path[(first+(flip ? PSize-n : n))%PSize];
  952.             Path=newPath;
  953.         }
  954.                
  955.         for (int n=0;n<Rings.size();n++)
  956.         {
  957.             int[] look=(int[])Rings.get(n);
  958.             boolean same=true;
  959.             for (int i=0;i<PSize;i++) if (look[i]!=Path[i]) {same=false; break;}
  960.             if (same) return;
  961.         }
  962.        
  963.         Rings.add(Path);
  964.     }
  965.    
  966.     // returns the angle maximally equidistant from Th1 and Th2
  967.     private double thetaObtuse(double Th1,double Th2)
  968.     {
  969.         double dth=Th2-Th1;
  970.         while (dth<-Math.PI) dth+=2*Math.PI;
  971.         while (dth>Math.PI) dth-=2*Math.PI;
  972.         return dth>0 ? Th1-0.5*(2*Math.PI-dth) : Th1+0.5*(2*Math.PI+dth);
  973.     }
  974. }
  975.