diff --git a/PG-PuReMD/Makefile.am b/PG-PuReMD/Makefile.am
index 5681e43a93fcb419f27e62ae25e3d72e5ee8740c..c5c84052219e835a12a5529ceff6223a524dd1a3 100644
--- a/PG-PuReMD/Makefile.am
+++ b/PG-PuReMD/Makefile.am
@@ -16,15 +16,15 @@ NVCCFLAGS += --compiler-options "$(DEFS) $(NVCC_OPT_CODE_DEFS) -O3 -funroll-loop
endif
-bin_PROGRAMS = bin/pg-puremd
-bin_pg_puremd_SOURCES = src/allocate.c src/basic_comm.c src/ffield.c src/grid.c src/list.c \
+lib_LTLIBRARIES = lib/libpuremd.la
+lib_libpuremd_la_SOURCES = src/allocate.c src/basic_comm.c src/ffield.c src/grid.c src/list.c \
src/lookup.c src/io_tools.c src/reset_tools.c src/restart.c src/random.c \
src/tool_box.c src/traj.c src/analyze.c src/box.c src/system_props.c \
src/control.c src/comm_tools.c src/geo_tools.c src/lin_alg.c src/neighbors.c \
src/charges.c src/bond_orders.c src/multi_body.c src/bonds.c src/valence_angles.c \
src/hydrogen_bonds.c src/torsion_angles.c src/nonbonded.c src/forces.c \
- src/integrate.c src/init_md.c src/parallelreax.c
-include_HEADERS = src/reax_types.h src/index_utils.h \
+ src/integrate.c src/init_md.c src/puremd.c
+lib_libpuremd_la_SOURCES += src/reax_types.h src/index_utils.h \
src/allocate.h src/basic_comm.h src/ffield.h src/grid.h src/list.h \
src/lookup.h src/io_tools.h src/reset_tools.h src/restart.h src/random.h src/vector.h \
src/tool_box.h src/traj.h src/analyze.h src/box.h src/system_props.h \
@@ -32,9 +32,10 @@ include_HEADERS = src/reax_types.h src/index_utils.h \
src/charges.h src/bond_orders.h src/multi_body.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
+include_HEADERS = src/puremd.h
if USE_CUDA
-bin_pg_puremd_SOURCES += src/cuda/cuda_utils.cu src/cuda/cuda_allocate.cu src/cuda/cuda_environment.cu \
+lib_libpuremd_la_SOURCES += src/cuda/cuda_utils.cu src/cuda/cuda_allocate.cu src/cuda/cuda_environment.cu \
src/cuda/cuda_system_props.cu src/cuda/cuda_reduction.cu src/cuda/cuda_box.cu src/cuda/cuda_list.cu \
src/cuda/cuda_copy.cu src/cuda/cuda_reset_tools.cu src/cuda/cuda_random.cu \
src/cuda/cuda_neighbors.cu src/cuda/cuda_bond_orders.cu src/cuda/cuda_bonds.cu \
@@ -43,7 +44,7 @@ bin_pg_puremd_SOURCES += src/cuda/cuda_utils.cu src/cuda/cuda_allocate.cu src/cu
src/cuda/cuda_charges.cu src/cuda/cuda_lin_alg.cu \
src/cuda/cuda_nonbonded.cu src/cuda/cuda_integrate.cu src/cuda/cuda_post_evolve.cu \
src/cuda/cuda_init_md.cu src/cuda/cuda_lookup.cu
-include_HEADERS += src/cuda/cuda_helpers.h src/cuda/cuda_shuffle.h \
+lib_libpuremd_la_SOURCES += src/cuda/cuda_helpers.h src/cuda/cuda_shuffle.h \
src/cuda/cuda_utils.h src/cuda/cuda_allocate.h src/cuda/cuda_environment.h \
src/cuda/cuda_system_props.h src/cuda/cuda_reduction.h src/cuda/cuda_box.h src/cuda/cuda_list.h \
src/cuda/cuda_copy.h src/cuda/cuda_reset_tools.h src/cuda/cuda_random.h src/cuda/cuda_vector.h \
@@ -55,20 +56,28 @@ include_HEADERS += src/cuda/cuda_helpers.h src/cuda/cuda_shuffle.h \
src/cuda/cuda_init_md.h src/cuda/cuda_lookup.h
if DEBUG
-bin_pg_puremd_SOURCES += src/cuda/cuda_validation.cu
-include_HEADERS += src/cuda/cuda_validation.h
+lib_libpuremd_la_SOURCES += src/cuda/cuda_validation.cu
+lib_libpuremd_la_SOURCES += src/cuda/cuda_validation.h
endif
# dummy source to cause C linking
-nodist_EXTRA_bin_pg_puremd_SOURCES = src/dummy.c
-
+nodist_EXTRA_lib_libpuremd_la_SOURCES = src/dummy.c
endif
+lib_libpuremd_la_CFLAGS = $(AM_CFLAGS) $(MPI_CFLAGS) $(CFLAGS)
+lib_libpuremd_la_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS)
+lib_libpuremd_la_LIBADD = $(AM_LDADD) $(MPI_LIBS) $(LDADD) -lstdc++
+lib_libpuremd_la_LDFLAGS = -version-info 1:0:0
+if USE_CUDA
+lib_libpuremd_la_CFLAGS += $(CUDA_CFLAGS)
+lib_libpuremd_la_LIBADD += $(CUDA_LIBS)
+endif
+bin_PROGRAMS = bin/pg-puremd
+bin_pg_puremd_SOURCES = src/driver.c
bin_pg_puremd_CFLAGS = $(AM_CFLAGS) $(MPI_CFLAGS) $(CFLAGS)
-bin_pg_puremd_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS)
-bin_pg_puremd_LDADD = $(AM_LDADD) $(MPI_LIBS) $(LDADD) -lstdc++
-
+bin_pg_puremd_CPPFLAGS = -I src $(AM_CPPFLAGS) $(CPPFLAGS)
+bin_pg_puremd_LDADD = lib/libpuremd.la $(AM_LDADD) $(MPI_LIBS) $(LDADD) -lstdc++
if USE_CUDA
bin_pg_puremd_CFLAGS += $(CUDA_CFLAGS)
bin_pg_puremd_LDADD += $(CUDA_LIBS)
diff --git a/PG-PuReMD/aclocal.m4 b/PG-PuReMD/aclocal.m4
index e1a7adedccfa9f18035be1d931ad439b6df5dc38..8206dbb643189be5c98f9d5fa31a85ad053315ff 100644
--- a/PG-PuReMD/aclocal.m4
+++ b/PG-PuReMD/aclocal.m4
@@ -56,6 +56,66 @@ m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])
+# Copyright (C) 2011-2017 Free Software Foundation, Inc.
+#
+# This file is free software; the Free Software Foundation
+# gives unlimited permission to copy and/or distribute it,
+# with or without modifications, as long as this notice is preserved.
+
+# AM_PROG_AR([ACT-IF-FAIL])
+# -------------------------
+# Try to determine the archiver interface, and trigger the ar-lib wrapper
+# if it is needed. If the detection of archiver interface fails, run
+# ACT-IF-FAIL (default is to abort configure with a proper error message).
+AC_DEFUN([AM_PROG_AR],
+[AC_BEFORE([$0], [LT_INIT])dnl
+AC_BEFORE([$0], [AC_PROG_LIBTOOL])dnl
+AC_REQUIRE([AM_AUX_DIR_EXPAND])dnl
+AC_REQUIRE_AUX_FILE([ar-lib])dnl
+AC_CHECK_TOOLS([AR], [ar lib "link -lib"], [false])
+: ${AR=ar}
+
+AC_CACHE_CHECK([the archiver ($AR) interface], [am_cv_ar_interface],
+ [AC_LANG_PUSH([C])
+ am_cv_ar_interface=ar
+ AC_COMPILE_IFELSE([AC_LANG_SOURCE([[int some_variable = 0;]])],
+ [am_ar_try='$AR cru libconftest.a conftest.$ac_objext >&AS_MESSAGE_LOG_FD'
+ AC_TRY_EVAL([am_ar_try])
+ if test "$ac_status" -eq 0; then
+ am_cv_ar_interface=ar
+ else
+ am_ar_try='$AR -NOLOGO -OUT:conftest.lib conftest.$ac_objext >&AS_MESSAGE_LOG_FD'
+ AC_TRY_EVAL([am_ar_try])
+ if test "$ac_status" -eq 0; then
+ am_cv_ar_interface=lib
+ else
+ am_cv_ar_interface=unknown
+ fi
+ fi
+ rm -f conftest.lib libconftest.a
+ ])
+ AC_LANG_POP([C])])
+
+case $am_cv_ar_interface in
+ar)
+ ;;
+lib)
+ # Microsoft lib, so override with the ar-lib wrapper script.
+ # FIXME: It is wrong to rewrite AR.
+ # But if we don't then we get into trouble of one sort or another.
+ # A longer-term fix would be to have automake use am__AR in this case,
+ # and then we could set am__AR="$am_aux_dir/ar-lib \$(AR)" or something
+ # similar.
+ AR="$am_aux_dir/ar-lib $AR"
+ ;;
+unknown)
+ m4_default([$1],
+ [AC_MSG_ERROR([could not determine $AR interface])])
+ ;;
+esac
+AC_SUBST([AR])dnl
+])
+
# AM_AUX_DIR_EXPAND -*- Autoconf -*-
# Copyright (C) 2001-2017 Free Software Foundation, Inc.
@@ -1153,3 +1213,8 @@ AC_SUBST([am__untar])
m4_include([../m4/acx_mpi.m4])
m4_include([../m4/ax_compiler_vendor.m4])
m4_include([../m4/ax_cuda.m4])
+m4_include([../m4/libtool.m4])
+m4_include([../m4/ltoptions.m4])
+m4_include([../m4/ltsugar.m4])
+m4_include([../m4/ltversion.m4])
+m4_include([../m4/lt~obsolete.m4])
diff --git a/PG-PuReMD/configure.ac b/PG-PuReMD/configure.ac
index bd79fae4f44e55a3bfef67f00beb8458d4094fd3..ff85bf947cc35b417af4c257b625a6365339aee1 100644
--- a/PG-PuReMD/configure.ac
+++ b/PG-PuReMD/configure.ac
@@ -10,8 +10,9 @@ sav_CFLAGS="$CFLAGS"
AM_INIT_AUTOMAKE([1.15 subdir-objects -Wall -Werror foreign])
# Enable silent build rules by default.
m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])], [AC_SUBST([AM_DEFAULT_VERBOSITY],[1])])
-#LT_PREREQ([2.2])
-#LT_INIT([dlopen])
+AM_PROG_AR
+LT_PREREQ([2.2])
+LT_INIT([dlopen])
AC_CONFIG_MACRO_DIR([../m4])
diff --git a/PG-PuReMD/src/driver.c b/PG-PuReMD/src/driver.c
new file mode 100644
index 0000000000000000000000000000000000000000..10e86b412d588607b69bbd51f4f02cfb1e7feb47
--- /dev/null
+++ b/PG-PuReMD/src/driver.c
@@ -0,0 +1,70 @@
+
+/*----------------------------------------------------------------------
+ SerialReax - Reax Force Field Simulator
+
+ Copyright (2010) Purdue University
+ Hasan Metin Aktulga, haktulga@cs.purdue.edu
+ Joseph Fogarty, jcfogart@mail.usf.edu
+ Sagar Pandit, pandit@usf.edu
+ Ananth Y Grama, ayg@cs.purdue.edu
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU General Public License as
+ published by the Free Software Foundation; either version 2 of
+ the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ See the GNU General Public License for more details:
+ <http://www.gnu.org/licenses/>.
+ ----------------------------------------------------------------------*/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "puremd.h"
+
+#define INVALID_INPUT (-1)
+
+
+static void usage( char * argv[] )
+{
+ fprintf( stderr, "usage: ./%s geometry_file force_field_params_file control_file\n", argv[0] );
+}
+
+
+int main( int argc, char* argv[] )
+{
+ void *handle;
+ int ret;
+
+ MPI_Init( &argc, &argv );
+
+ if ( argc != 4 )
+ {
+ usage( argv );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
+ }
+
+ handle = setup( argv[1], argv[2], argv[3] );
+ ret = PUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ ret = simulate( handle );
+ }
+
+ MPI_Finalized( &ret );
+ if ( !ret )
+ {
+ MPI_Finalize( );
+ }
+
+ if ( ret == PUREMD_SUCCESS )
+ {
+ ret = cleanup( handle );
+ }
+
+ return (ret == PUREMD_SUCCESS) ? 0 : (-1);
+}
diff --git a/PG-PuReMD/src/ffield.c b/PG-PuReMD/src/ffield.c
index 96d2e6b22de0fa56e1cfac417466d58905213cc2..0878882172cfcbe2adcfa775b2bdd905089ed8ca 100644
--- a/PG-PuReMD/src/ffield.c
+++ b/PG-PuReMD/src/ffield.c
@@ -30,7 +30,7 @@
#endif
-void Read_Force_Field( char *ffield_file, reax_interaction *reax,
+void Read_Force_Field_File( char *ffield_file, reax_interaction *reax,
reax_system *system, control_params *control )
{
FILE *fp;
diff --git a/PG-PuReMD/src/ffield.h b/PG-PuReMD/src/ffield.h
index 6b11374252f3922a59268071aa700ecfa6d61e5f..d847366fc115fec9df6c16de34eca3239ef09eef 100644
--- a/PG-PuReMD/src/ffield.h
+++ b/PG-PuReMD/src/ffield.h
@@ -29,7 +29,7 @@
extern "C" {
#endif
-void Read_Force_Field( char*, reax_interaction*, reax_system *, control_params* );
+void Read_Force_Field_File( char*, reax_interaction*, reax_system *, control_params* );
#ifdef __cplusplus
}
diff --git a/PG-PuReMD/src/geo_tools.c b/PG-PuReMD/src/geo_tools.c
index ecc1866cde7c9208dc5bc4e5ebac3aa7bcc83b28..24c5092aa0fe4090a40cc06e92a91ce62248357f 100644
--- a/PG-PuReMD/src/geo_tools.c
+++ b/PG-PuReMD/src/geo_tools.c
@@ -67,7 +67,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
}
-void Read_Geo( char* geo_file, reax_system* system, control_params *control,
+void Read_Geo_File( char* geo_file, reax_system* system, control_params *control,
simulation_data *data, storage *workspace, mpi_datatypes *mpi_data )
{
int i, j, serial, top;
@@ -79,7 +79,7 @@ void Read_Geo( char* geo_file, reax_system* system, control_params *control,
reax_atom *atom;
/* open the geometry file */
- geo = sfopen( geo_file, "r", "Read_Geo::geo" );
+ geo = sfopen( geo_file, "r", "Read_Geo_File::geo" );
/* read box information */
fscanf( geo, CUSTOM_BOXGEO_FORMAT,
@@ -136,7 +136,7 @@ void Read_Geo( char* geo_file, reax_system* system, control_params *control,
}
}
- sfclose( geo, "Read_Geo::geo" );
+ sfclose( geo, "Read_Geo_File::geo" );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: finished reading the geo file\n", system->my_rank );
@@ -247,7 +247,7 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
}
-void Read_PDB( char* pdb_file, reax_system* system, control_params *control,
+void Read_PDB_File( char* pdb_file, reax_system* system, control_params *control,
simulation_data *data, storage *workspace, mpi_datatypes *mpi_data )
{
FILE *pdb;
@@ -264,7 +264,7 @@ void Read_PDB( char* pdb_file, reax_system* system, control_params *control,
rvec x;
reax_atom *atom;
- pdb = sfopen( pdb_file, "r", "Read_PDB::pdb" );
+ pdb = sfopen( pdb_file, "r", "Read_PDB_File::pdb" );
/* allocate memory for tokenizing pdb lines */
Allocate_Tokenizer_Space( &s, &s1, &tmp );
@@ -451,14 +451,14 @@ void Read_PDB( char* pdb_file, reax_system* system, control_params *control,
}
}
- sfclose( pdb, "Read_PDB::pdb" );
+ sfclose( pdb, "Read_PDB_File::pdb" );
}
/* PDB serials are written without regard to the order, we'll see if this
* cause trouble, if so we'll have to rethink this approach
* Also, we do not write connect lines yet. */
-void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
+void Write_PDB_File( reax_system* system, reax_list* bonds, simulation_data *data,
control_params *control, mpi_datatypes *mpi_data, output_controls *out_control )
{
int i, cnt, me, np, buffer_req, buffer_len;
@@ -477,7 +477,7 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
np = control->nprocs;
/* Allocation*/
- line = smalloc( sizeof(char) * PDB_ATOM_FORMAT_O_LENGTH, "Write_PDB::line" );
+ line = smalloc( sizeof(char) * PDB_ATOM_FORMAT_O_LENGTH, "Write_PDB_File::line" );
if ( me == MASTER_NODE )
{
buffer_req = system->bigN * PDB_ATOM_FORMAT_O_LENGTH;
@@ -487,7 +487,7 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
buffer_req = system->n * PDB_ATOM_FORMAT_O_LENGTH;
}
- buffer = smalloc( sizeof(char) * buffer_req, "Write_PDB::buffer" );
+ buffer = smalloc( sizeof(char) * buffer_req, "Write_PDB_File::buffer" );
pdb = NULL;
line[0] = 0;
@@ -512,7 +512,7 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
sprintf( fname, "%s-%d.pdb", control->sim_name, data->step );
- pdb = sfopen( fname, "w", "Write_PDB::pdb" );
+ pdb = sfopen( fname, "w", "Write_PDB_File::pdb" );
fprintf( pdb, PDB_CRYST1_FORMAT_O,
"CRYST1",
system->big_box.box_norms[0], system->big_box.box_norms[1],
@@ -563,7 +563,7 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
if ( me == MASTER_NODE)
{
fprintf( pdb, "%s", buffer );
- sfclose( pdb, "Write_PDB::pdb" );
+ sfclose( pdb, "Write_PDB_File::pdb" );
}
/* Writing connect information */
@@ -585,6 +585,6 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
}
*/
- sfree( buffer, "Write_PDB::buffer" );
- sfree( line, "Write_PDB::line" );
+ sfree( buffer, "Write_PDB_File::buffer" );
+ sfree( line, "Write_PDB_File::line" );
}
diff --git a/PG-PuReMD/src/geo_tools.h b/PG-PuReMD/src/geo_tools.h
index 16f2b823e666929499aadd393e01f7ae108328b5..dc8f1adc9dc73f602fa5deebe5939d6d7d8b45ec 100644
--- a/PG-PuReMD/src/geo_tools.h
+++ b/PG-PuReMD/src/geo_tools.h
@@ -115,13 +115,13 @@ COLUMNS DATA TYPE FIELD DEFINITION
extern "C" {
#endif
-void Read_Geo( char*, reax_system*, control_params*,
+void Read_Geo_File( char*, reax_system*, control_params*,
simulation_data*, storage*, mpi_datatypes* );
-void Read_PDB( char*, reax_system*, control_params*,
+void Read_PDB_File( char*, reax_system*, control_params*,
simulation_data*, storage*, mpi_datatypes* );
-void Write_PDB( reax_system*, reax_list*, simulation_data*,
+void Write_PDB_File( reax_system*, reax_list*, simulation_data*,
control_params*, mpi_datatypes*, output_controls* );
#ifdef __cplusplus
diff --git a/PG-PuReMD/src/parallelreax.c b/PG-PuReMD/src/parallelreax.c
deleted file mode 100644
index 58a95d403746a145f5393cbed36d1a812d79e6f6..0000000000000000000000000000000000000000
--- a/PG-PuReMD/src/parallelreax.c
+++ /dev/null
@@ -1,563 +0,0 @@
-/*----------------------------------------------------------------------
- PuReMD - Purdue ReaxFF Molecular Dynamics Program
-
- Copyright (2010) Purdue University
- Hasan Metin Aktulga, haktulga@cs.purdue.edu
- Joseph Fogarty, jcfogart@mail.usf.edu
- Sagar Pandit, pandit@usf.edu
- Ananth Y Grama, ayg@cs.purdue.edu
-
- This program is free software; you can redistribute it and/or
- modify it under the terms of the GNU General Public License as
- published by the Free Software Foundation; either version 2 of
- the License, or (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
- See the GNU General Public License for more details:
- <http://www.gnu.org/licenses/>.
- ----------------------------------------------------------------------*/
-
-#include "reax_types.h"
-
-#include "allocate.h"
-#include "analyze.h"
-#include "comm_tools.h"
-#include "control.h"
-#include "ffield.h"
-#include "forces.h"
-#include "geo_tools.h"
-#include "init_md.h"
-#include "integrate.h"
-#include "io_tools.h"
-#include "neighbors.h"
-#include "reset_tools.h"
-#include "restart.h"
-#include "system_props.h"
-#include "tool_box.h"
-#include "traj.h"
-#include "vector.h"
-
-#ifdef HAVE_CUDA
- #include "cuda/cuda_copy.h"
- #include "cuda/cuda_environment.h"
- #include "cuda/cuda_forces.h"
- #include "cuda/cuda_init_md.h"
- #include "cuda/cuda_neighbors.h"
- #include "cuda/cuda_post_evolve.h"
- #include "cuda/cuda_reset_tools.h"
- #include "cuda/cuda_system_props.h"
- #include "cuda/cuda_utils.h"
- #if defined(DEBUG)
- #include "cuda/cuda_validation.h"
- #endif
-#endif
-
-/* CUDA-specific globals */
-reax_list **dev_lists;
-storage *dev_workspace;
-void *scratch;
-void *host_scratch;
-int BLOCKS, BLOCKS_POW_2, BLOCK_SIZE;
-int BLOCKS_N, BLOCKS_POW_2_N;
-int MATVEC_BLOCKS;
-
-
-void Read_Control_Files( char *geo_file, char *ffield_file, char *control_file,
- reax_system *system, control_params *control, simulation_data *data,
- storage *workspace, output_controls *out_control, mpi_datatypes *mpi_data )
-{
- Read_Force_Field( ffield_file, &(system->reax_param), system, control );
-
- Read_Control_File( control_file, control, out_control );
-
- if ( control->geo_format == CUSTOM )
- {
- Read_Geo( geo_file, system, control, data, workspace, mpi_data );
- }
- else if ( control->geo_format == PDB )
- {
- Read_PDB( geo_file, system, control, data, workspace, mpi_data );
- }
- else if ( control->geo_format == ASCII_RESTART )
- {
- Read_Restart( geo_file, system, control, data, workspace, mpi_data );
- control->restart = 1;
- }
- else if ( control->geo_format == BINARY_RESTART )
- {
- Read_Binary_Restart( geo_file, system, control, data, workspace, mpi_data );
- control->restart = 1;
- }
- else
- {
- fprintf( stderr, "unknown geo file format. terminating!\n" );
- MPI_Abort( MPI_COMM_WORLD, INVALID_GEO );
- }
-}
-
-
-void Post_Evolve( reax_system* system, control_params* control,
- simulation_data* data, storage* workspace, reax_list** lists,
- output_controls *out_control, mpi_datatypes *mpi_data )
-{
- int i;
- rvec diff, cross;
-
- /* remove translational and rotational velocity of the center of mass from system */
- if ( control->ensemble != NVE && control->remove_CoM_vel &&
- data->step % control->remove_CoM_vel == 0 )
- {
- /* compute velocity of the center of mass */
- Compute_Center_of_Mass( system, data, mpi_data, mpi_data->comm_mesh3D );
-
- for ( i = 0; i < system->n; i++ )
- {
- /* remove translational term */
- rvec_ScaledAdd( system->my_atoms[i].v, -1.0, data->vcm );
-
- /* remove rotational term */
- rvec_ScaledSum( diff, 1.0, system->my_atoms[i].x, -1.0, data->xcm );
- rvec_Cross( cross, data->avcm, diff );
- rvec_ScaledAdd( system->my_atoms[i].v, -1.0, cross );
- }
- }
-
- /* compute kinetic energy of system */
- Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
-
- if ( (out_control->energy_update_freq > 0
- && data->step % out_control->energy_update_freq == 0)
- || (out_control->write_steps > 0
- && data->step % out_control->write_steps == 0) )
- {
- Compute_Total_Energy( system, data, MPI_COMM_WORLD );
- }
-}
-
-
-#ifdef HAVE_CUDA
-int Cuda_Post_Evolve( reax_system* system, control_params* control,
- simulation_data* data, storage* workspace, reax_list** lists,
- output_controls *out_control, mpi_datatypes *mpi_data )
-{
- /* remove trans & rot velocity of the center of mass from system */
- if ( control->ensemble != NVE && control->remove_CoM_vel &&
- data->step % control->remove_CoM_vel == 0 )
- {
- /* compute velocity of the center of mass */
- Cuda_Compute_Center_of_Mass( system, data, mpi_data, mpi_data->comm_mesh3D );
-
- post_evolve_velocities( system, data );
- }
-
- /* compute kinetic energy of the system */
- Cuda_Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
-
- return SUCCESS;
-}
-#endif
-
-
-static void usage( char* argv[] )
-{
- fprintf( stderr, "usage: ./%s geometry_file force_field_params_file control_file\n", argv[0] );
-}
-
-
-int main( int argc, char* argv[] )
-{
- reax_system *system;
- control_params *control;
- simulation_data *data;
- storage *workspace;
- reax_list **lists;
- output_controls *out_control;
- mpi_datatypes *mpi_data;
- int i, ret, retries;
- real t_start = 0, t_elapsed;
-#if defined(DEBUG)
- real t_begin, t_end;
-#endif
-
- MPI_Init( &argc, &argv );
-
- if ( argc != 4 )
- {
- usage( argv );
- MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
- }
-
-#ifdef HAVE_CUDA
- system = smalloc( sizeof(reax_system), "main::system" );
- control = smalloc( sizeof(control_params), "main::control" );
- data = smalloc( sizeof(simulation_data), "main::data" );
- workspace = smalloc( sizeof(storage), "main::workspace" );
- lists = smalloc( sizeof(reax_list *) * LIST_N, "main::lists" );
- for ( i = 0; i < LIST_N; ++i )
- {
- lists[i] = smalloc( sizeof(reax_list), "main::lists[i]" );
- lists[i]->allocated = FALSE;
- }
- out_control = smalloc( sizeof(output_controls), "main::out_control" );
- mpi_data = smalloc( sizeof(mpi_datatypes), "main::mpi_data" );
-
- /* allocate auxiliary data structures (GPU) */
- dev_workspace = smalloc( sizeof(storage), "main::dev_workspace" );
- dev_lists = smalloc ( sizeof(reax_list *) * LIST_N, "main::dev_lists" );
- for ( i = 0; i < LIST_N; ++i )
- {
- dev_lists[i] = smalloc( sizeof(reax_list), "main::dev_lists[i]" );
- dev_lists[i]->allocated = FALSE;
- }
-
- /* setup MPI environment */
- MPI_Comm_size( MPI_COMM_WORLD, &(control->nprocs) );
- MPI_Comm_rank( MPI_COMM_WORLD, &(system->my_rank) );
-
- /* read system config files */
- Read_Control_Files( argv[1], argv[2], argv[3], system, control,
- data, workspace, out_control, mpi_data );
-
- /* setup the CUDA Device for this process */
- Setup_Cuda_Environment( system->my_rank, control->nprocs, control->gpus_per_node );
-
-#if defined(DEBUG)
- print_device_mem_usage( );
-#endif
-
- /* init blocks sizes */
- init_blocks( system );
-
- if ( system->my_rank == MASTER_NODE )
- {
- t_start = Get_Time( );
- }
-
- Cuda_Initialize( system, control, data, workspace, lists, out_control, mpi_data );
-
-#if defined(__CUDA_DEBUG__)
- Pure_Initialize( system, control, data, workspace, lists, out_control, mpi_data );
-#endif
-
-#if defined(DEBUG)
- print_device_mem_usage( );
-#endif
-
- /* init the blocks sizes for cuda kernels */
- init_blocks( system );
-
- /* compute f_0 */
- Comm_Atoms( system, control, data, workspace, mpi_data, TRUE );
- Sync_Atoms( system );
- Sync_Grid( &system->my_grid, &system->d_my_grid );
- init_blocks( system );
-
- Cuda_Reset( system, control, data, workspace, lists );
-
-#if defined(__CUDA_DEBUG__)
- Reset( system, control, data, workspace, lists );
-#endif
-
- Cuda_Generate_Neighbor_Lists( system, data, workspace, lists );
-
-#if defined(__CUDA_DEBUG__)
- Generate_Neighbor_Lists( system, data, workspace, lists );
-#endif
-
-#if defined(__CUDA_DEBUG__)
- Compute_Forces( system, control, data, workspace,
- lists, out_control, mpi_data );
-#endif
-
- Cuda_Compute_Forces( system, control, data, workspace, lists,
- out_control, mpi_data );
-
-#if defined (__CUDA_DEBUG__)
- Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
-#endif
-
- Cuda_Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
-
-#if defined(__CUDA_DEBUG__)
- validate_device( system, data, workspace, lists );
-#endif
-
-#if !defined(__CUDA_DEBUG__)
- Output_Results( system, control, data, lists, out_control, mpi_data );
-#endif
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d: step%d completed\n", system->my_rank, data->step );
- MPI_Barrier( MPI_COMM_WORLD );
-#endif
-
- ++data->step;
- retries = 0;
- while ( data->step <= control->nsteps && retries < MAX_RETRIES )
- {
- ret = SUCCESS;
-
- if ( control->T_mode && retries == 0 )
- {
- Temperature_Control( control, data );
- }
-
-#if defined(DEBUG)
- t_begin = Get_Time();
-#endif
-
-#if defined(__CUDA_DEBUG__)
- ret = control->Evolve( system, control, data, workspace,
- lists, out_control, mpi_data );
-#endif
-
- ret = control->Cuda_Evolve( system, control, data, workspace,
- lists, out_control, mpi_data );
-
-#if defined(DEBUG)
- t_end = Get_Timing_Info( t_begin );
- fprintf( stderr, " Evolve time: %f \n", t_end );
-#endif
-
-#if defined(DEBUG)
- t_begin = Get_Time();
-#endif
-
- if ( ret == SUCCESS )
- {
- ret = Cuda_Post_Evolve( system, control, data, workspace, lists,
- out_control, mpi_data );
- }
-
-#if defined(__CUDA_DEBUG__)
- Post_Evolve(system, control, data, workspace, lists, out_control, mpi_data);
-#endif
-
-#if defined(DEBUG)
- t_end = Get_Timing_Info( t_begin );
- fprintf( stderr, " Post Evolve time: %f \n", t_end );
-#endif
-
- if ( ret == SUCCESS )
- {
- data->timing.num_retries = retries;
-
-#if !defined(__CUDA_DEBUG__)
- Output_Results( system, control, data, lists, out_control, mpi_data );
-#endif
-
-// Analysis(system, control, data, workspace, lists, out_control, mpi_data);
-
- /* dump restart info */
-// if ( out_control->restart_freq &&
-// (data->step-data->prev_steps) % out_control->restart_freq == 0 )
-// {
-// if( out_control->restart_format == WRITE_ASCII )
-// {
-// Write_Restart( system, control, data, out_control, mpi_data );
-// }
-// else if( out_control->restart_format == WRITE_BINARY )
-// {
-// Write_Binary_Restart( system, control, data, out_control, mpi_data );
-// }
-// }
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d: step%d completed\n", system->my_rank, data->step );
- MPI_Barrier( MPI_COMM_WORLD );
-#endif
-
- ++data->step;
- retries = 0;
- }
- else
- {
- ++retries;
-
-#if defined(DEBUG)
- fprintf( stderr, "[INFO] p%d: retrying step %d...\n", system->my_rank, data->step );
-#endif
- }
- }
-
- if ( retries >= MAX_RETRIES )
- {
- fprintf( stderr, "[ERROR] Maximum retries reached for this step (%d). Terminating...\n",
- retries );
- MPI_Abort( MPI_COMM_WORLD, MAX_RETRIES_REACHED );
- }
-
-#if defined(__CUDA_DEBUG__)
- /* vaildate the results in debug mode */
- validate_device( system, data, workspace, lists );
-#endif
-
-#else
- system = smalloc( sizeof(reax_system), "main::system" );
- control = smalloc( sizeof(control_params), "main::control" );
- data = smalloc( sizeof(simulation_data), "main::data" );
- workspace = smalloc( sizeof(storage), "main::workspace" );
- lists = smalloc( sizeof(reax_list *) * LIST_N, "main::lists" );
- for ( i = 0; i < LIST_N; ++i )
- {
- lists[i] = smalloc( sizeof(reax_list), "main::lists[i]" );
- lists[i]->allocated = FALSE;
- }
- out_control = smalloc( sizeof(output_controls), "main::out_control" );
- mpi_data = smalloc( sizeof(mpi_datatypes), "main::mpi_data" );
-
- /* setup MPI environment */
- MPI_Comm_size( MPI_COMM_WORLD, &control->nprocs );
- MPI_Comm_rank( MPI_COMM_WORLD, &system->my_rank );
-
- /* read config files */
- Read_Control_Files( argv[1], argv[2], argv[3], system, control,
- data, workspace, out_control, mpi_data );
-
- if ( system->my_rank == MASTER_NODE )
- {
- t_start = Get_Time( );
- }
-
- Initialize( system, control, data, workspace, lists, out_control, mpi_data );
-
- /* compute f_0 */
- Comm_Atoms( system, control, data, workspace, mpi_data, TRUE );
-
- Reset( system, control, data, workspace, lists );
-
- if ( ret == FAILURE )
- {
- ReAllocate( system, control, data, workspace, lists, mpi_data );
- }
-
- ret = Generate_Neighbor_Lists( system, data, workspace, lists );
-
- if ( ret != SUCCESS )
- {
- fprintf( stderr, "[ERROR] cannot generate initial neighbor lists. Terminating...\n" );
- MPI_Abort( MPI_COMM_WORLD, CANNOT_INITIALIZE );
- }
-
- ret = Compute_Forces( system, control, data, workspace, lists, out_control, mpi_data );
-
- if ( ret != SUCCESS )
- {
- fprintf( stderr, "[ERROR] cannot compute initial forces. Terminating...\n" );
- MPI_Abort( MPI_COMM_WORLD, CANNOT_INITIALIZE );
- }
-
- Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
-
- Compute_Total_Energy( system, data, MPI_COMM_WORLD );
-
- Output_Results( system, control, data, lists, out_control, mpi_data );
-
- Check_Energy( data );
-
- retries = 0;
- ++data->step;
- while ( data->step <= control->nsteps && retries < MAX_RETRIES )
- {
- ret = SUCCESS;
-
- if ( control->T_mode && retries == 0 )
- {
- Temperature_Control( control, data );
- }
-
- ret = control->Evolve( system, control, data, workspace,
- lists, out_control, mpi_data );
-
- if ( ret == SUCCESS )
- {
- Post_Evolve(system, control, data, workspace, lists, out_control, mpi_data);
- }
-
- if ( ret == SUCCESS )
- {
- data->timing.num_retries = retries;
-
- Output_Results( system, control, data, lists, out_control, mpi_data );
-
-// Analysis(system, control, data, workspace, lists, out_control, mpi_data);
-
- /* dump restart info */
- if ( out_control->restart_freq &&
- (data->step - data->prev_steps) % out_control->restart_freq == 0 )
- {
- if ( out_control->restart_format == WRITE_ASCII )
- {
- Write_Restart( system, control, data, out_control, mpi_data );
- }
- else if ( out_control->restart_format == WRITE_BINARY )
- {
- Write_Binary_Restart( system, control, data, out_control, mpi_data );
- }
- }
-
- Check_Energy( data );
-
- ++data->step;
- retries = 0;
- }
- else
- {
- ++retries;
-
-#if defined(DEBUG)
- fprintf( stderr, "[INFO] p%d: retrying step %d...\n", system->my_rank, data->step );
-#endif
- }
- }
-
- if ( retries >= MAX_RETRIES )
- {
- fprintf( stderr, "[ERROR] Maximum retries reached for this step (%d). Terminating...\n",
- retries );
- MPI_Abort( MPI_COMM_WORLD, MAX_RETRIES_REACHED );
- }
-
-#endif
-
- /* end of the simulation, write total simulation time */
- if ( system->my_rank == MASTER_NODE )
- {
- t_elapsed = Get_Timing_Info( t_start );
- fprintf( out_control->out, "Total Simulation Time: %.2f secs\n", t_elapsed );
- }
-
-// Write_PDB( &system, &lists[BONDS], &out_control );
-
- Close_Output_Files( system, control, out_control, mpi_data );
-
-#if defined(TEST_ENERGY) || defined(TEST_FORCES)
-// Integrate_Results(control);
-#endif
-
-#if defined(DEBUG)
- fprintf( stderr, "p%d has reached the END\n", system->my_rank );
-#endif
-
- MPI_Finalized( &ret );
- if ( !ret )
- {
- MPI_Finalize( );
- }
-
- sfree( mpi_data, "main::mpi_data" );
- sfree( out_control, "main::out_control" );
- for ( i = 0; i < LIST_N; ++i )
- {
- sfree( lists[i], "main::lists[i]" );
- }
- sfree( lists, "main::lists" );
- sfree( workspace, "main::workspace" );
- sfree( data, "main::data" );
- sfree( control, "main::control" );
- sfree( system, "main::system" );
-
- return 0;
-}
diff --git a/PG-PuReMD/src/puremd.c b/PG-PuReMD/src/puremd.c
new file mode 100644
index 0000000000000000000000000000000000000000..9439ac9696ccbe3ad4de222dffb279ee304d18b8
--- /dev/null
+++ b/PG-PuReMD/src/puremd.c
@@ -0,0 +1,649 @@
+/*----------------------------------------------------------------------
+ PuReMD - Purdue ReaxFF Molecular Dynamics Program
+
+ Copyright (2010) Purdue University
+ Hasan Metin Aktulga, haktulga@cs.purdue.edu
+ Joseph Fogarty, jcfogart@mail.usf.edu
+ Sagar Pandit, pandit@usf.edu
+ Ananth Y Grama, ayg@cs.purdue.edu
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU General Public License as
+ published by the Free Software Foundation; either version 2 of
+ the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ See the GNU General Public License for more details:
+ <http://www.gnu.org/licenses/>.
+ ----------------------------------------------------------------------*/
+
+#include "puremd.h"
+
+#include "allocate.h"
+#include "analyze.h"
+#include "comm_tools.h"
+#include "control.h"
+#include "ffield.h"
+#include "forces.h"
+#include "geo_tools.h"
+#include "init_md.h"
+#include "integrate.h"
+#include "io_tools.h"
+#include "neighbors.h"
+#include "reset_tools.h"
+#include "restart.h"
+#include "system_props.h"
+#include "tool_box.h"
+#include "traj.h"
+#include "vector.h"
+
+#ifdef HAVE_CUDA
+ #include "cuda/cuda_copy.h"
+ #include "cuda/cuda_environment.h"
+ #include "cuda/cuda_forces.h"
+ #include "cuda/cuda_init_md.h"
+ #include "cuda/cuda_neighbors.h"
+ #include "cuda/cuda_post_evolve.h"
+ #include "cuda/cuda_reset_tools.h"
+ #include "cuda/cuda_system_props.h"
+ #include "cuda/cuda_utils.h"
+ #if defined(DEBUG)
+ #include "cuda/cuda_validation.h"
+ #endif
+#endif
+
+/* CUDA-specific globals */
+reax_list **dev_lists;
+storage *dev_workspace;
+void *scratch;
+void *host_scratch;
+int BLOCKS, BLOCKS_POW_2, BLOCK_SIZE;
+int BLOCKS_N, BLOCKS_POW_2_N;
+int MATVEC_BLOCKS;
+
+
+static void Read_Config_Files( const char * const geo_file, const char * const ffield_file,
+ const char * const control_file,
+ reax_system *system, control_params *control, simulation_data *data,
+ storage *workspace, output_controls *out_control, mpi_datatypes *mpi_data )
+{
+ Read_Force_Field_File( ffield_file, &system->reax_param, system, control );
+
+ Read_Control_File( control_file, control, out_control );
+
+ if ( control->geo_format == CUSTOM )
+ {
+ Read_Geo_File( geo_file, system, control, data, workspace, mpi_data );
+ }
+ else if ( control->geo_format == PDB )
+ {
+ Read_PDB_File( geo_file, system, control, data, workspace, mpi_data );
+ }
+ else if ( control->geo_format == ASCII_RESTART )
+ {
+ Read_Restart_File( geo_file, system, control, data, workspace, mpi_data );
+ control->restart = 1;
+ }
+ else if ( control->geo_format == BINARY_RESTART )
+ {
+ Read_Binary_Restart_File( geo_file, system, control, data, workspace, mpi_data );
+ control->restart = 1;
+ }
+ else
+ {
+ fprintf( stderr, "[ERROR] unknown geo file format. terminating!\n" );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_GEO );
+ }
+}
+
+
+static void Post_Evolve( reax_system* system, control_params* control,
+ simulation_data* data, storage* workspace, reax_list** lists,
+ output_controls *out_control, mpi_datatypes *mpi_data )
+{
+ int i;
+ rvec diff, cross;
+
+ /* remove translational and rotational velocity of the center of mass from system */
+ if ( control->ensemble != NVE && control->remove_CoM_vel &&
+ data->step % control->remove_CoM_vel == 0 )
+ {
+ /* compute velocity of the center of mass */
+ Compute_Center_of_Mass( system, data, mpi_data, mpi_data->comm_mesh3D );
+
+ for ( i = 0; i < system->n; i++ )
+ {
+ /* remove translational term */
+ rvec_ScaledAdd( system->my_atoms[i].v, -1.0, data->vcm );
+
+ /* remove rotational term */
+ rvec_ScaledSum( diff, 1.0, system->my_atoms[i].x, -1.0, data->xcm );
+ rvec_Cross( cross, data->avcm, diff );
+ rvec_ScaledAdd( system->my_atoms[i].v, -1.0, cross );
+ }
+ }
+
+ /* compute kinetic energy of system */
+ Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+
+ if ( (out_control->energy_update_freq > 0
+ && data->step % out_control->energy_update_freq == 0)
+ || (out_control->write_steps > 0
+ && data->step % out_control->write_steps == 0) )
+ {
+ Compute_Total_Energy( system, data, MPI_COMM_WORLD );
+ }
+}
+
+
+#ifdef HAVE_CUDA
+static void Cuda_Post_Evolve( reax_system* system, control_params* control,
+ simulation_data* data, storage* workspace, reax_list** lists,
+ output_controls *out_control, mpi_datatypes *mpi_data )
+{
+ /* remove trans & rot velocity of the center of mass from system */
+ if ( control->ensemble != NVE && control->remove_CoM_vel &&
+ data->step % control->remove_CoM_vel == 0 )
+ {
+ /* compute velocity of the center of mass */
+ Cuda_Compute_Center_of_Mass( system, data, mpi_data, mpi_data->comm_mesh3D );
+
+ post_evolve_velocities( system, data );
+ }
+
+ /* compute kinetic energy of the system */
+ Cuda_Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+}
+#endif
+
+
+void* setup( const char * const geo_file, const char * const ffield_file,
+ const char * const control_file )
+{
+ int i;
+ puremd_handle *pmd_handle;
+
+ /* top-level allocation */
+ pmd_handle = (puremd_handle*) smalloc( sizeof(puremd_handle),
+ "setup::pmd_handle" );
+
+ /* second-level allocations */
+ pmd_handle->system = smalloc( sizeof(reax_system),
+ "Setup::pmd_handle->system" );
+ pmd_handle->control = smalloc( sizeof(control_params),
+ "Setup::pmd_handle->control" );
+ pmd_handle->data = smalloc( sizeof(simulation_data),
+ "Setup::pmd_handle->data" );
+ pmd_handle->workspace = smalloc( sizeof(storage),
+ "Setup::pmd_handle->workspace" );
+ pmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
+ "Setup::pmd_handle->lists" );
+ for ( i = 0; i < LIST_N; ++i )
+ {
+ pmd_handle->lists[i] = smalloc( sizeof(reax_list),
+ "Setup::pmd_handle->lists[i]" );
+ pmd_handle->lists[i]->allocated = FALSE;
+ }
+ pmd_handle->out_control = smalloc( sizeof(output_controls),
+ "Setup::pmd_handle->out_control" );
+ pmd_handle->mpi_data = smalloc( sizeof(mpi_datatypes),
+ "Setup::pmd_handle->mpi_data" );
+
+#ifdef HAVE_CUDA
+ /* allocate auxiliary data structures (GPU) */
+ dev_workspace = smalloc( sizeof(storage), "Setup::dev_workspace" );
+ dev_lists = smalloc ( sizeof(reax_list *) * LIST_N, "Setup::dev_lists" );
+ for ( i = 0; i < LIST_N; ++i )
+ {
+ dev_lists[i] = smalloc( sizeof(reax_list), "Setup::dev_lists[i]" );
+ dev_lists[i]->allocated = FALSE;
+ }
+#endif
+
+ pmd_handle->output_enabled = TRUE;
+ pmd_handle->callback = NULL;
+
+ /* setup MPI environment */
+ MPI_Comm_size( MPI_COMM_WORLD, &pmd_handle->control->nprocs );
+ MPI_Comm_rank( MPI_COMM_WORLD, &pmd_handle->system->my_rank );
+
+ /* read system config files */
+ Read_Config_Files( geo_file, ffield_file, control_file,
+ pmd_handle->system, pmd_handle->control, pmd_handle->data,
+ pmd_handle->workspace, pmd_handle->out_control, pmd_handle->mpi_data );
+
+#ifdef HAVE_CUDA
+ /* setup the CUDA Device for this process */
+ Setup_Cuda_Environment( system->my_rank, control->nprocs, control->gpus_per_node );
+
+#if defined(DEBUG)
+ print_device_mem_usage( );
+#endif
+
+ /* init blocks sizes */
+ init_blocks( system );
+#endif
+
+ return (void*) pmd_handle;
+}
+
+
+int setup_callback( const void * const handle, const callback_function callback )
+{
+ int ret;
+ puremd_handle *pmd_handle;
+
+
+ ret = PUREMD_FAILURE;
+
+ if ( handle != NULL && callback != NULL )
+ {
+ pmd_handle = (puremd_handle*) handle;
+ pmd_handle->callback = callback;
+ ret = PUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+int simulate( const void * const handle )
+{
+ int ret, retries;
+ reax_system *system;
+ control_params *control;
+ simulation_data *data;
+ storage *workspace;
+ reax_list **lists;
+ output_controls *out_control;
+ mpi_datatypes *mpi_data;
+ puremd_handle *pmd_handle;
+ real t_start = 0, t_elapsed;
+#if defined(DEBUG)
+ real t_begin, t_end;
+#endif
+
+ ret = PUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ pmd_handle = (puremd_handle*) handle;
+
+ system = pmd_handle->system;
+ control = pmd_handle->control;
+ data = pmd_handle->data;
+ workspace = pmd_handle->workspace;
+ lists = pmd_handle->lists;
+ out_control = pmd_handle->out_control;
+ mpi_data = pmd_handle->mpi_data;
+
+#ifdef HAVE_CUDA
+ if ( system->my_rank == MASTER_NODE )
+ {
+ t_start = Get_Time( );
+ }
+
+ Cuda_Initialize( system, control, data, workspace, lists, out_control, mpi_data );
+
+#if defined(__CUDA_DEBUG__)
+ Pure_Initialize( system, control, data, workspace, lists, out_control, mpi_data );
+#endif
+
+#if defined(DEBUG)
+ print_device_mem_usage( );
+#endif
+
+ /* init the blocks sizes for cuda kernels */
+ init_blocks( system );
+
+ /* compute f_0 */
+ Comm_Atoms( system, control, data, workspace, mpi_data, TRUE );
+ Sync_Atoms( system );
+ Sync_Grid( &system->my_grid, &system->d_my_grid );
+ init_blocks( system );
+
+ Cuda_Reset( system, control, data, workspace, lists );
+
+#if defined(__CUDA_DEBUG__)
+ Reset( system, control, data, workspace, lists );
+#endif
+
+ Cuda_Generate_Neighbor_Lists( system, data, workspace, lists );
+
+#if defined(__CUDA_DEBUG__)
+ Generate_Neighbor_Lists( system, data, workspace, lists );
+#endif
+
+#if defined(__CUDA_DEBUG__)
+ Compute_Forces( system, control, data, workspace,
+ lists, out_control, mpi_data );
+#endif
+
+ Cuda_Compute_Forces( system, control, data, workspace, lists,
+ out_control, mpi_data );
+
+#if defined (__CUDA_DEBUG__)
+ Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+#endif
+
+ Cuda_Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+
+#if defined(__CUDA_DEBUG__)
+ validate_device( system, data, workspace, lists );
+#endif
+
+#if !defined(__CUDA_DEBUG__)
+ Output_Results( system, control, data, lists, out_control, mpi_data );
+#endif
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d: step%d completed\n", system->my_rank, data->step );
+ MPI_Barrier( MPI_COMM_WORLD );
+#endif
+
+ ++data->step;
+ retries = 0;
+ while ( data->step <= control->nsteps && retries < MAX_RETRIES )
+ {
+ ret = SUCCESS;
+
+ if ( control->T_mode && retries == 0 )
+ {
+ Temperature_Control( control, data );
+ }
+
+#if defined(DEBUG)
+ t_begin = Get_Time();
+#endif
+
+#if defined(__CUDA_DEBUG__)
+ ret = control->Evolve( system, control, data, workspace,
+ lists, out_control, mpi_data );
+#endif
+
+ ret = control->Cuda_Evolve( system, control, data, workspace,
+ lists, out_control, mpi_data );
+
+#if defined(DEBUG)
+ t_end = Get_Timing_Info( t_begin );
+ fprintf( stderr, " Evolve time: %f \n", t_end );
+#endif
+
+#if defined(DEBUG)
+ t_begin = Get_Time( );
+#endif
+
+ if ( ret == SUCCESS )
+ {
+ Cuda_Post_Evolve( system, control, data, workspace, lists,
+ out_control, mpi_data );
+
+#if defined(__CUDA_DEBUG__)
+ Post_Evolve( system, control, data, workspace, lists, out_control, mpi_data );
+#endif
+ }
+
+#if defined(DEBUG)
+ t_end = Get_Timing_Info( t_begin );
+ fprintf( stderr, " Post Evolve time: %f \n", t_end );
+#endif
+
+ if ( ret == SUCCESS )
+ {
+ data->timing.num_retries = retries;
+
+#if !defined(__CUDA_DEBUG__)
+ Output_Results( system, control, data, lists, out_control, mpi_data );
+#endif
+
+// Analysis( system, control, data, workspace, lists, out_control, mpi_data );
+
+ /* dump restart info */
+// if ( out_control->restart_freq &&
+// (data->step-data->prev_steps) % out_control->restart_freq == 0 )
+// {
+// if( out_control->restart_format == WRITE_ASCII )
+// {
+// Write_Restart_File( system, control, data, out_control, mpi_data );
+// }
+// else if( out_control->restart_format == WRITE_BINARY )
+// {
+// Write_Binary_Restart_File( system, control, data, out_control, mpi_data );
+// }
+// }
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d: step%d completed\n", system->my_rank, data->step );
+ MPI_Barrier( MPI_COMM_WORLD );
+#endif
+
+ ++data->step;
+ retries = 0;
+ }
+ else
+ {
+ ++retries;
+
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] p%d: retrying step %d...\n", system->my_rank, data->step );
+#endif
+ }
+ }
+
+ if ( retries >= MAX_RETRIES )
+ {
+ fprintf( stderr, "[ERROR] Maximum retries reached for this step (%d). Terminating...\n",
+ retries );
+ MPI_Abort( MPI_COMM_WORLD, MAX_RETRIES_REACHED );
+ }
+
+#if defined(__CUDA_DEBUG__)
+ /* vaildate the results in debug mode */
+ validate_device( system, data, workspace, lists );
+#endif
+
+#else
+ if ( system->my_rank == MASTER_NODE )
+ {
+ t_start = Get_Time( );
+ }
+
+ Initialize( system, control, data, workspace, lists, out_control, mpi_data );
+
+ /* compute f_0 */
+ Comm_Atoms( system, control, data, workspace, mpi_data, TRUE );
+
+ Reset( system, control, data, workspace, lists );
+
+ if ( ret == FAILURE )
+ {
+ ReAllocate( system, control, data, workspace, lists, mpi_data );
+ }
+
+ ret = Generate_Neighbor_Lists( system, data, workspace, lists );
+
+ if ( ret != SUCCESS )
+ {
+ fprintf( stderr, "[ERROR] cannot generate initial neighbor lists. Terminating...\n" );
+ MPI_Abort( MPI_COMM_WORLD, CANNOT_INITIALIZE );
+ }
+
+ ret = Compute_Forces( system, control, data, workspace, lists, out_control, mpi_data );
+
+ if ( ret != SUCCESS )
+ {
+ fprintf( stderr, "[ERROR] cannot compute initial forces. Terminating...\n" );
+ MPI_Abort( MPI_COMM_WORLD, CANNOT_INITIALIZE );
+ }
+
+ Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+
+ Compute_Total_Energy( system, data, MPI_COMM_WORLD );
+
+ Output_Results( system, control, data, lists, out_control, mpi_data );
+
+ Check_Energy( data );
+
+ retries = 0;
+ ++data->step;
+ while ( data->step <= control->nsteps && retries < MAX_RETRIES )
+ {
+ ret = SUCCESS;
+
+ if ( control->T_mode && retries == 0 )
+ {
+ Temperature_Control( control, data );
+ }
+
+ ret = control->Evolve( system, control, data, workspace,
+ lists, out_control, mpi_data );
+
+ if ( ret == SUCCESS )
+ {
+ Post_Evolve(system, control, data, workspace, lists, out_control, mpi_data);
+ }
+
+ if ( ret == SUCCESS )
+ {
+ data->timing.num_retries = retries;
+
+ Output_Results( system, control, data, lists, out_control, mpi_data );
+
+// Analysis( system, control, data, workspace, lists, out_control, mpi_data );
+
+ /* dump restart info */
+ if ( out_control->restart_freq &&
+ (data->step - data->prev_steps) % out_control->restart_freq == 0 )
+ {
+ if ( out_control->restart_format == WRITE_ASCII )
+ {
+ Write_Restart_File( system, control, data, out_control, mpi_data );
+ }
+ else if ( out_control->restart_format == WRITE_BINARY )
+ {
+ Write_Binary_Restart_File( system, control, data, out_control, mpi_data );
+ }
+ }
+
+ Check_Energy( data );
+
+ ++data->step;
+ retries = 0;
+ }
+ else
+ {
+ ++retries;
+
+#if defined(DEBUG)
+ fprintf( stderr, "[INFO] p%d: retrying step %d...\n", system->my_rank, data->step );
+#endif
+ }
+ }
+
+ if ( retries >= MAX_RETRIES )
+ {
+ fprintf( stderr, "[ERROR] Maximum retries reached for this step (%d). Terminating...\n",
+ retries );
+ MPI_Abort( MPI_COMM_WORLD, MAX_RETRIES_REACHED );
+ }
+
+#endif
+
+ /* end of the simulation, write total simulation time */
+ if ( system->my_rank == MASTER_NODE )
+ {
+ t_elapsed = Get_Timing_Info( t_start );
+ fprintf( out_control->out, "Total Simulation Time: %.2f secs\n", t_elapsed );
+ }
+
+// Write_PDB_File( &system, &lists[BONDS], &out_control );
+
+ Close_Output_Files( system, control, out_control, mpi_data );
+
+#if defined(TEST_ENERGY) || defined(TEST_FORCES)
+// Integrate_Results(control);
+#endif
+
+#if defined(DEBUG)
+ fprintf( stderr, "p%d has reached the END\n", system->my_rank );
+#endif
+
+ ret = PUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+int cleanup( const void * const handle )
+{
+ int i, ret;
+ puremd_handle *pmd_handle;
+
+ ret = PUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ pmd_handle = (puremd_handle*) handle;
+
+ //TODO
+// Finalize( pmd_handle->system, pmd_handle->control, pmd_handle->data,
+// pmd_handle->workspace, &pmd_handle->lists, pmd_handle->out_control,
+// pmd_handle->output_enabled );
+
+ sfree( pmd_handle->mpi_data, "cleanup::pmd_handle->mpi_data" );
+ sfree( pmd_handle->out_control, "cleanup::pmd_handle->out_control" );
+ for ( i = 0; i < LIST_N; ++i )
+ {
+ sfree( pmd_handle->lists[i], "cleanup::pmd_handle->lists[i]" );
+ }
+ sfree( pmd_handle->lists, "cleanup::pmd_handle->lists" );
+ sfree( pmd_handle->workspace, "cleanup::pmd_handle->workspace" );
+ sfree( pmd_handle->data, "cleanup::pmd_handle->data" );
+ sfree( pmd_handle->control, "cleanup::pmd_handle->control" );
+ sfree( pmd_handle->system, "cleanup::pmd_handle->system" );
+
+ sfree( pmd_handle, "cleanup::pmd_handle" );
+
+ ret = PUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+reax_atom* get_atoms( const void * const handle )
+{
+ puremd_handle *pmd_handle;
+ reax_atom *atoms;
+
+ atoms = NULL;
+
+ if ( handle != NULL )
+ {
+ pmd_handle = (puremd_handle*) handle;
+ atoms = pmd_handle->system->my_atoms;
+ }
+
+ return atoms;
+}
+
+
+int set_output_enabled( const void * const handle, const int enabled )
+{
+ int ret;
+ puremd_handle *pmd_handle;
+
+ ret = PUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ pmd_handle = (puremd_handle*) handle;
+ pmd_handle->output_enabled = enabled;
+ ret = PUREMD_SUCCESS;
+ }
+
+ return ret;
+}
diff --git a/PG-PuReMD/src/puremd.h b/PG-PuReMD/src/puremd.h
new file mode 100644
index 0000000000000000000000000000000000000000..b51f854c860cf02d427a80ff1287721f1d01758e
--- /dev/null
+++ b/PG-PuReMD/src/puremd.h
@@ -0,0 +1,54 @@
+/*----------------------------------------------------------------------
+ SerialReax - Reax Force Field Simulator
+
+ Copyright (2010) Purdue University
+ Hasan Metin Aktulga, haktulga@cs.purdue.edu
+ Joseph Fogarty, jcfogart@mail.usf.edu
+ Sagar Pandit, pandit@usf.edu
+ Ananth Y Grama, ayg@cs.purdue.edu
+
+ This program is free software; you can redistribute it and/or
+ modify it under the terms of the GNU General Public License as
+ published by the Free Software Foundation; either version 2 of
+ the License, or (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ See the GNU General Public License for more details:
+ <http://www.gnu.org/licenses/>.
+ ----------------------------------------------------------------------*/
+
+#ifndef __PUREMD_H_
+#define __PUREMD_H_
+
+#include "reax_types.h"
+
+
+#define PUREMD_SUCCESS (0)
+#define PUREMD_FAILURE (-1)
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void* setup( const char * const, const char * const,
+ const char * const );
+
+int setup_callback( const void * const, const callback_function );
+
+int simulate( const void * const );
+
+int cleanup( const void * const );
+
+reax_atom* get_atoms( const void * const );
+
+int set_output_enabled( const void * const, const int );
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif
diff --git a/PG-PuReMD/src/reax_types.h b/PG-PuReMD/src/reax_types.h
index e08d2b9b0ac8cbc67d78701d93bec03a219cbc48..57d7e414017c4882c80c594566965132289eeefa 100644
--- a/PG-PuReMD/src/reax_types.h
+++ b/PG-PuReMD/src/reax_types.h
@@ -548,6 +548,7 @@ typedef struct molecule molecule;
typedef struct LR_data LR_data;
typedef struct cubic_spline_coef cubic_spline_coef;
typedef struct LR_lookup_table LR_lookup_table;
+typedef struct puremd_handle puremd_handle;
/* function pointer definitions */
@@ -563,6 +564,8 @@ typedef real (*lookup_function)( real );
typedef void (*message_sorter)( reax_system*, int, int, int, mpi_out_data* );
/**/
typedef void (*unpacker)( reax_system*, int, void*, int, neighbor_proc*, int );
+/**/
+typedef void (*callback_function)(reax_atom*, simulation_data*, reax_list*);
/* struct definitions */
@@ -2458,6 +2461,30 @@ struct LR_lookup_table
};
+/* Handle for working with an instance of the PuReMD library */
+struct puremd_handle
+{
+ /* System info. struct pointer */
+ reax_system *system;
+ /* System struct pointer */
+ control_params *control;
+ /* Control parameters struct pointer */
+ simulation_data *data;
+ /* Internal workspace struct pointer */
+ storage *workspace;
+ /* Reax interaction list struct pointer */
+ reax_list **lists;
+ /* Output controls struct pointer */
+ output_controls *out_control;
+ /* MPI datatypes struct pointer */
+ mpi_datatypes *mpi_data;
+ /* TRUE if file I/O for simulation output enabled, FALSE otherwise */
+ int output_enabled;
+ /* Callback for getting simulation state at the end of each time step */
+ callback_function callback;
+};
+
+
/* CUDA-specific globals */
extern reax_list **dev_lists;
extern storage *dev_workspace;
diff --git a/PG-PuReMD/src/restart.c b/PG-PuReMD/src/restart.c
index 19e0e0a37caea16e728b97c202cd3641a074ca2e..73908ea80b122a1580ee7a0cbb8a41c7437b0cce 100644
--- a/PG-PuReMD/src/restart.c
+++ b/PG-PuReMD/src/restart.c
@@ -29,7 +29,7 @@
#include "vector.h"
-void Write_Binary_Restart( reax_system *system, control_params *control,
+void Write_Binary_Restart_File( reax_system *system, control_params *control,
simulation_data *data, output_controls *out_control,
mpi_datatypes *mpi_data )
{
@@ -49,7 +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 );
- fres = sfopen( fname, "wb", "Write_Binary_Restart::fres" );
+ fres = sfopen( fname, "wb", "Write_Binary_Restart_File::fres" );
/* master can write the header by itself */
res_header.step = data->step;
@@ -64,12 +64,12 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
/* master needs to allocate space for all atoms */
buffer = scalloc( system->bigN, sizeof(restart_atom),
- "Write_Binary_Restart::buffer" );
+ "Write_Binary_Restart_File::buffer" );
}
else
{
buffer = scalloc( system->n, sizeof(restart_atom),
- "Write_Binary_Restart::buffer" );
+ "Write_Binary_Restart_File::buffer" );
}
/* fill in the buffers */
@@ -108,14 +108,14 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
if ( me == MASTER_NODE )
{
fwrite( buffer, system->bigN, sizeof(restart_atom), fres );
- sfclose( fres, "Write_Binary_Restart::fres" );
+ sfclose( fres, "Write_Binary_Restart_File::fres" );
}
- sfree( buffer, "Write_Binary_Restart::buffer" );
+ sfree( buffer, "Write_Binary_Restart_File::buffer" );
}
-void Write_Restart( reax_system *system, control_params *control,
+void Write_Restart_File( reax_system *system, control_params *control,
simulation_data *data, output_controls *out_control,
mpi_datatypes *mpi_data )
{
@@ -135,7 +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 );
- fres = sfopen( fname, "w", "Write_Restart::fres" );
+ fres = sfopen( fname, "w", "Write_Restart_File::fres" );
/* write the header - only master writes it */
fprintf( fres, RESTART_HEADER,
@@ -156,7 +156,7 @@ void Write_Restart( reax_system *system, control_params *control,
buffer_req = system->n * RESTART_LINE_LEN + 1;
}
- buffer = smalloc( sizeof(char) * buffer_req, "Write_Restart::buffer" );
+ buffer = smalloc( sizeof(char) * buffer_req, "Write_Restart_File::buffer" );
line[0] = 0;
buffer[0] = 0;
@@ -200,10 +200,10 @@ void Write_Restart( reax_system *system, control_params *control,
if ( me == MASTER_NODE )
{
fprintf( fres, "%s", buffer );
- sfclose( fres, "Write_Restart::fres" );
+ sfclose( fres, "Write_Restart_File::fres" );
}
- sfree( buffer, "Write_Restart::buffer" );
- sfree( line, "Write_Restart::line" );
+ sfree( buffer, "Write_Restart_File::buffer" );
+ sfree( line, "Write_Restart_File::line" );
}
@@ -235,7 +235,7 @@ void Count_Binary_Restart_Atoms( FILE *fres, reax_system *system )
}
-void Read_Binary_Restart( char *res_file, reax_system *system,
+void Read_Binary_Restart_File( char *res_file, reax_system *system,
control_params *control, simulation_data *data,
storage *workspace, mpi_datatypes *mpi_data )
{
@@ -245,7 +245,7 @@ void Read_Binary_Restart( char *res_file, reax_system *system,
restart_atom res_atom;
reax_atom *p_atom;
- fres = sfopen( res_file, "rb", "Read_Binary_Restart::fres" );
+ fres = sfopen( res_file, "rb", "Read_Binary_Restart_File::fres" );
/* first read the header lines */
fread( &res_header, sizeof(restart_header), 1, fres );
@@ -305,7 +305,7 @@ void Read_Binary_Restart( char *res_file, reax_system *system,
}
}
- sfclose( fres, "Read_Binary_Restart::fres" );
+ sfclose( fres, "Read_Binary_Restart_File::fres" );
data->step = data->prev_steps;
// nsteps is updated based on the number of steps in the previous run
@@ -350,7 +350,7 @@ void Count_Restart_Atoms( FILE *fres, reax_system *system )
}
-void Read_Restart( char *res_file, reax_system *system,
+void Read_Restart_File( char *res_file, reax_system *system,
control_params *control, simulation_data *data,
storage *workspace, mpi_datatypes *mpi_data )
{
@@ -362,13 +362,13 @@ void Read_Restart( char *res_file, reax_system *system,
rvec x_temp, v_temp;
rtensor box;
- fres = sfopen( res_file, "r", "Read_Restart::fres" );
+ fres = sfopen( res_file, "r", "Read_Restart_File::fres" );
- s = smalloc( sizeof(char) * MAX_LINE, "Read_Restart::s" );
- tmp = smalloc( sizeof(char*) * MAX_TOKENS, "Read_Restart::tmp" );
+ s = smalloc( sizeof(char) * MAX_LINE, "Read_Restart_File::s" );
+ tmp = smalloc( sizeof(char*) * MAX_TOKENS, "Read_Restart_File::tmp" );
for (i = 0; i < MAX_TOKENS; i++)
{
- tmp[i] = smalloc( sizeof(char) * MAX_LINE, "Read_Restart::tmp[i]" );
+ tmp[i] = smalloc( sizeof(char) * MAX_LINE, "Read_Restart_File::tmp[i]" );
}
//read first header lines
@@ -486,15 +486,15 @@ void Read_Restart( char *res_file, reax_system *system,
top++;
}
}
- sfclose( fres, "Read_Restart::fres" );
+ sfclose( fres, "Read_Restart_File::fres" );
/* free memory allocations at the top */
for ( i = 0; i < MAX_TOKENS; i++ )
{
- sfree( tmp[i], "Read_Restart::tmp[i]" );
+ sfree( tmp[i], "Read_Restart_File::tmp[i]" );
}
- sfree( tmp, "Read_Restart::tmp" );
- sfree( s, "Read_Restart::s" );
+ sfree( tmp, "Read_Restart_File::tmp" );
+ sfree( s, "Read_Restart_File::s" );
data->step = data->prev_steps;
// nsteps is updated based on the number of steps in the previous run
diff --git a/PG-PuReMD/src/restart.h b/PG-PuReMD/src/restart.h
index 3d13a5a17256457829224a20f32c4d9f7305a06b..9f37e641c0c982dacaeb09e78f11afa2ec94c653 100644
--- a/PG-PuReMD/src/restart.h
+++ b/PG-PuReMD/src/restart.h
@@ -45,16 +45,16 @@
extern "C" {
#endif
-void Write_Binary_Restart( reax_system*, control_params*,
+void Write_Binary_Restart_File( reax_system*, control_params*,
simulation_data*, output_controls*, mpi_datatypes* );
-void Write_Restart( reax_system*, control_params*,
+void Write_Restart_File( reax_system*, control_params*,
simulation_data*, output_controls*, mpi_datatypes* );
-void Read_Binary_Restart( char*, reax_system*, control_params*,
+void Read_Binary_Restart_File( char*, reax_system*, control_params*,
simulation_data*, storage*, mpi_datatypes* );
-void Read_Restart( char*, reax_system*, control_params*,
+void Read_Restart_File( char*, reax_system*, control_params*,
simulation_data*, storage*, mpi_datatypes* );
#ifdef __cplusplus