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.algorithms.commands;
39
40 import ffx.algorithms.cli.AlgorithmsCommand;
41 import ffx.algorithms.cli.MinimizeOptions;
42 import ffx.crystal.Crystal;
43 import ffx.numerics.Potential;
44 import ffx.potential.ForceFieldEnergy;
45 import ffx.potential.MolecularAssembly;
46 import ffx.potential.bonded.LambdaInterface;
47 import ffx.potential.cli.AlchemicalOptions;
48 import ffx.potential.cli.AtomSelectionOptions;
49 import ffx.potential.cli.TopologyOptions;
50 import ffx.potential.parsers.PDBFilter;
51 import ffx.potential.parsers.SystemFilter;
52 import ffx.potential.parsers.XYZFilter;
53 import ffx.utilities.FFXBinding;
54 import ffx.utilities.FileUtils;
55 import org.apache.commons.io.FilenameUtils;
56 import picocli.CommandLine.Command;
57 import picocli.CommandLine.Mixin;
58 import picocli.CommandLine.Parameters;
59
60 import java.io.File;
61 import java.util.ArrayList;
62 import java.util.Collections;
63 import java.util.List;
64
65 import static java.lang.String.format;
66
67
68
69
70
71
72
73
74 @Command(description = " Run L-BFGS minimization on a system.", name = "Minimize")
75 public class Minimize extends AlgorithmsCommand {
76
77 @Mixin
78 private MinimizeOptions minimizeOptions;
79
80 @Mixin
81 private AtomSelectionOptions atomSelectionOptions;
82
83 @Mixin
84 private AlchemicalOptions alchemicalOptions;
85
86 @Mixin
87 private TopologyOptions topologyOptions;
88
89
90
91
92 @Parameters(arity = "1..*", paramLabel = "files", description = "Atomic coordinate files in PDB or XYZ format.")
93 private List<String> filenames = null;
94
95 private MolecularAssembly[] topologies;
96 private Potential potential;
97
98
99
100
101 private ffx.algorithms.optimize.Minimize minimize;
102
103
104
105
106 public Minimize() {
107 super();
108 }
109
110
111
112
113
114
115 public Minimize(FFXBinding binding) {
116 super(binding);
117 }
118
119
120
121
122
123
124 public Minimize(String[] args) {
125 super(args);
126 }
127
128
129
130
131 @Override
132 public Minimize run() {
133
134 if (!init()) {
135 return this;
136 }
137
138
139 int numTopologies = topologyOptions.getNumberOfTopologies(filenames);
140 int threadsPerTopology = topologyOptions.getThreadsPerTopology(numTopologies);
141 topologies = new MolecularAssembly[numTopologies];
142
143
144 alchemicalOptions.setAlchemicalProperties();
145 topologyOptions.setAlchemicalProperties(numTopologies);
146
147
148 if (filenames == null || filenames.isEmpty()) {
149 activeAssembly = getActiveAssembly(null);
150 if (activeAssembly == null) {
151 logger.info(helpString());
152 return this;
153 }
154 filenames = new ArrayList<>();
155 filenames.add(activeAssembly.getFile().getName());
156 topologies[0] = alchemicalOptions.processFile(topologyOptions, activeAssembly, 0);
157 } else {
158 logger.info(format(" Initializing %d topologies...", numTopologies));
159 for (int i = 0; i < numTopologies; i++) {
160 topologies[i] = alchemicalOptions.openFile(algorithmFunctions,
161 topologyOptions, threadsPerTopology, filenames.get(i), i);
162 }
163 activeAssembly = topologies[0];
164 }
165
166 if (topologies.length == 1) {
167 atomSelectionOptions.setActiveAtoms(topologies[0]);
168 }
169
170
171 StringBuilder sb = new StringBuilder("\n Minimizing energy of ");
172 potential = topologyOptions.assemblePotential(topologies, sb);
173 logger.info(sb.toString());
174
175 LambdaInterface linter = (potential instanceof LambdaInterface) ? (LambdaInterface) potential : null;
176 if (linter != null) {
177 linter.setLambda(alchemicalOptions.getInitialLambda());
178 }
179
180 SystemFilter systemFilter = algorithmFunctions.getFilter();
181
182 double[] x = new double[potential.getNumberOfVariables()];
183 potential.getCoordinates(x);
184 potential.energy(x, true);
185
186 minimize = new ffx.algorithms.optimize.Minimize(topologies[0], potential, algorithmListener);
187 minimize.minimize(minimizeOptions.getNBFGS(), minimizeOptions.getEps(), minimizeOptions.getIterations());
188
189 potential.getCoordinates(x);
190 activeAssembly = systemFilter.getActiveMolecularSystem();
191 updateTitle(potential.energy(x, true));
192 if (topologies.length > 1) {
193
194 for (MolecularAssembly molecularAssembly : topologies) {
195 String modelFilename = molecularAssembly.getFile().getAbsolutePath();
196
197 if (baseDir == null || !baseDir.exists() || !baseDir.isDirectory() || !baseDir.canWrite()) {
198 baseDir = new File(FilenameUtils.getFullPath(modelFilename));
199 }
200
201 String dirName = baseDir.toString() + File.separator;
202 String fileName = FilenameUtils.getName(modelFilename);
203 String ext = FilenameUtils.getExtension(fileName);
204 fileName = FilenameUtils.removeExtension(fileName);
205
206 if (ext.toUpperCase().contains("XYZ")) {
207 algorithmFunctions.saveAsXYZ(molecularAssembly, new File(dirName + fileName + ".xyz"));
208 } else {
209 algorithmFunctions.saveAsPDB(molecularAssembly, new File(dirName + fileName + ".pdb"));
210 }
211 }
212 } else {
213
214 setActiveAssembly(topologies[0]);
215 String modelFilename = activeAssembly.getFile().getAbsolutePath();
216 if (baseDir == null || !baseDir.exists() || !baseDir.isDirectory() || !baseDir.canWrite()) {
217 baseDir = new File(FilenameUtils.getFullPath(modelFilename));
218 }
219
220 String dirName = baseDir.toString() + File.separator;
221 String fileName = FilenameUtils.getName(modelFilename);
222 String ext = FilenameUtils.getExtension(fileName);
223 fileName = FilenameUtils.removeExtension(fileName);
224 File saveFile;
225 SystemFilter writeFilter;
226 if (ext.toUpperCase().contains("XYZ")) {
227 saveFile = new File(dirName + fileName + ".xyz");
228 writeFilter = new XYZFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
229 activeAssembly.getProperties());
230 algorithmFunctions.saveAsXYZ(activeAssembly, saveFile);
231 } else if (ext.toUpperCase().contains("ARC")) {
232 saveFile = new File(dirName + fileName + ".arc");
233 saveFile = algorithmFunctions.versionFile(saveFile);
234 writeFilter = new XYZFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
235 activeAssembly.getProperties());
236 algorithmFunctions.saveAsXYZ(activeAssembly, saveFile);
237 } else {
238 saveFile = new File(dirName + fileName + ".pdb");
239 saveFile = algorithmFunctions.versionFile(saveFile);
240 PDBFilter pdbFilter = new PDBFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
241 activeAssembly.getProperties());
242 writeFilter = pdbFilter;
243 int numModels = systemFilter.countNumModels();
244 if (numModels > 1) {
245 pdbFilter.setModelNumbering(0);
246 }
247 pdbFilter.writeFile(saveFile, true, false, false);
248 }
249
250 if (systemFilter instanceof XYZFilter || systemFilter instanceof PDBFilter) {
251 while (systemFilter.readNext()) {
252 Crystal crystal = activeAssembly.getCrystal();
253 ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy();
254 forceFieldEnergy.setCrystal(crystal);
255 if (systemFilter instanceof PDBFilter) {
256 FileUtils.append(saveFile, "ENDMDL\n");
257 minimize.minimize(minimizeOptions.getEps(), minimizeOptions.getIterations());
258 PDBFilter pdbFilter = (PDBFilter) systemFilter;
259 pdbFilter.writeFile(saveFile, true, false, false);
260 } else if (systemFilter instanceof XYZFilter) {
261 minimize.minimize(minimizeOptions.getEps(), minimizeOptions.getIterations());
262 writeFilter.writeFile(saveFile, true);
263 }
264 }
265 if (systemFilter instanceof PDBFilter) {
266 FileUtils.append(saveFile, "END\n");
267 }
268 }
269 }
270
271 return this;
272 }
273
274
275
276
277 @Override
278 public String toString() {
279 StringBuilder sb = new StringBuilder();
280 if (minimize != null) {
281 sb.append(minimize.toString());
282 } else {
283 sb.append("No minimization has been performed.");
284 }
285 return sb.toString();
286 }
287
288
289
290
291
292
293 public List<Double> getEnergyList() {
294 if (minimize != null) {
295 return minimize.getEnergyList();
296 }
297 return null;
298 }
299
300
301
302
303
304
305 public double getRMSGradient() {
306 if (minimize != null) {
307 return minimize.getRMSGradient();
308 }
309 return Double.NaN;
310 }
311
312
313
314
315
316
317 public double getEnergy() {
318 if (minimize != null) {
319 return minimize.getEnergy();
320 }
321 return Double.NaN;
322 }
323
324
325
326
327 @Override
328 public List<Potential> getPotentials() {
329 List<Potential> potentials;
330 if (potential == null) {
331 potentials = Collections.emptyList();
332 } else {
333 potentials = Collections.singletonList(potential);
334 }
335 return potentials;
336 }
337
338 }