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
159 System.setProperty("sor-scf-fallback", "false");
160 System.setProperty("direct-scf-fallback", "true");
161
162 xrayOptions.init();
163 double titrationPH = manyBodyOptions.getTitrationPH();
164 double inclusionCutoff = manyBodyOptions.getInclusionCutoff();
165 int mutatingResidue = manyBodyOptions.getInterestedResidue();
166 boolean onlyTitration = manyBodyOptions.getOnlyTitration();
167 double pHRestraint = manyBodyOptions.getPHRestraint();
168 if (manyBodyOptions.getTitration()) {
169 System.setProperty("manybody-titration", "true");
170 }
171
172
173 String nea = System.getProperty("native-environment-approximation", "true");
174 System.setProperty("native-environment-approximation", nea);
175
176 boolean lambdaTerm = alchemicalOptions.hasSoftcore();
177 if (lambdaTerm) {
178
179 System.setProperty("lambdaterm", "true");
180
181 System.setProperty("elec-lambdaterm", "false");
182
183 System.setProperty("intramolecular-softcore", "true");
184 }
185 System.setProperty("ro-ensembleEnergy", ensembleEnergy);
186
187 String modelFilename;
188 if (filenames != null && !filenames.isEmpty()) {
189 molecularAssemblies = algorithmFunctions.openAll(filenames.get(0));
190 activeAssembly = molecularAssemblies[0];
191 modelFilename = filenames.get(0);
192 } else if (activeAssembly == null) {
193 logger.info(helpString());
194 return this;
195 } else {
196 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
197 modelFilename = activeAssembly.getFile().getAbsolutePath();
198 }
199
200 CompositeConfiguration properties = activeAssembly.getProperties();
201 activeAssembly.getPotentialEnergy().setPrintOnFailure(false, false);
202 potentialEnergy = activeAssembly.getPotentialEnergy();
203 forceField = activeAssembly.getForceField();
204
205 if (forceField == null) {
206 logger.info("This force field is null");
207 }
208
209 String filename = filenames.get(0);
210
211 List<Residue> titrateResidues = new ArrayList<>();
212
213
214 Set<Atom> excludeAtoms = new HashSet<>();
215 boolean isTitrating = false;
216
217
218 if (xrayOptions.refinementMode != RefinementMode.COORDINATES) {
219 logger.info(" Refinement mode set to COORDINATES.");
220 xrayOptions.refinementMode = RefinementMode.COORDINATES;
221 }
222
223
224 List<Residue> residues = manyBodyOptions.collectResidues(activeAssembly);
225 if (residues == null || residues.isEmpty()) {
226 logger.info(" There are no residues in the active system to optimize.");
227 return this;
228 }
229
230 String[] titratableResidues = {"HIS", "HIE", "HID", "GLU", "GLH", "ASP", "ASH", "LYS", "LYD", "CYS", "CYD"};
231 List<String> titratableResiudesList = Arrays.asList(titratableResidues);
232 double boltzmannWeights;
233 double offsets;
234 double[][] populationArray = new double[1][55];
235 double[][] titrateBoltzmann;
236 double[] protonationBoltzmannSums;
237 double totalBoltzmann = 0;
238 List<Residue> residueList = activeAssembly.getResidueList();
239
240 String listResidues = "";
241
242
243 if (inclusionCutoff != -1 || onlyTitration) {
244 listResidues = manyBodyOptions.selectInclusionResidues(residueList, mutatingResidue, onlyTitration, inclusionCutoff);
245 }
246
247 List<Integer> residueNumber = new ArrayList<>();
248 for (Residue residue : residueList) {
249 residueNumber.add(residue.getResidueNumber());
250 }
251
252
253 if (properties.getBoolean("standardizeAtomNames", false)) {
254 renameAtomsToPDBStandard(activeAssembly);
255 }
256
257
258 if (manyBodyOptions.getTitration()) {
259 logger.info("\n Adding titration hydrogen to : " + filenames.get(0) + "\n");
260
261
262 List<Integer> resNumberList = new ArrayList<>();
263 for (Residue residue : residues) {
264 resNumberList.add(residue.getResidueNumber());
265 }
266
267
268 titrationManyBody = new TitrationManyBody(filename, activeAssembly.getForceField(),
269 resNumberList, titrationPH, manyBodyOptions);
270 activeAssembly = titrationManyBody.getProtonatedAssembly();
271 potentialEnergy = activeAssembly.getPotentialEnergy();
272 }
273
274
275 xrayOptions.setProperties(parseResult, properties);
276 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
277
278 DiffractionData diffractionData = xrayOptions.getDiffractionData(filenames, molecularAssemblies, properties);
279 refinementEnergy = xrayOptions.toXrayEnergy(diffractionData);
280 refinementEnergy.setScaling(null);
281
282
283 if (!pKa || onlyTitration) {
284 manyBodyOptions.setListResidues(listResidues);
285 }
286
287 if (lambdaTerm) {
288 alchemicalOptions.setFirstSystemAlchemistry(activeAssembly);
289 LambdaInterface lambdaInterface = (LambdaInterface) potentialEnergy;
290 double lambda = alchemicalOptions.getInitialLambda();
291 logger.info(format(" Setting ManyBody softcore lambda to: %5.3f", lambda));
292 lambdaInterface.setLambda(lambda);
293 }
294
295 RotamerOptimization rotamerOptimization = new RotamerOptimization(activeAssembly, refinementEnergy, algorithmListener);
296 manyBodyOptions.initRotamerOptimization(rotamerOptimization, activeAssembly);
297 rotamerOptimization.setPrintFiles(printFiles);
298 rotamerOptimization.setWriteEnergyRestart(printFiles);
299 rotamerOptimization.setPHRestraint(pHRestraint);
300 rotamerOptimization.setpH(titrationPH);
301 double[] x = new double[refinementEnergy.getNumberOfVariables()];
302 x = refinementEnergy.getCoordinates(x);
303 double e = refinementEnergy.energy(x, true);
304 logger.info(format("\n Initial target energy: %16.8f ", e));
305
306 selectedResidues = rotamerOptimization.getResidues();
307 rotamerOptimization.initFraction(selectedResidues);
308
309 RotamerLibrary.measureRotamers(selectedResidues, false);
310
311 rotamerOptimization.optimize(manyBodyOptions.getAlgorithm(selectedResidues.size()));
312
313 int[] currentRotamers = new int[selectedResidues.size()];
314
315
316 try {
317 rotamerOptimization.getFractions(selectedResidues.toArray(new Residue[0]), 0, currentRotamers);
318 if (pKa) {
319 rotamerOptimization.getProtonationPopulations(selectedResidues.toArray(new Residue[0]));
320 }
321 } catch (Exception ex) {
322 logger.severe(" Error calculating rotamer fractions: " + ex.getMessage());
323 return this;
324 }
325
326
327 boltzmannWeights = rotamerOptimization.getTotalBoltzmann();
328 offsets = rotamerOptimization.getRefEnergy();
329
330
331 populationArray = rotamerOptimization.getFraction();
332
333
334
335
336
337
338 try (FileWriter fileWriter = new FileWriter("populations.txt")) {
339 int residueIndex = 0;
340 for (Residue residue : selectedResidues) {
341 fileWriter.write("\n");
342 Rotamer[] rotamers = residue.getRotamers();
343 for (int rotIndex = 0; rotIndex < rotamers.length; rotIndex++) {
344 String rotPop = format("%.6f", populationArray[residueIndex][rotIndex]);
345 fileWriter.write(residue.getName() + residue.getResidueNumber() + "\t" +
346 rotamers[rotIndex].toString() + "\t" + rotPop + "\n");
347 }
348 residueIndex += 1;
349 }
350 } catch (IOException ioException) {
351 logger.warning("Error writing populations file: " + ioException.getMessage());
352 }
353
354 System.out.println("\n Successfully wrote to the populations file.");
355
356 int[][] conformers;
357 try {
358 conformers = rotamerOptimization.getConformers();
359 } catch (Exception ex) {
360 logger.severe(" Error getting conformers: " + ex.getMessage());
361 return this;
362 }
363
364 char[] altLocs = new char[]{'C', 'B', 'A'};
365 List<String> residueChainNum = new ArrayList<>();
366 for (Residue residue : activeAssembly.getResidueList()) {
367 char chain = residue.getChainID();
368 int resNum = residue.getResidueNumber();
369 residueChainNum.add(String.valueOf(chain) + resNum);
370 }
371
372 File structureFile = new File(filename);
373
374 int[] optimalRotamers = new int[selectedResidues.size()];
375 int assemblyIndex = 0;
376 String[] rotNames = new String[selectedResidues.size()];
377 for (int confIndex = 2; confIndex > -1; confIndex--) {
378 List<Residue> conformerResidueList = new ArrayList<>();
379 MolecularAssembly conformerAssembly = algorithmFunctions.open(filename);
380 if (manyBodyOptions.getTitration()) {
381 logger.info("\n Adding titration hydrogen to : " + filenames.get(0) + "\n");
382
383
384 List<Integer> resNumberList = new ArrayList<>();
385 for (Residue residue : residues) {
386 resNumberList.add(residue.getResidueNumber());
387 }
388
389
390 titrationManyBody = new TitrationManyBody(filename, activeAssembly.getForceField(),
391 resNumberList, titrationPH, manyBodyOptions);
392 conformerAssembly = titrationManyBody.getProtonatedAssembly();
393 potentialEnergy = conformerAssembly.getPotentialEnergy();
394 }
395 for (int resIndex = 0; resIndex < selectedResidues.size(); resIndex++) {
396 Residue residueSelect = selectedResidues.get(resIndex);
397 String resChainNum = String.valueOf(residueSelect.getChainID()) + residueSelect.getResidueNumber();
398 int index = residueChainNum.indexOf(resChainNum);
399 Residue residue = conformerAssembly.getResidueList().get(index);
400 conformerResidueList.add(residue);
401 residue.setRotamers(manyBodyOptions.getRotamerLibrary(true));
402 Rotamer[] rotamers = residue.getRotamers();
403 int rotIndex = conformers[resIndex][confIndex];
404 if (populationArray[resIndex][rotIndex] != 0) {
405 optimalRotamers[resIndex] = rotIndex;
406 RotamerLibrary.applyRotamer(residue, rotamers[rotIndex]);
407 double occupancy = populationArray[resIndex][rotIndex];
408 boolean diffStates = false;
409 for (int i = 2; i > -1; i--) {
410 int rotamerInd = conformers[resIndex][i];
411 String rotName = rotamers[rotamerInd].getName();
412 double occupancyTest = populationArray[resIndex][rotamerInd];
413 if (i == 2 && rotNames[resIndex] != null) {
414 rotNames[resIndex] = rotName;
415 } else if (i < 2 && occupancyTest != 0 && !rotNames[resIndex].contains(rotName)) {
416 diffStates = true;
417 String newString = rotNames[resIndex] + rotName;
418 logger.info(newString);
419 rotNames[resIndex] = newString;
420 } else if (i == 2) {
421 rotNames[resIndex] = rotName;
422 }
423 }
424 for (Atom atom : residue.getAtomList()) {
425 if (!residue.getBackboneAtoms().contains(atom) || diffStates) {
426 if (occupancy != 1) {
427 atom.setAltLoc(altLocs[confIndex]);
428 atom.setOccupancy(occupancy);
429 } else {
430 atom.setOccupancy(occupancy);
431 atom.setAltLoc(' ');
432 }
433 } else {
434 atom.setOccupancy(1.0);
435 atom.setAltLoc(' ');
436 }
437 }
438 }
439 }
440
441 conformerAssemblies[assemblyIndex] = conformerAssembly;
442 if (manyBodyOptions.getTitration()) {
443 titrationManyBody.excludeExcessAtoms(excludeAtoms, optimalRotamers, conformerAssemblies[assemblyIndex], conformerResidueList);
444 excludeAtoms.addAll(titrationManyBody.getExcludeAtoms());
445 }
446 assemblyIndex++;
447
448 if (pKa) {
449 rotamerOptimization.getProtonationPopulations(selectedResidues.toArray(new Residue[0]));
450 }
451 }
452
453
454 try (FileWriter popFileWriter = new FileWriter("populations.txt")) {
455 int residueIndex = 0;
456 for (Residue residue : selectedResidues) {
457 popFileWriter.write("\n");
458 protonationBoltzmannSums = new double[selectedResidues.size()];
459
460 Rotamer[] rotamers = residue.getRotamers();
461 for (int rotIndex = 0; rotIndex < rotamers.length; rotIndex++) {
462 String rotPop = format("%.6f", populationArray[residueIndex][rotIndex]);
463 popFileWriter.write(residue.getName() + residue.getResidueNumber() + "\t" +
464 rotamers[rotIndex].toString() + "\t" + rotPop + "\n");
465 }
466 residueIndex += 1;
467 }
468 } catch (IOException ioException) {
469 logger.warning("Error writing populations file: " + ioException.getMessage());
470 }
471
472 System.out.println("\n Successfully wrote to the populations file.");
473
474 PDBFilter pdbFilter = new PDBFilter(structureFile, Arrays.asList(conformerAssemblies), forceField, properties);
475 pdbFilter.writeFile(structureFile, false, excludeAtoms, true, true);
476
477 System.setProperty("standardizeAtomNames", "false");
478
479 return this;
480 }
481
482 @Override
483 public List<Potential> getPotentials() {
484 if (conformerAssemblies == null) {
485 return new ArrayList<>();
486 } else {
487 return Arrays.stream(conformerAssemblies)
488 .filter(a -> a != null)
489 .map(MolecularAssembly::getPotentialEnergy)
490 .filter(e -> e != null)
491 .collect(Collectors.toList());
492 }
493 }
494 }