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.DynamicsOptions;
42 import ffx.algorithms.cli.LambdaParticleOptions;
43 import ffx.algorithms.cli.MultiDynamicsOptions;
44 import ffx.algorithms.cli.OSTOptions;
45 import ffx.algorithms.cli.ThermodynamicsOptions;
46 import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering;
47 import ffx.numerics.Potential;
48 import ffx.potential.ForceFieldEnergy;
49 import ffx.potential.cli.AlchemicalOptions;
50 import ffx.potential.cli.AtomSelectionOptions;
51 import ffx.potential.cli.GradientOptions;
52 import ffx.potential.utils.GradientUtils;
53 import ffx.utilities.FFXBinding;
54 import org.apache.commons.io.FilenameUtils;
55 import picocli.CommandLine.Command;
56 import picocli.CommandLine.Mixin;
57 import picocli.CommandLine.Option;
58 import picocli.CommandLine.Parameters;
59
60 import java.io.File;
61 import java.util.Collections;
62 import java.util.List;
63
64 import static java.lang.String.format;
65
66
67
68
69
70
71
72
73 @Command(description = " OSTGradient script tests the Orthogonal Space Tempering Potential.", name = "test.OSTGradient")
74 public class OSTGradient extends AlgorithmsCommand {
75
76 @Mixin
77 private AtomSelectionOptions atomSelectionOptions;
78
79 @Mixin
80 private GradientOptions gradientOptions;
81
82 @Mixin
83 private AlchemicalOptions alchemicalOptions;
84
85
86
87
88 @Option(names = {"--meta", "--metaDynamics"},
89 defaultValue = "false", description = "Use a 1D metadynamics style bias.")
90 private boolean metaDynamics;
91
92
93
94
95 @Parameters(arity = "1", paramLabel = "file", description = "XYZ or PDB input file.")
96 private String filename;
97
98 private OrthogonalSpaceTempering orthogonalSpaceTempering;
99 private File saveDir = null;
100
101 public double dUdLError = 0.0;
102 public int nFailures = 0;
103
104
105
106
107 public OSTGradient() {
108 super();
109 }
110
111
112
113
114
115 public OSTGradient(FFXBinding binding) {
116 super(binding);
117 }
118
119
120
121
122
123 public OSTGradient(String[] args) {
124 super(args);
125 }
126
127
128
129
130 @Override
131 public OSTGradient run() {
132
133
134 if (!init()) {
135 return this;
136 }
137
138
139 System.setProperty("lambdaterm", "true");
140
141
142 activeAssembly = getActiveAssembly(filename);
143 if (activeAssembly == null) {
144 logger.info(helpString());
145 return this;
146 }
147
148
149 File structureFile = activeAssembly.getFile();
150 filename = structureFile.getAbsolutePath();
151
152 logger.info("\n Evaluating OST Gradient for " + filename);
153
154 String baseFilename = FilenameUtils.removeExtension(filename);
155 File histogramRestart = new File(baseFilename + ".his");
156 File lambdaRestart = null;
157
158 if (!histogramRestart.exists()) {
159 logger.warning("\n Histogram restart file does not exist: " + histogramRestart.toString());
160 } else if (!histogramRestart.canRead()) {
161 logger.warning("\n Histogram restart file can not be read." + histogramRestart.toString());
162 }
163
164
165 ForceFieldEnergy energy = activeAssembly.getPotentialEnergy();
166
167 if (saveDir == null || !saveDir.exists() || !saveDir.isDirectory() || !saveDir.canWrite()) {
168 saveDir = new File(FilenameUtils.getFullPath(filename));
169 }
170
171
172 DynamicsOptions dynamicsOptions = new DynamicsOptions();
173 LambdaParticleOptions lambdaParticleOptions = new LambdaParticleOptions();
174 MultiDynamicsOptions multiDynamicsOptions = new MultiDynamicsOptions();
175 ThermodynamicsOptions thermodynamicsOptions = new ThermodynamicsOptions();
176 OSTOptions ostOptions = new OSTOptions();
177
178 ostOptions.setMetaDynamics(metaDynamics);
179
180
181 orthogonalSpaceTempering = ostOptions.constructOST(energy, lambdaRestart, histogramRestart, activeAssembly, null,
182 dynamicsOptions, thermodynamicsOptions, lambdaParticleOptions,
183 null, !multiDynamicsOptions.isSynchronous());
184
185
186 double lambda = alchemicalOptions.getInitialLambda(false);
187 logger.info("\n Lambda value: " + lambda);
188
189
190 alchemicalOptions.setFirstSystemAlchemistry(activeAssembly);
191
192 alchemicalOptions.setFirstSystemUnchargedAtoms(activeAssembly);
193
194 atomSelectionOptions.setActiveAtoms(activeAssembly);
195
196
197
198
199
200
201 orthogonalSpaceTempering.setPropagateLambda(false);
202 orthogonalSpaceTempering.setLambda(lambda);
203 int n = orthogonalSpaceTempering.getNumberOfVariables();
204
205 if (n % 3 != 0) {
206 throw new IllegalStateException("Number of variables must be divisible by 3, got: " + n);
207 }
208 n /= 3;
209
210
211 double[] x = new double[3 * n];
212 orthogonalSpaceTempering.getCoordinates(x);
213
214 double[] g = new double[3 * n];
215 double e = orthogonalSpaceTempering.energyAndGradient(x, g);
216 double dEdLambda = orthogonalSpaceTempering.getTotaldEdLambda();
217
218
219 double step = gradientOptions.getDx();
220 orthogonalSpaceTempering.setLambda(lambda + step);
221 double lp = orthogonalSpaceTempering.energy(x);
222 orthogonalSpaceTempering.setLambda(lambda - step);
223 double lm = orthogonalSpaceTempering.energy(x);
224 double dedl = (lp - lm) / (2.0 * step);
225
226
227 logger.info(format(" Analytic dE/dL: %15.8f", dEdLambda));
228 logger.info(format(" Numeric dE/dL: %15.8f\n", dedl));
229 dUdLError = Math.abs(dEdLambda - dedl);
230
231
232 orthogonalSpaceTempering.setLambda(lambda);
233 GradientUtils gradientUtils = new GradientUtils(orthogonalSpaceTempering);
234 nFailures = gradientUtils.testGradient(gradientOptions);
235
236 return this;
237 }
238
239
240
241
242 @Override
243 public List<Potential> getPotentials() {
244 List<Potential> potentials;
245 if (orthogonalSpaceTempering == null) {
246 potentials = Collections.emptyList();
247 } else {
248 potentials = Collections.singletonList((Potential) orthogonalSpaceTempering);
249 }
250 return potentials;
251 }
252
253 }