helper function to determine whether a trace break should be introduced between two residues of the same type (amino acid, or nucleotides).  aaStretch: indicates whether the residues are to be treated as amino acids
425249
function shouldIntroduceTraceBreak(aaStretch, prevResidue, thisResidue) {
  // these checks are on purpose more relaxed than the checks we use in 
  // deriveConnectivity(). We don't really care about correctness of bond 
  // lengths here. The only thing that matters is that the residues are 
  // more or less close so that they could potentially be connected.
  var prevAtom, thisAtom;
  if (aaStretch) {
    prevAtom = prevResidue.atom('C');
    thisAtom = thisResidue.atom('N');
  } else {
    prevAtom = prevResidue.atom('O3\'');
    thisAtom = thisResidue.atom('P');
  }

  // in case there is a bond, we don't introduce a chain break
  if (prevAtom.isConnectedTo(thisAtom)) {
    return false;
  }
  var sqrDist = vec3.sqrDist(prevAtom.pos(), thisAtom.pos());
  return (Math.abs(sqrDist - 1.5*1.5) > 1);
}
785335
def check_structure_chain_quality(self, seqprop, structprop, chain_id,
                                      seq_ident_cutoff=0.5, allow_missing_on_termini=0.2,
                                      allow_mutants=True, allow_deletions=False,
                                      allow_insertions=False, allow_unresolved=True):
        """"""
        alignment = self._get_seqprop_to_structprop_alignment(seqprop=seqprop, structprop=structprop, chain_id=chain_id)

        # Compare sequence to structure's sequence using the alignment
        chain_passes_quality_check = ssbio.protein.structure.properties.quality.sequence_checker(reference_seq_aln=alignment[0],
                                                                   structure_seq_aln=alignment[1],
                                                                   seq_ident_cutoff=seq_ident_cutoff,
                                                                   allow_missing_on_termini=allow_missing_on_termini,
                                                                   allow_mutants=allow_mutants,
                                                                   allow_deletions=allow_deletions,
                                                                   allow_insertions=allow_insertions,
                                                                   allow_unresolved=allow_unresolved)
        return chain_passes_quality_check
877220
public static int[][] aaBondInfo() {

        if (aminoAcids == null) {
            createAAs();
        }

        int[][] info = new int[153][4];

        int counter = 0;
        int total = 0;
        for (int aa = 0; aa < aminoAcids.length; aa++) {
            AminoAcid acid = aminoAcids[aa];

            LOGGER.debug("#bonds for ", acid.getProperty(RESIDUE_NAME).toString(), " = " + acid.getBondCount());
            total += acid.getBondCount();
            LOGGER.debug("total #bonds: ", total);

            Iterator<IBond> bonds = acid.bonds().iterator();
            while (bonds.hasNext()) {
                IBond bond = (IBond) bonds.next();
                info[counter][0] = counter;
                info[counter][1] = acid.indexOf(bond.getBegin());
                info[counter][2] = acid.indexOf(bond.getEnd());
                info[counter][3] = bond.getOrder().numeric();
                counter++;
            }
        }

        if (counter > 153) {
            LOGGER.error("Error while creating AA info! Bond count is too large: ", counter);
            return null;
        }

        return info;
    }
678223
public static AFPChain fastaToAfpChain(SequencePair<Sequence<AminoAcidCompound>, AminoAcidCompound> alignment,
			Structure structure1, Structure structure2) throws StructureException {
		List<AlignedSequence<Sequence<AminoAcidCompound>, AminoAcidCompound>> seqs = alignment.getAlignedSequences();
		StringBuilder sb1 = new StringBuilder();
		for (AminoAcidCompound a : seqs.get(0)) {
			sb1.append(a.getBase());
		}
		try {
			ProteinSequence seq1 = new ProteinSequence(sb1.toString());
			StringBuilder sb2 = new StringBuilder();
			for (AminoAcidCompound a : seqs.get(1)) {
				sb1.append(a.getBase());
			}
			ProteinSequence seq2 = new ProteinSequence(sb2.toString());
			LinkedHashMap<String, ProteinSequence> map = new LinkedHashMap<String, ProteinSequence>();
			map.put(structure1.getName(), seq1);
			map.put(structure2.getName(), seq2);
			return fastaToAfpChain(map, structure1, structure2);
		} catch (CompoundNotFoundException e) {
			logger.error("Unexpected error while creating protein sequences: {}. This is most likely a bug.",e.getMessage());
			return null;
		}
	}
1295367
def findOrDie(s):
    
    aa = find(s)
    if aa:
        return aa
    else:
        print('Unknown amino acid or codon: %s' % s, file=sys.stderr)
        print('Valid arguments are: %s.' % list(CODONS.keys()),
              file=sys.stderr)
        sys.exit(1)
