| 1 | /*
|
|---|
| 2 | * Project: MoleCuilder
|
|---|
| 3 | * Description: creates and alters molecular systems
|
|---|
| 4 | * Copyright (C) 2010 University of Bonn. All rights reserved.
|
|---|
| 5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
|---|
| 6 | */
|
|---|
| 7 |
|
|---|
| 8 | /** \file config.cpp
|
|---|
| 9 | *
|
|---|
| 10 | * Function implementations for the class config.
|
|---|
| 11 | *
|
|---|
| 12 | */
|
|---|
| 13 |
|
|---|
| 14 | // include config.h
|
|---|
| 15 | #ifdef HAVE_CONFIG_H
|
|---|
| 16 | #include <config.h>
|
|---|
| 17 | #endif
|
|---|
| 18 |
|
|---|
| 19 | #include "Helpers/MemDebug.hpp"
|
|---|
| 20 |
|
|---|
| 21 | #include <stdio.h>
|
|---|
| 22 | #include <cstring>
|
|---|
| 23 |
|
|---|
| 24 | //#include "Actions/FragmentationAction/SubgraphDissectionAction.hpp"
|
|---|
| 25 | #include "atom.hpp"
|
|---|
| 26 | #include "bond.hpp"
|
|---|
| 27 | #include "bondgraph.hpp"
|
|---|
| 28 | #include "config.hpp"
|
|---|
| 29 | #include "ConfigFileBuffer.hpp"
|
|---|
| 30 | #include "element.hpp"
|
|---|
| 31 | #include "Helpers/helpers.hpp"
|
|---|
| 32 | #include "Helpers/Info.hpp"
|
|---|
| 33 | #include "lists.hpp"
|
|---|
| 34 | #include "Helpers/Verbose.hpp"
|
|---|
| 35 | #include "Helpers/Log.hpp"
|
|---|
| 36 | #include "molecule.hpp"
|
|---|
| 37 | #include "molecule.hpp"
|
|---|
| 38 | #include "periodentafel.hpp"
|
|---|
| 39 | #include "ThermoStatContainer.hpp"
|
|---|
| 40 | #include "World.hpp"
|
|---|
| 41 | #include "LinearAlgebra/Matrix.hpp"
|
|---|
| 42 | #include "Box.hpp"
|
|---|
| 43 |
|
|---|
| 44 | /************************************* Functions for class config ***************************/
|
|---|
| 45 |
|
|---|
| 46 | /** Constructor for config file class.
|
|---|
| 47 | */
|
|---|
| 48 | config::config() :
|
|---|
| 49 | BG(NULL),
|
|---|
| 50 | Thermostats(0),
|
|---|
| 51 | PsiType(0),
|
|---|
| 52 | MaxPsiDouble(0),
|
|---|
| 53 | PsiMaxNoUp(0),
|
|---|
| 54 | PsiMaxNoDown(0),
|
|---|
| 55 | MaxMinStopStep(1),
|
|---|
| 56 | InitMaxMinStopStep(1),
|
|---|
| 57 | ProcPEGamma(8),
|
|---|
| 58 | ProcPEPsi(1),
|
|---|
| 59 | configname(NULL),
|
|---|
| 60 | FastParsing(false),
|
|---|
| 61 | Deltat(0.01),
|
|---|
| 62 | basis(""),
|
|---|
| 63 | databasepath(NULL),
|
|---|
| 64 | DoConstrainedMD(0),
|
|---|
| 65 | MaxOuterStep(0),
|
|---|
| 66 | mainname(NULL),
|
|---|
| 67 | defaultpath(NULL),
|
|---|
| 68 | pseudopotpath(NULL),
|
|---|
| 69 | DoOutVis(0),
|
|---|
| 70 | DoOutMes(1),
|
|---|
| 71 | DoOutNICS(0),
|
|---|
| 72 | DoOutOrbitals(0),
|
|---|
| 73 | DoOutCurrent(0),
|
|---|
| 74 | DoFullCurrent(0),
|
|---|
| 75 | DoPerturbation(0),
|
|---|
| 76 | DoWannier(0),
|
|---|
| 77 | CommonWannier(0),
|
|---|
| 78 | SawtoothStart(0.01),
|
|---|
| 79 | VectorPlane(0),
|
|---|
| 80 | VectorCut(0.),
|
|---|
| 81 | UseAddGramSch(1),
|
|---|
| 82 | Seed(1),
|
|---|
| 83 | OutVisStep(10),
|
|---|
| 84 | OutSrcStep(5),
|
|---|
| 85 | MaxPsiStep(0),
|
|---|
| 86 | EpsWannier(1e-7),
|
|---|
| 87 | MaxMinStep(100),
|
|---|
| 88 | RelEpsTotalEnergy(1e-7),
|
|---|
| 89 | RelEpsKineticEnergy(1e-5),
|
|---|
| 90 | MaxMinGapStopStep(0),
|
|---|
| 91 | MaxInitMinStep(100),
|
|---|
| 92 | InitRelEpsTotalEnergy(1e-5),
|
|---|
| 93 | InitRelEpsKineticEnergy(1e-4),
|
|---|
| 94 | InitMaxMinGapStopStep(0),
|
|---|
| 95 | ECut(128.),
|
|---|
| 96 | MaxLevel(5),
|
|---|
| 97 | RiemannTensor(0),
|
|---|
| 98 | LevRFactor(0),
|
|---|
| 99 | RiemannLevel(0),
|
|---|
| 100 | Lev0Factor(2),
|
|---|
| 101 | RTActualUse(0),
|
|---|
| 102 | AddPsis(0),
|
|---|
| 103 | RCut(20.),
|
|---|
| 104 | StructOpt(0),
|
|---|
| 105 | IsAngstroem(1),
|
|---|
| 106 | RelativeCoord(0),
|
|---|
| 107 | MaxTypes(0)
|
|---|
| 108 | {
|
|---|
| 109 | mainname = new char[MAXSTRINGSIZE];
|
|---|
| 110 | defaultpath = new char[MAXSTRINGSIZE];
|
|---|
| 111 | pseudopotpath = new char[MAXSTRINGSIZE];
|
|---|
| 112 | databasepath = new char[MAXSTRINGSIZE];
|
|---|
| 113 | configname = new char[MAXSTRINGSIZE];
|
|---|
| 114 | Thermostats = new ThermoStatContainer();
|
|---|
| 115 | strcpy(mainname,"pcp");
|
|---|
| 116 | strcpy(defaultpath,"not specified");
|
|---|
| 117 | strcpy(pseudopotpath,"not specified");
|
|---|
| 118 | configname[0]='\0';
|
|---|
| 119 | basis = "3-21G";
|
|---|
| 120 | };
|
|---|
| 121 |
|
|---|
| 122 | /** Destructor for config file class.
|
|---|
| 123 | */
|
|---|
| 124 | config::~config()
|
|---|
| 125 | {
|
|---|
| 126 | delete[](mainname);
|
|---|
| 127 | delete[](defaultpath);
|
|---|
| 128 | delete[](pseudopotpath);
|
|---|
| 129 | delete[](databasepath);
|
|---|
| 130 | delete[](configname);
|
|---|
| 131 | if (Thermostats != NULL)
|
|---|
| 132 | delete(Thermostats);
|
|---|
| 133 |
|
|---|
| 134 | if (BG != NULL)
|
|---|
| 135 | delete(BG);
|
|---|
| 136 | };
|
|---|
| 137 |
|
|---|
| 138 | /** Displays menu for editing each entry of the config file.
|
|---|
| 139 | * Nothing fancy here, just lots of Log() << Verbose(0)s for the menu and a switch/case
|
|---|
| 140 | * for each entry of the config file structure.
|
|---|
| 141 | */
|
|---|
| 142 | void config::Edit()
|
|---|
| 143 | {
|
|---|
| 144 | char choice;
|
|---|
| 145 |
|
|---|
| 146 | do {
|
|---|
| 147 | DoLog(0) && (Log() << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl);
|
|---|
| 148 | DoLog(0) && (Log() << Verbose(0) << " A - mainname (prefix for all runtime files)" << endl);
|
|---|
| 149 | DoLog(0) && (Log() << Verbose(0) << " B - Default path (for runtime files)" << endl);
|
|---|
| 150 | DoLog(0) && (Log() << Verbose(0) << " C - Path of pseudopotential files" << endl);
|
|---|
| 151 | DoLog(0) && (Log() << Verbose(0) << " D - Number of coefficient sharing processes" << endl);
|
|---|
| 152 | DoLog(0) && (Log() << Verbose(0) << " E - Number of wave function sharing processes" << endl);
|
|---|
| 153 | DoLog(0) && (Log() << Verbose(0) << " F - 0: Don't output density for OpenDX, 1: do" << endl);
|
|---|
| 154 | DoLog(0) && (Log() << Verbose(0) << " G - 0: Don't output physical data, 1: do" << endl);
|
|---|
| 155 | DoLog(0) && (Log() << Verbose(0) << " H - 0: Don't output densities of each unperturbed orbital for OpenDX, 1: do" << endl);
|
|---|
| 156 | DoLog(0) && (Log() << Verbose(0) << " I - 0: Don't output current density for OpenDX, 1: do" << endl);
|
|---|
| 157 | DoLog(0) && (Log() << Verbose(0) << " J - 0: Don't do the full current calculation, 1: do" << endl);
|
|---|
| 158 | DoLog(0) && (Log() << Verbose(0) << " K - 0: Don't do perturbation calculation to obtain susceptibility and shielding, 1: do" << endl);
|
|---|
| 159 | DoLog(0) && (Log() << Verbose(0) << " L - 0: Wannier centres as calculated, 1: common centre for all, 2: unite centres according to spread, 3: cell centre, 4: shifted to nearest grid point" << endl);
|
|---|
| 160 | DoLog(0) && (Log() << Verbose(0) << " M - Absolute begin of unphysical sawtooth transfer for position operator within cell" << endl);
|
|---|
| 161 | DoLog(0) && (Log() << Verbose(0) << " N - (0,1,2) x,y,z-plane to do two-dimensional current vector cut" << endl);
|
|---|
| 162 | DoLog(0) && (Log() << Verbose(0) << " O - Absolute position along vector cut axis for cut plane" << endl);
|
|---|
| 163 | DoLog(0) && (Log() << Verbose(0) << " P - Additional Gram-Schmidt-Orthonormalization to stabilize numerics" << endl);
|
|---|
| 164 | DoLog(0) && (Log() << Verbose(0) << " Q - Initial integer value of random number generator" << endl);
|
|---|
| 165 | DoLog(0) && (Log() << Verbose(0) << " R - for perturbation 0, for structure optimization defines upper limit of iterations" << endl);
|
|---|
| 166 | DoLog(0) && (Log() << Verbose(0) << " T - Output visual after ...th step" << endl);
|
|---|
| 167 | DoLog(0) && (Log() << Verbose(0) << " U - Output source densities of wave functions after ...th step" << endl);
|
|---|
| 168 | DoLog(0) && (Log() << Verbose(0) << " X - minimization iterations per wave function, if unsure leave at default value 0" << endl);
|
|---|
| 169 | DoLog(0) && (Log() << Verbose(0) << " Y - tolerance value for total spread in iterative Jacobi diagonalization" << endl);
|
|---|
| 170 | DoLog(0) && (Log() << Verbose(0) << " Z - Maximum number of minimization iterations" << endl);
|
|---|
| 171 | DoLog(0) && (Log() << Verbose(0) << " a - Relative change in total energy to stop min. iteration" << endl);
|
|---|
| 172 | DoLog(0) && (Log() << Verbose(0) << " b - Relative change in kinetic energy to stop min. iteration" << endl);
|
|---|
| 173 | DoLog(0) && (Log() << Verbose(0) << " c - Check stop conditions every ..th step during min. iteration" << endl);
|
|---|
| 174 | DoLog(0) && (Log() << Verbose(0) << " e - Maximum number of minimization iterations during initial level" << endl);
|
|---|
| 175 | DoLog(0) && (Log() << Verbose(0) << " f - Relative change in total energy to stop min. iteration during initial level" << endl);
|
|---|
| 176 | DoLog(0) && (Log() << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl);
|
|---|
| 177 | DoLog(0) && (Log() << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl);
|
|---|
| 178 | // Log() << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
|
|---|
| 179 | DoLog(0) && (Log() << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl);
|
|---|
| 180 | DoLog(0) && (Log() << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl);
|
|---|
| 181 | DoLog(0) && (Log() << Verbose(0) << " m - Factor by which grid nodes increase between standard and upper level" << endl);
|
|---|
| 182 | DoLog(0) && (Log() << Verbose(0) << " n - 0: Don't use RiemannTensor, 1: Do" << endl);
|
|---|
| 183 | DoLog(0) && (Log() << Verbose(0) << " o - Factor by which grid nodes increase between Riemann and standard(?) level" << endl);
|
|---|
| 184 | DoLog(0) && (Log() << Verbose(0) << " p - Number of Riemann levels" << endl);
|
|---|
| 185 | DoLog(0) && (Log() << Verbose(0) << " r - 0: Don't Use RiemannTensor, 1: Do" << endl);
|
|---|
| 186 | DoLog(0) && (Log() << Verbose(0) << " s - 0: Doubly occupied orbitals, 1: Up-/Down-Orbitals" << endl);
|
|---|
| 187 | DoLog(0) && (Log() << Verbose(0) << " t - Number of orbitals (depends pn SpinType)" << endl);
|
|---|
| 188 | DoLog(0) && (Log() << Verbose(0) << " u - Number of SpinUp orbitals (depends on SpinType)" << endl);
|
|---|
| 189 | DoLog(0) && (Log() << Verbose(0) << " v - Number of SpinDown orbitals (depends on SpinType)" << endl);
|
|---|
| 190 | DoLog(0) && (Log() << Verbose(0) << " w - Number of additional, unoccupied orbitals" << endl);
|
|---|
| 191 | DoLog(0) && (Log() << Verbose(0) << " x - radial cutoff for ewald summation in Bohrradii" << endl);
|
|---|
| 192 | DoLog(0) && (Log() << Verbose(0) << " y - 0: Don't do structure optimization beforehand, 1: Do" << endl);
|
|---|
| 193 | DoLog(0) && (Log() << Verbose(0) << " z - 0: Units are in Bohr radii, 1: units are in Aengstrom" << endl);
|
|---|
| 194 | DoLog(0) && (Log() << Verbose(0) << " i - 0: Coordinates given in file are absolute, 1: ... are relative to unit cell" << endl);
|
|---|
| 195 | DoLog(0) && (Log() << Verbose(0) << "=========================================================" << endl);
|
|---|
| 196 | DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
|
|---|
| 197 | cin >> choice;
|
|---|
| 198 |
|
|---|
| 199 | switch (choice) {
|
|---|
| 200 | case 'A': // mainname
|
|---|
| 201 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::mainname << "\t new: ");
|
|---|
| 202 | cin >> config::mainname;
|
|---|
| 203 | break;
|
|---|
| 204 | case 'B': // defaultpath
|
|---|
| 205 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::defaultpath << "\t new: ");
|
|---|
| 206 | cin >> config::defaultpath;
|
|---|
| 207 | break;
|
|---|
| 208 | case 'C': // pseudopotpath
|
|---|
| 209 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::pseudopotpath << "\t new: ");
|
|---|
| 210 | cin >> config::pseudopotpath;
|
|---|
| 211 | break;
|
|---|
| 212 |
|
|---|
| 213 | case 'D': // ProcPEGamma
|
|---|
| 214 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::ProcPEGamma << "\t new: ");
|
|---|
| 215 | cin >> config::ProcPEGamma;
|
|---|
| 216 | break;
|
|---|
| 217 | case 'E': // ProcPEPsi
|
|---|
| 218 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::ProcPEPsi << "\t new: ");
|
|---|
| 219 | cin >> config::ProcPEPsi;
|
|---|
| 220 | break;
|
|---|
| 221 | case 'F': // DoOutVis
|
|---|
| 222 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoOutVis << "\t new: ");
|
|---|
| 223 | cin >> config::DoOutVis;
|
|---|
| 224 | break;
|
|---|
| 225 | case 'G': // DoOutMes
|
|---|
| 226 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoOutMes << "\t new: ");
|
|---|
| 227 | cin >> config::DoOutMes;
|
|---|
| 228 | break;
|
|---|
| 229 | case 'H': // DoOutOrbitals
|
|---|
| 230 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoOutOrbitals << "\t new: ");
|
|---|
| 231 | cin >> config::DoOutOrbitals;
|
|---|
| 232 | break;
|
|---|
| 233 | case 'I': // DoOutCurrent
|
|---|
| 234 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoOutCurrent << "\t new: ");
|
|---|
| 235 | cin >> config::DoOutCurrent;
|
|---|
| 236 | break;
|
|---|
| 237 | case 'J': // DoFullCurrent
|
|---|
| 238 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoFullCurrent << "\t new: ");
|
|---|
| 239 | cin >> config::DoFullCurrent;
|
|---|
| 240 | break;
|
|---|
| 241 | case 'K': // DoPerturbation
|
|---|
| 242 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoPerturbation << "\t new: ");
|
|---|
| 243 | cin >> config::DoPerturbation;
|
|---|
| 244 | break;
|
|---|
| 245 | case 'L': // CommonWannier
|
|---|
| 246 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::CommonWannier << "\t new: ");
|
|---|
| 247 | cin >> config::CommonWannier;
|
|---|
| 248 | break;
|
|---|
| 249 | case 'M': // SawtoothStart
|
|---|
| 250 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::SawtoothStart << "\t new: ");
|
|---|
| 251 | cin >> config::SawtoothStart;
|
|---|
| 252 | break;
|
|---|
| 253 | case 'N': // VectorPlane
|
|---|
| 254 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::VectorPlane << "\t new: ");
|
|---|
| 255 | cin >> config::VectorPlane;
|
|---|
| 256 | break;
|
|---|
| 257 | case 'O': // VectorCut
|
|---|
| 258 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::VectorCut << "\t new: ");
|
|---|
| 259 | cin >> config::VectorCut;
|
|---|
| 260 | break;
|
|---|
| 261 | case 'P': // UseAddGramSch
|
|---|
| 262 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::UseAddGramSch << "\t new: ");
|
|---|
| 263 | cin >> config::UseAddGramSch;
|
|---|
| 264 | break;
|
|---|
| 265 | case 'Q': // Seed
|
|---|
| 266 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::Seed << "\t new: ");
|
|---|
| 267 | cin >> config::Seed;
|
|---|
| 268 | break;
|
|---|
| 269 |
|
|---|
| 270 | case 'R': // MaxOuterStep
|
|---|
| 271 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxOuterStep << "\t new: ");
|
|---|
| 272 | cin >> config::MaxOuterStep;
|
|---|
| 273 | break;
|
|---|
| 274 | case 'T': // OutVisStep
|
|---|
| 275 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::OutVisStep << "\t new: ");
|
|---|
| 276 | cin >> config::OutVisStep;
|
|---|
| 277 | break;
|
|---|
| 278 | case 'U': // OutSrcStep
|
|---|
| 279 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::OutSrcStep << "\t new: ");
|
|---|
| 280 | cin >> config::OutSrcStep;
|
|---|
| 281 | break;
|
|---|
| 282 | case 'X': // MaxPsiStep
|
|---|
| 283 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxPsiStep << "\t new: ");
|
|---|
| 284 | cin >> config::MaxPsiStep;
|
|---|
| 285 | break;
|
|---|
| 286 | case 'Y': // EpsWannier
|
|---|
| 287 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::EpsWannier << "\t new: ");
|
|---|
| 288 | cin >> config::EpsWannier;
|
|---|
| 289 | break;
|
|---|
| 290 |
|
|---|
| 291 | case 'Z': // MaxMinStep
|
|---|
| 292 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxMinStep << "\t new: ");
|
|---|
| 293 | cin >> config::MaxMinStep;
|
|---|
| 294 | break;
|
|---|
| 295 | case 'a': // RelEpsTotalEnergy
|
|---|
| 296 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RelEpsTotalEnergy << "\t new: ");
|
|---|
| 297 | cin >> config::RelEpsTotalEnergy;
|
|---|
| 298 | break;
|
|---|
| 299 | case 'b': // RelEpsKineticEnergy
|
|---|
| 300 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RelEpsKineticEnergy << "\t new: ");
|
|---|
| 301 | cin >> config::RelEpsKineticEnergy;
|
|---|
| 302 | break;
|
|---|
| 303 | case 'c': // MaxMinStopStep
|
|---|
| 304 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxMinStopStep << "\t new: ");
|
|---|
| 305 | cin >> config::MaxMinStopStep;
|
|---|
| 306 | break;
|
|---|
| 307 | case 'e': // MaxInitMinStep
|
|---|
| 308 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxInitMinStep << "\t new: ");
|
|---|
| 309 | cin >> config::MaxInitMinStep;
|
|---|
| 310 | break;
|
|---|
| 311 | case 'f': // InitRelEpsTotalEnergy
|
|---|
| 312 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::InitRelEpsTotalEnergy << "\t new: ");
|
|---|
| 313 | cin >> config::InitRelEpsTotalEnergy;
|
|---|
| 314 | break;
|
|---|
| 315 | case 'g': // InitRelEpsKineticEnergy
|
|---|
| 316 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::InitRelEpsKineticEnergy << "\t new: ");
|
|---|
| 317 | cin >> config::InitRelEpsKineticEnergy;
|
|---|
| 318 | break;
|
|---|
| 319 | case 'h': // InitMaxMinStopStep
|
|---|
| 320 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::InitMaxMinStopStep << "\t new: ");
|
|---|
| 321 | cin >> config::InitMaxMinStopStep;
|
|---|
| 322 | break;
|
|---|
| 323 |
|
|---|
| 324 | // case 'j': // BoxLength
|
|---|
| 325 | // Log() << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
|
|---|
| 326 | // double * const cell_size = World::getInstance().getDomain();
|
|---|
| 327 | // for (int i=0;i<6;i++) {
|
|---|
| 328 | // Log() << Verbose(0) << "Cell size" << i << ": ";
|
|---|
| 329 | // cin >> cell_size[i];
|
|---|
| 330 | // }
|
|---|
| 331 | // break;
|
|---|
| 332 |
|
|---|
| 333 | case 'k': // ECut
|
|---|
| 334 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::ECut << "\t new: ");
|
|---|
| 335 | cin >> config::ECut;
|
|---|
| 336 | break;
|
|---|
| 337 | case 'l': // MaxLevel
|
|---|
| 338 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxLevel << "\t new: ");
|
|---|
| 339 | cin >> config::MaxLevel;
|
|---|
| 340 | break;
|
|---|
| 341 | case 'm': // RiemannTensor
|
|---|
| 342 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RiemannTensor << "\t new: ");
|
|---|
| 343 | cin >> config::RiemannTensor;
|
|---|
| 344 | break;
|
|---|
| 345 | case 'n': // LevRFactor
|
|---|
| 346 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::LevRFactor << "\t new: ");
|
|---|
| 347 | cin >> config::LevRFactor;
|
|---|
| 348 | break;
|
|---|
| 349 | case 'o': // RiemannLevel
|
|---|
| 350 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RiemannLevel << "\t new: ");
|
|---|
| 351 | cin >> config::RiemannLevel;
|
|---|
| 352 | break;
|
|---|
| 353 | case 'p': // Lev0Factor
|
|---|
| 354 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::Lev0Factor << "\t new: ");
|
|---|
| 355 | cin >> config::Lev0Factor;
|
|---|
| 356 | break;
|
|---|
| 357 | case 'r': // RTActualUse
|
|---|
| 358 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RTActualUse << "\t new: ");
|
|---|
| 359 | cin >> config::RTActualUse;
|
|---|
| 360 | break;
|
|---|
| 361 | case 's': // PsiType
|
|---|
| 362 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::PsiType << "\t new: ");
|
|---|
| 363 | cin >> config::PsiType;
|
|---|
| 364 | break;
|
|---|
| 365 | case 't': // MaxPsiDouble
|
|---|
| 366 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxPsiDouble << "\t new: ");
|
|---|
| 367 | cin >> config::MaxPsiDouble;
|
|---|
| 368 | break;
|
|---|
| 369 | case 'u': // PsiMaxNoUp
|
|---|
| 370 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::PsiMaxNoUp << "\t new: ");
|
|---|
| 371 | cin >> config::PsiMaxNoUp;
|
|---|
| 372 | break;
|
|---|
| 373 | case 'v': // PsiMaxNoDown
|
|---|
| 374 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::PsiMaxNoDown << "\t new: ");
|
|---|
| 375 | cin >> config::PsiMaxNoDown;
|
|---|
| 376 | break;
|
|---|
| 377 | case 'w': // AddPsis
|
|---|
| 378 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::AddPsis << "\t new: ");
|
|---|
| 379 | cin >> config::AddPsis;
|
|---|
| 380 | break;
|
|---|
| 381 |
|
|---|
| 382 | case 'x': // RCut
|
|---|
| 383 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RCut << "\t new: ");
|
|---|
| 384 | cin >> config::RCut;
|
|---|
| 385 | break;
|
|---|
| 386 | case 'y': // StructOpt
|
|---|
| 387 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::StructOpt << "\t new: ");
|
|---|
| 388 | cin >> config::StructOpt;
|
|---|
| 389 | break;
|
|---|
| 390 | case 'z': // IsAngstroem
|
|---|
| 391 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::IsAngstroem << "\t new: ");
|
|---|
| 392 | cin >> config::IsAngstroem;
|
|---|
| 393 | break;
|
|---|
| 394 | case 'i': // RelativeCoord
|
|---|
| 395 | DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RelativeCoord << "\t new: ");
|
|---|
| 396 | cin >> config::RelativeCoord;
|
|---|
| 397 | break;
|
|---|
| 398 | };
|
|---|
| 399 | } while (choice != 'q');
|
|---|
| 400 | };
|
|---|
| 401 |
|
|---|
| 402 | /** Tests whether a given configuration file adhears to old or new syntax.
|
|---|
| 403 | * \param *filename filename of config file to be tested
|
|---|
| 404 | * \param *periode pointer to a periodentafel class with all elements
|
|---|
| 405 | * \return 0 - old syntax, 1 - new syntax, -1 - unknown syntax
|
|---|
| 406 | */
|
|---|
| 407 | int config::TestSyntax(const char * const filename, const periodentafel * const periode) const
|
|---|
| 408 | {
|
|---|
| 409 | int test;
|
|---|
| 410 | ifstream file(filename);
|
|---|
| 411 |
|
|---|
| 412 | // search file for keyword: ProcPEGamma (new syntax)
|
|---|
| 413 | if (ParseForParameter(1,&file,"ProcPEGamma", 0, 1, 1, int_type, &test, 1, optional)) {
|
|---|
| 414 | file.close();
|
|---|
| 415 | return 1;
|
|---|
| 416 | }
|
|---|
| 417 | // search file for keyword: ProcsGammaPsi (old syntax)
|
|---|
| 418 | if (ParseForParameter(1,&file,"ProcsGammaPsi", 0, 1, 1, int_type, &test, 1, optional)) {
|
|---|
| 419 | file.close();
|
|---|
| 420 | return 0;
|
|---|
| 421 | }
|
|---|
| 422 | file.close();
|
|---|
| 423 | return -1;
|
|---|
| 424 | }
|
|---|
| 425 |
|
|---|
| 426 | /** Returns private config::IsAngstroem.
|
|---|
| 427 | * \return IsAngstroem
|
|---|
| 428 | */
|
|---|
| 429 | bool config::GetIsAngstroem() const
|
|---|
| 430 | {
|
|---|
| 431 | return (IsAngstroem == 1);
|
|---|
| 432 | };
|
|---|
| 433 |
|
|---|
| 434 | /** Returns private config::*defaultpath.
|
|---|
| 435 | * \return *defaultpath
|
|---|
| 436 | */
|
|---|
| 437 | char * config::GetDefaultPath() const
|
|---|
| 438 | {
|
|---|
| 439 | return defaultpath;
|
|---|
| 440 | };
|
|---|
| 441 |
|
|---|
| 442 |
|
|---|
| 443 | /** Returns private config::*defaultpath.
|
|---|
| 444 | * \return *defaultpath
|
|---|
| 445 | */
|
|---|
| 446 | void config::SetDefaultPath(const char * const path)
|
|---|
| 447 | {
|
|---|
| 448 | strcpy(defaultpath, path);
|
|---|
| 449 | };
|
|---|
| 450 |
|
|---|
| 451 | /** Loads a molecule from a ConfigFileBuffer.
|
|---|
| 452 | * \param *mol molecule to load
|
|---|
| 453 | * \param *FileBuffer ConfigFileBuffer to use
|
|---|
| 454 | * \param *periode periodentafel for finding elements
|
|---|
| 455 | * \param FastParsing whether to parse trajectories or not
|
|---|
| 456 | */
|
|---|
| 457 | void LoadMolecule(molecule * const &mol, struct ConfigFileBuffer * const &FileBuffer, const periodentafel * const periode, const bool FastParsing)
|
|---|
| 458 | {
|
|---|
| 459 | int MaxTypes = 0;
|
|---|
| 460 | const element *elementhash[MAX_ELEMENTS];
|
|---|
| 461 | char name[MAX_ELEMENTS];
|
|---|
| 462 | char keyword[MAX_ELEMENTS];
|
|---|
| 463 | int Z = -1;
|
|---|
| 464 | int No[MAX_ELEMENTS];
|
|---|
| 465 | int verbose = 0;
|
|---|
| 466 | double value[3];
|
|---|
| 467 |
|
|---|
| 468 | if (mol == NULL) {
|
|---|
| 469 | DoeLog(0) && (eLog()<< Verbose(0) << "Molecule is not allocated in LoadMolecule(), exit.");
|
|---|
| 470 | performCriticalExit();
|
|---|
| 471 | }
|
|---|
| 472 |
|
|---|
| 473 | ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(MaxTypes), 1, critical);
|
|---|
| 474 | if (MaxTypes == 0) {
|
|---|
| 475 | DoeLog(1) && (eLog()<< Verbose(1) << "There are no atoms according to MaxTypes in this config file." << endl);
|
|---|
| 476 | //performCriticalExit();
|
|---|
| 477 | } else {
|
|---|
| 478 | // prescan number of ions per type
|
|---|
| 479 | DoLog(0) && (Log() << Verbose(0) << "Prescanning ions per type: " << endl);
|
|---|
| 480 | int NoAtoms = 0;
|
|---|
| 481 | for (int i=0; i < MaxTypes; i++) {
|
|---|
| 482 | sprintf(name,"Ion_Type%i",i+1);
|
|---|
| 483 | ParseForParameter(verbose,FileBuffer, (const char*)name, 0, 1, 1, int_type, &No[i], 1, critical);
|
|---|
| 484 | ParseForParameter(verbose,FileBuffer, name, 0, 2, 1, int_type, &Z, 1, critical);
|
|---|
| 485 | elementhash[i] = periode->FindElement(Z);
|
|---|
| 486 | DoLog(1) && (Log() << Verbose(1) << i << ". Z = " << elementhash[i]->getAtomicNumber() << " with " << No[i] << " ions." << endl);
|
|---|
| 487 | NoAtoms += No[i];
|
|---|
| 488 | }
|
|---|
| 489 | int repetition = 0; // which repeated keyword shall be read
|
|---|
| 490 |
|
|---|
| 491 | // sort the lines via the LineMapping
|
|---|
| 492 | sprintf(name,"Ion_Type%i",MaxTypes);
|
|---|
| 493 | if (!ParseForParameter(verbose,FileBuffer, (const char*)name, 1, 1, 1, int_type, &value[0], 1, critical)) {
|
|---|
| 494 | DoeLog(0) && (eLog()<< Verbose(0) << "There are no atoms in the config file!" << endl);
|
|---|
| 495 | performCriticalExit();
|
|---|
| 496 | return;
|
|---|
| 497 | }
|
|---|
| 498 | FileBuffer->CurrentLine++;
|
|---|
| 499 | //Log() << Verbose(0) << FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine]];
|
|---|
| 500 | FileBuffer->MapIonTypesInBuffer(NoAtoms);
|
|---|
| 501 | //for (int i=0; i<(NoAtoms < 100 ? NoAtoms : 100 < 100 ? NoAtoms : 100);++i) {
|
|---|
| 502 | // Log() << Verbose(0) << FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine+i]];
|
|---|
| 503 | //}
|
|---|
| 504 |
|
|---|
| 505 | map<int, atom *> AtomList[MaxTypes];
|
|---|
| 506 | map<int, atom *> LinearList;
|
|---|
| 507 | atom *neues = NULL;
|
|---|
| 508 | Vector position;
|
|---|
| 509 | if (!FastParsing) {
|
|---|
| 510 | // parse in trajectories
|
|---|
| 511 | bool status = true;
|
|---|
| 512 | while (status) {
|
|---|
| 513 | DoLog(0) && (Log() << Verbose(0) << "Currently parsing MD step " << repetition << "." << endl);
|
|---|
| 514 | for (int i=0; i < MaxTypes; i++) {
|
|---|
| 515 | sprintf(name,"Ion_Type%i",i+1);
|
|---|
| 516 | for(int j=0;j<No[i];j++) {
|
|---|
| 517 | sprintf(keyword,"%s_%i",name, j+1);
|
|---|
| 518 | if (repetition == 0) {
|
|---|
| 519 | neues = World::getInstance().createAtom();
|
|---|
| 520 | AtomList[i][j] = neues;
|
|---|
| 521 | LinearList[ FileBuffer->LineMapping[FileBuffer->CurrentLine] ] = neues;
|
|---|
| 522 | neues->setType(elementhash[i]); // find element type
|
|---|
| 523 | } else
|
|---|
| 524 | neues = AtomList[i][j];
|
|---|
| 525 | status = (status &&
|
|---|
| 526 | ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &position[0], 1, (repetition == 0) ? critical : optional) &&
|
|---|
| 527 | ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &position[1], 1, (repetition == 0) ? critical : optional) &&
|
|---|
| 528 | ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &position[2], 1, (repetition == 0) ? critical : optional) &&
|
|---|
| 529 | ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
|
|---|
| 530 | if (!status)
|
|---|
| 531 | break;
|
|---|
| 532 | neues ->setPosition(position);
|
|---|
| 533 |
|
|---|
| 534 | // check size of vectors
|
|---|
| 535 | if (neues->Trajectory.R.size() <= (unsigned int)(repetition)) {
|
|---|
| 536 | //Log() << Verbose(0) << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
|
|---|
| 537 | neues->Trajectory.R.resize(repetition+10);
|
|---|
| 538 | neues->Trajectory.U.resize(repetition+10);
|
|---|
| 539 | neues->Trajectory.F.resize(repetition+10);
|
|---|
| 540 | }
|
|---|
| 541 |
|
|---|
| 542 | // put into trajectories list
|
|---|
| 543 | for (int d=0;d<NDIM;d++)
|
|---|
| 544 | neues->Trajectory.R.at(repetition)[d] = neues->at(d);
|
|---|
| 545 |
|
|---|
| 546 | // parse velocities if present
|
|---|
| 547 | if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->AtomicVelocity[0], 1,optional))
|
|---|
| 548 | neues->AtomicVelocity[0] = 0.;
|
|---|
| 549 | if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->AtomicVelocity[1], 1,optional))
|
|---|
| 550 | neues->AtomicVelocity[1] = 0.;
|
|---|
| 551 | if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->AtomicVelocity[2], 1,optional))
|
|---|
| 552 | neues->AtomicVelocity[2] = 0.;
|
|---|
| 553 | for (int d=0;d<NDIM;d++)
|
|---|
| 554 | neues->Trajectory.U.at(repetition)[d] = neues->AtomicVelocity[d];
|
|---|
| 555 |
|
|---|
| 556 | // parse forces if present
|
|---|
| 557 | if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
|
|---|
| 558 | value[0] = 0.;
|
|---|
| 559 | if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
|
|---|
| 560 | value[1] = 0.;
|
|---|
| 561 | if(!ParseForParameter(verbose,FileBuffer, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
|
|---|
| 562 | value[2] = 0.;
|
|---|
| 563 | for (int d=0;d<NDIM;d++)
|
|---|
| 564 | neues->Trajectory.F.at(repetition)[d] = value[d];
|
|---|
| 565 |
|
|---|
| 566 | // Log() << Verbose(0) << "Parsed position of step " << (repetition) << ": (";
|
|---|
| 567 | // for (int d=0;d<NDIM;d++)
|
|---|
| 568 | // Log() << Verbose(0) << neues->Trajectory.R.at(repetition).x[d] << " "; // next step
|
|---|
| 569 | // Log() << Verbose(0) << ")\t(";
|
|---|
| 570 | // for (int d=0;d<NDIM;d++)
|
|---|
| 571 | // Log() << Verbose(0) << neues->Trajectory.U.at(repetition).x[d] << " "; // next step
|
|---|
| 572 | // Log() << Verbose(0) << ")\t(";
|
|---|
| 573 | // for (int d=0;d<NDIM;d++)
|
|---|
| 574 | // Log() << Verbose(0) << neues->Trajectory.F.at(repetition).x[d] << " "; // next step
|
|---|
| 575 | // Log() << Verbose(0) << ")" << endl;
|
|---|
| 576 | }
|
|---|
| 577 | }
|
|---|
| 578 | repetition++;
|
|---|
| 579 | }
|
|---|
| 580 | repetition--;
|
|---|
| 581 | DoLog(0) && (Log() << Verbose(0) << "Found " << repetition << " trajectory steps." << endl);
|
|---|
| 582 | if (repetition <= 1) // if onyl one step, desactivate use of trajectories
|
|---|
| 583 | mol->MDSteps = 0;
|
|---|
| 584 | else
|
|---|
| 585 | mol->MDSteps = repetition;
|
|---|
| 586 | } else {
|
|---|
| 587 | // find the maximum number of MD steps so that we may parse last one (Ion_Type1_1 must always be present, because is the first atom)
|
|---|
| 588 | repetition = 0;
|
|---|
| 589 | while ( ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) &&
|
|---|
| 590 | ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) &&
|
|---|
| 591 | ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))
|
|---|
| 592 | repetition++;
|
|---|
| 593 | DoLog(0) && (Log() << Verbose(0) << "I found " << repetition << " times the keyword Ion_Type1_1." << endl);
|
|---|
| 594 | // parse in molecule coordinates
|
|---|
| 595 | for (int i=0; i < MaxTypes; i++) {
|
|---|
| 596 | sprintf(name,"Ion_Type%i",i+1);
|
|---|
| 597 | for(int j=0;j<No[i];j++) {
|
|---|
| 598 | sprintf(keyword,"%s_%i",name, j+1);
|
|---|
| 599 | if (repetition == 0) {
|
|---|
| 600 | neues = World::getInstance().createAtom();
|
|---|
| 601 | AtomList[i][j] = neues;
|
|---|
| 602 | LinearList[ FileBuffer->LineMapping[FileBuffer->CurrentLine] ] = neues;
|
|---|
| 603 | neues->setType(elementhash[i]); // find element type
|
|---|
| 604 | } else
|
|---|
| 605 | neues = AtomList[i][j];
|
|---|
| 606 | // then parse for each atom the coordinates as often as present
|
|---|
| 607 | ParseForParameter(verbose,FileBuffer, keyword, 0, 1, 1, double_type, &position[0], repetition,critical);
|
|---|
| 608 | ParseForParameter(verbose,FileBuffer, keyword, 0, 2, 1, double_type, &position[1], repetition,critical);
|
|---|
| 609 | ParseForParameter(verbose,FileBuffer, keyword, 0, 3, 1, double_type, &position[2], repetition,critical);
|
|---|
| 610 | neues->setPosition(position);
|
|---|
| 611 | ParseForParameter(verbose,FileBuffer, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
|
|---|
| 612 | if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 5, 1, double_type, &neues->AtomicVelocity[0], repetition,optional))
|
|---|
| 613 | neues->AtomicVelocity[0] = 0.;
|
|---|
| 614 | if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 6, 1, double_type, &neues->AtomicVelocity[1], repetition,optional))
|
|---|
| 615 | neues->AtomicVelocity[1] = 0.;
|
|---|
| 616 | if(!ParseForParameter(verbose,FileBuffer, keyword, 0, 7, 1, double_type, &neues->AtomicVelocity[2], repetition,optional))
|
|---|
| 617 | neues->AtomicVelocity[2] = 0.;
|
|---|
| 618 | // here we don't care if forces are present (last in trajectories is always equal to current position)
|
|---|
| 619 | neues->setType(elementhash[i]); // find element type
|
|---|
| 620 | mol->AddAtom(neues);
|
|---|
| 621 | }
|
|---|
| 622 | }
|
|---|
| 623 | }
|
|---|
| 624 | // put atoms into the molecule in their original order
|
|---|
| 625 | for(map<int, atom*>::iterator runner = LinearList.begin(); runner != LinearList.end(); ++runner) {
|
|---|
| 626 | mol->AddAtom(runner->second);
|
|---|
| 627 | }
|
|---|
| 628 | }
|
|---|
| 629 | };
|
|---|
| 630 |
|
|---|
| 631 | /** Stores all elements of config structure from which they can be re-read.
|
|---|
| 632 | * \param *filename name of file
|
|---|
| 633 | * \param *periode pointer to a periodentafel class with all elements
|
|---|
| 634 | * \param *mol pointer to molecule containing all atoms of the molecule
|
|---|
| 635 | */
|
|---|
| 636 | bool config::Save(const char * const filename, const periodentafel * const periode, molecule * const mol) const
|
|---|
| 637 | {
|
|---|
| 638 | bool result = true;
|
|---|
| 639 | const Matrix &domain = World::getInstance().getDomain().getM();
|
|---|
| 640 | ofstream * const output = new ofstream(filename, ios::out);
|
|---|
| 641 | if (output != NULL) {
|
|---|
| 642 | *output << "# ParallelCarParinello - main configuration file - created with molecuilder" << endl;
|
|---|
| 643 | *output << endl;
|
|---|
| 644 | *output << "mainname\t" << config::mainname << "\t# programm name (for runtime files)" << endl;
|
|---|
| 645 | *output << "defaultpath\t" << config::defaultpath << "\t# where to put files during runtime" << endl;
|
|---|
| 646 | *output << "pseudopotpath\t" << config::pseudopotpath << "\t# where to find pseudopotentials" << endl;
|
|---|
| 647 | *output << endl;
|
|---|
| 648 | *output << "ProcPEGamma\t" << config::ProcPEGamma << "\t# for parallel computing: share constants" << endl;
|
|---|
| 649 | *output << "ProcPEPsi\t" << config::ProcPEPsi << "\t# for parallel computing: share wave functions" << endl;
|
|---|
| 650 | *output << "DoOutVis\t" << config::DoOutVis << "\t# Output data for OpenDX" << endl;
|
|---|
| 651 | *output << "DoOutMes\t" << config::DoOutMes << "\t# Output data for measurements" << endl;
|
|---|
| 652 | *output << "DoOutOrbitals\t" << config::DoOutOrbitals << "\t# Output all Orbitals" << endl;
|
|---|
| 653 | *output << "DoOutCurr\t" << config::DoOutCurrent << "\t# Ouput current density for OpenDx" << endl;
|
|---|
| 654 | *output << "DoOutNICS\t" << config::DoOutNICS << "\t# Output Nucleus independent current shieldings" << endl;
|
|---|
| 655 | *output << "DoPerturbation\t" << config::DoPerturbation << "\t# Do perturbation calculate and determine susceptibility and shielding" << endl;
|
|---|
| 656 | *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
|
|---|
| 657 | *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
|
|---|
| 658 | *output << "Thermostat\t" << Thermostats->activeThermostat->name() << "\t";
|
|---|
| 659 | *output << Thermostats->activeThermostat->writeParams();
|
|---|
| 660 | *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
|
|---|
| 661 | *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
|
|---|
| 662 | *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;
|
|---|
| 663 | *output << "VectorPlane\t" << config::VectorPlane << "\t# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot" << endl;
|
|---|
| 664 | *output << "VectorCut\t" << config::VectorCut << "\t# Cut plane axis value" << endl;
|
|---|
| 665 | *output << "AddGramSch\t" << config::UseAddGramSch << "\t# Additional GramSchmidtOrtogonalization to be safe" << endl;
|
|---|
| 666 | *output << "Seed\t\t" << config::Seed << "\t# initial value for random seed for Psi coefficients" << endl;
|
|---|
| 667 | *output << endl;
|
|---|
| 668 | *output << "MaxOuterStep\t" << config::MaxOuterStep << "\t# number of MolecularDynamics/Structure optimization steps" << endl;
|
|---|
| 669 | *output << "Deltat\t" << config::Deltat << "\t# time per MD step" << endl;
|
|---|
| 670 | *output << "OutVisStep\t" << config::OutVisStep << "\t# Output visual data every ...th step" << endl;
|
|---|
| 671 | *output << "OutSrcStep\t" << config::OutSrcStep << "\t# Output \"restart\" data every ..th step" << endl;
|
|---|
| 672 | *output << "TargetTemp\t" << Thermostats->TargetTemp << "\t# Target temperature" << endl;
|
|---|
| 673 | *output << "MaxPsiStep\t" << config::MaxPsiStep << "\t# number of Minimisation steps per state (0 - default)" << endl;
|
|---|
| 674 | *output << "EpsWannier\t" << config::EpsWannier << "\t# tolerance value for spread minimisation of orbitals" << endl;
|
|---|
| 675 | *output << endl;
|
|---|
| 676 | *output << "# Values specifying when to stop" << endl;
|
|---|
| 677 | *output << "MaxMinStep\t" << config::MaxMinStep << "\t# Maximum number of steps" << endl;
|
|---|
| 678 | *output << "RelEpsTotalE\t" << config::RelEpsTotalEnergy << "\t# relative change in total energy" << endl;
|
|---|
| 679 | *output << "RelEpsKineticE\t" << config::RelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
|
|---|
| 680 | *output << "MaxMinStopStep\t" << config::MaxMinStopStep << "\t# check every ..th steps" << endl;
|
|---|
| 681 | *output << "MaxMinGapStopStep\t" << config::MaxMinGapStopStep << "\t# check every ..th steps" << endl;
|
|---|
| 682 | *output << endl;
|
|---|
| 683 | *output << "# Values specifying when to stop for INIT, otherwise same as above" << endl;
|
|---|
| 684 | *output << "MaxInitMinStep\t" << config::MaxInitMinStep << "\t# Maximum number of steps" << endl;
|
|---|
| 685 | *output << "InitRelEpsTotalE\t" << config::InitRelEpsTotalEnergy << "\t# relative change in total energy" << endl;
|
|---|
| 686 | *output << "InitRelEpsKineticE\t" << config::InitRelEpsKineticEnergy << "\t# relative change in kinetic energy" << endl;
|
|---|
| 687 | *output << "InitMaxMinStopStep\t" << config::InitMaxMinStopStep << "\t# check every ..th steps" << endl;
|
|---|
| 688 | *output << "InitMaxMinGapStopStep\t" << config::InitMaxMinGapStopStep << "\t# check every ..th steps" << endl;
|
|---|
| 689 | *output << endl;
|
|---|
| 690 | *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
|
|---|
| 691 | *output << domain.at(0,0) << "\t" << endl;
|
|---|
| 692 | *output << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
|
|---|
| 693 | *output << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
|
|---|
| 694 | // FIXME
|
|---|
| 695 | *output << endl;
|
|---|
| 696 | *output << "ECut\t\t" << config::ECut << "\t# energy cutoff for discretization in Hartrees" << endl;
|
|---|
| 697 | *output << "MaxLevel\t" << config::MaxLevel << "\t# number of different levels in the code, >=2" << endl;
|
|---|
| 698 | *output << "Level0Factor\t" << config::Lev0Factor << "\t# factor by which node number increases from S to 0 level" << endl;
|
|---|
| 699 | *output << "RiemannTensor\t" << config::RiemannTensor << "\t# (Use metric)" << endl;
|
|---|
| 700 | switch (config::RiemannTensor) {
|
|---|
| 701 | case 0: //UseNoRT
|
|---|
| 702 | break;
|
|---|
| 703 | case 1: // UseRT
|
|---|
| 704 | *output << "RiemannLevel\t" << config::RiemannLevel << "\t# Number of Riemann Levels" << endl;
|
|---|
| 705 | *output << "LevRFactor\t" << config::LevRFactor << "\t# factor by which node number increases from 0 to R level from" << endl;
|
|---|
| 706 | break;
|
|---|
| 707 | }
|
|---|
| 708 | *output << "PsiType\t\t" << config::PsiType << "\t# 0 - doubly occupied, 1 - SpinUp,SpinDown" << endl;
|
|---|
| 709 | // write out both types for easier changing afterwards
|
|---|
| 710 | // switch (PsiType) {
|
|---|
| 711 | // case 0:
|
|---|
| 712 | *output << "MaxPsiDouble\t" << config::MaxPsiDouble << "\t# here: specifying both maximum number of SpinUp- and -Down-states" << endl;
|
|---|
| 713 | // break;
|
|---|
| 714 | // case 1:
|
|---|
| 715 | *output << "PsiMaxNoUp\t" << config::PsiMaxNoUp << "\t# here: specifying maximum number of SpinUp-states" << endl;
|
|---|
| 716 | *output << "PsiMaxNoDown\t" << config::PsiMaxNoDown << "\t# here: specifying maximum number of SpinDown-states" << endl;
|
|---|
| 717 | // break;
|
|---|
| 718 | // }
|
|---|
| 719 | *output << "AddPsis\t\t" << config::AddPsis << "\t# Additional unoccupied Psis for bandgap determination" << endl;
|
|---|
| 720 | *output << endl;
|
|---|
| 721 | *output << "RCut\t\t" << config::RCut << "\t# R-cut for the ewald summation" << endl;
|
|---|
| 722 | *output << "StructOpt\t" << config::StructOpt << "\t# Do structure optimization beforehand" << endl;
|
|---|
| 723 | *output << "IsAngstroem\t" << config::IsAngstroem << "\t# 0 - Bohr, 1 - Angstroem" << endl;
|
|---|
| 724 | *output << "RelativeCoord\t" << config::RelativeCoord << "\t# whether ion coordinates are relative (1) or absolute (0)" << endl;
|
|---|
| 725 | *output << "MaxTypes\t" << mol->getElementCount() << "\t# maximum number of different ion types" << endl;
|
|---|
| 726 | *output << endl;
|
|---|
| 727 | result = result && mol->Checkout(output);
|
|---|
| 728 | if (mol->MDSteps <=1 )
|
|---|
| 729 | result = result && mol->Output(output);
|
|---|
| 730 | else
|
|---|
| 731 | result = result && mol->OutputTrajectories(output);
|
|---|
| 732 | output->close();
|
|---|
| 733 | output->clear();
|
|---|
| 734 | delete(output);
|
|---|
| 735 | return result;
|
|---|
| 736 | } else {
|
|---|
| 737 | DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open output file:" << filename << endl);
|
|---|
| 738 | return false;
|
|---|
| 739 | }
|
|---|
| 740 | };
|
|---|
| 741 |
|
|---|
| 742 | /** Stores all elements in a MPQC input file.
|
|---|
| 743 | * Note that this format cannot be parsed again.
|
|---|
| 744 | * \param *filename name of file (without ".in" suffix!)
|
|---|
| 745 | * \param *mol pointer to molecule containing all atoms of the molecule
|
|---|
| 746 | */
|
|---|
| 747 | bool config::SaveMPQC(const char * const filename, const molecule * const mol) const
|
|---|
| 748 | {
|
|---|
| 749 | Vector *center = NULL;
|
|---|
| 750 | ofstream *output = NULL;
|
|---|
| 751 |
|
|---|
| 752 | // first without hessian
|
|---|
| 753 | {
|
|---|
| 754 | stringstream * const fname = new stringstream;;
|
|---|
| 755 | *fname << filename << ".in";
|
|---|
| 756 | output = new ofstream(fname->str().c_str(), ios::out);
|
|---|
| 757 | if (output == NULL) {
|
|---|
| 758 | DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc output file:" << fname << endl);
|
|---|
| 759 | delete(fname);
|
|---|
| 760 | return false;
|
|---|
| 761 | }
|
|---|
| 762 | *output << "% Created by MoleCuilder" << endl;
|
|---|
| 763 | *output << "mpqc: (" << endl;
|
|---|
| 764 | *output << "\tsavestate = no" << endl;
|
|---|
| 765 | *output << "\tdo_gradient = yes" << endl;
|
|---|
| 766 | *output << "\tmole<MBPT2>: (" << endl;
|
|---|
| 767 | *output << "\t\tmaxiter = 200" << endl;
|
|---|
| 768 | *output << "\t\tbasis = $:basis" << endl;
|
|---|
| 769 | *output << "\t\tmolecule = $:molecule" << endl;
|
|---|
| 770 | *output << "\t\treference<CLHF>: (" << endl;
|
|---|
| 771 | *output << "\t\t\tbasis = $:basis" << endl;
|
|---|
| 772 | *output << "\t\t\tmolecule = $:molecule" << endl;
|
|---|
| 773 | *output << "\t\t)" << endl;
|
|---|
| 774 | *output << "\t)" << endl;
|
|---|
| 775 | *output << ")" << endl;
|
|---|
| 776 | *output << "molecule<Molecule>: (" << endl;
|
|---|
| 777 | *output << "\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
|
|---|
| 778 | *output << "\t{ atoms geometry } = {" << endl;
|
|---|
| 779 | center = mol->DetermineCenterOfAll();
|
|---|
| 780 | // output of atoms
|
|---|
| 781 | for(molecule::const_iterator iter = mol->begin(); iter!=mol->end();++iter){
|
|---|
| 782 | (*iter)->OutputMPQCLine(output,center);
|
|---|
| 783 | }
|
|---|
| 784 | delete(center);
|
|---|
| 785 | *output << "\t}" << endl;
|
|---|
| 786 | *output << ")" << endl;
|
|---|
| 787 | *output << "basis<GaussianBasisSet>: (" << endl;
|
|---|
| 788 | *output << "\tname = \"" << basis << "\"" << endl;
|
|---|
| 789 | *output << "\tmolecule = $:molecule" << endl;
|
|---|
| 790 | *output << ")" << endl;
|
|---|
| 791 | output->close();
|
|---|
| 792 | delete(output);
|
|---|
| 793 | delete(fname);
|
|---|
| 794 | }
|
|---|
| 795 |
|
|---|
| 796 | // second with hessian
|
|---|
| 797 | {
|
|---|
| 798 | stringstream * const fname = new stringstream;
|
|---|
| 799 | *fname << filename << ".hess.in";
|
|---|
| 800 | output = new ofstream(fname->str().c_str(), ios::out);
|
|---|
| 801 | if (output == NULL) {
|
|---|
| 802 | DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc hessian output file:" << fname << endl);
|
|---|
| 803 | delete(fname);
|
|---|
| 804 | return false;
|
|---|
| 805 | }
|
|---|
| 806 | *output << "% Created by MoleCuilder" << endl;
|
|---|
| 807 | *output << "mpqc: (" << endl;
|
|---|
| 808 | *output << "\tsavestate = no" << endl;
|
|---|
| 809 | *output << "\tdo_gradient = yes" << endl;
|
|---|
| 810 | *output << "\tmole<CLHF>: (" << endl;
|
|---|
| 811 | *output << "\t\tmaxiter = 200" << endl;
|
|---|
| 812 | *output << "\t\tbasis = $:basis" << endl;
|
|---|
| 813 | *output << "\t\tmolecule = $:molecule" << endl;
|
|---|
| 814 | *output << "\t)" << endl;
|
|---|
| 815 | *output << "\tfreq<MolecularFrequencies>: (" << endl;
|
|---|
| 816 | *output << "\t\tmolecule=$:molecule" << endl;
|
|---|
| 817 | *output << "\t)" << endl;
|
|---|
| 818 | *output << ")" << endl;
|
|---|
| 819 | *output << "molecule<Molecule>: (" << endl;
|
|---|
| 820 | *output << "\tunit = " << (IsAngstroem ? "angstrom" : "bohr" ) << endl;
|
|---|
| 821 | *output << "\t{ atoms geometry } = {" << endl;
|
|---|
| 822 | center = mol->DetermineCenterOfAll();
|
|---|
| 823 | // output of atoms
|
|---|
| 824 | for(molecule::const_iterator iter = mol->begin(); iter!=mol->end();++iter){
|
|---|
| 825 | (*iter)->OutputMPQCLine(output,center);
|
|---|
| 826 | }
|
|---|
| 827 | delete(center);
|
|---|
| 828 | *output << "\t}" << endl;
|
|---|
| 829 | *output << ")" << endl;
|
|---|
| 830 | *output << "basis<GaussianBasisSet>: (" << endl;
|
|---|
| 831 | *output << "\tname = \"3-21G\"" << endl;
|
|---|
| 832 | *output << "\tmolecule = $:molecule" << endl;
|
|---|
| 833 | *output << ")" << endl;
|
|---|
| 834 | output->close();
|
|---|
| 835 | delete(output);
|
|---|
| 836 | delete(fname);
|
|---|
| 837 | }
|
|---|
| 838 |
|
|---|
| 839 | return true;
|
|---|
| 840 | };
|
|---|
| 841 |
|
|---|
| 842 | /** Stores all atoms from all molecules in a PDB input file.
|
|---|
| 843 | * Note that this format cannot be parsed again.
|
|---|
| 844 | * \param *filename name of file (without ".in" suffix!)
|
|---|
| 845 | * \param *MolList pointer to MoleculeListClass containing all atoms
|
|---|
| 846 | */
|
|---|
| 847 | bool config::SavePDB(const char * const filename, const MoleculeListClass * const MolList) const
|
|---|
| 848 | {
|
|---|
| 849 | int AtomNo = -1;
|
|---|
| 850 | int MolNo = 0;
|
|---|
| 851 | FILE *f = NULL;
|
|---|
| 852 |
|
|---|
| 853 | char name[MAXSTRINGSIZE];
|
|---|
| 854 | strncpy(name, filename, MAXSTRINGSIZE-1);
|
|---|
| 855 | strncat(name, ".pdb", MAXSTRINGSIZE-(strlen(name)+1));
|
|---|
| 856 | f = fopen(name, "w" );
|
|---|
| 857 | if (f == NULL) {
|
|---|
| 858 | DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open pdb output file:" << name << endl);
|
|---|
| 859 | return false;
|
|---|
| 860 | }
|
|---|
| 861 | fprintf(f, "# Created by MoleCuilder\n");
|
|---|
| 862 |
|
|---|
| 863 | for (MoleculeList::const_iterator MolRunner = MolList->ListOfMolecules.begin(); MolRunner != MolList->ListOfMolecules.end(); MolRunner++) {
|
|---|
| 864 | int *elementNo = new int[MAX_ELEMENTS];
|
|---|
| 865 | for (int i=0;i<MAX_ELEMENTS;i++)
|
|---|
| 866 | elementNo[i] = 0;
|
|---|
| 867 | AtomNo = 0;
|
|---|
| 868 | for (molecule::const_iterator iter = (*MolRunner)->begin(); iter != (*MolRunner)->end(); ++iter) {
|
|---|
| 869 | sprintf(name, "%2s%2d",(*iter)->getType()->getSymbol().c_str(), elementNo[(*iter)->getType()->getAtomicNumber()]);
|
|---|
| 870 | elementNo[(*iter)->getType()->getAtomicNumber()] = (elementNo[(*iter)->getType()->getAtomicNumber()]+1) % 100; // confine to two digits
|
|---|
| 871 | fprintf(f,
|
|---|
| 872 | "ATOM %6u %-4s %4s%c%4u %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s\n",
|
|---|
| 873 | (*iter)->nr, /* atom serial number */
|
|---|
| 874 | name, /* atom name */
|
|---|
| 875 | (*MolRunner)->name, /* residue name */
|
|---|
| 876 | 'a'+(unsigned char)(AtomNo % 26), /* letter for chain */
|
|---|
| 877 | MolNo, /* residue sequence number */
|
|---|
| 878 | (*iter)->at(0), /* position X in Angstroem */
|
|---|
| 879 | (*iter)->at(1), /* position Y in Angstroem */
|
|---|
| 880 | (*iter)->at(2), /* position Z in Angstroem */
|
|---|
| 881 | (double)(*iter)->getType()->getValence(), /* occupancy */
|
|---|
| 882 | (double)(*iter)->getType()->getNoValenceOrbitals(), /* temperature factor */
|
|---|
| 883 | "0", /* segment identifier */
|
|---|
| 884 | (*iter)->getType()->getSymbol().c_str(), /* element symbol */
|
|---|
| 885 | "0"); /* charge */
|
|---|
| 886 | AtomNo++;
|
|---|
| 887 | }
|
|---|
| 888 | delete[](elementNo);
|
|---|
| 889 | MolNo++;
|
|---|
| 890 | }
|
|---|
| 891 | fclose(f);
|
|---|
| 892 |
|
|---|
| 893 | return true;
|
|---|
| 894 | };
|
|---|
| 895 |
|
|---|
| 896 | /** Stores all atoms in a PDB input file.
|
|---|
| 897 | * Note that this format cannot be parsed again.
|
|---|
| 898 | * \param *filename name of file (without ".in" suffix!)
|
|---|
| 899 | * \param *mol pointer to molecule
|
|---|
| 900 | */
|
|---|
| 901 | bool config::SavePDB(const char * const filename, const molecule * const mol) const
|
|---|
| 902 | {
|
|---|
| 903 | int AtomNo = -1;
|
|---|
| 904 | FILE *f = NULL;
|
|---|
| 905 |
|
|---|
| 906 | int *elementNo = new int[MAX_ELEMENTS];
|
|---|
| 907 | for (int i=0;i<MAX_ELEMENTS;i++)
|
|---|
| 908 | elementNo[i] = 0;
|
|---|
| 909 | char name[MAXSTRINGSIZE];
|
|---|
| 910 | strncpy(name, filename, MAXSTRINGSIZE-1);
|
|---|
| 911 | strncat(name, ".pdb", MAXSTRINGSIZE-(strlen(name)+1));
|
|---|
| 912 | f = fopen(name, "w" );
|
|---|
| 913 | if (f == NULL) {
|
|---|
| 914 | DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open pdb output file:" << name << endl);
|
|---|
| 915 | delete[](elementNo);
|
|---|
| 916 | return false;
|
|---|
| 917 | }
|
|---|
| 918 | fprintf(f, "# Created by MoleCuilder\n");
|
|---|
| 919 |
|
|---|
| 920 | AtomNo = 0;
|
|---|
| 921 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
|---|
| 922 | sprintf(name, "%2s%2d",(*iter)->getType()->getSymbol().c_str(), elementNo[(*iter)->getType()->getAtomicNumber()]);
|
|---|
| 923 | elementNo[(*iter)->getType()->getAtomicNumber()] = (elementNo[(*iter)->getType()->getAtomicNumber()]+1) % 100; // confine to two digits
|
|---|
| 924 | fprintf(f,
|
|---|
| 925 | "ATOM %6u %-4s %4s%c%4u %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s%2s\n",
|
|---|
| 926 | (*iter)->nr, /* atom serial number */
|
|---|
| 927 | name, /* atom name */
|
|---|
| 928 | mol->name, /* residue name */
|
|---|
| 929 | 'a'+(unsigned char)(AtomNo % 26), /* letter for chain */
|
|---|
| 930 | 0, /* residue sequence number */
|
|---|
| 931 | (*iter)->at(0), /* position X in Angstroem */
|
|---|
| 932 | (*iter)->at(1), /* position Y in Angstroem */
|
|---|
| 933 | (*iter)->at(2), /* position Z in Angstroem */
|
|---|
| 934 | (double)(*iter)->getType()->getValence(), /* occupancy */
|
|---|
| 935 | (double)(*iter)->getType()->getNoValenceOrbitals(), /* temperature factor */
|
|---|
| 936 | "0", /* segment identifier */
|
|---|
| 937 | (*iter)->getType()->getSymbol().c_str(), /* element symbol */
|
|---|
| 938 | "0"); /* charge */
|
|---|
| 939 | AtomNo++;
|
|---|
| 940 | }
|
|---|
| 941 | fclose(f);
|
|---|
| 942 | delete[](elementNo);
|
|---|
| 943 |
|
|---|
| 944 | return true;
|
|---|
| 945 | };
|
|---|
| 946 |
|
|---|
| 947 | /** Stores all atoms in a TREMOLO data input file.
|
|---|
| 948 | * Note that this format cannot be parsed again.
|
|---|
| 949 | * Note that TREMOLO does not like Id starting at 0, but at 1. Atoms with Id 0 are discarded!
|
|---|
| 950 | * \param *filename name of file (without ".in" suffix!)
|
|---|
| 951 | * \param *mol pointer to molecule
|
|---|
| 952 | */
|
|---|
| 953 | bool config::SaveTREMOLO(const char * const filename, const molecule * const mol) const
|
|---|
| 954 | {
|
|---|
| 955 | ofstream *output = NULL;
|
|---|
| 956 | stringstream * const fname = new stringstream;
|
|---|
| 957 |
|
|---|
| 958 | *fname << filename << ".data";
|
|---|
| 959 | output = new ofstream(fname->str().c_str(), ios::out);
|
|---|
| 960 | if (output == NULL) {
|
|---|
| 961 | DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open tremolo output file:" << fname << endl);
|
|---|
| 962 | delete(fname);
|
|---|
| 963 | return false;
|
|---|
| 964 | }
|
|---|
| 965 |
|
|---|
| 966 | // scan maximum number of neighbours
|
|---|
| 967 | int MaxNeighbours = 0;
|
|---|
| 968 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
|---|
| 969 | const int count = (*iter)->ListOfBonds.size();
|
|---|
| 970 | if (MaxNeighbours < count)
|
|---|
| 971 | MaxNeighbours = count;
|
|---|
| 972 | }
|
|---|
| 973 | *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;
|
|---|
| 974 |
|
|---|
| 975 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
|---|
| 976 | *output << (*iter)->nr << "\t";
|
|---|
| 977 | *output << (*iter)->getName() << "\t";
|
|---|
| 978 | *output << mol->name << "\t";
|
|---|
| 979 | *output << 0 << "\t";
|
|---|
| 980 | *output << (*iter)->at(0) << "\t" << (*iter)->at(1) << "\t" << (*iter)->at(2) << "\t";
|
|---|
| 981 | *output << static_cast<double>((*iter)->getType()->getValence()) << "\t";
|
|---|
| 982 | *output << (*iter)->getType()->getSymbol() << "\t";
|
|---|
| 983 | for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
|
|---|
| 984 | *output << (*runner)->GetOtherAtom(*iter)->nr << "\t";
|
|---|
| 985 | for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
|
|---|
| 986 | *output << "-\t";
|
|---|
| 987 | *output << endl;
|
|---|
| 988 | }
|
|---|
| 989 | output->flush();
|
|---|
| 990 | output->close();
|
|---|
| 991 | delete(output);
|
|---|
| 992 | delete(fname);
|
|---|
| 993 |
|
|---|
| 994 | return true;
|
|---|
| 995 | };
|
|---|
| 996 |
|
|---|
| 997 | /** Stores all atoms from all molecules in a TREMOLO data input file.
|
|---|
| 998 | * Note that this format cannot be parsed again.
|
|---|
| 999 | * Note that TREMOLO does not like Id starting at 0, but at 1. Atoms with Id 0 are discarded!
|
|---|
| 1000 | * \param *filename name of file (without ".in" suffix!)
|
|---|
| 1001 | * \param *MolList pointer to MoleculeListClass containing all atoms
|
|---|
| 1002 | */
|
|---|
| 1003 | bool config::SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const
|
|---|
| 1004 | {
|
|---|
| 1005 | Info FunctionInfo(__func__);
|
|---|
| 1006 | ofstream *output = NULL;
|
|---|
| 1007 | stringstream * const fname = new stringstream;
|
|---|
| 1008 |
|
|---|
| 1009 | *fname << filename << ".data";
|
|---|
| 1010 | output = new ofstream(fname->str().c_str(), ios::out);
|
|---|
| 1011 | if (output == NULL) {
|
|---|
| 1012 | DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open tremolo output file:" << fname << endl);
|
|---|
| 1013 | delete(fname);
|
|---|
| 1014 | return false;
|
|---|
| 1015 | }
|
|---|
| 1016 |
|
|---|
| 1017 | // scan maximum number of neighbours
|
|---|
| 1018 | int MaxNeighbours = 0;
|
|---|
| 1019 | for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
|
|---|
| 1020 | for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
|
|---|
| 1021 | const int count = (*iter)->ListOfBonds.size();
|
|---|
| 1022 | if (MaxNeighbours < count)
|
|---|
| 1023 | MaxNeighbours = count;
|
|---|
| 1024 | }
|
|---|
| 1025 | }
|
|---|
| 1026 | *output << "# ATOMDATA Id name resName resSeq x=3 Charge type neighbors=" << MaxNeighbours << endl;
|
|---|
| 1027 |
|
|---|
| 1028 | // create global to local id map
|
|---|
| 1029 | map<int, int> LocalNotoGlobalNoMap;
|
|---|
| 1030 | {
|
|---|
| 1031 | unsigned int MolCounter = 0;
|
|---|
| 1032 | int AtomNo = 1;
|
|---|
| 1033 | for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
|
|---|
| 1034 | for(molecule::iterator AtomRunner = (*MolWalker)->begin(); AtomRunner != (*MolWalker)->end(); ++AtomRunner) {
|
|---|
| 1035 | LocalNotoGlobalNoMap.insert( pair<int,int>((*AtomRunner)->getId(), AtomNo++) );
|
|---|
| 1036 | }
|
|---|
| 1037 | MolCounter++;
|
|---|
| 1038 | }
|
|---|
| 1039 | ASSERT(MolCounter == MolList->ListOfMolecules.size(), "SaveTREMOLO: LocalNotoGlobalNoMap[] has not been correctly initialized for each molecule");
|
|---|
| 1040 | }
|
|---|
| 1041 |
|
|---|
| 1042 | // write the file
|
|---|
| 1043 | {
|
|---|
| 1044 | int MolCounter = 0;
|
|---|
| 1045 | int AtomNo = 0;
|
|---|
| 1046 | for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
|
|---|
| 1047 | for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
|
|---|
| 1048 | *output << LocalNotoGlobalNoMap[ (*iter)->getId() ] << "\t";
|
|---|
| 1049 | *output << (*iter)->getName() << "\t";
|
|---|
| 1050 | *output << (*MolWalker)->name << "\t";
|
|---|
| 1051 | *output << MolCounter+1 << "\t";
|
|---|
| 1052 | *output << (*iter)->at(0) << "\t" << (*iter)->at(1) << "\t" << (*iter)->at(2) << "\t";
|
|---|
| 1053 | *output << (double)(*iter)->getType()->getValence() << "\t";
|
|---|
| 1054 | *output << (*iter)->getType()->getSymbol() << "\t";
|
|---|
| 1055 | for (BondList::iterator runner = (*iter)->ListOfBonds.begin(); runner != (*iter)->ListOfBonds.end(); runner++)
|
|---|
| 1056 | *output << LocalNotoGlobalNoMap[ (*runner)->GetOtherAtom((*iter))->getId() ] << "\t";
|
|---|
| 1057 | for(int i=(*iter)->ListOfBonds.size(); i < MaxNeighbours; i++)
|
|---|
| 1058 | *output << "-\t";
|
|---|
| 1059 | *output << endl;
|
|---|
| 1060 | AtomNo++;
|
|---|
| 1061 | }
|
|---|
| 1062 | MolCounter++;
|
|---|
| 1063 | }
|
|---|
| 1064 | }
|
|---|
| 1065 |
|
|---|
| 1066 | // store & free
|
|---|
| 1067 | output->flush();
|
|---|
| 1068 | output->close();
|
|---|
| 1069 | delete(output);
|
|---|
| 1070 | delete(fname);
|
|---|
| 1071 |
|
|---|
| 1072 | return true;
|
|---|
| 1073 | };
|
|---|
| 1074 |
|
|---|
| 1075 |
|
|---|
| 1076 | /** Tries given filename or standard on saving the config file.
|
|---|
| 1077 | * \param *ConfigFileName name of file
|
|---|
| 1078 | * \param *periode pointer to periodentafel structure with all the elements
|
|---|
| 1079 | * \param *molecules list of molecules structure with all the atoms and coordinates
|
|---|
| 1080 | */
|
|---|
| 1081 | void config::SaveAll(char *ConfigFileName, periodentafel *periode, MoleculeListClass *molecules)
|
|---|
| 1082 | {
|
|---|
| 1083 | char filename[MAXSTRINGSIZE];
|
|---|
| 1084 | ofstream output;
|
|---|
| 1085 | molecule *mol = NULL;
|
|---|
| 1086 |
|
|---|
| 1087 | // first save as PDB data
|
|---|
| 1088 | if (ConfigFileName != NULL)
|
|---|
| 1089 | strcpy(filename, ConfigFileName);
|
|---|
| 1090 | if (output == NULL)
|
|---|
| 1091 | strcpy(filename,"main_pcp_linux");
|
|---|
| 1092 | Log() << Verbose(0) << "Saving as pdb input ... " << endl;
|
|---|
| 1093 | if (SavePDB(filename, molecules))
|
|---|
| 1094 | Log() << Verbose(0) << "\t... done." << endl;
|
|---|
| 1095 | else
|
|---|
| 1096 | Log() << Verbose(0) << "\t... failed." << endl;
|
|---|
| 1097 |
|
|---|
| 1098 | // then save as tremolo data file
|
|---|
| 1099 | if (ConfigFileName != NULL)
|
|---|
| 1100 | strcpy(filename, ConfigFileName);
|
|---|
| 1101 | if (output == NULL)
|
|---|
| 1102 | strcpy(filename,"main_pcp_linux");
|
|---|
| 1103 | Log() << Verbose(0) << "Saving as tremolo data input ... " << endl;
|
|---|
| 1104 | if (SaveTREMOLO(filename, molecules))
|
|---|
| 1105 | Log() << Verbose(0) << "\t... done." << endl;
|
|---|
| 1106 | else
|
|---|
| 1107 | Log() << Verbose(0) << "\t... failed." << endl;
|
|---|
| 1108 |
|
|---|
| 1109 | // translate each to its center and merge all molecules in MoleculeListClass into this molecule
|
|---|
| 1110 | int N = molecules->ListOfMolecules.size();
|
|---|
| 1111 | if (N != 1) { // don't do anything in case of only one molecule (shifts mol ids otherwise)
|
|---|
| 1112 | int *src = new int[N];
|
|---|
| 1113 | N=0;
|
|---|
| 1114 | for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
|
|---|
| 1115 | src[N++] = (*ListRunner)->IndexNr;
|
|---|
| 1116 | }
|
|---|
| 1117 | mol = World::getInstance().createMolecule();
|
|---|
| 1118 | mol->SetNameFromFilename(ConfigFileName);
|
|---|
| 1119 | //mol->CalculateOrbitals(*this);
|
|---|
| 1120 | delete[](src);
|
|---|
| 1121 | } else {
|
|---|
| 1122 | if (!molecules->ListOfMolecules.empty()) {
|
|---|
| 1123 | mol = *(molecules->ListOfMolecules.begin());
|
|---|
| 1124 | mol->doCountAtoms();
|
|---|
| 1125 | //mol->CalculateOrbitals(*this);
|
|---|
| 1126 | } else {
|
|---|
| 1127 | DoeLog(1) && (eLog() << Verbose(1) << "There are no molecules to save!" << endl);
|
|---|
| 1128 | }
|
|---|
| 1129 | }
|
|---|
| 1130 |
|
|---|
| 1131 | Log() << Verbose(0) << "Storing configuration ... " << endl;
|
|---|
| 1132 | // get correct valence orbitals
|
|---|
| 1133 | if (ConfigFileName != NULL) { // test the file name
|
|---|
| 1134 | strcpy(filename, ConfigFileName);
|
|---|
| 1135 | output.open(filename, ios::trunc);
|
|---|
| 1136 | } else if (strlen(configname) != 0) {
|
|---|
| 1137 | strcpy(filename, configname);
|
|---|
| 1138 | output.open(configname, ios::trunc);
|
|---|
| 1139 | } else {
|
|---|
| 1140 | strcpy(filename, DEFAULTCONFIG);
|
|---|
| 1141 | output.open(DEFAULTCONFIG, ios::trunc);
|
|---|
| 1142 | }
|
|---|
| 1143 | output.close();
|
|---|
| 1144 | output.clear();
|
|---|
| 1145 | Log() << Verbose(0) << "Saving of config file ... " << endl;
|
|---|
| 1146 | if (Save(filename, periode, mol))
|
|---|
| 1147 | Log() << Verbose(0) << "\t... successful." << endl;
|
|---|
| 1148 | else
|
|---|
| 1149 | Log() << Verbose(0) << "\t... failed." << endl;
|
|---|
| 1150 |
|
|---|
| 1151 | // and save to xyz file
|
|---|
| 1152 | if (ConfigFileName != NULL) {
|
|---|
| 1153 | strcpy(filename, ConfigFileName);
|
|---|
| 1154 | strcat(filename, ".xyz");
|
|---|
| 1155 | output.open(filename, ios::trunc);
|
|---|
| 1156 | }
|
|---|
| 1157 | if (output == NULL) {
|
|---|
| 1158 | strcpy(filename,"main_pcp_linux");
|
|---|
| 1159 | strcat(filename, ".xyz");
|
|---|
| 1160 | output.open(filename, ios::trunc);
|
|---|
| 1161 | }
|
|---|
| 1162 | Log() << Verbose(0) << "Saving of XYZ file ... " << endl;
|
|---|
| 1163 | if (mol->MDSteps <= 1) {
|
|---|
| 1164 | if (mol->OutputXYZ(&output))
|
|---|
| 1165 | Log() << Verbose(0) << "\t... successful." << endl;
|
|---|
| 1166 | else
|
|---|
| 1167 | Log() << Verbose(0) << "\t... failed." << endl;
|
|---|
| 1168 | } else {
|
|---|
| 1169 | if (mol->OutputTrajectoriesXYZ(&output))
|
|---|
| 1170 | Log() << Verbose(0) << "\t... successful." << endl;
|
|---|
| 1171 | else
|
|---|
| 1172 | Log() << Verbose(0) << "\t... failed." << endl;
|
|---|
| 1173 | }
|
|---|
| 1174 | output.close();
|
|---|
| 1175 | output.clear();
|
|---|
| 1176 |
|
|---|
| 1177 | // and save as MPQC configuration
|
|---|
| 1178 | if (ConfigFileName != NULL)
|
|---|
| 1179 | strcpy(filename, ConfigFileName);
|
|---|
| 1180 | if (output == NULL)
|
|---|
| 1181 | strcpy(filename,"main_pcp_linux");
|
|---|
| 1182 | Log() << Verbose(0) << "Saving as mpqc input .. " << endl;
|
|---|
| 1183 | if (SaveMPQC(filename, mol))
|
|---|
| 1184 | Log() << Verbose(0) << "\t... done." << endl;
|
|---|
| 1185 | else
|
|---|
| 1186 | Log() << Verbose(0) << "\t... failed." << endl;
|
|---|
| 1187 |
|
|---|
| 1188 | // don't destroy molecule as it contains all our atoms
|
|---|
| 1189 | //World::getInstance().destroyMolecule(mol);
|
|---|
| 1190 | };
|
|---|
| 1191 |
|
|---|
| 1192 | /** Reads parameter from a parsed file.
|
|---|
| 1193 | * The file is either parsed for a certain keyword or if null is given for
|
|---|
| 1194 | * the value in row yth and column xth. If the keyword was necessity#critical,
|
|---|
| 1195 | * then an error is thrown and the programme aborted.
|
|---|
| 1196 | * \warning value is modified (both in contents and position)!
|
|---|
| 1197 | * \param verbose 1 - print found value to stderr, 0 - don't
|
|---|
| 1198 | * \param *file file to be parsed
|
|---|
| 1199 | * \param name Name of value in file (at least 3 chars!)
|
|---|
| 1200 | * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
|
|---|
| 1201 | * (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
|
|---|
| 1202 | * best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
|
|---|
| 1203 | * counted from this unresetted position!)
|
|---|
| 1204 | * \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!)
|
|---|
| 1205 | * \param yth In grid case specifying column number, otherwise the yth \a name matching line
|
|---|
| 1206 | * \param type Type of the Parameter to be read
|
|---|
| 1207 | * \param value address of the value to be read (must have been allocated)
|
|---|
| 1208 | * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply
|
|---|
| 1209 | * \param critical necessity of this keyword being specified (optional, critical)
|
|---|
| 1210 | * \return 1 - found, 0 - not found
|
|---|
| 1211 | * \note Routine is taken from the pcp project and hack-a-slack adapted to C++
|
|---|
| 1212 | */
|
|---|
| 1213 | int ParseForParameter(const int verbose, ifstream * const file, const char * const name, const int sequential, const int xth, const int yth, const int type, void * value, const int repetition, const int critical) {
|
|---|
| 1214 | int i = 0;
|
|---|
| 1215 | int j = 0; // loop variables
|
|---|
| 1216 | int length = 0;
|
|---|
| 1217 | int maxlength = -1;
|
|---|
| 1218 | long file_position = file->tellg(); // mark current position
|
|---|
| 1219 | char *dummy1 = NULL;
|
|---|
| 1220 | char *dummy = NULL;
|
|---|
| 1221 | char free_dummy[MAXSTRINGSIZE]; // pointers in the line that is read in per step
|
|---|
| 1222 | dummy1 = free_dummy;
|
|---|
| 1223 |
|
|---|
| 1224 | //fprintf(stderr,"Parsing for %s\n",name);
|
|---|
| 1225 | if (repetition == 0)
|
|---|
| 1226 | //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
|
|---|
| 1227 | return 0;
|
|---|
| 1228 |
|
|---|
| 1229 | int line = 0; // marks line where parameter was found
|
|---|
| 1230 | int found = (type >= grid) ? 0 : (-yth + 1); // marks if yth parameter name was found
|
|---|
| 1231 | while((found != repetition)) {
|
|---|
| 1232 | dummy1 = dummy = free_dummy;
|
|---|
| 1233 | do {
|
|---|
| 1234 | file->getline(dummy1, 256); // Read the whole line
|
|---|
| 1235 | if (file->eof()) {
|
|---|
| 1236 | if ((critical) && (found == 0)) {
|
|---|
| 1237 | //Error(InitReading, name);
|
|---|
| 1238 | fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
|
|---|
| 1239 | exit(255);
|
|---|
| 1240 | } else {
|
|---|
| 1241 | //if (!sequential)
|
|---|
| 1242 | file->clear();
|
|---|
| 1243 | file->seekg(file_position, ios::beg); // rewind to start position
|
|---|
| 1244 | return 0;
|
|---|
| 1245 | }
|
|---|
| 1246 | }
|
|---|
| 1247 | line++;
|
|---|
| 1248 | } while (dummy != NULL && dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines
|
|---|
| 1249 |
|
|---|
| 1250 | // C++ getline removes newline at end, thus re-add
|
|---|
| 1251 | if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
|
|---|
| 1252 | i = strlen(dummy1);
|
|---|
| 1253 | dummy1[i] = '\n';
|
|---|
| 1254 | dummy1[i+1] = '\0';
|
|---|
| 1255 | }
|
|---|
| 1256 | //fprintf(stderr,"line %i ends at %i, newline at %i\n", line, strlen(dummy1), strchr(dummy1,'\n')-free_dummy);
|
|---|
| 1257 |
|
|---|
| 1258 | if (dummy1 == NULL) {
|
|---|
| 1259 | if (verbose) fprintf(stderr,"Error reading line %i\n",line);
|
|---|
| 1260 | } else {
|
|---|
| 1261 | //fprintf(stderr,"Now parsing the line %i: %s\n", line, dummy1);
|
|---|
| 1262 | }
|
|---|
| 1263 | // Seek for possible end of keyword on line if given ...
|
|---|
| 1264 | if (name != NULL) {
|
|---|
| 1265 | dummy = strchr(dummy1,'\t'); // set dummy on first tab or space which ever's nearer
|
|---|
| 1266 | if (dummy == NULL) {
|
|---|
| 1267 | dummy = strchr(dummy1, ' '); // if not found seek for space
|
|---|
| 1268 | while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' '))) // skip some more tabs and spaces if necessary
|
|---|
| 1269 | dummy++;
|
|---|
| 1270 | }
|
|---|
| 1271 | if (dummy == NULL) {
|
|---|
| 1272 | dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
|
|---|
| 1273 | //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
|
|---|
| 1274 | //Error(FileOpenParams, NULL);
|
|---|
| 1275 | } else {
|
|---|
| 1276 | //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
|
|---|
| 1277 | }
|
|---|
| 1278 | } else dummy = dummy1;
|
|---|
| 1279 | // ... and check if it is the keyword!
|
|---|
| 1280 | //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
|
|---|
| 1281 | if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
|
|---|
| 1282 | found++; // found the parameter!
|
|---|
| 1283 | //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);
|
|---|
| 1284 |
|
|---|
| 1285 | if (found == repetition) {
|
|---|
| 1286 | for (i=0;i<xth;i++) { // i = rows
|
|---|
| 1287 | if (type >= grid) {
|
|---|
| 1288 | // grid structure means that grid starts on the next line, not right after keyword
|
|---|
| 1289 | dummy1 = dummy = free_dummy;
|
|---|
| 1290 | do {
|
|---|
| 1291 | file->getline(dummy1, 256); // Read the whole line, skip commentary and empty ones
|
|---|
| 1292 | if (file->eof()) {
|
|---|
| 1293 | if ((critical) && (found == 0)) {
|
|---|
| 1294 | //Error(InitReading, name);
|
|---|
| 1295 | fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
|
|---|
| 1296 | exit(255);
|
|---|
| 1297 | } else {
|
|---|
| 1298 | //if (!sequential)
|
|---|
| 1299 | file->clear();
|
|---|
| 1300 | file->seekg(file_position, ios::beg); // rewind to start position
|
|---|
| 1301 | return 0;
|
|---|
| 1302 | }
|
|---|
| 1303 | }
|
|---|
| 1304 | line++;
|
|---|
| 1305 | } while ((dummy1[0] == '#') || (dummy1[0] == '\n'));
|
|---|
| 1306 | if (dummy1 == NULL){
|
|---|
| 1307 | if (verbose) fprintf(stderr,"Error reading line %i\n", line);
|
|---|
| 1308 | } else {
|
|---|
| 1309 | //fprintf(stderr,"Reading next line %i: %s\n", line, dummy1);
|
|---|
| 1310 | }
|
|---|
| 1311 | } else { // simple int, strings or doubles start in the same line
|
|---|
| 1312 | while ((*dummy == '\t') || (*dummy == ' ')) // skip interjacent tabs and spaces
|
|---|
| 1313 | dummy++;
|
|---|
| 1314 | }
|
|---|
| 1315 | // C++ getline removes newline at end, thus re-add
|
|---|
| 1316 | if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
|
|---|
| 1317 | j = strlen(dummy1);
|
|---|
| 1318 | dummy1[j] = '\n';
|
|---|
| 1319 | dummy1[j+1] = '\0';
|
|---|
| 1320 | }
|
|---|
| 1321 |
|
|---|
| 1322 | int start = (type >= grid) ? 0 : yth-1 ;
|
|---|
| 1323 | for (j=start;j<yth;j++) { // j = columns
|
|---|
| 1324 | // check for lower triangular area and upper triangular area
|
|---|
| 1325 | if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
|
|---|
| 1326 | *((double *)value) = 0.0;
|
|---|
| 1327 | fprintf(stderr,"%f\t",*((double *)value));
|
|---|
| 1328 | value = (void *)((long)value + sizeof(double));
|
|---|
| 1329 | //value += sizeof(double);
|
|---|
| 1330 | } else {
|
|---|
| 1331 | // otherwise we must skip all interjacent tabs and spaces and find next value
|
|---|
| 1332 | dummy1 = dummy;
|
|---|
| 1333 | dummy = strchr(dummy1, '\t'); // seek for tab or space
|
|---|
| 1334 | if (dummy == NULL)
|
|---|
| 1335 | dummy = strchr(dummy1, ' '); // if not found seek for space
|
|---|
| 1336 | if (dummy == NULL) { // if still zero returned ...
|
|---|
| 1337 | dummy = strchr(dummy1, '\n'); // ... at line end then
|
|---|
| 1338 | if ((j < yth-1) && (type < 4)) { // check if xth value or not yet
|
|---|
| 1339 | if (critical) {
|
|---|
| 1340 | if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
|
|---|
| 1341 | //return 0;
|
|---|
| 1342 | exit(255);
|
|---|
| 1343 | //Error(FileOpenParams, NULL);
|
|---|
| 1344 | } else {
|
|---|
| 1345 | //if (!sequential)
|
|---|
| 1346 | file->clear();
|
|---|
| 1347 | file->seekg(file_position, ios::beg); // rewind to start position
|
|---|
| 1348 | return 0;
|
|---|
| 1349 | }
|
|---|
| 1350 | }
|
|---|
| 1351 | } else {
|
|---|
| 1352 | //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
|
|---|
| 1353 | }
|
|---|
| 1354 | if (*dummy1 == '#') {
|
|---|
| 1355 | // found comment, skipping rest of line
|
|---|
| 1356 | //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
|
|---|
| 1357 | if (!sequential) { // here we need it!
|
|---|
| 1358 | file->seekg(file_position, ios::beg); // rewind to start position
|
|---|
| 1359 | }
|
|---|
| 1360 | return 0;
|
|---|
| 1361 | }
|
|---|
| 1362 | //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
|
|---|
| 1363 | switch(type) {
|
|---|
| 1364 | case (row_int):
|
|---|
| 1365 | *((int *)value) = atoi(dummy1);
|
|---|
| 1366 | if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
|
|---|
| 1367 | if (verbose) fprintf(stderr,"%i\t",*((int *)value));
|
|---|
| 1368 | value = (void *)((long)value + sizeof(int));
|
|---|
| 1369 | //value += sizeof(int);
|
|---|
| 1370 | break;
|
|---|
| 1371 | case(row_double):
|
|---|
| 1372 | case(grid):
|
|---|
| 1373 | case(lower_trigrid):
|
|---|
| 1374 | case(upper_trigrid):
|
|---|
| 1375 | *((double *)value) = atof(dummy1);
|
|---|
| 1376 | if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
|
|---|
| 1377 | if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
|
|---|
| 1378 | value = (void *)((long)value + sizeof(double));
|
|---|
| 1379 | //value += sizeof(double);
|
|---|
| 1380 | break;
|
|---|
| 1381 | case(double_type):
|
|---|
| 1382 | *((double *)value) = atof(dummy1);
|
|---|
| 1383 | if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %lg\n", name, *((double *) value));
|
|---|
| 1384 | //value += sizeof(double);
|
|---|
| 1385 | break;
|
|---|
| 1386 | case(int_type):
|
|---|
| 1387 | *((int *)value) = atoi(dummy1);
|
|---|
| 1388 | if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %i\n", name, *((int *) value));
|
|---|
| 1389 | //value += sizeof(int);
|
|---|
| 1390 | break;
|
|---|
| 1391 | default:
|
|---|
| 1392 | case(string_type):
|
|---|
| 1393 | if (value != NULL) {
|
|---|
| 1394 | //if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
|
|---|
| 1395 | maxlength = MAXSTRINGSIZE;
|
|---|
| 1396 | length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
|
|---|
| 1397 | strncpy((char *)value, dummy1, length); // copy as much
|
|---|
| 1398 | ((char *)value)[length] = '\0'; // and set end marker
|
|---|
| 1399 | if ((verbose) && (i == xth-1)) fprintf(stderr,"%s is '%s' (%i chars)\n",name,((char *) value), length);
|
|---|
| 1400 | //value += sizeof(char);
|
|---|
| 1401 | } else {
|
|---|
| 1402 | }
|
|---|
| 1403 | break;
|
|---|
| 1404 | }
|
|---|
| 1405 | }
|
|---|
| 1406 | while (*dummy == '\t')
|
|---|
| 1407 | dummy++;
|
|---|
| 1408 | }
|
|---|
| 1409 | }
|
|---|
| 1410 | }
|
|---|
| 1411 | }
|
|---|
| 1412 | }
|
|---|
| 1413 | if ((type >= row_int) && (verbose))
|
|---|
| 1414 | fprintf(stderr,"\n");
|
|---|
| 1415 | if (!sequential) {
|
|---|
| 1416 | file->clear();
|
|---|
| 1417 | file->seekg(file_position, ios::beg); // rewind to start position
|
|---|
| 1418 | }
|
|---|
| 1419 | //fprintf(stderr, "End of Parsing\n\n");
|
|---|
| 1420 |
|
|---|
| 1421 | return (found); // true if found, false if not
|
|---|
| 1422 | }
|
|---|
| 1423 |
|
|---|
| 1424 |
|
|---|
| 1425 | /** Reads parameter from a parsed file.
|
|---|
| 1426 | * The file is either parsed for a certain keyword or if null is given for
|
|---|
| 1427 | * the value in row yth and column xth. If the keyword was necessity#critical,
|
|---|
| 1428 | * then an error is thrown and the programme aborted.
|
|---|
| 1429 | * \warning value is modified (both in contents and position)!
|
|---|
| 1430 | * \param verbose 1 - print found value to stderr, 0 - don't
|
|---|
| 1431 | * \param *FileBuffer pointer to buffer structure
|
|---|
| 1432 | * \param name Name of value in file (at least 3 chars!)
|
|---|
| 1433 | * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
|
|---|
| 1434 | * (if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
|
|---|
| 1435 | * best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
|
|---|
| 1436 | * counted from this unresetted position!)
|
|---|
| 1437 | * \param xth Which among a number of parameters it is (in grid case it's row number as grid is read as a whole!)
|
|---|
| 1438 | * \param yth In grid case specifying column number, otherwise the yth \a name matching line
|
|---|
| 1439 | * \param type Type of the Parameter to be read
|
|---|
| 1440 | * \param value address of the value to be read (must have been allocated)
|
|---|
| 1441 | * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply
|
|---|
| 1442 | * \param critical necessity of this keyword being specified (optional, critical)
|
|---|
| 1443 | * \return 1 - found, 0 - not found
|
|---|
| 1444 | * \note Routine is taken from the pcp project and hack-a-slack adapted to C++
|
|---|
| 1445 | */
|
|---|
| 1446 | int ParseForParameter(const int verbose, struct ConfigFileBuffer * const FileBuffer, const char * const name, const int sequential, const int xth, const int yth, const int type, void * value, const int repetition, const int critical) {
|
|---|
| 1447 | int i = 0;
|
|---|
| 1448 | int j = 0; // loop variables
|
|---|
| 1449 | int length = 0;
|
|---|
| 1450 | int maxlength = -1;
|
|---|
| 1451 | int OldCurrentLine = FileBuffer->CurrentLine;
|
|---|
| 1452 | char *dummy1 = NULL;
|
|---|
| 1453 | char *dummy = NULL; // pointers in the line that is read in per step
|
|---|
| 1454 |
|
|---|
| 1455 | //fprintf(stderr,"Parsing for %s\n",name);
|
|---|
| 1456 | if (repetition == 0)
|
|---|
| 1457 | //Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
|
|---|
| 1458 | return 0;
|
|---|
| 1459 |
|
|---|
| 1460 | int line = 0; // marks line where parameter was found
|
|---|
| 1461 | int found = (type >= grid) ? 0 : (-yth + 1); // marks if yth parameter name was found
|
|---|
| 1462 | while((found != repetition)) {
|
|---|
| 1463 | dummy1 = dummy = NULL;
|
|---|
| 1464 | do {
|
|---|
| 1465 | dummy1 = FileBuffer->buffer[ FileBuffer->LineMapping[FileBuffer->CurrentLine++] ];
|
|---|
| 1466 | if (FileBuffer->CurrentLine >= FileBuffer->NoLines) {
|
|---|
| 1467 | if ((critical) && (found == 0)) {
|
|---|
| 1468 | //Error(InitReading, name);
|
|---|
| 1469 | fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
|
|---|
| 1470 | exit(255);
|
|---|
| 1471 | } else {
|
|---|
| 1472 | FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
|
|---|
| 1473 | return 0;
|
|---|
| 1474 | }
|
|---|
| 1475 | }
|
|---|
| 1476 | if (dummy1 == NULL) {
|
|---|
| 1477 | if (verbose) fprintf(stderr,"Error reading line %i\n",line);
|
|---|
| 1478 | } else {
|
|---|
| 1479 | //fprintf(stderr,"Now parsing the line %i: %s\n", line, dummy1);
|
|---|
| 1480 | }
|
|---|
| 1481 | line++;
|
|---|
| 1482 | } while (dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines
|
|---|
| 1483 |
|
|---|
| 1484 | // Seek for possible end of keyword on line if given ...
|
|---|
| 1485 | if (name != NULL) {
|
|---|
| 1486 | dummy = strchr(dummy1,'\t'); // set dummy on first tab or space which ever's nearer
|
|---|
| 1487 | if (dummy == NULL) {
|
|---|
| 1488 | dummy = strchr(dummy1, ' '); // if not found seek for space
|
|---|
| 1489 | while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' '))) // skip some more tabs and spaces if necessary
|
|---|
| 1490 | dummy++;
|
|---|
| 1491 | }
|
|---|
| 1492 | if (dummy == NULL) {
|
|---|
| 1493 | dummy = strchr(dummy1, '\n'); // set on line end then (whole line = keyword)
|
|---|
| 1494 | //fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
|
|---|
| 1495 | //Error(FileOpenParams, NULL);
|
|---|
| 1496 | } else {
|
|---|
| 1497 | //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
|
|---|
| 1498 | }
|
|---|
| 1499 | } else dummy = dummy1;
|
|---|
| 1500 | // ... and check if it is the keyword!
|
|---|
| 1501 | //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
|
|---|
| 1502 | if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
|
|---|
| 1503 | found++; // found the parameter!
|
|---|
| 1504 | //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);
|
|---|
| 1505 |
|
|---|
| 1506 | if (found == repetition) {
|
|---|
| 1507 | for (i=0;i<xth;i++) { // i = rows
|
|---|
| 1508 | if (type >= grid) {
|
|---|
| 1509 | // grid structure means that grid starts on the next line, not right after keyword
|
|---|
| 1510 | dummy1 = dummy = NULL;
|
|---|
| 1511 | do {
|
|---|
| 1512 | dummy1 = FileBuffer->buffer[ FileBuffer->LineMapping[ FileBuffer->CurrentLine++] ];
|
|---|
| 1513 | if (FileBuffer->CurrentLine >= FileBuffer->NoLines) {
|
|---|
| 1514 | if ((critical) && (found == 0)) {
|
|---|
| 1515 | //Error(InitReading, name);
|
|---|
| 1516 | fprintf(stderr,"Error:InitReading, critical %s not found\n", name);
|
|---|
| 1517 | exit(255);
|
|---|
| 1518 | } else {
|
|---|
| 1519 | FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
|
|---|
| 1520 | return 0;
|
|---|
| 1521 | }
|
|---|
| 1522 | }
|
|---|
| 1523 | if (dummy1 == NULL) {
|
|---|
| 1524 | if (verbose) fprintf(stderr,"Error reading line %i\n", line);
|
|---|
| 1525 | } else {
|
|---|
| 1526 | //fprintf(stderr,"Reading next line %i: %s\n", line, dummy1);
|
|---|
| 1527 | }
|
|---|
| 1528 | line++;
|
|---|
| 1529 | } while ((dummy1 != NULL) && ((dummy1[0] == '#') || (dummy1[0] == '\n')));
|
|---|
| 1530 | dummy = dummy1;
|
|---|
| 1531 | } else { // simple int, strings or doubles start in the same line
|
|---|
| 1532 | while ((*dummy == '\t') || (*dummy == ' ')) // skip interjacent tabs and spaces
|
|---|
| 1533 | dummy++;
|
|---|
| 1534 | }
|
|---|
| 1535 |
|
|---|
| 1536 | for (j=((type >= grid) ? 0 : yth-1);j<yth;j++) { // j = columns
|
|---|
| 1537 | // check for lower triangular area and upper triangular area
|
|---|
| 1538 | if ( ((i > j) && (type == upper_trigrid)) || ((j > i) && (type == lower_trigrid))) {
|
|---|
| 1539 | *((double *)value) = 0.0;
|
|---|
| 1540 | fprintf(stderr,"%f\t",*((double *)value));
|
|---|
| 1541 | value = (void *)((long)value + sizeof(double));
|
|---|
| 1542 | //value += sizeof(double);
|
|---|
| 1543 | } else {
|
|---|
| 1544 | // otherwise we must skip all interjacent tabs and spaces and find next value
|
|---|
| 1545 | dummy1 = dummy;
|
|---|
| 1546 | dummy = strchr(dummy1, '\t'); // seek for tab or space
|
|---|
| 1547 | if (dummy == NULL)
|
|---|
| 1548 | dummy = strchr(dummy1, ' '); // if not found seek for space
|
|---|
| 1549 | if (dummy == NULL) { // if still zero returned ...
|
|---|
| 1550 | dummy = strchr(dummy1, '\n'); // ... at line end then
|
|---|
| 1551 | if ((j < yth-1) && (type < 4)) { // check if xth value or not yet
|
|---|
| 1552 | if (critical) {
|
|---|
| 1553 | if (verbose) fprintf(stderr,"Error: EoL at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
|
|---|
| 1554 | //return 0;
|
|---|
| 1555 | exit(255);
|
|---|
| 1556 | //Error(FileOpenParams, NULL);
|
|---|
| 1557 | } else {
|
|---|
| 1558 | if (!sequential) { // here we need it!
|
|---|
| 1559 | FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
|
|---|
| 1560 | }
|
|---|
| 1561 | return 0;
|
|---|
| 1562 | }
|
|---|
| 1563 | }
|
|---|
| 1564 | } else {
|
|---|
| 1565 | //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
|
|---|
| 1566 | }
|
|---|
| 1567 | if (*dummy1 == '#') {
|
|---|
| 1568 | // found comment, skipping rest of line
|
|---|
| 1569 | //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
|
|---|
| 1570 | if (!sequential) { // here we need it!
|
|---|
| 1571 | FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
|
|---|
| 1572 | }
|
|---|
| 1573 | return 0;
|
|---|
| 1574 | }
|
|---|
| 1575 | //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
|
|---|
| 1576 | switch(type) {
|
|---|
| 1577 | case (row_int):
|
|---|
| 1578 | *((int *)value) = atoi(dummy1);
|
|---|
| 1579 | if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
|
|---|
| 1580 | if (verbose) fprintf(stderr,"%i\t",*((int *)value));
|
|---|
| 1581 | value = (void *)((long)value + sizeof(int));
|
|---|
| 1582 | //value += sizeof(int);
|
|---|
| 1583 | break;
|
|---|
| 1584 | case(row_double):
|
|---|
| 1585 | case(grid):
|
|---|
| 1586 | case(lower_trigrid):
|
|---|
| 1587 | case(upper_trigrid):
|
|---|
| 1588 | *((double *)value) = atof(dummy1);
|
|---|
| 1589 | if ((verbose) && (i==0) && (j==0)) fprintf(stderr,"%s = ", name);
|
|---|
| 1590 | if (verbose) fprintf(stderr,"%lg\t",*((double *)value));
|
|---|
| 1591 | value = (void *)((long)value + sizeof(double));
|
|---|
| 1592 | //value += sizeof(double);
|
|---|
| 1593 | break;
|
|---|
| 1594 | case(double_type):
|
|---|
| 1595 | *((double *)value) = atof(dummy1);
|
|---|
| 1596 | if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %lg\n", name, *((double *) value));
|
|---|
| 1597 | //value += sizeof(double);
|
|---|
| 1598 | break;
|
|---|
| 1599 | case(int_type):
|
|---|
| 1600 | *((int *)value) = atoi(dummy1);
|
|---|
| 1601 | if ((verbose) && (i == xth-1)) fprintf(stderr,"%s = %i\n", name, *((int *) value));
|
|---|
| 1602 | //value += sizeof(int);
|
|---|
| 1603 | break;
|
|---|
| 1604 | default:
|
|---|
| 1605 | case(string_type):
|
|---|
| 1606 | if (value != NULL) {
|
|---|
| 1607 | //if (maxlength == -1) maxlength = strlen((char *)value); // get maximum size of string array
|
|---|
| 1608 | maxlength = MAXSTRINGSIZE;
|
|---|
| 1609 | length = maxlength > (dummy-dummy1) ? (dummy-dummy1) : maxlength; // cap at maximum
|
|---|
| 1610 | strncpy((char *)value, dummy1, length); // copy as much
|
|---|
| 1611 | ((char *)value)[length] = '\0'; // and set end marker
|
|---|
| 1612 | if ((verbose) && (i == xth-1)) fprintf(stderr,"%s is '%s' (%i chars)\n",name,((char *) value), length);
|
|---|
| 1613 | //value += sizeof(char);
|
|---|
| 1614 | } else {
|
|---|
| 1615 | }
|
|---|
| 1616 | break;
|
|---|
| 1617 | }
|
|---|
| 1618 | }
|
|---|
| 1619 | while (*dummy == '\t')
|
|---|
| 1620 | dummy++;
|
|---|
| 1621 | }
|
|---|
| 1622 | }
|
|---|
| 1623 | }
|
|---|
| 1624 | }
|
|---|
| 1625 | }
|
|---|
| 1626 | if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
|
|---|
| 1627 | if (!sequential) {
|
|---|
| 1628 | FileBuffer->CurrentLine = OldCurrentLine; // rewind to start position
|
|---|
| 1629 | }
|
|---|
| 1630 | //fprintf(stderr, "End of Parsing\n\n");
|
|---|
| 1631 |
|
|---|
| 1632 | return (found); // true if found, false if not
|
|---|
| 1633 | }
|
|---|
| 1634 |
|
|---|
| 1635 | /** Reading of Thermostat related values from parameter file.
|
|---|
| 1636 | * \param *fb file buffer containing the config file
|
|---|
| 1637 | */
|
|---|
| 1638 | void config::ParseThermostats(class ConfigFileBuffer * const fb)
|
|---|
| 1639 | {
|
|---|
| 1640 | char * const thermo = new char[12];
|
|---|
| 1641 | const int verbose = 0;
|
|---|
| 1642 |
|
|---|
| 1643 | // read desired Thermostat from file along with needed additional parameters
|
|---|
| 1644 | if (ParseForParameter(verbose,fb,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
|
|---|
| 1645 | Thermostats->makeActive(thermo,fb);
|
|---|
| 1646 | } else {
|
|---|
| 1647 | if ((Thermostats->TargetTemp != 0))
|
|---|
| 1648 | DoLog(2) && (Log() << Verbose(2) << "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
|
|---|
| 1649 | Thermostats->chooseNone();
|
|---|
| 1650 | }
|
|---|
| 1651 | delete[](thermo);
|
|---|
| 1652 | };
|
|---|
| 1653 |
|
|---|
