Changeset 66fd49 for src


Ignore:
Timestamp:
Feb 3, 2011, 9:51:18 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:
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)
Message:

Rewrite of FillVoidWithMoleculeAction.

FillVoidWithMoleculeAction:

  • new parameter MinDistance and default value of 0.
  • BUGFIX: filler is already created when parsing file, removed useless creation of it initially (also caused lots of confusion due to an "extra" molecule).
  • Undo implemented, regression test inserted.
  • Redo is somewhat hard to implement, as one would use performCall() if it only it would not retrieve its values from ValueStorage ...

FillVoidWithMolecule():

  • filler is now the zeroth not the last molecule, marked by firstInsertion and firstInserter. Filler is removed if no molecules are filled.
  • outsourced stuff into smaller functions
  • removed FillIt to through every atom despite only CurrentPosition, indepedent of atom position, is checked.

TESTFIXES:

  • Analysis/3: test.xyz changed because boundary is now 1.5 instead of 2.1 as 2.1 is not enough of molecules get filled in (and the filler already is).
  • Analysis/3: tensid.data was actually lacking water at (0,0,0) which is after the rewrite present.
Location:
src
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • src/Actions/MoleculeAction/FillVoidWithMoleculeAction.cpp

    r0b15bb r66fd49  
    2828#include "World.hpp"
    2929
     30#include "Descriptors/MoleculeIdDescriptor.hpp"
    3031#include "Parser/MpqcParser.hpp"
    3132#include "Parser/PcpParser.hpp"
     
    3536#include "Parser/FormatParserStorage.hpp"
    3637
     38#include <algorithm>
    3739#include <iostream>
    3840#include <string>
     
    4749/** =========== define the function ====================== */
    4850Action::state_ptr MoleculeFillVoidWithMoleculeAction::performCall() {
    49 
    5051  // obtain information
    5152  getParametersfromValueStorage();
    5253
    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);
    5460  // 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);
    5663  std::string FilenameSuffix = params.fillername.string().substr(params.fillername.string().find_last_of('.')+1, params.fillername.string().length());
    5764  ifstream input;
     
    9299      break;
    93100  }
    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)
    96107    filler = *iter; // get last molecule
    97108  filler->SetNameFromFilename(params.fillername.string().c_str());
     
    102113  for (int i=0;i<NDIM;i++)
    103114    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);
    105116
    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));
    107136}
    108137
    109138Action::state_ptr MoleculeFillVoidWithMoleculeAction::performUndo(Action::state_ptr _state) {
    110 //  MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get());
     139  MoleculeFillVoidWithMoleculeState *state = assert_cast<MoleculeFillVoidWithMoleculeState*>(_state.get());
    111140
    112 //  string newName = state->mol->getName();
    113 //  state->mol->setName(state->lastName);
     141  MoleculeListClass *MolList = World::getInstance().getMolecules();
    114142
    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));
    116160}
    117161
    118162Action::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);
    121167}
    122168
    123169bool MoleculeFillVoidWithMoleculeAction::canUndo() {
    124   return false;
     170  return true;
    125171}
    126172
    127173bool MoleculeFillVoidWithMoleculeAction::shouldUndo() {
    128   return false;
     174  return true;
    129175}
    130176/** =========== end of function ====================== */
  • src/Actions/MoleculeAction/FillVoidWithMoleculeAction.def

    r0b15bb r66fd49  
    1313// ValueStorage by the token "Z" -> first column: int, Z, "Z"
    1414// "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 paramdefaults
    19 #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)
    2020
    21 #undef statetypes
    22 #undef statereferences
     21#define statetypes (std::vector<molecule *>)
     22#define statereferences (fillermolecules)
    2323
    2424// some defines for all the names, you may use ACTION, STATE and PARAMS
  • src/Box.cpp

    r0b15bb r66fd49  
    3333
    3434#include "CodePatterns/Assert.hpp"
     35#include "CodePatterns/Log.hpp"
     36#include "CodePatterns/Verbose.hpp"
    3537
    3638using namespace std;
     
    220222}
    221223
     224double 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
    222239Shape Box::getShape() const{
    223240  return transform(Cuboid(Vector(0,0,0),Vector(1,1,1)),(*M));
    224241}
    225242
    226 const Box::Conditions_t Box::getConditions(){
     243const Box::Conditions_t Box::getConditions() const
     244{
    227245  return conditions;
    228246}
     
    232250}
    233251
    234 const vector<pair<Plane,Plane> >  Box::getBoundingPlanes(){
     252const vector<pair<Plane,Plane> >  Box::getBoundingPlanes() const
     253{
    235254  vector<pair<Plane,Plane> > res;
    236255  for(int i=0;i<NDIM;++i){
     
    248267    res.push_back(make_pair(p1,p2));
    249268  }
     269  ASSERT(res.size() == 3, "Box::getBoundingPlanes() - does not have three plane pairs!");
    250270  return res;
    251271}
  • src/Box.hpp

    r0b15bb r66fd49  
    106106  double periodicDistance(const Vector &point1,const Vector &point2) const;
    107107
     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
    108114  Shape getShape() const;
    109   const Conditions_t getConditions();
     115  const Conditions_t getConditions() const;
    110116  void setCondition(int,BoundaryCondition_t);
    111117
    112   const vector<pair<Plane,Plane> > getBoundingPlanes();
     118  const vector<pair<Plane,Plane> > getBoundingPlanes() const;
    113119
    114120  void setCuboid(const Vector&);
  • src/LinearAlgebra/RealSpaceMatrix.cpp

    r0b15bb r66fd49  
    121121  set(2,2, cos(x)*cos(y));
    122122}
     123
     124void 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
    123143
    124144RealSpaceMatrix &RealSpaceMatrix::operator=(const RealSpaceMatrix &src)
  • src/LinearAlgebra/RealSpaceMatrix.hpp

    r0b15bb r66fd49  
    7272   */
    7373  void setRotation(const double x, const double y, const double z);
     74  void setRandomRotation();
    7475
    7576  /**
  • src/boundary.cpp

    r0b15bb r66fd49  
    927927};
    928928
    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 */
     939void 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 */
     964void 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 */
     992bool 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.
    9321026 * \param *filler molecule which the box is to be filled with
    9331027 * \param configuration contains box dimensions
    9341028 * \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
    9371030 * \param RandAtomDisplacement maximum distance for random displacement per atom
    9381031 * \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 */
     1035void 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    )
    9421045{
    9431046  Info FunctionInfo(__func__);
     
    9491052  RealSpaceMatrix Rotations;
    9501053  const RealSpaceMatrix &MInverse = World::getInstance().getDomain().getMinv();
    951   Vector AtomTranslations;
    9521054  Vector FillerTranslations;
    9531055  Vector FillerDistance;
    9541056  Vector Inserter;
    9551057  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();
    9581061  map<molecule *, LinkedCell *> LCList;
    9591062  std::vector<molecule *> List = World::getInstance().getAllMolecules();
     
    9681071  // Center filler at origin
    9691072  filler->CenterEdge(&Inserter);
    970   const int FillerCount = filler->getAtomCount();
     1073  //const int FillerCount = filler->getAtomCount();
    9711074  DoLog(2) && (Log() << Verbose(2) << "INFO: Filler molecule has the following bonds:" << endl);
    9721075  for(molecule::iterator AtomRunner = filler->begin(); AtomRunner != filler->end(); ++AtomRunner)
     
    9961099
    9971100        // ... 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;
    10121119        }
    10131120
    1014         FillIt = true;
    1015         // go through all atoms
    1016         for(molecule::const_iterator iter = filler->begin(); iter !=filler->end();++iter){
    1017 
    1018           // Check whether there is anything too close by
    1019           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         }
    10281121        // insert into Filling
    10291122        if (FillIt) {
    10301123          Inserter = CurrentPosition + FillerTranslations;
    10311124          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);
    10431139          }
    1044           Filling->Translate(&Inserter);
    1045           for(molecule::iterator miter = Filling->begin(); miter != Filling->end(); ) {
    1046             // check whether each atom is inside box
    1047             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);
    10561140        } else {
    10571141          DoLog(1) && (Log() << Verbose(1) << "INFO: Position at " << Inserter << " is non-void point, within boundary or outside of MaxDistance." << endl);
     
    10591143        }
    10601144      }
    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);
    10701169
    10711170  for (std::map<molecule *, LinkedCell *>::iterator ListRunner = LCList.begin(); !LCList.empty(); ListRunner = LCList.begin()) {
  • src/boundary.hpp

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