Changeset eee966

Timestamp:

Aug 28, 2010, 6:08:52 AM (14 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:

ef3640e

Parents:

a42054

Message:

New Action FillVoidWithMoleculeAction.

• in contrast to FillWithMoleculeAction we fill all empty space, not just outside of surfaces.
• TEST: new tests in Filling/3 that fill with butane and a tenside micelle.
• new function FillVoidWithMolecule().

Files:

• src/Actions/Makefile.am (modified) (2 diffs)
• src/Actions/MapOfActions.cpp (modified) (6 diffs)
• src/Actions/MoleculeAction/FillVoidWithMoleculeAction.cpp (added)
• src/Actions/MoleculeAction/FillVoidWithMoleculeAction.def (added)
• src/Actions/MoleculeAction/FillVoidWithMoleculeAction.hpp (added)
• src/boundary.cpp (modified) (2 diffs)
• src/boundary.hpp (modified) (1 diff)
• tests/regression/Filling/3/post/tensid.data (added)
• tests/regression/Filling/3/post/test.conf (added)
• tests/regression/Filling/3/post/test.xyz (added)
• tests/regression/Filling/3/pre/tensid.data (added)
• tests/regression/Filling/3/pre/test.conf (added)
• tests/regression/Filling/3/pre/water.data (added)
• tests/regression/testsuite-filling.at (modified) (2 diffs)

src/Actions/Makefile.am

@@ -91 +91 @@
   MoleculeAction/CopyAction.cpp \
   MoleculeAction/FillWithMoleculeAction.cpp \
+  MoleculeAction/FillVoidWithMoleculeAction.cpp \
   MoleculeAction/LinearInterpolationofTrajectoriesAction.cpp \
   MoleculeAction/RotateAroundSelfByAngleAction.cpp \
@@ -104 +105 @@
   MoleculeAction/CopyAction.hpp \
   MoleculeAction/FillWithMoleculeAction.hpp \
+  MoleculeAction/FillVoidWithMoleculeAction.hpp \
   MoleculeAction/LinearInterpolationofTrajectoriesAction.hpp \
   MoleculeAction/RotateAroundSelfByAngleAction.hpp \

src/Actions/MapOfActions.cpp

@@ -71 +71 @@
 #include "Actions/MoleculeAction/CopyAction.hpp"
 #include "Actions/MoleculeAction/FillWithMoleculeAction.hpp"
+#include "Actions/MoleculeAction/FillVoidWithMoleculeAction.hpp"
 #include "Actions/MoleculeAction/LinearInterpolationofTrajectoriesAction.hpp"
 #include "Actions/MoleculeAction/RotateAroundSelfByAngleAction.hpp"
@@ -221 +222 @@
   DescriptionMap["element-db"] = "setting the path where the element databases can be found";
   DescriptionMap["fastparsing"] = "setting whether trajectories shall be parsed completely (n) or just first step (y)";
-  DescriptionMap["fill-molecule"] = "fill empty space of box with a filler molecule";
+  DescriptionMap["fill-molecule"] = "fill around molecules' surface with a filler molecule";
+  DescriptionMap["fill-void"] = "fill void space of box with a filler molecule";
   DescriptionMap["fragment-mol"] = "create for a given molecule into fragments up to given order";
   DescriptionMap["help"] = "Give this help screen";
@@ -361 +363 @@
   TypeMap["fastparsing"] = &typeid(bool);
   TypeMap["fill-molecule"] = &typeid(std::string);
+  TypeMap["fill-void"] = &typeid(std::string);
   TypeMap["fragment-mol"] = &typeid(std::string);
   TypeMap["input"] = &typeid(std::string);
@@ -498 +501 @@
   MenuContainsActionMap.insert( pair<std::string, std::string> ("molecule", "change-molname") );
   MenuContainsActionMap.insert( pair<std::string, std::string> ("molecule", "fill-molecule") );
+  MenuContainsActionMap.insert( pair<std::string, std::string> ("molecule", "fill-void") );
   MenuContainsActionMap.insert( pair<std::string, std::string> ("molecule", "linear-interpolate") );
   MenuContainsActionMap.insert( pair<std::string, std::string> ("molecule", "rotate-origin") );
@@ -571 +575 @@
         generic.insert("fastparsing");
   generic.insert("fill-molecule");
+  generic.insert("fill-void");
   generic.insert("fragment-mol");
   generic.insert("help");
@@ -954 +959 @@
   new MoleculeCopyAction();
   new MoleculeFillWithMoleculeAction();
+  new MoleculeFillVoidWithMoleculeAction();
   new MoleculeLinearInterpolationofTrajectoriesAction();
   new MoleculeRotateAroundSelfByAngleAction();

src/boundary.cpp

@@ -770 +770 @@
 
 
-/** Fills the empty space of the simulation box with water/
+/** Fills the empty space around other molecules' surface of the simulation box with a filler.
  * \param *out output stream for debugging
  * \param *List list of molecules already present in box
@@ -925 +925 @@
 };
 
+/** Fills the empty space of the simulation box with water/
+ * \param *out output stream for debugging
+ * \param *TesselStruct contains tesselated surface
+ * \param *filler molecule which the box is to be filled with
+ * \param configuration contains box dimensions
+ * \param distance[NDIM] distance between filling molecules in each direction
+ * \param boundary length of boundary zone between molecule and filling mollecules
+ * \param epsilon distance to surface which is not filled
+ * \param RandAtomDisplacement maximum distance for random displacement per atom
+ * \param RandMolDisplacement maximum distance for random displacement per filler molecule
+ * \param DoRandomRotation true - do random rotiations, false - don't
+ */
+void FillVoidWithMolecule(molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation)
+{
+  Info FunctionInfo(__func__);
+  molecule *Filling = NULL;
+  Vector CurrentPosition;
+  int N[NDIM];
+  int n[NDIM];
+  const Matrix &M = World::getInstance().getDomain().getM();
+  Matrix Rotations;
+  const Matrix &MInverse = World::getInstance().getDomain().getMinv();
+  Vector AtomTranslations;
+  Vector FillerTranslations;
+  Vector FillerDistance;
+  Vector Inserter;
+  double FillIt = false;
+  bond *Binder = NULL;
+  double phi[NDIM];
+  map<molecule *, LinkedCell *> LCList;
+  std::vector<molecule *> List = World::getInstance().getAllMolecules();
+  MoleculeListClass *MolList = World::getInstance().getMolecules();
+
+  for (std::vector<molecule *>::iterator ListRunner = List.begin(); ListRunner != List.end(); ListRunner++)
+    if ((*ListRunner)->getAtomCount() > 0) {
+      DoLog(1) && (Log() << Verbose(1) << "Pre-creating linked cell lists for molecule " << *ListRunner << "." << endl);
+      LCList[(*ListRunner)] = new LinkedCell(*(*ListRunner), 10.); // get linked cell list
+    }
+
+  // Center filler at origin
+  filler->CenterEdge(&Inserter);
+  const int FillerCount = filler->getAtomCount();
+  DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
+  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
+    for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
+      if ((*BondRunner)->leftatom == *AtomRunner)
+        DoLog(2) && (Log() << Verbose(2) << "  " << *(*BondRunner) << endl);
+
+  // calculate filler grid in [0,1]^3
+  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
+  for(int i=0;i<NDIM;i++)
+    N[i] = (int) ceil(1./FillerDistance[i]);
+  DoLog(1) && (Log() << Verbose(1) << "INFO: Grid steps are " << N[0] << ", " << N[1] << ", " << N[2] << "." << endl);
+
+  // initialize seed of random number generator to current time
+  srand ( time(NULL) );
+
+  // go over [0,1]^3 filler grid
+  for (n[0] = 0; n[0] < N[0]; n[0]++)
+    for (n[1] = 0; n[1] < N[1]; n[1]++)
+      for (n[2] = 0; n[2] < N[2]; n[2]++) {
+        // calculate position of current grid vector in untransformed box
+        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
+        // create molecule random translation vector ...
+        for (int i=0;i<NDIM;i++)
+          FillerTranslations[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+        DoLog(2) && (Log() << Verbose(2) << "INFO: Current Position is " << CurrentPosition << "+" << FillerTranslations << "." << endl);
+
+        // ... and rotation matrix
+        if (DoRandomRotation) {
+          for (int i=0;i<NDIM;i++) {
+            phi[i] = rand()/(RAND_MAX/(2.*M_PI));
+          }
+
+          Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
+          Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
+          Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
+          Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
+          Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
+          Rotations.set(1,2,              sin(phi[1])                                                    );
+          Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
+          Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
+          Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
+        }
+
+        FillIt = true;
+        // go through all atoms
+        for(molecule::const_iterator iter = filler->begin(); iter !=filler->end();++iter){
+
+          // Check whether there is anything too close by
+          LinkedCell::LinkedNodes* liste = NULL;
+          for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); ListRunner != LCList.end(); ++ListRunner) {
+            liste = ListRunner->second->GetPointsInsideSphere(boundary, &CurrentPosition);
+            FillIt = FillIt && (liste->size() == 0);
+            delete(liste);
+            if (!FillIt)
+              break;
+          }
+        }
+        // insert into Filling
+        if (FillIt) {
+          Inserter = CurrentPosition + FillerTranslations;
+          DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is void point." << endl);
+          Filling = filler->CopyMolecule();
+          for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ++miter) {
+            if (DoRandomRotation) {
+              Vector temp = (*miter)->getPosition();
+              temp *= Rotations;
+              (*miter)->setPosition(temp);
+            }
+            // create atomic random translation vector ...
+            for (int i=0;i<NDIM;i++)
+              AtomTranslations[i] = RandomAtomDisplacement*(rand()/(RAND_MAX/2.) - 1.);
+            (*miter)->setPosition((*miter)->getPosition() + AtomTranslations);
+          }
+          Filling->Translate(&Inserter);
+          for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {
+            // check whether each atom is inside box
+            if (!World::getInstance().getDomain().isInside((*miter)->getPosition())) {
+              atom *Walker = *miter;
+              ++miter;
+              World::getInstance().destroyAtom(Walker);
+            } else {
+              ++miter;
+            }
+          }
+          MolList->insert(Filling);
+        } else {
+          DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is non-void point, within boundary or outside of MaxDistance." << endl);
+          continue;
+        }
+      }
+  // last one is replaced by the filler, as we need the original atoms contained therein!
+  filler->Translate(&Inserter);
+  MolList->erase(Filling);
+  for (molecule::iterator iter = Filling->begin(); !Filling->empty(); iter = Filling->begin()) {
+    atom *Walker = *iter;
+    Filling->erase(iter);
+    World::getInstance().destroyAtom(Walker);
+  }
+  World::getInstance().destroyMolecule(Filling);
+
+  for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
+    delete ListRunner->second;
+    LCList.erase(ListRunner);
+  }
+};
 
 /** Tesselates the non convex boundary by rolling a virtual sphere along the surface of the molecule.

src/boundary.hpp

@@ -41 +41 @@
 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename);
 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);
-void FindConvexBorder(const molecule* const mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
+void FillVoidWithMolecule(molecule *filler, config &configuration, const double distance[NDIM], const double boundary, const double RandomAtomDisplacement, const double RandomMolDisplacement, const bool DoRandomRotation);
+    void FindConvexBorder(const molecule* const mol, Boundaries *BoundaryPts, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename);
 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol);
 bool FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LC, const double RADIUS, const char *tempbasename);

tests/regression/testsuite-filling.at

@@ -1 +1 @@
 AT_BANNER([MoleCuilder - Filling in molecules])
 # 1. filling box
-AT_SETUP([Filling - filling empty box])
+AT_SETUP([Filling - filling box around molecule's surface])
 AT_KEYWORDS([filling])
 AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/Filling/1/pre/test.conf .], 0)
@@ -22 +22 @@
 AT_CLEANUP
 
+# 3. filling box
+AT_SETUP([Filling - filling void space in box])
+AT_KEYWORDS([filling])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/Filling/3/pre/test.conf .], 0)
+AT_DATA([water.xyz], [[3
+ # test configuration, created by molecuilder test suite
+O       0.      0.      0.
+H       0.758602 0.     0.504284 
+H       0.758602 0.     -0.504284 
+]])
+AT_CHECK([../../molecuilder -i test.conf -e ${abs_top_srcdir}/src/ -o xyz --fill-void water.xyz --MaxDistance -1 --distances "3.1, 3.1, 3.1"  --lengths "2.1, 0., 0." --DoRotate 0], 0, [stdout], [stderr])
+AT_CHECK([file=test.xyz; diff -I '.*reated by molecuilder.*' $file ${abs_top_srcdir}/${AUTOTEST_PATH}/Filling/3/post/$file], 0, [ignore], [ignore])
+AT_CLEANUP
+
+AT_SETUP([Filling - filling void space besides tenside micelle in box])
+AT_KEYWORDS([filling])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/Filling/3/pre/*.data .], 0)
+AT_CHECK([../../molecuilder -i tensid.data -e ${abs_top_srcdir}/src/ -B "30,0,30,0,0,30" --fill-void water.data --MaxDistance -1 --distances "3.1, 3.1, 3.1"  --lengths "2.1, 0., 0." --DoRotate 0], 0, [stdout], [stderr])
+AT_CHECK([file=tensid.data; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/Filling/3/post/$file], 0, [ignore], [ignore])
+AT_CLEANUP