diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root new file mode 120000 index 00000000..945c9b46 --- /dev/null +++ b/_codeql_detected_source_root @@ -0,0 +1 @@ +. \ No newline at end of file diff --git a/src/Qgpu/topology.py b/src/Qgpu/topology.py index 50a85b74..5a91ead3 100644 --- a/src/Qgpu/topology.py +++ b/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: diff --git a/src/core/cuda/src/CudaNonbondedForce.cu b/src/core/cuda/src/CudaNonbondedForce.cu index 0966fb15..85f8690b 100644 --- a/src/core/cuda/src/CudaNonbondedForce.cu +++ b/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, diff --git a/src/core/nonbonded.cu b/src/core/nonbonded.cu index 9d49dab7..094b699c 100644 --- a/src/core/nonbonded.cu +++ b/src/core/nonbonded.cu @@ -1,5 +1,6 @@ #include "system.h" #include "nonbonded.h" +#include "vdw_rules.h" #include "utils.h" #include @@ -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; diff --git a/src/core/parse.cu b/src/core/parse.cu index a20dab3e..0fd6d3e1 100644 --- a/src/core/parse.cu +++ b/src/core/parse.cu @@ -1,5 +1,6 @@ #include "parse.h" #include "system.h" +#include "vdw_rules.h" #include #include @@ -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); } diff --git a/src/core/qatoms.cu b/src/core/qatoms.cu index b4a0307e..1bb59750 100644 --- a/src/core/qatoms.cu +++ b/src/core/qatoms.cu @@ -1,6 +1,7 @@ // TODO: Add impropers, bond pairs #include "system.h" #include "qatoms.h" +#include "vdw_rules.h" #include "utils.h" #include #include @@ -57,6 +58,7 @@ void calc_nonbonded_qp_forces() { r6 = r2 * r2 * r2; r2 = 1 / r2; r = sqrt(r2); + double r6inv = r2 * r2 * r2; // 1/r^6 for vdW calculation for (int state = 0; state < n_lambdas; state++) { qi_type = q_catypes[q_atypes[qi + n_qatoms * state].code - 1]; @@ -68,8 +70,11 @@ void calc_nonbonded_qp_forces() { aj_bii = bond14 ? aj_type.bii_1_4 : aj_type.bii_normal; Vel = topo.coulomb_constant * scaling * q_charges[qi + n_qatoms * state].q * ccharges[charges[j].code - 1].charge * r; - V_a = ai_aii * aj_aii / (r6 * r6); - V_b = ai_bii * aj_bii / r6; + if (topo.vdw_rule == VDW_GEOMETRIC) { + calc_vdw_geometric(ai_aii, aj_aii, ai_bii, aj_bii, r6inv, &V_a, &V_b); + } else { + calc_vdw_arithmetic(ai_aii, aj_aii, ai_bii, aj_bii, r6inv, &V_a, &V_b); + } dv = r2 * (-Vel - (12 * V_a - 6 * V_b)) * lambdas[state]; // Update forces @@ -296,6 +301,7 @@ void calc_nonbonded_qw_forces() { r6O = r2O * r2O * r2O; r2O = 1.0 / r2O; rO = sqrt(r2O); + double r6Oinv = r2O * r2O * r2O; // 1/r^6 for vdW calculation r2H1 = rH1 * rH1; r2H2 = rH2 * rH2; @@ -310,8 +316,11 @@ void calc_nonbonded_qw_forces() { ai_aii = qi_type.Ai; ai_bii = qi_type.Bi; - V_a = ai_aii * A_O / (r6O * r6O); - V_b = ai_bii * B_O / (r6O); + if (topo.vdw_rule == VDW_GEOMETRIC) { + calc_vdw_geometric(ai_aii, A_O, ai_bii, B_O, r6Oinv, &V_a, &V_b); + } else { + calc_vdw_arithmetic(ai_aii, A_O, ai_bii, B_O, r6Oinv, &V_a, &V_b); + } VelO = topo.coulomb_constant * crg_ow * q_charges[qi + n_qatoms * state].q * rO; VelH1 = topo.coulomb_constant * crg_hw * q_charges[qi + n_qatoms * state].q * rH1; @@ -493,8 +502,11 @@ void calc_nonbonded_qq_forces() { ai_bii = bond14 ? qi_type.Bi_14 : qi_type.Bi; aj_bii = bond14 ? qj_type.Bi_14 : qj_type.Bi; - 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) * lambdas[state]; dvelocities[ai].x -= dva * da.x; @@ -760,6 +772,7 @@ __device__ void calc_qw_dvel_matrix_incr(int row, int qi, int column, int n_lamb r6O = r2O * r2O * r2O; r2O = 1.0 / r2O; rO = sqrt(r2O); + double r6Oinv = r2O * r2O * r2O; // 1/r^6 for vdW calculation r2H1 = rH1 * rH1; r2H2 = rH2 * rH2; @@ -775,8 +788,11 @@ __device__ void calc_qw_dvel_matrix_incr(int row, int qi, int column, int n_lamb ai_aii = qi_type.Ai; ai_bii = qi_type.Bi; - V_a = ai_aii * A_O / (r6O * r6O); - V_b = ai_bii * B_O / (r6O); + if (D_topo.vdw_rule == VDW_GEOMETRIC) { + calc_vdw_geometric(ai_aii, A_O, ai_bii, B_O, r6Oinv, &V_a, &V_b); + } else { + calc_vdw_arithmetic(ai_aii, A_O, ai_bii, B_O, r6Oinv, &V_a, &V_b); + } VelO = D_topo.coulomb_constant * crg_ow * D_qcharges[qi + n_qatoms * state].q * rO; VelH1 = D_topo.coulomb_constant * crg_hw * D_qcharges[qi + n_qatoms * state].q * rH1; @@ -988,6 +1004,7 @@ __device__ void calc_qp_dvel_matrix_incr(int row, int qi, int column, int pj, in r6 = r2 * r2 * r2; r2 = 1 / r2; r = sqrt(r2); + double r6inv = r2 * r2 * r2; // 1/r^6 for vdW calculation for (int state = 0; state < n_lambdas; state++) { qi_type = D_qcatypes[D_qatypes[qi + n_qatoms * state].code - 1]; @@ -999,8 +1016,11 @@ __device__ void calc_qp_dvel_matrix_incr(int row, int qi, int column, int pj, in aj_bii = bond14 ? aj_type.bii_1_4 : aj_type.bii_normal; Vel = D_topo.coulomb_constant * scaling * D_qcharges[qi + n_qatoms * state].q * D_ccharges[D_charges[j].code - 1].charge * r; - V_a = ai_aii * aj_aii / (r6 * r6); - V_b = ai_bii * aj_bii / r6; + if (D_topo.vdw_rule == VDW_GEOMETRIC) { + calc_vdw_geometric(ai_aii, aj_aii, ai_bii, aj_bii, r6inv, &V_a, &V_b); + } else { + calc_vdw_arithmetic(ai_aii, aj_aii, ai_bii, aj_bii, r6inv, &V_a, &V_b); + } dv = r2 * (-Vel - (12 * V_a - 6 * V_b)) * D_lambdas[state]; // if (state == 0 && qi == 0 && pj == 1) { diff --git a/src/core/solvent.cu b/src/core/solvent.cu index 180152bc..c3b5c708 100644 --- a/src/core/solvent.cu +++ b/src/core/solvent.cu @@ -1,5 +1,6 @@ #include "system.h" #include "solvent.h" +#include "vdw_rules.h" #include "utils.h" #include @@ -40,8 +41,19 @@ void calc_nonbonded_ww_forces() { ccharge_ow = ccharges[charges[n_atoms_solute].code - 1]; ccharge_hw = ccharges[charges[n_atoms_solute+1].code - 1]; - A_OO = pow(catype_ow.aii_normal, 2); - B_OO = pow(catype_ow.bii_normal, 2); + if (topo.vdw_rule == VDW_GEOMETRIC) { + // Geometric: A_OO = sqrt(A) * sqrt(A) = A + A_OO = pow(catype_ow.aii_normal, 2); + B_OO = pow(catype_ow.bii_normal, 2); + } else { + // Arithmetic: R*_OO = R*_O + R*_O = 2*R*_O + double Rstar_OO = 2.0 * catype_ow.aii_normal; + double eps_OO = pow(catype_ow.bii_normal, 2); // sqrt(eps) * sqrt(eps) = eps + double R6_OO = pow(Rstar_OO, 6); + // Pre-compute coefficients for V_a = A_OO * r6^2 and V_b = B_OO * r6 + A_OO = eps_OO * R6_OO * R6_OO; + B_OO = 2.0 * eps_OO * R6_OO; + } crg_ow = ccharge_ow.charge; crg_hw = ccharge_hw.charge; @@ -538,8 +550,19 @@ void calc_nonbonded_ww_forces_host() { ccharge_ow = ccharges[charges[n_atoms_solute].code - 1]; ccharge_hw = ccharges[charges[n_atoms_solute + 1].code - 1]; - A_OO = pow(catype_ow.aii_normal, 2); - B_OO = pow(catype_ow.bii_normal, 2); + if (topo.vdw_rule == VDW_GEOMETRIC) { + // Geometric: A_OO = sqrt(A) * sqrt(A) = A + A_OO = pow(catype_ow.aii_normal, 2); + B_OO = pow(catype_ow.bii_normal, 2); + } else { + // Arithmetic: R*_OO = R*_O + R*_O = 2*R*_O + double Rstar_OO = 2.0 * catype_ow.aii_normal; + double eps_OO = pow(catype_ow.bii_normal, 2); // sqrt(eps) * sqrt(eps) = eps + double R6_OO = pow(Rstar_OO, 6); + // Pre-compute coefficients for V_a = A_OO * r6^2 and V_b = B_OO * r6 + A_OO = eps_OO * R6_OO * R6_OO; + B_OO = 2.0 * eps_OO * R6_OO; + } crg_ow = ccharge_ow.charge; crg_hw = ccharge_hw.charge; @@ -624,8 +647,11 @@ void calc_nonbonded_pw_forces() { ai_bii = ai_type.bii_normal; aj_bii = 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; @@ -1190,8 +1216,11 @@ __device__ void calc_pw_dvel_matrix_incr(int row, int pi, int column, int j, int ai_bii = ai_type.bii_normal; aj_bii = 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); diff --git a/src/core/system.h b/src/core/system.h index 6be0f9a8..67895107 100644 --- a/src/core/system.h +++ b/src/core/system.h @@ -215,6 +215,7 @@ struct topo_t { coord_t solvent_center; double el14_scale; double coulomb_constant; + int vdw_rule; // 1=geometric, 2=arithmetic }; struct cgrp_t { diff --git a/src/core/vdw_rules.h b/src/core/vdw_rules.h new file mode 100644 index 00000000..74776517 --- /dev/null +++ b/src/core/vdw_rules.h @@ -0,0 +1,52 @@ +// Defines vdW combination rules (geometric and arithmetic) and provides +// inline functions for computing vdW energy components for both rules. +#ifndef __VDW_RULES_H__ +#define __VDW_RULES_H__ + +#include + +// Handle CUDA keywords for non-CUDA compilation +#ifndef __CUDACC__ +#define __device__ +#define __host__ +#endif + +#define VDW_GEOMETRIC 1 +#define VDW_ARITHMETIC 2 + +// Geometric rule: A_ij = sqrt(A_i) * sqrt(A_j), B_ij = sqrt(B_i) * sqrt(B_j) +// Energy: V = A_ij * r^-12 - B_ij * r^-6 +// Parameters: ai_aii, aj_aii are sqrt(A_i), sqrt(A_j) +// ai_bii, aj_bii are sqrt(B_i), sqrt(B_j) +// r6 is 1/r^6 +__device__ __host__ inline void calc_vdw_geometric( + double ai_aii, double aj_aii, double ai_bii, double aj_bii, + double r6, double *V_a, double *V_b) +{ + *V_a = r6 * r6 * ai_aii * aj_aii; + *V_b = r6 * ai_bii * aj_bii; +} + +// Arithmetic rule: R*_ij = R*_i + R*_j, eps_ij = sqrt(eps_i * eps_j) +// Energy: V = sqrt(eps_ij) * R6^2 * r^-12 - 2*sqrt(eps_ij) * R6 * r^-6 +// where R6 = (R*_i + R*_j)^6 +// Parameters: For arithmetic rule, the same storage is reused: +// ai_aii, aj_aii store R*_i, R*_j (vdW radius) +// ai_bii, aj_bii store sqrt(eps_i), sqrt(eps_j) (after preprocessing) +// r6 is 1/r^6 +__device__ __host__ inline void calc_vdw_arithmetic( + double Rstar_i, double Rstar_j, double sqrt_eps_i, double sqrt_eps_j, + double r6, double *V_a, double *V_b) +{ + double Rstar_ij = Rstar_i + Rstar_j; // Arithmetic combination + double sqrt_eps_ij = sqrt_eps_i * sqrt_eps_j; // Geometric combination (already sqrt) + + // Compute R6 = (R*_ij)^6 + double R2 = Rstar_ij * Rstar_ij; + double R6 = R2 * R2 * R2; + + *V_a = sqrt_eps_ij * R6 * R6 * r6 * r6; // sqrt(eps_i * eps_j) * R^12 * r^-12 + *V_b = 2.0 * sqrt_eps_ij * R6 * r6; // 2 * sqrt(eps_i * eps_j) * R^6 * r^-6 +} + +#endif /* __VDW_RULES_H__ */