001/* This class is based on the original FATCAT implementation by
002 * <pre>
003 * Yuzhen Ye & Adam Godzik (2003)
004 * Flexible structure alignment by chaining aligned fragment pairs allowing twists.
005 * Bioinformatics vol.19 suppl. 2. ii246-ii255.
006 * https://www.ncbi.nlm.nih.gov/pubmed/14534198
007 * </pre>
008 *
009 * Thanks to Yuzhen Ye and A. Godzik for granting permission to freely use and redistribute this code.
010 *
011 * This code may be freely distributed and modified under the
012 * terms of the GNU Lesser General Public Licence.  This should
013 * be distributed with the code.  If you do not have a copy,
014 * see:
015 *
016 *      http://www.gnu.org/copyleft/lesser.html
017 *
018 * Copyright for this code is held jointly by the individual
019 * authors.  These should be listed in @author doc comments.
020 *
021 *
022 * Created on Jun 17, 2009
023 * Created by Andreas Prlic - RCSB PDB
024 *
025 */
026
027package org.biojava.nbio.structure.align.fatcat.calc;
028
029import org.biojava.nbio.structure.Atom;
030import org.biojava.nbio.structure.Calc;
031import org.biojava.nbio.structure.StructureException;
032import org.biojava.nbio.structure.align.AFPTwister;
033import org.biojava.nbio.structure.align.model.AFP;
034import org.biojava.nbio.structure.align.model.AFPChain;
035import org.biojava.nbio.structure.geometry.SuperPositions;
036import org.biojava.nbio.structure.jama.Matrix;
037
038import java.util.List;
039
040import javax.vecmath.Matrix4d;
041
042/** a class to chain AFPs to an alignment
043 *
044 * @author Andreas Prlic
045 *
046 */
047public class AFPChainer
048{
049        public static final boolean debug = FatCatAligner.debug;
050        // private static final boolean showAlig = false;
051
052        /**
053        // Key function: chain (assembly) the AFPs
054        // a AFP (k) is defined as (i, j, k), with i and j are staring points
055        // AFP extension (eg. AFP(k-1) -> AFP(k) ) requiring
056        // AFP(k-1) < AFP(k)(refer AFP.h definition),
057        // ie i(k-1) < i(k) and j(k-1) < j(k)
058        // in the figure, only (2) AFP can extend to AFP(k)
059        // Key features: a coordination transformation is allowed in the AFP extension
060        //                      gap penalties are also considered
061        //
062        //                                   protein1
063        //                 ---------------------------
064        //                 |        \                |
065        //                 |         \(1)            |
066        //                 |     \    \              |
067        //                 |      \(2) \             |
068        //              p  |       \                 |
069        //              r  |   \                     |
070        //              o  |    \(3)  \(i,j, k)      |
071        //              t  |     \     \             |
072        //              e  |            \            |
073        //              i  |                         |
074        //              n  |                         |
075        //              2  ---------------------------
076        //                  schematic of AFP chaining
077         */
078        public static void doChainAfp(FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[] ca2){
079                List<AFP> afpSet = afpChain.getAfpSet();
080
081                afpChain.setConn(0d);
082                afpChain.setDVar(0d);
083
084                int afpNum = afpSet.size();
085                if(afpNum <= 0) return;
086
087
088                int[] twi = afpChain.getTwi();
089                if(twi == null) {
090                        twi = new int[afpNum];
091                        afpChain.setTwi(twi);
092                }
093
094                //transformation, calculated at DoChainAfp, be used in List extraction
095
096                //forward: calculate the score matrix
097                boolean isConnected = false;
098                int     i, j, j0,  n;
099                double  stmp;
100
101                afpChain.setConn(0d);
102                afpChain.setDVar(0d);
103
104                double[]  sco = new double[afpNum]; //the score ending at an AFP
105                int[] pre = new int[afpNum];    //the previous AFP
106                double  maxsco = 0;
107                int     maxafp = 0;
108                int[] list = new int[afpNum];
109
110                int maxGap = params.getMaxGap();
111                int fragLen = params.getFragLen();
112                int maxTra = params.getMaxTra();
113
114
115                Matrix disTable1 = getDisTable(maxGap + 2 * fragLen + 1,ca1);
116                Matrix disTable2 = getDisTable(maxGap + 2 * fragLen + 1,ca2);
117
118                afpChain.setDisTable1(disTable1);
119                afpChain.setDisTable2(disTable2);
120
121                for(i = 0; i < afpNum; i ++)    {
122                        sco[i] = afpSet.get(i).getScore(); //start from itself
123                        pre[i] = -1;
124                        twi[i] = 0;
125                        if ( afpSet.get(i).getP1() < fragLen || afpSet.get(i).getP2() < fragLen)
126                                n = 0;
127                        else
128                                n = getCompatibleAfps(i, list, params, afpChain); //get a compatible list
129                        //printf("afp %d, compatible %d\n", i, n);
130                        for(j0 = 0; j0 < n; j0 ++)      {
131                                j = list[j0];
132                                isConnected = afpPairConn(j, i, params,afpChain); //note: j, i
133                                Double conn = afpChain.getConn();
134                                int t = 0;
135                                if ( isConnected)
136                                        t=1;
137                                if(twi[j] + t > maxTra) continue;
138                                //two many transformation are disfavored
139                                stmp = sco[j] + afpSet.get(i).getScore() + conn;
140                                if(stmp > sco[i])       { //considered all previous compatible AFPs
141                                        sco[i] = stmp;
142                                        twi[i] = twi[j] + t;
143                                        pre[i] = j;
144                                }
145                        }
146                        if(maxsco < sco[i])     {
147                                maxsco = sco[i];
148                                maxafp = i;
149                        }
150                }
151
152                int     currafp = maxafp;
153                if(debug)
154                        System.out.println(String.format("maximum score %f, %d\n", maxsco, twi[currafp]));
155
156                //trace-back from maxafp (maxsco)
157
158                afpChain.setAlignScore(maxsco);
159                afpChain.setAlignScoreUpdate(maxsco);
160                afpChain.setAfpChainTwiNum(0);
161
162                traceBack(pre, currafp, twi[currafp],params,afpChain,ca1,ca2);
163
164
165        }
166
167        private static Matrix getDisTable(int maxlen, Atom[]ca)
168
169        {
170                int length = ca.length;
171                Matrix dis = new Matrix(length,length);
172
173                int     i, j;
174                for(i = 0; i < length; i ++)    {
175                        dis.set(i,i,0);
176                        for(j = i + 1;( j < length) && (j <= i + maxlen); j ++)     {
177                                dis.set(i,j,0);
178
179                                double val = dis.get(i,j) + (Calc.getDistance(ca[i],ca[j])) * (Calc.getDistance(ca[i],ca[j]));
180                                dis.set(i,j,val);
181
182
183                                dis.set(i,j,Math.sqrt(dis.get(i,j)));
184                                dis.set(j,i,dis.get(i,j));
185                        }
186                }
187                return dis;
188
189        }
190
191        /*
192
193        derive the compabitle AFP lists for AFP-chaining
194        this is important for speeding up the process
195        for a given AFP(i1,j1), there are three regions that could be the starting
196        point for the compabitle AFPs of AFP(i1,j1)
197        //                 a1        a2   a3
198        //               i1-G    i1-f-c i1-f+1 i1
199        //                 |          |   |   |
200        //              ----------------------------
201        //              | B               |   |
202        //b1  j1-G  -|  ---------------|   |
203        //              |  |          |   |   |
204        //              |  |     C    | 3 |   |
205        //              |  |          |   |   |
206        //b2 j1-f-c -|  |--------------|   |
207        //              |  |     2    | 1 |   |
208        //b3 j1-f+1 -|------------------   |
209        //              |                   A |
210        //          j1 -|---------------------\
211        //              |                      \ (AFP(i1,j1))
212        //              -----------------------------
213        //
214        f: the length of AFPs (we use segments of same length)
215        G: g + f, where g is the maximum allowed gaps
216        c: the maximum allowed cross-over in AFP-connection,
217                 here we use c = f, and j1-f-c = j1-2f
218        incompatible region A: its AFPs overlap with given AFP(i1,j1)
219        incompatible region B: the gaps between its AFP with given AFP is larger than g
220        incompatible region C: its AFPs connect with given AFP but cross a given threshold.
221        compatible region 1: [i1-f-c,i1-f+1>,[j1-f-c,j1-f+1> or [a2,a3],[b2,b3]
222        compatible region 2: [i1-G,i1-f-c],[j1-f-c,j1-f] or [a1,a2],[b2,b3]
223        combine 1 and 2    : [i1-G,i1-f],[j1-f-c,j1-f]   or [a1,a3],[b2,b3]
224        compatible region 3: [i1-f-c,i1-f],[j1-G,j1-f-c] or [a2,a3],[b1,b2]
225        c->misCut
226        f->fragLen
227        G->fragLen+maxGap->maxGapFrag
228         *
229         *
230         */
231        private  static int getCompatibleAfps(int afp, int[] list, FatCatParameters params, AFPChain afpChain){
232
233                int     i, j, i1, j1, f, G, c, a1, a2, a3, b1, b2, b3, s1, s2;
234
235                int fragLen = params.getFragLen();
236                int maxGapFrag = params.getMaxGapFrag();
237                int misCut = params.getMisCut();
238                int maxTra = params.getMaxTra();
239                List<AFP> afpSet = afpChain.getAfpSet();
240
241                f = fragLen;
242                G = maxGapFrag;
243                c = misCut;
244
245                i1 = afpSet.get(afp).getP1();
246                j1 = afpSet.get(afp).getP2();
247                a3 = i1 - f;
248                a2 = a3 - c;
249                a1 = i1 - G;
250                a2 = a2>0?a2:0;
251                a1 = a1>0?a1:0;
252
253                b3 = j1 - f;
254                b2 = b3 - c;
255                b1 = j1 - G;
256                b2 = (b2 > 0)?b2:0;
257                b1 = (b1 > 0)?b1:0;
258
259                int[][] afpAftIndex = afpChain.getAfpAftIndex();
260                int[][] afpBefIndex = afpChain.getAfpBefIndex();
261                int[] twi = afpChain.getTwi();
262
263
264
265                int     n = 0;
266                //compatible region 1-2, [a1,a3][b2,b3]
267                                for(i = a1; i <= a3; i ++)      {//i <= a3 instead of i < a3
268                                        s1 = afpAftIndex[i][b2]; //note afpAftIndex, not afpIndex
269                                        if(s1 < 0)      continue;//no AFP for the given i with j > b2
270                                        s2 = afpBefIndex[i][b3]; //afps is sorted by j given a i,it's sparse matrix
271                                        if(s2 < 0)      continue;//no AFP for the given i with j < b3
272                                        for(j = s1; j <= s2; j ++)      { //j <= s2 instead of j < s2
273                                                if(twi[j] <= maxTra)    {
274                                                        list[n ++] = j;
275                                                }
276                                        }
277                                }
278
279                                //compatible region 3  [a2,a3][b1,b2]
280                                for(i = a2; i <= a3; i ++)      {
281                                        s1 = afpAftIndex[i][b1];
282                                        if(s1 < 0)      continue;
283                                        s2 = afpBefIndex[i][b2]; //afps is sorted by j given a i
284                                        if(s2 < 0)      continue;
285                                        //note j < s2, as the cases of j == s2 is alread considered in previous region
286                                        for(j = s1; j < s2; j ++)       {
287                                                if(twi[j] <= maxTra)    {
288                                                        list[n ++] = j;
289                                                }
290                                        }
291                                }
292
293                                return n;
294
295        }
296
297
298        /**
299        //Key function: calculate the connectivity of AFP pairs
300        //no compatibility criteria is executed
301        //note: afp1 is previous to afp2 in terms of the position
302                 //this module must be optimized
303         *
304         * @param afp1
305         * @param afp2
306         * @return flag if they are connected
307         */
308        public static boolean afpPairConn(int afp1, int afp2,  FatCatParameters params, AFPChain afpChain)
309
310        {
311
312                Double conn = afpChain.getConn();
313                Double dvar = afpChain.getDVar();
314
315                double misScore = params.getMisScore();
316                double maxPenalty = params.getMaxPenalty();
317                double disCut = params.getDisCut();
318                double gapExtend = params.getGapExtend();
319                double torsionPenalty = params.getTorsionPenalty();
320                double disSmooth = params.getDisSmooth();
321
322                List<AFP> afpSet = afpChain.getAfpSet();
323
324                int     m = calcGap(afpSet.get(afp2),afpSet.get(afp1));
325                int     g = calcMismatch(afpSet.get(afp2),afpSet.get(afp1));
326
327
328                double  gp = misScore * m;      //on average, penalty for a mismatch is misScore, no modification on score
329                if(g > 0)       {
330                        gp += gapExtend * g;
331                }
332                if(gp < maxPenalty)     gp = maxPenalty; //penalty cut-off
333                //note: use < (smaller) instead of >, because maxPenalty is a negative number
334
335                double  d;
336                d = calAfpDis(afp1, afp2,params, afpChain);
337                //note: the 'dis' value is numerically equivalent to the 'rms' with exceptions
338
339                boolean     ch = false;
340                double  tp = 0.0;
341                if(d >= disCut) {
342                        tp = torsionPenalty;
343                        ch = true;
344                } //use the variation of the distances between AFPs
345                else  if(d > disCut - disSmooth)        {
346                        double  wt = Math.sqrt((d - disCut + disSmooth) / disSmooth);
347                        //using sqrt: penalty increase with dis more quicker than linear function
348                        tp = torsionPenalty * wt;
349                }
350
351                dvar = d;
352                conn = tp + gp;
353
354                afpChain.setConn(conn);
355                afpChain.setDVar(dvar);
356                return ch;
357        }
358
359        /**
360         * return the gaps between this and afp
361         * requiring afp1 >  afp2
362         * ( operator % in C)
363         */
364
365        private static int  calcGap(AFP afp1 , AFP afp2)
366        {
367                int     g = (afp1.getP1() - afp2.getP1()) - (afp1.getP2() - afp2.getP2());
368                if(g < 0)       g = -g;
369                return g;
370        }
371
372        /**
373         * return the mis-matched between this and afp
374         *     requiring  AFP1 > afp2
375         *     (operator / in C)
376         * @param afp1
377         * @param afp2
378         * @return
379         */
380
381        //--------------------------------------------
382        private static int calcMismatch(AFP afp1, AFP afp2)
383        {
384                int     l1 = afp1.getP1() - afp2.getP1() - afp2.getFragLen();
385                int     l2 = afp1.getP2() - afp2.getP2() - afp2.getFragLen();
386                return (l1 > l2?l2:l1);
387        }
388
389        /**
390        //return the root mean square of the distance matrix between the residues
391        //from the segments that form the given AFP list
392        //this value can be a measurement (2) for the connectivity of the AFPs
393        //and its calculation is quicker than the measurement (1), rmsd
394        //currently only deal with AFP pair
395//
396//           |-d1--|
397//          |--d2---|
398//         |---d3----|
399        //-----------------------------------------------------------------------
400        //this module is optimized
401         *
402         * @param afp1
403         * @param afp2
404         * @return
405         */
406        private static double calAfpDis(int afp1, int afp2, FatCatParameters params, AFPChain afpChain)
407        {
408
409                List<AFP> afpSet = afpChain.getAfpSet();
410
411                Matrix disTable1 = afpChain.getDisTable1();
412                Matrix disTable2 = afpChain.getDisTable2();
413
414                int fragLen = params.getFragLen();
415                double afpDisCut = params.getAfpDisCut();
416                double disCut = params.getDisCut();
417                double fragLenSq = params.getFragLenSq();
418
419                int     i, j, ai, bi, aj, bj;
420                double  d;
421                double  rms = 0;
422                for(i = 0; i < fragLen; i ++)   {
423                        ai = afpSet.get(afp1).getP1() + i;
424                        bi = afpSet.get(afp1).getP2() + i;
425                        for(j = 0; j < fragLen; j ++)   {
426                                aj = afpSet.get(afp2).getP1() + j;
427                                bj = afpSet.get(afp2).getP2() + j;
428                                d = disTable1.get(aj,ai) - disTable2.get(bj,bi);
429                                rms += d * d;
430                                if(rms > afpDisCut)     { return (disCut); }
431                        }
432                }
433                return (Math.sqrt(rms / fragLenSq));
434        }
435
436
437        /**
438         * derive the optimal chaining of AFPs by trace-back
439         */
440        private static void traceBack(int[] pre, int currafp0, int twist, FatCatParameters params, AFPChain afpChain,Atom[] ca1, Atom[]ca2)
441        {
442
443                afpChain.setConn(0d);
444                afpChain.setDVar(0d);
445
446                int minLen = afpChain.getMinLen();
447                List<AFP> afpSet = afpChain.getAfpSet();
448
449                int afpChainLen = 0;
450
451                //trace-back from currafp (maxsco)
452                int[]     afpchain    = new int[minLen];
453                int[]     afptwibin   = new int[minLen];
454                double[]  afptwilist  = new double[minLen];
455                int       currafp = currafp0;
456                int       s = 0;
457
458                afptwibin[s] = 0;
459                afpchain[s ++] = currafp;
460
461                boolean isConnected = false;
462                int     prevafp;
463                // int     alnlen = afpSet.get(afpchain[s]).getFragLen();
464
465                //Double  conn = new Double(0) ;
466                //Double dvar =  new Double(0);
467
468                while((prevafp = pre[currafp]) != -1)   {
469
470                        isConnected = afpPairConn(prevafp, currafp, params,afpChain);
471
472                        if ( isConnected )
473                                afptwibin[s - 1] = 1;
474                        else
475                                afptwibin[s - 1] = 0;
476                        Double dvar = afpChain.getDVar();
477                        afptwilist[s - 1] = dvar;
478                        //note s - 1: the transformation of prevafp-currafp is recorded in currafp
479                        currafp = prevafp;
480                        // alnlen += afpSet.get(currafp).getFragLen();
481                        afpchain[s ++] = currafp;
482                }
483
484                afpChainLen = s;
485                afpChain.setAfpChainLen(afpChainLen);
486
487                //first afp without transformation
488                if ( isConnected)
489                        afptwibin[s - 1] = 1;
490                else
491                        afptwibin[s - 1] = 0;
492
493                //if(debug)
494                //   System.out.println(String.format("including %d AFPs, %d residues\n", afpChainLen, alnlen));
495
496                //record the optimal alignment in afpChainList (afpChainLen)
497
498                int[] afpChainList = afpChain.getAfpChainList();
499                double[] afpChainTwiBin = afpChain.getAfpChainTwiBin();
500                double[] afpChainTwiList = afpChain.getAfpChainTwiList();
501
502                if(afpChainList == null)        {
503                        afpChainList   = new int[s];
504                        afpChain.setAfpChainList(afpChainList);
505                        afpChainTwiBin     = new double[s];
506                        afpChain.setAfpChainTwiBin(afpChainTwiBin);
507                        afpChainTwiList = new double[s];
508                        afpChain.setAfpChainTwiList(afpChainTwiList);
509                }
510
511                int afpChainTwiNum = afpChain.getAfpChainTwiNum();
512
513                int     i;
514                for(i = 0; i < s; i ++) {
515                        afpChainList[i]     = afpchain[s - 1 - i];
516                        afpChainTwiBin[i]   = afptwibin[s - 1 - i];
517                        afpChainTwiList[i]  = afptwilist[s - 1 - i];
518                        afpChainTwiNum      += afptwibin[s - 1 - i];
519                }
520
521                if(afpChainTwiNum != twist)     {
522                        System.err.println(String.format("AFPChainer Warning: the twists number is not consistent %d %d\n", afpChainTwiNum, twist));
523                }
524
525                double alignScore = afpChain.getAlignScore();
526
527                double  checkscore = afpSet.get(afpChainList[0]).getScore();
528                for(i = 1; i < afpChainLen; i ++)       {
529                        isConnected = afpPairConn(afpChainList[i - 1], afpChainList[i], params,afpChain);
530                        checkscore = checkscore + afpSet.get(afpChainList[i]).getScore() + afpChain.getConn();
531                }
532                if(Math.abs(checkscore - alignScore) > 1e-4)        {
533                        System.err.println(String.format("AFPChainer Warning: fail in alignment score checking %.4f %.4f\n", alignScore, checkscore));
534                }
535
536                if ( debug )
537                        System.out.println("traceBack:" + afpChainLen + " " + afpChainList.length);
538                double  rmsd = calAfpRmsd(afpChainLen, afpChainList,0, afpChain,ca1,ca2);
539                afpChain.setChainRmsd(rmsd);
540                if (debug )
541                        System.out.println("Chain RMSD: " + rmsd);
542                int     b1 = 0;
543                int     bk = 0;
544                int     a, b;
545                afpChain.setChainLen( 0);
546                int chainLen       = afpChain.getChainLen();
547                int[] block2Afp = afpChain.getBlock2Afp();
548
549
550                double[] blockRmsd = afpChain.getBlockRmsd();
551                int[] blockSize = afpChain.getBlockSize();
552
553                block2Afp[0] = 0;
554                for(i = 0; i < afpChainLen; i ++)       {
555                        a = afpChainList[i];
556                        chainLen += afpSet.get(a).getFragLen();
557                        if(i > 0)       {
558                                b = afpChainList[i - 1];
559                                int misLen = afpChain.getMisLen();
560                                misLen += calcMismatch(afpSet.get(a),afpSet.get(b));
561                                afpChain.setMisLen(misLen);
562                                int gapLen = afpChain.getGapLen();
563                                gapLen += calcGap(afpSet.get(a),afpSet.get(b));
564                                afpChain.setGapLen(gapLen);
565
566                        }
567
568                        if(afpChainTwiBin[i] == 1)   {
569                                if (debug)
570                                        System.out.println(" ** calAfpTmsd : afpChainWtiBin == 1 : i: "+i+" i-b1: " + (i-b1) + " b1: " + b1 + " afpChainList.len: " + afpChainList.length);
571
572                                //int len = afpChainList.length - b1 +1;
573                                //int[] fakeList = new int[len];
574                                //int pos = -1;
575                                //for ( int fPos = b1 ; fPos< afpChainList.length ; fPos++){
576                                //   pos++;
577                                //   fakeList[pos] = afpChainList[fPos];
578                                //}
579                                if (debug )
580                                        System.err.println("calculation calAfpRmsd " + i + " " + b1 + " " );
581
582
583                                rmsd = calAfpRmsd(i - b1, afpChainList,b1, afpChain, ca1, ca2);
584
585
586                                blockRmsd[bk] = rmsd;
587                                blockSize[bk] = i - b1;
588                                b1 = i;
589
590                                //System.out.println("block2Afp.length:"+ block2Afp.length + " " + bk + " " + i + " " + afpChain.getMaxTra() );
591                                block2Afp[++bk] = i; //next-block
592                        }
593                }
594                afpChain.setBlock2Afp(block2Afp);
595                afpChain.setChainLen(chainLen);
596
597                if (debug)
598                        System.out.println("after loop over all afpChainList " + (i-b1) + " " + b1);
599
600                rmsd = calAfpRmsd(i - b1, afpChainList, b1, afpChain,ca1,ca2);
601                if (debug)
602                        System.out.println("*** i:" + i + " b1: " + b1 + " i-b1 " + (i-b1));
603
604
605
606                // argh this is the block RMSD, not the chain RMSD!
607                //afpChain.setChainRmsd(rmsd);
608                //rmsd = calAfpRmsd(i - b1, afpChainList[b1]);
609                blockSize[bk] = i - b1;
610                blockRmsd[bk] = rmsd;
611
612                afpChain.setBlockSize(blockSize);
613                afpChain.setBlockRmsd(blockRmsd);
614                int blockNum = afpChain.getBlockNum();
615                blockNum = ++bk;
616                if ( debug)
617                        System.err.println("AFPChainser setBlockNUm:" + blockNum);
618                afpChain.setBlockNum(blockNum);
619                afpChain.setAfpChainList(afpChainList);
620                afpChain.setAfpChainTwiList(afpChainTwiList);
621
622        }
623
624        /**
625        //return the rmsd of the residues from the segments that form the given AFP list
626        //this value can be a measurement (1) for the connectivity of the AFPs
627         *
628         * @param afpn
629         * @param afpPositions the positions of AFPs to work on.
630         * @param listStart the starting position in the list of AFPs
631         * @param afpChain
632         * @param ca1
633         * @param ca2
634         * @return rmsd
635         */
636
637        protected static double calAfpRmsd(int afpn, int[] afpPositions, int listStart,  AFPChain afpChain,Atom[] ca1, Atom[] ca2)
638        {
639
640
641                if (debug)
642                        System.out.println("XXX calling afp2res "+ afpn + " " + afpPositions.length);
643
644                int focusResn = AFPTwister.afp2Res(afpChain,afpn, afpPositions, listStart);
645
646
647                int[] focusRes1 = afpChain.getFocusRes1();
648                int[] focusRes2 = afpChain.getFocusRes2();
649
650                if (debug)
651                        System.out.println("XXX calculating calAfpRmsd: " + focusResn + " " + focusRes1.length + " " + focusRes2.length + " " + afpChain.getMinLen() + " " );
652                double rmsd = getRmsd(focusResn, focusRes1, focusRes2 , afpChain,ca1,  ca2);
653
654                return rmsd;
655        }
656
657
658
659        /** the the RMSD for the residues defined in the two arrays
660         *
661         * @param focusResn
662         * @param focusRes1
663         * @param focusRes2
664         * @return
665         */
666        private static double getRmsd(int focusResn, int[] focusRes1, int[] focusRes2, AFPChain afpChain, Atom[] ca1, Atom[] ca2){
667
668
669
670                Atom[] tmp1 = new Atom[focusResn];
671                Atom[] tmp2 = new Atom[focusResn];
672
673                for ( int i =0 ; i< focusResn;i++){
674                        tmp1[i] =       ca1[focusRes1[i]];
675                        tmp2[i] = (Atom)ca2[focusRes2[i]].clone();
676                        if (tmp1[i].getCoords() == null){
677                                System.err.println("tmp1 got null: " +i + " pos: " + focusRes1[i]);
678                        }
679                        if (tmp2[i].getCoords() == null){
680                                System.err.println("tmp1 got null: " +i + " pos: " + focusRes2[i]);
681                        }
682                        //XX
683                        //tmp2[i].setParent((Group) ca2[focusRes2[i]].getParent().clone());
684                }
685                double rmsd = 99;
686                try {
687                        rmsd = getRmsd(tmp1,tmp2);
688
689                } catch (Exception e){
690                        e.printStackTrace();
691                }
692                return rmsd;
693
694        }
695        /** Calculate the RMSD for two sets of atoms. Rotates the 2nd atom set so make sure this does not cause problems later
696         *
697         *
698         * @param catmp1
699         * @return
700         */
701        private static double getRmsd(Atom[] catmp1, Atom[] catmp2) throws StructureException{
702
703                Matrix4d trans = SuperPositions.superpose(Calc.atomsToPoints(catmp1), 
704                                Calc.atomsToPoints(catmp2));
705
706                Calc.transform(catmp2, trans);
707
708                //   if ( showAlig) {
709                //      StructureAlignmentJmol jmol = new StructureAlignmentJmol();
710                //      jmol.setTitle("AFPCHainer: getRmsd" + rmsd);
711                //
712                //      Chain c1 = new ChainImpl();
713                //      c1.setName("A");
714                //      for ( Atom a : catmp1){
715                //         c1.addGroup(a.getParent());
716                //      }
717                //
718                //      Chain c2 = new ChainImpl();
719                //      c2.setName("B");
720                //      for ( Atom a : catmp2){
721                //         c2.addGroup(a.getParent());
722                //      }
723                //
724                //      Structure fake = new StructureImpl();
725                //      fake.setPDBCode("AFPCHainer: getRmsd" + rmsd);
726                //      List<Chain> model1 = new ArrayList<Chain>();
727                //      model1.add(c1);
728                //      List<Chain> model2 = new ArrayList<Chain>();
729                //      model2.add(c2);
730                //      fake.addModel(model1);
731                //      fake.addModel(model2);
732                //      fake.setNmr(true);
733                //
734                //      jmol.setStructure(fake);
735                //      jmol.evalString("select *; backbone 0.4; wireframe off; spacefill off; " +
736                //      "select not protein and not solvent; spacefill on;");
737                //      jmol.evalString("select */1 ; color red; model 1; ");
738                //
739                //      // now show both models again.
740                //      jmol.evalString("model 0;");
741                //   }
742
743                return Calc.rmsd(catmp1,catmp2);
744        }
745
746}