Changeset 85bac0

Timestamp:

Sep 11, 2008, 1:28:25 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:

c30cda

Parents:

a1fe77

git-author:

Frederik Heber <heber@…> (09/09/08 15:29:58)

git-committer:

Frederik Heber <heber@…> (09/11/08 13:28:25)

Message:

Implemented molecule::LinearInterpolationBetweenConfiguration().

command line option "-L" with start and end step performs a linear interpolation between two atomic configurations. So far the mapping from initial atom labels to final labels is not yet finished, it is injective, but not yet minimal.

Location:

src

Files:

• builder.cpp (modified) (2 diffs)
• moleculelist.cpp (modified) (5 diffs)
• molecules.cpp (modified) (15 diffs)
• molecules.hpp (modified) (4 diffs)

src/builder.cpp

@@ -820 +820 @@
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
             cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
+            cout << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl; 
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
@@ -989 +990 @@
               argptr+=1;
               break;
+            case 'L':
+              ExitFlag = 1;
+              SaveFlag = true;
+              cout << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
+              if (!mol->LinearInterpolationBetweenConfiguration((ofstream *)&cout, atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration))
+                cout << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
+              else
+                cout << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
+              argptr+=3;
+              break;
             case 'P':
               ExitFlag = 1;

src/moleculelist.cpp

@@ -378 +378 @@
  * \param *configuration standard configuration to attach atoms in fragment molecule to.
  * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
+ * \param DoPeriodic true - call ScanForPeriodicCorrection, false - don't
+ * \param DoCentering true - call molecule::CenterEdge(), false - don't
  * \return true - success (each file was written), false - something went wrong.
  */
-bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
+bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, const char *fragmentprefix, config *configuration, int *SortIndex, bool DoPeriodic, bool DoCentering)
 {
   ofstream outputFragment;
@@ -404 +406 @@
 
     // correct periodic
-    ListOfMolecules[i]->ScanForPeriodicCorrection(out);
+    if (DoPeriodic)
+      ListOfMolecules[i]->ScanForPeriodicCorrection(out);
 
     // output xyz file
     FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
-    sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, fragmentprefix, FragmentNumber);
     outputFragment.open(FragmentName, ios::out);
-    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
+    *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
     if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
       *out << " done." << endl;
@@ -429 +432 @@
     
     // center on edge
-    ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
-    ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
-    int j = -1;
-    for (int k=0;k<NDIM;k++) {
-      j += k+1;
-      BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
-      ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
-    }
-    ListOfMolecules[i]->Translate(&BoxDimension);
+    if (DoCentering) {
+      ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
+      ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
+      int j = -1;
+      for (int k=0;k<NDIM;k++) {
+        j += k+1;
+        BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
+        ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
+      }
+      ListOfMolecules[i]->Translate(&BoxDimension);
+    }
 
     // also calculate necessary orbitals
@@ -445 +450 @@
     // change path in config
     //strcpy(PathBackup, configuration->configpath);
-    sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
+    sprintf(FragmentName, "%s/%s%s/", PathBackup, fragmentprefix, FragmentNumber);
     configuration->SetDefaultPath(FragmentName);
     
     // and save as config
-    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, fragmentprefix, FragmentNumber);
     outputFragment.open(FragmentName, ios::out);
