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 System.setProperty("sor-scf-fallback", "false");
129 System.setProperty("direct-scf-fallback", "true");
130
131
132
133
134
135
136 double titrationPH = manyBodyOptions.getTitrationPH();
137 if (titrationPH > 0) {
138 System.setProperty("manybody-titration", "true");
139 }
140
141 String filename;
142 if (filenames != null && !filenames.isEmpty()) {
143 activeAssembly = algorithmFunctions.open(filenames.get(0));
144 filename = filenames.get(0);
145 } else if (activeAssembly == null) {
146 logger.info(helpString());
147 return this;
148 } else {
149 filename = activeAssembly.getFile().getAbsolutePath();
150 }
151 MolecularAssembly[] molecularAssemblies = new MolecularAssembly[]{activeAssembly};
152
153 CompositeConfiguration properties = activeAssembly.getProperties();
154 activeAssembly.getPotentialEnergy().setPrintOnFailure(false, false);
155 potentialEnergy = activeAssembly.getPotentialEnergy();
156
157
158 List<Residue> residues = manyBodyOptions.collectResidues(activeAssembly);
159 if (residues == null || residues.isEmpty()) {
160 logger.info(" There are no residues in the active system to optimize.");
161 return this;
162 }
163
164
165 if (titrationPH > 0) {
166 logger.info("\n Adding titration hydrogen to : " + filenames.get(0) + "\n");
167
168 List<Integer> resNumberList = new ArrayList<>();
169 for (Residue residue : residues) {
170 resNumberList.add(residue.getResidueNumber());
171 }
172
173
174 titrationManyBody = new TitrationManyBody(filenames.getFirst(), activeAssembly.getForceField(),
175 resNumberList, titrationPH, manyBodyOptions);
176 activeAssembly = titrationManyBody.getProtonatedAssembly();
177 potentialEnergy = activeAssembly.getPotentialEnergy();
178 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
179 }
180
181 refinementEnergy = realSpaceOptions.toRealSpaceEnergy(filenames, molecularAssemblies);
182
183 RotamerOptimization rotamerOptimization = new RotamerOptimization(
184 activeAssembly, refinementEnergy, algorithmListener);
185 manyBodyOptions.initRotamerOptimization(rotamerOptimization, activeAssembly);
186
187 double[] x = new double[refinementEnergy.getNumberOfVariables()];
188 x = refinementEnergy.getCoordinates(x);
189 refinementEnergy.energy(x, true);
190
191 List<Residue> residueList = rotamerOptimization.getResidues();
192 RotamerLibrary.measureRotamers(residueList, false);
193
194 rotamerOptimization.optimize(manyBodyOptions.getAlgorithm(residueList.size()));
195
196 boolean isTitrating = false;
197 Set<Atom> excludeAtoms = new HashSet<>();
198 int[] optimalRotamers = rotamerOptimization.getOptimumRotamers();
199
200 if (titrationPH > 0) {
201 isTitrating = titrationManyBody.excludeExcessAtoms(excludeAtoms, optimalRotamers, residueList);
202 }
203
204 if (Comm.world().rank() == 0) {
205 logger.info(" Final Minimum Energy");
206
207 x = refinementEnergy.getCoordinates(x);
208 double energy = refinementEnergy.energy(x, true);
209
210 if (isTitrating) {
211 double phBias = rotamerOptimization.getEnergyExpansion().getTotalRotamerPhBias(residueList,
212 optimalRotamers, titrationPH, manyBodyOptions.getPHRestraint());
213 logger.info(format("\n Rotamer pH Bias %16.8f", phBias));
214 logger.info(format(" Potential with Bias%16.8f\n", phBias + energy));
215 } else {
216 logger.info(format("\n Real Space Target %16.8f\n", energy));
217 }
218 String ext = FilenameUtils.getExtension(filename);
219 filename = FilenameUtils.removeExtension(filename);
220 if (ext.toUpperCase().contains("XYZ")) {
221 algorithmFunctions.saveAsXYZ(molecularAssemblies[0], new File(filename + ".xyz"));
222 } else {
223 properties.setProperty("standardizeAtomNames", "false");
224 File modelFile = saveDirFile(activeAssembly.getFile());
225 PDBFilter pdbFilter = new PDBFilter(modelFile, activeAssembly,
226 activeAssembly.getForceField(),
227 properties);
228 if (titrationPH > 0) {
229 String remark = format("Titration pH: %6.3f", titrationPH);
230 if (!pdbFilter.writeFile(modelFile, false, excludeAtoms, true, true, new String[]{remark})) {
231 logger.info(format(" Save failed for %s", activeAssembly));
232 }
233 } else {
234 if (!pdbFilter.writeFile(modelFile, false, excludeAtoms, true, true)) {
235 logger.info(format(" Save failed for %s", activeAssembly));
236 }
237 }
238 }
239 }
240
241 return this;
242 }
243
244
245
246
247
248 public ManyBodyOptions getManyBodyOptions() {
249 return manyBodyOptions;
250 }
251
252 @Override
253 public List<Potential> getPotentials() {
254 return refinementEnergy == null ? Collections.emptyList() :
255 Collections.singletonList((Potential) refinementEnergy);
256 }
257 }