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.algorithms.misc;
39
40 import ffx.algorithms.AlgorithmListener;
41 import ffx.numerics.Potential;
42 import ffx.potential.MolecularAssembly;
43 import ffx.potential.bonded.AminoAcidUtils;
44 import ffx.potential.bonded.AminoAcidUtils.AminoAcid3;
45 import ffx.potential.bonded.Atom;
46 import ffx.potential.bonded.NucleicAcidUtils.NucleicAcid3;
47 import ffx.potential.bonded.Residue;
48 import ffx.potential.bonded.Rotamer;
49 import ffx.potential.bonded.RotamerLibrary;
50 import ffx.potential.parsers.PDBFilter;
51
52 import java.io.BufferedWriter;
53 import java.io.File;
54 import java.io.FileWriter;
55 import java.io.IOException;
56 import java.util.Arrays;
57 import java.util.logging.Level;
58 import java.util.logging.Logger;
59 import java.util.regex.Matcher;
60 import java.util.regex.Pattern;
61
62 import static java.lang.String.format;
63
64
65
66
67
68
69
70
71 public class GenerateRotamers {
72
73 private static final Logger logger = Logger.getLogger(GenerateRotamers.class.getName());
74 private static final Pattern atRangePatt = Pattern.compile("(\\d+)-(\\d+)");
75 private final Residue residue;
76 private final int nChi;
77 private final double[] currentChi;
78 private final File outFile;
79 private final MolecularAssembly molecularAssembly;
80 private final AlgorithmListener listener;
81 private final RotamerLibrary library;
82 private final Potential potential;
83 private double[] x;
84 private AminoAcid3 baselineAAres;
85 private Rotamer[] baselineRotamers;
86 private double incr = 10.0;
87 private boolean aroundLibrary = false;
88 private double width = 180.0;
89 private int startDepth = 0;
90 private int endDepth;
91 private boolean print = false;
92 private int nEval = 0;
93 private boolean writeVideo = false;
94 private File videoFile;
95 private PDBFilter videoFilter;
96
97
98
99
100
101
102
103
104
105
106
107 public GenerateRotamers(MolecularAssembly molecularAssembly, Potential potential, Residue residue,
108 File file, int nChi, AlgorithmListener listener) {
109 this(molecularAssembly, potential, residue, file, nChi, listener,
110 RotamerLibrary.getDefaultLibrary());
111 }
112
113
114
115
116
117
118
119
120
121
122
123
124 public GenerateRotamers(MolecularAssembly molecularAssembly, Potential potential, Residue residue,
125 File file, int nChi, AlgorithmListener listener, RotamerLibrary library) {
126 this.residue = residue;
127 this.nChi = nChi;
128 endDepth = nChi - 1;
129 this.potential = potential;
130 this.molecularAssembly = molecularAssembly;
131 this.listener = listener;
132 this.currentChi = new double[nChi];
133 this.library = library;
134 Arrays.fill(currentChi, 0.0);
135 File outputFile = file;
136 if (outputFile.exists()) {
137 String outName = outputFile.getName();
138 for (int i = 1; i < 1000; i++) {
139 outputFile = new File(format("%s_%d", outName, i));
140 if (!outputFile.exists()) {
141 break;
142 }
143 }
144 if (outputFile.exists()) {
145 logger.severe(format(" Could not version file %s", outName));
146 }
147 }
148 outFile = outputFile;
149 this.baselineAAres = residue.getAminoAcid3();
150 baselineRotamers = library.getRotamers(baselineAAres);
151 }
152
153
154
155
156
157
158 public void applyAndSaveTorsions(String[] torSets) {
159 for (String torSet : torSets) {
160 String[] torsions = torSet.split(",");
161 double[] values = new double[nChi * 2];
162 Arrays.fill(values, 0.0);
163 for (int i = 0; i < (Math.min(torsions.length, nChi)); i++) {
164 double chival = Double.parseDouble(torsions[i]);
165 currentChi[i] = chival;
166 values[2 * i] = chival;
167 }
168 Rotamer newRot = generateRotamer(values);
169 RotamerLibrary.applyRotamer(residue, newRot);
170 writeSnapshot();
171 }
172 }
173
174
175
176
177
178
179
180 public void setBaselineAARes(AminoAcid3 aa3) {
181 this.baselineAAres = aa3;
182 if (aa3 == AminoAcidUtils.AminoAcid3.UNK) {
183 boolean orig = library.getUsingOrigCoordsRotamer();
184 library.setUseOrigCoordsRotamer(false);
185 baselineRotamers = residue.setRotamers(library);
186 library.setUseOrigCoordsRotamer(orig);
187 } else {
188 baselineRotamers = library.getRotamers(aa3);
189 }
190 aroundLibrary = true;
191 }
192
193
194
195
196
197
198
199
200 public void setDepth(int start, int end) {
201 startDepth = start;
202 if (end < 1 || end > nChi) {
203 endDepth = nChi - 1;
204 } else {
205 endDepth = end;
206 }
207 }
208
209
210
211
212
213
214 public void setElectrostatics(String electrostatics) {
215 if (electrostatics != null) {
216 String[] tokens = electrostatics.split(",");
217 Atom[] atoms = molecularAssembly.getAtomArray();
218 for (String tok : tokens) {
219 Matcher m = atRangePatt.matcher(tok);
220 if (m.matches()) {
221 int begin = Integer.parseInt(m.group(1));
222 int end = Integer.parseInt(m.group(2));
223 logger.info(format(" Inactivating electrostatics for atoms %d-%d", begin, end));
224 for (int i = begin; i <= end; i++) {
225 Atom ai = atoms[i - 1];
226 ai.setElectrostatics(false);
227 ai.print(Level.FINE);
228 }
229 } else {
230 logger.info(format(" Discarding electrostatics input %s", tok));
231 }
232 }
233 }
234 }
235
236
237
238
239
240
241 public void setInactiveAtoms(String iatoms) {
242 if (iatoms != null) {
243 String[] tokens = iatoms.split(",");
244 Atom[] atoms = molecularAssembly.getAtomArray();
245 for (String tok : tokens) {
246 Matcher m = atRangePatt.matcher(tok);
247 if (m.matches()) {
248 int begin = Integer.parseInt(m.group(1));
249 int end = Integer.parseInt(m.group(2));
250 logger.info(format(" Inactivating atoms %d-%d", begin, end));
251 for (int i = begin; i <= end; i++) {
252 Atom ai = atoms[i - 1];
253 ai.setUse(false);
254 ai.print(Level.FINE);
255 }
256 } else {
257 logger.info(format(" Discarding inactive atoms input %s", tok));
258 }
259 }
260 }
261 }
262
263
264
265
266
267
268 public void setIncrement(double incr) {
269 this.incr = incr;
270 }
271
272
273
274
275
276
277 public void setPrint(boolean print) {
278 this.print = print;
279 }
280
281
282
283
284
285
286
287 public void setSearchWidth(double width) {
288 this.width = width;
289 }
290
291
292
293
294
295
296 public void setVideo(String videoFile) {
297 if (videoFile != null) {
298 File vidFile = new File(videoFile);
299 if (vidFile.exists()) {
300 for (int i = 0; i < 1000; i++) {
301 vidFile = new File(format("%s_%d", videoFile, i));
302 if (!vidFile.exists()) {
303 this.videoFile = vidFile;
304 writeVideo = true;
305 videoFilter = new PDBFilter(this.videoFile, molecularAssembly,
306 molecularAssembly.getForceField(), null);
307 videoFilter.setLogWrites(false);
308 break;
309 }
310 }
311 if (vidFile.exists()) {
312 logger.warning(format(" Could not version video file %s", videoFile));
313 }
314 } else {
315 this.videoFile = vidFile;
316 writeVideo = true;
317 videoFilter = new PDBFilter(this.videoFile, molecularAssembly,
318 molecularAssembly.getForceField(), null);
319 videoFilter.setLogWrites(false);
320 }
321 } else {
322 writeVideo = false;
323 }
324 }
325
326
327 public void tryRotamers() {
328 try (BufferedWriter bw = new BufferedWriter(new FileWriter(outFile))) {
329 if (aroundLibrary) {
330 for (Rotamer rotamer : baselineRotamers) {
331 applyLibRotamer(rotamer);
332 turnChi(startDepth, bw);
333 }
334 } else {
335 RotamerLibrary.measureRotamer(residue, currentChi, false);
336 turnChi(startDepth, bw);
337 }
338 } catch (IOException ex) {
339 logger.warning(format(" IO exception in rotamer generation: %s", ex));
340 }
341 }
342
343
344
345
346
347
348
349 private void applyLibRotamer(Rotamer rotamer) {
350 double[] rotValues = new double[nChi * 2];
351 Arrays.fill(rotValues, 0.0);
352 Arrays.fill(currentChi, 0.0);
353 int fillTo = Math.min(startDepth, (rotamer.length - 1));
354 for (int i = 0; i < fillTo; i++) {
355 int ii = 2 * i;
356 rotValues[ii] = rotamer.angles[i];
357 currentChi[i] = rotamer.angles[i];
358 rotValues[ii + 1] = rotamer.sigmas[i];
359 }
360 Rotamer newRot = generateRotamer(rotValues);
361 RotamerLibrary.applyRotamer(residue, newRot);
362 }
363
364
365
366
367
368
369
370 private Rotamer generateRotamer(double[] values) {
371 switch (residue.getResidueType()) {
372 case AA -> {
373 AminoAcid3 aa3 = residue.getAminoAcid3();
374 return new Rotamer(aa3, values);
375 }
376 case NA -> {
377 NucleicAcid3 na3 = residue.getNucleicAcid3();
378 return new Rotamer(na3, values);
379 }
380 default -> {
381 return new Rotamer(values);
382 }
383 }
384 }
385
386
387
388
389
390
391
392
393 private void turnChi(int depth, BufferedWriter bw) throws IOException {
394 double chi = currentChi[depth];
395 for (double i = chi - width; i <= chi + width; i += incr) {
396 currentChi[depth] = i;
397 if (depth == endDepth) {
398 evaluateChi(bw);
399 if (listener != null) {
400 listener.algorithmUpdate(molecularAssembly);
401 }
402 if (writeVideo) {
403 writeSnapshot();
404 }
405 } else {
406 turnChi(depth + 1, bw);
407 }
408 }
409 currentChi[depth] = chi;
410 }
411
412 private void writeSnapshot() {
413 try (BufferedWriter bw = new BufferedWriter(new FileWriter(videoFile, true))) {
414 bw.write(format("MODEL %d", nEval));
415 bw.newLine();
416 bw.write(format("REMARK 301 TORSIONS %s", formatChi()));
417 bw.newLine();
418 bw.flush();
419 videoFilter.writeFile(videoFile, true);
420 bw.write("ENDMDL");
421 bw.newLine();
422 bw.flush();
423 } catch (IOException ex) {
424 logger.warning(
425 format(" Exception writing to video file %s: " + "%s", videoFile.getName(), ex));
426 }
427 }
428
429
430
431
432
433
434
435 private void evaluateChi(BufferedWriter bw) throws IOException {
436 try {
437 applyChi();
438 double e = currentEnergy();
439 String result = format("%s,%12f", formatChi(), e);
440 ++nEval;
441 if (print) {
442 logger.info(format(" Evaluation %10d %s, energy %10.5f kcal/mol", nEval, formatChi(), e));
443 }
444 bw.write(result);
445 bw.newLine();
446 if (nEval % 1000 == 0) {
447 logger.info(format(" %12.7e states evaluated", (double) nEval));
448 }
449 } catch (ArithmeticException ex) {
450 logger.info(format(" Force field exception at chi %s", formatChi()));
451 }
452 }
453
454
455 private void applyChi() {
456 double[] rotValues = new double[nChi * 2];
457 for (int i = 0; i < nChi; i++) {
458 int ii = 2 * i;
459 rotValues[ii] = currentChi[i];
460 rotValues[ii + 1] = 0.0;
461 }
462 RotamerLibrary.applyRotamer(residue, generateRotamer(rotValues));
463 }
464
465
466
467
468
469
470 private String formatChi() {
471 StringBuilder sb = new StringBuilder(format("%8f", currentChi[0]));
472 for (int i = 1; i < nChi; i++) {
473
474 sb.append(format(",%8f", currentChi[i]));
475 }
476 return sb.toString();
477 }
478
479
480
481
482
483
484 private double currentEnergy() throws ArithmeticException {
485 if (x == null) {
486 int nVar = potential.getNumberOfVariables();
487 x = new double[nVar * 3];
488 }
489 potential.getCoordinates(x);
490 return potential.energy(x);
491 }
492 }