-    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
+    *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
     if ((intermediateResult = configuration->Save(&outputFragment, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
       *out << " done." << endl;
@@ -465 +470 @@
 
     // and save as mpqc input file
-    sprintf(FragmentName, "%s/%s%s.in", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    sprintf(FragmentName, "%s/%s%s.in", configuration->configpath, fragmentprefix, FragmentNumber);
     outputFragment.open(FragmentName, ios::out);
-    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
+    *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
     if ((intermediateResult = configuration->SaveMPQC(&outputFragment, ListOfMolecules[i])))
       *out << " done." << endl;

src/molecules.cpp

@@ -1008 +1008 @@
 
 /** Evaluates the potential energy used for constrained molecular dynamics.
- * \f$V_i^{con} = c^{bond} \cdot | r_{P(i)} - R_i | + sum_{i \neq j} C^{min} \cdot \tfrac{1}{C_{ij}} + C^{inj} \Bigl (1 - \theta \bigl (\prod_{i \neq j} (P(i) - P(j)) \bigr ) \Bigr )\f$
+ * \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$
  *     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}$ is minimum distance between 
  *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
@@ -1024 +1024 @@
  * \return potential energy
  * \note This routine is scaling quadratically which is not optimal.
+ * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
  */
 double molecule::ConstrainedPotential(ofstream *out, atom **PermutationMap, int startstep, int endstep, double *constants, bool IsAngstroem)
 {
-  double result = 0., tmp;
+  double result = 0., tmp, Norm1, Norm2;
   atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
   Vector trajectory1, trajectory2, normal, TestVector;
@@ -1042 +1043 @@
     tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     result += constants[0] * tmp;
+    //*out << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
     
     // second term: sum of distances to other trajectories
@@ -1047 +1049 @@
     while (Runner->next != end) {
       Runner = Runner->next;
+      if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
+        break;
       // determine normalized trajectories direction vector (n1, n2)
-      Sprinter = PermutationMap[Walker->nr];   // find target point
+      Sprinter = PermutationMap[Walker->nr];   // find first target point
       trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));
       trajectory1.SubtractVector(&Trajectories[Walker].R.at(startstep));
       trajectory1.Normalize();
-      Sprinter = PermutationMap[Runner->nr];   // find target point
+      Norm1 = trajectory1.Norm();
+      Sprinter = PermutationMap[Runner->nr];   // find second target point
       trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));
       trajectory2.SubtractVector(&Trajectories[Runner].R.at(startstep));
       trajectory2.Normalize();
-      // check whether they're linear dependent
-      if (fabs(trajectory1.ScalarProduct(&trajectory2)/trajectory1.Norm()/trajectory2.Norm() - 1.) < MYEPSILON) {
-        *out << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
-        *out << trajectory1;
-        *out << " and ";
-        *out << trajectory2 << endl;
+      Norm2 = trajectory1.Norm();
+      // check whether either is zero()
+      if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
         tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
+      } else if (Norm1 < MYEPSILON) {
+        Sprinter = PermutationMap[Walker->nr];   // find first target point
+        trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy first offset
+        trajectory1.SubtractVector(&Trajectories[Runner].R.at(startstep));  // subtract second offset
+        trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything 
+        trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
+        tmp = trajectory1.Norm();  // remaining norm is distance
+      } else if (Norm2 < MYEPSILON) {
+        Sprinter = PermutationMap[Runner->nr];   // find second target point
+        trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy second offset
+        trajectory2.SubtractVector(&Trajectories[Walker].R.at(startstep));  // subtract first offset
+        trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything 
+        trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
+        tmp = trajectory2.Norm();  // remaining norm is distance
+      } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
+//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
+//        *out << trajectory1;
+//        *out << " and ";
+//        *out << trajectory2;
+        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
+//        *out << " with distance " << tmp << "." << endl;
       } else { // determine distance by finding minimum distance
-        *out << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent." << endl;
+//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
+//        *out << endl;
+//        *out << "First Trajectory: ";
+//        *out << trajectory1 << endl;
+//        *out << "Second Trajectory: ";
+//        *out << trajectory2 << endl;
         // determine normal vector for both
         normal.MakeNormalVector(&trajectory1, &trajectory2);
+        // print all vectors for debugging
+//        *out << "Normal vector in between: ";
+//        *out << normal << endl;
         // setup matrix
         for (int i=NDIM;i--;) {
@@ -1077 +1108 @@
         gsl_linalg_HH_svx(A, x);
         // distance from last component
-        tmp = gsl_vector_get(x,2);   
+        tmp = gsl_vector_get(x,2);
+//        *out << " with distance " << tmp << "." << endl;
         // test whether we really have the intersection (by checking on c_1 and c_2)
         TestVector.CopyVector(&Trajectories[Runner].R.at(startstep));
@@ -1088 +1120 @@
         TestVector.SubtractVector(&trajectory1);
         if (TestVector.Norm() < MYEPSILON) {
-          *out << "Test: ok.\tDistance of " << tmp << " is correct." << endl;  
+//          *out << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;  
         } else {
-          *out << "Test: failed.\tIntersection is off by ";
-          *out << TestVector;
-          *out << "." << endl;  
+//          *out << Verbose(2) << "Test: failed.\tIntersection is off by ";
+//          *out << TestVector;
+//          *out << "." << endl;  
         }
       }
       // add up
       tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
-      result += constants[1] * 1./tmp;
+      if (fabs(tmp) > MYEPSILON) {
+        result += constants[1] * 1./tmp;
+        //*out << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
+      }
     }
       
@@ -1106 +1141 @@
       if ((PermutationMap[Walker->nr] == PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
         Sprinter = PermutationMap[Walker->nr];
-        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
-        *out << Trajectories[Sprinter].R.at(endstep);
-        *out << ", penalting." << endl;
+//        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
+//        *out << Trajectories[Sprinter].R.at(endstep);
+//        *out << ", penalting." << endl;
         result += constants[2];
+        //*out << Verbose(4) << "Adding " << constants[2] << "." << endl;
       }
     }
@@ -1123 +1159 @@
   zeile2 << "                ";
   for (int i=0;i<Nr;i++) {
-    zeile1 << i << "\t";
-    zeile2 << PermutationMap[i]->nr << "\t";
-  }
-  *out << zeile1 << endl << zeile2 << endl;
+    zeile1 << i << " ";
+    zeile2 << PermutationMap[i]->nr << " ";
+  }
+//  *out << zeile1.str() << endl << zeile2.str() << endl;
 };
 
