diff --git a/PuReMD/src/#basic_comm.c# b/PuReMD/src/#basic_comm.c#
new file mode 100644
index 0000000000000000000000000000000000000000..efa01e7b0aacc9ee7c63e51bdee9dad18556d31f
--- /dev/null
+++ b/PuReMD/src/#basic_comm.c#
@@ -0,0 +1,322 @@
+/*----------------------------------------------------------------------
+ PuReMD - Purdue ReaxFF Molecular Dynamics Program
+
+ Copyright (2010) Purdue University
+ Hasan Metin Aktulga, haktulga@cs.purdue.edu
+ Joseph Fogarty, jcfogart@mail.usf.edu
+ Sagar Pandit, pandit@usf.edu
+ Ananth Y Grama, ayg@cs.purdue.edu
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU General Public License as
+ published by the Free Software Foundation; either version 2 of
+ the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ See the GNU General Public License for more details:
+ <http://www.gnu.org/licenses/>.
+ ----------------------------------------------------------------------*/
+
+#include "reax_types.h"
+#if defined(PURE_REAX)
+#include "basic_comm.h"
+#include "vector.h"
+#elif defined(LAMMPS_REAX)
+#include "reax_basic_comm.h"
+#include "reax_vector.h"
+#endif
+
+#if defined(PURE_REAX)
+void real_packer( void *dummy, mpi_out_data *out_buf )
+{
+ int i;
+ real *buf = (real*) dummy;
+ real *out = (real*) out_buf->out_atoms;
+
+ for ( i = 0; i < out_buf->cnt; ++i )
+ out[i] = buf[ out_buf->index[i] ];
+}
+
+
+void rvec_packer( void *dummy, mpi_out_data *out_buf )
+{
+ int i;
+ rvec *buf = (rvec*) dummy;
+ rvec *out = (rvec*)out_buf->out_atoms;
+
+ for ( i = 0; i < out_buf->cnt; ++i )
+ memcpy( out[i], buf[ out_buf->index[i] ], sizeof(rvec) );
+}qsdfsdfsddsfsd
+
+
+VOIDuuuuuuuuu RVEC2_packer( void *dummy, mpi_out_data *out_buf )
+{
+ int i;
+ rvec2 *buf = (rvec2*) dummy;
+ rvec2 *out = (rvec2*) out_buf->out_atoms;
+
+ for ( i = 0; i < out_buf->cnt; ++i )
+ memcpy( out[i], buf[ out_buf->index[i] ], sizeof(rvec2) );
+}
+
+
+void Dist( reax_system* system, mpi_datatypes *mpi_data,
+ void *buf, MPI_Datatype type, int scale, dist_packer pack )
+{
+ int d;
+ mpi_out_data *out_bufs;
+ MPI_Comm comm;
+ MPI_Request req1, req2;
+q MPI_Status stat1, stat2;
+ neighbor_proc *nbr1, *nbr2;
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d dist: entered\n", system->my_rank );
+#endif
+ comm = mpi_data->comm_mesh3D;
+ out_bufs = mpi_data->out_buffers;
+
+ for ( d = 0; d < 3; ++d )
+ {
+ /* initiate recvs */
+ nbr1 = &(system->my_nbrs[2 * d]);
+ if ( nbr1->atoms_cnt )
+ MPI_Irecv( buf + nbr1->atoms_str * scale, nbr1->atoms_cnt, type,
+ nbr1->rank, 2 * d + 1, comm, &req1 );
+
+ nbr2 = &(system->my_nbrs[2 * d + 1]);
+ if ( nbr2->atoms_cnt )
+ MPI_Irecv( buf + nbr2->atoms_str * scale, nbr2->atoms_cnt, type,
+ nbr2->rank, 2 * d, comm, &req2 );
+
+ /* send both messages in dimension d */
+ if ( out_bufs[2 * d].cnt )
+ {
+ pack( buf, out_bufs + (2 * d) );
+ MPI_Send( out_bufs[2 * d].out_atoms, out_bufs[2 * d].cnt, type,
+ nbr1->rank, 2 * d, comm );
+ }
+
+ if ( out_bufs[2 * d + 1].cnt )
+ {
+ pack( buf, out_bufs + (2 * d + 1) );
+ MPI_Send( out_bufs[2 * d + 1].out_atoms, out_bufs[2 * d + 1].cnt, type,
+ nbr2->rank, 2 * d + 1, comm );
+ }
+
+ if ( nbr1->atoms_cnt ) MPI_Wait( &req1, &stat1 );
+ if ( nbr2->atoms_cnt ) MPI_Wait( &req2, &stat2 );
+ }
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d dist: done\n", system->my_rank );
+#endif
+}
+
+
+void real_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
+{
+ int i;
+ real *in = (real*) dummy_in;
+ real *buf = (real*) dummy_buf;
+
+ for ( i = 0; i < out_buf->cnt; ++i )
+ buf[ out_buf->index[i] ] += in[i];
+}
+
+
+void rvec_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
+{
+ int i;
+ rvec *in = (rvec*) dummy_in;
+ rvec *buf = (rvec*) dummy_buf;
+
+ for ( i = 0; i < out_buf->cnt; ++i )
+ {
+ rvec_Add( buf[ out_buf->index[i] ], in[i] );
+#if defined(DEBUG)
+ fprintf( stderr, "rvec_unpacker: cnt=%d i =%d index[i]=%d\n",
+ out_buf->cnt, i, out_buf->index[i] );
+#endif
+ }
+}
+
+
+void rvec2_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
+{
+ int i;
+ rvec2 *in = (rvec2*) dummy_in;
+ rvec2 *buf = (rvec2*) dummy_buf;
+
+ for ( i = 0; i < out_buf->cnt; ++i )
+ {
+ buf[ out_buf->index[i] ][0] += in[i][0];
+ buf[ out_buf->index[i] ][1] += in[i][1];
+ }
+}
+
+
+void Coll( reax_system* system, mpi_datatypes *mpi_data,
+ void *buf, MPI_Datatype type, int scale, coll_unpacker unpack )
+{
+ int d;
+ void *in1, *in2;
+ mpi_out_data *out_bufs;
+ MPI_Comm comm;
+ MPI_Request req1, req2;
+ MPI_Status stat1, stat2;
+ neighbor_proc *nbr1, *nbr2;
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d coll: entered\n", system->my_rank );
+#endif
+ comm = mpi_data->comm_mesh3D;
+ in1 = mpi_data->in1_buffer;
+ in2 = mpi_data->in2_buffer;
+ out_bufs = mpi_data->out_buffers;
+
+ for ( d = 2; d >= 0; --d )
+ {
+ /* initiate recvs */
+ nbr1 = &(system->my_nbrs[2 * d]);
+ if ( out_bufs[2 * d].cnt )
+ MPI_Irecv(in1, out_bufs[2 * d].cnt, type, nbr1->rank, 2 * d + 1, comm, &req1);
+
+ nbr2 = &(system->my_nbrs[2 * d + 1]);
+ if ( out_bufs[2 * d + 1].cnt )
+ MPI_Irecv(in2, out_bufs[2 * d + 1].cnt, type, nbr2->rank, 2 * d, comm, &req2);
+
+ /* send both messages in dimension d */
+ if ( nbr1->atoms_cnt )
+ MPI_Send( buf + nbr1->atoms_str * scale, nbr1->atoms_cnt, type,
+ nbr1->rank, 2 * d, comm );
+
+ if ( nbr2->atoms_cnt )
+ MPI_Send( buf + nbr2->atoms_str * scale, nbr2->atoms_cnt, type,
+ nbr2->rank, 2 * d + 1, comm );
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d coll[%d] nbr1: str=%d cnt=%d recv=%d\n",
+ system->my_rank, d, nbr1->atoms_str, nbr1->atoms_cnt,
+ out_bufs[2 * d].cnt );
+ fprintf( stderr, "p%d coll[%d] nbr2: str=%d cnt=%d recv=%d\n",
+ system->my_rank, d, nbr2->atoms_str, nbr2->atoms_cnt,
+ out_bufs[2 * d + 1].cnt );
+#endif
+
+ if ( out_bufs[2 * d].cnt )
+ {
+ MPI_Wait( &req1, &stat1 );
+ unpack( in1, buf, out_bufs + (2 * d) );
+ }
+
+ if ( out_bufs[2 * d + 1].cnt )
+ {
+ MPI_Wait( &req2, &stat2 );
+ unpack( in2, buf, out_bufs + (2 * d + 1) );
+ }
+ }
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d coll: done\n", system->my_rank );
+#endif
+}
+#endif /*PURE_REAX*/
+
+/*****************************************************************************/
+real Parallel_Norm( real *v, int n, MPI_Comm comm )
+{
+ int i;
+ real my_sum, norm_sqr;
+
+ my_sum = 0;
+ for ( i = 0; i < n; ++i )
+ my_sum += SQR( v[i] );
+
+ MPI_Allreduce( &my_sum, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, comm );
+
+ return sqrt( norm_sqr );
+}
+
+
+
+real Parallel_Dot( real *v1, real *v2, int n, MPI_Comm comm )
+{
+ int i;
+ real my_dot, res;
+
+ my_dot = 0;
+ for ( i = 0; i < n; ++i )
+ my_dot += v1[i] * v2[i];
+
+ MPI_Allreduce( &my_dot, &res, 1, MPI_DOUBLE, MPI_SUM, comm );
+
+ return res;
+}
+
+
+
+real Parallel_Vector_Acc( real *v, int n, MPI_Comm comm )
+{
+ int i;
+ real my_acc, res;
+
+ my_acc = 0;
+ for ( i = 0; i < n; ++i )
+ my_acc += v[i];
+
+ MPI_Allreduce( &my_acc, &res, 1, MPI_DOUBLE, MPI_SUM, comm );
+
+ return res;
+}
+
+
+/*****************************************************************************/
+#if defined(TEST_FORCES)
+void Coll_ids_at_Master( reax_system *system, storage *workspace,
+ mpi_datatypes *mpi_data )
+{
+ int i;
+ int *id_list;
+
+ MPI_Gather( &system->n, 1, MPI_INT, workspace->rcounts, 1, MPI_INT,
+ MASTER_NODE, mpi_data->world );
+
+ if ( system->my_rank == MASTER_NODE )
+ {
+ workspace->displs[0] = 0;
+ for ( i = 1; i < system->wsize; ++i )
+ workspace->displs[i] = workspace->displs[i - 1] + workspace->rcounts[i - 1];
+ }
+
+ id_list = (int*) malloc( system->n * sizeof(int) );
+ for ( i = 0; i < system->n; ++i )
+ id_list[i] = system->my_atoms[i].orig_id;
+
+ MPI_Gatherv( id_list, system->n, MPI_INT,
+ workspace->id_all, workspace->rcounts, workspace->displs,
+ MPI_INT, MASTER_NODE, mpi_data->world );
+
+ sfree( id_list, "id_list" );
+
+#if defined(DEBUG)
+ if ( system->my_rank == MASTER_NODE )
+ {
+ for ( i = 0 ; i < system->bigN; ++i )
+ fprintf( stderr, "id_all[%d]: %d\n", i, workspace->id_all[i] );
+ }
+#endif
+}
+
+
+void Coll_rvecs_at_Master( reax_system *system, storage *workspace,
+ mpi_datatypes *mpi_data, rvec* v )
+{
+ MPI_Gatherv( v, system->n, mpi_data->mpi_rvec,
+ workspace->f_all, workspace->rcounts, workspace->displs,
+ mpi_data->mpi_rvec, MASTER_NODE, mpi_data->world );
+}
+
+#endif
diff --git a/PuReMD/src/bond_orders.c b/PuReMD/src/bond_orders.c
index cf6c69911f889bb528ef32c39e6542960ba0d7e0..286e15c04ab667e99e8b0533675f3be03c01ff36 100644
--- a/PuReMD/src/bond_orders.c
+++ b/PuReMD/src/bond_orders.c
@@ -667,12 +667,17 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
single_body_parameters *sbp_i, single_body_parameters *sbp_j,
two_body_parameters *twbp )
{
- int j, btop_j;
+ int j;
real r2, C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
- bond_data *ibond, *jbond;
- bond_order_data *bo_ij, *bo_ji;
+ bond_data *ibond;
+ bond_order_data *bo_ij;
+#if defined(HALF_LIST)
+ int btop_j;
+ bond_data *jbond;
+ bond_order_data *bo_ji;
+#endif
j = nbr_pj->nbr;
r2 = SQR(nbr_pj->d);
@@ -705,29 +710,46 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
{
/****** bonds i-j and j-i ******/
ibond = &( bonds->bond_list[btop_i] );
+#if defined(HALF_LIST)
btop_j = End_Index( j, bonds );
jbond = &(bonds->bond_list[btop_j]);
+#endif
+#if defined(HALF_LIST)
ibond->nbr = j;
- jbond->nbr = i;
ibond->d = nbr_pj->d;
- jbond->d = nbr_pj->d;
rvec_Copy( ibond->dvec, nbr_pj->dvec );
- rvec_Scale( jbond->dvec, -1, nbr_pj->dvec );
ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
- ivec_Scale( jbond->rel_box, -1, nbr_pj->rel_box );
ibond->dbond_index = btop_i;
- jbond->dbond_index = btop_i;
ibond->sym_index = btop_j;
+ jbond->nbr = i;
+ jbond->d = nbr_pj->d;
+ rvec_Scale( jbond->dvec, -1.0, nbr_pj->dvec );
+ ivec_Scale( jbond->rel_box, -1.0, nbr_pj->rel_box );
+ jbond->dbond_index = btop_i;
jbond->sym_index = btop_i;
- Set_End_Index( j, btop_j + 1, bonds );
+ Set_End_Index( j, btop_j + 1, bonds );
+#else
+ ibond->nbr = j;
+ ibond->d = nbr_pj->d;
+ rvec_Copy( ibond->dvec, nbr_pj->dvec );
+ ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
+ ibond->dbond_index = btop_i;
+#endif
+
bo_ij = &( ibond->bo_data );
+ bo_ij->BO = BO;
+ bo_ij->BO_s = BO_s;
+ bo_ij->BO_pi = BO_pi;
+ bo_ij->BO_pi2 = BO_pi2;
+#if defined(HALF_LIST)
bo_ji = &( jbond->bo_data );
- bo_ji->BO = bo_ij->BO = BO;
- bo_ji->BO_s = bo_ij->BO_s = BO_s;
- bo_ji->BO_pi = bo_ij->BO_pi = BO_pi;
- bo_ji->BO_pi2 = bo_ij->BO_pi2 = BO_pi2;
+ bo_ji->BO = BO;
+ bo_ji->BO_s = BO_s;
+ bo_ji->BO_pi = BO_pi;
+ bo_ji->BO_pi2 = BO_pi2;
+#endif
/* Bond Order page2-3, derivative of total bond order prime */
Cln_BOp_s = twbp->p_bo2 * C12 / r2;
@@ -736,62 +758,51 @@ int BOp( storage *workspace, reax_list *bonds, real bo_cut,
/* Only dln_BOp_xx wrt. dr_i is stored here, note that
dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
- rvec_Scale(bo_ij->dln_BOp_s, -bo_ij->BO_s * Cln_BOp_s, ibond->dvec);
- rvec_Scale(bo_ij->dln_BOp_pi, -bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec);
- rvec_Scale(bo_ij->dln_BOp_pi2,
- -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec);
- rvec_Scale(bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s);
- rvec_Scale(bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
- rvec_Scale(bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
+#if defined(HALF_LIST)
+ rvec_Scale(bo_ij->dln_BOp_s, -1.0 * bo_ij->BO_s * Cln_BOp_s, ibond->dvec);
+ rvec_Scale(bo_ij->dln_BOp_pi, -1.0 * bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec);
+ rvec_Scale(bo_ij->dln_BOp_pi2, -1.0 * bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec);
+ rvec_Scale(bo_ji->dln_BOp_s, -1.0, bo_ij->dln_BOp_s);
+ rvec_Scale(bo_ji->dln_BOp_pi, -1.0, bo_ij->dln_BOp_pi );
+ rvec_Scale(bo_ji->dln_BOp_pi2, -1.0, bo_ij->dln_BOp_pi2 );
+#else
+ rvec_Scale(bo_ij->dln_BOp_s, -1.0 * bo_ij->BO_s * Cln_BOp_s, ibond->dvec);
+ rvec_Scale(bo_ij->dln_BOp_pi, -1.0 * bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec);
+ rvec_Scale(bo_ij->dln_BOp_pi2, -1.0 * bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec);
+#endif
/* Only dBOp wrt. dr_i is stored here, note that
dBOp/dr_i = -dBOp/dr_j and all others are 0 */
- rvec_Scale( bo_ij->dBOp,
- -(bo_ij->BO_s * Cln_BOp_s +
- bo_ij->BO_pi * Cln_BOp_pi +
- bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
- rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );
+#if defined(HALF_LIST)
+ rvec_Scale( bo_ij->dBOp, -1.0 * (bo_ij->BO_s * Cln_BOp_s
+ + bo_ij->BO_pi * Cln_BOp_pi
+ + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
+ rvec_Scale( bo_ji->dBOp, -1.0, bo_ij->dBOp );
+#else
+ rvec_Scale( bo_ij->dBOp, -1.0 * (bo_ij->BO_s * Cln_BOp_s
+ + bo_ij->BO_pi * Cln_BOp_pi
+ + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
+#endif
rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
+#if defined(HALF_LIST)
rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );
+#endif
bo_ij->BO_s -= bo_cut;
bo_ij->BO -= bo_cut;
+ workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
+ bo_ij->Cdbo = 0.0;
+ bo_ij->Cdbopi = 0.0;
+ bo_ij->Cdbopi2 = 0.0;
+#if defined(HALF_LIST)
bo_ji->BO_s -= bo_cut;
bo_ji->BO -= bo_cut;
- workspace->total_bond_order[i] += bo_ij->BO; //currently total_BOp
workspace->total_bond_order[j] += bo_ji->BO; //currently total_BOp
- bo_ij->Cdbo = bo_ij->Cdbopi = bo_ij->Cdbopi2 = 0.0;
- bo_ji->Cdbo = bo_ji->Cdbopi = bo_ji->Cdbopi2 = 0.0;
-
- /*fprintf( stderr, "%d %d %g %g %g\n",
- i+1, j+1, bo_ij->BO, bo_ij->BO_pi, bo_ij->BO_pi2 );*/
-
- /*fprintf( stderr, "Cln_BOp_s: %f, pbo2: %f, C12:%f\n",
- Cln_BOp_s, twbp->p_bo2, C12 );
- fprintf( stderr, "Cln_BOp_pi: %f, pbo4: %f, C34:%f\n",
- Cln_BOp_pi, twbp->p_bo4, C34 );
- fprintf( stderr, "Cln_BOp_pi2: %f, pbo6: %f, C56:%f\n",
- Cln_BOp_pi2, twbp->p_bo6, C56 );*/
- /*fprintf(stderr, "pbo1: %f, pbo2:%f\n", twbp->p_bo1, twbp->p_bo2);
- fprintf(stderr, "pbo3: %f, pbo4:%f\n", twbp->p_bo3, twbp->p_bo4);
- fprintf(stderr, "pbo5: %f, pbo6:%f\n", twbp->p_bo5, twbp->p_bo6);
- fprintf( stderr, "r_s: %f, r_p: %f, r_pp: %f\n",
- twbp->r_s, twbp->r_p, twbp->r_pp );
- fprintf( stderr, "C12: %g, C34:%g, C56:%g\n", C12, C34, C56 );*/
-
- /*fprintf( stderr, "\tfactors: %g %g %g\n",
- -(bo_ij->BO_s * Cln_BOp_s + bo_ij->BO_pi * Cln_BOp_pi +
- bo_ij->BO_pi2 * Cln_BOp_pp),
- -bo_ij->BO_pi * Cln_BOp_pi, -bo_ij->BO_pi2 * Cln_BOp_pi2 );*/
- /*fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
- bo_ij->dBOp[0], bo_ij->dBOp[1], bo_ij->dBOp[2] );
- fprintf( stderr, "dBOpi:\t[%g, %g, %g]\n",
- bo_ij->dln_BOp_pi[0], bo_ij->dln_BOp_pi[1],
- bo_ij->dln_BOp_pi[2] );
- fprintf( stderr, "dBOpi2:\t[%g, %g, %g]\n\n",
- bo_ij->dln_BOp_pi2[0], bo_ij->dln_BOp_pi2[1],
- bo_ij->dln_BOp_pi2[2] );*/
+ bo_ji->Cdbo = 0.0;
+ bo_ji->Cdbopi = 0.0;
+ bo_ji->Cdbopi2 = 0.0;
+#endif
return 1;
}
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index d438a0eba65d50ae11a8f6bd03fa5fe79492474b..d157eb6f913ed05aae59673e128834249cd52825 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -326,8 +326,8 @@ void Init_Forces( reax_system *system, control_params *control,
int i, j, pj;
int start_i, end_i;
int type_i, type_j;
- int Htop, btop_i, btop_j, num_bonds, num_hbonds;
- int ihb, jhb, ihb_top, jhb_top;
+ int Htop, btop_i, num_bonds, num_hbonds;
+ int ihb, jhb, ihb_top;
int local, flag, renbr;
real r_ij, cutoff;
sparse_matrix *H;
@@ -336,6 +336,12 @@ void Init_Forces( reax_system *system, control_params *control,
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
+#if defined(HALF_LIST)
+ int jhb_top;
+#else
+ int start_j, end_j;
+ int btop_j;
+#endif
far_nbrs = lists[FAR_NBRS];
bonds = lists[BONDS];
@@ -354,7 +360,7 @@ void Init_Forces( reax_system *system, control_params *control,
Htop = 0;
num_bonds = 0;
num_hbonds = 0;
- btop_i = btop_j = 0;
+ btop_i = 0;
renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
for ( i = 0; i < system->N; ++i )
@@ -363,7 +369,13 @@ void Init_Forces( reax_system *system, control_params *control,
type_i = atom_i->type;
start_i = Start_Index(i, far_nbrs);
end_i = End_Index(i, far_nbrs);
+#if defined(HALF_LIST)
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
btop_i = End_Index( i, bonds );
+#else
+ btop_i = Start_Index( i, bonds );
+#endif
sbp_i = &(system->reax_param.sbp[type_i]);
if ( i < system->n )
@@ -390,8 +402,19 @@ void Init_Forces( reax_system *system, control_params *control,
{
ihb = sbp_i->p_hbond;
if ( ihb == 1 )
+ {
+#if defined(HALF_LIST)
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
ihb_top = End_Index( atom_i->Hindex, hbonds );
- else ihb_top = -1;
+#else
+ ihb_top = Start_Index( atom_i->Hindex, hbonds );
+#endif
+ }
+ else
+ {
+ ihb_top = -1;
+ }
}
}
@@ -401,11 +424,6 @@ void Init_Forces( reax_system *system, control_params *control,
nbr_pj = &( far_nbrs->far_nbr_list[pj] );
j = nbr_pj->nbr;
atom_j = &(system->my_atoms[j]);
- //fprintf( stderr, "%d%d i=%d x_i: %f %f %f,j=%d x_j: %f %f %f, d=%f\n",
- // MIN(atom_i->orig_id, atom_j->orig_id),
- // MAX(atom_i->orig_id, atom_j->orig_id),
- // i, atom_i->x[0], atom_i->x[1], atom_i->x[2],
- // j, atom_j->x[0], atom_j->x[1], atom_j->x[2], nbr_pj->d );
if ( renbr )
{
if (nbr_pj->d <= cutoff)
@@ -439,14 +457,11 @@ void Init_Forces( reax_system *system, control_params *control,
if ( local )
{
/* H matrix entry */
+#if defined(HALF_LIST)
if ( j < system->n || atom_i->orig_id < atom_j->orig_id ) //tryQEq||1
+#endif
{
H->entries[Htop].j = j;
- //fprintf( stdout, "%d%d %d %d\n",
- // MIN(atom_i->orig_id, atom_j->orig_id),
- // MAX(atom_i->orig_id, atom_j->orig_id),
- // MIN(atom_i->orig_id, atom_j->orig_id),
- // MAX(atom_i->orig_id, atom_j->orig_id) );
if ( control->tabulate == 0 )
H->entries[Htop].val = Compute_H(r_ij, twbp->gamma, workspace->Tap);
else H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
@@ -467,6 +482,8 @@ void Init_Forces( reax_system *system, control_params *control,
++ihb_top;
++num_hbonds;
}
+#if defined(HALF_LIST)
+ /* only add to list for local j (far nbrs is half-list) */
else if ( j < system->n && ihb == 2 && jhb == 1 )
{
jhb_top = End_Index( atom_j->Hindex, hbonds );
@@ -476,6 +493,7 @@ void Init_Forces( reax_system *system, control_params *control,
Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
++num_hbonds;
}
+#endif
}
}
@@ -493,11 +511,7 @@ void Init_Forces( reax_system *system, control_params *control,
else if ( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 )
{
workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
- //if( workspace->bond_mark[i] == 1000 )
- // workspace->done_after[i] = pj;
}
- //fprintf( stdout, "%d%d - %d(%d) %d(%d)\n",
- // i , j, i, workspace->bond_mark[i], j, workspace->bond_mark[j] );
}
}
}
@@ -511,53 +525,30 @@ void Init_Forces( reax_system *system, control_params *control,
}
}
- //fprintf( stderr, "after the first init loop\n" );
- /*for( i = system->n; i < system->N; ++i )
- if( workspace->bond_mark[i] > 3 ) {
- start_i = Start_Index(i, bonds);
- end_i = End_Index(i, bonds);
- num_bonds -= (end_i - start_i);
- Set_End_Index(i, start_i, bonds );
- }*/
-
- /*for( i = system->n; i < system->N; ++i ) {
- start_i = Start_Index(i, far_nbrs);
- end_i = workspace->done_after[i];
+#if !defined(HALF_LIST)
+ /* set sym_index for bonds list (far_nbrs full list) */
+ for ( i = 0; i < system->N; ++i )
+ {
+ start_i = Start_Index( i, bonds );
+ end_i = End_Index( i, bonds );
- if( workspace->bond_mark[i] >= 2 && start_i < end_i ) {
- atom_i = &(system->my_atoms[i]);
- type_i = atom_i->type;
- btop_i = End_Index( i, bonds );
- sbp_i = &(system->reax_param.sbp[type_i]);
+ for ( btop_i = start_i; btop_i < end_i; ++btop_i )
+ {
+ j = bonds->bond_list[btop_i].nbr;
+ start_j = Start_Index( j, bonds );
+ end_j = End_Index( j, bonds );
- for( pj = start_i; pj < end_i; ++pj ) {
- nbr_pj = &( far_nbrs->far_nbr_list[pj] );
- j = nbr_pj->nbr;
-
- if( workspace->bond_mark[j] >= 2 && nbr_pj->d <= control->bond_cut ) {
- atom_j = &(system->my_atoms[j]);
- type_j = atom_j->type;
- sbp_j = &(system->reax_param.sbp[type_j]);
- twbp = &(system->reax_param.tbp[type_i][type_j]);
-
- if( BOp( workspace, bonds, control->bo_cut,
- i , btop_i, nbr_pj, sbp_i, sbp_j, twbp ) ) {
- num_bonds += 2;
- ++btop_i;
-
- if( workspace->bond_mark[j] > workspace->bond_mark[i] + 1 )
- workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
- else if( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 )
- workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
-
- //fprintf( stdout, "%d%d - %d(%d) %d(%d) new\n",
- // i , j, i, workspace->bond_mark[i], j, workspace->bond_mark[j] );
- }
- }
+ for ( btop_j = start_j; btop_j < end_j; ++btop_j )
+ {
+ if ( bonds->bond_list[btop_j].nbr == i )
+ {
+ bonds->bond_list[btop_i].sym_index = btop_j;
+ break;
+ }
+ }
}
- Set_End_Index( i, btop_i, bonds );
- }
- }*/
+ }
+#endif
workspace->realloc.Htop = Htop;
workspace->realloc.num_bonds = num_bonds;
@@ -595,8 +586,8 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
int i, j, pj;
int start_i, end_i;
int type_i, type_j;
- int btop_i, btop_j, num_bonds, num_hbonds;
- int ihb, jhb, ihb_top, jhb_top;
+ int btop_i, num_bonds, num_hbonds;
+ int ihb, jhb, ihb_top;
int local, flag, renbr;
real r_ij, cutoff;
reax_list *far_nbrs, *bonds, *hbonds;
@@ -604,6 +595,12 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
reax_atom *atom_i, *atom_j;
+#if defined(HALF_LIST)
+ int jhb_top;
+#else
+ int start_j, end_j;
+ int btop_j;
+#endif
far_nbrs = lists[FAR_NBRS];
bonds = lists[BONDS];
@@ -619,7 +616,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
num_bonds = 0;
num_hbonds = 0;
- btop_i = btop_j = 0;
+ btop_i = 0;
renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
for ( i = 0; i < system->N; ++i )
@@ -628,7 +625,13 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
type_i = atom_i->type;
start_i = Start_Index(i, far_nbrs);
end_i = End_Index(i, far_nbrs);
+#if defined(HALF_LIST)
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
btop_i = End_Index( i, bonds );
+#else
+ btop_i = Start_Index( i, bonds );
+#endif
sbp_i = &(system->reax_param.sbp[type_i]);
if ( i < system->n )
@@ -648,8 +651,19 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
{
ihb = sbp_i->p_hbond;
if ( ihb == 1 )
+ {
+#if defined(HALF_LIST)
+ /* start at end because other atoms
+ * can add to this atom's list (half-list) */
ihb_top = End_Index( atom_i->Hindex, hbonds );
- else ihb_top = -1;
+#else
+ ihb_top = Start_Index( atom_i->Hindex, hbonds );
+#endif
+ }
+ else
+ {
+ ihb_top = -1;
+ }
}
/* update i-j distance - check if j is within cutoff */
@@ -705,6 +719,8 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
++ihb_top;
++num_hbonds;
}
+#if defined(HALF_LIST)
+ /* only add to list for local j (far nbrs is half-list) */
else if ( j < system->n && ihb == 2 && jhb == 1 )
{
jhb_top = End_Index( atom_j->Hindex, hbonds );
@@ -714,6 +730,7 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
Set_End_Index( atom_j->Hindex, jhb_top + 1, hbonds );
++num_hbonds;
}
+#endif
}
}
@@ -746,13 +763,30 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
}
- /*for( i = system->n; i < system->N; ++i )
- if( workspace->bond_mark[i] > 3 ) {
- start_i = Start_Index(i, bonds);
- end_i = End_Index(i, bonds);
- num_bonds -= (end_i - start_i);
- Set_End_Index(i, start_i, bonds );
- }*/
+#if !defined(HALF_LIST)
+ /* set sym_index for bonds list (far_nbrs full list) */
+ for ( i = 0; i < system->N; ++i )
+ {
+ start_i = Start_Index( i, bonds );
+ end_i = End_Index( i, bonds );
+
+ for ( btop_i = start_i; btop_i < end_i; ++btop_i )
+ {
+ j = bonds->bond_list[btop_i].nbr;
+ start_j = Start_Index( j, bonds );
+ end_j = End_Index( j, bonds );
+
+ for ( btop_j = start_j; btop_j < end_j; ++btop_j )
+ {
+ if ( bonds->bond_list[btop_j].nbr == i )
+ {
+ bonds->bond_list[btop_i].sym_index = btop_j;
+ break;
+ }
+ }
+ }
+ }
+#endif
workspace->realloc.num_bonds = num_bonds;
workspace->realloc.num_hbonds = num_hbonds;
@@ -789,7 +823,10 @@ void Estimate_Storages( reax_system *system, control_params *control,
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
- reax_atom *atom_i, *atom_j;
+ reax_atom *atom_i;
+#if defined(HALF_LIST)
+ reax_atom *atom_j;
+#endif
far_nbrs = lists[FAR_NBRS];
*Htop = 0;
@@ -823,7 +860,9 @@ void Estimate_Storages( reax_system *system, control_params *control,
{
nbr_pj = &( far_nbrs->far_nbr_list[pj] );
j = nbr_pj->nbr;
+#if defined(HALF_LIST)
atom_j = &(system->my_atoms[j]);
+#endif
if (nbr_pj->d <= cutoff)
{
@@ -834,8 +873,12 @@ void Estimate_Storages( reax_system *system, control_params *control,
if ( local )
{
+#if defined(HALF_LIST)
if ( j < system->n || atom_i->orig_id < atom_j->orig_id ) //tryQEq ||1
+#endif
+ {
++(*Htop);
+ }
/* hydrogen bond lists */
if ( control->hbond_cut > 0.1 && (ihb == 1 || ihb == 2) &&
@@ -843,9 +886,16 @@ void Estimate_Storages( reax_system *system, control_params *control,
{
jhb = sbp_j->p_hbond;
if ( ihb == 1 && jhb == 2 )
+ {
++hb_top[i];
+ }
+ /* only add to list for local j (far nbrs is half-list) */
+#if defined(HALF_LIST)
else if ( j < system->n && ihb == 2 && jhb == 1 )
+ {
++hb_top[j];
+ }
+#endif
}
}
@@ -881,7 +931,9 @@ void Estimate_Storages( reax_system *system, control_params *control,
if ( BO >= control->bo_cut )
{
++bond_top[i];
+#if defined(HALF_LIST)
++bond_top[j];
+#endif
}
}
}
diff --git a/PuReMD/src/geo_tools.c b/PuReMD/src/geo_tools.c
index c1e3549fedf2039f96af84aa263cd114c7dee3cb..480b9ae7197ec53b839df411f2ba23262db5d0a8 100644
--- a/PuReMD/src/geo_tools.c
+++ b/PuReMD/src/geo_tools.c
@@ -57,7 +57,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d@count atoms:\n", system->my_rank );
- fprintf( stderr, "p%d: bigN = %d\n", system->my_rank, system->bigN );
+ fprintf( stderr, "p%d: bigNNN = %d\n", system->my_rank, system->bigN );
fprintf( stderr, "p%d: n = %d\n", system->my_rank, system->n );
fprintf( stderr, "p%d: N = %d\n\n", system->my_rank, system->N );
#endif
@@ -241,7 +241,7 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
//#if defined(DEBUG)
fprintf( stderr, "p%d@count atoms:\n", system->my_rank );
- fprintf( stderr, "p%d: bigN = %d\n", system->my_rank, system->bigN );
+ fprintf( stderr, "p%d: bigNNN = %d\n", system->my_rank, system->bigN );
fprintf( stderr, "p%d: n = %d\n", system->my_rank, system->n );
fprintf( stderr, "p%d: N = %d\n\n", system->my_rank, system->N );
//#endif
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index 541a132be7bc18354069022785af438030c2e286..1aaa18dffa127e7589c3bd8253aba478c3f1e277 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -44,10 +44,16 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
for ( i = 0; i < A->n; ++i )
{
si = A->start[i];
+#if defined(HALF_LIST)
b[i][0] += A->entries[si].val * x[i][0];
b[i][1] += A->entries[si].val * x[i][1];
+#endif
+#if defined(HALF_LIST)
for ( k = si + 1; k < A->end[i]; ++k )
+#else
+ for ( k = si; k < A->end[i]; ++k )
+#endif
{
j = A->entries[k].j;
H = A->entries[k].val;
@@ -55,11 +61,13 @@ void dual_Sparse_MatVec( sparse_matrix *A, rvec2 *x, rvec2 *b, int N )
b[i][0] += H * x[j][0];
b[i][1] += H * x[j][1];
+#if defined(HALF_LIST)
// comment out for tryQEq
//if( j < A->n ) {
b[j][0] += H * x[i][0];
b[j][1] += H * x[i][1];
//}
+#endif
}
}
}
@@ -91,8 +99,10 @@ int dual_CG( reax_system *system, storage *workspace, sparse_matrix *H,
Dist( system, mpi_data, x, mpi_data->mpi_rvec2, scale, rvec2_packer );
dual_Sparse_MatVec( H, x, workspace->q2, N );
+#if defined(HALF_LIST)
// tryQEq
Coll(system, mpi_data, workspace->q2, mpi_data->mpi_rvec2, scale, rvec2_unpacker);
+#endif
#if defined(CG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
@@ -142,8 +152,10 @@ int dual_CG( reax_system *system, storage *workspace, sparse_matrix *H,
{
Dist(system, mpi_data, workspace->d2, mpi_data->mpi_rvec2, scale, rvec2_packer);
dual_Sparse_MatVec( H, workspace->d2, workspace->q2, N );
+#if defined(HALF_LIST)
// tryQEq
Coll(system, mpi_data, workspace->q2, mpi_data->mpi_rvec2, scale, rvec2_unpacker);
+#endif
#if defined(CG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
@@ -265,14 +277,25 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
for ( i = 0; i < A->n; ++i )
{
si = A->start[i];
+
+#if defined(HALF_LIST)
b[i] += A->entries[si].val * x[i];
+#endif
+
+#if defined(HALF_LIST)
for ( k = si + 1; k < A->end[i]; ++k )
+#else
+ for ( k = si; k < A->end[i]; ++k )
+#endif
{
j = A->entries[k].j;
H = A->entries[k].val;
+
b[i] += H * x[j];
+#if defined(HALF_LIST)
//if( j < A->n ) // comment out for tryQEq
b[j] += H * x[i];
+#endif
}
}
}
@@ -322,8 +345,10 @@ int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
scale = sizeof(real) / sizeof(void);
Dist( system, mpi_data, x, MPI_DOUBLE, scale, real_packer );
Sparse_MatVec( H, x, workspace->q, system->N );
+#if defined(HALF_LIST)
// tryQEq
Coll( system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker );
+#endif
#if defined(CG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
@@ -356,8 +381,10 @@ int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
{
Dist( system, mpi_data, workspace->d, MPI_DOUBLE, scale, real_packer );
Sparse_MatVec( H, workspace->d, workspace->q, system->N );
+#if defined(HALF_LIST)
//tryQEq
Coll(system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker);
+#endif
#if defined(CG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
diff --git a/PuReMD/src/my* b/PuReMD/src/my*
new file mode 100644
index 0000000000000000000000000000000000000000..6ad2fd10f961011b33e2ef088343f431752b6da6
--- /dev/null
+++ b/PuReMD/src/my*
@@ -0,0 +1,67 @@
+./integrate.h:1:/*----------------------------------------------------------------------
+./reset_tools.h:1:/*----------------------------------------------------------------------
+./allocate.c:1:/*----------------------------------------------------------------------
+./lookup.h:1:/*----------------------------------------------------------------------
+./restart.c:1:/*----------------------------------------------------------------------
+./box.c:1:/*----------------------------------------------------------------------
+./torsion_angles.h:1:/*----------------------------------------------------------------------
+./bond_orders.h:1:/*----------------------------------------------------------------------
+./reax_types.h:1:/*----------------------------------------------------------------------
+./grid.c:1:/*----------------------------------------------------------------------
+./io_tools.h:1:/*----------------------------------------------------------------------
+./nonbonded.h:1:/*----------------------------------------------------------------------
+./reax_defs.h:1:/*----------------------------------------------------------------------
+./basic_comm.h:1:/*----------------------------------------------------------------------
+./multi_body.h:1:/*----------------------------------------------------------------------
+./linear_solvers.h:1:/*----------------------------------------------------------------------
+./neighbors.c:1:/*----------------------------------------------------------------------
+./valence_angles.h:1:/*----------------------------------------------------------------------
+./tool_box.c:1:/*----------------------------------------------------------------------
+./analyze.h:1:/*----------------------------------------------------------------------
+./qEq.c:1:/*----------------------------------------------------------------------
+./control.c:1:/*----------------------------------------------------------------------
+./vector.c:1:/*----------------------------------------------------------------------
+./bonds.h:1:/*----------------------------------------------------------------------
+./forces.c:1:/*----------------------------------------------------------------------
+./system_props.c:1:/*----------------------------------------------------------------------
+./traj.h:1:/*----------------------------------------------------------------------
+./comm_tools.h:1:/*----------------------------------------------------------------------
+./ffield.h:1:/*----------------------------------------------------------------------
+./geo_tools.c:1:/*----------------------------------------------------------------------
+./hydrogen_bonds.h:1:/*----------------------------------------------------------------------
+./init_md.h:1:/*----------------------------------------------------------------------
+./list.h:1:/*----------------------------------------------------------------------
+./random.c:1:/*----------------------------------------------------------------------
+./linear_solvers.c:1:/*----------------------------------------------------------------------
+./neighbors.h:1:/*----------------------------------------------------------------------
+./valence_angles.c:1:/*----------------------------------------------------------------------
+./tool_box.h:1:/*----------------------------------------------------------------------
+./analyze.c:1:/*----------------------------------------------------------------------
+./nonbonded.c:1:/*----------------------------------------------------------------------
+./basic_comm.c:1:/*----------------------------------------------------------------------
+./multi_body.c:1:/*----------------------------------------------------------------------
+./restart.h:1:/*----------------------------------------------------------------------
+./box.h:1:/*----------------------------------------------------------------------
+./parallelreax.c:1:/*----------------------------------------------------------------------
+./torsion_angles.c:1:/*----------------------------------------------------------------------
+./bond_orders.c:1:/*----------------------------------------------------------------------
+./grid.h:1:/*----------------------------------------------------------------------
+./io_tools.c:1:/*----------------------------------------------------------------------
+./integrate.c:1:/*----------------------------------------------------------------------
+./reset_tools.c:1:/*----------------------------------------------------------------------
+./allocate.h:1:/*----------------------------------------------------------------------
+./lookup.c:1:/*----------------------------------------------------------------------
+./init_md.c:1:/*----------------------------------------------------------------------
+./list.c:1:/*----------------------------------------------------------------------
+./random.h:1:/*----------------------------------------------------------------------
+./geo_tools.h:1:/*----------------------------------------------------------------------
+./hydrogen_bonds.c:1:/*----------------------------------------------------------------------
+./forces.h:1:/*----------------------------------------------------------------------
+./system_props.h:1:/*----------------------------------------------------------------------
+./traj.c:1:/*----------------------------------------------------------------------
+./comm_tools.c:1:/*----------------------------------------------------------------------
+./ffield.c:1:/*----------------------------------------------------------------------
+./qEq.h:1:/*----------------------------------------------------------------------
+./control.h:1:/*----------------------------------------------------------------------
+./vector.h:1:/*----------------------------------------------------------------------
+./bonds.c:1:/*----------------------------------------------------------------------
diff --git a/PuReMD/src/neighbors.c b/PuReMD/src/neighbors.c
index 5be0016eda64a6e4ec56f02d77695cd212fcdd55..fa86a7381404061d9b8a28a76c59fd4b492f0040 100644
--- a/PuReMD/src/neighbors.c
+++ b/PuReMD/src/neighbors.c
@@ -91,7 +91,9 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
/* first pick up a cell in the grid */
for ( i = 0; i < g->ncells[0]; i++ )
+ {
for ( j = 0; j < g->ncells[1]; j++ )
+ {
for ( k = 0; k < g->ncells[2]; k++ )
{
gci = &(g->cells[i][j][k]);
@@ -108,11 +110,23 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
itr = 0;
while ( (gcj = gci->nbrs[itr]) != NULL )
{
- if ( gci->str <= gcj->str &&
+ if (
+#if defined(HALF_LIST)
+ gci->str <= gcj->str &&
+#endif
(DistSqr_to_Special_Point(gci->nbrs_cp[itr], atom1->x) <= cutoff) )
+ {
/* pick up another atom from the neighbor cell */
for ( m = gcj->str; m < gcj->end; ++m )
- if ( l < m ) // prevent recounting same pairs within a gcell
+ {
+#if defined(HALF_LIST)
+ /* prevent recounting same pairs within a gcell and
+ * make half-list */
+ if ( l < m )
+#else
+ /* prevent recounting same pairs within a gcell */
+ if ( l != m )
+#endif
{
atom2 = &(system->my_atoms[m]);
dvec[0] = atom2->x[0] - atom1->x[0];
@@ -125,20 +139,22 @@ void Generate_Neighbor_Lists( reax_system *system, simulation_data *data,
nbr_data->nbr = m;
nbr_data->d = sqrt(d);
rvec_Copy( nbr_data->dvec, dvec );
- //ivec_Copy( nbr_data->rel_box, gcj->rel_box );
ivec_ScaledSum( nbr_data->rel_box,
1, gcj->rel_box, -1, gci->rel_box );
++num_far;
}
}
+ }
+ }
+
++itr;
}
+
Set_End_Index( l, num_far, far_nbrs );
- //fprintf(stderr, "i:%d, start: %d, end: %d - itr: %d\n",
- // atom1,Start_Index(atom1,far_nbrs),End_Index(atom1,far_nbrs),
- // itr);
}
}
+ }
+ }
workspace->realloc.num_far = num_far;
@@ -180,7 +196,9 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists )
/* first pick up a cell in the grid */
for ( i = 0; i < g->ncells[0]; i++ )
+ {
for ( j = 0; j < g->ncells[1]; j++ )
+ {
for ( k = 0; k < g->ncells[2]; k++ )
{
gci = &(g->cells[i][j][k]);
@@ -196,12 +214,23 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists )
itr = 0;
while ( (gcj = gci->nbrs[itr]) != NULL )
{
- if (gci->str <= gcj->str &&
+ if (
+#if defined(HALF_LIST)
+ gci->str <= gcj->str &&
+#endif
(DistSqr_to_Special_Point(gci->nbrs_cp[itr], atom1->x) <= cutoff))
- //fprintf( stderr, "\t\tgcell2: %d\n", itr );
+ {
/* pick up another atom from the neighbor cell */
for ( m = gcj->str; m < gcj->end; ++m )
+ {
+#if defined(HALF_LIST)
+ /* prevent recounting same pairs within a gcell and
+ * make half-list */
if ( l < m )
+#else
+ /* prevent recounting same pairs within a gcell */
+ if ( l != m )
+#endif
{
//fprintf( stderr, "\t\t\tatom2=%d\n", m );
atom2 = &(system->my_atoms[m]);
@@ -212,13 +241,15 @@ int Estimate_NumNeighbors( reax_system *system, reax_list **lists )
if ( d <= cutoff )
++num_far;
}
+ }
+ }
++itr;
}
- //fprintf( stderr, "itr: %d, tested: %d, num_nbrs: %d\n",
- // itr, tested, num_far-tmp );
}
}
+ }
+ }
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: estimate nbrs done - num_far=%d\n",
diff --git a/PuReMD/src/nonbonded.c b/PuReMD/src/nonbonded.c
index ab25b807d5e0a0263a4d8f89ecc475c8615b8bdd..bdcd80792597b02d120d5c587169d80e2efefcb4 100644
--- a/PuReMD/src/nonbonded.c
+++ b/PuReMD/src/nonbonded.c
@@ -71,7 +71,11 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
j = nbr_pj->nbr;
orig_j = system->my_atoms[j].orig_id;
+#if defined(HALF_LIST)
if ( nbr_pj->d <= control->nonb_cut && (j < natoms || orig_i < orig_j) )
+#else
+ if ( nbr_pj->d <= control->nonb_cut && orig_i < orig_j )
+#endif
{
r_ij = nbr_pj->d;
twbp = &(system->reax_param.tbp[ system->my_atoms[i].type ]
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 152c0c54c35638a06daa415813d34e7c97e8629c..0e58175e24aae6f9962e7f00529d7302ee9387d4 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -39,6 +39,7 @@
//#define DEBUG
//#define DEBUG_FOCUS
+//#define HALF_LIST
//#define TEST_ENERGY
//#define TEST_FORCES
//#define CG_PERFORMANCE