reaxc_ffield.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Nov 1, 22:02

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;
	}

reaxc_ffield.cppNo OneTemporaryActions

File Metadata

reaxc_ffield.cppView Options

Event Timeline

reaxc_ffield.cpp
No OneTemporary
Actions

reaxc_ffield.cpp
View Options