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:
20943b
Parents:
632508
git-author:
Frederik Heber <heber@…> (02/24/11 23:49:26)
git-committer:
Frederik Heber <heber@…> (03/01/11 13:17:09)
Message:

Rewrote MinimiseConstrainedPotential into a functor in Dynamics/.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecule_dynamics.cpp

    r632508 r9e23a3  
    2323#include "atom.hpp"
    2424#include "config.hpp"
    25 #include "element.hpp"
     25#include "Dynamics/MinimiseConstrainedPotential.hpp"
     26//#include "element.hpp"
    2627#include "CodePatterns/Info.hpp"
    27 #include "CodePatterns/Verbose.hpp"
    28 #include "CodePatterns/Log.hpp"
     28//#include "CodePatterns/Verbose.hpp"
     29//#include "CodePatterns/Log.hpp"
    2930#include "molecule.hpp"
    3031#include "parser.hpp"
    31 #include "LinearAlgebra/Plane.hpp"
     32//#include "LinearAlgebra/Plane.hpp"
    3233#include "ThermoStatContainer.hpp"
    3334#include "Thermostats/Berendsen.hpp"
    … …  
    3536#include "CodePatterns/enumeration.hpp"
    3637
    37 #include <gsl/gsl_matrix.h>
    38 #include <gsl/gsl_vector.h>
    39 #include <gsl/gsl_linalg.h>
    40 
    4138/************************************* Functions for class molecule *********************************/
    4239
    43 /** Penalizes long trajectories.
    44  * \param *Walker atom to check against others
    45  * \param *mol molecule with other atoms
    46  * \param &Params constraint potential parameters
    47  * \return penalty times each distance
    48  */
    49 double SumDistanceOfTrajectories(atom *Walker, molecule *mol, struct EvaluatePotential &Params)
    50 {
    51   gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
    52   gsl_vector *x = gsl_vector_alloc(NDIM);
    53   atom *Sprinter = NULL;
    54   Vector trajectory1, trajectory2, normal, TestVector;
    55   double Norm1, Norm2, tmp, result = 0.;
    56 
    57   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    58     if ((*iter) == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
    59       break;
    60     // determine normalized trajectories direction vector (n1, n2)
    61     Sprinter = Params.PermutationMap[Walker->getNr()];   // find first target point
    62     trajectory1 = Sprinter->getPositionAtStep(Params.endstep) - Walker->getPositionAtStep(Params.startstep);
    63     trajectory1.Normalize();
    64     Norm1 = trajectory1.Norm();
    65     Sprinter = Params.PermutationMap[(*iter)->getNr()];   // find second target point
    66     trajectory2 = Sprinter->getPositionAtStep(Params.endstep) - (*iter)->getPositionAtStep(Params.startstep);
    67     trajectory2.Normalize();
    68     Norm2 = trajectory1.Norm();
    69     // check whether either is zero()
    70     if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
    71       tmp = Walker->getPositionAtStep(Params.startstep).distance((*iter)->getPositionAtStep(Params.startstep));
    72     } else if (Norm1 < MYEPSILON) {
    73       Sprinter = Params.PermutationMap[Walker->getNr()];   // find first target point
    74       trajectory1 = Sprinter->getPositionAtStep(Params.endstep) - (*iter)->getPositionAtStep(Params.startstep);
    75       trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
    76       trajectory1 -= trajectory2;   // project the part in norm direction away
    77       tmp = trajectory1.Norm();  // remaining norm is distance
    78     } else if (Norm2 < MYEPSILON) {
    79       Sprinter = Params.PermutationMap[(*iter)->getNr()];   // find second target point
    80       trajectory2 = Sprinter->getPositionAtStep(Params.endstep) - Walker->getPositionAtStep(Params.startstep);  // copy second offset
    81       trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
    82       trajectory2 -= trajectory1;   // project the part in norm direction away
    83       tmp = trajectory2.Norm();  // remaining norm is distance
    84     } else if ((fabs(trajectory1.ScalarProduct(trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
    85   //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
    86   //        Log() << Verbose(0) << trajectory1;
    87   //        Log() << Verbose(0) << " and ";
    88   //        Log() << Verbose(0) << trajectory2;
    89       tmp = Walker->getPositionAtStep(Params.startstep).distance((*iter)->getPositionAtStep(Params.startstep));
    90   //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    91     } else { // determine distance by finding minimum distance
    92   //        Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear independent ";
    93   //        Log() << Verbose(0) << endl;
    94   //        Log() << Verbose(0) << "First Trajectory: ";
    95   //        Log() << Verbose(0) << trajectory1 << endl;
    96   //        Log() << Verbose(0) << "Second Trajectory: ";
    97   //        Log() << Verbose(0) << trajectory2 << endl;
    98       // determine normal vector for both
    99       normal = Plane(trajectory1, trajectory2,0).getNormal();
    100       // print all vectors for debugging
    101   //        Log() << Verbose(0) << "Normal vector in between: ";
    102   //        Log() << Verbose(0) << normal << endl;
    103       // setup matrix
    104       for (int i=NDIM;i--;) {
    105         gsl_matrix_set(A, 0, i, trajectory1[i]);
    106         gsl_matrix_set(A, 1, i, trajectory2[i]);
    107         gsl_matrix_set(A, 2, i, normal[i]);
    108         gsl_vector_set(x,i, (Walker->getPositionAtStep(Params.startstep)[i] - (*iter)->getPositionAtStep(Params.startstep)[i]));
    109       }
    110       // solve the linear system by Householder transformations
    111       gsl_linalg_HH_svx(A, x);
    112       // distance from last component
    113       tmp = gsl_vector_get(x,2);
    114   //        Log() << Verbose(0) << " with distance " << tmp << "." << endl;
    115       // test whether we really have the intersection (by checking on c_1 and c_2)
    116       trajectory1.Scale(gsl_vector_get(x,0));
    117       trajectory2.Scale(gsl_vector_get(x,1));
    118       normal.Scale(gsl_vector_get(x,2));
    119       TestVector = (*iter)->getPositionAtStep(Params.startstep) + trajectory2 + normal
    120                    - (Walker->getPositionAtStep(Params.startstep) + trajectory1);
    121       if (TestVector.Norm() < MYEPSILON) {
    122   //          Log() << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
    123       } else {
    124   //          Log() << Verbose(2) << "Test: failed.\tIntersection is off by ";
    125   //          Log() << Verbose(0) << TestVector;
    126   //          Log() << Verbose(0) << "." << endl;
    127       }
    128     }
    129     // add up
    130     tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    131     if (fabs(tmp) > MYEPSILON) {
    132       result += Params.PenaltyConstants[1] * 1./tmp;
    133       //Log() << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
    134     }
    135   }
    136   return result;
    137 };
    138 
    139 /** Penalizes atoms heading to same target.
    140  * \param *Walker atom to check against others
    141  * \param *mol molecule with other atoms
    142  * \param &Params constrained potential parameters
    143  * \return \a penalty times the number of equal targets
    144  */
    145 double PenalizeEqualTargets(atom *Walker, molecule *mol, struct EvaluatePotential &Params)
    146 {
    147   double result = 0.;
    148   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    149     if ((Params.PermutationMap[Walker->getNr()] == Params.PermutationMap[(*iter)->getNr()]) && (Walker->getNr() < (*iter)->getNr())) {
    150   //    atom *Sprinter = PermutationMap[Walker->getNr()];
    151   //        Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
    152   //        Log() << Verbose(0) << Sprinter->getPosition(endstep);
    153   //        Log() << Verbose(0) << ", penalting." << endl;
    154       result += Params.PenaltyConstants[2];
    155       //Log() << Verbose(4) << "Adding " << constants[2] << "." << endl;
    156     }
    157   }
    158   return result;
    159 };
    160 
    161 /** Evaluates the potential energy used for constrained molecular dynamics.
    162  * \f$V_i^{con} = c^{bond} \cdot | r_{P(i)} - R_i | + sum_{i \neq j} C^{min} \cdot \frac{1}{C_{ij}} + C^{inj} \Bigl (1 - \theta \bigl (\prod_{i \neq j} (P(i) - P(j)) \bigr ) \Bigr )\f$
    163  *     where the first term points to the target in minimum distance, the second is a penalty for trajectories lying too close to each other (\f$C_{ij}\f$ is minimum distance between
    164  *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
    165  * Note that for the second term we have to solve the following linear system:
    166  * \f$-c_1 \cdot n_1 + c_2 \cdot n_2 + C \cdot n_3 = - p_2 + p_1\f$, where \f$c_1\f$, \f$c_2\f$ and \f$C\f$ are constants,
    167  * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
    168  * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
    169  * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
    170  * \param *out output stream for debugging
    171  * \param &Params constrained potential parameters
    172  * \return potential energy
    173  * \note This routine is scaling quadratically which is not optimal.
    174  * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
    175  */
    176 double molecule::ConstrainedPotential(struct EvaluatePotential &Params)
    177 {
    178   double tmp = 0.;
    179   double result = 0.;
    180   // go through every atom
    181   atom *Runner = NULL;
    182   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
    183     // first term: distance to target
    184     Runner = Params.PermutationMap[(*iter)->getNr()];   // find target point
    185     tmp = ((*iter)->getPositionAtStep(Params.startstep).distance(Runner->getPositionAtStep(Params.endstep)));
    186     tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
    187     result += Params.PenaltyConstants[0] * tmp;
    188     //Log() << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
    189 
    190     // second term: sum of distances to other trajectories
    191     result += SumDistanceOfTrajectories((*iter), this, Params);
    192 
    193     // third term: penalty for equal targets
    194     result += PenalizeEqualTargets((*iter), this, Params);
    195   }
    196 
    197   return result;
    198 };
    199 
    200 /** print the current permutation map.
    201  * \param *out output stream for debugging
    202  * \param &Params constrained potential parameters
    203  * \param AtomCount number of atoms
    204  */
    205 void PrintPermutationMap(int AtomCount, struct EvaluatePotential &Params)
    206 {
    207   stringstream zeile1, zeile2;
    208   int *DoubleList = new int[AtomCount];
    209   for(int i=0;i<AtomCount;i++)
    210     DoubleList[i] = 0;
    211   int doubles = 0;
    212   zeile1 << "PermutationMap: ";
    213   zeile2 << "                ";
    214   for (int i=0;i<AtomCount;i++) {
    215     Params.DoubleList[Params.PermutationMap[i]->getNr()]++;
    216     zeile1 << i << " ";
    217     zeile2 << Params.PermutationMap[i]->getNr() << " ";
    218   }
    219   for (int i=0;i<AtomCount;i++)
    220     if (Params.DoubleList[i] > 1)
    221     doubles++;
    222   if (doubles >0)
    223     DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl);
    224   delete[](DoubleList);
    225 //  Log() << Verbose(2) << zeile1.str() << endl << zeile2.str() << endl;
    226 };
    227 
    228 /** \f$O(N^2)\f$ operation of calculation distance between each atom pair and putting into DistanceList.
    229  * \param *mol molecule to scan distances in
    230  * \param &Params constrained potential parameters
    231  */
    232 void FillDistanceList(molecule *mol, struct EvaluatePotential &Params)
    233 {
    234   for (int i=mol->getAtomCount(); i--;) {
    235     Params.DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
    236     Params.DistanceList[i]->clear();
    237   }
    238 
    239   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    240     for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
    241       Params.DistanceList[(*iter)->getNr()]->insert( DistancePair((*iter)->getPositionAtStep(Params.startstep).distance((*runner)->getPositionAtStep(Params.endstep)), (*runner)) );
    242     }
    243   }
    244 };
    245 
    246 /** initialize lists.
    247  * \param *out output stream for debugging
    248  * \param *mol molecule to scan distances in
    249  * \param &Params constrained potential parameters
    250  */
    251 void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params)
    252 {
    253   for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    254     Params.StepList[(*iter)->getNr()] = Params.DistanceList[(*iter)->getNr()]->begin();    // stores the step to the next iterator that could be a possible next target
    255     Params.PermutationMap[(*iter)->getNr()] = Params.DistanceList[(*iter)->getNr()]->begin()->second;   // always pick target with the smallest distance
    256     Params.DoubleList[Params.DistanceList[(*iter)->getNr()]->begin()->second->getNr()]++;            // increase this target's source count (>1? not injective)
    257     Params.DistanceIterators[(*iter)->getNr()] = Params.DistanceList[(*iter)->getNr()]->begin();    // and remember which one we picked
    258     DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << Params.DistanceList[(*iter)->getNr()]->begin()->first << "." << endl);
    259   }
    260 };
    261 
    262 /** Try the next nearest neighbour in order to make the permutation map injective.
    263  * \param *out output stream for debugging
    264  * \param *mol molecule
    265  * \param *Walker atom to change its target
    266  * \param &OldPotential old value of constraint potential to see if we do better with new target
    267  * \param &Params constrained potential parameters
    268  */
    269 double TryNextNearestNeighbourForInjectivePermutation(molecule *mol, atom *Walker, double &OldPotential, struct EvaluatePotential &Params)
    270 {
    271   double Potential = 0;
    272   DistanceMap::iterator NewBase = Params.DistanceIterators[Walker->getNr()];  // 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 ((Params.DoubleList[NewBase->second->getNr()] != 0) && (NewBase != Params.DistanceList[Walker->getNr()]->end()));
    276   if (NewBase != Params.DistanceList[Walker->getNr()]->end()) {
    277     Params.PermutationMap[Walker->getNr()] = NewBase->second;
    278     Potential = fabs(mol->ConstrainedPotential(Params));
    279     if (Potential > OldPotential) { // undo
    280       Params.PermutationMap[Walker->getNr()] = Params.DistanceIterators[Walker->getNr()]->second;
    281     } else {  // do
    282       Params.DoubleList[Params.DistanceIterators[Walker->getNr()]->second->getNr()]--;  // decrease the old entry in the doubles list
    283       Params.DoubleList[NewBase->second->getNr()]++;    // increase the old entry in the doubles list
    284       Params.DistanceIterators[Walker->getNr()] = NewBase;
    285       OldPotential = Potential;
    286       DoLog(3) && (Log() << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl);
    287     }
    288   }
    289   return Potential;
    290 };
    291 
    292 /** Permutes \a **&PermutationMap until the penalty is below constants[2].
    293  * \param *out output stream for debugging
    294  * \param *mol molecule to scan distances in
    295  * \param &Params constrained potential parameters
    296  */
    297 void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params)
    298 {
    299   molecule::const_iterator iter = mol->begin();
    300   DistanceMap::iterator NewBase;
    301   double Potential = fabs(mol->ConstrainedPotential(Params));
    302 
    303   if (mol->empty()) {
    304     eLog() << Verbose(1) << "Molecule is empty." << endl;
    305     return;
    306   }
    307   while ((Potential) > Params.PenaltyConstants[2]) {
    308     PrintPermutationMap(mol->getAtomCount(), Params);
    309     iter++;
    310     if (iter == mol->end()) // round-robin at the end
    311       iter = mol->begin();
    312     if (Params.DoubleList[Params.DistanceIterators[(*iter)->getNr()]->second->getNr()] <= 1)  // no need to make those injective that aren't
    313       continue;
    314     // now, try finding a new one
    315     Potential = TryNextNearestNeighbourForInjectivePermutation(mol, (*iter), Potential, Params);
    316   }
    317   for (int i=mol->getAtomCount(); i--;) // now each single entry in the DoubleList should be <=1
    318     if (Params.DoubleList[i] > 1) {
    319       DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
    320       performCriticalExit();
    321     }
    322   DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    323 };
    324 
    325 /** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
    326  * We do the following:
    327  *  -# Generate a distance list from all source to all target points
    328  *  -# Sort this per source point
    329  *  -# Take for each source point the target point with minimum distance, use this as initial permutation
    330  *  -# check whether molecule::ConstrainedPotential() is greater than injective penalty
    331  *     -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
    332  *  -# Next, we only apply transformations that keep the injectivity of the permutations list.
    333  *  -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target
    334  *     point and try to change it for one with lesser distance, or for the next one with greater distance, but only
    335  *     if this decreases the conditional potential.
    336  *  -# finished.
    337  *  -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
    338  *     that the total force is always pointing in direction of the constraint force (ensuring that we move in the
    339  *     right direction).
    340  *  -# Finally, we calculate the potential energy and return.
    341  * \param *out output stream for debugging
    342  * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
    343  * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
    344  * \param endstep step giving final position in constrained MD
    345  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
    346  * \sa molecule::VerletForceIntegration()
    347  * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
    348  * \todo The constrained potential's constants are set to fixed values right now, but they should scale based on checks of the system in order
    349  *       to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
    350  * \bug this all is not O(N log N) but O(N^2)
    351  */
    352 double molecule::MinimiseConstrainedPotential(atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem)
    353 {
    354   double Potential, OldPotential, OlderPotential;
    355   struct EvaluatePotential Params;
    356   Params.PermutationMap = new atom *[getAtomCount()];
    357   Params.DistanceList = new DistanceMap *[getAtomCount()];
    358   Params.DistanceIterators = new DistanceMap::iterator[getAtomCount()];
    359   Params.DoubleList = new int[getAtomCount()];
    360   Params.StepList = new DistanceMap::iterator[getAtomCount()];
    361   int round;
    362   atom *Sprinter = NULL;
    363   DistanceMap::iterator Rider, Strider;
    364 
    365   // set to zero
    366   for (int i=0;i<getAtomCount();i++) {
    367     Params.PermutationMap[i] = NULL;
    368     Params.DoubleList[i] = 0;
    369   }
    370 
    371   /// Minimise the potential
    372   // set Lagrange multiplier constants
    373   Params.PenaltyConstants[0] = 10.;
    374   Params.PenaltyConstants[1] = 1.;
    375   Params.PenaltyConstants[2] = 1e+7;    // just a huge penalty
    376   // generate the distance list
    377   DoLog(1) && (Log() << Verbose(1) << "Allocating, initializting and filling the distance list ... " << endl);
    378   FillDistanceList(this, Params);
    379 
    380   // create the initial PermutationMap (source -> target)
    381   CreateInitialLists(this, Params);
    382 
    383   // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
    384   DoLog(1) && (Log() << Verbose(1) << "Making the PermutationMap injective ... " << endl);
    385   MakeInjectivePermutation(this, Params);
    386   delete[](Params.DoubleList);
    387 
    388   // argument minimise the constrained potential in this injective PermutationMap
    389   DoLog(1) && (Log() << Verbose(1) << "Argument minimising the PermutationMap." << endl);
    390   OldPotential = 1e+10;
    391   round = 0;
    392   do {
    393     DoLog(2) && (Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl);
    394     OlderPotential = OldPotential;
    395     molecule::const_iterator iter;
    396     do {
    397       iter = begin();
    398       for (; iter != end(); ++iter) {
    399         PrintPermutationMap(getAtomCount(), Params);
    400         Sprinter = Params.DistanceIterators[(*iter)->getNr()]->second;   // store initial partner
    401         Strider = Params.DistanceIterators[(*iter)->getNr()];  //remember old iterator
    402         Params.DistanceIterators[(*iter)->getNr()] = Params.StepList[(*iter)->getNr()];
    403         if (Params.DistanceIterators[(*iter)->getNr()] == Params.DistanceList[(*iter)->getNr()]->end()) {// stop, before we run through the list and still on
    404           Params.DistanceIterators[(*iter)->getNr()] == Params.DistanceList[(*iter)->getNr()]->begin();
    405           break;
    406         }
    407         //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->getNr()]->second << "." << endl;
    408         // find source of the new target
    409         molecule::const_iterator runner = begin();
    410         for (; runner != end(); ++runner) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
    411           if (Params.PermutationMap[(*runner)->getNr()] == Params.DistanceIterators[(*iter)->getNr()]->second) {
    412             //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->getNr()] << "." << endl;
    413             break;
    414           }
    415         }
    416         if (runner != end()) { // we found the other source
    417           // then look in its distance list for Sprinter
    418           Rider = Params.DistanceList[(*runner)->getNr()]->begin();
    419           for (; Rider != Params.DistanceList[(*runner)->getNr()]->end(); Rider++)
    420             if (Rider->second == Sprinter)
    421               break;
    422           if (Rider != Params.DistanceList[(*runner)->getNr()]->end()) { // if we have found one
    423             //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->getNr()] << "/" << *Rider->second << "." << endl;
    424             // exchange both
    425             Params.PermutationMap[(*iter)->getNr()] = Params.DistanceIterators[(*iter)->getNr()]->second; // put next farther distance into PermutationMap
    426             Params.PermutationMap[(*runner)->getNr()] = Sprinter;  // and hand the old target to its respective owner
    427             PrintPermutationMap(getAtomCount(), Params);
    428             // calculate the new potential
    429             //Log() << Verbose(2) << "Checking new potential ..." << endl;
    430             Potential = ConstrainedPotential(Params);
    431             if (Potential > OldPotential) { // we made everything worse! Undo ...
    432               //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
    433               //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *Params.DistanceIterators[(*runner)->getNr()]->second << "." << endl;
    434               // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
    435               Params.PermutationMap[(*runner)->getNr()] = Params.DistanceIterators[(*runner)->getNr()]->second;
    436               // Undo for Walker
    437               Params.DistanceIterators[(*iter)->getNr()] = Strider;  // take next farther distance target
    438               //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *Params.DistanceIterators[(*iter)->getNr()]->second << "." << endl;
    439               Params.PermutationMap[(*iter)->getNr()] = Params.DistanceIterators[(*iter)->getNr()]->second;
    440             } else {
    441               Params.DistanceIterators[(*runner)->getNr()] = Rider;  // if successful also move the pointer in the iterator list
    442               DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl);
    443               OldPotential = Potential;
    444             }
    445             if (Potential > Params.PenaltyConstants[2]) {
    446               DoeLog(1) && (eLog()<< Verbose(1) << "The two-step permutation procedure did not maintain injectivity!" << endl);
    447               exit(255);
    448             }
    449             //Log() << Verbose(0) << endl;
    450           } else {
    451             DoeLog(1) && (eLog()<< Verbose(1) << **runner << " was not the owner of " << *Sprinter << "!" << endl);
    452             exit(255);
    453           }
    454         } else {
    455           Params.PermutationMap[(*iter)->getNr()] = Params.DistanceIterators[(*iter)->getNr()]->second; // new target has no source!
    456         }
    457         Params.StepList[(*iter)->getNr()]++; // take next farther distance target
    458       }
    459     } while (++iter != end());
    460   } while ((OlderPotential - OldPotential) > 1e-3);
    461   DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    462 
    463 
    464   /// free memory and return with evaluated potential
    465   for (int i=getAtomCount(); i--;)
    466     Params.DistanceList[i]->clear();
    467   delete[](Params.DistanceList);
    468   delete[](Params.DistanceIterators);
    469   return ConstrainedPotential(Params);
    470 };
    471 
    472 
    473 /** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
    474  * \param *out output stream for debugging
    475  * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
    476  * \param endstep step giving final position in constrained MD
    477  * \param **PermutationMap mapping between the atom label of the initial and the final configuration
    478  * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
    479  * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
    480  */
    481 void molecule::EvaluateConstrainedForces(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
    482 {
    483   double constant = 10.;
    484 
    485   /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
    486   DoLog(1) && (Log() << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl);
    487   BOOST_FOREACH(atom *_atom, atoms) {
    488     atom *Sprinter = PermutationMap[_atom->getNr()];
    489     // set forces
    490     for (int i=NDIM;i++;)
    491       Force->Matrix[0][_atom->getNr()][5+i] += 2.*constant*sqrt(_atom->getPositionAtStep(startstep).distance(Sprinter->getPositionAtStep(endstep)));
    492   }
    493   DoLog(1) && (Log() << Verbose(1) << "done." << endl);
    494 };
    49540
    49641/** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
    … …  
    51459  atom **PermutationMap = NULL;
    51560  atom *Sprinter = NULL;
    516   if (!MapByIdentity)
    517     MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
    518   else {
     61  if (!MapByIdentity) {
     62    MinimiseConstrainedPotential Minimiser(atoms);
     63    Minimiser(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
     64  } else {
    51965    // TODO: construction of enumeration goes here
    52066    PermutationMap = new atom *[getAtomCount()];
    … …  
    610156    // calculate forces and potential
    611157    atom **PermutationMap = NULL;
    612     ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
    613     EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force);
     158    MinimiseConstrainedPotential Minimiser(atoms);
     159    //double ConstrainedPotentialEnergy =
     160    Minimiser(PermutationMap, configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
     161    Minimiser.EvaluateConstrainedForces(&Force);
    614162    delete[](PermutationMap);
    615163  }
Note: See TracChangeset for help on using the changeset viewer.