source: pcp/src/ions.c@ 460e95

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

BOHRRADIUS changed to ANGSTROEMINBOHRRADIUS
  • Property mode set to 100644
File size: 29.0 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 int k;
64 struct Lattice *L = &P->Lat;
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);
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/** Calculated Ewald force.
193 * \f[
194 * 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|}
195 * \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 )
196 * \qquad (4.10)
197 * \f]
198 * Calculates the core-to-core force between the ions of all super cells.
199 * In order for this series to converge there must be a certain summation
200 * applied, which is the ewald summation and which is nothing else but to move
201 * in a circular motion from the current cell to the outside up to Ions::R_cut.
202 * \param *P Problem at hand
203 * \param first additional calculation beforehand if != 0
204 */
205void CalculateEwald(struct Problem *P, int first)
206{
207 int r,is,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec;
208 int MaxPar=P->Par.procs, ParMe=P->Par.me, MaxLocalVec, StartVec;
209 double rcut2,r2,erre2,addesr,addpre,arg,fac,gkl,esr,fxx,repand;
210 double R1[3];
211 double R2[3];
212 double R3[3];
213 double t[3];
214 struct Lattice *L = &P->Lat;
215 struct PseudoPot *PP = &P->PP;
216 struct Ions *I = &P->Ion;
217 if (I->R_cut != 0.0) {
218 if (first) {
219 SpeedMeasure(P, InitSimTime, StopTimeDo);
220 rcut2 = I->R_cut*I->R_cut;
221 ActualVec =0;
222 MaxCell = (int)5.*I->R_cut/pow(L->Volume,1./3.);
223 if (MaxCell < 2) MaxCell = 2;
224 for (i = -MaxCell; i <= MaxCell; i++) {
225 for (j = -MaxCell; j <= MaxCell; j++) {
226 for (k = -MaxCell; k <= MaxCell; k++) {
227 r2 = 0.0;
228 for (ir=0; ir <3; ir++) {
229 t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
230 r2 = t[ir]*t[ir];
231 }
232 Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
233 if ((r2 <= rcut2) || Is_Neighb_Cell) {
234 ActualVec++;
235 }
236 }
237 }
238 }
239 MaxVec = ActualVec;
240 I->MaxVec = ActualVec;
241 MaxLocalVec = MaxVec / MaxPar;
242 StartVec = ParMe*MaxLocalVec;
243 r = MaxVec % MaxPar;
244 if (ParMe >= r) {
245 StartVec += r;
246 } else {
247 StartVec += ParMe;
248 MaxLocalVec++;
249 }
250 I->MaxLocalVec = MaxLocalVec;
251 LocalActualVec = ActualVec = 0;
252 I->RLatticeVec = (double *)
253 Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald:");
254 for (i = -MaxCell; i <= MaxCell; i++) {
255 for (j = -MaxCell; j <= MaxCell; j++) {
256 for (k = -MaxCell; k <= MaxCell; k++) {
257 r2 = 0.0;
258 for (ir=0; ir <3; ir++) {
259 t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
260 r2 = t[ir]*t[ir];
261 }
262 Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
263 if ((r2 <= rcut2) || Is_Neighb_Cell) {
264 if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) {
265 I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
266 I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
267 I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
268 LocalActualVec++;
269 }
270 ActualVec++;
271 }
272 }
273 }
274 }
275 SpeedMeasure(P, InitSimTime, StartTimeDo);
276 }
277 SpeedMeasure(P, EwaldTime, StartTimeDo);
278 esr = 0.0;
279 for (is1=0;is1 < I->Max_Types; is1++)
280 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
281 for (i=0; i< NDIM; i++)
282 I->FTemp[i+NDIM*(ia1)] = 0;
283 for (is1=0;is1 < I->Max_Types; is1++) {
284 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
285 for (i=0;i<3;i++) {
286 R1[i]=I->I[is1].R[i+NDIM*ia1];
287 }
288 for (ir2=0;ir2<I->MaxLocalVec;ir2++) {
289 for (is2=0;is2 < I->Max_Types; is2++) {
290 gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
291 for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
292 for (i=0;i<3;i++) {
293 R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2];
294 }
295 erre2=0.0;
296 for (i=0;i<3;i++) {
297 R3[i] = R1[i]-R2[i];
298 erre2+=R3[i]*R3[i];
299 }
300 if (erre2 > MYEPSILON) {
301 arg=sqrt(erre2);
302 fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;
303
304 arg *= gkl;
305 addesr = derf(arg);
306 addesr = (1.0-addesr)*fac;
307 esr += addesr;
308 addpre=exp(-arg*arg)*gkl;
309 addpre=PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;
310 repand=(addesr+addpre)/erre2;
311 for (i=0;i<3;i++) {
312 fxx=repand*R3[i];
313 /*fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);*/
314 I->FTemp[i+NDIM*(ia1)] += fxx;
315 I->FTemp[i+NDIM*(ia2)] -= fxx;
316 }
317 }
318 }
319 }
320 }
321 }
322 }
323 } else {
324 esr = 0.0;
325 for (is1=0;is1 < I->Max_Types; is1++)
326 for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
327 for (i=0; i< NDIM; i++)
328 I->FTemp[i+NDIM*(ia1)] = 0.;
329 }
330 SpeedMeasure(P, EwaldTime, StopTimeDo);
331 for (is=0;is < I->Max_Types; is++) {
332 MPI_Allreduce (I->FTemp , I->I[is].FEwald, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);
333 }
334 MPI_Allreduce ( &esr, &L->E->AllTotalIonsEnergy[EwaldEnergy], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm);
335}
336
337/** Sum up all forces acting on Ions.
338 * IonType::FIon = IonType::FEwald + IonType::FIonL + IonType::FIonNL for all
339 * dimensions #NDIM and each Ion of IonType
340 * \param *P Problem at hand
341 */
342void CalculateIonForce(struct Problem *P)
343{
344 struct Ions *I = &P->Ion;
345 int i,j,d;
346 for (i=0; i < I->Max_Types; i++)
347 for (j=0; j < I->I[i].Max_IonsOfType; j++)
348 for (d=0; d<NDIM;d++)
349 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];
350}
351
352/** Write Forces to FileData::ForcesFile.
353 * goes through all Ions per IonType and each dimension of #NDIM and prints in one line:
354 * Position IonType::R, total force Ion IonType::FIon, local force IonType::FIonL,
355 * non-local IonType::FIonNL, ewald force IonType::FEwald
356 * \param *P Problem at hand
357 */
358void OutputIonForce(struct Problem *P)
359{
360 struct RunStruct *R = &P->R;
361 struct FileData *F = &P->Files;
362 struct Ions *I = &P->Ion;
363 FILE *fout = NULL;
364 int i,j;
365 if (F->MeOutMes != 1) return;
366 fout = F->ForcesFile;
367 for (i=0; i < I->Max_Types; i++)
368 for (j=0; j < I->I[i].Max_IonsOfType; j++)
369 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,
370 I->I[i].R[0+j*NDIM], I->I[i].R[1+j*NDIM], I->I[i].R[2+j*NDIM],
371 I->I[i].FIon[0+j*NDIM], I->I[i].FIon[1+j*NDIM], I->I[i].FIon[2+j*NDIM],
372 I->I[i].FIonL[0+j*NDIM], I->I[i].FIonL[1+j*NDIM], I->I[i].FIonL[2+j*NDIM],
373 I->I[i].FIonNL[0+j*NDIM], I->I[i].FIonNL[1+j*NDIM], I->I[i].FIonNL[2+j*NDIM],
374 I->I[i].FEwald[0+j*NDIM], I->I[i].FEwald[1+j*NDIM], I->I[i].FEwald[2+j*NDIM]);
375 fflush(fout);
376 fprintf(fout, "MeanForce:\t%e\n", R->MeanForce[0]);
377}
378
379/** Frees memory IonType.
380 * All pointers from IonType and the pointer on it are free'd
381 * \param *I Ions structure to be free'd
382 * \sa IonsInitRead where memory is allocated
383 */
384void RemoveIonsRead(struct Ions *I)
385{
386 int i, it;
387 for (i=0; i < I->Max_Types; i++) {
388 //fprintf(stderr,"Name\n");
389 Free(I->I[i].Name);
390 //fprintf(stderr,"Abbreviation\n");
391 Free(I->I[i].Symbol);
392 for (it=0;it<I->I[i].Max_IonsOfType;it++) {
393 //fprintf(stderr,"sigma[it]\n");
394 Free(I->I[i].sigma[it]);
395 //fprintf(stderr,"sigma_rezi[it]\n");
396 Free(I->I[i].sigma_rezi[it]);
397 //fprintf(stderr,"sigma_PAS[it]\n");
398 Free(I->I[i].sigma_PAS[it]);
399 //fprintf(stderr,"sigma_rezi_PAS[it]\n");
400 Free(I->I[i].sigma_rezi_PAS[it]);
401 }
402 //fprintf(stderr,"sigma\n");
403 Free(I->I[i].sigma);
404 //fprintf(stderr,"sigma_rezi\n");
405 Free(I->I[i].sigma_rezi);
406 //fprintf(stderr,"sigma_PAS\n");
407 Free(I->I[i].sigma_PAS);
408 //fprintf(stderr,"sigma_rezi_PAS\n");
409 Free(I->I[i].sigma_rezi_PAS);
410 //fprintf(stderr,"R\n");
411 Free(I->I[i].R);
412 //fprintf(stderr,"R_old\n");
413 Free(I->I[i].R_old);
414 //fprintf(stderr,"R_old_old\n");
415 Free(I->I[i].R_old_old);
416 //fprintf(stderr,"FIon\n");
417 Free(I->I[i].FIon);
418 //fprintf(stderr,"FIon_old\n");
419 Free(I->I[i].FIon_old);
420 //fprintf(stderr,"SearchDir\n");
421 Free(I->I[i].SearchDir);
422 //fprintf(stderr,"GammaA\n");
423 Free(I->I[i].GammaA);
424 //fprintf(stderr,"FIonL\n");
425 Free(I->I[i].FIonL);
426 //fprintf(stderr,"FIonNL\n");
427 Free(I->I[i].FIonNL);
428 //fprintf(stderr,"U\n");
429 Free(I->I[i].U);
430 //fprintf(stderr,"SFactor\n");
431 Free(I->I[i].SFactor);
432 //fprintf(stderr,"IMT\n");
433 Free(I->I[i].IMT);
434 //fprintf(stderr,"FEwald\n");
435 Free(I->I[i].FEwald);
436 Free(I->I[i].alpha);
437 }
438 //fprintf(stderr,"RLatticeVec\n");
439 if (I->R_cut == 0.0)
440 Free(I->RLatticeVec);
441 //fprintf(stderr,"FTemp\n");
442 Free(I->FTemp);
443 //fprintf(stderr,"I\n");
444 Free(I->I);
445}
446
447/** Shifts center of gravity of ion forces IonType::FIon.
448 * First sums up all forces of movable ions to a "force temperature",
449 * then reduces each force by this temp, so that all in all
450 * the net force is 0.
451 * \param *P Problem at hand
452 * \sa CorrectVelocity
453 * \note why is FTemp not divided by number of ions: is probably correct, but this
454 * function was not really meaningful from the beginning.
455 */
456void CorrectForces(struct Problem *P)
457{
458 struct Ions *I = &P->Ion;
459 double *FIon;
460 double FTemp[NDIM];
461 int is, ia, d;
462 for (d=0; d<NDIM; d++)
463 FTemp[d] = 0;
464 for (is=0; is < I->Max_Types; is++) { // calculate force temperature for each type ...
465 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of this type ...
466 FIon = &I->I[is].FIon[NDIM*ia];
467 if (I->I[is].IMT[ia] == MoveIon) { // .. if it's movable
468 for (d=0; d<NDIM; d++)
469 FTemp[d] += FIon[d];
470 }
471 }
472 }
473 for (is=0; is < I->Max_Types; is++) { // and then reduce each by this value
474 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
475 FIon = &I->I[is].FIon[NDIM*ia];
476 if (I->I[is].IMT[ia] == MoveIon) {
477 for (d=0; d<NDIM; d++)
478 FIon[d] -= FTemp[d]; // (?) why not -= FTemp[d]/I->Max_TotalIons
479 }
480 }
481 }
482}
483
484
485/** Shifts center of gravity of ion velocities (rather momentums).
486 * First sums up ion speed U times their IonMass (summed up for each
487 * dimension), then reduces the velocity by temp per Ions::TotalMass (so here
488 * total number of Ions is included)
489 * \param *P Problem at hand
490 * \sa CorrectForces
491 */
492void CorrectVelocity(struct Problem *P)
493{
494 struct Ions *I = &P->Ion;
495 double *U;
496 double UTemp[NDIM];
497 int is, ia, d;
498 for (d=0; d<NDIM; d++)
499 UTemp[d] = 0;
500 for (is=0; is < I->Max_Types; is++) { // all types ...
501 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // ... all ions per type ...
502 U = &I->I[is].U[NDIM*ia];
503 if (I->I[is].IMT[ia] == MoveIon) { // ... if it's movable
504 for (d=0; d<NDIM; d++)
505 UTemp[d] += I->I[is].IonMass*U[d];
506 }
507 }
508 }
509 for (is=0; is < I->Max_Types; is++) {
510 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
511 U = &I->I[is].U[NDIM*ia];
512 if (I->I[is].IMT[ia] == MoveIon) {
513 for (d=0; d<NDIM; d++)
514 U[d] -= UTemp[d]/I->TotalMass; // shift by mean velocity over mass and number of ions
515 }
516 }
517 }
518}
519
520/** Moves ions according to calculated force.
521 * Goes through each type of IonType and each ion therein, takes mass and
522 * Newton to move each ion to new coordinates IonType::R, while remembering the last
523 * two coordinates and the last force of the ion. Coordinates R are being
524 * transformed to inverted base, shifted by +-1.0 and back-transformed to
525 * real base.
526 * \param *P Problem at hand
527 * \sa CalculateIonForce
528 * \sa UpdateIonsU
529 * \warning U is not changed here, only used to move the ion.
530 */
531void UpdateIonsR(struct Problem *P)
532{
533 struct Ions *I = &P->Ion;
534 int is, ia, d;
535 double *R, *R_old, *R_old_old, *FIon, *FIon_old, *U, *FIonL, *FIonNL, *FEwald;
536 double IonMass, a;
537 double sR[NDIM], sRold[NDIM], sRoldold[NDIM];
538 const int delta_t = P->R.delta_t;
539 for (is=0; is < I->Max_Types; is++) { // go through all types ...
540 IonMass = I->I[is].IonMass;
541 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsR: IonMass = %lg",IonMass);
542 a = delta_t*0.5/IonMass; // F/m = a and thus: s = 0.5 * F/m * t^2 + v * t =: t * (F * a + v)
543 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
544 FIon = &I->I[is].FIon[NDIM*ia];
545 FIonL = &I->I[is].FIonL[NDIM*ia];
546 FIonNL = &I->I[is].FIonNL[NDIM*ia];
547 FEwald = &I->I[is].FEwald[NDIM*ia];
548 FIon_old = &I->I[is].FIon_old[NDIM*ia];
549 U = &I->I[is].U[NDIM*ia];
550 R = &I->I[is].R[NDIM*ia];
551 R_old = &I->I[is].R_old[NDIM*ia];
552 R_old_old = &I->I[is].R_old_old[NDIM*ia];
553 if (I->I[is].IMT[ia] == MoveIon) {
554 for (d=0; d<NDIM; d++) {
555 R_old_old[d] = R_old[d]; // shift old values
556 R_old[d] = R[d];
557 R[d] += delta_t*(U[d]+a*FIon[d]); // F = m * a and s = 0.5 * a * t^2
558 FIon_old[d] = FIon[d]; // store old force
559 FIon[d] = 0; // zero all as a sign that's been moved
560 FIonL[d] = 0;
561 FIonNL[d] = 0;
562 FEwald[d] = 0;
563 }
564 RMat33Vec3(sR, P->Lat.InvBasis, R);
565 RMat33Vec3(sRold, P->Lat.InvBasis, R_old);
566 RMat33Vec3(sRoldold, P->Lat.InvBasis, R_old_old);
567 for (d=0; d <NDIM; d++) {
568 while (sR[d] < 0) {
569 sR[d] += 1.0;
570 sRold[d] += 1.0;
571 sRoldold[d] += 1.0;
572 }
573 while (sR[d] >= 1.0) {
574 sR[d] -= 1.0;
575 sRold[d] -= 1.0;
576 sRoldold[d] -= 1.0;
577 }
578 }
579 RMat33Vec3(R, P->Lat.RealBasis, sR);
580 RMat33Vec3(R_old, P->Lat.RealBasis, sRold);
581 RMat33Vec3(R_old_old, P->Lat.RealBasis, sRoldold);
582 }
583 }
584 }
585}
586
587/** Resets the forces array.
588 * \param *P Problem at hand
589 */
590void ResetForces(struct Problem *P)
591{
592 struct Ions *I = &P->Ion;
593 double *FIon, *FIonL, *FIonNL, *FEwald;
594 int is, ia, d;
595 for (is=0; is < I->Max_Types; is++) { // go through all types ...
596 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) { // .. and each ion of the type
597 FIon = &I->I[is].FIon[NDIM*ia];
598 FIonL = &I->I[is].FIonL[NDIM*ia];
599 FIonNL = &I->I[is].FIonNL[NDIM*ia];
600 FEwald = &I->I[is].FEwald[NDIM*ia];
601 for (d=0; d<NDIM; d++) {
602 FIon[d] = 0.; // zero all as a sign that's been moved
603 FIonL[d] = 0.;
604 FIonNL[d] = 0.;
605 FEwald[d] = 0.;
606 }
607 }
608 }
609};
610
611/** Changes ion's velocity according to acting force.
612 * IonType::U is changed by the weighted force of actual step and last one
613 * according to Newton.
614 * \param *P Problem at hand
615 * \sa UpdateIonsR
616 * \sa CalculateIonForce
617 * \warning R is not changed here.
618 */
619void UpdateIonsU(struct Problem *P)
620{
621 struct Ions *I = &P->Ion;
622 int is, ia, d;
623 double *FIon, *FIon_old, *U;
624 double IonMass, a;
625 const int delta_t = P->R.delta_t;
626 for (is=0; is < I->Max_Types; is++) {
627 IonMass = I->I[is].IonMass;
628 if (IonMass < MYEPSILON) fprintf(stderr,"UpdateIonsU: IonMass = %lg",IonMass);
629 a = delta_t*0.5/IonMass; // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
630 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
631 FIon = &I->I[is].FIon[NDIM*ia];
632 FIon_old = &I->I[is].FIon_old[NDIM*ia];
633 U = &I->I[is].U[NDIM*ia];
634 if (I->I[is].IMT[ia] == MoveIon)
635 for (d=0; d<NDIM; d++) {
636 U[d] += a * (FIon[d]+FIon_old[d]);
637 }
638 }
639 }
640}
641
642/** CG process to optimise structure.
643 * \param *P Problem at hand
644 */
645void UpdateIons(struct Problem *P)
646{
647 struct Ions *I = &P->Ion;
648 int is, ia, d;
649 double *R, *R_old, *R_old_old, *FIon, *S, *GammaA;
650 double IonFac, GammaB; //, GammaT;
651 double FNorm, StepNorm;
652 for (is=0; is < I->Max_Types; is++) {
653 IonFac = I->I[is].IonFac;
654 GammaA = I->I[is].GammaA;
655 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
656 FIon = &I->I[is].FIon[NDIM*ia];
657 S = &I->I[is].SearchDir[NDIM*ia];
658 GammaB = GammaA[ia];
659 GammaA[ia] = RSP3(FIon,S);
660 FNorm = sqrt(RSP3(FIon,FIon));
661 StepNorm = log(1.+IonFac*FNorm); // Fix Hack
662 if (FNorm != 0)
663 StepNorm /= sqrt(RSP3(FIon,FIon));
664 else
665 StepNorm = 0;
666 //if (GammaB != 0.0) {
667 // GammaT = GammaA[ia]/GammaB;
668 // for (d=0; d>NDIM; d++)
669 // S[d] = FIon[d]+GammaT*S[d];
670 // } else {
671 for (d=0; d<NDIM; d++)
672 S[d] = StepNorm*FIon[d];
673 // }
674 R = &I->I[is].R[NDIM*ia];
675 R_old = &I->I[is].R_old[NDIM*ia];
676 R_old_old = &I->I[is].R_old_old[NDIM*ia];
677 if (I->I[is].IMT[ia] == MoveIon)
678 for (d=0; d<NDIM; d++) {
679 R_old_old[d] = R_old[d];
680 R_old[d] = R[d];
681 R[d] += S[d]; // FixHack*IonFac;
682 }
683 }
684 }
685}
686
687/** Print coordinates of all ions to stdout.
688 * \param *P Problem at hand
689 */
690void OutputIonCoordinates(struct Problem *P)
691{
692 struct Ions *I = &P->Ion;
693 int is, ia;
694 if (P->Par.me == 0) {
695 fprintf(stderr, "(%i) ======= Updated Ion Coordinates =========\n",P->Par.me);
696 for (is=0; is < I->Max_Types; is++)
697 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++)
698 //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]);
699 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]);
700 fprintf(stderr, "(%i) =========================================\n",P->Par.me);
701 }
702}
703
704/** Calculates kinetic energy (Ions::EKin) of movable Ions.
705 * Does 0.5 / IonType::IonMass * IonTye::U^2 for each Ion,
706 * also retrieves actual temperatur by 2/3 from just
707 * calculated Ions::EKin.
708 * \param *P Problem at hand
709 */
710void CalculateEnergyIonsU(struct Problem *P)
711{
712 struct Ions *I = &P->Ion;
713 int is, ia, d;
714 double *U;
715 double IonMass, a, ekin = 0;
716 for (is=0; is < I->Max_Types; is++) {
717 IonMass = I->I[is].IonMass;
718 a = 0.5*IonMass;
719 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
720 U = &I->I[is].U[NDIM*ia];
721 if (I->I[is].IMT[ia] == MoveIon)
722 for (d=0; d<NDIM; d++) {
723 ekin += a * U[d]*U[d];
724 }
725 }
726 }
727 I->EKin = ekin;
728 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
729}
730
731/** Scales ion velocities to match temperature.
732 * In order to match Ions::ActualTemp with desired RunStruct::TargetTemp
733 * each velocity in each dimension (for all types, all ions) is scaled
734 * with the ratio of these two. Ions::EKin is recalculated.
735 * \param *P Problem at hand
736 */
737void ScaleTemp(struct Problem *P)
738{
739 struct Ions *I = &P->Ion;
740 int is, ia, d;
741 double *U;
742 double IonMass, a, ekin = 0;
743 if (I->ActualTemp < MYEPSILON) fprintf(stderr,"ScaleTemp: I->ActualTemp = %lg",I->ActualTemp);
744 double ScaleTempFactor = sqrt(P->R.TargetTemp/I->ActualTemp);
745 for (is=0; is < I->Max_Types; is++) {
746 IonMass = I->I[is].IonMass;
747 a = 0.5*IonMass;
748 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
749 U = &I->I[is].U[NDIM*ia];
750 if (I->I[is].IMT[ia] == MoveIon)
751 for (d=0; d<NDIM; d++) {
752 U[d] *= ScaleTempFactor;
753 ekin += a * U[d]*U[d];
754 }
755 }
756 }
757 I->EKin = ekin;
758 I->ActualTemp = (2./(3.*I->Max_TotalIons)*I->EKin);
759}
760
761/** Calculates mean force vector as stop criteria in structure optimization.
762 * Calculates a mean force vector per ion as the usual euclidian norm over total
763 * number of ions and dimensions, being the sum of each ion force (for all type,
764 * all ions, all dims) squared.
765 * The mean force is archived in RunStruct::MeanForce and printed to screen.
766 * \param *P Problem at hand
767 */
768void GetOuterStop(struct Problem *P)
769{
770 struct RunStruct *R = &P->R;
771 struct Ions *I = &P->Ion;
772 int is, ia, IonNo=0, i, d;
773 double MeanForce = 0.0;
774 for (is=0; is < I->Max_Types; is++)
775 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
776 IonNo++;
777 for (d=0; d <NDIM; d++)
778 MeanForce += I->I[is].FIon[d+NDIM*ia]*I->I[is].FIon[d+NDIM*ia];
779 }
780 for (i=MAXOLD-1; i > 0; i--)
781 R->MeanForce[i] = R->MeanForce[i-1];
782 MeanForce = sqrt(MeanForce/(IonNo*NDIM));
783 R->MeanForce[0] = MeanForce;
784 if (P->Call.out[LeaderOut] && (P->Par.me == 0))
785 fprintf(stderr,"MeanForce = %e\n", R->MeanForce[0]);
786}
787
788
Note: See TracBrowser for help on using the repository browser.