Add support for the arithmetic vdw rule

_codeql_detected_source_root

@@ -0,0 +1 @@
+.

src/Qgpu/topology.py

@@ -7,6 +7,7 @@
 # Q-GPU libraries
 import IO
 
+
 class Topology():
     def __init__(self):
         self.data = {'header'           : {},
@@ -44,6 +45,7 @@ def __init__(self):
                      'topdir'           : None,
                      'coulomb'          : None,
                      '14scaling'        : None,
+                     'vdw_rule'         : None,
                     }
 
 class Read_Topology(object):
@@ -88,6 +90,18 @@ def Q(self):
         cnt = 0
         switch = 0
 
+        # Initialize vdW parameter lists to detect missing sections
+        Aii_normal = []
+        Bii_normal = []
+        Aii_polar = []
+        Bii_polar = []
+        Aii_14 = []
+        Bii_14 = []
+        Masses = []
+
+        # Track which vdW format was detected: 'geometric' or 'arithmetic'
+        vdw_format_detected = None
+
         with open(self.top) as infile:
             for line in infile:
 
@@ -152,6 +166,7 @@ def Q(self):
                     continue
 
                 if 'vdW combination rule' in line:
+                    self.data['vdw_rule'] = line.split()[0]
                     block = 13
                     continue
 
@@ -160,38 +175,71 @@ def Q(self):
 
                 if 'Masses' in line:
                     block = 15
-                    Masses = []
                     continue
-
+
+                # Geometric vdW format (vdw_rule=1): sqrt(Aii), sqrt(Bii)
                 if 'sqrt (Aii) normal' in line:
+                    vdw_format_detected = 'geometric'
                     Aii_normal = []
                     block = 16
                     continue
-                                                                                                                   
+
                 if 'sqrt (Bii) normal' in line:
                     Bii_normal = []
                     block = 17
                     continue
-                                                                                                                   
+
                 if 'sqrt (Aii) polar' in line:
                     Aii_polar = []
                     block = 18
                     continue
-                                                                                                                   
+
                 if 'sqrt (Bii) polar' in line:
                     Bii_polar = []
                     block = 19
                     continue
-                    
+
                 if 'sqrt (Aii) 1-4' in line:
                     Aii_14 = []
                     block = 20
                     continue
-                    
+
                 if 'sqrt (Bii) 1-4' in line:
                     Bii_14 = []
                     block = 21
                     continue
+
+                # Arithmetic vdW format (vdw_rule=2): R*, epsilon
+                if 'R* normal:' in line:
+                    vdw_format_detected = 'arithmetic'
+                    Aii_normal = []
+                    block = 16
+                    continue
+
+                if 'epsilon normal:' in line:
+                    Bii_normal = []
+                    block = 17
+                    continue
+
+                if 'R* polar:' in line:
+                    Aii_polar = []
+                    block = 18
+                    continue
+
+                if 'epsilon polar:' in line:
+                    Bii_polar = []
+                    block = 19
+                    continue
+
+                if 'R* 1-4:' in line:
+                    Aii_14 = []
+                    block = 20
+                    continue
+
+                if 'epsilon 1-4:' in line:
+                    Bii_14 = []
+                    block = 21
+                    continue
 
                 if 'No. of type-2 vdW interactions' in line:
                     block = 22
@@ -446,7 +494,33 @@ def Q(self):
                             l = '1'
 
                         self.data['excluded'].append(l)
-
+
+        # exit when vdW rule was not specified
+        if self.data['vdw_rule'] is None:
+            print("FATAL: vdW combination rule not specified in topology")
+            sys.exit()
+
+        if self.data['vdw_rule'] not in ('1', '2'):
+            print("FATAL: Invalid vdW combination rule '{}' (must be 1 or 2)".format(
+                self.data['vdw_rule']))
+            sys.exit()
+
+        # validate format matches declared rule
+        if vdw_format_detected is None:
+            print("FATAL: No vdW parameter sections found in topology")
+            sys.exit()
+
+        expected_format = 'geometric' if self.data['vdw_rule'] == '1' else 'arithmetic'
+        if vdw_format_detected != expected_format:
+            print("FATAL: vdw_rule={} but found {} format section headers".format(
+                self.data['vdw_rule'], vdw_format_detected))
+            sys.exit()
+
+        # validate all vdW parameter sections were populated
+        if not Aii_normal or not Bii_normal or not Aii_polar or not Bii_polar or not Aii_14 or not Bii_14:
+            print("FATAL: Missing required vdW parameter sections in topology")
+            sys.exit()
+
         # header construct
         self.data['header'] = header
 
