source: pcp/src/ions.c@ c510a7

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

char lengths of 255 and MAXDUMMYSTRING replaced with define MAXSTRINGSIZE in molecuilder and pcp
  • Property mode set to 100644
File size: 49.9 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(MAXSTRINGSIZE*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(MAXSTRINGSIZE, "IonsInitRead: Name");
300 I->I[i].Symbol = MallocString(MAXSTRINGSIZE, "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\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].FMagnetic[0+j*NDIM], I->I[i].FMagnetic[1+j*NDIM], I->I[i].FMagnetic[2+j*NDIM],
636 I->I[i].FEwald[0+j*NDIM], I->I[i].FEwald[1+j*NDIM], I->I[i].FEwald[2+j*NDIM]);
637 fflush(fout);
638 fprintf(fout, "MeanForce:\t%e\n", R->MeanForce[0]);
639}
640
641/** Reads Forces from CallOptions::ForcesFile.
642 * goes through all Ions per IonType and each dimension of #NDIM and parses in one line:
643 * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL,
644 * non-local IonType::FIonNL, ewald force IonType::FEwald
645 * \param *P Problem at hand
646 */
647void ParseIonForce(struct Problem *P)
648{
649 //struct RunStruct *Run = &P->R;
650 struct Ions *I = &P->Ion;
651 FILE *finput = NULL;
652 char line[1024];
653 int i,j;
654 double i1, j1;
655 double R[NDIM];
656
657 if (P->Call.ForcesFile == NULL) return;
658 finput = fopen(P->Call.ForcesFile, "r");
659 //fprintf(stderr, "File pointer says %p\n", finput);
660 if (finput == NULL)
661 Error(InitReading, "ForcesFile does not exist.");
662 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
663 for (i=0; i < I->Max_Types; i++)
664 for (j=0; j < I->I[i].Max_IonsOfType; j++)
665 if (feof(finput)) {
666 Error(InitReading, "ForcesFile does not contain enough lines.");
667 } else {
668 //fscanf(finput, "%s\n", line);
669 //if (strlen(line) < 100) {
670 //fprintf(stderr, "i %d, j %d, length of line %ld, line %s.\n", i,j, strlen(line), line);
671 //Error(InitReading, "ForcesFile does not contain enough lines or line is faulty.");
672 //} else
673 //fprintf(stderr, "Parsed line: '%s'\n", line);
674 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,
675 &R[0], &R[1], &R[2],
676 &I->I[i].FIon[0+j*NDIM], &I->I[i].FIon[1+j*NDIM], &I->I[i].FIon[2+j*NDIM],
677 &I->I[i].FIonL[0+j*NDIM], &I->I[i].FIonL[1+j*NDIM], &I->I[i].FIonL[2+j*NDIM],
678 &I->I[i].FIonNL[0+j*NDIM], &I->I[i].FIonNL[1+j*NDIM], &I->I[i].FIonNL[2+j*NDIM],
679 &I->I[i].FMagnetic[0+j*NDIM], &I->I[i].FMagnetic[1+j*NDIM], &I->I[i].FMagnetic[2+j*NDIM],
680 &I->I[i].FEwald[0+j*NDIM], &I->I[i].FEwald[1+j*NDIM], &I->I[i].FEwald[2+j*NDIM]);
681 if ((i != (int)i1) || (j != (int)j1))
682 Error(InitReading, "Line does not match to desired ion.");
683 }
684 fclose(finput);
685 //fprintf(stderr, "MeanForce:\t%e\n", Run->MeanForce[0]);
686};
687
688/** Frees memory IonType.
689 * All pointers from IonType and the pointer on it are free'd
690 * \param *I Ions structure to be free'd
691 * \sa IonsInitRead where memory is allocated
692 */
693void RemoveIonsRead(struct Ions *I)
694{
695 int i, it;
696 for (i=0; i < I->Max_Types; i++) {
697 Free(I->I[i].Name, "RemoveIonsRead: I->I[i].Name");
698 Free(I->I[i].Symbol, "RemoveIonsRead: I->I[i].Symbol");
699 for (it=0;it<I->I[i].Max_IonsOfType;it++) {
700 Free(I->I[i].moment[it], "RemoveIonsRead: I->I[i].moment[it]");
701 Free(I->I[i].sigma[it], "RemoveIonsRead: I->I[i].sigma[it]");
702 Free(I->I[i].sigma_rezi[it], "RemoveIonsRead: I->I[i].sigma_rezi[it]");
703 Free(I->I[i].sigma_PAS[it], "RemoveIonsRead: I->I[i].sigma_PAS[it]");
704 Free(I->I[i].sigma_rezi_PAS[it], "RemoveIonsRead: I->I[i].sigma_rezi_PAS[it]");
705 while (RemoveConstraintItem(&I->I[i], it)); // empty the constrained item list
706 }
707 Free(I->I[i].sigma, "RemoveIonsRead: I->I[i].sigma");
708 Free(I->I[i].sigma_rezi, "RemoveIonsRead: I->I[i].sigma_rezi");
709 Free(I->I[i].sigma_PAS, "RemoveIonsRead: I->I[i].sigma_PAS");
710 Free(I->I[i].sigma_rezi_PAS, "RemoveIonsRead: I->I[i].sigma_rezi_PAS");
711 Free(I->I[i].R, "RemoveIonsRead: I->I[i].R");
712 Free(I->I[i].R_old, "RemoveIonsRead: I->I[i].R_old");
713 Free(I->I[i].R_old_old, "RemoveIonsRead: I->I[i].R_old_old");
714 Free(I->I[i].Constraints, "RemoveIonsRead: I->I[i].Constraints");
715 Free(I->I[i].FIon, "RemoveIonsRead: I->I[i].FIon");
716 Free(I->I[i].FIon_old, "RemoveIonsRead: I->I[i].FIon_old");
717 Free(I->I[i].SearchDir, "RemoveIonsRead: I->I[i].SearchDir");
718 Free(I->I[i].GammaA, "RemoveIonsRead: I->I[i].GammaA");
719 Free(I->I[i].FIonL, "RemoveIonsRead: I->I[i].FIonL");
720 Free(I->I[i].FIonNL, "RemoveIonsRead: I->I[i].FIonNL");
721 Free(I->I[i].FMagnetic, "RemoveIonsRead: I->I[i].FMagnetic");
722 Free(I->I[i].U, "RemoveIonsRead: I->I[i].U");
723 Free(I->I[i].SFactor, "RemoveIonsRead: I->I[i].SFactor");
724 Free(I->I[i].IMT, "RemoveIonsRead: I->I[i].IMT");
725 Free(I->I[i].FEwald, "RemoveIonsRead: I->I[i].FEwald");
726 Free(I->I[i].alpha, "RemoveIonsRead: I->I[i].alpha");
727 }
728 if (I->R_cut == 0.0)
729 Free(I->RLatticeVec, "RemoveIonsRead: I->RLatticeVec");
730 Free(I->FTemp, "RemoveIonsRead: I->FTemp");
731 Free(I->I, "RemoveIonsRead: I->I");
732}
733
734/** Shifts center of gravity of ion forces IonType::FIon.
735 * First sums up all forces of movable ions for net center of gravity force,
736 * then reduces each force by this temp over Ions#Max_TotalIons, so that all in all
737 * the net force is 0.
738 * \param *P Problem at hand
739 * \sa CorrectVelocity()
740 */
741void CorrectForces(struct Problem *P)
742{
743 struct Ions *I = &P->Ion;
744 double *FIon;
745 double FTemp[NDIM];
746 int is, ia, d, NumIons = 0;
747 for (d=0; d<NDIM; d++)
748 FTemp[d] = 0.;
749 for (is=0; is < I->Max_Types; is++) { // calculate force temperature for each type ...
750 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of this type ...
751 FIon = &I->I[is].FIon[NDIM*ia];
752 if (I->I[is].IMT[ia] == MoveIon) { // .. if it's influenced
753 NumIons++;
754 for (d=0; d<NDIM; d++)
755 FTemp[d] += FIon[d];
756 }
757 }
758 }
759 //if (P->Files.MeOutMes != 1) return;
760 // fprintf(stderr, "TotalForce before: %e\t %e\t %e\n", FTemp[0], FTemp[1], FTemp[2]);
761 for (is=0; is < I->Max_Types; is++) { // and then reduce each by this value
762 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
763 FIon = &I->I[is].FIon[NDIM*ia];
764 if (I->I[is].IMT[ia] == MoveIon) {
765 for (d=0; d<NDIM; d++)
766 FIon[d] -= FTemp[d]/NumIons;
767 }
768 }
769 }
770/* for (d=0; d<NDIM; d++)
771 FTemp[d] = 0.;
772 for (is=0; is < I->Max_Types; is++) { // calculate force temperature for each type ...
773 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of this type ...
774 FIon = &I->I[is].FIon[NDIM*ia];
775 if (I->I[is].IMT[ia] == MoveIon) { // .. if it's influenced
776 for (d=0; d<NDIM; d++)
777 FTemp[d] += FIon[d];
778 }
779 }
780 }
781 if (P->Files.MeOutMes != 1) return;
782 fprintf(stderr, "TotalForce after: %e\t %e\t %e\n", FTemp[0], FTemp[1], FTemp[2]);
783 */
784}
785
786
787/** Shifts center of gravity of ion velocities (rather momentums).
788 * First sums up ion speed U times their IonMass (summed up for each
789 * dimension) in temp, then reduces the velocity by temp per Ions::TotalMass (so here
790 * total number of Ions is included)
791 * \param *P Problem at hand
792 * \sa CorrectForces()
793 */
794void CorrectVelocity(struct Problem *P)
795{
796 struct Ions *I = &P->Ion;
797 double *U;
798 double UTemp[NDIM];
799 int is, ia, d;
800 for (d=0; d<NDIM; d++)
801 UTemp[d] = 0;
802 for (is=0; is < I->Max_Types; is++) { // all types ...
803 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // ... all ions per type ...
804 U = &I->I[is].U[NDIM*ia];
805 //if (I->I[is].IMT[ia] == MoveIon) { // ... even FixedIon moves, only not by other's forces
806 for (d=0; d<NDIM; d++)
807 UTemp[d] += I->I[is].IonMass*U[d];
808 //}
809 }
810 }
811 for (is=0; is < I->Max_Types; is++) {
812 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
813 U = &I->I[is].U[NDIM*ia];
814 //if (I->I[is].IMT[ia] == MoveIon) {
815 for (d=0; d<NDIM; d++)
816 U[d] -= UTemp[d]/I->TotalMass; // shift by mean velocity over mass and number of ions
817 //}
818 }
819 }
820}
821
822/** Moves ions according to calculated force.
823 * Goes through each type of IonType and each ion therein, takes mass and
824 * Newton to move each ion to new coordinates IonType::R, while remembering the last
825 * two coordinates and the last force of the ion. Coordinates R are being
826 * transformed to inverted base, shifted by +-1.0 and back-transformed to
827 * real base.
828 * \param *P Problem at hand
829 * \sa CalculateIonForce
830 * \sa UpdateIonsU
831 * \warning U is not changed here, only used to move the ion.
832 */
833void UpdateIonsR(struct Problem *P)
834{
835 struct Ions *I = &P->Ion;
836 int is, ia, d;
837 double *R, *R_old, *R_old_old, *FIon, *FIon_old, *U, *FIonL, *FIonNL, *FEwald;
838 double IonMass, a;
839 double sR[NDIM], sRold[NDIM], sRoldold[NDIM];
840 const int delta_t = P->R.delta_t;
841 for (is=0; is < I->Max_Types; is++) { // go through all types ...
842 IonMass = I->I[is].IonMass;
843 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsR: IonMass = %lg",IonMass);
844 a = delta_t*0.5/IonMass; // F/m = a and thus: s = 0.5 * F/m * t^2 + v * t =: t * (F * a + v)
845 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
846 FIon = &I->I[is].FIon[NDIM*ia];
847 FIonL = &I->I[is].FIonL[NDIM*ia];
848 FIonNL = &I->I[is].FIonNL[NDIM*ia];
849 FEwald = &I->I[is].FEwald[NDIM*ia];
850 FIon_old = &I->I[is].FIon_old[NDIM*ia];
851 U = &I->I[is].U[NDIM*ia];
852 R = &I->I[is].R[NDIM*ia];
853 R_old = &I->I[is].R_old[NDIM*ia];
854 R_old_old = &I->I[is].R_old_old[NDIM*ia];
855 //if (I->I[is].IMT[ia] == MoveIon) { // even FixedIon moves, only not by other's forces
856 for (d=0; d<NDIM; d++) {
857 R_old_old[d] = R_old[d]; // shift old values
858 R_old[d] = R[d];
859 R[d] += delta_t*(U[d]+a*FIon[d]); // F = m * a and s = 0.5 * a * t^2
860 FIon_old[d] = FIon[d]; // store old force
861// FIon[d] = 0.; // zero all as a sign that's been moved
862// FIonL[d] = 0.;
863// FIonNL[d] = 0.;
864// FEwald[d] = 0.;
865// don't reset forces here anymore, as we OutputVis() after this and before next CalculateForce() (rendering them null in the cute graphs)
866 }
867 if (I->I[is].Constraints[ia] != NULL) { // override current resulting R if constrained is given
868 fprintf(stderr, "(%i) Using constrained coordinates R of step %d on ion (%i,%i)\n", P->Par.me, I->I[is].Constraints[ia]->step, is, ia);
869 if ((I->I[is].Constraints[ia]->step != P->R.OuterStep) && (P->Par.me == 0))
870 fprintf(stderr, "(%i) WARNING: constraint step %d does not match Outerstep %d for ion (%d,%d)\n", P->Par.me, I->I[is].Constraints[ia]->step, P->R.OuterStep, is, ia);
871 for (d=0; d<NDIM; d++)
872 R[d] = I->I[is].Constraints[ia]->R[d];
873 I->I[is].IMT[ia] = I->I[is].Constraints[ia]->IMT;
874 //RemoveConstraintItem(&I->I[is],ia); // and remove this first item from the list: // is done now at UpdateIonsU()
875 }
876 RMat33Vec3(sR, P->Lat.InvBasis, R);
877 RMat33Vec3(sRold, P->Lat.InvBasis, R_old);
878 RMat33Vec3(sRoldold, P->Lat.InvBasis, R_old_old);
879 for (d=0; d <NDIM; d++) {
880 while (sR[d] < 0) {
881 sR[d] += 1.0;
882 sRold[d] += 1.0;
883 sRoldold[d] += 1.0;
884 }
885 while (sR[d] >= 1.0) {
886 sR[d] -= 1.0;
887 sRold[d] -= 1.0;
888 sRoldold[d] -= 1.0;
889 }
890 }
891 RMat33Vec3(R, P->Lat.RealBasis, sR);
892 RMat33Vec3(R_old, P->Lat.RealBasis, sRold);
893 RMat33Vec3(R_old_old, P->Lat.RealBasis, sRoldold);
894 //}
895 }
896 }
897}
898
899/** Resets the forces array.
900 * \param *P Problem at hand
901 */
902void ResetForces(struct Problem *P)
903{
904 struct Ions *I = &P->Ion;
905 double *FIon, *FIonL, *FIonNL, *FEwald;
906 int is, ia, d;
907 for (is=0; is < I->Max_Types; is++) { // go through all types ...
908 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
909 FIon = &I->I[is].FIon[NDIM*ia];
910 FIonL = &I->I[is].FIonL[NDIM*ia];
911 FIonNL = &I->I[is].FIonNL[NDIM*ia];
912 FEwald = &I->I[is].FEwald[NDIM*ia];
913 for (d=0; d<NDIM; d++) {
914 FIon[d] = 0.; // zero all as a sign that's been moved
915 FIonL[d] = 0.;
916 FIonNL[d] = 0.;
917 FEwald[d] = 0.;
918 }
919 }
920 }
921};
922
923/** Changes ion's velocity according to acting force.
924 * IonType::U is changed by the weighted force of actual step and last one
925 * according to Newton.
926 * \param *P Problem at hand
927 * \sa UpdateIonsR
928 * \sa CalculateIonForce
929 * \warning R is not changed here.
930 */
931void UpdateIonsU(struct Problem *P)
932{
933 struct Ions *I = &P->Ion;
934 int is, ia, d;
935 double *FIon, *FIon_old, *U;
936 double IonMass, a;
937 const int delta_t = P->R.delta_t;
938 for (is=0; is < I->Max_Types; is++) {
939 IonMass = I->I[is].IonMass;
940 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsU: IonMass = %lg",IonMass);
941 a = delta_t*0.5/IonMass; // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
942 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
943 U = &I->I[is].U[NDIM*ia];
944 if (I->I[is].Constraints[ia] != NULL) { // override current resulting R and U if constrained is given
945 fprintf(stderr, "(%i) Setting constrained motion U at step %d on ion (%i,%i)\n", P->Par.me, I->I[is].Constraints[ia]->step, is, ia);
946 if ((I->I[is].Constraints[ia]->step != P->R.OuterStep) && (P->Par.me == 0))
947 fprintf(stderr, "(%i) WARNING: constraint step %d does not match Outerstep %d for ion (%d,%d)\n", P->Par.me, I->I[is].Constraints[ia]->step, P->R.OuterStep, is, ia);
948 for (d=0; d<NDIM; d++)
949 U[d] = I->I[is].Constraints[ia]->U[d];
950 RemoveConstraintItem(&I->I[is],ia); // and remove this first item from the list
951 } else {
952 FIon = &I->I[is].FIon[NDIM*ia];
953 FIon_old = &I->I[is].FIon_old[NDIM*ia];
954 //if (I->I[is].IMT[ia] == MoveIon)
955 for (d=0; d<NDIM; d++) {
956 U[d] += a * (FIon[d]+FIon_old[d]);
957 }
958 }
959 }
960 }
961}
962
963/** CG process to optimise structure.
964 * \param *P Problem at hand
965 */
966void UpdateIons(struct Problem *P)
967{
968 struct Ions *I = &P->Ion;
969 int is, ia, d;
970 double *R, *R_old, *R_old_old, *FIon, *S, *GammaA;
971 double IonFac, GammaB; //, GammaT;
972 double FNorm, StepNorm;
973 for (is=0; is < I->Max_Types; is++) {
974 IonFac = I->I[is].IonFac;
975 GammaA = I->I[is].GammaA;
976 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
977 FIon = &I->I[is].FIon[NDIM*ia];
978 S = &I->I[is].SearchDir[NDIM*ia];
979 GammaB = GammaA[ia];
980 GammaA[ia] = RSP3(FIon,S);
981 FNorm = sqrt(RSP3(FIon,FIon));
982 StepNorm = log(1.+IonFac*FNorm); // Fix Hack
983 if (FNorm != 0)
984 StepNorm /= sqrt(RSP3(FIon,FIon));
985 else
986 StepNorm = 0;
987 //if (GammaB != 0.0) {
988 // GammaT = GammaA[ia]/GammaB;
989 // for (d=0; d>NDIM; d++)
990 // S[d] = FIon[d]+GammaT*S[d];
991 // } else {
992 for (d=0; d<NDIM; d++)
993 S[d] = StepNorm*FIon[d];
994 // }
995 R = &I->I[is].R[NDIM*ia];
996 R_old = &I->I[is].R_old[NDIM*ia];
997 R_old_old = &I->I[is].R_old_old[NDIM*ia];
998 if (I->I[is].IMT[ia] == MoveIon)
999 for (d=0; d<NDIM; d++) {
1000 R_old_old[d] = R_old[d];
1001 R_old[d] = R[d];
1002 R[d] += S[d]; // FixHack*IonFac;
1003 }
1004 }
1005 }
1006}
1007
1008/** Print coordinates of all ions to stdout.
1009 * \param *P Problem at hand
1010 * \param actual (1 - don't at current RunStruct#OuterStep, 0 - do)
1011 */
1012void OutputIonCoordinates(struct Problem *P, int actual)
1013{
1014 //struct RunStruct *R = &P->R;
1015 struct Ions *I = &P->Ion;
1016 char filename[MAXSTRINGSIZE];
1017 FILE *output;
1018 int is, ia, nr = 0;
1019 double Bohr = (I->IsAngstroem) ? 1./ANGSTROEMINBOHRRADIUS : 1.;
1020/* if (P->Par.me == 0 && P->Call.out[ReadOut]) {
1021 fprintf(stderr, "(%i) ======= Differential Ion Coordinates =========\n",P->Par.me);
1022 for (is=0; is < I->Max_Types; is++)
1023 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++)
1024 //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]);
1025 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]);
1026 fprintf(stderr, "(%i) =========================================\n",P->Par.me);
1027 }*/
1028 if (actual)
1029 sprintf(filename, "%s.MD", P->Call.MainParameterFile);
1030 else
1031 sprintf(filename, "%s.opt", P->Call.MainParameterFile);
1032 if (((actual) && (P->R.OuterStep <= 0)) || ((!actual) && (P->R.StructOptStep == 1))) { // on first step write complete config file
1033 WriteParameters(P, filename);
1034 } else { // afterwards simply add the current ion coordinates as constrained steps
1035 if ((P->Par.me == 0) && (output = fopen(filename,"a"))) {
1036 fprintf(output,"# ====== MD step %d ========= \n", (actual) ? P->R.OuterStep : P->R.StructOptStep);
1037 for (is=0; is < I->Max_Types; is++)
1038 for (ia=0;ia<I->I[is].Max_IonsOfType;ia++) {
1039 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);
1040 }
1041 fflush(output);
1042 }
1043 }
1044}
1045
1046/** Calculates kinetic energy (Ions::EKin) of movable Ions.
1047 * Does 0.5 / IonType::IonMass * IonTye::U^2 for each Ion,
1048 * also retrieves actual temperatur by 2/3 from just
1049 * calculated Ions::EKin.
1050 * \param *P Problem at hand
1051 */
1052void CalculateEnergyIonsU(struct Problem *P)
1053{
1054 struct Ions *I = &P->Ion;
1055 int is, ia, d;
1056 double *U;
1057 double IonMass, a, ekin = 0;
1058 for (is=0; is < I->Max_Types; is++) {
1059 IonMass = I->I[is].IonMass;
1060 a = 0.5*IonMass;
1061 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1062 U = &I->I[is].U[NDIM*ia];
1063 //if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
1064 for (d=0; d<NDIM; d++) {
1065 ekin += a * U[d]*U[d];
1066 }
1067 }
1068 }
1069 I->EKin = ekin;
1070 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
1071}
1072
1073/** Scales ion velocities to match temperature.
1074 * In order to match Ions::ActualTemp with desired RunStruct::TargetTemp
1075 * each velocity in each dimension (for all types, all ions) is scaled
1076 * with the ratio of these two. Ions::EKin is recalculated.
1077 * \param *P Problem at hand
1078 */
1079void ScaleTemp(struct Problem *P)
1080{
1081 struct Ions *I = &P->Ion;
1082 int is, ia, d;
1083 double *U;
1084 double IonMass, a, ekin = 0;
1085 if (I->ActualTemp < MYEPSILON) fprintf(stderr,"ScaleTemp: I->ActualTemp = %lg",I->ActualTemp);
1086 double ScaleTempFactor = sqrt(P->R.TargetTemp/I->ActualTemp);
1087 for (is=0; is < I->Max_Types; is++) {
1088 IonMass = I->I[is].IonMass;
1089 a = 0.5*IonMass;
1090 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1091 U = &I->I[is].U[NDIM*ia];
1092 //if (I->I[is].IMT[ia] == MoveIon) // even FixedIon moves, only not by other's forces
1093 for (d=0; d<NDIM; d++) {
1094 U[d] *= ScaleTempFactor;
1095 ekin += a * U[d]*U[d];
1096 }
1097 }
1098 }
1099 I->EKin = ekin;
1100 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
1101}
1102
1103/** Calculates mean force vector as stop criteria in structure optimization.
1104 * Calculates a mean force vector per ion as the usual euclidian norm over total
1105 * number of ions and dimensions, being the sum of each ion force (for all type,
1106 * all ions, all dims) squared.
1107 * The mean force is archived in RunStruct::MeanForce and printed to screen.
1108 * \param *P Problem at hand
1109 */
1110void GetOuterStop(struct Problem *P)
1111{
1112 struct RunStruct *R = &P->R;
1113 struct Ions *I = &P->Ion;
1114 int is, ia, IonNo=0, i, d;
1115 double MeanForce = 0.0;
1116 for (is=0; is < I->Max_Types; is++)
1117 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1118 IonNo++;
1119 for (d=0; d <NDIM; d++)
1120 MeanForce += I->I[is].FIon[d+NDIM*ia]*I->I[is].FIon[d+NDIM*ia];
1121 }
1122 for (i=MAXOLD-1; i > 0; i--)
1123 R->MeanForce[i] = R->MeanForce[i-1];
1124 MeanForce = sqrt(MeanForce/(IonNo*NDIM));
1125 R->MeanForce[0] = MeanForce;
1126 if (P->Call.out[LeaderOut] && (P->Par.me == 0))
1127 fprintf(stderr,"MeanForce = %e\n", R->MeanForce[0]);
1128}
1129
1130
Note: See TracBrowser for help on using the repository browser.