Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/builder.cpp

    r9907e7 rf65e1f  
    5353
    5454#include <cstring>
    55 #include <cstdlib>
    5655
    5756#include "analysis_bonds.hpp"
     
    865864
    866865        mol->CountAtoms(); // recount atoms
    867         if (mol->getAtomCount() != 0) {  // if there is more than none
    868           count = mol->getAtomCount();  // is changed becausing of adding, thus has to be stored away beforehand
     866        if (mol->AtomCount != 0) {  // if there is more than none
     867          count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
    869868          Elements = new element *[count];
    870869          vectors = new Vector *[count];
     
    12961295  // generate some KeySets
    12971296  DoLog(0) && (Log() << Verbose(0) << "Generating KeySets." << endl);
    1298   KeySet TestSets[mol->getAtomCount()+1];
     1297  KeySet TestSets[mol->AtomCount+1];
    12991298  i=1;
    13001299  while (Walker->next != mol->end) {
     
    13071306  DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl);
    13081307  KeySetTestPair test;
    1309   test = TestSets[mol->getAtomCount()-1].insert(Walker->nr);
     1308  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
    13101309  if (test.second) {
    13111310    DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl);
     
    13131312    DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl);
    13141313  }
    1315   TestSets[mol->getAtomCount()].insert(mol->end->previous->nr);
    1316   TestSets[mol->getAtomCount()].insert(mol->end->previous->previous->previous->nr);
     1314  TestSets[mol->AtomCount].insert(mol->end->previous->nr);
     1315  TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
    13171316
    13181317  // constructing Graph structure
     
    13221321  // insert KeySets into Subgraphs
    13231322  DoLog(0) && (Log() << Verbose(0) << "Inserting KeySets into Subgraph class." << endl);
    1324   for (int j=0;j<mol->getAtomCount();j++) {
     1323  for (int j=0;j<mol->AtomCount;j++) {
    13251324    Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
    13261325  }
    13271326  DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl);
    13281327  GraphTestPair test2;
    1329   test2 = Subgraphs.insert(GraphPair (TestSets[mol->getAtomCount()],pair<int, double>(counter++, 1.)));
     1328  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
    13301329  if (test2.second) {
    13311330    DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl);
     
    14951494 */
    14961495static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode,
    1497                                    config& configuration, char **ConfigFileName, set<int> &ArgcList)
     1496                                   config& configuration, char *&ConfigFileName, set<int> &ArgcList)
    14981497{
    14991498  Vector x,y,z,n;  // coordinates for absolute point in cell volume
     
    15131512  molecule *mol = NULL;
    15141513  string BondGraphFileName("\n");
    1515   bool DatabasePathGiven = false;
     1514  strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1);
    15161515
    15171516  if (argc > 1) { // config file specified as option
     
    15301529            break;
    15311530          case 'v':
    1532             setVerbosity(atoi(argv[argptr]));
     1531            ArgcList.insert(argptr-1);
     1532            return(1);
     1533            break;
     1534          case 'V':
    15331535            ArgcList.insert(argptr-1);
    15341536            ArgcList.insert(argptr);
    15351537            argptr++;
    15361538            break;
    1537           case 'V':
    1538             ArgcList.insert(argptr-1);
    1539             return(1);
    1540             break;
    15411539          case 'B':
    1542             ArgcList.insert(argptr-1);
    1543             ArgcList.insert(argptr);
    1544             ArgcList.insert(argptr+1);
    1545             ArgcList.insert(argptr+2);
    1546             ArgcList.insert(argptr+3);
    1547             ArgcList.insert(argptr+4);
    1548             ArgcList.insert(argptr+5);
    1549             argptr+=6;
    15501540            if (ExitFlag == 0) ExitFlag = 1;
     1541            if ((argptr+5 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) ) {
     1542              ExitFlag = 255;
     1543              DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for bounding in box: -B <xx> <xy> <xz> <yy> <yz> <zz>" << endl);
     1544              performCriticalExit();
     1545            } else {
     1546              SaveFlag = true;
     1547              j = -1;
     1548              DoLog(1) && (Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl);
     1549              double * const cell_size = World::getInstance().getDomain();
     1550              for (int i=0;i<6;i++) {
     1551                cell_size[i] = atof(argv[argptr+i]);
     1552              }
     1553              argptr+=6;
     1554            }
    15511555            break;
    15521556          case 'e':
     
    15571561              DoLog(0) && (Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl);
    15581562              strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1);
    1559               DatabasePathGiven = true;
    15601563              argptr+=1;
    15611564            }
     
    15711574            }
    15721575            break;
    1573           case 'M':
    1574             if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    1575               ExitFlag = 255;
    1576               DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -M <basis name>" << endl);
    1577               performCriticalExit();
    1578             } else {
    1579               configuration.basis = argv[argptr];
    1580               DoLog(1) && (Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl);
    1581               argptr+=1;
    1582             }
    1583             break;
    15841576          case 'n':
    15851577            DoLog(0) && (Log() << Verbose(0) << "I won't parse trajectories." << endl);
     
    15891581            {
    15901582              World::getInstance().setDefaultName(argv[argptr]);
    1591               DoLog(0) && (Log() << Verbose(0) << "Default name of new molecules set to " << World::getInstance().getDefaultName() << "." << endl);
     1583              DoLog(0) && (Log() << Verbose(0) << "Default name of new molecules set to " << *World::getInstance().getDefaultName() << "." << endl);
    15921584            }
    15931585            break;
     
    16011593
    16021594    // 3a. Parse the element database
    1603     if (DatabasePathGiven)
    1604       if (periode->LoadPeriodentafel(configuration.databasepath)) {
    1605         DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl);
    1606         //periode->Output();
    1607       } else {
    1608         DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl);
    1609         return 1;
    1610       }
     1595    if (periode->LoadPeriodentafel(configuration.databasepath)) {
     1596      DoLog(0) && (Log() << Verbose(0) << "Element list loaded successfully." << endl);
     1597      //periode->Output();
     1598    } else {
     1599      DoLog(0) && (Log() << Verbose(0) << "Element list loading failed." << endl);
     1600      return 1;
     1601    }
    16111602    // 3b. Find config file name and parse if possible, also BondGraphFileName
    16121603    if (argv[1][0] != '-') {
     
    16221613        } else {
    16231614          DoLog(0) && (Log() << Verbose(0) << "Empty configuration file." << endl);
    1624           strcpy(*ConfigFileName, argv[1]);
     1615          ConfigFileName = argv[1];
    16251616          configPresent = empty;
    16261617          output.close();
     
    16281619      } else {
    16291620        test.close();
    1630         strcpy(*ConfigFileName, argv[1]);
     1621        ConfigFileName = argv[1];
    16311622        DoLog(1) && (Log() << Verbose(1) << "Specified config file found, parsing ... ");
    1632         switch (configuration.TestSyntax(*ConfigFileName, periode)) {
     1623        switch (configuration.TestSyntax(ConfigFileName, periode)) {
    16331624          case 1:
    16341625            DoLog(0) && (Log() << Verbose(0) << "new syntax." << endl);
    1635             configuration.Load(*ConfigFileName, BondGraphFileName, periode, molecules);
     1626            configuration.Load(ConfigFileName, BondGraphFileName, periode, molecules);
    16361627            configPresent = present;
    16371628            break;
    16381629          case 0:
    16391630            DoLog(0) && (Log() << Verbose(0) << "old syntax." << endl);
    1640             configuration.LoadOld(*ConfigFileName, BondGraphFileName, periode, molecules);
     1631            configuration.LoadOld(ConfigFileName, BondGraphFileName, periode, molecules);
    16411632            configPresent = present;
    16421633            break;
     
    16591650       mol = World::getInstance().createMolecule();
    16601651       mol->ActiveFlag = true;
    1661        if (*ConfigFileName != NULL)
    1662          mol->SetNameFromFilename(*ConfigFileName);
     1652       if (ConfigFileName != NULL)
     1653         mol->SetNameFromFilename(ConfigFileName);
    16631654       molecules->insert(mol);
    16641655     }
     
    17141705                if (first->type != NULL) {
    17151706                  mol->AddAtom(first);  // add to molecule
    1716                   if ((configPresent == empty) && (mol->getAtomCount() != 0))
     1707                  if ((configPresent == empty) && (mol->AtomCount != 0))
    17171708                    configPresent = present;
    17181709                } else
     
    17271718        if (configPresent == present) {
    17281719          switch(argv[argptr-1][1]) {
     1720            case 'M':
     1721              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1722                ExitFlag = 255;
     1723                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl);
     1724                performCriticalExit();
     1725              } else {
     1726                configuration.basis = argv[argptr];
     1727                DoLog(1) && (Log() << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl);
     1728                argptr+=1;
     1729              }
     1730              break;
    17291731            case 'D':
    17301732              if (ExitFlag == 0) ExitFlag = 1;
     
    17321734                DoLog(1) && (Log() << Verbose(1) << "Depth-First-Search Analysis." << endl);
    17331735                MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
    1734                 int *MinimumRingSize = new int[mol->getAtomCount()];
     1736                int *MinimumRingSize = new int[mol->AtomCount];
    17351737                atom ***ListOfLocalAtoms = NULL;
    17361738                class StackClass<bond *> *BackEdgeStack = NULL;
     
    17521754                  delete(Subgraphs);
    17531755                  for (int i=0;i<FragmentCounter;i++)
    1754                     delete[](ListOfLocalAtoms[i]);
    1755                   delete[](ListOfLocalAtoms);
     1756                    Free(&ListOfLocalAtoms[i]);
     1757                  Free(&ListOfLocalAtoms);
    17561758                }
    17571759                delete(BackEdgeStack);
     
    17931795                        if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+2])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-')) {
    17941796                          ExitFlag = 255;
    1795                           DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C E <Z1> <Z2> <output> <bin output> <binstart> <binend>" << endl);
     1797                          DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C E <Z1> <Z2> <output> <bin output>" << endl);
    17961798                          performCriticalExit();
    17971799                        } else {
     
    18081810                          else
    18091811                            correlationmap = PairCorrelation(molecules, elemental, elemental2);
    1810                           OutputPairCorrelation(&output, correlationmap);
     1812                          //OutputCorrelationToSurface(&output, correlationmap);
    18111813                          BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
    18121814                          OutputCorrelation ( &binoutput, binmap );
     
    18241826                        if ((argptr+8 >= argc) || (!IsValidNumber(argv[argptr+1])) ||  (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+7])) || (!IsValidNumber(argv[argptr+8])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-') || (argv[argptr+4][0] == '-') || (argv[argptr+5][0] == '-') || (argv[argptr+6][0] == '-')) {
    18251827                          ExitFlag = 255;
    1826                           DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C P <Z1> <x> <y> <z> <output> <bin output> <binstart> <binend>" << endl);
     1828                          DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C P <Z1> <x> <y> <z> <output> <bin output>" << endl);
    18271829                          performCriticalExit();
    18281830                        } else {
     
    18391841                          else
    18401842                            correlationmap = CorrelationToPoint(molecules, elemental, Point);
    1841                           OutputCorrelationToPoint(&output, correlationmap);
     1843                          //OutputCorrelationToSurface(&output, correlationmap);
    18421844                          BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
    18431845                          OutputCorrelation ( &binoutput, binmap );
     
    18801882                          int counter  = 0;
    18811883                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    1882                             if ((Boundary == NULL) || (Boundary->getAtomCount() < (*BigFinder)->getAtomCount())) {
     1884                            if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
    18831885                              Boundary = *BigFinder;
    18841886                            }
    18851887                            counter++;
    18861888                          }
    1887                           bool *Actives = new bool[counter];
     1889                          bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
    18881890                          counter = 0;
    18891891                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
     
    19111913                          output.close();
    19121914                          binoutput.close();
    1913                           counter = 0;
    19141915                          for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
    19151916                            (*BigFinder)->ActiveFlag = Actives[counter++];
    1916                           delete[](Actives);
     1917                          Free(&Actives);
    19171918                          delete(LCList);
    19181919                          delete(TesselStruct);
     
    19401941                performCriticalExit();
    19411942              } else {
    1942                 mol->getAtomCount();
    19431943                SaveFlag = true;
    19441944                DoLog(1) && (Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl);
     
    19501950            case 'F':
    19511951              if (ExitFlag == 0) ExitFlag = 1;
    1952               ArgcList.insert(argptr-1);
    1953               ArgcList.insert(argptr);
    1954               ArgcList.insert(argptr+1);
    1955               ArgcList.insert(argptr+2);
    1956               ArgcList.insert(argptr+3);
    1957               ArgcList.insert(argptr+4);
    1958               ArgcList.insert(argptr+5);
    1959               ArgcList.insert(argptr+6);
    1960               ArgcList.insert(argptr+7);
    1961               ArgcList.insert(argptr+8);
    1962               ArgcList.insert(argptr+9);
    1963               ArgcList.insert(argptr+10);
    1964               ArgcList.insert(argptr+11);
    1965               ArgcList.insert(argptr+12);
    1966               argptr+=13;
     1952              MaxDistance = -1;
     1953              if (argv[argptr-1][2] == 'F') { // option is -FF?
     1954                // fetch first argument as max distance to surface
     1955                MaxDistance = atof(argv[argptr++]);
     1956                DoLog(0) && (Log() << Verbose(0) << "Filling with maximum layer distance of " << MaxDistance << "." << endl);
     1957              }
     1958              if ((argptr+7 >=argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (!IsValidNumber(argv[argptr+7]))) {
     1959                ExitFlag = 255;
     1960                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <xyz of filler> <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl);
     1961                performCriticalExit();
     1962              } else {
     1963                SaveFlag = true;
     1964                DoLog(1) && (Log() << Verbose(1) << "Filling Box with water molecules." << endl);
     1965                // construct water molecule
     1966                molecule *filler = World::getInstance().createMolecule();
     1967                if (!filler->AddXYZFile(argv[argptr])) {
     1968                  DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse filler molecule from " << argv[argptr] << "." << endl);
     1969                }
     1970                filler->SetNameFromFilename(argv[argptr]);
     1971                configuration.BG->ConstructBondGraph(filler);
     1972                molecule *Filling = NULL;
     1973                atom *second = NULL, *third = NULL;
     1974                first = World::getInstance().createAtom();
     1975                first->type = periode->FindElement(1);
     1976                first->x = Vector(0.441, -0.143, 0.);
     1977                filler->AddAtom(first);
     1978                second = World::getInstance().createAtom();
     1979                second->type = periode->FindElement(1);
     1980                second->x = Vector(-0.464, 1.137, 0.0);
     1981                filler->AddAtom(second);
     1982                third = World::getInstance().createAtom();
     1983                third->type = periode->FindElement(8);
     1984                third->x = Vector(-0.464, 0.177, 0.);
     1985                filler->AddAtom(third);
     1986                filler->AddBond(first, third, 1);
     1987                filler->AddBond(second, third, 1);
     1988                // call routine
     1989                double distance[NDIM];
     1990                for (int i=0;i<NDIM;i++)
     1991                  distance[i] = atof(argv[argptr+i+1]);
     1992                Filling = FillBoxWithMolecule(molecules, filler, configuration, MaxDistance, distance, atof(argv[argptr+4]), atof(argv[argptr+5]), atof(argv[argptr+6]), atoi(argv[argptr+7]));
     1993                if (Filling != NULL) {
     1994                  Filling->ActiveFlag = false;
     1995                  molecules->insert(Filling);
     1996                }
     1997                World::getInstance().destroyMolecule(filler);
     1998                argptr+=6;
     1999              }
    19672000              break;
    19682001            case 'A':
     
    19742007              } else {
    19752008                DoLog(0) && (Log() << Verbose(0) << "Parsing bonds from " << argv[argptr] << "." << endl);
    1976                 ifstream input(argv[argptr]);
    1977                 mol->CreateAdjacencyListFromDbondFile(&input);
    1978                 input.close();
     2009                ifstream *input = new ifstream(argv[argptr]);
     2010                mol->CreateAdjacencyListFromDbondFile(input);
     2011                input->close();
    19792012                argptr+=1;
    19802013              }
     
    20262059                int counter  = 0;
    20272060                for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    2028                   if ((Boundary == NULL) || (Boundary->getAtomCount() < (*BigFinder)->getAtomCount())) {
     2061                  (*BigFinder)->CountAtoms();
     2062                  if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
    20292063                    Boundary = *BigFinder;
    20302064                  }
    20312065                  counter++;
    20322066                }
    2033                 DoLog(1) && (Log() << Verbose(1) << "Biggest molecule has " << Boundary->getAtomCount() << " atoms." << endl);
     2067                DoLog(1) && (Log() << Verbose(1) << "Biggest molecule has " << Boundary->AtomCount << " atoms." << endl);
    20342068                start = clock();
    20352069                LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.);
     
    20662100              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
    20672101                ExitFlag = 255;
    2068                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for linear interpolation: -L <step0> <step1> <prefix> <identity mapping?>" << endl);
     2102                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for storing tempature: -L <step0> <step1> <prefix> <identity mapping?>" << endl);
    20692103                performCriticalExit();
    20702104              } else {
     
    20982132            case 'R':
    20992133              if (ExitFlag == 0) ExitFlag = 1;
    2100               if ((argptr+3 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) || (!IsValidNumber(argv[argptr+3])))  {
     2134              if ((argptr+1 >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])))  {
    21012135                ExitFlag = 255;
    2102                 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <x> <y> <z> <distance>" << endl);
     2136                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for removing atoms: -R <id> <distance>" << endl);
    21032137                performCriticalExit();
    21042138              } else {
    21052139                SaveFlag = true;
    2106                 const double radius = atof(argv[argptr+3]);
    2107                 Vector point(atof(argv[argptr]),atof(argv[argptr+1]),atof(argv[argptr+2]));
    2108                 DoLog(1) && (Log() << Verbose(1) << "Removing atoms around " << point << " with radius " << radius << "." << endl);
    2109                 atom *Walker = NULL;
    2110                 molecule::iterator advancer = mol->begin();
    2111                 for(molecule::iterator iter = advancer; advancer != mol->end();) {
    2112                   iter = advancer++;
    2113                   if ((*iter)->x.DistanceSquared(point) > radius*radius){ // distance to first above radius ...
    2114                     Walker = (*iter);
    2115                     DoLog(1) && (Log() << Verbose(1) << "Removing atom " << *Walker << "." << endl);
    2116                     mol->RemoveAtom(*(iter));
    2117                     World::getInstance().destroyAtom(Walker);
     2140                DoLog(1) && (Log() << Verbose(1) << "Removing atoms around " << argv[argptr] << " with radius " << argv[argptr+1] << "." << endl);
     2141                double tmp1 = atof(argv[argptr+1]);
     2142                atom *third = mol->FindAtom(atoi(argv[argptr]));
     2143                atom *first = mol->start;
     2144                if ((third != NULL) && (first != mol->end)) {
     2145                  atom *second = first->next;
     2146                  while(second != mol->end) {
     2147                    first = second;
     2148                    second = first->next;
     2149                    if (first->x.DistanceSquared(third->x) > tmp1*tmp1) // distance to first above radius ...
     2150                      mol->RemoveAtom(first);
    21182151                  }
     2152                } else {
     2153                  DoeLog(1) && (eLog()<< Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl);
    21192154                }
    2120                 argptr+=4;
     2155                argptr+=2;
    21212156              }
    21222157              break;
     
    21552190            case 's':
    21562191              if (ExitFlag == 0) ExitFlag = 1;
    2157               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
     2192              if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
    21582193                ExitFlag = 255;
    21592194                DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl);
    21602195                performCriticalExit();
    21612196              } else {
    2162                 ArgcList.insert(argptr-1);
    2163                 ArgcList.insert(argptr);
    2164                 ArgcList.insert(argptr+1);
    2165                 ArgcList.insert(argptr+2);
     2197                SaveFlag = true;
     2198                j = -1;
     2199                DoLog(1) && (Log() << Verbose(1) << "Scaling all ion positions by factor." << endl);
     2200                factor = new double[NDIM];
     2201                factor[0] = atof(argv[argptr]);
     2202                factor[1] = atof(argv[argptr+1]);
     2203                factor[2] = atof(argv[argptr+2]);
     2204                mol->Scale((const double ** const)&factor);
     2205                double * const cell_size = World::getInstance().getDomain();
     2206                for (int i=0;i<NDIM;i++) {
     2207                  j += i+1;
     2208                  x[i] = atof(argv[NDIM+i]);
     2209                  cell_size[j]*=factor[i];
     2210                }
     2211                delete[](factor);
    21662212                argptr+=3;
    21672213              }
     
    21742220                performCriticalExit();
    21752221              } else {
    2176                 ArgcList.insert(argptr-1);
    2177                 ArgcList.insert(argptr);
    2178                 ArgcList.insert(argptr+1);
    2179                 ArgcList.insert(argptr+2);
    2180                 ArgcList.insert(argptr+3);
    2181                 ArgcList.insert(argptr+4);
    2182                 ArgcList.insert(argptr+5);
     2222                SaveFlag = true;
     2223                j = -1;
     2224                DoLog(1) && (Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl);
     2225                double * const cell_size = World::getInstance().getDomain();
     2226                for (int i=0;i<6;i++) {
     2227                  cell_size[i] = atof(argv[argptr+i]);
     2228                }
     2229                // center
     2230                mol->CenterInBox();
    21832231                argptr+=6;
    21842232              }
     
    21912239                performCriticalExit();
    21922240              } else {
    2193                 ArgcList.insert(argptr-1);
    2194                 ArgcList.insert(argptr);
    2195                 ArgcList.insert(argptr+1);
    2196                 ArgcList.insert(argptr+2);
    2197                 ArgcList.insert(argptr+3);
    2198                 ArgcList.insert(argptr+4);
    2199                 ArgcList.insert(argptr+5);
     2241                SaveFlag = true;
     2242                j = -1;
     2243                DoLog(1) && (Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl);
     2244                double * const cell_size = World::getInstance().getDomain();
     2245                for (int i=0;i<6;i++) {
     2246                  cell_size[i] = atof(argv[argptr+i]);
     2247                }
     2248                // center
     2249                mol->BoundInBox();
    22002250                argptr+=6;
    22012251              }
     
    22082258                performCriticalExit();
    22092259              } else {
    2210                 ArgcList.insert(argptr-1);
    2211                 ArgcList.insert(argptr);
    2212                 ArgcList.insert(argptr+1);
    2213                 ArgcList.insert(argptr+2);
     2260                SaveFlag = true;
     2261                j = -1;
     2262                DoLog(1) && (Log() << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl);
     2263                // make every coordinate positive
     2264                mol->CenterEdge(&x);
     2265                // update Box of atoms by boundary
     2266                mol->SetBoxDimension(&x);
     2267                // translate each coordinate by boundary
     2268                double * const cell_size = World::getInstance().getDomain();
     2269                j=-1;
     2270                for (int i=0;i<NDIM;i++) {
     2271                  j += i+1;
     2272                  x[i] = atof(argv[argptr+i]);
     2273                  cell_size[j] += x[i]*2.;
     2274                }
     2275                mol->Translate((const Vector *)&x);
    22142276                argptr+=3;
    22152277              }
     
    22172279            case 'O':
    22182280              if (ExitFlag == 0) ExitFlag = 1;
    2219               ArgcList.insert(argptr-1);
     2281              SaveFlag = true;
     2282              DoLog(1) && (Log() << Verbose(1) << "Centering atoms on edge and setting box dimensions." << endl);
     2283              x.Zero();
     2284              mol->CenterEdge(&x);
     2285              mol->SetBoxDimension(&x);
    22202286              argptr+=0;
    22212287              break;
     
    22272293                performCriticalExit();
    22282294              } else {
    2229                 mol->getAtomCount();
    22302295                SaveFlag = true;
    22312296                DoLog(1) && (Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl);
     
    22472312                mol->CreateAdjacencyList(atof(argv[argptr]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
    22482313                DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl);
    2249                 if (mol->hasBondStructure()) {
     2314                if (mol->first->next != mol->last) {
    22502315                  ExitFlag = mol->FragmentMolecule(atoi(argv[argptr+1]), &configuration);
    22512316                }
     
    23362401                performCriticalExit();
    23372402              } else {
    2338                 ArgcList.insert(argptr-1);
    2339                 ArgcList.insert(argptr);
    2340                 ArgcList.insert(argptr+1);
    2341                 ArgcList.insert(argptr+2);
    2342                 argptr+=3;
     2403                SaveFlag = true;
     2404                double * const cell_size = World::getInstance().getDomain();
     2405                for (int axis = 1; axis <= NDIM; axis++) {
     2406                  int faktor = atoi(argv[argptr++]);
     2407                  int count;
     2408                  const element ** Elements;
     2409                  Vector ** vectors;
     2410                  if (faktor < 1) {
     2411                    DoeLog(1) && (eLog()<< Verbose(1) << "Repetition factor mus be greater than 1!" << endl);
     2412                    faktor = 1;
     2413                  }
     2414                  mol->CountAtoms();  // recount atoms
     2415                  if (mol->AtomCount != 0) {  // if there is more than none
     2416                    count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
     2417                    Elements = new const element *[count];
     2418                    vectors = new Vector *[count];
     2419                    j = 0;
     2420                    first = mol->start;
     2421                    while (first->next != mol->end) {  // make a list of all atoms with coordinates and element
     2422                      first = first->next;
     2423                      Elements[j] = first->type;
     2424                      vectors[j] = &first->x;
     2425                      j++;
     2426                    }
     2427                    if (count != j)
     2428                      DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
     2429                    x.Zero();
     2430                    y.Zero();
     2431                    y[abs(axis)-1] = cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
     2432                    for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
     2433                      x += y; // per factor one cell width further
     2434                      for (int k=count;k--;) { // go through every atom of the original cell
     2435                        first = World::getInstance().createAtom(); // create a new body
     2436                        first->x = (*vectors[k]) + x;
     2437                        first->type = Elements[k];  // insert original element
     2438                        mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
     2439                      }
     2440                    }
     2441                    // free memory
     2442                    delete[](Elements);
     2443                    delete[](vectors);
     2444                    // correct cell size
     2445                    if (axis < 0) { // if sign was negative, we have to translate everything
     2446                      x =(-(faktor-1)) * y;
     2447                      mol->Translate(&x);
     2448                    }
     2449                    cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     2450                  }
     2451                }
    23432452              }
    23442453              break;
     
    23522461    } while (argptr < argc);
    23532462    if (SaveFlag)
    2354       configuration.SaveAll(*ConfigFileName, periode, molecules);
     2463      configuration.SaveAll(ConfigFileName, periode, molecules);
    23552464  } else {  // no arguments, hence scan the elements db
    23562465    if (periode->LoadPeriodentafel(configuration.databasepath))
     
    23672476void cleanUp(){
    23682477  World::purgeInstance();
     2478  Log() << Verbose(0) <<  "Maximum of allocated memory: "
     2479    << MemoryUsageObserver::getInstance()->getMaximumUsedMemory() << endl;
     2480  Log() << Verbose(0) <<  "Remaining non-freed memory: "
     2481    << MemoryUsageObserver::getInstance()->getUsedMemorySize() << endl;
     2482  MemoryUsageObserver::purgeInstance();
    23692483  logger::purgeInstance();
    23702484  errorLogger::purgeInstance();
     
    23742488  ActionRegistry::purgeInstance();
    23752489  ActionHistory::purgeInstance();
    2376   Memory::getState();
    23772490}
    23782491
     
    23802493{
    23812494    config *configuration = World::getInstance().getConfig();
    2382     // while we are non interactive, we want to abort from asserts
    2383     //ASSERT_DO(Assert::Abort);
    2384     molecule *mol = NULL;
    23852495    Vector x, y, z, n;
    23862496    ifstream test;
     
    23892499    char **Arguments = NULL;
    23902500    int ArgcSize = 0;
    2391     int ExitFlag = 0;
    23922501    bool ArgumentsCopied = false;
    2393     char *ConfigFileName = new char[MAXSTRINGSIZE];
    2394 
    2395     // print version check whether arguments are present at all
     2502
    23962503    cout << ESPACKVersion << endl;
    2397     if (argc < 2) {
    2398       cout << "Obtain help with " << argv[0] << " -h." << endl;
    2399       cleanUp();
    2400       Memory::getState();
    2401       return(1);
    2402     }
    2403 
    24042504
    24052505    setVerbosity(0);
    24062506    // need to init the history before any action is created
    24072507    ActionHistory::init();
    2408 
    2409     // In the interactive mode, we can leave the user the choice in case of error
    2410     ASSERT_DO(Assert::Ask);
    2411 
    2412     // from this moment on, we need to be sure to deeinitialize in the correct order
    2413     // this is handled by the cleanup function
    2414     atexit(cleanUp);
    24152508
    24162509    // Parse command line options and if present create respective UI
     
    24192512      ArgcList.insert(0); // push back program!
    24202513      ArgcList.insert(1); // push back config file name
     2514      char ConfigFileName[MAXSTRINGSIZE];
    24212515      // handle arguments by ParseCommandLineOptions()
    2422       ExitFlag = ParseCommandLineOptions(argc,argv,World::getInstance().getMolecules(),World::getInstance().getPeriode(),*World::getInstance().getConfig(), &ConfigFileName, ArgcList);
     2516      ParseCommandLineOptions(argc,argv,World::getInstance().getMolecules(),World::getInstance().getPeriode(),*World::getInstance().getConfig(), (char *&)ConfigFileName, ArgcList);
    24232517      // copy all remaining arguments to a new argv
    2424       Arguments = new char *[ArgcList.size()];
     2518      Arguments = Malloc<char *>(ArgcList.size(), "main - **Arguments");
    24252519      cout << "The following arguments are handled by CommandLineParser: ";
    24262520      for (set<int>::iterator ArgcRunner = ArgcList.begin(); ArgcRunner != ArgcList.end(); ++ArgcRunner) {
    2427         Arguments[ArgcSize] = new char[strlen(argv[*ArgcRunner])+2];
     2521        Arguments[ArgcSize] = Malloc<char>(strlen(argv[*ArgcRunner])+2, "main - *Arguments[]");
    24282522        strcpy(Arguments[ArgcSize], argv[*ArgcRunner]);
    24292523        cout << " " << argv[*ArgcRunner];
     
    24502544    }
    24512545
    2452     Log() << Verbose(0) << "Saving to " << ConfigFileName << "." << endl;
    2453     World::getInstance().getConfig()->SaveAll(ConfigFileName, World::getInstance().getPeriode(), World::getInstance().getMolecules());
     2546    if(World::getInstance().getPeriode()->StorePeriodentafel(configuration->databasepath))
     2547        Log() << Verbose(0) << "Saving of elements.db successful." << endl;
     2548
     2549    else
     2550        Log() << Verbose(0) << "Saving of elements.db failed." << endl;
    24542551
    24552552  // free the new argv
    24562553  if (ArgumentsCopied) {
    24572554    for (int i=0; i<ArgcSize;i++)
    2458       delete[](Arguments[i]);
    2459     delete[](Arguments);
     2555      Free(&Arguments[i]);
     2556    Free(&Arguments);
    24602557  }
    2461   delete[](ConfigFileName);
    2462 
    2463   return (ExitFlag == 1 ? 0 : ExitFlag);
     2558
     2559  cleanUp();
     2560  Memory::getState();
     2561  return (0);
    24642562}
    24652563
Note: See TracChangeset for help on using the changeset viewer.