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 edu.rit.pj.Comm;
41 import ffx.algorithms.cli.AlgorithmsCommand;
42 import ffx.algorithms.cli.ManyBodyOptions;
43 import ffx.algorithms.optimize.RotamerOptimization;
44 import ffx.algorithms.optimize.TitrationManyBody;
45 import ffx.numerics.Potential;
46 import ffx.potential.ForceFieldEnergy;
47 import ffx.potential.MolecularAssembly;
48 import ffx.potential.bonded.Atom;
49 import ffx.potential.bonded.Residue;
50 import ffx.potential.bonded.RotamerLibrary;
51 import ffx.potential.parsers.PDBFilter;
52 import ffx.realspace.cli.RealSpaceOptions;
53 import ffx.utilities.FFXBinding;
54 import ffx.xray.RefinementEnergy;
55 import org.apache.commons.configuration2.CompositeConfiguration;
56 import org.apache.commons.io.FilenameUtils;
57 import picocli.CommandLine.Command;
58 import picocli.CommandLine.Mixin;
59 import picocli.CommandLine.Parameters;
60
61 import java.io.File;
62 import java.util.ArrayList;
63 import java.util.Collections;
64 import java.util.HashSet;
65 import java.util.List;
66 import java.util.Set;
67
68 import static java.lang.String.format;
69
70
71
72
73
74
75
76
77 @Command(description = " Discrete optimization using a many-body expansion and elimination expressions.", name = "realspace.ManyBody")
78 public class ManyBody extends AlgorithmsCommand {
79
80 @Mixin
81 private RealSpaceOptions realSpaceOptions;
82
83 @Mixin
84 private ManyBodyOptions manyBodyOptions;
85
86
87
88
89 @Parameters(arity = "1..*", paramLabel = "files", description = "PDB and Real Space input files.")
90 private List<String> filenames;
91
92 private RefinementEnergy refinementEnergy;
93
94 ForceFieldEnergy potentialEnergy;
95 TitrationManyBody titrationManyBody;
96
97
98
99
100 public ManyBody() {
101 super();
102 }
103
104
105
106
107
108 public ManyBody(String[] args) {
109 super(args);
110 }
111
112
113
114
115
116 public ManyBody(FFXBinding binding) {
117 super(binding);
118 }
119
120 @Override
121 public ManyBody run() {
122
123 if (!init()) {
124 return this;
125 }
126
127
128
129
130
131
132 double titrationPH = manyBodyOptions.getTitrationPH();
133 if (titrationPH > 0) {
134 System.setProperty("manybody-titration", "true");
135 }
136
137 String filename;
138 if (filenames != null && !filenames.isEmpty()) {
139 activeAssembly = algorithmFunctions.open(filenames.get(0));
140 filename = filenames.get(0);
141 } else if (activeAssembly == null) {
142 logger.info(helpString());
143 return this;
144 } else {
145 filename = activeAssembly.getFile().getAbsolutePath();
146 }
147 MolecularAssembly[] molecularAssemblies = new MolecularAssembly[]{activeAssembly};
148
149 CompositeConfiguration properties = activeAssembly.getProperties();
150 activeAssembly.getPotentialEnergy().setPrintOnFailure(false, false);
151 potentialEnergy = activeAssembly.getPotentialEnergy();
152
153
154 List<Residue> residues = manyBodyOptions.collectResidues(activeAssembly);
155 if (residues == null || residues.isEmpty()) {
156 logger.info(" There are no residues in the active system to optimize.");
157 return this;
158 }
159
160
161 if (titrationPH > 0) {
162 logger.info("\n Adding titration hydrogen to : " + filenames.get(0) + "\n");
163
164 List<Integer> resNumberList = new ArrayList<>();
165 for (Residue residue : residues) {
166 resNumberList.add(residue.getResidueNumber());
167 }
168
169
170 titrationManyBody = new TitrationManyBody(filenames.getFirst(), activeAssembly.getForceField(),
171 resNumberList, titrationPH, manyBodyOptions);
172 activeAssembly = titrationManyBody.getProtonatedAssembly();
173 potentialEnergy = activeAssembly.getPotentialEnergy();
174 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
175 }
176
177 refinementEnergy = realSpaceOptions.toRealSpaceEnergy(filenames, molecularAssemblies);
178
179 RotamerOptimization rotamerOptimization = new RotamerOptimization(
180 activeAssembly, refinementEnergy, algorithmListener);
181 manyBodyOptions.initRotamerOptimization(rotamerOptimization, activeAssembly);
182
183 double[] x = new double[refinementEnergy.getNumberOfVariables()];
184 x = refinementEnergy.getCoordinates(x);
185 refinementEnergy.energy(x, true);
186
187 List<Residue> residueList = rotamerOptimization.getResidues();
188 RotamerLibrary.measureRotamers(residueList, false);
189
190 rotamerOptimization.optimize(manyBodyOptions.getAlgorithm(residueList.size()));
191
192 boolean isTitrating = false;
193 Set<Atom> excludeAtoms = new HashSet<>();
194 int[] optimalRotamers = rotamerOptimization.getOptimumRotamers();
195
196 if (titrationPH > 0) {
197 isTitrating = titrationManyBody.excludeExcessAtoms(excludeAtoms, optimalRotamers, residueList);
198 }
199
200 if (Comm.world().rank() == 0) {
201 logger.info(" Final Minimum Energy");
202
203 x = refinementEnergy.getCoordinates(x);
204 double energy = refinementEnergy.energy(x, true);
205
206 if (isTitrating) {
207 double phBias = rotamerOptimization.getEnergyExpansion().getTotalRotamerPhBias(residueList,
208 optimalRotamers, titrationPH, manyBodyOptions.getPHRestraint());
209 logger.info(format("\n Rotamer pH Bias %16.8f", phBias));
210 logger.info(format(" Potential with Bias%16.8f\n", phBias + energy));
211 } else {
212 logger.info(format("\n Real Space Target %16.8f\n", energy));
213 }
214 String ext = FilenameUtils.getExtension(filename);
215 filename = FilenameUtils.removeExtension(filename);
216 if (ext.toUpperCase().contains("XYZ")) {
217 algorithmFunctions.saveAsXYZ(molecularAssemblies[0], new File(filename + ".xyz"));
218 } else {
219 properties.setProperty("standardizeAtomNames", "false");
220 File modelFile = saveDirFile(activeAssembly.getFile());
221 PDBFilter pdbFilter = new PDBFilter(modelFile, activeAssembly,
222 activeAssembly.getForceField(),
223 properties);
224 if (titrationPH > 0) {
225 String remark = format("Titration pH: %6.3f", titrationPH);
226 if (!pdbFilter.writeFile(modelFile, false, excludeAtoms, true, true, new String[]{remark})) {
227 logger.info(format(" Save failed for %s", activeAssembly));
228 }
229 } else {
230 if (!pdbFilter.writeFile(modelFile, false, excludeAtoms, true, true)) {
231 logger.info(format(" Save failed for %s", activeAssembly));
232 }
233 }
234 }
235 }
236
237 return this;
238 }
239
240
241
242
243
244 public ManyBodyOptions getManyBodyOptions() {
245 return manyBodyOptions;
246 }
247
248 @Override
249 public List<Potential> getPotentials() {
250 return refinementEnergy == null ? Collections.emptyList() :
251 Collections.singletonList((Potential) refinementEnergy);
252 }
253 }