-
Kurt A. O'Hearn authored
PG-PuReMD: merge memory management and step retry logic from MPI+GPU code into MPI-only code. Fix and backport numerous bugs (void pointer arithmetic and dereferences, etc.). Update code to match new data structures. Begin merging changes from OpenMP code (sPuReMD) for preconditioning work.
Kurt A. O'Hearn authoredPG-PuReMD: merge memory management and step retry logic from MPI+GPU code into MPI-only code. Fix and backport numerous bugs (void pointer arithmetic and dereferences, etc.). Update code to match new data structures. Begin merging changes from OpenMP code (sPuReMD) for preconditioning work.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
ffield.c 25.51 KiB
/*----------------------------------------------------------------------
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"
#if defined(PURE_REAX)
#include "ffield.h"
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reax_ffield.h"
#include "reax_tool_box.h"
#endif
void Read_Force_Field( char *ffield_file, reax_interaction *reax,
reax_system *system, control_params *control )
{
FILE *fp;
char *s;
char **tmp;
char *tor_flag;
int c, i, j, k, l, m, n, o, p, cnt;
real val;
int __N;
int index1, index2;
/* open force field file */
if ( (fp = fopen( ffield_file, "r" ) ) == NULL )
{
fprintf( stderr, "[ERROR] p%d: cannot open force field file! terminating...\n",
system->my_rank );
MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
}
s = smalloc( sizeof(char) * MAX_LINE, "Read_Force_Field::s" );
tmp = smalloc( sizeof(char *) * MAX_TOKENS, "Read_Force_Field::tmp");
for (i = 0; i < MAX_TOKENS; i++)
{
tmp[i] = smalloc( sizeof(char) * MAX_TOKEN_LEN, "Read_Force_Field::tmp[i]" );
}
/* reading first header comment */
fgets( s, MAX_LINE, fp );
/* line 2 is number of global parameters */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
/* reading the number of global parameters */
n = atoi(tmp[0]);
if ( n < 1 )
{ fprintf( stderr, "[WARNING] p%d: number of globals in ffield file is 0!\n",
system->my_rank );
return;
}
reax->gp.n_global = n;
reax->gp.l = smalloc( sizeof(real) * n, "Read_Force_Field::reax->gp.l" );
/* see reax_types.h for mapping between l[i] and the lambdas used in ff */
for (i = 0; i < n; i++)
{
fgets(s, MAX_LINE, fp);
c = Tokenize(s, &tmp);
val = (real) atof(tmp[0]);
reax->gp.l[i] = val;
}
control->bo_cut = 0.01 * reax->gp.l[29];
control->nonb_low = reax->gp.l[11];
control->nonb_cut = reax->gp.l[12];
/* next line is number of atom types and some comments */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
reax->num_atom_types = atoi(tmp[0]);
/* 3 lines of comments */
fgets(s, MAX_LINE, fp);
fgets(s, MAX_LINE, fp);
fgets(s, MAX_LINE, fp);
/* Allocating structures in reax_interaction */
__N = reax->num_atom_types;
reax->sbp = scalloc( reax->num_atom_types, sizeof(single_body_parameters),
"Read_Force_Field::reax->sbp" );
reax->tbp = scalloc( POW(reax->num_atom_types, 2.0), sizeof(two_body_parameters),
"Read_Force_Field::reax->tbp" );
reax->thbp = scalloc( POW(reax->num_atom_types, 3.0), sizeof(three_body_header),
"Read_Force_Field::reax->thbp" );
reax->hbp = scalloc( POW(reax->num_atom_types, 3.0), sizeof(hbond_parameters),
"Read_Force_Field::reax->hbp" );
reax->fbp = scalloc( POW(reax->num_atom_types, 4.0), sizeof(four_body_header),
"Read_Force_Field::reax->fbp" );
tor_flag = scalloc( POW(reax->num_atom_types, 4.0), sizeof(char),
"Read_Force_Field::tor_flag" );
/* vdWaals type:
* 1: Shielded Morse, no inner-wall
* 2: inner wall, no shielding
* 3: inner wall+shielding */
reax->gp.vdw_type = 0;
/* reading single atom parameters */
/* there are 4 lines of each single atom parameters in ff files. these
* parameters later determine some of the pair and triplet parameters using
* combination rules. */
for ( i = 0; i < reax->num_atom_types; i++ )
{
/* line one */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
for ( j = 0; j < (int)(strlen(tmp[0])); ++j ) {
reax->sbp[i].name[j] = toupper( tmp[0][j] );
}
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: Atom Name in the force field : %s \n",
system->my_rank, reax->sbp[i].name );
#endif
val = atof(tmp[1]);
reax->sbp[i].r_s = val;
val = atof(tmp[2]);
reax->sbp[i].valency = val;
val = atof(tmp[3]);
reax->sbp[i].mass = val;
val = atof(tmp[4]);
reax->sbp[i].r_vdw = val;
val = atof(tmp[5]);
reax->sbp[i].epsilon = val;
val = atof(tmp[6]);
reax->sbp[i].gamma = val;
val = atof(tmp[7]);
reax->sbp[i].r_pi = val;
val = atof(tmp[8]);
reax->sbp[i].valency_e = val;
reax->sbp[i].nlp_opt = 0.5 * (reax->sbp[i].valency_e - reax->sbp[i].valency);
/* line two */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
val = atof(tmp[0]);
reax->sbp[i].alpha = val;
val = atof(tmp[1]);
reax->sbp[i].gamma_w = val;
val = atof(tmp[2]);
reax->sbp[i].valency_boc = val;
val = atof(tmp[3]);
reax->sbp[i].p_ovun5 = val;
val = atof(tmp[4]);
val = atof(tmp[5]);
reax->sbp[i].chi = val;
val = atof(tmp[6]);
reax->sbp[i].eta = 2.0 * val;
val = atof(tmp[7]);
reax->sbp[i].p_hbond = (int)val;
/* line 3 */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
val = atof(tmp[0]);
reax->sbp[i].r_pi_pi = val;
val = atof(tmp[1]);
reax->sbp[i].p_lp2 = val;
val = atof(tmp[2]);
val = atof(tmp[3]);
reax->sbp[i].b_o_131 = val;
val = atof(tmp[4]);
reax->sbp[i].b_o_132 = val;
val = atof(tmp[5]);
reax->sbp[i].b_o_133 = val;
val = atof(tmp[6]);
val = atof(tmp[7]);
/* line 4 */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
val = atof(tmp[0]); reax->sbp[i].p_ovun2 = val;
val = atof(tmp[1]);
reax->sbp[i].p_val3 = val;
val = atof(tmp[2]);
val = atof(tmp[3]);
reax->sbp[i].valency_val = val;
val = atof(tmp[4]);
reax->sbp[i].p_val5 = val;
val = atof(tmp[5]);
reax->sbp[i].rcore2 = val;
val = atof(tmp[6]);
reax->sbp[i].ecore2 = val;
val = atof(tmp[7]);
reax->sbp[i].acore2 = val;
/* Inner-wall */
if ( reax->sbp[i].rcore2 > 0.01 && reax->sbp[i].acore2 > 0.01 )
{
/* Shielding vdWaals */
if ( reax->sbp[i].gamma_w > 0.5 )
{
if ( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3 )
{
fprintf( stderr, "[WARNING] p%d: inconsistent vdWaals-parameters\n"
"Force field parameters for element %s\n"
"indicate inner wall+shielding, but earlier\n"
"atoms indicate different vdWaals-method.\n"
"This may cause division-by-zero errors.\n"
"Keeping vdWaals-setting for earlier atoms.\n",
system->my_rank, reax->sbp[i].name );
}
else
{
reax->gp.vdw_type = 3;
#if defined(DEBUG)
fprintf( stderr, "p%d: vdWaals type for element %s: Shielding+inner-wall",
system->my_rank, reax->sbp[i].name );
#endif
}
}
/* No shielding vdWaals parameters present */
else
{
if ( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 2 )
{
fprintf( stderr, "[WARNING] p%d: inconsistent vdWaals-parameters\n",
system->my_rank );
fprintf( stderr, " [INFO] Force field parameters for element %s\n", reax->sbp[i].name );
fprintf( stderr, " [INFO] indicate inner wall without shielding, but earlier\n" );
fprintf( stderr, " [INFO] atoms indicate different vdWaals-method.\n" );
fprintf( stderr, " [INFO] This may cause division-by-zero errors.\n" );
fprintf( stderr, " [INFO] Keeping vdWaals-setting for earlier atoms.\n" );
}
else
{
reax->gp.vdw_type = 2;
#if defined(DEBUG)
fprintf( stderr, "p%d: vdWaals type for element%s: No Shielding,inner-wall",
system->my_rank, reax->sbp[i].name );
#endif
}
}
}
/* No Inner wall parameters present */
else
{
/* Shielding vdWaals */
if ( reax->sbp[i].gamma_w > 0.5 )
{
if ( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 1 ) fprintf( stderr, "[WARNING] p%d: inconsistent vdWaals-parameters\n" \
" [INFO] Force field parameters for element %s\n" \
" [INFO] indicate shielding without inner wall, but earlier\n" \
" [INFO] atoms indicate different vdWaals-method.\n" \
" [INFO] This may cause division-by-zero errors.\n" \
" [INFO] Keeping vdWaals-setting for earlier atoms.\n",
system->my_rank, reax->sbp[i].name );
else
{
reax->gp.vdw_type = 1;
#if defined(DEBUG)
fprintf( stderr, "p%d, vdWaals type for element%s: Shielding,no inner-wall",
system->my_rank, reax->sbp[i].name );
#endif
}
}
else
{
fprintf( stderr, "[ERROR] p%d: inconsistent vdWaals-parameters\n" \
" [INFO] No shielding or inner-wall set for element %s\n",
system->my_rank, reax->sbp[i].name );
MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
}
}
}
#if defined(DEBUG)
fprintf( stderr, "p%d: vdWaals type: %d\n", system->my_rank, reax->gp.vdw_type );
#endif
/* Equate vval3 to valf for first-row elements (25/10/2004) */
for ( i = 0; i < reax->num_atom_types; i++ )
{
if ( reax->sbp[i].mass < 21 &&
reax->sbp[i].valency_val != reax->sbp[i].valency_boc )
{
fprintf( stderr, "[WARNING] p%d: changed valency_val to valency_boc for atom type %s\n",
system->my_rank, reax->sbp[i].name );
reax->sbp[i].valency_val = reax->sbp[i].valency_boc;
}
}
/* next line is number of two body combination and some comments */
fgets(s, MAX_LINE, fp);
c = Tokenize(s, &tmp);
l = atoi(tmp[0]);
/* a line of comments */
fgets(s, MAX_LINE, fp);
for (i = 0; i < l; i++)
{
/* line 1 */
fgets(s, MAX_LINE, fp);
c = Tokenize(s, &tmp);
j = atoi(tmp[0]) - 1;
k = atoi(tmp[1]) - 1;
index1 = j * __N + k;
index2 = k * __N + j;
if (j < reax->num_atom_types && k < reax->num_atom_types)
{
val = atof(tmp[2]);
reax->tbp[ index1 ].De_s = val;
reax->tbp[ index2 ].De_s = val;
val = atof(tmp[3]);
reax->tbp[ index1 ].De_p = val;
reax->tbp[ index2 ].De_p = val;
val = atof(tmp[4]); reax->tbp[ index1 ].De_pp = val;
reax->tbp[ index2 ].De_pp = val;
val = atof(tmp[5]);
reax->tbp[ index1 ].p_be1 = val;
reax->tbp[ index2 ].p_be1 = val;
val = atof(tmp[6]);
reax->tbp[ index1 ].p_bo5 = val;
reax->tbp[ index2 ].p_bo5 = val;
val = atof(tmp[7]);
reax->tbp[ index1 ].v13cor = val;
reax->tbp[ index2 ].v13cor = val;
val = atof(tmp[8]);
reax->tbp[ index1 ].p_bo6 = val;
reax->tbp[ index2 ].p_bo6 = val;
val = atof(tmp[9]);
reax->tbp[ index1 ].p_ovun1 = val;
reax->tbp[ index2 ].p_ovun1 = val;
/* line 2 */
fgets(s, MAX_LINE, fp);
c = Tokenize(s, &tmp);
val = atof(tmp[0]);
reax->tbp[ index1 ].p_be2 = val;
reax->tbp[ index2 ].p_be2 = val;
val = atof(tmp[1]);
reax->tbp[ index1 ].p_bo3 = val;
reax->tbp[ index2 ].p_bo3 = val;
val = atof(tmp[2]);
reax->tbp[ index1 ].p_bo4 = val;
reax->tbp[ index2 ].p_bo4 = val;
val = atof(tmp[3]);
val = atof(tmp[4]);
reax->tbp[ index1 ].p_bo1 = val;
reax->tbp[ index2 ].p_bo1 = val;
val = atof(tmp[5]);
reax->tbp[ index1 ].p_bo2 = val;
reax->tbp[ index2 ].p_bo2 = val;
val = atof(tmp[6]);
reax->tbp[ index1 ].ovc = val;
reax->tbp[ index2 ].ovc = val;
val = atof(tmp[7]);
}
}
/* calculating combination rules and filling up remaining fields. */
for (i = 0; i < reax->num_atom_types; i++)
{
for (j = i; j < reax->num_atom_types; j++)
{
index1 = i * __N + j;
index2 = j * __N + i;
reax->tbp[index1].r_s =
0.5 * (reax->sbp[i].r_s + reax->sbp[j].r_s);
reax->tbp[index2].r_s =
0.5 * (reax->sbp[j].r_s + reax->sbp[i].r_s);
reax->tbp[index1].r_p =
0.5 * (reax->sbp[i].r_pi + reax->sbp[j].r_pi);
reax->tbp[index2].r_p =
0.5 * (reax->sbp[j].r_pi + reax->sbp[i].r_pi);
reax->tbp[index1].r_pp =
0.5 * (reax->sbp[i].r_pi_pi + reax->sbp[j].r_pi_pi);
reax->tbp[index2].r_pp =
0.5 * (reax->sbp[j].r_pi_pi + reax->sbp[i].r_pi_pi);
reax->tbp[index1].p_boc3 =
SQRT(reax->sbp[i].b_o_132 * reax->sbp[j].b_o_132);
reax->tbp[index2].p_boc3 =
SQRT(reax->sbp[j].b_o_132 * reax->sbp[i].b_o_132);
reax->tbp[index1].p_boc4 =
SQRT(reax->sbp[i].b_o_131 * reax->sbp[j].b_o_131);
reax->tbp[index2].p_boc4 =
SQRT(reax->sbp[j].b_o_131 * reax->sbp[i].b_o_131);
reax->tbp[index1].p_boc5 =
SQRT(reax->sbp[i].b_o_133 * reax->sbp[j].b_o_133);
reax->tbp[index2].p_boc5 =
SQRT(reax->sbp[j].b_o_133 * reax->sbp[i].b_o_133);
reax->tbp[index1].D =
SQRT(reax->sbp[i].epsilon * reax->sbp[j].epsilon);
reax->tbp[index2].D =
SQRT(reax->sbp[j].epsilon * reax->sbp[i].epsilon);
reax->tbp[index1].alpha =
SQRT(reax->sbp[i].alpha * reax->sbp[j].alpha);
reax->tbp[index2].alpha =
SQRT(reax->sbp[j].alpha * reax->sbp[i].alpha);
reax->tbp[index1].r_vdW =
2.0 * SQRT(reax->sbp[i].r_vdw * reax->sbp[j].r_vdw);
reax->tbp[index2].r_vdW =
2.0 * SQRT(reax->sbp[j].r_vdw * reax->sbp[i].r_vdw);
reax->tbp[index1].gamma_w =
SQRT(reax->sbp[i].gamma_w * reax->sbp[j].gamma_w);
reax->tbp[index2].gamma_w =
SQRT(reax->sbp[j].gamma_w * reax->sbp[i].gamma_w);
reax->tbp[index1].gamma =
POW(reax->sbp[i].gamma * reax->sbp[j].gamma, -1.5);
reax->tbp[index2].gamma =
POW(reax->sbp[j].gamma * reax->sbp[i].gamma, -1.5);
/* additions for additional vdWaals interaction types - inner core */
reax->tbp[index1].rcore = reax->tbp[index2].rcore =
SQRT( reax->sbp[i].rcore2 * reax->sbp[j].rcore2 );
reax->tbp[index1].ecore = reax->tbp[index2].ecore =
SQRT( reax->sbp[i].ecore2 * reax->sbp[j].ecore2 );
reax->tbp[index1].acore = reax->tbp[index2].acore =
SQRT( reax->sbp[i].acore2 * reax->sbp[j].acore2 );
}
}
/* next line is number of two body offdiagonal combinations and comments */
/* these are two body offdiagonal terms that are different from the
combination rules defined above */
fgets(s, MAX_LINE, fp);
c = Tokenize(s, &tmp);
l = atoi(tmp[0]);
for (i = 0; i < l; i++)
{
fgets(s, MAX_LINE, fp);
c = Tokenize(s, &tmp);
j = atoi(tmp[0]) - 1;
k = atoi(tmp[1]) - 1;
index1 = j * __N + k;
index2 = k * __N + j;
if (j < reax->num_atom_types && k < reax->num_atom_types)
{
val = atof(tmp[2]);
if (val > 0.0)
{
reax->tbp[index1].D = val;
reax->tbp[index2].D = val;
}
val = atof(tmp[3]);
if (val > 0.0)
{
reax->tbp[index1].r_vdW = 2 * val;
reax->tbp[index2].r_vdW = 2 * val;
}
val = atof(tmp[4]);
if (val > 0.0)
{
reax->tbp[index1].alpha = val;
reax->tbp[index2].alpha = val;
}
val = atof(tmp[5]);
if (val > 0.0)
{
reax->tbp[index1].r_s = val;
reax->tbp[index2].r_s = val;
}
val = atof(tmp[6]);
if (val > 0.0)
{
reax->tbp[index1].r_p = val;
reax->tbp[index2].r_p = val;
}
val = atof(tmp[7]);
if (val > 0.0)
{
reax->tbp[index1].r_pp = val;
reax->tbp[index2].r_pp = val;
}
}
}
/* 3-body parameters -
supports multi-well potentials (upto MAX_3BODY_PARAM in mytypes.h) */
/* clear entries first */
for ( i = 0; i < reax->num_atom_types; ++i )
{
for ( j = 0; j < reax->num_atom_types; ++j )
{
for ( k = 0; k < reax->num_atom_types; ++k )
{
reax->thbp[i * __N * __N + j * __N + k].cnt = 0;
}
}
}
/* next line is number of 3-body params and some comments */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
l = atoi( tmp[0] );
for ( i = 0; i < l; i++ )
{
fgets(s, MAX_LINE, fp);
c = Tokenize(s, &tmp);
j = atoi(tmp[0]) - 1;
k = atoi(tmp[1]) - 1;
m = atoi(tmp[2]) - 1;
index1 = j * __N * __N + k * __N + m;
index2 = m * __N * __N + k * __N + j;
if (j < reax->num_atom_types && k < reax->num_atom_types &&
m < reax->num_atom_types)
{
cnt = reax->thbp[index1].cnt;
reax->thbp[index1].cnt++;
reax->thbp[index2].cnt++;
val = atof(tmp[3]);
reax->thbp[index1].prm[cnt].theta_00 = val;
reax->thbp[index2].prm[cnt].theta_00 = val;
val = atof(tmp[4]);
reax->thbp[index1].prm[cnt].p_val1 = val;
reax->thbp[index2].prm[cnt].p_val1 = val;
val = atof(tmp[5]);
reax->thbp[index1].prm[cnt].p_val2 = val;
reax->thbp[index2].prm[cnt].p_val2 = val;
val = atof(tmp[6]);
reax->thbp[index1].prm[cnt].p_coa1 = val;
reax->thbp[index2].prm[cnt].p_coa1 = val;
val = atof(tmp[7]);
reax->thbp[index1].prm[cnt].p_val7 = val;
reax->thbp[index2].prm[cnt].p_val7 = val;
val = atof(tmp[8]);
reax->thbp[index1].prm[cnt].p_pen1 = val;
reax->thbp[index2].prm[cnt].p_pen1 = val;
val = atof(tmp[9]);
reax->thbp[index1].prm[cnt].p_val4 = val;
reax->thbp[index2].prm[cnt].p_val4 = val;
}
}
/* 4-body parameters are entered in compact form. i.e. 0-X-Y-0
* correspond to any type of pair of atoms in 1 and 4
* position. However, explicit X-Y-Z-W takes precedence over the
* default description.
* supports multi-well potentials (upto MAX_4BODY_PARAM in mytypes.h)
* IMPORTANT: for now, directions on how to read multi-entries from ffield
* is not clear */
/* clear all entries first */
for ( i = 0; i < reax->num_atom_types; ++i )
{
for ( j = 0; j < reax->num_atom_types; ++j )
{
for ( k = 0; k < reax->num_atom_types; ++k )
{
for ( m = 0; m < reax->num_atom_types; ++m )
{
reax->fbp[i * __N * __N * __N + j * __N * __N + k * __N + m].cnt = 0;
tor_flag[i * __N * __N * __N + j * __N * __N + k * __N + m] = 0;
} }
}
}
/* next line is number of 4-body params and some comments */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
l = atoi( tmp[0] );
for ( i = 0; i < l; i++ )
{
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
j = atoi(tmp[0]) - 1;
k = atoi(tmp[1]) - 1;
m = atoi(tmp[2]) - 1;
n = atoi(tmp[3]) - 1;
index1 = j * __N * __N * __N + k * __N * __N + m * __N + n;
index2 = n * __N * __N * __N + m * __N * __N + k * __N + j;
/* this means the entry is not in compact form */
if ( j >= 0 && n >= 0 )
{
if ( j < reax->num_atom_types && k < reax->num_atom_types &&
m < reax->num_atom_types && n < reax->num_atom_types )
{
/* these flags ensure that this entry take precedence
over the compact form entries */
tor_flag[index1] = 1;
tor_flag[index2] = 1;
reax->fbp[index1].cnt = 1;
reax->fbp[index2].cnt = 1;
val = atof(tmp[4]);
reax->fbp[index1].prm[0].V1 = val;
reax->fbp[index2].prm[0].V1 = val;
val = atof(tmp[5]);
reax->fbp[index1].prm[0].V2 = val;
reax->fbp[index2].prm[0].V2 = val;
val = atof(tmp[6]);
reax->fbp[index1].prm[0].V3 = val;
reax->fbp[index2].prm[0].V3 = val;
val = atof(tmp[7]);
reax->fbp[index1].prm[0].p_tor1 = val;
reax->fbp[index2].prm[0].p_tor1 = val;
val = atof(tmp[8]);
reax->fbp[index1].prm[0].p_cot1 = val;
reax->fbp[index2].prm[0].p_cot1 = val;
}
}
/* This means the entry is of the form 0-X-Y-0 */
else
{
if ( k < reax->num_atom_types && m < reax->num_atom_types )
{
for ( p = 0; p < reax->num_atom_types; p++ )
{
for ( o = 0; o < reax->num_atom_types; o++ )
{
index1 = p * __N * __N * __N + k * __N * __N + m * __N + o;
index2 = o * __N * __N * __N + m * __N * __N + k * __N + p;
reax->fbp[index1].cnt = 1;
reax->fbp[index2].cnt = 1;
if (tor_flag[index1] == 0)
{
reax->fbp[index1].prm[0].V1 = atof(tmp[4]);
reax->fbp[index1].prm[0].V2 = atof(tmp[5]);
reax->fbp[index1].prm[0].V3 = atof(tmp[6]);
reax->fbp[index1].prm[0].p_tor1 = atof(tmp[7]);
reax->fbp[index1].prm[0].p_cot1 = atof(tmp[8]);
}
if (tor_flag[index2] == 0)
{
reax->fbp[index2].prm[0].V1 = atof(tmp[4]);
reax->fbp[index2].prm[0].V2 = atof(tmp[5]);
reax->fbp[index2].prm[0].V3 = atof(tmp[6]);
reax->fbp[index2].prm[0].p_tor1 = atof(tmp[7]);
reax->fbp[index2].prm[0].p_cot1 = atof(tmp[8]);
}
}
}
}
}
}
/* next line is number of hydrogen bond params and some comments */
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
l = atoi( tmp[0] );
for ( i = 0; i < l; i++ )
{
fgets( s, MAX_LINE, fp );
c = Tokenize( s, &tmp );
j = atoi(tmp[0]) - 1;
k = atoi(tmp[1]) - 1;
m = atoi(tmp[2]) - 1;
index1 = j * __N * __N + k * __N + m;
if ( j < reax->num_atom_types && m < reax->num_atom_types )
{
val = atof(tmp[3]);
reax->hbp[index1].r0_hb = val;
val = atof(tmp[4]);
reax->hbp[index1].p_hb1 = val;
val = atof(tmp[5]);
reax->hbp[index1].p_hb2 = val;
val = atof(tmp[6]);
reax->hbp[index1].p_hb3 = val;
}
}
/* deallocate helper storage */
for ( i = 0; i < MAX_TOKENS; i++ )
{
sfree( tmp[i], "Read_Force_Field::tmp[i]" );
}
sfree( tmp, "Read_Force_Field::tmp" );
sfree( s, "Read_Force_Field::s" );
sfree( tor_flag, "Read_Force_Field::tor_flag" );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: force field read\n", system->my_rank );
#endif
}