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 edu.rit.pj.Comm;
41 import ffx.algorithms.cli.AlgorithmsCommand;
42 import ffx.algorithms.cli.BarostatOptions;
43 import ffx.algorithms.cli.DynamicsOptions;
44 import ffx.algorithms.cli.LambdaParticleOptions;
45 import ffx.algorithms.cli.MultiDynamicsOptions;
46 import ffx.algorithms.cli.OSTOptions;
47 import ffx.algorithms.cli.RandomUnitCellOptions;
48 import ffx.algorithms.cli.ThermodynamicsOptions;
49 import ffx.algorithms.cli.ThermodynamicsOptions.ThermodynamicsAlgorithm;
50 import ffx.algorithms.thermodynamics.MonteCarloOST;
51 import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering;
52 import ffx.crystal.CrystalPotential;
53 import ffx.numerics.Potential;
54 import ffx.potential.MolecularAssembly;
55 import ffx.potential.bonded.LambdaInterface;
56 import ffx.potential.cli.AlchemicalOptions;
57 import ffx.potential.cli.TopologyOptions;
58 import ffx.potential.cli.WriteoutOptions;
59 import ffx.utilities.FFXBinding;
60 import org.apache.commons.configuration2.Configuration;
61 import org.apache.commons.io.FilenameUtils;
62 import picocli.CommandLine.Command;
63 import picocli.CommandLine.Mixin;
64 import picocli.CommandLine.Option;
65 import picocli.CommandLine.Parameters;
66
67 import java.io.File;
68 import java.util.ArrayList;
69 import java.util.Collections;
70 import java.util.List;
71
72 import static java.lang.String.format;
73
74
75
76
77
78
79
80
81
82 @Command(description = " Use the Transition-Tempered Orthogonal Space Random Walk algorithm to estimate a free energy.", name = "Thermodynamics")
83 public class Thermodynamics extends AlgorithmsCommand {
84
85 @Mixin
86 public DynamicsOptions dynamicsOptions;
87
88 @Mixin
89 public BarostatOptions barostatOptions;
90
91 @Mixin
92 public RandomUnitCellOptions randomSymopOptions;
93
94 @Mixin
95 public AlchemicalOptions alchemicalOptions;
96
97 @Mixin
98 public TopologyOptions topologyOptions;
99
100 @Mixin
101 public WriteoutOptions writeoutOptions;
102
103 @Mixin
104 public ThermodynamicsOptions thermodynamicsOptions;
105
106 @Mixin
107 public OSTOptions ostOptions;
108
109 @Mixin
110 public LambdaParticleOptions lambdaParticleOptions;
111
112 @Mixin
113 public MultiDynamicsOptions multiDynamicsOptions;
114
115
116
117
118 @Option(names = {"-v", "--verbose"},
119 description = "Log additional information (primarily for MC-OST).")
120 public boolean verbose = false;
121
122
123
124
125 @Parameters(arity = "1..*", paramLabel = "files", description = "The atomic coordinate file in PDB or XYZ format.")
126 public List<String> filenames = null;
127
128 public MolecularAssembly[] topologies;
129 public CrystalPotential potential;
130 public OrthogonalSpaceTempering orthogonalSpaceTempering = null;
131 public Configuration additionalProperties;
132
133
134
135
136
137
138 public void setProperties(Configuration additionalProps) {
139 this.additionalProperties = additionalProps;
140 }
141
142
143
144
145 public Thermodynamics() {
146 super();
147 }
148
149
150
151
152
153
154 public Thermodynamics(FFXBinding binding) {
155 super(binding);
156 }
157
158
159
160
161
162
163 public Thermodynamics(String[] args) {
164 super(args);
165 }
166
167
168
169
170 @Override
171 public Thermodynamics run() {
172
173
174 if (!init()) {
175 return this;
176 }
177
178
179 int numTopologies = topologyOptions.getNumberOfTopologies(filenames);
180 int threadsPerTopology = topologyOptions.getThreadsPerTopology(numTopologies);
181 topologies = new MolecularAssembly[numTopologies];
182
183
184 alchemicalOptions.setAlchemicalProperties();
185 topologyOptions.setAlchemicalProperties(numTopologies);
186
187 Comm world = Comm.world();
188 int size = world.size();
189 int rank = (size > 1) ? world.rank() : 0;
190
191
192 List<File> structureFiles = new ArrayList<>();
193 for (String filename : filenames) {
194 File file = new File(FilenameUtils.normalize(filename));
195 structureFiles.add(file);
196 }
197
198 File firstStructure = structureFiles.get(0);
199 String filePathNoExtension = firstStructure.getAbsolutePath().replaceFirst("\\.[^.]+$", "");
200 File histogramRestart = new File(filePathNoExtension + ".his");
201
202
203 String withRankName = filePathNoExtension;
204
205 if (size > 1) {
206 List<File> rankedFiles = new ArrayList<>(numTopologies);
207 String rankDirName = FilenameUtils.getFullPath(filePathNoExtension);
208 rankDirName = format("%s%d", rankDirName, rank + multiDynamicsOptions.getFirstDir());
209 File rankDirectory = new File(rankDirName);
210 if (!rankDirectory.exists()) {
211 rankDirectory.mkdir();
212 }
213 rankDirName = rankDirName + File.separator;
214 withRankName = format("%s%s", rankDirName, FilenameUtils.getName(filePathNoExtension));
215
216 for (File structureFile : structureFiles) {
217 rankedFiles.add(new File(format("%s%s", rankDirName,
218 FilenameUtils.getName(structureFile.getName()))));
219 }
220 structureFiles = rankedFiles;
221 }
222
223 File lambdaRestart = new File(withRankName + ".lam");
224 File dyn = new File(withRankName + ".dyn");
225 if (ostOptions.getIndependentWalkers()) {
226 histogramRestart = new File(withRankName + ".his");
227 }
228
229
230 if (filenames == null || filenames.isEmpty()) {
231 activeAssembly = getActiveAssembly(null);
232 if (activeAssembly == null) {
233 logger.info(helpString());
234 return this;
235 }
236 filenames = new ArrayList<>();
237 filenames.add(activeAssembly.getFile().getName());
238 topologies[0] = alchemicalOptions.processFile(topologyOptions, activeAssembly, 0);
239 } else {
240 logger.info(format(" Initializing %d topologies...", numTopologies));
241 for (int i = 0; i < numTopologies; i++) {
242 topologies[i] = multiDynamicsOptions.openFile(algorithmFunctions, topologyOptions,
243 threadsPerTopology, filenames.get(i), i, alchemicalOptions, structureFiles.get(i), rank);
244 }
245 }
246
247 StringBuilder sb = new StringBuilder("\n Running ");
248
249 ThermodynamicsAlgorithm algorithm = thermodynamicsOptions.getAlgorithm();
250 double initLambda = alchemicalOptions.getInitialLambda(size, rank, true);
251 if (algorithm == ThermodynamicsAlgorithm.OST) {
252 sb.append("Orthogonal Space Tempering");
253 } else if (algorithm == ThermodynamicsAlgorithm.FIXED) {
254 sb.append("Fixed Lambda Sampling at Window L=").append(format("%5.3f ", initLambda));
255 } else if (algorithm == ThermodynamicsAlgorithm.NEQ) {
256 sb.append("Non-Equilibrium Sampling");
257 } else {
258 logger.severe(" Unknown Thermodynamics Algorithm " + algorithm);
259 }
260 sb.append(" for ");
261
262 potential = (CrystalPotential) topologyOptions.assemblePotential(topologies, sb);
263 logger.info(sb.toString());
264
265 LambdaInterface lambdaInterface = (LambdaInterface) potential;
266 boolean lamExists = lambdaRestart.exists();
267 double[] x = new double[potential.getNumberOfVariables()];
268 potential.getCoordinates(x);
269 lambdaInterface.setLambda(initLambda);
270 potential.energy(x, true);
271
272 if (numTopologies == 1) {
273 randomSymopOptions.randomize(topologies[0]);
274 }
275
276 multiDynamicsOptions.distribute(topologies, potential, algorithmFunctions, rank, size);
277
278 if (algorithm == ThermodynamicsAlgorithm.OST) {
279 orthogonalSpaceTempering =
280 ostOptions.constructOST(potential, lambdaRestart, histogramRestart, topologies[0],
281 additionalProperties, dynamicsOptions, thermodynamicsOptions, lambdaParticleOptions,
282 algorithmListener, !multiDynamicsOptions.isSynchronous());
283 if (!lamExists) {
284 orthogonalSpaceTempering.setLambda(initLambda);
285 }
286
287 CrystalPotential ostPotential = ostOptions.applyAllOSTOptions(orthogonalSpaceTempering, topologies[0],
288 dynamicsOptions, barostatOptions);
289 if (ostOptions.isMonteCarlo()) {
290 MonteCarloOST mcOST = ostOptions.setupMCOST(orthogonalSpaceTempering, topologies,
291 dynamicsOptions, thermodynamicsOptions, verbose, dyn, algorithmListener);
292 ostOptions.beginMCOST(mcOST, dynamicsOptions, thermodynamicsOptions);
293 } else {
294 ostOptions.beginMDOST(orthogonalSpaceTempering, topologies, ostPotential, dynamicsOptions,
295 writeoutOptions, thermodynamicsOptions, dyn, algorithmListener);
296 }
297 logger.info(" Done running OST sampling.");
298 } else if (algorithm == ThermodynamicsAlgorithm.FIXED) {
299 orthogonalSpaceTempering = null;
300 potential = barostatOptions.checkNPT(topologies[0], potential);
301 thermodynamicsOptions.runFixedAlchemy(topologies, potential, dynamicsOptions, writeoutOptions, dyn, algorithmListener);
302 logger.info(" Done running fixed lambda sampling.");
303 } else if (algorithm == ThermodynamicsAlgorithm.NEQ) {
304 orthogonalSpaceTempering = null;
305 potential = barostatOptions.checkNPT(topologies[0], potential);
306 thermodynamicsOptions.runNEQ(topologies, potential, dynamicsOptions, writeoutOptions, dyn, algorithmListener);
307 logger.info(" Done running non-equilibrium sampling.");
308 } else {
309 logger.severe(" Unknown Thermodynamics Algorithm " + algorithm);
310 }
311
312 return this;
313 }
314
315 public OrthogonalSpaceTempering getOST() {
316 return orthogonalSpaceTempering;
317 }
318
319 public CrystalPotential getPotential() {
320 return potential;
321 }
322
323 @Override
324 public List<Potential> getPotentials() {
325 List<Potential> potentials;
326 if (orthogonalSpaceTempering == null) {
327 if (potential == null) {
328 potentials = Collections.emptyList();
329 } else {
330 potentials = Collections.singletonList((Potential) potential);
331 }
332 } else {
333 potentials = Collections.singletonList((Potential) orthogonalSpaceTempering);
334 }
335 return potentials;
336 }
337
338 @Override
339 public boolean destroyPotentials() {
340 return getPotentials().stream().allMatch(potential -> potential.destroy());
341 }
342 }