jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8355 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2012-06-20 15:11:43 +00:00
parent 0cd34fb8a2
commit f3eb64fd0f
6 changed files with 314 additions and 199 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

122
src/USER-REAXC/fix_qeq_reax.cpp
View File

@ -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_j<A->firstnbr[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];

205
src/USER-REAXC/pair_reax_c.cpp
View File

@ -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<int> (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 ) {

134
src/neigh_half_bin.cpp
View File

@ -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

23
src/neigh_stencil.cpp
View File

@ -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)
{

27
src/neighbor.cpp
View File

@ -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) {

2
src/neighbor.h
View File

@ -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);