Java Code Examples for org.biojava.nbio.structure.align.model.AFPChain#getOptLen()

/**
 * Creates a Map specifying the alignment as a mapping between residue indices
 * of protein 1 and residue indices of protein 2.
 *
 * <p>For example,<pre>
 * 1234
 * 5678</pre>
 * becomes<pre>
 * 1->5
 * 2->6
 * 3->7
 * 4->8</pre>
 *
 * @param afpChain An alignment
 * @return A mapping from aligned residues of protein 1 to their partners in protein 2.
 * @throws StructureException If afpChain is not one-to-one
 */
public static Map<Integer, Integer> alignmentAsMap(AFPChain afpChain) throws StructureException {
	Map<Integer,Integer> map = new HashMap<Integer,Integer>();

	if( afpChain.getAlnLength() < 1 ) {
		return map;
	}
	int[][][] optAln = afpChain.getOptAln();
	int[] optLen = afpChain.getOptLen();
	for(int block = 0; block < afpChain.getBlockNum(); block++) {
		for(int pos = 0; pos < optLen[block]; pos++) {
			int res1 = optAln[block][0][pos];
			int res2 = optAln[block][1][pos];
			if(map.containsKey(res1)) {
				throw new StructureException(String.format("Residue %d aligned to both %d and %d.", res1,map.get(res1),res2));
			}
			map.put(res1,res2);
		}
	}
	return map;
}
 

/**
 * Retrieves the optimum alignment from an AFPChain and returns it as a
 * java collection. The result is indexed in the same way as
 * {@link AFPChain#getOptAln()}, but has the correct size().
 * <pre>
 * List<List<List<Integer>>> aln = getOptAlnAsList(AFPChain afpChain);
 * aln.get(blockNum).get(structureNum={0,1}).get(pos)</pre>
 *
 * @param afpChain
 * @return
 */
public static List<List<List<Integer>>> getOptAlnAsList(AFPChain afpChain) {
	int[][][] optAln = afpChain.getOptAln();
	int[] optLen = afpChain.getOptLen();
	List<List<List<Integer>>> blocks = new ArrayList<List<List<Integer>>>(afpChain.getBlockNum());
	for(int blockNum=0;blockNum<afpChain.getBlockNum();blockNum++) {
		//TODO could improve speed an memory by wrapping the arrays with
		// an unmodifiable list, similar to Arrays.asList(...) but with the
		// correct size parameter.
		List<Integer> align1 = new ArrayList<Integer>(optLen[blockNum]);
		List<Integer> align2 = new ArrayList<Integer>(optLen[blockNum]);
		for(int pos=0;pos<optLen[blockNum];pos++) {
			align1.add(optAln[blockNum][0][pos]);
			align2.add(optAln[blockNum][1][pos]);
		}
		List<List<Integer>> block = new ArrayList<List<Integer>>(2);
		block.add(align1);
		block.add(align2);
		blocks.add(block);
	}

	return blocks;
}
 

/**
 * Print an AFPChain manually for debugging
 * @param cpAlignment
 */
@SuppressWarnings("unused")
private static void printOptAln(AFPChain cpAlignment) {
	int[] optLen = cpAlignment.getOptLen();
	int[][][] optAln = cpAlignment.getOptAln();

	for(int block=0;block<cpAlignment.getBlockNum();block++) {
		for(int pos=0;pos<optLen[block]; pos++) {
			System.out.format("%s\t%s\n", optAln[block][0][pos],
					optAln[block][1][pos]);
		}
		System.out.println();
	}
}
 

private static void printXMLEQRKnownPositions(PrettyXMLWriter xml,
		AFPChain afpChain, int blockNr) throws IOException, StructureException{
	int[] optLen       = afpChain.getOptLen();


	String[][][] pdbAln = afpChain.getPdbAln();
	if ( pdbAln == null){
		throw new StructureException("Can't convert to XML without known the PDB coordinates. Please provide Ca atoms and call toXML(afpChain,ca1, ca2)");
	}

	for ( int eqrNr = 0 ; eqrNr < optLen[blockNr] ; eqrNr++ ){

		String pdbResnum1 = pdbAln[blockNr][0][eqrNr];
		String pdbResnum2 = pdbAln[blockNr][1][eqrNr];

		//System.out.println(eqrNr + " got resnum: " + pdbResnum1 + " " + pdbResnum2);

		String[] spl1 = pdbResnum1.split(":");
		String[] spl2 = pdbResnum2.split(":");

		String chain1 = spl1[0];
		String pdbres1 = spl1[1];

		String chain2 = spl2[0];
		String pdbres2 = spl2[1];

		xml.openTag("eqr");
		xml.attribute("eqrNr", String.valueOf(eqrNr));
		xml.attribute("pdbres1",pdbres1);
		xml.attribute("chain1", chain1);
		xml.attribute("pdbres2",pdbres2);
		xml.attribute("chain2", chain2);
		xml.closeTag("eqr");
	}
}
 