1348302
public void identify(final List<Chain> chains,
			final Set<ProteinModification> potentialModifications) {

		if (chains==null) {
			throw new IllegalArgumentException("Null structure.");
		}

		if (potentialModifications==null) {
			throw new IllegalArgumentException("Null potentialModifications.");
		}


		reset();

		if (potentialModifications.isEmpty()) {
			return;
		}


		residues = new ArrayList<Group>();
		List<Group> ligands = new ArrayList<Group>();
		Map<Component, Set<Group>> mapCompGroups = new HashMap<Component, Set<Group>>();

		for (Chain chain : chains) {

			List<Group> ress = StructureUtil.getAminoAcids(chain);

			//List<Group> ligs = chain.getAtomLigands();
			List<Group> ligs = StructureTools.filterLigands(chain.getAtomGroups());
			residues.addAll(ress);
			residues.removeAll(ligs);
			ligands.addAll(ligs);
			addModificationGroups(potentialModifications, ress, ligs, mapCompGroups);
		}

		if (residues.isEmpty()) {
			String pdbId = "?";
			if ( chains.size() > 0) {
				Structure struc = chains.get(0).getStructure();
				if ( struc != null)
					pdbId = struc.getPDBCode();
			}
			logger.warn("No amino acids found for {}. Either you did not parse the PDB file with alignSEQRES records, or this record does not contain any amino acids.", pdbId);
		}
		List<ModifiedCompound> modComps = new ArrayList<ModifiedCompound>();

		for (ProteinModification mod : potentialModifications) {
			ModificationCondition condition = mod.getCondition();
			List<Component> components = condition.getComponents();
			if (!mapCompGroups.keySet().containsAll(components)) {
				// not all components exist for this mod.
				continue;
			}

			int sizeComps = components.size();
			if (sizeComps==1) {

				processCrosslink1(mapCompGroups, modComps, mod, components);

			} else {

				processMultiCrosslink(mapCompGroups, modComps, mod, condition);
			}
		}

		if (recordAdditionalAttachments) {
			// identify additional groups that are not directly attached to amino acids.
			for (ModifiedCompound mc : modComps) {
				identifyAdditionalAttachments(mc, ligands, chains);
			}
		}

		mergeModComps(modComps);

		identifiedModifiedCompounds.addAll(modComps);


		// record unidentifiable linkage
		if (recordUnidentifiableModifiedCompounds) {
			recordUnidentifiableAtomLinkages(modComps, ligands);
			recordUnidentifiableModifiedResidues(modComps);
		}
	}
735422
def _validate_residue_types(self):
        ''''''

        for chain_id, sequence_map in self.rosetta_to_atom_sequence_maps.iteritems():
            rosetta_sequence = self.rosetta_sequences[chain_id]
            atom_sequence = self.atom_sequences[chain_id]
            for rosetta_id, atom_id, _ in sequence_map:
                assert(rosetta_sequence[rosetta_id].ResidueAA == atom_sequence[atom_id].ResidueAA)

        for chain_id, sequence_map in self.atom_to_seqres_sequence_maps.iteritems():
            atom_sequence = self.atom_sequences[chain_id]
            seqres_sequence = self.seqres_sequences[chain_id]
            for atom_id, seqres_id, _ in sorted(sequence_map):
                assert(atom_sequence[atom_id].ResidueAA == seqres_sequence[seqres_id].ResidueAA)

        for chain_id, sequence_map in self.seqres_to_uniparc_sequence_maps.iteritems():
            if self.pdb_chain_to_uniparc_chain_mapping.get(chain_id):
                seqres_sequence = self.seqres_sequences[chain_id]
                uniparc_sequence = self.uniparc_sequences[self.pdb_chain_to_uniparc_chain_mapping[chain_id]]
                for seqres_id, uniparc_id_resid_pair, substitution_match in sequence_map:
                    uniparc_id = uniparc_id_resid_pair[1]
                    # Some of the matches may not be identical but all the '*' Clustal Omega matches should be identical
                    if substitution_match and substitution_match.clustal == 1:
                        assert(seqres_sequence[seqres_id].ResidueAA == uniparc_sequence[uniparc_id].ResidueAA)
