/*
 *                  BioJava development code
 *
 * This code may be freely distributed and modified under the
 * terms of the GNU Lesser General Public Licence.  This should
 * be distributed with the code.  If you do not have a copy,
 * see:
 *
 *      http://www.gnu.org/copyleft/lesser.html
 *
 * Copyright for this code is held jointly by the individual
 * authors.  These should be listed in @author doc comments.
 *
 * For more information on the BioJava project and its aims,
 * or to join the biojava-l mailing list, visit the home page
 * at:
 *
 *      http://www.biojava.org/
 *
 * Created on Jan 4, 2006
 *
 */
package org.biojava.bio.structure;

import java.io.StringWriter;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.HashSet;
import java.util.Iterator;
import java.util.LinkedHashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;

import java.util.logging.Level;
import java.util.logging.Logger;
import java.util.regex.Matcher;
import java.util.regex.Pattern;

import org.biojava.bio.structure.align.util.AtomCache;


/**
 * A class that provides some tool methods.
 *
 * @author Andreas Prlic, Jules Jacobsen
 * @since 1.0
 * @version %I% %G%
 */
public class StructureTools {

        /** The Atom name of C-alpha atoms.
         *
         */
        public static final String caAtomName = " CA ";

        public static final String nAtomName = "N";

        public static final String oAtomName = "O";

        public static final String cbAtomName = "CB";


        /** The names of the Atoms that form the backbone.
         *
         */
        public static final String[] backboneAtomNames = {nAtomName,caAtomName,"C",oAtomName, cbAtomName};

        public static final Character UNKNOWN_GROUP_LABEL = new Character('x');;


        //private static final String insertionCodeRegExp = "([0-9]+)([a-zA-Z]*)";
        //private static final Pattern insertionCodePattern = Pattern.compile(insertionCodeRegExp);


        // there is a file format change in PDB 3.0 and nucleotides are being renamed
        static private Map<String, Integer> nucleotides30 ;
        static private Map<String, Integer> nucleotides23 ;

        //amino acid 3 and 1 letter code definitions
        private static final Map<String, Character> aminoAcids;

        private static final Set<Element> hBondDonorAcceptors;
        //      // for conversion 3code 1code
        //      private static  SymbolTokenization threeLetter ;
        //      private static  SymbolTokenization oneLetter ;

        public static Logger logger =  Logger.getLogger("org.biojava.bio.structure");

