main.cpp
Go to the documentation of this file.
Code
1
2//
3// d-SEAMS molecular dynamics analysis engine code
4// Copyright (C) <2018--present> Amrita Goswami, Rohit Goswami
5// amrita16thaug646[at]gmail.com, r95g10[at]gmail.com
6//
7// This program is free software: you can redistribute it and/or modify
8// it under the terms of the MIT License as published by
9// the Open Source Initiative.
10//
11// A copy of the MIT License is included in the LICENSE file of this repository.
12// You should have received a copy of the MIT License along with this program.
13// If not, see <https://opensource.org/licenses/MIT>.
14//
16
17// Standard Library
18#include <array>
19#include <cstdlib>
20#include <ctime>
21#include <sstream>
22#include <string>
23
24// Internal Libraries
25#include "opt_parser.h"
26
27// Newer pointCloud
28#include <bond.hpp>
29#include <bop.hpp>
30#include <bulkTUM.hpp>
31#include <cluster.hpp>
32#include <franzblau.hpp>
33#include <generic.hpp>
34#include <mol_sys.hpp>
35#include <neighbours.hpp>
36#include <rdf2d.hpp>
37#include <ring.hpp>
38#include <seams_input.hpp>
39#include <seams_output.hpp>
40#include <topo_bulk.hpp>
41#include <topo_one_dim.hpp>
42#include <topo_two_dim.hpp>
43#include <selection.hpp>
44
45// Externally bundled-input libraries
46// #include <cxxopts.hpp>
47#include <rang.hpp>
48#include <sol/sol.hpp>
49
50#include <fmt/core.h>
51#include <yaml-cpp/yaml.h>
52
53int main(int argc, char *argv[]) {
54 // Parse Things
55 auto result = parse(argc, argv);
56 auto &arguments = result.arguments();
57 // Initialize yaml config
58 YAML::Node config = YAML::LoadFile(result["c"].as<std::string>());
59 // This is a dummy used to figure out the order of options (cmd > yml)
60 std::string script, tFile;
61 // Initialize Lua
62 sol::state lua;
63 // Use all libraries
64 lua.open_libraries();
65 // Get the trajectory string
66 if (config["trajectory"]) {
67 tFile = config["trajectory"].as<std::string>();
68 } // end of getting the trajectory
69 // Get variable file string
70 std::string vars = config["variables"].as<std::string>();
71 // --------------------------------------
72 // Structure determination block for TWO-DIMENSIONAL ICE
73 if (config["topoTwoDim"]["use"].as<bool>()) {
74 // Use the variables script
75 lua.script_file(vars);
76 // -----------------
77 // Variables which must be declared in C++
78 //
79 // Newer pointCloud (rescloud -> ice structure, solcloud -> largest cluster)
81 // Some neighbor lists
82 std::vector<std::vector<int>> nList, hbnList;
83 // For the list of all rings (of all sizes)
84 std::vector<std::vector<int>> ringsAllSizes;
85 std::vector<std::vector<int>> rings;
86 // RDF stuff
87 std::vector<double> rdfValues; // RDF vector
88 // -----------------
89 // This section basically only registers functions and handles the rest in
90 // lua Use the functions defined here
91 auto lscript = lua.get<std::string>("functionScript");
92 // Transfer variables to lua
93 lua["doBOP"] = config["bulk"]["use"].as<bool>();
94 lua["topoOneDim"] = config["topoOneDim"]["use"].as<bool>();
95 lua["topoTwoDim"] = config["topoTwoDim"]["use"].as<bool>();
96 lua["topoBulk"] = config["bulk"]["use"].as<bool>();
97 //
98 lua["nList"] = &nList;
99 lua["hbnList"] = &hbnList;
100 lua["resCloud"] = &resCloud;
101 lua["trajectory"] = tFile;
102 // Confined ice stuff
103 lua["ringsAllSizes"] = &rings;
104 // RDF stuff
105 lua["rdf"] = &rdfValues;
106 // -----------------
107 // Register functions
108 //
109 // Writing stuff
110 // Generic requirements
111 lua.set_function("readFrameOnlyOne", sinp::readLammpsTrjreduced);
112 lua.set_function("readFrameOnlyOneAllAtoms", sinp::readLammpsTrj); // reads in all atoms regardless of type
113 lua.set_function("neighborList", nneigh::neighListO);
114 // -----------------
115 // Topological Network Method Specific Functions
116 // Generic requirements (read in only inside the slice)
117 lua.set_function("getHbondNetwork", bond::populateHbonds);
118 lua.set_function("bondNetworkByIndex", nneigh::neighbourListByIndex);
119 // -----------------
120 // Primitive rings
121 lua.set_function("getPrimitiveRings", primitive::ringNetwork);
122 // -----------------
123 // Quasi-two-dimensional ice
124 lua.set_function("ringAnalysis", ring::polygonRingAnalysis);
125 // --------------------------
126 // RDF functions
127 lua.set_function("calcRDF", rdf2::rdf2Danalysis_AA);
128 // --------------------------
129 // Use the script
130 lua.script_file(lscript);
131 // --------------------------
132
133 } // end of two-dimensional ice block
134 // --------------------------------------
135 // Structure determination block for ONE-DIMENSIONAL ICE
136 if (config["topoOneDim"]["use"].as<bool>()) {
137 // Use the script
138 lua.script_file(vars);
139 // -----------------
140 // Variables which must be declared in C++
141 //
142 // Newer pointCloud (rescloud -> ice structure, solcloud -> largest cluster)
144 molSys::PointCloud<molSys::Point<double>, double> oCloud; // O atom pointCloud
145 molSys::PointCloud<molSys::Point<double>, double> hCloud; // H atom pointCloud
146 // Some neighbor
147 std::vector<std::vector<int>> nList, hbnList;
148 // For the list of all rings (of all sizes)
149 std::vector<std::vector<int>> ringsAllSizes;
150 std::vector<std::vector<int>> rings;
151 int atomID;
152 // -----------------
153 // This section basically only registers functions and handles the rest in
154 // lua Use the functions defined here
155 auto lscript = lua.get<std::string>("functionScript");
156 // Transfer variables to lua
157 lua["doBOP"] = config["bulk"]["use"].as<bool>();
158 lua["topoOneDim"] = config["topoOneDim"]["use"].as<bool>();
159 lua["topoTwoDim"] = config["topoTwoDim"]["use"].as<bool>();
160 lua["topoBulk"] = config["bulk"]["use"].as<bool>();
161 //
162 lua["nList"] = &nList;
163 lua["hbnList"] = &hbnList;
164 lua["resCloud"] = &resCloud;
165 lua["oCloud"] = &oCloud;
166 lua["hCloud"] = &hCloud;
167 lua["trajectory"] = tFile;
168 // Confined ice stuff
169 lua["ringsAllSizes"] = &rings;
170 lua["lowestAtomID"] = &atomID;
171 // Register functions
172 //
173 // Writing stuff
174 // Generic requirements
175 lua.set_function("readFrameOnlyOne", sinp::readLammpsTrjreduced);
176 lua.set_function("readFrameOnlyOneAllAtoms", sinp::readLammpsTrj); // reads in all atoms regardless of type
177 lua.set_function("neighborList", nneigh::neighListO);
178 // -----------------
179 // Topological Network Method Specific Functions
180 // Generic requirements (read in only inside the slice)
181 lua.set_function("getHbondNetwork", bond::populateHbonds);
182 lua.set_function("bondNetworkByIndex", nneigh::neighbourListByIndex);
183 // -----------------
184 // Primitive rings
185 lua.set_function("getPrimitiveRings", primitive::ringNetwork);
186 // -----------------
187 // Quasi-one-dimensional ice
188 lua.set_function("prismAnalysis", ring::prismAnalysis);
189 // --------------------------
190 // Use the script
191 lua.script_file(lscript);
192 // --------------------------
193
194 } // end of one-dimensional ice block
195 // --------------------------------------
196 // Ice Structure Determination for BULK ICE
197 if (config["bulk"]["use"].as<bool>()) {
198 // Use the variables script
199 lua.script_file(vars);
200 // Variables which must be declared in C++
201 //
202 // Newer pointCloud (rescloud -> ice structure, solcloud -> largest
203 molSys::PointCloud<molSys::Point<double>, double> resCloud, solCloud;
204 // PointCloud for O atoms and H atoms separately
205 molSys::PointCloud<molSys::Point<double>, double> oCloud, hCloud;
206 // Some neighbor
207 std::vector<std::vector<int>> nList,
208 hbnList; // Neighbour lists (by cutoff and hydrogen-bonded neighbour
209 // lists)
210 std::vector<std::vector<int>>
211 iceList; // Neighbour list for the largest ice cluster
212 // For averaged q6
213 std::vector<double> avgQ6;
214 // For the list of all rings (of all sizes)
215 std::vector<std::vector<int>> ringsAllSizes;
216 std::vector<std::vector<int>> rings;
217 // -----------------
218 // Variables defined in C++ specific to confined systems
219
220 // This section basically only registers functions and handles the rest in
221 // lua lua Use the functions defined here
222 auto lscript = lua.get<std::string>("functionScript");
223 // Transfer variables to lua
224 lua["doBOP"] = config["bulk"]["bondOrderParameters"].as<bool>();
225 lua["topoOneDim"] = config["topoOneDim"]["use"].as<bool>();
226 lua["topoTwoDim"] = config["topoTwoDim"]["use"].as<bool>();
227 lua["topoBulk"] = config["bulk"]["topologicalNetworkCriterion"].as<bool>();
228 //
229 lua["nList"] = &nList;
230 lua["hbnList"] = &hbnList;
231 lua["iceNeighbourList"] = &iceList;
232 lua["resCloud"] = &resCloud;
233 lua["oCloud"] = &oCloud;
234 lua["hCloud"] = &hCloud;
235 lua["clusterCloud"] = &solCloud;
236 lua["avgQ6"] = &avgQ6;
237 lua["trajectory"] = tFile;
238 // Confined ice stuff
239 lua["ringsAllSizes"] = &rings;
240 // Register functions
241 //
242 // Writing stuff
243 lua.set_function("writeDump", sout::writeDump);
244 lua.set_function("writeHistogram", sout::writeHisto);
245 // Generic requirements
246 lua.set_function("readFrame", sinp::readLammpsTrjO);
247 lua.set_function("neighborList", nneigh::neighListO);
248 // CHILL+ and modifications
249 lua.set_function("chillPlus_cij", chill::getCorrelPlus);
250 lua.set_function("chillPlus_iceType", chill::getIceTypePlus);
251 // CHILL functions
252 lua.set_function("chill_cij", chill::getCorrel);
253 lua.set_function("chill_iceType", chill::getIceType);
254 // Reclassify using q6
255 lua.set_function("averageQ6", chill::getq6);
256 lua.set_function("modifyChill", chill::reclassifyWater);
257 lua.set_function("percentage_Ice", chill::printIceType);
258 // Largest ice cluster
259 lua.set_function("clusterAnalysis", clump::clusterAnalysis);
260 lua.set_function("recenterCluster", clump::recenterClusterCloud);
261 // -----------------
262 // Selection Functions
263 // Function for getting an output PointCloud of a particular atom type from an existing PointCloud
264 lua.set_function("getPointCloudAtomsOfOneAtomType", gen::getPointCloudOneAtomType);
265 lua.set_function("selectInSingleSlice", gen::moleculesInSingleSlice);
266 lua.set_function("selectEdgeAtomsInRingsWithinSlice", ring::getEdgeMoleculesInRings);
267 lua.set_function("selectAtomsInSliceWithRingEdgeAtoms", ring::printSliceGetEdgeMoleculesInRings);
268 // -----------------
269 // Topological Network Methods
270 // Generic requirements (read in only inside the slice)
271 lua.set_function("readFrameOnlyOne", sinp::readLammpsTrjreduced);
272 lua.set_function("readFrameOnlyOneAllAtoms", sinp::readLammpsTrj); // reads in all atoms regardless of type
273 lua.set_function("getHbondNetwork", bond::populateHbonds);
274 lua.set_function("getHbondNetworkFromClouds", bond::populateHbondsWithInputClouds);
275 lua.set_function("bondNetworkByIndex", nneigh::neighbourListByIndex);
276 // -----------------
277 // Primitive rings
278 lua.set_function("getPrimitiveRings", primitive::ringNetwork);
279 // Function for just getting and writing out the ring numbers
280 lua.set_function("bulkRingNumberAnalysis", ring::bulkPolygonRingAnalysis);
281 // -----------------
282 // Bulk ice, using the topological network criterion
283 lua.set_function("bulkTopologicalNetworkCriterion", ring::topoBulkAnalysis);
284 // --------------------------
285 // Bulk ice, using the TUM (Topological Unit Matching Criterion). No need to
286 // use bulkTopologicalNetworkCriterion if you use this function
287 lua.set_function("bulkTopoUnitMatching", tum3::topoUnitMatchingBulk);
288 // --------------------------
289 // Use the script
290 lua.script_file(lscript);
291 // --------------------------
292
293 } // end of bulk ice structure determination block
294 // --------------------------------------
295
296 std::cout << rang::style::bold
297 << fmt::format("Welcome to the Black Parade.\nYou ran:-\n")
298 << rang::style::reset
299 << fmt::format("\nBulk Ice Analysis: {}",
300 config["bulk"]["use"].as<bool>())
301 << fmt::format("\nQuasi-one-dimensional Ice Analysis: {}",
302 config["topoOneDim"]["use"].as<bool>())
303 << fmt::format("\nQuasi-two-dimensional Ice Analysis: {}",
304 config["topoTwoDim"]["use"].as<bool>())
305 << "\n";
306
307 return 0;
308}
File for bond-related analyses (hydrogen bonds, bonded atoms for data file write-outs etc....
File for the bond order parameter analysis.
File containing functions used specific to bulk topological unit matching (TUM) criterion.
File for the bond order parameter analysis.
File for generating shortest-path rings according to the Franzblau algorithm.
File for containing generic or common functions.
std::vector< std::vector< int > > populateHbonds(std::string filename, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList, int targetFrame, int Htype)
Definition bond.cpp:174
std::vector< std::vector< int > > populateHbondsWithInputClouds(molSys::PointCloud< molSys::Point< double >, double > *yCloud, molSys::PointCloud< molSys::Point< double >, double > *hCloud, std::vector< std::vector< int > > nList)
Definition bond.cpp:343
molSys::PointCloud< molSys::Point< double >, double > getIceType(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="chill.txt")
Classifies each atom according to the CHILL algorithm.
Definition bop.cpp:505
int printIceType(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="superChill.txt")
Prints out the iceType for a particular frame onto the terminal.
Definition bop.cpp:1060
molSys::PointCloud< molSys::Point< double >, double > getCorrelPlus(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList, bool isSlice=false)
Gets c_ij and then classifies bond types according to the CHILL+ algorithm.
Definition bop.cpp:602
molSys::PointCloud< molSys::Point< double >, double > getCorrel(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList, bool isSlice=false)
Definition bop.cpp:299
molSys::PointCloud< molSys::Point< double >, double > reclassifyWater(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< double > *q6)
Definition bop.cpp:1011
std::vector< double > getq6(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList, bool isSlice=false)
Definition bop.cpp:868
molSys::PointCloud< molSys::Point< double >, double > getIceTypePlus(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList, std::string path, int firstFrame, bool isSlice=false, std::string outputFileName="chillPlus.txt")
Classifies each atom according to the CHILL+ algorithm.
Definition bop.cpp:755
int recenterClusterCloud(molSys::PointCloud< molSys::Point< double >, double > *iceCloud, std::vector< std::vector< int > > nList)
Recenters the coordinates of a pointCloud.
Definition cluster.cpp:398
int clusterAnalysis(std::string path, molSys::PointCloud< molSys::Point< double >, double > *iceCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList, std::vector< std::vector< int > > &iceNeighbourList, double cutoff, int firstFrame, std::string bopAnalysis="q6")
Definition cluster.cpp:312
std::vector< std::vector< int > > neighbourListByIndex(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList)
std::vector< std::vector< int > > neighListO(double rcutoff, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeI)
std::vector< std::vector< int > > ringNetwork(std::vector< std::vector< int > > nList, int maxDepth)
Definition franzblau.cpp:32
int topoBulkAnalysis(std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, bool onlyTetrahedral=true)
int bulkPolygonRingAnalysis(std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int maxDepth, int firstFrame)
Definition topo_bulk.cpp:44
int rdf2Danalysis_AA(std::string path, std::vector< double > *rdfValues, molSys::PointCloud< molSys::Point< double >, double > *yCloud, double cutoff, double binwidth, int firstFrame, int finalFrame)
Definition rdf2d.cpp:45
int topoUnitMatchingBulk(std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int firstFrame, bool printClusters, bool onlyTetrahedral)
Definition bulkTUM.cpp:52
molSys::PointCloud< molSys::Point< double >, double > getPointCloudOneAtomType(molSys::PointCloud< molSys::Point< double >, double > *yCloud, molSys::PointCloud< molSys::Point< double >, double > *outCloud, int atomTypeI, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
Definition selection.cpp:37
int polygonRingAnalysis(std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int maxDepth, double sheetArea, int firstFrame)
void getEdgeMoleculesInRings(std::vector< std::vector< int > > rings, molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh, bool identicalCloud=false)
int prismAnalysis(std::string path, std::vector< std::vector< int > > rings, std::vector< std::vector< int > > nList, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int maxDepth, int *atomID, int firstFrame, int currentFrame, bool doShapeMatching=false)
void moleculesInSingleSlice(molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool clearPreviousSliceSelection=true, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
void printSliceGetEdgeMoleculesInRings(std::string path, std::vector< std::vector< int > > rings, molSys::PointCloud< molSys::Point< double >, double > *oCloud, molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::array< double, 3 > coordLow, std::array< double, 3 > coordHigh, bool identicalCloud=false)
molSys::PointCloud< molSys::Point< double >, double > readLammpsTrj(std::string filename, int targetFrame, molSys::PointCloud< molSys::Point< double >, double > *yCloud, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
molSys::PointCloud< molSys::Point< double >, double > readLammpsTrjreduced(std::string filename, int targetFrame, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeI, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
molSys::PointCloud< molSys::Point< double >, double > readLammpsTrjO(std::string filename, int targetFrame, molSys::PointCloud< molSys::Point< double >, double > *yCloud, int typeO, bool isSlice=false, std::array< double, 3 > coordLow=std::array< double, 3 >{0, 0, 0}, std::array< double, 3 > coordHigh=std::array< double, 3 >{0, 0, 0})
int main(int argc, char *argv[])
Definition main.cpp:53
The main molecular system handler.
int writeHisto(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::vector< std::vector< int > > nList, std::vector< double > avgQ6)
int writeDump(molSys::PointCloud< molSys::Point< double >, double > *yCloud, std::string path, std::string outFile)
Generic function for writing out to a dump file.
Header file for neighbour list generation.
cxxopts::ParseResult parse(int argc, char *argv[])
File containing functions used to calculate the in-plane radial distribution functions.
File containing common functions used by bulk and confined topological network critera.
File for functions that read in files).
File containing functions used to define 'selections' in a given range, using ring information.
This contains a collection of points; contains information for a particular frame.
Definition mol_sys.hpp:170
File containing functions used specific to bulk topological network critera.
File containing functions used specific to quasi-one-dimensional topological network critera (the pri...