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.MolecularAssembly;
43 import ffx.potential.cli.PotentialCommand;
44 import ffx.potential.cli.WriteoutOptions;
45 import ffx.potential.parsers.SystemFilter;
46 import ffx.utilities.FFXBinding;
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
68
69 @Option(names = {"-s", "--start"}, paramLabel = "1", defaultValue = "1",
70 description = "First frame to evaluate (1-indexed).")
71 private int start = 1;
72
73
74
75
76 @Option(names = {"-f", "--final"}, paramLabel = "all frames", defaultValue = "0",
77 description = "Last frame to evaluate (1-indexed); values less than 1 evaluate to end of trajectory.")
78 private int finish = 0;
79
80
81
82
83 @Option(names = {"--st", "--stride"}, paramLabel = "1", defaultValue = "1",
84 description = "Stride: evaluate density every N frames. Must be positive.")
85 private int stride = 1;
86
87
88
89
90 @Option(names = {"-p", "--printout"}, defaultValue = "false",
91 description = "Print out a file with density adjusted to match mean calculated density.")
92 private boolean doPrint = false;
93
94
95
96
97 @Parameters(arity = "1", paramLabel = "file",
98 description = "An atomic coordinate file in PDB or XYZ format.")
99 private String filename = null;
100
101 public Density() {
102 super();
103 }
104
105 public Density(FFXBinding binding) {
106 super(binding);
107 }
108
109 public Density(String[] args) {
110 super(args);
111 }
112
113 @Override
114 public Density run() {
115
116 if (!init()) {
117 return this;
118 }
119
120
121 activeAssembly = getActiveAssembly(filename);
122 if (activeAssembly == null) {
123 logger.info(helpString());
124 return this;
125 }
126
127
128 filename = activeAssembly.getFile().getAbsolutePath();
129
130 SystemFilter openFilter = potentialFunctions.getFilter();
131 double totMass = activeAssembly.getMass();
132 Crystal crystal = activeAssembly.getCrystal().getUnitCell();
133
134 if (crystal.aperiodic()) {
135 logger.info(format(" System %s is aperiodic: total mass %16.7g g/mol", filename, totMass));
136 } else {
137 double volume = crystal.volume / crystal.getNumSymOps();
138 double density = crystal.getDensity(totMass);
139 int nFrames = openFilter.countNumModels();
140 double[] densities = new double[nFrames];
141 int lastFrame = (finish < 1) ? nFrames : finish;
142
143 densities[0] = density;
144 logger.info(format("\n Evaluating density for system %s with mass %16.7g (g/mol).", filename, totMass));
145 logger.info(format(" Density at frame %9d is %16.7g (g/mL) from a volume of %16.7g (A^3)", 1, density, volume));
146
147 int ctr = 1;
148 while (openFilter.readNext(false, false)) {
149 volume = crystal.volume / crystal.getNumSymOps();
150 density = crystal.getDensity(totMass);
151 logger.info(format(" Density at frame %9d is %16.7g (g/mL) from a volume of %16.7g (A^3)", ++ctr, density, volume));
152 int i = ctr - 1;
153 densities[i] = density;
154 }
155
156 SummaryStatistics densStats = new SummaryStatistics(densities, start - 1, lastFrame, stride);
157 logger.info(" Summary statistics for density:");
158 logger.info(densStats.toString());
159
160 if (doPrint) {
161 activeAssembly.setFractionalMode(MolecularAssembly.FractionalMode.MOLECULE);
162 activeAssembly.computeFractionalCoordinates();
163 crystal.setDensity(densStats.mean, totMass);
164 activeAssembly.moveToFractionalCoordinates();
165 saveByExtension(activeAssembly, filename, writeoutOptions.fileType);
166
167 filename = activeAssembly.getFile().getAbsolutePath();
168 logger.info("\n Saved density adjusted coordinates: " + filename);
169 logger.info("\n Molecule based fractional coordinates were maintained.");
170 }
171 }
172
173 return this;
174 }
175 }