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.potential.cli;
39
40 import static java.lang.Double.parseDouble;
41 import static java.lang.String.format;
42
43 import ffx.numerics.Potential;
44 import ffx.numerics.switching.MultiplicativeSwitch;
45 import ffx.numerics.switching.PowerSwitch;
46 import ffx.numerics.switching.SquaredTrigSwitch;
47 import ffx.numerics.switching.UnivariateSwitchingFunction;
48 import ffx.potential.DualTopologyEnergy;
49 import ffx.potential.MolecularAssembly;
50 import ffx.potential.QuadTopologyEnergy;
51 import ffx.potential.bonded.Atom;
52 import java.util.Arrays;
53 import java.util.Collections;
54 import java.util.HashSet;
55 import java.util.List;
56 import java.util.Set;
57 import java.util.logging.Level;
58 import java.util.logging.Logger;
59 import java.util.regex.Matcher;
60 import java.util.stream.Collectors;
61 import picocli.CommandLine.ArgGroup;
62 import picocli.CommandLine.Option;
63
64
65
66
67
68
69
70
71 public class TopologyOptions {
72
73
74 public static final Logger logger = Logger.getLogger(TopologyOptions.class.getName());
75
76
77
78
79 @ArgGroup(heading = "%n Alchemical Options for Dual and Quad Topologies%n", validate = false)
80 private final TopologyOptionGroup group = new TopologyOptionGroup();
81
82
83
84
85
86
87 public String getAlchemicalAtoms2() {
88 return group.alchemicalAtoms2;
89 }
90
91
92
93
94
95
96 public String getUnchargedAtoms2() {
97 return group.unchargedAtoms2;
98 }
99
100
101
102
103
104
105 public int getNPar() {
106 return group.nPar;
107 }
108
109
110
111
112
113
114
115 public String getUnsharedA() {
116 return group.unsharedA;
117 }
118
119
120
121
122
123
124
125 public String getUnsharedB() {
126 return group.unsharedB;
127 }
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148 public String getLambdaFunction() {
149 return group.lambdaFunction;
150 }
151
152
153
154
155
156
157
158
159 public int getNumParallel(int nThreads, int nTopologies) {
160 int numParallel = getNPar();
161 if (nThreads % numParallel != 0) {
162 logger.warning(format(
163 " Number of threads available %d not evenly divisible by np %d; reverting to sequential.",
164 nThreads, numParallel));
165 numParallel = 1;
166 } else if (nTopologies % numParallel != 0) {
167 logger.warning(
168 format(" Number of topologies %d not evenly divisible by np %d; reverting to sequential.",
169 nTopologies, numParallel));
170 numParallel = 1;
171 }
172 return numParallel;
173 }
174
175
176
177
178
179
180
181
182
183
184
185
186 public Potential assemblePotential(MolecularAssembly[] assemblies, int threadsAvail,
187 StringBuilder sb) {
188 int nargs = assemblies.length;
189 int numPar = getNumParallel(threadsAvail, nargs);
190 UnivariateSwitchingFunction sf = nargs > 1 ? getSwitchingFunction() : null;
191 List<Integer> uniqueA;
192 List<Integer> uniqueB;
193 if (assemblies.length >= 4) {
194 uniqueA = getUniqueAtomsA(assemblies[0]);
195 uniqueB = getUniqueAtomsB(assemblies[2]);
196 } else {
197 uniqueA = Collections.emptyList();
198 uniqueB = Collections.emptyList();
199 }
200 return getTopology(assemblies, sf, uniqueA, uniqueB, numPar, sb);
201 }
202
203
204
205
206
207
208 public UnivariateSwitchingFunction getSwitchingFunction() {
209 UnivariateSwitchingFunction sf;
210 String lambdaFunction = getLambdaFunction();
211 if (!lambdaFunction.equalsIgnoreCase("1.0")) {
212 String lf = lambdaFunction.toUpperCase();
213 switch (lf) {
214 case "TRIG" -> sf = new SquaredTrigSwitch(false);
215 case "MULT" -> sf = new MultiplicativeSwitch(1.0, 0.0);
216 default -> {
217 try {
218 double beta = parseDouble(lf);
219 sf = new PowerSwitch(1.0, beta);
220 } catch (NumberFormatException ex) {
221 logger.warning(format(
222 "Argument to option -sf %s could not be properly parsed; using default linear switch",
223 lambdaFunction));
224 sf = new PowerSwitch(1.0, 1.0);
225 }
226 }
227 }
228 } else {
229 sf = new PowerSwitch(1.0, 1.0);
230 }
231 return sf;
232 }
233
234
235
236
237
238
239
240 public List<Integer> getUniqueAtomsA(MolecularAssembly topology) {
241 return getUniqueAtoms(topology, "A", getUnsharedA());
242 }
243
244
245
246
247
248
249
250
251
252
253
254
255 public Potential getTopology(MolecularAssembly[] topologies, UnivariateSwitchingFunction sf,
256 List<Integer> uniqueA, List<Integer> uniqueB, int numParallel, StringBuilder sb) {
257 Potential potential = null;
258 switch (topologies.length) {
259 case 1 -> {
260 sb.append("Single Topology ");
261 potential = topologies[0].getPotentialEnergy();
262 }
263 case 2 -> {
264 sb.append("Dual Topology ");
265 DualTopologyEnergy dte = new DualTopologyEnergy(topologies[0], topologies[1], sf);
266 if (numParallel == 2) {
267 dte.setParallel(true);
268 }
269 potential = dte;
270 }
271 case 4 -> {
272 sb.append("Quad Topology ");
273 DualTopologyEnergy dta = new DualTopologyEnergy(topologies[0], topologies[1], sf);
274 DualTopologyEnergy dtb = new DualTopologyEnergy(topologies[3], topologies[2], sf);
275 QuadTopologyEnergy qte = new QuadTopologyEnergy(dta, dtb, uniqueA, uniqueB);
276 if (numParallel >= 2) {
277 qte.setParallel(true);
278 if (numParallel == 4) {
279 dta.setParallel(true);
280 dtb.setParallel(true);
281 }
282 }
283 potential = qte;
284 }
285 default -> logger.severe(" Must have 2, 4, or 8 topologies!");
286 }
287 sb.append(Arrays.stream(topologies).map(MolecularAssembly::toString)
288 .collect(Collectors.joining(", ", " [", "] ")));
289 return potential;
290 }
291
292
293
294
295
296
297
298
299
300 public static List<Integer> getUniqueAtoms(MolecularAssembly assembly, String label,
301 String unshared) {
302 if (!unshared.isEmpty()) {
303 logger.info(" Finding unique atoms for dual topology " + label);
304 Set<Integer> indices = new HashSet<>();
305 String[] toks = unshared.split("\\.");
306 Atom[] atoms1 = assembly.getAtomArray();
307 for (String range : toks) {
308 Matcher m = AlchemicalOptions.rangeRegEx.matcher(range);
309 if (m.find()) {
310 int rangeStart = Integer.parseInt(m.group(1));
311 int rangeEnd = (m.group(2) != null) ? Integer.parseInt(m.group(2)) : rangeStart;
312 if (rangeStart > rangeEnd) {
313 logger.severe(format(" Range %s was invalid start was greater than end", range));
314 }
315 logger.info(
316 format(" Range %s for %s, start %d end %d", range, label, rangeStart, rangeEnd));
317 if(logger.isLoggable(Level.FINE)) {
318 logger.fine(format(" First atom in range: %s", atoms1[rangeStart - 1]));
319 if (rangeEnd > rangeStart) {
320 logger.fine(format(" Last atom in range: %s", atoms1[rangeEnd - 1]));
321 }
322 }
323 for (int i = rangeStart; i <= rangeEnd; i++) {
324 indices.add(i - 1);
325 }
326 }
327 }
328 int counter = 0;
329 Set<Integer> adjustedIndices = new HashSet<>();
330 int atomLength = atoms1.length;
331 for (int i = 0; i < atomLength; i++) {
332 Atom ai = atoms1[i];
333 if (indices.contains(i)) {
334 if (ai.applyLambda()) {
335 logger.warning(format(
336 " Ranges defined in %s should not overlap with ligand atoms they are assumed to not be shared.",
337 label));
338 } else {
339 logger.fine(format(" Unshared %s: %d variables %d-%d", label, i, counter, counter + 2));
340 for (int j = 0; j < 3; j++) {
341 adjustedIndices.add(counter + j);
342 }
343 }
344 }
345 if (!ai.applyLambda()) {
346 counter += 3;
347 }
348 }
349 return adjustedIndices.stream().sorted().collect(Collectors.toList());
350 } else {
351 return Collections.emptyList();
352 }
353 }
354
355
356
357
358
359
360
361 public List<Integer> getUniqueAtomsB(MolecularAssembly topology) {
362 return getUniqueAtoms(topology, "B", getUnsharedB());
363 }
364
365
366
367
368
369
370 public boolean hasSoftcore() {
371 String alchemicalAtoms2 = getAlchemicalAtoms2();
372 return (alchemicalAtoms2 != null && !alchemicalAtoms2.equalsIgnoreCase("NONE")
373 && !alchemicalAtoms2.equalsIgnoreCase(""));
374 }
375
376
377
378
379
380
381 public void setSecondSystemAlchemistry(MolecularAssembly topology) {
382 String alchemicalAtoms2 = getAlchemicalAtoms2();
383 AlchemicalOptions.setAlchemicalAtoms(topology, alchemicalAtoms2);
384 }
385
386
387
388
389
390
391 public void setSecondSystemUnchargedAtoms(MolecularAssembly topology) {
392 String unchargedAtoms2 = getUnchargedAtoms2();
393 AlchemicalOptions.setUnchargedAtoms(topology, unchargedAtoms2);
394 }
395
396
397
398
399 private static class TopologyOptionGroup {
400
401
402 @Option(names = {"--ac2", "--alchemicalAtoms2"}, paramLabel = "<selection>", defaultValue = "",
403 description = "Specify alchemical atoms for the 2nd topology [ALL, NONE, Range(s): 1-3,6-N].")
404 String alchemicalAtoms2 = "";
405
406
407
408
409
410 @Option(names = {"--uc2", "--unchargedAtoms2"}, paramLabel = "<selection>", defaultValue = "",
411 description = "Specify atoms without electrostatics for the 2nd topology [ALL, NONE, Range(s): 1-3,6-N].")
412 String unchargedAtoms2 = "";
413
414
415
416
417
418 @Option(names = {"--np", "--nParallel"}, paramLabel = "1",
419 description = "Number of topologies to evaluate in parallel")
420 int nPar = 1;
421
422
423
424
425
426 @Option(names = {"--uaA", "--unsharedA"}, paramLabel = "-1",
427 description = "Unshared atoms in the A dual topology (e.g. 1-24.32-65).")
428 String unsharedA = null;
429
430
431
432
433
434 @Option(names = {"--uaB", "--unsharedB"}, paramLabel = "-1",
435 description = "Unshared atoms in the B dual topology (e.g. 1-24.32-65).")
436 String unsharedB = null;
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463 @Option(names = {"--sf", "--switchingFunction"}, paramLabel = "1.0",
464 description = "Switching function to use for dual topology: options are TRIG, MULT, or a number (original behavior with specified lambda exponent)")
465 String lambdaFunction = "1.0";
466 }
467 }