@@ -711,20 +785,21 @@ def CSV(self,wd):
 
         #Topo.csv
         with open(self.wd + '/topo.csv','w') as outfile:
-            outfile.write('7\n')
+            outfile.write('8\n')
             outfile.write(self.data['solvtype'] + '\n')
             outfile.write(self.data['exclusion'] + '\n')
             outfile.write(self.data['radii'] + '\n')
             outfile.write('{};{};{}\n'.format(self.data['solucenter'][0],
                                               self.data['solucenter'][1],
                                               self.data['solucenter'][2],))
-            
+
             outfile.write('{};{};{}\n'.format(self.data['solvcenter'][0],
                                               self.data['solvcenter'][1],
                                               self.data['solvcenter'][2],))
-            
+
             outfile.write(self.data['14scaling'] + '\n')
             outfile.write(self.data['coulomb'] + '\n')
+            outfile.write(self.data['vdw_rule'] + '\n')
 
         # Charge groups
         with open(self.wd + '/charge_groups.csv','w') as outfile:

src/core/cuda/src/CudaNonbondedForce.cu

@@ -3,6 +3,7 @@
 #include "cuda/include/CudaContext.cuh"
 #include "cuda/include/CudaNonbondedForce.cuh"
 #include "system.h"
+#include "vdw_rules.h"
 
 namespace CudaNonbondedForce {
 bool is_initialized = false;
@@ -45,9 +46,6 @@ __device__ __forceinline__ catype_t shfl_catype(catype_t v, int srcLane, unsigne
 }
 
 __device__ void calculate_unforce_bound(
-    // const int x_idx,
-    // const int y_idx,
-
     const coord_t& x,
     const coord_t& y,
 
@@ -64,6 +62,7 @@ __device__ void calculate_unforce_bound(
 
     const double scaling,
 
+    const int vdw_rule,
     const double lambda,
 
     double& evdw,
@@ -81,8 +80,12 @@ __device__ void calculate_unforce_bound(
 
     ecoul = scaling * coulomb_constant * x_charge * y_charge * r * lambda;
 
-    double v_a = r6 * r6 * x_aii * y_aii;
-    double v_b = r6 * x_bii * y_bii;
+    double v_a, v_b;
+    if (vdw_rule == VDW_GEOMETRIC) {
+        calc_vdw_geometric(x_aii, y_aii, x_bii, y_bii, r6, &v_a, &v_b);
+    } else {
+        calc_vdw_arithmetic(x_aii, y_aii, x_bii, y_bii, r6, &v_a, &v_b);
+    }
     v_a *= lambda;
     v_b *= lambda;
     evdw = v_a - v_b;
@@ -263,6 +266,7 @@ __global__ void calc_nonbonded_force_kernel(
                 bj_bii,
                 d_topo.coulomb_constant,
                 scaling,
+                d_topo.vdw_rule,
                 lambda,
                 evdw,
                 ecoul,

src/core/nonbonded.cu

@@ -1,5 +1,6 @@
 #include "system.h"
 #include "nonbonded.h"
+#include "vdw_rules.h"
 #include "utils.h"
 #include <stdio.h>
 
@@ -69,8 +70,11 @@ void calc_nonbonded_pp_forces() {
             ai_bii = bond14 ? ai_type.bii_1_4 : ai_type.bii_normal;
             aj_bii = bond14 ? aj_type.bii_1_4 : aj_type.bii_normal;
 
-            V_a = r6a * r6a * ai_aii * aj_aii;
-            V_b = r6a * ai_bii * aj_bii;
+            if (topo.vdw_rule == VDW_GEOMETRIC) {
+                calc_vdw_geometric(ai_aii, aj_aii, ai_bii, aj_bii, r6a, &V_a, &V_b);
+            } else {
+                calc_vdw_arithmetic(ai_aii, aj_aii, ai_bii, aj_bii, r6a, &V_a, &V_b);
+            }
             dva = r2a * ( -Vela -12 * V_a + 6 * V_b);
 
             dvelocities[i].x -= dva * da.x;
@@ -202,8 +206,11 @@ __device__ void calc_pp_dvel_matrix_incr(int row, int pi, int column, int pj,
     ai_bii = bond14 ? ai_type.bii_1_4 : ai_type.bii_normal;
     aj_bii = bond14 ? aj_type.bii_1_4 : aj_type.bii_normal;
 
-    V_a = r6a * r6a * ai_aii * aj_aii;
-    V_b = r6a * ai_bii * aj_bii;
+    if (D_topo.vdw_rule == VDW_GEOMETRIC) {
+        calc_vdw_geometric(ai_aii, aj_aii, ai_bii, aj_bii, r6a, &V_a, &V_b);
+    } else {
+        calc_vdw_arithmetic(ai_aii, aj_aii, ai_bii, aj_bii, r6a, &V_a, &V_b);
+    }
     dva = r2a * ( -Vela -12 * V_a + 6 * V_b);
 
     patom_a->x = -dva * da.x;

src/core/parse.cu

@@ -1,5 +1,6 @@
 #include "parse.h"
 #include "system.h"
+#include "vdw_rules.h"
 #include <stdio.h>
 #include <unistd.h>
 
@@ -350,6 +351,15 @@ void init_topo(const char *filename) {
 
     topo.el14_scale = strtod(file.buffer[6][0], &eptr);
     topo.coulomb_constant = strtod(file.buffer[7][0], &eptr);
+    topo.vdw_rule = atoi(file.buffer[8][0]);
+    #ifdef VERBOSE
+    printf("read %d into vdw_rule (1=geometric, 2=arithmetic)\n", topo.vdw_rule);
+    #endif
+
+    if (topo.vdw_rule != VDW_GEOMETRIC && topo.vdw_rule != VDW_ARITHMETIC) {
+        printf(">>> FATAL: Invalid vdw_rule %d. Must be 1 (geometric) or 2 (arithmetic). Exiting...\n", topo.vdw_rule);
+        exit(EXIT_FAILURE);
+    }
 
     clean_csv(file);
 }
@@ -881,6 +891,19 @@ void init_catypes(const char* filename) {
         catypes[i] = catype;
     }
 
+    // Preprocess bii parameters for arithmetic rule: convert ε to √ε
+    // This matches Fortran md.f90:14747-14752 preprocessing
+    if (topo.vdw_rule == VDW_ARITHMETIC) {
+        for (int i = 0; i < n_catypes; i++) {
+            catypes[i].bii_normal = sqrt(fabs(catypes[i].bii_normal));
+            catypes[i].bii_polar = sqrt(fabs(catypes[i].bii_polar));
+            catypes[i].bii_1_4 = sqrt(fabs(catypes[i].bii_1_4));
+        }
+        #ifdef VERBOSE
+        printf("Preprocessed catypes bii parameters for arithmetic vdW rule\n");
+        #endif
+    }
+
     clean_csv(file);
 }
 
@@ -1056,6 +1079,17 @@ void init_qcatypes(const char*filename) {
         q_catypes[i].m = strtod(file.buffer[i+1][7], &eptr);
     }
 
+    // Preprocess Bi parameters for arithmetic rule: convert ε to √ε
+    if (topo.vdw_rule == VDW_ARITHMETIC) {
+        for (int i = 0; i < n_qcatypes; i++) {
+            q_catypes[i].Bi = sqrt(fabs(q_catypes[i].Bi));
+            q_catypes[i].Bi_14 = sqrt(fabs(q_catypes[i].Bi_14));
+        }
+        #ifdef VERBOSE
+        printf("Preprocessed q_catypes Bi parameters for arithmetic vdW rule\n");
+        #endif
+    }
+
     clean_csv(file);
 }