source: pcp/src/ions.c@ 64ca279

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

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