From b1e52a8b836ed2a15d743ab536e6283ac6fdd4d3 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Thu, 7 Oct 2021 00:03:50 -0400
Subject: [PATCH] sPuReMD: rework memory management logic to better match
PG-PuReMD code. Adapt new memory management logic for multiple simulation
logic (as motivated by memory issues observed with QM/MM code for AMBER
integration work). Other code clean-up.
---
sPuReMD/src/allocate.c | 196 ++---
sPuReMD/src/allocate.h | 9 +-
sPuReMD/src/box.c | 153 ++--
sPuReMD/src/box.h | 22 +-
sPuReMD/src/forces.c | 1372 +++++++++++++++-------------------
sPuReMD/src/forces.h | 11 +-
sPuReMD/src/init_md.c | 266 +++----
sPuReMD/src/integrate.c | 757 ++++++++++++-------
sPuReMD/src/integrate.h | 14 +-
sPuReMD/src/list.c | 25 +-
sPuReMD/src/list.h | 27 +-
sPuReMD/src/neighbors.c | 57 +-
sPuReMD/src/neighbors.h | 8 +-
sPuReMD/src/reax_types.h | 98 ++-
sPuReMD/src/reset_tools.c | 3 +-
sPuReMD/src/spuremd.c | 89 ++-
sPuReMD/src/system_props.c | 9 +-
sPuReMD/src/system_props.h | 6 +-
sPuReMD/src/valence_angles.c | 14 +-
19 files changed, 1617 insertions(+), 1519 deletions(-)
diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index 4f980bc8..c235913a 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -21,6 +21,7 @@
#include "allocate.h"
+#include "grid.h"
#include "list.h"
#include "tool_box.h"
@@ -97,7 +98,8 @@ void PreAllocate_Space( reax_system * const system,
}
-static void Reallocate_Neighbor_List( reax_list *far_nbr_list, int n, int n_max, int num_intrs )
+static void Reallocate_Neighbor_List( reax_list * const far_nbr_list, int n,
+ int n_max, int num_intrs )
{
if ( far_nbr_list->allocated == TRUE )
{
@@ -146,7 +148,6 @@ void Deallocate_Matrix( sparse_matrix * const H )
static void Reallocate_Matrix( sparse_matrix *H, int n, int n_max, int m )
{
Deallocate_Matrix( H );
-
Allocate_Matrix( H, n, n_max, m );
}
@@ -178,177 +179,104 @@ void Initialize_HBond_List( int n, int const * const h_index,
}
-static void Reallocate_Initialize_HBond_List( int n, int num_h, int num_h_max,
- int *h_index, reax_list *hbond_list )
+static void Reallocate_List( reax_list * const list, int n, int n_max,
+ int max_intrs, int type )
{
- int i, num_hbonds, *hb_top;
-
- hb_top = scalloc( n, sizeof(int), __FILE__, __LINE__ );
- num_hbonds = 0;
-
- for ( i = 0; i < n; ++i )
- {
- if ( h_index[i] >= 0 )
- {
- hb_top[i] = MAX( Num_Entries( h_index[i], hbond_list ) * SAFE_HBONDS,
- MIN_HBONDS );
- num_hbonds += hb_top[i];
- }
- }
-
- if ( hbond_list->allocated == TRUE )
+ if ( list->allocated == TRUE )
{
- Delete_List( TYP_HBOND, hbond_list );
+ Delete_List( type, list );
}
- Make_List( num_h, num_h_max, num_hbonds, TYP_HBOND, hbond_list );
-
- Initialize_HBond_List( n, h_index, hb_top, hbond_list );
-
- sfree( hb_top, __FILE__, __LINE__ );
+ Make_List( n, n_max, max_intrs, type, list );
}
-void Initialize_Bond_List( int * const bond_top,
- reax_list * const bond_list )
+void Reallocate_Part1( reax_system * const system, control_params const * const control,
+ static_storage * const workspace, reax_list ** const lists )
{
- int i;
+ int i, j, k;
+ reallocate_data *realloc;
+ grid *g;
- /* find starting indexes for each atom and the total number of bonds */
- for ( i = 1; i < bond_list->n; ++i )
- {
- bond_top[i] += bond_top[i - 1];
- }
+ realloc = &workspace->realloc;
+ g = &system->g;
- Set_Start_Index( 0, 0, bond_list );
- Set_End_Index( 0, 0, bond_list );
- for ( i = 1; i < bond_list->n; ++i )
+ if ( realloc->gcell_atoms > -1 )
{
- Set_Start_Index( i, bond_top[i - 1], bond_list );
- Set_End_Index( i, bond_top[i - 1], bond_list );
- }
-}
-
-
-static void Reallocate_Initialize_Bond_List( int n, int n_max,
- reax_list *bond_list, int *num_bonds, int *est_3body )
-{
- int i;
- int *bond_top;
-
- bond_top = scalloc( n, sizeof(int), __FILE__, __LINE__ );
- *num_bonds = 0;
- *est_3body = 0;
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr, "[INFO] reallocating gcell: g->max_atoms: %d\n", g->max_atoms );
+#endif
- for ( i = 0; i < n; ++i )
- {
- *est_3body += SQR( Num_Entries( i, bond_list ) );
- bond_top[i] = MAX( Num_Entries( i, bond_list ) * 2, MIN_BONDS );
- *num_bonds += bond_top[i];
- }
+ for ( i = 0; i < g->ncell_max[0]; i++ )
+ {
+ for ( j = 0; j < g->ncell_max[1]; j++ )
+ {
+ for ( k = 0; k < g->ncell_max[2]; k++ )
+ {
+ sfree( g->atoms[i][j][k], __FILE__, __LINE__ );
+ g->atoms[i][j][k] = scalloc( workspace->realloc.gcell_atoms,
+ sizeof(int), __FILE__, __LINE__ );
+ }
+ }
+ }
- if ( bond_list->allocated == TRUE )
- {
- Delete_List( TYP_BOND, bond_list );
+ realloc->gcell_atoms = -1;
}
- Make_List( n, n_max, (int) CEIL( *num_bonds * SAFE_ZONE ),
- TYP_BOND, bond_list );
-
- Initialize_Bond_List( bond_top, bond_list );
-
- sfree( bond_top, __FILE__, __LINE__ );
}
-void Reallocate( reax_system * const system, control_params const * const control,
- static_storage * const workspace, reax_list ** const lists,
- int nbr_flag )
+void Reallocate_Part2( reax_system const * const system,
+ control_params const * const control, simulation_data const * const data,
+ static_storage * const workspace, reax_list ** const lists )
{
- int i, j, k;
- int num_bonds, est_3body;
+ int renbr;
reallocate_data *realloc;
- grid *g;
realloc = &workspace->realloc;
- g = &system->g;
+ renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
- if ( realloc->num_far > 0 && nbr_flag == TRUE )
+ if ( renbr == TRUE && realloc->far_nbrs == TRUE )
{
Reallocate_Neighbor_List( lists[FAR_NBRS],
- system->N, system->N_max, realloc->num_far * SAFE_ZONE );
- realloc->num_far = -1;
+ system->N, system->N_max, realloc->total_far_nbrs );
+ Init_List_Indices( lists[FAR_NBRS] );
+
+ realloc->far_nbrs = FALSE;
}
- if ( realloc->Htop > 0 )
+ if ( realloc->cm == TRUE )
{
Reallocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max,
- realloc->Htop * SAFE_ZONE );
- realloc->Htop = -1;
+ realloc->total_cm_entries );
- Deallocate_Matrix( &workspace->L );
- Deallocate_Matrix( &workspace->U );
+ realloc->cm = FALSE;
}
- if ( control->hbond_cut > 0.0 && realloc->hbonds > 0 )
+ if ( realloc->bonds == TRUE )
{
- Reallocate_Initialize_HBond_List( system->N, workspace->num_H,
- workspace->num_H_max, workspace->hbond_index, lists[HBONDS] );
- realloc->hbonds = -1;
- }
+ Reallocate_List( lists[BONDS], system->N, system->N_max,
+ realloc->total_bonds, TYP_BOND );
+ Init_List_Indices( lists[BONDS] );
- num_bonds = est_3body = -1;
- if ( realloc->bonds > 0 )
- {
- Reallocate_Initialize_Bond_List( system->N, system->N_max, lists[BONDS], &num_bonds, &est_3body );
- realloc->bonds = -1;
- realloc->num_3body = MAX( realloc->num_3body, est_3body );
+ realloc->bonds = FALSE;
+ realloc->thbody = TRUE;
}
- if ( realloc->num_3body > 0 )
+ if ( control->hbond_cut > 0.0 && workspace->num_H > 0 && realloc->hbonds == TRUE )
{
- if ( lists[THREE_BODIES]->allocated == TRUE )
- {
- Delete_List( TYP_THREE_BODY, lists[THREE_BODIES] );
- }
-
- if ( num_bonds == -1 )
- {
- num_bonds = lists[BONDS]->total_intrs;
- }
- realloc->num_3body *= SAFE_ZONE;
+ Reallocate_List( lists[HBONDS], workspace->num_H, workspace->num_H_max,
+ realloc->total_hbonds, TYP_HBOND );
+ Init_List_Indices( lists[HBONDS] );
- Make_List( num_bonds, num_bonds, realloc->num_3body,
- TYP_THREE_BODY, lists[THREE_BODIES] );
- realloc->num_3body = -1;
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "[INFO] reallocating 3 bodies\n" );
- fprintf( stderr, "[INFO] reallocated - num_bonds: %d\n", num_bonds );
- fprintf( stderr, "[INFO] reallocated - num_3body: %d\n", realloc->num_3body );
- fprintf( stderr, "[INFO] reallocated 3body memory: %ldMB\n",
- realloc->num_3body * sizeof(three_body_interaction_data) /
- (1024 * 1024) );
-#endif
+ realloc->hbonds = FALSE;
}
- if ( realloc->gcell_atoms > -1 )
+ if ( realloc->thbody == TRUE )
{
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "[INFO] reallocating gcell: g->max_atoms: %d\n", g->max_atoms );
-#endif
+ Reallocate_List( lists[THREE_BODIES], realloc->total_bonds,
+ (int) CEIL( realloc->total_bonds * SAFE_ZONE ),
+ realloc->total_thbodies, TYP_THREE_BODY );
+ Init_List_Indices( lists[THREE_BODIES] );
- for ( i = 0; i < g->ncell_max[0]; i++ )
- {
- for ( j = 0; j < g->ncell_max[1]; j++ )
- {
- for ( k = 0; k < g->ncell_max[2]; k++ )
- {
- sfree( g->atoms[i][j][k], __FILE__, __LINE__ );
- g->atoms[i][j][k] = scalloc( workspace->realloc.gcell_atoms, sizeof(int),
- __FILE__, __LINE__ );
- }
- }
- }
-
- realloc->gcell_atoms = -1;
+ realloc->thbody = FALSE;
}
}
diff --git a/sPuReMD/src/allocate.h b/sPuReMD/src/allocate.h
index 05ee7fdc..f94ded87 100644
--- a/sPuReMD/src/allocate.h
+++ b/sPuReMD/src/allocate.h
@@ -28,9 +28,6 @@
void PreAllocate_Space( reax_system * const, control_params const * const,
static_storage * const, int );
-void Reallocate( reax_system * const, control_params const * const,
- static_storage * const, reax_list ** const, int );
-
void Allocate_Matrix( sparse_matrix * const , int, int, int );
void Deallocate_Matrix( sparse_matrix * const );
@@ -40,5 +37,11 @@ void Initialize_HBond_List( int, int const * const, int * const,
void Initialize_Bond_List( int * const, reax_list * const );
+void Reallocate_Part1( reax_system * const, control_params const * const,
+ static_storage * const, reax_list ** const );
+
+void Reallocate_Part2( reax_system const * const, control_params const * const,
+ simulation_data const * const, static_storage * const, reax_list ** const );
+
#endif
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index 7a20fa9f..bdc558dd 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -243,7 +243,7 @@ void Update_Box_Semi_Isotropic( simulation_box *box, rvec mu )
* the simulation box that this atom was re-mapped into
* during the position update where (0,0,0) represents no mapping
* */
-void Update_Atom_Position_Periodic( rvec x, rvec dx, ivec rel_map, simulation_box *box )
+void Update_Atom_Position_Periodic( rvec x, rvec dx, ivec rel_map, simulation_box const * const box )
{
int i, remapped;
real tmp;
@@ -302,7 +302,7 @@ void Update_Atom_Position_Periodic( rvec x, rvec dx, ivec rel_map, simulation_bo
* x: updated atom position
* rel_map: unused
* */
-void Update_Atom_Position_Non_Periodic( rvec x, rvec dx, ivec rel_map, simulation_box *box )
+void Update_Atom_Position_Non_Periodic( rvec x, rvec dx, ivec rel_map, simulation_box const * const box )
{
int i;
real tmp;
@@ -344,8 +344,9 @@ void Update_Atom_Position_Non_Periodic( rvec x, rvec dx, ivec rel_map, simulatio
* r: displacement vector betweeon the atoms
* return value: distance
*/
-real Compute_Atom_Distance_Periodic( simulation_box* box, rvec x1, rvec x2,
- ivec x1_rel_map, ivec x2_rel_map, ivec x2_rel_box, rvec r )
+real Compute_Atom_Distance_Periodic( simulation_box const * const box,
+ rvec x1, rvec x2, ivec x1_rel_map, ivec x2_rel_map, ivec x2_rel_box,
+ rvec r )
{
int i;
real norm;
@@ -385,7 +386,7 @@ real Compute_Atom_Distance_Periodic( simulation_box* box, rvec x1, rvec x2,
* r: displacement vector betweeon the atoms
* return value: distance
*/
-real Compute_Atom_Distance_Non_Periodic( simulation_box* box, rvec x1, rvec x2,
+real Compute_Atom_Distance_Non_Periodic( simulation_box const * const box, rvec x1, rvec x2,
ivec x1_rel_map, ivec x2_rel_map, ivec x2_rel_box, rvec r )
{
int i;
@@ -425,7 +426,7 @@ void Compute_Atom_Distance_Triclinic( control_params *control,
}
-void Update_Atom_Position_Triclinic( control_params *control, simulation_box *box,
+void Update_Atom_Position_Triclinic( control_params *control, simulation_box * const box,
rvec x, rvec dx, ivec rel_map )
{
rvec xa, dxa;
@@ -468,7 +469,8 @@ 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 vlist_cut, far_neighbor_data *data )
+ simulation_box const * const box, real vlist_cut,
+ far_neighbor_data * const data, int max )
{
int count;
real norm_sqr;
@@ -479,11 +481,14 @@ int Find_Non_Periodic_Far_Neighbors( rvec x1, rvec x2, int atom, int nbr_atom,
if ( norm_sqr <= SQR( vlist_cut ) && norm_sqr >= MIN_NBR_DIST )
{
- data->nbr = nbr_atom;
- ivec_MakeZero( data->rel_box );
-// rvec_MakeZero( data->ext_factor );
- data->d = SQRT( norm_sqr );
- rvec_Copy( data->dvec, dvec );
+ if ( max > 0 )
+ {
+ 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
@@ -498,7 +503,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 vlist_cut )
+ simulation_box const * const box, real vlist_cut )
{
int count;
real norm_sqr;
@@ -525,7 +530,8 @@ 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 vlist_cut, far_neighbor_data *data )
+ simulation_box const * const box, real vlist_cut,
+ far_neighbor_data * const data, int max )
{
int i, count;
real norm_sqr, tmp, sqr_vlist_cut;
@@ -568,11 +574,14 @@ int Find_Periodic_Far_Neighbors_Big_Box( rvec x1, rvec x2, int atom, int nbr_ato
if ( norm_sqr <= sqr_vlist_cut && norm_sqr >= MIN_NBR_DIST )
{
- 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 );
+ if ( max > 0 )
+ {
+ 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;
}
else
@@ -587,7 +596,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 vlist_cut )
+ simulation_box const * const box, real vlist_cut )
{
int i, count;
real norm_sqr, d, tmp;
@@ -640,7 +649,8 @@ 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 vlist_cut, far_neighbor_data *data )
+ simulation_box const * const box, real vlist_cut,
+ far_neighbor_data * const data, int max )
{
int i, j, k, count;
int imax, jmax, kmax;
@@ -676,54 +686,57 @@ int Find_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, int atom, int nbr_a
if ( sqr_norm <= sqr_vlist_cut && sqr_norm >= MIN_NBR_DIST )
{
- 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 );
- assert( data[count].d > 0.0 );
-
- data[count].dvec[0] = d_i;
- data[count].dvec[1] = d_j;
- data[count].dvec[2] = d_k;
+ if ( count < max )
+ {
+ 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 );
+ assert( data[count].d > 0.0 );
+
+ data[count].dvec[0] = d_i;
+ data[count].dvec[1] = d_j;
+ data[count].dvec[2] = d_k;
+ }
++count;
}
@@ -738,7 +751,7 @@ 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 vlist_cut )
+ simulation_box const * const box, real vlist_cut )
{
int i, j, k, count;
int imax, jmax, kmax;
diff --git a/sPuReMD/src/box.h b/sPuReMD/src/box.h
index e9b63d84..d7d81a2a 100644
--- a/sPuReMD/src/box.h
+++ b/sPuReMD/src/box.h
@@ -37,38 +37,38 @@ 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 const * const, real, far_neighbor_data * const, int );
int Count_Non_Periodic_Far_Neighbors( rvec, rvec, int, int,
- simulation_box*, real );
+ simulation_box const * const , real );
int Find_Periodic_Far_Neighbors_Big_Box( rvec, rvec, int, int,
- simulation_box*, real, far_neighbor_data* );
+ simulation_box const * const, real, far_neighbor_data * const, int );
int Count_Periodic_Far_Neighbors_Big_Box( rvec, rvec, int, int,
- simulation_box*, real );
+ simulation_box const * const, real );
int Find_Periodic_Far_Neighbors_Small_Box( rvec, rvec, int, int,
- simulation_box*, real, far_neighbor_data* );
+ simulation_box const * const, real, far_neighbor_data * const, int );
int Count_Periodic_Far_Neighbors_Small_Box( rvec, rvec, int, int,
- simulation_box*, real );
+ simulation_box const * const, real );
void Compute_Atom_Distance_Triclinic( control_params *,
simulation_box *, rvec, rvec, ivec, ivec, ivec, rvec );
void Update_Atom_Position_Triclinic( control_params *,
- simulation_box *, rvec, rvec, ivec );
+ simulation_box * const, rvec, rvec, ivec );
/* These functions assume that the coordinates are in triclinic system */
/* this function returns cartesian norm but triclinic distance vector */
-real Compute_Atom_Distance_Periodic( simulation_box*, rvec, rvec, ivec, ivec, ivec, rvec );
+real Compute_Atom_Distance_Periodic( simulation_box const * const, rvec, rvec, ivec, ivec, ivec, rvec );
-real Compute_Atom_Distance_Non_Periodic( simulation_box*, rvec, rvec, ivec, ivec, ivec, rvec );
+real Compute_Atom_Distance_Non_Periodic( simulation_box const * const, rvec, rvec, ivec, ivec, ivec, rvec );
-void Update_Atom_Position_Periodic( rvec, rvec, ivec, simulation_box* );
+void Update_Atom_Position_Periodic( rvec, rvec, ivec, simulation_box const * const );
-void Update_Atom_Position_Non_Periodic( rvec, rvec, ivec, simulation_box* );
+void Update_Atom_Position_Non_Periodic( rvec, rvec, ivec, simulation_box const * const );
real Metric_Product( rvec, rvec, simulation_box* );
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index dd59f912..dfe47bb0 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -54,7 +54,7 @@ static int compare_bonds( const void *p1, const void *p2 )
#endif
-void Init_Bonded_Force_Functions( control_params *control )
+void Init_Bonded_Force_Functions( control_params * const control )
{
control->intr_funcs[0] = &BO;
control->intr_funcs[1] = &Bonds;
@@ -268,106 +268,8 @@ static void Compute_Total_Force( reax_system *system, control_params *control,
}
-static void Validate_Lists( static_storage *workspace, reax_list **lists, int step, int n,
- int Hmax, int Htop, int num_bonds, int num_hbonds )
-{
- int i, flag;
- reax_list *bonds, *hbonds;
-
- bonds = lists[BONDS];
- hbonds = lists[HBONDS];
-
- /* far neighbors */
- if ( Htop > Hmax * DANGER_ZONE )
- {
- workspace->realloc.Htop = Htop;
- if ( Htop > Hmax )
- {
- fprintf( stderr,
- "[ERROR] step%d - ran out of space on H matrix: Htop=%d, max = %d",
- step, Htop, Hmax );
- exit( INSUFFICIENT_MEMORY );
- }
- }
-
- /* bond list */
- flag = -1;
- workspace->realloc.num_bonds = num_bonds;
- for ( i = 0; i < n - 1; ++i )
- {
- if ( End_Index(i, bonds) >= Start_Index(i + 1, bonds) - 2 )
- {
- workspace->realloc.bonds = 1;
- if ( End_Index(i, bonds) > Start_Index(i + 1, bonds) )
- {
- flag = i;
- }
- }
- }
-
- if ( flag > -1 )
- {
- fprintf( stderr, "[ERROR] step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
- step, flag, End_Index(flag, bonds), Start_Index(flag + 1, bonds) );
- exit( INSUFFICIENT_MEMORY );
- }
-
- if ( End_Index(i, bonds) >= bonds->total_intrs - 2 )
- {
- workspace->realloc.bonds = 1;
-
- if ( End_Index(i, bonds) > bonds->total_intrs )
- {
- fprintf( stderr, "[ERROR] step%d-bondchk failed: i=%d end(i)=%d bond_end=%d\n",
- step, flag, End_Index(i, bonds), bonds->total_intrs );
- exit( INSUFFICIENT_MEMORY );
- }
- }
-
-
- /* hbonds list */
- if ( workspace->num_H > 0 )
- {
- flag = -1;
- workspace->realloc.num_hbonds = num_hbonds;
- for ( i = 0; i < workspace->num_H - 1; ++i )
- {
- if ( Num_Entries(i, hbonds) >=
- (Start_Index(i + 1, hbonds) - Start_Index(i, hbonds)) * DANGER_ZONE )
- {
- workspace->realloc.hbonds = 1;
- if ( End_Index(i, hbonds) > Start_Index(i + 1, hbonds) )
- {
- flag = i;
- }
- }
- }
-
- if ( flag > -1 )
- {
- fprintf( stderr, "[ERROR] step%d-hbondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
- step, flag, End_Index(flag, hbonds), Start_Index(flag + 1, hbonds) );
- exit( INSUFFICIENT_MEMORY );
- }
-
- if ( Num_Entries(i, hbonds) >=
- (hbonds->total_intrs - Start_Index(i, hbonds)) * DANGER_ZONE )
- {
- workspace->realloc.hbonds = 1;
-
- if ( End_Index(i, hbonds) > hbonds->total_intrs )
- {
- fprintf( stderr, "[ERROR] step%d-hbondchk failed: i=%d end(i)=%d hbondend=%d\n",
- step, flag, End_Index(i, hbonds), hbonds->total_intrs );
- exit( INSUFFICIENT_MEMORY );
- }
- }
- }
-}
-
-
-static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
- control_params *control, LR_lookup_table **LR, int i, int j,
+static inline real Init_Charge_Matrix_Entry_Tab( reax_system const * const system,
+ control_params const * const control, LR_lookup_table ** const LR, int i, int j,
real r_ij, MATRIX_ENTRY_POSITION pos )
{
int r;
@@ -501,23 +403,34 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
for ( i = 0; i < system->N; ++i )
{
/* total charge constraint on atoms */
- H->j[*Htop] = i;
- H->val[*Htop] = 1.0;
-
- H_sp->j[*H_sp_top] = i;
- H_sp->val[*H_sp_top] = 1.0;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = i;
+ H->val[*Htop] = 1.0;
+ }
+ ++(*Htop);
- *Htop = *Htop + 1;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = i;
+ H_sp->val[*H_sp_top] = 1.0;
+ }
+ ++(*H_sp_top);
}
- H->j[*Htop] = system->N;
- H->val[*Htop] = 0.0;
- *Htop = *Htop + 1;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N;
+ H->val[*Htop] = 0.0;
+ }
+ ++(*Htop);
- H_sp->j[*H_sp_top] = system->N;
- H_sp->val[*H_sp_top] = 0.0;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
+ ++(*H_sp_top);
}
else
{
@@ -530,24 +443,35 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
j <= system->molec_charge_constraint_ranges[2 * i + 1]; ++j )
{
/* molecule charge constraint on atoms */
- H->j[*Htop] = j - 1;
- H->val[*Htop] = 1.0;
-
- H_sp->j[*H_sp_top] = j - 1;
- H_sp->val[*H_sp_top] = 1.0;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = j - 1;
+ H->val[*Htop] = 1.0;
+ }
+ ++(*Htop);
- *Htop = *Htop + 1;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = j - 1;
+ H_sp->val[*H_sp_top] = 1.0;
+ }
+ ++(*H_sp_top);
}
/* explicit zeros on diagonals */
- H->j[*Htop] = system->N + i;
- H->val[*Htop] = 0.0;
- *Htop = *Htop + 1;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N + i;
+ H->val[*Htop] = 0.0;
+ }
+ ++(*Htop);
- H_sp->j[*H_sp_top] = system->N + i;
- H_sp->val[*H_sp_top] = 0.0;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N + i;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
+ ++(*H_sp_top);
}
}
break;
@@ -566,13 +490,19 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
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;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = i;
+ H->val[*Htop] = 1.0;
+ }
+ ++(*Htop);
- H_sp->j[*H_sp_top] = i;
- H_sp->val[*H_sp_top] = 1.0;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = i;
+ H_sp->val[*H_sp_top] = 1.0;
+ }
+ ++(*H_sp_top);
/* kinetic energy terms */
for ( pj = Start_Index(i, far_nbr_list); pj < End_Index(i, far_nbr_list); ++pj )
@@ -609,8 +539,11 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
if ( val_flag == FALSE )
{
- H->j[*Htop] = system->N + j;
- H->val[*Htop] = bond_softness;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N + j;
+ H->val[*Htop] = bond_softness;
+ }
++(*Htop);
}
@@ -628,8 +561,11 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
if ( val_flag == FALSE )
{
- H_sp->j[*H_sp_top] = system->N + j;
- H_sp->val[*H_sp_top] = bond_softness;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N + j;
+ H_sp->val[*H_sp_top] = bond_softness;
+ }
++(*H_sp_top);
}
@@ -641,13 +577,19 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
}
/* placeholders for diagonal entries, to be replaced below */
- H->j[*Htop] = system->N + i;
- H->val[*Htop] = 0.0;
- *Htop = *Htop + 1;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N + i;
+ H->val[*Htop] = 0.0;
+ }
+ ++(*Htop);
- H_sp->j[*H_sp_top] = system->N + i;
- H_sp->val[*H_sp_top] = 0.0;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N + i;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
+ ++(*H_sp_top);
}
/* second to last row */
@@ -679,23 +621,35 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
/* 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;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N + i;
+ H->val[*Htop] = 1.0;
+ }
+ ++(*Htop);
- H_sp->j[*H_sp_top] = system->N + i;
- H_sp->val[*H_sp_top] = 1.0;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N + i;
+ H_sp->val[*H_sp_top] = 1.0;
+ }
+ ++(*H_sp_top);
}
/* explicitly store zero on diagonal */
- H->j[*Htop] = system->N_cm - 2;
- H->val[*Htop] = 0.0;
- *Htop = *Htop + 1;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N_cm - 2;
+ H->val[*Htop] = 0.0;
+ }
+ ++(*Htop);
- 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;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N_cm - 2;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
+ ++(*H_sp_top);
/* last row */
H->start[system->N_cm - 1] = *Htop;
@@ -703,23 +657,35 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
for ( i = 0; i < system->N; ++i )
{
- H->j[*Htop] = i;
- H->val[*Htop] = 1.0;
- *Htop = *Htop + 1;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = i;
+ H->val[*Htop] = 1.0;
+ }
+ ++(*Htop);
- H_sp->j[*H_sp_top] = i;
- H_sp->val[*H_sp_top] = 1.0;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = i;
+ H_sp->val[*H_sp_top] = 1.0;
+ }
+ ++(*H_sp_top);
}
/* explicitly store zero on diagonal */
- H->j[*Htop] = system->N_cm - 1;
- H->val[*Htop] = 0.0;
- *Htop = *Htop + 1;
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N_cm - 1;
+ H->val[*Htop] = 0.0;
+ }
+ ++(*Htop);
- H_sp->j[*H_sp_top] = system->N_cm - 1;
- H_sp->val[*H_sp_top] = 0.0;
- *H_sp_top = *H_sp_top + 1;
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N_cm - 1;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
+ ++(*H_sp_top);
sfree( X_diag, __FILE__, __LINE__ );
break;
@@ -737,26 +703,26 @@ static void Init_Distance( reax_system const * const system,
{
int i, j, pj;
int start_i, end_i;
- reax_list *far_nbr_list;
+ reax_list *far_nbrs;
- far_nbr_list = lists[FAR_NBRS];
+ far_nbrs = lists[FAR_NBRS];
- for ( i = 0; i < system->N; ++i )
+ for ( i = 0; i < far_nbrs->n; ++i )
{
- start_i = Start_Index( i, far_nbr_list );
- end_i = End_Index( i, far_nbr_list );
+ start_i = Start_Index( i, far_nbrs );
+ end_i = End_Index( i, far_nbrs );
/* update distance and displacement vector between atoms i and j (i-j)
* for the j atom entry in the far nbr list */
for ( pj = start_i; pj < end_i; ++pj )
{
- j = far_nbr_list->far_nbr_list[pj].nbr;
+ j = far_nbrs->far_nbr_list[pj].nbr;
- far_nbr_list->far_nbr_list[pj].d = control->compute_atom_distance(
+ far_nbrs->far_nbr_list[pj].d = control->compute_atom_distance(
&system->box, system->atoms[i].x, system->atoms[j].x,
system->atoms[i].rel_map, system->atoms[j].rel_map,
- far_nbr_list->far_nbr_list[pj].rel_box,
- far_nbr_list->far_nbr_list[pj].dvec );
+ far_nbrs->far_nbr_list[pj].rel_box,
+ far_nbrs->far_nbr_list[pj].dvec );
}
}
}
@@ -775,37 +741,32 @@ static void Init_CM_Half( reax_system const * const system,
int flag, flag_sp, val_flag;
real val;
sparse_matrix *H, *H_sp;
- reax_list *far_nbr_list;
+ reax_list *far_nbrs;
- far_nbr_list = lists[FAR_NBRS];
+ far_nbrs = lists[FAR_NBRS];
H = &workspace->H;
H_sp = &workspace->H_sp;
Htop = 0;
H_sp_top = 0;
- for ( i = 0; i < system->N; ++i )
+ for ( i = 0; i < far_nbrs->n; ++i )
{
- start_i = Start_Index( i, far_nbr_list );
- end_i = End_Index( i, far_nbr_list );
+ start_i = Start_Index( i, far_nbrs );
+ end_i = End_Index( i, far_nbrs );
H->start[i] = Htop;
H_sp->start[i] = H_sp_top;
for ( pj = start_i; pj < end_i; ++pj )
{
- j = far_nbr_list->far_nbr_list[pj].nbr;
+ j = far_nbrs->far_nbr_list[pj].nbr;
flag = FALSE;
flag_sp = FALSE;
-#if defined(QMMM)
-// if ( system->atoms[i].qmmm_mask == TRUE
-// || system->atoms[j].qmmm_mask == TRUE )
-// {
-#endif
- if ( far_nbr_list->far_nbr_list[pj].d <= control->nonb_cut )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
{
flag = TRUE;
- if ( far_nbr_list->far_nbr_list[pj].d <= control->nonb_sp_cut )
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_sp_cut )
{
flag_sp = TRUE;
}
@@ -814,7 +775,7 @@ static void Init_CM_Half( reax_system const * const system,
if ( flag == TRUE )
{
val = Init_Charge_Matrix_Entry( system, control,
- workspace, i, j, far_nbr_list->far_nbr_list[pj].d,
+ workspace, i, j, far_nbrs->far_nbr_list[pj].d,
OFF_DIAGONAL );
val_flag = FALSE;
@@ -830,8 +791,11 @@ static void Init_CM_Half( reax_system const * const system,
if ( val_flag == FALSE )
{
- H->j[Htop] = j;
- H->val[Htop] = val;
+ if ( Htop < H->m )
+ {
+ H->j[Htop] = j;
+ H->val[Htop] = val;
+ }
++Htop;
}
@@ -858,15 +822,15 @@ static void Init_CM_Half( reax_system const * const system,
}
}
}
-#if defined(QMMM)
-// }
-#endif
}
/* diagonal entry */
- H->j[Htop] = i;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
- workspace, i, i, far_nbr_list->far_nbr_list[pj].d, DIAGONAL );
+ if ( Htop < H->m )
+ {
+ H->j[Htop] = i;
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
+ workspace, i, i, 0.0, DIAGONAL );
+ }
++Htop;
H_sp->j[H_sp_top] = i;
@@ -874,403 +838,187 @@ static void Init_CM_Half( reax_system const * const system,
++H_sp_top;
}
- Init_Charge_Matrix_Remaining_Entries( system, control, far_nbr_list,
+ Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
H, H_sp, &Htop, &H_sp_top );
H->start[system->N_cm] = Htop;
H_sp->start[system->N_cm] = H_sp_top;
+
+
+ if ( Htop > H->m )
+ {
+ workspace->realloc.cm = TRUE;
+ }
}
-/* Compute entries of the bonds/hbonds lists and store the lists in full format
- * using the far neighbors list (stored in full format) */
-static void Init_Bond_Full( reax_system const * const system,
+/* Compute the tabulated charge matrix entries and store the matrix in half format
+ * using the far neighbors list (stored in half format)
+ */
+static void Init_CM_Tab_Half( reax_system const * const system,
control_params const * const control,
- static_storage * const workspace, reax_list ** const lists,
- int * const num_bonds, int * const num_hbonds )
+ static_storage * const workspace, reax_list ** const lists )
{
- int i, j, pj;
+ int i, j, pj, target;
int start_i, end_i;
- int type_i, type_j;
- int btop_i, btop_j;
- int ihb, jhb, ihb_top, jhb_top;
- real r_ij, r2;
- real C12, C34, C56;
- real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
- real BO, BO_s, BO_pi, BO_pi2;
- reax_list *far_nbrs, *bonds, *hbonds;
- single_body_parameters *sbp_i, *sbp_j;
- two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
- bond_data *ibond, *jbond;
- bond_order_data *bo_ij, *bo_ji;
+ int Htop, H_sp_top;
+ int flag, flag_sp, val_flag;
+ real val;
+ sparse_matrix *H, *H_sp;
+ reax_list *far_nbrs;
far_nbrs = lists[FAR_NBRS];
- bonds = lists[BONDS];
- hbonds = lists[HBONDS];
- *num_bonds = 0;
- *num_hbonds = 0;
- btop_i = 0;
- btop_j = 0;
+ H = &workspace->H;
+ H_sp = &workspace->H_sp;
+ Htop = 0;
+ H_sp_top = 0;
- for ( i = 0; i < system->N; ++i )
+ for ( i = 0; i < far_nbrs->n; ++i )
{
- type_i = system->atoms[i].type;
start_i = Start_Index( i, far_nbrs );
end_i = End_Index( i, far_nbrs );
- btop_i = End_Index( i, bonds );
- sbp_i = &system->reax_param.sbp[type_i];
+ H->start[i] = Htop;
+ H_sp->start[i] = H_sp_top;
- if ( control->hbond_cut > 0.0 )
+ for ( pj = start_i; pj < end_i; ++pj )
{
- ihb = sbp_i->p_hbond;
+ j = far_nbrs->far_nbr_list[pj].nbr;
+ flag = FALSE;
+ flag_sp = FALSE;
- if ( ihb == H_ATOM )
- {
- ihb_top = End_Index( workspace->hbond_index[i], hbonds );
- }
- else
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
{
- ihb_top = -1;
- }
- }
- else
- {
- ihb = NON_H_BONDING_ATOM;
- ihb_top = -1;
- }
+ flag = TRUE;
- for ( pj = start_i; pj < end_i; ++pj )
- {
- nbr_pj = &far_nbrs->far_nbr_list[pj];
- j = nbr_pj->nbr;
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_sp_cut )
+ {
+ flag_sp = TRUE;
+ }
+ }
-#if defined(QMMM)
- if ( system->atoms[i].qmmm_mask == TRUE
- || system->atoms[j].qmmm_mask == TRUE )
- {
-#endif
- if ( nbr_pj->d <= control->nonb_cut )
+ if ( flag == TRUE )
{
- type_j = system->atoms[j].type;
- sbp_j = &system->reax_param.sbp[type_j];
- twbp = &system->reax_param.tbp[type_i][type_j];
- r_ij = nbr_pj->d;
+ val = Init_Charge_Matrix_Entry_Tab( system, control,
+ workspace->LR, i, j, far_nbrs->far_nbr_list[pj].d,
+ OFF_DIAGONAL );
+ val_flag = FALSE;
-#if defined(QMMM)
- if ( system->atoms[i].qmmm_mask == TRUE
- && system->atoms[j].qmmm_mask == TRUE )
- {
-#endif
- /* Only non-dummy atoms can form bonds */
- if ( system->atoms[i].is_dummy == FALSE
- && system->atoms[j].is_dummy == FALSE )
+ for ( target = H->start[i]; target < Htop; ++target )
{
-
- /* hydrogen bond lists */
- if ( control->hbond_cut > 0.0
- && (ihb == H_ATOM || ihb == H_BONDING_ATOM)
- && nbr_pj->d <= control->hbond_cut )
+ if ( H->j[target] == j )
{
- jhb = sbp_j->p_hbond;
-
- if ( ihb == H_ATOM && jhb == H_BONDING_ATOM )
- {
- hbonds->hbond_list[ihb_top].nbr = j;
- hbonds->hbond_list[ihb_top].scl = 1;
- hbonds->hbond_list[ihb_top].ptr = nbr_pj;
- ++ihb_top;
- ++(*num_hbonds);
- }
- else if ( ihb == H_BONDING_ATOM && jhb == H_ATOM )
- {
- jhb_top = End_Index( workspace->hbond_index[j], hbonds );
- hbonds->hbond_list[jhb_top].nbr = i;
- hbonds->hbond_list[jhb_top].scl = -1;
- hbonds->hbond_list[jhb_top].ptr = nbr_pj;
- Set_End_Index( workspace->hbond_index[j], jhb_top + 1, hbonds );
- ++(*num_hbonds);
- }
+ H->val[target] += val;
+ val_flag = TRUE;
+ break;
}
+ }
- /* uncorrected bond orders */
- if ( nbr_pj->d <= control->bond_cut )
+ if ( val_flag == FALSE )
+ {
+ if ( Htop < H->m )
{
- r2 = SQR( r_ij );
-
- if ( sbp_i->r_s > 0.0 && sbp_j->r_s > 0.0 )
- {
- C12 = twbp->p_bo1 * POW( r_ij / twbp->r_s, twbp->p_bo2 );
- BO_s = (1.0 + control->bo_cut) * EXP( C12 );
- }
- else
- {
- C12 = 0.0;
- BO_s = 0.0;
- }
+ H->j[Htop] = j;
+ H->val[Htop] = val;
+ }
+ ++Htop;
+ }
- if ( sbp_i->r_pi > 0.0 && sbp_j->r_pi > 0.0 )
- {
- C34 = twbp->p_bo3 * POW( r_ij / twbp->r_p, twbp->p_bo4 );
- BO_pi = EXP( C34 );
- }
- else
- {
- C34 = 0.0;
- BO_pi = 0.0;
- }
+ /* H_sp matrix entry */
+ if ( flag_sp == TRUE )
+ {
+ val_flag = FALSE;
- if ( sbp_i->r_pi_pi > 0.0 && sbp_j->r_pi_pi > 0.0 )
- {
- C56 = twbp->p_bo5 * POW( r_ij / twbp->r_pp, twbp->p_bo6 );
- BO_pi2 = EXP( C56 );
- }
- else
+ for ( target = H_sp->start[i]; target < H_sp_top; ++target )
+ {
+ if ( H_sp->j[target] == j )
{
- C56 = 0.0;
- BO_pi2 = 0.0;
+ H_sp->val[target] += val;
+ val_flag = TRUE;
+ break;
}
+ }
- /* Initially BO values are the uncorrected ones, page 1 */
- BO = BO_s + BO_pi + BO_pi2;
-
- if ( BO >= control->bo_cut )
- {
- *num_bonds += 2;
- /****** bonds i-j and j-i ******/
- ibond = &bonds->bond_list[btop_i];
- btop_j = End_Index( j, bonds );
- jbond = &bonds->bond_list[btop_j];
-
- ibond->nbr = j;
- jbond->nbr = i;
- ibond->d = r_ij;
- jbond->d = r_ij;
- rvec_Copy( ibond->dvec, nbr_pj->dvec );
- rvec_Scale( jbond->dvec, -1, nbr_pj->dvec );
- ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
- ivec_Scale( jbond->rel_box, -1, nbr_pj->rel_box );
- ibond->dbond_index = btop_i;
- jbond->dbond_index = btop_i;
- ibond->sym_index = btop_j;
- jbond->sym_index = btop_i;
- ++btop_i;
- Set_End_Index( j, btop_j + 1, bonds );
-
- bo_ij = &ibond->bo_data;
- bo_ij->BO = BO;
- bo_ij->BO_s = BO_s;
- bo_ij->BO_pi = BO_pi;
- bo_ij->BO_pi2 = BO_pi2;
- bo_ji = &jbond->bo_data;
- bo_ji->BO = BO;
- bo_ji->BO_s = BO_s;
- bo_ji->BO_pi = BO_pi;
- bo_ji->BO_pi2 = BO_pi2;
-
- /* Bond Order page2-3, derivative of total bond order prime */
- Cln_BOp_s = twbp->p_bo2 * C12 / r2;
- Cln_BOp_pi = twbp->p_bo4 * C34 / r2;
- Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2;
-
- /* Only dln_BOp_xx wrt. dr_i is stored here, note that
- dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
- rvec_Scale( bo_ij->dln_BOp_s, -bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
- rvec_Scale( bo_ij->dln_BOp_pi, -bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
- rvec_Scale( bo_ij->dln_BOp_pi2,
- -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
- rvec_Scale( bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s );
- rvec_Scale( bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
- rvec_Scale( bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
-
- /* Only dBOp wrt. dr_i is stored here, note that
- dBOp/dr_i = -dBOp/dr_j and all others are 0 */
- rvec_Scale( bo_ij->dBOp, -(bo_ij->BO_s * Cln_BOp_s
- + bo_ij->BO_pi * Cln_BOp_pi
- + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
- rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );
-
- rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
- rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );
-
- bo_ij->BO_s -= control->bo_cut;
- bo_ij->BO -= control->bo_cut;
- bo_ji->BO_s -= control->bo_cut;
- bo_ji->BO -= control->bo_cut;
- workspace->total_bond_order[i] += bo_ij->BO;
- workspace->total_bond_order[j] += bo_ji->BO;
- bo_ij->Cdbo = 0.0;
- bo_ij->Cdbopi = 0.0;
- bo_ij->Cdbopi2 = 0.0;
- bo_ji->Cdbo = 0.0;
- bo_ji->Cdbopi = 0.0;
- bo_ji->Cdbopi2 = 0.0;
-
- Set_End_Index( j, btop_j + 1, bonds );
- }
+ if ( val_flag == FALSE )
+ {
+ H_sp->j[H_sp_top] = j;
+ H_sp->val[H_sp_top] = val;
+ ++H_sp_top;
}
}
-#if defined(QMMM)
- }
-#endif
- }
-#if defined(QMMM)
}
-#endif
}
- Set_End_Index( i, btop_i, bonds );
- if ( ihb == H_ATOM )
+ /* diagonal entry */
+ if ( Htop < H->m )
{
- Set_End_Index( workspace->hbond_index[i], ihb_top, hbonds );
+ H->j[Htop] = i;
+ H->val[Htop] = Init_Charge_Matrix_Entry_Tab( system, control,
+ workspace->LR, i, i, 0.0, DIAGONAL );
}
- }
-}
-
-
-/* Generate bond list (full format), hydrogen bond list (full format),
- * and charge matrix (half symmetric format)
- * from the far neighbors list (with distance updates, if necessary)
- * */
-static void Init_Forces( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
-{
- int renbr;
- int num_bonds, num_hbonds;
- static int dist_done = FALSE, cm_done = FALSE, bonds_done = FALSE;
-
- renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
- num_bonds = 0;
- num_hbonds = 0;
-
- if ( renbr == FALSE && dist_done == FALSE )
- {
- Init_Distance( system, control, lists );
+ ++Htop;
- dist_done = TRUE;
+ H_sp->j[H_sp_top] = i;
+ H_sp->val[H_sp_top] = H->val[Htop - 1];
+ ++H_sp_top;
}
- if ( cm_done == FALSE )
- {
-// if ( workspace->H.format == SYM_HALF_MATRIX )
-// {
- Init_CM_Half( system, control, workspace, lists );
-// }
-// else
-// {
-// Init_CM_Full( system, control, data, workspace, lists, out_control );
-// }
- }
+ Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
+ H, H_sp, &Htop, &H_sp_top );
- if ( bonds_done == FALSE )
+ H->start[system->N_cm] = Htop;
+ H_sp->start[system->N_cm] = H_sp_top;
+
+
+ if ( Htop > H->m )
{
-// if ( lists[FAR_NBRS]->format == HALF_LIST )
-// {
-// Init_Bond_Half( system, control, workspace, lists, &num_bonds, &num_hbonds );
-// }
-// else
-// {
- Init_Bond_Full( system, control, workspace, lists, &num_bonds, &num_hbonds );
-// }
+ workspace->realloc.cm = TRUE;
}
-
-// ret = (workspace->realloc.cm == FALSE
-// && workspace->realloc.bonds == FALSE
-// && workspace->realloc.hbonds == FALSE
-// ? SUCCESS : FAILURE);
-//
-// if ( workspace->realloc.cm == FALSE )
-// {
-// cm_done = TRUE;
-// }
-// if ( workspace->realloc.bonds == FALSE && workspace->realloc.hbonds == FALSE )
-// {
-// bonds_done = TRUE;
-// }
-//
-// if ( ret == SUCCESS )
-// {
- dist_done = FALSE;
- cm_done = FALSE;
- bonds_done = FALSE;
-// }
-
- /* validate lists - decide if reallocation is required! */
- Validate_Lists( workspace, lists,
- data->step, system->N, workspace->H.m, workspace->H.start[system->N_cm],
- num_bonds, num_hbonds );
-
-#if defined(TEST_FORCES)
- /* Calculate_dBO requires a sorted bonds list */
-// for ( i = 0; i < bonds->n; ++i )
-// {
-// if ( Num_Entries(i, bonds) > 0 )
-// {
-// qsort( &bonds->bond_list[Start_Index(i, bonds)],
-// Num_Entries(i, bonds), sizeof(bond_data), compare_bonds );
-// }
-// }
-#endif
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "step%d: Htop = %d, num_bonds = %d, num_hbonds = %d\n",
- data->step, Htop, num_bonds, num_hbonds );
-
-#endif
}
-static void Init_Forces_Tab( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+/* Compute entries of the bonds/hbonds lists and store the lists in full format
+ * using the far neighbors list (stored in full format) */
+static void Init_Bond_Full( reax_system const * const system,
+ control_params const * const control,
+ static_storage * const workspace, reax_list ** const lists )
{
- int i, j, pj, target;
+ int i, j, pj;
int start_i, end_i;
int type_i, type_j;
- int Htop, H_sp_top, btop_i, btop_j, num_bonds, num_hbonds;
+ int btop_i, btop_j;
int ihb, jhb, ihb_top, jhb_top;
- int flag, flag_sp, val_flag, renbr;
- real r_ij, r2, val;
+ int num_bonds, num_hbonds;
+ real r_ij, r2;
real C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
- sparse_matrix *H, *H_sp;
reax_list *far_nbrs, *bonds, *hbonds;
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
- reax_atom *atom_i, *atom_j;
bond_data *ibond, *jbond;
bond_order_data *bo_ij, *bo_ji;
far_nbrs = lists[FAR_NBRS];
bonds = lists[BONDS];
hbonds = lists[HBONDS];
- H = &workspace->H;
- H_sp = &workspace->H_sp;
- Htop = 0;
- H_sp_top = 0;
num_bonds = 0;
num_hbonds = 0;
btop_i = 0;
btop_j = 0;
- renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
- for ( i = 0; i < system->N; ++i )
+ for ( i = 0; i < far_nbrs->n; ++i )
{
- atom_i = &system->atoms[i];
- type_i = atom_i->type;
+ type_i = system->atoms[i].type;
start_i = Start_Index( i, far_nbrs );
end_i = End_Index( i, far_nbrs );
- H->start[i] = Htop;
- H_sp->start[i] = H_sp_top;
btop_i = End_Index( i, bonds );
sbp_i = &system->reax_param.sbp[type_i];
- if ( control->hbond_cut > 0.0 )
+ if ( control->hbond_cut > 0.0 && workspace->num_H > 0 )
{
ihb = sbp_i->p_hbond;
+
if ( ihb == H_ATOM )
{
ihb_top = End_Index( workspace->hbond_index[i], hbonds );
@@ -1286,116 +1034,34 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
ihb_top = -1;
}
- /* update i-j distance - check if j is within cutoff */
for ( pj = start_i; pj < end_i; ++pj )
{
nbr_pj = &far_nbrs->far_nbr_list[pj];
j = nbr_pj->nbr;
- flag = FALSE;
- flag_sp = FALSE;
#if defined(QMMM)
if ( system->atoms[i].qmmm_mask == TRUE
|| system->atoms[j].qmmm_mask == TRUE )
{
#endif
- /* check if reneighboring step --
- * atomic distances just computed via
- * Verlet list, so use current distances */
- if ( renbr == TRUE )
- {
- if ( nbr_pj->d <= control->nonb_cut )
- {
- flag = TRUE;
-
- if ( nbr_pj->d <= control->nonb_sp_cut )
- {
- flag_sp = TRUE;
- }
- }
- }
- /* update atomic distances */
- else
- {
- atom_j = &system->atoms[j];
- nbr_pj->d = control->compute_atom_distance( &system->box,
- atom_i->x, atom_j->x, atom_i->rel_map,
- atom_j->rel_map, nbr_pj->rel_box,
- nbr_pj->dvec );
-
- if ( nbr_pj->d <= control->nonb_cut )
- {
- flag = TRUE;
-
- if ( nbr_pj->d <= control->nonb_sp_cut )
- {
- flag_sp = TRUE;
- }
- }
- }
-
- if ( flag == TRUE )
+ if ( nbr_pj->d <= control->nonb_cut )
{
type_j = system->atoms[j].type;
sbp_j = &system->reax_param.sbp[type_j];
twbp = &system->reax_param.tbp[type_i][type_j];
r_ij = nbr_pj->d;
- 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 )
- {
- 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;
- }
- }
-
#if defined(QMMM)
if ( system->atoms[i].qmmm_mask == TRUE
&& system->atoms[j].qmmm_mask == TRUE )
{
#endif
- /* Only non-dummy atoms can form bonds */
if ( system->atoms[i].is_dummy == FALSE
&& system->atoms[j].is_dummy == FALSE )
{
+
/* hydrogen bond lists */
- if ( control->hbond_cut > 0.0
+ if ( control->hbond_cut > 0.0 && workspace->num_H > 0
&& (ihb == H_ATOM || ihb == H_BONDING_ATOM)
&& nbr_pj->d <= control->hbond_cut )
{
@@ -1403,19 +1069,25 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
if ( ihb == H_ATOM && jhb == H_BONDING_ATOM )
{
- hbonds->hbond_list[ihb_top].nbr = j;
- hbonds->hbond_list[ihb_top].scl = 1;
- hbonds->hbond_list[ihb_top].ptr = nbr_pj;
- ++ihb_top;
+ if ( num_hbonds < hbonds->total_intrs )
+ {
+ hbonds->hbond_list[ihb_top].nbr = j;
+ hbonds->hbond_list[ihb_top].scl = 1;
+ hbonds->hbond_list[ihb_top].ptr = nbr_pj;
+ ++ihb_top;
+ }
++num_hbonds;
}
else if ( ihb == H_BONDING_ATOM && jhb == H_ATOM )
{
- jhb_top = End_Index( workspace->hbond_index[j], hbonds );
- hbonds->hbond_list[jhb_top].nbr = i;
- hbonds->hbond_list[jhb_top].scl = -1;
- hbonds->hbond_list[jhb_top].ptr = nbr_pj;
- Set_End_Index( workspace->hbond_index[j], jhb_top + 1, hbonds );
+ if ( num_hbonds < hbonds->total_intrs )
+ {
+ jhb_top = End_Index( workspace->hbond_index[j], hbonds );
+ hbonds->hbond_list[jhb_top].nbr = i;
+ hbonds->hbond_list[jhb_top].scl = -1;
+ hbonds->hbond_list[jhb_top].ptr = nbr_pj;
+ Set_End_Index( workspace->hbond_index[j], jhb_top + 1, hbonds );
+ }
++num_hbonds;
}
}
@@ -1463,77 +1135,80 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
if ( BO >= control->bo_cut )
{
+ if ( num_bonds < bonds->total_intrs )
+ {
+ /****** bonds i-j and j-i ******/
+ ibond = &bonds->bond_list[btop_i];
+ btop_j = End_Index( j, bonds );
+ jbond = &bonds->bond_list[btop_j];
+
+ ibond->nbr = j;
+ jbond->nbr = i;
+ ibond->d = r_ij;
+ jbond->d = r_ij;
+ rvec_Copy( ibond->dvec, nbr_pj->dvec );
+ rvec_Scale( jbond->dvec, -1, nbr_pj->dvec );
+ ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
+ ivec_Scale( jbond->rel_box, -1, nbr_pj->rel_box );
+ ibond->dbond_index = btop_i;
+ jbond->dbond_index = btop_i;
+ ibond->sym_index = btop_j;
+ jbond->sym_index = btop_i;
+ ++btop_i;
+ Set_End_Index( j, btop_j + 1, bonds );
+
+ bo_ij = &ibond->bo_data;
+ bo_ij->BO = BO;
+ bo_ij->BO_s = BO_s;
+ bo_ij->BO_pi = BO_pi;
+ bo_ij->BO_pi2 = BO_pi2;
+ bo_ji = &jbond->bo_data;
+ bo_ji->BO = BO;
+ bo_ji->BO_s = BO_s;
+ bo_ji->BO_pi = BO_pi;
+ bo_ji->BO_pi2 = BO_pi2;
+
+ /* Bond Order page2-3, derivative of total bond order prime */
+ Cln_BOp_s = twbp->p_bo2 * C12 / r2;
+ Cln_BOp_pi = twbp->p_bo4 * C34 / r2;
+ Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2;
+
+ /* Only dln_BOp_xx wrt. dr_i is stored here, note that
+ dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
+ rvec_Scale( bo_ij->dln_BOp_s, -bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
+ rvec_Scale( bo_ij->dln_BOp_pi, -bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
+ rvec_Scale( bo_ij->dln_BOp_pi2,
+ -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
+ rvec_Scale( bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s );
+ rvec_Scale( bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
+ rvec_Scale( bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
+
+ /* Only dBOp wrt. dr_i is stored here, note that
+ dBOp/dr_i = -dBOp/dr_j and all others are 0 */
+ rvec_Scale( bo_ij->dBOp, -(bo_ij->BO_s * Cln_BOp_s
+ + bo_ij->BO_pi * Cln_BOp_pi
+ + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
+ rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );
+
+ rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
+ rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );
+
+ bo_ij->BO_s -= control->bo_cut;
+ bo_ij->BO -= control->bo_cut;
+ bo_ji->BO_s -= control->bo_cut;
+ bo_ji->BO -= control->bo_cut;
+ workspace->total_bond_order[i] += bo_ij->BO;
+ workspace->total_bond_order[j] += bo_ji->BO;
+ bo_ij->Cdbo = 0.0;
+ bo_ij->Cdbopi = 0.0;
+ bo_ij->Cdbopi2 = 0.0;
+ bo_ji->Cdbo = 0.0;
+ bo_ji->Cdbopi = 0.0;
+ bo_ji->Cdbopi2 = 0.0;
+
+ Set_End_Index( j, btop_j + 1, bonds );
+ }
num_bonds += 2;
- /****** bonds i-j and j-i ******/
- ibond = &bonds->bond_list[btop_i];
- btop_j = End_Index( j, bonds );
- jbond = &bonds->bond_list[btop_j];
-
- ibond->nbr = j;
- jbond->nbr = i;
- ibond->d = r_ij;
- jbond->d = r_ij;
- rvec_Copy( ibond->dvec, nbr_pj->dvec );
- rvec_Scale( jbond->dvec, -1, nbr_pj->dvec );
- ivec_Copy( ibond->rel_box, nbr_pj->rel_box );
- ivec_Scale( jbond->rel_box, -1, nbr_pj->rel_box );
- ibond->dbond_index = btop_i;
- jbond->dbond_index = btop_i;
- ibond->sym_index = btop_j;
- jbond->sym_index = btop_i;
- ++btop_i;
- Set_End_Index( j, btop_j + 1, bonds );
-
- bo_ij = &ibond->bo_data;
- bo_ij->BO = BO;
- bo_ij->BO_s = BO_s;
- bo_ij->BO_pi = BO_pi;
- bo_ij->BO_pi2 = BO_pi2;
- bo_ji = &jbond->bo_data;
- bo_ji->BO = BO;
- bo_ji->BO_s = BO_s;
- bo_ji->BO_pi = BO_pi;
- bo_ji->BO_pi2 = BO_pi2;
-
- /* Bond Order page2-3, derivative of total bond order prime */
- Cln_BOp_s = twbp->p_bo2 * C12 / r2;
- Cln_BOp_pi = twbp->p_bo4 * C34 / r2;
- Cln_BOp_pi2 = twbp->p_bo6 * C56 / r2;
-
- /* Only dln_BOp_xx wrt. dr_i is stored here, note that
- dln_BOp_xx/dr_i = -dln_BOp_xx/dr_j and all others are 0 */
- rvec_Scale( bo_ij->dln_BOp_s, -bo_ij->BO_s * Cln_BOp_s, ibond->dvec );
- rvec_Scale( bo_ij->dln_BOp_pi, -bo_ij->BO_pi * Cln_BOp_pi, ibond->dvec );
- rvec_Scale( bo_ij->dln_BOp_pi2,
- -bo_ij->BO_pi2 * Cln_BOp_pi2, ibond->dvec );
- rvec_Scale( bo_ji->dln_BOp_s, -1., bo_ij->dln_BOp_s );
- rvec_Scale( bo_ji->dln_BOp_pi, -1., bo_ij->dln_BOp_pi );
- rvec_Scale( bo_ji->dln_BOp_pi2, -1., bo_ij->dln_BOp_pi2 );
-
- /* Only dBOp wrt. dr_i is stored here, note that
- dBOp/dr_i = -dBOp/dr_j and all others are 0 */
- rvec_Scale( bo_ij->dBOp, -(bo_ij->BO_s * Cln_BOp_s
- + bo_ij->BO_pi * Cln_BOp_pi
- + bo_ij->BO_pi2 * Cln_BOp_pi2), ibond->dvec );
- rvec_Scale( bo_ji->dBOp, -1., bo_ij->dBOp );
-
- rvec_Add( workspace->dDeltap_self[i], bo_ij->dBOp );
- rvec_Add( workspace->dDeltap_self[j], bo_ji->dBOp );
-
- bo_ij->BO_s -= control->bo_cut;
- bo_ij->BO -= control->bo_cut;
- bo_ji->BO_s -= control->bo_cut;
- bo_ji->BO -= control->bo_cut;
- workspace->total_bond_order[i] += bo_ij->BO;
- workspace->total_bond_order[j] += bo_ji->BO;
- bo_ij->Cdbo = 0.0;
- bo_ij->Cdbopi = 0.0;
- bo_ij->Cdbopi2 = 0.0;
- bo_ji->Cdbo = 0.0;
- bo_ji->Cdbopi = 0.0;
- bo_ji->Cdbopi2 = 0.0;
-
- Set_End_Index( j, btop_j + 1, bonds );
}
}
}
@@ -1546,16 +1221,6 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
#endif
}
- /* diagonal entry */
- H->j[Htop] = i;
- H->val[Htop] = Init_Charge_Matrix_Entry( system, control,
- workspace, i, i, r_ij, DIAGONAL );
- ++Htop;
-
- H_sp->j[H_sp_top] = i;
- H_sp->val[H_sp_top] = H->val[Htop - 1];
- ++H_sp_top;
-
Set_End_Index( i, btop_i, bonds );
if ( ihb == H_ATOM )
{
@@ -1563,28 +1228,190 @@ static void Init_Forces_Tab( reax_system *system, control_params *control,
}
}
- Init_Charge_Matrix_Remaining_Entries( system, control, far_nbrs,
- H, H_sp, &Htop, &H_sp_top );
+ if ( num_bonds > bonds->total_intrs )
+ {
+ workspace->realloc.bonds = TRUE;
+ }
- H->start[system->N_cm] = Htop;
- H_sp->start[system->N_cm] = H_sp_top;
+ if ( control->hbond_cut > 0.0 && workspace->num_H > 0
+ && num_hbonds > hbonds->total_intrs )
+ {
+ workspace->realloc.hbonds = TRUE;
+ }
+}
+
+
+/* Generate bond list (full format), hydrogen bond list (full format),
+ * and charge matrix (half symmetric format)
+ * from the far neighbors list (with distance updates, if necessary)
+ * */
+static int Init_Forces( reax_system * const system,
+ control_params const * const control,
+ simulation_data * const data, static_storage * const workspace,
+ reax_list ** const lists, output_controls * const out_control )
+{
+ int i, renbr, ret;
+ static int dist_done = FALSE, cm_done = FALSE, bonds_done = FALSE;
+
+ renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
+
+ if ( renbr == FALSE && dist_done == FALSE )
+ {
+ Init_Distance( system, control, lists );
+
+ dist_done = TRUE;
+ }
+
+ if ( cm_done == FALSE )
+ {
+ if ( control->tabulate <= 0 )
+ {
+// if ( workspace->H.format == SYM_HALF_MATRIX )
+// {
+ Init_CM_Half( system, control, workspace, lists );
+// }
+// else
+// {
+// Init_CM_Full( system, control, data, workspace, lists, out_control );
+// }
+ }
+ else
+ {
+// if ( workspace->H.format == SYM_HALF_MATRIX )
+// {
+ Init_CM_Tab_Half( system, control, workspace, lists );
+// }
+// else
+// {
+// Init_CM_Tab_Full( system, control, data, workspace, lists, out_control );
+// }
+ }
+ }
+
+ if ( bonds_done == FALSE )
+ {
+ for ( i = 0; i < system->N; ++i )
+ {
+ workspace->total_bond_order[i] = 0.0;
+ }
+ for ( i = 0; i < system->N; ++i )
+ {
+ rvec_MakeZero( workspace->dDeltap_self[i] );
+ }
+
+// if ( lists[FAR_NBRS]->format == HALF_LIST )
+// {
+// Init_Bond_Half( system, control, workspace, lists );
+// }
+// else
+// {
+ Init_Bond_Full( system, control, workspace, lists );
+// }
+ }
+
+ ret = (workspace->realloc.cm == FALSE
+ && workspace->realloc.bonds == FALSE
+ && workspace->realloc.hbonds == FALSE ? SUCCESS : FAILURE);
+
+ if ( workspace->realloc.cm == FALSE )
+ {
+ cm_done = TRUE;
+ }
+ if ( workspace->realloc.bonds == FALSE && workspace->realloc.hbonds == FALSE )
+ {
+ bonds_done = TRUE;
+ }
- /* validate lists - decide if reallocation is required! */
- Validate_Lists( workspace, lists,
- data->step, system->N, H->m, Htop, num_bonds, num_hbonds );
+ if ( ret == SUCCESS )
+ {
+ dist_done = FALSE;
+ cm_done = FALSE;
+ bonds_done = FALSE;
+ }
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "step%d: Htop = %d, num_bonds = %d, num_hbonds = %d\n",
- data->step, Htop, num_bonds, num_hbonds );
- //Print_Bonds( system, bonds, "sbonds.out" );
- //Print_Bond_List2( system, bonds, "sbonds.out" );
- //Print_Sparse_Matrix2( H, "H.out", NULL );
+#if defined(TEST_FORCES)
+ /* Calculate_dBO requires a sorted bonds list */
+// for ( i = 0; i < bonds->n; ++i )
+// {
+// if ( Num_Entries(i, bonds) > 0 )
+// {
+// qsort( &bonds->bond_list[Start_Index(i, bonds)],
+// Num_Entries(i, bonds), sizeof(bond_data), compare_bonds );
+// }
+// }
#endif
+
+ return ret;
}
-void Estimate_Storage_Sizes( reax_system *system, control_params *control,
- reax_list **lists, int *Htop, int *hb_top, int *bond_top, int *num_3body )
+static void Estimate_Storages_CM( reax_system const * const system,
+ control_params const * const control, static_storage * const workspace,
+ reax_list ** const lists )
+{
+ int i, pj;
+ int start_i, end_i;
+ int Htop;
+ reax_list *far_nbrs;
+
+ Htop = 0;
+ far_nbrs = lists[FAR_NBRS];
+
+ for ( i = 0; i < far_nbrs->n; ++i )
+ {
+ start_i = Start_Index( i, far_nbrs );
+ end_i = End_Index( i, far_nbrs );
+
+ /* update i-j distance - check if j is within cutoff */
+ for ( pj = start_i; pj < end_i; ++pj )
+ {
+ if ( far_nbrs->far_nbr_list[pj].d <= control->nonb_cut )
+ {
+ ++Htop;
+ }
+ }
+
+ /* diagonal entry */
+ ++Htop;
+ }
+
+ switch ( control->charge_method )
+ {
+ case QEQ_CM:
+ break;
+
+ case EE_CM:
+ if ( system->num_molec_charge_constraints == 0 )
+ {
+ Htop += system->N_cm;
+ }
+ else
+ {
+ for ( i = 0; i < system->num_molec_charge_constraints; ++i )
+ {
+ Htop += system->molec_charge_constraint_ranges[2 * i + 1]
+ - system->molec_charge_constraint_ranges[2 * i] + 1;
+ }
+ }
+ break;
+
+ case ACKS2_CM:
+ Htop = 2 * Htop + 3 * system->N + 2;
+ break;
+
+ default:
+ fprintf( stderr, "[ERROR] Unknown charge method type. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+
+ workspace->realloc.total_cm_entries = (int) CEIL( Htop * SAFE_ZONE );
+}
+
+
+static void Estimate_Storages_Bonds( reax_system const * const system,
+ control_params const * const control, static_storage * const workspace,
+ reax_list ** const lists )
{
int i, j, pj;
int start_i, end_i;
@@ -1601,6 +1428,15 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
far_nbrs = lists[FAR_NBRS];
+ for ( i = 0; i < system->N_max; ++i )
+ {
+ system->bonds[i] = 0;
+ }
+ for ( i = 0; i < system->N_max; ++i )
+ {
+ system->hbonds[i] = 0;
+ }
+
for ( i = 0; i < far_nbrs->n; ++i )
{
atom_i = &system->atoms[i];
@@ -1622,7 +1458,6 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
type_j = atom_j->type;
sbp_j = &system->reax_param.sbp[type_j];
twbp = &system->reax_param.tbp[type_i][type_j];
- ++(*Htop);
/* hydrogen bond lists */
if ( control->hbond_cut > 0.0
@@ -1633,11 +1468,11 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
if ( ihb == H_ATOM && jhb == H_BONDING_ATOM )
{
- ++hb_top[i];
+ ++system->hbonds[i];
}
else if ( ihb == H_BONDING_ATOM && jhb == H_ATOM )
{
- ++hb_top[j];
+ ++system->hbonds[j];
}
}
@@ -1684,93 +1519,100 @@ void Estimate_Storage_Sizes( reax_system *system, control_params *control,
if ( BO >= control->bo_cut )
{
- ++bond_top[i];
- ++bond_top[j];
+ ++system->bonds[i];
+ ++system->bonds[j];
}
}
}
}
}
- switch ( control->charge_method )
+ workspace->realloc.total_thbodies = 0;
+ for ( i = 0; i < system->N; ++i )
{
- case QEQ_CM:
- break;
+ //TODO: the factor of 2 depends on the bond list format (and also effects thbody list), revisit later
+ system->bonds[i] = MAX( (int) CEIL( system->bonds[i] * 2.0 * SAFE_ZONE ), MIN_BONDS );
+ }
- case EE_CM:
- if ( system->num_molec_charge_constraints == 0 )
- {
- *Htop += system->N_cm;
- }
- else
- {
- for ( i = 0; i < system->num_molec_charge_constraints; ++i )
- {
- *Htop += system->molec_charge_constraint_ranges[2 * i + 1]
- - system->molec_charge_constraint_ranges[2 * i] + 1;
- }
- }
- break;
+ for ( i = 0; i < system->N; ++i )
+ {
+ system->hbonds[i] = MAX( (int) CEIL( system->hbonds[i] * SAFE_HBONDS ), MIN_HBONDS );
+ }
- case ACKS2_CM:
- *Htop = 2 * *Htop + 3 * system->N + 2;
- break;
+ /* reductions to get totals */
+ workspace->realloc.total_bonds = 0;
+ workspace->realloc.total_hbonds = 0;
+ workspace->realloc.total_thbodies = 0;
- default:
- fprintf( stderr, "[ERROR] Unknown charge method type. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
+ for ( i = 0; i < system->N_max; ++i )
+ {
+ workspace->realloc.total_bonds += system->bonds[i];
+ }
+ for ( i = 0; i < system->N_max; ++i )
+ {
+ workspace->realloc.total_hbonds += system->hbonds[i];
}
+ for ( i = 0; i < system->N_max; ++i )
+ {
+ workspace->realloc.total_thbodies += SQR( system->bonds[i] );
+ }
+}
- *Htop *= SAFE_ZONE;
- for ( i = 0; i < system->N; ++i )
+void Estimate_Storages( reax_system const * const system,
+ control_params const * const control, static_storage * const workspace,
+ reax_list ** const lists, int realloc_cm, int realloc_bonds )
+{
+ if ( realloc_cm == TRUE )
+ {
+ Estimate_Storages_CM( system, control, workspace, lists );
+ }
+
+ if ( realloc_bonds == TRUE )
{
- hb_top[i] = MAX( hb_top[i] * SAFE_HBONDS, MIN_HBONDS );
- *num_3body += SQR( bond_top[i] );
- bond_top[i] = MAX( bond_top[i] * 2, MIN_BONDS );
+ Estimate_Storages_Bonds( system, control, workspace, lists );
}
- *num_3body *= SAFE_ZONE;
}
-void Compute_Forces( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list** lists, output_controls *out_control, int realloc )
+int Compute_Forces( reax_system * const system, control_params * const control,
+ simulation_data * const data, static_storage * const workspace,
+ reax_list ** const lists, output_controls * const out_control, int realloc )
{
+ int ret;
real t_start, t_elapsed;
t_start = Get_Time( );
- if ( control->tabulate <= 0 )
- {
- Init_Forces( system, control, data, workspace, lists, out_control );
- }
- else
- {
- Init_Forces_Tab( system, control, data, workspace, lists, out_control );
- }
+ ret = Init_Forces( system, control, data, workspace, lists, out_control );
t_elapsed = Get_Timing_Info( t_start );
data->timing.init_forces += t_elapsed;
- t_start = Get_Time( );
- Compute_Bonded_Forces( system, control, data, workspace, lists, out_control );
- t_elapsed = Get_Timing_Info( t_start );
- data->timing.bonded += t_elapsed;
+ if ( ret != SUCCESS )
+ {
+ Estimate_Storages( system, control, workspace, lists, workspace->realloc.cm,
+ workspace->realloc.bonds || workspace->realloc.hbonds );
+ }
- t_start = Get_Time( );
- Compute_NonBonded_Forces( system, control, data, workspace,
- lists, out_control, realloc );
- t_elapsed = Get_Timing_Info( t_start );
- data->timing.nonb += t_elapsed;
+ if ( ret == SUCCESS )
+ {
+ t_start = Get_Time( );
+ Compute_Bonded_Forces( system, control, data, workspace, lists, out_control );
+ t_elapsed = Get_Timing_Info( t_start );
+ data->timing.bonded += t_elapsed;
- Compute_Total_Force( system, control, data, workspace, lists );
+ t_start = Get_Time( );
+ Compute_NonBonded_Forces( system, control, data, workspace,
+ lists, out_control, realloc );
+ t_elapsed = Get_Timing_Info( t_start );
+ data->timing.nonb += t_elapsed;
-#if defined(DEBUG_FOCUS)
- //Print_Total_Force( system, control, data, workspace, lists, out_control );
-#endif
+ Compute_Total_Force( system, control, data, workspace, lists );
#if defined(TEST_FORCES)
- Print_Total_Force( system, control, data, workspace, lists, out_control );
- Compare_Total_Forces( system, control, data, workspace, lists, out_control );
+ Print_Total_Force( system, control, data, workspace, lists, out_control );
+ Compare_Total_Forces( system, control, data, workspace, lists, out_control );
#endif
+ }
+
+ return ret;
}
diff --git a/sPuReMD/src/forces.h b/sPuReMD/src/forces.h
index 982ba286..e34b304c 100644
--- a/sPuReMD/src/forces.h
+++ b/sPuReMD/src/forces.h
@@ -25,13 +25,14 @@
#include "reax_types.h"
-void Init_Bonded_Force_Functions( control_params* );
+void Init_Bonded_Force_Functions( control_params * const );
-void Compute_Forces( reax_system *, control_params *, simulation_data *,
- static_storage *, reax_list **, output_controls *, int );
+void Estimate_Storages( reax_system const * const, control_params const * const,
+ static_storage * const, reax_list ** const, int, int );
-void Estimate_Storage_Sizes( reax_system *, control_params *, reax_list **,
- int *, int *, int *, int * );
+int Compute_Forces( reax_system * const, control_params * const,
+ simulation_data * const, static_storage * const, reax_list ** const,
+ output_controls * const, int );
#endif
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 7b43cb96..cd9d6075 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -38,8 +38,8 @@
#include "vector.h"
-static void Generate_Initial_Velocities( reax_system *system,
- control_params *control, real T )
+static void Generate_Initial_Velocities( reax_system * const system,
+ control_params const * const control, real T )
{
int i;
real scale, norm;
@@ -98,8 +98,8 @@ static void Generate_Initial_Velocities( reax_system *system,
}
-static void Init_System( reax_system *system, control_params *control,
- simulation_data *data )
+static void Init_System( reax_system * const system, control_params * const control,
+ simulation_data * const data, static_storage * const workspace, int realloc )
{
int i;
rvec dx;
@@ -155,12 +155,26 @@ static void Init_System( reax_system *system, control_params *control,
}
Setup_Grid( system );
+
+ Bin_Atoms( system, workspace );
+
+#if defined(REORDER_ATOMS)
+ Reorder_Atoms( system, workspace, control );
+#endif
+
+ if ( realloc == TRUE )
+ {
+ /* list management */
+ system->bonds = smalloc( sizeof(int) * system->N_max, __FILE__, __LINE__ );
+
+ system->hbonds = smalloc( sizeof(int) * system->N_max, __FILE__, __LINE__ );
+ }
}
-static void Init_Simulation_Data( reax_system *system, control_params *control,
- simulation_data *data, output_controls *out_control,
- evolve_function *Evolve, int realloc )
+static void Init_Simulation_Data( reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ evolve_function * const Evolve, int realloc )
{
#if defined(_OPENMP)
if ( realloc == TRUE && (control->ensemble == sNPT || control->ensemble == iNPT
@@ -278,8 +292,9 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
}
-/* Initialize Taper params */
-static void Init_Taper( control_params *control, static_storage *workspace )
+/* Initialize taper function applied to van der Waals and Coulomb interactions */
+static void Init_Taper( control_params const * const control,
+ static_storage * const workspace )
{
real d1, d7;
real swa, swa2, swa3;
@@ -322,8 +337,9 @@ static void Init_Taper( control_params *control, static_storage *workspace )
}
-static void Init_Workspace( reax_system *system, control_params *control,
- static_storage *workspace, int realloc )
+static void Init_Workspace( reax_system * const system,
+ control_params const * const control, static_storage * const workspace,
+ int realloc )
{
int i;
@@ -770,71 +786,82 @@ static void Init_Workspace( reax_system *system, control_params *control,
__FILE__, __LINE__ );
#endif
- workspace->realloc.num_far = -1;
- workspace->realloc.Htop = -1;
- workspace->realloc.hbonds = -1;
- workspace->realloc.bonds = -1;
- workspace->realloc.num_3body = -1;
+ workspace->realloc.far_nbrs = FALSE;
+ workspace->realloc.cm = FALSE;
+ workspace->realloc.hbonds = FALSE;
+ workspace->realloc.bonds = FALSE;
+ workspace->realloc.thbody = FALSE;
workspace->realloc.gcell_atoms = -1;
Reset_Workspace( system, workspace );
- /* Initialize Taper function */
Init_Taper( control, workspace );
}
-static void Init_Lists( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control, int realloc )
+static void Init_Lists( reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ static_storage * const workspace, reax_list ** const lists, int realloc )
{
- int i, num_nbrs, num_bonds, num_hbonds, num_3body, Htop;
- int *hb_top, *bond_top;
-
- num_nbrs = Estimate_Num_Neighbors( system, control, workspace, lists );
+ int i, ret;
- if ( lists[FAR_NBRS]->allocated == FALSE )
- {
- Make_List( system->N, system->N_max, num_nbrs, TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
- }
- else if ( realloc == TRUE || lists[FAR_NBRS]->total_intrs < num_nbrs )
+ if ( realloc == TRUE )
{
- if ( lists[FAR_NBRS]->allocated == TRUE )
+ Estimate_Num_Neighbors( system, control, workspace, lists );
+
+ if ( lists[FAR_NBRS]->allocated == FALSE )
{
- Delete_List( TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
+ Make_List( system->N, system->N_max, workspace->realloc.total_far_nbrs,
+ TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
+ }
+ else if ( realloc == TRUE
+ || lists[FAR_NBRS]->total_intrs < workspace->realloc.total_far_nbrs )
+ {
+ if ( lists[FAR_NBRS]->allocated == TRUE )
+ {
+ Delete_List( TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
+ }
+ Make_List( system->N, system->N_max,
+ MAX( workspace->realloc.total_far_nbrs, lists[FAR_NBRS]->total_intrs ),
+ TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
}
- Make_List( system->N, system->N_max,
- MAX( num_nbrs, lists[FAR_NBRS]->total_intrs ),
- TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
}
- else
+
+ lists[FAR_NBRS]->n = system->N;
+
+ ret = Generate_Neighbor_Lists( system, control, data, workspace, lists );
+
+ if ( ret != SUCCESS )
{
- lists[FAR_NBRS]->n = system->N;
- }
+ Estimate_Num_Neighbors( system, control, workspace, lists );
- Generate_Neighbor_Lists( system, control, data, workspace, lists );
+ Delete_List( TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
+ Make_List( system->N, system->N_max,
+ MAX( workspace->realloc.total_far_nbrs, lists[FAR_NBRS]->total_intrs ),
+ TYP_FAR_NEIGHBOR, lists[FAR_NBRS] );
- Htop = 0;
- hb_top = scalloc( system->N, sizeof(int), __FILE__, __LINE__ );
- bond_top = scalloc( system->N, sizeof(int), __FILE__, __LINE__ );
- num_3body = 0;
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
+ }
- Estimate_Storage_Sizes( system, control, lists, &Htop,
- hb_top, bond_top, &num_3body );
- num_3body = MAX( num_3body, MIN_BONDS );
+ if ( realloc == TRUE )
+ {
+ Estimate_Storages( system, control, workspace, lists, TRUE, TRUE );
+ }
if ( workspace->H.allocated == FALSE )
{
- Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max, Htop );
+ Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max,
+ workspace->realloc.total_cm_entries );
}
- else if ( realloc == TRUE || workspace->H.m < Htop
+ else if ( realloc == TRUE || workspace->H.m < workspace->realloc.total_cm_entries
|| workspace->H.n_max < system->N_cm_max )
{
if ( workspace->H.allocated == TRUE )
{
Deallocate_Matrix( &workspace->H );
}
- Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max, Htop );
+ Allocate_Matrix( &workspace->H, system->N_cm, system->N_cm_max,
+ workspace->realloc.total_cm_entries );
}
else
{
@@ -844,12 +871,13 @@ static void Init_Lists( reax_system *system, control_params *control,
if ( workspace->H_sp.allocated == FALSE )
{
/* TODO: better estimate for H_sp?
- * If so, need to refactor Estimate_Storage_Sizes
+ * If so, need to refactor Estimate_Storages
* to use various cut-off distances as parameters
* (non-bonded, hydrogen, 3body, etc.) */
- Allocate_Matrix( &workspace->H_sp, system->N_cm, system->N_cm_max, Htop );
+ Allocate_Matrix( &workspace->H_sp, system->N_cm, system->N_cm_max,
+ workspace->realloc.total_cm_entries );
}
- else if ( realloc == TRUE || workspace->H_sp.m < Htop
+ else if ( realloc == TRUE || workspace->H_sp.m < workspace->realloc.total_cm_entries
|| workspace->H.n_max < system->N_cm_max )
{
if ( workspace->H_sp.allocated == TRUE )
@@ -857,10 +885,11 @@ static void Init_Lists( reax_system *system, control_params *control,
Deallocate_Matrix( &workspace->H_sp );
}
/* TODO: better estimate for H_sp?
- * If so, need to refactor Estimate_Storage_Sizes
+ * If so, need to refactor Estimate_Storages
* to use various cut-off distances as parameters
* (non-bonded, hydrogen, 3body, etc.) */
- Allocate_Matrix( &workspace->H_sp, system->N_cm, system->N_cm_max, Htop );
+ Allocate_Matrix( &workspace->H_sp, system->N_cm, system->N_cm_max,
+ workspace->realloc.total_cm_entries );
}
else
{
@@ -870,7 +899,6 @@ static void Init_Lists( reax_system *system, control_params *control,
workspace->num_H = 0;
if ( control->hbond_cut > 0.0 )
{
- /* init hydrogen atom indexes */
for ( i = 0; i < system->N; ++i )
{
if ( system->reax_param.sbp[ system->atoms[i].type ].p_hbond == H_ATOM )
@@ -882,78 +910,60 @@ static void Init_Lists( reax_system *system, control_params *control,
workspace->hbond_index[i] = -1;
}
}
-
- if ( workspace->num_H == 0 )
+ for ( i = system->N; i < system->N_max; ++i )
{
- control->hbond_cut = 0.0;
+ workspace->hbond_index[i] = -1;
}
- else
+ }
+
+ if ( control->hbond_cut > 0.0 && workspace->num_H > 0 )
+ {
+ if ( lists[HBONDS]->allocated == FALSE )
{
- num_hbonds = 0;
- for ( i = 0; i < system->N; ++i )
- {
- num_hbonds += hb_top[i];
- }
+ workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
- if ( lists[HBONDS]->allocated == FALSE )
+ Make_List( workspace->num_H, workspace->num_H_max,
+ (int) CEIL( SAFE_ZONE * workspace->realloc.total_hbonds ),
+ TYP_HBOND, lists[HBONDS] );
+ }
+ else if ( workspace->num_H_max < workspace->num_H
+ || lists[HBONDS]->total_intrs < workspace->realloc.total_hbonds )
+ {
+ if ( workspace->num_H_max < workspace->num_H )
{
workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
-
- Make_List( workspace->num_H, workspace->num_H_max,
- (int) CEIL( SAFE_ZONE * num_hbonds ),
- TYP_HBOND, lists[HBONDS] );
}
- else if ( workspace->num_H_max < workspace->num_H
- || lists[HBONDS]->total_intrs < num_hbonds )
- {
- if ( workspace->num_H_max < workspace->num_H )
- {
- workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
- }
- if ( lists[HBONDS]->allocated == TRUE )
- {
- Delete_List( TYP_HBOND, lists[HBONDS] );
- }
- Make_List( workspace->num_H, workspace->num_H_max,
- MAX( num_hbonds, lists[HBONDS]->total_intrs ),
- TYP_HBOND, lists[HBONDS] );
- }
- else
+ if ( lists[HBONDS]->allocated == TRUE )
{
- lists[HBONDS]->n = workspace->num_H;
+ Delete_List( TYP_HBOND, lists[HBONDS] );
}
-
- Initialize_HBond_List( system->N, workspace->hbond_index, hb_top, lists[HBONDS] );
+ Make_List( workspace->num_H, workspace->num_H_max,
+ MAX( workspace->realloc.total_hbonds, lists[HBONDS]->total_intrs ),
+ TYP_HBOND, lists[HBONDS] );
+ }
+ else
+ {
+ lists[HBONDS]->n = workspace->num_H;
}
- }
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "estimated storage - num_hbonds: %d\n", num_hbonds );
- fprintf( stderr, "memory allocated: hbonds = %ldMB\n",
- num_hbonds * sizeof(hbond_data) / (1024 * 1024) );
-#endif
- num_bonds = 0;
- for ( i = 0; i < system->N; ++i )
- {
- num_bonds += bond_top[i];
+ Init_List_Indices( lists[HBONDS] );
}
/* bonds list */
if ( lists[BONDS]->allocated == FALSE )
{
- Make_List( system->N, system->N_max, (int) CEIL( num_bonds * SAFE_ZONE ),
+ Make_List( system->N, system->N_max, (int) CEIL( workspace->realloc.total_bonds * SAFE_ZONE ),
TYP_BOND, lists[BONDS] );
}
- else if ( realloc == TRUE || lists[BONDS]->total_intrs < num_bonds )
+ else if ( realloc == TRUE || lists[BONDS]->total_intrs < workspace->realloc.total_bonds )
{
if ( lists[BONDS]->allocated == TRUE )
{
Delete_List( TYP_BOND, lists[BONDS] );
}
Make_List( system->N, system->N_max,
- MAX( num_bonds, lists[BONDS]->total_intrs ),
+ MAX( workspace->realloc.total_bonds, lists[BONDS]->total_intrs ),
TYP_BOND, lists[BONDS] );
}
else
@@ -961,45 +971,34 @@ static void Init_Lists( reax_system *system, control_params *control,
lists[BONDS]->n = system->N;
}
- Initialize_Bond_List( bond_top, lists[BONDS] );
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "estimated storage - num_bonds: %d\n", num_bonds );
- fprintf( stderr, "memory allocated: bonds = %ldMB\n",
- num_bonds * sizeof(bond_data) / (1024 * 1024) );
-#endif
+ Init_List_Indices( lists[BONDS] );
/* 3bodies list */
if ( lists[THREE_BODIES]->allocated == FALSE )
{
- Make_List( num_bonds, num_bonds, num_3body, TYP_THREE_BODY, lists[THREE_BODIES] );
+ Make_List( workspace->realloc.total_bonds, (int) CEIL( workspace->realloc.total_bonds * SAFE_ZONE ),
+ workspace->realloc.total_thbodies, TYP_THREE_BODY, lists[THREE_BODIES] );
}
- else if ( lists[THREE_BODIES]->n_max < num_bonds
- || lists[THREE_BODIES]->total_intrs < num_3body )
+ else if ( lists[THREE_BODIES]->n_max < workspace->realloc.total_bonds
+ || lists[THREE_BODIES]->total_intrs < workspace->realloc.total_thbodies )
{
if ( lists[THREE_BODIES]->allocated == TRUE )
{
Delete_List( TYP_THREE_BODY, lists[THREE_BODIES] );
}
- Make_List( MAX( num_bonds, lists[THREE_BODIES]->n_max),
- MAX( num_bonds, lists[THREE_BODIES]->n_max),
- MAX( num_3body, lists[THREE_BODIES]->total_intrs ),
+ Make_List( workspace->realloc.total_bonds,
+ MAX( workspace->realloc.total_bonds, lists[THREE_BODIES]->n_max),
+ MAX( workspace->realloc.total_thbodies, lists[THREE_BODIES]->total_intrs ),
TYP_THREE_BODY, lists[THREE_BODIES] );
}
else
{
- lists[THREE_BODIES]->n = num_bonds;
+ lists[THREE_BODIES]->n = workspace->realloc.total_bonds;
}
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "estimated storage - num_3body: %d\n", num_3body );
- fprintf( stderr, "memory allocated: 3-body = %ldMB\n",
- num_3body * sizeof(three_body_interaction_data) / (1024 * 1024) );
-#endif
-
#if defined(TEST_FORCES)
//TODO: increased num. of DDELTA list elements, find a better count later
- Make_List( system->N, num_bonds * 20, TYP_DDELTA, lists[DDELTA] );
+ Make_List( system->N, workspace->realloc.total_bonds * 20, TYP_DDELTA, lists[DDELTA] );
for ( i = 0; i < lists[DDELTA]->n; ++i )
{
@@ -1007,7 +1006,8 @@ static void Init_Lists( reax_system *system, control_params *control,
Set_End_Index( i, 0, lists[DDELTA] );
}
- Make_List( num_bonds, num_bonds * MAX_BONDS * 3, TYP_DBO, lists[DBO] );
+ Make_List( workspace->realloc.total_bonds, workspace->realloc.total_bonds * MAX_BONDS * 3,
+ TYP_DBO, lists[DBO] );
for ( i = 0; i < lists[DBO]->n; ++i )
{
@@ -1015,9 +1015,6 @@ static void Init_Lists( reax_system *system, control_params *control,
Set_End_Index( i, 0, lists[DBO] );
}
#endif
-
- sfree( hb_top, __FILE__, __LINE__ );
- sfree( bond_top, __FILE__, __LINE__ );
}
@@ -1282,10 +1279,10 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
}
-void Initialize( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace, reax_list **lists,
- output_controls *out_control, evolve_function *Evolve,
- int output_enabled, int realloc )
+void Initialize( reax_system * const system, control_params * const control,
+ simulation_data * const data, static_storage * const workspace,
+ reax_list ** const lists, output_controls * const out_control,
+ evolve_function * const Evolve, int output_enabled, int realloc )
{
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(control)
@@ -1308,13 +1305,13 @@ void Initialize( reax_system *system, control_params *control,
Randomize( );
- Init_System( system, control, data );
+ Init_System( system, control, data, workspace, realloc );
- Init_Simulation_Data( system, control, data, out_control, Evolve, realloc );
+ Init_Simulation_Data( system, control, data, Evolve, realloc );
Init_Workspace( system, control, workspace, realloc );
- Init_Lists( system, control, data, workspace, lists, out_control, realloc );
+ Init_Lists( system, control, data, workspace, lists, realloc );
Init_Out_Controls( system, control, workspace, out_control, output_enabled );
@@ -1386,6 +1383,9 @@ static void Finalize_System( reax_system *system, control_params *control,
sfree( system->atoms, __FILE__, __LINE__ );
}
+
+ sfree( system->bonds, __FILE__, __LINE__ );
+ sfree( system->hbonds, __FILE__, __LINE__ );
}
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index fc1a95f1..05893a8e 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -32,152 +32,234 @@
#include "vector.h"
-/* Perform simulation using the microcanonical ensemble with
- * the Velocity Verlet integrator. */
-void Velocity_Verlet_NVE( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+/* Velocity Verlet integrator for microcanonical ensemble. */
+int Velocity_Verlet_NVE( reax_system * const system, control_params * const control,
+ simulation_data * const data, static_storage * const workspace,
+ reax_list ** const lists, output_controls * const out_control )
{
- int i, renbr;
+ int i, renbr, ret;
+ static int verlet_part1_done = FALSE, gen_nbr_list = FALSE;
real inv_m, scalar1, scalar2;
rvec dx;
+ ret = SUCCESS;
renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
scalar1 = -0.5 * control->dt * F_CONV;
scalar2 = -0.5 * SQR( control->dt ) * F_CONV;
- for ( i = 0; i < system->N; i++ )
+ if ( verlet_part1_done == FALSE )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
- /* Compute x(t + dt) */
- rvec_ScaledSum( dx, control->dt, system->atoms[i].v,
- scalar2 * inv_m, system->atoms[i].f );
+ /* Compute x(t + dt) */
+ rvec_ScaledSum( dx, control->dt, system->atoms[i].v,
+ scalar2 * inv_m, system->atoms[i].f );
- control->update_atom_position( system->atoms[i].x, dx,
- system->atoms[i].rel_map, &system->box );
+ control->update_atom_position( system->atoms[i].x, dx,
+ system->atoms[i].rel_map, &system->box );
- /* Compute v(t + dt/2) */
- rvec_ScaledAdd( system->atoms[i].v,
- scalar1 * inv_m, system->atoms[i].f );
+ /* Compute v(t + dt/2) */
+ rvec_ScaledAdd( system->atoms[i].v,
+ scalar1 * inv_m, system->atoms[i].f );
+ }
+
+ if ( renbr == TRUE )
+ {
+ Reallocate_Part1( system, control, workspace, lists );
+
+ Bin_Atoms( system, workspace );
+
+#if defined(REORDER_ATOMS)
+ Reorder_Atoms( system, workspace, control );
+#endif
+ }
+
+ verlet_part1_done = TRUE;
}
- Reallocate( system, control, workspace, lists, renbr );
+ Reallocate_Part2( system, control, data, workspace, lists );
Reset( system, control, data, workspace, lists );
- if ( renbr == TRUE )
+ if ( renbr == TRUE && gen_nbr_list == FALSE )
{
- Generate_Neighbor_Lists( system, control, data, workspace, lists );
+ ret = Generate_Neighbor_Lists( system, control, data, workspace, lists );
+
+ if ( ret == SUCCESS )
+ {
+ gen_nbr_list = TRUE;
+ }
+ else
+ {
+ Estimate_Num_Neighbors( system, control, workspace, lists );
+ }
}
- Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
+ if ( ret == SUCCESS )
+ {
+ ret = Compute_Forces( system, control, data, workspace, lists,
+ out_control, FALSE );
+ }
- for ( i = 0; i < system->N; i++ )
+ if ( ret == SUCCESS )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
- /* Compute v(t + dt) */
- rvec_ScaledAdd( system->atoms[i].v,
- scalar1 * inv_m, system->atoms[i].f );
+ /* Compute v(t + dt) */
+ rvec_ScaledAdd( system->atoms[i].v,
+ scalar1 * inv_m, system->atoms[i].f );
+ }
+
+ verlet_part1_done = FALSE;
+ gen_nbr_list = FALSE;
}
+
+ return ret;
}
-/* Perform simulation using constant volume and temperature (NVT ensemble)
- * with a Berendsen thermostat and the Velocity Verlet integrator. */
-void Velocity_Verlet_Berendsen_NVT( reax_system* system,
- control_params* control, simulation_data *data,
- static_storage *workspace, reax_list **lists,
- output_controls *out_control )
+/* Velocity Verlet integrator for constant volume and temperature
+ * with Berendsen thermostat.
+ *
+ * NOTE: All box dimensions are scaled by the same amount, and
+ * there is no change in the angles between axes. */
+int Velocity_Verlet_Berendsen_NVT( reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ static_storage * const workspace, reax_list ** const lists,
+ output_controls * const out_control )
{
- int i, renbr;
+ int i, renbr, ret;
+ static int verlet_part1_done = FALSE, gen_nbr_list = FALSE;
real inv_m, scalar1, scalar2, lambda;
rvec dx;
+ ret = SUCCESS;
renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
scalar1 = -0.5 * control->dt * F_CONV;
scalar2 = -0.5 * SQR( control->dt ) * F_CONV;
- /* velocity verlet, 1st part */
- for ( i = 0; i < system->N; i++ )
+ if ( verlet_part1_done == FALSE )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ /* velocity verlet, 1st part */
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
- /* Compute x(t + dt) */
- rvec_ScaledSum( dx, control->dt, system->atoms[i].v, scalar2 * inv_m, system->atoms[i].f );
+ /* Compute x(t + dt) */
+ rvec_ScaledSum( dx, control->dt, system->atoms[i].v, scalar2 * inv_m, system->atoms[i].f );
- control->update_atom_position( system->atoms[i].x, dx, system->atoms[i].rel_map, &system->box );
+ control->update_atom_position( system->atoms[i].x, dx, system->atoms[i].rel_map, &system->box );
- /* Compute v(t + dt/2) */
- rvec_ScaledAdd( system->atoms[i].v, scalar1 * inv_m, system->atoms[i].f );
+ /* Compute v(t + dt/2) */
+ rvec_ScaledAdd( system->atoms[i].v, scalar1 * inv_m, system->atoms[i].f );
+ }
+
+ Reallocate_Part1( system, control, workspace, lists );
+
+ if ( renbr == TRUE )
+ {
+ Bin_Atoms( system, workspace );
+
+#if defined(REORDER_ATOMS)
+ Reorder_Atoms( system, workspace, control );
+#endif
+ }
+
+ verlet_part1_done = TRUE;
}
- Reallocate( system, control, workspace, lists, renbr );
+ Reallocate_Part2( system, control, data, workspace, lists );
Reset( system, control, data, workspace, lists );
- if ( renbr == TRUE )
+ if ( renbr == TRUE && gen_nbr_list == FALSE )
{
- Generate_Neighbor_Lists( system, control, data, workspace, lists );
+ ret = Generate_Neighbor_Lists( system, control, data, workspace, lists );
+
+ if ( ret == SUCCESS )
+ {
+ gen_nbr_list = TRUE;
+ }
+ else
+ {
+ Estimate_Num_Neighbors( system, control, workspace, lists );
+ }
}
- Compute_Forces( system, control, data, workspace,
- lists, out_control, FALSE );
+ if ( ret == SUCCESS )
+ {
+ ret = Compute_Forces( system, control, data, workspace, lists,
+ out_control, FALSE );
+ }
- /* velocity verlet, 2nd part */
- for ( i = 0; i < system->N; i++ )
+ if ( ret == SUCCESS )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ /* velocity verlet, 2nd part */
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
- /* Compute v(t + dt) */
- rvec_ScaledAdd( system->atoms[i].v, scalar1 * inv_m, system->atoms[i].f );
- }
+ /* Compute v(t + dt) */
+ rvec_ScaledAdd( system->atoms[i].v, scalar1 * inv_m, system->atoms[i].f );
+ }
- Compute_Kinetic_Energy( system, data );
+ Compute_Kinetic_Energy( system, data );
- /* temperature scaler */
- lambda = 1.0 + ((control->dt * 1.0e-12) / control->Tau_T)
- * (control->T / data->therm.T - 1.0);
+ /* temperature scaler */
+ lambda = 1.0 + ((control->dt * 1.0e-12) / control->Tau_T)
+ * (control->T / data->therm.T - 1.0);
- if ( lambda < MIN_dT )
- {
- lambda = MIN_dT;
- }
+ if ( lambda < MIN_dT )
+ {
+ lambda = MIN_dT;
+ }
- lambda = SQRT( lambda );
+ lambda = SQRT( lambda );
- if ( lambda > MAX_dT )
- {
- lambda = MAX_dT;
- }
+ if ( lambda > MAX_dT )
+ {
+ lambda = MAX_dT;
+ }
- /* Scale velocities and positions at t+dt */
- for ( i = 0; i < system->N; ++i )
- {
- rvec_Scale( system->atoms[i].v, lambda, system->atoms[i].v );
+ /* Scale velocities and positions at t+dt */
+ for ( i = 0; i < system->N; ++i )
+ {
+ rvec_Scale( system->atoms[i].v, lambda, system->atoms[i].v );
+ }
+
+ /* update kinetic energy and temperature based on new velocities */
+ Compute_Kinetic_Energy( system, data );
+
+ verlet_part1_done = FALSE;
+ gen_nbr_list = FALSE;
}
- /* update kinetic energy and temperature based on new velocities */
- Compute_Kinetic_Energy( system, data );
+ return ret;
}
-/* Perform simulation using constant volume and temperature (NVT ensemble)
- * with a Nose-Hoover thermostat and the Velocity Verlet integrator.
+/* Velocity Verlet integrator for constant volume and constant temperature
+ * with Nose-Hoover thermostat.
*
* Reference: Understanding Molecular Simulation, Frenkel and Smit
* Academic Press Inc. San Diego, 1996 p. 388-391 */
-void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system* system, control_params* control,
- simulation_data *data, static_storage *workspace, reax_list **lists,
- output_controls *out_control )
+int Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ static_storage * const workspace, reax_list ** const lists,
+ output_controls * const out_control )
{
- int i, itr, renbr;
+ int i, itr, renbr, ret;
+ static int verlet_part1_done = FALSE, gen_nbr_list = FALSE;
real inv_m, scalar1, scalar2, scalar3, scalar4, coef_v;
real E_kin_new, G_xi_new, v_xi_new, v_xi_old;
rvec dx;
thermostat *therm;
+ ret = SUCCESS;
renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
scalar1 = -0.5 * control->dt * F_CONV;
scalar2 = -0.5 * SQR( control->dt ) * F_CONV;
@@ -185,296 +267,413 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system* system, control_params*
scalar4 = -0.5 * SQR( control->dt );
therm = &data->therm;
- /* Compute x(t + dt) and copy old forces */
- for ( i = 0; i < system->N; i++ )
+ if ( verlet_part1_done == FALSE )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ /* Compute x(t + dt) and copy old forces */
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
- rvec_ScaledSum( dx, control->dt + scalar4 * therm->v_xi, system->atoms[i].v,
- scalar2 * inv_m, system->atoms[i].f );
+ rvec_ScaledSum( dx, control->dt + scalar4 * therm->v_xi, system->atoms[i].v,
+ scalar2 * inv_m, system->atoms[i].f );
- control->update_atom_position( system->atoms[i].x, dx,
- system->atoms[i].rel_map, &system->box );
+ control->update_atom_position( system->atoms[i].x, dx,
+ system->atoms[i].rel_map, &system->box );
- rvec_Copy( workspace->f_old[i], system->atoms[i].f );
- }
+ rvec_Copy( workspace->f_old[i], system->atoms[i].f );
+ }
- /* Compute xi(t + dt) */
- therm->xi += therm->v_xi * control->dt + 0.5 * SQR( control->dt ) * therm->G_xi;
+ /* Compute xi(t + dt) */
+ therm->xi += therm->v_xi * control->dt + 0.5 * SQR( control->dt ) * therm->G_xi;
- Reallocate( system, control, workspace, lists, renbr );
+ if ( renbr == TRUE )
+ {
+ Reallocate_Part1( system, control, workspace, lists );
- Reset( system, control, data, workspace, lists );
+ Bin_Atoms( system, workspace );
- if ( renbr == TRUE )
- {
- Generate_Neighbor_Lists( system, control, data, workspace, lists );
+#if defined(REORDER_ATOMS)
+ Reorder_Atoms( system, workspace, control );
+#endif
+ }
+
+ verlet_part1_done = TRUE;
}
- /* Calculate Forces at time (t + dt) */
- Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
+ Reallocate_Part2( system, control, data, workspace, lists );
- /* Compute iteration constants for each atom's velocity */
- for ( i = 0; i < system->N; ++i )
+ Reset( system, control, data, workspace, lists );
+
+ if ( renbr == TRUE && gen_nbr_list == FALSE )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ ret = Generate_Neighbor_Lists( system, control, data, workspace, lists );
- rvec_Scale( workspace->v_const[i],
- 1.0 + scalar3 * therm->v_xi, system->atoms[i].v );
- rvec_ScaledAdd( workspace->v_const[i],
- scalar1 * inv_m, workspace->f_old[i] );
- rvec_ScaledAdd( workspace->v_const[i],
- scalar1 * inv_m, system->atoms[i].f );
+ if ( ret == SUCCESS )
+ {
+ gen_nbr_list = TRUE;
+ }
+ else
+ {
+ Estimate_Num_Neighbors( system, control, workspace, lists );
+ }
}
- v_xi_new = therm->v_xi_old + 2.0 * control->dt * therm->G_xi;
- E_kin_new = 0.0;
- G_xi_new = 0.0;
- v_xi_old = 0.0;
- itr = 0;
- do
+ if ( ret == SUCCESS )
{
- itr++;
- v_xi_old = v_xi_new;
- coef_v = 1.0 / (1.0 + 0.5 * control->dt * v_xi_old);
- E_kin_new = 0.0;
+ /* Calculate Forces at time (t + dt) */
+ ret = Compute_Forces( system, control, data, workspace, lists,
+ out_control, FALSE );
+ }
+ if ( ret == SUCCESS )
+ {
+ /* Compute iteration constants for each atom's velocity */
for ( i = 0; i < system->N; ++i )
{
- rvec_Scale( system->atoms[i].v, coef_v, workspace->v_const[i] );
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+
+ rvec_Scale( workspace->v_const[i],
+ 1.0 + scalar3 * therm->v_xi, system->atoms[i].v );
+ rvec_ScaledAdd( workspace->v_const[i],
+ scalar1 * inv_m, workspace->f_old[i] );
+ rvec_ScaledAdd( workspace->v_const[i],
+ scalar1 * inv_m, system->atoms[i].f );
+ }
- E_kin_new += 0.5 * system->reax_param.sbp[system->atoms[i].type].mass
- * rvec_Dot( system->atoms[i].v, system->atoms[i].v );
+ v_xi_new = therm->v_xi_old + 2.0 * control->dt * therm->G_xi;
+ E_kin_new = 0.0;
+ G_xi_new = 0.0;
+ v_xi_old = 0.0;
+ itr = 0;
+ do
+ {
+ itr++;
+ v_xi_old = v_xi_new;
+ coef_v = 1.0 / (1.0 + 0.5 * control->dt * v_xi_old);
+ E_kin_new = 0.0;
+
+ for ( i = 0; i < system->N; ++i )
+ {
+ rvec_Scale( system->atoms[i].v, coef_v, workspace->v_const[i] );
+
+ E_kin_new += 0.5 * system->reax_param.sbp[system->atoms[i].type].mass
+ * rvec_Dot( system->atoms[i].v, system->atoms[i].v );
+ }
+
+ G_xi_new = control->Tau_T * (2.0 * E_kin_new
+ - data->N_f * K_B * control->T);
+ v_xi_new = therm->v_xi + 0.5 * control->dt * (therm->G_xi + G_xi_new);
}
+ while ( FABS( v_xi_new - v_xi_old ) > 1.0e-5 );
- G_xi_new = control->Tau_T * (2.0 * E_kin_new
- - data->N_f * K_B * control->T);
- v_xi_new = therm->v_xi + 0.5 * control->dt * (therm->G_xi + G_xi_new);
+ therm->v_xi_old = therm->v_xi;
+ therm->v_xi = v_xi_new;
+ therm->G_xi = G_xi_new;
+
+ verlet_part1_done = FALSE;
+ gen_nbr_list = FALSE;
}
- while ( FABS( v_xi_new - v_xi_old ) > 1.0e-5 );
- therm->v_xi_old = therm->v_xi;
- therm->v_xi = v_xi_new;
- therm->G_xi = G_xi_new;
+ return ret;
}
-/* Perform simulation using constant pressure and temperature (NPT ensemble)
- * with a Berendsen thermostat and the Velocity Verlet integrator.
+/* Velocity Verlet integrator for constant pressure and constant temperature
+ * with a Berendsen thermostat.
*
* NOTE: All box dimensions are scaled by the same amount, and
* there is no change in the angles between axes. */
-void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
- control_params* control, simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+int Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ static_storage * const workspace, reax_list ** const lists,
+ output_controls * const out_control )
{
- int i, renbr;
+ int i, renbr, ret;
+ static int verlet_part1_done = FALSE, gen_nbr_list = FALSE;
real inv_m, dt, lambda, mu;
rvec dx, mu_3;
+ ret = SUCCESS;
dt = control->dt;
renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
- /* velocity verlet, 1st part */
- for ( i = 0; i < system->N; i++ )
+ if ( verlet_part1_done == FALSE )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ /* velocity verlet, 1st part */
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
- /* Compute x(t + dt) */
- rvec_ScaledSum( dx, dt, system->atoms[i].v,
- F_CONV * inv_m * -0.5 * SQR(dt), system->atoms[i].f );
+ /* Compute x(t + dt) */
+ rvec_ScaledSum( dx, dt, system->atoms[i].v,
+ F_CONV * inv_m * -0.5 * SQR(dt), system->atoms[i].f );
- control->update_atom_position( system->atoms[i].x, dx, system->atoms[i].rel_map, &system->box );
+ control->update_atom_position( system->atoms[i].x, dx, system->atoms[i].rel_map, &system->box );
- /* Compute v(t + dt/2) */
- rvec_ScaledAdd( system->atoms[i].v,
- F_CONV * inv_m * -0.5 * dt, system->atoms[i].f );
+ /* Compute v(t + dt/2) */
+ rvec_ScaledAdd( system->atoms[i].v,
+ F_CONV * inv_m * -0.5 * dt, system->atoms[i].f );
+ }
+
+ if ( renbr == TRUE )
+ {
+ Update_Grid( system );
+
+ Reallocate_Part1( system, control, workspace, lists );
+
+ Bin_Atoms( system, workspace );
+
+#if defined(REORDER_ATOMS)
+ Reorder_Atoms( system, workspace, control );
+#endif
+ }
+
+ verlet_part1_done = TRUE;
}
- Reallocate( system, control, workspace, lists, renbr );
+ Reallocate_Part2( system, control, data, workspace, lists );
Reset( system, control, data, workspace, lists );
- if ( renbr == TRUE )
+ if ( renbr == TRUE && gen_nbr_list == FALSE )
{
- Update_Grid( system );
+ ret = Generate_Neighbor_Lists( system, control, data, workspace, lists );
- Generate_Neighbor_Lists( system, control, data, workspace, lists );
+ if ( ret == SUCCESS )
+ {
+ gen_nbr_list = TRUE;
+ }
+ else
+ {
+ Estimate_Num_Neighbors( system, control, workspace, lists );
+ }
}
- Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
-
- /* velocity verlet, 2nd part */
- for ( i = 0; i < system->N; i++ )
+ if ( ret == SUCCESS )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
-
- /* Compute v(t + dt) */
- rvec_ScaledAdd( system->atoms[i].v,
- -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
+ ret = Compute_Forces( system, control, data, workspace, lists,
+ out_control, FALSE );
}
- Compute_Kinetic_Energy( system, data );
+ if ( ret == SUCCESS )
+ {
+ /* velocity verlet, 2nd part */
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+
+ /* Compute v(t + dt) */
+ rvec_ScaledAdd( system->atoms[i].v,
+ -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
+ }
- Compute_Pressure_Isotropic( system, control, data, out_control );
+ Compute_Kinetic_Energy( system, data );
- /* pressure scaler */
- for ( i = 0; i < 3; ++i )
- {
- mu_3[i] = POW( 1.0 + (dt / control->Tau_P[i])
- * (data->tot_press[i] - control->P[i]), 1.0 / 3.0 );
+ Compute_Pressure_Isotropic( system, control, data, out_control );
- if ( mu_3[i] < MIN_dV )
+ /* pressure scaler */
+ for ( i = 0; i < 3; ++i )
{
- mu_3[i] = MIN_dV;
+ mu_3[i] = POW( 1.0 + (dt / control->Tau_P[i])
+ * (data->tot_press[i] - control->P[i]), 1.0 / 3.0 );
+
+ if ( mu_3[i] < MIN_dV )
+ {
+ mu_3[i] = MIN_dV;
+ }
+ else if ( mu_3[i] > MAX_dV )
+ {
+ mu_3[i] = MAX_dV;
+ }
}
- else if ( mu_3[i] > MAX_dV )
+ mu = (mu_3[0] + mu_3[1] + mu_3[2]) / 3.0;
+
+ /* temperature scaler */
+ lambda = 1.0 + ((dt * 1.0e-12) / control->Tau_T)
+ * (control->T / data->therm.T - 1.0);
+
+ if ( lambda < MIN_dT )
{
- mu_3[i] = MAX_dV;
+ lambda = MIN_dT;
}
- }
- mu = (mu_3[0] + mu_3[1] + mu_3[2]) / 3.0;
- /* temperature scaler */
- lambda = 1.0 + ((dt * 1.0e-12) / control->Tau_T)
- * (control->T / data->therm.T - 1.0);
+ lambda = SQRT( lambda );
- if ( lambda < MIN_dT )
- {
- lambda = MIN_dT;
- }
+ if ( lambda > MAX_dT )
+ {
+ lambda = MAX_dT;
+ }
- lambda = SQRT( lambda );
+ /* Scale velocities and positions at t+dt */
+ for ( i = 0; i < system->N; ++i )
+ {
+ rvec_Scale( system->atoms[i].v, lambda, system->atoms[i].v );
+
+ /* IMPORTANT: What Adri does with scaling positions first to
+ * unit coordinates and then back to cartesian coordinates essentially
+ * is scaling the coordinates with mu^2. However, this causes unphysical
+ * modifications on the system because box dimensions
+ * are being scaled with mu! We need to discuss this with Adri! */
+ rvec_Scale( system->atoms[i].x, mu, system->atoms[i].x );
+ }
- if ( lambda > MAX_dT )
- {
- lambda = MAX_dT;
- }
+ /* update kinetic energy and temperature based on new velocities */
+ Compute_Kinetic_Energy( system, data );
- /* Scale velocities and positions at t+dt */
- for ( i = 0; i < system->N; ++i )
- {
- rvec_Scale( system->atoms[i].v, lambda, system->atoms[i].v );
-
- /* IMPORTANT: What Adri does with scaling positions first to
- * unit coordinates and then back to cartesian coordinates essentially
- * is scaling the coordinates with mu^2. However, this causes unphysical
- * modifications on the system because box dimensions
- * are being scaled with mu! We need to discuss this with Adri! */
- rvec_Scale( system->atoms[i].x, mu, system->atoms[i].x );
- }
+ Update_Box_Isotropic( &system->box, mu );
- /* update kinetic energy and temperature based on new velocities */
- Compute_Kinetic_Energy( system, data );
+ verlet_part1_done = FALSE;
+ gen_nbr_list = FALSE;
+ }
- Update_Box_Isotropic( &system->box, mu );
+ return ret;
}
/* Perform simulation using constant pressure and temperature (NPT ensemble)
* with a Berendsen thermostat and the Velocity Verlet integrator. */
-void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
- control_params* control, simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+int Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ static_storage * const workspace, reax_list ** const lists,
+ output_controls * const out_control )
{
- int i, renbr;
+ int i, renbr, ret;
+ static int verlet_part1_done = FALSE, gen_nbr_list = FALSE;
real dt, inv_m, lambda;
rvec dx, mu;
+ ret = SUCCESS;
dt = control->dt;
renbr = ((data->step - data->prev_steps) % control->reneighbor) == 0 ? TRUE : FALSE;
- /* velocity verlet, 1st part */
- for ( i = 0; i < system->N; i++ )
+ if ( verlet_part1_done == FALSE )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ /* velocity verlet, 1st part */
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
- /* Compute x(t + dt) */
- rvec_ScaledSum( dx, dt, system->atoms[i].v,
- -0.5 * F_CONV * inv_m * SQR(dt), system->atoms[i].f );
+ /* Compute x(t + dt) */
+ rvec_ScaledSum( dx, dt, system->atoms[i].v,
+ -0.5 * F_CONV * inv_m * SQR(dt), system->atoms[i].f );
- control->update_atom_position( system->atoms[i].x, dx, system->atoms[i].rel_map, &system->box );
+ control->update_atom_position( system->atoms[i].x, dx, system->atoms[i].rel_map, &system->box );
- /* Compute v(t + dt/2) */
- rvec_ScaledAdd( system->atoms[i].v,
- -0.5 * F_CONV * inv_m * dt, system->atoms[i].f );
+ /* Compute v(t + dt/2) */
+ rvec_ScaledAdd( system->atoms[i].v,
+ -0.5 * F_CONV * inv_m * dt, system->atoms[i].f );
+ }
+
+ if ( renbr == TRUE )
+ {
+ Update_Grid( system );
+
+ Reallocate_Part1( system, control, workspace, lists );
+
+ Bin_Atoms( system, workspace );
+
+#if defined(REORDER_ATOMS)
+ Reorder_Atoms( system, workspace, control );
+#endif
+ }
+
+ verlet_part1_done = TRUE;
}
- Reallocate( system, control, workspace, lists, renbr );
+ Reallocate_Part2( system, control, data, workspace, lists );
Reset( system, control, data, workspace, lists );
- if ( renbr == TRUE )
+ if ( renbr == TRUE && gen_nbr_list == FALSE )
{
- Update_Grid( system );
+ ret = Generate_Neighbor_Lists( system, control, data, workspace, lists );
- Generate_Neighbor_Lists( system, control, data, workspace, lists );
+ if ( ret == SUCCESS )
+ {
+ gen_nbr_list = TRUE;
+ }
+ else
+ {
+ Estimate_Num_Neighbors( system, control, workspace, lists );
+ }
}
- Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
-
- /* velocity verlet, 2nd part */
- for ( i = 0; i < system->N; i++ )
+ if ( ret == SUCCESS )
{
- inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
- /* Compute v(t + dt) */
- rvec_ScaledAdd( system->atoms[i].v,
- -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
+ ret = Compute_Forces( system, control, data, workspace, lists,
+ out_control, FALSE );
}
- Compute_Kinetic_Energy( system, data );
+ if ( ret == SUCCESS )
+ {
+ /* velocity verlet, 2nd part */
+ for ( i = 0; i < system->N; i++ )
+ {
+ inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+ /* Compute v(t + dt) */
+ rvec_ScaledAdd( system->atoms[i].v,
+ -0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
+ }
- Compute_Pressure_Isotropic( system, control, data, out_control );
+ Compute_Kinetic_Energy( system, data );
- /* pressure scaler */
- for ( i = 0; i < 3; ++i )
- {
- mu[i] = POW( 1.0 + (dt / control->Tau_P[i])
- * (data->tot_press[i] - control->P[i]), 1.0 / 3.0 );
+ Compute_Pressure_Isotropic( system, control, data, out_control );
- if ( mu[i] < MIN_dV )
+ /* pressure scaler */
+ for ( i = 0; i < 3; ++i )
{
- mu[i] = MIN_dV;
+ mu[i] = POW( 1.0 + (dt / control->Tau_P[i])
+ * (data->tot_press[i] - control->P[i]), 1.0 / 3.0 );
+
+ if ( mu[i] < MIN_dV )
+ {
+ mu[i] = MIN_dV;
+ }
+ else if ( mu[i] > MAX_dV )
+ {
+ mu[i] = MAX_dV;
+ }
}
- else if ( mu[i] > MAX_dV )
+
+ /* temperature scaler */
+ lambda = 1.0 + ((dt * 1.0e-12) / control->Tau_T)
+ * (control->T / data->therm.T - 1.0);
+
+ if ( lambda < MIN_dT )
{
- mu[i] = MAX_dV;
+ lambda = MIN_dT;
}
- }
- /* temperature scaler */
- lambda = 1.0 + ((dt * 1.0e-12) / control->Tau_T)
- * (control->T / data->therm.T - 1.0);
+ lambda = SQRT( lambda );
- if ( lambda < MIN_dT )
- {
- lambda = MIN_dT;
- }
+ if ( lambda > MAX_dT )
+ {
+ lambda = MAX_dT;
+ }
- lambda = SQRT( lambda );
+ /* Scale velocities and positions at t+dt */
+ for ( i = 0; i < system->N; ++i )
+ {
+ rvec_Scale( system->atoms[i].v, lambda, system->atoms[i].v );
+
+ /* IMPORTANT: What Adri does with scaling positions first to
+ * unit coordinates and then back to cartesian coordinates essentially
+ * is scaling the coordinates with mu^2. However, this causes unphysical
+ * modifications on the system because box dimensions
+ * are being scaled with mu! We need to discuss this with Adri! */
+ rvec_Multiply( system->atoms[i].x, mu, system->atoms[i].x );
+ }
- if ( lambda > MAX_dT )
- {
- lambda = MAX_dT;
- }
+ /* update kinetic energy and temperature based on new velocities */
+ Compute_Kinetic_Energy( system, data );
- /* Scale velocities and positions at t+dt */
- for ( i = 0; i < system->N; ++i )
- {
- rvec_Scale( system->atoms[i].v, lambda, system->atoms[i].v );
-
- /* IMPORTANT: What Adri does with scaling positions first to
- * unit coordinates and then back to cartesian coordinates essentially
- * is scaling the coordinates with mu^2. However, this causes unphysical
- * modifications on the system because box dimensions
- * are being scaled with mu! We need to discuss this with Adri! */
- rvec_Multiply( system->atoms[i].x, mu, system->atoms[i].x );
- }
+ Update_Box_Semi_Isotropic( &system->box, mu );
- /* update kinetic energy and temperature based on new velocities */
- Compute_Kinetic_Energy( system, data );
+ verlet_part1_done = FALSE;
+ gen_nbr_list = FALSE;
+ }
- Update_Box_Semi_Isotropic( &system->box, mu );
+ return ret;
}
@@ -482,16 +681,19 @@ void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
/* BELOW FUNCTIONS ARE NOT BEING USED ANYMORE! */
/************************************************/
#if defined(ANISOTROPIC)
-void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* control,
- simulation_data *data, static_storage *workspace, reax_list **lists,
- output_controls *out_control )
+int Velocity_Verlet_Nose_Hoover_NVT( reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ static_storage * const workspace, reax_list ** const lists,
+ output_controls * const out_control )
{
- int i;
+ int i, ret;
real inv_m;
real dt = control->dt;
real dt_sqr = SQR(dt);
rvec dx;
+ ret = SUCCESS;
+
for ( i = 0; i < system->N; i++ )
{
inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
@@ -515,6 +717,15 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* contr
data->therm.xi += 0.5 * dt * control->Tau_T
* (2.0 * data->E_Kin - data->N_f * K_B / F_CONV * control->T);
+ if ( renbr == TRUE )
+ {
+ Bin_Atoms( system, workspace );
+
+#if defined(REORDER_ATOMS)
+ Reorder_Atoms( system, workspace, control );
+#endif
+ }
+
Reset( system, control, data, workspace );
fprintf(out_control->log, "reset-");
fflush( out_control->log );
@@ -552,14 +763,17 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* contr
fprintf( out_control->log, "Xi: %8.3f %8.3f %8.3f\n",
data->therm.xi, data->E_Kin, data->N_f * K_B / F_CONV * control->T );
fflush( out_control->log );
+
+ return ret;
}
-void Velocity_Verlet_Isotropic_NPT( reax_system* system, control_params* control,
- simulation_data *data, static_storage *workspace, reax_list **lists,
- output_controls *out_control )
+int Velocity_Verlet_Isotropic_NPT( reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ static_storage * const workspace, reax_list ** const lists,
+ output_controls * const out_control )
{
- int i, itr;
+ int i, itr, ret;
real deps, v_eps_new = 0, v_eps_old = 0, G_xi_new;
real dxi, v_xi_new = 0, v_xi_old = 0, a_eps_new;
real inv_m, exp_deps, inv_3V;
@@ -572,6 +786,8 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system, control_params* control
simulation_box *box = &system->box;
rvec dx, dv;
+ ret = SUCCESS;
+
// Here we just calculate how much to increment eps, xi, v_eps, v_xi.
// Commits are done after positions and velocities of atoms are updated
// because position, velocity updates uses v_eps, v_xi terms;
@@ -618,6 +834,15 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system, control_params* control
iso_bar->eps += deps;
Update_Box_Isotropic( &system->box, EXP( 3.0 * iso_bar->eps ) );
+ if ( renbr == TRUE )
+ {
+ Bin_Atoms( system, workspace );
+
+#if defined(REORDER_ATOMS)
+ Reorder_Atoms( system, workspace, control );
+#endif
+ }
+
/* Calculate new forces, f(t + dt) */
Reset( system, control, data, workspace );
fprintf(out_control->log, "reset-");
@@ -718,5 +943,7 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system, control_params* control
iso_bar->a_eps, iso_bar->v_eps, iso_bar->eps);
fprintf(out_control->log, "xi: \tG- %8.3f v- %8.3f xi - %8.3f\n",
therm->G_xi, therm->v_xi, therm->xi);
+
+ return ret;
}
#endif
diff --git a/sPuReMD/src/integrate.h b/sPuReMD/src/integrate.h
index 7207fa8f..b7f6dd3b 100644
--- a/sPuReMD/src/integrate.h
+++ b/sPuReMD/src/integrate.h
@@ -26,26 +26,26 @@
#include "reax_types.h"
-void Velocity_Verlet_NVE( reax_system *, control_params *, simulation_data *,
+int Velocity_Verlet_NVE( reax_system *, control_params *, simulation_data *,
static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Berendsen_NVT( reax_system *, control_params *,
+int Velocity_Verlet_Berendsen_NVT( reax_system *, control_params *,
simulation_data *, static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system *, control_params *,
+int Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system *, control_params *,
simulation_data *, static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system *, control_params *,
+int Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system *, control_params *,
simulation_data *, static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system *, control_params *,
+int Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system *, control_params *,
simulation_data *, static_storage *, reax_list **, output_controls * );
#if defined(ANISOTROPIC)
-void Velocity_Verlet_Nose_Hoover_NVT( reax_system *, control_params *,
+int Velocity_Verlet_Nose_Hoover_NVT( reax_system *, control_params *,
simulation_data *, static_storage *, reax_list **, output_controls * );
-void Velocity_Verlet_Isotropic_NPT( reax_system *, control_params *,
+int Velocity_Verlet_Isotropic_NPT( reax_system *, control_params *,
simulation_data *, static_storage *, reax_list **, output_controls * );
#endif
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index fd9087e4..cd20c089 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -24,7 +24,8 @@
#include "tool_box.h"
-void Make_List( int n, int n_max, int total_intrs, int type, reax_list* l )
+/* Allocates a new instance of reax_list with internal lists of specified type */
+void Make_List( int n, int n_max, int total_intrs, int type, reax_list * const l )
{
assert( n > 0 );
assert( n_max > 0 );
@@ -141,7 +142,8 @@ void Make_List( int n, int n_max, int total_intrs, int type, reax_list* l )
}
-void Delete_List( int type, reax_list* l )
+/* Deallocates any space allocated for a reax_list instance */
+void Delete_List( int type, reax_list * const l )
{
assert( l != NULL );
@@ -217,3 +219,22 @@ void Delete_List( int type, reax_list* l )
break;
}
}
+
+
+/* Initialize list indices
+ *
+ * l: pointer to list
+ * */
+void Init_List_Indices( reax_list * const l )
+{
+ int i;
+
+ assert( l != NULL );
+ assert( l->n > 0 );
+
+ for ( i = 0; i < l->n; ++i )
+ {
+ Set_Start_Index( i, 0, l );
+ Set_End_Index( i, 0, l );
+ }
+}
diff --git a/sPuReMD/src/list.h b/sPuReMD/src/list.h
index 6c46c1e9..43b5dca6 100644
--- a/sPuReMD/src/list.h
+++ b/sPuReMD/src/list.h
@@ -25,11 +25,28 @@
#include "reax_types.h"
-void Make_List( int, int, int, int, reax_list* );
-
-void Delete_List( int, reax_list* );
-
-
+/* Allocates a new instance of reax_list with internal lists of specified type */
+void Make_List( int, int, int, int, reax_list * const );
+
+/* Deallocates any space allocated for a reax_list instance */
+void Delete_List( int, reax_list * const );
+
+void Init_List_Indices( reax_list * const );
+
+
+/* Return the size for the i-th internal list
+ *
+ * Inputs:
+ * l: list of lists (reax_list)
+ * i: list to return size for
+ *
+ * Returns:
+ * size of i-th internal list
+ *
+ * Unchecked runtime exceptions:
+ * l is not a NULL pointer
+ * i is a valid internal list index
+ * */
static inline int Num_Entries( int i, reax_list const * const l )
{
assert( l != NULL );
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index a673c46d..c75db365 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -42,10 +42,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 const * const, real );
typedef int (*find_far_neighbors_function)( rvec, rvec, int, int,
- simulation_box*, real, far_neighbor_data* );
+ simulation_box const * const, real, far_neighbor_data * const, int );
static void Choose_Neighbor_Counter( reax_system const * const system,
@@ -132,7 +132,7 @@ static inline real DistSqr_to_CP( rvec cp, rvec x )
/* Estimate the storage requirements for the far neighbor list */
-int Estimate_Num_Neighbors( reax_system * const system,
+void Estimate_Num_Neighbors( reax_system const * const system,
control_params const * const control, static_storage * const workspace,
reax_list ** const lists )
{
@@ -143,16 +143,13 @@ int Estimate_Num_Neighbors( reax_system * const system,
int *nbr_atoms;
ivec *nbrs;
rvec *nbrs_cp;
- grid *g;
+ grid const * const g = &system->g;
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 );
-
/* for each cell in the grid along the 3
* Cartesian directions: (i, j, k) => (x, y, z) */
for ( i = 0; i < g->ncell[0]; i++ )
@@ -211,46 +208,34 @@ int Estimate_Num_Neighbors( reax_system * const system,
}
}
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "[INFO] Estimate_Num_Neighbors: num_far = %d\n",
- (int) CEIL( num_far * SAFE_ZONE ) );
-#endif
-
- return (int) CEIL( num_far * SAFE_ZONE );
+ workspace->realloc.total_far_nbrs = (int) CEIL( num_far * SAFE_ZONE );
}
/* Generate the far neighbor list */
-void Generate_Neighbor_Lists( reax_system * const system,
+int Generate_Neighbor_Lists( reax_system * const system,
control_params const * const control, simulation_data * const data,
static_storage * const workspace, reax_list ** const lists )
{
int i, j, k, l, m, itr;
int x, y, z;
int atom1, atom2, max;
- int num_far, count;
+ int num_far, count, ret;
int *nbr_atoms;
ivec *nbrs;
rvec *nbrs_cp;
- grid *g;
+ grid const * const g = &system->g;
reax_list *far_nbrs;
- far_neighbor_data *nbr_data;
find_far_neighbors_function Find_Far_Neighbors;
real t_start, t_elapsed;
t_start = Get_Time( );
- g = &system->g;
num_far = 0;
far_nbrs = lists[FAR_NBRS];
+ ret = SUCCESS;
Choose_Neighbor_Finder( system, control, &Find_Far_Neighbors );
- Bin_Atoms( system, workspace );
-
-#if defined(REORDER_ATOMS)
- Reorder_Atoms( system, workspace, control );
-#endif
-
for ( i = 0; i < far_nbrs->n; ++i )
{
Set_Start_Index( i, 0, far_nbrs );
@@ -300,11 +285,11 @@ void Generate_Neighbor_Lists( reax_system * const system,
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 );
+ &system->box, control->vlist_cut,
+ &far_nbrs->far_nbr_list[num_far],
+ far_nbrs->total_intrs - num_far );
num_far += count;
}
@@ -315,11 +300,6 @@ void Generate_Neighbor_Lists( reax_system * const system,
}
Set_End_Index( atom1, num_far, far_nbrs );
-
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "[INFO] Generate_Neighbor_Lists: i = %d, start = %d, end = %d, itr = %d\n",
- atom1, Start_Index(atom1,far_nbrs), End_Index(atom1,far_nbrs), itr );
-#endif
}
}
}
@@ -331,16 +311,9 @@ void Generate_Neighbor_Lists( reax_system * const system,
ivec_MakeZero( system->atoms[i].rel_map );
}
- if ( num_far > far_nbrs->total_intrs * DANGER_ZONE )
+ if ( num_far > far_nbrs->total_intrs )
{
- workspace->realloc.num_far = num_far;
-
- if ( num_far > far_nbrs->total_intrs )
- {
- fprintf( stderr, "[ERROR] Generate_Neighbor_Lists: step%d-ran out of space on far_nbrs: top=%d, max=%d",
- data->step, num_far, far_nbrs->total_intrs );
- exit( INSUFFICIENT_MEMORY );
- }
+ ret = FAILURE;
}
#if defined(DEBUG_FOCUS)
@@ -358,4 +331,6 @@ void Generate_Neighbor_Lists( reax_system * const system,
t_elapsed = Get_Timing_Info( t_start );
data->timing.nbrs += t_elapsed;
+
+ return ret;
}
diff --git a/sPuReMD/src/neighbors.h b/sPuReMD/src/neighbors.h
index 45292a9d..24047cac 100644
--- a/sPuReMD/src/neighbors.h
+++ b/sPuReMD/src/neighbors.h
@@ -25,11 +25,11 @@
#include "reax_types.h"
-void Generate_Neighbor_Lists( reax_system * const, control_params const * const,
- simulation_data * const, static_storage * const, reax_list ** const );
-
-int Estimate_Num_Neighbors( reax_system * const, control_params const * const,
+void Estimate_Num_Neighbors( reax_system const * const, control_params const * const,
static_storage * const, reax_list ** const );
+int Generate_Neighbor_Lists( reax_system * const, control_params const * const,
+ simulation_data * const, static_storage * const, reax_list ** const );
+
#endif
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 286a0d29..7cf2fabf 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -70,12 +70,12 @@
#if defined(USE_REF_FORTRAN_REAXFF_CONSTANTS)
/* transcendental constant pi */
#define PI (3.14159265)
- /* unit conversion from ??? to kcal / mol */
+ /* unit conversion, ??? to kcal / mol */
#define C_ELE (332.0638)
- /* Boltzmann constant, AMU * A^2 / (ps^2 * K) */
+ /* Boltzmann constant, in J / mol / K */
// #define K_B (0.831687)
#define K_B (0.8314510)
- /* unit conversion for atomic force to AMU * A / ps^2 */
+ /* unit conversion for atomic force, kcal to J */
#define F_CONV (4.184e2)
/* energy conversion constant from electron volts to kilo-calories per mole */
#define KCALpMOL_to_EV (23.02)
@@ -110,13 +110,13 @@
#if !defined(K_B)
/* in ??? */
// #define K_B (503.398008)
- /* Boltzmann constant, AMU * A^2 / (ps^2 * K) */
+ /* in J / mol / K */
// #define K_B (0.831687)
#define K_B (0.8314510)
#endif
/* unit conversion for atomic force */
#if !defined(F_CONV)
- /* to AMU * A / ps^2 */
+ /* ??? to AMU * A / ps^2 */
#define F_CONV (1.0e6 / 48.88821291 / 48.88821291)
#endif
/* unit conversion for atomic energy */
@@ -196,6 +196,10 @@
/* min. temperature scaler for atomic positions and velocities in NPT ensembles */
#define MIN_dT (0.0)
+/* max. num. of main simulation loop retries;
+ * retries occur when memory allocation checks determine more memory is needed */
+#define MAX_RETRIES (5)
+
#define ALMOST_ZERO (1.0e-10)
#define NEG_INF (-1.0e10)
#define NO_BOND (1.0e-3)
@@ -433,31 +437,32 @@ typedef struct spuremd_handle spuremd_handle;
/* function pointer for calculating a bonded interaction */
-typedef void (*interaction_function)( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls* );
+typedef void (*interaction_function)( reax_system *, control_params *,
+ simulation_data *, static_storage *, reax_list **, output_controls * );
/* function pointer for calculating pairwise atom distance */
-typedef real (*atom_distance_function)( simulation_box*,
+typedef real (*atom_distance_function)( simulation_box const * const,
rvec, rvec, ivec, ivec, ivec, rvec );
/* function pointer for updating atom positions */
-typedef void (*update_atom_position_function)( rvec, rvec, ivec, simulation_box* );
+typedef void (*update_atom_position_function)( rvec, rvec, ivec, simulation_box const * const );
#if defined(TEST_FORCES)
/* function pointers for printed bonded interactions */
-typedef void (*print_interaction)(reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+typedef void (*print_interaction)( reax_system *, control_params *, simulation_data *,
+ static_storage *, reax_list **, output_controls * );
#endif
/* function pointer for evolving the atomic system (i.e., updating the positions)
* given the pre-computed forces from the prescribed interactions */
-typedef void (*evolve_function)( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list**, output_controls* );
+typedef int (*evolve_function)( reax_system * const, control_params * const,
+ simulation_data * const, static_storage * const, reax_list ** const,
+ output_controls * const );
/* function pointer for a callback function to be triggered after
* completion of a simulation step -- useful for, e.g., the Python wrapper */
-typedef void (*callback_function)( int, reax_atom*, simulation_data* );
+typedef void (*callback_function)( int, reax_atom *, simulation_data * );
/* function pointer for writing trajectory file header */
-typedef int (*write_header_function)( reax_system*, control_params*,
- static_storage*, output_controls* );
+typedef int (*write_header_function)( reax_system *, control_params *,
+ static_storage *, output_controls * );
/* function pointer for apending a frame to the trajectory file */
-typedef int (*append_traj_frame_function)( reax_system*, control_params*,
- simulation_data*, static_storage*, reax_list **, output_controls* );
+typedef int (*append_traj_frame_function)( reax_system *, control_params *,
+ simulation_data *, static_storage *, reax_list **, output_controls * );
/* function pointer for writing to the trajectory file */
typedef int (*write_function)( FILE *, const char *, ... );
@@ -891,14 +896,21 @@ struct reax_system
unsigned int num_molec_charge_constraints;
/* max. num. of charge constraints on groups of atoms (i.e., molecules) */
unsigned int max_num_molec_charge_constraints;
- /* atom info */
- reax_atom *atoms;
/* atomic interaction parameters */
reax_interaction reax_param;
- /* simulation space (a.k.a. box) parameters */
+ /* grid specifying domain (i.e., spatial) decompisition
+ * of atoms within simulation box */
simulation_box box;
/* grid structure used for binning atoms and tracking neighboring bins */
grid g;
+ /* collection of atomic info. */
+ reax_atom *atoms;
+ /* num. bonds per atom */
+ int *bonds;
+ /* num. hydrogen bonds per atom */
+ int *hbonds;
+ /* num. matrix entries per row */
+ int *cm_entries;
};
@@ -1439,15 +1451,35 @@ struct sparse_matrix
};
+/* used to determine if and how much space should be reallocated */
struct reallocate_data
{
- int num_far;
- int Htop;
+ /* TRUE if far neighbor list needs
+ * to be reallocated, FALSE otherwise */
+ int far_nbrs;
+ /* total num. of (max) far neighbors across all atoms */
+ int total_far_nbrs;
+ /* TRUE if charge matrix needs
+ * to be reallocated, FALSE otherwise */
+ int cm;
+ /* total num. matrix entries (sum over max) */
+ int total_cm_entries;
+ /* TRUE if hbond list needs
+ * to be reallocated, FALSE otherwise */
int hbonds;
- int num_hbonds;
+ /* total num. hydrogen bonds (sum over max) */
+ int total_hbonds;
+ /* TRUE if bond list needs
+ * to be reallocated, FALSE otherwise */
int bonds;
- int num_bonds;
- int num_3body;
+ /* total num. bonds (sum over max) */
+ int total_bonds;
+ /* TRUE if three body list needs
+ * to be reallocated, FALSE otherwise */
+ int thbody;
+ /* total num. three body interactions */
+ int total_thbodies;
+ /**/
int gcell_atoms;
};
@@ -1664,23 +1696,21 @@ struct static_storage
};
-/* interaction lists */
+/* 1-D flattened list of lists used for various interaction types */
struct reax_list
{
- /* 0 if struct members are NOT allocated, 1 otherwise */
+ /* FALSE if struct members are NOT allocated, TRUE otherwise */
int allocated;
- /* num. active entries in list for this simulation */
+ /* num. active entries (i.e., nested lists) in the list for this simulation */
int n;
- /* max. num. entries in list across all simulations */
+ /* max. num. entries (i.e., nested lists) in the list across all simulations */
int n_max;
/* total interactions across all entries which can be stored in the list */
int total_intrs;
- /* starting position of atom's interactions */
+ /* starting position of interactions for an inner list */
int *index;
- /* ending position of atom's interactions */
+ /* ending position of interactions for an inner list */
int *end_index;
- /* max. num. of interactions per list entry */
- int *max_intrs;
/* interaction list (polymorphic via union dereference) */
// union
// {
diff --git a/sPuReMD/src/reset_tools.c b/sPuReMD/src/reset_tools.c
index 57ececef..6756bb87 100644
--- a/sPuReMD/src/reset_tools.c
+++ b/sPuReMD/src/reset_tools.c
@@ -200,7 +200,7 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
Set_End_Index( i, tmp, bonds );
}
- if ( control->hbond_cut > 0.0 )
+ if ( control->hbond_cut > 0.0 && workspace->num_H > 0 )
{
for ( i = 0; i < system->N; ++i )
{
@@ -268,7 +268,6 @@ void Reset_Grid( grid *g )
}
-
void Reset_Marks( grid *g, ivec *grid_stack, int grid_top )
{
int i;
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 913d3050..28214f53 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -353,7 +353,7 @@ int setup_callback( const void * const handle, const callback_function callback
*/
int simulate( const void * const handle )
{
- int steps, ret;
+ int steps, ret, ret_forces, retries;
evolve_function Evolve;
spuremd_handle *spmd_handle;
@@ -374,10 +374,23 @@ int simulate( const void * const handle )
Reset( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists );
- Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ ret_forces = Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->realloc );
+ if ( ret_forces != SUCCESS )
+ {
+ ret_forces = Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
+ spmd_handle->realloc );
+
+ if ( ret_forces != SUCCESS )
+ {
+ fprintf( stderr, "[ERROR] Unrecoverable memory allocation issue (Compute_Forces). Terminating...\n" );
+ exit( INVALID_GEO );
+ }
+ }
+
Compute_Kinetic_Energy( spmd_handle->system, spmd_handle->data );
if ( spmd_handle->control->compute_pressure == TRUE && spmd_handle->control->ensemble != sNPT
@@ -414,54 +427,76 @@ int simulate( const void * const handle )
}
//}
- for ( ++spmd_handle->data->step; spmd_handle->data->step <= spmd_handle->control->nsteps; spmd_handle->data->step++ )
+ ++spmd_handle->data->step;
+ retries = 0;
+ while ( spmd_handle->data->step <= spmd_handle->control->nsteps
+ && retries < MAX_RETRIES )
{
+ ret = SUCCESS;
+
if ( spmd_handle->control->T_mode != 0 )
{
Temperature_Control( spmd_handle->control, spmd_handle->data,
spmd_handle->out_control );
}
- Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
-
- Post_Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ ret = Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
- if ( spmd_handle->output_enabled == TRUE )
+ if ( ret == SUCCESS )
{
- Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ Post_Evolve( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
}
- Check_Energy( spmd_handle->data );
-
- if ( spmd_handle->output_enabled == TRUE )
+ if ( ret == SUCCESS )
{
- steps = spmd_handle->data->step - spmd_handle->data->prev_steps;
+ if ( spmd_handle->output_enabled == TRUE )
+ {
+ Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
+ }
- Analysis( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
+ Check_Energy( spmd_handle->data );
- if ( spmd_handle->out_control->restart_freq > 0
- && steps % spmd_handle->out_control->restart_freq == 0 )
+ if ( spmd_handle->output_enabled == TRUE )
{
- Write_Restart( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, spmd_handle->out_control );
+ steps = spmd_handle->data->step - spmd_handle->data->prev_steps;
+
+ Analysis( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
+
+ if ( spmd_handle->out_control->restart_freq > 0
+ && steps % spmd_handle->out_control->restart_freq == 0 )
+ {
+ Write_Restart( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ spmd_handle->workspace, spmd_handle->out_control );
+ }
+
+ if ( spmd_handle->out_control->write_steps > 0
+ && steps % spmd_handle->out_control->write_steps == 0 )
+ {
+ Write_PDB( spmd_handle->system, spmd_handle->lists[BONDS], spmd_handle->data,
+ spmd_handle->control, spmd_handle->workspace, spmd_handle->out_control );
+ }
}
- if ( spmd_handle->out_control->write_steps > 0
- && steps % spmd_handle->out_control->write_steps == 0 )
+ if ( spmd_handle->callback != NULL )
{
- Write_PDB( spmd_handle->system, spmd_handle->lists[BONDS], spmd_handle->data,
- spmd_handle->control, spmd_handle->workspace, spmd_handle->out_control );
+ spmd_handle->callback( spmd_handle->system->N, spmd_handle->system->atoms,
+ spmd_handle->data );
}
- }
- if ( spmd_handle->callback != NULL )
+ ++spmd_handle->data->step;
+ retries = 0;
+ }
+ else
{
- spmd_handle->callback( spmd_handle->system->N, spmd_handle->system->atoms,
- spmd_handle->data );
+ ++retries;
+
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr, "[INFO] p%d: retrying step %d...\n", system->my_rank, data->step );
+#endif
}
}
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index 57a060a8..0829259b 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -59,7 +59,8 @@ void Temperature_Control( control_params *control, simulation_data *data,
}
-void Compute_Total_Mass( reax_system *system, simulation_data *data )
+void Compute_Total_Mass( reax_system const * const system,
+ simulation_data * const data )
{
int i;
@@ -74,7 +75,8 @@ void Compute_Total_Mass( reax_system *system, simulation_data *data )
}
-void Compute_Center_of_Mass( reax_system *system, simulation_data *data )
+void Compute_Center_of_Mass( reax_system const * const system,
+ simulation_data * const data )
{
int i;
real m, xx, xy, xz, yy, yz, zz, det;
@@ -181,7 +183,8 @@ void Compute_Center_of_Mass( reax_system *system, simulation_data *data )
}
-void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
+void Compute_Kinetic_Energy( reax_system const * const system,
+ simulation_data * const data )
{
int i;
real m;
diff --git a/sPuReMD/src/system_props.h b/sPuReMD/src/system_props.h
index 26e1594b..6fccbc4a 100644
--- a/sPuReMD/src/system_props.h
+++ b/sPuReMD/src/system_props.h
@@ -27,11 +27,11 @@
void Temperature_Control( control_params*, simulation_data*, output_controls* );
-void Compute_Total_Mass( reax_system*, simulation_data* );
+void Compute_Total_Mass( reax_system const * const, simulation_data * const );
-void Compute_Center_of_Mass( reax_system*, simulation_data* );
+void Compute_Center_of_Mass( reax_system const * const, simulation_data * const );
-void Compute_Kinetic_Energy( reax_system*, simulation_data* );
+void Compute_Kinetic_Energy( reax_system const * const , simulation_data * const );
void Compute_Total_Energy( simulation_data* );
diff --git a/sPuReMD/src/valence_angles.c b/sPuReMD/src/valence_angles.c
index 5f094610..1a85372f 100644
--- a/sPuReMD/src/valence_angles.c
+++ b/sPuReMD/src/valence_angles.c
@@ -718,22 +718,26 @@ void Valence_Angles( reax_system *system, control_params *control,
if ( num_thb_intrs >= thb_intrs->total_intrs * DANGER_ZONE )
{
- workspace->realloc.num_3body = num_thb_intrs;
+ workspace->realloc.total_thbodies = (int) CEIL( num_thb_intrs * SAFE_ZONE );
+ workspace->realloc.thbody = TRUE;
+
+ /* retry functionality is not implemented as valence angle list
+ * computation is not an atomic transaction, so use older reallocation strategy for now */
if ( num_thb_intrs > thb_intrs->total_intrs )
{
- fprintf( stderr, "step%d-ran out of space on angle_list: top=%d, max=%d",
+ fprintf( stderr, "[ERROR] step%d-ran out of space on angle_list: top=%d, max=%d",
data->step, num_thb_intrs, thb_intrs->total_intrs );
exit( INSUFFICIENT_MEMORY );
}
}
#if defined(TEST_ENERGY)
- fprintf( stderr, "Number of angle interactions: %d\n", num_thb_intrs );
+ fprintf( stderr, "[INFO] Number of angle interactions: %d\n", num_thb_intrs );
- fprintf( stderr, "Angle Energy: %g\t Penalty Energy: %g\t Coalition Energy: %g\n",
+ fprintf( stderr, "[INFO] Angle Energy: %g\t Penalty Energy: %g\t Coalition Energy: %g\n",
data->E_Ang, data->E_Pen, data->E_Coa );
- fprintf( stderr, "3body: press (%23.15e %23.15e %23.15e)\n",
+ fprintf( stderr, "[INFO] 3body: press (%23.15e %23.15e %23.15e)\n",
data->press[0][0], data->press[1][1], data->press[2][2] );
#endif
}
--
GitLab