250662
private List<List<Atom[]>> getMatchedAtomsOfLinkages(
			ModificationCondition condition, Map<Component, Set<Group>> mapCompGroups) {
		List<ModificationLinkage> linkages = condition.getLinkages();
		int nLink = linkages.size();

		List<List<Atom[]>> matchedAtomsOfLinkages =
				new ArrayList<List<Atom[]>>(nLink);

		for (int iLink=0; iLink<nLink; iLink++) {
			ModificationLinkage linkage = linkages.get(iLink);
			Component comp1 = linkage.getComponent1();
			Component comp2 = linkage.getComponent2();

//			boolean isAA1 = comp1.;
//			boolean isAA2 = comp2.getType()==true;

			Set<Group> groups1 = mapCompGroups.get(comp1);
			Set<Group> groups2 = mapCompGroups.get(comp2);

			List<Atom[]> list = new ArrayList<Atom[]>();

			List<String> potentialNamesOfAtomOnGroup1 = linkage.getPDBNameOfPotentialAtomsOnComponent1();
			for (String name : potentialNamesOfAtomOnGroup1) {
				if (name.equals("*")) {
					// wildcard
					potentialNamesOfAtomOnGroup1 = null; // search all atoms
					break;
				}
			}

			List<String> potentialNamesOfAtomOnGroup2 = linkage.getPDBNameOfPotentialAtomsOnComponent2();
			for (String name : potentialNamesOfAtomOnGroup2) {
				if (name.equals("*")) {
					// wildcard
					potentialNamesOfAtomOnGroup2 = null; // search all atoms
					break;
				}
			}

			for (Group g1 : groups1) {
				for (Group g2 : groups2) {
					if (g1.equals(g2)) {
						continue;
					}

					// only for wildcard match of two residues
					boolean ignoreNCLinkage =
						potentialNamesOfAtomOnGroup1 == null &&
						potentialNamesOfAtomOnGroup2 == null &&
						residues.contains(g1) &&
						residues.contains(g2);

					Atom[] atoms = StructureUtil.findNearestAtomLinkage(
							g1, g2,
							potentialNamesOfAtomOnGroup1,
							potentialNamesOfAtomOnGroup2,
							ignoreNCLinkage,
							bondLengthTolerance);
					if (atoms!=null) {
						list.add(atoms);
					}
				}
			}

			if (list.isEmpty()) {
				// broken linkage
				break;
			}

			matchedAtomsOfLinkages.add(list);
		}

		return matchedAtomsOfLinkages;
	}
