source: pcp/src/ions.c@ 3716b2

Last change on this file since 3716b2 was 3716b2, checked in by Frederik Heber <heber@…>, 17 years ago

commented-out lines
  • Property mode set to 100644
File size: 45.1 KB
Line 
1/** \file ions.c
2 * Ionic Force, speed, coordinate and energy calculation.
3 * Herein are all the routines for the molecular dynamics calculations such as
4 * reading of configuration files IonsInitRead() and
5 * free'ing RemoveIonsRead(),
6 * summing up forces CalculateIonForce(),
7 * correcting for temperature CorrectVelocities(),
8 * or for fixation center of balance CorrectForces(),
9 * outputting them to file OutputIonForce(),
10 * calculating kinetic CalculateEnergyIonsU(), ewald CalculateEwald() energies,
11 * moving ion UpdateIonsR() (position) UpdateIonsU() (speed),
12 * scaling to match some temperature ScaleTemp()
13 * are gathered.
14 * Finally, needed for structure optimization GetOuterStop() evaluates the stop
15 * condition of a minimal mean force acting on each ion.
16 *
17 Project: ParallelCarParrinello
18 \author Jan Hamaekers
19 \date 2000
20
21 File: ions.c
22 $Id: ions.c,v 1.34 2007/02/05 16:15:27 foo Exp $
23*/
24
25#include <stdlib.h>
26#include <stdio.h>
27#include <stddef.h>
28#include <math.h>
29#include <string.h>
30#include <unistd.h>
31#include <gsl/gsl_randist.h>
32#include "mymath.h"
33
34
35// use double precision fft when we have it
36#ifdef HAVE_CONFIG_H
37#include <config.h>
38#endif
39
40#ifdef HAVE_DFFTW_H
41#include "dfftw.h"
42#else
43#include "fftw.h"
44#endif
45
46#include "data.h"
47#include "errors.h"
48#include "helpers.h"
49#include "mymath.h"
50#include "ions.h"
51#include "init.h"
52
53/** Readin of Thermostat related values from parameter file.
54 * \param *P Problem at hand
55 * \param *source parameter file
56 */
57void InitThermostats(struct Problem *P, FILE *source)
58{
59 struct FileData *F = &P->Files;
60 struct Ions *I = &P->Ion;
61 char *thermo = MallocString(12, "IonsInitRead: thermo");
62 int j;
63
64 // initialize the naming stuff and all
65 I->ThermostatImplemented = (int *) Malloc((MaxThermostats)*(sizeof(int)), "IonsInitRead: *ThermostatImplemented");
66 I->ThermostatNames = (char **) Malloc((MaxThermostats)*(sizeof(char *)), "IonsInitRead: *ThermostatNames");
67 for (j=0;j<MaxThermostats;j++)
68 I->ThermostatNames[j] = (char *) MallocString(12*(sizeof(char)), "IonsInitRead: ThermostatNames[]");
69 strcpy(I->ThermostatNames[0],"None");
70 I->ThermostatImplemented[0] = 1;
71 strcpy(I->ThermostatNames[1],"Woodcock");
72 I->ThermostatImplemented[1] = 1;
73 strcpy(I->ThermostatNames[2],"Gaussian");
74 I->ThermostatImplemented[2] = 1;
75 strcpy(I->ThermostatNames[3],"Langevin");
76 I->ThermostatImplemented[3] = 1;
77 strcpy(I->ThermostatNames[4],"Berendsen");
78 I->ThermostatImplemented[4] = 1;
79 strcpy(I->ThermostatNames[5],"NoseHoover");
80 I->ThermostatImplemented[5] = 1;
81 I->Thermostat = -1;
82 if (F->MeOutMes) fprintf(F->TemperatureFile, "# %s\t%s\t%s --- ","Step", "ActualTemp(1)", "ActualTemp(2)");
83
84 // read desired Thermostat from file along with needed additional parameters
85 if (ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 1, 1, string_type, thermo, 1, optional)) {
86 if (strcmp(thermo, I->ThermostatNames[0]) == 0) { // None
87 if (I->ThermostatImplemented[0] == 1) {
88 I->Thermostat = None;
89 } else {
90 if (P->Par.me == 0)
91 fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[0]);
92 I->Thermostat = None;
93 }
94 } else if (strcmp(thermo, I->ThermostatNames[1]) == 0) { // Woodcock
95 if (I->ThermostatImplemented[1] == 1) {
96 I->Thermostat = Woodcock;
97 ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, int_type, &P->R.ScaleTempStep, 1, critical); // read scaling frequency
98 } else {
99 if (P->Par.me == 0)
100 fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[1]);
101 I->Thermostat = None;
102 }
103 } else if (strcmp(thermo, I->ThermostatNames[2]) == 0) { // Gaussian
104 if (I->ThermostatImplemented[2] == 1) {
105 I->Thermostat = Gaussian;
106 ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, int_type, &P->R.ScaleTempStep, 1, critical); // read collision rate
107 } else {
108 if (P->Par.me == 0)
109 fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[2]);
110 I->Thermostat = None;
111 }
112 } else if (strcmp(thermo, I->ThermostatNames[3]) == 0) { // Langevin
113 if (I->ThermostatImplemented[3] == 1) {
114 I->Thermostat = Langevin;
115 ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.TempFrequency, 1, critical); // read gamma
116 if (ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 3, 1, double_type, &P->R.alpha, 1, optional)) {
117 if (P->Par.me == 0)
118 fprintf(stderr,"(%i) Extended Stochastic Thermostat detected with interpolation coefficient %lg\n", P->Par.me, P->R.alpha);
119 } else {
120 P->R.alpha = 1.;
121 }
122 } else {
123 if (P->Par.me == 0)
124 fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[3]);
125 I->Thermostat = None;
126 }
127 } else if (strcmp(thermo, I->ThermostatNames[4]) == 0) { // Berendsen
128 if (I->ThermostatImplemented[4] == 1) {
129 I->Thermostat = Berendsen;
130 ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.TempFrequency, 1, critical); // read \tau_T
131 } else {
132 if (P->Par.me == 0)
133 fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[4]);
134 I->Thermostat = None;
135 }
136 } else if (strcmp(thermo, I->ThermostatNames[5]) == 0) { // Nose-Hoover
137 if (I->ThermostatImplemented[5] == 1) {
138 I->Thermostat = NoseHoover;
139 ParseForParameter(P->Call.out[ReadOut],source,"Thermostat", 0, 2, 1, double_type, &P->R.HooverMass, 1, critical); // read Hoovermass
140 P->R.alpha = 0.;
141 } else {
142 if (P->Par.me == 0)
143 fprintf(stderr, "(%i) Warning: %s thermostat not implemented, falling back to None.", P->Par.me, I->ThermostatNames[5]);
144 I->Thermostat = None;
145 }
146 }
147 if (I->Thermostat == -1) {
148 if (P->Par.me == 0)
149 fprintf(stderr,"(%i) Warning: thermostat name was not understood!\n", P->Par.me);
150 I->Thermostat = None;
151 }
152 } else {
153 if ((P->R.MaxOuterStep > 0) && (P->R.TargetTemp != 0))
154 debug(P, "No thermostat chosen despite finite temperature MD, falling back to None.");
155 I->Thermostat = None;
156 }
157 if (F->MeOutMes) fprintf(F->TemperatureFile, "%s Thermostat\n",I->ThermostatNames[(int)I->Thermostat]);
158 Free(thermo, "InitThermostats: thermo");
159}
160
161/** Parses for a given keyword the ion coordinates and velocites.
162 * \param *P Problem at hand, contains pointer to Ion and IonType structure
163 * \param *source
164 * \param *keyword keyword string
165 * \param repetition which repeated version of the keyword shall be read by ReadParameters()
166 * \param relative whether atom coordinates are relative to unit cell size or absolute
167 * \param Bohr conversion factor to atomiclength (e.g. 1.0 if coordinates in bohr radius, 0.529... if given in Angstrom)
168 * \param sigma variance of maxwell-boltzmann distribution
169 * \param *R where to put the read coordinates
170 * \param *U where to put the read velocities
171 * \param *IMT whether its MoveType#FixedIon or not
172 * \return if 1 - success, 0 - failure
173 */
174int ParseIonsCoordinatesAndVelocities(struct Problem *P, FILE *source, char *keyword, int repetition, int relative, double sigma, double *R, double *U, int *IMT)
175{
176 //struct RunStruct *Run = &P->R;
177 struct Ions *I = &P->Ion;
178 struct Lattice *L = &P->Lat;
179 int k;
180 double R_trafo[3];
181 gsl_rng * r;
182 const gsl_rng_type * T;
183 double Bohr = (I->IsAngstroem) ? ANGSTROEMINBOHRRADIUS : 1.;
184
185 gsl_rng_env_setup();
186 T = gsl_rng_default;
187 r = gsl_rng_alloc (T);
188 if (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 1, 1, double_type, &R_trafo[0], repetition, optional) &&
189 ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 2, 1, double_type, &R_trafo[1], repetition, optional) &&
190 ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 3, 1, double_type, &R_trafo[2], repetition, optional) &&
191 ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 4, 1, int_type, IMT, repetition, optional)) {
192 // if read successful, then try also parsing velocities
193 if ((ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 5, 1, double_type, &U[0], repetition, optional)) &&
194 (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 6, 1, double_type, &U[1], repetition, optional)) &&
195 (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 7, 1, double_type, &U[2], repetition, optional))) {
196 //nothing
197 } else { // if no velocities specified, set to maxwell-boltzmann distributed
198 if ((I->Thermostat != None)) {
199 if (P->Par.me == 0) fprintf(stderr, "(%i) Missing velocities of ion %s: Maxwell-Boltzmann with %lg\n", P->Par.me, keyword, sigma);
200 for (k=0;k<NDIM;k++)
201 U[k] = gsl_ran_gaussian (r, sigma);
202 } else {
203 for (k=0;k<NDIM;k++)
204 U[k] = 0.;
205 }
206 }
207 // change if coordinates were relative
208 if (relative) {
209 //fprintf(stderr,"(%i)Ion coordinates are relative %i ... \n",P->Par.me, relative);
210 RMat33Vec3(R, L->RealBasis, R_trafo); // multiply with real basis
211 } else {
212 for (k=0;k<NDIM;k++) // otherweise just copy to vector
213 R[k] = R_trafo[k];
214 }
215 // check if given MoveType is valid
216 if (*IMT < 0 || *IMT >= MaxIonMoveType) {
217 fprintf(stderr,"Bad Ion MoveType set to MoveIon for %s\n", keyword);
218 *IMT = 0;
219 }
220 if (!relative) {
221 //fprintf(stderr,"converting with %lg\n", Bohr);
222 SM(R, Bohr, NDIM); // convert angstroem to bohrradius
223 SM(U, Bohr, NDIM); // convert angstroem to bohrradius
224 }
225 // finally periodicly correct coordinates to remain in unit cell
226 RMat33Vec3(R_trafo, L->InvBasis, R);
227 for (k=0; k <NDIM; k++) { // periodic correction until within unit cell
228 while (R_trafo[k] < 0)
229 R_trafo[k] += 1.0;
230 while (R_trafo[k] >= 1.0)
231 R_trafo[k] -= 1.0;
232 }
233 RMat33Vec3(R, L->RealBasis, R_trafo);
234 // print coordinates if desired
235 if ((P->Call.out[ReadOut])) {
236 fprintf(stderr,"(%i) coordinates of Ion %s: (x,y,z) = (",P->Par.me, keyword);
237 for(k=0;k<NDIM;k++) {
238 fprintf(stderr,"%lg ",R[k]);
239 if (k != NDIM-1) fprintf(stderr,", ");
240 else fprintf(stderr,")\n");
241 }
242 }
243 gsl_rng_free (r);
244 return 1;
245 } else {
246 for (k=0;k<NDIM;k++) {
247 U[k] = R[k] = 0.;
248 }
249 *IMT = 0;
250 gsl_rng_free (r);
251 return 0;
252 }
253}
254
255/** Readin of the ion section of parameter file.
256 * Among others the following paramaters are read from file:
257 * Ions::MaxTypes, IonType::Max_IonsOfType, IonType::Z, IonType::R,
258 * IonType:IMT, ...
259 * \param P Problem at hand (containing Ions and Lattice structures)
260 * \param *source file pointer
261 */
262void IonsInitRead(struct Problem *P, FILE *source) {
263 struct Ions *I = &P->Ion;
264 //struct RunStruct *Run = &P->R;
265 int i,it,j,IMT,d,me, count;
266 int relative; // whether read-in coordinate are relative (1) or absolute
267 double Bohr = ANGSTROEMINBOHRRADIUS; /* Angstroem */
268 double sigma, R[3], U[3];
269 struct IonConstrained *ptr = NULL;
270
271 I->Max_TotalIons = 0;
272 I->TotalZval = 0;
273 MPI_Comm_rank(MPI_COMM_WORLD, &me);
274 ParseForParameter(P->Call.out[ReadOut],source,"RCut", 0, 1, 1, double_type, &I->R_cut, 1, critical);
275 ParseForParameter(P->Call.out[ReadOut],source,"IsAngstroem", 0, 1, 1, int_type, &I->IsAngstroem, 1, critical);
276 if (!I->IsAngstroem) {
277 Bohr = 1.0;
278 } else { // adapt RealBasis (is in Angstroem as well)
279 SM(P->Lat.RealBasis, Bohr, NDIM_NDIM);
280 RMatReci3(P->Lat.InvBasis, P->Lat.RealBasis);
281 }
282 ParseForParameter(P->Call.out[ReadOut],source,"MaxTypes", 0, 1, 1, int_type, &I->Max_Types, 1, critical);
283 I->I = (struct IonType *) Malloc(sizeof(struct IonType)*I->Max_Types, "IonsInitRead: IonType");
284 if (!ParseForParameter(P->Call.out[ReadOut],source,"RelativeCoord", 0, 1, 1, int_type, &relative, 1, optional))
285 relative = 0;
286 if (!ParseForParameter(P->Call.out[ReadOut],source,"StructOpt", 0, 1, 1, int_type, &I->StructOpt, 1, optional))
287 I->StructOpt = 0;
288
289 InitThermostats(P, source);
290
291 /* Ions Data */
292 I->RLatticeVec = NULL;
293 I->TotalMass = 0;
294 char *free_name, *name;
295 name = free_name = Malloc(255*sizeof(char),"IonsInitRead: Name");
296 for (i=0; i < I->Max_Types; i++) {
297 sprintf(name,"Ion_Type%i",i+1);
298 I->I[i].corecorr = NotCoreCorrected;
299 I->I[i].Name = MallocString(255, "IonsInitRead: Name");
300 I->I[i].Symbol = MallocString(255, "IonsInitRead: Symbol");
301 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 1, 1, int_type, &I->I[i].Max_IonsOfType, 1, critical);
302 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, int_type, &I->I[i].Z, 1, critical);
303 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &I->I[i].rgauss, 1, critical);
304 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &I->I[i].l_max, 1, critical);
305 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 5, 1, int_type, &I->I[i].l_loc, 1, critical);
306 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 6, 1, double_type, &I->I[i].IonMass, 1, critical);
307 if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 7, 1, string_type, &I->I[i].Name[0], 1, optional))
308 sprintf(&I->I[i].Name[0],"type%i",i);
309 if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 8, 1, string_type, &I->I[i].Symbol[0], 1, optional))
310 sprintf(&I->I[i].Symbol[0],"t%i",i);
311 if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) Element name: %s, symbol: %s\n",P->Par.me, I->I[i].Name, I->I[i].Symbol);
312 I->I[i].moment = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: moment");
313 I->I[i].sigma = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma");
314 I->I[i].sigma_rezi = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_rezi");
315 I->I[i].sigma_PAS = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_PAS");
316 I->I[i].sigma_rezi_PAS = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_rezi_PAS");
317 for (it=0;it<I->I[i].Max_IonsOfType;it++) {
318 I->I[i].moment[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: moment");
319 I->I[i].sigma[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: sigma");
320 I->I[i].sigma_rezi[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: sigma_rezi");
321 I->I[i].sigma_PAS[it] = (double *) Malloc(sizeof(double) * NDIM, "IonsInitRead: sigma_PAS");
322 I->I[i].sigma_rezi_PAS[it] = (double *) Malloc(sizeof(double) * NDIM, "IonsInitRead: sigma_rezi_PAS");
323 }
324 //I->I[i].IonFac = I->I[i].IonMass;
325 I->I[i].IonMass *= Units2Electronmass; // in config given in amu not in "atomic units" (electron mass, m_e = 1)
326 // proportional factor between electron and proton mass
327 I->I[i].IonFac = Units2Electronmass/I->I[i].IonMass;
328 I->TotalMass += I->I[i].Max_IonsOfType*I->I[i].IonMass;
329 sigma = sqrt(P->R.TargetTemp/I->I[i].IonMass);
330 I->I[i].ZFactor = -I->I[i].Z*4.*PI;
331 I->I[i].alpha = (double *) Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: alpha");
332 I->I[i].R = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: R");
333 I->I[i].R_old = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: R_old");
334 I->I[i].R_old_old = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: R_old_old");
335 I->I[i].Constraints = (struct IonConstrained **) Malloc(sizeof(struct IonConstrained *)*I->I[i].Max_IonsOfType, "IonsInitRead: **Constraints");
336 for (it=0;it<I->I[i].Max_IonsOfType;it++)
337 I->I[i].Constraints[it] = NULL;
338 I->I[i].FIon = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIon");
339 I->I[i].FIon_old = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIon_old");
340 SetArrayToDouble0(I->I[i].FIon_old, NDIM*I->I[i].Max_IonsOfType);
341 I->I[i].SearchDir = Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: SearchDir");
342 I->I[i].GammaA = Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: GammaA");
343 SetArrayToDouble0(I->I[i].GammaA, I->I[i].Max_IonsOfType);
344 I->I[i].U = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: U");
345 SetArrayToDouble0(I->I[i].U, NDIM*I->I[i].Max_IonsOfType);
346 I->I[i].SFactor = NULL;
347 I->I[i].IMT = (enum IonMoveType *) Malloc(sizeof(enum IonMoveType) * I->I[i].Max_IonsOfType, "IonsInitRead: IMT");
348 I->I[i].FIonL = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIonL");
349 I->I[i].FIonNL = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIonNL");
350 I->I[i].FMagnetic = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FMagnetic");
351 I->I[i].FConstraint = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FConstraint");
352 I->I[i].FEwald = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FEwald");
353 I->I[i].alpha = (double *) Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: alpha");
354 // now parse ion coordination
355 for (j=0; j < I->I[i].Max_IonsOfType; j++) {
356 I->I[i].alpha[j] = 2.;
357 sprintf(name,"Ion_Type%i_%i",i+1,j+1);
358 for (d=0; d <NDIM; d++) // transfor to old and old_old coordinates as well
359 I->I[i].FMagnetic[d+j*NDIM] = 0.;
360 // first ones are starting positions
361 if ((count = ParseIonsCoordinatesAndVelocities(P, source, name, 1, relative, sigma, &I->I[i].R[NDIM*j], &I->I[i].U[NDIM*j], &IMT))) {
362 for (d=0; d <NDIM; d++) // transfor to old and old_old coordinates as well
363 I->I[i].R_old_old[d+NDIM*j] = I->I[i].R_old[d+NDIM*j] = I->I[i].R[d+NDIM*j];
364 I->I[i].IMT[j] = IMT;
365 if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) R[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,I->I[i].R[0+NDIM*j],I->I[i].R[1+NDIM*j],I->I[i].R[2+NDIM*j]);
366 if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) U[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,I->I[i].U[0+NDIM*j],I->I[i].U[1+NDIM*j],I->I[i].U[2+NDIM*j]);
367 if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) IMT[%d,%d] = %i\n", P->Par.me, i,j,I->I[i].IMT[j]);
368 } else {
369 Error(SomeError, "Could not find or parse coordinates of an ion");
370 }
371 // following ones are constrained motions
372 while (ParseIonsCoordinatesAndVelocities(P, source, name, ++count, relative, sigma, R, U, &IMT)) {
373 if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) Constrained read %d ...\n", P->Par.me, count-1);
374 ptr = AddConstraintItem(&I->I[i], j); // allocate memory
375 for (d=0; d< NDIM; d++) { // fill the constraint item
376 ptr->R[d] = R[d];
377 ptr->U[d] = U[d];
378 }
379 ptr->IMT = IMT;
380 ptr->step = count-1;
381 if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) R[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,ptr->R[0],ptr->R[1],ptr->R[2]);
382 if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) U[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,ptr->U[0],ptr->U[1],ptr->U[2]);
383 if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) IMT[%d,%d] = %i\n", P->Par.me, i,j,ptr->IMT);
384 }
385 //if (P->Call.out[ReadOut])
386 fprintf(stderr,"(%i) A total of %d additional constrained moves for ion (%d,%d) was parsed.\n", P->Par.me, count-2, i,j);
387
388 I->Max_TotalIons++;
389 }
390 }
391 Free(free_name, "IonsInitRead: free_name");
392 I->Max_Max_IonsOfType = 0;
393 for (i=0; i < I->Max_Types; i++)
394 I->Max_Max_IonsOfType = MAX(I->Max_Max_IonsOfType, I->I[i].Max_IonsOfType);
395 I->FTemp = (double *) Malloc(sizeof(double) * NDIM * I->Max_Max_IonsOfType, "IonsInitRead:");
396}
397
398/** Allocates memory for a IonConstrained item and adds to end of list.
399 * \param *I IonType structure
400 * \param ion which ion of the type
401 * \return pointer to created item
402 */
403struct IonConstrained * AddConstraintItem(struct IonType *I, int ion)
404{
405 struct IonConstrained **ptr = &(I->Constraints[ion]);
406 while (*ptr != NULL) { // advance to end of list
407 ptr = &((*ptr)->next);
408 }
409 // allocated memory for structure and items within
410 *ptr = (struct IonConstrained *) Malloc(sizeof(struct IonConstrained), "CreateConstraintItem: IonConstrained");
411 (*ptr)->R = (double *) Malloc(sizeof(double)*NDIM, "AddConstraintItem: *R");
412 (*ptr)->U = (double *) Malloc(sizeof(double)*NDIM, "AddConstraintItem: *U");
413 (*ptr)->next = NULL;
414 return *ptr;
415}
416
417/** Removes first of item of IonConstrained list for given \a ion.
418 * Frees memory of items within and then
419 * \param *I IonType structure
420 * \param ion which ion of the type
421 * \return 1 - success, 0 - failure (list is already empty, start pointer was NULL)
422 */
423int RemoveConstraintItem(struct IonType *I, int ion)
424{
425 struct IonConstrained **ptr = &(I->Constraints[ion]);
426 struct IonConstrained *next_ptr;
427 if (*ptr != NULL) {
428 next_ptr = (*ptr)->next;
429 Free((*ptr)->R, "RemoveConstraintItem: R");
430 Free((*ptr)->U, "RemoveConstraintItem: U");
431 Free((*ptr), "RemoveConstraintItem: IonConstrained structure");
432 *ptr = next_ptr; // set start of list to next item (automatically is null if ptr was last item)
433 return 1;
434 } else {
435 return 0;
436 }
437}
438
439/** Calculated Ewald force.
440 * \f[
441 * E_{K-K} = \frac{1}{2} \sum_L \sum_{K=\{(I_s,I_a,J_s,J_a)|R_{I_s,I_a} - R_{J_s,J_a} - L\neq 0\}} \frac{Z_{I_s} Z_{J_s}} {|R_{I_s,I_a} - R_{J_s,J_a} - L|}
442 * \textnormal{erfc} \Bigl ( \frac{|R_{I_s,I_a} - R_{J_s,J_a} - L|} {\sqrt{r_{I_s}^{Gauss}+r_{J_s}^{Gauss}}}\Bigr )
443 * \qquad (4.10)
444 * \f]
445 * Calculates the core-to-core force between the ions of all super cells.
446 * In order for this series to converge there must be a certain summation
447 * applied, which is the ewald summation and which is nothing else but to move
448 * in a circular motion from the current cell to the outside up to Ions::R_cut.
449 * \param *P Problem at hand
450 * \param first additional calculation beforehand if != 0
451 */
452void CalculateEwald(struct Problem *P, int first)
453{
454 int r,is,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec;
455 int MaxPar=P->Par.procs, ParMe=P->Par.me, MaxLocalVec, StartVec;
456 double rcut2,r2,erre2,addesr,addpre,arg,fac,gkl,esr,fxx,repand;
457 double R1[3];
458 double R2[3];
459 double R3[3];
460 double t[3];
461 struct Lattice *L = &P->Lat;
462 struct PseudoPot *PP = &P->PP;
463 struct Ions *I = &P->Ion;
464 if (I->R_cut != 0.0) {
465 if (first) {
466 SpeedMeasure(P, InitSimTime, StopTimeDo);
467 rcut2 = I->R_cut*I->R_cut;
468 ActualVec =0;
469 MaxCell = (int)5.*I->R_cut/pow(L->Volume,1./3.);
470 if (MaxCell < 2) MaxCell = 2;
471 for (i = -MaxCell; i <= MaxCell; i++) {
472 for (j = -MaxCell; j <= MaxCell; j++) {
473 for (k = -MaxCell; k <= MaxCell; k++) {
474 r2 = 0.0;
475 for (ir=0; ir <3; ir++) {
476 t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
477 r2 = t[ir]*t[ir];
478 }
479 Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
480 if ((r2 <= rcut2) || Is_Neighb_Cell) {
481 ActualVec++;
482 }
483 }
484 }
485 }
486 MaxVec = ActualVec;
487 I->MaxVec = ActualVec;
488 MaxLocalVec = MaxVec / MaxPar;
489 StartVec = ParMe*MaxLocalVec;
490 r = MaxVec % MaxPar;
491 if (ParMe >= r) {
492 StartVec += r;
493 } else {
494 StartVec += ParMe;
495 MaxLocalVec++;
496 }
497 I->MaxLocalVec = MaxLocalVec;
498 LocalActualVec = ActualVec = 0;
499 I->RLatticeVec = (double *)
500 Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald:");
501 for (i = -MaxCell; i <= MaxCell; i++) {
502 for (j = -MaxCell; j <= MaxCell; j++) {
503 for (k = -MaxCell; k <= MaxCell; k++) {
504 r2 = 0.0;
505 for (ir=0; ir <3; ir++) {
506 t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
507 r2 = t[ir]*t[ir];
508 }
509 Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
510 if ((r2 <= rcut2) || Is_Neighb_Cell) {
511 if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) {
512 I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
513 I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
514 I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
515 LocalActualVec++;
516 }
517 ActualVec++;
518 }
519 }
520 }
521 }
522 SpeedMeasure(P, InitSimTime, StartTimeDo);
523 }
524 SpeedMeasure(P, EwaldTime, StartTimeDo);
525 esr = 0.0;
526 for (is1=0;is1 < I->Max_Types; is1++)
527 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
528 for (i=0; i< NDIM; i++)
529 I->FTemp[i+NDIM*(ia1)] = 0;
530 for (is1=0;is1 < I->Max_Types; is1++) {
531 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
532 for (i=0;i<3;i++) {
533 R1[i]=I->I[is1].R[i+NDIM*ia1];
534 }
535 for (ir2=0;ir2<I->MaxLocalVec;ir2++) {
536 for (is2=0;is2 < I->Max_Types; is2++) {
537 gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
538 for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
539 for (i=0;i<3;i++) {
540 R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2];
541 }
542 erre2=0.0;
543 for (i=0;i<3;i++) {
544 R3[i] = R1[i]-R2[i];
545 erre2+=R3[i]*R3[i];
546 }
547 if (erre2 > MYEPSILON) {
548 arg=sqrt(erre2);
549 fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;
550
551 arg *= gkl;
552 addesr = derf(arg);
553 addesr = (1.0-addesr)*fac;
554 esr += addesr;
555 addpre=exp(-arg*arg)*gkl;
556 addpre=PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;
557 repand=(addesr+addpre)/erre2;
558 for (i=0;i<3;i++) {
559 fxx=repand*R3[i];
560 /*fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);*/
561 I->FTemp[i+NDIM*(ia1)] += fxx;
562 I->FTemp[i+NDIM*(ia2)] -= fxx;
563 }
564 }
565 }
566 }
567 }
568 }
569 }
570 } else {
571 esr = 0.0;
572 for (is1=0;is1 < I->Max_Types; is1++)
573 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
574 for (i=0; i< NDIM; i++)
575 I->FTemp[i+NDIM*(ia1)] = 0.;
576 }
577 SpeedMeasure(P, EwaldTime, StopTimeDo);
578 for (is=0;is < I->Max_Types; is++) {
579 MPI_Allreduce (I->FTemp , I->I[is].FEwald, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);
580 }
581 MPI_Allreduce ( &esr, &L->E->AllTotalIonsEnergy[EwaldEnergy], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm);
582}
583
584/** Sum up all forces acting on Ions.
585 * IonType::FIon = IonType::FEwald + IonType::FIonL + IonType::FIonNL for all
586 * dimensions #NDIM and each Ion of IonType
587 * \param *P Problem at hand
588 */
589void CalculateIonForce(struct Problem *P)
590{
591 struct Ions *I = &P->Ion;
592 int i,j,d;
593 for (i=0; i < I->Max_Types; i++)
594 for (j=0; j < I->I[i].Max_IonsOfType; j++)
595 for (d=0; d<NDIM;d++)
596 I->I[i].FIon[d+j*NDIM] = I->I[i].FEwald[d+j*NDIM] + I->I[i].FIonL[d+j*NDIM] + I->I[i].FIonNL[d+j*NDIM];
597}
598
599/** Write Forces to FileData::ForcesFile.
600 * goes through all Ions per IonType and each dimension of #NDIM and prints in one line:
601 * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL,
602 * non-local IonType::FIonNL, ewald force IonType::FEwald
603 * \param *P Problem at hand
604 */
605void OutputIonForce(struct Problem *P)
606{
607 struct RunStruct *R = &P->R;
608 struct FileData *F = &P->Files;
609 struct Ions *I = &P->Ion;
610 FILE *fout = NULL;
611 int i,j;
612 if (F->MeOutMes != 1) return;
613 fout = F->ForcesFile;
614 for (i=0; i < I->Max_Types; i++)
615 for (j=0; j < I->I[i].Max_IonsOfType; j++)
616 fprintf(fout, "%i\t%i\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n", i,j,
617 I->I[i].R[0+j*NDIM], I->I[i].R[1+j*NDIM], I->I[i].R[2+j*NDIM],
618 I->I[i].FIon[0+j*NDIM], I->I[i].FIon[1+j*NDIM], I->I[i].FIon[2+j*NDIM],
619 I->I[i].FIonL[0+j*NDIM], I->I[i].FIonL[1+j*NDIM], I->I[i].FIonL[2+j*NDIM],
620 I->I[i].FIonNL[0+j*NDIM], I->I[i].FIonNL[1+j*NDIM], I->I[i].FIonNL[2+j*NDIM],
621 I->I[i].FEwald[0+j*NDIM], I->I[i].FEwald[1+j*NDIM], I->I[i].FEwald[2+j*NDIM]);
622 fflush(fout);
623 fprintf(fout, "MeanForce:\t%e\n", R->MeanForce[0]);
624}
625
626/** Reads Forces from CallOptions::ForcesFile.
627 * goes through all Ions per IonType and each dimension of #NDIM and parses in one line:
628 * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL,
629 * non-local IonType::FIonNL, ewald force IonType::FEwald
630 * \param *P Problem at hand
631 */
632void ParseIonForce(struct Problem *P)
633{
634 //struct RunStruct *Run = &P->R;
635 struct Ions *I = &P->Ion;
636 FILE *finput = NULL;
637 char line[1024];
638 int i,j;
639 double i1, j1;
640 double R[NDIM];
641
642 if (P->Call.ForcesFile == NULL) return;
643 finput = fopen(P->Call.ForcesFile, "r");
644 //fprintf(stderr, "File pointer says %p\n", finput);
645 if (finput == NULL)
646 Error(InitReading, "ForcesFile does not exist.");
647 fscanf(finput, "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n", line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line,line); // skip header line
648 for (i=0; i < I->Max_Types; i++)
649 for (j=0; j < I->I[i].Max_IonsOfType; j++)
650 if (feof(finput)) {
651 Error(InitReading, "ForcesFile does not contain enough lines.");
652 } else {
653 //fscanf(finput, "%s\n", line);
654 //if (strlen(line) < 100) {
655 //fprintf(stderr, "i %d, j %d, length of line %ld, line %s.\n", i,j, strlen(line), line);
656 //Error(InitReading, "ForcesFile does not contain enough lines or line is faulty.");
657 //} else
658 //fprintf(stderr, "Parsed line: '%s'\n", line);
659 fscanf(finput, "%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n", &i1,&j1,
660 &R[0], &R[1], &R[2],
661 &I->I[i].FIon[0+j*NDIM], &I->I[i].FIon[1+j*NDIM], &I->I[i].FIon[2+j*NDIM],
662 &I->I[i].FIonL[0+j*NDIM], &I->I[i].FIonL[1+j*NDIM], &I->I[i].FIonL[2+j*NDIM],
663 &I->I[i].FIonNL[0+j*NDIM], &I->I[i].FIonNL[1+j*NDIM], &I->I[i].FIonNL[2+j*NDIM],
664 &I->I[i].FMagnetic[0+j*NDIM], &I->I[i].FMagnetic[1+j*NDIM], &I->I[i].FMagnetic[2+j*NDIM],
665 &I->I[i].FEwald[0+j*NDIM], &I->I[i].FEwald[1+j*NDIM], &I->I[i].FEwald[2+j*NDIM]);
666 if ((i != (int)i1) || (j != (int)j1))
667 Error(InitReading, "Line does not match to desired ion.");
668 }
669 fclose(finput);
670 //fprintf(stderr, "MeanForce:\t%e\n", Run->MeanForce[0]);
671};
672
673/** Frees memory IonType.
674 * All pointers from IonType and the pointer on it are free'd
675 * \param *I Ions structure to be free'd
676 * \sa IonsInitRead where memory is allocated
677 */
678void RemoveIonsRead(struct Ions *I)
679{
680 int i, it;
681 for (i=0; i < I->Max_Types; i++) {
682 Free(I->I[i].Name, "RemoveIonsRead: I->I[i].Name");
683 Free(I->I[i].Symbol, "RemoveIonsRead: I->I[i].Symbol");
684 for (it=0;it<I->I[i].Max_IonsOfType;it++) {
685 Free(I->I[i].moment[it], "RemoveIonsRead: I->I[i].moment[it]");
686 Free(I->I[i].sigma[it], "RemoveIonsRead: I->I[i].sigma[it]");
687 Free(I->I[i].sigma_rezi[it], "RemoveIonsRead: I->I[i].sigma_rezi[it]");
688 Free(I->I[i].sigma_PAS[it], "RemoveIonsRead: I->I[i].sigma_PAS[it]");
689 Free(I->I[i].sigma_rezi_PAS[it], "RemoveIonsRead: I->I[i].sigma_rezi_PAS[it]");
690 while (RemoveConstraintItem(&I->I[i], it)); // empty the constrained item list
691 }
692 Free(I->I[i].sigma, "RemoveIonsRead: I->I[i].sigma");
693 Free(I->I[i].sigma_rezi, "RemoveIonsRead: I->I[i].sigma_rezi");
694 Free(I->I[i].sigma_PAS, "RemoveIonsRead: I->I[i].sigma_PAS");
695 Free(I->I[i].sigma_rezi_PAS, "RemoveIonsRead: I->I[i].sigma_rezi_PAS");
696 Free(I->I[i].R, "RemoveIonsRead: I->I[i].R");
697 Free(I->I[i].R_old, "RemoveIonsRead: I->I[i].R_old");
698 Free(I->I[i].R_old_old, "RemoveIonsRead: I->I[i].R_old_old");
699 Free(I->I[i].Constraints, "RemoveIonsRead: I->I[i].Constraints");
700 Free(I->I[i].FIon, "RemoveIonsRead: I->I[i].FIon");
701 Free(I->I[i].FIon_old, "RemoveIonsRead: I->I[i].FIon_old");
702 Free(I->I[i].SearchDir, "RemoveIonsRead: I->I[i].SearchDir");
703 Free(I->I[i].GammaA, "RemoveIonsRead: I->I[i].GammaA");
704 Free(I->I[i].FIonL, "RemoveIonsRead: I->I[i].FIonL");
705 Free(I->I[i].FIonNL, "RemoveIonsRead: I->I[i].FIonNL");
706 Free(I->I[i].FMagnetic, "RemoveIonsRead: I->I[i].FMagnetic");
707 Free(I->I[i].U, "RemoveIonsRead: I->I[i].U");
708 Free(I->I[i].SFactor, "RemoveIonsRead: I->I[i].SFactor");
709 Free(I->I[i].IMT, "RemoveIonsRead: I->I[i].IMT");
710 Free(I->I[i].FEwald, "RemoveIonsRead: I->I[i].FEwald");
711 Free(I->I[i].alpha, "RemoveIonsRead: I->I[i].alpha");
712 }
713 if (I->R_cut == 0.0)
714 Free(I->RLatticeVec, "RemoveIonsRead: I->RLatticeVec");
715 Free(I->FTemp, "RemoveIonsRead: I->FTemp");
716 Free(I->I, "RemoveIonsRead: I->I");
717}
718
719/** Shifts center of gravity of ion forces IonType::FIon.
720 * First sums up all forces of movable ions to a "force temperature",
721 * then reduces each force by this temp, so that all in all
722 * the net force is 0.
723 * \param *P Problem at hand
724 * \sa CorrectVelocity
725 * \note why is FTemp not divided by number of ions: is probably correct, but this
726 * function was not really meaningful from the beginning.
727 */
728void CorrectForces(struct Problem *P)
729{
730 struct Ions *I = &P->Ion;
731 double *FIon;
732 double FTemp[NDIM];
733 int is, ia, d;
734 for (d=0; d<NDIM; d++)
735 FTemp[d] = 0;
736 for (is=0; is < I->Max_Types; is++) { // calculate force temperature for each type ...
737 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of this type ...
738 FIon = &I->I[is].FIon[NDIM*ia];
739 if (I->I[is].IMT[ia] == MoveIon) { // .. if it's movable
740 for (d=0; d<NDIM; d++)
741 FTemp[d] += FIon[d];
742 }
743 }
744 }
745 for (is=0; is < I->Max_Types; is++) { // and then reduce each by this value
746 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
747 FIon = &I->I[is].FIon[NDIM*ia];
748 if (I->I[is].IMT[ia] == MoveIon) {
749 for (d=0; d<NDIM; d++)
750 FIon[d] -= FTemp[d]; // (?) why not -= FTemp[d]/I->Max_TotalIons
751 }
752 }
753 }
754}
755
756
757/** Shifts center of gravity of ion velocities (rather momentums).
758 * First sums up ion speed U times their IonMass (summed up for each
759 * dimension), then reduces the velocity by temp per Ions::TotalMass (so here
760 * total number of Ions is included)
761 * \param *P Problem at hand
762 * \sa CorrectForces
763 */
764void CorrectVelocity(struct Problem *P)
765{
766 struct Ions *I = &P->Ion;
767 double *U;
768 double UTemp[NDIM];
769 int is, ia, d;
770 for (d=0; d<NDIM; d++)
771 UTemp[d] = 0;
772 for (is=0; is < I->Max_Types; is++) { // all types ...
773 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // ... all ions per type ...
774 U = &I->I[is].U[NDIM*ia];
775 if (I->I[is].IMT[ia] == MoveIon) { // ... if it's movable
776 for (d=0; d<NDIM; d++)
777 UTemp[d] += I->I[is].IonMass*U[d];
778 }
779 }
780 }
781 for (is=0; is < I->Max_Types; is++) {
782 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
783 U = &I->I[is].U[NDIM*ia];
784 if (I->I[is].IMT[ia] == MoveIon) {
785 for (d=0; d<NDIM; d++)
786 U[d] -= UTemp[d]/I->TotalMass; // shift by mean velocity over mass and number of ions
787 }
788 }
789 }
790}
791
792/** Moves ions according to calculated force.
793 * Goes through each type of IonType and each ion therein, takes mass and
794 * Newton to move each ion to new coordinates IonType::R, while remembering the last
795 * two coordinates and the last force of the ion. Coordinates R are being
796 * transformed to inverted base, shifted by +-1.0 and back-transformed to
797 * real base.
798 * \param *P Problem at hand
799 * \sa CalculateIonForce
800 * \sa UpdateIonsU
801 * \warning U is not changed here, only used to move the ion.
802 */
803void UpdateIonsR(struct Problem *P)
804{
805 struct Ions *I = &P->Ion;
806 int is, ia, d;
807 double *R, *R_old, *R_old_old, *FIon, *FIon_old, *U, *FIonL, *FIonNL, *FEwald;
808 double IonMass, a;
809 double sR[NDIM], sRold[NDIM], sRoldold[NDIM];
810 const int delta_t = P->R.delta_t;
811 for (is=0; is < I->Max_Types; is++) { // go through all types ...
812 IonMass = I->I[is].IonMass;
813 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsR: IonMass = %lg",IonMass);
814 a = delta_t*0.5/IonMass; // F/m = a and thus: s = 0.5 * F/m * t^2 + v * t =: t * (F * a + v)
815 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
816 FIon = &I->I[is].FIon[NDIM*ia];
817 FIonL = &I->I[is].FIonL[NDIM*ia];
818 FIonNL = &I->I[is].FIonNL[NDIM*ia];
819 FEwald = &I->I[is].FEwald[NDIM*ia];
820 FIon_old = &I->I[is].FIon_old[NDIM*ia];
821 U = &I->I[is].U[NDIM*ia];
822 R = &I->I[is].R[NDIM*ia];
823 R_old = &I->I[is].R_old[NDIM*ia];
824 R_old_old = &I->I[is].R_old_old[NDIM*ia];
825 if (I->I[is].IMT[ia] == MoveIon) {
826 for (d=0; d<NDIM; d++) {
827 R_old_old[d] = R_old[d]; // shift old values
828 R_old[d] = R[d];
829 R[d] += delta_t*(U[d]+a*FIon[d]); // F = m * a and s = 0.5 * a * t^2
830 FIon_old[d] = FIon[d]; // store old force
831 FIon[d] = 0; // zero all as a sign that's been moved
832 FIonL[d] = 0;
833 FIonNL[d] = 0;
834 FEwald[d] = 0;
835 }
836 RMat33Vec3(sR, P->Lat.InvBasis, R);
837 RMat33Vec3(sRold, P->Lat.InvBasis, R_old);
838 RMat33Vec3(sRoldold, P->Lat.InvBasis, R_old_old);
839 for (d=0; d <NDIM; d++) {
840 while (sR[d] < 0) {
841 sR[d] += 1.0;
842 sRold[d] += 1.0;
843 sRoldold[d] += 1.0;
844 }
845 while (sR[d] >= 1.0) {
846 sR[d] -= 1.0;
847 sRold[d] -= 1.0;
848 sRoldold[d] -= 1.0;
849 }
850 }
851 RMat33Vec3(R, P->Lat.RealBasis, sR);
852 RMat33Vec3(R_old, P->Lat.RealBasis, sRold);
853 RMat33Vec3(R_old_old, P->Lat.RealBasis, sRoldold);
854 }
855 }
856 }
857}
858
859/** Resets the forces array.
860 * \param *P Problem at hand
861 */
862void ResetForces(struct Problem *P)
863{
864 struct Ions *I = &P->Ion;
865 double *FIon, *FIonL, *FIonNL, *FEwald;
866 int is, ia, d;
867 for (is=0; is < I->Max_Types; is++) { // go through all types ...
868 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
869 FIon = &I->I[is].FIon[NDIM*ia];
870 FIonL = &I->I[is].FIonL[NDIM*ia];
871 FIonNL = &I->I[is].FIonNL[NDIM*ia];
872 FEwald = &I->I[is].FEwald[NDIM*ia];
873 for (d=0; d<NDIM; d++) {
874 FIon[d] = 0.; // zero all as a sign that's been moved
875 FIonL[d] = 0.;
876 FIonNL[d] = 0.;
877 FEwald[d] = 0.;
878 }
879 }
880 }
881};
882
883/** Changes ion's velocity according to acting force.
884 * IonType::U is changed by the weighted force of actual step and last one
885 * according to Newton.
886 * \param *P Problem at hand
887 * \sa UpdateIonsR
888 * \sa CalculateIonForce
889 * \warning R is not changed here.
890 */
891void UpdateIonsU(struct Problem *P)
892{
893 struct Ions *I = &P->Ion;
894 int is, ia, d;
895 double *FIon, *FIon_old, *U;
896 double IonMass, a;
897 const int delta_t = P->R.delta_t;
898 for (is=0; is < I->Max_Types; is++) {
899 IonMass = I->I[is].IonMass;
900 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsU: IonMass = %lg",IonMass);
901 a = delta_t*0.5/IonMass; // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
902 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
903 FIon = &I->I[is].FIon[NDIM*ia];
904 FIon_old = &I->I[is].FIon_old[NDIM*ia];
905 U = &I->I[is].U[NDIM*ia];
906 if (I->I[is].IMT[ia] == MoveIon)
907 for (d=0; d<NDIM; d++) {
908 U[d] += a * (FIon[d]+FIon_old[d]);
909 }
910 }
911 }
912}
913
914/** CG process to optimise structure.
915 * \param *P Problem at hand
916 */
917void UpdateIons(struct Problem *P)
918{
919 struct Ions *I = &P->Ion;
920 int is, ia, d;
921 double *R, *R_old, *R_old_old, *FIon, *S, *GammaA;
922 double IonFac, GammaB; //, GammaT;
923 double FNorm, StepNorm;
924 for (is=0; is < I->Max_Types; is++) {
925 IonFac = I->I[is].IonFac;
926 GammaA = I->I[is].GammaA;
927 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
928 FIon = &I->I[is].FIon[NDIM*ia];
929 S = &I->I[is].SearchDir[NDIM*ia];
930 GammaB = GammaA[ia];
931 GammaA[ia] = RSP3(FIon,S);
932 FNorm = sqrt(RSP3(FIon,FIon));
933 StepNorm = log(1.+IonFac*FNorm); // Fix Hack
934 if (FNorm != 0)
935 StepNorm /= sqrt(RSP3(FIon,FIon));
936 else
937 StepNorm = 0;
938 //if (GammaB != 0.0) {
939 // GammaT = GammaA[ia]/GammaB;
940 // for (d=0; d>NDIM; d++)
941 // S[d] = FIon[d]+GammaT*S[d];
942 // } else {
943 for (d=0; d<NDIM; d++)
944 S[d] = StepNorm*FIon[d];
945 // }
946 R = &I->I[is].R[NDIM*ia];
947 R_old = &I->I[is].R_old[NDIM*ia];
948 R_old_old = &I->I[is].R_old_old[NDIM*ia];
949 if (I->I[is].IMT[ia] == MoveIon)
950 for (d=0; d<NDIM; d++) {
951 R_old_old[d] = R_old[d];
952 R_old[d] = R[d];
953 R[d] += S[d]; // FixHack*IonFac;
954 }
955 }
956 }
957}
958
959/** Print coordinates of all ions to stdout.
960 * \param *P Problem at hand
961 * \param actual (1 - don't at current RunStruct#OuterStep, 0 - do)
962 */
963void OutputIonCoordinates(struct Problem *P, int actual)
964{
965 //struct RunStruct *R = &P->R;
966 struct Ions *I = &P->Ion;
967 char filename[255];
968 FILE *output;
969 int is, ia, nr = 0;
970 double Bohr = (I->IsAngstroem) ? 1./ANGSTROEMINBOHRRADIUS : 1.;
971/* if (P->Par.me == 0 && P->Call.out[ReadOut]) {
972 fprintf(stderr, "(%i) ======= Differential Ion Coordinates =========\n",P->Par.me);
973 for (is=0; is < I->Max_Types; is++)
974 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++)
975 //fprintf(stderr, "(%i) R[%i/%i][%i/%i] = (%e,%e,%e)\n", P->Par.me, is, I->Max_Types, ia, I->I[is].Max_IonsOfType, I->I[is].R[NDIM*ia+0],I->I[is].R[NDIM*ia+1],I->I[is].R[NDIM*ia+2]);
976 fprintf(stderr, "Ion_Type%i_%i\t%.6f\t%.6f\t%.6f\t0\t# Atom from StructOpt\n", is+1, ia+1, I->I[is].R[NDIM*ia+0]-I->I[is].R_old[NDIM*ia+0],I->I[is].R[NDIM*ia+1]-I->I[is].R_old[NDIM*ia+1],I->I[is].R[NDIM*ia+2]-I->I[is].R_old[NDIM*ia+2]);
977 fprintf(stderr, "(%i) =========================================\n",P->Par.me);
978 }*/
979 if (actual)
980 sprintf(filename, "%s.MD", P->Call.MainParameterFile);
981 else
982 sprintf(filename, "%s.opt", P->Call.MainParameterFile);
983 if (((actual) && (P->R.OuterStep <= 0)) || ((!actual) && (P->R.StructOptStep == 1))) { // on first step write complete config file
984 WriteParameters(P, filename);
985 } else { // afterwards simply add the current ion coordinates as constrained steps
986 if ((P->Par.me == 0) && (output = fopen(filename,"a"))) {
987 fprintf(output,"# ====== MD step %d ========= \n", (actual) ? P->R.OuterStep : P->R.StructOptStep);
988 for (is=0; is < I->Max_Types; is++)
989 for (ia=0;ia<I->I[is].Max_IonsOfType;ia++) {
990 fprintf(output,"Ion_Type%i_%i\t%2.9f\t%2.9f\t%2.9f\t%i\t%2.9f\t%2.9f\t%2.9f\t# Number in molecule %i\n", is+1, ia+1, I->I[is].R[0+NDIM*ia]*Bohr, I->I[is].R[1+NDIM*ia]*Bohr, I->I[is].R[2+NDIM*ia]*Bohr, I->I[is].IMT[ia], I->I[is].U[0+NDIM*ia]*Bohr, I->I[is].U[1+NDIM*ia]*Bohr, I->I[is].U[2+NDIM*ia]*Bohr, ++nr);
991 }
992 fflush(output);
993 }
994 }
995}
996
997/** Calculates kinetic energy (Ions::EKin) of movable Ions.
998 * Does 0.5 / IonType::IonMass * IonTye::U^2 for each Ion,
999 * also retrieves actual temperatur by 2/3 from just
1000 * calculated Ions::EKin.
1001 * \param *P Problem at hand
1002 */
1003void CalculateEnergyIonsU(struct Problem *P)
1004{
1005 struct Ions *I = &P->Ion;
1006 int is, ia, d;
1007 double *U;
1008 double IonMass, a, ekin = 0;
1009 for (is=0; is < I->Max_Types; is++) {
1010 IonMass = I->I[is].IonMass;
1011 a = 0.5*IonMass;
1012 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1013 U = &I->I[is].U[NDIM*ia];
1014 //if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
1015 for (d=0; d<NDIM; d++) {
1016 ekin += a * U[d]*U[d];
1017 }
1018 }
1019 }
1020 I->EKin = ekin;
1021 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
1022}
1023
1024/** Scales ion velocities to match temperature.
1025 * In order to match Ions::ActualTemp with desired RunStruct::TargetTemp
1026 * each velocity in each dimension (for all types, all ions) is scaled
1027 * with the ratio of these two. Ions::EKin is recalculated.
1028 * \param *P Problem at hand
1029 */
1030void ScaleTemp(struct Problem *P)
1031{
1032 struct Ions *I = &P->Ion;
1033 int is, ia, d;
1034 double *U;
1035 double IonMass, a, ekin = 0;
1036 if (I->ActualTemp < MYEPSILON) fprintf(stderr,"ScaleTemp: I->ActualTemp = %lg",I->ActualTemp);
1037 double ScaleTempFactor = sqrt(P->R.TargetTemp/I->ActualTemp);
1038 for (is=0; is < I->Max_Types; is++) {
1039 IonMass = I->I[is].IonMass;
1040 a = 0.5*IonMass;
1041 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1042 U = &I->I[is].U[NDIM*ia];
1043 //if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
1044 for (d=0; d<NDIM; d++) {
1045 U[d] *= ScaleTempFactor;
1046 ekin += a * U[d]*U[d];
1047 }
1048 }
1049 }
1050 I->EKin = ekin;
1051 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
1052}
1053
1054/** Calculates mean force vector as stop criteria in structure optimization.
1055 * Calculates a mean force vector per ion as the usual euclidian norm over total
1056 * number of ions and dimensions, being the sum of each ion force (for all type,
1057 * all ions, all dims) squared.
1058 * The mean force is archived in RunStruct::MeanForce and printed to screen.
1059 * \param *P Problem at hand
1060 */
1061void GetOuterStop(struct Problem *P)
1062{
1063 struct RunStruct *R = &P->R;
1064 struct Ions *I = &P->Ion;
1065 int is, ia, IonNo=0, i, d;
1066 double MeanForce = 0.0;
1067 for (is=0; is < I->Max_Types; is++)
1068 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1069 IonNo++;
1070 for (d=0; d <NDIM; d++)
1071 MeanForce += I->I[is].FIon[d+NDIM*ia]*I->I[is].FIon[d+NDIM*ia];
1072 }
1073 for (i=MAXOLD-1; i > 0; i--)
1074 R->MeanForce[i] = R->MeanForce[i-1];
1075 MeanForce = sqrt(MeanForce/(IonNo*NDIM));
1076 R->MeanForce[0] = MeanForce;
1077 if (P->Call.out[LeaderOut] && (P->Par.me == 0))
1078 fprintf(stderr,"MeanForce = %e\n", R->MeanForce[0]);
1079}
1080
1081
Note: See TracBrowser for help on using the repository browser.