public static void printXMLEQRInferPositions(PrettyXMLWriter xml,
		AFPChain afpChain, int bk, Atom[] ca1, Atom[] ca2)  throws IOException{

	int[] optLen       = afpChain.getOptLen();

	if ( optLen == null)
		return;

	int[][][] optAln   = afpChain.getOptAln();


	for ( int pos=0;pos< optLen[bk];pos++){
		int pos1 = optAln[bk][0][pos];
		int pos2 = optAln[bk][1][pos];
		xml.openTag("eqr");
		xml.attribute("eqrNr", String.valueOf(pos));
		xml.attribute("pdbres1",ca1[pos1].getGroup().getResidueNumber().toString());
		xml.attribute("chain1", ca1[pos1].getGroup().getChain().getName());
		xml.attribute("pdbres2",ca2[pos2].getGroup().getResidueNumber().toString());
		xml.attribute("chain2", ca2[pos2].getGroup().getChain().getName());

		xml.closeTag("eqr");
		//System.out.println("aligned position: " + pos1  + ":" + pos2 +
		//" pdbresnum " + ca1[pos1].getGroup().getResidueNumber().toString() + " " +
		//ca1[pos1].getParent().getPDBName()+":" +
		//ca2[pos2].getGroup().getResidueNumber().toString() + " " + ca2[pos2].getParent().getPDBName());

	}

}
 

/**
 * Return a list of pdb Strings corresponding to the aligned positions of the molecule.
 * Only supports a pairwise alignment with the AFPChain DS.
 *
 * @param aligPos
 * @param afpChain
 * @param ca
 */
public static final List<String> getPDBresnum(int aligPos, AFPChain afpChain, Atom[] ca){
	List<String> lst = new ArrayList<String>();
	if ( aligPos > 1) {
		System.err.println("multiple alignments not supported yet!");
		return lst;
	}

	int blockNum = afpChain.getBlockNum();
	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	if ( optLen == null)
		return lst;

	for(int bk = 0; bk < blockNum; bk ++)       {

		for ( int i=0;i< optLen[bk];i++){

			int pos = optAln[bk][aligPos][i];
			if ( pos < ca.length) {
				String pdbInfo = JmolTools.getPdbInfo(ca[pos]);
				//lst.add(ca1[pos].getParent().getPDBCode());
				lst.add(pdbInfo);
			}
		}

	}
	return lst;
}
 

private static void reverseOptAln(AFPChain afpChain) {
	int ca2len = afpChain.getCa2Length();

	int[][][] optAln = afpChain.getOptAln();
	int[] optLen = afpChain.getOptLen();

	for (int block = 0; block < afpChain.getBlockNum(); block++) {
		for (int pos = 0; pos < optLen[block]; pos++) {
			optAln[block][1][pos] = ca2len - 1 - optAln[block][1][pos];
		}
	}

	afpChain.setOptAln(optAln);
}
 

/**
 * Fill the aligned Atom arrays with the equivalent residues in the afpChain.
 * @param afpChain
 * @param ca1
 * @param ca2
 * @param ca1aligned
 * @param ca2aligned
 */
public static void fillAlignedAtomArrays(AFPChain afpChain, Atom[] ca1,
		Atom[] ca2, Atom[] ca1aligned, Atom[] ca2aligned) {

	int pos=0;
	int[] blockLens = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();
	assert(afpChain.getBlockNum() <= optAln.length);

	for (int block=0; block < afpChain.getBlockNum(); block++) {
		for(int i=0;i<blockLens[block];i++) {
			int pos1 = optAln[block][0][i];
			int pos2 = optAln[block][1][i];
			Atom a1 = ca1[pos1];
			Atom a2 = (Atom) ca2[pos2].clone();
			ca1aligned[pos] = a1;
			ca2aligned[pos] = a2;
			pos++;
		}
	}

	// this can happen when we load an old XML serialization which did not support modern ChemComp representation of modified residues.
	if (pos != afpChain.getOptLength()){
		logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + afpChain.getOptLength());
		// we need to resize the array, because we allocated too many atoms earlier on.
		ca1aligned = (Atom[]) resizeArray(ca1aligned, pos);
		ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
	}

}
 

