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