Changes in / [b38b64:631dcb]


Ignore:
Location:
src
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • src/Makefile.am

    rb38b64 r631dcb  
    1010
    1111
    12 #EXTRA_DIST = ${molecuilder_DATA}
     12EXTRA_DIST = ${molecuilder_DATA}

  • src/analyzer.cpp

    rb38b64 r631dcb  
    5858  bool NoHessian = false;
    5959  bool NoTime = false;
    60   double norm;
    6160  int counter;
    6261 
    … …  
    215214    output << endl << "Total Chis" << endl << "===============" << endl << Chi.Header[Chi.MatrixCounter] << endl;
    216215    for(int j=0;j<Chi.RowCounter[Chi.MatrixCounter];j++) {
    217       for(int k=0;k<Chi.ColumnCounter;k++)
     216      for(int k=0;k<Chi.ColumnCounter[Chi.MatrixCounter];k++)
    218217        output << scientific << Chi.Matrix[ Chi.MatrixCounter ][j][k] << "\t";
    219218      output << endl;
    … …  
    223222    output << endl << "Total Chis PAS" << endl << "===============" << endl << ChiPAS.Header[ChiPAS.MatrixCounter] << endl;
    224223    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
    225       for(int k=0;k<ChiPAS.ColumnCounter;k++)
     224      for(int k=0;k<ChiPAS.ColumnCounter[Chi.MatrixCounter];k++)
    226225        output << scientific << ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t";
    227226      output << endl;
    … …  
    306305
    307306  // +++++++++++++++++++++++++++++++++++++++ Plotting shieldings
    308 
    309307  if (periode != NULL) { // also look for PAS values
    310308    if (!CreateDataDeltaForcesOrderPerAtom(ShieldingPAS, ShieldingPASFragments, KeySet, argv[3], "DeltaShieldingsPAS-Order", "Plot of error between approximated shieldings and full shieldings versus the Bond Order", datum)) return 1;
    … …  
    325323    for(int j=0;j<ChiPAS.RowCounter[ChiPAS.MatrixCounter];j++) {
    326324      output << j << "\t";
    327       for(int k=0;k<ChiPAS.ColumnCounter;k++)
     325      for(int k=0;k<ChiPAS.ColumnCounter[ChiPAS.MatrixCounter];k++)
    328326        output << scientific <<  ChiPAS.Matrix[ ChiPAS.MatrixCounter ][j][k] << "\t"; //*(((k>1) && (k<6))? 1.e6 : 1.) << "\t";
    329327      output << endl;
    … …  
    507505    }
    508506    output << "'ShieldingsPAS-Order.dat' index " << KeySet.Order << " title 'Full' using ($1+" << step*(double)KeySet.Order << "):7 with boxes" << endl;
    509     output.close(); 
    510507    output2.close(); 
    511508

  • src/atom.cpp

    rb38b64 r631dcb  
    8585};
    8686
    87 ostream & operator << (ostream &ost, atom &a)
     87ostream & operator << (ostream &ost, const atom &a)
    8888{
    8989  ost << "[" << a.Name << "|" << &a << "]";

  • src/bond.cpp

    rb38b64 r631dcb  
    8282};
    8383
    84 ostream & operator << (ostream &ost, bond &b)
     84ostream & operator << (ostream &ost, const bond &b)
    8585{
    8686  ost << "[" << b.leftatom->Name << " <" << b.BondDegree << "(H" << b.HydrogenBond << ")>" << b.rightatom->Name << "]";
    … …  
    9898  if(rightatom == Atom)
    9999    return leftatom;
     100  cerr << "Bond " << *this << " does not contain atom " << *Atom << "!" << endl;
    100101  return NULL;
    101102};

  • src/boundary.cpp

    rb38b64 r631dcb  
    14691469
    14701470  cout << Verbose(1) << "Adding triangle to its lines" << endl;
    1471   int i = 0;
    14721471  TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
    14731472  TrianglesOnBoundaryCount++;
    … …  
    14761475   * this is apparently done when constructing triangle
    14771476
    1478    for (i=0; i<3; i++)
     1477   for (int i=0; i<3; i++)
    14791478   {
    14801479   BLS[i]->AddTriangle(BTS);
    … …  
    16021601        Vector Mittelpunkt;
    16031602
    1604         double CurrentEpsilon = 0.1;
    1605         double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
     1603        double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius;
    16061604        double BallAngle, AlternativeSign;
    16071605        atom *Walker; // variable atom point
    … …  
    17791777
    17801778        double CurrentEpsilon = 0.1;
    1781         double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius, Distance;
     1779        double alpha, beta, gamma, SideA, SideB, SideC, sign, Umkreisradius, Restradius;
    17821780        double BallAngle;
    17831781        atom *Walker; // variable atom point

  • src/builder.cpp

    rb38b64 r631dcb  
    822822            cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
    823823            cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
     824            cout << "\t-L <step0> <step1> <prefix>\tStore a linear interpolation between two configurations <step0> and <step1> into single config files with prefix <prefix> and as Trajectories into the current config file." << endl;
    824825            cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
    825826            cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
    … …  
    10651066              }
    10661067              break;
     1068            case 'L':
     1069              ExitFlag = 1;
     1070              SaveFlag = true;
     1071              cout << Verbose(1) << "Linear interpolation between configuration " << argv[argptr] << " and " << argv[argptr+1] << "." << endl;
     1072              if (!mol->LinearInterpolationBetweenConfiguration((ofstream *)&cout, atoi(argv[argptr]), atoi(argv[argptr+1]), argv[argptr+2], configuration))
     1073                cout << Verbose(2) << "Could not store " << argv[argptr+2] << " files." << endl;
     1074              else
     1075                cout << Verbose(2) << "Steps created and " << argv[argptr+2] << " files stored." << endl;
     1076              argptr+=3;
     1077              break;
    10671078            case 'P':
    10681079              ExitFlag = 1;
    … …  
    10731084                SaveFlag = true;
    10741085                cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
    1075                 if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
     1086                if (!mol->VerletForceIntegration((ofstream *)&cout, argv[argptr], configuration))
    10761087                  cout << Verbose(2) << "File not found." << endl;
    10771088                else

  • src/config.cpp

    rb38b64 r631dcb  
    1818  configpath = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
    1919  configname = (char *) MallocString(sizeof(char)*MAXSTRINGSIZE,"config constructor: mainname");
     20  ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "IonsInitRead: *ThermostatImplemented");
     21  ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "IonsInitRead: *ThermostatNames");
     22  for (int j=0;j<MaxThermostats;j++)
     23    ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "IonsInitRead: ThermostatNames[]");
     24  Thermostat = 4;
     25  alpha = 0.;
     26  ScaleTempStep = 25;
     27  TempFrequency = 2.5;
    2028  strcpy(mainname,"pcp");
    2129  strcpy(defaultpath,"not specified");
    … …  
    2533  basis="3-21G";
    2634 
     35
     36  strcpy(ThermostatNames[0],"None");
     37  ThermostatImplemented[0] = 1;
     38  strcpy(ThermostatNames[1],"Woodcock");
     39  ThermostatImplemented[1] = 1;
     40  strcpy(ThermostatNames[2],"Gaussian");
     41  ThermostatImplemented[2] = 1;
     42  strcpy(ThermostatNames[3],"Langevin");
     43  ThermostatImplemented[3] = 1;
     44  strcpy(ThermostatNames[4],"Berendsen");
     45  ThermostatImplemented[4] = 1;
     46  strcpy(ThermostatNames[5],"NoseHoover");
     47  ThermostatImplemented[5] = 1;
     48
    2749  FastParsing = false;
    2850  ProcPEGamma=8;
    … …  
    3658  DoFullCurrent=0;
    3759  DoWannier=0;
     60  DoConstrainedMD=0;
    3861  CommonWannier=0;
    3962  SawtoothStart=0.01;
    … …  
    4265  UseAddGramSch=1;
    4366  Seed=1;
    44  
    4567  MaxOuterStep=0;
    46   Deltat=1;
     68  Deltat=0.01;
    4769  OutVisStep=10;
    4870  OutSrcStep=5;
    … …  
    96118  Free((void **)&configname, "config::~config: *configname");
    97119};
     120
     121/** Readin of Thermostat related values from parameter file.
     122 * \param *source parameter file
     123 */
     124void config::InitThermostats(ifstream *source)
     125{
     126  char *thermo = MallocString(12, "IonsInitRead: thermo");
     127  int verbose = 0;
     128
     129  // read desired Thermostat from file along with needed additional parameters
     130  if (ParseForParameter(verbose,source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
     131    if (strcmp(thermo, ThermostatNames[0]) == 0) { // None
     132      if (ThermostatImplemented[0] == 1) {
     133        Thermostat = None;
     134      } else {
     135        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     136        Thermostat = None;
     137      }
     138    } else if (strcmp(thermo, ThermostatNames[1]) == 0) { // Woodcock
     139      if (ThermostatImplemented[1] == 1) {
     140        Thermostat = Woodcock;
     141        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read scaling frequency
     142      } else {
     143        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     144        Thermostat = None;
     145      }
     146    } else if (strcmp(thermo, ThermostatNames[2]) == 0) { // Gaussian
     147      if (ThermostatImplemented[2] == 1) {
     148        Thermostat = Gaussian;
     149        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, int_type, &ScaleTempStep, 1, critical); // read collision rate
     150      } else {
     151        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     152        Thermostat = None;
     153      }
     154    } else if (strcmp(thermo, ThermostatNames[3]) == 0) { // Langevin
     155      if (ThermostatImplemented[3] == 1) {
     156        Thermostat = Langevin;
     157        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read gamma
     158        if (ParseForParameter(verbose,source,"Thermostat", 0, 3, 1, double_type, &alpha, 1, optional)) {
     159          cout << Verbose(2) << "Extended Stochastic Thermostat detected with interpolation coefficient " << alpha << "." << endl;
     160        } else {
     161          alpha = 1.;
     162        }
     163      } else {
     164        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     165        Thermostat = None;
     166      }
     167    } else if (strcmp(thermo, ThermostatNames[4]) == 0) { // Berendsen
     168      if (ThermostatImplemented[4] == 1) {
     169        Thermostat = Berendsen;
     170        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &TempFrequency, 1, critical); // read \tau_T
     171      } else {
     172        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     173        Thermostat = None;
     174      }
     175    } else if (strcmp(thermo, ThermostatNames[5]) == 0) { // Nose-Hoover
     176      if (ThermostatImplemented[5] == 1) {
     177        Thermostat = NoseHoover;
     178        ParseForParameter(verbose,source,"Thermostat", 0, 2, 1, double_type, &HooverMass, 1, critical); // read Hoovermass
     179        alpha = 0.;
     180      } else {
     181        cout << Verbose(1) << "Warning: " << ThermostatNames[0] << " thermostat not implemented, falling back to None." << endl;
     182        Thermostat = None;
     183      }
     184    } else {
     185      cout << Verbose(1) << " Warning: thermostat name was not understood!" << endl;
     186      Thermostat = None;
     187    }
     188  } else {
     189    if ((MaxOuterStep > 0) && (TargetTemp != 0))
     190      cout << Verbose(2) <<  "No thermostat chosen despite finite temperature MD, falling back to None." << endl;
     191    Thermostat = None;
     192  }
     193  Free((void **)&thermo, "InitThermostats: thermo");
     194};
     195
    98196
    99197/** Displays menu for editing each entry of the config file.
    … …  
    465563  double value[3];
    466564 
     565  InitThermostats(file);
     566 
    467567  /* Namen einlesen */
    468568
    … …  
    518618  }
    519619 
     620  if (ParseForParameter(verbose,file,"DoConstrainedMD", 0, 1, 1, int_type, &(config::DoConstrainedMD), 1, optional))
     621    if (config::DoConstrainedMD < 0)
     622      config::DoConstrainedMD = 0;
    520623  ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
    521624  if (!ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
    … …  
    9991102    *output << "DoPerturbation\t" << config::DoPerturbation << "\t# Do perturbation calculate and determine susceptibility and shielding" << endl;
    10001103    *output << "DoFullCurrent\t" << config::DoFullCurrent << "\t# Do full perturbation" << endl;
     1104    *output << "DoConstrainedMD\t" << config::DoConstrainedMD << "\t# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD" << endl;
     1105    *output << "Thermostat\t" << ThermostatNames[Thermostat] << "\t";
     1106    switch(Thermostat) {
     1107      default:
     1108      case None:
     1109        break;
     1110      case Woodcock:
     1111        *output << ScaleTempStep;
     1112        break;
     1113      case Gaussian:
     1114        *output << ScaleTempStep;
     1115        break;
     1116      case Langevin:
     1117        *output << TempFrequency << "\t" << alpha;
     1118        break;
     1119      case Berendsen:
     1120        *output << TempFrequency;
     1121        break;
     1122      case NoseHoover:
     1123        *output << HooverMass;
     1124        break;
     1125    };
     1126    *output << "\t# Which Thermostat and its parameters to use in MD case." << endl;
    10011127    *output << "CommonWannier\t" << config::CommonWannier << "\t# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center" << endl;
    10021128    *output << "SawtoothStart\t" << config::SawtoothStart << "\t# Absolute value for smooth transition at cell border " << endl;

  • src/datacreator.cpp

    rb38b64 r631dcb  
    276276  stringstream filename;
    277277  ofstream output;
    278   double norm = 0.;
    279278
    280279  filename << prefix << ".dat";

  • src/defs.hpp

    rb38b64 r631dcb  
    1616#define BONDTHRESHOLD 0.5   //!< CSD threshold in bond check which is the width of the interval whose center is the sum of the covalent radii
    1717#define AtomicEnergyToKelvin 315774.67  //!< conversion factor from atomic energy to kelvin via boltzmann factor
     18#define KelvinToAtomicTemperature 3.1668152e-06    //!< conversion factor for Kelvin to atomic temperature (Hartree over k_B)
     19#define KelvinToeV 8.6173422e-05                   //!< conversion factor for Kelvin to Ht (k_B * T = energy), thus Boltzmann constant in eV/K
     20#define AtomicMassUnitsToeV 931494088.        //!< conversion factor for atomic weight in units to mass in eV
     21#define AtomicMassUnitsToHt 34480864.        //!< conversion factor for atomic weight in units to mass in Ht (protonmass/electronmass * electron_mass_in_Ht
     22#define ElectronMass_Ht 18778.865            //!< electron mass in Ht
     23#define ElectronMass_eV 510998.903           //!< electron mass in eV
     24#define Units2Electronmass (AtomicMassUnitsToeV/ElectronMass_eV) //!< atomic mass unit in eV/ electron mass in eV = 1 822.88863
     25#define Atomictime2Femtoseconds 0.024188843     //!< Atomictime in fs
     26
    1827#define VERSIONSTRING "v1.0"
    1928

  • src/joiner.cpp

    rb38b64 r631dcb  
    104104    if(!Shielding.ParseIndices(argv[1])) return 1;
    105105    if(!ShieldingPAS.ParseIndices(argv[1])) return 1;
    106     if(!Chi.ParseIndices()) return 1;
    107     if(!ChiPAS.ParseIndices()) return 1;
     106    //if(!Chi.ParseIndices()) return 1;
     107    //if(!ChiPAS.ParseIndices()) return 1;
    108108  }
    109109

  • src/moleculelist.cpp

    rb38b64 r631dcb  
    378378 * \param *configuration standard configuration to attach atoms in fragment molecule to.
    379379 * \param *SortIndex Index to map from the BFS labeling to the sequence how of Ion_Type in the config
     380 * \param DoPeriodic true - call ScanForPeriodicCorrection, false - don't
     381 * \param DoCentering true - call molecule::CenterEdge(), false - don't
    380382 * \return true - success (each file was written), false - something went wrong.
    381383 */
    382 bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
     384bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, const char *fragmentprefix, config *configuration, int *SortIndex, bool DoPeriodic, bool DoCentering)
    383385{
    384386  ofstream outputFragment;
    … …  
    405407
    406408    // correct periodic
    407     //ListOfMolecules[i]->ScanForPeriodicCorrection(out);
     409    if (DoPeriodic)
     410      ListOfMolecules[i]->ScanForPeriodicCorrection(out);
    408411
    409412    // output xyz file
    410413    FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
    411     sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     414    sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, fragmentprefix, FragmentNumber);
    412415    outputFragment.open(FragmentName, ios::out);
    413     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
     416    *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
    414417    if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
    415418      *out << " done." << endl;
    … …  
    430433   
    431434    // center on edge
    432     ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
    433     ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
    434     int j = -1;
    435     for (int k=0;k<NDIM;k++) {
    436       j += k+1;
    437       BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
    438       ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
    439     }
    440     ListOfMolecules[i]->Translate(&BoxDimension);
     435    if (DoCentering) {
     436      ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
     437      ListOfMolecules[i]->SetBoxDimension(&BoxDimension);  // update Box of atoms by boundary
     438      int j = -1;
     439      for (int k=0;k<NDIM;k++) {
     440        j += k+1;
     441        BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
     442        ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
     443      }
     444      ListOfMolecules[i]->Translate(&BoxDimension);
     445    }
    441446
    442447    // also calculate necessary orbitals
    … …  
    446451    // change path in config
    447452    //strcpy(PathBackup, configuration->configpath);
    448     sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
     453    sprintf(FragmentName, "%s/%s%s/", PathBackup, fragmentprefix, FragmentNumber);
    449454    configuration->SetDefaultPath(FragmentName);
    450455   
    451456    // and save as config
    452     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    453     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
     457    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, fragmentprefix, FragmentNumber);
     458    *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
    454459    if ((intermediateResult = configuration->Save(FragmentName, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
    455460      *out << " done." << endl;
    … …  
    463468
    464469    // and save as mpqc input file
    465     sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
    466     *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
     470    sprintf(FragmentName, "%s/%s%s.in", configuration->configpath, fragmentprefix, FragmentNumber);
     471    *out << Verbose(2) << "Saving " << fragmentprefix << " No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
    467472    if ((intermediateResult = configuration->SaveMPQC(FragmentName, ListOfMolecules[i])))
    468473      *out << " done." << endl;

  • 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;

  • src/molecules.hpp

    rb38b64 r631dcb  
    1010
    1111// GSL headers
     12#include <gsl/gsl_eigen.h>
     13#include <gsl/gsl_heapsort.h>
     14#include <gsl/gsl_linalg.h>
     15#include <gsl/gsl_matrix.h>
    1216#include <gsl/gsl_multimin.h>
    1317#include <gsl/gsl_vector.h>
    14 #include <gsl/gsl_matrix.h>
    15 #include <gsl/gsl_eigen.h>
    16 #include <gsl/gsl_heapsort.h>
     18#include <gsl/gsl_randist.h>
    1719
    1820// STL headers
    … …  
    152154};
    153155
    154 ostream & operator << (ostream &ost, atom &a);
     156ostream & operator << (ostream &ost, const atom &a);
    155157
    156158/** Bonds between atoms.
    … …  
    191193};
    192194
    193 ostream & operator << (ostream &ost, bond &b);
     195
     196ostream & operator << (ostream &ost, const bond &b);
    194197
    195198class MoleculeLeafClass;
     199
     200
     201#define MaxThermostats 6      //!< maximum number of thermostat entries in Ions#ThermostatNames and Ions#ThermostatImplemented
     202enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
     203
    196204
    197205/** The complete molecule.
    … …  
    259267        void PrincipalAxisSystem(ofstream *out, bool DoRotate);
    260268        double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
    261         bool VerletForceIntegration(char *file, double delta_t, bool IsAngstroem);
    262 
     269       
     270  bool VerletForceIntegration(ofstream *out, char *file, config &configuration);
     271  void Thermostats(config &configuration, double ActualTemp, int Thermostat);
     272
     273  double ConstrainedPotential(ofstream *out, atom **permutation, int start, int end, double *constants, bool IsAngstroem);
     274  double MinimiseConstrainedPotential(ofstream *out, atom **&permutation, int startstep, int endstep, bool IsAngstroem);
     275  void EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force);
     276  bool LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration);
     277       
    263278  bool CheckBounds(const Vector *x) const;
    264279  void GetAlignvector(struct lsq_params * par) const;
    … …  
    339354  bool AddHydrogenCorrection(ofstream *out, char *path);
    340355  bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
    341   bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
     356  bool OutputConfigForListOfFragments(ofstream *out, const char *fragmentprefix, config *configuration, int *SortIndex, bool DoPeriodic, bool DoCentering);
    342357  void Output(ofstream *out);
    343358
    … …  
    389404    string basis;
    390405
     406    int DoConstrainedMD;
     407    int MaxOuterStep;
     408    int Thermostat;
     409    int *ThermostatImplemented;
     410    char **ThermostatNames;
     411    double TempFrequency;
     412    double alpha;
     413    double HooverMass;
     414    double TargetTemp;
     415    int ScaleTempStep;
     416   
    391417    private:
    392418    char *mainname;
    … …  
    409435    int Seed;
    410436
    411     int MaxOuterStep;
    412437    int OutVisStep;
    413438    int OutSrcStep;
    414     double TargetTemp;
    415     int ScaleTempStep;
    416439    int MaxPsiStep;
    417440    double EpsWannier;
    … …  
    460483  char *GetDefaultPath() const;
    461484  void SetDefaultPath(const char *path);
     485  void InitThermostats(ifstream *source);
    462486};
    463487
Note: See TracChangeset for help on using the changeset viewer.