private static final int getUngappedFatCatPos(AFPChain afpChain, int chainNr, int aligPos){
	char[] aseq;
	if ( chainNr == 0 )
		aseq = afpChain.getAlnseq1();
	else
		aseq = afpChain.getAlnseq2();

	if ( aligPos > aseq.length)
		return -1;
	if ( aligPos < 0)
		return -1;

	int blockNum = afpChain.getBlockNum();
	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	int p1, p2;
	int p1b = 0;
	int p2b = 0;
	int len = 0;


	for(int i = 0; i < blockNum; i ++)  {
		for(int j = 0; j < optLen[i]; j ++) {

			p1 = optAln[i][0][j];
			p2 = optAln[i][1][j];


			if(len > 0)     {

				int lmax = (p1 - p1b - 1)>(p2 - p2b - 1)?(p1 - p1b - 1):(p2 - p2b - 1);

				// lmax gives the length of an alignment gap

				//System.out.println("   p1-p2: " + p1 + " - " + p2 + " lmax: " + lmax + " p1b-p2b:"+p1b + " " + p2b);
				for(int k = 0; k < lmax; k ++)      {

					if(k >= (p1 - p1b - 1)) {
						// a gap position in chain 0
						if ( aligPos == len && chainNr == 0){
							return -1;
						}
					}
					else {
						if ( aligPos == len && chainNr == 0)
							return p1b+1+k;


					}
					if(k >= (p2 - p2b - 1)) {
						// a gap position in chain 1
						if ( aligPos == len && chainNr == 1){
							return -1;
						}
					}
					else  {
						if ( aligPos == len && chainNr == 1) {
							return p2b+1+k;
						}


					}
					len++;

				}
			}

			if ( aligPos == len && chainNr == 0)
				return p1;
			if ( aligPos == len && chainNr == 1)
				return p2;

			len++;
			p1b = p1;
			p2b = p2;


		}
	}


	// we did not find an aligned position
	return -1;
}
 

public static boolean[][] blankOutBreakFlag(AFPChain afpChain, Atom[] ca2,
		int rows, int cols, CECalculator calculator, boolean[][] breakFlag,
		int blankWindowSize) {

	int[][][] optAln = afpChain.getOptAln();
	int blockNum = afpChain.getBlockNum();

	int[] optLen = afpChain.getOptLen();

	// ca2 is circularly permutated at this point.
	int breakPoint = ca2.length / 2;

	for (int bk = 0; bk < blockNum; bk++) {

		// Matrix m= afpChain.getBlockRotationMatrix()[bk];
		// Atom shift = afpChain.getBlockShiftVector()[bk];
		for (int i = 0; i < optLen[bk]; i++) {
			int pos1 = optAln[bk][0][i];
			int pos2 = optAln[bk][1][i];
			// blank out area around these positions...

			int dist = blankWindowSize;
			int start1 = Math.max(pos1 - dist, 0);
			int start2 = Math.max(pos2 - dist, 0);
			int end1 = Math.min(pos1 + dist, rows - 1);
			int end2 = Math.min(pos2 + dist, cols - 1);

			for (int i1 = start1; i1 < end1; i1++) {

				for (int j2 = start2; j2 < end2; j2++) {
					// System.out.println(i1 + " " + j2 + " (***)");
					breakFlag[i1][j2] = true;
					if (j2 < breakPoint) {
						breakFlag[i1][j2 + breakPoint] = true;
					}
				}
			}

		}
	}

	return breakFlag;
}
 

public static void rotateAtoms2(AFPChain afpChain,Atom[] ca2){



		int blockNum = afpChain.getBlockNum();

		int[] optLen = afpChain.getOptLen();
		int[][][] optAln = afpChain.getOptAln();

		for(int bk = 0; bk < blockNum; bk ++)       {

			Matrix m= afpChain.getBlockRotationMatrix()[bk];
			Atom shift = afpChain.getBlockShiftVector()[bk];
			for ( int i=0;i< optLen[bk];i++){
				int pos = optAln[bk][1][i];
				Atom a = ca2[pos];

				Calc.rotate(a, m);
				Calc.shift(a, shift);

				//atoms.add(ca2[pos]);
			}

		}

	}
 