        /**
         * Pattern to describe subranges. Matches "A", "A:", "A:7-53","A_7-53", etc.
         * @see #getSubRanges(Structure, String)
         */
        public static final Pattern pdbNumRangeRegex = Pattern.compile(
                        "^\\s*(\\w)" + //chain ID
                        "(?:" + //begin range, this is a "non-capturing group"
                          "(?::|_|:$|_$|$)" + //colon or underscore, could be at the end of a line, another non-capt. group.
                                "(?:"+ // another non capturing group for the residue range
                                 "([-+]?[0-9]+[A-Za-z]?)" + // first residue
                                 "\\s*-\\s*" + // -
                                 "([-+]?[0-9]+[A-Za-z]?)" + // second residue
                                ")?+"+
                        ")?" + //end range
                        "\\s*");
        
        
        static {
                nucleotides30 = new HashMap<String,Integer>();
                nucleotides30.put("DA",1);
                nucleotides30.put("DC",1);
                nucleotides30.put("DG",1);
                nucleotides30.put("DT",1);
                nucleotides30.put("DI",1);
                nucleotides30.put("A",1);
                nucleotides30.put("G",1);
                nucleotides30.put("C",1);
                nucleotides30.put("U",1);
                nucleotides30.put("I",1);

                //TODO: check if they are always HETATMs, in that case this will not be necessary
                // the DNA linkers - the +C , +G, +A  +T +U and +I have been replaced with these:
                nucleotides30.put("TAF",1); // 2'-DEOXY-2'-FLUORO-ARABINO-FURANOSYL THYMINE-5'-PHOSPHATE
                nucleotides30.put("TC1",1); // 3-(5-PHOSPHO-2-DEOXY-BETA-D-RIBOFURANOSYL)-2-OXO-1,3-DIAZA-PHENOTHIAZINE
                nucleotides30.put("TFE",1); // 2'-O-[2-(TRIFLUORO)ETHYL] THYMIDINE-5'-MONOPHOSPHATE
                nucleotides30.put("TFO",1); // [2-(6-AMINO-9H-PURIN-9-YL)-1-METHYLETHOXY]METHYLPHOSPHONIC ACID"
                nucleotides30.put("TGP",1); // 5'-THIO-2'-DEOXY-GUANOSINE PHOSPHONIC ACID
                nucleotides30.put("THX",1); // PHOSPHONIC ACID 6-({6-[6-(6-CARBAMOYL-3,6,7,8-TETRAHYDRO-3,6-DIAZA-AS-INDACENE-2-CARBONYL)-3,6,7,8-TETRAHYDRO-3,6-DIAZA-AS-INDOCENE-2-CARBONYL]-3,6,7,8-TETRAHYDRO-3,6-DIAZA-AS-INDACENE-2-CARBONL}-AMINO)-HEXYL ESTER 5-(5-METHYL-2,4-DIOXO-3,4-DIHYDRO-2H-PYRIMIDIN-1-YL)-TETRAHYDRO-FURAN-2-YLMETHYL ESTER
                nucleotides30.put("TLC",1); // 2-O,3-ETHDIYL-ARABINOFURANOSYL-THYMINE-5'-MONOPHOSPHATE
                nucleotides30.put("TLN",1); //  [(1R,3R,4R,7S)-7-HYDROXY-3-(THYMIN-1-YL)-2,5-DIOXABICYCLO[2.2.1]HEPT-1-YL]METHYL DIHYDROGEN PHOSPHATE"
                nucleotides30.put("TP1",1); // 2-(METHYLAMINO)-ETHYLGLYCINE-CARBONYLMETHYLENE-THYMINE
                nucleotides30.put("TPC",1); // 5'-THIO-2'-DEOXY-CYTOSINE PHOSPHONIC ACID
                nucleotides30.put("TPN",1); // 2-AMINOETHYLGLYCINE-CARBONYLMETHYLENE-THYMINE



                // store nucleic acids (C, G, A, T, U, and I), and
                // the modified versions of nucleic acids (+C, +G, +A, +T, +U, and +I), and
                nucleotides23  = new HashMap<String,Integer>();
                String[] names = {"C","G","A","T","U","I","+C","+G","+A","+T","+U","+I"};
                for (int i = 0; i < names.length; i++) {
                        String n = names[i];
                        nucleotides23.put(n,1);
                }

                aminoAcids = new HashMap<String, Character>();
                aminoAcids.put("GLY", new Character('G'));
                aminoAcids.put("ALA", new Character('A'));
                aminoAcids.put("VAL", new Character('V'));
                aminoAcids.put("LEU", new Character('L'));
                aminoAcids.put("ILE", new Character('I'));
                aminoAcids.put("PHE", new Character('F'));
                aminoAcids.put("TYR", new Character('Y'));
                aminoAcids.put("TRP", new Character('W'));
                aminoAcids.put("PRO", new Character('P'));
                aminoAcids.put("HIS", new Character('H'));
                aminoAcids.put("LYS", new Character('K'));
                aminoAcids.put("ARG", new Character('R'));
                aminoAcids.put("SER", new Character('S'));
                aminoAcids.put("THR", new Character('T'));
                aminoAcids.put("GLU", new Character('E'));
                aminoAcids.put("GLN", new Character('Q'));
                aminoAcids.put("ASP", new Character('D'));
                aminoAcids.put("ASN", new Character('N'));
                aminoAcids.put("CYS", new Character('C'));
                aminoAcids.put("MET", new Character('M'));
                //MSE is only found as a molecular replacement for MET
                aminoAcids.put("MSE", new Character('M'));
                //'non-standard', genetically encoded
                //http://www.chem.qmul.ac.uk/iubmb/newsletter/1999/item3.html
                //IUBMB recommended name is 'SEC' but the wwPDB currently use 'CSE'
                //likewise 'PYL' (IUBMB) and 'PYH' (PDB)
                aminoAcids.put("CSE", new Character('U'));
                aminoAcids.put("SEC", new Character('U'));
                aminoAcids.put("PYH", new Character('O'));
                aminoAcids.put("PYL", new Character('O'));

                hBondDonorAcceptors = new HashSet<Element>();
                hBondDonorAcceptors.add(Element.N);
                hBondDonorAcceptors.add(Element.O);
                hBondDonorAcceptors.add(Element.S);

        }


