| 1 | /** \file builder.cpp
|
|---|
| 2 | *
|
|---|
| 3 | * By stating absolute positions or binding angles and distances atomic positions of a molecule can be constructed.
|
|---|
| 4 | * The output is the complete configuration file for PCP for direct use.
|
|---|
| 5 | * Features:
|
|---|
| 6 | * -# Atomic data is retrieved from a file, if not found requested and stored there for later re-use
|
|---|
| 7 | * -# step-by-step construction of the molecule beginning either at a centre of with a certain atom
|
|---|
| 8 | *
|
|---|
| 9 | */
|
|---|
| 10 |
|
|---|
| 11 | /*! \mainpage Molecuilder - a molecular set builder
|
|---|
| 12 | *
|
|---|
| 13 | * This introductory shall briefly make aquainted with the program, helping in installing and a first run.
|
|---|
| 14 | *
|
|---|
| 15 | * \section about About the Program
|
|---|
| 16 | *
|
|---|
| 17 | * Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
|
|---|
| 18 | * atoms making up an molecule by the successive statement of binding angles and distances and referencing to
|
|---|
| 19 | * already constructed atoms.
|
|---|
| 20 | *
|
|---|
| 21 | * A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
|
|---|
| 22 | * molecular dynamics implementation.
|
|---|
| 23 | *
|
|---|
| 24 | * \section install Installation
|
|---|
| 25 | *
|
|---|
| 26 | * Installation should without problems succeed as follows:
|
|---|
| 27 | * -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
|
|---|
| 28 | * -# make
|
|---|
| 29 | * -# make install
|
|---|
| 30 | *
|
|---|
| 31 | * Further useful commands are
|
|---|
| 32 | * -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
|
|---|
| 33 | * -# make doxygen-doc: Creates these html pages out of the documented source
|
|---|
| 34 | *
|
|---|
| 35 | * \section run Running
|
|---|
| 36 | *
|
|---|
| 37 | * The program can be executed by running: ./molecuilder
|
|---|
| 38 | *
|
|---|
| 39 | * Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
|
|---|
| 40 | * it is created and any given data on elements of the periodic table will be stored therein and re-used on
|
|---|
| 41 | * later re-execution.
|
|---|
| 42 | *
|
|---|
| 43 | * \section ref References
|
|---|
| 44 | *
|
|---|
| 45 | * For the special configuration file format, see the documentation of pcp.
|
|---|
| 46 | *
|
|---|
| 47 | */
|
|---|
| 48 |
|
|---|
| 49 |
|
|---|
| 50 | #include <boost/bind.hpp>
|
|---|
| 51 |
|
|---|
| 52 | using namespace std;
|
|---|
| 53 |
|
|---|
| 54 | #include "analysis_correlation.hpp"
|
|---|
| 55 | #include "atom.hpp"
|
|---|
| 56 | #include "bond.hpp"
|
|---|
| 57 | #include "bondgraph.hpp"
|
|---|
| 58 | #include "boundary.hpp"
|
|---|
| 59 | #include "config.hpp"
|
|---|
| 60 | #include "element.hpp"
|
|---|
| 61 | #include "ellipsoid.hpp"
|
|---|
| 62 | #include "helpers.hpp"
|
|---|
| 63 | #include "leastsquaremin.hpp"
|
|---|
| 64 | #include "linkedcell.hpp"
|
|---|
| 65 | #include "log.hpp"
|
|---|
| 66 | #include "memoryusageobserverunittest.hpp"
|
|---|
| 67 | #include "molecule.hpp"
|
|---|
| 68 | #include "periodentafel.hpp"
|
|---|
| 69 | #include "UIElements/UIFactory.hpp"
|
|---|
| 70 | #include "UIElements/MainWindow.hpp"
|
|---|
| 71 | #include "UIElements/Dialog.hpp"
|
|---|
| 72 | #include "Menu/ActionMenuItem.hpp"
|
|---|
| 73 | #include "Actions/ActionRegistry.hpp"
|
|---|
| 74 | #include "Actions/MethodAction.hpp"
|
|---|
| 75 | #include "Actions/small_actions.hpp"
|
|---|
| 76 |
|
|---|
| 77 |
|
|---|
| 78 | /** Parses the command line options.
|
|---|
| 79 | * \param argc argument count
|
|---|
| 80 | * \param **argv arguments array
|
|---|
| 81 | * \param *molecules list of molecules structure
|
|---|
| 82 | * \param *periode elements structure
|
|---|
| 83 | * \param configuration config file structure
|
|---|
| 84 | * \param *ConfigFileName pointer to config file name in **argv
|
|---|
| 85 | * \param *PathToDatabases pointer to db's path in **argv
|
|---|
| 86 | * \return exit code (0 - successful, all else - something's wrong)
|
|---|
| 87 | */
|
|---|
| 88 | static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode,\
|
|---|
| 89 | config& configuration, char *&ConfigFileName)
|
|---|
| 90 | {
|
|---|
| 91 | Vector x,y,z,n; // coordinates for absolute point in cell volume
|
|---|
| 92 | double *factor; // unit factor if desired
|
|---|
| 93 | ifstream test;
|
|---|
| 94 | ofstream output;
|
|---|
| 95 | string line;
|
|---|
| 96 | atom *first;
|
|---|
| 97 | bool SaveFlag = false;
|
|---|
| 98 | int ExitFlag = 0;
|
|---|
| 99 | int j;
|
|---|
| 100 | double volume = 0.;
|
|---|
| 101 | enum ConfigStatus configPresent = absent;
|
|---|
| 102 | clock_t start,end;
|
|---|
| 103 | int argptr;
|
|---|
| 104 | molecule *mol = NULL;
|
|---|
| 105 | string BondGraphFileName("");
|
|---|
| 106 | int verbosity = 0;
|
|---|
| 107 | strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1);
|
|---|
| 108 |
|
|---|
| 109 | if (argc > 1) { // config file specified as option
|
|---|
| 110 | // 1. : Parse options that just set variables or print help
|
|---|
| 111 | argptr = 1;
|
|---|
| 112 | do {
|
|---|
| 113 | if (argv[argptr][0] == '-') {
|
|---|
| 114 | Log() << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";
|
|---|
| 115 | argptr++;
|
|---|
| 116 | switch(argv[argptr-1][1]) {
|
|---|
| 117 | case 'h':
|
|---|
| 118 | case 'H':
|
|---|
| 119 | case '?':
|
|---|
| 120 | Log() << Verbose(0) << "MoleCuilder suite" << endl << "==================" << endl << endl;
|
|---|
| 121 | Log() << Verbose(0) << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
|
|---|
| 122 | Log() << Verbose(0) << "or simply " << argv[0] << " without arguments for interactive session." << endl;
|
|---|
| 123 | Log() << Verbose(0) << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
|
|---|
| 124 | Log() << Verbose(0) << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
|
|---|
| 125 | Log() << Verbose(0) << "\t-b xx xy xz yy yz zz\tCenter atoms in domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
|
|---|
| 126 | Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
|
|---|
| 127 | Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
|
|---|
| 128 | Log() << Verbose(0) << "\t-C\tPair Correlation analysis." << endl;
|
|---|
| 129 | Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
|
|---|
| 130 | Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
|
|---|
| 131 | Log() << Verbose(0) << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
|
|---|
| 132 | Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
|
|---|
| 133 | Log() << Verbose(0) << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
|
|---|
| 134 | Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
|
|---|
| 135 | Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
|
|---|
| 136 | Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
|
|---|
| 137 | Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
|
|---|
| 138 | Log() << Verbose(0) << "\t-M <basis>\tSetting basis to store to MPQC config files." << endl;
|
|---|
| 139 | Log() << Verbose(0) << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
|
|---|
| 140 | Log() << Verbose(0) << "\t-N <radius> <file>\tGet non-convex-envelope." << endl;
|
|---|
| 141 | Log() << Verbose(0) << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
|
|---|
| 142 | Log() << Verbose(0) << "\t-O\tCenter atoms in origin." << endl;
|
|---|
| 143 | Log() << Verbose(0) << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
|
|---|
| 144 | Log() << Verbose(0) << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
|
|---|
| 145 | Log() << Verbose(0) << "\t-r <id>\t\tRemove an atom with given id." << endl;
|
|---|
| 146 | Log() << Verbose(0) << "\t-R <id> <radius>\t\tRemove all atoms out of sphere around a given one." << endl;
|
|---|
| 147 | Log() << Verbose(0) << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
|
|---|
| 148 | Log() << Verbose(0) << "\t-S <file> Store temperatures from the config file in <file>." << endl;
|
|---|
| 149 | Log() << Verbose(0) << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
|
|---|
| 150 | Log() << Verbose(0) << "\t-T x1 x2 x3\tTranslate periodically all atoms by this vector (x1,x2,x3)." << endl;
|
|---|
| 151 | Log() << Verbose(0) << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
|
|---|
| 152 | Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl;
|
|---|
| 153 | Log() << Verbose(0) << "\t-V\t\tGives version information." << endl;
|
|---|
| 154 | Log() << Verbose(0) << "Note: config files must not begin with '-' !" << endl;
|
|---|
| 155 | return (1);
|
|---|
| 156 | break;
|
|---|
| 157 | case 'v':
|
|---|
| 158 | while (argv[argptr-1][verbosity+1] == 'v') {
|
|---|
| 159 | verbosity++;
|
|---|
| 160 | }
|
|---|
| 161 | setVerbosity(verbosity);
|
|---|
| 162 | Log() << Verbose(0) << "Setting verbosity to " << verbosity << "." << endl;
|
|---|
| 163 | break;
|
|---|
| 164 | case 'V':
|
|---|
| 165 | Log() << Verbose(0) << argv[0] << " " << VERSIONSTRING << endl;
|
|---|
| 166 | Log() << Verbose(0) << "Build your own molecule position set." << endl;
|
|---|
| 167 | return (1);
|
|---|
| 168 | break;
|
|---|
| 169 | case 'e':
|
|---|
| 170 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
|---|
| 171 | eLog() << Verbose(0) << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
|
|---|
| 172 | performCriticalExit();
|
|---|
| 173 | } else {
|
|---|
| 174 | Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl;
|
|---|
| 175 | strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1);
|
|---|
| 176 | argptr+=1;
|
|---|
| 177 | }
|
|---|
| 178 | break;
|
|---|
| 179 | case 'g':
|
|---|
| 180 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
|---|
| 181 | eLog() << Verbose(0) << "Not enough or invalid arguments for specifying bond length table: -g <table file>" << endl;
|
|---|
| 182 | performCriticalExit();
|
|---|
| 183 | } else {
|
|---|
| 184 | BondGraphFileName = argv[argptr];
|
|---|
| 185 | Log() << Verbose(0) << "Using " << BondGraphFileName << " as bond length table." << endl;
|
|---|
| 186 | argptr+=1;
|
|---|
| 187 | }
|
|---|
| 188 | break;
|
|---|
| 189 | case 'n':
|
|---|
| 190 | Log() << Verbose(0) << "I won't parse trajectories." << endl;
|
|---|
| 191 | configuration.FastParsing = true;
|
|---|
| 192 | break;
|
|---|
| 193 | default: // no match? Step on
|
|---|
| 194 | argptr++;
|
|---|
| 195 | break;
|
|---|
| 196 | }
|
|---|
| 197 | } else
|
|---|
| 198 | argptr++;
|
|---|
| 199 | } while (argptr < argc);
|
|---|
| 200 |
|
|---|
| 201 | // 3a. Parse the element database
|
|---|
| 202 | if (periode->LoadPeriodentafel(configuration.databasepath)) {
|
|---|
| 203 | Log() << Verbose(0) << "Element list loaded successfully." << endl;
|
|---|
| 204 | //periode->Output();
|
|---|
| 205 | } else {
|
|---|
| 206 | Log() << Verbose(0) << "Element list loading failed." << endl;
|
|---|
| 207 | return 1;
|
|---|
| 208 | }
|
|---|
| 209 | // 3b. Find config file name and parse if possible, also BondGraphFileName
|
|---|
| 210 | if (argv[1][0] != '-') {
|
|---|
| 211 | // simply create a new molecule, wherein the config file is loaded and the manipulation takes place
|
|---|
| 212 | Log() << Verbose(0) << "Config file given." << endl;
|
|---|
| 213 | test.open(argv[1], ios::in);
|
|---|
| 214 | if (test == NULL) {
|
|---|
| 215 | //return (1);
|
|---|
| 216 | output.open(argv[1], ios::out);
|
|---|
| 217 | if (output == NULL) {
|
|---|
| 218 | Log() << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
|
|---|
| 219 | configPresent = absent;
|
|---|
| 220 | } else {
|
|---|
| 221 | Log() << Verbose(0) << "Empty configuration file." << endl;
|
|---|
| 222 | ConfigFileName = argv[1];
|
|---|
| 223 | configPresent = empty;
|
|---|
| 224 | output.close();
|
|---|
| 225 | }
|
|---|
| 226 | } else {
|
|---|
| 227 | test.close();
|
|---|
| 228 | ConfigFileName = argv[1];
|
|---|
| 229 | Log() << Verbose(1) << "Specified config file found, parsing ... ";
|
|---|
| 230 | switch (configuration.TestSyntax(ConfigFileName, periode)) {
|
|---|
| 231 | case 1:
|
|---|
| 232 | Log() << Verbose(0) << "new syntax." << endl;
|
|---|
| 233 | configuration.Load(ConfigFileName, BondGraphFileName, periode, molecules);
|
|---|
| 234 | configPresent = present;
|
|---|
| 235 | break;
|
|---|
| 236 | case 0:
|
|---|
| 237 | Log() << Verbose(0) << "old syntax." << endl;
|
|---|
| 238 | configuration.LoadOld(ConfigFileName, BondGraphFileName, periode, molecules);
|
|---|
| 239 | configPresent = present;
|
|---|
| 240 | break;
|
|---|
| 241 | default:
|
|---|
| 242 | Log() << Verbose(0) << "Unknown syntax or empty, yet present file." << endl;
|
|---|
| 243 | configPresent = empty;
|
|---|
| 244 | }
|
|---|
| 245 | }
|
|---|
| 246 | } else
|
|---|
| 247 | configPresent = absent;
|
|---|
| 248 | // set mol to first active molecule
|
|---|
| 249 | if (molecules->ListOfMolecules.size() != 0) {
|
|---|
| 250 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
|
|---|
| 251 | if ((*ListRunner)->ActiveFlag) {
|
|---|
| 252 | mol = *ListRunner;
|
|---|
| 253 | break;
|
|---|
| 254 | }
|
|---|
| 255 | }
|
|---|
| 256 | if (mol == NULL) {
|
|---|
| 257 | mol = new molecule(periode);
|
|---|
| 258 | mol->ActiveFlag = true;
|
|---|
| 259 | molecules->insert(mol);
|
|---|
| 260 | }
|
|---|
| 261 |
|
|---|
| 262 | // 4. parse again through options, now for those depending on elements db and config presence
|
|---|
| 263 | argptr = 1;
|
|---|
| 264 | do {
|
|---|
| 265 | Log() << Verbose(0) << "Current Command line argument: " << argv[argptr] << "." << endl;
|
|---|
| 266 | if (argv[argptr][0] == '-') {
|
|---|
| 267 | argptr++;
|
|---|
| 268 | if ((configPresent == present) || (configPresent == empty)) {
|
|---|
| 269 | switch(argv[argptr-1][1]) {
|
|---|
| 270 | case 'p':
|
|---|
| 271 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 272 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
|---|
| 273 | ExitFlag = 255;
|
|---|
| 274 | eLog() << Verbose(0) << "Not enough arguments for parsing: -p <xyz file>" << endl;
|
|---|
| 275 | performCriticalExit();
|
|---|
| 276 | } else {
|
|---|
| 277 | SaveFlag = true;
|
|---|
| 278 | Log() << Verbose(1) << "Parsing xyz file for new atoms." << endl;
|
|---|
| 279 | if (!mol->AddXYZFile(argv[argptr]))
|
|---|
| 280 | Log() << Verbose(2) << "File not found." << endl;
|
|---|
| 281 | else {
|
|---|
| 282 | Log() << Verbose(2) << "File found and parsed." << endl;
|
|---|
| 283 | configPresent = present;
|
|---|
| 284 | }
|
|---|
| 285 | }
|
|---|
| 286 | break;
|
|---|
| 287 | case 'a':
|
|---|
| 288 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 289 | if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3]))) {
|
|---|
| 290 | ExitFlag = 255;
|
|---|
| 291 | eLog() << Verbose(0) << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
|
|---|
| 292 | performCriticalExit();
|
|---|
| 293 | } else {
|
|---|
| 294 | SaveFlag = true;
|
|---|
| 295 | Log() << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
|
|---|
| 296 | first = new atom;
|
|---|
| 297 | first->type = periode->FindElement(atoi(argv[argptr]));
|
|---|
| 298 | if (first->type != NULL)
|
|---|
| 299 | Log() << Verbose(2) << "found element " << first->type->name << endl;
|
|---|
| 300 | for (int i=NDIM;i--;)
|
|---|
| 301 | first->x.x[i] = atof(argv[argptr+1+i]);
|
|---|
| 302 | if (first->type != NULL) {
|
|---|
| 303 | mol->AddAtom(first); // add to molecule
|
|---|
| 304 | if ((configPresent == empty) && (mol->AtomCount != 0))
|
|---|
| 305 | configPresent = present;
|
|---|
| 306 | } else
|
|---|
| 307 | eLog() << Verbose(1) << "Could not find the specified element." << endl;
|
|---|
| 308 | argptr+=4;
|
|---|
| 309 | }
|
|---|
| 310 | break;
|
|---|
| 311 | default: // no match? Don't step on (this is done in next switch's default)
|
|---|
| 312 | break;
|
|---|
| 313 | }
|
|---|
| 314 | }
|
|---|
| 315 | if (configPresent == present) {
|
|---|
| 316 | switch(argv[argptr-1][1]) {
|
|---|
| 317 | case 'M':
|
|---|
| 318 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
|---|
| 319 | ExitFlag = 255;
|
|---|
| 320 | eLog() << Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl;
|
|---|
| 321 | performCriticalExit();
|
|---|
| 322 | } else {
|
|---|
| 323 | configuration.basis = argv[argptr];
|
|---|
| 324 | Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;
|
|---|
| 325 | argptr+=1;
|
|---|
| 326 | }
|
|---|
| 327 | break;
|
|---|
| 328 | case 'D':
|
|---|
| 329 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 330 | {
|
|---|
| 331 | Log() << Verbose(1) << "Depth-First-Search Analysis." << endl;
|
|---|
| 332 | MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis
|
|---|
| 333 | int *MinimumRingSize = new int[mol->AtomCount];
|
|---|
| 334 | atom ***ListOfLocalAtoms = NULL;
|
|---|
| 335 | class StackClass<bond *> *BackEdgeStack = NULL;
|
|---|
| 336 | class StackClass<bond *> *LocalBackEdgeStack = NULL;
|
|---|
| 337 | mol->CreateAdjacencyList(atof(argv[argptr]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
|
|---|
| 338 | Subgraphs = mol->DepthFirstSearchAnalysis(BackEdgeStack);
|
|---|
| 339 | if (Subgraphs != NULL) {
|
|---|
| 340 | int FragmentCounter = 0;
|
|---|
| 341 | while (Subgraphs->next != NULL) {
|
|---|
| 342 | Subgraphs = Subgraphs->next;
|
|---|
| 343 | Subgraphs->FillBondStructureFromReference(mol, FragmentCounter, ListOfLocalAtoms, false); // we want to keep the created ListOfLocalAtoms
|
|---|
| 344 | LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
|
|---|
| 345 | Subgraphs->Leaf->PickLocalBackEdges(ListOfLocalAtoms[FragmentCounter], BackEdgeStack, LocalBackEdgeStack);
|
|---|
| 346 | Subgraphs->Leaf->CyclicStructureAnalysis(LocalBackEdgeStack, MinimumRingSize);
|
|---|
| 347 | delete(LocalBackEdgeStack);
|
|---|
| 348 | delete(Subgraphs->previous);
|
|---|
| 349 | FragmentCounter++;
|
|---|
| 350 | }
|
|---|
| 351 | delete(Subgraphs);
|
|---|
| 352 | for (int i=0;i<FragmentCounter;i++)
|
|---|
| 353 | Free(&ListOfLocalAtoms[i]);
|
|---|
| 354 | Free(&ListOfLocalAtoms);
|
|---|
| 355 | }
|
|---|
| 356 | delete(BackEdgeStack);
|
|---|
| 357 | delete[](MinimumRingSize);
|
|---|
| 358 | }
|
|---|
| 359 | //argptr+=1;
|
|---|
| 360 | break;
|
|---|
| 361 | case 'C':
|
|---|
| 362 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 363 | if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-') || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-')) {
|
|---|
| 364 | ExitFlag = 255;
|
|---|
| 365 | eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <Z> <output> <bin output>" << endl;
|
|---|
| 366 | performCriticalExit();
|
|---|
| 367 | } else {
|
|---|
| 368 | SaveFlag = false;
|
|---|
| 369 | ofstream output(argv[argptr+1]);
|
|---|
| 370 | ofstream binoutput(argv[argptr+2]);
|
|---|
| 371 | const double radius = 5.;
|
|---|
| 372 |
|
|---|
| 373 | // get the boundary
|
|---|
| 374 | class molecule *Boundary = NULL;
|
|---|
| 375 | class Tesselation *TesselStruct = NULL;
|
|---|
| 376 | const LinkedCell *LCList = NULL;
|
|---|
| 377 | // find biggest molecule
|
|---|
| 378 | int counter = 0;
|
|---|
| 379 | for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
|
|---|
| 380 | if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
|
|---|
| 381 | Boundary = *BigFinder;
|
|---|
| 382 | }
|
|---|
| 383 | counter++;
|
|---|
| 384 | }
|
|---|
| 385 | bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
|
|---|
| 386 | counter = 0;
|
|---|
| 387 | for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
|
|---|
| 388 | Actives[counter] = (*BigFinder)->ActiveFlag;
|
|---|
| 389 | (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
|
|---|
| 390 | }
|
|---|
| 391 | LCList = new LinkedCell(Boundary, 2.*radius);
|
|---|
| 392 | element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
|
|---|
| 393 | FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
|
|---|
| 394 | int ranges[NDIM] = {1,1,1};
|
|---|
| 395 | CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
|
|---|
| 396 | BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
|
|---|
| 397 | OutputCorrelation ( &binoutput, binmap );
|
|---|
| 398 | output.close();
|
|---|
| 399 | binoutput.close();
|
|---|
| 400 | for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
|
|---|
| 401 | (*BigFinder)->ActiveFlag = Actives[counter];
|
|---|
| 402 | Free(&Actives);
|
|---|
| 403 | delete(LCList);
|
|---|
| 404 | delete(TesselStruct);
|
|---|
| 405 | argptr+=3;
|
|---|
| 406 | }
|
|---|
| 407 | break;
|
|---|
| 408 | case 'E':
|
|---|
| 409 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 410 | if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
|
|---|
| 411 | ExitFlag = 255;
|
|---|
| 412 | eLog() << Verbose(0) << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
|
|---|
| 413 | performCriticalExit();
|
|---|
| 414 | } else {
|
|---|
| 415 | SaveFlag = true;
|
|---|
| 416 | Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
|
|---|
| 417 | first = mol->FindAtom(atoi(argv[argptr]));
|
|---|
| 418 | first->type = periode->FindElement(atoi(argv[argptr+1]));
|
|---|
| 419 | argptr+=2;
|
|---|
| 420 | }
|
|---|
| 421 | break;
|
|---|
| 422 | case 'F':
|
|---|
| 423 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 424 | if (argptr+5 >=argc) {
|
|---|
| 425 | ExitFlag = 255;
|
|---|
| 426 | eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <randatom> <randmol> <DoRotate>" << endl;
|
|---|
| 427 | performCriticalExit();
|
|---|
| 428 | } else {
|
|---|
| 429 | SaveFlag = true;
|
|---|
| 430 | Log() << Verbose(1) << "Filling Box with water molecules." << endl;
|
|---|
| 431 | // construct water molecule
|
|---|
| 432 | molecule *filler = new molecule(periode);;
|
|---|
| 433 | molecule *Filling = NULL;
|
|---|
| 434 | atom *second = NULL, *third = NULL;
|
|---|
| 435 | first = new atom();
|
|---|
| 436 | first->type = periode->FindElement(1);
|
|---|
| 437 | first->x.Init(0.441, -0.143, 0.);
|
|---|
| 438 | filler->AddAtom(first);
|
|---|
| 439 | second = new atom();
|
|---|
| 440 | second->type = periode->FindElement(1);
|
|---|
| 441 | second->x.Init(-0.464, 1.137, 0.0);
|
|---|
| 442 | filler->AddAtom(second);
|
|---|
| 443 | third = new atom();
|
|---|
| 444 | third->type = periode->FindElement(8);
|
|---|
| 445 | third->x.Init(-0.464, 0.177, 0.);
|
|---|
| 446 | filler->AddAtom(third);
|
|---|
| 447 | filler->AddBond(first, third, 1);
|
|---|
| 448 | filler->AddBond(second, third, 1);
|
|---|
| 449 | // call routine
|
|---|
| 450 | double distance[NDIM];
|
|---|
| 451 | for (int i=0;i<NDIM;i++)
|
|---|
| 452 | distance[i] = atof(argv[argptr+i]);
|
|---|
| 453 | Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atoi(argv[argptr+5]));
|
|---|
| 454 | if (Filling != NULL) {
|
|---|
| 455 | molecules->insert(Filling);
|
|---|
| 456 | }
|
|---|
| 457 | delete(filler);
|
|---|
| 458 | argptr+=6;
|
|---|
| 459 | }
|
|---|
| 460 | break;
|
|---|
| 461 | case 'A':
|
|---|
| 462 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 463 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
|---|
| 464 | ExitFlag =255;
|
|---|
| 465 | eLog() << Verbose(0) << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
|
|---|
| 466 | performCriticalExit();
|
|---|
| 467 | } else {
|
|---|
| 468 | Log() << Verbose(0) << "Parsing bonds from " << argv[argptr] << "." << endl;
|
|---|
| 469 | ifstream *input = new ifstream(argv[argptr]);
|
|---|
| 470 | mol->CreateAdjacencyListFromDbondFile(input);
|
|---|
| 471 | input->close();
|
|---|
| 472 | argptr+=1;
|
|---|
| 473 | }
|
|---|
| 474 | break;
|
|---|
| 475 | case 'N':
|
|---|
| 476 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 477 | if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
|
|---|
| 478 | ExitFlag = 255;
|
|---|
| 479 | eLog() << Verbose(0) << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
|
|---|
| 480 | performCriticalExit();
|
|---|
| 481 | } else {
|
|---|
| 482 | class Tesselation *T = NULL;
|
|---|
| 483 | const LinkedCell *LCList = NULL;
|
|---|
| 484 | molecule * Boundary = NULL;
|
|---|
| 485 | //string filename(argv[argptr+1]);
|
|---|
| 486 | //filename.append(".csv");
|
|---|
| 487 | Log() << Verbose(0) << "Evaluating non-convex envelope of biggest molecule.";
|
|---|
| 488 | Log() << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
|
|---|
| 489 | // find biggest molecule
|
|---|
| 490 | int counter = 0;
|
|---|
| 491 | for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
|
|---|
| 492 | (*BigFinder)->CountAtoms();
|
|---|
| 493 | if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
|
|---|
| 494 | Boundary = *BigFinder;
|
|---|
| 495 | }
|
|---|
| 496 | counter++;
|
|---|
| 497 | }
|
|---|
| 498 | Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl;
|
|---|
| 499 | start = clock();
|
|---|
| 500 | LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.);
|
|---|
| 501 | FindNonConvexBorder(Boundary, T, LCList, atof(argv[argptr]), argv[argptr+1]);
|
|---|
| 502 | //FindDistributionOfEllipsoids(T, &LCList, N, number, filename.c_str());
|
|---|
| 503 | end = clock();
|
|---|
| 504 | Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
|
|---|
| 505 | delete(LCList);
|
|---|
| 506 | argptr+=2;
|
|---|
| 507 | }
|
|---|
| 508 | break;
|
|---|
| 509 | case 'S':
|
|---|
| 510 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 511 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
|---|
| 512 | ExitFlag = 255;
|
|---|
| 513 | eLog() << Verbose(0) << "Not enough or invalid arguments given for storing tempature: -S <temperature file>" << endl;
|
|---|
| 514 | performCriticalExit();
|
|---|
| 515 | } else {
|
|---|
| 516 | Log() << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
|
|---|
| 517 | ofstream *output = new ofstream(argv[argptr], ios::trunc);
|
|---|
| 518 | if (!mol->OutputTemperatureFromTrajectories(output, 0, mol->MDSteps))
|
|---|
| 519 | Log() << Verbose(2) << "File could not be written." << endl;
|
|---|
| 520 | else
|
|---|
| 521 | Log() << Verbose(2) << "File stored." << endl;
|
|---|
| 522 | output->close();
|
|---|
| 523 | delete(output);
|
|---|
| 524 | argptr+=1;
|
|---|
| 525 | }
|
|---|
| 526 | break;
|
|---|
| 527 | case 'L':
|
|---|
| 528 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 529 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
|---|
| 530 | ExitFlag = 255;
|
|---|
| 531 | eLog() << Verbose(0) << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl;
|
|---|
| 532 | performCriticalExit();
|
|---|
| 533 | } else {
|
|---|
| 534 | SaveFlag = true;
|
|---|
| 535 | Log() << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
|
|---|
| 536 | if (atoi(argv[argptr+3]) == 1)
|
|---|
| 537 | Log() << Verbose(1) << "Using Identity for the permutation map." << endl;
|
|---|
| 538 | if (!mol->LinearInterpolationBetweenConfiguration(atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration, atoi(argv[argptr+3])) == 1 ? true : false)
|
|---|
| 539 | Log() << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
|
|---|
| 540 | else
|
|---|
| 541 | Log() << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
|
|---|
| 542 | argptr+=4;
|
|---|
| 543 | }
|
|---|
| 544 | break;
|
|---|
| 545 | case 'P':
|
|---|
| 546 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 547 | if ((argptr >= argc) || (argv[argptr][0] == '-')) {
|
|---|
| 548 | ExitFlag = 255;
|
|---|
| 549 | eLog() << Verbose(0) << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
|
|---|
| 550 | performCriticalExit();
|
|---|
| 551 | } else {
|
|---|
| 552 | SaveFlag = true;
|
|---|
| 553 | Log() << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
|
|---|
| 554 | if (!mol->VerletForceIntegration(argv[argptr], configuration))
|
|---|
| 555 | Log() << Verbose(2) << "File not found." << endl;
|
|---|
| 556 | else
|
|---|
| 557 | Log() << Verbose(2) << "File found and parsed." << endl;
|
|---|
| 558 | argptr+=1;
|
|---|
| 559 | }
|
|---|
| 560 | break;
|
|---|
| 561 | case 'R':
|
|---|
| 562 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 563 | if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
|
|---|
| 564 | ExitFlag = 255;
|
|---|
| 565 | eLog() << Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl;
|
|---|
| 566 | performCriticalExit();
|
|---|
| 567 | } else {
|
|---|
| 568 | SaveFlag = true;
|
|---|
| 569 | Log() << Verbose(1) << "Removing atoms around " << argv[argptr] << " with radius " << argv[argptr+1] << "." << endl;
|
|---|
| 570 | double tmp1 = atof(argv[argptr+1]);
|
|---|
| 571 | atom *third = mol->FindAtom(atoi(argv[argptr]));
|
|---|
| 572 | atom *first = mol->start;
|
|---|
| 573 | if ((third != NULL) && (first != mol->end)) {
|
|---|
| 574 | atom *second = first->next;
|
|---|
| 575 | while(second != mol->end) {
|
|---|
| 576 | first = second;
|
|---|
| 577 | second = first->next;
|
|---|
| 578 | if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
|
|---|
| 579 | mol->RemoveAtom(first);
|
|---|
| 580 | }
|
|---|
| 581 | } else {
|
|---|
| 582 | eLog() << Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl;
|
|---|
| 583 | }
|
|---|
| 584 | argptr+=2;
|
|---|
| 585 | }
|
|---|
| 586 | break;
|
|---|
| 587 | case 't':
|
|---|
| 588 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 589 | if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
|---|
| 590 | ExitFlag = 255;
|
|---|
| 591 | eLog() << Verbose(0) << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
|
|---|
| 592 | performCriticalExit();
|
|---|
| 593 | } else {
|
|---|
| 594 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 595 | SaveFlag = true;
|
|---|
| 596 | Log() << Verbose(1) << "Translating all ions by given vector." << endl;
|
|---|
| 597 | for (int i=NDIM;i--;)
|
|---|
| 598 | x.x[i] = atof(argv[argptr+i]);
|
|---|
| 599 | mol->Translate((const Vector *)&x);
|
|---|
| 600 | argptr+=3;
|
|---|
| 601 | }
|
|---|
| 602 | break;
|
|---|
| 603 | case 'T':
|
|---|
| 604 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 605 | if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
|---|
| 606 | ExitFlag = 255;
|
|---|
| 607 | eLog() << Verbose(0) << "Not enough or invalid arguments given for periodic translation: -T <x> <y> <z>" << endl;
|
|---|
| 608 | performCriticalExit();
|
|---|
| 609 | } else {
|
|---|
| 610 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 611 | SaveFlag = true;
|
|---|
| 612 | Log() << Verbose(1) << "Translating all ions periodically by given vector." << endl;
|
|---|
| 613 | for (int i=NDIM;i--;)
|
|---|
| 614 | x.x[i] = atof(argv[argptr+i]);
|
|---|
| 615 | mol->TranslatePeriodically((const Vector *)&x);
|
|---|
| 616 | argptr+=3;
|
|---|
| 617 | }
|
|---|
| 618 | break;
|
|---|
| 619 | case 's':
|
|---|
| 620 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 621 | if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
|---|
| 622 | ExitFlag = 255;
|
|---|
| 623 | eLog() << Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl;
|
|---|
| 624 | performCriticalExit();
|
|---|
| 625 | } else {
|
|---|
| 626 | SaveFlag = true;
|
|---|
| 627 | j = -1;
|
|---|
| 628 | Log() << Verbose(1) << "Scaling all ion positions by factor." << endl;
|
|---|
| 629 | factor = new double[NDIM];
|
|---|
| 630 | factor[0] = atof(argv[argptr]);
|
|---|
| 631 | factor[1] = atof(argv[argptr+1]);
|
|---|
| 632 | factor[2] = atof(argv[argptr+2]);
|
|---|
| 633 | mol->Scale((const double ** const)&factor);
|
|---|
| 634 | for (int i=0;i<NDIM;i++) {
|
|---|
| 635 | j += i+1;
|
|---|
| 636 | x.x[i] = atof(argv[NDIM+i]);
|
|---|
| 637 | mol->cell_size[j]*=factor[i];
|
|---|
| 638 | }
|
|---|
| 639 | delete[](factor);
|
|---|
| 640 | argptr+=3;
|
|---|
| 641 | }
|
|---|
| 642 | break;
|
|---|
| 643 | case 'b':
|
|---|
| 644 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 645 | if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
|
|---|
| 646 | ExitFlag = 255;
|
|---|
| 647 | eLog() << Verbose(0) << "Not enough or invalid arguments given for centering in box: -b <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
|
|---|
| 648 | performCriticalExit();
|
|---|
| 649 | } else {
|
|---|
| 650 | SaveFlag = true;
|
|---|
| 651 | j = -1;
|
|---|
| 652 | Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
|
|---|
| 653 | for (int i=0;i<6;i++) {
|
|---|
| 654 | mol->cell_size[i] = atof(argv[argptr+i]);
|
|---|
| 655 | }
|
|---|
| 656 | // center
|
|---|
| 657 | mol->CenterInBox();
|
|---|
| 658 | argptr+=6;
|
|---|
| 659 | }
|
|---|
| 660 | break;
|
|---|
| 661 | case 'B':
|
|---|
| 662 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 663 | if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
|
|---|
| 664 | ExitFlag = 255;
|
|---|
| 665 | eLog() << Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl;
|
|---|
| 666 | performCriticalExit();
|
|---|
| 667 | } else {
|
|---|
| 668 | SaveFlag = true;
|
|---|
| 669 | j = -1;
|
|---|
| 670 | Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
|
|---|
| 671 | for (int i=0;i<6;i++) {
|
|---|
| 672 | mol->cell_size[i] = atof(argv[argptr+i]);
|
|---|
| 673 | }
|
|---|
| 674 | // center
|
|---|
| 675 | mol->BoundInBox();
|
|---|
| 676 | argptr+=6;
|
|---|
| 677 | }
|
|---|
| 678 | break;
|
|---|
| 679 | case 'c':
|
|---|
| 680 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 681 | if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
|---|
| 682 | ExitFlag = 255;
|
|---|
| 683 | eLog() << Verbose(0) << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
|
|---|
| 684 | performCriticalExit();
|
|---|
| 685 | } else {
|
|---|
| 686 | SaveFlag = true;
|
|---|
| 687 | j = -1;
|
|---|
| 688 | Log() << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
|
|---|
| 689 | // make every coordinate positive
|
|---|
| 690 | mol->CenterEdge(&x);
|
|---|
| 691 | // update Box of atoms by boundary
|
|---|
| 692 | mol->SetBoxDimension(&x);
|
|---|
| 693 | // translate each coordinate by boundary
|
|---|
| 694 | j=-1;
|
|---|
| 695 | for (int i=0;i<NDIM;i++) {
|
|---|
| 696 | j += i+1;
|
|---|
| 697 | x.x[i] = atof(argv[argptr+i]);
|
|---|
| 698 | mol->cell_size[j] += x.x[i]*2.;
|
|---|
| 699 | }
|
|---|
| 700 | mol->Translate((const Vector *)&x);
|
|---|
| 701 | argptr+=3;
|
|---|
| 702 | }
|
|---|
| 703 | break;
|
|---|
| 704 | case 'O':
|
|---|
| 705 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 706 | SaveFlag = true;
|
|---|
| 707 | Log() << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl;
|
|---|
| 708 | x.Zero();
|
|---|
| 709 | mol->CenterEdge(&x);
|
|---|
| 710 | mol->SetBoxDimension(&x);
|
|---|
| 711 | argptr+=0;
|
|---|
| 712 | break;
|
|---|
| 713 | case 'r':
|
|---|
| 714 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 715 | if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr]))) {
|
|---|
| 716 | ExitFlag = 255;
|
|---|
| 717 | eLog() << Verbose(0) << "Not enough or invalid arguments given for removing atoms: -r <id>" << endl;
|
|---|
| 718 | performCriticalExit();
|
|---|
| 719 | } else {
|
|---|
| 720 | SaveFlag = true;
|
|---|
| 721 | Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl;
|
|---|
| 722 | atom *first = mol->FindAtom(atoi(argv[argptr]));
|
|---|
| 723 | mol->RemoveAtom(first);
|
|---|
| 724 | argptr+=1;
|
|---|
| 725 | }
|
|---|
| 726 | break;
|
|---|
| 727 | case 'f':
|
|---|
| 728 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 729 | if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
|
|---|
| 730 | ExitFlag = 255;
|
|---|
| 731 | eLog() << Verbose(0) << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
|
|---|
| 732 | performCriticalExit();
|
|---|
| 733 | } else {
|
|---|
| 734 | Log() << Verbose(0) << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
|
|---|
| 735 | Log() << Verbose(0) << "Creating connection matrix..." << endl;
|
|---|
| 736 | start = clock();
|
|---|
| 737 | mol->CreateAdjacencyList(atof(argv[argptr++]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
|
|---|
| 738 | Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
|
|---|
| 739 | if (mol->first->next != mol->last) {
|
|---|
| 740 | ExitFlag = mol->FragmentMolecule(atoi(argv[argptr]), &configuration);
|
|---|
| 741 | }
|
|---|
| 742 | end = clock();
|
|---|
| 743 | Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
|
|---|
| 744 | argptr+=2;
|
|---|
| 745 | }
|
|---|
| 746 | break;
|
|---|
| 747 | case 'm':
|
|---|
| 748 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 749 | j = atoi(argv[argptr++]);
|
|---|
| 750 | if ((j<0) || (j>1)) {
|
|---|
| 751 | eLog() << Verbose(1) << "Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
|
|---|
| 752 | j = 0;
|
|---|
| 753 | }
|
|---|
| 754 | if (j) {
|
|---|
| 755 | SaveFlag = true;
|
|---|
| 756 | Log() << Verbose(0) << "Converting to prinicipal axis system." << endl;
|
|---|
| 757 | } else
|
|---|
| 758 | Log() << Verbose(0) << "Evaluating prinicipal axis." << endl;
|
|---|
| 759 | mol->PrincipalAxisSystem((bool)j);
|
|---|
| 760 | break;
|
|---|
| 761 | case 'o':
|
|---|
| 762 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 763 | if ((argptr+1 >= argc) || (argv[argptr][0] == '-')){
|
|---|
| 764 | ExitFlag = 255;
|
|---|
| 765 | eLog() << Verbose(0) << "Not enough or invalid arguments given for convex envelope: -o <convex output file> <non-convex output file>" << endl;
|
|---|
| 766 | performCriticalExit();
|
|---|
| 767 | } else {
|
|---|
| 768 | class Tesselation *TesselStruct = NULL;
|
|---|
| 769 | const LinkedCell *LCList = NULL;
|
|---|
| 770 | Log() << Verbose(0) << "Evaluating volume of the convex envelope.";
|
|---|
| 771 | Log() << Verbose(1) << "Storing tecplot convex data in " << argv[argptr] << "." << endl;
|
|---|
| 772 | Log() << Verbose(1) << "Storing tecplot non-convex data in " << argv[argptr+1] << "." << endl;
|
|---|
| 773 | LCList = new LinkedCell(mol, 10.);
|
|---|
| 774 | //FindConvexBorder(mol, LCList, argv[argptr]);
|
|---|
| 775 | FindNonConvexBorder(mol, TesselStruct, LCList, 5., argv[argptr+1]);
|
|---|
| 776 | // RemoveAllBoundaryPoints(TesselStruct, mol, argv[argptr]);
|
|---|
| 777 | double volumedifference = ConvexizeNonconvexEnvelope(TesselStruct, mol, argv[argptr]);
|
|---|
| 778 | double clustervolume = VolumeOfConvexEnvelope(TesselStruct, &configuration);
|
|---|
| 779 | Log() << Verbose(0) << "The tesselated volume area is " << clustervolume << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
|
|---|
| 780 | Log() << Verbose(0) << "The non-convex tesselated volume area is " << clustervolume-volumedifference << " " << (configuration.GetIsAngstroem() ? "angstrom" : "atomiclength") << "^3." << endl;
|
|---|
| 781 | delete(TesselStruct);
|
|---|
| 782 | delete(LCList);
|
|---|
| 783 | argptr+=2;
|
|---|
| 784 | }
|
|---|
| 785 | break;
|
|---|
| 786 | case 'U':
|
|---|
| 787 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 788 | if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
|
|---|
| 789 | ExitFlag = 255;
|
|---|
| 790 | eLog() << Verbose(0) << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
|
|---|
| 791 | performCriticalExit();
|
|---|
| 792 | } else {
|
|---|
| 793 | volume = atof(argv[argptr++]);
|
|---|
| 794 | Log() << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
|
|---|
| 795 | }
|
|---|
| 796 | case 'u':
|
|---|
| 797 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 798 | if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) ) {
|
|---|
| 799 | if (volume != -1)
|
|---|
| 800 | ExitFlag = 255;
|
|---|
| 801 | eLog() << Verbose(0) << "Not enough arguments given for suspension: -u <density>" << endl;
|
|---|
| 802 | performCriticalExit();
|
|---|
| 803 | } else {
|
|---|
| 804 | double density;
|
|---|
| 805 | SaveFlag = true;
|
|---|
| 806 | Log() << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
|
|---|
| 807 | density = atof(argv[argptr++]);
|
|---|
| 808 | if (density < 1.0) {
|
|---|
| 809 | eLog() << Verbose(1) << "Density must be greater than 1.0g/cm^3 !" << endl;
|
|---|
| 810 | density = 1.3;
|
|---|
| 811 | }
|
|---|
| 812 | // for(int i=0;i<NDIM;i++) {
|
|---|
| 813 | // repetition[i] = atoi(argv[argptr++]);
|
|---|
| 814 | // if (repetition[i] < 1)
|
|---|
| 815 | // eLog() << Verbose(1) << "repetition value must be greater 1!" << endl;
|
|---|
| 816 | // repetition[i] = 1;
|
|---|
| 817 | // }
|
|---|
| 818 | PrepareClustersinWater(&configuration, mol, volume, density); // if volume == 0, will calculate from ConvexEnvelope
|
|---|
| 819 | }
|
|---|
| 820 | break;
|
|---|
| 821 | case 'd':
|
|---|
| 822 | if (ExitFlag == 0) ExitFlag = 1;
|
|---|
| 823 | if ((argptr+2 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
|
|---|
| 824 | ExitFlag = 255;
|
|---|
| 825 | eLog() << Verbose(0) << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
|
|---|
| 826 | performCriticalExit();
|
|---|
| 827 | } else {
|
|---|
| 828 | SaveFlag = true;
|
|---|
| 829 | for (int axis = 1; axis <= NDIM; axis++) {
|
|---|
| 830 | int faktor = atoi(argv[argptr++]);
|
|---|
| 831 | int count;
|
|---|
| 832 | element ** Elements;
|
|---|
| 833 | Vector ** vectors;
|
|---|
| 834 | if (faktor < 1) {
|
|---|
| 835 | eLog() << Verbose(1) << "Repetition factor mus be greater than 1!" << endl;
|
|---|
| 836 | faktor = 1;
|
|---|
| 837 | }
|
|---|
| 838 | mol->CountAtoms(); // recount atoms
|
|---|
| 839 | if (mol->AtomCount != 0) { // if there is more than none
|
|---|
| 840 | count = mol->AtomCount; // is changed becausing of adding, thus has to be stored away beforehand
|
|---|
| 841 | Elements = new element *[count];
|
|---|
| 842 | vectors = new Vector *[count];
|
|---|
| 843 | j = 0;
|
|---|
| 844 | first = mol->start;
|
|---|
| 845 | while (first->next != mol->end) { // make a list of all atoms with coordinates and element
|
|---|
| 846 | first = first->next;
|
|---|
| 847 | Elements[j] = first->type;
|
|---|
| 848 | vectors[j] = &first->x;
|
|---|
| 849 | j++;
|
|---|
| 850 | }
|
|---|
| 851 | if (count != j)
|
|---|
| 852 | eLog() << Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
|
|---|
| 853 | x.Zero();
|
|---|
| 854 | y.Zero();
|
|---|
| 855 | y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
|
|---|
| 856 | for (int i=1;i<faktor;i++) { // then add this list with respective translation factor times
|
|---|
| 857 | x.AddVector(&y); // per factor one cell width further
|
|---|
| 858 | for (int k=count;k--;) { // go through every atom of the original cell
|
|---|
| 859 | first = new atom(); // create a new body
|
|---|
| 860 | first->x.CopyVector(vectors[k]); // use coordinate of original atom
|
|---|
| 861 | first->x.AddVector(&x); // translate the coordinates
|
|---|
| 862 | first->type = Elements[k]; // insert original element
|
|---|
| 863 | mol->AddAtom(first); // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
|
|---|
| 864 | }
|
|---|
| 865 | }
|
|---|
| 866 | // free memory
|
|---|
| 867 | delete[](Elements);
|
|---|
| 868 | delete[](vectors);
|
|---|
| 869 | // correct cell size
|
|---|
| 870 | if (axis < 0) { // if sign was negative, we have to translate everything
|
|---|
| 871 | x.Zero();
|
|---|
| 872 | x.AddVector(&y);
|
|---|
| 873 | x.Scale(-(faktor-1));
|
|---|
| 874 | mol->Translate(&x);
|
|---|
| 875 | }
|
|---|
| 876 | mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
|
|---|
| 877 | }
|
|---|
| 878 | }
|
|---|
| 879 | }
|
|---|
| 880 | break;
|
|---|
| 881 | default: // no match? Step on
|
|---|
| 882 | if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
|
|---|
| 883 | argptr++;
|
|---|
| 884 | break;
|
|---|
| 885 | }
|
|---|
| 886 | }
|
|---|
| 887 | } else argptr++;
|
|---|
| 888 | } while (argptr < argc);
|
|---|
| 889 | if (SaveFlag)
|
|---|
| 890 | configuration.SaveAll(ConfigFileName, periode, molecules);
|
|---|
| 891 | } else { // no arguments, hence scan the elements db
|
|---|
| 892 | if (periode->LoadPeriodentafel(configuration.databasepath))
|
|---|
| 893 | Log() << Verbose(0) << "Element list loaded successfully." << endl;
|
|---|
| 894 | else
|
|---|
| 895 | Log() << Verbose(0) << "Element list loading failed." << endl;
|
|---|
| 896 | configuration.RetrieveConfigPathAndName("main_pcp_linux");
|
|---|
| 897 | }
|
|---|
| 898 | return(ExitFlag);
|
|---|
| 899 | };
|
|---|
| 900 |
|
|---|
| 901 | /***************************************** Functions used to build all menus **********************/
|
|---|
| 902 |
|
|---|
| 903 | void populateEditMoleculesMenu(Menu* editMoleculesMenu,MoleculeListClass *molecules, config *configuration, periodentafel *periode){
|
|---|
| 904 | // build the EditMoleculesMenu
|
|---|
| 905 | Action *createMoleculeAction = new MethodAction("createMoleculeAction",boost::bind(&MoleculeListClass::createNewMolecule,molecules,periode));
|
|---|
| 906 | new ActionMenuItem('c',"create new molecule",editMoleculesMenu,createMoleculeAction);
|
|---|
| 907 |
|
|---|
| 908 | Action *loadMoleculeAction = new MethodAction("loadMoleculeAction",boost::bind(&MoleculeListClass::loadFromXYZ,molecules,periode));
|
|---|
| 909 | new ActionMenuItem('l',"load molecule from xyz file",editMoleculesMenu,loadMoleculeAction);
|
|---|
| 910 |
|
|---|
| 911 | Action *changeFilenameAction = new ChangeMoleculeNameAction(molecules);
|
|---|
| 912 | new ActionMenuItem('n',"change molecule's name",editMoleculesMenu,changeFilenameAction);
|
|---|
| 913 |
|
|---|
| 914 | Action *giveFilenameAction = new MethodAction("giveFilenameAction",boost::bind(&MoleculeListClass::setMoleculeFilename,molecules));
|
|---|
| 915 | new ActionMenuItem('N',"give molecules filename",editMoleculesMenu,giveFilenameAction);
|
|---|
| 916 |
|
|---|
| 917 | Action *parseAtomsAction = new MethodAction("parseAtomsAction",boost::bind(&MoleculeListClass::parseXYZIntoMolecule,molecules));
|
|---|
| 918 | new ActionMenuItem('p',"parse atoms in xyz file into molecule",editMoleculesMenu,parseAtomsAction);
|
|---|
| 919 |
|
|---|
| 920 | Action *eraseMoleculeAction = new MethodAction("eraseMoleculeAction",boost::bind(&MoleculeListClass::eraseMolecule,molecules));
|
|---|
| 921 | new ActionMenuItem('r',"remove a molecule",editMoleculesMenu,eraseMoleculeAction);
|
|---|
| 922 | }
|
|---|
| 923 |
|
|---|
| 924 |
|
|---|
| 925 | /********************************************** Main routine **************************************/
|
|---|
| 926 |
|
|---|
| 927 | int main(int argc, char **argv)
|
|---|
| 928 | {
|
|---|
| 929 | periodentafel *periode = new periodentafel;
|
|---|
| 930 | MoleculeListClass *molecules = new MoleculeListClass;
|
|---|
| 931 | molecule *mol = NULL;
|
|---|
| 932 | config *configuration = new config;
|
|---|
| 933 | Vector x, y, z, n;
|
|---|
| 934 | ifstream test;
|
|---|
| 935 | ofstream output;
|
|---|
| 936 | string line;
|
|---|
| 937 | char *ConfigFileName = NULL;
|
|---|
| 938 | int j;
|
|---|
| 939 | setVerbosity(0);
|
|---|
| 940 | /* structure of ParseCommandLineOptions will be refactored later */
|
|---|
| 941 | j = ParseCommandLineOptions(argc, argv, molecules, periode, *configuration, ConfigFileName);
|
|---|
| 942 | switch (j){
|
|---|
| 943 | case 255:
|
|---|
| 944 | case 2:
|
|---|
| 945 | case 1:
|
|---|
| 946 | delete (molecules);
|
|---|
| 947 | delete (periode);
|
|---|
| 948 | delete (configuration);
|
|---|
| 949 | Log() << Verbose(0) << "Maximum of allocated memory: " << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
|
|---|
| 950 | Log() << Verbose(0) << "Remaining non-freed memory: " << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
|
|---|
| 951 | MemoryUsageObserver::getInstance()->purgeInstance();
|
|---|
| 952 | logger::purgeInstance();
|
|---|
| 953 | errorLogger::purgeInstance();
|
|---|
| 954 | return (j == 1 ? 0 : j);
|
|---|
| 955 | default:
|
|---|
| 956 | break;
|
|---|
| 957 | }
|
|---|
| 958 | if(molecules->ListOfMolecules.size() == 0){
|
|---|
| 959 | mol = new molecule(periode);
|
|---|
| 960 | if(mol->cell_size[0] == 0.){
|
|---|
| 961 | Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
|
|---|
| 962 | for(int i = 0;i < 6;i++){
|
|---|
| 963 | Log() << Verbose(1) << "Cell size" << i << ": ";
|
|---|
| 964 | cin >> mol->cell_size[i];
|
|---|
| 965 | }
|
|---|
| 966 | }
|
|---|
| 967 |
|
|---|
| 968 | mol->ActiveFlag = true;
|
|---|
| 969 | molecules->insert(mol);
|
|---|
| 970 | }
|
|---|
| 971 |
|
|---|
| 972 | {
|
|---|
| 973 | menuPopulaters populaters;
|
|---|
| 974 | populaters.MakeEditMoleculesMenu = populateEditMoleculesMenu;
|
|---|
| 975 |
|
|---|
| 976 | #ifdef USE_GUI_QT
|
|---|
| 977 | UIFactory::makeUserInterface(UIFactory::QT4);
|
|---|
| 978 | #else
|
|---|
| 979 | UIFactory::makeUserInterface(UIFactory::Text);
|
|---|
| 980 | #endif
|
|---|
| 981 | MainWindow *mainWindow = UIFactory::get()->makeMainWindow(populaters,molecules, configuration, periode, ConfigFileName);
|
|---|
| 982 | mainWindow->display();
|
|---|
| 983 | delete mainWindow;
|
|---|
| 984 | }
|
|---|
| 985 |
|
|---|
| 986 | if(periode->StorePeriodentafel(configuration->databasepath))
|
|---|
| 987 | Log() << Verbose(0) << "Saving of elements.db successful." << endl;
|
|---|
| 988 |
|
|---|
| 989 | else
|
|---|
| 990 | Log() << Verbose(0) << "Saving of elements.db failed." << endl;
|
|---|
| 991 |
|
|---|
| 992 | delete (molecules);
|
|---|
| 993 | delete(periode);
|
|---|
| 994 | delete(configuration);
|
|---|
| 995 |
|
|---|
| 996 |
|
|---|
| 997 |
|
|---|
| 998 | Log() << Verbose(0) << "Maximum of allocated memory: "
|
|---|
| 999 | << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
|
|---|
| 1000 | Log() << Verbose(0) << "Remaining non-freed memory: "
|
|---|
| 1001 | << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
|
|---|
| 1002 | MemoryUsageObserver::purgeInstance();
|
|---|
| 1003 | logger::purgeInstance();
|
|---|
| 1004 | errorLogger::purgeInstance();
|
|---|
| 1005 | UIFactory::purgeInstance();
|
|---|
| 1006 | ActionRegistry::purgeRegistry();
|
|---|
| 1007 | return (0);
|
|---|
| 1008 | }
|
|---|
| 1009 |
|
|---|
| 1010 | /********************************************** E N D **************************************************/
|
|---|