public static String getJmolScript4Block(AFPChain afpChain, Atom[] ca1, Atom[] ca2, int blockNr) {
	int blockNum = afpChain.getBlockNum();

	if (blockNr >= blockNum)
		return DEFAULT_SCRIPT;

	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	if (optLen == null)
		return DEFAULT_SCRIPT;

	StringWriter jmol = new StringWriter();
	jmol.append(DEFAULT_SCRIPT);

	jmol.append("select */2; color lightgrey; model 2; ");

	printJmolScript4Block(ca1, ca2, blockNum, optLen, optAln, jmol, blockNr);

	jmol.append("model 0;  ");
	jmol.append(LIGAND_DISPLAY_SCRIPT);
	// System.out.println(jmol);
	return jmol.toString();

}
 

@Test
	public void testFilterDuplicateAFPs() throws Exception {
		int[][][] dupAlign = new int[1][2][];

		int ca2len = 12;
		dupAlign[0][0] = new int[] { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13 };
		dupAlign[0][1] = new int[] { 3, 5, 6, 7, 8, 9,10,11, 0+ca2len, 1+ca2len, 2+ca2len, 3+ca2len, 4+ca2len, 7+ca2len };
		Atom[] ca1,ca2;

		ca1 = makeDummyCA(dupAlign[0][0].length);
		ca2 = makeDummyCA(ca2len);
		ca2 = StructureTools.duplicateCA2(ca2);
		AFPChain afp = makeDummyAFPChain(dupAlign, ca1, ca2);

		CECPParameters params = new CECPParameters();
		params.setMinCPLength(0);
//		AFPChain newAFP = (AFPChain) filterDuplicateAFPs.invoke(null, afp, new CECalculator(null), ca1, ca2);
		AFPChain newAFP = CeCPMain.filterDuplicateAFPs(afp, new CECalculator(null), ca1, ca2, params);

		int[][][] align = newAFP.getOptAln();
		int[] blkLen = newAFP.getOptLen();
		// optimal alignment should be
		//  1  2  3  4  5  6  7 | 8  9 10 11 12
		//  5  6  7  8  9 10 11 | 0  1  2  3  4
		int[][][] expected = new int[][][] {
				new int[][] {
						new int[] { 1,  2,  3,  4,  5,  6,  7, },
						new int[] { 5,  6,  7,  8,  9, 10, 11, },
				},
				new int[][] {
						new int[] { 8,  9, 10, 11, 12, },
						new int[] { 0,  1,  2,  3,  4, },
				},
		};

		int[] expectedLen = new int[] { expected[0][0].length, expected[1][0].length };

		Assert.assertTrue(Arrays.deepEquals(expected, align));
		Assert.assertTrue(Arrays.equals(expectedLen, blkLen));

	}
 

/**
 * superimposing according to the optimized alignment
 *
 * @param afpChain
 * @param ca1
 * @param ca2
 * @return Group array twisted.
 * @throws StructureException
 */