        /** Count how many number of Atoms are contained within a Structure object.
         *
         * @param s the structure object
         * @return the number of Atoms in this Structure
         */
        public static final int getNrAtoms(Structure s){

                int nrAtoms = 0;

                Iterator<Group> iter = new GroupIterator(s);

                while ( iter.hasNext()){
                        Group g = (Group) iter.next();
                        nrAtoms += g.size();
                }

                return nrAtoms;
        }


        /** Count how many groups are contained within a structure object.
         *
         * @param s the structure object
         * @return the number of groups in the structure
         */
        public static final int getNrGroups(Structure s){
                int nrGroups = 0;

                List<Chain> chains = s.getChains(0);
                Iterator<Chain> iter = chains.iterator();
                while (iter.hasNext()){
                        Chain c = (Chain) iter.next();
                        nrGroups += c.getAtomLength();
                }
                return nrGroups;
        }


        /** Returns an array of the requested Atoms from the Structure object. Iterates over all groups
         * and checks if the requested atoms are in this group, no matter if this is a  {@link AminoAcid} or {@link HetatomImpl} group.
         * For structures with more than one model, only model 0 will be used.
         *
         * @param s the structure to get the atoms from
         *
         * @param atomNames  contains the atom names to be used.
         * @return an Atom[] array
         */
        public static final Atom[] getAtomArray(Structure s, String[] atomNames){
                List<Chain> chains = s.getModel(0);

                List<Atom> atoms = new ArrayList<Atom>();

                extractCAatoms(atomNames, chains, atoms);
                return (Atom[]) atoms.toArray(new Atom[atoms.size()]);

        }
        
        /** Returns an array of the requested Atoms from the Structure object. 
         * In contrast to {@link #getAtomArray(Structure, String[])} this method iterates over all chains.
         * Iterates over all chains and groups
         * and checks if the requested atoms are in this group, no matter if this is a {@link AminoAcid} or {@link HetatomImpl} group.
         * For structures with more than one model, only model 0 will be used.
         *
         * @param s the structure to get the atoms from
         *
         * @param atomNames  contains the atom names to be used.
         * @return an Atom[] array
         */
        public static final Atom[] getAtomArrayAllModels(Structure s, String[] atomNames){
                
                List<Atom> atoms = new ArrayList<Atom>();

                for (int i =0 ; i < s.nrModels(); i++ ) {
                        List<Chain> chains = s.getModel(i);
                        extractCAatoms(atomNames, chains, atoms);
                }
                return (Atom[]) atoms.toArray(new Atom[atoms.size()]);

        }



        private static void extractCAatoms(String[] atomNames, List<Chain> chains,
                        List<Atom> atoms) {
                for ( Chain c : chains) {
                        
                        for ( Group g : c.getAtomGroups()) {

                                // a temp container for the atoms of this group
                                List<Atom> thisGroupAtoms = new ArrayList<Atom>();
                                // flag to check if this group contains all the requested atoms.
                                boolean thisGroupAllAtoms = true;
                                for ( int i = 0 ; i < atomNames.length; i++){
                                        String atomName = atomNames[i];
                                        try {
                                                Atom a = g.getAtom(atomName);
                                                thisGroupAtoms.add(a);
                                        } catch (StructureException e){
                                                // this group does not have a required atom, skip it...
                                                thisGroupAllAtoms = false;
                                                break;
                                        }
                                }
                                if ( thisGroupAllAtoms){
                                        // add the atoms of this group to the array.
                                        Iterator<Atom> aIter = thisGroupAtoms.iterator();
                                        while(aIter.hasNext()){
                                                Atom a = (Atom) aIter.next();
                                                atoms.add(a);
                                        }
                                }

                        }
                }
        }

        /** Returns an array of the requested Atoms from the Structure object. Iterates over all groups
         * and checks if the requested atoms are in this group, no matter if this is a AminoAcid or Hetatom group.
         *
         *
         * @param c the Chain to get the atoms from
         *
         * @param atomNames  contains the atom names to be used.
         * @return an Atom[] array
         */
        public static final Atom[] getAtomArray(Chain c, String[] atomNames){

                List<Group> groups = c.getAtomGroups();

                List<Atom> atoms = new ArrayList<Atom>();

                for (Group g : groups){

                        // a temp container for the atoms of this group
                        List<Atom> thisGroupAtoms = new ArrayList<Atom>();
                        // flag to check if this group contains all the requested atoms.
                        boolean thisGroupAllAtoms = true;
                        for ( int i = 0 ; i < atomNames.length; i++){
                                String atomName = atomNames[i];
                                try {
                                        Atom a = g.getAtom(atomName);
                                        thisGroupAtoms.add(a);
                                } catch (StructureException e){
                                        // this group does not have a required atom, skip it...
                                        thisGroupAllAtoms = false;
                                        break;
                                }
                        }
                        if ( thisGroupAllAtoms){
                                // add the atoms of this group to the array.
                                Iterator<Atom> aIter = thisGroupAtoms.iterator();
                                while(aIter.hasNext()){
                                        Atom a = (Atom) aIter.next();
                                        atoms.add(a);
                                }
                        }

                }
                return (Atom[]) atoms.toArray(new Atom[atoms.size()]);

        }

