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.RepExOptions;
45 import ffx.algorithms.dynamics.Barostat;
46 import ffx.algorithms.dynamics.MolecularDynamics;
47 import ffx.algorithms.dynamics.ReplicaExchange;
48 import ffx.crystal.CrystalPotential;
49 import ffx.numerics.Potential;
50 import ffx.potential.cli.AtomSelectionOptions;
51 import ffx.potential.cli.WriteoutOptions;
52 import ffx.utilities.FFXBinding;
53 import org.apache.commons.io.FilenameUtils;
54 import picocli.CommandLine.Command;
55 import picocli.CommandLine.Mixin;
56 import picocli.CommandLine.Parameters;
57
58 import java.io.File;
59 import java.util.Collections;
60 import java.util.List;
61
62
63
64
65
66
67
68
69 @Command(description = " Run dynamics on a system.", name = "Dynamics")
70 public class Dynamics extends AlgorithmsCommand {
71
72 @Mixin
73 private AtomSelectionOptions atomSelectionOptions;
74
75 @Mixin
76 private DynamicsOptions dynamicsOptions;
77
78 @Mixin
79 private BarostatOptions barostatOptions;
80
81 @Mixin
82 private WriteoutOptions writeOut;
83
84 @Mixin
85 private RepExOptions repEx;
86
87
88
89
90 @Parameters(arity = "1", paramLabel = "file",
91 description = "XYZ or PDB input file.")
92 private String filename;
93
94
95
96
97
98 public Potential potential = null;
99 public MolecularDynamics molDyn = null;
100
101 public MolecularDynamics getMolecularDynamics() {
102 return molDyn;
103 }
104
105 public Potential getPotentialObject() {
106 return potential;
107 }
108
109
110
111
112 public Dynamics() {
113 super();
114 }
115
116
117
118
119
120
121 public Dynamics(FFXBinding binding) {
122 super(binding);
123 }
124
125
126
127
128
129
130 public Dynamics(String[] args) {
131 super(args);
132 }
133
134
135
136
137 @Override
138 public Dynamics run() {
139
140
141 if (!init()) {
142 return this;
143 }
144
145
146 dynamicsOptions.init();
147
148
149 activeAssembly = getActiveAssembly(filename);
150 if (activeAssembly == null) {
151 logger.info(helpString());
152 return this;
153 }
154
155
156 filename = activeAssembly.getFile().getAbsolutePath();
157
158
159 atomSelectionOptions.setActiveAtoms(activeAssembly);
160
161 potential = activeAssembly.getPotentialEnergy();
162 double[] x = new double[potential.getNumberOfVariables()];
163 potential.getCoordinates(x);
164 potential.energy(x, true);
165
166 if (barostatOptions.getPressure() > 0) {
167 CrystalPotential crystalPotential = (CrystalPotential) potential;
168 potential = barostatOptions.createBarostat(activeAssembly, crystalPotential);
169 }
170
171 Comm world = Comm.world();
172 int size = world.size();
173
174 if (!repEx.getRepEx() || size < 2) {
175 logger.info("\n Running molecular dynamics on " + filename);
176
177 File dyn = new File(FilenameUtils.removeExtension(filename) + ".dyn");
178 if (!dyn.exists()) {
179 dyn = null;
180 }
181
182 molDyn = dynamicsOptions.getDynamics(writeOut, potential, activeAssembly, algorithmListener);
183
184 molDyn.dynamic(dynamicsOptions.getSteps(), dynamicsOptions.getDt(),
185 dynamicsOptions.getReport(), dynamicsOptions.getWrite(), dynamicsOptions.getTemperature(), true, dyn);
186
187 } else {
188 logger.info("\n Running replica exchange molecular dynamics on " + filename);
189 int rank = (size > 1) ? world.rank() : 0;
190 logger.info("Rank:" + rank);
191
192 File structureFile = new File(filename);
193 String baseFilename = FilenameUtils.removeExtension(structureFile.getName());
194 File rankDirectory = new File(structureFile.getParent() + File.separator + rank);
195 if (!rankDirectory.exists()) {
196 rankDirectory.mkdir();
197 }
198
199
200 String withRankName = rankDirectory.getPath() + File.separator + baseFilename;
201 File dyn = new File(withRankName + ".dyn");
202 if (!dyn.exists()) {
203 dyn = null;
204 }
205
206 final String newMolAssemblyFile = rankDirectory.getPath() + File.separator + structureFile.getName();
207 logger.info("Set activeAssembly filename: " + newMolAssemblyFile);
208 activeAssembly.setFile(new File(newMolAssemblyFile));
209 molDyn = dynamicsOptions.getDynamics(writeOut, potential, activeAssembly, algorithmListener);
210 ReplicaExchange replicaExchange = new ReplicaExchange(molDyn, algorithmListener,
211 dynamicsOptions.getTemperature(), repEx.getExponent(), repEx.getMonteCarlo());
212
213 long totalSteps = dynamicsOptions.getSteps();
214 int nSteps = repEx.getReplicaSteps();
215 int cycles = (int) (totalSteps / nSteps);
216 if (cycles <= 0) {
217 cycles = 1;
218 }
219
220 replicaExchange.sample(cycles, nSteps, dynamicsOptions.getDt(), dynamicsOptions.getReport(), dynamicsOptions.getWrite());
221 }
222
223 return this;
224 }
225
226
227
228
229 @Override
230 public List<Potential> getPotentials() {
231 List<Potential> potentials;
232 if (potential == null) {
233 potentials = Collections.emptyList();
234 } else {
235 potentials = Collections.singletonList(potential);
236 }
237 return potentials;
238 }
239
240 }