525429
private static void traceBack(int[] pre, int currafp0, int twist, FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[]ca2)
	{

		afpChain.setConn(0d);
		afpChain.setDVar(0d);

		int minLen = afpChain.getMinLen();
		List<AFP> afpSet = afpChain.getAfpSet();

		int afpChainLen = 0;

		//trace-back from currafp (maxsco)
		int[]     afpchain    = new int[minLen];
		int[]     afptwibin   = new int[minLen];
		double[]  afptwilist  = new double[minLen];
		int       currafp = currafp0;
		int       s = 0;

		afptwibin[s] = 0;
		afpchain[s ++] = currafp;

		boolean isConnected = false;
		int     prevafp;
		// int     alnlen = afpSet.get(afpchain[s]).getFragLen();

		//Double  conn = new Double(0) ;
		//Double dvar =  new Double(0);

		while((prevafp = pre[currafp]) != -1)   {

			isConnected = afpPairConn(prevafp, currafp, params,afpChain);

			if ( isConnected )
				afptwibin[s - 1] = 1;
			else
				afptwibin[s - 1] = 0;
			Double dvar = afpChain.getDVar();
			afptwilist[s - 1] = dvar;
			//note s - 1: the transformation of prevafp-currafp is recorded in currafp
			currafp = prevafp;
			// alnlen += afpSet.get(currafp).getFragLen();
			afpchain[s ++] = currafp;
		}

		afpChainLen = s;
		afpChain.setAfpChainLen(afpChainLen);

		//first afp without transformation
		if ( isConnected)
			afptwibin[s - 1] = 1;
		else
			afptwibin[s - 1] = 0;

		//if(debug)
		//   System.out.println(String.format("including %d AFPs, %d residues\n", afpChainLen, alnlen));

		//record the optimal alignment in afpChainList (afpChainLen)

		int[] afpChainList = afpChain.getAfpChainList();
		double[] afpChainTwiBin = afpChain.getAfpChainTwiBin();
		double[] afpChainTwiList = afpChain.getAfpChainTwiList();

		if(afpChainList == null)        {
			afpChainList   = new int[s];
			afpChain.setAfpChainList(afpChainList);
			afpChainTwiBin     = new double[s];
			afpChain.setAfpChainTwiBin(afpChainTwiBin);
			afpChainTwiList = new double[s];
			afpChain.setAfpChainTwiList(afpChainTwiList);
		}

		int afpChainTwiNum = afpChain.getAfpChainTwiNum();

		int     i;
		for(i = 0; i < s; i ++) {
			afpChainList[i]     = afpchain[s - 1 - i];
			afpChainTwiBin[i]   = afptwibin[s - 1 - i];
			afpChainTwiList[i]  = afptwilist[s - 1 - i];
			afpChainTwiNum      += afptwibin[s - 1 - i];
		}

		if(afpChainTwiNum != twist)     {
			System.err.println(String.format("AFPChainer Warning: the twists number is not consistent %d %d\n", afpChainTwiNum, twist));
		}

		double alignScore = afpChain.getAlignScore();

		double  checkscore = afpSet.get(afpChainList[0]).getScore();
		for(i = 1; i < afpChainLen; i ++)       {
			isConnected = afpPairConn(afpChainList[i - 1], afpChainList[i], params,afpChain);
			checkscore = checkscore + afpSet.get(afpChainList[i]).getScore() + afpChain.getConn();
		}
		if(Math.abs(checkscore - alignScore) > 1e-4)        {
			System.err.println(String.format("AFPChainer Warning: fail in alignment score checking %.4f %.4f\n", alignScore, checkscore));
		}

		if ( debug )
			System.out.println("traceBack:" + afpChainLen + " " + afpChainList.length);
		double  rmsd = calAfpRmsd(afpChainLen, afpChainList,0, afpChain,ca1,ca2);
		afpChain.setChainRmsd(rmsd);
		if (debug )
			System.out.println("Chain RMSD: " + rmsd);
		int     b1 = 0;
		int     bk = 0;
		int     a, b;
		afpChain.setChainLen( 0);
		int chainLen       = afpChain.getChainLen();
		int[] block2Afp = afpChain.getBlock2Afp();


		double[] blockRmsd = afpChain.getBlockRmsd();
		int[] blockSize = afpChain.getBlockSize();

		block2Afp[0] = 0;
		for(i = 0; i < afpChainLen; i ++)       {
			a = afpChainList[i];
			chainLen += afpSet.get(a).getFragLen();
			if(i > 0)       {
				b = afpChainList[i - 1];
				int misLen = afpChain.getMisLen();
				misLen += calcMismatch(afpSet.get(a),afpSet.get(b));
				afpChain.setMisLen(misLen);
				int gapLen = afpChain.getGapLen();
				gapLen += calcGap(afpSet.get(a),afpSet.get(b));
				afpChain.setGapLen(gapLen);

			}

			if(afpChainTwiBin[i] == 1)   {
				if (debug)
					System.out.println(" ** calAfpTmsd : afpChainWtiBin == 1 : i: "+i+" i-b1: " + (i-b1) + " b1: " + b1 + " afpChainList.len: " + afpChainList.length);

				//int len = afpChainList.length - b1 +1;
				//int[] fakeList = new int[len];
				//int pos = -1;
				//for ( int fPos = b1 ; fPos< afpChainList.length ; fPos++){
				//   pos++;
				//   fakeList[pos] = afpChainList[fPos];
				//}
				if (debug )
					System.err.println("calculation calAfpRmsd " + i + " " + b1 + " " );


				rmsd = calAfpRmsd(i - b1, afpChainList,b1, afpChain, ca1, ca2);


				blockRmsd[bk] = rmsd;
				blockSize[bk] = i - b1;
				b1 = i;

				//System.out.println("block2Afp.length:"+ block2Afp.length + " " + bk + " " + i + " " + afpChain.getMaxTra() );
				block2Afp[++bk] = i; //next-block
			}
		}
		afpChain.setBlock2Afp(block2Afp);
		afpChain.setChainLen(chainLen);

		if (debug)
			System.out.println("after loop over all afpChainList " + (i-b1) + " " + b1);

		rmsd = calAfpRmsd(i - b1, afpChainList, b1, afpChain,ca1,ca2);
		if (debug)
			System.out.println("*** i:" + i + " b1: " + b1 + " i-b1 " + (i-b1));



		// argh this is the block RMSD, not the chain RMSD!
		//afpChain.setChainRmsd(rmsd);
		//rmsd = calAfpRmsd(i - b1, afpChainList[b1]);
		blockSize[bk] = i - b1;
		blockRmsd[bk] = rmsd;

		afpChain.setBlockSize(blockSize);
		afpChain.setBlockRmsd(blockRmsd);
		int blockNum = afpChain.getBlockNum();
		blockNum = ++bk;
		if ( debug)
			System.err.println("AFPChainser setBlockNUm:" + blockNum);
		afpChain.setBlockNum(blockNum);
		afpChain.setAfpChainList(afpChainList);
		afpChain.setAfpChainTwiList(afpChainTwiList);

	}
