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