From da7b9380bb9be55a5988d13f1f2d314f45f8fd54 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Wed, 15 Aug 2018 12:53:32 -0400
Subject: [PATCH] sPuReMD: backport efforts for SMALL_BOX_SUPPORT from OGOLEM
version of code and conditionally compile. Remove unused code. Other code
clean-up.
---
sPuReMD/src/analyze.c | 27 +-
sPuReMD/src/bond_orders.c | 36 --
sPuReMD/src/box.c | 51 +--
sPuReMD/src/lookup.c | 114 ++-----
sPuReMD/src/lookup.h | 7 -
sPuReMD/src/neighbors.c | 688 +++++++++++++++++---------------------
sPuReMD/src/print_utils.c | 25 --
sPuReMD/src/reax_types.h | 171 +++++-----
sPuReMD/src/reset_utils.c | 1 +
9 files changed, 458 insertions(+), 662 deletions(-)
diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index eba15b97..f51212a5 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -29,7 +29,7 @@
#include "vector.h"
-#define MAX_FRAGMENT_TYPES 100
+#define MAX_FRAGMENT_TYPES (100)
enum atoms
@@ -54,9 +54,9 @@ typedef struct
} molecule;
-// copy bond list into old bond list
+/* copy bond list into old bond list */
static void Copy_Bond_List( reax_system *system, control_params *control,
- reax_list **lists )
+ reax_list **lists )
{
int i, j, top_old;
reax_list *new_bonds = lists[BONDS];
@@ -111,13 +111,17 @@ static int Compare_Bond_Lists( int atom, control_params *control, reax_list **li
oldp = MIN( oldp + 1, End_Index( atom, old_bonds ) ),
newp = MIN( newp + 1, End_Index( atom, new_bonds ) ) )
{
- while ( oldp < End_Index( atom, old_bonds ) &&
- old_bonds->select.bond_list[oldp].bo_data.BO < control->bg_cut )
+ while ( oldp < End_Index( atom, old_bonds )
+ && old_bonds->select.bond_list[oldp].bo_data.BO < control->bg_cut )
+ {
++oldp;
+ }
- while ( newp < End_Index( atom, new_bonds ) &&
- new_bonds->select.bond_list[newp].bo_data.BO < control->bg_cut )
+ while ( newp < End_Index( atom, new_bonds )
+ && new_bonds->select.bond_list[newp].bo_data.BO < control->bg_cut )
+ {
++newp;
+ }
/*fprintf( fout, "%d, oldp: %d - %d, newp: %d - %d",
atom, oldp, old_bonds->select.bond_list[oldp].nbr,
@@ -165,8 +169,7 @@ static int Compare_Bond_Lists( int atom, control_params *control, reax_list **li
static void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
- control_params *control, reax_list *bonds, int print,
- FILE *fout )
+ control_params *control, reax_list *bonds, int print, FILE *fout )
{
int i, start, end;
@@ -174,17 +177,23 @@ static void Get_Molecule( int atom, molecule *m, int *mark, reax_system *system,
end = End_Index( atom, bonds );
if ( print )
+ {
fprintf( fout, "%5d(%2s)",
atom + 1, system->reaxprm.sbp[ system->atoms[atom].type ].name );
+ }
mark[atom] = 1;
m->atom_list[ m->atom_count++ ] = atom;
m->mtypes[ system->atoms[ atom ].type ]++;
for ( i = start; i < end; ++i )
+ {
if ( bonds->select.bond_list[i].bo_data.BO >= control->bg_cut &&
!mark[bonds->select.bond_list[i].nbr] )
+ {
Get_Molecule( bonds->select.bond_list[i].nbr, m, mark,
system, control, bonds, print, fout );
+ }
+ }
}
diff --git a/sPuReMD/src/bond_orders.c b/sPuReMD/src/bond_orders.c
index 98d26e04..0b02a6dd 100644
--- a/sPuReMD/src/bond_orders.c
+++ b/sPuReMD/src/bond_orders.c
@@ -669,36 +669,6 @@ void Add_dBond_to_Forces( int i, int pj, reax_system *system,
}
-/* Locate j on i's list.
- This function assumes that j is there for sure!
- And this is the case given our method of neighbor generation*/
-static int Locate_Symmetric_Bond( reax_list *bonds, int i, int j )
-{
- int start = Start_Index(i, bonds);
- int end = End_Index(i, bonds);
- int mid = (start + end) / 2;
- int mid_nbr;
-
- while ( (mid_nbr = bonds->select.bond_list[mid].nbr) != j )
- {
- /*fprintf( stderr, "\tstart: %d end: %d mid: %d\n",
- start, end, mid );*/
- if ( mid_nbr < j )
- {
- start = mid + 1;
- }
- else
- {
- end = mid - 1;
- }
-
- mid = (start + end) / 2;
- }
-
- return mid;
-}
-
-
static inline void Copy_Neighbor_Data( bond_data *dest, near_neighbor_data *src )
{
dest->nbr = src->nbr;
@@ -723,12 +693,6 @@ static inline void Copy_Bond_Order_Data( bond_order_data *dest, bond_order_data
}
-static int compare_bonds( const void *p1, const void *p2 )
-{
- return ((bond_data *)p1)->nbr - ((bond_data *)p2)->nbr;
-}
-
-
/* A very important and crucial assumption here is that each segment
* belonging to a different atom in nbrhoods->nbr_list is sorted in its own.
* This can either be done in the general coordinator function or here */
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index fba8bd32..39e15d77 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -398,9 +398,9 @@ int Are_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
}
-/* Determines if the distance between x1 and x2 is < vlist_cut.
- If so, this neighborhood is added to the list of far neighbors.
- Periodic boundary conditions do not apply. */
+/* Determines if the distance between atoms x1 and x2 is strictly less than
+ * vlist_cut. If so, this neighborhood is added to the list of far neighbors.
+ * Note: Periodic boundary conditions do not apply. */
void Get_NonPeriodic_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
control_params *control, far_neighbor_data *new_nbrs, int *count )
{
@@ -424,10 +424,10 @@ void Get_NonPeriodic_Far_Neighbors( rvec x1, rvec x2, simulation_box *box,
}
-/* Finds periodic neighbors in a 'big_box'. Here 'big_box' means:
- the current simulation box has all dimensions > 2 *vlist_cut.
- If the periodic distance between x1 and x2 is than vlist_cut, this
- neighborhood is added to the list of far neighbors. */
+/* Finds periodic neighbors in a 'big_box'. Here 'big_box' means the current
+ * simulation box has all dimensions strictly greater than twice of vlist_cut.
+ * If the periodic distance between x1 and x2 is than vlist_cut, this
+ * neighborhood is added to the list of far neighbors. */
void Get_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, simulation_box *box,
control_params *control, far_neighbor_data *periodic_nbrs, int *count )
{
@@ -439,9 +439,7 @@ void Get_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, simulation_box *box,
for ( i = 0; i < 3; i++ )
{
d = x2[i] - x1[i];
- tmp = SQR(d);
- // fprintf(out,"Inside Sq_Distance_on_T3, %d, %lf, %lf\n",
- // i,tmp,SQR(box->box_norms[i]/2.0));
+ tmp = SQR( d );
if ( tmp >= SQR( box->box_norms[i] / 2.0 ) )
{
@@ -482,13 +480,14 @@ void Get_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, simulation_box *box,
}
-/* Finds all periodic far neighborhoods between x1 and x2
- ((dist(x1, x2') < vlist_cut, periodic images of x2 are also considered).
- Here the box is 'small' meaning that at least one dimension is < 2*vlist_cut.
- IMPORTANT: This part might need some improvement. In NPT, the simulation box
- might get too small (such as <5 A!). In this case we have to consider the
- periodic images of x2 that are two boxs away!!!
-*/
+/* Finds all periodic far neighborhoods between atoms x1 and x2
+ * ((dist(x1, x2') < 2 * vlist_cut, periodic images of x2 are also considered).
+ * Here the box is 'small' meaning that at least one dimension is strictly
+ * less than twice of vlist_cut.
+ *
+ * NOTE: This part might need some improvement. In NPT, the simulation box
+ * might get too small (such as <5 A!). In this case we have to consider the
+ * periodic images of x2 that are two boxs away!!! */
void Get_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, simulation_box *box,
control_params *control, far_neighbor_data *periodic_nbrs, int *count )
{
@@ -498,29 +497,36 @@ void Get_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, simulation_box *box
*count = 0;
/* determine the max stretch of imaginary boxs in each direction
- to handle periodic boundary conditions correctly. */
+ * to handle periodic boundary conditions correctly */
imax = (int)(control->vlist_cut / box->box_norms[0] + 1);
jmax = (int)(control->vlist_cut / box->box_norms[1] + 1);
kmax = (int)(control->vlist_cut / box->box_norms[2] + 1);
+
/*if( imax > 1 || jmax > 1 || kmax > 1 )
fprintf( stderr, "box %8.3f x %8.3f x %8.3f --> %2d %2d %2d\n",
box->box_norms[0], box->box_norms[1], box->box_norms[2],
imax, jmax, kmax ); */
-
for ( i = -imax; i <= imax; ++i )
{
- if (FABS(d_i = ((x2[0] + i * box->box_norms[0]) - x1[0])) <= control->vlist_cut)
+ d_i = (x2[0] + i * box->box_norms[0]) - x1[0];
+
+ if ( FABS(d_i) <= control->vlist_cut )
{
for ( j = -jmax; j <= jmax; ++j )
{
- if (FABS(d_j = ((x2[1] + j * box->box_norms[1]) - x1[1])) <= control->vlist_cut)
+ d_j = (x2[1] + j * box->box_norms[1]) - x1[1];
+
+ if ( FABS(d_j) <= control->vlist_cut )
{
for ( k = -kmax; k <= kmax; ++k )
{
- if (FABS(d_k = ((x2[2] + k * box->box_norms[2]) - x1[2])) <= control->vlist_cut)
+ d_k = (x2[2] + k * box->box_norms[2]) - x1[2];
+
+ if ( FABS(d_k) <= control->vlist_cut )
{
sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
+
if ( sqr_norm <= SQR(control->vlist_cut) )
{
periodic_nbrs[ *count ].d = SQRT( sqr_norm );
@@ -553,6 +559,7 @@ void Get_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, simulation_box *box
* periodic_nbrs[ *count ].imaginary = 0;
* else periodic_nbrs[ *count ].imaginary = 1;
*/
+
++(*count);
}
}
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index 267c3542..6af97014 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -25,32 +25,6 @@
#include "two_body_interactions.h"
-void Make_Lookup_Table( real xmin, real xmax, int n,
- lookup_function f, lookup_table* t )
-{
- int i;
-
- t->xmin = xmin;
- t->xmax = xmax;
- t->n = n;
- t->dx = (xmax - xmin) / (n - 1);
- t->inv_dx = 1.0 / t->dx;
- t->a = (n - 1) / (xmax - xmin);
- t->y = (real*) smalloc( n * sizeof(real),
- "Make_Lookup_Table::t->y" );
-
- for (i = 0; i < n; i++)
- {
- t->y[i] = f(i * t->dx + t->xmin);
- }
-
- // fprintf(stdout,"dx = %lf\n",t->dx);
- // for(i=0; i < n; i++)
- // fprintf( stdout,"%d %lf %lf %lf\n",
- // i, i/t->a+t->xmin, t->y[i], exp(i/t->a+t->xmin) );
-}
-
-
/* Fills solution into x. Warning: will modify c and d! */
static void Tridiagonal_Solve( const real *a, const real *b,
real *c, real *d, real *x, unsigned int n)
@@ -68,7 +42,7 @@ static void Tridiagonal_Solve( const real *a, const real *b,
d[i] = (d[i] - d[i - 1] * a[i]) / id;
}
- /* Now back substitute. */
+ /* solve via back substitution */
x[n - 1] = d[n - 1];
for (i = n - 2; i >= 0; i--)
{
@@ -83,19 +57,19 @@ static void Natural_Cubic_Spline( const real *h, const real *f,
int i;
real *a, *b, *c, *d, *v;
- /* allocate space for the linear system */
- a = (real*) smalloc( n * sizeof(real),
+ /* allocate space for linear system */
+ a = smalloc( n * sizeof(real),
"Natural_Cubic_Spline::a" );
- b = (real*) smalloc( n * sizeof(real),
+ b = smalloc( n * sizeof(real),
"Natural_Cubic_Spline::b" );
- c = (real*) smalloc( n * sizeof(real),
+ c = smalloc( n * sizeof(real),
"Natural_Cubic_Spline::c" );
- d = (real*) smalloc( n * sizeof(real),
+ d = smalloc( n * sizeof(real),
"Natural_Cubic_Spline::d" );
- v = (real*) smalloc( n * sizeof(real),
+ v = smalloc( n * sizeof(real),
"Natural_Cubic_Spline::v" );
- /* build the linear system */
+ /* build linear system */
a[0] = 0.0;
a[1] = 0.0;
a[n - 1] = 0.0;
@@ -130,6 +104,7 @@ static void Natural_Cubic_Spline( const real *h, const real *f,
/*fprintf( stderr, "i a b c d\n" );
for( i = 0; i < n; ++i )
fprintf( stderr, "%d %f %f %f %f\n", i, a[i], b[i], c[i], d[i] );*/
+
v[0] = 0.0;
v[n - 1] = 0.0;
Tridiagonal_Solve( a + 1, b + 1, c + 1, d + 1, v + 1, n - 2 );
@@ -158,15 +133,15 @@ static void Complete_Cubic_Spline( const real *h, const real *f, real v0, real v
real *a, *b, *c, *d, *v;
/* allocate space for the linear system */
- a = (real*) smalloc( n * sizeof(real),
+ a = smalloc( n * sizeof(real),
"Complete_Cubic_Spline::a" );
- b = (real*) smalloc( n * sizeof(real),
+ b = smalloc( n * sizeof(real),
"Complete_Cubic_Spline::b" );
- c = (real*) smalloc( n * sizeof(real),
+ c = smalloc( n * sizeof(real),
"Complete_Cubic_Spline::c" );
- d = (real*) smalloc( n * sizeof(real),
+ d = smalloc( n * sizeof(real),
"Complete_Cubic_Spline::d" );
- v = (real*) smalloc( n * sizeof(real),
+ v = smalloc( n * sizeof(real),
"Complete_Cubic_Spline::v" );
/* build the linear system */
@@ -271,17 +246,17 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control,
num_atom_types = system->reaxprm.num_atom_types;
dr = control->r_cut / control->tabulate;
- h = (real*) scalloc( (control->tabulate + 1), sizeof(real),
+ h = scalloc( (control->tabulate + 1), sizeof(real),
"Make_LR_Lookup_Table::h" );
- fh = (real*) scalloc( (control->tabulate + 1), sizeof(real),
+ fh = scalloc( (control->tabulate + 1), sizeof(real),
"Make_LR_Lookup_Table::fh" );
- fvdw = (real*) scalloc( (control->tabulate + 1), sizeof(real),
+ fvdw = scalloc( (control->tabulate + 1), sizeof(real),
"Make_LR_Lookup_Table::fvdw" );
- fCEvd = (real*) scalloc( (control->tabulate + 1), sizeof(real),
+ fCEvd = scalloc( (control->tabulate + 1), sizeof(real),
"Make_LR_Lookup_Table::fCEvd" );
- fele = (real*) scalloc( (control->tabulate + 1), sizeof(real),
+ fele = scalloc( (control->tabulate + 1), sizeof(real),
"Make_LR_Lookup_Table::fele" );
- fCEclmb = (real*) scalloc( (control->tabulate + 1), sizeof(real),
+ fCEclmb = scalloc( (control->tabulate + 1), sizeof(real),
"Make_LR_Lookup_Table::fCEclmb" );
/* allocate Long-Range LookUp Table space based on
@@ -321,22 +296,22 @@ void Make_LR_Lookup_Table( reax_system *system, control_params *control,
workspace->LR[i][j].n = control->tabulate + 1;
workspace->LR[i][j].dx = dr;
workspace->LR[i][j].inv_dx = control->tabulate / control->r_cut;
- workspace->LR[i][j].y = (LR_data*)
+ workspace->LR[i][j].y =
smalloc( workspace->LR[i][j].n * sizeof(LR_data),
"Make_LR_Lookup_Table::LR[i][j].y" );
- workspace->LR[i][j].H = (cubic_spline_coef*)
+ workspace->LR[i][j].H =
smalloc( workspace->LR[i][j].n * sizeof(cubic_spline_coef),
"Make_LR_Lookup_Table::LR[i][j].H" );
- workspace->LR[i][j].vdW = (cubic_spline_coef*)
+ workspace->LR[i][j].vdW =
smalloc( workspace->LR[i][j].n * sizeof(cubic_spline_coef),
"Make_LR_Lookup_Table::LR[i][j].vdW" );
- workspace->LR[i][j].CEvd = (cubic_spline_coef*)
+ workspace->LR[i][j].CEvd =
smalloc( workspace->LR[i][j].n * sizeof(cubic_spline_coef),
"Make_LR_Lookup_Table::LR[i][j].CEvd" );
- workspace->LR[i][j].ele = (cubic_spline_coef*)
+ workspace->LR[i][j].ele =
smalloc( workspace->LR[i][j].n * sizeof(cubic_spline_coef),
"Make_LR_Lookup_Table::LR[i][j].ele" );
- workspace->LR[i][j].CEclmb = (cubic_spline_coef*)
+ workspace->LR[i][j].CEclmb =
smalloc( workspace->LR[i][j].n * sizeof(cubic_spline_coef),
"Make_LR_Lookup_Table::LR[i][j].CEclmb" );
@@ -499,40 +474,3 @@ void Finalize_LR_Lookup_Table( reax_system *system, control_params *control,
sfree( workspace->LR, "Finalize_LR_Lookup_Table::LR" );
}
-
-
-static int Lookup_Index_Of( real x, lookup_table* t )
-{
- return (int)( t->a * ( x - t->xmin ) );
-}
-
-
-real Lookup( real x, lookup_table* t )
-{
- real x1, x2;
- real b;
- int i;
-
- /*
- if ( x < t->xmin)
- {
- fprintf(stderr,"Domain check %lf > %lf\n",t->xmin,x);
- exit(0);
- }
- if ( x > t->xmax)
- {
- fprintf(stderr,"Domain check %lf < %lf\n",t->xmax,x);
- exit(0);
- }
- */
-
- i = Lookup_Index_Of( x, t );
- x1 = i * t->dx + t->xmin;
- x2 = (i + 1) * t->dx + t->xmin;
-
- b = ( x2 * t->y[i] - x1 * t->y[i + 1] ) * t->inv_dx;
- // fprintf( stdout,"SLookup_Entry: %d, %lf, %lf, %lf, %lf: %lf, %lf\n",
- // i,x1,x2,x,b,t->one_over_dx*(t->y[i+1]-t->y[i])*x+b,exp(x));
-
- return t->inv_dx * ( t->y[i + 1] - t->y[i] ) * x + b;
-}
diff --git a/sPuReMD/src/lookup.h b/sPuReMD/src/lookup.h
index 8a63a1a9..109b54d1 100644
--- a/sPuReMD/src/lookup.h
+++ b/sPuReMD/src/lookup.h
@@ -25,13 +25,6 @@
#include "reax_types.h"
-typedef real (*lookup_function)(real);
-
-
-void Make_Lookup_Table( real, real, int, lookup_function, lookup_table* );
-
-real Lookup( real, lookup_table* );
-
void Make_LR_Lookup_Table( reax_system*, control_params*,
static_storage* );
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index 441cb238..b2937ea9 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -30,15 +30,10 @@
#include "vector.h"
-/* Function pointer definitions */
-typedef void (*get_far_neighbors_function)(rvec, rvec, simulation_box*,
- control_params*, far_neighbor_data*, int*);
-
-
static inline real DistSqr_to_CP( rvec cp, rvec x )
{
- int i;
- real d_sqr = 0;
+ int i;
+ real d_sqr = 0.0;
for ( i = 0; i < 3; ++i )
{
@@ -52,9 +47,8 @@ static inline real DistSqr_to_CP( rvec cp, rvec x )
}
-void Generate_Neighbor_Lists( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+int Estimate_NumNeighbors( reax_system *system, control_params *control,
+ static_storage *workspace, reax_list **lists )
{
int i, j, k, l, m, itr;
int x, y, z;
@@ -64,17 +58,11 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
ivec *nbrs;
rvec *nbrs_cp;
grid *g;
- reax_list *far_nbrs;
- far_neighbor_data *nbr_data;
- real t_start, t_elapsed;
+ far_neighbor_data nbr_data;
- t_start = Get_Time( );
- // fprintf( stderr, "\n\tentered nbrs - " );
- g = &( system->g );
- far_nbrs = lists[FAR_NBRS];
+ g = &system->g;
Bin_Atoms( system, workspace );
- // fprintf( stderr, "atoms sorted - " );
num_far = 0;
@@ -87,14 +75,11 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
{
nbrs = g->nbrs[i][j][k];
nbrs_cp = g->nbrs_cp[i][j][k];
- //fprintf( stderr, "gridcell %d %d %d\n", i, j, k );
/* pick up an atom from the current cell */
- for (l = 0; l < g->top[i][j][k]; ++l )
+ for ( l = 0; l < g->top[i][j][k]; ++l )
{
atom1 = g->atoms[i][j][k][l];
- Set_Start_Index( atom1, num_far, far_nbrs );
- //fprintf( stderr, "\tatom %d\n", atom1 );
itr = 0;
while ( nbrs[itr][0] >= 0 )
@@ -102,29 +87,27 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
x = nbrs[itr][0];
y = nbrs[itr][1];
z = nbrs[itr][2];
- //fprintf( stderr, "\t\tgridcell %d %d %d\n", x, y, z );
- if ( DistSqr_to_CP(nbrs_cp[itr], system->atoms[atom1].x ) <=
- SQR(control->vlist_cut) )
+ if ( DistSqr_to_CP(nbrs_cp[itr], system->atoms[atom1].x )
+ <= SQR(control->vlist_cut) )
{
nbr_atoms = g->atoms[x][y][z];
max = g->top[x][y][z];
- //fprintf( stderr, "\t\tmax: %d\n", max );
- /* pick up another atom from the neighbor cell */
+ /* pick up another atom from the neighbor cell -
+ * we have to compare atom1 with its own periodic images as well,
+ * that's why there is also equality in the if stmt below */
for ( m = 0; m < max; ++m )
{
atom2 = nbr_atoms[m];
+
+ //if( nbrs[itr+1][0] >= 0 || atom1 > atom2 ) {
if ( atom1 > atom2 )
{
- nbr_data = &(far_nbrs->select.far_nbr_list[num_far]);
- //fprintf (stderr, " %f %f %f \n", nbr_data->dvec[0], nbr_data->dvec[1], nbr_data->dvec[2]);
if ( Are_Far_Neighbors(system->atoms[atom1].x,
system->atoms[atom2].x,
- &(system->box), control->vlist_cut,
- nbr_data) )
+ &system->box, control->vlist_cut, &nbr_data) )
{
- nbr_data->nbr = atom2;
++num_far;
}
}
@@ -133,66 +116,109 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
++itr;
}
-
- Set_End_Index( atom1, 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);
}
}
}
}
- if ( num_far > far_nbrs->total_intrs * DANGER_ZONE )
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr, "estimate nbrs done, num_far: %d\n", num_far );
+#endif
+
+ return num_far * SAFE_ZONE;
+}
+
+
+/* If the code is not compiled to handle small periodic boxes (i.e. a
+ * simulation box with any dimension less than twice the Verlet list cutoff
+ * distance, vlist_cut), it will use the optimized Generate_Neighbor_Lists
+ * function. Otherwise it will execute the neighbor routine with small
+ * periodic box support.
+ *
+ * Define the preprocessor definition SMALL_BOX_SUPPORT to enable (in
+ * reax_types.h). */
+#if defined SMALL_BOX_SUPPORT
+typedef void (*get_far_neighbors_function)( rvec, rvec, simulation_box*,
+ control_params*, far_neighbor_data*, int* );
+
+
+void Choose_Neighbor_Finder( reax_system *system, control_params *control,
+ get_far_neighbors_function *Get_Far_Neighbors )
+{
+ if ( control->periodic_boundaries )
{
- workspace->realloc.num_far = num_far;
- if ( num_far > far_nbrs->total_intrs )
+ if ( system->box.box_norms[0] > 2.0 * control->vlist_cut
+ && system->box.box_norms[1] > 2.0 * control->vlist_cut
+ && system->box.box_norms[2] > 2.0 * control->vlist_cut )
{
- fprintf( stderr, "step%d-ran out of space on far_nbrs: top=%d, max=%d",
- data->step, num_far, far_nbrs->total_intrs );
- exit( INSUFFICIENT_MEMORY );
+ *Get_Far_Neighbors = &Get_Periodic_Far_Neighbors_Big_Box;
+ }
+ else
+ {
+ *Get_Far_Neighbors = &Get_Periodic_Far_Neighbors_Small_Box;
}
}
-
- t_elapsed = Get_Timing_Info( t_start );
- data->timing.nbrs += t_elapsed;
-
-#if defined(DEBUG)
- for ( i = 0; i < system->N; ++i )
+ else
{
- qsort( &(far_nbrs->select.far_nbr_list[ Start_Index(i, far_nbrs) ]),
- Num_Entries(i, far_nbrs), sizeof(far_neighbor_data),
- compare_far_nbrs );
+ *Get_Far_Neighbors = &Get_NonPeriodic_Far_Neighbors;
}
-#endif
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "nbrs - ");
- fprintf( stderr, "nbrs done, num_far: %d\n", num_far );
-#endif
-#if defined(TEST_ENERGY)
- //Print_Far_Neighbors( system, control, workspace, lists );
-#endif
}
-int Estimate_NumNeighbors( reax_system *system, control_params *control,
- static_storage *workspace, reax_list **lists )
+int compare_far_nbrs(const void *v1, const void *v2)
{
- int i, j, k, l, m, itr;
- int x, y, z;
- int atom1, atom2, max;
- int num_far;
- int *nbr_atoms;
+ return ((*(far_neighbor_data *)v1).nbr - (*(far_neighbor_data *)v2).nbr);
+}
+
+
+static inline void Set_Far_Neighbor( far_neighbor_data *dest, int nbr, real d, real C,
+ rvec dvec, ivec rel_box/*, rvec ext_factor*/ )
+{
+ dest->nbr = nbr;
+ dest->d = d;
+ rvec_Scale( dest->dvec, C, dvec );
+ ivec_Copy( dest->rel_box, rel_box );
+ // rvec_Scale( dest->ext_factor, C, ext_factor );
+}
+
+
+void Generate_Neighbor_Lists( reax_system *system, control_params *control,
+ simulation_data *data, static_storage *workspace,
+ reax_list **lists, output_controls *out_control )
+{
+ int i, j, k, l, m, itr;
+ int x, y, z;
+ int atom1, atom2, max;
+ int num_far, c, count;
+ int *nbr_atoms;
ivec *nbrs;
rvec *nbrs_cp;
grid *g;
- far_neighbor_data nbr_data;
+ reax_list *far_nbrs;
+ get_far_neighbors_function Get_Far_Neighbors;
+ far_neighbor_data new_nbrs[125];
+ real t_start, t_elapsed;
+
+ t_start = Get_Time( );
+ g = &system->g;
+ far_nbrs = lists[FAR_NBRS];
+
+ if ( control->ensemble == iNPT || control->ensemble == sNPT
+ || control->ensemble == NPT )
+ {
+ Update_Grid( system );
+ }
+
+ Bin_Atoms( system, out_control );
+
+#ifdef REORDER_ATOMS
+ Cluster_Atoms( system, workspace, control );
+#endif
+
+ Choose_Neighbor_Finder( system, control, &Get_Far_Neighbors );
- // fprintf( stderr, "\n\tentered nbrs - " );
- g = &( system->g );
- Bin_Atoms( system, workspace );
- // fprintf( stderr, "atoms sorted - " );
num_far = 0;
+ c = 0;
/* first pick up a cell in the grid */
for ( i = 0; i < g->ncell[0]; i++ )
@@ -203,152 +229,206 @@ int Estimate_NumNeighbors( reax_system *system, control_params *control,
{
nbrs = g->nbrs[i][j][k];
nbrs_cp = g->nbrs_cp[i][j][k];
- //fprintf( stderr, "gridcell %d %d %d\n", i, j, k );
/* pick up an atom from the current cell */
- for (l = 0; l < g->top[i][j][k]; ++l )
+ for ( l = 0; l < g->top[i][j][k]; ++l )
{
atom1 = g->atoms[i][j][k][l];
+ Set_Start_Index( atom1, num_far, far_nbrs );
itr = 0;
- while ( nbrs[itr][0] >= 0 )
+ while ( nbrs[itr][0] > 0 )
{
x = nbrs[itr][0];
y = nbrs[itr][1];
z = nbrs[itr][2];
- //fprintf( stderr, "\t\tgridcell %d %d %d\n", x, y, z );
- if ( DistSqr_to_CP(nbrs_cp[itr], system->atoms[atom1].x ) <=
- SQR(control->vlist_cut) )
+ // if( DistSqr_to_CP(nbrs_cp[itr], system->atoms[atom1].x ) <=
+ // SQR(control->r_cut))
+ nbr_atoms = g->atoms[x][y][z];
+ max = g->top[x][y][z];
+
+ /* pick up another atom from the neighbor cell -
+ * we have to compare atom1 with its own periodic images as well,
+ * that's why there is also equality in the if stmt below */
+ for ( m = 0; m < max; ++m )
{
- nbr_atoms = g->atoms[x][y][z];
- max = g->top[x][y][z];
- //fprintf( stderr, "\t\tmax: %d\n", max );
+ atom2 = nbr_atoms[m];
- /* pick up another atom from the neighbor cell -
- we have to compare atom1 with its own periodic images as well,
- that's why there is also equality in the if stmt below */
- for ( m = 0; m < max; ++m )
+ if ( atom1 >= atom2 )
{
- atom2 = nbr_atoms[m];
- //if( nbrs[itr+1][0] >= 0 || atom1 > atom2 ) {
- if ( atom1 > atom2 )
+ Get_Far_Neighbors( system->atoms[atom1].x, system->atoms[atom2].x,
+ &system->box, control, new_nbrs, &count );
+
+ for ( c = 0; c < count; ++c )
{
- if ( Are_Far_Neighbors(system->atoms[atom1].x,
- system->atoms[atom2].x,
- &(system->box), control->vlist_cut,
- &nbr_data) )
+ if ( atom1 != atom2 || (atom1 == atom2 && new_nbrs[c].d >= 0.1) )
+ {
+ Set_Far_Neighbor( &far_nbrs->select.far_nbr_list[num_far],
+ atom2, new_nbrs[c].d, 1.0,
+ new_nbrs[c].dvec, new_nbrs[c].rel_box );
++num_far;
+ }
}
}
}
++itr;
}
+
+ Set_End_Index( atom1, num_far, far_nbrs );
}
}
}
}
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "estimate nbrs done, num_far: %d\n", num_far );
-#endif
- return num_far * SAFE_ZONE;
-}
-
-
+ far_nbrs->total_intrs = num_far;
-#if defined DONE
-void Choose_Neighbor_Finder( reax_system *system, control_params *control,
- get_far_neighbors_function *Get_Far_Neighbors )
-{
- if ( control->periodic_boundaries )
- {
- if ( system->box.box_norms[0] > 2.0 * control->vlist_cut &&
- system->box.box_norms[1] > 2.0 * control->vlist_cut &&
- system->box.box_norms[2] > 2.0 * control->vlist_cut )
- {
- (*Get_Far_Neighbors) = Get_Periodic_Far_Neighbors_Big_Box;
- }
- else
- {
- (*Get_Far_Neighbors) = Get_Periodic_Far_Neighbors_Small_Box;
- }
- }
- else
+#if defined(DEBUG)
+ for ( i = 0; i < system->N; ++i )
{
- (*Get_Far_Neighbors) = Get_NonPeriodic_Far_Neighbors;
+ qsort( &far_nbrs->select.far_nbr_list[ Start_Index(i, far_nbrs) ],
+ Num_Entries(i, far_nbrs), sizeof(far_neighbor_data),
+ compare_far_nbrs );
}
-}
+ fprintf( stderr, "step%d: num of farnbrs=%6d\n", data->step, num_far );
+ fprintf( stderr, "\tallocated farnbrs: %6d\n",
+ system->N * far_nbrs->intrs_per_unit );
+#endif
-int compare_near_nbrs(const void *v1, const void *v2)
-{
- return ((*(near_neighbor_data *)v1).nbr - (*(near_neighbor_data *)v2).nbr);
+ t_elapsed = Get_Timing_Info( t_start );
+ data->timing.nbrs += t_elapsed;
}
-int compare_far_nbrs(const void *v1, const void *v2)
+#else
+void Generate_Neighbor_Lists( reax_system *system, control_params *control,
+ simulation_data *data, static_storage *workspace,
+ reax_list **lists, output_controls *out_control )
{
- return ((*(far_neighbor_data *)v1).nbr - (*(far_neighbor_data *)v2).nbr);
-}
+ int i, j, k, l, m, itr;
+ int x, y, z;
+ int atom1, atom2, max;
+ int num_far;
+ int *nbr_atoms;
+ ivec *nbrs;
+ rvec *nbrs_cp;
+ grid *g;
+ reax_list *far_nbrs;
+ far_neighbor_data *nbr_data;
+ real t_start, t_elapsed;
+ t_start = Get_Time( );
+ g = &system->g;
+ far_nbrs = lists[FAR_NBRS];
-static inline void Set_Far_Neighbor( far_neighbor_data *dest, int nbr, real d, real C,
- rvec dvec, ivec rel_box/*, rvec ext_factor*/ )
-{
- dest->nbr = nbr;
- dest->d = d;
- rvec_Scale( dest->dvec, C, dvec );
- ivec_Copy( dest->rel_box, rel_box );
- // rvec_Scale( dest->ext_factor, C, ext_factor );
-}
+ Bin_Atoms( system, workspace );
+#ifdef REORDER_ATOMS
+ Cluster_Atoms( system, workspace, control );
+#endif
-static inline void Set_Near_Neighbor(near_neighbor_data *dest, int nbr, real d, real C,
- rvec dvec, ivec rel_box/*, rvec ext_factor*/)
-{
- dest->nbr = nbr;
- dest->d = d;
- rvec_Scale( dest->dvec, C, dvec );
- ivec_Scale( dest->rel_box, C, rel_box );
- // rvec_Scale( dest->ext_factor, C, ext_factor );
-}
+ num_far = 0;
-/* In case bond restrictions are applied, this method checks if
- atom1 and atom2 are allowed to bond with each other */
-static inline int can_Bond( static_storage *workspace, int atom1, int atom2 )
-{
- int i;
+ /* first pick up a cell in the grid */
+ for ( i = 0; i < g->ncell[0]; i++ )
+ {
+ for ( j = 0; j < g->ncell[1]; j++ )
+ {
+ for ( k = 0; k < g->ncell[2]; k++ )
+ {
+ nbrs = g->nbrs[i][j][k];
+ nbrs_cp = g->nbrs_cp[i][j][k];
- // fprintf( stderr, "can bond %6d %6d?\n", atom1, atom2 );
+ /* pick up an atom from the current cell */
+ for ( l = 0; l < g->top[i][j][k]; ++l )
+ {
+ atom1 = g->atoms[i][j][k][l];
+ Set_Start_Index( atom1, num_far, far_nbrs );
- if ( !workspace->restricted[ atom1 ] && !workspace->restricted[ atom2 ] )
- {
- return FALSE;
+ itr = 0;
+ while ( nbrs[itr][0] >= 0 )
+ {
+ x = nbrs[itr][0];
+ y = nbrs[itr][1];
+ z = nbrs[itr][2];
+
+ if ( DistSqr_to_CP(nbrs_cp[itr], system->atoms[atom1].x ) <=
+ SQR(control->vlist_cut) )
+ {
+ nbr_atoms = g->atoms[x][y][z];
+ max = g->top[x][y][z];
+
+ /* pick up another atom from the neighbor cell */
+ for ( m = 0; m < max; ++m )
+ {
+ atom2 = nbr_atoms[m];
+
+ if ( atom1 > atom2 )
+ {
+ nbr_data = &far_nbrs->select.far_nbr_list[num_far];
+
+ if ( Are_Far_Neighbors(system->atoms[atom1].x,
+ system->atoms[atom2].x, &system->box,
+ control->vlist_cut, nbr_data) )
+ {
+ nbr_data->nbr = atom2;
+ ++num_far;
+ }
+ }
+ }
+ }
+
+ ++itr;
+ }
+
+ Set_End_Index( atom1, 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);
+ }
+ }
+ }
}
- for ( i = 0; i < workspace->restricted[ atom1 ]; ++i )
+ if ( num_far > far_nbrs->total_intrs * DANGER_ZONE )
{
- if ( workspace->restricted_list[ atom1 ][i] == atom2 )
+ workspace->realloc.num_far = num_far;
+ if ( num_far > far_nbrs->total_intrs )
{
- return FALSE;
+ fprintf( stderr, "step%d-ran out of space on far_nbrs: top=%d, max=%d",
+ data->step, num_far, far_nbrs->total_intrs );
+ exit( INSUFFICIENT_MEMORY );
}
}
- for ( i = 0; i < workspace->restricted[ atom2 ]; ++i )
+#if defined(DEBUG)
+ for ( i = 0; i < system->N; ++i )
{
- if ( workspace->restricted_list[ atom2 ][i] == atom1 )
- {
- return FALSE;
- }
+ qsort( &far_nbrs->select.far_nbr_list[ Start_Index(i, far_nbrs) ],
+ Num_Entries(i, far_nbrs), sizeof(far_neighbor_data),
+ compare_far_nbrs );
}
+#endif
- return TRUE;
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr, "nbrs - ");
+ fprintf( stderr, "nbrs done, num_far: %d\n", num_far );
+#endif
+
+#if defined(TEST_ENERGY)
+ //Print_Far_Neighbors( system, control, workspace, lists );
+#endif
+
+ t_elapsed = Get_Timing_Info( t_start );
+ data->timing.nbrs += t_elapsed;
}
+#endif
+#if defined(LEGACY)
/* check if atom2 is on atom1's near neighbor list */
static inline int is_Near_Neighbor( reax_list *near_nbrs, int atom1, int atom2 )
{
@@ -358,7 +438,6 @@ static inline int is_Near_Neighbor( reax_list *near_nbrs, int atom1, int atom2 )
{
if ( near_nbrs->select.near_nbr_list[i].nbr == atom2 )
{
- // fprintf( stderr, "near neighbors %6d %6d\n", atom1, atom2 );
return FALSE;
}
}
@@ -367,9 +446,26 @@ static inline int is_Near_Neighbor( reax_list *near_nbrs, int atom1, int atom2 )
}
+int compare_near_nbrs(const void *v1, const void *v2)
+{
+ return ((*(near_neighbor_data *)v1).nbr - (*(near_neighbor_data *)v2).nbr);
+}
+
+
+static inline void Set_Near_Neighbor( near_neighbor_data *dest, int nbr, real d, real C,
+ rvec dvec, ivec rel_box/*, rvec ext_factor*/ )
+{
+ dest->nbr = nbr;
+ dest->d = d;
+ rvec_Scale( dest->dvec, C, dvec );
+ ivec_Scale( dest->rel_box, C, rel_box );
+ // rvec_Scale( dest->ext_factor, C, ext_factor );
+}
+
+
void Generate_Neighbor_Lists( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ reax_list **lists, output_controls *out_control )
{
int i, j, k;
int x, y, z;
@@ -380,36 +476,27 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
int grid_top;
grid *g;
reax_list *far_nbrs;
- //int hb_type1, hb_type2;
- //reax_list *hbonds = &(*lists)[HBOND];
- //int top_hbond1, top_hbond2;
get_far_neighbors_function Get_Far_Neighbors;
far_neighbor_data new_nbrs[125];
+ real t_start, t_elapsed;
- g = &( system->g );
+ t_start = Get_Time( );
+ g = &system->g;
far_nbrs = lists[FAR_NBRS];
- // fprintf( stderr, "\n\tentered nbrs - " );
- if ( control->ensemble == iNPT || control->ensemble == sNPT ||
- control->ensemble == NPT )
+ if ( control->ensemble == iNPT || control->ensemble == sNPT
+ || control->ensemble == NPT )
{
Update_Grid( system );
}
- // fprintf( stderr, "grid updated - " );
Bin_Atoms( system, out_control );
- // fprintf( stderr, "atoms sorted - " );
#ifdef REORDER_ATOMS
Cluster_Atoms( system, workspace, control );
- // fprintf( stderr, "atoms clustered - " );
#endif
Choose_Neighbor_Finder( system, control, &Get_Far_Neighbors );
- // fprintf( stderr, "function chosen - " );
-
- Reset_Neighbor_Lists( system, workspace, lists );
- // fprintf( stderr, "lists cleared - " );
num_far = 0;
num_near = 0;
@@ -426,14 +513,10 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
nbrs_cp = g->nbrs_cp[i][j][k];
/* pick up an atom from the current cell */
- //#ifdef REORDER_ATOMS
- // for(atom1 = g->start[i][j][k]; atom1 < g->end[i][j][k]; atom1++)
- //#else
- for (l = 0; l < g->top[i][j][k]; ++l )
+ for ( l = 0; l < g->top[i][j][k]; ++l )
{
atom1 = g->atoms[i][j][k][l];
Set_End_Index( atom1, num_far, far_nbrs );
- // fprintf( stderr, "atom %d:\n", atom1 );
itr = 0;
while ( nbrs[itr][0] > 0 )
@@ -445,63 +528,30 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
// if( DistSqr_to_CP(nbrs_cp[itr], system->atoms[atom1].x ) <=
// SQR(control->r_cut))
nbr_atoms = g->atoms[x][y][z];
- max_atoms = g->top[x][y][z];
+ max = g->top[x][y][z];
/* pick up another atom from the neighbor cell -
- we have to compare atom1 with its own periodic images as well,
- that's why there is also equality in the if stmt below */
- //#ifdef REORDER_ATOMS
- //for(atom2=g->start[x][y][z]; atom2<g->end[x][y][z]; atom2++)
- //#else
- for ( m = 0, atom2 = nbr_atoms[m]; m < max; ++m, atom2 = nbr_atoms[m] )
+ * we have to compare atom1 with its own periodic images as well,
+ * that's why there is also equality in the if stmt below */
+ for ( m = 0; m < max; ++m )
{
+ atom2 = nbr_atoms[m];
+
if ( atom1 >= atom2 )
{
- //fprintf( stderr, "\tatom2 %d", atom2 );
//top_near1 = End_Index( atom1, near_nbrs );
//Set_Start_Index( atom1, num_far, far_nbrs );
- //hb_type1=system->reaxprm.sbp[system->atoms[atom1].type].p_hbond;
- Get_Far_Neighbors( system->atoms[atom1].x,
- system->atoms[atom2].x,
- &(system->box), control, new_nbrs, &count );
- fprintf( stderr, "\t%d count:%d\n", atom2, count );
+ Get_Far_Neighbors( system->atoms[atom1].x, system->atoms[atom2].x,
+ &system->box, control, new_nbrs, &count );
for ( c = 0; c < count; ++c )
{
- if (atom1 != atom2 || (atom1 == atom2 && new_nbrs[c].d >= 0.1))
+ if ( atom1 != atom2 || (atom1 == atom2 && new_nbrs[c].d >= 0.1) )
{
- Set_Far_Neighbor(&(far_nbrs->select.far_nbr_list[num_far]),
- atom2, new_nbrs[c].d, 1.0,
- new_nbrs[c].dvec, new_nbrs[c].rel_box );
+ Set_Far_Neighbor( &far_nbrs->select.far_nbr_list[num_far],
+ atom2, new_nbrs[c].d, 1.0,
+ new_nbrs[c].dvec, new_nbrs[c].rel_box );
++num_far;
-
- /*fprintf(stderr,"FARNBR:%6d%6d%8.3f[%8.3f%8.3f%8.3f]\n",
- atom1, atom2, new_nbrs[c].d,
- new_nbrs[c].dvec[0], new_nbrs[c].dvec[1],
- new_nbrs[c].dvec[2] ); */
-
-
- /* hydrogen bond lists */
- /*if( control->hb_cut > 0.1 &&
- new_nbrs[c].d <= control->hb_cut ) {
- // fprintf( stderr, "%d %d\n", atom1, atom2 );
- hb_type2=system->reaxprm.sbp[system->atoms[atom2].type].p_hbond;
- if( hb_type1 == 1 && hb_type2 == 2 ) {
- top_hbond1=End_Index(workspace->hbond_index[atom1],hbonds);
- Set_Near_Neighbor(&(hbonds->select.hbond_list[top_hbond1]),
- atom2, new_nbrs[c].d, 1.0, new_nbrs[c].dvec,
- new_nbrs[c].rel_box );
- Set_End_Index( workspace->hbond_index[atom1],
- top_hbond1 + 1, hbonds );
- }
- else if( hb_type1 == 2 && hb_type2 == 1 ) {
- top_hbond2 = End_Index( workspace->hbond_index[atom2], hbonds );
- Set_Near_Neighbor(&(hbonds->select.hbond_list[top_hbond2]),
- atom1, new_nbrs[c].d, -1.0, new_nbrs[c].dvec,
- new_nbrs[c].rel_box );
- Set_End_Index( workspace->hbond_index[atom2],
- top_hbond2 + 1, hbonds );
- }*/
}
}
}
@@ -514,9 +564,15 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
}
}
- fprintf( stderr, "nbrs done-" );
-
-
+ /* Below is some piece of legacy code. It tries to force bonded
+ * interactions between atoms which are given to be bonded in the
+ * BGF file through the CONECT lines. The aim is to ensure that
+ * molecules stay intact during an energy minimization procudure.
+ * Energy minimization is not actually not implemented as of now
+ * (July 7, 2014), but we mimic it by running at very low temperature
+ * and small dt (timestep length) NVT simulation.
+ * However, I (HMA) am not sure if it can possibly achieve the desired
+ * effect. Therefore, this functionality is currently disabled. */
/* apply restrictions on near neighbors only */
if ( (data->step - data->prev_steps) < control->restrict_bonds )
{
@@ -524,35 +580,31 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
{
if ( workspace->restricted[ atom1 ] )
{
- // fprintf( stderr, "atom1: %d\n", atom1 );
-
top_near1 = End_Index( atom1, near_nbrs );
for ( j = 0; j < workspace->restricted[ atom1 ]; ++j )
{
- if (is_Near_Neighbor(near_nbrs, atom1,
- atom2 = workspace->restricted_list[atom1][j]) == FALSE)
+ atom2 = workspace->restricted_list[atom1][j];
+
+ if ( is_Near_Neighbor(near_nbrs, atom1, atom2) == FALSE )
{
fprintf( stderr, "%3d-%3d: added bond by applying restrictions!\n",
- atom1, atom2 );
+ atom1, atom2 );
top_near2 = End_Index( atom2, near_nbrs );
/* we just would like to get the nearest image, so a call to
- Get_Periodic_Far_Neighbors_Big_Box is good enough. */
+ * Get_Periodic_Far_Neighbors_Big_Box is good enough. */
Get_Periodic_Far_Neighbors_Big_Box( system->atoms[ atom1 ].x,
- system->atoms[ atom2 ].x,
- &(system->box), control,
- new_nbrs, &count );
+ system->atoms[ atom2 ].x, &system->box, control,
+ new_nbrs, &count );
- Set_Near_Neighbor( &(near_nbrs->select.near_nbr_list[ top_near1 ]),
- atom2, new_nbrs[c].d, 1.0,
- new_nbrs[c].dvec, new_nbrs[c].rel_box );
+ Set_Near_Neighbor( &near_nbrs->select.near_nbr_list[ top_near1 ], atom2,
+ new_nbrs[c].d, 1.0, new_nbrs[c].dvec, new_nbrs[c].rel_box );
++top_near1;
- Set_Near_Neighbor( &(near_nbrs->select.near_nbr_list[ top_near2 ]),
- atom1, new_nbrs[c].d, -1.0,
- new_nbrs[c].dvec, new_nbrs[c].rel_box );
+ Set_Near_Neighbor( &near_nbrs->select.near_nbr_list[ top_near2 ], atom1,
+ new_nbrs[c].d, -1.0, new_nbrs[c].dvec, new_nbrs[c].rel_box );
Set_End_Index( atom2, top_near2 + 1, near_nbrs );
}
}
@@ -561,8 +613,6 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
}
}
}
- // fprintf( stderr, "restrictions applied-" );
-
/* verify nbrlists, count total_intrs, sort nearnbrs */
near_nbrs->total_intrs = 0;
@@ -577,7 +627,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
exit( RUNTIME_ERROR );
}
- near_nbrs->total_intrs += Num_Entries(i, near_nbrs);
+ near_nbrs->total_intrs += Num_Entries( i, near_nbrs );
if ( End_Index(i, far_nbrs) > Start_Index(i + 1, far_nbrs) )
{
@@ -587,17 +637,15 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
exit( RUNTIME_ERROR );
}
- far_nbrs->total_intrs += Num_Entries(i, far_nbrs);
+ far_nbrs->total_intrs += Num_Entries( i, far_nbrs );
}
for ( i = 0; i < system->N; ++i )
{
- qsort( &(near_nbrs->select.near_nbr_list[ Start_Index(i, near_nbrs) ]),
+ qsort( &near_nbrs->select.near_nbr_list[ Start_Index(i, near_nbrs) ],
Num_Entries(i, near_nbrs), sizeof(near_neighbor_data),
compare_near_nbrs );
}
- // fprintf( stderr, "near nbrs sorted\n" );
-
#ifdef TEST_ENERGY
/* for( i = 0; i < system->N; ++i ) {
@@ -610,137 +658,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
num_near / system->N );
#endif
- //fprintf( stderr, "step%d: num of nearnbrs = %6d num of farnbrs: %6d\n",
- // data->step, num_near, num_far );
-
- //fprintf( stderr, "\talloc nearnbrs = %6d alloc farnbrs: %6d\n",
- // system->N * near_nbrs->intrs_per_unit,
- // system->N * far_nbrs->intrs_per_unit );
-}
-
-
-void Generate_Neighbor_Lists( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
-{
- int i, j, k, l, m, itr;
- int x, y, z;
- int atom1, atom2, max;
- int num_far, c, count;
- int *nbr_atoms;
- ivec *nbrs;
- rvec *nbrs_cp;
- grid *g;
- reax_list *far_nbrs;
- get_far_neighbors_function Get_Far_Neighbors;
- far_neighbor_data new_nbrs[125];
-
- g = &( system->g );
- far_nbrs = lists[FAR_NBRS];
-
- // fprintf( stderr, "\n\tentered nbrs - " );
- if ( control->ensemble == iNPT ||
- control->ensemble == sNPT ||
- control->ensemble == NPT )
- {
- Update_Grid( system );
- }
- // fprintf( stderr, "grid updated - " );
-
- Bin_Atoms( system, out_control );
- // fprintf( stderr, "atoms sorted - " );
- Choose_Neighbor_Finder( system, control, &Get_Far_Neighbors );
- // fprintf( stderr, "function chosen - " );
- Reset_Neighbor_Lists( system, workspace, lists );
- // fprintf( stderr, "lists cleared - " );
-
- num_far = 0;
- c = 0;
-
- /* first pick up a cell in the grid */
- for ( i = 0; i < g->ncell[0]; i++ )
- {
- for ( j = 0; j < g->ncell[1]; j++ )
- {
- for ( k = 0; k < g->ncell[2]; k++ )
- {
- nbrs = g->nbrs[i][j][k];
- nbrs_cp = g->nbrs_cp[i][j][k];
- fprintf( stderr, "gridcell %d %d %d\n", i, j, k );
-
- /* pick up an atom from the current cell */
- for (l = 0; l < g->top[i][j][k]; ++l )
- {
- atom1 = g->atoms[i][j][k][l];
- Set_Start_Index( atom1, num_far, far_nbrs );
- fprintf( stderr, "\tatom %d\n", atom1 );
-
- itr = 0;
- while ( nbrs[itr][0] > 0 )
- {
- x = nbrs[itr][0];
- y = nbrs[itr][1];
- z = nbrs[itr][2];
- fprintf( stderr, "\t\tgridcell %d %d %d\n", x, y, z );
-
- // if( DistSqr_to_CP(nbrs_cp[itr], system->atoms[atom1].x ) <=
- // SQR(control->r_cut))
- nbr_atoms = g->atoms[x][y][z];
- max = g->top[x][y][z];
- fprintf( stderr, "\t\tmax: %d\n", max );
-
-
- /* pick up another atom from the neighbor cell -
- we have to compare atom1 with its own periodic images as well,
- that's why there is also equality in the if stmt below */
- for ( m = 0, atom2 = nbr_atoms[m]; m < max; ++m, atom2 = nbr_atoms[m] )
- {
- if ( atom1 >= atom2 )
- {
- Get_Far_Neighbors( system->atoms[atom1].x,
- system->atoms[atom2].x,
- &(system->box), control, new_nbrs, &count );
- fprintf( stderr, "\t\t\t%d count:%d\n", atom2, count );
-
- for ( c = 0; c < count; ++c )
- if (atom1 != atom2 || (atom1 == atom2 && new_nbrs[c].d >= 0.1))
- {
- Set_Far_Neighbor(&(far_nbrs->select.far_nbr_list[num_far]),
- atom2, new_nbrs[c].d, 1.0,
- new_nbrs[c].dvec, new_nbrs[c].rel_box );
- ++num_far;
-
- /*fprintf(stderr,"FARNBR:%6d%6d%8.3f[%8.3f%8.3f%8.3f]\n",
- atom1, atom2, new_nbrs[c].d,
- new_nbrs[c].dvec[0], new_nbrs[c].dvec[1],
- new_nbrs[c].dvec[2] ); */
- }
- }
- }
-
- ++itr;
- }
-
- Set_End_Index( atom1, num_far, far_nbrs );
- }
- }
- }
- }
-
- far_nbrs->total_intrs = num_far;
- fprintf( stderr, "nbrs done, num_far: %d\n", num_far );
-
-#if defined(DEBUG)
- for ( i = 0; i < system->N; ++i )
- {
- qsort( &(far_nbrs->select.far_nbr_list[ Start_Index(i, far_nbrs) ]),
- Num_Entries(i, far_nbrs), sizeof(far_neighbor_data),
- compare_far_nbrs );
- }
-
- fprintf( stderr, "step%d: num of farnbrs=%6d\n", data->step, num_far );
- fprintf( stderr, "\tallocated farnbrs: %6d\n",
- system->N * far_nbrs->intrs_per_unit );
-#endif
+ t_elapsed = Get_Timing_Info( t_start );
+ data->timing.nbrs += t_elapsed;
}
#endif
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index 65e62467..c837063d 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -959,28 +959,3 @@ void Print_Bond_List2( reax_system *system, reax_list *bonds, char *fname )
fprintf( f, "\n" );
}
}
-
-
-#ifdef LGJ
-Print_XYZ_Serial( reax_system* system, static_storage *workspace )
-{
- rvec p;
- char fname[100];
- FILE *fout;
- int i;
-
- snprintf( fname, 100, "READ_PDB.0" );
- fout = sfopen( fname, "w" );
-
- for ( i = 0; i < system->N; i++ )
- {
- fprintf( fout, "%6d%24.15e%24.15e%24.15e\n",
- workspace->orig_id[i],
- p[0] = system->atoms[i].x[0],
- p[1] = system->atoms[i].x[1],
- p[2] = system->atoms[i].x[2] );
- }
-
- sfclose( fout, "Print_XYZ_Serial::fout" );
-}
-#endif
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 6793dfb9..2be51420 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -36,36 +36,78 @@
#include <omp.h>
#endif
+/* enables debugging code */
//#define DEBUG_FOCUS
+/* enables test forces code */
//#define TEST_FORCES
+/* enables test energy code */
//#define TEST_ENERGY
-#define REORDER_ATOMS /* turns on nbrgen opt by re-ordering atoms */
-//#define LGJ
+/* enables reordering atoms after neighbor list generation (optimization) */
+#define REORDER_ATOMS
+/* enables support for small simulation boxes (i.e. a simulation box with any
+ * dimension less than twice the Verlet list cutoff distance, vlist_cut) */
+//#define SMALL_BOX_SUPPORT
#define SUCCESS (1)
#define FAILURE (0)
#define TRUE (1)
#define FALSE (0)
-#define LOG log
-#define EXP exp
-#define SQRT sqrt
-#define POW pow
-#define ACOS acos
-#define COS cos
-#define SIN sin
-#define TAN tan
-#define CEIL ceil
-#define FLOOR floor
-#define FABS fabs
-#define FMOD fmod
-
-#define SQR(x) ((x)*(x))
-#define CUBE(x) ((x)*(x)*(x))
-#define DEG2RAD(a) ((a)*PI/180.0)
-#define RAD2DEG(a) ((a)*180.0/PI)
-#define MAX( x, y ) (((x) > (y)) ? (x) : (y))
-#define MIN( x, y ) (((x) < (y)) ? (x) : (y))
+/* transcendental constant pi */
+#if defined(M_PI)
+ /* GNU C library (libc), defined in math.h */
+ #define PI (M_PI)
+#else
+ #define PI (3.14159265)
+#endif
+#define C_ele (332.06371)
+//#define K_B (503.398008) // kcal/mol/K
+#define K_B (0.831687) // amu A^2 / ps^2 / K
+#define F_CONV (1e6 / 48.88821291 / 48.88821291) // --> amu A / ps^2
+#define E_CONV (0.002391) // amu A^2 / ps^2 --> kcal/mol
+#define EV_to_KCALpMOL (14.400000) // ElectronVolt --> KCAL per MOLe
+#define KCALpMOL_to_EV (23.060549) // 23.020000//KCAL per MOLe --> ElectronVolt
+#define ECxA_to_DEBYE (4.803204) // elem. charge * angstrom -> debye conv
+#define CAL_to_JOULES (4.184000) // CALories --> JOULES
+#define JOULES_to_CAL (1.0 / 4.184000) // JOULES --> CALories
+#define AMU_to_GRAM (1.6605e-24)
+#define ANG_to_CM (1.0e-8)
+#define AVOGNR (6.0221367e23)
+#define P_CONV (1.0e-24 * AVOGNR * JOULES_to_CAL)
+
+#define MAX_STR (1024)
+#define MAX_LINE (1024)
+#define MAX_TOKENS (1024)
+#define MAX_TOKEN_LEN (1024)
+
+#define MAX_ATOM_ID (100000)
+#define MAX_RESTRICT (15)
+#define MAX_MOLECULE_SIZE (20)
+#define MAX_ATOM_TYPES (25)
+
+#define MAX_GRID (50)
+#define MAX_3BODY_PARAM (5)
+#define MAX_4BODY_PARAM (5)
+#define NO_OF_INTERACTIONS (10)
+
+#define MAX_dV (1.01)
+#define MIN_dV (0.99)
+#define MAX_dT (4.00)
+#define MIN_dT (0.00)
+
+#define ZERO (0.000000000000000e+00)
+#define ALMOST_ZERO (1e-10)
+#define NEG_INF (-1e10)
+#define NO_BOND (1e-3)
+#define HB_THRESHOLD (1e-2)
+#define MAX_BONDS (40)
+#define MIN_BONDS (15)
+#define MIN_HBONDS (50)
+#define SAFE_HBONDS (1.4)
+#define MIN_GCELL_POPL (50)
+#define SAFE_ZONE (1.2)
+#define DANGER_ZONE (0.95)
+#define LOOSE_ZONE (0.75)
/* NaN IEEE 754 representation for C99 in math.h
* Note: function choice must match REAL typedef below */
@@ -75,58 +117,24 @@
#warn "No support for NaN"
#define IS_NAN_REAL(a) (0)
#endif
-
-#if !defined(PI)
- #define PI (3.14159265)
-#endif
-#define C_ele (332.06371)
-//#define K_B (503.398008) // kcal/mol/K
-#define K_B (0.831687) // amu A^2 / ps^2 / K
-#define F_CONV (1e6 / 48.88821291 / 48.88821291) // --> amu A / ps^2
-#define E_CONV (0.002391) // amu A^2 / ps^2 --> kcal/mol
-#define EV_to_KCALpMOL (14.400000) // ElectronVolt --> KCAL per MOLe
-#define KCALpMOL_to_EV (23.060549) // 23.020000//KCAL per MOLe --> ElectronVolt
-#define ECxA_to_DEBYE (4.803204) // elem. charge * angstrom -> debye conv
-#define CAL_to_JOULES (4.184000) // CALories --> JOULES
-#define JOULES_to_CAL (1.0 / 4.184000) // JOULES --> CALories
-#define AMU_to_GRAM (1.6605e-24)
-#define ANG_to_CM (1.0e-8)
-#define AVOGNR (6.0221367e23)
-#define P_CONV (1.0e-24 * AVOGNR * JOULES_to_CAL)
-
-#define MAX_STR (1024)
-#define MAX_LINE (1024)
-#define MAX_TOKENS (1024)
-#define MAX_TOKEN_LEN (1024)
-
-#define MAX_ATOM_ID (100000)
-#define MAX_RESTRICT (15)
-#define MAX_MOLECULE_SIZE (20)
-#define MAX_ATOM_TYPES (25)
-
-#define MAX_GRID (50)
-#define MAX_3BODY_PARAM (5)
-#define MAX_4BODY_PARAM (5)
-#define NO_OF_INTERACTIONS (10)
-
-#define MAX_dV (1.01)
-#define MIN_dV (0.99)
-#define MAX_dT (4.00)
-#define MIN_dT (0.00)
-
-#define ZERO (0.000000000000000e+00)
-#define ALMOST_ZERO (1e-10)
-#define NEG_INF (-1e10)
-#define NO_BOND (1e-3)
-#define HB_THRESHOLD (1e-2)
-#define MAX_BONDS (40)
-#define MIN_BONDS (15)
-#define MIN_HBONDS (50)
-#define SAFE_HBONDS (1.4)
-#define MIN_GCELL_POPL (50)
-#define SAFE_ZONE (1.2)
-#define DANGER_ZONE (0.95)
-#define LOOSE_ZONE (0.75)
+#define LOG (log)
+#define EXP (exp)
+#define SQRT (sqrt)
+#define POW (pow)
+#define ACOS (acos)
+#define COS (cos)
+#define SIN (sin)
+#define TAN (tan)
+#define CEIL (ceil)
+#define FLOOR (floor)
+#define FABS (fabs)
+#define FMOD (fmod)
+#define SQR(x) ((x)*(x))
+#define CUBE(x) ((x)*(x)*(x))
+#define DEG2RAD(a) ((a)*PI/180.0)
+#define RAD2DEG(a) ((a)*180.0/PI)
+#define MAX(x,y) (((x) > (y)) ? (x) : (y))
+#define MIN(x,y) (((x) < (y)) ? (x) : (y))
/* ensemble type */
@@ -295,7 +303,6 @@ typedef struct sparse_matrix sparse_matrix;
typedef struct reallocate_data reallocate_data;
typedef struct LR_data LR_data;
typedef struct cubic_spline_coef cubic_spline_coef;
-typedef struct lookup_table lookup_table;
typedef struct LR_lookup_table LR_lookup_table;
typedef struct static_storage static_storage;
typedef struct reax_list reax_list;
@@ -1010,22 +1017,6 @@ struct cubic_spline_coef
};
-struct lookup_table
-{
- real xmin;
- real xmax;
- int n;
- real dx;
- real inv_dx;
- real a;
-
- real m;
- real c;
-
- real *y;
-};
-
-
struct LR_lookup_table
{
real xmin;
diff --git a/sPuReMD/src/reset_utils.c b/sPuReMD/src/reset_utils.c
index a97ae246..aa58850f 100644
--- a/sPuReMD/src/reset_utils.c
+++ b/sPuReMD/src/reset_utils.c
@@ -138,6 +138,7 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
{
tmp = Start_Index( workspace->hbond_index[i], hbonds );
Set_End_Index( workspace->hbond_index[i], tmp, hbonds );
+
/* fprintf( stderr, "i:%d, hbond: %d-%d\n",
i, Start_Index( workspace->hbond_index[i], hbonds ),
End_Index( workspace->hbond_index[i], hbonds ) );*/
--
GitLab