source: pcp/src/ions.c@ e00f47

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

Free(): now takes a debug string to know where free error occured
AddConstraintItem(), RemoveConstraintItem() and IonConstrainted structure to take constrained motion info
  • Property mode set to 100644
File size: 30.8 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 "mymath.h"
32
33
34// use double precision fft when we have it
35#ifdef HAVE_CONFIG_H
36#include <config.h>
37#endif
38
39#ifdef HAVE_DFFTW_H
40#include "dfftw.h"
41#else
42#include "fftw.h"
43#endif
44
45#include "data.h"
46#include "errors.h"
47#include "helpers.h"
48#include "mymath.h"
49#include "ions.h"
50#include "init.h"
51
52/** Readin of the ion section of parameter file.
53 * Among others the following paramaters are read from file:
54 * Ions::MaxTypes, IonType::Max_IonsOfType, IonType::Z, IonType::R,
55 * IonType:IMT, ...
56 * \param P Problem at hand (containing Ions and Lattice structures)
57 * \param *source file pointer
58 */
59void IonsInitRead(struct Problem *P, FILE *source) {
60 struct Ions *I = &P->Ion;
61 //struct RunStruct *Run = &P->R;
62 int i,it,j,IMT,BorAng, d,me, count;
63 struct Lattice *L = &P->Lat;
64 int k;
65 int relative; // whether read-in coordinate are relative (1) or absolute
66 double Bohr = ANGSTROEMINBOHRRADIUS; /* Angstroem */
67 double sigma, R[3], U[3];
68
69 I->Max_TotalIons = 0;
70 I->TotalZval = 0;
71 MPI_Comm_rank(MPI_COMM_WORLD, &me);
72 ParseForParameter(P->Call.out[ReadOut],source,"RCut", 0, 1, 1, double_type, &I->R_cut, critical);
73 ParseForParameter(P->Call.out[ReadOut],source,"IsAngstroem", 0, 1, 1, int_type, &BorAng, critical);
74 if (!I->IsAngstroem) {
75 Bohr = 1.0;
76 } else { // adapt RealBasis (is in Angstroem as well)
77 SM(P->Lat.RealBasis, Bohr, NDIM_NDIM);
78 RMatReci3(P->Lat.InvBasis, P->Lat.RealBasis);
79 }
80 ParseForParameter(P->Call.out[ReadOut],source,"MaxTypes", 0, 1, 1, int_type, &I->Max_Types, critical);
81 I->I = (struct IonType *) Malloc(sizeof(struct IonType)*I->Max_Types, "IonsInitRead: IonType");
82 if (!ParseForParameter(P->Call.out[ReadOut],source,"RelativeCoord", 0, 1, 1, int_type, &relative, optional))
83 relative = 0;
84 if (!ParseForParameter(P->Call.out[ReadOut],source,"StructOpt", 0, 1, 1, int_type, &I->StructOpt, optional))
85 I->StructOpt = 0;
86
87 /* Ions Data */
88 I->RLatticeVec = NULL;
89 I->TotalMass = 0;
90 char *free_name, *name;
91 name = free_name = Malloc(255*sizeof(char),"IonsInitRead: Name");
92 for (i=0; i < I->Max_Types; i++) {
93 sprintf(name,"Ion_Type%i",i+1);
94 I->I[i].corecorr = NotCoreCorrected;
95 I->I[i].Name = MallocString(255, "IonsInitRead: Name");
96 I->I[i].Symbol = MallocString(255, "IonsInitRead: Symbol");
97 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 1, 1, int_type, &I->I[i].Max_IonsOfType, critical);
98 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, int_type, &I->I[i].Z, critical);
99 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &I->I[i].rgauss, critical);
100 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &I->I[i].l_max, critical);
101 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 5, 1, int_type, &I->I[i].l_loc, critical);
102 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 6, 1, double_type, &I->I[i].IonMass, critical);
103 if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 7, 1, string_type, &I->I[i].Name[0], optional))
104 sprintf(&I->I[i].Name[0],"type%i",i);
105 if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 8, 1, string_type, &I->I[i].Symbol[0], optional))
106 sprintf(&I->I[i].Symbol[0],"t%i",i);
107 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);
108 I->I[i].sigma = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma");
109 I->I[i].sigma_rezi = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_rezi");
110 I->I[i].sigma_PAS = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_PAS");
111 I->I[i].sigma_rezi_PAS = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_rezi_PAS");
112 for (it=0;it<I->I[i].Max_IonsOfType;it++) {
113 I->I[i].sigma[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: sigma");
114 I->I[i].sigma_rezi[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: sigma_rezi");
115 I->I[i].sigma_PAS[it] = (double *) Malloc(sizeof(double) * NDIM, "IonsInitRead: sigma_PAS");
116 I->I[i].sigma_rezi_PAS[it] = (double *) Malloc(sizeof(double) * NDIM, "IonsInitRead: sigma_rezi_PAS");
117 }
118 //I->I[i].IonFac = I->I[i].IonMass;
119 I->I[i].IonFac = 1/I->I[i].IonMass;
120 I->TotalMass += I->I[i].Max_IonsOfType*I->I[i].IonMass;
121 I->I[i].ZFactor = -I->I[i].Z*4.*PI;
122 I->I[i].alpha = (double *) Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: alpha");
123 I->I[i].R = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: R");
124 I->I[i].R_old = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: R_old");
125 I->I[i].R_old_old = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: R_old_old");
126 I->I[i].FIon = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIon");
127 I->I[i].FIon_old = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIon_old");
128 SetArrayToDouble0(I->I[i].FIon_old, NDIM*I->I[i].Max_IonsOfType);
129 I->I[i].SearchDir = Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: SearchDir");
130 I->I[i].GammaA = Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: GammaA");
131 SetArrayToDouble0(I->I[i].GammaA, I->I[i].Max_IonsOfType);
132 I->I[i].U = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: U");
133 SetArrayToDouble0(I->I[i].U, NDIM*I->I[i].Max_IonsOfType);
134 I->I[i].SFactor = NULL;
135 I->I[i].IMT = (enum IonMoveType *) Malloc(sizeof(enum IonMoveType) * I->I[i].Max_IonsOfType, "IonsInitRead: IMT");
136 I->I[i].FIonL = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIonL");
137 I->I[i].FIonNL = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIonNL");
138 I->I[i].FEwald = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FEwald");
139 I->I[i].alpha = (double *) Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: alpha");
140 for (j=0; j < I->I[i].Max_IonsOfType; j++) {
141 I->I[i].alpha[j] = 2.;
142 sprintf(name,"Ion_Type%i_%i",i+1,j+1);
143 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 1, 1, double_type, &R[0], critical);
144 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, double_type, &R[1], critical);
145 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &R[2], critical);
146 ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &IMT, critical);
147 // change if coordinates were relative
148 if (relative) {
149 //fprintf(stderr,"(%i)Ion coordinates are relative %i ... \n",P->Par.me, relative);
150 RMat33Vec3(&I->I[i].R[0+NDIM*j], L->RealBasis, R); // multiply with real basis
151 } else {
152 for (k=0;k<NDIM;k++) // and copy to old vector
153 I->I[i].R[k+NDIM*j] = R[k];
154 }
155 if (IMT < 0 || IMT >= MaxIonMoveType) {
156 fprintf(stderr,"Bad Ion MoveType set to MoveIon for Ion (%i,%i)\n",i,j);
157 IMT = 0;
158 }
159
160 if ((P->Call.out[ReadOut])) { // rotate coordinates
161 fprintf(stderr,"(%i) coordinates of Ion %i: (x,y,z) = (",P->Par.me,j);
162 for(k=0;k<NDIM;k++) {
163 fprintf(stderr,"%lg ",I->I[i].R[k+NDIM*j]);
164 if (k != NDIM-1) fprintf(stderr,", ");
165 else fprintf(stderr,")\n");
166 }
167 }
168
169 I->I[i].IMT[j] = (enum IonMoveType)IMT;
170 SM(&I->I[i].R[NDIM*j], 1./Bohr, NDIM);
171 RMat33Vec3(R, L->InvBasis, &I->I[i].R[NDIM*j]);
172 for (d=0; d <NDIM; d++) {
173 while (R[d] < 0)
174 R[d] += 1.0;
175 while (R[d] >= 1.0)
176 R[d] -= 1.0;
177 }
178 RMat33Vec3(&I->I[i].R[NDIM*j], L->RealBasis, R);
179 for (d=0; d <NDIM; d++) {
180 I->I[i].R_old_old[d+NDIM*j] = I->I[i].R_old[d+NDIM*j] = I->I[i].R[d+NDIM*j];
181 }
182 I->Max_TotalIons++;
183 }
184 }
185 Free(free_name, "IonsInitRead: free_name");
186 I->Max_Max_IonsOfType = 0;
187 for (i=0; i < I->Max_Types; i++)
188 I->Max_Max_IonsOfType = MAX(I->Max_Max_IonsOfType, I->I[i].Max_IonsOfType);
189 I->FTemp = (double *) Malloc(sizeof(double) * NDIM * I->Max_Max_IonsOfType, "IonsInitRead:");
190}
191
192/** Allocates memory for a IonConstrained item and adds to end of list.
193 * \param *I IonType structure
194 * \param ion which ion of the type
195 * \return pointer to created item
196 */
197struct IonConstrained * AddConstraintItem(struct IonType *I, int ion)
198{
199 struct IonConstrained **ptr = &(I->Constraints[ion]);
200 while (*ptr != NULL) { // advance to end of list
201 ptr = &((*ptr)->next);
202 }
203 // allocated memory for structure and items within
204 *ptr = (struct IonConstrained *) Malloc(sizeof(struct IonConstrained), "CreateConstraintItem: IonConstrained");
205 (*ptr)->R = (double *) Malloc(sizeof(double)*NDIM, "AddConstraintItem: *R");
206 (*ptr)->U = (double *) Malloc(sizeof(double)*NDIM, "AddConstraintItem: *U");
207 (*ptr)->next = NULL;
208 return *ptr;
209}
210
211/** Removes first of item of IonConstrained list for given \a ion.
212 * Frees memory of items within and then
213 * \param *I IonType structure
214 * \param ion which ion of the type
215 * \return 1 - success, 0 - failure (list is already empty, start pointer was NULL)
216 */
217int RemoveConstraintItem(struct IonType *I, int ion)
218{
219 struct IonConstrained **ptr = &(I->Constraints[ion]);
220 struct IonConstrained *next_ptr;
221 if (*ptr != NULL) {
222 next_ptr = (*ptr)->next;
223 Free((*ptr)->R, "RemoveConstraintItem: R");
224 Free((*ptr)->U, "RemoveConstraintItem: U");
225 Free((*ptr), "RemoveConstraintItem: IonConstrained structure");
226 *ptr = next_ptr; // set start of list to next item (automatically is null if ptr was last item)
227 return 1;
228 } else {
229 return 0;
230 }
231}
232
233/** Calculated Ewald force.
234 * \f[
235 * 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|}
236 * \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 )
237 * \qquad (4.10)
238 * \f]
239 * Calculates the core-to-core force between the ions of all super cells.
240 * In order for this series to converge there must be a certain summation
241 * applied, which is the ewald summation and which is nothing else but to move
242 * in a circular motion from the current cell to the outside up to Ions::R_cut.
243 * \param *P Problem at hand
244 * \param first additional calculation beforehand if != 0
245 */
246void CalculateEwald(struct Problem *P, int first)
247{
248 int r,is,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec;
249 int MaxPar=P->Par.procs, ParMe=P->Par.me, MaxLocalVec, StartVec;
250 double rcut2,r2,erre2,addesr,addpre,arg,fac,gkl,esr,fxx,repand;
251 double R1[3];
252 double R2[3];
253 double R3[3];
254 double t[3];
255 struct Lattice *L = &P->Lat;
256 struct PseudoPot *PP = &P->PP;
257 struct Ions *I = &P->Ion;
258 if (I->R_cut != 0.0) {
259 if (first) {
260 SpeedMeasure(P, InitSimTime, StopTimeDo);
261 rcut2 = I->R_cut*I->R_cut;
262 ActualVec =0;
263 MaxCell = (int)5.*I->R_cut/pow(L->Volume,1./3.);
264 if (MaxCell < 2) MaxCell = 2;
265 for (i = -MaxCell; i <= MaxCell; i++) {
266 for (j = -MaxCell; j <= MaxCell; j++) {
267 for (k = -MaxCell; k <= MaxCell; k++) {
268 r2 = 0.0;
269 for (ir=0; ir <3; ir++) {
270 t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
271 r2 = t[ir]*t[ir];
272 }
273 Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
274 if ((r2 <= rcut2) || Is_Neighb_Cell) {
275 ActualVec++;
276 }
277 }
278 }
279 }
280 MaxVec = ActualVec;
281 I->MaxVec = ActualVec;
282 MaxLocalVec = MaxVec / MaxPar;
283 StartVec = ParMe*MaxLocalVec;
284 r = MaxVec % MaxPar;
285 if (ParMe >= r) {
286 StartVec += r;
287 } else {
288 StartVec += ParMe;
289 MaxLocalVec++;
290 }
291 I->MaxLocalVec = MaxLocalVec;
292 LocalActualVec = ActualVec = 0;
293 I->RLatticeVec = (double *)
294 Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald:");
295 for (i = -MaxCell; i <= MaxCell; i++) {
296 for (j = -MaxCell; j <= MaxCell; j++) {
297 for (k = -MaxCell; k <= MaxCell; k++) {
298 r2 = 0.0;
299 for (ir=0; ir <3; ir++) {
300 t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
301 r2 = t[ir]*t[ir];
302 }
303 Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
304 if ((r2 <= rcut2) || Is_Neighb_Cell) {
305 if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) {
306 I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
307 I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
308 I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
309 LocalActualVec++;
310 }
311 ActualVec++;
312 }
313 }
314 }
315 }
316 SpeedMeasure(P, InitSimTime, StartTimeDo);
317 }
318 SpeedMeasure(P, EwaldTime, StartTimeDo);
319 esr = 0.0;
320 for (is1=0;is1 < I->Max_Types; is1++)
321 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
322 for (i=0; i< NDIM; i++)
323 I->FTemp[i+NDIM*(ia1)] = 0;
324 for (is1=0;is1 < I->Max_Types; is1++) {
325 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
326 for (i=0;i<3;i++) {
327 R1[i]=I->I[is1].R[i+NDIM*ia1];
328 }
329 for (ir2=0;ir2<I->MaxLocalVec;ir2++) {
330 for (is2=0;is2 < I->Max_Types; is2++) {
331 gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
332 for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
333 for (i=0;i<3;i++) {
334 R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2];
335 }
336 erre2=0.0;
337 for (i=0;i<3;i++) {
338 R3[i] = R1[i]-R2[i];
339 erre2+=R3[i]*R3[i];
340 }
341 if (erre2 > MYEPSILON) {
342 arg=sqrt(erre2);
343 fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;
344
345 arg *= gkl;
346 addesr = derf(arg);
347 addesr = (1.0-addesr)*fac;
348 esr += addesr;
349 addpre=exp(-arg*arg)*gkl;
350 addpre=PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;
351 repand=(addesr+addpre)/erre2;
352 for (i=0;i<3;i++) {
353 fxx=repand*R3[i];
354 /*fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);*/
355 I->FTemp[i+NDIM*(ia1)] += fxx;
356 I->FTemp[i+NDIM*(ia2)] -= fxx;
357 }
358 }
359 }
360 }
361 }
362 }
363 }
364 } else {
365 esr = 0.0;
366 for (is1=0;is1 < I->Max_Types; is1++)
367 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
368 for (i=0; i< NDIM; i++)
369 I->FTemp[i+NDIM*(ia1)] = 0.;
370 }
371 SpeedMeasure(P, EwaldTime, StopTimeDo);
372 for (is=0;is < I->Max_Types; is++) {
373 MPI_Allreduce (I->FTemp , I->I[is].FEwald, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);
374 }
375 MPI_Allreduce ( &esr, &L->E->AllTotalIonsEnergy[EwaldEnergy], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm);
376}
377
378/** Sum up all forces acting on Ions.
379 * IonType::FIon = IonType::FEwald + IonType::FIonL + IonType::FIonNL for all
380 * dimensions #NDIM and each Ion of IonType
381 * \param *P Problem at hand
382 */
383void CalculateIonForce(struct Problem *P)
384{
385 struct Ions *I = &P->Ion;
386 int i,j,d;
387 for (i=0; i < I->Max_Types; i++)
388 for (j=0; j < I->I[i].Max_IonsOfType; j++)
389 for (d=0; d<NDIM;d++)
390 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];
391}
392
393/** Write Forces to FileData::ForcesFile.
394 * goes through all Ions per IonType and each dimension of #NDIM and prints in one line:
395 * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL,
396 * non-local IonType::FIonNL, ewald force IonType::FEwald
397 * \param *P Problem at hand
398 */
399void OutputIonForce(struct Problem *P)
400{
401 struct RunStruct *R = &P->R;
402 struct FileData *F = &P->Files;
403 struct Ions *I = &P->Ion;
404 FILE *fout = NULL;
405 int i,j;
406 if (F->MeOutMes != 1) return;
407 fout = F->ForcesFile;
408 for (i=0; i < I->Max_Types; i++)
409 for (j=0; j < I->I[i].Max_IonsOfType; j++)
410 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,
411 I->I[i].R[0+j*NDIM], I->I[i].R[1+j*NDIM], I->I[i].R[2+j*NDIM],
412 I->I[i].FIon[0+j*NDIM], I->I[i].FIon[1+j*NDIM], I->I[i].FIon[2+j*NDIM],
413 I->I[i].FIonL[0+j*NDIM], I->I[i].FIonL[1+j*NDIM], I->I[i].FIonL[2+j*NDIM],
414 I->I[i].FIonNL[0+j*NDIM], I->I[i].FIonNL[1+j*NDIM], I->I[i].FIonNL[2+j*NDIM],
415 I->I[i].FEwald[0+j*NDIM], I->I[i].FEwald[1+j*NDIM], I->I[i].FEwald[2+j*NDIM]);
416 fflush(fout);
417 fprintf(fout, "MeanForce:\t%e\n", R->MeanForce[0]);
418}
419
420/** Frees memory IonType.
421 * All pointers from IonType and the pointer on it are free'd
422 * \param *I Ions structure to be free'd
423 * \sa IonsInitRead where memory is allocated
424 */
425void RemoveIonsRead(struct Ions *I)
426{
427 int i, it;
428 for (i=0; i < I->Max_Types; i++) {
429 Free(I->I[i].Name, "RemoveIonsRead: I->I[i].Name");
430 Free(I->I[i].Symbol, "RemoveIonsRead: I->I[i].Symbol");
431 for (it=0;it<I->I[i].Max_IonsOfType;it++) {
432 Free(I->I[i].moment[it], "RemoveIonsRead: I->I[i].moment[it]");
433 Free(I->I[i].sigma[it], "RemoveIonsRead: I->I[i].sigma[it]");
434 Free(I->I[i].sigma_rezi[it], "RemoveIonsRead: I->I[i].sigma_rezi[it]");
435 Free(I->I[i].sigma_PAS[it], "RemoveIonsRead: I->I[i].sigma_PAS[it]");
436 Free(I->I[i].sigma_rezi_PAS[it], "RemoveIonsRead: I->I[i].sigma_rezi_PAS[it]");
437 while (RemoveConstraintItem(&I->I[i], it)); // empty the constrained item list
438 }
439 Free(I->I[i].sigma, "RemoveIonsRead: I->I[i].sigma");
440 Free(I->I[i].sigma_rezi, "RemoveIonsRead: I->I[i].sigma_rezi");
441 Free(I->I[i].sigma_PAS, "RemoveIonsRead: I->I[i].sigma_PAS");
442 Free(I->I[i].sigma_rezi_PAS, "RemoveIonsRead: I->I[i].sigma_rezi_PAS");
443 Free(I->I[i].R, "RemoveIonsRead: I->I[i].R");
444 Free(I->I[i].R_old, "RemoveIonsRead: I->I[i].R_old");
445 Free(I->I[i].R_old_old, "RemoveIonsRead: I->I[i].R_old_old");
446 Free(I->I[i].Constraints, "RemoveIonsRead: I->I[i].Constraints");
447 Free(I->I[i].FIon, "RemoveIonsRead: I->I[i].FIon");
448 Free(I->I[i].FIon_old, "RemoveIonsRead: I->I[i].FIon_old");
449 Free(I->I[i].SearchDir, "RemoveIonsRead: I->I[i].SearchDir");
450 Free(I->I[i].GammaA, "RemoveIonsRead: I->I[i].GammaA");
451 Free(I->I[i].FIonL, "RemoveIonsRead: I->I[i].FIonL");
452 Free(I->I[i].FIonNL, "RemoveIonsRead: I->I[i].FIonNL");
453 Free(I->I[i].FMagnetic, "RemoveIonsRead: I->I[i].FMagnetic");
454 Free(I->I[i].U, "RemoveIonsRead: I->I[i].U");
455 Free(I->I[i].SFactor, "RemoveIonsRead: I->I[i].SFactor");
456 Free(I->I[i].IMT, "RemoveIonsRead: I->I[i].IMT");
457 Free(I->I[i].FEwald, "RemoveIonsRead: I->I[i].FEwald");
458 Free(I->I[i].alpha, "RemoveIonsRead: I->I[i].alpha");
459 }
460 if (I->R_cut == 0.0)
461 Free(I->RLatticeVec, "RemoveIonsRead: I->RLatticeVec");
462 Free(I->FTemp, "RemoveIonsRead: I->FTemp");
463 Free(I->I, "RemoveIonsRead: I->I");
464}
465
466/** Shifts center of gravity of ion forces IonType::FIon.
467 * First sums up all forces of movable ions to a "force temperature",
468 * then reduces each force by this temp, so that all in all
469 * the net force is 0.
470 * \param *P Problem at hand
471 * \sa CorrectVelocity
472 * \note why is FTemp not divided by number of ions: is probably correct, but this
473 * function was not really meaningful from the beginning.
474 */
475void CorrectForces(struct Problem *P)
476{
477 struct Ions *I = &P->Ion;
478 double *FIon;
479 double FTemp[NDIM];
480 int is, ia, d;
481 for (d=0; d<NDIM; d++)
482 FTemp[d] = 0;
483 for (is=0; is < I->Max_Types; is++) { // calculate force temperature for each type ...
484 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of this type ...
485 FIon = &I->I[is].FIon[NDIM*ia];
486 if (I->I[is].IMT[ia] == MoveIon) { // .. if it's movable
487 for (d=0; d<NDIM; d++)
488 FTemp[d] += FIon[d];
489 }
490 }
491 }
492 for (is=0; is < I->Max_Types; is++) { // and then reduce each by this value
493 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
494 FIon = &I->I[is].FIon[NDIM*ia];
495 if (I->I[is].IMT[ia] == MoveIon) {
496 for (d=0; d<NDIM; d++)
497 FIon[d] -= FTemp[d]; // (?) why not -= FTemp[d]/I->Max_TotalIons
498 }
499 }
500 }
501}
502
503
504/** Shifts center of gravity of ion velocities (rather momentums).
505 * First sums up ion speed U times their IonMass (summed up for each
506 * dimension), then reduces the velocity by temp per Ions::TotalMass (so here
507 * total number of Ions is included)
508 * \param *P Problem at hand
509 * \sa CorrectForces
510 */
511void CorrectVelocity(struct Problem *P)
512{
513 struct Ions *I = &P->Ion;
514 double *U;
515 double UTemp[NDIM];
516 int is, ia, d;
517 for (d=0; d<NDIM; d++)
518 UTemp[d] = 0;
519 for (is=0; is < I->Max_Types; is++) { // all types ...
520 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // ... all ions per type ...
521 U = &I->I[is].U[NDIM*ia];
522 if (I->I[is].IMT[ia] == MoveIon) { // ... if it's movable
523 for (d=0; d<NDIM; d++)
524 UTemp[d] += I->I[is].IonMass*U[d];
525 }
526 }
527 }
528 for (is=0; is < I->Max_Types; is++) {
529 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
530 U = &I->I[is].U[NDIM*ia];
531 if (I->I[is].IMT[ia] == MoveIon) {
532 for (d=0; d<NDIM; d++)
533 U[d] -= UTemp[d]/I->TotalMass; // shift by mean velocity over mass and number of ions
534 }
535 }
536 }
537}
538
539/** Moves ions according to calculated force.
540 * Goes through each type of IonType and each ion therein, takes mass and
541 * Newton to move each ion to new coordinates IonType::R, while remembering the last
542 * two coordinates and the last force of the ion. Coordinates R are being
543 * transformed to inverted base, shifted by +-1.0 and back-transformed to
544 * real base.
545 * \param *P Problem at hand
546 * \sa CalculateIonForce
547 * \sa UpdateIonsU
548 * \warning U is not changed here, only used to move the ion.
549 */
550void UpdateIonsR(struct Problem *P)
551{
552 struct Ions *I = &P->Ion;
553 int is, ia, d;
554 double *R, *R_old, *R_old_old, *FIon, *FIon_old, *U, *FIonL, *FIonNL, *FEwald;
555 double IonMass, a;
556 double sR[NDIM], sRold[NDIM], sRoldold[NDIM];
557 const int delta_t = P->R.delta_t;
558 for (is=0; is < I->Max_Types; is++) { // go through all types ...
559 IonMass = I->I[is].IonMass;
560 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsR: IonMass = %lg",IonMass);
561 a = delta_t*0.5/IonMass; // F/m = a and thus: s = 0.5 * F/m * t^2 + v * t =: t * (F * a + v)
562 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
563 FIon = &I->I[is].FIon[NDIM*ia];
564 FIonL = &I->I[is].FIonL[NDIM*ia];
565 FIonNL = &I->I[is].FIonNL[NDIM*ia];
566 FEwald = &I->I[is].FEwald[NDIM*ia];
567 FIon_old = &I->I[is].FIon_old[NDIM*ia];
568 U = &I->I[is].U[NDIM*ia];
569 R = &I->I[is].R[NDIM*ia];
570 R_old = &I->I[is].R_old[NDIM*ia];
571 R_old_old = &I->I[is].R_old_old[NDIM*ia];
572 if (I->I[is].IMT[ia] == MoveIon) {
573 for (d=0; d<NDIM; d++) {
574 R_old_old[d] = R_old[d]; // shift old values
575 R_old[d] = R[d];
576 R[d] += delta_t*(U[d]+a*FIon[d]); // F = m * a and s = 0.5 * a * t^2
577 FIon_old[d] = FIon[d]; // store old force
578 FIon[d] = 0; // zero all as a sign that's been moved
579 FIonL[d] = 0;
580 FIonNL[d] = 0;
581 FEwald[d] = 0;
582 }
583 RMat33Vec3(sR, P->Lat.InvBasis, R);
584 RMat33Vec3(sRold, P->Lat.InvBasis, R_old);
585 RMat33Vec3(sRoldold, P->Lat.InvBasis, R_old_old);
586 for (d=0; d <NDIM; d++) {
587 while (sR[d] < 0) {
588 sR[d] += 1.0;
589 sRold[d] += 1.0;
590 sRoldold[d] += 1.0;
591 }
592 while (sR[d] >= 1.0) {
593 sR[d] -= 1.0;
594 sRold[d] -= 1.0;
595 sRoldold[d] -= 1.0;
596 }
597 }
598 RMat33Vec3(R, P->Lat.RealBasis, sR);
599 RMat33Vec3(R_old, P->Lat.RealBasis, sRold);
600 RMat33Vec3(R_old_old, P->Lat.RealBasis, sRoldold);
601 }
602 }
603 }
604}
605
606/** Resets the forces array.
607 * \param *P Problem at hand
608 */
609void ResetForces(struct Problem *P)
610{
611 struct Ions *I = &P->Ion;
612 double *FIon, *FIonL, *FIonNL, *FEwald;
613 int is, ia, d;
614 for (is=0; is < I->Max_Types; is++) { // go through all types ...
615 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
616 FIon = &I->I[is].FIon[NDIM*ia];
617 FIonL = &I->I[is].FIonL[NDIM*ia];
618 FIonNL = &I->I[is].FIonNL[NDIM*ia];
619 FEwald = &I->I[is].FEwald[NDIM*ia];
620 for (d=0; d<NDIM; d++) {
621 FIon[d] = 0.; // zero all as a sign that's been moved
622 FIonL[d] = 0.;
623 FIonNL[d] = 0.;
624 FEwald[d] = 0.;
625 }
626 }
627 }
628};
629
630/** Changes ion's velocity according to acting force.
631 * IonType::U is changed by the weighted force of actual step and last one
632 * according to Newton.
633 * \param *P Problem at hand
634 * \sa UpdateIonsR
635 * \sa CalculateIonForce
636 * \warning R is not changed here.
637 */
638void UpdateIonsU(struct Problem *P)
639{
640 struct Ions *I = &P->Ion;
641 int is, ia, d;
642 double *FIon, *FIon_old, *U;
643 double IonMass, a;
644 const int delta_t = P->R.delta_t;
645 for (is=0; is < I->Max_Types; is++) {
646 IonMass = I->I[is].IonMass;
647 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsU: IonMass = %lg",IonMass);
648 a = delta_t*0.5/IonMass; // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
649 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
650 FIon = &I->I[is].FIon[NDIM*ia];
651 FIon_old = &I->I[is].FIon_old[NDIM*ia];
652 U = &I->I[is].U[NDIM*ia];
653 if (I->I[is].IMT[ia] == MoveIon)
654 for (d=0; d<NDIM; d++) {
655 U[d] += a * (FIon[d]+FIon_old[d]);
656 }
657 }
658 }
659}
660
661/** CG process to optimise structure.
662 * \param *P Problem at hand
663 */
664void UpdateIons(struct Problem *P)
665{
666 struct Ions *I = &P->Ion;
667 int is, ia, d;
668 double *R, *R_old, *R_old_old, *FIon, *S, *GammaA;
669 double IonFac, GammaB; //, GammaT;
670 double FNorm, StepNorm;
671 for (is=0; is < I->Max_Types; is++) {
672 IonFac = I->I[is].IonFac;
673 GammaA = I->I[is].GammaA;
674 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
675 FIon = &I->I[is].FIon[NDIM*ia];
676 S = &I->I[is].SearchDir[NDIM*ia];
677 GammaB = GammaA[ia];
678 GammaA[ia] = RSP3(FIon,S);
679 FNorm = sqrt(RSP3(FIon,FIon));
680 StepNorm = log(1.+IonFac*FNorm); // Fix Hack
681 if (FNorm != 0)
682 StepNorm /= sqrt(RSP3(FIon,FIon));
683 else
684 StepNorm = 0;
685 //if (GammaB != 0.0) {
686 // GammaT = GammaA[ia]/GammaB;
687 // for (d=0; d>NDIM; d++)
688 // S[d] = FIon[d]+GammaT*S[d];
689 // } else {
690 for (d=0; d<NDIM; d++)
691 S[d] = StepNorm*FIon[d];
692 // }
693 R = &I->I[is].R[NDIM*ia];
694 R_old = &I->I[is].R_old[NDIM*ia];
695 R_old_old = &I->I[is].R_old_old[NDIM*ia];
696 if (I->I[is].IMT[ia] == MoveIon)
697 for (d=0; d<NDIM; d++) {
698 R_old_old[d] = R_old[d];
699 R_old[d] = R[d];
700 R[d] += S[d]; // FixHack*IonFac;
701 }
702 }
703 }
704}
705
706/** Print coordinates of all ions to stdout.
707 * \param *P Problem at hand
708 */
709void OutputIonCoordinates(struct Problem *P)
710{
711 struct Ions *I = &P->Ion;
712 int is, ia;
713 if (P->Par.me == 0) {
714 fprintf(stderr, "(%i) ======= Updated Ion Coordinates =========\n",P->Par.me);
715 for (is=0; is < I->Max_Types; is++)
716 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++)
717 //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]);
718 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[NDIM*ia+1],I->I[is].R[NDIM*ia+2]);
719 fprintf(stderr, "(%i) =========================================\n",P->Par.me);
720 }
721}
722
723/** Calculates kinetic energy (Ions::EKin) of movable Ions.
724 * Does 0.5 / IonType::IonMass * IonTye::U^2 for each Ion,
725 * also retrieves actual temperatur by 2/3 from just
726 * calculated Ions::EKin.
727 * \param *P Problem at hand
728 */
729void CalculateEnergyIonsU(struct Problem *P)
730{
731 struct Ions *I = &P->Ion;
732 int is, ia, d;
733 double *U;
734 double IonMass, a, ekin = 0;
735 for (is=0; is < I->Max_Types; is++) {
736 IonMass = I->I[is].IonMass;
737 a = 0.5*IonMass;
738 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
739 U = &I->I[is].U[NDIM*ia];
740 if (I->I[is].IMT[ia] == MoveIon)
741 for (d=0; d<NDIM; d++) {
742 ekin += a * U[d]*U[d];
743 }
744 }
745 }
746 I->EKin = ekin;
747 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
748}
749
750/** Scales ion velocities to match temperature.
751 * In order to match Ions::ActualTemp with desired RunStruct::TargetTemp
752 * each velocity in each dimension (for all types, all ions) is scaled
753 * with the ratio of these two. Ions::EKin is recalculated.
754 * \param *P Problem at hand
755 */
756void ScaleTemp(struct Problem *P)
757{
758 struct Ions *I = &P->Ion;
759 int is, ia, d;
760 double *U;
761 double IonMass, a, ekin = 0;
762 if (I->ActualTemp < MYEPSILON) fprintf(stderr,"ScaleTemp: I->ActualTemp = %lg",I->ActualTemp);
763 double ScaleTempFactor = sqrt(P->R.TargetTemp/I->ActualTemp);
764 for (is=0; is < I->Max_Types; is++) {
765 IonMass = I->I[is].IonMass;
766 a = 0.5*IonMass;
767 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
768 U = &I->I[is].U[NDIM*ia];
769 if (I->I[is].IMT[ia] == MoveIon)
770 for (d=0; d<NDIM; d++) {
771 U[d] *= ScaleTempFactor;
772 ekin += a * U[d]*U[d];
773 }
774 }
775 }
776 I->EKin = ekin;
777 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
778}
779
780/** Calculates mean force vector as stop criteria in structure optimization.
781 * Calculates a mean force vector per ion as the usual euclidian norm over total
782 * number of ions and dimensions, being the sum of each ion force (for all type,
783 * all ions, all dims) squared.
784 * The mean force is archived in RunStruct::MeanForce and printed to screen.
785 * \param *P Problem at hand
786 */
787void GetOuterStop(struct Problem *P)
788{
789 struct RunStruct *R = &P->R;
790 struct Ions *I = &P->Ion;
791 int is, ia, IonNo=0, i, d;
792 double MeanForce = 0.0;
793 for (is=0; is < I->Max_Types; is++)
794 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
795 IonNo++;
796 for (d=0; d <NDIM; d++)
797 MeanForce += I->I[is].FIon[d+NDIM*ia]*I->I[is].FIon[d+NDIM*ia];
798 }
799 for (i=MAXOLD-1; i > 0; i--)
800 R->MeanForce[i] = R->MeanForce[i-1];
801 MeanForce = sqrt(MeanForce/(IonNo*NDIM));
802 R->MeanForce[0] = MeanForce;
803 if (P->Call.out[LeaderOut] && (P->Par.me == 0))
804 fprintf(stderr,"MeanForce = %e\n", R->MeanForce[0]);
805}
806
807
Note: See TracBrowser for help on using the repository browser.