source: src/Filling/Inserter/RandomInserter.cpp@ 1be100

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 2 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * RandomInserter.cpp
 *
 *  Created on: Feb 21, 2012
 *      Author: heber
 */


// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

//#include "CodePatterns/MemDebug.hpp"

#include "RandomInserter.hpp"

#include <algorithm>

#include "Atom/atom.hpp"
#include "CodePatterns/Log.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "RandomNumbers/RandomNumberGeneratorFactory.hpp"
#include "RandomNumbers/RandomNumberGenerator.hpp"


size_t RandomInserter::Max_Attempts = 10;

/** Sets given 3x3 matrix to a random rotation matrix.
 *
 * @param a matrix to set
 */
inline void setRandomRotation(RealSpaceMatrix &a)
{
  double phi[NDIM];
  RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
  const double rng_min = random.min();
  const double rng_max = random.max();

  for (int i=0;i<NDIM;i++) {
    phi[i] = (random()/(rng_max-rng_min))*(2.*M_PI);
    LOG(4, "DEBUG: Random angle is " << phi[i] << ".");
  }

  a.setRotation(phi);
}

/** Constructor for class RandomInserter.
 *
 * @param _MaxAtomComponent maximum component for random atom translations
 * @param _MaxMoleculeComponent maximum component for random molecule translations
 * @param  _DoRandomRotation whether to do random rotations
 */
RandomInserter::RandomInserter(
    const double _MaxAtomComponent,
    const double _MaxMoleculeComponent,
    const bool _DoRandomRotation) :
   random(RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator()),
   rng_min(random.min()),
   rng_max(random.max()),
   MaxAtomComponent(_MaxAtomComponent),
   MaxMoleculeComponent(_MaxMoleculeComponent),
   DoRandomRotation(_DoRandomRotation)
{}

/** Destructor for class RandomInserter.
 *
 */
RandomInserter::~RandomInserter()
{}

/** Checks whether all atoms currently are inside the cluster's shape.
 *
 * @param cluster cluster to check
 * @return true - all atoms are inside cluster's shape, false - else
 */
bool RandomInserter::AreClustersAtomsInside(ClusterInterface::Cluster_impl cluster) const
{
  ClusterInterface::atomIdSet atoms = cluster->getAtomIds();
  bool status = true;

  for (ClusterInterface::atomIdSet::const_iterator iter = atoms.begin();
      iter != atoms.end();
      ++iter)
    status = status && cluster->isInside(*iter);

  return status;
}

/** Perform the given random translations and rotations on a cluster.
 *
 * @param cluster cluster to translate and rotate
 * @param Rotations random rotation matrix
 * @param RandomAtomTranslations vector with random translation for each atom in cluster
 * @param RandomMoleculeTranslation vector with random translation for cluster
 * @param offset vector with offset for cluster
 */
void RandomInserter::doTranslation(
    ClusterInterface::Cluster_impl cluster,
    const RealSpaceMatrix &Rotations,
    const std::vector<Vector> &RandomAtomTranslations,
    const Vector &RandomMoleculeTranslation) const
{
  AtomIdSet atoms = cluster->getAtoms();

  ASSERT( atoms.size() <=  RandomAtomTranslations.size(),
      "RandomInserter::doTranslation() - insufficient random atom translations given.");

  cluster->transform(Rotations);
  cluster->translate(RandomMoleculeTranslation);
  AtomIdSet::iterator miter = atoms.begin();
  std::vector<Vector>::const_iterator aiter = RandomAtomTranslations.begin();
  for(;miter != atoms.end(); ++miter, ++aiter)
    (*miter)->setPosition((*miter)->getPosition() + *aiter);
  ASSERT( aiter == RandomAtomTranslations.end(),
      "RandomInserter::doTranslation() - there are more RandomAtomTranslations than atoms?");
}

/** Undos a given random translations and rotations on a cluster.
 *
 * @param cluster cluster to translate and rotate
 * @param Rotations random rotation matrix
 * @param RandomAtomTranslations vector with random translation for each atom in cluster
 * @param RandomMoleculeTranslations vector with random translation for cluster
 * @param offset vector with offset for cluster
 */
