diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp index c79c58fe98..0aca805f7d 100644 --- a/src/USER-REAXC/pair_reax_c.cpp +++ b/src/USER-REAXC/pair_reax_c.cpp @@ -777,7 +777,7 @@ double PairReaxC::memory_usage() // From reaxc_allocate: BO bytes += 1.0 * system->total_cap * sizeof(reax_atom); - bytes += 19.0 * system->total_cap * sizeof(real); + bytes += 19.0 * system->total_cap * sizeof(double); bytes += 3.0 * system->total_cap * sizeof(int); // From reaxc_lists diff --git a/src/USER-REAXC/reaxc_allocate.cpp b/src/USER-REAXC/reaxc_allocate.cpp index 474b60d060..dc8545e006 100644 --- a/src/USER-REAXC/reaxc_allocate.cpp +++ b/src/USER-REAXC/reaxc_allocate.cpp @@ -191,14 +191,14 @@ int Allocate_Workspace( reax_system *system, control_params *control, int i, total_real, total_rvec, local_rvec; workspace->allocated = 1; - total_real = total_cap * sizeof(real); + total_real = total_cap * sizeof(double); total_rvec = total_cap * sizeof(rvec); local_rvec = local_cap * sizeof(rvec); /* communication storage */ for( i = 0; i < MAX_NBRS; ++i ) { - workspace->tmp_dbl[i] = (real*) - scalloc( total_cap, sizeof(real), "tmp_dbl", comm ); + workspace->tmp_dbl[i] = (double*) + scalloc( total_cap, sizeof(double), "tmp_dbl", comm ); workspace->tmp_rvec[i] = (rvec*) scalloc( total_cap, sizeof(rvec), "tmp_rvec", comm ); workspace->tmp_rvec2[i] = (rvec2*) @@ -208,62 +208,62 @@ int Allocate_Workspace( reax_system *system, control_params *control, /* bond order related storage */ workspace->within_bond_box = (int*) scalloc( total_cap, sizeof(int), "skin", comm ); - workspace->total_bond_order = (real*) smalloc( total_real, "total_bo", comm ); - workspace->Deltap = (real*) smalloc( total_real, "Deltap", comm ); - workspace->Deltap_boc = (real*) smalloc( total_real, "Deltap_boc", comm ); + workspace->total_bond_order = (double*) smalloc( total_real, "total_bo", comm ); + workspace->Deltap = (double*) smalloc( total_real, "Deltap", comm ); + workspace->Deltap_boc = (double*) smalloc( total_real, "Deltap_boc", comm ); workspace->dDeltap_self = (rvec*) smalloc( total_rvec, "dDeltap_self", comm ); - workspace->Delta = (real*) smalloc( total_real, "Delta", comm ); - workspace->Delta_lp = (real*) smalloc( total_real, "Delta_lp", comm ); - workspace->Delta_lp_temp = (real*) + workspace->Delta = (double*) smalloc( total_real, "Delta", comm ); + workspace->Delta_lp = (double*) smalloc( total_real, "Delta_lp", comm ); + workspace->Delta_lp_temp = (double*) smalloc( total_real, "Delta_lp_temp", comm ); - workspace->dDelta_lp = (real*) smalloc( total_real, "dDelta_lp", comm ); - workspace->dDelta_lp_temp = (real*) + workspace->dDelta_lp = (double*) smalloc( total_real, "dDelta_lp", comm ); + workspace->dDelta_lp_temp = (double*) smalloc( total_real, "dDelta_lp_temp", comm ); - workspace->Delta_e = (real*) smalloc( total_real, "Delta_e", comm ); - workspace->Delta_boc = (real*) smalloc( total_real, "Delta_boc", comm ); - workspace->Delta_val = (real*) smalloc( total_real, "Delta_val", comm ); - workspace->nlp = (real*) smalloc( total_real, "nlp", comm ); - workspace->nlp_temp = (real*) smalloc( total_real, "nlp_temp", comm ); - workspace->Clp = (real*) smalloc( total_real, "Clp", comm ); - workspace->vlpex = (real*) smalloc( total_real, "vlpex", comm ); + workspace->Delta_e = (double*) smalloc( total_real, "Delta_e", comm ); + workspace->Delta_boc = (double*) smalloc( total_real, "Delta_boc", comm ); + workspace->Delta_val = (double*) smalloc( total_real, "Delta_val", comm ); + workspace->nlp = (double*) smalloc( total_real, "nlp", comm ); + workspace->nlp_temp = (double*) smalloc( total_real, "nlp_temp", comm ); + workspace->Clp = (double*) smalloc( total_real, "Clp", comm ); + workspace->vlpex = (double*) smalloc( total_real, "vlpex", comm ); workspace->bond_mark = (int*) scalloc( total_cap, sizeof(int), "bond_mark", comm ); workspace->done_after = (int*) scalloc( total_cap, sizeof(int), "done_after", comm ); /* QEq storage */ - workspace->Hdia_inv = (real*) - scalloc( total_cap, sizeof(real), "Hdia_inv", comm ); - workspace->b_s = (real*) scalloc( total_cap, sizeof(real), "b_s", comm ); - workspace->b_t = (real*) scalloc( total_cap, sizeof(real), "b_t", comm ); - workspace->b_prc = (real*) scalloc( total_cap, sizeof(real), "b_prc", comm ); - workspace->b_prm = (real*) scalloc( total_cap, sizeof(real), "b_prm", comm ); - workspace->s = (real*) scalloc( total_cap, sizeof(real), "s", comm ); - workspace->t = (real*) scalloc( total_cap, sizeof(real), "t", comm ); - workspace->droptol = (real*) - scalloc( total_cap, sizeof(real), "droptol", comm ); + workspace->Hdia_inv = (double*) + scalloc( total_cap, sizeof(double), "Hdia_inv", comm ); + workspace->b_s = (double*) scalloc( total_cap, sizeof(double), "b_s", comm ); + workspace->b_t = (double*) scalloc( total_cap, sizeof(double), "b_t", comm ); + workspace->b_prc = (double*) scalloc( total_cap, sizeof(double), "b_prc", comm ); + workspace->b_prm = (double*) scalloc( total_cap, sizeof(double), "b_prm", comm ); + workspace->s = (double*) scalloc( total_cap, sizeof(double), "s", comm ); + workspace->t = (double*) scalloc( total_cap, sizeof(double), "t", comm ); + workspace->droptol = (double*) + scalloc( total_cap, sizeof(double), "droptol", comm ); workspace->b = (rvec2*) scalloc( total_cap, sizeof(rvec2), "b", comm ); workspace->x = (rvec2*) scalloc( total_cap, sizeof(rvec2), "x", comm ); /* GMRES storage */ - workspace->y = (real*) scalloc( RESTART+1, sizeof(real), "y", comm ); - workspace->z = (real*) scalloc( RESTART+1, sizeof(real), "z", comm ); - workspace->g = (real*) scalloc( RESTART+1, sizeof(real), "g", comm ); - workspace->h = (real**) scalloc( RESTART+1, sizeof(real*), "h", comm ); - workspace->hs = (real*) scalloc( RESTART+1, sizeof(real), "hs", comm ); - workspace->hc = (real*) scalloc( RESTART+1, sizeof(real), "hc", comm ); - workspace->v = (real**) scalloc( RESTART+1, sizeof(real*), "v", comm ); + workspace->y = (double*) scalloc( RESTART+1, sizeof(double), "y", comm ); + workspace->z = (double*) scalloc( RESTART+1, sizeof(double), "z", comm ); + workspace->g = (double*) scalloc( RESTART+1, sizeof(double), "g", comm ); + workspace->h = (double**) scalloc( RESTART+1, sizeof(double*), "h", comm ); + workspace->hs = (double*) scalloc( RESTART+1, sizeof(double), "hs", comm ); + workspace->hc = (double*) scalloc( RESTART+1, sizeof(double), "hc", comm ); + workspace->v = (double**) scalloc( RESTART+1, sizeof(double*), "v", comm ); for( i = 0; i < RESTART+1; ++i ) { - workspace->h[i] = (real*) scalloc( RESTART+1, sizeof(real), "h[i]", comm ); - workspace->v[i] = (real*) scalloc( total_cap, sizeof(real), "v[i]", comm ); + workspace->h[i] = (double*) scalloc( RESTART+1, sizeof(double), "h[i]", comm ); + workspace->v[i] = (double*) scalloc( total_cap, sizeof(double), "v[i]", comm ); } /* CG storage */ - workspace->r = (real*) scalloc( total_cap, sizeof(real), "r", comm ); - workspace->d = (real*) scalloc( total_cap, sizeof(real), "d", comm ); - workspace->q = (real*) scalloc( total_cap, sizeof(real), "q", comm ); - workspace->p = (real*) scalloc( total_cap, sizeof(real), "p", comm ); + workspace->r = (double*) scalloc( total_cap, sizeof(double), "r", comm ); + workspace->d = (double*) scalloc( total_cap, sizeof(double), "d", comm ); + workspace->q = (double*) scalloc( total_cap, sizeof(double), "q", comm ); + workspace->p = (double*) scalloc( total_cap, sizeof(double), "p", comm ); workspace->r2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "r2", comm ); workspace->d2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "d2", comm ); workspace->q2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "q2", comm ); @@ -274,8 +274,8 @@ int Allocate_Workspace( reax_system *system, control_params *control, // /* force related storage */ workspace->f = (rvec*) scalloc( total_cap, sizeof(rvec), "f", comm ); - workspace->CdDelta = (real*) - scalloc( total_cap, sizeof(real), "CdDelta", comm ); + workspace->CdDelta = (double*) + scalloc( total_cap, sizeof(double), "CdDelta", comm ); return SUCCESS; } diff --git a/src/USER-REAXC/reaxc_bond_orders.cpp b/src/USER-REAXC/reaxc_bond_orders.cpp index 7103f48299..e524c75e87 100644 --- a/src/USER-REAXC/reaxc_bond_orders.cpp +++ b/src/USER-REAXC/reaxc_bond_orders.cpp @@ -257,14 +257,14 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj, } -int BOp( storage *workspace, reax_list *bonds, real bo_cut, +int BOp( storage *workspace, reax_list *bonds, double bo_cut, int i, int btop_i, far_neighbor_data *nbr_pj, single_body_parameters *sbp_i, single_body_parameters *sbp_j, two_body_parameters *twbp ) { int j, btop_j; - real r2, C12, C34, C56; - real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2; - real BO, BO_s, BO_pi, BO_pi2; + double r2, C12, C34, C56; + double Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2; + double BO, BO_s, BO_pi, BO_pi2; bond_data *ibond, *jbond; bond_order_data *bo_ij, *bo_ji; @@ -370,14 +370,14 @@ void BO( reax_system *system, control_params *control, simulation_data *data, { int i, j, pj, type_i, type_j; int start_i, end_i, sym_index; - real val_i, Deltap_i, Deltap_boc_i; - real val_j, Deltap_j, Deltap_boc_j; - real f1, f2, f3, f4, f5, f4f5, exp_f4, exp_f5; - real exp_p1i, exp_p2i, exp_p1j, exp_p2j; - real temp, u1_ij, u1_ji, Cf1A_ij, Cf1B_ij, Cf1_ij, Cf1_ji; - real Cf45_ij, Cf45_ji, p_lp1; //u_ij, u_ji - real A0_ij, A1_ij, A2_ij, A2_ji, A3_ij, A3_ji; - real explp1, p_boc1, p_boc2; + double val_i, Deltap_i, Deltap_boc_i; + double val_j, Deltap_j, Deltap_boc_j; + double f1, f2, f3, f4, f5, f4f5, exp_f4, exp_f5; + double exp_p1i, exp_p2i, exp_p1j, exp_p2j; + double temp, u1_ij, u1_ji, Cf1A_ij, Cf1B_ij, Cf1_ij, Cf1_ji; + double Cf45_ij, Cf45_ji, p_lp1; //u_ij, u_ji + double A0_ij, A1_ij, A2_ij, A2_ji, A3_ij, A3_ji; + double explp1, p_boc1, p_boc2; single_body_parameters *sbp_i, *sbp_j; two_body_parameters *twbp; bond_order_data *bo_ij, *bo_ji; diff --git a/src/USER-REAXC/reaxc_bond_orders.h b/src/USER-REAXC/reaxc_bond_orders.h index 0912a936d6..8f1e182aa2 100644 --- a/src/USER-REAXC/reaxc_bond_orders.h +++ b/src/USER-REAXC/reaxc_bond_orders.h @@ -30,33 +30,33 @@ #include "reaxc_types.h" typedef struct{ - real C1dbo, C2dbo, C3dbo; - real C1dbopi, C2dbopi, C3dbopi, C4dbopi; - real C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2; - real C1dDelta, C2dDelta, C3dDelta; + double C1dbo, C2dbo, C3dbo; + double C1dbopi, C2dbopi, C3dbopi, C4dbopi; + double C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2; + double C1dDelta, C2dDelta, C3dDelta; }dbond_coefficients; #ifdef TEST_FORCES -void Get_dBO( reax_system*, reax_list**, int, int, real, rvec* ); +void Get_dBO( reax_system*, reax_list**, int, int, double, rvec* ); void Get_dBOpinpi2( reax_system*, reax_list**, - int, int, real, real, rvec*, rvec* ); + int, int, double, double, rvec*, rvec* ); -void Add_dBO( reax_system*, reax_list**, int, int, real, rvec* ); +void Add_dBO( reax_system*, reax_list**, int, int, double, rvec* ); void Add_dBOpinpi2( reax_system*, reax_list**, - int, int, real, real, rvec*, rvec* ); + int, int, double, double, rvec*, rvec* ); -void Add_dBO_to_Forces( reax_system*, reax_list**, int, int, real ); +void Add_dBO_to_Forces( reax_system*, reax_list**, int, int, double ); void Add_dBOpinpi2_to_Forces( reax_system*, reax_list**, - int, int, real, real ); + int, int, double, double ); -void Add_dDelta( reax_system*, reax_list**, int, real, rvec* ); -void Add_dDelta_to_Forces( reax_system *, reax_list**, int, real ); +void Add_dDelta( reax_system*, reax_list**, int, double, rvec* ); +void Add_dDelta_to_Forces( reax_system *, reax_list**, int, double ); #endif void Add_dBond_to_Forces( reax_system*, int, int, storage*, reax_list** ); void Add_dBond_to_Forces_NPT( int, int, simulation_data*, storage*, reax_list** ); -int BOp(storage*, reax_list*, real, int, int, far_neighbor_data*, +int BOp(storage*, reax_list*, double, int, int, far_neighbor_data*, single_body_parameters*, single_body_parameters*, two_body_parameters*); void BO( reax_system*, control_params*, simulation_data*, storage*, reax_list**, output_controls* ); diff --git a/src/USER-REAXC/reaxc_bonds.cpp b/src/USER-REAXC/reaxc_bonds.cpp index 857ee40701..6e3fa75978 100644 --- a/src/USER-REAXC/reaxc_bonds.cpp +++ b/src/USER-REAXC/reaxc_bonds.cpp @@ -38,10 +38,10 @@ void Bonds( reax_system *system, control_params *control, int i, j, pj, natoms; int start_i, end_i; int type_i, type_j; - real ebond, pow_BOs_be2, exp_be12, CEbo; - real gp3, gp4, gp7, gp10, gp37; - real exphu, exphua1, exphub1, exphuov, hulpov, estriph; - real decobdbo, decobdboua, decobdboub; + double ebond, pow_BOs_be2, exp_be12, CEbo; + double gp3, gp4, gp7, gp10, gp37; + double exphu, exphua1, exphub1, exphuov, hulpov, estriph; + double decobdbo, decobdboua, decobdboub; single_body_parameters *sbp_i, *sbp_j; two_body_parameters *twbp; bond_order_data *bo_ij; diff --git a/src/USER-REAXC/reaxc_control.cpp b/src/USER-REAXC/reaxc_control.cpp index 9d067c16a5..3753360c68 100644 --- a/src/USER-REAXC/reaxc_control.cpp +++ b/src/USER-REAXC/reaxc_control.cpp @@ -34,7 +34,7 @@ char Read_Control_File( char *control_file, control_params* control, FILE *fp; char *s, **tmp; int i,ival; - real val; + double val; /* open control file */ if ( (fp = fopen( control_file, "r" ) ) == NULL ) { diff --git a/src/USER-REAXC/reaxc_ffield.cpp b/src/USER-REAXC/reaxc_ffield.cpp index 0292e76320..74791d47fd 100644 --- a/src/USER-REAXC/reaxc_ffield.cpp +++ b/src/USER-REAXC/reaxc_ffield.cpp @@ -38,7 +38,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, int c, i, j, k, l, m, n, o, p, cnt; int lgflag = control->lgflag; int errorflag = 1; - real val; + double val; MPI_Comm comm; comm = MPI_COMM_WORLD; @@ -64,14 +64,14 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax, } reax->gp.n_global = n; - reax->gp.l = (real*) malloc(sizeof(real)*n); + reax->gp.l = (double*) malloc(sizeof(double)*n); /* see reax_types.h for mapping between l[i] and the lambdas used in ff */ for (i=0; i < n; i++) { fgets(s,MAX_LINE,fp); c = Tokenize(s,&tmp); - val = (real) atof(tmp[0]); + val = (double) atof(tmp[0]); reax->gp.l[i] = val; } diff --git a/src/USER-REAXC/reaxc_forces.cpp b/src/USER-REAXC/reaxc_forces.cpp index b9fd78dd76..84bb86e5d2 100644 --- a/src/USER-REAXC/reaxc_forces.cpp +++ b/src/USER-REAXC/reaxc_forces.cpp @@ -171,9 +171,9 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists, } -real Compute_H( real r, real gamma, real *ctap ) +double Compute_H( double r, double gamma, double *ctap ) { - real taper, dr3gamij_1, dr3gamij_3; + double taper, dr3gamij_1, dr3gamij_3; taper = ctap[7] * r + ctap[6]; taper = taper * r + ctap[5]; @@ -189,10 +189,10 @@ real Compute_H( real r, real gamma, real *ctap ) } -real Compute_tabH( real r_ij, int ti, int tj ) +double Compute_tabH( double r_ij, int ti, int tj ) { int r, tmin, tmax; - real val, dif, base; + double val, dif, base; LR_lookup_table *t; tmin = MIN( ti, tj ); @@ -202,7 +202,7 @@ real Compute_tabH( real r_ij, int ti, int tj ) /* cubic spline interpolation */ r = (int)(r_ij * t->inv_dx); if( r == 0 ) ++r; - base = (real)(r+1) * t->dx; + base = (double)(r+1) * t->dx; dif = r_ij - base; val = ((t->ele[r].d*dif + t->ele[r].c)*dif + t->ele[r].b)*dif + t->ele[r].a; @@ -221,7 +221,7 @@ void Init_Forces( reax_system *system, control_params *control, int Htop, btop_i, btop_j, num_bonds, num_hbonds; int ihb, jhb, ihb_top, jhb_top; int local, flag, renbr; - real r_ij, cutoff; + double r_ij, cutoff; sparse_matrix *H; reax_list *far_nbrs, *bonds, *hbonds; single_body_parameters *sbp_i, *sbp_j; @@ -388,7 +388,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control, int btop_i, btop_j, num_bonds, num_hbonds; int ihb, jhb, ihb_top, jhb_top; int local, flag, renbr; - real cutoff; + double cutoff; reax_list *far_nbrs, *bonds, *hbonds; single_body_parameters *sbp_i, *sbp_j; two_body_parameters *twbp; @@ -531,10 +531,10 @@ void Estimate_Storages( reax_system *system, control_params *control, int type_i, type_j; int ihb, jhb; int local; - real cutoff; - real r_ij; - real C12, C34, C56; - real BO, BO_s, BO_pi, BO_pi2; + double cutoff; + double r_ij; + double C12, C34, C56; + double BO, BO_s, BO_pi, BO_pi2; reax_list *far_nbrs; single_body_parameters *sbp_i, *sbp_j; two_body_parameters *twbp; diff --git a/src/USER-REAXC/reaxc_hydrogen_bonds.cpp b/src/USER-REAXC/reaxc_hydrogen_bonds.cpp index 7cbd0bb0a9..f58f9c5842 100644 --- a/src/USER-REAXC/reaxc_hydrogen_bonds.cpp +++ b/src/USER-REAXC/reaxc_hydrogen_bonds.cpp @@ -42,8 +42,8 @@ void Hydrogen_Bonds( reax_system *system, control_params *control, int itr, top; int num_hb_intrs = 0; ivec rel_jk; - real r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; - real e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3; + double r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; + double e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3; rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk; rvec dvec_jk, force, ext_press; hbond_parameters *hbp; @@ -55,7 +55,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control, hbond_data *hbond_list; // tally variables - real fi_tmp[3], fk_tmp[3], delij[3], delkj[3]; + double fi_tmp[3], fk_tmp[3], delij[3], delkj[3]; bonds = (*lists) + BONDS; bond_list = bonds->select.bond_list; diff --git a/src/USER-REAXC/reaxc_init_md.cpp b/src/USER-REAXC/reaxc_init_md.cpp index 505f3b5a33..f912c95ea5 100644 --- a/src/USER-REAXC/reaxc_init_md.cpp +++ b/src/USER-REAXC/reaxc_init_md.cpp @@ -82,9 +82,9 @@ int Init_Simulation_Data( reax_system *system, control_params *control, void Init_Taper( control_params *control, storage *workspace, MPI_Comm comm ) { - real d1, d7; - real swa, swa2, swa3; - real swb, swb2, swb3; + double d1, d7; + double swa, swa2, swa3; + double swb, swb2, swb3; swa = control->nonb_low; swb = control->nonb_cut; diff --git a/src/USER-REAXC/reaxc_io_tools.cpp b/src/USER-REAXC/reaxc_io_tools.cpp index 513b59c9f6..1a978d0073 100644 --- a/src/USER-REAXC/reaxc_io_tools.cpp +++ b/src/USER-REAXC/reaxc_io_tools.cpp @@ -430,7 +430,7 @@ void Print_Linear_System( reax_system *system, control_params *control, } -void Print_LinSys_Soln( reax_system *system, real *x, real *b_prm, real *b ) +void Print_LinSys_Soln( reax_system *system, double *x, double *b_prm, double *b ) { int i; char fname[100]; diff --git a/src/USER-REAXC/reaxc_io_tools.h b/src/USER-REAXC/reaxc_io_tools.h index 8d25a85d12..12312d344e 100644 --- a/src/USER-REAXC/reaxc_io_tools.h +++ b/src/USER-REAXC/reaxc_io_tools.h @@ -49,7 +49,7 @@ void Print_Far_Neighbors( reax_system*, reax_list**, control_params *); void Print_Sparse_Matrix( reax_system*, sparse_matrix* ); void Print_Sparse_Matrix2( reax_system*, sparse_matrix*, char* ); void Print_Linear_System( reax_system*, control_params*, storage*, int ); -void Print_LinSys_Soln( reax_system*, real*, real*, real* ); +void Print_LinSys_Soln( reax_system*, double*, double*, double* ); void Print_Charges( reax_system* ); void Print_Bonds( reax_system*, reax_list*, char* ); void Print_Bond_List2( reax_system*, reax_list*, char* ); diff --git a/src/USER-REAXC/reaxc_lookup.cpp b/src/USER-REAXC/reaxc_lookup.cpp index fe01d4bb46..7497a5c096 100644 --- a/src/USER-REAXC/reaxc_lookup.cpp +++ b/src/USER-REAXC/reaxc_lookup.cpp @@ -31,10 +31,10 @@ LR_lookup_table **LR; -void Tridiagonal_Solve( const real *a, const real *b, - real *c, real *d, real *x, unsigned int n){ +void Tridiagonal_Solve( const double *a, const double *b, + double *c, double *d, double *x, unsigned int n){ int i; - real id; + double id; c[0] /= b[0]; /* Division by zero risk. */ d[0] /= b[0]; /* Division by zero would imply a singular matrix. */ @@ -50,19 +50,19 @@ void Tridiagonal_Solve( const real *a, const real *b, } -void Natural_Cubic_Spline( const real *h, const real *f, +void Natural_Cubic_Spline( const double *h, const double *f, cubic_spline_coef *coef, unsigned int n, MPI_Comm comm ) { int i; - real *a, *b, *c, *d, *v; + double *a, *b, *c, *d, *v; /* allocate space for the linear system */ - a = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); - b = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); - c = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); - d = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); - v = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); + a = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); + b = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); + c = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); + d = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); + v = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); /* build the linear system */ a[0] = a[1] = a[n-1] = 0; @@ -101,19 +101,19 @@ void Natural_Cubic_Spline( const real *h, const real *f, -void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast, +void Complete_Cubic_Spline( const double *h, const double *f, double v0, double vlast, cubic_spline_coef *coef, unsigned int n, MPI_Comm comm ) { int i; - real *a, *b, *c, *d, *v; + double *a, *b, *c, *d, *v; /* allocate space for the linear system */ - a = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); - b = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); - c = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); - d = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); - v = (real*) smalloc( n * sizeof(real), "cubic_spline:a", comm ); + a = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); + b = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); + c = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); + d = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); + v = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm ); /* build the linear system */ a[0] = 0; @@ -150,14 +150,14 @@ void Complete_Cubic_Spline( const real *h, const real *f, real v0, real vlast, } -void LR_Lookup( LR_lookup_table *t, real r, LR_data *y ) +void LR_Lookup( LR_lookup_table *t, double r, LR_data *y ) { int i; - real base, dif; + double base, dif; i = (int)(r * t->inv_dx); if( i == 0 ) ++i; - base = (real)(i+1) * t->dx; + base = (double)(i+1) * t->dx; dif = r - base; y->e_vdW = ((t->vdW[i].d*dif + t->vdW[i].c)*dif + t->vdW[i].b)*dif + @@ -180,9 +180,9 @@ int Init_Lookup_Tables( reax_system *system, control_params *control, int i, j, r; int num_atom_types; int existing_types[MAX_ATOM_TYPES], aggregated[MAX_ATOM_TYPES]; - real dr; - real *h, *fh, *fvdw, *fele, *fCEvd, *fCEclmb; - real v0_vdw, v0_ele, vlast_vdw, vlast_ele; + double dr; + double *h, *fh, *fvdw, *fele, *fCEvd, *fCEclmb; + double v0_vdw, v0_ele, vlast_vdw, vlast_ele; MPI_Comm comm; /* initializations */ @@ -194,18 +194,18 @@ int Init_Lookup_Tables( reax_system *system, control_params *control, num_atom_types = system->reax_param.num_atom_types; dr = control->nonb_cut / control->tabulate; - h = (real*) - smalloc( (control->tabulate+2) * sizeof(real), "lookup:h", comm ); - fh = (real*) - smalloc( (control->tabulate+2) * sizeof(real), "lookup:fh", comm ); - fvdw = (real*) - smalloc( (control->tabulate+2) * sizeof(real), "lookup:fvdw", comm ); - fCEvd = (real*) - smalloc( (control->tabulate+2) * sizeof(real), "lookup:fCEvd", comm ); - fele = (real*) - smalloc( (control->tabulate+2) * sizeof(real), "lookup:fele", comm ); - fCEclmb = (real*) - smalloc( (control->tabulate+2) * sizeof(real), "lookup:fCEclmb", comm ); + h = (double*) + smalloc( (control->tabulate+2) * sizeof(double), "lookup:h", comm ); + fh = (double*) + smalloc( (control->tabulate+2) * sizeof(double), "lookup:fh", comm ); + fvdw = (double*) + smalloc( (control->tabulate+2) * sizeof(double), "lookup:fvdw", comm ); + fCEvd = (double*) + smalloc( (control->tabulate+2) * sizeof(double), "lookup:fCEvd", comm ); + fele = (double*) + smalloc( (control->tabulate+2) * sizeof(double), "lookup:fele", comm ); + fCEclmb = (double*) + smalloc( (control->tabulate+2) * sizeof(double), "lookup:fCEclmb", comm ); LR = (LR_lookup_table**) scalloc( num_atom_types, sizeof(LR_lookup_table*), "lookup:LR", comm ); diff --git a/src/USER-REAXC/reaxc_multi_body.cpp b/src/USER-REAXC/reaxc_multi_body.cpp index f1b56ef5de..1923668e89 100644 --- a/src/USER-REAXC/reaxc_multi_body.cpp +++ b/src/USER-REAXC/reaxc_multi_body.cpp @@ -35,17 +35,17 @@ void Atom_Energy( reax_system *system, control_params *control, output_controls *out_control ) { int i, j, pj, type_i, type_j; - real Delta_lpcorr, dfvl; - real e_lp, expvd2, inv_expvd2, dElp, CElp, DlpVi; - real e_lph, Di, vov3, deahu2dbo, deahu2dsbo; - real e_ov, CEover1, CEover2, CEover3, CEover4; - real exp_ovun1, exp_ovun2, sum_ovun1, sum_ovun2; - real exp_ovun2n, exp_ovun6, exp_ovun8; - real inv_exp_ovun1, inv_exp_ovun2, inv_exp_ovun2n, inv_exp_ovun8; - real e_un, CEunder1, CEunder2, CEunder3, CEunder4; - real p_lp2, p_lp3; - real p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8; - real eng_tmp; + double Delta_lpcorr, dfvl; + double e_lp, expvd2, inv_expvd2, dElp, CElp, DlpVi; + double e_lph, Di, vov3, deahu2dbo, deahu2dsbo; + double e_ov, CEover1, CEover2, CEover3, CEover4; + double exp_ovun1, exp_ovun2, sum_ovun1, sum_ovun2; + double exp_ovun2n, exp_ovun6, exp_ovun8; + double inv_exp_ovun1, inv_exp_ovun2, inv_exp_ovun2n, inv_exp_ovun8; + double e_un, CEunder1, CEunder2, CEunder3, CEunder4; + double p_lp2, p_lp3; + double p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8; + double eng_tmp; int numbonds; single_body_parameters *sbp_i; diff --git a/src/USER-REAXC/reaxc_nonbonded.cpp b/src/USER-REAXC/reaxc_nonbonded.cpp index 5cb9c9ee43..fff046ba7c 100644 --- a/src/USER-REAXC/reaxc_nonbonded.cpp +++ b/src/USER-REAXC/reaxc_nonbonded.cpp @@ -38,20 +38,20 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control, int i, j, pj, natoms; int start_i, end_i, flag; rc_tagint orig_i, orig_j; - real p_vdW1, p_vdW1i; - real powr_vdW1, powgi_vdW1; - real tmp, r_ij, fn13, exp1, exp2; - real Tap, dTap, dfn13, CEvd, CEclmb, de_core; - real dr3gamij_1, dr3gamij_3; - real e_ele, e_vdW, e_core, SMALL = 0.0001; - real e_lg, de_lg, r_ij5, r_ij6, re6; + double p_vdW1, p_vdW1i; + double powr_vdW1, powgi_vdW1; + double tmp, r_ij, fn13, exp1, exp2; + double Tap, dTap, dfn13, CEvd, CEclmb, de_core; + double dr3gamij_1, dr3gamij_3; + double e_ele, e_vdW, e_core, SMALL = 0.0001; + double e_lg, de_lg, r_ij5, r_ij6, re6; rvec temp, ext_press; two_body_parameters *twbp; far_neighbor_data *nbr_pj; reax_list *far_nbrs; // Tallying variables: - real pe_vdw, f_tmp, delij[3]; + double pe_vdw, f_tmp, delij[3]; natoms = system->n; far_nbrs = (*lists) + FAR_NBRS; @@ -212,10 +212,10 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control, int type_i, type_j, tmin, tmax; int start_i, end_i, flag; rc_tagint orig_i, orig_j; - real r_ij, base, dif; - real e_vdW, e_ele; - real CEvd, CEclmb, SMALL = 0.0001; - real f_tmp, delij[3]; + double r_ij, base, dif; + double e_vdW, e_ele; + double CEvd, CEclmb, SMALL = 0.0001; + double f_tmp, delij[3]; rvec temp, ext_press; far_neighbor_data *nbr_pj; @@ -265,7 +265,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control, /* Cubic Spline Interpolation */ r = (int)(r_ij * t->inv_dx); if( r == 0 ) ++r; - base = (real)(r+1) * t->dx; + base = (double)(r+1) * t->dx; dif = r_ij - base; e_vdW = ((t->vdW[r].d*dif + t->vdW[r].c)*dif + t->vdW[r].b)*dif + @@ -319,7 +319,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system,control_params *control, void Compute_Polarization_Energy( reax_system *system, simulation_data *data ) { int i, type_i; - real q, en_tmp; + double q, en_tmp; data->my_en.e_pol = 0.0; for( i = 0; i < system->n; i++ ) { @@ -338,16 +338,16 @@ void Compute_Polarization_Energy( reax_system *system, simulation_data *data ) } void LR_vdW_Coulomb( reax_system *system, storage *workspace, - control_params *control, int i, int j, real r_ij, LR_data *lr ) + control_params *control, int i, int j, double r_ij, LR_data *lr ) { - real p_vdW1 = system->reax_param.gp.l[28]; - real p_vdW1i = 1.0 / p_vdW1; - real powr_vdW1, powgi_vdW1; - real tmp, fn13, exp1, exp2; - real Tap, dTap, dfn13; - real dr3gamij_1, dr3gamij_3; - real e_core, de_core; - real e_lg, de_lg, r_ij5, r_ij6, re6; + double p_vdW1 = system->reax_param.gp.l[28]; + double p_vdW1i = 1.0 / p_vdW1; + double powr_vdW1, powgi_vdW1; + double tmp, fn13, exp1, exp2; + double Tap, dTap, dfn13; + double dr3gamij_1, dr3gamij_3; + double e_core, de_core; + double e_lg, de_lg, r_ij5, r_ij6, re6; two_body_parameters *twbp; twbp = &(system->reax_param.tbp[i][j]); diff --git a/src/USER-REAXC/reaxc_nonbonded.h b/src/USER-REAXC/reaxc_nonbonded.h index ecb07494c6..9a29d4d8c1 100644 --- a/src/USER-REAXC/reaxc_nonbonded.h +++ b/src/USER-REAXC/reaxc_nonbonded.h @@ -39,5 +39,5 @@ void Tabulated_vdW_Coulomb_Energy( reax_system*, control_params*, void Compute_Polarization_Energy( reax_system*, simulation_data* ); void LR_vdW_Coulomb( reax_system*, storage*, control_params*, - int, int, real, LR_data* ); + int, int, double, LR_data* ); #endif diff --git a/src/USER-REAXC/reaxc_reset_tools.cpp b/src/USER-REAXC/reaxc_reset_tools.cpp index 7c3ba779ab..7d9e9deb76 100644 --- a/src/USER-REAXC/reaxc_reset_tools.cpp +++ b/src/USER-REAXC/reaxc_reset_tools.cpp @@ -111,9 +111,9 @@ void Reset_Timing( reax_timing *rt ) void Reset_Workspace( reax_system *system, storage *workspace ) { - memset( workspace->total_bond_order, 0, system->total_cap * sizeof( real ) ); + memset( workspace->total_bond_order, 0, system->total_cap * sizeof( double ) ); memset( workspace->dDeltap_self, 0, system->total_cap * sizeof( rvec ) ); - memset( workspace->CdDelta, 0, system->total_cap * sizeof( real ) ); + memset( workspace->CdDelta, 0, system->total_cap * sizeof( double ) ); memset( workspace->f, 0, system->total_cap * sizeof( rvec ) ); } diff --git a/src/USER-REAXC/reaxc_system_props.cpp b/src/USER-REAXC/reaxc_system_props.cpp index f9abcc770a..554151ba82 100644 --- a/src/USER-REAXC/reaxc_system_props.cpp +++ b/src/USER-REAXC/reaxc_system_props.cpp @@ -31,7 +31,7 @@ void Temperature_Control( control_params *control, simulation_data *data ) { - real tmp; + double tmp; if( control->T_mode == 1 ) {// step-wise temperature control if((data->step-data->prev_steps) % ((int)(control->T_freq/control->dt))==0){ @@ -55,7 +55,7 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data, { int i; rvec p; - real m; + double m; data->my_en.e_kin = 0.0; data->sys_en.e_kin = 0.0; @@ -82,7 +82,7 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data, void Compute_System_Energy( reax_system *system, simulation_data *data, MPI_Comm comm ) { - real my_en[15], sys_en[15]; + double my_en[15], sys_en[15]; my_en[0] = data->my_en.e_bond; my_en[1] = data->my_en.e_ov; @@ -141,7 +141,7 @@ void Compute_Total_Mass( reax_system *system, simulation_data *data, MPI_Comm comm ) { int i; - real tmp; + double tmp; tmp = 0; for( i = 0; i < system->n; i++ ) @@ -158,8 +158,8 @@ void Compute_Center_of_Mass( reax_system *system, simulation_data *data, mpi_datatypes *mpi_data, MPI_Comm comm ) { int i; - real m, det; //xx, xy, xz, yy, yz, zz; - real tmp_mat[6], tot_mat[6]; + double m, det; //xx, xy, xz, yy, yz, zz; + double tmp_mat[6], tot_mat[6]; rvec my_xcm, my_vcm, my_amcm, my_avcm; rvec tvec, diff; rtensor mat, inv; diff --git a/src/USER-REAXC/reaxc_tool_box.cpp b/src/USER-REAXC/reaxc_tool_box.cpp index 84a654fecf..4096c75118 100644 --- a/src/USER-REAXC/reaxc_tool_box.cpp +++ b/src/USER-REAXC/reaxc_tool_box.cpp @@ -30,7 +30,7 @@ void Transform( rvec x1, simulation_box *box, char flag, rvec x2 ) { int i, j; - real tmp; + double tmp; if (flag > 0) { for (i=0; i < 3; i++) { @@ -79,16 +79,16 @@ void Fit_to_Periodic_Box( simulation_box *box, rvec *p ) } struct timeval tim; -real t_end; +double t_end; -real Get_Time( ) +double Get_Time( ) { gettimeofday(&tim, NULL ); return( tim.tv_sec + (tim.tv_usec / 1000000.0) ); } -real Get_Timing_Info( real t_start ) +double Get_Timing_Info( double t_start ) { gettimeofday(&tim, NULL ); t_end = tim.tv_sec + (tim.tv_usec / 1000000.0); @@ -96,7 +96,7 @@ real Get_Timing_Info( real t_start ) } -void Update_Timing_Info( real *t_start, real *timing ) +void Update_Timing_Info( double *t_start, double *timing ) { gettimeofday(&tim, NULL ); t_end = tim.tv_sec + (tim.tv_usec / 1000000.0); diff --git a/src/USER-REAXC/reaxc_tool_box.h b/src/USER-REAXC/reaxc_tool_box.h index e45b07d000..2a5ea1d6db 100644 --- a/src/USER-REAXC/reaxc_tool_box.h +++ b/src/USER-REAXC/reaxc_tool_box.h @@ -44,14 +44,14 @@ int iown_midpoint( simulation_box*, rvec, rvec ); /* from grid.h */ void GridCell_Closest_Point( grid_cell*, grid_cell*, ivec, ivec, rvec ); void GridCell_to_Box_Points( grid_cell*, ivec, rvec, rvec ); -real DistSqr_between_Special_Points( rvec, rvec ); -real DistSqr_to_Special_Point( rvec, rvec ); +double DistSqr_between_Special_Points( rvec, rvec ); +double DistSqr_to_Special_Point( rvec, rvec ); int Relative_Coord_Encoding( ivec ); /* from system_props.h */ -real Get_Time( ); -real Get_Timing_Info( real ); -void Update_Timing_Info( real*, real* ); +double Get_Time( ); +double Get_Timing_Info( double ); +void Update_Timing_Info( double*, double* ); /* from io_tools.h */ int Get_Atom_Type( reax_interaction*, char*, MPI_Comm ); diff --git a/src/USER-REAXC/reaxc_torsion_angles.cpp b/src/USER-REAXC/reaxc_torsion_angles.cpp index 7609c66022..bf73fc6ae8 100644 --- a/src/USER-REAXC/reaxc_torsion_angles.cpp +++ b/src/USER-REAXC/reaxc_torsion_angles.cpp @@ -33,20 +33,20 @@ #define MIN_SINE 1e-10 -real Calculate_Omega( rvec dvec_ij, real r_ij, - rvec dvec_jk, real r_jk, - rvec dvec_kl, real r_kl, - rvec dvec_li, real r_li, +double Calculate_Omega( rvec dvec_ij, double r_ij, + rvec dvec_jk, double r_jk, + rvec dvec_kl, double r_kl, + rvec dvec_li, double r_li, three_body_interaction_data *p_ijk, three_body_interaction_data *p_jkl, rvec dcos_omega_di, rvec dcos_omega_dj, rvec dcos_omega_dk, rvec dcos_omega_dl, output_controls *out_control ) { - real unnorm_cos_omega, unnorm_sin_omega, omega; - real sin_ijk, cos_ijk, sin_jkl, cos_jkl; - real htra, htrb, htrc, hthd, hthe, hnra, hnrc, hnhd, hnhe; - real arg, poem, tel; + double unnorm_cos_omega, unnorm_sin_omega, omega; + double sin_ijk, cos_ijk, sin_jkl, cos_jkl; + double htra, htrb, htrc, hthd, hthe, hnra, hnrc, hnhd, hnhe; + double arg, poem, tel; rvec cross_jk_kl; sin_ijk = sin( p_ijk->theta ); @@ -128,25 +128,25 @@ void Torsion_Angles( reax_system *system, control_params *control, int start_pj, end_pj, start_pk, end_pk; int num_frb_intrs = 0; - real Delta_j, Delta_k; - real r_ij, r_jk, r_kl, r_li; - real BOA_ij, BOA_jk, BOA_kl; + double Delta_j, Delta_k; + double r_ij, r_jk, r_kl, r_li; + double BOA_ij, BOA_jk, BOA_kl; - real exp_tor2_ij, exp_tor2_jk, exp_tor2_kl; - real exp_tor1, exp_tor3_DjDk, exp_tor4_DjDk, exp_tor34_inv; - real exp_cot2_jk, exp_cot2_ij, exp_cot2_kl; - real fn10, f11_DjDk, dfn11, fn12; - real theta_ijk, theta_jkl; - real sin_ijk, sin_jkl; - real cos_ijk, cos_jkl; - real tan_ijk_i, tan_jkl_i; - real omega, cos_omega, cos2omega, cos3omega; + double exp_tor2_ij, exp_tor2_jk, exp_tor2_kl; + double exp_tor1, exp_tor3_DjDk, exp_tor4_DjDk, exp_tor34_inv; + double exp_cot2_jk, exp_cot2_ij, exp_cot2_kl; + double fn10, f11_DjDk, dfn11, fn12; + double theta_ijk, theta_jkl; + double sin_ijk, sin_jkl; + double cos_ijk, cos_jkl; + double tan_ijk_i, tan_jkl_i; + double omega, cos_omega, cos2omega, cos3omega; rvec dcos_omega_di, dcos_omega_dj, dcos_omega_dk, dcos_omega_dl; - real CV, cmn, CEtors1, CEtors2, CEtors3, CEtors4; - real CEtors5, CEtors6, CEtors7, CEtors8, CEtors9; - real Cconj, CEconj1, CEconj2, CEconj3; - real CEconj4, CEconj5, CEconj6; - real e_tor, e_con; + double CV, cmn, CEtors1, CEtors2, CEtors3, CEtors4; + double CEtors5, CEtors6, CEtors7, CEtors8, CEtors9; + double Cconj, CEconj1, CEconj2, CEconj3; + double CEconj4, CEconj5, CEconj6; + double e_tor, e_con; rvec dvec_li; rvec force, ext_press; ivec rel_box_jl; @@ -156,16 +156,16 @@ void Torsion_Angles( reax_system *system, control_params *control, bond_data *pbond_ij, *pbond_jk, *pbond_kl; bond_order_data *bo_ij, *bo_jk, *bo_kl; three_body_interaction_data *p_ijk, *p_jkl; - real p_tor2 = system->reax_param.gp.l[23]; - real p_tor3 = system->reax_param.gp.l[24]; - real p_tor4 = system->reax_param.gp.l[25]; - real p_cot2 = system->reax_param.gp.l[27]; + double p_tor2 = system->reax_param.gp.l[23]; + double p_tor3 = system->reax_param.gp.l[24]; + double p_tor4 = system->reax_param.gp.l[25]; + double p_cot2 = system->reax_param.gp.l[27]; reax_list *bonds = (*lists) + BONDS; reax_list *thb_intrs = (*lists) + THREE_BODIES; // Virial tallying variables - real delil[3], deljl[3], delkl[3]; - real eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; + double delil[3], deljl[3], delkl[3]; + double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; natoms = system->n; diff --git a/src/USER-REAXC/reaxc_types.h b/src/USER-REAXC/reaxc_types.h index b05e86f16f..5d6d5d8024 100644 --- a/src/USER-REAXC/reaxc_types.h +++ b/src/USER-REAXC/reaxc_types.h @@ -62,11 +62,10 @@ /********************** TYPE DEFINITIONS ********************/ typedef int ivec[3]; -typedef double real; -typedef real rvec[3]; -typedef real rtensor[3][3]; -typedef real rvec2[2]; -typedef real rvec4[4]; +typedef double rvec[3]; +typedef double rtensor[3][3]; +typedef double rvec2[2]; +typedef double rvec4[4]; // import LAMMPS' definition of tagint @@ -111,7 +110,7 @@ typedef struct typedef struct { int n_global; - real* l; + double* l; int vdw_type; } global_parameters; @@ -122,44 +121,44 @@ typedef struct /* Line one in field file */ char name[15]; // Two character atom name - real r_s; - real valency; // Valency of the atom - real mass; // Mass of atom - real r_vdw; - real epsilon; - real gamma; - real r_pi; - real valency_e; - real nlp_opt; + double r_s; + double valency; // Valency of the atom + double mass; // Mass of atom + double r_vdw; + double epsilon; + double gamma; + double r_pi; + double valency_e; + double nlp_opt; /* Line two in field file */ - real alpha; - real gamma_w; - real valency_boc; - real p_ovun5; - real chi; - real eta; + double alpha; + double gamma_w; + double valency_boc; + double p_ovun5; + double chi; + double eta; int p_hbond; // 1 for H, 2 for hbonding atoms (O,S,P,N), 0 for others /* Line three in field file */ - real r_pi_pi; - real p_lp2; - real b_o_131; - real b_o_132; - real b_o_133; + double r_pi_pi; + double p_lp2; + double b_o_131; + double b_o_132; + double b_o_133; /* Line four in the field file */ - real p_ovun2; - real p_val3; - real valency_val; - real p_val5; - real rcore2; - real ecore2; - real acore2; + double p_ovun2; + double p_val3; + double valency_val; + double p_val5; + double rcore2; + double ecore2; + double acore2; /* Line five in the ffield file, only for lgvdw yes */ - real lgcij; - real lgre; + double lgcij; + double lgre; } single_body_parameters; @@ -168,29 +167,29 @@ typedef struct /* Two Body Parameters */ typedef struct { /* Bond Order parameters */ - real p_bo1,p_bo2,p_bo3,p_bo4,p_bo5,p_bo6; - real r_s, r_p, r_pp; // r_o distances in BO formula - real p_boc3, p_boc4, p_boc5; + double p_bo1,p_bo2,p_bo3,p_bo4,p_bo5,p_bo6; + double r_s, r_p, r_pp; // r_o distances in BO formula + double p_boc3, p_boc4, p_boc5; /* Bond Energy parameters */ - real p_be1, p_be2; - real De_s, De_p, De_pp; + double p_be1, p_be2; + double De_s, De_p, De_pp; /* Over/Under coordination parameters */ - real p_ovun1; + double p_ovun1; /* Van der Waal interaction parameters */ - real D; - real alpha; - real r_vdW; - real gamma_w; - real rcore, ecore, acore; - real lgcij, lgre; + double D; + double alpha; + double r_vdW; + double gamma_w; + double rcore, ecore, acore; + double lgcij, lgre; /* electrostatic parameters */ - real gamma; // note: this parameter is gamma^-3 and not gamma. + double gamma; // note: this parameter is gamma^-3 and not gamma. - real v13cor, ovc; + double v13cor, ovc; } two_body_parameters; @@ -198,14 +197,14 @@ typedef struct { /* 3-body parameters */ typedef struct { /* valence angle */ - real theta_00; - real p_val1, p_val2, p_val4, p_val7; + double theta_00; + double p_val1, p_val2, p_val4, p_val7; /* penalty */ - real p_pen1; + double p_pen1; /* 3-body conjugation */ - real p_coa1; + double p_coa1; } three_body_parameters; @@ -218,20 +217,20 @@ typedef struct{ /* hydrogen-bond parameters */ typedef struct{ - real r0_hb, p_hb1, p_hb2, p_hb3; + double r0_hb, p_hb1, p_hb2, p_hb3; } hbond_parameters; /* 4-body parameters */ typedef struct { - real V1, V2, V3; + double V1, V2, V3; /* torsion angle */ - real p_tor1; + double p_tor1; /* 4-body conjugation */ - real p_cot1; + double p_cot1; } four_body_parameters; @@ -267,7 +266,7 @@ struct _reax_atom rvec f; // force rvec f_old; - real q; // charge + double q; // charge rvec4 s; // they take part in rvec4 t; // computing q @@ -287,7 +286,7 @@ typedef _reax_atom reax_atom; typedef struct { - real V; + double V; rvec min, max, box_norms; rtensor box, box_inv; @@ -299,7 +298,7 @@ typedef struct struct grid_cell { - real cutoff; + double cutoff; rvec min, max; ivec rel_box; @@ -332,7 +331,7 @@ typedef struct ivec native_str; ivec native_end; - real ghost_cut; + double ghost_cut; ivec ghost_span; ivec ghost_nonb_span; ivec ghost_hbond_span; @@ -372,10 +371,10 @@ typedef struct typedef struct { - real ghost_nonb; - real ghost_hbond; - real ghost_bond; - real ghost_cutoff; + double ghost_nonb; + double ghost_hbond; + double ghost_bond; + double ghost_cutoff; } boundary_cutoff; using LAMMPS_NS::Pair; @@ -399,7 +398,7 @@ struct _reax_system class Pair *pair_ptr; int my_bonds; int mincap; - real safezone, saferzone; + double safezone, saferzone; }; typedef _reax_system reax_system; @@ -421,7 +420,7 @@ typedef struct 5 : NPT (Parrinello-Rehman-Nose-Hoover) Anisotropic*/ int ensemble; int nsteps; - real dt; + double dt; int geo_format; int restart; @@ -431,33 +430,33 @@ typedef struct int reposition_atoms; int reneighbor; - real vlist_cut; - real bond_cut; - real nonb_cut, nonb_low; - real hbond_cut; - real user_ghost_cut; + double vlist_cut; + double bond_cut; + double nonb_cut, nonb_low; + double hbond_cut; + double user_ghost_cut; - real bg_cut; - real bo_cut; - real thb_cut; - real thb_cutsq; + double bg_cut; + double bo_cut; + double thb_cut; + double thb_cutsq; int tabulate; int qeq_freq; - real q_err; + double q_err; int refactor; - real droptol; + double droptol; - real T_init, T_final, T; - real Tau_T; + double T_init, T_final, T; + double Tau_T; int T_mode; - real T_rate, T_freq; + double T_rate, T_freq; int virial; rvec P, Tau_P, Tau_PT; int press_mode; - real compressibility; + double compressibility; int molecular_analysis; int num_ignored; @@ -476,22 +475,22 @@ typedef struct typedef struct { - real T; - real xi; - real v_xi; - real v_xi_old; - real G_xi; + double T; + double xi; + double v_xi; + double v_xi_old; + double G_xi; } thermostat; typedef struct { - real P; - real eps; - real v_eps; - real v_eps_old; - real a_eps; + double P; + double eps; + double v_eps; + double v_eps_old; + double a_eps; } isotropic_barostat; @@ -499,12 +498,12 @@ typedef struct typedef struct { rtensor P; - real P_scalar; + double P_scalar; - real eps; - real v_eps; - real v_eps_old; - real a_eps; + double eps; + double v_eps; + double v_eps_old; + double a_eps; rtensor h0; rtensor v_g0; @@ -516,17 +515,17 @@ typedef struct typedef struct { - real start; - real end; - real elapsed; + double start; + double end; + double elapsed; - real total; - real comm; - real nbrs; - real init_forces; - real bonded; - real nonb; - real qEq; + double total; + double comm; + double nbrs; + double init_forces; + double bonded; + double nonb; + double qEq; int s_matvecs; int t_matvecs; } reax_timing; @@ -534,41 +533,41 @@ typedef struct typedef struct { - real e_tot; - real e_kin; // Total kinetic energy - real e_pot; + double e_tot; + double e_kin; // Total kinetic energy + double e_pot; - real e_bond; // Total bond energy - real e_ov; // Total over coordination - real e_un; // Total under coordination energy - real e_lp; // Total under coordination energy - real e_ang; // Total valance angle energy - real e_pen; // Total penalty energy - real e_coa; // Total three body conjgation energy - real e_hb; // Total Hydrogen bond energy - real e_tor; // Total torsional energy - real e_con; // Total four body conjugation energy - real e_vdW; // Total van der Waals energy - real e_ele; // Total electrostatics energy - real e_pol; // Polarization energy + double e_bond; // Total bond energy + double e_ov; // Total over coordination + double e_un; // Total under coordination energy + double e_lp; // Total under coordination energy + double e_ang; // Total valance angle energy + double e_pen; // Total penalty energy + double e_coa; // Total three body conjgation energy + double e_hb; // Total Hydrogen bond energy + double e_tor; // Total torsional energy + double e_con; // Total four body conjugation energy + double e_vdW; // Total van der Waals energy + double e_ele; // Total electrostatics energy + double e_pol; // Polarization energy } energy_data; typedef struct { int step; int prev_steps; - real time; + double time; - real M; // Total Mass - real inv_M; // 1 / Total Mass + double M; // Total Mass + double inv_M; // 1 / Total Mass rvec xcm; // Center of mass rvec vcm; // Center of mass velocity rvec fcm; // Center of mass force rvec amcm; // Angular momentum of CoM rvec avcm; // Angular velocity of CoM - real etran_cm; // Translational kinetic energy of CoM - real erot_cm; // Rotational kinetic energy of CoM + double etran_cm; // Translational kinetic energy of CoM + double erot_cm; // Rotational kinetic energy of CoM rtensor kinetic; // Kinetic energy tensor rtensor virial; // Hydrodynamic virial @@ -576,15 +575,15 @@ typedef struct energy_data my_en; energy_data sys_en; - real N_f; //Number of degrees of freedom + double N_f; //Number of degrees of freedom rvec t_scale; rtensor p_scale; thermostat therm; // Used in Nose_Hoover method isotropic_barostat iso_bar; flexible_barostat flex_bar; - real inv_W; + double inv_W; - real kin_press; + double kin_press; rvec int_press; rvec my_ext_press; rvec ext_press; @@ -597,7 +596,7 @@ typedef struct typedef struct{ int thb; int pthb; // pointer to the third body on the central atom's nbrlist - real theta, cos_theta; + double theta, cos_theta; rvec dcos_di, dcos_dj, dcos_dk; } three_body_interaction_data; @@ -605,7 +604,7 @@ typedef struct{ typedef struct { int nbr; ivec rel_box; - real d; + double d; rvec dvec; } far_neighbor_data; @@ -629,11 +628,11 @@ typedef struct{ } dbond_data; typedef struct{ - real BO, BO_s, BO_pi, BO_pi2; - real Cdbo, Cdbopi, Cdbopi2; - real C1dbo, C2dbo, C3dbo; - real C1dbopi, C2dbopi, C3dbopi, C4dbopi; - real C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2; + double BO, BO_s, BO_pi, BO_pi2; + double Cdbo, Cdbopi, Cdbopi2; + double C1dbo, C2dbo, C3dbo; + double C1dbopi, C2dbopi, C3dbopi, C4dbopi; + double C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2; rvec dBOp, dln_BOp_s, dln_BOp_pi, dln_BOp_pi2; } bond_order_data; @@ -643,7 +642,7 @@ typedef struct { int dbond_index; ivec rel_box; // rvec ext_factor; - real d; + double d; rvec dvec; bond_order_data bo_data; } bond_data; @@ -651,7 +650,7 @@ typedef struct { typedef struct { int j; - real val; + double val; } sparse_matrix_entry; typedef struct { @@ -676,35 +675,35 @@ typedef struct int allocated; /* communication storage */ - real *tmp_dbl[REAX_MAX_NBRS]; + double *tmp_dbl[REAX_MAX_NBRS]; rvec *tmp_rvec[REAX_MAX_NBRS]; rvec2 *tmp_rvec2[REAX_MAX_NBRS]; int *within_bond_box; /* bond order related storage */ - real *total_bond_order; - real *Deltap, *Deltap_boc; - real *Delta, *Delta_lp, *Delta_lp_temp, *Delta_e, *Delta_boc, *Delta_val; - real *dDelta_lp, *dDelta_lp_temp; - real *nlp, *nlp_temp, *Clp, *vlpex; + double *total_bond_order; + double *Deltap, *Deltap_boc; + double *Delta, *Delta_lp, *Delta_lp_temp, *Delta_e, *Delta_boc, *Delta_val; + double *dDelta_lp, *dDelta_lp_temp; + double *nlp, *nlp_temp, *Clp, *vlpex; rvec *dDeltap_self; int *bond_mark, *done_after; /* QEq storage */ sparse_matrix *H, *L, *U; - real *Hdia_inv, *b_s, *b_t, *b_prc, *b_prm, *s, *t; - real *droptol; + double *Hdia_inv, *b_s, *b_t, *b_prc, *b_prm, *s, *t; + double *droptol; rvec2 *b, *x; /* GMRES storage */ - real *y, *z, *g; - real *hc, *hs; - real **h, **v; + double *y, *z, *g; + double *hc, *hs; + double **h, **v; /* CG storage */ - real *r, *d, *q, *p; + double *r, *d, *q, *p; rvec2 *r2, *d2, *q2, *p2; /* Taper */ - real Tap[8]; //Tap7, Tap6, Tap5, Tap4, Tap3, Tap2, Tap1, Tap0; + double Tap[8]; //Tap7, Tap6, Tap5, Tap4, Tap3, Tap2, Tap1, Tap0; /* storage for analysis */ int *mark, *old_mark; @@ -718,7 +717,7 @@ typedef struct rvec *v_const; /* force calculations */ - real *CdDelta; // coefficient of dDelta + double *CdDelta; // coefficient of dDelta rvec *f; reallocate_data realloc; @@ -802,27 +801,27 @@ typedef struct typedef struct { - real H; - real e_vdW, CEvd; - real e_ele, CEclmb; + double H; + double e_vdW, CEvd; + double e_ele, CEclmb; } LR_data; typedef struct { - real a, b, c, d; + double a, b, c, d; } cubic_spline_coef; typedef struct { - real xmin, xmax; + double xmin, xmax; int n; - real dx, inv_dx; - real a; - real m; - real c; + double dx, inv_dx; + double a; + double m; + double c; LR_data *y; cubic_spline_coef *H; @@ -844,7 +843,7 @@ typedef void (*print_interaction)(reax_system*, control_params*, simulation_data*, storage*, reax_list**, output_controls*); -typedef real (*lookup_function)(real); +typedef double (*lookup_function)(double); typedef void (*message_sorter) (reax_system*, int, int, int, mpi_out_data*); typedef void (*unpacker) ( reax_system*, int, void*, int, neighbor_proc*, int ); diff --git a/src/USER-REAXC/reaxc_valence_angles.cpp b/src/USER-REAXC/reaxc_valence_angles.cpp index 60f985598d..c2b3287be5 100644 --- a/src/USER-REAXC/reaxc_valence_angles.cpp +++ b/src/USER-REAXC/reaxc_valence_angles.cpp @@ -30,9 +30,9 @@ #include "reaxc_list.h" #include "reaxc_vector.h" -static real Dot( real* v1, real* v2, int k ) +static double Dot( double* v1, double* v2, int k ) { - real ret = 0.0; + double ret = 0.0; for( int i=0; i < k; ++i ) ret += v1[i] * v2[i]; @@ -40,8 +40,8 @@ static real Dot( real* v1, real* v2, int k ) return ret; } -void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk, - real *theta, real *cos_theta ) +void Calculate_Theta( rvec dvec_ji, double d_ji, rvec dvec_jk, double d_jk, + double *theta, double *cos_theta ) { (*cos_theta) = Dot( dvec_ji, dvec_jk, 3 ) / ( d_ji * d_jk ); if( *cos_theta > 1. ) *cos_theta = 1.0; @@ -50,18 +50,18 @@ void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk, (*theta) = acos( *cos_theta ); } -void Calculate_dCos_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk, +void Calculate_dCos_Theta( rvec dvec_ji, double d_ji, rvec dvec_jk, double d_jk, rvec* dcos_theta_di, rvec* dcos_theta_dj, rvec* dcos_theta_dk ) { int t; - real sqr_d_ji = SQR(d_ji); - real sqr_d_jk = SQR(d_jk); - real inv_dists = 1.0 / (d_ji * d_jk); - real inv_dists3 = pow( inv_dists, 3.0 ); - real dot_dvecs = Dot( dvec_ji, dvec_jk, 3 ); - real Cdot_inv3 = dot_dvecs * inv_dists3; + double sqr_d_ji = SQR(d_ji); + double sqr_d_jk = SQR(d_jk); + double inv_dists = 1.0 / (d_ji * d_jk); + double inv_dists3 = pow( inv_dists, 3.0 ); + double dot_dvecs = Dot( dvec_ji, dvec_jk, 3 ); + double Cdot_inv3 = dot_dvecs * inv_dists3; for( t = 0; t < 3; ++t ) { (*dcos_theta_di)[t] = dvec_jk[t] * inv_dists - @@ -83,27 +83,27 @@ void Valence_Angles( reax_system *system, control_params *control, int start_j, end_j, start_pk, end_pk; int cnt, num_thb_intrs; - real temp, temp_bo_jt, pBOjt7; - real p_val1, p_val2, p_val3, p_val4, p_val5; - real p_val6, p_val7, p_val8, p_val9, p_val10; - real p_pen1, p_pen2, p_pen3, p_pen4; - real p_coa1, p_coa2, p_coa3, p_coa4; - real trm8, expval6, expval7, expval2theta, expval12theta, exp3ij, exp3jk; - real exp_pen2ij, exp_pen2jk, exp_pen3, exp_pen4, trm_pen34, exp_coa2; - real dSBO1, dSBO2, SBO, SBO2, CSBO2, SBOp, prod_SBO, vlpadj; - real CEval1, CEval2, CEval3, CEval4, CEval5, CEval6, CEval7, CEval8; - real CEpen1, CEpen2, CEpen3; - real e_ang, e_coa, e_pen; - real CEcoa1, CEcoa2, CEcoa3, CEcoa4, CEcoa5; - real Cf7ij, Cf7jk, Cf8j, Cf9j; - real f7_ij, f7_jk, f8_Dj, f9_Dj; - real Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta; - real BOA_ij, BOA_jk; + double temp, temp_bo_jt, pBOjt7; + double p_val1, p_val2, p_val3, p_val4, p_val5; + double p_val6, p_val7, p_val8, p_val9, p_val10; + double p_pen1, p_pen2, p_pen3, p_pen4; + double p_coa1, p_coa2, p_coa3, p_coa4; + double trm8, expval6, expval7, expval2theta, expval12theta, exp3ij, exp3jk; + double exp_pen2ij, exp_pen2jk, exp_pen3, exp_pen4, trm_pen34, exp_coa2; + double dSBO1, dSBO2, SBO, SBO2, CSBO2, SBOp, prod_SBO, vlpadj; + double CEval1, CEval2, CEval3, CEval4, CEval5, CEval6, CEval7, CEval8; + double CEpen1, CEpen2, CEpen3; + double e_ang, e_coa, e_pen; + double CEcoa1, CEcoa2, CEcoa3, CEcoa4, CEcoa5; + double Cf7ij, Cf7jk, Cf8j, Cf9j; + double f7_ij, f7_jk, f8_Dj, f9_Dj; + double Ctheta_0, theta_0, theta_00, theta, cos_theta, sin_theta; + double BOA_ij, BOA_jk; rvec force, ext_press; // Tallying variables - real eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; - real delij[3], delkj[3]; + double eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; + double delij[3], delkj[3]; three_body_header *thbh; three_body_parameters *thbp; diff --git a/src/USER-REAXC/reaxc_valence_angles.h b/src/USER-REAXC/reaxc_valence_angles.h index d7b79d7fc7..31936ba190 100644 --- a/src/USER-REAXC/reaxc_valence_angles.h +++ b/src/USER-REAXC/reaxc_valence_angles.h @@ -32,8 +32,8 @@ void Valence_Angles( reax_system*, control_params*, simulation_data*, storage*, reax_list**, output_controls* ); -void Calculate_Theta( rvec, real, rvec, real, real*, real* ); +void Calculate_Theta( rvec, double, rvec, double, double*, double* ); -void Calculate_dCos_Theta( rvec, real, rvec, real, rvec*, rvec*, rvec* ); +void Calculate_dCos_Theta( rvec, double, rvec, double, rvec*, rvec*, rvec* ); #endif diff --git a/src/USER-REAXC/reaxc_vector.cpp b/src/USER-REAXC/reaxc_vector.cpp index 2565064217..17a2851b00 100644 --- a/src/USER-REAXC/reaxc_vector.cpp +++ b/src/USER-REAXC/reaxc_vector.cpp @@ -34,7 +34,7 @@ void rvec_Copy( rvec dest, rvec src ) } -void rvec_Scale( rvec ret, real c, rvec v ) +void rvec_Scale( rvec ret, double c, rvec v ) { ret[0] = c * v[0], ret[1] = c * v[1], ret[2] = c * v[2]; } @@ -46,7 +46,7 @@ void rvec_Add( rvec ret, rvec v ) } -void rvec_ScaledAdd( rvec ret, real c, rvec v ) +void rvec_ScaledAdd( rvec ret, double c, rvec v ) { ret[0] += c * v[0], ret[1] += c * v[1], ret[2] += c * v[2]; } @@ -60,7 +60,7 @@ void rvec_Sum( rvec ret, rvec v1 ,rvec v2 ) } -void rvec_ScaledSum( rvec ret, real c1, rvec v1 ,real c2, rvec v2 ) +void rvec_ScaledSum( rvec ret, double c1, rvec v1 ,double c2, rvec v2 ) { ret[0] = c1 * v1[0] + c2 * v2[0]; ret[1] = c1 * v1[1] + c2 * v2[1]; @@ -68,13 +68,13 @@ void rvec_ScaledSum( rvec ret, real c1, rvec v1 ,real c2, rvec v2 ) } -real rvec_Dot( rvec v1, rvec v2 ) +double rvec_Dot( rvec v1, rvec v2 ) { return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]; } -real rvec_ScaledDot( real c1, rvec v1, real c2, rvec v2 ) +double rvec_ScaledDot( double c1, rvec v1, double c2, rvec v2 ) { return (c1*c2) * (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); } @@ -138,13 +138,13 @@ void rvec_OuterProduct( rtensor r, rvec v1, rvec v2 ) } -real rvec_Norm_Sqr( rvec v ) +double rvec_Norm_Sqr( rvec v ) { return SQR(v[0]) + SQR(v[1]) + SQR(v[2]); } -real rvec_Norm( rvec v ) +double rvec_Norm( rvec v ) { return sqrt( SQR(v[0]) + SQR(v[1]) + SQR(v[2]) ); } @@ -187,7 +187,7 @@ void rtensor_MatVec( rvec ret, rtensor m, rvec v ) } -void rtensor_Scale( rtensor ret, real c, rtensor m ) +void rtensor_Scale( rtensor ret, double c, rtensor m ) { int i, j; @@ -217,7 +217,7 @@ void ivec_Copy( ivec dest, ivec src ) } -void ivec_Scale( ivec dest, real C, ivec src ) +void ivec_Scale( ivec dest, double C, ivec src ) { dest[0] = (int)(C * src[0]); dest[1] = (int)(C * src[1]); diff --git a/src/USER-REAXC/reaxc_vector.h b/src/USER-REAXC/reaxc_vector.h index 189b866dc0..4aafff87e9 100644 --- a/src/USER-REAXC/reaxc_vector.h +++ b/src/USER-REAXC/reaxc_vector.h @@ -32,13 +32,13 @@ #include "reaxc_defs.h" void rvec_Copy( rvec, rvec ); -void rvec_Scale( rvec, real, rvec ); +void rvec_Scale( rvec, double, rvec ); void rvec_Add( rvec, rvec ); -void rvec_ScaledAdd( rvec, real, rvec ); +void rvec_ScaledAdd( rvec, double, rvec ); void rvec_Sum( rvec, rvec, rvec ); -void rvec_ScaledSum( rvec, real, rvec, real, rvec ); -real rvec_Dot( rvec, rvec ); -real rvec_ScaledDot( real, rvec, real, rvec ); +void rvec_ScaledSum( rvec, double, rvec, double, rvec ); +double rvec_Dot( rvec, rvec ); +double rvec_ScaledDot( double, rvec, double, rvec ); void rvec_Multiply( rvec, rvec, rvec ); void rvec_iMultiply( rvec, ivec, rvec ); void rvec_Divide( rvec, rvec, rvec ); @@ -46,19 +46,19 @@ void rvec_iDivide( rvec, rvec, ivec ); void rvec_Invert( rvec, rvec ); void rvec_Cross( rvec, rvec, rvec ); void rvec_OuterProduct( rtensor, rvec, rvec ); -real rvec_Norm_Sqr( rvec ); -real rvec_Norm( rvec ); +double rvec_Norm_Sqr( rvec ); +double rvec_Norm( rvec ); int rvec_isZero( rvec ); void rvec_MakeZero( rvec ); void rvec_Random( rvec ); void rtensor_MakeZero( rtensor ); void rtensor_MatVec( rvec, rtensor, rvec ); -void rtensor_Scale( rtensor, real, rtensor ); +void rtensor_Scale( rtensor, double, rtensor ); void ivec_MakeZero( ivec ); void ivec_Copy( ivec, ivec ); -void ivec_Scale( ivec, real, ivec ); +void ivec_Scale( ivec, double, ivec ); void ivec_Sum( ivec, ivec, ivec ); #endif diff --git a/src/USER-SMD/compute_smd_tlsph_strain_rate.cpp b/src/USER-SMD/compute_smd_tlsph_strain_rate.cpp index a9e909db10..4298793aef 100644 --- a/src/USER-SMD/compute_smd_tlsph_strain_rate.cpp +++ b/src/USER-SMD/compute_smd_tlsph_strain_rate.cpp @@ -37,6 +37,7 @@ using namespace Eigen; using namespace LAMMPS_NS; + /* ---------------------------------------------------------------------- */ ComputeSMDTLSPHStrainRate::ComputeSMDTLSPHStrainRate(LAMMPS *lmp, int narg, char **arg) : @@ -73,7 +74,6 @@ void ComputeSMDTLSPHStrainRate::init() { void ComputeSMDTLSPHStrainRate::compute_peratom() { invoked_peratom = update->ntimestep; - int *mol = atom->molecule; // grow vector array if necessary @@ -84,8 +84,6 @@ void ComputeSMDTLSPHStrainRate::compute_peratom() { array_atom = strain_rate_array; } - int ulsph_flag = 0; - int tlsph_flag = 0; int itmp = 0; Matrix3d *D = (Matrix3d *) force->pair->extract("smd/tlsph/strain_rate_ptr", itmp); if (D == NULL) { @@ -114,13 +112,3 @@ double ComputeSMDTLSPHStrainRate::memory_usage() { double bytes = size_peratom_cols * nmax * sizeof(double); return bytes; } - -/* - * deviator of a tensor - */ -Matrix3d ComputeSMDTLSPHStrainRate::Deviator(Matrix3d M) { - Matrix3d eye; - eye.setIdentity(); - eye *= M.trace() / 3.0; - return M - eye; -} diff --git a/src/USER-SMD/compute_smd_tlsph_strain_rate.h b/src/USER-SMD/compute_smd_tlsph_strain_rate.h index 6332f03764..86da271b45 100644 --- a/src/USER-SMD/compute_smd_tlsph_strain_rate.h +++ b/src/USER-SMD/compute_smd_tlsph_strain_rate.h @@ -33,8 +33,6 @@ ComputeStyle(smd/tlsph_strain_rate,ComputeSMDTLSPHStrainRate) #define LMP_COMPUTE_SMD_TLSPH_STRAIN_RATE_H #include "compute.h" -#include -using namespace Eigen; namespace LAMMPS_NS { @@ -45,7 +43,6 @@ class ComputeSMDTLSPHStrainRate : public Compute { void init(); void compute_peratom(); double memory_usage(); - Matrix3d Deviator(Matrix3d); private: int nmax; diff --git a/src/USER-SMD/compute_smd_tlsph_stress.cpp b/src/USER-SMD/compute_smd_tlsph_stress.cpp index bd38b6af33..e5eccd9aeb 100644 --- a/src/USER-SMD/compute_smd_tlsph_stress.cpp +++ b/src/USER-SMD/compute_smd_tlsph_stress.cpp @@ -36,6 +36,17 @@ using namespace Eigen; using namespace LAMMPS_NS; + +/* + * deviator of a tensor + */ +static Matrix3d Deviator(Matrix3d M) { + Matrix3d eye; + eye.setIdentity(); + eye *= M.trace() / 3.0; + return M - eye; +} + /* ---------------------------------------------------------------------- */ ComputeSMDTLSPHStress::ComputeSMDTLSPHStress(LAMMPS *lmp, int narg, char **arg) : @@ -119,13 +130,3 @@ double ComputeSMDTLSPHStress::memory_usage() { double bytes = size_peratom_cols * nmax * sizeof(double); return bytes; } - -/* - * deviator of a tensor - */ -Matrix3d ComputeSMDTLSPHStress::Deviator(Matrix3d M) { - Matrix3d eye; - eye.setIdentity(); - eye *= M.trace() / 3.0; - return M - eye; -} diff --git a/src/USER-SMD/compute_smd_tlsph_stress.h b/src/USER-SMD/compute_smd_tlsph_stress.h index 5348584d04..ed3ce7e7a9 100644 --- a/src/USER-SMD/compute_smd_tlsph_stress.h +++ b/src/USER-SMD/compute_smd_tlsph_stress.h @@ -33,8 +33,6 @@ ComputeStyle(smd/tlsph_stress, ComputeSMDTLSPHStress) #define LMP_COMPUTE_SMD_TLSPH_STRESS_H #include "compute.h" -#include -using namespace Eigen; namespace LAMMPS_NS { @@ -45,7 +43,6 @@ class ComputeSMDTLSPHStress : public Compute { void init(); void compute_peratom(); double memory_usage(); - Matrix3d Deviator(Matrix3d); private: int nmax; diff --git a/src/USER-SMD/compute_smd_ulsph_strain_rate.cpp b/src/USER-SMD/compute_smd_ulsph_strain_rate.cpp index 542c87bc32..b341fdf31a 100644 --- a/src/USER-SMD/compute_smd_ulsph_strain_rate.cpp +++ b/src/USER-SMD/compute_smd_ulsph_strain_rate.cpp @@ -36,6 +36,7 @@ using namespace Eigen; using namespace LAMMPS_NS; + /* ---------------------------------------------------------------------- */ ComputeSMDULSPHStrainRate::ComputeSMDULSPHStrainRate(LAMMPS *lmp, int narg, char **arg) : @@ -120,13 +121,3 @@ double ComputeSMDULSPHStrainRate::memory_usage() { double bytes = size_peratom_cols * nmax * sizeof(double); return bytes; } - -/* - * deviator of a tensor - */ -Matrix3d ComputeSMDULSPHStrainRate::Deviator(Matrix3d M) { - Matrix3d eye; - eye.setIdentity(); - eye *= M.trace() / 3.0; - return M - eye; -} diff --git a/src/USER-SMD/compute_smd_ulsph_strain_rate.h b/src/USER-SMD/compute_smd_ulsph_strain_rate.h index b200e15722..b712ab3c3e 100644 --- a/src/USER-SMD/compute_smd_ulsph_strain_rate.h +++ b/src/USER-SMD/compute_smd_ulsph_strain_rate.h @@ -33,8 +33,6 @@ ComputeStyle(smd/ulsph_strain_rate,ComputeSMDULSPHStrainRate) #define LMP_COMPUTE_SMD_ULSPH_STRAIN_RATE_H #include "compute.h" -#include -using namespace Eigen; namespace LAMMPS_NS { @@ -45,7 +43,6 @@ class ComputeSMDULSPHStrainRate : public Compute { void init(); void compute_peratom(); double memory_usage(); - Matrix3d Deviator(Matrix3d); private: int nmax; diff --git a/src/USER-SMD/compute_smd_ulsph_stress.cpp b/src/USER-SMD/compute_smd_ulsph_stress.cpp index 26b323e7ee..dc723fc639 100644 --- a/src/USER-SMD/compute_smd_ulsph_stress.cpp +++ b/src/USER-SMD/compute_smd_ulsph_stress.cpp @@ -36,6 +36,16 @@ using namespace Eigen; using namespace LAMMPS_NS; +/* + * deviator of a tensor + */ +static Matrix3d Deviator(Matrix3d M) { + Matrix3d eye; + eye.setIdentity(); + eye *= M.trace() / 3.0; + return M - eye; +} + /* ---------------------------------------------------------------------- */ ComputeSMDULSPHStress::ComputeSMDULSPHStress(LAMMPS *lmp, int narg, char **arg) : @@ -122,12 +132,3 @@ double ComputeSMDULSPHStress::memory_usage() { return bytes; } -/* - * deviator of a tensor - */ -Matrix3d ComputeSMDULSPHStress::Deviator(Matrix3d M) { - Matrix3d eye; - eye.setIdentity(); - eye *= M.trace() / 3.0; - return M - eye; -} diff --git a/src/USER-SMD/compute_smd_ulsph_stress.h b/src/USER-SMD/compute_smd_ulsph_stress.h index 5b31a111a4..4656076e48 100644 --- a/src/USER-SMD/compute_smd_ulsph_stress.h +++ b/src/USER-SMD/compute_smd_ulsph_stress.h @@ -33,8 +33,6 @@ ComputeStyle(smd/ulsph_stress, ComputeSMDULSPHStress) #define LMP_COMPUTE_SMD_ULSPH_STRESS_H #include "compute.h" -#include -using namespace Eigen; namespace LAMMPS_NS { @@ -45,7 +43,6 @@ class ComputeSMDULSPHStress : public Compute { void init(); void compute_peratom(); double memory_usage(); - Matrix3d Deviator(Matrix3d); private: int nmax; diff --git a/src/USER-SMD/fix_smd_move_triangulated_surface.h b/src/USER-SMD/fix_smd_move_triangulated_surface.h index e6593ea8b0..ce4eaed88e 100644 --- a/src/USER-SMD/fix_smd_move_triangulated_surface.h +++ b/src/USER-SMD/fix_smd_move_triangulated_surface.h @@ -31,10 +31,8 @@ FixStyle(smd/move_tri_surf,FixSMDMoveTriSurf) #ifndef LMP_FIX_SMD_INTEGRATE_TRIANGULAR_SURFACE_H #define LMP_FIX_SMD_INTEGRATE_TRIANGULAR_SURFACE_H -#include "fix.h" #include - -using namespace Eigen; +#include "fix.h" namespace LAMMPS_NS { @@ -53,9 +51,9 @@ protected: double dtv; bool linearFlag, wiggleFlag, rotateFlag; double vx, vy, vz; - Vector3d rotation_axis, origin; + Eigen::Vector3d rotation_axis, origin; double rotation_period; - Matrix3d u_cross, uxu; + Eigen::Matrix3d u_cross, uxu; double wiggle_travel, wiggle_max_travel, wiggle_direction; }; diff --git a/src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp b/src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp index bbfd1ec22b..5c20b803e8 100644 --- a/src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp +++ b/src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp @@ -202,7 +202,6 @@ void FixSMD_TLSPH_ReferenceConfiguration::pre_exchange() { void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) { int i, j, ii, jj, n, inum, jnum; int *ilist, *jlist, *numneigh, **firstneigh; - int itype, jtype; double r, h, wf, wfd; Vector3d dx; @@ -219,7 +218,6 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) { double **x0 = atom->x; double *radius = atom->radius; - int *type = atom->type; int *mask = atom->mask; tagint *tag = atom->tag; NeighList *list = pair->list; @@ -234,14 +232,12 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) { for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; - jtype = type[j]; if (INSERT_PREDEFINED_CRACKS) { if (!crack_exclude(i, j)) @@ -284,7 +280,6 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) { for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; @@ -296,7 +291,6 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) { dx(1) = x0[i][1] - x0[j][1]; dx(2) = x0[i][2] - x0[j][2]; r = dx.norm(); - jtype = type[j]; h = radius[i] + radius[j]; if (INSERT_PREDEFINED_CRACKS) { diff --git a/src/USER-SMD/fix_smd_wall_surface.cpp b/src/USER-SMD/fix_smd_wall_surface.cpp index 16bd1c7809..776acff7a8 100644 --- a/src/USER-SMD/fix_smd_wall_surface.cpp +++ b/src/USER-SMD/fix_smd_wall_surface.cpp @@ -31,8 +31,8 @@ #include "error.h" #include #include -#include #include "atom_vec.h" +#include using namespace LAMMPS_NS; using namespace FixConst; @@ -60,7 +60,7 @@ FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) : if (narg != 6) error->all(FLERR, "Illegal number of arguments for fix smd/wall_surface"); - filename.assign(arg[3]); + filename = strdup(arg[3]); wall_particle_type = force->inumeric(FLERR, arg[4]); wall_molecule_id = force->inumeric(FLERR, arg[5]); if (wall_molecule_id < 65535) { @@ -69,7 +69,7 @@ FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) : if (comm->me == 0) { printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n"); - printf("fix smd/wall_surface reads trianglulated surface from file: %s\n", filename.c_str()); + printf("fix smd/wall_surface reads trianglulated surface from file: %s\n", filename); printf("fix smd/wall_surface has particle type %d \n", wall_particle_type); printf("fix smd/wall_surface has molecule id %d \n", wall_molecule_id); printf(">>========>>========>>========>>========>>========>>========>>========>>========\n"); @@ -79,6 +79,8 @@ FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ FixSMDWallSurface::~FixSMDWallSurface() { + free(filename); + filename = NULL; // unregister this fix so atom class doesn't invoke it any more //atom->delete_callback(id, 0); @@ -237,10 +239,10 @@ void FixSMDWallSurface::read_triangles(int pass) { vert = new Vector3d[3]; Vector3d normal, center; - FILE *fp = fopen(filename.c_str(), "r"); + FILE *fp = fopen(filename, "r"); if (fp == NULL) { char str[128]; - sprintf(str, "Cannot open file %s", filename.c_str()); + sprintf(str, "Cannot open file %s", filename); error->one(FLERR, str); } diff --git a/src/USER-SMD/fix_smd_wall_surface.h b/src/USER-SMD/fix_smd_wall_surface.h index 32fb934d63..8bd7002a9e 100644 --- a/src/USER-SMD/fix_smd_wall_surface.h +++ b/src/USER-SMD/fix_smd_wall_surface.h @@ -21,8 +21,6 @@ FixStyle(smd/wall_surface,FixSMDWallSurface) #define LMP_FIX_SMD_WALL_SURFACE_H #include "fix.h" -#include -using namespace std; namespace LAMMPS_NS { @@ -42,7 +40,7 @@ public: private: int first; // flag for first time initialization double sublo[3], subhi[3]; // epsilon-extended proc sub-box for adding atoms; - std::string filename; + char *filename; int wall_particle_type; int wall_molecule_id; }; diff --git a/src/USER-SMD/pair_smd_tlsph.cpp b/src/USER-SMD/pair_smd_tlsph.cpp index 73adf2ed31..77a45ee754 100644 --- a/src/USER-SMD/pair_smd_tlsph.cpp +++ b/src/USER-SMD/pair_smd_tlsph.cpp @@ -2102,7 +2102,6 @@ void PairTlsph::ComputeDamage(const int i, const Matrix3d strain, const Matrix3d double *damage = atom->damage; int *type = atom->type; int itype = type[i]; - bool damage_flag = false; double jc_failure_strain; //double damage_gap, damage_rate; Matrix3d eye, stress_deviator; @@ -2116,16 +2115,15 @@ void PairTlsph::ComputeDamage(const int i, const Matrix3d strain, const Matrix3d /* * maximum stress failure criterion: */ - damage_flag = IsotropicMaxStressDamage(stress, Lookup[FAILURE_MAX_PRINCIPAL_STRESS_THRESHOLD][itype]); + IsotropicMaxStressDamage(stress, Lookup[FAILURE_MAX_PRINCIPAL_STRESS_THRESHOLD][itype]); } else if (failureModel[itype].failure_max_principal_strain) { error->one(FLERR, "not yet implemented"); /* * maximum strain failure criterion: */ - damage_flag = IsotropicMaxStrainDamage(strain, Lookup[FAILURE_MAX_PRINCIPAL_STRAIN_THRESHOLD][itype]); + IsotropicMaxStrainDamage(strain, Lookup[FAILURE_MAX_PRINCIPAL_STRAIN_THRESHOLD][itype]); } else if (failureModel[itype].failure_max_plastic_strain) { if (eff_plastic_strain[i] >= Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype]) { - damage_flag = true; damage[i] = 1.0; //double damage_gap = 0.5 * Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype]; //damage[i] = (eff_plastic_strain[i] - Lookup[FAILURE_MAX_PLASTIC_STRAIN_THRESHOLD][itype]) / damage_gap; @@ -2142,7 +2140,6 @@ void PairTlsph::ComputeDamage(const int i, const Matrix3d strain, const Matrix3d //printf("JC failure strain is: %f\n", jc_failure_strain); if (eff_plastic_strain[i] >= jc_failure_strain) { - damage_flag = true; double damage_rate = Lookup[SIGNAL_VELOCITY][itype] / (100.0 * radius[i]); damage[i] += damage_rate * update->dt; //damage[i] = 1.0; diff --git a/src/USER-SMD/pair_smd_tlsph.h b/src/USER-SMD/pair_smd_tlsph.h index b02a1c585e..17db11b816 100644 --- a/src/USER-SMD/pair_smd_tlsph.h +++ b/src/USER-SMD/pair_smd_tlsph.h @@ -33,12 +33,7 @@ PairStyle(smd/tlsph,PairTlsph) #include "pair.h" #include -#include -#include -#include -using namespace std; -using namespace Eigen; namespace LAMMPS_NS { class PairTlsph: public Pair { @@ -67,13 +62,13 @@ public: void PreCompute(); void ComputeForces(int eflag, int vflag); void effective_longitudinal_modulus(const int itype, const double dt, const double d_iso, const double p_rate, - const Matrix3d d_dev, const Matrix3d sigma_dev_rate, const double damage, double &K_eff, double &mu_eff, double &M_eff); + const Eigen::Matrix3d d_dev, const Eigen::Matrix3d sigma_dev_rate, const double damage, double &K_eff, double &mu_eff, double &M_eff); void ComputePressure(const int i, const double rho, const double mass_specific_energy, const double vol_specific_energy, const double pInitial, const double d_iso, double &pFinal, double &p_rate); - void ComputeStressDeviator(const int i, const Matrix3d sigmaInitial_dev, const Matrix3d d_dev, Matrix3d &sigmaFinal_dev, - Matrix3d &sigma_dev_rate, double &plastic_strain_increment); - void ComputeDamage(const int i, const Matrix3d strain, const Matrix3d sigmaFinal, Matrix3d &sigma_damaged); + void ComputeStressDeviator(const int i, const Eigen::Matrix3d sigmaInitial_dev, const Eigen::Matrix3d d_dev, Eigen::Matrix3d &sigmaFinal_dev, + Eigen::Matrix3d &sigma_dev_rate, double &plastic_strain_increment); + void ComputeDamage(const int i, const Eigen::Matrix3d strain, const Eigen::Matrix3d sigmaFinal, Eigen::Matrix3d &sigma_damaged); protected: @@ -89,12 +84,12 @@ protected: /* * per atom arrays */ - Matrix3d *K, *PK1, *Fdot, *Fincr; - Matrix3d *R; // rotation matrix - Matrix3d *FincrInv; - Matrix3d *D, *W; // strain rate and spin tensor - Vector3d *smoothVelDifference; - Matrix3d *CauchyStress; + Eigen::Matrix3d *K, *PK1, *Fdot, *Fincr; + Eigen::Matrix3d *R; // rotation matrix + Eigen::Matrix3d *FincrInv; + Eigen::Matrix3d *D, *W; // strain rate and spin tensor + Eigen::Vector3d *smoothVelDifference; + Eigen::Matrix3d *CauchyStress; double *detF, *particle_dt; double *hourglass_error; int *numNeighsRefConfig; diff --git a/src/USER-SMD/pair_smd_triangulated_surface.h b/src/USER-SMD/pair_smd_triangulated_surface.h index 1cfd88f1fa..9f92e05ca4 100644 --- a/src/USER-SMD/pair_smd_triangulated_surface.h +++ b/src/USER-SMD/pair_smd_triangulated_surface.h @@ -34,7 +34,6 @@ PairStyle(smd/tri_surface,PairTriSurf) #include "pair.h" #include -using namespace Eigen; namespace LAMMPS_NS { @@ -50,8 +49,8 @@ class PairTriSurf : public Pair { void init_style(); void init_list(int, class NeighList *); virtual double memory_usage(); - void PointTriangleDistance(const Vector3d P, const Vector3d TRI1, const Vector3d TRI2, const Vector3d TRI3, - Vector3d &CP, double &dist); + void PointTriangleDistance(const Eigen::Vector3d P, const Eigen::Vector3d TRI1, const Eigen::Vector3d TRI2, const Eigen::Vector3d TRI3, + Eigen::Vector3d &CP, double &dist); double clamp(const double a, const double min, const double max); void *extract(const char *, int &); diff --git a/src/USER-SMD/pair_smd_ulsph.cpp b/src/USER-SMD/pair_smd_ulsph.cpp index 92b6a6c0fd..9cb3bd2788 100644 --- a/src/USER-SMD/pair_smd_ulsph.cpp +++ b/src/USER-SMD/pair_smd_ulsph.cpp @@ -222,7 +222,7 @@ void PairULSPH::PreCompute() { int *ilist, *jlist, *numneigh; int **firstneigh; int nlocal = atom->nlocal; - int inum, jnum, ii, jj, i, itype, jtype, j, idim; + int inum, jnum, ii, jj, i, itype, j, idim; double wfd, h, irad, r, rSq, wf, ivol, jvol; Vector3d dx, dv, g, du; Matrix3d Ktmp, Ltmp, Ftmp, K3di, D; @@ -281,7 +281,6 @@ void PairULSPH::PreCompute() { if (rSq < h * h) { r = sqrt(rSq); - jtype = type[j]; jvol = vfrac[j]; // distance vectors in current and reference configuration, velocity difference @@ -649,7 +648,6 @@ void PairULSPH::AssembleStressTensor() { double **tlsph_stress = atom->smd_stress; double *e = atom->e; int *type = atom->type; - double pFinal; int i, itype; int nlocal = atom->nlocal; Matrix3d D, Ddev, W, V, sigma_diag; @@ -686,7 +684,6 @@ void PairULSPH::AssembleStressTensor() { error->one(FLERR, "unknown EOS."); break; case NONE: - pFinal = 0.0; c0[i] = 1.0; break; case EOS_TAIT: diff --git a/src/USER-SMD/pair_smd_ulsph.h b/src/USER-SMD/pair_smd_ulsph.h index 81919afdb0..40ccc37e93 100644 --- a/src/USER-SMD/pair_smd_ulsph.h +++ b/src/USER-SMD/pair_smd_ulsph.h @@ -33,11 +33,7 @@ PairStyle(smd/ulsph,PairULSPH) #include "pair.h" #include -#include -#include -using namespace Eigen; -using namespace std; namespace LAMMPS_NS { class PairULSPH: public Pair { @@ -57,10 +53,10 @@ public: void *extract(const char *, int &); void PreCompute(); void PreCompute_DensitySummation(); - double effective_shear_modulus(const Matrix3d d_dev, const Matrix3d deltaStressDev, const double dt, const int itype); + double effective_shear_modulus(const Eigen::Matrix3d d_dev, const Eigen::Matrix3d deltaStressDev, const double dt, const int itype); - Vector3d ComputeHourglassForce(const int i, const int itype, const int j, const int jtype, const Vector3d dv, - const Vector3d xij, const Vector3d g, const double c_ij, const double mu_ij, const double rho_ij); + Eigen::Vector3d ComputeHourglassForce(const int i, const int itype, const int j, const int jtype, const Eigen::Vector3d dv, + const Eigen::Vector3d xij, const Eigen::Vector3d g, const double c_ij, const double mu_ij, const double rho_ij); protected: @@ -78,10 +74,10 @@ protected: int nmax; // max number of atoms on this proc int *numNeighs; - Matrix3d *K; + Eigen::Matrix3d *K; double *shepardWeight, *c0, *rho; - Vector3d *smoothVel; - Matrix3d *stressTensor, *L, *F; + Eigen::Vector3d *smoothVel; + Eigen::Matrix3d *stressTensor, *L, *F; double dtCFL;