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