public static Group[] twistOptimized(AFPChain afpChain, Atom[] ca1,
		Atom[] ca2) throws StructureException {

	Atom[] optTwistPdb = new Atom[ca2.length];

	int gPos = -1;
	for (Atom a : ca2) {
		gPos++;
		optTwistPdb[gPos] = a;
	}

	int blockNum = afpChain.getBlockNum();

	int b2 = 0;
	int e2 = 0;
	int focusResn = 0;
	int[] focusRes1 = afpChain.getFocusRes1();
	int[] focusRes2 = afpChain.getFocusRes2();

	if (focusRes1 == null) {
		focusRes1 = new int[afpChain.getCa1Length()];
		afpChain.setFocusRes1(focusRes1);
	}
	if (focusRes2 == null) {
		focusRes2 = new int[afpChain.getCa2Length()];
		afpChain.setFocusRes2(focusRes2);
	}

	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	for (int bk = 0; bk < blockNum; bk++) {
		// THIS IS TRANSFORMING THE ORIGINAL ca2 COORDINATES, NO CLONING...
		// copies the atoms over to iniTwistPdb later on in modifyCod
		transformOrigPDB(optLen[bk], optAln[bk][0], optAln[bk][1], ca1,
				ca2, afpChain, bk);

		// transform pro2 according to comparison of pro1 and pro2 at give
		// residues
		if (bk > 0) {
			b2 = e2;
		}
		if (bk < blockNum - 1) { // bend at the middle of two consecutive
									// blocks
			e2 = optAln[bk][1][optLen[bk] - 1];
			e2 = (optAln[bk + 1][1][0] - e2) / 2 + e2;
		} else {
			e2 = ca2.length;
		}
		cloneAtomRange(optTwistPdb, ca2, b2, e2);
		for (int i = 0; i < optLen[bk]; i++) {
			focusRes1[focusResn] = optAln[bk][0][i];
			focusRes2[focusResn] = optAln[bk][1][i];
			focusResn++;
		}
	}
	int totalLenOpt = focusResn;
	logger.debug("calrmsdopt for {} residues", focusResn);
	double totalRmsdOpt = calCaRmsd(ca1, optTwistPdb, focusResn, focusRes1,
			focusRes2);
	logger.debug("got opt RMSD: {}", totalRmsdOpt);
	int optLength = afpChain.getOptLength();

	if (totalLenOpt != optLength) {
		logger.warn("Final alignment length is different {} {}",
				totalLenOpt, optLength);
	}
	logger.debug("final alignment length {}, rmsd {}", focusResn,
			totalRmsdOpt);

	afpChain.setTotalLenOpt(totalLenOpt);
	afpChain.setTotalRmsdOpt(totalRmsdOpt);

	return StructureTools.cloneGroups(optTwistPdb);

}
 

public  static double getTMScore(AFPChain align, Atom[] ca1, Atom[] ca2, boolean normalizeMin) throws StructureException
{
	if ( align.getNrEQR() == 0)
		return -1;


	// Create new arrays for the subset of atoms in the alignment.
	Atom[] ca1aligned = new Atom[align.getOptLength()];
	Atom[] ca2aligned = new Atom[align.getOptLength()];
	int pos=0;
	int[] blockLens = align.getOptLen();
	int[][][] optAln = align.getOptAln();
	assert(align.getBlockNum() <= optAln.length);

	for(int block=0;block< align.getBlockNum();block++) {

		if ( ! ( blockLens[block] <= optAln[block][0].length)) {
			logger.warn("AFPChainScorer getTMScore: errors reconstructing alignment block [" + block + "]. Length is " + blockLens[block] + " but should be <=" + optAln[block][0].length);
		}

		for(int i=0;i<blockLens[block];i++) {
			int pos1 = optAln[block][0][i];
			int pos2 = optAln[block][1][i];
			Atom a1 = ca1[pos1];
			Atom a2 = (Atom) ca2[pos2].clone();

			ca1aligned[pos] = a1;
			ca2aligned[pos] = a2;
			pos++;
		}
	}

	// this can happen when we load an old XML serialization which did not support modern ChemComp representation of modified residues.
	if ( pos != align.getOptLength()){
		logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + align.getOptLength());
		// we need to resize the array, because we allocated too many atoms earlier on.
		ca1aligned = (Atom[]) resizeArray(ca1aligned, pos);
		ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
	}
	//Superimpose
	Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(ca1aligned),
			Calc.atomsToPoints(ca2aligned));

	Calc.transform(ca2aligned, trans);

	return Calc.getTMScore(ca1aligned, ca2aligned, ca1.length, ca2.length, normalizeMin);
}
 

/** get the block number for an aligned position
 *
 * @param afpChain
 * @param aligPos
 * @return
 */
public static int getBlockNrForAlignPos(AFPChain afpChain, int aligPos){
	// moved here from DisplayAFP;

	int blockNum = afpChain.getBlockNum();

	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();

	int len = 0;
	int p1b=0;
	int p2b=0;

	for(int i = 0; i < blockNum; i ++)  {

		for(int j = 0; j < optLen[i]; j ++) {

			int p1 = optAln[i][0][j];
			int p2 = optAln[i][1][j];

			if (len != 0) {
				// check for gapped region
				int lmax = (p1 - p1b - 1)>(p2 - p2b - 1)?(p1 - p1b - 1):(p2 - p2b - 1);
				for(int k = 0; k < lmax; k ++)      {
					len++;
				}
			}


			p1b = p1;
			p2b = p2;
			if ( len >= aligPos) {

				return i;
			}
			len++;
		}
	}

	return blockNum;

}
 

