第2章 GPU计算基础与CUDA编程简介

好，我们直接进入正题。做分子对接加速，GPU 是你的火箭引擎，CUDA 就是你的发射台控制器。我见过不少做计算的同行，拿到 Vina 源码就想往 GPU 上怼，结果连最基础的线程模型都没搞明白，跑出来的结果还不如 CPU 快。今天咱们把底子打牢。

2.1 GPU 为什么能加速分子对接？

CPU 和 GPU 的设计哲学完全不同。CPU 是聪明能干的总经理，擅长处理复杂的逻辑；GPU 是成千上万的流水线工人，专干重复的体力活。

分子对接的本质是什么？说白了就是算——算配体在受体口袋里每个位置、每个扭转角的打分，然后挑出最低能量的那个。这中间有海量的

独立的

你想想看，Vina 打分时，每个配体构象的计算互相独立。这不就是典型的"数据并行"吗？把数万个构象同时丢给 GPU，一次就出结果——这种场景是我们 GPU 加速的红利区。

2.2 CUDA 编程模型：从 "Hello World" 到线程网格

CUDA 的编程模型，核心就两句话：主机（Host）控制流程，设备（Device）执行计算。

我个人习惯先画一张图来理解执行模型。你看下面这个流程，从调用 kernel 到线程实际执行，中间经过了网格、线程块、线程三级抽象。

这段代码是一个典型的 CUDA kernel 调用模板，你可以对照着上面的图看——这其实是 Vina 打分函数中常用的模式。

// 定义 kernel：每个线程计算一个配体构象的能量
__global__ void score_ligands_kernel(
    float* d_scores,       // 输出：每个构象的分数
    float* d_coords,       // 输入：所有构象的原子坐标
    float* d_receptor,     // 输入：受体原子坐标和参数
    int num_conformations  // 构象总数
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if (idx >= num_conformations) return;

    // 当前线程计算第 idx 个构象的分数
    float score = 0.0f;
    for (int i = 0; i < NUM_RECEPTOR_ATOMS; i++) {
        score += compute_interaction(d_coords + idx * NUM_LIGAND_ATOMS * 3, d_receptor + i);
    }
    d_scores[idx] = score;
}

// 主机端调用
int blocks = (N + 255) / 256; // 向上取整，每块 256 线程
score_ligands_kernel<<<blocks, 256>>>(d_scores, d_coords, d_receptor, N);

key 点：

• blockIdx.x + threadIdx.x 构成全局唯一线程 ID，这是寻址的基础
• 每个线程独立访问不同的数据区，不存在写冲突——这就是完美的数据并行
• block 大小一般为 128 或 256 个线程，这是经验值，Vina 场景用 256 比较稳妥

2.3 内存层次：让数据流动起来

GPU 的内存不是一块平的大仓库，它有严格的层级和速度差异。我刚开始搞 CUDA 时，在这上面栽过跟头——总想着把所有数据都塞到显存里，结果忽略了数据在 Host 和 Device 之间的搬运开销。

我个人在优化 Vina 的路子上，最常用的组合是：全局内存存大数组 + 共享内存做中间缓存 + 寄存器算局部变量。 你想想看，计算每个构象时，受体的坐标是不变的，如果每个线程每次都从全局内存读一遍受体坐标，那带宽就炸了。所以我通常把这些数据先搬到共享内存里，让一个 block 的所有线程共用。

2.4 开发环境配置：别卡在第一步

好，理论撸完了，咱们来点实战的。配置工具链这一步其实没啥难度，但容易踩错。我建议你从 CUDA 12.x + GCC 11.x 开始，兼容性最稳。

2.4.1 检查 CUDA 和 GPU 设备

# 检查 CUDA 版本
nvcc --version

# 查看 GPU 设备信息
nvidia-smi

# 显示 GPU 计算能力（Compute Capability）
nvidia-smi --query-gpu=compute_cap --format=csv

2.4.2 编译你的第一个 CUDA 程序

# 单文件编译
nvcc -o test test.cu -arch=sm_80 -O3

# CMake 集成（推荐，Vina 项目也用 CMake）
cmake_minimum_required(VERSION 3.20)
project(VinaGPU)
set(CMAKE_CUDA_STANDARD 17)
enable_language(CUDA)

find_package(CUDAToolkit REQUIRED)

add_executable(vina_gpu main.cu)
target_compile_options(vina_gpu PRIVATE -arch=sm_80 -O3 -lineinfo)

注意 -arch=sm_80 对应的是 Ampere 架构（RTX 30 系列/A100），如果你用 RTX 40 系列（Ada Lovelace）就换成 sm_89，用 RTX 20 系列（Turing）换 sm_75。这个参数编译出来是给哪个架构用的，跑在别的架构上虽然也能跑，但性能会打折扣。我习惯查表确认 GPU 的计算能力，再写对应的 arch flags。

