1 /* MolecularDynamics.c */
2 /**********************************************************************************************************
3 Copyright (c) 2002-2013 Abdul-Rahman Allouche. All rights reserved
4
5 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
6 documentation files (the Gabedit), to deal in the Software without restriction, including without limitation
7 the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
8 and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
9
10 The above copyright notice and this permission notice shall be included in all copies or substantial portions
11 of the Software.
12
13 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
14 TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
15 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
16 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
17 DEALINGS IN THE SOFTWARE.
18 ************************************************************************************************************/
19 #include "../../Config.h"
20 #include <stdlib.h>
21 #include <math.h>
22 #include <time.h>
23
24 #include "../Common/Global.h"
25 #include "../Utils/AtomsProp.h"
26 #include "../Geometry/Fragments.h"
27 #include "../Geometry/DrawGeom.h"
28 #include "../Utils/Utils.h"
29 #include "Atom.h"
30 #include "Molecule.h"
31 #include "ForceField.h"
32 #include "MolecularDynamics.h"
33
34
35 /*********************************************************************************/
36 static void initMD(MolecularDynamics* molecularDynamics, gdouble temperature, gdouble stepSize, MDIntegratorType integratorType, MDThermostatType thermostat, gdouble friction, gdouble collide, gchar* fileNameTraj, gchar* fileNameProp, gint numberOfRunSteps);
37 static void berendsen(MolecularDynamics* molecularDynamics);
38 static void andersen(MolecularDynamics* molecularDynamics);
39 static void bussi(MolecularDynamics* molecularDynamics);
40 static void rescaleVelocities(MolecularDynamics* molecularDynamics);
41 static void computeEnergies(MolecularDynamics* molecularDynamics);
42 static void applyOneStep(MolecularDynamics* molecularDynamics);
43 static void applyVerlet(MolecularDynamics* molecularDynamics);
44 static void applyBeeman(MolecularDynamics* molecularDynamics);
45 static void applyStochastic(MolecularDynamics* molecularDynamics);
46 static gdouble maxwel(gdouble masse, gdouble temperature);
47 static void newProperties(MolecularDynamics* molecularDynamics, gchar* comments);
48 static void saveProperties(MolecularDynamics* molecularDynamics, gint iStep0, gint iStep, gchar* comments);
49 static void saveTrajectory(MolecularDynamics* molecularDynamics, gint iStep);
50 static void getRandVect(gdouble len, gdouble V[]);
51 static gdouble getEKin(MolecularDynamics* molecularDynamics);
52 static gdouble getKelvin(MolecularDynamics* molecularDynamics);
53 static gdouble drandom();
54 static gdouble normal();
55 /**********************************************************************/
runMolecularDynamicsConfo(MolecularDynamics * molecularDynamics,ForceField * forceField,gint updateFrequency,gdouble heatTime,gdouble equiTime,gdouble runTime,gdouble heatTemperature,gdouble equiTemperature,gdouble runTemperature,gdouble stepSize,MDIntegratorType integratorType,MDThermostatType thermostat,gdouble friction,gdouble collide,gint numberOfGeometries,gchar * fileNameTraj,gchar * fileNameProp)56 ForceField** runMolecularDynamicsConfo(
57 MolecularDynamics* molecularDynamics, ForceField* forceField,
58 gint updateFrequency,
59 gdouble heatTime, gdouble equiTime, gdouble runTime,
60 gdouble heatTemperature, gdouble equiTemperature, gdouble runTemperature,
61 gdouble stepSize,
62 MDIntegratorType integratorType,
63 MDThermostatType thermostat,
64 gdouble friction,
65 gdouble collide,
66 gint numberOfGeometries,
67 gchar* fileNameTraj,
68 gchar* fileNameProp
69 )
70 {
71 gint i;
72 gint j;
73 gchar* str = NULL;
74 gdouble gradientNorm = 0;
75 gint numberOfHeatSteps = 0;
76 gint numberOfEquiSteps = 0;
77 gint numberOfRunSteps = 0;
78 gdouble currentTemp;
79 gint updateNumber = 0;
80 gint n0 = 0;
81 ForceField** geometries = NULL;
82 gint iSel = 0;
83 gint stepSel = 1;
84 /*
85 * physical constants in SI units
86 * ------------------------------
87 * Kb = 1.380662 E-23 J/K
88 * Na = 6.022045 E23 1/mol
89 * e = 1.6021892 E-19 C
90 * eps = 8.85418782 E-12 F/m
91 *
92 * 1 Kcal = 4184.0 J
93 * 1 amu = 1.6605655 E-27 Kg
94 * 1 A = 1.0 E-10 m
95 *
96 * Internally, AKMA units are used:
97 *
98 * timeFactor = SQRT ( ( 1A )**2 * 1amu * Na / 1Kcal )
99 * kBoltzmann = Na *Kb / 1 Kcal
100 */
101
102 /* printf("basname = %s\n",g_path_get_basename(fileNameTraj));*/
103
104 if(forceField->molecule.nAtoms<1) return NULL;
105 if(numberOfGeometries<2) return NULL;
106 geometries = g_malloc(numberOfGeometries*sizeof(ForceField*));
107 for(i=0;i<numberOfGeometries;i++) geometries[i] = NULL;
108
109 molecularDynamics->forceField = forceField;
110 molecularDynamics->numberOfAtoms = forceField->molecule.nAtoms;
111 molecularDynamics->updateFrequency = updateFrequency;
112
113 currentTemp = heatTemperature/2;
114
115 numberOfHeatSteps = heatTime/stepSize*1000;
116 numberOfEquiSteps = equiTime/stepSize*1000;;
117 numberOfRunSteps = runTime/stepSize*1000;;
118
119
120 currentTemp = heatTemperature;
121 if(numberOfHeatSteps==0) currentTemp = equiTemperature;
122 if(numberOfHeatSteps==0 && numberOfEquiSteps==0 ) currentTemp = runTemperature;
123
124 initMD(molecularDynamics,currentTemp,stepSize, integratorType, thermostat, friction, collide, fileNameTraj, fileNameProp, numberOfRunSteps);
125 molecularDynamics->forceField->klass->calculateGradient(molecularDynamics->forceField);
126 computeEnergies(molecularDynamics);
127
128 iSel = -1;
129 {
130 if(str) g_free(str);
131 str = g_strdup_printf("Geometry selected Potential energy = %0.4f", molecularDynamics->potentialEnergy);
132 redrawMolecule(&molecularDynamics->forceField->molecule,str);
133 iSel++;
134 geometries[iSel] = g_malloc(sizeof(ForceField));
135 *geometries[iSel] = copyForceField(molecularDynamics->forceField);
136 Waiting(1);
137 }
138
139 molecularDynamics->temperature = heatTemperature;
140 rescaleVelocities(molecularDynamics);
141
142 currentTemp = heatTemperature;
143 n0 = 0;
144 newProperties(molecularDynamics," ");
145 /*newProperties(molecularDynamics," ----> Heating");*/
146 for (i = 0; i < numberOfHeatSteps; i++ )
147 {
148 molecularDynamics->temperature = currentTemp;
149 applyOneStep(molecularDynamics);
150 currentTemp = heatTemperature + ( runTemperature - heatTemperature ) *
151 ( ( gdouble )( i + 1 )/ numberOfHeatSteps );
152 molecularDynamics->temperature = currentTemp;
153 rescaleVelocities(molecularDynamics);
154 if(StopCalcul) break;
155 if (++updateNumber >= molecularDynamics->updateFrequency )
156 {
157 if(str) g_free(str);
158 str = g_strdup_printf(_("MD Heating: %0.2f fs, T = %0.2f K T(t) = %0.2f Kin = %0.4f Pot = %0.4f Tot = %0.4f"),
159 i*stepSize, currentTemp,
160 molecularDynamics->kelvin,
161 molecularDynamics->kineticEnergy,
162 molecularDynamics->potentialEnergy,
163 molecularDynamics->totalEnergy
164 );
165 redrawMolecule(&molecularDynamics->forceField->molecule,str);
166 updateNumber = 0;
167 }
168 saveProperties(molecularDynamics, n0+i+1, i+1," Heating");
169 }
170
171 currentTemp = equiTemperature;
172 molecularDynamics->temperature = currentTemp;
173 rescaleVelocities(molecularDynamics);
174 if(StopCalcul) numberOfEquiSteps =0;
175 if(StopCalcul) numberOfRunSteps =0;
176 updateNumber = molecularDynamics->updateFrequency;
177 n0 += numberOfHeatSteps;
178 /* newProperties(molecularDynamics," ----> Equilibrium");*/
179 for (i = 0; i < numberOfEquiSteps; i++ )
180 {
181 molecularDynamics->temperature = currentTemp;
182 applyOneStep(molecularDynamics);
183 molecularDynamics->temperature = currentTemp;
184 rescaleVelocities(molecularDynamics);
185 if(StopCalcul) break;
186 if (++updateNumber >= molecularDynamics->updateFrequency )
187 {
188 if(str) g_free(str);
189 str = g_strdup_printf(_("MD Equilibrium: %0.2f fs, T = %0.2f K T(t) = %0.2f K Kin = %0.4f Pot = %0.4f Tot = %0.4f"),
190 i*stepSize, currentTemp,
191 molecularDynamics->kelvin,
192 molecularDynamics->kineticEnergy,
193 molecularDynamics->potentialEnergy,
194 molecularDynamics->totalEnergy
195 );
196 redrawMolecule(&molecularDynamics->forceField->molecule,str);
197 updateNumber = 0;
198 }
199 saveProperties(molecularDynamics, n0+i+1, i+1, " Equilibrium");
200 }
201 updateNumber = molecularDynamics->updateFrequency;
202
203 currentTemp = runTemperature;
204 molecularDynamics->temperature = currentTemp;
205 rescaleVelocities(molecularDynamics);
206 if(StopCalcul) numberOfRunSteps =0;
207 updateNumber = molecularDynamics->updateFrequency;
208 n0 += numberOfEquiSteps;
209 /* newProperties(molecularDynamics," ----> Runing");*/
210 if(str) g_free(str);
211 str = g_strdup_printf(_("Geometry selected Potential energy = %0.4f"), molecularDynamics->potentialEnergy);
212 redrawMolecule(&molecularDynamics->forceField->molecule,str);
213 if(numberOfGeometries>2) stepSel = numberOfRunSteps/(numberOfGeometries-1);
214 else stepSel = numberOfRunSteps;
215 /* printf("Isel = %d\n",stepSel);*/
216 for (i = 0; i < numberOfRunSteps; i++ )
217 {
218 molecularDynamics->temperature = currentTemp;
219 applyOneStep(molecularDynamics);
220 if(molecularDynamics->thermostat == ANDERSEN) andersen(molecularDynamics);
221 if(molecularDynamics->thermostat == BERENDSEN) berendsen(molecularDynamics);
222 if(molecularDynamics->thermostat == BUSSI) bussi(molecularDynamics);
223 if(StopCalcul) break;
224 if (++updateNumber >= molecularDynamics->updateFrequency )
225 {
226 if(str) g_free(str);
227 str = g_strdup_printf(_("MD Running: %0.2f fs, T = %0.2f K T(t) = %0.2f K Kin = %0.4f Pot = %0.4f Tot = %0.4f"),
228 i*stepSize, currentTemp,
229 molecularDynamics->kelvin,
230 molecularDynamics->kineticEnergy,
231 molecularDynamics->potentialEnergy,
232 molecularDynamics->totalEnergy
233 );
234 redrawMolecule(&molecularDynamics->forceField->molecule,str);
235 updateNumber = 0;
236 saveTrajectory(molecularDynamics, i+1);
237 }
238 if((i+1)%stepSel==0 && (iSel+1)<numberOfGeometries)
239 {
240 if(str) g_free(str);
241 str = g_strdup_printf(_("Geometry selected Potential energy = %0.4f"), molecularDynamics->potentialEnergy);
242 redrawMolecule(&molecularDynamics->forceField->molecule,str);
243 iSel++;
244 geometries[iSel] = g_malloc(sizeof(ForceField));
245 *geometries[iSel] = copyForceField(molecularDynamics->forceField);
246 Waiting(1);
247 }
248 saveProperties(molecularDynamics, n0+i+1, i+1," Running");
249 }
250 if(iSel<numberOfGeometries-1)
251 {
252 if(str) g_free(str);
253 str = g_strdup_printf(_("Geometry selected Potential energy = %0.4f"), molecularDynamics->potentialEnergy);
254 redrawMolecule(&molecularDynamics->forceField->molecule,str);
255 iSel++;
256 geometries[iSel] = g_malloc(sizeof(ForceField));
257 *geometries[iSel] = copyForceField(molecularDynamics->forceField);
258 Waiting(1);
259 }
260
261 updateNumber = molecularDynamics->updateFrequency;
262 n0 += numberOfRunSteps;
263
264 molecularDynamics->forceField->klass->calculateGradient(molecularDynamics->forceField);
265 gradientNorm = 0;
266 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
267 for ( j = 0; j < 3; j++)
268 gradientNorm +=
269 molecularDynamics->forceField->molecule.gradient[j][i] *
270 molecularDynamics->forceField->molecule.gradient[j][i];
271
272 gradientNorm = sqrt( gradientNorm );
273 if(str) g_free(str);
274 str = g_strdup_printf(_("End of MD Simulation. Gradient = %f Ekin = %f (Kcal/mol) EPot = %0.4f ETot = %0.4f T(t) = %0.2f"),
275 (gdouble)gradientNorm,
276 molecularDynamics->kineticEnergy,
277 molecularDynamics->potentialEnergy,
278 molecularDynamics->totalEnergy,
279 molecularDynamics->kelvin
280 );
281 redrawMolecule(&molecularDynamics->forceField->molecule,str);
282 g_free(str);
283 if(molecularDynamics->fileTraj)fclose(molecularDynamics->fileTraj);
284 if(molecularDynamics->fileProp)fclose(molecularDynamics->fileProp);
285 freeMolecularDynamics(molecularDynamics);
286 return geometries;
287 }
288 /**********************************************************************/
runMolecularDynamics(MolecularDynamics * molecularDynamics,ForceField * forceField,gint updateFrequency,gdouble heatTime,gdouble equiTime,gdouble runTime,gdouble coolTime,gdouble heatTemperature,gdouble equiTemperature,gdouble runTemperature,gdouble coolTemperature,gdouble stepSize,MDIntegratorType integratorType,MDThermostatType thermostat,gdouble friction,gdouble collide,gchar * fileNameTraj,gchar * fileNameProp)289 void runMolecularDynamics(
290 MolecularDynamics* molecularDynamics, ForceField* forceField,
291 gint updateFrequency,
292 gdouble heatTime, gdouble equiTime, gdouble runTime, gdouble coolTime,
293 gdouble heatTemperature, gdouble equiTemperature, gdouble runTemperature, gdouble coolTemperature,
294 gdouble stepSize,
295 MDIntegratorType integratorType,
296 MDThermostatType thermostat,
297 gdouble friction,
298 gdouble collide,
299 gchar* fileNameTraj,
300 gchar* fileNameProp
301 )
302 {
303 gint i;
304 gint j;
305 gchar* str = NULL;
306 gdouble gradientNorm = 0;
307 gint numberOfHeatSteps = 0;
308 gint numberOfEquiSteps = 0;
309 gint numberOfRunSteps = 0;
310 gint numberOfCoolSteps = 0;
311 gdouble currentTemp;
312 gint updateNumber = 0;
313 gint n0 = 0;
314 /*
315 * physical constants in SI units
316 * ------------------------------
317 * Kb = 1.380662 E-23 J/K
318 * Na = 6.022045 E23 1/mol
319 * e = 1.6021892 E-19 C
320 * eps = 8.85418782 E-12 F/m
321 *
322 * 1 Kcal = 4184.0 J
323 * 1 amu = 1.6605655 E-27 Kg
324 * 1 A = 1.0 E-10 m
325 *
326 * Internally, AKMA units are used:
327 *
328 * timeFactor = SQRT ( ( 1A )**2 * 1amu * Na / 1Kcal )
329 * kBoltzmann = Na *Kb / 1 Kcal
330 */
331
332 /* printf("basname = %s\n",g_path_get_basename(fileNameTraj));*/
333
334 if(forceField->molecule.nAtoms<1) return;
335
336 molecularDynamics->forceField = forceField;
337 molecularDynamics->numberOfAtoms = forceField->molecule.nAtoms;
338 molecularDynamics->updateFrequency = updateFrequency;
339
340 currentTemp = heatTemperature/2;
341
342 numberOfHeatSteps = heatTime/stepSize*1000;
343 numberOfEquiSteps = equiTime/stepSize*1000;;
344 numberOfRunSteps = runTime/stepSize*1000;;
345 numberOfCoolSteps = coolTime/stepSize*1000;;
346
347
348 currentTemp = heatTemperature;
349 if(numberOfHeatSteps==0) currentTemp = equiTemperature;
350 if(numberOfHeatSteps==0 && numberOfEquiSteps==0 ) currentTemp = runTemperature;
351 if(numberOfHeatSteps==0 && numberOfEquiSteps==0 && numberOfRunSteps==0 ) currentTemp = coolTemperature;
352
353 initMD(molecularDynamics,currentTemp,stepSize, integratorType, thermostat, friction, collide, fileNameTraj, fileNameProp, numberOfRunSteps);
354 molecularDynamics->forceField->klass->calculateGradient(molecularDynamics->forceField);
355
356 molecularDynamics->temperature = heatTemperature;
357 rescaleVelocities(molecularDynamics);
358
359 currentTemp = heatTemperature;
360 n0 = 0;
361 newProperties(molecularDynamics," ");
362 /*newProperties(molecularDynamics," ----> Heating");*/
363 for (i = 0; i < numberOfHeatSteps; i++ )
364 {
365 molecularDynamics->temperature = currentTemp;
366 applyOneStep(molecularDynamics);
367 currentTemp = heatTemperature + ( runTemperature - heatTemperature ) *
368 ( ( gdouble )( i + 1 )/ numberOfHeatSteps );
369 molecularDynamics->temperature = currentTemp;
370 rescaleVelocities(molecularDynamics);
371 if(StopCalcul) break;
372 if (++updateNumber >= molecularDynamics->updateFrequency )
373 {
374 if(str) g_free(str);
375 str = g_strdup_printf(_("MD Heating: %0.2f fs, T = %0.2f K T(t) = %0.2f Kin = %0.4f Pot = %0.4f Tot = %0.4f"),
376 i*stepSize, currentTemp,
377 molecularDynamics->kelvin,
378 molecularDynamics->kineticEnergy,
379 molecularDynamics->potentialEnergy,
380 molecularDynamics->totalEnergy
381 );
382 redrawMolecule(&molecularDynamics->forceField->molecule,str);
383 updateNumber = 0;
384 }
385 saveProperties(molecularDynamics, n0+i+1, i+1," Heating");
386 }
387
388 currentTemp = equiTemperature;
389 molecularDynamics->temperature = currentTemp;
390 rescaleVelocities(molecularDynamics);
391 if(StopCalcul) numberOfEquiSteps =0;
392 if(StopCalcul) numberOfRunSteps =0;
393 if(StopCalcul) numberOfCoolSteps =0;
394 updateNumber = molecularDynamics->updateFrequency;
395 n0 += numberOfHeatSteps;
396 /* newProperties(molecularDynamics," ----> Equilibrium");*/
397 for (i = 0; i < numberOfEquiSteps; i++ )
398 {
399 molecularDynamics->temperature = currentTemp;
400 applyOneStep(molecularDynamics);
401 molecularDynamics->temperature = currentTemp;
402 rescaleVelocities(molecularDynamics);
403 if(StopCalcul) break;
404 if (++updateNumber >= molecularDynamics->updateFrequency )
405 {
406 if(str) g_free(str);
407 str = g_strdup_printf(_("MD Equilibrium: %0.2f fs, T = %0.2f K T(t) = %0.2f K Kin = %0.4f Pot = %0.4f Tot = %0.4f"),
408 i*stepSize, currentTemp,
409 molecularDynamics->kelvin,
410 molecularDynamics->kineticEnergy,
411 molecularDynamics->potentialEnergy,
412 molecularDynamics->totalEnergy
413 );
414 redrawMolecule(&molecularDynamics->forceField->molecule,str);
415 updateNumber = 0;
416 }
417 saveProperties(molecularDynamics, n0+i+1, i+1, " Equilibrium");
418 }
419 updateNumber = molecularDynamics->updateFrequency;
420
421 currentTemp = runTemperature;
422 molecularDynamics->temperature = currentTemp;
423 rescaleVelocities(molecularDynamics);
424 if(StopCalcul) numberOfRunSteps =0;
425 if(StopCalcul) numberOfCoolSteps =0;
426 updateNumber = molecularDynamics->updateFrequency;
427 n0 += numberOfEquiSteps;
428 /* newProperties(molecularDynamics," ----> Runing");*/
429 for (i = 0; i < numberOfRunSteps; i++ )
430 {
431 molecularDynamics->temperature = currentTemp;
432 applyOneStep(molecularDynamics);
433 if(molecularDynamics->thermostat == ANDERSEN) andersen(molecularDynamics);
434 if(molecularDynamics->thermostat == BERENDSEN) berendsen(molecularDynamics);
435 if(molecularDynamics->thermostat == BUSSI) bussi(molecularDynamics);
436 if(StopCalcul) break;
437 if (++updateNumber >= molecularDynamics->updateFrequency )
438 {
439 if(str) g_free(str);
440 str = g_strdup_printf(_("MD Running: %0.2f fs, T = %0.2f K T(t) = %0.2f K Kin = %0.4f Pot = %0.4f Tot = %0.4f"),
441 i*stepSize, currentTemp,
442 molecularDynamics->kelvin,
443 molecularDynamics->kineticEnergy,
444 molecularDynamics->potentialEnergy,
445 molecularDynamics->totalEnergy
446 );
447 redrawMolecule(&molecularDynamics->forceField->molecule,str);
448 updateNumber = 0;
449 saveTrajectory(molecularDynamics, i+1);
450 }
451 saveProperties(molecularDynamics, n0+i+1, i+1," Running");
452 }
453
454 updateNumber = molecularDynamics->updateFrequency;
455 if(StopCalcul) numberOfCoolSteps =0;
456 n0 += numberOfRunSteps;
457 /* newProperties(molecularDynamics," ----> Cooling");*/
458 for (i = 0; i < numberOfCoolSteps; i++ )
459 {
460 currentTemp = runTemperature - ( runTemperature - coolTemperature ) *
461 ( ( gdouble )( i + 1 )/ numberOfCoolSteps );
462 molecularDynamics->temperature = currentTemp;
463 rescaleVelocities(molecularDynamics);
464 molecularDynamics->temperature = currentTemp;
465 applyOneStep(molecularDynamics);
466 if(StopCalcul) break;
467 if (++updateNumber >= molecularDynamics->updateFrequency )
468 {
469 if(str) g_free(str);
470 str = g_strdup_printf(_("MD Cooling: %0.2f fs, T = %0.2f K T(t) = %0.2f K Kin = %0.4f Pot = %0.4f Tot = %0.4f"),
471 i*stepSize, currentTemp,
472 molecularDynamics->kelvin,
473 molecularDynamics->kineticEnergy,
474 molecularDynamics->potentialEnergy,
475 molecularDynamics->totalEnergy
476 );
477 redrawMolecule(&molecularDynamics->forceField->molecule,str);
478 updateNumber = 0;
479 }
480 saveProperties(molecularDynamics, n0+i+1, i+1," Cooling");
481 }
482 molecularDynamics->forceField->klass->calculateGradient(molecularDynamics->forceField);
483 gradientNorm = 0;
484 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
485 for ( j = 0; j < 3; j++)
486 gradientNorm +=
487 molecularDynamics->forceField->molecule.gradient[j][i] *
488 molecularDynamics->forceField->molecule.gradient[j][i];
489
490 gradientNorm = sqrt( gradientNorm );
491 if(str) g_free(str);
492 str = g_strdup_printf(_("End of MD Simulation. Gradient = %f Ekin = %f (Kcal/mol) EPot = %0.4f ETot = %0.4f T(t) = %0.2f"),
493 (gdouble)gradientNorm,
494 molecularDynamics->kineticEnergy,
495 molecularDynamics->potentialEnergy,
496 molecularDynamics->totalEnergy,
497 molecularDynamics->kelvin
498 );
499 redrawMolecule(&molecularDynamics->forceField->molecule,str);
500 g_free(str);
501 if(molecularDynamics->fileTraj)fclose(molecularDynamics->fileTraj);
502 if(molecularDynamics->fileProp)fclose(molecularDynamics->fileProp);
503 freeMolecularDynamics(molecularDynamics);
504 }
505 /*********************************************************************************/
removeTranslation(MolecularDynamics * molecularDynamics)506 static void removeTranslation(MolecularDynamics* molecularDynamics)
507 {
508 gdouble vtot[3] = {0,0,0};
509 gint i;
510 gint j;
511 gdouble mass = 1.0;
512 gdouble totMass = 0.0;
513 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
514 {
515 mass = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
516 totMass += mass;
517 for ( j = 0; j < 3; j++)
518 {
519 vtot[j] += mass*molecularDynamics->velocity[i][j];
520 }
521 }
522
523 for ( j = 0; j < 3; j++)
524 vtot[j] /= totMass;
525
526 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
527 for ( j = 0; j < 3; j++)
528 molecularDynamics->velocity[i][j] -= vtot[j];
529 /* check */
530 /*
531 for ( j = 0; j < 3; j++)
532 vtot[j] = 0;
533 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
534 {
535 mass = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
536 for ( j = 0; j < 3; j++)
537 {
538 vtot[j] += mass*molecularDynamics->velocity[i][j];
539 }
540 }
541 printf("Trans velocity = %f %f %f\n",vtot[0], vtot[1], vtot[2]);
542 */
543 }
544 /**************************************************/
InverseTensor(gdouble mat[3][3],gdouble invmat[3][3])545 static gboolean InverseTensor(gdouble mat[3][3],gdouble invmat[3][3])
546 {
547 gdouble t4,t6,t8,t10,t12,t14,t17;
548 gdouble d = 0;
549
550 t4 = mat[0][0]*mat[1][1];
551 t6 = mat[0][0]*mat[1][2];
552 t8 = mat[0][1]*mat[1][0];
553 t10 = mat[0][2]*mat[1][0];
554 t12 = mat[0][1]*mat[2][0];
555 t14 = mat[0][2]*mat[2][0];
556 d =(t4*mat[2][2]-t6*mat[2][1]-t8*mat[2][2]+t10*mat[2][1]+t12*mat[1][2]-t14*mat[1][1]);
557 if(d == 0)
558 {
559 invmat[0][0] = 0;
560 invmat[0][1] = 0;
561 invmat[0][2] = 0;
562 invmat[1][0] = 0;
563 invmat[1][1] = 0;
564 invmat[1][2] = 0;
565 invmat[2][0] = 0;
566 invmat[2][1] = 0;
567 invmat[2][2] = 0;
568 return FALSE;
569 }
570 t17 = 1/d;
571 invmat[0][0] = (mat[1][1]*mat[2][2]-mat[1][2]*mat[2][1])*t17;
572 invmat[0][1] = -(mat[0][1]*mat[2][2]-mat[0][2]*mat[2][1])*t17;
573 invmat[0][2] = -(-mat[0][1]*mat[1][2]+mat[0][2]*mat[1][1])*t17;
574 invmat[1][0] = -(mat[1][0]*mat[2][2]-mat[1][2]*mat[2][0])*t17;
575 invmat[1][1] = (mat[0][0]*mat[2][2]-t14)*t17;
576 invmat[1][2] = -(t6-t10)*t17;
577 invmat[2][0] = -(-mat[1][0]*mat[2][1]+mat[1][1]*mat[2][0])*t17;
578 invmat[2][1] = -(mat[0][0]*mat[2][1]-t12)*t17;
579 invmat[2][2] = (t4-t8)*t17;
580
581 return TRUE;
582 }
583 /*********************************************************************************/
removeRotation(MolecularDynamics * molecularDynamics)584 static void removeRotation(MolecularDynamics* molecularDynamics)
585 {
586 gdouble vtot[3] = {0,0,0};
587 gdouble cm[3] = {0,0,0};
588 gdouble L[3] = {0,0,0};
589 gint i;
590 gint j;
591 gint k;
592 gdouble mass = 1.0;
593 gdouble totMass = 0.0;
594 gdouble cdel[3];
595 gdouble vAng[3]={0,0,0};
596 gdouble tensor[3][3];
597 gdouble invTensor[3][3];
598 gdouble xx, xy,xz,yy,yz,zz;
599 /* find the center of mass coordinates and total velocity*/
600 AtomMol* atoms = molecularDynamics->forceField->molecule.atoms;
601
602
603 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
604 {
605 mass = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
606 totMass += mass;
607 for ( j = 0; j < 3; j++)
608 cm[j] += mass*atoms[i].coordinates[j];
609 for ( j = 0; j < 3; j++)
610 vtot[j] += mass*molecularDynamics->velocity[i][j];
611 }
612
613
614 for ( j = 0; j < 3; j++)
615 cm[j] /= totMass;
616 for ( j = 0; j < 3; j++)
617 vtot[j] /= totMass;
618
619 /* compute the angular momentum */
620 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
621 {
622 mass = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
623 for ( j = 0; j < 3; j++)
624 L[j] += (
625 atoms[i].coordinates[(j+1)%3]*molecularDynamics->velocity[i][(j+2)%3]
626 - atoms[i].coordinates[(j+2)%3]*molecularDynamics->velocity[i][(j+1)%3]
627 )*mass;
628 }
629 for ( j = 0; j < 3; j++)
630 L[j] -= (
631 cm[(j+1)%3]*vtot[(j+2)%3]
632 - cm[(j+2)%3]*vtot[(j+1)%3]
633 )*totMass;
634
635 /* calculate and invert the inertia tensor */
636 for ( k = 0; k < 3; k++)
637 for ( j = 0; j < 3; j++)
638 tensor[k][j] = 0;
639 xx = 0;
640 yy = 0;
641 zz = 0;
642 xy = 0;
643 xz = 0;
644 yz = 0;
645 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
646 {
647 mass = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
648 for ( j = 0; j < 3; j++)
649 cdel[j] = atoms[i].coordinates[j]-cm[j];
650 xx += cdel[0]*cdel[0]*mass;
651 xy += cdel[0]*cdel[1]*mass;
652 xz += cdel[0]*cdel[2]*mass;
653 yy += cdel[1]*cdel[1]*mass;
654 yz += cdel[1]*cdel[2]*mass;
655 zz += cdel[2]*cdel[2]*mass;
656 }
657 tensor[0][0] = yy+zz;
658 tensor[1][0] = -xy;
659 tensor[2][0] = -xz;
660 tensor[0][1] = -xy;
661 tensor[1][1] = xx+zz;
662 tensor[2][1] = -yz;
663 tensor[0][2] = -xz;
664 tensor[1][2] = -yz;
665 tensor[2][2] = xx+yy;
666 if(InverseTensor(tensor,invTensor))
667 {
668 for ( j = 0; j < 3; j++)
669 {
670 vAng[j] = 0;
671 for ( k = 0; k < 3; k++)
672 vAng[j] += invTensor[j][k]*L[k];
673 }
674 /* eliminate any rotation about the system center of mass */
675 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
676 {
677 for ( j = 0; j < 3; j++)
678 cdel[j] = atoms[i].coordinates[j]-cm[j];
679 for ( j = 0; j < 3; j++)
680 molecularDynamics->velocity[i][j] +=
681 cdel[(j+1)%3]*vAng[(j+2)%3]-
682 cdel[(j+2)%3]*vAng[(j+1)%3];
683 }
684 }
685 /* check */
686 /*
687 for ( j = 0; j < 3; j++)
688 L[j] = 0;
689 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
690 {
691 mass = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
692 for ( j = 0; j < 3; j++)
693 L[j] += (
694 atoms[i].coordinates[(j+1)%3]*molecularDynamics->velocity[i][(j+2)%3]
695 - atoms[i].coordinates[(j+2)%3]*molecularDynamics->velocity[i][(j+1)%3]
696 )*mass;
697 }
698 for ( j = 0; j < 3; j++)
699 {
700 vAng[j] = 0;
701 for ( k = 0; k < 3; k++)
702 vAng[j] += invTensor[j][k]*L[k];
703 }
704 printf("Angular velocity = %f %f %f\n",vAng[0], vAng[1], vAng[2]);
705 */
706
707 }
708 /*********************************************************************************/
removeTranslationAndRotation(MolecularDynamics * molecularDynamics)709 static void removeTranslationAndRotation(MolecularDynamics* molecularDynamics)
710 {
711 removeTranslation(molecularDynamics);
712 removeRotation(molecularDynamics);
713 }
714 /*********************************************************************************/
initSD(MolecularDynamics * molecularDynamics,gdouble friction)715 static void initSD(MolecularDynamics* molecularDynamics, gdouble friction)
716 {
717 /* gdouble fsInAKMA = 1.0/sqrt(1e-10*1e-10*1.6605655e-27*6.022045e23/4184.0)/1e15;*/
718 static gdouble fsInAKMA = 0.020454828110640;
719 gint i;
720
721
722 molecularDynamics->friction = friction/(fsInAKMA)/1000;
723
724 molecularDynamics->positionFriction = NULL;
725 molecularDynamics->velocityFriction = NULL;
726 molecularDynamics->accelarationFriction = NULL;
727 molecularDynamics->gamma = NULL;
728 molecularDynamics->positionRandom = NULL;
729 molecularDynamics->velocityRandom = NULL;
730
731 if(molecularDynamics->integratorType != STOCHASTIC) return;
732
733 molecularDynamics->positionFriction = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble));
734 molecularDynamics->velocityFriction = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble));
735 molecularDynamics->accelarationFriction = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble));
736 molecularDynamics->gamma = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble));
737
738 molecularDynamics->positionRandom = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble*));
739 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
740 molecularDynamics->positionRandom[i] = g_malloc(3*sizeof(gdouble));
741
742 molecularDynamics->velocityRandom = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble*));
743 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
744 molecularDynamics->velocityRandom[i] = g_malloc(3*sizeof(gdouble));
745
746 }
747 /*********************************************************************************/
initMD(MolecularDynamics * molecularDynamics,gdouble temperature,gdouble stepSize,MDIntegratorType integratorType,MDThermostatType thermostat,gdouble friction,gdouble collide,gchar * fileNameTraj,gchar * fileNameProp,gint numberOfRunSteps)748 static void initMD(MolecularDynamics* molecularDynamics, gdouble temperature, gdouble stepSize, MDIntegratorType integratorType, MDThermostatType thermostat, gdouble friction, gdouble collide, gchar* fileNameTraj, gchar* fileNameProp, gint numberOfRunSteps)
749 {
750 gint i;
751 gint j;
752 /* gdouble fsInAKMA = 1.0/sqrt(1e-10*1e-10*1.6605655e-27*6.022045e23/4184.0)/1e15;*/
753 static gdouble fsInAKMA = 0.020454828110640;
754 gdouble dt = stepSize * fsInAKMA;
755
756 molecularDynamics->collide = collide;
757 molecularDynamics->potentialEnergy = 0;
758 molecularDynamics->kineticEnergy = 0;
759 molecularDynamics->totalEnergy = 0;
760 molecularDynamics->kelvin = 0;
761 molecularDynamics->temperature = temperature;
762 molecularDynamics->thermostat = NONE;
763
764 molecularDynamics->integratorType = integratorType;
765 molecularDynamics->thermostat = thermostat;
766 molecularDynamics->fileTraj = NULL;
767 molecularDynamics->fileProp = NULL;
768
769 molecularDynamics->velocity = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble*));
770 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
771 molecularDynamics->velocity[i] = g_malloc(3*sizeof(gdouble));
772
773 molecularDynamics->a = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble*));
774 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
775 molecularDynamics->a[i] = g_malloc(3*sizeof(gdouble));
776
777 molecularDynamics->aold = NULL;
778 if(molecularDynamics->integratorType==BEEMAN)
779 {
780 molecularDynamics->aold = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble*));
781 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
782 molecularDynamics->aold[i] = g_malloc(3*sizeof(gdouble));
783 }
784 molecularDynamics->coordinatesOld = NULL;
785 molecularDynamics->moved = NULL;
786 molecularDynamics->update = NULL;
787 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS)
788 {
789 molecularDynamics->coordinatesOld = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gdouble*));
790 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
791 molecularDynamics->coordinatesOld[i] = g_malloc(3*sizeof(gdouble));
792 molecularDynamics->moved = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gboolean));
793 molecularDynamics->update = g_malloc(molecularDynamics->numberOfAtoms *sizeof(gboolean));
794
795 }
796
797 if(fileNameTraj)
798 {
799 molecularDynamics->fileTraj = fopen(fileNameTraj, "w");
800 if(molecularDynamics->fileTraj != NULL)
801 {
802 fprintf(molecularDynamics->fileTraj,"[Gabedit Format]\n");
803 fprintf(molecularDynamics->fileTraj,"\n");
804 fprintf(molecularDynamics->fileTraj,"[MD]\n");
805 if(molecularDynamics->updateFrequency>0) numberOfRunSteps/=molecularDynamics->updateFrequency;
806 fprintf(molecularDynamics->fileTraj," %d\n",numberOfRunSteps);
807 }
808 }
809 if(fileNameProp)
810 {
811 molecularDynamics->fileProp = fopen(fileNameProp, "w");
812 }
813
814 srand ( (unsigned)time (NULL));
815
816 molecularDynamics->dt = dt;
817 molecularDynamics->dt_2 = dt/2.0;
818 molecularDynamics->dt2_2 = dt*dt/2;;
819 molecularDynamics->dt_8 = dt/8.0;
820 molecularDynamics->dt2_8 = dt*dt/8.0;
821
822 initSD(molecularDynamics, friction);
823
824
825 molecularDynamics->forceField->klass->calculateGradient(molecularDynamics->forceField);
826 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
827 {
828 gdouble m = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
829 for ( j = 0; j < 3; j++) molecularDynamics->a[i][j] = -molecularDynamics->forceField->molecule.gradient[j][i]/m;
830 if(molecularDynamics->aold)
831 for ( j = 0; j < 3; j++) molecularDynamics->aold[i][j] = molecularDynamics->a[i][j];
832 if(molecularDynamics->coordinatesOld)
833 for ( j = 0; j < 3; j++) molecularDynamics->coordinatesOld[i][j] = molecularDynamics->forceField->molecule.atoms[i].coordinates[j];
834 if(molecularDynamics->moved) molecularDynamics->moved[i] = FALSE;
835 if(molecularDynamics->update) molecularDynamics->update[i] = FALSE;
836 }
837
838 if(temperature<=0)
839 {
840 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
841 for ( j = 0; j < 3; j++)
842 molecularDynamics->velocity[i][j] = 0.0;
843 }
844 else
845 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
846 {
847 if(molecularDynamics->forceField->molecule.atoms[i].variable)
848 {
849 gdouble speed = maxwel(molecularDynamics->forceField->molecule.atoms[i].prop.masse,temperature);
850 getRandVect(speed, molecularDynamics->velocity[i]);
851 }
852 else
853 for( j = 0; j < 3; j++) molecularDynamics->velocity[i][j] = 0.0;
854 }
855 molecularDynamics->nvariables = 0;
856 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
857 if(molecularDynamics->forceField->molecule.atoms[i].variable) molecularDynamics->nvariables +=1;
858 if(molecularDynamics->nvariables==0)
859 {
860 molecularDynamics->nvariables= molecularDynamics->numberOfAtoms;
861 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
862 {
863 molecularDynamics->forceField->molecule.atoms[i].variable = TRUE;
864 if(temperature>0)
865 {
866 gdouble speed = maxwel(molecularDynamics->forceField->molecule.atoms[i].prop.masse,temperature);
867 getRandVect(speed, molecularDynamics->velocity[i]);
868 }
869 }
870 }
871 molecularDynamics->nfree = 3* molecularDynamics->nvariables-molecularDynamics->forceField->numberOfRattleConstraintsTerms;
872 removeTranslationAndRotation(molecularDynamics);
873 if(molecularDynamics->nvariables==molecularDynamics->numberOfAtoms) molecularDynamics->nfree -=6;
874 if(molecularDynamics->nvariables==molecularDynamics->numberOfAtoms-1) molecularDynamics->nfree -=3;
875 if(molecularDynamics->nvariables==molecularDynamics->numberOfAtoms-2) molecularDynamics->nfree -=1;
876 printf("nfree =%d\n",molecularDynamics->nfree);
877 if( molecularDynamics->nfree<1) {
878 StopCalcul=TRUE;
879 }
880 if( molecularDynamics->nfree<1) molecularDynamics->nfree = 1;
881
882 removeTranslationAndRotation(molecularDynamics);
883 }
884 /*********************************************************************************/
rescaleVelocities(MolecularDynamics * molecularDynamics)885 static void rescaleVelocities(MolecularDynamics* molecularDynamics)
886 {
887 berendsen(molecularDynamics);
888 }
889 /*********************************************************************************/
berendsen(MolecularDynamics * molecularDynamics)890 static void berendsen(MolecularDynamics* molecularDynamics)
891 {
892 gint i;
893 gint j;
894 static gdouble fsInAKMA = 0.020454828110640;
895 gdouble ekin = 0;
896 gdouble kelvin = 0;
897 gint nfree = molecularDynamics->nfree;
898 static gdouble Kb = 1.9871914e-3;
899 gdouble scale = 1.0;
900 gdouble dt = molecularDynamics->dt;
901 gdouble tautemp = 1.0/(molecularDynamics->collide)*1000*fsInAKMA;
902 gdouble masse = 1.0;
903 if(molecularDynamics->temperature<=0) return;
904 if(nfree<1) return;
905 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
906 {
907 masse = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
908 for ( j = 0; j < 3; j++)
909 ekin += molecularDynamics->velocity[i][j]*molecularDynamics->velocity[i][j]*
910 masse;
911 }
912 /*
913 ekin /= 2;
914 kelvin = 2* ekin / ( nfree * Kb);
915 */
916 kelvin = ekin / ( nfree * Kb);
917 /* if(tautemp>dt) tautemp = dt;*/
918 scale = sqrt(1.0 + (dt/tautemp)*(molecularDynamics->temperature/kelvin-1.0));
919 /* printf("temp = %f kelvin = %f scale = %f\n",molecularDynamics->temperature, kelvin, scale);*/
920 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
921 {
922 if(molecularDynamics->forceField->molecule.atoms[i].variable)
923 for ( j = 0; j < 3; j++)
924 molecularDynamics->velocity[i][j] *= scale;
925 }
926 removeTranslationAndRotation(molecularDynamics);
927 }
928 /*********************************************************************************/
andersen(MolecularDynamics * molecularDynamics)929 static void andersen(MolecularDynamics* molecularDynamics)
930 {
931 gint i;
932 /* gdouble fsInAKMA = 1.0/sqrt(1e-10*1e-10*1.6605655e-27*6.022045e23/4184.0)/1e15;*/
933 static gdouble fsInAKMA = 0.020454828110640;
934 gdouble tau = 1.0/molecularDynamics->collide*1000*fsInAKMA; /* in fs */
935 gdouble rate;
936 static gdouble Kb = 1.9871914e-3;
937 if(molecularDynamics->temperature<=0) return;
938 if(molecularDynamics->numberOfAtoms<1) return;
939
940 rate = molecularDynamics->dt / tau;
941 rate /= pow(molecularDynamics->nvariables,2.0/3.0);
942
943 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
944 {
945 gdouble trial = drandom();
946 if(molecularDynamics->forceField->molecule.atoms[i].variable)
947 if(trial<rate)
948 {
949 /*
950 gdouble speed = maxwel(
951 molecularDynamics->forceField->molecule.atoms[i].prop.masse,
952 molecularDynamics->temperature
953 );
954 getRandVect(speed, molecularDynamics->velocity[i]);
955 */
956 double speed = sqrt(Kb* molecularDynamics->temperature/molecularDynamics->forceField->molecule.atoms[i].prop.masse);
957 double pnorm = normal();
958 molecularDynamics->velocity[i][0] = pnorm*speed;
959 pnorm = normal();
960 molecularDynamics->velocity[i][1] = pnorm*speed;
961 pnorm = normal();
962 molecularDynamics->velocity[i][2] = pnorm*speed;
963 }
964 }
965 }
966 /*********************************************************************************/
bussi(MolecularDynamics * molecularDynamics)967 static void bussi(MolecularDynamics* molecularDynamics)
968 {
969 static gdouble fsInAKMA = 0.020454828110640;
970 gint nfree = molecularDynamics->nfree;
971 static gdouble Kb = 1.9871914e-3;
972 gdouble scale = 1.0;
973 gdouble dt = molecularDynamics->dt;
974 gdouble tautemp = 1.0/(molecularDynamics->collide)*1000*fsInAKMA;
975 gdouble c = exp(-dt/tautemp);
976 gdouble ekin = getEKin(molecularDynamics);
977 gdouble kelvin = 2*ekin / ( nfree * Kb);
978 gdouble d = (1.0-c) * (molecularDynamics->temperature/kelvin) / (nfree);
979 gdouble r = normal ();
980 gdouble si = 0.0;
981 gdouble s = 0.0;
982 gint i,j;
983 if(molecularDynamics->temperature<=0) return;
984 if(nfree<1) return;
985 for(i=0;i<nfree-1;i++)
986 {
987 si = normal ();
988 s += si*si;
989 }
990 scale = c + (s+r*r)*d + 2.0*r*sqrt(c*d);
991 scale = sqrt(scale);
992 if (r+sqrt(c/d)<0) scale = -scale;
993 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
994 if(molecularDynamics->forceField->molecule.atoms[i].variable)
995 for ( j = 0; j < 3; j++)
996 molecularDynamics->velocity[i][j] *= scale;
997 removeTranslationAndRotation(molecularDynamics);
998 }
999 /*********************************************************************************/
newAccelaration(MolecularDynamics * molecularDynamics)1000 static void newAccelaration(MolecularDynamics* molecularDynamics)
1001 {
1002 gint i;
1003 gint j;
1004 molecularDynamics->forceField->klass->calculateGradient(molecularDynamics->forceField);
1005 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
1006 {
1007 gdouble m = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
1008 if(molecularDynamics->aold)
1009 for ( j = 0; j < 3; j++)
1010 molecularDynamics->aold[i][j] = molecularDynamics->a[i][j];
1011
1012 for ( j = 0; j < 3; j++)
1013 molecularDynamics->a[i][j] = -molecularDynamics->forceField->molecule.gradient[j][i]/m;
1014 }
1015 }
1016 /*********************************************************************************/
computeEnergies(MolecularDynamics * molecularDynamics)1017 static void computeEnergies(MolecularDynamics* molecularDynamics)
1018 {
1019 molecularDynamics->kineticEnergy = getEKin(molecularDynamics);
1020 molecularDynamics->potentialEnergy = molecularDynamics->forceField->klass->calculateEnergyTmp(
1021 molecularDynamics->forceField,&molecularDynamics->forceField->molecule);
1022 molecularDynamics->totalEnergy = molecularDynamics->kineticEnergy + molecularDynamics->potentialEnergy;
1023 molecularDynamics->kelvin = getKelvin(molecularDynamics);
1024 }
1025 /*********************************************************************************/
applyOneStep(MolecularDynamics * molecularDynamics)1026 static void applyOneStep(MolecularDynamics* molecularDynamics)
1027 {
1028 if(molecularDynamics->integratorType==VERLET) applyVerlet(molecularDynamics);
1029 else if(molecularDynamics->integratorType==BEEMAN) applyBeeman(molecularDynamics);
1030 else applyStochastic(molecularDynamics);
1031 computeEnergies(molecularDynamics);
1032 }
1033 /*********************************************************************************/
applyRattleFirstPortion(MolecularDynamics * molecularDynamics)1034 static void applyRattleFirstPortion(MolecularDynamics* molecularDynamics)
1035 {
1036 gint i;
1037 gint k;
1038 gint maxIter = 1000;
1039 gdouble omega = 1.0;
1040 gdouble tolerance = 1e-6;
1041 gboolean done = FALSE;
1042 gint nIter = 0;
1043 gint a1 = 0;
1044 gint a2 = 0;
1045 gdouble r2ij;
1046 gdouble dot;
1047 gdouble invMass1;
1048 gdouble invMass2;
1049 gdouble delta;
1050 gdouble term = 0;
1051 gdouble terms[3];
1052 gdouble d;
1053 Molecule* m = &molecularDynamics->forceField->molecule;
1054 ForceField* forceField = molecularDynamics->forceField;
1055 gdouble deltaMax = 0;
1056
1057 if(forceField->options.rattleConstraints==NOCONSTRAINTS) return;
1058 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1059 {
1060 molecularDynamics->moved[i] = molecularDynamics->forceField->molecule.atoms[i].variable;
1061 molecularDynamics->update[i] = FALSE;
1062 }
1063 /* maxIter *= molecularDynamics->forceField->numberOfRattleConstraintsTerms;*/
1064 do{
1065 nIter++;
1066 done=TRUE;
1067 deltaMax = 0;
1068 for (i = 0; i < molecularDynamics->forceField->numberOfRattleConstraintsTerms; i++)
1069 {
1070 a1 = (gint)molecularDynamics->forceField->rattleConstraintsTerms[0][i];
1071 a2 = (gint)molecularDynamics->forceField->rattleConstraintsTerms[1][i];
1072 if( !molecularDynamics->moved[a1] && !molecularDynamics->moved[a2] ) continue;
1073 r2ij = 0;
1074 for (k=0;k<3;k++)
1075 {
1076 d = m->atoms[a2].coordinates[k]-m->atoms[a1].coordinates[k];
1077 r2ij +=d*d;
1078 }
1079 delta = molecularDynamics->forceField->rattleConstraintsTerms[2][i]-r2ij;
1080 /* if(fabs(delta)<=tolerance) continue;*/
1081 if(r2ij>0 && fabs(delta/r2ij)<=tolerance) continue;
1082 if(deltaMax<fabs(delta)) deltaMax = fabs(delta);
1083 done = FALSE;
1084 molecularDynamics->update[a1] = TRUE;
1085 molecularDynamics->update[a2] = TRUE;
1086 /* here : rattle image for PBC, not yet implemented */
1087 dot = 0;
1088 for (k=0;k<3;k++)
1089 {
1090 d = m->atoms[a2].coordinates[k]-m->atoms[a1].coordinates[k];
1091 dot +=d*(molecularDynamics->coordinatesOld[a2][k]-molecularDynamics->coordinatesOld[a1][k]);
1092 }
1093 invMass1 = 1/m->atoms[a1].prop.masse;
1094 invMass2 = 1/m->atoms[a2].prop.masse;
1095 term = omega*delta / (2.0*(invMass1+invMass2)*dot);
1096 for (k=0;k<3;k++)
1097 {
1098 terms[k] = (molecularDynamics->coordinatesOld[a2][k]-molecularDynamics->coordinatesOld[a1][k])*term;
1099 }
1100 for (k=0;k<3;k++) m->atoms[a1].coordinates[k] -= terms[k]*invMass1;
1101 for (k=0;k<3;k++) m->atoms[a2].coordinates[k] += terms[k]*invMass2;
1102
1103 invMass1 /= molecularDynamics->dt;
1104 invMass2 /= molecularDynamics->dt;
1105 for (k=0;k<3;k++) molecularDynamics->velocity[a1][k] -= terms[k]*invMass1;
1106 for (k=0;k<3;k++) molecularDynamics->velocity[a2][k] += terms[k]*invMass2;
1107 }
1108 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1109 {
1110 molecularDynamics->moved[i] = molecularDynamics->update[i];
1111 molecularDynamics->update[i] = FALSE;
1112 }
1113 }while(!done && nIter<maxIter);
1114 if(nIter>=maxIter && deltaMax>tolerance*10)
1115 {
1116 printf(_("Rattle first portion : Warning, distance constraints not satisfied\n"));
1117 /*
1118 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1119 {
1120 printf("atom#%d\n",i);
1121 printf("Old coord\n");
1122 for (k=0;k<3;k++) printf("%f ",molecularDynamics->coordinatesOld[i][k]);
1123 printf("\n");
1124 printf("New coord\n");
1125 for (k=0;k<3;k++) printf("%f ",m->atoms[i].coordinates[k]);
1126 printf("\n");
1127 }
1128 exit(1);
1129 */
1130 for (i = 0; i < molecularDynamics->forceField->numberOfRattleConstraintsTerms; i++)
1131 {
1132 a1 = (gint)molecularDynamics->forceField->rattleConstraintsTerms[0][i];
1133 a2 = (gint)molecularDynamics->forceField->rattleConstraintsTerms[1][i];
1134 r2ij = 0;
1135 for (k=0;k<3;k++)
1136 {
1137 d = m->atoms[a2].coordinates[k]-m->atoms[a1].coordinates[k];
1138 r2ij +=d*d;
1139 }
1140 delta = molecularDynamics->forceField->rattleConstraintsTerms[2][i]-r2ij;
1141 printf("%d %d %s %s r2ij=%f r2Old=%f delta=%f\n",
1142 a1,a2,
1143 molecularDynamics->forceField->molecule.atoms[a1].mmType,
1144 molecularDynamics->forceField->molecule.atoms[a2].mmType,
1145 r2ij, molecularDynamics->forceField->rattleConstraintsTerms[2][i],delta);
1146 }
1147 StopCalcul=TRUE;
1148 }
1149 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1150 if(!m->atoms[i].variable)
1151 {
1152 for (k=0;k<3;k++) m->atoms[i].coordinates[k] = molecularDynamics->coordinatesOld[i][k];
1153 for (k=0;k<3;k++) molecularDynamics->velocity[i][k] = 0;
1154 }
1155
1156 }
1157 /*********************************************************************************/
applyRattleSecondPortion(MolecularDynamics * molecularDynamics)1158 static void applyRattleSecondPortion(MolecularDynamics* molecularDynamics)
1159 {
1160 gint i;
1161 gint k;
1162 gint maxIter = 1000;
1163 gdouble omega = 1.0;
1164 gdouble tolerance = 1e-6;
1165 gboolean done = FALSE;
1166 gint nIter = 0;
1167 gint a1 = 0;
1168 gint a2 = 0;
1169 gdouble r2ij;
1170 gdouble dot;
1171 gdouble invMass1;
1172 gdouble invMass2;
1173 gdouble term = 0;
1174 gdouble terms[3];
1175 gdouble d;
1176 Molecule* m = &molecularDynamics->forceField->molecule;
1177 ForceField* forceField = molecularDynamics->forceField;
1178 gdouble deltaMax = 0;
1179
1180 if(forceField->options.rattleConstraints==NOCONSTRAINTS) return;
1181 tolerance /= molecularDynamics->dt;
1182 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1183 {
1184 molecularDynamics->moved[i] = molecularDynamics->forceField->molecule.atoms[i].variable;
1185 molecularDynamics->update[i] = FALSE;
1186 }
1187 /* maxIter *= molecularDynamics->forceField->numberOfRattleConstraintsTerms;*/
1188 do{
1189 nIter++;
1190 done=TRUE;
1191 deltaMax = 0;
1192 for (i = 0; i < molecularDynamics->forceField->numberOfRattleConstraintsTerms; i++)
1193 {
1194 a1 = (gint)molecularDynamics->forceField->rattleConstraintsTerms[0][i];
1195 a2 = (gint)molecularDynamics->forceField->rattleConstraintsTerms[1][i];
1196 r2ij = molecularDynamics->forceField->rattleConstraintsTerms[2][i];
1197 if( !molecularDynamics->moved[a1] && !molecularDynamics->moved[a2] ) continue;
1198 /* here : rattle image for PBC, not yet implemented */
1199 dot = 0;
1200 for (k=0;k<3;k++)
1201 {
1202 d = m->atoms[a2].coordinates[k]-m->atoms[a1].coordinates[k];
1203 dot +=d*(molecularDynamics->velocity[a2][k]-molecularDynamics->velocity[a1][k]);
1204 }
1205 invMass1 = 1/molecularDynamics->forceField->molecule.atoms[a1].prop.masse;
1206 invMass2 = 1/molecularDynamics->forceField->molecule.atoms[a2].prop.masse;
1207 term = -dot / ((invMass1+invMass2)*r2ij);
1208 if(fabs(term)<=tolerance) continue;
1209 /* if(fabs(dot/r2ij)<=tolerance) continue;*/
1210 if(deltaMax<fabs(term)) deltaMax = fabs(term);
1211
1212 done = FALSE;
1213 molecularDynamics->update[a1] = TRUE;
1214 molecularDynamics->update[a2] = TRUE;
1215 term *= omega;
1216
1217 for (k=0;k<3;k++)
1218 {
1219 d = m->atoms[a2].coordinates[k]-m->atoms[a1].coordinates[k];
1220 terms[k] = d*term;
1221 }
1222 for (k=0;k<3;k++) molecularDynamics->velocity[a1][k] -= terms[k]*invMass1;
1223 for (k=0;k<3;k++) molecularDynamics->velocity[a2][k] += terms[k]*invMass2;
1224 }
1225 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1226 {
1227 molecularDynamics->moved[i] = molecularDynamics->update[i];
1228 molecularDynamics->update[i] = FALSE;
1229 }
1230 }while(!done && nIter<maxIter);
1231 if(nIter>=maxIter && deltaMax>tolerance*10)
1232 {
1233 printf(_("Rattle second portion : Warning, velocity constraints not satisfied\n"));
1234 for (i = 0; i < molecularDynamics->forceField->numberOfRattleConstraintsTerms; i++)
1235 {
1236 a1 = (gint)molecularDynamics->forceField->rattleConstraintsTerms[0][i];
1237 a2 = (gint)molecularDynamics->forceField->rattleConstraintsTerms[1][i];
1238 r2ij = 0;
1239 for (k=0;k<3;k++)
1240 {
1241 d = m->atoms[a2].coordinates[k]-m->atoms[a1].coordinates[k];
1242 r2ij +=d*d;
1243 }
1244 dot = 0;
1245 for (k=0;k<3;k++)
1246 {
1247 d = m->atoms[a2].coordinates[k]-m->atoms[a1].coordinates[k];
1248 dot +=d*(molecularDynamics->velocity[a2][k]-molecularDynamics->velocity[a1][k]);
1249 }
1250 invMass1 = 1/molecularDynamics->forceField->molecule.atoms[a1].prop.masse;
1251 invMass2 = 1/molecularDynamics->forceField->molecule.atoms[a2].prop.masse;
1252 term = -dot / ((invMass1+invMass2)*r2ij);
1253 printf("%d %d %s %s r2ij=%f r2Old=%f term=%f\n",
1254 a1,a2,
1255 molecularDynamics->forceField->molecule.atoms[a1].mmType,
1256 molecularDynamics->forceField->molecule.atoms[a2].mmType,
1257 r2ij, molecularDynamics->forceField->rattleConstraintsTerms[2][i],term);
1258 }
1259 StopCalcul=TRUE;
1260 }
1261 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1262 if(!m->atoms[i].variable) for (k=0;k<3;k++) molecularDynamics->velocity[i][k] = 0;
1263 }
1264 /*********************************************************************************/
applyVerlet(MolecularDynamics * molecularDynamics)1265 static void applyVerlet(MolecularDynamics* molecularDynamics)
1266 {
1267 gint i;
1268 gint j;
1269
1270 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS)
1271 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1272 for ( j = 0; j < 3; j++)
1273 molecularDynamics->coordinatesOld[i][j]= molecularDynamics->forceField->molecule.atoms[i].coordinates[j];
1274
1275 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1276 {
1277 if(!molecularDynamics->forceField->molecule.atoms[i].variable) continue;
1278
1279 for ( j = 0; j < 3; j++)
1280 {
1281 molecularDynamics->forceField->molecule.atoms[i].coordinates[j] +=
1282 molecularDynamics->velocity[i][j] * molecularDynamics->dt +
1283 molecularDynamics->a[i][j]*molecularDynamics->dt2_2;
1284 }
1285 for ( j = 0; j < 3; j++)
1286 molecularDynamics->velocity[i][j] += molecularDynamics->a[i][j] * molecularDynamics->dt_2;
1287 }
1288
1289 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS) applyRattleFirstPortion(molecularDynamics);
1290
1291 newAccelaration(molecularDynamics);
1292
1293 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1294 {
1295 if(!molecularDynamics->forceField->molecule.atoms[i].variable) continue;
1296 for ( j = 0; j < 3; j++)
1297 molecularDynamics->velocity[i][j] += molecularDynamics->a[i][j] * molecularDynamics->dt_2;
1298 }
1299 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS) applyRattleSecondPortion(molecularDynamics);
1300 }
1301 /*********************************************************************************/
applyBeeman(MolecularDynamics * molecularDynamics)1302 static void applyBeeman(MolecularDynamics* molecularDynamics)
1303 {
1304 gint i;
1305 gint j;
1306 gdouble terms[3];
1307
1308 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS)
1309 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1310 for ( j = 0; j < 3; j++)
1311 molecularDynamics->coordinatesOld[i][j]= molecularDynamics->forceField->molecule.atoms[i].coordinates[j];
1312
1313
1314 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1315 {
1316 if(!molecularDynamics->forceField->molecule.atoms[i].variable) continue;
1317 for ( j = 0; j < 3; j++)
1318 terms[j] = 5.0*molecularDynamics->a[i][j]-molecularDynamics->aold[i][j];
1319
1320 for ( j = 0; j < 3; j++)
1321 {
1322 molecularDynamics->forceField->molecule.atoms[i].coordinates[j] +=
1323 molecularDynamics->velocity[i][j] * molecularDynamics->dt +
1324 terms[j]*molecularDynamics->dt2_8;
1325 }
1326 for ( j = 0; j < 3; j++)
1327 molecularDynamics->velocity[i][j] += terms[j] * molecularDynamics->dt_8;
1328 }
1329
1330 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS) applyRattleFirstPortion(molecularDynamics);
1331
1332 newAccelaration(molecularDynamics);
1333
1334 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1335 {
1336 if(!molecularDynamics->forceField->molecule.atoms[i].variable) continue;
1337 for ( j = 0; j < 3; j++)
1338 molecularDynamics->velocity[i][j] += (3.0*molecularDynamics->a[i][j]+molecularDynamics->aold[i][j]) * molecularDynamics->dt_8;
1339 }
1340
1341 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS) applyRattleSecondPortion(molecularDynamics);
1342 }
1343 /**********************************************************************/
erfinv(gdouble y)1344 static gdouble erfinv( gdouble y )
1345 {
1346 static gdouble a[] = {0, 0.886226899, -1.645349621, 0.914624893, -0.140543331 };
1347 static gdouble b[] = {0, -2.118377725, 1.442710462, -0.329097515, 0.012229801 };
1348 static gdouble c[] = {0, -1.970840454, -1.624906493, 3.429567803, 1.641345311 };
1349 static gdouble d[] = {0, 3.543889200, 1.637067800 };
1350 gdouble x=1e100, z;
1351
1352 if ( y < -1. ) return x;
1353 if ( y > 1. ) return x;
1354 if ( y >= -.7 )
1355 {
1356 if ( y <= .7 )
1357 {
1358 z = y*y;
1359 x = y * (((a[4]*z+a[3])*z+a[2])*z+a[1]) /
1360 ((((b[4]*z+b[3])*z+b[2])*z+b[1])*z+1);
1361 }
1362 else if ( y < 1 )
1363 {
1364 z = sqrt(-log((1-y)/2));
1365 x = (((c[4]*z+c[3])*z+c[2])*z+c[1]) / ((d[2]*z+d[1])*z+1);
1366 }
1367 }
1368 else
1369 {
1370 z = sqrt(-log((1+y)/2));
1371 x = -(((c[4]*z+c[3])*z+c[2])*z+c[1]) / ((d[2]*z+d[1])*z+1);
1372 }
1373 return x;
1374 }
1375 /**********************************************************************/
getRandVect(gdouble len,gdouble V[])1376 static void getRandVect(gdouble len, gdouble V[])
1377 {
1378 gdouble l = 0;
1379 gint j;
1380 for(j=0;j<3;j++)
1381 {
1382 V [j] = drandom();
1383 l += V[j]*V[j];
1384 }
1385
1386 if(l<=0) return;
1387 l = sqrt(l);
1388 for(j=0;j<3;j++)
1389 V [j] *= len/l;
1390 }
1391 /**********************************************************************/
maxwel(gdouble masse,gdouble temperature)1392 static gdouble maxwel(gdouble masse, gdouble temperature)
1393 {
1394 /*
1395 * physical constants in SI units
1396 * ------------------------------
1397 * Kb = 1.380662 E-23 J/K
1398 * Na = 6.022045 E23 1/mol
1399 * e = 1.6021892 E-19 C
1400 * eps = 8.85418782 E-12 F/m
1401 *
1402 * 1 Kcal = 4184.0 J
1403 * 1 amu = 1.6605655 E-27 Kg
1404 * 1 A = 1.0 E-10 m
1405 *
1406 * Internally, AKMA units are used:
1407 * KBOLTZ = Na *Kb / 1 Kcal
1408 */
1409 /* gdouble Kb = 6.022045e23*1.380662e-23/4184.0;*/
1410 gdouble Kb = 1.9871914e-3;
1411 gdouble beta = sqrt(masse / (2.0*Kb*temperature));
1412 gdouble rho;
1413 gdouble xs, ys, zs;
1414 rho = drandom();
1415 xs = erfinv(rho)/beta;
1416 rho = drandom();
1417 ys = erfinv(rho)/beta;
1418 rho = drandom();
1419 zs = erfinv(rho)/beta;
1420
1421 return sqrt(xs*xs+ys*ys+zs*zs);
1422
1423 }
1424 /*********************************************************************************/
newProperties(MolecularDynamics * molecularDynamics,gchar * comments)1425 static void newProperties(MolecularDynamics* molecularDynamics, gchar* comments)
1426 {
1427 if( molecularDynamics->fileProp == NULL) return;
1428 fprintf(molecularDynamics->fileProp,"time0(fs)\ttime(fs)\tTotal Energy(Kcal/mol)\tPotential Energy(kcal/mol) Kinetic Energy(Kcal/mol)\tT(t) (K)\tTaver(K)\tsigma(T)(K)");
1429 if(comments) fprintf(molecularDynamics->fileProp,"%s\n", comments);
1430 else fprintf(molecularDynamics->fileProp,"\n");
1431 }
1432 /*********************************************************************************/
saveProperties(MolecularDynamics * molecularDynamics,gint iStep0,gint iStep,gchar * comments)1433 static void saveProperties(MolecularDynamics* molecularDynamics, gint iStep0, gint iStep, gchar* comments)
1434 {
1435 /* gdouble fsInAKMA = 1.0/sqrt(1e-10*1e-10*1.6605655e-27*6.022045e23/4184.0)/1e15;*/
1436 static gdouble fsInAKMA = 0.020454828110640;
1437 gdouble dt = molecularDynamics->dt/(fsInAKMA);
1438 static gdouble Ttot = 0;
1439 static gdouble T2tot = 0;
1440 gdouble Taver = 0;
1441 gdouble T2aver = 0;
1442
1443
1444 if( molecularDynamics->fileProp == NULL) return;
1445 if(iStep==1)
1446 {
1447 Ttot = 0;
1448 T2tot = 0;
1449 }
1450 Ttot += molecularDynamics->kelvin;
1451 T2tot += molecularDynamics->kelvin*molecularDynamics->kelvin;
1452 Taver = Ttot/iStep;
1453 T2aver = T2tot/iStep;
1454
1455
1456 fprintf(molecularDynamics->fileProp,"%f\t%f\t%f\t\t%f\t\t%f\t%f\t%f\t%f",
1457 (iStep0)*dt,
1458 (iStep)*dt,
1459 molecularDynamics->totalEnergy,
1460 molecularDynamics->potentialEnergy,
1461 molecularDynamics->kineticEnergy,
1462 molecularDynamics->kelvin,
1463 Taver,
1464 sqrt(fabs(T2aver-Taver*Taver))
1465 );
1466 if(comments) fprintf(molecularDynamics->fileProp,"%s\n", comments);
1467 else fprintf(molecularDynamics->fileProp,"\n");
1468 }
1469 /*********************************************************************************/
saveTrajectory(MolecularDynamics * molecularDynamics,gint iStep)1470 static void saveTrajectory(MolecularDynamics* molecularDynamics, gint iStep)
1471 {
1472 /* gdouble fsInAKMA = 1.0/sqrt(1e-10*1e-10*1.6605655e-27*6.022045e23/4184.0)/1e15;*/
1473 static gdouble fsInAKMA = 0.020454828110640;
1474 gdouble dt = molecularDynamics->dt/(fsInAKMA);
1475 gint i;
1476 if( molecularDynamics->fileTraj == NULL) return;
1477
1478 fprintf(molecularDynamics->fileTraj," %d %f %f %f %f nAtoms time(fs) TotalEnery(Kcal/mol) Kinetic Potential\n",
1479 molecularDynamics->numberOfAtoms,
1480 (iStep)*dt,
1481 molecularDynamics->totalEnergy,
1482 molecularDynamics->kineticEnergy,
1483 molecularDynamics->potentialEnergy
1484 );
1485 fprintf(molecularDynamics->fileTraj," %s\n", "Coord in Ang, Velocity in AKMA, time in fs");
1486
1487 for (i = 0; i < molecularDynamics->numberOfAtoms; i++)
1488 {
1489 fprintf(molecularDynamics->fileTraj," %s %f %f %f %f %f %f %f %s %s %s %d %d\n",
1490 molecularDynamics->forceField->molecule.atoms[i].prop.symbol,
1491 molecularDynamics->forceField->molecule.atoms[i].coordinates[0],
1492 molecularDynamics->forceField->molecule.atoms[i].coordinates[1],
1493 molecularDynamics->forceField->molecule.atoms[i].coordinates[2],
1494 molecularDynamics->velocity[i][0],
1495 molecularDynamics->velocity[i][1],
1496 molecularDynamics->velocity[i][2],
1497 molecularDynamics->forceField->molecule.atoms[i].charge,
1498 molecularDynamics->forceField->molecule.atoms[i].mmType,
1499 molecularDynamics->forceField->molecule.atoms[i].pdbType,
1500 molecularDynamics->forceField->molecule.atoms[i].residueName,
1501 molecularDynamics->forceField->molecule.atoms[i].residueNumber,
1502 molecularDynamics->forceField->molecule.atoms[i].variable
1503 );
1504 }
1505 }
1506
1507 /**********************************************************************/
freeMolecularDynamics(MolecularDynamics * molecularDynamics)1508 void freeMolecularDynamics(MolecularDynamics* molecularDynamics)
1509 {
1510
1511 molecularDynamics->forceField = NULL;
1512 molecularDynamics->numberOfAtoms = 0;
1513 molecularDynamics->updateFrequency = 0;
1514 if(molecularDynamics->velocity)
1515 {
1516 gint i;
1517 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
1518 if(molecularDynamics->velocity[i]) g_free(molecularDynamics->velocity[i]);
1519 g_free(molecularDynamics->velocity);
1520 }
1521 if(molecularDynamics->a)
1522 {
1523 gint i;
1524 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
1525 if(molecularDynamics->a[i]) g_free(molecularDynamics->a[i]);
1526 g_free(molecularDynamics->a);
1527 }
1528 if(molecularDynamics->aold)
1529 {
1530 gint i;
1531 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
1532 if(molecularDynamics->aold[i]) g_free(molecularDynamics->aold[i]);
1533 g_free(molecularDynamics->aold);
1534 }
1535 if(molecularDynamics->coordinatesOld)
1536 {
1537 gint i;
1538 for(i=0;i<molecularDynamics->numberOfAtoms;i++)
1539 if(molecularDynamics->coordinatesOld[i]) g_free(molecularDynamics->coordinatesOld[i]);
1540 g_free(molecularDynamics->coordinatesOld);
1541 }
1542 if(molecularDynamics->moved) g_free(molecularDynamics->moved);
1543 if(molecularDynamics->update) g_free(molecularDynamics->update);
1544 }
1545 /********************************************************************************/
getEKin(MolecularDynamics * molecularDynamics)1546 static gdouble getEKin(MolecularDynamics* molecularDynamics)
1547 {
1548 gdouble ekin = 0;
1549 gint i;
1550 gint j;
1551 gdouble masse;
1552 for ( i = 0; i < molecularDynamics->numberOfAtoms; i++)
1553 {
1554 masse = molecularDynamics->forceField->molecule.atoms[i].prop.masse;
1555 for ( j = 0; j < 3; j++)
1556 ekin += molecularDynamics->velocity[i][j]*molecularDynamics->velocity[i][j]*
1557 masse;
1558 }
1559 return ekin/2;
1560 }
1561 /********************************************************************************/
getKelvin(MolecularDynamics * molecularDynamics)1562 static gdouble getKelvin(MolecularDynamics* molecularDynamics)
1563 {
1564 gint nfree = molecularDynamics->nfree;
1565 static gdouble Kb = 1.9871914e-3;
1566 if(nfree<1) return 0;
1567 return 2*getEKin(molecularDynamics) / ( nfree * Kb);
1568 }
1569 /********************************************************************************/
1570 /*
1571 literature references:
1572
1573 M. P. Allen, "Brownian Dynamics Simulation of a Chemical
1574 Reaction in Solution", Molecular Physics, 40, 1073-1087 (1980)
1575
1576 F. Guarnieri and W. C. Still, "A Rapidly Convergent Simulation
1577 Method: Mixed Monte Carlo / Stochastic Dynamics", Journal of
1578 Computational Chemistry, 15, 1302-1310 (1994)
1579 */
1580 /*********************************************************************************/
getsFrictionalAndRandomForce(MolecularDynamics * molecularDynamics)1581 static void getsFrictionalAndRandomForce(MolecularDynamics* molecularDynamics)
1582 {
1583 gdouble* gamma = molecularDynamics->gamma;
1584 gdouble* positionFriction = molecularDynamics->positionFriction;
1585 gdouble* velocityFriction = molecularDynamics->velocityFriction;
1586 gdouble* accelarationFriction = molecularDynamics->accelarationFriction;
1587 gdouble** positionRandom = molecularDynamics->positionRandom;
1588 gdouble** velocityRandom = molecularDynamics->velocityRandom;
1589 gdouble dt = molecularDynamics->dt;
1590
1591 gint n = molecularDynamics->numberOfAtoms;
1592
1593 gint i;
1594 gint j;
1595 gdouble gdt;
1596 gdouble egdt;
1597 gdouble ktm = 0;
1598 gdouble pterm;
1599 gdouble vterm;
1600 gdouble psig;
1601 gdouble vsig;
1602 gdouble rho;
1603 gdouble rhoc;
1604 gdouble pnorm;
1605 gdouble vnorm;
1606 static gdouble Kb = 1.9871914e-3;
1607
1608 for(i=0;i<n;i++)
1609 gamma[i] = molecularDynamics->friction;
1610
1611 /* printf(" friction = %f\n", molecularDynamics->friction);*/
1612 for(i=0;i<n;i++)
1613 {
1614 gdt = gamma[i] * dt;
1615 /* printf("gdt = %f\n",gdt);*/
1616 if (gdt <= 0.0)
1617 {
1618 positionFriction[i] = 1.0;
1619 velocityFriction[i] = dt;
1620 accelarationFriction[i] = 0.5 * dt * dt;
1621 for(j=0;j<3;j++)
1622 {
1623 positionRandom[i][j] = 0.0;
1624 velocityRandom[i][j] = 0.0;
1625 }
1626 }
1627 else
1628 {
1629 /* analytical expressions when friction coefficient is large */
1630 if (gdt>=0.05)
1631 {
1632 egdt = exp(-gdt);
1633 positionFriction[i] = egdt;
1634 velocityFriction[i] = (1.0-egdt) / gamma[i];
1635 accelarationFriction[i] = (dt-velocityFriction[i]) / gamma[i];
1636 pterm = 2.0*gdt - 3.0 + (4.0-egdt)*egdt;
1637 vterm = 1.0 - egdt*egdt;
1638 rho = (1.0-egdt)*(1.0-egdt) / sqrt(pterm*vterm);
1639 }
1640 /* use series expansions when friction coefficient is small */
1641 else
1642 {
1643 gdouble gdt2 = gdt * gdt;
1644 gdouble gdt3 = gdt * gdt2;
1645 gdouble gdt4 = gdt2 * gdt2;
1646 gdouble gdt5 = gdt2 * gdt3;
1647 gdouble gdt6 = gdt3 * gdt3;
1648 gdouble gdt7 = gdt3 * gdt4;
1649 gdouble gdt8 = gdt4 * gdt4;
1650 gdouble gdt9 = gdt4 * gdt5;
1651 accelarationFriction[i] = (gdt2/2.0 - gdt3/6.0 + gdt4/24.0
1652 - gdt5/120.0 + gdt6/720.0
1653 - gdt7/5040.0 + gdt8/40320.0
1654 - gdt9/362880.0) / gamma[i]/gamma[i];
1655 velocityFriction[i] = dt - gamma[i]*accelarationFriction[i];
1656 positionFriction[i] = 1.0 - gamma[i]*velocityFriction[i];
1657 pterm = 2.0*gdt3/3.0 - gdt4/2.0
1658 + 7.0*gdt5/30.0 - gdt6/12.0
1659 + 31.0*gdt7/1260.0 - gdt8/160.0
1660 + 127.0*gdt9/90720.0;
1661 vterm = 2.0*gdt - 2.0*gdt2 + 4.0*gdt3/3.0
1662 - 2.0*gdt4/3.0 + 4.0*gdt5/15.0
1663 - 4.0*gdt6/45.0 + 8.0*gdt7/315.0
1664 - 2.0*gdt8/315.0 + 4.0*gdt9/2835.0;
1665 rho = sqrt(3.0) * (0.5 - 3.0*gdt/16.0
1666 - 17.0*gdt2/1280.0
1667 + 17.0*gdt3/6144.0
1668 + 40967.0*gdt4/34406400.0
1669 - 57203.0*gdt5/275251200.0
1670 - 1429487.0*gdt6/13212057600.0);
1671 }
1672 ktm = Kb * molecularDynamics->temperature / molecularDynamics->forceField->molecule.atoms[i].prop.masse;
1673 psig = sqrt(ktm*pterm) / gamma[i];
1674 vsig = sqrt(ktm*vterm);
1675 rhoc = sqrt(1.0 - rho*rho);
1676 for(j=0;j<3;j++)
1677 {
1678 pnorm = normal();
1679 vnorm = normal ();
1680 positionRandom[i][j] = psig * pnorm;
1681 velocityRandom[i][j] = vsig * (rho*pnorm+rhoc*vnorm);
1682 }
1683 }
1684 }
1685 }
1686 /*********************************************************************************/
applyStochastic(MolecularDynamics * molecularDynamics)1687 static void applyStochastic(MolecularDynamics* molecularDynamics)
1688 {
1689 gdouble* positionFriction = molecularDynamics->positionFriction;
1690 gdouble* velocityFriction = molecularDynamics->velocityFriction;
1691 gdouble* accelarationFriction = molecularDynamics->accelarationFriction;
1692 gdouble** positionRandom = molecularDynamics->positionRandom;
1693 gdouble** velocityRandom = molecularDynamics->velocityRandom;
1694 gdouble**v = molecularDynamics->velocity;
1695 gdouble**a = molecularDynamics->a;
1696
1697 gint n = molecularDynamics->numberOfAtoms;
1698 gint i;
1699 gint j;
1700 AtomMol* atoms = molecularDynamics->forceField->molecule.atoms;
1701
1702 getsFrictionalAndRandomForce(molecularDynamics);
1703
1704 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS)
1705 for (i = 0; i < n; i++)
1706 for ( j = 0; j < 3; j++)
1707 molecularDynamics->coordinatesOld[i][j]= molecularDynamics->forceField->molecule.atoms[i].coordinates[j];
1708
1709
1710 for(i=0;i<n;i++)
1711 {
1712 if(!molecularDynamics->forceField->molecule.atoms[i].variable) continue;
1713 for(j=0;j<3;j++)
1714 atoms[i].coordinates[j] += v[i][j]*velocityFriction[i] + a[i][j]*accelarationFriction[i] + positionRandom[i][j];
1715 for(j=0;j<3;j++)
1716 v[i][j] = v[i][j]*positionFriction[i] + 0.5*a[i][j]*velocityFriction[i];
1717 }
1718 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS) applyRattleFirstPortion(molecularDynamics);
1719 newAccelaration(molecularDynamics);
1720
1721 for (i = 0; i < n; i++)
1722 {
1723 if(!molecularDynamics->forceField->molecule.atoms[i].variable) continue;
1724 for ( j = 0; j < 3; j++)
1725 v[i][j] += 0.5*a[i][j]*velocityFriction[i] + velocityRandom[i][j];
1726 }
1727 if(molecularDynamics->forceField->options.rattleConstraints!=NOCONSTRAINTS) applyRattleSecondPortion(molecularDynamics);
1728 computeEnergies(molecularDynamics);
1729 }
1730 /*********************************************************************************/
drandom()1731 static gdouble drandom()
1732 {
1733 return (rand()/(gdouble)RAND_MAX);
1734 }
1735 /*********************************************************************************/
1736 /* "normal" generates a random number from a normal Gaussian
1737 distribution with a mean of zero and a variance of one
1738 */
normal()1739 static gdouble normal()
1740 {
1741 gdouble v1,v2,rsq;
1742 gdouble factor;
1743 static gdouble store;
1744 static gboolean compute = TRUE;
1745
1746 if (compute)
1747 {
1748 do{
1749 v1 = 2.0 * drandom() - 1.0;
1750 v2 = 2.0 * drandom () - 1.0;
1751 rsq = v1*v1 + v2*v2;
1752 }while(rsq >= 1.0);
1753 compute = FALSE;
1754 factor = sqrt(-2.0*log(rsq)/rsq);
1755 store = v1 * factor;
1756 return v2 * factor;
1757 }
1758 /* use the second random value computed at the last call */
1759 else
1760 {
1761 compute = TRUE;
1762 return store;
1763 }
1764 }
1765 /*********************************************************************************/
1766