diff --git a/src/USER-REAXC/reaxc_allocate.cpp b/src/USER-REAXC/reaxc_allocate.cpp
index 7df5241a9..2d571269d 100644
--- a/src/USER-REAXC/reaxc_allocate.cpp
+++ b/src/USER-REAXC/reaxc_allocate.cpp
@@ -1,1026 +1,884 @@
/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Copyright (2010) Purdue University
Hasan Metin Aktulga, hmaktulga@lbl.gov
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
Please cite the related publication:
H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama,
"Parallel Reactive Molecular Dynamics: Numerical Methods and
Algorithmic Techniques", Parallel Computing, in press.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details:
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "pair_reax_c.h"
-#if defined(PURE_REAX)
-#include "allocate.h"
-#include "list.h"
-#include "reset_tools.h"
-#include "tool_box.h"
-#include "vector.h"
-#elif defined(LAMMPS_REAX)
#include "reaxc_allocate.h"
#include "reaxc_list.h"
#include "reaxc_reset_tools.h"
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
-#endif
/* allocate space for my_atoms
important: we cannot know the exact number of atoms that will fall into a
process's box throughout the whole simulation. therefore
we need to make upper bound estimates for various data structures */
int PreAllocate_Space( reax_system *system, control_params *control,
storage *workspace, MPI_Comm comm )
{
int mincap = system->mincap;
double safezone = system->safezone;
// determine the local and total capacity
system->local_cap = MAX( (int)(system->n * safezone), mincap );
system->total_cap = MAX( (int)(system->N * safezone), mincap );
#if defined(DEBUG)
fprintf( stderr, "p%d: local_cap=%d total_cap=%d\n",
system->my_rank, system->local_cap, system->total_cap );
#endif
system->my_atoms = (reax_atom*)
scalloc( system->total_cap, sizeof(reax_atom), "my_atoms", comm );
/* space for keeping restriction info, if any */
// not yet implemented in the parallel version!!!
// if( control->restrict_bonds ) {
// workspace->restricted = (int*)
// scalloc( system->local_cap, sizeof(int), "restricted_atoms", comm );
// workspace->restricted_list = (int**)
// scalloc( system->local_cap, sizeof(int*), "restricted_list", comm );
// for( i = 0; i < system->local_cap; ++i )
// workspace->restricted_list[i] = (int*)
// scalloc( MAX_RESTRICT, sizeof(int), "restricted_list[i]", comm );
// }
return SUCCESS;
}
/************* system *************/
inline void reax_atom_Copy( reax_atom *dest, reax_atom *src )
{
dest->orig_id = src->orig_id;
dest->type = src->type;
strcpy( dest->name, src->name );
rvec_Copy( dest->x, src->x );
rvec_Copy( dest->v, src->v );
rvec_Copy( dest->f_old, src->f_old );
rvec_Copy( dest->s, src->s );
rvec_Copy( dest->t, src->t );
dest->Hindex = src->Hindex;
dest->num_bonds = src->num_bonds;
dest->num_hbonds = src->num_hbonds;
dest->numbonds = src->numbonds;
}
void Copy_Atom_List( reax_atom *dest, reax_atom *src, int n )
{
int i;
for( i = 0; i < n; ++i )
memcpy( dest+i, src+i, sizeof(reax_atom) );
}
int Allocate_System( reax_system *system, int local_cap, int total_cap,
char *msg )
{
system->my_atoms = (reax_atom*)
realloc( system->my_atoms, total_cap*sizeof(reax_atom) );
return SUCCESS;
}
void DeAllocate_System( reax_system *system )
{
int i, j, k;
int ntypes;
reax_interaction *ff_params;
// dealloocate the atom list
sfree( system->my_atoms, "system->my_atoms" );
// deallocate the ffield parameters storage
ff_params = &(system->reax_param);
ntypes = ff_params->num_atom_types;
sfree( ff_params->gp.l, "ff:globals" );
for( i = 0; i < ntypes; ++i ) {
for( j = 0; j < ntypes; ++j ) {
for( k = 0; k < ntypes; ++k ) {
sfree( ff_params->fbp[i][j][k], "ff:fbp[i,j,k]" );
}
sfree( ff_params->fbp[i][j], "ff:fbp[i,j]" );
sfree( ff_params->thbp[i][j], "ff:thbp[i,j]" );
sfree( ff_params->hbp[i][j], "ff:hbp[i,j]" );
}
sfree( ff_params->fbp[i], "ff:fbp[i]" );
sfree( ff_params->thbp[i], "ff:thbp[i]" );
sfree( ff_params->hbp[i], "ff:hbp[i]" );
sfree( ff_params->tbp[i], "ff:tbp[i]" );
}
sfree( ff_params->fbp, "ff:fbp" );
sfree( ff_params->thbp, "ff:thbp" );
sfree( ff_params->hbp, "ff:hbp" );
sfree( ff_params->tbp, "ff:tbp" );
sfree( ff_params->sbp, "ff:sbp" );
}
/************* workspace *************/
void DeAllocate_Workspace( control_params *control, storage *workspace )
{
int i;
if( !workspace->allocated )
return;
workspace->allocated = 0;
/* communication storage */
for( i = 0; i < MAX_NBRS; ++i ) {
sfree( workspace->tmp_dbl[i], "tmp_dbl[i]" );
sfree( workspace->tmp_rvec[i], "tmp_rvec[i]" );
sfree( workspace->tmp_rvec2[i], "tmp_rvec2[i]" );
}
/* bond order storage */
sfree( workspace->within_bond_box, "skin" );
sfree( workspace->total_bond_order, "total_bo" );
sfree( workspace->Deltap, "Deltap" );
sfree( workspace->Deltap_boc, "Deltap_boc" );
sfree( workspace->dDeltap_self, "dDeltap_self" );
sfree( workspace->Delta, "Delta" );
sfree( workspace->Delta_lp, "Delta_lp" );
sfree( workspace->Delta_lp_temp, "Delta_lp_temp" );
sfree( workspace->dDelta_lp, "dDelta_lp" );
sfree( workspace->dDelta_lp_temp, "dDelta_lp_temp" );
sfree( workspace->Delta_e, "Delta_e" );
sfree( workspace->Delta_boc, "Delta_boc" );
sfree( workspace->Delta_val, "Delta_val" );
sfree( workspace->nlp, "nlp" );
sfree( workspace->nlp_temp, "nlp_temp" );
sfree( workspace->Clp, "Clp" );
sfree( workspace->vlpex, "vlpex" );
sfree( workspace->bond_mark, "bond_mark" );
sfree( workspace->done_after, "done_after" );
/* QEq storage */
sfree( workspace->Hdia_inv, "Hdia_inv" );
sfree( workspace->b_s, "b_s" );
sfree( workspace->b_t, "b_t" );
sfree( workspace->b_prc, "b_prc" );
sfree( workspace->b_prm, "b_prm" );
sfree( workspace->s, "s" );
sfree( workspace->t, "t" );
sfree( workspace->droptol, "droptol" );
sfree( workspace->b, "b" );
sfree( workspace->x, "x" );
/* GMRES storage */
for( i = 0; i < RESTART+1; ++i ) {
sfree( workspace->h[i], "h[i]" );
sfree( workspace->v[i], "v[i]" );
}
sfree( workspace->h, "h" );
sfree( workspace->v, "v" );
sfree( workspace->y, "y" );
sfree( workspace->z, "z" );
sfree( workspace->g, "g" );
sfree( workspace->hs, "hs" );
sfree( workspace->hc, "hc" );
/* CG storage */
sfree( workspace->r, "r" );
sfree( workspace->d, "d" );
sfree( workspace->q, "q" );
sfree( workspace->p, "p" );
sfree( workspace->r2, "r2" );
sfree( workspace->d2, "d2" );
sfree( workspace->q2, "q2" );
sfree( workspace->p2, "p2" );
/* integrator */
// sfree( workspace->f_old );
sfree( workspace->v_const, "v_const" );
/*workspace->realloc.num_far = -1;
workspace->realloc.Htop = -1;
workspace->realloc.hbonds = -1;
workspace->realloc.bonds = -1;
workspace->realloc.num_3body = -1;
workspace->realloc.gcell_atoms = -1;*/
/* storage for analysis */
// if( control->molecular_analysis || control->diffusion_coef ) {
// sfree( workspace->mark, "mark" );
// sfree( workspace->old_mark, "old_mark" );
// }
// if( control->diffusion_coef )
// sfree( workspace->x_old, "x_old" );
/* force related storage */
sfree( workspace->f, "f" );
sfree( workspace->CdDelta, "CdDelta" );
#ifdef TEST_FORCES
sfree(workspace->dDelta, "dDelta" );
sfree( workspace->f_ele, "f_ele" );
sfree( workspace->f_vdw, "f_vdw" );
sfree( workspace->f_bo, "f_bo" );
sfree( workspace->f_be, "f_be" );
sfree( workspace->f_lp, "f_lp" );
sfree( workspace->f_ov, "f_ov" );
sfree( workspace->f_un, "f_un" );
sfree( workspace->f_ang, "f_ang" );
sfree( workspace->f_coa, "f_coa" );
sfree( workspace->f_pen, "f_pen" );
sfree( workspace->f_hb, "f_hb" );
sfree( workspace->f_tor, "f_tor" );
sfree( workspace->f_con, "f_con" );
sfree( workspace->f_tot, "f_tot" );
sfree( workspace->rcounts, "rcounts" );
sfree( workspace->displs, "displs" );
sfree( workspace->id_all, "id_all" );
sfree( workspace->f_all, "f_all" );
#endif
/* hbond storage */
//sfree( workspace->Hindex, "Hindex" );
//sfree( workspace->num_bonds );
//sfree( workspace->num_hbonds );
//sfree( workspace->hash, "hash" );
//sfree( workspace->rev_hash, "rev_hash" );
}
int Allocate_Workspace( reax_system *system, control_params *control,
storage *workspace, int local_cap, int total_cap,
MPI_Comm comm, char *msg )
{
int i, total_real, total_rvec, local_rvec;
workspace->allocated = 1;
total_real = total_cap * sizeof(real);
total_rvec = total_cap * sizeof(rvec);
local_rvec = local_cap * sizeof(rvec);
/* communication storage */
for( i = 0; i < MAX_NBRS; ++i ) {
workspace->tmp_dbl[i] = (real*)
scalloc( total_cap, sizeof(real), "tmp_dbl", comm );
workspace->tmp_rvec[i] = (rvec*)
scalloc( total_cap, sizeof(rvec), "tmp_rvec", comm );
workspace->tmp_rvec2[i] = (rvec2*)
scalloc( total_cap, sizeof(rvec2), "tmp_rvec2", comm );
}
/* bond order related storage */
workspace->within_bond_box = (int*)
scalloc( total_cap, sizeof(int), "skin", comm );
workspace->total_bond_order = (real*) smalloc( total_real, "total_bo", comm );
workspace->Deltap = (real*) smalloc( total_real, "Deltap", comm );
workspace->Deltap_boc = (real*) smalloc( total_real, "Deltap_boc", comm );
workspace->dDeltap_self = (rvec*) smalloc( total_rvec, "dDeltap_self", comm );
workspace->Delta = (real*) smalloc( total_real, "Delta", comm );
workspace->Delta_lp = (real*) smalloc( total_real, "Delta_lp", comm );
workspace->Delta_lp_temp = (real*)
smalloc( total_real, "Delta_lp_temp", comm );
workspace->dDelta_lp = (real*) smalloc( total_real, "dDelta_lp", comm );
workspace->dDelta_lp_temp = (real*)
smalloc( total_real, "dDelta_lp_temp", comm );
workspace->Delta_e = (real*) smalloc( total_real, "Delta_e", comm );
workspace->Delta_boc = (real*) smalloc( total_real, "Delta_boc", comm );
workspace->Delta_val = (real*) smalloc( total_real, "Delta_val", comm );
workspace->nlp = (real*) smalloc( total_real, "nlp", comm );
workspace->nlp_temp = (real*) smalloc( total_real, "nlp_temp", comm );
workspace->Clp = (real*) smalloc( total_real, "Clp", comm );
workspace->vlpex = (real*) smalloc( total_real, "vlpex", comm );
workspace->bond_mark = (int*)
scalloc( total_cap, sizeof(int), "bond_mark", comm );
workspace->done_after = (int*)
scalloc( total_cap, sizeof(int), "done_after", comm );
// fprintf( stderr, "p%d: bond order storage\n", system->my_rank );
/* QEq storage */
workspace->Hdia_inv = (real*)
scalloc( total_cap, sizeof(real), "Hdia_inv", comm );
workspace->b_s = (real*) scalloc( total_cap, sizeof(real), "b_s", comm );
workspace->b_t = (real*) scalloc( total_cap, sizeof(real), "b_t", comm );
workspace->b_prc = (real*) scalloc( total_cap, sizeof(real), "b_prc", comm );
workspace->b_prm = (real*) scalloc( total_cap, sizeof(real), "b_prm", comm );
workspace->s = (real*) scalloc( total_cap, sizeof(real), "s", comm );
workspace->t = (real*) scalloc( total_cap, sizeof(real), "t", comm );
workspace->droptol = (real*)
scalloc( total_cap, sizeof(real), "droptol", comm );
workspace->b = (rvec2*) scalloc( total_cap, sizeof(rvec2), "b", comm );
workspace->x = (rvec2*) scalloc( total_cap, sizeof(rvec2), "x", comm );
/* GMRES storage */
workspace->y = (real*) scalloc( RESTART+1, sizeof(real), "y", comm );
workspace->z = (real*) scalloc( RESTART+1, sizeof(real), "z", comm );
workspace->g = (real*) scalloc( RESTART+1, sizeof(real), "g", comm );
workspace->h = (real**) scalloc( RESTART+1, sizeof(real*), "h", comm );
workspace->hs = (real*) scalloc( RESTART+1, sizeof(real), "hs", comm );
workspace->hc = (real*) scalloc( RESTART+1, sizeof(real), "hc", comm );
workspace->v = (real**) scalloc( RESTART+1, sizeof(real*), "v", comm );
for( i = 0; i < RESTART+1; ++i ) {
workspace->h[i] = (real*) scalloc( RESTART+1, sizeof(real), "h[i]", comm );
workspace->v[i] = (real*) scalloc( total_cap, sizeof(real), "v[i]", comm );
}
/* CG storage */
workspace->r = (real*) scalloc( total_cap, sizeof(real), "r", comm );
workspace->d = (real*) scalloc( total_cap, sizeof(real), "d", comm );
workspace->q = (real*) scalloc( total_cap, sizeof(real), "q", comm );
workspace->p = (real*) scalloc( total_cap, sizeof(real), "p", comm );
workspace->r2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "r2", comm );
workspace->d2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "d2", comm );
workspace->q2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "q2", comm );
workspace->p2 = (rvec2*) scalloc( total_cap, sizeof(rvec2), "p2", comm );
/* integrator storage */
workspace->v_const = (rvec*) smalloc( local_rvec, "v_const", comm );
/* storage for analysis */
// not yet implemented in the parallel version!!!
// if( control->molecular_analysis || control->diffusion_coef ) {
// workspace->mark = (int*) scalloc( local_cap, sizeof(int), "mark", comm );
// workspace->old_mark = (int*)
// scalloc( local_cap, sizeof(int), "old_mark", comm );
// }
// else
// workspace->mark = workspace->old_mark = NULL;
// if( control->diffusion_coef )
// workspace->x_old = (rvec*)
// scalloc( local_cap, sizeof(rvec), "x_old", comm );
// else workspace->x_old = NULL;
// /* force related storage */
workspace->f = (rvec*) scalloc( total_cap, sizeof(rvec), "f", comm );
workspace->CdDelta = (real*)
scalloc( total_cap, sizeof(real), "CdDelta", comm );
#ifdef TEST_FORCES
workspace->dDelta=(rvec*) smalloc( total_rvec, "dDelta", comm );
workspace->f_ele =(rvec*) smalloc( total_rvec, "f_ele", comm );
workspace->f_vdw =(rvec*) smalloc( total_rvec, "f_vdw", comm );
workspace->f_bo =(rvec*) smalloc( total_rvec, "f_bo", comm );
workspace->f_be =(rvec*) smalloc( total_rvec, "f_be", comm );
workspace->f_lp =(rvec*) smalloc( total_rvec, "f_lp", comm );
workspace->f_ov =(rvec*) smalloc( total_rvec, "f_ov", comm );
workspace->f_un =(rvec*) smalloc( total_rvec, "f_un", comm );
workspace->f_ang =(rvec*) smalloc( total_rvec, "f_ang", comm );
workspace->f_coa =(rvec*) smalloc( total_rvec, "f_coa", comm );
workspace->f_pen =(rvec*) smalloc( total_rvec, "f_pen", comm );
workspace->f_hb =(rvec*) smalloc( total_rvec, "f_hb", comm );
workspace->f_tor =(rvec*) smalloc( total_rvec, "f_tor", comm );
workspace->f_con =(rvec*) smalloc( total_rvec, "f_con", comm );
workspace->f_tot =(rvec*) smalloc( total_rvec, "f_tot", comm );
if( system->my_rank == MASTER_NODE ) {
workspace->rcounts = (int*)
smalloc( system->wsize*sizeof(int), "rcount", comm );
workspace->displs = (int*)
smalloc( system->wsize*sizeof(int), "displs", comm );
workspace->id_all = (int*)
smalloc( system->bigN*sizeof(int), "id_all", comm );
workspace->f_all = (rvec*)
smalloc( system->bigN*sizeof(rvec), "f_all", comm );
}
else{
workspace->rcounts = NULL;
workspace->displs = NULL;
workspace->id_all = NULL;
workspace->f_all = NULL;
}
#endif
return SUCCESS;
}
void Reallocate_Neighbor_List( reax_list *far_nbrs, int n, int num_intrs,
MPI_Comm comm )
{
Delete_List( far_nbrs, comm );
if(!Make_List( n, num_intrs, TYP_FAR_NEIGHBOR, far_nbrs, comm )){
fprintf(stderr, "Problem in initializing far nbrs list. Terminating!\n");
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
}
int Allocate_Matrix( sparse_matrix **pH, int cap, int m, MPI_Comm comm )
{
sparse_matrix *H;
*pH = (sparse_matrix*)
smalloc( sizeof(sparse_matrix), "sparse_matrix", comm );
H = *pH;
H->cap = cap;
H->m = m;
H->start = (int*) smalloc( sizeof(int) * cap, "matrix_start", comm );
H->end = (int*) smalloc( sizeof(int) * cap, "matrix_end", comm );
H->entries = (sparse_matrix_entry*)
smalloc( sizeof(sparse_matrix_entry)*m, "matrix_entries", comm );
return SUCCESS;
}
void Deallocate_Matrix( sparse_matrix *H )
{
sfree(H->start, "H->start");
sfree(H->end, "H->end");
sfree(H->entries, "H->entries");
sfree(H, "H");
}
int Reallocate_Matrix( sparse_matrix **H, int n, int m, char *name,
MPI_Comm comm )
{
Deallocate_Matrix( *H );
if( !Allocate_Matrix( H, n, m, comm ) ) {
fprintf(stderr, "not enough space for %s matrix. terminating!\n", name);
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "reallocating %s matrix, n = %d, m = %d\n", name, n, m );
fprintf( stderr, "memory allocated: %s = %dMB\n",
name, (int)(m * sizeof(sparse_matrix_entry) / (1024*1024)) );
#endif
return SUCCESS;
}
int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
MPI_Comm comm )
{
int i, id, total_hbonds;
int mincap = system->mincap;
double saferzone = system->saferzone;
total_hbonds = 0;
for( i = 0; i < system->n; ++i )
if( (id = system->my_atoms[i].Hindex) >= 0 ) {
// commented out - already updated in validate_lists in forces.c
// system->my_atoms[i].num_hbonds = MAX(Num_Entries(id,hbonds)*SAFER_ZONE,
// MIN_HBONDS);
total_hbonds += system->my_atoms[i].num_hbonds;
}
total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
Delete_List( hbonds, comm );
if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND, hbonds, comm ) ) {
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
return total_hbonds;
}
int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
int *total_bonds, int *est_3body, MPI_Comm comm )
{
int i;
int mincap = system->mincap;
double safezone = system->safezone;
*total_bonds = 0;
*est_3body = 0;
for( i = 0; i < system->N; ++i ){
*est_3body += SQR(system->my_atoms[i].num_bonds);
// commented out - already updated in validate_lists in forces.c
// system->my_atoms[i].num_bonds = MAX( Num_Entries(i,bonds)*2, MIN_BONDS );
*total_bonds += system->my_atoms[i].num_bonds;
}
*total_bonds = (int)(MAX( *total_bonds * safezone, mincap*MIN_BONDS ));
Delete_List( bonds, comm );
if(!Make_List(system->total_cap, *total_bonds, TYP_BOND, bonds, comm)) {
fprintf( stderr, "not enough space for bonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
return SUCCESS;
}
/************* grid *************/
int Estimate_GCell_Population( reax_system* system, MPI_Comm comm )
{
int d, i, j, k, l, max_atoms, my_max, all_max;
ivec c;
grid *g;
grid_cell *gc;
simulation_box *my_ext_box;
reax_atom *atoms;
my_ext_box = &(system->my_ext_box);
g = &(system->my_grid);
atoms = system->my_atoms;
Reset_Grid( g );
for( l = 0; l < system->n; l++ ) {
for( d = 0; d < 3; ++d ) {
c[d] = (int)((atoms[l].x[d]-my_ext_box->min[d])*g->inv_len[d]);
if( c[d] >= g->native_end[d] )
c[d] = g->native_end[d] - 1;
else if( c[d] < g->native_str[d] )
c[d] = g->native_str[d];
}
#if defined(DEBUG)
fprintf( stderr, "p%d bin_my_atoms: l:%d - atom%d @ %.5f %.5f %.5f" \
"--> cell: %d %d %d\n",
system->my_rank, l, atoms[l].orig_id,
atoms[l].x[0], atoms[l].x[1], atoms[l].x[2],
c[0], c[1], c[2] );
#endif
gc = &( g->cells[c[0]][c[1]][c[2]] );
gc->top++;
}
max_atoms = 0;
for( i = 0; i < g->ncells[0]; i++ )
for( j = 0; j < g->ncells[1]; j++ )
for( k = 0; k < g->ncells[2]; k++ ) {
gc = &(g->cells[i][j][k]);
if( max_atoms < gc->top )
max_atoms = gc->top;
#if defined(DEBUG)
fprintf( stderr, "p%d gc[%d,%d,%d]->top=%d\n",
system->my_rank, i, j, k, gc->top );
#endif
}
my_max = (int)(MAX(max_atoms*SAFE_ZONE, MIN_GCELL_POPL));
MPI_Allreduce( &my_max, &all_max, 1, MPI_INT, MPI_MAX, comm );
#if defined(DEBUG)
fprintf( stderr, "p%d max_atoms=%d, my_max=%d, all_max=%d\n",
system->my_rank, max_atoms, my_max, all_max );
#endif
return all_max;
}
void Allocate_Grid( reax_system *system, MPI_Comm comm )
{
int i, j, k, l;
grid *g;
grid_cell *gc;
g = &( system->my_grid );
/* allocate gcell reordering space */
g->order = (ivec*) scalloc( g->total+1, sizeof(ivec), "g:order", comm );
/* allocate the gcells for the new grid */
g->max_nbrs = (2*g->vlist_span[0]+1)*(2*g->vlist_span[1]+1)*
(2*g->vlist_span[2]+1)+3;
g->cells = (grid_cell***)
scalloc( g->ncells[0], sizeof(grid_cell**), "gcells", comm );
for( i = 0; i < g->ncells[0]; i++ ) {
g->cells[i] = (grid_cell**)
scalloc( g->ncells[1], sizeof(grid_cell*),"gcells[i]", comm );
for( j = 0; j < g->ncells[1]; ++j ) {
g->cells[i][j] = (grid_cell*)
scalloc( g->ncells[2], sizeof(grid_cell), "gcells[i][j]", comm );
for( k = 0; k < g->ncells[2]; k++ ) {
gc = &(g->cells[i][j][k]);
gc->top = gc->mark = gc->str = gc->end = 0;
gc->nbrs = (grid_cell**)
scalloc( g->max_nbrs, sizeof(grid_cell*), "g:nbrs", comm );
gc->nbrs_x = (ivec*)
scalloc( g->max_nbrs, sizeof(ivec), "g:nbrs_x", comm );
gc->nbrs_cp = (rvec*)
scalloc( g->max_nbrs, sizeof(rvec), "g:nbrs_cp", comm );
for( l = 0; l < g->max_nbrs; ++l )
gc->nbrs[l] = NULL;
}
}
}
/* allocate atom id storage in gcells */
g->max_atoms = Estimate_GCell_Population( system, comm );
/* space for storing atom id's is required only for native cells */
for( i = g->native_str[0]; i < g->native_end[0]; ++i )
for( j = g->native_str[1]; j < g->native_end[1]; ++j )
for( k = g->native_str[2]; k < g->native_end[2]; ++k )
g->cells[i][j][k].atoms = (int*) scalloc( g->max_atoms, sizeof(int),
"g:atoms", comm );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d-allocated %dx%dx%d grid: nbrs=%d atoms=%d space=%dMB\n",
system->my_rank, g->ncells[0], g->ncells[1], g->ncells[2],
g->max_nbrs, g->max_atoms,
(int)
((g->total*sizeof(grid_cell)+g->total*g->max_nbrs*sizeof(int*) +
g->total*g->max_nbrs*sizeof(rvec) +
(g->native_end[0]-g->native_str[0])*
(g->native_end[1]-g->native_str[1])*
(g->native_end[2]-g->native_str[2])*g->max_atoms*sizeof(int))/
(1024*1024)) );
#endif
}
void Deallocate_Grid( grid *g )
{
int i, j, k;
grid_cell *gc;
sfree( g->order, "g:order" );
/* deallocate the grid cells */
for( i = 0; i < g->ncells[0]; i++ ) {
for( j = 0; j < g->ncells[1]; j++ ) {
for( k = 0; k < g->ncells[2]; k++ ) {
gc = &(g->cells[i][j][k]);
sfree( gc->nbrs, "g:nbrs" );
sfree( gc->nbrs_x, "g:nbrs_x" );
sfree( gc->nbrs_cp, "g:nbrs_cp" );
if(gc->atoms != NULL )
sfree( gc->atoms, "g:atoms" );
}
sfree( g->cells[i][j], "g:cells[i][j]" );
}
sfree( g->cells[i], "g:cells[i]" );
}
sfree( g->cells, "g:cells" );
}
-/************ mpi buffers ********************/
- /* make the allocations based on the largest need by the three comms:
- 1- transfer an atom who has moved into other proc's domain (mpi_atom)
- 2- exchange boundary atoms (boundary_atom)
- 3- update position info for boundary atoms (mpi_rvec)
-
- the largest space by far is required for the 2nd comm operation above.
- buffers are void*, type cast to the correct pointer type to access
- the allocated buffers */
-int Allocate_MPI_Buffers( mpi_datatypes *mpi_data, int est_recv,
- neighbor_proc *my_nbrs, char *msg )
-{
- int i;
- mpi_out_data *mpi_buf;
- MPI_Comm comm;
-
- comm = mpi_data->world;
-
- /* in buffers */
- mpi_data->in1_buffer = (void*)
- scalloc( est_recv, sizeof(boundary_atom), "in1_buffer", comm );
- mpi_data->in2_buffer = (void*)
- scalloc( est_recv, sizeof(boundary_atom), "in2_buffer", comm );
-
- /* out buffers */
- for( i = 0; i < MAX_NBRS; ++i ) {
- mpi_buf = &( mpi_data->out_buffers[i] );
- /* allocate storage for the neighbor processor i */
- mpi_buf->index = (int*)
- scalloc( my_nbrs[i].est_send, sizeof(int), "mpibuf:index", comm );
- mpi_buf->out_atoms = (void*)
- scalloc( my_nbrs[i].est_send, sizeof(boundary_atom), "mpibuf:out_atoms",
- comm );
- }
-
- return SUCCESS;
-}
-
-
-void Deallocate_MPI_Buffers( mpi_datatypes *mpi_data )
-{
- int i;
- mpi_out_data *mpi_buf;
-
- sfree( mpi_data->in1_buffer, "in1_buffer" );
- sfree( mpi_data->in2_buffer, "in2_buffer" );
-
- for( i = 0; i < MAX_NBRS; ++i ) {
- mpi_buf = &( mpi_data->out_buffers[i] );
- sfree( mpi_buf->index, "mpibuf:index" );
- sfree( mpi_buf->out_atoms, "mpibuf:out_atoms" );
- }
-}
-
-
void ReAllocate( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
mpi_datatypes *mpi_data )
{
int num_bonds, est_3body, Hflag, ret;
int renbr, newsize;
reallocate_data *realloc;
reax_list *far_nbrs;
grid *g;
MPI_Comm comm;
char msg[200];
int mincap = system->mincap;
double safezone = system->safezone;
double saferzone = system->saferzone;
realloc = &(workspace->realloc);
g = &(system->my_grid);
comm = mpi_data->world;
#if defined(DEBUG)
fprintf( stderr, "p%d@reallocate: n: %d, N: %d, numH: %d\n",
system->my_rank, system->n, system->N, system->numH );
fprintf( stderr, "p%d@reallocate: local_cap: %d, total_cap: %d, Hcap: %d\n",
system->my_rank, system->local_cap, system->total_cap,
system->Hcap);
fprintf( stderr, "p%d: realloc.num_far: %d\n",
system->my_rank, realloc->num_far );
fprintf( stderr, "p%d: realloc.H: %d, realloc.Htop: %d\n",
system->my_rank, realloc->H, realloc->Htop );
fprintf( stderr, "p%d: realloc.Hbonds: %d, realloc.num_hbonds: %d\n",
system->my_rank, realloc->hbonds, realloc->num_hbonds );
fprintf( stderr, "p%d: realloc.bonds: %d, num_bonds: %d\n",
system->my_rank, realloc->bonds, realloc->num_bonds );
fprintf( stderr, "p%d: realloc.num_3body: %d\n",
system->my_rank, realloc->num_3body );
#endif
// IMPORTANT: LOOSE ZONES CHECKS ARE DISABLED FOR NOW BY &&'ing with 0!!!
int nflag = 0;
if( system->n >= DANGER_ZONE * system->local_cap ||
(0 && system->n <= LOOSE_ZONE * system->local_cap) ) {
nflag = 1;
system->local_cap = MAX( (int)(system->n * safezone), mincap );
}
int Nflag = 0;
if( system->N >= DANGER_ZONE * system->total_cap ||
(0 && system->N <= LOOSE_ZONE * system->total_cap) ) {
Nflag = 1;
system->total_cap = MAX( (int)(system->N * safezone), mincap );
}
if( Nflag ) {
/* system */
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating system and workspace -"\
"n=%d N=%d local_cap=%d total_cap=%d\n",
system->my_rank, system->n, system->N,
system->local_cap, system->total_cap );
#endif
ret = Allocate_System( system, system->local_cap, system->total_cap, msg );
if( ret != SUCCESS ) {
fprintf( stderr, "not enough space for atom_list: total_cap=%d",
system->total_cap );
fprintf( stderr, "terminating...\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
/* workspace */
DeAllocate_Workspace( control, workspace );
ret = Allocate_Workspace( system, control, workspace, system->local_cap,
system->total_cap, comm, msg );
if( ret != SUCCESS ) {
fprintf( stderr, "no space for workspace: local_cap=%d total_cap=%d",
system->local_cap, system->total_cap );
fprintf( stderr, "terminating...\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
}
renbr = (data->step - data->prev_steps) % control->reneighbor == 0;
/* far neighbors */
if( renbr ) {
far_nbrs = *lists + FAR_NBRS;
if( Nflag || realloc->num_far >= far_nbrs->num_intrs * DANGER_ZONE ) {
if( realloc->num_far > far_nbrs->num_intrs ) {
fprintf( stderr, "step%d-ran out of space on far_nbrs: top=%d, max=%d",
data->step, realloc->num_far, far_nbrs->num_intrs );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
newsize = static_cast<int>
(MAX( realloc->num_far*safezone, mincap*MIN_NBRS ));
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating far_nbrs: num_fars=%d, space=%dMB\n",
system->my_rank, (int)(realloc->num_far*SAFE_ZONE),
(newsize*sizeof(far_neighbor_data)/(1024*1024)) );
#endif
Reallocate_Neighbor_List( far_nbrs, system->total_cap, newsize, comm );
realloc->num_far = 0;
}
}
#if defined(PURE_REAX)
/* qeq coef matrix */
H = workspace->H;
if( nflag || realloc->Htop >= H->m * DANGER_ZONE ) {
if( realloc->Htop > H->m ) {
fprintf( stderr,
"step%d - ran out of space on H matrix: Htop=%d, max = %d",
data->step, realloc->Htop, H->m );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating H matrix: Htop=%d, space=%dMB\n",
system->my_rank, (int)(realloc->Htop*SAFE_ZONE),
(int)(realloc->Htop * SAFE_ZONE * sizeof(sparse_matrix_entry) /
(1024*1024)) );
#endif
newsize = static_cast<int>
(MAX( realloc->Htop*safezone, mincap*MIN_NBRS ));
Reallocate_Matrix( &(workspace->H), system->local_cap,
newsize, "H", comm );
//Deallocate_Matrix( workspace->L );
//Deallocate_Matrix( workspace->U );
workspace->L = NULL;
workspace->U = NULL;
realloc->Htop = 0;
}
#endif /*PURE_REAX*/
/* hydrogen bonds list */
if( control->hbond_cut > 0 ) {
Hflag = 0;
if( system->numH >= DANGER_ZONE * system->Hcap ||
(0 && system->numH <= LOOSE_ZONE * system->Hcap) ) {
Hflag = 1;
system->Hcap = int(MAX( system->numH * saferzone, mincap ));
}
if( Hflag || realloc->hbonds ) {
ret = Reallocate_HBonds_List( system, (*lists)+HBONDS, comm );
realloc->hbonds = 0;
#if defined(DEBUG_FOCUS)
fprintf(stderr, "p%d: reallocating hbonds: total_hbonds=%d space=%dMB\n",
system->my_rank, ret, (int)(ret*sizeof(hbond_data)/(1024*1024)));
#endif
}
}
/* bonds list */
num_bonds = est_3body = -1;
if( Nflag || realloc->bonds ){
Reallocate_Bonds_List( system, (*lists)+BONDS, &num_bonds,
&est_3body, comm );
realloc->bonds = 0;
realloc->num_3body = MAX( realloc->num_3body, est_3body );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating bonds: total_bonds=%d, space=%dMB\n",
system->my_rank, num_bonds,
(int)(num_bonds*sizeof(bond_data)/(1024*1024)) );
#endif
}
/* 3-body list */
if( realloc->num_3body > 0 ) {
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating 3body list: num_3body=%d, space=%dMB\n",
system->my_rank, realloc->num_3body,
(int)(realloc->num_3body * sizeof(three_body_interaction_data) /
(1024*1024)) );
#endif
Delete_List( (*lists)+THREE_BODIES, comm );
if( num_bonds == -1 )
num_bonds = ((*lists)+BONDS)->num_intrs;
realloc->num_3body = (int)(MAX(realloc->num_3body*safezone, MIN_3BODIES));
if( !Make_List( num_bonds, realloc->num_3body, TYP_THREE_BODY,
(*lists)+THREE_BODIES, comm ) ) {
fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
MPI_Abort( comm, CANNOT_INITIALIZE );
}
realloc->num_3body = -1;
}
-#if defined(PURE_REAX)
- /* grid */
- if( renbr && realloc->gcell_atoms > -1 ) {
-#if defined(DEBUG_FOCUS)
- fprintf(stderr, "reallocating gcell: g->max_atoms: %d\n", g->max_atoms);
-#endif
- for( i = g->native_str[0]; i < g->native_end[0]; i++ )
- for( j = g->native_str[1]; j < g->native_end[1]; j++ )
- for( k = g->native_str[2]; k < g->native_end[2]; k++ ) {
- // reallocate g->atoms
- sfree( g->cells[i][j][k].atoms, "g:atoms" );
- g->cells[i][j][k].atoms = (int*)
- scalloc( realloc->gcell_atoms, sizeof(int), "g:atoms", comm);
- }
- realloc->gcell_atoms = -1;
- }
-
- /* mpi buffers */
- // we have to be at a renbring step -
- // to ensure correct values at mpi_buffers for update_boundary_positions
- if( !renbr )
- mpi_flag = 0;
- // check whether in_buffer capacity is enough
- else if( system->max_recved >= system->est_recv * 0.90 )
- mpi_flag = 1;
- else {
- // otherwise check individual outgoing buffers
- mpi_flag = 0;
- for( p = 0; p < MAX_NBRS; ++p ) {
- nbr_pr = &( system->my_nbrs[p] );
- nbr_data = &( mpi_data->out_buffers[p] );
- if( nbr_data->cnt >= nbr_pr->est_send * 0.90 ) {
- mpi_flag = 1;
- break;
- }
- }
- }
-
- if( mpi_flag ) {
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "p%d: reallocating mpi_buf: old_recv=%d\n",
- system->my_rank, system->est_recv );
- for( p = 0; p < MAX_NBRS; ++p )
- fprintf( stderr, "p%d: nbr%d old_send=%d\n",
- system->my_rank, p, system->my_nbrs[p].est_send );
-#endif
- /* update mpi buffer estimates based on last comm */
- system->est_recv = MAX( system->max_recved*SAFER_ZONE, MIN_SEND );
- system->est_trans =
- (system->est_recv * sizeof(boundary_atom)) / sizeof(mpi_atom);
- total_send = 0;
- for( p = 0; p < MAX_NBRS; ++p ) {
- nbr_pr = &( system->my_nbrs[p] );
- nbr_data = &( mpi_data->out_buffers[p] );
- nbr_pr->est_send = MAX( nbr_data->cnt*SAFER_ZONE, MIN_SEND );
- total_send += nbr_pr->est_send;
- }
-#if defined(DEBUG_FOCUS)
- fprintf( stderr, "p%d: reallocating mpi_buf: recv=%d send=%d total=%dMB\n",
- system->my_rank, system->est_recv, total_send,
- (int)((system->est_recv+total_send)*sizeof(boundary_atom)/
- (1024*1024)));
- for( p = 0; p < MAX_NBRS; ++p )
- fprintf( stderr, "p%d: nbr%d new_send=%d\n",
- system->my_rank, p, system->my_nbrs[p].est_send );
-#endif
-
- /* reallocate mpi buffers */
- Deallocate_MPI_Buffers( mpi_data );
- ret = Allocate_MPI_Buffers( mpi_data, system->est_recv,
- system->my_nbrs, msg );
- if( ret != SUCCESS ) {
- fprintf( stderr, "%s", msg );
- fprintf( stderr, "terminating...\n" );
- MPI_Abort( comm, INSUFFICIENT_MEMORY );
- }
- }
-#endif /*PURE_REAX*/
-
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d @ step%d: reallocate done\n",
system->my_rank, data->step );
MPI_Barrier( comm );
#endif
}
diff --git a/src/USER-REAXC/reaxc_types.h b/src/USER-REAXC/reaxc_types.h
index c7580bdbd..df8bb666a 100644
--- a/src/USER-REAXC/reaxc_types.h
+++ b/src/USER-REAXC/reaxc_types.h
@@ -1,1016 +1,976 @@
/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Copyright (2010) Purdue University
Hasan Metin Aktulga, hmaktulga@lbl.gov
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
Please cite the related publication:
H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama,
"Parallel Reactive Molecular Dynamics: Numerical Methods and
Algorithmic Techniques", Parallel Computing, in press.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details:
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#ifndef __REAX_TYPES_H_
#define __REAX_TYPES_H_
#include "ctype.h"
#include "math.h"
#include "mpi.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "sys/time.h"
#include "time.h"
/************* SOME DEFS - crucial for reax_types.h *********/
//#define PURE_REAX
#define LAMMPS_REAX
//#define DEBUG
//#define DEBUG_FOCUS
//#define TEST_ENERGY
//#define TEST_FORCES
//#define CG_PERFORMANCE
#define LOG_PERFORMANCE
#define STANDARD_BOUNDARIES
//#define OLD_BOUNDARIES
//#define MIDPOINT_BOUNDARIES
#define REAX_MAX_STR 1024
#define REAX_MAX_NBRS 6
#define REAX_MAX_3BODY_PARAM 5
#define REAX_MAX_4BODY_PARAM 5
#define REAX_MAX_ATOM_TYPES 25
#define REAX_MAX_MOLECULE_SIZE 20
#define MAX_BOND 20 // same as reaxc_defs.h
/********************** TYPE DEFINITIONS ********************/
typedef int ivec[3];
typedef double real;
typedef real rvec[3];
typedef real rtensor[3][3];
typedef real rvec2[2];
typedef real rvec4[4];
#ifdef LAMMPS_SMALLSMALL
typedef int tagint;
#endif
#ifdef LAMMPS_SMALLBIG
typedef int tagint;
#endif
#ifdef LAMMPS_BIGBIG
typedef int64_t tagint;
#endif
-typedef struct {
- int step, bigN;
- real T, xi, v_xi, v_xi_old, G_xi;
- rtensor box;
-} restart_header;
-
-typedef struct {
- tagint orig_id, type;
- char name[8];
- rvec x, v;
-} restart_atom;
-
-typedef struct
-{
- tagint orig_id;
- int imprt_id;
- int type;
- int num_bonds;
- int num_hbonds;
- //int pad; // pad to 8-byte address boundary
- char name[8];
- rvec x; // position
- rvec v; // velocity
- rvec f_old; // old force
- rvec4 s, t; // for calculating q
-} mpi_atom;
-
-
-typedef struct
-{
- tagint orig_id;
- int imprt_id;
- int type;
- int num_bonds;
- int num_hbonds;
- //int pad;
- rvec x; // position
-} boundary_atom;
-
-
typedef struct
{
//int ncells;
//int *cnt_by_gcell;
int cnt;
//int *block;
int *index;
//MPI_Datatype out_dtype;
void *out_atoms;
} mpi_out_data;
typedef struct
{
MPI_Comm world;
MPI_Comm comm_mesh3D;
MPI_Datatype sys_info;
MPI_Datatype mpi_atom_type;
MPI_Datatype boundary_atom_type;
MPI_Datatype mpi_rvec, mpi_rvec2;
MPI_Datatype restart_atom_type;
MPI_Datatype header_line;
MPI_Datatype header_view;
MPI_Datatype init_desc_line;
MPI_Datatype init_desc_view;
MPI_Datatype atom_line;
MPI_Datatype atom_view;
MPI_Datatype bond_line;
MPI_Datatype bond_view;
MPI_Datatype angle_line;
MPI_Datatype angle_view;
//MPI_Request send_req1[REAX_MAX_NBRS];
//MPI_Request send_req2[REAX_MAX_NBRS];
//MPI_Status send_stat1[REAX_MAX_NBRS];
//MPI_Status send_stat2[REAX_MAX_NBRS];
//MPI_Status recv_stat1[REAX_MAX_NBRS];
//MPI_Status recv_stat2[REAX_MAX_NBRS];
mpi_out_data out_buffers[REAX_MAX_NBRS];
void *in1_buffer;
void *in2_buffer;
} mpi_datatypes;
/* Global params mapping */
/*
l[0] = p_boc1
l[1] = p_boc2
l[2] = p_coa2
l[3] = N/A
l[4] = N/A
l[5] = N/A
l[6] = p_ovun6
l[7] = N/A
l[8] = p_ovun7
l[9] = p_ovun8
l[10] = N/A
l[11] = swa
l[12] = swb
l[13] = N/A
l[14] = p_val6
l[15] = p_lp1
l[16] = p_val9
l[17] = p_val10
l[18] = N/A
l[19] = p_pen2
l[20] = p_pen3
l[21] = p_pen4
l[22] = N/A
l[23] = p_tor2
l[24] = p_tor3
l[25] = p_tor4
l[26] = N/A
l[27] = p_cot2
l[28] = p_vdW1
l[29] = v_par30
l[30] = p_coa4
l[31] = p_ovun4
l[32] = p_ovun3
l[33] = p_val8
l[34] = N/A
l[35] = N/A
l[36] = N/A
l[37] = version number
l[38] = p_coa3
*/
typedef struct
{
int n_global;
real* l;
int vdw_type;
} global_parameters;
typedef struct
{
/* Line one in field file */
char name[15]; // Two character atom name
real r_s;
real valency; // Valency of the atom
real mass; // Mass of atom
real r_vdw;
real epsilon;
real gamma;
real r_pi;
real valency_e;
real nlp_opt;
/* Line two in field file */
real alpha;
real gamma_w;
real valency_boc;
real p_ovun5;
real chi;
real eta;
int p_hbond; // 1 for H, 2 for hbonding atoms (O,S,P,N), 0 for others
/* Line three in field file */
real r_pi_pi;
real p_lp2;
real b_o_131;
real b_o_132;
real b_o_133;
/* Line four in the field file */
real p_ovun2;
real p_val3;
real valency_val;
real p_val5;
real rcore2;
real ecore2;
real acore2;
/* Line five in the ffield file, only for lgvdw yes */
real lgcij;
real lgre;
} single_body_parameters;
/* Two Body Parameters */
typedef struct {
/* Bond Order parameters */
real p_bo1,p_bo2,p_bo3,p_bo4,p_bo5,p_bo6;
real r_s, r_p, r_pp; // r_o distances in BO formula
real p_boc3, p_boc4, p_boc5;
/* Bond Energy parameters */
real p_be1, p_be2;
real De_s, De_p, De_pp;
/* Over/Under coordination parameters */
real p_ovun1;
/* Van der Waal interaction parameters */
real D;
real alpha;
real r_vdW;
real gamma_w;
real rcore, ecore, acore;
real lgcij, lgre;
/* electrostatic parameters */
real gamma; // note: this parameter is gamma^-3 and not gamma.
real v13cor, ovc;
} two_body_parameters;
/* 3-body parameters */
typedef struct {
/* valence angle */
real theta_00;
real p_val1, p_val2, p_val4, p_val7;
/* penalty */
real p_pen1;
/* 3-body conjugation */
real p_coa1;
} three_body_parameters;
typedef struct{
int cnt;
three_body_parameters prm[REAX_MAX_3BODY_PARAM];
} three_body_header;
/* hydrogen-bond parameters */
typedef struct{
real r0_hb, p_hb1, p_hb2, p_hb3;
} hbond_parameters;
/* 4-body parameters */
typedef struct {
real V1, V2, V3;
/* torsion angle */
real p_tor1;
/* 4-body conjugation */
real p_cot1;
} four_body_parameters;
typedef struct
{
int cnt;
four_body_parameters prm[REAX_MAX_4BODY_PARAM];
} four_body_header;
typedef struct
{
int num_atom_types;
global_parameters gp;
single_body_parameters *sbp;
two_body_parameters **tbp;
three_body_header ***thbp;
hbond_parameters ***hbp;
four_body_header ****fbp;
} reax_interaction;
struct _reax_atom
{
tagint orig_id;
int imprt_id;
int type;
char name[8];
rvec x; // position
rvec v; // velocity
rvec f; // force
rvec f_old;
real q; // charge
rvec4 s; // they take part in
rvec4 t; // computing q
int Hindex;
int num_bonds;
int num_hbonds;
int renumber;
int numbonds; // true number of bonds around atoms
int nbr_id[MAX_BOND]; // ids of neighbors around atoms
double nbr_bo[MAX_BOND]; // BO values of bond between i and nbr
double sum_bo, no_lp; // sum of BO values and no. of lone pairs
};
typedef _reax_atom reax_atom;
typedef struct
{
real V;
rvec min, max, box_norms;
rtensor box, box_inv;
rtensor trans, trans_inv;
rtensor g;
} simulation_box;
struct grid_cell
{
real cutoff;
rvec min, max;
ivec rel_box;
int mark;
int type;
int str;
int end;
int top;
int* atoms;
struct grid_cell** nbrs;
ivec* nbrs_x;
rvec* nbrs_cp;
};
typedef struct grid_cell grid_cell;
typedef struct
{
int total, max_atoms, max_nbrs;
ivec ncells;
rvec cell_len;
rvec inv_len;
ivec bond_span;
ivec nonb_span;
ivec vlist_span;
ivec native_cells;
ivec native_str;
ivec native_end;
real ghost_cut;
ivec ghost_span;
ivec ghost_nonb_span;
ivec ghost_hbond_span;
ivec ghost_bond_span;
grid_cell*** cells;
ivec *order;
} grid;
typedef struct
{
int rank;
int est_send, est_recv;
int atoms_str, atoms_cnt;
ivec rltv, prdc;
rvec bndry_min, bndry_max;
int send_type;
int recv_type;
ivec str_send;
ivec end_send;
ivec str_recv;
ivec end_recv;
} neighbor_proc;
typedef struct
{
int N;
int exc_gcells;
int exc_atoms;
} bound_estimate;
typedef struct
{
real ghost_nonb;
real ghost_hbond;
real ghost_bond;
real ghost_cutoff;
} boundary_cutoff;
using LAMMPS_NS::Pair;
struct _reax_system
{
reax_interaction reax_param;
int n, N, bigN, numH;
int local_cap, total_cap, gcell_cap, Hcap;
int est_recv, est_trans, max_recved;
int wsize, my_rank, num_nbrs;
ivec my_coords;
neighbor_proc my_nbrs[REAX_MAX_NBRS];
int *global_offset;
simulation_box big_box, my_box, my_ext_box;
grid my_grid;
boundary_cutoff bndry_cuts;
reax_atom *my_atoms;
class Pair *pair_ptr;
int my_bonds;
int mincap;
real safezone, saferzone;
};
typedef _reax_system reax_system;
/* system control parameters */
typedef struct
{
char sim_name[REAX_MAX_STR];
int nprocs;
ivec procs_by_dim;
/* ensemble values:
0 : NVE
1 : bNVT (Berendsen)
2 : nhNVT (Nose-Hoover)
3 : sNPT (Parrinello-Rehman-Nose-Hoover) semiisotropic
4 : iNPT (Parrinello-Rehman-Nose-Hoover) isotropic
5 : NPT (Parrinello-Rehman-Nose-Hoover) Anisotropic*/
int ensemble;
int nsteps;
real dt;
int geo_format;
int restart;
int restrict_bonds;
int remove_CoM_vel;
int random_vel;
int reposition_atoms;
int reneighbor;
real vlist_cut;
real bond_cut;
real nonb_cut, nonb_low;
real hbond_cut;
real user_ghost_cut;
real bg_cut;
real bo_cut;
real thb_cut;
real thb_cutsq;
int tabulate;
int qeq_freq;
real q_err;
int refactor;
real droptol;
real T_init, T_final, T;
real Tau_T;
int T_mode;
real T_rate, T_freq;
int virial;
rvec P, Tau_P, Tau_PT;
int press_mode;
real compressibility;
int molecular_analysis;
int num_ignored;
int ignore[REAX_MAX_ATOM_TYPES];
int dipole_anal;
int freq_dipole_anal;
int diffusion_coef;
int freq_diffusion_coef;
int restrict_type;
int lgflag;
} control_params;
typedef struct
{
real T;
real xi;
real v_xi;
real v_xi_old;
real G_xi;
} thermostat;
typedef struct
{
real P;
real eps;
real v_eps;
real v_eps_old;
real a_eps;
} isotropic_barostat;
typedef struct
{
rtensor P;
real P_scalar;
real eps;
real v_eps;
real v_eps_old;
real a_eps;
rtensor h0;
rtensor v_g0;
rtensor v_g0_old;
rtensor a_g0;
} flexible_barostat;
typedef struct
{
real start;
real end;
real elapsed;
real total;
real comm;
real nbrs;
real init_forces;
real bonded;
real nonb;
real qEq;
int s_matvecs;
int t_matvecs;
} reax_timing;
typedef struct
{
real e_tot;
real e_kin; // Total kinetic energy
real e_pot;
real e_bond; // Total bond energy
real e_ov; // Total over coordination
real e_un; // Total under coordination energy
real e_lp; // Total under coordination energy
real e_ang; // Total valance angle energy
real e_pen; // Total penalty energy
real e_coa; // Total three body conjgation energy
real e_hb; // Total Hydrogen bond energy
real e_tor; // Total torsional energy
real e_con; // Total four body conjugation energy
real e_vdW; // Total van der Waals energy
real e_ele; // Total electrostatics energy
real e_pol; // Polarization energy
} energy_data;
typedef struct
{
int step;
int prev_steps;
real time;
real M; // Total Mass
real inv_M; // 1 / Total Mass
rvec xcm; // Center of mass
rvec vcm; // Center of mass velocity
rvec fcm; // Center of mass force
rvec amcm; // Angular momentum of CoM
rvec avcm; // Angular velocity of CoM
real etran_cm; // Translational kinetic energy of CoM
real erot_cm; // Rotational kinetic energy of CoM
rtensor kinetic; // Kinetic energy tensor
rtensor virial; // Hydrodynamic virial
energy_data my_en;
energy_data sys_en;
real N_f; //Number of degrees of freedom
rvec t_scale;
rtensor p_scale;
thermostat therm; // Used in Nose_Hoover method
isotropic_barostat iso_bar;
flexible_barostat flex_bar;
real inv_W;
real kin_press;
rvec int_press;
rvec my_ext_press;
rvec ext_press;
rvec tot_press;
reax_timing timing;
} simulation_data;
typedef struct{
int thb;
int pthb; // pointer to the third body on the central atom's nbrlist
real theta, cos_theta;
rvec dcos_di, dcos_dj, dcos_dk;
} three_body_interaction_data;
typedef struct {
int nbr;
ivec rel_box;
real d;
rvec dvec;
} far_neighbor_data;
typedef struct {
int nbr;
int scl;
far_neighbor_data *ptr;
} hbond_data;
typedef struct{
int wrt;
rvec dVal;
} dDelta_data;
typedef struct{
int wrt;
rvec dBO, dBOpi, dBOpi2;
} dbond_data;
typedef struct{
real BO, BO_s, BO_pi, BO_pi2;
real Cdbo, Cdbopi, Cdbopi2;
real C1dbo, C2dbo, C3dbo;
real C1dbopi, C2dbopi, C3dbopi, C4dbopi;
real C1dbopi2, C2dbopi2, C3dbopi2, C4dbopi2;
rvec dBOp, dln_BOp_s, dln_BOp_pi, dln_BOp_pi2;
} bond_order_data;
typedef struct {
int nbr;
int sym_index;
int dbond_index;
ivec rel_box;
// rvec ext_factor;
real d;
rvec dvec;
bond_order_data bo_data;
} bond_data;
typedef struct {
int j;
real val;
} sparse_matrix_entry;
typedef struct {
int cap, n, m;
int *start, *end;
sparse_matrix_entry *entries;
} sparse_matrix;
typedef struct {
int num_far;
int H, Htop;
int hbonds, num_hbonds;
int bonds, num_bonds;
int num_3body;
int gcell_atoms;
} reallocate_data;
typedef struct
{
int allocated;
/* communication storage */
real *tmp_dbl[REAX_MAX_NBRS];
rvec *tmp_rvec[REAX_MAX_NBRS];
rvec2 *tmp_rvec2[REAX_MAX_NBRS];
int *within_bond_box;
/* bond order related storage */
real *total_bond_order;
real *Deltap, *Deltap_boc;
real *Delta, *Delta_lp, *Delta_lp_temp, *Delta_e, *Delta_boc, *Delta_val;
real *dDelta_lp, *dDelta_lp_temp;
real *nlp, *nlp_temp, *Clp, *vlpex;
rvec *dDeltap_self;
int *bond_mark, *done_after;
/* QEq storage */
sparse_matrix *H, *L, *U;
real *Hdia_inv, *b_s, *b_t, *b_prc, *b_prm, *s, *t;
real *droptol;
rvec2 *b, *x;
/* GMRES storage */
real *y, *z, *g;
real *hc, *hs;
real **h, **v;
/* CG storage */
real *r, *d, *q, *p;
rvec2 *r2, *d2, *q2, *p2;
/* Taper */
real Tap[8]; //Tap7, Tap6, Tap5, Tap4, Tap3, Tap2, Tap1, Tap0;
/* storage for analysis */
int *mark, *old_mark;
rvec *x_old;
/* storage space for bond restrictions */
int *restricted;
int **restricted_list;
/* integrator */
rvec *v_const;
/* force calculations */
real *CdDelta; // coefficient of dDelta
rvec *f;
#ifdef TEST_FORCES
rvec *f_ele;
rvec *f_vdw;
rvec *f_bo;
rvec *f_be;
rvec *f_lp;
rvec *f_ov;
rvec *f_un;
rvec *f_ang;
rvec *f_coa;
rvec *f_pen;
rvec *f_hb;
rvec *f_tor;
rvec *f_con;
rvec *f_tot;
rvec *dDelta; // calculated on the fly in bond_orders.c together with bo'
int *rcounts;
int *displs;
int *id_all;
rvec *f_all;
#endif
reallocate_data realloc;
//int *num_bonds;
/* hydrogen bonds */
//int num_H, Hcap;
//int *Hindex;
//int *num_hbonds;
//int *hash;
//int *rev_hash;
} storage;
typedef union
{
void *v;
three_body_interaction_data *three_body_list;
bond_data *bond_list;
dbond_data *dbo_list;
dDelta_data *dDelta_list;
far_neighbor_data *far_nbr_list;
hbond_data *hbond_list;
} list_type;
struct _reax_list
{
int allocated;
int n;
int num_intrs;
int *index;
int *end_index;
int type;
list_type select;
};
typedef _reax_list reax_list;
typedef struct
{
#if defined(PURE_REAX)
MPI_File trj;
#endif
FILE *strj;
int trj_offset;
int atom_line_len;
int bond_line_len;
int angle_line_len;
int write_atoms;
int write_bonds;
int write_angles;
char *line;
int buffer_len;
char *buffer;
FILE *out;
FILE *pot;
FILE *log;
FILE *mol, *ign;
FILE *dpl;
FILE *drft;
FILE *pdb;
FILE *prs;
int write_steps;
int traj_compress;
int traj_method;
char traj_title[81];
int atom_info;
int bond_info;
int angle_info;
int restart_format;
int restart_freq;
int debug_level;
int energy_update_freq;
#ifdef TEST_ENERGY
FILE *ebond;
FILE *elp, *eov, *eun;
FILE *eval, *epen, *ecoa;
FILE *ehb;
FILE *etor, *econ;
FILE *evdw, *ecou;
#endif
#ifdef TEST_FORCES
FILE *fbo, *fdbo;
FILE *fbond;
FILE *flp, *fov, *fun;
FILE *fang, *fcoa, *fpen;
FILE *fhb;
FILE *ftor, *fcon;
FILE *fvdw, *fele;
FILE *ftot, *fcomp;
#endif
#if defined(TEST_ENERGY) || defined(TEST_FORCES)
FILE *flist; // far neighbor list
FILE *blist; // bond list
FILE *nlist; // near neighbor list
#endif
} output_controls;
typedef struct
{
int atom_count;
int atom_list[REAX_MAX_MOLECULE_SIZE];
int mtypes[REAX_MAX_ATOM_TYPES];
} molecule;
typedef struct
{
real H;
real e_vdW, CEvd;
real e_ele, CEclmb;
} LR_data;
typedef struct
{
real a, b, c, d;
} cubic_spline_coef;
typedef struct
{
real xmin, xmax;
int n;
real dx, inv_dx;
real a;
real m;
real c;
LR_data *y;
cubic_spline_coef *H;
cubic_spline_coef *vdW, *CEvd;
cubic_spline_coef *ele, *CEclmb;
} LR_lookup_table;
extern LR_lookup_table **LR;
/* function pointer defs */
typedef void (*evolve_function)(reax_system*, control_params*,
simulation_data*, storage*, reax_list**,
output_controls*, mpi_datatypes* );
#if defined(PURE_REAX)
evolve_function Evolve;
#endif
typedef void (*interaction_function) (reax_system*, control_params*,
simulation_data*, storage*,
reax_list**, output_controls*);
typedef void (*print_interaction)(reax_system*, control_params*,
simulation_data*, storage*,
reax_list**, output_controls*);
typedef real (*lookup_function)(real);
typedef void (*message_sorter) (reax_system*, int, int, int, mpi_out_data*);
typedef void (*unpacker) ( reax_system*, int, void*, int, neighbor_proc*, int );
typedef void (*dist_packer) (void*, mpi_out_data*);
typedef void (*coll_unpacker) (void*, void*, mpi_out_data*);
#endif