jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1367 from mkanski/reax_better_errors 

Removing calls to MPI_Abort from reax/c + a fix for a memory leak
This commit is contained in:
Axel Kohlmeyer 2019-04-04 12:47:02 -04:00 committed by GitHub
parent f73ccc3025 6ff1fee6d6
commit 968ac3d8d7
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
25 changed files with 576 additions and 610 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
src/KOKKOS/pair_reaxc_kokkos.cpp
View File

@ -447,23 +447,23 @@ int PairReaxCKokkos<DeviceType>::Init_Lookup_Tables()
num_atom_types = atom->ntypes;
dr = control->nonb_cut / control->tabulate;
h = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:h", world );
smalloc( control->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:h");
fh = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fh", world );
smalloc( control->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fh");
fvdw = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fvdw", world );
smalloc( control->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fvdw");
fCEvd = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fCEvd", world );
smalloc( control->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fCEvd");
fele = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fele", world );
smalloc( control->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fele");
fCEclmb = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fCEclmb", world );
smalloc( control->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fCEclmb");
LR = (LR_lookup_table**)
scalloc( num_atom_types+1, sizeof(LR_lookup_table*), "lookup:LR", world );
scalloc( control->error_ptr, num_atom_types+1, sizeof(LR_lookup_table*), "lookup:LR");
for( i = 0; i < num_atom_types+1; ++i )
LR[i] = (LR_lookup_table*)
scalloc( num_atom_types+1, sizeof(LR_lookup_table), "lookup:LR[i]", world );
scalloc( control->error_ptr, num_atom_types+1, sizeof(LR_lookup_table), "lookup:LR[i]");
for( i = 1; i <= num_atom_types; ++i ) {
for( j = i; j <= num_atom_types; ++j ) {
@ -473,22 +473,18 @@ int PairReaxCKokkos<DeviceType>::Init_Lookup_Tables()
LR[i][j].dx = dr;
LR[i][j].inv_dx = control->tabulate / control->nonb_cut;
LR[i][j].y = (LR_data*)
smalloc( LR[i][j].n * sizeof(LR_data), "lookup:LR[i,j].y", world );
smalloc( control->error_ptr, LR[i][j].n * sizeof(LR_data), "lookup:LR[i,j].y");
LR[i][j].H = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].H" ,
world );
smalloc( control->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].H");
LR[i][j].vdW = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].vdW",
world);
smalloc( control->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].vdW");
LR[i][j].CEvd = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].CEvd",
world);
smalloc( control->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].CEvd");
LR[i][j].ele = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].ele",
world );
smalloc( control->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].ele");
LR[i][j].CEclmb = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),
"lookup:LR[i,j].CEclmb", world );
smalloc( control->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),
"lookup:LR[i,j].CEclmb");
for( r = 1; r <= control->tabulate; ++r ) {
LR_vdW_Coulomb(i, j, r * dr, &(LR[i][j].y[r]) );
@ -512,24 +508,20 @@ int PairReaxCKokkos<DeviceType>::Init_Lookup_Tables()
vlast_vdw = fCEvd[r-1];
vlast_ele = fele[r-1];
Natural_Cubic_Spline( &h[1], &fh[1],
&(LR[i][j].H[1]), control->tabulate+1, world );
Natural_Cubic_Spline( control->error_ptr, &h[1], &fh[1],
&(LR[i][j].H[1]), control->tabulate+1 );
Complete_Cubic_Spline( &h[1], &fvdw[1], v0_vdw, vlast_vdw,
&(LR[i][j].vdW[1]), control->tabulate+1,
world );
Complete_Cubic_Spline( control->error_ptr, &h[1], &fvdw[1], v0_vdw, vlast_vdw,
&(LR[i][j].vdW[1]), control->tabulate+1 );
Natural_Cubic_Spline( &h[1], &fCEvd[1],
&(LR[i][j].CEvd[1]), control->tabulate+1,
world );
Natural_Cubic_Spline( control->error_ptr, &h[1], &fCEvd[1],
&(LR[i][j].CEvd[1]), control->tabulate+1 );
Complete_Cubic_Spline( &h[1], &fele[1], v0_ele, vlast_ele,
&(LR[i][j].ele[1]), control->tabulate+1,
world );
Complete_Cubic_Spline( control->error_ptr, &h[1], &fele[1], v0_ele, vlast_ele,
&(LR[i][j].ele[1]), control->tabulate+1 );
Natural_Cubic_Spline( &h[1], &fCEclmb[1],
&(LR[i][j].CEclmb[1]), control->tabulate+1,
world );
Natural_Cubic_Spline( control->error_ptr, &h[1], &fCEclmb[1],
&(LR[i][j].CEclmb[1]), control->tabulate+1 );
}
}
free(h);
@ -555,16 +547,16 @@ void PairReaxCKokkos<DeviceType>::Deallocate_Lookup_Tables()
for( i = 0; i <= ntypes; ++i ) {
for( j = i; j <= ntypes; ++j )
if (LR[i][j].n) {
sfree( LR[i][j].y, "LR[i,j].y" );
sfree( LR[i][j].H, "LR[i,j].H" );
sfree( LR[i][j].vdW, "LR[i,j].vdW" );
sfree( LR[i][j].CEvd, "LR[i,j].CEvd" );
sfree( LR[i][j].ele, "LR[i,j].ele" );
sfree( LR[i][j].CEclmb, "LR[i,j].CEclmb" );
sfree( control->error_ptr, LR[i][j].y, "LR[i,j].y" );
sfree( control->error_ptr, LR[i][j].H, "LR[i,j].H" );
sfree( control->error_ptr, LR[i][j].vdW, "LR[i,j].vdW" );
sfree( control->error_ptr, LR[i][j].CEvd, "LR[i,j].CEvd" );
sfree( control->error_ptr, LR[i][j].ele, "LR[i,j].ele" );
sfree( control->error_ptr, LR[i][j].CEclmb, "LR[i,j].CEclmb" );
}
sfree( LR[i], "LR[i]" );
sfree( control->error_ptr, LR[i], "LR[i]" );
}
sfree( LR, "LR" );
sfree( control->error_ptr, LR, "LR" );
}
/* ---------------------------------------------------------------------- */

10
src/USER-OMP/pair_reaxc_omp.cpp
View File

@ -113,7 +113,7 @@ PairReaxCOMP::~PairReaxCOMP()
if (setup_flag) {
reax_list * bonds = lists+BONDS;
for (int i=0; i<bonds->num_intrs; ++i)
sfree(bonds->select.bond_list[i].bo_data.CdboReduction, "CdboReduction");
sfree(error, bonds->select.bond_list[i].bo_data.CdboReduction, "CdboReduction");
}
memory->destroy(num_nbrs_offset);
@ -209,7 +209,7 @@ void PairReaxCOMP::compute(int eflag, int vflag)
setup();
Reset( system, control, data, workspace, &lists, world );
Reset( system, control, data, workspace, &lists );
// Why not update workspace like in MPI-only code?
// Using the MPI-only way messes up the hb energy
@ -410,12 +410,12 @@ void PairReaxCOMP::setup( )
// initialize my data structures
PreAllocate_Space( system, control, workspace, world );
PreAllocate_Space( system, control, workspace );
write_reax_atoms();
int num_nbrs = estimate_reax_lists();
if(!Make_List(system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR,
lists+FAR_NBRS, world))
lists+FAR_NBRS))
error->all(FLERR,"Pair reax/c problem in far neighbor list");
write_reax_lists();
@ -445,7 +445,7 @@ void PairReaxCOMP::setup( )
// check if I need to shrink/extend my data-structs
ReAllocate( system, control, data, workspace, &lists, mpi_data );
ReAllocate( system, control, data, workspace, &lists );
}
}

12
src/USER-OMP/reaxc_forces_omp.cpp
View File

@ -289,9 +289,10 @@ void Validate_ListsOMP( reax_system *system, storage * /*workspace */, reax_list
else comp = bonds->num_intrs;
if (End_Index(i, bonds) > comp) {
fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
step, i, End_Index(i,bonds), comp );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
char errmsg[256];
snprintf(errmsg, 256, "step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
step, i, End_Index(i,bonds), comp );
system->error_ptr->one(FLERR,errmsg);
}
}
}
@ -315,9 +316,10 @@ void Validate_ListsOMP( reax_system *system, storage * /*workspace */, reax_list
else comp = hbonds->num_intrs;
if (End_Index(Hindex, hbonds) > comp) {
fprintf(stderr,"step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
char errmsg[256];
snprintf(errmsg, 256, "step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
step, Hindex, End_Index(Hindex,hbonds), comp );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
system->error_ptr->one(FLERR, errmsg);
}
}
}

73
src/USER-OMP/reaxc_init_md_omp.cpp
View File

