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