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.test;
39
40 import ffx.algorithms.cli.AlgorithmsCommand;
41 import ffx.algorithms.cli.MinimizeOptions;
42 import ffx.algorithms.optimize.PhMinimize;
43 import ffx.crystal.Crystal;
44 import ffx.potential.ForceFieldEnergy;
45 import ffx.potential.extended.ExtendedSystem;
46 import ffx.potential.parsers.PDBFilter;
47 import ffx.potential.parsers.SystemFilter;
48 import ffx.potential.parsers.XPHFilter;
49 import ffx.potential.parsers.XYZFilter;
50 import ffx.utilities.FFXBinding;
51 import ffx.utilities.FileUtils;
52 import org.apache.commons.io.FilenameUtils;
53 import picocli.CommandLine.Command;
54 import picocli.CommandLine.Mixin;
55 import picocli.CommandLine.Option;
56 import picocli.CommandLine.Parameters;
57
58 import java.io.File;
59
60 import static java.lang.String.format;
61 import static org.apache.commons.math3.util.FastMath.abs;
62
63
64
65
66
67
68
69
70 @Command(description = " Run L-BFGS minimization on a CpHMD system.", name = "test.MinimizePh")
71 public class MinimizePh extends AlgorithmsCommand {
72
73 @Mixin
74 private MinimizeOptions minimizeOptions;
75
76
77
78
79 @Option(names = {"--pH", "--constantPH"}, paramLabel = "7.4", defaultValue = "7.4",
80 description = "pH value for the energy evaluation. (Only applies when esvTerm is true)")
81 private double pH;
82
83
84
85
86 @Option(names = {"--coords"}, paramLabel = "false", defaultValue = "false",
87 description = "Minimize spatial coordinates along with titration")
88 private boolean coords;
89
90
91
92
93 @Parameters(arity = "1..*", paramLabel = "files", description = "Atomic coordinate files in PDB or XYZ format.")
94 private String filename;
95
96 private ForceFieldEnergy forceFieldEnergy;
97
98
99
100
101 public MinimizePh() {
102 super();
103 }
104
105
106
107
108
109
110 public MinimizePh(FFXBinding binding) {
111 super(binding);
112 }
113
114
115
116
117
118
119 public MinimizePh(String[] args) {
120 super(args);
121 }
122
123
124
125
126 @Override
127 public MinimizePh run() {
128
129
130 if (!init()) {
131 return this;
132 }
133
134
135 activeAssembly = getActiveAssembly(filename);
136 if (activeAssembly == null) {
137 logger.info(helpString());
138 return this;
139 }
140
141
142 filename = activeAssembly.getFile().getAbsolutePath();
143
144 logger.info("\n Running MinimizePh on " + filename);
145 forceFieldEnergy = activeAssembly.getPotentialEnergy();
146
147 ExtendedSystem esvSystem = new ExtendedSystem(activeAssembly, pH, null);
148 esvSystem.setConstantPh(pH);
149 int numESVs = esvSystem.getExtendedResidueList().size();
150 forceFieldEnergy.attachExtendedSystem(esvSystem);
151 logger.info(format(" Attached extended system with %d residues.", numESVs));
152
153 double[] x = new double[forceFieldEnergy.getNumberOfVariables()];
154 forceFieldEnergy.getCoordinates(x);
155 forceFieldEnergy.energy(x, true);
156
157 SystemFilter systemFilter = algorithmFunctions.getFilter();
158 if (systemFilter instanceof XYZFilter) {
159 XPHFilter xphFilter = new XPHFilter(activeAssembly.getFile(), activeAssembly,
160 activeAssembly.getForceField(), activeAssembly.getProperties(), esvSystem);
161 xphFilter.readFile();
162 logger.info(" Reading ESV lambdas from XPH file");
163 forceFieldEnergy.getCoordinates(x);
164 forceFieldEnergy.energy(x, true);
165 }
166 PhMinimize minimize = new PhMinimize(activeAssembly, forceFieldEnergy, algorithmListener,
167 esvSystem);
168
169 double energy = minimize.getEnergy();
170 double tolerance = 1.0e-10;
171
172 if (coords) {
173 while (true) {
174
175 minimize.minimizeCoordinates(minimizeOptions.getNBFGS(),
176 minimizeOptions.getEps(), minimizeOptions.getIterations());
177 double newEnergy = minimize.getEnergy();
178 int status = minimize.getStatus();
179 if (status != 0) {
180 break;
181 }
182 if (abs(newEnergy - energy) <= tolerance) {
183 break;
184 }
185 energy = newEnergy;
186
187
188 minimize.minimizeTitration(minimizeOptions.getNBFGS(),
189 minimizeOptions.getEps(), minimizeOptions.getIterations());
190 newEnergy = minimize.getEnergy();
191 status = minimize.getStatus();
192 if (status != 0) {
193 break;
194 }
195 if (abs(newEnergy - energy) <= tolerance) {
196 break;
197 }
198 energy = newEnergy;
199 }
200 } else {
201 minimize.minimizeTitration(minimizeOptions.getNBFGS(),
202 minimizeOptions.getEps(), minimizeOptions.getIterations());
203 }
204
205 forceFieldEnergy.getCoordinates(x);
206 forceFieldEnergy.energy(x, true);
207 esvSystem.writeRestart();
208
209
210 String modelFilename = activeAssembly.getFile().getAbsolutePath();
211 if (baseDir == null || !baseDir.exists() || !baseDir.isDirectory() || !baseDir.canWrite()) {
212 baseDir = new File(FilenameUtils.getFullPath(modelFilename));
213 }
214
215 String dirName = baseDir.toString() + File.separator;
216 String fileName = FilenameUtils.getName(modelFilename);
217 String ext = FilenameUtils.getExtension(fileName);
218 fileName = FilenameUtils.removeExtension(fileName);
219 File saveFile = null;
220 SystemFilter writeFilter = null;
221 if (ext.toUpperCase().contains("XYZ")) {
222 saveFile = new File(dirName + fileName + ".xyz");
223 writeFilter = new XYZFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
224 activeAssembly.getProperties());
225 algorithmFunctions.saveAsXYZ(activeAssembly, saveFile);
226 } else if (ext.toUpperCase().contains("ARC")) {
227 saveFile = new File(dirName + fileName + ".arc");
228 saveFile = algorithmFunctions.versionFile(saveFile);
229 writeFilter = new XYZFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
230 activeAssembly.getProperties());
231 algorithmFunctions.saveAsXYZ(activeAssembly, saveFile);
232 } else {
233 saveFile = new File(dirName + fileName + ".pdb");
234 saveFile = algorithmFunctions.versionFile(saveFile);
235
236 PDBFilter pdbFilter = new PDBFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
237 activeAssembly.getProperties());
238 writeFilter = pdbFilter;
239 int numModels = systemFilter.countNumModels();
240 if (numModels > 1) {
241 pdbFilter.setModelNumbering(0);
242 }
243 pdbFilter.writeFile(saveFile, true, false, false);
244 }
245
246 if (systemFilter instanceof XYZFilter || systemFilter instanceof PDBFilter) {
247 while (systemFilter.readNext()) {
248 Crystal crystal = activeAssembly.getCrystal();
249 ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy();
250 forceFieldEnergy.setCrystal(crystal);
251 if (systemFilter instanceof PDBFilter) {
252 FileUtils.append(saveFile, "ENDMDL\n");
253 minimize.minimizeTitration(minimizeOptions.getNBFGS(), minimizeOptions.getEps(), minimizeOptions.getIterations());
254 PDBFilter pdbFilter = (PDBFilter) systemFilter;
255 pdbFilter.writeFile(saveFile, true, false, false);
256 } else {
257 minimize.minimizeTitration(minimizeOptions.getNBFGS(), minimizeOptions.getEps(), minimizeOptions.getIterations());
258 writeFilter.writeFile(saveFile, true);
259 }
260 }
261 if (systemFilter instanceof PDBFilter) {
262 FileUtils.append(saveFile, "END\n");
263 }
264 }
265
266 return this;
267 }
268 }