Homec4science
Diffusion lammps (lammm-master) Browse

lammps/src/USER-REAXC/reaxc_ffield.cpp
3d36d38a5183lammm-master

Last Change
Raw File

reaxc_ffield.cpp
View Options

/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Copyright (2010) Purdue University
Hasan Metin Aktulga, hmaktulga@lbl.gov
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
Please cite the related publication:
H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama,
"Parallel Reactive Molecular Dynamics: Numerical Methods and
Algorithmic Techniques", Parallel Computing, in press.
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 "pair_reax_c.h"
#if defined(PURE_REAX)
#include "ffield.h"
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reaxc_ffield.h"
#include "reaxc_tool_box.h"
#endif
char Read_Force_Field( char *ffield_file, reax_interaction *reax,
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;
MPI_Comm comm;
comm = MPI_COMM_WORLD;
/* open force field file */
if ( (fp = fopen( ffield_file, "r" ) ) == NULL ) {
fprintf( stderr, "error opening the force filed file! terminating...\n" );
MPI_Abort( comm, FILE_NOT_FOUND );
}
s = (char*) malloc(sizeof(char)*MAX_LINE);
tmp = (char**) malloc(sizeof(char*)*MAX_TOKENS);
for (i=0; i < MAX_TOKENS; i++)
tmp[i] = (char*) malloc(sizeof(char)*MAX_TOKEN_LEN);
/* 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: number of globals in ffield file is 0!\n" );
fclose(fp);
return 1;
}
reax->gp.n_global = n;
reax->gp.l = (real*) malloc(sizeof(real)*n);
/* 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 */
reax->sbp = (single_body_parameters*)
scalloc( reax->num_atom_types, sizeof(single_body_parameters), "sbp",
comm );
reax->tbp = (two_body_parameters**)
scalloc( reax->num_atom_types, sizeof(two_body_parameters*), "tbp", comm );
reax->thbp= (three_body_header***)
scalloc( reax->num_atom_types, sizeof(three_body_header**), "thbp", comm );
reax->hbp = (hbond_parameters***)
scalloc( reax->num_atom_types, sizeof(hbond_parameters**), "hbp", comm );
reax->fbp = (four_body_header****)
scalloc( reax->num_atom_types, sizeof(four_body_header***), "fbp", comm );
tor_flag = (char****)
scalloc( reax->num_atom_types, sizeof(char***), "tor_flag", comm );
for( i = 0; i < reax->num_atom_types; i++ ) {
reax->tbp[i] = (two_body_parameters*)
scalloc( reax->num_atom_types, sizeof(two_body_parameters), "tbp[i]",
comm );
reax->thbp[i]= (three_body_header**)
scalloc( reax->num_atom_types, sizeof(three_body_header*), "thbp[i]",
comm );
reax->hbp[i] = (hbond_parameters**)
scalloc( reax->num_atom_types, sizeof(hbond_parameters*), "hbp[i]",
comm );
reax->fbp[i] = (four_body_header***)
scalloc( reax->num_atom_types, sizeof(four_body_header**), "fbp[i]",
comm );
tor_flag[i] = (char***)
scalloc( reax->num_atom_types, sizeof(char**), "tor_flag[i]", comm );
for( j = 0; j < reax->num_atom_types; j++ ) {
reax->thbp[i][j]= (three_body_header*)
scalloc( reax->num_atom_types, sizeof(three_body_header), "thbp[i,j]",
comm );
reax->hbp[i][j] = (hbond_parameters*)
scalloc( reax->num_atom_types, sizeof(hbond_parameters), "hbp[i,j]",
comm );
reax->fbp[i][j] = (four_body_header**)
scalloc( reax->num_atom_types, sizeof(four_body_header*), "fbp[i,j]",
comm );
tor_flag[i][j] = (char**)
scalloc( reax->num_atom_types, sizeof(char*), "tor_flag[i,j]", comm );
for (k=0; k < reax->num_atom_types; k++) {
reax->fbp[i][j][k] = (four_body_header*)
scalloc( reax->num_atom_types, sizeof(four_body_header), "fbp[i,j,k]",
comm );
tor_flag[i][j][k] = (char*)
scalloc( reax->num_atom_types, sizeof(char), "tor_flag[i,j,k]",
comm );
}
}
}
// 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] );
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;
if( reax->sbp[i].rcore2>0.01 && reax->sbp[i].acore2>0.01 ){ // Inner-wall
if( reax->sbp[i].gamma_w>0.5 ){ // Shielding vdWaals
if( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 3 )
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
"Force field parameters for element %s\n" \
"indicate inner wall+shielding, but earlier\n" \
"atoms indicate different vdWaals-method.\n" \
"This may cause division-by-zero errors.\n" \
"Keeping vdWaals-setting for earlier atoms.\n",
reax->sbp[i].name );
else{
reax->gp.vdw_type = 3;
#if defined(DEBUG)
fprintf( stderr, "vdWaals type for element %s: Shielding+inner-wall",
reax->sbp[i].name );
#endif
}
}
else { // No shielding vdWaals parameters present
if( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 2 )
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
"Force field parameters for element %s\n" \
"indicate inner wall without shielding, but earlier\n" \
"atoms indicate different vdWaals-method.\n" \
"This may cause division-by-zero errors.\n" \
"Keeping vdWaals-setting for earlier atoms.\n",
reax->sbp[i].name );
else{
reax->gp.vdw_type = 2;
#if defined(DEBUG)
fprintf( stderr,"vdWaals type for element%s: No Shielding,inner-wall",
reax->sbp[i].name );
#endif
}
}
}
else{ // No Inner wall parameters present
if( reax->sbp[i].gamma_w>0.5 ){ // Shielding vdWaals
if( reax->gp.vdw_type != 0 && reax->gp.vdw_type != 1 )
fprintf( stderr, "Warning: inconsistent vdWaals-parameters\n" \
"Force field parameters for element %s\n" \
"indicate shielding without inner wall, but earlier\n" \
"atoms indicate different vdWaals-method.\n" \
"This may cause division-by-zero errors.\n" \
"Keeping vdWaals-setting for earlier atoms.\n",
reax->sbp[i].name );
else{
reax->gp.vdw_type = 1;
#if defined(DEBUG)
fprintf( stderr,"vdWaals type for element%s: Shielding,no inner-wall",
reax->sbp[i].name );
#endif
}
}
else{
fprintf( stderr, "Error: inconsistent vdWaals-parameters\n"\
"No shielding or inner-wall set for element %s\n",
reax->sbp[i].name );
MPI_Abort( comm, INVALID_INPUT );
}
}
}
#if defined(DEBUG)
fprintf( stderr, "vdWaals type: %d\n", 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: changed valency_val to valency_boc for %s\n",
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;
if (j < reax->num_atom_types && k < reax->num_atom_types) {
val = atof(tmp[2]); reax->tbp[j][k].De_s = val;
reax->tbp[k][j].De_s = val;
val = atof(tmp[3]); reax->tbp[j][k].De_p = val;
reax->tbp[k][j].De_p = val;
val = atof(tmp[4]); reax->tbp[j][k].De_pp = val;
reax->tbp[k][j].De_pp = val;
val = atof(tmp[5]); reax->tbp[j][k].p_be1 = val;
reax->tbp[k][j].p_be1 = val;
val = atof(tmp[6]); reax->tbp[j][k].p_bo5 = val;
reax->tbp[k][j].p_bo5 = val;
val = atof(tmp[7]); reax->tbp[j][k].v13cor = val;
reax->tbp[k][j].v13cor = val;
val = atof(tmp[8]); reax->tbp[j][k].p_bo6 = val;
reax->tbp[k][j].p_bo6 = val;
val = atof(tmp[9]); reax->tbp[j][k].p_ovun1 = val;
reax->tbp[k][j].p_ovun1 = val;
/* line 2 */
fgets(s,MAX_LINE,fp);
c=Tokenize(s,&tmp);
val = atof(tmp[0]); reax->tbp[j][k].p_be2 = val;
reax->tbp[k][j].p_be2 = val;
val = atof(tmp[1]); reax->tbp[j][k].p_bo3 = val;
reax->tbp[k][j].p_bo3 = val;
val = atof(tmp[2]); reax->tbp[j][k].p_bo4 = val;
reax->tbp[k][j].p_bo4 = val;
val = atof(tmp[3]);
val = atof(tmp[4]); reax->tbp[j][k].p_bo1 = val;
reax->tbp[k][j].p_bo1 = val;
val = atof(tmp[5]); reax->tbp[j][k].p_bo2 = val;
reax->tbp[k][j].p_bo2 = val;
val = atof(tmp[6]); reax->tbp[j][k].ovc = val;
reax->tbp[k][j].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++) {
reax->tbp[i][j].r_s = 0.5 *
(reax->sbp[i].r_s + reax->sbp[j].r_s);
reax->tbp[j][i].r_s = 0.5 *
(reax->sbp[j].r_s + reax->sbp[i].r_s);
reax->tbp[i][j].r_p = 0.5 *
(reax->sbp[i].r_pi + reax->sbp[j].r_pi);
reax->tbp[j][i].r_p = 0.5 *
(reax->sbp[j].r_pi + reax->sbp[i].r_pi);
reax->tbp[i][j].r_pp = 0.5 *
(reax->sbp[i].r_pi_pi + reax->sbp[j].r_pi_pi);
reax->tbp[j][i].r_pp = 0.5 *
(reax->sbp[j].r_pi_pi + reax->sbp[i].r_pi_pi);
reax->tbp[i][j].p_boc3 =
sqrt(reax->sbp[i].b_o_132 *
reax->sbp[j].b_o_132);
reax->tbp[j][i].p_boc3 =
sqrt(reax->sbp[j].b_o_132 *
reax->sbp[i].b_o_132);
reax->tbp[i][j].p_boc4 =
sqrt(reax->sbp[i].b_o_131 *
reax->sbp[j].b_o_131);
reax->tbp[j][i].p_boc4 =
sqrt(reax->sbp[j].b_o_131 *
reax->sbp[i].b_o_131);
reax->tbp[i][j].p_boc5 =
sqrt(reax->sbp[i].b_o_133 *
reax->sbp[j].b_o_133);
reax->tbp[j][i].p_boc5 =
sqrt(reax->sbp[j].b_o_133 *
reax->sbp[i].b_o_133);
reax->tbp[i][j].D =
sqrt(reax->sbp[i].epsilon *
reax->sbp[j].epsilon);
reax->tbp[j][i].D =
sqrt(reax->sbp[j].epsilon *
reax->sbp[i].epsilon);
reax->tbp[i][j].alpha =
sqrt(reax->sbp[i].alpha *
reax->sbp[j].alpha);
reax->tbp[j][i].alpha =
sqrt(reax->sbp[j].alpha *
reax->sbp[i].alpha);
reax->tbp[i][j].r_vdW =
2.0 * sqrt(reax->sbp[i].r_vdw * reax->sbp[j].r_vdw);
reax->tbp[j][i].r_vdW =
2.0 * sqrt(reax->sbp[j].r_vdw * reax->sbp[i].r_vdw);
reax->tbp[i][j].gamma_w =
sqrt(reax->sbp[i].gamma_w *
reax->sbp[j].gamma_w);
reax->tbp[j][i].gamma_w =
sqrt(reax->sbp[j].gamma_w *
reax->sbp[i].gamma_w);
reax->tbp[i][j].gamma =
pow(reax->sbp[i].gamma *
reax->sbp[j].gamma,-1.5);
reax->tbp[j][i].gamma =
pow(reax->sbp[j].gamma *
reax->sbp[i].gamma,-1.5);
// additions for additional vdWaals interaction types - inner core
reax->tbp[i][j].rcore = reax->tbp[j][i].rcore =
sqrt( reax->sbp[i].rcore2 * reax->sbp[j].rcore2 );
reax->tbp[i][j].ecore = reax->tbp[j][i].ecore =
sqrt( reax->sbp[i].ecore2 * reax->sbp[j].ecore2 );
reax->tbp[i][j].acore = reax->tbp[j][i].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;
if (j < reax->num_atom_types && k < reax->num_atom_types) {
val = atof(tmp[2]);
if (val > 0.0) {
reax->tbp[j][k].D = val;
reax->tbp[k][j].D = val;
}
val = atof(tmp[3]);
if (val > 0.0) {
reax->tbp[j][k].r_vdW = 2 * val;
reax->tbp[k][j].r_vdW = 2 * val;
}
val = atof(tmp[4]);
if (val > 0.0) {
reax->tbp[j][k].alpha = val;
reax->tbp[k][j].alpha = val;
}
val = atof(tmp[5]);
if (val > 0.0) {
reax->tbp[j][k].r_s = val;
reax->tbp[k][j].r_s = val;
}
val = atof(tmp[6]);
if (val > 0.0) {
reax->tbp[j][k].r_p = val;
reax->tbp[k][j].r_p = val;
}
val = atof(tmp[7]);
if (val > 0.0) {
reax->tbp[j][k].r_pp = val;
reax->tbp[k][j].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][j][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;
if (j < reax->num_atom_types && k < reax->num_atom_types &&
m < reax->num_atom_types) {
cnt = reax->thbp[j][k][m].cnt;
reax->thbp[j][k][m].cnt++;
reax->thbp[m][k][j].cnt++;
val = atof(tmp[3]);
reax->thbp[j][k][m].prm[cnt].theta_00 = val;
reax->thbp[m][k][j].prm[cnt].theta_00 = val;
val = atof(tmp[4]);
reax->thbp[j][k][m].prm[cnt].p_val1 = val;
reax->thbp[m][k][j].prm[cnt].p_val1 = val;
val = atof(tmp[5]);
reax->thbp[j][k][m].prm[cnt].p_val2 = val;
reax->thbp[m][k][j].prm[cnt].p_val2 = val;
val = atof(tmp[6]);
reax->thbp[j][k][m].prm[cnt].p_coa1 = val;
reax->thbp[m][k][j].prm[cnt].p_coa1 = val;
val = atof(tmp[7]);
reax->thbp[j][k][m].prm[cnt].p_val7 = val;
reax->thbp[m][k][j].prm[cnt].p_val7 = val;
val = atof(tmp[8]);
reax->thbp[j][k][m].prm[cnt].p_pen1 = val;
reax->thbp[m][k][j].prm[cnt].p_pen1 = val;
val = atof(tmp[9]);
reax->thbp[j][k][m].prm[cnt].p_val4 = val;
reax->thbp[m][k][j].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][j][k][m].cnt = 0;
tor_flag[i][j][k][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;
if (j >= 0 && n >= 0) { // this means the entry is not in compact form
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[j][k][m][n] = 1;
tor_flag[n][m][k][j] = 1;
reax->fbp[j][k][m][n].cnt = 1;
reax->fbp[n][m][k][j].cnt = 1;
/* cnt = reax->fbp[j][k][m][n].cnt;
reax->fbp[j][k][m][n].cnt++;
reax->fbp[n][m][k][j].cnt++; */
val = atof(tmp[4]);
reax->fbp[j][k][m][n].prm[0].V1 = val;
reax->fbp[n][m][k][j].prm[0].V1 = val;
val = atof(tmp[5]);
reax->fbp[j][k][m][n].prm[0].V2 = val;
reax->fbp[n][m][k][j].prm[0].V2 = val;
val = atof(tmp[6]);
reax->fbp[j][k][m][n].prm[0].V3 = val;
reax->fbp[n][m][k][j].prm[0].V3 = val;
val = atof(tmp[7]);
reax->fbp[j][k][m][n].prm[0].p_tor1 = val;
reax->fbp[n][m][k][j].prm[0].p_tor1 = val;
val = atof(tmp[8]);
reax->fbp[j][k][m][n].prm[0].p_cot1 = val;
reax->fbp[n][m][k][j].prm[0].p_cot1 = val;
}
}
else { /* This means the entry is of the form 0-X-Y-0 */
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++ ) {
reax->fbp[p][k][m][o].cnt = 1;
reax->fbp[o][m][k][p].cnt = 1;
/* cnt = reax->fbp[p][k][m][o].cnt;
reax->fbp[p][k][m][o].cnt++;
reax->fbp[o][m][k][p].cnt++; */
if (tor_flag[p][k][m][o] == 0) {
reax->fbp[p][k][m][o].prm[0].V1 = atof(tmp[4]);
reax->fbp[p][k][m][o].prm[0].V2 = atof(tmp[5]);
reax->fbp[p][k][m][o].prm[0].V3 = atof(tmp[6]);
reax->fbp[p][k][m][o].prm[0].p_tor1 = atof(tmp[7]);
reax->fbp[p][k][m][o].prm[0].p_cot1 = atof(tmp[8]);
}
if (tor_flag[o][m][k][p] == 0) {
reax->fbp[o][m][k][p].prm[0].V1 = atof(tmp[4]);
reax->fbp[o][m][k][p].prm[0].V2 = atof(tmp[5]);
reax->fbp[o][m][k][p].prm[0].V3 = atof(tmp[6]);
reax->fbp[o][m][k][p].prm[0].p_tor1 = atof(tmp[7]);
reax->fbp[o][m][k][p].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;
if( j < reax->num_atom_types && m < reax->num_atom_types ) {
val = atof(tmp[3]);
reax->hbp[j][k][m].r0_hb = val;
val = atof(tmp[4]);
reax->hbp[j][k][m].p_hb1 = val;
val = atof(tmp[5]);
reax->hbp[j][k][m].p_hb2 = val;
val = atof(tmp[6]);
reax->hbp[j][k][m].p_hb3 = val;
}
}
/* deallocate helper storage */
for( i = 0; i < MAX_TOKENS; i++ )
free( tmp[i] );
free( tmp );
free( s );
/* deallocate tor_flag */
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++ ) {
free( tor_flag[i][j][k] );
}
free( tor_flag[i][j] );
}
free( tor_flag[i] );
}
free( tor_flag );
// close file
fclose(fp);
#if defined(DEBUG_FOCUS)
fprintf( stderr, "force field read\n" );
#endif
return SUCCESS;
}
c4science · Help