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