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.RefinementMinimize;
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 != RefinementMinimize.RefinementMode.COORDINATES) {
215 logger.info(" Refinement mode set to COORDINATES.");
216 xrayOptions.refinementMode = RefinementMinimize.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 MolecularAssembly protonatedAssembly = titrationManyBody.getProtonatedAssembly();
267 setActiveAssembly(protonatedAssembly);
268 potentialEnergy = protonatedAssembly.getPotentialEnergy();
269 }
270
271
272 xrayOptions.setProperties(parseResult, properties);
273 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
274
275 DiffractionData diffractionData = xrayOptions.getDiffractionData(filenames, molecularAssemblies, properties);
276 refinementEnergy = xrayOptions.toXrayEnergy(diffractionData);
277 refinementEnergy.setScaling(null);
278
279
280 if (!pKa || onlyTitration) {
281 manyBodyOptions.setListResidues(listResidues);
282 }
283
284 if (lambdaTerm) {
285 alchemicalOptions.setFirstSystemAlchemistry(activeAssembly);
286 LambdaInterface lambdaInterface = (LambdaInterface) potentialEnergy;
287 double lambda = alchemicalOptions.getInitialLambda();
288 logger.info(format(" Setting ManyBody softcore lambda to: %5.3f", lambda));
289 lambdaInterface.setLambda(lambda);
290 }
291
292 RotamerOptimization rotamerOptimization = new RotamerOptimization(activeAssembly, refinementEnergy, algorithmListener);
293 manyBodyOptions.initRotamerOptimization(rotamerOptimization, activeAssembly);
294 rotamerOptimization.setPrintFiles(printFiles);
295 rotamerOptimization.setWriteEnergyRestart(printFiles);
296 rotamerOptimization.setPHRestraint(pHRestraint);
297 rotamerOptimization.setpH(titrationPH);
298 double[] x = new double[refinementEnergy.getNumberOfVariables()];
299 x = refinementEnergy.getCoordinates(x);
300 double e = refinementEnergy.energy(x, true);
301 logger.info(format("\n Initial target energy: %16.8f ", e));
302
303 selectedResidues = rotamerOptimization.getResidues();
304 rotamerOptimization.initFraction(selectedResidues);
305
306 RotamerLibrary.measureRotamers(selectedResidues, false);
307
308 rotamerOptimization.optimize(manyBodyOptions.getAlgorithm(selectedResidues.size()));
309
310 int[] currentRotamers = new int[selectedResidues.size()];
311
312
313 try {
314 rotamerOptimization.getFractions(selectedResidues.toArray(new Residue[0]), 0, currentRotamers);
315 if (pKa) {
316 rotamerOptimization.getProtonationPopulations(selectedResidues.toArray(new Residue[0]));
317 }
318 } catch (Exception ex) {
319 logger.severe(" Error calculating rotamer fractions: " + ex.getMessage());
320 return this;
321 }
322
323
324 boltzmannWeights = rotamerOptimization.getTotalBoltzmann();
325 offsets = rotamerOptimization.getRefEnergy();
326
327
328 populationArray = rotamerOptimization.getFraction();
329
330
331
332
333
334
335 try (FileWriter fileWriter = new FileWriter("populations.txt")) {
336 int residueIndex = 0;
337 for (Residue residue : selectedResidues) {
338 fileWriter.write("\n");
339 Rotamer[] rotamers = residue.getRotamers();
340 for (int rotIndex = 0; rotIndex < rotamers.length; rotIndex++) {
341 String rotPop = format("%.6f", populationArray[residueIndex][rotIndex]);
342 fileWriter.write(residue.getName() + residue.getResidueNumber() + "\t" +
343 rotamers[rotIndex].toString() + "\t" + rotPop + "\n");
344 }
345 residueIndex += 1;
346 }
347 } catch (IOException ioException) {
348 logger.warning("Error writing populations file: " + ioException.getMessage());
349 }
350
351 System.out.println("\n Successfully wrote to the populations file.");
352
353 int[][] conformers;
354 try {
355 conformers = rotamerOptimization.getConformers();
356 } catch (Exception ex) {
357 logger.severe(" Error getting conformers: " + ex.getMessage());
358 return this;
359 }
360
361 char[] altLocs = new char[]{'C', 'B', 'A'};
362 List<String> residueChainNum = new ArrayList<>();
363 for (Residue residue : activeAssembly.getResidueList()) {
364 char chain = residue.getChainID();
365 int resNum = residue.getResidueNumber();
366 residueChainNum.add(String.valueOf(chain) + resNum);
367 }
368
369 File structureFile = new File(filename);
370
371 int[] optimalRotamers = new int[selectedResidues.size()];
372 int assemblyIndex = 0;
373 String[] rotNames = new String[selectedResidues.size()];
374 for (int confIndex = 2; confIndex > -1; confIndex--) {
375 List<Residue> conformerResidueList = new ArrayList<>();
376 MolecularAssembly conformerAssembly = algorithmFunctions.open(filename);
377 if (manyBodyOptions.getTitration()) {
378 logger.info("\n Adding titration hydrogen to : " + filenames.get(0) + "\n");
379
380
381 List<Integer> resNumberList = new ArrayList<>();
382 for (Residue residue : residues) {
383 resNumberList.add(residue.getResidueNumber());
384 }
385
386
387 titrationManyBody = new TitrationManyBody(filename, activeAssembly.getForceField(),
388 resNumberList, titrationPH, manyBodyOptions);
389 conformerAssembly = titrationManyBody.getProtonatedAssembly();
390 potentialEnergy = conformerAssembly.getPotentialEnergy();
391 }
392 for (int resIndex = 0; resIndex < selectedResidues.size(); resIndex++) {
393 Residue residueSelect = selectedResidues.get(resIndex);
394 String resChainNum = String.valueOf(residueSelect.getChainID()) + residueSelect.getResidueNumber();
395 int index = residueChainNum.indexOf(resChainNum);
396 Residue residue = conformerAssembly.getResidueList().get(index);
397 conformerResidueList.add(residue);
398 residue.setRotamers(manyBodyOptions.getRotamerLibrary(true));
399 Rotamer[] rotamers = residue.getRotamers();
400 int rotIndex = conformers[resIndex][confIndex];
401 if (populationArray[resIndex][rotIndex] != 0) {
402 optimalRotamers[resIndex] = rotIndex;
403 RotamerLibrary.applyRotamer(residue, rotamers[rotIndex]);
404 double occupancy = populationArray[resIndex][rotIndex];
405 boolean diffStates = false;
406 for (int i = 2; i > -1; i--) {
407 int rotamerInd = conformers[resIndex][i];
408 String rotName = rotamers[rotamerInd].getName();
409 double occupancyTest = populationArray[resIndex][rotamerInd];
410 if (i == 2 && rotNames[resIndex] != null) {
411 rotNames[resIndex] = rotName;
412 } else if (i < 2 && occupancyTest != 0 && !rotNames[resIndex].contains(rotName)) {
413 diffStates = true;
414 String newString = rotNames[resIndex] + rotName;
415 logger.info(newString);
416 rotNames[resIndex] = newString;
417 } else if (i == 2) {
418 rotNames[resIndex] = rotName;
419 }
420 }
421 for (Atom atom : residue.getAtomList()) {
422 if (!residue.getBackboneAtoms().contains(atom) || diffStates) {
423 if (occupancy != 1) {
424 atom.setAltLoc(altLocs[confIndex]);
425 atom.setOccupancy(occupancy);
426 } else {
427 atom.setOccupancy(occupancy);
428 atom.setAltLoc(' ');
429 }
430 } else {
431 atom.setOccupancy(1.0);
432 atom.setAltLoc(' ');
433 }
434 }
435 }
436 }
437
438 conformerAssemblies[assemblyIndex] = conformerAssembly;
439 if (manyBodyOptions.getTitration()) {
440 titrationManyBody.excludeExcessAtoms(excludeAtoms, optimalRotamers, conformerAssemblies[assemblyIndex], conformerResidueList);
441 excludeAtoms.addAll(titrationManyBody.getExcludeAtoms());
442 }
443 assemblyIndex++;
444
445 if (pKa) {
446 rotamerOptimization.getProtonationPopulations(selectedResidues.toArray(new Residue[0]));
447 }
448 }
449
450
451 try (FileWriter popFileWriter = new FileWriter("populations.txt")) {
452 int residueIndex = 0;
453 for (Residue residue : selectedResidues) {
454 popFileWriter.write("\n");
455 protonationBoltzmannSums = new double[selectedResidues.size()];
456
457 Rotamer[] rotamers = residue.getRotamers();
458 for (int rotIndex = 0; rotIndex < rotamers.length; rotIndex++) {
459 String rotPop = format("%.6f", populationArray[residueIndex][rotIndex]);
460 popFileWriter.write(residue.getName() + residue.getResidueNumber() + "\t" +
461 rotamers[rotIndex].toString() + "\t" + rotPop + "\n");
462 }
463 residueIndex += 1;
464 }
465 } catch (IOException ioException) {
466 logger.warning("Error writing populations file: " + ioException.getMessage());
467 }
468
469 System.out.println("\n Successfully wrote to the populations file.");
470
471 PDBFilter pdbFilter = new PDBFilter(structureFile, Arrays.asList(conformerAssemblies), forceField, properties);
472 pdbFilter.writeFile(structureFile, false, excludeAtoms, true, true);
473
474 System.setProperty("standardizeAtomNames", "false");
475
476 return this;
477 }
478
479 @Override
480 public List<Potential> getPotentials() {
481 if (conformerAssemblies == null) {
482 return new ArrayList<>();
483 } else {
484 return Arrays.stream(conformerAssemblies)
485 .filter(a -> a != null)
486 .map(MolecularAssembly::getPotentialEnergy)
487 .filter(e -> e != null)
488 .collect(Collectors.toList());
489 }
490 }
491 }