sPuReMD: corrections to ACKS2 matrix to match reference Fortran ReaxFF code.

sPuReMD/src/ffield.c

@@ -426,8 +426,8 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
             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, -3.0 );
-            reax->tbp[j][i].gamma = POW( reax->sbp[j].gamma * reax->sbp[i].gamma, -3.0 );
+            reax->tbp[i][j].gamma = SQRT( reax->sbp[i].gamma * reax->sbp[j].gamma );
+            reax->tbp[j][i].gamma = SQRT( reax->sbp[j].gamma * reax->sbp[i].gamma );
 
             reax->tbp[i][j].acore = SQRT( reax->sbp[i].acore2 * reax->sbp[j].acore2 );
             reax->tbp[j][i].acore = SQRT( reax->sbp[j].acore2 * reax->sbp[i].acore2 );

sPuReMD/src/forces.c

@@ -432,8 +432,8 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
                 Tap = Tap * r_ij + workspace->Tap[0];
 
                 /* shielding */
-                dr3gamij_1 = ( r_ij * r_ij * r_ij
-                        + system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma );
+                dr3gamij_1 = r_ij * r_ij * r_ij
+                        + POW( system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma, -3.0 );
                 dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
 
                 /* i == j: non-periodic self-interaction term
@@ -456,7 +456,8 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
         switch ( pos )
         {
             case OFF_DIAGONAL:
-                Tap = workspace->Tap[7] * r_ij + workspace->Tap[6];
+                Tap = workspace->Tap[7] * r_ij
+                    + workspace->Tap[6];
                 Tap = Tap * r_ij + workspace->Tap[5];
                 Tap = Tap * r_ij + workspace->Tap[4];
                 Tap = Tap * r_ij + workspace->Tap[3];
@@ -465,18 +466,17 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
                 Tap = Tap * r_ij + workspace->Tap[0];
 
                 /* shielding */
-                dr3gamij_1 = ( r_ij * r_ij * r_ij
-                        + system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma );
+                dr3gamij_1 = r_ij * r_ij * r_ij
+                        + POW( system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma, -3.0 );
                 dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
 
                 /* i == j: non-periodic self-interaction term
                  * i != j: periodic self-interaction term */
-                ret = ((i == j) ? 0.5 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
+                ret = ((i == j) ? 0.0 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
             break;
 
             case DIAGONAL:
-                /* EE parameters in electron-volts */
-//                ret = 2.0 * system->reaxprm.sbp[system->atoms[i].type].eta;
+                /* parameters in electron-volts */
                 ret = system->reaxprm.sbp[system->atoms[i].type].eta;
             break;