diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index e142f24cc127248ff40903443fc09eab59f5c9b9..ceb87644dfdc18527e5bda19781768d2b12fb146 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -353,70 +353,115 @@ real Metric_Product( rvec x1, rvec x2, simulation_box* box )
}
-/* Similar to Find_Periodic_Far_Neighbors_Big_Box but does not
+/* 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. */
+int Find_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
+ simulation_box *box, real cutoff, far_neighbor_data *data )
+{
+ int count;
+ real norm_sqr;
+ rvec dvec;
+
+ rvec_ScaledSum( dvec, 1.0, x2, -1.0, x1 );
+ norm_sqr = rvec_Norm_Sqr( dvec );
+
+ if ( norm_sqr <= SQR( cutoff ) )
+ {
+ data->nbr = nbr_atom;
+ ivec_MakeZero( data->rel_box );
+// rvec_MakeZero( data->ext_factor );
+ data->d = SQRT( norm_sqr );
+ rvec_Copy( data->dvec, dvec );
+ count = 1;
+ }
+ else
+ {
+ count = 0;
+ }
+
+ return count;
+}
+
+
+/* Similar to Find_Non_Periodic_Far_Neighbors but does not
* update the far neighbors list */
-int Count_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2,
+int Count_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
+ simulation_box *box, real cutoff )
+{
+ real norm_sqr;
+ rvec d;
+ int count;
+
+ rvec_ScaledSum( d, 1.0, x2, -1.0, x1 );
+ norm_sqr = rvec_Norm_Sqr( d );
+
+ if ( norm_sqr <= SQR( cutoff ) )
+ {
+ count = 1;
+ }
+ else
+ {
+ count = 0;
+ }
+
+ return count;
+}
+
+
+/* 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 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 )
{
- real norm_sqr, d, tmp, ret;
int i;
+ real norm_sqr, tmp, count;
+ ivec rel_box;
+// rvec ext_factor;
+ rvec dvec;
norm_sqr = 0.0;
- ret = FALSE;
for ( i = 0; i < 3; i++ )
{
- d = x2[i] - x1[i];
- tmp = SQR( d );
+ dvec[i] = x2[i] - x1[i];
+ tmp = SQR( dvec[i] );
if ( tmp >= SQR( box->box_norms[i] / 2.0 ) )
{
if ( x2[i] > x1[i] )
{
- d -= box->box_norms[i];
+ rel_box[i] = -1;
+// ext_factor[i] = +1.0;
+ dvec[i] -= box->box_norms[i];
}
else
{
- d += box->box_norms[i];
+ rel_box[i] = +1;
+// ext_factor[i] = -1.0;
+ dvec[i] += box->box_norms[i];
}
- norm_sqr += SQR( d );
+ norm_sqr += SQR( dvec[i] );
}
else
{
norm_sqr += tmp;
+ rel_box[i] = 0;
+// ext_factor[i] = 0.0;
}
}
if ( norm_sqr <= SQR( cutoff ) )
{
- ret = TRUE;
- }
-
- return ret;
-}
-
-
-/* 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. */
-int Find_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
- simulation_box *box, real cutoff, far_neighbor_data *data )
-{
- real norm_sqr;
- int count;
-
- rvec_ScaledSum( data[0].dvec, 1.0, x2, -1.0, x1 );
- norm_sqr = rvec_Norm_Sqr( data[0].dvec );
-
- if ( norm_sqr <= SQR( cutoff ) )
- {
+ data->nbr = nbr_atom;
+ ivec_Copy( data->rel_box, rel_box );
+// rvec_Copy( data->ext_factor, ext_factor );
+ data->d = SQRT( norm_sqr );
+ rvec_Copy( data->dvec, dvec );
count = 1;
- data[0].d = SQRT( norm_sqr );
- data[0].nbr = nbr_atom;
-
- ivec_MakeZero( data[0].rel_box );
- // rvec_MakeZero( data[0].ext_factor );
}
else
{
@@ -427,12 +472,10 @@ int Find_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
}
-/* 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 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 )
+/* 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 )
{
real norm_sqr, d, tmp, count;
int i;
@@ -449,33 +492,23 @@ int Find_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, int atom, int nbr_ato
if ( x2[i] > x1[i] )
{
d -= box->box_norms[i];
- data[0].rel_box[i] = -1;
-// data[0].ext_factor[i] = +1;
}
else
{
d += box->box_norms[i];
- data[0].rel_box[i] = +1;
-// data[0].ext_factor[i] = -1;
}
- data[0].dvec[i] = d;
- norm_sqr += SQR(d);
+ norm_sqr += SQR( d );
}
else
{
- data[0].dvec[i] = d;
norm_sqr += tmp;
- data[0].rel_box[i] = 0;
-// data[0].ext_factor[i] = 0;
}
}
if ( norm_sqr <= SQR( cutoff ) )
{
count = 1;
- data[0].d = SQRT( norm_sqr );
- data[0].nbr = nbr_atom;
}
else
{
@@ -537,23 +570,12 @@ int Find_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, int atom, int nbr_a
if ( sqr_norm <= SQR(cutoff) && (atom != nbr_atom
|| (atom == nbr_atom && SQRT( sqr_norm ) >= 0.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;
+ 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 || j || k )
-// {
-// fprintf( stderr, "x1: %.2f %.2f %.2f\n", x1[0], x1[1], x1[2] );
-// fprintf( stderr, "x2: %.2f %.2f %.2f\n", x2[0], x2[1], x2[2] );
-// fprintf( stderr, "d : %8.2f%8.2f%8.2f\n\n", d_i, d_j, d_k );
-// }
-//
// if ( i )
// {
// data[count].ext_factor[0] = (real)i / -abs(i);
@@ -590,6 +612,65 @@ int Find_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, int atom, int nbr_a
// 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;
+}
+
+
+/* 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 )
+{
+ int i, j, k, count;
+ int imax, jmax, kmax;
+ real sqr_norm, d_i, d_j, d_k;
+
+ count = 0;
+
+ /* 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 )
+ {
+ d_i = (x2[0] + i * box->box_norms[0]) - x1[0];
+
+ if ( FABS(d_i) <= cutoff )
+ {
+ for ( j = -jmax; j <= jmax; ++j )
+ {
+ d_j = (x2[1] + j * box->box_norms[1]) - x1[1];
+
+ if ( FABS(d_j) <= cutoff )
+ {
+ for ( k = -kmax; k <= kmax; ++k )
+ {
+ 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;
}
}
diff --git a/sPuReMD/src/box.h b/sPuReMD/src/box.h
index d6ab4551bf75327f58bcb36ef9b9a218dd4a8161..1f110266198f726f29094503a9ee7986c0e1a4aa 100644
--- a/sPuReMD/src/box.h
+++ b/sPuReMD/src/box.h
@@ -36,18 +36,24 @@ void Update_Box( rtensor, simulation_box* );
void Update_Box_Isotropic( simulation_box*, real );
void Update_Box_Semi_Isotropic( simulation_box*, rvec );
-int Count_Periodic_Far_Neighbors_Big_Box( rvec, rvec, simulation_box*, real,
- far_neighbor_data* );
-
int Find_Non_Periodic_Far_Neighbors( rvec, rvec, int, int,
simulation_box*, real, far_neighbor_data* );
+int Count_Non_Periodic_Far_Neighbors( rvec, rvec, int, int,
+ simulation_box*, real );
+
int Find_Periodic_Far_Neighbors_Big_Box( rvec, rvec, int, int,
simulation_box*, real, far_neighbor_data* );
+int Count_Periodic_Far_Neighbors_Big_Box( rvec, rvec, int, int,
+ simulation_box*, real );
+
int Find_Periodic_Far_Neighbors_Small_Box( rvec, rvec, int, int,
simulation_box*, real, far_neighbor_data* );
+int Count_Periodic_Far_Neighbors_Small_Box( rvec, rvec, int, int,
+ simulation_box*, real );
+
void Distance_on_T3_Gen( rvec, rvec, simulation_box*, rvec );
void Inc_on_T3_Gen( rvec, rvec, simulation_box* );
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 61ec6bc3f82f6e157c3c6efe64797e2695d5e2e4..08a5fab0a3cc9389e77f541f66199cc39b25ca99 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -1108,6 +1108,26 @@ static void Calculate_Charges_EE( const reax_system * const system,
}
+/* Get atomic charge q for ACKS2 method
+ */
+static void Calculate_Charges_ACKS2( const reax_system * const system,
+ static_storage * const workspace )
+{
+ int i;
+
+ for ( i = 0; i < system->N; ++i )
+ {
+ system->atoms[i].q = workspace->s[0][i];
+
+#if defined(DEBUG_FOCUS)
+ printf( "atom %4d: %f\n", i, system->atoms[i].q );
+ printf( " x[0]: %10.5f, x[1]: %10.5f, x[2]: %10.5f\n",
+ system->atoms[i].x[0], system->atoms[i].x[1], system->atoms[i].x[2] );
+#endif
+ }
+}
+
+
/* Main driver method for QEq kernel
* 1) init / setup routines for preconditioning of linear solver
* 2) compute preconditioner
@@ -1356,7 +1376,7 @@ static void ACKS2( reax_system * const system, control_params * const control,
data->timing.cm_solver_iters += iters;
- Calculate_Charges_EE( system, workspace );
+ Calculate_Charges_ACKS2( system, workspace );
}
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index f561da8ff57aa430f50e3e04b7b25d52ccd74ccd..c426b836425837973fde112da1834df8d48bde67 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -240,7 +240,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
if ( reax->sbp[i].gamma_w > 0.5 )
{
if ( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3 )
- fprintf( stderr, "[WARNING] Inconsistent vdWaals-parameters\n" \
+ fprintf( stderr, "[WARNING] Inconsistent vdWaals-parameters!\n" \
" Force field parameters for element %s\n" \
" indicate inner wall+shielding, but earlier\n" \
" atoms indicate different vdWaals-method.\n" \
@@ -262,7 +262,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
{
if ( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 2 )
{
- fprintf( stderr, "[WARNING] Inconsistent vdWaals-parameters\n" \
+ fprintf( stderr, "[WARNING] Inconsistent vdWaals-parameters!\n" \
" Force field parameters for element %s\n" \
" indicate inner wall without shielding, but earlier\n" \
" atoms indicate different vdWaals-method.\n" \
@@ -288,9 +288,9 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
if ( reax->sbp[i].gamma_w > 0.5 )
{
if ( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 1 )
- fprintf( stderr, "[WARNING] Inconsistent vdWaals-parameters\n" \
+ fprintf( stderr, "[WARNING] Inconsistent vdWaals-parameters!\n" \
" Force field parameters for element %s\n" \
- " indicate shielding without inner wall, but earlier\n" \
+ " indicate shielding without inner wall, but earlier\n" \
" atoms indicate different vdWaals-method.\n" \
" This may cause division-by-zero errors.\n" \
" Keeping vdWaals-setting for earlier atoms.\n",
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 1f29e7300096e8937382612f7b003597353dd0e3..f8ad6eb1c7b62a77c6437dbba59e85209f59db82 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -340,7 +340,7 @@ static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
dif = r_ij - base;
val = ((t->ele[r].d * dif + t->ele[r].c) * dif + t->ele[r].b) * dif +
t->ele[r].a;
- val *= EV_to_KCALpMOL / C_ele;
+ val *= EV_to_KCALpMOL / C_ELE;
ret = ((i == j) ? 0.5 : 1.0) * val;
break;
@@ -378,6 +378,38 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
{
case QEQ_CM:
case EE_CM:
+ switch ( pos )
+ {
+ case OFF_DIAGONAL:
+ Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
+ Tap = Tap * r_ij + workspace->Tap[5];
+ Tap = Tap * r_ij + workspace->Tap[4];
+ Tap = Tap * r_ij + workspace->Tap[3];
+ Tap = Tap * r_ij + workspace->Tap[2];
+ Tap = Tap * r_ij + workspace->Tap[1];
+ Tap = Tap * r_ij + workspace->Tap[0];
+
+ /* shielding */
+ dr3gamij_1 = ( r_ij * r_ij * r_ij +
+ system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma );
+ dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
+
+ /* i == j: non-periodic self-interaction term
+ * i != j: periodic self-interaction term */
+ ret = ((i == j) ? 0.5 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
+ break;
+
+ case DIAGONAL:
+ ret = system->reaxprm.sbp[system->atoms[i].type].eta;
+ break;
+
+ default:
+ fprintf( stderr, "[ERROR] Invalid matrix position. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+ break;
+
case ACKS2_CM:
switch ( pos )
{
@@ -395,10 +427,14 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma );
dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
+ /* i == j: non-periodic self-interaction term
+ * i != j: periodic self-interaction term */
ret = ((i == j) ? 0.5 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
break;
case DIAGONAL:
+ /* EE parameters in electron-volts */
+// ret = 2.0 * system->reaxprm.sbp[system->atoms[i].type].eta;
ret = system->reaxprm.sbp[system->atoms[i].type].eta;
break;
@@ -424,7 +460,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
sparse_matrix * H, sparse_matrix * H_sp,
int * Htop, int * H_sp_top )
{
- int i, j, pj;
+ int i, j, pj, target, val_flag;
real d, xcut, bond_softness, * X_diag;
switch ( control->charge_method )
@@ -457,36 +493,20 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
break;
case ACKS2_CM:
- X_diag = (real*) scalloc( system->N, sizeof(real),
+ X_diag = smalloc( sizeof(real) * system->N,
"Init_Charge_Matrix_Remaining_Entries::X_diag" );
- H->start[system->N] = *Htop;
- H_sp->start[system->N] = *H_sp_top;
-
for ( i = 0; i < system->N; ++i )
{
- H->j[*Htop] = i;
- H->val[*Htop] = 1.0;
- *Htop = *Htop + 1;
-
- H_sp->j[*H_sp_top] = i;
- H_sp->val[*H_sp_top] = 1.0;
- *H_sp_top = *H_sp_top + 1;
+ X_diag[i] = 0.0;
}
- H->j[*Htop] = system->N;
- H->val[*Htop] = 0.0;
- *Htop = *Htop + 1;
-
- H_sp->j[*H_sp_top] = system->N;
- H_sp->val[*H_sp_top] = 0.0;
- *H_sp_top = *H_sp_top + 1;
-
for ( i = 0; i < system->N; ++i )
{
- H->start[system->N + i + 1] = *Htop;
- H_sp->start[system->N + i + 1] = *H_sp_top;
+ H->start[system->N + i] = *Htop;
+ H_sp->start[system->N + i] = *H_sp_top;
+ /* constraint on ref. value for kinetic energy potential */
H->j[*Htop] = i;
H->val[*Htop] = 1.0;
*Htop = *Htop + 1;
@@ -495,55 +515,94 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
H_sp->val[*H_sp_top] = 1.0;
*H_sp_top = *H_sp_top + 1;
+ /* kinetic energy terms */
for ( pj = Start_Index(i, far_nbrs); pj < End_Index(i, far_nbrs); ++pj )
{
+ /* exclude self-periodic images of atoms for
+ * kinetic energy term because the effective
+ * potential is the same on an atom and its periodic image */
if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
{
j = far_nbrs->far_nbr_list[pj].nbr;
- xcut = ( system->reaxprm.sbp[ system->atoms[i].type ].b_s_acks2
- + system->reaxprm.sbp[ system->atoms[j].type ].b_s_acks2 )
- / 2.0;
+ xcut = 0.5 * ( system->reaxprm.sbp[ system->atoms[i].type ].b_s_acks2
+ + system->reaxprm.sbp[ system->atoms[j].type ].b_s_acks2 );
- if ( far_nbrs->far_nbr_list[pj].d < xcut &&
- far_nbrs->far_nbr_list[pj].d > 0.001 )
+ if ( far_nbrs->far_nbr_list[pj].d < xcut )
{
d = far_nbrs->far_nbr_list[pj].d / xcut;
- bond_softness = system->reaxprm.gp.l[34] * POW( d, 3.0 ) * POW( 1.0 - d, 6.0 );
-
- H->j[*Htop] = system->N + j + 1;
- H->val[*Htop] = MAX( 0.0, bond_softness );
- *Htop = *Htop + 1;
-
- H_sp->j[*H_sp_top] = system->N + j + 1;
- H_sp->val[*H_sp_top] = MAX( 0.0, bond_softness );
- *H_sp_top = *H_sp_top + 1;
-
- X_diag[i] -= bond_softness;
- X_diag[j] -= bond_softness;
+ bond_softness = system->reaxprm.gp.l[34] * POW( d, 3.0 )
+ * POW( 1.0 - d, 6.0 );
+
+ if ( bond_softness > 0.0 )
+ {
+ val_flag = FALSE;
+
+ for ( target = H->start[system->N + i]; target < *Htop; ++target )
+ {
+ if ( H->j[target] == system->N + j )
+ {
+ H->val[target] += bond_softness;
+ val_flag = TRUE;
+ break;
+ }
+ }
+
+ if ( val_flag == FALSE )
+ {
+ H->j[*Htop] = system->N + j;
+ H->val[*Htop] = bond_softness;
+ ++(*Htop);
+ }
+
+ val_flag = FALSE;
+
+ for ( target = H_sp->start[system->N + i]; target < *H_sp_top; ++target )
+ {
+ if ( H_sp->j[target] == system->N + j )
+ {
+ H_sp->val[target] += bond_softness;
+ val_flag = TRUE;
+ break;
+ }
+ }
+
+ if ( val_flag == FALSE )
+ {
+ H_sp->j[*H_sp_top] = system->N + j;
+ H_sp->val[*H_sp_top] = bond_softness;
+ ++(*H_sp_top);
+ }
+
+ X_diag[i] -= bond_softness;
+ X_diag[j] -= bond_softness;
+ }
}
}
}
- H->j[*Htop] = system->N + i + 1;
+ /* placeholders for diagonal entries, to be replaced below */
+ H->j[*Htop] = system->N + i;
H->val[*Htop] = 0.0;
*Htop = *Htop + 1;
- H_sp->j[*H_sp_top] = system->N + i + 1;
+ H_sp->j[*H_sp_top] = system->N + i;
H_sp->val[*H_sp_top] = 0.0;
*H_sp_top = *H_sp_top + 1;
}
- H->start[system->N_cm - 1] = *Htop;
- H_sp->start[system->N_cm - 1] = *H_sp_top;
+ /* second to last row */
+ H->start[system->N_cm - 2] = *Htop;
+ H_sp->start[system->N_cm - 2] = *H_sp_top;
- for ( i = system->N + 1; i < system->N_cm - 1; ++i )
+ /* place accumulated diagonal entries (needed second to last row marker above before this code) */
+ for ( i = system->N; i < 2 * system->N; ++i )
{
for ( pj = H->start[i]; pj < H->start[i + 1]; ++pj )
{
if ( H->j[pj] == i )
{
- H->val[pj] = X_diag[i - system->N - 1];
+ H->val[pj] = X_diag[i - system->N];
break;
}
}
@@ -552,13 +611,38 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
{
if ( H_sp->j[pj] == i )
{
- H_sp->val[pj] = X_diag[i - system->N - 1];
+ H_sp->val[pj] = X_diag[i - system->N];
break;
}
}
}
- for ( i = system->N + 1; i < system->N_cm - 1; ++i )
+ /* coupling with the kinetic energy potential */
+ for ( i = 0; i < system->N; ++i )
+ {
+ H->j[*Htop] = system->N + i;
+ H->val[*Htop] = 1.0;
+ *Htop = *Htop + 1;
+
+ H_sp->j[*H_sp_top] = system->N + i;
+ H_sp->val[*H_sp_top] = 1.0;
+ *H_sp_top = *H_sp_top + 1;
+ }
+
+ /* explicitly store zero on diagonal */
+ H->j[*Htop] = system->N_cm - 2;
+ H->val[*Htop] = 0.0;
+ *Htop = *Htop + 1;
+
+ H_sp->j[*H_sp_top] = system->N_cm - 2;
+ H_sp->val[*H_sp_top] = 0.0;
+ *H_sp_top = *H_sp_top + 1;
+
+ /* last row */
+ H->start[system->N_cm - 1] = *Htop;
+ H_sp->start[system->N_cm - 1] = *H_sp_top;
+
+ for ( i = 0; i < system->N; ++i )
{
H->j[*Htop] = i;
H->val[*Htop] = 1.0;
@@ -569,6 +653,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
*H_sp_top = *H_sp_top + 1;
}
+ /* explicitly store zero on diagonal */
H->j[*Htop] = system->N_cm - 1;
H->val[*Htop] = 0.0;
*Htop = *Htop + 1;
@@ -577,6 +662,7 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
H_sp->val[*H_sp_top] = 0.0;
*H_sp_top = *H_sp_top + 1;
+
sfree( X_diag, "Init_Charge_Matrix_Remaining_Entries::X_diag" );
break;
@@ -591,12 +677,13 @@ static void Init_Forces( reax_system *system, control_params *control,
reax_list **lists, output_controls *out_control )
{
int i, j, pj;
+ int target;
int start_i, end_i;
int type_i, type_j;
int Htop, H_sp_top, btop_i, btop_j, num_bonds, num_hbonds;
int ihb, jhb, ihb_top, jhb_top;
int flag, flag_sp;
- real r_ij, r2;
+ real r_ij, r2, val, val_flag;
real C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
@@ -689,17 +776,48 @@ static void Init_Forces( reax_system *system, control_params *control,
sbp_j = &system->reaxprm.sbp[type_j];
twbp = &system->reaxprm.tbp[type_i][type_j];
- H->j[Htop] = j;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
- workspace, i, j, r_ij, OFF_DIAGONAL );
- ++Htop;
+ val = Init_Charge_Matrix_Entry( system, control,
+ workspace, i, j, r_ij, OFF_DIAGONAL );
+ val_flag = FALSE;
+
+ for ( target = H->start[i]; target < Htop; ++target )
+ {
+ if ( H->j[target] == j )
+ {
+ H->val[target] += val;
+ val_flag = TRUE;
+ break;
+ }
+ }
+
+ if ( val_flag == FALSE )
+ {
+ H->j[Htop] = j;
+ H->val[Htop] = val;
+ ++Htop;
+ }
/* H_sp matrix entry */
if ( flag_sp == TRUE )
{
- H_sp->j[H_sp_top] = j;
- H_sp->val[H_sp_top] = H->val[Htop - 1];
- ++H_sp_top;
+ val_flag = FALSE;
+
+ for ( target = H_sp->start[i]; target < H_sp_top; ++target )
+ {
+ if ( H_sp->j[target] == j )
+ {
+ H_sp->val[target] += val;
+ val_flag = TRUE;
+ break;
+ }
+ }
+
+ if ( val_flag == FALSE )
+ {
+ H_sp->j[H_sp_top] = j;
+ H_sp->val[H_sp_top] = val;
+ ++H_sp_top;
+ }
}
/* hydrogen bond lists */
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index af1d351798c6fb82dcf8f53b66977c04f61fb998..e0899c7084397b2c06ba674558fa821ccda865aa 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -416,16 +416,22 @@ static void Init_Workspace( reax_system *system, control_params *control,
workspace->b[i] = -system->reaxprm.sbp[ system->atoms[i].type ].chi;
}
- workspace->b_s[system->N] = control->cm_q_net;
- workspace->b[system->N] = control->cm_q_net;
-
- for ( i = system->N + 1; i < system->N_cm; ++i )
+ /* Non-zero total charge can lead to unphysical results.
+ * As such, set the ACKS2 reference charge of every atom
+ * to the total charge divided by the number of atoms.
+ * Except for trivial cases, this leads to fractional
+ * reference charges, which is usually not desirable. */
+ for ( i = 0; i < system->N; ++i )
{
- workspace->b_s[i] = 0.0;
+ workspace->b_s[system->N + i] = control->cm_q_net / system->N;
//TODO: check if unused (redundant)
- workspace->b[i] = 0.0;
+ workspace->b[system->N + i] = control->cm_q_net / system->N;
}
+
+ /* system charge defines the total charge constraint */
+ workspace->b_s[system->N_cm] = control->cm_q_net;
+ workspace->b[system->N_cm] = control->cm_q_net;
break;
default:
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index 332b1e14b3e0bdcf0bd54f630360511eae91b3c9..6f6f7139a4309683aee182dcf4373eb60a7d7ad3 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -218,7 +218,7 @@ static void LR_Lookup( LR_lookup_table *t, real r, LR_data *y )
y->CEclmb = ((t->CEclmb[i].d * dif + t->CEclmb[i].c) * dif + t->CEclmb[i].b) * dif +
t->CEclmb[i].a;
- y->H = y->e_ele * EV_to_KCALpMOL / C_ele;
+ y->H = y->e_ele * EV_to_KCALpMOL / C_ELE;
//y->H = ((t->H[i].d*dif + t->H[i].c)*dif + t->H[i].b)*dif + t->H[i].a;
}
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index aab95c343be611fca388e265a1225222311598ed..60efdd59820da98dcdffe4c3c2a673ecacab78a6 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -29,6 +29,85 @@
#include "tool_box.h"
#include "vector.h"
+#include "print_utils.h"
+
+
+/* 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). */
+typedef int (*count_far_neighbors_function)( rvec, rvec, int, int,
+ simulation_box*, real );
+
+typedef int (*find_far_neighbors_function)( rvec, rvec, int, int,
+ simulation_box*, real, far_neighbor_data* );
+
+
+static void Choose_Neighbor_Counter( reax_system *system, control_params *control,
+ count_far_neighbors_function *Count_Far_Neighbors )
+{
+ if ( control->periodic_boundaries )
+ {
+#ifdef SMALL_BOX_SUPPORT
+ 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 )
+ {
+ *Count_Far_Neighbors = &Count_Periodic_Far_Neighbors_Big_Box;
+ }
+ else
+ {
+ *Count_Far_Neighbors = &Count_Periodic_Far_Neighbors_Small_Box;
+ }
+#else
+ *Count_Far_Neighbors = &Count_Periodic_Far_Neighbors_Big_Box;
+#endif
+ }
+ else
+ {
+ *Count_Far_Neighbors = &Count_Non_Periodic_Far_Neighbors;
+ }
+}
+
+
+static void Choose_Neighbor_Finder( reax_system *system, control_params *control,
+ find_far_neighbors_function *Find_Far_Neighbors )
+{
+ if ( control->periodic_boundaries )
+ {
+#ifdef SMALL_BOX_SUPPORT
+ 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 )
+ {
+ *Find_Far_Neighbors = &Find_Periodic_Far_Neighbors_Big_Box;
+ }
+ else
+ {
+ *Find_Far_Neighbors = &Find_Periodic_Far_Neighbors_Small_Box;
+ }
+#else
+ *Find_Far_Neighbors = &Find_Periodic_Far_Neighbors_Big_Box;
+#endif
+ }
+ else
+ {
+ *Find_Far_Neighbors = &Find_Non_Periodic_Far_Neighbors;
+ }
+}
+
+
+#ifdef DEBUG
+static int compare_far_nbrs(const void *v1, const void *v2)
+{
+ return ((*(far_neighbor_data *)v1).nbr - (*(far_neighbor_data *)v2).nbr);
+}
+#endif
+
static inline real DistSqr_to_CP( rvec cp, rvec x )
{
@@ -53,16 +132,18 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
int i, j, k, l, m, itr;
int x, y, z;
int atom1, atom2, max;
- int num_far;
+ int num_far, count;
int *nbr_atoms;
ivec *nbrs;
rvec *nbrs_cp;
grid *g;
- far_neighbor_data nbr_data;
+ count_far_neighbors_function Count_Far_Neighbors;
g = &system->g;
num_far = 0;
+ Choose_Neighbor_Counter( system, control, &Count_Far_Neighbors );
+
Bin_Atoms( system, workspace );
/* first pick up a cell in the grid */
@@ -79,8 +160,8 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
for ( l = 0; l < g->top[i][j][k]; ++l )
{
atom1 = g->atoms[i][j][k][l];
-
itr = 0;
+
while ( nbrs[itr][0] >= 0 )
{
x = nbrs[itr][0];
@@ -93,24 +174,21 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
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 */
+ /* pick up another atom from the neighbor cell;
+ * we have to compare atom1 with its own periodic images as well
+ * in the case of periodic boundary conditions,
+ * hence the 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 )
+ if ( atom1 >= atom2 )
{
- /* assume periodic boundary conditions since it is
- * safe to over-estimate for the non-periodic case */
- if ( Count_Periodic_Far_Neighbors_Big_Box(system->atoms[atom1].x,
- system->atoms[atom2].x,
- &system->box, control->vlist_cut, &nbr_data) )
- {
- ++num_far;
- }
+ count = Count_Far_Neighbors( system->atoms[atom1].x,
+ system->atoms[atom2].x, atom1, atom2,
+ &system->box, control->vlist_cut );
+
+ num_far += count;
}
}
}
@@ -130,53 +208,6 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
}
-/* 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). */
-typedef int (*find_far_neighbors_function)( rvec, rvec, int, int,
- simulation_box*, real, far_neighbor_data* );
-
-
-void Choose_Neighbor_Finder( reax_system *system, control_params *control,
- find_far_neighbors_function *Find_Far_Neighbors )
-{
- if ( control->periodic_boundaries )
- {
-#ifdef SMALL_BOX_SUPPORT
- 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 )
- {
- *Find_Far_Neighbors = &Find_Periodic_Far_Neighbors_Big_Box;
- }
- else
- {
- *Find_Far_Neighbors = &Find_Periodic_Far_Neighbors_Small_Box;
- }
-#else
- *Find_Far_Neighbors = &Find_Periodic_Far_Neighbors_Big_Box;
-#endif
- }
- else
- {
- *Find_Far_Neighbors = &Find_Non_Periodic_Far_Neighbors;
- }
-}
-
-
-#ifdef DEBUG
-int compare_far_nbrs(const void *v1, const void *v2)
-{
- return ((*(far_neighbor_data *)v1).nbr - (*(far_neighbor_data *)v2).nbr);
-}
-#endif
-
-
void Generate_Neighbor_Lists( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
@@ -196,6 +227,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
t_start = Get_Time( );
g = &system->g;
+ num_far = 0;
far_nbrs = lists[FAR_NBRS];
Bin_Atoms( system, workspace );
@@ -206,8 +238,6 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
Choose_Neighbor_Finder( system, control, &Find_Far_Neighbors );
- num_far = 0;
-
/* first pick up a cell in the grid */
for ( i = 0; i < g->ncell[0]; i++ )
{
@@ -237,17 +267,21 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
nbr_atoms = g->atoms[x][y][z];
max = g->top[x][y][z];
- /* 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
+ * in the case of periodic boundary conditions,
+ * hence the equality in the if stmt below */
for ( m = 0; m < max; ++m )
{
atom2 = nbr_atoms[m];
- if ( atom1 > atom2 )
+ if ( atom1 >= atom2 )
{
nbr_data = &far_nbrs->far_nbr_list[num_far];
- count = Find_Far_Neighbors( system->atoms[atom1].x, system->atoms[atom2].x,
- atom1, atom2, &system->box, control->vlist_cut, nbr_data );
+ count = Find_Far_Neighbors( system->atoms[atom1].x,
+ system->atoms[atom2].x, atom1, atom2,
+ &system->box, control->vlist_cut, nbr_data );
num_far += count;
}
@@ -269,6 +303,7 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
if ( num_far > far_nbrs->total_intrs * DANGER_ZONE )
{
workspace->realloc.num_far = num_far;
+
if ( num_far > far_nbrs->total_intrs )
{
fprintf( stderr, "step%d-ran out of space on far_nbrs: top=%d, max=%d",
@@ -287,239 +322,9 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
#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;
-}
-
-
-#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 )
-{
- int i;
-
- for ( i = Start_Index(atom1, near_nbrs); i < End_Index(atom1, near_nbrs); ++i )
- {
- if ( near_nbrs->near_nbr_list[i].nbr == atom2 )
- {
- return FALSE;
- }
- }
-
- return TRUE;
-}
-
-
-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 )
-{
- int i, j, k;
- int x, y, z;
- int *nbr_atoms;
- int atom1, atom2, max;
- int num_far;
- int c, count;
- int grid_top;
- grid *g;
- reax_list *far_nbrs;
- find_far_neighbors_function Find_Far_Neighbors;
- far_neighbor_data *nbr_data;
- 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 == aNPT )
- {
- Update_Grid( system );
- }
-
- Bin_Atoms( system, workspace );
-
-#ifdef REORDER_ATOMS
- //Cluster_Atoms( system, workspace, control );
-#endif
-
- Choose_Neighbor_Finder( system, control, &Find_Far_Neighbors );
-
- num_far = 0;
- num_near = 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];
-
- /* 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_End_Index( atom1, num_far, far_nbrs );
-
- 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->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 )
- {
- atom2 = nbr_atoms[m];
-
- if ( atom1 >= atom2 )
- {
- nbr_data = &far_nbrs->far_nbr_list[num_far];
-
- //top_near1 = End_Index( atom1, near_nbrs );
- //Set_Start_Index( atom1, num_far, far_nbrs );
- count = Find_Far_Neighbors( system->atoms[atom1].x, system->atoms[atom2].x,
- atom1, atom2, &system->box, control->vlist_cut, nbr_data );
-
- num_far += count;
- }
- }
- }
-
- Set_End_Index( atom1, top_far1, far_nbrs );
- }
- }
- }
- }
-
- /* 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 )
- {
- for ( atom1 = 0; atom1 < system->N; ++atom1 )
- {
- if ( workspace->restricted[ atom1 ] )
- {
- top_near1 = End_Index( atom1, near_nbrs );
-
- for ( j = 0; j < workspace->restricted[ atom1 ]; ++j )
- {
- 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 );
-
- 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( system->atoms[ atom1 ].x,
- system->atoms[ atom2 ].x, &system->box, control,
- new_nbrs, &count );
-
- Set_Near_Neighbor( &near_nbrs->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->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 );
- }
- }
-
- Set_End_Index( atom1, top_near1, near_nbrs );
- }
- }
- }
-
- /* verify nbrlists, count total_intrs, sort nearnbrs */
- near_nbrs->total_intrs = 0;
- far_nbrs->total_intrs = 0;
- for ( i = 0; i < system->N - 1; ++i )
- {
- if ( End_Index(i, near_nbrs) > Start_Index(i + 1, near_nbrs) )
- {
- fprintf( stderr,
- "step%3d: nearnbr list of atom%d is overwritten by atom%d\n",
- data->step, i + 1, i );
- exit( RUNTIME_ERROR );
- }
-
- near_nbrs->total_intrs += Num_Entries( i, near_nbrs );
-
- if ( End_Index(i, far_nbrs) > Start_Index(i + 1, far_nbrs) )
- {
- fprintf( stderr,
- "step%3d: farnbr list of atom%d is overwritten by atom%d\n",
- data->step, i + 1, i );
- exit( RUNTIME_ERROR );
- }
-
- far_nbrs->total_intrs += Num_Entries( i, far_nbrs );
- }
-
- for ( i = 0; i < system->N; ++i )
- {
- qsort( &near_nbrs->near_nbr_list[ Start_Index(i, near_nbrs) ],
- Num_Entries(i, near_nbrs), sizeof(near_neighbor_data),
- compare_near_nbrs );
- }
-
-#ifdef TEST_ENERGY
-// for ( i = 0; i < system->N; ++i )
-// {
-// qsort( &far_nbrs->far_nbr_list[ Start_Index(i, far_nbrs) ],
-// Num_Entries(i, far_nbrs), sizeof(far_neighbor_data), compare_far_nbrs );
-// }
-
- fprintf( stderr, "Near neighbors/atom: %d (compare to 150)\n",
- num_near / system->N );
+ Print_Far_Neighbors( system, control, workspace, lists );
#endif
t_elapsed = Get_Timing_Info( t_start );
data->timing.nbrs += t_elapsed;
}
-#endif
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 3355a85c809fe0a42b86ae608a54be30d6611649..2412eee442d0b571b59991bab099bb84779c30df 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -45,8 +45,9 @@
/* 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
+ * dimension less than twice the Verlet list cutoff distance, vlist_cut),
+ * which means multiple periodic images must be searched for interactions */
+#define SMALL_BOX_SUPPORT
#define SUCCESS (1)
#define FAILURE (0)
@@ -60,19 +61,34 @@
#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 C_ELE (332.06371)
+/* kcal/mol/K */
+//#define K_B (503.398008)
+/* amu A^2 / ps^2 / K */
+#define K_B (0.831687)
+/* --> amu A / ps^2 */
+#define F_CONV (1e6 / 48.88821291 / 48.88821291)
+/* amu A^2 / ps^2 --> kcal/mol */
+#define E_CONV (0.002391)
+/* conversion constant from electron volts to kilo-calories per mole */
+#define EV_to_KCALpMOL (14.400000)
+/* conversion constant from kilo-calories per mole to electron volts */
+//#define KCALpMOL_to_EV (23.060549) // value used in LAMMPS
+#define KCALpMOL_to_EV (23.0200000) // value used in ReaxFF Fortran code
+/* elem. charge * angstrom -> debye conv */
+#define ECxA_to_DEBYE (4.803204)
+/* CALories --> JOULES */
+#define CAL_to_JOULES (4.184000)
+/* JOULES --> CALories */
+#define JOULES_to_CAL (1.0 / 4.184000)
+/* */
#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)
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index ab7b7e4cd3c2f487696afc018fa8b6b11340f0b7..337e853e423c99faa350af37478f84e9e2ff262a 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -69,7 +69,7 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
|| (out_control->write_steps > 0
&& data->step % out_control->write_steps == 0) )
{
- Compute_Total_Energy( system, data );
+ Compute_Total_Energy( system, control, data, workspace );
}
}
@@ -220,7 +220,6 @@ int simulate( const void * const handle )
Generate_Neighbor_Lists( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
- //fprintf( stderr, "total: %.2f secs\n", data.timing.nbrs);
Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
@@ -233,7 +232,8 @@ int simulate( const void * const handle )
|| (spmd_handle->out_control->write_steps > 0
&& spmd_handle->data->step % spmd_handle->out_control->write_steps == 0) )
{
- Compute_Total_Energy( spmd_handle->system, spmd_handle->data );
+ Compute_Total_Energy( spmd_handle->system, spmd_handle->control,
+ spmd_handle->data, spmd_handle->workspace );
}
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index 412843d4416e940b9b48b7c3587792b2991f1eb0..8969c7adf03a910cd6f51f7e752485a2b7de4153 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -209,7 +209,8 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
}
-void Compute_Total_Energy( reax_system* system, simulation_data* data )
+void Compute_Total_Energy( reax_system* system, control_params *control,
+ simulation_data* data, static_storage *workspace )
{
int i, type_i;
real e_pol, q;
@@ -217,18 +218,41 @@ void Compute_Total_Energy( reax_system* system, simulation_data* data )
/* Compute Polarization Energy */
e_pol = 0.0;
+ if ( control->charge_method == QEQ_CM
+ || control->charge_method == EE_CM )
+ {
#ifdef _OPENMP
- #pragma omp parallel for default(none) private(q, type_i) shared(system) \
- reduction(+: e_pol) schedule(static)
+ #pragma omp parallel for default(none) private(q, type_i) shared(system) \
+ reduction(+: e_pol) schedule(static)
#endif
- for ( i = 0; i < system->N; i++ )
+ for ( i = 0; i < system->N; i++ )
+ {
+ q = system->atoms[i].q;
+ type_i = system->atoms[i].type;
+
+ e_pol += ( system->reaxprm.sbp[ type_i ].chi * q
+ + (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) )
+ * KCALpMOL_to_EV;
+ }
+ }
+ else if ( control->charge_method == ACKS2_CM )
{
- q = system->atoms[i].q;
- type_i = system->atoms[i].type;
+#ifdef _OPENMP
+ #pragma omp parallel for default(none) private(q, type_i) shared(system) \
+ reduction(+: e_pol) schedule(static)
+#endif
+ for ( i = 0; i < system->N; i++ )
+ {
+ q = system->atoms[i].q;
+ type_i = system->atoms[i].type;
- e_pol += ( system->reaxprm.sbp[ type_i ].chi * q
- + (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) )
- * KCALpMOL_to_EV;
+ /* energy due to first and second order EE parameters */
+ e_pol += KCALpMOL_to_EV * ( system->reaxprm.sbp[ type_i ].chi
+ + system->reaxprm.sbp[ type_i ].eta / 2.0 * q) * q;
+
+ /* energy due to coupling with kinetic energy potential */
+ e_pol += KCALpMOL_to_EV * system->atoms[i].q * workspace->s[0][ system->N + i ];
+ }
}
data->E_Pol = e_pol;
diff --git a/sPuReMD/src/system_props.h b/sPuReMD/src/system_props.h
index 50a81b1b04171fae64cf3ac87061179312270dba..a7842a7cd216ae5307521c9c4d76b3af88f4cab3 100644
--- a/sPuReMD/src/system_props.h
+++ b/sPuReMD/src/system_props.h
@@ -33,7 +33,7 @@ void Compute_Center_of_Mass( reax_system*, simulation_data*, FILE* );
void Compute_Kinetic_Energy( reax_system*, simulation_data* );
-void Compute_Total_Energy( reax_system*, simulation_data* );
+void Compute_Total_Energy( reax_system*, control_params*, simulation_data*, static_storage* );
void Check_Energy( simulation_data* );
diff --git a/sPuReMD/src/two_body_interactions.c b/sPuReMD/src/two_body_interactions.c
index e6316f813285bb5b066556aad0042b99dd11d0ba..8c9c307c9f4aebfad49feceba4a825639002add2 100644
--- a/sPuReMD/src/two_body_interactions.c
+++ b/sPuReMD/src/two_body_interactions.c
@@ -112,27 +112,27 @@ void Bond_Energy( reax_system *system, control_params *control,
/* Stabilisation terminal triple bond */
if ( bo_ij->BO >= 1.00 )
{
- if ( gp37 == 2 ||
- (sbp_i->mass == 12.0000 && sbp_j->mass == 15.9990) ||
- (sbp_j->mass == 12.0000 && sbp_i->mass == 15.9990) )
+ if ( gp37 == 2
+ || (sbp_i->mass == 12.0000 && sbp_j->mass == 15.9990)
+ || (sbp_j->mass == 12.0000 && sbp_i->mass == 15.9990) )
{
- //ba = SQR(bo_ij->BO - 2.50);
- exphu = EXP( -gp7 * SQR(bo_ij->BO - 2.50) );
- //oboa=abo(j1)-boa;
- //obob=abo(j2)-boa;
+ //ba = SQR( bo_ij->BO - 2.5 );
+ exphu = EXP( -gp7 * SQR(bo_ij->BO - 2.5) );
+ //oboa = abo(j1) - boa;
+ //obob = abo(j2) - boa;
exphua1 = EXP(-gp3 * (workspace->total_bond_order[i] - bo_ij->BO));
exphub1 = EXP(-gp3 * (workspace->total_bond_order[j] - bo_ij->BO));
- //ovoab=abo(j1)-aval(it1)+abo(j2)-aval(it2);
+ //ovoab = abo(j1) - aval(it1) + abo(j2) - aval(it2);
exphuov = EXP(gp4 * (workspace->Delta[i] + workspace->Delta[j]));
hulpov = 1.0 / (1.0 + 25.0 * exphuov);
estriph = gp10 * exphu * hulpov * (exphua1 + exphub1);
- //estrain(j1) = estrain(j1) + 0.50*estriph;
- //estrain(j2) = estrain(j2) + 0.50*estriph;
+ //estrain(j1) = estrain(j1) + 0.5 * estriph;
+ //estrain(j2) = estrain(j2) + 0.5 * estriph;
ebond_total += estriph;
decobdbo = gp10 * exphu * hulpov * (exphua1 + exphub1) *
- ( gp3 - 2.0 * gp7 * (bo_ij->BO - 2.50) );
+ ( gp3 - 2.0 * gp7 * (bo_ij->BO - 2.5) );
decobdboua = -gp10 * exphu * hulpov *
(gp3 * exphua1 + 25.0 * gp4 * exphuov * hulpov * (exphua1 + exphub1));
decobdboub = -gp10 * exphu * hulpov *
@@ -146,7 +146,7 @@ void Bond_Energy( reax_system *system, control_params *control,
fprintf( out_control->ebond,
"%6d%6d%24.15e%24.15e%24.15e%24.15e\n",
workspace->orig_id[i], workspace->orig_id[j],
- //i+1, j+1,
+ //i + 1, j + 1,
estriph, decobdbo, decobdboua, decobdboub );
#endif
@@ -193,6 +193,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
real Tap, dTap, dfn13, CEvd, CEclmb;
real dr3gamij_1, dr3gamij_3;
real e_ele, e_vdW, e_core, de_core;
+ real d, xcut, bond_softness, d_bond_softness, effpot_diff;
rvec temp, ext_press;
//rtensor temp_rtensor, total_rtensor;
two_body_parameters *twbp;
@@ -295,10 +296,10 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
tmp = Tap / dr3gamij_3;
- e_ele = self_coef * C_ele * system->atoms[i].q * system->atoms[j].q * tmp;
+ e_ele = self_coef * C_ELE * system->atoms[i].q * system->atoms[j].q * tmp;
e_ele_total += e_ele;
- CEclmb = self_coef * C_ele * system->atoms[i].q * system->atoms[j].q *
+ CEclmb = self_coef * C_ELE * system->atoms[i].q * system->atoms[j].q *
( dTap - Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
if ( control->ensemble == NVE || control->ensemble == nhNVT
@@ -339,30 +340,30 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
rvec_Add( data->ext_press, ext_press );
}
- /*fprintf( stderr, "nonbonded(%d,%d): rel_box (%f %f %f)",
- i,j,nbr_pj->rel_box[0],nbr_pj->rel_box[1],nbr_pj->rel_box[2] );
-
- fprintf( stderr, "force(%f %f %f)", temp[0], temp[1], temp[2] );
-
- fprintf( stderr, "ext_press (%12.6f %12.6f %12.6f)\n",
- data->ext_press[0], data->ext_press[1], data->ext_press[2] );*/
-
- /* This part is intended for a fully-flexible box */
- /* rvec_OuterProduct( temp_rtensor, nbr_pj->dvec,
- system->atoms[i].x );
- rtensor_Scale( total_rtensor,
- F_C * -(CEvd + CEclmb), temp_rtensor );
- rvec_OuterProduct( temp_rtensor,
- nbr_pj->dvec, system->atoms[j].x );
- rtensor_ScaledAdd( total_rtensor,
- F_C * +(CEvd + CEclmb), temp_rtensor );
-
- if( nbr_pj->imaginary )
- // This is an external force due to an imaginary nbr
- rtensor_ScaledAdd( data->flex_bar.P, -1.0, total_rtensor );
- else
- // This interaction is completely internal
- rtensor_Add( data->flex_bar.P, total_rtensor ); */
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr, "nonbonded(%d,%d): rel_box (%f %f %f)",
+ i, j, nbr_pj->rel_box[0],
+ nbr_pj->rel_box[1], nbr_pj->rel_box[2] );
+
+ fprintf( stderr, "force(%f %f %f)", temp[0], temp[1], temp[2] );
+
+ fprintf( stderr, "ext_press (%12.6f %12.6f %12.6f)\n",
+ data->ext_press[0], data->ext_press[1],
+ data->ext_press[2] );
+#endif
+
+// /* This part is intended for a fully-flexible box */
+// rvec_OuterProduct( temp_rtensor, nbr_pj->dvec, system->atoms[i].x );
+// rtensor_Scale( total_rtensor, F_C * -(CEvd + CEclmb), temp_rtensor );
+// rvec_OuterProduct( temp_rtensor, nbr_pj->dvec, system->atoms[j].x );
+// rtensor_ScaledAdd( total_rtensor, F_C * +(CEvd + CEclmb), temp_rtensor );
+//
+// /* This is an external force due to an imaginary nbr */
+// if ( nbr_pj->imaginary )
+// rtensor_ScaledAdd( data->flex_bar.P, -1.0, total_rtensor );
+// /* This interaction is completely internal */
+// else
+// rtensor_Add( data->flex_bar.P, total_rtensor );
}
#ifdef TEST_ENERGY
@@ -391,15 +392,80 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
}
}
}
+
+ //TODO: better integration of ACKS2 code below with above code for performance
+#ifdef _OPENMP
+ #pragma omp barrier
+#endif
+
+ /* contribution to energy and gradients (atoms and cell)
+ * due to geometry-dependent terms in the ACKS2
+ * kinetic energy */
+ if ( control->charge_method == ACKS2_CM )
+ {
+ for ( i = 0; i < system->N; ++i )
+ {
+ for ( pj = Start_Index(i, far_nbrs); pj < End_Index(i, far_nbrs); ++pj )
+ {
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
+ {
+ nbr_pj = &far_nbrs->far_nbr_list[pj];
+ j = nbr_pj->nbr;
+
+ /* kinetic energy terms */
+ xcut = 0.5 * ( system->reaxprm.sbp[ system->atoms[i].type ].b_s_acks2
+ + system->reaxprm.sbp[ system->atoms[j].type ].b_s_acks2 );
+
+ if ( far_nbrs->far_nbr_list[pj].d < xcut )
+ {
+ d = far_nbrs->far_nbr_list[pj].d / xcut;
+ bond_softness = system->reaxprm.gp.l[34] * POW( d, 3.0 )
+ * POW( 1.0 - d, 6.0 );
+
+ if ( bond_softness > 0.0 )
+ {
+ /* Coulombic energy contribution */
+ effpot_diff = workspace->s[0][system->N + i]
+ - workspace->s[0][system->N + j];
+ e_ele = -0.5 * KCALpMOL_to_EV * bond_softness
+ * SQR( effpot_diff );
+ e_ele_total += e_ele;
+
+ /* forces contribution */
+ d_bond_softness = system->reaxprm.gp.l[34]
+ * 3.0 / xcut * POW( d, 2.0 )
+ * POW( 1.0 - d, 5.0 ) * (1.0 - 3.0 * d);
+ d_bond_softness = -0.5 * d_bond_softness
+ * SQR( effpot_diff );
+ d_bond_softness = KCALpMOL_to_EV * d_bond_softness
+ / far_nbrs->far_nbr_list[pj].d;
+
+#ifndef _OPENMP
+ rvec_ScaledAdd( system->atoms[i].f,
+ -d_bond_softness, nbr_pj->dvec );
+ rvec_ScaledAdd( system->atoms[j].f,
+ d_bond_softness, nbr_pj->dvec );
+#else
+ rvec_ScaledAdd( workspace->f_local[tid * system->N + i],
+ -d_bond_softness, nbr_pj->dvec );
+ rvec_ScaledAdd( workspace->f_local[tid * system->N + j],
+ d_bond_softness, nbr_pj->dvec );
+#endif
+ }
+ }
+ }
+ }
+ }
+ }
}
data->E_vdW = e_vdW_total;
data->E_Ele = e_ele_total;
- // sfclose( fout, "vdW_Coulomb_Energy::fout" );
-
- // fprintf( stderr, "nonbonded: ext_press (%24.15e %24.15e %24.15e)\n",
- // data->ext_press[0], data->ext_press[1], data->ext_press[2] );
+#if defined(DEBUG)
+ fprintf( stderr, "nonbonded: ext_press (%24.15e %24.15e %24.15e)\n",
+ data->ext_press[0], data->ext_press[1], data->ext_press[2] );
+#endif
}
@@ -435,7 +501,6 @@ void LR_vdW_Coulomb( reax_system *system, control_params *control,
dTap = dTap * r_ij + 2 * workspace->Tap[2];
dTap += workspace->Tap[1] / r_ij;
-
/* vdWaals calculations */
powr_vdW1 = POW( r_ij, p_vdW1 );
powgi_vdW1 = POW( 1.0 / twbp->gamma_w, p_vdW1 );
@@ -503,14 +568,14 @@ void LR_vdW_Coulomb( reax_system *system, control_params *control,
tmp = Tap / dr3gamij_3;
lr->H = EV_to_KCALpMOL * tmp;
- lr->e_ele = C_ele * tmp;
+ lr->e_ele = C_ELE * tmp;
/* fprintf( stderr,"i:%d(%d), j:%d(%d), gamma:%f,\
Tap:%f, dr3gamij_3:%f, qi: %f, qj: %f\n",
i, system->atoms[i].type, j, system->atoms[j].type,
twbp->gamma, Tap, dr3gamij_3,
system->atoms[i].q, system->atoms[j].q ); */
- lr->CEclmb = C_ele * ( dTap - Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
+ lr->CEclmb = C_ELE * ( dTap - Tap * r_ij / dr3gamij_1 ) / dr3gamij_3;
/* fprintf( stdout, "%d %d\t%g\t%g %g\t%g %g\t%g %g\n",
i+1, j+1, r_ij, e_vdW, CEvd * r_ij,
system->atoms[i].q, system->atoms[j].q, e_ele, CEclmb * r_ij ); */
@@ -575,36 +640,39 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
t = &workspace->LR[tmin][tmax];
/* Cubic Spline Interpolation */
- r = (int)(r_ij * t->inv_dx);
+ r = (int) (r_ij * t->inv_dx);
if ( r == 0 )
{
++r;
}
- base = (real)(r + 1) * t->dx;
+ base = (real) (r + 1) * t->dx;
dif = r_ij - base;
- //fprintf(stderr, "r: %f, i: %d, base: %f, dif: %f\n", r, i, base, dif);
+
+#if defined(DEBUG)
+ fprintf( stderr, "r: %f, i: %d, base: %f, dif: %f\n", r, i, base, dif );
+#endif
if ( update_energies )
{
- e_vdW = ((t->vdW[r].d * dif + t->vdW[r].c) * dif + t->vdW[r].b) * dif +
- t->vdW[r].a;
+ e_vdW = ((t->vdW[r].d * dif + t->vdW[r].c) * dif + t->vdW[r].b)
+ * dif + t->vdW[r].a;
e_vdW *= self_coef;
- e_ele = ((t->ele[r].d * dif + t->ele[r].c) * dif + t->ele[r].b) * dif +
- t->ele[r].a;
+ e_ele = ((t->ele[r].d * dif + t->ele[r].c) * dif + t->ele[r].b)
+ * dif + t->ele[r].a;
e_ele *= self_coef * system->atoms[i].q * system->atoms[j].q;
e_vdW_total += e_vdW;
e_ele_total += e_ele;
}
- CEvd = ((t->CEvd[r].d * dif + t->CEvd[r].c) * dif + t->CEvd[r].b) * dif +
- t->CEvd[r].a;
+ CEvd = ((t->CEvd[r].d * dif + t->CEvd[r].c) * dif + t->CEvd[r].b)
+ * dif + t->CEvd[r].a;
CEvd *= self_coef;
- //CEvd = (3*t->vdW[r].d*dif + 2*t->vdW[r].c)*dif + t->vdW[r].b;
+// CEvd = (3 * t->vdW[r].d * dif + 2 * t->vdW[r].c) * dif + t->vdW[r].b;
- CEclmb = ((t->CEclmb[r].d * dif + t->CEclmb[r].c) * dif + t->CEclmb[r].b) * dif +
- t->CEclmb[r].a;
+ CEclmb = ((t->CEclmb[r].d * dif + t->CEclmb[r].c) * dif + t->CEclmb[r].b)
+ * dif + t->CEclmb[r].a;
CEclmb *= self_coef * system->atoms[i].q * system->atoms[j].q;
if ( control->ensemble == NVE || control->ensemble == nhNVT
@@ -666,21 +734,3 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
data->E_vdW += e_vdW_total;
data->E_Ele += e_ele_total;
}
-
-
-#if defined(OLD)
-/* Linear extrapolation */
-/*p = (r_ij * t->inv_dx;
- r = (int) p;
- prev = &( t->y[r] );
- next = &( t->y[r+1] );
-
- tmp = p - r;
- e_vdW = self_coef * (prev->e_vdW + tmp*(next->e_vdW - prev->e_vdW ));
- CEvd = self_coef * (prev->CEvd + tmp*(next->CEvd - prev->CEvd ));
-
- e_ele = self_coef * (prev->e_ele + tmp*(next->e_ele - prev->e_ele ));
- e_ele = e_ele * system->atoms[i].q * system->atoms[j].q;
- CEclmb = self_coef * (prev->CEclmb+tmp*(next->CEclmb - prev->CEclmb));
- CEclmb = CEclmb * system->atoms[i].q * system->atoms[j].q;*/
-#endif