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