Changeset 2528d8 for src


Ignore:
Timestamp:
Sep 10, 2014, 7:10:26 PM (10 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
aa55d0
Parents:
16c7dd
git-author:
Frederik Heber <heber@…> (09/01/14 15:57:23)
git-committer:
Frederik Heber <heber@…> (09/10/14 19:10:26)
Message:

Removed FillBoxWithMolecule(), FillVoidWithMolecule, and helper functions.

  • these are no longer needed as Actions have been removed.
Location:
src/Tesselation
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • src/Tesselation/boundary.cpp

    r16c7dd r2528d8  
    767767};
    768768
    769 
    770 /** Fills the empty space around other molecules' surface of the simulation box with a filler.
    771  * \param *out output stream for debugging
    772  * \param *List list of molecules already present in box
    773  * \param *TesselStruct contains tesselated surface
    774  * \param *filler molecule which the box is to be filled with
    775  * \param configuration contains box dimensions
    776  * \param MaxDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
    777  * \param distance[NDIM] distance between filling molecules in each direction
    778  * \param boundary length of boundary zone between molecule and filling mollecules
    779  * \param epsilon distance to surface which is not filled
    780  * \param RandAtomDisplacement maximum distance for random displacement per atom
    781  * \param RandMolDisplacement maximum distance for random displacement per filler molecule
    782  * \param DoRandomRotation true - do random rotiations, false - don't
    783  */
    784 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)
    785 {
    786         //Info FunctionInfo(__func__);
    787   molecule *Filling = World::getInstance().createMolecule();
    788   Vector CurrentPosition;
    789   int N[NDIM];
    790   int n[NDIM];
    791   const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
    792   RealSpaceMatrix Rotations;
    793   const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
    794   Vector AtomTranslations;
    795   Vector FillerTranslations;
    796   Vector FillerDistance;
    797   Vector Inserter;
    798   double FillIt = false;
    799   bond::ptr Binder;
    800   double phi[NDIM];
    801   map<molecule *, Tesselation *> TesselStruct;
    802   map<molecule *, LinkedCell_deprecated *> LCList;
    803 
    804   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++)
    805     if ((*ListRunner)->getAtomCount() > 0) {
    806       LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
    807       PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
    808       LCList[(*ListRunner)] = new LinkedCell_deprecated(cloud, 10.); // get linked cell list
    809       LOG(1, "Pre-creating tesselation for molecule " << *ListRunner << ".");
    810       TesselStruct[(*ListRunner)] = NULL;
    811       FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell_deprecated *&)LCList[(*ListRunner)], 5., NULL);
    812     }
    813 
    814   // Center filler at origin
    815   filler->CenterEdge();
    816   const int FillerCount = filler->getAtomCount();
    817   LOG(2, "INFO: Filler molecule has the following bonds:");
    818   for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
    819     const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    820     for(BondList::const_iterator BondRunner = ListOfBonds.begin();
    821         BondRunner != ListOfBonds.end();
    822         ++BondRunner) {
    823       if ((*BondRunner)->leftatom == *AtomRunner)
    824         LOG(2, "  " << *(*BondRunner));
    825     }
    826   }
    827 
    828   atom * CopyAtoms[FillerCount];
    829 
    830   // calculate filler grid in [0,1]^3
    831   FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
    832   for(int i=0;i<NDIM;i++)
    833     N[i] = (int) ceil(1./FillerDistance[i]);
    834   LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");
    835 
    836   // initialize seed of random number generator to current time
    837   RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
    838   const double rng_min = random.min();
    839   const double rng_max = random.max();
    840   //srand ( time(NULL) );
    841 
    842   // go over [0,1]^3 filler grid
    843   for (n[0] = 0; n[0] < N[0]; n[0]++)
    844     for (n[1] = 0; n[1] < N[1]; n[1]++)
    845       for (n[2] = 0; n[2] < N[2]; n[2]++) {
    846         // calculate position of current grid vector in untransformed box
    847         CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    848         // create molecule random translation vector ...
    849         for (int i=0;i<NDIM;i++)
    850           FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
    851         LOG(2, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");
    852 
    853         // go through all atoms
    854         for (int i=0;i<FillerCount;i++)
    855           CopyAtoms[i] = NULL;
    856 
    857         // have same rotation angles for all molecule's atoms
    858         if (DoRandomRotation)
    859           for (int i=0;i<NDIM;i++)
    860             phi[i] = (random()/(rng_max-rng_min))*(2.*M_PI);
    861 
    862         // atom::clone is not const member function, hence we need iterator here
    863         for(molecule::iterator iter = filler->begin(); iter !=filler->end();++iter){
    864 
    865           // create atomic random translation vector ...
    866           for (int i=0;i<NDIM;i++)
    867             AtomTranslations[i] = RandomAtomDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
    868 
    869           // ... and rotation matrix
    870           if (DoRandomRotation) {
    871             Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
    872             Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
    873             Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
    874             Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
    875             Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
    876             Rotations.set(1,2,              sin(phi[1])                                                    );
    877             Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
    878             Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
    879             Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
    880           }
    881 
    882           // ... and put at new position
    883           Inserter = (*iter)->getPosition();
    884           if (DoRandomRotation)
    885             Inserter *= Rotations;
    886           Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
    887 
    888           // check whether inserter is inside box
    889           Inserter *= MInverse;
    890           FillIt = true;
    891           for (int i=0;i<NDIM;i++)
    892             FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
    893           Inserter *= M;
    894 
    895           // Check whether point is in- or outside
    896           for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    897             // get linked cell list
    898             if (TesselStruct[(*ListRunner)] != NULL) {
    899               const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
    900               FillIt = FillIt && (distance > boundary) && ((MaxDistance < 0) || (MaxDistance > distance));
    901             }
    902           }
    903           // insert into Filling
    904           if (FillIt) {
    905             LOG(1, "INFO: Position at " << Inserter << " is outer point.");
    906             // copy atom ...
    907             CopyAtoms[(*iter)->getNr()] = (*iter)->clone();
    908             (*CopyAtoms[(*iter)->getNr()]).setPosition(Inserter);
    909             Filling->AddAtom(CopyAtoms[(*iter)->getNr()]);
    910             LOG(1, "Filling atom " << **iter << ", translated to " << AtomTranslations << ", at final position is " << (CopyAtoms[(*iter)->getNr()]->getPosition()) << ".");
    911           } else {
    912             LOG(1, "INFO: Position at " << Inserter << " is inner point, within boundary or outside of MaxDistance.");
    913             CopyAtoms[(*iter)->getNr()] = NULL;
    914             continue;
    915           }
    916         }
    917         // go through all bonds and add as well
    918         for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
    919           const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    920           for(BondList::const_iterator BondRunner = ListOfBonds.begin();
    921               BondRunner != ListOfBonds.end();
    922               ++BondRunner)
    923             if ((*BondRunner)->leftatom == *AtomRunner) {
    924               Binder = (*BondRunner);
    925               if ((CopyAtoms[Binder->leftatom->getNr()] != NULL) && (CopyAtoms[Binder->rightatom->getNr()] != NULL)) {
    926                 LOG(3, "Adding Bond between " << *CopyAtoms[Binder->leftatom->getNr()] << " and " << *CopyAtoms[Binder->rightatom->getNr()]<< ".");
    927                 Filling->AddBond(CopyAtoms[Binder->leftatom->getNr()], CopyAtoms[Binder->rightatom->getNr()], Binder->getDegree());
    928               }
    929             }
    930         }
    931       }
    932   for (MoleculeList::iterator ListRunner = List->ListOfMolecules.begin(); ListRunner != List->ListOfMolecules.end(); ListRunner++) {
    933     delete LCList[*ListRunner];
    934     delete TesselStruct[(*ListRunner)];
    935   }
    936 };
    937 
    938 /** Rotates given molecule \a Filling and moves its atoms according to given
    939  *  \a RandomAtomDisplacement.
    940  *
    941  *  Note that for rotation to be sensible, the molecule should be centered at
    942  *  the origin. This is not done here!
    943  *
    944  *  \param &Filling molecule whose atoms to displace
    945  *  \param RandomAtomDisplacement magnitude of random displacement
    946  *  \param &Rotations 3D rotation matrix (or unity if no rotation desired)
    947  */
    948 void RandomizeMoleculePositions(
    949     molecule *&Filling,
    950     double RandomAtomDisplacement,
    951     RealSpaceMatrix &Rotations,
    952     RandomNumberGenerator &random
    953     )
    954 {
    955   const double rng_min = random.min();
    956   const double rng_max = random.max();
    957 
    958   Vector AtomTranslations;
    959   for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) {
    960     Vector temp = (*miter)->getPosition();
    961     temp *= Rotations;
    962     (*miter)->setPosition(temp);
    963     // create atomic random translation vector ...
    964     for (int i=0;i<NDIM;i++)
    965       AtomTranslations[i] = RandomAtomDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
    966     (*miter)->setPosition((*miter)->getPosition() + AtomTranslations);
    967   }
    968 }
    969 
    970 /** Removes all atoms of a molecule outside.
    971  *
    972  * If the molecule is empty, it is removed as well.
    973  *
    974  * @param Filling molecule whose atoms to check, removed if eventually left
    975  *        empty.
    976  * @return true - atoms had to be removed, false - no atoms have been removed
    977  */
    978 bool RemoveAtomsOutsideDomain(molecule *&Filling)
    979 {
    980   bool status = false;
    981   Box &Domain = World::getInstance().getDomain();
    982   // check if all is still inside domain
    983   for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {
    984     // check whether each atom is inside box
    985     if (!Domain.isInside((*miter)->getPosition())) {
    986       status = true;
    987       atom *Walker = *miter;
    988       ++miter;
    989       World::getInstance().destroyAtom(Walker);
    990     } else {
    991       ++miter;
    992     }
    993   }
    994   if (Filling->empty()) {
    995     LOG(0, "Removing molecule " << Filling->getName() << ", all atoms have been removed.");
    996     World::getInstance().destroyMolecule(Filling);
    997     Filling = NULL;
    998   }
    999   return status;
    1000 }
    1001 
    1002 /** Checks whether there are no atoms inside a sphere around \a CurrentPosition
    1003  *  except those atoms present in \a *filler.
    1004  *  If filler is NULL, then we just call LinkedCell_deprecated::GetPointsInsideSphere() and
    1005  *  check whether the return list is empty.
    1006  * @param *filler
    1007  * @param boundary
    1008  * @param CurrentPosition
    1009  */
    1010 bool isSpaceAroundPointVoid(
    1011     LinkedCell_deprecated *LC,
    1012     molecule *filler,
    1013     const double boundary,
    1014     Vector &CurrentPosition)
    1015 {
    1016   size_t compareTo = 0;
    1017   TesselPointSTLList* liste = LC->GetPointsInsideSphere(boundary == 0. ? MYEPSILON : boundary, &CurrentPosition);
    1018   if (filler != NULL) {
    1019     for (TesselPointSTLList::const_iterator iter = liste->begin();
    1020         iter != liste->end();
    1021         ++iter) {
    1022       for (molecule::iterator miter = filler->begin();
    1023           miter != filler->end();
    1024           ++miter) {
    1025         if (*iter == *miter)
    1026           ++compareTo;
    1027       }
    1028     }
    1029   }
    1030   const bool result = (liste->size() == compareTo);
    1031   if (!result) {
    1032     LOG(0, "Skipping because of the following atoms:");
    1033     for (TesselPointSTLList::const_iterator iter = liste->begin();
    1034         iter != liste->end();
    1035         ++iter) {
    1036       LOG(0, **iter);
    1037     }
    1038   }
    1039   delete(liste);
    1040   return result;
    1041 }
    1042 
    1043 /** Sets given 3x3 matrix to a random rotation matrix.
    1044  *
    1045  * @param a matrix to set
    1046  */
    1047 inline void setRandomRotation(RealSpaceMatrix &a)
    1048 {
    1049   double phi[NDIM];
    1050   RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
    1051   const double rng_min = random.min();
    1052   const double rng_max = random.max();
    1053 
    1054   for (int i=0;i<NDIM;i++) {
    1055     phi[i] = (random()/(rng_max-rng_min))*(2.*M_PI);
    1056     LOG(4, "DEBUG: Random angle is " << phi[i] << ".");
    1057   }
    1058 
    1059   a.setRotation(phi);
    1060 }
    1061 
    1062 /** Fills the empty space of the simulation box with water.
    1063  * \param *filler molecule which the box is to be filled with
    1064  * \param configuration contains box dimensions
    1065  * \param distance[NDIM] distance between filling molecules in each direction
    1066  * \param boundary length of boundary zone between molecule and filling molecules
    1067  * \param RandAtomDisplacement maximum distance for random displacement per atom
    1068  * \param RandMolDisplacement maximum distance for random displacement per filler molecule
    1069  * \param MinDistance minimum distance to boundary of domain and present molecules
    1070  * \param DoRandomRotation true - do random rotations, false - don't
    1071  */
    1072 void FillVoidWithMolecule(
    1073     molecule *&filler,
    1074     config &configuration,
    1075     const double distance[NDIM],
    1076     const double boundary,
    1077     const double RandomAtomDisplacement,
    1078     const double RandomMolDisplacement,
    1079     const double MinDistance,
    1080     const bool DoRandomRotation
    1081     )
    1082 {
    1083   //Info FunctionInfo(__func__);
    1084   molecule *Filling = NULL;
    1085   Vector CurrentPosition;
    1086   int N[NDIM];
    1087   int n[NDIM];
    1088   const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
    1089   RealSpaceMatrix Rotations;
    1090   const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
    1091   Vector FillerTranslations;
    1092   Vector FillerDistance;
    1093   Vector Inserter;
    1094   double FillIt = false;
    1095   Vector firstInserter;
    1096   bool firstInsertion = true;
    1097   const Box &Domain = World::getInstance().getDomain();
    1098   map<molecule *, LinkedCell_deprecated *> LCList;
    1099   std::vector<molecule *> List = World::getInstance().getAllMolecules();
    1100   MoleculeListClass *MolList = World::getInstance().getMolecules();
    1101 
    1102   for (std::vector<molecule *>::iterator ListRunner = List.begin(); ListRunner != List.end(); ListRunner++)
    1103     if ((*ListRunner)->getAtomCount() > 0) {
    1104       LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
    1105       PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
    1106       LCList[(*ListRunner)] = new LinkedCell_deprecated(cloud, 10.); // get linked cell list
    1107     }
    1108 
    1109   // Center filler at its center of gravity
    1110   const Vector gravity = filler->DetermineCenterOfGravity();
    1111   filler->CenterAtVector(gravity);
    1112   //const int FillerCount = filler->getAtomCount();
    1113   LOG(2, "INFO: Filler molecule has the following bonds:");
    1114   for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
    1115     const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    1116     for(BondList::const_iterator BondRunner = ListOfBonds.begin();
    1117         BondRunner != ListOfBonds.end();
    1118         ++BondRunner)
    1119       if ((*BondRunner)->leftatom == *AtomRunner)
    1120         LOG(2, "  " << *(*BondRunner));
    1121   }
    1122 
    1123   // calculate filler grid in [0,1]^3
    1124   FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
    1125   for(int i=0;i<NDIM;i++)
    1126     N[i] = (int) ceil(1./FillerDistance[i]);
    1127   LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");
    1128 
    1129   // initialize seed of random number generator to current time
    1130   RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
    1131   const double rng_min = random.min();
    1132   const double rng_max = random.max();
    1133   //srand ( time(NULL) );
    1134 
    1135   // go over [0,1]^3 filler grid
    1136   for (n[0] = 0; n[0] < N[0]; n[0]++)
    1137     for (n[1] = 0; n[1] < N[1]; n[1]++)
    1138       for (n[2] = 0; n[2] < N[2]; n[2]++) {
    1139         // calculate position of current grid vector in untransformed box
    1140         CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    1141         // create molecule random translation vector ...
    1142         for (int i=0;i<NDIM;i++) // have the random values [-1,1]*RandomMolDisplacement
    1143           FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
    1144         LOG(2, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");
    1145 
    1146         // ... and rotation matrix
    1147         if (DoRandomRotation)
    1148           setRandomRotation(Rotations);
    1149         else
    1150           Rotations.setIdentity();
    1151 
    1152 
    1153         // Check whether there is anything too close by and whether atom is outside of domain
    1154         FillIt = true;
    1155         for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
    1156           FillIt = FillIt && isSpaceAroundPointVoid(
    1157               ListRunner->second,
    1158               (firstInsertion ? filler : NULL),
    1159               boundary,
    1160               CurrentPosition);
    1161           FillIt = FillIt && (Domain.isValid(CurrentPosition))
    1162               && ((Domain.DistanceToBoundary(CurrentPosition) - MinDistance) > -MYEPSILON);
    1163           if (!FillIt)
    1164             break;
    1165         }
    1166 
    1167         // insert into Filling
    1168         if (FillIt) {
    1169           Inserter = CurrentPosition + FillerTranslations;
    1170           LOG(1, "INFO: Position at " << Inserter << " is void point.");
    1171           // fill!
    1172           Filling = filler->CopyMolecule();
    1173           RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations, random);
    1174           // translation
    1175           Filling->Translate(Inserter);
    1176           // remove out-of-bounds atoms
    1177           const bool status = RemoveAtomsOutsideDomain(Filling);
    1178           if ((firstInsertion) && (!status)) { // no atom has been removed
    1179             // remove copied atoms and molecule again
    1180             Filling->removeAtomsinMolecule();
    1181             World::getInstance().destroyMolecule(Filling);
    1182             // and mark is final filler position
    1183             Filling = filler;
    1184             firstInsertion = false;
    1185             firstInserter = Inserter;
    1186           } else {
    1187             // TODO: Remove when World has no MoleculeListClass anymore
    1188             if (Filling)
    1189               MolList->insert(Filling);
    1190           }
    1191         } else {
    1192          LOG(1, "INFO: Position at " << Inserter << " is non-void point, within boundary or outside of MaxDistance.");
    1193           continue;
    1194         }
    1195       }
    1196 
    1197   // have we inserted any molecules?
    1198   if (firstInsertion) {
    1199     // If not remove filler
    1200     for(molecule::iterator miter = filler->begin(); !filler->empty(); miter = filler->begin()) {
    1201       atom *Walker = *miter;
    1202       World::getInstance().destroyAtom(Walker);
    1203     }
    1204     World::getInstance().destroyMolecule(filler);
    1205   } else {
    1206     // otherwise translate and randomize to final position
    1207     if (DoRandomRotation)
    1208       setRandomRotation(Rotations);
    1209     else
    1210       Rotations.setIdentity();
    1211     RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations, random);
    1212     // translation
    1213     filler->Translate(firstInserter);
    1214     // remove out-of-bounds atoms
    1215     RemoveAtomsOutsideDomain(filler);
    1216   }
    1217 
    1218   LOG(0, MolList->ListOfMolecules.size() << " molecules have been inserted.");
    1219 
    1220   for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
    1221     delete ListRunner->second;
    1222     LCList.erase(ListRunner);
    1223   }
    1224 };
    1225 
    1226 
    1227 /** Fills the empty space around other molecules' surface of the simulation box with a filler.
    1228  *
    1229  * Note that we use \a FindNonConvexBorder to determine the surface of the found molecules.
    1230  * There, we use a radius of twice the given \a boundary.
    1231  *
    1232  * \param *out output stream for debugging
    1233  * \param *List list of molecules already present in box
    1234  * \param *TesselStruct contains tesselated surface
    1235  * \param *filler molecule which the box is to be filled with
    1236  * \param configuration contains box dimensions
    1237  * \param MaxSurfaceDistance fills in molecules only up to this distance (set to -1 if whole of the domain)
    1238  * \param distance[NDIM] distance between filling molecules in each direction
    1239  * \param boundary length of boundary zone between molecule and filling molecules
    1240  * \param MinDistance distance to boundary of domain which is not filled
    1241  * \param RandAtomDisplacement maximum distance for random displacement per atom
    1242  * \param RandMolDisplacement maximum distance for random displacement per filler molecule
    1243  * \param DoRandomRotation true - do random rotations, false - don't
    1244  */
    1245 void FillBoxWithMolecule(
    1246     MoleculeListClass *MolList,
    1247     molecule *filler,
    1248     config &configuration,
    1249     const double MaxSurfaceDistance,
    1250     const double distance[NDIM],
    1251     const double boundary,
    1252     const double MinDistance,
    1253     const double RandomAtomDisplacement,
    1254     const double RandomMolDisplacement,
    1255     const bool DoRandomRotation)
    1256 {
    1257   //Info FunctionInfo(__func__);
    1258   molecule *Filling = NULL;
    1259   Vector CurrentPosition;
    1260   int N[NDIM];
    1261   int n[NDIM];
    1262   const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
    1263   RealSpaceMatrix Rotations;
    1264   const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
    1265   Vector FillerTranslations;
    1266   Vector FillerDistance;
    1267   Vector Inserter;
    1268   double FillIt = false;
    1269   Vector firstInserter;
    1270   bool firstInsertion = true;
    1271   const Box &Domain = World::getInstance().getDomain();
    1272   map<molecule *, LinkedCell_deprecated *> LCList;
    1273   std::vector<molecule *> List = World::getInstance().getAllMolecules();
    1274   map<molecule *, Tesselation *> TesselStruct;
    1275 
    1276   for (MoleculeList::iterator ListRunner = MolList->ListOfMolecules.begin(); ListRunner != MolList->ListOfMolecules.end(); ListRunner++)
    1277     if ((*ListRunner)->getAtomCount() > 0) {
    1278       LOG(1, "Pre-creating linked cell lists for molecule " << *ListRunner << ".");
    1279       PointCloudAdaptor< molecule > cloud(*ListRunner, (*ListRunner)->name);
    1280       LCList[(*ListRunner)] = new LinkedCell_deprecated(cloud, 4.*boundary); // get linked cell list
    1281       LOG(1, "Pre-creating tesselation for molecule " << *ListRunner << ".");
    1282       TesselStruct[(*ListRunner)] = NULL;
    1283       FindNonConvexBorder((*ListRunner), TesselStruct[(*ListRunner)], (const LinkedCell_deprecated *&)LCList[(*ListRunner)], 2.*boundary, NULL);
    1284     }
    1285 
    1286   // Center filler at origin
    1287   filler->CenterEdge();
    1288 //  const int FillerCount = filler->getAtomCount();
    1289   LOG(2, "INFO: Filler molecule has the following bonds:");
    1290   for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) {
    1291     const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    1292     for(BondList::const_iterator BondRunner = ListOfBonds.begin();
    1293         BondRunner != ListOfBonds.end();
    1294         ++BondRunner) {
    1295       if ((*BondRunner)->leftatom == *AtomRunner)
    1296         LOG(2, "  " << *(*BondRunner));
    1297     }
    1298   }
    1299 
    1300 //  atom * CopyAtoms[FillerCount];
    1301 
    1302   setVerbosity(4);
    1303 
    1304   // calculate filler grid in [0,1]^3
    1305   FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
    1306   for(int i=0;i<NDIM;i++)
    1307     N[i] = (int) ceil(1./FillerDistance[i]);
    1308   LOG(1, "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << ".");
    1309 
    1310   // initialize seed of random number generator to current time
    1311   RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
    1312   const double rng_min = random.min();
    1313   const double rng_max = random.max();
    1314   //srand ( time(NULL) );
    1315 
    1316   // go over [0,1]^3 filler grid
    1317   for (n[0] = 0; n[0] < N[0]; n[0]++)
    1318     for (n[1] = 0; n[1] < N[1]; n[1]++)
    1319       for (n[2] = 0; n[2] < N[2]; n[2]++) {
    1320         // calculate position of current grid vector in untransformed box
    1321         CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    1322         // create molecule random translation vector ...
    1323         for (int i=0;i<NDIM;i++)
    1324           FillerTranslations[i] = RandomMolDisplacement*(random()/((rng_max-rng_min)/2.) - 1.);
    1325         LOG(1, "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << ".");
    1326 
    1327         // ... and rotation matrix
    1328         if (DoRandomRotation)
    1329           setRandomRotation(Rotations);
    1330         else
    1331           Rotations.setIdentity();
    1332 
    1333 
    1334         // Check whether there is anything too close by and whether atom is outside of domain
    1335         FillIt = true;
    1336         for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
    1337           // check whether its void
    1338           FillIt = FillIt && isSpaceAroundPointVoid(
    1339               ListRunner->second,
    1340               (firstInsertion ? filler : NULL),
    1341               boundary,
    1342               CurrentPosition);
    1343           if (!FillIt) {
    1344             LOG(2, "REJECT: Position at " << Inserter << " is non-void.");
    1345             break;
    1346           }
    1347           // check whether inside domain
    1348           FillIt = FillIt && (Domain.isValid(CurrentPosition));
    1349           if (!FillIt) {
    1350             LOG(2, "REJECT: Position at " << CurrentPosition << " is "
    1351                 << distance << ", hence outside of domain.");
    1352             break;
    1353           }
    1354           // check minimum distance to boundary
    1355           const double distance = (Domain.DistanceToBoundary(CurrentPosition) - MinDistance);
    1356           FillIt = FillIt && (distance > -MYEPSILON);
    1357           if (!FillIt) {
    1358             LOG(2, "REJECT: Position at " << CurrentPosition << " is " << distance << ", less than "
    1359                 << MinDistance << " hence, too close to boundary.");
    1360             break;
    1361           }
    1362         }
    1363         // Check whether point is in- or outside of tesselations
    1364         if (FillIt) {
    1365           for (MoleculeList::iterator ListRunner = MolList->ListOfMolecules.begin(); ListRunner != MolList->ListOfMolecules.end(); ListRunner++) {
    1366             // get linked cell list
    1367             if (TesselStruct[(*ListRunner)] != NULL) {
    1368               const double distance = (TesselStruct[(*ListRunner)]->GetDistanceToSurface(Inserter, LCList[(*ListRunner)]));
    1369               LOG(2, "INFO: Distance to surface is " << distance << ".");
    1370               FillIt = FillIt && ((distance == -1.) || (distance > boundary)) && ((MaxSurfaceDistance < 0) || (MaxSurfaceDistance > distance));
    1371               if (!FillIt) {
    1372                 LOG(2, "REJECT: Position at " << CurrentPosition << " is in distance of " << distance
    1373                     << " to a surface, less than " << MaxSurfaceDistance  << " hence, too close.");
    1374                 break;
    1375               }
    1376             }
    1377           }
    1378         }
    1379 
    1380         // insert into Filling
    1381         if (FillIt) {
    1382           Inserter = CurrentPosition + FillerTranslations;
    1383           LOG(2, "ACCEPT: Position at " << CurrentPosition << " is void point.");
    1384           // fill!
    1385           Filling = filler->CopyMolecule();
    1386           RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations, random);
    1387           // translation
    1388           Filling->Translate(Inserter);
    1389           // remove out-of-bounds atoms
    1390           const bool status = RemoveAtomsOutsideDomain(Filling);
    1391           if ((firstInsertion) && (!status)) { // no atom has been removed
    1392             // remove copied atoms and molecule again
    1393             Filling->removeAtomsinMolecule();
    1394             World::getInstance().destroyMolecule(Filling);
    1395             // and mark is final filler position
    1396             Filling = filler;
    1397             firstInsertion = false;
    1398             firstInserter = Inserter;
    1399           } else {
    1400             // TODO: Remove when World has no MoleculeListClass anymore
    1401             if (Filling)
    1402               MolList->insert(Filling);
    1403           }
    1404         } else {
    1405           LOG(2, "REJECT: Position at " << CurrentPosition << " is non-void point, within boundary or outside of MaxSurfaceDistance.");
    1406           continue;
    1407         }
    1408       }
    1409 
    1410   // have we inserted any molecules?
    1411   if (firstInsertion) {
    1412     // If not remove filler
    1413     for(molecule::iterator miter = filler->begin(); !filler->empty(); miter = filler->begin()) {
    1414       atom *Walker = *miter;
    1415       World::getInstance().destroyAtom(Walker);
    1416     }
    1417     World::getInstance().destroyMolecule(filler);
    1418   } else {
    1419     // otherwise translate and randomize to final position
    1420     if (DoRandomRotation)
    1421       setRandomRotation(Rotations);
    1422     else
    1423       Rotations.setIdentity();
    1424     RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations, random);
    1425     // translation
    1426     filler->Translate(firstInserter);
    1427     // remove out-of-bounds atoms
    1428     RemoveAtomsOutsideDomain(filler);
    1429   }
    1430 
    1431   LOG(0, MolList->ListOfMolecules.size() << " molecules have been inserted.");
    1432 
    1433   for (std::map<molecule *, LinkedCell_deprecated *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
    1434     delete ListRunner->second;
    1435     LCList.erase(ListRunner);
    1436   }
    1437 };
    1438 
    1439769/** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.
    1440770 * \param *out output stream for debugging

  • src/Tesselation/boundary.hpp

    r16c7dd r2528d8  
    4040
    4141double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
    42 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);
    43 void FillVoidWithMolecule(molecule *&filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const double MinDistance, const bool DoRandomRotation);
    4442void FindConvexBorder(const molecule* const mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell_deprecated *LCList, const char *filename);
    4543Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);
Note: See TracChangeset for help on using the changeset viewer.