1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.commands;
39
40 import ffx.potential.bonded.Atom;
41 import ffx.potential.bonded.Residue;
42 import ffx.potential.bonded.Residue.ResidueType;
43 import ffx.potential.cli.PotentialCommand;
44 import ffx.potential.utils.PotentialsUtils;
45 import ffx.utilities.FFXBinding;
46 import ffx.utilities.FilePathInfo;
47 import picocli.CommandLine.Command;
48 import picocli.CommandLine.Parameters;
49
50 import java.util.ArrayList;
51 import java.util.List;
52
53 import static ffx.numerics.math.DoubleMath.dist;
54 import static java.lang.String.format;
55 import static org.apache.commons.math3.util.FastMath.abs;
56
57
58
59
60
61
62
63 @Command(name = "ChainBreaks", description = " Fix chain breaks in a pdb file.")
64 public class ChainBreaks extends PotentialCommand {
65
66
67
68
69 @Parameters(arity = "1", paramLabel = "file",
70 description = "The atomic coordinate file in PDB or XYZ format.")
71 private String filename = null;
72
73 private List<double[]> newCoordinates;
74 private final PotentialsUtils potentialsUtils = new PotentialsUtils();
75
76 public ChainBreaks() {
77 super();
78 }
79
80 public ChainBreaks(FFXBinding binding) {
81 super(binding);
82 }
83
84 public ChainBreaks(String[] args) {
85 super(args);
86 }
87
88 @Override
89 public ChainBreaks run() {
90 if (!init()) {
91 return null;
92 }
93
94
95 activeAssembly = getActiveAssembly(filename);
96 if (activeAssembly == null) {
97 logger.info(helpString());
98 return null;
99 }
100
101
102 filename = activeAssembly.getFile().getAbsolutePath();
103
104
105 FilePathInfo fileInfo = parseFilePath(filename);
106
107 if (!fileInfo.extension().toLowerCase().contains("pdb")) {
108 logger.info(format(" The file extension does not include 'pdb': %s", filename));
109 return null;
110 }
111
112 logger.info(format("\n Fixing Chain Breaks for %s", filename));
113 List<Residue> residues = activeAssembly.getResidueList();
114 List<String> chainBreaks = findChainBreaks(residues, 3.0);
115 saveByOriginalExtension(activeAssembly, filename);
116 return this;
117 }
118
119 private List<String> findChainBreaks(List<Residue> residues, double cutoff) {
120 List<List<Residue>> subChains = new ArrayList<>();
121 List<String> chainBreaks = new ArrayList<>();
122
123
124
125 ResidueType rType = residues.get(0).getResidueType();
126 String startAtName;
127 String endAtName;
128 switch (rType) {
129 case AA:
130 startAtName = "N";
131 endAtName = "C";
132 break;
133 case NA:
134 boolean namedStar = residues.stream()
135 .flatMap((Residue r) -> r.getAtomList().stream())
136 .anyMatch((Atom a) -> a.getName().equals("O5*"));
137 if (namedStar) {
138 startAtName = "O5*";
139 endAtName = "O3*";
140 } else {
141 startAtName = "O5'";
142 endAtName = "O3'";
143 }
144 break;
145 case UNK:
146 default:
147 logger.fine(" Not attempting to find chain breaks for chain with residue " + residues.get(0));
148 return null;
149 }
150
151 List<Residue> subChain = null;
152 Residue previousResidue = null;
153 Atom priorEndAtom = null;
154
155 for (Residue residue : residues) {
156 List<Atom> resAtoms = residue.getAtomList();
157 if (priorEndAtom == null) {
158
159 subChain = new ArrayList<>();
160 subChain.add(residue);
161 subChains.add(subChain);
162 } else {
163
164 Atom startAtom = null;
165 for (Atom a : resAtoms) {
166 if (a.getName().equalsIgnoreCase(startAtName)) {
167 startAtom = a;
168 break;
169 }
170 }
171 if (startAtom == null) {
172 subChain.add(residue);
173 continue;
174 }
175
176 double r = dist(priorEndAtom.getXYZ(null), startAtom.getXYZ(null));
177
178 if (r > cutoff) {
179
180 subChain = new ArrayList<>();
181 subChain.add(residue);
182 subChains.add(subChain);
183 chainBreaks.add("C " + previousResidue + " N " + residue);
184 fixChainBreaks(priorEndAtom.getXYZ(null), startAtom.getXYZ(null));
185 priorEndAtom.setXYZ(newCoordinates.get(0));
186 startAtom.setXYZ(newCoordinates.get(1));
187 } else {
188
189 subChain.add(residue);
190 }
191 }
192
193
194 for (Atom a : resAtoms) {
195 if (a.getName().equalsIgnoreCase(endAtName)) {
196 priorEndAtom = a;
197 break;
198 }
199 }
200 previousResidue = residue;
201 }
202
203 return chainBreaks;
204 }
205
206 private void fixChainBreaks(double[] cCoordinates, double[] nCoordinates) {
207 logger.info(" Generating new coordinates.");
208 newCoordinates = new ArrayList<>();
209 double distance = dist(cCoordinates, nCoordinates);
210 while (distance > 3.0) {
211 double dx = abs((cCoordinates[0] - nCoordinates[0]) / 4.0);
212 double dy = abs((cCoordinates[1] - nCoordinates[1]) / 4.0);
213 double dz = abs((cCoordinates[2] - nCoordinates[2]) / 4.0);
214
215 if (cCoordinates[0] > nCoordinates[0]) {
216 cCoordinates[0] = cCoordinates[0] - dx;
217 nCoordinates[0] = nCoordinates[0] + dx;
218 } else if (cCoordinates[0] < nCoordinates[0]) {
219 cCoordinates[0] = cCoordinates[0] + dx;
220 nCoordinates[0] = nCoordinates[0] - dx;
221 }
222
223 if (cCoordinates[1] > nCoordinates[1]) {
224 cCoordinates[1] = cCoordinates[1] - dy;
225 nCoordinates[1] = nCoordinates[1] + dy;
226 } else if (cCoordinates[1] < nCoordinates[1]) {
227 cCoordinates[1] = cCoordinates[1] + dy;
228 nCoordinates[1] = nCoordinates[1] - dy;
229 }
230
231 if (cCoordinates[2] > nCoordinates[2]) {
232 cCoordinates[2] = cCoordinates[2] - dz;
233 nCoordinates[2] = nCoordinates[2] + dz;
234 } else if (cCoordinates[2] < nCoordinates[0]) {
235 cCoordinates[2] = cCoordinates[2] + dz;
236 nCoordinates[2] = nCoordinates[2] - dz;
237 }
238
239 distance = dist(cCoordinates, nCoordinates);
240 }
241 newCoordinates.add(cCoordinates);
242 newCoordinates.add(nCoordinates);
243 }
244 }