001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.structure.xtal;
022
023
024import org.biojava.nbio.structure.*;
025import org.biojava.nbio.structure.contact.AtomContactSet;
026import org.biojava.nbio.structure.contact.StructureInterface;
027import org.biojava.nbio.structure.contact.StructureInterfaceList;
028import org.slf4j.Logger;
029import org.slf4j.LoggerFactory;
030
031import javax.vecmath.Matrix4d;
032import javax.vecmath.Point3i;
033import javax.vecmath.Vector3d;
034import java.util.*;
035
036
037/**
038 * A class containing methods to find interfaces in a given crystallographic Structure by
039 * reconstructing the crystal lattice through application of symmetry operators
040 *
041 * @author Jose Duarte
042 *
043 */
044
045public class CrystalBuilder {
046
047        public static final String NCS_CHAINID_SUFFIX_CHAR = "n";
048
049        // Default number of cell neighbors to try in interface search (in 3 directions of space).
050        // In the search, only bounding box overlaps are tried, thus there's not so much overhead in adding
051        // more cells. We actually tested it and using numCells from 1 to 10 didn't change runtimes at all.
052        // Examples with interfaces in distant neighbor cells:
053        //   2nd neighbors: 3hz3, 1wqj, 2de3, 1jcd
054        //   3rd neighbors: 3bd3, 1men, 2gkp, 1wui
055        //   5th neighbors: 2ahf, 2h2z
056        //   6th neighbors: 1was (in fact interfaces appear only at 5th neighbors for it)
057        // Maybe this could be avoided by previously translating the given molecule to the first cell,
058        // BUT! some bona fide cases exist, e.g. 2d3e: it is properly placed at the origin but the molecule
059        // is enormously long in comparison with the dimensions of the unit cell, some interfaces come at the 7th neighbor.
060        // After a scan of the whole PDB (Oct 2013) using numCells=50, the highest one was 4jgc with
061        // interfaces up to the 11th neighbor. Other high ones (9th neighbors) are 4jbm and 4k3t.
062        // We set the default value to 12 based on that (having not seen any difference in runtime)
063        public static final int DEF_NUM_CELLS = 12;
064
065        /**
066         * Default maximum distance between two chains to be considered an interface.
067         * @see #getUniqueInterfaces(double)
068         */
069        public static final double DEFAULT_INTERFACE_DISTANCE_CUTOFF = 5.5;
070
071        public static final Matrix4d IDENTITY = new Matrix4d(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1);
072
073
074        /**
075         * Whether to consider HETATOMs in contact calculations
076         */
077        private static final boolean INCLUDE_HETATOMS = true;
078
079        private Structure structure;
080        private PDBCrystallographicInfo crystallographicInfo;
081        private int numPolyChainsAu;
082        private int numOperatorsSg;
083        private Map<String,Matrix4d> chainNcsOps = null;
084        private Map<String,String> chainOrigNames = null;
085
086        private static final Logger logger = LoggerFactory.getLogger(CrystalBuilder.class);
087
088        private int numCells;
089
090        private ArrayList<CrystalTransform> visitedCrystalTransforms;
091        private Map<String,Map<Matrix4d,StructureInterface>> visitedNcsChainPairs = null;
092
093        private boolean searchBeyondAU;
094        private Matrix4d[] ops;
095
096        /**
097         * Special constructor for NCS-aware CrystalBuilder.
098         * The output list of interfaces will be pre-clustered by NCS-equivalence.
099         * Run {@link CrystalBuilder#expandNcsOps(Structure, Map, Map)} first to extend the AU
100         * and get the equivalence information.
101         * @param structure
102         *          NCS-extended structure
103         * @param chainOrigNames
104         *          chain names mapped to the original chain names (pre-NCS extension)
105         * @param chainNcsOps
106         *          chain names mapped to the ncs operators that was used to generate them
107         * @since 5.0.0
108         */
109        public CrystalBuilder(Structure structure, Map<String,String> chainOrigNames, Map<String,Matrix4d> chainNcsOps) {
110                this(structure);
111                this.chainOrigNames = chainOrigNames;
112                this.chainNcsOps = chainNcsOps;
113        }
114
115        public CrystalBuilder(Structure structure) {
116                this.structure = structure;
117                this.crystallographicInfo = structure.getCrystallographicInfo();
118                this.numPolyChainsAu = structure.getPolyChains().size();
119
120                this.searchBeyondAU = false;
121                if (structure.isCrystallographic()) {
122
123                        this.searchBeyondAU = true;
124
125                        // we need to check space group not null for the cases where the entry is crystallographic but
126                        // the space group is not a standard one recognized by biojava, e.g. 1mnk (SG: 'I 21')
127                        if (this.crystallographicInfo.isNonStandardSg()) {
128                                logger.warn("Space group is non-standard, will only calculate asymmetric unit interfaces.");
129                                this.searchBeyondAU = false;
130                        }
131
132                        // just in case we still check for space group null (a user pdb file could potentially be crystallographic and have no space group)
133                        if (this.crystallographicInfo.getSpaceGroup() == null) {
134                                logger.warn("Space group is null, will only calculate asymmetric unit interfaces.");
135                                this.searchBeyondAU = false;
136                        }
137
138                        // we need to check crystal cell not null for the rare cases where the entry is crystallographic but
139                        // the crystal cell is not given, e.g. 2i68, 2xkm, 4bpq
140                        if (this.crystallographicInfo.getCrystalCell() == null) {
141                                logger.warn("Could not find a crystal cell definition, will only calculate asymmetric unit interfaces.");
142                                this.searchBeyondAU = false;
143                        }
144
145                        // check for cases like 4hhb that are in a non-standard coordinate frame convention, see https://github.com/eppic-team/owl/issues/4
146                        if (this.crystallographicInfo.isNonStandardCoordFrameConvention()) {
147                                logger.warn("Non-standard coordinate frame convention, will only calculate asymmetric unit interfaces.");
148                                this.searchBeyondAU = false;
149                        }
150                }
151
152                if (this.searchBeyondAU) {
153                        // explore the crystal
154                        this.numOperatorsSg = this.crystallographicInfo.getSpaceGroup().getMultiplicity();
155                        this.ops = this.crystallographicInfo.getTransformationsOrthonormal();
156                } else {
157                        // look for contacts within structure as given
158                        this.numOperatorsSg = 1;
159                        this.ops = new Matrix4d[1];
160                        this.ops[0] = new Matrix4d(IDENTITY);
161                }
162
163                this.numCells = DEF_NUM_CELLS;
164
165        }
166
167
168        /**
169         * @return true if this CrystalBuilder is NCS-aware.
170         * @since 5.0.0
171         */
172        public boolean hasNcsOps() {
173                return chainNcsOps != null;
174        }
175
176        /**
177         * Set the number of neighboring crystal cells that will be used in the search for contacts
178         * @param numCells
179         */
180        public void setNumCells(int numCells) {
181                this.numCells = numCells;
182        }
183
184        private void initialiseVisited() {
185                visitedCrystalTransforms = new ArrayList<>();
186                if(this.hasNcsOps()) {
187                        visitedNcsChainPairs = new HashMap<>();
188                }
189        }
190
191        /**
192         * Returns the list of unique interfaces that the given Structure has upon
193         * generation of all crystal symmetry mates. An interface is defined as any pair of chains
194         * that contact, i.e. for which there is at least a pair of atoms (one from each chain) within
195         * the default cutoff distance.
196         * @return
197         * @see #DEFAULT_INTERFACE_DISTANCE_CUTOFF
198         */
199        public StructureInterfaceList getUniqueInterfaces() {
200                return getUniqueInterfaces(DEFAULT_INTERFACE_DISTANCE_CUTOFF);
201        }
202
203        /**
204         * Returns the list of unique interfaces that the given Structure has upon
205         * generation of all crystal symmetry mates. An interface is defined as any pair of chains
206         * that contact, i.e. for which there is at least a pair of atoms (one from each chain) within
207         * the given cutoff distance.
208         * @param cutoff the distance cutoff for 2 chains to be considered in contact
209         * @return
210         */
211        public StructureInterfaceList getUniqueInterfaces(double cutoff) {
212
213
214                StructureInterfaceList set = new StructureInterfaceList();
215
216                // certain structures in the PDB are not macromolecules (contain no polymeric chains at all), e.g. 1ao2
217                // with the current mmCIF parsing, those will be empty since purely non-polymeric chains are removed
218                // see commit e9562781f23da0ebf3547146a307d7edd5741090
219                if (numPolyChainsAu==0) {
220                        logger.warn("No chains present in the structure! No interfaces will be calculated");
221                        return set;
222                }
223
224                // pass the chainOrigNames map in NCS case so that StructureInterfaceList can deal with original to NCS chain names conversion
225                if (chainOrigNames!=null) {
226                        set.setChainOrigNamesMap(chainOrigNames);
227                }
228
229                // initialising the visited ArrayList for keeping track of symmetry redundancy
230                initialiseVisited();
231
232
233
234                // the isCrystallographic() condition covers 3 cases:
235                // a) entries with expMethod X-RAY/other diffraction and defined crystalCell (most usual case)
236                // b) entries with expMethod null but defined crystalCell (e.g. PDB file with CRYST1 record but no expMethod annotation)
237                // c) entries with expMethod not X-RAY (e.g. NMR) and defined crystalCell (NMR entries do have a dummy CRYST1 record "1 1 1 90 90 90 P1")
238                // d) isCrystallographic will be false if the structure is crystallographic but the space group was not recognized
239
240
241                calcInterfacesCrystal(set, cutoff);
242
243                return set;
244        }
245
246        /**
247         * Calculate interfaces between original asymmetric unit and neighboring
248         * whole unit cells, including the original full unit cell i.e. i=0,j=0,k=0
249         * @param set
250         * @param cutoff
251         */
252        private void calcInterfacesCrystal(StructureInterfaceList set, double cutoff) {
253
254
255                // initialising debugging vars
256                long start = -1;
257                long end = -1;
258                int trialCount = 0;
259                int skippedRedundant = 0;
260                int skippedAUsNoOverlap = 0;
261                int skippedChainsNoOverlap = 0;
262                int skippedSelfEquivalent = 0;
263
264                // The bounding boxes of all AUs of the unit cell
265                UnitCellBoundingBox bbGrid = new UnitCellBoundingBox(numOperatorsSg, numPolyChainsAu);;
266                // we calculate all the bounds of each of the asym units, those will then be reused and translated
267                bbGrid.setBbs(structure, ops, INCLUDE_HETATOMS);
268
269
270                // if not crystallographic there's no search to do in other cells, only chains within "AU" will be checked
271                if (!searchBeyondAU) numCells = 0;
272
273                boolean verbose = logger.isDebugEnabled();
274
275                if (verbose) {
276                        trialCount = 0;
277                        start= System.currentTimeMillis();
278                        int neighbors = (2*numCells+1)*(2*numCells+1)*(2*numCells+1)-1;
279                        int auTrials = (numPolyChainsAu*(numPolyChainsAu-1))/2;
280                        int trials = numPolyChainsAu*numOperatorsSg*numPolyChainsAu*neighbors;
281                        logger.debug("Chain clash trials within original AU: "+auTrials);
282                        logger.debug(
283                                        "Chain clash trials between the original AU and the neighbouring "+neighbors+
284                                        " whole unit cells ("+numCells+" neighbours)" +
285                                        "(2x"+numPolyChainsAu+"chains x "+numOperatorsSg+"AUs x "+neighbors+"cells) : "+trials);
286                        logger.debug("Total trials: "+(auTrials+trials));
287                }
288
289                List<Chain> polyChains = structure.getPolyChains();
290
291                for (int a=-numCells;a<=numCells;a++) {
292                        for (int b=-numCells;b<=numCells;b++) {
293                                for (int c=-numCells;c<=numCells;c++) {
294
295                                        Point3i trans = new Point3i(a,b,c);
296                                        Vector3d transOrth = new Vector3d(a,b,c);
297                                        if (a!=0 || b!=0 || c!=0) {
298                                                // we avoid doing the transformation for 0,0,0 (in case it's not crystallographic)
299                                                this.crystallographicInfo.getCrystalCell().transfToOrthonormal(transOrth);
300                                        }
301
302                                        UnitCellBoundingBox bbGridTrans = bbGrid.getTranslatedBbs(transOrth);
303
304                                        for (int n=0;n<numOperatorsSg;n++) {
305
306                                                // short-cut strategies
307                                                // 1) we skip first of all if the bounding boxes of the AUs don't overlap
308                                                if (!bbGrid.getAuBoundingBox(0).overlaps(bbGridTrans.getAuBoundingBox(n), cutoff)) {
309                                                        skippedAUsNoOverlap++;
310                                                        continue;
311                                                }
312
313                                                // 2) we check if we didn't already see its equivalent symmetry operator partner
314                                                CrystalTransform tt = new CrystalTransform(this.crystallographicInfo.getSpaceGroup(), n);
315                                                tt.translate(trans);
316                                                if (isRedundantTransform(tt)) {
317                                                        skippedRedundant++;
318                                                        continue;
319                                                }
320                                                addVisitedTransform(tt);
321
322
323                                                boolean selfEquivalent = false;
324
325                                                // 3) an operator can be "self redundant" if it is the inverse of itself (involutory, e.g. all pure 2-folds with no translation)
326                                                if (tt.isEquivalent(tt)) {
327                                                        logger.debug("Transform {} is equivalent to itself, will skip half of i-chains to j-chains comparisons", tt.toString());
328                                                        // in this case we can't skip the operator, but we can skip half of the matrix comparisons e.g. j>i
329                                                        // we set a flag and do that within the loop below
330                                                        selfEquivalent = true;
331                                                }
332
333                                                StringBuilder builder = null;
334                                                if (verbose) builder = new StringBuilder(String.valueOf(tt)).append(" ");
335
336                                                // Now that we know that boxes overlap and operator is not redundant, we have to go to the details
337                                                int contactsFound = 0;
338
339                                                for (int j=0;j<numPolyChainsAu;j++) {
340
341                                                        for (int i=0;i<numPolyChainsAu;i++) { // we only have to compare the original asymmetric unit to every full cell around
342
343                                                                if(selfEquivalent && (j>i)) {
344                                                                        // in case of self equivalency of the operator we can safely skip half of the matrix
345                                                                        skippedSelfEquivalent++;
346                                                                        continue;
347                                                                }
348                                                                // special case of original AU, we don't compare a chain to itself
349                                                                if (n==0 && a==0 && b==0 && c==0 && i==j) continue;
350
351                                                                // before calculating the AtomContactSet we check for overlap, then we save putting atoms into the grid
352                                                                if (!bbGrid.getChainBoundingBox(0,i).overlaps(bbGridTrans.getChainBoundingBox(n,j),cutoff)) {
353                                                                        skippedChainsNoOverlap++;
354                                                                        if (verbose) {
355                                                                                builder.append(".");
356                                                                        }
357                                                                        continue;
358                                                                }
359
360                                                                trialCount++;
361
362                                                                // finally we've gone through all short-cuts and the 2 chains seem to be close enough:
363                                                                // we do the calculation of contacts
364                                                                Chain chaini = polyChains.get(i);
365                                                                Chain chainj = polyChains.get(j);
366
367                                                                if (n!=0 || a!=0 || b!=0 || c!=0) {
368                                                                        Matrix4d mJCryst = new Matrix4d(ops[n]);
369                                                                        translate(mJCryst, transOrth);
370                                                                        chainj = (Chain)chainj.clone();
371                                                                        Calc.transform(chainj,mJCryst);
372                                                                }
373
374                                                                StructureInterface interf = calcContacts(chaini, chainj, cutoff, tt, builder);
375                                                                if (interf == null) {
376                                                                        continue;
377                                                                }
378
379                                                                contactsFound++;
380                                                                if(this.hasNcsOps()) {
381                                                                        StructureInterface interfNcsRef = findNcsRef(interf);
382                                                                        set.addNcsEquivalent(interf,interfNcsRef);
383                                                                } else {
384                                                                        set.add(interf);
385                                                                }
386                                                        }
387                                                }
388
389                                                if( verbose ) {
390                                                        if (a==0 && b==0 && c==0 && n==0)
391                                                                builder.append(" "+contactsFound+"("+(numPolyChainsAu*(numPolyChainsAu-1))/2+")");
392                                                        else if (selfEquivalent)
393                                                                builder.append(" "+contactsFound+"("+(numPolyChainsAu*(numPolyChainsAu+1))/2+")");
394                                                        else
395                                                                builder.append(" "+contactsFound+"("+numPolyChainsAu*numPolyChainsAu+")");
396
397                                                        logger.debug(builder.toString());
398                                                }
399                                        }
400                                }
401                        }
402                }
403
404                end = System.currentTimeMillis();
405                logger.debug("\n"+trialCount+" chain-chain clash trials done. Time "+(end-start)/1000+"s");
406                logger.debug("  skipped (not overlapping AUs)       : "+skippedAUsNoOverlap);
407                logger.debug("  skipped (not overlapping chains)    : "+skippedChainsNoOverlap);
408                logger.debug("  skipped (sym redundant op pairs)    : "+skippedRedundant);
409                logger.debug("  skipped (sym redundant self op)     : "+skippedSelfEquivalent);
410                logger.debug("Found "+set.size()+" interfaces.");
411        }
412
413
414        /**
415         * Checks whether given interface is NCS-redundant, i.e., an identical interface between NCS copies of
416         * these molecules has already been seen, and returns this (reference) interface.
417         *
418         * @param interf
419         *          StructureInterface
420         * @return  already seen interface that is NCS-equivalent to interf,
421         *          null if such interface is not found.
422         */
423        private StructureInterface findNcsRef(StructureInterface interf) {
424                if (!this.hasNcsOps()) {
425                        return null;
426                }
427                String chainIName = interf.getMoleculeIds().getFirst();
428                String iOrigName = chainOrigNames.get(chainIName);
429
430                String chainJName = interf.getMoleculeIds().getSecond();
431                String jOrigName = chainOrigNames.get(chainJName);
432
433                Matrix4d mJCryst;
434                if(this.searchBeyondAU) {
435                        mJCryst = interf.getTransforms().getSecond().getMatTransform();
436                        mJCryst = crystallographicInfo.getCrystalCell().transfToOrthonormal(mJCryst);
437                } else {
438                        mJCryst = IDENTITY;
439                }
440
441                // Let X1,...Xn be the original coords, before NCS transforms (M1...Mk)
442                // current chain i: M_i * X_i
443                // current chain j: Cn * M_j * X_j
444
445                // transformation to bring chain j near X_i: M_i^(-1) * Cn * M_j
446                // transformation to bring chain i near X_j: (Cn * M_j)^(-1) * M_i = (M_i^(-1) * Cn * M_j)^(-1)
447
448                Matrix4d mChainIInv = new Matrix4d(chainNcsOps.get(chainIName));
449                mChainIInv.invert();
450
451                Matrix4d mJNcs = new Matrix4d(chainNcsOps.get(chainJName));
452
453                Matrix4d j2iNcsOrigin = new Matrix4d(mChainIInv);
454                j2iNcsOrigin.mul(mJCryst);
455                //overall transformation to bring current chainj from its NCS origin to i's
456                j2iNcsOrigin.mul(mJNcs);
457
458                //overall transformation to bring current chaini from its NCS origin to j's
459                Matrix4d i2jNcsOrigin = new Matrix4d(j2iNcsOrigin);
460                i2jNcsOrigin.invert();
461
462                String matchChainIdsIJ = iOrigName + jOrigName;
463                String matchChainIdsJI = jOrigName + iOrigName;
464
465                // same original chain names
466                Optional<Matrix4d> matchDirect =
467                                visitedNcsChainPairs.computeIfAbsent(matchChainIdsIJ, k-> new HashMap<>()).entrySet().stream().
468                                        map(r->r.getKey()).
469                                        filter(r->r.epsilonEquals(j2iNcsOrigin,0.01)).
470                                        findFirst();
471
472                Matrix4d matchMatrix = matchDirect.orElse(null);
473                String matchChainIds = matchChainIdsIJ;
474
475                if(matchMatrix == null) {
476                        // reversed original chain names with inverted transform
477                        Optional<Matrix4d> matchInverse =
478                                        visitedNcsChainPairs.computeIfAbsent(matchChainIdsJI, k-> new HashMap<>()).entrySet().stream().
479                                        map(r->r.getKey()).
480                                        filter(r->r.epsilonEquals(i2jNcsOrigin,0.01)).
481                                        findFirst();
482                        matchMatrix = matchInverse.orElse(null);
483                        matchChainIds = matchChainIdsJI;
484                }
485
486                StructureInterface matchInterface = null;
487
488                if (matchMatrix == null) {
489                        visitedNcsChainPairs.get(matchChainIdsIJ).put(j2iNcsOrigin,interf);
490                } else {
491                        matchInterface = visitedNcsChainPairs.get(matchChainIds).get(matchMatrix);
492                }
493
494                return matchInterface;
495        }
496
497        private StructureInterface calcContacts(Chain chaini, Chain chainj, double cutoff, CrystalTransform tt, StringBuilder builder) {
498                // note that we don't consider hydrogens when calculating contacts
499                AtomContactSet graph = StructureTools.getAtomsInContact(chaini, chainj, cutoff, INCLUDE_HETATOMS);
500
501                if (graph.size()>0) {
502                        if (builder != null) builder.append("x");
503
504                        CrystalTransform transf = new CrystalTransform(this.crystallographicInfo.getSpaceGroup());
505                        StructureInterface interf = new StructureInterface(
506                                        StructureTools.getAllAtomArray(chaini), StructureTools.getAllAtomArray(chainj),
507                                        chaini.getName(), chainj.getName(),
508                                        graph,
509                                        transf, tt);
510
511                        return interf;
512
513                } else {
514                        if (builder != null) builder.append("o");
515                        return null;
516                }
517        }
518
519        private void addVisitedTransform(CrystalTransform tt) {
520                visitedCrystalTransforms.add(tt);
521        }
522
523        /**
524         * Checks whether given transformId/translation is symmetry redundant
525         * Two transformations are symmetry redundant if their matrices (4d) multiplication gives the identity, i.e.
526         * if one is the inverse of the other.
527         * @param tt
528         * @return
529         */
530        private boolean isRedundantTransform(CrystalTransform tt) {
531
532                Iterator<CrystalTransform> it = visitedCrystalTransforms.iterator();
533                while (it.hasNext()) {
534                        CrystalTransform v = it.next();
535
536                        if (tt.isEquivalent(v)) {
537
538                                logger.debug("Skipping redundant transformation: "+tt+", equivalent to "+v);
539
540                                // there's only 1 possible equivalent partner for each visited matrix
541                                // (since the equivalent is its inverse matrix and the inverse matrix is unique)
542                                // thus once the partner has been seen, we don't need to check it ever again
543                                it.remove();
544
545                                return true;
546                        }
547                }
548
549                return false;
550        }
551
552        public void translate(Matrix4d m, Vector3d translation) {
553                m.m03 = m.m03+translation.x;
554                m.m13 = m.m13+translation.y;
555                m.m23 = m.m23+translation.z;
556        }
557
558        /**
559         * Apply the NCS operators in the given Structure adding new chains as needed.
560         * All chains are (re)assigned ids of the form: original_chain_id+ncs_operator_index+{@value #NCS_CHAINID_SUFFIX_CHAR}.
561         * @param structure
562         *          the structure to expand
563         * @param chainOrigNames
564         *          new chain names mapped to the original chain names
565         * @param chainNcsOps
566         *          new chain names mapped to the ncs operators that was used to generate them
567         * @since 5.0.0
568         */
569        public static void expandNcsOps(Structure structure, Map<String,String> chainOrigNames, Map<String,Matrix4d> chainNcsOps) {
570                PDBCrystallographicInfo xtalInfo = structure.getCrystallographicInfo();
571                if (xtalInfo ==null) return;
572
573                if (xtalInfo.getNcsOperators()==null || xtalInfo.getNcsOperators().length==0)
574                        return;
575
576                List<Chain> chainsToAdd = new ArrayList<>();
577
578                Matrix4d identity = new Matrix4d();
579                identity.setIdentity();
580
581                Matrix4d[] ncsOps = xtalInfo.getNcsOperators();
582
583                for (Chain c:structure.getChains()) {
584                        String cOrigId = c.getId();
585                        String cOrigName = c.getName();
586
587                        for (int iOperator = 0; iOperator < ncsOps.length; iOperator++) {
588                                Matrix4d m = ncsOps[iOperator];
589
590                                Chain clonedChain = (Chain)c.clone();
591                                String newChainId = cOrigId+(iOperator+1)+NCS_CHAINID_SUFFIX_CHAR;
592                                String newChainName = cOrigName+(iOperator+1)+NCS_CHAINID_SUFFIX_CHAR;
593                                clonedChain.setId(newChainId);
594                                clonedChain.setName(newChainName);
595
596                                setChainIdsInResidueNumbers(clonedChain, newChainName);
597                                Calc.transform(clonedChain, m);
598
599                                chainsToAdd.add(clonedChain);
600                                c.getEntityInfo().addChain(clonedChain);
601
602                                chainOrigNames.put(newChainName,cOrigName);
603                                chainNcsOps.put(newChainName,m);
604                        }
605
606                        chainNcsOps.put(cOrigName,identity);
607                        chainOrigNames.put(cOrigName,cOrigName);
608                }
609
610                chainsToAdd.forEach(structure::addChain);
611        }
612
613        /**
614         * Auxiliary method to reset chain ids of residue numbers in a chain.
615         * Used when cloning chains and resetting their ids: one needs to take care of
616         * resetting the ids within residue numbers too.
617         * @param c
618         * @param newChainName
619         */
620        private static void setChainIdsInResidueNumbers(Chain c, String newChainName) {
621                for (Group g:c.getAtomGroups()) {
622                        g.setResidueNumber(newChainName, g.getResidueNumber().getSeqNum(), g.getResidueNumber().getInsCode());
623                }
624                for (Group g:c.getSeqResGroups()) {
625                        if (g.getResidueNumber()==null) continue;
626                        g.setResidueNumber(newChainName, g.getResidueNumber().getSeqNum(), g.getResidueNumber().getInsCode());
627                }
628        }
629
630}