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.xray.commands;
39
40 import ffx.algorithms.cli.AlgorithmsCommand;
41 import ffx.algorithms.cli.ManyBodyOptions;
42 import ffx.algorithms.optimize.RotamerOptimization;
43 import ffx.algorithms.optimize.TitrationManyBody;
44 import ffx.numerics.Potential;
45 import ffx.potential.ForceFieldEnergy;
46 import ffx.potential.MolecularAssembly;
47 import ffx.potential.bonded.Atom;
48 import ffx.potential.bonded.LambdaInterface;
49 import ffx.potential.bonded.Residue;
50 import ffx.potential.bonded.Rotamer;
51 import ffx.potential.bonded.RotamerLibrary;
52 import ffx.potential.cli.AlchemicalOptions;
53 import ffx.potential.parameters.ForceField;
54 import ffx.potential.parsers.PDBFilter;
55 import ffx.utilities.FFXBinding;
56 import ffx.xray.DiffractionData;
57 import ffx.xray.RefinementEnergy;
58 import ffx.xray.refine.RefinementMode;
59 import ffx.xray.cli.XrayOptions;
60 import org.apache.commons.configuration2.CompositeConfiguration;
61 import picocli.CommandLine.Command;
62 import picocli.CommandLine.Mixin;
63 import picocli.CommandLine.Option;
64 import picocli.CommandLine.Parameters;
65
66 import java.io.File;
67 import java.io.FileWriter;
68 import java.io.IOException;
69 import java.util.ArrayList;
70 import java.util.Arrays;
71 import java.util.HashSet;
72 import java.util.List;
73 import java.util.Set;
74 import java.util.stream.Collectors;
75
76 import static ffx.potential.bonded.NamingUtils.renameAtomsToPDBStandard;
77 import static java.lang.String.format;
78
79
80
81
82
83
84
85
86 @Command(description = " Run GenZ function for free energy change.", name = "xray.GenZ")
87 public class GenZ extends AlgorithmsCommand {
88
89 @Mixin
90 private ManyBodyOptions manyBodyOptions;
91
92 @Mixin
93 private AlchemicalOptions alchemicalOptions;
94
95 @Mixin
96 private XrayOptions xrayOptions;
97
98 @Option(names = {"--rEE", "--ro-ensembleEnergy"}, paramLabel = "0.0",
99 description = "Keep permutations within ensemble Energy kcal/mol from the GMEC.")
100 private String ensembleEnergy = "0.0";
101
102 @Option(names = {"--pF", "--printFiles"}, paramLabel = "false",
103 description = "Write to an energy restart file and ensemble file.")
104 private boolean printFiles = false;
105
106 @Option(names = {"--pKa"}, paramLabel = "false",
107 description = "Calculating protonation populations for pKa shift.")
108 private boolean pKa = false;
109
110
111
112
113 @Parameters(arity = "1..*", paramLabel = "files", description = "PDB and Real Space input files.")
114 private List<String> filenames;
115
116 private RefinementEnergy refinementEnergy;
117 private ForceFieldEnergy potentialEnergy;
118 private MolecularAssembly[] molecularAssemblies;
119 private ForceField forceField;
120 private TitrationManyBody titrationManyBody;
121 private List<Residue> selectedResidues;
122 private MolecularAssembly[] conformerAssemblies = new MolecularAssembly[3];
123
124
125
126
127 public GenZ() {
128 super();
129 }
130
131
132
133
134
135
136 public GenZ(String[] args) {
137 super(args);
138 }
139
140
141
142
143
144
145 public GenZ(FFXBinding binding) {
146 super(binding);
147 }
148
149
150
151
152 @Override
153 public GenZ run() {
154 if (!init()) {
155 return this;
156 }
157
158 xrayOptions.init();
159 double titrationPH = manyBodyOptions.getTitrationPH();
160 double inclusionCutoff = manyBodyOptions.getInclusionCutoff();
161 int mutatingResidue = manyBodyOptions.getInterestedResidue();
162 boolean onlyTitration = manyBodyOptions.getOnlyTitration();
163 double pHRestraint = manyBodyOptions.getPHRestraint();
164 if (manyBodyOptions.getTitration()) {
165 System.setProperty("manybody-titration", "true");
166 }
167
168
169 String nea = System.getProperty("native-environment-approximation", "true");
170 System.setProperty("native-environment-approximation", nea);
171
172 boolean lambdaTerm = alchemicalOptions.hasSoftcore();
173 if (lambdaTerm) {
174
175 System.setProperty("lambdaterm", "true");
176
177 System.setProperty("elec-lambdaterm", "false");
178
179 System.setProperty("intramolecular-softcore", "true");
180 }
181 System.setProperty("ro-ensembleEnergy", ensembleEnergy);
182
183 String modelFilename;
184 if (filenames != null && !filenames.isEmpty()) {
185 molecularAssemblies = algorithmFunctions.openAll(filenames.get(0));
186 activeAssembly = molecularAssemblies[0];
187 modelFilename = filenames.get(0);
188 } else if (activeAssembly == null) {
189 logger.info(helpString());
190 return this;
191 } else {
192 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
193 modelFilename = activeAssembly.getFile().getAbsolutePath();
194 }
195
196 CompositeConfiguration properties = activeAssembly.getProperties();
197 activeAssembly.getPotentialEnergy().setPrintOnFailure(false, false);
198 potentialEnergy = activeAssembly.getPotentialEnergy();
199 forceField = activeAssembly.getForceField();
200
201 if (forceField == null) {
202 logger.info("This force field is null");
203 }
204
205 String filename = filenames.get(0);
206
207 List<Residue> titrateResidues = new ArrayList<>();
208
209
210 Set<Atom> excludeAtoms = new HashSet<>();
211 boolean isTitrating = false;
212
213
214 if (xrayOptions.refinementMode != RefinementMode.COORDINATES) {
215 logger.info(" Refinement mode set to COORDINATES.");
216 xrayOptions.refinementMode = RefinementMode.COORDINATES;
217 }
218
219
220 List<Residue> residues = manyBodyOptions.collectResidues(activeAssembly);
221 if (residues == null || residues.isEmpty()) {
222 logger.info(" There are no residues in the active system to optimize.");
223 return this;
224 }
225
226 String[] titratableResidues = {"HIS", "HIE", "HID", "GLU", "GLH", "ASP", "ASH", "LYS", "LYD", "CYS", "CYD"};
227 List<String> titratableResiudesList = Arrays.asList(titratableResidues);
228 double boltzmannWeights;
229 double offsets;
230 double[][] populationArray = new double[1][55];
231 double[][] titrateBoltzmann;
232 double[] protonationBoltzmannSums;
233 double totalBoltzmann = 0;
234 List<Residue> residueList = activeAssembly.getResidueList();
235
236 String listResidues = "";
237
238
239 if (inclusionCutoff != -1 || onlyTitration) {
240 listResidues = manyBodyOptions.selectInclusionResidues(residueList, mutatingResidue, onlyTitration, inclusionCutoff);
241 }
242
243 List<Integer> residueNumber = new ArrayList<>();
244 for (Residue residue : residueList) {
245 residueNumber.add(residue.getResidueNumber());
246 }
247
248
249 if (properties.getBoolean("standardizeAtomNames", false)) {
250 renameAtomsToPDBStandard(activeAssembly);
251 }
252
253
254 if (manyBodyOptions.getTitration()) {
255 logger.info("\n Adding titration hydrogen to : " + filenames.get(0) + "\n");
256
257
258 List<Integer> resNumberList = new ArrayList<>();
259 for (Residue residue : residues) {
260 resNumberList.add(residue.getResidueNumber());
261 }
262
263
264 titrationManyBody = new TitrationManyBody(filename, activeAssembly.getForceField(),
265 resNumberList, titrationPH, manyBodyOptions);
266 activeAssembly = titrationManyBody.getProtonatedAssembly();
267 potentialEnergy = activeAssembly.getPotentialEnergy();
268 }
269
270
271 xrayOptions.setProperties(parseResult, properties);
272 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
273
274 DiffractionData diffractionData = xrayOptions.getDiffractionData(filenames, molecularAssemblies, properties);
275 refinementEnergy = xrayOptions.toXrayEnergy(diffractionData);
276 refinementEnergy.setScaling(null);
277
278
279 if (!pKa || onlyTitration) {
280 manyBodyOptions.setListResidues(listResidues);
281 }
282
283 if (lambdaTerm) {
284 alchemicalOptions.setFirstSystemAlchemistry(activeAssembly);
285 LambdaInterface lambdaInterface = (LambdaInterface) potentialEnergy;
286 double lambda = alchemicalOptions.getInitialLambda();
287 logger.info(format(" Setting ManyBody softcore lambda to: %5.3f", lambda));
288 lambdaInterface.setLambda(lambda);
289 }
290
291 RotamerOptimization rotamerOptimization = new RotamerOptimization(activeAssembly, refinementEnergy, algorithmListener);
292 manyBodyOptions.initRotamerOptimization(rotamerOptimization, activeAssembly);
293 rotamerOptimization.setPrintFiles(printFiles);
294 rotamerOptimization.setWriteEnergyRestart(printFiles);
295 rotamerOptimization.setPHRestraint(pHRestraint);
296 rotamerOptimization.setpH(titrationPH);
297 double[] x = new double[refinementEnergy.getNumberOfVariables()];
298 x = refinementEnergy.getCoordinates(x);
299 double e = refinementEnergy.energy(x, true);
300 logger.info(format("\n Initial target energy: %16.8f ", e));
301
302 selectedResidues = rotamerOptimization.getResidues();
303 rotamerOptimization.initFraction(selectedResidues);
304
305 RotamerLibrary.measureRotamers(selectedResidues, false);
306
307 rotamerOptimization.optimize(manyBodyOptions.getAlgorithm(selectedResidues.size()));
308
309 int[] currentRotamers = new int[selectedResidues.size()];
310
311
312 try {
313 rotamerOptimization.getFractions(selectedResidues.toArray(new Residue[0]), 0, currentRotamers);
314 if (pKa) {
315 rotamerOptimization.getProtonationPopulations(selectedResidues.toArray(new Residue[0]));
316 }
317 } catch (Exception ex) {
318 logger.severe(" Error calculating rotamer fractions: " + ex.getMessage());
319 return this;
320 }
321
322
323 boltzmannWeights = rotamerOptimization.getTotalBoltzmann();
324 offsets = rotamerOptimization.getRefEnergy();
325
326
327 populationArray = rotamerOptimization.getFraction();
328
329
330
331
332
333
334 try (FileWriter fileWriter = new FileWriter("populations.txt")) {
335 int residueIndex = 0;
336 for (Residue residue : selectedResidues) {
337 fileWriter.write("\n");
338 Rotamer[] rotamers = residue.getRotamers();
339 for (int rotIndex = 0; rotIndex < rotamers.length; rotIndex++) {
340 String rotPop = format("%.6f", populationArray[residueIndex][rotIndex]);
341 fileWriter.write(residue.getName() + residue.getResidueNumber() + "\t" +
342 rotamers[rotIndex].toString() + "\t" + rotPop + "\n");
343 }
344 residueIndex += 1;
345 }
346 } catch (IOException ioException) {
347 logger.warning("Error writing populations file: " + ioException.getMessage());
348 }
349
350 System.out.println("\n Successfully wrote to the populations file.");
351
352 int[][] conformers;
353 try {
354 conformers = rotamerOptimization.getConformers();
355 } catch (Exception ex) {
356 logger.severe(" Error getting conformers: " + ex.getMessage());
357 return this;
358 }
359
360 char[] altLocs = new char[]{'C', 'B', 'A'};
361 List<String> residueChainNum = new ArrayList<>();
362 for (Residue residue : activeAssembly.getResidueList()) {
363 char chain = residue.getChainID();
364 int resNum = residue.getResidueNumber();
365 residueChainNum.add(String.valueOf(chain) + resNum);
366 }
367
368 File structureFile = new File(filename);
369
370 int[] optimalRotamers = new int[selectedResidues.size()];
371 int assemblyIndex = 0;
372 String[] rotNames = new String[selectedResidues.size()];
373 for (int confIndex = 2; confIndex > -1; confIndex--) {
374 List<Residue> conformerResidueList = new ArrayList<>();
375 MolecularAssembly conformerAssembly = algorithmFunctions.open(filename);
376 if (manyBodyOptions.getTitration()) {
377 logger.info("\n Adding titration hydrogen to : " + filenames.get(0) + "\n");
378
379
380 List<Integer> resNumberList = new ArrayList<>();
381 for (Residue residue : residues) {
382 resNumberList.add(residue.getResidueNumber());
383 }
384
385
386 titrationManyBody = new TitrationManyBody(filename, activeAssembly.getForceField(),
387 resNumberList, titrationPH, manyBodyOptions);
388 conformerAssembly = titrationManyBody.getProtonatedAssembly();
389 potentialEnergy = conformerAssembly.getPotentialEnergy();
390 }
391 for (int resIndex = 0; resIndex < selectedResidues.size(); resIndex++) {
392 Residue residueSelect = selectedResidues.get(resIndex);
393 String resChainNum = String.valueOf(residueSelect.getChainID()) + residueSelect.getResidueNumber();
394 int index = residueChainNum.indexOf(resChainNum);
395 Residue residue = conformerAssembly.getResidueList().get(index);
396 conformerResidueList.add(residue);
397 residue.setRotamers(manyBodyOptions.getRotamerLibrary(true));
398 Rotamer[] rotamers = residue.getRotamers();
399 int rotIndex = conformers[resIndex][confIndex];
400 if (populationArray[resIndex][rotIndex] != 0) {
401 optimalRotamers[resIndex] = rotIndex;
402 RotamerLibrary.applyRotamer(residue, rotamers[rotIndex]);
403 double occupancy = populationArray[resIndex][rotIndex];
404 boolean diffStates = false;
405 for (int i = 2; i > -1; i--) {
406 int rotamerInd = conformers[resIndex][i];
407 String rotName = rotamers[rotamerInd].getName();
408 double occupancyTest = populationArray[resIndex][rotamerInd];
409 if (i == 2 && rotNames[resIndex] != null) {
410 rotNames[resIndex] = rotName;
411 } else if (i < 2 && occupancyTest != 0 && !rotNames[resIndex].contains(rotName)) {
412 diffStates = true;
413 String newString = rotNames[resIndex] + rotName;
414 logger.info(newString);
415 rotNames[resIndex] = newString;
416 } else if (i == 2) {
417 rotNames[resIndex] = rotName;
418 }
419 }
420 for (Atom atom : residue.getAtomList()) {
421 if (!residue.getBackboneAtoms().contains(atom) || diffStates) {
422 if (occupancy != 1) {
423 atom.setAltLoc(altLocs[confIndex]);
424 atom.setOccupancy(occupancy);
425 } else {
426 atom.setOccupancy(occupancy);
427 atom.setAltLoc(' ');
428 }
429 } else {
430 atom.setOccupancy(1.0);
431 atom.setAltLoc(' ');
432 }
433 }
434 }
435 }
436
437 conformerAssemblies[assemblyIndex] = conformerAssembly;
438 if (manyBodyOptions.getTitration()) {
439 titrationManyBody.excludeExcessAtoms(excludeAtoms, optimalRotamers, conformerAssemblies[assemblyIndex], conformerResidueList);
440 excludeAtoms.addAll(titrationManyBody.getExcludeAtoms());
441 }
442 assemblyIndex++;
443
444 if (pKa) {
445 rotamerOptimization.getProtonationPopulations(selectedResidues.toArray(new Residue[0]));
446 }
447 }
448
449
450 try (FileWriter popFileWriter = new FileWriter("populations.txt")) {
451 int residueIndex = 0;
452 for (Residue residue : selectedResidues) {
453 popFileWriter.write("\n");
454 protonationBoltzmannSums = new double[selectedResidues.size()];
455
456 Rotamer[] rotamers = residue.getRotamers();
457 for (int rotIndex = 0; rotIndex < rotamers.length; rotIndex++) {
458 String rotPop = format("%.6f", populationArray[residueIndex][rotIndex]);
459 popFileWriter.write(residue.getName() + residue.getResidueNumber() + "\t" +
460 rotamers[rotIndex].toString() + "\t" + rotPop + "\n");
461 }
462 residueIndex += 1;
463 }
464 } catch (IOException ioException) {
465 logger.warning("Error writing populations file: " + ioException.getMessage());
466 }
467
468 System.out.println("\n Successfully wrote to the populations file.");
469
470 PDBFilter pdbFilter = new PDBFilter(structureFile, Arrays.asList(conformerAssemblies), forceField, properties);
471 pdbFilter.writeFile(structureFile, false, excludeAtoms, true, true);
472
473 System.setProperty("standardizeAtomNames", "false");
474
475 return this;
476 }
477
478 @Override
479 public List<Potential> getPotentials() {
480 if (conformerAssemblies == null) {
481 return new ArrayList<>();
482 } else {
483 return Arrays.stream(conformerAssemblies)
484 .filter(a -> a != null)
485 .map(MolecularAssembly::getPotentialEnergy)
486 .filter(e -> e != null)
487 .collect(Collectors.toList());
488 }
489 }
490 }