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 edu.rit.pj.ParallelTeam;
41 import ffx.potential.bonded.Atom;
42 import ffx.potential.cli.PotentialCommand;
43 import ffx.potential.nonbonded.implicit.ConnollyRegion;
44 import ffx.potential.nonbonded.implicit.GaussVol;
45 import ffx.potential.parameters.ForceField;
46 import ffx.utilities.FFXBinding;
47 import picocli.CommandLine.Command;
48 import picocli.CommandLine.Option;
49 import picocli.CommandLine.Parameters;
50
51 import static ffx.utilities.Constants.NS2SEC;
52 import static java.lang.String.format;
53 import static org.apache.commons.math3.util.FastMath.pow;
54
55
56
57
58
59
60
61 @Command(name = "Volume", description = " Calculate the surface area and volume using the GaussVol (default) or Connolly algorithm.")
62 public class Volume extends PotentialCommand {
63
64 private static final double rminToSigma = 1.0 / pow(2.0, 1.0 / 6.0);
65
66
67 @Option(names = {"-c", "--connolly"}, paramLabel = "false",
68 description = "Use the Connolly algorithm to compute solvent excluded volume and solvent accessible surface area.")
69 private boolean connolly = false;
70
71
72 @Option(names = {"-m", "--molecular"}, paramLabel = "false",
73 description = "For Connolly, compute molecular volume and surface area (instead of SEV/SASA).")
74 private boolean molecular = false;
75
76
77 @Option(names = {"--vdW", "--vanDerWaals"}, paramLabel = "false",
78 description = "For Connolly, compute van der Waals volume and surface area (instead of SEV/SASA)")
79 private boolean vdW = false;
80
81
82 @Option(names = {"-p", "--probe"}, paramLabel = "1.4",
83 description = "For Connolly, set the exclude radius (SASA) or probe radius (molecular surface). Ignored for vdW.")
84 private double probe = 1.4;
85
86
87 @Option(names = {"-y", "--includeHydrogen"}, paramLabel = "false",
88 description = "Include Hydrogen in calculation volume and surface area.")
89 private boolean includeHydrogen = false;
90
91
92 @Option(names = {"-s", "--sigma"}, paramLabel = "false",
93 description = "Use sigma radii instead of Rmin.")
94 private boolean sigma = false;
95
96
97 @Option(names = {"-o", "--offset"}, paramLabel = "0.0",
98 description = "Add an offset to all atomic radii for GaussVol volume and surface area.")
99 private double offset = 0.0;
100
101
102 @Option(names = {"-v", "--verbose"}, paramLabel = "false",
103 description = "Print out all components of volume of molecule and offset.")
104 private boolean verbose = false;
105
106
107 @Parameters(arity = "1", paramLabel = "file",
108 description = "The atomic coordinate file in PDB or XYZ format.")
109 private String filename = null;
110
111
112 public double totalVolume = 0.0;
113 public double totalSurfaceArea = 0.0;
114
115 public Volume() { super(); }
116 public Volume(FFXBinding binding) { super(binding); }
117 public Volume(String[] args) { super(args); }
118
119 @Override
120 public Volume run() {
121
122 if (!init()) {
123 return this;
124 }
125
126
127 activeAssembly = getActiveAssembly(filename);
128 if (activeAssembly == null) {
129 logger.info(helpString());
130 return this;
131 }
132
133
134 filename = activeAssembly.getFile().getAbsolutePath();
135
136 logger.info("\n Calculating volume and surface area for " + filename);
137
138 Atom[] atoms = activeAssembly.getAtomArray();
139 int nAtoms = atoms.length;
140
141 if (!connolly) {
142
143 double[][] positions = new double[nAtoms][3];
144 int index = 0;
145 for (Atom atom : atoms) {
146 positions[index][0] = atom.getX();
147 positions[index][1] = atom.getY();
148 positions[index][2] = atom.getZ();
149 index++;
150 }
151
152 ForceField forceField = activeAssembly.getForceField();
153 if (includeHydrogen) {
154 forceField.addProperty("GAUSSVOL_HYDROGEN", "true");
155 }
156 if (sigma) {
157 forceField.addProperty("GAUSSVOL_USE_SIGMA", "true");
158 }
159 if (offset > 0.0) {
160 forceField.addProperty("GAUSSVOL_RADII_OFFSET", Double.toString(offset));
161 }
162 if (!forceField.hasProperty("GAUSSVOL_RADII_SCALE")) {
163 forceField.addProperty("GAUSSVOL_RADII_SCALE", Double.toString(1.0));
164 }
165
166
167 ParallelTeam parallelTeam = new ParallelTeam();
168 GaussVol gaussVol = new GaussVol(atoms, activeAssembly.getForceField(), parallelTeam);
169 long time = -System.nanoTime();
170 gaussVol.computeVolumeAndSA(positions);
171 time += System.nanoTime();
172
173 if (verbose) {
174 logger.info(format("\n Maximum depth of overlaps in tree: %d", gaussVol.getMaximumDepth()));
175 logger.info(format(" Total number of overlaps in tree: %d", gaussVol.getTotalNumberOfOverlaps()));
176 double[] radii = gaussVol.getRadii();
177 index = 0;
178 for (Atom atom : atoms) {
179 logger.info(format(" %5d %4s: %6.4f", index, atom.getName(), radii[index]));
180 index++;
181 }
182 }
183
184 logger.info("\n GaussVol Surface Area and Volume\n");
185 if (sigma) {
186 logger.info(format(" Radii: Sigma"));
187 } else {
188 logger.info(format(" Radii: Rmin"));
189 }
190 logger.info(format(" Radii offset: %8.4f (Ang)", offset));
191 logger.info(format(" Include hydrogen: %8b", includeHydrogen));
192 logger.info(format(" Volume: %10.4f (Ang^3)", gaussVol.getVolume()));
193 logger.info(format(" Surface Area: %10.4f (Ang^2)", gaussVol.getSurfaceArea()));
194 logger.info(format(" Time: %10.4f (sec)", time * NS2SEC));
195
196
197 totalVolume = gaussVol.getVolume();
198 totalSurfaceArea = gaussVol.getSurfaceArea();
199 } else {
200
201 double exclude = 0.0;
202
203 if (vdW) {
204
205 exclude = 0.0;
206 probe = 0.0;
207 } else if (!molecular) {
208
209 exclude = probe;
210 probe = 0.0;
211 }
212
213 double[] radii = new double[nAtoms];
214 int index = 0;
215 for (Atom atom : atoms) {
216 radii[index] = atom.getVDWType().radius / 2.0;
217 if (sigma) {
218 radii[index] *= rminToSigma;
219 }
220 boolean hydrogen = atom.isHydrogen();
221 if (!includeHydrogen && hydrogen) {
222 radii[index] = 0.0;
223 }
224 index++;
225 }
226
227
228 int nThreads = 1;
229 ConnollyRegion connollyRegion = new ConnollyRegion(atoms, radii, nThreads);
230
231 connollyRegion.setExclude(exclude);
232
233 connollyRegion.setProbe(probe);
234 long time = -System.nanoTime();
235 connollyRegion.runVolume();
236 time += System.nanoTime();
237
238 double zero = 0.0;
239 if (vdW || (probe == zero && exclude == zero)) {
240 logger.info("\n Connolly van der Waals Surface Area and Volume\n");
241 } else if (!molecular) {
242 logger.info("\n Connolly Solvent Accessible Surface Area and Solvent Excluded Volume\n");
243 logger.info(format(" Exclude radius: %8.4f (Ang)", exclude));
244 } else {
245 logger.info("\n Connolly Molecular Surface Area and Volume");
246 logger.info(format(" Probe radius: %8.4f (Ang)", probe));
247 }
248 if (sigma) {
249 logger.info(format(" Radii: Sigma"));
250 } else {
251 logger.info(format(" Radii: Rmin"));
252 }
253 logger.info(format(" Include hydrogen: %8b", includeHydrogen));
254 logger.info(format(" Volume: %10.4f (Ang^3)", connollyRegion.getVolume()));
255 logger.info(format(" Surface Area: %10.4f (Ang^2)", connollyRegion.getSurfaceArea()));
256 logger.info(format(" Time: %10.4f (sec)", time * NS2SEC));
257
258
259 totalVolume = connollyRegion.getVolume();
260 totalSurfaceArea = connollyRegion.getSurfaceArea();
261 }
262
263 return this;
264 }
265 }