Changeset 36ec71 for src/molecules.cpp


Ignore:
Timestamp:
Jul 24, 2009, 10:38:32 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:
d30402
Parents:
042f82 (diff), 51c910 (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.
git-author:
Frederik Heber <heber@…> (07/23/09 14:23:32)
git-committer:
Frederik Heber <heber@…> (07/24/09 10:38:32)
Message:

Merge branch 'master' into ConcaveHull

Conflicts:

molecuilder/src/analyzer.cpp
molecuilder/src/bond.cpp
molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/datacreator.cpp
molecuilder/src/datacreator.hpp
molecuilder/src/defs.hpp
molecuilder/src/ellipsoid.cpp
molecuilder/src/joiner.cpp
molecuilder/src/molecules.cpp
molecuilder/src/molecules.hpp
molecuilder/src/parser.cpp
molecuilder/src/parser.hpp

merges:

compilation fixes:

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/molecules.cpp

    r042f82 r36ec71  
    6363  cell_size[1] = cell_size[3] = cell_size[4]= 0.;
    6464  strcpy(name,"none");
     65  IndexNr  = -1;
     66  ActiveFlag = false;
    6567};
    6668
    … …  
    602604 * \param *filename filename
    603605 */
    604 void molecule::SetNameFromFilename(char *filename)
     606void molecule::SetNameFromFilename(const char *filename)
    605607{
    606608  int length = 0;
    607609  char *molname = strrchr(filename, '/')+sizeof(char);  // search for filename without dirs
    608   char *endname = strrchr(filename, '.');
     610  char *endname = strchr(molname, '.');
    609611  if ((endname == NULL) || (endname < molname))
    610612    length = strlen(molname);
    … …  
    612614    length = strlen(molname) - strlen(endname);
    613615  strncpy(name, molname, length);
     616  name[length]='\0';
    614617};
    615618
    … …  
    669672        }
    670673      }
    671       // matrix multiply
    672       x.MatrixMultiplication(M);
    673       ptr->x.CopyVector(&x);
    674674    }
    675675  }
    … …  
    710710    max->AddVector(min);
    711711    Translate(min);
     712    Center.Zero();
    712713  }
    713714  delete(min);
    … …  
    719720 * \param *center return vector for translation vector
    720721 */
    721 void molecule::CenterOrigin(ofstream *out, Vector *center)
     722void molecule::CenterOrigin(ofstream *out)
    722723{
    723724  int Num = 0;
    724725  atom *ptr = start->next;  // start at first in list
    725726
    726   for(int i=NDIM;i--;) // zero center vector
    727     center->x[i] = 0.;
     727  Center.Zero();
    728728
    729729  if (ptr != end) {   //list not empty?
    … …  
    731731      ptr = ptr->next;
    732732      Num++;
    733       center->AddVector(&ptr->x);
    734     }
    735     center->Scale(-1./Num); // divide through total number (and sign for direction)
    736     Translate(center);
     733      Center.AddVector(&ptr->x);
     734    }
     735    Center.Scale(-1./Num); // divide through total number (and sign for direction)
     736    Translate(&Center);
     737    Center.Zero();
    737738  }
    738739};
    … …  
    799800 * \param *center return vector for translation vector
    800801 */
    801 void molecule::CenterGravity(ofstream *out, Vector *center)
    802 {
    803   if (center == NULL) {
    804     DetermineCenter(*center);
    805     Translate(center);
    806     delete(center);
    807   } else {
    808     Translate(center);
    809   }
    810 };
     802void molecule::CenterPeriodic(ofstream *out)
     803{
     804  DeterminePeriodicCenter(Center);
     805};
     806
     807
     808/** Centers the center of gravity of the atoms at (0,0,0).
     809 * \param *out output stream for debugging
     810 * \param *center return vector for translation vector
     811 */
     812void molecule::CenterAtVector(ofstream *out, Vector *newcenter)
     813{
     814  Center.CopyVector(newcenter);
     815};
     816
    811817
    812818/** Scales all atoms by \a *factor.
    … …  
    911917
    912918/** Determines center of molecule (yet not considering atom masses).
    913  * \param Center reference to return vector
    914  */
    915 void molecule::DetermineCenter(Vector &Center)
     919 * \param center reference to return vector
     920 */
     921void molecule::DeterminePeriodicCenter(Vector &center)
    916922{
    917923  atom *Walker = start;
    … …  
    987993  Vector *CenterOfGravity = DetermineCenterOfGravity(out);
    988994
    989   CenterGravity(out, CenterOfGravity);
     995  CenterPeriodic(out);
    990996
    991997  // reset inertia tensor
    … …  
    10831089};
    10841090
     1091/** Evaluates the potential energy used for constrained molecular dynamics.
     1092 * \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$
     1093 *     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
     1094 *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
     1095 * Note that for the second term we have to solve the following linear system:
     1096 * \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,
     1097 * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
     1098 * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
     1099 * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
     1100 * \param *out output stream for debugging
     1101 * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
     1102 * \param startstep start configuration (MDStep in molecule::trajectories)
     1103 * \param endstep end configuration (MDStep in molecule::trajectories)
     1104 * \param *constants constant in front of each term
     1105 * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
     1106 * \return potential energy
     1107 * \note This routine is scaling quadratically which is not optimal.
     1108 * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
     1109 */
     1110double molecule::ConstrainedPotential(ofstream *out, atom **PermutationMap, int startstep, int endstep, double *constants, bool IsAngstroem)
     1111{
     1112  double result = 0., tmp, Norm1, Norm2;
     1113  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
     1114  Vector trajectory1, trajectory2, normal, TestVector;
     1115  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
     1116  gsl_vector *x = gsl_vector_alloc(NDIM);
     1117
     1118  // go through every atom
     1119  Walker = start;
     1120  while (Walker->next != end) {
     1121    Walker = Walker->next;
     1122    // first term: distance to target
     1123    Runner = PermutationMap[Walker->nr];   // find target point
     1124    tmp = (Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)));
     1125    tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     1126    result += constants[0] * tmp;
     1127    //*out << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
     1128   
     1129    // second term: sum of distances to other trajectories
     1130    Runner = start;
     1131    while (Runner->next != end) {
     1132      Runner = Runner->next;
     1133      if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
     1134        break;
     1135      // determine normalized trajectories direction vector (n1, n2)
     1136      Sprinter = PermutationMap[Walker->nr];   // find first target point
     1137      trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));
     1138      trajectory1.SubtractVector(&Trajectories[Walker].R.at(startstep));
     1139      trajectory1.Normalize();
     1140      Norm1 = trajectory1.Norm();
     1141      Sprinter = PermutationMap[Runner->nr];   // find second target point
     1142      trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));
     1143      trajectory2.SubtractVector(&Trajectories[Runner].R.at(startstep));
     1144      trajectory2.Normalize();
     1145      Norm2 = trajectory1.Norm();
     1146      // check whether either is zero()
     1147      if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
     1148        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
     1149      } else if (Norm1 < MYEPSILON) {
     1150        Sprinter = PermutationMap[Walker->nr];   // find first target point
     1151        trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy first offset
     1152        trajectory1.SubtractVector(&Trajectories[Runner].R.at(startstep));  // subtract second offset
     1153        trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything
     1154        trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
     1155        tmp = trajectory1.Norm();  // remaining norm is distance
     1156      } else if (Norm2 < MYEPSILON) {
     1157        Sprinter = PermutationMap[Runner->nr];   // find second target point
     1158        trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy second offset
     1159        trajectory2.SubtractVector(&Trajectories[Walker].R.at(startstep));  // subtract first offset
     1160        trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything
     1161        trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
     1162        tmp = trajectory2.Norm();  // remaining norm is distance
     1163      } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
     1164//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
     1165//        *out << trajectory1;
     1166//        *out << " and ";
     1167//        *out << trajectory2;
     1168        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
     1169//        *out << " with distance " << tmp << "." << endl;
     1170      } else { // determine distance by finding minimum distance
     1171//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
     1172//        *out << endl;
     1173//        *out << "First Trajectory: ";
     1174//        *out << trajectory1 << endl;
     1175//        *out << "Second Trajectory: ";
     1176//        *out << trajectory2 << endl;
     1177        // determine normal vector for both
     1178        normal.MakeNormalVector(&trajectory1, &trajectory2);
     1179        // print all vectors for debugging
     1180//        *out << "Normal vector in between: ";
     1181//        *out << normal << endl;
     1182        // setup matrix
     1183        for (int i=NDIM;i--;) {
     1184          gsl_matrix_set(A, 0, i, trajectory1.x[i]);
     1185          gsl_matrix_set(A, 1, i, trajectory2.x[i]);
     1186          gsl_matrix_set(A, 2, i, normal.x[i]);
     1187          gsl_vector_set(x,i, (Trajectories[Walker].R.at(startstep).x[i] - Trajectories[Runner].R.at(startstep).x[i]));
     1188        }
     1189        // solve the linear system by Householder transformations
     1190        gsl_linalg_HH_svx(A, x);
     1191        // distance from last component
     1192        tmp = gsl_vector_get(x,2);
     1193//        *out << " with distance " << tmp << "." << endl;
     1194        // test whether we really have the intersection (by checking on c_1 and c_2)
     1195        TestVector.CopyVector(&Trajectories[Runner].R.at(startstep));
     1196        trajectory2.Scale(gsl_vector_get(x,1));
     1197        TestVector.AddVector(&trajectory2);
     1198        normal.Scale(gsl_vector_get(x,2));
     1199        TestVector.AddVector(&normal);
     1200        TestVector.SubtractVector(&Trajectories[Walker].R.at(startstep));
     1201        trajectory1.Scale(gsl_vector_get(x,0));
     1202        TestVector.SubtractVector(&trajectory1);
     1203        if (TestVector.Norm() < MYEPSILON) {
     1204//          *out << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl; 
     1205        } else {
     1206//          *out << Verbose(2) << "Test: failed.\tIntersection is off by ";
     1207//          *out << TestVector;
     1208//          *out << "." << endl; 
     1209        }
     1210      }
     1211      // add up
     1212      tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
     1213      if (fabs(tmp) > MYEPSILON) {
     1214        result += constants[1] * 1./tmp;
     1215        //*out << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
     1216      }
     1217    }
     1218     
     1219    // third term: penalty for equal targets
     1220    Runner = start;
     1221    while (Runner->next != end) {
     1222      Runner = Runner->next;
     1223      if ((PermutationMap[Walker->nr] == PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
     1224        Sprinter = PermutationMap[Walker->nr];
     1225//        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
     1226//        *out << Trajectories[Sprinter].R.at(endstep);
     1227//        *out << ", penalting." << endl;
     1228        result += constants[2];
     1229        //*out << Verbose(4) << "Adding " << constants[2] << "." << endl;
     1230      }
     1231    }
     1232  }
     1233 
     1234  return result;
     1235};
     1236
     1237void PrintPermutationMap(ofstream *out, atom **PermutationMap, int Nr)
     1238{
     1239  stringstream zeile1, zeile2;
     1240  int *DoubleList = (int *) Malloc(Nr*sizeof(int), "PrintPermutationMap: *DoubleList");
     1241  int doubles = 0;
     1242  for (int i=0;i<Nr;i++)
     1243    DoubleList[i] = 0;
     1244  zeile1 << "PermutationMap: ";
     1245  zeile2 << "                ";
     1246  for (int i=0;i<Nr;i++) {
     1247    DoubleList[PermutationMap[i]->nr]++;
     1248    zeile1 << i << " ";
     1249    zeile2 << PermutationMap[i]->nr << " ";
     1250  }
     1251  for (int i=0;i<Nr;i++)
     1252    if (DoubleList[i] > 1)
     1253    doubles++;
     1254 // *out << "Found " << doubles << " Doubles." << endl;
     1255  Free((void **)&DoubleList, "PrintPermutationMap: *DoubleList");
     1256//  *out << zeile1.str() << endl << zeile2.str() << endl;
     1257};
     1258
     1259/** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
     1260 * We do the following:
     1261 *  -# Generate a distance list from all source to all target points
     1262 *  -# Sort this per source point
     1263 *  -# Take for each source point the target point with minimum distance, use this as initial permutation
     1264 *  -# check whether molecule::ConstrainedPotential() is greater than injective penalty
     1265 *     -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
     1266 *  -# Next, we only apply transformations that keep the injectivity of the permutations list.
     1267 *  -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target
     1268 *     point and try to change it for one with lesser distance, or for the next one with greater distance, but only
     1269 *     if this decreases the conditional potential.
     1270 *  -# finished.
     1271 *  -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
     1272 *     that the total force is always pointing in direction of the constraint force (ensuring that we move in the
     1273 *     right direction).
     1274 *  -# Finally, we calculate the potential energy and return.
     1275 * \param *out output stream for debugging
     1276 * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
     1277 * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
     1278 * \param endstep step giving final position in constrained MD
     1279 * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
     1280 * \sa molecule::VerletForceIntegration()
     1281 * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
     1282 * \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
     1283 *       to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
     1284 * \bug this all is not O(N log N) but O(N^2)
     1285 */
     1286double molecule::MinimiseConstrainedPotential(ofstream *out, atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem)
     1287{
     1288  double Potential, OldPotential, OlderPotential;
     1289  PermutationMap = (atom **) Malloc(AtomCount*sizeof(atom *), "molecule::MinimiseConstrainedPotential: **PermutationMap");
     1290  DistanceMap **DistanceList = (DistanceMap **) Malloc(AtomCount*sizeof(DistanceMap *), "molecule::MinimiseConstrainedPotential: **DistanceList");
     1291  DistanceMap::iterator *DistanceIterators = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *DistanceIterators");
     1292  int *DoubleList = (int *) Malloc(AtomCount*sizeof(int), "molecule::MinimiseConstrainedPotential: *DoubleList");
     1293  DistanceMap::iterator *StepList = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *StepList");
     1294  double constants[3];
     1295  int round;
     1296  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
     1297  DistanceMap::iterator Rider, Strider;
     1298 
     1299  /// Minimise the potential
     1300  // set Lagrange multiplier constants
     1301  constants[0] = 10.;
     1302  constants[1] = 1.;
     1303  constants[2] = 1e+7;    // just a huge penalty
     1304  // generate the distance list
     1305  *out << Verbose(1) << "Creating the distance list ... " << endl;
     1306  for (int i=AtomCount; i--;) {
     1307    DoubleList[i] = 0;  // stores for how many atoms in startstep this atom is a target in endstep
     1308    DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
     1309    DistanceList[i]->clear();
     1310  }
     1311  *out << Verbose(1) << "Filling the distance list ... " << endl;
     1312  Walker = start;
     1313  while (Walker->next != end) {
     1314    Walker = Walker->next;
     1315    Runner = start;
     1316    while(Runner->next != end) {
     1317      Runner = Runner->next;
     1318      DistanceList[Walker->nr]->insert( DistancePair(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)), Runner) );
     1319    }
     1320  }
     1321  // create the initial PermutationMap (source -> target)
     1322  Walker = start;
     1323  while (Walker->next != end) {
     1324    Walker = Walker->next;
     1325    StepList[Walker->nr] = DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
     1326    PermutationMap[Walker->nr] = DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
     1327    DoubleList[DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
     1328    DistanceIterators[Walker->nr] = DistanceList[Walker->nr]->begin();    // and remember which one we picked
     1329    *out << *Walker << " starts with distance " << DistanceList[Walker->nr]->begin()->first << "." << endl;
     1330  }
     1331  *out << Verbose(1) << "done." << endl;
     1332  // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
     1333  *out << Verbose(1) << "Making the PermutationMap injective ... " << endl;
     1334  Walker = start;
     1335  DistanceMap::iterator NewBase;
     1336  OldPotential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
     1337  while ((OldPotential) > constants[2]) {
     1338    PrintPermutationMap(out, PermutationMap, AtomCount);
     1339    Walker = Walker->next;
     1340    if (Walker == end) // round-robin at the end
     1341      Walker = start->next;
     1342    if (DoubleList[DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
     1343      continue;
     1344    // now, try finding a new one
     1345    NewBase = DistanceIterators[Walker->nr];  // store old base
     1346    do {
     1347      NewBase++;  // take next further distance in distance to targets list that's a target of no one
     1348    } while ((DoubleList[NewBase->second->nr] != 0) && (NewBase != DistanceList[Walker->nr]->end()));
     1349    if (NewBase != DistanceList[Walker->nr]->end()) {
     1350      PermutationMap[Walker->nr] = NewBase->second;
     1351      Potential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
     1352      if (Potential > OldPotential) { // undo
     1353        PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
     1354      } else {  // do
     1355        DoubleList[DistanceIterators[Walker->nr]->second->nr]--;  // decrease the old entry in the doubles list
     1356        DoubleList[NewBase->second->nr]++;    // increase the old entry in the doubles list
     1357        DistanceIterators[Walker->nr] = NewBase;
     1358        OldPotential = Potential;
     1359        *out << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl;
     1360      }
     1361    }
     1362  }
     1363  for (int i=AtomCount; i--;) // now each single entry in the DoubleList should be <=1
     1364    if (DoubleList[i] > 1) {
     1365      cerr << "Failed to create an injective PermutationMap!" << endl;
     1366      exit(1);
     1367    }
     1368  *out << Verbose(1) << "done." << endl;
     1369  Free((void **)&DoubleList, "molecule::MinimiseConstrainedPotential: *DoubleList");
     1370  // argument minimise the constrained potential in this injective PermutationMap
     1371  *out << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
     1372  OldPotential = 1e+10;
     1373  round = 0;
     1374  do {
     1375    *out << "Starting round " << ++round << " ... " << endl;
     1376    OlderPotential = OldPotential;
     1377    do {
     1378      Walker = start;
     1379      while (Walker->next != end) { // pick one
     1380        Walker = Walker->next;
     1381        PrintPermutationMap(out, PermutationMap, AtomCount);
     1382        Sprinter = DistanceIterators[Walker->nr]->second;   // store initial partner
     1383        Strider = DistanceIterators[Walker->nr];  //remember old iterator
     1384        DistanceIterators[Walker->nr] = StepList[Walker->nr]; 
     1385        if (DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
     1386          DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->begin();
     1387          break;
     1388        }
     1389        //*out << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
     1390        // find source of the new target
     1391        Runner = start->next;
     1392        while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
     1393          if (PermutationMap[Runner->nr] == DistanceIterators[Walker->nr]->second) {
     1394            //*out << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
     1395            break;
     1396          }
     1397          Runner = Runner->next;
     1398        }
     1399        if (Runner != end) { // we found the other source
     1400          // then look in its distance list for Sprinter
     1401          Rider = DistanceList[Runner->nr]->begin();
     1402          for (; Rider != DistanceList[Runner->nr]->end(); Rider++)
     1403            if (Rider->second == Sprinter)
     1404              break;
     1405          if (Rider != DistanceList[Runner->nr]->end()) { // if we have found one
     1406            //*out << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl;
     1407            // exchange both
     1408            PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap
     1409            PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
     1410            PrintPermutationMap(out, PermutationMap, AtomCount);
     1411            // calculate the new potential
     1412            //*out << Verbose(2) << "Checking new potential ..." << endl;
     1413            Potential = ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
     1414            if (Potential > OldPotential) { // we made everything worse! Undo ...
     1415              //*out << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
     1416              //*out << Verbose(3) << "Setting " << *Runner << "'s source to " << *DistanceIterators[Runner->nr]->second << "." << endl;
     1417              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
     1418              PermutationMap[Runner->nr] = DistanceIterators[Runner->nr]->second;
     1419              // Undo for Walker
     1420              DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
     1421              //*out << Verbose(3) << "Setting " << *Walker << "'s source to " << *DistanceIterators[Walker->nr]->second << "." << endl;
     1422              PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
     1423            } else {
     1424              DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
     1425              *out << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl;
     1426              OldPotential = Potential;
     1427            }
     1428            if (Potential > constants[2]) {
     1429              cerr << "ERROR: The two-step permutation procedure did not maintain injectivity!" << endl;
     1430              exit(255);
     1431            }
     1432            //*out << endl;
     1433          } else {
     1434            cerr << "ERROR: " << *Runner << " was not the owner of " << *Sprinter << "!" << endl;
     1435            exit(255);
     1436          }
     1437        } else {
     1438          PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // new target has no source!
     1439        }
     1440        StepList[Walker->nr]++; // take next farther distance target
     1441      }
     1442    } while (Walker->next != end);
     1443  } while ((OlderPotential - OldPotential) > 1e-3);
     1444  *out << Verbose(1) << "done." << endl;
     1445
     1446 
     1447  /// free memory and return with evaluated potential
     1448  for (int i=AtomCount; i--;)
     1449    DistanceList[i]->clear();
     1450  Free((void **)&DistanceList, "molecule::MinimiseConstrainedPotential: **DistanceList");
     1451  Free((void **)&DistanceIterators, "molecule::MinimiseConstrainedPotential: *DistanceIterators");
     1452  return ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
     1453};
     1454
     1455/** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
     1456 * \param *out output stream for debugging
     1457 * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
     1458 * \param endstep step giving final position in constrained MD
     1459 * \param **PermutationMap mapping between the atom label of the initial and the final configuration
     1460 * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
     1461 * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
     1462 */
     1463void molecule::EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
     1464{
     1465  double constant = 10.;
     1466  atom *Walker = NULL, *Sprinter = NULL;
     1467 
     1468  /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
     1469  *out << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl;
     1470  Walker = start;
     1471  while (Walker->next != NULL) {
     1472    Walker = Walker->next;
     1473    Sprinter = PermutationMap[Walker->nr];
     1474    // set forces
     1475    for (int i=NDIM;i++;)
     1476      Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
     1477  } 
     1478  *out << Verbose(1) << "done." << endl;
     1479};
     1480
     1481/** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
     1482 * Note, step number is config::MaxOuterStep
     1483 * \param *out output stream for debugging
     1484 * \param startstep stating initial configuration in molecule::Trajectories
     1485 * \param endstep stating final configuration in molecule::Trajectories
     1486 * \param &config configuration structure
     1487 * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
     1488 */
     1489bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration)
     1490{
     1491  molecule *mol = NULL;
     1492  bool status = true;
     1493  int MaxSteps = configuration.MaxOuterStep;
     1494  MoleculeListClass *MoleculePerStep = new MoleculeListClass();
     1495  // Get the Permutation Map by MinimiseConstrainedPotential
     1496  atom **PermutationMap = NULL;
     1497  atom *Walker = NULL, *Sprinter = NULL;
     1498  MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
     1499
     1500  // check whether we have sufficient space in Trajectories for each atom
     1501  Walker = start;
     1502  while (Walker->next != end) {
     1503    Walker = Walker->next;
     1504    if (Trajectories[Walker].R.size() <= (unsigned int)(MaxSteps)) {
     1505      //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
     1506      Trajectories[Walker].R.resize(MaxSteps+1);
     1507      Trajectories[Walker].U.resize(MaxSteps+1);
     1508      Trajectories[Walker].F.resize(MaxSteps+1);
     1509    }
     1510  }
     1511  // push endstep to last one
     1512  Walker = start;
     1513  while (Walker->next != end) { // remove the endstep (is now the very last one)
     1514    Walker = Walker->next;
     1515    for (int n=NDIM;n--;) {
     1516      Trajectories[Walker].R.at(MaxSteps).x[n] = Trajectories[Walker].R.at(endstep).x[n];
     1517      Trajectories[Walker].U.at(MaxSteps).x[n] = Trajectories[Walker].U.at(endstep).x[n];
     1518      Trajectories[Walker].F.at(MaxSteps).x[n] = Trajectories[Walker].F.at(endstep).x[n];
     1519    }
     1520  }
     1521  endstep = MaxSteps;
     1522 
     1523  // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
     1524  *out << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl;
     1525  for (int step = 0; step <= MaxSteps; step++) {
     1526    mol = new molecule(elemente);
     1527    MoleculePerStep->insert(mol);
     1528    Walker = start;
     1529    while (Walker->next != end) {
     1530      Walker = Walker->next;
     1531      // add to molecule list
     1532      Sprinter = mol->AddCopyAtom(Walker);
     1533      for (int n=NDIM;n--;) {
     1534        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);
     1535        // add to Trajectories
     1536        //*out << Verbose(3) << step << ">=" << MDSteps-1 << endl;
     1537        if (step < MaxSteps) {
     1538          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);
     1539          Trajectories[Walker].U.at(step).x[n] = 0.;
     1540          Trajectories[Walker].F.at(step).x[n] = 0.;
     1541        }
     1542      }
     1543    }
     1544  }
     1545  MDSteps = MaxSteps+1;   // otherwise new Trajectories' points aren't stored on save&exit
     1546 
     1547  // store the list to single step files
     1548  int *SortIndex = (int *) Malloc(AtomCount*sizeof(int), "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
     1549  for (int i=AtomCount; i--; )
     1550    SortIndex[i] = i;
     1551  status = MoleculePerStep->OutputConfigForListOfFragments(out, &configuration, SortIndex);
     1552 
     1553  // free and return
     1554  Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
     1555  delete(MoleculePerStep);
     1556  return status;
     1557};
     1558
    10851559/** Parses nuclear forces from file and performs Verlet integration.
    10861560 * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
    10871561 * have to transform them).
    10881562 * This adds a new MD step to the config file.
     1563 * \param *out output stream for debugging
    10891564 * \param *file filename
     1565 * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
    10901566 * \param delta_t time step width in atomic units
    10911567 * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
     1568 * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential()
    10921569 * \return true - file found and parsed, false - file not found or imparsable
    1093  */
    1094 bool molecule::VerletForceIntegration(char *file, double delta_t, bool IsAngstroem)
    1095 {
    1096   element *runner = elemente->start;
     1570 * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
     1571 */
     1572bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
     1573{
    10971574  atom *walker = NULL;
    1098   int AtomNo;
    10991575  ifstream input(file);
    11001576  string token;
    11011577  stringstream item;
    1102   double a, IonMass;
     1578  double IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp;
    11031579  ForceMatrix Force;
    1104   Vector tmpvector;
    11051580
    11061581  CountElements();  // make sure ElementsInMolecule is up to date
    … …  
    11201595    }
    11211596    // correct Forces
    1122 //    for(int d=0;d<NDIM;d++)
    1123 //      tmpvector.x[d] = 0.;
    1124 //    for(int i=0;i<AtomCount;i++)
    1125 //      for(int d=0;d<NDIM;d++) {
    1126 //        tmpvector.x[d] += Force.Matrix[0][i][d+5];
    1127 //      }
    1128 //    for(int i=0;i<AtomCount;i++)
    1129 //      for(int d=0;d<NDIM;d++) {
    1130 //        Force.Matrix[0][i][d+5] -= tmpvector.x[d]/(double)AtomCount;
    1131 //      }
     1597    for(int d=0;d<NDIM;d++)
     1598      Vector[d] = 0.;
     1599    for(int i=0;i<AtomCount;i++)
     1600      for(int d=0;d<NDIM;d++) {
     1601        Vector[d] += Force.Matrix[0][i][d+5];
     1602      }
     1603    for(int i=0;i<AtomCount;i++)
     1604      for(int d=0;d<NDIM;d++) {
     1605        Force.Matrix[0][i][d+5] -= Vector[d]/(double)AtomCount;
     1606      }
     1607    // solve a constrained potential if we are meant to
     1608    if (configuration.DoConstrainedMD) {
     1609      // calculate forces and potential
     1610      atom **PermutationMap = NULL;
     1611      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
     1612      EvaluateConstrainedForces(out, configuration.DoConstrainedMD, 0, PermutationMap, &Force);
     1613      Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
     1614    }
     1615   
    11321616    // and perform Verlet integration for each atom with position, velocity and force vector
    1133     runner = elemente->start;
    1134     while (runner->next != elemente->end) { // go through every element
    1135       runner = runner->next;
    1136       IonMass = runner->mass;
    1137       a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    1138       if (ElementsInMolecule[runner->Z]) { // if this element got atoms
    1139         AtomNo = 0;
     1617    walker = start;
     1618    while (walker->next != end) { // go through every atom of this element
     1619      walker = walker->next;
     1620      //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
     1621      // check size of vectors
     1622      if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
     1623        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
     1624        Trajectories[walker].R.resize(MDSteps+10);
     1625        Trajectories[walker].U.resize(MDSteps+10);
     1626        Trajectories[walker].F.resize(MDSteps+10);
     1627      }
     1628     
     1629      // Update R (and F)
     1630      for (int d=0; d<NDIM; d++) {
     1631        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
     1632        Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
     1633        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
     1634        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
     1635      }
     1636      // Update U
     1637      for (int d=0; d<NDIM; d++) {
     1638        Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
     1639        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
     1640      }
     1641//      out << "Integrated position&velocity of step " << (MDSteps) << ": (";
     1642//      for (int d=0;d<NDIM;d++)
     1643//        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
     1644//      out << ")\t(";
     1645//      for (int d=0;d<NDIM;d++)
     1646//        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
     1647//      out << ")" << endl;
     1648            // next atom
     1649    }
     1650  }
     1651  // correct velocities (rather momenta) so that center of mass remains motionless
     1652  for(int d=0;d<NDIM;d++)
     1653    Vector[d] = 0.;
     1654  IonMass = 0.;
     1655  walker = start;
     1656  while (walker->next != end) { // go through every atom
     1657    walker = walker->next;
     1658    IonMass += walker->type->mass;  // sum up total mass
     1659    for(int d=0;d<NDIM;d++) {
     1660      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
     1661    }
     1662  }
     1663  // correct velocities (rather momenta) so that center of mass remains motionless
     1664  for(int d=0;d<NDIM;d++)
     1665    Vector[d] /= IonMass;
     1666  ActualTemp = 0.;
     1667  walker = start;
     1668  while (walker->next != end) { // go through every atom of this element
     1669    walker = walker->next;
     1670    for(int d=0;d<NDIM;d++) {
     1671      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d];
     1672      ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
     1673    }
     1674  }
     1675  Thermostats(configuration, ActualTemp, Berendsen);
     1676  MDSteps++;
     1677
     1678
     1679  // exit
     1680  return true;
     1681};
     1682
     1683/** Implementation of various thermostats.
     1684 * All these thermostats apply an additional force which has the following forms:
     1685 * -# Woodcock
     1686 *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
     1687 * -# Gaussian
     1688 *  \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$
     1689 * -# Langevin
     1690 *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$ 
     1691 * -# Berendsen
     1692 *  \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$
     1693 * -# Nose-Hoover
     1694 *  \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$
     1695 * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
     1696 * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
     1697 * belatedly and the constraint force set.
     1698 * \param *P Problem at hand
     1699 * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
     1700 * \sa InitThermostat()
     1701 */
     1702void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
     1703{
     1704  double ekin = 0.;
     1705  double E = 0., G = 0.;
     1706  double delta_alpha = 0.;
     1707  double ScaleTempFactor;
     1708  double sigma;
     1709  double IonMass;
     1710  int d;
     1711  gsl_rng * r;
     1712  const gsl_rng_type * T;
     1713  double *U = NULL, *F = NULL, FConstraint[NDIM];
     1714  atom *walker = NULL;
     1715 
     1716  // calculate scale configuration
     1717  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
     1718   
     1719  // differentating between the various thermostats
     1720  switch(Thermostat) {
     1721     case None:
     1722      cout << Verbose(2) <<  "Applying no thermostat..." << endl;
     1723      break;
     1724     case Woodcock:
     1725      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
     1726        cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
    11401727        walker = start;
    11411728        while (walker->next != end) { // go through every atom of this element
    11421729          walker = walker->next;
    1143           if (walker->type == runner) { // if this atom fits to element
    1144             // check size of vectors
    1145             if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
    1146               //cout << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
    1147               Trajectories[walker].R.resize(MDSteps+10);
    1148               Trajectories[walker].U.resize(MDSteps+10);
    1149               Trajectories[walker].F.resize(MDSteps+10);
     1730          IonMass = walker->type->mass;
     1731          U = Trajectories[walker].U.at(MDSteps).x;
     1732          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1733            for (d=0; d<NDIM; d++) {
     1734              U[d] *= sqrt(ScaleTempFactor);
     1735              ekin += 0.5*IonMass * U[d]*U[d];
    11501736            }
    1151             // 1. calculate x(t+\delta t)
    1152             for (int d=0; d<NDIM; d++) {
    1153               Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5];
    1154               Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
    1155               Trajectories[walker].R.at(MDSteps).x[d] += delta_t*(Trajectories[walker].U.at(MDSteps-1).x[d]);
    1156               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
     1737        }
     1738      }
     1739      break;
     1740     case Gaussian:
     1741      cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
     1742      walker = start;
     1743      while (walker->next != end) { // go through every atom of this element
     1744        walker = walker->next;
     1745        IonMass = walker->type->mass;
     1746        U = Trajectories[walker].U.at(MDSteps).x;
     1747        F = Trajectories[walker].F.at(MDSteps).x;
     1748        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1749          for (d=0; d<NDIM; d++) {
     1750            G += U[d] * F[d];
     1751            E += U[d]*U[d]*IonMass;
     1752          }
     1753      }
     1754      cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
     1755      walker = start;
     1756      while (walker->next != end) { // go through every atom of this element
     1757        walker = walker->next;
     1758        IonMass = walker->type->mass;
     1759        U = Trajectories[walker].U.at(MDSteps).x;
     1760        F = Trajectories[walker].F.at(MDSteps).x;
     1761        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
     1762          for (d=0; d<NDIM; d++) {
     1763            FConstraint[d] = (G/E) * (U[d]*IonMass);
     1764            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
     1765            ekin += IonMass * U[d]*U[d];
     1766          }
     1767      }
     1768      break;
     1769     case Langevin:
     1770      cout << Verbose(2) <<  "Applying Langevin thermostat..." << endl;
     1771      // init random number generator
     1772      gsl_rng_env_setup();
     1773      T = gsl_rng_default;
     1774      r = gsl_rng_alloc (T);
     1775      // Go through each ion
     1776      walker = start;
     1777      while (walker->next != end) { // go through every atom of this element
     1778        walker = walker->next;
     1779        IonMass = walker->type->mass;
     1780        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
     1781        U = Trajectories[walker].U.at(MDSteps).x;
     1782        F = Trajectories[walker].F.at(MDSteps).x;
     1783        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1784          // throw a dice to determine whether it gets hit by a heat bath particle
     1785          if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) { 
     1786            cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
     1787            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
     1788            for (d=0; d<NDIM; d++) {
     1789              U[d] = gsl_ran_gaussian (r, sigma);
    11571790            }
    1158             // 2. Calculate v(t+\delta t)
    1159             for (int d=0; d<NDIM; d++) {
    1160               Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d];
    1161               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;
    1162             }
    1163 //            cout << "Integrated position&velocity of step " << (MDSteps) << ": (";
    1164 //            for (int d=0;d<NDIM;d++)
    1165 //              cout << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
    1166 //            cout << ")\t(";
    1167 //            for (int d=0;d<NDIM;d++)
    1168 //              cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
    1169 //            cout << ")" << endl;
    1170             // next atom
    1171             AtomNo++;
     1791            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
     1792          }
     1793          for (d=0; d<NDIM; d++)
     1794            ekin += 0.5*IonMass * U[d]*U[d];
     1795        }
     1796      }
     1797      break;
     1798     case Berendsen:
     1799      cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
     1800      walker = start;
     1801      while (walker->next != end) { // go through every atom of this element
     1802        walker = walker->next;
     1803        IonMass = walker->type->mass;
     1804        U = Trajectories[walker].U.at(MDSteps).x;
     1805        F = Trajectories[walker].F.at(MDSteps).x;
     1806        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1807          for (d=0; d<NDIM; d++) {
     1808            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
     1809            ekin += 0.5*IonMass * U[d]*U[d];
    11721810          }
    11731811        }
    11741812      }
    1175     }
    1176   }
    1177 //  // correct velocities (rather momenta) so that center of mass remains motionless
    1178 //  tmpvector.zero()
    1179 //  IonMass = 0.;
    1180 //  walker = start;
    1181 //  while (walker->next != end) { // go through every atom
    1182 //    walker = walker->next;
    1183 //    IonMass += walker->type->mass;  // sum up total mass
    1184 //    for(int d=0;d<NDIM;d++) {
    1185 //      tmpvector.x[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
    1186 //    }
    1187 //  }
    1188 //  walker = start;
    1189 //  while (walker->next != end) { // go through every atom of this element
    1190 //    walker = walker->next;
    1191 //    for(int d=0;d<NDIM;d++) {
    1192 //      Trajectories[walker].U.at(MDSteps).x[d] -= tmpvector.x[d]*walker->type->mass/IonMass;
    1193 //    }
    1194 //  }
    1195   MDSteps++;
    1196 
    1197 
    1198   // exit
    1199   return true;
    1200 };
    1201 
    1202 /** Align all atoms in such a manner that given vector \a *n is along z axis.
     1813      break;
     1814     case NoseHoover:
     1815      cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
     1816      // dynamically evolve alpha (the additional degree of freedom)
     1817      delta_alpha = 0.;
     1818      walker = start;
     1819      while (walker->next != end) { // go through every atom of this element
     1820        walker = walker->next;
     1821        IonMass = walker->type->mass;
     1822        U = Trajectories[walker].U.at(MDSteps).x;
     1823        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1824          for (d=0; d<NDIM; d++) {
     1825            delta_alpha += U[d]*U[d]*IonMass;
     1826          }
     1827        }
     1828      }
     1829      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
     1830      configuration.alpha += delta_alpha*configuration.Deltat;
     1831      cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
     1832      // apply updated alpha as additional force
     1833      walker = start;
     1834      while (walker->next != end) { // go through every atom of this element
     1835        walker = walker->next;
     1836        IonMass = walker->type->mass;
     1837        U = Trajectories[walker].U.at(MDSteps).x;
     1838        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
     1839          for (d=0; d<NDIM; d++) {
     1840              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
     1841              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
     1842              ekin += (0.5*IonMass) * U[d]*U[d];
     1843            }
     1844        }
     1845      }
     1846      break;
     1847  }   
     1848  cout << Verbose(1) << "Kinetic energy is " << ekin << "." << endl;
     1849};
     1850
     1851/** Align all atoms in such a manner that given vector \a *n is along z axis.
    12031852 * \param n[] alignment vector.
    12041853 */
    … …  
    12671916bool molecule::RemoveAtom(atom *pointer)
    12681917{
    1269   if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
     1918  if (ElementsInMolecule[pointer->type->Z] != 0)  { // this would indicate an error
    12701919    ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
    1271   else
     1920    AtomCount--;
     1921  } else
    12721922    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
    12731923  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
    … …  
    30933743  while (MolecularWalker->next != NULL) {
    30943744    MolecularWalker = MolecularWalker->next;
     3745    *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
    30953746    LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
    30963747//    // check the list of local atoms for debugging
    … …  
    31083759    delete(LocalBackEdgeStack);
    31093760  }
    3110 
     3761 
    31113762  // ===== 3. if structure still valid, parse key set file and others =====
    31123763  FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
    … …  
    31143765  // ===== 4. check globally whether there's something to do actually (first adaptivity check)
    31153766  FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
    3116 
    3117   // =================================== Begin of FRAGMENTATION ===============================
    3118   // ===== 6a. assign each keyset to its respective subgraph =====
     3767 
     3768  // =================================== Begin of FRAGMENTATION =============================== 
     3769  // ===== 6a. assign each keyset to its respective subgraph ===== 
    31193770  Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);
    31203771
    … …  
    31513802  delete(ParsedFragmentList);
    31523803  delete[](MinimumRingSize);
    3153 
     3804 
    31543805
    31553806  // ==================================== End of FRAGMENTATION ============================================
    … …  
    32443895  atom *Walker = NULL, *OtherAtom = NULL;
    32453896  ReferenceStack->Push(Binder);
    3246 
     3897 
    32473898  do {  // go through all bonds and push local ones
    32483899    Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
    … …  
    37144365};
    37154366
     4367/** Creates \a MoleculeListClass of all unique fragments of the \a molecule containing \a Order atoms or vertices.
     4368 * The picture to have in mind is that of a DFS "snake" of a certain length \a Order, i.e. as in the infamous
     4369 * computer game, that winds through the connected graph representing the molecule. Color (white,
     4370 * lightgray, darkgray, black) indicates whether a vertex has been discovered so far or not. Labels will help in
     4371 * creating only unique fragments and not additional ones with vertices simply in different sequence.
     4372 * The Predecessor is always the one that came before in discovering, needed on backstepping. And
     4373 * finally, the ShortestPath is needed for removing vertices from the snake stack during the back-
     4374 * stepping.
     4375 * \param *out output stream for debugging
     4376 * \param Order number of atoms in each fragment
     4377 * \param *configuration configuration for writing config files for each fragment
     4378 * \return List of all unique fragments with \a Order atoms
     4379 */
     4380/*
     4381MoleculeListClass * molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, config *configuration)
     4382{
     4383  atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
     4384  int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
     4385  int *Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
     4386  enum Shading *ColorVertexList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
     4387  enum Shading *ColorEdgeList = (enum Shading *) Malloc(sizeof(enum Shading)*BondCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorBondList");
     4388  StackClass<atom *> *RootStack = new StackClass<atom *>(AtomCount);
     4389  StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
     4390  StackClass<atom *> *SnakeStack = new StackClass<atom *>(Order+1); // equal to Order is not possible, as then the StackClass<atom *> cannot discern between full and empty stack!
     4391  MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
     4392  MoleculeListClass *FragmentList = NULL;
     4393  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
     4394  bond *Binder = NULL;
     4395  int RunningIndex = 0, FragmentCounter = 0;
     4396
     4397  *out << Verbose(1) << "Begin of CreateListOfUniqueFragmentsOfOrder." << endl;
     4398
     4399  // reset parent list
     4400  *out << Verbose(3) << "Resetting labels, parent, predecessor, color and shortest path lists." << endl;
     4401  for (int i=0;i<AtomCount;i++) { // reset all atom labels
     4402    // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
     4403    Labels[i] = -1;
     4404    SonList[i] = NULL;
     4405    PredecessorList[i] = NULL;
     4406    ColorVertexList[i] = white;
     4407    ShortestPathList[i] = -1;
     4408  }
     4409  for (int i=0;i<BondCount;i++)
     4410    ColorEdgeList[i] = white;
     4411  RootStack->ClearStack();  // clearstack and push first atom if exists
     4412  TouchedStack->ClearStack();
     4413  Walker = start->next;
     4414  while ((Walker != end)
     4415#ifdef ADDHYDROGEN
     4416   && (Walker->type->Z == 1)
     4417        }
     4418      }
     4419      *out << ", SP of " << ShortestPathList[Walker->nr]  << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
     4420
     4421      // then check the stack for a newly stumbled upon fragment
     4422      if (SnakeStack->ItemCount() == Order) { // is stack full?
     4423        // store the fragment if it is one and get a removal candidate
     4424        Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
     4425        // remove the candidate if one was found
     4426        if (Removal != NULL) {
     4427          *out << Verbose(2) << "Removing item " << Removal->Name << " with SP of " << ShortestPathList[Removal->nr] << " from snake stack." << endl;
     4428          SnakeStack->RemoveItem(Removal);
     4429          ColorVertexList[Removal->nr] = lightgray; // return back to not on snake stack but explored marking
     4430          if (Walker == Removal) { // if the current atom is to be removed, we also have to take a step back
     4431            Walker = PredecessorList[Removal->nr];
     4432            *out << Verbose(2) << "Stepping back to " << Walker->Name << "." << endl;
     4433          }
     4434        }
     4435      } else
     4436        Removal = NULL;
     4437
     4438      // finally, look for a white neighbour as the next Walker
     4439      Binder = NULL;
     4440      if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) {  // don't look, if a new walker has been set above
     4441        *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
     4442        OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
     4443        if (ShortestPathList[Walker->nr] < Order) {
     4444          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
     4445            Binder = ListOfBondsPerAtom[Walker->nr][i];
     4446            *out << Verbose(2) << "Current bond is " << *Binder << ": ";
     4447            OtherAtom = Binder->GetOtherAtom(Walker);
     4448            if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
     4449              *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
     4450              //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
     4451            } else { // otherwise check its colour and element
     4452              if (
     4453              (OtherAtom->type->Z != 1) &&
     4454#endif
     4455                    (ColorEdgeList[Binder->nr] == white)) {  // skip hydrogen, look for unexplored vertices
     4456                *out << "Moving along " << GetColor(ColorEdgeList[Binder->nr]) << " bond " << Binder << " to " << ((ColorVertexList[OtherAtom->nr] == white) ? "unexplored" : "explored") << " item: " << OtherAtom->Name << "." << endl;
     4457                // i find it currently rather sensible to always set the predecessor in order to find one's way back
     4458                //if (PredecessorList[OtherAtom->nr] == NULL) {
     4459                PredecessorList[OtherAtom->nr] = Walker;
     4460                *out << Verbose(3) << "Setting Predecessor of " << OtherAtom->Name << " to " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
     4461                //} else {
     4462                //  *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
     4463                //}
     4464                Walker = OtherAtom;
     4465                break;
     4466              } else {
     4467                if (OtherAtom->type->Z == 1)
     4468                  *out << "Links to a hydrogen atom." << endl;
     4469                else
     4470                  *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
     4471              }
     4472            }
     4473          }
     4474        } else {  // means we have stepped beyond the horizon: Return!
     4475          Walker = PredecessorList[Walker->nr];
     4476          OtherAtom = Walker;
     4477          *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
     4478        }
     4479        if (Walker != OtherAtom) {  // if no white neighbours anymore, color it black
     4480          *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
     4481          ColorVertexList[Walker->nr] = black;
     4482          Walker = PredecessorList[Walker->nr];
     4483        }
     4484      }
     4485    } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
     4486    *out << Verbose(2) << "Inner Looping is finished." << endl;
     4487
     4488    // if we reset all AtomCount atoms, we have again technically O(N^2) ...
     4489    *out << Verbose(2) << "Resetting lists." << endl;
     4490    Walker = NULL;
     4491    Binder = NULL;
     4492    while (!TouchedStack->IsEmpty()) {
     4493      Walker = TouchedStack->PopLast();
     4494      *out << Verbose(3) << "Re-initialising entries of " << *Walker << "." << endl;
     4495      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
     4496        ColorEdgeList[ListOfBondsPerAtom[Walker->nr][i]->nr] = white;
     4497      PredecessorList[Walker->nr] = NULL;
     4498      ColorVertexList[Walker->nr] = white;
     4499      ShortestPathList[Walker->nr] = -1;
     4500    }
     4501  }
     4502  *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
     4503
     4504  // copy together
     4505  *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
     4506  FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
     4507  RunningIndex = 0;
     4508  while ((Leaflet != NULL) && (RunningIndex < FragmentCounter))  {
     4509    FragmentList->ListOfMolecules[RunningIndex++] = Leaflet->Leaf;
     4510    Leaflet->Leaf = NULL; // prevent molecule from being removed
     4511    TempLeaf = Leaflet;
     4512    Leaflet = Leaflet->previous;
     4513    delete(TempLeaf);
     4514  };
     4515
     4516  // free memory and exit
     4517  Free((void **)&PredecessorList, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
     4518  Free((void **)&ShortestPathList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
     4519  Free((void **)&Labels, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
     4520  Free((void **)&ColorVertexList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
     4521  delete(RootStack);
     4522  delete(TouchedStack);
     4523  delete(SnakeStack);
     4524
     4525  *out << Verbose(1) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
     4526  return FragmentList;
     4527};
     4528*/
     4529
    37164530/** Structure containing all values in power set combination generation.
    37174531 */
    … …  
    45455359  if (result) {
    45465360    *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
    4547     DetermineCenter(CenterOfGravity);
    4548     OtherMolecule->DetermineCenter(OtherCenterOfGravity);
     5361    DeterminePeriodicCenter(CenterOfGravity);
     5362    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
    45495363    *out << Verbose(5) << "Center of Gravity: ";
    45505364    CenterOfGravity.Output(out);
Note: See TracChangeset for help on using the changeset viewer.