diff --git a/src/BODY/fix_nh_body.cpp b/src/BODY/fix_nh_body.cpp index a6255a486d..a1865dfd18 100644 --- a/src/BODY/fix_nh_body.cpp +++ b/src/BODY/fix_nh_body.cpp @@ -104,7 +104,6 @@ void FixNHBody::nve_x() AtomVecBody::Bonus *bonus = avec->bonus; int *body = atom->body; double **angmom = atom->angmom; - double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; if (igroup == atom->firstgroup) nlocal = atom->nfirst; diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 12f39bf931..49728d6aa6 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -3184,12 +3184,11 @@ double PairAIREBO::piRCSpline(double Nij, double Nji, double Nijconj, if (done==0) { for (i=0; i=(double) i && Nij<=(double) i+1 || Nij==(double) i) x=i; + if (Nij>=(double) i && Nij<=(double) i+1) x=i; for (i=0; i=(double) i && Nji<=(double) i+1 || Nji==(double) i) y=i; + if (Nji>=(double) i && Nji<=(double) i+1) y=i; for (i=0; i=(double) i && Nijconj<=(double) i+1 || - Nijconj==(double) i) z=i; + if (Nijconj>=(double) i && Nijconj<=(double) i+1) z=i; for (i=0; i<64; i++) coeffs[i]=piCC[x][y][z][i]; piRC=Sptricubic(Nij,Nji,Nijconj,coeffs,dN3); @@ -3199,7 +3198,7 @@ double PairAIREBO::piRCSpline(double Nij, double Nji, double Nijconj, // CH interaction - if (typei==0 && typej==1 || typei==1 && typej==0) { + if ((typei==0 && typej==1) || (typei==1 && typej==0)) { // if the inputs are out of bounds set them back to a point in bounds if (NijpiCHdom[0][1] || diff --git a/src/RIGID/fix_rigid_nh.cpp b/src/RIGID/fix_rigid_nh.cpp index 6a1476ba35..5db799e9f0 100644 --- a/src/RIGID/fix_rigid_nh.cpp +++ b/src/RIGID/fix_rigid_nh.cpp @@ -84,7 +84,7 @@ FixRigidNH::FixRigidNH(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Cannot use fix rigid npt/nph on a non-periodic dimension"); - if (pcouple == XYZ && dimension == 3 && + if (pcouple == XYZ && dimension == 3 && (p_start[0] != p_start[1] || p_start[0] != p_start[2] || p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] || p_period[0] != p_period[1] || p_period[0] != p_period[2])) diff --git a/src/USER-MISC/fix_ttm_mod.cpp b/src/USER-MISC/fix_ttm_mod.cpp index f54c21532b..5a58dd8489 100644 --- a/src/USER-MISC/fix_ttm_mod.cpp +++ b/src/USER-MISC/fix_ttm_mod.cpp @@ -613,110 +613,110 @@ void FixTTMMod::end_of_step() double stability_criterion = 0.0; for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - T_electron_first[ixnode][iynode][iznode] = - T_electron[ixnode][iynode][iznode]; + for (int iynode = 0; iynode < nynodes; iynode++) + for (int iznode = 0; iznode < nznodes; iznode++) + T_electron_first[ixnode][iynode][iznode] = + T_electron[ixnode][iynode][iznode]; do { - for (int ixnode = 0; ixnode < nxnodes; ixnode++) + for (int ixnode = 0; ixnode < nxnodes; ixnode++) for (int iynode = 0; iynode < nynodes; iynode++) for (int iznode = 0; iznode < nznodes; iznode++) T_electron[ixnode][iynode][iznode] = T_electron_first[ixnode][iynode][iznode]; - stability_criterion = 1.0 - - 2.0*inner_dt/el_specific_heat * - (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz)); - if (stability_criterion < 0.0) { - inner_dt = 0.25*el_specific_heat / - (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz)); - } - num_inner_timesteps = static_cast(update->dt/inner_dt) + 1; - inner_dt = update->dt/double(num_inner_timesteps); - if (num_inner_timesteps > 1000000) - error->warning(FLERR,"Too many inner timesteps in fix ttm",0); - for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps; - ith_inner_timestep++) { - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) - T_electron_old[ixnode][iynode][iznode] = - T_electron[ixnode][iynode][iznode]; - // compute new electron T profile - duration = duration + inner_dt; - for (int ixnode = 0; ixnode < nxnodes; ixnode++) - for (int iynode = 0; iynode < nynodes; iynode++) - for (int iznode = 0; iznode < nznodes; iznode++) { - int right_xnode = ixnode + 1; - int right_ynode = iynode + 1; - int right_znode = iznode + 1; - if (right_xnode == nxnodes) right_xnode = 0; - if (right_ynode == nynodes) right_ynode = 0; - if (right_znode == nznodes) right_znode = 0; - int left_xnode = ixnode - 1; - int left_ynode = iynode - 1; - int left_znode = iznode - 1; - if (left_xnode == -1) left_xnode = nxnodes - 1; - if (left_ynode == -1) left_ynode = nynodes - 1; - if (left_znode == -1) left_znode = nznodes - 1; - double skin_layer_d = double(skin_layer); - double ixnode_d = double(ixnode); - double surface_d = double(t_surface_l); - mult_factor = 0.0; - if (duration < width){ - if (ixnode >= t_surface_l) mult_factor = (intensity/(dx*skin_layer_d))*exp((-1.0)*(ixnode_d - surface_d)/skin_layer_d); - } - if (ixnode < t_surface_l) net_energy_transfer_all[ixnode][iynode][iznode] = 0.0; - double cr_vac = 1; - if (T_electron_old[ixnode][iynode][iznode] == 0) cr_vac = 0; - double cr_v_l_x = 1; - if (T_electron_old[left_xnode][iynode][iznode] == 0) cr_v_l_x = 0; - double cr_v_r_x = 1; - if (T_electron_old[right_xnode][iynode][iznode] == 0) cr_v_r_x = 0; - double cr_v_l_y = 1; - if (T_electron_old[ixnode][left_ynode][iznode] == 0) cr_v_l_y = 0; - double cr_v_r_y = 1; - if (T_electron_old[ixnode][right_ynode][iznode] == 0) cr_v_r_y = 0; - double cr_v_l_z = 1; - if (T_electron_old[ixnode][iynode][left_znode] == 0) cr_v_l_z = 0; - double cr_v_r_z = 1; - if (T_electron_old[ixnode][iynode][right_znode] == 0) cr_v_r_z = 0; - if (cr_vac != 0) { - T_electron[ixnode][iynode][iznode] = - T_electron_old[ixnode][iynode][iznode] + - inner_dt/el_properties(T_electron_old[ixnode][iynode][iznode]).el_heat_capacity * - ((cr_v_r_x*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[right_xnode][iynode][iznode]/2.0).el_thermal_conductivity* - (T_electron_old[right_xnode][iynode][iznode]-T_electron_old[ixnode][iynode][iznode])/dx - - cr_v_l_x*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[left_xnode][iynode][iznode]/2.0).el_thermal_conductivity* - (T_electron_old[ixnode][iynode][iznode]-T_electron_old[left_xnode][iynode][iznode])/dx)/dx + - (cr_v_r_y*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][right_ynode][iznode]/2.0).el_thermal_conductivity* - (T_electron_old[ixnode][right_ynode][iznode]-T_electron_old[ixnode][iynode][iznode])/dy - - cr_v_l_y*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][left_ynode][iznode]/2.0).el_thermal_conductivity* - (T_electron_old[ixnode][iynode][iznode]-T_electron_old[ixnode][left_ynode][iznode])/dy)/dy + - (cr_v_r_z*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][iynode][right_znode]/2.0).el_thermal_conductivity* - (T_electron_old[ixnode][iynode][right_znode]-T_electron_old[ixnode][iynode][iznode])/dz - - cr_v_l_z*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][iynode][left_znode]/2.0).el_thermal_conductivity* - (T_electron_old[ixnode][iynode][iznode]-T_electron_old[ixnode][iynode][left_znode])/dz)/dz); - T_electron[ixnode][iynode][iznode]+=inner_dt/el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity* - (mult_factor - - net_energy_transfer_all[ixnode][iynode][iznode]/del_vol); - } - else T_electron[ixnode][iynode][iznode] = - T_electron_old[ixnode][iynode][iznode]; - if ((T_electron[ixnode][iynode][iznode] > 0.0) && (T_electron[ixnode][iynode][iznode] < electron_temperature_min)) - T_electron[ixnode][iynode][iznode] = T_electron[ixnode][iynode][iznode] + 0.5*(electron_temperature_min - T_electron[ixnode][iynode][iznode]); + stability_criterion = 1.0 - + 2.0*inner_dt/el_specific_heat * + (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz)); + if (stability_criterion < 0.0) { + inner_dt = 0.25*el_specific_heat / + (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz)); + } + num_inner_timesteps = static_cast(update->dt/inner_dt) + 1; + inner_dt = update->dt/double(num_inner_timesteps); + if (num_inner_timesteps > 1000000) + error->warning(FLERR,"Too many inner timesteps in fix ttm",0); + for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps; + ith_inner_timestep++) { + for (int ixnode = 0; ixnode < nxnodes; ixnode++) + for (int iynode = 0; iynode < nynodes; iynode++) + for (int iznode = 0; iznode < nznodes; iznode++) + T_electron_old[ixnode][iynode][iznode] = + T_electron[ixnode][iynode][iznode]; + // compute new electron T profile + duration = duration + inner_dt; + for (int ixnode = 0; ixnode < nxnodes; ixnode++) + for (int iynode = 0; iynode < nynodes; iynode++) + for (int iznode = 0; iznode < nznodes; iznode++) { + int right_xnode = ixnode + 1; + int right_ynode = iynode + 1; + int right_znode = iznode + 1; + if (right_xnode == nxnodes) right_xnode = 0; + if (right_ynode == nynodes) right_ynode = 0; + if (right_znode == nznodes) right_znode = 0; + int left_xnode = ixnode - 1; + int left_ynode = iynode - 1; + int left_znode = iznode - 1; + if (left_xnode == -1) left_xnode = nxnodes - 1; + if (left_ynode == -1) left_ynode = nynodes - 1; + if (left_znode == -1) left_znode = nznodes - 1; + double skin_layer_d = double(skin_layer); + double ixnode_d = double(ixnode); + double surface_d = double(t_surface_l); + mult_factor = 0.0; + if (duration < width){ + if (ixnode >= t_surface_l) mult_factor = (intensity/(dx*skin_layer_d))*exp((-1.0)*(ixnode_d - surface_d)/skin_layer_d); + } + if (ixnode < t_surface_l) net_energy_transfer_all[ixnode][iynode][iznode] = 0.0; + double cr_vac = 1; + if (T_electron_old[ixnode][iynode][iznode] == 0) cr_vac = 0; + double cr_v_l_x = 1; + if (T_electron_old[left_xnode][iynode][iznode] == 0) cr_v_l_x = 0; + double cr_v_r_x = 1; + if (T_electron_old[right_xnode][iynode][iznode] == 0) cr_v_r_x = 0; + double cr_v_l_y = 1; + if (T_electron_old[ixnode][left_ynode][iznode] == 0) cr_v_l_y = 0; + double cr_v_r_y = 1; + if (T_electron_old[ixnode][right_ynode][iznode] == 0) cr_v_r_y = 0; + double cr_v_l_z = 1; + if (T_electron_old[ixnode][iynode][left_znode] == 0) cr_v_l_z = 0; + double cr_v_r_z = 1; + if (T_electron_old[ixnode][iynode][right_znode] == 0) cr_v_r_z = 0; + if (cr_vac != 0) { + T_electron[ixnode][iynode][iznode] = + T_electron_old[ixnode][iynode][iznode] + + inner_dt/el_properties(T_electron_old[ixnode][iynode][iznode]).el_heat_capacity * + ((cr_v_r_x*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[right_xnode][iynode][iznode]/2.0).el_thermal_conductivity* + (T_electron_old[right_xnode][iynode][iznode]-T_electron_old[ixnode][iynode][iznode])/dx - + cr_v_l_x*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[left_xnode][iynode][iznode]/2.0).el_thermal_conductivity* + (T_electron_old[ixnode][iynode][iznode]-T_electron_old[left_xnode][iynode][iznode])/dx)/dx + + (cr_v_r_y*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][right_ynode][iznode]/2.0).el_thermal_conductivity* + (T_electron_old[ixnode][right_ynode][iznode]-T_electron_old[ixnode][iynode][iznode])/dy - + cr_v_l_y*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][left_ynode][iznode]/2.0).el_thermal_conductivity* + (T_electron_old[ixnode][iynode][iznode]-T_electron_old[ixnode][left_ynode][iznode])/dy)/dy + + (cr_v_r_z*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][iynode][right_znode]/2.0).el_thermal_conductivity* + (T_electron_old[ixnode][iynode][right_znode]-T_electron_old[ixnode][iynode][iznode])/dz - + cr_v_l_z*el_properties(T_electron_old[ixnode][iynode][iznode]/2.0+T_electron_old[ixnode][iynode][left_znode]/2.0).el_thermal_conductivity* + (T_electron_old[ixnode][iynode][iznode]-T_electron_old[ixnode][iynode][left_znode])/dz)/dz); + T_electron[ixnode][iynode][iznode]+=inner_dt/el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity* + (mult_factor - + net_energy_transfer_all[ixnode][iynode][iznode]/del_vol); + } + else T_electron[ixnode][iynode][iznode] = + T_electron_old[ixnode][iynode][iznode]; + if ((T_electron[ixnode][iynode][iznode] > 0.0) && (T_electron[ixnode][iynode][iznode] < electron_temperature_min)) + T_electron[ixnode][iynode][iznode] = T_electron[ixnode][iynode][iznode] + 0.5*(electron_temperature_min - T_electron[ixnode][iynode][iznode]); - if (el_properties(T_electron[ixnode][iynode][iznode]).el_thermal_conductivity > el_thermal_conductivity) - el_thermal_conductivity = el_properties(T_electron[ixnode][iynode][iznode]).el_thermal_conductivity; - if ((T_electron[ixnode][iynode][iznode] > 0.0) && (el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity < el_specific_heat)) + if (el_properties(T_electron[ixnode][iynode][iznode]).el_thermal_conductivity > el_thermal_conductivity) + el_thermal_conductivity = el_properties(T_electron[ixnode][iynode][iznode]).el_thermal_conductivity; + if ((T_electron[ixnode][iynode][iznode] > 0.0) && (el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity < el_specific_heat)) el_specific_heat = el_properties(T_electron[ixnode][iynode][iznode]).el_heat_capacity; - } - } - stability_criterion = 1.0 - - 2.0*inner_dt/el_specific_heat * - (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz)); + } + } + stability_criterion = 1.0 - + 2.0*inner_dt/el_specific_heat * + (el_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz)); - } while (stability_criterion < 0.0); + } while (stability_criterion < 0.0); // output nodal temperatures for current timestep if ((nfileevery) && !(update->ntimestep % nfileevery)) { // compute atomic Ta for each grid point diff --git a/src/USER-REAXC/reaxc_tool_box.cpp b/src/USER-REAXC/reaxc_tool_box.cpp index fa1f8b8a75..22576e9f3b 100644 --- a/src/USER-REAXC/reaxc_tool_box.cpp +++ b/src/USER-REAXC/reaxc_tool_box.cpp @@ -54,7 +54,7 @@ int Tokenize( char* s, char*** tok ) } /* safe malloc */ -void *smalloc( long n, const char *name, MPI_Comm comm ) +void *smalloc( rc_bigint n, const char *name, MPI_Comm comm ) { void *ptr; @@ -77,19 +77,19 @@ void *smalloc( long n, const char *name, MPI_Comm comm ) /* safe calloc */ -void *scalloc( int n, int size, const char *name, MPI_Comm comm ) +void *scalloc( rc_bigint n, rc_bigint size, const char *name, MPI_Comm comm ) { void *ptr; if( n <= 0 ) { - fprintf( stderr, "WARNING: trying to allocate %d elements for array %s. ", + fprintf( stderr, "WARNING: trying to allocate %ld elements for array %s. ", n, name ); fprintf( stderr, "returning NULL.\n" ); return NULL; } if( size <= 0 ) { - fprintf( stderr, "WARNING: elements size for array %s is %d. ", + fprintf( stderr, "WARNING: elements size for array %s is %ld. ", name, size ); fprintf( stderr, "returning NULL.\n" ); return NULL; @@ -97,7 +97,7 @@ void *scalloc( int n, int size, const char *name, MPI_Comm comm ) ptr = calloc( n, size ); if( ptr == NULL ) { - fprintf( stderr, "ERROR: failed to allocate %d bytes for array %s", + fprintf( stderr, "ERROR: failed to allocate %ld bytes for array %s", n*size, name ); MPI_Abort( comm, INSUFFICIENT_MEMORY ); } diff --git a/src/USER-REAXC/reaxc_tool_box.h b/src/USER-REAXC/reaxc_tool_box.h index 90a711fdb7..7f8dd455ca 100644 --- a/src/USER-REAXC/reaxc_tool_box.h +++ b/src/USER-REAXC/reaxc_tool_box.h @@ -37,7 +37,7 @@ double Get_Time( ); int Tokenize( char*, char*** ); /* from lammps */ -void *smalloc( long, const char*, MPI_Comm ); -void *scalloc( int, int, const char*, MPI_Comm ); +void *smalloc( rc_bigint, const char*, MPI_Comm ); +void *scalloc( rc_bigint, rc_bigint, const char*, MPI_Comm ); void sfree( void*, const char* ); #endif diff --git a/src/USER-REAXC/reaxc_traj.cpp b/src/USER-REAXC/reaxc_traj.cpp index 1dcc2f5232..9d4fa73524 100644 --- a/src/USER-REAXC/reaxc_traj.cpp +++ b/src/USER-REAXC/reaxc_traj.cpp @@ -101,7 +101,7 @@ int Write_Header( reax_system *system, control_params *control, out_control->traj_title ); strncat( out_control->buffer, out_control->line, HEADER_LINE_LEN+1 ); - sprintf( out_control->line, INT_LINE, "number_of_atoms:", system->bigN ); + sprintf( out_control->line, BIGINT_LINE, "number_of_atoms:", system->bigN ); strncat( out_control->buffer, out_control->line, HEADER_LINE_LEN+1 ); sprintf( out_control->line, STR_LINE, "ensemble_type:", diff --git a/src/USER-REAXC/reaxc_traj.h b/src/USER-REAXC/reaxc_traj.h index e00341e5e1..72c56637eb 100644 --- a/src/USER-REAXC/reaxc_traj.h +++ b/src/USER-REAXC/reaxc_traj.h @@ -36,6 +36,7 @@ #define HEADER_LINE_LEN 62 #define STR_LINE "%-37s%-24s\n" #define INT_LINE "%-37s%-24d\n" +#define BIGINT_LINE "%-37s%-24ld\n" #define INT2_LINE "%-36s%-12d,%-12d\n" #define REAL_LINE "%-37s%-24.3f\n" #define SCI_LINE "%-37s%-24g\n" diff --git a/src/USER-REAXC/reaxc_types.h b/src/USER-REAXC/reaxc_types.h index 0779acaf39..9737b55d27 100644 --- a/src/USER-REAXC/reaxc_types.h +++ b/src/USER-REAXC/reaxc_types.h @@ -68,8 +68,9 @@ typedef double rvec2[2]; typedef double rvec4[4]; -// import LAMMPS' definition of tagint +// import LAMMPS' definition of tagint and bigint typedef LAMMPS_NS::tagint rc_tagint; +typedef LAMMPS_NS::bigint rc_bigint; typedef struct { @@ -383,7 +384,8 @@ struct _reax_system { reax_interaction reax_param; - int n, N, bigN, numH; + rc_bigint bigN; + int n, N, numH; int local_cap, total_cap, gcell_cap, Hcap; int est_recv, est_trans, max_recved; int wsize, my_rank, num_nbrs;