        /** Returns an Atom array of the CA atoms.
         * @param c the structure object
         * @return an Atom[] array
         */
        public static final Atom[] getAtomCAArray(Chain c){
                String[] atomNames = {" CA " };
                return getAtomArray(c,atomNames);
        }

        /** Provides an equivalent copy of Atoms in a new array. Clones everything, starting with parent 
         * groups and chains. The chain will only contain groups that are part of the CA array.
         * 
         * @param ca array of CA atoms
         * @return Atom array
         */
        public static final Atom[] cloneCAArray(Atom[] ca) throws StructureException{
                Atom[] newCA = new Atom[ca.length];

                List<Chain> model = new ArrayList<Chain>();
                int apos = -1;
                for(Atom a: ca){
                        apos++;
                        Group parentG = a.getGroup();
                        Chain parentC = parentG.getChain();

                        Chain newChain = null;
                        for ( Chain c : model){
                                if ( c.getChainID().equals(parentC.getChainID())){
                                        newChain = c;
                                        break;
                                }
                        }
                        if ( newChain == null){
                                newChain = new ChainImpl();
                                newChain.setChainID(parentC.getChainID());
                                model.add(newChain);
                        }

                        Group parentN = (Group)parentG.clone();

                        newCA[apos] = parentN.getAtom(" CA ");
                        newChain.addGroup(parentN);
                }
                return newCA;
        }

        /** Clone a set of CA Atoms, but returns the parent groups
         *  
         * @param ca Atom array
         * @return Group array
         */
        public static Group[] cloneGroups(Atom[] ca) {
                Group[] newGroup = new Group[ca.length]; 

                List<Chain> model = new ArrayList<Chain>();
                int apos = -1;
                for(Atom a: ca){
                        apos++;
                        Group parentG = a.getGroup();
                        Chain parentC = parentG.getChain();

                        Chain newChain = null;
                        for ( Chain c : model){
                                if ( c.getChainID().equals(parentC.getChainID())){
                                        newChain = c;
                                        break;
                                }
                        }
                        if ( newChain == null){
                                newChain = new ChainImpl();
                                newChain.setChainID(parentC.getChainID());
                                model.add(newChain);
                        }

                        Group ng = (Group)parentG.clone();
                        newGroup[apos] = ng;
                        newChain.addGroup(ng);
                }
                return newGroup;
        }
        
        /** Utility method for working with circular permutations. Creates a duplicated and cloned set of Calpha atoms from the input array.
         * 
         * @param ca2 atom array
         * @return cloned and duplicated set of input array
         * @throws StructureException
         */
        public static Atom[] duplicateCA2(Atom[] ca2) throws StructureException{
                // we don't want to rotate input atoms, do we?
                Atom[] ca2clone = new Atom[ca2.length*2];

                int pos = 0;

                Chain c = null;
                String prevChainId = "";
                for (Atom a : ca2){                     
                        Group g = (Group) a.getGroup().clone(); // works because each group has only a CA atom
                        
                        if (c == null ) {
                                c = new ChainImpl();
                                Chain orig= a.getGroup().getChain();
                                c.setChainID(orig.getChainID());
                        } else {
                                Chain orig= a.getGroup().getChain();
                                if ( ! orig.getChainID().equals(prevChainId)){
                                        c = new ChainImpl();
                                        c.setChainID(orig.getChainID());
                                }
                        }
                        
                        c.addGroup(g);
                        ca2clone[pos] = g.getAtom(StructureTools.caAtomName);

                        pos++;
                }

                // Duplicate ca2!
                c = null;
                prevChainId = "";
                for (Atom a : ca2){
                        Group g = (Group)a.getGroup().clone();
                        
                        if (c == null ) {
                                c = new ChainImpl();
                                Chain orig= a.getGroup().getChain();
                                c.setChainID(orig.getChainID());
                        } else {
                                Chain orig= a.getGroup().getChain();
                                if ( ! orig.getChainID().equals(prevChainId)){
                                        c = new ChainImpl();
                                        c.setChainID(orig.getChainID());
                                }
                        }
                        
                        c.addGroup(g);
                        ca2clone[pos] = g.getAtom(StructureTools.caAtomName);

                        pos++;
                }

                return ca2clone;

        }

        

