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