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.get(0), activeAssembly.getForceField(),
171 resNumberList, titrationPH, manyBodyOptions);
172 MolecularAssembly protonatedAssembly = titrationManyBody.getProtonatedAssembly();
173 setActiveAssembly(protonatedAssembly);
174 potentialEnergy = protonatedAssembly.getPotentialEnergy();
175 molecularAssemblies = new MolecularAssembly[]{activeAssembly};
176 }
177
178 refinementEnergy = realSpaceOptions.toRealSpaceEnergy(filenames, molecularAssemblies);
179
180 RotamerOptimization rotamerOptimization = new RotamerOptimization(
181 activeAssembly, refinementEnergy, algorithmListener);
182 manyBodyOptions.initRotamerOptimization(rotamerOptimization, activeAssembly);
183
184 double[] x = new double[refinementEnergy.getNumberOfVariables()];
185 x = refinementEnergy.getCoordinates(x);
186 refinementEnergy.energy(x, true);
187
188 List<Residue> residueList = rotamerOptimization.getResidues();
189 RotamerLibrary.measureRotamers(residueList, false);
190
191 rotamerOptimization.optimize(manyBodyOptions.getAlgorithm(residueList.size()));
192
193 boolean isTitrating = false;
194 Set<Atom> excludeAtoms = new HashSet<>();
195 int[] optimalRotamers = rotamerOptimization.getOptimumRotamers();
196
197 if (titrationPH > 0) {
198 isTitrating = titrationManyBody.excludeExcessAtoms(excludeAtoms, optimalRotamers, residueList);
199 }
200
201 if (Comm.world().rank() == 0) {
202 logger.info(" Final Minimum Energy");
203
204 x = refinementEnergy.getCoordinates(x);
205 double energy = refinementEnergy.energy(x, true);
206
207 if (isTitrating) {
208 double phBias = rotamerOptimization.getEnergyExpansion().getTotalRotamerPhBias(residueList,
209 optimalRotamers, titrationPH, manyBodyOptions.getPHRestraint());
210 logger.info(format("\n Rotamer pH Bias %16.8f", phBias));
211 logger.info(format(" Potential with Bias%16.8f\n", phBias + energy));
212 } else {
213 logger.info(format("\n Real Space Target %16.8f\n", energy));
214 }
215 String ext = FilenameUtils.getExtension(filename);
216 filename = FilenameUtils.removeExtension(filename);
217 if (ext.toUpperCase().contains("XYZ")) {
218 algorithmFunctions.saveAsXYZ(molecularAssemblies[0], new File(filename + ".xyz"));
219 } else {
220 properties.setProperty("standardizeAtomNames", "false");
221 File modelFile = saveDirFile(activeAssembly.getFile());
222 PDBFilter pdbFilter = new PDBFilter(modelFile, activeAssembly,
223 activeAssembly.getForceField(),
224 properties);
225 if (titrationPH > 0) {
226 String remark = format("Titration pH: %6.3f", titrationPH);
227 if (!pdbFilter.writeFile(modelFile, false, excludeAtoms, true, true, new String[]{remark})) {
228 logger.info(format(" Save failed for %s", activeAssembly));
229 }
230 } else {
231 if (!pdbFilter.writeFile(modelFile, false, excludeAtoms, true, true)) {
232 logger.info(format(" Save failed for %s", activeAssembly));
233 }
234 }
235 }
236 }
237
238 return this;
239 }
240
241
242
243
244
245 public ManyBodyOptions getManyBodyOptions() {
246 return manyBodyOptions;
247 }
248
249 @Override
250 public List<Potential> getPotentials() {
251 return refinementEnergy == null ? Collections.emptyList() :
252 Collections.singletonList((Potential) refinementEnergy);
253 }
254 }