diff --git a/PuReMD/Makefile.am b/PuReMD/Makefile.am
index 192ad0e9a1457c9e781fb94c941ead49eb54c1c2..718c453fa84b45f83c4c3890eda1735fbb0b7f20 100644
--- a/PuReMD/Makefile.am
+++ b/PuReMD/Makefile.am
@@ -16,5 +16,5 @@ include_HEADERS = src/reax_defs.h src/reax_types.h \
src/bonds.h src/valence_angles.h src/hydrogen_bonds.h src/torsion_angles.h src/nonbonded.h src/forces.h \
src/integrate.h src/init_md.h
-bin_puremd_CFLAGS = $(AM_CFLAGS) $(MPI_CFLAGS) $(CFLAGS)
-bin_puremd_LDADD = $(AM_LDADD) $(MPI_LIBS) $(LDADD)
+bin_puremd_CFLAGS = $(CFLAGS) $(MPI_CFLAGS)
+bin_puremd_LDADD = $(LDADD) $(MPI_LIBS)
diff --git a/PuReMD/configure.ac b/PuReMD/configure.ac
index 0168249747115f2cd2ce99926f61a21c7d179a4c..8833fc248a64e49a684106feeaeb4682c54bf715 100644
--- a/PuReMD/configure.ac
+++ b/PuReMD/configure.ac
@@ -5,7 +5,7 @@ AC_PREREQ([2.69])
AC_INIT([PuReMD], [1.0], [ohearnku@msu.edu hma@msu.edu])
# Do not allow AC_PROG_CC to set CFLAGS (this line must be after AC_INIT but before AC_PROG_CC)
-sav_CFLAGS="${CFLAGS}"
+save_CFLAGS="${CFLAGS}"
: ${CFLAGS=""}
AM_INIT_AUTOMAKE([1.15 subdir-objects -Wall -Werror foreign])
# Enable silent build rules by default.
@@ -26,7 +26,7 @@ AC_LANG([C])
# Checks for programs.
AC_PROG_CC([cc icc gcc])
AC_PROG_CPP
-CFLAGS="${sav_CFLAGS}"
+CFLAGS="${save_CFLAGS}"
AC_CONFIG_SRCDIR([src/torsion_angles.h])
AC_CONFIG_HEADERS([src/config.h])
@@ -50,32 +50,37 @@ AC_FUNC_REALLOC
AC_FUNC_STRTOD
AC_CHECK_FUNCS([gettimeofday memset pow sqrt])
-# Check for compiler vendor
+# Check for compiler vendor. If the compiler is recognized,
+# the variable ax_cv_c_compiler_vendor is set accordingly.
AX_COMPILER_VENDOR
if test "x${ax_cv_c_compiler_vendor}" = "xgnu"; then
if test "x${BUILD_DEBUG}" = "x"; then
- CFLAGS="${CFLAGS} -Wall -O3 -funroll-loops"
+ CFLAGS="-Wall -O2 -funroll-loops"
else
- CFLAGS="${CFLAGS} -Wall -O0 -g"
+ CFLAGS="-Wall -O0 -g2 -D_GLIBCXX_DEBUG"
fi
-fi
-if test "x${ax_cv_c_compiler_vendor}" = "xintel"; then
+elif test "x${ax_cv_c_compiler_vendor}" = "xintel"; then
if test "x${BUILD_DEBUG}" = "x"; then
- CFLAGS="$CFLAGS -fast"
+ CFLAGS="-Wall -O2 -funroll-loops"
+ else
+ CFLAGS="-Wall -O0"
fi
fi
-# Check for MPI support.
+# Check for MPI support. If found, sets the following variables:
+# Compiler-related: MPICC, MPICXX, MPIF77
+# Library linking: MPILIBS
+# Compilation support (via preprocessor def): HAVE_MPI
CONFIGURE_HEADLINE([ MPI compiler ])
ACX_MPI([], [AC_MSG_ERROR([could not find mpi library])])
AC_CHECK_PROG(MPIRUN, mpirun, mpirun)
AC_SUBST(MPIRUN)
# try to find if we are using OpenMPI / MPICH by looking inside mpi.h
-sav_CC="${CC}"
-sav_CFLAGS="${CFLAGS}"
+save_CC="${CC}"
+save_CFLAGS="${CFLAGS}"
CC="${MPICC}"
-CFLAGS="${CFLAGS}"
+CFLAGS=""
AC_CHECK_DECL([OPEN_MPI], [mpi_vendor="OpenMPI"],
[], [#include "mpi.h"])
# MPICH v2
@@ -84,8 +89,8 @@ AC_CHECK_DECL([MPICH2], [mpi_vendor="MPICH2"],
# MPICH v3
AC_CHECK_DECL([MPICH_VERSION], [mpi_vendor="MPICH3"],
[], [#include "mpi.h"])
-CC="${sav_CC}"
-CFLAGS="${sav_CFLAGS}"
+CC="${save_CC}"
+CFLAGS="${save_CFLAGS}"
# try to set MPI_CFLAGS and MPI_LIBS
MPI_CFLAGS=
@@ -104,7 +109,7 @@ then
MPI_CFLAGS=
for i in $tmp
do
- case $i in
+ case ${i} in
-[[DIUbi]]*)
MPI_CFLAGS="${MPI_CFLAGS} ${i}"
;;
@@ -114,7 +119,7 @@ then
tmp=`${MPICC} -link-info | awk '{$1=""; print $0 }'`
for i in $tmp
do
- case $i in
+ case ${i} in
[[\\/]]*.a | ?:[[\\/]]*.a | -[[lLRu]]* | -Wl* )
MPI_LIBS="${MPI_LIBS} ${i}"
;;
@@ -138,7 +143,10 @@ then
[-lmkl_sequential -lmkl_core -lpthread -lm -ldl])
AS_IF([test "x${MKL_FOUND_LIBS}" != "xyes"],
[AC_MSG_ERROR([Unable to find MKL LAPACKE library.])])
- LIBS="${LIBS} -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl"
+ # dynamic linking
+ #LIBS="${LIBS} -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl"
+ # static linking
+ LIBS="${LIBS} -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_ilp64.a ${MKLROOT}/lib/intel64/libmkl_sequential.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl"
AC_DEFINE([HAVE_LAPACKE_MKL], [1], [Define to 1 if you have MKL LAPACKE support enabled.])
else
AC_CHECK_HEADERS([lapacke.h], [LAPACKE_FOUND_HEADERS="yes"])
diff --git a/PuReMD/src/allocate.c b/PuReMD/src/allocate.c
index 08f455866404bbd5d4f8ab93f5a513127f684e15..2227656be245149b527a6d7b0968c0975b150fde 100644
--- a/PuReMD/src/allocate.c
+++ b/PuReMD/src/allocate.c
@@ -709,7 +709,8 @@ void Deallocate_Grid( grid *g )
buffers are void*, type cast to the correct pointer type to access
the allocated buffers */
int Allocate_MPI_Buffers( mpi_datatypes *mpi_data, int est_recv,
- neighbor_proc *my_nbrs, char *msg )
+ neighbor_proc *my_nbrs,
+ char *msg )
{
int i;
mpi_out_data *mpi_buf;
@@ -734,6 +735,19 @@ int Allocate_MPI_Buffers( mpi_datatypes *mpi_data, int est_recv,
scalloc( my_nbrs[i].est_send, sizeof(boundary_atom), "mpibuf:out_atoms",
comm );
}
+#if defined(NEUTRAL_TERRITORY)
+ /* Neutral Territory out buffers */
+ for ( i = 0; i < MAX_NT_NBRS; ++i )
+ {
+ mpi_buf = &( mpi_data->out_nt_buffers[i] );
+ /* allocate storage for the neighbor processor i */
+ mpi_buf->index = (int*)
+ scalloc( my_nbrs[i].est_send, sizeof(int), "mpibuf:nt_index", comm );
+ mpi_buf->out_atoms = (void*)
+ scalloc( my_nbrs[i].est_send, sizeof(boundary_atom), "mpibuf:nt_out_atoms",
+ comm );
+ }
+#endif
return SUCCESS;
}
@@ -753,6 +767,14 @@ void Deallocate_MPI_Buffers( mpi_datatypes *mpi_data )
sfree( mpi_buf->index, "mpibuf:index" );
sfree( mpi_buf->out_atoms, "mpibuf:out_atoms" );
}
+#if defined(NEUTRAL_TERRITORY)
+ for ( i = 0; i < MAX_NT_NBRS; ++i )
+ {
+ mpi_buf = &( mpi_data->out_nt_buffers[i] );
+ sfree( mpi_buf->index, "mpibuf:nt_index" );
+ sfree( mpi_buf->out_atoms, "mpibuf:nt_out_atoms" );
+ }
+#endif
}
@@ -1004,6 +1026,20 @@ void ReAllocate( reax_system *system, control_params *control,
break;
}
}
+#if defined(NEUTRAL_TERRITORY)
+ // otherwise check individual outgoing Neutral Territory buffers
+ mpi_flag = 0;
+ for ( p = 0; p < MAX_NT_NBRS; ++p )
+ {
+ nbr_pr = &( system->my_nt_nbrs[p] );
+ nbr_data = &( mpi_data->out_nt_buffers[p] );
+ if ( nbr_data->cnt >= nbr_pr->est_send * 0.90 )
+ {
+ mpi_flag = 1;
+ break;
+ }
+ }
+#endif
}
if ( mpi_flag )
@@ -1020,6 +1056,7 @@ void ReAllocate( reax_system *system, control_params *control,
system->est_trans =
(system->est_recv * sizeof(boundary_atom)) / sizeof(mpi_atom);
total_send = 0;
+
for ( p = 0; p < MAX_NBRS; ++p )
{
nbr_pr = &( system->my_nbrs[p] );
@@ -1027,6 +1064,15 @@ void ReAllocate( reax_system *system, control_params *control,
nbr_pr->est_send = MAX( nbr_data->cnt * SAFER_ZONE, MIN_SEND );
total_send += nbr_pr->est_send;
}
+#if defined(NEUTRAL_TERRITORY)
+ for ( p = 0; p < MAX_NT_NBRS; ++p )
+ {
+ nbr_pr = &( system->my_nt_nbrs[p] );
+ nbr_data = &( mpi_data->out_nt_buffers[p] );
+ nbr_pr->est_send = MAX( nbr_data->cnt * SAFER_ZONE, MIN_SEND );
+ }
+#endif
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: reallocating mpi_buf: recv=%d send=%d total=%dMB\n",
system->my_rank, system->est_recv, total_send,
diff --git a/PuReMD/src/basic_comm.c b/PuReMD/src/basic_comm.c
index 3014299478da41545f9bdd0da2efc5f4e379d7db..135abb6e49404c8305938ea510e5c0bb2aafd9c6 100644
--- a/PuReMD/src/basic_comm.c
+++ b/PuReMD/src/basic_comm.c
@@ -129,6 +129,60 @@ void Dist( reax_system* system, mpi_datatypes *mpi_data,
}
+void Dist_NT( reax_system* system, mpi_datatypes *mpi_data,
+ void *buf, MPI_Datatype type, int scale, dist_packer pack )
+{
+ int d;
+ mpi_out_data *out_bufs;
+ MPI_Comm comm;
+ MPI_Request req[6];
+ MPI_Status stat[6];
+ neighbor_proc *nbr;
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d dist: entered\n", system->my_rank );
+#endif
+ comm = mpi_data->comm_mesh3D;
+ out_bufs = mpi_data->out_nt_buffers;
+
+ /* initiate recvs */
+ for( d = 0; d < 6; ++d )
+ {
+ nbr = &(system->my_nt_nbrs[d]);
+ if ( nbr->atoms_cnt )
+ {
+ MPI_Irecv( buf + nbr->atoms_str * scale, nbr->atoms_cnt, type,
+ nbr->receive_rank, d, comm, &(req[d]) );
+ }
+ }
+
+ for( d = 0; d < 6; ++d)
+ {
+ /* send both messages in dimension d */
+ if ( out_bufs[d].cnt )
+ {
+ pack( buf, out_bufs + d );
+ MPI_Send( out_bufs[d].out_atoms, out_bufs[d].cnt, type,
+ nbr->rank, d, comm );
+ }
+ }
+
+ for( d = 0; d < 6; ++d )
+ {
+ nbr = &(system->my_nbrs[d]);
+ if ( nbr->atoms_cnt )
+ {
+ // TODO: Wait(MPI_WAITANY)
+ MPI_Wait( &(req[d]), &(stat[d]) );
+ }
+ }
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d dist: done\n", system->my_rank );
+#endif
+}
+
+
void real_unpacker( void *dummy_in, void *dummy_buf, mpi_out_data *out_buf )
{
int i;
@@ -236,6 +290,63 @@ void Coll( reax_system* system, mpi_datatypes *mpi_data,
fprintf( stderr, "p%d coll: done\n", system->my_rank );
#endif
}
+
+
+void Coll_NT( reax_system* system, mpi_datatypes *mpi_data,
+ void *buf, MPI_Datatype type, int scale, coll_unpacker unpack )
+{
+ int d;
+ void *in[6];
+ mpi_out_data *out_bufs;
+ MPI_Comm comm;
+ MPI_Request req[6];
+ MPI_Status stat[6];
+ neighbor_proc *nbr;
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d coll: entered\n", system->my_rank );
+#endif
+ comm = mpi_data->comm_mesh3D;
+ out_bufs = mpi_data->out_buffers;
+
+ for ( d = 0; d < 6; ++d )
+ {
+ /* initiate recvs */
+ nbr = &(system->my_nbrs[d]);
+ in[d] = mpi_data->in_nt_buffer[d];
+ if ( out_bufs[d].cnt )
+ {
+ MPI_Irecv(in[d], out_bufs[d].cnt, type, nbr->receive_rank, d, comm, &(req[d]));
+ }
+ }
+
+ for( d = 0; d < 6; ++d )
+ {
+ /* send both messages in direction d */
+ nbr = &(system->my_nbrs[d]);
+ if ( nbr->atoms_cnt )
+ {
+ MPI_Send( buf + nbr->atoms_str * scale, nbr->atoms_cnt, type,
+ nbr->rank, d, comm );
+ }
+ }
+
+ for( d = 0; d < 6; ++d )
+ {
+ if ( out_bufs[d].cnt )
+ {
+ //TODO: WAITANY
+ MPI_Wait( &(req[d]), &(stat[d]));
+ unpack( in[d], buf, out_bufs + d );
+ }
+ }
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d coll: done\n", system->my_rank );
+#endif
+}
+
+
#endif /*PURE_REAX*/
/*****************************************************************************/
diff --git a/PuReMD/src/comm_tools.c b/PuReMD/src/comm_tools.c
index 8419e3efe44377e7d1680fe09b2bc06ef5756ca6..0b399825a58054fdfa6aa07f5351d4f221f45c55 100644
--- a/PuReMD/src/comm_tools.c
+++ b/PuReMD/src/comm_tools.c
@@ -26,6 +26,143 @@
#include "vector.h"
+void Setup_NT_Comm( reax_system* system, control_params* control,
+ mpi_datatypes *mpi_data )
+{
+ int i, d;
+ real bndry_cut;
+ neighbor_proc *nbr_pr;
+ simulation_box *my_box;
+ ivec nbr_coords;
+ ivec r[12] = {
+ {0, 0, -1}, // -z
+ {0, 0, +1}, // +z
+ {0, -1, 0}, // -y
+ {-1, -1, 0}, // -x-y
+ {-1, 0, 0}, // -x
+ {-1, +1, 0}, // -x+y
+
+ {0, 0, +1}, // +z
+ {0, 0, -1}, // -z
+ {0, +1, 0}, // +y
+ {+1, +1, 0}, // +x+y
+ {+1, 0, 0}, // -x
+ {+1, -1, 0} // -x+y
+ };
+ my_box = &(system->my_box);
+ bndry_cut = system->bndry_cuts.ghost_cutoff;
+
+ /* identify my neighbors */
+ system->num_nt_nbrs = MAX_NT_NBRS;
+ for ( i = 0; i < system->num_nt_nbrs; ++i )
+ {
+ ivec_Sum( nbr_coords, system->my_coords, r[i] ); /* actual nbr coords */
+ nbr_pr = &(system->my_nt_nbrs[i]);
+ MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_coords, &(nbr_pr->rank) );
+
+ /* set the rank of the neighbor processor in the receiving direction */
+ ivec_Sum( nbr_coords, system->my_coords, r[i + 6] ); /* actual nbr coords */
+ MPI_Cart_rank( mpi_data->comm_mesh3D, nbr_coords, &(nbr_pr->receive_rank) );
+
+ for ( d = 0; d < 3; ++d )
+ {
+ /* determine the boundary area with this nbr */
+ if ( r[i][d] < 0 )
+ {
+ nbr_pr->bndry_min[d] = my_box->min[d];
+ nbr_pr->bndry_max[d] = my_box->min[d] + bndry_cut;
+ }
+ else if ( r[i][d] > 0 )
+ {
+ nbr_pr->bndry_min[d] = my_box->max[d] - bndry_cut;
+ nbr_pr->bndry_max[d] = my_box->max[d];
+ }
+ else
+ {
+ nbr_pr->bndry_min[d] = my_box->min[d];
+ nbr_pr->bndry_max[d] = my_box->max[d];
+ }
+
+ /* determine if it is a periodic neighbor */
+ if ( nbr_coords[d] < 0 )
+ nbr_pr->prdc[d] = -1;
+ else if ( nbr_coords[d] >= control->procs_by_dim[d] )
+ nbr_pr->prdc[d] = +1;
+ else
+ nbr_pr->prdc[d] = 0;
+ }
+
+ }
+}
+
+
+void Sort_Neutral_Territory( reax_system *system, int dir, mpi_out_data *out_bufs )
+{
+ int i, d, p, out_cnt;
+ reax_atom *atoms;
+ simulation_box *my_box;
+ boundary_atom *out_buf;
+ neighbor_proc *nbr_pr;
+
+ atoms = system->my_atoms;
+ my_box = &( system->my_box );
+
+ /* place each atom into the appropriate outgoing list */
+ nbr_pr = &( system->my_nt_nbrs[dir] );
+ for ( i = 0; i < system->n; ++i )
+ {
+ if ( nbr_pr->bndry_min[0] <= atoms[i].x[0] &&
+ atoms[i].x[0] < nbr_pr->bndry_max[0] &&
+ nbr_pr->bndry_min[1] <= atoms[i].x[1] &&
+ atoms[i].x[1] < nbr_pr->bndry_max[1] &&
+ nbr_pr->bndry_min[2] <= atoms[i].x[2] &&
+ atoms[i].x[2] < nbr_pr->bndry_max[2] )
+ {
+ out_bufs[dir].index[out_cnt] = i;
+ out_bufs[dir].cnt++;
+ }
+ }
+}
+
+
+void Init_Neutral_Territory( reax_system* system, mpi_datatypes *mpi_data )
+{
+ int d, start, end, cnt;
+ mpi_out_data *out_bufs;
+ void *in;
+ MPI_Comm comm;
+ MPI_Request req;
+ MPI_Status stat;
+ neighbor_proc *nbr;
+
+ Reset_Out_Buffers( mpi_data->out_nt_buffers, system->num_nt_nbrs );
+ comm = mpi_data->comm_mesh3D;
+ out_bufs = mpi_data->out_nt_buffers;
+ cnt = 0;
+ end = system->n;
+
+ for ( d = 0; d < 6; ++d )
+ {
+ Sort_Neutral_Territory( system, d, out_bufs );
+ start = end;
+
+ MPI_Irecv( &cnt, 1, MPI_INT, nbr->receive_rank, d, comm, &req );
+ MPI_Send( &(out_bufs[d].cnt), 1, MPI_INT, nbr->rank, d, comm );
+ MPI_Wait( &req, &stat );
+
+ if( mpi_data->in_nt_buffer[d] == NULL )
+ {
+ mpi_data->in_nt_buffer[d] = (void *) smalloc( SAFER_ZONE * out_bufs[d].cnt * sizeof(real), "in", comm );
+ }
+
+ nbr = &(system->my_nt_nbrs[d]);
+ nbr->atoms_str = end;
+ nbr->atoms_cnt = cnt;
+ end += cnt;
+ }
+}
+
+
void Setup_Comm( reax_system* system, control_params* control,
mpi_datatypes *mpi_data )
{
@@ -653,6 +790,10 @@ void Comm_Atoms( reax_system *system, control_params *control,
#endif
Bin_Boundary_Atoms( system );
+
+#if defined(NEUTRAL_TERRITORY)
+ Init_Neutral_Territory( system, mpi_data );
+#endif
}
else
{
diff --git a/PuReMD/src/ffield.c b/PuReMD/src/ffield.c
index b05216bdcb19f8002bfe02293a3cf6e28dd4faa7..b25b8db44ab7d20d1380069533e9a65bc82db55b 100644
--- a/PuReMD/src/ffield.c
+++ b/PuReMD/src/ffield.c
@@ -43,11 +43,7 @@ char Read_Force_Field( char *ffield_file, reax_interaction *reax,
comm = MPI_COMM_WORLD;
/* open force field file */
- if ( (fp = fopen( ffield_file, "r" ) ) == NULL )
- {
- fprintf( stderr, "error opening the force filed file! terminating...\n" );
- MPI_Abort( comm, FILE_NOT_FOUND );
- }
+ fp = sfopen( ffield_file, "r", "Read_Force_Field::fp" );
s = (char*) malloc(sizeof(char) * MAX_LINE);
tmp = (char**) malloc(sizeof(char*)*MAX_TOKENS);
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index 69ab40e79b1fbfb948fc95dd84a1998229aff3e8..30a169bcdbc1402cad9c4e65d17382b2836835cf 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -378,6 +378,7 @@ void Init_Forces( reax_system *system, control_params *control,
}
sbp_i = &(system->reax_param.sbp[type_i]);
+ //TODO: edit this part to include NT atoms
if ( i < system->n )
{
local = 1;
diff --git a/PuReMD/src/geo_tools.c b/PuReMD/src/geo_tools.c
index c01f65cbb63c3be0e2ad1c1f54844e01bfc93db3..77e1f95b32b7bf8e12b5de46578af8364fe505e3 100644
--- a/PuReMD/src/geo_tools.c
+++ b/PuReMD/src/geo_tools.c
@@ -81,11 +81,7 @@ char Read_Geo( char* geo_file, reax_system* system, control_params *control,
comm = MPI_COMM_WORLD;
/* open the geometry file */
- if ( (geo = fopen(geo_file, "r")) == NULL )
- {
- fprintf( stderr, "fopen: error opening the geo file! terminating...\n" );
- MPI_Abort( comm, FILE_NOT_FOUND );
- }
+ geo = sfopen( geo_file, "r", "Read_Geo::geo" );
/* read box information */
fscanf( geo, CUSTOM_BOXGEO_FORMAT,
@@ -140,7 +136,7 @@ char Read_Geo( char* geo_file, reax_system* system, control_params *control,
}
}
- fclose( geo );
+ sfclose( geo, "Read_Geo::geo" );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: finished reading the geo file\n", system->my_rank );
@@ -271,11 +267,7 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
comm = MPI_COMM_WORLD;
/* open pdb file */
- if ( (pdb = fopen(pdb_file, "r")) == NULL )
- {
- fprintf( stderr, "fopen: error opening the pdb file! terminating...\n" );
- MPI_Abort( comm, FILE_NOT_FOUND );
- }
+ pdb = sfopen( pdb_file, "r", "Read_PDB::pdb" );
/* allocate memory for tokenizing pdb lines */
if ( Allocate_Tokenizer_Space( &s, &s1, &tmp ) == FAILURE )
@@ -481,7 +473,7 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
return FAILURE;
}
- fclose( pdb );
+ sfclose( pdb, "Read_PDB::pdb" );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: finished reading the pdb file\n", system->my_rank );
@@ -550,7 +542,7 @@ char Write_PDB(reax_system* system, reax_list** bonds, simulation_data *data,
sprintf(fname, "%s-%d.pdb", control->sim_name, data->step);
- pdb = fopen(fname, "w");
+ pdb = sfopen( fname, "w", "Write_PDB::pdb" );
/*fprintf( pdb, PDB_CRYST1_FORMAT_O,
"CRYST1",
system->big_box.box_norms[0], system->big_box.box_norms[1],
@@ -601,7 +593,7 @@ char Write_PDB(reax_system* system, reax_list** bonds, simulation_data *data,
if ( me == MASTER_NODE)
{
fprintf( pdb, "%s", buffer );
- fclose( pdb );
+ sfclose( pdb, "Write_PDB::pdb" );
}
/* Writing connect information */
diff --git a/PuReMD/src/io_tools.c b/PuReMD/src/io_tools.c
index ef7a76c73d24496aad99cef2d7b60693fa9c98b8..b22ef6bf18e03c4e24f0e6905c6e546ead49d838 100644
--- a/PuReMD/src/io_tools.c
+++ b/PuReMD/src/io_tools.c
@@ -64,66 +64,45 @@ int Init_Output_Files( reax_system *system, control_params *control,
{
/* init out file */
sprintf( temp, "%s.out", control->sim_name );
- if ( (out_control->out = fopen( temp, "w" )) != NULL )
- {
+ out_control->out = sfopen( temp, "w", "Init_Output_Files" );
#if !defined(DEBUG) && !defined(DEBUG_FOCUS)
- fprintf( out_control->out, "%-6s%14s%14s%14s%11s%13s%13s\n",
- "step", "total energy", "potential", "kinetic",
- "T(K)", "V(A^3)", "P(Gpa)" );
+ fprintf( out_control->out, "%-6s%14s%14s%14s%11s%13s%13s\n",
+ "step", "total energy", "potential", "kinetic",
+ "T(K)", "V(A^3)", "P(Gpa)" );
#else
- fprintf( out_control->out, "%-6s%24s%24s%24s%13s%16s%13s\n",
- "step", "total energy", "potential", "kinetic",
- "T(K)", "V(A^3)", "P(GPa)" );
+ fprintf( out_control->out, "%-6s%24s%24s%24s%13s%16s%13s\n",
+ "step", "total energy", "potential", "kinetic",
+ "T(K)", "V(A^3)", "P(GPa)" );
#endif
- fflush( out_control->out );
- }
- else
- {
- strcpy( msg, "init_out_controls: .out file could not be opened\n" );
- return FAILURE;
- }
+ fflush( out_control->out );
/* init potentials file */
sprintf( temp, "%s.pot", control->sim_name );
- if ( (out_control->pot = fopen( temp, "w" )) != NULL )
- {
+ out_control->pot = sfopen( temp, "w", "Init_Output_Files" );
#if !defined(DEBUG) && !defined(DEBUG_FOCUS)
- fprintf( out_control->pot,
- "%-6s%14s%14s%14s%14s%14s%14s%14s%14s%14s%14s%14s\n",
- "step", "ebond", "eatom", "elp",
- "eang", "ecoa", "ehb", "etor", "econj",
- "evdw", "ecoul", "epol" );
+ fprintf( out_control->pot,
+ "%-6s%14s%14s%14s%14s%14s%14s%14s%14s%14s%14s%14s\n",
+ "step", "ebond", "eatom", "elp",
+ "eang", "ecoa", "ehb", "etor", "econj",
+ "evdw", "ecoul", "epol" );
#else
- fprintf( out_control->pot,
- "%-6s%24s%24s%24s%24s%24s%24s%24s%24s%24s%24s%24s\n",
- "step", "ebond", "eatom", "elp",
- "eang", "ecoa", "ehb", "etor", "econj",
- "evdw", "ecoul", "epol" );
+ fprintf( out_control->pot,
+ "%-6s%24s%24s%24s%24s%24s%24s%24s%24s%24s%24s%24s\n",
+ "step", "ebond", "eatom", "elp",
+ "eang", "ecoa", "ehb", "etor", "econj",
+ "evdw", "ecoul", "epol" );
#endif
- fflush( out_control->pot );
- }
- else
- {
- strcpy( msg, "init_out_controls: .pot file could not be opened\n" );
- return FAILURE;
- }
+ fflush( out_control->pot );
/* init log file */
#if defined(LOG_PERFORMANCE)
sprintf( temp, "%s.log", control->sim_name );
- if ( (out_control->log = fopen( temp, "w" )) != NULL )
- {
- fprintf( out_control->log, "%-6s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
- "step", "total", "comm", "neighbors", "init", "bonded", "nonbonded",
- "CM", "CM Sort", "S iters", "Pre Comp", "Pre App", "S comm", "S allr",
- "S spmv", "S vec ops", "S orthog", "S tsolve" );
- fflush( out_control->log );
- }
- else
- {
- strcpy( msg, "init_out_controls: .log file could not be opened\n" );
- return FAILURE;
- }
+ out_control->log = sfopen( temp, "w", "Init_Output_Files" );
+ fprintf( out_control->log, "%-6s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
+ "step", "total", "comm", "neighbors", "init", "bonded", "nonbonded",
+ "CM", "CM Sort", "S iters", "Pre Comp", "Pre App", "S comm", "S allr",
+ "S spmv", "S vec ops", "S orthog", "S tsolve" );
+ fflush( out_control->log );
#endif
}
@@ -133,52 +112,31 @@ int Init_Output_Files( reax_system *system, control_params *control,
control->ensemble == sNPT )
{
sprintf( temp, "%s.prs", control->sim_name );
- if ( (out_control->prs = fopen( temp, "w" )) != NULL )
- {
- fprintf(out_control->prs, "%8s%13s%13s%13s%13s%13s%13s%13s\n",
- "step", "Pint/norm[x]", "Pint/norm[y]", "Pint/norm[z]",
- "Pext/Ptot[x]", "Pext/Ptot[y]", "Pext/Ptot[z]", "Pkin/V" );
- fflush( out_control->prs );
- }
- else
- {
- strcpy(msg, "init_out_controls: .prs file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->prs = sfopen( temp, "w", "Init_Output_Files" );
+ fprintf(out_control->prs, "%8s%13s%13s%13s%13s%13s%13s%13s\n",
+ "step", "Pint/norm[x]", "Pint/norm[y]", "Pint/norm[z]",
+ "Pext/Ptot[x]", "Pext/Ptot[y]", "Pext/Ptot[z]", "Pkin/V" );
+ fflush( out_control->prs );
}
/* init electric dipole moment analysis file */
if ( control->dipole_anal )
{
sprintf( temp, "%s.dpl", control->sim_name );
- if ( (out_control->dpl = fopen( temp, "w" )) != NULL )
- {
- fprintf( out_control->dpl, "%6s%20s%30s",
- "step", "molecule count", "avg dipole moment norm" );
- fflush( out_control->dpl );
- }
- else
- {
- strcpy(msg, "init_out_controls: .dpl file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->dpl = sfopen( temp, "w", "Init_Output_Files" );
+ fprintf( out_control->dpl, "%6s%20s%30s",
+ "step", "molecule count", "avg dipole moment norm" );
+ fflush( out_control->dpl );
}
/* init diffusion coef analysis file */
if ( control->diffusion_coef )
{
sprintf( temp, "%s.drft", control->sim_name );
- if ( (out_control->drft = fopen( temp, "w" )) != NULL )
- {
- fprintf( out_control->drft, "%7s%20s%20s\n",
- "step", "type count", "avg disp^2" );
- fflush( out_control->drft );
- }
- else
- {
- strcpy(msg, "init_out_controls: .drft file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->drft = sfopen( temp, "w", "Init_Output_Files" );
+ fprintf( out_control->drft, "%7s%20s%20s\n",
+ "step", "type count", "avg disp^2" );
+ fflush( out_control->drft );
}
}
@@ -190,9 +148,7 @@ int Init_Output_Files( reax_system *system, control_params *control,
/*if( control->molecular_analysis ) {
if( system->my_rank == MASTER_NODE ) {
sprintf( temp, "%s.mol", control->sim_name );
- if( (out_control->mol = fopen( temp, "w" )) == NULL ) {
- strcpy(msg,"init_out_controls: .mol file could not be opened\n");
- return FAILURE;
+ out_control->mol = sfopen( temp, "w", "Init_Output_Files" );
}
}
@@ -203,233 +159,121 @@ int Init_Output_Files( reax_system *system, control_params *control,
#ifdef TEST_ENERGY
/* open bond energy file */
sprintf( temp, "%s.ebond.%d", control->sim_name, system->my_rank );
- if ( (out_control->ebond = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .ebond file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->ebond = sfopen( temp, "w", "Init_Output_Files" );
/* open lone-pair energy file */
sprintf( temp, "%s.elp.%d", control->sim_name, system->my_rank );
- if ( (out_control->elp = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .elp file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->elp = sfopen( temp, "w", "Init_Output_Files" );
/* open overcoordination energy file */
sprintf( temp, "%s.eov.%d", control->sim_name, system->my_rank );
- if ( (out_control->eov = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .eov file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->eov = sfopen( temp, "w", "Init_Output_Files" );
/* open undercoordination energy file */
sprintf( temp, "%s.eun.%d", control->sim_name, system->my_rank );
- if ( (out_control->eun = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .eun file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->eun = sfopen( temp, "w", "Init_Output_Files" );
/* open angle energy file */
sprintf( temp, "%s.eval.%d", control->sim_name, system->my_rank );
- if ( (out_control->eval = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .eval file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->eval = sfopen( temp, "w", "Init_Output_Files" );
/* open coalition energy file */
sprintf( temp, "%s.ecoa.%d", control->sim_name, system->my_rank );
- if ( (out_control->ecoa = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .ecoa file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->ecoa = sfopen( temp, "w", "Init_Output_Files" );
/* open penalty energy file */
sprintf( temp, "%s.epen.%d", control->sim_name, system->my_rank );
- if ( (out_control->epen = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .epen file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->epen = sfopen( temp, "w", "Init_Output_Files" );
/* open torsion energy file */
sprintf( temp, "%s.etor.%d", control->sim_name, system->my_rank );
- if ( (out_control->etor = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .etor file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->etor = sfopen( temp, "w", "Init_Output_Files" );
/* open conjugation energy file */
sprintf( temp, "%s.econ.%d", control->sim_name, system->my_rank );
- if ( (out_control->econ = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .econ file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->econ = sfopen( temp, "w", "Init_Output_Files" );
/* open hydrogen bond energy file */
sprintf( temp, "%s.ehb.%d", control->sim_name, system->my_rank );
- if ( (out_control->ehb = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .ehb file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->ehb = sfopen( temp, "w", "Init_Output_Files" );
/* open vdWaals energy file */
sprintf( temp, "%s.evdw.%d", control->sim_name, system->my_rank );
- if ( (out_control->evdw = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .evdw file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->evdw = sfopen( temp, "w", "Init_Output_Files" );
/* open coulomb energy file */
sprintf( temp, "%s.ecou.%d", control->sim_name, system->my_rank );
- if ( (out_control->ecou = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .ecou file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->ecou = sfopen( temp, "w", "Init_Output_Files" );
#endif
#ifdef TEST_FORCES
/* open bond orders file */
sprintf( temp, "%s.fbo.%d", control->sim_name, system->my_rank );
- if ( (out_control->fbo = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fbo file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fbo = sfopen( temp, "w", "Init_Output_Files" );
/* open bond orders derivatives file */
sprintf( temp, "%s.fdbo.%d", control->sim_name, system->my_rank );
- if ( (out_control->fdbo = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fdbo file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fdbo = sfopen( temp, "w", "Init_Output_Files" );
/* produce a single force file - to be written by p0 */
if ( system->my_rank == MASTER_NODE )
{
/* open bond forces file */
sprintf( temp, "%s.fbond", control->sim_name );
- if ( (out_control->fbond = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fbond file couldn't be opened\n");
- return FAILURE;
- }
+ Output_Files" )) == NULL )
/* open lone-pair forces file */
sprintf( temp, "%s.flp", control->sim_name );
- if ( (out_control->flp = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .flp file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->flp = sfopen( temp, "w", "Init_Output_Files" );
/* open overcoordination forces file */
sprintf( temp, "%s.fov", control->sim_name );
- if ( (out_control->fov = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fov file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fov = sfopen( temp, "w", "Init_Output_Files" );
/* open undercoordination forces file */
sprintf( temp, "%s.fun", control->sim_name );
- if ( (out_control->fun = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fun file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fun = sfopen( temp, "w", "Init_Output_Files" );
/* open angle forces file */
sprintf( temp, "%s.fang", control->sim_name );
- if ( (out_control->fang = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fang file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fang = sfopen( temp, "w", "Init_Output_Files" );
/* open coalition forces file */
sprintf( temp, "%s.fcoa", control->sim_name );
- if ( (out_control->fcoa = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fcoa file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fcoa = sfopen( temp, "w", "Init_Output_Files" );
/* open penalty forces file */
sprintf( temp, "%s.fpen", control->sim_name );
- if ( (out_control->fpen = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fpen file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fpen = sfopen( temp, "w", "Init_Output_Files" );
/* open torsion forces file */
sprintf( temp, "%s.ftor", control->sim_name );
- if ( (out_control->ftor = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .ftor file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->ftor = sfopen( temp, "w", "Init_Output_Files" );
/* open conjugation forces file */
sprintf( temp, "%s.fcon", control->sim_name );
- if ( (out_control->fcon = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fcon file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fcon = sfopen( temp, "w", "Init_Output_Files" );
/* open hydrogen bond forces file */
sprintf( temp, "%s.fhb", control->sim_name );
- if ( (out_control->fhb = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fhb file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fhb = sfopen( temp, "w", "Init_Output_Files" );
/* open vdw forces file */
sprintf( temp, "%s.fvdw", control->sim_name );
- if ( (out_control->fvdw = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fvdw file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fvdw = sfopen( temp, "w", "Init_Output_Files" );
/* open nonbonded forces file */
sprintf( temp, "%s.fele", control->sim_name );
- if ( (out_control->fele = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fele file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fele = sfopen( temp, "w", "Init_Output_Files" );
/* open total force file */
sprintf( temp, "%s.ftot", control->sim_name );
- if ( (out_control->ftot = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .ftot file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->ftot = sfopen( temp, "w", "Init_Output_Files" );
/* open force comprison file */
sprintf( temp, "%s.fcomp", control->sim_name );
- if ( (out_control->fcomp = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .fcomp file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->fcomp = sfopen( temp, "w", "Init_Output_Files" );
}
#endif
@@ -437,27 +281,15 @@ int Init_Output_Files( reax_system *system, control_params *control,
#if defined(TEST_FORCES) || defined(TEST_ENERGY)
/* open far neighbor list file */
sprintf( temp, "%s.far_nbrs_list.%d", control->sim_name, system->my_rank );
- if ( (out_control->flist = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .far_nbrs_list file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->flist = sfopen( temp, "w", "Init_Output_Files" );
/* open bond list file */
sprintf( temp, "%s.bond_list.%d", control->sim_name, system->my_rank );
- if ( (out_control->blist = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .bond_list file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->blist = sfopen( temp, "w", "Init_Output_Files" );
/* open near neighbor list file */
sprintf( temp, "%s.near_nbrs_list.%d", control->sim_name, system->my_rank );
- if ( (out_control->nlist = fopen( temp, "w" )) == NULL )
- {
- strcpy(msg, "Init_Out_Files: .near_nbrs_list file couldn't be opened\n");
- return FAILURE;
- }
+ out_control->nlist = sfopen( temp, "w", "Init_Output_Files" );
#endif
#endif
@@ -476,65 +308,68 @@ int Close_Output_Files( reax_system *system, control_params *control,
{
if ( out_control->energy_update_freq > 0 )
{
- fclose( out_control->out );
- fclose( out_control->pot );
+ sfclose( out_control->out, "Close_Output_Files" );
+ sfclose( out_control->pot, "Close_Output_Files" );
#if defined(LOG_PERFORMANCE)
- fclose( out_control->log );
+ sfclose( out_control->log, "Close_Output_Files" );
#endif
}
if ( control->ensemble == NPT || control->ensemble == iNPT ||
control->ensemble == sNPT )
- fclose( out_control->prs );
+ sfclose( out_control->prs, "Close_Output_Files" );
- if ( control->dipole_anal ) fclose( out_control->dpl );
- if ( control->diffusion_coef ) fclose( out_control->drft );
- if ( control->molecular_analysis ) fclose( out_control->mol );
+ if ( control->dipole_anal )
+ sfclose( out_control->dpl, "Close_Output_Files" );
+ if ( control->diffusion_coef )
+ sfclose( out_control->drft, "Close_Output_Files" );
+ if ( control->molecular_analysis )
+ sfclose( out_control->mol, "Close_Output_Files" );
}
#ifdef TEST_ENERGY
- fclose( out_control->ebond );
- fclose( out_control->elp );
- fclose( out_control->eov );
- fclose( out_control->eun );
- fclose( out_control->eval );
- fclose( out_control->epen );
- fclose( out_control->ecoa );
- fclose( out_control->ehb );
- fclose( out_control->etor );
- fclose( out_control->econ );
- fclose( out_control->evdw );
- fclose( out_control->ecou );
+ sfclose( out_control->ebond, "Close_Output_Files" );
+ sfclose( out_control->elp, "Close_Output_Files" );
+ sfclose( out_control->eov, "Close_Output_Files" );
+ sfclose( out_control->eun, "Close_Output_Files" );
+ sfclose( out_control->eval, "Close_Output_Files" );
+ sfclose( out_control->epen, "Close_Output_Files" );
+ sfclose( out_control->ecoa, "Close_Output_Files" );
+ sfclose( out_control->ehb, "Close_Output_Files" );
+ sfclose( out_control->etor, "Close_Output_Files" );
+ sfclose( out_control->econ, "Close_Output_Files" );
+ sfclose( out_control->evdw, "Close_Output_Files" );
+ sfclose( out_control->ecou, "Close_Output_Files" );
#endif
#ifdef TEST_FORCES
- fclose( out_control->fbo );
- fclose( out_control->fdbo );
+ sfclose( out_control->fbo, "Close_Output_Files" );
+ sfclose( out_control->fdbo, "Close_Output_Files" );
if ( system->my_rank == MASTER_NODE )
{
- fclose( out_control->fbond );
- fclose( out_control->flp );
- fclose( out_control->fov );
- fclose( out_control->fun );
- fclose( out_control->fang );
- fclose( out_control->fcoa );
- fclose( out_control->fpen );
- fclose( out_control->ftor );
- fclose( out_control->fcon );
- fclose( out_control->fhb );
- fclose( out_control->fvdw );
- fclose( out_control->fele );
- fclose( out_control->ftot );
- fclose( out_control->fcomp );
+ sfclose( out_control->fbond, "Close_Output_Files" );
+ sfclose( out_control->flp, "Close_Output_Files" );
+ sfclose( out_control->fov, "Close_Output_Files" );
+ sfclose( out_control->fun, "Close_Output_Files" );
+ sfclose( out_control->fang, "Close_Output_Files" );
+ sfclose( out_control->fcoa, "Close_Output_Files" );
+ sfclose( out_control->fpen, "Close_Output_Files" );
+ sfclose( out_control->ftor, "Close_Output_Files" );
+ sfclose( out_control->fcon, "Close_Output_Files" );
+ sfclose( out_control->fhb, "Close_Output_Files" );
+ sfclose( out_control->fvdw, "Close_Output_Files" );
+ sfclose( out_control->fele, "Close_Output_Files" );
+ sfclose( out_control->ftot, "Close_Output_Files" );
+ sfclose( out_control->fcomp, "Close_Output_Files" );
}
#endif
#if defined(PURE_REAX)
#if defined(TEST_FORCES) || defined(TEST_ENERGY)
- fclose( out_control->flist );
- fclose( out_control->blist );
- fclose( out_control->nlist );
+ sfclose( out_control->flist, "Close_Output_Files" );
+ sfclose( out_control->blist, "Close_Output_Files" );
+ sfclose( out_control->nlist, "Close_Output_Files" );
#endif
#endif
@@ -668,7 +503,7 @@ void Print_GCell_Exchange_Bounds( int my_rank, neighbor_proc *my_nbrs )
char exch[3][10] = { "NONE", "NEAR_EXCH", "FULL_EXCH" };
sprintf( fname, "gcell_exchange_bounds%d", my_rank );
- f = fopen( fname, "w" );
+ f = sfopen( fname, "w", "Print_GCell_Exchange_Bounds" );
/* loop over neighbor processes */
for ( r[0] = -1; r[0] <= 1; ++r[0])
@@ -696,7 +531,7 @@ void Print_GCell_Exchange_Bounds( int my_rank, neighbor_proc *my_nbrs )
nbr_pr->end_recv[2] );
}
- fclose(f);
+ sfclose( f, "Print_GCell_Exchange_Bounds" );
}
@@ -715,7 +550,7 @@ void Print_Native_GCells( reax_system *system )
};
sprintf( fname, "native_gcells.%d", system->my_rank );
- f = fopen( fname, "w" );
+ f = sfopen( fname, "w", "Print_Native_GCells" );
g = &(system->my_grid);
for ( i = g->native_str[0]; i < g->native_end[0]; i++ )
@@ -735,7 +570,7 @@ void Print_Native_GCells( reax_system *system )
fprintf( f, "\n" );
}
- fclose(f);
+ sfclose( f, "Print_Native_GCells" );
}
@@ -754,7 +589,7 @@ void Print_All_GCells( reax_system *system )
};
sprintf( fname, "all_gcells.%d", system->my_rank );
- f = fopen( fname, "w" );
+ f = sfopen( fname, "w", "Print_All_GCells" );
g = &(system->my_grid);
for ( i = 0; i < g->ncells[0]; i++ )
@@ -774,7 +609,7 @@ void Print_All_GCells( reax_system *system )
fprintf( f, "\n" );
}
- fclose(f);
+ sfclose( f, "Print_All_GCells" );
}
@@ -786,11 +621,7 @@ void Print_My_Atoms( reax_system *system, control_params *control, int step )
FILE *fh;
sprintf( fname, "%s.my_atoms.%d.%d", control->sim_name, step, system->my_rank );
- if ( (fh = fopen( fname, "w" )) == NULL )
- {
- fprintf( stderr, "error in opening my_atoms file" );
- MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
- }
+ fh = sfopen( fname, "w", "Print_My_Atoms" );
// fprintf( stderr, "p%d had %d atoms\n",
// system->my_rank, system->n );
@@ -803,7 +634,7 @@ void Print_My_Atoms( reax_system *system, control_params *control, int step )
system->my_atoms[i].x[1],
system->my_atoms[i].x[2] );
- fclose( fh );
+ sfclose( fh, "Print_My_Atoms" );
}
@@ -814,11 +645,7 @@ void Print_My_Ext_Atoms( reax_system *system )
FILE *fh;
sprintf( fname, "my_ext_atoms.%d", system->my_rank );
- if ( (fh = fopen( fname, "w" )) == NULL )
- {
- fprintf( stderr, "error in opening my_ext_atoms file" );
- MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
- }
+ fh = sfopen( fname, "w", "Print_My_Ext_Atoms" );
// fprintf( stderr, "p%d had %d atoms\n",
// system->my_rank, system->n );
@@ -831,7 +658,7 @@ void Print_My_Ext_Atoms( reax_system *system )
system->my_atoms[i].x[1],
system->my_atoms[i].x[2] );
- fclose( fh );
+ sfclose( fh, "Print_My_Ext_Atoms" );
}
@@ -844,7 +671,7 @@ void Print_Far_Neighbors( reax_system *system, reax_list **lists,
reax_list *far_nbrs;
sprintf( fname, "%s.far_nbrs.%d", control->sim_name, system->my_rank );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w", "Print_Far_Neighbors" );
far_nbrs = lists[FAR_NBRS];
natoms = system->N;
@@ -871,7 +698,7 @@ void Print_Far_Neighbors( reax_system *system, reax_list **lists,
}
}
- fclose( fout );
+ sfclose( fout, "Print_Far_Neighbors" );
}
@@ -891,7 +718,7 @@ void Print_Sparse_Matrix( reax_system *system, sparse_matrix *A )
void Print_Sparse_Matrix2( reax_system *system, sparse_matrix *A, char *fname )
{
int i, j;
- FILE *f = fopen( fname, "w" );
+ FILE *f = sfopen( fname, "w", "Print_Sparse_Matrix2" );
for ( i = 0; i < A->n; ++i )
for ( j = A->start[i]; j < A->end[i]; ++j )
@@ -900,7 +727,7 @@ void Print_Sparse_Matrix2( reax_system *system, sparse_matrix *A, char *fname )
system->my_atoms[A->entries[j].j].orig_id,
A->entries[j].val );
- fclose(f);
+ sfclose( f, "Print_Sparse_Matrix2" );
}
@@ -908,7 +735,7 @@ void Print_Symmetric_Sparse(reax_system *system, sparse_matrix *A, char *fname)
{
int i, j;
reax_atom *ai, *aj;
- FILE *f = fopen( fname, "w" );
+ FILE *f = sfopen( fname, "w", "Print_Symmetric_Sparse" );
for ( i = 0; i < A->n; ++i )
{
@@ -924,7 +751,7 @@ void Print_Symmetric_Sparse(reax_system *system, sparse_matrix *A, char *fname)
}
}
- fclose(f);
+ sfclose( f, "Print_Symmetric_Sparse" );
}
@@ -939,7 +766,7 @@ void Print_Linear_System( reax_system *system, control_params *control,
// print rhs and init guesses for QEq
sprintf( fname, "%s.p%dstate%d", control->sim_name, system->my_rank, step );
- out = fopen( fname, "w" );
+ out = sfopen( fname, "w", "Print_Linear_System" );
for ( i = 0; i < system->n; i++ )
{
ai = &(system->my_atoms[i]);
@@ -948,7 +775,7 @@ void Print_Linear_System( reax_system *system, control_params *control,
workspace->s[i], workspace->b_s[i],
workspace->t[i], workspace->b_t[i] );
}
- fclose( out );
+ sfclose( out, "Print_Linear_System" );
// print QEq coef matrix
sprintf( fname, "%s.p%dH%d", control->sim_name, system->my_rank, step );
@@ -956,7 +783,7 @@ void Print_Linear_System( reax_system *system, control_params *control,
// print the incomplete H matrix
/*sprintf( fname, "%s.p%dHinc%d", control->sim_name, system->my_rank, step );
- out = fopen( fname, "w" );
+ out = sfopen( fname, "w", "Print_Linear_System" );
H = workspace->H;
for( i = 0; i < H->n; ++i ) {
ai = &(system->my_atoms[i]);
@@ -970,7 +797,7 @@ void Print_Linear_System( reax_system *system, control_params *control,
aj->orig_id, ai->orig_id, H->entries[j].val );
}
}
- fclose( out );*/
+ sfclose( out, "Print_Linear_System" );*/
// print the L from incomplete cholesky decomposition
/*sprintf( fname, "%s.p%dL%d", control->sim_name, system->my_rank, step );
@@ -985,13 +812,13 @@ void Print_LinSys_Soln( reax_system *system, real *x, real *b_prm, real *b )
FILE *fout;
sprintf( fname, "qeq.%d.out", system->my_rank );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w", "Print_LinSys_Soln" );
for ( i = 0; i < system->n; ++i )
fprintf( fout, "%6d%10.4f%10.4f%10.4f\n",
system->my_atoms[i].orig_id, x[i], b_prm[i], b[i] );
- fclose( fout );
+ sfclose( fout, "Print_LinSys_Soln" );
}
@@ -1002,7 +829,7 @@ void Print_Charges( reax_system *system )
FILE *fout;
sprintf( fname, "q.%d.out", system->my_rank );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w", "Print_Charges" );
for ( i = 0; i < system->n; ++i )
fprintf( fout, "%6d %10.7f %10.7f %10.7f\n",
@@ -1011,7 +838,7 @@ void Print_Charges( reax_system *system )
system->my_atoms[i].t[0],
system->my_atoms[i].q );
- fclose( fout );
+ sfclose( fout, "Print_Charges" );
}
@@ -1025,7 +852,7 @@ void Print_HBonds( reax_system *system, reax_list **lists,
reax_list *hbonds = lists[HBONDS];
sprintf( fname, "%s.hbonds.%d.%d", control->sim_name, step, system->my_rank );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w", "Print_HBonds" );
for ( i = 0; i < system->numH; ++i )
{
@@ -1040,7 +867,7 @@ void Print_HBonds( reax_system *system, reax_list **lists,
}
}
- fclose( fout );
+ sfclose( fout, "Print_HBonds" );
}
@@ -1053,7 +880,7 @@ void Print_HBond_Indices( reax_system *system, reax_list **lists,
reax_list *hbonds = lists[HBONDS];
sprintf( fname, "%s.hbonds_indices.%d.%d", control->sim_name, step, system->my_rank );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w", "Print_HBond_Indices" );
for ( i = 0; i < system->N; ++i )
{
@@ -1061,7 +888,7 @@ void Print_HBond_Indices( reax_system *system, reax_list **lists,
i, Start_Index(i, hbonds), End_Index(i, hbonds) );
}
- fclose( fout );
+ sfclose( fout, "Print_HBond_Indices" );
}
@@ -1076,7 +903,7 @@ void Print_Bonds( reax_system *system, reax_list **lists,
reax_list *bonds = lists[BONDS];
sprintf( fname, "%s.bonds.%d.%d", control->sim_name, step, system->my_rank );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w", "Print_Bonds" );
for ( i = 0; i < system->N; ++i )
{
@@ -1093,7 +920,7 @@ void Print_Bonds( reax_system *system, reax_list **lists,
}
}
- fclose( fout );
+ sfclose( fout, "Print_Bonds" );
}
@@ -1105,7 +932,7 @@ int fn_qsort_intcmp( const void *a, const void *b )
void Print_Bond_List2( reax_system *system, reax_list *bonds, char *fname )
{
int i, j, id_i, id_j, nbr, pj;
- FILE *f = fopen( fname, "w" );
+ FILE *f = sfopen( fname, "w", "Print_Bond_List2" );
int temp[500];
int num = 0;
@@ -1156,7 +983,7 @@ void Print_Far_Neighbors_List_Adj_Format( reax_system *system,
FILE *fout;
sprintf( fname, "%s.far.%d.%d", control->sim_name, step, system->my_rank );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w", "Print_Far_Neighbors_Adj_Format" );
num_intrs = 0;
intrs = NULL;
@@ -1203,7 +1030,7 @@ void Print_Far_Neighbors_List_Adj_Format( reax_system *system,
free( intrs );
}
- fclose( fout );
+ sfclose( fout, "Print_Far_Neighbors_List_Adj_Format" );
}
void Output_Results( reax_system *system, control_params *control,
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index a404bdb113d3ee539ba15c0a1363f7476eac9a92..98bc9b2923baf2259db4dc3062801b2202e3bcd1 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -1268,7 +1268,6 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
t_start = Get_Time( );
tmp = Parallel_Dot(workspace->d, workspace->q, system->n, mpi_data->world);
- //TODO: all_Reduce time
t_allreduce += Get_Timing_Info ( t_start );
t_start = Get_Time( );
@@ -1302,7 +1301,6 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
t_start = Get_Time( );
sig_old = sig_new;
sig_new = Parallel_Dot(workspace->r, workspace->p, system->n, mpi_data->world);
- //TODO all_reduce time
t_allreduce += Get_Timing_Info( t_start );
t_start = Get_Time( );
diff --git a/PuReMD/src/reax_defs.h b/PuReMD/src/reax_defs.h
index dc2ac1b925cc9f60daca52d0520d5e2e173e19bb..310a46f769ba98f67f6ee4ce22d9fba6b1e30944 100644
--- a/PuReMD/src/reax_defs.h
+++ b/PuReMD/src/reax_defs.h
@@ -100,7 +100,6 @@
#define RESTART (30)
-
/******************* ENUMERATIONS *************************/
enum geo_formats { CUSTOM = 0, PDB = 1, ASCII_RESTART = 2, BINARY_RESTART = 3, GF_N = 4 };
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 25749323d464d581388a486cb466ed20a179fb52..2b08d03e15726c1893cf236dbe22685b1313edea 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -54,6 +54,7 @@
#define ZERO (0.000000000000000e+00)
#define REAX_MAX_STR 1024
#define REAX_MAX_NBRS 6
+#define REAX_MAX_NT_NBRS 6
#define REAX_MAX_3BODY_PARAM 5
#define REAX_MAX_4BODY_PARAM 5
#define REAX_MAX_ATOM_TYPES 25
@@ -224,8 +225,10 @@ typedef struct
//MPI_Status recv_stat2[REAX_MAX_NBRS];
mpi_out_data out_buffers[REAX_MAX_NBRS];
+ mpi_out_data out_nt_buffers[REAX_MAX_NT_NBRS];
void *in1_buffer;
void *in2_buffer;
+ void *in_nt_buffer[REAX_MAX_NT_NBRS];
} mpi_datatypes;
@@ -503,6 +506,7 @@ typedef struct
typedef struct
{
int rank;
+ int receive_rank;
int est_send, est_recv;
int atoms_str, atoms_cnt;
ivec rltv, prdc;
@@ -544,9 +548,10 @@ typedef struct
int n, N, bigN, numH;
int local_cap, total_cap, gcell_cap, Hcap;
int est_recv, est_trans, max_recved;
- int wsize, my_rank, num_nbrs;
+ int wsize, my_rank, num_nbrs, num_nt_nbrs;
ivec my_coords;
neighbor_proc my_nbrs[REAX_MAX_NBRS];
+ neighbor_proc my_nt_nbrs[REAX_MAX_NT_NBRS];
int *global_offset;
simulation_box big_box, my_box, my_ext_box;
grid my_grid;
diff --git a/PuReMD/src/restart.c b/PuReMD/src/restart.c
index b687d82abc0a2618b0cc7d496f39b8630d51d450..bdd2476f74fdbe797c009643b021ef69f1fd78a6 100644
--- a/PuReMD/src/restart.c
+++ b/PuReMD/src/restart.c
@@ -49,11 +49,7 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
{
/* master handles the restart file */
sprintf( fname, "%s.res%d", control->sim_name, data->step );
- if ( (fres = fopen( fname, "wb" )) == NULL )
- {
- fprintf( stderr, "ERROR: can't open the restart file! terminating...\n" );
- MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
- }
+ fres = sfopen( fname, "wb", "Write_Binary_Restart" );
/* master can write the header by itself */
res_header.step = data->step;
@@ -108,7 +104,7 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
if ( me == MASTER_NODE )
{
fwrite( buffer, system->bigN, sizeof(restart_atom), fres );
- fclose( fres );
+ sfclose( fres, "Write_Binary_Restart" );
}
sfree(buffer, "buffer");
@@ -139,11 +135,7 @@ void Write_Restart( reax_system *system, control_params *control,
if ( me == MASTER_NODE )
{
sprintf( fname, "%s.res%d", control->sim_name, data->step );
- if ( (fres = fopen( fname, "w" )) == NULL )
- {
- fprintf( stderr, "ERROR: can't open the restart file! terminating...\n" );
- MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
- }
+ fres = sfopen( fname, "w", "Write_Restart" );
/* write the header - only master writes it */
fprintf( fres, RESTART_HEADER,
@@ -204,7 +196,7 @@ void Write_Restart( reax_system *system, control_params *control,
if ( me == MASTER_NODE )
{
fprintf( fres, "%s", buffer );
- fclose( fres );
+ sfclose( fres, "Write_Restart" );
}
sfree(buffer, "buffer");
sfree(line, "line");
@@ -250,11 +242,7 @@ void Read_Binary_Restart( char *res_file, reax_system *system,
comm = MPI_COMM_WORLD;
- if ( (fres = fopen(res_file, "rb")) == NULL )
- {
- fprintf( stderr, "ERROR: cannot open the restart file! terminating...\n" );
- MPI_Abort( comm, FILE_NOT_FOUND );
- }
+ fres = sfopen( res_file, "rb", "Read_Binary_Restart" );
/* first read the header lines */
fread(&res_header, sizeof(restart_header), 1, fres);
@@ -313,7 +301,7 @@ void Read_Binary_Restart( char *res_file, reax_system *system,
}
}
- fclose( fres );
+ sfclose( fres, "Read_Binary_Restart" );
data->step = data->prev_steps;
// nsteps is updated based on the number of steps in the previous run
@@ -368,11 +356,7 @@ void Read_Restart( char *res_file, reax_system *system,
comm = MPI_COMM_WORLD;
- if ( (fres = fopen(res_file, "r")) == NULL )
- {
- fprintf( stderr, "ERROR: cannot open the restart file! terminating...\n" );
- MPI_Abort( comm, FILE_NOT_FOUND );
- }
+ fres = sfopen( res_file, "r", "Read_Binary_Restart" );
s = (char*) malloc(sizeof(char) * MAX_LINE);
tmp = (char**) malloc(sizeof(char*)*MAX_TOKENS);
@@ -464,7 +448,7 @@ void Read_Restart( char *res_file, reax_system *system,
top++;
}
}
- fclose( fres );
+ sfclose( fres, "Read_Restart" );
/* free memory allocations at the top */
for ( i = 0; i < MAX_TOKENS; i++ )
sfree( tmp[i], "tmp[i]" );
diff --git a/PuReMD/src/torsion_angles.c b/PuReMD/src/torsion_angles.c
index f915fdb8325de91f5c5ea72f47e36ae0ab306992..8d9949dadee37483c189e16ac79aad2c27e0b03e 100644
--- a/PuReMD/src/torsion_angles.c
+++ b/PuReMD/src/torsion_angles.c
@@ -197,7 +197,7 @@ void Torsion_Angles( reax_system *system, control_params *control,
// FILE *ftor;
// sprintf( fname, "tor%d.out", system->my_rank );
- // ftor = fopen( fname, "w" );
+ // ftor = sfopen( fname, "w", "Torsion_Angles" );
natoms = system->n;
diff --git a/PuReMD/src/traj.c b/PuReMD/src/traj.c
index d03a8289f22639477db5c9c7b3af919d37adba9e..e189c4fc2c479c6bb0a5b7a7180e267c820d29d2 100644
--- a/PuReMD/src/traj.c
+++ b/PuReMD/src/traj.c
@@ -527,7 +527,7 @@ int Init_Traj( reax_system *system, control_params *control,
else if ( out_control->traj_method == REG_TRAJ)
{
if ( system->my_rank == MASTER_NODE )
- out_control->strj = fopen( fname, "w" );
+ out_control->strj = sfopen( fname, "w", "Init_Traj" );
}
else
{
@@ -538,7 +538,7 @@ int Init_Traj( reax_system *system, control_params *control,
if ( out_control->traj_method == REG_TRAJ)
{
if ( system->my_rank == MASTER_NODE )
- out_control->strj = fopen( fname, "w" );
+ out_control->strj = sfopen( fname, "w", "Init_Traj" );
}
else
{
@@ -1116,10 +1116,10 @@ int End_Traj( int my_rank, output_controls *out_control )
if ( out_control->traj_method == MPI_TRAJ )
MPI_File_close( &(out_control->trj) );
else if ( my_rank == MASTER_NODE )
- fclose( out_control->strj );
+ sfclose( out_control->strj, "End_Traj" );
#elif defined(LAMMPS_REAX)
if ( my_rank == MASTER_NODE )
- fclose( out_control->strj );
+ sfclose( out_control->strj, "End_Traj" );
#endif
sfree( out_control->buffer, "out_control->buffer" );