From f3eb64fd0f7adb24b5f1b39accc04a880cab1c86 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 20 Jun 2012 15:11:43 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8355 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-REAXC/fix_qeq_reax.cpp | 122 ++++++++++--------- src/USER-REAXC/pair_reax_c.cpp | 205 ++++++++++++-------------------- src/neigh_half_bin.cpp | 134 +++++++++++++++++++++ src/neigh_stencil.cpp | 23 ++++ src/neighbor.cpp | 27 +++-- src/neighbor.h | 2 + 6 files changed, 314 insertions(+), 199 deletions(-) diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index b61c149002..c79b79448f 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. @@ -53,11 +53,11 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : +FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - if (narg != 8) error->all(FLERR,"Illegal fix qeq/reax command"); - + if (narg != 8) error->all(FLERR,"Illegal fix qeq/reax command"); + nevery = atoi(arg[3]); swa = atof(arg[4]); swb = atof(arg[5]); @@ -73,7 +73,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : s = NULL; t = NULL; nprev = 5; - + Hdia_inv = NULL; b_s = NULL; b_t = NULL; @@ -111,7 +111,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : FixQEqReax::~FixQEqReax() { // unregister callbacks to this fix from Atom class - + atom->delete_callback(id,0); memory->destroy(s_hist); @@ -153,7 +153,7 @@ void FixQEqReax::pertype_parameters(char *arg) gamma = (double *) pair->extract("gamma",tmp); if (chi == NULL || eta == NULL || gamma == NULL) error->all(FLERR, - "Fix qeq/reax could not extract params from pair reax/c"); + "Fix qeq/reax could not extract params from pair reax/c"); return; } @@ -171,11 +171,11 @@ void FixQEqReax::pertype_parameters(char *arg) if (comm->me == 0) { if ((pf = fopen(arg,"r")) == NULL) error->one(FLERR,"Fix qeq/reax parameter file could not be found"); - + for (i = 1; i <= ntypes && !feof(pf); i++) { fscanf(pf,"%d %lg %lg %lg",&itype,&v1,&v2,&v3); if (itype < 1 || itype > ntypes) - error->one(FLERR,"Fix qeq/reax invalid atom type in param file"); + error->one(FLERR,"Fix qeq/reax invalid atom type in param file"); chi[itype] = v1; eta[itype] = v2; gamma[itype] = v3; @@ -255,7 +255,7 @@ void FixQEqReax::allocate_matrix() m += list->numneigh[i]; } m_cap = MAX( (int)(m * SAFE_ZONE), MIN_CAP * MIN_NBRS ); - + H.n = n_cap; H.m = m_cap; memory->create(H.firstnbr,n_cap,"qeq:H.firstnbr"); @@ -271,7 +271,7 @@ void FixQEqReax::deallocate_matrix() memory->destroy( H.firstnbr ); memory->destroy( H.numnbrs ); memory->destroy( H.jlist ); - memory->destroy( H.val ); + memory->destroy( H.val ); } /* ---------------------------------------------------------------------- */ @@ -287,14 +287,21 @@ void FixQEqReax::reallocate_matrix() void FixQEqReax::init() { if (!atom->q_flag) error->all(FLERR,"Fix qeq/reax requires atom attribute q"); - + // need a half neighbor list w/ Newton off // built whenever re-neighboring occurs int irequest = neighbor->request(this); neighbor->requests[irequest]->pair = 0; neighbor->requests[irequest]->fix = 1; + + // Ray + //neighbor->requests[irequest]->half = 0; + //neighbor->requests[irequest]->full = 1; + //neighbor->requests[irequest]->ghost = 2; // 2 for newton off + neighbor->requests[irequest]->newton = 2; + neighbor->requests[irequest]->ghost = 1; init_shielding(); init_taper(); @@ -319,7 +326,7 @@ void FixQEqReax::init_shielding() ntypes = atom->ntypes; memory->create(shld,ntypes+1,ntypes+1,"qeq:shileding"); - + for( i = 1; i <= ntypes; ++i ) for( j = 1; j <= ntypes; ++j ) shld[i][j] = pow( gamma[i] * gamma[j], -1.5 ); @@ -338,11 +345,12 @@ void FixQEqReax::init_taper() else if (swb < 5 && comm->me == 0) error->warning(FLERR,"Fix qeq/reax has very low Taper radius cutoff"); - d7 = pow( swb - swa, 7.0 ); + d7 = pow( swb - swa, 7 ); swa2 = SQR( swa ); swa3 = CUBE( swa ); swb2 = SQR( swb ); swb3 = CUBE( swb ); + Tap[7] = 20.0 / d7; Tap[6] = -70.0 * (swa + swb) / d7; Tap[5] = 84.0 * (swa2 + 3.0*swa*swb + swb2) / d7; @@ -351,7 +359,7 @@ void FixQEqReax::init_taper() Tap[2] =-210.0 * (swa3*swb2 + swa2*swb3) / d7; Tap[1] = 140.0 * swa3 * swb3 / d7; Tap[0] = (-35.0*swa3*swb2*swb2 + 21.0*swa2*swb3*swb2 + - 7.0*swa*swb3*swb3 + swb3*swb3*swb ) / d7; + 7.0*swa*swb3*swb3 + swb3*swb3*swb ) / d7; } /* ---------------------------------------------------------------------- */ @@ -364,7 +372,7 @@ void FixQEqReax::setup_pre_force(int vflag) allocate_matrix(); pre_force(vflag); -} +} /* ---------------------------------------------------------------------- */ @@ -372,14 +380,14 @@ void FixQEqReax::setup_pre_force_respa(int vflag, int ilevel) { if (ilevel < nlevels_respa-1) return; setup_pre_force(vflag); -} +} /* ---------------------------------------------------------------------- */ void FixQEqReax::min_setup_pre_force(int vflag) { setup_pre_force(vflag); -} +} /* ---------------------------------------------------------------------- */ @@ -387,7 +395,7 @@ void FixQEqReax::init_storage() { N = atom->nlocal + atom->nghost; for( int i = 0; i < N; i++ ) { - Hdia_inv[i] = 1. / eta[atom->type[i]]; + Hdia_inv[i] = 1. / eta[atom->type[i]]; b_s[i] = -chi[atom->type[i]]; b_t[i] = -1.0; b_prc[i] = 0; @@ -412,9 +420,9 @@ void FixQEqReax::pre_force(int vflag) if( atom->nmax > nmax ) reallocate_storage(); if( n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE ) reallocate_matrix(); - + init_matvec(); - matvecs = CG(b_s, s); // CG on s - parallel + matvecs = CG(b_s, s); // CG on s - parallel matvecs += CG(b_t, t); // CG on t - parallel calculate_Q(); @@ -456,7 +464,7 @@ void FixQEqReax::init_matvec() //t[i] = 2 * t_hist[i][0] - t_hist[i][1]; /* quadratic extrapolation for s & t from previous solutions */ - //s[i] = s_hist[i][2] + 3 * ( s_hist[i][0] - s_hist[i][1] ); + //s[i] = s_hist[i][2] + 3 * ( s_hist[i][0] - s_hist[i][1] ); t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1] ); /* cubic extrapolation for s & t from previous solutions */ @@ -488,7 +496,7 @@ void FixQEqReax::compute_H() ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - + // fill in the H matrix m_fill = 0; r_sqr = 0; @@ -497,20 +505,20 @@ void FixQEqReax::compute_H() jlist = firstneigh[i]; jnum = numneigh[i]; H.firstnbr[i] = m_fill; - + for( jj = 0; jj < jnum; jj++ ) { j = jlist[jj]; - + dx = x[j][0] - x[i][0]; dy = x[j][1] - x[i][1]; dz = x[j][2] - x[i][2]; r_sqr = SQR(dx) + SQR(dy) + SQR(dz); - + flag = 0; if (r_sqr <= SQR(swb)) { if (j < n) flag = 1; else if (tag[i] < tag[j]) flag = 1; - else if (tag[i] == tag[j]) { + else if (tag[i] == tag[j]) { if (dz > SMALL) flag = 1; else if (fabs(dz) < SMALL) { if (dy > SMALL) flag = 1; @@ -519,21 +527,21 @@ void FixQEqReax::compute_H() } } } - + if( flag ) { - H.jlist[m_fill] = j; - H.val[m_fill] = calculate_H( sqrt(r_sqr), shld[type[i]][type[j]] ); - m_fill++; + H.jlist[m_fill] = j; + H.val[m_fill] = calculate_H( sqrt(r_sqr), shld[type[i]][type[j]] ); + m_fill++; } } - + H.numnbrs[i] = m_fill - H.firstnbr[i]; } if (m_fill >= H.m) { char str[128]; sprintf(str,"H matrix size has been exceeded: m_fill=%d H.m=%d\n", - m_fill, H.m ); + m_fill, H.m ); error->warning(FLERR,str); error->all(FLERR,"Fix qeq/reax has insufficient QEq matrix size"); } @@ -583,21 +591,21 @@ int FixQEqReax::CG( double *b, double *x ) comm->forward_comm_fix(this); //Dist_vector( d ); sparse_matvec( &H, d, q ); comm->reverse_comm_fix(this); //Coll_vector( q ); - + tmp = parallel_dot( d, q, n ); alpha = sig_new / tmp; // comm->me, i, parallel_norm( d, n ), parallel_norm( q, n ), tmp ); - + vector_add( x, alpha, d, n ); vector_add( r, -alpha, q, n ); - + // pre-conditioning for( j = 0; j < n; ++j ) p[j] = r[j] * Hdia_inv[j]; - + sig_old = sig_new; sig_new = parallel_dot( r, p, n ); - + beta = sig_new / sig_old; vector_sum( d, 1., p, beta, d, n ); @@ -619,7 +627,7 @@ void FixQEqReax::sparse_matvec( sparse_matrix *A, double *x, double *b ) b[i] = eta[ atom->type[i] ] * x[i]; for( i = n; i < N; ++i ) b[i] = 0; - + for( i = 0; i < n; ++i ) { for( itr_j=A->firstnbr[i]; itr_jfirstnbr[i]+A->numnbrs[i]; itr_j++) { j = A->jlist[itr_j]; @@ -640,10 +648,10 @@ void FixQEqReax::calculate_Q() s_sum = parallel_vector_acc( s, n ); t_sum = parallel_vector_acc( t, n); u = s_sum / t_sum; - + for( i = 0; i < n; ++i ) { q[i] = s[i] - u * t[i]; - + /* backup s & t */ for( k = 4; k > 0; --k ) { s_hist[i][k] = s_hist[i][k-1]; @@ -659,12 +667,12 @@ void FixQEqReax::calculate_Q() /* ---------------------------------------------------------------------- */ -int FixQEqReax::pack_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) +int FixQEqReax::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) { int m; - if( pack_flag == 1) + if( pack_flag == 1) for(m = 0; m < n; m++) buf[m] = d[list[m]]; else if( pack_flag == 2 ) for(m = 0; m < n; m++) buf[m] = s[list[m]]; @@ -675,23 +683,23 @@ int FixQEqReax::pack_comm(int n, int *list, double *buf, return 1; } - + /* ---------------------------------------------------------------------- */ void FixQEqReax::unpack_comm(int n, int first, double *buf) { int i, m; - - if( pack_flag == 1) + + if( pack_flag == 1) for(m = 0, i = first; m < n; m++, i++) d[i] = buf[m]; - else if( pack_flag == 2) + else if( pack_flag == 2) for(m = 0, i = first; m < n; m++, i++) s[i] = buf[m]; - else if( pack_flag == 3) + else if( pack_flag == 3) for(m = 0, i = first; m < n; m++, i++) t[i] = buf[m]; - else if( pack_flag == 4) + else if( pack_flag == 4) for(m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m]; } - + /* ---------------------------------------------------------------------- */ int FixQEqReax::pack_reverse_comm(int n, int first, double *buf) @@ -700,7 +708,7 @@ int FixQEqReax::pack_reverse_comm(int n, int first, double *buf) for(m = 0, i = first; m < n; m++, i++) buf[m] = q[i]; return 1; } - + /* ---------------------------------------------------------------------- */ void FixQEqReax::unpack_reverse_comm(int n, int *list, double *buf) @@ -720,7 +728,7 @@ double FixQEqReax::memory_usage() bytes += atom->nmax*11 * sizeof(double); // storage bytes += n_cap*2 * sizeof(int); // matrix... bytes += m_cap * sizeof(int); - bytes += m_cap * sizeof(double); + bytes += m_cap * sizeof(double); return bytes; } @@ -823,7 +831,7 @@ double FixQEqReax::parallel_vector_acc( double *v, int n ) double FixQEqReax::norm( double* v1, int k ) { double ret = 0; - + for( --k; k>=0; --k ) ret += ( v1[k] * v1[k] ); @@ -832,8 +840,8 @@ double FixQEqReax::norm( double* v1, int k ) /* ---------------------------------------------------------------------- */ -void FixQEqReax::vector_sum( double* dest, double c, double* v, - double d, double* y, int k ) +void FixQEqReax::vector_sum( double* dest, double c, double* v, + double d, double* y, int k ) { for( --k; k>=0; --k ) dest[k] = c * v[k] + d * y[k]; @@ -852,7 +860,7 @@ void FixQEqReax::vector_scale( double* dest, double c, double* v, int k ) double FixQEqReax::dot( double* v1, double* v2, int k ) { double ret = 0; - + for( --k; k>=0; --k ) ret += v1[k] * v2[k]; diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp index 2195132a0f..eba3a6d89e 100644 --- a/src/USER-REAXC/pair_reax_c.cpp +++ b/src/USER-REAXC/pair_reax_c.cpp @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. @@ -59,16 +59,17 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp) single_enable = 0; restartinfo = 0; one_coeff = 1; + ghostneigh = 1; system = (reax_system *) memory->smalloc(sizeof(reax_system),"reax:system"); - control = (control_params *) + control = (control_params *) memory->smalloc(sizeof(control_params),"reax:control"); data = (simulation_data *) memory->smalloc(sizeof(simulation_data),"reax:data"); workspace = (storage *) memory->smalloc(sizeof(storage),"reax:storage"); - lists = (reax_list *) + lists = (reax_list *) memory->smalloc(LIST_N * sizeof(reax_list),"reax:lists"); out_control = (output_controls *) memory->smalloc(sizeof(output_controls),"reax:out_control"); @@ -128,7 +129,7 @@ PairReaxC::~PairReaxC() DeAllocate_System( system ); } //fprintf( stderr, "4\n" ); - + memory->destroy( system ); memory->destroy( control ); memory->destroy( data ); @@ -142,6 +143,7 @@ PairReaxC::~PairReaxC() if( allocated ) { memory->destroy(setflag); memory->destroy(cutsq); + memory->destroy(cutghost); delete [] map; delete [] chi; @@ -163,6 +165,7 @@ void PairReaxC::allocate( ) memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(cutghost,n+1,n+1,"pair:cutghost"); map = new int[n+1]; chi = new double[n+1]; @@ -190,7 +193,7 @@ void PairReaxC::settings(int narg, char **arg) control->hbond_cut = 7.50; control->thb_cut = 0.001; control->thb_cutsq = 0.00001; - + out_control->write_steps = 0; out_control->traj_method = 0; strcpy( out_control->traj_title, "default_title" ); @@ -246,7 +249,7 @@ void PairReaxC::coeff( int nargs, char **args ) // read ffield file Read_Force_Field(args[2], &(system->reax_param), control); - + // read args that map atom types to elements in potential file // map[i] = which element the Ith atom type is, -1 if NULL @@ -276,7 +279,7 @@ void PairReaxC::coeff( int nargs, char **args ) setflag[i][j] = 1; count++; } - + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -290,7 +293,7 @@ void PairReaxC::init_style( ) int iqeq; for (iqeq = 0; iqeq < modify->nfix; iqeq++) if (strcmp(modify->fix[iqeq]->style,"qeq/reax") == 0) break; - if (iqeq == modify->nfix && qeqflag == 1) + if (iqeq == modify->nfix && qeqflag == 1) error->all(FLERR,"Pair reax/c requires use of fix qeq/reax"); system->n = atom->nlocal; // my atoms @@ -299,8 +302,8 @@ void PairReaxC::init_style( ) system->wsize = comm->nprocs; system->big_box.V = 0; - system->big_box.box_norms[0] = 0; - system->big_box.box_norms[1] = 0; + system->big_box.box_norms[0] = 0; + system->big_box.box_norms[1] = 0; system->big_box.box_norms[2] = 0; if (atom->tag_enable == 0) @@ -308,11 +311,28 @@ void PairReaxC::init_style( ) if (force->newton_pair == 0) error->all(FLERR,"Pair style reax/c requires newton pair on"); + // original // need a half neighbor list w/ Newton off // built whenever re-neighboring occurs + //int irequest = neighbor->request(this); + //neighbor->requests[irequest]->newton = 2; + + // Ray + // need a full+ghost neighbor list w/ Newton off + // built whenever re-neighboring occurs + + //int irequest = neighbor->request(this); + //neighbor->requests[irequest]->half = 0; + //neighbor->requests[irequest]->full = 1; + //neighbor->requests[irequest]->ghost = 2; // 2 for newton off + + // need a half neighbor list with ghosts w/ Newton off + // built whenever re-neighboring occurs + int irequest = neighbor->request(this); neighbor->requests[irequest]->newton = 2; + neighbor->requests[irequest]->ghost = 1; cutmax = MAX3(control->nonb_cut, control->hbond_cut, 2*control->bond_cut); @@ -344,7 +364,7 @@ void PairReaxC::setup( ) if (setup_flag == 0) { setup_flag = 1; - + int *num_bonds = fix_reax->num_bonds; int *num_hbonds = fix_reax->num_hbonds; @@ -359,15 +379,15 @@ void PairReaxC::setup( ) PreAllocate_Space( system, control, workspace, world ); write_reax_atoms(); - + int num_nbrs = estimate_reax_lists(); - if(!Make_List(system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR, - lists+FAR_NBRS, world)) + if(!Make_List(system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR, + lists+FAR_NBRS, world)) error->all(FLERR,"Pair reax/c problem in far neighbor list"); - + write_reax_lists(); - Initialize( system, control, data, workspace, &lists, out_control, - mpi_data, world ); + Initialize( system, control, data, workspace, &lists, out_control, + mpi_data, world ); for( int k = 0; k < system->N; ++k ) { num_bonds[k] = system->my_atoms[k].num_bonds; num_hbonds[k] = system->my_atoms[k].num_hbonds; @@ -383,7 +403,7 @@ void PairReaxC::setup( ) for(int k = oldN; k < system->N; ++k) Set_End_Index( k, Start_Index( k, lists+BONDS ), lists+BONDS ); - + // check if I need to shrink/extend my data-structs ReAllocate( system, control, data, workspace, &lists, mpi_data ); @@ -394,6 +414,7 @@ void PairReaxC::setup( ) double PairReaxC::init_one(int i, int j) { + cutghost[i][j] = cutghost[j][i] = cutmax; return cutmax; } @@ -417,9 +438,9 @@ void PairReaxC::compute(int eflag, int vflag) /* if ((eflag_atom || vflag_atom) && firstwarn) { firstwarn = 0; - if (comm->me == 0) + if (comm->me == 0) error->warning(FLERR,"Pair reax/c cannot yet compute " - "per-atom energy or stress"); + "per-atom energy or stress"); } */ if (vflag_global) control->virial = 1; @@ -430,15 +451,15 @@ void PairReaxC::compute(int eflag, int vflag) system->bigN = static_cast (atom->natoms); // all atoms in the system system->big_box.V = 0; - system->big_box.box_norms[0] = 0; - system->big_box.box_norms[1] = 0; + system->big_box.box_norms[0] = 0; + system->big_box.box_norms[1] = 0; system->big_box.box_norms[2] = 0; if( comm->me == 0 ) t_start = MPI_Wtime(); // setup data structures setup(); - + Reset( system, control, data, workspace, &lists, world ); workspace->realloc.num_far = write_reax_lists(); // timing for filling in the reax lists @@ -481,7 +502,7 @@ void PairReaxC::compute(int eflag, int vflag) // Store the different parts of the energy // in a list for output by compute pair command - pvector[0] = data->my_en.e_bond; + pvector[0] = data->my_en.e_bond; pvector[1] = data->my_en.e_ov + data->my_en.e_un; pvector[2] = data->my_en.e_lp; pvector[3] = 0.0; @@ -512,11 +533,11 @@ void PairReaxC::compute(int eflag, int vflag) Output_Results( system, control, data, &lists, out_control, mpi_data ); - if(fixbond_flag) - fixbond( system, control, data, &lists, out_control, mpi_data ); + if(fixbond_flag) + fixbond( system, control, data, &lists, out_control, mpi_data ); - if(fixspecies_flag) - fixspecies( system, control, data, &lists, out_control, mpi_data ); + if(fixspecies_flag) + fixspecies( system, control, data, &lists, out_control, mpi_data ); } @@ -526,7 +547,7 @@ void PairReaxC::write_reax_atoms() { int *num_bonds = fix_reax->num_bonds; int *num_hbonds = fix_reax->num_hbonds; - + for( int i = 0; i < system->N; ++i ){ system->my_atoms[i].orig_id = atom->tag[i]; system->my_atoms[i].type = map[atom->type[i]]; @@ -551,8 +572,8 @@ void PairReaxC::get_distance( rvec xj, rvec xi, double *d_sqr, rvec *dvec ) /* ---------------------------------------------------------------------- */ -void PairReaxC::set_far_nbr( far_neighbor_data *fdest, - int j, double d, rvec dvec ) +void PairReaxC::set_far_nbr( far_neighbor_data *fdest, + int j, double d, rvec dvec ) { fdest->nbr = j; fdest->d = d; @@ -588,7 +609,11 @@ int PairReaxC::estimate_reax_lists() marked = (int*) calloc( system->N, sizeof(int) ); dist = (double*) calloc( system->N, sizeof(double) ); - for( itr_i = 0; itr_i < list->inum; ++itr_i ){ + int inum = list->inum; + int gnum = list->gnum; + int numall = inum + gnum; + + for( itr_i = 0; itr_i < inum+gnum; ++itr_i ){ i = ilist[itr_i]; marked[i] = 1; ++num_marked; @@ -598,54 +623,12 @@ int PairReaxC::estimate_reax_lists() j = jlist[itr_j]; j &= NEIGHMASK; get_distance( x[j], x[i], &d_sqr, &dvec ); - dist[j] = sqrt(d_sqr); - - if( dist[j] <= control->nonb_cut ) - ++num_nbrs; - } - - // compute the nbrs among ghost atoms - for( itr_j = 0; itr_j < numneigh[i]; ++itr_j ){ - j = jlist[itr_j]; - j &= NEIGHMASK; - - if( j >= nlocal && !marked[j] && - dist[j] <= (control->vlist_cut - control->bond_cut) ){ - marked[j] = 1; - ++num_marked; - - for( itr_g = 0; itr_g < numneigh[i]; ++itr_g ){ - g = jlist[itr_g]; - g &= NEIGHMASK; - - if( g >= nlocal && !marked[g] ){ - get_distance( x[g], x[j], &g_d_sqr, &g_dvec ); - //g_dvec[0] = x[g][0] - x[j][0]; - //g_dvec[1] = x[g][1] - x[j][1]; - //g_dvec[2] = x[g][2] - x[j][2]; - //g_d_sqr = SQR(g_dvec[0]) + SQR(g_dvec[1]) + SQR(g_dvec[2]); - - if( g_d_sqr <= SQR(control->bond_cut) ) - ++num_nbrs; - } - } - } + + if( d_sqr <= SQR(control->nonb_cut) ) + ++num_nbrs; } } - for( i = 0; i < system->N; ++i ) - if( !marked[i] ) { - marked[i] = 1; - ++num_marked; - - for( j = i+1; j < system->N; ++j ) - if( !marked[j] ) { - get_distance( x[j], x[i], &d_sqr, &dvec ); - if( d_sqr <= SQR(control->bond_cut) ) - ++num_nbrs; - } - } - free( marked ); free( dist ); @@ -676,11 +659,15 @@ int PairReaxC::write_reax_lists() far_nbrs = lists + FAR_NBRS; far_list = far_nbrs->select.far_nbr_list; - num_nbrs = 0; + num_nbrs = 0; marked = (int*) calloc( system->N, sizeof(int) ); dist = (double*) calloc( system->N, sizeof(double) ); - for( itr_i = 0; itr_i < list->inum; ++itr_i ){ + int inum = list->inum; + int gnum = list->gnum; + int numall = inum + gnum; + + for( itr_i = 0; itr_i < inum+gnum; ++itr_i ){ i = ilist[itr_i]; marked[i] = 1; jlist = firstneigh[i]; @@ -690,70 +677,24 @@ int PairReaxC::write_reax_lists() j = jlist[itr_j]; j &= NEIGHMASK; get_distance( x[j], x[i], &d_sqr, &dvec ); - dist[j] = sqrt( d_sqr ); - if( dist[j] <= control->nonb_cut ){ - set_far_nbr( &far_list[num_nbrs], j, dist[j], dvec ); - ++num_nbrs; + if( d_sqr <= (control->nonb_cut*control->nonb_cut) ){ + dist[j] = sqrt( d_sqr ); + set_far_nbr( &far_list[num_nbrs], j, dist[j], dvec ); + ++num_nbrs; } } Set_End_Index( i, num_nbrs, far_nbrs ); - - // compute the nbrs among ghost atoms - for( itr_j = 0; itr_j < numneigh[i]; ++itr_j ){ - j = jlist[itr_j]; - j &= NEIGHMASK; - - if( j >= nlocal && !marked[j] && - dist[j] <= (control->vlist_cut - control->bond_cut) ){ - marked[j] = 1; - Set_Start_Index( j, num_nbrs, far_nbrs ); - - for( itr_g = 0; itr_g < numneigh[i]; ++itr_g ){ - g = jlist[itr_g]; - g &= NEIGHMASK; - - if( g >= nlocal && !marked[g] ){ - get_distance( x[g], x[j], &g_d_sqr, &g_dvec ); - - if( g_d_sqr <= SQR(control->bond_cut) ){ - g_d = sqrt( g_d_sqr ); - - set_far_nbr( &far_list[num_nbrs], g, g_d, g_dvec ); - ++num_nbrs; - } - } - } - Set_End_Index( j, num_nbrs, far_nbrs ); - } - } } - for( i = 0; i < system->N; ++i ) - if( !marked[i] ) { - marked[i] = 1; - Set_Start_Index( i, num_nbrs, far_nbrs ); - - for( j = i+1; j < system->N; ++j ) - if( !marked[j] ) { - get_distance( x[j], x[i], &d_sqr, &dvec ); - if( d_sqr <= SQR(control->bond_cut) ) { - set_far_nbr( &far_list[num_nbrs], j, sqrt(d_sqr), dvec ); - ++num_nbrs; - } - } - - Set_End_Index( i, num_nbrs, far_nbrs ); - } - free( marked ); free( dist ); - + return num_nbrs; } - + /* ---------------------------------------------------------------------- */ - + void PairReaxC::read_reax_forces() { for( int i = 0; i < system->N; ++i ) { diff --git a/src/neigh_half_bin.cpp b/src/neigh_half_bin.cpp index aa748c5f7e..dfe8456aea 100644 --- a/src/neigh_half_bin.cpp +++ b/src/neigh_half_bin.cpp @@ -119,6 +119,140 @@ void Neighbor::half_bin_no_newton(NeighList *list) list->inum = inum; } +/* ---------------------------------------------------------------------- + binned neighbor list construction with partial Newton's 3rd law + include neighbors of ghost atoms, but no "special neighbors" for ghosts + owned and ghost atoms check own bin and other bins in stencil + pair stored once if i,j are both owned and i < j + pair stored by me if i owned and j ghost (also stored by proc owning j) + pair stored once if i,j are both ghost and i < j +------------------------------------------------------------------------- */ + +void Neighbor::half_bin_no_newton_ghost(NeighList *list) +{ + int i,j,k,n,itype,jtype,ibin,which; + double xtmp,ytmp,ztmp,delx,dely,delz,rsq; + int xbin,ybin,zbin,xbin2,ybin2,zbin2; + int *neighptr; + + // bin local & ghost atoms + + bin_atoms(); + + // loop over each atom, storing neighbors + + int **special = atom->special; + int **nspecial = atom->nspecial; + int *tag = atom->tag; + + double **x = atom->x; + int *type = atom->type; + int *mask = atom->mask; + int *molecule = atom->molecule; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int molecular = atom->molecular; + + int *ilist = list->ilist; + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + int **pages = list->pages; + int nstencil = list->nstencil; + int *stencil = list->stencil; + int **stencilxyz = list->stencilxyz; + + int inum = 0; + int npage = 0; + int npnt = 0; + + for (i = 0; i < nall; i++) { + + if (pgsize - npnt < oneatom) { + npnt = 0; + npage++; + if (npage == list->maxpage) pages = list->add_pages(); + } + + neighptr = &pages[npage][npnt]; + n = 0; + + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // loop over all atoms in other bins in stencil including self + // when i is a ghost atom, must check if stencil bin is out of bounds + // only store pair if i < j + // stores own/own pairs only once + // stores own/ghost pairs with owned atom only, on both procs + // stores ghost/ghost pairs only once + + if (i < nlocal) { + ibin = coord2bin(x[i]); + + for (k = 0; k < nstencil; k++) { + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighsq[itype][jtype]) { + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; + } + } + } + + } else { + ibin = coord2bin(x[i],xbin,ybin,zbin); + for (k = 0; k < nstencil; k++) { + xbin2 = xbin + stencilxyz[k][0]; + ybin2 = ybin + stencilxyz[k][1]; + zbin2 = zbin + stencilxyz[k][2]; + if (xbin2 < 0 || xbin2 >= mbinx || + ybin2 < 0 || ybin2 >= mbiny || + zbin2 < 0 || zbin2 >= mbinz) continue; + for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + if (j <= i) continue; + + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if (rsq <= cutneighghostsq[itype][jtype]) + neighptr[n++] = j; + } + } + } + + ilist[inum++] = i; + firstneigh[i] = neighptr; + numneigh[i] = n; + npnt += n; + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); + } + + list->inum = atom->nlocal; + list->gnum = inum - atom->nlocal; +} + /* ---------------------------------------------------------------------- binned neighbor list construction with full Newton's 3rd law each owned atom i checks its own bin and other bins in Newton stencil diff --git a/src/neigh_stencil.cpp b/src/neigh_stencil.cpp index 34a5ddc305..53c71286da 100644 --- a/src/neigh_stencil.cpp +++ b/src/neigh_stencil.cpp @@ -79,6 +79,29 @@ void Neighbor::stencil_half_bin_3d_no_newton(NeighList *list, /* ---------------------------------------------------------------------- */ +void Neighbor::stencil_half_ghost_bin_3d_no_newton(NeighList *list, + int sx, int sy, int sz) +{ + int i,j,k; + int *stencil = list->stencil; + int **stencilxyz = list->stencilxyz; + int nstencil = 0; + + for (k = -sz; k <= sz; k++) + for (j = -sy; j <= sy; j++) + for (i = -sx; i <= sx; i++) + if (bin_distance(i,j,k) < cutneighmaxsq) { + stencilxyz[nstencil][0] = i; + stencilxyz[nstencil][1] = j; + stencilxyz[nstencil][2] = k; + stencil[nstencil++] = k*mbiny*mbinx + j*mbinx + i; + } + + list->nstencil = nstencil; +} + +/* ---------------------------------------------------------------------- */ + void Neighbor::stencil_half_bin_2d_newton(NeighList *list, int sx, int sy, int sz) { diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 8c8183bf2d..d3b06d1409 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -493,8 +493,10 @@ void Neighbor::init() lists[i]->listskip = lists[requests[i]->otherlist]; lists[i]->copy_skip_info(requests[i]->iskip,requests[i]->ijskip); - } else if (requests[i]->half_from_full) + } else if (requests[i]->half_from_full) { lists[i]->listfull = lists[i-1]; + printf("AAA\n"); + } else if (requests[i]->granhistory) { lists[i-1]->listgranhistory = lists[i]; @@ -520,6 +522,7 @@ void Neighbor::init() requests[i]->half = 0; requests[i]->half_from_full = 1; lists[i]->listfull = lists[j]; + printf("BBB %d %d\n",nlist,j); } } else if (requests[i]->fix || requests[i]->compute) { @@ -547,6 +550,7 @@ void Neighbor::init() requests[i]->half = 0; requests[i]->half_from_full = 1; lists[i]->listfull = lists[j]; + printf("CCC\n"); } } } @@ -829,9 +833,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq) else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton; } else if (rq->newton == 1) { pb = &Neighbor::half_nsq_newton; - } else if (rq->newton == 2) { - pb = &Neighbor::half_nsq_no_newton; - } + } else if (rq->newton == 2) pb = &Neighbor::half_nsq_no_newton; } else if (style == BIN) { if (rq->newton == 0) { if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton; @@ -840,7 +842,10 @@ void Neighbor::choose_build(int index, NeighRequest *rq) } else if (rq->newton == 1) { if (triclinic == 0) pb = &Neighbor::half_bin_newton; else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri; - } else if (rq->newton == 2) pb = &Neighbor::half_bin_no_newton; + } else if (rq->newton == 2) { + if (rq->ghost) pb = &Neighbor::half_bin_no_newton_ghost; + else pb = &Neighbor::half_bin_no_newton; + } } else if (style == MULTI) { if (rq->newton == 0) { if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton; @@ -982,9 +987,9 @@ void Neighbor::choose_build(int index, NeighRequest *rq) // general error check - if (rq->ghost && !rq->full) - error->all(FLERR, - "Neighbors of ghost atoms only allowed for full neighbor lists"); + //if (rq->ghost && !rq->full) + // error->all(FLERR, + // "Neighbors of ghost atoms only allowed for full neighbor lists"); pair_build[index] = pb; } @@ -1038,8 +1043,10 @@ void Neighbor::choose_stencil(int index, NeighRequest *rq) } else if (rq->newton == 2) { if (dimension == 2) sc = &Neighbor::stencil_half_bin_2d_no_newton; - else if (dimension == 3) - sc = &Neighbor::stencil_half_bin_3d_no_newton; + else if (dimension == 3) { + if (rq->ghost) sc = &Neighbor::stencil_half_ghost_bin_3d_no_newton; + else sc = &Neighbor::stencil_half_bin_3d_no_newton; + } } } else if (style == MULTI) { diff --git a/src/neighbor.h b/src/neighbor.h index 3cffa42bf8..b93e2fab9a 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -180,6 +180,7 @@ class Neighbor : protected Pointers { void half_nsq_newton(class NeighList *); void half_bin_no_newton(class NeighList *); + void half_bin_no_newton_ghost(class NeighList *); void half_bin_newton(class NeighList *); void half_bin_newton_tri(class NeighList *); @@ -223,6 +224,7 @@ class Neighbor : protected Pointers { void stencil_half_bin_2d_no_newton(class NeighList *, int, int, int); void stencil_half_bin_3d_no_newton(class NeighList *, int, int, int); + void stencil_half_ghost_bin_3d_no_newton(class NeighList *, int, int, int); void stencil_half_bin_2d_newton(class NeighList *, int, int, int); void stencil_half_bin_3d_newton(class NeighList *, int, int, int); void stencil_half_bin_2d_newton_tri(class NeighList *, int, int, int);