@ -43,7 +43,7 @@
extern int Init_MPI_Datatypes(reax_system*, storage*, mpi_datatypes*, MPI_Comm, char*);
extern int Init_System(reax_system*, control_params*, char*);
extern int Init_Simulation_Data(reax_system*, control_params*, simulation_data*, char*);
extern int Init_Workspace(reax_system*, control_params*, storage*, MPI_Comm, char*);
extern int Init_Workspace(reax_system*, control_params*, storage*, char*);
/* ---------------------------------------------------------------------- */
@ -63,7 +63,7 @@ int Init_ListsOMP( reax_system *system, control_params *control,
bond_top = (int*) calloc( system->total_cap, sizeof(int) );
hb_top = (int*) calloc( system->local_cap, sizeof(int) );
Estimate_Storages( system, control, lists,
&Htop, hb_top, bond_top, &num_3body, comm );
&Htop, hb_top, bond_top, &num_3body );
if (control->hbond_cut > 0) {
/* init H indexes */
@ -75,9 +75,8 @@ int Init_ListsOMP( reax_system *system, control_params *control,
total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
*lists+HBONDS, comm ) ) {
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
*lists+HBONDS ) ) {
system->error_ptr->one( FLERR, "Not enough space for hbonds list. Terminating!" );
}
}
@ -89,9 +88,8 @@ int Init_ListsOMP( reax_system *system, control_params *control,
bond_cap = (int)(MAX( total_bonds*safezone, mincap*MIN_BONDS ));
if( !Make_List( system->total_cap, bond_cap, TYP_BOND,
*lists+BONDS, comm ) ) {
fprintf( stderr, "not enough space for bonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
*lists+BONDS ) ) {
system->error_ptr->one( FLERR, "Not enough space for bonds list. Terminating!\n" );
}
int nthreads = control->nthreads;
@ -99,15 +97,14 @@ int Init_ListsOMP( reax_system *system, control_params *control,
for (i = 0; i < bonds->num_intrs; ++i)
bonds->select.bond_list[i].bo_data.CdboReduction =
(double*) smalloc(sizeof(double)*nthreads, "CdboReduction", comm);
(double*) smalloc(system->error_ptr, sizeof(double)*nthreads, "CdboReduction");
/* 3bodies list */
cap_3body = (int)(MAX( num_3body*safezone, MIN_3BODIES ));
if( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
*lists+THREE_BODIES, comm ) ){
*lists+THREE_BODIES ) ){
fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
system->error_ptr->one( FLERR, "Problem in initializing angles list. Terminating!" );
}
free( hb_top );
@ -125,60 +122,50 @@ void InitializeOMP( reax_system *system, control_params *control,
mpi_datatypes *mpi_data, MPI_Comm comm )
{
char msg[MAX_STR];
char errmsg[512];
if (Init_MPI_Datatypes(system, workspace, mpi_data, comm, msg) == FAILURE) {
fprintf( stderr, "p%d: init_mpi_datatypes: could not create datatypes\n",
system->my_rank );
fprintf( stderr, "p%d: mpi_data couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
system->error_ptr->one( FLERR, "init_mpi_datatypes: could not create datatypes. "
"Mpi_data couldn't be initialized! Terminating.");
}
if (Init_System(system, control, msg) == FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: system could not be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
snprintf( errmsg, 512, "Error on: %s. "
"System could not be initialized! Terminating.", msg );
system->error_ptr->one(FLERR, errmsg);
}
if (Init_Simulation_Data( system, control, data, msg ) == FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: sim_data couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
snprintf( errmsg, 512, "Error on: %s. "
"Sim_data couldn't be initialized! Terminating.", msg );
system->error_ptr->one(FLERR, errmsg);
}
if (Init_Workspace( system, control, workspace, mpi_data->world, msg ) ==
if (Init_Workspace( system, control, workspace, msg ) ==
FAILURE) {
fprintf( stderr, "p%d:init_workspace: not enough memory\n",
system->my_rank );
fprintf( stderr, "p%d:workspace couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
system->error_ptr->one(FLERR, "init_workspace: not enough memory. "
"Workspace couldn't be initialized! Terminating.");
}
if (Init_ListsOMP( system, control, data, workspace, lists, mpi_data, msg ) ==
FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: system could not be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
snprintf( errmsg, 512, "Error on: %s. "
"System could not be initialized! Terminating.", msg );
system->error_ptr->one(FLERR, errmsg);
}
if (Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: could not open output files! terminating...\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
snprintf( errmsg, 512, "Error on: %s"
"Could not open output files! Terminating.", msg );
system->error_ptr->one(FLERR, errmsg);
}
if (control->tabulate) {
if (Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: couldn't create lookup table! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
snprintf( errmsg, 512, "Error on: %s."
" Couldn't create lookup table! Terminating.", msg );
system->error_ptr->one(FLERR, errmsg);
}
}

19
src/USER-OMP/reaxc_valence_angles_omp.cpp
View File

@ -237,12 +237,12 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
// Confirm that thb_intrs->num_intrs / nthreads is enough to hold all angles from a single atom
if(my_offset >= (tid+1)*per_thread) {
int me;
MPI_Comm_rank(MPI_COMM_WORLD,&me);
fprintf( stderr, "step%d-ran out of space on angle_list on proc %i for atom %i:", data->step, me, j);
fprintf( stderr, " nthreads= %d, tid=%d, my_offset=%d, per_thread=%d\n", nthreads, tid, my_offset, per_thread);
fprintf( stderr, " num_intrs= %i N= %i\n",thb_intrs->num_intrs , system->N);
MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
char errmsg[512];
snprintf( errmsg, 512, "step%d-ran out of space on angle_list for atom %i:\n"
" nthreads= %d, tid=%d, my_offset=%d, per_thread=%d\n"
" num_intrs= %i N= %i\n"
, data->step, j, nthreads, tid, my_offset, per_thread,thb_intrs->num_intrs , system->N);
control->error_ptr->one(FLERR, errmsg);
}
// Number of angles owned by this atom
@ -601,9 +601,10 @@ void Valence_AnglesOMP( reax_system *system, control_params *control,
if (num_thb_intrs >= thb_intrs->num_intrs * DANGER_ZONE) {
workspace->realloc.num_3body = num_thb_intrs * TWICE;
if (num_thb_intrs > thb_intrs->num_intrs) {
fprintf( stderr, "step%d-ran out of space on angle_list: top=%d, max=%d",
data->step, num_thb_intrs, thb_intrs->num_intrs );
MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
char errmsg[128];
snprintf(errmsg, 128, "step%d-ran out of space on angle_list: top=%d, max=%d",
data->step, num_thb_intrs, thb_intrs->num_intrs);
control->error_ptr->one(FLERR, errmsg);
}
}

25
src/USER-REAXC/pair_reaxc.cpp
View File

@ -91,6 +91,7 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp)
memory->smalloc(sizeof(mpi_datatypes),"reax:mpi");
MPI_Comm_rank(world,&system->my_rank);
control->me = system->my_rank;
system->my_coords[0] = 0;
system->my_coords[1] = 0;
@ -108,6 +109,8 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp)
system->bndry_cuts.ghost_cutoff = 0;
system->my_atoms = NULL;
system->pair_ptr = this;
system->error_ptr = error;
control->error_ptr = error;
system->omp_active = 0;
@ -139,10 +142,10 @@ PairReaxC::~PairReaxC()
if (control->tabulate ) Deallocate_Lookup_Tables( system);
if (control->hbond_cut > 0 ) Delete_List( lists+HBONDS, world);
Delete_List( lists+BONDS, world );
Delete_List( lists+THREE_BODIES, world );
Delete_List( lists+FAR_NBRS, world );
if (control->hbond_cut > 0 ) Delete_List( lists+HBONDS );
Delete_List( lists+BONDS );
Delete_List( lists+THREE_BODIES );
Delete_List( lists+FAR_NBRS );
DeAllocate_Workspace( control, workspace );
DeAllocate_System( system );
@ -394,7 +397,8 @@ void PairReaxC::init_style( )
"increased neighbor list skin.");
for( int i = 0; i < LIST_N; ++i )
lists[i].allocated = 0;
if (lists[i].allocated != 1)
lists[i].allocated = 0;
if (fix_reax == NULL) {
char **fixarg = new char*[3];
@ -436,13 +440,14 @@ void PairReaxC::setup( )
// initialize my data structures
PreAllocate_Space( system, control, workspace, world );
PreAllocate_Space( system, control, workspace );
write_reax_atoms();
int num_nbrs = estimate_reax_lists();
if(!Make_List(system->total_cap, num_nbrs, TYP_FAR_NEIGHBOR,
lists+FAR_NBRS, world))
error->all(FLERR,"Pair reax/c problem in far neighbor list");
lists+FAR_NBRS))
error->one(FLERR,"Pair reax/c problem in far neighbor list");
(lists+FAR_NBRS)->error_ptr=error;
write_reax_lists();
Initialize( system, control, data, workspace, &lists, out_control,
@ -465,7 +470,7 @@ void PairReaxC::setup( )
// check if I need to shrink/extend my data-structs
ReAllocate( system, control, data, workspace, &lists, mpi_data );
ReAllocate( system, control, data, workspace, &lists );
}
bigint local_ngroup = list->inum;
@ -516,7 +521,7 @@ void PairReaxC::compute(int eflag, int vflag)
setup();
Reset( system, control, data, workspace, &lists, world );
Reset( system, control, data, workspace, &lists );
workspace->realloc.num_far = write_reax_lists();
// timing for filling in the reax lists
if (comm->me == 0) {

341
src/USER-REAXC/reaxc_allocate.cpp
View File

@ -35,12 +35,14 @@
#include <omp.h>
#endif
#include "error.h"
/* 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 )
storage * workspace )
{
int mincap = system->mincap;
double safezone = system->safezone;
@ -51,7 +53,7 @@ int PreAllocate_Space( reax_system *system, control_params * /*control*/,
system->total_cap = MAX( (int)(system->N * safezone), mincap );
system->my_atoms = (reax_atom*)
scalloc( system->total_cap, sizeof(reax_atom), "my_atoms", comm );
scalloc(system->error_ptr, system->total_cap, sizeof(reax_atom), "my_atoms");
// Nullify some arrays only used in omp styles
// Should be safe to do here since called in pair->setup();
@ -86,39 +88,39 @@ void DeAllocate_System( reax_system *system )
int ntypes;
reax_interaction *ff_params;
// dealloocate the atom list
sfree( system->my_atoms, "system->my_atoms" );
// deallocate the atom list
sfree(system->error_ptr, 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" );
sfree(system->error_ptr, 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(system->error_ptr, 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(system->error_ptr, ff_params->fbp[i][j], "ff:fbp[i,j]" );
sfree(system->error_ptr, ff_params->thbp[i][j], "ff:thbp[i,j]" );
sfree(system->error_ptr, 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(system->error_ptr, ff_params->fbp[i], "ff:fbp[i]" );
sfree(system->error_ptr, ff_params->thbp[i], "ff:thbp[i]" );
sfree(system->error_ptr, ff_params->hbp[i], "ff:hbp[i]" );
sfree(system->error_ptr, 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" );
sfree(system->error_ptr, ff_params->fbp, "ff:fbp" );
sfree(system->error_ptr, ff_params->thbp, "ff:thbp" );
sfree(system->error_ptr, ff_params->hbp, "ff:hbp" );
sfree(system->error_ptr, ff_params->tbp, "ff:tbp" );
sfree(system->error_ptr, ff_params->sbp, "ff:sbp" );
}
/************* workspace *************/
void DeAllocate_Workspace( control_params * /*control*/, storage *workspace )
void DeAllocate_Workspace( control_params * control, storage *workspace )
{
int i;
@ -129,86 +131,86 @@ void DeAllocate_Workspace( control_params * /*control*/, storage *workspace )
/* 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]" );
sfree(control->error_ptr, workspace->tmp_dbl[i], "tmp_dbl[i]" );
sfree(control->error_ptr, workspace->tmp_rvec[i], "tmp_rvec[i]" );
sfree(control->error_ptr, 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" );
sfree(control->error_ptr, workspace->within_bond_box, "skin" );
sfree(control->error_ptr, workspace->total_bond_order, "total_bo" );
sfree(control->error_ptr, workspace->Deltap, "Deltap" );
sfree(control->error_ptr, workspace->Deltap_boc, "Deltap_boc" );
sfree(control->error_ptr, workspace->dDeltap_self, "dDeltap_self" );
sfree(control->error_ptr, workspace->Delta, "Delta" );
sfree(control->error_ptr, workspace->Delta_lp, "Delta_lp" );
sfree(control->error_ptr, workspace->Delta_lp_temp, "Delta_lp_temp" );
sfree(control->error_ptr, workspace->dDelta_lp, "dDelta_lp" );
sfree(control->error_ptr, workspace->dDelta_lp_temp, "dDelta_lp_temp" );
sfree(control->error_ptr, workspace->Delta_e, "Delta_e" );
sfree(control->error_ptr, workspace->Delta_boc, "Delta_boc" );
sfree(control->error_ptr, workspace->Delta_val, "Delta_val" );
sfree(control->error_ptr, workspace->nlp, "nlp" );
sfree(control->error_ptr, workspace->nlp_temp, "nlp_temp" );
sfree(control->error_ptr, workspace->Clp, "Clp" );
sfree(control->error_ptr, workspace->vlpex, "vlpex" );
sfree(control->error_ptr, workspace->bond_mark, "bond_mark" );
sfree(control->error_ptr, 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" );
sfree(control->error_ptr, workspace->Hdia_inv, "Hdia_inv" );
sfree(control->error_ptr, workspace->b_s, "b_s" );
sfree(control->error_ptr, workspace->b_t, "b_t" );
sfree(control->error_ptr, workspace->b_prc, "b_prc" );
sfree(control->error_ptr, workspace->b_prm, "b_prm" );
sfree(control->error_ptr, workspace->s, "s" );
sfree(control->error_ptr, workspace->t, "t" );
sfree(control->error_ptr, workspace->droptol, "droptol" );
sfree(control->error_ptr, workspace->b, "b" );
sfree(control->error_ptr, 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(control->error_ptr, workspace->h[i], "h[i]" );
sfree(control->error_ptr, 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" );
sfree(control->error_ptr, workspace->h, "h" );
sfree(control->error_ptr, workspace->v, "v" );
sfree(control->error_ptr, workspace->y, "y" );
sfree(control->error_ptr, workspace->z, "z" );
sfree(control->error_ptr, workspace->g, "g" );
sfree(control->error_ptr, workspace->hs, "hs" );
sfree(control->error_ptr, 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" );
sfree(control->error_ptr, workspace->r, "r" );
sfree(control->error_ptr, workspace->d, "d" );
sfree(control->error_ptr, workspace->q, "q" );
sfree(control->error_ptr, workspace->p, "p" );
sfree(control->error_ptr, workspace->r2, "r2" );
sfree(control->error_ptr, workspace->d2, "d2" );
sfree(control->error_ptr, workspace->q2, "q2" );
sfree(control->error_ptr, workspace->p2, "p2" );
/* integrator storage */
sfree( workspace->v_const, "v_const" );
sfree(control->error_ptr, workspace->v_const, "v_const" );
/* force related storage */
sfree( workspace->f, "f" );
sfree( workspace->CdDelta, "CdDelta" );
sfree(control->error_ptr, workspace->f, "f" );
sfree(control->error_ptr, workspace->CdDelta, "CdDelta" );
/* reductions */
#ifdef LMP_USER_OMP
if (workspace->CdDeltaReduction) sfree( workspace->CdDeltaReduction, "cddelta_reduce" );
if (workspace->forceReduction) sfree( workspace->forceReduction, "f_reduce" );
if (workspace->valence_angle_atom_myoffset) sfree( workspace->valence_angle_atom_myoffset, "valence_angle_atom_myoffset");
if (workspace->my_ext_pressReduction) sfree( workspace->my_ext_pressReduction, "ext_press_reduce");
if (workspace->CdDeltaReduction) sfree(control->error_ptr, workspace->CdDeltaReduction, "cddelta_reduce" );
if (workspace->forceReduction) sfree(control->error_ptr, workspace->forceReduction, "f_reduce" );
if (workspace->valence_angle_atom_myoffset) sfree(control->error_ptr, workspace->valence_angle_atom_myoffset, "valence_angle_atom_myoffset");
if (workspace->my_ext_pressReduction) sfree(control->error_ptr, workspace->my_ext_pressReduction, "ext_press_reduce");
#endif
}
int Allocate_Workspace( reax_system * /*system*/, control_params * control,
storage *workspace, int local_cap, int total_cap,
MPI_Comm comm, char * /*msg*/ )
char * /*msg*/ )
{
int i, total_real, total_rvec, local_rvec;
@ -220,94 +222,94 @@ int Allocate_Workspace( reax_system * /*system*/, control_params * control,
/* communication storage */
for( i = 0; i < MAX_NBRS; ++i ) {
workspace->tmp_dbl[i] = (double*)
scalloc( total_cap, sizeof(double), "tmp_dbl", comm );
scalloc(control->error_ptr, total_cap, sizeof(double), "tmp_dbl");
workspace->tmp_rvec[i] = (rvec*)
scalloc( total_cap, sizeof(rvec), "tmp_rvec", comm );
scalloc(control->error_ptr, total_cap, sizeof(rvec), "tmp_rvec");
workspace->tmp_rvec2[i] = (rvec2*)
scalloc( total_cap, sizeof(rvec2), "tmp_rvec2", comm );
scalloc(control->error_ptr, total_cap, sizeof(rvec2), "tmp_rvec2");
}
/* bond order related storage */
workspace->within_bond_box = (int*)
scalloc( total_cap, sizeof(int), "skin", comm );
workspace->total_bond_order = (double*) smalloc( total_real, "total_bo", comm );
workspace->Deltap = (double*) smalloc( total_real, "Deltap", comm );
workspace->Deltap_boc = (double*) smalloc( total_real, "Deltap_boc", comm );
workspace->dDeltap_self = (rvec*) smalloc( total_rvec, "dDeltap_self", comm );
workspace->Delta = (double*) smalloc( total_real, "Delta", comm );
workspace->Delta_lp = (double*) smalloc( total_real, "Delta_lp", comm );
scalloc(control->error_ptr, total_cap, sizeof(int), "skin");
workspace->total_bond_order = (double*) smalloc(control->error_ptr, total_real, "total_bo");
workspace->Deltap = (double*) smalloc(control->error_ptr, total_real, "Deltap");
workspace->Deltap_boc = (double*) smalloc(control->error_ptr, total_real, "Deltap_boc");
workspace->dDeltap_self = (rvec*) smalloc(control->error_ptr, total_rvec, "dDeltap_self");
workspace->Delta = (double*) smalloc(control->error_ptr, total_real, "Delta");
workspace->Delta_lp = (double*) smalloc(control->error_ptr, total_real, "Delta_lp");
workspace->Delta_lp_temp = (double*)
smalloc( total_real, "Delta_lp_temp", comm );
workspace->dDelta_lp = (double*) smalloc( total_real, "dDelta_lp", comm );
smalloc(control->error_ptr, total_real, "Delta_lp_temp");
workspace->dDelta_lp = (double*) smalloc(control->error_ptr, total_real, "dDelta_lp");
workspace->dDelta_lp_temp = (double*)
smalloc( total_real, "dDelta_lp_temp", comm );
workspace->Delta_e = (double*) smalloc( total_real, "Delta_e", comm );
workspace->Delta_boc = (double*) smalloc( total_real, "Delta_boc", comm );
workspace->Delta_val = (double*) smalloc( total_real, "Delta_val", comm );
workspace->nlp = (double*) smalloc( total_real, "nlp", comm );
workspace->nlp_temp = (double*) smalloc( total_real, "nlp_temp", comm );
workspace->Clp = (double*) smalloc( total_real, "Clp", comm );
workspace->vlpex = (double*) smalloc( total_real, "vlpex", comm );
smalloc(control->error_ptr, total_real, "dDelta_lp_temp");
workspace->Delta_e = (double*) smalloc(control->error_ptr, total_real, "Delta_e");
workspace->Delta_boc = (double*) smalloc(control->error_ptr, total_real, "Delta_boc");
workspace->Delta_val = (double*) smalloc(control->error_ptr, total_real, "Delta_val");
workspace->nlp = (double*) smalloc(control->error_ptr, total_real, "nlp");
workspace->nlp_temp = (double*) smalloc(control->error_ptr, total_real, "nlp_temp");
workspace->Clp = (double*) smalloc(control->error_ptr, total_real, "Clp");
workspace->vlpex = (double*) smalloc(control->error_ptr, total_real, "vlpex");
workspace->bond_mark = (int*)
scalloc( total_cap, sizeof(int), "bond_mark", comm );
scalloc(control->error_ptr, total_cap, sizeof(int), "bond_mark");
workspace->done_after = (int*)
scalloc( total_cap, sizeof(int), "done_after", comm );
scalloc(control->error_ptr, total_cap, sizeof(int), "done_after");
/* QEq storage */
workspace->Hdia_inv = (double*)
scalloc( total_cap, sizeof(double), "Hdia_inv", comm );
workspace->b_s = (double*) scalloc( total_cap, sizeof(double), "b_s", comm );
workspace->b_t = (double*) scalloc( total_cap, sizeof(double), "b_t", comm );
workspace->b_prc = (double*) scalloc( total_cap, sizeof(double), "b_prc", comm );
workspace->b_prm = (double*) scalloc( total_cap, sizeof(double), "b_prm", comm );
workspace->s = (double*) scalloc( total_cap, sizeof(double), "s", comm );
workspace->t = (double*) scalloc( total_cap, sizeof(double), "t", comm );
scalloc(control->error_ptr, total_cap, sizeof(double), "Hdia_inv");
workspace->b_s = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "b_s");
workspace->b_t = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "b_t");
workspace->b_prc = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "b_prc");
workspace->b_prm = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "b_prm");
workspace->s = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "s");
workspace->t = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "t");
workspace->droptol = (double*)
scalloc( total_cap, sizeof(double), "droptol", comm );
workspace->b = (rvec2*) scalloc( total_cap, sizeof(rvec2), "b", comm );
workspace->x = (rvec2*) scalloc( total_cap, sizeof(rvec2), "x", comm );
scalloc(control->error_ptr, total_cap, sizeof(double), "droptol");
workspace->b = (rvec2*) scalloc(control->error_ptr, total_cap, sizeof(rvec2), "b");
workspace->x = (rvec2*) scalloc(control->error_ptr, total_cap, sizeof(rvec2), "x");
/* GMRES storage */
workspace->y = (double*) scalloc( RESTART+1, sizeof(double), "y", comm );
workspace->z = (double*) scalloc( RESTART+1, sizeof(double), "z", comm );
workspace->g = (double*) scalloc( RESTART+1, sizeof(double), "g", comm );
workspace->h = (double**) scalloc( RESTART+1, sizeof(double*), "h", comm );
workspace->hs = (double*) scalloc( RESTART+1, sizeof(double), "hs", comm );
workspace->hc = (double*) scalloc( RESTART+1, sizeof(double), "hc", comm );
workspace->v = (double**) scalloc( RESTART+1, sizeof(double*), "v", comm );
workspace->y = (double*) scalloc(control->error_ptr, RESTART+1, sizeof(double), "y");
workspace->z = (double*) scalloc(control->error_ptr, RESTART+1, sizeof(double), "z");
workspace->g = (double*) scalloc(control->error_ptr, RESTART+1, sizeof(double), "g");
workspace->h = (double**) scalloc(control->error_ptr, RESTART+1, sizeof(double*), "h");
workspace->hs = (double*) scalloc(control->error_ptr, RESTART+1, sizeof(double), "hs");
workspace->hc = (double*) scalloc(control->error_ptr, RESTART+1, sizeof(double), "hc");
workspace->v = (double**) scalloc(control->error_ptr, RESTART+1, sizeof(double*), "v");
for( i = 0; i < RESTART+1; ++i ) {
workspace->h[i] = (double*) scalloc( RESTART+1, sizeof(double), "h[i]", comm );
workspace->v[i] = (double*) scalloc( total_cap, sizeof(double), "v[i]", comm );
workspace->h[i] = (double*) scalloc(control->error_ptr, RESTART+1, sizeof(double), "h[i]");
workspace->v[i] = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "v[i]");
}
/* CG storage */
workspace->r = (double*) scalloc( total_cap, sizeof(double), "r", comm );
workspace->d = (double*) scalloc( total_cap, sizeof(double), "d", comm );
workspace->q = (double*) scalloc( total_cap, sizeof(double), "q", comm );
workspace->p = (double*) scalloc( total_cap, sizeof(double), "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 );
workspace->r = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "r");
workspace->d = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "d");
workspace->q = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "q");
workspace->p = (double*) scalloc(control->error_ptr, total_cap, sizeof(double), "p");
workspace->r2 = (rvec2*) scalloc(control->error_ptr, total_cap, sizeof(rvec2), "r2");
workspace->d2 = (rvec2*) scalloc(control->error_ptr, total_cap, sizeof(rvec2), "d2");
workspace->q2 = (rvec2*) scalloc(control->error_ptr, total_cap, sizeof(rvec2), "q2");
workspace->p2 = (rvec2*) scalloc(control->error_ptr, total_cap, sizeof(rvec2), "p2");
/* integrator storage */
workspace->v_const = (rvec*) smalloc( local_rvec, "v_const", comm );
workspace->v_const = (rvec*) smalloc(control->error_ptr, local_rvec, "v_const");
/* force related storage */
workspace->f = (rvec*) scalloc( total_cap, sizeof(rvec), "f", comm );
workspace->f = (rvec*) scalloc(control->error_ptr, total_cap, sizeof(rvec), "f");
workspace->CdDelta = (double*)
scalloc( total_cap, sizeof(double), "CdDelta", comm );
scalloc(control->error_ptr, total_cap, sizeof(double), "CdDelta");
// storage for reductions with multiple threads
#ifdef LMP_USER_OMP
workspace->CdDeltaReduction = (double *) scalloc(sizeof(double), total_cap*control->nthreads,
"cddelta_reduce", comm);
workspace->CdDeltaReduction = (double *) scalloc(control->error_ptr, sizeof(double), total_cap*control->nthreads,
"cddelta_reduce");
workspace->forceReduction = (rvec *) scalloc(sizeof(rvec), total_cap*control->nthreads,
"forceReduction", comm);
workspace->forceReduction = (rvec *) scalloc(control->error_ptr, sizeof(rvec), total_cap*control->nthreads,
"forceReduction");
workspace->valence_angle_atom_myoffset = (int *) scalloc(sizeof(int), total_cap, "valence_angle_atom_myoffset", comm);
workspace->valence_angle_atom_myoffset = (int *) scalloc(control->error_ptr, sizeof(int), total_cap, "valence_angle_atom_myoffset");
workspace->my_ext_pressReduction = (rvec *) calloc(sizeof(rvec), control->nthreads);
#else
LMP_UNUSED_PARAM(control);
@ -318,18 +320,16 @@ int Allocate_Workspace( reax_system * /*system*/, control_params * control,
static void Reallocate_Neighbor_List( reax_list *far_nbrs, int n,
int num_intrs, MPI_Comm comm )
int num_intrs )
{
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 );
Delete_List( far_nbrs);
if(!Make_List( n, num_intrs, TYP_FAR_NEIGHBOR, far_nbrs )){
far_nbrs->error_ptr->one(FLERR,"Problem in initializing far neighbors list");
}
}
static int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
MPI_Comm comm )
static int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds )
{
int i, total_hbonds;
@ -343,10 +343,9 @@ static int Reallocate_HBonds_List( reax_system *system, reax_list *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 );
Delete_List( hbonds);
if (!Make_List( system->Hcap, total_hbonds, TYP_HBOND, hbonds )) {
hbonds->error_ptr->one(FLERR, "Not enough space for hydrogen bonds list");
}
return total_hbonds;
@ -354,8 +353,7 @@ static int Reallocate_HBonds_List( reax_system *system, reax_list *hbonds,
static int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
int *total_bonds, int *est_3body,
MPI_Comm comm )
int *total_bonds, int *est_3body )
{
int i;
@ -373,13 +371,12 @@ static int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
#ifdef LMP_USER_OMP
if (system->omp_active)
for (i = 0; i < bonds->num_intrs; ++i)
sfree(bonds->select.bond_list[i].bo_data.CdboReduction, "CdboReduction");
sfree(system->error_ptr, bonds->select.bond_list[i].bo_data.CdboReduction, "CdboReduction");
#endif
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 );
Delete_List( bonds);
if(!Make_List(system->total_cap, *total_bonds, TYP_BOND, bonds)) {
bonds->error_ptr->one(FLERR, "Not enough space for bonds list");
}
#ifdef LMP_USER_OMP
@ -392,7 +389,7 @@ static int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
if (system->omp_active)
for (i = 0; i < bonds->num_intrs; ++i)
bonds->select.bond_list[i].bo_data.CdboReduction =
(double*) smalloc(sizeof(double)*nthreads, "CdboReduction", comm);
(double*) smalloc(system->error_ptr, sizeof(double)*nthreads, "CdboReduction");
#endif
return SUCCESS;
@ -400,14 +397,12 @@ static int Reallocate_Bonds_List( reax_system *system, reax_list *bonds,
void ReAllocate( reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists,
mpi_datatypes *mpi_data )
simulation_data *data, storage *workspace, reax_list **lists )
{
int num_bonds, est_3body, Hflag, ret;
int renbr, newsize;
reallocate_data *realloc;
reax_list *far_nbrs;
MPI_Comm comm;
char msg[200];
int mincap = system->mincap;
@ -415,7 +410,6 @@ void ReAllocate( reax_system *system, control_params *control,
double saferzone = system->saferzone;
realloc = &(workspace->realloc);
comm = mpi_data->world;
if( system->n >= DANGER_ZONE * system->local_cap ||
(0 && system->n <= LOOSE_ZONE * system->local_cap) ) {
@ -433,21 +427,19 @@ void ReAllocate( reax_system *system, control_params *control,
/* system */
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 );
char errmsg[256];
snprintf(errmsg, 256, "Not enough space for atom_list: total_cap=%d", system->total_cap);
system->error_ptr->one(FLERR, errmsg);
}
/* workspace */
DeAllocate_Workspace( control, workspace );
ret = Allocate_Workspace( system, control, workspace, system->local_cap,
system->total_cap, comm, msg );
system->total_cap, 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 );
char errmsg[256];
snprintf(errmsg, 256, "Not enough space for workspace: local_cap=%d total_cap=%d", system->local_cap, system->total_cap);
system->error_ptr->one(FLERR, errmsg);
}
}
@ -459,15 +451,15 @@ void ReAllocate( reax_system *system, control_params *control,
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 );
char errmsg[256];
snprintf(errmsg, 256, "step%d-ran out of space on far_nbrs: top=%d, max=%d", data->step, realloc->num_far, far_nbrs->num_intrs);
system->error_ptr->one(FLERR, errmsg);
}
newsize = static_cast<int>
(MAX( realloc->num_far*safezone, mincap*MIN_NBRS ));
Reallocate_Neighbor_List( far_nbrs, system->total_cap, newsize, comm );
Reallocate_Neighbor_List( far_nbrs, system->total_cap, newsize);
realloc->num_far = 0;
}
}
@ -482,7 +474,7 @@ void ReAllocate( reax_system *system, control_params *control,
}
if (Hflag || realloc->hbonds) {
ret = Reallocate_HBonds_List( system, (*lists)+HBONDS, comm );
ret = Reallocate_HBonds_List( system, (*lists)+HBONDS);
realloc->hbonds = 0;
}
}
@ -491,14 +483,14 @@ void ReAllocate( reax_system *system, control_params *control,
num_bonds = est_3body = -1;
if (Nflag || realloc->bonds) {
Reallocate_Bonds_List( system, (*lists)+BONDS, &num_bonds,
&est_3body, comm );
&est_3body);
realloc->bonds = 0;
realloc->num_3body = MAX( realloc->num_3body, est_3body ) * 2;
}
/* 3-body list */
if (realloc->num_3body > 0) {
Delete_List( (*lists)+THREE_BODIES, comm );
Delete_List( (*lists)+THREE_BODIES);
if (num_bonds == -1)
num_bonds = ((*lists)+BONDS)->num_intrs;
@ -506,9 +498,8 @@ void ReAllocate( reax_system *system, control_params *control,
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 );
(*lists)+THREE_BODIES ) ) {
system->error_ptr->one(FLERR, "Problem in initializing angles list");
}
realloc->num_3body = -1;
}