2.5 实战：Vina 中 CUDA 代码的第一步走读

我们把刚才的理论拼起来，看看 Vina 实际的 CUDA 代码长什么样。下面这个代码片段来自 Vina-GPU 的 scoring 核心，不要觉得长——拆开看，每一行你都能对上号。

// 实际 Vina-GPU 源码风格（简化版）
#define BLOCK_SIZE 256

__global__ void vina_score_kernel(
    float* d_scores,
    const float* d_ligand_coords,     // [num_conformations][num_lig_atoms][3]
    const float* d_receptor_coords,   // [num_rec_atoms][3]
    const float* d_vdw_params,
    const float* d_hbond_params,
    int num_conformations,
    int num_lig_atoms,
    int num_rec_atoms
) {
    int conf_idx = blockIdx.x * BLOCK_SIZE + threadIdx.x;
    if (conf_idx >= num_conformations) return;

    // 当前线程处理第 conf_idx 个构象的总分
    float total_score = 0.0f;

    // 注册变量保存当前配体构象的原子坐标偏移
    int lig_offset = conf_idx * num_lig_atoms * 3;

    // 遍历受体原子，计算所有相互作用
    for (int i = threadIdx.x; i < num_rec_atoms; i += BLOCK_SIZE) {
        float rx = d_receptor_coords[i * 3 + 0];
        float ry = d_receptor_coords[i * 3 + 1];
        float rz = d_receptor_coords[i * 3 + 2];

        // 内联打分函数：范德华 + 氢键 + 静电
        for (int j = 0; j < num_lig_atoms; j++) {
            float lx = d_ligand_coords[lig_offset + j * 3 + 0];
            float ly = d_ligand_coords[lig_offset + j * 3 + 1];
            float lz = d_ligand_coords[lig_offset + j * 3 + 2];

            float dx = rx - lx;
            float dy = ry - ly;
            float dz = rz - lz;
            float dist2 = dx*dx + dy*dy + dz*dz;

            total_score += compute_vdw(dist2, d_vdw_params[j], d_vdw_params[i]);
            total_score += compute_hbond(dist2, d_hbond_params[j], d_hbond_params[i]);
        }
    }

    // 归约：将同一 block 的线程结果累加
    __shared__ float partial_sum[BLOCK_SIZE];
    partial_sum[threadIdx.x] = total_score;
    __syncthreads();

    for (int s = BLOCK_SIZE/2; s > 0; s >>= 1) {
        if (threadIdx.x < s) {
            partial_sum[threadIdx.x] += partial_sum[threadIdx.x + s];
        }
        __syncthreads();
    }

    if (threadIdx.x == 0) {
        d_scores[conf_idx] = partial_sum[0];
    }
}

这段代码有四个关键点值得你记住：

1. 线程-构象 1:1 映射：每个线程生命周期内只处理一个构象，这样不会出现同步撕扯
2. 共享内存做归约：受体原子的遍历循环里，每个线程都算自己的累加值，最后在 block 内做一次 tree-based 求和——这是经典并行归约模式
3. 循环分块：for (int i = threadIdx.x; i < num_rec_atoms; i += BLOCK_SIZE) 保证同一个 block 的线程对受体原子连续存取，合并访存，带宽利用率高
4. __syncthreads()：共享内存写完后必须同步，否则多线程会读到对方还没写进去的值，算出来的结果就是错的

2.6 性能分析的第一个抓手：占用率

代码写完了，跑起来后怎么知道它发挥了多少 GPU 的功力？我最先看的指标永远是 占用率（Occupancy）：

# 使用 nvprof 或 Nsight Compute 获取占用率
ncu --set full ./vina_gpu

# 简化版：查看 kernel 的占用率
ncu --metrics sm__occupancy ./vina_gpu

占用率 = (活跃的线程束数 / 理论最大线程束数) × 100%。如果低于 50%，你大概率有这些毛病：

• 每个线程用的寄存器太多，导致每个 SM 能容纳的线程束数减少（用 __launch_bounds__ 可以限制寄存器用量）
• 共享内存用太多，超过每个 block 的配额
• block 维度配比太小，SM 的线程束不够满

嗯，到这基本把 GPU 计算最核心的东西讲清楚了。你从理解 GPU 为什么快，到看懂 Vina 的 CUDA kernel，再到怎么配置环境、怎么初步分析性能——这条路就算是铺出来了。后面几章我们会深入每一个细节，但基础不牢地动山摇，这一章我讲得慢了点，希望你吃透了再往下走。