/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010 University of Bonn. All rights reserved. * Please see the LICENSE file or "Copyright notice" in builder.cpp for details. */ /* * molecule_dynamics.cpp * * Created on: Oct 5, 2009 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "CodePatterns/MemDebug.hpp" #include "World.hpp" #include "atom.hpp" #include "config.hpp" #include "element.hpp" #include "CodePatterns/Info.hpp" #include "CodePatterns/Verbose.hpp" #include "CodePatterns/Log.hpp" #include "molecule.hpp" #include "parser.hpp" #include "LinearAlgebra/Plane.hpp" #include "ThermoStatContainer.hpp" #include "Thermostats/Berendsen.hpp" #include "CodePatterns/enumeration.hpp" #include #include #include /************************************* Functions for class molecule *********************************/ /** Penalizes long trajectories. * \param *Walker atom to check against others * \param *mol molecule with other atoms * \param &Params constraint potential parameters * \return penalty times each distance */ double SumDistanceOfTrajectories(atom *Walker, molecule *mol, struct EvaluatePotential &Params) { gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM); gsl_vector *x = gsl_vector_alloc(NDIM); atom *Sprinter = NULL; Vector trajectory1, trajectory2, normal, TestVector; double Norm1, Norm2, tmp, result = 0.; for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) { if ((*iter) == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; igetNr()]; // find first target point trajectory1 = Sprinter->getPositionAtStep(Params.endstep) - Walker->getPositionAtStep(Params.startstep); trajectory1.Normalize(); Norm1 = trajectory1.Norm(); Sprinter = Params.PermutationMap[(*iter)->getNr()]; // find second target point trajectory2 = Sprinter->getPositionAtStep(Params.endstep) - (*iter)->getPositionAtStep(Params.startstep); trajectory2.Normalize(); Norm2 = trajectory1.Norm(); // check whether either is zero() if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) { tmp = Walker->getPositionAtStep(Params.startstep).distance((*iter)->getPositionAtStep(Params.startstep)); } else if (Norm1 < MYEPSILON) { Sprinter = Params.PermutationMap[Walker->getNr()]; // find first target point trajectory1 = Sprinter->getPositionAtStep(Params.endstep) - (*iter)->getPositionAtStep(Params.startstep); trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything trajectory1 -= trajectory2; // project the part in norm direction away tmp = trajectory1.Norm(); // remaining norm is distance } else if (Norm2 < MYEPSILON) { Sprinter = Params.PermutationMap[(*iter)->getNr()]; // find second target point trajectory2 = Sprinter->getPositionAtStep(Params.endstep) - Walker->getPositionAtStep(Params.startstep); // copy second offset trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything trajectory2 -= trajectory1; // project the part in norm direction away tmp = trajectory2.Norm(); // remaining norm is distance } else if ((fabs(trajectory1.ScalarProduct(trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent // Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: "; // Log() << Verbose(0) << trajectory1; // Log() << Verbose(0) << " and "; // Log() << Verbose(0) << trajectory2; tmp = Walker->getPositionAtStep(Params.startstep).distance((*iter)->getPositionAtStep(Params.startstep)); // Log() << Verbose(0) << " with distance " << tmp << "." << endl; } else { // determine distance by finding minimum distance // Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear independent "; // Log() << Verbose(0) << endl; // Log() << Verbose(0) << "First Trajectory: "; // Log() << Verbose(0) << trajectory1 << endl; // Log() << Verbose(0) << "Second Trajectory: "; // Log() << Verbose(0) << trajectory2 << endl; // determine normal vector for both normal = Plane(trajectory1, trajectory2,0).getNormal(); // print all vectors for debugging // Log() << Verbose(0) << "Normal vector in between: "; // Log() << Verbose(0) << normal << endl; // setup matrix for (int i=NDIM;i--;) { gsl_matrix_set(A, 0, i, trajectory1[i]); gsl_matrix_set(A, 1, i, trajectory2[i]); gsl_matrix_set(A, 2, i, normal[i]); gsl_vector_set(x,i, (Walker->getPositionAtStep(Params.startstep)[i] - (*iter)->getPositionAtStep(Params.startstep)[i])); } // solve the linear system by Householder transformations gsl_linalg_HH_svx(A, x); // distance from last component tmp = gsl_vector_get(x,2); // Log() << Verbose(0) << " with distance " << tmp << "." << endl; // test whether we really have the intersection (by checking on c_1 and c_2) trajectory1.Scale(gsl_vector_get(x,0)); trajectory2.Scale(gsl_vector_get(x,1)); normal.Scale(gsl_vector_get(x,2)); TestVector = (*iter)->getPositionAtStep(Params.startstep) + trajectory2 + normal - (Walker->getPositionAtStep(Params.startstep) + trajectory1); if (TestVector.Norm() < MYEPSILON) { // Log() << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl; } else { // Log() << Verbose(2) << "Test: failed.\tIntersection is off by "; // Log() << Verbose(0) << TestVector; // Log() << Verbose(0) << "." << endl; } } // add up tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem; if (fabs(tmp) > MYEPSILON) { result += Params.PenaltyConstants[1] * 1./tmp; //Log() << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl; } } return result; }; /** Penalizes atoms heading to same target. * \param *Walker atom to check against others * \param *mol molecule with other atoms * \param &Params constrained potential parameters * \return \a penalty times the number of equal targets */ double PenalizeEqualTargets(atom *Walker, molecule *mol, struct EvaluatePotential &Params) { double result = 0.; for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) { if ((Params.PermutationMap[Walker->getNr()] == Params.PermutationMap[(*iter)->getNr()]) && (Walker->getNr() < (*iter)->getNr())) { // atom *Sprinter = PermutationMap[Walker->getNr()]; // Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at "; // Log() << Verbose(0) << Sprinter->getPosition(endstep); // Log() << Verbose(0) << ", penalting." << endl; result += Params.PenaltyConstants[2]; //Log() << Verbose(4) << "Adding " << constants[2] << "." << endl; } } return result; }; /** Evaluates the potential energy used for constrained molecular dynamics. * \f$V_i^{con} = c^{bond} \cdot | r_{P(i)} - R_i | + sum_{i \neq j} C^{min} \cdot \frac{1}{C_{ij}} + C^{inj} \Bigl (1 - \theta \bigl (\prod_{i \neq j} (P(i) - P(j)) \bigr ) \Bigr )\f$ * where the first term points to the target in minimum distance, the second is a penalty for trajectories lying too close to each other (\f$C_{ij}\f$ is minimum distance between * trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point. * Note that for the second term we have to solve the following linear system: * \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, * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$, * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines. * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration() * \param *out output stream for debugging * \param &Params constrained potential parameters * \return potential energy * \note This routine is scaling quadratically which is not optimal. * \todo There's a bit double counting going on for the first time, bu nothing to worry really about. */ double molecule::ConstrainedPotential(struct EvaluatePotential &Params) { double tmp = 0.; double result = 0.; // go through every atom atom *Runner = NULL; for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { // first term: distance to target Runner = Params.PermutationMap[(*iter)->getNr()]; // find target point tmp = ((*iter)->getPositionAtStep(Params.startstep).distance(Runner->getPositionAtStep(Params.endstep))); tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem; result += Params.PenaltyConstants[0] * tmp; //Log() << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl; // second term: sum of distances to other trajectories result += SumDistanceOfTrajectories((*iter), this, Params); // third term: penalty for equal targets result += PenalizeEqualTargets((*iter), this, Params); } return result; }; /** print the current permutation map. * \param *out output stream for debugging * \param &Params constrained potential parameters * \param AtomCount number of atoms */ void PrintPermutationMap(int AtomCount, struct EvaluatePotential &Params) { stringstream zeile1, zeile2; int *DoubleList = new int[AtomCount]; for(int i=0;igetNr()]++; zeile1 << i << " "; zeile2 << Params.PermutationMap[i]->getNr() << " "; } for (int i=0;i 1) doubles++; if (doubles >0) DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl); delete[](DoubleList); // Log() << Verbose(2) << zeile1.str() << endl << zeile2.str() << endl; }; /** \f$O(N^2)\f$ operation of calculation distance between each atom pair and putting into DistanceList. * \param *mol molecule to scan distances in * \param &Params constrained potential parameters */ void FillDistanceList(molecule *mol, struct EvaluatePotential &Params) { for (int i=mol->getAtomCount(); i--;) { Params.DistanceList[i] = new DistanceMap; // is the distance sorted target list per atom Params.DistanceList[i]->clear(); } for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) { for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) { Params.DistanceList[(*iter)->getNr()]->insert( DistancePair((*iter)->getPositionAtStep(Params.startstep).distance((*runner)->getPositionAtStep(Params.endstep)), (*runner)) ); } } }; /** initialize lists. * \param *out output stream for debugging * \param *mol molecule to scan distances in * \param &Params constrained potential parameters */ void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params) { for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) { Params.StepList[(*iter)->getNr()] = Params.DistanceList[(*iter)->getNr()]->begin(); // stores the step to the next iterator that could be a possible next target Params.PermutationMap[(*iter)->getNr()] = Params.DistanceList[(*iter)->getNr()]->begin()->second; // always pick target with the smallest distance Params.DoubleList[Params.DistanceList[(*iter)->getNr()]->begin()->second->getNr()]++; // increase this target's source count (>1? not injective) Params.DistanceIterators[(*iter)->getNr()] = Params.DistanceList[(*iter)->getNr()]->begin(); // and remember which one we picked DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << Params.DistanceList[(*iter)->getNr()]->begin()->first << "." << endl); } }; /** Try the next nearest neighbour in order to make the permutation map injective. * \param *out output stream for debugging * \param *mol molecule * \param *Walker atom to change its target * \param &OldPotential old value of constraint potential to see if we do better with new target * \param &Params constrained potential parameters */ double TryNextNearestNeighbourForInjectivePermutation(molecule *mol, atom *Walker, double &OldPotential, struct EvaluatePotential &Params) { double Potential = 0; DistanceMap::iterator NewBase = Params.DistanceIterators[Walker->getNr()]; // store old base do { NewBase++; // take next further distance in distance to targets list that's a target of no one } while ((Params.DoubleList[NewBase->second->getNr()] != 0) && (NewBase != Params.DistanceList[Walker->getNr()]->end())); if (NewBase != Params.DistanceList[Walker->getNr()]->end()) { Params.PermutationMap[Walker->getNr()] = NewBase->second; Potential = fabs(mol->ConstrainedPotential(Params)); if (Potential > OldPotential) { // undo Params.PermutationMap[Walker->getNr()] = Params.DistanceIterators[Walker->getNr()]->second; } else { // do Params.DoubleList[Params.DistanceIterators[Walker->getNr()]->second->getNr()]--; // decrease the old entry in the doubles list Params.DoubleList[NewBase->second->getNr()]++; // increase the old entry in the doubles list Params.DistanceIterators[Walker->getNr()] = NewBase; OldPotential = Potential; DoLog(3) && (Log() << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl); } } return Potential; }; /** Permutes \a **&PermutationMap until the penalty is below constants[2]. * \param *out output stream for debugging * \param *mol molecule to scan distances in * \param &Params constrained potential parameters */ void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params) { molecule::const_iterator iter = mol->begin(); DistanceMap::iterator NewBase; double Potential = fabs(mol->ConstrainedPotential(Params)); if (mol->empty()) { eLog() << Verbose(1) << "Molecule is empty." << endl; return; } while ((Potential) > Params.PenaltyConstants[2]) { PrintPermutationMap(mol->getAtomCount(), Params); iter++; if (iter == mol->end()) // round-robin at the end iter = mol->begin(); if (Params.DoubleList[Params.DistanceIterators[(*iter)->getNr()]->second->getNr()] <= 1) // no need to make those injective that aren't continue; // now, try finding a new one Potential = TryNextNearestNeighbourForInjectivePermutation(mol, (*iter), Potential, Params); } for (int i=mol->getAtomCount(); i--;) // now each single entry in the DoubleList should be <=1 if (Params.DoubleList[i] > 1) { DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl); performCriticalExit(); } DoLog(1) && (Log() << Verbose(1) << "done." << endl); }; /** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy. * We do the following: * -# Generate a distance list from all source to all target points * -# Sort this per source point * -# Take for each source point the target point with minimum distance, use this as initial permutation * -# check whether molecule::ConstrainedPotential() is greater than injective penalty * -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential. * -# Next, we only apply transformations that keep the injectivity of the permutations list. * -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target * point and try to change it for one with lesser distance, or for the next one with greater distance, but only * if this decreases the conditional potential. * -# finished. * -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree, * that the total force is always pointing in direction of the constraint force (ensuring that we move in the * right direction). * -# Finally, we calculate the potential energy and return. * \param *out output stream for debugging * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated) * \param endstep step giving final position in constrained MD * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false) * \sa molecule::VerletForceIntegration() * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2) * \todo The constrained potential's constants are set to fixed values right now, but they should scale based on checks of the system in order * to ensure they're properties (e.g. constants[2] always greater than the energy of the system). * \bug this all is not O(N log N) but O(N^2) */ double molecule::MinimiseConstrainedPotential(atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem) { double Potential, OldPotential, OlderPotential; struct EvaluatePotential Params; Params.PermutationMap = new atom *[getAtomCount()]; Params.DistanceList = new DistanceMap *[getAtomCount()]; Params.DistanceIterators = new DistanceMap::iterator[getAtomCount()]; Params.DoubleList = new int[getAtomCount()]; Params.StepList = new DistanceMap::iterator[getAtomCount()]; int round; atom *Sprinter = NULL; DistanceMap::iterator Rider, Strider; // set to zero for (int i=0;i target) CreateInitialLists(this, Params); // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it DoLog(1) && (Log() << Verbose(1) << "Making the PermutationMap injective ... " << endl); MakeInjectivePermutation(this, Params); delete[](Params.DoubleList); // argument minimise the constrained potential in this injective PermutationMap DoLog(1) && (Log() << Verbose(1) << "Argument minimising the PermutationMap." << endl); OldPotential = 1e+10; round = 0; do { DoLog(2) && (Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl); OlderPotential = OldPotential; molecule::const_iterator iter; do { iter = begin(); for (; iter != end(); ++iter) { PrintPermutationMap(getAtomCount(), Params); Sprinter = Params.DistanceIterators[(*iter)->getNr()]->second; // store initial partner Strider = Params.DistanceIterators[(*iter)->getNr()]; //remember old iterator Params.DistanceIterators[(*iter)->getNr()] = Params.StepList[(*iter)->getNr()]; if (Params.DistanceIterators[(*iter)->getNr()] == Params.DistanceList[(*iter)->getNr()]->end()) {// stop, before we run through the list and still on Params.DistanceIterators[(*iter)->getNr()] == Params.DistanceList[(*iter)->getNr()]->begin(); break; } //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->getNr()]->second << "." << endl; // find source of the new target molecule::const_iterator runner = begin(); for (; runner != end(); ++runner) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already) if (Params.PermutationMap[(*runner)->getNr()] == Params.DistanceIterators[(*iter)->getNr()]->second) { //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->getNr()] << "." << endl; break; } } if (runner != end()) { // we found the other source // then look in its distance list for Sprinter Rider = Params.DistanceList[(*runner)->getNr()]->begin(); for (; Rider != Params.DistanceList[(*runner)->getNr()]->end(); Rider++) if (Rider->second == Sprinter) break; if (Rider != Params.DistanceList[(*runner)->getNr()]->end()) { // if we have found one //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->getNr()] << "/" << *Rider->second << "." << endl; // exchange both Params.PermutationMap[(*iter)->getNr()] = Params.DistanceIterators[(*iter)->getNr()]->second; // put next farther distance into PermutationMap Params.PermutationMap[(*runner)->getNr()] = Sprinter; // and hand the old target to its respective owner PrintPermutationMap(getAtomCount(), Params); // calculate the new potential //Log() << Verbose(2) << "Checking new potential ..." << endl; Potential = ConstrainedPotential(Params); if (Potential > OldPotential) { // we made everything worse! Undo ... //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl; //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *Params.DistanceIterators[(*runner)->getNr()]->second << "." << endl; // Undo for Runner (note, we haven't moved the iteration yet, we may use this) Params.PermutationMap[(*runner)->getNr()] = Params.DistanceIterators[(*runner)->getNr()]->second; // Undo for Walker Params.DistanceIterators[(*iter)->getNr()] = Strider; // take next farther distance target //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *Params.DistanceIterators[(*iter)->getNr()]->second << "." << endl; Params.PermutationMap[(*iter)->getNr()] = Params.DistanceIterators[(*iter)->getNr()]->second; } else { Params.DistanceIterators[(*runner)->getNr()] = Rider; // if successful also move the pointer in the iterator list DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl); OldPotential = Potential; } if (Potential > Params.PenaltyConstants[2]) { DoeLog(1) && (eLog()<< Verbose(1) << "The two-step permutation procedure did not maintain injectivity!" << endl); exit(255); } //Log() << Verbose(0) << endl; } else { DoeLog(1) && (eLog()<< Verbose(1) << **runner << " was not the owner of " << *Sprinter << "!" << endl); exit(255); } } else { Params.PermutationMap[(*iter)->getNr()] = Params.DistanceIterators[(*iter)->getNr()]->second; // new target has no source! } Params.StepList[(*iter)->getNr()]++; // take next farther distance target } } while (++iter != end()); } while ((OlderPotential - OldPotential) > 1e-3); DoLog(1) && (Log() << Verbose(1) << "done." << endl); /// free memory and return with evaluated potential for (int i=getAtomCount(); i--;) Params.DistanceList[i]->clear(); delete[](Params.DistanceList); delete[](Params.DistanceIterators); return ConstrainedPotential(Params); }; /** Evaluates the (distance-related part) of the constrained potential for the constrained forces. * \param *out output stream for debugging * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated) * \param endstep step giving final position in constrained MD * \param **PermutationMap mapping between the atom label of the initial and the final configuration * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation. * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential() */ void molecule::EvaluateConstrainedForces(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force) { double constant = 10.; /// evaluate forces (only the distance to target dependent part) with the final PermutationMap DoLog(1) && (Log() << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl); BOOST_FOREACH(atom *_atom, atoms) { atom *Sprinter = PermutationMap[_atom->getNr()]; // set forces for (int i=NDIM;i++;) Force->Matrix[0][_atom->getNr()][5+i] += 2.*constant*sqrt(_atom->getPositionAtStep(startstep).distance(Sprinter->getPositionAtStep(endstep))); } DoLog(1) && (Log() << Verbose(1) << "done." << endl); }; /** Performs a linear interpolation between two desired atomic configurations with a given number of steps. * Note, step number is config::MaxOuterStep * \param *out output stream for debugging * \param startstep stating initial configuration in molecule::Trajectories * \param endstep stating final configuration in molecule::Trajectories * \param &prefix path and prefix * \param &config configuration structure * \param MapByIdentity if true we just use the identity to map atoms in start config to end config, if not we find mapping by \sa MinimiseConstrainedPotential() * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories */ bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, std::string prefix, config &configuration, bool MapByIdentity) { // TODO: rewrite permutationMaps using enumeration objects molecule *mol = NULL; bool status = true; int MaxSteps = configuration.MaxOuterStep; MoleculeListClass *MoleculePerStep = new MoleculeListClass(World::getPointer()); // Get the Permutation Map by MinimiseConstrainedPotential atom **PermutationMap = NULL; atom *Sprinter = NULL; if (!MapByIdentity) MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem()); else { // TODO: construction of enumeration goes here PermutationMap = new atom *[getAtomCount()]; for(internal_iterator iter = atoms.begin(); iter != atoms.end();++iter){ PermutationMap[(*iter)->getNr()] = (*iter); } } // check whether we have sufficient space in Trajectories for each atom LOG(1, "STATUS: Extending each trajectory size to " << MaxSteps+1 << "."); for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::ResizeTrajectory),MaxSteps+1)); // push endstep to last one for_each(atoms.begin(),atoms.end(),boost::bind(&atom::CopyStepOnStep,_1,MaxSteps,endstep)); endstep = MaxSteps; // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule DoLog(1) && (Log() << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl); for (int step = 0; step <= MaxSteps; step++) { mol = World::getInstance().createMolecule(); MoleculePerStep->insert(mol); for (molecule::const_iterator iter = begin(); iter != end(); ++iter) { // add to molecule list Sprinter = mol->AddCopyAtom((*iter)); // add to Trajectories Vector temp = (*iter)->getPositionAtStep(startstep) + (PermutationMap[(*iter)->getNr()]->getPositionAtStep(endstep) - (*iter)->getPositionAtStep(startstep))*((double)step/(double)MaxSteps); Sprinter->setPosition(temp); (*iter)->setAtomicVelocityAtStep(step, zeroVec); (*iter)->setAtomicForceAtStep(step, zeroVec); //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl; } } MDSteps = MaxSteps+1; // otherwise new Trajectories' points aren't stored on save&exit // store the list to single step files int *SortIndex = new int[getAtomCount()]; for (int i=getAtomCount(); i--; ) SortIndex[i] = i; status = MoleculePerStep->OutputConfigForListOfFragments(prefix, SortIndex); delete[](SortIndex); // free and return delete[](PermutationMap); delete(MoleculePerStep); return status; }; /** Parses nuclear forces from file and performs Verlet integration. * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we * have to transform them). * This adds a new MD step to the config file. * \param *file filename * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained * \param offset offset in matrix file to the first force component * \return true - file found and parsed, false - file not found or imparsable * \todo This is not yet checked if it is correctly working with DoConstrained set to true. */ bool molecule::VerletForceIntegration(char *file, config &configuration, const size_t offset) { Info FunctionInfo(__func__); string token; stringstream item; double IonMass, ConstrainedPotentialEnergy, ActualTemp; Vector Velocity; ForceMatrix Force; const int AtomCount = getAtomCount(); // parse file into ForceMatrix std::ifstream input(file); if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) { DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl); performCriticalExit(); return false; } input.close(); if (Force.RowCounter[0] != AtomCount) { DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl); performCriticalExit(); return false; } // correct Forces Velocity.Zero(); for(int i=0;i(AtomCount); } // solve a constrained potential if we are meant to if (configuration.DoConstrainedMD) { // calculate forces and potential atom **PermutationMap = NULL; ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem()); EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force); delete[](PermutationMap); } // and perform Verlet integration for each atom with position, velocity and force vector // check size of vectors BOOST_FOREACH(atom *_atom, atoms) { _atom->VelocityVerletUpdate(_atom->getNr(), MDSteps+1, &configuration, &Force, (const size_t) 0); } // correct velocities (rather momenta) so that center of mass remains motionless Velocity = atoms.totalMomentumAtStep(MDSteps+1); IonMass = atoms.totalMass(); // correct velocities (rather momenta) so that center of mass remains motionless Velocity *= 1./IonMass; atoms.addVelocityAtStep(-1*Velocity,MDSteps+1); ActualTemp = atoms.totalTemperatureAtStep(MDSteps+1); Berendsen berendsen = Berendsen(); berendsen.addToContainer(World::getInstance().getThermostats()); double ekin = berendsen.scaleAtoms(MDSteps,ActualTemp,atoms); DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl); MDSteps++; // exit return true; };