1 /**********************************************************************
2 neb.c:
3
4 neb is a program to perform the nudged elastic band (NEB) method
5 for finding a minimum energy path (MEP) connecting two structures,
6 which is based on JCP 113, 9978 (2000).
7
8 Log of neb.c:
9
10 Apr./06/2011 Released by T. Ozaki
11 Feb./01/2013 Modified by Y. Kubota (supervised by Prof. F. Ishii)
12 Feb./17/2015 Modified by T. Ozaki
13
14 ***********************************************************************/
15
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <math.h>
19 #include <string.h>
20 #include <ctype.h>
21 #include <time.h>
22 #include <sys/types.h>
23 #include <sys/stat.h>
24 #include <unistd.h>
25 #include "openmx_common.h"
26 #include "Inputtools.h"
27 #include "lapack_prototypes.h"
28
29 #define MAXBUF 1024
30 #define Criterion_Max_Step 0.10
31
32 #ifdef MAX
33 #undef MAX
34 #endif
35 #define MAX(a,b) (((a)>(b))? (a):(b))
36
37 #ifdef MIN
38 #undef MIN
39 #endif
40 #define MIN(a,b) (((a)<(b))? (a):(b))
41
42
43 static void allocate_arrays();
44 static void free_arrays();
45 static void read_input(char *file);
46 static void generate_input_files(char *file, int iter);
47 static void Calc_NEB_Gradients(double ***neb_atom_coordinates);
48 static void Update_Coordinates(int iter, double ***neb_atom_coordinates);
49 static void Make_XYZ_File(char fname[YOUSO10], double ***neb_atom_coordinates);
50 static void Generate_Restart_File(char fname[YOUSO10], double ***neb_atom_coordinates);
51 static void Steepest_Decent(int iter, double ***neb_atom_coordinates, double norm, double Max_Force);
52 static double DIIS_BFGS(int iter, double ***neb_atom_coordinates, double norm, double Max_Force);
53
54 double **All_Grid_Origin;
55 double **Tmp_Grid_Origin;
56 double ***neb_atom_coordinates;
57 double ***tmp_neb_atom_coordinates;
58 double ****His_neb_atom_coordinates;
59 double ******InvHess;
60 double ***Vec0;
61 int **atomFixedXYZ;
62
63
64 static char fileXYZ[YOUSO10] = ".neb.xyz";
65 char system_name[YOUSO10];
66 int *WhatSpecies_NEB;
67 int *Spe_WhatAtom_NEB;
68 char **SpeName_NEB;
69
70 int PN;
71
72
73
74
75
neb(int argc,char * argv[])76 void neb(int argc, char *argv[])
77 {
78 int iter,index_images;
79 int po,i,j,k,p,h;
80 int myid,numprocs;
81 int myid1,numprocs1;
82 int myid2,numprocs2;
83 int myid3,numprocs3;
84 int parallel1_flag;
85 int myworld1,ID0,ID1;
86 int Num_Comm_World1;
87 int *NPROCS_ID1;
88 int *Comm_World1;
89 int *NPROCS_WD1;
90 int *Comm_World_StartID1;
91 int parallel2_flag;
92 int myworld2,myworld3;
93 int Num_Comm_World2,Num_Comm_World3;
94 int *NPROCS_ID2;
95 int *NPROCS_ID3;
96 int *Comm_World2;
97 int *Comm_World3;
98 int *NPROCS_WD2;
99 int *NPROCS_WD3;
100 int *Comm_World_StartID2;
101 int *Comm_WOrld_StartID3;
102 char fname_original[YOUSO10];
103 char fname1[YOUSO10];
104 MPI_Comm *MPI_CommWD1;
105 MPI_Comm *MPI_CommWD2;
106 MPI_Comm *MPI_CommWD3;
107 char file_neb_utot[YOUSO10];
108 double TStime,TEtime,a0time,a1time,f0time;
109 double b0time,b1time,c0time,c1time,flatstime,flatetime,sumtime=0.0;
110 FILE *fp;
111
112 dtime(&TStime);
113
114 MPI_Comm_size(MPI_COMM_WORLD1,&numprocs);
115 MPI_Comm_rank(MPI_COMM_WORLD1,&myid);
116
117 /* check argv */
118
119 if (argc==1){
120 printf("\nCould not find an input file.\n\n");
121 MPI_Finalize();
122 exit(0);
123 }
124
125 /****************************************************
126 show a message by the NEB code
127 ****************************************************/
128
129 if (myid==Host_ID){
130 printf("\n*******************************************************\n");
131 printf("*******************************************************\n");
132 printf(" Welcome to the NEB extension of OpenMX \n");
133 printf(" Copyright (C), 2002-2011, T. Ozaki \n");
134 printf(" OpenMX comes with ABSOLUTELY NO WARRANTY. \n");
135 printf(" This is free software, and you are welcome to \n");
136 printf(" redistribute it under the constitution of the GNU-GPL.\n");
137 printf("*******************************************************\n");
138 printf("*******************************************************\n\n");
139 }
140
141 /****************************************************
142 read the input file
143 ****************************************************/
144
145 read_input(argv[1]);
146 fnjoint(filepath,filename,fileXYZ);
147 sprintf(fname_original,"%s",argv[1]);
148
149 /****************************************************
150 compare PN with numprocs and NEB_Num_Images
151 ****************************************************/
152
153 if ((PN > numprocs)||(PN > NEB_Num_Images)){
154
155 PN=0;
156 if (myid==Host_ID){
157 printf("The keyword MD.NEB.Parallel.Number will be ignored if the following conditions are satisfied:\n");
158 printf("MD.NEB.Parallel.Number is equal to or smaller than the number of MPI processes.\n");
159 printf("MD.NEB.Parallel.Number is equal to or smaller than the number of MD.NEB.Number.Images.\n");
160 }
161 }
162
163 /*******************************************************************************
164 If MD.NEB.Parallel.Number is not specifed,
165 *******************************************************************************/
166
167 if (PN==0){
168
169 /****************************************************
170 Two level parallelization is performed.
171 Outer loop: images
172 Inner loop: calculation in each image
173 ****************************************************/
174
175 if ( NEB_Num_Images<=numprocs || numprocs==1 ){
176
177 /****************************************************
178 allocate processes to each image in the World1
179 ****************************************************/
180
181 Num_Comm_World1 = NEB_Num_Images;
182
183 if ( Num_Comm_World1<=numprocs ){
184
185 parallel1_flag = 1;
186
187 NPROCS_ID1 = (int*)malloc(sizeof(int)*numprocs);
188 Comm_World1 = (int*)malloc(sizeof(int)*numprocs);
189 NPROCS_WD1 = (int*)malloc(sizeof(int)*Num_Comm_World1);
190 Comm_World_StartID1 = (int*)malloc(sizeof(int)*Num_Comm_World1);
191 MPI_CommWD1 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World1);
192
193 Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World1, &myworld1, MPI_CommWD1,
194 NPROCS_ID1, Comm_World1, NPROCS_WD1, Comm_World_StartID1);
195
196 MPI_Comm_size(MPI_CommWD1[myworld1],&numprocs1);
197 MPI_Comm_rank(MPI_CommWD1[myworld1],&myid1);
198
199 }
200 else {
201 parallel1_flag = 0;
202 myworld1 = 0;
203 }
204
205 /****************************************************
206 allocate processes to each image in the World2
207 ****************************************************/
208
209 Num_Comm_World2 = 2;
210
211 if ( Num_Comm_World2<=numprocs ){
212
213 parallel2_flag = 1;
214
215 NPROCS_ID2 = (int*)malloc(sizeof(int)*numprocs);
216 Comm_World2 = (int*)malloc(sizeof(int)*numprocs);
217 NPROCS_WD2 = (int*)malloc(sizeof(int)*Num_Comm_World2);
218 Comm_World_StartID2 = (int*)malloc(sizeof(int)*Num_Comm_World2);
219 MPI_CommWD2 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World2);
220
221 Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World2, &myworld2, MPI_CommWD2,
222 NPROCS_ID2, Comm_World2, NPROCS_WD2, Comm_World_StartID2);
223
224 MPI_Comm_size(MPI_CommWD2[myworld2],&numprocs2);
225 MPI_Comm_rank(MPI_CommWD2[myworld2],&myid2);
226 }
227 else {
228 parallel2_flag = 0;
229 myworld2 = 0;
230 }
231
232 /****************************************************
233 SCF calculations for the two terminal structures
234 ****************************************************/
235
236 /* check whether the restart is performed or not */
237
238 sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name);
239
240 if (Scf_RestartFromFile==1 && (fp = fopen(file_neb_utot,"rb"))!=NULL){
241
242 if (myid==Host_ID){
243 fread(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp);
244 fread(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp);
245 }
246
247 MPI_Bcast(&neb_atom_coordinates[0][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1);
248 MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1);
249
250 fclose(fp);
251 }
252
253 /* if the restart is not performed, perform the SCF calulations for the two terminals */
254
255 else {
256
257 /* generate input files */
258
259 generate_input_files(fname_original,1);
260
261 /* In case of parallel2_flag==1 */
262
263 if (parallel2_flag==1){
264
265 if (myworld2==0)
266 index_images = 0;
267 else
268 index_images = NEB_Num_Images + 1;
269
270 sprintf(fname1,"%s_%i",fname_original,index_images);
271 argv[1] = fname1;
272
273 neb_run( argv,MPI_CommWD2[myworld2],index_images,neb_atom_coordinates,
274 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
275
276 /* find a representative ID in the top world for ID=0 in the world2 */
277
278 i = 0; po = 0;
279 do {
280
281 if (Comm_World2[i]==0){
282 ID0 = i;
283 po = 1;
284 }
285 i++;
286 } while (po==0);
287
288 /* find a representative ID in the top world for ID=1 in the world2 */
289
290 i = 0; po = 0;
291 do {
292
293 if (Comm_World2[i]==1){
294 ID1 = i;
295 po = 1;
296 }
297 i++;
298 } while (po==0);
299
300 /* MPI broadcast of neb_atom_coordinates[0] and neb_atom_coordinates[NEB_Num_Images+1] */
301
302 MPI_Barrier(MPI_COMM_WORLD);
303 for (i=0; i<=atomnum; i++){
304 MPI_Bcast(&neb_atom_coordinates[0][i][0], 20, MPI_DOUBLE, ID0, MPI_COMM_WORLD);
305 MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][i][0], 20, MPI_DOUBLE, ID1, MPI_COMM_WORLD);
306 }
307
308 } /* if (parallel2_flag==1) */
309
310 /* In case of parallel2_flag==0 */
311
312 else {
313
314 /* SCF calculation of a terminal with the index of 0 */
315 index_images = 0;
316 sprintf(fname1,"%s_%i",fname_original,index_images);
317 argv[1] = fname1;
318 neb_run( argv,MPI_COMM_WORLD1,index_images,neb_atom_coordinates,
319 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
320
321 /* SCF calculation of a terminal with the index of (NEB_Num_Images + 1) */
322 index_images = NEB_Num_Images + 1;
323 sprintf(fname1,"%s_%i",fname_original,index_images);
324 argv[1] = fname1;
325 neb_run( argv,MPI_COMM_WORLD1,index_images,neb_atom_coordinates,
326 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
327 }
328
329 /* save *_0_rst/*.neb.utot in binary mode */
330
331 if (myid==Host_ID){
332
333 sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name);
334
335 if ((fp = fopen(file_neb_utot,"wb")) != NULL){
336
337 fwrite(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp);
338 fwrite(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp);
339 fclose(fp);
340 }
341 else{
342 printf("Could not open a file %s in neb\n",file_neb_utot);
343 }
344 }
345
346 } /* else */
347
348 /****************************************************
349 optimiziation for finding a minimum energy path
350 connecting the two terminal structures.
351 ****************************************************/
352
353 iter = 1;
354 MD_Opt_OK = 0;
355 dtime(&f0time);
356 do {
357
358 /* if iter==1, generate an input file for restarting */
359
360 if (iter==1) Generate_Restart_File(fname_original,neb_atom_coordinates);
361
362 /* generate input files */
363
364 generate_input_files(fname_original,iter);
365
366 /* In case of parallel1_mode==1 */
367
368 if (parallel1_flag==1){
369
370 sprintf(fname1,"%s_%i",fname_original,myworld1+1);
371 argv[1] = fname1;
372 dtime(&flatstime);
373 neb_run( argv,MPI_CommWD1[myworld1],myworld1+1,neb_atom_coordinates,
374 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
375 dtime(&flatetime);
376 sumtime += (flatetime-flatstime);
377 /* MPI: All_Grid_Origin */
378
379 for (i=0; i<=(NEB_Num_Images+1); i++){
380 for (j=0; j<=3; j++) Tmp_Grid_Origin[i][j] = 0.0;
381 }
382
383 if (myid1==Host_ID){
384 Tmp_Grid_Origin[myworld1+1][1] = Grid_Origin[1];
385 Tmp_Grid_Origin[myworld1+1][2] = Grid_Origin[2];
386 Tmp_Grid_Origin[myworld1+1][3] = Grid_Origin[3];
387 }
388
389 MPI_Barrier(MPI_COMM_WORLD);
390
391 for (p=1; p<=NEB_Num_Images; p++){
392 MPI_Allreduce( &Tmp_Grid_Origin[p][0], &All_Grid_Origin[p][0],
393 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
394 }
395
396 /* MPI: neb_atom_coordinates */
397
398 for (p=1; p<=NEB_Num_Images; p++){
399 for (i=0; i<=atomnum; i++){
400 for (j=0; j<20; j++){
401 tmp_neb_atom_coordinates[p][i][j] = 0.0;
402 }
403 }
404 }
405
406 if (myid1==Host_ID){
407 for (i=0; i<=atomnum; i++){
408 for (j=0; j<20; j++){
409 tmp_neb_atom_coordinates[myworld1+1][i][j] = neb_atom_coordinates[myworld1+1][i][j];;
410 }
411 }
412 }
413
414 MPI_Barrier(MPI_COMM_WORLD);
415
416 for (p=1; p<=NEB_Num_Images; p++){
417 for (i=0; i<=atomnum; i++){
418
419 MPI_Allreduce( &tmp_neb_atom_coordinates[p][i][0], &neb_atom_coordinates[p][i][0],
420 20, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
421
422 }
423 }
424
425 } /* if (parallel1_flag==1) */
426
427 /* In case of parallel1_flag==0 */
428 else {
429
430 for (p=1; p<=NEB_Num_Images; p++){
431
432 sprintf(fname1,"%s_%i",fname_original,p);
433 argv[1] = fname1;
434 neb_run( argv,MPI_COMM_WORLD1,p,neb_atom_coordinates,
435 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
436
437 /* store Grid_Origin */
438 All_Grid_Origin[p][1] = Grid_Origin[1];
439 All_Grid_Origin[p][2] = Grid_Origin[2];
440 All_Grid_Origin[p][3] = Grid_Origin[3];
441 }
442 }
443
444 MPI_Barrier(MPI_COMM_WORLD);
445
446 /* calculate the gradients defined by the NEB method */
447
448 Calc_NEB_Gradients(neb_atom_coordinates);
449
450 /* make a xyz file storing a set of structures of images at the current step */
451
452 Make_XYZ_File(fileXYZ,neb_atom_coordinates);
453
454 /* update the atomic structures of the images */
455
456 Update_Coordinates(iter,neb_atom_coordinates);
457
458 /* generate an input file for restarting */
459
460 Generate_Restart_File(fname_original,neb_atom_coordinates);
461
462 /* increment of iter */
463
464 iter++;
465
466 } while (MD_Opt_OK==0 && iter<=MD_IterNumber);
467 MPI_Barrier(MPI_COMM_WORLD);
468
469 /* freeing of arrays for the World1 */
470
471 if ( Num_Comm_World1<=numprocs ){
472
473 MPI_Comm_free(&MPI_CommWD1[myworld1]);
474 free(MPI_CommWD1);
475 free(Comm_World_StartID1);
476 free(NPROCS_WD1);
477 free(Comm_World1);
478 free(NPROCS_ID1);
479 }
480
481 /* freeing of arrays for the World2 */
482
483 if ( Num_Comm_World2<=numprocs ){
484
485 MPI_Comm_free(&MPI_CommWD2[myworld2]);
486 free(MPI_CommWD2);
487 free(Comm_World_StartID2);
488 free(NPROCS_WD2);
489 free(Comm_World2);
490 free(NPROCS_ID2);
491 }
492
493 /* freeing of arrays */
494
495 free_arrays();
496
497 /* show a final message */
498
499
500 dtime(&TEtime);
501
502 printf("\nThe calculation was normally finished. (proc=%3d) TIME=%lf (s) flat time=%lf \n",
503 myid,(TEtime-TStime),sumtime/(iter-1));
504
505 MPI_Finalize();
506 exit(0);
507
508 } /* if ( NEB_Num_Images<=numprocs || numprocs==1 ) */
509
510 /****************************************************
511 One level parallelization is performed.
512 Outer loop: images
513 ****************************************************/
514
515 else{
516
517 dtime(&a0time);
518 int syou,amari,g,bb;
519
520 syou = NEB_Num_Images/numprocs;
521 amari = NEB_Num_Images%numprocs;
522
523 if(amari==0){
524 g=syou;
525 }
526 else{
527 g=syou+1;
528 }
529
530 Num_Comm_World1=numprocs;
531
532 parallel1_flag = 1;
533
534 NPROCS_ID1 = (int*)malloc(sizeof(int)*numprocs);
535 Comm_World1 = (int*)malloc(sizeof(int)*numprocs);
536 NPROCS_WD1 = (int*)malloc(sizeof(int)*Num_Comm_World1);
537 Comm_World_StartID1 = (int*)malloc(sizeof(int)*Num_Comm_World1);
538 MPI_CommWD1 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World1);
539
540 Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World1, &myworld1, MPI_CommWD1,
541 NPROCS_ID1, Comm_World1, NPROCS_WD1, Comm_World_StartID1);
542
543 MPI_Comm_size(MPI_CommWD1[myworld1],&numprocs1);
544 MPI_Comm_rank(MPI_CommWD1[myworld1],&myid1);
545
546 /****************************************************
547 allocate processes to each image in the World2
548 ****************************************************/
549
550 Num_Comm_World2 = 2;
551
552 if ( Num_Comm_World2<=numprocs ){
553
554 parallel2_flag = 1;
555
556 NPROCS_ID2 = (int*)malloc(sizeof(int)*numprocs);
557 Comm_World2 = (int*)malloc(sizeof(int)*numprocs);
558 NPROCS_WD2 = (int*)malloc(sizeof(int)*Num_Comm_World2);
559 Comm_World_StartID2 = (int*)malloc(sizeof(int)*Num_Comm_World2);
560 MPI_CommWD2 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World2);
561
562 Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World2, &myworld2, MPI_CommWD2,
563 NPROCS_ID2, Comm_World2, NPROCS_WD2, Comm_World_StartID2);
564
565 MPI_Comm_size(MPI_CommWD2[myworld2],&numprocs2);
566 MPI_Comm_rank(MPI_CommWD2[myworld2],&myid2);
567 }
568
569 /****************************************************
570 SCF calculations for the two terminal structures
571 ****************************************************/
572
573 /* check whether the restart is performed or not */
574
575 sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name);
576
577 if (Scf_RestartFromFile==1 && (fp = fopen(file_neb_utot,"rb"))!=NULL){
578
579 if (myid==Host_ID){
580 fread(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp);
581 fread(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp);
582 }
583
584 MPI_Bcast(&neb_atom_coordinates[0][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1);
585 MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1);
586
587 fclose(fp);
588 }
589
590 /* if the restart is not performed, perform the SCF calulations for the two terminals */
591
592 else {
593
594 /* generate input files */
595
596 generate_input_files(fname_original,1);
597
598 /* In case of parallel2_flag==1 */
599
600 if (parallel2_flag==1){
601
602 if (myworld2==0)
603 index_images = 0;
604 else
605 index_images = NEB_Num_Images + 1;
606
607 sprintf(fname1,"%s_%i",fname_original,index_images);
608 argv[1] = fname1;
609
610 neb_run( argv,MPI_CommWD2[myworld2],index_images,neb_atom_coordinates,
611 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
612
613 /* find a representative ID in the top world for ID=0 in the world2 */
614
615 i = 0; po = 0;
616 do {
617
618 if (Comm_World2[i]==0){
619 ID0 = i;
620 po = 1;
621 }
622 i++;
623 } while (po==0);
624
625 /* find a representative ID in the top world for ID=1 in the world2 */
626
627 i = 0; po = 0;
628 do {
629
630 if (Comm_World2[i]==1){
631 ID1 = i;
632 po = 1;
633 }
634 i++;
635 } while (po==0);
636
637 /* MPI broadcast of neb_atom_coordinates[0] and neb_atom_coordinates[NEB_Num_Images+1] */
638
639 MPI_Barrier(MPI_COMM_WORLD);
640 for (i=0; i<=atomnum; i++){
641 MPI_Bcast(&neb_atom_coordinates[0][i][0], 20, MPI_DOUBLE, ID0, MPI_COMM_WORLD);
642 MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][i][0], 20, MPI_DOUBLE, ID1, MPI_COMM_WORLD);
643 }
644
645 } /* if (parallel2_flag==1) */
646
647 /* save *_0_rst/*.neb.utot in binary mode */
648
649 if (myid==Host_ID){
650
651 sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name);
652
653 if ((fp = fopen(file_neb_utot,"wb")) != NULL){
654
655 fwrite(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp);
656 fwrite(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp);
657 fclose(fp);
658 }
659 else{
660 printf("Could not open a file %s in neb\n",file_neb_utot);
661 }
662 }
663
664 } /* else */
665
666 /* amari world */
667
668 if(amari!=0){
669 Num_Comm_World3=amari;
670
671 NPROCS_ID3 = (int*)malloc(sizeof(int)*numprocs);
672 Comm_World3 = (int*)malloc(sizeof(int)*numprocs);
673 NPROCS_WD3 = (int*)malloc(sizeof(int)*Num_Comm_World3);
674 Comm_WOrld_StartID3 = (int*)malloc(sizeof(int)*Num_Comm_World3);
675 MPI_CommWD3 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World3);
676
677 Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World3, &myworld3, MPI_CommWD3,
678 NPROCS_ID3, Comm_World3, NPROCS_WD3, Comm_WOrld_StartID3);
679
680 MPI_Comm_size(MPI_CommWD3[myworld3],&numprocs3);
681 MPI_Comm_rank(MPI_CommWD3[myworld3],&myid3);
682 }
683 /* amari world */
684
685 dtime(&a1time);
686
687 /****************************************************
688 optimization for finding a minimum energy path
689 connecting the two terminal structures.
690 ****************************************************/
691
692 iter = 1;
693 MD_Opt_OK = 0;
694 dtime(&b0time);
695
696 do {
697
698 /* if iter==1, generate an input file for restarting */
699
700 if (iter==1) Generate_Restart_File(fname_original,neb_atom_coordinates);
701
702 /* generate input files */
703
704 generate_input_files(fname_original,iter);
705
706 for(h=1;h<=g;h++){
707 if(h==g && g==syou+1){
708 sprintf(fname1,"%s_%i",fname_original,myworld3+1+(h-1)*numprocs);
709 }
710 else{
711 sprintf(fname1,"%s_%i",fname_original,myworld1+1+(h-1)*numprocs);
712 }
713 argv[1] = fname1;
714 if((h==g) && (g==syou+1)){
715 neb_run( argv,MPI_CommWD3[myworld3],myworld3+1+(h-1)*numprocs,neb_atom_coordinates,
716 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
717
718 if (myid3==Host_ID){
719 Tmp_Grid_Origin[myworld3+1+(h-1)*numprocs][1] = Grid_Origin[1];
720 Tmp_Grid_Origin[myworld3+1+(h-1)*numprocs][2] = Grid_Origin[2];
721 Tmp_Grid_Origin[myworld3+1+(h-1)*numprocs][3] = Grid_Origin[3];
722 }
723 }
724 else{
725 neb_run( argv,MPI_CommWD1[myworld1],myworld1+1+(h-1)*numprocs,neb_atom_coordinates,
726 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
727
728 if(h==1){
729 for (i=0; i<=(NEB_Num_Images+1); i++){
730 for (j=0; j<=3; j++) Tmp_Grid_Origin[i][j] = 0.0;
731 }
732 }
733 if (myid1==Host_ID){
734 Tmp_Grid_Origin[myworld1+1+(h-1)*numprocs][1] = Grid_Origin[1];
735 Tmp_Grid_Origin[myworld1+1+(h-1)*numprocs][2] = Grid_Origin[2];
736 Tmp_Grid_Origin[myworld1+1+(h-1)*numprocs][3] = Grid_Origin[3];
737 }
738 }
739
740 /* MPI: All_Grid_Origin */
741
742 MPI_Barrier(MPI_COMM_WORLD);
743
744 }
745
746 MPI_Barrier(MPI_COMM_WORLD);
747
748 for (p=1; p<=NEB_Num_Images; p++){
749 MPI_Allreduce( &Tmp_Grid_Origin[p][0], &All_Grid_Origin[p][0],
750 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
751 }
752
753 /* MPI: neb_atom_coordinates */
754
755 for (p=1; p<=NEB_Num_Images; p++){
756 for (i=0; i<=atomnum; i++){
757 for (j=0; j<20; j++){
758 tmp_neb_atom_coordinates[p][i][j] = 0.0;
759 }
760 }
761 }
762
763 for(h=1;h<=g;h++){
764 if(h==g && g==syou+1){
765 if (myid3==Host_ID){
766 for (i=0; i<=atomnum; i++){
767 for (j=0; j<20; j++){
768 tmp_neb_atom_coordinates[myworld3+1+(h-1)*numprocs][i][j] = neb_atom_coordinates[myworld3+1+(h-1)*numprocs][i][j];;
769 }
770 }
771 }
772 }
773
774 else{
775 if (myid1==Host_ID){
776 for (i=0; i<=atomnum; i++){
777 for (j=0; j<20; j++){
778 tmp_neb_atom_coordinates[myworld1+1+(h-1)*numprocs][i][j] = neb_atom_coordinates[myworld1+1+(h-1)*numprocs][i][j];;
779 }
780 }
781 }
782 }
783 }
784
785 MPI_Barrier(MPI_COMM_WORLD);
786
787 for (p=1; p<=NEB_Num_Images; p++){
788 for (i=0; i<=atomnum; i++){
789
790 MPI_Allreduce( &tmp_neb_atom_coordinates[p][i][0], &neb_atom_coordinates[p][i][0],
791 20, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
792
793 }
794 }
795
796 MPI_Barrier(MPI_COMM_WORLD);
797
798 /* calculate the gradients defined by the NEB method */
799
800 Calc_NEB_Gradients(neb_atom_coordinates);
801
802 /* make a xyz file storing a set of structures of images at the current step */
803
804 Make_XYZ_File(fileXYZ,neb_atom_coordinates);
805
806 /* update the atomic structures of the images */
807
808 Update_Coordinates(iter,neb_atom_coordinates);
809
810 /* generate an input file for restarting */
811
812 Generate_Restart_File(fname_original,neb_atom_coordinates);
813
814 /* increment of iter */
815
816 iter++;
817
818 } while (MD_Opt_OK==0 && iter<=MD_IterNumber);
819
820 dtime(&b1time);
821
822 MPI_Barrier(MPI_COMM_WORLD);
823
824 /* freeing of arrays for the World1 */
825
826 dtime(&c0time);
827
828 MPI_Comm_free(&MPI_CommWD1[myworld1]);
829 free(MPI_CommWD1);
830 free(Comm_World_StartID1);
831 free(NPROCS_WD1);
832 free(Comm_World1);
833 free(NPROCS_ID1);
834
835 /* freeing of arrays for the World2 */
836
837 MPI_Comm_free(&MPI_CommWD2[myworld2]);
838 free(MPI_CommWD2);
839 free(Comm_World_StartID2);
840 free(NPROCS_WD2);
841 free(Comm_World2);
842 free(NPROCS_ID2);
843
844 /* freeing of arrays for the World3 */
845
846 if(amari!=0){
847
848 MPI_Comm_free(&MPI_CommWD3[myworld3]);
849 free(MPI_CommWD3);
850 free(Comm_WOrld_StartID3);
851 free(NPROCS_WD3);
852 free(Comm_World3);
853 free(NPROCS_ID3);
854
855 }
856
857 /* freeing of arrays */
858
859 free_arrays();
860
861 /* show a final message */
862 dtime(&c1time);
863
864 dtime(&TEtime);
865
866 printf("\nThe calculation was normally finished. (proc=%3d) TIME=%lf (s)\n",myid,(TEtime-TStime));
867
868 MPI_Finalize();
869 exit(0);
870 }
871
872 } /* if (PN==0) */
873
874 /****************************************************************************
875 If MD.NEB.Parallel.Number (PN) is specifed,
876 corresponding to 1<=MD.NEB.Parallel.Number (PN)
877 *****************************************************************************/
878
879 else{
880
881 dtime(&a0time);
882 int syou,amari,g,bb;
883
884 syou = NEB_Num_Images/PN;
885 amari = NEB_Num_Images%PN;
886
887 if (amari==0) g = syou;
888 else g = syou + 1;
889
890 Num_Comm_World1 = PN;
891
892 if(PN>numprocs){
893
894 if (myid==Host_ID) printf("PN is larger than the number of MPI processes.\n");
895
896 MPI_Finalize();
897 exit(0);
898 }
899
900 parallel1_flag = 1;
901
902 NPROCS_ID1 = (int*)malloc(sizeof(int)*numprocs);
903 Comm_World1 = (int*)malloc(sizeof(int)*numprocs);
904 NPROCS_WD1 = (int*)malloc(sizeof(int)*Num_Comm_World1);
905 Comm_World_StartID1 = (int*)malloc(sizeof(int)*Num_Comm_World1);
906 MPI_CommWD1 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World1);
907
908 Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World1, &myworld1, MPI_CommWD1,
909 NPROCS_ID1, Comm_World1, NPROCS_WD1, Comm_World_StartID1);
910
911 MPI_Comm_size(MPI_CommWD1[myworld1],&numprocs1);
912 MPI_Comm_rank(MPI_CommWD1[myworld1],&myid1);
913
914 /****************************************************
915 allocate processes to each image in the World2
916 ****************************************************/
917
918 Num_Comm_World2 = 2;
919
920 if ( Num_Comm_World2<=numprocs ){
921
922 parallel2_flag = 1;
923
924 NPROCS_ID2 = (int*)malloc(sizeof(int)*numprocs);
925 Comm_World2 = (int*)malloc(sizeof(int)*numprocs);
926 NPROCS_WD2 = (int*)malloc(sizeof(int)*Num_Comm_World2);
927 Comm_World_StartID2 = (int*)malloc(sizeof(int)*Num_Comm_World2);
928 MPI_CommWD2 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World2);
929
930 Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World2, &myworld2, MPI_CommWD2,
931 NPROCS_ID2, Comm_World2, NPROCS_WD2, Comm_World_StartID2);
932
933 MPI_Comm_size(MPI_CommWD2[myworld2],&numprocs2);
934 MPI_Comm_rank(MPI_CommWD2[myworld2],&myid2);
935 }
936 else{
937 parallel2_flag = 0;
938 myworld2 = 0;
939 }
940
941 /****************************************************
942 SCF calculations for the two terminal structures
943 ****************************************************/
944
945 /* check whether the restart is performed or not */
946
947 sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name);
948
949 if (Scf_RestartFromFile==1 && (fp = fopen(file_neb_utot,"rb"))!=NULL){
950
951 if (myid==Host_ID){
952 fread(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp);
953 fread(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp);
954 }
955
956 MPI_Bcast(&neb_atom_coordinates[0][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1);
957 MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1);
958
959 fclose(fp);
960 }
961
962 /* if the restart is not performed, perform the SCF calulations for the two terminals */
963
964 else {
965
966 /* generate input files */
967
968 generate_input_files(fname_original,1);
969
970 /***********************************
971 if (PN==1)
972 ***********************************/
973
974 if (PN==1){
975
976 for (i=0; i<=1; i++){
977
978 if (i==0) index_images = 0;
979 else index_images = NEB_Num_Images + 1;
980
981 sprintf(fname1,"%s_%i",fname_original,index_images);
982 argv[1] = fname1;
983
984 neb_run( argv,MPI_CommWD2[myworld2],index_images,neb_atom_coordinates,
985 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
986 }
987 }
988
989 /***********************************
990 if (2<=PN)
991 ***********************************/
992
993 else {
994
995 if (myworld2==0)
996 index_images = 0;
997 else
998 index_images = NEB_Num_Images + 1;
999
1000 sprintf(fname1,"%s_%i",fname_original,index_images);
1001 argv[1] = fname1;
1002
1003 neb_run( argv,MPI_CommWD2[myworld2],index_images,neb_atom_coordinates,
1004 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
1005
1006 /* find a representative ID in the top world for ID=0 in the world2 */
1007
1008 i = 0; po = 0;
1009 do {
1010
1011 if (Comm_World2[i]==0){
1012 ID0 = i;
1013 po = 1;
1014 }
1015 i++;
1016 } while (po==0);
1017
1018 /* find a representative ID in the top world for ID=1 in the world2 */
1019
1020 i = 0; po = 0;
1021 do {
1022
1023 if (Comm_World2[i]==1){
1024 ID1 = i;
1025 po = 1;
1026 }
1027 i++;
1028 } while (po==0);
1029
1030 /* MPI broadcast of neb_atom_coordinates[0] and neb_atom_coordinates[NEB_Num_Images+1] */
1031
1032 MPI_Barrier(MPI_COMM_WORLD);
1033 for (i=0; i<=atomnum; i++){
1034 MPI_Bcast(&neb_atom_coordinates[0][i][0], 20, MPI_DOUBLE, ID0, MPI_COMM_WORLD);
1035 MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][i][0], 20, MPI_DOUBLE, ID1, MPI_COMM_WORLD);
1036 }
1037
1038 } /* else which corresponds to if (2<=PN) */
1039
1040 /* save *_0_rst/*.neb.utot in binary mode */
1041
1042 if (myid==Host_ID){
1043
1044 sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name);
1045
1046 if ((fp = fopen(file_neb_utot,"wb")) != NULL){
1047
1048 fwrite(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp);
1049 fwrite(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp);
1050 fclose(fp);
1051 }
1052 else{
1053 printf("Could not open a file %s in neb\n",file_neb_utot);
1054 }
1055 }
1056
1057 } /* else */
1058
1059 /****************************************************
1060 make amari world
1061 ****************************************************/
1062
1063 if (amari!=0){
1064
1065 Num_Comm_World3 = amari;
1066
1067 NPROCS_ID3 = (int*)malloc(sizeof(int)*numprocs);
1068 Comm_World3 = (int*)malloc(sizeof(int)*numprocs);
1069 NPROCS_WD3 = (int*)malloc(sizeof(int)*Num_Comm_World3);
1070 Comm_WOrld_StartID3 = (int*)malloc(sizeof(int)*Num_Comm_World3);
1071 MPI_CommWD3 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World3);
1072
1073 Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World3, &myworld3, MPI_CommWD3,
1074 NPROCS_ID3, Comm_World3, NPROCS_WD3, Comm_WOrld_StartID3);
1075
1076 MPI_Comm_size(MPI_CommWD3[myworld3],&numprocs3);
1077 MPI_Comm_rank(MPI_CommWD3[myworld3],&myid3);
1078 }
1079 /* amari world */
1080
1081 dtime(&a1time);
1082
1083 /****************************************************
1084 optimization for finding a minimum energy path
1085 connecting the two terminal structures.
1086 ****************************************************/
1087
1088 iter = 1;
1089 MD_Opt_OK = 0;
1090 dtime(&b0time);
1091
1092 do {
1093
1094 /* if iter==1, generate an input file for restarting */
1095
1096 if (iter==1) Generate_Restart_File(fname_original,neb_atom_coordinates);
1097
1098 /* generate input files */
1099
1100 generate_input_files(fname_original,iter);
1101
1102 for (h=1; h<=g; h++){
1103
1104 if( h==g && g==(syou+1) ){
1105 sprintf(fname1,"%s_%i",fname_original,myworld3+1+(h-1)*PN);
1106 }
1107 else{
1108 sprintf(fname1,"%s_%i",fname_original,myworld1+1+(h-1)*PN);
1109 }
1110
1111 argv[1] = fname1;
1112
1113 if( h==g && g==(syou+1) ){
1114
1115 neb_run( argv,MPI_CommWD3[myworld3],myworld3+1+(h-1)*PN,neb_atom_coordinates,
1116 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
1117
1118 if (myid3==Host_ID){
1119 Tmp_Grid_Origin[myworld3+1+(h-1)*PN][1] = Grid_Origin[1];
1120 Tmp_Grid_Origin[myworld3+1+(h-1)*PN][2] = Grid_Origin[2];
1121 Tmp_Grid_Origin[myworld3+1+(h-1)*PN][3] = Grid_Origin[3];
1122 }
1123 }
1124
1125 else{
1126
1127 neb_run( argv,MPI_CommWD1[myworld1],myworld1+1+(h-1)*PN,neb_atom_coordinates,
1128 WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB );
1129
1130 if(h==1){
1131 for (i=0; i<=(NEB_Num_Images+1); i++){
1132 for (j=0; j<=3; j++) Tmp_Grid_Origin[i][j] = 0.0;
1133 }
1134 }
1135 if (myid1==Host_ID){
1136 Tmp_Grid_Origin[myworld1+1+(h-1)*PN][1] = Grid_Origin[1];
1137 Tmp_Grid_Origin[myworld1+1+(h-1)*PN][2] = Grid_Origin[2];
1138 Tmp_Grid_Origin[myworld1+1+(h-1)*PN][3] = Grid_Origin[3];
1139 }
1140 }
1141
1142 /* MPI: All_Grid_Origin */
1143 MPI_Barrier(MPI_COMM_WORLD);
1144
1145 } /* for (h=1; h<=g; h++) */
1146
1147 MPI_Barrier(MPI_COMM_WORLD);
1148
1149 for (p=1; p<=NEB_Num_Images; p++){
1150 MPI_Allreduce( &Tmp_Grid_Origin[p][0], &All_Grid_Origin[p][0],
1151 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1152 }
1153
1154 /* MPI: neb_atom_coordinates */
1155
1156 for (p=1; p<=NEB_Num_Images; p++){
1157 for (i=0; i<=atomnum; i++){
1158 for (j=0; j<20; j++){
1159 tmp_neb_atom_coordinates[p][i][j] = 0.0;
1160 }
1161 }
1162 }
1163
1164 for(h=1;h<=g;h++){
1165 if(h==g && g==syou+1){
1166 if (myid3==Host_ID){
1167 for (i=0; i<=atomnum; i++){
1168 for (j=0; j<20; j++){
1169 tmp_neb_atom_coordinates[myworld3+1+(h-1)*PN][i][j] = neb_atom_coordinates[myworld3+1+(h-1)*PN][i][j];;
1170 }
1171 }
1172 }
1173 }
1174
1175 else{
1176 if (myid1==Host_ID){
1177 for (i=0; i<=atomnum; i++){
1178 for (j=0; j<20; j++){
1179 tmp_neb_atom_coordinates[myworld1+1+(h-1)*PN][i][j] = neb_atom_coordinates[myworld1+1+(h-1)*PN][i][j];;
1180 }
1181 }
1182 }
1183 }
1184 }
1185
1186 MPI_Barrier(MPI_COMM_WORLD);
1187
1188 for (p=1; p<=NEB_Num_Images; p++){
1189 for (i=0; i<=atomnum; i++){
1190
1191 MPI_Allreduce( &tmp_neb_atom_coordinates[p][i][0], &neb_atom_coordinates[p][i][0],
1192 20, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1193
1194 }
1195 }
1196
1197 MPI_Barrier(MPI_COMM_WORLD);
1198
1199 /* calculate the gradients defined by the NEB method */
1200
1201 Calc_NEB_Gradients(neb_atom_coordinates);
1202
1203 /* make a xyz file storing a set of structures of images at the current step */
1204
1205 Make_XYZ_File(fileXYZ,neb_atom_coordinates);
1206
1207 /* update the atomic structures of the images */
1208
1209 Update_Coordinates(iter,neb_atom_coordinates);
1210
1211 /* generate an input file for restarting */
1212
1213 Generate_Restart_File(fname_original,neb_atom_coordinates);
1214
1215 /* increment of iter */
1216
1217 iter++;
1218
1219 } while (MD_Opt_OK==0 && iter<=MD_IterNumber);
1220
1221 dtime(&b1time);
1222
1223 MPI_Barrier(MPI_COMM_WORLD);
1224
1225 /* freeing of arrays for the World1 */
1226
1227 dtime(&c0time);
1228
1229 MPI_Comm_free(&MPI_CommWD1[myworld1]);
1230 free(MPI_CommWD1);
1231 free(Comm_World_StartID1);
1232 free(NPROCS_WD1);
1233 free(Comm_World1);
1234 free(NPROCS_ID1);
1235
1236 /* freeing of arrays for the World2 */
1237
1238 MPI_Comm_free(&MPI_CommWD2[myworld2]);
1239 free(MPI_CommWD2);
1240 free(Comm_World_StartID2);
1241 free(NPROCS_WD2);
1242 free(Comm_World2);
1243 free(NPROCS_ID2);
1244
1245 /* freeing of arrays for the World3 */
1246
1247 if(amari!=0){
1248
1249 MPI_Comm_free(&MPI_CommWD3[myworld3]);
1250 free(MPI_CommWD3);
1251 free(Comm_WOrld_StartID3);
1252 free(NPROCS_WD3);
1253 free(Comm_World3);
1254 free(NPROCS_ID3);
1255
1256 }
1257
1258 /* freeing of arrays */
1259
1260 free_arrays();
1261
1262 /* show a final message */
1263 dtime(&c1time);
1264
1265 dtime(&TEtime);
1266
1267 printf("\nThe calculation was normally finished. (proc=%3d) TIME=%lf (s)\n",myid,(TEtime-TStime));
1268
1269 MPI_Finalize();
1270 exit(0);
1271
1272
1273 }
1274
1275 }
1276
1277
1278
Generate_Restart_File(char fname_original[YOUSO10],double *** neb_atom_coordinates)1279 void Generate_Restart_File(char fname_original[YOUSO10], double ***neb_atom_coordinates)
1280 {
1281 int i,p,Gc_AN,c,n1,k,po;
1282 int restart_flag;
1283 int unit_flag;
1284 double c1,c2,c3;
1285 double tmpxyz[4];
1286 char st[800];
1287 char st1[800];
1288 char rm_operate[YOUSO10];
1289 char fname1[YOUSO10];
1290 char fname2[YOUSO10];
1291 char keyword[YOUSO10];
1292 FILE *fp1,*fp2;
1293 char *tp;
1294 int numprocs,myid;
1295
1296 /* MPI */
1297 MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
1298 MPI_Comm_rank(MPI_COMM_WORLD,&myid);
1299
1300 /* generate the input files */
1301
1302 if (myid==Host_ID){
1303
1304 /* initialize */
1305
1306 restart_flag = 0;
1307 unit_flag = 0;
1308
1309 /* the new input file */
1310
1311 sprintf(fname1,"%s#",fname_original);
1312 fp1 = fopen(fname1,"w");
1313 fseek(fp1,0,SEEK_END);
1314
1315 /* the original input file */
1316
1317 fp2 = fopen(fname_original,"r");
1318
1319 if (fp2!=NULL){
1320
1321 while (fgets(st,800,fp2)!=NULL){
1322
1323 string_tolower(st,st1);
1324
1325 /* find the specification of <atoms.speciesandcoordinates */
1326
1327 if (strncmp(st1,"<atoms.speciesandcoordinates",28)==0){
1328
1329 fprintf(fp1,"%s",st);
1330
1331 /* replace the atomic coordinates */
1332
1333 for (i=1; i<=atomnum; i++){
1334
1335 fgets(st,800,fp2);
1336 string_tolower(st,st1);
1337
1338 /* serial number */
1339 tp = strtok(st, " ");
1340 if (tp!=NULL) fprintf(fp1,"%4s",tp);
1341
1342 /* name of species */
1343 tp =strtok(NULL, " ");
1344 if (tp!=NULL) fprintf(fp1," %4s",tp);
1345
1346 /* x-coordinate */
1347 tp =strtok(NULL, " ");
1348 fprintf(fp1," %18.14f",neb_atom_coordinates[0][i][1]);
1349
1350 /* y-coordinate */
1351 tp =strtok(NULL, " ");
1352 fprintf(fp1," %18.14f",neb_atom_coordinates[0][i][2]);
1353
1354 /* z-coordinate */
1355 tp =strtok(NULL, " ");
1356 fprintf(fp1," %18.14f",neb_atom_coordinates[0][i][3]);
1357
1358 while (tp!=NULL){
1359 tp =strtok(NULL, " ");
1360 if (tp!=NULL) fprintf(fp1," %s",tp);
1361 }
1362
1363 }
1364 }
1365
1366 /* find the specification of <atoms.speciesandcoordinates */
1367
1368 else if (strncmp(st1,"<neb.atoms.speciesandcoordinates",32)==0){
1369
1370 fprintf(fp1,"%s",st);
1371
1372 /* replace the atomic coordinates */
1373
1374 for (i=1; i<=atomnum; i++){
1375
1376 fgets(st,800,fp2);
1377 string_tolower(st,st1);
1378
1379 /* serial number */
1380 tp = strtok(st, " ");
1381 if (tp!=NULL) fprintf(fp1,"%4s",tp);
1382
1383 /* name of species */
1384 tp =strtok(NULL, " ");
1385 if (tp!=NULL) fprintf(fp1," %4s",tp);
1386
1387 /* x-coordinate */
1388 tp =strtok(NULL, " ");
1389 fprintf(fp1," %18.14f",neb_atom_coordinates[NEB_Num_Images+1][i][1]);
1390
1391 /* y-coordinate */
1392 tp =strtok(NULL, " ");
1393 fprintf(fp1," %18.14f",neb_atom_coordinates[NEB_Num_Images+1][i][2]);
1394
1395 /* z-coordinate */
1396 tp =strtok(NULL, " ");
1397 fprintf(fp1," %18.14f",neb_atom_coordinates[NEB_Num_Images+1][i][3]);
1398
1399 while (tp!=NULL){
1400 tp =strtok(NULL, " ");
1401 if (tp!=NULL) fprintf(fp1," %s",tp);
1402 }
1403
1404 }
1405 }
1406
1407 /* Atoms.SpeciesAndCoordinates.Unit */
1408
1409 else if (strncmp(st1,"atoms.speciesandcoordinates.unit",32)==0){
1410 fprintf(fp1,"Atoms.SpeciesAndCoordinates.Unit AU\n");
1411 unit_flag = 1;
1412 }
1413
1414 /* sc.restart */
1415
1416 else if (strncmp(st1,"scf.restart",11)==0){
1417 fprintf(fp1,"scf.restart on\n");
1418 restart_flag = 1;
1419 }
1420
1421 else {
1422
1423 po = 0;
1424
1425 for (p=1; p<=NEB_Num_Images; p++){
1426
1427 /* find the specification of <neb%d.atoms.speciesandcoordinates and forget it */
1428
1429 k = (int)log10(p) + 1;
1430 sprintf(keyword,"<neb%i.atoms.speciesandcoordinates",p);
1431
1432 if (strncmp(st1,keyword,32+k)==0){
1433
1434 po = 1;
1435
1436 /* get atomic coordinates and forget them */
1437 for (i=1; i<=atomnum; i++) fgets(st,800,fp2);
1438
1439 /* get neb%i.atoms.speciesandcoordinates> and forget it */
1440 fgets(st,800,fp2);
1441
1442 }
1443 } /* p */
1444
1445 /* other cases */
1446
1447 if (po==0) fprintf(fp1,"%s",st);
1448
1449 } /* else */
1450 } /* while */
1451
1452 /* close fp2 */
1453 fclose(fp2);
1454 }
1455
1456 /* add the restart flag if it was not found. */
1457
1458 if (restart_flag==0){
1459 fprintf(fp1,"\n\nscf.restart on\n");
1460 }
1461
1462 /* add Atoms.SpeciesAndCoordinates.Unit if it was not found. */
1463
1464 if (unit_flag==0){
1465 fprintf(fp1,"Atoms.SpeciesAndCoordinates.Unit AU\n");
1466 }
1467
1468 /* add atomic coordinates of the images */
1469
1470 for (p=1; p<=NEB_Num_Images; p++){
1471
1472 fp2 = fopen(fname_original,"r");
1473
1474 if (fp2!=NULL){
1475
1476 while (fgets(st,800,fp2)!=NULL){
1477
1478 string_tolower(st,st1);
1479
1480 if (strncmp(st1,"<atoms.speciesandcoordinates",28)==0){
1481
1482 fprintf(fp1,"\n<NEB%d.Atoms.SpeciesAndCoordinates\n",p);
1483
1484 for (k=1; k<=atomnum; k++){
1485
1486 fgets(st,800,fp2);
1487 string_tolower(st,st1);
1488
1489 /* serial number */
1490 tp = strtok(st, " ");
1491 if (tp!=NULL) fprintf(fp1,"%4s",tp);
1492
1493 /* name of species */
1494 tp =strtok(NULL, " ");
1495 if (tp!=NULL) fprintf(fp1," %4s",tp);
1496
1497 /* x-coordinate */
1498 tp =strtok(NULL, " ");
1499 fprintf(fp1," %18.14f",neb_atom_coordinates[p][k][1]);
1500
1501 /* y-coordinate */
1502 tp =strtok(NULL, " ");
1503 fprintf(fp1," %18.14f",neb_atom_coordinates[p][k][2]);
1504
1505 /* z-coordinate */
1506 tp =strtok(NULL, " ");
1507 fprintf(fp1," %18.14f",neb_atom_coordinates[p][k][3]);
1508
1509 while (tp!=NULL){
1510 tp =strtok(NULL, " ");
1511 if (tp!=NULL) fprintf(fp1," %s",tp);
1512 }
1513
1514 } /* k */
1515
1516 }
1517
1518 else if (strncmp(st1,"atoms.speciesandcoordinates>",28)==0){
1519 fprintf(fp1,"NEB%d.Atoms.SpeciesAndCoordinates>\n",p);
1520 }
1521 }
1522
1523 /* close fp2 */
1524 fclose(fp2);
1525
1526 } /* if (fp2!=NULL){ */
1527 } /* p */
1528
1529 /* fclose */
1530 fclose(fp1);
1531
1532 } /* if (myid==Host_ID) */
1533
1534 }
1535
1536
1537
1538
1539
Make_XYZ_File(char fname[YOUSO10],double *** neb_atom_coordinates)1540 void Make_XYZ_File(char fname[YOUSO10], double ***neb_atom_coordinates)
1541 {
1542 FILE *fp;
1543 int p,k,i,j;
1544 int myid;
1545 char file_neb_utot[YOUSO10];
1546
1547 MPI_Comm_rank(MPI_COMM_WORLD,&myid);
1548
1549 if (myid==Host_ID){
1550
1551 /* save *.neb.xyz */
1552
1553 if ((fp = fopen(fname,"w")) != NULL){
1554
1555 for (p=0; p<=(NEB_Num_Images+1); p++){
1556
1557 fprintf(fp,"%i \n",atomnum);
1558 fprintf(fp," Image index = %3d Energy= %8.5f (Hatree)\n",p,neb_atom_coordinates[p][0][0]);
1559
1560 for (k=1; k<=atomnum; k++){
1561
1562 i = WhatSpecies_NEB[k];
1563 j = Spe_WhatAtom_NEB[i];
1564
1565 fprintf(fp,"%4s %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f\n",
1566 Atom_Symbol[j],
1567 neb_atom_coordinates[p][k][1]*BohrR,
1568 neb_atom_coordinates[p][k][2]*BohrR,
1569 neb_atom_coordinates[p][k][3]*BohrR,
1570 -neb_atom_coordinates[p][k][4],
1571 -neb_atom_coordinates[p][k][5],
1572 -neb_atom_coordinates[p][k][6]);
1573 }
1574 }
1575 fclose(fp);
1576 }
1577 else{
1578 printf("failure of saving the xyz file.\n");
1579 fclose(fp);
1580 }
1581
1582 } /* if (myid==Host_ID) */
1583
1584 }
1585
1586
1587
1588
1589
Update_Coordinates(int iter,double *** neb_atom_coordinates)1590 void Update_Coordinates(int iter, double ***neb_atom_coordinates)
1591 {
1592 int p,j,k;
1593 double norm,dis,dif;
1594 double Max_Force,Max_Step;
1595 double tmp0,sum,sum_E,sum_dis;
1596 int numprocs,myid;
1597 char fileOPT[YOUSO10];
1598 FILE *fp_OPT;
1599 char fileE[YOUSO10];
1600 FILE *fp_E;
1601
1602 MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
1603 MPI_Comm_rank(MPI_COMM_WORLD,&myid);
1604
1605 /*********************************************
1606 calculate the norm of gradient
1607 *********************************************/
1608
1609 norm = 0.0;
1610 for (p=1; p<=NEB_Num_Images; p++){
1611 for (j=1; j<=atomnum; j++){
1612 for (k=1; k<=3; k++){
1613 norm += neb_atom_coordinates[p][j][k+3]*neb_atom_coordinates[p][j][k+3];
1614 }
1615 }
1616 }
1617
1618 norm = sqrt(norm);
1619
1620 /****************************************************
1621 find the maximum value of force
1622 ****************************************************/
1623
1624 Max_Force = 0.0;
1625 for (p=1; p<=NEB_Num_Images; p++){
1626 for (j=1; j<=atomnum; j++){
1627
1628 sum = 0.0;
1629 for (k=1; k<=3; k++){
1630 tmp0 = neb_atom_coordinates[p][j][k+3];
1631 sum += tmp0*tmp0;
1632 }
1633 sum = sqrt(sum);
1634 if (Max_Force<sum) Max_Force = sum;
1635 }
1636 }
1637
1638 if (Max_Force<MD_Opt_criterion) MD_Opt_OK = 1;
1639
1640 /*********************************************
1641 update atomic coordinates
1642 *********************************************/
1643
1644 Max_Step = DIIS_BFGS(iter, neb_atom_coordinates,norm,Max_Force);
1645
1646 /*
1647 Steepest_Decent(iter, neb_atom_coordinates,norm,Max_Force);
1648 */
1649
1650 /*********************************************
1651 save the history of optimization
1652 *********************************************/
1653
1654 if (myid==Host_ID){
1655
1656 sum_E = 0.0;
1657 for (p=1; p<=NEB_Num_Images; p++){
1658 sum_E += neb_atom_coordinates[p][0][0];
1659 }
1660
1661 sprintf(fileOPT,"%s%s.neb.opt",filepath,system_name);
1662
1663 if ((fp_OPT = fopen(fileOPT,"a")) != NULL){
1664
1665 if (iter==1){
1666 fprintf(fp_OPT,"\n");
1667 fprintf(fp_OPT,"#***********************************************************\n");
1668 fprintf(fp_OPT,"#***********************************************************\n");
1669 fprintf(fp_OPT,"# History of optimization by the NEB method \n");
1670 fprintf(fp_OPT,"#***********************************************************\n");
1671 fprintf(fp_OPT,"#***********************************************************\n");
1672 fprintf(fp_OPT,"#\n");
1673 fprintf(fp_OPT,"# iter SD_scaling |Maximum force| Maximum step Norm Sum of Total Energy of Images\n");
1674 fprintf(fp_OPT,"# (Hartree/Bohr) (Ang) (Hartree/Bohr) (Hartree) \n\n");
1675 }
1676
1677 fprintf(fp_OPT," %3d %15.8f %15.8f %15.8f %15.8f %15.8f\n",
1678 iter,SD_scaling,Max_Force,Max_Step,norm,sum_E);
1679
1680 fclose(fp_OPT);
1681 }
1682 else
1683 printf("Error in saving the neb.opt file\n");
1684
1685 } /* if (myid==Host_ID) */
1686
1687 /*********************************************
1688 save the total DFT energy of each image
1689 and iterdistance between them.
1690 *********************************************/
1691
1692 if (myid==Host_ID){
1693
1694 sprintf(fileE,"%s%s.neb.ene",filepath,system_name);
1695
1696 if ((fp_E = fopen(fileE,"w")) != NULL){
1697
1698 fprintf(fp_E,"#\n");
1699 fprintf(fp_E,"# 1st column: index of images, where 0 and MD.NEB.Number.Images+1 are the terminals\n");
1700 fprintf(fp_E,"# 2nd column: Total energy (Hartree) of each image\n");
1701 fprintf(fp_E,"# 3rd column: distance (Bohr) between neighbors\n");
1702 fprintf(fp_E,"# 4th column: distance (Bohr) from the image of the index 0\n");
1703 fprintf(fp_E,"# 5th column: x distance\n " );
1704 fprintf(fp_E,"#\n");
1705
1706 sum_dis = 0.0;
1707
1708 for (p=0; p<=(NEB_Num_Images+1); p++){
1709
1710 if (p!=0){
1711 dis = 0.0;
1712 for (j=1; j<=atomnum; j++){
1713 for (k=1; k<=3; k++){
1714 dif = neb_atom_coordinates[p][j][k] - neb_atom_coordinates[p-1][j][k];
1715 dis += dif*dif;
1716 }
1717 }
1718 }
1719 else{
1720 dis = 0.0;
1721 }
1722
1723 sum_dis += dis;
1724
1725 fprintf(fp_E," %3d %15.8f %15.8f %15.8f %15.8f\n",
1726 p, neb_atom_coordinates[p][0][0], dis, sum_dis, neb_atom_coordinates[p][2][1]*BohrR-neb_atom_coordinates[p][1][1]*BohrR);
1727
1728 }
1729
1730 fclose(fp_E);
1731 }
1732 else
1733 printf("Error in saving the neb.ene file\n");
1734
1735 }
1736
1737 }
1738
1739
1740
1741
DIIS_BFGS(int iter,double *** neb_atom_coordinates,double norm,double Max_Force)1742 double DIIS_BFGS(int iter, double ***neb_atom_coordinates, double norm, double Max_Force)
1743 {
1744 int p,i,j,k,m,n,dim;
1745 int p1,i1,j1,p2,i2,j2;
1746 double sum,tmp1,tmp2,tmp,c0,c1;
1747 double Max_Step;
1748 double SD_min,SD_max,SD_init,dt;
1749 double dif_g,dif_x,dif_x1,dif_x2;
1750 double sum0,sum1;
1751 char *JOBB="L";
1752 double *A,*B,max_A,RR;
1753 double *work;
1754 int *ipiv;
1755 INTEGER LDA,LWORK,info,N;
1756
1757 /****************************************************
1758 SD method
1759 ****************************************************/
1760
1761 if (iter<OptStartDIIS){
1762
1763 /*************************************************************
1764 store coordinates and gradients
1765 *************************************************************/
1766
1767 for (p=1; p<=NEB_Num_Images; p++){
1768 for (j=1; j<=atomnum; j++){
1769 for (k=1; k<=6; k++){
1770 His_neb_atom_coordinates[1][p][j][k] = His_neb_atom_coordinates[0][p][j][k];
1771 His_neb_atom_coordinates[0][p][j][k] = neb_atom_coordinates[p][j][k];
1772 }
1773 }
1774 }
1775
1776 /****************************************************
1777 set up SD_scaling
1778 ****************************************************/
1779
1780 dt = 41.3411*4.0;
1781 SD_init = dt*dt/1836.1526;
1782 SD_max = SD_init*10.0; /* default 10 */
1783 SD_min = SD_init*0.005; /* default 0.2 */
1784
1785 if (iter==1){
1786
1787 SD_scaling_user = Max_Force/BohrR/5.0;
1788 SD_scaling = SD_scaling_user/(Max_Force+1.0e-10);
1789
1790 if (SD_max<SD_scaling) SD_scaling = SD_max;
1791 if (SD_scaling<SD_min) SD_scaling = SD_min;
1792
1793 Past_Norm[1] = norm;
1794 }
1795 else {
1796
1797 if (Past_Norm[1]<norm && iter%4==1){
1798 SD_scaling = SD_scaling/4.0;
1799 }
1800 else if (Past_Norm[1]<Past_Norm[2] && norm<Past_Norm[1] && iter%4==1){
1801 SD_scaling = SD_scaling*1.2;
1802 }
1803
1804 if (SD_max<SD_scaling) SD_scaling = SD_max;
1805 if (SD_scaling<SD_min) SD_scaling = SD_min;
1806
1807 Past_Norm[5] = Past_Norm[4];
1808 Past_Norm[4] = Past_Norm[3];
1809 Past_Norm[3] = Past_Norm[2];
1810 Past_Norm[2] = Past_Norm[1];
1811 Past_Norm[1] = norm;
1812 }
1813
1814 /*************************************************************
1815 update coordinate by the simple steepest decent (SD) method
1816 *************************************************************/
1817
1818 Max_Step = 0.0;
1819 for (p=1; p<=NEB_Num_Images; p++){
1820 for (j=1; j<=atomnum; j++){
1821 for (k=1; k<=3; k++){
1822
1823 tmp = SD_scaling*neb_atom_coordinates[p][j][k+3];
1824 neb_atom_coordinates[p][j][k] -= tmp;
1825 if (Max_Step<fabs(tmp)) Max_Step = fabs(tmp);
1826 }
1827 }
1828 }
1829
1830 } /* if (iter<OptStartDIIS) */
1831
1832 /****************************************************
1833 DIIS + BFGS method
1834 ****************************************************/
1835
1836 else {
1837
1838 /*************************************************************
1839 store coordinates and gradients
1840 *************************************************************/
1841
1842 for (p=1; p<=NEB_Num_Images; p++){
1843 for (j=1; j<=atomnum; j++){
1844 for (k=1; k<=6; k++){
1845
1846 for (m=(M_GDIIS_HISTORY-1); 0<m; m--){
1847 His_neb_atom_coordinates[m][p][j][k] = His_neb_atom_coordinates[m-1][p][j][k];
1848 }
1849
1850 His_neb_atom_coordinates[0][p][j][k] = neb_atom_coordinates[p][j][k];
1851 }
1852 }
1853 }
1854
1855 /*************************************************************
1856 estimate an optimum coordinates by the DIIS method
1857 *************************************************************/
1858
1859 dim = iter - OptStartDIIS + 1;
1860 if (M_GDIIS_HISTORY<dim) dim = M_GDIIS_HISTORY;
1861
1862 /* allocation of arrays */
1863
1864 A = (double*)malloc(sizeof(double)*(dim+1)*(dim+1));
1865 B = (double*)malloc(sizeof(double)*(dim+1));
1866
1867 /* construct the A matrix */
1868
1869 for (m=0; m<dim; m++){
1870 for (n=0; n<dim; n++){
1871
1872 sum = 0.0;
1873
1874 for (p=1; p<=NEB_Num_Images; p++){
1875 for (j=1; j<=atomnum; j++){
1876 for (k=4; k<=6; k++){
1877 sum += His_neb_atom_coordinates[m][p][j][k]*His_neb_atom_coordinates[n][p][j][k];
1878 }
1879 }
1880 }
1881
1882 A[m*(dim+1)+n] = sum;
1883 }
1884 }
1885
1886 /* find max of A and scale elements by it. */
1887
1888 max_A = 0.0;
1889 for (m=0; m<dim; m++){
1890 for (n=0; n<dim; n++){
1891 RR = fabs(A[m*(dim+1)+n]) ;
1892 if (max_A<RR ) max_A = RR;
1893 }
1894 }
1895
1896 max_A = 1.0/max_A;
1897
1898 for (m=0; m<dim; m++){
1899 for (n=0; n<dim; n++){
1900 A[m*(dim+1)+n] *= max_A;
1901 }
1902 }
1903
1904 for (m=0; m<dim; m++){
1905 A[m*(dim+1)+dim] = 1.0;
1906 A[dim*(dim+1)+m] = 1.0;
1907 }
1908 A[dim*(dim+1)+dim] = 0.0;
1909
1910 /*
1911 printf("Mat A\n");
1912 for (m=0; m<=dim; m++){
1913 for (n=0; n<=dim; n++){
1914 printf("%10.5f ",A[m*(dim+1)+n]);
1915 }
1916 printf("\n");
1917 }
1918 */
1919
1920 for (m=0; m<dim; m++) B[m] = 0.0;
1921 B[dim] = 1.0;
1922
1923 /* solve Ax = B */
1924
1925 N = dim + 1;
1926 LDA = dim + 1;
1927 LWORK = dim + 1;
1928 work = (double*)malloc(sizeof(double)*LWORK);
1929 ipiv = (int*)malloc(sizeof(int)*(dim+1));
1930 i = 1;
1931
1932 F77_NAME(dsysv,DSYSV)( JOBB, &N, &i, A, &LDA, ipiv, B, &LDA, work, &LWORK, &info);
1933
1934 if (info!=0) {
1935 printf(" ERROR in dsysv_, info=%d\n",info);
1936 MD_Opt_OK =1;
1937 goto Last_Step;
1938 }
1939
1940 /*
1941 for (m=0; m<dim; m++){
1942 printf("m=%2d B=%15.12f\n",m,B[m]);
1943 }
1944 */
1945
1946 /* update the coordinates by the DIIS method */
1947
1948 for (p=1; p<=NEB_Num_Images; p++){
1949 for (j=1; j<=atomnum; j++){
1950 for (k=1; k<=3; k++){
1951
1952 sum0 = 0.0;
1953 sum1 = 0.0;
1954 for (m=0; m<dim; m++){
1955 sum0 += B[m]*His_neb_atom_coordinates[m][p][j][k ];
1956 sum1 += B[m]*His_neb_atom_coordinates[m][p][j][k+3];
1957 }
1958
1959 neb_atom_coordinates[p][j][k] = sum0;
1960 neb_atom_coordinates[p][j][k+3] = sum1;
1961 }
1962 }
1963 }
1964
1965 Last_Step:
1966
1967 /* freeing of arrays */
1968
1969 free(work);
1970 free(ipiv);
1971 free(A);
1972 free(B);
1973
1974 /*************************************************************
1975 update an approximate Hessian by the BFGS method
1976 *************************************************************/
1977
1978 if (iter==OptStartDIIS){
1979
1980 for (p1=1; p1<=NEB_Num_Images; p1++){
1981 for (i1=1; i1<=atomnum; i1++){
1982 for (j1=1; j1<=3; j1++){
1983 for (p2=1; p2<=NEB_Num_Images; p2++){
1984 for (i2=1; i2<=atomnum; i2++){
1985 for (j2=1; j2<=3; j2++){
1986 InvHess[p1][i1][j1][p2][i2][j2] = 0.0;
1987 }
1988 }
1989 }
1990 InvHess[p1][i1][j1][p1][i1][j1] = 1.0;
1991 }
1992 }
1993 }
1994 }
1995
1996 else {
1997
1998 /* (H^(n-1))^{-1} Delta g^(n) */
1999
2000 for (p1=1; p1<=NEB_Num_Images; p1++){
2001 for (i1=1; i1<=atomnum; i1++){
2002 for (j1=1; j1<=3; j1++){
2003
2004 sum = 0.0;
2005 for (p2=1; p2<=NEB_Num_Images; p2++){
2006 for (i2=1; i2<=atomnum; i2++){
2007 for (j2=1; j2<=3; j2++){
2008 dif_g = His_neb_atom_coordinates[0][p2][i2][j2+3]-His_neb_atom_coordinates[1][p2][i2][j2+3];
2009 sum += InvHess[p1][i1][j1][p2][i2][j2]*dif_g;
2010 }
2011 }
2012 }
2013
2014 Vec0[p1][i1][j1] = sum;
2015 }
2016 }
2017 }
2018
2019 /* < Delta x^(n) | Delta g^(n) > */
2020
2021 tmp1 = 0.0;
2022 for (p2=1; p2<=NEB_Num_Images; p2++){
2023 for (i2=1; i2<=atomnum; i2++){
2024 for (j2=1; j2<=3; j2++){
2025
2026 dif_x = His_neb_atom_coordinates[0][p2][i2][j2 ]-His_neb_atom_coordinates[1][p2][i2][j2 ];
2027 dif_g = His_neb_atom_coordinates[0][p2][i2][j2+3]-His_neb_atom_coordinates[1][p2][i2][j2+3];
2028 tmp1 += dif_x*dif_g;
2029 }
2030 }
2031 }
2032
2033 /* < Delta g^(n) | Vec0 > */
2034
2035 tmp2 = 0.0;
2036 for (p2=1; p2<=NEB_Num_Images; p2++){
2037 for (i2=1; i2<=atomnum; i2++){
2038 for (j2=1; j2<=3; j2++){
2039
2040 dif_g = His_neb_atom_coordinates[0][p2][i2][j2+3]-His_neb_atom_coordinates[1][p2][i2][j2+3];
2041 tmp2 += dif_g*Vec0[p2][i2][j2];
2042 }
2043 }
2044 }
2045
2046 /* update InvHess */
2047
2048 c0 = (tmp1 + tmp2)/(tmp1*tmp1);
2049 c1 = -1.0/tmp1;
2050
2051 for (p1=1; p1<=NEB_Num_Images; p1++){
2052 for (i1=1; i1<=atomnum; i1++){
2053 for (j1=1; j1<=3; j1++){
2054 for (p2=1; p2<=NEB_Num_Images; p2++){
2055 for (i2=1; i2<=atomnum; i2++){
2056 for (j2=1; j2<=3; j2++){
2057
2058 dif_x1 = His_neb_atom_coordinates[0][p1][i1][j1]
2059 -His_neb_atom_coordinates[1][p1][i1][j1];
2060
2061 dif_x2 = His_neb_atom_coordinates[0][p2][i2][j2]
2062 -His_neb_atom_coordinates[1][p2][i2][j2];
2063
2064 InvHess[p1][i1][j1][p2][i2][j2] += c0*dif_x1*dif_x2;
2065
2066 dif_x1 = Vec0[p1][i1][j1];
2067 dif_x2 = His_neb_atom_coordinates[0][p2][i2][j2]
2068 -His_neb_atom_coordinates[1][p2][i2][j2];
2069
2070 InvHess[p1][i1][j1][p2][i2][j2] += c1*dif_x1*dif_x2;
2071
2072 dif_x1 = His_neb_atom_coordinates[0][p1][i1][j1]
2073 -His_neb_atom_coordinates[1][p1][i1][j1];
2074 dif_x2 = Vec0[p2][i2][j2];
2075
2076 InvHess[p1][i1][j1][p2][i2][j2] += c1*dif_x1*dif_x2;
2077
2078 }
2079 }
2080 }
2081 }
2082 }
2083 }
2084
2085 } /* else */
2086
2087 /* update atomic coordinates */
2088
2089 for (p1=1; p1<=NEB_Num_Images; p1++){
2090 for (i1=1; i1<=atomnum; i1++){
2091 for (j1=1; j1<=3; j1++){
2092
2093 sum1 = 0.0;
2094 for (p2=1; p2<=NEB_Num_Images; p2++){
2095 for (i2=1; i2<=atomnum; i2++){
2096 for (j2=1; j2<=3; j2++){
2097
2098 tmp = neb_atom_coordinates[p2][i2][j2+3];
2099 sum1 += InvHess[p1][i1][j1][p2][i2][j2]*tmp;
2100 }
2101 }
2102 }
2103
2104 neb_atom_coordinates[p1][i1][j1] -= sum1;
2105 }
2106 }
2107 }
2108
2109 /* In case of a too large updating, do a modest updating */
2110
2111 Max_Step = 0.0;
2112 for (p=1; p<=NEB_Num_Images; p++){
2113 for (j=1; j<=atomnum; j++){
2114 for (k=1; k<=3; k++){
2115
2116 dif_x = fabs(neb_atom_coordinates[p][j][k] - His_neb_atom_coordinates[0][p][j][k]);
2117 if (Max_Step<dif_x) Max_Step = dif_x;
2118 }
2119 }
2120 }
2121
2122 if (Criterion_Max_Step<Max_Step){
2123
2124 for (p=1; p<=NEB_Num_Images; p++){
2125
2126 for (j=1; j<=atomnum; j++){
2127 for (k=1; k<=3; k++){
2128
2129 dif_x = neb_atom_coordinates[p][j][k] - His_neb_atom_coordinates[0][p][j][k];
2130 neb_atom_coordinates[p][j][k] = His_neb_atom_coordinates[0][p][j][k]
2131 + dif_x/Max_Step*Criterion_Max_Step;
2132 }
2133 }
2134 }
2135
2136 Max_Step = Criterion_Max_Step;
2137 }
2138
2139 } /* else */
2140
2141 return Max_Step;
2142 }
2143
2144
2145
2146
2147
2148
2149
Steepest_Decent(int iter,double *** neb_atom_coordinates,double norm,double Max_Force)2150 void Steepest_Decent(int iter, double ***neb_atom_coordinates, double norm, double Max_Force)
2151 {
2152 int p,j,k;
2153 double SD_min,SD_max,SD_init;
2154
2155 /****************************************************
2156 set up SD_scaling
2157 ****************************************************/
2158
2159 SD_init = 0.1/(Max_Force+1.0e-10);
2160 SD_max = SD_init*20.0;
2161 SD_min = SD_init*0.05;
2162
2163 if (iter==1){
2164 SD_scaling = 0.1/(Max_Force+1.0e-10);
2165 }
2166 else {
2167
2168 if (Past_Norm[1]<norm && iter%4==1){
2169 SD_scaling = SD_scaling/4.0;
2170 }
2171 else if (Past_Norm[1]<Past_Norm[2] && norm<Past_Norm[1] && iter%4==1){
2172 SD_scaling = SD_scaling*1.2;
2173 }
2174
2175 Past_Norm[5] = Past_Norm[4];
2176 Past_Norm[4] = Past_Norm[3];
2177 Past_Norm[3] = Past_Norm[2];
2178 Past_Norm[2] = Past_Norm[1];
2179 Past_Norm[1] = norm;
2180 }
2181
2182 if (SD_max<SD_scaling) SD_scaling = SD_max;
2183 if (SD_scaling<SD_min) SD_scaling = SD_min;
2184
2185 /*************************************************************
2186 update coordinate by the simple steepest decent (SD) method
2187 *************************************************************/
2188
2189 for (p=1; p<=NEB_Num_Images; p++){
2190 for (j=1; j<=atomnum; j++){
2191 for (k=1; k<=3; k++){
2192 neb_atom_coordinates[p][j][k] -= SD_scaling*neb_atom_coordinates[p][j][k+3];
2193 }
2194 }
2195 }
2196 }
2197
2198
2199
2200
2201
Calc_NEB_Gradients(double *** neb_atom_coordinates)2202 void Calc_NEB_Gradients(double ***neb_atom_coordinates)
2203 {
2204 int i,j,k,p;
2205 double ***tangents;
2206 double dtmp1,dtmp2;
2207 double sum,prod1;
2208 double sconst;
2209 double Dvmax,Dvmin,vm,vi,vp;
2210 double dis1,dis2,dif;
2211
2212 /* allocation of arrays */
2213
2214 tangents = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2));
2215 for (i=0; i<(NEB_Num_Images+2); i++){
2216 tangents[i] = (double**)malloc(sizeof(double*)*(atomnum+1));
2217 for (j=0; j<(atomnum+1); j++){
2218 tangents[i][j] = (double*)malloc(sizeof(double)*4);
2219 for (k=0; k<4; k++) tangents[i][j][k] = 0.0;
2220 }
2221 }
2222
2223 /***********************************************************************
2224 calculate tangents based on Eqs. (8)-(11) in JCP 113, 9978 (2000)
2225 ***********************************************************************/
2226
2227 for (p=1; p<=NEB_Num_Images; p++){
2228
2229 vm = neb_atom_coordinates[p-1][0][0];
2230 vi = neb_atom_coordinates[p+0][0][0];
2231 vp = neb_atom_coordinates[p+1][0][0];
2232
2233 if ( vm<vi && vi<vp ){
2234
2235 for (j=1; j<=atomnum; j++){
2236 for (k=1; k<=3; k++){
2237 tangents[p][j][k] = neb_atom_coordinates[p+1][j][k] - neb_atom_coordinates[p][j][k];
2238 }
2239 }
2240 }
2241
2242 else if ( vp<vi && vi<vm ){
2243
2244 for (j=1; j<=atomnum; j++){
2245 for (k=1; k<=3; k++){
2246 tangents[p][j][k] = neb_atom_coordinates[p][j][k] - neb_atom_coordinates[p-1][j][k];
2247 }
2248 }
2249 }
2250
2251 else if ( (vp>=vi && vi<=vm) || (vp<vi && vi>vm) ){
2252
2253 Dvmax = MAX(fabs(vp-vi),fabs(vm-vi));
2254 Dvmin = MIN(fabs(vp-vi),fabs(vm-vi));
2255
2256 if (vm<=vp){
2257
2258 for (j=1; j<=atomnum; j++){
2259 for (k=1; k<=3; k++){
2260 tangents[p][j][k] = Dvmax*(neb_atom_coordinates[p+1][j][k] - neb_atom_coordinates[p ][j][k])
2261 + Dvmin*(neb_atom_coordinates[p ][j][k] - neb_atom_coordinates[p-1][j][k]);
2262 }
2263 }
2264 }
2265 else {
2266
2267 for (j=1; j<=atomnum; j++){
2268 for (k=1; k<=3; k++){
2269 tangents[p][j][k] = Dvmin*(neb_atom_coordinates[p+1][j][k] - neb_atom_coordinates[p ][j][k])
2270 + Dvmax*(neb_atom_coordinates[p ][j][k] - neb_atom_coordinates[p-1][j][k]);
2271 }
2272 }
2273 }
2274
2275 }
2276 else {
2277 printf("unknown situation\n");
2278 }
2279
2280 /* normalization of tangent */
2281
2282 sum = 0.0;
2283 for (j=1; j<=atomnum; j++){
2284 for (k=1; k<=3; k++){
2285 sum += tangents[p][j][k]*tangents[p][j][k];
2286 }
2287 }
2288
2289 sum = 1.0/sqrt(sum);
2290 for (j=1; j<=atomnum; j++){
2291 for (k=1; k<=3; k++){
2292 tangents[p][j][k] *= sum;
2293 }
2294 }
2295 }
2296
2297 /* calculate forces defined by the NEB method */
2298
2299 for (p=1; p<=NEB_Num_Images; p++){
2300
2301 /****************************************************************
2302 Perpendicular force
2303 ****************************************************************/
2304
2305 /* inner product of gradient[p] and tangents[p] */
2306
2307 prod1 = 0.0;
2308 for (j=1; j<=atomnum; j++){
2309 for (k=1; k<=3; k++){
2310 prod1 += neb_atom_coordinates[p][j][k+16]*tangents[p][j][k];
2311 }
2312 }
2313
2314 /* calculate the perpendicular gradient */
2315
2316 for (j=1; j<=atomnum; j++){
2317 for (k=1; k<=3; k++){
2318 neb_atom_coordinates[p][j][k+3] = neb_atom_coordinates[p][j][k+16] - prod1*tangents[p][j][k];
2319 }
2320 }
2321
2322 /****************************************************************
2323 Parallel force
2324 ****************************************************************/
2325
2326 /* force constant */
2327
2328 sconst = NEB_Spring_Const;
2329
2330 /* calculate the distance between Ri+1 and Ri */
2331
2332 dis1 = 0.0;
2333 for (j=1; j<=atomnum; j++){
2334 for (k=1; k<=3; k++){
2335 dif = neb_atom_coordinates[p+1][j][k] - neb_atom_coordinates[p][j][k];
2336 dis1 += dif*dif;
2337 }
2338 }
2339 dis1 = sqrt(dis1);
2340
2341 /* calculate the distance between Ri and Ri-1 */
2342
2343 dis2 = 0.0;
2344 for (j=1; j<=atomnum; j++){
2345 for (k=1; k<=3; k++){
2346 dif = neb_atom_coordinates[p][j][k] - neb_atom_coordinates[p-1][j][k];
2347 dis2 += dif*dif;
2348 }
2349 }
2350 dis2 = sqrt(dis2);
2351
2352 /* calculate the parallel gradient by springs, and add the contribution
2353 based on Eq. (12) in JCP 113, 9978 (2000) */
2354
2355 for (j=1; j<=atomnum; j++){
2356 for (k=1; k<=3; k++){
2357 neb_atom_coordinates[p][j][k+3] += -sconst*(dis1 - dis2)*tangents[p][j][k];
2358 }
2359 }
2360
2361 } /* p */
2362
2363 /****************************************************************
2364 introduce the contraints defined by user
2365 by setting gradients zero
2366 1: fixed
2367 0: relaxed
2368 ****************************************************************/
2369
2370 for (p=1; p<=NEB_Num_Images; p++){
2371 for (j=1; j<=atomnum; j++){
2372 for (k=1; k<=3; k++){
2373 if (atomFixedXYZ[j][k]==1){
2374 neb_atom_coordinates[p][j][k+3] = 0.0;
2375 }
2376 }
2377 }
2378 }
2379
2380 /* freeing of arrays */
2381
2382 for (i=0; i<(NEB_Num_Images+2); i++){
2383 for (j=0; j<(atomnum+1); j++){
2384 free(tangents[i][j]);
2385 }
2386 free(tangents[i]);
2387 }
2388 free(tangents);
2389
2390 }
2391
2392
2393
2394
2395
generate_input_files(char * file,int iter)2396 void generate_input_files(char *file, int iter)
2397 {
2398 int i,p,Gc_AN,c,n1,k;
2399 int restart_flag,fixed_flag;
2400 int unit_flag,level_flag;
2401 double c1,c2,c3;
2402 double tmpxyz[4];
2403 char st[800];
2404 char st1[800];
2405 char rm_operate[YOUSO10];
2406 char fname[YOUSO10];
2407 char fname1[YOUSO10];
2408 char fname2[YOUSO10];
2409 FILE *fp1,*fp2;
2410 char *tp;
2411 int numprocs,myid;
2412
2413 /* MPI */
2414 MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
2415 MPI_Comm_rank(MPI_COMM_WORLD,&myid);
2416
2417 /* get system.name */
2418
2419 if (input_open(file)==0){
2420 MPI_Finalize();
2421 exit(0);
2422 }
2423
2424 input_string("System.CurrrentDirectory",filepath,"./");
2425 input_string("System.Name",filename,"default");
2426
2427 input_close();
2428
2429 /* generate the input files */
2430
2431 if (myid==Host_ID){
2432
2433 for (p=0; p<=(NEB_Num_Images+1); p++){
2434
2435 /* initialize */
2436
2437 restart_flag = 0;
2438 fixed_flag = 0;
2439 unit_flag = 0;
2440 level_flag = 0;
2441
2442 /* the new input file */
2443
2444 sprintf(fname1,"%s_%i",file,p);
2445 fp1 = fopen(fname1,"w");
2446 fseek(fp1,0,SEEK_END);
2447
2448 /* the original input file */
2449
2450 fp2 = fopen(file,"r");
2451
2452 if (fp2!=NULL){
2453
2454 while (fgets(st,800,fp2)!=NULL){
2455
2456 string_tolower(st,st1);
2457
2458 /* find the specification of <atoms.speciesandcoordinates */
2459
2460 if (strncmp(st1,"<atoms.speciesandcoordinates",28)==0){
2461
2462 fprintf(fp1,"%s",st);
2463
2464 /* replace the atomic coordinates */
2465
2466 for (i=1; i<=atomnum; i++){
2467
2468 fgets(st,800,fp2);
2469 string_tolower(st,st1);
2470
2471 /* serial number */
2472 tp = strtok(st, " ");
2473 if (tp!=NULL) fprintf(fp1,"%4s",tp);
2474
2475 /* name of species */
2476 tp =strtok(NULL, " ");
2477 if (tp!=NULL) fprintf(fp1," %4s",tp);
2478
2479 /* x-coordinate */
2480 tp =strtok(NULL, " ");
2481 fprintf(fp1," %18.14f",neb_atom_coordinates[p][i][1]);
2482
2483 /* y-coordinate */
2484 tp =strtok(NULL, " ");
2485 fprintf(fp1," %18.14f",neb_atom_coordinates[p][i][2]);
2486
2487 /* z-coordinate */
2488 tp =strtok(NULL, " ");
2489 fprintf(fp1," %18.14f",neb_atom_coordinates[p][i][3]);
2490
2491 while (tp!=NULL){
2492 tp =strtok(NULL, " ");
2493 if (tp!=NULL) fprintf(fp1," %s",tp);
2494 }
2495
2496 }
2497 }
2498
2499 /* Atoms.SpeciesAndCoordinates.Unit */
2500
2501 else if (strncmp(st1,"atoms.speciesandcoordinates.unit",32)==0){
2502 fprintf(fp1,"Atoms.SpeciesAndCoordinates.Unit AU\n");
2503 unit_flag = 1;
2504 }
2505
2506 /* scf.restart for iter==1 */
2507
2508 else if (strncmp(st1,"scf.restart",11)==0 && iter==1 && Scf_RestartFromFile==0){
2509 fprintf(fp1,"scf.restart off\n");
2510 restart_flag = 1;
2511 }
2512
2513 else if (strncmp(st1,"scf.restart",11)==0 && iter==1 && Scf_RestartFromFile==1){
2514 fprintf(fp1,"scf.restart on\n");
2515 restart_flag = 1;
2516 }
2517
2518 /* scf.restart for iter!=1 */
2519
2520 else if (strncmp(st1,"scf.restart",11)==0 && iter!=1){
2521 fprintf(fp1,"scf.restart on\n");
2522 restart_flag = 1;
2523 }
2524
2525 /* scf.fixed.grid */
2526
2527 else if (strncmp(st1,"scf.fixed.grid",14)==0){
2528 fprintf(fp1,"%s",st);
2529 fixed_flag = 1;
2530 }
2531
2532 /* System.Name */
2533
2534 else if (strncmp(st1,"system.name",11)==0){
2535 fprintf(fp1,"System.Name %s_%i\n",filename,p);
2536 }
2537
2538 /* level.of.stdout */
2539
2540 else if (strncmp(st1,"level.of.stdout",15)==0 && p<=1 ){
2541 fprintf(fp1,"level.of.stdout 1\n");
2542 level_flag = 1;
2543 }
2544
2545 else if (strncmp(st1,"level.of.stdout",15)==0 && 1<p ){
2546 fprintf(fp1,"level.of.stdout 0\n");
2547 level_flag = 1;
2548 }
2549
2550 else{
2551 fprintf(fp1,"%s",st);
2552 }
2553 }
2554
2555 fclose(fp2);
2556 }
2557
2558 /* add the unit flag if it was not found. */
2559
2560 if (unit_flag==0){
2561 fprintf(fp1,"Atoms.SpeciesAndCoordinates.Unit AU\n");
2562 }
2563
2564 /* add the restart flag if it was not found. */
2565
2566 if (restart_flag==0 && iter==1){
2567 fprintf(fp1,"\n\nscf.restart off\n");
2568 }
2569 else if (restart_flag==0 && iter!=1){
2570 fprintf(fp1,"\n\nscf.restart on\n");
2571 }
2572
2573 /* add scf.fixed.grid if it was not found. */
2574
2575 if (fixed_flag==0 && 1<=p && p<=NEB_Num_Images && 1<iter){
2576 fprintf(fp1,"\n\nscf.fixed.grid %18.14f %18.14f %18.14f\n",
2577 All_Grid_Origin[p][1],All_Grid_Origin[p][2],All_Grid_Origin[p][3]);
2578 }
2579
2580 /* add level.of.stdout if it was not found. */
2581
2582 if (level_flag==0 && p<=1){
2583 fprintf(fp1,"\n\nlevel.of.stdout 1\n");
2584 }
2585 else if (level_flag==0 && 1<p){
2586 fprintf(fp1,"\n\nlevel.of.stdout 0\n");
2587 }
2588
2589 /* fclose */
2590 fclose(fp1);
2591
2592 } /* p */
2593 } /* if (myid==Host_ID) */
2594
2595 MPI_Barrier(MPI_COMM_WORLD);
2596
2597 }
2598
2599
2600
read_input(char * file)2601 void read_input(char *file)
2602 {
2603 int i,j,k,p;
2604 int SpinP_switch;
2605 int po,itmp;
2606 int unitvector_unit;
2607 int unitvectors_flag;
2608 int itmp1,itmp2,itmp3,itmp4;
2609 int myid,numprocs;
2610 double tv[4][4];
2611 double dtmp1,dtmp2,dtmp3,dtmp4,dtmp5;
2612 double dtmp6,dtmp7,dtmp8,dtmp9;
2613 double tmpx,tmpy,tmpz,x,y,z;
2614 double dif,sx,sy,sz;
2615 double xc0,yc0,zc0,xc1,yc1,zc1;
2616 double xc,yc,zc;
2617 char Species[100];
2618 char OrbPol[100];
2619 char buf[MAXBUF];
2620 char keyword[YOUSO10];
2621 int i_vec[40];
2622 char *s_vec[40];
2623 FILE *fp;
2624
2625 MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
2626 MPI_Comm_rank(MPI_COMM_WORLD,&myid);
2627
2628 /* open the input file */
2629
2630 if (input_open(file)==0){
2631 MPI_Finalize();
2632 exit(0);
2633 }
2634
2635 input_string("System.CurrrentDirectory",filepath,"./");
2636 input_string("System.Name",filename,"default");
2637 sprintf(system_name,"%s",filename);
2638
2639 /**************************************************************
2640 NEB_Num_Images is the number of images (p) excluding
2641 the end points (0 and p+1).
2642 *************************************************************/
2643
2644 input_int("MD.NEB.Number.Images",&NEB_Num_Images,10);
2645 input_int("MD.NEB.Spring.Const",&NEB_Spring_Const,0.1); /* hartree/bohr^2 */
2646 input_int("MD.maxIter",&MD_IterNumber,1);
2647
2648 input_int("MD.NEB.Parallel.Number",&PN,0);
2649
2650 s_vec[0]="Off"; s_vec[1]="On"; s_vec[2]="NC";
2651 i_vec[0]=0 ; i_vec[1]=1 ; i_vec[2]=3;
2652 input_string2int("scf.SpinPolarization", &SpinP_switch, 3, s_vec,i_vec);
2653
2654 input_int("Atoms.Number",&atomnum,0);
2655 input_int("Species.Number",&SpeciesNum,0);
2656
2657 input_int("MD.Opt.DIIS.History",&M_GDIIS_HISTORY,3);
2658
2659 /* allocate arrays */
2660 allocate_arrays();
2661
2662 /*********************************************************
2663 read the unit cell
2664 *********************************************************/
2665
2666 s_vec[0]="Ang"; s_vec[1]="AU";
2667 i_vec[0]=0; i_vec[1]=1;
2668 input_string2int("Atoms.UnitVectors.Unit",&unitvector_unit,2,s_vec,i_vec);
2669
2670 unitvectors_flag = 0;
2671
2672 if (fp=input_find("<Atoms.Unitvectors")) {
2673
2674 unitvectors_flag = 1;
2675
2676 for (i=1; i<=3; i++){
2677 fscanf(fp,"%lf %lf %lf",&tv[i][1],&tv[i][2],&tv[i][3]);
2678 }
2679 if ( ! input_last("Atoms.Unitvectors>") ) {
2680 /* format error */
2681 if (myid==Host_ID){
2682 printf("Format error for Atoms.Unitvectors\n");
2683 }
2684 MPI_Finalize();
2685 exit(0);
2686 }
2687
2688 /* Ang to AU */
2689 if (unitvector_unit==0){
2690 for (i=1; i<=3; i++){
2691 tv[i][1] = tv[i][1]/BohrR;
2692 tv[i][2] = tv[i][2]/BohrR;
2693 tv[i][3] = tv[i][3]/BohrR;
2694 }
2695 }
2696 }
2697
2698 if (unitvectors_flag==0){
2699 if (myid==Host_ID){
2700 printf("A common unit cell for two terminal structures has to be specified.\n");fflush(stdout);
2701 }
2702 MPI_Finalize();
2703 exit(0);
2704 }
2705
2706 /**************************************************************
2707 read atomic coodinates given by Atoms.SpeciesAndCoordinates.
2708 **************************************************************/
2709
2710 s_vec[0]="Ang"; s_vec[1]="AU"; s_vec[2]="FRAC";
2711 i_vec[0]= 0; i_vec[1]= 1; i_vec[2]= 2;
2712 input_string2int("Atoms.SpeciesAndCoordinates.Unit",
2713 &coordinates_unit,3,s_vec,i_vec);
2714
2715 if (fp=input_find("<Atoms.SpeciesAndCoordinates") ) {
2716
2717 for (i=1; i<=atomnum; i++){
2718
2719 fgets(buf,MAXBUF,fp);
2720
2721 /* spin non-collinear */
2722 if (SpinP_switch==3){
2723
2724 /*******************************************************
2725 (1) spin non-collinear
2726 *******************************************************/
2727
2728 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %s",
2729 &j, Species,
2730 &neb_atom_coordinates[0][i][1],
2731 &neb_atom_coordinates[0][i][2],
2732 &neb_atom_coordinates[0][i][3],
2733 &dtmp1,&dtmp2,
2734 &dtmp3,&dtmp4,
2735 &dtmp5,&dtmp6,
2736 &itmp2,
2737 OrbPol );
2738
2739 }
2740
2741 /**************************************************
2742 (2) spin collinear
2743 **************************************************/
2744
2745 else{
2746
2747 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %s",
2748 &j, Species,
2749 &neb_atom_coordinates[0][i][1],
2750 &neb_atom_coordinates[0][i][2],
2751 &neb_atom_coordinates[0][i][3],
2752 &dtmp1,&dtmp2, OrbPol );
2753 }
2754
2755 if (i!=j){
2756 if (myid==Host_ID){
2757 printf("Format error of the sequential number %i in <Atoms.SpeciesAndCoordinates\n",j);
2758 }
2759 MPI_Finalize();
2760 exit(0);
2761 }
2762
2763 }
2764
2765 ungetc('\n',fp);
2766 if (!input_last("Atoms.SpeciesAndCoordinates>")) {
2767 /* format error */
2768 if (myid==Host_ID){
2769 printf("Format error for Atoms.SpeciesAndCoordinates\n");fflush(stdout);
2770 }
2771
2772 MPI_Finalize();
2773 exit(0);
2774 }
2775
2776 /* Ang to AU */
2777
2778 if (coordinates_unit==0){
2779 for (i=1; i<=atomnum; i++){
2780 neb_atom_coordinates[0][i][1] /= BohrR;
2781 neb_atom_coordinates[0][i][2] /= BohrR;
2782 neb_atom_coordinates[0][i][3] /= BohrR;
2783 }
2784 }
2785 }
2786
2787 /* FRAC to AU */
2788 if (coordinates_unit==2){
2789
2790 /* The fractional coordinates should be kept within 0 to 1. */
2791
2792 for (i=1; i<=atomnum; i++){
2793 for (k=1; k<=3; k++){
2794
2795 itmp = (int)neb_atom_coordinates[0][i][k];
2796
2797 if (1.0<neb_atom_coordinates[0][i][k]){
2798
2799 /* ended by T.Ohwaki */
2800
2801 neb_atom_coordinates[0][i][k] = neb_atom_coordinates[0][i][k] - (double)itmp;
2802
2803 if (myid==Host_ID){
2804 if (k==1) printf("The fractional coordinate of a-axis for atom %d was translated within the range (0 to 1).\n",i);
2805 if (k==2) printf("The fractional coordinate of b-axis for atom %d was translated within the range (0 to 1).\n",i);
2806 if (k==3) printf("The fractional coordinate of c-axis for atom %d was translated within the range (0 to 1).\n",i);
2807 }
2808 }
2809 else if (neb_atom_coordinates[0][i][k]<0.0){
2810
2811 neb_atom_coordinates[0][i][k] = neb_atom_coordinates[0][i][k] + (double)(abs(itmp)+1);
2812
2813 if (myid==Host_ID){
2814 if (k==1) printf("The fractional coordinate of a-axis for atom %d was translated within the range (0 to 1).\n",i);
2815 if (k==2) printf("The fractional coordinate of b-axis for atom %d was translated within the range (0 to 1).\n",i);
2816 if (k==3) printf("The fractional coordinate of c-axis for atom %d was translated within the range (0 to 1).\n",i);
2817 }
2818 }
2819 }
2820 }
2821
2822 /* calculation of xyz-coordinate in A.U. The grid origin is zero. */
2823
2824 tmpx = 0.0;
2825 tmpy = 0.0;
2826 tmpz = 0.0;
2827
2828 for (i=1; i<=atomnum; i++){
2829
2830 x = neb_atom_coordinates[0][i][1]*tv[1][1]
2831 + neb_atom_coordinates[0][i][2]*tv[2][1]
2832 + neb_atom_coordinates[0][i][3]*tv[3][1] + tmpx;
2833
2834 y = neb_atom_coordinates[0][i][1]*tv[1][2]
2835 + neb_atom_coordinates[0][i][2]*tv[2][2]
2836 + neb_atom_coordinates[0][i][3]*tv[3][2] + tmpy;
2837
2838 z = neb_atom_coordinates[0][i][1]*tv[1][3]
2839 + neb_atom_coordinates[0][i][2]*tv[2][3]
2840 + neb_atom_coordinates[0][i][3]*tv[3][3] + tmpz;
2841
2842 neb_atom_coordinates[0][i][1] = x;
2843 neb_atom_coordinates[0][i][2] = y;
2844 neb_atom_coordinates[0][i][3] = z;
2845 }
2846 }
2847
2848 /*****************************************************************
2849 read atomic coodinates given by NEB.Atoms.SpeciesAndCoordinates
2850 *****************************************************************/
2851
2852 if (fp=input_find("<NEB.Atoms.SpeciesAndCoordinates") ) {
2853
2854 for (i=1; i<=atomnum; i++){
2855
2856 fgets(buf,MAXBUF,fp);
2857
2858 /* spin non-collinear */
2859 if (SpinP_switch==3){
2860
2861 /*******************************************************
2862 (1) spin non-collinear
2863 *******************************************************/
2864
2865 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %d %s",
2866 &j, Species,
2867 &neb_atom_coordinates[NEB_Num_Images+1][i][1],
2868 &neb_atom_coordinates[NEB_Num_Images+1][i][2],
2869 &neb_atom_coordinates[NEB_Num_Images+1][i][3],
2870 &dtmp1,&dtmp2,
2871 &dtmp3,&dtmp4,
2872 &dtmp5,&dtmp6,
2873 &itmp2,
2874 OrbPol );
2875
2876 }
2877
2878 /**************************************************
2879 (2) spin collinear
2880 **************************************************/
2881
2882 else{
2883
2884 sscanf(buf,"%i %s %lf %lf %lf %lf %lf %s",
2885 &j, Species,
2886 &neb_atom_coordinates[NEB_Num_Images+1][i][1],
2887 &neb_atom_coordinates[NEB_Num_Images+1][i][2],
2888 &neb_atom_coordinates[NEB_Num_Images+1][i][3],
2889 &dtmp1,&dtmp2, OrbPol );
2890 }
2891
2892 if (i!=j){
2893 if (myid==Host_ID){
2894 printf("Format error of the sequential number %i in <NEB.Atoms.SpeciesAndCoordinates\n",j);
2895 }
2896 MPI_Finalize();
2897 exit(0);
2898 }
2899
2900 }
2901
2902 ungetc('\n',fp);
2903 if (!input_last("NEB.Atoms.SpeciesAndCoordinates>")) {
2904 /* format error */
2905 if (myid==Host_ID){
2906 printf("Format error for NEB.Atoms.SpeciesAndCoordinates\n");
2907 }
2908 MPI_Finalize();
2909 exit(0);
2910 }
2911
2912 /* Ang to AU */
2913
2914 if (coordinates_unit==0){
2915 for (i=1; i<=atomnum; i++){
2916 neb_atom_coordinates[NEB_Num_Images+1][i][1] /= BohrR;
2917 neb_atom_coordinates[NEB_Num_Images+1][i][2] /= BohrR;
2918 neb_atom_coordinates[NEB_Num_Images+1][i][3] /= BohrR;
2919 }
2920 }
2921 }
2922
2923 /* FRAC to AU */
2924 if (coordinates_unit==2){
2925
2926 /* The fractional coordinates should be kept within 0 to 1. */
2927
2928 for (i=1; i<=atomnum; i++){
2929 for (k=1; k<=3; k++){
2930
2931 itmp = (int)neb_atom_coordinates[NEB_Num_Images+1][i][k];
2932
2933 if (1.0<neb_atom_coordinates[NEB_Num_Images+1][i][k]){
2934
2935 /* ended by T.Ohwaki */
2936
2937 neb_atom_coordinates[NEB_Num_Images+1][i][k] = neb_atom_coordinates[NEB_Num_Images+1][i][k] - (double)itmp;
2938
2939 if (myid==Host_ID){
2940 if (k==1) printf("The fractional coordinate of a-axis for atom %d was translated within the range (0 to 1).\n",i);
2941 if (k==2) printf("The fractional coordinate of b-axis for atom %d was translated within the range (0 to 1).\n",i);
2942 if (k==3) printf("The fractional coordinate of c-axis for atom %d was translated within the range (0 to 1).\n",i);
2943 }
2944 }
2945 else if (neb_atom_coordinates[NEB_Num_Images+1][i][k]<0.0){
2946
2947 neb_atom_coordinates[NEB_Num_Images+1][i][k] = neb_atom_coordinates[NEB_Num_Images+1][i][k] + (double)(abs(itmp)+1);
2948
2949 if (myid==Host_ID){
2950 if (k==1) printf("The fractional coordinate of a-axis for atom %d was translated within the range (0 to 1).\n",i);
2951 if (k==2) printf("The fractional coordinate of b-axis for atom %d was translated within the range (0 to 1).\n",i);
2952 if (k==3) printf("The fractional coordinate of c-axis for atom %d was translated within the range (0 to 1).\n",i);
2953 }
2954 }
2955 }
2956 }
2957
2958 /* calculation of xyz-coordinate in A.U. The grid origin is zero. */
2959
2960 tmpx = 0.0;
2961 tmpy = 0.0;
2962 tmpz = 0.0;
2963
2964 for (i=1; i<=atomnum; i++){
2965
2966 x = neb_atom_coordinates[NEB_Num_Images+1][i][1]*tv[1][1]
2967 + neb_atom_coordinates[NEB_Num_Images+1][i][2]*tv[2][1]
2968 + neb_atom_coordinates[NEB_Num_Images+1][i][3]*tv[3][1] + tmpx;
2969
2970 y = neb_atom_coordinates[NEB_Num_Images+1][i][1]*tv[1][2]
2971 + neb_atom_coordinates[NEB_Num_Images+1][i][2]*tv[2][2]
2972 + neb_atom_coordinates[NEB_Num_Images+1][i][3]*tv[3][2] + tmpy;
2973
2974 z = neb_atom_coordinates[NEB_Num_Images+1][i][1]*tv[1][3]
2975 + neb_atom_coordinates[NEB_Num_Images+1][i][2]*tv[2][3]
2976 + neb_atom_coordinates[NEB_Num_Images+1][i][3]*tv[3][3] + tmpz;
2977
2978 neb_atom_coordinates[NEB_Num_Images+1][i][1] = x;
2979 neb_atom_coordinates[NEB_Num_Images+1][i][2] = y;
2980 neb_atom_coordinates[NEB_Num_Images+1][i][3] = z;
2981 }
2982 }
2983
2984 /*****************************************************************
2985 if scf.restart==on, then read the coordinates of images
2986 and read total energies of the two terminal structures.
2987 *****************************************************************/
2988
2989 input_logical("scf.restart",&Scf_RestartFromFile, 0);
2990
2991 if (Scf_RestartFromFile){
2992
2993 /* read the structures */
2994
2995 for (p=1; p<=NEB_Num_Images; p++){
2996
2997 sprintf(keyword,"<NEB%d.Atoms.SpeciesAndCoordinates",p);
2998
2999 if (fp=input_find(keyword) ) {
3000
3001 for (i=1; i<=atomnum; i++){
3002
3003 fscanf(fp,"%i %s %lf %lf %lf %lf %lf",
3004 &j, Species,
3005 &neb_atom_coordinates[p][i][1],
3006 &neb_atom_coordinates[p][i][2],
3007 &neb_atom_coordinates[p][i][3],
3008 &dtmp1,&dtmp2);
3009 }
3010
3011 sprintf(keyword,"NEB%d.Atoms.SpeciesAndCoordinates>",p);
3012
3013 if (!input_last(keyword)) {
3014 /* format error */
3015 if (myid==Host_ID){
3016 printf("Format error for NEB%d.Atoms.SpeciesAndCoordinates\n",p);
3017 }
3018 MPI_Finalize();
3019 exit(0);
3020 }
3021
3022 }
3023 }
3024
3025 }
3026
3027 else {
3028
3029 /*****************************************************************
3030 Making the centroid of two terminal coordinates equivalent.
3031 This treatment allows us to hold the equivalence of the centroid
3032 of all the images. Thus, Grid_Origin also becomes equivalent to
3033 eath other.
3034 *****************************************************************/
3035
3036 xc0 = 0.0;
3037 yc0 = 0.0;
3038 zc0 = 0.0;
3039
3040 xc1 = 0.0;
3041 yc1 = 0.0;
3042 zc1 = 0.0;
3043
3044 for (i=1; i<=atomnum; i++){
3045
3046 xc0 += neb_atom_coordinates[0][i][1];
3047 yc0 += neb_atom_coordinates[0][i][2];
3048 zc0 += neb_atom_coordinates[0][i][3];
3049
3050 xc1 += neb_atom_coordinates[NEB_Num_Images+1][i][1];
3051 yc1 += neb_atom_coordinates[NEB_Num_Images+1][i][2];
3052 zc1 += neb_atom_coordinates[NEB_Num_Images+1][i][3];
3053 }
3054
3055 xc0 = xc0/(double)atomnum;
3056 yc0 = yc0/(double)atomnum;
3057 zc0 = zc0/(double)atomnum;
3058
3059 xc1 = xc1/(double)atomnum;
3060 yc1 = yc1/(double)atomnum;
3061 zc1 = zc1/(double)atomnum;
3062
3063 xc = 0.5*(xc0 + xc1);
3064 yc = 0.5*(yc0 + yc1);
3065 zc = 0.5*(zc0 + zc1);
3066
3067 for (i=1; i<=atomnum; i++){
3068
3069 neb_atom_coordinates[0][i][1] += -xc0 + xc;
3070 neb_atom_coordinates[0][i][2] += -yc0 + yc;
3071 neb_atom_coordinates[0][i][3] += -zc0 + zc;
3072
3073 neb_atom_coordinates[NEB_Num_Images+1][i][1] += -xc1 + xc;
3074 neb_atom_coordinates[NEB_Num_Images+1][i][2] += -yc1 + yc;
3075 neb_atom_coordinates[NEB_Num_Images+1][i][3] += -zc1 + zc;
3076 }
3077
3078 /*****************************************************************
3079 generate atomic coordinates for images
3080 *****************************************************************/
3081
3082 for (p=1; p<=NEB_Num_Images; p++){
3083 for (i=1; i<=atomnum; i++){
3084 for (j=1; j<=3; j++){
3085 dif = neb_atom_coordinates[NEB_Num_Images+1][i][j]-neb_atom_coordinates[0][i][j];
3086 neb_atom_coordinates[p][i][j] = neb_atom_coordinates[0][i][j]
3087 + dif*(double)p/(double)(NEB_Num_Images+1);
3088 }
3089 }
3090 }
3091 }
3092
3093 /*****************************************************************
3094 read atomFixedXYZ
3095 set fixed atomic position in geometry optimization
3096 and MD:
3097
3098 1: fixed
3099 0: relaxed
3100 *****************************************************************/
3101
3102 if (fp=input_find("<MD.Fixed.XYZ")) {
3103
3104 for (i=1; i<=atomnum; i++){
3105 fscanf(fp,"%d %d %d %d",
3106 &j,&atomFixedXYZ[i][1],&atomFixedXYZ[i][2],&atomFixedXYZ[i][3]);
3107 }
3108
3109 if ( ! input_last("MD.Fixed.XYZ>") ) {
3110 /* format error */
3111
3112 if (myid==Host_ID){
3113 printf("Format error for MD.Fixed.XYZ\n");
3114 }
3115 MPI_Finalize();
3116 exit(0);
3117 }
3118 }
3119
3120 /* close the input file */
3121
3122 input_close();
3123
3124 }
3125
3126
allocate_arrays()3127 void allocate_arrays()
3128 {
3129 int i,j,k,m;
3130 int p1,i1,j1,p2,i2,j2;
3131
3132 neb_atom_coordinates = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2));
3133 for (i=0; i<(NEB_Num_Images+2); i++){
3134 neb_atom_coordinates[i] = (double**)malloc(sizeof(double*)*(atomnum+1));
3135 for (j=0; j<(atomnum+1); j++){
3136 neb_atom_coordinates[i][j] = (double*)malloc(sizeof(double)*20);
3137 for (k=0; k<20; k++) neb_atom_coordinates[i][j][k] = 0.0;
3138 }
3139 }
3140
3141 tmp_neb_atom_coordinates = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2));
3142 for (i=0; i<(NEB_Num_Images+2); i++){
3143 tmp_neb_atom_coordinates[i] = (double**)malloc(sizeof(double*)*(atomnum+1));
3144 for (j=0; j<(atomnum+1); j++){
3145 tmp_neb_atom_coordinates[i][j] = (double*)malloc(sizeof(double)*20);
3146 for (k=0; k<20; k++) neb_atom_coordinates[i][j][k] = 0.0;
3147 }
3148 }
3149
3150 His_neb_atom_coordinates = (double****)malloc(sizeof(double***)*(M_GDIIS_HISTORY+2));
3151 for (m=0; m<(M_GDIIS_HISTORY+2); m++){
3152 His_neb_atom_coordinates[m] = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2));
3153 for (i=0; i<(NEB_Num_Images+2); i++){
3154 His_neb_atom_coordinates[m][i] = (double**)malloc(sizeof(double*)*(atomnum+1));
3155 for (j=0; j<(atomnum+1); j++){
3156 His_neb_atom_coordinates[m][i][j] = (double*)malloc(sizeof(double)*20);
3157 for (k=0; k<20; k++) His_neb_atom_coordinates[m][i][j][k] = 0.0;
3158 }
3159 }
3160 }
3161
3162 All_Grid_Origin = (double**)malloc(sizeof(double*)*(NEB_Num_Images+2));
3163 for (i=0; i<(NEB_Num_Images+2); i++){
3164 All_Grid_Origin[i] = (double*)malloc(sizeof(double)*4);
3165 for (j=0; j<4; j++) All_Grid_Origin[i][j] = 0.0;
3166 }
3167
3168 Tmp_Grid_Origin = (double**)malloc(sizeof(double*)*(NEB_Num_Images+2));
3169 for (i=0; i<(NEB_Num_Images+2); i++){
3170 Tmp_Grid_Origin[i] = (double*)malloc(sizeof(double)*4);
3171 for (j=0; j<4; j++) Tmp_Grid_Origin[i][j] = 0.0;
3172 }
3173
3174 WhatSpecies_NEB = (int*)malloc(sizeof(int)*(atomnum+1));
3175 Spe_WhatAtom_NEB = (int*)malloc(sizeof(int)*SpeciesNum);
3176
3177 SpeName_NEB = (char**)malloc(sizeof(char*)*SpeciesNum);
3178 for (i=0; i<SpeciesNum; i++){
3179 SpeName_NEB[i] = (char*)malloc(sizeof(char)*YOUSO10);
3180 }
3181
3182 InvHess = (double******)malloc(sizeof(double*****)*(NEB_Num_Images+2));
3183 for (p1=0; p1<(NEB_Num_Images+2); p1++){
3184 InvHess[p1] = (double*****)malloc(sizeof(double****)*(atomnum+1));
3185 for (i1=0; i1<(atomnum+1); i1++){
3186 InvHess[p1][i1] = (double****)malloc(sizeof(double***)*4);
3187 for (j1=0; j1<4; j1++){
3188 InvHess[p1][i1][j1] = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2));
3189 for (p2=0; p2<(NEB_Num_Images+2); p2++){
3190 InvHess[p1][i1][j1][p2] = (double**)malloc(sizeof(double*)*(atomnum+1));
3191 for (i2=0; i2<(atomnum+1); i2++){
3192 InvHess[p1][i1][j1][p2][i2] = (double*)malloc(sizeof(double)*4);
3193 for (j2=0; j2<4; j2++){
3194 InvHess[p1][i1][j1][p2][i2][j2] = 0.0;
3195 }
3196 }
3197 }
3198 }
3199 }
3200 }
3201
3202 /* initialize InvHess */
3203 for (p1=1; p1<=NEB_Num_Images; p1++){
3204 for (i1=1; i1<=atomnum; i1++){
3205 for (j1=1; j1<=3; j1++){
3206 InvHess[p1][i1][j1][p1][i1][j1] = 1.0;
3207 }
3208 }
3209 }
3210
3211 Vec0 = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2));
3212 for (p1=0; p1<(NEB_Num_Images+2); p1++){
3213 Vec0[p1] = (double**)malloc(sizeof(double*)*(atomnum+1));
3214 for (i1=0; i1<(atomnum+1); i1++){
3215 Vec0[p1][i1] = (double*)malloc(sizeof(double)*4);
3216 for (j1=0; j1<4; j1++) Vec0[p1][i1][j1] = 0.0;
3217 }
3218 }
3219
3220 atomFixedXYZ = (int**)malloc(sizeof(int*)*(atomnum+1));
3221 for(i=0; i<=atomnum; i++){
3222 atomFixedXYZ[i] = (int*)malloc(sizeof(int)*4);
3223 /* default='relaxed' */
3224 atomFixedXYZ[i][1] = 0;
3225 atomFixedXYZ[i][2] = 0;
3226 atomFixedXYZ[i][3] = 0;
3227 }
3228
3229 }
3230
3231
3232
free_arrays()3233 void free_arrays()
3234 {
3235 int m,i,j,k;
3236 int p1,i1,j1,p2,i2,j2;
3237
3238 for (i=0; i<(NEB_Num_Images+2); i++){
3239 for (j=0; j<(atomnum+1); j++){
3240 free(neb_atom_coordinates[i][j]);
3241 }
3242 free(neb_atom_coordinates[i]);
3243 }
3244 free(neb_atom_coordinates);
3245
3246 for (i=0; i<(NEB_Num_Images+2); i++){
3247 for (j=0; j<(atomnum+1); j++){
3248 free(tmp_neb_atom_coordinates[i][j]);
3249 }
3250 free(tmp_neb_atom_coordinates[i]);
3251 }
3252 free(tmp_neb_atom_coordinates);
3253
3254 for (m=0; m<(M_GDIIS_HISTORY+2); m++){
3255 for (i=0; i<(NEB_Num_Images+2); i++){
3256 for (j=0; j<(atomnum+1); j++){
3257 free(His_neb_atom_coordinates[m][i][j]);
3258 }
3259 free(His_neb_atom_coordinates[m][i]);
3260 }
3261 free(His_neb_atom_coordinates[m]);
3262 }
3263 free(His_neb_atom_coordinates);
3264
3265 for (i=0; i<(NEB_Num_Images+2); i++){
3266 free(All_Grid_Origin[i]);
3267 }
3268 free(All_Grid_Origin);
3269
3270 for (i=0; i<(NEB_Num_Images+2); i++){
3271 free(Tmp_Grid_Origin[i]);
3272 }
3273 free(Tmp_Grid_Origin);
3274
3275 free(WhatSpecies_NEB);
3276 free(Spe_WhatAtom_NEB);
3277
3278 for (i=0; i<SpeciesNum; i++){
3279 free(SpeName_NEB[i]);
3280 }
3281 free(SpeName_NEB);
3282
3283 for (p1=0; p1<(NEB_Num_Images+2); p1++){
3284 for (i1=0; i1<(atomnum+1); i1++){
3285 for (j1=0; j1<4; j1++){
3286 for (p2=0; p2<(NEB_Num_Images+2); p2++){
3287 for (i2=0; i2<(atomnum+1); i2++){
3288 free(InvHess[p1][i1][j1][p2][i2]);
3289 }
3290 free(InvHess[p1][i1][j1][p2]);
3291 }
3292 free(InvHess[p1][i1][j1]);
3293 }
3294 free(InvHess[p1][i1]);
3295 }
3296 free(InvHess[p1]);
3297 }
3298 free(InvHess);
3299
3300 for (p1=0; p1<(NEB_Num_Images+2); p1++){
3301 for (i1=0; i1<(atomnum+1); i1++){
3302 free(Vec0[p1][i1]);
3303 }
3304 free(Vec0[p1]);
3305 }
3306 free(Vec0);
3307
3308 for(i=0; i<=atomnum; i++){
3309 free(atomFixedXYZ[i]);
3310 }
3311 free(atomFixedXYZ);
3312
3313 }
3314