1 /**********************************************************************
2 Make_InputFile_with_FinalCoord.c:
3
4 Make_InputFile_with_FinalCoord.c is a subrutine to make an input file
5 with the final coordinate of the system.
6
7 Log of Make_InputFile_with_FinalCoord.c:
8
9 19/Sep./2007 Released by T. Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <time.h>
16 #include <ctype.h>
17 #include <string.h>
18
19 /* stat section */
20 #include <sys/types.h>
21 #include <sys/stat.h>
22 #include <unistd.h>
23 /* end stat section */
24 #include "openmx_common.h"
25 #include "tran_variables.h"
26 #include "mpi.h"
27
28
29 void Make_InputFile_with_FinalCoord_Normal(char *file, int MD_iter);
30 void Make_InputFile_with_FinalCoord_NEGF(char *file, int MD_iter);
31
32
33
Make_InputFile_with_FinalCoord(char * file,int MD_iter)34 void Make_InputFile_with_FinalCoord(char *file, int MD_iter)
35 {
36
37 if (Solver!=4)
38 Make_InputFile_with_FinalCoord_Normal(file, MD_iter);
39 else
40 Make_InputFile_with_FinalCoord_NEGF(file, MD_iter);
41
42 }
43
44
45
46
Make_InputFile_with_FinalCoord_Normal(char * file,int MD_iter)47 void Make_InputFile_with_FinalCoord_Normal(char *file, int MD_iter)
48 {
49 int i,j,Gc_AN,c,n1,k;
50 int restart_flag,fixed_flag;
51 int geoopt_restart_flag;
52 int velocity_flag;
53 int MD_Current_Iter_flag;
54 int Atoms_UnitVectors_flag;
55 int Atoms_UnitVectors_Unit_flag;
56 int rstfile_num;
57 double c1,c2,c3;
58 double vx,vy,vz;
59 double tmpxyz[4];
60 char st[800];
61 char st1[800];
62 char rm_operate[YOUSO10];
63 char fname[YOUSO10];
64 char fname1[YOUSO10];
65 char fname2[YOUSO10];
66 FILE *fp1,*fp2;
67 char *tp;
68 int numprocs,myid;
69
70 /* MPI */
71 MPI_Comm_size(mpi_comm_level1,&numprocs);
72 MPI_Comm_rank(mpi_comm_level1,&myid);
73
74 if (myid==Host_ID){
75
76 /* initialize */
77
78 restart_flag = 0;
79 geoopt_restart_flag = 0;
80 fixed_flag = 0;
81 velocity_flag = 0;
82 MD_Current_Iter_flag = 0;
83 Atoms_UnitVectors_flag = 0;
84 Atoms_UnitVectors_Unit_flag = 0;
85
86 /* the new input file */
87
88 sprintf(fname1,"%s#",file);
89 fp1 = fopen(fname1,"w");
90 fseek(fp1,0,SEEK_END);
91
92 /* the original input file */
93
94 fp2 = fopen(file,"r");
95
96 if (fp2!=NULL){
97
98 while (fgets(st,800,fp2)!=NULL){
99
100 string_tolower(st,st1);
101
102 /* find the specification of <atoms.speciesandcoordinates */
103
104 if (strncmp(st1,"<atoms.speciesandcoordinates",28)==0){
105
106 fprintf(fp1,"%s",st);
107
108 /* replace the atomic coordinates */
109
110 for (i=1; i<=atomnum; i++){
111
112 fgets(st,800,fp2);
113 string_tolower(st,st1);
114
115 /* serial number */
116 tp = strtok(st, " ");
117 if (tp!=NULL) fprintf(fp1,"%4s",tp);
118
119 /* name of species */
120 tp =strtok(NULL, " ");
121 if (tp!=NULL) fprintf(fp1," %4s",tp);
122
123 /* "Ang" */
124 if (coordinates_unit==0){
125
126 /* x-coordinate */
127 tp =strtok(NULL, " ");
128 fprintf(fp1," %18.14f",Gxyz[i][1]*BohrR);
129
130 /* y-coordinate */
131 tp =strtok(NULL, " ");
132 fprintf(fp1," %18.14f",Gxyz[i][2]*BohrR);
133
134 /* z-coordinate */
135 tp =strtok(NULL, " ");
136 fprintf(fp1," %18.14f",Gxyz[i][3]*BohrR);
137 }
138
139 /* AU */
140 else if (coordinates_unit==1){
141
142 /* x-coordinate */
143 tp =strtok(NULL, " ");
144 fprintf(fp1," %18.14f",Gxyz[i][1]);
145
146 /* y-coordinate */
147 tp =strtok(NULL, " ");
148 fprintf(fp1," %18.14f",Gxyz[i][2]);
149
150 /* z-coordinate */
151 tp =strtok(NULL, " ");
152 fprintf(fp1," %18.14f",Gxyz[i][3]);
153 }
154
155 /* FRAC */
156 else if (coordinates_unit==2){
157
158 /* The zero is taken as the origin of the unit cell. */
159
160 tmpxyz[1] = Gxyz[i][1] - Grid_Origin[1];
161 tmpxyz[2] = Gxyz[i][2] - Grid_Origin[2];
162 tmpxyz[3] = Gxyz[i][3] - Grid_Origin[3];
163
164 c1 = Dot_Product(tmpxyz,rtv[1])*0.5/PI;
165 c2 = Dot_Product(tmpxyz,rtv[2])*0.5/PI;
166 c3 = Dot_Product(tmpxyz,rtv[3])*0.5/PI;
167
168 /* a-axis */
169 tp =strtok(NULL, " ");
170 fprintf(fp1," %18.14f",c1);
171
172 /* b-axis */
173 tp =strtok(NULL, " ");
174 fprintf(fp1," %18.14f",c2);
175
176 /* c-axis */
177 tp =strtok(NULL, " ");
178 fprintf(fp1," %18.14f",c3);
179 }
180
181 while (tp!=NULL){
182 tp =strtok(NULL, " ");
183 if (tp!=NULL) fprintf(fp1," %s",tp);
184 }
185
186 }
187
188 }
189
190 /* scf.restart */
191
192 else if (strncmp(st1,"scf.restart",11)==0){
193 fprintf(fp1,"scf.restart on\n");
194 restart_flag = 1;
195 }
196
197 /* geoopt.restart */
198
199 else if (strncmp(st1,"geoopt.restart",11)==0){
200 fprintf(fp1,"geoopt.restart on\n");
201 geoopt_restart_flag = 1;
202 }
203
204 /* scf.fixed.grid */
205
206 else if (strncmp(st1,"scf.fixed.grid",14)==0){
207 fprintf(fp1,"scf.fixed.grid %18.14f %18.14f %18.14f\n",
208 Grid_Origin[1],Grid_Origin[2],Grid_Origin[3]);
209 fixed_flag = 1;
210 }
211
212 /* MD.Init.Velocity && VerletXYZ (NEV) */
213
214 else if (strncmp(st1,"<md.init.velocity",16)==0
215 && (MD_switch==1 || MD_switch==2 || MD_switch==9 || MD_switch==11
216 || MD_switch==14 || MD_switch==15)){
217
218 for (i=1; i<=(atomnum+1); i++){
219 fgets(st,800,fp2);
220 }
221
222 fprintf(fp1,"<MD.Init.Velocity\n");
223
224 for (i=1; i<=atomnum; i++){
225
226 vx = Gxyz[i][24]/(0.4571028*0.000001);
227 vy = Gxyz[i][25]/(0.4571028*0.000001);
228 vz = Gxyz[i][26]/(0.4571028*0.000001);
229 fprintf(fp1," %4d %20.14f %20.14f %20.14f\n", i, vx, vy, vz);
230 }
231
232 fprintf(fp1,"MD.Init.Velocity>\n");
233
234 velocity_flag = 1;
235 }
236
237 /* MD.Current.Iter */
238
239 else if (strncmp(st1,"md.current.iter",15)==0){
240 fprintf(fp1,"\n\nMD.Current.Iter %2d\n",MD_iter);
241 MD_Current_Iter_flag = 1;
242 }
243
244 /* Atoms.UnitVectors */
245
246 else if (strncmp(st1,"<atoms.unitvectors",18)==0){
247
248 for (i=1; i<=4; i++){
249 fgets(st,800,fp2);
250 }
251
252 fprintf(fp1,"<Atoms.UnitVectors\n");
253
254 for (i=1; i<=3; i++){
255 for (j=1; j<=3; j++){
256 fprintf(fp1," %18.15f ",tv[i][j]);
257 }
258 fprintf(fp1,"\n");
259 }
260
261 fprintf(fp1,"Atoms.UnitVectors>\n");
262
263 Atoms_UnitVectors_flag = 1;
264 }
265
266 /* Atoms.UnitVectors.Unit */
267
268 else if (strncmp(st1,"atoms.unitvectors.unit",22)==0){
269
270 fprintf(fp1,"Atoms.UnitVectors.Unit AU\n");
271 Atoms_UnitVectors_Unit_flag = 1;
272 }
273
274 else{
275 fprintf(fp1,"%s",st);
276 }
277 }
278
279 fclose(fp2);
280 }
281
282 /* add the restart flag if it was not found. */
283
284 if (restart_flag==0){
285 fprintf(fp1,"\n\nscf.restart on\n");
286 }
287
288 /* add the restart flag for geometry optimization if it was not found. */
289
290 if (geoopt_restart_flag==0 &&
291 (MD_switch==3 || MD_switch==4 || MD_switch==5 || MD_switch==6 || MD_switch==7) ){
292 fprintf(fp1,"\n\ngeoopt.restart on\n");
293 }
294
295 /* add scf.fixed.grid if it was not found. */
296
297 if (fixed_flag==0){
298 fprintf(fp1,"\n\nscf.fixed.grid %18.14f %18.14f %18.14f\n",
299 Grid_Origin[1],Grid_Origin[2],Grid_Origin[3]);
300 }
301
302 /* add velocity frag if it was not found. */
303
304 if (velocity_flag==0 && (MD_switch==1 || MD_switch==2 || MD_switch==9 || MD_switch==11
305 || MD_switch==14 || MD_switch==15)){
306
307 fprintf(fp1,"\n\n<MD.Init.Velocity\n");
308
309 for (i=1; i<=atomnum; i++){
310 fprintf(fp1," %4d %20.14f %20.14f %20.14f\n",
311 i,
312 Gxyz[i][24]/(0.4571028*0.000001),
313 Gxyz[i][25]/(0.4571028*0.000001),
314 Gxyz[i][26]/(0.4571028*0.000001));
315 }
316
317 fprintf(fp1,"MD.Init.Velocity>\n");
318
319 }
320
321 /* add MD.Current.Iter if it was not found. */
322
323 if (MD_Current_Iter_flag==0){
324 fprintf(fp1,"\n\nMD.Current.Iter %2d\n",MD_iter);
325 }
326
327 /* add Atoms.UnitVectors if it was not found. */
328
329 if (Atoms_UnitVectors_flag==0){
330
331 fprintf(fp1,"\n\n");
332 fprintf(fp1,"<Atoms.UnitVectors\n");
333
334 for (i=1; i<=3; i++){
335 for (j=1; j<=3; j++){
336 fprintf(fp1," %18.15f ",tv[i][j]);
337 }
338 fprintf(fp1,"\n");
339 }
340
341 fprintf(fp1,"Atoms.UnitVectors>\n");
342 }
343
344 /* add Atoms.UnitVectors if it was not found. */
345
346 if (Atoms_UnitVectors_Unit_flag==0){
347
348 fprintf(fp1,"\n\n");
349 fprintf(fp1,"Atoms.UnitVectors.Unit AU\n");
350 }
351
352 /* fclose */
353 fclose(fp1);
354
355 } /* if (myid==Host_ID) */
356
357 }
358
359
360
Make_InputFile_with_FinalCoord_NEGF(char * file,int MD_iter)361 void Make_InputFile_with_FinalCoord_NEGF(char *file, int MD_iter)
362 {
363 int i,Gc_AN,c,n1,k;
364 int restart_flag,fixed_flag;
365 int velocity_flag;
366 int rstfile_num;
367 double c1,c2,c3;
368 char st[800];
369 char st1[800];
370 char rm_operate[YOUSO10];
371 char fname[YOUSO10];
372 char fname1[YOUSO10];
373 char fname2[YOUSO10];
374 FILE *fp1,*fp2;
375 char *tp;
376 int numprocs,myid;
377
378 /* MPI */
379 MPI_Comm_size(mpi_comm_level1,&numprocs);
380 MPI_Comm_rank(mpi_comm_level1,&myid);
381
382 if (myid==Host_ID){
383
384 /* initialize */
385
386 restart_flag = 0;
387 fixed_flag = 0;
388 velocity_flag = 0;
389
390 /* the new input file */
391
392 sprintf(fname1,"%s#",file);
393 fp1 = fopen(fname1,"w");
394 fseek(fp1,0,SEEK_END);
395
396 /* the original input file */
397
398 fp2 = fopen(file,"r");
399
400 if (fp2!=NULL){
401
402 while (fgets(st,800,fp2)!=NULL){
403
404 string_tolower(st,st1);
405
406 /* find the specification of <leftleadatoms.speciesandcoordinates */
407
408 if (strncmp(st1,"<leftleadatoms.speciesandcoordinates",35)==0){
409
410 fprintf(fp1,"%s",st);
411
412 /* replace the atomic coordinates */
413
414 for (i=1; i<=Latomnum; i++){
415
416 fgets(st,800,fp2);
417 string_tolower(st,st1);
418
419 /* serial number */
420 tp = strtok(st, " ");
421 if (tp!=NULL) fprintf(fp1,"%4s",tp);
422
423 /* name of species */
424 tp =strtok(NULL, " ");
425 if (tp!=NULL) fprintf(fp1," %4s",tp);
426
427 /* "Ang" */
428 if (coordinates_unit==0){
429
430 /* x-coordinate */
431 tp =strtok(NULL, " ");
432 fprintf(fp1," %18.14f",Gxyz[i][1]*BohrR);
433
434 /* y-coordinate */
435 tp =strtok(NULL, " ");
436 fprintf(fp1," %18.14f",Gxyz[i][2]*BohrR);
437
438 /* z-coordinate */
439 tp =strtok(NULL, " ");
440 fprintf(fp1," %18.14f",Gxyz[i][3]*BohrR);
441 }
442
443 /* AU */
444 else if (coordinates_unit==1){
445
446 /* x-coordinate */
447 tp =strtok(NULL, " ");
448 fprintf(fp1," %18.14f",Gxyz[i][1]);
449
450 /* y-coordinate */
451 tp =strtok(NULL, " ");
452 fprintf(fp1," %18.14f",Gxyz[i][2]);
453
454 /* z-coordinate */
455 tp =strtok(NULL, " ");
456 fprintf(fp1," %18.14f",Gxyz[i][3]);
457 }
458
459 while (tp!=NULL){
460 tp =strtok(NULL, " ");
461 if (tp!=NULL) fprintf(fp1," %s",tp);
462 }
463
464 }
465 }
466
467 /* find the specification of <atoms.speciesandcoordinates */
468
469 else if (strncmp(st1,"<atoms.speciesandcoordinates",28)==0){
470
471 fprintf(fp1,"%s",st);
472
473 /* replace the atomic coordinates */
474
475 for (i=1; i<=Catomnum; i++){
476
477 fgets(st,800,fp2);
478 string_tolower(st,st1);
479
480 /* serial number */
481 tp = strtok(st, " ");
482 if (tp!=NULL) fprintf(fp1,"%4s",tp);
483
484 /* name of species */
485 tp =strtok(NULL, " ");
486 if (tp!=NULL) fprintf(fp1," %4s",tp);
487
488 /* "Ang" */
489 if (coordinates_unit==0){
490
491 /* x-coordinate */
492 tp =strtok(NULL, " ");
493 fprintf(fp1," %18.14f",Gxyz[i+Latomnum][1]*BohrR);
494
495 /* y-coordinate */
496 tp =strtok(NULL, " ");
497 fprintf(fp1," %18.14f",Gxyz[i+Latomnum][2]*BohrR);
498
499 /* z-coordinate */
500 tp =strtok(NULL, " ");
501 fprintf(fp1," %18.14f",Gxyz[i+Latomnum][3]*BohrR);
502 }
503
504 /* AU */
505 else if (coordinates_unit==1){
506
507 /* x-coordinate */
508 tp =strtok(NULL, " ");
509 fprintf(fp1," %18.14f",Gxyz[i+Latomnum][1]);
510
511 /* y-coordinate */
512 tp =strtok(NULL, " ");
513 fprintf(fp1," %18.14f",Gxyz[i+Latomnum][2]);
514
515 /* z-coordinate */
516 tp =strtok(NULL, " ");
517 fprintf(fp1," %18.14f",Gxyz[i+Latomnum][3]);
518 }
519
520 while (tp!=NULL){
521 tp =strtok(NULL, " ");
522 if (tp!=NULL) fprintf(fp1," %s",tp);
523 }
524 }
525 }
526
527 /* find the specification of <rightleadatoms.speciesandcoordinates */
528
529 else if (strncmp(st1,"<rightleadatoms.speciesandcoordinates",36)==0){
530
531 fprintf(fp1,"%s",st);
532
533 /* replace the atomic coordinates */
534
535 for (i=1; i<=Ratomnum; i++){
536
537 fgets(st,800,fp2);
538 string_tolower(st,st1);
539
540 /* serial number */
541 tp = strtok(st, " ");
542 if (tp!=NULL) fprintf(fp1,"%4s",tp);
543
544 /* name of species */
545 tp =strtok(NULL, " ");
546 if (tp!=NULL) fprintf(fp1," %4s",tp);
547
548 /* "Ang" */
549 if (coordinates_unit==0){
550
551 /* x-coordinate */
552 tp =strtok(NULL, " ");
553 fprintf(fp1," %18.14f",Gxyz[i+Latomnum+Catomnum][1]*BohrR);
554
555 /* y-coordinate */
556 tp =strtok(NULL, " ");
557 fprintf(fp1," %18.14f",Gxyz[i+Latomnum+Catomnum][2]*BohrR);
558
559 /* z-coordinate */
560 tp =strtok(NULL, " ");
561 fprintf(fp1," %18.14f",Gxyz[i+Latomnum+Catomnum][3]*BohrR);
562 }
563
564 /* AU */
565 else if (coordinates_unit==1){
566
567 /* x-coordinate */
568 tp =strtok(NULL, " ");
569 fprintf(fp1," %18.14f",Gxyz[i+Latomnum+Catomnum][1]);
570
571 /* y-coordinate */
572 tp =strtok(NULL, " ");
573 fprintf(fp1," %18.14f",Gxyz[i+Latomnum+Catomnum][2]);
574
575 /* z-coordinate */
576 tp =strtok(NULL, " ");
577 fprintf(fp1," %18.14f",Gxyz[i+Latomnum+Catomnum][3]);
578 }
579
580 while (tp!=NULL){
581 tp =strtok(NULL, " ");
582 if (tp!=NULL) fprintf(fp1," %s",tp);
583 }
584
585 }
586
587 }
588
589 /* scf.restart */
590
591 else if (strncmp(st1,"scf.restart",11)==0){
592 fprintf(fp1,"scf.restart on\n");
593 restart_flag = 1;
594 }
595
596 /* scf.fixed.grid */
597
598 else if (strncmp(st1,"scf.fixed.grid",14)==0){
599
600 fprintf(fp1,"scf.fixed.grid %18.14f %18.14f %18.14f\n",
601 Grid_Origin[1],Grid_Origin[2],Grid_Origin[3]);
602 fixed_flag = 1;
603 }
604
605 /* MD.Init.Velocity && VerletXYZ (NEV) */
606
607 else if (strncmp(st1,"<md.init.velocity",16)==0
608 && (MD_switch==1 || MD_switch==2 || MD_switch==9 || MD_switch==11
609 || MD_switch==14 || MD_switch==15)){
610
611 for (i=1; i<=(atomnum+1); i++){
612 fgets(st,800,fp2);
613 }
614
615 fprintf(fp1,"<MD.Init.Velocity\n");
616
617 for (i=1; i<=atomnum; i++){
618 fprintf(fp1," %4d %20.14f %20.14f %20.14f\n",
619 i,
620 Gxyz[i][24]/(0.4571028*0.000001),
621 Gxyz[i][25]/(0.4571028*0.000001),
622 Gxyz[i][26]/(0.4571028*0.000001));
623 }
624
625 fprintf(fp1,"MD.Init.Velocity>\n");
626
627 velocity_flag = 1;
628 }
629
630 else{
631 fprintf(fp1,"%s",st);
632 }
633 }
634
635 fclose(fp2);
636 }
637
638 /* add the restart flag if it was not found. */
639
640 if (restart_flag==0){
641 fprintf(fp1,"\n\nscf.restart on\n");
642 }
643
644 /* add scf.fixed.grid if it was not found. */
645
646 if (fixed_flag==0){
647 fprintf(fp1,"\n\nscf.fixed.grid %18.14f %18.14f %18.14f\n",
648 Grid_Origin[1],Grid_Origin[2],Grid_Origin[3]);
649 }
650
651 /* add velocity frag if it was not found. */
652
653 if (velocity_flag==0 && (MD_switch==1 || MD_switch==2 || MD_switch==9 || MD_switch==11)){
654
655 fprintf(fp1,"\n\n<MD.Init.Velocity\n");
656
657 for (i=1; i<=atomnum; i++){
658 fprintf(fp1," %4d %20.14f %20.14f %20.14f\n",
659 i,
660 Gxyz[i][24]/(0.4571028*0.000001),
661 Gxyz[i][25]/(0.4571028*0.000001),
662 Gxyz[i][26]/(0.4571028*0.000001));
663 }
664
665 fprintf(fp1,"MD.Init.Velocity>\n");
666 }
667
668 /* fclose */
669 fclose(fp1);
670
671 } /* if (myid==Host_ID) */
672
673 }
674
675
676
677
678
679
680
681
682
683
684