From de461af222350da2cc054eeaede54bc24746f92a Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Mon, 17 Dec 2018 22:45:12 -0500
Subject: [PATCH] sPuReMD: default to vlist_buffer of 2.5 angstroms. Fixes to
cutoffs used in far neighbor list generation.
---
doc/src/manual.tex | 6 +-
sPuReMD/src/box.c | 310 ++++++++++++++++++++++++++--------------
sPuReMD/src/box.h | 12 +-
sPuReMD/src/control.c | 2 +-
sPuReMD/src/neighbors.c | 10 +-
5 files changed, 217 insertions(+), 123 deletions(-)
diff --git a/doc/src/manual.tex b/doc/src/manual.tex
index 069adabc..25b9f517 100644
--- a/doc/src/manual.tex
+++ b/doc/src/manual.tex
@@ -377,14 +377,14 @@ Defaults: 4.0, 0.001, 0.0
\label{sec:vlist_buffer}
\begin{verbatim}
- reneighbor 10
- vlist_buffer 2
+ reneighbor 1
+ vlist_buffer 2.5
\end{verbatim}
PuReMD features delayed neighbor generation by using Verlet lists.
{\tt reneighbor} controls the reneighboring frequency and {\tt vlist\_buffer}
controls the buffer space beyond the maximum ReaxFF interaction cutoff.
-Defaults: 1 (reneighbor every step), 0
+Defaults: 1 (reneighbor every step), 2.5
\subsubsection{charge\_method}
\label{sec:charge_method}
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index ceb87644..aa0ed4eb 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -160,7 +160,7 @@ void Setup_Box( real a, real b, real c, real alpha, real beta, real gamma,
if ( IS_NAN_REAL(a) || IS_NAN_REAL(b) || IS_NAN_REAL(c)
|| IS_NAN_REAL(alpha) || IS_NAN_REAL(beta) || IS_NAN_REAL(gamma) )
{
- fprintf( stderr, "Invalid simulation box boundaries for big box (NaN). Terminating...\n" );
+ fprintf( stderr, "[ERROR] Invalid simulation box boundaries for big box (NaN). Terminating...\n" );
exit( INVALID_INPUT );
}
@@ -357,7 +357,7 @@ real Metric_Product( rvec x1, rvec x2, simulation_box* box )
* vlist_cut. If so, this neighborhood is added to the list of far neighbors.
* Note: Periodic boundary conditions do not apply. */
int Find_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
- simulation_box *box, real cutoff, far_neighbor_data *data )
+ simulation_box *box, real nonb_cut, real vlist_cut, far_neighbor_data *data )
{
int count;
real norm_sqr;
@@ -366,7 +366,7 @@ int Find_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
rvec_ScaledSum( dvec, 1.0, x2, -1.0, x1 );
norm_sqr = rvec_Norm_Sqr( dvec );
- if ( norm_sqr <= SQR( cutoff ) )
+ if ( norm_sqr <= SQR( vlist_cut ) )
{
data->nbr = nbr_atom;
ivec_MakeZero( data->rel_box );
@@ -387,7 +387,7 @@ int Find_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
/* Similar to Find_Non_Periodic_Far_Neighbors but does not
* update the far neighbors list */
int Count_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
- simulation_box *box, real cutoff )
+ simulation_box *box, real nonb_cut, real vlist_cut )
{
real norm_sqr;
rvec d;
@@ -396,7 +396,7 @@ int Count_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
rvec_ScaledSum( d, 1.0, x2, -1.0, x1 );
norm_sqr = rvec_Norm_Sqr( d );
- if ( norm_sqr <= SQR( cutoff ) )
+ if ( norm_sqr <= SQR( vlist_cut ) )
{
count = 1;
}
@@ -414,7 +414,7 @@ int Count_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
* If the periodic distance between x1 and x2 is less than vlist_cut, this
* neighborhood is added to the list of far neighbors. */
int Find_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, int atom, int nbr_atom,
- simulation_box *box, real cutoff, far_neighbor_data *data )
+ simulation_box *box, real nonb_cut, real vlist_cut, far_neighbor_data *data )
{
int i;
real norm_sqr, tmp, count;
@@ -454,7 +454,7 @@ int Find_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, int atom, int nbr_ato
}
}
- if ( norm_sqr <= SQR( cutoff ) )
+ if ( norm_sqr <= SQR( vlist_cut ) )
{
data->nbr = nbr_atom;
ivec_Copy( data->rel_box, rel_box );
@@ -475,7 +475,7 @@ int Find_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, int atom, int nbr_ato
/* Similar to Find_Periodic_Far_Neighbors_Big_Box but does not
* update the far neighbors list */
int Count_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, int atom, int nbr_atom,
- simulation_box *box, real cutoff )
+ simulation_box *box, real nonb_cut, real vlist_cut )
{
real norm_sqr, d, tmp, count;
int i;
@@ -506,7 +506,7 @@ int Count_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, int atom, int nbr_at
}
}
- if ( norm_sqr <= SQR( cutoff ) )
+ if ( norm_sqr <= SQR( vlist_cut ) )
{
count = 1;
}
@@ -528,103 +528,172 @@ int Count_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, int atom, int nbr_at
* 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!!! */
int Find_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, int atom, int nbr_atom,
- simulation_box *box, real cutoff, far_neighbor_data *data )
+ simulation_box *box, real nonb_cut, real vlist_cut, far_neighbor_data *data )
{
int i, j, k, count;
int imax, jmax, kmax;
- real sqr_norm, d_i, d_j, d_k;
+ real sqr_norm, d_i, d_j, d_k, sqr_cutoff2;
count = 0;
+ sqr_cutoff2 = SQR( vlist_cut );
/* determine the max stretch of imaginary boxs in each direction
* to handle periodic boundary conditions correctly */
- imax = (int)(cutoff / box->box_norms[0] + 1);
- jmax = (int)(cutoff / box->box_norms[1] + 1);
- kmax = (int)(cutoff / box->box_norms[2] + 1);
+// imax = (int)(nonb_cut / box->box_norms[0] + 1);
+// jmax = (int)(nonb_cut / box->box_norms[1] + 1);
+// kmax = (int)(nonb_cut / box->box_norms[2] + 1);
+ imax = (int)(2.0 * nonb_cut / box->box_norms[0]);
+ jmax = (int)(2.0 * nonb_cut / box->box_norms[1]);
+ kmax = (int)(2.0 * nonb_cut / box->box_norms[2]);
+
+ /* non-self interactions */
+ if ( atom != nbr_atom )
+ {
+ /* assumption: orthogonal coordinates */
+ for ( i = -imax; i <= imax; ++i )
+ {
+ d_i = x2[0] - x1[0] + i * box->box_norms[0];
- /*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 ( j = -jmax; j <= jmax; ++j )
+ {
+ d_j = x2[1] - x1[1] + j * box->box_norms[1];
- for ( i = -imax; i <= imax; ++i )
- {
- d_i = (x2[0] + i * box->box_norms[0]) - x1[0];
+ for ( k = -kmax; k <= kmax; ++k )
+ {
+ d_k = x2[2] - x1[2] + k * box->box_norms[2];
+
+ sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
- if ( FABS(d_i) <= cutoff )
+ if ( sqr_norm <= sqr_cutoff2 )
+ {
+ data[count].nbr = nbr_atom;
+
+ data[count].rel_box[0] = i;
+ data[count].rel_box[1] = j;
+ data[count].rel_box[2] = k;
+
+// if ( i )
+// {
+// data[count].ext_factor[0] = (real)i / -abs(i);
+// }
+// else
+// {
+// data[count].ext_factor[0] = 0;
+// }
+//
+// if ( j )
+// {
+// data[count].ext_factor[1] = (real)j / -abs(j);
+// }
+// else
+// {
+// data[count].ext_factor[1] = 0;
+// }
+//
+// if ( k )
+// {
+// data[count].ext_factor[2] = (real)k / -abs(k);
+// }
+// else
+// {
+// data[count].ext_factor[2] = 0;
+// }
+//
+// if ( i == 0 && j == 0 && k == 0 )
+// {
+// data[count].imaginary = 0;
+// }
+// else
+// {
+// data[count].imaginary = 1;
+// }
+
+ data[count].d = SQRT( sqr_norm );
+
+ data[count].dvec[0] = d_i;
+ data[count].dvec[1] = d_j;
+ data[count].dvec[2] = d_k;
+
+ ++count;
+ }
+ }
+ }
+ }
+ }
+ /* self interactions */
+ else
+ {
+ /* assumption: orthogonal coordinates */
+ for ( i = -imax; i <= imax; ++i )
{
+ d_i = x2[0] - x1[0] + i * box->box_norms[0];
+
for ( j = -jmax; j <= jmax; ++j )
{
- d_j = (x2[1] + j * box->box_norms[1]) - x1[1];
+ d_j = x2[1] - x1[1] + j * box->box_norms[1];
- if ( FABS(d_j) <= cutoff )
+ for ( k = -kmax; k <= kmax; ++k )
{
- for ( k = -kmax; k <= kmax; ++k )
+ d_k = x2[2] - x1[2] + k * box->box_norms[2];
+
+ sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
+
+ if ( sqr_norm <= sqr_cutoff2 && sqr_norm >= 0.01 )
{
- d_k = (x2[2] + k * box->box_norms[2]) - x1[2];
-
- if ( FABS(d_k) <= cutoff )
- {
- sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
-
- if ( sqr_norm <= SQR(cutoff) && (atom != nbr_atom
- || (atom == nbr_atom && SQRT( sqr_norm ) >= 0.1)) )
- {
- data[count].nbr = nbr_atom;
-
- data[count].rel_box[0] = i;
- data[count].rel_box[1] = j;
- data[count].rel_box[2] = k;
-
-// if ( i )
-// {
-// data[count].ext_factor[0] = (real)i / -abs(i);
-// }
-// else
-// {
-// data[count].ext_factor[0] = 0;
-// }
-//
-// if ( j )
-// {
-// data[count].ext_factor[1] = (real)j / -abs(j);
-// }
-// else
-// {
-// data[count].ext_factor[1] = 0;
-// }
-//
-// if ( k )
-// {
-// data[count].ext_factor[2] = (real)k / -abs(k);
-// }
-// else
-// {
-// data[count].ext_factor[2] = 0;
-// }
-//
-// if ( i == 0 && j == 0 && k == 0 )
-// {
-// data[count].imaginary = 0;
-// }
-// else
-// {
-// data[count].imaginary = 1;
-// }
-
- data[count].d = SQRT( sqr_norm );
-
- data[count].dvec[0] = d_i;
- data[count].dvec[1] = d_j;
- data[count].dvec[2] = d_k;
-
- ++count;
- }
- }
+ data[count].nbr = nbr_atom;
+
+ data[count].rel_box[0] = i;
+ data[count].rel_box[1] = j;
+ data[count].rel_box[2] = k;
+
+// if ( i )
+// {
+// data[count].ext_factor[0] = (real)i / -abs(i);
+// }
+// else
+// {
+// data[count].ext_factor[0] = 0;
+// }
+//
+// if ( j )
+// {
+// data[count].ext_factor[1] = (real)j / -abs(j);
+// }
+// else
+// {
+// data[count].ext_factor[1] = 0;
+// }
+//
+// if ( k )
+// {
+// data[count].ext_factor[2] = (real)k / -abs(k);
+// }
+// else
+// {
+// data[count].ext_factor[2] = 0;
+// }
+//
+// if ( i == 0 && j == 0 && k == 0 )
+// {
+// data[count].imaginary = 0;
+// }
+// else
+// {
+// data[count].imaginary = 1;
+// }
+
+ data[count].d = SQRT( sqr_norm );
+
+ data[count].dvec[0] = d_i;
+ data[count].dvec[1] = d_j;
+ data[count].dvec[2] = d_k;
+
+ ++count;
}
}
}
}
+
}
return count;
@@ -634,46 +703,69 @@ int Find_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, int atom, int nbr_a
/* Similar to Find_Periodic_Far_Neighbors_Small_Box but does not
* update the far neighbors list */
int Count_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, int atom, int nbr_atom,
- simulation_box *box, real cutoff )
+ simulation_box *box, real nonb_cut, real vlist_cut )
{
int i, j, k, count;
int imax, jmax, kmax;
- real sqr_norm, d_i, d_j, d_k;
+ real sqr_norm, d_i, d_j, d_k, sqr_cutoff2;
count = 0;
+ sqr_cutoff2 = SQR( vlist_cut );
/* determine the max stretch of imaginary boxs in each direction
* to handle periodic boundary conditions correctly */
- imax = (int)(cutoff / box->box_norms[0] + 1);
- jmax = (int)(cutoff / box->box_norms[1] + 1);
- kmax = (int)(cutoff / box->box_norms[2] + 1);
-
- for ( i = -imax; i <= imax; ++i )
+// imax = (int)(nonb_cut / box->box_norms[0] + 1);
+// jmax = (int)(nonb_cut / box->box_norms[1] + 1);
+// kmax = (int)(nonb_cut / box->box_norms[2] + 1);
+ imax = (int)(2.0 * nonb_cut / box->box_norms[0]);
+ jmax = (int)(2.0 * nonb_cut / box->box_norms[1]);
+ kmax = (int)(2.0 * nonb_cut / box->box_norms[2]);
+
+ /* non-self interactions */
+ if ( atom != nbr_atom )
{
- d_i = (x2[0] + i * box->box_norms[0]) - x1[0];
+ for ( i = -imax; i <= imax; ++i )
+ {
+ d_i = x2[0] - x1[0] + i * box->box_norms[0];
+
+ for ( j = -jmax; j <= jmax; ++j )
+ {
+ d_j = x2[1] - x1[1] + j * box->box_norms[1];
- if ( FABS(d_i) <= cutoff )
+ for ( k = -kmax; k <= kmax; ++k )
+ {
+ d_k = x2[2] - x1[2] + k * box->box_norms[2];
+
+ sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
+
+ if ( sqr_norm <= sqr_cutoff2 )
+ {
+ ++count;
+ }
+ }
+ }
+ }
+ }
+ /* self interactions */
+ else
+ {
+ for ( i = -imax; i <= imax; ++i )
{
+ d_i = x2[0] - x1[0] + i * box->box_norms[0];
+
for ( j = -jmax; j <= jmax; ++j )
{
- d_j = (x2[1] + j * box->box_norms[1]) - x1[1];
+ d_j = x2[1] - x1[1] + j * box->box_norms[1];
- if ( FABS(d_j) <= cutoff )
+ for ( k = -kmax; k <= kmax; ++k )
{
- for ( k = -kmax; k <= kmax; ++k )
+ d_k = x2[2] - x1[2] + k * box->box_norms[2];
+
+ sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
+
+ if ( sqr_norm <= sqr_cutoff2 && sqr_norm >= 0.01 )
{
- d_k = (x2[2] + k * box->box_norms[2]) - x1[2];
-
- if ( FABS(d_k) <= cutoff )
- {
- sqr_norm = SQR(d_i) + SQR(d_j) + SQR(d_k);
-
- if ( sqr_norm <= SQR(cutoff) && (atom != nbr_atom
- || (atom == nbr_atom && SQRT( sqr_norm ) >= 0.1)) )
- {
- ++count;
- }
- }
+ ++count;
}
}
}
diff --git a/sPuReMD/src/box.h b/sPuReMD/src/box.h
index 1f110266..b4a3a99c 100644
--- a/sPuReMD/src/box.h
+++ b/sPuReMD/src/box.h
@@ -37,22 +37,22 @@ void Update_Box_Isotropic( simulation_box*, real );
void Update_Box_Semi_Isotropic( simulation_box*, rvec );
int Find_Non_Periodic_Far_Neighbors( rvec, rvec, int, int,
- simulation_box*, real, far_neighbor_data* );
+ simulation_box*, real, real, far_neighbor_data* );
int Count_Non_Periodic_Far_Neighbors( rvec, rvec, int, int,
- simulation_box*, real );
+ simulation_box*, real, real );
int Find_Periodic_Far_Neighbors_Big_Box( rvec, rvec, int, int,
- simulation_box*, real, far_neighbor_data* );
+ simulation_box*, real, real, far_neighbor_data* );
int Count_Periodic_Far_Neighbors_Big_Box( rvec, rvec, int, int,
- simulation_box*, real );
+ simulation_box*, real, real );
int Find_Periodic_Far_Neighbors_Small_Box( rvec, rvec, int, int,
- simulation_box*, real, far_neighbor_data* );
+ simulation_box*, real, real, far_neighbor_data* );
int Count_Periodic_Far_Neighbors_Small_Box( rvec, rvec, int, int,
- simulation_box*, real );
+ simulation_box*, real, real );
void Distance_on_T3_Gen( rvec, rvec, simulation_box*, rvec );
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index fd1a8e21..6e68b5f3 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -64,7 +64,7 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->nonb_cut = system->reaxprm.gp.l[12];
/* defaults values for other cutoffs */
- control->vlist_cut = control->nonb_cut;
+ control->vlist_cut = control->nonb_cut + 2.5;
control->bond_cut = 5.0;
control->bg_cut = 0.3;
control->thb_cut = 0.001;
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index 60efdd59..fafcfa34 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -41,10 +41,10 @@
* Define the preprocessor definition SMALL_BOX_SUPPORT to enable (in
* reax_types.h). */
typedef int (*count_far_neighbors_function)( rvec, rvec, int, int,
- simulation_box*, real );
+ simulation_box*, real, real );
typedef int (*find_far_neighbors_function)( rvec, rvec, int, int,
- simulation_box*, real, far_neighbor_data* );
+ simulation_box*, real, real, far_neighbor_data* );
static void Choose_Neighbor_Counter( reax_system *system, control_params *control,
@@ -186,7 +186,8 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
{
count = Count_Far_Neighbors( system->atoms[atom1].x,
system->atoms[atom2].x, atom1, atom2,
- &system->box, control->vlist_cut );
+ &system->box, control->nonb_cut,
+ control->vlist_cut );
num_far += count;
}
@@ -281,7 +282,8 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
count = Find_Far_Neighbors( system->atoms[atom1].x,
system->atoms[atom2].x, atom1, atom2,
- &system->box, control->vlist_cut, nbr_data );
+ &system->box, control->nonb_cut,
+ control->vlist_cut, nbr_data );
num_far += count;
}
--
GitLab