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 MolecularAssembly protonatedAssembly = titrationManyBody.getProtonatedAssembly();
242 setActiveAssembly(protonatedAssembly);
243 potentialEnergy = protonatedAssembly.getPotentialEnergy();
244 }
245 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
246 refinementEnergy = realSpaceOptions.toRealSpaceEnergy(filenames, molecularAssemblies);
247
248 if (lambdaTerm) {
249 alchemicalOptions.setFirstSystemAlchemistry(activeAssembly);
250 LambdaInterface lambdaInterface = (LambdaInterface) potentialEnergy;
251 double lambda = alchemicalOptions.getInitialLambda();
252 logger.info(format(" Setting ManyBody softcore lambda to: %5.3f", lambda));
253 lambdaInterface.setLambda(lambda);
254 }
255
256 RotamerOptimization rotamerOptimization = new RotamerOptimization(activeAssembly, refinementEnergy, algorithmListener);
257 rotamerOptimization.setPrintFiles(printFiles);
258 rotamerOptimization.setWriteEnergyRestart(printFiles);
259 rotamerOptimization.setPHRestraint(pHRestraint);
260 rotamerOptimization.setpH(titrationPH);
261
262 manyBodyOptions.initRotamerOptimization(rotamerOptimization, activeAssembly);
263
264 double[] x = new double[refinementEnergy.getNumberOfVariables()];
265 x = refinementEnergy.getCoordinates(x);
266 double e = refinementEnergy.energy(x, true);
267 logger.info(format("\n Initial target energy: %16.8f ", e));
268
269 selectedResidues = rotamerOptimization.getResidues();
270 rotamerOptimization.initFraction(selectedResidues);
271
272 RotamerLibrary.measureRotamers(residueList, false);
273
274 rotamerOptimization.optimize(manyBodyOptions.getAlgorithm(residueList.size()));
275
276 int[] currentRotamers = new int[selectedResidues.size()];
277
278
279 try {
280 rotamerOptimization.getFractions(selectedResidues.toArray(new Residue[0]), 0, currentRotamers);
281 if (pKa) {
282 rotamerOptimization.getProtonationPopulations(selectedResidues.toArray(new Residue[0]));
283 }
284 } catch (Exception ex) {
285 logger.severe(" Error calculating rotamer fractions: " + ex.getMessage());
286 return this;
287 }
288
289
290 boltzmannWeights = rotamerOptimization.getTotalBoltzmann();
291 offsets = rotamerOptimization.getRefEnergy();
292
293
294 populationArray = rotamerOptimization.getFraction();
295
296 try (FileWriter fileWriter = new FileWriter("populations.txt")) {
297 int residueIndex = 0;
298 for (Residue residue : selectedResidues) {
299 fileWriter.write("\n");
300 Rotamer[] rotamers = residue.getRotamers();
301 for (int rotIndex = 0; rotIndex < rotamers.length; rotIndex++) {
302 String rotPop = format("%.6f", populationArray[residueIndex][rotIndex]);
303 fileWriter.write(residue.getName() + residue.getResidueNumber() + "\t" +
304 rotamers[rotIndex].toString() + "\t" + rotPop + "\n");
305 }
306 residueIndex += 1;
307 }
308 System.out.println("\n Successfully wrote to the populations file.");
309 } catch (IOException ioException) {
310 logger.severe(" Error writing populations file: " + ioException.getMessage());
311 }
312
313 int[][] conformers;
314 try {
315 conformers = rotamerOptimization.getConformers();
316 } catch (Exception ex) {
317 logger.severe(" Error getting conformers: " + ex.getMessage());
318 return this;
319 }
320
321 char[] altLocs = new char[]{'C', 'B', 'A'};
322 List<String> residueChainNum = new ArrayList<>();
323 for (Residue residue : activeAssembly.getResidueList()) {
324 char chain = residue.getChainID();
325 int resNum = residue.getResidueNumber();
326 residueChainNum.add(String.valueOf(chain) + String.valueOf(resNum));
327 }
328
329 File structureFile = new File(filename);
330 int[] optimalRotamers = new int[selectedResidues.size()];
331 int assemblyIndex = 0;
332 String[] rotNames = new String[selectedResidues.size()];
333 for (int confIndex = 2; confIndex > -1; confIndex--) {
334 List<Residue> conformerResidueList = new ArrayList<>();
335 MolecularAssembly conformerAssembly = algorithmFunctions.open(filename);
336 if (manyBodyOptions.getTitration()) {
337 logger.info("\n Adding titration hydrogen to : " + filenames.get(0) + "\n");
338
339
340 List<Integer> resNumberList = new ArrayList<>();
341 for (Residue residue : residues) {
342 resNumberList.add(residue.getResidueNumber());
343 }
344
345
346 titrationManyBody = new TitrationManyBody(filename, activeAssembly.getForceField(),
347 resNumberList, titrationPH, manyBodyOptions);
348 conformerAssembly = titrationManyBody.getProtonatedAssembly();
349 potentialEnergy = conformerAssembly.getPotentialEnergy();
350 }
351 for (int resIndex = 0; resIndex < selectedResidues.size(); resIndex++) {
352 Residue residueSelect = selectedResidues.get(resIndex);
353 String resChainNum = String.valueOf(residueSelect.getChainID()) + String.valueOf(residueSelect.getResidueNumber());
354 int index = residueChainNum.indexOf(resChainNum);
355 Residue residue = conformerAssembly.getResidueList().get(index);
356 conformerResidueList.add(residue);
357 residue.setRotamers(manyBodyOptions.getRotamerLibrary(true));
358 Rotamer[] rotamers = residue.getRotamers();
359 int rotIndex = conformers[resIndex][confIndex];
360 if (populationArray[resIndex][rotIndex] != 0) {
361 optimalRotamers[resIndex] = rotIndex;
362 RotamerLibrary.applyRotamer(residue, rotamers[rotIndex]);
363 double occupancy = populationArray[resIndex][rotIndex];
364 boolean diffStates = false;
365 for (int i = 2; i > -1; i--) {
366 int rotamerInd = conformers[resIndex][i];
367 String rotName = rotamers[rotamerInd].getName();
368 double occupancyTest = populationArray[resIndex][rotamerInd];
369 if (i == 2 && rotNames[resIndex] != null) {
370 rotNames[resIndex] = rotName;
371 } else if (i < 2 && occupancyTest != 0 && !rotNames[resIndex].contains(rotName)) {
372 diffStates = true;
373 String newString = rotNames[resIndex] + rotName;
374 logger.info(newString);
375 rotNames[resIndex] = newString;
376 } else if (i == 2) {
377 rotNames[resIndex] = rotName;
378 }
379 }
380 for (Atom atom : residue.getAtomList()) {
381 if (!residue.getBackboneAtoms().contains(atom) || diffStates) {
382 if (occupancy != 1) {
383 atom.setAltLoc(altLocs[confIndex]);
384 atom.setOccupancy(occupancy);
385 } else {
386 atom.setOccupancy(occupancy);
387 atom.setAltLoc(' ');
388 }
389 } else {
390 atom.setOccupancy(1.0);
391 atom.setAltLoc(' ');
392 }
393 }
394 }
395 }
396
397 conformerAssemblies[assemblyIndex] = conformerAssembly;
398 titrationManyBody.excludeExcessAtoms(excludeAtoms, optimalRotamers, conformerAssemblies[assemblyIndex], conformerResidueList);
399 excludeAtoms.addAll(titrationManyBody.getExcludeAtoms());
400 assemblyIndex++;
401 }
402 PDBFilter pdbFilter = new PDBFilter(structureFile, Arrays.asList(conformerAssemblies), forceField, properties);
403 pdbFilter.writeFile(structureFile, false, excludeAtoms, true, true);
404 System.setProperty("standardizeAtomNames", "false");
405
406 return this;
407 }
408 }