From a27584a5710cd0772105d7b5f57ad4edd62b5e19 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@cse.msu.edu>
Date: Sat, 11 Feb 2017 13:21:10 -0500
Subject: [PATCH] Begin refactoring code to facilate other charge methods (EEM,
ACKS2, etc.).
---
environ/param.gpu.water | 19 ++--
sPuReMD/src/QEq.c | 70 +++++++-------
sPuReMD/src/allocate.c | 2 +
sPuReMD/src/control.c | 66 +++++++------
sPuReMD/src/forces.c | 210 ++++++++++++++++++++++++++++++++--------
sPuReMD/src/lin_alg.c | 26 ++---
sPuReMD/src/mytypes.h | 26 +++--
sPuReMD/src/traj.c | 2 +-
8 files changed, 283 insertions(+), 138 deletions(-)
diff --git a/environ/param.gpu.water b/environ/param.gpu.water
index 6712d0b6..0b251cf4 100644
--- a/environ/param.gpu.water
+++ b/environ/param.gpu.water
@@ -14,15 +14,16 @@ bond_graph_cutoff 0.3 ! bond strength cutoff for bond
thb_cutoff 0.001 ! cutoff value for three body interactions (Angstroms)
hbond_cutoff 7.50 ! cutoff distance for hydrogen bond interactions (Angstroms)
-qeq_solver_type 0 ! iterative linear solver used for equilibration kernel (QEq)
-qeq_solver_q_err 1e-6 ! relative residual norm threshold used in solver
-qeq_domain_sparsity 1.0 ! scalar for scaling cut-off distance, used to sparsify QEq matrix (between 0.0 and 1.0)
-pre_comp_type 1 ! method used to compute QEq preconditioner, if applicable
-pre_comp_refactor 100 ! nsteps to recompute preconditioner
-pre_comp_droptol 0.0 ! threshold tolerance for dropping values in preconditioner computation, if applicable
-pre_comp_sweeps 3 ! sweeps to compute preconditioner (ILU_PAR)
-pre_app_type 1 ! method used to apply QEq preconditioner
-pre_app_jacobi_iters 50 ! number of Jacobi iterations used for applying QEq precondition, if applicable
+charge_method 0 ! charge method: 0 = QEq, 1 = EEM, 2 = ACKS2
+cm_solver_type 0 ! iterative linear solver used in charge method
+cm_solver_q_err 1e-6 ! 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 100 ! 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.0 ! desired initial temperature of the simulated system
temp_final 300.0 ! desired final temperature of the simulated system
diff --git a/sPuReMD/src/QEq.c b/sPuReMD/src/QEq.c
index 026a3ae1..972886fd 100644
--- a/sPuReMD/src/QEq.c
+++ b/sPuReMD/src/QEq.c
@@ -1365,7 +1365,7 @@ static void Init_MatVec( const reax_system * const system, const control_params
sparse_matrix *Hptr;
// char fname[100];
- if (control->qeq_domain_sparsify_enabled)
+ if (control->cm_domain_sparsify_enabled)
{
Hptr = workspace->H_sp;
}
@@ -1378,23 +1378,23 @@ static void Init_MatVec( const reax_system * const system, const control_params
Hptr = create_test_mat( );
#endif
- if (control->pre_comp_refactor > 0 &&
- ((data->step - data->prev_steps) % control->pre_comp_refactor == 0 || workspace->L == NULL))
+ if (control->cm_solver_pre_comp_refactor > 0 &&
+ ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0 || workspace->L == NULL))
{
//Print_Linear_System( system, control, workspace, data->step );
time = Get_Time( );
- if ( control->pre_comp_type != DIAG_PC )
+ if ( control->cm_solver_pre_comp_type != DIAG_PC )
{
Sort_Matrix_Rows( workspace->H );
- if ( control->qeq_domain_sparsify_enabled == TRUE )
+ if ( control->cm_domain_sparsify_enabled == TRUE )
{
Sort_Matrix_Rows( workspace->H_sp );
}
- if ( control->pre_app_type == TRI_SOLVE_GC_PA )
+ if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
- if ( control->qeq_domain_sparsify_enabled == TRUE )
+ if ( control->cm_domain_sparsify_enabled == TRUE )
{
Hptr = setup_graph_coloring( workspace->H_sp );
}
@@ -1412,7 +1412,7 @@ static void Init_MatVec( const reax_system * const system, const control_params
fprintf( stderr, "H matrix sorted\n" );
#endif
- switch ( control->pre_comp_type )
+ switch ( control->cm_solver_pre_comp_type )
{
case DIAG_PC:
if ( workspace->Hdia_inv == NULL )
@@ -1427,7 +1427,7 @@ static void Init_MatVec( const reax_system * const system, const control_params
break;
case ICHOLT_PC:
- Calculate_Droptol( Hptr, workspace->droptol, control->pre_comp_droptol );
+ Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "drop tolerances calculated\n" );
@@ -1465,11 +1465,11 @@ static void Init_MatVec( const reax_system * const system, const control_params
}
}
- data->timing.pre_comp += ILU_PAR( Hptr, control->pre_comp_sweeps, workspace->L, workspace->U );
+ data->timing.pre_comp += ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L, workspace->U );
break;
case ILUT_PAR_PC:
- Calculate_Droptol( Hptr, workspace->droptol, control->pre_comp_droptol );
+ Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "drop tolerances calculated\n" );
#endif
@@ -1485,7 +1485,7 @@ static void Init_MatVec( const reax_system * const system, const control_params
}
}
- data->timing.pre_comp += ILUT_PAR( Hptr, workspace->droptol, control->pre_comp_sweeps,
+ data->timing.pre_comp += ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
workspace->L, workspace->U );
break;
@@ -1619,44 +1619,44 @@ void QEq( reax_system * const system, control_params * const control, simulation
// Print_Linear_System( system, control, workspace, data->step );
// }
- switch ( control->qeq_solver_type )
+ switch ( control->cm_solver_type )
{
case GMRES_S:
- iters = GMRES( workspace, control, data, workspace->H, workspace->b_s, control->qeq_solver_q_err,
+ iters = GMRES( workspace, control, data, workspace->H, workspace->b_s, control->cm_solver_q_err,
workspace->s[0], out_control->log,
- ((data->step - data->prev_steps) % control->pre_comp_refactor == 0) ? TRUE : FALSE );
- iters += GMRES( workspace, control, data, workspace->H, workspace->b_t, control->qeq_solver_q_err,
+ ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) ? TRUE : FALSE );
+ iters += GMRES( workspace, control, data, workspace->H, workspace->b_t, control->cm_solver_q_err,
workspace->t[0], out_control->log, FALSE );
break;
case GMRES_H_S:
- iters = GMRES_HouseHolder( workspace, control, data, workspace->H, workspace->b_s, control->qeq_solver_q_err,
- workspace->s[0], out_control->log, (data->step - data->prev_steps) % control->pre_comp_refactor == 0 );
- iters += GMRES_HouseHolder( workspace, control, data, workspace->H, workspace->b_t, control->qeq_solver_q_err,
+ iters = GMRES_HouseHolder( workspace, control, data, workspace->H, workspace->b_s, control->cm_solver_q_err,
+ workspace->s[0], out_control->log, (data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0 );
+ iters += GMRES_HouseHolder( workspace, control, data, workspace->H, workspace->b_t, control->cm_solver_q_err,
workspace->t[0], out_control->log, 0 );
break;
case CG_S:
- iters = CG( workspace, workspace->H, workspace->b_s, control->qeq_solver_q_err,
+ iters = CG( workspace, workspace->H, workspace->b_s, control->cm_solver_q_err,
workspace->s[0], out_control->log ) + 1;
- iters += CG( workspace, workspace->H, workspace->b_t, control->qeq_solver_q_err,
+ iters += CG( workspace, workspace->H, workspace->b_t, control->cm_solver_q_err,
workspace->t[0], out_control->log ) + 1;
-// iters = CG( workspace, workspace->H, workspace->b_s, control->qeq_solver_q_err,
-// workspace->L, workspace->U, workspace->s[0], control->pre_app_type,
-// control->pre_app_jacobi_iters, out_control->log ) + 1;
-// iters += CG( workspace, workspace->H, workspace->b_t, control->qeq_solver_q_err,
-// workspace->L, workspace->U, workspace->t[0], control->pre_app_type,
-// control->pre_app_jacobi_iters, out_control->log ) + 1;
+// iters = CG( workspace, workspace->H, workspace->b_s, control->cm_solver_q_err,
+// workspace->L, workspace->U, workspace->s[0], control->cm_solver_pre_app_type,
+// control->cm_solver_pre_app_jacobi_iters, out_control->log ) + 1;
+// iters += CG( workspace, workspace->H, workspace->b_t, control->cm_solver_q_err,
+// workspace->L, workspace->U, workspace->t[0], control->cm_solver_pre_app_type,
+// control->cm_solver_pre_app_jacobi_iters, out_control->log ) + 1;
break;
case SDM_S:
- iters = SDM( workspace, workspace->H, workspace->b_s, control->qeq_solver_q_err,
+ iters = SDM( workspace, workspace->H, workspace->b_s, control->cm_solver_q_err,
workspace->s[0], out_control->log ) + 1;
- iters += SDM( workspace, workspace->H, workspace->b_t, control->qeq_solver_q_err,
+ iters += SDM( workspace, workspace->H, workspace->b_t, control->cm_solver_q_err,
workspace->t[0], out_control->log ) + 1;
-// iters = SDM( workspace, workspace->H, workspace->b_s, control->qeq_solver_q_err,
-// workspace->L, workspace->U, workspace->s[0], control->pre_app_type,
-// control->pre_app_jacobi_iters, out_control->log ) + 1;
-// iters += SDM( workspace, workspace->H, workspace->b_t, control->qeq_solver_q_err,
-// workspace->L, workspace->U, workspace->t[0], control->pre_app_type,
-// control->pre_app_jacobi_iters, out_control->log ) + 1;
+// iters = SDM( workspace, workspace->H, workspace->b_s, control->cm_solver_q_err,
+// workspace->L, workspace->U, workspace->s[0], control->cm_solver_pre_app_type,
+// control->cm_solver_pre_app_jacobi_iters, out_control->log ) + 1;
+// iters += SDM( workspace, workspace->H, workspace->b_t, control->cm_solver_q_err,
+// workspace->L, workspace->U, workspace->t[0], control->cm_solver_pre_app_type,
+// control->cm_solver_pre_app_jacobi_iters, out_control->log ) + 1;
break;
default:
fprintf( stderr, "Unrecognized QEq solver selection. Terminating...\n" );
diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index c9c5321c..c8bdddb0 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -73,6 +73,7 @@ void Reallocate_Neighbor_List( list *far_nbrs, int n, int num_intrs )
}
+/* dynamic allocation of memory for matrix in CSR format */
int Allocate_Matrix( sparse_matrix **pH, int n, int m )
{
sparse_matrix *H;
@@ -97,6 +98,7 @@ int Allocate_Matrix( sparse_matrix **pH, int n, int m )
}
+/* deallocate memory for matrix in CSR format */
void Deallocate_Matrix( sparse_matrix *H )
{
free(H->start);
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 41f74496..a10ee4f2 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -70,16 +70,17 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->tabulate = 0;
- control->qeq_solver_type = GMRES_S;
- control->qeq_solver_q_err = 0.000001;
- control->qeq_domain_sparsify_enabled = FALSE;
- control->qeq_domain_sparsity = 1.0;
- control->pre_comp_type = ICHOLT_PC;
- control->pre_comp_sweeps = 3;
- control->pre_comp_refactor = 100;
- control->pre_comp_droptol = 0.01;
- control->pre_app_type = TRI_SOLVE_PA;
- control->pre_app_jacobi_iters = 50;
+ control->charge_method = QEQ_CM;
+ control->cm_solver_type = GMRES_S;
+ control->cm_solver_q_err = 0.000001;
+ control->cm_domain_sparsify_enabled = FALSE;
+ control->cm_domain_sparsity = 1.0;
+ control->cm_solver_pre_comp_type = ICHOLT_PC;
+ control->cm_solver_pre_comp_sweeps = 3;
+ control->cm_solver_pre_comp_refactor = 100;
+ 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;
control->T_init = 0.;
control->T_final = 300.;
@@ -240,51 +241,56 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
val = atof( tmp[1] );
control->hb_cut = val;
}
- else if ( strcmp(tmp[0], "qeq_solver_type") == 0 )
+ else if ( strcmp(tmp[0], "charge_method") == 0 )
{
ival = atoi( tmp[1] );
- control->qeq_solver_type = ival;
+ control->charge_method = ival;
}
- else if ( strcmp(tmp[0], "qeq_solver_q_err") == 0 )
+ else if ( strcmp(tmp[0], "cm_solver_type") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_type = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_q_err") == 0 )
{
val = atof( tmp[1] );
- control->qeq_solver_q_err = val;
+ control->cm_solver_q_err = val;
}
- else if ( strcmp(tmp[0], "qeq_domain_sparsity") == 0 )
+ else if ( strcmp(tmp[0], "cm_domain_sparsity") == 0 )
{
val = atof( tmp[1] );
- control->qeq_domain_sparsity = val;
- control->qeq_domain_sparsify_enabled = TRUE;
+ control->cm_domain_sparsity = val;
+ control->cm_domain_sparsify_enabled = TRUE;
}
- else if ( strcmp(tmp[0], "pre_comp_type") == 0 )
+ else if ( strcmp(tmp[0], "cm_solver_pre_comp_type") == 0 )
{
ival = atoi( tmp[1] );
- control->pre_comp_type = ival;
+ control->cm_solver_pre_comp_type = ival;
}
- else if ( strcmp(tmp[0], "pre_comp_refactor") == 0 )
+ else if ( strcmp(tmp[0], "cm_solver_pre_comp_refactor") == 0 )
{
ival = atoi( tmp[1] );
- control->pre_comp_refactor = ival;
+ control->cm_solver_pre_comp_refactor = ival;
}
- else if ( strcmp(tmp[0], "pre_comp_droptol") == 0 )
+ else if ( strcmp(tmp[0], "cm_solver_pre_comp_droptol") == 0 )
{
val = atof( tmp[1] );
- control->pre_comp_droptol = val;
+ control->cm_solver_pre_comp_droptol = val;
}
- else if ( strcmp(tmp[0], "pre_comp_sweeps") == 0 )
+ else if ( strcmp(tmp[0], "cm_solver_pre_comp_sweeps") == 0 )
{
ival = atoi( tmp[1] );
- control->pre_comp_sweeps = ival;
+ control->cm_solver_pre_comp_sweeps = ival;
}
- else if ( strcmp(tmp[0], "pre_app_type") == 0 )
+ else if ( strcmp(tmp[0], "cm_solver_pre_app_type") == 0 )
{
ival = atoi( tmp[1] );
- control->pre_app_type = ival;
+ control->cm_solver_pre_app_type = ival;
}
- else if ( strcmp(tmp[0], "pre_app_jacobi_iters") == 0 )
+ else if ( strcmp(tmp[0], "cm_solver_pre_app_jacobi_iters") == 0 )
{
ival = atoi( tmp[1] );
- control->pre_app_jacobi_iters = ival;
+ control->cm_solver_pre_app_jacobi_iters = ival;
}
else if ( strcmp(tmp[0], "temp_init") == 0 )
{
@@ -530,7 +536,7 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->bo_cut = 0.01 * system->reaxprm.gp.l[29];
control->r_low = system->reaxprm.gp.l[11];
control->r_cut = system->reaxprm.gp.l[12];
- control->r_sp_cut = control->r_cut * control->qeq_domain_sparsity;
+ control->r_sp_cut = control->r_cut * control->cm_domain_sparsity;
control->vlist_cut += control->r_cut;
system->g.cell_size = control->vlist_cut / 2.;
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 9108a8dd..8c21776a 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -33,6 +33,13 @@
#include "vector.h"
+typedef enum
+{
+ DIAGONAL = 0,
+ OFF_DIAGONAL = 1,
+} MATRIX_ENTRY_POSITION;
+
+
void Dummy_Interaction( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
list **lists, output_controls *out_control )
@@ -263,6 +270,159 @@ void Validate_Lists( static_storage *workspace, list **lists, int step, int n,
}
+static inline real Init_Charge_Matrix_Entry_Tab( reax_system *system,
+ control_params *control, int i, int j,
+ real r_ij, MATRIX_ENTRY_POSITION pos )
+{
+ int r;
+ real base, dif, val, ret = 0.0;
+ LR_lookup_table *t;
+
+ switch ( control->charge_method )
+ {
+ case QEQ_CM:
+ switch ( pos )
+ {
+ case OFF_DIAGONAL:
+ t = &( LR
+ [MIN( system->atoms[i].type, system->atoms[j].type )]
+ [MAX( system->atoms[i].type, system->atoms[j].type )] );
+
+ /* cubic spline interpolation */
+ r = (int)(r_ij * t->inv_dx);
+ if ( r == 0 ) ++r;
+ base = (real)(r + 1) * t->dx;
+ dif = r_ij - base;
+ val = ((t->ele[r].d * dif + t->ele[r].c) * dif + t->ele[r].b) * dif +
+ t->ele[r].a;
+ val *= EV_to_KCALpMOL / C_ele;
+
+ ret = ((i == j) ? 0.5 : 1.0) * val;
+ break;
+ case DIAGONAL:
+ ret = system->reaxprm.sbp[system->atoms[i].type].eta;
+ break;
+ default:
+ fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+ break;
+
+ case EEM_CM:
+ switch ( pos )
+ {
+ case OFF_DIAGONAL:
+ break;
+ case DIAGONAL:
+ break;
+ default:
+ fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+ break;
+
+ case ACKS2_CM:
+ //TODO
+ switch ( pos )
+ {
+ case OFF_DIAGONAL:
+ break;
+ case DIAGONAL:
+ break;
+ default:
+ fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+ break;
+
+ default:
+ fprintf( stderr, "Invalid charge method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+
+ return ret;
+}
+
+
+static inline real Init_Charge_Matrix_Entry( reax_system *system,
+ control_params *control, int i, int j,
+ real r_ij, MATRIX_ENTRY_POSITION pos )
+{
+ real Tap, dr3gamij_1, dr3gamij_3, ret = 0.0;
+
+ switch ( control->charge_method )
+ {
+ case QEQ_CM:
+ switch ( pos )
+ {
+ case OFF_DIAGONAL:
+ Tap = control->Tap7 * r_ij + control->Tap6;
+ Tap = Tap * r_ij + control->Tap5;
+ Tap = Tap * r_ij + control->Tap4;
+ Tap = Tap * r_ij + control->Tap3;
+ Tap = Tap * r_ij + control->Tap2;
+ Tap = Tap * r_ij + control->Tap1;
+ Tap = Tap * r_ij + control->Tap0;
+
+ dr3gamij_1 = ( r_ij * r_ij * r_ij +
+ system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma );
+ dr3gamij_3 = POW( dr3gamij_1 , 0.33333333333333 );
+
+ ret = ((i == j) ? 0.5 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
+ break;
+ case DIAGONAL:
+ ret = system->reaxprm.sbp[system->atoms[i].type].eta;
+ break;
+ default:
+ fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+ break;
+
+ case EEM_CM:
+ switch ( pos )
+ {
+ case OFF_DIAGONAL:
+ break;
+ case DIAGONAL:
+ break;
+ default:
+ fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+ break;
+
+ case ACKS2_CM:
+ //TODO
+ switch ( pos )
+ {
+ case OFF_DIAGONAL:
+ break;
+ case DIAGONAL:
+ break;
+ default:
+ fprintf( stderr, "[Init_forces] Invalid matrix position. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+ break;
+
+ default:
+ fprintf( stderr, "Invalid charge method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
+
+ return ret;
+}
+
+
void Init_Forces( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
list **lists, output_controls *out_control )
@@ -273,9 +433,7 @@ void Init_Forces( reax_system *system, control_params *control,
int Htop, H_sp_top, btop_i, btop_j, num_bonds, num_hbonds;
int ihb, jhb, ihb_top, jhb_top;
int flag, flag_sp;
- real r_ij, r2, self_coef;
- real dr3gamij_1, dr3gamij_3, Tap;
- //real val, dif, base;
+ real r_ij, r2;
real C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
@@ -285,7 +443,6 @@ void Init_Forces( reax_system *system, control_params *control,
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
- //LR_lookup_table *t;
reax_atom *atom_i, *atom_j;
bond_data *ibond, *jbond;
bond_order_data *bo_ij, *bo_ji;
@@ -359,22 +516,10 @@ void Init_Forces( reax_system *system, control_params *control,
r_ij = nbr_pj->d;
sbp_j = &(system->reaxprm.sbp[type_j]);
twbp = &(system->reaxprm.tbp[type_i][type_j]);
- self_coef = (i == j) ? 0.5 : 1.0;
-
- /* H matrix entry */
- Tap = control->Tap7 * r_ij + control->Tap6;
- Tap = Tap * r_ij + control->Tap5;
- Tap = Tap * r_ij + control->Tap4;
- Tap = Tap * r_ij + control->Tap3;
- Tap = Tap * r_ij + control->Tap2;
- Tap = Tap * r_ij + control->Tap1;
- Tap = Tap * r_ij + control->Tap0;
-
- dr3gamij_1 = ( r_ij * r_ij * r_ij + twbp->gamma );
- dr3gamij_3 = POW( dr3gamij_1 , 0.33333333333333 );
H->j[Htop] = j;
- H->val[Htop] = self_coef * Tap * EV_to_KCALpMOL / dr3gamij_3;
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, j,
+ r_ij, OFF_DIAGONAL );
++Htop;
/* H_sp matrix entry */
@@ -541,10 +686,10 @@ void Init_Forces( reax_system *system, control_params *control,
/* diagonal entry */
H->j[Htop] = i;
- H->val[Htop] = system->reaxprm.sbp[type_i].eta;
+ H->val[Htop] = Init_Charge_Matrix_Entry( system, control, i, j,
+ r_ij, DIAGONAL );
++Htop;
- /* diagonal entry */
H_sp->j[H_sp_top] = i;
H_sp->val[H_sp_top] = H->val[Htop - 1];
++H_sp_top;
@@ -582,11 +727,9 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
int start_i, end_i;
int type_i, type_j;
int Htop, H_sp_top, btop_i, btop_j, num_bonds, num_hbonds;
- int tmin, tmax, r;
int ihb, jhb, ihb_top, jhb_top;
int flag, flag_sp;
- real r_ij, r2, self_coef;
- real val, dif, base;
+ real r_ij, r2;
real C12, C34, C56;
real Cln_BOp_s, Cln_BOp_pi, Cln_BOp_pi2;
real BO, BO_s, BO_pi, BO_pi2;
@@ -596,7 +739,6 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
single_body_parameters *sbp_i, *sbp_j;
two_body_parameters *twbp;
far_neighbor_data *nbr_pj;
- LR_lookup_table *t;
reax_atom *atom_i, *atom_j;
bond_data *ibond, *jbond;
bond_order_data *bo_ij, *bo_ji;
@@ -670,22 +812,10 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
r_ij = nbr_pj->d;
sbp_j = &(system->reaxprm.sbp[type_j]);
twbp = &(system->reaxprm.tbp[type_i][type_j]);
- self_coef = (i == j) ? 0.5 : 1.0;
- tmin = MIN( type_i, type_j );
- tmax = MAX( type_i, type_j );
- t = &( LR[tmin][tmax] );
-
- /* cubic spline interpolation */
- r = (int)(r_ij * t->inv_dx);
- if ( r == 0 ) ++r;
- base = (real)(r + 1) * t->dx;
- dif = r_ij - base;
- val = ((t->ele[r].d * dif + t->ele[r].c) * dif + t->ele[r].b) * dif +
- t->ele[r].a;
- val *= EV_to_KCALpMOL / C_ele;
H->j[Htop] = j;
- H->val[Htop] = self_coef * val;
+ H->val[Htop] = Init_Charge_Matrix_Entry_Tab( system, control, i, j,
+ r_ij, OFF_DIAGONAL );
++Htop;
/* H_sp matrix entry */
@@ -824,10 +954,10 @@ void Init_Forces_Tab( reax_system *system, control_params *control,
/* diagonal entry */
H->j[Htop] = i;
- H->val[Htop] = system->reaxprm.sbp[type_i].eta;
+ H->val[Htop] = Init_Charge_Matrix_Entry_Tab( system, control, i, j,
+ r_ij, DIAGONAL );
++Htop;
- /* diagonal entry */
H_sp->j[H_sp_top] = i;
H_sp->val[H_sp_top] = H->val[Htop - 1];
++H_sp_top;
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 1759b41c..4d37ae6f 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -1030,12 +1030,12 @@ static void apply_preconditioner( const static_storage * const workspace, const
{
int i, si;
- switch ( control->pre_app_type )
+ switch ( control->cm_solver_pre_app_type )
{
case NONE_PA:
break;
case TRI_SOLVE_PA:
- switch ( control->pre_comp_type )
+ switch ( control->cm_solver_pre_comp_type )
{
case DIAG_PC:
diag_pre_app( workspace->Hdia_inv, y, x, workspace->H->n );
@@ -1053,7 +1053,7 @@ static void apply_preconditioner( const static_storage * const workspace, const
}
break;
case TRI_SOLVE_LEVEL_SCHED_PA:
- switch ( control->pre_comp_type )
+ switch ( control->cm_solver_pre_comp_type )
{
case DIAG_PC:
diag_pre_app( workspace->Hdia_inv, y, x, workspace->H->n );
@@ -1071,7 +1071,7 @@ static void apply_preconditioner( const static_storage * const workspace, const
}
break;
case TRI_SOLVE_GC_PA:
- switch ( control->pre_comp_type )
+ switch ( control->cm_solver_pre_comp_type )
{
case DIAG_PC:
fprintf( stderr, "Unsupported preconditioner computation/application method combination. Terminating...\n" );
@@ -1099,7 +1099,7 @@ static void apply_preconditioner( const static_storage * const workspace, const
}
break;
case JACOBI_ITER_PA:
- switch ( control->pre_comp_type )
+ switch ( control->cm_solver_pre_comp_type )
{
case DIAG_PC:
fprintf( stderr, "Unsupported preconditioner computation/application method combination. Terminating...\n" );
@@ -1133,7 +1133,7 @@ static void apply_preconditioner( const static_storage * const workspace, const
}
}
- jacobi_iter( workspace->L, Dinv_L, y, x, LOWER, control->pre_app_jacobi_iters );
+ jacobi_iter( workspace->L, Dinv_L, y, x, LOWER, control->cm_solver_pre_app_jacobi_iters );
#pragma omp master
{
@@ -1160,7 +1160,7 @@ static void apply_preconditioner( const static_storage * const workspace, const
}
}
- jacobi_iter( workspace->U, Dinv_U, y, x, UPPER, control->pre_app_jacobi_iters );
+ jacobi_iter( workspace->U, Dinv_U, y, x, UPPER, control->cm_solver_pre_app_jacobi_iters );
break;
default:
fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
@@ -1203,7 +1203,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
data->timing.solver_vector_ops += Get_Timing_Info( time_start );
}
- if ( control->pre_comp_type == DIAG_PC )
+ if ( control->cm_solver_pre_comp_type == DIAG_PC )
{
/* apply preconditioner to RHS */
#pragma omp master
@@ -1231,7 +1231,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
data->timing.solver_spmv += Get_Timing_Info( time_start );
}
- if ( control->pre_comp_type == DIAG_PC )
+ if ( control->cm_solver_pre_comp_type == DIAG_PC )
{
#pragma omp master
{
@@ -1244,7 +1244,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
}
}
- if ( control->pre_comp_type == DIAG_PC )
+ if ( control->cm_solver_pre_comp_type == DIAG_PC )
{
#pragma omp master
{
@@ -1269,7 +1269,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
}
}
- if ( control->pre_comp_type != DIAG_PC )
+ if ( control->cm_solver_pre_comp_type != DIAG_PC )
{
#pragma omp master
{
@@ -1322,7 +1322,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
data->timing.pre_app += Get_Timing_Info( time_start );
}
- if ( control->pre_comp_type == DIAG_PC )
+ if ( control->cm_solver_pre_comp_type == DIAG_PC )
{
/* apply modified Gram-Schmidt to orthogonalize the new residual */
#pragma omp master
@@ -1389,7 +1389,7 @@ int GMRES( const static_storage * const workspace, const control_params * const
#pragma omp master
{
time_start = Get_Time( );
- if ( control->pre_comp_type == DIAG_PC )
+ if ( control->cm_solver_pre_comp_type == DIAG_PC )
{
/* Givens rotations on the upper-Hessenberg matrix to make it U */
for ( i = 0; i <= j; i++ )
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 69441bb7..23653bb3 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -194,6 +194,11 @@ enum geo_formats
CUSTOM = 0, PDB = 1, BGF = 2, ASCII_RESTART = 3, BINARY_RESTART = 4, GF_N = 5,
};
+enum charge_method
+{
+ QEQ_CM = 0, EEM_CM = 1, ACKS2_CM = 2,
+};
+
enum solver
{
GMRES_S = 0, GMRES_H_S = 1, CG_S = 2, SDM_S = 3,
@@ -514,16 +519,17 @@ typedef struct
int freq_diffusion_coef;
int restrict_type;
- unsigned int qeq_solver_type;
- real qeq_solver_q_err;
- real qeq_domain_sparsity;
- unsigned int qeq_domain_sparsify_enabled;
- unsigned int pre_comp_type;
- unsigned int pre_comp_refactor;
- real pre_comp_droptol;
- unsigned int pre_comp_sweeps;
- unsigned int pre_app_type;
- unsigned int pre_app_jacobi_iters;
+ unsigned int charge_method;
+ unsigned int cm_solver_type;
+ real cm_solver_q_err;
+ real cm_domain_sparsity;
+ unsigned int cm_domain_sparsify_enabled;
+ unsigned int cm_solver_pre_comp_type;
+ unsigned int cm_solver_pre_comp_refactor;
+ real cm_solver_pre_comp_droptol;
+ unsigned int cm_solver_pre_comp_sweeps;
+ unsigned int cm_solver_pre_app_type;
+ unsigned int cm_solver_pre_app_jacobi_iters;
int molec_anal;
int freq_molec_anal;
diff --git a/sPuReMD/src/traj.c b/sPuReMD/src/traj.c
index 81207a00..f679186d 100644
--- a/sPuReMD/src/traj.c
+++ b/sPuReMD/src/traj.c
@@ -53,7 +53,7 @@ int Write_Custom_Header(reax_system *system, control_params *control,
control->bo_cut,
control->thb_cut,
control->hb_cut,
- control->qeq_solver_q_err,
+ control->cm_solver_q_err,
control->T_init,
control->T_final,
control->Tau_T,
--
GitLab