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 Barostat barostat = barostatOptions.createBarostat(activeAssembly, crystalPotential);
169 potential = barostat;
170 }
171
172 Comm world = Comm.world();
173 int size = world.size();
174
175 if (!repEx.getRepEx() || size < 2) {
176 logger.info("\n Running molecular dynamics on " + filename);
177
178 File dyn = new File(FilenameUtils.removeExtension(filename) + ".dyn");
179 if (!dyn.exists()) {
180 dyn = null;
181 }
182
183 molDyn = dynamicsOptions.getDynamics(writeOut, potential, activeAssembly, algorithmListener);
184
185 molDyn.dynamic(dynamicsOptions.getSteps(), dynamicsOptions.getDt(),
186 dynamicsOptions.getReport(), dynamicsOptions.getWrite(), dynamicsOptions.getTemperature(), true, dyn);
187
188 } else {
189 logger.info("\n Running replica exchange molecular dynamics on " + filename);
190 int rank = (size > 1) ? world.rank() : 0;
191 logger.info("Rank:" + Integer.toString(rank));
192
193 File structureFile = new File(filename);
194 String baseFilename = FilenameUtils.removeExtension(structureFile.getName());
195 File rankDirectory = new File(structureFile.getParent() + File.separator + Integer.toString(rank));
196 if (!rankDirectory.exists()) {
197 rankDirectory.mkdir();
198 }
199
200
201 String withRankName = rankDirectory.getPath() + File.separator + baseFilename;
202 File dyn = new File(withRankName + ".dyn");
203 if (!dyn.exists()) {
204 dyn = null;
205 }
206
207 final String newMolAssemblyFile = rankDirectory.getPath() + File.separator + structureFile.getName();
208 logger.info("Set activeAssembly filename: " + newMolAssemblyFile);
209 activeAssembly.setFile(new File(newMolAssemblyFile));
210 molDyn = dynamicsOptions.getDynamics(writeOut, potential, activeAssembly, algorithmListener);
211 ReplicaExchange replicaExchange = new ReplicaExchange(molDyn, algorithmListener,
212 dynamicsOptions.getTemperature(), repEx.getExponent(), repEx.getMonteCarlo());
213
214 long totalSteps = dynamicsOptions.getSteps();
215 int nSteps = repEx.getReplicaSteps();
216 int cycles = (int) (totalSteps / nSteps);
217 if (cycles <= 0) {
218 cycles = 1;
219 }
220
221 replicaExchange.sample(cycles, nSteps, dynamicsOptions.getDt(), dynamicsOptions.getReport(), dynamicsOptions.getWrite());
222 }
223
224 return this;
225 }
226
227
228
229
230 @Override
231 public List<Potential> getPotentials() {
232 List<Potential> potentials;
233 if (potential == null) {
234 potentials = Collections.emptyList();
235 } else {
236 potentials = Collections.singletonList(potential);
237 }
238 return potentials;
239 }
240
241 }