        /** Returns an Atom array of the CA atoms.
         * @param s the structure object
         * @return an Atom[] array
         */
        public static Atom[] getAtomCAArray(Structure s){
                String[] atomNames = {" CA "};
                return getAtomArray(s,atomNames);
        }

        /** Returns an Atom array of the MainChain atoms.

         * @param s the structure object
         * @return an Atom[] array
         */
        public static Atom[] getBackboneAtomArray(Structure s){
                String[] atomNames = backboneAtomNames;
                return getAtomArray(s,atomNames);
        }


        /** convert three character amino acid codes into single character
         *  e.g. convert CYS to C
         *  @return a character
         *  @param code3 a three character amino acid representation String
         *  @throws IllegalSymbolException
         */

        public static final Character convert_3code_1code(String code3)
        throws UnknownPdbAminoAcidException {
                //      {
                //              Symbol sym   =  threeLetter.parseToken(code3) ;
                //              String code1 =  oneLetter.tokenizeSymbol(sym);
                //
                //              return new Character(code1.charAt(0)) ;
                Character code1 = null;
                code1 = aminoAcids.get(code3);

                if (code1 == null) {
                        throw new UnknownPdbAminoAcidException(code3 + " not a standard amino acid");
                } else {
                        return code1;
                }

        }

        /** convert a three letter code into single character.
         * catches for unusual characters
         *
         * @param groupCode3 three letter representation
         * @return null if group is a nucleotide code
         */
        public static final Character get1LetterCode(String groupCode3){

                Character aminoCode1 = null;
                try {
                        // is it a standard amino acid ?
                        aminoCode1 = convert_3code_1code(groupCode3);
                } catch (UnknownPdbAminoAcidException e){
                        // hm groupCode3 is not standard
                        // perhaps it is an nucleotide?
                        if ( isNucleotide(groupCode3) ) {
                                //System.out.println("nucleotide, aminoCode1:"+aminoCode1);
                                aminoCode1= null;
                        } else {
                                // does not seem to be so let's assume it is
                                //  nonstandard aminoacid and label it "X"
                                //logger.warning("unknown group name "+groupCode3 );
                                aminoCode1 = UNKNOWN_GROUP_LABEL;
                        }
                }

                return aminoCode1;

        }


        /* Test if the threelettercode of an ATOM entry corresponds to a
         * nucleotide or to an aminoacid.
         * @param a 3-character code for a group.
         *
         */
        public static final boolean isNucleotide(String groupCode3){

                String code = groupCode3.trim();
                if ( nucleotides30.containsKey(code)){
                        return true;
                }

                if ( nucleotides23.containsKey(code)){
                        return true;
                }

                return false ;
        }

        /** Reduce a structure to provide a smaller representation . Only takes the first model of the structure. If chainId is provided only return a structure containing that Chain ID. 
         * Converts lower case chain IDs to upper case if structure does not contain a chain with that ID. 
         * 
         * @param s
         * @param chainId
         * @return Structure
         * @since 3.0
         */
        @SuppressWarnings("deprecation")
        public static final Structure getReducedStructure(Structure s, String chainId) throws StructureException{
                // since we deal here with structure alignments,
                // only use Model 1...

                Structure newS = new StructureImpl();
                newS.setHeader(s.getHeader());
                newS.setPDBCode(s.getPDBCode());
                newS.setPDBHeader(s.getPDBHeader());
                newS.setName(s.getName());
                newS.setSSBonds(s.getSSBonds());
                newS.setDBRefs(s.getDBRefs());
                newS.setSites(s.getSites());
                newS.setNmr(s.isNmr());
                newS.setBiologicalAssembly(s.isBiologicalAssembly());
                newS.setCompounds(s.getCompounds());
                newS.setConnections(s.getConnections());
                newS.setSSBonds(s.getSSBonds());
                newS.setSites(s.getSites());
                
                if ( chainId != null)
                        chainId = chainId.trim();

                if ( chainId == null || chainId.equals("")){
                        // only get model 0
                        List<Chain> model0 = s.getModel(0);
                        for (Chain c : model0){
                                newS.addChain(c);
                        }
                        return newS;

                }

                Chain c =  null;
                try {
                        c = s.getChainByPDB(chainId);
                } catch (StructureException e){
                        System.err.println(e.getMessage() + " trying upper case Chain id...");
                        c = s.getChainByPDB(chainId.toUpperCase());
                        
                        
                }
                if ( c != null) {
                        newS.addChain(c);
                        for ( Compound comp : s.getCompounds()){
                                if ( comp.getChainId().contains(c.getChainID())){
                                        // found matching compound. set description...
                                        newS.getPDBHeader().setDescription("Chain " + c.getChainID() + " of " + s.getPDBCode() + " " + comp.getMolName());
                                }
                        }
                }


                return newS;
        }


