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.commands;
39
40 import ffx.crystal.Crystal;
41 import ffx.numerics.math.SummaryStatistics;
42 import ffx.potential.cli.PotentialCommand;
43 import ffx.potential.cli.WriteoutOptions;
44 import ffx.potential.parsers.SystemFilter;
45 import ffx.utilities.FFXBinding;
46 import org.apache.commons.io.FilenameUtils;
47 import picocli.CommandLine.Command;
48 import picocli.CommandLine.Mixin;
49 import picocli.CommandLine.Option;
50 import picocli.CommandLine.Parameters;
51
52 import static java.lang.String.format;
53
54
55
56
57
58
59
60 @Command(name = "Density", description = " Calculate the mean system density for single topology systems.")
61 public class Density extends PotentialCommand {
62
63 @Mixin
64 private WriteoutOptions writeoutOptions = new WriteoutOptions();
65
66
67 @Option(names = {"-s", "--start"}, paramLabel = "1", defaultValue = "1",
68 description = "First frame to evaluate (1-indexed).")
69 private int start = 1;
70
71
72 @Option(names = {"-f", "--final"}, paramLabel = "all frames", defaultValue = "0",
73 description = "Last frame to evaluate (1-indexed); values less than 1 evaluate to end of trajectory.")
74 private int finish = 0;
75
76
77 @Option(names = {"--st", "--stride"}, paramLabel = "1", defaultValue = "1",
78 description = "Stride: evaluate density every N frames. Must be positive.")
79 private int stride = 1;
80
81
82 @Option(names = {"-p", "--printout"}, defaultValue = "false",
83 description = "Print out a file with density adjusted to match mean calculated density.")
84 private boolean doPrint = false;
85
86
87 @Parameters(arity = "1", paramLabel = "file",
88 description = "An atomic coordinate file in PDB or XYZ format.")
89 private String filename = null;
90
91 public Density() {
92 super();
93 }
94
95 public Density(FFXBinding binding) {
96 super(binding);
97 }
98
99 public Density(String[] args) {
100 super(args);
101 }
102
103 @Override
104 public Density run() {
105
106 if (!init()) {
107 return this;
108 }
109
110
111 activeAssembly = getActiveAssembly(filename);
112 if (activeAssembly == null) {
113 logger.info(helpString());
114 return this;
115 }
116
117
118 filename = activeAssembly.getFile().getAbsolutePath();
119
120 SystemFilter openFilter = potentialFunctions.getFilter();
121 double totMass = activeAssembly.getMass();
122 Crystal crystal = activeAssembly.getCrystal().getUnitCell();
123
124 if (crystal.aperiodic()) {
125 logger.info(format(" System %s is aperiodic: total mass %16.7g g/mol", filename, totMass));
126 } else {
127 double volume = crystal.volume / crystal.getNumSymOps();
128 double density = crystal.getDensity(totMass);
129 int nFrames = openFilter.countNumModels();
130 double[] densities = new double[nFrames];
131 int lastFrame = (finish < 1) ? nFrames : finish;
132
133 densities[0] = density;
134 logger.info(format(" Evaluating density for system %s with mass %16.7g (g/mol).", filename, totMass));
135 logger.info(format(" Density at frame %9d is %16.7g (g/mL) from a volume of %16.7g (A^3)", 1, density, volume));
136
137 int ctr = 1;
138 while (openFilter.readNext(false, false)) {
139 volume = crystal.volume / crystal.getNumSymOps();
140 density = crystal.getDensity(totMass);
141 logger.info(format(" Density at frame %9d is %16.7g (g/mL) from a volume of %16.7g (A^3)", ++ctr, density, volume));
142 int i = ctr - 1;
143 densities[i] = density;
144 }
145
146 SummaryStatistics densStats = new SummaryStatistics(densities, start - 1, lastFrame, stride);
147 logger.info(" Summary statistics for density:");
148 logger.info(densStats.toString());
149
150 if (doPrint) {
151 crystal.setDensity(densStats.mean, totMass);
152 String outFileName = FilenameUtils.removeExtension(filename);
153 writeoutOptions.saveFile(outFileName, potentialFunctions, activeAssembly);
154 }
155 }
156
157 return this;
158 }
159 }