1331152
def prune_structure_according_to_loops_file(self, loops_file_content, search_radius, expected_loop_length = None, generate_pymol_session = True):
        ''''''
        lf = LoopsFile(loops_file_content)
        assert(len(lf.data) == 1) # todo: remove

        # Extract the list of Residues defined in the loops file
        loop_residues = []
        for loop in lf.data:
            # Identify the Residues corresponding to the loops file definition
            start_id = str(loop['start'])
            end_id = str(loop['end'])
            start = []
            end = []
            for chain, chain_residues in self.residues.iteritems():
                for res, residue_object in chain_residues.iteritems():
                    if res.strip() == start_id:
                        start.append((chain, res, residue_object))
                    if res.strip() == end_id:
                        end.append((chain, res, residue_object))
            if len(start) != len(end) != 1:
                raise Exception('The PDB is ambiguous with respect to the loops file i.e. more than one PDB residue corresponds to the start or end residue defined in the loops file.')
            start, end = start[0], end[0]

            # We assume that the loops file identifies a range in document order
            loop_start = False
            for l in self.indexed_lines:
                if l[0] == 'ATOM' or l[0] == 'HETATM':
                    atom_residue = l[3].residue
                    if not loop_start:
                        if atom_residue == start[2]:
                            loop_start = True
                            if atom_residue not in loop_residues:
                                loop_residues.append(atom_residue)
                    else:
                        if atom_residue not in loop_residues:
                            loop_residues.append(atom_residue)
                        if atom_residue == end[2]:
                            break
        if expected_loop_length != None and expected_loop_length != len(loop_residues):
            raise Exception('Expected to identify {0} residues but {1} were identified.'.format(expected_loop_length, len(loop_residues)))

        # Determine the sidechain atoms to be removed
        sidechain_atom_serial_numbers = self.find_sidechain_atoms_within_radius_of_residue_objects(loop_residues, search_radius)

        # Determine the loop residue atoms to be removed
        loop_residues_ids = set()
        loop_atom_serial_numbers = set()
        for loop_residue in loop_residues:
            for aatom in loop_residue.get('ATOM'):
                loop_atom_serial_numbers.add(aatom.serial_number)
        assert(len(sidechain_atom_serial_numbers.intersection(loop_atom_serial_numbers)) == 0)

        return self.prune(loop_atom_serial_numbers, sidechain_atom_serial_numbers, generate_pymol_session = generate_pymol_session)
269440
private boolean updateTraceSpec(TraceSpecification ts) {
        List<TraceElement> traceSpecs = ts.getSpecs();
        Integer minimumLevel = null;
        if (ts.isSensitiveTraceSuppressed()) {
            minimumLevel = findMinimumSafeLevel(ts.getSafeLevelsIndex());
        }

        int newFineLevel = fineLevel;
        int newFineLevelsEnabled = 0;
        int newSpecTraceLevel = TrLevelConstants.SPEC_TRACE_LEVEL_OFF;

        for (TraceElement spec : traceSpecs) {
            String clazz = spec.groupName;
            int traceElementFineLevel = spec.fineLevel;
            int specTraceLevel = spec.specTraceLevel;
            boolean setValue = spec.action;

            // Do we have a match with our package name?
            boolean process = false;

            if (clazz.endsWith("*")) // packages can end with wildcard
            {
                if (1 == clazz.length()) {
                    process = true;
                } else {
                    clazz = clazz.substring(0, clazz.length() - 1);
                    for (String group : groups) {
                        if (group.startsWith(clazz)) {
                            process = true;
                            break;
                        }
                    }
                    process = process || name.startsWith(clazz);
                }
            } else {
                // or can be a complete package or parent package
                // - look for a full or partial match with the TC package name
                int lastDot = name.lastIndexOf('.');
                if (lastDot > 0) {
                    String packageName = name.substring(0, lastDot);
                    if (packageName.startsWith(clazz))
                        process = true;
                }
                // groups may be class names (eg WsLogger impls)
                for (String group : groups) {
                    lastDot = group.lastIndexOf('.');
                    if (lastDot > 0) {
                        String packageName = group.substring(0, lastDot);
                        if (packageName.startsWith(clazz))
                            process = true;
                        break;
                    }
                }

                // could be a straight group name match
                for (String group : groups) {
                    if (group.equalsIgnoreCase(clazz)) {
                        process = true;
                        break;
                    }
                }
                process = process || name.equalsIgnoreCase(clazz);
            }

            if (process) {
                newFineLevel = traceElementFineLevel;
                newSpecTraceLevel = specTraceLevel;
                if (minimumLevel != null && newFineLevel < minimumLevel) {
                    newFineLevel = minimumLevel;
                }

                for (int level = newFineLevel; level < TrLevelConstants.TRACE_LEVEL_OFF; level++) {
                    if (setValue) {
                        newFineLevelsEnabled |= 1 << level;
                    } else {
                        newFineLevelsEnabled &= ~(1 << level);
                    }
                }
                
                if (newFineLevel == TrLevelConstants.TRACE_LEVEL_OFF) {
                	isTraceOff = true;
                } else {
                	isTraceOff = false;
                }
                // Indicate that the trace spec matched something
                spec.setMatched(true);
            }
        } // end for each spec

        boolean updated = false;
        if (newFineLevel != fineLevel) {
            fineLevel = newFineLevel;
            updated = true;
        }
        if (newFineLevelsEnabled != fineLevelsEnabled) {
            fineLevelsEnabled = newFineLevelsEnabled;
            updated = true;
        }
        if (newSpecTraceLevel != specTraceLevel) {
            specTraceLevel = newSpecTraceLevel;
            updated = true;
        }

        return updated;
    }
