source: pcp/src/ions.c@ 69eca8

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

CalculateIonForce(): some change
  • Property mode set to 100644
File size: 47.5 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 * The basic idea of the ewald method is to add a countercharge - here of gaussian form -
441 * which shields the long-range charges making them effectively short-ranged, while the
442 * countercharges themselves can be analytically (here by (hidden) Fouriertransform: it's
443 * added to the electron density) integrated and withdrawn.
444 * \f[
445 * 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|}
446 * \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 )
447 * \qquad (4.10)
448 * \f]
449 * Calculates the core-to-core force between the ions of all super cells.
450 * In order for this series to converge there must be a certain summation
451 * applied, which is the ewald summation and which is nothing else but to move
452 * in a circular motion from the current cell to the outside up to Ions::R_cut.
453 * \param *P Problem at hand
454 * \param first if != 0 table mirror cells to take into (local!) account of ewald summation
455 */
456void CalculateEwald(struct Problem *P, int first)
457{
458 int r,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec; //,is
459 int MaxPar=P->Par.procs, ParMe=P->Par.me, MaxLocalVec, StartVec;
460 double rcut2,r2,erre2,addesr,addpre,arg,fac,gkl,esr,fxx,repand;
461 double R1[3];
462 double R2[3];
463 double R3[3];
464 double t[3];
465 struct Lattice *L = &P->Lat;
466 struct PseudoPot *PP = &P->PP;
467 struct Ions *I = &P->Ion;
468 if (first) {
469 SpeedMeasure(P, InitSimTime, StopTimeDo);
470 rcut2 = I->R_cut*I->R_cut;
471 ActualVec = 0; // counts number of cells taken into account
472 MaxCell = (int)(5.*I->R_cut/pow(L->Volume,1./3.)); // the number of cells in each direction is prop. to axis length over cutoff
473 if (MaxCell < 2) MaxCell = 2; // take at least 2
474 for (i = -MaxCell; i <= MaxCell; i++) {
475 for (j = -MaxCell; j <= MaxCell; j++) {
476 for (k = -MaxCell; k <= MaxCell; k++) {
477 r2 = 0.;
478 for (ir=0; ir <NDIM; ir++) { // t is the offset vector pointing to distant cell (only diagonal entries)
479 t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
480 r2 = t[ir]*t[ir]; // is squared distance to other cell
481 }
482 Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1)); // whether it's directly adjacent
483 if ((r2 <= rcut2) || Is_Neighb_Cell) { // if it's either adjacent or within cutoff, count it in
484 ActualVec++;
485 }
486 }
487 }
488 }
489 MaxVec = ActualVec; // store number of cells
490 I->MaxVec = ActualVec;
491 MaxLocalVec = MaxVec / MaxPar; // share among processes
492 StartVec = ParMe*MaxLocalVec; // for this process start at its patch
493 r = MaxVec % MaxPar; // rest not fitting...
494 if (ParMe >= r) { // ones who don't get extra cells have to shift their start vector
495 StartVec += r;
496 } else { // others get extra cell and have to shift also by a bit
497 StartVec += ParMe;
498 MaxLocalVec++;
499 }
500 I->MaxLocalVec = MaxLocalVec;
501 LocalActualVec = ActualVec = 0;
502 // now go through the same loop again and note down every offset vector for cells in this process' patch
503 if (MaxLocalVec != 0) {
504 I->RLatticeVec = (double *)
505 Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald: I->RLatticeVec");
506 } else {
507 fprintf(stderr,"(%i) Warning in Ewald summation: Got MaxLocalVec == 0\n", P->Par.me);
508 }
509 for (i = -MaxCell; i <= MaxCell; i++) {
510 for (j = -MaxCell; j <= MaxCell; j++) {
511 for (k = -MaxCell; k <= MaxCell; k++) {
512 r2 = 0.;
513 for (ir=0; ir <NDIM; ir++) { // create offset vector from diagonal entries
514 t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
515 r2 = t[ir]*t[ir];
516 }
517 Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
518 if ((r2 <= rcut2) || Is_Neighb_Cell) { // if adjacent or within cutoff ...
519 if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) { // ... and belongs to patch, store 'em
520 I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
521 I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
522 I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
523 LocalActualVec++;
524 }
525 ActualVec++;
526 }
527 }
528 }
529 }
530 SpeedMeasure(P, InitSimTime, StartTimeDo);
531 }
532 SpeedMeasure(P, EwaldTime, StartTimeDo);
533
534 // first set everything to zero
535 esr = 0.0;
536 //for (is1=0;is1 < I->Max_Types; is1++)
537 // then take each ion of each type once
538 for (is1=0;is1 < I->Max_Types; is1++) {
539 // reset FTemp for IonType
540 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
541 for (i=0; i< NDIM; i++)
542 I->FTemp[i+NDIM*(ia1)] = 0.;
543 // then sum for each on of the IonType
544 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
545 for (i=0;i<NDIM;i++) {
546 R1[i]=I->I[is1].R[i+NDIM*ia1]; // R1 is local coordinate of current first ion
547 }
548 for (ir2=0;ir2<I->MaxLocalVec;ir2++) { // for every cell of local patch (each with the same atoms)
549 for (is2=0;is2 < I->Max_Types; is2++) { // and for each ion of each type a second time
550 // gkl = 1./(sqrt(r_is1,gauss^2 + r_is2,gauss^2))
551 gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
552 for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
553 for (i=0;i<3;i++) {
554 R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2]; // R2 is coordinate of current second ion in the distant cell!
555 }
556
557 erre2=0.0;
558 for (i=0;i<NDIM;i++) {
559 R3[i] = R1[i]-R2[i]; // R3 = R_is1,ia1 - R_is2,ia2 - L
560 erre2+=R3[i]*R3[i]; // erre2 = |R_is1,ia1 - R_is2,ia2 - L|^2
561 }
562 if (erre2 > MYEPSILON) { // appears as one over ... in fac, thus check if zero
563 arg = sqrt(erre2); // arg = |R_is1,ia1 - R_is2,ia2 - L|
564 fac=PP->zval[is1]*PP->zval[is2]/arg*0.5; // fac = -1/2 * Z_is1 * Z_is2 / arg
565
566 arg *= gkl; // arg = |R_is1,ia1 - R_is2,ia2 - L|/(sqrt(r_is1,gauss^2 + r_is2,gauss^2))
567 addesr = derf(arg);
568 addesr = (1.-addesr)*fac; // complementary error function: addesr = 1/2*erfc(arg)/|R_is1,ia1 - R_is2,ia2 - L|
569 esr += addesr;
570 addpre=exp(-arg*arg)*gkl;
571 addpre = PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre; // addpre = exp(-|R_is1,ia1 - R_is2,ia2 - L|^2/(r_is1,gauss^2 + r_is2,gauss^2))/(sqrt(pi)*sqrt(r_is1,gauss^2 + r_is2,gauss^2))
572 repand = (addesr+addpre)/erre2;
573// if ((fabs(repand) > MYEPSILON)) {
574// fprintf(stderr,"(%i) Ewald[is%d,ia%d/is%d,ia%d,(%d)]: Z1 %lg, Z2 %lg\tpre %lg\t esr %lg\t fac %lg\t arg %lg\tR3 (%lg,%lg,%lg)\n", P->Par.me, is1, ia1, is2, ia2, ir2, PP->zval[is1], PP->zval[is2],addpre,addesr, fac, arg, R3[0], R3[1], R3[2]);
575// }
576 for (i=0;i<NDIM;i++) {
577 fxx=repand*R3[i];
578// if ((fabs(repand) > MYEPSILON)) {
579// fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);
580// }
581 I->FTemp[i+NDIM*(ia1)] += fxx*2.0; // actio = reactio?
582 //I->FTemp[i+NDIM*(ia2)] -= fxx;
583 }
584 }
585 }
586 }
587 }
588 }
589 MPI_Allreduce (I->FTemp , I->I[is1].FEwald, NDIM*I->I[is1].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);
590 }
591 SpeedMeasure(P, EwaldTime, StopTimeDo);
592 //for (is=0;is < I->Max_Types; is++) {
593 //}
594 MPI_Allreduce ( &esr, &L->E->AllTotalIonsEnergy[EwaldEnergy], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm);
595// fprintf(stderr,"\n");
596}
597
598/** Sum up all forces acting on Ions.
599 * IonType::FIon = IonType::FEwald + IonType::FIonL + IonType::FIonNL + IonType::FMagnetic for all
600 * dimensions #NDIM and each Ion of IonType
601 * \param *P Problem at hand
602 */
603void CalculateIonForce(struct Problem *P)
604{
605 struct Ions *I = &P->Ion;
606 int i,j,d;
607 for (i=0; i < I->Max_Types; i++)
608 for (j=0; j < I->I[i].Max_IonsOfType; j++)
609 for (d=0; d<NDIM;d++)
610 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] + I->I[i].FMagnetic[d+j*NDIM]);
611}
612
613/** Write Forces to FileData::ForcesFile.
614 * goes through all Ions per IonType and each dimension of #NDIM and prints in one line:
615 * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL,
616 * non-local IonType::FIonNL, ewald force IonType::FEwald
617 * \param *P Problem at hand
618 */
619void OutputIonForce(struct Problem *P)
620{
621 struct RunStruct *R = &P->R;
622 struct FileData *F = &P->Files;
623 struct Ions *I = &P->Ion;
624 FILE *fout = NULL;
625 int i,j;
626 if (F->MeOutMes != 1) return;
627 fout = F->ForcesFile;
628 for (i=0; i < I->Max_Types; i++)
629 for (j=0; j < I->I[i].Max_IonsOfType; j++)
630 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,
631 I->I[i].R[0+j*NDIM], I->I[i].R[1+j*NDIM], I->I[i].R[2+j*NDIM],
632 I->I[i].FIon[0+j*NDIM], I->I[i].FIon[1+j*NDIM], I->I[i].FIon[2+j*NDIM],
633 I->I[i].FIonL[0+j*NDIM], I->I[i].FIonL[1+j*NDIM], I->I[i].FIonL[2+j*NDIM],
634 I->I[i].FIonNL[0+j*NDIM], I->I[i].FIonNL[1+j*NDIM], I->I[i].FIonNL[2+j*NDIM],
635 I->I[i].FEwald[0+j*NDIM], I->I[i].FEwald[1+j*NDIM], I->I[i].FEwald[2+j*NDIM]);
636 fflush(fout);
637 fprintf(fout, "MeanForce:\t%e\n", R->MeanForce[0]);
638}
639
640/** Reads Forces from CallOptions::ForcesFile.
641 * goes through all Ions per IonType and each dimension of #NDIM and parses in one line:
642 * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL,
643 * non-local IonType::FIonNL, ewald force IonType::FEwald
644 * \param *P Problem at hand
645 */
646void ParseIonForce(struct Problem *P)
647{
648 //struct RunStruct *Run = &P->R;
649 struct Ions *I = &P->Ion;
650 FILE *finput = NULL;
651 char line[1024];
652 int i,j;
653 double i1, j1;
654 double R[NDIM];
655
656 if (P->Call.ForcesFile == NULL) return;
657 finput = fopen(P->Call.ForcesFile, "r");
658 //fprintf(stderr, "File pointer says %p\n", finput);
659 if (finput == NULL)
660 Error(InitReading, "ForcesFile does not exist.");
661 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
662 for (i=0; i < I->Max_Types; i++)
663 for (j=0; j < I->I[i].Max_IonsOfType; j++)
664 if (feof(finput)) {
665 Error(InitReading, "ForcesFile does not contain enough lines.");
666 } else {
667 //fscanf(finput, "%s\n", line);
668 //if (strlen(line) < 100) {
669 //fprintf(stderr, "i %d, j %d, length of line %ld, line %s.\n", i,j, strlen(line), line);
670 //Error(InitReading, "ForcesFile does not contain enough lines or line is faulty.");
671 //} else
672 //fprintf(stderr, "Parsed line: '%s'\n", line);
673 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,
674 &R[0], &R[1], &R[2],
675 &I->I[i].FIon[0+j*NDIM], &I->I[i].FIon[1+j*NDIM], &I->I[i].FIon[2+j*NDIM],
676 &I->I[i].FIonL[0+j*NDIM], &I->I[i].FIonL[1+j*NDIM], &I->I[i].FIonL[2+j*NDIM],
677 &I->I[i].FIonNL[0+j*NDIM], &I->I[i].FIonNL[1+j*NDIM], &I->I[i].FIonNL[2+j*NDIM],
678 &I->I[i].FMagnetic[0+j*NDIM], &I->I[i].FMagnetic[1+j*NDIM], &I->I[i].FMagnetic[2+j*NDIM],
679 &I->I[i].FEwald[0+j*NDIM], &I->I[i].FEwald[1+j*NDIM], &I->I[i].FEwald[2+j*NDIM]);
680 if ((i != (int)i1) || (j != (int)j1))
681 Error(InitReading, "Line does not match to desired ion.");
682 }
683 fclose(finput);
684 //fprintf(stderr, "MeanForce:\t%e\n", Run->MeanForce[0]);
685};
686
687/** Frees memory IonType.
688 * All pointers from IonType and the pointer on it are free'd
689 * \param *I Ions structure to be free'd
690 * \sa IonsInitRead where memory is allocated
691 */
692void RemoveIonsRead(struct Ions *I)
693{
694 int i, it;
695 for (i=0; i < I->Max_Types; i++) {
696 Free(I->I[i].Name, "RemoveIonsRead: I->I[i].Name");
697 Free(I->I[i].Symbol, "RemoveIonsRead: I->I[i].Symbol");
698 for (it=0;it<I->I[i].Max_IonsOfType;it++) {
699 Free(I->I[i].moment[it], "RemoveIonsRead: I->I[i].moment[it]");
700 Free(I->I[i].sigma[it], "RemoveIonsRead: I->I[i].sigma[it]");
701 Free(I->I[i].sigma_rezi[it], "RemoveIonsRead: I->I[i].sigma_rezi[it]");
702 Free(I->I[i].sigma_PAS[it], "RemoveIonsRead: I->I[i].sigma_PAS[it]");
703 Free(I->I[i].sigma_rezi_PAS[it], "RemoveIonsRead: I->I[i].sigma_rezi_PAS[it]");
704 while (RemoveConstraintItem(&I->I[i], it)); // empty the constrained item list
705 }
706 Free(I->I[i].sigma, "RemoveIonsRead: I->I[i].sigma");
707 Free(I->I[i].sigma_rezi, "RemoveIonsRead: I->I[i].sigma_rezi");
708 Free(I->I[i].sigma_PAS, "RemoveIonsRead: I->I[i].sigma_PAS");
709 Free(I->I[i].sigma_rezi_PAS, "RemoveIonsRead: I->I[i].sigma_rezi_PAS");
710 Free(I->I[i].R, "RemoveIonsRead: I->I[i].R");
711 Free(I->I[i].R_old, "RemoveIonsRead: I->I[i].R_old");
712 Free(I->I[i].R_old_old, "RemoveIonsRead: I->I[i].R_old_old");
713 Free(I->I[i].Constraints, "RemoveIonsRead: I->I[i].Constraints");
714 Free(I->I[i].FIon, "RemoveIonsRead: I->I[i].FIon");
715 Free(I->I[i].FIon_old, "RemoveIonsRead: I->I[i].FIon_old");
716 Free(I->I[i].SearchDir, "RemoveIonsRead: I->I[i].SearchDir");
717 Free(I->I[i].GammaA, "RemoveIonsRead: I->I[i].GammaA");
718 Free(I->I[i].FIonL, "RemoveIonsRead: I->I[i].FIonL");
719 Free(I->I[i].FIonNL, "RemoveIonsRead: I->I[i].FIonNL");
720 Free(I->I[i].FMagnetic, "RemoveIonsRead: I->I[i].FMagnetic");
721 Free(I->I[i].U, "RemoveIonsRead: I->I[i].U");
722 Free(I->I[i].SFactor, "RemoveIonsRead: I->I[i].SFactor");
723 Free(I->I[i].IMT, "RemoveIonsRead: I->I[i].IMT");
724 Free(I->I[i].FEwald, "RemoveIonsRead: I->I[i].FEwald");
725 Free(I->I[i].alpha, "RemoveIonsRead: I->I[i].alpha");
726 }
727 if (I->R_cut == 0.0)
728 Free(I->RLatticeVec, "RemoveIonsRead: I->RLatticeVec");
729 Free(I->FTemp, "RemoveIonsRead: I->FTemp");
730 Free(I->I, "RemoveIonsRead: I->I");
731}
732
733/** Shifts center of gravity of ion forces IonType::FIon.
734 * First sums up all forces of movable ions to a "force temperature",
735 * then reduces each force by this temp, so that all in all
736 * the net force is 0.
737 * \param *P Problem at hand
738 * \sa CorrectVelocity
739 * \note why is FTemp not divided by number of ions: is probably correct, but this
740 * function was not really meaningful from the beginning.
741 */
742void CorrectForces(struct Problem *P)
743{
744 struct Ions *I = &P->Ion;
745 double *FIon;
746 double FTemp[NDIM];
747 int is, ia, d;
748 for (d=0; d<NDIM; d++)
749 FTemp[d] = 0;
750 for (is=0; is < I->Max_Types; is++) { // calculate force temperature for each type ...
751 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of this type ...
752 FIon = &I->I[is].FIon[NDIM*ia];
753 if (I->I[is].IMT[ia] == MoveIon) { // .. if it's movable
754 for (d=0; d<NDIM; d++)
755 FTemp[d] += FIon[d];
756 }
757 }
758 }
759 for (is=0; is < I->Max_Types; is++) { // and then reduce each by this value
760 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
761 FIon = &I->I[is].FIon[NDIM*ia];
762 if (I->I[is].IMT[ia] == MoveIon) {
763 for (d=0; d<NDIM; d++)
764 FIon[d] -= FTemp[d]; // (?) why not -= FTemp[d]/I->Max_TotalIons
765 }
766 }
767 }
768}
769
770
771/** Shifts center of gravity of ion velocities (rather momentums).
772 * First sums up ion speed U times their IonMass (summed up for each
773 * dimension), then reduces the velocity by temp per Ions::TotalMass (so here
774 * total number of Ions is included)
775 * \param *P Problem at hand
776 * \sa CorrectForces
777 */
778void CorrectVelocity(struct Problem *P)
779{
780 struct Ions *I = &P->Ion;
781 double *U;
782 double UTemp[NDIM];
783 int is, ia, d;
784 for (d=0; d<NDIM; d++)
785 UTemp[d] = 0;
786 for (is=0; is < I->Max_Types; is++) { // all types ...
787 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // ... all ions per type ...
788 U = &I->I[is].U[NDIM*ia];
789 if (I->I[is].IMT[ia] == MoveIon) { // ... if it's movable
790 for (d=0; d<NDIM; d++)
791 UTemp[d] += I->I[is].IonMass*U[d];
792 }
793 }
794 }
795 for (is=0; is < I->Max_Types; is++) {
796 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
797 U = &I->I[is].U[NDIM*ia];
798 if (I->I[is].IMT[ia] == MoveIon) {
799 for (d=0; d<NDIM; d++)
800 U[d] -= UTemp[d]/I->TotalMass; // shift by mean velocity over mass and number of ions
801 }
802 }
803 }
804}
805
806/** Moves ions according to calculated force.
807 * Goes through each type of IonType and each ion therein, takes mass and
808 * Newton to move each ion to new coordinates IonType::R, while remembering the last
809 * two coordinates and the last force of the ion. Coordinates R are being
810 * transformed to inverted base, shifted by +-1.0 and back-transformed to
811 * real base.
812 * \param *P Problem at hand
813 * \sa CalculateIonForce
814 * \sa UpdateIonsU
815 * \warning U is not changed here, only used to move the ion.
816 */
817void UpdateIonsR(struct Problem *P)
818{
819 struct Ions *I = &P->Ion;
820 int is, ia, d;
821 double *R, *R_old, *R_old_old, *FIon, *FIon_old, *U, *FIonL, *FIonNL, *FEwald;
822 double IonMass, a;
823 double sR[NDIM], sRold[NDIM], sRoldold[NDIM];
824 const int delta_t = P->R.delta_t;
825 for (is=0; is < I->Max_Types; is++) { // go through all types ...
826 IonMass = I->I[is].IonMass;
827 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsR: IonMass = %lg",IonMass);
828 a = delta_t*0.5/IonMass; // F/m = a and thus: s = 0.5 * F/m * t^2 + v * t =: t * (F * a + v)
829 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
830 FIon = &I->I[is].FIon[NDIM*ia];
831 FIonL = &I->I[is].FIonL[NDIM*ia];
832 FIonNL = &I->I[is].FIonNL[NDIM*ia];
833 FEwald = &I->I[is].FEwald[NDIM*ia];
834 FIon_old = &I->I[is].FIon_old[NDIM*ia];
835 U = &I->I[is].U[NDIM*ia];
836 R = &I->I[is].R[NDIM*ia];
837 R_old = &I->I[is].R_old[NDIM*ia];
838 R_old_old = &I->I[is].R_old_old[NDIM*ia];
839 if (I->I[is].IMT[ia] == MoveIon) {
840 for (d=0; d<NDIM; d++) {
841 R_old_old[d] = R_old[d]; // shift old values
842 R_old[d] = R[d];
843 R[d] += delta_t*(U[d]+a*FIon[d]); // F = m * a and s = 0.5 * a * t^2
844 FIon_old[d] = FIon[d]; // store old force
845 FIon[d] = 0; // zero all as a sign that's been moved
846 FIonL[d] = 0;
847 FIonNL[d] = 0;
848 FEwald[d] = 0;
849 }
850 RMat33Vec3(sR, P->Lat.InvBasis, R);
851 RMat33Vec3(sRold, P->Lat.InvBasis, R_old);
852 RMat33Vec3(sRoldold, P->Lat.InvBasis, R_old_old);
853 for (d=0; d <NDIM; d++) {
854 while (sR[d] < 0) {
855 sR[d] += 1.0;
856 sRold[d] += 1.0;
857 sRoldold[d] += 1.0;
858 }
859 while (sR[d] >= 1.0) {
860 sR[d] -= 1.0;
861 sRold[d] -= 1.0;
862 sRoldold[d] -= 1.0;
863 }
864 }
865 RMat33Vec3(R, P->Lat.RealBasis, sR);
866 RMat33Vec3(R_old, P->Lat.RealBasis, sRold);
867 RMat33Vec3(R_old_old, P->Lat.RealBasis, sRoldold);
868 }
869 }
870 }
871}
872
873/** Resets the forces array.
874 * \param *P Problem at hand
875 */
876void ResetForces(struct Problem *P)
877{
878 struct Ions *I = &P->Ion;
879 double *FIon, *FIonL, *FIonNL, *FEwald;
880 int is, ia, d;
881 for (is=0; is < I->Max_Types; is++) { // go through all types ...
882 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
883 FIon = &I->I[is].FIon[NDIM*ia];
884 FIonL = &I->I[is].FIonL[NDIM*ia];
885 FIonNL = &I->I[is].FIonNL[NDIM*ia];
886 FEwald = &I->I[is].FEwald[NDIM*ia];
887 for (d=0; d<NDIM; d++) {
888 FIon[d] = 0.; // zero all as a sign that's been moved
889 FIonL[d] = 0.;
890 FIonNL[d] = 0.;
891 FEwald[d] = 0.;
892 }
893 }
894 }
895};
896
897/** Changes ion's velocity according to acting force.
898 * IonType::U is changed by the weighted force of actual step and last one
899 * according to Newton.
900 * \param *P Problem at hand
901 * \sa UpdateIonsR
902 * \sa CalculateIonForce
903 * \warning R is not changed here.
904 */
905void UpdateIonsU(struct Problem *P)
906{
907 struct Ions *I = &P->Ion;
908 int is, ia, d;
909 double *FIon, *FIon_old, *U;
910 double IonMass, a;
911 const int delta_t = P->R.delta_t;
912 for (is=0; is < I->Max_Types; is++) {
913 IonMass = I->I[is].IonMass;
914 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsU: IonMass = %lg",IonMass);
915 a = delta_t*0.5/IonMass; // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
916 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
917 FIon = &I->I[is].FIon[NDIM*ia];
918 FIon_old = &I->I[is].FIon_old[NDIM*ia];
919 U = &I->I[is].U[NDIM*ia];
920 if (I->I[is].IMT[ia] == MoveIon)
921 for (d=0; d<NDIM; d++) {
922 U[d] += a * (FIon[d]+FIon_old[d]);
923 }
924 }
925 }
926}
927
928/** CG process to optimise structure.
929 * \param *P Problem at hand
930 */
931void UpdateIons(struct Problem *P)
932{
933 struct Ions *I = &P->Ion;
934 int is, ia, d;
935 double *R, *R_old, *R_old_old, *FIon, *S, *GammaA;
936 double IonFac, GammaB; //, GammaT;
937 double FNorm, StepNorm;
938 for (is=0; is < I->Max_Types; is++) {
939 IonFac = I->I[is].IonFac;
940 GammaA = I->I[is].GammaA;
941 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
942 FIon = &I->I[is].FIon[NDIM*ia];
943 S = &I->I[is].SearchDir[NDIM*ia];
944 GammaB = GammaA[ia];
945 GammaA[ia] = RSP3(FIon,S);
946 FNorm = sqrt(RSP3(FIon,FIon));
947 StepNorm = log(1.+IonFac*FNorm); // Fix Hack
948 if (FNorm != 0)
949 StepNorm /= sqrt(RSP3(FIon,FIon));
950 else
951 StepNorm = 0;
952 //if (GammaB != 0.0) {
953 // GammaT = GammaA[ia]/GammaB;
954 // for (d=0; d>NDIM; d++)
955 // S[d] = FIon[d]+GammaT*S[d];
956 // } else {
957 for (d=0; d<NDIM; d++)
958 S[d] = StepNorm*FIon[d];
959 // }
960 R = &I->I[is].R[NDIM*ia];
961 R_old = &I->I[is].R_old[NDIM*ia];
962 R_old_old = &I->I[is].R_old_old[NDIM*ia];
963 if (I->I[is].IMT[ia] == MoveIon)
964 for (d=0; d<NDIM; d++) {
965 R_old_old[d] = R_old[d];
966 R_old[d] = R[d];
967 R[d] += S[d]; // FixHack*IonFac;
968 }
969 }
970 }
971}
972
973/** Print coordinates of all ions to stdout.
974 * \param *P Problem at hand
975 * \param actual (1 - don't at current RunStruct#OuterStep, 0 - do)
976 */
977void OutputIonCoordinates(struct Problem *P, int actual)
978{
979 //struct RunStruct *R = &P->R;
980 struct Ions *I = &P->Ion;
981 char filename[255];
982 FILE *output;
983 int is, ia, nr = 0;
984 double Bohr = (I->IsAngstroem) ? 1./ANGSTROEMINBOHRRADIUS : 1.;
985/* if (P->Par.me == 0 && P->Call.out[ReadOut]) {
986 fprintf(stderr, "(%i) ======= Differential Ion Coordinates =========\n",P->Par.me);
987 for (is=0; is < I->Max_Types; is++)
988 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++)
989 //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]);
990 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]);
991 fprintf(stderr, "(%i) =========================================\n",P->Par.me);
992 }*/
993 if (actual)
994 sprintf(filename, "%s.MD", P->Call.MainParameterFile);
995 else
996 sprintf(filename, "%s.opt", P->Call.MainParameterFile);
997 if (((actual) && (P->R.OuterStep <= 0)) || ((!actual) && (P->R.StructOptStep == 1))) { // on first step write complete config file
998 WriteParameters(P, filename);
999 } else { // afterwards simply add the current ion coordinates as constrained steps
1000 if ((P->Par.me == 0) && (output = fopen(filename,"a"))) {
1001 fprintf(output,"# ====== MD step %d ========= \n", (actual) ? P->R.OuterStep : P->R.StructOptStep);
1002 for (is=0; is < I->Max_Types; is++)
1003 for (ia=0;ia<I->I[is].Max_IonsOfType;ia++) {
1004 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);
1005 }
1006 fflush(output);
1007 }
1008 }
1009}
1010
1011/** Calculates kinetic energy (Ions::EKin) of movable Ions.
1012 * Does 0.5 / IonType::IonMass * IonTye::U^2 for each Ion,
1013 * also retrieves actual temperatur by 2/3 from just
1014 * calculated Ions::EKin.
1015 * \param *P Problem at hand
1016 */
1017void CalculateEnergyIonsU(struct Problem *P)
1018{
1019 struct Ions *I = &P->Ion;
1020 int is, ia, d;
1021 double *U;
1022 double IonMass, a, ekin = 0;
1023 for (is=0; is < I->Max_Types; is++) {
1024 IonMass = I->I[is].IonMass;
1025 a = 0.5*IonMass;
1026 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1027 U = &I->I[is].U[NDIM*ia];
1028 //if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
1029 for (d=0; d<NDIM; d++) {
1030 ekin += a * U[d]*U[d];
1031 }
1032 }
1033 }
1034 I->EKin = ekin;
1035 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
1036}
1037
1038/** Scales ion velocities to match temperature.
1039 * In order to match Ions::ActualTemp with desired RunStruct::TargetTemp
1040 * each velocity in each dimension (for all types, all ions) is scaled
1041 * with the ratio of these two. Ions::EKin is recalculated.
1042 * \param *P Problem at hand
1043 */
1044void ScaleTemp(struct Problem *P)
1045{
1046 struct Ions *I = &P->Ion;
1047 int is, ia, d;
1048 double *U;
1049 double IonMass, a, ekin = 0;
1050 if (I->ActualTemp < MYEPSILON) fprintf(stderr,"ScaleTemp: I->ActualTemp = %lg",I->ActualTemp);
1051 double ScaleTempFactor = sqrt(P->R.TargetTemp/I->ActualTemp);
1052 for (is=0; is < I->Max_Types; is++) {
1053 IonMass = I->I[is].IonMass;
1054 a = 0.5*IonMass;
1055 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1056 U = &I->I[is].U[NDIM*ia];
1057 //if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
1058 for (d=0; d<NDIM; d++) {
1059 U[d] *= ScaleTempFactor;
1060 ekin += a * U[d]*U[d];
1061 }
1062 }
1063 }
1064 I->EKin = ekin;
1065 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
1066}
1067
1068/** Calculates mean force vector as stop criteria in structure optimization.
1069 * Calculates a mean force vector per ion as the usual euclidian norm over total
1070 * number of ions and dimensions, being the sum of each ion force (for all type,
1071 * all ions, all dims) squared.
1072 * The mean force is archived in RunStruct::MeanForce and printed to screen.
1073 * \param *P Problem at hand
1074 */
1075void GetOuterStop(struct Problem *P)
1076{
1077 struct RunStruct *R = &P->R;
1078 struct Ions *I = &P->Ion;
1079 int is, ia, IonNo=0, i, d;
1080 double MeanForce = 0.0;
1081 for (is=0; is < I->Max_Types; is++)
1082 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1083 IonNo++;
1084 for (d=0; d <NDIM; d++)
1085 MeanForce += I->I[is].FIon[d+NDIM*ia]*I->I[is].FIon[d+NDIM*ia];
1086 }
1087 for (i=MAXOLD-1; i > 0; i--)
1088 R->MeanForce[i] = R->MeanForce[i-1];
1089 MeanForce = sqrt(MeanForce/(IonNo*NDIM));
1090 R->MeanForce[0] = MeanForce;
1091 if (P->Call.out[LeaderOut] && (P->Par.me == 0))
1092 fprintf(stderr,"MeanForce = %e\n", R->MeanForce[0]);
1093}
1094
1095
Note: See TracBrowser for help on using the repository browser.