1 #include "dynmat.h"
2 #include "math.h"
3 #include "version.h"
4 #include "global.h"
5
6 /* ----------------------------------------------------------------------------
7 * Class DynMat stores the Dynamic Matrix read from the binary file from
8 * fix-phonon and incorporates some operations w.r.t the DM.
9 * ---------------------------------------------------------------------------- */
DynMat(int narg,char ** arg)10 DynMat::DynMat(int narg, char **arg)
11 {
12 attyp = NULL;
13 memory = NULL;
14 M_inv_sqrt = NULL;
15 interpolate = NULL;
16 DM_q = DM_all = NULL;
17 binfile = funit = dmfile = NULL;
18
19 attyp = NULL;
20 basis = NULL;
21 flag_reset_gamma = flag_skip = 0;
22
23 // analyze the command line options
24 int iarg = 1;
25 while (narg > iarg){
26 if (strcmp(arg[iarg], "-s") == 0){
27 flag_reset_gamma = flag_skip = 1;
28
29 } else if (strcmp(arg[iarg], "-r") == 0){
30 flag_reset_gamma = 1;
31
32 } else if (strcmp(arg[iarg], "-h") == 0){
33 help();
34
35 } else {
36 if (binfile) delete []binfile;
37 int n = strlen(arg[iarg]) + 1;
38 binfile = new char[n];
39 strcpy(binfile, arg[iarg]);
40 }
41
42 iarg++;
43 }
44
45 ShowVersion();
46 // get the binary file name from user input if not found in command line
47 char str[MAXLINE];
48 if (binfile == NULL) {
49 char *ptr = NULL;
50 printf("\n");
51 do {
52 printf("Please input the binary file name from fix_phonon: ");
53 fgets(str,MAXLINE,stdin);
54 ptr = strtok(str, " \n\t\r\f");
55 } while (ptr == NULL);
56
57 int n = strlen(ptr) + 1;
58 binfile = new char[n];
59 strcpy(binfile, ptr);
60 }
61
62 // open the binary file
63 FILE *fp = fopen(binfile, "rb");
64 if (fp == NULL) {
65 printf("\nFile %s not found! Programe terminated.\n", binfile);
66 help();
67 }
68
69 // read header info from the binary file
70 if (fread(&sysdim, sizeof(int), 1, fp) != 1){
71 printf("\nError while reading sysdim from file: %s\n", binfile);
72 fclose(fp);
73 exit(2);
74 }
75 if (fread(&nx, sizeof(int), 1, fp) != 1){
76 printf("\nError while reading nx from file: %s\n", binfile);
77 fclose(fp);
78 exit(2);
79 }
80 if (fread(&ny, sizeof(int), 1, fp) != 1){
81 printf("\nError while reading ny from file: %s\n", binfile);
82 fclose(fp);
83 exit(2);
84 }
85 if (fread(&nz, sizeof(int), 1, fp) != 1){
86 printf("\nError while reading nz from file: %s\n", binfile);
87 fclose(fp);
88 exit(2);
89 }
90 if (fread(&nucell, sizeof(int), 1, fp) != 1){
91 printf("\nError while reading nucell from file: %s\n", binfile);
92 fclose(fp);
93 exit(2);
94 }
95 if (fread(&boltz, sizeof(double), 1, fp) != 1){
96 printf("\nError while reading boltz from file: %s\n", binfile);
97 fclose(fp);
98 exit(2);
99 }
100
101 fftdim = sysdim * nucell;
102 fftdim2 = fftdim * fftdim;
103 npt = nx * ny * nz;
104
105 // display info related to the read file
106 ShowInfo();
107 if (sysdim < 1||sysdim > 3||nx < 1||ny < 1||nz < 1||nucell < 1){
108 printf("Wrong values read from header of file: %s, please check the binary file!\n", binfile);
109 fclose(fp); exit(3);
110 }
111 Define_Conversion_Factor();
112
113 // now to allocate memory for DM
114 memory = new Memory();
115 memory->create(DM_all, npt, fftdim2, "DynMat:DM_all");
116 memory->create(DM_q, fftdim,fftdim,"DynMat:DM_q");
117
118 // read all dynamical matrix info into DM_all
119 if (fread(DM_all[0], sizeof(doublecomplex), npt*fftdim2, fp) != size_t(npt*fftdim2)){
120 printf("\nError while reading the DM from file: %s\n", binfile);
121 fclose(fp);
122 exit(1);
123 }
124
125 // now try to read unit cell info from the binary file
126 memory->create(basis, nucell, sysdim, "DynMat:basis");
127 memory->create(attyp, nucell, "DynMat:attyp");
128 memory->create(M_inv_sqrt, nucell, "DynMat:M_inv_sqrt");
129
130 if (fread(&Tmeasure, sizeof(double), 1, fp) != 1){
131 printf("\nError while reading temperature from file: %s\n", binfile);
132 fclose(fp);
133 exit(3);
134 }
135 if (fread(&basevec[0], sizeof(double), 9, fp) != 9){
136 printf("\nError while reading lattice info from file: %s\n", binfile);
137 fclose(fp);
138 exit(3);
139 }
140 if (fread(basis[0], sizeof(double), fftdim, fp) != fftdim){
141 printf("\nError while reading basis info from file: %s\n", binfile);
142 fclose(fp);
143 exit(3);
144 }
145 if (fread(&attyp[0], sizeof(int), nucell, fp) != nucell){
146 printf("\nError while reading atom types from file: %s\n", binfile);
147 fclose(fp);
148 exit(3);
149 }
150 if (fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != nucell){
151 printf("\nError while reading atomic masses from file: %s\n", binfile);
152 fclose(fp);
153 exit(3);
154 }
155 fclose(fp);
156
157 car2dir();
158 real2rec();
159
160 // initialize interpolation
161 interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all);
162 if (flag_reset_gamma) interpolate->reset_gamma();
163
164 // Enforcing Austic Sum Rule
165 EnforceASR();
166
167 // get the dynamical matrix from force constant matrix: D = 1/M x Phi
168 for (int idq = 0; idq < npt; ++idq){
169 int ndim =0;
170 for (int idim = 0; idim < fftdim; ++idim)
171 for (int jdim = 0; jdim < fftdim; ++jdim){
172 double inv_mass = M_inv_sqrt[idim/sysdim] * M_inv_sqrt[jdim/sysdim];
173 DM_all[idq][ndim].r *= inv_mass;
174 DM_all[idq][ndim].i *= inv_mass;
175 ndim++;
176 }
177 }
178
179 // ask for the interpolation method
180 interpolate->set_method();
181
182 return;
183 }
184
185
186 /* ----------------------------------------------------------------------------
187 * Free the memories used.
188 * ---------------------------------------------------------------------------- */
~DynMat()189 DynMat::~DynMat()
190 {
191 // destroy all memory allocated
192 if (funit) delete []funit;
193 if (dmfile) delete []dmfile;
194 if (binfile) delete []binfile;
195 if (interpolate) delete interpolate;
196
197 memory->destroy(DM_q);
198 memory->destroy(attyp);
199 memory->destroy(basis);
200 memory->destroy(DM_all);
201 memory->destroy(M_inv_sqrt);
202 delete memory;
203
204 return;
205 }
206
207 /* ----------------------------------------------------------------------------
208 * method to write DM_q to file, single point
209 * ---------------------------------------------------------------------------- */
writeDMq(double * q)210 void DynMat::writeDMq(double *q)
211 {
212 FILE *fp;
213 // only ask for file name for the first time
214 // other calls will append the result to the file.
215 if (dmfile == NULL){
216 char str[MAXLINE], *ptr;
217 printf("\n");
218 while ( 1 ){
219 printf("Please input the filename to output the DM at selected q: ");
220 fgets(str,MAXLINE,stdin);
221 ptr = strtok(str, " \r\t\n\f");
222 if (ptr) break;
223 }
224
225 int n = strlen(ptr) + 1;
226 dmfile = new char[n];
227 strcpy(dmfile, ptr);
228 fp = fopen(dmfile,"w");
229
230 } else {
231 fp = fopen(dmfile,"a");
232 }
233 fprintf(fp,"# q = [%lg %lg %lg]\n", q[0], q[1], q[2]);
234
235 for (int i = 0; i < fftdim; ++i){
236 for (int j = 0; j < fftdim; ++j) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i);
237 fprintf(fp,"\n");
238 }
239 fprintf(fp,"\n");
240 fclose(fp);
241
242 return;
243 }
244
245 /* ----------------------------------------------------------------------------
246 * method to write DM_q to file, dispersion-like
247 * ---------------------------------------------------------------------------- */
writeDMq(double * q,const double qr,FILE * fp)248 void DynMat::writeDMq(double *q, const double qr, FILE *fp)
249 {
250 fprintf(fp, "%lg %lg %lg %lg ", q[0], q[1], q[2], qr);
251
252 for (int i = 0; i < fftdim; ++i)
253 for (int j = 0; j < fftdim; ++j) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i);
254
255 fprintf(fp,"\n");
256 return;
257 }
258
259 /* ----------------------------------------------------------------------------
260 * method to evaluate the eigenvalues of current q-point;
261 * return the eigenvalues in egv.
262 * cLapack subroutine zheevd is employed.
263 * ---------------------------------------------------------------------------- */
geteigen(double * egv,int flag)264 int DynMat::geteigen(double *egv, int flag)
265 {
266 char jobz, uplo;
267 integer n, lda, lwork, lrwork, *iwork, liwork, info;
268 doublecomplex *work;
269 doublereal *w = &egv[0], *rwork;
270
271 n = fftdim;
272 if (flag) jobz = 'V';
273 else jobz = 'N';
274
275 uplo = 'U';
276 lwork = (n+2)*n;
277 lrwork = 1 + (5+n+n)*n;
278 liwork = 3 + 5*n;
279 lda = n;
280
281 memory->create(work, lwork, "geteigen:work");
282 memory->create(rwork, lrwork, "geteigen:rwork");
283 memory->create(iwork, liwork, "geteigen:iwork");
284
285 zheevd_(&jobz, &uplo, &n, DM_q[0], &lda, w, work, &lwork, rwork, &lrwork, iwork, &liwork, &info);
286
287 // to get w instead of w^2; and convert w into v (THz hopefully)
288 for (int i = 0; i < n; ++i){
289 if (w[i]>= 0.) w[i] = sqrt(w[i]);
290 else w[i] = -sqrt(-w[i]);
291
292 w[i] *= eml2f;
293 }
294
295 memory->destroy(work);
296 memory->destroy(rwork);
297 memory->destroy(iwork);
298
299 return info;
300 }
301
302 /* ----------------------------------------------------------------------------
303 * method to get the Dynamical Matrix at q
304 * ---------------------------------------------------------------------------- */
getDMq(double * q)305 void DynMat::getDMq(double *q)
306 {
307 interpolate->execute(q, DM_q[0]);
308
309 return;
310 }
311
312 /* ----------------------------------------------------------------------------
313 * method to get the Dynamical Matrix at q
314 * ---------------------------------------------------------------------------- */
getDMq(double * q,double * wt)315 void DynMat::getDMq(double *q, double *wt)
316 {
317 interpolate->execute(q, DM_q[0]);
318
319 if (flag_skip && interpolate->UseGamma ) wt[0] = 0.;
320 return;
321 }
322
323 /* ----------------------------------------------------------------------------
324 * private method to convert the cartisan coordinate of basis into fractional
325 * ---------------------------------------------------------------------------- */
car2dir()326 void DynMat::car2dir()
327 {
328 double mat[9];
329 for (int idim = 0; idim < 9; ++idim) mat[idim] = basevec[idim];
330 GaussJordan(3, mat);
331
332 for (int i = 0; i < nucell; ++i){
333 double x[3];
334 x[0] = x[1] = x[2] = 0.;
335 for (int idim = 0; idim < sysdim; idim++) x[idim] = basis[i][idim];
336 for (int idim = 0; idim < sysdim; idim++)
337 basis[i][idim] = x[0]*mat[idim] + x[1]*mat[3+idim] + x[2]*mat[6+idim];
338 }
339
340 return;
341 }
342
343 /* ----------------------------------------------------------------------------
344 * private method to enforce the acoustic sum rule on force constant matrix at G
345 * ---------------------------------------------------------------------------- */
EnforceASR()346 void DynMat::EnforceASR()
347 {
348 char str[MAXLINE];
349 int nasr = 20;
350 if (nucell <= 1) nasr = 1;
351 printf("\n"); for (int i = 0; i < 80; ++i) printf("=");
352
353 // compute and display eigenvalues of Phi at gamma before ASR
354 if (nucell > 100){
355 printf("\nYour unit cell is rather large, eigenvalue evaluation takes some time...");
356 fflush(stdout);
357 }
358
359 double egvs[fftdim];
360 for (int i = 0; i < fftdim; ++i)
361 for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j];
362 geteigen(egvs, 0);
363 printf("\nEigenvalues of Phi at gamma before enforcing ASR:\n");
364 for (int i = 0; i < fftdim; ++i){
365 printf("%lg ", egvs[i]);
366 if (i%10 == 9) printf("\n");
367 if (i == 99){ printf("...... (%d more skipped)\n", fftdim-100); break;}
368 }
369 printf("\n\n");
370
371 // ask for iterations to enforce ASR
372 printf("Please input the # of iterations to enforce ASR [%d]: ", nasr);
373 fgets(str,MAXLINE,stdin);
374 char *ptr = strtok(str," \t\n\r\f");
375 if (ptr) nasr = atoi(ptr);
376 if (nasr < 1){
377 for (int i=0; i<80; i++) printf("="); printf("\n");
378 return;
379 }
380
381 for (int iit = 0; iit < nasr; ++iit){
382 // simple ASR; the resultant matrix might not be symmetric
383 for (int a = 0; a < sysdim; ++a)
384 for (int b = 0; b < sysdim; ++b){
385 for (int k = 0; k < nucell; ++k){
386 double sum = 0.;
387 for (int kp = 0; kp < nucell; ++kp){
388 int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
389 sum += DM_all[0][idx].r;
390 }
391 sum /= double(nucell);
392 for (int kp = 0; kp < nucell; ++kp){
393 int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
394 DM_all[0][idx].r -= sum;
395 }
396 }
397 }
398
399 // symmetrize
400 for (int k = 0; k < nucell; ++k)
401 for (int kp = k; kp < nucell; ++kp){
402 double csum = 0.;
403 for (int a = 0; a < sysdim; ++a)
404 for (int b = 0; b < sysdim; ++b){
405 int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
406 int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a;
407 csum = (DM_all[0][idx].r + DM_all[0][jdx].r )*0.5;
408 DM_all[0][idx].r = DM_all[0][jdx].r = csum;
409 }
410 }
411 }
412
413 // symmetric ASR
414 for (int a = 0; a < sysdim; ++a)
415 for (int b = 0; b < sysdim; ++b){
416 for (int k = 0; k < nucell; ++k){
417 double sum = 0.;
418 for (int kp = 0; kp < nucell; ++kp){
419 int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
420 sum += DM_all[0][idx].r;
421 }
422 sum /= double(nucell-k);
423 for (int kp = k; kp < nucell; ++kp){
424 int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
425 int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a;
426 DM_all[0][idx].r -= sum;
427 DM_all[0][jdx].r = DM_all[0][idx].r;
428 }
429 }
430 }
431
432 // compute and display eigenvalues of Phi at gamma after ASR
433 for (int i = 0; i < fftdim; ++i)
434 for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j];
435 geteigen(egvs, 0);
436 printf("Eigenvalues of Phi at gamma after enforcing ASR:\n");
437 for (int i = 0; i < fftdim; ++i){
438 printf("%lg ", egvs[i]);
439 if (i%10 == 9) printf("\n");
440 if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;}
441 }
442 printf("\n");
443 for (int i = 0; i < 80; ++i) printf("="); printf("\n\n");
444
445 return;
446 }
447
448 /* ----------------------------------------------------------------------------
449 * private method to get the reciprocal lattice vectors from the real space ones
450 * ---------------------------------------------------------------------------- */
real2rec()451 void DynMat::real2rec()
452 {
453 ibasevec[0] = basevec[4]*basevec[8] - basevec[5]*basevec[7];
454 ibasevec[1] = basevec[5]*basevec[6] - basevec[3]*basevec[8];
455 ibasevec[2] = basevec[3]*basevec[7] - basevec[4]*basevec[6];
456
457 ibasevec[3] = basevec[7]*basevec[2] - basevec[8]*basevec[1];
458 ibasevec[4] = basevec[8]*basevec[0] - basevec[6]*basevec[2];
459 ibasevec[5] = basevec[6]*basevec[1] - basevec[7]*basevec[0];
460
461 ibasevec[6] = basevec[1]*basevec[5] - basevec[2]*basevec[4];
462 ibasevec[7] = basevec[2]*basevec[3] - basevec[0]*basevec[5];
463 ibasevec[8] = basevec[0]*basevec[4] - basevec[1]*basevec[3];
464
465 double vol = 0.;
466 for (int i = 0; i < sysdim; ++i) vol += ibasevec[i] * basevec[i];
467 vol = 8.*atan(1.)/vol;
468
469 for (int i = 0; i < 9; ++i) ibasevec[i] *= vol;
470
471 printf("\n"); for (int i = 0; i < 80; ++i) printf("=");
472 printf("\nBasis vectors of the unit cell in real space:");
473 for (int i = 0; i < sysdim; ++i){
474 printf("\n A%d: ", i+1);
475 for (int j = 0; j < sysdim; ++j) printf("%8.4f ", basevec[i*3+j]);
476 }
477 printf("\nBasis vectors of the corresponding reciprocal cell:");
478 for (int i = 0; i < sysdim; ++i){
479 printf("\n B%d: ", i+1);
480 for (int j = 0; j < sysdim; ++j) printf("%8.4f ", ibasevec[i*3+j]);
481 }
482 printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
483
484 return;
485 }
486
487 /* ----------------------------------------------------------------------
488 * private method, to get the inverse of a double matrix by means of
489 * Gaussian-Jordan Elimination with full pivoting; square matrix required.
490 *
491 * Adapted from the Numerical Recipes in Fortran.
492 * --------------------------------------------------------------------*/
GaussJordan(int n,double * Mat)493 void DynMat::GaussJordan(int n, double *Mat)
494 {
495 int i,icol,irow,j,k,l,ll,idr,idc;
496 int *indxc,*indxr,*ipiv;
497 double big, nmjk;
498 double dum, pivinv;
499
500 indxc = new int[n];
501 indxr = new int[n];
502 ipiv = new int[n];
503
504 for (i = 0; i < n; ++i) ipiv[i] = 0;
505 for (i = 0; i < n; ++i){
506 big = 0.;
507 for (j = 0; j < n; ++j){
508 if (ipiv[j] != 1){
509 for (k = 0; k < n; ++k){
510 if (ipiv[k] == 0){
511 idr = j * n + k;
512 nmjk = abs(Mat[idr]);
513 if (nmjk >= big){
514 big = nmjk;
515 irow = j;
516 icol = k;
517 }
518 } else if (ipiv[k] > 1){
519 printf("DynMat: Singular matrix in double GaussJordan!\n"); exit(1);
520 }
521 }
522 }
523 }
524 ipiv[icol] += 1;
525 if (irow != icol){
526 for (l = 0; l < n; ++l){
527 idr = irow*n+l;
528 idc = icol*n+l;
529 dum = Mat[idr];
530 Mat[idr] = Mat[idc];
531 Mat[idc] = dum;
532 }
533 }
534 indxr[i] = irow;
535 indxc[i] = icol;
536 idr = icol * n + icol;
537 if (Mat[idr] == 0.){
538 printf("DynMat: Singular matrix in double GaussJordan!");
539 exit(1);
540 }
541
542 pivinv = 1./ Mat[idr];
543 Mat[idr] = 1.;
544 idr = icol*n;
545 for (l = 0; l < n; ++l) Mat[idr+l] *= pivinv;
546 for (ll = 0; ll < n; ++ll){
547 if (ll != icol){
548 idc = ll * n + icol;
549 dum = Mat[idc];
550 Mat[idc] = 0.;
551 idc -= icol;
552 for (l = 0; l < n; ++l) Mat[idc+l] -= Mat[idr+l]*dum;
553 }
554 }
555 }
556 for (l = n-1; l >= 0; --l){
557 int rl = indxr[l];
558 int cl = indxc[l];
559 if (rl != cl){
560 for (k = 0; k < n; ++k){
561 idr = k * n + rl;
562 idc = k * n + cl;
563 dum = Mat[idr];
564 Mat[idr] = Mat[idc];
565 Mat[idc] = dum;
566 }
567 }
568 }
569 delete []indxr;
570 delete []indxc;
571 delete []ipiv;
572
573 return;
574 }
575
576 /* ----------------------------------------------------------------------------
577 * Public method to reset the interpolation method
578 * ---------------------------------------------------------------------------- */
reset_interp_method()579 void DynMat::reset_interp_method()
580 {
581 interpolate->set_method();
582
583 return;
584 }
585
586 /* ----------------------------------------------------------------------------
587 * Private method to display help info
588 * ---------------------------------------------------------------------------- */
help()589 void DynMat::help()
590 {
591 ShowVersion();
592 printf("\nUsage:\n phana [options] [file]\n\n");
593 printf("Available options:\n");
594 printf(" -r To reset the dynamical matrix at the gamma point by a 4th order\n");
595 printf(" polynomial interpolation along the [100] direction; this might be\n");
596 printf(" useful for systems with charges. As for charged system, the dynamical\n");
597 printf(" matrix at Gamma is far from accurate because of the long range nature\n");
598 printf(" of Coulombic interaction. By reset it by interpolation, will partially\n");
599 printf(" elliminate the unwanted behavior, but the result is still inaccurate.\n");
600 printf(" By default, this is not set; and not expected for uncharged systems.\n\n");
601 printf(" -s This will reset the dynamical matrix at the gamma point, too, but it\n");
602 printf(" will also inform the code to skip all q-points that is in the vicinity\n");
603 printf(" of the gamma point when evaluating phonon DOS and/or phonon dispersion.\n\n");
604 printf(" By default, this is not set; and not expected for uncharged systems.\n\n");
605 printf(" -h To print out this help info.\n\n");
606 printf(" file To define the filename that carries the binary dynamical matrice generated\n");
607 printf(" by fix-phonon. If not provided, the code will ask for it.\n");
608 printf("\n\n");
609 exit(0);
610 }
611
612 /* ----------------------------------------------------------------------------
613 * Private method to display the version info
614 * ---------------------------------------------------------------------------- */
ShowVersion()615 void DynMat::ShowVersion()
616 {
617 printf(" ____ _ _ __ _ _ __ \n");
618 printf(" ( _ \\( )_( ) /__\\ ( \\( ) /__\\ \n");
619 printf(" )___/ ) _ ( /(__)\\ ) ( /(__)\\ \n");
620 printf(" (__) (_) (_)(__)(__)(_)\\_)(__)(__)\n");
621 printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%02d, compiled on %s.\n", VERSION, __DATE__);
622 printf("Reference: https://doi.org/10.1016/j.cpc.2011.04.019\n");
623
624 return;
625 }
626
627 /* ----------------------------------------------------------------------------
628 * Private method to define the conversion factor from the DM measured to THZ
629 * for the phonon frequency and force constants.
630 * ---------------------------------------------------------------------------- */
Define_Conversion_Factor()631 void DynMat::Define_Conversion_Factor()
632 {
633 funit = new char[4];
634 strcpy(funit, "THz");
635
636 if (fabs(boltz - 1.) <= ZERO){ // LJ Unit
637 eml2f = eml2fc = 1.;
638 delete funit;
639 funit = new char[27];
640 strcpy(funit,"sqrt(epsilon/(m.sigma^2))");
641
642 } else if (fabs(boltz - 0.0019872067) <= ZERO){ // real
643 eml2f = 3.255487031;
644 eml2fc = 0.0433641042418;
645
646 } else if (fabs(boltz*1.e3 - 8.617343e-2) <= ZERO){ // metal
647 eml2f = 15.633304237154924;
648 eml2fc = 1.;
649
650 } else if (fabs(boltz*1.e20 - 1.3806504e-3) <= ZERO){ // si
651 eml2f = 1.591549431e-13;
652 eml2fc = 0.06241509074460763;
653
654 } else if (fabs(boltz*1.e13 - 1.3806504e-3) <= ZERO){ // cgs
655 eml2f = 1.591549431e-13;
656 eml2fc = 6.241509074460763e-05;
657
658 } else if (fabs(boltz*1.e3 - 3.16681534e-3) <= ZERO){ // electron
659 eml2f = 154.10792761319672;
660 eml2fc = 97.1736242922823;
661
662 } else if (fabs(boltz*1.e5 - 1.3806504e-3) <= ZERO){ // micro
663 eml2f = 1.5915494309189532e-07;
664 eml2fc = 6.241509074460763e-05;
665
666 } else if (fabs(boltz - 0.013806504) <= ZERO){ // nano
667 eml2f = 0.0001591549431;
668 eml2fc = 6.241509074460763e-05;
669
670 } else {
671 printf("WARNING: Perhaps because of float precision, I cannot get the factor to convert\n");
672 printf("sqrt(E/ML^2)/(2*pi) into THz, instead, I set it to 1; you should check the unit\nused by LAMMPS.\n");
673 eml2f = eml2fc = 1.;
674 }
675 return;
676 }
677
678 /* ----------------------------------------------------------------------------
679 * Private method to output the information read
680 * ---------------------------------------------------------------------------- */
ShowInfo()681 void DynMat::ShowInfo()
682 {
683 printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
684 printf("Dynamical matrix is read from file: %s\n", binfile);
685 printf("The system size in three dimension: %d x %d x %d\n", nx, ny, nz);
686 printf("Number of atoms per unit cell : %d\n", nucell);
687 printf("System dimension : %d\n", sysdim);
688 printf("Boltzmann constant in used units : %g\n", boltz);
689 for (int i = 0; i < 80; ++i) printf("="); printf("\n");
690 return;
691 }
692 /* --------------------------------------------------------------------*/
693