void RandomInserter::undoTranslation(
    ClusterInterface::Cluster_impl cluster,
    const RealSpaceMatrix &Rotations,
    const std::vector<Vector> &RandomAtomTranslations,
    const Vector &RandomMoleculeTranslation) const
{
  AtomIdSet atoms = cluster->getAtoms();
  ASSERT( atoms.size() <=  RandomAtomTranslations.size(),
      "RandomInserter::doTranslation() - insufficient random atom translations given.");

  // get inverse rotation
  RealSpaceMatrix inverseRotations = Rotations.invert();

  AtomIdSet::iterator miter = atoms.begin();
  std::vector<Vector>::const_iterator aiter = RandomAtomTranslations.begin();
  for(;miter != atoms.end(); ++miter, ++aiter) {
    (*miter)->setPosition((*miter)->getPosition() - *aiter);
  }
  ASSERT( aiter == RandomAtomTranslations.end(),
      "RandomInserter::undoTranslation() - there are more RandomAtomTranslations than atoms?");
  cluster->translate(zeroVec-RandomMoleculeTranslation);
  cluster->transform(inverseRotations);
}

/** Creates a random vector
 *
 * @param range range of components, i.e. \f$ v[i] \in [0,range)\f$
 * @param offset offset for each component
 * @return \a range * rnd() + \a offset
 */
Vector RandomInserter::getRandomVector(const double range, const double offset) const
{
  Vector returnVector;
  for (size_t i=0; i<NDIM; ++i)
    returnVector[i] = range*(random()/((rng_max-rng_min)/2.)) + offset;
  LOG(2, "DEBUG: Generated random vector is " << returnVector << ".");
  return returnVector;
}

/** Inserter operator that randomly translates and rotates the Cluster.
 *
 * \note we assume that clusters are always cloned at origin.
 *
 * @param offset
 * @return true - random translations and rotations did not violate bounding shape, false - else
 */
bool RandomInserter::operator()(ClusterInterface::Cluster_impl cluster, const Vector &offset) const
{
  // calculate center
  AtomIdSet atoms = cluster->getAtoms();
  Vector center;
  center.Zero();
  for(AtomIdSet::iterator miter = atoms.begin();miter != atoms.end(); ++miter)
    center += (*miter)->getPosition();
  center *= 1./atoms.size();

  // shift cluster to center
  cluster->translate(zeroVec-center);

  size_t attempt = 0;
  for (;attempt < Max_Attempts; ++attempt) {
    // generate random rotation matrix
    RealSpaceMatrix Rotations;
    if (DoRandomRotation)
      setRandomRotation(Rotations);
    else
      Rotations.setIdentity();

    // generate random molecule translation vector
    Vector RandomMoleculeTranslation(getRandomVector(MaxMoleculeComponent, -MaxMoleculeComponent));

    // generate random atom translation vector
    std::vector<Vector> RandomAtomTranslations(atoms.size(), zeroVec);
    std::generate_n(RandomAtomTranslations.begin(), atoms.size(),
        boost::bind(&RandomInserter::getRandomVector, this, MaxAtomComponent, -MaxAtomComponent));

    // apply!
//    LOG(2, "DEBUG: cluster before doTranslation(): " << *cluster << ".");
    doTranslation(cluster, Rotations, RandomAtomTranslations, RandomMoleculeTranslation);
//    LOG(2, "DEBUG: cluster after doTranslation(): " << *cluster << ".");

    // ... and check
    if (!AreClustersAtomsInside(cluster)) {
      undoTranslation(cluster, Rotations, RandomAtomTranslations, RandomMoleculeTranslation);
//      LOG(2, "DEBUG: cluster after undoTranslation(): " << *cluster << ".");
    } else {
      break;
    }
  }

  // and move to final position
  cluster->translate(offset+center);

  return attempt != Max_Attempts;
}