1 | /*
|
---|
2 | * molecule_dynamics.cpp
|
---|
3 | *
|
---|
4 | * Created on: Oct 5, 2009
|
---|
5 | * Author: heber
|
---|
6 | */
|
---|
7 |
|
---|
8 | #include "World.hpp"
|
---|
9 | #include "atom.hpp"
|
---|
10 | #include "config.hpp"
|
---|
11 | #include "element.hpp"
|
---|
12 | #include "info.hpp"
|
---|
13 | #include "log.hpp"
|
---|
14 | #include "memoryallocator.hpp"
|
---|
15 | #include "molecule.hpp"
|
---|
16 | #include "parser.hpp"
|
---|
17 | #include "Plane.hpp"
|
---|
18 |
|
---|
19 | /************************************* Functions for class molecule *********************************/
|
---|
20 |
|
---|
21 | /** Penalizes long trajectories.
|
---|
22 | * \param *Walker atom to check against others
|
---|
23 | * \param *mol molecule with other atoms
|
---|
24 | * \param &Params constraint potential parameters
|
---|
25 | * \return penalty times each distance
|
---|
26 | */
|
---|
27 | double SumDistanceOfTrajectories(atom *Walker, molecule *mol, struct EvaluatePotential &Params)
|
---|
28 | {
|
---|
29 | gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
|
---|
30 | gsl_vector *x = gsl_vector_alloc(NDIM);
|
---|
31 | atom *Sprinter = NULL;
|
---|
32 | Vector trajectory1, trajectory2, normal, TestVector;
|
---|
33 | double Norm1, Norm2, tmp, result = 0.;
|
---|
34 |
|
---|
35 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
36 | if ((*iter) == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
|
---|
37 | break;
|
---|
38 | // determine normalized trajectories direction vector (n1, n2)
|
---|
39 | Sprinter = Params.PermutationMap[Walker->nr]; // find first target point
|
---|
40 | trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep);
|
---|
41 | trajectory1.Normalize();
|
---|
42 | Norm1 = trajectory1.Norm();
|
---|
43 | Sprinter = Params.PermutationMap[(*iter)->nr]; // find second target point
|
---|
44 | trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
|
---|
45 | trajectory2.Normalize();
|
---|
46 | Norm2 = trajectory1.Norm();
|
---|
47 | // check whether either is zero()
|
---|
48 | if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
|
---|
49 | tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
|
---|
50 | } else if (Norm1 < MYEPSILON) {
|
---|
51 | Sprinter = Params.PermutationMap[Walker->nr]; // find first target point
|
---|
52 | trajectory1 = Sprinter->Trajectory.R.at(Params.endstep) - (*iter)->Trajectory.R.at(Params.startstep);
|
---|
53 | trajectory2 *= trajectory1.ScalarProduct(trajectory2); // trajectory2 is scaled to unity, hence we don't need to divide by anything
|
---|
54 | trajectory1 -= trajectory2; // project the part in norm direction away
|
---|
55 | tmp = trajectory1.Norm(); // remaining norm is distance
|
---|
56 | } else if (Norm2 < MYEPSILON) {
|
---|
57 | Sprinter = Params.PermutationMap[(*iter)->nr]; // find second target point
|
---|
58 | trajectory2 = Sprinter->Trajectory.R.at(Params.endstep) - Walker->Trajectory.R.at(Params.startstep); // copy second offset
|
---|
59 | trajectory1 *= trajectory2.ScalarProduct(trajectory1); // trajectory1 is scaled to unity, hence we don't need to divide by anything
|
---|
60 | trajectory2 -= trajectory1; // project the part in norm direction away
|
---|
61 | tmp = trajectory2.Norm(); // remaining norm is distance
|
---|
62 | } else if ((fabs(trajectory1.ScalarProduct(trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
|
---|
63 | // Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
|
---|
64 | // Log() << Verbose(0) << trajectory1;
|
---|
65 | // Log() << Verbose(0) << " and ";
|
---|
66 | // Log() << Verbose(0) << trajectory2;
|
---|
67 | tmp = Walker->Trajectory.R.at(Params.startstep).distance((*iter)->Trajectory.R.at(Params.startstep));
|
---|
68 | // Log() << Verbose(0) << " with distance " << tmp << "." << endl;
|
---|
69 | } else { // determine distance by finding minimum distance
|
---|
70 | // Log() << Verbose(3) << "Both trajectories of " << *Walker << " and " << *(*iter) << " are linear independent ";
|
---|
71 | // Log() << Verbose(0) << endl;
|
---|
72 | // Log() << Verbose(0) << "First Trajectory: ";
|
---|
73 | // Log() << Verbose(0) << trajectory1 << endl;
|
---|
74 | // Log() << Verbose(0) << "Second Trajectory: ";
|
---|
75 | // Log() << Verbose(0) << trajectory2 << endl;
|
---|
76 | // determine normal vector for both
|
---|
77 | normal = Plane(trajectory1, trajectory2,0).getNormal();
|
---|
78 | // print all vectors for debugging
|
---|
79 | // Log() << Verbose(0) << "Normal vector in between: ";
|
---|
80 | // Log() << Verbose(0) << normal << endl;
|
---|
81 | // setup matrix
|
---|
82 | for (int i=NDIM;i--;) {
|
---|
83 | gsl_matrix_set(A, 0, i, trajectory1[i]);
|
---|
84 | gsl_matrix_set(A, 1, i, trajectory2[i]);
|
---|
85 | gsl_matrix_set(A, 2, i, normal[i]);
|
---|
86 | gsl_vector_set(x,i, (Walker->Trajectory.R.at(Params.startstep)[i] - (*iter)->Trajectory.R.at(Params.startstep)[i]));
|
---|
87 | }
|
---|
88 | // solve the linear system by Householder transformations
|
---|
89 | gsl_linalg_HH_svx(A, x);
|
---|
90 | // distance from last component
|
---|
91 | tmp = gsl_vector_get(x,2);
|
---|
92 | // Log() << Verbose(0) << " with distance " << tmp << "." << endl;
|
---|
93 | // test whether we really have the intersection (by checking on c_1 and c_2)
|
---|
94 | trajectory1.Scale(gsl_vector_get(x,0));
|
---|
95 | trajectory2.Scale(gsl_vector_get(x,1));
|
---|
96 | normal.Scale(gsl_vector_get(x,2));
|
---|
97 | TestVector = (*iter)->Trajectory.R.at(Params.startstep) + trajectory2 + normal
|
---|
98 | - (Walker->Trajectory.R.at(Params.startstep) + trajectory1);
|
---|
99 | if (TestVector.Norm() < MYEPSILON) {
|
---|
100 | // Log() << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;
|
---|
101 | } else {
|
---|
102 | // Log() << Verbose(2) << "Test: failed.\tIntersection is off by ";
|
---|
103 | // Log() << Verbose(0) << TestVector;
|
---|
104 | // Log() << Verbose(0) << "." << endl;
|
---|
105 | }
|
---|
106 | }
|
---|
107 | // add up
|
---|
108 | tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
|
---|
109 | if (fabs(tmp) > MYEPSILON) {
|
---|
110 | result += Params.PenaltyConstants[1] * 1./tmp;
|
---|
111 | //Log() << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
|
---|
112 | }
|
---|
113 | }
|
---|
114 | return result;
|
---|
115 | };
|
---|
116 |
|
---|
117 | /** Penalizes atoms heading to same target.
|
---|
118 | * \param *Walker atom to check against others
|
---|
119 | * \param *mol molecule with other atoms
|
---|
120 | * \param &Params constrained potential parameters
|
---|
121 | * \return \a penalty times the number of equal targets
|
---|
122 | */
|
---|
123 | double PenalizeEqualTargets(atom *Walker, molecule *mol, struct EvaluatePotential &Params)
|
---|
124 | {
|
---|
125 | double result = 0.;
|
---|
126 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
127 | if ((Params.PermutationMap[Walker->nr] == Params.PermutationMap[(*iter)->nr]) && (Walker->nr < (*iter)->nr)) {
|
---|
128 | // atom *Sprinter = PermutationMap[Walker->nr];
|
---|
129 | // Log() << Verbose(0) << *Walker << " and " << *(*iter) << " are heading to the same target at ";
|
---|
130 | // Log() << Verbose(0) << Sprinter->Trajectory.R.at(endstep);
|
---|
131 | // Log() << Verbose(0) << ", penalting." << endl;
|
---|
132 | result += Params.PenaltyConstants[2];
|
---|
133 | //Log() << Verbose(4) << "Adding " << constants[2] << "." << endl;
|
---|
134 | }
|
---|
135 | }
|
---|
136 | return result;
|
---|
137 | };
|
---|
138 |
|
---|
139 | /** Evaluates the potential energy used for constrained molecular dynamics.
|
---|
140 | * \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$
|
---|
141 | * 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
|
---|
142 | * trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
|
---|
143 | * Note that for the second term we have to solve the following linear system:
|
---|
144 | * \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,
|
---|
145 | * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
|
---|
146 | * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
|
---|
147 | * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
|
---|
148 | * \param *out output stream for debugging
|
---|
149 | * \param &Params constrained potential parameters
|
---|
150 | * \return potential energy
|
---|
151 | * \note This routine is scaling quadratically which is not optimal.
|
---|
152 | * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
|
---|
153 | */
|
---|
154 | double molecule::ConstrainedPotential(struct EvaluatePotential &Params)
|
---|
155 | {
|
---|
156 | double tmp = 0.;
|
---|
157 | double result = 0.;
|
---|
158 | // go through every atom
|
---|
159 | atom *Runner = NULL;
|
---|
160 | for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
|
---|
161 | // first term: distance to target
|
---|
162 | Runner = Params.PermutationMap[(*iter)->nr]; // find target point
|
---|
163 | tmp = ((*iter)->Trajectory.R.at(Params.startstep).distance(Runner->Trajectory.R.at(Params.endstep)));
|
---|
164 | tmp *= Params.IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
|
---|
165 | result += Params.PenaltyConstants[0] * tmp;
|
---|
166 | //Log() << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
|
---|
167 |
|
---|
168 | // second term: sum of distances to other trajectories
|
---|
169 | result += SumDistanceOfTrajectories((*iter), this, Params);
|
---|
170 |
|
---|
171 | // third term: penalty for equal targets
|
---|
172 | result += PenalizeEqualTargets((*iter), this, Params);
|
---|
173 | }
|
---|
174 |
|
---|
175 | return result;
|
---|
176 | };
|
---|
177 |
|
---|
178 | /** print the current permutation map.
|
---|
179 | * \param *out output stream for debugging
|
---|
180 | * \param &Params constrained potential parameters
|
---|
181 | * \param AtomCount number of atoms
|
---|
182 | */
|
---|
183 | void PrintPermutationMap(int AtomCount, struct EvaluatePotential &Params)
|
---|
184 | {
|
---|
185 | stringstream zeile1, zeile2;
|
---|
186 | int *DoubleList = new int[AtomCount];
|
---|
187 | for(int i=0;i<AtomCount;i++)
|
---|
188 | DoubleList[i] = 0;
|
---|
189 | int doubles = 0;
|
---|
190 | zeile1 << "PermutationMap: ";
|
---|
191 | zeile2 << " ";
|
---|
192 | for (int i=0;i<AtomCount;i++) {
|
---|
193 | Params.DoubleList[Params.PermutationMap[i]->nr]++;
|
---|
194 | zeile1 << i << " ";
|
---|
195 | zeile2 << Params.PermutationMap[i]->nr << " ";
|
---|
196 | }
|
---|
197 | for (int i=0;i<AtomCount;i++)
|
---|
198 | if (Params.DoubleList[i] > 1)
|
---|
199 | doubles++;
|
---|
200 | if (doubles >0)
|
---|
201 | DoLog(2) && (Log() << Verbose(2) << "Found " << doubles << " Doubles." << endl);
|
---|
202 | delete[](DoubleList);
|
---|
203 | // Log() << Verbose(2) << zeile1.str() << endl << zeile2.str() << endl;
|
---|
204 | };
|
---|
205 |
|
---|
206 | /** \f$O(N^2)\f$ operation of calculation distance between each atom pair and putting into DistanceList.
|
---|
207 | * \param *mol molecule to scan distances in
|
---|
208 | * \param &Params constrained potential parameters
|
---|
209 | */
|
---|
210 | void FillDistanceList(molecule *mol, struct EvaluatePotential &Params)
|
---|
211 | {
|
---|
212 | for (int i=mol->getAtomCount(); i--;) {
|
---|
213 | Params.DistanceList[i] = new DistanceMap; // is the distance sorted target list per atom
|
---|
214 | Params.DistanceList[i]->clear();
|
---|
215 | }
|
---|
216 |
|
---|
217 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
218 | for (molecule::const_iterator runner = mol->begin(); runner != mol->end(); ++runner) {
|
---|
219 | Params.DistanceList[(*iter)->nr]->insert( DistancePair((*iter)->Trajectory.R.at(Params.startstep).distance((*runner)->Trajectory.R.at(Params.endstep)), (*runner)) );
|
---|
220 | }
|
---|
221 | }
|
---|
222 | };
|
---|
223 |
|
---|
224 | /** initialize lists.
|
---|
225 | * \param *out output stream for debugging
|
---|
226 | * \param *mol molecule to scan distances in
|
---|
227 | * \param &Params constrained potential parameters
|
---|
228 | */
|
---|
229 | void CreateInitialLists(molecule *mol, struct EvaluatePotential &Params)
|
---|
230 | {
|
---|
231 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
232 | Params.StepList[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin(); // stores the step to the next iterator that could be a possible next target
|
---|
233 | Params.PermutationMap[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin()->second; // always pick target with the smallest distance
|
---|
234 | Params.DoubleList[Params.DistanceList[(*iter)->nr]->begin()->second->nr]++; // increase this target's source count (>1? not injective)
|
---|
235 | Params.DistanceIterators[(*iter)->nr] = Params.DistanceList[(*iter)->nr]->begin(); // and remember which one we picked
|
---|
236 | DoLog(2) && (Log() << Verbose(2) << **iter << " starts with distance " << Params.DistanceList[(*iter)->nr]->begin()->first << "." << endl);
|
---|
237 | }
|
---|
238 | };
|
---|
239 |
|
---|
240 | /** Try the next nearest neighbour in order to make the permutation map injective.
|
---|
241 | * \param *out output stream for debugging
|
---|
242 | * \param *mol molecule
|
---|
243 | * \param *Walker atom to change its target
|
---|
244 | * \param &OldPotential old value of constraint potential to see if we do better with new target
|
---|
245 | * \param &Params constrained potential parameters
|
---|
246 | */
|
---|
247 | double TryNextNearestNeighbourForInjectivePermutation(molecule *mol, atom *Walker, double &OldPotential, struct EvaluatePotential &Params)
|
---|
248 | {
|
---|
249 | double Potential = 0;
|
---|
250 | DistanceMap::iterator NewBase = Params.DistanceIterators[Walker->nr]; // store old base
|
---|
251 | do {
|
---|
252 | NewBase++; // take next further distance in distance to targets list that's a target of no one
|
---|
253 | } while ((Params.DoubleList[NewBase->second->nr] != 0) && (NewBase != Params.DistanceList[Walker->nr]->end()));
|
---|
254 | if (NewBase != Params.DistanceList[Walker->nr]->end()) {
|
---|
255 | Params.PermutationMap[Walker->nr] = NewBase->second;
|
---|
256 | Potential = fabs(mol->ConstrainedPotential(Params));
|
---|
257 | if (Potential > OldPotential) { // undo
|
---|
258 | Params.PermutationMap[Walker->nr] = Params.DistanceIterators[Walker->nr]->second;
|
---|
259 | } else { // do
|
---|
260 | Params.DoubleList[Params.DistanceIterators[Walker->nr]->second->nr]--; // decrease the old entry in the doubles list
|
---|
261 | Params.DoubleList[NewBase->second->nr]++; // increase the old entry in the doubles list
|
---|
262 | Params.DistanceIterators[Walker->nr] = NewBase;
|
---|
263 | OldPotential = Potential;
|
---|
264 | DoLog(3) && (Log() << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl);
|
---|
265 | }
|
---|
266 | }
|
---|
267 | return Potential;
|
---|
268 | };
|
---|
269 |
|
---|
270 | /** Permutes \a **&PermutationMap until the penalty is below constants[2].
|
---|
271 | * \param *out output stream for debugging
|
---|
272 | * \param *mol molecule to scan distances in
|
---|
273 | * \param &Params constrained potential parameters
|
---|
274 | */
|
---|
275 | void MakeInjectivePermutation(molecule *mol, struct EvaluatePotential &Params)
|
---|
276 | {
|
---|
277 | molecule::const_iterator iter = mol->begin();
|
---|
278 | DistanceMap::iterator NewBase;
|
---|
279 | double Potential = fabs(mol->ConstrainedPotential(Params));
|
---|
280 |
|
---|
281 | if (mol->empty()) {
|
---|
282 | eLog() << Verbose(1) << "Molecule is empty." << endl;
|
---|
283 | return;
|
---|
284 | }
|
---|
285 | while ((Potential) > Params.PenaltyConstants[2]) {
|
---|
286 | PrintPermutationMap(mol->getAtomCount(), Params);
|
---|
287 | iter++;
|
---|
288 | if (iter == mol->end()) // round-robin at the end
|
---|
289 | iter = mol->begin();
|
---|
290 | if (Params.DoubleList[Params.DistanceIterators[(*iter)->nr]->second->nr] <= 1) // no need to make those injective that aren't
|
---|
291 | continue;
|
---|
292 | // now, try finding a new one
|
---|
293 | Potential = TryNextNearestNeighbourForInjectivePermutation(mol, (*iter), Potential, Params);
|
---|
294 | }
|
---|
295 | for (int i=mol->getAtomCount(); i--;) // now each single entry in the DoubleList should be <=1
|
---|
296 | if (Params.DoubleList[i] > 1) {
|
---|
297 | DoeLog(0) && (eLog()<< Verbose(0) << "Failed to create an injective PermutationMap!" << endl);
|
---|
298 | performCriticalExit();
|
---|
299 | }
|
---|
300 | DoLog(1) && (Log() << Verbose(1) << "done." << endl);
|
---|
301 | };
|
---|
302 |
|
---|
303 | /** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
|
---|
304 | * We do the following:
|
---|
305 | * -# Generate a distance list from all source to all target points
|
---|
306 | * -# Sort this per source point
|
---|
307 | * -# Take for each source point the target point with minimum distance, use this as initial permutation
|
---|
308 | * -# check whether molecule::ConstrainedPotential() is greater than injective penalty
|
---|
309 | * -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
|
---|
310 | * -# Next, we only apply transformations that keep the injectivity of the permutations list.
|
---|
311 | * -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target
|
---|
312 | * point and try to change it for one with lesser distance, or for the next one with greater distance, but only
|
---|
313 | * if this decreases the conditional potential.
|
---|
314 | * -# finished.
|
---|
315 | * -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
|
---|
316 | * that the total force is always pointing in direction of the constraint force (ensuring that we move in the
|
---|
317 | * right direction).
|
---|
318 | * -# Finally, we calculate the potential energy and return.
|
---|
319 | * \param *out output stream for debugging
|
---|
320 | * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
|
---|
321 | * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
|
---|
322 | * \param endstep step giving final position in constrained MD
|
---|
323 | * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
|
---|
324 | * \sa molecule::VerletForceIntegration()
|
---|
325 | * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
|
---|
326 | * \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
|
---|
327 | * to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
|
---|
328 | * \bug this all is not O(N log N) but O(N^2)
|
---|
329 | */
|
---|
330 | double molecule::MinimiseConstrainedPotential(atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem)
|
---|
331 | {
|
---|
332 | double Potential, OldPotential, OlderPotential;
|
---|
333 | struct EvaluatePotential Params;
|
---|
334 | Params.PermutationMap = new atom *[getAtomCount()];
|
---|
335 | Params.DistanceList = new DistanceMap *[getAtomCount()];
|
---|
336 | Params.DistanceIterators = new DistanceMap::iterator[getAtomCount()];
|
---|
337 | Params.DoubleList = new int[getAtomCount()];
|
---|
338 | Params.StepList = new DistanceMap::iterator[getAtomCount()];
|
---|
339 | int round;
|
---|
340 | atom *Sprinter = NULL;
|
---|
341 | DistanceMap::iterator Rider, Strider;
|
---|
342 |
|
---|
343 | // set to zero
|
---|
344 | for (int i=0;i<getAtomCount();i++) {
|
---|
345 | Params.PermutationMap[i] = NULL;
|
---|
346 | Params.DoubleList[i] = 0;
|
---|
347 | }
|
---|
348 |
|
---|
349 | /// Minimise the potential
|
---|
350 | // set Lagrange multiplier constants
|
---|
351 | Params.PenaltyConstants[0] = 10.;
|
---|
352 | Params.PenaltyConstants[1] = 1.;
|
---|
353 | Params.PenaltyConstants[2] = 1e+7; // just a huge penalty
|
---|
354 | // generate the distance list
|
---|
355 | DoLog(1) && (Log() << Verbose(1) << "Allocating, initializting and filling the distance list ... " << endl);
|
---|
356 | FillDistanceList(this, Params);
|
---|
357 |
|
---|
358 | // create the initial PermutationMap (source -> target)
|
---|
359 | CreateInitialLists(this, Params);
|
---|
360 |
|
---|
361 | // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
|
---|
362 | DoLog(1) && (Log() << Verbose(1) << "Making the PermutationMap injective ... " << endl);
|
---|
363 | MakeInjectivePermutation(this, Params);
|
---|
364 | delete[](Params.DoubleList);
|
---|
365 |
|
---|
366 | // argument minimise the constrained potential in this injective PermutationMap
|
---|
367 | DoLog(1) && (Log() << Verbose(1) << "Argument minimising the PermutationMap." << endl);
|
---|
368 | OldPotential = 1e+10;
|
---|
369 | round = 0;
|
---|
370 | do {
|
---|
371 | DoLog(2) && (Log() << Verbose(2) << "Starting round " << ++round << ", at current potential " << OldPotential << " ... " << endl);
|
---|
372 | OlderPotential = OldPotential;
|
---|
373 | molecule::const_iterator iter;
|
---|
374 | do {
|
---|
375 | iter = begin();
|
---|
376 | for (; iter != end(); ++iter) {
|
---|
377 | PrintPermutationMap(getAtomCount(), Params);
|
---|
378 | Sprinter = Params.DistanceIterators[(*iter)->nr]->second; // store initial partner
|
---|
379 | Strider = Params.DistanceIterators[(*iter)->nr]; //remember old iterator
|
---|
380 | Params.DistanceIterators[(*iter)->nr] = Params.StepList[(*iter)->nr];
|
---|
381 | if (Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->end()) {// stop, before we run through the list and still on
|
---|
382 | Params.DistanceIterators[(*iter)->nr] == Params.DistanceList[(*iter)->nr]->begin();
|
---|
383 | break;
|
---|
384 | }
|
---|
385 | //Log() << Verbose(2) << "Current Walker: " << *(*iter) << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[(*iter)->nr]->second << "." << endl;
|
---|
386 | // find source of the new target
|
---|
387 | molecule::const_iterator runner = begin();
|
---|
388 | 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)
|
---|
389 | if (Params.PermutationMap[(*runner)->nr] == Params.DistanceIterators[(*iter)->nr]->second) {
|
---|
390 | //Log() << Verbose(2) << "Found the corresponding owner " << *(*runner) << " to " << *PermutationMap[(*runner)->nr] << "." << endl;
|
---|
391 | break;
|
---|
392 | }
|
---|
393 | }
|
---|
394 | if (runner != end()) { // we found the other source
|
---|
395 | // then look in its distance list for Sprinter
|
---|
396 | Rider = Params.DistanceList[(*runner)->nr]->begin();
|
---|
397 | for (; Rider != Params.DistanceList[(*runner)->nr]->end(); Rider++)
|
---|
398 | if (Rider->second == Sprinter)
|
---|
399 | break;
|
---|
400 | if (Rider != Params.DistanceList[(*runner)->nr]->end()) { // if we have found one
|
---|
401 | //Log() << Verbose(2) << "Current Other: " << *(*runner) << " with old/next candidate " << *PermutationMap[(*runner)->nr] << "/" << *Rider->second << "." << endl;
|
---|
402 | // exchange both
|
---|
403 | Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // put next farther distance into PermutationMap
|
---|
404 | Params.PermutationMap[(*runner)->nr] = Sprinter; // and hand the old target to its respective owner
|
---|
405 | PrintPermutationMap(getAtomCount(), Params);
|
---|
406 | // calculate the new potential
|
---|
407 | //Log() << Verbose(2) << "Checking new potential ..." << endl;
|
---|
408 | Potential = ConstrainedPotential(Params);
|
---|
409 | if (Potential > OldPotential) { // we made everything worse! Undo ...
|
---|
410 | //Log() << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
|
---|
411 | //Log() << Verbose(3) << "Setting " << *(*runner) << "'s source to " << *Params.DistanceIterators[(*runner)->nr]->second << "." << endl;
|
---|
412 | // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
|
---|
413 | Params.PermutationMap[(*runner)->nr] = Params.DistanceIterators[(*runner)->nr]->second;
|
---|
414 | // Undo for Walker
|
---|
415 | Params.DistanceIterators[(*iter)->nr] = Strider; // take next farther distance target
|
---|
416 | //Log() << Verbose(3) << "Setting " << *(*iter) << "'s source to " << *Params.DistanceIterators[(*iter)->nr]->second << "." << endl;
|
---|
417 | Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second;
|
---|
418 | } else {
|
---|
419 | Params.DistanceIterators[(*runner)->nr] = Rider; // if successful also move the pointer in the iterator list
|
---|
420 | DoLog(3) && (Log() << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl);
|
---|
421 | OldPotential = Potential;
|
---|
422 | }
|
---|
423 | if (Potential > Params.PenaltyConstants[2]) {
|
---|
424 | DoeLog(1) && (eLog()<< Verbose(1) << "The two-step permutation procedure did not maintain injectivity!" << endl);
|
---|
425 | exit(255);
|
---|
426 | }
|
---|
427 | //Log() << Verbose(0) << endl;
|
---|
428 | } else {
|
---|
429 | DoeLog(1) && (eLog()<< Verbose(1) << **runner << " was not the owner of " << *Sprinter << "!" << endl);
|
---|
430 | exit(255);
|
---|
431 | }
|
---|
432 | } else {
|
---|
433 | Params.PermutationMap[(*iter)->nr] = Params.DistanceIterators[(*iter)->nr]->second; // new target has no source!
|
---|
434 | }
|
---|
435 | Params.StepList[(*iter)->nr]++; // take next farther distance target
|
---|
436 | }
|
---|
437 | } while (++iter != end());
|
---|
438 | } while ((OlderPotential - OldPotential) > 1e-3);
|
---|
439 | DoLog(1) && (Log() << Verbose(1) << "done." << endl);
|
---|
440 |
|
---|
441 |
|
---|
442 | /// free memory and return with evaluated potential
|
---|
443 | for (int i=getAtomCount(); i--;)
|
---|
444 | Params.DistanceList[i]->clear();
|
---|
445 | delete[](Params.DistanceList);
|
---|
446 | delete[](Params.DistanceIterators);
|
---|
447 | return ConstrainedPotential(Params);
|
---|
448 | };
|
---|
449 |
|
---|
450 |
|
---|
451 | /** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
|
---|
452 | * \param *out output stream for debugging
|
---|
453 | * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
|
---|
454 | * \param endstep step giving final position in constrained MD
|
---|
455 | * \param **PermutationMap mapping between the atom label of the initial and the final configuration
|
---|
456 | * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
|
---|
457 | * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
|
---|
458 | */
|
---|
459 | void molecule::EvaluateConstrainedForces(int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
|
---|
460 | {
|
---|
461 | /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
|
---|
462 | DoLog(1) && (Log() << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl);
|
---|
463 | ActOnAllAtoms( &atom::EvaluateConstrainedForce, startstep, endstep, PermutationMap, Force );
|
---|
464 | DoLog(1) && (Log() << Verbose(1) << "done." << endl);
|
---|
465 | };
|
---|
466 |
|
---|
467 | /** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
|
---|
468 | * Note, step number is config::MaxOuterStep
|
---|
469 | * \param *out output stream for debugging
|
---|
470 | * \param startstep stating initial configuration in molecule::Trajectories
|
---|
471 | * \param endstep stating final configuration in molecule::Trajectories
|
---|
472 | * \param &config configuration structure
|
---|
473 | * \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()
|
---|
474 | * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
|
---|
475 | */
|
---|
476 | bool molecule::LinearInterpolationBetweenConfiguration(int startstep, int endstep, const char *prefix, config &configuration, bool MapByIdentity)
|
---|
477 | {
|
---|
478 | molecule *mol = NULL;
|
---|
479 | bool status = true;
|
---|
480 | int MaxSteps = configuration.MaxOuterStep;
|
---|
481 | MoleculeListClass *MoleculePerStep = new MoleculeListClass(World::getPointer());
|
---|
482 | // Get the Permutation Map by MinimiseConstrainedPotential
|
---|
483 | atom **PermutationMap = NULL;
|
---|
484 | atom *Sprinter = NULL;
|
---|
485 | if (!MapByIdentity)
|
---|
486 | MinimiseConstrainedPotential(PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
|
---|
487 | else {
|
---|
488 | PermutationMap = new atom *[getAtomCount()];
|
---|
489 | SetIndexedArrayForEachAtomTo( PermutationMap, &atom::nr );
|
---|
490 | }
|
---|
491 |
|
---|
492 | // check whether we have sufficient space in Trajectories for each atom
|
---|
493 | ActOnAllAtoms( &atom::ResizeTrajectory, MaxSteps );
|
---|
494 | // push endstep to last one
|
---|
495 | ActOnAllAtoms( &atom::CopyStepOnStep, MaxSteps, endstep );
|
---|
496 | endstep = MaxSteps;
|
---|
497 |
|
---|
498 | // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
|
---|
499 | DoLog(1) && (Log() << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl);
|
---|
500 | for (int step = 0; step <= MaxSteps; step++) {
|
---|
501 | mol = World::getInstance().createMolecule();
|
---|
502 | MoleculePerStep->insert(mol);
|
---|
503 | for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
|
---|
504 | // add to molecule list
|
---|
505 | Sprinter = mol->AddCopyAtom((*iter));
|
---|
506 | for (int n=NDIM;n--;) {
|
---|
507 | Sprinter->x[n] = (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
|
---|
508 | // add to Trajectories
|
---|
509 | //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
|
---|
510 | if (step < MaxSteps) {
|
---|
511 | (*iter)->Trajectory.R.at(step)[n] = (*iter)->Trajectory.R.at(startstep)[n] + (PermutationMap[(*iter)->nr]->Trajectory.R.at(endstep)[n] - (*iter)->Trajectory.R.at(startstep)[n])*((double)step/(double)MaxSteps);
|
---|
512 | (*iter)->Trajectory.U.at(step)[n] = 0.;
|
---|
513 | (*iter)->Trajectory.F.at(step)[n] = 0.;
|
---|
514 | }
|
---|
515 | }
|
---|
516 | }
|
---|
517 | }
|
---|
518 | MDSteps = MaxSteps+1; // otherwise new Trajectories' points aren't stored on save&exit
|
---|
519 |
|
---|
520 | // store the list to single step files
|
---|
521 | int *SortIndex = new int[getAtomCount()];
|
---|
522 | for (int i=getAtomCount(); i--; )
|
---|
523 | SortIndex[i] = i;
|
---|
524 | status = MoleculePerStep->OutputConfigForListOfFragments(&configuration, SortIndex);
|
---|
525 | delete[](SortIndex);
|
---|
526 |
|
---|
527 | // free and return
|
---|
528 | delete[](PermutationMap);
|
---|
529 | delete(MoleculePerStep);
|
---|
530 | return status;
|
---|
531 | };
|
---|
532 |
|
---|
533 | /** Parses nuclear forces from file and performs Verlet integration.
|
---|
534 | * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
|
---|
535 | * have to transform them).
|
---|
536 | * This adds a new MD step to the config file.
|
---|
537 | * \param *out output stream for debugging
|
---|
538 | * \param *file filename
|
---|
539 | * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
|
---|
540 | * \param delta_t time step width in atomic units
|
---|
541 | * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
|
---|
542 | * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential()
|
---|
543 | * \return true - file found and parsed, false - file not found or imparsable
|
---|
544 | * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
|
---|
545 | */
|
---|
546 | bool molecule::VerletForceIntegration(char *file, config &configuration)
|
---|
547 | {
|
---|
548 | Info FunctionInfo(__func__);
|
---|
549 | ifstream input(file);
|
---|
550 | string token;
|
---|
551 | stringstream item;
|
---|
552 | double IonMass, ConstrainedPotentialEnergy, ActualTemp;
|
---|
553 | Vector Velocity;
|
---|
554 | ForceMatrix Force;
|
---|
555 |
|
---|
556 | CountElements(); // make sure ElementsInMolecule is up to date
|
---|
557 |
|
---|
558 | // check file
|
---|
559 | if (input == NULL) {
|
---|
560 | return false;
|
---|
561 | } else {
|
---|
562 | // parse file into ForceMatrix
|
---|
563 | if (!Force.ParseMatrix(file, 0,0,0)) {
|
---|
564 | DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl);
|
---|
565 | performCriticalExit();
|
---|
566 | return false;
|
---|
567 | }
|
---|
568 | if (Force.RowCounter[0] != getAtomCount()) {
|
---|
569 | DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
|
---|
570 | performCriticalExit();
|
---|
571 | return false;
|
---|
572 | }
|
---|
573 | // correct Forces
|
---|
574 | Velocity.Zero();
|
---|
575 | for(int i=0;i<getAtomCount();i++)
|
---|
576 | for(int d=0;d<NDIM;d++) {
|
---|
577 | Velocity[d] += Force.Matrix[0][i][d+5];
|
---|
578 | }
|
---|
579 | for(int i=0;i<getAtomCount();i++)
|
---|
580 | for(int d=0;d<NDIM;d++) {
|
---|
581 | Force.Matrix[0][i][d+5] -= Velocity[d]/static_cast<double>(getAtomCount());
|
---|
582 | }
|
---|
583 | // solve a constrained potential if we are meant to
|
---|
584 | if (configuration.DoConstrainedMD) {
|
---|
585 | // calculate forces and potential
|
---|
586 | atom **PermutationMap = NULL;
|
---|
587 | ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
|
---|
588 | EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force);
|
---|
589 | delete[](PermutationMap);
|
---|
590 | }
|
---|
591 |
|
---|
592 | // and perform Verlet integration for each atom with position, velocity and force vector
|
---|
593 | // check size of vectors
|
---|
594 | //ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
|
---|
595 |
|
---|
596 | ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps+1, &configuration, &Force);
|
---|
597 | }
|
---|
598 | // correct velocities (rather momenta) so that center of mass remains motionless
|
---|
599 | Velocity.Zero();
|
---|
600 | IonMass = 0.;
|
---|
601 | ActOnAllAtoms ( &atom::SumUpKineticEnergy, MDSteps+1, &IonMass, &Velocity );
|
---|
602 |
|
---|
603 | // correct velocities (rather momenta) so that center of mass remains motionless
|
---|
604 | Velocity.Scale(1./IonMass);
|
---|
605 | ActualTemp = 0.;
|
---|
606 | ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps+1, &Velocity );
|
---|
607 | Thermostats(configuration, ActualTemp, Berendsen);
|
---|
608 | MDSteps++;
|
---|
609 |
|
---|
610 | // exit
|
---|
611 | return true;
|
---|
612 | };
|
---|
613 |
|
---|
614 | /** Implementation of various thermostats.
|
---|
615 | * All these thermostats apply an additional force which has the following forms:
|
---|
616 | * -# Woodcock
|
---|
617 | * \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
|
---|
618 | * -# Gaussian
|
---|
619 | * \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$
|
---|
620 | * -# Langevin
|
---|
621 | * \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$
|
---|
622 | * -# Berendsen
|
---|
623 | * \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$
|
---|
624 | * -# Nose-Hoover
|
---|
625 | * \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$
|
---|
626 | * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
|
---|
627 | * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
|
---|
628 | * belatedly and the constraint force set.
|
---|
629 | * \param *P Problem at hand
|
---|
630 | * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
|
---|
631 | * \sa InitThermostat()
|
---|
632 | */
|
---|
633 | void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
|
---|
634 | {
|
---|
635 | double ekin = 0.;
|
---|
636 | double E = 0., G = 0.;
|
---|
637 | double delta_alpha = 0.;
|
---|
638 | double ScaleTempFactor;
|
---|
639 | gsl_rng * r;
|
---|
640 | const gsl_rng_type * T;
|
---|
641 |
|
---|
642 | // calculate scale configuration
|
---|
643 | ScaleTempFactor = configuration.TargetTemp/ActualTemp;
|
---|
644 |
|
---|
645 | // differentating between the various thermostats
|
---|
646 | switch(Thermostat) {
|
---|
647 | case None:
|
---|
648 | DoLog(2) && (Log() << Verbose(2) << "Applying no thermostat..." << endl);
|
---|
649 | break;
|
---|
650 | case Woodcock:
|
---|
651 | if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
|
---|
652 | DoLog(2) && (Log() << Verbose(2) << "Applying Woodcock thermostat..." << endl);
|
---|
653 | ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
|
---|
654 | }
|
---|
655 | break;
|
---|
656 | case Gaussian:
|
---|
657 | DoLog(2) && (Log() << Verbose(2) << "Applying Gaussian thermostat..." << endl);
|
---|
658 | ActOnAllAtoms( &atom::Thermostat_Gaussian_init, MDSteps, &G, &E );
|
---|
659 |
|
---|
660 | DoLog(1) && (Log() << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl);
|
---|
661 | ActOnAllAtoms( &atom::Thermostat_Gaussian_least_constraint, MDSteps, G/E, &ekin, &configuration);
|
---|
662 |
|
---|
663 | break;
|
---|
664 | case Langevin:
|
---|
665 | DoLog(2) && (Log() << Verbose(2) << "Applying Langevin thermostat..." << endl);
|
---|
666 | // init random number generator
|
---|
667 | gsl_rng_env_setup();
|
---|
668 | T = gsl_rng_default;
|
---|
669 | r = gsl_rng_alloc (T);
|
---|
670 | // Go through each ion
|
---|
671 | ActOnAllAtoms( &atom::Thermostat_Langevin, MDSteps, r, &ekin, &configuration );
|
---|
672 | break;
|
---|
673 |
|
---|
674 | case Berendsen:
|
---|
675 | DoLog(2) && (Log() << Verbose(2) << "Applying Berendsen-VanGunsteren thermostat..." << endl);
|
---|
676 | ActOnAllAtoms( &atom::Thermostat_Berendsen, MDSteps, ScaleTempFactor, &ekin, &configuration );
|
---|
677 | break;
|
---|
678 |
|
---|
679 | case NoseHoover:
|
---|
680 | DoLog(2) && (Log() << Verbose(2) << "Applying Nose-Hoover thermostat..." << endl);
|
---|
681 | // dynamically evolve alpha (the additional degree of freedom)
|
---|
682 | delta_alpha = 0.;
|
---|
683 | ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
|
---|
684 | delta_alpha = (delta_alpha - (3.*getAtomCount()+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
|
---|
685 | configuration.alpha += delta_alpha*configuration.Deltat;
|
---|
686 | DoLog(3) && (Log() << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl);
|
---|
687 | // apply updated alpha as additional force
|
---|
688 | ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );
|
---|
689 | break;
|
---|
690 | }
|
---|
691 | DoLog(1) && (Log() << Verbose(1) << "Kinetic energy is " << ekin << "." << endl);
|
---|
692 | };
|
---|