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.optimize.TorsionSearch;
43 import ffx.numerics.Potential;
44 import ffx.potential.AssemblyState;
45 import ffx.potential.ForceFieldEnergy;
46 import ffx.potential.MolecularAssembly;
47 import ffx.potential.parsers.SystemFilter;
48 import ffx.potential.parsers.XYZFilter;
49 import ffx.utilities.FFXBinding;
50 import org.apache.commons.io.FileUtils;
51 import org.apache.commons.io.FilenameUtils;
52 import org.apache.commons.math3.util.FastMath;
53 import picocli.CommandLine.Command;
54 import picocli.CommandLine.Option;
55 import picocli.CommandLine.Parameters;
56
57 import java.io.File;
58 import java.io.FileInputStream;
59 import java.io.FileOutputStream;
60 import java.util.Collections;
61 import java.util.List;
62
63 import static java.lang.String.format;
64
65
66
67
68
69
70
71
72
73
74
75
76 @Command(description = " The TorsionSearch command enumerates conformations of a molecule using torsional scans around rotatable bonds.",
77 name = "TorsionSearch")
78 public class TorsionScan extends AlgorithmsCommand {
79
80 @Option(names = {"--th", "--theta"}, paramLabel = "60.0", defaultValue = "60.0",
81 description = "Step size for bond rotations (in Degrees).")
82 private double increment = 60.0;
83
84 @Option(names = {"--saveNumStates", "--sns"}, paramLabel = "10", defaultValue = "10",
85 description = "Save this many of the lowest energy states per worker. This is the default.")
86 private int saveNumStates = 10;
87
88 @Option(names = {"--elimMax", "--em"}, paramLabel = "0", defaultValue = "0",
89 description = "Eliminate bonds where one torsion causes high energies. Reduces the complexity of the search.")
90 private int elimMax = 0;
91
92 @Option(names = {"--saveAll"}, paramLabel = "false", defaultValue = "false",
93 description = "Save out all states. Not recommended for large systems.")
94 private boolean saveAll = false;
95
96 @Option(names = {"--sc", "--staticComparison"}, paramLabel = "false", defaultValue = "false",
97 description = "If set, each bond is rotated independently (faster, but fewer permutations).")
98 private boolean staticCompare;
99
100 @Option(names = {"--eliminationThreshold", "--et"}, paramLabel = "50.0", defaultValue = "50.0",
101 description = "Remove bonds that cause > this energy change during static analysis (kcal/mol).")
102 private double eliminationThreshold = 50.0;
103
104 @Option(names = {"--startIndex"}, paramLabel = "0", defaultValue = "0",
105 description = "Start at this hilbert index.")
106 private long startIndex = 0;
107
108 @Option(names = {"--endIndex"}, paramLabel = "0", defaultValue = "0",
109 description = "End at this hilbert index.")
110 private long endIndex = 0;
111
112 @Parameters(arity = "1", paramLabel = "files",
113 description = "Atomic coordinate file(s) to permute in XYZ format.")
114 String filename = null;
115
116 public TorsionScan() {
117 super();
118 }
119
120 public TorsionScan(FFXBinding binding) {
121 super(binding);
122 }
123
124 public TorsionScan(String[] args) {
125 super(args);
126 }
127
128 @Override
129 public TorsionScan run() {
130 System.setProperty("polarization", "direct");
131
132 if (!init()) {
133 return this;
134 }
135
136 if (filename == null) {
137 logger.info(helpString());
138 return this;
139 }
140
141 activeAssembly = getActiveAssembly(filename);
142 SystemFilter systemFilter = algorithmFunctions.getFilter();
143 MolecularAssembly ma = systemFilter.getActiveMolecularSystem();
144 logger.info("");
145 ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy();
146 double[] x = new double[forceFieldEnergy.getNumberOfVariables()];
147 Comm world = Comm.world();
148 int rank = world.rank();
149 int size = world.size();
150 if (saveAll) {
151 saveNumStates = -1;
152 }
153
154 int numTorsions = (int) FastMath.ceil(360 / increment);
155 TorsionSearch torsionSearch = new TorsionSearch(ma, ma.getMoleculeArray()[0], numTorsions, saveNumStates);
156 if (size > 1 && rank == 0) {
157 torsionSearch.staticAnalysis(elimMax, eliminationThreshold);
158 } else if (size == 1) {
159 torsionSearch.staticAnalysis(elimMax, eliminationThreshold);
160 }
161
162 if (world.size() > 1 && !staticCompare) {
163 torsionSearch.buildWorker(rank, size);
164 torsionSearch.runWorker();
165 } else if (!staticCompare) {
166 endIndex = endIndex == 0 ? torsionSearch.getEnd() : endIndex;
167 torsionSearch.spinTorsions(startIndex, endIndex);
168 }
169
170 List<AssemblyState> states = torsionSearch.getStates();
171 List<Double> energies = torsionSearch.getEnergies();
172 List<Long> hilbertIndices = torsionSearch.getHilbertIndices();
173
174 int count = 0;
175 logger.info("\n Saving " + states.size() + " states.");
176 String extension = "_rot.arc";
177 if (world.size() > 1) {
178 extension = "_rank" + world.rank() + extension;
179 }
180 File saveLocation = new File(FilenameUtils.removeExtension(filename) + extension);
181 logger.info(" Logging structures into: " + saveLocation);
182 XYZFilter xyzFilter = new XYZFilter(saveLocation,
183 activeAssembly,
184 activeAssembly.getForceField(),
185 activeAssembly.getProperties());
186 while (!states.isEmpty()) {
187 AssemblyState assembly = states.get(0);
188 states.remove(0);
189 assembly.revertState();
190 forceFieldEnergy.getCoordinates(x);
191 double e = energies.get(0);
192 energies.remove(0);
193 long hilbertIndex = hilbertIndices.get(0);
194 hilbertIndices.remove(0);
195 logger.info(format(" Writing to file. Configuration #%-6d energy: %-12.5f Hilbert index: %-15d",
196 (count + 1),
197 e,
198 hilbertIndex));
199 xyzFilter.writeFile(saveLocation, true);
200 count++;
201 }
202 logger.info("\n -1 indices come from static torsion scan.");
203
204 File key = new File(FilenameUtils.removeExtension(filename) + ".key");
205 File properties = new File(FilenameUtils.removeExtension(filename) + ".properties");
206 try {
207 if (key.exists()) {
208 File keyComparison = new File(FilenameUtils.removeExtension(filename) + "_rot.key");
209 if (keyComparison.createNewFile()) {
210 FileUtils.copyFile(key, keyComparison);
211 }
212 } else if (properties.exists()) {
213 File propertiesComparison = new File(FilenameUtils.removeExtension(filename) + "_rot.properties");
214 if (propertiesComparison.createNewFile()) {
215 FileUtils.copyFile(properties, propertiesComparison);
216 }
217 } else {
218 logger.info(" No key or properties file found.");
219 }
220 } catch (Exception e) {
221 e.printStackTrace();
222 }
223
224 if (world.rank() == 0 && world.size() > 1) {
225 logger.info("\n Combining all rank files into one file.");
226 File combined = new File(FilenameUtils.removeExtension(filename) + "_rot.arc");
227 try {
228 if (combined.createNewFile()) {
229 try (FileOutputStream fos = new FileOutputStream(combined)) {
230 for (int i = 0; i < world.size(); i++) {
231 File rankFile = new File(FilenameUtils.removeExtension(filename) + "_rank" + i + "_rot.arc");
232 if (rankFile.exists()) {
233 try (FileInputStream fis = new FileInputStream(rankFile)) {
234 byte[] buffer = new byte[1024];
235 int length;
236 while ((length = fis.read(buffer)) > 0) {
237 fos.write(buffer, 0, length);
238 }
239 }
240 }
241 }
242 }
243 }
244 } catch (Exception e) {
245 e.printStackTrace();
246 }
247 }
248 return this;
249 }
250
251 @Override
252 public List<Potential> getPotentials() {
253 return Collections.emptyList();
254 }
255 }