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