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 edu.rit.pj.Comm;
41 import ffx.algorithms.cli.RepexOSTOptions;
42 import ffx.algorithms.cli.ThermodynamicsOptions;
43 import ffx.algorithms.commands.Thermodynamics;
44 import ffx.algorithms.dynamics.MolecularDynamics;
45 import ffx.algorithms.thermodynamics.HistogramData;
46 import ffx.algorithms.thermodynamics.LambdaData;
47 import ffx.algorithms.thermodynamics.MonteCarloOST;
48 import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering;
49 import ffx.algorithms.thermodynamics.RepExOST;
50 import ffx.crystal.CrystalPotential;
51 import ffx.numerics.Potential;
52 import ffx.potential.MolecularAssembly;
53 import ffx.potential.bonded.LambdaInterface;
54 import ffx.utilities.FFXBinding;
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
60 import java.io.File;
61 import java.io.IOException;
62 import java.util.ArrayList;
63 import java.util.Collections;
64 import java.util.List;
65 import java.util.stream.Collectors;
66
67
68
69
70
71
72
73
74 @Command(description = " Use Orthogonal Space Tempering with histogram replica exchange to estimate a free energy difference.", name = "test.RepexThermo")
75 public class RepexThermo extends Thermodynamics {
76
77 @Mixin
78 private RepexOSTOptions repex;
79
80 private RepExOST repExOST;
81 private CrystalPotential finalPotential;
82
83
84
85
86 public RepexThermo() {
87 super();
88 }
89
90
91
92
93
94
95 public RepexThermo(FFXBinding binding) {
96 super(binding);
97 }
98
99
100
101
102
103
104 public RepexThermo(String[] args) {
105 super(args);
106 }
107
108
109
110
111 @Override
112 public RepexThermo run() {
113
114
115 if (!init()) {
116 return this;
117 }
118
119 boolean fromActive;
120 List<String> arguments;
121 if (filenames != null && !filenames.isEmpty()) {
122 arguments = filenames;
123 fromActive = false;
124 } else {
125 logger.warning(" Untested: use of active assembly instead of provided filenames!");
126 MolecularAssembly mola = algorithmFunctions.getActiveAssembly();
127 if (mola == null) {
128 logger.info(helpString());
129 return this;
130 }
131 arguments = Collections.singletonList(mola.getFile().getName());
132 fromActive = true;
133 }
134
135 int nArgs = arguments.size();
136
137 topologies = new MolecularAssembly[nArgs];
138
139
140 int numTopologies = topologyOptions.getNumberOfTopologies(filenames);
141 int threadsPer = topologyOptions.getThreadsPerTopology(numTopologies);
142 topologies = new MolecularAssembly[numTopologies];
143
144
145
146
147 boolean lambdaTerm = (
148 nArgs == 1 || alchemicalOptions.hasSoftcore() || topologyOptions.hasSoftcore());
149
150 if (lambdaTerm) {
151 System.setProperty("lambdaterm", "true");
152 }
153
154
155
156 if (nArgs >= 2) {
157
158
159 System.setProperty("ligand-vapor-elec", "false");
160 }
161
162 List<MolecularAssembly> topologyList = new ArrayList<>(nArgs);
163
164 Comm world = Comm.world();
165 int size = world.size();
166 if (size < 2) {
167 logger.severe(" RepexThermo requires multiple processes, found only one!");
168 }
169 int rank = (size > 1) ? world.rank() : 0;
170
171 double initLambda = alchemicalOptions.getInitialLambda(size, rank);
172
173
174 List<File> structureFiles = arguments.stream()
175 .map(fn -> new File(new File(FilenameUtils.normalize(fn)).getAbsolutePath()))
176 .collect(Collectors.toList());
177
178 File firstStructure = structureFiles.get(0);
179 String filePathNoExtension = firstStructure.getAbsolutePath().replaceFirst("\\.[^.]+$", "");
180
181
182
183 String filepath = FilenameUtils.getFullPath(filePathNoExtension);
184 String fileBase = FilenameUtils.getBaseName(FilenameUtils.getName(filePathNoExtension));
185 String rankDirName = String.format("%s%d", filepath, rank);
186 File rankDirectory = new File(rankDirName);
187 if (!rankDirectory.exists()) {
188 rankDirectory.mkdir();
189 }
190 rankDirName = rankDirName + File.separator;
191 String withRankName = rankDirName + fileBase;
192 File lambdaRestart = new File(withRankName + ".lam");
193
194 boolean lamExists = lambdaRestart.exists();
195
196
197 logger.info(String.format(" Initializing %d topologies...", nArgs));
198 if (fromActive) {
199 topologyList.add(alchemicalOptions.processFile(topologyOptions, activeAssembly, 0));
200 } else {
201 for (int i = 0; i < nArgs; i++) {
202 topologyList.add(multiDynamicsOptions.openFile(algorithmFunctions, topologyOptions,
203 threadsPer, arguments.get(i), i, alchemicalOptions, structureFiles.get(i), rank));
204 }
205 }
206
207 MolecularAssembly[] topologies =
208 topologyList.toArray(new MolecularAssembly[0]);
209
210 StringBuilder sb = new StringBuilder("\n Running ");
211 switch (thermodynamicsOptions.getAlgorithm()) {
212 case OST:
213 sb.append("Orthogonal Space Tempering");
214 break;
215 default:
216 throw new IllegalArgumentException(
217 " RepexThermo currently does not support fixed-lambda alchemy!");
218 }
219 sb.append(" with histogram replica exchange for ");
220
221 potential = (CrystalPotential) topologyOptions.assemblePotential(topologies, sb);
222
223 LambdaInterface linter = (LambdaInterface) potential;
224 logger.info(sb.toString());
225
226 double[] x = new double[potential.getNumberOfVariables()];
227 potential.getCoordinates(x);
228 linter.setLambda(initLambda);
229 potential.energy(x, true);
230
231 if (nArgs == 1) {
232 randomSymopOptions.randomize(topologies[0]);
233 }
234
235 multiDynamicsOptions.distribute(topologies, potential, algorithmFunctions, rank, size);
236
237 boolean isMC = ostOptions.isMonteCarlo();
238 boolean twoStep = ostOptions.isTwoStep();
239 MonteCarloOST mcOST = null;
240 MolecularDynamics md;
241
242 if (thermodynamicsOptions.getAlgorithm() == ThermodynamicsOptions.ThermodynamicsAlgorithm.OST) {
243 File firstHisto = new File(filepath + "0" + File.separator + fileBase + ".his");
244
245 orthogonalSpaceTempering = ostOptions.constructOST(potential, lambdaRestart, firstHisto, topologies[0],
246 additionalProperties, dynamicsOptions, thermodynamicsOptions, lambdaParticleOptions,
247 algorithmListener,
248 false);
249 finalPotential = ostOptions.applyAllOSTOptions(orthogonalSpaceTempering, topologies[0],
250 dynamicsOptions, barostatOptions);
251
252 if (isMC) {
253 mcOST = ostOptions.setupMCOST(orthogonalSpaceTempering, topologies, dynamicsOptions, thermodynamicsOptions,
254 verbose, null, algorithmListener);
255 md = mcOST.getMD();
256 } else {
257 md = ostOptions.assembleMolecularDynamics(topologies, finalPotential, dynamicsOptions, algorithmListener);
258 }
259 if (!lamExists) {
260 if (finalPotential instanceof LambdaInterface) {
261 ((LambdaInterface) finalPotential).setLambda(initLambda);
262 } else {
263 orthogonalSpaceTempering.setLambda(initLambda);
264 }
265 }
266
267 CompositeConfiguration allProperties = new CompositeConfiguration(
268 topologies[0].getProperties());
269 if (additionalProperties != null) {
270 allProperties.addConfiguration(additionalProperties);
271 }
272
273 for (int i = 1; i < size; i++) {
274 File histogramFile = new File(filepath + i + File.separator + fileBase + ".his");
275 File lambdaFile = new File(filepath + i + File.separator + fileBase + ".lam");
276 HistogramData histogramData = HistogramData.readHistogram(histogramFile);
277 LambdaData lambdaData = LambdaData.readLambdaData(lambdaFile);
278 orthogonalSpaceTempering.addHistogram(histogramData, lambdaData);
279 }
280
281 if (isMC) {
282 try {
283 repExOST = RepExOST.repexMC(orthogonalSpaceTempering, mcOST, dynamicsOptions, ostOptions,
284 topologies[0].getProperties(), writeoutOptions.getFileType(), twoStep,
285 repex.getRepexFrequency());
286 } catch (IOException e) {
287 throw new RuntimeException(e);
288 }
289 } else {
290 try {
291 repExOST = RepExOST.repexMD(orthogonalSpaceTempering, md, dynamicsOptions, ostOptions,
292 topologies[0].getProperties(), writeoutOptions.getFileType(), repex.getRepexFrequency());
293 } catch (IOException e) {
294 throw new RuntimeException(e);
295 }
296 }
297
298 long eSteps = thermodynamicsOptions.getEquilSteps();
299 if (eSteps > 0) {
300 try {
301 repExOST.mainLoop(eSteps, true);
302 } catch (IOException e) {
303 throw new RuntimeException(e);
304 }
305 }
306 try {
307 repExOST.mainLoop(dynamicsOptions.getNumSteps(), false);
308 } catch (IOException e) {
309 throw new RuntimeException(e);
310 }
311 } else {
312 logger.severe(" RepexThermo currently does not support fixed-lambda alchemy!");
313 }
314
315 logger.info(" " + thermodynamicsOptions.getAlgorithm() + " with Histogram Replica Exchange Done.");
316
317 return this;
318 }
319
320 @Override
321 public OrthogonalSpaceTempering getOST() {
322 return repExOST == null ? null : repExOST.getOST();
323 }
324
325 @Override
326 public CrystalPotential getPotential() {
327 return (repExOST == null) ? potential : repExOST.getOST();
328 }
329
330 @Override
331 public List<Potential> getPotentials() {
332 if (repExOST == null) {
333 if (potential == null) {
334 return Collections.emptyList();
335 } else {
336 return Collections.singletonList(potential);
337 }
338 }
339 return Collections.singletonList((Potential) repExOST.getOST());
340 }
341 }