Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/config.cpp

    re138de ra67d19  
    55 */
    66
     7#include <stdio.h>
     8#include <cstring>
     9
    710#include "atom.hpp"
     11#include "bond.hpp"
    812#include "config.hpp"
    913#include "element.hpp"
    … …  
    1519#include "molecule.hpp"
    1620#include "periodentafel.hpp"
     21#include "World.hpp"
    1722
    1823/******************************** Functions for class ConfigFileBuffer **********************/
    … …  
    2429    char number1[8];
    2530    char number2[8];
    26     char *dummy1, *dummy2;
     31    const char *dummy1, *dummy2;
    2732    //Log() << Verbose(0) << s1 << "  " << s2 << endl;
    2833    dummy1 = strchr(s1, '_')+sizeof(char)*5;  // go just after "Ion_Type"
    … …  
    6974  file= new ifstream(filename);
    7075  if (file == NULL) {
    71     eLog() << Verbose(0) << "ERROR: config file " << filename << " missing!" << endl;
     76    DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
    7277    return;
    7378  }
    … …  
    8085  file->clear();
    8186  file->seekg(file_position, ios::beg);
    82   Log() << Verbose(1) << NoLines-1 << " lines were recognized." << endl;
     87  DoLog(1) && (Log() << Verbose(1) << NoLines-1 << " lines were recognized." << endl);
    8388
    8489  // allocate buffer's 1st dimension
    8590  if (buffer != NULL) {
    86     eLog() << Verbose(0) << "ERROR: FileBuffer->buffer is not NULL!" << endl;
     91    DoeLog(1) && (eLog()<< Verbose(1) << "FileBuffer->buffer is not NULL!" << endl);
    8792    return;
    8893  } else
    … …  
    100105    lines++;
    101106  } while((!file->eof()) && (lines < NoLines));
    102   Log() << Verbose(1) << lines-1 << " lines were read into the buffer." << endl;
     107  DoLog(1) && (Log() << Verbose(1) << lines-1 << " lines were read into the buffer." << endl);
    103108
    104109  // close and exit
    … …  
    137142void ConfigFileBuffer::MapIonTypesInBuffer(const int NoAtoms)
    138143{
    139   map<const char *, int, IonTypeCompare> LineList;
     144  map<const char *, int, IonTypeCompare> IonTypeLineMap;
    140145  if (LineMapping == NULL) {
    141     eLog() << Verbose(0) << "map pointer is NULL: " << LineMapping << endl;
     146    DoeLog(0) && (eLog()<< Verbose(0) << "map pointer is NULL: " << LineMapping << endl);
     147    performCriticalExit();
    142148    return;
    143149  }
    … …  
    145151  // put all into hashed map
    146152  for (int i=0; i<NoAtoms; ++i) {
    147     LineList.insert(pair<const char *, int> (buffer[CurrentLine+i], CurrentLine+i));
     153    IonTypeLineMap.insert(pair<const char *, int> (buffer[CurrentLine+i], CurrentLine+i));
    148154  }
    149155
    150156  // fill map
    151157  int nr=0;
    152   for (map<const char *, int, IonTypeCompare>::iterator runner = LineList.begin(); runner != LineList.end(); ++runner) {
     158  for (map<const char *, int, IonTypeCompare>::iterator runner = IonTypeLineMap.begin(); runner != IonTypeLineMap.end(); ++runner) {
    153159    if (CurrentLine+nr < NoLines)
    154160      LineMapping[CurrentLine+(nr++)] = runner->second;
    155     else
    156       eLog() << Verbose(0) << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl;
     161    else {
     162      DoeLog(0) && (eLog()<< Verbose(0) << "config::MapIonTypesInBuffer - NoAtoms is wrong: We are past the end of the file!" << endl);
     163      performCriticalExit();
     164    }
    157165  }
    158166}
    … …  
    200208    Free(&ThermostatNames[j]);
    201209  Free(&ThermostatNames);
     210
     211  if (BG != NULL)
     212    delete(BG);
    202213};
    203214
    … …  
    239250        Thermostat = None;
    240251      } else {
    241         Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     252        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    242253        Thermostat = None;
    243254      }
    … …  
    247258        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
    248259      } else {
    249         Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     260        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    250261        Thermostat = None;
    251262      }
    … …  
    255266        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
    256267      } else {
    257         Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     268        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    258269        Thermostat = None;
    259270      }
    … …  
    263274        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
    264275        if (ParseForParameter(verbose,fb,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
    265           Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
     276          DoLog(2) && (Log() << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl);
    266277        } else {
    267278          alpha = 1.;
    268279        }
    269280      } else {
    270         Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     281        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    271282        Thermostat = None;
    272283      }
    … …  
    276287        ParseForParameter(verbose,fb,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
    277288      } else {
    278         Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     289        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    279290        Thermostat = None;
    280291      }
    … …  
    285296        alpha = 0.;
    286297      } else {
    287         Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     298        DoLog(1) && (Log() << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl);
    288299        Thermostat = None;
    289300      }
    290301    } else {
    291       Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl;
     302      DoLog(1) && (Log() << Verbose(1) << " Warning: thermostat name was not understood!" << endl);
    292303      Thermostat = None;
    293304    }
    294305  } else {
    295306    if ((MaxOuterStep > 0) && (TargetTemp != 0))
    296       Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl;
     307      DoLog(2) && (Log() << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl);
    297308    Thermostat = None;
    298309  }
    … …  
    310321
    311322  do {
    312     Log() << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl;
    313     Log() << Verbose(0) << " A - mainname (prefix for all runtime files)" << endl;
    314     Log() << Verbose(0) << " B - Default path (for runtime files)" << endl;
    315     Log() << Verbose(0) << " C - Path of pseudopotential files" << endl;
    316     Log() << Verbose(0) << " D - Number of coefficient sharing processes" << endl;
    317     Log() << Verbose(0) << " E - Number of wave function sharing processes" << endl;
    318     Log() << Verbose(0) << " F - 0: Don't output density for OpenDX, 1: do" << endl;
    319     Log() << Verbose(0) << " G - 0: Don't output physical data, 1: do" << endl;
    320     Log() << Verbose(0) << " H - 0: Don't output densities of each unperturbed orbital for OpenDX, 1: do" << endl;
    321     Log() << Verbose(0) << " I - 0: Don't output current density for OpenDX, 1: do" << endl;
    322     Log() << Verbose(0) << " J - 0: Don't do the full current calculation, 1: do" << endl;
    323     Log() << Verbose(0) << " K - 0: Don't do perturbation calculation to obtain susceptibility and shielding, 1: do" << endl;
    324     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;
    325     Log() << Verbose(0) << " M - Absolute begin of unphysical sawtooth transfer for position operator within cell" << endl;
    326     Log() << Verbose(0) << " N - (0,1,2) x,y,z-plane to do two-dimensional current vector cut" << endl;
    327     Log() << Verbose(0) << " O - Absolute position along vector cut axis for cut plane" << endl;
    328     Log() << Verbose(0) << " P - Additional Gram-Schmidt-Orthonormalization to stabilize numerics" << endl;
    329     Log() << Verbose(0) << " Q - Initial integer value of random number generator" << endl;
    330     Log() << Verbose(0) << " R - for perturbation 0, for structure optimization defines upper limit of iterations" << endl;
    331     Log() << Verbose(0) << " T - Output visual after ...th step" << endl;
    332     Log() << Verbose(0) << " U - Output source densities of wave functions after ...th step" << endl;
    333     Log() << Verbose(0) << " X - minimization iterations per wave function, if unsure leave at default value 0" << endl;
    334     Log() << Verbose(0) << " Y - tolerance value for total spread in iterative Jacobi diagonalization" << endl;
    335     Log() << Verbose(0) << " Z - Maximum number of minimization iterations" << endl;
    336     Log() << Verbose(0) << " a - Relative change in total energy to stop min. iteration" << endl;
    337     Log() << Verbose(0) << " b - Relative change in kinetic energy to stop min. iteration" << endl;
    338     Log() << Verbose(0) << " c - Check stop conditions every ..th step during min. iteration" << endl;
    339     Log() << Verbose(0) << " e - Maximum number of minimization iterations during initial level" << endl;
    340     Log() << Verbose(0) << " f - Relative change in total energy to stop min. iteration during initial level" << endl;
    341     Log() << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl;
    342     Log() << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl;
     323    DoLog(0) && (Log() << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl);
     324    DoLog(0) && (Log() << Verbose(0) << " A - mainname (prefix for all runtime files)" << endl);
     325    DoLog(0) && (Log() << Verbose(0) << " B - Default path (for runtime files)" << endl);
     326    DoLog(0) && (Log() << Verbose(0) << " C - Path of pseudopotential files" << endl);
     327    DoLog(0) && (Log() << Verbose(0) << " D - Number of coefficient sharing processes" << endl);
     328    DoLog(0) && (Log() << Verbose(0) << " E - Number of wave function sharing processes" << endl);
     329    DoLog(0) && (Log() << Verbose(0) << " F - 0: Don't output density for OpenDX, 1: do" << endl);
     330    DoLog(0) && (Log() << Verbose(0) << " G - 0: Don't output physical data, 1: do" << endl);
     331    DoLog(0) && (Log() << Verbose(0) << " H - 0: Don't output densities of each unperturbed orbital for OpenDX, 1: do" << endl);
     332    DoLog(0) && (Log() << Verbose(0) << " I - 0: Don't output current density for OpenDX, 1: do" << endl);
     333    DoLog(0) && (Log() << Verbose(0) << " J - 0: Don't do the full current calculation, 1: do" << endl);
     334    DoLog(0) && (Log() << Verbose(0) << " K - 0: Don't do perturbation calculation to obtain susceptibility and shielding, 1: do" << endl);
     335    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);
     336    DoLog(0) && (Log() << Verbose(0) << " M - Absolute begin of unphysical sawtooth transfer for position operator within cell" << endl);
     337    DoLog(0) && (Log() << Verbose(0) << " N - (0,1,2) x,y,z-plane to do two-dimensional current vector cut" << endl);
     338    DoLog(0) && (Log() << Verbose(0) << " O - Absolute position along vector cut axis for cut plane" << endl);
     339    DoLog(0) && (Log() << Verbose(0) << " P - Additional Gram-Schmidt-Orthonormalization to stabilize numerics" << endl);
     340    DoLog(0) && (Log() << Verbose(0) << " Q - Initial integer value of random number generator" << endl);
     341    DoLog(0) && (Log() << Verbose(0) << " R - for perturbation 0, for structure optimization defines upper limit of iterations" << endl);
     342    DoLog(0) && (Log() << Verbose(0) << " T - Output visual after ...th step" << endl);
     343    DoLog(0) && (Log() << Verbose(0) << " U - Output source densities of wave functions after ...th step" << endl);
     344    DoLog(0) && (Log() << Verbose(0) << " X - minimization iterations per wave function, if unsure leave at default value 0" << endl);
     345    DoLog(0) && (Log() << Verbose(0) << " Y - tolerance value for total spread in iterative Jacobi diagonalization" << endl);
     346    DoLog(0) && (Log() << Verbose(0) << " Z - Maximum number of minimization iterations" << endl);
     347    DoLog(0) && (Log() << Verbose(0) << " a - Relative change in total energy to stop min. iteration" << endl);
     348    DoLog(0) && (Log() << Verbose(0) << " b - Relative change in kinetic energy to stop min. iteration" << endl);
     349    DoLog(0) && (Log() << Verbose(0) << " c - Check stop conditions every ..th step during min. iteration" << endl);
     350    DoLog(0) && (Log() << Verbose(0) << " e - Maximum number of minimization iterations during initial level" << endl);
     351    DoLog(0) && (Log() << Verbose(0) << " f - Relative change in total energy to stop min. iteration during initial level" << endl);
     352    DoLog(0) && (Log() << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl);
     353    DoLog(0) && (Log() << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl);
    343354//    Log() << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
    344     Log() << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl;
    345     Log() << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl;
    346     Log() << Verbose(0) << " m - Factor by which grid nodes increase between standard and upper level" << endl;
    347     Log() << Verbose(0) << " n - 0: Don't use RiemannTensor, 1: Do" << endl;
    348     Log() << Verbose(0) << " o - Factor by which grid nodes increase between Riemann and standard(?) level" << endl;
    349     Log() << Verbose(0) << " p - Number of Riemann levels" << endl;
    350     Log() << Verbose(0) << " r - 0: Don't Use RiemannTensor, 1: Do" << endl;
    351     Log() << Verbose(0) << " s - 0: Doubly occupied orbitals, 1: Up-/Down-Orbitals" << endl;
    352     Log() << Verbose(0) << " t - Number of orbitals (depends pn SpinType)" << endl;
    353     Log() << Verbose(0) << " u - Number of SpinUp orbitals (depends on SpinType)" << endl;
    354     Log() << Verbose(0) << " v - Number of SpinDown orbitals (depends on SpinType)" << endl;
    355     Log() << Verbose(0) << " w - Number of additional, unoccupied orbitals" << endl;
    356     Log() << Verbose(0) << " x - radial cutoff for ewald summation in Bohrradii" << endl;
    357     Log() << Verbose(0) << " y - 0: Don't do structure optimization beforehand, 1: Do" << endl;
    358     Log() << Verbose(0) << " z - 0: Units are in Bohr radii, 1: units are in Aengstrom" << endl;
    359     Log() << Verbose(0) << " i - 0: Coordinates given in file are absolute, 1: ... are relative to unit cell" << endl;
    360     Log() << Verbose(0) << "=========================================================" << endl;
    361     Log() << Verbose(0) << "INPUT: ";
     355    DoLog(0) && (Log() << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl);
     356    DoLog(0) && (Log() << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl);
     357    DoLog(0) && (Log() << Verbose(0) << " m - Factor by which grid nodes increase between standard and upper level" << endl);
     358    DoLog(0) && (Log() << Verbose(0) << " n - 0: Don't use RiemannTensor, 1: Do" << endl);
     359    DoLog(0) && (Log() << Verbose(0) << " o - Factor by which grid nodes increase between Riemann and standard(?) level" << endl);
     360    DoLog(0) && (Log() << Verbose(0) << " p - Number of Riemann levels" << endl);
     361    DoLog(0) && (Log() << Verbose(0) << " r - 0: Don't Use RiemannTensor, 1: Do" << endl);
     362    DoLog(0) && (Log() << Verbose(0) << " s - 0: Doubly occupied orbitals, 1: Up-/Down-Orbitals" << endl);
     363    DoLog(0) && (Log() << Verbose(0) << " t - Number of orbitals (depends pn SpinType)" << endl);
     364    DoLog(0) && (Log() << Verbose(0) << " u - Number of SpinUp orbitals (depends on SpinType)" << endl);
     365    DoLog(0) && (Log() << Verbose(0) << " v - Number of SpinDown orbitals (depends on SpinType)" << endl);
     366    DoLog(0) && (Log() << Verbose(0) << " w - Number of additional, unoccupied orbitals" << endl);
     367    DoLog(0) && (Log() << Verbose(0) << " x - radial cutoff for ewald summation in Bohrradii" << endl);
     368    DoLog(0) && (Log() << Verbose(0) << " y - 0: Don't do structure optimization beforehand, 1: Do" << endl);
     369    DoLog(0) && (Log() << Verbose(0) << " z - 0: Units are in Bohr radii, 1: units are in Aengstrom" << endl);
     370    DoLog(0) && (Log() << Verbose(0) << " i - 0: Coordinates given in file are absolute, 1: ... are relative to unit cell" << endl);
     371    DoLog(0) && (Log() << Verbose(0) << "=========================================================" << endl);
     372    DoLog(0) && (Log() << Verbose(0) << "INPUT: ");
    362373    cin >> choice;
    363374
    364375    switch (choice) {
    365376        case 'A': // mainname
    366           Log() << Verbose(0) << "Old: " << config::mainname << "\t new: ";
     377          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::mainname << "\t new: ");
    367378          cin >> config::mainname;
    368379          break;
    369380        case 'B': // defaultpath
    370           Log() << Verbose(0) << "Old: " << config::defaultpath << "\t new: ";
     381          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::defaultpath << "\t new: ");
    371382          cin >> config::defaultpath;
    372383          break;
    373384        case 'C': // pseudopotpath
    374           Log() << Verbose(0) << "Old: " << config::pseudopotpath << "\t new: ";
     385          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::pseudopotpath << "\t new: ");
    375386          cin >> config::pseudopotpath;
    376387          break;
    377388
    378389        case 'D': // ProcPEGamma
    379           Log() << Verbose(0) << "Old: " << config::ProcPEGamma << "\t new: ";
     390          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::ProcPEGamma << "\t new: ");
    380391          cin >> config::ProcPEGamma;
    381392          break;
    382393        case 'E': // ProcPEPsi
    383           Log() << Verbose(0) << "Old: " << config::ProcPEPsi << "\t new: ";
     394          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::ProcPEPsi << "\t new: ");
    384395          cin >> config::ProcPEPsi;
    385396          break;
    386397        case 'F': // DoOutVis
    387           Log() << Verbose(0) << "Old: " << config::DoOutVis << "\t new: ";
     398          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoOutVis << "\t new: ");
    388399          cin >> config::DoOutVis;
    389400          break;
    390401        case 'G': // DoOutMes
    391           Log() << Verbose(0) << "Old: " << config::DoOutMes << "\t new: ";
     402          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoOutMes << "\t new: ");
    392403          cin >> config::DoOutMes;
    393404          break;
    394405        case 'H': // DoOutOrbitals
    395           Log() << Verbose(0) << "Old: " << config::DoOutOrbitals << "\t new: ";
     406          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoOutOrbitals << "\t new: ");
    396407          cin >> config::DoOutOrbitals;
    397408          break;
    398409        case 'I': // DoOutCurrent
    399           Log() << Verbose(0) << "Old: " << config::DoOutCurrent << "\t new: ";
     410          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoOutCurrent << "\t new: ");
    400411          cin >> config::DoOutCurrent;
    401412          break;
    402413        case 'J': // DoFullCurrent
    403           Log() << Verbose(0) << "Old: " << config::DoFullCurrent << "\t new: ";
     414          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoFullCurrent << "\t new: ");
    404415          cin >> config::DoFullCurrent;
    405416          break;
    406417        case 'K': // DoPerturbation
    407           Log() << Verbose(0) << "Old: " << config::DoPerturbation << "\t new: ";
     418          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::DoPerturbation << "\t new: ");
    408419          cin >> config::DoPerturbation;
    409420          break;
    410421        case 'L': // CommonWannier
    411           Log() << Verbose(0) << "Old: " << config::CommonWannier << "\t new: ";
     422          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::CommonWannier << "\t new: ");
    412423          cin >> config::CommonWannier;
    413424          break;
    414425        case 'M': // SawtoothStart
    415           Log() << Verbose(0) << "Old: " << config::SawtoothStart << "\t new: ";
     426          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::SawtoothStart << "\t new: ");
    416427          cin >> config::SawtoothStart;
    417428          break;
    418429        case 'N': // VectorPlane
    419           Log() << Verbose(0) << "Old: " << config::VectorPlane << "\t new: ";
     430          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::VectorPlane << "\t new: ");
    420431          cin >> config::VectorPlane;
    421432          break;
    422433        case 'O': // VectorCut
    423           Log() << Verbose(0) << "Old: " << config::VectorCut << "\t new: ";
     434          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::VectorCut << "\t new: ");
    424435          cin >> config::VectorCut;
    425436          break;
    426437        case 'P': // UseAddGramSch
    427           Log() << Verbose(0) << "Old: " << config::UseAddGramSch << "\t new: ";
     438          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::UseAddGramSch << "\t new: ");
    428439          cin >> config::UseAddGramSch;
    429440          break;
    430441        case 'Q': // Seed
    431           Log() << Verbose(0) << "Old: " << config::Seed << "\t new: ";
     442          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::Seed << "\t new: ");
    432443          cin >> config::Seed;
    433444          break;
    434445
    435446        case 'R': // MaxOuterStep
    436           Log() << Verbose(0) << "Old: " << config::MaxOuterStep << "\t new: ";
     447          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxOuterStep << "\t new: ");
    437448          cin >> config::MaxOuterStep;
    438449          break;
    439450        case 'T': // OutVisStep
    440           Log() << Verbose(0) << "Old: " << config::OutVisStep << "\t new: ";
     451          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::OutVisStep << "\t new: ");
    441452          cin >> config::OutVisStep;
    442453          break;
    443454        case 'U': // OutSrcStep
    444           Log() << Verbose(0) << "Old: " << config::OutSrcStep << "\t new: ";
     455          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::OutSrcStep << "\t new: ");
    445456          cin >> config::OutSrcStep;
    446457          break;
    447458        case 'X': // MaxPsiStep
    448           Log() << Verbose(0) << "Old: " << config::MaxPsiStep << "\t new: ";
     459          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxPsiStep << "\t new: ");
    449460          cin >> config::MaxPsiStep;
    450461          break;
    451462        case 'Y': // EpsWannier
    452           Log() << Verbose(0) << "Old: " << config::EpsWannier << "\t new: ";
     463          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::EpsWannier << "\t new: ");
    453464          cin >> config::EpsWannier;
    454465          break;
    455466
    456467        case 'Z': // MaxMinStep
    457           Log() << Verbose(0) << "Old: " << config::MaxMinStep << "\t new: ";
     468          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxMinStep << "\t new: ");
    458469          cin >> config::MaxMinStep;
    459470          break;
    460471        case 'a': // RelEpsTotalEnergy
    461           Log() << Verbose(0) << "Old: " << config::RelEpsTotalEnergy << "\t new: ";
     472          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RelEpsTotalEnergy << "\t new: ");
    462473          cin >> config::RelEpsTotalEnergy;
    463474          break;
    464475        case 'b': // RelEpsKineticEnergy
    465           Log() << Verbose(0) << "Old: " << config::RelEpsKineticEnergy << "\t new: ";
     476          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RelEpsKineticEnergy << "\t new: ");
    466477          cin >> config::RelEpsKineticEnergy;
    467478          break;
    468479        case 'c': // MaxMinStopStep
    469           Log() << Verbose(0) << "Old: " << config::MaxMinStopStep << "\t new: ";
     480          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxMinStopStep << "\t new: ");
    470481          cin >> config::MaxMinStopStep;
    471482          break;
    472483        case 'e': // MaxInitMinStep
    473           Log() << Verbose(0) << "Old: " << config::MaxInitMinStep << "\t new: ";
     484          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxInitMinStep << "\t new: ");
    474485          cin >> config::MaxInitMinStep;
    475486          break;
    476487        case 'f': // InitRelEpsTotalEnergy
    477           Log() << Verbose(0) << "Old: " << config::InitRelEpsTotalEnergy << "\t new: ";
     488          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::InitRelEpsTotalEnergy << "\t new: ");
    478489          cin >> config::InitRelEpsTotalEnergy;
    479490          break;
    480491        case 'g': // InitRelEpsKineticEnergy
    481           Log() << Verbose(0) << "Old: " << config::InitRelEpsKineticEnergy << "\t new: ";
     492          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::InitRelEpsKineticEnergy << "\t new: ");
    482493          cin >> config::InitRelEpsKineticEnergy;
    483494          break;
    484495        case 'h': // InitMaxMinStopStep
    485           Log() << Verbose(0) << "Old: " << config::InitMaxMinStopStep << "\t new: ";
     496          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::InitMaxMinStopStep << "\t new: ");
    486497          cin >> config::InitMaxMinStopStep;
    487498          break;
    … …  
    489500//        case 'j': // BoxLength
    490501//          Log() << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
     502//          double * const cell_size = World::get()->cell_size;
    491503//          for (int i=0;i<6;i++) {
    492504//            Log() << Verbose(0) << "Cell size" << i << ": ";
    493 //            cin >> mol->cell_size[i];
     505//            cin >> cell_size[i];
    494506//          }
    495507//          break;
    496508
    497509        case 'k': // ECut
    498           Log() << Verbose(0) << "Old: " << config::ECut << "\t new: ";
     510          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::ECut << "\t new: ");
    499511          cin >> config::ECut;
    500512          break;
    501513        case 'l': // MaxLevel
    502           Log() << Verbose(0) << "Old: " << config::MaxLevel << "\t new: ";
     514          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxLevel << "\t new: ");
    503515          cin >> config::MaxLevel;
    504516          break;
    505517        case 'm': // RiemannTensor
    506           Log() << Verbose(0) << "Old: " << config::RiemannTensor << "\t new: ";
     518          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RiemannTensor << "\t new: ");
    507519          cin >> config::RiemannTensor;
    508520          break;
    509521        case 'n': // LevRFactor
    510           Log() << Verbose(0) << "Old: " << config::LevRFactor << "\t new: ";
     522          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::LevRFactor << "\t new: ");
    511523          cin >> config::LevRFactor;
    512524          break;
    513525        case 'o': // RiemannLevel
    514           Log() << Verbose(0) << "Old: " << config::RiemannLevel << "\t new: ";
     526          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RiemannLevel << "\t new: ");
    515527          cin >> config::RiemannLevel;
    516528          break;
    517529        case 'p': // Lev0Factor
    518           Log() << Verbose(0) << "Old: " << config::Lev0Factor << "\t new: ";
     530          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::Lev0Factor << "\t new: ");
    519531          cin >> config::Lev0Factor;
    520532          break;
    521533        case 'r': // RTActualUse
    522           Log() << Verbose(0) << "Old: " << config::RTActualUse << "\t new: ";
     534          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RTActualUse << "\t new: ");
    523535          cin >> config::RTActualUse;
    524536          break;
    525537        case 's': // PsiType
    526           Log() << Verbose(0) << "Old: " << config::PsiType << "\t new: ";
     538          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::PsiType << "\t new: ");
    527539          cin >> config::PsiType;
    528540          break;
    529541        case 't': // MaxPsiDouble
    530           Log() << Verbose(0) << "Old: " << config::MaxPsiDouble << "\t new: ";
     542          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::MaxPsiDouble << "\t new: ");
    531543          cin >> config::MaxPsiDouble;
    532544          break;
    533545        case 'u': // PsiMaxNoUp
    534           Log() << Verbose(0) << "Old: " << config::PsiMaxNoUp << "\t new: ";
     546          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::PsiMaxNoUp << "\t new: ");
    535547          cin >> config::PsiMaxNoUp;
    536548          break;
    537549        case 'v': // PsiMaxNoDown
    538           Log() << Verbose(0) << "Old: " << config::PsiMaxNoDown << "\t new: ";
     550          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::PsiMaxNoDown << "\t new: ");
    539551          cin >> config::PsiMaxNoDown;
    540552          break;
    541553        case 'w': // AddPsis
    542           Log() << Verbose(0) << "Old: " << config::AddPsis << "\t new: ";
     554          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::AddPsis << "\t new: ");
    543555          cin >> config::AddPsis;
    544556          break;
    545557
    546558        case 'x': // RCut
    547           Log() << Verbose(0) << "Old: " << config::RCut << "\t new: ";
     559          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RCut << "\t new: ");
    548560          cin >> config::RCut;
    549561          break;
    550562        case 'y': // StructOpt
    551           Log() << Verbose(0) << "Old: " << config::StructOpt << "\t new: ";
     563          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::StructOpt << "\t new: ");
    552564          cin >> config::StructOpt;
    553565          break;
    554566        case 'z': // IsAngstroem
    555           Log() << Verbose(0) << "Old: " << config::IsAngstroem << "\t new: ";
     567          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::IsAngstroem << "\t new: ");
    556568          cin >> config::IsAngstroem;
    557569          break;
    558570        case 'i': // RelativeCoord
    559           Log() << Verbose(0) << "Old: " << config::RelativeCoord << "\t new: ";
     571          DoLog(0) && (Log() << Verbose(0) << "Old: " << config::RelativeCoord << "\t new: ");
    560572          cin >> config::RelativeCoord;
    561573          break;
    … …  
    636648  }
    637649  strcpy(configname, ptr);
    638   Log() << Verbose(0) << "Found configpath: " << configpath << ", dir slash was found at " << last << ", config name is " << configname << "." << endl;
     650  DoLog(0) && (Log() << Verbose(0) << "Found configpath: " << configpath << ", dir slash was found at " << last << ", config name is " << configname << "." << endl);
    639651  delete[](buffer);
    640652};
    … …  
    647659{
    648660  if (FileBuffer != NULL) {
    649     eLog() << Verbose(1) << "WARNING: deleting present FileBuffer in PrepareFileBuffer()." << endl;
     661    DoeLog(2) && (eLog()<< Verbose(2) << "deleting present FileBuffer in PrepareFileBuffer()." << endl);
    650662    delete(FileBuffer);
    651663  }
    … …  
    673685
    674686  if (mol == NULL) {
    675     eLog() << Verbose(0) << "Molecule is not allocated in LoadMolecule(), exit.";
     687    DoeLog(0) && (eLog()<< Verbose(0) << "Molecule is not allocated in LoadMolecule(), exit.");
    676688    performCriticalExit();
    677689  }
    … …  
    679691  ParseForParameter(verbose,FileBuffer,"MaxTypes", 0, 1, 1, int_type, &(MaxTypes), 1, critical);
    680692  if (MaxTypes == 0) {
    681     eLog() << Verbose(0) << "There are no atoms according to MaxTypes in this config file." << endl;
     693    DoeLog(1) && (eLog()<< Verbose(1) << "There are no atoms according to MaxTypes in this config file." << endl);
     694    //performCriticalExit();
    682695  } else {
    683696    // prescan number of ions per type
    684     Log() << Verbose(0) << "Prescanning ions per type: " << endl;
     697    DoLog(0) && (Log() << Verbose(0) << "Prescanning ions per type: " << endl);
    685698    int NoAtoms = 0;
    686699    for (int i=0; i < MaxTypes; i++) {
    … …  
    689702      ParseForParameter(verbose,FileBuffer, name, 0, 2, 1, int_type, &Z, 1, critical);
    690703      elementhash[i] = periode->FindElement(Z);
    691       Log() << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl;
     704      DoLog(1) && (Log() << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl);
    692705      NoAtoms += No[i];
    693706    }
    … …  
    697710    sprintf(name,"Ion_Type%i",MaxTypes);
    698711    if (!ParseForParameter(verbose,FileBuffer, (const char*)name, 1, 1, 1, int_type, &value[0], 1, critical)) {
    699       eLog() << Verbose(0) << "There are no atoms in the config file!" << endl;
     712      DoeLog(0) && (eLog()<< Verbose(0) << "There are no atoms in the config file!" << endl);
     713      performCriticalExit();
    700714      return;
    701715    }
    … …  
    714728      bool status = true;
    715729      while (status) {
    716         Log() << Verbose(0) << "Currently parsing MD step " << repetition << "." << endl;
     730        DoLog(0) && (Log() << Verbose(0) << "Currently parsing MD step " << repetition << "." << endl);
    717731        for (int i=0; i < MaxTypes; i++) {
    718732          sprintf(name,"Ion_Type%i",i+1);
    … …  
    780794      }
    781795      repetition--;
    782       Log() << Verbose(0) << "Found " << repetition << " trajectory steps." << endl;
     796      DoLog(0) && (Log() << Verbose(0) << "Found " << repetition << " trajectory steps." << endl);
    783797      if (repetition <= 1)  // if onyl one step, desactivate use of trajectories
    784798        mol->MDSteps = 0;
    … …  
    792806              ParseForParameter(verbose,FileBuffer, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))
    793807        repetition++;
    794       Log() << Verbose(0) << "I found " << repetition << " times the keyword Ion_Type1_1." << endl;
     808      DoLog(0) && (Log() << Verbose(0) << "I found " << repetition << " times the keyword Ion_Type1_1." << endl);
    795809      // parse in molecule coordinates
    796810      for (int i=0; i < MaxTypes; i++) {
    … …  
    841855  ifstream *file = new ifstream(filename);
    842856  if (file == NULL) {
    843     eLog() << Verbose(0) << "ERROR: config file " << filename << " missing!" << endl;
     857    DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
    844858    return;
    845859  }
    … …  
    953967  // Unit cell and magnetic field
    954968  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    955   mol->cell_size[0] = BoxLength[0];
    956   mol->cell_size[1] = BoxLength[3];
    957   mol->cell_size[2] = BoxLength[4];
    958   mol->cell_size[3] = BoxLength[6];
    959   mol->cell_size[4] = BoxLength[7];
    960   mol->cell_size[5] = BoxLength[8];
     969  double * const cell_size = World::get()->cell_size;
     970  cell_size[0] = BoxLength[0];
     971  cell_size[1] = BoxLength[3];
     972  cell_size[2] = BoxLength[4];
     973  cell_size[3] = BoxLength[6];
     974  cell_size[4] = BoxLength[7];
     975  cell_size[5] = BoxLength[8];
    961976  //if (1) fprintf(stderr,"\n");
    962977
    … …  
    10451060
    10461061  // 2. parse the bond graph file if given
    1047   BG = new BondGraph(IsAngstroem);
    1048   if (BG->LoadBondLengthTable(BondGraphFileName)) {
    1049     Log() << Verbose(0) << "Bond length table loaded successfully." << endl;
    1050   } else {
    1051     Log() << Verbose(0) << "Bond length table loading failed." << endl;
     1062  if (BG == NULL) {
     1063    BG = new BondGraph(IsAngstroem);
     1064    if (BG->LoadBondLengthTable(BondGraphFileName)) {
     1065      DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl);
     1066    } else {
     1067      DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
     1068    }
    10521069  }
    10531070
    10541071  // 3. parse the molecule in
    10551072  LoadMolecule(mol, FileBuffer, periode, FastParsing);
     1073  mol->SetNameFromFilename(filename);
    10561074  mol->ActiveFlag = true;
    10571075  MolList->insert(mol);
    10581076
    10591077  // 4. dissect the molecule into connected subgraphs
    1060   BG->ConstructBondGraph(mol);
     1078  // don't do this here ...
     1079  //MolList->DissectMoleculeIntoConnectedSubgraphs(mol,this);
     1080  //delete(mol);
    10611081
    10621082  delete(FileBuffer);
    … …  
    10741094  ifstream *file = new ifstream(filename);
    10751095  if (file == NULL) {
    1076     eLog() << Verbose(0) << "ERROR: config file " << filename << " missing!" << endl;
     1096    DoeLog(1) && (eLog()<< Verbose(1) << "config file " << filename << " missing!" << endl);
    10771097    return;
    10781098  }
    … …  
    11521172
    11531173  ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    1154   mol->cell_size[0] = BoxLength[0];
    1155   mol->cell_size[1] = BoxLength[3];
    1156   mol->cell_size[2] = BoxLength[4];
    1157   mol->cell_size[3] = BoxLength[6];
    1158   mol->cell_size[4] = BoxLength[7];
    1159   mol->cell_size[5] = BoxLength[8];
     1174  double * const cell_size = World::get()->cell_size;
     1175  cell_size[0] = BoxLength[0];
     1176  cell_size[1] = BoxLength[3];
     1177  cell_size[2] = BoxLength[4];
     1178  cell_size[3] = BoxLength[6];
     1179  cell_size[4] = BoxLength[7];
     1180  cell_size[5] = BoxLength[8];
    11601181  if (1) fprintf(stderr,"\n");
    11611182  config::DoPerturbation = 0;
    … …  
    12331254  BG = new BondGraph(IsAngstroem);
    12341255  if (BG->LoadBondLengthTable(BondGraphFileName)) {
    1235     Log() << Verbose(0) << "Bond length table loaded successfully." << endl;
     1256    DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl);
    12361257  } else {
    1237     Log() << Verbose(0) << "Bond length table loading failed." << endl;
     1258    DoLog(0) && (Log() << Verbose(0) << "Bond length table loading failed." << endl);
    12381259  }
    12391260
    … …  
    12421263  for (i=MAX_ELEMENTS;i--;)
    12431264    elementhash[i] = NULL;
    1244   Log() << Verbose(0) << "Parsing Ions ..." << endl;
     1265  DoLog(0) && (Log() << Verbose(0) << "Parsing Ions ..." << endl);
    12451266  No=0;
    12461267  found = 0;
    12471268  while (getline(*file,zeile,'\n')) {
    12481269    if (zeile.find("Ions_Data") == 0) {
    1249       Log() << Verbose(1) << "found Ions_Data...begin parsing" << endl;
     1270      DoLog(1) && (Log() << Verbose(1) << "found Ions_Data...begin parsing" << endl);
    12501271      found ++;
    12511272    }
    … …  
    12611282      input >> b;     // element mass
    12621283      elementhash[No] = periode->FindElement(Z);
    1263       Log() << Verbose(1) << "AtomNo: " << AtomNo << "\tZ: " << Z << "\ta:" << a << "\tl:"  << l << "\b:" << b << "\tElement:" << elementhash[No] << "\t:" << endl;
     1284      DoLog(1) && (Log() << Verbose(1) << "AtomNo: " << AtomNo << "\tZ: " << Z << "\ta:" << a << "\tl:"  << l << "\b:" << b << "\tElement:" << elementhash[No] << "\t:" << endl);
    12641285      for(i=0;i<AtomNo;i++) {
    12651286        if (!getline(*file,zeile,'\n')) {// parse on and on
    1266           Log() << Verbose(2) << "Error: Too few items in ionic list of element" << elementhash[No] << "." << endl << "Exiting." << endl;
     1287          DoLog(2) && (Log() << Verbose(2) << "Error: Too few items in ionic list of element" << elementhash[No] << "." << endl << "Exiting." << endl);
    12671288          // return 1;
    12681289        } else {
    … …  
    12951316  // bring MaxTypes up to date
    12961317  mol->CountElements();
     1318  const double * const cell_size = World::get()->cell_size;
    12971319  ofstream * const output = new ofstream(filename, ios::out);
    12981320  if (output != NULL) {
    … …  
    13651387    *output << endl;
    13661388    *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
    1367     *output << mol->cell_size[0] << "\t" << endl;
    1368     *output << mol->cell_size[1] << "\t" << mol->cell_size[2] << "\t" << endl;
    1369     *output << mol->cell_size[3] << "\t" << mol->cell_size[4] << "\t" << mol->cell_size[5] << "\t" << endl;
     1389    *output << cell_size[0] << "\t" << endl;
     1390    *output << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
     1391    *output << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
    13701392    // FIXME
    13711393    *output << endl;
    … …  
    14101432    delete(output);
    14111433    return result;
    1412   } else
     1434  } else {
     1435    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open output file:" << filename << endl);
    14131436    return false;
     1437  }
    14141438};
    14151439
    … …  
    14301454    *fname << filename << ".in";
    14311455    output = new ofstream(fname->str().c_str(), ios::out);
     1456    if (output == NULL) {
     1457      DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc output file:" << fname << endl);
     1458      delete(fname);
     1459      return false;
     1460    }
    14321461    *output << "% Created by MoleCuilder" << endl;
    14331462    *output << "mpqc: (" << endl;
    … …  
    14681497    *fname << filename << ".hess.in";
    14691498    output = new ofstream(fname->str().c_str(), ios::out);
     1499    if (output == NULL) {
     1500      DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open mpqc hessian output file:" << fname << endl);
     1501      delete(fname);
     1502      return false;
     1503    }
    14701504    *output << "% Created by MoleCuilder" << endl;
    14711505    *output << "mpqc: (" << endl;
    … …  
    14991533    delete(fname);
    15001534  }
     1535
     1536  return true;
     1537};
     1538
     1539/** Stores all atoms from all molecules in a PDB input file.
     1540 * Note that this format cannot be parsed again.
     1541 * \param *filename name of file (without ".in" suffix!)
     1542 * \param *MolList pointer to MoleculeListClass containing all atoms
     1543 */
     1544bool config::SavePDB(const char * const filename, const MoleculeListClass * const MolList) const
     1545{
     1546  int AtomNo = -1;
     1547  int MolNo = 0;
     1548  atom *Walker = NULL;
     1549  FILE *f = NULL;
     1550
     1551  char name[MAXSTRINGSIZE];
     1552  strncpy(name, filename, MAXSTRINGSIZE-1);
     1553  strncat(name, ".pdb", MAXSTRINGSIZE-(strlen(name)+1));
     1554  f = fopen(name, "w" );
     1555  if (f == NULL) {
     1556    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open pdb output file:" << name << endl);
     1557    return false;
     1558  }
     1559  fprintf(f, "# Created by MoleCuilder\n");
     1560
     1561  for (MoleculeList::const_iterator Runner = MolList->ListOfMolecules.begin(); Runner != MolList->ListOfMolecules.end(); Runner++) {
     1562    Walker = (*Runner)->start;
     1563    int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo");
     1564    AtomNo = 0;
     1565    while (Walker->next != (*Runner)->end) {
     1566      Walker = Walker->next;
     1567      sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
     1568      elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
     1569      fprintf(f,
     1570             "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
     1571             Walker->nr,                /* atom serial number */
     1572             name,         /* atom name */
     1573             (*Runner)->name,      /* residue name */
     1574             'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
     1575             MolNo,         /* residue sequence number */
     1576             Walker->node->x[0],                 /* position X in Angstroem */
     1577             Walker->node->x[1],                 /* position Y in Angstroem */
     1578             Walker->node->x[2],                 /* position Z in Angstroem */
     1579             (double)Walker->type->Valence,         /* occupancy */
     1580             (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
     1581             "0",            /* segment identifier */
     1582             Walker->type->symbol,    /* element symbol */
     1583             "0");           /* charge */
     1584      AtomNo++;
     1585    }
     1586    Free(&elementNo);
     1587    MolNo++;
     1588  }
     1589  fclose(f);
     1590
     1591  return true;
     1592};
     1593
     1594/** Stores all atoms in a PDB input file.
     1595 * Note that this format cannot be parsed again.
     1596 * \param *filename name of file (without ".in" suffix!)
     1597 * \param *mol pointer to molecule
     1598 */
     1599bool config::SavePDB(const char * const filename, const molecule * const mol) const
     1600{
     1601  int AtomNo = -1;
     1602  atom *Walker = NULL;
     1603  FILE *f = NULL;
     1604
     1605  int *elementNo = Calloc<int>(MAX_ELEMENTS, "config::SavePDB - elementNo");
     1606  char name[MAXSTRINGSIZE];
     1607  strncpy(name, filename, MAXSTRINGSIZE-1);
     1608  strncat(name, ".pdb", MAXSTRINGSIZE-(strlen(name)+1));
     1609  f = fopen(name, "w" );
     1610  if (f == NULL) {
     1611    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open pdb output file:" << name << endl);
     1612    Free(&elementNo);
     1613    return false;
     1614  }
     1615  fprintf(f, "# Created by MoleCuilder\n");
     1616
     1617  Walker = mol->start;
     1618  AtomNo = 0;
     1619  while (Walker->next != mol->end) {
     1620    Walker = Walker->next;
     1621    sprintf(name, "%2s%2d",Walker->type->symbol, elementNo[Walker->type->Z]);
     1622    elementNo[Walker->type->Z] = (elementNo[Walker->type->Z]+1) % 100;   // confine to two digits
     1623    fprintf(f,
     1624           "ATOM %6u %-4s %4s%c%4u    %8.3f%8.3f%8.3f%6.2f%6.2f      %4s%2s%2s\n",
     1625           Walker->nr,                /* atom serial number */
     1626           name,         /* atom name */
     1627           mol->name,      /* residue name */
     1628           'a'+(unsigned char)(AtomNo % 26),           /* letter for chain */
     1629           0,         /* residue sequence number */
     1630           Walker->node->x[0],                 /* position X in Angstroem */
     1631           Walker->node->x[1],                 /* position Y in Angstroem */
     1632           Walker->node->x[2],                 /* position Z in Angstroem */
     1633           (double)Walker->type->Valence,         /* occupancy */
     1634           (double)Walker->type->NoValenceOrbitals,          /* temperature factor */
     1635           "0",            /* segment identifier */
     1636           Walker->type->symbol,    /* element symbol */
     1637           "0");           /* charge */
     1638    AtomNo++;
     1639  }
     1640  fclose(f);
     1641  Free(&elementNo);
     1642
     1643  return true;
     1644};
     1645
     1646/** Stores all atoms in a TREMOLO data input file.
     1647 * Note that this format cannot be parsed again.
     1648 * Note that TREMOLO does not like Id starting at 0, but at 1. Atoms with Id 0 are discarded!
     1649 * \param *filename name of file (without ".in" suffix!)
     1650 * \param *mol pointer to molecule
     1651 */
     1652bool config::SaveTREMOLO(const char * const filename, const molecule * const mol) const
     1653{
     1654  atom *Walker = NULL;
     1655  ofstream *output = NULL;
     1656  stringstream * const fname = new stringstream;
     1657
     1658  *fname << filename << ".data";
     1659  output = new ofstream(fname->str().c_str(), ios::out);
     1660  if (output == NULL) {
     1661    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open tremolo output file:" << fname << endl);
     1662    delete(fname);
     1663    return false;
     1664  }
     1665
     1666  // scan maximum number of neighbours
     1667  Walker = mol->start;
     1668  int MaxNeighbours = 0;
     1669  while (Walker->next != mol->end) {
     1670    Walker = Walker->next;
     1671    const int count = Walker->ListOfBonds.size();
     1672    if (MaxNeighbours < count)
     1673      MaxNeighbours = count;
     1674  }
     1675  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
     1676
     1677  Walker = mol->start;
     1678  while (Walker->next != mol->end) {
     1679    Walker = Walker->next;
     1680    *output << Walker->nr << "\t";
     1681    *output << Walker->Name << "\t";
     1682    *output << mol->name << "\t";
     1683    *output << 0 << "\t";
     1684    *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
     1685    *output << (double)Walker->type->Valence << "\t";
     1686    *output << Walker->type->symbol << "\t";
     1687    for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
     1688      *output << (*runner)->GetOtherAtom(Walker)->nr << "\t";
     1689    for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
     1690      *output << "-\t";
     1691    *output << endl;
     1692  }
     1693  output->flush();
     1694  output->close();
     1695  delete(output);
     1696  delete(fname);
     1697
     1698  return true;
     1699};
     1700
     1701/** Stores all atoms from all molecules in a TREMOLO data input file.
     1702 * Note that this format cannot be parsed again.
     1703 * Note that TREMOLO does not like Id starting at 0, but at 1. Atoms with Id 0 are discarded!
     1704 * \param *filename name of file (without ".in" suffix!)
     1705 * \param *MolList pointer to MoleculeListClass containing all atoms
     1706 */
     1707bool config::SaveTREMOLO(const char * const filename, const MoleculeListClass * const MolList) const
     1708{
     1709  atom *Walker = NULL;
     1710  ofstream *output = NULL;
     1711  stringstream * const fname = new stringstream;
     1712
     1713  *fname << filename << ".data";
     1714  output = new ofstream(fname->str().c_str(), ios::out);
     1715  if (output == NULL) {
     1716    DoeLog(1) && (eLog()<< Verbose(1) << "Cannot open tremolo output file:" << fname << endl);
     1717    delete(fname);
     1718    return false;
     1719  }
     1720
     1721  // scan maximum number of neighbours
     1722  int MaxNeighbours = 0;
     1723  for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
     1724    Walker = (*MolWalker)->start;
     1725    while (Walker->next != (*MolWalker)->end) {
     1726      Walker = Walker->next;
     1727      const int count = Walker->ListOfBonds.size();
     1728      if (MaxNeighbours < count)
     1729        MaxNeighbours = count;
     1730    }
     1731  }
     1732  *output << "# ATOMDATA Id name resName resSeq x=3 charge type neighbors=" << MaxNeighbours << endl;
     1733
     1734  // create global to local id map
     1735  int **LocalNotoGlobalNoMap = Calloc<int *>(MolList->ListOfMolecules.size(), "config::SaveTREMOLO - **LocalNotoGlobalNoMap");
     1736  {
     1737    int MolCounter = 0;
     1738    int AtomNo = 0;
     1739    for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
     1740      LocalNotoGlobalNoMap[MolCounter] = Calloc<int>(MolList->CountAllAtoms(), "config::SaveTREMOLO - *LocalNotoGlobalNoMap[]");
     1741
     1742      (*MolWalker)->SetIndexedArrayForEachAtomTo( LocalNotoGlobalNoMap[MolCounter], &atom::nr, IncrementalAbsoluteValue, &AtomNo);
     1743
     1744      MolCounter++;
     1745    }
     1746  }
     1747
     1748  // write the file
     1749  {
     1750    int MolCounter = 0;
     1751    int AtomNo = 0;
     1752    for (MoleculeList::const_iterator MolWalker = MolList->ListOfMolecules.begin(); MolWalker != MolList->ListOfMolecules.end(); MolWalker++) {
     1753      Walker = (*MolWalker)->start;
     1754      while (Walker->next != (*MolWalker)->end) {
     1755        Walker = Walker->next;
     1756        *output << AtomNo+1 << "\t";
     1757        *output << Walker->Name << "\t";
     1758        *output << (*MolWalker)->name << "\t";
     1759        *output << MolCounter+1 << "\t";
     1760        *output << Walker->node->x[0] << "\t" << Walker->node->x[1] << "\t" << Walker->node->x[2] << "\t";
     1761        *output << (double)Walker->type->Valence << "\t";
     1762        *output << Walker->type->symbol << "\t";
     1763        for (BondList::iterator runner = Walker->ListOfBonds.begin(); runner != Walker->ListOfBonds.end(); runner++)
     1764          *output << LocalNotoGlobalNoMap[MolCounter][ (*runner)->GetOtherAtom(Walker)->nr ]+1 << "\t";
     1765        for(int i=Walker->ListOfBonds.size(); i < MaxNeighbours; i++)
     1766          *output << "-\t";
     1767        *output << endl;
     1768        AtomNo++;
     1769      }
     1770      MolCounter++;
     1771    }
     1772  }
     1773
     1774  // store & free
     1775  output->flush();
     1776  output->close();
     1777  delete(output);
     1778  delete(fname);
     1779  for(size_t i=0;i<MolList->ListOfMolecules.size(); i++)
     1780    Free(&LocalNotoGlobalNoMap[i]);
     1781  Free(&LocalNotoGlobalNoMap);
    15011782
    15021783  return true;
    … …  
    18502131              }
    18512132              line++;
    1852             } while (dummy1 != NULL && (dummy1[0] == '#') || (dummy1[0] == '\n'));
     2133            } while ((dummy1 != NULL) && ((dummy1[0] == '#') || (dummy1[0] == '\n')));
    18532134            dummy = dummy1;
    18542135          } else { // simple int, strings or doubles start in the same line
Note: See TracChangeset for help on using the changeset viewer.