- Timestamp:
- Feb 3, 2011, 9:51:18 AM (14 years ago)
- 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:
- d6f886
- Parents:
- 0b15bb
- git-author:
- Frederik Heber <heber@…> (12/30/10 20:52:17)
- git-committer:
- Frederik Heber <heber@…> (02/03/11 09:51:18)
- Location:
- src
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
src/Actions/MoleculeAction/FillVoidWithMoleculeAction.cpp
r0b15bb r66fd49 28 28 #include "World.hpp" 29 29 30 #include "Descriptors/MoleculeIdDescriptor.hpp" 30 31 #include "Parser/MpqcParser.hpp" 31 32 #include "Parser/PcpParser.hpp" … … 35 36 #include "Parser/FormatParserStorage.hpp" 36 37 38 #include <algorithm> 37 39 #include <iostream> 38 40 #include <string> … … 47 49 /** =========== define the function ====================== */ 48 50 Action::state_ptr MoleculeFillVoidWithMoleculeAction::performCall() { 49 50 51 // obtain information 51 52 getParametersfromValueStorage(); 52 53 53 DoLog(1) && (Log() << Verbose(1) << "Filling Box with water molecules, lengths(" << params.lengths[0] << "," << params.lengths[1] << "," << params.lengths[2] << "), distances (" << params.distances[0] << "," << params.distances[1] << "," << params.distances[2] << "), DoRotate " << params.DoRotate << "." << endl); 54 if (!boost::filesystem::exists(params.fillername)) { 55 DoeLog(1) && (eLog() << Verbose(1) << "File with filler molecule " << params.fillername << " does not exist!" << endl); 56 return Action::failure; 57 } 58 59 DoLog(1) && (Log() << Verbose(1) << "Filling Box with water molecules, lengths(" << params.lengths[0] << "," << params.lengths[1] << "," << params.lengths[2] << "), distances (" << params.distances[0] << "," << params.distances[1] << "," << params.distances[2] << "), MinDistance " << params.MinDistance << ", DoRotate " << params.DoRotate << "." << endl); 54 60 // construct water molecule 55 molecule *filler = World::getInstance().createMolecule(); 61 std::vector<molecule *> presentmolecules = World::getInstance().getAllMolecules(); 62 // DoLog(0) && (Log() << Verbose(0) << presentmolecules.size() << " molecules initially are present." << std::endl); 56 63 std::string FilenameSuffix = params.fillername.string().substr(params.fillername.string().find_last_of('.')+1, params.fillername.string().length()); 57 64 ifstream input; … … 92 99 break; 93 100 } 94 World::MoleculeIterator iter = World::getInstance().getMoleculeIter(); 95 for (; iter != World::getInstance().moleculeEnd(); ++iter) 101 102 // search the filler molecule that has been just parsed 103 molecule *filler = NULL; 104 for (World::MoleculeIterator iter = World::getInstance().getMoleculeIter(); 105 iter != World::getInstance().moleculeEnd(); 106 ++iter) 96 107 filler = *iter; // get last molecule 97 108 filler->SetNameFromFilename(params.fillername.string().c_str()); … … 102 113 for (int i=0;i<NDIM;i++) 103 114 distance[i] = params.distances[i]; 104 FillVoidWithMolecule(filler, *(World::getInstance().getConfig()), distance, params.lengths[0], params.lengths[1], params.lengths[2], params. DoRotate);115 FillVoidWithMolecule(filler, *(World::getInstance().getConfig()), distance, params.lengths[0], params.lengths[1], params.lengths[2], params.MinDistance, params.DoRotate); 105 116 106 return Action::success; 117 // generate list of newly created molecules 118 // (we can in general remove more quickly from a list than a vector) 119 std::vector<molecule *> fillermolecules = World::getInstance().getAllMolecules(); 120 // DoLog(0) && (Log() << Verbose(0) << fillermolecules.size() << " molecules are present." << std::endl); 121 std::list<molecule *> fillermolecules_list; 122 std::copy( fillermolecules.begin(), fillermolecules.end(), std::back_inserter( fillermolecules_list )); 123 // DoLog(0) && (Log() << Verbose(0) << fillermolecules_list.size() << " molecules have been copied." << std::endl); 124 for (std::vector<molecule *>::const_iterator iter = presentmolecules.begin(); 125 iter != presentmolecules.end(); 126 ++iter) { 127 fillermolecules_list.remove(*iter); 128 } 129 // DoLog(0) && (Log() << Verbose(0) << fillermolecules_list.size() << " molecules left after removal." << std::endl); 130 fillermolecules.clear(); 131 std::copy(fillermolecules_list.begin(), fillermolecules_list.end(), std::back_inserter( fillermolecules )); 132 133 // DoLog(0) && (Log() << Verbose(0) << fillermolecules.size() << " molecules have been inserted." << std::endl); 134 135 return Action::state_ptr(new MoleculeFillVoidWithMoleculeState(fillermolecules,params)); 107 136 } 108 137 109 138 Action::state_ptr MoleculeFillVoidWithMoleculeAction::performUndo(Action::state_ptr _state) { 110 //MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get());139 MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get()); 111 140 112 // string newName = state->mol->getName(); 113 // state->mol->setName(state->lastName); 141 MoleculeListClass *MolList = World::getInstance().getMolecules(); 114 142 115 return Action::failure; 143 BOOST_FOREACH(molecule *_mol, state->fillermolecules) { 144 MolList->erase(_mol); 145 if ((_mol != NULL) && (!(World::getInstance().getAllMolecules(MoleculeById(_mol->getId()))).empty())) { 146 for(molecule::iterator iter = _mol->begin(); 147 !_mol->empty(); 148 iter = _mol->begin()) { 149 atom *Walker = *iter; 150 _mol->erase(iter); 151 World::getInstance().destroyAtom(Walker); 152 } 153 World::getInstance().destroyMolecule(_mol); 154 } 155 } 156 157 // as molecules and atoms from state are removed, we have to create a new one 158 std::vector<molecule *> fillermolecules; 159 return Action::state_ptr(new MoleculeFillVoidWithMoleculeState(fillermolecules,state->params)); 116 160 } 117 161 118 162 Action::state_ptr MoleculeFillVoidWithMoleculeAction::performRedo(Action::state_ptr _state){ 119 // Undo and redo have to do the same for this action 120 return performUndo(_state); 163 //MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get()); 164 165 return Action::failure; 166 //return Action::state_ptr(_state); 121 167 } 122 168 123 169 bool MoleculeFillVoidWithMoleculeAction::canUndo() { 124 return false;170 return true; 125 171 } 126 172 127 173 bool MoleculeFillVoidWithMoleculeAction::shouldUndo() { 128 return false;174 return true; 129 175 } 130 176 /** =========== end of function ====================== */ -
src/Actions/MoleculeAction/FillVoidWithMoleculeAction.def
r0b15bb r66fd49 13 13 // ValueStorage by the token "Z" -> first column: int, Z, "Z" 14 14 // "undefine" if no parameters are required, use (NODEFAULT) for each (undefined) default value 15 #define paramtypes (boost::filesystem::path)(Vector)(Vector)( bool)16 #define paramtokens ("fill-void")("distances")("lengths")(" DoRotate")17 #define paramdescriptions ("name of xyz file of filler molecule")("list of three of distances in space, one for each axis direction")("list of three of lengths in space, one for each axis direction")(" whether to rotate or not")18 # undef paramdefaults19 #define paramreferences (fillername)(distances)(lengths)( DoRotate)15 #define paramtypes (boost::filesystem::path)(Vector)(Vector)(double)(bool) 16 #define paramtokens ("fill-void")("distances")("lengths")("MinDistance")("DoRotate") 17 #define paramdescriptions ("name of xyz file of filler molecule")("list of three of distances in space, one for each axis direction")("list of three of lengths in space, one for each axis direction")("minimum distance to boundary")("whether to rotate or not") 18 #define paramdefaults (NODEFAULT)(NODEFAULT)(NODEFAULT)("0.")("0") 19 #define paramreferences (fillername)(distances)(lengths)(MinDistance)(DoRotate) 20 20 21 # undef statetypes22 # undef statereferences21 #define statetypes (std::vector<molecule *>) 22 #define statereferences (fillermolecules) 23 23 24 24 // some defines for all the names, you may use ACTION, STATE and PARAMS -
src/Box.cpp
r0b15bb r66fd49 33 33 34 34 #include "CodePatterns/Assert.hpp" 35 #include "CodePatterns/Log.hpp" 36 #include "CodePatterns/Verbose.hpp" 35 37 36 38 using namespace std; … … 220 222 } 221 223 224 double Box::DistanceToBoundary(const Vector &point) const 225 { 226 std::map<double, Plane> DistanceSet; 227 std::vector<std::pair<Plane,Plane> > Boundaries = getBoundingPlanes(); 228 for (int i=0;i<NDIM;++i) { 229 const double tempres1 = Boundaries[i].first.distance(point); 230 const double tempres2 = Boundaries[i].second.distance(point); 231 DistanceSet.insert( make_pair(tempres1, Boundaries[i].first) ); 232 DoLog(1) && (Log() << Verbose(1) << "Inserting distance " << tempres1 << " and " << tempres2 << "." << std::endl); 233 DistanceSet.insert( make_pair(tempres2, Boundaries[i].second) ); 234 } 235 ASSERT(!DistanceSet.empty(), "Box::DistanceToBoundary() - no distances in map!"); 236 return (DistanceSet.begin())->first; 237 } 238 222 239 Shape Box::getShape() const{ 223 240 return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M)); 224 241 } 225 242 226 const Box::Conditions_t Box::getConditions(){ 243 const Box::Conditions_t Box::getConditions() const 244 { 227 245 return conditions; 228 246 } … … 232 250 } 233 251 234 const vector<pair<Plane,Plane> > Box::getBoundingPlanes(){ 252 const vector<pair<Plane,Plane> > Box::getBoundingPlanes() const 253 { 235 254 vector<pair<Plane,Plane> > res; 236 255 for(int i=0;i<NDIM;++i){ … … 248 267 res.push_back(make_pair(p1,p2)); 249 268 } 269 ASSERT(res.size() == 3, "Box::getBoundingPlanes() - does not have three plane pairs!"); 250 270 return res; 251 271 } -
src/Box.hpp
r0b15bb r66fd49 106 106 double periodicDistance(const Vector &point1,const Vector &point2) const; 107 107 108 /** 109 * Calculates the minimum distance to the boundary of the periodic 110 * space defined by this box. 111 */ 112 double DistanceToBoundary(const Vector &point) const; 113 108 114 Shape getShape() const; 109 const Conditions_t getConditions() ;115 const Conditions_t getConditions() const; 110 116 void setCondition(int,BoundaryCondition_t); 111 117 112 const vector<pair<Plane,Plane> > getBoundingPlanes() ;118 const vector<pair<Plane,Plane> > getBoundingPlanes() const; 113 119 114 120 void setCuboid(const Vector&); -
src/LinearAlgebra/RealSpaceMatrix.cpp
r0b15bb r66fd49 121 121 set(2,2, cos(x)*cos(y)); 122 122 } 123 124 void RealSpaceMatrix::setRandomRotation() 125 { 126 double phi[NDIM]; 127 128 for (int i=0;i<NDIM;i++) { 129 phi[i] = rand()/(RAND_MAX/(2.*M_PI)); 130 } 131 132 set(0,0, cos(phi[0]) *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]))); 133 set(0,1, sin(phi[0]) *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]))); 134 set(0,2, cos(phi[1])*sin(phi[2]) ); 135 set(1,0, -sin(phi[0])*cos(phi[1]) ); 136 set(1,1, cos(phi[0])*cos(phi[1]) ); 137 set(1,2, sin(phi[1]) ); 138 set(2,0, -cos(phi[0]) *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]))); 139 set(2,1, -sin(phi[0]) *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]))); 140 set(2,2, cos(phi[1])*cos(phi[2]) ); 141 } 142 123 143 124 144 RealSpaceMatrix &RealSpaceMatrix::operator=(const RealSpaceMatrix &src) -
src/LinearAlgebra/RealSpaceMatrix.hpp
r0b15bb r66fd49 72 72 */ 73 73 void setRotation(const double x, const double y, const double z); 74 void setRandomRotation(); 74 75 75 76 /** -
src/boundary.cpp
r0b15bb r66fd49 927 927 }; 928 928 929 /** Fills the empty space of the simulation box with water/ 930 * \param *out output stream for debugging 931 * \param *TesselStruct contains tesselated surface 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 ) 944 { 945 Vector AtomTranslations; 946 for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) { 947 Vector temp = (*miter)->getPosition(); 948 temp *= Rotations; 949 (*miter)->setPosition(temp); 950 // create atomic random translation vector ... 951 for (int i=0;i<NDIM;i++) 952 AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.); 953 (*miter)->setPosition((*miter)->getPosition() + AtomTranslations); 954 } 955 } 956 957 /** Removes all atoms of a molecule outside. 958 * 959 * If the molecule is empty, it is removed as well. 960 * 961 * @param Filling molecule whose atoms to check, removed if eventually left 962 * empty. 963 */ 964 void RemoveAtomsOutsideDomain(molecule *&Filling) 965 { 966 Box &Domain = World::getInstance().getDomain(); 967 // check if all is still inside domain 968 for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) { 969 // check whether each atom is inside box 970 if (!Domain.isInside((*miter)->getPosition())) { 971 atom *Walker = *miter; 972 ++miter; 973 World::getInstance().destroyAtom(Walker); 974 } else { 975 ++miter; 976 } 977 } 978 if (Filling->empty()) { 979 DoLog(0) && (Log() << Verbose(0) << "Removing molecule " << Filling->getName() << ", all atoms have been removed." << std::endl); 980 World::getInstance().destroyMolecule(Filling); 981 } 982 } 983 984 /** Checks whether there are no atoms inside a sphere around \a CurrentPosition 985 * except those atoms present in \a *filler. 986 * If filler is NULL, then we just call LinkedCell::GetPointsInsideSphere() and 987 * check whether the return list is empty. 988 * @param *filler 989 * @param boundary 990 * @param CurrentPosition 991 */ 992 bool isSpaceAroundPointVoid( 993 LinkedCell *LC, 994 molecule *filler, 995 const double boundary, 996 Vector &CurrentPosition) 997 { 998 size_t compareTo = 0; 999 LinkedCell::LinkedNodes* liste = LC->GetPointsInsideSphere(boundary == 0. ? MYEPSILON : boundary, &CurrentPosition); 1000 if (filler != NULL) { 1001 for (LinkedCell::LinkedNodes::const_iterator iter = liste->begin(); 1002 iter != liste->end(); 1003 ++iter) { 1004 for (molecule::iterator miter = filler->begin(); 1005 miter != filler->end(); 1006 ++miter) { 1007 if (*iter == *miter) 1008 ++compareTo; 1009 } 1010 } 1011 } 1012 const bool result = (liste->size() == compareTo); 1013 if (!result) { 1014 DoLog(0) && (Log() << Verbose(0) << "Skipping because of the following atoms:" << std::endl); 1015 for (LinkedCell::LinkedNodes::const_iterator iter = liste->begin(); 1016 iter != liste->end(); 1017 ++iter) { 1018 DoLog(0) && (Log() << Verbose(0) << **iter << std::endl); 1019 } 1020 } 1021 delete(liste); 1022 return result; 1023 } 1024 1025 /** Fills the empty space of the simulation box with water. 932 1026 * \param *filler molecule which the box is to be filled with 933 1027 * \param configuration contains box dimensions 934 1028 * \param distance[NDIM] distance between filling molecules in each direction 935 * \param boundary length of boundary zone between molecule and filling mollecules 936 * \param epsilon distance to surface which is not filled 1029 * \param boundary length of boundary zone between molecule and filling molecules 937 1030 * \param RandAtomDisplacement maximum distance for random displacement per atom 938 1031 * \param RandMolDisplacement maximum distance for random displacement per filler molecule 939 * \param DoRandomRotation true - do random rotiations, false - don't 940 */ 941 void FillVoidWithMolecule(molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation) 1032 * \param MinDistance minimum distance to boundary of domain and present molecules 1033 * \param DoRandomRotation true - do random rotations, false - don't 1034 */ 1035 void FillVoidWithMolecule( 1036 molecule *&filler, 1037 config &configuration, 1038 const double distance[NDIM], 1039 const double boundary, 1040 const double RandomAtomDisplacement, 1041 const double RandomMolDisplacement, 1042 const double MinDistance, 1043 const bool DoRandomRotation 1044 ) 942 1045 { 943 1046 Info FunctionInfo(__func__); … … 949 1052 RealSpaceMatrix Rotations; 950 1053 const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv(); 951 Vector AtomTranslations;952 1054 Vector FillerTranslations; 953 1055 Vector FillerDistance; 954 1056 Vector Inserter; 955 1057 double FillIt = false; 956 bond *Binder = NULL; 957 double phi[NDIM]; 1058 Vector firstInserter; 1059 bool firstInsertion = true; 1060 const Box &Domain = World::getInstance().getDomain(); 958 1061 map<molecule *, LinkedCell *> LCList; 959 1062 std::vector<molecule *> List = World::getInstance().getAllMolecules(); … … 968 1071 // Center filler at origin 969 1072 filler->CenterEdge(&Inserter); 970 const int FillerCount = filler->getAtomCount();1073 //const int FillerCount = filler->getAtomCount(); 971 1074 DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl); 972 1075 for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner) … … 996 1099 997 1100 // ... and rotation matrix 998 if (DoRandomRotation) { 999 for (int i=0;i<NDIM;i++) { 1000 phi[i] = rand()/(RAND_MAX/(2.*M_PI)); 1001 } 1002 1003 Rotations.set(0,0, cos(phi[0]) *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]))); 1004 Rotations.set(0,1, sin(phi[0]) *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]))); 1005 Rotations.set(0,2, cos(phi[1])*sin(phi[2]) ); 1006 Rotations.set(1,0, -sin(phi[0])*cos(phi[1]) ); 1007 Rotations.set(1,1, cos(phi[0])*cos(phi[1]) ); 1008 Rotations.set(1,2, sin(phi[1]) ); 1009 Rotations.set(2,0, -cos(phi[0]) *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]))); 1010 Rotations.set(2,1, -sin(phi[0]) *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]))); 1011 Rotations.set(2,2, cos(phi[1])*cos(phi[2]) ); 1101 if (DoRandomRotation) 1102 Rotations.setRandomRotation(); 1103 else 1104 Rotations.setIdentity(); 1105 1106 1107 // Check whether there is anything too close by and whether atom is outside of domain 1108 FillIt = true; 1109 for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) { 1110 FillIt = FillIt && isSpaceAroundPointVoid( 1111 ListRunner->second, 1112 (firstInsertion ? filler : NULL), 1113 boundary, 1114 CurrentPosition); 1115 FillIt = FillIt && (Domain.isInside(CurrentPosition)) 1116 && ((Domain.DistanceToBoundary(CurrentPosition) - MinDistance) > -MYEPSILON); 1117 if (!FillIt) 1118 break; 1012 1119 } 1013 1120 1014 FillIt = true;1015 // go through all atoms1016 for(molecule::const_iterator iter = filler->begin(); iter !=filler->end();++iter){1017 1018 // Check whether there is anything too close by1019 LinkedCell::LinkedNodes* liste = NULL;1020 for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {1021 liste = ListRunner->second->GetPointsInsideSphere(boundary, &CurrentPosition);1022 FillIt = FillIt && (liste->size() == 0);1023 delete(liste);1024 if (!FillIt)1025 break;1026 }1027 }1028 1121 // insert into Filling 1029 1122 if (FillIt) { 1030 1123 Inserter = CurrentPosition + FillerTranslations; 1031 1124 DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is void point." << endl); 1032 Filling = filler->CopyMolecule(); 1033 for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) { 1034 if (DoRandomRotation) { 1035 Vector temp = (*miter)->getPosition(); 1036 temp *= Rotations; 1037 (*miter)->setPosition(temp); 1038 } 1039 // create atomic random translation vector ... 1040 for (int i=0;i<NDIM;i++) 1041 AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.); 1042 (*miter)->setPosition((*miter)->getPosition() + AtomTranslations); 1125 // fill! 1126 if (firstInsertion) { // use filler as first molecule 1127 Filling = filler; 1128 firstInsertion = false; 1129 firstInserter = Inserter; 1130 } else { // copy from filler molecule 1131 Filling = filler->CopyMolecule(); 1132 RandomizeMoleculePositions(Filling, RandomAtomDisplacement, Rotations); 1133 // translation 1134 Filling->Translate(&Inserter); 1135 // remove out-of-bounds atoms 1136 RemoveAtomsOutsideDomain(Filling); 1137 // TODO: Remove when World has no MoleculeListClass anymore 1138 MolList->insert(Filling); 1043 1139 } 1044 Filling->Translate(&Inserter);1045 for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {1046 // check whether each atom is inside box1047 if (!World::getInstance().getDomain().isInside((*miter)->getPosition())) {1048 atom *Walker = *miter;1049 ++miter;1050 World::getInstance().destroyAtom(Walker);1051 } else {1052 ++miter;1053 }1054 }1055 MolList->insert(Filling);1056 1140 } else { 1057 1141 DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is non-void point, within boundary or outside of MaxDistance." << endl); … … 1059 1143 } 1060 1144 } 1061 // last one is replaced by the filler, as we need the original atoms contained therein! 1062 filler->Translate(&Inserter); 1063 MolList->erase(Filling); 1064 for (molecule::iterator iter = Filling->begin(); !Filling->empty(); iter = Filling->begin()) { 1065 atom *Walker = *iter; 1066 Filling->erase(iter); 1067 World::getInstance().destroyAtom(Walker); 1068 } 1069 World::getInstance().destroyMolecule(Filling); 1145 1146 // have we inserted any molecules? 1147 if (firstInsertion) { 1148 // If not remove filler 1149 for(molecule::iterator miter = filler->begin(); !filler->empty(); miter = filler->begin()) { 1150 atom *Walker = *miter; 1151 filler->erase(Walker); 1152 World::getInstance().destroyAtom(Walker); 1153 } 1154 World::getInstance().destroyMolecule(filler); 1155 } else { 1156 // otherwise translate and randomize to final position 1157 if (DoRandomRotation) 1158 Rotations.setRandomRotation(); 1159 else 1160 Rotations.setIdentity(); 1161 RandomizeMoleculePositions(filler, RandomAtomDisplacement, Rotations); 1162 // translation 1163 filler->Translate(&firstInserter); 1164 // remove out-of-bounds atoms 1165 RemoveAtomsOutsideDomain(filler); 1166 } 1167 1168 DoLog(0) && (Log() << Verbose(0) << MolList->ListOfMolecules.size() << " molecules have been inserted." << std::endl); 1070 1169 1071 1170 for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) { -
src/boundary.hpp
r0b15bb r66fd49 41 41 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename); 42 42 molecule * 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 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); 44 44 void FindConvexBorder(const molecule* const mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename); 45 45 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);
Note:
See TracChangeset
for help on using the changeset viewer.