Changeset 631dcb for src/molecules.cpp


Ignore:
Timestamp:
Jul 23, 2009, 11:23:59 AM (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:
71e7c7
Parents:
b38b64 (diff), fcbfc8 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'Thermostat'

Conflicts:

.gitignore
Makefile.am
molecuilder/src/analyzer.cpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/molecules.cpp
molecuilder/src/molecules.hpp

  • config::SaveMPQC() has different call parameters
  • analyzer and joiner had conflicts due to Chi and ChiPAS values
  • molecule::VerletForceIntegration() is slightly different too, but Thermostat supposedly is old version
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecules.cpp

    rb38b64 r631dcb  
    10091009};
    10101010
     1011/** Evaluates the potential energy used for constrained molecular dynamics.
     1012 * \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$
     1013 *     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
     1014 *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
     1015 * Note that for the second term we have to solve the following linear system:
     1016 * \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,
     1017 * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
     1018 * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
     1019 * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
     1020 * \param *out output stream for debugging
     1021 * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
     1022 * \param startstep start configuration (MDStep in molecule::trajectories)
     1023 * \param endstep end configuration (MDStep in molecule::trajectories)
     1024 * \param *constants constant in front of each term
     1025 * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
     1026 * \return potential energy
     1027 * \note This routine is scaling quadratically which is not optimal.
     1028 * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
     1029 */
     1030double molecule::ConstrainedPotential(ofstream *out, atom **PermutationMap, int startstep, int endstep, double *constants, bool IsAngstroem)
     1031{
     1032  double result = 0., tmp, Norm1, Norm2;
     1033  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
     1034  Vector trajectory1, trajectory2, normal, TestVector;
     1035  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
     1036  gsl_vector *x = gsl_vector_alloc(NDIM);
     1037
     1038  // go through every atom
     1039  Walker = start;
     1040  while (Walker->next != end) {
     1041    Walker = Walker->next;
     1042    // first term: distance to target
     1043    Runner = PermutationMap[Walker->nr];   // find target point
     1044    tmp = (Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)));
     1045    tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     1046    result += constants[0] * tmp;
     1047    //*out << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
     1048   
     1049    // second term: sum of distances to other trajectories
     1050    Runner = start;
     1051    while (Runner->next != end) {
     1052      Runner = Runner->next;
     1053      if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
     1054        break;
     1055      // determine normalized trajectories direction vector (n1, n2)
     1056      Sprinter = PermutationMap[Walker->nr];   // find first target point
     1057      trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));
     1058      trajectory1.SubtractVector(&Trajectories[Walker].R.at(startstep));
     1059      trajectory1.Normalize();
     1060      Norm1 = trajectory1.Norm();
     1061      Sprinter = PermutationMap[Runner->nr];   // find second target point
     1062      trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));
     1063      trajectory2.SubtractVector(&Trajectories[Runner].R.at(startstep));
     1064      trajectory2.Normalize();
     1065      Norm2 = trajectory1.Norm();
     1066      // check whether either is zero()
     1067      if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
     1068        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
     1069      } else if (Norm1 < MYEPSILON) {
     1070        Sprinter = PermutationMap[Walker->nr];   // find first target point
     1071        trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy first offset
     1072        trajectory1.SubtractVector(&Trajectories[Runner].R.at(startstep));  // subtract second offset
     1073        trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
     1074        trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
     1075        tmp = trajectory1.Norm();  // remaining norm is distance
     1076      } else if (Norm2 < MYEPSILON) {
     1077        Sprinter = PermutationMap[Runner->nr];   // find second target point
     1078        trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy second offset
     1079        trajectory2.SubtractVector(&Trajectories[Walker].R.at(startstep));  // subtract first offset
     1080        trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
     1081        trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
     1082        tmp = trajectory2.Norm();  // remaining norm is distance
     1083      } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
     1084//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
     1085//        *out << trajectory1;
     1086//        *out << " and ";
     1087//        *out << trajectory2;
     1088        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
     1089//        *out << " with distance " << tmp << "." << endl;
     1090      } else { // determine distance by finding minimum distance
     1091//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
     1092//        *out << endl;
     1093//        *out << "First Trajectory: ";
     1094//        *out << trajectory1 << endl;
     1095//        *out << "Second Trajectory: ";
     1096//        *out << trajectory2 << endl;
     1097        // determine normal vector for both
     1098        normal.MakeNormalVector(&trajectory1, &trajectory2);
     1099        // print all vectors for debugging
     1100//        *out << "Normal vector in between: ";
     1101//        *out << normal << endl;
     1102        // setup matrix
     1103        for (int i=NDIM;i--;) {
     1104          gsl_matrix_set(A, 0, i, trajectory1.x[i]);
     1105          gsl_matrix_set(A, 1, i, trajectory2.x[i]);
     1106          gsl_matrix_set(A, 2, i, normal.x[i]);
     1107          gsl_vector_set(x,i, (Trajectories[Walker].R.at(startstep).x[i] - Trajectories[Runner].R.at(startstep).x[i]));
     1108        }
     1109        // solve the linear system by Householder transformations
     1110        gsl_linalg_HH_svx(A, x);
     1111        // distance from last component
     1112        tmp = gsl_vector_get(x,2);
     1113//        *out << " with distance " << tmp << "." << endl;
     1114        // test whether we really have the intersection (by checking on c_1 and c_2)
     1115        TestVector.CopyVector(&Trajectories[Runner].R.at(startstep));
     1116        trajectory2.Scale(gsl_vector_get(x,1));
     1117        TestVector.AddVector(&trajectory2);
     1118        normal.Scale(gsl_vector_get(x,2));
     1119        TestVector.AddVector(&normal);
     1120        TestVector.SubtractVector(&Trajectories[Walker].R.at(startstep));
     1121        trajectory1.Scale(gsl_vector_get(x,0));
     1122        TestVector.SubtractVector(&trajectory1);
     1123        if (TestVector.Norm() < MYEPSILON) {
     1124//          *out << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl; 
     1125        } else {
     1126//          *out << Verbose(2) << "Test: failed.\tIntersection is off by ";
     1127//          *out << TestVector;
     1128//          *out << "." << endl; 
     1129        }
     1130      }
     1131      // add up
     1132      tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     1133      if (fabs(tmp) > MYEPSILON) {
     1134        result += constants[1] * 1./tmp;
     1135        //*out << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
     1136      }
     1137    }
     1138     
     1139    // third term: penalty for equal targets
     1140    Runner = start;
     1141    while (Runner->next != end) {
     1142      Runner = Runner->next;
     1143      if ((PermutationMap[Walker->nr] == PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
     1144        Sprinter = PermutationMap[Walker->nr];
     1145//        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
     1146//        *out << Trajectories[Sprinter].R.at(endstep);
     1147//        *out << ", penalting." << endl;
     1148        result += constants[2];
     1149        //*out << Verbose(4) << "Adding " << constants[2] << "." << endl;
     1150      }
     1151    }
     1152  }
     1153 
     1154  return result;
     1155};
     1156
     1157void PrintPermutationMap(ofstream *out, atom **PermutationMap, int Nr)
     1158{
     1159  stringstream zeile1, zeile2;
     1160  int *DoubleList = (int *) Malloc(Nr*sizeof(int), "PrintPermutationMap: *DoubleList");
     1161  int doubles = 0;
     1162  for (int i=0;i<Nr;i++)
     1163    DoubleList[i] = 0;
     1164  zeile1 << "PermutationMap: ";
     1165  zeile2 << "                ";
     1166  for (int i=0;i<Nr;i++) {
     1167    DoubleList[PermutationMap[i]->nr]++;
     1168    zeile1 << i << " ";
     1169    zeile2 << PermutationMap[i]->nr << " ";
     1170  }
     1171  for (int i=0;i<Nr;i++)
     1172    if (DoubleList[i] > 1)
     1173    doubles++;
     1174 // *out << "Found " << doubles << " Doubles." << endl;
     1175  Free((void **)&DoubleList, "PrintPermutationMap: *DoubleList");
     1176//  *out << zeile1.str() << endl << zeile2.str() << endl;
     1177};
     1178
     1179/** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
     1180 * We do the following:
     1181 *  -# Generate a distance list from all source to all target points
     1182 *  -# Sort this per source point
     1183 *  -# Take for each source point the target point with minimum distance, use this as initial permutation
     1184 *  -# check whether molecule::ConstrainedPotential() is greater than injective penalty
     1185 *     -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
     1186 *  -# Next, we only apply transformations that keep the injectivity of the permutations list.
     1187 *  -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target
     1188 *     point and try to change it for one with lesser distance, or for the next one with greater distance, but only
     1189 *     if this decreases the conditional potential.
     1190 *  -# finished.
     1191 *  -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
     1192 *     that the total force is always pointing in direction of the constraint force (ensuring that we move in the
     1193 *     right direction).
     1194 *  -# Finally, we calculate the potential energy and return.
     1195 * \param *out output stream for debugging
     1196 * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
     1197 * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
     1198 * \param endstep step giving final position in constrained MD
     1199 * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
     1200 * \sa molecule::VerletForceIntegration()
     1201 * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
     1202 * \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
     1203 *       to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
     1204 * \bug this all is not O(N log N) but O(N^2)
     1205 */
     1206double molecule::MinimiseConstrainedPotential(ofstream *out, atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem)
     1207{
     1208  double Potential, OldPotential, OlderPotential;
     1209  PermutationMap = (atom **) Malloc(AtomCount*sizeof(atom *), "molecule::MinimiseConstrainedPotential: **PermutationMap");
     1210  DistanceMap **DistanceList = (DistanceMap **) Malloc(AtomCount*sizeof(DistanceMap *), "molecule::MinimiseConstrainedPotential: **DistanceList");
     1211  DistanceMap::iterator *DistanceIterators = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *DistanceIterators");
     1212  int *DoubleList = (int *) Malloc(AtomCount*sizeof(int), "molecule::MinimiseConstrainedPotential: *DoubleList");
     1213  DistanceMap::iterator *StepList = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *StepList");
     1214  double constants[3];
     1215  int round;
     1216  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
     1217  DistanceMap::iterator Rider, Strider;
     1218 
     1219  /// Minimise the potential
     1220  // set Lagrange multiplier constants
     1221  constants[0] = 10.;
     1222  constants[1] = 1.;
     1223  constants[2] = 1e+7;    // just a huge penalty
     1224  // generate the distance list
     1225  *out << Verbose(1) << "Creating the distance list ... " << endl;
     1226  for (int i=AtomCount; i--;) {
     1227    DoubleList[i] = 0;  // stores for how many atoms in startstep this atom is a target in endstep
     1228    DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
     1229    DistanceList[i]->clear();
     1230  }
     1231  *out << Verbose(1) << "Filling the distance list ... " << endl;
     1232  Walker = start;
     1233  while (Walker->next != end) {
     1234    Walker = Walker->next;
     1235    Runner = start;
     1236    while(Runner->next != end) {
     1237      Runner = Runner->next;
     1238      DistanceList[Walker->nr]->insert( DistancePair(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)), Runner) );
     1239    }
     1240  }
     1241  // create the initial PermutationMap (source -> target)
     1242  Walker = start;
     1243  while (Walker->next != end) {
     1244    Walker = Walker->next;
     1245    StepList[Walker->nr] = DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
     1246    PermutationMap[Walker->nr] = DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
     1247    DoubleList[DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
     1248    DistanceIterators[Walker->nr] = DistanceList[Walker->nr]->begin();    // and remember which one we picked
     1249    *out << *Walker << " starts with distance " << DistanceList[Walker->nr]->begin()->first << "." << endl;
     1250  }
     1251  *out << Verbose(1) << "done." << endl;
     1252  // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
     1253  *out << Verbose(1) << "Making the PermutationMap injective ... " << endl;
     1254  Walker = start;
     1255  DistanceMap::iterator NewBase;
     1256  OldPotential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
     1257  while ((OldPotential) > constants[2]) {
     1258    PrintPermutationMap(out, PermutationMap, AtomCount);
     1259    Walker = Walker->next;
     1260    if (Walker == end) // round-robin at the end
     1261      Walker = start->next;
     1262    if (DoubleList[DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
     1263      continue;
     1264    // now, try finding a new one
     1265    NewBase = DistanceIterators[Walker->nr];  // store old base
     1266    do {
     1267      NewBase++;  // take next further distance in distance to targets list that's a target of no one
     1268    } while ((DoubleList[NewBase->second->nr] != 0) && (NewBase != DistanceList[Walker->nr]->end()));
     1269    if (NewBase != DistanceList[Walker->nr]->end()) {
     1270      PermutationMap[Walker->nr] = NewBase->second;
     1271      Potential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
     1272      if (Potential > OldPotential) { // undo
     1273        PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
     1274      } else {  // do
     1275        DoubleList[DistanceIterators[Walker->nr]->second->nr]--;  // decrease the old entry in the doubles list
     1276        DoubleList[NewBase->second->nr]++;    // increase the old entry in the doubles list
     1277        DistanceIterators[Walker->nr] = NewBase;
     1278        OldPotential = Potential;
     1279        *out << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl;
     1280      }
     1281    }
     1282  }
     1283  for (int i=AtomCount; i--;) // now each single entry in the DoubleList should be <=1
     1284    if (DoubleList[i] > 1) {
     1285      cerr << "Failed to create an injective PermutationMap!" << endl;
     1286      exit(1);
     1287    }
     1288  *out << Verbose(1) << "done." << endl;
     1289  Free((void **)&DoubleList, "molecule::MinimiseConstrainedPotential: *DoubleList");
     1290  // argument minimise the constrained potential in this injective PermutationMap
     1291  *out << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
     1292  OldPotential = 1e+10;
     1293  round = 0;
     1294  do {
     1295    *out << "Starting round " << ++round << " ... " << endl;
     1296    OlderPotential = OldPotential;
     1297    do {
     1298      Walker = start;
     1299      while (Walker->next != end) { // pick one
     1300        Walker = Walker->next;
     1301        PrintPermutationMap(out, PermutationMap, AtomCount);
     1302        Sprinter = DistanceIterators[Walker->nr]->second;   // store initial partner
     1303        Strider = DistanceIterators[Walker->nr];  //remember old iterator
     1304        DistanceIterators[Walker->nr] = StepList[Walker->nr]; 
     1305        if (DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
     1306          DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->begin();
     1307          break;
     1308        }
     1309        //*out << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
     1310        // find source of the new target
     1311        Runner = start->next;
     1312        while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
     1313          if (PermutationMap[Runner->nr] == DistanceIterators[Walker->nr]->second) {
     1314            //*out << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
     1315            break;
     1316          }
     1317          Runner = Runner->next;
     1318        }
     1319        if (Runner != end) { // we found the other source
     1320          // then look in its distance list for Sprinter
     1321          Rider = DistanceList[Runner->nr]->begin();
     1322          for (; Rider != DistanceList[Runner->nr]->end(); Rider++)
     1323            if (Rider->second == Sprinter)
     1324              break;
     1325          if (Rider != DistanceList[Runner->nr]->end()) { // if we have found one
     1326            //*out << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl;
     1327            // exchange both
     1328            PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
     1329            PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
     1330            PrintPermutationMap(out, PermutationMap, AtomCount);
     1331            // calculate the new potential
     1332            //*out << Verbose(2) << "Checking new potential ..." << endl;
     1333            Potential = ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
     1334            if (Potential > OldPotential) { // we made everything worse! Undo ...
     1335              //*out << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
     1336              //*out << Verbose(3) << "Setting " << *Runner << "'s source to " << *DistanceIterators[Runner->nr]->second << "." << endl;
     1337              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
     1338              PermutationMap[Runner->nr] = DistanceIterators[Runner->nr]->second;
     1339              // Undo for Walker
     1340              DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
     1341              //*out << Verbose(3) << "Setting " << *Walker << "'s source to " << *DistanceIterators[Walker->nr]->second << "." << endl;
     1342              PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
     1343            } else {
     1344              DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
     1345              *out << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl;
     1346              OldPotential = Potential;
     1347            }
     1348            if (Potential > constants[2]) {
     1349              cerr << "ERROR: The two-step permutation procedure did not maintain injectivity!" << endl;
     1350              exit(255);
     1351            }
     1352            //*out << endl;
     1353          } else {
     1354            cerr << "ERROR: " << *Runner << " was not the owner of " << *Sprinter << "!" << endl;
     1355            exit(255);
     1356          }
     1357        } else {
     1358          PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // new target has no source!
     1359        }
     1360        StepList[Walker->nr]++; // take next farther distance target
     1361      }
     1362    } while (Walker->next != end);
     1363  } while ((OlderPotential - OldPotential) > 1e-3);
     1364  *out << Verbose(1) << "done." << endl;
     1365
     1366 
     1367  /// free memory and return with evaluated potential
     1368  for (int i=AtomCount; i--;)
     1369    DistanceList[i]->clear();
     1370  Free((void **)&DistanceList, "molecule::MinimiseConstrainedPotential: **DistanceList");
     1371  Free((void **)&DistanceIterators, "molecule::MinimiseConstrainedPotential: *DistanceIterators");
     1372  return ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
     1373};
     1374
     1375/** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
     1376 * \param *out output stream for debugging
     1377 * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
     1378 * \param endstep step giving final position in constrained MD
     1379 * \param **PermutationMap mapping between the atom label of the initial and the final configuration
     1380 * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
     1381 * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
     1382 */
     1383void molecule::EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
     1384{
     1385  double constant = 10.;
     1386  atom *Walker = NULL, *Sprinter = NULL;
     1387 
     1388  /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
     1389  *out << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl;
     1390  Walker = start;
     1391  while (Walker->next != NULL) {
     1392    Walker = Walker->next;
     1393    Sprinter = PermutationMap[Walker->nr];
     1394    // set forces
     1395    for (int i=NDIM;i++;)
     1396      Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
     1397  } 
     1398  *out << Verbose(1) << "done." << endl;
     1399};
     1400
     1401/** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
     1402 * Note, step number is config::MaxOuterStep
     1403 * \param *out output stream for debugging
     1404 * \param startstep stating initial configuration in molecule::Trajectories
     1405 * \param endstep stating final configuration in molecule::Trajectories
     1406 * \param &config configuration structure
     1407 * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
     1408 */
     1409bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration)
     1410{
     1411  bool status = true;
     1412  int MaxSteps = configuration.MaxOuterStep;
     1413  MoleculeListClass *MoleculePerStep = new MoleculeListClass(MaxSteps+1, AtomCount);
     1414  // Get the Permutation Map by MinimiseConstrainedPotential
     1415  atom **PermutationMap = NULL;
     1416  atom *Walker = NULL, *Sprinter = NULL;
     1417  MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
     1418
     1419  // check whether we have sufficient space in Trajectories for each atom
     1420  Walker = start;
     1421  while (Walker->next != end) {
     1422    Walker = Walker->next;
     1423    if (Trajectories[Walker].R.size() <= (unsigned int)(MaxSteps)) {
     1424      //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
     1425      Trajectories[Walker].R.resize(MaxSteps+1);
     1426      Trajectories[Walker].U.resize(MaxSteps+1);
     1427      Trajectories[Walker].F.resize(MaxSteps+1);
     1428    }
     1429  }
     1430  // push endstep to last one
     1431  Walker = start;
     1432  while (Walker->next != end) { // remove the endstep (is now the very last one)
     1433    Walker = Walker->next;
     1434    for (int n=NDIM;n--;) {
     1435      Trajectories[Walker].R.at(MaxSteps).x[n] = Trajectories[Walker].R.at(endstep).x[n];
     1436      Trajectories[Walker].U.at(MaxSteps).x[n] = Trajectories[Walker].U.at(endstep).x[n];
     1437      Trajectories[Walker].F.at(MaxSteps).x[n] = Trajectories[Walker].F.at(endstep).x[n];
     1438    }
     1439  }
     1440  endstep = MaxSteps;
     1441 
     1442  // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
     1443  *out << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl;
     1444  for (int step = 0; step <= MaxSteps; step++) {
     1445    MoleculePerStep->ListOfMolecules[step] = new molecule(elemente);
     1446    Walker = start;
     1447    while (Walker->next != end) {
     1448      Walker = Walker->next;
     1449      // add to molecule list
     1450      Sprinter = MoleculePerStep->ListOfMolecules[step]->AddCopyAtom(Walker);
     1451      for (int n=NDIM;n--;) {
     1452        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);
     1453        // add to Trajectories
     1454        //*out << Verbose(3) << step << ">=" << MDSteps-1 << endl;
     1455        if (step < MaxSteps) {
     1456          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);
     1457          Trajectories[Walker].U.at(step).x[n] = 0.;
     1458          Trajectories[Walker].F.at(step).x[n] = 0.;
     1459        }
     1460      }
     1461    }
     1462  }
     1463  MDSteps = MaxSteps+1;   // otherwise new Trajectories' points aren't stored on save&exit
     1464 
     1465  // store the list to single step files
     1466  int *SortIndex = (int *) Malloc(AtomCount*sizeof(int), "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
     1467  for (int i=AtomCount; i--; )
     1468    SortIndex[i] = i;
     1469  status = MoleculePerStep->OutputConfigForListOfFragments(out, "ConstrainedStep", &configuration, SortIndex, false, false);
     1470 
     1471  // free and return
     1472  Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
     1473  delete(MoleculePerStep);
     1474  return status;
     1475};
     1476
    10111477/** Parses nuclear forces from file and performs Verlet integration.
    10121478 * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
    10131479 * have to transform them).
    10141480 * This adds a new MD step to the config file.
     1481 * \param *out output stream for debugging
    10151482 * \param *file filename
     1483 * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
    10161484 * \param delta_t time step width in atomic units
    10171485 * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
     1486 * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential()
    10181487 * \return true - file found and parsed, false - file not found or imparsable
    1019  */
    1020 bool molecule::VerletForceIntegration(char *file, double delta_t, bool IsAngstroem)
    1021 {
    1022   element *runner = elemente->start;
     1488 * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
     1489 */
     1490bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
     1491{
    10231492  atom *walker = NULL;
    1024   int AtomNo;
    10251493  ifstream input(file);
    10261494  string token;
    10271495  stringstream item;
    1028   double a, IonMass;
     1496  double IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp;
    10291497  ForceMatrix Force;
    1030   Vector tmpvector;
    10311498
    10321499  CountElements();  // make sure ElementsInMolecule is up to date
    … …  
    10461513    }
    10471514    // correct Forces
    1048 //    for(int d=0;d<NDIM;d++)
    1049 //      tmpvector.x[d] = 0.;
    1050 //    for(int i=0;i<AtomCount;i++)
    1051 //      for(int d=0;d<NDIM;d++) {
    1052 //        tmpvector.x[d] += Force.Matrix[0][i][d+5];
    1053 //      }
    1054 //    for(int i=0;i<AtomCount;i++)
    1055 //      for(int d=0;d<NDIM;d++) {
    1056 //        Force.Matrix[0][i][d+5] -= tmpvector.x[d]/(double)AtomCount;
    1057 //      }
     1515    for(int d=0;d<NDIM;d++)
     1516      Vector[d] = 0.;
     1517    for(int i=0;i<AtomCount;i++)
     1518      for(int d=0;d<NDIM;d++) {
     1519        Vector[d] += Force.Matrix[0][i][d+5];
     1520      }
     1521    for(int i=0;i<AtomCount;i++)
     1522      for(int d=0;d<NDIM;d++) {
     1523        Force.Matrix[0][i][d+5] -= Vector[d]/(double)AtomCount;
     1524      }
     1525    // solve a constrained potential if we are meant to
     1526    if (configuration.DoConstrainedMD) {
     1527      // calculate forces and potential
     1528      atom **PermutationMap = NULL;
     1529      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
     1530      EvaluateConstrainedForces(out, configuration.DoConstrainedMD, 0, PermutationMap, &Force);
     1531      Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
     1532    }
     1533   
    10581534    // and perform Verlet integration for each atom with position, velocity and force vector
    1059     runner = elemente->start;
    1060     while (runner->next != elemente->end) { // go through every element
    1061       runner = runner->next;
    1062       IonMass = runner->mass;
    1063       a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    1064       if (ElementsInMolecule[runner->Z]) { // if this element got atoms
    1065         AtomNo = 0;
     1535    walker = start;
     1536    while (walker->next != end) { // go through every atom of this element
     1537      walker = walker->next;
     1538      //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
     1539      // check size of vectors
     1540      if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
     1541        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
     1542        Trajectories[walker].R.resize(MDSteps+10);
     1543        Trajectories[walker].U.resize(MDSteps+10);
     1544        Trajectories[walker].F.resize(MDSteps+10);
     1545      }
     1546     
     1547      // Update R (and F)
     1548      for (int d=0; d<NDIM; d++) {
     1549        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
     1550        Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
     1551        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
     1552        Trajectories[walker].R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(Trajectories[walker].F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
     1553      }
     1554      // Update U
     1555      for (int d=0; d<NDIM; d++) {
     1556        Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
     1557        Trajectories[walker].U.at(MDSteps).x[d] += configuration.Deltat * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
     1558      }
     1559//      out << "Integrated position&velocity of step " << (MDSteps) << ": (";
     1560//      for (int d=0;d<NDIM;d++)
     1561//        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
     1562//      out << ")\t(";
     1563//      for (int d=0;d<NDIM;d++)
     1564//        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
     1565//      out << ")" << endl;
     1566            // next atom
     1567    }
     1568  }
     1569  // correct velocities (rather momenta) so that center of mass remains motionless
     1570  for(int d=0;d<NDIM;d++)
     1571    Vector[d] = 0.;
     1572  IonMass = 0.;
     1573  walker = start;
     1574  while (walker->next != end) { // go through every atom
     1575    walker = walker->next;
     1576    IonMass += walker->type->mass;  // sum up total mass
     1577    for(int d=0;d<NDIM;d++) {
     1578      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
     1579    }
     1580  }
     1581  // correct velocities (rather momenta) so that center of mass remains motionless
     1582  for(int d=0;d<NDIM;d++)
     1583    Vector[d] /= IonMass;
     1584  ActualTemp = 0.;
     1585  walker = start;
     1586  while (walker->next != end) { // go through every atom of this element
     1587    walker = walker->next;
     1588    for(int d=0;d<NDIM;d++) {
     1589      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d];
     1590      ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
     1591    }
     1592  }
     1593  Thermostats(configuration, ActualTemp, Berendsen);
     1594  MDSteps++;
     1595
     1596
     1597  // exit
     1598  return true;
     1599};
     1600
     1601/** Implementation of various thermostats.
     1602 * All these thermostats apply an additional force which has the following forms:
     1603 * -# Woodcock
     1604 *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
     1605 * -# Gaussian
     1606 *  \f$ \frac{ \sum_i \frac{p_i}{m_i} \frac{\partial V}{\partial q_i}} {\sum_i \frac{p^2_i}{m_i}} \cdot p_i\f$
     1607 * -# Langevin
     1608 *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$ 
     1609 * -# Berendsen
     1610 *  \f$p_i \rightarrow \left [ 1+ \frac{\delta t}{\tau_T} \left ( \frac{T_0}{T} \right ) \right ]^{\frac{1}{2}} \cdot p_i\f$
     1611 * -# Nose-Hoover
     1612 *  \f$\zeta p_i \f$ with \f$\frac{\partial \zeta}{\partial t} = \frac{1}{M_s} \left ( \sum^N_{i=1} \frac{p_i^2}{m_i} - g k_B T \right )\f$
     1613 * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
     1614 * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
     1615 * belatedly and the constraint force set.
     1616 * \param *P Problem at hand
     1617 * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
     1618 * \sa InitThermostat()
     1619 */
     1620void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
     1621{
     1622  double ekin = 0.;
     1623  double E = 0., G = 0.;
     1624  double delta_alpha = 0.;
     1625  double ScaleTempFactor;
     1626  double sigma;
     1627  double IonMass;
     1628  int d;
     1629  gsl_rng * r;
     1630  const gsl_rng_type * T;
     1631  double *U = NULL, *F = NULL, FConstraint[NDIM];
     1632  atom *walker = NULL;
     1633 
     1634  // calculate scale configuration
     1635  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
     1636   
     1637  // differentating between the various thermostats
     1638  switch(Thermostat) {
     1639     case None:
     1640      cout << Verbose(2) <<  "Applying no thermostat..." << endl;
     1641      break;
     1642     case Woodcock:
     1643      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
     1644        cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
    10661645        walker = start;
    10671646        while (walker->next != end) { // go through every atom of this element
    10681647          walker = walker->next;
    1069           if (walker->type == runner) { // if this atom fits to element
    1070             // check size of vectors
    1071             if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
    1072               //cout << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
    1073               Trajectories[walker].R.resize(MDSteps+10);
    1074               Trajectories[walker].U.resize(MDSteps+10);
    1075               Trajectories[walker].F.resize(MDSteps+10);
     1648          IonMass = walker->type->mass;
     1649          U = Trajectories[walker].U.at(MDSteps).x;
     1650          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1651            for (d=0; d<NDIM; d++) {
     1652              U[d] *= sqrt(ScaleTempFactor);
     1653              ekin += 0.5*IonMass * U[d]*U[d];
    10761654            }
    1077             // 1. calculate x(t+\delta t)
    1078             for (int d=0; d<NDIM; d++) {
    1079               Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5];
    1080               Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
    1081               Trajectories[walker].R.at(MDSteps).x[d] += delta_t*(Trajectories[walker].U.at(MDSteps-1).x[d]);
    1082               Trajectories[walker].R.at(MDSteps).x[d] += 0.5*delta_t*delta_t*(Trajectories[walker].F.at(MDSteps-1).x[d])/IonMass;     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
     1655        }
     1656      }
     1657      break;
     1658     case Gaussian:
     1659      cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
     1660      walker = start;
     1661      while (walker->next != end) { // go through every atom of this element
     1662        walker = walker->next;
     1663        IonMass = walker->type->mass;
     1664        U = Trajectories[walker].U.at(MDSteps).x;
     1665        F = Trajectories[walker].F.at(MDSteps).x;
     1666        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1667          for (d=0; d<NDIM; d++) {
     1668            G += U[d] * F[d];
     1669            E += U[d]*U[d]*IonMass;
     1670          }
     1671      }
     1672      cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
     1673      walker = start;
     1674      while (walker->next != end) { // go through every atom of this element
     1675        walker = walker->next;
     1676        IonMass = walker->type->mass;
     1677        U = Trajectories[walker].U.at(MDSteps).x;
     1678        F = Trajectories[walker].F.at(MDSteps).x;
     1679        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1680          for (d=0; d<NDIM; d++) {
     1681            FConstraint[d] = (G/E) * (U[d]*IonMass);
     1682            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
     1683            ekin += IonMass * U[d]*U[d];
     1684          }
     1685      }
     1686      break;
     1687     case Langevin:
     1688      cout << Verbose(2) <<  "Applying Langevin thermostat..." << endl;
     1689      // init random number generator
     1690      gsl_rng_env_setup();
     1691      T = gsl_rng_default;
     1692      r = gsl_rng_alloc (T);
     1693      // Go through each ion
     1694      walker = start;
     1695      while (walker->next != end) { // go through every atom of this element
     1696        walker = walker->next;
     1697        IonMass = walker->type->mass;
     1698        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
     1699        U = Trajectories[walker].U.at(MDSteps).x;
     1700        F = Trajectories[walker].F.at(MDSteps).x;
     1701        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1702          // throw a dice to determine whether it gets hit by a heat bath particle
     1703          if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) { 
     1704            cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
     1705            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
     1706            for (d=0; d<NDIM; d++) {
     1707              U[d] = gsl_ran_gaussian (r, sigma);
    10831708            }
    1084             // 2. Calculate v(t+\delta t)
    1085             for (int d=0; d<NDIM; d++) {
    1086               Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
    1087               Trajectories[walker].U.at(MDSteps).x[d] += 0.5*delta_t*(Trajectories[walker].F.at(MDSteps-1).x[d]+Trajectories[walker].F.at(MDSteps).x[d])/IonMass;
    1088             }
    1089 //            cout << "Integrated position&velocity of step " << (MDSteps) << ": (";
    1090 //            for (int d=0;d<NDIM;d++)
    1091 //              cout << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
    1092 //            cout << ")\t(";
    1093 //            for (int d=0;d<NDIM;d++)
    1094 //              cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
    1095 //            cout << ")" << endl;
    1096             // next atom
    1097             AtomNo++;
     1709            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
     1710          }
     1711          for (d=0; d<NDIM; d++)
     1712            ekin += 0.5*IonMass * U[d]*U[d];
     1713        }
     1714      }
     1715      break;
     1716     case Berendsen:
     1717      cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
     1718      walker = start;
     1719      while (walker->next != end) { // go through every atom of this element
     1720        walker = walker->next;
     1721        IonMass = walker->type->mass;
     1722        U = Trajectories[walker].U.at(MDSteps).x;
     1723        F = Trajectories[walker].F.at(MDSteps).x;
     1724        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1725          for (d=0; d<NDIM; d++) {
     1726            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
     1727            ekin += 0.5*IonMass * U[d]*U[d];
    10981728          }
    10991729        }
    11001730      }
    1101     }
    1102   }
    1103 //  // correct velocities (rather momenta) so that center of mass remains motionless
    1104 //  tmpvector.zero()
    1105 //  IonMass = 0.;
    1106 //  walker = start;
    1107 //  while (walker->next != end) { // go through every atom
    1108 //    walker = walker->next;
    1109 //    IonMass += walker->type->mass;  // sum up total mass
    1110 //    for(int d=0;d<NDIM;d++) {
    1111 //      tmpvector.x[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
    1112 //    }
    1113 //  }
    1114 //  walker = start;
    1115 //  while (walker->next != end) { // go through every atom of this element
    1116 //    walker = walker->next;
    1117 //    for(int d=0;d<NDIM;d++) {
    1118 //      Trajectories[walker].U.at(MDSteps).x[d] -= tmpvector.x[d]*walker->type->mass/IonMass;
    1119 //    }
    1120 //  }
    1121   MDSteps++;
    1122 
    1123 
    1124   // exit
    1125   return true;
    1126 };
    1127 
    1128 /** Align all atoms in such a manner that given vector \a *n is along z axis.
     1731      break;
     1732     case NoseHoover:
     1733      cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
     1734      // dynamically evolve alpha (the additional degree of freedom)
     1735      delta_alpha = 0.;
     1736      walker = start;
     1737      while (walker->next != end) { // go through every atom of this element
     1738        walker = walker->next;
     1739        IonMass = walker->type->mass;
     1740        U = Trajectories[walker].U.at(MDSteps).x;
     1741        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1742          for (d=0; d<NDIM; d++) {
     1743            delta_alpha += U[d]*U[d]*IonMass;
     1744          }
     1745        }
     1746      }
     1747      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
     1748      configuration.alpha += delta_alpha*configuration.Deltat;
     1749      cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
     1750      // apply updated alpha as additional force
     1751      walker = start;
     1752      while (walker->next != end) { // go through every atom of this element
     1753        walker = walker->next;
     1754        IonMass = walker->type->mass;
     1755        U = Trajectories[walker].U.at(MDSteps).x;
     1756        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1757          for (d=0; d<NDIM; d++) {
     1758              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
     1759              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
     1760              ekin += (0.5*IonMass) * U[d]*U[d];
     1761            }
     1762        }
     1763      }
     1764      break;
     1765  }   
     1766  cout << Verbose(1) << "Kinetic energy is " << ekin << "." << endl;
     1767};
     1768
     1769/** Align all atoms in such a manner that given vector \a *n is along z axis.
    11291770 * \param n[] alignment vector.
    11301771 */
    … …  
    17642405  Vector x;
    17652406  int FalseBondDegree = 0;
    1766 
     2407 
    17672408  BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
    17682409  *out << Verbose(0) << "Begin of CreateAdjacencyList." << endl;
    … …  
    19232564                        *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
    19242565          *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
    1925 
     2566               
    19262567          // output bonds for debugging (if bond chain list was correctly installed)
    19272568          *out << Verbose(1) << endl << "From contents of bond chain list:";
    … …  
    21602801    ColorList[i] = white;
    21612802  }
    2162 
     2803 
    21632804  *out << Verbose(1) << "Back edge list - ";
    21642805  BackEdgeStack->Output(out);
    … …  
    30093650  while (MolecularWalker->next != NULL) {
    30103651    MolecularWalker = MolecularWalker->next;
     3652    *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    30113653    LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
    30123654//    // check the list of local atoms for debugging
    … …  
    30243666    delete(LocalBackEdgeStack);
    30253667  }
    3026 
     3668 
    30273669  // ===== 3. if structure still valid, parse key set file and others =====
    30283670  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
    … …  
    30303672  // ===== 4. check globally whether there's something to do actually (first adaptivity check)
    30313673  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
    3032 
    3033   // =================================== Begin of FRAGMENTATION ===============================
    3034   // ===== 6a. assign each keyset to its respective subgraph =====
     3674 
     3675  // =================================== Begin of FRAGMENTATION =============================== 
     3676  // ===== 6a. assign each keyset to its respective subgraph ===== 
    30353677  Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);
    30363678
    … …  
    30673709  delete(ParsedFragmentList);
    30683710  delete[](MinimumRingSize);
    3069 
     3711 
    30703712
    30713713  // ==================================== End of FRAGMENTATION ============================================
    … …  
    31053747
    31063748      *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
    3107       if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
     3749      if (BondFragments->OutputConfigForListOfFragments(out, FRAGMENTPREFIX, configuration, SortIndex, true, true))
    31083750        *out << Verbose(1) << "All configs written." << endl;
    31093751      else
    … …  
    31603802  atom *Walker = NULL, *OtherAtom = NULL;
    31613803  ReferenceStack->Push(Binder);
    3162 
     3804 
    31633805  do {  // go through all bonds and push local ones
    31643806    Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
    … …  
    31733815        }
    31743816    Binder = ReferenceStack->PopFirst();  // loop the stack for next item
    3175                 *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
     3817    *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
    31763818    ReferenceStack->Push(Binder);
    31773819  } while (FirstBond != Binder);
    3178 
     3820 
    31793821  return status;
    31803822};
    … …  
    33213963  Walker = start;
    33223964  while (Walker->next != end) {
    3323     Walker = Walker->next;
     3965    Walker = Walker->next; 
    33243966    *out << Verbose(4) << "Atom " << Walker->Name << "/" << Walker->nr << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
    33253967    TotalDegree = 0;
Note: See TracChangeset for help on using the changeset viewer.