562936
public static boolean isAmbiguosProtein(String sequence) {
	sequence = SequenceUtil.cleanSequence(sequence);
	if (SequenceUtil.isNonAmbNucleotideSequence(sequence)) {
		return false;
	}
	if (SequenceUtil.DIGIT.matcher(sequence).find()) {
		return false;
	}
	if (SequenceUtil.NON_AA.matcher(sequence).find()) {
		return false;
	}
	if (SequenceUtil.AA.matcher(sequence).find()) {
		return false;
	}
	final Matcher amb_prot = SequenceUtil.AMBIGUOUS_AA.matcher(sequence);
	return amb_prot.find();
	}
278376
boolean matchNext() {
        if (numAtoms == 0)
            return false;
        if (sptr > 1)
            backtrack();
        main:
        while (sptr != 0) {
            final int bidx = currBondIdx();

            if (bidx == numBonds) {

                // done
                if (numMapped == numAtoms)
                    return true;

                // handle disconnected atoms
                for (IAtom qatom : query.atoms()) {
                    if (amap[qatom.getIndex()] == UNMAPPED) {
                        Iterator<IAtom> iter = atoms();
                        while (iter.hasNext()) {
                            IAtom atom = iter.next();
                            if (feasible(bidx, (IQueryAtom) qatom, atom))
                                continue main;
                        }
                        break;
                    }
                }
                backtrack();
                continue;
            }

            IQueryBond qbond = qbonds[bidx];
            IQueryAtom qbeg  = (IQueryAtom) qbond.getBegin();
            IQueryAtom qend  = (IQueryAtom) qbond.getEnd();

            int begIdx = amap[qbeg.getIndex()];
            int endIdx = amap[qend.getIndex()];

            // both atoms matched, there must be a bond between them
            if (begIdx != UNMAPPED && endIdx != UNMAPPED) {
                IBond bond = mol.getAtom(begIdx).getBond(mol.getAtom(endIdx));
                if (feasible(qbond, bond))
                    continue;
            }
            // 'beg' is mapped, find a feasible 'end' from it's neighbor list
            else if (begIdx != UNMAPPED) {
                IAtom           beg   = mol.getAtom(begIdx);
                Iterator<IBond> biter = bonds(beg);
                while (biter.hasNext()) {
                    IBond bond = biter.next();
                    IAtom end  = bond.getOther(beg);
                    if (qbond.matches(bond) && feasible(bidx + 1, qend, end))
                        continue main;
                }
            }
            // 'end' is mapped, find a feasible 'beg' from it's neighbor list
            else if (endIdx != UNMAPPED) {
                IAtom           end   = mol.getAtom(endIdx);
                Iterator<IBond> biter = bonds(end);
                while (biter.hasNext()) {
                    IBond bond = biter.next();
                    IAtom beg  = bond.getOther(end);
                    if (qbond.matches(bond) && feasible(bidx + 1, qbeg, beg))
                        continue main;
                }
            }
            // 'beg' nor 'end' matched, find a feasible mapping from
            // any atom in the molecule
            else {
                Iterator<IAtom> aiter = atoms();
                while (aiter.hasNext()) {
                    if (feasible(bidx, qbeg, aiter.next()))
                        continue main;
                }
            }
            backtrack();
        }
        return false;
    }
