Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/Tesselation/boundary.cpp

    r2b7375 r052c10  
    183183  LineMap LinesOnBoundary;
    184184  TriangleMap TrianglesOnBoundary;
    185   Vector MolCenter = mol->DetermineCenterOfAll();
     185  Vector *MolCenter = mol->DetermineCenterOfAll();
    186186  Vector helper;
    187187  BoundariesTestPair BoundaryTestPair;
    … …  
    207207    // Boundaries stores non-const TesselPoint ref, hence we need iterator here
    208208    for (molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    209       ProjectedVector = (*iter)->getPosition() - (MolCenter);
     209      ProjectedVector = (*iter)->getPosition() - (*MolCenter);
    210210      ProjectedVector.ProjectOntoPlane(AxisVector);
    211211
    … …  
    233233          LOG(2, "Keeping new vector due to larger projected distance " << ProjectedVectorNorm << ".");
    234234        } else if (fabs(ProjectedVectorNorm - BoundaryTestPair.first->second.first) < MYEPSILON) {
    235           helper = (*iter)->getPosition() - (MolCenter);
     235          helper = (*iter)->getPosition() - (*MolCenter);
    236236          const double oldhelperNorm = helper.NormSquared();
    237           helper = BoundaryTestPair.first->second.second->getPosition() - (MolCenter);
     237          helper = BoundaryTestPair.first->second.second->getPosition() - (*MolCenter);
    238238          if (helper.NormSquared() < oldhelperNorm) {
    239239            BoundaryTestPair.first->second.second = (*iter);
    … …  
    289289        {
    290290          Vector SideA, SideB, SideC, SideH;
    291           SideA = left->second.second->getPosition() - (MolCenter);
     291          SideA = left->second.second->getPosition() - (*MolCenter);
    292292          SideA.ProjectOntoPlane(AxisVector);
    293293          //          LOG(1, "SideA: " << SideA);
    294294
    295           SideB = right->second.second->getPosition() -(MolCenter);
     295          SideB = right->second.second->getPosition() -(*MolCenter);
    296296          SideB.ProjectOntoPlane(AxisVector);
    297297          //          LOG(1, "SideB: " << SideB);
    … …  
    301301          //          LOG(1, "SideC: " << SideC);
    302302
    303           SideH = runner->second.second->getPosition() -(MolCenter);
     303          SideH = runner->second.second->getPosition() -(*MolCenter);
    304304          SideH.ProjectOntoPlane(AxisVector);
    305305          //          LOG(1, "SideH: " << SideH);
    … …  
    329329    } while (flag);
    330330  }
     331  delete(MolCenter);
    331332  return BoundaryPoints;
    332333};
    … …  
    654655};
    655656
     657/** Creates multiples of the by \a *mol given cluster and suspends them in water with a given final density.
     658 * We get cluster volume by Tesselation::getVolumeOfConvexEnvelope() and its diameters by GetDiametersOfCluster()
     659 * TODO: Here, we need a VolumeOfGeneralEnvelope (i.e. non-convex one)
     660 * \param *out output stream for debugging
     661 * \param *configuration needed for path to store convex envelope file
     662 * \param *mol molecule structure representing the cluster
     663 * \param *&TesselStruct Tesselation structure with triangles on return
     664 * \param ClusterVolume guesstimated cluster volume, if equal 0 we used Tesselation::getVolumeOfConvexEnvelope() instead.
     665 * \param celldensity desired average density in final cell
     666 */
     667void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity)
     668{
     669        //Info FunctionInfo(__func__);
     670  bool IsAngstroem = true;
     671  double *GreatestDiameter = NULL;
     672  Boundaries *BoundaryPoints = NULL;
     673  class Tesselation *TesselStruct = NULL;
     674  Vector BoxLengths;
     675  int repetition[NDIM] = { 1, 1, 1 };
     676  int TotalNoClusters = 1;
     677  double totalmass = 0.;
     678  double clustervolume = 0.;
     679  double cellvolume = 0.;
     680
     681  // transform to PAS by Action
     682  Vector MainAxis(0.,0.,1.);
     683  mol->RotateToPrincipalAxisSystem(MainAxis);
     684
     685  IsAngstroem = configuration->GetIsAngstroem();
     686  BoundaryPoints = GetBoundaryPoints(mol, TesselStruct);
     687  GreatestDiameter = GetDiametersOfCluster(BoundaryPoints, mol, TesselStruct, IsAngstroem);
     688  PointCloudAdaptor< molecule > cloud(mol, mol->name);
     689  LinkedCell_deprecated *LCList = new LinkedCell_deprecated(cloud, 10.);
     690  FindConvexBorder(mol, BoundaryPoints, TesselStruct, (const LinkedCell_deprecated *&)LCList, NULL);
     691  delete (LCList);
     692  delete[] BoundaryPoints;
     693
     694
     695  // some preparations beforehand
     696  if (ClusterVolume == 0)
     697    clustervolume = TesselStruct->getVolumeOfConvexEnvelope(configuration->GetIsAngstroem());
     698  else
     699    clustervolume = ClusterVolume;
     700
     701  delete TesselStruct;
     702
     703  for (int i = 0; i < NDIM; i++)
     704    TotalNoClusters *= repetition[i];
     705
     706  // sum up the atomic masses
     707  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
     708      totalmass += (*iter)->getType()->getMass();
     709  }
     710  LOG(0, "RESULT: The summed mass is " << setprecision(10) << totalmass << " atomicmassunit.");
     711  LOG(0, "RESULT: The average density is " << setprecision(10) << totalmass / clustervolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
     712
     713  // solve cubic polynomial
     714  LOG(1, "Solving equidistant suspension in water problem ...");
     715  if (IsAngstroem)
     716    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_A - (totalmass / clustervolume)) / (celldensity - 1);
     717  else
     718    cellvolume = (TotalNoClusters * totalmass / SOLVENTDENSITY_a0 - (totalmass / clustervolume)) / (celldensity - 1);
     719  LOG(1, "Cellvolume needed for a density of " << celldensity << " g/cm^3 is " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
     720
     721  double minimumvolume = TotalNoClusters * (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
     722  LOG(1, "Minimum volume of the convex envelope contained in a rectangular box is " << minimumvolume << " atomicmassunit/" << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
     723  if (minimumvolume > cellvolume) {
     724    ELOG(1, "the containing box already has a greater volume than the envisaged cell volume!");
     725    LOG(0, "Setting Box dimensions to minimum possible, the greatest diameters.");
     726    for (int i = 0; i < NDIM; i++)
     727      BoxLengths[i] = GreatestDiameter[i];
     728    mol->CenterEdge(&BoxLengths);
     729  } else {
     730    BoxLengths[0] = (repetition[0] * GreatestDiameter[0] + repetition[1] * GreatestDiameter[1] + repetition[2] * GreatestDiameter[2]);
     731    BoxLengths[1] = (repetition[0] * repetition[1] * GreatestDiameter[0] * GreatestDiameter[1] + repetition[0] * repetition[2] * GreatestDiameter[0] * GreatestDiameter[2] + repetition[1] * repetition[2] * GreatestDiameter[1] * GreatestDiameter[2]);
     732    BoxLengths[2] = minimumvolume - cellvolume;
     733    double x0 = 0.;
     734    double x1 = 0.;
     735    double x2 = 0.;
     736    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x0, &x1, &x2) == 1) // either 1 or 3 on return
     737      LOG(0, "RESULT: The resulting spacing is: " << x0 << " .");
     738    else {
     739      LOG(0, "RESULT: The resulting spacings are: " << x0 << " and " << x1 << " and " << x2 << " .");
     740      x0 = x2; // sorted in ascending order
     741    }
     742
     743    cellvolume = 1.;
     744    for (int i = 0; i < NDIM; i++) {
     745      BoxLengths[i] = repetition[i] * (x0 + GreatestDiameter[i]);
     746      cellvolume *= BoxLengths[i];
     747    }
     748
     749    // set new box dimensions
     750    LOG(0, "Translating to box with these boundaries.");
     751    mol->SetBoxDimension(&BoxLengths);
     752    mol->CenterInBox();
     753  }
     754  delete[] GreatestDiameter;
     755  // update Box of atoms by boundary
     756  mol->SetBoxDimension(&BoxLengths);
     757  LOG(0, "RESULT: The resulting cell dimensions are: " << BoxLengths[0] << " and " << BoxLengths[1] << " and " << BoxLengths[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");
     758};
     759
     760
     761/** Fills the empty space around other molecules' surface of the simulation box with a filler.
     762 * \param *out output stream for debugging
     763 * \param *List list of molecules already present in box
     764 * \param *TesselStruct contains tesselated surface
     765 * \param *filler molecule which the box is to be filled with
     766 * \param configuration contains box dimensions
     767 * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
     768 * \param distance[NDIM] distance between filling molecules in each direction
     769 * \param boundary length of boundary zone between molecule and filling mollecules
     770 * \param epsilon distance to surface which is not filled
     771 * \param RandAtomDisplacement maximum distance for random displacement per atom
     772 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
     773 * \param DoRandomRotation true - do random rotiations, false - don't
     774 */
     775void FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, const double MaxDistance, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
     776{
     777        //Info FunctionInfo(__func__);
     778  molecule *Filling = World::getInstance().createMolecule();
     779  Vector CurrentPosition;
     780  int N[NDIM];
     781  int n[NDIM];
     782  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
     783  RealSpaceMatrix Rotations;
     784  const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
     785  Vector AtomTranslations;
     786  Vector FillerTranslations;
     787  Vector FillerDistance;
     788  Vector Inserter;
     789  double FillIt = false;
     790  bond::ptr Binder;
     791  double phi[NDIM];
     792  map<molecule *, Tesselation *> TesselStruct;
     793  map<molecule *, LinkedCell_deprecated *> LCList;
     794
     795  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
     796    if ((*ListRunner)->getAtomCount() > 0) {
     797      LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
     798      PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
     799      LCList[(*ListRunner)] = new LinkedCell_deprecated(cloud, 10.); // get linked cell list
     800      LOG(1, "Pre-creating tesselation for molecule " << *ListRunner << ".");
     801      TesselStruct[(*ListRunner)] = NULL;
     802      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell_deprecated *&)LCList[(*ListRunner)], 5., NULL);
     803    }
     804
     805  // Center filler at origin
     806  filler->CenterEdge(&Inserter);
     807  const int FillerCount = filler->getAtomCount();
     808  LOG(2, "INFO: Filler molecule has the following bonds:");
     809  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
     810    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
     811    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
     812        BondRunner != ListOfBonds.end();
     813        ++BondRunner) {
     814      if ((*BondRunner)->leftatom == *AtomRunner)
     815        LOG(2, "  " << *(*BondRunner));
     816    }
     817  }
     818
     819  atom * CopyAtoms[FillerCount];
     820
     821  // calculate filler grid in [0,1]^3
     822  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
     823  for(int i=0;i<NDIM;i++)
     824    N[i] = (int) ceil(1./FillerDistance[i]);
     825  LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");
     826
     827  // initialize seed of random number generator to current time
     828  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
     829  const double rng_min = random.min();
     830  const double rng_max = random.max();
     831  //srand ( time(NULL) );
     832
     833  // go over [0,1]^3 filler grid
     834  for (n[0] = 0; n[0] < N[0]; n[0]++)
     835    for (n[1] = 0; n[1] < N[1]; n[1]++)
     836      for (n[2] = 0; n[2] < N[2]; n[2]++) {
     837        // calculate position of current grid vector in untransformed box
     838        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
     839        // create molecule random translation vector ...
     840        for (int i=0;i<NDIM;i++)
     841          FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
     842        LOG(2, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");
     843
     844        // go through all atoms
     845        for (int i=0;i<FillerCount;i++)
     846          CopyAtoms[i] = NULL;
     847
     848        // have same rotation angles for all molecule's atoms
     849        if (DoRandomRotation)
     850          for (int i=0;i<NDIM;i++)
     851            phi[i] = (random()/(rng_max-rng_min))*(2.*M_PI);
     852
     853        // atom::clone is not const member function, hence we need iterator here
     854        for(molecule::iterator iter = filler->begin(); iter !=filler->end();++iter){
     855
     856          // create atomic random translation vector ...
     857          for (int i=0;i<NDIM;i++)
     858            AtomTranslations[i] = RandomAtomDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
     859
     860          // ... and rotation matrix
     861          if (DoRandomRotation) {
     862            Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
     863            Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
     864            Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
     865            Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
     866            Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
     867            Rotations.set(1,2,              sin(phi[1])                                                    );
     868            Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
     869            Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
     870            Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
     871          }
     872
     873          // ... and put at new position
     874          Inserter = (*iter)->getPosition();
     875          if (DoRandomRotation)
     876            Inserter *= Rotations;
     877          Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
     878
     879          // check whether inserter is inside box
     880          Inserter *= MInverse;
     881          FillIt = true;
     882          for (int i=0;i<NDIM;i++)
     883            FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
     884          Inserter *= M;
     885
     886          // Check whether point is in- or outside
     887          for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
     888            // get linked cell list
     889            if (TesselStruct[(*ListRunner)] != NULL) {
     890              const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
     891              FillIt = FillIt && (distance > boundary) && ((MaxDistance < 0) || (MaxDistance > distance));
     892            }
     893          }
     894          // insert into Filling
     895          if (FillIt) {
     896            LOG(1, "INFO: Position at " << Inserter << " is outer point.");
     897            // copy atom ...
     898            CopyAtoms[(*iter)->getNr()] = (*iter)->clone();
     899            (*CopyAtoms[(*iter)->getNr()]).setPosition(Inserter);
     900            Filling->AddAtom(CopyAtoms[(*iter)->getNr()]);
     901            LOG(1, "Filling atom " << **iter << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->getNr()]->getPosition()) << ".");
     902          } else {
     903            LOG(1, "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance.");
     904            CopyAtoms[(*iter)->getNr()] = NULL;
     905            continue;
     906          }
     907        }
     908        // go through all bonds and add as well
     909        for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
     910          const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
     911          for(BondList::const_iterator BondRunner = ListOfBonds.begin();
     912              BondRunner != ListOfBonds.end();
     913              ++BondRunner)
     914            if ((*BondRunner)->leftatom == *AtomRunner) {
     915              Binder = (*BondRunner);
     916              if ((CopyAtoms[Binder->leftatom->getNr()] != NULL) && (CopyAtoms[Binder->rightatom->getNr()] != NULL)) {
     917                LOG(3, "Adding Bond between " << *CopyAtoms[Binder->leftatom->getNr()] << " and " << *CopyAtoms[Binder->rightatom->getNr()]<< ".");
     918                Filling->AddBond(CopyAtoms[Binder->leftatom->getNr()], CopyAtoms[Binder->rightatom->getNr()], Binder->getDegree());
     919              }
     920            }
     921        }
     922      }
     923  for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
     924    delete LCList[*ListRunner];
     925    delete TesselStruct[(*ListRunner)];
     926  }
     927};
     928
     929/** Rotates given molecule \a Filling and moves its atoms according to given
     930 *  \a RandomAtomDisplacement.
     931 *
     932 *  Note that for rotation to be sensible, the molecule should be centered at
     933 *  the origin. This is not done here!
     934 *
     935 *  \param &Filling molecule whose atoms to displace
     936 *  \param RandomAtomDisplacement magnitude of random displacement
     937 *  \param &Rotations 3D rotation matrix (or unity if no rotation desired)
     938 */
     939void RandomizeMoleculePositions(
     940    molecule *&Filling,
     941    double RandomAtomDisplacement,
     942    RealSpaceMatrix &Rotations,
     943    RandomNumberGenerator &random
     944    )
     945{
     946  const double rng_min = random.min();
     947  const double rng_max = random.max();
     948
     949  Vector AtomTranslations;
     950  for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) {
     951    Vector temp = (*miter)->getPosition();
     952    temp *= Rotations;
     953    (*miter)->setPosition(temp);
     954    // create atomic random translation vector ...
     955    for (int i=0;i<NDIM;i++)
     956      AtomTranslations[i] = RandomAtomDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
     957    (*miter)->setPosition((*miter)->getPosition() + AtomTranslations);
     958  }
     959}
     960
     961/** Removes all atoms of a molecule outside.
     962 *
     963 * If the molecule is empty, it is removed as well.
     964 *
     965 * @param Filling molecule whose atoms to check, removed if eventually left
     966 *        empty.
     967 * @return true - atoms had to be removed, false - no atoms have been removed
     968 */
     969bool RemoveAtomsOutsideDomain(molecule *&Filling)
     970{
     971  bool status = false;
     972  Box &Domain = World::getInstance().getDomain();
     973  // check if all is still inside domain
     974  for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {
     975    // check whether each atom is inside box
     976    if (!Domain.isInside((*miter)->getPosition())) {
     977      status = true;
     978      atom *Walker = *miter;
     979      ++miter;
     980      World::getInstance().destroyAtom(Walker);
     981    } else {
     982      ++miter;
     983    }
     984  }
     985  if (Filling->empty()) {
     986    LOG(0, "Removing molecule " << Filling->getName() << ", all atoms have been removed.");
     987    World::getInstance().destroyMolecule(Filling);
     988    Filling = NULL;
     989  }
     990  return status;
     991}
     992
     993/** Checks whether there are no atoms inside a sphere around \a CurrentPosition
     994 *  except those atoms present in \a *filler.
     995 *  If filler is NULL, then we just call LinkedCell_deprecated::GetPointsInsideSphere() and
     996 *  check whether the return list is empty.
     997 * @param *filler
     998 * @param boundary
     999 * @param CurrentPosition
     1000 */
     1001bool isSpaceAroundPointVoid(
     1002    LinkedCell_deprecated *LC,
     1003    molecule *filler,
     1004    const double boundary,
     1005    Vector &CurrentPosition)
     1006{
     1007  size_t compareTo = 0;
     1008  TesselPointSTLList* liste = LC->GetPointsInsideSphere(boundary == 0. ? MYEPSILON : boundary, &CurrentPosition);
     1009  if (filler != NULL) {
     1010    for (TesselPointSTLList::const_iterator iter = liste->begin();
     1011        iter != liste->end();
     1012        ++iter) {
     1013      for (molecule::iterator miter = filler->begin();
     1014          miter != filler->end();
     1015          ++miter) {
     1016        if (*iter == *miter)
     1017          ++compareTo;
     1018      }
     1019    }
     1020  }
     1021  const bool result = (liste->size() == compareTo);
     1022  if (!result) {
     1023    LOG(0, "Skipping because of the following atoms:");
     1024    for (TesselPointSTLList::const_iterator iter = liste->begin();
     1025        iter != liste->end();
     1026        ++iter) {
     1027      LOG(0, **iter);
     1028    }
     1029  }
     1030  delete liste;
     1031  return result;
     1032}
     1033
     1034/** Sets given 3x3 matrix to a random rotation matrix.
     1035 *
     1036 * @param a matrix to set
     1037 */
     1038inline void setRandomRotation(RealSpaceMatrix &a)
     1039{
     1040  double phi[NDIM];
     1041  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
     1042  const double rng_min = random.min();
     1043  const double rng_max = random.max();
     1044
     1045  for (int i=0;i<NDIM;i++) {
     1046    phi[i] = (random()/(rng_max-rng_min))*(2.*M_PI);
     1047    LOG(4, "DEBUG: Random angle is " << phi[i] << ".");
     1048  }
     1049
     1050  a.setRotation(phi);
     1051}
     1052
     1053/** Fills the empty space of the simulation box with water.
     1054 * \param *filler molecule which the box is to be filled with
     1055 * \param configuration contains box dimensions
     1056 * \param distance[NDIM] distance between filling molecules in each direction
     1057 * \param boundary length of boundary zone between molecule and filling molecules
     1058 * \param RandAtomDisplacement maximum distance for random displacement per atom
     1059 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
     1060 * \param MinDistance minimum distance to boundary of domain and present molecules
     1061 * \param DoRandomRotation true - do random rotations, false - don't
     1062 */
     1063void FillVoidWithMolecule(
     1064    molecule *&filler,
     1065    config &configuration,
     1066    const double distance[NDIM],
     1067    const double boundary,
     1068    const double RandomAtomDisplacement,
     1069    const double RandomMolDisplacement,
     1070    const double MinDistance,
     1071    const bool DoRandomRotation
     1072    )
     1073{
     1074  //Info FunctionInfo(__func__);
     1075  molecule *Filling = NULL;
     1076  Vector CurrentPosition;
     1077  int N[NDIM];
     1078  int n[NDIM];
     1079  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
     1080  RealSpaceMatrix Rotations;
     1081  const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
     1082  Vector FillerTranslations;
     1083  Vector FillerDistance;
     1084  Vector Inserter;
     1085  double FillIt = false;
     1086  Vector firstInserter;
     1087  bool firstInsertion = true;
     1088  const Box &Domain = World::getInstance().getDomain();
     1089  map<molecule *, LinkedCell_deprecated *> LCList;
     1090  std::vector<molecule *> List = World::getInstance().getAllMolecules();
     1091  MoleculeListClass *MolList = World::getInstance().getMolecules();
     1092
     1093  for (std::vector<molecule *>::iterator ListRunner = List.begin(); ListRunner != List.end(); ListRunner++)
     1094    if ((*ListRunner)->getAtomCount() > 0) {
     1095      LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
     1096      PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
     1097      LCList[(*ListRunner)] = new LinkedCell_deprecated(cloud, 10.); // get linked cell list
     1098    }
     1099
     1100  // Center filler at its center of gravity
     1101  Vector *gravity = filler->DetermineCenterOfGravity();
     1102  filler->CenterAtVector(gravity);
     1103  delete gravity;
     1104  //const int FillerCount = filler->getAtomCount();
     1105  LOG(2, "INFO: Filler molecule has the following bonds:");
     1106  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
     1107    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
     1108    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
     1109        BondRunner != ListOfBonds.end();
     1110        ++BondRunner)
     1111      if ((*BondRunner)->leftatom == *AtomRunner)
     1112        LOG(2, "  " << *(*BondRunner));
     1113  }
     1114
     1115  // calculate filler grid in [0,1]^3
     1116  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
     1117  for(int i=0;i<NDIM;i++)
     1118    N[i] = (int) ceil(1./FillerDistance[i]);
     1119  LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");
     1120
     1121  // initialize seed of random number generator to current time
     1122  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
     1123  const double rng_min = random.min();
     1124  const double rng_max = random.max();
     1125  //srand ( time(NULL) );
     1126
     1127  // go over [0,1]^3 filler grid
     1128  for (n[0] = 0; n[0] < N[0]; n[0]++)
     1129    for (n[1] = 0; n[1] < N[1]; n[1]++)
     1130      for (n[2] = 0; n[2] < N[2]; n[2]++) {
     1131        // calculate position of current grid vector in untransformed box
     1132        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
     1133        // create molecule random translation vector ...
     1134        for (int i=0;i<NDIM;i++) // have the random values [-1,1]*RandomMolDisplacement
     1135          FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
     1136        LOG(2, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");
     1137
     1138        // ... and rotation matrix
     1139        if (DoRandomRotation)
     1140          setRandomRotation(Rotations);
     1141        else
     1142          Rotations.setIdentity();
     1143
     1144
     1145        // Check whether there is anything too close by and whether atom is outside of domain
     1146        FillIt = true;
     1147        for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
     1148          FillIt = FillIt && isSpaceAroundPointVoid(
     1149              ListRunner->second,
     1150              (firstInsertion ? filler : NULL),
     1151              boundary,
     1152              CurrentPosition);
     1153          FillIt = FillIt && (Domain.isValid(CurrentPosition))
     1154              && ((Domain.DistanceToBoundary(CurrentPosition) - MinDistance) > -MYEPSILON);
     1155          if (!FillIt)
     1156            break;
     1157        }
     1158
     1159        // insert into Filling
     1160        if (FillIt) {
     1161          Inserter = CurrentPosition + FillerTranslations;
     1162          LOG(1, "INFO: Position at " << Inserter << " is void point.");
     1163          // fill!
     1164          Filling = filler->CopyMolecule();
     1165          RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations, random);
     1166          // translation
     1167          Filling->Translate(&Inserter);
     1168          // remove out-of-bounds atoms
     1169          const bool status = RemoveAtomsOutsideDomain(Filling);
     1170          if ((firstInsertion) && (!status)) { // no atom has been removed
     1171            // remove copied atoms and molecule again
     1172            Filling->removeAtomsinMolecule();
     1173            World::getInstance().destroyMolecule(Filling);
     1174            // and mark is final filler position
     1175            Filling = filler;
     1176            firstInsertion = false;
     1177            firstInserter = Inserter;
     1178          } else {
     1179            // TODO: Remove when World has no MoleculeListClass anymore
     1180            if (Filling)
     1181              MolList->insert(Filling);
     1182          }
     1183        } else {
     1184         LOG(1, "INFO: Position at " << Inserter << " is non-void point, within boundary or outside of MaxDistance.");
     1185          continue;
     1186        }
     1187      }
     1188
     1189  // have we inserted any molecules?
     1190  if (firstInsertion) {
     1191    // If not remove filler
     1192    for(molecule::iterator miter = filler->begin(); !filler->empty(); miter = filler->begin()) {
     1193      atom *Walker = *miter;
     1194      World::getInstance().destroyAtom(Walker);
     1195    }
     1196    World::getInstance().destroyMolecule(filler);
     1197  } else {
     1198    // otherwise translate and randomize to final position
     1199    if (DoRandomRotation)
     1200      setRandomRotation(Rotations);
     1201    else
     1202      Rotations.setIdentity();
     1203    RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations, random);
     1204    // translation
     1205    filler->Translate(&firstInserter);
     1206    // remove out-of-bounds atoms
     1207    RemoveAtomsOutsideDomain(filler);
     1208  }
     1209
     1210  LOG(0, MolList->ListOfMolecules.size() << " molecules have been inserted.");
     1211
     1212  for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
     1213    delete ListRunner->second;
     1214    LCList.erase(ListRunner);
     1215  }
     1216};
     1217
     1218
     1219/** Fills the empty space around other molecules' surface of the simulation box with a filler.
     1220 *
     1221 * Note that we use \a FindNonConvexBorder to determine the surface of the found molecules.
     1222 * There, we use a radius of twice the given \a boundary.
     1223 *
     1224 * \param *out output stream for debugging
     1225 * \param *List list of molecules already present in box
     1226 * \param *TesselStruct contains tesselated surface
     1227 * \param *filler molecule which the box is to be filled with
     1228 * \param configuration contains box dimensions
     1229 * \param MaxSurfaceDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
     1230 * \param distance[NDIM] distance between filling molecules in each direction
     1231 * \param boundary length of boundary zone between molecule and filling molecules
     1232 * \param MinDistance distance to boundary of domain which is not filled
     1233 * \param RandAtomDisplacement maximum distance for random displacement per atom
     1234 * \param RandMolDisplacement maximum distance for random displacement per filler molecule
     1235 * \param DoRandomRotation true - do random rotations, false - don't
     1236 */
     1237void FillBoxWithMolecule(
     1238    MoleculeListClass *MolList,
     1239    molecule *filler,
     1240    config &configuration,
     1241    const double MaxSurfaceDistance,
     1242    const double distance[NDIM],
     1243    const double boundary,
     1244    const double MinDistance,
     1245    const double RandomAtomDisplacement,
     1246    const double RandomMolDisplacement,
     1247    const bool DoRandomRotation)
     1248{
     1249  //Info FunctionInfo(__func__);
     1250  molecule *Filling = NULL;
     1251  Vector CurrentPosition;
     1252  int N[NDIM];
     1253  int n[NDIM];
     1254  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
     1255  RealSpaceMatrix Rotations;
     1256  const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
     1257  Vector FillerTranslations;
     1258  Vector FillerDistance;
     1259  Vector Inserter;
     1260  double FillIt = false;
     1261  Vector firstInserter;
     1262  bool firstInsertion = true;
     1263  const Box &Domain = World::getInstance().getDomain();
     1264  map<molecule *, LinkedCell_deprecated *> LCList;
     1265  std::vector<molecule *> List = World::getInstance().getAllMolecules();
     1266  map<molecule *, Tesselation *> TesselStruct;
     1267
     1268  for (MoleculeList::iterator ListRunner = MolList->ListOfMolecules.begin(); ListRunner != MolList->ListOfMolecules.end(); ListRunner++)
     1269    if ((*ListRunner)->getAtomCount() > 0) {
     1270      LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
     1271      PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
     1272      LCList[(*ListRunner)] = new LinkedCell_deprecated(cloud, 4.*boundary); // get linked cell list
     1273      LOG(1, "Pre-creating tesselation for molecule " << *ListRunner << ".");
     1274      TesselStruct[(*ListRunner)] = NULL;
     1275      FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell_deprecated *&)LCList[(*ListRunner)], 2.*boundary, NULL);
     1276    }
     1277
     1278  // Center filler at origin
     1279  filler->CenterEdge(&Inserter);
     1280//  const int FillerCount = filler->getAtomCount();
     1281  LOG(2, "INFO: Filler molecule has the following bonds:");
     1282  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
     1283    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
     1284    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
     1285        BondRunner != ListOfBonds.end();
     1286        ++BondRunner) {
     1287      if ((*BondRunner)->leftatom == *AtomRunner)
     1288        LOG(2, "  " << *(*BondRunner));
     1289    }
     1290  }
     1291
     1292//  atom * CopyAtoms[FillerCount];
     1293
     1294  setVerbosity(4);
     1295
     1296  // calculate filler grid in [0,1]^3
     1297  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
     1298  for(int i=0;i<NDIM;i++)
     1299    N[i] = (int) ceil(1./FillerDistance[i]);
     1300  LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");
     1301
     1302  // initialize seed of random number generator to current time
     1303  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
     1304  const double rng_min = random.min();
     1305  const double rng_max = random.max();
     1306  //srand ( time(NULL) );
     1307
     1308  // go over [0,1]^3 filler grid
     1309  for (n[0] = 0; n[0] < N[0]; n[0]++)
     1310    for (n[1] = 0; n[1] < N[1]; n[1]++)
     1311      for (n[2] = 0; n[2] < N[2]; n[2]++) {
     1312        // calculate position of current grid vector in untransformed box
     1313        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
     1314        // create molecule random translation vector ...
     1315        for (int i=0;i<NDIM;i++)
     1316          FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
     1317        LOG(1, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");
     1318
     1319        // ... and rotation matrix
     1320        if (DoRandomRotation)
     1321          setRandomRotation(Rotations);
     1322        else
     1323          Rotations.setIdentity();
     1324
     1325
     1326        // Check whether there is anything too close by and whether atom is outside of domain
     1327        FillIt = true;
     1328        for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
     1329          // check whether its void
     1330          FillIt = FillIt && isSpaceAroundPointVoid(
     1331              ListRunner->second,
     1332              (firstInsertion ? filler : NULL),
     1333              boundary,
     1334              CurrentPosition);
     1335          if (!FillIt) {
     1336            LOG(2, "REJECT: Position at " << Inserter << " is non-void.");
     1337            break;
     1338          }
     1339          // check whether inside domain
     1340          FillIt = FillIt && (Domain.isValid(CurrentPosition));
     1341          if (!FillIt) {
     1342            LOG(2, "REJECT: Position at " << CurrentPosition << " is "
     1343                << distance << ", hence outside of domain.");
     1344            break;
     1345          }
     1346          // check minimum distance to boundary
     1347          const double distance = (Domain.DistanceToBoundary(CurrentPosition) - MinDistance);
     1348          FillIt = FillIt && (distance > -MYEPSILON);
     1349          if (!FillIt) {
     1350            LOG(2, "REJECT: Position at " << CurrentPosition << " is " << distance << ", less than "
     1351                << MinDistance << " hence, too close to boundary.");
     1352            break;
     1353          }
     1354        }
     1355        // Check whether point is in- or outside of tesselations
     1356        if (FillIt) {
     1357          for (MoleculeList::iterator ListRunner = MolList->ListOfMolecules.begin(); ListRunner != MolList->ListOfMolecules.end(); ListRunner++) {
     1358            // get linked cell list
     1359            if (TesselStruct[(*ListRunner)] != NULL) {
     1360              const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
     1361              LOG(2, "INFO: Distance to surface is " << distance << ".");
     1362              FillIt = FillIt && ((distance == -1.) || (distance > boundary)) && ((MaxSurfaceDistance < 0) || (MaxSurfaceDistance > distance));
     1363              if (!FillIt) {
     1364                LOG(2, "REJECT: Position at " << CurrentPosition << " is in distance of " << distance
     1365                    << " to a surface, less than " << MaxSurfaceDistance  << " hence, too close.");
     1366                break;
     1367              }
     1368            }
     1369          }
     1370        }
     1371
     1372        // insert into Filling
     1373        if (FillIt) {
     1374          Inserter = CurrentPosition + FillerTranslations;
     1375          LOG(2, "ACCEPT: Position at " << CurrentPosition << " is void point.");
     1376          // fill!
     1377          Filling = filler->CopyMolecule();
     1378          RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations, random);
     1379          // translation
     1380          Filling->Translate(&Inserter);
     1381          // remove out-of-bounds atoms
     1382          const bool status = RemoveAtomsOutsideDomain(Filling);
     1383          if ((firstInsertion) && (!status)) { // no atom has been removed
     1384            // remove copied atoms and molecule again
     1385            Filling->removeAtomsinMolecule();
     1386            World::getInstance().destroyMolecule(Filling);
     1387            // and mark is final filler position
     1388            Filling = filler;
     1389            firstInsertion = false;
     1390            firstInserter = Inserter;
     1391          } else {
     1392            // TODO: Remove when World has no MoleculeListClass anymore
     1393            if (Filling)
     1394              MolList->insert(Filling);
     1395          }
     1396        } else {
     1397          LOG(2, "REJECT: Position at " << CurrentPosition << " is non-void point, within boundary or outside of MaxSurfaceDistance.");
     1398          continue;
     1399        }
     1400      }
     1401
     1402  // have we inserted any molecules?
     1403  if (firstInsertion) {
     1404    // If not remove filler
     1405    for(molecule::iterator miter = filler->begin(); !filler->empty(); miter = filler->begin()) {
     1406      atom *Walker = *miter;
     1407      World::getInstance().destroyAtom(Walker);
     1408    }
     1409    World::getInstance().destroyMolecule(filler);
     1410  } else {
     1411    // otherwise translate and randomize to final position
     1412    if (DoRandomRotation)
     1413      setRandomRotation(Rotations);
     1414    else
     1415      Rotations.setIdentity();
     1416    RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations, random);
     1417    // translation
     1418    filler->Translate(&firstInserter);
     1419    // remove out-of-bounds atoms
     1420    RemoveAtomsOutsideDomain(filler);
     1421  }
     1422
     1423  LOG(0, MolList->ListOfMolecules.size() << " molecules have been inserted.");
     1424
     1425  for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
     1426    delete ListRunner->second;
     1427    LCList.erase(ListRunner);
     1428  }
     1429};
     1430
    6561431/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
    6571432 * \param *out output stream for debugging
Note: See TracChangeset for help on using the changeset viewer.