source: src/Actions/FillAction/SuspendInMoleculeAction.cpp@ ffa69c

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2014 Frederik Heber. 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/>.
 */

/*
 * SuspendInMoleculeAction.cpp
 *
 *  Created on: Sep 05, 2014
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include "Actions/UndoRedoHelpers.hpp"
#include "Atom/atom.hpp"
#include "Atom/AtomicInfo.hpp"
#include "Atom/CopyAtoms/CopyAtoms_withBonds.hpp"
#include "Bond/BondInfo.hpp"
#include "CodePatterns/Log.hpp"
#include "Descriptors/MoleculeOrderDescriptor.hpp"
#include "Element/element.hpp"
#include "Filling/Cluster.hpp"
#include "Filling/Filler.hpp"
#include "Filling/Preparators/BoxFillerPreparator.hpp"
#include "LinkedCell/linkedcell.hpp"
#include "LinkedCell/PointCloudAdaptor.hpp"
#include "molecule.hpp"
#include "Parser/FormatParserInterface.hpp"
#include "Parser/FormatParserStorage.hpp"
#include "Tesselation/boundary.hpp"
#include "Tesselation/tesselation.hpp"
#include "World.hpp"

#include <algorithm>
#include<gsl/gsl_poly.h>
#include <iostream>
#include <string>
#include <vector>

#include "Actions/FillAction/SuspendInMoleculeAction.hpp"

using namespace MoleCuilder;

static double calculateMass(const molecule &_mol)
{
  // sum up the atomic masses
  const double mass = _mol.getAtomSet().totalMass();
  LOG(2, "DEBUG: Molecule "+_mol.getName()+"'s summed mass is "
      << setprecision(10) << mass << " atomicmassunit.");
  return mass;
}

static double calculateEnvelopeVolume(
    molecule &_mol,
    std::vector<double> &_diameters)
{
  const bool IsAngstroem = true;
  class Tesselation *TesselStruct = NULL;

  Boundaries *BoundaryPoints = GetBoundaryPoints(&_mol, TesselStruct);
  const double * diameters =
      GetDiametersOfCluster(BoundaryPoints, &_mol, TesselStruct, IsAngstroem);
  std::copy(&diameters[0], &diameters[3], _diameters.begin());
  delete[] diameters;
  PointCloudAdaptor< molecule > cloud(&_mol, _mol.getName());
  LinkedCell_deprecated *LCList = new LinkedCell_deprecated(cloud, 10.);
  FindConvexBorder(&_mol, BoundaryPoints, TesselStruct, (const LinkedCell_deprecated *&)LCList, NULL);
  delete (LCList);
  delete[] BoundaryPoints;

  // some preparations beforehand
  const double volume = TesselStruct->getVolumeOfConvexEnvelope(IsAngstroem);

  delete TesselStruct;

  LOG(2, "DEBUG: Molecule "+_mol.getName()+"'s volume is "
      << setprecision(10) << volume << " angstrom^3.");

  return volume;
}

// and construct the stuff
#include "SuspendInMoleculeAction.def"
#include "Action_impl_pre.hpp"
/** =========== define the function ====================== */
ActionState::ptr FillSuspendInMoleculeAction::performCall() {
  typedef std::vector<atom*> AtomVector;

  // get the filler molecule
  std::vector<AtomicInfo> movedatoms;
  molecule *filler = NULL;
  {
    const std::vector< molecule *> molecules = World::getInstance().getSelectedMolecules();
    if (molecules.size() != 1) {
      STATUS("No exactly one molecule selected, aborting,");
      return Action::failure;
    }
    filler = *(molecules.begin());
  }
  for(molecule::const_iterator iter = const_cast<const molecule *>(filler)->begin();
      iter != const_cast<const molecule *>(filler)->end(); ++iter)
    movedatoms.push_back( AtomicInfo(*(*iter)) );
  LOG(1, "INFO: Chosen molecule has " << filler->size() << " atoms.");

  // center filler's tip at origin
  filler->CenterEdge();

  std::vector<molecule *> molecules = World::getInstance().getAllMolecules();
  if (molecules.size() < 2) {
    STATUS("There must be at least two molecules: filler and to be suspended.");
    return Action::failure;
  }

  /// first we need to calculate some volumes and masses
  double totalmass = 0.;
  const bool IsAngstroem = true;
  Vector BoxLengths;
  double clustervolume = 0.;
  std::vector<double> GreatestDiameter(NDIM, 0.);
  for (std::vector<molecule *>::const_iterator iter = molecules.begin();
      iter != molecules.end(); ++iter)
  {
    // skip the filler
    if (*iter == filler)
      continue;
    molecule & mol = **iter;
    const double mass = calculateMass(mol);
    totalmass += mass;
    std::vector<double> diameters(NDIM, 0.);
    const double volume = calculateEnvelopeVolume(mol, diameters);
    clustervolume += volume;
    for (size_t i=0;i<NDIM;++i)
      GreatestDiameter[i] = std::max(GreatestDiameter[i], diameters[i]);
  }
  LOG(1, "INFO: The summed mass is " << setprecision(10)
      << totalmass << " atomicmassunit.");
  LOG(1, "INFO: The average density is " << setprecision(10)
      << totalmass / clustervolume << " atomicmassunit/"
      << (IsAngstroem ? " angstrom" : " atomiclength") << "^3.");
  if ( ((totalmass / clustervolume < 1.) && (params.density.get() > 1.))
      || ((totalmass / clustervolume > 1.) && (params.density.get() < 1.))) {
    STATUS("Desired and present molecular densities must both be either in [0,1) or in (1, inf).");
    return Action::failure;
  }

  // calculate maximum solvent density
  std::vector<double> fillerdiameters(NDIM, 0.);
  const double fillervolume = calculateEnvelopeVolume(*filler, fillerdiameters);
  const double fillermass = calculateMass(*filler);
  LOG(1, "INFO: The filler's mass is " << setprecision(10)
      << fillermass << " atomicmassunit, and it's volume is "
      << fillervolume << (IsAngstroem ? " angstrom" : " atomiclength") << "^3.");
//  const double solventdensity = fillermass / fillervolume;

  /// solve cubic polynomial
  double cellvolume = 0.;
  LOG(1, "Solving equidistant suspension in water problem ...");
  // s = solvent, f = filler, 0 = initial molecules/cluster
  // v_s = v_0 + v_f, m_s = m_0 + rho_f * v_f --> rho_s = m_s/v_s ==> v_f = (m_0 - rho_s * v_o) / (rho_s - rho_f)
  cellvolume = (totalmass - params.density.get() * clustervolume) / (params.density.get() - 1.) + clustervolume;
  LOG(1, "Cellvolume needed for a density of " << params.density.get()
      << " g/cm^3 is " << cellvolume << " angstroem^3.");

  const double minimumvolume =
      (GreatestDiameter[0] * GreatestDiameter[1] * GreatestDiameter[2]);
  LOG(1, "Minimum volume of the convex envelope contained in a rectangular box is "
      << minimumvolume << " angstrom^3.");

  if (minimumvolume > cellvolume) {
    ELOG(1, "The containing box already has a greater volume than the envisaged cell volume!");
    LOG(0, "Setting Box dimensions to minimum possible, the greatest diameters.");
    for (int i = 0; i < NDIM; i++)
      BoxLengths[i] = GreatestDiameter[i];
//    mol->CenterEdge();
  } else {
    BoxLengths[0] = GreatestDiameter[0] + GreatestDiameter[1] + GreatestDiameter[2];
    BoxLengths[1] = GreatestDiameter[0] * GreatestDiameter[1]
        + GreatestDiameter[0] * GreatestDiameter[2]
        + GreatestDiameter[1] * GreatestDiameter[2];
    BoxLengths[2] = minimumvolume - cellvolume;
    std::vector<double> x(3, 0.);
    // for cubic polynomial there are either 1 or 3 unique solutions
    if (gsl_poly_solve_cubic(BoxLengths[0], BoxLengths[1], BoxLengths[2], &x[0], &x[1], &x[2]) == 1) {
      x[1] = x[0];
      x[2] = x[0];
    } else {
      std::swap(x[0], x[2]); // sorted in ascending order
    }
    LOG(0, "RESULT: The resulting spacing is: " << x << " .");

    cellvolume = 1.;
    for (size_t i = 0; i < NDIM; ++i) {
      BoxLengths[i] = x[i] + GreatestDiameter[i];
      cellvolume *= BoxLengths[i];
    }
  }

  // TODO: Determine counts from resulting mass correctly (hard problem due to integers)
  std::vector<unsigned int> counts(3, 0);
  const unsigned int totalcounts = round(params.density.get() * cellvolume - totalmass) / fillermass;
  if (totalcounts > 0) {
    counts[0] = ceil(BoxLengths[0]/3.1);
    counts[1] = ceil(BoxLengths[1]/3.1);
    counts[2] = ceil(BoxLengths[2]/3.1);
  }

  // update Box of atoms by boundary
  {
    RealSpaceMatrix domain;
    for(size_t i =0; i<NDIM;++i)
      domain.at(i,i) = BoxLengths[i];
    World::getInstance().setDomain(domain);
  }
  LOG(0, "RESULT: The resulting cell dimensions are: " << BoxLengths[0] << " and " << BoxLengths[1] << " and " << BoxLengths[2] << " with total volume of " << cellvolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3.");

  // prepare the filler preparator
  BoxFillerPreparator filler_preparator(filler);
  filler_preparator.addVoidPredicate(params.mindistance.get());
  filler_preparator.addRandomInserter(
      params.RandAtomDisplacement.get(),
      params.RandMoleculeDisplacement.get(),
      params.DoRotate.get());
  Vector offset(.5,.5,.5);
  filler_preparator.addCubeMesh(
      counts,
      offset,
      World::getInstance().getDomain().getM());
  if (!filler_preparator()) {
    STATUS("Filler was not fully constructed.");
    return Action::failure;
  }

  // use filler
  bool successflag = false;
  FillSuspendInMoleculeState *UndoState = NULL;
  {
    // fill
    Filler *fillerFunction = filler_preparator.obtainFiller();
    // TODO: When molecule::getBoundingSphere() does not use a sphere anymore,
    // we need to check whether we rotate the molecule randomly. For this to
    // work we need a sphere!
    const Shape s = filler->getBoundingSphere(params.RandAtomDisplacement.get());
    ClusterInterface::Cluster_impl cluster( new Cluster(filler->getAtomIds(), s) );
    CopyAtoms_withBonds copyMethod;
    Filler::ClusterVector_t ClonedClusters;
    successflag = (*fillerFunction)(copyMethod, cluster, ClonedClusters);
    delete fillerFunction;

    // append each cluster's atoms to clonedatoms (however not selected ones)
    std::vector<const atom *> clonedatoms;
    std::vector<AtomicInfo> clonedatominfos;
    for (Filler::ClusterVector_t::const_iterator iter = ClonedClusters.begin();
        iter != ClonedClusters.end(); ++iter) {
      const AtomIdSet &atoms = (*iter)->getAtomIds();
      clonedatoms.reserve(clonedatoms.size()+atoms.size());
      for (AtomIdSet::const_iterator atomiter = atoms.begin(); atomiter != atoms.end(); ++atomiter)
        if (!filler->containsAtom(*atomiter)) {
          clonedatoms.push_back( *atomiter );
          clonedatominfos.push_back( AtomicInfo(*(*atomiter)) );
        }
    }
    std::vector< BondInfo > clonedbonds;
    StoreBondInformationFromAtoms(clonedatoms, clonedbonds);
    LOG(2, "DEBUG: There are " << clonedatominfos.size() << " newly created atoms.");

    if (!successflag) {
      STATUS("Insertion failed, removing inserted clusters, translating original one back");
      RemoveAtomsFromAtomicInfo(clonedatominfos);
      clonedatoms.clear();
      SetAtomsFromAtomicInfo(movedatoms);
    } else {
      std::vector<Vector> MovedToVector(filler->size(), zeroVec);
      std::transform(filler->begin(), filler->end(), MovedToVector.begin(),
          boost::bind(&AtomInfo::getPosition, _1) );
      UndoState = new FillSuspendInMoleculeState(clonedatominfos,clonedbonds,movedatoms,MovedToVector,params);
    }
  }

  if (successflag)
    return ActionState::ptr(UndoState);
  else {
    return Action::failure;
  }
}

ActionState::ptr FillSuspendInMoleculeAction::performUndo(ActionState::ptr _state) {
  FillSuspendInMoleculeState *state = assert_cast<FillSuspendInMoleculeState*>(_state.get());

  // remove all created atoms
  RemoveAtomsFromAtomicInfo(state->clonedatoms);
  // add the original cluster
  SetAtomsFromAtomicInfo(state->movedatoms);

  return ActionState::ptr(_state);
}

ActionState::ptr FillSuspendInMoleculeAction::performRedo(ActionState::ptr _state){
  FillSuspendInMoleculeState *state = assert_cast<FillSuspendInMoleculeState*>(_state.get());

  // place filler cluster again at new spot
  ResetAtomPosition(state->movedatoms, state->MovedToVector);

  // re-create all clusters
  bool statusflag = AddAtomsFromAtomicInfo(state->clonedatoms);

  // re-create the bonds
  if (statusflag)
    AddBondsFromBondInfo(state->clonedbonds);
  if (statusflag)
    return ActionState::ptr(_state);
  else {
    STATUS("Failed re-adding filled in atoms.");
    return Action::failure;
  }
}

bool FillSuspendInMoleculeAction::canUndo() {
  return false;
}

bool FillSuspendInMoleculeAction::shouldUndo() {
  return false;
}
/** =========== end of function ====================== */