Changeset e355762 for src


Ignore:
Timestamp:
Mar 1, 2011, 1:17:09 PM (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:
eb33c4
Parents:
fe6f20
git-author:
Frederik Heber <heber@…> (02/25/11 11:41:34)
git-committer:
Frederik Heber <heber@…> (03/01/11 13:17:09)
Message:

Rewrote MinimiseConstrainedPotential to just use std::map instead of arrays.

  • MinimiseConstrainedPotential::operator() needs std::Map<atom*,atom*> as parameter for PermutationMap.
Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/Dynamics/LinearInterpolationBetweenSteps.hpp

    rfe6f20 re355762  
    5959    atom *Sprinter = NULL;
    6060    if (!MapByIdentity) {
     61      std::map<atom *, atom*> PermutationMap;
    6162      molecule::atomSet atoms_list;
    6263      copy(atoms.begin(), atoms.end(), atoms_list.begin());
    63       MinimiseConstrainedPotential Minimiser(atoms_list);
    64       Minimiser(PermutationMap, startstep, endstep, IsAngstroem);
     64      MinimiseConstrainedPotential Minimiser(atoms_list, PermutationMap);
     65      Minimiser(startstep, endstep, IsAngstroem);
    6566    } else {
    6667      // TODO: construction of enumeration goes here

  • src/Dynamics/MinimiseConstrainedPotential.cpp

    rfe6f20 re355762  
    3939
    4040
    41 MinimiseConstrainedPotential::MinimiseConstrainedPotential(molecule::atomSet &_atoms) :
    42   atoms(_atoms)
     41MinimiseConstrainedPotential::MinimiseConstrainedPotential(
     42    molecule::atomSet &_atoms,
     43    std::map<atom*, atom *> &_PermutationMap) :
     44  atoms(_atoms),
     45  PermutationMap(_PermutationMap)
    4346{}
    4447
    … …  
    4649{}
    4750
    48 double MinimiseConstrainedPotential::operator()(atom **&PermutationMap, int _startstep, int _endstep, bool IsAngstroem)
     51double MinimiseConstrainedPotential::operator()(int _startstep, int _endstep, bool IsAngstroem)
    4952{
    5053  double Potential, OldPotential, OlderPotential;
    51   PermutationMap = new atom *[atoms.size()];
    52   DistanceList = new DistanceMap *[atoms.size()];
    53   DistanceIterators = new DistanceMap::iterator[atoms.size()];
    54   DoubleList = new int[atoms.size()];
    55   StepList = new DistanceMap::iterator[atoms.size()];
    5654  int round;
    5755  atom *Sprinter = NULL;
    … …  
    5957
    6058  // set to zero
    61   for (size_t i=0;i<atoms.size();i++) {
    62     PermutationMap[i] = NULL;
    63     DoubleList[i] = 0;
    64   }
     59  PermutationMap.clear();
     60  DoubleList.clear();
     61  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     62    DistanceList[*iter].clear();
     63  }
     64  DistanceList.clear();
     65  DistanceIterators.clear();
     66  DistanceIterators.clear();
    6567
    6668  /// Minimise the potential
    … …  
    7981  DoLog(1) && (Log() << Verbose(1) << "Making the PermutationMap injective ... " << endl);
    8082  MakeInjectivePermutation();
    81   delete[](DoubleList);
     83  DoubleList.clear();
    8284
    8385  // argument minimise the constrained potential in this injective PermutationMap
    … …  
    9294      iter = atoms.begin();
    9395      for (; iter != atoms.end(); ++iter) {
     96        CalculateDoubleList();
    9497        PrintPermutationMap();
    95         Sprinter = DistanceIterators[(*iter)->getNr()]->second;   // store initial partner
    96         Strider = DistanceIterators[(*iter)->getNr()];  //remember old iterator
    97         DistanceIterators[(*iter)->getNr()] = StepList[(*iter)->getNr()];
    98         if (DistanceIterators[(*iter)->getNr()] == DistanceList[(*iter)->getNr()]->end()) {// stop, before we run through the list and still on
    99           DistanceIterators[(*iter)->getNr()] == DistanceList[(*iter)->getNr()]->begin();
     98        Sprinter = DistanceIterators[(*iter)]->second;   // store initial partner
     99        Strider = DistanceIterators[(*iter)];  //remember old iterator
     100        DistanceIterators[(*iter)] = StepList[(*iter)];
     101        if (DistanceIterators[(*iter)] == DistanceList[(*iter)].end()) {// stop, before we run through the list and still on
     102          DistanceIterators[(*iter)] == DistanceList[(*iter)].begin();
    100103          break;
    101104        }
    102         //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->getNr()]->second << "." << endl;
     105        //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)]->second << "." << endl;
    103106        // find source of the new target
    104107        molecule::atomSet::const_iterator runner = atoms.begin();
    105108        for (; runner != atoms.end(); ++runner) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
    106           if (PermutationMap[(*runner)->getNr()] == DistanceIterators[(*iter)->getNr()]->second) {
    107             //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->getNr()] << "." << endl;
     109          if (PermutationMap[(*runner)] == DistanceIterators[(*iter)]->second) {
     110            //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)] << "." << endl;
    108111            break;
    109112          }
    … …  
    111114        if (runner != atoms.end()) { // we found the other source
    112115          // then look in its distance list for Sprinter
    113           Rider = DistanceList[(*runner)->getNr()]->begin();
    114           for (; Rider != DistanceList[(*runner)->getNr()]->end(); Rider++)
     116          Rider = DistanceList[(*runner)].begin();
     117          for (; Rider != DistanceList[(*runner)].end(); Rider++)
    115118            if (Rider->second == Sprinter)
    116119              break;
    117           if (Rider != DistanceList[(*runner)->getNr()]->end()) { // if we have found one
    118             //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->getNr()] << "/" << *Rider->second << "." << endl;
     120          if (Rider != DistanceList[(*runner)].end()) { // if we have found one
     121            //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)] << "/" << *Rider->second << "." << endl;
    119122            // exchange both
    120             PermutationMap[(*iter)->getNr()] = DistanceIterators[(*iter)->getNr()]->second; // put next farther distance into PermutationMap
    121             PermutationMap[(*runner)->getNr()] = Sprinter;  // and hand the old target to its respective owner
     123            PermutationMap[(*iter)] = DistanceIterators[(*iter)]->second; // put next farther distance into PermutationMap
     124            PermutationMap[(*runner)] = Sprinter;  // and hand the old target to its respective owner
     125            CalculateDoubleList();
    122126            PrintPermutationMap();
    123127            // calculate the new potential
    … …  
    126130            if (Potential > OldPotential) { // we made everything worse! Undo ...
    127131              //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
    128               //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *DistanceIterators[(*runner)->getNr()]->second << "." << endl;
     132              //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *DistanceIterators[(*runner)]->second << "." << endl;
    129133              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
    130               PermutationMap[(*runner)->getNr()] = DistanceIterators[(*runner)->getNr()]->second;
     134              PermutationMap[(*runner)] = DistanceIterators[(*runner)]->second;
    131135              // Undo for Walker
    132               DistanceIterators[(*iter)->getNr()] = Strider;  // take next farther distance target
    133               //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *DistanceIterators[(*iter)->getNr()]->second << "." << endl;
    134               PermutationMap[(*iter)->getNr()] = DistanceIterators[(*iter)->getNr()]->second;
     136              DistanceIterators[(*iter)] = Strider;  // take next farther distance target
     137              //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *DistanceIterators[(*iter)]->second << "." << endl;
     138              PermutationMap[(*iter)] = DistanceIterators[(*iter)]->second;
    135139            } else {
    136               DistanceIterators[(*runner)->getNr()] = Rider;  // if successful also move the pointer in the iterator list
     140              DistanceIterators[(*runner)] = Rider;  // if successful also move the pointer in the iterator list
    137141              DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl);
    138142              OldPotential = Potential;
    … …  
    148152          }
    149153        } else {
    150           PermutationMap[(*iter)->getNr()] = DistanceIterators[(*iter)->getNr()]->second; // new target has no source!
     154          PermutationMap[(*iter)] = DistanceIterators[(*iter)]->second; // new target has no source!
    151155        }
    152         StepList[(*iter)->getNr()]++; // take next farther distance target
     156        StepList[(*iter)]++; // take next farther distance target
    153157      }
    154158    } while (++iter != atoms.end());
    … …  
    157161
    158162
    159   /// free memory and return with evaluated potential
    160   for (int i=atoms.size(); i--;)
    161     DistanceList[i]->clear();
    162   delete[](DistanceList);
    163   delete[](DistanceIterators);
    164163  return ConstrainedPotential();
    165164};
    … …  
    167166void MinimiseConstrainedPotential::FillDistanceList()
    168167{
    169   for (int i=atoms.size(); i--;) {
    170     DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
    171     DistanceList[i]->clear();
    172   }
    173 
    174168  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    175169    for (molecule::atomSet::const_iterator runner = atoms.begin(); runner != atoms.end(); ++runner) {
    176       DistanceList[(*iter)->getNr()]->insert( DistancePair((*iter)->getPositionAtStep(startstep).distance((*runner)->getPositionAtStep(endstep)), (*runner)) );
     170      DistanceList[(*iter)].insert( DistancePair((*iter)->getPositionAtStep(startstep).distance((*runner)->getPositionAtStep(endstep)), (*runner)) );
    177171    }
    178172  }
    … …  
    182176{
    183177  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    184     StepList[(*iter)->getNr()] = DistanceList[(*iter)->getNr()]->begin();    // stores the step to the next iterator that could be a possible next target
    185     PermutationMap[(*iter)->getNr()] = DistanceList[(*iter)->getNr()]->begin()->second;   // always pick target with the smallest distance
    186     DoubleList[DistanceList[(*iter)->getNr()]->begin()->second->getNr()]++;            // increase this target's source count (>1? not injective)
    187     DistanceIterators[(*iter)->getNr()] = DistanceList[(*iter)->getNr()]->begin();    // and remember which one we picked
    188     DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << DistanceList[(*iter)->getNr()]->begin()->first << "." << endl);
     178    StepList[(*iter)] = DistanceList[(*iter)].begin();    // stores the step to the next iterator that could be a possible next target
     179    PermutationMap[(*iter)] = DistanceList[(*iter)].begin()->second;   // always pick target with the smallest distance
     180    DoubleList[DistanceList[(*iter)].begin()->second]++;            // increase this target's source count (>1? not injective)
     181    DistanceIterators[(*iter)] = DistanceList[(*iter)].begin();    // and remember which one we picked
     182    DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << DistanceList[(*iter)].begin()->first << "." << endl);
    189183  }
    190184};
    … …  
    201195  }
    202196  while ((Potential) > PenaltyConstants[2]) {
     197    CalculateDoubleList();
    203198    PrintPermutationMap();
    204199    iter++;
    205200    if (iter == atoms.end()) // round-robin at the end
    206201      iter = atoms.begin();
    207     if (DoubleList[DistanceIterators[(*iter)->getNr()]->second->getNr()] <= 1)  // no need to make those injective that aren't
     202    if (DoubleList[DistanceIterators[(*iter)]->second] <= 1)  // no need to make those injective that aren't
    208203      continue;
    209204    // now, try finding a new one
    210205    Potential = TryNextNearestNeighbourForInjectivePermutation((*iter), Potential);
    211206  }
    212   for (int i=atoms.size(); i--;) // now each single entry in the DoubleList should be <=1
    213     if (DoubleList[i] > 1) {
     207  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     208    // now each single entry in the DoubleList should be <=1
     209    if (DoubleList[*iter] > 1) {
    214210      DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
    215211      performCriticalExit();
    216212    }
     213  }
    217214  DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    218215};
    219216
     217unsigned int MinimiseConstrainedPotential::CalculateDoubleList()
     218{
     219  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter)
     220    DoubleList[*iter] = 0;
     221  unsigned int doubles = 0;
     222  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter)
     223      DoubleList[ PermutationMap[*iter] ]++;
     224  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter)
     225    if (DoubleList[*iter] > 1)
     226      doubles++;
     227  if (doubles >0)
     228    DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl);
     229  return doubles;
     230};
     231
    220232void MinimiseConstrainedPotential::PrintPermutationMap() const
    221233{
    222234  stringstream zeile1, zeile2;
    223   const int AtomCount = atoms.size();
    224   int *DoubleList = new int[AtomCount];
    225   for(int i=0;i<AtomCount;i++)
    226     DoubleList[i] = 0;
    227235  int doubles = 0;
    228236  zeile1 << "PermutationMap: ";
    229237  zeile2 << "                ";
    230   for (int i=0;i<AtomCount;i++) {
    231     DoubleList[PermutationMap[i]->getNr()]++;
    232     zeile1 << i << " ";
    233     zeile2 << PermutationMap[i]->getNr() << " ";
    234   }
    235   for (int i=0;i<AtomCount;i++)
    236     if (DoubleList[i] > 1)
    237     doubles++;
     238  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     239    zeile1 << (*iter)->getName() << " ";
     240    zeile2 << (PermutationMap[*iter])->getName() << " ";
     241  }
     242  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     243    std::map<atom *, unsigned int>::const_iterator value_iter = DoubleList.find(*iter);
     244    if (value_iter->second > (unsigned int)1)
     245      doubles++;
     246  }
    238247  if (doubles >0)
    239248    DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl);
    240   delete[](DoubleList);
    241249//  Log() << Verbose(2) << zeile1.str() << endl << zeile2.str() << endl;
    242250};
    … …  
    250258  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    251259    // first term: distance to target
    252     Runner = PermutationMap[(*iter)->getNr()];   // find target point
     260    Runner = PermutationMap[(*iter)];   // find target point
    253261    tmp = ((*iter)->getPositionAtStep(startstep).distance(Runner->getPositionAtStep(endstep)));
    254262    tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    … …  
    269277{
    270278  double Potential = 0;
    271   DistanceMap::iterator NewBase = DistanceIterators[Walker->getNr()];  // store old base
     279  DistanceMap::iterator NewBase = DistanceIterators[Walker];  // store old base
    272280  do {
    273281    NewBase++;  // take next further distance in distance to targets list that's a target of no one
    274   } while ((DoubleList[NewBase->second->getNr()] != 0) && (NewBase != DistanceList[Walker->getNr()]->end()));
    275   if (NewBase != DistanceList[Walker->getNr()]->end()) {
    276     PermutationMap[Walker->getNr()] = NewBase->second;
     282  } while ((DoubleList[NewBase->second] != 0) && (NewBase != DistanceList[Walker].end()));
     283  if (NewBase != DistanceList[Walker].end()) {
     284    PermutationMap[Walker] = NewBase->second;
    277285    Potential = fabs(ConstrainedPotential());
    278286    if (Potential > OldPotential) { // undo
    279       PermutationMap[Walker->getNr()] = DistanceIterators[Walker->getNr()]->second;
     287      PermutationMap[Walker] = DistanceIterators[Walker]->second;
    280288    } else {  // do
    281       DoubleList[DistanceIterators[Walker->getNr()]->second->getNr()]--;  // decrease the old entry in the doubles list
    282       DoubleList[NewBase->second->getNr()]++;    // increase the old entry in the doubles list
    283       DistanceIterators[Walker->getNr()] = NewBase;
     289      DoubleList[DistanceIterators[Walker]->second]--;  // decrease the old entry in the doubles list
     290      DoubleList[NewBase->second]++;    // increase the old entry in the doubles list
     291      DistanceIterators[Walker] = NewBase;
    284292      OldPotential = Potential;
    285293      DoLog(3) && (Log() << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl);
    … …  
    293301  double result = 0.;
    294302  for (molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    295     if ((PermutationMap[Walker->getNr()] == PermutationMap[(*iter)->getNr()]) && (Walker->getNr() < (*iter)->getNr())) {
     303    if ((PermutationMap[Walker] == PermutationMap[(*iter)]) && (Walker < (*iter))) {
    296304  //    atom *Sprinter = PermutationMap[Walker->nr];
    297305  //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
    … …  
    317325      break;
    318326    // determine normalized trajectories direction vector (n1, n2)
    319     Sprinter = PermutationMap[Walker->getNr()];   // find first target point
     327    Sprinter = PermutationMap[Walker];   // find first target point
    320328    trajectory1 = Sprinter->getPositionAtStep(endstep) - Walker->getPositionAtStep(startstep);
    321329    trajectory1.Normalize();
    322330    Norm1 = trajectory1.Norm();
    323     Sprinter = PermutationMap[(*iter)->getNr()];   // find second target point
     331    Sprinter = PermutationMap[(*iter)];   // find second target point
    324332    trajectory2 = Sprinter->getPositionAtStep(endstep) - (*iter)->getPositionAtStep(startstep);
    325333    trajectory2.Normalize();
    … …  
    329337      tmp = Walker->getPositionAtStep(startstep).distance((*iter)->getPositionAtStep(startstep));
    330338    } else if (Norm1 < MYEPSILON) {
    331       Sprinter = PermutationMap[Walker->getNr()];   // find first target point
     339      Sprinter = PermutationMap[Walker];   // find first target point
    332340      trajectory1 = Sprinter->getPositionAtStep(endstep) - (*iter)->getPositionAtStep(startstep);
    333341      trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    … …  
    335343      tmp = trajectory1.Norm();  // remaining norm is distance
    336344    } else if (Norm2 < MYEPSILON) {
    337       Sprinter = PermutationMap[(*iter)->getNr()];   // find second target point
     345      Sprinter = PermutationMap[(*iter)];   // find second target point
    338346      trajectory2 = Sprinter->getPositionAtStep(endstep) - Walker->getPositionAtStep(startstep);  // copy second offset
    339347      trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    … …  
    402410  DoLog(1) && (Log() << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl);
    403411  for(molecule::atomSet::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    404     atom *Sprinter = PermutationMap[(*iter)->getNr()];
     412    atom *Sprinter = PermutationMap[(*iter)];
    405413    // set forces
    406414    for (int i=NDIM;i++;)

  • src/Dynamics/MinimiseConstrainedPotential.hpp

    rfe6f20 re355762  
    2727{
    2828public:
    29   MinimiseConstrainedPotential(molecule::atomSet &_atoms);
     29  /** Constructor.
     30   *
     31   * @param _atoms set of atoms to operate on
     32   * \param _PermutationMap on return: mapping between the atom label of the initial and the final configuration
     33   * @return
     34   */
     35  MinimiseConstrainedPotential(molecule::atomSet &_atoms, std::map<atom*, atom *> &_PermutationMap);
     36
     37  /** Destructor.
     38   *
     39   * @return
     40   */
    3041  ~MinimiseConstrainedPotential();
    3142
    … …  
    4657   *     right direction).
    4758   *  -# Finally, we calculate the potential energy and return.
    48    * \param *out output stream for debugging
    49    * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
    5059   * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
    5160   * \param endstep step giving final position in constrained MD
    … …  
    5766   * \bug this all is not O(N log N) but O(N^2)
    5867   */
    59   double operator()(atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem);
     68  double operator()(int startstep, int endstep, bool IsAngstroem);
    6069
    6170  /** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
    62    * \param *out output stream for debugging
    63    * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
    64    * \param endstep step giving final position in constrained MD
    65    * \param **PermutationMap mapping between the atom label of the initial and the final configuration
    6671   * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
    6772   * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
    … …  
    7782  int startstep; //!< start configuration (MDStep in atom::trajectory)
    7883  int endstep; //!< end configuration (MDStep in atom::trajectory)
    79   atom **PermutationMap; //!< gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
    80   DistanceMap **DistanceList; //!< distance list of each atom to each atom
    81   DistanceMap::iterator *StepList; //!< iterator to ascend through NearestNeighbours \a **DistanceList
    82   int *DoubleList; //!< count of which sources want to move to this target, basically the injective measure (>1 -> not injective)
    83   DistanceMap::iterator *DistanceIterators; //!< marks which was the last picked target as injective candidate with smallest distance
     84  std::map<atom*, atom *> &PermutationMap; //!< gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
     85  std::map<atom *, DistanceMap> DistanceList; //!< distance list of each atom to each atom
     86  std::map<atom *, DistanceMap::iterator> StepList; //!< iterator to ascend through NearestNeighbours \a **DistanceList
     87  std::map<atom *, unsigned int> DoubleList; //!< count of which sources want to move to this target, basically the injective measure (>1 -> not injective)
     88  std::map<atom *, DistanceMap::iterator> DistanceIterators; //!< marks which was the last picked target as injective candidate with smallest distance
    8489  bool IsAngstroem; //!< whether coordinates are in angstroem (true) or bohrradius (false)
    8590  double *PenaltyConstants; //!<  penalty constant in front of each term
    8691
    8792  /** \f$O(N^2)\f$ operation of calculation distance between each atom pair and putting into DistanceList.
    88    * \param *mol molecule to scan distances in
    89    * \param &Params constrained potential parameters
    9093   */
    9194  void FillDistanceList();
    9295
    9396  /** Initialize lists.
    94    * \param *out output stream for debugging
    95    * \param *mol molecule to scan distances in
    96    * \param &Params constrained potential parameters
    9797   */
    9898  void CreateInitialLists();
    9999
    100100  /** Permutes \a **&PermutationMap until the penalty is below constants[2].
    101    * \param *out output stream for debugging
    102    * \param *mol molecule to scan distances in
    103    * \param &Params constrained potential parameters
    104101   */
    105102  void MakeInjectivePermutation();
    106103
     104  /** Calculates the number of doubles in PermutationMap.
     105   */
     106  unsigned int CalculateDoubleList();
     107
    107108  /** Print the current permutation map.
    108    * \param *out output stream for debugging
    109    * \param &Params constrained potential parameters
    110    * \param AtomCount number of atoms
    111109   */
    112110  void PrintPermutationMap() const;
    … …  
    121119   * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
    122120   * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
    123    * \param *out output stream for debugging
    124    * \param &Params constrained potential parameters
    125121   * \return potential energy
    126122   * \note This routine is scaling quadratically which is not optimal.
    … …  
    130126
    131127  /** Try the next nearest neighbour in order to make the permutation map injective.
    132    * \param *out output stream for debugging
    133    * \param *mol molecule
    134128   * \param *Walker atom to change its target
    135129   * \param &OldPotential old value of constraint potential to see if we do better with new target
    136    * \param &Params constrained potential parameters
    137130   */
    138131  double TryNextNearestNeighbourForInjectivePermutation(atom *Walker, double &OldPotential);
    … …  
    140133  /** Penalizes atoms heading to same target.
    141134   * \param *Walker atom to check against others
    142    * \param *mol molecule with other atoms
    143    * \param &Params constrained potential parameters
    144135   * \return \a penalty times the number of equal targets
    145136   */
    … …  
    148139  /** Penalizes long trajectories.
    149140   * \param *Walker atom to check against others
    150    * \param *mol molecule with other atoms
    151    * \param &Params constraint potential parameters
    152141   * \return penalty times each distance
    153142   */

  • src/molecule_dynamics.cpp

    rfe6f20 re355762  
    8585  if (configuration.DoConstrainedMD) {
    8686    // calculate forces and potential
    87     atom **PermutationMap = NULL;
    88     MinimiseConstrainedPotential Minimiser(atoms);
     87    std::map<atom*, atom *> PermutationMap;
     88    MinimiseConstrainedPotential Minimiser(atoms, PermutationMap);
    8989    //double ConstrainedPotentialEnergy =
    90     Minimiser(PermutationMap, configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
     90    Minimiser(configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
    9191    Minimiser.EvaluateConstrainedForces(&Force);
    92     delete[](PermutationMap);
    9392  }
    9493
Note: See TracChangeset for help on using the changeset viewer.