diff --git a/PG-PuReMD/src/box.c b/PG-PuReMD/src/box.c
index a223eebed2dfb7c798e4eaae3ae17cde20c1aa9a..26dc4640bdcb05330bec43a5a4aa291238179257 100644
--- a/PG-PuReMD/src/box.c
+++ b/PG-PuReMD/src/box.c
@@ -182,7 +182,7 @@ void Setup_Big_Box( real a, real b, real c, real alpha, real beta, real gamma,
box->box[2][2] = c * SQRT(1.0 - SQR(c_beta) - SQR(zi));
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "box is %8.2f x %8.2f x %8.2f\n",
+ fprintf( stderr, "[INFO] box = (%8.4f, %8.4f, %8.4f)\n",
box->box[0][0], box->box[1][1], box->box[2][2] );
#endif
@@ -378,26 +378,34 @@ void Scale_Box( reax_system * const system, control_params * const control,
}
/* temperature scaler */
- lambda = 1.0 + (dt / control->Tau_T) * (control->T / data->therm.T - 1.0);
+ lambda = 1.0 + ((dt * 1.0e-12) / control->Tau_T)
+ * (control->T / data->therm.T - 1.0);
+
if ( lambda < MIN_dT )
{
lambda = MIN_dT;
}
- else if ( lambda > MAX_dT )
+
+ lambda = SQRT( lambda );
+
+ if ( lambda > MAX_dT )
{
lambda = MAX_dT;
}
- lambda = SQRT( lambda );
/* Scale velocities and positions at t+dt */
for ( i = 0; i < system->n; ++i )
{
atom = &system->my_atoms[i];
+
rvec_Scale( atom->v, lambda, atom->v );
+
atom->x[0] = mu[0] * atom->x[0];
atom->x[1] = mu[1] * atom->x[1];
atom->x[2] = mu[2] * atom->x[2];
}
+
+ /* update kinetic energy and temperature based on new positions and velocities */
Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
/* update box & grid */
diff --git a/PG-PuReMD/src/control.c b/PG-PuReMD/src/control.c
index 7c2d478cd31ee18784e0752218b3ef31ff6d4785..718372023d7ba6ebcac3d96991532b1019dfcc28 100644
--- a/PG-PuReMD/src/control.c
+++ b/PG-PuReMD/src/control.c
@@ -91,7 +91,7 @@ void Read_Control_File( const char * const control_file, control_params * const
control->cm_domain_sparsity = 1.0;
control->cm_solver_pre_comp_type = JACOBI_PC;
control->cm_solver_pre_comp_sweeps = 3;
- control->cm_solver_pre_comp_refactor = 100;
+ control->cm_solver_pre_comp_refactor = 1;
control->cm_solver_pre_comp_droptol = 0.01;
control->cm_solver_pre_app_type = TRI_SOLVE_PA;
control->cm_solver_pre_app_jacobi_iters = 50;
@@ -378,9 +378,9 @@ void Read_Control_File( const char * const control_file, control_params * const
val = atof(tmp[1]);
control->T_init = val;
- if ( control->T_init < 0.1 )
+ if ( control->T_init < 0.001 )
{
- control->T_init = 0.1;
+ control->T_init = 0.001;
}
}
else if ( strncmp(tmp[0], "temp_final", MAX_LINE) == 0 )
@@ -396,7 +396,8 @@ void Read_Control_File( const char * const control_file, control_params * const
else if ( strncmp(tmp[0], "t_mass", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
- control->Tau_T = val * 1.e-3; // convert t_mass from fs to ps
+ /* convert from fs to s */
+ control->Tau_T = val * 1.0e-15;
}
else if ( strncmp(tmp[0], "t_mode", MAX_LINE) == 0 )
{
diff --git a/PG-PuReMD/src/cuda/cuda_integrate.cu b/PG-PuReMD/src/cuda/cuda_integrate.cu
index 9889297426f7e09498bc7d504b9b312423043c4b..8a32018a6239edbad37646f46e53f89f71afc3b3 100644
--- a/PG-PuReMD/src/cuda/cuda_integrate.cu
+++ b/PG-PuReMD/src/cuda/cuda_integrate.cu
@@ -43,6 +43,7 @@ CUDA_GLOBAL void k_velocity_verlet_part1( reax_atom *my_atoms,
/* Compute x(t + dt) */
rvec_ScaledSum( dx, dt, atom->v, -0.5 * F_CONV * inv_m * SQR(dt), atom->f );
rvec_Add( atom->x, dx );
+
/* Compute v(t + dt/2) */
rvec_ScaledAdd( atom->v, -0.5 * F_CONV * inv_m * dt, atom->f );
}
diff --git a/PG-PuReMD/src/geo_tools.c b/PG-PuReMD/src/geo_tools.c
index 3b3b0c1a2cef4f9eb947fcf5b88925e233deffbb..f504a1134db8ab32091de1307db59a1edbef2217 100644
--- a/PG-PuReMD/src/geo_tools.c
+++ b/PG-PuReMD/src/geo_tools.c
@@ -87,8 +87,10 @@ void Read_Geo_File( const char * const geo_file, reax_system * const system,
/* read box information */
fscanf( geo, CUSTOM_BOXGEO_FORMAT,
descriptor, &box_x, &box_y, &box_z, &alpha, &beta, &gamma );
+
/* initialize the box */
- Setup_Big_Box( box_x, box_y, box_z, alpha, beta, gamma, &(system->big_box) );
+ Setup_Big_Box( box_x, box_y, box_z, alpha, beta, gamma, &system->big_box );
+
/* initialize the simulation environment */
Setup_Environment( system, control, mpi_data );
@@ -122,7 +124,7 @@ void Read_Geo_File( const char * const geo_file, reax_system * const system,
{
element[j] = toupper( element[j] );
}
- atom->type = Get_Atom_Type( &system->reax_param, element );
+ atom->type = Get_Atom_Type( &system->reax_param, element, sizeof(element) );
strncpy( atom->name, name, sizeof(atom->name) - 1 );
atom->name[sizeof(atom->name) - 1] = '\0';
rvec_Copy( atom->x, x );
@@ -160,7 +162,7 @@ static int Read_Box_Info( reax_system * const system, FILE *geo, int geo_format
char s_group[12], s_zValue[12];
int ret;
- ret = SUCCESS;
+ ret = FAILURE;
fseek( geo, 0, SEEK_SET );
@@ -169,8 +171,12 @@ static int Read_Box_Info( reax_system * const system, FILE *geo, int geo_format
case PDB:
cryst = "CRYST1";
break;
+ case BGF:
+ cryst = "CRYSTX";
+ break;
default:
cryst = "BOX";
+ break;
}
/* locate the cryst line in the geo file, read it and
@@ -181,18 +187,23 @@ static int Read_Box_Info( reax_system * const system, FILE *geo, int geo_format
{
if ( geo_format == PDB )
{
- sscanf( line, PDB_CRYST1_FORMAT,
- &descriptor[0],
- &s_a[0], &s_b[0], &s_c[0],
- &s_alpha[0], &s_beta[0], &s_gamma[0],
- &s_group[0], &s_zValue[0] );
+ sscanf( line, PDB_CRYST1_FORMAT, descriptor,
+ s_a, s_b, s_c, s_alpha, s_beta, s_gamma,
+ s_group, s_zValue );
+ }
+ else if ( geo_format == BGF )
+ {
+ sscanf( line, BGF_CRYSTX_FORMAT, descriptor,
+ s_a, s_b, s_c, s_alpha, s_beta, s_gamma );
}
/* compute full volume tensor from the angles */
- Setup_Big_Box( atof(s_a), atof(s_b), atof(s_c),
+ Setup_Big_Box( atof(s_a), atof(s_b), atof(s_c),
atof(s_alpha), atof(s_beta), atof(s_gamma),
&system->big_box );
+ ret = SUCCESS;
+
break;
}
}
@@ -200,6 +211,10 @@ static int Read_Box_Info( reax_system * const system, FILE *geo, int geo_format
{
ret = FAILURE;
}
+ else
+ {
+ fseek( geo, 0, SEEK_SET );
+ }
return ret;
}
@@ -274,7 +289,8 @@ void Read_PDB_File( const char * const pdb_file, reax_system * const system,
pdb = sfopen( pdb_file, "r", "Read_PDB_File::pdb" );
- Allocate_Tokenizer_Space( &s, &s1, &tmp );
+ Allocate_Tokenizer_Space( &s, MAX_LINE, &s1, MAX_LINE,
+ &tmp, MAX_TOKENS, MAX_TOKEN_LEN );
if ( Read_Box_Info( system, pdb, PDB ) == FAILURE )
{
@@ -283,7 +299,9 @@ void Read_PDB_File( const char * const pdb_file, reax_system * const system,
MPI_Abort( MPI_COMM_WORLD, INVALID_GEO );
}
+ /* initialize the simulation environment */
Setup_Environment( system, control, mpi_data );
+
Count_PDB_Atoms( pdb, system );
PreAllocate_Space( system, control, workspace );
@@ -311,70 +329,71 @@ void Read_PDB_File( const char * const pdb_file, reax_system * const system,
{
if ( strncmp(tmp[0], "ATOM", 4) == 0 )
{
- strncpy( &descriptor[0], s1, 6 );
- descriptor[6] = 0;
- strncpy( &serial[0], &s1[6], 5 );
- serial[5] = 0;
- strncpy( &atom_name[0], &s1[12], 4 );
- atom_name[4] = 0;
- //strncpy( &serial[0], &s1[6], 7 ); serial[7] = 0;
- //strncpy( &atom_name[0], &s1[13], 3 ); atom_name[3] = 0;
+ strncpy( descriptor, s1, 6 );
+ descriptor[6] = '\0';
+ strncpy( serial, &s1[6], 5 );
+ serial[5] = '\0';
+ strncpy( atom_name, &s1[12], 4 );
+ atom_name[4] = '\0';
+ //strncpy( serial, &s1[6], 7 ); serial[7] = 0;
+ //strncpy( atom_name, &s1[13], 3 ); atom_name[3] = 0;
alt_loc = s1[16];
- strncpy( &res_name[0], &s1[17], 3 );
- res_name[3] = 0;
+ strncpy( res_name, &s1[17], 3 );
+ res_name[3] = '\0';
chain_id = s1[21];
- strncpy( &res_seq[0], &s1[22], 4 );
- res_seq[4] = 0;
+ strncpy( res_seq, &s1[22], 4 );
+ res_seq[4] = '\0';
icode = s1[26];
- strncpy( &s_x[0], &s1[30], 8 );
- s_x[8] = 0;
- strncpy( &s_y[0], &s1[38], 8 );
- s_y[8] = 0;
- strncpy( &s_z[0], &s1[46], 8 );
- s_z[8] = 0;
- strncpy( &occupancy[0], &s1[54], 6 );
- occupancy[6] = 0;
- strncpy( &temp_factor[0], &s1[60], 6 );
- temp_factor[6] = 0;
- strncpy( &seg_id[0], &s1[72], 4 );
- seg_id[4] = 0;
- strncpy( &element[0], &s1[76], 2 );
- element[2] = 0;
- strncpy( &charge[0], &s1[78], 2 );
- charge[2] = 0;
+ strncpy( s_x, &s1[30], 8 );
+ s_x[8] = '\0';
+ strncpy( s_y, &s1[38], 8 );
+ s_y[8] = '\0';
+ strncpy( s_z, &s1[46], 8 );
+ s_z[8] = '\0';
+ strncpy( occupancy, &s1[54], 6 );
+ occupancy[6] = '\0';
+ strncpy( temp_factor, &s1[60], 6 );
+ temp_factor[6] = '\0';
+ strncpy( seg_id, &s1[72], 4 );
+ seg_id[4] = '\0';
+ strncpy( element, &s1[76], 2 );
+ element[2] = '\0';
+ strncpy( charge, &s1[78], 2 );
+ charge[2] = '\0';
}
else if (strncmp(tmp[0], "HETATM", 6) == 0)
{
- strncpy( &descriptor[0], s1, 6 );
- descriptor[6] = 0;
- strncpy( &serial[0], &s1[6], 5 );
- serial[5] = 0;
- strncpy( &atom_name[0], &s1[12], 4 );
- atom_name[4] = 0;
- //strncpy( &serial[0], &s1[6], 7 ); serial[7] = 0;
- //strncpy( &atom_name[0], &s1[13], 3 ); atom_name[3] = 0;
+ strncpy( descriptor, s1, 6 );
+ descriptor[6] = '\0';
+ strncpy( serial, &s1[6], 5 );
+ serial[5] = '\0';
+ strncpy( atom_name, &s1[12], 4 );
+ atom_name[4] = '\0';
+ //strncpy( serial, &s1[6], 7 ); serial[7] = 0;
+ //strncpy( atom_name, &s1[13], 3 ); atom_name[3] = 0;
alt_loc = s1[16];
- strncpy( &res_name[0], &s1[17], 3 );
- res_name[3] = 0;
+ strncpy( res_name, &s1[17], 3 );
+ res_name[3] = '\0';
chain_id = s1[21];
- strncpy( &res_seq[0], &s1[22], 4 );
- res_seq[4] = 0;
+ strncpy( res_seq, &s1[22], 4 );
+ res_seq[4] = '\0';
icode = s1[26];
- strncpy( &s_x[0], &s1[30], 8 );
- s_x[8] = 0;
- strncpy( &s_y[0], &s1[38], 8 );
- s_y[8] = 0;
- strncpy( &s_z[0], &s1[46], 8 );
- s_z[8] = 0;
- strncpy( &occupancy[0], &s1[54], 6 );
- occupancy[6] = 0;
- strncpy( &temp_factor[0], &s1[60], 6 );
- temp_factor[6] = 0;
- //strncpy( &seg_id[0], &s1[72], 4 ); seg_id[4] = 0;
- strncpy( &element[0], &s1[76], 2 );
- element[2] = 0;
- strncpy( &charge[0], &s1[78], 2 );
- charge[2] = 0;
+ strncpy( s_x, &s1[30], 8 );
+ s_x[8] = '\0';
+ strncpy( s_y, &s1[38], 8 );
+ s_y[8] = '\0';
+ strncpy( s_z, &s1[46], 8 );
+ s_z[8] = '\0';
+ strncpy( occupancy, &s1[54], 6 );
+ occupancy[6] = '\0';
+ strncpy( temp_factor, &s1[60], 6 );
+ temp_factor[6] = '\0';
+ strncpy( seg_id, &s1[72], 4 );
+ seg_id[4] = '\0';
+ strncpy( element, &s1[76], 2 );
+ element[2] = '\0';
+ strncpy( charge, &s1[78], 2 );
+ charge[2] = '\0';
}
/* if the point is inside my_box, add it to my lists */
@@ -391,8 +410,8 @@ void Read_PDB_File( const char * const pdb_file, reax_system * const system,
pdb_serial = (int) strtod( &serial[0], &endptr );
atom->orig_id = pdb_serial;
- Trim_Spaces( element );
- atom->type = Get_Atom_Type( &system->reax_param, element );
+ Trim_Spaces( element, sizeof(element) );
+ atom->type = Get_Atom_Type( &system->reax_param, element, sizeof(element) );
strncpy( atom->name, atom_name, sizeof(atom->name) - 1 );
atom->name[sizeof(atom->name) - 1] = '\0';
@@ -462,7 +481,7 @@ void Read_PDB_File( const char * const pdb_file, reax_system * const system,
MPI_Abort( MPI_COMM_WORLD, INVALID_GEO );
}
- Deallocate_Tokenizer_Space( s, s1, tmp );
+ Deallocate_Tokenizer_Space( &s, &s1, &tmp, MAX_TOKENS );
sfclose( pdb, "Read_PDB_File::pdb" );
}
@@ -545,7 +564,7 @@ void Write_PDB_File( reax_system * const system, reax_list * const bond_list,
strncpy( name, p_atom->name, sizeof(name) - 1 );
name[sizeof(name) - 1] = '\0';
- Trim_Spaces( name );
+ Trim_Spaces( name, sizeof(name) );
snprintf( line, sizeof(line) - 1, PDB_ATOM_FORMAT_O,
"ATOM ", p_atom->orig_id, p_atom->name, ' ', "REX", ' ', 1, ' ',
@@ -611,3 +630,215 @@ void Write_PDB_File( reax_system * const system, reax_list * const bond_list,
sfree( buffer, "Write_PDB_File::buffer" );
sfree( line, "Write_PDB_File::line" );
}
+
+
+void Read_BGF( const char * const bgf_file, reax_system * const system,
+ control_params * const control, simulation_data * const data,
+ storage * const workspace, mpi_datatypes * const mpi_data )
+{
+ FILE *bgf;
+ char **tokens;
+ char *line, *backup;
+ char descriptor[7], serial[6];
+ char atom_name[6], res_name[4], res_seq[6];
+ char s_x[11], s_y[11], s_z[11];
+ char occupancy[4], temp_factor[3];
+ char element[6], charge[9];
+ char chain_id;
+ char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
+ char *endptr;
+ int i, n, atom_cnt, token_cnt, bgf_serial, ratom;
+ rvec x;
+
+ endptr = NULL;
+ ratom = 0;
+
+ bgf = sfopen( bgf_file, "r", "Read_BGF::bgf" );
+
+ Allocate_Tokenizer_Space( &line, MAX_LINE, &backup, MAX_LINE,
+ &tokens, MAX_TOKENS, MAX_TOKEN_LEN );
+
+ /* count number of atoms in the BGF file */
+ n = 0;
+ line[0] = '\0';
+
+ if ( Read_Box_Info( system, bgf, BGF ) == FAILURE )
+ {
+ fprintf( stderr, "[ERROR] Read_Box_Info: no CRYSTX line in the BGF file!" );
+ fprintf( stderr, " Terminating...\n" );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_GEO );
+ }
+
+ /* initialize the simulation environment */
+ Setup_Environment( system, control, mpi_data );
+
+ while ( fgets( line, MAX_LINE, bgf ) )
+ {
+ tokens[0][0] = '\0';
+ token_cnt = Tokenize( line, &tokens, MAX_TOKEN_LEN );
+
+ if ( strncmp( tokens[0], "ATOM", 4 ) == 0
+ || strncmp( tokens[0], "HETATM", 6 ) == 0 )
+ {
+ ++n;
+ }
+
+ line[0] = '\0';
+ }
+ if ( ferror( bgf ) )
+ {
+ fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_GEO );
+ }
+
+ sfclose( bgf, "Read_BGF::bgf" );
+
+ system->n = n;
+ system->N = n;
+ system->bigN = n;
+ PreAllocate_Space( system, control, workspace );
+
+// workspace->map_serials = scalloc( MAX_ATOM_ID, sizeof(int),
+// "Read_BGF::workspace->map_serials" );
+// for ( i = 0; i < MAX_ATOM_ID; ++i )
+// {
+// workspace->map_serials[i] = -1;
+// }
+
+ bgf = sfopen( bgf_file, "r", "Read_BGF::bgf" );
+ atom_cnt = 0;
+ token_cnt = 0;
+
+ while ( fgets( line, MAX_LINE, bgf ) )
+ {
+ /* read new line and tokenize it */
+ strncpy( backup, line, MAX_LINE - 1 );
+ backup[MAX_LINE - 1] = '\0';
+ token_cnt = Tokenize( line, &tokens, MAX_TOKEN_LEN );
+
+ /* process lines with atom info (i.e., begin with keywords HETATM or ATOM),
+ * format: "%6s %5d %-5s %3s %c %5s%10.5f%10.5f%10.5f %-5s%3d%2d %8.5f"
+ *
+ * also, it's common to see the atom info format
+ * inlined in MSI BGF files as follows:
+ * FORMAT ATOM (a6,1x,i5,1x,a5,1x,a3,1x,a1,1x,a5,3f10.5,1x,a5,i3,i2,1x,f8.5)
+ * */
+ if ( strncmp( tokens[0], "ATOM", 4 ) == 0
+ || strncmp( tokens[0], "HETATM", 6 ) == 0 )
+ {
+ strncpy( descriptor, backup, sizeof(descriptor) - 1 );
+ descriptor[sizeof(descriptor) - 1] = '\0';
+ strncpy( serial, backup + 7, sizeof(serial) - 1 );
+ serial[sizeof(serial) - 1] = '\0';
+ strncpy( atom_name, backup + 13, sizeof(atom_name) - 1 );
+ atom_name[sizeof(atom_name) - 1] = '\0';
+ strncpy( res_name, backup + 19, sizeof(res_name) - 1 );
+ res_name[sizeof(res_name) - 1] = '\0';
+ chain_id = backup[23];
+ strncpy( res_seq, backup + 25, sizeof(res_seq) - 1 );
+ res_seq[sizeof(res_seq) - 1] = '\0';
+ strncpy( s_x, backup + 30, sizeof(s_x) - 1 );
+ s_x[sizeof(s_x) - 1] = '\0';
+ strncpy( s_y, backup + 40, sizeof(s_y) - 1 );
+ s_y[sizeof(s_y) - 1] = '\0';
+ strncpy( s_z, backup + 50, sizeof(s_x) - 1 );
+ s_z[sizeof(s_x) - 1] = '\0';
+ strncpy( element, backup + 61, sizeof(element) - 1 );
+ element[sizeof(element) - 1] = '\0';
+ strncpy( occupancy, backup + 66, sizeof(occupancy) - 1 );
+ occupancy[sizeof(occupancy) - 1] = '\0';
+ strncpy( temp_factor, backup + 69, sizeof(temp_factor) - 1 );
+ temp_factor[sizeof(temp_factor) - 1] = '\0';
+ strncpy( charge, backup + 72, sizeof(charge) - 1 );
+ charge[sizeof(charge) - 1] = '\0';
+
+// sscanf( backup, BGF_ATOM_FORMAT, descriptor, serial, atom_name,
+// res_name, &chain_id, res_seq, s_x, s_y, s_z, element,
+// occupancy, temp_factor, charge );
+
+ Make_Point( strtod( s_x, &endptr ), strtod( s_y, &endptr ),
+ strtod( s_z, &endptr ), &x );
+
+ Fit_to_Periodic_Box( &system->big_box, &x );
+
+ /* if the point is inside my_box, add it to my lists */
+ if ( is_Inside_Box( &system->my_box, x ) == TRUE )
+ {
+ /* add to mapping */
+ bgf_serial = strtod( serial, &endptr );
+ Check_Input_Range( bgf_serial, 0, MAX_ATOM_ID, "Invalid bgf serial" );
+// workspace->map_serials[bgf_serial] = atom_cnt;
+ system->my_atoms[atom_cnt].orig_id = bgf_serial;
+
+ /* copy atomic positions */
+ system->my_atoms[atom_cnt].x[0] = x[0];
+ system->my_atoms[atom_cnt].x[1] = x[1];
+ system->my_atoms[atom_cnt].x[2] = x[2];
+
+ /* atom name and type */
+ strncpy( system->my_atoms[atom_cnt].name, atom_name,
+ sizeof(system->my_atoms[atom_cnt].name) - 1 );
+ system->my_atoms[atom_cnt].name[sizeof(system->my_atoms[atom_cnt].name) - 1] = '\0';
+ Trim_Spaces( element, sizeof(element) );
+ system->my_atoms[atom_cnt].type =
+ Get_Atom_Type( &system->reax_param, element, sizeof(element) );
+
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr,
+ "[INFO] atom_cnt = %5d, atom_type = %3d, x = (%10.5f,%10.5f,%10.5f), q = %10.5f, occ = %s, temp = %s, res_name = %4s, element = %s\n",
+ atom_cnt, system->my_atoms[atom_cnt].type,
+ system->my_atoms[atom_cnt].x[0],
+ system->my_atoms[atom_cnt].x[1],
+ system->my_atoms[atom_cnt].x[2],
+ system->my_atoms[atom_cnt].q, occupancy, temp_factor,
+ res_name, element );
+#endif
+
+ atom_cnt++;
+ }
+ }
+ else if ( strncmp( tokens[0], "CONECT", 6 ) == 0 )
+ {
+// if ( control->restrict_bonds )
+// {
+// /* check number of restrictions */
+// Check_Input_Range( token_cnt - 2, 0, MAX_RESTRICT,
+// "CONECT line exceeds max restrictions allowed.\n" );
+//
+// /* read bond restrictions */
+// bgf_serial = atoi(tokens[1]);
+// if ( is_Valid_Serial( workspace, bgf_serial ) )
+// {
+// ratom = workspace->map_serials[ bgf_serial ];
+// }
+//
+// workspace->restricted[ ratom ] = token_cnt - 2;
+// for ( i = 2; i < token_cnt; ++i )
+// {
+// if ( is_Valid_Serial( workspace, bgf_serial = atoi(tokens[i]) ) )
+// {
+// workspace->restricted_list[ratom][i - 2] =
+// workspace->map_serials[bgf_serial];
+// }
+// }
+// }
+ }
+
+ /* clear previous input line */
+ line[0] = '\0';
+
+ for ( i = 0; i < token_cnt; ++i )
+ {
+ tokens[i][0] = '\0';
+ }
+ }
+ if ( ferror( bgf ) )
+ {
+ fprintf( stderr, "[ERROR] Unable to read BGF file. Terminating...\n" );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_GEO );
+ }
+
+ Deallocate_Tokenizer_Space( &line, &backup, &tokens, MAX_TOKENS );
+
+ sfclose( bgf, "Read_BGF::bgf" );
+}
diff --git a/PG-PuReMD/src/geo_tools.h b/PG-PuReMD/src/geo_tools.h
index d6b332d3ba3c3f40eb9709969fa8e562660f92bc..3369d64c48a2ac5c291295b34f4cf72d805e266b 100644
--- a/PG-PuReMD/src/geo_tools.h
+++ b/PG-PuReMD/src/geo_tools.h
@@ -110,6 +110,11 @@ COLUMNS DATA TYPE FIELD DEFINITION
#define PDB_ATOM_FORMAT_O_LENGTH 81
#define PDB_CRYST1_FORMAT_O "%6s%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f%11s%4d\n"
+#define BGF_ATOM_FORMAT "%6s %5s %5s %3s %c %5s%10s%10s%10s %5s%3s%2s %8s"
+#define BGF_CRYSTX_FORMAT "%8s%11s%11s%11s%11s%11s%11s"
+
+#define BGF_ATOM_FORMAT_O "%6s %5d %-5s %3s %c %5s%10.5f%10.5f%10.5f %-5s%3d%2d %8.5f"
+
#ifdef __cplusplus
extern "C" {
@@ -121,6 +126,9 @@ void Read_Geo_File( const char * const, reax_system * const, control_params * co
void Read_PDB_File( const char * const, reax_system * const, control_params * const,
simulation_data * const, storage * const, mpi_datatypes * const );
+void Read_BGF( const char * const, reax_system * const, control_params * const,
+ simulation_data * const, storage * const, mpi_datatypes * const );
+
void Write_PDB_File( reax_system * const, reax_list * const, simulation_data * const,
control_params * const, mpi_datatypes * const, output_controls * const );
diff --git a/PG-PuReMD/src/init_md.c b/PG-PuReMD/src/init_md.c
index 0aa321f416628674341c2be894c1466575c151c8..6ccc3bf948b0d44f02736179e1c824f944447fd8 100644
--- a/PG-PuReMD/src/init_md.c
+++ b/PG-PuReMD/src/init_md.c
@@ -100,13 +100,30 @@ static void Reposition_Atoms( reax_system * const system, control_params * const
-void Generate_Initial_Velocities( reax_system * const system, real T )
+void Generate_Initial_Velocities( reax_system * const system,
+ control_params * const control, real T )
{
int i;
real m, scale, norm;
- if ( T <= 0.1 )
+ if ( T <= 0.1 || control->random_vel == FALSE )
{
+ /* warnings if conflicts between initial temperature and control file parameter */
+ if ( control->random_vel == TRUE )
+ {
+ fprintf( stderr, "[ERROR] conflicting control file parameters\n" );
+ fprintf( stderr, "[INFO] random_vel = 1 and small initial temperature (t_init = %f)\n", T );
+ fprintf( stderr, "[INFO] set random_vel = 0 to resolve this (atom initial velocites set to zero)\n" );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
+ }
+ else if ( T > 0.1 )
+ {
+ fprintf( stderr, "[ERROR] conflicting control file paramters\n" );
+ fprintf( stderr, "[INFO] random_vel = 0 and large initial temperature (t_init = %f)\n", T );
+ fprintf( stderr, "[INFO] set random_vel = 1 to resolve this (random atom initial velocites according to t_init)\n" );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
+ }
+
for ( i = 0; i < system->n; i++ )
{
rvec_MakeZero( system->my_atoms[i].v );
@@ -114,6 +131,13 @@ void Generate_Initial_Velocities( reax_system * const system, real T )
}
else
{
+ if ( T <= 0.0 )
+ {
+ fprintf( stderr, "[ERROR] random atom initial velocities specified with invalid temperature (%f). Terminating...\n",
+ T );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
+ }
+
Randomize( );
for ( i = 0; i < system->n; i++ )
@@ -211,9 +235,10 @@ void Init_System( reax_system * const system, control_params * const control,
// Reposition_Atoms( system, control, data, mpi_data );
/* initialize velocities so that desired init T can be attained */
- if ( !control->restart || (control->restart && control->random_vel) )
+ if ( control->restart == FALSE
+ || (control->restart == TRUE && control->random_vel == TRUE) )
{
- Generate_Initial_Velocities( system, control->T_init );
+ Generate_Initial_Velocities( system, control, control->T_init );
}
Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
@@ -242,7 +267,7 @@ void Init_Simulation_Data( reax_system * const system, control_params * const co
break;
case bNVT:
- data->N_f = 3 * system->bigN + 1;
+ data->N_f = 3 * system->bigN;
control->Evolve = &Velocity_Verlet_Berendsen_NVT;
control->virial = 0;
break;
diff --git a/PG-PuReMD/src/init_md.h b/PG-PuReMD/src/init_md.h
index 805a1e8c398e5431207a62142ee1a825e1efc51f..f4fa2a323f21b76303c3fb995edb58e352263ea7 100644
--- a/PG-PuReMD/src/init_md.h
+++ b/PG-PuReMD/src/init_md.h
@@ -29,7 +29,7 @@
extern "C" {
#endif
-void Generate_Initial_Velocities( reax_system * const, real );
+void Generate_Initial_Velocities( reax_system * const, control_params * const, real );
void Init_MPI_Datatypes( reax_system * const, storage * const, mpi_datatypes * const );
diff --git a/PG-PuReMD/src/integrate.c b/PG-PuReMD/src/integrate.c
index a230d74ad80fb2fa62a1a687d78447004beb82c3..704c437afff82d73633c0b96533829efcdb22e77 100644
--- a/PG-PuReMD/src/integrate.c
+++ b/PG-PuReMD/src/integrate.c
@@ -34,35 +34,37 @@
#include "vector.h"
+/* Velocity Verlet integrator for microcanonical ensemble. */
int Velocity_Verlet_NVE( reax_system * const system, control_params * const control,
simulation_data * const data, storage * const workspace, reax_list ** const lists,
output_controls * const out_control, mpi_datatypes * const mpi_data )
{
int i, steps, renbr, ret;
static int verlet_part1_done = FALSE, gen_nbr_list = FALSE;
- real inv_m, dt, dt_sqr;
+ real inv_m, scalar1, scalar2;
rvec dx;
reax_atom *atom;
ret = SUCCESS;
- dt = control->dt;
steps = data->step - data->prev_steps;
renbr = steps % control->reneighbor == 0 ? TRUE : FALSE;
+ scalar1 = -0.5 * control->dt * F_CONV;
+ scalar2 = -0.5 * SQR( control->dt ) * F_CONV;
if ( verlet_part1_done == FALSE )
{
- dt_sqr = SQR(dt);
-
/* velocity verlet, 1st part */
for ( i = 0; i < system->n; i++ )
{
atom = &system->my_atoms[i];
inv_m = 1.0 / system->reax_param.sbp[atom->type].mass;
+
/* Compute x(t + dt) */
- rvec_ScaledSum( dx, dt, atom->v, -0.5 * dt_sqr * F_CONV * inv_m, atom->f );
+ rvec_ScaledSum( dx, control->dt, atom->v, scalar2 * inv_m, atom->f );
rvec_Add( system->my_atoms[i].x, dx );
+
/* Compute v(t + dt/2) */
- rvec_ScaledAdd( atom->v, -0.5 * dt * F_CONV * inv_m, atom->f );
+ rvec_ScaledAdd( atom->v, scalar1 * inv_m, atom->f );
}
Reallocate_Part1( system, control, data, workspace, lists, mpi_data );
@@ -102,7 +104,9 @@ int Velocity_Verlet_NVE( reax_system * const system, control_params * const cont
{
atom = &system->my_atoms[i];
inv_m = 1.0 / system->reax_param.sbp[ atom->type ].mass;
- rvec_ScaledAdd( atom->v, -0.5 * dt * F_CONV * inv_m, atom->f );
+
+ /* Compute v(t + dt) */
+ rvec_ScaledAdd( atom->v, scalar1 * inv_m, atom->f );
}
verlet_part1_done = FALSE;
@@ -113,52 +117,57 @@ int Velocity_Verlet_NVE( reax_system * const system, control_params * const cont
}
-int Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system * const system,
- control_params * const control, simulation_data * const data, storage * const workspace,
- reax_list ** const lists, output_controls * const out_control, mpi_datatypes * const mpi_data )
+/* Velocity Verlet integrator for constant volume and temperature
+ * with Berendsen thermostat.
+ *
+ * NOTE: All box dimensions are scaled by the same amount, and
+ * there is no change in the angles between axes. */
+int Velocity_Verlet_Berendsen_NVT( reax_system * const system, control_params * const control,
+ simulation_data * const data, storage * const workspace, reax_list ** const lists,
+ output_controls * const out_control, mpi_datatypes * const mpi_data )
{
- int i, itr, steps, renbr, ret, ret_mpi;
+ int i, steps, renbr, ret;
static int verlet_part1_done = FALSE, gen_nbr_list = FALSE;
- real inv_m, coef_v;
- real dt, dt_sqr;
- real my_ekin, new_ekin;
- real G_xi_new, v_xi_new, v_xi_old;
+ real inv_m, scalar1, scalar2, lambda;
rvec dx;
- thermostat *therm;
reax_atom *atom;
ret = SUCCESS;
- dt = control->dt;
- therm = &data->therm;
steps = data->step - data->prev_steps;
renbr = steps % control->reneighbor == 0 ? TRUE : FALSE;
+ scalar1 = -0.5 * control->dt * F_CONV;
+ scalar2 = -0.5 * SQR( control->dt ) * F_CONV;
if ( verlet_part1_done == FALSE )
{
- dt_sqr = SQR(dt);
-
/* velocity verlet, 1st part */
for ( i = 0; i < system->n; i++ )
{
atom = &system->my_atoms[i];
inv_m = 1.0 / system->reax_param.sbp[atom->type].mass;
- rvec_ScaledSum( dx, dt, atom->v, -0.5 * dt_sqr * F_CONV * inv_m, atom->f );
- rvec_Add( system->my_atoms[i].x, dx );
- rvec_Copy( atom->f_old, atom->f );
+
+ /* Compute x(t + dt) */
+ rvec_ScaledSum( dx, control->dt, atom->v, scalar2 * inv_m, atom->f );
+ rvec_Add( atom->x, dx );
+
+ /* Compute v(t + dt/2) */
+ rvec_ScaledAdd( atom->v, scalar1 * inv_m, atom->f );
}
-
- /* Compute xi(t + dt) */
- therm->xi += therm->v_xi * dt + 0.5 * dt_sqr * therm->G_xi;
Reallocate_Part1( system, control, data, workspace, lists, mpi_data );
+ if ( renbr == TRUE )
+ {
+ Update_Grid( system, control, MPI_COMM_WORLD );
+ }
+
Comm_Atoms( system, control, data, workspace, mpi_data, renbr );
verlet_part1_done = TRUE;
}
Reallocate_Part2( system, control, data, workspace, lists, mpi_data );
-
+
Reset( system, control, data, workspace, lists );
if ( renbr == TRUE && gen_nbr_list == FALSE )
@@ -183,47 +192,44 @@ int Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system * const system,
if ( ret == SUCCESS )
{
- /* Compute iteration constants for each atom's velocity */
- for ( i = 0; i < system->n; ++i )
+ /* velocity verlet, 2nd part */
+ for ( i = 0; i < system->n; i++ )
{
atom = &system->my_atoms[i];
inv_m = 1.0 / system->reax_param.sbp[atom->type].mass;
- rvec_Scale( workspace->v_const[i], 1.0 - 0.5 * dt * therm->v_xi, atom->v );
- rvec_ScaledAdd( workspace->v_const[i], -0.5 * dt * inv_m * F_CONV, atom->f_old );
- rvec_ScaledAdd( workspace->v_const[i], -0.5 * dt * inv_m * F_CONV, atom->f );
+
+ /* Compute v(t + dt) */
+ rvec_ScaledAdd( atom->v, scalar1 * inv_m, atom->f );
}
-
- v_xi_new = therm->v_xi_old + 2.0 * dt * therm->G_xi;
- my_ekin = G_xi_new = v_xi_old = 0;
- itr = 0;
- do
+
+ Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+
+ /* temperature scaler */
+ lambda = 1.0 + ((control->dt * 1.0e-12) / control->Tau_T)
+ * (control->T / data->therm.T - 1.0);
+
+ if ( lambda < MIN_dT )
{
- itr++;
-
- /* new values become old in this iteration */
- v_xi_old = v_xi_new;
-
- my_ekin = 0;
- for ( i = 0; i < system->n; ++i )
- {
- atom = &system->my_atoms[i];
- coef_v = 1.0 / (1.0 + 0.5 * dt * v_xi_old);
- rvec_Scale( atom->v, coef_v, workspace->v_const[i] );
- my_ekin += 0.5 * system->reax_param.sbp[atom->type].mass
- * rvec_Dot( atom->v, atom->v );
- }
-
- ret_mpi = MPI_Allreduce( &my_ekin, &new_ekin, 1, MPI_DOUBLE,
- MPI_SUM, mpi_data->comm_mesh3D );
- Check_MPI_Error( ret_mpi, __FILE__, __LINE__ );
-
- G_xi_new = control->Tau_T * ( 2.0 * new_ekin - data->N_f * K_B * control->T );
- v_xi_new = therm->v_xi + 0.5 * dt * ( therm->G_xi + G_xi_new );
- } while ( FABS(v_xi_new - v_xi_old) > 1e-5 );
+ lambda = MIN_dT;
+ }
- therm->v_xi_old = therm->v_xi;
- therm->v_xi = v_xi_new;
- therm->G_xi = G_xi_new;
+ lambda = SQRT( lambda );
+
+ if ( lambda > MAX_dT )
+ {
+ lambda = MAX_dT;
+ }
+
+ /* Scale velocities and positions at t+dt */
+ for ( i = 0; i < system->n; ++i )
+ {
+ atom = &system->my_atoms[i];
+
+ rvec_Scale( atom->v, lambda, atom->v );
+ }
+
+ /* update kinetic energy and temperature based on new velocities */
+ Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
verlet_part1_done = FALSE;
gen_nbr_list = FALSE;
@@ -233,54 +239,57 @@ int Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system * const system,
}
-/* uses Berendsen-type coupling for both T and P.
- * All box dimensions are scaled by the same amount,
- * there is no change in the angles between axes. */
-int Velocity_Verlet_Berendsen_NVT( reax_system * const system, control_params * const control,
- simulation_data * const data, storage * const workspace, reax_list ** const lists,
- output_controls * const out_control, mpi_datatypes * const mpi_data )
+/* Velocity Verlet integrator for constant volume and constant temperature
+ * with Nose-Hoover thermostat.
+ *
+ * Reference: Understanding Molecular Simulation, Frenkel and Smit
+ * Academic Press Inc. San Diego, 1996 p. 388-391 */
+int Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system * const system,
+ control_params * const control, simulation_data * const data, storage * const workspace,
+ reax_list ** const lists, output_controls * const out_control, mpi_datatypes * const mpi_data )
{
- int i, steps, renbr, ret;
+ int i, itr, steps, renbr, ret, ret_mpi;
static int verlet_part1_done = FALSE, gen_nbr_list = FALSE;
- real inv_m, dt, dt_sqr, lambda;
+ real inv_m, coef_v;
+ real dt, dt_sqr;
+ real my_ekin, new_ekin;
+ real G_xi_new, v_xi_new, v_xi_old;
rvec dx;
+ thermostat *therm;
reax_atom *atom;
ret = SUCCESS;
dt = control->dt;
+ therm = &data->therm;
steps = data->step - data->prev_steps;
renbr = steps % control->reneighbor == 0 ? TRUE : FALSE;
if ( verlet_part1_done == FALSE )
{
- dt_sqr = SQR( dt );
+ dt_sqr = SQR(dt);
/* velocity verlet, 1st part */
for ( i = 0; i < system->n; i++ )
{
atom = &system->my_atoms[i];
inv_m = 1.0 / system->reax_param.sbp[atom->type].mass;
- /* Compute x(t + dt) */
- rvec_ScaledSum( dx, dt, atom->v, -0.5 * F_CONV * inv_m * dt_sqr, atom->f );
- rvec_Add( atom->x, dx );
- /* Compute v(t + dt/2) */
- rvec_ScaledAdd( atom->v, -0.5 * F_CONV * inv_m * dt, atom->f );
+ rvec_ScaledSum( dx, dt, atom->v, -0.5 * dt_sqr * F_CONV * inv_m, atom->f );
+ rvec_Add( system->my_atoms[i].x, dx );
+ rvec_Copy( atom->f_old, atom->f );
}
+
+ /* Compute xi(t + dt) */
+ therm->xi += therm->v_xi * dt + 0.5 * dt_sqr * therm->G_xi;
Reallocate_Part1( system, control, data, workspace, lists, mpi_data );
- if ( renbr == TRUE )
- {
- Update_Grid( system, control, MPI_COMM_WORLD );
- }
-
Comm_Atoms( system, control, data, workspace, mpi_data, renbr );
verlet_part1_done = TRUE;
}
Reallocate_Part2( system, control, data, workspace, lists, mpi_data );
-
+
Reset( system, control, data, workspace, lists );
if ( renbr == TRUE && gen_nbr_list == FALSE )
@@ -305,39 +314,47 @@ int Velocity_Verlet_Berendsen_NVT( reax_system * const system, control_params *
if ( ret == SUCCESS )
{
- /* velocity verlet, 2nd part */
- for ( i = 0; i < system->n; i++ )
+ /* Compute iteration constants for each atom's velocity */
+ for ( i = 0; i < system->n; ++i )
{
atom = &system->my_atoms[i];
inv_m = 1.0 / system->reax_param.sbp[atom->type].mass;
- /* Compute v(t + dt) */
- rvec_ScaledAdd( atom->v, -0.5 * dt * F_CONV * inv_m, atom->f );
- }
-
- Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
-
- /* temperature scaler */
- lambda = 1.0 + (dt / control->Tau_T) * (control->T / data->therm.T - 1.0);
-
- if ( lambda < MIN_dT )
- {
- lambda = MIN_dT;
- }
- else if (lambda > MAX_dT )
- {
- lambda = MAX_dT;
+ rvec_Scale( workspace->v_const[i], 1.0 - 0.5 * dt * therm->v_xi, atom->v );
+ rvec_ScaledAdd( workspace->v_const[i], -0.5 * dt * inv_m * F_CONV, atom->f_old );
+ rvec_ScaledAdd( workspace->v_const[i], -0.5 * dt * inv_m * F_CONV, atom->f );
}
-
- lambda = SQRT( lambda );
-
- /* Scale velocities and positions at t+dt */
- for ( i = 0; i < system->n; ++i )
+
+ v_xi_new = therm->v_xi_old + 2.0 * dt * therm->G_xi;
+ my_ekin = G_xi_new = v_xi_old = 0;
+ itr = 0;
+ do
{
- atom = &system->my_atoms[i];
- rvec_Scale( atom->v, lambda, atom->v );
- }
+ itr++;
+
+ /* new values become old in this iteration */
+ v_xi_old = v_xi_new;
+
+ my_ekin = 0;
+ for ( i = 0; i < system->n; ++i )
+ {
+ atom = &system->my_atoms[i];
+ coef_v = 1.0 / (1.0 + 0.5 * dt * v_xi_old);
+ rvec_Scale( atom->v, coef_v, workspace->v_const[i] );
+ my_ekin += 0.5 * system->reax_param.sbp[atom->type].mass
+ * rvec_Dot( atom->v, atom->v );
+ }
+
+ ret_mpi = MPI_Allreduce( &my_ekin, &new_ekin, 1, MPI_DOUBLE,
+ MPI_SUM, mpi_data->comm_mesh3D );
+ Check_MPI_Error( ret_mpi, __FILE__, __LINE__ );
+
+ G_xi_new = control->Tau_T * ( 2.0 * new_ekin - data->N_f * K_B * control->T );
+ v_xi_new = therm->v_xi + 0.5 * dt * ( therm->G_xi + G_xi_new );
+ } while ( FABS(v_xi_new - v_xi_old) > 1e-5 );
- Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+ therm->v_xi_old = therm->v_xi;
+ therm->v_xi = v_xi_new;
+ therm->G_xi = G_xi_new;
verlet_part1_done = FALSE;
gen_nbr_list = FALSE;
@@ -347,8 +364,9 @@ int Velocity_Verlet_Berendsen_NVT( reax_system * const system, control_params *
}
-/* uses Berendsen-type coupling for both T and P.
- * All box dimensions are scaled by the same amount,
+/* Velocity Verlet integrator for constant pressure and constant temperature.
+ *
+ * NOTE: All box dimensions are scaled by the same amount, and
* there is no change in the angles between axes. */
int Velocity_Verlet_Berendsen_NPT( reax_system * const system, control_params * const control,
simulation_data * const data, storage * const workspace, reax_list ** const lists,
@@ -372,9 +390,11 @@ int Velocity_Verlet_Berendsen_NPT( reax_system * const system, control_params *
{
atom = &system->my_atoms[i];
inv_m = 1.0 / system->reax_param.sbp[atom->type].mass;
+
/* Compute x(t + dt) */
rvec_ScaledSum( dx, dt, atom->v, -0.5 * F_CONV * inv_m * SQR(dt), atom->f );
rvec_Add( atom->x, dx );
+
/* Compute v(t + dt/2) */
rvec_ScaledAdd( atom->v, -0.5 * F_CONV * inv_m * dt, atom->f );
}
@@ -420,6 +440,7 @@ int Velocity_Verlet_Berendsen_NPT( reax_system * const system, control_params *
{
atom = &system->my_atoms[i];
inv_m = 1.0 / system->reax_param.sbp[atom->type].mass;
+
/* Compute v(t + dt) */
rvec_ScaledAdd( atom->v, -0.5 * dt * F_CONV * inv_m, atom->f );
}
diff --git a/PG-PuReMD/src/puremd.c b/PG-PuReMD/src/puremd.c
index 7ec1f5186811dfdea00a47c019915422330d1a00..0d1b8f579d046e9aeee016696e58f830106cb588 100644
--- a/PG-PuReMD/src/puremd.c
+++ b/PG-PuReMD/src/puremd.c
@@ -70,6 +70,10 @@ static void Read_Config_Files( const char * const geo_file,
{
Read_PDB_File( geo_file, system, control, data, workspace, mpi_data );
}
+ else if ( control->geo_format == BGF )
+ {
+ Read_BGF( 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 );
@@ -104,9 +108,12 @@ static void Cuda_Post_Evolve( reax_system * const system, control_params * const
Cuda_Remove_CoM_Velocities( system, control, data );
}
- /* compute kinetic energy of the system */
- Cuda_Compute_Kinetic_Energy( system, control, workspace,
- data, mpi_data->comm_mesh3D );
+ if ( control->ensemble == NVE )
+ {
+ /* compute kinetic energy of the system */
+ Cuda_Compute_Kinetic_Energy( system, control, workspace,
+ data, mpi_data->comm_mesh3D );
+ }
if ( (out_control->energy_update_freq > 0
&& (data->step - data->prev_steps) % out_control->energy_update_freq == 0)
@@ -145,8 +152,11 @@ static void Post_Evolve( reax_system * const system, control_params * const cont
}
}
- /* compute kinetic energy of system */
- Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+ if ( control->ensemble == NVE )
+ {
+ /* compute kinetic energy of system */
+ Compute_Kinetic_Energy( system, data, mpi_data->comm_mesh3D );
+ }
if ( (out_control->energy_update_freq > 0
&& (data->step - data->prev_steps) % out_control->energy_update_freq == 0)
diff --git a/PG-PuReMD/src/reax_types.h b/PG-PuReMD/src/reax_types.h
index 3aaeac92e0956134419f656658456d8748509c4c..0d1895d114529d7a0be8917228fae84e063a79b3 100644
--- a/PG-PuReMD/src/reax_types.h
+++ b/PG-PuReMD/src/reax_types.h
@@ -69,6 +69,12 @@
//#define TEST_ENERGY
/* compile force calculation test code */
//#define TEST_FORCES
+/* constants defined in reference Fortran ReaxFF code (useful for comparisons) */
+#define USE_REF_FORTRAN_REAXFF_CONSTANTS
+/* constants defined in reference Fortran eReaxFF code (useful for comparisons) */
+//#define USE_REF_FORTRAN_EREAXFF_CONSTANTS
+/* constants defined in LAMMPS ReaxFF code (useful for comparisons) */
+//#define USE_LAMMPS_REAXFF_CONSTANTS
/* pack/unpack MPI buffers on device and use CUDA-aware MPI for messaging */
//TODO: rewrite this to use configure option to include
/* OpenMPI extensions for CUDA-aware support */
@@ -121,35 +127,92 @@
#define MIN(x,y) (((x) < (y)) ? (x) : (y))
#define MAX3(x,y,z) MAX( MAX((x),(y)), (z))
-/* ??? */
-#define C_ELE (332.06371)
-/* kcal/mol/K */
-//#define K_B (503.398008)
-/* amu A^2 / ps^2 / K */
-#define K_B (0.831687)
-/* --> amu A / ps^2 */
-#define F_CONV (1.0e6 / 48.88821291 / 48.88821291)
-/* amu A^2 / ps^2 --> kcal/mol */
-#define E_CONV (0.002391)
-/* conversion factor from electron volts to kilo calories per mole */
-#define EV_to_KCALpMOL (14.400000)
-/* conversion factor from kilo calories per mode to electron volts */
-//#define KCALpMOL_to_EV (23.060549) // (23.020000)
-#define KCALpMOL_to_EV (23.020000)
-/* conversion factor from (elemental charge * angstroms) to debye */
-#define ECxA_to_DEBYE (4.803204)
-/* conversion factor from calories to joules */
-#define CAL_to_JOULES (4.184000)
-/* conversion factor from joules to calories */
-#define JOULES_to_CAL (1.0 / 4.184000)
-/* conversion factor from (unified) atomic mass units to grams */
-#define AMU_to_GRAM (1.6605e-24)
-/* conversion factor from angstroms to centimenters */
-#define ANG_to_CM (1.0e-8)
-/* Avogadro's constant */
-#define AVOGNR (6.0221367e23)
-/* ??? */
-#define P_CONV (1.0e-24 * AVOGNR * JOULES_to_CAL)
+#if defined(USE_REF_FORTRAN_REAXFF_CONSTANTS)
+ /* transcendental constant pi */
+ #define PI (3.14159265)
+ /* unit conversion from ??? to kcal / mol */
+ #define C_ELE (332.0638)
+ /* Boltzmann constant, AMU * A^2 / (ps^2 * K) */
+ #define K_B (0.831687)
+// #define K_B (0.8314510)
+ /* unit conversion for atomic force to AMU * A / ps^2 */
+ #define F_CONV (4.184e2)
+ /* energy conversion constant from electron volts to kilo-calories per mole */
+ #define KCALpMOL_to_EV (23.02)
+ /* electric dipole moment conversion constant from elementary charge * angstrom to debye */
+ #define ECxA_to_DEBYE (4.80320679913)
+#elif defined(USE_REF_FORTRAN_EREAXFF_CONSTANTS)
+ //TODO
+ /* energy conversion constant from electron volts to kilo-calories per mole */
+ #define KCALpMOL_to_EV (23.02)
+ /* electric dipole moment conversion constant from elementary charge * angstrom to debye */
+ #define ECxA_to_DEBYE (4.80320679913)
+#elif defined(USE_LAMMPS_REAXFF_CONSTANTS)
+ //TODO
+ /* energy conversion constant from electron volts to kilo-calories per mole */
+ #define KCALpMOL_to_EV (23.060549)
+#endif
+
+/* Coulomb energy conversion */
+#if !defined(C_ELE)
+ #define C_ELE (332.06371)
+#endif
+/* Boltzmann constant */
+#if !defined(K_B)
+ /* in ??? */
+// #define K_B (503.398008)
+ /* Boltzmann constant, AMU * A^2 / (ps^2 * K) */
+// #define K_B (0.831687)
+ #define K_B (0.8314510)
+#endif
+/* unit conversion for atomic force */
+#if !defined(F_CONV)
+ /* to AMU * A / ps^2 */
+ #define F_CONV (1.0e6 / 48.88821291 / 48.88821291)
+#endif
+/* unit conversion for atomic energy */
+#if !defined(E_CONV)
+ /* AMU * Angstroms^2 / ps^2 --> kcal / mol */
+ #define E_CONV (0.002391)
+#endif
+/* energy conversion constant from electron volts to kilo-calories per mole */
+#if !defined(EV_to_KCALpMOL)
+ #define EV_to_KCALpMOL (14.40)
+#endif
+/* energy conversion constant from electron volts to kilo-calories per mole */
+#if !defined(KCALpMOL_to_EV)
+ #define KCALpMOL_to_EV (23.0408)
+#endif
+/* electric dipole moment conversion constant from elementary charge * angstrom to debye */
+#if !defined(ECxA_to_DEBYE)
+ #define ECxA_to_DEBYE (4.803204)
+#endif
+/* energy conversion constant from (gram) calories to Joules (SI) */
+#if !defined(CAL_to_JOULES)
+ #define CAL_to_JOULES (4.1840)
+#endif
+/* energy conversion constant from Joules (SI) to (gram) calories */
+#if !defined(JOULES_to_CAL)
+ #define JOULES_to_CAL (1.0 / 4.1840)
+#endif
+/* mass conversion constant from unified atomic mass units (daltons) to grams */
+#if !defined(AMU_to_GRAM)
+ #define AMU_to_GRAM (1.6605e-24)
+#endif
+/* distance conversion constant from angstroms to centimeters */
+#if !defined(ANG_to_CM)
+ #define ANG_to_CM (1.0e-8)
+#endif
+/* Avogradro's constant */
+#if !defined(AVOGNR)
+ #define AVOGNR (6.0221367e23)
+#endif
+/* unit conversion for pressure:
+ * (1 s / 10^12 ps) * (1 m / 10^10 Angstroms) * (6.0221367^23 atoms / mole) * (0.2390057 J / cal)
+ * ps * Angstroms * moles * cals => s * m * atoms * J */
+#if !defined(P_CONV)
+ #define P_CONV (1.0e-24 * AVOGNR * JOULES_to_CAL)
+#endif
/* max. num. of characters for string buffers */
#define MAX_STR (1024)
@@ -162,6 +225,8 @@
/* max. num. of characters in atom name */
#define MAX_ATOM_NAME_LEN (8)
+/* max. atom ID in geo file */
+#define MAX_ATOM_ID (100000)
/* ??? */
#define MAX_RESTRICT (15)
/* max. num. atoms per molecule */
@@ -222,9 +287,9 @@
/* min. pressure scaler for simulation box dimenion in NPT ensembles */
#define MIN_dV (0.99)
/* max. temperature scaler for atomic positions and velocities in NPT ensembles */
-#define MAX_dT (4.00)
+#define MAX_dT (2.0)
/* min. temperature scaler for atomic positions and velocities in NPT ensembles */
-#define MIN_dT (0.00)
+#define MIN_dT (0.0)
/* ??? */
#define MASTER_NODE (0)
@@ -403,9 +468,10 @@ enum geo_formats
{
CUSTOM = 0,
PDB = 1,
- ASCII_RESTART = 2,
- BINARY_RESTART = 3,
- GF_N = 4,
+ BGF = 2,
+ ASCII_RESTART = 3,
+ BINARY_RESTART = 4,
+ GF_N = 5,
};
/* method for computing atomic charges */
diff --git a/PG-PuReMD/src/system_props.c b/PG-PuReMD/src/system_props.c
index f42f5a908f25a780ec412ef86b639a26ff92136d..32365607fb3993c073df3c42d47de5ae8764bb61 100644
--- a/PG-PuReMD/src/system_props.c
+++ b/PG-PuReMD/src/system_props.c
@@ -89,8 +89,9 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data,
m = system->reax_param.sbp[system->my_atoms[i].type].mass;
rvec_Scale( p, m, system->my_atoms[i].v );
- data->my_en.e_kin += 0.5 * rvec_Dot( p, system->my_atoms[i].v );
+ data->my_en.e_kin += rvec_Dot( p, system->my_atoms[i].v );
}
+ data->my_en.e_kin *= 0.5;
ret = MPI_Allreduce( &data->my_en.e_kin, &data->sys_en.e_kin, 1, MPI_DOUBLE,
MPI_SUM, comm );
diff --git a/PG-PuReMD/src/tool_box.c b/PG-PuReMD/src/tool_box.c
index ea446a6ab59e6357a0bde7c69a21bf4adc4ce936..c9f0cf8db5cceecefad74b9bcb30ec86ca971f3b 100644
--- a/PG-PuReMD/src/tool_box.c
+++ b/PG-PuReMD/src/tool_box.c
@@ -34,6 +34,9 @@
#include "reax_comm_tools.h"
#endif
+#include <ctype.h>
+#include <string.h>
+
/************** taken from comm_tools.c **************/
int SumScan( int n, int me, int root, MPI_Comm comm )
@@ -155,22 +158,31 @@ int Check_Input_Range( int val, int lo, int hi, char *message )
}
-/* Note: element must be NULL-terminated before calling */
-void Trim_Spaces( char *element )
+void Trim_Spaces( char * const element, const size_t size )
{
- int i, j;
+ int i, j, n;
+
+ element[size - 1] = '\0';
+ n = strlen( element );
+
+ /* buffer not NULL-terminated, abort */
+ if ( n == size )
+ {
+ fprintf( stderr, "[ERROR] buffer not NULL-terminated (Trim_Spaces). Terminating...\n" );
+ exit( RUNTIME_ERROR );
+ }
/* skip initial space chars */
for ( i = 0; element[i] == ' '; ++i )
;
- /* strlen safe here only if element is NULL-terminated before calling Trim_Spaces */
- for ( j = i; j < (int) strlen(element) && element[j] != ' '; ++j )
+ /* make uppercase, offset to 0 */
+ for ( j = i; j < n && element[j] != ' '; ++j )
{
- /* make uppercase, offset to 0 */
element[j - i] = toupper( element[j] );
}
- /* NULL-terminate the string */
+
+ /* NULL terminate */
element[j - i] = '\0';
}
@@ -210,16 +222,16 @@ void Update_Timing_Info( real *t_start, real *timing )
/*********** from io_tools.c **************/
-int Get_Atom_Type( reax_interaction *reax_param, char *s )
+int Get_Atom_Type( reax_interaction *reax_param, char *s, size_t n )
{
int i, ret, flag;
flag = FAILURE;
- ret = -1;
for ( i = 0; i < reax_param->num_atom_types; ++i )
{
- if ( strncmp( reax_param->sbp[i].name, s, 15 ) == 0 )
+ if ( strncmp( reax_param->sbp[i].name, s,
+ MIN( sizeof(reax_param->sbp[i].name), n ) ) == 0 )
{
ret = i;
flag = SUCCESS;
@@ -249,33 +261,37 @@ char *Get_Atom_Name( reax_system *system, int i )
}
-void Deallocate_Tokenizer_Space( char *line, char *backup, char **tokens )
+void Allocate_Tokenizer_Space( char **line, size_t line_size,
+ char **backup, size_t backup_size,
+ char ***tokens, size_t num_tokens, size_t token_size )
{
int i;
- for ( i = 0; i < MAX_TOKENS; i++ )
+ *line = smalloc( sizeof(char) * line_size, "Allocate_Tokenizer_Space::*line" );
+ *backup = smalloc( sizeof(char) * backup_size, "Allocate_Tokenizer_Space::*backup" );
+ *tokens = smalloc( sizeof(char*) * num_tokens, "Allocate_Tokenizer_Space::*tokens" );
+
+ for ( i = 0; i < num_tokens; i++ )
{
- sfree( tokens[i], "Deallocate_Tokenizer_Space::tokens[i]" );
+ (*tokens)[i] = smalloc( sizeof(char) * token_size,
+ "Allocate_Tokenizer_Space::(*tokens)[i]" );
}
-
- sfree( line, "Deallocate_Tokenizer_Space::line" );
- sfree( backup, "Deallocate_Tokenizer_Space::backup" );
- sfree( tokens, "Deallocate_Tokenizer_Space::tokens" );
}
-void Allocate_Tokenizer_Space( char **line, char **backup, char ***tokens )
+void Deallocate_Tokenizer_Space( char **line, char **backup,
+ char ***tokens, size_t num_tokens )
{
int i;
- *line = smalloc( sizeof(char) * MAX_LINE, "Allocate_Tokenizer_Space::line" );
- *backup = smalloc( sizeof(char) * MAX_LINE, "Allocate_Tokenizer_Space::backup" );
- *tokens = smalloc( sizeof(char*) * MAX_TOKENS, "Allocate_Tokenizer_Space::tokens" );
-
- for ( i = 0; i < MAX_TOKENS; i++ )
+ for ( i = 0; i < num_tokens; i++ )
{
- (*tokens)[i] = smalloc( sizeof(char) * MAX_TOKEN_LEN, "Allocate_Tokenizer_Space::tokens[i]" );
+ sfree( (*tokens)[i], "Deallocate_Tokenizer_Space::tokens[i]" );
}
+
+ sfree( *line, "Deallocate_Tokenizer_Space::line" );
+ sfree( *backup, "Deallocate_Tokenizer_Space::backup" );
+ sfree( *tokens, "Deallocate_Tokenizer_Space::tokens" );
}
diff --git a/PG-PuReMD/src/tool_box.h b/PG-PuReMD/src/tool_box.h
index 9ae34084ef43fedccdb4bfe940687fe85cf023a1..6b80cfc0b162379d477d00e439a9d829f10444c0 100644
--- a/PG-PuReMD/src/tool_box.h
+++ b/PG-PuReMD/src/tool_box.h
@@ -43,7 +43,7 @@ int is_Valid_Serial( storage*, int );
int Check_Input_Range( int, int, int, char* );
-void Trim_Spaces( char* );
+void Trim_Spaces( char * const, const size_t );
/* from system_props.h */
real Get_Time( );
@@ -53,15 +53,17 @@ real Get_Elapsed_Time( real );
void Update_Timing_Info( real*, real* );
/* from io_tools.h */
-int Get_Atom_Type( reax_interaction*, char* );
+int Get_Atom_Type( reax_interaction*, char*, size_t );
char *Get_Element( reax_system*, int );
char *Get_Atom_Name( reax_system*, int );
-void Deallocate_Tokenizer_Space( char *, char *, char ** );
+void Allocate_Tokenizer_Space( char**, size_t, char**, size_t, char***,
+ size_t, size_t );
-void Allocate_Tokenizer_Space( char **, char **, char *** );
+void Deallocate_Tokenizer_Space( char **, char **, char ***,
+ size_t );
int Tokenize( char*, char***, size_t );
diff --git a/environ/control_water b/environ/control_water
index 9cd2d3fbed13901ecf7e48d5d2b20489ad2bc6b3..1961d0bac8f120a45532596e9193440cac9ea7f7 100644
--- a/environ/control_water
+++ b/environ/control_water
@@ -1,4 +1,4 @@
-simulation_name water_6540_notab_qeq ! output files will carry this name + their specific extension
+simulation_name petn_48256_notab_qeq ! output files will carry this name + their specific extension
ensemble_type 0 ! 0: NVE, 1: Berendsen NVT, 2: nose-Hoover NVT, 3: semi-isotropic NPT, 4: isotropic NPT, 5: anisotropic NPT
nsteps 100 ! number of simulation steps
dt 0.25 ! time step in fs
@@ -44,7 +44,7 @@ p_mass 5000.00 ! in fs, pressure inertia parame
compress 0.008134 ! in ps^2 * A / amu ( 4.5X10^(-5) bar^(-1) )
press_mode 0 ! 0: internal + external pressure, 1: ext only, 2: int only
-geo_format 1 ! 0: custom, 1: pdb, 2: bgf
+geo_format 1 ! 0: custom, 1: pdb, 2: bgf 3: ASCII restart 3: binary restart
write_freq 0 ! write trajectory after so many steps
traj_compress 0 ! 0: no compression 1: uses zlib to compress trajectory output
traj_format 0 ! 0: our own format (below options apply to this only), 1: xyz, 2: bgf, 3: pdb
diff --git a/environ/parallel_control b/environ/parallel_control
index 4f6d75b2e1a4f572dc3e8dc03cc63eeb4896780c..fe851fe6000268153f07738b2a9b1d2b429040d8 100644
--- a/environ/parallel_control
+++ b/environ/parallel_control
@@ -1,6 +1,6 @@
-simulation_name water.6.notab.111 ! output files will carry this name + their specific ext.
-ensemble_type 1 ! 0: NVE, 1: Berendsen NVT, 2: Nose-Hoover NVT(under testing), 3: semi-isotropic NPT, 4: isotropic NPT, 5: anisotropic NPT (under development)
-nsteps 100 ! number of simulation steps
+simulation_name water_6540_notab_nve_qeq_111_mpi_gpu_probe ! output files will carry this name + their specific ext.
+ensemble_type 0 ! 0: NVE, 1: Berendsen NVT, 2: Nose-Hoover NVT(under testing), 3: semi-isotropic NPT, 4: isotropic NPT, 5: anisotropic NPT (under development)
+nsteps 1000 ! number of simulation steps
dt 0.25 ! time step in fs
proc_by_dim 1 1 1 ! distribution of processors by dimensions
gpus_per_node 1 ! GPUs per node
@@ -9,14 +9,14 @@ reposition_atoms 0 ! 0: just fit to periodic bounda
restrict_bonds 0 ! enforce the bonds given in CONECT lines of pdb file for this many steps
tabulate_long_range 0 ! number of sampling points for cubic spline interpolation, 0 no interpolation
energy_update_freq 1
-remove_CoM_vel 500 ! remove the translational and rotational vel around the center of mass at every 'this many' steps
+remove_CoM_vel 0 ! remove the translational and rotational vel around the center of mass at every 'this many' steps
reneighbor 1
-vlist_buffer 0
+vlist_buffer 0.0
nbrhood_cutoff 5.0 ! near neighbors cutoff for bond calculations (Angstroms)
bond_graph_cutoff 0.3 ! bond strength cutoff for bond graphs (Angstroms)
thb_cutoff 0.001 ! cutoff value for three body interactions (Angstroms)
-hbond_cutoff 7.50 ! cutoff distance for hydrogen bond interactions (Angstroms)
+hbond_cutoff 7.5 ! cutoff distance for hydrogen bond interactions (Angstroms)
charge_method 0 ! charge method: 0 = QEq, 1 = EEM, 2 = ACKS2
charge_freq 1 ! frequency (sim step) at which atomic charges are computed
@@ -24,19 +24,19 @@ cm_q_net 0.0 ! net system charge
cm_solver_type 2 ! iterative linear solver for charge method: 0 = GMRES, 1 = GMRES_H, 2 = CG, 3 = SDM
cm_solver_max_iters 1000 ! max solver iterations
cm_solver_restart 100 ! inner iterations of GMRES before restarting
-cm_solver_q_err 1e-6 ! relative residual norm threshold used in solver
+cm_solver_q_err 1.0e-14 ! relative residual norm threshold used in solver
cm_domain_sparsity 1.0 ! scalar for scaling cut-off distance, used to sparsify charge matrix (between 0.0 and 1.0)
cm_solver_pre_comp_type 1 ! method used to compute preconditioner, if applicable
-cm_solver_pre_comp_refactor 1000 ! number of steps before recomputing preconditioner
+cm_solver_pre_comp_refactor 1 ! number of steps before recomputing preconditioner
cm_solver_pre_comp_droptol 0.0 ! threshold tolerance for dropping values in preconditioner computation, if applicable
cm_solver_pre_comp_sweeps 3 ! number of sweeps used to compute preconditioner (ILU_PAR)
cm_solver_pre_app_type 1 ! method used to apply preconditioner
cm_solver_pre_app_jacobi_iters 50 ! number of Jacobi iterations used for applying precondition, if applicable
-temp_init 0.01 ! desired initial temperature of the simulated system
-temp_final 300.0 ! desired final temperature of the simulated system
-t_mass 500.0 ! thermal inertia parameter (fs): Nose-Hoover-NVT: 0.16666, NVP: 100.0, bNVT/iNPT/sNPT: 500.0
-t_mode 2 ! 0: T-coupling only, 1: step-wise, 2: constant slope
+temp_init 0.0 ! desired initial temperature of the simulated system
+temp_final 0.0 ! desired final temperature of the simulated system
+t_mass 0.16666 ! thermal inertia parameter (fs): Nose-Hoover-NVT: 0.16666, NVP: 100.0, bNVT/iNPT/sNPT: 500.0
+t_mode 0 ! 0: T-coupling only, 1: step-wise, 2: constant slope
t_rate 5.0 ! in K
t_freq 1.0 ! in ps
@@ -45,20 +45,20 @@ p_mass 10000.00 10000.00 10000.00 ! in fs, pressure
compress 0.008134 ! in ps^2 * A / amu ( 4.5X10^(-5) bar^(-1) )
press_mode 0 ! 0: internal + external pressure, 1: ext only, 2: int only
-geo_format 1 ! 0: custom 1: pdb (only if natoms < 100000) 2: ASCII restart 3: binary restart
+geo_format 1 ! 0: custom 1: pdb (only if natoms < 100000) 2: bgf 3: ASCII restart 3: binary restart
write_freq 0 ! write trajectory after so many steps
-traj_method 1 ! 0: simple parallel I/O, 1: MPI I/O
-traj_title WATER_NVE ! (no white spaces)
+traj_method 0 ! 0: simple parallel I/O, 1: MPI I/O
+traj_title water_6540_notab_nve_qeq_111_mpi_gpu_probe ! (no white spaces)
atom_info 1 ! 0: no atom info, 1: print basic atom info in the trajectory file
atom_forces 0 ! 0: basic atom format, 1: print force on each atom in the trajectory file
atom_velocities 0 ! 0: basic atom format, 1: print the velocity of each atom in the trajectory file
-bond_info 1 ! 0: do not print bonds, 1: print bonds in the trajectory file
-angle_info 1 ! 0: do not print angles, 1: print angles in the trajectory file
+bond_info 0 ! 0: do not print bonds, 1: print bonds in the trajectory file
+angle_info 0 ! 0: do not print angles, 1: print angles in the trajectory file
dipole_anal 0 ! 1: calculate a electric dipole moment of the system
-freq_dipole_anal 1 ! calculate electric dipole moment at every 'this many' steps
+freq_dipole_anal 0 ! calculate electric dipole moment at every 'this many' steps
diffusion_coef 0 ! 1: calculate diffusion coefficient of the system
-freq_diffusion_coef 1 ! calculate diffusion coefficient at every 'this many' steps
+freq_diffusion_coef 0 ! calculate diffusion coefficient at every 'this many' steps
restrict_type 2 ! -1: all types of atoms, 0 and up: only this type of atoms
restart_format 1 ! 0: restarts in ASCII 1: restarts in binary
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index 766ea766c7fd495267db7e17b6b8e1e1b617c970..c4fc219b64ba54aa3416ee7b95b333ee4729ed68 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -53,7 +53,8 @@ void Velocity_Verlet_NVE( reax_system *system, control_params *control,
rvec_ScaledSum( dx, control->dt, system->atoms[i].v,
scalar2 * inv_m, system->atoms[i].f );
- control->update_atom_position( system->atoms[i].x, dx, system->atoms[i].rel_map, &system->box );
+ control->update_atom_position( system->atoms[i].x, dx,
+ system->atoms[i].rel_map, &system->box );
/* Compute v(t + dt/2) */
rvec_ScaledAdd( system->atoms[i].v,
@@ -307,6 +308,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
for ( i = 0; i < system->N; i++ )
{
inv_m = 1.0 / system->reax_param.sbp[system->atoms[i].type].mass;
+
/* Compute v(t + dt) */
rvec_ScaledAdd( system->atoms[i].v,
-0.5 * dt * F_CONV * inv_m, system->atoms[i].f );
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index b030ac672603bf492cb22a24337ed6bcbecabd38..58decd7519e0d3875ce5cde9bf9c3a355883f8e7 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -38,10 +38,10 @@
/* enables debugging code */
//#define DEBUG_FOCUS
-/* enables test forces code */
-//#define TEST_FORCES
/* enables test energy code */
//#define TEST_ENERGY
+/* enables test forces code */
+//#define TEST_FORCES
/* constants defined in reference Fortran ReaxFF code (useful for comparisons) */
#define USE_REF_FORTRAN_REAXFF_CONSTANTS
/* constants defined in reference Fortran eReaxFF code (useful for comparisons) */
@@ -70,7 +70,7 @@
#if defined(USE_REF_FORTRAN_REAXFF_CONSTANTS)
/* transcendental constant pi */
#define PI (3.14159265)
- /* uni conversion from ??? to kcal / mol */
+ /* unit conversion from ??? to kcal / mol */
#define C_ELE (332.0638)
/* Boltzmann constant, AMU * A^2 / (ps^2 * K) */
#define K_B (0.831687)
@@ -163,14 +163,22 @@
#define P_CONV (1.0e-24 * AVOGNR * JOULES_to_CAL)
#endif
+/* max. num. of characters for string buffers */
#define MAX_STR (1024)
+/* max. num. of characters for a line in files */
#define MAX_LINE (1024)
+/* max. num. of tokens per line */
#define MAX_TOKENS (1024)
+/* max. num. of characters per token */
#define MAX_TOKEN_LEN (1024)
+/* max. atom ID in geo file */
#define MAX_ATOM_ID (100000)
+/* ??? */
#define MAX_RESTRICT (15)
+/* max. num. atoms per molecule */
#define MAX_MOLECULE_SIZE (20)
+/* max. num. atom types defined in the force field parameter file */
#define MAX_ATOM_TYPES (25)
#define MAX_GRID (50)
@@ -178,9 +186,13 @@
#define MAX_4BODY_PARAM (5)
#define NUM_INTRS (10)
+/* max. pressure scaler for simulation box dimenion in NPT ensembles */
#define MAX_dV (1.01)
+/* min. pressure scaler for simulation box dimenion in NPT ensembles */
#define MIN_dV (1.0 / MAX_dV)
+/* max. temperature scaler for atomic positions and velocities in NPT ensembles */
#define MAX_dT (2.0)
+/* min. temperature scaler for atomic positions and velocities in NPT ensembles */
#define MIN_dT (0.0)
#define ZERO (0.000000000000000e+00)
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 02e1b1bc51f3b684f319d4e0c6c124fa4182c37c..83eb2a24b607e0a653dccb4e7bb775a079c0567f 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -238,21 +238,28 @@ real Get_Timing_Info( real t_start )
/*********** from io_tools.c **************/
int Get_Atom_Type( reax_interaction *reax_param, char *s, size_t n )
{
- int i;
+ int i, ret, flag;
+
+ flag = FAILURE;
for ( i = 0; i < reax_param->num_atom_types; ++i )
{
if ( strncmp( reax_param->sbp[i].name, s,
MIN( sizeof(reax_param->sbp[i].name), n ) ) == 0 )
{
- return i;
+ ret = i;
+ flag = SUCCESS;
+ break;
}
}
- fprintf( stderr, "[ERROR] Unknown atom type: %s. Terminating...\n", s );
- exit( UNKNOWN_ATOM_TYPE );
+ if ( flag == FAILURE )
+ {
+ fprintf( stderr, "[ERROR] Unknown atom type: %s. Terminating...\n", s );
+ exit( UNKNOWN_ATOM_TYPE );
+ }
- return FAILURE;
+ return ret;
}