StructureTools.java

/*
 *                  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.nbio.structure;

import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
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 javax.vecmath.Point3d;

import org.biojava.nbio.structure.align.util.AtomCache;
import org.biojava.nbio.structure.contact.AtomContactSet;
import org.biojava.nbio.structure.contact.Grid;
import org.biojava.nbio.structure.io.FileParsingParameters;
import org.biojava.nbio.structure.io.PDBFileParser;
import org.biojava.nbio.core.util.FileDownloadUtils;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

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

	private static final Logger logger = LoggerFactory
			.getLogger(StructureTools.class);

	// Amino Acid backbone
	/**
	 * The atom name of the backbone C-alpha atom. Note that this can be
	 * ambiguous depending on the context since Calcium atoms use the same name
	 * in PDB.
	 */
	public static final String CA_ATOM_NAME = "CA";

	/**
	 * The atom name for the backbone amide nitrogen
	 */
	public static final String N_ATOM_NAME = "N";

	/**
	 * The atom name for the backbone carbonyl
	 */
	public static final String C_ATOM_NAME = "C";

	/**
	 * The atom name for the backbone carbonyl oxygen
	 */
	public static final String O_ATOM_NAME = "O";

	/**
	 * The atom name of the side-chain C-beta atom
	 */
	public static final String CB_ATOM_NAME = "CB";

	// Nucleotide backbone
	/**
	 * The atom name of the backbone C1' in RNA
	 */
	public static final String C1_ATOM_NAME = "C1'";
	/**
	 * The atom name of the backbone C2' in RNA
	 */
	public static final String C2_ATOM_NAME = "C2'";
	/**
	 * The atom name of the backbone C3' in RNA
	 */
	public static final String C3_ATOM_NAME = "C3'";
	/**
	 * The atom name of the backbone C4' in RNA
	 */
	public static final String C4_ATOM_NAME = "C4'";
	/**
	 * The atom name of the backbone O2' in RNA
	 */
	public static final String O2_ATOM_NAME = "O2'";
	/**
	 * The atom name of the backbone O3' in RNA
	 */
	public static final String O3_ATOM_NAME = "O3'";
	/**
	 * The atom name of the backbone O4' in RNA
	 */
	public static final String O4_ATOM_NAME = "O4'";
	/**
	 * The atom name of the backbone O4' in RNA
	 */
	public static final String O5_ATOM_NAME = "O5'";
	/**
	 * The atom name of the backbone O4' in RNA
	 */
	public static final String OP1_ATOM_NAME = "OP1";
	/**
	 * The atom name of the backbone O4' in RNA
	 */
	public static final String OP2_ATOM_NAME = "OP2";
	/**
	 * The atom name of the backbone phosphate in RNA
	 */
	public static final String P_ATOM_NAME = "P";

	/**
	 * The atom used as representative for nucleotides, equivalent to
	 * {@link #CA_ATOM_NAME} for proteins
	 */
	public static final String NUCLEOTIDE_REPRESENTATIVE = P_ATOM_NAME;

	/**
	 * The character to use for unknown compounds in sequence strings
	 */
	public static final char UNKNOWN_GROUP_LABEL = 'X';

	/**
	 * Below this ratio of aminoacid/nucleotide residues to the sequence total,
	 * we use simple majority of aminoacid/nucleotide residues to decide the
	 * character of the chain (protein/nucleotide)
	 */
	public static final double RATIO_RESIDUES_TO_TOTAL = 0.95;

	/**
	 * Threshold for plausible binding of a ligand to the selected substructure
	 */
	public static final double DEFAULT_LIGAND_PROXIMITY_CUTOFF = 5;

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

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

	private static final Set<Element> hBondDonorAcceptors;

	static {
		nucleotides30 = new HashMap<String, Character>();
		nucleotides30.put("DA", 'A');
		nucleotides30.put("DC", 'C');
		nucleotides30.put("DG", 'G');
		nucleotides30.put("DT", 'T');
		nucleotides30.put("DI", 'I');
		nucleotides30.put("A", 'A');
		nucleotides30.put("G", 'G');
		nucleotides30.put("C", 'C');
		nucleotides30.put("U", 'U');
		nucleotides30.put("I", 'I');

		// the DNA linkers - the +C , +G, +A +T +U and +I have been replaced
		// with these:
		nucleotides30.put("TAF", UNKNOWN_GROUP_LABEL); // Fluorinated Thymine
		nucleotides30.put("TC1", UNKNOWN_GROUP_LABEL); // Furanosyl
		nucleotides30.put("TFE", UNKNOWN_GROUP_LABEL); // Fluorinated Thymine
		nucleotides30.put("TFO", UNKNOWN_GROUP_LABEL); // Tenofovir (3'
		// terminator)
		nucleotides30.put("TGP", UNKNOWN_GROUP_LABEL); // Guanine variant
		nucleotides30.put("THX", UNKNOWN_GROUP_LABEL); // 5' terminator
		nucleotides30.put("TLC", UNKNOWN_GROUP_LABEL); // Thymine with dicyclic
		// sugar
		nucleotides30.put("TLN", UNKNOWN_GROUP_LABEL); // locked Thymine
		nucleotides30.put("LCG", UNKNOWN_GROUP_LABEL); // locked Guanine
		nucleotides30.put("TP1", UNKNOWN_GROUP_LABEL); // Thymine peptide
		// nucleic acid, with
		// added methyl
		nucleotides30.put("CP1", UNKNOWN_GROUP_LABEL); // Cytidine peptide
		// nucleic acid, with
		// added methyl
		nucleotides30.put("TPN", UNKNOWN_GROUP_LABEL); // Thymine peptide
		// nucleic acid
		nucleotides30.put("CPN", UNKNOWN_GROUP_LABEL); // Cytidine peptide
		// nucleic acid
		nucleotides30.put("GPN", UNKNOWN_GROUP_LABEL); // Guanine peptide
		// nucleic acid
		nucleotides30.put("APN", UNKNOWN_GROUP_LABEL); // Adenosine peptide
		// nucleic acid
		nucleotides30.put("TPC", UNKNOWN_GROUP_LABEL); // Thymine variant

		// 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, Character>();
		String[] names = { "C", "G", "A", "T", "U", "I", "+C", "+G", "+A",
				"+T", "+U", "+I" };
		for (String n : names) {
			nucleotides23.put(n, n.charAt(n.length() - 1));
		}

		aminoAcids = new HashMap<String, Character>();
		aminoAcids.put("GLY", 'G');
		aminoAcids.put("ALA", 'A');
		aminoAcids.put("VAL", 'V');
		aminoAcids.put("LEU", 'L');
		aminoAcids.put("ILE", 'I');
		aminoAcids.put("PHE", 'F');
		aminoAcids.put("TYR", 'Y');
		aminoAcids.put("TRP", 'W');
		aminoAcids.put("PRO", 'P');
		aminoAcids.put("HIS", 'H');
		aminoAcids.put("LYS", 'K');
		aminoAcids.put("ARG", 'R');
		aminoAcids.put("SER", 'S');
		aminoAcids.put("THR", 'T');
		aminoAcids.put("GLU", 'E');
		aminoAcids.put("GLN", 'Q');
		aminoAcids.put("ASP", 'D');
		aminoAcids.put("ASN", 'N');
		aminoAcids.put("CYS", 'C');
		aminoAcids.put("MET", 'M');
		// MSE is only found as a molecular replacement for MET
		aminoAcids.put("MSE", '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", 'U');
		aminoAcids.put("SEC", 'U');
		aminoAcids.put("PYH", 'O');
		aminoAcids.put("PYL", 'O');

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

	}

	/**
	 * Count how many 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 = 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);
		for (Chain c : chains) {
			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. If the group does not contain all requested atoms then no atoms
	 * are added for that 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>();

		extractAtoms(atomNames, chains, atoms);
		return 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. If the group does not
	 * contain all requested atoms then no atoms are added for that 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);
			extractAtoms(atomNames, chains, atoms);
		}
		return atoms.toArray(new Atom[atoms.size()]);

	}

	/**
	 * Convert all atoms of the structure (all models) into an Atom array
	 *
	 * @param s
	 *            input structure
	 * @return all atom array
	 */
	public static final Atom[] getAllAtomArray(Structure s) {
		List<Atom> atoms = new ArrayList<Atom>();

		AtomIterator iter = new AtomIterator(s);
		while (iter.hasNext()) {
			Atom a = iter.next();
			atoms.add(a);
		}
		return atoms.toArray(new Atom[atoms.size()]);
	}
	/**
	 * Convert all atoms of the structure (specified model) into an Atom array
	 *
	 * @param s
	 *            input structure
	 * @return all atom array
	 */
	public static final Atom[] getAllAtomArray(Structure s, int model) {
		List<Atom> atoms = new ArrayList<Atom>();

		AtomIterator iter = new AtomIterator(s,model);
		while (iter.hasNext()) {
			Atom a = iter.next();
			atoms.add(a);
		}
		return atoms.toArray(new Atom[atoms.size()]);

	}

	/**
	 * Returns and array of all atoms of the chain, including
	 * Hydrogens (if present) and all HETATOMs. Waters are not included.
	 *
	 * @param c
	 *            input chain
	 * @return all atom array
	 */
	public static final Atom[] getAllAtomArray(Chain c) {
		List<Atom> atoms = new ArrayList<Atom>();

		for (Group g : c.getAtomGroups()) {
			if (g.isWater())
				continue;
			for (Atom a : g.getAtoms()) {
				atoms.add(a);
			}
		}
		return atoms.toArray(new Atom[atoms.size()]);
	}

	/**
	 * List of groups from the structure not included in ca (e.g. ligands).
	 *
	 * Unaligned groups are searched from all chains referenced in ca, as well
	 * as any chains in the first model of the structure from ca[0], if any.
	 *
	 * @param ca an array of atoms
	 * @return
	 */
	public static List<Group> getUnalignedGroups(Atom[] ca) {
		Set<Chain> chains = new HashSet<Chain>();
		Set<Group> caGroups = new HashSet<Group>();

		// Create list of all chains in this structure
		Structure s = null;
		if (ca.length > 0) {
			Group g = ca[0].getGroup();
			if (g != null) {
				Chain c = g.getChain();
				if (c != null) {
					s = c.getStructure();
				}
			}
		}
		if (s != null) {
			// Add all chains from the structure
			for (Chain c : s.getChains(0)) {
				chains.add(c);
			}
		}

		// Add groups and chains from ca
		for (Atom a : ca) {
			Group g = a.getGroup();
			if (g != null) {
				caGroups.add(g);

				Chain c = g.getChain();
				if (c != null) {
					chains.add(c);
				}
			}
		}

		// Iterate through all chains, finding groups not in ca
		List<Group> unadded = new ArrayList<Group>();
		for (Chain c : chains) {
			for (Group g : c.getAtomGroups()) {
				if (!caGroups.contains(g)) {
					unadded.add(g);
				}
			}
		}
		return unadded;
	}

	/**
	 * Finds all ligand groups from the target which fall within the cutoff distance
	 * of some atom from the query set.
	 *
	 * @param target Set of groups including the ligands
	 * @param query Atom selection
	 * @param cutoff Distance from query atoms to consider, in angstroms
	 * @return All groups from the target with at least one atom within cutoff of a query atom
	 * @see StructureTools#DEFAULT_LIGAND_PROXIMITY_CUTOFF
	 */
	public static List<Group> getLigandsByProximity(Collection<Group> target, Atom[] query, double cutoff) {
		// Spatial hashing of the reduced structure
		Grid grid = new Grid(cutoff);
		grid.addAtoms(query);

		List<Group> ligands = new ArrayList<>();
		for(Group g :target ) {
			// don't worry about waters
			if(g.isWater()) {
				continue;
			}

			if(g.isPolymeric() ) {
				// Polymers aren't ligands
				continue;
			}

			// It is a ligand!

			// Check that it's within cutoff of something in reduced
			List<Atom> groupAtoms = g.getAtoms();
			if( ! grid.hasAnyContact(Calc.atomsToPoints(groupAtoms))) {
				continue;
			}

			ligands.add(g);
		}
		return ligands;
	}

	/**
	 * Adds a particular group to a structure. A new chain will be created if necessary.
	 *
	 * <p>When adding multiple groups, pass the return value of one call as the
	 * chainGuess parameter of the next call for efficiency.
	 * <pre>
	 * Chain guess = null;
	 * for(Group g : groups) {
	 *     guess = addGroupToStructure(s, g, guess );
	 * }
	 * </pre>
	 * @param s structure to receive the group
	 * @param g group to add
	 * @param chainGuess (optional) If not null, should be a chain from s. Used
	 *  to improve performance when adding many groups from the same chain
	 * @param clone Indicates whether the input group should be cloned before
	 *  being added to the new chain
	 * @return the chain g was added to
	 */
	public static Chain addGroupToStructure(Structure s, Group g, int model, Chain chainGuess, boolean clone ) {
		synchronized(s) {
			// Find or create the chain
			String chainId = g.getChainId();
			assert !chainId.isEmpty();
			Chain chain;
			if(chainGuess != null && chainGuess.getId() == chainId) {
				// previously guessed chain
				chain = chainGuess;
			} else {
				// Try to guess
				chain = s.getChain(chainId, model);
				if(chain == null) {
					// no chain found
					chain = new ChainImpl();
					chain.setId(chainId);

					Chain oldChain = g.getChain();
					chain.setName(oldChain.getName());

					EntityInfo oldEntityInfo = oldChain.getEntityInfo();

					EntityInfo newEntityInfo;
					if(oldEntityInfo == null) {
						newEntityInfo = new EntityInfo();
						s.addEntityInfo(newEntityInfo);
					} else {
						newEntityInfo = s.getEntityById(oldEntityInfo.getMolId());
						if( newEntityInfo == null ) {
							newEntityInfo = new EntityInfo(oldEntityInfo);
							s.addEntityInfo(newEntityInfo);
						}
					}
					newEntityInfo.addChain(chain);
					chain.setEntityInfo(newEntityInfo);

					// TODO Do the seqres need to be cloned too? -SB 2016-10-7
					chain.setSeqResGroups(oldChain.getSeqResGroups());
					chain.setSeqMisMatches(oldChain.getSeqMisMatches());

					s.addChain(chain,model);
				}
			}

			// Add cloned group
			if(clone) {
				g = (Group)g.clone();
			}
			chain.addGroup(g);

			return chain;
		}
	}

	/**
	 * Add a list of groups to a new structure. Chains will be automatically
	 * created in the new structure as needed.
	 * @param s structure to receive the group
	 * @param g group to add
	 * @param clone Indicates whether the input groups should be cloned before
	 *  being added to the new chain
	 */
	public static void addGroupsToStructure(Structure s, Collection<Group> groups, int model, boolean clone) {
		Chain chainGuess = null;
		for(Group g : groups) {
			chainGuess = addGroupToStructure(s, g, model, chainGuess, clone);
		}
	}

	/**
	 * Expand a set of atoms into all groups from the same structure.
	 *
	 * If the structure is set, only the first atom is used (assuming all
	 * atoms come from the same original structure).
	 * If the atoms aren't linked to a structure (for instance, for cloned atoms),
	 * searches all chains of all atoms for groups.
	 * @param atoms Sample of atoms
	 * @return All groups from all chains accessible from the input atoms
	 */
	public static Set<Group> getAllGroupsFromSubset(Atom[] atoms) {
		return getAllGroupsFromSubset(atoms,null);
	}
	/**
	 * Expand a set of atoms into all groups from the same structure.
	 *
	 * If the structure is set, only the first atom is used (assuming all
	 * atoms come from the same original structure).
	 * If the atoms aren't linked to a structure (for instance, for cloned atoms),
	 * searches all chains of all atoms for groups.
	 * @param atoms Sample of atoms
	 * @param types Type of groups to return (useful for getting only ligands, for instance).
	 *  Null gets all groups.
	 * @return All groups from all chains accessible from the input atoms
	 */
	public static Set<Group> getAllGroupsFromSubset(Atom[] atoms,GroupType types) {
		// Get the full structure
		Structure s = null;
		if (atoms.length > 0) {
			Group g = atoms[0].getGroup();
			if (g != null) {
				Chain c = g.getChain();
				if (c != null) {
					s = c.getStructure();
				}
			}
		}
		// Collect all groups from the structure
		Set<Chain> allChains = new HashSet<>();
		if( s != null ) {
			allChains.addAll(s.getChains());
		}
		// In case the structure wasn't set, need to use ca chains too
		for(Atom a : atoms) {
			Group g = a.getGroup();
			if(g != null) {
				Chain c = g.getChain();
				if( c != null ) {
					allChains.add(c);
				}
			}
		}

		if(allChains.isEmpty() ) {
			return Collections.emptySet();
		}

		// Extract all ligand groups
		Set<Group> full = new HashSet<>();
		for(Chain c : allChains) {
			if(types == null) {
				full.addAll(c.getAtomGroups());
			} else {
				full.addAll(c.getAtomGroups(types));
			}
		}

		return full;
	}


	/**
	 * Returns and array of all non-Hydrogen atoms in the given Structure,
	 * optionally including HET atoms or not. Waters are not included.
	 *
	 * @param s
	 * @param hetAtoms
	 *            if true HET atoms are included in array, if false they are not
	 * @return
	 */
	public static final Atom[] getAllNonHAtomArray(Structure s, boolean hetAtoms) {
		AtomIterator iter = new AtomIterator(s);
		return getAllNonHAtomArray(s, hetAtoms, iter);
	}
	/**
	 * Returns and array of all non-Hydrogen atoms in the given Structure,
	 * optionally including HET atoms or not. Waters are not included.
	 *
	 * @param s
	 * @param hetAtoms
	 *            if true HET atoms are included in array, if false they are not
	 * @param modelNr Model number to draw atoms from
	 * @return
	 */
	public static final Atom[] getAllNonHAtomArray(Structure s, boolean hetAtoms, int modelNr) {
		AtomIterator iter = new AtomIterator(s,modelNr);
		return getAllNonHAtomArray(s, hetAtoms, iter);
	}
	private static final Atom[] getAllNonHAtomArray(Structure s, boolean hetAtoms, AtomIterator iter) {
		List<Atom> atoms = new ArrayList<Atom>();

		while (iter.hasNext()) {
			Atom a = iter.next();
			if (a.getElement() == Element.H)
				continue;

			Group g = a.getGroup();

			if (g.isWater())
				continue;

			if (!hetAtoms && g.getType().equals(GroupType.HETATM))
				continue;

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

	/**
	 * Returns and array of all non-Hydrogen atoms in the given Chain,
	 * optionally including HET atoms or not Waters are not included.
	 *
	 * @param c
	 * @param hetAtoms
	 *            if true HET atoms are included in array, if false they are not
	 * @return
	 */
	public static final Atom[] getAllNonHAtomArray(Chain c, boolean hetAtoms) {
		List<Atom> atoms = new ArrayList<Atom>();

		for (Group g : c.getAtomGroups()) {
			if (g.isWater())
				continue;
			for (Atom a : g.getAtoms()) {

				if (a.getElement() == Element.H)
					continue;

				if (!hetAtoms && g.getType().equals(GroupType.HETATM))
					continue;

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

	/**
	 * Returns and array of all non-Hydrogen atoms coordinates in the given Chain,
	 * optionally including HET atoms or not Waters are not included.
	 *
	 * @param c
	 * @param hetAtoms
	 *            if true HET atoms are included in array, if false they are not
	 * @return
	 */
	public static final Point3d[] getAllNonHCoordsArray(Chain c, boolean hetAtoms) {
		List<Point3d> atoms = new ArrayList<Point3d>();

		for (Group g : c.getAtomGroups()) {
			if (g.isWater())
				continue;
			for (Atom a : g.getAtoms()) {

				if (a.getElement() == Element.H)
					continue;

				if (!hetAtoms && g.getType().equals(GroupType.HETATM))
					continue;

				atoms.add(a.getCoordsAsPoint3d());
			}
		}
		return atoms.toArray(new Point3d[atoms.size()]);
	}

	/**
	 * Adds to the given atoms list, all atoms of groups that contained all
	 * requested atomNames, i.e. if a group does not contain all of the
	 * requested atom names, its atoms won't be added.
	 *
	 * @param atomNames
	 * @param chains
	 * @param atoms
	 */
	private static void extractAtoms(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 (String atomName : atomNames) {
					Atom a = g.getAtom(atomName);

					if (a == null) {
						// this group does not have a required atom, skip it...
						thisGroupAllAtoms = false;
						break;
					}
					thisGroupAtoms.add(a);
				}
				if (thisGroupAllAtoms) {
					// add the atoms of this group to the array.
					for (Atom a : thisGroupAtoms) {
						atoms.add(a);
					}
				}

			}
		}
	}

	/**
	 * Returns an array of the requested Atoms from the Chain 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. If the group does not
	 * contain all requested atoms then no atoms are added for that 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<Atom> atoms = new ArrayList<Atom>();

		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 (String atomName : atomNames) {
				Atom a = g.getAtom(atomName);
				if (a == null) {
					logger.debug("Group " + g.getResidueNumber() + " ("
							+ g.getPDBName()
							+ ") does not have the required atom '" + atomName
							+ "'");
					// this group does not have a required atom, skip it...
					thisGroupAllAtoms = false;
					break;
				}
				thisGroupAtoms.add(a);
			}

			if (thisGroupAllAtoms) {
				// add the atoms of this group to the array.
				for (Atom a : thisGroupAtoms) {
					atoms.add(a);
				}
			}

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

	}

	/**
	 * Returns an Atom array of the C-alpha atoms. Any atom that is a carbon and
	 * has CA name will be returned.
	 *
	 * @param c
	 *            the structure object
	 * @return an Atom[] array
	 * @see #getRepresentativeAtomArray(Chain)
	 */
	public static final Atom[] getAtomCAArray(Chain c) {
		List<Atom> atoms = new ArrayList<Atom>();

		for (Group g : c.getAtomGroups()) {
			if (g.hasAtom(CA_ATOM_NAME)
					&& g.getAtom(CA_ATOM_NAME).getElement() == Element.C) {
				atoms.add(g.getAtom(CA_ATOM_NAME));
			}
		}

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

	/**
	 * Gets a representative atom for each group that is part of the chain
	 * backbone. Note that modified aminoacids won't be returned as part of the
	 * backbone if the {@link org.biojava.nbio.structure.io.mmcif.ReducedChemCompProvider} was used to load the
	 * structure.
	 *
	 * For amino acids, the representative is a CA carbon. For nucleotides, the
	 * representative is the {@value #NUCLEOTIDE_REPRESENTATIVE}. Other group
	 * types will be ignored.
	 *
	 * @param c
	 * @return representative Atoms of the chain backbone
	 * @since Biojava 4.1.0
	 */
	public static final Atom[] getRepresentativeAtomArray(Chain c) {
		List<Atom> atoms = new ArrayList<Atom>();

		for (Group g : c.getAtomGroups()) {

			switch (g.getType()) {
			case AMINOACID:
				if (g.hasAtom(CA_ATOM_NAME)
						&& g.getAtom(CA_ATOM_NAME).getElement() == Element.C) {
					atoms.add(g.getAtom(CA_ATOM_NAME));
				}
				break;
			case NUCLEOTIDE:
				if (g.hasAtom(NUCLEOTIDE_REPRESENTATIVE)) {
					atoms.add(g.getAtom(NUCLEOTIDE_REPRESENTATIVE));
				}
				break;
			default:
				// don't add
			}
		}

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

	}

	/**
	 * 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 input array.
	 *
	 * @param ca
	 *            array of representative atoms, e.g. CA atoms
	 * @return Atom array
	 * @since Biojava 4.1.0
	 */
	public static final Atom[] cloneAtomArray(Atom[] ca) {
		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.getName().equals(parentC.getName())) {
					newChain = c;
					break;
				}
			}
			if (newChain == null) {
				newChain = new ChainImpl();
				newChain.setId(parentC.getId());
				newChain.setName(parentC.getName());
				model.add(newChain);
			}

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

			newCA[apos] = parentN.getAtom(a.getName());
			try {
				// if the group doesn't exist yet, this produces a StructureException
				newChain.getGroupByPDB(parentN.getResidueNumber());
			} catch (StructureException e) {
				// the group doesn't exist yet in the newChain, let's add it
				newChain.addGroup(parentN);
			}

		}
		return newCA;
	}

	/**
	 * Clone a set of representative 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.getName().equals(parentC.getName())) {
					newChain = c;
					break;
				}
			}
			if (newChain == null) {
				newChain = new ChainImpl();
				newChain.setName(parentC.getName());
				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
	 */
	public static Atom[] duplicateCA2(Atom[] ca2) {
		// 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 single atom

			if (c == null) {
				c = new ChainImpl();
				Chain orig = a.getGroup().getChain();
				c.setId(orig.getId());
				c.setName(orig.getName());
			} else {
				Chain orig = a.getGroup().getChain();
				if (!orig.getId().equals(prevChainId)) {
					c = new ChainImpl();
					c.setId(orig.getId());
					c.setName(orig.getName());
				}
			}

			c.addGroup(g);
			ca2clone[pos] = g.getAtom(a.getName());

			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.setId(orig.getId());
				c.setName(orig.getName());
			} else {
				Chain orig = a.getGroup().getChain();
				if (!orig.getId().equals(prevChainId)) {
					c = new ChainImpl();
					c.setId(orig.getId());
					c.setName(orig.getName());
				}
			}

			c.addGroup(g);
			ca2clone[pos] = g.getAtom(a.getName());

			pos++;
		}

		return ca2clone;

	}

	/**
	 * Return an Atom array of the C-alpha atoms. Any atom that is a carbon and
	 * has CA name will be returned.
	 *
	 * @param s
	 *            the structure object
	 * @return an Atom[] array
	 * @see #getRepresentativeAtomArray(Structure)
	 */
	public static Atom[] getAtomCAArray(Structure s) {

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

		for (Chain c : s.getChains()) {
			for (Group g : c.getAtomGroups()) {
				if (g.hasAtom(CA_ATOM_NAME)
						&& g.getAtom(CA_ATOM_NAME).getElement() == Element.C) {
					atoms.add(g.getAtom(CA_ATOM_NAME));
				}
			}
		}

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

	/**
	 * Gets a representative atom for each group that is part of the chain
	 * backbone. Note that modified aminoacids won't be returned as part of the
	 * backbone if the {@link org.biojava.nbio.structure.io.mmcif.ReducedChemCompProvider} was used to load the
	 * structure.
	 *
	 * For amino acids, the representative is a CA carbon. For nucleotides, the
	 * representative is the {@value #NUCLEOTIDE_REPRESENTATIVE}. Other group
	 * types will be ignored.
	 *
	 * @param s
	 *            Input structure
	 * @return representative Atoms of the structure backbone
	 * @since Biojava 4.1.0
	 */
	public static Atom[] getRepresentativeAtomArray(Structure s) {

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

		for (Chain c : s.getChains()) {
			Atom[] chainAtoms = getRepresentativeAtomArray(c);
			for (Atom a : chainAtoms) {
				atoms.add(a);
			}
		}

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

	/**
	 * Return an Atom array of the main chain atoms: CA, C, N, O Any group that
	 * contains those atoms will be included, be it a standard aminoacid or not
	 *
	 * @param s
	 *            the structure object
	 * @return an Atom[] array
	 */
	public static Atom[] getBackboneAtomArray(Structure s) {

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

		for (Chain c : s.getChains()) {
			for (Group g : c.getAtomGroups()) {
				if (g.hasAminoAtoms()) {
					// this means we will only take atoms grom groups that have
					// complete backbones
					for (Atom a : g.getAtoms()) {
						switch (g.getType()) {
						case NUCLEOTIDE:
							// Nucleotide backbone
							if (a.getName().equals(C1_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(C2_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(C3_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(C4_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(O2_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(O3_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(O4_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(O5_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(OP1_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(OP2_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(P_ATOM_NAME))
								atoms.add(a);
							// TODO Allow C4* names as well as C4'? -SB 3/2015
							break;
						case AMINOACID:
						default:
							// we do it this way instead of with g.getAtom() to
							// be sure we always use the same order as original
							if (a.getName().equals(CA_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(C_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(N_ATOM_NAME))
								atoms.add(a);
							if (a.getName().equals(O_ATOM_NAME))
								atoms.add(a);
							break;
						}
					}
				}
			}

		}

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

	/**
	 * Convert three character amino acid codes into single character e.g.
	 * convert CYS to C. Valid 3-letter codes will be those of the standard 20
	 * amino acids plus MSE, CSE, SEC, PYH, PYL (see the {@link #aminoAcids}
	 * map)
	 *
	 * @return the 1 letter code, or null if the given 3 letter code does not
	 *         correspond to an amino acid code
	 * @param groupCode3
	 *            a three character amino acid representation String
	 * @see {@link #get1LetterCode(String)}
	 */
	public static final Character get1LetterCodeAmino(String groupCode3) {
		return aminoAcids.get(groupCode3);
	}

	/**
	 * Convert a three letter amino acid or nucleotide code into a single
	 * character code. If the code does not correspond to an amino acid or
	 * nucleotide, returns {@link #UNKNOWN_GROUP_LABEL}.
	 *
	 * Returned null for nucleotides prior to version 4.0.1.
	 *
	 * @param groupCode3
	 *            three letter representation
	 * @return The 1-letter abbreviation
	 */
	public static final Character get1LetterCode(String groupCode3) {

		Character code1;

		// is it a standard amino acid ?
		code1 = get1LetterCodeAmino(groupCode3);

		if (code1 == null) {
			// hm groupCode3 is not standard
			// perhaps it is a nucleotide?
			groupCode3 = groupCode3.trim();
			if (isNucleotide(groupCode3)) {
				code1 = nucleotides30.get(groupCode3);
				if (code1 == null) {
					code1 = nucleotides23.get(groupCode3);
				}
				if (code1 == null) {
					code1 = UNKNOWN_GROUP_LABEL;
				}
			} else {
				// does not seem to be so let's assume it is
				// nonstandard aminoacid and label it "X"
				// logger.warning("unknown group name "+groupCode3 );
				code1 = UNKNOWN_GROUP_LABEL;
			}
		}

		return code1;

	}

	/**
	 * Test if the three-letter code of an ATOM entry corresponds to a
	 * nucleotide or to an aminoacid.
	 *
	 * @param groupCode3
	 *            3-character code for a group.
	 *
	 */
	public static boolean isNucleotide(String groupCode3) {
		String code = groupCode3.trim();
		return nucleotides30.containsKey(code)
				|| nucleotides23.containsKey(code);
	}

	/**
	 * Reduce a structure to provide a smaller representation . Only takes the
	 * first model of the structure. If chainName 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
	 * @deprecated Use {@link StructureIdentifier#reduce(Structure)} instead (v. 4.2.0)
	 */
	@Deprecated
	public static Structure getReducedStructure(Structure s,
			String chainId) throws StructureException {
		// since we deal here with structure alignments,
		// only use Model 1...

		Structure newS = new StructureImpl();
		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.setBiologicalAssembly(s.isBiologicalAssembly());
		newS.setEntityInfos(s.getEntityInfos());
		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) {
			logger.warn(e.getMessage() + ". Chain id " + chainId
					+ " did not match, trying upper case Chain id.");
			c = s.getChainByPDB(chainId.toUpperCase());

		}
		if (c != null) {
			newS.addChain(c);
			for (EntityInfo comp : s.getEntityInfos()) {
				if (comp.getChainIds() != null
						&& comp.getChainIds().contains(c.getChainID())) {
					// found matching entity info. set description...
					newS.getPDBHeader().setDescription(
							"Chain " + c.getChainID() + " of " + s.getPDBCode()
							+ " " + comp.getDescription());
				}
			}
		}

		return newS;
	}

	public static String convertAtomsToSeq(Atom[] atoms) {

		StringBuilder buf = new StringBuilder();
		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();
			Character code1 = get1LetterCodeAmino(code3);
			if (code1 == null)
				code1 = UNKNOWN_GROUP_LABEL;

			buf.append(code1);

			prevGroup = g;

		}
		return buf.toString();
	}

	/**
	 * 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 Group getGroupByPDBResidueNumber(Structure struc,
			ResidueNumber pdbResNum) throws StructureException {
		if (struc == null || pdbResNum == null) {
			throw new IllegalArgumentException("Null argument(s).");
		}

		Chain chain = struc.getPolyChainByPDB(pdbResNum.getChainName());

		return chain.getGroupByPDB(pdbResNum);
	}

	/**
	 * Returns the set of intra-chain contacts for the given chain for given
	 * atom names, i.e. the contact map. Uses a spatial hashing algorithm that
	 * speeds up the calculation without need of full distance matrix. The
	 * parsing mode {@link FileParsingParameters#setAlignSeqRes(boolean)} needs
	 * to be set to true for this to work.
	 *
	 * @param chain
	 * @param atomNames
	 *            the array with atom names to be used. Beware: CA will do both
	 *            C-alphas an Calciums if null all non-H atoms of non-hetatoms
	 *            will be used
	 * @param cutoff
	 * @return
	 */
	public static AtomContactSet getAtomsInContact(Chain chain,
			String[] atomNames, double cutoff) {
		Grid grid = new Grid(cutoff);

		Atom[] atoms = null;
		if (atomNames == null) {
			atoms = getAllNonHAtomArray(chain, false);
		} else {
			atoms = getAtomArray(chain, atomNames);
		}
		// If tha
		if(atoms.length==0){
			logger.warn("No atoms found for buidling grid!");
			return new AtomContactSet(cutoff);
		}
		grid.addAtoms(atoms);

		return grid.getAtomContacts();
	}

	/**
	 * Returns the set of intra-chain contacts for the given chain for all non-H
	 * atoms of non-hetatoms, i.e. the contact map. Uses a spatial hashing
	 * algorithm that speeds up the calculation without need of full distance
	 * matrix. The parsing mode
	 * {@link FileParsingParameters#setAlignSeqRes(boolean)} needs to be set to
	 * true for this to work.
	 *
	 * @param chain
	 * @param cutoff
	 * @return
	 */
	public static AtomContactSet getAtomsInContact(Chain chain, double cutoff) {
		return getAtomsInContact(chain, (String[]) null, cutoff);
	}

	/**
	 * Returns the set of intra-chain contacts for the given chain for C-alpha
	 * atoms (including non-standard aminoacids appearing as HETATM groups),
	 * i.e. the contact map. Uses a spatial hashing algorithm that speeds up
	 * the calculation without need of full distance matrix. The parsing mode
	 * {@link FileParsingParameters#setAlignSeqRes(boolean)} needs to be set to
	 * true for this to work.
	 *
	 * @param chain
	 * @param cutoff
	 * @return
	 * @see #getRepresentativeAtomsInContact(Chain, double)
	 */
	public static AtomContactSet getAtomsCAInContact(Chain chain, double cutoff) {
		Grid grid = new Grid(cutoff);

		Atom[] atoms = getAtomCAArray(chain);

		grid.addAtoms(atoms);

		return grid.getAtomContacts();
	}

	/**
	 * Returns the set of intra-chain contacts for the given chain for C-alpha
	 * or C3' atoms (including non-standard aminoacids appearing as HETATM
	 * groups), i.e. the contact map. Uses a spatial hashing algorithm that
	 * speeds up the calculation without need of full distance matrix.
	 *
	 * @param chain
	 * @param cutoff
	 * @return
	 * @since Biojava 4.1.0
	 */
	public static AtomContactSet getRepresentativeAtomsInContact(Chain chain,
			double cutoff) {
		Grid grid = new Grid(cutoff);

		Atom[] atoms = getRepresentativeAtomArray(chain);

		grid.addAtoms(atoms);

		return grid.getAtomContacts();
	}

	/**
	 * Returns the set of inter-chain contacts between the two given chains for
	 * the given atom names. Uses a spatial hashing algorithm that speeds up
	 * the calculation without need of full distance matrix. The parsing mode
	 * {@link FileParsingParameters#setAlignSeqRes(boolean)} needs to be set to
	 * true for this to work.
	 *
	 * @param chain1
	 * @param chain2
	 * @param atomNames
	 *            the array with atom names to be used. For Calphas use {"CA"},
	 *            if null all non-H atoms will be used. Note HET atoms are
	 *            ignored unless this parameter is null.
	 * @param cutoff
	 * @param hetAtoms
	 *            if true HET atoms are included, if false they are not
	 * @return
	 */
	public static AtomContactSet getAtomsInContact(Chain chain1, Chain chain2,
			String[] atomNames, double cutoff, boolean hetAtoms) {
		Grid grid = new Grid(cutoff);
		Atom[] atoms1 = null;
		Atom[] atoms2 = null;
		if (atomNames == null) {
			atoms1 = getAllNonHAtomArray(chain1, hetAtoms);
			atoms2 = getAllNonHAtomArray(chain2, hetAtoms);
		} else {
			atoms1 = getAtomArray(chain1, atomNames);
			atoms2 = getAtomArray(chain2, atomNames);
		}
		grid.addAtoms(atoms1, atoms2);

		return grid.getAtomContacts();
	}

	/**
	 * Returns the set of inter-chain contacts between the two given chains for
	 * all non-H atoms. Uses a spatial hashing algorithm that speeds up the
	 * calculation without need of full distance matrix. The parsing mode
	 * {@link FileParsingParameters#setAlignSeqRes(boolean)} needs to be set to
	 * true for this to work.
	 *
	 * @param chain1
	 * @param chain2
	 * @param cutoff
	 * @param hetAtoms
	 *            if true HET atoms are included, if false they are not
	 * @return
	 */
	public static AtomContactSet getAtomsInContact(Chain chain1, Chain chain2,
			double cutoff, boolean hetAtoms) {
		return getAtomsInContact(chain1, chain2, null, cutoff, hetAtoms);
	}

	/**
	 * Finds Groups in {@code structure} that contain at least one Atom that is
	 * within {@code radius} Angstroms of {@code centroid}.
	 *
	 * @param structure
	 *            The structure from which to find Groups
	 * @param centroid
	 *            The centroid of the shell
	 * @param excludeResidues
	 *            A list of ResidueNumbers to exclude
	 * @param radius
	 *            The radius from {@code centroid}, in Angstroms
	 * @param includeWater
	 *            Whether to include Groups whose <em>only</em> atoms are water
	 * @param useAverageDistance
	 *            When set to true, distances are the arithmetic mean (1-norm)
	 *            of the distances of atoms that belong to the group and that
	 *            are within the shell; otherwise, distances are the minimum of
	 *            these values
	 * @return A map of Groups within (or partially within) the shell, to their
	 *         distances in Angstroms
	 */
	public static Map<Group, Double> getGroupDistancesWithinShell(
			Structure structure, Atom centroid,
			Set<ResidueNumber> excludeResidues, double radius,
			boolean includeWater, boolean useAverageDistance) {

		// for speed, we avoid calculating square roots
		radius = radius * radius;

		Map<Group, Double> distances = new HashMap<Group, Double>();

		// we only need this if we're averaging distances
		// note that we can't use group.getAtoms().size() because some the
		// group's atoms be outside the shell
		Map<Group, Integer> atomCounts = new HashMap<Group, Integer>();

		for (Chain chain : structure.getChains()) {
			groupLoop: for (Group chainGroup : chain.getAtomGroups()) {

				// exclude water
				if (!includeWater && chainGroup.isWater())
					continue;

				// check blacklist of residue numbers
				for (ResidueNumber rn : excludeResidues) {
					if (rn.equals(chainGroup.getResidueNumber()))
						continue groupLoop;
				}

				for (Atom testAtom : chainGroup.getAtoms()) {

					// use getDistanceFast as we are doing a lot of comparisons
					double dist = Calc.getDistanceFast(centroid, testAtom);

					// if we're the shell
					if (dist <= radius) {
						if (!distances.containsKey(chainGroup))
							distances.put(chainGroup, Double.POSITIVE_INFINITY);
						if (useAverageDistance) {
							// sum the distance; we'll divide by the total
							// number later
							// here, we CANNOT use fastDistance (distance
							// squared) because we want the arithmetic mean
							distances.put(chainGroup, distances.get(chainGroup)
									+ Math.sqrt(dist));
							if (!atomCounts.containsKey(chainGroup))
								atomCounts.put(chainGroup, 0);
							atomCounts.put(chainGroup,
									atomCounts.get(chainGroup) + 1);
						} else {
							// take the minimum distance among all atoms of
							// chainGroup
							// note that we can't break here because we might
							// find a smaller distance
							if (dist < distances.get(chainGroup)) {
								distances.put(chainGroup, dist);
							}
						}
					}

				}
			}
		}

		if (useAverageDistance) {
			for (Map.Entry<Group, Double> entry : distances.entrySet()) {
				int count = atomCounts.get(entry.getKey());
				distances.put(entry.getKey(), entry.getValue() / count);
			}
		} else {
			// in this case we used getDistanceFast
			for (Map.Entry<Group, Double> entry : distances.entrySet()) {
				distances.put(entry.getKey(), Math.sqrt(entry.getValue()));
			}
		}

		return distances;

	}

	public static Set<Group> getGroupsWithinShell(Structure structure,
			Atom atom, Set<ResidueNumber> excludeResidues, double distance,
			boolean includeWater) {

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

		Set<Group> returnSet = new LinkedHashSet<Group>();
		for (Chain chain : structure.getChains()) {
			groupLoop: for (Group chainGroup : chain.getAtomGroups()) {
				if (!includeWater && chainGroup.isWater())
					continue;
				for (ResidueNumber rn : excludeResidues) {
					if (rn.equals(chainGroup.getResidueNumber()))
						continue groupLoop;
				}
				for (Atom atomB : chainGroup.getAtoms()) {

					// use getDistanceFast as we are doing a lot of comparisons
					double dist = Calc.getDistanceFast(atom, atomB);
					if (dist <= distance) {
						returnSet.add(chainGroup);
						break;
					}

				}
			}
		}
		return returnSet;
	}

	/**
	 * <p>
	 * Returns a Set of Groups in a structure within the distance specified of a
	 * given group.
	 * </p>
	 * <p>
	 * Updated 18-Sep-2015 sroughley to return a Set so only a unique set of
	 * Groups returned
	 *
	 * @param structure
	 *            The structure to work with
	 * @param group
	 *            The 'query' group
	 * @param distance
	 *            The cutoff distance
	 * @param includeWater
	 *            Should water residues be included in the output?
	 * @return {@link LinkedHashSet} of {@link Group}s within at least one atom
	 *         with {@code distance} of at least one atom in {@code group}
	 */
	public static Set<Group> getGroupsWithinShell(Structure structure,
			Group group, double distance, boolean includeWater) {

		Set<Group> returnList = new LinkedHashSet<Group>();

		Set<ResidueNumber> excludeGroups = new HashSet<ResidueNumber>();
		excludeGroups.add(group.getResidueNumber());
		for (Atom atom : group.getAtoms()) {
			Set<Group> set = getGroupsWithinShell(structure, atom,
					excludeGroups, distance, includeWater);
			returnList.addAll(set);
		}

		return returnList;
	}

	/**
	 * Remove all models from a Structure and keep only the first
	 *
	 * @param s
	 *            original Structure
	 * @return a structure that contains only the first model
	 * @since 3.0.5
	 */
	public static Structure removeModels(Structure s) {
		if (s.nrModels() == 1)
			return s;

		Structure n = new StructureImpl();
		// go through whole substructure and clone ...

		// copy structure data

		n.setPDBCode(s.getPDBCode());
		n.setName(s.getName());

		// TODO: do deep copying of data!
		n.setPDBHeader(s.getPDBHeader());
		n.setDBRefs(s.getDBRefs());

		n.setSites(s.getSites());

		n.setChains(s.getModel(0));

		return n;

	}

	/**
	 * Removes all polymeric and solvent groups from a list of groups
	 *
	 */
	public static List<Group> filterLigands(List<Group> allGroups) {

		List<Group> groups = new ArrayList<Group>();
		for (Group g : allGroups) {

			if ( g.isPolymeric())
				continue;

			if (!g.isWater()) {
				groups.add(g);
			}
		}

		return groups;
	}

	/**
	 * Short version of {@link #getStructure(String, PDBFileParser, AtomCache)}
	 * which creates new parsers when needed
	 *
	 * @param name
	 * @return
	 * @throws IOException
	 * @throws StructureException
	 */
	public static Structure getStructure(String name) throws IOException,
	StructureException {
		return StructureTools.getStructure(name, null, null);
	}

	/**
	 * Flexibly get a structure from an input String. The intent of this method
	 * is to allow any reasonable string which could refer to a structure to be
	 * correctly parsed. The following are currently supported:
	 * <ol>
	 * <li>Filename (if name refers to an existing file)
	 * <li>PDB ID
	 * <li>SCOP domains
	 * <li>PDP domains
	 * <li>Residue ranges
	 * <li>Other formats supported by AtomCache
	 * </ol>
	 *
	 * @param name
	 *            Some reference to the protein structure
	 * @param parser
	 *            A clean PDBFileParser to use if it is a file. If null, a
	 *            PDBFileParser will be instantiated if needed.
	 * @param cache
	 *            An AtomCache to use if the structure can be fetched from the
	 *            PDB. If null, a AtomCache will be instantiated if needed.
	 * @return A Structure object
	 * @throws IOException
	 *             if name is an existing file, but doesn't parse correctly
	 * @throws StructureException
	 *             if the format is unknown, or if AtomCache throws an
	 *             exception.
	 */
	public static Structure getStructure(String name, PDBFileParser parser,
			AtomCache cache) throws IOException, StructureException {
		File f = new File(FileDownloadUtils.expandUserHome(name));
		if (f.exists()) {
			if (parser == null) {
				parser = new PDBFileParser();
			}
			InputStream inStream = new FileInputStream(f);
			return parser.parsePDBFile(inStream);
		} else {
			if (cache == null) {
				cache = new AtomCache();
			}
			return cache.getStructure(name);
		}
	}

	/**
	 * @deprecated  use {@link Chain#isProtein()} instead.
	 */
	@Deprecated
	public static boolean isProtein(Chain c) {

		return c.isProtein();
	}

	/**
	 * @deprecated use {@link Chain#isNucleicAcid()} instead.
	 */
	@Deprecated
	public static boolean isNucleicAcid(Chain c) {
		return c.isNucleicAcid();
	}

	/**
	 * @deprecated use {@link Chain#getPredominantGroupType()} instead.
	 */
	@Deprecated
	public static GroupType getPredominantGroupType(Chain c) {
		return c.getPredominantGroupType();
	}

	/**
	 * @deprecated use {@link Chain#isWaterOnly()} instead.
	 */
	@Deprecated
	public static boolean isChainWaterOnly(Chain c) {
		return c.isWaterOnly();
	}

	/**
	 * @deprecated  use {@link Chain#isPureNonPolymer()} instead.
	 */
	@Deprecated
	public static boolean isChainPureNonPolymer(Chain c) {

		return c.isPureNonPolymer();
	}

	/**
	 * Cleans up the structure's alternate location (altloc) groups. All alternate location groups should have all atoms (except
	 * in the case of microheterogenity) or when a deuterium exists.
	 * Ensure that all the alt loc groups have all the atoms in the main group.
	 * @param structure The Structure to be cleaned up
	 */
	public static void cleanUpAltLocs(Structure structure) {
		for (int i =0; i< structure.nrModels() ; i++){
			for (Chain chain : structure.getModel(i)) {
				for (Group group : chain.getAtomGroups()) {
					for (Group altLocGroup : group.getAltLocs()) {
						for ( Atom groupAtom : group.getAtoms()) {
							// If this alt loc doesn't have this atom
							if (! altLocGroup.hasAtom(groupAtom.getName())) {
								// Fix for microheterogenity
								if (altLocGroup.getPDBName().equals(group.getPDBName())) {
									// If it's a Hydrogen then we check for it's Deuterated brother
									if(!hasDeuteratedEquiv(groupAtom, altLocGroup)){
										altLocGroup.addAtom(groupAtom);
									}
								}
							}
						}
					}
				}
			}
		}
	}

	/**
	 * Check to see if an Deuterated atom has a non deuterated brother in the group.
	 * @param atom the input atom that is putatively deuterium
	 * @param currentGroup the group the atom is in
	 * @return true if the atom is deuterated and it's hydrogen equive exists.
	 */
	public static boolean hasNonDeuteratedEquiv(Atom atom, Group currentGroup) {
		if(atom.getElement()==Element.D && currentGroup.hasAtom(replaceFirstChar(atom.getName(),'D', 'H'))) {
			// If it's deuterated and has a non-deuterated brother
			return true;
		}
		return false;
	}

	/**
	 * Check to see if a Hydrogen has a  Deuterated brother in the group.
	 * @param atom the input atom that is putatively hydorgen
	 * @param currentGroup the group the atom is in
	 * @return true if the atom is hydrogen and it's Deuterium equiv exists.
	 */
	public static boolean hasDeuteratedEquiv(Atom atom, Group currentGroup) {
		if(atom.getElement()==Element.H && currentGroup.hasAtom(replaceFirstChar(atom.getName(),'H', 'D'))) {
			// If it's hydrogen and has a deuterated brother
			return true;
		}
		return false;
	}

	private static String replaceFirstChar(String name, char c, char d) {
		if(name.charAt(0)==c){
			return name.replaceFirst(String.valueOf(c), String.valueOf(d));
		}
		return name;
	}

	/**
	 * Remove all atoms but the representative atoms (C alphas or phosphates) from the given structure.
	 * @param structure the structure
	 * @since 5.4.0
	 */
	public static void reduceToRepresentativeAtoms(Structure structure) {
		for (int modelIdx = 0; modelIdx<structure.nrModels(); modelIdx++) {
			for (Chain c : structure.getPolyChains(modelIdx)) {
				for (Group g : c.getAtomGroups()) {
					List<Atom> atoms = g.getAtoms();
					if (g.isAminoAcid()) {
						atoms.removeIf(a->!a.getName().equals(CA_ATOM_NAME));
					} else if (g.isNucleotide()) {
						atoms.removeIf(a->!a.getName().equals(NUCLEOTIDE_REPRESENTATIVE));
					}
					// else we keep all other atoms. We are concerned only about aminoacids and nucleotides that make up the bulk of the structures
				}
			}
		}
	}
}