From c954418dad3230619cec1ddb12fbaf517b7e094f Mon Sep 17 00:00:00 2001 From: ltimmerman3 Date: Tue, 19 Nov 2024 12:09:21 -0500 Subject: [PATCH] Updating socket commit with modified mlff code --- ChangeLog | 7 +- src/include/isddft.h | 4 +- src/initialization.c | 10 +- src/mlff/covariance_matrix.c | 29 +- src/mlff/covariance_matrix.h | 14 +- src/mlff/linearsys.c | 1589 ----------------- src/mlff/linearsys.h | 59 - src/mlff/soap_descriptor.c | 51 - src/mlff/soap_descriptor.h | 5 - src/mlff/sparc_mlff_interface.c | 1207 ++++++------- src/mlff/sparc_mlff_interface.h | 113 +- src/readfiles.c | 5 +- .../high_accuracy/AlSi_NPTNH_mlff.inpt | 2 +- .../standard/AlSi_NPTNH_mlff.inpt | 2 +- .../high_accuracy/AlSi_mlff_nonorth.inpt | 2 +- .../standard/AlSi_mlff_nonorth.inpt | 2 +- .../high_accuracy/AlSi_mlffflag1.inpt | 2 +- .../standard/AlSi_mlffflag1.inpt | 2 +- .../mlff/C_sqmlff/high_accuracy/C_sqmlff.inpt | 5 +- tests/mlff/C_sqmlff/standard/C_sqmlff.inpt | 3 +- tests/samplescript_cluster | 8 +- 21 files changed, 612 insertions(+), 2509 deletions(-) delete mode 100644 src/mlff/linearsys.c delete mode 100644 src/mlff/linearsys.h diff --git a/ChangeLog b/ChangeLog index 585fab0d..8899d221 100644 --- a/ChangeLog +++ b/ChangeLog @@ -5,14 +5,17 @@ -------------- Nov 18, 2024 Name: Tian Tian, Lucas Timmerman -Changes: (initialization.c, initialization.h, electronicGroundState.c, - include/isddft.h, main.c, main_socket.c, makefile, md.c, relax.c, +Changes: (initialization.c, readfiles.c, initialization.h, electronicGroundState.c, + include/isddft.h, main.c, main_socket.c, makefile, md.c, relax.c, tests/, socket/*.h, socket/*.c, tests/Socket/*, CI workflow, doc) 1. Add socket communication layer and socket submodule 2. Unifying MLFF and DFT single point calculations to Calculate_Properties function 3. Update CI workflow for compiling socket code 4. Update socket mode simple tests 5. Update documentation for socket mode +6. Rewrite of sparc_mlff_interface to allow mlff compatibility with socket and relaxation +7. Removed extraneous variable hnl_file_name +8. Modified .inpt files for mlff tests to account for src changes -------------- Nov 13, 2024 diff --git a/src/include/isddft.h b/src/include/isddft.h index de0fb238..73a322fb 100644 --- a/src/include/isddft.h +++ b/src/include/isddft.h @@ -670,7 +670,7 @@ typedef struct _SPARC_OBJ{ int n_str_max_mlff; int n_train_max_mlff; int mlff_flag; - int last_train_MD_iter; + int last_train_iter; int kernel_typ_MLFF; int descriptor_typ_MLFF; int N_rgrid_MLFF; @@ -688,7 +688,6 @@ typedef struct _SPARC_OBJ{ double xi_3_SOAP; double F_tol_SOAP; double F_rel_scale; - char hnl_file_name[L_STRING]; char mlff_data_folder[L_STRING]; double stress_rel_scale[6]; int MLFF_DFT_fq; @@ -1153,7 +1152,6 @@ typedef struct _SPARC_INPUT_OBJ{ double F_tol_SOAP; double F_rel_scale; double stress_rel_scale[6]; - char hnl_file_name[L_STRING]; char mlff_data_folder[L_STRING]; int MLFF_DFT_fq; diff --git a/src/initialization.c b/src/initialization.c index d56b49cb..82aad807 100644 --- a/src/initialization.c +++ b/src/initialization.c @@ -51,7 +51,8 @@ #define min(x,y) ((x)<(y)?(x):(y)) #define max(x,y) ((x)>(y)?(x):(y)) -#define N_MEMBR 209 +//#define N_MEMBR 209 +#define N_MEMBR 208 /** @@ -1310,8 +1311,7 @@ void SPARC_copy_input(SPARC_OBJ *pSPARC, SPARC_INPUT_OBJ *pSPARC_Input) { pSPARC->xi_3_SOAP = pSPARC_Input->xi_3_SOAP; pSPARC->F_tol_SOAP = pSPARC_Input->F_tol_SOAP; pSPARC->F_rel_scale = pSPARC_Input->F_rel_scale; - strncpy(pSPARC->hnl_file_name, pSPARC_Input->hnl_file_name,sizeof(pSPARC->hnl_file_name)); - strncpy(pSPARC->mlff_data_folder, pSPARC_Input->mlff_data_folder,sizeof(pSPARC->mlff_data_folder)); + strncpy(pSPARC->mlff_data_folder, pSPARC_Input->mlff_data_folder,sizeof(pSPARC->mlff_data_folder)); pSPARC->stress_rel_scale[0] = pSPARC_Input->stress_rel_scale[0]; pSPARC->stress_rel_scale[1] = pSPARC_Input->stress_rel_scale[1]; pSPARC->stress_rel_scale[2] = pSPARC_Input->stress_rel_scale[2]; @@ -3957,7 +3957,6 @@ void write_output_init(SPARC_OBJ *pSPARC) { fprintf(output_fp, "MLFF_FLAG: %d\n", pSPARC->mlff_flag); fprintf(output_fp, "MLFF_INITIAL_STEPS_TRAIN: %d\n", pSPARC->begin_train_steps); - // fprintf(output_fp, "MLFF_hnl_FILE: %s\n", pSPARC->hnl_file_name); if (pSPARC->mlff_flag>1){ fprintf(output_fp, "MLFF_IF_ATOM_DATA_AVAILABLE: %d\n", pSPARC->if_atom_data_available); fprintf(output_fp, "MLFF_MODEL_FOLDER: %s\n", pSPARC->mlff_data_folder); @@ -4326,7 +4325,7 @@ void SPARC_Input_MPI_create(MPI_Datatype *pSPARC_INPUT_MPI) { 1, 1, 1, 1, 1, 1, /* double */ 32, 32, 32, L_STRING, L_STRING, /* char */ - L_STRING, L_STRING, L_STRING, L_STRING, L_STRING, + L_STRING, L_STRING, L_STRING, L_STRING, L_STRING}; // calculating offsets in an architecture independent manner @@ -4551,7 +4550,6 @@ void SPARC_Input_MPI_create(MPI_Datatype *pSPARC_INPUT_MPI) { MPI_Get_address(&sparc_input_tmp.InDensUCubFilename, addr + i++); MPI_Get_address(&sparc_input_tmp.InDensDCubFilename, addr + i++); - MPI_Get_address(&sparc_input_tmp.hnl_file_name, addr + i++); MPI_Get_address(&sparc_input_tmp.mlff_data_folder, addr + i++); for (i = 0; i < N_MEMBR; i++) { disps[i] = addr[i] - base; diff --git a/src/mlff/covariance_matrix.c b/src/mlff/covariance_matrix.c index 546c0bbe..867346f8 100644 --- a/src/mlff/covariance_matrix.c +++ b/src/mlff/covariance_matrix.c @@ -215,7 +215,7 @@ void copy_descriptors(DescriptorObj *desc_str_MLFF, DescriptorObj *desc_str){ /* -add_firstMD function updates the MLFF_Obj by updating design matrix, b vector etc. for the first MD +add_firstDFT function updates the MLFF_Obj by updating design matrix, b vector etc. for the first MD [Input] 1. desc_str: DescriptorObj structure of the first MD @@ -229,7 +229,7 @@ add_firstMD function updates the MLFF_Obj by updating design matrix, b vector et */ -void add_firstMD(DescriptorObj *desc_str, NeighList *nlist, MLFF_Obj *mlff_str, double E, double *F, double *stress_sparc) { +void add_firstDFT(DescriptorObj *desc_str, NeighList *nlist, MLFF_Obj *mlff_str, double E, double *F, double *stress_sparc) { int rank, nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); @@ -1426,31 +1426,6 @@ void remove_train_cols(MLFF_Obj *mlff_str, int col_ID){ } -/* -get_N_r_hnl function calculates the number of grid points in hnl_file - -[Input] -1. pSPARC: SPARC structure -[Output] -1. pSPARC: SPARC structure -*/ -void get_N_r_hnl(SPARC_OBJ *pSPARC){ - int i, j, info; - FILE *fp; - char line[512]; - char a1[512], a2[512], a3[512], a4[512]; - int count1=0, count2=0; - - // fp = fopen("hnl.txt","r"); - fp = fopen(pSPARC->hnl_file_name,"r"); - - fgets(line, sizeof (line), fp); - sscanf(line, "%s%s%s%s", a1, a2, a3, a4); - int N_r = atoi(a4); - pSPARC->N_rgrid_MLFF = N_r; - fclose(fp); -} - // The functions below are outdated. // /* // mlff_kernel function computes the SOAP kernel between two descriptors {X2_str, X3_str} and {X2_tr, X3_tr} diff --git a/src/mlff/covariance_matrix.h b/src/mlff/covariance_matrix.h index c2c0ad18..b19ad59d 100644 --- a/src/mlff/covariance_matrix.h +++ b/src/mlff/covariance_matrix.h @@ -47,7 +47,7 @@ void copy_descriptors(DescriptorObj *desc_str_MLFF, DescriptorObj *desc_str); /* -add_firstMD function updates the MLFF_Obj by updating design matrix, b vector etc. for the first MD +add_firstDFT function updates the MLFF_Obj by updating design matrix, b vector etc. for the first MD [Input] 1. desc_str: DescriptorObj structure of the first MD @@ -59,7 +59,7 @@ add_firstMD function updates the MLFF_Obj by updating design matrix, b vector et [Output] 1. mlff_str: MLFF_Obj structure */ -void add_firstMD(DescriptorObj *desc_str, NeighList *nlist, MLFF_Obj *mlff_str, double E, double* F, double* stress_sparc); +void add_firstDFT(DescriptorObj *desc_str, NeighList *nlist, MLFF_Obj *mlff_str, double E, double* F, double* stress_sparc); /* @@ -128,16 +128,6 @@ remove_train_cols function removes a given local confiugration from the training */ void remove_train_cols(MLFF_Obj *mlff_str, int col_ID); -/* -get_N_r_hnl function calculates the number of grid points in hnl_file - -[Input] -1. pSPARC: SPARC structure -[Output] -1. pSPARC: SPARC structure -*/ -void get_N_r_hnl(SPARC_OBJ *pSPARC); - // double mlff_kernel(int kernel_typ, double **X2_str, double **X3_str, double **X2_tr, double **X3_tr, int atom, int atom_tr, // double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); diff --git a/src/mlff/linearsys.c b/src/mlff/linearsys.c deleted file mode 100644 index 3b87ecc8..00000000 --- a/src/mlff/linearsys.c +++ /dev/null @@ -1,1589 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#ifdef USE_MKL - #define MKL_Complex16 double _Complex - #include -#else - #include - #include -#endif - -#include "tools_mlff.h" -#include "spherical_harmonics.h" -#include "soap_descriptor.h" -#include "descriptor.h" -#include "mlff_types.h" -#include "sparsification.h" -#include "regression.h" -#include "isddft.h" -#include "ddbp_tools.h" -#include "linearsys.h" -#include "sparc_mlff_interface.h" - -#define max(a,b) ((a)>(b)?(a):(b)) -#define min(a,b) ((a)<(b)?(a):(b)) -//#define au2GPa 29421.02648438959 - -/* -mlff_kernel function computes the SOAP kernel between two descriptors {X2_str, X3_str} and {X2_tr, X3_tr} - -[Input] -1. X2_str: pointer to the first descriptor (2-body) -2. X3_str: pointer to the first descriptor (3-body) -3. X2_tr: pointer to the second descriptor (2-body) -4. X3_tr: pointer to the second descriptor (3-body) -5. beta_2: weight to the 2-body term in the kernel -6. beta_3: weight to the 3-body term in the kernel -7. xi_3: exponent in the kernel -8. size_X2: length of the 2-body kernel -9. size_X3: length of the 3-body kernel -[Output] -1. kernel_val: Value of the kernel -*/ - - -double mlff_kernel(int kernel_typ, double **X2_str, double **X3_str, double **X2_tr, double **X3_tr, int atom, int atom_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3) { - double kernel_val; - if (size_X3 > 0){ - if (kernel_typ ==0){// polynomial kernel implemented in VASP MLFF scheme - kernel_val = soap_kernel_polynomial(X2_str[atom], X3_str[atom], X2_tr[atom_tr], X3_tr[atom_tr], beta_2, beta_3, xi_3, size_X2, size_X3); - } else if (kernel_typ==1){// Gaussiam Kernel - kernel_val = soap_kernel_Gaussian(X2_str[atom], X3_str[atom], X2_tr[atom_tr], X3_tr[atom_tr], beta_2, beta_3, xi_3, size_X2, size_X3); - } else{// Laplacian Kernel - kernel_val = soap_kernel_Laplacian(X2_str[atom], X3_str[atom], X2_tr[atom_tr], X3_tr[atom_tr], beta_2, beta_3, xi_3, size_X2, size_X3); - } - if (kernel_val>1.0+1e-5){ - printf("Error in soap kernel evaluation(>1) error %f\n", kernel_val); - exit(1); - } - } else { - if (kernel_typ ==0){// polynomial kernel implemented in VASP MLFF scheme - kernel_val = GMP_kernel_polynomial(X2_str[atom], X2_tr[atom_tr], xi_3, size_X2); - } else if (kernel_typ==1){// Gaussiam Kernel - kernel_val = GMP_kernel_Gaussian(X2_str[atom], X2_tr[atom_tr], xi_3, size_X2); - } else{// Laplacian Kernel - printf("Laplacian kernel not implemented with GMP\n"); - exit(1); - } - if (kernel_val>1.0+1e-5){ - printf("Error in soap kernel evaluation(>1) error %f\n", kernel_val); - exit(1); - } - } - return kernel_val; - -} - -double der_mlff_kernel(int kernel_typ, double ***dX2_str, double ***dX3_str, double **X2_str, double **X3_str, double **X2_tr, double **X3_tr, int atom, int atom_tr, int neighbor, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3) { - if (size_X3 > 0){ - if (kernel_typ ==0) - return der_soap_kernel_polynomial(dX2_str[atom][neighbor], dX3_str[atom][neighbor], X2_str[atom], X3_str[atom], X2_tr[atom_tr], X3_tr[atom_tr], beta_2, beta_3, xi_3, size_X2, size_X3); - else if (kernel_typ==1) - return der_soap_kernel_Gaussian(dX2_str[atom][neighbor], dX3_str[atom][neighbor], X2_str[atom], X3_str[atom], X2_tr[atom_tr], X3_tr[atom_tr], beta_2, beta_3, xi_3, size_X2, size_X3); - else - return der_soap_kernel_Laplacian(dX2_str[atom][neighbor], dX3_str[atom][neighbor], X2_str[atom], X3_str[atom], X2_tr[atom_tr], X3_tr[atom_tr], beta_2, beta_3, xi_3, size_X2, size_X3); - } else { - if (kernel_typ ==0) - return der_GMP_kernel_polynomial(dX2_str[atom][neighbor], X2_str[atom], X2_tr[atom_tr], xi_3, size_X2); - else if (kernel_typ==1) - return der_GMP_kernel_Gaussian(dX2_str[atom][neighbor], X2_str[atom], X2_tr[atom_tr], xi_3, size_X2); - else{ - printf("Laplacian kernel not implemented with GMP\n"); - exit(1); - } - } -} - -double soap_kernel_polynomial(double *X2_str, double *X3_str, double *X2_tr, double *X3_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3) { - - double norm_X3_str, norm_X3_tr, X3_str_temp[size_X3], X3_tr_temp[size_X3], kernel_val; - - norm_X3_str = sqrt(dotProduct(X3_str, X3_str, size_X3)); - norm_X3_tr = sqrt(dotProduct(X3_tr, X3_tr, size_X3)); - for (int i = 0; i0) - temp = -0.5; - else if ((X2_str[i]-X2_tr[i])<0) - temp = 0.5; - else - temp = 0.0; - temp2[i] = temp*exp_temp2; - } - for (int i=0; i 0) - temp = -0.5; - else if ((X3_str[i]-X3_tr[i])<0) - temp = 0.5; - else - temp = 0.0; - temp3[i] = temp*exp_temp3; - } - - der_val = beta_2 * dotProduct(temp2, dX2_str, size_X2) + beta_3 * dotProduct(temp3, dX3_str, size_X3); - - - - // temp = pow(norm_X3_str, -1-xi_3); - - // temp0 = dotProduct(X3_str, X3_tr_temp, size_X3); - - // temp1 = pow(temp0, xi_3-1); - - // const1 = -1 * beta_3* xi_3 * temp * temp1*temp0; - // const2 = beta_3 * temp * norm_X3_str * xi_3 * temp1; - - // der_val = beta_2 * dotProduct(X2_tr, dX2_str, size_X2) + - // const1*dotProduct(X3_str_temp, dX3_str, size_X3) + const2*dotProduct(X3_tr_temp, dX3_str, size_X3); - - return der_val; -} - -double GMP_kernel_polynomial(double *X2_str, double *X2_tr, - double xi_3, int size_X2) { - - double norm_X2_str, norm_X2_tr, X2_str_temp[size_X2], X2_tr_temp[size_X2], kernel_val; - int i; - norm_X2_str = sqrt(dotProduct(X2_str, X2_str, size_X2)); - norm_X2_tr = sqrt(dotProduct(X2_tr, X2_tr, size_X2)); - for (i = 0; iatom_idx_domain, desc_str->atom_idx_domain, desc_str->natom_domain*sizeof(int)); - memcpy(desc_str_MLFF->el_idx_domain, desc_str->el_idx_domain, desc_str->natom_domain*sizeof(int)); - memcpy(desc_str_MLFF->Nneighbors, desc_str->Nneighbors, desc_str->natom_domain*sizeof(int)); - memcpy(desc_str_MLFF->unique_Nneighbors, desc_str->unique_Nneighbors, desc_str->natom_domain*sizeof(int)); - memcpy(desc_str_MLFF->natom_elem, desc_str->natom_elem, desc_str->nelem*sizeof(int)); - for (int i = 0; i < desc_str->natom_domain; i++){ - memcpy(desc_str_MLFF->unique_Nneighbors_elemWise[i], desc_str->unique_Nneighbors_elemWise[i], desc_str->nelem*sizeof(int)); - } - - for (int i = 0; i < desc_str->natom_domain; i++){ - free(desc_str_MLFF->unique_neighborList[i].array); - desc_str_MLFF->unique_neighborList[i].capacity = desc_str->unique_neighborList[i].capacity; - desc_str_MLFF->unique_neighborList[i].array = (int *)malloc(desc_str->unique_neighborList[i].capacity*sizeof(int)); - desc_str_MLFF->unique_neighborList[i].len = desc_str->unique_neighborList[i].len; - for (int k =0; kunique_neighborList[i].len; k++){ - desc_str_MLFF->unique_neighborList[i].array[k] = desc_str->unique_neighborList[i].array[k]; - } - } - - for (int i = 0; i < desc_str->natom_domain; i++){ - memcpy(desc_str_MLFF->X2[i], desc_str->X2[i], desc_str->size_X2*sizeof(double)); - int uniq_natms = desc_str->unique_Nneighbors[i]; - // Large rcuts result in self image in neighbor list - if (uniq_natms == desc_str->natom) uniq_natms -= 1; - for (int j = 0; j < 1+uniq_natms; j++){ - memcpy(desc_str_MLFF->dX2_dX[i][j], desc_str->dX2_dX[i][j], desc_str->size_X2*sizeof(double)); - memcpy(desc_str_MLFF->dX2_dY[i][j], desc_str->dX2_dY[i][j], desc_str->size_X2*sizeof(double)); - memcpy(desc_str_MLFF->dX2_dZ[i][j], desc_str->dX2_dZ[i][j], desc_str->size_X2*sizeof(double)); - } - for (int j = 0; j < 6; j++){ - memcpy(desc_str_MLFF->dX2_dF[i][j], desc_str->dX2_dF[i][j], desc_str->size_X2*sizeof(double)); - } - } - if (desc_str->descriptor_typ < 2) { - for (int i = 0; i < desc_str->natom_domain; i++){ - memcpy(desc_str_MLFF->X3[i], desc_str->X3[i], desc_str->size_X3*sizeof(double)); - int uniq_natms = desc_str->unique_Nneighbors[i]; - // Large rcuts result in self image in neighbor list - if (uniq_natms == desc_str->natom) uniq_natms -= 1; - for (int j = 0; j < 1+uniq_natms; j++){ - memcpy(desc_str_MLFF->dX3_dX[i][j], desc_str->dX3_dX[i][j], desc_str->size_X3*sizeof(double)); - memcpy(desc_str_MLFF->dX3_dY[i][j], desc_str->dX3_dY[i][j], desc_str->size_X3*sizeof(double)); - memcpy(desc_str_MLFF->dX3_dZ[i][j], desc_str->dX3_dZ[i][j], desc_str->size_X3*sizeof(double)); - } - for (int j = 0; j < 6; j++){ - memcpy(desc_str_MLFF->dX3_dF[i][j], desc_str->dX3_dF[i][j], desc_str->size_X3*sizeof(double)); - } - } - } -} - -/* -add_firstMD function updates the MLFF_Obj by updating design matrix, b vector etc. for the first MD - -[Input] -1. desc_str: DescriptorObj structure of the first MD -2. nlist: NeighList strcuture of the first MD -3. mlff_str: MLFF_Obj structure -4. E: energy per atom of first MD structure (Ha/atom) -5. F: atomic foces of first MD structure (Ha/bohr) [ColMajor] -6. stress: stress of first MD structure (GPa) -[Output] -1. mlff_str: MLFF_Obj structure -*/ - - -void add_firstMD(DescriptorObj *desc_str, NeighList *nlist, MLFF_Obj *mlff_str, double E, double *F, double *stress_sparc) { - - int rank, nprocs; - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - - - - int natom = desc_str->natom; - int nelem = desc_str->nelem; - double xi_3 = desc_str->xi_3; - - int size_X2, size_X3; - double beta_2, beta_3; - size_X2 = desc_str->size_X2; - beta_2 = desc_str->beta_2; - beta_3 = desc_str->beta_3; - size_X3 = desc_str->size_X3; - - - int row_idx, col_idx, count, atom_idx, idx; - int iatm_str, istress; - - // calculate mean and std deviation to do the normalization - mlff_str->E_store[mlff_str->E_store_counter] = E; - for (int i = 0; i < mlff_str->stress_len; i++) - mlff_str->stress_store[i][mlff_str->E_store_counter] = stress_sparc[i]; - - for (int i=0; i < 3*natom; i++){ - mlff_str->F_store[mlff_str->F_store_counter+i] = F[i]; - } - mlff_str->E_store_counter += 1; - mlff_str->F_store_counter += 3*natom; - - int kernel_typ = mlff_str->kernel_typ; - mlff_str->mu_E = E; - mlff_str->std_E = 1; - mlff_str->std_F = sqrt(get_variance(F, 3*natom)); - for (int i = 0; i < mlff_str->stress_len; i++){ - mlff_str->std_stress[i] = 1; - } - - // populate bvec - if (rank==0){ - mlff_str->b_no_norm[0] = E; - for (int istress=0; istress < mlff_str->stress_len; istress++){ - mlff_str->b_no_norm[1+istress] = stress_sparc[istress]; - } - - for (int i = 0; i < mlff_str->natom_domain; i++){ - idx = mlff_str->atom_idx_domain[i]; - mlff_str->b_no_norm[3*i+1+mlff_str->stress_len] = F[3*idx]; - mlff_str->b_no_norm[3*i+2+mlff_str->stress_len] = F[3*idx+1]; - mlff_str->b_no_norm[3*i+3+mlff_str->stress_len] = F[3*idx+2]; - - } - } else{ - for (int i = 0; i < mlff_str->natom_domain; i++){ - idx = mlff_str->atom_idx_domain[i]; - mlff_str->b_no_norm[3*i] = F[3*idx]; - mlff_str->b_no_norm[3*i+1] = F[3*idx+1]; - mlff_str->b_no_norm[3*i+2] = F[3*idx+2]; - } - } - - int *cum_natm_elem = (int *)malloc(sizeof(int)*nelem); - - //Stores the starting index wrt SPARC for each element type - sequential - cum_natm_elem[0] = 0; - for (int i = 1; i < nelem; i++){ - cum_natm_elem[i] = cum_natm_elem[i-1]+desc_str->natom_elem[i-1]; - } - - double *X2_gathered, *X3_gathered; - - X2_gathered = (double *) malloc(sizeof(double)*size_X2*desc_str->natom); - X3_gathered = (double *) malloc(sizeof(double)*size_X3*desc_str->natom); - - double *X2_local, *X3_local; - - X2_local = (double *) malloc(sizeof(double)*size_X2*mlff_str->natom_domain); - X3_local = (double *) malloc(sizeof(double)*size_X3*mlff_str->natom_domain); - - for (int i=0; i < mlff_str->natom_domain; i++){ - for (int j=0; j < size_X2; j++){ - X2_local[i*size_X2+j] = desc_str->X2[i][j]; - } - for (int j=0; j < size_X3; j++){ - X3_local[i*size_X3+j] = desc_str->X3[i][j]; - } - } - - int local_natoms[nprocs]; - MPI_Allgather(&mlff_str->natom_domain, 1, MPI_INT, local_natoms, 1, MPI_INT, MPI_COMM_WORLD); - - int recvcounts_X2[nprocs], recvcounts_X3[nprocs], displs_X2[nprocs], displs_X3[nprocs]; - displs_X2[0] = 0; - displs_X3[0] = 0; - for (int i=0; i < nprocs; i++){ - recvcounts_X2[i] = local_natoms[i]*size_X2; - recvcounts_X3[i] = local_natoms[i]*size_X3; - if (i>0){ - displs_X2[i] = displs_X2[i-1]+local_natoms[i-1]*size_X2; - displs_X3[i] = displs_X3[i-1]+local_natoms[i-1]*size_X3; - } - } - - - MPI_Allgatherv(X2_local, size_X2*mlff_str->natom_domain, MPI_DOUBLE, X2_gathered, recvcounts_X2, displs_X2, MPI_DOUBLE, MPI_COMM_WORLD); - MPI_Allgatherv(X3_local, size_X3*mlff_str->natom_domain, MPI_DOUBLE, X3_gathered, recvcounts_X3, displs_X3, MPI_DOUBLE, MPI_COMM_WORLD); - - double **X2_gathered_2D = (double **) malloc(sizeof(double*)*natom); - double **X3_gathered_2D = (double **) malloc(sizeof(double*)*natom); - - for (int i=0; i < natom; i++){ - X2_gathered_2D[i] = (double *) malloc(sizeof(double)*size_X2); - X3_gathered_2D[i] = (double *) malloc(sizeof(double)*size_X3); - for (int j=0; j < size_X2; j++){ - X2_gathered_2D[i][j] = X2_gathered[i*size_X2+j]; - } - for (int j=0; j < size_X3; j++){ - X3_gathered_2D[i][j] = X3_gathered[i*size_X3+j]; - } - } - - // if (rank == 0) { - // fprintf(mlff_str->fp_mlff, "X2\n"); - // for (int i=0; i < natom; i++){ - // for (int j=0; j < size_X2; j++){ - // fprintf(mlff_str->fp_mlff, "%lf ", X2_gathered_2D[i][j]); - // } - // fprintf(mlff_str->fp_mlff, "\n"); - // } - // fprintf(mlff_str->fp_mlff, "X3\n"); - // for (int i=0; i < natom; i++){ - // for (int j=0; j < size_X3; j++){ - // fprintf(mlff_str->fp_mlff, "%lf ", X3_gathered_2D[i][j]); - // } - // fprintf(mlff_str->fp_mlff, "\n"); - // } - // } - - dyArray *highrank_ID_descriptors = (dyArray *) malloc(sizeof(dyArray)*nelem); - for (int i=0; i natom_elem[i] - 500; - init_dyarray(&highrank_ID_descriptors[i]); - mlff_CUR_sparsify(kernel_typ, &X2_gathered_2D[cum_natm_elem[i]], &X3_gathered_2D[cum_natm_elem[i]], - desc_str->natom_elem[i], size_X2, size_X3, beta_2, beta_3, xi_3, &highrank_ID_descriptors[i], N_low_min); - } - - for (int i = 0; i < nelem; i++){ - mlff_str->natm_train_elemwise[i] = (highrank_ID_descriptors[i]).len; - } - - count=0; - for (int i = 0; i < nelem; i++){ - for (int j = 0; j < (highrank_ID_descriptors[i]).len; j++){ - mlff_str->natm_typ_train[count] = i; - for(int jj = 0; jj < size_X2; jj++){ - mlff_str->X2_traindataset[count][jj] = X2_gathered_2D[cum_natm_elem[i]+(highrank_ID_descriptors[i]).array[j]][jj]; - } - for(int jj = 0; jj < size_X3; jj++){ - mlff_str->X3_traindataset[count][jj] = X3_gathered_2D[cum_natm_elem[i]+(highrank_ID_descriptors[i]).array[j]][jj]; - } - count++; - } - } - - initialize_descriptor(mlff_str->descriptor_strdataset, mlff_str, nlist); - copy_descriptors(mlff_str->descriptor_strdataset, desc_str); - - - mlff_str->n_str = 1; - mlff_str->natm_train_total = count; - if (rank == 0) fprintf(mlff_str->fp_mlff, "natm_train_total: %d\n", mlff_str->natm_train_total); - if (rank==0){ - mlff_str->n_rows = 3*desc_str->natom_domain + 1 + mlff_str->stress_len; - } else{ - mlff_str->n_rows = 3*desc_str->natom_domain; - } - - - mlff_str->n_cols = count; - - - int *cum_natm_elem1 = (int *)malloc(sizeof(int)*nelem); - cum_natm_elem1[0] = 0; - for (int i = 1; i < nelem; i++){ - cum_natm_elem1[i] = cum_natm_elem1[i-1] + mlff_str->natm_train_elemwise[i-1]; - } - - - double *K_train_local = (double *) malloc(sizeof(double)*mlff_str->n_cols*(3*natom+1+mlff_str->stress_len)); // row major; - for (int i=0; i < mlff_str->n_cols*(3*natom+1+mlff_str->stress_len); i++){ - K_train_local[i] = 0.0; - } - - // Energy term to populate K_train matrix - int el_type; - double temp = (1.0/ (double) desc_str->natom); - for (int i=0; i < desc_str->natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < (highrank_ID_descriptors[el_type]).len; iatm_el_tr++){ - col_idx = cum_natm_elem1[el_type] + iatm_el_tr; - K_train_local[col_idx] += temp * mlff_kernel(kernel_typ, desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, - beta_2, beta_3, xi_3, size_X2, size_X3); - } - } - - // Force term to populate K_train matrix - int atm_idx; - for (int i=0; i natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - atm_idx = desc_str->atom_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < (highrank_ID_descriptors[el_type]).len; iatm_el_tr++){ - col_idx = cum_natm_elem1[el_type] + iatm_el_tr; - row_idx = 3*atm_idx+1; - - K_train_local[row_idx*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dX, desc_str->dX3_dX, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+1)*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dY, desc_str->dX3_dY, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+2)*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dZ, desc_str->dX3_dZ, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - // unique neighbor list never contains self image - for (int neighs =0; neighs < desc_str->unique_Nneighbors[i]; neighs++){ - row_idx = 3*nlist->unique_neighborList[i].array[neighs]+1; - K_train_local[row_idx*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dX, desc_str->dX3_dX, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+1)*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dY, desc_str->dX3_dY, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+2)*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dZ, desc_str->dX3_dZ, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - } - - } - } - - double volume = desc_str->cell_measure; - - - - for (int i=0; i < desc_str->natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - atm_idx = desc_str->atom_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < (highrank_ID_descriptors[el_type]).len; iatm_el_tr++){ - col_idx = cum_natm_elem1[el_type] + iatm_el_tr; - for (int istress = 0; istress < mlff_str->stress_len; istress++){ - row_idx = 3*desc_str->natom+1+istress; - K_train_local[(row_idx)*mlff_str->n_cols + col_idx] += (1.0/volume)* - der_mlff_kernel(kernel_typ, desc_str->dX2_dF, desc_str->dX3_dF, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, istress, - beta_2, beta_3, xi_3, size_X2, size_X3); - } - } - } - - - double *K_train_assembled; - K_train_assembled = (double *)malloc(sizeof(double)*mlff_str->n_cols*(3*natom+1+mlff_str->stress_len)); - - MPI_Allreduce(K_train_local, K_train_assembled, mlff_str->n_cols*(3*natom+1+mlff_str->stress_len), MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - - int r_idx; - if (rank==0){ - for (int i=0; i < mlff_str->n_cols; i++){ - mlff_str->K_train[0][i] = K_train_assembled[i]; - } - for (int istress =0; istress < mlff_str->stress_len; istress++){ - r_idx = 3*desc_str->natom+1+istress; - for (int i=0; i < mlff_str->n_cols; i++){ - mlff_str->K_train[1+istress][i] = K_train_assembled[r_idx*mlff_str->n_cols + i]; - } - } - } - - - - for (int i = 0; i < desc_str->natom_domain; i++){ - atm_idx = desc_str->atom_idx_domain[i]; - r_idx = 3*atm_idx+1; - - for (int j=0; j < mlff_str->n_cols; j++){ - if (rank==0){ - mlff_str->K_train[3*i+1+mlff_str->stress_len][j] = K_train_assembled[r_idx*mlff_str->n_cols + j]; - mlff_str->K_train[3*i+2+mlff_str->stress_len][j] = K_train_assembled[(1+r_idx)*mlff_str->n_cols + j]; - mlff_str->K_train[3*i+3+mlff_str->stress_len][j] = K_train_assembled[(2+r_idx)*mlff_str->n_cols + j]; - } else{ - mlff_str->K_train[3*i][j] = K_train_assembled[r_idx*mlff_str->n_cols + j]; - mlff_str->K_train[3*i+1][j] = K_train_assembled[(1+r_idx)*mlff_str->n_cols + j]; - mlff_str->K_train[3*i+2][j] = K_train_assembled[(2+r_idx)*mlff_str->n_cols + j]; - } - - } - } - - for (int i=0; i natom; - int nelem = desc_str->nelem ; - int size_X2 = desc_str->size_X2; - int size_X3 = desc_str->size_X3; - double beta_2 = desc_str->beta_2; - double beta_3 = desc_str->beta_3; - double xi_3 = desc_str->xi_3; - - int kernel_typ = mlff_str->kernel_typ; - - int rank, nprocs; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); - - mlff_str->E_store[mlff_str->E_store_counter] = E; - for (int i = 0; i < mlff_str->stress_len; i++) - mlff_str->stress_store[i][mlff_str->E_store_counter] = stress_sparc[i]; - for (int i=0; i < 3*natom; i++){ - mlff_str->F_store[mlff_str->F_store_counter+i] = F[i]; - } - mlff_str->E_store_counter += 1; - mlff_str->F_store_counter += 3*natom; - - if (mlff_str->n_str > 1) { - mlff_str->mu_E = get_mean(mlff_str->E_store, mlff_str->E_store_counter); - mlff_str->std_E = sqrt(get_variance(mlff_str->E_store, mlff_str->E_store_counter)); - mlff_str->std_F = sqrt(get_variance(mlff_str->F_store, mlff_str->F_store_counter)); - for (int i = 0; i < mlff_str->stress_len; i++){ - mlff_str->std_stress[i] = sqrt(get_variance(mlff_str->stress_store[i], mlff_str->E_store_counter)); - } - } else { - mlff_str->mu_E = mlff_str->E_store[0]; - mlff_str->std_E = 1.0; - mlff_str->std_F = sqrt(get_variance(mlff_str->F_store, mlff_str->F_store_counter)); - for (int i = 0; i < mlff_str->stress_len; i++){ - mlff_str->std_stress[i] = 1.0; - } - } - - int idx; - if (rank==0){ - mlff_str->b_no_norm[mlff_str->n_rows] = E; - for (int istress=0; istress < mlff_str->stress_len; istress++){ - mlff_str->b_no_norm[mlff_str->n_rows+1+istress] = stress_sparc[istress]; - } - for (int i = 0; i < desc_str->natom_domain; i++){ - idx = mlff_str->atom_idx_domain[i]; - mlff_str->b_no_norm[mlff_str->n_rows+3*i+1+mlff_str->stress_len] = F[3*idx]; - mlff_str->b_no_norm[mlff_str->n_rows+3*i+2+mlff_str->stress_len] = F[3*idx+1]; - mlff_str->b_no_norm[mlff_str->n_rows+3*i+3+mlff_str->stress_len] = F[3*idx+2]; - } - } else{ - for (int i = 0; i < desc_str->natom_domain; i++){ - idx = mlff_str->atom_idx_domain[i]; - mlff_str->b_no_norm[mlff_str->n_rows+3*i] = F[3*idx]; - mlff_str->b_no_norm[mlff_str->n_rows+3*i+1] = F[3*idx+1]; - mlff_str->b_no_norm[mlff_str->n_rows+3*i+2] = F[3*idx+2]; - } - } - - initialize_descriptor(mlff_str->descriptor_strdataset+mlff_str->n_str, mlff_str, nlist); - copy_descriptors(mlff_str->descriptor_strdataset+mlff_str->n_str, desc_str); - - int *cum_natm_elem = (int *)malloc(sizeof(int)*nelem); - cum_natm_elem[0] = 0; - for (int i = 1; i < nelem; i++){ - cum_natm_elem[i] = cum_natm_elem[i-1] + desc_str->natom_elem[i-1]; - } - - int **cum_natm_ele_cols = (int **)malloc(sizeof(int*)*nelem); - for (int i = 0; i natm_train_elemwise[i]); - } - - for (int i = 0; i < nelem; i++){ - int count=0; - for (int j=0; j < mlff_str->natm_train_total; j++){ - if (mlff_str->natm_typ_train[j] == i){ - cum_natm_ele_cols[i][count] = j; - count++; - } - } - } - - double temp = (1.0/ (double) natom); - double *K_train_local = (double *) malloc(sizeof(double)*mlff_str->n_cols*(3*natom+1+mlff_str->stress_len)); // row major - for (int i=0; i < mlff_str->n_cols*(3*natom+1+mlff_str->stress_len); i++){ - K_train_local[i] = 0.0; - } - - int el_type; - for (int i=0; i < desc_str->natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < mlff_str->natm_train_elemwise[el_type]; iatm_el_tr++){ - col_idx = cum_natm_ele_cols[el_type][iatm_el_tr]; - K_train_local[col_idx] += temp * mlff_kernel(kernel_typ, desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, - beta_2, beta_3, xi_3, size_X2, size_X3); - } - } - - int atm_idx; - for (int i=0; i natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - atm_idx = desc_str->atom_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < mlff_str->natm_train_elemwise[el_type]; iatm_el_tr++){ - col_idx = cum_natm_ele_cols[el_type][iatm_el_tr]; - row_idx = 3*atm_idx+1; - - K_train_local[row_idx*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dX, desc_str->dX3_dX, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+1)*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dY, desc_str->dX3_dY, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+2)*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dZ, desc_str->dX3_dZ, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - for (int neighs =0; neighs < desc_str->unique_Nneighbors[i]; neighs++){ - row_idx = 3*nlist->unique_neighborList[i].array[neighs]+1; - K_train_local[row_idx*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dX, desc_str->dX3_dX, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+1)*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dY, desc_str->dX3_dY, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+2)*mlff_str->n_cols + col_idx] += - der_mlff_kernel(kernel_typ, desc_str->dX2_dZ, desc_str->dX3_dZ, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - } - - } - } - - - double volume = desc_str->cell_measure; - - for (int i=0; i < desc_str->natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - atm_idx = desc_str->atom_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < mlff_str->natm_train_elemwise[el_type]; iatm_el_tr++){ - col_idx = cum_natm_ele_cols[el_type][iatm_el_tr]; - for (int istress = 0; istress < mlff_str->stress_len; istress++){ - row_idx = 3*desc_str->natom+1+istress; - K_train_local[(row_idx)*mlff_str->n_cols + col_idx] += (1.0/volume)* - der_mlff_kernel(kernel_typ, desc_str->dX2_dF, desc_str->dX3_dF, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, istress, - beta_2, beta_3, xi_3, size_X2, size_X3); - - } - } - } - - - double *K_train_assembled; - K_train_assembled = (double *)malloc(sizeof(double)*mlff_str->n_cols*(3*natom+1+mlff_str->stress_len)); - - MPI_Allreduce(K_train_local, K_train_assembled, mlff_str->n_cols*(3*natom+1+mlff_str->stress_len), MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - - int r_idx; - - if (rank==0){ - for (int i=0; i < mlff_str->n_cols; i++){ - mlff_str->K_train[mlff_str->n_rows][i] = K_train_assembled[i]; - } - for (int istress =0; istress < mlff_str->stress_len; istress++){ - r_idx = 3*desc_str->natom+1+istress; - for (int i=0; i < mlff_str->n_cols; i++){ - mlff_str->K_train[mlff_str->n_rows+1+istress][i] = K_train_assembled[r_idx*mlff_str->n_cols + i]; - } - } - } - - - for (int i = 0; i < desc_str->natom_domain; i++){ - atm_idx = desc_str->atom_idx_domain[i]; - r_idx = 3*atm_idx+1; - - for (int j=0; j < mlff_str->n_cols; j++){ - if (rank==0){ - mlff_str->K_train[mlff_str->n_rows+3*i+1+mlff_str->stress_len][j] = K_train_assembled[r_idx*mlff_str->n_cols + j]; - mlff_str->K_train[mlff_str->n_rows+3*i+2+mlff_str->stress_len][j] = K_train_assembled[(1+r_idx)*mlff_str->n_cols + j]; - mlff_str->K_train[mlff_str->n_rows+3*i+3+mlff_str->stress_len][j] = K_train_assembled[(2+r_idx)*mlff_str->n_cols + j]; - } else{ - mlff_str->K_train[mlff_str->n_rows+3*i][j] = K_train_assembled[r_idx*mlff_str->n_cols + j]; - mlff_str->K_train[mlff_str->n_rows+3*i+1][j] = K_train_assembled[(1+r_idx)*mlff_str->n_cols + j]; - mlff_str->K_train[mlff_str->n_rows+3*i+2][j] = K_train_assembled[(2+r_idx)*mlff_str->n_cols + j]; - } - - } - } - free(K_train_local); - free(K_train_assembled); - - - // updating other MLFF parameters such as number of structures, number of training environment, element typ of training env - mlff_str->n_str += 1; - if (rank==0){ - mlff_str->n_rows += 3*desc_str->natom_domain + 1 + mlff_str->stress_len; - } else { - mlff_str->n_rows += 3*desc_str->natom_domain; - } - free(cum_natm_elem); - for (int i = 0; i natom; - int nelem = desc_str->nelem; - int size_X2 = desc_str->size_X2; - int size_X3 = desc_str->size_X3; - double beta_2 = desc_str->beta_2; - double beta_3 = desc_str->beta_3; - double xi_3 = desc_str->xi_3; - stress_scale = (double*) malloc(mlff_str->stress_len*sizeof(double)); - - - int rank, nproc; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &nproc); - - int kernel_typ = mlff_str->kernel_typ; - - E_scale = mlff_str->E_scale; - F_scale = mlff_str->F_scale * mlff_str->relative_scale_F; - - for (int i=0; i < mlff_str->stress_len; i++){ - stress_scale[i] = mlff_str->stress_scale[i] * mlff_str->relative_scale_stress[i]; - - } - - int *cum_natm_elem = (int *)malloc(sizeof(int)*nelem); - cum_natm_elem[0] = 0; - for (int i = 1; i < nelem; i++){ - cum_natm_elem[i] = cum_natm_elem[i-1] + desc_str->natom_elem[i-1]; - } - - int **cum_natm_ele_cols = (int **)malloc(sizeof(int*)*nelem); - for (int i = 0; i natm_train_elemwise[i]); - } - - for (int i = 0; i < nelem; i++){ - int count=0; - for (int j=0; j < mlff_str->natm_train_total; j++){ - if (mlff_str->natm_typ_train[j] == i){ - cum_natm_ele_cols[i][count] = j; - count++; - } - } - } - - double temp = (1.0/ (double) desc_str->natom); - double *K_train_local = (double *) malloc(sizeof(double)*mlff_str->n_cols*(3*natom+1+mlff_str->stress_len)); // row major; - for (int i=0; i < mlff_str->n_cols*(3*natom+1+mlff_str->stress_len); i++){ - K_train_local[i] = 0.0; - } - - int el_type; - for (int i=0; i natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < mlff_str->natm_train_elemwise[el_type]; iatm_el_tr++){ - col_idx = cum_natm_ele_cols[el_type][iatm_el_tr]; - K_train_local[col_idx] += temp * mlff_kernel(kernel_typ, desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, - beta_2, beta_3, xi_3, size_X2, size_X3); - } - } - - int atm_idx; - for (int i=0; i natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - atm_idx = desc_str->atom_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < mlff_str->natm_train_elemwise[el_type]; iatm_el_tr++){ - col_idx = cum_natm_ele_cols[el_type][iatm_el_tr]; - row_idx = 3*atm_idx+1; - - K_train_local[row_idx*mlff_str->n_cols + col_idx] += F_scale* - der_mlff_kernel(kernel_typ, desc_str->dX2_dX, desc_str->dX3_dX, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+1)*mlff_str->n_cols + col_idx] += F_scale* - der_mlff_kernel(kernel_typ, desc_str->dX2_dY, desc_str->dX3_dY, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+2)*mlff_str->n_cols + col_idx] += F_scale* - der_mlff_kernel(kernel_typ, desc_str->dX2_dZ, desc_str->dX3_dZ, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 0, - beta_2, beta_3, xi_3, size_X2, size_X3); - for (int neighs =0; neighs < desc_str->unique_Nneighbors[i]; neighs++){ - row_idx = 3*nlist->unique_neighborList[i].array[neighs]+1; - K_train_local[row_idx*mlff_str->n_cols + col_idx] += F_scale* - der_mlff_kernel(kernel_typ, desc_str->dX2_dX, desc_str->dX3_dX, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+1)*mlff_str->n_cols + col_idx] += F_scale* - der_mlff_kernel(kernel_typ, desc_str->dX2_dY, desc_str->dX3_dY, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - - K_train_local[(row_idx+2)*mlff_str->n_cols + col_idx] += F_scale* - der_mlff_kernel(kernel_typ, desc_str->dX2_dZ, desc_str->dX3_dZ, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, 1+neighs, - beta_2, beta_3, xi_3, size_X2, size_X3); - } - - } - } - - - double volume = desc_str->cell_measure; - - for (int i=0; i < desc_str->natom_domain; i++){ - el_type = desc_str->el_idx_domain[i]; - atm_idx = desc_str->atom_idx_domain[i]; - for (int iatm_el_tr = 0; iatm_el_tr < mlff_str->natm_train_elemwise[el_type]; iatm_el_tr++){ - col_idx = cum_natm_ele_cols[el_type][iatm_el_tr]; - for (int istress = 0; istress < mlff_str->stress_len; istress++){ - row_idx = 3*desc_str->natom+1+istress; - K_train_local[(row_idx)*mlff_str->n_cols + col_idx] += (1.0/volume)*stress_scale[istress]* - der_mlff_kernel(kernel_typ, desc_str->dX2_dF, desc_str->dX3_dF, - desc_str->X2, desc_str->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, col_idx, istress, - beta_2, beta_3, xi_3, size_X2, size_X3); - - } - } - } - - - - - double *K_train_assembled; - K_train_assembled = (double *)malloc(sizeof(double)*mlff_str->n_cols*(3*natom+1+mlff_str->stress_len)); - - - MPI_Allreduce(K_train_local, K_train_assembled, mlff_str->n_cols*(3*natom+1+mlff_str->stress_len), MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - - int r_idx; - if (rank==0){ - for (int i=0; i < mlff_str->n_cols; i++){ - K_predict[0][i] = K_train_assembled[i]; - } - for (int istress =0; istress < mlff_str->stress_len; istress++){ - r_idx = 3*desc_str->natom+1+istress; - for (int i=0; i < mlff_str->n_cols; i++){ - K_predict[1+istress][i] = K_train_assembled[r_idx*mlff_str->n_cols + i]; - } - } - } - - - - for (int i = 0; i < desc_str->natom_domain; i++){ - atm_idx = desc_str->atom_idx_domain[i]; - r_idx = 3*atm_idx+1; - - for (int j=0; j < mlff_str->n_cols; j++){ - if (rank==0){ - K_predict[3*i+1+mlff_str->stress_len][j] = K_train_assembled[r_idx*mlff_str->n_cols + j]; - K_predict[3*i+2+mlff_str->stress_len][j] = K_train_assembled[(1+r_idx)*mlff_str->n_cols + j]; - K_predict[3*i+3+mlff_str->stress_len][j] = K_train_assembled[(2+r_idx)*mlff_str->n_cols + j]; - } else{ - K_predict[3*i][j] = K_train_assembled[r_idx*mlff_str->n_cols + j]; - K_predict[3*i+1][j] = K_train_assembled[(1+r_idx)*mlff_str->n_cols + j]; - K_predict[3*i+2][j] = K_train_assembled[(2+r_idx)*mlff_str->n_cols + j]; - } - - } - } - - - - - - free(K_train_local); - free(K_train_assembled); - - - - free(cum_natm_elem); - for (int i = 0; i nelem; - int natom = mlff_str->natom; - int size_X2 = mlff_str->size_X2; - int size_X3 = mlff_str->size_X3; - int kernel_typ = mlff_str->kernel_typ; - double beta_2 = mlff_str->beta_2; - double beta_3 = mlff_str->beta_3; - double xi_3 = mlff_str->xi_3; - - int rank, nprocs; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &nprocs); - - // copy the X2 and X3 into train history data - - for(int j = 0; j < size_X2; j++){ - mlff_str->X2_traindataset[mlff_str->n_cols][j] = X2[j]; - } - for(int j = 0; j < size_X3; j++){ - mlff_str->X3_traindataset[mlff_str->n_cols][j] = X3[j]; - } - - double *k_local; - int total_rows = 0; - for (int istr = 0; istr < mlff_str->n_str; istr++){ - total_rows = total_rows + 1 + 3*(mlff_str->descriptor_strdataset+istr)->natom + mlff_str->stress_len; - } - k_local = (double *)malloc(sizeof(double)*total_rows); - for (int i=0; i< total_rows; i++){ - k_local[i] = 0.0; - } - - row_idx=0; - int el_type; - double temp; - for (int istr = 0; istr < mlff_str->n_str; istr++){ - temp = 1.0/ ((double) (mlff_str->descriptor_strdataset+istr)->natom); - for (int i=0; i <(mlff_str->descriptor_strdataset+istr)->natom_domain; i++){ - el_type = (mlff_str->descriptor_strdataset+istr)->el_idx_domain[i]; - if (el_type==elem_typ){ - k_local[row_idx] += temp* - mlff_kernel(kernel_typ, (mlff_str->descriptor_strdataset+istr)->X2, (mlff_str->descriptor_strdataset+istr)->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, i, mlff_str->n_cols, beta_2, beta_3, xi_3, size_X2, size_X3); - } - } - row_idx += 1+3*(mlff_str->descriptor_strdataset+istr)->natom+mlff_str->stress_len; - } - - // Force term to populate K_train matrix - row_idx = 0; - int row_index_F; - for (int istr = 0; istr < mlff_str->n_str; istr++){ - for (int iatm_str = 0; iatm_str < (mlff_str->descriptor_strdataset+istr)->natom_domain; iatm_str++){ - el_type = (mlff_str->descriptor_strdataset+istr)->el_idx_domain[iatm_str]; - if (el_type==elem_typ){ - atom_idx = (mlff_str->descriptor_strdataset+istr)->atom_idx_domain[iatm_str]; - row_index_F = row_idx + 3*atom_idx+1; - // x-component (w.r.t itself) "because an atom is not considered it's neighbour hence dealt outside the neighs loop" - k_local[row_index_F] += - der_mlff_kernel(kernel_typ, (mlff_str->descriptor_strdataset+istr)->dX2_dX, (mlff_str->descriptor_strdataset+istr)->dX3_dX, - (mlff_str->descriptor_strdataset+istr)->X2, (mlff_str->descriptor_strdataset+istr)->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, iatm_str, mlff_str->n_cols, 0, beta_2, beta_3, xi_3, size_X2, size_X3); - // y-component (w.r.t itself) "because an atom is not considered it's neighbour hence dealt outside the neighs loop" - k_local[row_index_F+1] += - der_mlff_kernel(kernel_typ, (mlff_str->descriptor_strdataset+istr)->dX2_dY, (mlff_str->descriptor_strdataset+istr)->dX3_dY, - (mlff_str->descriptor_strdataset+istr)->X2, (mlff_str->descriptor_strdataset+istr)->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, iatm_str, mlff_str->n_cols, 0, beta_2, beta_3, xi_3, size_X2, size_X3); - // z-component (w.r.t itself) "because an atom is not considered it's neighbour hence dealt outside the neighs loop" - k_local[row_index_F+2] += - der_mlff_kernel(kernel_typ, (mlff_str->descriptor_strdataset+istr)->dX2_dZ, (mlff_str->descriptor_strdataset+istr)->dX3_dZ, - (mlff_str->descriptor_strdataset+istr)->X2, (mlff_str->descriptor_strdataset+istr)->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, iatm_str, mlff_str->n_cols, 0, beta_2, beta_3, xi_3, size_X2, size_X3); - for (int neighs = 0; neighs < (mlff_str->descriptor_strdataset+istr)->unique_Nneighbors[iatm_str]; neighs++){ - row_index_F = row_idx + 3*(mlff_str->descriptor_strdataset+istr)->unique_neighborList[iatm_str].array[neighs]+1; - k_local[row_index_F] += - der_mlff_kernel(kernel_typ, (mlff_str->descriptor_strdataset+istr)->dX2_dX, (mlff_str->descriptor_strdataset+istr)->dX3_dX, - (mlff_str->descriptor_strdataset+istr)->X2, (mlff_str->descriptor_strdataset+istr)->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, iatm_str, mlff_str->n_cols, 1+neighs, beta_2, beta_3, xi_3, size_X2, size_X3); - // y-component (w.r.t neighs neighbour) - k_local[row_index_F+1] += - der_mlff_kernel(kernel_typ, (mlff_str->descriptor_strdataset+istr)->dX2_dY, (mlff_str->descriptor_strdataset+istr)->dX3_dY, - (mlff_str->descriptor_strdataset+istr)->X2, (mlff_str->descriptor_strdataset+istr)->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, iatm_str, mlff_str->n_cols, 1+neighs, beta_2, beta_3, xi_3, size_X2, size_X3); - // z-component (w.r.t neighs neighbour) - k_local[row_index_F+2] += - der_mlff_kernel(kernel_typ, (mlff_str->descriptor_strdataset+istr)->dX2_dZ, (mlff_str->descriptor_strdataset+istr)->dX3_dZ, - (mlff_str->descriptor_strdataset+istr)->X2, (mlff_str->descriptor_strdataset+istr)->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, iatm_str, mlff_str->n_cols, 1+neighs, beta_2, beta_3, xi_3, size_X2, size_X3); - } - } - } - row_idx += 1+3*(mlff_str->descriptor_strdataset+istr)->natom+mlff_str->stress_len; - } - - int row_index_stress; - double volume; - row_idx = 0; - for (int istr = 0; istr < mlff_str->n_str; istr++){ - volume = (mlff_str->descriptor_strdataset+istr)->cell_measure; - for (int iatm_str = 0; iatm_str < (mlff_str->descriptor_strdataset+istr)->natom_domain; iatm_str++){ - el_type = (mlff_str->descriptor_strdataset+istr)->el_idx_domain[iatm_str]; - if (el_type==elem_typ){ - for (int istress = 0; istress < mlff_str->stress_len; istress++){ - row_index_stress = row_idx + 3*(mlff_str->descriptor_strdataset+istr)->natom + 1 + istress; - k_local[row_index_stress] += (1.0/volume)* - der_mlff_kernel(kernel_typ, (mlff_str->descriptor_strdataset+istr)->dX2_dF, (mlff_str->descriptor_strdataset+istr)->dX3_dF, - (mlff_str->descriptor_strdataset+istr)->X2, (mlff_str->descriptor_strdataset+istr)->X3, - mlff_str->X2_traindataset, mlff_str->X3_traindataset, iatm_str, mlff_str->n_cols, istress, beta_2, beta_3, xi_3, size_X2, size_X3); - } - } - } - row_idx += 1+3*(mlff_str->descriptor_strdataset+istr)->natom+mlff_str->stress_len; - } - - double *K_train_assembled; - K_train_assembled = (double *)malloc(sizeof(double)*total_rows); - - MPI_Allreduce(k_local, K_train_assembled, total_rows, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - - int temp_idx1 = 0, temp_idx2 = 0; - if (rank==0){ - for (int i=0; i < mlff_str->n_str; i++){ - - mlff_str->K_train[temp_idx1][mlff_str->n_cols] = K_train_assembled[temp_idx2]; - for (int istress =0; istress < mlff_str->stress_len; istress++){ - mlff_str->K_train[temp_idx1+1+istress][mlff_str->n_cols] = K_train_assembled[temp_idx2 + 3*((mlff_str->descriptor_strdataset+i)->natom)+1+istress]; - } - temp_idx1 = temp_idx1 + 1 + mlff_str->stress_len + 3*((mlff_str->descriptor_strdataset+i)->natom_domain); - temp_idx2 = temp_idx2 + 1 + mlff_str->stress_len + 3*((mlff_str->descriptor_strdataset+i)->natom); - } - } - - int r_idx, temp_idx, atm_idx; - temp_idx=0; - temp_idx1 = 0; - for (int s=0; s < mlff_str->n_str; s++){ - for (int i = 0; i < (mlff_str->descriptor_strdataset+s)->natom_domain; i++){ - atm_idx = (mlff_str->descriptor_strdataset+s)->atom_idx_domain[i]; - r_idx = 3*atm_idx+1; - if (rank==0){ - mlff_str->K_train[temp_idx+3*i+1+mlff_str->stress_len][mlff_str->n_cols] = K_train_assembled[temp_idx1+r_idx]; - mlff_str->K_train[temp_idx+3*i+2+mlff_str->stress_len][mlff_str->n_cols] = K_train_assembled[temp_idx1+r_idx+1]; - mlff_str->K_train[temp_idx+3*i+3+mlff_str->stress_len][mlff_str->n_cols] = K_train_assembled[temp_idx1+r_idx+2]; - } else { - mlff_str->K_train[temp_idx+3*i][mlff_str->n_cols] = K_train_assembled[temp_idx1+r_idx]; - mlff_str->K_train[temp_idx+3*i+1][mlff_str->n_cols] = K_train_assembled[temp_idx1+r_idx+1]; - mlff_str->K_train[temp_idx+3*i+2][mlff_str->n_cols] = K_train_assembled[temp_idx1+r_idx+2]; - } - } - if (rank==0){ - temp_idx = temp_idx + 1 + 3*((mlff_str->descriptor_strdataset+s)->natom_domain) +mlff_str->stress_len; - } else { - temp_idx = temp_idx + 3*((mlff_str->descriptor_strdataset+s)->natom_domain); - } - temp_idx1 = temp_idx1 + 1 + mlff_str->stress_len + 3*((mlff_str->descriptor_strdataset+s)->natom); - } - - // updating other MLFF parameters such as number of structures, number of training environment, element typ of training env - mlff_str->natm_typ_train[mlff_str->n_cols] = elem_typ; - mlff_str->natm_train_total += 1; - mlff_str->n_cols += 1; - mlff_str->natm_train_elemwise[elem_typ] += 1; - free(K_train_assembled); - free(k_local); - -} - -/* -remove_str_rows function removes a given reference structure from the training dataset - -[Input] -1. mlff_str: MLFF_Obj structure -2. str_ID: ID of the reference structure in training dataset -[Output] -1. mlff_str: MLFF_Obj structure -*/ - -void remove_str_rows(MLFF_Obj *mlff_str, int str_ID){ - int start_idx = 0, end_idx = 0, i, j, istr, rows_to_delete; - - for (i = 0; i < str_ID; i++){ - start_idx += 3*(mlff_str->descriptor_strdataset+i)->natom + 1; - } - end_idx = start_idx + 3*(mlff_str->descriptor_strdataset+str_ID)->natom + 1; - - rows_to_delete = end_idx - start_idx; - - - for (istr = str_ID + 1; istr < mlff_str->n_str; istr++){ - copy_descriptors(mlff_str->descriptor_strdataset + istr - 1, mlff_str->descriptor_strdataset + istr ); - } - delete_descriptor(mlff_str->descriptor_strdataset + mlff_str->n_str); - mlff_str->n_str = mlff_str->n_str - 1; - - for (i = start_idx; i < end_idx; i++) { - mlff_str->b_no_norm[i] = mlff_str->b_no_norm[i + rows_to_delete]; - for (j = 0; j < mlff_str->n_cols; j++){ - mlff_str->K_train[i][j] = mlff_str->K_train[i + rows_to_delete][j]; - } - } - - for (i = mlff_str->n_rows - rows_to_delete; i < mlff_str->n_rows; i++) { - mlff_str->b_no_norm[i] = 0.0; - for (j = 0; j < mlff_str->n_cols; j++){ - mlff_str->K_train[i][j] = 0.0; - } - } - mlff_str->n_rows = mlff_str->n_rows - rows_to_delete; -} - -/* -remove_train_cols function removes a given local confiugration from the training dataset - -[Input] -1. mlff_str: MLFF_Obj structure -2. col_ID: ID of the local confiugration in training dataset -[Output] -1. mlff_str: MLFF_Obj structure -*/ -void remove_train_cols(MLFF_Obj *mlff_str, int col_ID){ - int i, j; - for (i = col_ID; i < mlff_str->n_cols-1; i++){ - for (j = 0; j < mlff_str->n_rows; j++){ - mlff_str->K_train[j][i] = mlff_str->K_train[j][i+1]; - } - } - - for (j = 0; j < mlff_str->n_rows; j++){ - mlff_str->K_train[j][mlff_str->n_cols-1] = 0.0; - } - - for (i =col_ID; i < mlff_str->n_cols-1; i++){ - for (j=0; j < mlff_str->size_X3; j++){ - mlff_str->X3_traindataset[i][j] = mlff_str->X3_traindataset[i+1][j]; - } - for (j=0; j < mlff_str->size_X2; j++){ - mlff_str->X2_traindataset[i][j] = mlff_str->X2_traindataset[i+1][j]; - } - } - - for (j=0; j < mlff_str->size_X3; j++){ - mlff_str->X3_traindataset[mlff_str->n_cols-1][j] = 0.0; - } - for (j=0; j < mlff_str->size_X2; j++){ - mlff_str->X2_traindataset[mlff_str->n_cols-1][j] = 0.0; - } - - int atom_typ = mlff_str->natm_typ_train[col_ID]; - mlff_str->natm_train_elemwise[atom_typ] = mlff_str->natm_train_elemwise[atom_typ] -1; - mlff_str->natm_train_total = mlff_str->natm_train_total -1; - - for (i =col_ID; i < mlff_str->n_cols-1; i++){ - mlff_str->natm_typ_train[i] = mlff_str->natm_typ_train[i+1]; - } - - mlff_str->n_cols = mlff_str->n_cols - 1; -} - - -void get_N_r_hnl(SPARC_OBJ *pSPARC){ - int i, j, info; - FILE *fp; - char line[512]; - char a1[512], a2[512], a3[512], a4[512]; - int count1=0, count2=0; - - // fp = fopen("hnl.txt","r"); - fp = fopen(pSPARC->hnl_file_name,"r"); - - fgets(line, sizeof (line), fp); - sscanf(line, "%s%s%s%s", a1, a2, a3, a4); - int N_r = atoi(a4); - pSPARC->N_rgrid_MLFF = N_r; - fclose(fp); -} diff --git a/src/mlff/linearsys.h b/src/mlff/linearsys.h deleted file mode 100644 index 14e5bc7b..00000000 --- a/src/mlff/linearsys.h +++ /dev/null @@ -1,59 +0,0 @@ -#ifndef LINEARSYS_H -#define LINEARSYS_H - -#include "mlff_types.h" -#include "isddft.h" - -double mlff_kernel(int kernel_typ, double **X2_str, double **X3_str, double **X2_tr, double **X3_tr, int atom, int atom_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); - -double der_mlff_kernel(int kernel_typ, double ***dX2_str, double ***dX3_str, double **X2_str, double **X3_str, double **X2_tr, double **X3_tr, int atom, int atom_tr, int neighbor, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); - -double soap_kernel_polynomial(double *X2_str, double *X3_str, double *X2_tr, double *X3_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); - -double soap_kernel_Gaussian(double *X2_str, double *X3_str, double *X2_tr, double *X3_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); - -double soap_kernel_Laplacian(double *X2_str, double *X3_str, double *X2_tr, double *X3_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); - -double der_soap_kernel_polynomial(double *dX2_str, double *dX3_str, double *X2_str, double *X3_str, double *X2_tr, double *X3_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); - -double der_soap_kernel_Gaussian(double *dX2_str, double *dX3_str, double *X2_str, double *X3_str, double *X2_tr, double *X3_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); - -double der_soap_kernel_Laplacian(double *dX2_str, double *dX3_str, double *X2_str, double *X3_str, double *X2_tr, double *X3_tr, - double beta_2, double beta_3, double xi_3, int size_X2, int size_X3); - -double GMP_kernel_polynomial(double *X2_str, double *X2_tr, - double xi_3, int size_X2); - -double GMP_kernel_Gaussian(double *X2_str, double *X2_tr, - double xi_3, int size_X2); - -double der_GMP_kernel_polynomial(double *dX2_str, double *X2_str, double *X2_tr, - double xi_3, int size_X2); - -double der_GMP_kernel_Gaussian(double *dX2_str, double *X2_str, double *X2_tr, - double xi_3, int size_X2); - - -void copy_descriptors(DescriptorObj *desc_str_MLFF, DescriptorObj *desc_str); - -void add_firstMD(DescriptorObj *desc_str, NeighList *nlist, MLFF_Obj *mlff_str, double E, double* F, double* stress_sparc); - -void add_newstr_rows(DescriptorObj *desc_str, NeighList *nlist, MLFF_Obj *mlff_str, double E, double *F, double* stress_sparc); - -void calculate_Kpredict(DescriptorObj *desc_str, NeighList *nlist, MLFF_Obj *mlff_str, double **K_predict); - -void add_newtrain_cols(double *X2, double *X3, int elem_typ, MLFF_Obj *mlff_str); - -void remove_str_rows(MLFF_Obj *mlff_str, int str_ID); - -void remove_train_cols(MLFF_Obj *mlff_str, int col_ID); - -void get_N_r_hnl(SPARC_OBJ *pSPARC); -#endif diff --git a/src/mlff/soap_descriptor.c b/src/mlff/soap_descriptor.c index 87ae424e..2da8795c 100644 --- a/src/mlff/soap_descriptor.c +++ b/src/mlff/soap_descriptor.c @@ -17,57 +17,6 @@ #define min(a,b) ((a)<(b)?(a):(b)) #define temp_tol 1.0E-10 - - - - - - - -void read_h_nl(const int N, const int L, double *rgrid, double *h_nl, double *dh_nl, SPARC_OBJ *pSPARC){ - int i, j, info; - FILE *fp; - char line[512]; - char a1[512], a2[512], a3[512], a4[512]; - int count1=0, count2=0; - - // fp = fopen("hnl.txt","r"); - fp = fopen(pSPARC->hnl_file_name,"r"); - - fgets(line, sizeof (line), fp); - sscanf(line, "%s%s%s%s", a1, a2, a3, a4); - int N_r = atoi(a4); - pSPARC->N_rgrid_MLFF = N_r; - - for (i = 0; i < N_r; i++){ - info = fscanf(fp, "%s", a1); - rgrid[i] = atof(a1); - - } - - for (i =0; i < N*(L+1); i++){ - info = fscanf(fp, "%s %s", a1, a2); - - for (j = 0; j < N_r; j++){ - info = fscanf(fp, "%s", a1); - h_nl[count1] = atof(a1); - count1++; - } - for (j=0; j0) [Input] -1. soap_str: SOAP structure to the first MD structure -2. mlff_str: MLFF_Obj structure to be initialized -3. n_str_max: max reference structure in training dataset -4. n_train_max: max local descriptor per element type in training dataset +1. pSPARC: SPARC object [Output] -1. mlff_str: MLFF_Obj structure initialized +1. pSPARC: SPARC object with energies, forces, and stresses predicted by MLFF */ -void init_MLFF(SPARC_OBJ *pSPARC) { +void MLFF_call(SPARC_OBJ* pSPARC){ + // TODO: Account for socket restarts where elecgs_Count may be reset + int calc_count = pSPARC->elecgs_Count; + #ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) calc_count = pSPARC->SocketSCFCount - 1; + #endif + if (calc_count < 1) { + init_MLFF_structure(pSPARC); + init_MLFF_simulation(pSPARC); + } else { + // pSPARC->mlff_flag == 1 means on-the-fly MD from scratch, pSPARC->mlff_flag == 21 means only prediction using a known model, pSPARC->mlff_flag == 22 on-the-fly starting from a known model + if (pSPARC->mlff_flag == 1 || pSPARC->mlff_flag == 22) { + MLFF_train_predict(pSPARC, pSPARC->mlff_str); + } else if (pSPARC->mlff_flag == 21) { + MLFF_only_predict(pSPARC, pSPARC->mlff_str); + } + if (pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1) { + if(pSPARC->cell_typ != 0){ + for(int i = 0; i < pSPARC->n_atom; i++) + nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); + } + } + } +} - +/* +init_MLFF_structure function initializes and dynamically allocates memory to various objects in MLFF_Obj +[Input] +1. pSPARC: SPARC_OBJ structure +[Output] +1. pSPARC: SPARC_OBJ with MLFF_Obj mlff_str structure initialized +*/ + +void init_MLFF_structure(SPARC_OBJ *pSPARC) { + double t1, t2; +t1 = MPI_Wtime(); // initialized the mlff structure within SPARC pSPARC->mlff_str = (MLFF_Obj *) malloc(sizeof(MLFF_Obj)*1); MLFF_Obj* mlff_str = pSPARC->mlff_str; - int rank, nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - - // formatting the hnl file name start + // begin formatting the model folder path char *INPUT_DIR = malloc(L_PSD * sizeof(char)); char *inpt_path = pSPARC->filename; char *pch = strrchr(inpt_path,'/'); // find last occurrence of '/' @@ -70,75 +98,37 @@ void init_MLFF(SPARC_OBJ *pSPARC) { memcpy(INPUT_DIR, inpt_path, pch-inpt_path); INPUT_DIR[(int)(pch-inpt_path)] = '\0'; } - - char *merged_string = malloc(L_PSD * sizeof(char)); - snprintf(merged_string, L_PSD, "%s%s%s", INPUT_DIR, "/",pSPARC->hnl_file_name); - snprintf(pSPARC->hnl_file_name, L_STRING, "%s", merged_string); - free(merged_string); merged_string = malloc(L_PSD * sizeof(char)); snprintf(merged_string, L_PSD, "%s%s%s", INPUT_DIR, "/",pSPARC->mlff_data_folder); snprintf(pSPARC->mlff_data_folder, L_STRING, "%s", merged_string); free(merged_string); - // formatting the hnl file name end - - + // end formatting the model folder path // initialized the MLFF output print names int nelem = pSPARC->Ntypes; int natom = pSPARC->n_atom; - - mlff_str->xi_3 = pSPARC->xi_3_SOAP; mlff_str->Nmax = pSPARC->N_max_SOAP; mlff_str->Lmax = pSPARC->L_max_SOAP; mlff_str->rcut = pSPARC->rcut_SOAP; mlff_str->sigma_atom = pSPARC->sigma_atom_SOAP; - - if (pSPARC->descriptor_typ_MLFF < 2) { compute_hnl_soap(pSPARC, mlff_str); - // read the hnl files - // if (rank==0){ - // get_N_r_hnl(pSPARC); - - // mlff_str->rgrid = (double *) malloc(sizeof(double)* pSPARC->N_rgrid_MLFF); - // mlff_str->h_nl = (double *) malloc(sizeof(double)* pSPARC->N_rgrid_MLFF* pSPARC->N_max_SOAP*(pSPARC->L_max_SOAP+1)); - // mlff_str->dh_nl = (double *) malloc(sizeof(double)* pSPARC->N_rgrid_MLFF* pSPARC->N_max_SOAP*(pSPARC->L_max_SOAP+1)); - // read_h_nl(pSPARC->N_max_SOAP, pSPARC->L_max_SOAP, mlff_str->rgrid, mlff_str->h_nl, mlff_str->dh_nl, pSPARC); - - // MPI_Bcast(&pSPARC->N_rgrid_MLFF, 1, MPI_INT, 0, MPI_COMM_WORLD); - - // MPI_Bcast(mlff_str->rgrid, pSPARC->N_rgrid_MLFF, MPI_DOUBLE, 0, MPI_COMM_WORLD); - // MPI_Bcast(mlff_str->h_nl, pSPARC->N_rgrid_MLFF* pSPARC->N_max_SOAP*(pSPARC->L_max_SOAP+1), MPI_DOUBLE, 0, MPI_COMM_WORLD); - // MPI_Bcast(mlff_str->dh_nl, pSPARC->N_rgrid_MLFF* pSPARC->N_max_SOAP*(pSPARC->L_max_SOAP+1), MPI_DOUBLE, 0, MPI_COMM_WORLD); - // } else { - // MPI_Bcast(&pSPARC->N_rgrid_MLFF, 1, MPI_INT, 0, MPI_COMM_WORLD); - // mlff_str->rgrid = (double *) malloc(sizeof(double)* pSPARC->N_rgrid_MLFF); - // mlff_str->h_nl = (double *) malloc(sizeof(double)* pSPARC->N_rgrid_MLFF* pSPARC->N_max_SOAP*(pSPARC->L_max_SOAP+1)); - // mlff_str->dh_nl = (double *) malloc(sizeof(double)* pSPARC->N_rgrid_MLFF* pSPARC->N_max_SOAP*(pSPARC->L_max_SOAP+1)); - - // MPI_Bcast(mlff_str->rgrid, pSPARC->N_rgrid_MLFF, MPI_DOUBLE, 0, MPI_COMM_WORLD); - // MPI_Bcast(mlff_str->h_nl, pSPARC->N_rgrid_MLFF* pSPARC->N_max_SOAP*(pSPARC->L_max_SOAP+1), MPI_DOUBLE, 0, MPI_COMM_WORLD); - // MPI_Bcast(mlff_str->dh_nl, pSPARC->N_rgrid_MLFF* pSPARC->N_max_SOAP*(pSPARC->L_max_SOAP+1), MPI_DOUBLE, 0, MPI_COMM_WORLD); - // } - // size of the descriptor - int size_X3 = 0; if (pSPARC->descriptor_typ_MLFF==0){ size_X3 = ((nelem * pSPARC->N_max_SOAP+1)*(nelem * pSPARC->N_max_SOAP))/2 * (pSPARC->L_max_SOAP+1); } - mlff_str->size_X3 = size_X3; // rgrid of the hnl for spline interpolation mlff_str->N_rgrid = pSPARC->N_rgrid_MLFF; - mlff_str->X3_traindataset = (double **) malloc(sizeof(double*)*pSPARC->n_train_max_mlff*nelem); for (int i = 0; i < pSPARC->n_train_max_mlff*nelem; i++){ mlff_str->X3_traindataset[i] = (double *) malloc(sizeof(double)*size_X3); } - } // put GMP in else here - - + } else { + printf("GMP Descriptor not implemented yet\n"); + exit(1); + } mlff_str->descriptor_strdataset = (DescriptorObj *) malloc(sizeof(DescriptorObj)*pSPARC->n_str_max_mlff); // Number of stress components depedning on the BC int stress_len = (1-pSPARC->BCx) + (1-pSPARC->BCy) + (1-pSPARC->BCz) + ( (1-pSPARC->BCx-pSPARC->BCy) > 0) + ( (1-pSPARC->BCx-pSPARC->BCz) > 0) + ( (1-pSPARC->BCy-pSPARC->BCz) > 0); @@ -146,18 +136,12 @@ void init_MLFF(SPARC_OBJ *pSPARC) { if (pSPARC->mlff_pressure_train_flag){ stress_len = 1; } - - - - #ifdef DEBUG if(!rank) printf("No. of physical stress components = %d\n", stress_len); #endif - mlff_str->stress_len = stress_len; mlff_str->mlff_flag = pSPARC->mlff_flag; - int el_idx_global[natom]; // element ID of all atoms int count = 0; @@ -167,17 +151,13 @@ void init_MLFF(SPARC_OBJ *pSPARC) { count++; } } - - mlff_str->Znucl = (int *)malloc(pSPARC->Ntypes*sizeof(int)); for (int i=0; i Ntypes; i++){ mlff_str->Znucl[i] = pSPARC->Znucl[i]; } - // Atoms distributed to processors for MLFF parallelization int natom_per_procs_quot = natom/nprocs; int natom_per_procs_rem = natom%nprocs; - int natom_procs[nprocs]; for (int i=0; i < nprocs; i++){ if (i < natom_per_procs_rem){ @@ -186,24 +166,19 @@ void init_MLFF(SPARC_OBJ *pSPARC) { natom_procs[i] = natom_per_procs_quot; } } - mlff_str->natom_domain = natom_procs[rank]; mlff_str->atom_idx_domain = (int *)malloc(mlff_str->natom_domain*sizeof(int)); mlff_str->el_idx_domain = (int *)malloc(mlff_str->natom_domain*sizeof(int)); mlff_str->print_mlff_flag = pSPARC->print_mlff_flag; - char str1[100] = "MLFF_data_reference_atoms.txt"; snprintf(mlff_str->ref_atom_name, L_STRING, "%s%s%s", pSPARC->mlff_data_folder, "/",str1); // strcpy(mlff_str->ref_atom_name, str1); - char str2[100] = "MLFF_data_reference_structures.txt"; snprintf(mlff_str->ref_str_name, L_STRING, "%s%s%s", pSPARC->mlff_data_folder, "/",str2); // strcpy(mlff_str->ref_str_name, str2); - char str3[100] = "MLFF_RESTART.txt"; snprintf(mlff_str->restart_name, L_STRING, "%s%s%s", pSPARC->mlff_data_folder, "/",str3); // strcpy(mlff_str->restart_name, str3); - int strt = 0; if (rank==0){ strt = 0; @@ -217,14 +192,11 @@ void init_MLFF(SPARC_OBJ *pSPARC) { mlff_str->atom_idx_domain[i] = strt +i; mlff_str->el_idx_domain[i] = el_idx_global[mlff_str->atom_idx_domain[i]]; } - mlff_str->descriptor_typ = pSPARC->descriptor_typ_MLFF; - // This stores the scaling factors and std dev of stress components mlff_str->relative_scale_stress = (double*) calloc(stress_len, sizeof(double)); mlff_str->stress_scale = (double*) calloc(stress_len, sizeof(double)); mlff_str->std_stress = (double*) calloc(stress_len, sizeof(double)); - mlff_str->condK_min = pSPARC->condK_min; mlff_str->if_sparsify_before_train = pSPARC->if_sparsify_before_train; mlff_str->n_str_max = pSPARC->n_str_max_mlff; @@ -235,44 +207,33 @@ void init_MLFF(SPARC_OBJ *pSPARC) { mlff_str->natm_train_total = 0; mlff_str->nelem = nelem; mlff_str->natom = natom; - - mlff_str->F_tol = pSPARC->F_tol_SOAP; mlff_str->relative_scale_F = pSPARC->F_rel_scale; int BC[] = {pSPARC->BCx,pSPARC->BCy,pSPARC->BCz}; int index[] = {0,1,2,3,4,5}; - reshape_stress(pSPARC->cell_typ,BC,index); for(int i = 0; i < stress_len; i++){ mlff_str->relative_scale_stress[i] = pSPARC->stress_rel_scale[index[i]]; } - mlff_str->sigma_w = 100.0; mlff_str->sigma_v = 0.1; mlff_str->kernel_typ = pSPARC->kernel_typ_MLFF; mlff_str->error_train_E = 0.0; mlff_str->error_train_F = 0.0; mlff_str->error_train_stress = (double*) calloc(stress_len, sizeof(double)); - mlff_str->E_store_counter = 0; mlff_str->F_store_counter = 0; - mlff_str->F_store = (double *) calloc(pSPARC->n_str_max_mlff*natom*3, sizeof(double)); // hardcoded here (use realloc in future) mlff_str->E_store = (double *) calloc(pSPARC->n_str_max_mlff, sizeof(double)); // hardcoded here mlff_str->stress_store = (double **) calloc(stress_len, sizeof(double*)); - mlff_str->internal_energy_DFT = (double *) calloc(pSPARC->n_str_max_mlff, sizeof(double)); mlff_str->free_energy_DFT = (double *) calloc(pSPARC->n_str_max_mlff, sizeof(double)); mlff_str->internal_energy_DFT_count = 0; mlff_str->mlff_internal_energy_flag = pSPARC->mlff_internal_energy_flag; mlff_str->mlff_pressure_train_flag = pSPARC->mlff_pressure_train_flag; - - - for (int i = 0; i < stress_len; i++){ mlff_str->stress_store[i] = (double *) calloc(pSPARC->n_str_max_mlff, sizeof(double)); } - // initialized the arrays to be used in regression mlff_str->cov_train = (double *) malloc(1*sizeof(double)); mlff_str->AtA_SVD_U = (double *) malloc(1*sizeof(double)); @@ -283,39 +244,30 @@ void init_MLFF(SPARC_OBJ *pSPARC) { for (int i=0; inatm_train_elemwise[i] = 0; } - mlff_str->natm_typ_train = (int *)malloc(sizeof(int)*nelem * pSPARC->n_train_max_mlff); - int K_size_row = pSPARC->n_str_max_mlff * (3*mlff_str->natom_domain + 1+stress_len); int K_size_column = nelem * pSPARC->n_train_max_mlff; int b_size = pSPARC->n_str_max_mlff*(3*mlff_str->natom_domain + 1+stress_len); int w_size = nelem * pSPARC->n_train_max_mlff; - mlff_str->K_train = (double **) malloc(sizeof(double*)*K_size_row); for (int i =0; i < K_size_row; i++){ mlff_str->K_train[i] = (double *) malloc(sizeof(double)*K_size_column); } - mlff_str->b_no_norm = (double *) malloc(sizeof(double)*b_size); mlff_str->weights = (double *) malloc(sizeof(double)*w_size); - for (int i = 0; i < K_size_row; i++){ for (int j=0; j < K_size_column; j++){ mlff_str->K_train[i][j] = 0.0; } } - for (int i = 0; i < b_size; i++){ mlff_str->b_no_norm[i] = 0.0; } - for (int i = 0; i < w_size; i++){ mlff_str->weights[i] = 0.0; } - mlff_str->n_train_max = pSPARC->n_train_max_mlff; mlff_str->n_str_max = pSPARC->n_str_max_mlff; - char fname_mlff_print[100] = "mlff.log"; char fname_mlff_print_loc[L_STRING]; if (rank==0 && mlff_str->print_mlff_flag==1){ @@ -324,34 +276,312 @@ void init_MLFF(SPARC_OBJ *pSPARC) { mlff_str->fp_mlff = fopen(fname_mlff_print_loc,"w"); } free(INPUT_DIR); + if (pSPARC->descriptor_typ_MLFF>1 && pSPARC->elec_T > 2000){ + if (rank==0){ + printf("WARNING: GMP for high temepratures calculations has not been tested yet! Please use SOAP.\n"); + exit(1); + } + } +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "MLFF object intialized. Time taken: %.3f s\n", t2-t1); + } +} + +/** + * @brief Initializes the simulations that will be using the MLFF structure. + */ + +void init_MLFF_simulation(SPARC_OBJ *pSPARC) { + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MLFF_Obj* mlff_str = pSPARC->mlff_str; + FILE *fp_mlff; + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fp_mlff = mlff_str->fp_mlff; + fprintf(fp_mlff,"--------------------------------------------------------------------------------------------\n"); + fprintf(fp_mlff, "Simulation step: 0\n"); + } + switch (pSPARC->mlff_flag) { + case 1: + init_new_MLFF(pSPARC, mlff_str); + break; + case 22: + init_existing_trainable_MLFF(pSPARC, mlff_str); + break; + case 21: + init_static_MLFF(pSPARC, mlff_str); + break; + default: + if (rank == 0) { + printf("Error: MLFF flag not recognized\n"); + exit(1); + } + } +} - +/** + * @brief Initializes a new MLFF simulation from scratch. + */ - +void init_new_MLFF(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str) { + int rank; + double t1, t2; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + intialize_print_MLFF(mlff_str, pSPARC); + pSPARC->last_train_iter = -1; + // Puts coords into cart if MD or Relax flag = 1 +t1 = MPI_Wtime(); + Calculate_electronicGroundState(pSPARC); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "DFT call done for the first step. Time taken: %.3f s\n", t2-t1); + } + MPI_Bcast(pSPARC->forces, 3*pSPARC->n_atom, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(pSPARC->stress, 6, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&pSPARC->pres, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count] = pSPARC->Etot - pSPARC->Entropy; + mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count] = pSPARC->Etot; + if (pSPARC->print_mlff_flag == 1 && rank ==0 && pSPARC->mlff_internal_energy_flag){ + fprintf(mlff_str->fp_mlff, "Internal energy DFT data added! Free energy: %.6E Ha, Internal energy: %.6E Ha, Entropy: %.6E Ha\n", + mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count], mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count], pSPARC->Entropy); + } + mlff_str->internal_energy_DFT_count += 1; + if(pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1){ + if(pSPARC->cell_typ != 0){ + for(int i = 0; i < pSPARC->n_atom; i++) + Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); + } + } +t1 = MPI_Wtime(); + sparc_mlff_interface_initialDFT(pSPARC, mlff_str); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "Covariance matrices formed from first DFT data. Time taken: %.3f s\n", t2-t1); + } + if(pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1){ + if(pSPARC->cell_typ != 0){ + for(int i = 0; i < pSPARC->n_atom; i++) + nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); + } + } } -/* -free_MLFF function frees the memory allocated to various objects in MLFF_Obj +/** + * @brief Initializes a new MLFF simulation from a known model that can be updated. + */ -[Input] -1. mlff_str: MLFF_Obj structure -[Output] -1. mlff_str: MLFF_Obj structure -*/ +void init_existing_trainable_MLFF(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str) { + int rank; + double t1, t2; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + if (!rank){ + // TODO: Probably need to edit or debug + read_MLFF_files(mlff_str, pSPARC); + } +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "Existing MLFF ref-atom file read. Time taken: %.3f s\n", t2-t1); + } +t1 = MPI_Wtime(); + MPI_Bcast(&mlff_str->n_cols, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->n_str, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->weights, mlff_str->n_cols, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->mu_E, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->std_E, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->std_F, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->std_stress, mlff_str->stress_len, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->E_scale, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->F_scale, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->stress_scale, mlff_str->stress_len, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->n_str, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->natm_train_elemwise, pSPARC->Ntypes, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->n_rows, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->natm_train_total, 1, MPI_INT, 0, MPI_COMM_WORLD); + for (int i=0; i < mlff_str->n_cols; i++){ + MPI_Bcast(mlff_str->X3_traindataset[i], mlff_str->size_X3, MPI_DOUBLE, 0, MPI_COMM_WORLD); + } + MPI_Bcast(&mlff_str->relative_scale_F, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->natm_typ_train, mlff_str->n_cols, MPI_INT, 0, MPI_COMM_WORLD); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "MLFF ref-atom read data broadcasted. Time taken: %.3f s\n", t2-t1); + } + char *fname_str; + fname_str = (char *) malloc(sizeof(char)*512); + char str1[512] = "MLFF_data_reference_structures.txt"; + strcpy(fname_str, str1); + double **cell_data; + double **LatUVec_data; + double **apos_data; + double *Etot_data; + double **F_data; + double **stress_data; + int *natom_data; + int **natom_elem_data; + natom_data = (int *) malloc(sizeof(int)*mlff_str->n_str); + Etot_data = (double *) malloc(sizeof(double)*mlff_str->n_str); + cell_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); + LatUVec_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); + stress_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); + natom_elem_data = (int **) malloc(sizeof(int*)*mlff_str->n_str); + for (int i = 0; i < mlff_str->n_str; i++){ + cell_data[i] = (double *) malloc(sizeof(double)*3); + LatUVec_data[i] = (double *) malloc(sizeof(double)*9); + stress_data[i] = (double *) malloc(sizeof(double)*6); + natom_elem_data[i] = (int *) malloc(sizeof(int)*pSPARC->Ntypes); + } +t1 = MPI_Wtime(); + if (rank==0){ + apos_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); + F_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); + // TODO: Edit and debug + read_structures_MLFF_data(mlff_str->ref_str_name, mlff_str->n_str, pSPARC->Ntypes, cell_data, LatUVec_data, apos_data, Etot_data, F_data, stress_data, natom_data, natom_elem_data); + if(pSPARC->cell_typ != 0){ + for (int istr = 0; istr < mlff_str->n_str; istr++){ + // TODO: handling training on a mix of orthogonal and non-orthogonal structures + for (int i = 0; i < natom_data[istr]; i++) { + Cart2nonCart_coord(pSPARC, &apos_data[istr][3*i], &apos_data[istr][3*i+1], &apos_data[istr][3*i+2]); + } + } + } + } +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "Read structures data from the exisiting MLFF model. Time taken: %.3f s\n", t2-t1); + } +t1 = MPI_Wtime(); + pretrain_MLFF_model(mlff_str, pSPARC, cell_data, LatUVec_data, apos_data, Etot_data, F_data, stress_data, natom_data, natom_elem_data); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "Pretraining from the exisiting data done. Time taken: %.3f s\n", t2-t1); + } +t1 = MPI_Wtime(); + Calculate_electronicGroundState(pSPARC); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "DFT call made after pretraining. Time taken: %.3f s\n", t2-t1); + } + MPI_Bcast(pSPARC->forces, 3*pSPARC->n_atom, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(pSPARC->stress, 6, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&pSPARC->pres, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + if (pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1) { + if(pSPARC->cell_typ != 0){ + for(int i = 0; i < pSPARC->n_atom; i++) + Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); + } + } + init_dyarray(&mlff_str->atom_idx_addtrain); + for (int i = 0; i < pSPARC->n_atom; i++){ + append_dyarray(&(mlff_str->atom_idx_addtrain),i); + } +t1 = MPI_Wtime(); + sparc_mlff_interface_addDFT(pSPARC, mlff_str); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "Covariance matrices formed after the DFT call. Time taken: %.3f s\n", t2-t1); + } +t1 = MPI_Wtime(); + mlff_train_Bayesian(mlff_str); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "MLFF model trained with Bayesian regression. Time taken: %.3f s\n", t2-t1); + } + pSPARC->last_train_iter = 0; + delete_dyarray(&mlff_str->atom_idx_addtrain); + for (int i = 0; i < mlff_str->n_str-1; i++){ + if(rank == 0){ + free(apos_data[i]); // Check this for memory leak + free(F_data[i]); // Check this for memory leak + } + free(cell_data[i]); + free(LatUVec_data[i]); + free(stress_data[i]); + free(natom_elem_data[i]); + } + if (rank == 0) { + free(apos_data); // Check this for memory leak + free(F_data); // Check this for memory leak + } + free(cell_data); + free(LatUVec_data); + free(stress_data); + free(natom_elem_data); + free(Etot_data); + free(natom_data); + if (pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1) { + if(pSPARC->cell_typ != 0){ + for(int i = 0; i < pSPARC->n_atom; i++) + nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); + } + } +} + +/** + * @brief Initializes a new MLFF simulation from a known model that remains static. + */ + +void init_static_MLFF(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str) { + int rank; + double t1, t2; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + if (!rank){ + // TODO: Probably need to edit or debug + read_MLFF_files(mlff_str, pSPARC); + } +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "Existing MLFF ref-atom file read. Time taken: %.3f s\n", t2-t1); + } +t1 = MPI_Wtime(); + MPI_Bcast(&mlff_str->n_cols, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->n_str, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->weights, mlff_str->n_cols, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->mu_E, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->std_E, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->std_F, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->std_stress, mlff_str->stress_len, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->E_scale, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->F_scale, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->stress_scale, mlff_str->stress_len, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->n_str, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->natm_train_elemwise, pSPARC->Ntypes, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->n_rows, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&mlff_str->natm_train_total, 1, MPI_INT, 0, MPI_COMM_WORLD); + for (int i=0; i < mlff_str->n_cols; i++){ + MPI_Bcast(mlff_str->X3_traindataset[i], mlff_str->size_X3, MPI_DOUBLE, 0, MPI_COMM_WORLD); + } + MPI_Bcast(&mlff_str->relative_scale_F, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast(mlff_str->natm_typ_train, mlff_str->n_cols, MPI_INT, 0, MPI_COMM_WORLD); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "MLFF ref-atom read data broadcasted. Time taken: %.3f s\n", t2-t1); + } + MLFF_only_predict(pSPARC, mlff_str); +t2 = MPI_Wtime(); + if (pSPARC->print_mlff_flag == 1 && rank ==0){ + fprintf(mlff_str->fp_mlff, "MLFF model prediction for E, F, stress done. Time taken: %.3f s\n", t2-t1); + } + if(pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1){ + if(pSPARC->cell_typ != 0){ + for(int i = 0; i < pSPARC->n_atom; i++) + nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); + } + } +} + +/** + * @brief Frees memory allocated for MLFF structure. + */ void free_MLFF(MLFF_Obj *mlff_str){ - int rank, nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - - if (rank==0 && mlff_str->print_mlff_flag==1){ fclose(mlff_str->fp_mlff); } - int natom = mlff_str->natom, K_size_row = mlff_str->n_str_max*(3*mlff_str->natom_domain + 1 + mlff_str->stress_len), nelem = mlff_str->nelem; - free(mlff_str->Znucl); free(mlff_str->atom_idx_domain); free(mlff_str->el_idx_domain); @@ -360,7 +590,6 @@ void free_MLFF(MLFF_Obj *mlff_str){ for(int i = 0; i < mlff_str->stress_len; i++){ free(mlff_str->stress_store[i]); } - free(mlff_str->stress_store); free(mlff_str->cov_train); free(mlff_str->AtA_SVD_U); @@ -372,7 +601,6 @@ void free_MLFF(MLFF_Obj *mlff_str){ for (int i = 0; i < (mlff_str->n_str_max * (3*mlff_str->natom_domain + 1+mlff_str->stress_len)); i++){ free(mlff_str->K_train[i]); } - free(mlff_str->K_train); free(mlff_str->b_no_norm); free(mlff_str->weights); @@ -381,8 +609,7 @@ void free_MLFF(MLFF_Obj *mlff_str){ delete_descriptor(mlff_str->descriptor_strdataset+i); } free(mlff_str->descriptor_strdataset); - } - + } if (mlff_str->descriptor_typ < 2) { for (int i=0; i < mlff_str->n_train_max*nelem; i++){ free(mlff_str->X3_traindataset[i]); @@ -390,305 +617,35 @@ void free_MLFF(MLFF_Obj *mlff_str){ free(mlff_str->X3_traindataset); free(mlff_str->rgrid); free(mlff_str->h_nl); - free(mlff_str->dh_nl); - } else { - for (int i=0; i < mlff_str->size_X3; i++) { - free(mlff_str->params_i[i]); - } - free(mlff_str->params_i); - for (int i = 0; i < mlff_str->size_X3; i++){ - free(mlff_str->params_d[i]); - } - free(mlff_str->params_d); - for (int i=0; i < mlff_str->nelem; i++){ - free(mlff_str->atom_gaussians_p[i]); - } - free(mlff_str->atom_gaussians_p); - free(mlff_str->ngaussians); - free(mlff_str->element_index_to_order_p); - free(mlff_str->atom_indices_p); - free(mlff_str->atom_type_to_indices_p); - free(mlff_str->sigmas); - } - - free(mlff_str->relative_scale_stress); - free(mlff_str->stress_scale); - free(mlff_str->std_stress); - free(mlff_str->error_train_stress); - -} - - -/* -MLFF_main function initializes the MLFF calculations, also sets up MLFF restarts - -[Input] -1. pSPARC: SPARC structure -[Output] -1. pSPARC: SPARC structure - -Important variables: -pSPARC->mlff_flag==1 : Start on-the-fly from scratch (don't use any existing model) -pSPARC->mlff_flag==21 : Only use the trained MLFF model and do not train -pSPARC->mlff_flag==22 : Continue training the existing MLFF model -*/ - -void MLFF_main(SPARC_OBJ *pSPARC){ - int rank; - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - double t1, t2, t3, t4; - -t1 = MPI_Wtime(); - init_MLFF(pSPARC); -t2 = MPI_Wtime(); - - MLFF_Obj* mlff_str = pSPARC->mlff_str; - FILE *fp_mlff; - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - // fp_mlff = fopen("mlff.log","a"); - fp_mlff = mlff_str->fp_mlff; - fprintf(fp_mlff,"--------------------------------------------------------------------------------------------\n"); - fprintf(fp_mlff, "MD step: 0\n"); - } - - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "MLFF object intialized. Time taken: %.3f s\n", t2-t1); - } - - if (pSPARC->mlff_flag==1) { - intialize_print_MLFF(mlff_str, pSPARC); -t1 = MPI_Wtime(); - Calculate_electronicGroundState(pSPARC); -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "DFT call done for the first MD. Time taken: %.3f s\n", t2-t1); - } - - - MPI_Bcast(pSPARC->forces, 3*pSPARC->n_atom, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(pSPARC->stress, 6, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&pSPARC->pres, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - - mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count] = pSPARC->Etot - pSPARC->Entropy; - - - mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count] = pSPARC->Etot; - - if (pSPARC->print_mlff_flag == 1 && rank ==0 && pSPARC->mlff_internal_energy_flag){ - fprintf(fp_mlff, "Internal energy DFT data added! Free energy: %.6E Ha, Internal energy: %.6E Ha, Entropy: %.6E Ha\n", - mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count], mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count], pSPARC->Entropy); - } - - mlff_str->internal_energy_DFT_count += 1; - - Initialize_MD(pSPARC); - if(pSPARC->cell_typ != 0){ - coordinatetransform_map(pSPARC, pSPARC->n_atom, pSPARC->atom_pos); - } - -t1 = MPI_Wtime(); - sparc_mlff_interface_firstMD(pSPARC, mlff_str); -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "Covariance matrices formed from first MD DFT data. Time taken: %.3f s\n", t2-t1); - } - - if(pSPARC->cell_typ != 0){ - for(int i = 0; i < pSPARC->n_atom; i++) - nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); - } - - } else { -t1 = MPI_Wtime(); - if (!rank){ - // TODO: Probably need to edit or debug - read_MLFF_files(mlff_str, pSPARC); - } -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "Existing MLFF ref-atom file read. Time taken: %.3f s\n", t2-t1); - } -t1 = MPI_Wtime(); - MPI_Bcast(&mlff_str->n_cols, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->n_str, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(mlff_str->weights, mlff_str->n_cols, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->mu_E, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->std_E, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->std_F, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(mlff_str->std_stress, mlff_str->stress_len, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->E_scale, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->F_scale, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(mlff_str->stress_scale, mlff_str->stress_len, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->n_str, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(mlff_str->natm_train_elemwise, pSPARC->Ntypes, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->n_rows, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(&mlff_str->natm_train_total, 1, MPI_INT, 0, MPI_COMM_WORLD); - for (int i=0; i < mlff_str->n_cols; i++){ - MPI_Bcast(mlff_str->X3_traindataset[i], mlff_str->size_X3, MPI_DOUBLE, 0, MPI_COMM_WORLD); - } - MPI_Bcast(&mlff_str->relative_scale_F, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(mlff_str->natm_typ_train, mlff_str->n_cols, MPI_INT, 0, MPI_COMM_WORLD); -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "MLFF ref-atom read data broadcasted. Time taken: %.3f s\n", t2-t1); - } - if (pSPARC->mlff_flag==21){ -t1 = MPI_Wtime(); - MLFF_call_from_MD_only_predict(pSPARC, mlff_str); -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "MLFF model prediction for E, F, stress done. Time taken: %.3f s\n", t2-t1); - } - if(pSPARC->cell_typ != 0){ - for(int i = 0; i < pSPARC->n_atom; i++) - nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); - } - - Initialize_MD(pSPARC); - } else if (pSPARC->mlff_flag ==22){ - char *fname_str; - fname_str = (char *) malloc(sizeof(char)*512); - char str1[512] = "MLFF_data_reference_structures.txt"; - strcpy(fname_str, str1); - - double **cell_data; - double **LatUVec_data; - double **apos_data; - double *Etot_data; - double **F_data; - double **stress_data; - int *natom_data; - int **natom_elem_data; - - natom_data = (int *) malloc(sizeof(int)*mlff_str->n_str); - Etot_data = (double *) malloc(sizeof(double)*mlff_str->n_str); - cell_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); - LatUVec_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); - stress_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); - natom_elem_data = (int **) malloc(sizeof(int*)*mlff_str->n_str); - for (int i = 0; i < mlff_str->n_str; i++){ - cell_data[i] = (double *) malloc(sizeof(double)*3); - LatUVec_data[i] = (double *) malloc(sizeof(double)*9); - stress_data[i] = (double *) malloc(sizeof(double)*6); - natom_elem_data[i] = (int *) malloc(sizeof(int)*pSPARC->Ntypes); - } - -t1 = MPI_Wtime(); - if (rank==0){ - apos_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); - F_data = (double **) malloc(sizeof(double*)*mlff_str->n_str); - // TODO: Edit and debug - read_structures_MLFF_data(mlff_str->ref_str_name, mlff_str->n_str, pSPARC->Ntypes, cell_data, LatUVec_data, apos_data, Etot_data, F_data, stress_data, natom_data, natom_elem_data); - if(pSPARC->cell_typ != 0){ - for (int istr = 0; istr < mlff_str->n_str; istr++){ - coordinatetransform_map(pSPARC, natom_data[istr], apos_data[istr]); - } - } - } - -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "Read structures data from the exisiting MLFF model. Time taken: %.3f s\n", t2-t1); - } - -t1 = MPI_Wtime(); - pretrain_MLFF_model(mlff_str, pSPARC, cell_data, LatUVec_data, apos_data, Etot_data, F_data, stress_data, natom_data, natom_elem_data); -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "Pretraining from the exisiting data done. Time taken: %.3f s\n", t2-t1); - } - -t1 = MPI_Wtime(); - Calculate_electronicGroundState(pSPARC); -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "DFT call made after pretraining. Time taken: %.3f s\n", t2-t1); - } - MPI_Bcast(pSPARC->forces, 3*pSPARC->n_atom, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(pSPARC->stress, 6, MPI_DOUBLE, 0, MPI_COMM_WORLD); - MPI_Bcast(&pSPARC->pres, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - - Initialize_MD(pSPARC); - if(pSPARC->cell_typ != 0){ - coordinatetransform_map(pSPARC, pSPARC->n_atom, pSPARC->atom_pos); - } - init_dyarray(&mlff_str->atom_idx_addtrain); - for (int i = 0; i < pSPARC->n_atom; i++){ - append_dyarray(&(mlff_str->atom_idx_addtrain),i); - } - -t1 = MPI_Wtime(); - sparc_mlff_interface_addMD(pSPARC, mlff_str); -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "Covariance matrices formed after the DFT call. Time taken: %.3f s\n", t2-t1); - } - -t1 = MPI_Wtime(); - mlff_train_Bayesian(mlff_str); -t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "MLFF model trained with Bayesian regression. Time taken: %.3f s\n", t2-t1); - } - - pSPARC->last_train_MD_iter = 0; - delete_dyarray(&mlff_str->atom_idx_addtrain); - - // MPI_Barrier(MPI_COMM_WORLD); - // exit(3); - - for (int i = 0; i < mlff_str->n_str-1; i++){ - if(rank == 0){ - free(apos_data[i]); // Check this for memory leak - free(F_data[i]); // Check this for memory leak - } - free(cell_data[i]); - free(LatUVec_data[i]); - free(stress_data[i]); - free(natom_elem_data[i]); - - } - - if (rank == 0) { - free(apos_data); // Check this for memory leak - free(F_data); // Check this for memory leak - } - - free(cell_data); - free(LatUVec_data); - free(stress_data); - free(natom_elem_data); - - free(Etot_data); - free(natom_data); - - if(pSPARC->cell_typ != 0){ - for(int i = 0; i < pSPARC->n_atom; i++) - nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); - } - + free(mlff_str->dh_nl); + } else { + for (int i=0; i < mlff_str->size_X3; i++) { + free(mlff_str->params_i[i]); + } + free(mlff_str->params_i); + for (int i = 0; i < mlff_str->size_X3; i++){ + free(mlff_str->params_d[i]); } - } + free(mlff_str->params_d); + for (int i=0; i < mlff_str->nelem; i++){ + free(mlff_str->atom_gaussians_p[i]); + } + free(mlff_str->atom_gaussians_p); + free(mlff_str->ngaussians); + free(mlff_str->element_index_to_order_p); + free(mlff_str->atom_indices_p); + free(mlff_str->atom_type_to_indices_p); + free(mlff_str->sigmas); + } + free(mlff_str->relative_scale_stress); + free(mlff_str->stress_scale); + free(mlff_str->std_stress); + free(mlff_str->error_train_stress); } -/* -pretrain_MLFF_model function trains the MLFF model from the available SPARC-DFT data of energy, force and stress stored from previous on-the-fly run - -[Input] -1. mlff_str: MLFF structure -2. pSPARC: SPARC structure -3. cell_data: lattice constant of all structures -4. apos_data: atom positions of all structures -5. Etot_data: Total energy of all structures -6. F_data: Forces in all structures -7. stress_data: stresss on al structures -8. natom_data: number of atoms on all structures -9. natom_elem_data: Number of atoms of each eleement type of all structures -[Output] -1. mlff_str: MLFF structure - -*/ - +/** + * @brief Constructs design matrices from saved reference structures and reference atomic descriptors. + */ void pretrain_MLFF_model( MLFF_Obj *mlff_str, @@ -705,30 +662,23 @@ void pretrain_MLFF_model( int rank, nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - double t1, t2, t3, t4; - FILE *fp_mlff; if (pSPARC->print_mlff_flag == 1 && rank ==0){ fp_mlff = mlff_str->fp_mlff; } - double cell[3], LatUVec[9], *apos, Etot, *F, *stress, *stress1; int *natom_elem, *atomtyp, count, natom_domain, *atom_idx_domain, *el_idx_domain; NeighList *nlist; DescriptorObj *desc_str; - int n_str = mlff_str->n_str; mlff_str->n_str = 0; stress1 = (double *) malloc(sizeof(double)*mlff_str->stress_len); - int *BC = (int*)malloc(3*sizeof(int)); BC[0] = pSPARC->BCx; BC[1] = pSPARC->BCy; BC[2] = pSPARC->BCz; - MPI_Bcast(natom_data, n_str, MPI_INT, 0, MPI_COMM_WORLD); - mlff_str->n_rows = 0; for (int i = 0; i < n_str; i++){ apos = (double *) malloc(sizeof(double)*3*natom_data[i]); @@ -737,12 +687,10 @@ void pretrain_MLFF_model( nlist = (NeighList *) malloc(sizeof(NeighList)*1); desc_str = (DescriptorObj *) malloc(sizeof(DescriptorObj)*1); natom_elem = (int *) malloc(sizeof(int)*pSPARC->Ntypes); - if (rank==0){ cell[0] = cell_data[i][0]; cell[1] = cell_data[i][1]; cell[2] = cell_data[i][2]; - for (int j = 0; j < 9; j++){ LatUVec[j] = LatUVec_data[i][j]; } @@ -765,15 +713,12 @@ void pretrain_MLFF_model( MPI_Bcast(stress, 6, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&Etot, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(natom_elem, pSPARC->Ntypes, MPI_INT, 0, MPI_COMM_WORLD); - int index[6] = {0,1,2,3,4,5}; reshape_stress(pSPARC->cell_typ, BC, index); for(int j = 0; j < 6; j++){ stress[j] = stress[index[j]]; } - atomtyp = (int *) malloc(natom_data[i]*sizeof(int)); - count = 0; for (int ii=0; ii < pSPARC->Ntypes; ii++){ for (int j=0; j < natom_elem[ii]; j++){ @@ -781,28 +726,21 @@ void pretrain_MLFF_model( count++; } } - get_domain_decompose_mlff_natom(natom_data[i], pSPARC->Ntypes, pSPARC->nAtomv, nprocs, rank, &natom_domain); - atom_idx_domain = (int *)malloc(sizeof(int)*natom_domain); el_idx_domain = (int *)malloc(sizeof(int)*natom_domain); - get_domain_decompose_mlff_idx(natom_data[i], pSPARC->Ntypes, pSPARC->nAtomv, nprocs, rank, natom_domain, atom_idx_domain, el_idx_domain); - double *geometric_ratio = (double*) malloc(2 * sizeof(double)); geometric_ratio[0] = pSPARC->CUTOFF_y[0]/pSPARC->CUTOFF_x[0]; geometric_ratio[1] = pSPARC->CUTOFF_z[0]/pSPARC->CUTOFF_x[0]; - t1 = MPI_Wtime(); build_nlist(mlff_str->rcut, pSPARC->Ntypes, natom_data[i], apos, atomtyp, pSPARC->cell_typ, BC, cell, LatUVec, pSPARC->twist, geometric_ratio, nlist, natom_domain, atom_idx_domain, el_idx_domain); t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Neighbor list done. Time taken: %.3f s\n", t2-t1); } - t1 = MPI_Wtime(); build_descriptor(desc_str, nlist, mlff_str, apos); - t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "SOAP descriptor done. Time taken: %.3f s\n", t2-t1); @@ -810,15 +748,12 @@ t2 = MPI_Wtime(); for (int j = 0; j < mlff_str->stress_len; j++){ stress1[j] = stress[j]/au2GPa; // check this } - t1 = MPI_Wtime(); add_newstr_rows(desc_str, nlist, mlff_str, Etot/natom_data[i], F, stress1); t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Rows added for structure no. %d. Time taken: %.3f s\n", i+1, t2-t1); } - free(apos); free(F); free(stress); @@ -836,25 +771,14 @@ t2 = MPI_Wtime(); printf("Added new structure # %d\n",i+1); } #endif - } free(stress1); free(BC); } -/* -get_domain_decompose_mlff_natom function obtains the number of atoms for parallelization in MLFF - -[Input] -1. natom: number of atom in the structure -2. nelem: number of element types in the structure -3. nAtomv: number of atom of different elements in the structure -4. nprocs: total number of processors -5. rank: rank for which the decomposition to be outputted -[Output] -1. natom_domain: number of atoms to be handled by the processor rank - -*/ +/** + * @brief Gets the number of atoms on each processor + */ void get_domain_decompose_mlff_natom( int natom, @@ -866,9 +790,7 @@ void get_domain_decompose_mlff_natom( { int natom_per_procs_quot = natom/nprocs; int natom_per_procs_rem = natom%nprocs; - int el_idx_global[natom]; - // element ID of all atoms int count = 0; for (int i=0; i < nelem; i++){ @@ -877,7 +799,6 @@ void get_domain_decompose_mlff_natom( count++; } } - int natom_procs[nprocs]; for (int i=0; i < nprocs; i++){ if (i < natom_per_procs_rem){ @@ -886,24 +807,12 @@ void get_domain_decompose_mlff_natom( natom_procs[i] = natom_per_procs_quot; } } - *natom_domain = natom_procs[rank]; } -/* -get_domain_decompose_mlff_idx function obtains the idx of those atoms for parallelization in MLFF - -[Input] -1. natom: number of atom in the structure -2. nelem: number of element types in the structure -3. nAtomv: number of atom of different elements in the structure -4. nprocs: total number of processors -5. rank: rank for which the decomposition to be outputted -6. natom_domain: number of atoms to be handled by the processor rank -[Output] -1. atom_idx_domain: index of atoms to be handled by the processor rank -2. el_idx_domain: index of element types of the atoms to be handled by the processor rank -*/ +/** + * @brief Returns the atom and element indices for parallelization in MLFF + */ void get_domain_decompose_mlff_idx( int natom, @@ -917,9 +826,7 @@ void get_domain_decompose_mlff_idx( { int natom_per_procs_quot = natom/nprocs; int natom_per_procs_rem = natom%nprocs; - int el_idx_global[natom]; - // element ID of all atoms int count = 0; for (int i=0; i < nelem; i++){ @@ -928,7 +835,6 @@ void get_domain_decompose_mlff_idx( count++; } } - int natom_procs[nprocs]; for (int i=0; i < nprocs; i++){ if (i < natom_per_procs_rem){ @@ -937,7 +843,6 @@ void get_domain_decompose_mlff_idx( natom_procs[i] = natom_per_procs_quot; } } - int strt = 0; if (rank==0){ strt = 0; @@ -953,82 +858,54 @@ void get_domain_decompose_mlff_idx( } } -/* -MLFF_call function acts as the interface in the md.c file. This function is called if (mlff_idx>0) whenever the electronicgroundstate was called before - -[Input] -1. pSPARC: SPARC object -[Output] -1. pSPARC: SPARC object -*/ - -void MLFF_call(SPARC_OBJ* pSPARC){ - // pSPARC->mlff_flag == 1 means on-the-fly MD from scratch, pSPARC->mlff_flag == 21 means only prediction using a known model, pSPARC->mlff_flag == 22 on-the-fly MD starting from a known model - if (pSPARC->mlff_flag == 1 || pSPARC->mlff_flag == 22) { - MLFF_call_from_MD(pSPARC, pSPARC->mlff_str); - } else if (pSPARC->mlff_flag == 21) { - MLFF_call_from_MD_only_predict(pSPARC, pSPARC->mlff_str); - } - - if(pSPARC->cell_typ != 0){ - for(int i = 0; i < pSPARC->n_atom; i++) - nonCart2Cart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); - } - -} - -/* -MLFF_call_from_MD function performs the on-the-fly MLFF (it has the logic implemented to skip DFT steps on the basis of UQ) - -[Input] -1. pSPARC: SPARC object -2. mlff_str: MLFF object -[Output] -1. pSPARC: SPARC object -2. mlff_str: MLFF object -*/ - -void MLFF_call_from_MD(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str){ +/** + * @brief Returns predictions from MLFF or DFT depending on size of training set or UQ. Updates MLFF model as necessary. + */ +void MLFF_train_predict(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str){ double t1, t2, t3, t4; t3 = MPI_Wtime(); - int nprocs, rank; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - FILE *fp_mlff; if (pSPARC->print_mlff_flag == 1 && rank ==0){ // fp_mlff = fopen("mlff.log","a"); fp_mlff = mlff_str->fp_mlff; } - int doMLFF = pSPARC->mlff_flag; int initial_train; if (pSPARC->mlff_flag==1){ initial_train= pSPARC->begin_train_steps; } + // TODO: Allow restart from untrained MLFF model with existing structures and atoms files if (pSPARC->mlff_flag==22){ initial_train = -1; - + if (mlff_str->n_str < 1) { + printf("Error: Attempting to restart a mlff simulation from pretrained model, but no reference structures found. Exiting...\n"); + exit(1); + } } - if (pSPARC->print_mlff_flag==1 && rank==0){ fprintf(fp_mlff,"--------------------------------------------------------------------------------------------\n"); - fprintf(fp_mlff, "MD step: %d\n", pSPARC->MDCount); + int print_count = pSPARC->elecgs_Count; + #ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) print_count = pSPARC->SocketSCFCount - 1; + #endif + fprintf(fp_mlff, "Calculation: %d\n", print_count); } - int BC[] = {pSPARC->BCx, pSPARC->BCy, pSPARC->BCz}; int index[] = {0,1,2,3,4,5}; reshape_stress(pSPARC->cell_typ, BC, index); - - if (pSPARC->MDCount <= initial_train){ + if (mlff_str->n_str <= initial_train && pSPARC->last_train_iter < 0){ t1 = MPI_Wtime(); Calculate_electronicGroundState(pSPARC); - if(pSPARC->cell_typ != 0){ - coordinatetransform_map(pSPARC, pSPARC->n_atom, pSPARC->atom_pos); + if (pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1) { + if(pSPARC->cell_typ != 0){ + for(int i = 0; i < pSPARC->n_atom; i++) + Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); + } } - t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "DFT call made. Time taken: %.3f s\n", t2-t1); @@ -1036,7 +913,6 @@ t2 = MPI_Wtime(); MPI_Bcast(pSPARC->forces, 3*pSPARC->n_atom, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(pSPARC->stress, 6, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&pSPARC->pres, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count] = pSPARC->Etot-pSPARC->Entropy; mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count] = pSPARC->Etot; if (pSPARC->print_mlff_flag == 1 && rank ==0 && mlff_str->mlff_internal_energy_flag){ @@ -1044,29 +920,30 @@ t2 = MPI_Wtime(); mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count], mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count], pSPARC->Entropy); } mlff_str->internal_energy_DFT_count += 1; - init_dyarray(&mlff_str->atom_idx_addtrain); for (int i = 0; i < pSPARC->n_atom; i++){ append_dyarray(&(mlff_str->atom_idx_addtrain),i); } t1 = MPI_Wtime(); - sparc_mlff_interface_addMD(pSPARC, mlff_str); + sparc_mlff_interface_addDFT(pSPARC, mlff_str); t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "The covariance matrix updated. Time taken: %.3f s\n", t2-t1); } - if (pSPARC->MDCount == initial_train){ + if (mlff_str->n_str == initial_train){ t1 = MPI_Wtime(); mlff_train_Bayesian(mlff_str); if (pSPARC->mlff_internal_energy_flag){ train_internal_energy_model(mlff_str); } - t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Bayesian linear regression done. Time taken: %.3f s\n", t2-t1); } - pSPARC->last_train_MD_iter = pSPARC->MDCount; + pSPARC->last_train_iter = pSPARC->elecgs_Count; + #ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) pSPARC->last_train_iter = pSPARC->SocketSCFCount - 1; + #endif } delete_dyarray(&mlff_str->atom_idx_addtrain); } else { @@ -1079,13 +956,16 @@ t1 = MPI_Wtime(); sparc_mlff_interface_predict(pSPARC, mlff_str, E_predict, F_predict, stress_predict, bayesian_error); t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "MLFF prediction made for the current MD. Time taken: %.3f s\n", t2-t1); + fprintf(fp_mlff, "MLFF prediction made for the current step. Time taken: %.3f s\n", t2-t1); } double max_pred_error = fabs(largest(&bayesian_error[1+mlff_str->stress_len], 3*pSPARC->n_atom)); - if (pSPARC->MDCount == pSPARC->last_train_MD_iter + 1){ + int train_iter_compare = pSPARC->elecgs_Count; + #ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) train_iter_compare = pSPARC->SocketSCFCount - 1; + #endif + if (train_iter_compare == pSPARC->last_train_iter + 1) { pSPARC->F_tol_SOAP = pSPARC->factor_multiply_sigma_tol*max_pred_error; } - if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "max_pred_error: %.9E, F_tol_SOAP: %.9E.\n", max_pred_error, pSPARC->F_tol_SOAP); } @@ -1094,17 +974,27 @@ t2 = MPI_Wtime(); printf("max_pred_error: %10.9f, pSPARC->F_tol_SOAP: %10.9f\n",max_pred_error,pSPARC->F_tol_SOAP); } #endif - int Count_MD = pSPARC->MDCount + pSPARC->restartCount + (pSPARC->RestartFlag == 0); - if (max_pred_error >= pSPARC->F_tol_SOAP || !(Count_MD % pSPARC->MLFF_DFT_fq)){ + // TODO: Confirm this logic works + int Count_Sim = pSPARC->elecgs_Count + pSPARC->restartCount + (pSPARC->RestartFlag == 0); + #ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) Count_Sim = pSPARC->SocketSCFCount - 1 + pSPARC->restartCount + (pSPARC->RestartFlag == 0); + #endif + if (max_pred_error >= pSPARC->F_tol_SOAP || !(Count_Sim % pSPARC->MLFF_DFT_fq)){ if (pSPARC->print_mlff_flag == 1 && rank ==0){ - fprintf(fp_mlff, "Error prediction in ML Model is too large!\n"); + if (max_pred_error >= pSPARC->F_tol_SOAP){ + fprintf(fp_mlff, "Error prediction in ML Model is too large!\n"); + } else { + fprintf(fp_mlff, "Max consecutive MLFF steps reached. Calling DFT.\n"); + } } t1 = MPI_Wtime(); Calculate_electronicGroundState(pSPARC); - if(pSPARC->cell_typ != 0){ - coordinatetransform_map(pSPARC, pSPARC->n_atom, pSPARC->atom_pos); + if (pSPARC->MDFlag == 1 || pSPARC->RelaxFlag == 1) { + if(pSPARC->cell_typ != 0){ + for(int i = 0; i < pSPARC->n_atom; i++) + Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); + } } - t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "DFT call made. Time taken: %.3f s\n", t2-t1); @@ -1112,22 +1002,17 @@ t2 = MPI_Wtime(); MPI_Bcast(pSPARC->forces, 3*pSPARC->n_atom, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(pSPARC->stress, 6, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(&pSPARC->pres, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); - mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count] = pSPARC->Etot-pSPARC->Entropy; mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count] = pSPARC->Etot; if (pSPARC->print_mlff_flag == 1 && rank ==0 && pSPARC->mlff_internal_energy_flag){ fprintf(fp_mlff, "Internal energy DFT data added! Free energy: %.6E Ha, Internal energy: %.6E Ha, Entropy: %.6E Ha\n", - mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count], mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count], pSPARC->Entropy); + mlff_str->free_energy_DFT[mlff_str->internal_energy_DFT_count], mlff_str->internal_energy_DFT[mlff_str->internal_energy_DFT_count], pSPARC->Entropy); } - mlff_str->internal_energy_DFT_count += 1; - double sum_square_error = 0.0; double max_force_error = 0.0; - if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Error Energy: %10.9f (Ha/atom)\n", fabs((1.0/pSPARC->n_atom)*pSPARC->Etot - E_predict[0])); - if(pSPARC->BC == 2){ if (pSPARC->mlff_pressure_train_flag==0){ for(int i = 0; i< mlff_str->stress_len; i++){ @@ -1136,14 +1021,10 @@ t2 = MPI_Wtime(); } else { fprintf(fp_mlff, "Pressure: Error (GPa) = %10.9f, Error(percent) = %10.9f\n", fabs(stress_predict[0]- pSPARC->pres)*au2GPa, fabs(stress_predict[0]- pSPARC->pres)/fabs(pSPARC->pres) *100); } - - - } else { for(int i = 0; i< mlff_str->stress_len; i++) fprintf(fp_mlff, "Stress[%d]: Error (Ha/Bohr^n) = %10.9f, Error(percent) = %10.9f\n", index[i], fabs(stress_predict[i]-pSPARC->stress[index[i]]), fabs(stress_predict[i]-pSPARC->stress[index[i]])/fabs(pSPARC->stress[index[i]]) *100); } - fprintf(fp_mlff, "Force Error atom-wise: \n"); sum_square_error = 0.0; max_force_error = 0.0; @@ -1151,7 +1032,6 @@ t2 = MPI_Wtime(); sum_square_error += (pSPARC->forces[3*i]+F_predict[3*i])*(pSPARC->forces[3*i]+F_predict[3*i]); sum_square_error += (pSPARC->forces[3*i+1]+F_predict[3*i+1])*(pSPARC->forces[3*i+1]+F_predict[3*i+1]); sum_square_error += (pSPARC->forces[3*i+2]+F_predict[3*i+2])*(pSPARC->forces[3*i+2]+F_predict[3*i+2]); - fprintf(fp_mlff, "%10.9f %10.9f %10.9f\n",pSPARC->forces[3*i]+F_predict[3*i], pSPARC->forces[3*i+1]+F_predict[3*i+1], pSPARC->forces[3*i+2]+F_predict[3*i+2]); @@ -1167,11 +1047,9 @@ t2 = MPI_Wtime(); } fprintf(fp_mlff, "Error Force (max): %10.9f (Ha/Bohr)\n", max_force_error); fprintf(fp_mlff, "Error Force (RMS): %10.9f (Ha/Bohr)\n", sqrt(sum_square_error/(3*pSPARC->n_atom))); - fprintf(fp_mlff, "Error Force Train (max): %10.9f (Ha/Bohr)\n", mlff_str->error_train_F); fprintf(fp_mlff, "Error Energy Train (max): %10.9f (Ha/atom)\n", mlff_str->error_train_E); } - init_dyarray(&mlff_str->atom_idx_addtrain); for (int i = 0; i < pSPARC->n_atom; i++){ if (bayesian_error[1+mlff_str->stress_len+3*i]>pSPARC->F_tol_SOAP || bayesian_error[1+mlff_str->stress_len+3*i+1]>pSPARC->F_tol_SOAP || bayesian_error[1+mlff_str->stress_len+3*i+2]>pSPARC->F_tol_SOAP){ @@ -1179,7 +1057,7 @@ t2 = MPI_Wtime(); } } t1 = MPI_Wtime(); - sparc_mlff_interface_addMD(pSPARC, mlff_str); + sparc_mlff_interface_addDFT(pSPARC, mlff_str); t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Covariance matrices updated! Time taken: %.3f s\n", t2-t1); @@ -1193,7 +1071,10 @@ t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Bayesian linear regression done! Time taken: %.3f s\n", t2-t1); } - pSPARC->last_train_MD_iter = pSPARC->MDCount; + pSPARC->last_train_iter = pSPARC->elecgs_Count; + #ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) pSPARC->last_train_iter = pSPARC->SocketSCFCount - 1; + #endif delete_dyarray(&mlff_str->atom_idx_addtrain); } else { #ifdef DEBUG @@ -1214,10 +1095,7 @@ t2 = MPI_Wtime(); E_internal = pSPARC->Etot*mlff_str->internal_energy_model_weights[1] + mlff_str->internal_energy_model_weights[0]; entropy = pSPARC->Etot - E_internal; } - - pSPARC->Entropy = entropy; - if (pSPARC->mlff_pressure_train_flag==0){ for(int i = 0; i < mlff_str->stress_len; i++){ pSPARC->stress[index[i]] = stress_predict[i]; @@ -1227,8 +1105,6 @@ t2 = MPI_Wtime(); } else { pSPARC->pres = stress_predict[0]; } - - if (rank==0){ write_MLFF_results(pSPARC); } @@ -1238,44 +1114,36 @@ t2 = MPI_Wtime(); free(stress_predict); free(bayesian_error); } - t4 = MPI_Wtime(); - } -/* -MLFF_call_from_MD_only_predict function performs MLFF caclulation only - -[Input] -1. pSPARC: SPARC object -2. mlff_str: MLFF object -[Output] -1. pSPARC: SPARC object -2. mlff_str: MLFF object +/** + * @brief Performs MLFF calculation without UQ for DFT check */ -void MLFF_call_from_MD_only_predict(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str){ + +void MLFF_only_predict(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str){ int nprocs, rank; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - double t1, t2; FILE *fp_mlff; if (pSPARC->print_mlff_flag==1 && rank==0){ // fp_mlff = fopen("mlff.log","a"); fp_mlff = mlff_str->fp_mlff; } - if (pSPARC->print_mlff_flag==1 && rank==0){ fprintf(fp_mlff,"--------------------------------------------------------------------------------------------\n"); - fprintf(fp_mlff, "MD step: %d\n", pSPARC->MDCount); + int print_count = pSPARC->elecgs_Count; + #ifdef USE_SOCKET + if (pSPARC->SocketFlag == 1) print_count = pSPARC->SocketSCFCount - 1; + #endif + fprintf(fp_mlff, "Calculation: %d\n", print_count); } - double *E_predict, *F_predict, *bayesian_error, *stress_predict; E_predict = (double *)malloc(1*sizeof(double)); F_predict = (double *)malloc(3*pSPARC->n_atom*sizeof(double)); stress_predict = (double *)malloc(mlff_str->stress_len*sizeof(double)); bayesian_error = (double *)malloc((3*pSPARC->n_atom+1+mlff_str->stress_len)*sizeof(double)); - t1 = MPI_Wtime(); sparc_mlff_interface_predict(pSPARC, mlff_str, E_predict, F_predict, stress_predict, bayesian_error); t2 = MPI_Wtime(); @@ -1286,11 +1154,9 @@ t2 = MPI_Wtime(); for (int i=0; i < pSPARC->n_atom*3; i++){ pSPARC->forces[i] = -1.0*F_predict[i]; } - int index[] = {0,1,2,3,4,5}; int BC[] = {pSPARC->BCx, pSPARC->BCy, pSPARC->BCz}; reshape_stress(pSPARC->cell_typ, BC, index); - if (pSPARC->mlff_pressure_train_flag){ for(int i = 0; i < mlff_str->stress_len; i++){ pSPARC->stress[index[i]] = stress_predict[i]; @@ -1300,13 +1166,10 @@ t2 = MPI_Wtime(); } else { pSPARC->pres = stress_predict[0]; } - - free(E_predict); free(F_predict); free(bayesian_error); free(stress_predict); - if (rank==0){ write_MLFF_results(pSPARC); } @@ -1315,45 +1178,38 @@ t2 = MPI_Wtime(); } } -void sparc_mlff_interface_firstMD(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str){ - +/** + * @brief This function is used to add the first DFT step to the MLFF training set + */ +void sparc_mlff_interface_initialDFT(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str){ double t1, t2, t3, t4; t3 = MPI_Wtime(); - int rank, nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - FILE *fp_mlff; if (pSPARC->print_mlff_flag==1 && rank==0){ - // fp_mlff = fopen("mlff.log","a"); fp_mlff = mlff_str->fp_mlff; } - int *Z = (int *)malloc(pSPARC->Ntypes*sizeof(int)); for (int i=0; i Ntypes; i++){ Z[i] = pSPARC->Znucl[i]; } - double F_estimate_max; double **K_predict, *K_predict_col_major; - double *geometric_ratio = (double*) malloc(2 * sizeof(double)); geometric_ratio[0] = pSPARC->CUTOFF_y[0]/pSPARC->CUTOFF_x[0]; geometric_ratio[1] = pSPARC->CUTOFF_z[0]/pSPARC->CUTOFF_x[0]; - int *BC = (int*)malloc(3*sizeof(int)); double *cell =(double *) malloc(3*sizeof(double)); int *atomtyp = (int *) malloc(pSPARC->n_atom*sizeof(int)); - BC[0] = pSPARC->BCx; BC[1] = pSPARC->BCy; BC[2] = pSPARC->BCz; cell[0] = pSPARC->range_x; cell[1] = pSPARC->range_y; cell[2] = pSPARC->range_z; - int count = 0; for (int i=0; i < pSPARC->Ntypes; i++){ for (int j=0; j < pSPARC->nAtomv[i]; j++){ @@ -1361,22 +1217,16 @@ t3 = MPI_Wtime(); count++; } } - t1 = MPI_Wtime(); NeighList *nlist = (NeighList *) malloc(sizeof(NeighList)*1); build_nlist(mlff_str->rcut, pSPARC->Ntypes, pSPARC->n_atom, pSPARC->atom_pos, atomtyp, pSPARC->cell_typ, BC, cell, pSPARC->LatUVec, pSPARC->twist, geometric_ratio, nlist, mlff_str->natom_domain, mlff_str->atom_idx_domain, mlff_str->el_idx_domain); DescriptorObj *desc_str = (DescriptorObj *) malloc(sizeof(DescriptorObj)*1); build_descriptor(desc_str, nlist, mlff_str, pSPARC->atom_pos); - - - t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "neighbor list and descriptor done! Time taken: %.3f s\n", t2-t1); } - t1 = MPI_Wtime(); - double *force = (double *)malloc(3*pSPARC->n_atom*sizeof(double)); for (int i=0; i < 3*pSPARC->n_atom; i++){ force[i] = -1.0*pSPARC->forces[i]; @@ -1386,17 +1236,13 @@ t1 = MPI_Wtime(); reshape_stress(pSPARC->cell_typ, BC, index); for(int i = 0; i < 6; i++) stress[i] = pSPARC->stress[index[i]]; - if (pSPARC->mlff_pressure_train_flag){ for (int i = 0; i < 6; i++){ index[i] = i; } stress[0]=pSPARC->pres; } - - - add_firstMD(desc_str, nlist, mlff_str, pSPARC->Etot/pSPARC->n_atom, force, stress); - + add_firstDFT(desc_str, nlist, mlff_str, pSPARC->Etot/pSPARC->n_atom, force, stress); if (rank==0){ if(pSPARC->cell_typ != 0){ for(int i = 0; i < pSPARC->n_atom; i++) @@ -1408,18 +1254,15 @@ t1 = MPI_Wtime(); } print_new_ref_structure_MLFF(mlff_str, mlff_str->n_str, nlist, pSPARC->atom_pos, pSPARC->Etot, force, pSPARC->stress); if(pSPARC->cell_typ != 0){ - coordinatetransform_map(pSPARC, pSPARC->n_atom, pSPARC->atom_pos); + for(int i = 0; i < pSPARC->n_atom; i++) + Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); } - } - t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Covariance matrix updated! Time taken: %.3f s\n", t2-t1); // fclose(fp_mlff); } - free(geometric_ratio); free(BC); free(cell); @@ -1431,61 +1274,48 @@ t2 = MPI_Wtime(); free(force); free(stress); free(Z); -t4 = MPI_Wtime(); - - - } +/** + * @brief Add DFT data to the MLFF model + */ -void sparc_mlff_interface_addMD(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str){ - +void sparc_mlff_interface_addDFT(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str){ double t1, t2, t3, t4; t3 = MPI_Wtime(); - int rank, nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - FILE *fp_mlff; if (pSPARC->print_mlff_flag==1 && rank==0){ // fp_mlff = fopen("mlff.log","a"); fp_mlff = mlff_str->fp_mlff; } - NeighList *nlist; DescriptorObj *desc_str; - // printf("sparc_mlff_interface_initialMD calld in pSPARC->MDCount: %d\n", pSPARC->MDCount); int *Z; Z = (int *)malloc(pSPARC->Ntypes*sizeof(int)); for (int i=0; i Ntypes; i++){ Z[i] = pSPARC->Znucl[i]; } - nlist = (NeighList *) malloc(sizeof(NeighList)*1); desc_str = (DescriptorObj *) malloc(sizeof(DescriptorObj)*1); - int *BC, *atomtyp; double *cell, F_estimate_max; double **K_predict, *K_predict_col_major; int nelem = mlff_str->nelem; int size_X3 = mlff_str->size_X3; double xi_3 = mlff_str->xi_3; - double *geometric_ratio = (double*) malloc(2 * sizeof(double)); geometric_ratio[0] = pSPARC->CUTOFF_y[0]/pSPARC->CUTOFF_x[0]; geometric_ratio[1] = pSPARC->CUTOFF_z[0]/pSPARC->CUTOFF_x[0]; - BC = (int*)malloc(3*sizeof(int)); cell =(double *) malloc(3*sizeof(double)); atomtyp = (int *) malloc(pSPARC->n_atom*sizeof(int)); - BC[0] = pSPARC->BCx; BC[1] = pSPARC->BCy; BC[2] = pSPARC->BCz; // periodic cell[0] = pSPARC->range_x; cell[1] = pSPARC->range_y; cell[2] = pSPARC->range_z; - - int count = 0; for (int i=0; i < pSPARC->Ntypes; i++){ for (int j=0; j < pSPARC->nAtomv[i]; j++){ @@ -1493,20 +1323,14 @@ t3 = MPI_Wtime(); count++; } } - - t1 = MPI_Wtime(); - build_nlist(mlff_str->rcut, pSPARC->Ntypes, pSPARC->n_atom, pSPARC->atom_pos, atomtyp, pSPARC->cell_typ, BC, cell, pSPARC->LatUVec, pSPARC->twist, geometric_ratio, nlist, mlff_str->natom_domain, mlff_str->atom_idx_domain, mlff_str->el_idx_domain); - build_descriptor(desc_str, nlist, mlff_str, pSPARC->atom_pos); - t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "neighbor list and soap descriptor done! Time taken: %.3f s\n", t2-t1); } t1 = MPI_Wtime(); - double *force = (double *)malloc(3*pSPARC->n_atom*sizeof(double)); for (int i=0; i < 3*pSPARC->n_atom; i++){ force[i] = -1.0*pSPARC->forces[i]; @@ -1516,39 +1340,29 @@ t1 = MPI_Wtime(); reshape_stress(pSPARC->cell_typ, BC, index); for(int i = 0; i < 6; i++) stress[i] = pSPARC->stress[index[i]]; - if (pSPARC->mlff_pressure_train_flag){ for (int i = 0; i < 6; i++){ index[i] = i; } stress[0]=pSPARC->pres; } - add_newstr_rows(desc_str, nlist, mlff_str, pSPARC->Etot/pSPARC->n_atom, force, stress); t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Covariance matrix rows updated! Time taken: %.3f s\n", t2-t1); } - t1 = MPI_Wtime(); - double *X3_gathered; - X3_gathered = (double *) malloc(sizeof(double)*size_X3*mlff_str->natom); - double *X3_local; - X3_local = (double *) malloc(sizeof(double)*size_X3*mlff_str->natom_domain); - for (int i=0; i < mlff_str->natom_domain; i++){ for (int j=0; j < size_X3; j++){ X3_local[i*size_X3+j] = desc_str->X3[i][j]; } } - int local_natoms[nprocs]; MPI_Allgather(&mlff_str->natom_domain, 1, MPI_INT, local_natoms, 1, MPI_INT, MPI_COMM_WORLD); - int recvcounts_X3[nprocs], displs_X3[nprocs]; displs_X3[0] = 0; for (int i=0; i < nprocs; i++){ @@ -1557,9 +1371,7 @@ t1 = MPI_Wtime(); displs_X3[i] = displs_X3[i-1]+local_natoms[i-1]*size_X3; } } - MPI_Allgatherv(X3_local, size_X3*mlff_str->natom_domain, MPI_DOUBLE, X3_gathered, recvcounts_X3, displs_X3, MPI_DOUBLE, MPI_COMM_WORLD); - double **X3_gathered_2D; X3_gathered_2D = (double **) malloc(sizeof(double*)*mlff_str->natom); for (int i=0; i < mlff_str->natom; i++){ @@ -1568,24 +1380,20 @@ t1 = MPI_Wtime(); X3_gathered_2D[i][j] = X3_gathered[i*size_X3+j]; } } - double **X3_toadd; int no_desc_toadd[nelem]; for (int i=0; i n_atom; i++){ if (lin_search(((&mlff_str->atom_idx_addtrain))->array, ((&mlff_str->atom_idx_addtrain))->len, i) != -1){ no_desc_toadd[atomtyp[i]] += 1; } } - // Only atoms whose force error was larger than the threshold X3_toadd = (double **) malloc(sizeof(double *)* (((&mlff_str->atom_idx_addtrain))->len)); int atom_typ1[((&mlff_str->atom_idx_addtrain))->len]; - for (int i =0; i<(((&mlff_str->atom_idx_addtrain))->len); i++){ X3_toadd[i] = (double *) malloc(sizeof(double)*size_X3); atom_typ1[i] = atomtyp[((&mlff_str->atom_idx_addtrain))->array[i]]; @@ -1593,17 +1401,14 @@ t1 = MPI_Wtime(); X3_toadd[i][j] = X3_gathered_2D[((&mlff_str->atom_idx_addtrain))->array[i]][j]; } } - int *cum_natm_elem; cum_natm_elem = (int *)malloc(sizeof(int)*nelem); cum_natm_elem[0] = 0; for (int i = 1; i < nelem; i++){ cum_natm_elem[i] = cum_natm_elem[i-1]+no_desc_toadd[i-1]; } - dyArray *highrank_ID_descriptors; highrank_ID_descriptors = (dyArray *) malloc(sizeof(dyArray)*nelem); - int total_descriptors_toadd = 0; for (int i=0; i print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Sparsification done! Time taken: %.3f s\n", t2-t1); } t1 = MPI_Wtime(); - int temp_idx; for (int i=0; i < nelem; i++){ for (int j=0; j<(highrank_ID_descriptors[i]).len; j++){ @@ -1628,7 +1431,6 @@ t1 = MPI_Wtime(); add_newtrain_cols(X3_toadd[temp_idx], i, mlff_str); } } - t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Covariance matrix columns updated! Time taken: %.3f s\n", t2-t1); @@ -1644,18 +1446,16 @@ t2 = MPI_Wtime(); } print_new_ref_structure_MLFF(mlff_str, mlff_str->n_str, nlist, pSPARC->atom_pos, pSPARC->Etot, force, pSPARC->stress); if(pSPARC->cell_typ != 0){ - coordinatetransform_map(pSPARC, pSPARC->n_atom, pSPARC->atom_pos); + for(int i = 0; i < pSPARC->n_atom; i++) + Cart2nonCart_coord(pSPARC, &pSPARC->atom_pos[3*i], &pSPARC->atom_pos[3*i+1], &pSPARC->atom_pos[3*i+2]); } } - for (int i=0; i atom_idx_addtrain))->len); i++){ free(X3_toadd[i]); } - free(X3_toadd); free(cum_natm_elem); free(X3_local); @@ -1668,9 +1468,7 @@ t2 = MPI_Wtime(); free(BC); free(cell); free(atomtyp); - delete_descriptor(desc_str); - clear_nlist(nlist, mlff_str->natom_domain); free(nlist); free(desc_str); @@ -1678,54 +1476,46 @@ t2 = MPI_Wtime(); free(stress); free(highrank_ID_descriptors); free(Z); - t4 = MPI_Wtime(); } +/** + * @brief Predict energies, forces, stresses and Bayesian error from MLFF + */ + void sparc_mlff_interface_predict(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str, double *E_predict, double *F_predict, double *stress_predict, double *bayesian_error){ double t1, t2, t3, t4; - t3 = MPI_Wtime(); - int nprocs, rank; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - FILE *fp_mlff; if (pSPARC->print_mlff_flag==1 && rank==0){ // fp_mlff = fopen("mlff.log","a"); fp_mlff = mlff_str->fp_mlff; } - NeighList *nlist; DescriptorObj *desc_str; - nlist = (NeighList *) malloc(sizeof(NeighList)*1); desc_str = (DescriptorObj *) malloc(sizeof(DescriptorObj)*1); - int *BC, *atomtyp; double *cell, F_estimate_max; double **K_predict, *K_predict_col_major; - int *Z; Z = (int *)malloc(pSPARC->Ntypes*sizeof(int)); for (int i=0; i Ntypes; i++){ Z[i] = pSPARC->Znucl[i]; } - double *geometric_ratio = (double*) malloc(2 * sizeof(double)); geometric_ratio[0] = pSPARC->CUTOFF_y[0]/pSPARC->CUTOFF_x[0]; geometric_ratio[1] = pSPARC->CUTOFF_z[0]/pSPARC->CUTOFF_x[0]; - BC = (int*)malloc(3*sizeof(int)); cell =(double *) malloc(3*sizeof(double)); atomtyp = (int *) malloc(pSPARC->n_atom*sizeof(int)); - BC[0] = pSPARC->BCx; BC[1] = pSPARC->BCy; BC[2] = pSPARC->BCz; cell[0] = pSPARC->range_x; cell[1] = pSPARC->range_y; cell[2] = pSPARC->range_z; - int count = 0; for (int i=0; i < pSPARC->Ntypes; i++){ for (int j=0; j < pSPARC->nAtomv[i]; j++){ @@ -1733,18 +1523,14 @@ t3 = MPI_Wtime(); count++; } } - t1 = MPI_Wtime(); build_nlist(mlff_str->rcut, pSPARC->Ntypes, pSPARC->n_atom, pSPARC->atom_pos, atomtyp, pSPARC->cell_typ, BC, cell, pSPARC->LatUVec, pSPARC->twist, geometric_ratio, nlist, mlff_str->natom_domain, mlff_str->atom_idx_domain, mlff_str->el_idx_domain); build_descriptor(desc_str, nlist, mlff_str, pSPARC->atom_pos); t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Neighbor list and soap descriptor done! Time taken: %.3f s\n", t2-t1); } - t1 = MPI_Wtime(); - if (rank==0){ K_predict = (double **)malloc(sizeof(double*)*(1+3*mlff_str->natom_domain+mlff_str->stress_len)); K_predict_col_major = (double *)malloc(sizeof(double)*(1+3*mlff_str->natom_domain+mlff_str->stress_len)*mlff_str->n_cols); @@ -1765,9 +1551,7 @@ t1 = MPI_Wtime(); } } - calculate_Kpredict(desc_str, nlist, mlff_str, K_predict); - if (rank==0){ for (int j=0; j < mlff_str->n_cols; j++){ for (int i=0; i<(1+3*mlff_str->natom_domain+mlff_str->stress_len); i++){ @@ -1781,33 +1565,26 @@ t1 = MPI_Wtime(); } } } - double *F_predict_local, *bayesian_error_local; - F_predict_local = (double *) malloc(sizeof(double)*3*mlff_str->natom_domain); if (rank==0){ bayesian_error_local = (double *) malloc(sizeof(double)*(1+3*mlff_str->natom_domain+mlff_str->stress_len)); } else{ bayesian_error_local = (double *) malloc(sizeof(double)*(3*mlff_str->natom_domain)); } - t2 = MPI_Wtime(); if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "Kpredict calculation done done! Time taken: %.3f s\n", t2-t1); } t1 = MPI_Wtime(); - mlff_predict(K_predict_col_major, mlff_str, E_predict, F_predict_local, stress_predict, bayesian_error_local, pSPARC->n_atom); MPI_Bcast(E_predict, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(stress_predict, mlff_str->stress_len, MPI_DOUBLE, 0, MPI_COMM_WORLD); - int local_natoms[nprocs]; MPI_Allgather(&mlff_str->natom_domain, 1, MPI_INT, local_natoms, 1, MPI_INT, MPI_COMM_WORLD); - int recvcounts[nprocs], displs[nprocs], recvcounts_bayesian[nprocs], displs_bayesian[nprocs]; displs[0] = 0; displs_bayesian[0] = 0; - for (int i=0; i < nprocs; i++){ recvcounts[i] = local_natoms[i]*3; if (i==0){ @@ -1824,17 +1601,14 @@ t1 = MPI_Wtime(); } } } - int send_bayesian; if (rank==0){ send_bayesian = 3*mlff_str->natom_domain+1+mlff_str->stress_len; } else { send_bayesian = 3*mlff_str->natom_domain; } - MPI_Allgatherv(F_predict_local, 3*mlff_str->natom_domain, MPI_DOUBLE, F_predict, recvcounts, displs, MPI_DOUBLE, MPI_COMM_WORLD); MPI_Allgatherv(bayesian_error_local, send_bayesian, MPI_DOUBLE, bayesian_error, recvcounts_bayesian, displs_bayesian, MPI_DOUBLE, MPI_COMM_WORLD); - if (rank==0){ for (int i=0; i<(1+mlff_str->stress_len+3*mlff_str->natom_domain); i++){ free(K_predict[i]); @@ -1845,12 +1619,10 @@ t1 = MPI_Wtime(); } } t2 = MPI_Wtime(); - if (pSPARC->print_mlff_flag == 1 && rank ==0){ fprintf(fp_mlff, "prediction and gather done! Time taken: %.3f s\n", t2-t1); // fclose(fp_mlff); } - free(bayesian_error_local); free(F_predict_local); free(K_predict); @@ -1860,7 +1632,6 @@ t2 = MPI_Wtime(); free(cell); free(atomtyp); delete_descriptor(desc_str); - clear_nlist(nlist, mlff_str->natom_domain); free(nlist); free(desc_str); @@ -1868,15 +1639,17 @@ t2 = MPI_Wtime(); t4 = MPI_Wtime(); } +/** + * @brief Write the MLFF results to the output and static files + */ + +// TODO: Implement static printing void write_MLFF_results(SPARC_OBJ *pSPARC){ - int rank, nproc, i; FILE *output_fp, *static_fp; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &nproc); - // write energies into output file - if (!rank && pSPARC->Verbosity) { output_fp = fopen(pSPARC->OutFilename,"a"); if (output_fp == NULL) { @@ -1884,7 +1657,16 @@ void write_MLFF_results(SPARC_OBJ *pSPARC){ exit(EXIT_FAILURE); } fprintf(output_fp,"====================================================================\n"); - fprintf(output_fp," Energy and force calculation (MLFF #%d) \n", pSPARC->MDCount + pSPARC->restartCount + (pSPARC->RestartFlag == 0)); + if (pSPARC->MDFlag == 1) + fprintf(output_fp," Energy and force calculation (MLFF #%d) \n", pSPARC->MDCount + pSPARC->restartCount + (pSPARC->RestartFlag == 0)); + else if (pSPARC->RelaxFlag >= 1) + fprintf(output_fp," Energy and force calculation (MLFF #%d) \n", pSPARC->RelaxCount + pSPARC->restartCount + (pSPARC->RestartFlag == 0)); +#ifdef USE_SOCKET + else if (pSPARC->SocketFlag == 1) + fprintf(output_fp," Energy and force calculation (MLFF #%d) \n", pSPARC->SocketSCFCount - 1 + pSPARC->restartCount + (pSPARC->RestartFlag == 0)); +#endif + else + fprintf(output_fp," Energy and force calculation (MLFF #%d) \n", 1); fprintf(output_fp,"====================================================================\n"); fprintf(output_fp,"Free energy per atom :%18.10E (Ha/atom)\n", pSPARC->Etot / pSPARC->n_atom); fprintf(output_fp,"Total free energy :%18.10E (Ha)\n", pSPARC->Etot); @@ -1901,8 +1683,35 @@ void write_MLFF_results(SPARC_OBJ *pSPARC){ fclose(static_fp); } } + // write forces into .static file if required + if (pSPARC->PrintForceFlag == 1 && pSPARC->MDFlag == 0 && pSPARC->RelaxFlag == 0) { + static_fp = fopen(pSPARC->StaticFilename,"a"); + if (static_fp == NULL) { + printf("\nCannot open file \"%s\"\n",pSPARC->StaticFilename); + exit(EXIT_FAILURE); + } + fprintf(static_fp, "Atomic forces (Ha/Bohr):\n"); + for (i = 0; i < pSPARC->n_atom; i++) { + fprintf(static_fp,"%18.10E %18.10E %18.10E\n", + pSPARC->forces[3*i],pSPARC->forces[3*i+1],pSPARC->forces[3*i+2]); + } + fclose(static_fp); + } + // Calculate Stress and pressure + if(pSPARC->Calc_stress == 1){ + // write stress to .static file + if (pSPARC->MDFlag == 0 && pSPARC->RelaxFlag == 0) { + static_fp = fopen(pSPARC->StaticFilename,"a"); + if (static_fp == NULL) { + printf("\nCannot open file \"%s\"\n",pSPARC->StaticFilename); + exit(EXIT_FAILURE); + } + fprintf(static_fp,"Stress"); + PrintStress(pSPARC, pSPARC->stress, static_fp); + fclose(static_fp); + } + } } - if(!rank && pSPARC->Verbosity) { output_fp = fopen(pSPARC->OutFilename,"a"); if (output_fp == NULL) { @@ -1922,7 +1731,6 @@ void write_MLFF_results(SPARC_OBJ *pSPARC){ fprintf(output_fp,"Maximum force :%18.10E (Ha/Bohr)\n",maxF); fclose(output_fp); } - output_fp = fopen(pSPARC->OutFilename,"a"); int stress_len = (1-pSPARC->BCx) + (1-pSPARC->BCy) + (1-pSPARC->BCz) + ( (1-pSPARC->BCx-pSPARC->BCy) > 0) + ( (1-pSPARC->BCx-pSPARC->BCz) > 0) + ( (1-pSPARC->BCy-pSPARC->BCz) > 0); stress_len = (pSPARC->cell_typ > 20 && pSPARC->cell_typ < 30) ? 1 : stress_len; @@ -1932,7 +1740,6 @@ void write_MLFF_results(SPARC_OBJ *pSPARC){ reshape_stress(pSPARC->cell_typ, BC, index); for(int i = 0; i < 6; i++) stress[i] = pSPARC->stress[index[i]]; - double maxS = 0.0, temp; for (i = 0; i < stress_len; i++){ temp = fabs(stress[i]); @@ -1947,27 +1754,24 @@ void write_MLFF_results(SPARC_OBJ *pSPARC){ fprintf(output_fp,"Maximum Stress :%18.10E (Ha/Bohr^n) (n = periodicity order)\n",maxS); } } - if(pSPARC->Calc_pres == 1){ fprintf(output_fp,"Pressure :%18.10E (GPa)\n",pSPARC->pres*CONST_HA_BOHR3_GPA); } - if (pSPARC->mlff_internal_energy_flag){ fprintf(output_fp,"-Entropy*kb*T :%18.10E (Ha)\n", pSPARC->Entropy); } fclose(output_fp); - free(stress); } +/** + * @brief Reshape the stress tensor based on the cell type and boundary conditions for MLFF training + */ void reshape_stress(int cell_typ, int *BC, int *index) { - int rank, nprocs; MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &rank); - - int count = 0, skip_count=0; if(cell_typ < 20){ for(int i = 0; i < 3; i++){ @@ -1981,21 +1785,4 @@ void reshape_stress(int cell_typ, int *BC, int *index) { } else{ index[0] = 5; } - - // if (mlff_pressure_train_flag){ - // for (int i = 0; i < 6; i++){ - // index[i] = i; - // } - // } -} - - - -void coordinatetransform_map(SPARC_OBJ *pSPARC, int natom, double *coord) { - // Convert Cart to nonCart coordinates for non orthogonal cell - if(pSPARC->cell_typ != 0){ - for(int atm = 0; atm < natom; atm++) - Cart2nonCart_coord(pSPARC, &coord[3*atm], &coord[3*atm+1], &coord[3*atm+2]); - } - -} +} \ No newline at end of file diff --git a/src/mlff/sparc_mlff_interface.h b/src/mlff/sparc_mlff_interface.h index fca8b777..98b5c5be 100644 --- a/src/mlff/sparc_mlff_interface.h +++ b/src/mlff/sparc_mlff_interface.h @@ -3,21 +3,80 @@ #include "isddft.h" #include "mlff_types.h" + +/** + * @brief Governing function that calls appropriate subroutines based on the MLFF flag and state of simulation + */ + +void MLFF_call(SPARC_OBJ* pSPARC); + +/** + * @brief Initializes the MLFF structure. + */ // Updated -void init_MLFF(SPARC_OBJ *pSPARC); +void init_MLFF_structure(SPARC_OBJ *pSPARC); + +/** + * @brief Initializes a new MLFF simulation depending on the flag. + */ +// New +void init_MLFF_simulation(SPARC_OBJ *pSPARC); + +/** + * @brief Initializes a new MLFF simulation from scratch. + */ + +void init_new_MLFF(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); + +/** + * @brief Initializes a new MLFF simulation from a known model that can be updated. + */ + +void init_existing_trainable_MLFF(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); + +/** + * @brief Initializes a new MLFF simulation from a known model that remains static. + */ + +void init_static_MLFF(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); + +/** + * @brief Frees memory allocated for MLFF structure. + */ // Updated void free_MLFF(MLFF_Obj *mlff_str); -// Updated -void MLFF_main(SPARC_OBJ *pSPARC); + /** + * @brief Constructs design matrices from saved reference structures and reference atomic descriptors. + */ -void get_domain_decompose_mlff_natom(int natom, +void pretrain_MLFF_model( + MLFF_Obj *mlff_str, + SPARC_OBJ *pSPARC, + double **cell_data, + double **LatUVec_data, + double **apos_data, + double *Etot_data, + double **F_data, + double **stress_data, + int *natom_data, + int **natom_elem_data); + + /** + * @brief Gets the number of atoms on each processor + */ + + void get_domain_decompose_mlff_natom(int natom, int nelem, int *nAtomv, int nprocs, int rank, int *natom_domain); + /** + * @brief Returns the atom and element indices for parallelization in MLFF + */ + void get_domain_decompose_mlff_idx( int natom, int nelem, @@ -28,42 +87,46 @@ void get_domain_decompose_mlff_idx( int *atom_idx_domain, int *el_idx_domain); -void pretrain_MLFF_model( - MLFF_Obj *mlff_str, - SPARC_OBJ *pSPARC, - double **cell_data, - double **LatUVec_data, - double **apos_data, - double *Etot_data, - double **F_data, - double **stress_data, - int *natom_data, - int **natom_elem_data); +/** + * @brief Returns predictions from MLFF or DFT depending on size of training set or UQ. Updates MLFF model as necessary. + */ +void MLFF_train_predict(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); -void MLFF_call(SPARC_OBJ* pSPARC); - - -void MLFF_call_from_MD(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); +/** + * @brief Performs MLFF calculation without UQ for DFT check +*/ +void MLFF_only_predict(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); -void MLFF_call_from_MD_only_predict(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); +/** + * @brief This function is used to add the first MD step to the MLFF training set + */ +void sparc_mlff_interface_initialDFT(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); -void sparc_mlff_interface_firstMD(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); +/** + * @brief Add DFT data to the MLFF model + */ +void sparc_mlff_interface_addDFT(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); -void sparc_mlff_interface_addMD(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str); - +/** + * @brief Predict energies, forces, stresses and Bayesian error from MLFF + */ void sparc_mlff_interface_predict(SPARC_OBJ *pSPARC, MLFF_Obj *mlff_str, double *E_predict, double *F_predict, double *stress_predict, double *bayesian_error); +/** + * @brief Write the MLFF results to the output and static files + */ void write_MLFF_results(SPARC_OBJ *pSPARC); +/** + * @brief Reshape the stress tensor based on the cell type and boundary conditions for MLFF training + */ void reshape_stress(int cell_typ, int *BC, int *index); -void coordinatetransform_map(SPARC_OBJ *pSPARC, int natom, double *coord); - #endif \ No newline at end of file diff --git a/src/readfiles.c b/src/readfiles.c index 946db6c3..a0947fc4 100644 --- a/src/readfiles.c +++ b/src/readfiles.c @@ -647,10 +647,7 @@ void read_input(SPARC_INPUT_OBJ *pSPARC_Input, SPARC_OBJ *pSPARC) { } else if(strcmpi(str,"MLFF_IF_ATOM_DATA_AVAILABLE:") == 0) { fscanf(input_fp,"%d",&pSPARC_Input->if_atom_data_available); fscanf(input_fp, "%*[^\n]\n"); - } else if (strcmpi(str,"MLFF_hnl_FILE:") == 0) { - fscanf(input_fp,"%s",pSPARC_Input->hnl_file_name); - fscanf(input_fp, "%*[^\n]\n"); - } else if (strcmpi(str,"MLFF_MODEL_FOLDER:") == 0) { // MLFF start + } else if (strcmpi(str,"MLFF_MODEL_FOLDER:") == 0) { // MLFF start fscanf(input_fp,"%s",pSPARC_Input->mlff_data_folder); fscanf(input_fp, "%*[^\n]\n"); } else if(strcmpi(str,"MLFF_DESCRIPTOR_TYPE:") == 0) { diff --git a/tests/mlff/AlSi_NPTNH_mlff/high_accuracy/AlSi_NPTNH_mlff.inpt b/tests/mlff/AlSi_NPTNH_mlff/high_accuracy/AlSi_NPTNH_mlff.inpt index e53a887a..4e3e5abf 100644 --- a/tests/mlff/AlSi_NPTNH_mlff/high_accuracy/AlSi_NPTNH_mlff.inpt +++ b/tests/mlff/AlSi_NPTNH_mlff/high_accuracy/AlSi_NPTNH_mlff.inpt @@ -31,7 +31,7 @@ MLFF_RADIAL_BASIS: 8 # N_R (Number of raddial basis functions in SOAP) MLFF_RCUT_SOAP: 10 # (Cutoff radius for the SOAP descriptor) MLFF_MODEL_FOLDER: ./ -MLFF_INITIAL_STEPS_TRAIN: 3 # Initial number of MD steps where DFT calculation is imposed +MLFF_INITIAL_STEPS_TRAIN: 4 # Initial number of MD steps where DFT calculation is imposed MLFF_MAX_STR_STORE: 500 # This is the maximum number of structure for which memory is allocated to store MLFF data MLFF_MAX_CONFIG_STORE: 5000 # THis is the maximum number of atomic descriptors per element type that is stored MLFF_SCALE_FORCE: 1.0 # Scaling of force relative to energy diff --git a/tests/mlff/AlSi_NPTNH_mlff/standard/AlSi_NPTNH_mlff.inpt b/tests/mlff/AlSi_NPTNH_mlff/standard/AlSi_NPTNH_mlff.inpt index d1b698c8..99dc4e8c 100644 --- a/tests/mlff/AlSi_NPTNH_mlff/standard/AlSi_NPTNH_mlff.inpt +++ b/tests/mlff/AlSi_NPTNH_mlff/standard/AlSi_NPTNH_mlff.inpt @@ -27,7 +27,7 @@ MLFF_RADIAL_BASIS: 8 # N_R (Number of raddial basis functions in SOAP) MLFF_RCUT_SOAP: 10 # (Cutoff radius for the SOAP descriptor) MLFF_MODEL_FOLDER: ./ -MLFF_INITIAL_STEPS_TRAIN: 3 # Initial number of MD steps where DFT calculation is imposed +MLFF_INITIAL_STEPS_TRAIN: 4 # Initial number of MD steps where DFT calculation is imposed MLFF_MAX_STR_STORE: 500 # This is the maximum number of structure for which memory is allocated to store MLFF data MLFF_MAX_CONFIG_STORE: 5000 # THis is the maximum number of atomic descriptors per element type that is stored MLFF_SCALE_FORCE: 1.0 # Scaling of force relative to energy diff --git a/tests/mlff/AlSi_mlff_nonorth/high_accuracy/AlSi_mlff_nonorth.inpt b/tests/mlff/AlSi_mlff_nonorth/high_accuracy/AlSi_mlff_nonorth.inpt index a77bc5d4..ef88b339 100644 --- a/tests/mlff/AlSi_mlff_nonorth/high_accuracy/AlSi_mlff_nonorth.inpt +++ b/tests/mlff/AlSi_mlff_nonorth/high_accuracy/AlSi_mlff_nonorth.inpt @@ -32,7 +32,7 @@ MLFF_RADIAL_BASIS: 8 # N_R (Number of raddial basis functions in SOAP) MLFF_RCUT_SOAP: 10 # (Cutoff radius for the SOAP descriptor) MLFF_MODEL_FOLDER: ./ -MLFF_INITIAL_STEPS_TRAIN: 3 # Initial number of MD steps where DFT calculation is imposed +MLFF_INITIAL_STEPS_TRAIN: 4 # Initial number of MD steps where DFT calculation is imposed MLFF_MAX_STR_STORE: 500 # This is the maximum number of structure for which memory is allocated to store MLFF data MLFF_MAX_CONFIG_STORE: 5000 # THis is the maximum number of atomic descriptors per element type that is stored MLFF_SCALE_FORCE: 1.0 # Scaling of force relative to energy diff --git a/tests/mlff/AlSi_mlff_nonorth/standard/AlSi_mlff_nonorth.inpt b/tests/mlff/AlSi_mlff_nonorth/standard/AlSi_mlff_nonorth.inpt index 28a6ce45..77ed46fa 100644 --- a/tests/mlff/AlSi_mlff_nonorth/standard/AlSi_mlff_nonorth.inpt +++ b/tests/mlff/AlSi_mlff_nonorth/standard/AlSi_mlff_nonorth.inpt @@ -31,7 +31,7 @@ MLFF_RADIAL_BASIS: 8 # N_R (Number of raddial basis functions in SOAP) MLFF_RCUT_SOAP: 10 # (Cutoff radius for the SOAP descriptor) MLFF_MODEL_FOLDER: ./ -MLFF_INITIAL_STEPS_TRAIN: 3 # Initial number of MD steps where DFT calculation is imposed +MLFF_INITIAL_STEPS_TRAIN: 4 # Initial number of MD steps where DFT calculation is imposed MLFF_MAX_STR_STORE: 500 # This is the maximum number of structure for which memory is allocated to store MLFF data MLFF_MAX_CONFIG_STORE: 5000 # THis is the maximum number of atomic descriptors per element type that is stored MLFF_SCALE_FORCE: 1.0 # Scaling of force relative to energy diff --git a/tests/mlff/AlSi_mlffflag1/high_accuracy/AlSi_mlffflag1.inpt b/tests/mlff/AlSi_mlffflag1/high_accuracy/AlSi_mlffflag1.inpt index 0c506738..d5e36bd5 100644 --- a/tests/mlff/AlSi_mlffflag1/high_accuracy/AlSi_mlffflag1.inpt +++ b/tests/mlff/AlSi_mlffflag1/high_accuracy/AlSi_mlffflag1.inpt @@ -28,7 +28,7 @@ MLFF_RADIAL_BASIS: 8 # N_R (Number of raddial basis functions in SOAP) MLFF_RCUT_SOAP: 10 # (Cutoff radius for the SOAP descriptor) MLFF_MODEL_FOLDER: ./ -MLFF_INITIAL_STEPS_TRAIN: 3 # Initial number of MD steps where DFT calculation is imposed +MLFF_INITIAL_STEPS_TRAIN: 4 # Initial number of MD steps where DFT calculation is imposed MLFF_MAX_STR_STORE: 500 # This is the maximum number of structure for which memory is allocated to store MLFF data MLFF_MAX_CONFIG_STORE: 5000 # THis is the maximum number of atomic descriptors per element type that is stored MLFF_SCALE_FORCE: 1.0 # Scaling of force relative to energy diff --git a/tests/mlff/AlSi_mlffflag1/standard/AlSi_mlffflag1.inpt b/tests/mlff/AlSi_mlffflag1/standard/AlSi_mlffflag1.inpt index eabbe311..9d77f10e 100644 --- a/tests/mlff/AlSi_mlffflag1/standard/AlSi_mlffflag1.inpt +++ b/tests/mlff/AlSi_mlffflag1/standard/AlSi_mlffflag1.inpt @@ -27,7 +27,7 @@ MLFF_RADIAL_BASIS: 8 # N_R (Number of raddial basis functions in SOAP) MLFF_RCUT_SOAP: 10 # (Cutoff radius for the SOAP descriptor) MLFF_MODEL_FOLDER: ./ -MLFF_INITIAL_STEPS_TRAIN: 3 # Initial number of MD steps where DFT calculation is imposed +MLFF_INITIAL_STEPS_TRAIN: 4 # Initial number of MD steps where DFT calculation is imposed MLFF_MAX_STR_STORE: 500 # This is the maximum number of structure for which memory is allocated to store MLFF data MLFF_MAX_CONFIG_STORE: 5000 # THis is the maximum number of atomic descriptors per element type that is stored MLFF_SCALE_FORCE: 1.0 # Scaling of force relative to energy diff --git a/tests/mlff/C_sqmlff/high_accuracy/C_sqmlff.inpt b/tests/mlff/C_sqmlff/high_accuracy/C_sqmlff.inpt index 46f7757d..aa7eaefc 100644 --- a/tests/mlff/C_sqmlff/high_accuracy/C_sqmlff.inpt +++ b/tests/mlff/C_sqmlff/high_accuracy/C_sqmlff.inpt @@ -56,13 +56,12 @@ MLFF_FLAG: 1 MLFF_ANGULAR_BASIS: 6 # Lmax (maximum order of spherical harmonics in SOAP) MLFF_RADIAL_BASIS: 8 # N_R (Number of raddial basis functions in SOAP) MLFF_RCUT_SOAP: 10 # (Cutoff radius for the SOAP descriptor) -MLFF_hnl_FILE: ../hnl_N_8.txt # (name of hnl file which has h_nl and d h_nl/ dr) - (Lmax, N_R and RCUT_SOAP has to be consistent with this file) MLFF_MODEL_FOLDER: ./ -MLFF_INITIAL_STEPS_TRAIN: 8 # Initial number of MD steps where DFT calculation is imposed +MLFF_INITIAL_STEPS_TRAIN: 9 # Initial number of MD steps where DFT calculation is imposed MLFF_MAX_STR_STORE: 500 # This is the maximum number of structure for which memory is allocated to store MLFF data MLFF_MAX_CONFIG_STORE: 5000 # THis is the maximum number of atomic descriptors per element type that is stored MLFF_SCALE_FORCE: 1.0 # Scaling of force relative to energy MLFF_SCALE_STRESS: 1.0 1.0 1.0 1.0 1.0 1.0 # Scaling of stress relative to energy MLFF_REGUL_MIN: 1.00E-10 # regularization hyperparameters are chosen such that condition number of matrix to be inverted is less than 1/REGUL_MIN_MLFF MLFF_FACTOR_MULTIPLY_SIGMATOL: 1.001 # Parameter to control sigma_tol (always use >1) closer to will result in more DFT calulations -MLFF_IF_SPARSIFY_BEFORE_TRAIN: 1 # choose 0 if there is no sparsification just before training \ No newline at end of file +MLFF_IF_SPARSIFY_BEFORE_TRAIN: 1 # choose 0 if there is no sparsification just before training diff --git a/tests/mlff/C_sqmlff/standard/C_sqmlff.inpt b/tests/mlff/C_sqmlff/standard/C_sqmlff.inpt index 253b8f96..a031789a 100644 --- a/tests/mlff/C_sqmlff/standard/C_sqmlff.inpt +++ b/tests/mlff/C_sqmlff/standard/C_sqmlff.inpt @@ -55,9 +55,8 @@ MLFF_FLAG: 1 MLFF_ANGULAR_BASIS: 6 # Lmax (maximum order of spherical harmonics in SOAP) MLFF_RADIAL_BASIS: 8 # N_R (Number of raddial basis functions in SOAP) MLFF_RCUT_SOAP: 10 # (Cutoff radius for the SOAP descriptor) -MLFF_hnl_FILE: ../hnl_N_8.txt # (name of hnl file which has h_nl and d h_nl/ dr) - (Lmax, N_R and RCUT_SOAP has to be consistent with this file) MLFF_MODEL_FOLDER: ./ -MLFF_INITIAL_STEPS_TRAIN: 8 # Initial number of MD steps where DFT calculation is imposed +MLFF_INITIAL_STEPS_TRAIN: 9 # Initial number of MD steps where DFT calculation is imposed MLFF_MAX_STR_STORE: 500 # This is the maximum number of structure for which memory is allocated to store MLFF data MLFF_MAX_CONFIG_STORE: 5000 # THis is the maximum number of atomic descriptors per element type that is stored MLFF_SCALE_FORCE: 1.0 # Scaling of force relative to energy diff --git a/tests/samplescript_cluster b/tests/samplescript_cluster index c854eab7..a78c8c9c 100644 --- a/tests/samplescript_cluster +++ b/tests/samplescript_cluster @@ -1,16 +1,14 @@ #!/bin/bash #SBATCH -J SPARC_testsuite -#SBATCH -A gts-phanish6 # Account name -#SBATCH -p inferno # Partition +#SBATCH -A gts-amedford6 # Account name +#SBATCH -p inferno #SBATCH -N 2 --ntasks-per-node=24 #SBATCH --mem-per-cpu=7G #SBATCH -t1:00:00 cd $SLURM_SUBMIT_DIR module purge -module load gcc -module load mvapich2 -module load intel-oneapi-mkl/2023.1.0 +ml intel-oneapi-compilers intel-oneapi-mkl intel-oneapi-mpi echo $PWD