diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index c0ccb01cad..c885140711 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -553,7 +553,7 @@ void FixQEqReax::init_matvec() void FixQEqReax::compute_H() { int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh; - int i, j, ii, jj, temp, newnbr, flag; + int i, j, ii, jj, flag; double **x, SMALL = 0.0001; double dx, dy, dz, r_sqr; @@ -651,7 +651,7 @@ int FixQEqReax::CG( double *b, double *x ) { int i, j, imax; double tmp, alpha, beta, b_norm; - double sig_old, sig_new, sig0; + double sig_old, sig_new; int nn, jj; int *ilist; @@ -677,17 +677,14 @@ int FixQEqReax::CG( double *b, double *x ) d[j] = r[j] * Hdia_inv[j]; //pre-condition } - int ttype = 1; b_norm = parallel_norm( b, nn ); sig_new = parallel_dot( r, d, nn); - sig0 = sig_new; for( i = 1; i < imax && sqrt(sig_new) / b_norm > tolerance; ++i ) { comm->forward_comm_fix(this); //Dist_vector( d ); sparse_matvec( &H, d, q ); comm->reverse_comm_fix(this); //Coll_vector( q ); - ttype = 2; tmp = parallel_dot( d, q, nn); alpha = sig_new / tmp; diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp index 47c9ea24a6..cb8fa0be57 100644 --- a/src/USER-REAXC/pair_reax_c.cpp +++ b/src/USER-REAXC/pair_reax_c.cpp @@ -630,37 +630,28 @@ void PairReaxC::set_far_nbr( far_neighbor_data *fdest, int PairReaxC::estimate_reax_lists() { - int itr_i, itr_j, itr_g, i, j, g; - int nlocal, nghost, num_nbrs, num_marked; + int itr_i, itr_j, i, j; + int num_nbrs, num_marked; int *ilist, *jlist, *numneigh, **firstneigh, *marked; - double d_sqr, g_d_sqr; - rvec dvec, g_dvec; + double d_sqr; + rvec dvec; double **x; - reax_list *far_nbrs; - far_neighbor_data *far_list; int mincap = system->mincap; double safezone = system->safezone; x = atom->x; - nlocal = atom->nlocal; - nghost = atom->nghost; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - far_nbrs = lists + FAR_NBRS; - far_list = far_nbrs->select.far_nbr_list; - num_nbrs = 0; num_marked = 0; marked = (int*) calloc( system->N, sizeof(int) ); - int inum = list->inum; - int gnum = list->gnum; - int numall = inum + gnum; + int numall = list->inum + list->gnum; - for( itr_i = 0; itr_i < inum+gnum; ++itr_i ){ + for( itr_i = 0; itr_i < numall; ++itr_i ){ i = ilist[itr_i]; marked[i] = 1; ++num_marked; @@ -685,18 +676,16 @@ int PairReaxC::estimate_reax_lists() int PairReaxC::write_reax_lists() { - int itr_i, itr_j, itr_g, i, j, g, flag; - int nlocal, nghost, num_nbrs; + int itr_i, itr_j, i, j; + int num_nbrs; int *ilist, *jlist, *numneigh, **firstneigh; - double d_sqr, g_d, g_d_sqr; - rvec dvec, g_dvec; - double *dist, **x, SMALL = 0.0001; + double d_sqr; + rvec dvec; + double *dist, **x; reax_list *far_nbrs; far_neighbor_data *far_list; x = atom->x; - nlocal = atom->nlocal; - nghost = atom->nghost; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; @@ -707,11 +696,9 @@ int PairReaxC::write_reax_lists() num_nbrs = 0; dist = (double*) calloc( system->N, sizeof(double) ); - int inum = list->inum; - int gnum = list->gnum; - int numall = inum + gnum; + int numall = list->inum + list->gnum; - for( itr_i = 0; itr_i < inum+gnum; ++itr_i ){ + for( itr_i = 0; itr_i < numall; ++itr_i ){ i = ilist[itr_i]; jlist = firstneigh[i]; Set_Start_Index( i, num_nbrs, far_nbrs ); @@ -792,8 +779,6 @@ double PairReaxC::memory_usage() bytes += 19.0 * system->total_cap * sizeof(real); bytes += 3.0 * system->total_cap * sizeof(int); - double mem1 = bytes; - // From reaxc_lists bytes += 2.0 * lists->n * sizeof(int); bytes += lists->num_intrs * sizeof(three_body_interaction_data); @@ -813,8 +798,8 @@ double PairReaxC::memory_usage() void PairReaxC::FindBond() { - int i, ii, j, pj, nj, jtmp, jj; - double bo_tmp, bo_cut, rij, rsq, r_tmp; + int i, j, pj, nj; + double bo_tmp, bo_cut; bond_data *bo_ij; bo_cut = 0.10; @@ -827,7 +812,6 @@ void PairReaxC::FindBond() if (j < i) continue; bo_tmp = bo_ij->bo_data.BO; - r_tmp = bo_ij->d; if (bo_tmp >= bo_cut ) { tmpid[i][nj] = j; diff --git a/src/USER-REAXC/reaxc_allocate.cpp b/src/USER-REAXC/reaxc_allocate.cpp index 1bf3d85750..474b60d060 100644 --- a/src/USER-REAXC/reaxc_allocate.cpp +++ b/src/USER-REAXC/reaxc_allocate.cpp @@ -54,31 +54,6 @@ int PreAllocate_Space( reax_system *system, control_params *control, /************* system *************/ -inline void reax_atom_Copy( reax_atom *dest, reax_atom *src ) -{ - dest->orig_id = src->orig_id; - dest->type = src->type; - strcpy( dest->name, src->name ); - rvec_Copy( dest->x, src->x ); - rvec_Copy( dest->v, src->v ); - rvec_Copy( dest->f_old, src->f_old ); - rvec_Copy( dest->s, src->s ); - rvec_Copy( dest->t, src->t ); - dest->Hindex = src->Hindex; - dest->num_bonds = src->num_bonds; - dest->num_hbonds = src->num_hbonds; - dest->numbonds = src->numbonds; -} - - -void Copy_Atom_List( reax_atom *dest, reax_atom *src, int n ) -{ - int i; - - for( i = 0; i < n; ++i ) - memcpy( dest+i, src+i, sizeof(reax_atom) ); -} - int Allocate_System( reax_system *system, int local_cap, int total_cap, char *msg ) @@ -306,8 +281,8 @@ int Allocate_Workspace( reax_system *system, control_params *control, } -void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs, - MPI_Comm comm ) +static void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, + int num_intrs, MPI_Comm comm ) { Delete_List( far_nbrs, comm ); if(!Make_List( n, num_intrs, TYP_FAR_NEIGHBOR, far_nbrs, comm )){ @@ -317,48 +292,8 @@ void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs, } -int Allocate_Matrix( sparse_matrix **pH, int cap, int m, MPI_Comm comm ) -{ - sparse_matrix *H; - - *pH = (sparse_matrix*) - smalloc( sizeof(sparse_matrix), "sparse_matrix", comm ); - H = *pH; - H->cap = cap; - H->m = m; - H->start = (int*) smalloc( sizeof(int) * cap, "matrix_start", comm ); - H->end = (int*) smalloc( sizeof(int) * cap, "matrix_end", comm ); - H->entries = (sparse_matrix_entry*) - smalloc( sizeof(sparse_matrix_entry)*m, "matrix_entries", comm ); - - return SUCCESS; -} - - -void Deallocate_Matrix( sparse_matrix *H ) -{ - sfree(H->start, "H->start"); - sfree(H->end, "H->end"); - sfree(H->entries, "H->entries"); - sfree(H, "H"); -} - - -int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name, - MPI_Comm comm ) -{ - Deallocate_Matrix( *H ); - if( !Allocate_Matrix( H, n, m, comm ) ) { - fprintf(stderr, "not enough space for %s matrix. terminating!\n", name); - MPI_Abort( comm, INSUFFICIENT_MEMORY ); - } - - return SUCCESS; -} - - -int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds, - MPI_Comm comm ) +static int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds, + MPI_Comm comm ) { int i, id, total_hbonds; @@ -382,8 +317,9 @@ int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds, } -int Reallocate_Bonds_List( reax_system *system, reax_list *bonds, - int *total_bonds, int *est_3body, MPI_Comm comm ) +static int Reallocate_Bonds_List( reax_system *system, reax_list *bonds, + int *total_bonds, int *est_3body, + MPI_Comm comm ) { int i; @@ -408,174 +344,6 @@ int Reallocate_Bonds_List( reax_system *system, reax_list *bonds, } -/************* grid *************/ -int Estimate_GCell_Population( reax_system* system, MPI_Comm comm ) -{ - int d, i, j, k, l, max_atoms, my_max, all_max; - ivec c; - grid *g; - grid_cell *gc; - simulation_box *my_ext_box; - reax_atom *atoms; - - my_ext_box = &(system->my_ext_box); - g = &(system->my_grid); - atoms = system->my_atoms; - Reset_Grid( g ); - - for( l = 0; l < system->n; l++ ) { - for( d = 0; d < 3; ++d ) { - - c[d] = (int)((atoms[l].x[d]-my_ext_box->min[d])*g->inv_len[d]); - - if( c[d] >= g->native_end[d] ) - c[d] = g->native_end[d] - 1; - else if( c[d] < g->native_str[d] ) - c[d] = g->native_str[d]; - } - gc = &( g->cells[c[0]][c[1]][c[2]] ); - gc->top++; - } - - max_atoms = 0; - for( i = 0; i < g->ncells[0]; i++ ) - for( j = 0; j < g->ncells[1]; j++ ) - for( k = 0; k < g->ncells[2]; k++ ) { - gc = &(g->cells[i][j][k]); - if( max_atoms < gc->top ) - max_atoms = gc->top; - } - - my_max = (int)(MAX(max_atoms*SAFE_ZONE, MIN_GCELL_POPL)); - MPI_Allreduce( &my_max, &all_max, 1, MPI_INT, MPI_MAX, comm ); - - return all_max; -} - - -void Allocate_Grid( reax_system *system, MPI_Comm comm ) -{ - int i, j, k, l; - grid *g; - grid_cell *gc; - - g = &( system->my_grid ); - - /* allocate gcell reordering space */ - g->order = (ivec*) scalloc( g->total+1, sizeof(ivec), "g:order", comm ); - - /* allocate the gcells for the new grid */ - g->max_nbrs = (2*g->vlist_span[0]+1)*(2*g->vlist_span[1]+1)* - (2*g->vlist_span[2]+1)+3; - - g->cells = (grid_cell***) - scalloc( g->ncells[0], sizeof(grid_cell**), "gcells", comm ); - for( i = 0; i < g->ncells[0]; i++ ) { - g->cells[i] = (grid_cell**) - scalloc( g->ncells[1], sizeof(grid_cell*),"gcells[i]", comm ); - - for( j = 0; j < g->ncells[1]; ++j ) { - g->cells[i][j] = (grid_cell*) - scalloc( g->ncells[2], sizeof(grid_cell), "gcells[i][j]", comm ); - - for( k = 0; k < g->ncells[2]; k++ ) { - gc = &(g->cells[i][j][k]); - gc->top = gc->mark = gc->str = gc->end = 0; - gc->nbrs = (grid_cell**) - scalloc( g->max_nbrs, sizeof(grid_cell*), "g:nbrs", comm ); - gc->nbrs_x = (ivec*) - scalloc( g->max_nbrs, sizeof(ivec), "g:nbrs_x", comm ); - gc->nbrs_cp = (rvec*) - scalloc( g->max_nbrs, sizeof(rvec), "g:nbrs_cp", comm ); - for( l = 0; l < g->max_nbrs; ++l ) - gc->nbrs[l] = NULL; - } - } - } - - /* allocate atom id storage in gcells */ - g->max_atoms = Estimate_GCell_Population( system, comm ); - /* space for storing atom id's is required only for native cells */ - for( i = g->native_str[0]; i < g->native_end[0]; ++i ) - for( j = g->native_str[1]; j < g->native_end[1]; ++j ) - for( k = g->native_str[2]; k < g->native_end[2]; ++k ) - g->cells[i][j][k].atoms = (int*) scalloc( g->max_atoms, sizeof(int), - "g:atoms", comm ); -} - - -void Deallocate_Grid( grid *g ) -{ - int i, j, k; - grid_cell *gc; - - sfree( g->order, "g:order" ); - - /* deallocate the grid cells */ - for( i = 0; i < g->ncells[0]; i++ ) { - for( j = 0; j < g->ncells[1]; j++ ) { - for( k = 0; k < g->ncells[2]; k++ ) { - gc = &(g->cells[i][j][k]); - sfree( gc->nbrs, "g:nbrs" ); - sfree( gc->nbrs_x, "g:nbrs_x" ); - sfree( gc->nbrs_cp, "g:nbrs_cp" ); - if(gc->atoms != NULL ) - sfree( gc->atoms, "g:atoms" ); - } - sfree( g->cells[i][j], "g:cells[i][j]" ); - } - sfree( g->cells[i], "g:cells[i]" ); - } - sfree( g->cells, "g:cells" ); -} - - -int Allocate_MPI_Buffers( mpi_datatypes *mpi_data, int est_recv, - neighbor_proc *my_nbrs, char *msg ) -{ - int i; - mpi_out_data *mpi_buf; - MPI_Comm comm; - - comm = mpi_data->world; - - /* in buffers */ - mpi_data->in1_buffer = (void*) - scalloc( est_recv, sizeof(boundary_atom), "in1_buffer", comm ); - mpi_data->in2_buffer = (void*) - scalloc( est_recv, sizeof(boundary_atom), "in2_buffer", comm ); - - /* out buffers */ - for( i = 0; i < MAX_NBRS; ++i ) { - mpi_buf = &( mpi_data->out_buffers[i] ); - /* allocate storage for the neighbor processor i */ - mpi_buf->index = (int*) - scalloc( my_nbrs[i].est_send, sizeof(int), "mpibuf:index", comm ); - mpi_buf->out_atoms = (void*) - scalloc( my_nbrs[i].est_send, sizeof(boundary_atom), "mpibuf:out_atoms", - comm ); - } - - return SUCCESS; -} - - -void Deallocate_MPI_Buffers( mpi_datatypes *mpi_data ) -{ - int i; - mpi_out_data *mpi_buf; - - sfree( mpi_data->in1_buffer, "in1_buffer" ); - sfree( mpi_data->in2_buffer, "in2_buffer" ); - - for( i = 0; i < MAX_NBRS; ++i ) { - mpi_buf = &( mpi_data->out_buffers[i] ); - sfree( mpi_buf->index, "mpibuf:index" ); - sfree( mpi_buf->out_atoms, "mpibuf:out_atoms" ); - } -} - - void ReAllocate( reax_system *system, control_params *control, simulation_data *data, storage *workspace, reax_list **lists, mpi_datatypes *mpi_data ) @@ -584,7 +352,6 @@ void ReAllocate( reax_system *system, control_params *control, int renbr, newsize; reallocate_data *realloc; reax_list *far_nbrs; - grid *g; MPI_Comm comm; char msg[200]; @@ -593,13 +360,10 @@ void ReAllocate( reax_system *system, control_params *control, double saferzone = system->saferzone; realloc = &(workspace->realloc); - g = &(system->my_grid); comm = mpi_data->world; - int nflag = 0; if( system->n >= DANGER_ZONE * system->local_cap || (0 && system->n <= LOOSE_ZONE * system->local_cap) ) { - nflag = 1; system->local_cap = MAX( (int)(system->n * safezone), mincap ); } diff --git a/src/USER-REAXC/reaxc_allocate.h b/src/USER-REAXC/reaxc_allocate.h index 0ad5f1a710..f8814859af 100644 --- a/src/USER-REAXC/reaxc_allocate.h +++ b/src/USER-REAXC/reaxc_allocate.h @@ -28,9 +28,8 @@ #define __ALLOCATE_H_ #include "reaxc_types.h" -int PreAllocate_Space( reax_system*, control_params*, storage*, MPI_Comm ); +int PreAllocate_Space( reax_system*, control_params*, storage*, MPI_Comm ); -void reax_atom_Copy( reax_atom*, reax_atom* ); int Allocate_System( reax_system*, int, int, char* ); void DeAllocate_System( reax_system* ); @@ -38,17 +37,6 @@ int Allocate_Workspace( reax_system*, control_params*, storage*, int, int, MPI_Comm, char* ); void DeAllocate_Workspace( control_params*, storage* ); -void Allocate_Grid( reax_system*, MPI_Comm ); -void Deallocate_Grid( grid* ); - -int Allocate_MPI_Buffers( mpi_datatypes*, int, neighbor_proc*, char* ); - -int Allocate_Matrix( sparse_matrix**, int, int, MPI_Comm ); - -int Allocate_HBond_List( int, int, int*, int*, reax_list* ); - -int Allocate_Bond_List( int, int*, reax_list* ); - void ReAllocate( reax_system*, control_params*, simulation_data*, storage*, reax_list**, mpi_datatypes* ); #endif diff --git a/src/USER-REAXC/reaxc_bond_orders.cpp b/src/USER-REAXC/reaxc_bond_orders.cpp index 11f862f64b..7103f48299 100644 --- a/src/USER-REAXC/reaxc_bond_orders.cpp +++ b/src/USER-REAXC/reaxc_bond_orders.cpp @@ -150,7 +150,6 @@ void Add_dBond_to_Forces( reax_system *system, int i, int pj, int pk, k, j; /* Virial Tallying variables */ - real f_scaler; rvec fi_tmp, fj_tmp, fk_tmp, delij, delji, delki, delkj, temp; /* Initializations */ @@ -370,7 +369,7 @@ void BO( reax_system *system, control_params *control, simulation_data *data, storage *workspace, reax_list **lists, output_controls *out_control ) { int i, j, pj, type_i, type_j; - int start_i, end_i, sym_index, num_bonds; + 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; @@ -384,7 +383,6 @@ void BO( reax_system *system, control_params *control, simulation_data *data, bond_order_data *bo_ij, *bo_ji; reax_list *bonds = (*lists) + BONDS; - num_bonds = 0; p_boc1 = system->reax_param.gp.l[0]; p_boc2 = system->reax_param.gp.l[1]; diff --git a/src/USER-REAXC/reaxc_control.cpp b/src/USER-REAXC/reaxc_control.cpp index eadce53a0c..9d067c16a5 100644 --- a/src/USER-REAXC/reaxc_control.cpp +++ b/src/USER-REAXC/reaxc_control.cpp @@ -33,7 +33,7 @@ char Read_Control_File( char *control_file, control_params* control, { FILE *fp; char *s, **tmp; - int c,i,ival; + int i,ival; real val; /* open control file */ @@ -115,7 +115,7 @@ char Read_Control_File( char *control_file, control_params* control, /* read control parameters file */ while (!feof(fp)) { fgets( s, MAX_LINE, fp ); - c = Tokenize( s, &tmp ); + Tokenize( s, &tmp ); if( strcmp(tmp[0], "simulation_name") == 0 ) { strcpy( control->sim_name, tmp[1] ); diff --git a/src/USER-REAXC/reaxc_forces.cpp b/src/USER-REAXC/reaxc_forces.cpp index d581cbebbd..b9fd78dd76 100644 --- a/src/USER-REAXC/reaxc_forces.cpp +++ b/src/USER-REAXC/reaxc_forces.cpp @@ -28,7 +28,6 @@ #include "reaxc_forces.h" #include "reaxc_bond_orders.h" #include "reaxc_bonds.h" -#include "reaxc_basic_comm.h" #include "reaxc_hydrogen_bonds.h" #include "reaxc_io_tools.h" #include "reaxc_list.h" @@ -120,8 +119,6 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists, { int i, comp, Hindex; reax_list *bonds, *hbonds; - reallocate_data *realloc; - realloc = &(workspace->realloc); double saferzone = system->saferzone; @@ -391,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 r_ij, cutoff; + real cutoff; reax_list *far_nbrs, *bonds, *hbonds; single_body_parameters *sbp_i, *sbp_j; two_body_parameters *twbp; @@ -468,7 +465,6 @@ void Init_Forces_noQEq( reax_system *system, control_params *control, if( flag ) { type_j = atom_j->type; if (type_j < 0) continue; - r_ij = nbr_pj->d; sbp_j = &(system->reax_param.sbp[type_j]); twbp = &(system->reax_param.tbp[type_i][type_j]); @@ -536,7 +532,7 @@ void Estimate_Storages( reax_system *system, control_params *control, int ihb, jhb; int local; real cutoff; - real r_ij, r2; + real r_ij; real C12, C34, C56; real BO, BO_s, BO_pi, BO_pi2; reax_list *far_nbrs; @@ -604,8 +600,6 @@ void Estimate_Storages( reax_system *system, control_params *control, /* uncorrected bond orders */ if( nbr_pj->d <= control->bond_cut ) { - r2 = SQR(r_ij); - if( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0) { C12 = twbp->p_bo1 * pow( r_ij / twbp->r_s, twbp->p_bo2 ); BO_s = (1.0 + control->bo_cut) * exp( C12 ); diff --git a/src/USER-REAXC/reaxc_hydrogen_bonds.cpp b/src/USER-REAXC/reaxc_hydrogen_bonds.cpp index b366872b7a..b51ac6a183 100644 --- a/src/USER-REAXC/reaxc_hydrogen_bonds.cpp +++ b/src/USER-REAXC/reaxc_hydrogen_bonds.cpp @@ -42,7 +42,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control, int itr, top; int num_hb_intrs = 0; ivec rel_jk; - real r_ij, r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; + real r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; real 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; @@ -102,7 +102,6 @@ void Hydrogen_Bonds( reax_system *system, control_params *control, bo_ij = &(pbond_ij->bo_data); type_i = system->my_atoms[i].type; if (type_i < 0) continue; - r_ij = pbond_ij->d; hbp = &(system->reax_param.hbp[ type_i ][ type_j ][ type_k ]); ++num_hb_intrs; diff --git a/src/USER-REAXC/reaxc_init_md.cpp b/src/USER-REAXC/reaxc_init_md.cpp index df1a857028..a674b2e39f 100644 --- a/src/USER-REAXC/reaxc_init_md.cpp +++ b/src/USER-REAXC/reaxc_init_md.cpp @@ -155,10 +155,8 @@ int Init_Lists( reax_system *system, control_params *control, simulation_data *data, storage *workspace, reax_list **lists, mpi_datatypes *mpi_data, char *msg ) { - int i, num_nbrs; - int total_hbonds, total_bonds, bond_cap, num_3body, cap_3body, Htop; + int i, total_hbonds, total_bonds, bond_cap, num_3body, cap_3body, Htop; int *hb_top, *bond_top; - int nrecv[MAX_NBRS]; MPI_Comm comm; int mincap = system->mincap; diff --git a/src/USER-REAXC/reaxc_io_tools.cpp b/src/USER-REAXC/reaxc_io_tools.cpp index 9645234bd4..513b59c9f6 100644 --- a/src/USER-REAXC/reaxc_io_tools.cpp +++ b/src/USER-REAXC/reaxc_io_tools.cpp @@ -27,7 +27,6 @@ #include "pair_reax_c.h" #include "update.h" #include "reaxc_io_tools.h" -#include "reaxc_basic_comm.h" #include "reaxc_list.h" #include "reaxc_reset_tools.h" #include "reaxc_system_props.h" diff --git a/src/USER-REAXC/reaxc_multi_body.cpp b/src/USER-REAXC/reaxc_multi_body.cpp index 220fcb072c..f1b56ef5de 100644 --- a/src/USER-REAXC/reaxc_multi_body.cpp +++ b/src/USER-REAXC/reaxc_multi_body.cpp @@ -43,19 +43,18 @@ void Atom_Energy( reax_system *system, control_params *control, 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_lp1, p_lp2, p_lp3; + real p_lp2, p_lp3; real p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8; - real eng_tmp, f_tmp; + real eng_tmp; int numbonds; - single_body_parameters *sbp_i, *sbp_j; + single_body_parameters *sbp_i; two_body_parameters *twbp; bond_data *pbond; bond_order_data *bo_ij; reax_list *bonds = (*lists) + BONDS; /* Initialize parameters */ - p_lp1 = system->reax_param.gp.l[15]; p_lp3 = system->reax_param.gp.l[5]; p_ovun3 = system->reax_param.gp.l[32]; p_ovun4 = system->reax_param.gp.l[31]; @@ -143,7 +142,6 @@ void Atom_Energy( reax_system *system, control_params *control, type_j = system->my_atoms[j].type; if (type_j < 0) continue; bo_ij = &(bonds->select.bond_list[pj].bo_data); - sbp_j = &(system->reax_param.sbp[ type_j ]); twbp = &(system->reax_param.tbp[ type_i ][ type_j ]); sum_ovun1 += twbp->p_ovun1 * twbp->De_s * bo_ij->BO; @@ -205,7 +203,6 @@ void Atom_Energy( reax_system *system, control_params *control, if( system->pair_ptr->evflag) { eng_tmp = e_ov; if (numbonds > 0) eng_tmp += e_un; - f_tmp = CEover3 + CEunder3; system->pair_ptr->ev_tally(i,i,system->n,1,eng_tmp,0.0,0.0,0.0,0.0,0.0); } diff --git a/src/USER-REAXC/reaxc_tool_box.cpp b/src/USER-REAXC/reaxc_tool_box.cpp index 55feac209f..7fccbb3427 100644 --- a/src/USER-REAXC/reaxc_tool_box.cpp +++ b/src/USER-REAXC/reaxc_tool_box.cpp @@ -78,46 +78,6 @@ void Fit_to_Periodic_Box( simulation_box *box, rvec *p ) } } -/************** from geo_tools.c *****************/ -void Make_Point( real x, real y, real z, rvec* p ) -{ - (*p)[0] = x; - (*p)[1] = y; - (*p)[2] = z; -} - - - -int is_Valid_Serial( storage *workspace, int serial ) -{ - return SUCCESS; -} - - - -int Check_Input_Range( int val, int lo, int hi, char *message, MPI_Comm comm ) -{ - if( val < lo || val > hi ) { - fprintf( stderr, "%s\nInput %d - Out of range %d-%d. Terminating...\n", - message, val, lo, hi ); - MPI_Abort( comm, INVALID_INPUT ); - } - - return 1; -} - - -void Trim_Spaces( char *element ) -{ - int i, j; - - for( i = 0; element[i] == ' '; ++i ); // skip initial space chars - - for( j = i; j < (int)(strlen(element)) && element[j] != ' '; ++j ) - element[j-i] = toupper( element[j] ); // make uppercase, offset to 0 - element[j-i] = 0; // finalize the string -} - struct timeval tim; real t_end; diff --git a/src/USER-REAXC/reaxc_tool_box.h b/src/USER-REAXC/reaxc_tool_box.h index 2ef27dd6da..d7784fc0d3 100644 --- a/src/USER-REAXC/reaxc_tool_box.h +++ b/src/USER-REAXC/reaxc_tool_box.h @@ -48,12 +48,6 @@ real DistSqr_between_Special_Points( rvec, rvec ); real DistSqr_to_Special_Point( rvec, rvec ); int Relative_Coord_Encoding( ivec ); -/* from geo_tools.h */ -void Make_Point( real, real, real, rvec* ); -int is_Valid_Serial( storage*, int ); -int Check_Input_Range( int, int, int, char*, MPI_Comm ); -void Trim_Spaces( char* ); - /* from system_props.h */ real Get_Time( ); real Get_Timing_Info( real ); diff --git a/src/USER-REAXC/reaxc_torsion_angles.cpp b/src/USER-REAXC/reaxc_torsion_angles.cpp index 16f6a4539f..223fd1d40c 100644 --- a/src/USER-REAXC/reaxc_torsion_angles.cpp +++ b/src/USER-REAXC/reaxc_torsion_angles.cpp @@ -124,7 +124,7 @@ void Torsion_Angles( reax_system *system, control_params *control, { int i, j, k, l, pi, pj, pk, pl, pij, plk, natoms; int type_i, type_j, type_k, type_l; - int start_j, end_j, start_k, end_k; + int start_j, end_j; int start_pj, end_pj, start_pk, end_pk; int num_frb_intrs = 0; @@ -165,7 +165,7 @@ void Torsion_Angles( reax_system *system, control_params *control, // Virial tallying variables real delil[3], deljl[3], delkl[3]; - real eng_tmp, f_scaler, fi_tmp[3], fj_tmp[3], fk_tmp[3]; + real eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; natoms = system->n; @@ -183,8 +183,6 @@ void Torsion_Angles( reax_system *system, control_params *control, if( system->my_atoms[j].orig_id < system->my_atoms[k].orig_id && bo_jk->BO > control->thb_cut/*0*/ && Num_Entries(pk, thb_intrs) ) { - start_k = Start_Index(k, bonds); - end_k = End_Index(k, bonds); pj = pbond_jk->sym_index; // pj points to j on k's list if( Num_Entries(pj, thb_intrs) ) { diff --git a/src/USER-REAXC/reaxc_traj.cpp b/src/USER-REAXC/reaxc_traj.cpp index 3a94ed45ef..05bf73a29f 100644 --- a/src/USER-REAXC/reaxc_traj.cpp +++ b/src/USER-REAXC/reaxc_traj.cpp @@ -54,7 +54,6 @@ void Write_Skip_Line( output_controls *out_control, mpi_datatypes *mpi_data, if( my_rank == MASTER_NODE ) fprintf( out_control->strj, INT2_LINE, "chars_to_skip_section:", skip, num_section ); - } diff --git a/src/USER-REAXC/reaxc_types.h b/src/USER-REAXC/reaxc_types.h index 816906ff3e..b05e86f16f 100644 --- a/src/USER-REAXC/reaxc_types.h +++ b/src/USER-REAXC/reaxc_types.h @@ -68,45 +68,10 @@ typedef real rtensor[3][3]; typedef real rvec2[2]; typedef real rvec4[4]; + +// import LAMMPS' definition of tagint typedef LAMMPS_NS::tagint rc_tagint; -typedef struct { - int step, bigN; - real T, xi, v_xi, v_xi_old, G_xi; - rtensor box; -} restart_header; - -typedef struct { - rc_tagint orig_id, type; - char name[8]; - rvec x, v; -} restart_atom; - -typedef struct -{ - rc_tagint orig_id; - int imprt_id; - int type; - int num_bonds; - int num_hbonds; - char name[8]; - rvec x; // position - rvec v; // velocity - rvec f_old; // old force - rvec4 s, t; // for calculating q -} mpi_atom; - -typedef struct -{ - rc_tagint orig_id; - int imprt_id; - int type; - int num_bonds; - int num_hbonds; - rvec x; // position -} boundary_atom; - - typedef struct { int cnt; diff --git a/src/USER-REAXC/reaxc_valence_angles.cpp b/src/USER-REAXC/reaxc_valence_angles.cpp index 532d9a6a7e..60f985598d 100644 --- a/src/USER-REAXC/reaxc_valence_angles.cpp +++ b/src/USER-REAXC/reaxc_valence_angles.cpp @@ -30,6 +30,16 @@ #include "reaxc_list.h" #include "reaxc_vector.h" +static real Dot( real* v1, real* v2, int k ) +{ + real ret = 0.0; + + for( int i=0; i < k; ++i ) + ret += v1[i] * v2[i]; + + return ret; +} + void Calculate_Theta( rvec dvec_ji, real d_ji, rvec dvec_jk, real d_jk, real *theta, real *cos_theta ) { @@ -88,12 +98,11 @@ void Valence_Angles( reax_system *system, control_params *control, 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 r_ij, r_jk; real BOA_ij, BOA_jk; rvec force, ext_press; // Tallying variables - real eng_tmp, f_scaler, fi_tmp[3], fj_tmp[3], fk_tmp[3]; + real eng_tmp, fi_tmp[3], fj_tmp[3], fk_tmp[3]; real delij[3], delkj[3]; three_body_header *thbh; @@ -168,7 +177,6 @@ void Valence_Angles( reax_system *system, control_params *control, if( BOA_ij/*bo_ij->BO*/ > 0.0 && ( j < system->n || pbond_ij->nbr < system->n ) ) { i = pbond_ij->nbr; - r_ij = pbond_ij->d; type_i = system->my_atoms[i].type; for( pk = start_j; pk < pi; ++pk ) { @@ -223,7 +231,6 @@ void Valence_Angles( reax_system *system, control_params *control, (bo_ij->BO > control->thb_cut) && (bo_jk->BO > control->thb_cut) && (bo_ij->BO * bo_jk->BO > control->thb_cutsq) ) { - r_jk = pbond_jk->d; thbh = &( system->reax_param.thbp[ type_i ][ type_j ][ type_k ] ); for( cnt = 0; cnt < thbh->cnt; ++cnt ) { diff --git a/src/USER-REAXC/reaxc_vector.cpp b/src/USER-REAXC/reaxc_vector.cpp index 4d56442172..2565064217 100644 --- a/src/USER-REAXC/reaxc_vector.cpp +++ b/src/USER-REAXC/reaxc_vector.cpp @@ -27,89 +27,13 @@ #include "pair_reax_c.h" #include "reaxc_vector.h" -int Vector_isZero( real* v, int k ) -{ - for( --k; k>=0; --k ) - if( fabs( v[k] ) > ALMOST_ZERO ) - return 0; - - return 1; -} - - -void Vector_MakeZero( real *v, int k ) -{ - for( --k; k>=0; --k ) - v[k] = 0; -} - - -void Vector_Copy( real* dest, real* v, int k ) -{ - for( --k; k>=0; --k ) - dest[k] = v[k]; -} - - -void Vector_Scale( real* dest, real c, real* v, int k ) -{ - for( --k; k>=0; --k ) - dest[k] = c * v[k]; -} - - -void Vector_Sum( real* dest, real c, real* v, real d, real* y, int k ) -{ - for( --k; k>=0; --k ) - dest[k] = c * v[k] + d * y[k]; -} - - -void Vector_Add( real* dest, real c, real* v, int k ) -{ - for( --k; k>=0; --k ) - dest[k] += c * v[k]; -} - - -real Dot( real* v1, real* v2, int k ) -{ - real ret = 0; - - for( --k; k>=0; --k ) - ret += v1[k] * v2[k]; - - return ret; -} - - -real Norm( real* v1, int k ) -{ - real ret = 0; - - for( --k; k>=0; --k ) - ret += SQR( v1[k] ); - - return sqrt( ret ); -} - - -void Vector_Print( FILE *fout, char *vname, real *v, int k ) -{ - int i; - - fprintf( fout, "%s:", vname ); - for( i = 0; i < k; ++i ) - fprintf( fout, "%24.15e\n", v[i] ); - fprintf( fout, "\n" ); -} - void rvec_Copy( rvec dest, rvec src ) { dest[0] = src[0], dest[1] = src[1], dest[2] = src[2]; } + void rvec_Scale( rvec ret, real c, rvec v ) { ret[0] = c * v[0], ret[1] = c * v[1], ret[2] = c * v[2]; @@ -241,37 +165,6 @@ void rvec_MakeZero( rvec v ) v[0] = v[1] = v[2] = 0.000000000000000e+00; } -void rtensor_Multiply( rtensor ret, rtensor m1, rtensor m2 ) -{ - int i, j, k; - rtensor temp; - - if( ret == m1 || ret == m2 ) - { - for( i = 0; i < 3; ++i ) - for( j = 0; j < 3; ++j ) - { - temp[i][j] = 0; - for( k = 0; k < 3; ++k ) - temp[i][j] += m1[i][k] * m2[k][j]; - } - - for( i = 0; i < 3; ++i ) - for( j = 0; j < 3; ++j ) - ret[i][j] = temp[i][j]; - } - else - { - for( i = 0; i < 3; ++i ) - for( j = 0; j < 3; ++j ) - { - ret[i][j] = 0; - for( k = 0; k < 3; ++k ) - ret[i][j] += m1[i][k] * m2[k][j]; - } - } -} - void rtensor_MatVec( rvec ret, rtensor m, rvec v ) { @@ -304,64 +197,6 @@ void rtensor_Scale( rtensor ret, real c, rtensor m ) } -void rtensor_Add( rtensor ret, rtensor t ) -{ - int i, j; - - for( i = 0; i < 3; ++i ) - for( j = 0; j < 3; ++j ) - ret[i][j] += t[i][j]; -} - - -void rtensor_ScaledAdd( rtensor ret, real c, rtensor t ) -{ - int i, j; - - for( i = 0; i < 3; ++i ) - for( j = 0; j < 3; ++j ) - ret[i][j] += c * t[i][j]; -} - - -void rtensor_Sum( rtensor ret, rtensor t1, rtensor t2 ) -{ - int i, j; - - for( i = 0; i < 3; ++i ) - for( j = 0; j < 3; ++j ) - ret[i][j] = t1[i][j] + t2[i][j]; -} - - -void rtensor_ScaledSum( rtensor ret, real c1, rtensor t1, - real c2, rtensor t2 ) -{ - int i, j; - - for( i = 0; i < 3; ++i ) - for( j = 0; j < 3; ++j ) - ret[i][j] = c1 * t1[i][j] + c2 * t2[i][j]; -} - - -void rtensor_Copy( rtensor ret, rtensor t ) -{ - int i, j; - - for( i = 0; i < 3; ++i ) - for( j = 0; j < 3; ++j ) - ret[i][j] = t[i][j]; -} - - -void rtensor_Identity( rtensor t ) -{ - t[0][0] = t[1][1] = t[2][2] = 1; - t[0][1] = t[0][2] = t[1][0] = t[1][2] = t[2][0] = t[2][1] = 0; -} - - void rtensor_MakeZero( rtensor t ) { t[0][0] = t[0][1] = t[0][2] = 0; @@ -370,43 +205,6 @@ void rtensor_MakeZero( rtensor t ) } -void rtensor_Transpose( rtensor ret, rtensor t ) -{ - ret[0][0] = t[0][0], ret[1][1] = t[1][1], ret[2][2] = t[2][2]; - ret[0][1] = t[1][0], ret[0][2] = t[2][0]; - ret[1][0] = t[0][1], ret[1][2] = t[2][1]; - ret[2][0] = t[0][2], ret[2][1] = t[1][2]; -} - - -real rtensor_Det( rtensor t ) -{ - return ( t[0][0] * (t[1][1] * t[2][2] - t[1][2] * t[2][1] ) + - t[0][1] * (t[1][2] * t[2][0] - t[1][0] * t[2][2] ) + - t[0][2] * (t[1][0] * t[2][1] - t[1][1] * t[2][0] ) ); -} - - -real rtensor_Trace( rtensor t ) -{ - return (t[0][0] + t[1][1] + t[2][2]); -} - - -void Print_rTensor(FILE* fp, rtensor t) -{ - int i, j; - - for (i=0; i < 3; i++) - { - fprintf(fp,"["); - for (j=0; j < 3; j++) - fprintf(fp,"%8.3f,\t",t[i][j]); - fprintf(fp,"]\n"); - } -} - - void ivec_MakeZero( ivec v ) { v[0] = v[1] = v[2] = 0; @@ -427,30 +225,6 @@ void ivec_Scale( ivec dest, real C, ivec src ) } -void ivec_rScale( ivec dest, real C, rvec src ) -{ - dest[0] = (int)(C * src[0]); - dest[1] = (int)(C * src[1]); - dest[2] = (int)(C * src[2]); -} - - -int ivec_isZero( ivec v ) -{ - if( v[0]==0 && v[1]==0 && v[2]==0 ) - return 1; - return 0; -} - - -int ivec_isEqual( ivec v1, ivec v2 ) -{ - if( v1[0]==v2[0] && v1[1]==v2[1] && v1[2]==v2[2] ) - return 1; - return 0; -} - - void ivec_Sum( ivec dest, ivec v1, ivec v2 ) { dest[0] = v1[0] + v2[0]; @@ -459,42 +233,3 @@ void ivec_Sum( ivec dest, ivec v1, ivec v2 ) } -void ivec_ScaledSum( ivec dest, int k1, ivec v1, int k2, ivec v2 ) -{ - dest[0] = k1*v1[0] + k2*v2[0]; - dest[1] = k1*v1[1] + k2*v2[1]; - dest[2] = k1*v1[2] + k2*v2[2]; -} - - -void ivec_Add( ivec dest, ivec v ) -{ - dest[0] += v[0]; - dest[1] += v[1]; - dest[2] += v[2]; -} - - -void ivec_ScaledAdd( ivec dest, int k, ivec v ) -{ - dest[0] += k * v[0]; - dest[1] += k * v[1]; - dest[2] += k * v[2]; -} - - - -void ivec_Max( ivec res, ivec v1, ivec v2 ) -{ - res[0] = MAX( v1[0], v2[0] ); - res[1] = MAX( v1[1], v2[1] ); - res[2] = MAX( v1[2], v2[2] ); -} - - -void ivec_Max3( ivec res, ivec v1, ivec v2, ivec v3 ) -{ - res[0] = MAX3( v1[0], v2[0], v3[0] ); - res[1] = MAX3( v1[1], v2[1], v3[1] ); - res[2] = MAX3( v1[2], v2[2], v3[2] ); -} diff --git a/src/USER-REAXC/reaxc_vector.h b/src/USER-REAXC/reaxc_vector.h index a33e249bcd..189b866dc0 100644 --- a/src/USER-REAXC/reaxc_vector.h +++ b/src/USER-REAXC/reaxc_vector.h @@ -31,16 +31,6 @@ #include "reaxc_types.h" #include "reaxc_defs.h" -int Vector_isZero( real*, int ); -void Vector_MakeZero( real*, int ); -void Vector_Copy( real*, real*, int ); -void Vector_Scale( real*, real, real*, int ); -void Vector_Sum( real*, real, real*, real, real*, int ); -void Vector_Add( real*, real, real*, int ); -real Dot( real*, real*, int ); -real Norm( real*, int ); -void Vector_Print( FILE*, char*, real*, int ); - void rvec_Copy( rvec, rvec ); void rvec_Scale( rvec, real, rvec ); void rvec_Add( rvec, rvec ); @@ -63,33 +53,12 @@ void rvec_MakeZero( rvec ); void rvec_Random( rvec ); void rtensor_MakeZero( rtensor ); -void rtensor_Multiply( rtensor, rtensor, rtensor ); void rtensor_MatVec( rvec, rtensor, rvec ); void rtensor_Scale( rtensor, real, rtensor ); -void rtensor_Add( rtensor, rtensor ); -void rtensor_ScaledAdd( rtensor, real, rtensor ); -void rtensor_Sum( rtensor, rtensor, rtensor ); -void rtensor_ScaledSum( rtensor, real, rtensor, real, rtensor ); -void rtensor_Scale( rtensor, real, rtensor ); -void rtensor_Copy( rtensor, rtensor ); -void rtensor_Identity( rtensor ); -void rtensor_Transpose( rtensor, rtensor ); -real rtensor_Det( rtensor ); -real rtensor_Trace( rtensor ); -void Print_rTensor(FILE*,rtensor); - -int ivec_isZero( ivec ); -int ivec_isEqual( ivec, ivec ); void ivec_MakeZero( ivec ); void ivec_Copy( ivec, ivec ); void ivec_Scale( ivec, real, ivec ); -void ivec_rScale( ivec, real, rvec ); void ivec_Sum( ivec, ivec, ivec ); -void ivec_ScaledSum( ivec, int, ivec, int, ivec ); -void ivec_Add( ivec, ivec ); -void ivec_ScaledAdd( ivec, int, ivec ); -void ivec_Max( ivec, ivec, ivec ); -void ivec_Max3( ivec, ivec, ivec, ivec ); #endif