Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/builder.cpp

    r481601 r046783  
    6868#include "periodentafel.hpp"
    6969#include "version.h"
     70#include "World.hpp"
    7071
    7172/********************************************* Subsubmenu routine ************************************/
    … …  
    103104        Log() << Verbose(0) << "Enter absolute coordinates." << endl;
    104105        first = new atom;
    105         first->x.AskPosition(mol->cell_size, false);
     106        first->x.AskPosition(World::get()->cell_size, false);
    106107        first->type = periode->AskElement();  // give type
    107108        mol->AddAtom(first);  // add to molecule
    … …  
    114115          if (!valid) eLog() << Verbose(2) << "Resulting position out of cell." << endl;
    115116          Log() << Verbose(0) << "Enter reference coordinates." << endl;
    116           x.AskPosition(mol->cell_size, true);
     117          x.AskPosition(World::get()->cell_size, true);
    117118          Log() << Verbose(0) << "Enter relative coordinates." << endl;
    118           first->x.AskPosition(mol->cell_size, false);
     119          first->x.AskPosition(World::get()->cell_size, false);
    119120          first->x.AddVector((const Vector *)&x);
    120121          Log() << Verbose(0) << "\n";
    … …  
    131132          second = mol->AskAtom("Enter atom number: ");
    132133          Log() << Verbose(0) << "Enter relative coordinates." << endl;
    133           first->x.AskPosition(mol->cell_size, false);
     134          first->x.AskPosition(World::get()->cell_size, false);
    134135          for (int i=NDIM;i--;) {
    135136            first->x.x[i] += second->x.x[i];
    … …  
    358359    case 'b': // normal vector of mirror plane
    359360      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
    360       n.AskPosition(mol->cell_size,0);
     361      n.AskPosition(World::get()->cell_size,0);
    361362      n.Normalize();
    362363      break;
    … …  
    425426    case 'b': // normal vector of mirror plane
    426427      Log() << Verbose(0) << "Enter normal vector of mirror plane." << endl;
    427       n.AskPosition(mol->cell_size,0);
     428      n.AskPosition(World::get()->cell_size,0);
    428429      n.Normalize();
    429430      break;
    … …  
    655656static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
    656657{
    657   atom *first, *second;
     658  atom *first, *second, *third;
    658659  molecule *mol = NULL;
    659660  Vector x,y,z,n; // coordinates for absolute point in cell volume
    … …  
    667668  Log() << Verbose(0) << "r - remove an atom" << endl;
    668669  Log() << Verbose(0) << "b - scale a bond between atoms" << endl;
     670  Log() << Verbose(0) << "t - turn an atom round another bond" << endl;
    669671  Log() << Verbose(0) << "u - change an atoms element" << endl;
    670672  Log() << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
    … …  
    748750      break;
    749751
     752    case 't': // turn/rotate atom
     753      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
     754        if ((*ListRunner)->ActiveFlag) {
     755          mol = *ListRunner;
     756          Log() << Verbose(0) << "Turning atom around another bond - first is atom to turn, second (central) and third specify bond" << endl;
     757          first = mol->AskAtom("Enter turning atom: ");
     758          second = mol->AskAtom("Enter central atom: ");
     759          third  = mol->AskAtom("Enter bond atom: ");
     760          cout << Verbose(0) << "Enter new angle in degrees: ";
     761          double tmp = 0.;
     762          cin >> tmp;
     763          // calculate old angle
     764          x.CopyVector((const Vector *)&first->x);
     765          x.SubtractVector((const Vector *)&second->x);
     766          y.CopyVector((const Vector *)&third->x);
     767          y.SubtractVector((const Vector *)&second->x);
     768          double alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
     769          cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
     770          cout << Verbose(0) << alpha << " degrees" << endl;
     771          // rotate
     772          z.MakeNormalVector(&x,&y);
     773          x.RotateVector(&z,(alpha-tmp)*M_PI/180.);
     774          x.AddVector(&second->x);
     775          first->x.CopyVector(&x);
     776          // check new angle
     777          x.CopyVector((const Vector *)&first->x);
     778          x.SubtractVector((const Vector *)&second->x);
     779          alpha = (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.);
     780          cout << Verbose(0) << "new Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
     781          cout << Verbose(0) << alpha << " degrees" << endl;
     782        }
     783      break;
     784
    750785    case 'u': // change an atom's element
    751786      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
    … …  
    831866          x.Zero();
    832867          y.Zero();
    833           y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
     868          y.x[abs(axis)-1] = World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
    834869          for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
    835870            x.AddVector(&y); // per factor one cell width further
    … …  
    854889            mol->Translate(&x);
    855890          }
    856           mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     891          World::get()->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    857892        }
    858893      }
    … …  
    911946        Log() << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
    912947        Log() << Verbose(0) << "Enter translation vector." << endl;
    913         x.AskPosition(mol->cell_size,0);
     948        x.AskPosition(World::get()->cell_size,0);
    914949        mol->Center.AddVector((const Vector *)&x);
    915950     }
    … …  
    9711006        // center at set box dimensions
    9721007        mol->CenterEdge(&center);
    973         mol->cell_size[0] = center.x[0];
    974         mol->cell_size[1] = 0;
    975         mol->cell_size[2] = center.x[1];
    976         mol->cell_size[3] = 0;
    977         mol->cell_size[4] = 0;
    978         mol->cell_size[5] = center.x[2];
     1008        double * const cell_size = World::get()->cell_size;
     1009        cell_size[0] = center.x[0];
     1010        cell_size[1] = 0;
     1011        cell_size[2] = center.x[1];
     1012        cell_size[3] = 0;
     1013        cell_size[4] = 0;
     1014        cell_size[5] = center.x[2];
    9791015        molecules->insert(mol);
    9801016      }
    … …  
    13871423  enum ConfigStatus configPresent = absent;
    13881424  clock_t start,end;
     1425  double MaxDistance = -1;
    13891426  int argptr;
    13901427  molecule *mol = NULL;
    … …  
    14121449            Log() << Verbose(0) << "\t-B xx xy xz yy yz zz\tBound atoms by domain with given symmetric matrix of (xx,xy,xz,yy,yz,zz)." << endl;
    14131450            Log() << Verbose(0) << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
    1414             Log() << Verbose(0) << "\t-C <Z> <output> <bin output>\tPair Correlation analysis." << endl;
     1451            Log() << Verbose(0) << "\t-C <type> [params] <output> <bin output> <BinWidth> <BinStart> <BinEnd>\tPair Correlation analysis." << endl;
    14151452            Log() << Verbose(0) << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
    14161453            Log() << Verbose(0) << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
    … …  
    14181455            Log() << Verbose(0) << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
    14191456            Log() << Verbose(0) << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
    1420             Log() << Verbose(0) << "\t-F <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
     1457            Log() << Verbose(0) << "\t-F <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
     1458            Log() << Verbose(0) << "\t-FF <MaxDistance> <xyz of filler> <dist_x> <dist_y> <dist_z> <epsilon> <randatom> <randmol> <DoRotate>\tFilling Box with water molecules." << endl;
    14211459            Log() << Verbose(0) << "\t-g <file>\tParses a bond length table from the given file." << endl;
    14221460            Log() << Verbose(0) << "\t-h/-H/-?\tGive this help screen." << endl;
    14231461            Log() << Verbose(0) << "\t-I\t Dissect current system of molecules into a set of disconnected (subgraphs of) molecules." << endl;
     1462            Log() << Verbose(0) << "\t-j\t<path> Store all bonds to file." << endl;
     1463            Log() << Verbose(0) << "\t-J\t<path> Store adjacency per atom to file." << endl;
    14241464            Log() << Verbose(0) << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
    14251465            Log() << Verbose(0) << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
    … …  
    14401480            Log() << Verbose(0) << "\t-v\t\tsets verbosity (more is more)." << endl;
    14411481            Log() << Verbose(0) << "\t-V\t\tGives version information." << endl;
     1482            Log() << Verbose(0) << "\t-X\t\tset default name of a molecule." << endl;
    14421483            Log() << Verbose(0) << "Note: config files must not begin with '-' !" << endl;
    14431484            return (1);
    … …  
    14781519            Log() << Verbose(0) << "I won't parse trajectories." << endl;
    14791520            configuration.FastParsing = true;
     1521            break;
     1522          case 'X':
     1523            {
     1524              char **name = &(World::get()->DefaultName);
     1525              delete[](*name);
     1526              const int length = strlen(argv[argptr]);
     1527              *name = new char[length+2];
     1528              strncpy(*name, argv[argptr], length);
     1529              Log() << Verbose(0) << "Default name of new molecules set to " << *name << "." << endl;
     1530            }
    14801531            break;
    14811532          default:   // no match? Step on
    … …  
    16691720                  }
    16701721              }
    1671               if (mol == NULL) {
     1722              if ((mol == NULL) && (!molecules->ListOfMolecules.empty())) {
    16721723                mol = *(molecules->ListOfMolecules.begin());
    1673                 mol->ActiveFlag = true;
     1724                if (mol != NULL)
     1725                  mol->ActiveFlag = true;
    16741726              }
    16751727              break;
    16761728            case 'C':
    16771729              if (ExitFlag == 0) ExitFlag = 1;
    1678               if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr][0] == '-') || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-')) {
     1730              if ((argptr >= argc)) {
    16791731                ExitFlag = 255;
    1680                 eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <Z> <output> <bin output>" << endl;
     1732                eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output> <BinStart> <BinEnd>" << endl;
    16811733                performCriticalExit();
    16821734              } else {
    1683                 ofstream output(argv[argptr+1]);
    1684                 ofstream binoutput(argv[argptr+2]);
    1685                 const double radius = 5.;
    1686 
    1687                 // get the boundary
    1688                 class molecule *Boundary = NULL;
    1689                 class Tesselation *TesselStruct = NULL;
    1690                 const LinkedCell *LCList = NULL;
    1691                 // find biggest molecule
    1692                 int counter  = 0;
    1693                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    1694                   if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
    1695                     Boundary = *BigFinder;
    1696                   }
    1697                   counter++;
     1735                switch(argv[argptr][0]) {
     1736                  case 'E':
     1737                    {
     1738                      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] == '-')) {
     1739                        ExitFlag = 255;
     1740                        eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C E <Z1> <Z2> <output> <bin output>" << endl;
     1741                        performCriticalExit();
     1742                      } else {
     1743                        ofstream output(argv[argptr+3]);
     1744                        ofstream binoutput(argv[argptr+4]);
     1745                        const double BinStart = atof(argv[argptr+5]);
     1746                        const double BinEnd = atof(argv[argptr+6]);
     1747
     1748                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
     1749                        element *elemental2 = periode->FindElement((const int) atoi(argv[argptr+2]));
     1750                        PairCorrelationMap *correlationmap = PairCorrelation(molecules, elemental, elemental2);
     1751                        //OutputCorrelationToSurface(&output, correlationmap);
     1752                        BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
     1753                        OutputCorrelation ( &binoutput, binmap );
     1754                        output.close();
     1755                        binoutput.close();
     1756                        delete(binmap);
     1757                        delete(correlationmap);
     1758                        argptr+=7;
     1759                      }
     1760                    }
     1761                    break;
     1762
     1763                  case 'P':
     1764                    {
     1765                      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] == '-')) {
     1766                        ExitFlag = 255;
     1767                        eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C P <Z1> <x> <y> <z> <output> <bin output>" << endl;
     1768                        performCriticalExit();
     1769                      } else {
     1770                        ofstream output(argv[argptr+5]);
     1771                        ofstream binoutput(argv[argptr+6]);
     1772                        const double BinStart = atof(argv[argptr+7]);
     1773                        const double BinEnd = atof(argv[argptr+8]);
     1774
     1775                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
     1776                        Vector *Point = new Vector((const double) atof(argv[argptr+1]),(const double) atof(argv[argptr+2]),(const double) atof(argv[argptr+3]));
     1777                        CorrelationToPointMap *correlationmap = CorrelationToPoint(molecules, elemental, Point);
     1778                        //OutputCorrelationToSurface(&output, correlationmap);
     1779                        BinPairMap *binmap = BinData( correlationmap, 0.5, BinStart, BinEnd );
     1780                        OutputCorrelation ( &binoutput, binmap );
     1781                        output.close();
     1782                        binoutput.close();
     1783                        delete(Point);
     1784                        delete(binmap);
     1785                        delete(correlationmap);
     1786                        argptr+=9;
     1787                      }
     1788                    }
     1789                    break;
     1790
     1791                  case 'S':
     1792                    {
     1793                      if ((argptr+6 >= argc) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+4])) || (!IsValidNumber(argv[argptr+5])) || (!IsValidNumber(argv[argptr+6])) || (argv[argptr+1][0] == '-') || (argv[argptr+2][0] == '-') || (argv[argptr+3][0] == '-')) {
     1794                        ExitFlag = 255;
     1795                        eLog() << Verbose(0) << "Not enough or invalid arguments given for pair correlation analysis: -C S <Z> <output> <bin output> <BinWidth> <BinStart> <BinEnd>" << endl;
     1796                        performCriticalExit();
     1797                      } else {
     1798                        ofstream output(argv[argptr+2]);
     1799                        ofstream binoutput(argv[argptr+3]);
     1800                        const double radius = 4.;
     1801                        const double BinWidth = atof(argv[argptr+4]);
     1802                        const double BinStart = atof(argv[argptr+5]);
     1803                        const double BinEnd = atof(argv[argptr+6]);
     1804                        double LCWidth = 20.;
     1805                        if (BinEnd > 0) {
     1806                          if (BinEnd > 2.*radius)
     1807                              LCWidth = BinEnd;
     1808                          else
     1809                            LCWidth = 2.*radius;
     1810                        }
     1811
     1812                        // get the boundary
     1813                        class molecule *Boundary = NULL;
     1814                        class Tesselation *TesselStruct = NULL;
     1815                        const LinkedCell *LCList = NULL;
     1816                        // find biggest molecule
     1817                        int counter  = 0;
     1818                        for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
     1819                          if ((Boundary == NULL) || (Boundary->AtomCount < (*BigFinder)->AtomCount)) {
     1820                            Boundary = *BigFinder;
     1821                          }
     1822                          counter++;
     1823                        }
     1824                        bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
     1825                        counter = 0;
     1826                        for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
     1827                          Actives[counter++] = (*BigFinder)->ActiveFlag;
     1828                          (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
     1829                        }
     1830                        LCList = new LinkedCell(Boundary, LCWidth);
     1831                        element *elemental = periode->FindElement((const int) atoi(argv[argptr+1]));
     1832                        FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
     1833                        //int ranges[NDIM] = {1,1,1};
     1834                        CorrelationToSurfaceMap *surfacemap = CorrelationToSurface( molecules, elemental, TesselStruct, LCList); // for Periodic..(): ..., ranges );
     1835                        OutputCorrelationToSurface(&output, surfacemap);
     1836                        // check whether radius was appropriate
     1837                        {
     1838                        double start; double end;
     1839                        GetMinMax( surfacemap, start, end);
     1840                        if (LCWidth < end)
     1841                          eLog() << Verbose(1) << "Linked Cell width is smaller than the found range of values! Bins can only be correct up to: " << radius << "." << endl;
     1842                        }
     1843                        BinPairMap *binmap = BinData( surfacemap, BinWidth, BinStart, BinEnd );
     1844                        OutputCorrelation ( &binoutput, binmap );
     1845                        output.close();
     1846                        binoutput.close();
     1847                        for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
     1848                          (*BigFinder)->ActiveFlag = Actives[counter++];
     1849                        Free(&Actives);
     1850                        delete(LCList);
     1851                        delete(TesselStruct);
     1852                        delete(binmap);
     1853                        delete(surfacemap);
     1854                        argptr+=7;
     1855                      }
     1856                    }
     1857                    break;
     1858
     1859                  default:
     1860                    ExitFlag = 255;
     1861                    eLog() << Verbose(0) << "Invalid type given for pair correlation analysis: -C <type: E/P/S> [more params] <output> <bin output>" << endl;
     1862                    performCriticalExit();
     1863                    break;
    16981864                }
    1699                 bool *Actives = Malloc<bool>(counter, "ParseCommandLineOptions() - case C -- *Actives");
    1700                 counter = 0;
    1701                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++) {
    1702                   Actives[counter++] = (*BigFinder)->ActiveFlag;
    1703                   (*BigFinder)->ActiveFlag = (*BigFinder == Boundary) ? false : true;
    1704                 }
    1705                 LCList = new LinkedCell(Boundary, 2.*radius);
    1706                 element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
    1707                 FindNonConvexBorder(Boundary, TesselStruct, LCList, radius, NULL);
    1708                 int ranges[NDIM] = {1,1,1};
    1709                 CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( molecules, elemental, TesselStruct, LCList, ranges );
    1710                 OutputCorrelationToSurface(&output, surfacemap);
    1711                 BinPairMap *binmap = BinData( surfacemap, 0.5, 0., 0. );
    1712                 OutputCorrelation ( &binoutput, binmap );
    1713                 output.close();
    1714                 binoutput.close();
    1715                 for (MoleculeList::iterator BigFinder = molecules->ListOfMolecules.begin(); BigFinder != molecules->ListOfMolecules.end(); BigFinder++)
    1716                   (*BigFinder)->ActiveFlag = Actives[counter++];
    1717                 Free(&Actives);
    1718                 delete(LCList);
    1719                 delete(TesselStruct);
    1720                 argptr+=3;
    17211865              }
    17221866              break;
    … …  
    17371881            case 'F':
    17381882              if (ExitFlag == 0) ExitFlag = 1;
    1739               if (argptr+6 >=argc) {
     1883              MaxDistance = -1;
     1884              if (argv[argptr-1][2] == 'F') {
     1885                // fetch first argument as max distance to surface
     1886                MaxDistance = atof(argv[argptr++]);
     1887                Log() << Verbose(0) << "Filling with maximum layer distance of " << MaxDistance << "." << endl;
     1888              }
     1889              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]))) {
    17401890                ExitFlag = 255;
    1741                 eLog() << Verbose(0) << "Not enough or invalid arguments given for filling box with water: -F <dist_x> <dist_y> <dist_z> <boundary> <randatom> <randmol> <DoRotate>" << endl;
     1891                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;
    17421892                performCriticalExit();
    17431893              } else {
    … …  
    17461896                // construct water molecule
    17471897                molecule *filler = new molecule(periode);
     1898                if (!filler->AddXYZFile(argv[argptr])) {
     1899                  eLog() << Verbose(0) << "Could not parse filler molecule from " << argv[argptr] << "." << endl;
     1900                }
     1901                filler->SetNameFromFilename(argv[argptr]);
     1902                configuration.BG->ConstructBondGraph(filler);
    17481903                molecule *Filling = NULL;
    1749                 atom *second = NULL, *third = NULL;
    1750 //                first = new atom();
    1751 //                first->type = periode->FindElement(5);
    1752 //                first->x.Zero();
    1753 //                filler->AddAtom(first);
    1754                 first = new atom();
    1755                 first->type = periode->FindElement(1);
    1756                 first->x.Init(0.441, -0.143, 0.);
    1757                 filler->AddAtom(first);
    1758                 second = new atom();
    1759                 second->type = periode->FindElement(1);
    1760                 second->x.Init(-0.464, 1.137, 0.0);
    1761                 filler->AddAtom(second);
    1762                 third = new atom();
    1763                 third->type = periode->FindElement(8);
    1764                 third->x.Init(-0.464, 0.177, 0.);
    1765                 filler->AddAtom(third);
    1766                 filler->AddBond(first, third, 1);
    1767                 filler->AddBond(second, third, 1);
    17681904                // call routine
    17691905                double distance[NDIM];
    17701906                for (int i=0;i<NDIM;i++)
    1771                   distance[i] = atof(argv[argptr+i]);
    1772                 Filling = FillBoxWithMolecule(molecules, filler, configuration, distance, atof(argv[argptr+3]), atof(argv[argptr+4]), atof(argv[argptr+5]), atoi(argv[argptr+6]));
     1907                  distance[i] = atof(argv[argptr+i+1]);
     1908                Filling = FillBoxWithMolecule(molecules, filler, configuration, MaxDistance, distance, atof(argv[argptr+4]), atof(argv[argptr+5]), atof(argv[argptr+6]), atoi(argv[argptr+7]));
    17731909                if (Filling != NULL) {
    17741910                  Filling->ActiveFlag = false;
    … …  
    17931929              }
    17941930              break;
     1931
     1932            case 'J':
     1933              if (ExitFlag == 0) ExitFlag = 1;
     1934              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1935                ExitFlag =255;
     1936                eLog() << Verbose(0) << "Missing path of adjacency file: -j <path>" << endl;
     1937                performCriticalExit();
     1938              } else {
     1939                Log() << Verbose(0) << "Storing adjacency to path " << argv[argptr] << "." << endl;
     1940                configuration.BG->ConstructBondGraph(mol);
     1941                mol->StoreAdjacencyToFile(argv[argptr]);
     1942                argptr+=1;
     1943              }
     1944              break;
     1945
     1946            case 'j':
     1947              if (ExitFlag == 0) ExitFlag = 1;
     1948              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
     1949                ExitFlag =255;
     1950                eLog() << Verbose(0) << "Missing path of bonds file: -j <path>" << endl;
     1951                performCriticalExit();
     1952              } else {
     1953                Log() << Verbose(0) << "Storing bonds to path " << argv[argptr] << "." << endl;
     1954                configuration.BG->ConstructBondGraph(mol);
     1955                mol->StoreBondsToFile(argv[argptr]);
     1956                argptr+=1;
     1957              }
     1958              break;
     1959
    17951960            case 'N':
    17961961              if (ExitFlag == 0) ExitFlag = 1;
    … …  
    19542119                factor[2] = atof(argv[argptr+2]);
    19552120                mol->Scale((const double ** const)&factor);
     2121                double * const cell_size = World::get()->cell_size;
    19562122                for (int i=0;i<NDIM;i++) {
    19572123                  j += i+1;
    19582124                  x.x[i] = atof(argv[NDIM+i]);
    1959                   mol->cell_size[j]*=factor[i];
     2125                  cell_size[j]*=factor[i];
    19602126                }
    19612127                delete[](factor);
    … …  
    19732139                j = -1;
    19742140                Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
     2141                double * const cell_size = World::get()->cell_size;
    19752142                for (int i=0;i<6;i++) {
    1976                   mol->cell_size[i] = atof(argv[argptr+i]);
     2143                  cell_size[i] = atof(argv[argptr+i]);
    19772144                }
    19782145                // center
    … …  
    19912158                j = -1;
    19922159                Log() << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
     2160                double * const cell_size = World::get()->cell_size;
    19932161                for (int i=0;i<6;i++) {
    1994                   mol->cell_size[i] = atof(argv[argptr+i]);
     2162                  cell_size[i] = atof(argv[argptr+i]);
    19952163                }
    19962164                // center
    … …  
    20142182                mol->SetBoxDimension(&x);
    20152183                // translate each coordinate by boundary
     2184                double * const cell_size = World::get()->cell_size;
    20162185                j=-1;
    20172186                for (int i=0;i<NDIM;i++) {
    20182187                  j += i+1;
    20192188                  x.x[i] = atof(argv[argptr+i]);
    2020                   mol->cell_size[j] += x.x[i]*2.;
     2189                  cell_size[j] += x.x[i]*2.;
    20212190                }
    20222191                mol->Translate((const Vector *)&x);
    … …  
    21492318              } else {
    21502319                SaveFlag = true;
     2320                double * const cell_size = World::get()->cell_size;
    21512321                for (int axis = 1; axis <= NDIM; axis++) {
    21522322                  int faktor = atoi(argv[argptr++]);
    … …  
    21752345                    x.Zero();
    21762346                    y.Zero();
    2177                     y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
     2347                    y.x[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
    21782348                    for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
    21792349                      x.AddVector(&y); // per factor one cell width further
    … …  
    21962366                      mol->Translate(&x);
    21972367                    }
    2198                     mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
     2368                    cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
    21992369                  }
    22002370                }
    … …  
    22652435  if (molecules->ListOfMolecules.size() == 0) {
    22662436    mol = new molecule(periode);
    2267     if (mol->cell_size[0] == 0.) {
     2437    double * const cell_size = World::get()->cell_size;
     2438    if (cell_size[0] == 0.) {
    22682439      Log() << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
    22692440      for (int i=0;i<6;i++) {
    22702441        Log() << Verbose(1) << "Cell size" << i << ": ";
    2271         cin >> mol->cell_size[i];
     2442        cin >> cell_size[i];
    22722443      }
    22732444    }
Note: See TracChangeset for help on using the changeset viewer.