Changeset 0c42f2 for src


Ignore:
Timestamp:
Sep 12, 2016, 2:03:11 PM (9 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, 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_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable
Children:
6cb1cd
Parents:
64cafb2
git-author:
Frederik Heber <heber@…> (05/10/16 20:04:47)
git-committer:
Frederik Heber <heber@…> (09/12/16 14:03:11)
Message:

Extended SaturateAction to be able to deal with bonded atoms.

Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/AtomAction/SaturateAction.cpp

    r64cafb2 r0c42f2  
    7575    }
    7676    ASSERT( mol != NULL,
    77         "ChangeElementAction::performCall() - mol is still NULL." );
     77        "AtomSaturateAction::performCall() - mol is still NULL." );
    7878
    79 //    // add the hydrogens
    80 //    const Vector AtomsPosition = _atom->getPosition();
    81 //    double typical_distance = _atom->getType()->getHBondDistance(0);
    82 //    if (typical_distance == -1.)
    83 //      typical_distance = 1.;
    84 //    SphericalPointDistribution PointSphere(typical_distance);
    85 //    PointSphere.initSelf(_atom->getType()->getNoValenceOrbitals());
    86 //    for (SphericalPointDistribution::Polygon_t::const_iterator iter = PointSphere.points.begin();
    87 //        iter != PointSphere.points.end(); ++iter) {
    88 //      // for every Vector add a point relative to atom's position
    89 //      atom * hydrogenAtom = World::getInstance().createAtom();
    90 //      hydrogenAtom->setType(hydrogen);
    91 //      hydrogenAtom->setPosition(AtomsPosition + *iter);
    92 //      mol->AddAtom(hydrogenAtom);
    93 //      // add bond
    94 //      _atom->addBond(WorldTime::getTime(), hydrogenAtom);
    95 //      // mark down for undo state
    96 //      addedHydrogens.push_back(AtomicInfo(*(hydrogenAtom)));
    97 //    }
     79    // radius of polygon on sphere
     80    double typical_distance = _atom->getType()->getHBondDistance(0);
     81    if (typical_distance == -1.)
     82      typical_distance = 1.;
     83
     84    // check for any bonds and get vacant positions
     85    SphericalPointDistribution::Polygon_t vacant_positions;
     86    const BondList& ListOfBonds = _atom->getListOfBonds();
     87    SphericalPointDistribution PointSphere(typical_distance);
     88    if (ListOfBonds.size() == 0) {
     89      vacant_positions = PointSphere.getSimplePolygon(_atom->getType()->getNoValenceOrbitals());
     90      LOG(3, "DEBUG: Using ideal positions as " << vacant_positions);
     91    } else {
     92      // get ideal polygon and currently occupied positions
     93      const SphericalPointDistribution::Polygon_t ideal_positions =
     94          PointSphere.getSimplePolygon(_atom->getType()->getNoValenceOrbitals());
     95      LOG(3, "DEBUG: ideal positions are " << ideal_positions);
     96      SphericalPointDistribution::Polygon_t current_positions;
     97      for (BondList::const_iterator bonditer = ListOfBonds.begin();
     98          bonditer != ListOfBonds.end(); ++bonditer) {
     99        const Vector position =
     100            (*bonditer)->GetOtherAtom(_atom)->getPosition().getVectorToPoint(_atom->getPosition());
     101        current_positions.push_back((1./position.Norm())*position);
     102      }
     103      LOG(3, "DEBUG: current occupied positions are " << current_positions);
     104
     105      // find the best matching rotated polygon
     106      vacant_positions = SphericalPointDistribution::matchSphericalPointDistributions(
     107          current_positions,
     108          ideal_positions);
     109      LOG(3, "DEBUG: Resulting vacant positions are " << vacant_positions);
     110
     111      // scale vacant positions to typical_distance
     112      std::for_each(
     113          vacant_positions.begin(), vacant_positions.end(),
     114          boost::bind(&Vector::Scale, _1, boost::cref(typical_distance)));
     115    }
     116
     117    // add the hydrogens
     118    const Vector AtomsPosition = _atom->getPosition();
     119    for (SphericalPointDistribution::Polygon_t::const_iterator iter = vacant_positions.begin();
     120        iter != vacant_positions.end(); ++iter) {
     121      // for every Vector add a point relative to atom's position
     122      atom * hydrogenAtom = World::getInstance().createAtom();
     123      hydrogenAtom->setType(hydrogen);
     124      hydrogenAtom->setPosition(AtomsPosition + *iter);
     125      mol->AddAtom(hydrogenAtom);
     126      // add bond
     127      _atom->addBond(WorldTime::getTime(), hydrogenAtom);
     128      // mark down for undo state
     129      addedHydrogens.push_back(AtomicInfo(*(hydrogenAtom)));
     130    }
    98131  }
    99132

  • src/Fragmentation/Exporters/SphericalPointDistribution.cpp

    r64cafb2 r0c42f2  
    3939#include "CodePatterns/Assert.hpp"
    4040#include "CodePatterns/IteratorAdaptors.hpp"
     41#include "CodePatterns/Log.hpp"
    4142#include "CodePatterns/toString.hpp"
    4243
    … …  
    5758typedef std::vector<double> DistanceArray_t;
    5859
     60// static instances
     61const double SphericalPointDistribution::SQRT_3(sqrt(3.0));
     62
    5963DistanceArray_t calculatePairwiseDistances(
    60     const std::vector<Vector> &_points,
     64    const std::vector<Vector> &_returnpolygon,
    6165    const IndexList_t &_indices
    6266    )
    … …  
    6973      if (firstiter == seconditer)
    7074        continue;
    71       const double distance = (_points[*firstiter] - _points[*seconditer]).NormSquared();
     75      const double distance = (_returnpolygon[*firstiter] - _returnpolygon[*seconditer]).NormSquared();
    7276      result.push_back(distance);
    7377    }
    … …  
    8993 * \a _new.
    9094 *
    91  * \param _old first set of points (fewer or equal to \a _new)
    92  * \param _new second set of points
     95 * \param _old first set of returnpolygon (fewer or equal to \a _new)
     96 * \param _new second set of returnpolygon
    9397 * \param _Matching matching between the two sets
    9498 * \return pair with L1 and squared L2 error
    … …  
    129133
    130134SphericalPointDistribution::Polygon_t removeMatchingPoints(
    131     const SphericalPointDistribution::Polygon_t &_points,
     135    const SphericalPointDistribution::Polygon_t &_returnpolygon,
    132136    const IndexList_t &_matchingindices
    133137    )
    134138{
    135   SphericalPointDistribution::Polygon_t remainingpoints;
     139  SphericalPointDistribution::Polygon_t remainingreturnpolygon;
    136140  IndexArray_t indices(_matchingindices.begin(), _matchingindices.end());
    137141  std::sort(indices.begin(), indices.end());
    138142  IndexArray_t::const_iterator valueiter = indices.begin();
    139143  SphericalPointDistribution::Polygon_t::const_iterator pointiter =
    140       _points.begin();
    141   for (unsigned int i=0; i< _points.size(); ++i, ++pointiter) {
     144      _returnpolygon.begin();
     145  for (unsigned int i=0; i< _returnpolygon.size(); ++i, ++pointiter) {
    142146    // skip all those in values
    143147    if (*valueiter == i)
    144148      ++valueiter;
    145149    else
    146       remainingpoints.push_back(*pointiter);
    147   }
    148 
    149   return remainingpoints;
     150      remainingreturnpolygon.push_back(*pointiter);
     151  }
     152
     153  return remainingreturnpolygon;
    150154}
    151155
    152156/** Rotates a given polygon around x, y, and z axis by the given angles.
    153157 *
    154  * Essentially, we concentrate on the three points of the polygon to rotate
     158 * Essentially, we concentrate on the three returnpolygon of the polygon to rotate
    155159 * to the correct position. First, we rotate its center via \a angles,
    156160 * then we rotate the "triangle" around itself/\a _RotationAxis by
    157161 * \a _RotationAngle.
    158162 *
    159  * \param _polygon polygon whose points to rotate
     163 * \param _polygon polygon whose returnpolygon to rotate
    160164 * \param _angles vector with rotation angles for x,y,z axis
    161165 * \param _RotationAxis
    … …  
    188192  double bestL2;
    189193  IndexList_t bestmatching;
    190   VectorArray_t oldpoints;
    191   VectorArray_t newpoints;
     194  VectorArray_t oldreturnpolygon;
     195  VectorArray_t newreturnpolygon;
    192196};
    193197
    … …  
    229233      // calculate errors
    230234      std::pair<double, double> errors = calculateErrorOfMatching(
    231           _MCS.oldpoints, _MCS.newpoints, _matching);
     235          _MCS.oldreturnpolygon, _MCS.newreturnpolygon, _matching);
    232236      if (errors.first < L1THRESHOLD) {
    233237        _MCS.bestmatching = _matching;
    … …  
    248252    )
    249253{
    250   SphericalPointDistribution::Polygon_t remainingpoints;
     254  SphericalPointDistribution::Polygon_t remainingreturnpolygon;
    251255  VectorArray_t remainingold(_polygon.begin(), _polygon.end());
    252256  VectorArray_t remainingnew(_newpolygon.begin(), _newpolygon.end());
    … …  
    256260    MCS.foundflag = false;
    257261    MCS.bestL2 = std::numeric_limits<double>::max();
    258     MCS.oldpoints.insert(MCS.oldpoints.begin(), _polygon.begin(),_polygon.end() );
    259     MCS.newpoints.insert(MCS.newpoints.begin(), _newpolygon.begin(),_newpolygon.end() );
     262    MCS.oldreturnpolygon.insert(MCS.oldreturnpolygon.begin(), _polygon.begin(),_polygon.end() );
     263    MCS.newreturnpolygon.insert(MCS.newreturnpolygon.begin(), _newpolygon.begin(),_newpolygon.end() );
    260264
    261265    // search for bestmatching combinatorially
    … …  
    269273      const unsigned int matchingsize = _polygon.size();
    270274      ASSERT( matchingsize <= indices.size(),
    271           "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new points to choose for matching to old ones.");
     275          "SphericalPointDistribution::matchSphericalPointDistributions() - not enough new returnpolygon to choose for matching to old ones.");
    272276      recurseMatchings(MCS, matching, indices, matchingsize);
    273277    }
    … …  
    278282    Vector newCenter;
    279283    {
    280       // calculate center of triangle/line/point consisting of first points of matching
     284      // calculate center of triangle/line/point consisting of first returnpolygon of matching
    281285      Vector oldCenter;
    282286      IndexList_t::const_iterator iter = MCS.bestmatching.begin();
    … …  
    313317        rotatePolygon(_newpolygon, angles, RotationAxis, RotationAngle);
    314318
    315     // remove all points in matching and return remaining ones
     319    // remove all returnpolygon in matching and return remaining ones
    316320    return removeMatchingPoints(rotated_newpolygon, MCS.bestmatching);
    317321  } else
    … …  
    319323}
    320324
    321 
    322 
     325SphericalPointDistribution::Polygon_t
     326SphericalPointDistribution::getSimplePolygon(const int _NumberOfPoints) const
     327{
     328  Polygon_t returnpolygon;
     329
     330        switch (_NumberOfPoints)
     331        {
     332          case 0:
     333            returnpolygon = get<0>();
     334            break;
     335          case 1:
     336            returnpolygon = get<1>();
     337            break;
     338          case 2:
     339            returnpolygon = get<2>();
     340            break;
     341          case 3:
     342            returnpolygon = get<3>();
     343            break;
     344          case 4:
     345            returnpolygon = get<4>();
     346            break;
     347          case 5:
     348            returnpolygon = get<5>();
     349            break;
     350          case 6:
     351            returnpolygon = get<6>();
     352            break;
     353          case 7:
     354            returnpolygon = get<7>();
     355            break;
     356          case 8:
     357            returnpolygon = get<8>();
     358            break;
     359          case 9:
     360            returnpolygon = get<9>();
     361            break;
     362          case 10:
     363            returnpolygon = get<10>();
     364            break;
     365          case 11:
     366            returnpolygon = get<11>();
     367            break;
     368          case 12:
     369            returnpolygon = get<12>();
     370            break;
     371          case 14:
     372            returnpolygon = get<14>();
     373            break;
     374          default:
     375            ASSERT(0, "SphericalPointDistribution::initSelf() - cannot deal with the case "
     376                +toString(_NumberOfPoints)+".");
     377        }
     378
     379        return returnpolygon;
     380}
     381

  • src/Fragmentation/Exporters/SphericalPointDistribution.hpp

    r64cafb2 r0c42f2  
    3939   */
    4040  SphericalPointDistribution(const double _Bondlength = 1.) :
    41     Bondlength(_Bondlength),
    42     SQRT_3(sqrt(3.0))
     41    Bondlength(_Bondlength)
    4342  {}
    4443
    … …  
    5251   * \return Polygon_t with points on the surface centered at (0,0,0)
    5352   */
    54   template <int N> Polygon_t get()
     53  template <int N> Polygon_t get() const
    5554  {
    5655    ASSERT(0, "SphericalPointDistribution::get() - not specialized for "+toString(N)+".");
     56    return Polygon_t();
    5757  }
    5858
     59  /** Initializes the polygon with the given \a _NumberOfPoints.
     60   *
     61   * \param _NumberOfPoints number of points
     62   */
     63  Polygon_t getSimplePolygon(const int _NumberOfPoints) const;
    5964
    6065  /** Matches a given spherical distribution with another containing more
    … …  
    7883  const double Bondlength;
    7984  //!> precalculated value for root of 3
    80   const double SQRT_3;
     85  static const double SQRT_3;
    8186};
    8287
    8388// declare specializations
    8489
    85 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<0>();
    86 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<1>();
    87 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<2>();
    88 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<3>();
    89 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<4>();
    90 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<5>();
    91 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<6>();
    92 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<7>();
    93 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<8>();
    94 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<9>();
    95 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<10>();
    96 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<11>();
    97 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<12>();
    98 template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<14>();
     90template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<0>() const;
     91template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<1>() const;
     92template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<2>() const;
     93template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<3>() const;
     94template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<4>() const;
     95template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<5>() const;
     96template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<6>() const;
     97template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<7>() const;
     98template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<8>() const;
     99template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<9>() const;
     100template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<10>() const;
     101template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<11>() const;
     102template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<12>() const;
     103template <> SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<14>() const;
    99104
    100105#endif /* SPHERICALPOINTDISTRIBUTION_HPP_ */

  • src/Fragmentation/Exporters/SphericalPointDistribution_getPoints.cpp

    r64cafb2 r0c42f2  
    4242
    4343template <>
    44 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<0>()
     44SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<0>() const
    4545{
    4646  Polygon_t polygon;
    … …  
    5151
    5252template <>
    53 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<1>()
     53SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<1>() const
    5454{
    5555  Polygon_t polygon;
    … …  
    6161
    6262template <>
    63 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<2>()
     63SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<2>() const
    6464{
    6565  Polygon_t polygon;
    … …  
    7272
    7373template <>
    74 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<3>()
     74SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<3>() const
    7575{
    7676  Polygon_t polygon;
    … …  
    8484
    8585template <>
    86 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<4>()
     86SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<4>() const
    8787{
    8888  Polygon_t polygon;
    … …  
    9797
    9898template <>
    99 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<5>()
     99SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<5>() const
    100100{
    101101  Polygon_t polygon;
    … …  
    111111
    112112template <>
    113 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<6>()
     113SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<6>() const
    114114{
    115115  Polygon_t polygon;
    … …  
    126126
    127127template <>
    128 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<7>()
     128SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<7>() const
    129129{
    130130  Polygon_t polygon;
    … …  
    142142
    143143template <>
    144 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<8>()
     144SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<8>() const
    145145{
    146146  Polygon_t polygon;
    … …  
    159159
    160160template <>
    161 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<9>()
     161SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<9>() const
    162162{
    163163  Polygon_t polygon;
    … …  
    177177
    178178template <>
    179 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<10>()
     179SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<10>() const
    180180{
    181181  Polygon_t polygon;
    … …  
    197197
    198198template <>
    199 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<11>()
     199SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<11>() const
    200200{
    201201  Polygon_t polygon;
    … …  
    223223
    224224template <>
    225 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<12>()
     225SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<12>() const
    226226{
    227227  Polygon_t polygon;
    … …  
    245245
    246246template <>
    247 SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<14>()
     247SphericalPointDistribution::Polygon_t SphericalPointDistribution::get<14>() const
    248248{
    249249  Polygon_t polygon;
Note: See TracChangeset for help on using the changeset viewer.