6
src/USER-REAXC/reaxc_allocate.h
View File

@ -28,15 +28,15 @@
#define __ALLOCATE_H_
#include "reaxc_types.h"
int PreAllocate_Space( reax_system*, control_params*, storage*, MPI_Comm );
int PreAllocate_Space( reax_system*, control_params*, storage* );
int Allocate_System( reax_system*, int, int, char* );
void DeAllocate_System( reax_system* );
int Allocate_Workspace( reax_system*, control_params*, storage*,
int, int, MPI_Comm, char* );
int, int, char* );
void DeAllocate_Workspace( control_params*, storage* );
void ReAllocate( reax_system*, control_params*, simulation_data*, storage*,
reax_list**, mpi_datatypes* );
reax_list** );
#endif

10
src/USER-REAXC/reaxc_control.cpp
View File

@ -28,6 +28,8 @@
#include "reaxc_control.h"
#include "reaxc_tool_box.h"
#include "error.h"
char Read_Control_File( char *control_file, control_params* control,
output_controls *out_control )
{
@ -38,8 +40,7 @@ char Read_Control_File( char *control_file, control_params* control,
/* open control file */
if ( (fp = fopen( control_file, "r" ) ) == NULL ) {
fprintf( stderr, "error opening the control file! terminating...\n" );
MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
control->error_ptr->all(FLERR, "The control file cannot be opened");
}
/* assign default values */
@ -364,8 +365,9 @@ char Read_Control_File( char *control_file, control_params* control,
control->restrict_type = ival;
}
else {
fprintf( stderr, "WARNING: unknown parameter %s\n", tmp[0] );
MPI_Abort( MPI_COMM_WORLD, 15 );
char errmsg[128];
snprintf(errmsg,128,"Unknown parameter %s in the control file", tmp[0]);
control->error_ptr->all(FLERR, errmsg);
}
}

159
src/USER-REAXC/reaxc_ffield.cpp
View File

@ -29,6 +29,7 @@
#include "reaxc_ffield.h"
#include "reaxc_tool_box.h"
char Read_Force_Field( FILE *fp, reax_interaction *reax,
control_params *control )
{
@ -41,7 +42,6 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
double val;
MPI_Comm comm;
int me;
comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm, &me);
@ -61,7 +61,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
n = atoi(tmp[0]);
if (n < 1) {
if (me == 0)
fprintf( stderr, "WARNING: number of globals in ffield file is 0!\n" );
control->error_ptr->warning( FLERR, "Number of globals in ffield file is 0. The file will not be read." );
fclose(fp);
free(s);
free(tmp);
@ -96,61 +96,52 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
/* Allocating structures in reax_interaction */
reax->sbp = (single_body_parameters*)
scalloc( reax->num_atom_types, sizeof(single_body_parameters), "sbp",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(single_body_parameters), "sbp");
reax->tbp = (two_body_parameters**)
scalloc( reax->num_atom_types, sizeof(two_body_parameters*), "tbp", comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(two_body_parameters*), "tbp");
reax->thbp= (three_body_header***)
scalloc( reax->num_atom_types, sizeof(three_body_header**), "thbp", comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(three_body_header**), "thbp");
reax->hbp = (hbond_parameters***)
scalloc( reax->num_atom_types, sizeof(hbond_parameters**), "hbp", comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(hbond_parameters**), "hbp");
reax->fbp = (four_body_header****)
scalloc( reax->num_atom_types, sizeof(four_body_header***), "fbp", comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(four_body_header***), "fbp");
tor_flag = (char****)
scalloc( reax->num_atom_types, sizeof(char***), "tor_flag", comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(char***), "tor_flag");
for( i = 0; i < reax->num_atom_types; i++ ) {
reax->tbp[i] = (two_body_parameters*)
scalloc( reax->num_atom_types, sizeof(two_body_parameters), "tbp[i]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(two_body_parameters), "tbp[i]");
reax->thbp[i]= (three_body_header**)
scalloc( reax->num_atom_types, sizeof(three_body_header*), "thbp[i]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(three_body_header*), "thbp[i]");
reax->hbp[i] = (hbond_parameters**)
scalloc( reax->num_atom_types, sizeof(hbond_parameters*), "hbp[i]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(hbond_parameters*), "hbp[i]");
reax->fbp[i] = (four_body_header***)
scalloc( reax->num_atom_types, sizeof(four_body_header**), "fbp[i]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(four_body_header**), "fbp[i]");
tor_flag[i] = (char***)
scalloc( reax->num_atom_types, sizeof(char**), "tor_flag[i]", comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(char**), "tor_flag[i]");
for( j = 0; j < reax->num_atom_types; j++ ) {
reax->thbp[i][j]= (three_body_header*)
scalloc( reax->num_atom_types, sizeof(three_body_header), "thbp[i,j]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(three_body_header), "thbp[i,j]");
reax->hbp[i][j] = (hbond_parameters*)
scalloc( reax->num_atom_types, sizeof(hbond_parameters), "hbp[i,j]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(hbond_parameters), "hbp[i,j]");
reax->fbp[i][j] = (four_body_header**)
scalloc( reax->num_atom_types, sizeof(four_body_header*), "fbp[i,j]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(four_body_header*), "fbp[i,j]");
tor_flag[i][j] = (char**)
scalloc( reax->num_atom_types, sizeof(char*), "tor_flag[i,j]", comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(char*), "tor_flag[i,j]");
for (k=0; k < reax->num_atom_types; k++) {
reax->fbp[i][j][k] = (four_body_header*)
scalloc( reax->num_atom_types, sizeof(four_body_header), "fbp[i,j,k]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(four_body_header), "fbp[i,j,k]");
tor_flag[i][j][k] = (char*)
scalloc( reax->num_atom_types, sizeof(char), "tor_flag[i,j,k]",
comm );
scalloc(control->error_ptr, reax->num_atom_types, sizeof(char), "tor_flag[i,j,k]");
}
}
}
reax->gp.vdw_type = 0;
char errmsg[1024];
for( i = 0; i < reax->num_atom_types; i++ ) {
/* line one */
@ -158,15 +149,12 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
c = Tokenize( s, &tmp );
/* Sanity checks */
if (c == 2 && !lgflag) {
if (me == 0)
fprintf(stderr, "Force field file requires using 'lgvdw yes'\n");
MPI_Abort( comm, FILE_NOT_FOUND );
}
if (c == 2 && !lgflag)
control->error_ptr->all(FLERR, "Force field file requires using 'lgvdw yes'");
if (c < 9) {
if (me == 0)
fprintf(stderr, "Inconsistent ffield file (reaxc_ffield.cpp) \n");
MPI_Abort( comm, FILE_NOT_FOUND );
snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);
control->error_ptr->all(FLERR, errmsg);
}
for( j = 0; j < (int)(strlen(tmp[0])); ++j )
@ -188,9 +176,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
/* Sanity check */
if (c < 8) {
if (me == 0)
fprintf(stderr, "Inconsistent ffield file (reaxc_ffield.cpp) \n");
MPI_Abort( comm, FILE_NOT_FOUND );
snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);
control->error_ptr->all(FLERR, errmsg);
}
val = atof(tmp[0]); reax->sbp[i].alpha = val;
@ -208,9 +195,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
/* Sanity check */
if (c < 8) {
if (me == 0)
fprintf(stderr, "Inconsistent ffield file (reaxc_ffield.cpp) \n");
MPI_Abort( comm, FILE_NOT_FOUND );
snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);
control->error_ptr->all(FLERR, errmsg);
}
val = atof(tmp[0]); reax->sbp[i].r_pi_pi = val;
@ -228,9 +214,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
/* Sanity check */
if (c < 8) {
if (me == 0)
fprintf(stderr, "Inconsistent ffield file (reaxc_ffield.cpp) \n");
MPI_Abort( comm, FILE_NOT_FOUND );
snprintf (errmsg, 1024, "Missing parameter(s) in line %s", s);
control->error_ptr->all(FLERR, errmsg);
}
val = atof(tmp[0]); reax->sbp[i].p_ovun2 = val;
@ -249,9 +234,7 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
/* Sanity check */
if (c > 2) {
if (me == 0)
fprintf(stderr, "Force field file incompatible with 'lgvdw yes'\n");
MPI_Abort( comm, FILE_NOT_FOUND );
control->error_ptr->all(FLERR,"Force field file incompatible with 'lgvdw yes'");
}
val = atof(tmp[0]); reax->sbp[i].lgcij = val;
@ -261,28 +244,32 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
if (reax->sbp[i].rcore2>0.01 && reax->sbp[i].acore2>0.01) { // Inner-wall
if (reax->sbp[i].gamma_w>0.5) { // Shielding vdWaals
if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3) {
if (errorflag && (me == 0))
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
"Force field parameters for element %s\n" \
"indicate inner wall+shielding, but earlier\n" \
"atoms indicate different vdWaals-method.\n" \
"This may cause division-by-zero errors.\n" \
"Keeping vdWaals-setting for earlier atoms.\n",
reax->sbp[i].name );
if (errorflag && (me == 0)) {
char errmsg[512];
snprintf(errmsg, 512, "VdWaals-parameters for element %s "
"indicate inner wall+shielding, but earlier "
"atoms indicate different vdWaals-method. "
"This may cause division-by-zero errors. "
"Keeping vdWaals-setting for earlier atoms.",
reax->sbp[i].name);
control->error_ptr->warning(FLERR,errmsg);
}
errorflag = 0;
} else {
reax->gp.vdw_type = 3;
}
} else { // No shielding vdWaals parameters present
if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 2) {
if (me == 0)
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
"Force field parameters for element %s\n" \
"indicate inner wall without shielding, but earlier\n" \
"atoms indicate different vdWaals-method.\n" \
"This may cause division-by-zero errors.\n" \
"Keeping vdWaals-setting for earlier atoms.\n",
reax->sbp[i].name );
if (me == 0) {
char errmsg[512];
snprintf(errmsg, 512, "VdWaals-parameters for element %s "
"indicate inner wall without shielding, but earlier "
"atoms indicate different vdWaals-method. "
"This may cause division-by-zero errors. "
"Keeping vdWaals-setting for earlier atoms.",
reax->sbp[i].name);
control->error_ptr->warning(FLERR,errmsg);
}
} else {
reax->gp.vdw_type = 2;
}
@ -290,23 +277,25 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
} else { // No Inner wall parameters present
if (reax->sbp[i].gamma_w>0.5) { // Shielding vdWaals
if (reax->gp.vdw_type != 0 && reax->gp.vdw_type != 1) {
if (me == 0)
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
"Force field parameters for element %s\n" \
"indicate shielding without inner wall, but earlier\n" \
"atoms indicate different vdWaals-method.\n" \
"This may cause division-by-zero errors.\n" \
"Keeping vdWaals-setting for earlier atoms.\n",
reax->sbp[i].name );
if (me == 0) {
char errmsg[512];
snprintf(errmsg, 512, "VdWaals parameters for element %s "
"indicate shielding without inner wall, but earlier "
"elements indicate different vdWaals-method. "
"This may cause division-by-zero errors. "
"Keeping vdWaals-setting for earlier atoms.",
reax->sbp[i].name);
control->error_ptr->warning(FLERR,errmsg);
}
} else {
reax->gp.vdw_type = 1;
}
} else {
if (me == 0)
fprintf( stderr, "Error: inconsistent vdWaals-parameters\n" \
"No shielding or inner-wall set for element %s\n",
reax->sbp[i].name );
MPI_Abort( comm, INVALID_INPUT );
char errmsg[256];
snprintf(errmsg, 256, "Inconsistent vdWaals-parameters: "
"No shielding or inner-wall set for element %s",
reax->sbp[i].name);
control->error_ptr->all(FLERR, errmsg);
}
}
}
@ -315,15 +304,23 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
for( i = 0; i < reax->num_atom_types; i++ )
if( reax->sbp[i].mass < 21 &&
reax->sbp[i].valency_val != reax->sbp[i].valency_boc ) {
if (me == 0)
fprintf(stderr,"Warning: changed valency_val to valency_boc for %s\n",
reax->sbp[i].name );
if (me == 0) {
char errmsg[256];
snprintf(errmsg, 256, "Changed valency_val to valency_boc for %s",
reax->sbp[i].name);
control->error_ptr->warning(FLERR,errmsg);
}
reax->sbp[i].valency_val = reax->sbp[i].valency_boc;
}
/* next line is number of two body combination and some comments */
fgets(s,MAX_LINE,fp);
c=Tokenize(s,&tmp);
if (c == 2 && !lgflag) {
control->error_ptr->all(FLERR, "Force field file requires using 'lgvdw yes'");
}
l = atoi(tmp[0]);
/* a line of comments */

27
src/USER-REAXC/reaxc_forces.cpp
View File

@ -39,11 +39,14 @@
#include "reaxc_valence_angles.h"
#include "reaxc_vector.h"
#include "error.h"
interaction_function Interaction_Functions[NUM_INTRS];
void Dummy_Interaction( reax_system * /*system*/, control_params * /*control*/,
simulation_data * /*data*/, storage * /*workspace*/,
reax_list **/*lists*/, output_controls * /*out_control*/ )
reax_list ** /*lists*/, output_controls * /*out_control*/ )
{
}
@ -115,7 +118,7 @@ void Compute_Total_Force( reax_system *system, control_params *control,
}
void Validate_Lists( reax_system *system, storage * /*workspace*/, reax_list **lists,
int step, int /*n*/, int N, int numH, MPI_Comm comm )
int step, int /*n*/, int N, int numH )
{
int i, comp, Hindex;
reax_list *bonds, *hbonds;
@ -134,9 +137,10 @@ void Validate_Lists( reax_system *system, storage * /*workspace*/, reax_list **l
else comp = bonds->num_intrs;
if (End_Index(i, bonds) > comp) {
fprintf( stderr, "step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
char errmsg[256];
snprintf(errmsg, 256, "step%d-bondchk failed: i=%d end(i)=%d str(i+1)=%d\n",
step, i, End_Index(i,bonds), comp );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
system->error_ptr->one(FLERR,errmsg);
}
}
}
@ -161,9 +165,10 @@ void Validate_Lists( reax_system *system, storage * /*workspace*/, reax_list **l
else comp = hbonds->num_intrs;
if (End_Index(Hindex, hbonds) > comp) {
fprintf(stderr,"step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
char errmsg[256];
snprintf(errmsg, 256, "step%d-hbondchk failed: H=%d end(H)=%d str(H+1)=%d\n",
step, Hindex, End_Index(Hindex,hbonds), comp );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
system->error_ptr->one(FLERR, errmsg);
}
}
}
@ -173,8 +178,7 @@ void Validate_Lists( reax_system *system, storage * /*workspace*/, reax_list **l
void Init_Forces_noQEq( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls * /*out_control*/,
MPI_Comm comm ) {
reax_list **lists, output_controls * /*out_control*/ ) {
int i, j, pj;
int start_i, end_i;
int type_i, type_j;
@ -308,13 +312,13 @@ void Init_Forces_noQEq( reax_system *system, control_params *control,
workspace->realloc.num_hbonds = num_hbonds;
Validate_Lists( system, workspace, lists, data->step,
system->n, system->N, system->numH, comm );
system->n, system->N, system->numH);
}
void Estimate_Storages( reax_system *system, control_params *control,
reax_list **lists, int *Htop, int *hb_top,
int *bond_top, int *num_3body, MPI_Comm /*comm*/ )
int *bond_top, int *num_3body )
{
int i, j, pj;
int start_i, end_i;
@ -436,10 +440,9 @@ void Compute_Forces( reax_system *system, control_params *control,
reax_list **lists, output_controls *out_control,
mpi_datatypes *mpi_data )
{
MPI_Comm comm = mpi_data->world;
Init_Forces_noQEq( system, control, data, workspace,
lists, out_control, comm );
lists, out_control);
/********* bonded interactions ************/
Compute_Bonded_Forces( system, control, data, workspace,

2
src/USER-REAXC/reaxc_forces.h
View File

@ -36,5 +36,5 @@ void Init_Force_Functions( control_params* );
void Compute_Forces( reax_system*, control_params*, simulation_data*,
storage*, reax_list**, output_controls*, mpi_datatypes* );
void Estimate_Storages( reax_system*, control_params*, reax_list**,
int*, int*, int*, int*, MPI_Comm );
int*, int*, int*, int* );
#endif

86
src/USER-REAXC/reaxc_init_md.cpp
View File

@ -36,6 +36,8 @@
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
#include "error.h"
int Init_System( reax_system *system, control_params *control, char * /*msg*/ )
{
int i;
@ -80,7 +82,7 @@ int Init_Simulation_Data( reax_system *system, control_params *control,
return SUCCESS;
}
void Init_Taper( control_params *control, storage *workspace, MPI_Comm comm )
void Init_Taper( control_params *control, storage *workspace )
{
double d1, d7;
double swa, swa2, swa3;
@ -89,15 +91,17 @@ void Init_Taper( control_params *control, storage *workspace, MPI_Comm comm )
swa = control->nonb_low;
swb = control->nonb_cut;
if (fabs( swa ) > 0.01)
fprintf( stderr, "Warning: non-zero lower Taper-radius cutoff\n" );
if (fabs( swa ) > 0.01 && control->me == 0)
control->error_ptr->warning( FLERR, "Non-zero lower Taper-radius cutoff" );
if (swb < 0) {
fprintf( stderr, "Negative upper Taper-radius cutoff\n" );
MPI_Abort( comm, INVALID_INPUT );
control->error_ptr->all(FLERR,"Negative upper Taper-radius cutoff");
}
else if( swb < 5 && control->me == 0) {
char errmsg[256];
snprintf(errmsg, 256, "Very low Taper-radius cutoff: %f", swb );
control->error_ptr->warning( FLERR, errmsg );
}
else if( swb < 5 )
fprintf( stderr, "Warning: very low Taper-radius cutoff: %f\n", swb );
d1 = swb - swa;
d7 = pow( d1, 7.0 );
@ -119,12 +123,12 @@ void Init_Taper( control_params *control, storage *workspace, MPI_Comm comm )
int Init_Workspace( reax_system *system, control_params *control,
storage *workspace, MPI_Comm comm, char *msg )
storage *workspace, char *msg )
{
int ret;
ret = Allocate_Workspace( system, control, workspace,
system->local_cap, system->total_cap, comm, msg );
system->local_cap, system->total_cap, msg );
if (ret != SUCCESS)
return ret;
@ -132,7 +136,7 @@ int Init_Workspace( reax_system *system, control_params *control,
Reset_Workspace( system, workspace );
/* Initialize the Taper function */
Init_Taper( control, workspace, comm );
Init_Taper( control, workspace);
return SUCCESS;
}
@ -156,17 +160,15 @@ int Init_Lists( reax_system *system, control_params *control,
{
int i, total_hbonds, total_bonds, bond_cap, num_3body, cap_3body, Htop;
int *hb_top, *bond_top;
MPI_Comm comm;
int mincap = system->mincap;
double safezone = system->safezone;
double saferzone = system->saferzone;
comm = mpi_data->world;
bond_top = (int*) calloc( system->total_cap, sizeof(int) );
hb_top = (int*) calloc( system->local_cap, sizeof(int) );
Estimate_Storages( system, control, lists,
&Htop, hb_top, bond_top, &num_3body, comm );
&Htop, hb_top, bond_top, &num_3body);
if (control->hbond_cut > 0) {
/* init H indexes */
@ -178,10 +180,10 @@ int Init_Lists( reax_system *system, control_params *control,
total_hbonds = (int)(MAX( total_hbonds*saferzone, mincap*MIN_HBONDS ));
if( !Make_List( system->Hcap, total_hbonds, TYP_HBOND,
*lists+HBONDS, comm ) ) {
fprintf( stderr, "not enough space for hbonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
*lists+HBONDS ) ) {
control->error_ptr->one(FLERR, "Not enough space for hbonds list.");
}
(*lists+HBONDS)->error_ptr = system->error_ptr;
}
total_bonds = 0;
@ -192,18 +194,18 @@ int Init_Lists( reax_system *system, control_params *control,
bond_cap = (int)(MAX( total_bonds*safezone, mincap*MIN_BONDS ));
if( !Make_List( system->total_cap, bond_cap, TYP_BOND,
*lists+BONDS, comm ) ) {
fprintf( stderr, "not enough space for bonds list. terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
*lists+BONDS ) ) {
control->error_ptr->one(FLERR, "Not enough space for bonds list.");
}
(*lists+BONDS)->error_ptr = system->error_ptr;
/* 3bodies list */
cap_3body = (int)(MAX( num_3body*safezone, MIN_3BODIES ));
if( !Make_List( bond_cap, cap_3body, TYP_THREE_BODY,
*lists+THREE_BODIES, comm ) ){
fprintf( stderr, "Problem in initializing angles list. Terminating!\n" );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
*lists+THREE_BODIES ) ){
control->error_ptr->one(FLERR,"Problem in initializing angles list.");
}
(*lists+THREE_BODIES)->error_ptr = system->error_ptr;
free( hb_top );
free( bond_top );
@ -217,60 +219,38 @@ void Initialize( reax_system *system, control_params *control,
mpi_datatypes *mpi_data, MPI_Comm comm )
{
char msg[MAX_STR];
char errmsg[128];
if (Init_MPI_Datatypes(system, workspace, mpi_data, comm, msg) == FAILURE) {
fprintf( stderr, "p%d: init_mpi_datatypes: could not create datatypes\n",
system->my_rank );
fprintf( stderr, "p%d: mpi_data couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
control->error_ptr->one(FLERR,"Could not create datatypes");
}
if (Init_System(system, control, msg) == FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: system could not be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
control->error_ptr->one(FLERR,"System could not be initialized");
}
if (Init_Simulation_Data( system, control, data, msg ) == FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: sim_data couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
control->error_ptr->one(FLERR,"Sim_data could not be initialized");
}
if (Init_Workspace( system, control, workspace, mpi_data->world, msg ) ==
if (Init_Workspace( system, control, workspace, msg ) ==
FAILURE) {
fprintf( stderr, "p%d:init_workspace: not enough memory\n",
system->my_rank );
fprintf( stderr, "p%d:workspace couldn't be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
control->error_ptr->one(FLERR,"Workspace could not be initialized");
}
if (Init_Lists( system, control, data, workspace, lists, mpi_data, msg ) ==
FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: system could not be initialized! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
control->error_ptr->one(FLERR,"Lists could not be initialized");
}
if (Init_Output_Files(system,control,out_control,mpi_data,msg)== FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: could not open output files! terminating...\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
control->error_ptr->one(FLERR,"Could not open output files");
}
if (control->tabulate) {
if (Init_Lookup_Tables( system, control, workspace, mpi_data, msg ) == FAILURE) {
fprintf( stderr, "p%d: %s\n", system->my_rank, msg );
fprintf( stderr, "p%d: couldn't create lookup table! terminating.\n",
system->my_rank );
MPI_Abort( mpi_data->world, CANNOT_INITIALIZE );
control->error_ptr->one(FLERR,"Lookup table could not be created");
}
}

1
src/USER-REAXC/reaxc_init_md.h
View File

@ -29,6 +29,7 @@
#include "reaxc_types.h"
void Initialize( reax_system*, control_params*, simulation_data*, storage*,
reax_list**, output_controls*, mpi_datatypes*, MPI_Comm );
#endif

72
src/USER-REAXC/reaxc_list.cpp
View File

@ -28,117 +28,121 @@
#include "reaxc_list.h"
#include "reaxc_tool_box.h"
#include "error.h"
/************* allocate list space ******************/
int Make_List(int n, int num_intrs, int type, reax_list *l, MPI_Comm comm)
int Make_List(int n, int num_intrs, int type, reax_list *l )
{
l->allocated = 1;
l->n = n;
l->num_intrs = num_intrs;
if (l->index) sfree(l->index, "list:index");
if (l->end_index) sfree(l->end_index, "list:end_index");
l->index = (int*) smalloc( n * sizeof(int), "list:index", comm );
l->end_index = (int*) smalloc( n * sizeof(int), "list:end_index", comm );
if (l->index) sfree(l->error_ptr, l->index, "list:index");
if (l->end_index) sfree(l->error_ptr, l->end_index, "list:end_index");
l->index = (int*) smalloc(l->error_ptr, n * sizeof(int), "list:index");
l->end_index = (int*) smalloc(l->error_ptr, n * sizeof(int), "list:end_index");
l->type = type;
switch(l->type) {
case TYP_VOID:
if (l->select.v) sfree(l->select.v, "list:v");
l->select.v = (void*) smalloc(l->num_intrs * sizeof(void*), "list:v", comm);
if (l->select.v) sfree(l->error_ptr, l->select.v, "list:v");
l->select.v = (void*) smalloc(l->error_ptr, l->num_intrs * sizeof(void*), "list:v");
break;
case TYP_THREE_BODY:
if (l->select.three_body_list) sfree(l->select.three_body_list,"list:three_bodies");
if (l->select.three_body_list) sfree(l->error_ptr, l->select.three_body_list,"list:three_bodies");
l->select.three_body_list = (three_body_interaction_data*)
smalloc( l->num_intrs * sizeof(three_body_interaction_data),
"list:three_bodies", comm );
smalloc(l->error_ptr, l->num_intrs * sizeof(three_body_interaction_data),
"list:three_bodies");
break;
case TYP_BOND:
if (l->select.bond_list) sfree(l->select.bond_list,"list:bonds");
if (l->select.bond_list) sfree(l->error_ptr, l->select.bond_list,"list:bonds");
l->select.bond_list = (bond_data*)
smalloc( l->num_intrs * sizeof(bond_data), "list:bonds", comm );
smalloc(l->error_ptr, l->num_intrs * sizeof(bond_data), "list:bonds");
break;
case TYP_DBO:
if (l->select.dbo_list) sfree(l->select.dbo_list,"list:dbonds");
if (l->select.dbo_list) sfree(l->error_ptr, l->select.dbo_list,"list:dbonds");
l->select.dbo_list = (dbond_data*)
smalloc( l->num_intrs * sizeof(dbond_data), "list:dbonds", comm );
smalloc(l->error_ptr, l->num_intrs * sizeof(dbond_data), "list:dbonds");
break;
case TYP_DDELTA:
if (l->select.dDelta_list) sfree(l->select.dDelta_list,"list:dDeltas");
if (l->select.dDelta_list) sfree(l->error_ptr, l->select.dDelta_list,"list:dDeltas");
l->select.dDelta_list = (dDelta_data*)
smalloc( l->num_intrs * sizeof(dDelta_data), "list:dDeltas", comm );
smalloc(l->error_ptr, l->num_intrs * sizeof(dDelta_data), "list:dDeltas");
break;
case TYP_FAR_NEIGHBOR:
if (l->select.far_nbr_list) sfree(l->select.far_nbr_list,"list:far_nbrs");
if (l->select.far_nbr_list) sfree(l->error_ptr, l->select.far_nbr_list,"list:far_nbrs");
l->select.far_nbr_list = (far_neighbor_data*)
smalloc(l->num_intrs * sizeof(far_neighbor_data), "list:far_nbrs", comm);
smalloc(l->error_ptr, l->num_intrs * sizeof(far_neighbor_data), "list:far_nbrs");
break;
case TYP_HBOND:
if (l->select.hbond_list) sfree(l->select.hbond_list,"list:hbonds");
if (l->select.hbond_list) sfree(l->error_ptr, l->select.hbond_list,"list:hbonds");
l->select.hbond_list = (hbond_data*)
smalloc( l->num_intrs * sizeof(hbond_data), "list:hbonds", comm );
smalloc(l->error_ptr, l->num_intrs * sizeof(hbond_data), "list:hbonds");
break;
default:
fprintf( stderr, "ERROR: no %d list type defined!\n", l->type );
MPI_Abort( comm, INVALID_INPUT );
char errmsg[128];
snprintf(errmsg, 128, "No %d list type defined", l->type);
l->error_ptr->one(FLERR,errmsg);
}
return SUCCESS;
}
void Delete_List( reax_list *l, MPI_Comm comm )
void Delete_List( reax_list *l )
{
if (l->allocated == 0)
return;
l->allocated = 0;
sfree( l->index, "list:index" );
sfree( l->end_index, "list:end_index" );
sfree(l->error_ptr, l->index, "list:index" );
sfree(l->error_ptr, l->end_index, "list:end_index" );
l->index = NULL;
l->end_index = NULL;
switch(l->type) {
case TYP_VOID:
sfree( l->select.v, "list:v" );
sfree(l->error_ptr, l->select.v, "list:v" );
l->select.v = NULL;
break;
case TYP_HBOND:
sfree( l->select.hbond_list, "list:hbonds" );
sfree(l->error_ptr, l->select.hbond_list, "list:hbonds" );
l->select.hbond_list = NULL;
break;
case TYP_FAR_NEIGHBOR:
sfree( l->select.far_nbr_list, "list:far_nbrs" );
sfree(l->error_ptr, l->select.far_nbr_list, "list:far_nbrs" );
l->select.far_nbr_list = NULL;
break;
case TYP_BOND:
sfree( l->select.bond_list, "list:bonds" );
sfree(l->error_ptr, l->select.bond_list, "list:bonds" );
l->select.bond_list = NULL;
break;
case TYP_DBO:
sfree( l->select.dbo_list, "list:dbos" );
sfree(l->error_ptr, l->select.dbo_list, "list:dbos" );
l->select.dbo_list = NULL;
break;
case TYP_DDELTA:
sfree( l->select.dDelta_list, "list:dDeltas" );
sfree(l->error_ptr, l->select.dDelta_list, "list:dDeltas" );
l->select.dDelta_list = NULL;
break;
case TYP_THREE_BODY:
sfree( l->select.three_body_list, "list:three_bodies" );
sfree(l->error_ptr, l->select.three_body_list, "list:three_bodies" );
l->select.three_body_list = NULL;
break;
default:
fprintf( stderr, "ERROR: no %d list type defined!\n", l->type );
MPI_Abort( comm, INVALID_INPUT );
char errmsg[128];
snprintf(errmsg, 128, "No %d list type defined", l->type);
l->error_ptr->all(FLERR,errmsg);
}
}

4
src/USER-REAXC/reaxc_list.h
View File

@ -29,8 +29,8 @@
#include "reaxc_types.h"
int Make_List( int, int, int, reax_list*, MPI_Comm );
void Delete_List( reax_list*, MPI_Comm );
int Make_List( int, int, int, reax_list* );
void Delete_List( reax_list* );
inline int Num_Entries(int,reax_list*);
inline int Start_Index( int, reax_list* );

126
src/USER-REAXC/reaxc_lookup.cpp
View File

@ -50,19 +50,18 @@ void Tridiagonal_Solve( const double *a, const double *b,
}
void Natural_Cubic_Spline( const double *h, const double *f,
cubic_spline_coef *coef, unsigned int n,
MPI_Comm comm )
void Natural_Cubic_Spline( LAMMPS_NS::Error* error_ptr, const double *h, const double *f,
cubic_spline_coef *coef, unsigned int n )
{
int i;
double *a, *b, *c, *d, *v;
/* allocate space for the linear system */
a = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
b = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
c = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
d = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
v = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
a = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
b = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
c = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
d = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
v = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
/* build the linear system */
a[0] = a[1] = a[n-1] = 0;
@ -92,28 +91,27 @@ void Natural_Cubic_Spline( const double *h, const double *f,
coef[i-1].a = f[i];
}
sfree( a, "cubic_spline:a" );
sfree( b, "cubic_spline:b" );
sfree( c, "cubic_spline:c" );
sfree( d, "cubic_spline:d" );
sfree( v, "cubic_spline:v" );
sfree(error_ptr, a, "cubic_spline:a" );
sfree(error_ptr, b, "cubic_spline:b" );
sfree(error_ptr, c, "cubic_spline:c" );
sfree(error_ptr, d, "cubic_spline:d" );
sfree(error_ptr, v, "cubic_spline:v" );
}
void Complete_Cubic_Spline( const double *h, const double *f, double v0, double vlast,
cubic_spline_coef *coef, unsigned int n,
MPI_Comm comm )
void Complete_Cubic_Spline( LAMMPS_NS::Error* error_ptr, const double *h, const double *f, double v0, double vlast,
cubic_spline_coef *coef, unsigned int n )
{
int i;
double *a, *b, *c, *d, *v;
/* allocate space for the linear system */
a = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
b = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
c = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
d = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
v = (double*) smalloc( n * sizeof(double), "cubic_spline:a", comm );
a = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
b = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
c = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
d = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
v = (double*) smalloc(error_ptr, n * sizeof(double), "cubic_spline:a");
/* build the linear system */
a[0] = 0;
@ -142,11 +140,11 @@ void Complete_Cubic_Spline( const double *h, const double *f, double v0, double
coef[i-1].a = f[i];
}
sfree( a, "cubic_spline:a" );
sfree( b, "cubic_spline:b" );
sfree( c, "cubic_spline:c" );
sfree( d, "cubic_spline:d" );
sfree( v, "cubic_spline:v" );
sfree(error_ptr, a, "cubic_spline:a" );
sfree(error_ptr, b, "cubic_spline:b" );
sfree(error_ptr, c, "cubic_spline:c" );
sfree(error_ptr, d, "cubic_spline:d" );
sfree(error_ptr, v, "cubic_spline:v" );
}
@ -159,35 +157,33 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
double dr;
double *h, *fh, *fvdw, *fele, *fCEvd, *fCEclmb;
double v0_vdw, v0_ele, vlast_vdw, vlast_ele;
MPI_Comm comm;
/* initializations */
v0_vdw = 0;
v0_ele = 0;
vlast_vdw = 0;
vlast_ele = 0;
comm = mpi_data->world;
num_atom_types = system->reax_param.num_atom_types;
dr = control->nonb_cut / control->tabulate;
h = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:h", comm );
smalloc(system->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:h");
fh = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fh", comm );
smalloc(system->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fh");
fvdw = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fvdw", comm );
smalloc(system->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fvdw");
fCEvd = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fCEvd", comm );
smalloc(system->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fCEvd");
fele = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fele", comm );
smalloc(system->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fele");
fCEclmb = (double*)
smalloc( (control->tabulate+2) * sizeof(double), "lookup:fCEclmb", comm );
smalloc(system->error_ptr, (control->tabulate+2) * sizeof(double), "lookup:fCEclmb");
LR = (LR_lookup_table**)
scalloc( num_atom_types, sizeof(LR_lookup_table*), "lookup:LR", comm );
scalloc(system->error_ptr, num_atom_types, sizeof(LR_lookup_table*), "lookup:LR");
for( i = 0; i < num_atom_types; ++i )
LR[i] = (LR_lookup_table*)
scalloc( num_atom_types, sizeof(LR_lookup_table), "lookup:LR[i]", comm );
scalloc(system->error_ptr, num_atom_types, sizeof(LR_lookup_table), "lookup:LR[i]");
for( i = 0; i < MAX_ATOM_TYPES; ++i )
existing_types[i] = 0;
@ -207,22 +203,18 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
LR[i][j].dx = dr;
LR[i][j].inv_dx = control->tabulate / control->nonb_cut;
LR[i][j].y = (LR_data*)
smalloc( LR[i][j].n * sizeof(LR_data), "lookup:LR[i,j].y", comm );
smalloc(system->error_ptr, LR[i][j].n * sizeof(LR_data), "lookup:LR[i,j].y");
LR[i][j].H = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].H" ,
comm );
smalloc(system->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].H");
LR[i][j].vdW = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].vdW",
comm);
smalloc(system->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].vdW");
LR[i][j].CEvd = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].CEvd",
comm);
smalloc(system->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].CEvd");
LR[i][j].ele = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].ele",
comm );
smalloc(system->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),"lookup:LR[i,j].ele");
LR[i][j].CEclmb = (cubic_spline_coef*)
smalloc( LR[i][j].n*sizeof(cubic_spline_coef),
"lookup:LR[i,j].CEclmb", comm );
smalloc(system->error_ptr, LR[i][j].n*sizeof(cubic_spline_coef),
"lookup:LR[i,j].CEclmb");
for( r = 1; r <= control->tabulate; ++r ) {
LR_vdW_Coulomb( system, workspace, control, i, j, r * dr, &(LR[i][j].y[r]) );
@ -246,24 +238,20 @@ int Init_Lookup_Tables( reax_system *system, control_params *control,
vlast_vdw = fCEvd[r-1];
vlast_ele = fele[r-1];
Natural_Cubic_Spline( &h[1], &fh[1],
&(LR[i][j].H[1]), control->tabulate+1, comm );
Natural_Cubic_Spline( control->error_ptr, &h[1], &fh[1],
&(LR[i][j].H[1]), control->tabulate+1);
Complete_Cubic_Spline( &h[1], &fvdw[1], v0_vdw, vlast_vdw,
&(LR[i][j].vdW[1]), control->tabulate+1,
comm );
Complete_Cubic_Spline( control->error_ptr, &h[1], &fvdw[1], v0_vdw, vlast_vdw,
&(LR[i][j].vdW[1]), control->tabulate+1);
Natural_Cubic_Spline( &h[1], &fCEvd[1],
&(LR[i][j].CEvd[1]), control->tabulate+1,
comm );
Natural_Cubic_Spline( control->error_ptr, &h[1], &fCEvd[1],
&(LR[i][j].CEvd[1]), control->tabulate+1);
Complete_Cubic_Spline( &h[1], &fele[1], v0_ele, vlast_ele,
&(LR[i][j].ele[1]), control->tabulate+1,
comm );
Complete_Cubic_Spline( control->error_ptr, &h[1], &fele[1], v0_ele, vlast_ele,
&(LR[i][j].ele[1]), control->tabulate+1);
Natural_Cubic_Spline( &h[1], &fCEclmb[1],
&(LR[i][j].CEclmb[1]), control->tabulate+1,
comm );
Natural_Cubic_Spline( control->error_ptr, &h[1], &fCEclmb[1],
&(LR[i][j].CEclmb[1]), control->tabulate+1);
} else {
LR[i][j].n = 0;
}
@ -291,14 +279,14 @@ void Deallocate_Lookup_Tables( reax_system *system )
for( i = 0; i < ntypes; ++i ) {
for( j = i; j < ntypes; ++j )
if (LR[i][j].n) {
sfree( LR[i][j].y, "LR[i,j].y" );
sfree( LR[i][j].H, "LR[i,j].H" );
sfree( LR[i][j].vdW, "LR[i,j].vdW" );
sfree( LR[i][j].CEvd, "LR[i,j].CEvd" );
sfree( LR[i][j].ele, "LR[i,j].ele" );
sfree( LR[i][j].CEclmb, "LR[i,j].CEclmb" );
sfree(system->error_ptr, LR[i][j].y, "LR[i,j].y" );
sfree(system->error_ptr, LR[i][j].H, "LR[i,j].H" );
sfree(system->error_ptr, LR[i][j].vdW, "LR[i,j].vdW" );
sfree(system->error_ptr, LR[i][j].CEvd, "LR[i,j].CEvd" );
sfree(system->error_ptr, LR[i][j].ele, "LR[i,j].ele" );
sfree(system->error_ptr, LR[i][j].CEclmb, "LR[i,j].CEclmb" );
}
sfree( LR[i], "LR[i]" );
sfree(system->error_ptr, LR[i], "LR[i]" );
}
sfree( LR, "LR" );
sfree(system->error_ptr, LR, "LR" );
}

10
src/USER-REAXC/reaxc_lookup.h
View File

@ -32,13 +32,11 @@
void Tridiagonal_Solve( const double *a, const double *b,
double *c, double *d, double *x, unsigned int n);
void Natural_Cubic_Spline( const double *h, const double *f,
cubic_spline_coef *coef, unsigned int n,
MPI_Comm comm );
void Natural_Cubic_Spline( LAMMPS_NS::Error*, const double *h, const double *f,
cubic_spline_coef *coef, unsigned int n );
void Complete_Cubic_Spline( const double *h, const double *f, double v0, double vlast,
cubic_spline_coef *coef, unsigned int n,
MPI_Comm comm );
void Complete_Cubic_Spline( LAMMPS_NS::Error*, const double *h, const double *f, double v0, double vlast,
cubic_spline_coef *coef, unsigned int n );
int Init_Lookup_Tables( reax_system*, control_params*, storage*,
mpi_datatypes*, char* );

26
src/USER-REAXC/reaxc_reset_tools.cpp
View File

@ -30,6 +30,9 @@
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
#include "error.h"
void Reset_Atoms( reax_system* system, control_params *control )
{
int i;
@ -120,8 +123,7 @@ void Reset_Workspace( reax_system *system, storage *workspace )
void Reset_Neighbor_Lists( reax_system *system, control_params *control,
storage *workspace, reax_list **lists,
MPI_Comm comm )
storage *workspace, reax_list **lists )
{
int i, total_bonds, Hindex, total_hbonds;
reax_list *bonds, *hbonds;
@ -142,10 +144,10 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
if (total_bonds >= bonds->num_intrs * DANGER_ZONE) {
workspace->realloc.bonds = 1;
if (total_bonds >= bonds->num_intrs) {
fprintf(stderr,
"p%d: not enough space for bonds! total=%d allocated=%d\n",
system->my_rank, total_bonds, bonds->num_intrs );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
char errmsg[256];
snprintf(errmsg, 256, "Not enough space for bonds! total=%d allocated=%d\n",
total_bonds, bonds->num_intrs);
control->error_ptr->one(FLERR, errmsg);
}
}
}
@ -168,10 +170,10 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
if (total_hbonds >= hbonds->num_intrs * 0.90/*DANGER_ZONE*/) {
workspace->realloc.hbonds = 1;
if (total_hbonds >= hbonds->num_intrs) {
fprintf(stderr,
"p%d: not enough space for hbonds! total=%d allocated=%d\n",
system->my_rank, total_hbonds, hbonds->num_intrs );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
char errmsg[256];
snprintf(errmsg, 256, "Not enough space for hbonds! total=%d allocated=%d\n",
total_hbonds, hbonds->num_intrs);
control->error_ptr->one(FLERR, errmsg);
}
}
}
@ -179,7 +181,7 @@ void Reset_Neighbor_Lists( reax_system *system, control_params *control,
void Reset( reax_system *system, control_params *control, simulation_data *data,
storage *workspace, reax_list **lists, MPI_Comm comm )
storage *workspace, reax_list **lists )
{
Reset_Atoms( system, control );
@ -187,6 +189,6 @@ void Reset( reax_system *system, control_params *control, simulation_data *data,
Reset_Workspace( system, workspace );
Reset_Neighbor_Lists( system, control, workspace, lists, comm );
Reset_Neighbor_Lists( system, control, workspace, lists );
}

4
src/USER-REAXC/reaxc_reset_tools.h
View File

@ -34,7 +34,7 @@ void Reset_Simulation_Data( simulation_data*, int );
void Reset_Timing( reax_timing* );
void Reset_Workspace( reax_system*, storage* );
void Reset_Neighbor_Lists( reax_system*, control_params*, storage*,
reax_list**, MPI_Comm );
reax_list** );
void Reset( reax_system*, control_params*, simulation_data*, storage*,
reax_list**, MPI_Comm );
reax_list** );
#endif

46
src/USER-REAXC/reaxc_tool_box.cpp
View File

@ -27,6 +27,8 @@
#include "pair_reaxc.h"
#include "reaxc_tool_box.h"
#include "error.h"
struct timeval tim;
double t_end;
@ -53,23 +55,25 @@ int Tokenize( char* s, char*** tok )
return count;
}
/* safe malloc */
void *smalloc( rc_bigint n, const char *name, MPI_Comm comm )
void *smalloc( LAMMPS_NS::Error *error_ptr, rc_bigint n, const char *name )
{
void *ptr;
char errmsg[256];
if (n <= 0) {
fprintf( stderr, "WARNING: trying to allocate %ld bytes for array %s. ",
n, name );
fprintf( stderr, "returning NULL.\n" );
snprintf(errmsg, 256, "Trying to allocate %ld bytes for array %s. "
"returning NULL.", n, name);
error_ptr->one(FLERR,errmsg);
return NULL;
}
ptr = malloc( n );
if (ptr == NULL) {
fprintf( stderr, "ERROR: failed to allocate %ld bytes for array %s",
n, name );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
snprintf(errmsg, 256, "Failed to allocate %ld bytes for array %s", n, name);
error_ptr->one(FLERR,errmsg);
}
return ptr;
@ -77,29 +81,30 @@ void *smalloc( rc_bigint n, const char *name, MPI_Comm comm )
/* safe calloc */
void *scalloc( rc_bigint n, rc_bigint size, const char *name, MPI_Comm comm )
void *scalloc( LAMMPS_NS::Error *error_ptr, rc_bigint n, rc_bigint size, const char *name )
{
void *ptr;
char errmsg[256];
if (n <= 0) {
fprintf( stderr, "WARNING: trying to allocate %ld elements for array %s. ",
n, name );
fprintf( stderr, "returning NULL.\n" );
snprintf(errmsg, 256, "Trying to allocate %ld elements for array %s. "
"returning NULL.\n", n, name );
error_ptr->one(FLERR,errmsg);
return NULL;
}
if (size <= 0) {
fprintf( stderr, "WARNING: elements size for array %s is %ld. ",
name, size );
fprintf( stderr, "returning NULL.\n" );
snprintf(errmsg, 256, "Elements size for array %s is %ld. "
"returning NULL", name, size );
error_ptr->one(FLERR,errmsg);
return NULL;
}
ptr = calloc( n, size );
if (ptr == NULL) {
fprintf( stderr, "ERROR: failed to allocate %ld bytes for array %s",
n*size, name );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
char errmsg[256];
snprintf(errmsg, 256, "Failed to allocate %ld bytes for array %s", n*size, name);
error_ptr->one(FLERR,errmsg);
}
return ptr;
@ -107,11 +112,12 @@ void *scalloc( rc_bigint n, rc_bigint size, const char *name, MPI_Comm comm )
/* safe free */
void sfree( void *ptr, const char *name )
void sfree( LAMMPS_NS::Error* error_ptr, void *ptr, const char *name )
{
if (ptr == NULL) {
fprintf( stderr, "WARNING: trying to free the already NULL pointer %s!\n",
name );
char errmsg[256];
snprintf(errmsg, 256, "Trying to free the already NULL pointer %s", name );
error_ptr->one(FLERR,errmsg);
return;
}

6
src/USER-REAXC/reaxc_tool_box.h
View File

@ -37,7 +37,7 @@ double Get_Time( );
int Tokenize( char*, char*** );
/* from lammps */
void *smalloc( rc_bigint, const char*, MPI_Comm );
void *scalloc( rc_bigint, rc_bigint, const char*, MPI_Comm );
void sfree( void*, const char* );
void *smalloc( LAMMPS_NS::Error*, rc_bigint, const char* );
void *scalloc( LAMMPS_NS::Error*, rc_bigint, rc_bigint, const char* );
void sfree( LAMMPS_NS::Error*, void*, const char* );
#endif

31
src/USER-REAXC/reaxc_traj.cpp
View File

@ -29,8 +29,9 @@
#include "reaxc_list.h"
#include "reaxc_tool_box.h"
int Reallocate_Output_Buffer( output_controls *out_control, int req_space,
MPI_Comm comm )
#include "error.h"
int Reallocate_Output_Buffer( LAMMPS_NS::Error *error_ptr, output_controls *out_control, int req_space )
{
if (out_control->buffer_len > 0)
free( out_control->buffer );
@ -38,10 +39,9 @@ int Reallocate_Output_Buffer( output_controls *out_control, int req_space,
out_control->buffer_len = (int)(req_space*SAFE_ZONE);
out_control->buffer = (char*) malloc(out_control->buffer_len*sizeof(char));
if (out_control->buffer == NULL) {
fprintf( stderr,
"insufficient memory for required buffer size %d. terminating!\n",
(int) (req_space*SAFE_ZONE) );
MPI_Abort( comm, INSUFFICIENT_MEMORY );
char errmsg[256];
snprintf(errmsg, 256, "Insufficient memory for required buffer size %d", (int) (req_space*SAFE_ZONE));
error_ptr->one(FLERR,errmsg);
}
return SUCCESS;
@ -83,7 +83,7 @@ int Write_Header( reax_system *system, control_params *control,
my_hdr_lines = num_hdr_lines * ( system->my_rank == MASTER_NODE );
buffer_req = my_hdr_lines * HEADER_LINE_LEN;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
Reallocate_Output_Buffer( control->error_ptr, out_control, buffer_req );
/* only the master node writes into trajectory header */
if (system->my_rank == MASTER_NODE) {
@ -278,7 +278,7 @@ int Write_Init_Desc( reax_system *system, control_params * /*control*/,
else buffer_req = system->n * INIT_DESC_LEN + 1;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
Reallocate_Output_Buffer( system->error_ptr, out_control, buffer_req );
out_control->line[0] = 0;
out_control->buffer[0] = 0;
@ -367,7 +367,7 @@ int Write_Frame_Header( reax_system *system, control_params *control,
my_frm_hdr_lines = num_frm_hdr_lines * ( me == MASTER_NODE );
buffer_req = my_frm_hdr_lines * HEADER_LINE_LEN;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
Reallocate_Output_Buffer( control->error_ptr, out_control, buffer_req );
/* only the master node writes into trajectory header */
if (me == MASTER_NODE) {
@ -500,7 +500,7 @@ int Write_Atoms( reax_system *system, control_params * /*control*/,
else buffer_req = system->n * line_len + 1;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
Reallocate_Output_Buffer( system->error_ptr, out_control, buffer_req );
/* fill in buffer */
out_control->line[0] = 0;
@ -531,9 +531,7 @@ int Write_Atoms( reax_system *system, control_params * /*control*/,
p_atom->f[0], p_atom->f[1], p_atom->f[2], p_atom->q );
break;
default:
fprintf( stderr,
"write_traj_atoms: unknown atom trajectroy format!\n");
MPI_Abort( mpi_data->world, UNKNOWN_OPTION );
system->error_ptr->one(FLERR,"Write_traj_atoms: unknown atom trajectory format");
}
strncpy( out_control->buffer + i*line_len, out_control->line, line_len+1 );
@ -592,7 +590,7 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
else buffer_req = my_bonds * line_len + 1;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
Reallocate_Output_Buffer( system->error_ptr, out_control, buffer_req );
/* fill in the buffer */
out_control->line[0] = 0;
@ -619,8 +617,7 @@ int Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
bo_ij->bo_data.BO_pi, bo_ij->bo_data.BO_pi2 );
break;
default:
fprintf(stderr, "write_traj_bonds: FATAL! invalid bond_info option");
MPI_Abort( mpi_data->world, UNKNOWN_OPTION );
system->error_ptr->one(FLERR, "Write_traj_bonds: FATAL! invalid bond_info option");
}
strncpy( out_control->buffer + my_bonds*line_len,
out_control->line, line_len+1 );
@ -693,7 +690,7 @@ int Write_Angles( reax_system *system, control_params *control,
else buffer_req = my_angles * line_len + 1;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer( out_control, buffer_req, mpi_data->world );
Reallocate_Output_Buffer( system->error_ptr, out_control, buffer_req );
/* fill in the buffer */
my_angles = 0;

7
src/USER-REAXC/reaxc_types.h
View File

@ -39,6 +39,9 @@
#include <sys/time.h>
#include "accelerator_kokkos.h"
namespace LAMMPS_NS { class Error;}
#if defined LMP_USER_OMP
#define OMP_TIMING 0
@ -409,6 +412,7 @@ struct _reax_system
boundary_cutoff bndry_cuts;
reax_atom *my_atoms;
class LAMMPS_NS::Error *error_ptr;
class LAMMPS_NS::Pair *pair_ptr;
int my_bonds;
int mincap;
@ -486,6 +490,8 @@ typedef struct
int lgflag;
int enobondsflag;
class LAMMPS_NS::Error *error_ptr;
int me;
} control_params;
@ -772,6 +778,7 @@ struct _reax_list
int type;
list_type select;
class LAMMPS_NS::Error *error_ptr;
};
typedef _reax_list reax_list;

9
src/USER-REAXC/reaxc_valence_angles.cpp
View File

@ -30,6 +30,8 @@
#include "reaxc_list.h"
#include "reaxc_vector.h"
#include "error.h"
static double Dot( double* v1, double* v2, int k )
{
double ret = 0.0;
@ -405,9 +407,10 @@ void Valence_Angles( reax_system *system, control_params *control,
if (num_thb_intrs >= thb_intrs->num_intrs * DANGER_ZONE) {
workspace->realloc.num_3body = num_thb_intrs;
if (num_thb_intrs > thb_intrs->num_intrs) {
fprintf( stderr, "step%d-ran out of space on angle_list: top=%d, max=%d",
data->step, num_thb_intrs, thb_intrs->num_intrs );
MPI_Abort( MPI_COMM_WORLD, INSUFFICIENT_MEMORY );
char errmsg[128];
snprintf(errmsg, 128, "step%d-ran out of space on angle_list: top=%d, max=%d",
data->step, num_thb_intrs, thb_intrs->num_intrs);
control->error_ptr->one(FLERR, errmsg);
}
}