        /** Reduce a structure to provide a smaller representation.
         * Only takes the first model of the structure. If chainNr >=0 only takes
         * the chain at that position into account.
         * 
         * @param s
         * @param chainNr can be -1 to request all chains of model 0, otherwise will only add chain at this position 
         * @return Structure object
         * @since 3.0
         */
        @SuppressWarnings("deprecation")
        public static final Structure getReducedStructure(Structure s, int chainNr) throws StructureException{
                // since we deal here with structure alignments,
                // only use Model 1...

                Structure newS = new StructureImpl();
                newS.setHeader(s.getHeader());
                newS.setPDBCode(s.getPDBCode());
                newS.setPDBHeader(s.getPDBHeader());
                newS.setName(s.getName());
                newS.setSSBonds(s.getSSBonds());
                newS.setDBRefs(s.getDBRefs());
                newS.setSites(s.getSites());
                newS.setNmr(s.isNmr());
                newS.setBiologicalAssembly(s.isBiologicalAssembly());
                newS.setCompounds(s.getCompounds());
                newS.setConnections(s.getConnections());
                newS.setSSBonds(s.getSSBonds());
                newS.setSites(s.getSites());
                newS.getPDBHeader().setDescription("subset of " + s.getPDBCode() + " " + s.getPDBHeader().getDescription() );
                
                if ( chainNr < 0 ) {

                        // only get model 0
                        List<Chain> model0 = s.getModel(0);
                        for (Chain c : model0){
                                newS.addChain(c);
                        }
                        return newS;
                }
                Chain c =  null;

                c = s.getChain(0, chainNr);

                newS.addChain(c);

                return newS;
        }

        
        
        /** In addition to the functionality provided by getReducedStructure also provides a way to specify sub-regions of a structure with the following 
         * specification:
         * 
         * 
         * ranges can be surrounded by ( and ). (but will be removed).
         * ranges are specified as
         * PDBresnum1 : PDBresnum2
         * 
         *  a list of ranges is separated by ,
         *  
         *  Example
         *  4GCR (A:1-83)
         *  1CDG (A:407-495,A:582-686)
         *  1CDG (A_407-495,A_582-686)
         * 
         * @param s The full structure
         * @param ranges A comma-seperated list of ranges, optionally surrounded by parentheses
         * @return Substructure of s specified by ranges
         */
        