public static Atom[] getAlignedAtoms2(AFPChain afpChain,Atom[] ca2){

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

		int blockNum = afpChain.getBlockNum();

		int[] optLen = afpChain.getOptLen();
		int[][][] optAln = afpChain.getOptAln();


		for(int bk = 0; bk < blockNum; bk ++)       {


			for ( int i=0;i< optLen[bk];i++){
				int pos = optAln[bk][1][i];
				atoms.add(ca2[pos]);
			}

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

public static Atom[] getAlignedAtoms1(AFPChain afpChain,Atom[] ca1){
	List<Atom> atoms = new ArrayList<Atom>();

	int blockNum = afpChain.getBlockNum();

	int[] optLen = afpChain.getOptLen();
	int[][][] optAln = afpChain.getOptAln();


	for(int bk = 0; bk < blockNum; bk ++)       {


		for ( int i=0;i< optLen[bk];i++){
			int pos = optAln[bk][0][i];
			atoms.add(ca1[pos]);
		}

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

/**
 * Modifies the {@link AFPChain#setOptAln(int[][][]) optAln} of an AFPChain
 * by permuting the second protein.
 *
 * Sets residue numbers in the second protein to <i>(i-cp)%len</i>
 *
 * @param afpChain
 * @param cp Amount leftwards (or rightward, if negative) to shift the
 */
private static void permuteOptAln(AFPChain afpChain, int cp)
{
	int ca2len = afpChain.getCa2Length();

	if( ca2len <= 0) {
		throw new IllegalArgumentException("No Ca2Length specified in "+afpChain);
	}

	// Allow negative cp points for convenience.
	if(cp == 0) {
		return;
	}
	if(cp <= -ca2len || cp >= ca2len) {
		// could just take cp%ca2len, but probably its a bug if abs(cp)>=ca2len
		throw new ArrayIndexOutOfBoundsException( String.format(
				"Permutation point %d must be between %d and %d for %s",
				cp, 1-ca2len,ca2len-1, afpChain.getName2() ) );
	}
	if(cp < 0) {
		cp = cp + ca2len;
	}

	// the unprocessed alignment
	int[][][] optAln = afpChain.getOptAln();
	int[] optLen = afpChain.getOptLen();

	// the processed alignment
	List<List<List<Integer>>> blocks = new ArrayList<List<List<Integer>>>(afpChain.getBlockNum()*2);

	//Update residue indices
	// newi = (oldi-cp) % N
	for(int block = 0; block < afpChain.getBlockNum(); block++) {
		if(optLen[block]<1)
			continue;

		// set up storage for the current block
		List<List<Integer>> currBlock = new ArrayList<List<Integer>>(2);
		currBlock.add( new ArrayList<Integer>());
		currBlock.add( new ArrayList<Integer>());
		blocks.add(currBlock);

		// pos = 0 case
		currBlock.get(0).add( optAln[block][0][0] );
		currBlock.get(1).add( (optAln[block][1][0]+cp ) % ca2len);

		for(int pos = 1; pos < optLen[block]; pos++) {
			//check if we need to start a new block
			//this happens when the new alignment crosses the protein terminus
			if( optAln[block][1][pos-1]+cp<ca2len &&
					optAln[block][1][pos]+cp >= ca2len) {
				currBlock = new ArrayList<List<Integer>>(2);
				currBlock.add( new ArrayList<Integer>());
				currBlock.add( new ArrayList<Integer>());
				blocks.add(currBlock);
			}
			currBlock.get(0).add( optAln[block][0][pos] );
			currBlock.get(1).add( (optAln[block][1][pos]+cp ) % ca2len);
		}
	}

	// save permuted blocks to afpChain
	assignOptAln(afpChain,blocks);
}
 

public static void rotateAtoms2(AFPChain afpChain,Atom[] ca2){



		int blockNum = afpChain.getBlockNum();

		int[] optLen = afpChain.getOptLen();
		int[][][] optAln = afpChain.getOptAln();

		for(int bk = 0; bk < blockNum; bk ++)       {

			Matrix m= afpChain.getBlockRotationMatrix()[bk];
			Atom shift = afpChain.getBlockShiftVector()[bk];
			for ( int i=0;i< optLen[bk];i++){
				int pos = optAln[bk][1][i];
				Atom a = ca2[pos];

				Calc.rotate(a, m);
				Calc.shift(a, shift);

				//atoms.add(ca2[pos]);
			}

		}

	}