Changes in src/builder.cpp [9907e7:f65e1f]
- File:
-
- 1 edited
-
src/builder.cpp (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/builder.cpp
r9907e7 rf65e1f 53 53 54 54 #include <cstring> 55 #include <cstdlib>56 55 57 56 #include "analysis_bonds.hpp" … … 865 864 866 865 mol->CountAtoms(); // recount atoms 867 if (mol-> getAtomCount()!= 0) { // if there is more than none868 count = mol-> getAtomCount(); // is changed becausing of adding, thus has to be stored away beforehand866 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 869 868 Elements = new element *[count]; 870 869 vectors = new Vector *[count]; … … 1296 1295 // generate some KeySets 1297 1296 DoLog(0) && (Log() << Verbose(0) << "Generating KeySets." << endl); 1298 KeySet TestSets[mol-> getAtomCount()+1];1297 KeySet TestSets[mol->AtomCount+1]; 1299 1298 i=1; 1300 1299 while (Walker->next != mol->end) { … … 1307 1306 DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in KeySets." << endl); 1308 1307 KeySetTestPair test; 1309 test = TestSets[mol-> getAtomCount()-1].insert(Walker->nr);1308 test = TestSets[mol->AtomCount-1].insert(Walker->nr); 1310 1309 if (test.second) { 1311 1310 DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl); … … 1313 1312 DoLog(1) && (Log() << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl); 1314 1313 } 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); 1317 1316 1318 1317 // constructing Graph structure … … 1322 1321 // insert KeySets into Subgraphs 1323 1322 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++) { 1325 1324 Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.))); 1326 1325 } 1327 1326 DoLog(0) && (Log() << Verbose(0) << "Testing insertion of already present item in Subgraph." << endl); 1328 1327 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.))); 1330 1329 if (test2.second) { 1331 1330 DoLog(1) && (Log() << Verbose(1) << "Insertion worked?!" << endl); … … 1495 1494 */ 1496 1495 static 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) 1498 1497 { 1499 1498 Vector x,y,z,n; // coordinates for absolute point in cell volume … … 1513 1512 molecule *mol = NULL; 1514 1513 string BondGraphFileName("\n"); 1515 bool DatabasePathGiven = false;1514 strncpy(configuration.databasepath, LocalPath, MAXSTRINGSIZE-1); 1516 1515 1517 1516 if (argc > 1) { // config file specified as option … … 1530 1529 break; 1531 1530 case 'v': 1532 setVerbosity(atoi(argv[argptr])); 1531 ArgcList.insert(argptr-1); 1532 return(1); 1533 break; 1534 case 'V': 1533 1535 ArgcList.insert(argptr-1); 1534 1536 ArgcList.insert(argptr); 1535 1537 argptr++; 1536 1538 break; 1537 case 'V':1538 ArgcList.insert(argptr-1);1539 return(1);1540 break;1541 1539 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;1550 1540 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 } 1551 1555 break; 1552 1556 case 'e': … … 1557 1561 DoLog(0) && (Log() << Verbose(0) << "Using " << argv[argptr] << " as elements database." << endl); 1558 1562 strncpy (configuration.databasepath, argv[argptr], MAXSTRINGSIZE-1); 1559 DatabasePathGiven = true;1560 1563 argptr+=1; 1561 1564 } … … 1571 1574 } 1572 1575 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;1584 1576 case 'n': 1585 1577 DoLog(0) && (Log() << Verbose(0) << "I won't parse trajectories." << endl); … … 1589 1581 { 1590 1582 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); 1592 1584 } 1593 1585 break; … … 1601 1593 1602 1594 // 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 } 1611 1602 // 3b. Find config file name and parse if possible, also BondGraphFileName 1612 1603 if (argv[1][0] != '-') { … … 1622 1613 } else { 1623 1614 DoLog(0) && (Log() << Verbose(0) << "Empty configuration file." << endl); 1624 strcpy(*ConfigFileName, argv[1]);1615 ConfigFileName = argv[1]; 1625 1616 configPresent = empty; 1626 1617 output.close(); … … 1628 1619 } else { 1629 1620 test.close(); 1630 strcpy(*ConfigFileName, argv[1]);1621 ConfigFileName = argv[1]; 1631 1622 DoLog(1) && (Log() << Verbose(1) << "Specified config file found, parsing ... "); 1632 switch (configuration.TestSyntax( *ConfigFileName, periode)) {1623 switch (configuration.TestSyntax(ConfigFileName, periode)) { 1633 1624 case 1: 1634 1625 DoLog(0) && (Log() << Verbose(0) << "new syntax." << endl); 1635 configuration.Load( *ConfigFileName, BondGraphFileName, periode, molecules);1626 configuration.Load(ConfigFileName, BondGraphFileName, periode, molecules); 1636 1627 configPresent = present; 1637 1628 break; 1638 1629 case 0: 1639 1630 DoLog(0) && (Log() << Verbose(0) << "old syntax." << endl); 1640 configuration.LoadOld( *ConfigFileName, BondGraphFileName, periode, molecules);1631 configuration.LoadOld(ConfigFileName, BondGraphFileName, periode, molecules); 1641 1632 configPresent = present; 1642 1633 break; … … 1659 1650 mol = World::getInstance().createMolecule(); 1660 1651 mol->ActiveFlag = true; 1661 if ( *ConfigFileName != NULL)1662 mol->SetNameFromFilename( *ConfigFileName);1652 if (ConfigFileName != NULL) 1653 mol->SetNameFromFilename(ConfigFileName); 1663 1654 molecules->insert(mol); 1664 1655 } … … 1714 1705 if (first->type != NULL) { 1715 1706 mol->AddAtom(first); // add to molecule 1716 if ((configPresent == empty) && (mol-> getAtomCount()!= 0))1707 if ((configPresent == empty) && (mol->AtomCount != 0)) 1717 1708 configPresent = present; 1718 1709 } else … … 1727 1718 if (configPresent == present) { 1728 1719 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; 1729 1731 case 'D': 1730 1732 if (ExitFlag == 0) ExitFlag = 1; … … 1732 1734 DoLog(1) && (Log() << Verbose(1) << "Depth-First-Search Analysis." << endl); 1733 1735 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis 1734 int *MinimumRingSize = new int[mol-> getAtomCount()];1736 int *MinimumRingSize = new int[mol->AtomCount]; 1735 1737 atom ***ListOfLocalAtoms = NULL; 1736 1738 class StackClass<bond *> *BackEdgeStack = NULL; … … 1752 1754 delete(Subgraphs); 1753 1755 for (int i=0;i<FragmentCounter;i++) 1754 delete[](ListOfLocalAtoms[i]);1755 delete[](ListOfLocalAtoms);1756 Free(&ListOfLocalAtoms[i]); 1757 Free(&ListOfLocalAtoms); 1756 1758 } 1757 1759 delete(BackEdgeStack); … … 1793 1795 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] == '-')) { 1794 1796 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); 1796 1798 performCriticalExit(); 1797 1799 } else { … … 1808 1810 else 1809 1811 correlationmap = PairCorrelation(molecules, elemental, elemental2); 1810 OutputPairCorrelation(&output, correlationmap);1812 //OutputCorrelationToSurface(&output, correlationmap); 1811 1813 BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd ); 1812 1814 OutputCorrelation ( &binoutput, binmap ); … … 1824 1826 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] == '-')) { 1825 1827 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); 1827 1829 performCriticalExit(); 1828 1830 } else { … … 1839 1841 else 1840 1842 correlationmap = CorrelationToPoint(molecules, elemental, Point); 1841 OutputCorrelationToPoint(&output, correlationmap);1843 //OutputCorrelationToSurface(&output, correlationmap); 1842 1844 BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd ); 1843 1845 OutputCorrelation ( &binoutput, binmap ); … … 1880 1882 int counter = 0; 1881 1883 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)) { 1883 1885 Boundary = *BigFinder; 1884 1886 } 1885 1887 counter++; 1886 1888 } 1887 bool *Actives = new bool[counter];1889 bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives"); 1888 1890 counter = 0; 1889 1891 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) { … … 1911 1913 output.close(); 1912 1914 binoutput.close(); 1913 counter = 0;1914 1915 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) 1915 1916 (*BigFinder)->ActiveFlag = Actives[counter++]; 1916 delete[](Actives);1917 Free(&Actives); 1917 1918 delete(LCList); 1918 1919 delete(TesselStruct); … … 1940 1941 performCriticalExit(); 1941 1942 } else { 1942 mol->getAtomCount();1943 1943 SaveFlag = true; 1944 1944 DoLog(1) && (Log() << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl); … … 1950 1950 case 'F': 1951 1951 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 } 1967 2000 break; 1968 2001 case 'A': … … 1974 2007 } else { 1975 2008 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(); 1979 2012 argptr+=1; 1980 2013 } … … 2026 2059 int counter = 0; 2027 2060 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)) { 2029 2063 Boundary = *BigFinder; 2030 2064 } 2031 2065 counter++; 2032 2066 } 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); 2034 2068 start = clock(); 2035 2069 LCList = new LinkedCell(Boundary, atof(argv[argptr])*2.); … … 2066 2100 if ((argptr >= argc) || (argv[argptr][0] == '-')) { 2067 2101 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); 2069 2103 performCriticalExit(); 2070 2104 } else { … … 2098 2132 case 'R': 2099 2133 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]))) { 2101 2135 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); 2103 2137 performCriticalExit(); 2104 2138 } else { 2105 2139 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); 2118 2151 } 2152 } else { 2153 DoeLog(1) && (eLog()<< Verbose(1) << "Removal failed due to missing atoms on molecule or wrong id." << endl); 2119 2154 } 2120 argptr+= 4;2155 argptr+=2; 2121 2156 } 2122 2157 break; … … 2155 2190 case 's': 2156 2191 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])) ) { 2158 2193 ExitFlag = 255; 2159 2194 DoeLog(0) && (eLog()<< Verbose(0) << "Not enough or invalid arguments given for scaling: -s <factor_x> [factor_y] [factor_z]" << endl); 2160 2195 performCriticalExit(); 2161 2196 } 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); 2166 2212 argptr+=3; 2167 2213 } … … 2174 2220 performCriticalExit(); 2175 2221 } 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(); 2183 2231 argptr+=6; 2184 2232 } … … 2191 2239 performCriticalExit(); 2192 2240 } 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(); 2200 2250 argptr+=6; 2201 2251 } … … 2208 2258 performCriticalExit(); 2209 2259 } 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); 2214 2276 argptr+=3; 2215 2277 } … … 2217 2279 case 'O': 2218 2280 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); 2220 2286 argptr+=0; 2221 2287 break; … … 2227 2293 performCriticalExit(); 2228 2294 } else { 2229 mol->getAtomCount();2230 2295 SaveFlag = true; 2231 2296 DoLog(1) && (Log() << Verbose(1) << "Removing atom " << argv[argptr] << "." << endl); … … 2247 2312 mol->CreateAdjacencyList(atof(argv[argptr]), configuration.GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL); 2248 2313 DoLog(0) && (Log() << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl); 2249 if (mol-> hasBondStructure()) {2314 if (mol->first->next != mol->last) { 2250 2315 ExitFlag = mol->FragmentMolecule(atoi(argv[argptr+1]), &configuration); 2251 2316 } … … 2336 2401 performCriticalExit(); 2337 2402 } 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 } 2343 2452 } 2344 2453 break; … … 2352 2461 } while (argptr < argc); 2353 2462 if (SaveFlag) 2354 configuration.SaveAll( *ConfigFileName, periode, molecules);2463 configuration.SaveAll(ConfigFileName, periode, molecules); 2355 2464 } else { // no arguments, hence scan the elements db 2356 2465 if (periode->LoadPeriodentafel(configuration.databasepath)) … … 2367 2476 void cleanUp(){ 2368 2477 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(); 2369 2483 logger::purgeInstance(); 2370 2484 errorLogger::purgeInstance(); … … 2374 2488 ActionRegistry::purgeInstance(); 2375 2489 ActionHistory::purgeInstance(); 2376 Memory::getState();2377 2490 } 2378 2491 … … 2380 2493 { 2381 2494 config *configuration = World::getInstance().getConfig(); 2382 // while we are non interactive, we want to abort from asserts2383 //ASSERT_DO(Assert::Abort);2384 molecule *mol = NULL;2385 2495 Vector x, y, z, n; 2386 2496 ifstream test; … … 2389 2499 char **Arguments = NULL; 2390 2500 int ArgcSize = 0; 2391 int ExitFlag = 0;2392 2501 bool ArgumentsCopied = false; 2393 char *ConfigFileName = new char[MAXSTRINGSIZE]; 2394 2395 // print version check whether arguments are present at all 2502 2396 2503 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 2404 2504 2405 2505 setVerbosity(0); 2406 2506 // need to init the history before any action is created 2407 2507 ActionHistory::init(); 2408 2409 // In the interactive mode, we can leave the user the choice in case of error2410 ASSERT_DO(Assert::Ask);2411 2412 // from this moment on, we need to be sure to deeinitialize in the correct order2413 // this is handled by the cleanup function2414 atexit(cleanUp);2415 2508 2416 2509 // Parse command line options and if present create respective UI … … 2419 2512 ArgcList.insert(0); // push back program! 2420 2513 ArgcList.insert(1); // push back config file name 2514 char ConfigFileName[MAXSTRINGSIZE]; 2421 2515 // 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); 2423 2517 // copy all remaining arguments to a new argv 2424 Arguments = new char *[ArgcList.size()];2518 Arguments = Malloc<char *>(ArgcList.size(), "main - **Arguments"); 2425 2519 cout << "The following arguments are handled by CommandLineParser: "; 2426 2520 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[]"); 2428 2522 strcpy(Arguments[ArgcSize], argv[*ArgcRunner]); 2429 2523 cout << " " << argv[*ArgcRunner]; … … 2450 2544 } 2451 2545 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; 2454 2551 2455 2552 // free the new argv 2456 2553 if (ArgumentsCopied) { 2457 2554 for (int i=0; i<ArgcSize;i++) 2458 delete[](Arguments[i]);2459 delete[](Arguments);2555 Free(&Arguments[i]); 2556 Free(&Arguments); 2460 2557 } 2461 delete[](ConfigFileName); 2462 2463 return (ExitFlag == 1 ? 0 : ExitFlag); 2558 2559 cleanUp(); 2560 Memory::getState(); 2561 return (0); 2464 2562 } 2465 2563
Note:
See TracChangeset
for help on using the changeset viewer.