194137
@Override
	public List<Sequence<AminoAcidCompound>> createSequences(
			Sequence<NucleotideCompound> originalSequence) {

		List<List<AminoAcidCompound>> workingList = new ArrayList<List<AminoAcidCompound>>();

		Iterable<SequenceView<NucleotideCompound>> iter = new WindowedSequence<NucleotideCompound>(
				originalSequence, 3);

		boolean first = true;

		// If not waiting for a start codon, start translating immediately
		boolean doTranslate = !waitForStartCodon;

		for (SequenceView<NucleotideCompound> element : iter) {
			AminoAcidCompound aminoAcid = null;

			int i = 1;
			Table.CaseInsensitiveTriplet triplet = new Table.CaseInsensitiveTriplet(
					element.getCompoundAt(i++), element.getCompoundAt(i++),
					element.getCompoundAt(i++));

			Codon target = null;

			target = quickLookup.get(triplet);

			// Check for a start
			if (!doTranslate && target.isStart()) {
				doTranslate = true;
			}

			if (doTranslate) {
				if (target != null)
					aminoAcid = target.getAminoAcid();
				if (aminoAcid == null && translateNCodons()) {
					aminoAcid = unknownAminoAcidCompound;
				} else {
					if (first && initMetOnly && target.isStart()) {
						aminoAcid = methionineAminoAcidCompound;
					}
				}

				addCompoundsToList(Arrays.asList(aminoAcid), workingList);
			}

			if (doTranslate && stopAtStopCodons && target.isStop()) {
				// Check if we need to stop, but dont stop until started!
				break;
			}

			first = false;
		}
		postProcessCompoundLists(workingList);

		return workingListToSequences(workingList);
	}
265508
private List<Group> trySimpleMatch(List<Group> seqResGroups,List<Group> atmResGroups) {
		// by default first ATOM position is 1
		//

		@SuppressWarnings("unchecked")
		List<Group> newSeqResGroups = (ArrayList<Group>)((ArrayList<Group>)seqResGroups).clone();

		boolean startAt1 = true;

		for ( int atomResPos = 0 ; atomResPos < atmResGroups.size() ; atomResPos++){

			// let's try to match this case
			Group atomResGroup = atmResGroups.get(atomResPos);

			// let's ignore waters
			if ( atomResGroup.isWater()){
				continue;
			}

			ResidueNumber atomResNum = atomResGroup.getResidueNumber();

			int seqResPos = atomResNum.getSeqNum();



			if ( seqResPos < 0) {
				logger.debug("ATOM residue number < 0 : {}", seqResPos);
				return null;
			}

			if ( seqResPos == 0){
				// make sure the first SEQRES is matching.
				Group seqResGroup = seqResGroups.get(0);
				if (  seqResGroup.getPDBName().equals(atomResGroup.getPDBName())){
					// they match! seems in this case the numbering starts with 0...
					startAt1 = false;
				} else {

					logger.debug("SEQRES position 1 ({}) does not match ATOM PDB res num 0 ({})",
							seqResGroup.getPDBName(), atomResGroup.getPDBName());


					return null;

				}
			}


			if ( startAt1 )
				seqResPos--;

			// another check that would require the alignment approach
			if ( startAt1 && seqResPos >=  seqResGroups.size()  )
			{

				// could be a HETATOM...
				if ( atomResGroup instanceof AminoAcid) {
					logger.debug(" ATOM residue nr: " + seqResPos + " > seqres! " + seqResGroups.size() + " " + atomResGroup);
					return null;
				} else if ( atomResGroup instanceof NucleotideImpl) {
					logger.debug(" NUCLEOTIDE residue nr: " + seqResPos + " > seqres! " + seqResGroups.size() + " " + atomResGroup);
					return null;
				} else {
					// we won't map HETATOM groups...
					continue;
				}
			}


			//			if ( seqResPos < 0){
			//
			//				System.err.println("What is going on??? " + atomRes.getChainID() + " " + atomResGroup);
			//			}

			if ( seqResPos >= seqResGroups.size()){
				logger.debug("seqres groups don't match atom indices " + seqResPos);
				if ( atomResGroup instanceof AminoAcid )
					return null;
				else
					continue;
			}

			Group seqResGroup = seqResGroups.get(seqResPos );

			if ( ! seqResGroup.getPDBName().trim().equals(atomResGroup.getPDBName().trim())){
				// a mismatch! something is wrong in the mapping and we need to do an alignment
				logger.debug("Mismatch of SEQRES pos " + seqResPos + " and ATOM record: " + atomResGroup + " | " + seqResGroup);
				return null;
			}

			// the two groups are identical and we can merge them
			// replace the SEQRES group with the ATOM group...

			Group replacedGroup = newSeqResGroups.set(seqResPos, atomResGroup);
			logger.debug("Merging index {}: replaced seqres group {} ({}) with atom group {} ({})",
					seqResPos,
					replacedGroup.getResidueNumber(), replacedGroup.getPDBName(),
					atomResGroup.getResidueNumber(), atomResGroup.getPDBName());

		}

		// all went ok. copy over the updated list to the original one.
		// note: if something went wrong, we did not modifiy the original list.
		//seqResGroups = newSeqResGroups;


		//			int pos = -1;
		//			for (Group g: seqResGroups){
		//				pos++;
		//				logger.debug(pos + " " + g);
		//			}

		//System.out.println("I:" + seqResGroups);
		// all atom records could get matched correctly!
		return newSeqResGroups;

	}
