/********************************************************************** neb.c: neb is a program to perform the nudged elastic band (NEB) method for finding a minimum energy path (MEP) connecting two structures, which is based on JCP 113, 9978 (2000). Log of neb.c: Apr./06/2011 Released by T. Ozaki Feb./01/2013 Modified by Y. Kubota (supervised by Prof. F. Ishii) Feb./17/2015 Modified by T. Ozaki ***********************************************************************/ #include #include #include #include #include #include #include #include #include #include "openmx_common.h" #include "Inputtools.h" #include "lapack_prototypes.h" #define MAXBUF 1024 #define Criterion_Max_Step 0.10 #ifdef MAX #undef MAX #endif #define MAX(a,b) (((a)>(b))? (a):(b)) #ifdef MIN #undef MIN #endif #define MIN(a,b) (((a)<(b))? (a):(b)) static void allocate_arrays(); static void free_arrays(); static void read_input(char *file); static void generate_input_files(char *file, int iter); static void Calc_NEB_Gradients(double ***neb_atom_coordinates); static void Update_Coordinates(int iter, double ***neb_atom_coordinates); static void Make_XYZ_File(char fname[YOUSO10], double ***neb_atom_coordinates); static void Generate_Restart_File(char fname[YOUSO10], double ***neb_atom_coordinates); static void Steepest_Decent(int iter, double ***neb_atom_coordinates, double norm, double Max_Force); static double DIIS_BFGS(int iter, double ***neb_atom_coordinates, double norm, double Max_Force); double **All_Grid_Origin; double **Tmp_Grid_Origin; double ***neb_atom_coordinates; double ***tmp_neb_atom_coordinates; double ****His_neb_atom_coordinates; double ******InvHess; double ***Vec0; int **atomFixedXYZ; static char fileXYZ[YOUSO10] = ".neb.xyz"; char system_name[YOUSO10]; int *WhatSpecies_NEB; int *Spe_WhatAtom_NEB; char **SpeName_NEB; int PN; void neb(int argc, char *argv[]) { int iter,index_images; int po,i,j,k,p,h; int myid,numprocs; int myid1,numprocs1; int myid2,numprocs2; int myid3,numprocs3; int parallel1_flag; int myworld1,ID0,ID1; int Num_Comm_World1; int *NPROCS_ID1; int *Comm_World1; int *NPROCS_WD1; int *Comm_World_StartID1; int parallel2_flag; int myworld2,myworld3; int Num_Comm_World2,Num_Comm_World3; int *NPROCS_ID2; int *NPROCS_ID3; int *Comm_World2; int *Comm_World3; int *NPROCS_WD2; int *NPROCS_WD3; int *Comm_World_StartID2; int *Comm_WOrld_StartID3; char fname_original[YOUSO10]; char fname1[YOUSO10]; MPI_Comm *MPI_CommWD1; MPI_Comm *MPI_CommWD2; MPI_Comm *MPI_CommWD3; char file_neb_utot[YOUSO10]; double TStime,TEtime,a0time,a1time,f0time; double b0time,b1time,c0time,c1time,flatstime,flatetime,sumtime=0.0; FILE *fp; dtime(&TStime); MPI_Comm_size(MPI_COMM_WORLD1,&numprocs); MPI_Comm_rank(MPI_COMM_WORLD1,&myid); /* check argv */ if (argc==1){ printf("\nCould not find an input file.\n\n"); MPI_Finalize(); exit(0); } /**************************************************** show a message by the NEB code ****************************************************/ if (myid==Host_ID){ printf("\n*******************************************************\n"); printf("*******************************************************\n"); printf(" Welcome to the NEB extension of OpenMX \n"); printf(" Copyright (C), 2002-2011, T. Ozaki \n"); printf(" OpenMX comes with ABSOLUTELY NO WARRANTY. \n"); printf(" This is free software, and you are welcome to \n"); printf(" redistribute it under the constitution of the GNU-GPL.\n"); printf("*******************************************************\n"); printf("*******************************************************\n\n"); } /**************************************************** read the input file ****************************************************/ read_input(argv[1]); fnjoint(filepath,filename,fileXYZ); sprintf(fname_original,"%s",argv[1]); /**************************************************** compare PN with numprocs and NEB_Num_Images ****************************************************/ if ((PN > numprocs)||(PN > NEB_Num_Images)){ PN=0; if (myid==Host_ID){ printf("The keyword MD.NEB.Parallel.Number will be ignored if the following conditions are satisfied:\n"); printf("MD.NEB.Parallel.Number is equal to or smaller than the number of MPI processes.\n"); printf("MD.NEB.Parallel.Number is equal to or smaller than the number of MD.NEB.Number.Images.\n"); } } /******************************************************************************* If MD.NEB.Parallel.Number is not specifed, *******************************************************************************/ if (PN==0){ /**************************************************** Two level parallelization is performed. Outer loop: images Inner loop: calculation in each image ****************************************************/ if ( NEB_Num_Images<=numprocs || numprocs==1 ){ /**************************************************** allocate processes to each image in the World1 ****************************************************/ Num_Comm_World1 = NEB_Num_Images; if ( Num_Comm_World1<=numprocs ){ parallel1_flag = 1; NPROCS_ID1 = (int*)malloc(sizeof(int)*numprocs); Comm_World1 = (int*)malloc(sizeof(int)*numprocs); NPROCS_WD1 = (int*)malloc(sizeof(int)*Num_Comm_World1); Comm_World_StartID1 = (int*)malloc(sizeof(int)*Num_Comm_World1); MPI_CommWD1 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World1); Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World1, &myworld1, MPI_CommWD1, NPROCS_ID1, Comm_World1, NPROCS_WD1, Comm_World_StartID1); MPI_Comm_size(MPI_CommWD1[myworld1],&numprocs1); MPI_Comm_rank(MPI_CommWD1[myworld1],&myid1); } else { parallel1_flag = 0; myworld1 = 0; } /**************************************************** allocate processes to each image in the World2 ****************************************************/ Num_Comm_World2 = 2; if ( Num_Comm_World2<=numprocs ){ parallel2_flag = 1; NPROCS_ID2 = (int*)malloc(sizeof(int)*numprocs); Comm_World2 = (int*)malloc(sizeof(int)*numprocs); NPROCS_WD2 = (int*)malloc(sizeof(int)*Num_Comm_World2); Comm_World_StartID2 = (int*)malloc(sizeof(int)*Num_Comm_World2); MPI_CommWD2 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World2); Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World2, &myworld2, MPI_CommWD2, NPROCS_ID2, Comm_World2, NPROCS_WD2, Comm_World_StartID2); MPI_Comm_size(MPI_CommWD2[myworld2],&numprocs2); MPI_Comm_rank(MPI_CommWD2[myworld2],&myid2); } else { parallel2_flag = 0; myworld2 = 0; } /**************************************************** SCF calculations for the two terminal structures ****************************************************/ /* check whether the restart is performed or not */ sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name); if (Scf_RestartFromFile==1 && (fp = fopen(file_neb_utot,"rb"))!=NULL){ if (myid==Host_ID){ fread(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp); fread(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp); } MPI_Bcast(&neb_atom_coordinates[0][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1); MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1); fclose(fp); } /* if the restart is not performed, perform the SCF calulations for the two terminals */ else { /* generate input files */ generate_input_files(fname_original,1); /* In case of parallel2_flag==1 */ if (parallel2_flag==1){ if (myworld2==0) index_images = 0; else index_images = NEB_Num_Images + 1; sprintf(fname1,"%s_%i",fname_original,index_images); argv[1] = fname1; neb_run( argv,MPI_CommWD2[myworld2],index_images,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); /* find a representative ID in the top world for ID=0 in the world2 */ i = 0; po = 0; do { if (Comm_World2[i]==0){ ID0 = i; po = 1; } i++; } while (po==0); /* find a representative ID in the top world for ID=1 in the world2 */ i = 0; po = 0; do { if (Comm_World2[i]==1){ ID1 = i; po = 1; } i++; } while (po==0); /* MPI broadcast of neb_atom_coordinates[0] and neb_atom_coordinates[NEB_Num_Images+1] */ MPI_Barrier(MPI_COMM_WORLD); for (i=0; i<=atomnum; i++){ MPI_Bcast(&neb_atom_coordinates[0][i][0], 20, MPI_DOUBLE, ID0, MPI_COMM_WORLD); MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][i][0], 20, MPI_DOUBLE, ID1, MPI_COMM_WORLD); } } /* if (parallel2_flag==1) */ /* In case of parallel2_flag==0 */ else { /* SCF calculation of a terminal with the index of 0 */ index_images = 0; sprintf(fname1,"%s_%i",fname_original,index_images); argv[1] = fname1; neb_run( argv,MPI_COMM_WORLD1,index_images,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); /* SCF calculation of a terminal with the index of (NEB_Num_Images + 1) */ index_images = NEB_Num_Images + 1; sprintf(fname1,"%s_%i",fname_original,index_images); argv[1] = fname1; neb_run( argv,MPI_COMM_WORLD1,index_images,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); } /* save *_0_rst/*.neb.utot in binary mode */ if (myid==Host_ID){ sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name); if ((fp = fopen(file_neb_utot,"wb")) != NULL){ fwrite(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp); fwrite(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp); fclose(fp); } else{ printf("Could not open a file %s in neb\n",file_neb_utot); } } } /* else */ /**************************************************** optimiziation for finding a minimum energy path connecting the two terminal structures. ****************************************************/ iter = 1; MD_Opt_OK = 0; dtime(&f0time); do { /* if iter==1, generate an input file for restarting */ if (iter==1) Generate_Restart_File(fname_original,neb_atom_coordinates); /* generate input files */ generate_input_files(fname_original,iter); /* In case of parallel1_mode==1 */ if (parallel1_flag==1){ sprintf(fname1,"%s_%i",fname_original,myworld1+1); argv[1] = fname1; dtime(&flatstime); neb_run( argv,MPI_CommWD1[myworld1],myworld1+1,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); dtime(&flatetime); sumtime += (flatetime-flatstime); /* MPI: All_Grid_Origin */ for (i=0; i<=(NEB_Num_Images+1); i++){ for (j=0; j<=3; j++) Tmp_Grid_Origin[i][j] = 0.0; } if (myid1==Host_ID){ Tmp_Grid_Origin[myworld1+1][1] = Grid_Origin[1]; Tmp_Grid_Origin[myworld1+1][2] = Grid_Origin[2]; Tmp_Grid_Origin[myworld1+1][3] = Grid_Origin[3]; } MPI_Barrier(MPI_COMM_WORLD); for (p=1; p<=NEB_Num_Images; p++){ MPI_Allreduce( &Tmp_Grid_Origin[p][0], &All_Grid_Origin[p][0], 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); } /* MPI: neb_atom_coordinates */ for (p=1; p<=NEB_Num_Images; p++){ for (i=0; i<=atomnum; i++){ for (j=0; j<20; j++){ tmp_neb_atom_coordinates[p][i][j] = 0.0; } } } if (myid1==Host_ID){ for (i=0; i<=atomnum; i++){ for (j=0; j<20; j++){ tmp_neb_atom_coordinates[myworld1+1][i][j] = neb_atom_coordinates[myworld1+1][i][j];; } } } MPI_Barrier(MPI_COMM_WORLD); for (p=1; p<=NEB_Num_Images; p++){ for (i=0; i<=atomnum; i++){ MPI_Allreduce( &tmp_neb_atom_coordinates[p][i][0], &neb_atom_coordinates[p][i][0], 20, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); } } } /* if (parallel1_flag==1) */ /* In case of parallel1_flag==0 */ else { for (p=1; p<=NEB_Num_Images; p++){ sprintf(fname1,"%s_%i",fname_original,p); argv[1] = fname1; neb_run( argv,MPI_COMM_WORLD1,p,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); /* store Grid_Origin */ All_Grid_Origin[p][1] = Grid_Origin[1]; All_Grid_Origin[p][2] = Grid_Origin[2]; All_Grid_Origin[p][3] = Grid_Origin[3]; } } MPI_Barrier(MPI_COMM_WORLD); /* calculate the gradients defined by the NEB method */ Calc_NEB_Gradients(neb_atom_coordinates); /* make a xyz file storing a set of structures of images at the current step */ Make_XYZ_File(fileXYZ,neb_atom_coordinates); /* update the atomic structures of the images */ Update_Coordinates(iter,neb_atom_coordinates); /* generate an input file for restarting */ Generate_Restart_File(fname_original,neb_atom_coordinates); /* increment of iter */ iter++; } while (MD_Opt_OK==0 && iter<=MD_IterNumber); MPI_Barrier(MPI_COMM_WORLD); /* freeing of arrays for the World1 */ if ( Num_Comm_World1<=numprocs ){ MPI_Comm_free(&MPI_CommWD1[myworld1]); free(MPI_CommWD1); free(Comm_World_StartID1); free(NPROCS_WD1); free(Comm_World1); free(NPROCS_ID1); } /* freeing of arrays for the World2 */ if ( Num_Comm_World2<=numprocs ){ MPI_Comm_free(&MPI_CommWD2[myworld2]); free(MPI_CommWD2); free(Comm_World_StartID2); free(NPROCS_WD2); free(Comm_World2); free(NPROCS_ID2); } /* freeing of arrays */ free_arrays(); /* show a final message */ dtime(&TEtime); printf("\nThe calculation was normally finished. (proc=%3d) TIME=%lf (s) flat time=%lf \n", myid,(TEtime-TStime),sumtime/(iter-1)); MPI_Finalize(); exit(0); } /* if ( NEB_Num_Images<=numprocs || numprocs==1 ) */ /**************************************************** One level parallelization is performed. Outer loop: images ****************************************************/ else{ dtime(&a0time); int syou,amari,g,bb; syou = NEB_Num_Images/numprocs; amari = NEB_Num_Images%numprocs; if(amari==0){ g=syou; } else{ g=syou+1; } Num_Comm_World1=numprocs; parallel1_flag = 1; NPROCS_ID1 = (int*)malloc(sizeof(int)*numprocs); Comm_World1 = (int*)malloc(sizeof(int)*numprocs); NPROCS_WD1 = (int*)malloc(sizeof(int)*Num_Comm_World1); Comm_World_StartID1 = (int*)malloc(sizeof(int)*Num_Comm_World1); MPI_CommWD1 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World1); Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World1, &myworld1, MPI_CommWD1, NPROCS_ID1, Comm_World1, NPROCS_WD1, Comm_World_StartID1); MPI_Comm_size(MPI_CommWD1[myworld1],&numprocs1); MPI_Comm_rank(MPI_CommWD1[myworld1],&myid1); /**************************************************** allocate processes to each image in the World2 ****************************************************/ Num_Comm_World2 = 2; if ( Num_Comm_World2<=numprocs ){ parallel2_flag = 1; NPROCS_ID2 = (int*)malloc(sizeof(int)*numprocs); Comm_World2 = (int*)malloc(sizeof(int)*numprocs); NPROCS_WD2 = (int*)malloc(sizeof(int)*Num_Comm_World2); Comm_World_StartID2 = (int*)malloc(sizeof(int)*Num_Comm_World2); MPI_CommWD2 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World2); Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World2, &myworld2, MPI_CommWD2, NPROCS_ID2, Comm_World2, NPROCS_WD2, Comm_World_StartID2); MPI_Comm_size(MPI_CommWD2[myworld2],&numprocs2); MPI_Comm_rank(MPI_CommWD2[myworld2],&myid2); } /**************************************************** SCF calculations for the two terminal structures ****************************************************/ /* check whether the restart is performed or not */ sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name); if (Scf_RestartFromFile==1 && (fp = fopen(file_neb_utot,"rb"))!=NULL){ if (myid==Host_ID){ fread(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp); fread(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp); } MPI_Bcast(&neb_atom_coordinates[0][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1); MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1); fclose(fp); } /* if the restart is not performed, perform the SCF calulations for the two terminals */ else { /* generate input files */ generate_input_files(fname_original,1); /* In case of parallel2_flag==1 */ if (parallel2_flag==1){ if (myworld2==0) index_images = 0; else index_images = NEB_Num_Images + 1; sprintf(fname1,"%s_%i",fname_original,index_images); argv[1] = fname1; neb_run( argv,MPI_CommWD2[myworld2],index_images,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); /* find a representative ID in the top world for ID=0 in the world2 */ i = 0; po = 0; do { if (Comm_World2[i]==0){ ID0 = i; po = 1; } i++; } while (po==0); /* find a representative ID in the top world for ID=1 in the world2 */ i = 0; po = 0; do { if (Comm_World2[i]==1){ ID1 = i; po = 1; } i++; } while (po==0); /* MPI broadcast of neb_atom_coordinates[0] and neb_atom_coordinates[NEB_Num_Images+1] */ MPI_Barrier(MPI_COMM_WORLD); for (i=0; i<=atomnum; i++){ MPI_Bcast(&neb_atom_coordinates[0][i][0], 20, MPI_DOUBLE, ID0, MPI_COMM_WORLD); MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][i][0], 20, MPI_DOUBLE, ID1, MPI_COMM_WORLD); } } /* if (parallel2_flag==1) */ /* save *_0_rst/*.neb.utot in binary mode */ if (myid==Host_ID){ sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name); if ((fp = fopen(file_neb_utot,"wb")) != NULL){ fwrite(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp); fwrite(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp); fclose(fp); } else{ printf("Could not open a file %s in neb\n",file_neb_utot); } } } /* else */ /* amari world */ if(amari!=0){ Num_Comm_World3=amari; NPROCS_ID3 = (int*)malloc(sizeof(int)*numprocs); Comm_World3 = (int*)malloc(sizeof(int)*numprocs); NPROCS_WD3 = (int*)malloc(sizeof(int)*Num_Comm_World3); Comm_WOrld_StartID3 = (int*)malloc(sizeof(int)*Num_Comm_World3); MPI_CommWD3 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World3); Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World3, &myworld3, MPI_CommWD3, NPROCS_ID3, Comm_World3, NPROCS_WD3, Comm_WOrld_StartID3); MPI_Comm_size(MPI_CommWD3[myworld3],&numprocs3); MPI_Comm_rank(MPI_CommWD3[myworld3],&myid3); } /* amari world */ dtime(&a1time); /**************************************************** optimization for finding a minimum energy path connecting the two terminal structures. ****************************************************/ iter = 1; MD_Opt_OK = 0; dtime(&b0time); do { /* if iter==1, generate an input file for restarting */ if (iter==1) Generate_Restart_File(fname_original,neb_atom_coordinates); /* generate input files */ generate_input_files(fname_original,iter); for(h=1;h<=g;h++){ if(h==g && g==syou+1){ sprintf(fname1,"%s_%i",fname_original,myworld3+1+(h-1)*numprocs); } else{ sprintf(fname1,"%s_%i",fname_original,myworld1+1+(h-1)*numprocs); } argv[1] = fname1; if((h==g) && (g==syou+1)){ neb_run( argv,MPI_CommWD3[myworld3],myworld3+1+(h-1)*numprocs,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); if (myid3==Host_ID){ Tmp_Grid_Origin[myworld3+1+(h-1)*numprocs][1] = Grid_Origin[1]; Tmp_Grid_Origin[myworld3+1+(h-1)*numprocs][2] = Grid_Origin[2]; Tmp_Grid_Origin[myworld3+1+(h-1)*numprocs][3] = Grid_Origin[3]; } } else{ neb_run( argv,MPI_CommWD1[myworld1],myworld1+1+(h-1)*numprocs,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); if(h==1){ for (i=0; i<=(NEB_Num_Images+1); i++){ for (j=0; j<=3; j++) Tmp_Grid_Origin[i][j] = 0.0; } } if (myid1==Host_ID){ Tmp_Grid_Origin[myworld1+1+(h-1)*numprocs][1] = Grid_Origin[1]; Tmp_Grid_Origin[myworld1+1+(h-1)*numprocs][2] = Grid_Origin[2]; Tmp_Grid_Origin[myworld1+1+(h-1)*numprocs][3] = Grid_Origin[3]; } } /* MPI: All_Grid_Origin */ MPI_Barrier(MPI_COMM_WORLD); } MPI_Barrier(MPI_COMM_WORLD); for (p=1; p<=NEB_Num_Images; p++){ MPI_Allreduce( &Tmp_Grid_Origin[p][0], &All_Grid_Origin[p][0], 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); } /* MPI: neb_atom_coordinates */ for (p=1; p<=NEB_Num_Images; p++){ for (i=0; i<=atomnum; i++){ for (j=0; j<20; j++){ tmp_neb_atom_coordinates[p][i][j] = 0.0; } } } for(h=1;h<=g;h++){ if(h==g && g==syou+1){ if (myid3==Host_ID){ for (i=0; i<=atomnum; i++){ for (j=0; j<20; j++){ tmp_neb_atom_coordinates[myworld3+1+(h-1)*numprocs][i][j] = neb_atom_coordinates[myworld3+1+(h-1)*numprocs][i][j];; } } } } else{ if (myid1==Host_ID){ for (i=0; i<=atomnum; i++){ for (j=0; j<20; j++){ tmp_neb_atom_coordinates[myworld1+1+(h-1)*numprocs][i][j] = neb_atom_coordinates[myworld1+1+(h-1)*numprocs][i][j];; } } } } } MPI_Barrier(MPI_COMM_WORLD); for (p=1; p<=NEB_Num_Images; p++){ for (i=0; i<=atomnum; i++){ MPI_Allreduce( &tmp_neb_atom_coordinates[p][i][0], &neb_atom_coordinates[p][i][0], 20, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); } } MPI_Barrier(MPI_COMM_WORLD); /* calculate the gradients defined by the NEB method */ Calc_NEB_Gradients(neb_atom_coordinates); /* make a xyz file storing a set of structures of images at the current step */ Make_XYZ_File(fileXYZ,neb_atom_coordinates); /* update the atomic structures of the images */ Update_Coordinates(iter,neb_atom_coordinates); /* generate an input file for restarting */ Generate_Restart_File(fname_original,neb_atom_coordinates); /* increment of iter */ iter++; } while (MD_Opt_OK==0 && iter<=MD_IterNumber); dtime(&b1time); MPI_Barrier(MPI_COMM_WORLD); /* freeing of arrays for the World1 */ dtime(&c0time); MPI_Comm_free(&MPI_CommWD1[myworld1]); free(MPI_CommWD1); free(Comm_World_StartID1); free(NPROCS_WD1); free(Comm_World1); free(NPROCS_ID1); /* freeing of arrays for the World2 */ MPI_Comm_free(&MPI_CommWD2[myworld2]); free(MPI_CommWD2); free(Comm_World_StartID2); free(NPROCS_WD2); free(Comm_World2); free(NPROCS_ID2); /* freeing of arrays for the World3 */ if(amari!=0){ MPI_Comm_free(&MPI_CommWD3[myworld3]); free(MPI_CommWD3); free(Comm_WOrld_StartID3); free(NPROCS_WD3); free(Comm_World3); free(NPROCS_ID3); } /* freeing of arrays */ free_arrays(); /* show a final message */ dtime(&c1time); dtime(&TEtime); printf("\nThe calculation was normally finished. (proc=%3d) TIME=%lf (s)\n",myid,(TEtime-TStime)); MPI_Finalize(); exit(0); } } /* if (PN==0) */ /**************************************************************************** If MD.NEB.Parallel.Number (PN) is specifed, corresponding to 1<=MD.NEB.Parallel.Number (PN) *****************************************************************************/ else{ dtime(&a0time); int syou,amari,g,bb; syou = NEB_Num_Images/PN; amari = NEB_Num_Images%PN; if (amari==0) g = syou; else g = syou + 1; Num_Comm_World1 = PN; if(PN>numprocs){ if (myid==Host_ID) printf("PN is larger than the number of MPI processes.\n"); MPI_Finalize(); exit(0); } parallel1_flag = 1; NPROCS_ID1 = (int*)malloc(sizeof(int)*numprocs); Comm_World1 = (int*)malloc(sizeof(int)*numprocs); NPROCS_WD1 = (int*)malloc(sizeof(int)*Num_Comm_World1); Comm_World_StartID1 = (int*)malloc(sizeof(int)*Num_Comm_World1); MPI_CommWD1 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World1); Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World1, &myworld1, MPI_CommWD1, NPROCS_ID1, Comm_World1, NPROCS_WD1, Comm_World_StartID1); MPI_Comm_size(MPI_CommWD1[myworld1],&numprocs1); MPI_Comm_rank(MPI_CommWD1[myworld1],&myid1); /**************************************************** allocate processes to each image in the World2 ****************************************************/ Num_Comm_World2 = 2; if ( Num_Comm_World2<=numprocs ){ parallel2_flag = 1; NPROCS_ID2 = (int*)malloc(sizeof(int)*numprocs); Comm_World2 = (int*)malloc(sizeof(int)*numprocs); NPROCS_WD2 = (int*)malloc(sizeof(int)*Num_Comm_World2); Comm_World_StartID2 = (int*)malloc(sizeof(int)*Num_Comm_World2); MPI_CommWD2 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World2); Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World2, &myworld2, MPI_CommWD2, NPROCS_ID2, Comm_World2, NPROCS_WD2, Comm_World_StartID2); MPI_Comm_size(MPI_CommWD2[myworld2],&numprocs2); MPI_Comm_rank(MPI_CommWD2[myworld2],&myid2); } else{ parallel2_flag = 0; myworld2 = 0; } /**************************************************** SCF calculations for the two terminal structures ****************************************************/ /* check whether the restart is performed or not */ sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name); if (Scf_RestartFromFile==1 && (fp = fopen(file_neb_utot,"rb"))!=NULL){ if (myid==Host_ID){ fread(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp); fread(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp); } MPI_Bcast(&neb_atom_coordinates[0][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1); MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][0][0], 1, MPI_DOUBLE, Host_ID, MPI_COMM_WORLD1); fclose(fp); } /* if the restart is not performed, perform the SCF calulations for the two terminals */ else { /* generate input files */ generate_input_files(fname_original,1); /*********************************** if (PN==1) ***********************************/ if (PN==1){ for (i=0; i<=1; i++){ if (i==0) index_images = 0; else index_images = NEB_Num_Images + 1; sprintf(fname1,"%s_%i",fname_original,index_images); argv[1] = fname1; neb_run( argv,MPI_CommWD2[myworld2],index_images,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); } } /*********************************** if (2<=PN) ***********************************/ else { if (myworld2==0) index_images = 0; else index_images = NEB_Num_Images + 1; sprintf(fname1,"%s_%i",fname_original,index_images); argv[1] = fname1; neb_run( argv,MPI_CommWD2[myworld2],index_images,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); /* find a representative ID in the top world for ID=0 in the world2 */ i = 0; po = 0; do { if (Comm_World2[i]==0){ ID0 = i; po = 1; } i++; } while (po==0); /* find a representative ID in the top world for ID=1 in the world2 */ i = 0; po = 0; do { if (Comm_World2[i]==1){ ID1 = i; po = 1; } i++; } while (po==0); /* MPI broadcast of neb_atom_coordinates[0] and neb_atom_coordinates[NEB_Num_Images+1] */ MPI_Barrier(MPI_COMM_WORLD); for (i=0; i<=atomnum; i++){ MPI_Bcast(&neb_atom_coordinates[0][i][0], 20, MPI_DOUBLE, ID0, MPI_COMM_WORLD); MPI_Bcast(&neb_atom_coordinates[NEB_Num_Images+1][i][0], 20, MPI_DOUBLE, ID1, MPI_COMM_WORLD); } } /* else which corresponds to if (2<=PN) */ /* save *_0_rst/*.neb.utot in binary mode */ if (myid==Host_ID){ sprintf(file_neb_utot,"%s%s_0_rst/%s.neb.utot",filepath,system_name,system_name); if ((fp = fopen(file_neb_utot,"wb")) != NULL){ fwrite(&neb_atom_coordinates[0][0][0],sizeof(double),1,fp); fwrite(&neb_atom_coordinates[NEB_Num_Images+1][0][0],sizeof(double),1,fp); fclose(fp); } else{ printf("Could not open a file %s in neb\n",file_neb_utot); } } } /* else */ /**************************************************** make amari world ****************************************************/ if (amari!=0){ Num_Comm_World3 = amari; NPROCS_ID3 = (int*)malloc(sizeof(int)*numprocs); Comm_World3 = (int*)malloc(sizeof(int)*numprocs); NPROCS_WD3 = (int*)malloc(sizeof(int)*Num_Comm_World3); Comm_WOrld_StartID3 = (int*)malloc(sizeof(int)*Num_Comm_World3); MPI_CommWD3 = (MPI_Comm*)malloc(sizeof(MPI_Comm)*Num_Comm_World3); Make_Comm_Worlds(MPI_COMM_WORLD, myid, numprocs, Num_Comm_World3, &myworld3, MPI_CommWD3, NPROCS_ID3, Comm_World3, NPROCS_WD3, Comm_WOrld_StartID3); MPI_Comm_size(MPI_CommWD3[myworld3],&numprocs3); MPI_Comm_rank(MPI_CommWD3[myworld3],&myid3); } /* amari world */ dtime(&a1time); /**************************************************** optimization for finding a minimum energy path connecting the two terminal structures. ****************************************************/ iter = 1; MD_Opt_OK = 0; dtime(&b0time); do { /* if iter==1, generate an input file for restarting */ if (iter==1) Generate_Restart_File(fname_original,neb_atom_coordinates); /* generate input files */ generate_input_files(fname_original,iter); for (h=1; h<=g; h++){ if( h==g && g==(syou+1) ){ sprintf(fname1,"%s_%i",fname_original,myworld3+1+(h-1)*PN); } else{ sprintf(fname1,"%s_%i",fname_original,myworld1+1+(h-1)*PN); } argv[1] = fname1; if( h==g && g==(syou+1) ){ neb_run( argv,MPI_CommWD3[myworld3],myworld3+1+(h-1)*PN,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); if (myid3==Host_ID){ Tmp_Grid_Origin[myworld3+1+(h-1)*PN][1] = Grid_Origin[1]; Tmp_Grid_Origin[myworld3+1+(h-1)*PN][2] = Grid_Origin[2]; Tmp_Grid_Origin[myworld3+1+(h-1)*PN][3] = Grid_Origin[3]; } } else{ neb_run( argv,MPI_CommWD1[myworld1],myworld1+1+(h-1)*PN,neb_atom_coordinates, WhatSpecies_NEB,Spe_WhatAtom_NEB,SpeName_NEB ); if(h==1){ for (i=0; i<=(NEB_Num_Images+1); i++){ for (j=0; j<=3; j++) Tmp_Grid_Origin[i][j] = 0.0; } } if (myid1==Host_ID){ Tmp_Grid_Origin[myworld1+1+(h-1)*PN][1] = Grid_Origin[1]; Tmp_Grid_Origin[myworld1+1+(h-1)*PN][2] = Grid_Origin[2]; Tmp_Grid_Origin[myworld1+1+(h-1)*PN][3] = Grid_Origin[3]; } } /* MPI: All_Grid_Origin */ MPI_Barrier(MPI_COMM_WORLD); } /* for (h=1; h<=g; h++) */ MPI_Barrier(MPI_COMM_WORLD); for (p=1; p<=NEB_Num_Images; p++){ MPI_Allreduce( &Tmp_Grid_Origin[p][0], &All_Grid_Origin[p][0], 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); } /* MPI: neb_atom_coordinates */ for (p=1; p<=NEB_Num_Images; p++){ for (i=0; i<=atomnum; i++){ for (j=0; j<20; j++){ tmp_neb_atom_coordinates[p][i][j] = 0.0; } } } for(h=1;h<=g;h++){ if(h==g && g==syou+1){ if (myid3==Host_ID){ for (i=0; i<=atomnum; i++){ for (j=0; j<20; j++){ tmp_neb_atom_coordinates[myworld3+1+(h-1)*PN][i][j] = neb_atom_coordinates[myworld3+1+(h-1)*PN][i][j];; } } } } else{ if (myid1==Host_ID){ for (i=0; i<=atomnum; i++){ for (j=0; j<20; j++){ tmp_neb_atom_coordinates[myworld1+1+(h-1)*PN][i][j] = neb_atom_coordinates[myworld1+1+(h-1)*PN][i][j];; } } } } } MPI_Barrier(MPI_COMM_WORLD); for (p=1; p<=NEB_Num_Images; p++){ for (i=0; i<=atomnum; i++){ MPI_Allreduce( &tmp_neb_atom_coordinates[p][i][0], &neb_atom_coordinates[p][i][0], 20, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD ); } } MPI_Barrier(MPI_COMM_WORLD); /* calculate the gradients defined by the NEB method */ Calc_NEB_Gradients(neb_atom_coordinates); /* make a xyz file storing a set of structures of images at the current step */ Make_XYZ_File(fileXYZ,neb_atom_coordinates); /* update the atomic structures of the images */ Update_Coordinates(iter,neb_atom_coordinates); /* generate an input file for restarting */ Generate_Restart_File(fname_original,neb_atom_coordinates); /* increment of iter */ iter++; } while (MD_Opt_OK==0 && iter<=MD_IterNumber); dtime(&b1time); MPI_Barrier(MPI_COMM_WORLD); /* freeing of arrays for the World1 */ dtime(&c0time); MPI_Comm_free(&MPI_CommWD1[myworld1]); free(MPI_CommWD1); free(Comm_World_StartID1); free(NPROCS_WD1); free(Comm_World1); free(NPROCS_ID1); /* freeing of arrays for the World2 */ MPI_Comm_free(&MPI_CommWD2[myworld2]); free(MPI_CommWD2); free(Comm_World_StartID2); free(NPROCS_WD2); free(Comm_World2); free(NPROCS_ID2); /* freeing of arrays for the World3 */ if(amari!=0){ MPI_Comm_free(&MPI_CommWD3[myworld3]); free(MPI_CommWD3); free(Comm_WOrld_StartID3); free(NPROCS_WD3); free(Comm_World3); free(NPROCS_ID3); } /* freeing of arrays */ free_arrays(); /* show a final message */ dtime(&c1time); dtime(&TEtime); printf("\nThe calculation was normally finished. (proc=%3d) TIME=%lf (s)\n",myid,(TEtime-TStime)); MPI_Finalize(); exit(0); } } void Generate_Restart_File(char fname_original[YOUSO10], double ***neb_atom_coordinates) { int i,p,Gc_AN,c,n1,k,po; int restart_flag; int unit_flag; double c1,c2,c3; double tmpxyz[4]; char st[800]; char st1[800]; char rm_operate[YOUSO10]; char fname1[YOUSO10]; char fname2[YOUSO10]; char keyword[YOUSO10]; FILE *fp1,*fp2; char *tp; int numprocs,myid; /* MPI */ MPI_Comm_size(MPI_COMM_WORLD,&numprocs); MPI_Comm_rank(MPI_COMM_WORLD,&myid); /* generate the input files */ if (myid==Host_ID){ /* initialize */ restart_flag = 0; unit_flag = 0; /* the new input file */ sprintf(fname1,"%s#",fname_original); fp1 = fopen(fname1,"w"); fseek(fp1,0,SEEK_END); /* the original input file */ fp2 = fopen(fname_original,"r"); if (fp2!=NULL){ while (fgets(st,800,fp2)!=NULL){ string_tolower(st,st1); /* find the specification of and forget it */ fgets(st,800,fp2); } } /* p */ /* other cases */ if (po==0) fprintf(fp1,"%s",st); } /* else */ } /* while */ /* close fp2 */ fclose(fp2); } /* add the restart flag if it was not found. */ if (restart_flag==0){ fprintf(fp1,"\n\nscf.restart on\n"); } /* add Atoms.SpeciesAndCoordinates.Unit if it was not found. */ if (unit_flag==0){ fprintf(fp1,"Atoms.SpeciesAndCoordinates.Unit AU\n"); } /* add atomic coordinates of the images */ for (p=1; p<=NEB_Num_Images; p++){ fp2 = fopen(fname_original,"r"); if (fp2!=NULL){ while (fgets(st,800,fp2)!=NULL){ string_tolower(st,st1); if (strncmp(st1,"",28)==0){ fprintf(fp1,"NEB%d.Atoms.SpeciesAndCoordinates>\n",p); } } /* close fp2 */ fclose(fp2); } /* if (fp2!=NULL){ */ } /* p */ /* fclose */ fclose(fp1); } /* if (myid==Host_ID) */ } void Make_XYZ_File(char fname[YOUSO10], double ***neb_atom_coordinates) { FILE *fp; int p,k,i,j; int myid; char file_neb_utot[YOUSO10]; MPI_Comm_rank(MPI_COMM_WORLD,&myid); if (myid==Host_ID){ /* save *.neb.xyz */ if ((fp = fopen(fname,"w")) != NULL){ for (p=0; p<=(NEB_Num_Images+1); p++){ fprintf(fp,"%i \n",atomnum); fprintf(fp," Image index = %3d Energy= %8.5f (Hatree)\n",p,neb_atom_coordinates[p][0][0]); for (k=1; k<=atomnum; k++){ i = WhatSpecies_NEB[k]; j = Spe_WhatAtom_NEB[i]; fprintf(fp,"%4s %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f\n", Atom_Symbol[j], neb_atom_coordinates[p][k][1]*BohrR, neb_atom_coordinates[p][k][2]*BohrR, neb_atom_coordinates[p][k][3]*BohrR, -neb_atom_coordinates[p][k][4], -neb_atom_coordinates[p][k][5], -neb_atom_coordinates[p][k][6]); } } fclose(fp); } else{ printf("failure of saving the xyz file.\n"); fclose(fp); } } /* if (myid==Host_ID) */ } void Update_Coordinates(int iter, double ***neb_atom_coordinates) { int p,j,k; double norm,dis,dif; double Max_Force,Max_Step; double tmp0,sum,sum_E,sum_dis; int numprocs,myid; char fileOPT[YOUSO10]; FILE *fp_OPT; char fileE[YOUSO10]; FILE *fp_E; MPI_Comm_size(MPI_COMM_WORLD,&numprocs); MPI_Comm_rank(MPI_COMM_WORLD,&myid); /********************************************* calculate the norm of gradient *********************************************/ norm = 0.0; for (p=1; p<=NEB_Num_Images; p++){ for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ norm += neb_atom_coordinates[p][j][k+3]*neb_atom_coordinates[p][j][k+3]; } } } norm = sqrt(norm); /**************************************************** find the maximum value of force ****************************************************/ Max_Force = 0.0; for (p=1; p<=NEB_Num_Images; p++){ for (j=1; j<=atomnum; j++){ sum = 0.0; for (k=1; k<=3; k++){ tmp0 = neb_atom_coordinates[p][j][k+3]; sum += tmp0*tmp0; } sum = sqrt(sum); if (Max_Force */ tmp1 = 0.0; for (p2=1; p2<=NEB_Num_Images; p2++){ for (i2=1; i2<=atomnum; i2++){ for (j2=1; j2<=3; j2++){ dif_x = His_neb_atom_coordinates[0][p2][i2][j2 ]-His_neb_atom_coordinates[1][p2][i2][j2 ]; dif_g = His_neb_atom_coordinates[0][p2][i2][j2+3]-His_neb_atom_coordinates[1][p2][i2][j2+3]; tmp1 += dif_x*dif_g; } } } /* < Delta g^(n) | Vec0 > */ tmp2 = 0.0; for (p2=1; p2<=NEB_Num_Images; p2++){ for (i2=1; i2<=atomnum; i2++){ for (j2=1; j2<=3; j2++){ dif_g = His_neb_atom_coordinates[0][p2][i2][j2+3]-His_neb_atom_coordinates[1][p2][i2][j2+3]; tmp2 += dif_g*Vec0[p2][i2][j2]; } } } /* update InvHess */ c0 = (tmp1 + tmp2)/(tmp1*tmp1); c1 = -1.0/tmp1; for (p1=1; p1<=NEB_Num_Images; p1++){ for (i1=1; i1<=atomnum; i1++){ for (j1=1; j1<=3; j1++){ for (p2=1; p2<=NEB_Num_Images; p2++){ for (i2=1; i2<=atomnum; i2++){ for (j2=1; j2<=3; j2++){ dif_x1 = His_neb_atom_coordinates[0][p1][i1][j1] -His_neb_atom_coordinates[1][p1][i1][j1]; dif_x2 = His_neb_atom_coordinates[0][p2][i2][j2] -His_neb_atom_coordinates[1][p2][i2][j2]; InvHess[p1][i1][j1][p2][i2][j2] += c0*dif_x1*dif_x2; dif_x1 = Vec0[p1][i1][j1]; dif_x2 = His_neb_atom_coordinates[0][p2][i2][j2] -His_neb_atom_coordinates[1][p2][i2][j2]; InvHess[p1][i1][j1][p2][i2][j2] += c1*dif_x1*dif_x2; dif_x1 = His_neb_atom_coordinates[0][p1][i1][j1] -His_neb_atom_coordinates[1][p1][i1][j1]; dif_x2 = Vec0[p2][i2][j2]; InvHess[p1][i1][j1][p2][i2][j2] += c1*dif_x1*dif_x2; } } } } } } } /* else */ /* update atomic coordinates */ for (p1=1; p1<=NEB_Num_Images; p1++){ for (i1=1; i1<=atomnum; i1++){ for (j1=1; j1<=3; j1++){ sum1 = 0.0; for (p2=1; p2<=NEB_Num_Images; p2++){ for (i2=1; i2<=atomnum; i2++){ for (j2=1; j2<=3; j2++){ tmp = neb_atom_coordinates[p2][i2][j2+3]; sum1 += InvHess[p1][i1][j1][p2][i2][j2]*tmp; } } } neb_atom_coordinates[p1][i1][j1] -= sum1; } } } /* In case of a too large updating, do a modest updating */ Max_Step = 0.0; for (p=1; p<=NEB_Num_Images; p++){ for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ dif_x = fabs(neb_atom_coordinates[p][j][k] - His_neb_atom_coordinates[0][p][j][k]); if (Max_Step=vi && vi<=vm) || (vpvm) ){ Dvmax = MAX(fabs(vp-vi),fabs(vm-vi)); Dvmin = MIN(fabs(vp-vi),fabs(vm-vi)); if (vm<=vp){ for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ tangents[p][j][k] = Dvmax*(neb_atom_coordinates[p+1][j][k] - neb_atom_coordinates[p ][j][k]) + Dvmin*(neb_atom_coordinates[p ][j][k] - neb_atom_coordinates[p-1][j][k]); } } } else { for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ tangents[p][j][k] = Dvmin*(neb_atom_coordinates[p+1][j][k] - neb_atom_coordinates[p ][j][k]) + Dvmax*(neb_atom_coordinates[p ][j][k] - neb_atom_coordinates[p-1][j][k]); } } } } else { printf("unknown situation\n"); } /* normalization of tangent */ sum = 0.0; for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ sum += tangents[p][j][k]*tangents[p][j][k]; } } sum = 1.0/sqrt(sum); for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ tangents[p][j][k] *= sum; } } } /* calculate forces defined by the NEB method */ for (p=1; p<=NEB_Num_Images; p++){ /**************************************************************** Perpendicular force ****************************************************************/ /* inner product of gradient[p] and tangents[p] */ prod1 = 0.0; for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ prod1 += neb_atom_coordinates[p][j][k+16]*tangents[p][j][k]; } } /* calculate the perpendicular gradient */ for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ neb_atom_coordinates[p][j][k+3] = neb_atom_coordinates[p][j][k+16] - prod1*tangents[p][j][k]; } } /**************************************************************** Parallel force ****************************************************************/ /* force constant */ sconst = NEB_Spring_Const; /* calculate the distance between Ri+1 and Ri */ dis1 = 0.0; for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ dif = neb_atom_coordinates[p+1][j][k] - neb_atom_coordinates[p][j][k]; dis1 += dif*dif; } } dis1 = sqrt(dis1); /* calculate the distance between Ri and Ri-1 */ dis2 = 0.0; for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ dif = neb_atom_coordinates[p][j][k] - neb_atom_coordinates[p-1][j][k]; dis2 += dif*dif; } } dis2 = sqrt(dis2); /* calculate the parallel gradient by springs, and add the contribution based on Eq. (12) in JCP 113, 9978 (2000) */ for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ neb_atom_coordinates[p][j][k+3] += -sconst*(dis1 - dis2)*tangents[p][j][k]; } } } /* p */ /**************************************************************** introduce the contraints defined by user by setting gradients zero 1: fixed 0: relaxed ****************************************************************/ for (p=1; p<=NEB_Num_Images; p++){ for (j=1; j<=atomnum; j++){ for (k=1; k<=3; k++){ if (atomFixedXYZ[j][k]==1){ neb_atom_coordinates[p][j][k+3] = 0.0; } } } } /* freeing of arrays */ for (i=0; i<(NEB_Num_Images+2); i++){ for (j=0; j<(atomnum+1); j++){ free(tangents[i][j]); } free(tangents[i]); } free(tangents); } void generate_input_files(char *file, int iter) { int i,p,Gc_AN,c,n1,k; int restart_flag,fixed_flag; int unit_flag,level_flag; double c1,c2,c3; double tmpxyz[4]; char st[800]; char st1[800]; char rm_operate[YOUSO10]; char fname[YOUSO10]; char fname1[YOUSO10]; char fname2[YOUSO10]; FILE *fp1,*fp2; char *tp; int numprocs,myid; /* MPI */ MPI_Comm_size(MPI_COMM_WORLD,&numprocs); MPI_Comm_rank(MPI_COMM_WORLD,&myid); /* get system.name */ if (input_open(file)==0){ MPI_Finalize(); exit(0); } input_string("System.CurrrentDirectory",filepath,"./"); input_string("System.Name",filename,"default"); input_close(); /* generate the input files */ if (myid==Host_ID){ for (p=0; p<=(NEB_Num_Images+1); p++){ /* initialize */ restart_flag = 0; fixed_flag = 0; unit_flag = 0; level_flag = 0; /* the new input file */ sprintf(fname1,"%s_%i",file,p); fp1 = fopen(fname1,"w"); fseek(fp1,0,SEEK_END); /* the original input file */ fp2 = fopen(file,"r"); if (fp2!=NULL){ while (fgets(st,800,fp2)!=NULL){ string_tolower(st,st1); /* find the specification of ") ) { /* format error */ if (myid==Host_ID){ printf("Format error for Atoms.Unitvectors\n"); } MPI_Finalize(); exit(0); } /* Ang to AU */ if (unitvector_unit==0){ for (i=1; i<=3; i++){ tv[i][1] = tv[i][1]/BohrR; tv[i][2] = tv[i][2]/BohrR; tv[i][3] = tv[i][3]/BohrR; } } } if (unitvectors_flag==0){ if (myid==Host_ID){ printf("A common unit cell for two terminal structures has to be specified.\n");fflush(stdout); } MPI_Finalize(); exit(0); } /************************************************************** read atomic coodinates given by Atoms.SpeciesAndCoordinates. **************************************************************/ s_vec[0]="Ang"; s_vec[1]="AU"; s_vec[2]="FRAC"; i_vec[0]= 0; i_vec[1]= 1; i_vec[2]= 2; input_string2int("Atoms.SpeciesAndCoordinates.Unit", &coordinates_unit,3,s_vec,i_vec); if (fp=input_find("")) { /* format error */ if (myid==Host_ID){ printf("Format error for Atoms.SpeciesAndCoordinates\n");fflush(stdout); } MPI_Finalize(); exit(0); } /* Ang to AU */ if (coordinates_unit==0){ for (i=1; i<=atomnum; i++){ neb_atom_coordinates[0][i][1] /= BohrR; neb_atom_coordinates[0][i][2] /= BohrR; neb_atom_coordinates[0][i][3] /= BohrR; } } } /* FRAC to AU */ if (coordinates_unit==2){ /* The fractional coordinates should be kept within 0 to 1. */ for (i=1; i<=atomnum; i++){ for (k=1; k<=3; k++){ itmp = (int)neb_atom_coordinates[0][i][k]; if (1.0")) { /* format error */ if (myid==Host_ID){ printf("Format error for NEB.Atoms.SpeciesAndCoordinates\n"); } MPI_Finalize(); exit(0); } /* Ang to AU */ if (coordinates_unit==0){ for (i=1; i<=atomnum; i++){ neb_atom_coordinates[NEB_Num_Images+1][i][1] /= BohrR; neb_atom_coordinates[NEB_Num_Images+1][i][2] /= BohrR; neb_atom_coordinates[NEB_Num_Images+1][i][3] /= BohrR; } } } /* FRAC to AU */ if (coordinates_unit==2){ /* The fractional coordinates should be kept within 0 to 1. */ for (i=1; i<=atomnum; i++){ for (k=1; k<=3; k++){ itmp = (int)neb_atom_coordinates[NEB_Num_Images+1][i][k]; if (1.0",p); if (!input_last(keyword)) { /* format error */ if (myid==Host_ID){ printf("Format error for NEB%d.Atoms.SpeciesAndCoordinates\n",p); } MPI_Finalize(); exit(0); } } } } else { /***************************************************************** Making the centroid of two terminal coordinates equivalent. This treatment allows us to hold the equivalence of the centroid of all the images. Thus, Grid_Origin also becomes equivalent to eath other. *****************************************************************/ xc0 = 0.0; yc0 = 0.0; zc0 = 0.0; xc1 = 0.0; yc1 = 0.0; zc1 = 0.0; for (i=1; i<=atomnum; i++){ xc0 += neb_atom_coordinates[0][i][1]; yc0 += neb_atom_coordinates[0][i][2]; zc0 += neb_atom_coordinates[0][i][3]; xc1 += neb_atom_coordinates[NEB_Num_Images+1][i][1]; yc1 += neb_atom_coordinates[NEB_Num_Images+1][i][2]; zc1 += neb_atom_coordinates[NEB_Num_Images+1][i][3]; } xc0 = xc0/(double)atomnum; yc0 = yc0/(double)atomnum; zc0 = zc0/(double)atomnum; xc1 = xc1/(double)atomnum; yc1 = yc1/(double)atomnum; zc1 = zc1/(double)atomnum; xc = 0.5*(xc0 + xc1); yc = 0.5*(yc0 + yc1); zc = 0.5*(zc0 + zc1); for (i=1; i<=atomnum; i++){ neb_atom_coordinates[0][i][1] += -xc0 + xc; neb_atom_coordinates[0][i][2] += -yc0 + yc; neb_atom_coordinates[0][i][3] += -zc0 + zc; neb_atom_coordinates[NEB_Num_Images+1][i][1] += -xc1 + xc; neb_atom_coordinates[NEB_Num_Images+1][i][2] += -yc1 + yc; neb_atom_coordinates[NEB_Num_Images+1][i][3] += -zc1 + zc; } /***************************************************************** generate atomic coordinates for images *****************************************************************/ for (p=1; p<=NEB_Num_Images; p++){ for (i=1; i<=atomnum; i++){ for (j=1; j<=3; j++){ dif = neb_atom_coordinates[NEB_Num_Images+1][i][j]-neb_atom_coordinates[0][i][j]; neb_atom_coordinates[p][i][j] = neb_atom_coordinates[0][i][j] + dif*(double)p/(double)(NEB_Num_Images+1); } } } } /***************************************************************** read atomFixedXYZ set fixed atomic position in geometry optimization and MD: 1: fixed 0: relaxed *****************************************************************/ if (fp=input_find("") ) { /* format error */ if (myid==Host_ID){ printf("Format error for MD.Fixed.XYZ\n"); } MPI_Finalize(); exit(0); } } /* close the input file */ input_close(); } void allocate_arrays() { int i,j,k,m; int p1,i1,j1,p2,i2,j2; neb_atom_coordinates = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2)); for (i=0; i<(NEB_Num_Images+2); i++){ neb_atom_coordinates[i] = (double**)malloc(sizeof(double*)*(atomnum+1)); for (j=0; j<(atomnum+1); j++){ neb_atom_coordinates[i][j] = (double*)malloc(sizeof(double)*20); for (k=0; k<20; k++) neb_atom_coordinates[i][j][k] = 0.0; } } tmp_neb_atom_coordinates = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2)); for (i=0; i<(NEB_Num_Images+2); i++){ tmp_neb_atom_coordinates[i] = (double**)malloc(sizeof(double*)*(atomnum+1)); for (j=0; j<(atomnum+1); j++){ tmp_neb_atom_coordinates[i][j] = (double*)malloc(sizeof(double)*20); for (k=0; k<20; k++) neb_atom_coordinates[i][j][k] = 0.0; } } His_neb_atom_coordinates = (double****)malloc(sizeof(double***)*(M_GDIIS_HISTORY+2)); for (m=0; m<(M_GDIIS_HISTORY+2); m++){ His_neb_atom_coordinates[m] = (double***)malloc(sizeof(double**)*(NEB_Num_Images+2)); for (i=0; i<(NEB_Num_Images+2); i++){ His_neb_atom_coordinates[m][i] = (double**)malloc(sizeof(double*)*(atomnum+1)); for (j=0; j<(atomnum+1); j++){ His_neb_atom_coordinates[m][i][j] = (double*)malloc(sizeof(double)*20); for (k=0; k<20; k++) His_neb_atom_coordinates[m][i][j][k] = 0.0; } } } All_Grid_Origin = (double**)malloc(sizeof(double*)*(NEB_Num_Images+2)); for (i=0; i<(NEB_Num_Images+2); i++){ All_Grid_Origin[i] = (double*)malloc(sizeof(double)*4); for (j=0; j<4; j++) All_Grid_Origin[i][j] = 0.0; } Tmp_Grid_Origin = (double**)malloc(sizeof(double*)*(NEB_Num_Images+2)); for (i=0; i<(NEB_Num_Images+2); i++){ Tmp_Grid_Origin[i] = (double*)malloc(sizeof(double)*4); for (j=0; j<4; j++) Tmp_Grid_Origin[i][j] = 0.0; } WhatSpecies_NEB = (int*)malloc(sizeof(int)*(atomnum+1)); Spe_WhatAtom_NEB = (int*)malloc(sizeof(int)*SpeciesNum); SpeName_NEB = (char**)malloc(sizeof(char*)*SpeciesNum); for (i=0; i