        @SuppressWarnings("deprecation")
        public static final Structure getSubRanges(Structure s, String ranges ) 
        throws StructureException
        {
                Structure struc = getReducedStructure(s, null);

                if ( ranges == null || ranges.equals(""))
                        throw new IllegalArgumentException("ranges can't be null or empty");

                ranges = ranges.trim();

                if ( ranges.startsWith("("))
                        ranges = ranges.substring(1);
                if ( ranges.endsWith(")")) {
                        ranges = ranges.substring(0,ranges.length()-1);
                }

                //special case: '-' means 'everything'
                if ( ranges.equals("-") ) {
                        return s;
                }

                Structure newS = new StructureImpl();
                
                newS.setHeader(s.getHeader());
                newS.setPDBCode(s.getPDBCode());
                newS.setPDBHeader(s.getPDBHeader());
                newS.setName(s.getName());
                newS.setDBRefs(s.getDBRefs());
                newS.setNmr(s.isNmr());
                newS.setBiologicalAssembly(s.isBiologicalAssembly());
                newS.getPDBHeader().setDescription("sub-range " + ranges + " of "  + newS.getPDBCode() + " " + s.getPDBHeader().getDescription());
                
                // TODO The following should be only copied for atoms which are present in the range.
                //newS.setCompounds(s.getCompounds());
                //newS.setConnections(s.getConnections());
                //newS.setSSBonds(s.getSSBonds());
                //newS.setSites(s.getSites());
                
                String[] rangS =ranges.split(",");

                StringWriter name = new StringWriter();
                name.append(s.getName());
                boolean firstRange = true;
                String prevChainId = null;

                // parse the ranges, adding the specified residues to newS
                for ( String r: rangS){
                        //System.out.println(">"+r+"<");
                        // Match a single range, eg "A_4-27"
                        
                        Matcher matcher = pdbNumRangeRegex.matcher(r);
                        if( ! matcher.matches() ){
                                throw new StructureException("wrong range specification, should be provided as chainID_pdbResnum1-pdbRensum2: "+ranges);
                        }
                        String chainId = matcher.group(1);
                        Chain chain;
                        
                        if(chainId.equals("_") && struc.size() == 1) {
                                // Handle special case of "_" chain for single-chain proteins
                                chain = struc.getChain(0);
                        } else {
                                // Explicit chain
                                chain = struc.getChainByPDB(chainId);
                        }
                        
                        Group[] groups;
                        
                        String pdbresnumStart = matcher.group(2);
                        String pdbresnumEnd   = matcher.group(3);
                        
                        if ( ! firstRange){
                                name.append( ",");
                        } else {
                                name.append(AtomCache.CHAIN_SPLIT_SYMBOL);
                        }
                        if( pdbresnumStart != null && pdbresnumEnd != null) {
                                // not a full chain
                                //since Java doesn't allow '+' before integers, fix this up.
                                if(pdbresnumStart.charAt(0) == '+')
                                        pdbresnumStart = pdbresnumStart.substring(1);
                                if(pdbresnumEnd.charAt(0) == '+')
                                        pdbresnumEnd = pdbresnumEnd.substring(1);
                                groups = chain.getGroupsByPDB(pdbresnumStart, pdbresnumEnd);
                                
                                name.append( chainId + AtomCache.UNDERSCORE + pdbresnumStart+"-" + pdbresnumEnd);
                                
                        } else {
                                // full chain
                                groups = chain.getAtomGroups().toArray(new Group[chain.getAtomGroups().size()]);
                                name.append(chainId);
                        }
                        firstRange = true;
                        
                        // Create new chain, if needed
                        Chain c = null;
                        if ( prevChainId == null) {
                                // first chain...
                                c = new ChainImpl();
                                c.setChainID(chain.getChainID());
                                newS.addChain(c);
                        } else if ( prevChainId.equals(chain.getChainID())) {
                                c = newS.getChainByPDB(prevChainId);

                        } else {
                                try {
                                        c = newS.getChainByPDB(chain.getChainID());
                                } catch (StructureException e){
                                        // chain not in structure yet...
                                        c = new ChainImpl();
                                        c.setChainID(chain.getChainID());
                                        newS.addChain(c);
                                }
                        }

                        // add the groups to the chain:
                        for ( Group g: groups) {
                                c.addGroup(g);
                        }

                        prevChainId = c.getChainID();
                }
                
                newS.setName(name.toString());
                
                return newS;
        }

        public static final String convertAtomsToSeq(Atom[] atoms) {

                StringBuffer buf = new StringBuffer();
                Group prevGroup  = null;
                for (Atom a : atoms){
                        Group g = a.getGroup();
                        if ( prevGroup != null) {
                                if ( prevGroup.equals(g)) {
                                        // we add each group only once.
                                        continue;
                                }
                        }
                        String code3 = g.getPDBName();
                        try {
                                buf.append(convert_3code_1code(code3) );
                        } catch (UnknownPdbAminoAcidException e){
                                buf.append('X');
                        }
                        prevGroup = g;

                }
                return buf.toString();
        }

        /** get a PDB residue number object for this group
         * 
         * @param g Group object
         * @return a ResidueNumber object
         * @deprecated replaced by  Group.getResidueNumber()
         */
        public static final ResidueNumber getPDBResidueNumber(Group g){

                return g.getResidueNumber();

        }

        /** Get a group represented by a ResidueNumber.
         * 
         * @param struc a {@link Structure}
         * @param pdbResNum a {@link ResidueNumber}
         * @return a group in the structure that is represented by the pdbResNum. 
         * @throws StructureException if the group cannot be found.
         */
        public static final Group getGroupByPDBResidueNumber(Structure struc, 
                        ResidueNumber pdbResNum) throws StructureException {
                if (struc == null || pdbResNum==null) {
                        throw new IllegalArgumentException("Null argument(s).");
                }

                Chain chain = struc.findChain(pdbResNum.getChainId());

//              String numIns = "" + pdbResNum.getSeqNum();
//              if (pdbResNum.getInsCode() != null) {
//                      numIns += pdbResNum.getInsCode();
//              }

                
                return chain.getGroupByPDB(pdbResNum);
        }

