Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/Tesselation/boundary.cpp

    r052c10 r2b7375  
    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);
    332331  return BoundaryPoints;
    333332};
    … …  
    655654};
    656655
    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  */
    667 void 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  */
    775 void 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  */
    939 void 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  */
    969 bool 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  */
    1001 bool 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  */
    1038 inline 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  */
    1063 void 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  */
    1237 void 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 
    1431656/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
    1432657 * \param *out output stream for debugging
Note: See TracChangeset for help on using the changeset viewer.