697579
function checkPair( bodyA, bodyB ){

        // filter out bodies that don't collide with each other
        if (
            ( bodyA.treatment === 'static' || bodyA.treatment === 'kinematic' ) &&
            ( bodyB.treatment === 'static' || bodyB.treatment === 'kinematic' )
        ){
            return false;
        }

        if ( bodyA.geometry.name === 'circle' && bodyB.geometry.name === 'circle' ){

            return checkCircles( bodyA, bodyB );

        } else if ( bodyA.geometry.name === 'compound' || bodyB.geometry.name === 'compound' ){
            // compound bodies are special. We can't use gjk because
            // they could have concavities. so we do the pieces individually
            var test = (bodyA.geometry.name === 'compound')
                ,compound = test ? bodyA : bodyB
                ,other = test ? bodyB : bodyA
                ,cols
                ,ch
                ,ret = []
                ,scratch = Physics.scratchpad()
                ,vec = scratch.vector()
                ,oldPos = scratch.vector()
                ,otherAABB = other.aabb()
                ,i
                ,l
                ;

            for ( i = 0, l = compound.children.length; i < l; i++ ){

                ch = compound.children[ i ];
                // move body to fake position
                oldPos.clone( ch.state.pos );
                ch.offset.vadd( oldPos.vadd( compound.offset ).rotate( -ch.state.angular.pos ) );
                ch.state.pos.clone( compound.state.pos );
                ch.state.angular.pos += compound.state.angular.pos;

                // check it if the aabbs overlap
                if ( Physics.aabb.overlap(otherAABB, ch.aabb()) ){

                    cols = checkPair( other, ch );

                    if ( cols instanceof Array ){
                        for ( var j = 0, c, ll = cols.length; j < ll; j++ ){
                            c = cols[j];
                            // set body to be the compound body
                            if ( c.bodyA === ch ){
                                c.bodyA = compound;
                            } else {
                                c.bodyB = compound;
                            }
                            ret.push( c );
                        }

                    } else if ( cols ) {
                        // set body to be the compound body
                        if ( cols.bodyA === ch ){
                            cols.bodyA = compound;
                        } else {
                            cols.bodyB = compound;
                        }
                        ret.push( cols );
                    }
                }

                // transform it back
                ch.state.angular.pos -= compound.state.angular.pos;
                ch.offset.vsub( oldPos );
                ch.state.pos.clone( oldPos.rotate( ch.state.angular.pos ).vsub( compound.offset ) );
            }

            return scratch.done( ret );

        } else {

            return checkGJK( bodyA, bodyB );
        }
    }
822974
private Group getNewGroup(String recordName,Character aminoCode1, long seq_id,String groupCode3) {

		Group g = ChemCompGroupFactory.getGroupFromChemCompDictionary(groupCode3);
		if ( g != null && !g.getChemComp().isEmpty()) {
			if ( g instanceof AminoAcidImpl) {
				AminoAcidImpl aa = (AminoAcidImpl) g;
				aa.setId(seq_id);
			} else if ( g instanceof NucleotideImpl) {
				NucleotideImpl nuc =  (NucleotideImpl) g;
				nuc.setId(seq_id);
			} else if ( g instanceof HetatomImpl) {
				HetatomImpl het = (HetatomImpl)g;
				het.setId(seq_id);
			}
			return g;
		}



		Group group;
		if ( recordName.equals("ATOM") ) {
			if (StructureTools.isNucleotide(groupCode3))  {
				// it is a nucleotide
				NucleotideImpl nu = new NucleotideImpl();
				group = nu;
				nu.setId(seq_id);

			} else if (aminoCode1==null || aminoCode1 == StructureTools.UNKNOWN_GROUP_LABEL){
				HetatomImpl h = new HetatomImpl();
				h.setId(seq_id);
				group = h;

			} else {
				AminoAcidImpl aa = new AminoAcidImpl() ;
				aa.setAminoType(aminoCode1);
				aa.setId(seq_id);
				group = aa ;
			}
		}
		else {
			if (StructureTools.isNucleotide(groupCode3))  {
				// it is a nucleotide
				NucleotideImpl nu = new NucleotideImpl();
				group = nu;
				nu.setId(seq_id);
			}
			else if (aminoCode1 != null ) {
				AminoAcidImpl aa = new AminoAcidImpl() ;
				aa.setAminoType(aminoCode1);
				aa.setId(seq_id);
				group = aa ;
			} else {
				HetatomImpl h = new HetatomImpl();
				h.setId(seq_id);
				group = h;
			}
		}
		return  group ;
	}