@@ -1137 +1173 @@
  *     -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
  *  -# Next, we only apply transformations that keep the injectivity of the permutations list.
- *  -# Hence, we for one source point we step down the ladder and seek the corresponding owner of this new target 
+ *  -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target 
  *     point and try to change it for one with lesser distance, or for the next one with greater distance, but only
  *     if this decreases the conditional potential.
@@ -1146 +1182 @@
  *  -# Finally, we calculate the potential energy and return.
  * \param *out output stream for debugging
- * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
+ * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
  * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
  * \param endstep step giving final position in constrained MD 
  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
  * \sa molecule::VerletForceIntegration()
- * \return potential energy
- */
-double molecule::MinimiseConstrainedPotential(ofstream *out, ForceMatrix *Force, int startstep, int endstep, bool IsAngstroem)
+ * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
+ * \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
+ *       to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
+ * \bug this all is not O(N log N) but O(N^2)
+ */
+double molecule::MinimiseConstrainedPotential(ofstream *out, atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem)
 {
   double Potential, OldPotential;
-  atom **PermutationMap = (atom **) Malloc(AtomCount*sizeof(atom *), "molecule::MinimiseConstrainedPotential: *PermutationMap");
+  PermutationMap = (atom **) Malloc(AtomCount*sizeof(atom *), "molecule::MinimiseConstrainedPotential: **PermutationMap");
   DistanceMap **DistanceList = (DistanceMap **) Malloc(AtomCount*sizeof(DistanceMap *), "molecule::MinimiseConstrainedPotential: **DistanceList");
-  DistanceMap::iterator *DistanceIterators = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *DistanceIterators"); 
+  DistanceMap::iterator *DistanceIterators = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *DistanceIterators");
+  int *DoubleList = (int *) Malloc(AtomCount*sizeof(int), "molecule::MinimiseConstrainedPotential: *DoubleList");
   double constants[3];
   atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
-
+  DistanceMap::iterator Rider;
+  
   /// Minimise the potential
-  // set constants (TODO)
+  // set Lagrange multiplier constants
   constants[0] = 10.;
   constants[1] = 1.;
   constants[2] = 1e+7;    // just a huge penalty
   // generate the distance list
-  *out << Verbose(1) << "Generating the distance list ... " << endl;
+  *out << Verbose(1) << "Creating the distance list ... " << endl;
   for (int i=AtomCount; i--;) {
-    DistanceList[i] = new DistanceMap;
+    DoubleList[i] = 0;  // store's for how many atoms in startstep this atom is a target in endstep
+    DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
     DistanceList[i]->clear();
   }
+  *out << Verbose(1) << "Filling the distance list ... " << endl;
   Walker = start;
-  while (Walker->next != NULL) {
+  while (Walker->next != end) {
     Walker = Walker->next;
     Runner = start;
-    while(Runner->next != NULL) {
+    while(Runner->next != end) {
       Runner = Runner->next;
       DistanceList[Walker->nr]->insert( DistancePair(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)), Runner) );
     }
   }
-  // create the initial PermutationMap
+  // create the initial PermutationMap (source -> target)
   Walker = start;
-  while (Walker->next != NULL) {
+  while (Walker->next != end) {
     Walker = Walker->next;
-    PermutationMap[Walker->nr] = DistanceList[Walker->nr]->begin()->second;
-    DistanceIterators[Walker->nr] = DistanceList[Walker->nr]->begin();
+    PermutationMap[Walker->nr] = DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
+    DoubleList[DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
+    DistanceIterators[Walker->nr] = DistanceList[Walker->nr]->begin();    // and remember which one we picked
     *out << *Walker << " starts with distance " << DistanceList[Walker->nr]->begin()->first << "." << endl;
   }
-  // make the PermutationMap injective 
+  *out << Verbose(1) << "done." << endl;
+  // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
   *out << Verbose(1) << "Making the PermutationMap injective ... " << endl;
   Walker = start;
-  while ((OldPotential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem)) > constants[2])) {
+  DistanceMap::iterator NewBase;
+  OldPotential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
+  while ((OldPotential) > constants[2]) {
     PrintPermutationMap(out, PermutationMap, AtomCount);
     Walker = Walker->next;
     if (Walker == end) // round-robin at the end
       Walker = start->next;
-    DistanceIterators[Walker->nr]++;  // take next further distance in distance to targets list
-    PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
-  }
+    if (DoubleList[DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
+      continue;
+    // now, try finding a new one
+    NewBase = DistanceIterators[Walker->nr];  // store old base
+    do {
+      NewBase++;  // take next further distance in distance to targets list that's a target of no one
+    } while ((DoubleList[NewBase->second->nr] != 0) && (NewBase != DistanceList[Walker->nr]->end()));
+    if (NewBase != DistanceList[Walker->nr]->end()) {
+      DoubleList[DistanceIterators[Walker->nr]->second->nr]--;  // decrease the old entry in the doubles list
+      DoubleList[NewBase->second->nr]++;    // increase the old entry in the doubles list
+      PermutationMap[Walker->nr] = NewBase->second;
+      Potential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
+      if (Potential > OldPotential) { // undo
+        DoubleList[NewBase->second->nr]--;
+        DistanceIterators[Walker->nr] = NewBase;
+        PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
+        DoubleList[DistanceIterators[Walker->nr]->second->nr]++;
+      } else {
+        OldPotential = Potential;
+        *out << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl;
+      }
+    }
+  } 
+  for (int i=AtomCount; i--;) // now each single entry in the DoubleList should be <=1
+    if (DoubleList[i] > 1) { 
+      cerr << "Failed to create an injective PermutationMap!" << endl;
+      exit(1);
+    }
+  *out << Verbose(1) << "done." << endl;
+  Free((void **)&DoubleList, "molecule::MinimiseConstrainedPotential: *DoubleList");
   // argument minimise the constrained potential in this injective PermutationMap
-  *out << Verbose(1) << "Argument minimising the PermutationMap ... " << endl;
+  *out << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
   do {
-    PrintPermutationMap(out, PermutationMap, AtomCount);
     Walker = start;
     while (Walker->next != end) { // pick one
       Walker = Walker->next;
+      PrintPermutationMap(out, PermutationMap, AtomCount);
       Sprinter = DistanceIterators[Walker->nr]->second;   // store initial partner
       DistanceIterators[Walker->nr]++;  // take next farther distance target
-      PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
-      Runner = start;
-      while(Runner->next != end) { // find the other whose toes were stepping on (this target should be used by another already
+      if (DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->end()) // stop, before we run through the list and still on
+        break;
+      *out << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl; 
+      Runner = start->next;
+      while(Runner != end) { // find the other one whose toes we might be stepping on (this target should be used by another already)
+        if (PermutationMap[Runner->nr] == DistanceIterators[Walker->nr]->second) {
+          *out << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *DistanceIterators[Walker->nr]->second << "." << endl;
+          break;
+        }
         Runner = Runner->next;
-        if (PermutationMap[Runner->nr] == DistanceIterators[Walker->nr]->second)
-          break;
-      }
-      if (Runner->next != end) { // we found one a tripel
+      }
+      if (Runner != end) { // we found one a tripel
         // then look in distance list for Sprinter
-        DistanceMap::iterator Rider = DistanceList[Walker->nr]->begin();
-        for (; Rider != DistanceList[Walker->nr]->end(); Rider++)
+        Rider = DistanceList[Runner->nr]->begin();
+        for (; Rider != DistanceList[Runner->nr]->end(); Rider++)
           if (Rider->second == Sprinter)
             break;
-        if (Rider != DistanceList[Walker->nr]->end()) { // if we have found one
+        if (Rider != DistanceList[Runner->nr]->end()) { // if we have found one
+          *out << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl; 
           // exchange the second also 
-          PermutationMap[Runner->nr] = Rider->second;
+          PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap 
+          PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
           // calculate the new potential
+          PrintPermutationMap(out, PermutationMap, AtomCount);
+          *out << Verbose(2) << "Checking new potential ..." << endl;
           Potential = ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
           if (Potential > OldPotential) { // we made everything worse! Undo ...
+            *out << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl; 
             // Undo for Walker
             DistanceIterators[Walker->nr]--;  // take next farther distance target
@@ -1234 +1317 @@
             // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
             PermutationMap[Runner->nr] = DistanceIterators[Runner->nr]->second;
+            Potential = OldPotential;
           } else {
             DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
+            *out << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl;
             OldPotential = Potential;
-            *out << "Found a better permutation, new potential is " << OldPotential << "." << endl;
             break;
           }
-          if (Potential > constants[2])
+          if (Potential > constants[2]) {
             cerr << "ERROR: The two-step permutation procedure did not maintain injectivity!" << endl;
+            exit(255);
+          }
+          *out << endl;
         }
       } else {
-        DistanceIterators[Walker->nr]--;  // take next farther distance target
-        PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
+        DistanceIterators[Walker->nr]--;  // take old distance target
       }
     }
   } while (Walker->next != end);
-
+  *out << Verbose(1) << "done." << endl;
+
+  
+  /// free memory and return with evaluated potential 
+  for (int i=AtomCount; i--;)
+    DistanceList[i]->clear();
+  Free((void **)&DistanceList, "molecule::MinimiseConstrainedPotential: **DistanceList");
+  Free((void **)&DistanceIterators, "molecule::MinimiseConstrainedPotential: *DistanceIterators");
+  return ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
+};
+
+/** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
+ * \param *out output stream for debugging
+ * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
+ * \param endstep step giving final position in constrained MD 
+ * \param **PermutationMap mapping between the atom label of the initial and the final configuration
+ * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
+ * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
+ */
+void molecule::EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
+{
+  double constant = 10.;
+  atom *Walker = NULL, *Sprinter = NULL; 
+  
   /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
   *out << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl;
@@ -1258 +1367 @@
     // set forces
     for (int i=NDIM;i++;)
-      Force->Matrix[0][Walker->nr][5+i] += 2.*constants[0]*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
+      Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
   }  
-  
-  /// evaluate potential, free memory and return
-  for (int i=AtomCount; i--;)
-    DistanceList[i]->clear();
-  Free((void **)&DistanceList, "molecule::MinimiseConstrainedPotential: **DistanceList");
-  Free((void **)&DistanceIterators, "molecule::MinimiseConstrainedPotential: *DistanceIterators");
+  *out << Verbose(1) << "done." << endl;
+};
+
+/** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
+ * Note, step number is config::MaxOuterStep
+ * \param *out output stream for debugging
+ * \param startstep stating initial configuration in molecule::Trajectories
+ * \param endstep stating final configuration in molecule::Trajectories
+ * \param &config configuration structure
+ * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
+ */
+bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration) 
+{
+  bool status = true;
+  int MaxSteps = configuration.MaxOuterStep; 
+  MoleculeListClass *MoleculePerStep = new MoleculeListClass(MaxSteps+1, AtomCount);
+  // Get the Permutation Map by MinimiseConstrainedPotential
+  atom **PermutationMap = NULL;
+  atom *Walker = NULL, *Sprinter = NULL;
+  MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
+
+  // check whether we have sufficient space in Trajectories for each atom
+  Walker = start;
+  while (Walker->next != end) {
+    Walker = Walker->next;
+    if (Trajectories[Walker].R.size() <= (unsigned int)(MaxSteps)) {
+      //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
+      Trajectories[Walker].R.resize(MaxSteps+1);
+      Trajectories[Walker].U.resize(MaxSteps+1);
+      Trajectories[Walker].F.resize(MaxSteps+1);
+    }
+  }
+  // push endstep to last one
+  Walker = start;
+  while (Walker->next != end) { // remove the endstep (is now the very last one)
+    Walker = Walker->next;
+    for (int n=NDIM;n--;) {
+      Trajectories[Walker].R.at(MaxSteps).x[n] = Trajectories[Walker].R.at(endstep).x[n];
+      Trajectories[Walker].U.at(MaxSteps).x[n] = Trajectories[Walker].U.at(endstep).x[n];
+      Trajectories[Walker].F.at(MaxSteps).x[n] = Trajectories[Walker].F.at(endstep).x[n];
+    }
+  }
+  endstep = MaxSteps;
+  
+  // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
+  *out << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl;
+  for (int step = 0; step <= MaxSteps; step++) {
+    MoleculePerStep->ListOfMolecules[step] = new molecule(elemente);
+    Walker = start;
+    while (Walker->next != end) {
+      Walker = Walker->next;
+      // add to molecule list
+      Sprinter = MoleculePerStep->ListOfMolecules[step]->AddCopyAtom(Walker);
+      for (int n=NDIM;n--;) {
+        Sprinter->x.x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+        // add to Trajectories
+        //*out << Verbose(3) << step << ">=" << MDSteps-1 << endl;
+        //if (step >= MDSteps-1) {
+          Trajectories[Walker].R.at(step).x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
+          Trajectories[Walker].U.at(step).x[n] = 0.;
+          Trajectories[Walker].F.at(step).x[n] = 0.;
+        //}
+      }
+    }
+  }
+  MDSteps = MaxSteps+1;   // otherwise new Trajectories' points aren't stored on save&exit
+  
+  // store the list to single step files
+  int *SortIndex = (int *) Malloc(AtomCount*sizeof(int), "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
+  for (int i=AtomCount; i--; )
+    SortIndex[i] = i;
+  status = MoleculePerStep->OutputConfigForListOfFragments(out, "ConstrainedStep", &configuration, SortIndex, false, false);
+  
+  // free and return
   Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
-  return ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
+  delete(MoleculePerStep);
+  return status;
 };
 
@@ -1280 +1458 @@
  * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential() 
  * \return true - file found and parsed, false - file not found or imparsable
+ * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
  */
 bool molecule::VerletForceIntegration(ofstream *out, char *file, double delta_t, bool IsAngstroem, int DoConstrained)
@@ -1321 +1500 @@
     if (DoConstrained) {
       // calculate forces and potential
-      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, &Force, DoConstrained, 0, IsAngstroem);
+      atom **PermutationMap = NULL;
+      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap, DoConstrained, 0, IsAngstroem);
+      EvaluateConstrainedForces(out, DoConstrained, 0, PermutationMap, &Force);
+      Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
     }
     
@@ -3292 +3474 @@
       
       *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
-      if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
+      if (BondFragments->OutputConfigForListOfFragments(out, FRAGMENTPREFIX, configuration, SortIndex, true, true))
         *out << Verbose(1) << "All configs written." << endl;
       else

src/molecules.hpp

@@ -262 +262 @@
         bool VerletForceIntegration(ofstream *out, char *file, double delta_t, bool IsAngstroem, int DoConstrained);
         double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
-        double MinimiseConstrainedPotential(ofstream *out, ForceMatrix *Force, int startstep, int endstep, bool IsAngstroem);
+        double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
+        void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
+        bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration); 
         
   bool CheckBounds(const Vector *x) const;
@@ -340 +342 @@
   bool AddHydrogenCorrection(ofstream *out, char *path);
   bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
-  bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
+  bool OutputConfigForListOfFragments(ofstream *out, const char *fragmentprefix, config *configuration, int *SortIndex, bool DoPeriodic, bool DoCentering);
   void Output(ofstream *out);
   
@@ -389 +391 @@
     double Deltat;
     int DoConstrainedMD;
+    int MaxOuterStep;
     
     private:
@@ -410 +413 @@
     int Seed;
     
-    int MaxOuterStep;
     int OutVisStep;
     int OutSrcStep;