Changeset ccd9f5


Ignore:
Timestamp:
Oct 9, 2009, 2:35:14 PM (16 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:
4a7776a
Parents:
c111db
Message:

First half of molecule_dynamics.cpp is refactored into smaller parts.

Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/atom.cpp

    rc111db rccd9f5  
    99#include "element.hpp"
    1010#include "memoryallocator.hpp"
     11#include "parser.hpp"
    1112#include "vector.hpp"
    1213
    … …  
    297298};
    298299
     300/** Evaluates some constraint potential if atom moves from \a startstep at once to \endstep in trajectory.
     301 * \param startstep trajectory begins at
     302 * \param endstep trajectory ends at
     303 * \param **PermutationMap if atom switches places with some other atom, there is no translation but a permutaton noted here (not in the trajectories of each).
     304 * \param *Force Force matrix to store result in
     305 */
     306void atom::EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
     307{
     308  double constant = 10.;
     309  atom *Sprinter = PermutationMap[nr];
     310  // set forces
     311  for (int i=NDIM;i++;)
     312    Force->Matrix[0][nr][5+i] += 2.*constant*sqrt(Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
     313};

  • src/atom.hpp

    rc111db rccd9f5  
    2727class bond;
    2828class element;
     29class ForceMatrix;
    2930class Vector;
    3031
    … …  
    8384
    8485  void AddKineticToTemperature(double *temperature, int step) const;
     86  void EvaluateConstrainedForce(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
    8587
    8688  bool IsInParallelepiped(Vector offset, double *parallelepiped);

  • src/molecule.hpp

    rc111db rccd9f5  
    5656/************************************* Class definitions ****************************************/
    5757
    58 
     58/** Structure to contain parameters needed for evaluation of constraint potential.
     59 */
     60struct EvaluatePotential
     61{
     62  int startstep;              //!< start configuration (MDStep in atom::trajectory)
     63  int endstep;                //!< end configuration (MDStep in atom::trajectory)
     64  atom **PermutationMap;      //!< gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
     65  DistanceMap **DistanceList; //!< distance list of each atom to each atom
     66  DistanceMap::iterator *StepList; //!< iterator to ascend through NearestNeighbours \a **DistanceList
     67  int *DoubleList;            //!< count of which sources want to move to this target, basically the injective measure (>1 -> not injective)
     68  DistanceMap::iterator *DistanceIterators; //!< marks which was the last picked target as injective candidate with smallest distance
     69  bool IsAngstroem;           //!< whether coordinates are in angstroem (true) or bohrradius (false)
     70  double *PenaltyConstants;   //!<  penalty constant in front of each term
     71};
    5972
    6073#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
    … …  
    216229
    217230
    218   double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
     231  double ConstrainedPotential(ofstream *out, struct EvaluatePotential &Params);
    219232  double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
    220233  void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);

  • src/molecule_dynamics.cpp

    rc111db rccd9f5  
    1515/************************************* Functions for class molecule *********************************/
    1616
     17/** Penalizes long trajectories.
     18 * \param *Walker atom to check against others
     19 * \param *mol molecule with other atoms
     20 * \param &Params constraint potential parameters
     21 * \return penalty times each distance
     22 */
     23double SumDistanceOfTrajectories(atom *Walker, molecule *mol, struct EvaluatePotential &Params)
     24{
     25  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
     26  gsl_vector *x = gsl_vector_alloc(NDIM);
     27  atom * Runner = mol->start;
     28  atom *Sprinter = NULL;
     29  Vector trajectory1, trajectory2, normal, TestVector;
     30  double Norm1, Norm2, tmp, result = 0.;
     31
     32  while (Runner->next != mol->end) {
     33    Runner = Runner->next;
     34    if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
     35      break;
     36    // determine normalized trajectories direction vector (n1, n2)
     37    Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
     38    trajectory1.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));
     39    trajectory1.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));
     40    trajectory1.Normalize();
     41    Norm1 = trajectory1.Norm();
     42    Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
     43    trajectory2.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));
     44    trajectory2.SubtractVector(&Runner->Trajectory.R.at(Params.startstep));
     45    trajectory2.Normalize();
     46    Norm2 = trajectory1.Norm();
     47    // check whether either is zero()
     48    if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
     49      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.startstep));
     50    } else if (Norm1 < MYEPSILON) {
     51      Sprinter = Params.PermutationMap[Walker->nr];   // find first target point
     52      trajectory1.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));  // copy first offset
     53      trajectory1.SubtractVector(&Runner->Trajectory.R.at(Params.startstep));  // subtract second offset
     54      trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
     55      trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
     56      tmp = trajectory1.Norm();  // remaining norm is distance
     57    } else if (Norm2 < MYEPSILON) {
     58      Sprinter = Params.PermutationMap[Runner->nr];   // find second target point
     59      trajectory2.CopyVector(&Sprinter->Trajectory.R.at(Params.endstep));  // copy second offset
     60      trajectory2.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));  // subtract first offset
     61      trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
     62      trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
     63      tmp = trajectory2.Norm();  // remaining norm is distance
     64    } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
     65  //        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
     66  //        *out << trajectory1;
     67  //        *out << " and ";
     68  //        *out << trajectory2;
     69      tmp = Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.startstep));
     70  //        *out << " with distance " << tmp << "." << endl;
     71    } else { // determine distance by finding minimum distance
     72  //        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
     73  //        *out << endl;
     74  //        *out << "First Trajectory: ";
     75  //        *out << trajectory1 << endl;
     76  //        *out << "Second Trajectory: ";
     77  //        *out << trajectory2 << endl;
     78      // determine normal vector for both
     79      normal.MakeNormalVector(&trajectory1, &trajectory2);
     80      // print all vectors for debugging
     81  //        *out << "Normal vector in between: ";
     82  //        *out << normal << endl;
     83      // setup matrix
     84      for (int i=NDIM;i--;) {
     85        gsl_matrix_set(A, 0, i, trajectory1.x[i]);
     86        gsl_matrix_set(A, 1, i, trajectory2.x[i]);
     87        gsl_matrix_set(A, 2, i, normal.x[i]);
     88        gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep).x[i] - Runner->Trajectory.R.at(Params.startstep).x[i]));
     89      }
     90      // solve the linear system by Householder transformations
     91      gsl_linalg_HH_svx(A, x);
     92      // distance from last component
     93      tmp = gsl_vector_get(x,2);
     94  //        *out << " with distance " << tmp << "." << endl;
     95      // test whether we really have the intersection (by checking on c_1 and c_2)
     96      TestVector.CopyVector(&Runner->Trajectory.R.at(Params.startstep));
     97      trajectory2.Scale(gsl_vector_get(x,1));
     98      TestVector.AddVector(&trajectory2);
     99      normal.Scale(gsl_vector_get(x,2));
     100      TestVector.AddVector(&normal);
     101      TestVector.SubtractVector(&Walker->Trajectory.R.at(Params.startstep));
     102      trajectory1.Scale(gsl_vector_get(x,0));
     103      TestVector.SubtractVector(&trajectory1);
     104      if (TestVector.Norm() < MYEPSILON) {
     105  //          *out << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
     106      } else {
     107  //          *out << Verbose(2) << "Test: failed.\tIntersection is off by ";
     108  //          *out << TestVector;
     109  //          *out << "." << endl;
     110      }
     111    }
     112    // add up
     113    tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     114    if (fabs(tmp) > MYEPSILON) {
     115      result += Params.PenaltyConstants[1] * 1./tmp;
     116      //*out << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
     117    }
     118  }
     119  return result;
     120};
     121
     122/** Penalizes atoms heading to same target.
     123 * \param *Walker atom to check against others
     124 * \param *mol molecule with other atoms
     125 * \param &Params constrained potential parameters
     126 * \return \a penalty times the number of equal targets
     127 */
     128double PenalizeEqualTargets(atom *Walker, molecule *mol, struct EvaluatePotential &Params)
     129{
     130  double result = 0.;
     131  atom * Runner = mol->start;
     132  while (Runner->next != mol->end) {
     133    Runner = Runner->next;
     134    if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
     135  //    atom *Sprinter = PermutationMap[Walker->nr];
     136  //        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
     137  //        *out << Sprinter->Trajectory.R.at(endstep);
     138  //        *out << ", penalting." << endl;
     139      result += Params.PenaltyConstants[2];
     140      //*out << Verbose(4) << "Adding " << constants[2] << "." << endl;
     141    }
     142  }
     143  return result;
     144};
    17145
    18146/** Evaluates the potential energy used for constrained molecular dynamics.
    … …  
    26154 * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
    27155 * \param *out output stream for debugging
    28  * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
    29  * \param startstep start configuration (MDStep in molecule::trajectories)
    30  * \param endstep end configuration (MDStep in molecule::trajectories)
    31  * \param *constants constant in front of each term
    32  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
     156 * \param &Params constrained potential parameters
    33157 * \return potential energy
    34158 * \note This routine is scaling quadratically which is not optimal.
    35159 * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
    36160 */
    37 double molecule::ConstrainedPotential(ofstream *out, atom **PermutationMap, int startstep, int endstep, double *constants, bool IsAngstroem)
    38 {
    39   double result = 0., tmp, Norm1, Norm2;
    40   atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
    41   Vector trajectory1, trajectory2, normal, TestVector;
    42   gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
    43   gsl_vector *x = gsl_vector_alloc(NDIM);
     161double molecule::ConstrainedPotential(ofstream *out, struct EvaluatePotential &Params)
     162{
     163  double tmp, result;
    44164
    45165  // go through every atom
    46   Walker = start;
     166  atom *Runner = NULL;
     167  atom *Walker = start;
    47168  while (Walker->next != end) {
    48169    Walker = Walker->next;
    49170    // first term: distance to target
    50     Runner = PermutationMap[Walker->nr];   // find target point
    51     tmp = (Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(endstep)));
    52     tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    53     result += constants[0] * tmp;
     171    Runner = Params.PermutationMap[Walker->nr];   // find target point
     172    tmp = (Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)));
     173    tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     174    result += Params.PenaltyConstants[0] * tmp;
    54175    //*out << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
    55176
    56177    // second term: sum of distances to other trajectories
    57     Runner = start;
    58     while (Runner->next != end) {
    59       Runner = Runner->next;
    60       if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
    61         break;
    62       // determine normalized trajectories direction vector (n1, n2)
    63       Sprinter = PermutationMap[Walker->nr];   // find first target point
    64       trajectory1.CopyVector(&Sprinter->Trajectory.R.at(endstep));
    65       trajectory1.SubtractVector(&Walker->Trajectory.R.at(startstep));
    66       trajectory1.Normalize();
    67       Norm1 = trajectory1.Norm();
    68       Sprinter = PermutationMap[Runner->nr];   // find second target point
    69       trajectory2.CopyVector(&Sprinter->Trajectory.R.at(endstep));
    70       trajectory2.SubtractVector(&Runner->Trajectory.R.at(startstep));
    71       trajectory2.Normalize();
    72       Norm2 = trajectory1.Norm();
    73       // check whether either is zero()
    74       if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    75         tmp = Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(startstep));
    76       } else if (Norm1 < MYEPSILON) {
    77         Sprinter = PermutationMap[Walker->nr];   // find first target point
    78         trajectory1.CopyVector(&Sprinter->Trajectory.R.at(endstep));  // copy first offset
    79         trajectory1.SubtractVector(&Runner->Trajectory.R.at(startstep));  // subtract second offset
    80         trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    81         trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
    82         tmp = trajectory1.Norm();  // remaining norm is distance
    83       } else if (Norm2 < MYEPSILON) {
    84         Sprinter = PermutationMap[Runner->nr];   // find second target point
    85         trajectory2.CopyVector(&Sprinter->Trajectory.R.at(endstep));  // copy second offset
    86         trajectory2.SubtractVector(&Walker->Trajectory.R.at(startstep));  // subtract first offset
    87         trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    88         trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
    89         tmp = trajectory2.Norm();  // remaining norm is distance
    90       } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
    91 //        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
    92 //        *out << trajectory1;
    93 //        *out << " and ";
    94 //        *out << trajectory2;
    95         tmp = Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(startstep));
    96 //        *out << " with distance " << tmp << "." << endl;
    97       } else { // determine distance by finding minimum distance
    98 //        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
    99 //        *out << endl;
    100 //        *out << "First Trajectory: ";
    101 //        *out << trajectory1 << endl;
    102 //        *out << "Second Trajectory: ";
    103 //        *out << trajectory2 << endl;
    104         // determine normal vector for both
    105         normal.MakeNormalVector(&trajectory1, &trajectory2);
    106         // print all vectors for debugging
    107 //        *out << "Normal vector in between: ";
    108 //        *out << normal << endl;
    109         // setup matrix
    110         for (int i=NDIM;i--;) {
    111           gsl_matrix_set(A, 0, i, trajectory1.x[i]);
    112           gsl_matrix_set(A, 1, i, trajectory2.x[i]);
    113           gsl_matrix_set(A, 2, i, normal.x[i]);
    114           gsl_vector_set(x,i, (Walker->Trajectory.R.at(startstep).x[i] - Runner->Trajectory.R.at(startstep).x[i]));
    115         }
    116         // solve the linear system by Householder transformations
    117         gsl_linalg_HH_svx(A, x);
    118         // distance from last component
    119         tmp = gsl_vector_get(x,2);
    120 //        *out << " with distance " << tmp << "." << endl;
    121         // test whether we really have the intersection (by checking on c_1 and c_2)
    122         TestVector.CopyVector(&Runner->Trajectory.R.at(startstep));
    123         trajectory2.Scale(gsl_vector_get(x,1));
    124         TestVector.AddVector(&trajectory2);
    125         normal.Scale(gsl_vector_get(x,2));
    126         TestVector.AddVector(&normal);
    127         TestVector.SubtractVector(&Walker->Trajectory.R.at(startstep));
    128         trajectory1.Scale(gsl_vector_get(x,0));
    129         TestVector.SubtractVector(&trajectory1);
    130         if (TestVector.Norm() < MYEPSILON) {
    131 //          *out << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
    132         } else {
    133 //          *out << Verbose(2) << "Test: failed.\tIntersection is off by ";
    134 //          *out << TestVector;
    135 //          *out << "." << endl;
    136         }
    137       }
    138       // add up
    139       tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    140       if (fabs(tmp) > MYEPSILON) {
    141         result += constants[1] * 1./tmp;
    142         //*out << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
    143       }
    144     }
     178    result += SumDistanceOfTrajectories(Walker, this, Params);
    145179
    146180    // third term: penalty for equal targets
    147     Runner = start;
    148     while (Runner->next != end) {
    149       Runner = Runner->next;
    150       if ((PermutationMap[Walker->nr] == PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
    151         Sprinter = PermutationMap[Walker->nr];
    152 //        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
    153 //        *out << Sprinter->Trajectory.R.at(endstep);
    154 //        *out << ", penalting." << endl;
    155         result += constants[2];
    156         //*out << Verbose(4) << "Adding " << constants[2] << "." << endl;
    157       }
    158     }
     181    result += PenalizeEqualTargets(Walker, this, Params);
    159182  }
    160183
    … …  
    162185};
    163186
    164 void PrintPermutationMap(ofstream *out, atom **PermutationMap, int Nr)
     187/** print the current permutation map.
     188 * \param *out output stream for debugging
     189 * \param &Params constrained potential parameters
     190 * \param AtomCount number of atoms
     191 */
     192void PrintPermutationMap(ofstream *out, int AtomCount, struct EvaluatePotential &Params)
    165193{
    166194  stringstream zeile1, zeile2;
    167   int *DoubleList = Malloc<int>(Nr, "PrintPermutationMap: *DoubleList");
     195  int *DoubleList = Malloc<int>(AtomCount, "PrintPermutationMap: *DoubleList");
    168196  int doubles = 0;
    169   for (int i=0;i<Nr;i++)
     197  for (int i=0;i<AtomCount;i++)
    170198    DoubleList[i] = 0;
    171199  zeile1 << "PermutationMap: ";
    172200  zeile2 << "                ";
    173   for (int i=0;i<Nr;i++) {
    174     DoubleList[PermutationMap[i]->nr]++;
     201  for (int i=0;i<AtomCount;i++) {
     202    Params.DoubleList[Params.PermutationMap[i]->nr]++;
    175203    zeile1 << i << " ";
    176     zeile2 << PermutationMap[i]->nr << " ";
    177   }
    178   for (int i=0;i<Nr;i++)
    179     if (DoubleList[i] > 1)
     204    zeile2 << Params.PermutationMap[i]->nr << " ";
     205  }
     206  for (int i=0;i<AtomCount;i++)
     207    if (Params.DoubleList[i] > 1)
    180208    doubles++;
    181  // *out << "Found " << doubles << " Doubles." << endl;
     209  if (doubles >0)
     210    *out << "Found " << doubles << " Doubles." << endl;
    182211  Free(&DoubleList);
    183212//  *out << zeile1.str() << endl << zeile2.str() << endl;
     213};
     214
     215/** \f$O(N^2)\f$ operation of calculation distance between each atom pair and putting into DistanceList.
     216 * \param *mol molecule to scan distances in
     217 * \param &Params constrained potential parameters
     218 */
     219void FillDistanceList(molecule *mol, struct EvaluatePotential &Params)
     220{
     221  for (int i=mol->AtomCount; i--;) {
     222    Params.DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
     223    Params.DistanceList[i]->clear();
     224  }
     225
     226  atom *Runner = NULL;
     227  atom *Walker = mol->start;
     228  while (Walker->next != mol->end) {
     229    Walker = Walker->next;
     230    Runner = mol->start;
     231    while(Runner->next != mol->end) {
     232      Runner = Runner->next;
     233      Params.DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(Params.startstep).Distance(&Runner->Trajectory.R.at(Params.endstep)), Runner) );
     234    }
     235  }
     236};
     237
     238/** initialize lists.
     239 * \param *out output stream for debugging
     240 * \param *mol molecule to scan distances in
     241 * \param &Params constrained potential parameters
     242 */
     243void CreateInitialLists(ofstream *out, molecule *mol, struct EvaluatePotential &Params)
     244{
     245  for (int i=mol->AtomCount; i--;)
     246    Params.DoubleList[i] = 0;  // stores for how many atoms in startstep this atom is a target in endstep
     247
     248  atom *Walker = mol->start;
     249  while (Walker->next != mol->end) {
     250    Walker = Walker->next;
     251    Params.StepList[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
     252    Params.PermutationMap[Walker->nr] = Params.DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
     253    Params.DoubleList[Params.DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
     254    Params.DistanceIterators[Walker->nr] = Params.DistanceList[Walker->nr]->begin();    // and remember which one we picked
     255    *out << *Walker << " starts with distance " << Params.DistanceList[Walker->nr]->begin()->first << "." << endl;
     256  }
     257};
     258
     259/** Try the next nearest neighbour in order to make the permutation map injective.
     260 * \param *out output stream for debugging
     261 * \param *mol molecule
     262 * \param *Walker atom to change its target
     263 * \param &OldPotential old value of constraint potential to see if we do better with new target
     264 * \param &Params constrained potential parameters
     265 */
     266double TryNextNearestNeighbourForInjectivePermutation(ofstream *out, molecule *mol, atom *Walker, double &OldPotential, struct EvaluatePotential &Params)
     267{
     268  double Potential = 0;
     269  DistanceMap::iterator NewBase = Params.DistanceIterators[Walker->nr];  // store old base
     270  do {
     271    NewBase++;  // take next further distance in distance to targets list that's a target of no one
     272  } while ((Params.DoubleList[NewBase->second->nr] != 0) && (NewBase != Params.DistanceList[Walker->nr]->end()));
     273  if (NewBase != Params.DistanceList[Walker->nr]->end()) {
     274    Params.PermutationMap[Walker->nr] = NewBase->second;
     275    Potential = fabs(mol->ConstrainedPotential(out, Params));
     276    if (Potential > OldPotential) { // undo
     277      Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second;
     278    } else {  // do
     279      Params.DoubleList[Params.DistanceIterators[Walker->nr]->second->nr]--;  // decrease the old entry in the doubles list
     280      Params.DoubleList[NewBase->second->nr]++;    // increase the old entry in the doubles list
     281      Params.DistanceIterators[Walker->nr] = NewBase;
     282      OldPotential = Potential;
     283      *out << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl;
     284    }
     285  }
     286  return Potential;
     287};
     288
     289/** Permutes \a **&PermutationMap until the penalty is below constants[2].
     290 * \param *out output stream for debugging
     291 * \param *mol molecule to scan distances in
     292 * \param &Params constrained potential parameters
     293 */
     294void MakeInjectivePermutation(ofstream *out, molecule *mol, struct EvaluatePotential &Params)
     295{
     296  atom *Walker = mol->start;
     297  DistanceMap::iterator NewBase;
     298  double Potential = fabs(mol->ConstrainedPotential(out, Params));
     299
     300  while ((Potential) > Params.PenaltyConstants[2]) {
     301    PrintPermutationMap(out, mol->AtomCount, Params);
     302    Walker = Walker->next;
     303    if (Walker == mol->end) // round-robin at the end
     304      Walker = mol->start->next;
     305    if (Params.DoubleList[Params.DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
     306      continue;
     307    // now, try finding a new one
     308    Potential = TryNextNearestNeighbourForInjectivePermutation(out, mol, Walker, Potential, Params);
     309  }
     310  for (int i=mol->AtomCount; i--;) // now each single entry in the DoubleList should be <=1
     311    if (Params.DoubleList[i] > 1) {
     312      cerr << "Failed to create an injective PermutationMap!" << endl;
     313      exit(1);
     314    }
     315  *out << Verbose(1) << "done." << endl;
    184316};
    185317
    … …  
    214346{
    215347  double Potential, OldPotential, OlderPotential;
    216   PermutationMap = Malloc<atom*>(AtomCount, "molecule::MinimiseConstrainedPotential: **PermutationMap");
    217   DistanceMap **DistanceList = Malloc<DistanceMap*>(AtomCount, "molecule::MinimiseConstrainedPotential: **DistanceList");
    218   DistanceMap::iterator *DistanceIterators = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: *DistanceIterators");
    219   int *DoubleList = Malloc<int>(AtomCount, "molecule::MinimiseConstrainedPotential: *DoubleList");
    220   DistanceMap::iterator *StepList = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: *StepList");
    221   double constants[3];
     348  struct EvaluatePotential Params;
     349  Params.PermutationMap = Malloc<atom*>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.**PermutationMap");
     350  Params.DistanceList = Malloc<DistanceMap*>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.**DistanceList");
     351  Params.DistanceIterators = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*DistanceIterators");
     352  Params.DoubleList = Malloc<int>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*DoubleList");
     353  Params.StepList = Malloc<DistanceMap::iterator>(AtomCount, "molecule::MinimiseConstrainedPotential: Params.*StepList");
    222354  int round;
    223355  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
    … …  
    226358  /// Minimise the potential
    227359  // set Lagrange multiplier constants
    228   constants[0] = 10.;
    229   constants[1] = 1.;
    230   constants[2] = 1e+7;    // just a huge penalty
     360  Params.PenaltyConstants[0] = 10.;
     361  Params.PenaltyConstants[1] = 1.;
     362  Params.PenaltyConstants[2] = 1e+7;    // just a huge penalty
    231363  // generate the distance list
    232   *out << Verbose(1) << "Creating the distance list ... " << endl;
    233   for (int i=AtomCount; i--;) {
    234     DoubleList[i] = 0;  // stores for how many atoms in startstep this atom is a target in endstep
    235     DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
    236     DistanceList[i]->clear();
    237   }
    238   *out << Verbose(1) << "Filling the distance list ... " << endl;
    239   Walker = start;
    240   while (Walker->next != end) {
    241     Walker = Walker->next;
    242     Runner = start;
    243     while(Runner->next != end) {
    244       Runner = Runner->next;
    245       DistanceList[Walker->nr]->insert( DistancePair(Walker->Trajectory.R.at(startstep).Distance(&Runner->Trajectory.R.at(endstep)), Runner) );
    246     }
    247   }
     364  *out << Verbose(1) << "Allocating, initializting and filling the distance list ... " << endl;
     365  FillDistanceList(this, Params);
     366
    248367  // create the initial PermutationMap (source -> target)
    249   Walker = start;
    250   while (Walker->next != end) {
    251     Walker = Walker->next;
    252     StepList[Walker->nr] = DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
    253     PermutationMap[Walker->nr] = DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
    254     DoubleList[DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
    255     DistanceIterators[Walker->nr] = DistanceList[Walker->nr]->begin();    // and remember which one we picked
    256     *out << *Walker << " starts with distance " << DistanceList[Walker->nr]->begin()->first << "." << endl;
    257   }
    258   *out << Verbose(1) << "done." << endl;
     368  CreateInitialLists(out, this, Params);
     369
    259370  // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
    260371  *out << Verbose(1) << "Making the PermutationMap injective ... " << endl;
    261   Walker = start;
    262   DistanceMap::iterator NewBase;
    263   OldPotential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
    264   while ((OldPotential) > constants[2]) {
    265     PrintPermutationMap(out, PermutationMap, AtomCount);
    266     Walker = Walker->next;
    267     if (Walker == end) // round-robin at the end
    268       Walker = start->next;
    269     if (DoubleList[DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
    270       continue;
    271     // now, try finding a new one
    272     NewBase = DistanceIterators[Walker->nr];  // store old base
    273     do {
    274       NewBase++;  // take next further distance in distance to targets list that's a target of no one
    275     } while ((DoubleList[NewBase->second->nr] != 0) && (NewBase != DistanceList[Walker->nr]->end()));
    276     if (NewBase != DistanceList[Walker->nr]->end()) {
    277       PermutationMap[Walker->nr] = NewBase->second;
    278       Potential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
    279       if (Potential > OldPotential) { // undo
    280         PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
    281       } else {  // do
    282         DoubleList[DistanceIterators[Walker->nr]->second->nr]--;  // decrease the old entry in the doubles list
    283         DoubleList[NewBase->second->nr]++;    // increase the old entry in the doubles list
    284         DistanceIterators[Walker->nr] = NewBase;
    285         OldPotential = Potential;
    286         *out << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl;
    287       }
    288     }
    289   }
    290   for (int i=AtomCount; i--;) // now each single entry in the DoubleList should be <=1
    291     if (DoubleList[i] > 1) {
    292       cerr << "Failed to create an injective PermutationMap!" << endl;
    293       exit(1);
    294     }
    295   *out << Verbose(1) << "done." << endl;
    296   Free(&DoubleList);
     372  MakeInjectivePermutation(out, this, Params);
     373  Free(&Params.DoubleList);
     374
    297375  // argument minimise the constrained potential in this injective PermutationMap
    298376  *out << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
    … …  
    306384      while (Walker->next != end) { // pick one
    307385        Walker = Walker->next;
    308         PrintPermutationMap(out, PermutationMap, AtomCount);
    309         Sprinter = DistanceIterators[Walker->nr]->second;   // store initial partner
    310         Strider = DistanceIterators[Walker->nr];  //remember old iterator
    311         DistanceIterators[Walker->nr] = StepList[Walker->nr];
    312         if (DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
    313           DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->begin();
     386        PrintPermutationMap(out, AtomCount, Params);
     387        Sprinter = Params.DistanceIterators[Walker->nr]->second;   // store initial partner
     388        Strider = Params.DistanceIterators[Walker->nr];  //remember old iterator
     389        Params.DistanceIterators[Walker->nr] = Params.StepList[Walker->nr];
     390        if (Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
     391          Params.DistanceIterators[Walker->nr] == Params.DistanceList[Walker->nr]->begin();
    314392          break;
    315393        }
    … …  
    318396        Runner = start->next;
    319397        while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
    320           if (PermutationMap[Runner->nr] == DistanceIterators[Walker->nr]->second) {
     398          if (Params.PermutationMap[Runner->nr] == Params.DistanceIterators[Walker->nr]->second) {
    321399            //*out << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
    322400            break;
    … …  
    326404        if (Runner != end) { // we found the other source
    327405          // then look in its distance list for Sprinter
    328           Rider = DistanceList[Runner->nr]->begin();
    329           for (; Rider != DistanceList[Runner->nr]->end(); Rider++)
     406          Rider = Params.DistanceList[Runner->nr]->begin();
     407          for (; Rider != Params.DistanceList[Runner->nr]->end(); Rider++)
    330408            if (Rider->second == Sprinter)
    331409              break;
    332           if (Rider != DistanceList[Runner->nr]->end()) { // if we have found one
     410          if (Rider != Params.DistanceList[Runner->nr]->end()) { // if we have found one
    333411            //*out << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl;
    334412            // exchange both
    335             PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
    336             PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
    337             PrintPermutationMap(out, PermutationMap, AtomCount);
     413            Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
     414            Params.PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
     415            PrintPermutationMap(out, AtomCount, Params);
    338416            // calculate the new potential
    339417            //*out << Verbose(2) << "Checking new potential ..." << endl;
    340             Potential = ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
     418            Potential = ConstrainedPotential(out, Params);
    341419            if (Potential > OldPotential) { // we made everything worse! Undo ...
    342420              //*out << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
    343               //*out << Verbose(3) << "Setting " << *Runner << "'s source to " << *DistanceIterators[Runner->nr]->second << "." << endl;
     421              //*out << Verbose(3) << "Setting " << *Runner << "'s source to " << *Params.DistanceIterators[Runner->nr]->second << "." << endl;
    344422              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
    345               PermutationMap[Runner->nr] = DistanceIterators[Runner->nr]->second;
     423              Params.PermutationMap[Runner->nr] = Params.DistanceIterators[Runner->nr]->second;
    346424              // Undo for Walker
    347               DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
    348               //*out << Verbose(3) << "Setting " << *Walker << "'s source to " << *DistanceIterators[Walker->nr]->second << "." << endl;
    349               PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
     425              Params.DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
     426              //*out << Verbose(3) << "Setting " << *Walker << "'s source to " << *Params.DistanceIterators[Walker->nr]->second << "." << endl;
     427              Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second;
    350428            } else {
    351               DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
     429              Params.DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
    352430              *out << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl;
    353431              OldPotential = Potential;
    354432            }
    355             if (Potential > constants[2]) {
     433            if (Potential > Params.PenaltyConstants[2]) {
    356434              cerr << "ERROR: The two-step permutation procedure did not maintain injectivity!" << endl;
    357435              exit(255);
    … …  
    363441          }
    364442        } else {
    365           PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // new target has no source!
     443          Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second; // new target has no source!
    366444        }
    367         StepList[Walker->nr]++; // take next farther distance target
     445        Params.StepList[Walker->nr]++; // take next farther distance target
    368446      }
    369447    } while (Walker->next != end);
    … …  
    374452  /// free memory and return with evaluated potential
    375453  for (int i=AtomCount; i--;)
    376     DistanceList[i]->clear();
    377   Free(&DistanceList);
    378   Free(&DistanceIterators);
    379   return ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
    380 };
     454    Params.DistanceList[i]->clear();
     455  Free(&Params.DistanceList);
     456  Free(&Params.DistanceIterators);
     457  return ConstrainedPotential(out, Params);
     458};
     459
    381460
    382461/** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
    … …  
    390469void molecule::EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
    391470{
    392   double constant = 10.;
    393   atom *Walker = NULL, *Sprinter = NULL;
    394 
    395471  /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
    396472  *out << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl;
    397   Walker = start;
    398   while (Walker->next != NULL) {
    399     Walker = Walker->next;
    400     Sprinter = PermutationMap[Walker->nr];
    401     // set forces
    402     for (int i=NDIM;i++;)
    403       Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Walker->Trajectory.R.at(startstep).Distance(&Sprinter->Trajectory.R.at(endstep)));
    404   }
     473  ActOnAllAtoms( &atom::EvaluateConstrainedForce, startstep, endstep, PermutationMap, Force );
    405474  *out << Verbose(1) << "done." << endl;
    406475};
Note: See TracChangeset for help on using the changeset viewer.