        /*
         * Returns a List of Groups in a structure within the distance specified of a given group.
         */
        public static List<Group> getGroupsWithinShell(Structure structure, Group group, double distance, boolean includeWater) {
            Set<Group> returnSet = new LinkedHashSet<Group>();

            //square the distance to use as a comparison against getDistanceFast which returns the square of a distance.
            distance = distance * distance;

            for (Atom atomA : group.getAtoms()) {
                for (Chain chain : structure.getChains()) {
                    for (Group chainGroup : chain.getAtomGroups()) {
//                        System.out.println("Checking group: " + chainGroup);
                        if (chainGroup.getResidueNumber().equals(group.getResidueNumber())) {
                            continue;
                        }
                        else if (!includeWater && chainGroup.getPDBName().equals("HOH")) {
                            continue;
                        }
                        else {
                            for (Atom atomB : chainGroup.getAtoms()) {
                            try {
                                //use getDistanceFast as we are doing a lot of comparisons
                                double dist = Calc.getDistanceFast(atomA, atomB);
                                if (dist <= distance) {
                                    returnSet.add(chainGroup);
//                                    System.out.println(String.format("%s within %s of %s", atomB, dist, atomA));
                                    break;
                                }
                                
                            } catch (StructureException ex) {
                                Logger.getLogger(StructureTools.class.getName()).log(Level.SEVERE, null, ex);
                            }
                                
                            }
                        }
                    }
                }
            }
            List<Group> returnList = new ArrayList<Group>();
            returnList.addAll(returnSet);
            return returnList;
        }

        /*
         * Very simple distance-based bond calculator. Will give approximations,
         * but do not rely on this to be chemically correct.
         */
        public static List<Bond> findBonds(Group group, List<Group> groups) {
            List<Bond> bondList = new ArrayList<Bond>();
            for (Atom atomA : group.getAtoms()) {
                for (Group groupB : groups) {
                    if (groupB.getType().equals(GroupType.HETATM)) {
                        continue;
                    }
                    for (Atom atomB : groupB.getAtoms()) {
                        try {
                            double dist = Calc.getDistance(atomA, atomB);
                            BondType bondType = BondType.UNDEFINED;
                            if (dist <= 2) {
                                bondType = BondType.COVALENT;
                                Bond bond = new Bond(dist, bondType, group, atomA, groupB, atomB);
                                bondList.add(bond);
//                                    System.out.println(String.format("%s within %s of %s", atomB, dist, atomA));
                            }
                            else if (dist <= 3.25) {
                                
                                if (isHbondDonorAcceptor(atomA) && isHbondDonorAcceptor(atomB)) {
                                    bondType = BondType.HBOND;
                                }
                                else if (atomA.getElement().isMetal() && isHbondDonorAcceptor(atomB)) {
                                    bondType = BondType.METAL;
                                }
                                else if (atomA.getElement().equals(Element.C) && atomB.getElement().equals(Element.C)) {
                                    bondType = BondType.HYDROPHOBIC;
                                }
                                //not really interested in 'undefined' types
                                if (bondType != BondType.UNDEFINED) {
                                    Bond bond = new Bond(dist, bondType, group, atomA, groupB, atomB);
                                    bondList.add(bond);
                                }
//                                    System.out.println(String.format("%s within %s of %s", atomB, dist, atomA));
                            } else if (dist <= 3.9) {
                                if (atomA.getElement().equals(Element.C) && atomB.getElement().equals(Element.C)) {
                                    bondType = BondType.HYDROPHOBIC;
                                }
                                //not really interested in 'undefined' types
                                if (bondType != BondType.UNDEFINED) {
                                    Bond bond = new Bond(dist, bondType, group, atomA, groupB, atomB);
                                    bondList.add(bond);
                                }
                            }

                        } catch (StructureException ex) {
                            Logger.getLogger(StructureTools.class.getName()).log(Level.SEVERE, null, ex);
                        }

                    }
                }
            }


            return bondList;
        }

        private static boolean isHbondDonorAcceptor(Atom atom) {
            if (hBondDonorAcceptors.contains(atom.getElement())) {
                return true;
            }
            return false;
        }
}