생성형 AI를 활용한 GPU 가속 기반의 분자 동역학, 화학 반응 시뮬레이션

0. 누구세요?

나는 화학공학과 학생이다. 그냥 아주 평범한.

얼마 전까지는 2학년이라고 할 수 있었는데, 이제 3학년 1학기를 눈앞에 둔. 그냥 그런 학생이다.

그리고 교내 공모전에서 입상에 실패한. 후술할 내용들은 교내 공모전에 출품했던 내용인데, 음. 떨어져서 슬퍼서 네 그냥 네 그래서 써봐요...

그냥누구한명이라도내가이런걸했다는걸알아줬으면했어.............

개발자들로 가득한 공간에 비전공자, 그것도 디지털이 아닌 실물을 다루는 사람이 글을 쓴다니. 뭔가 언어, 문화 그 아무것도 모른 채로 외국 어딘가에 갑자기 떨어진 기분이라 솔직히 좀 불안하긴하다...

뭔가 이상해도... 너그럽게 용서해주셨으면 하는 약간의 이기적인 바램을 가져봅니다...

1. 그래서 뭘 했나요?

내가 했던건 가장 원초적인 형태의 분자동역학(molecular dynamics) 시뮬레이션을 만드는것이었다.

화학공학이라는 학문은 물질세계를 구성하는, 어떠한 성질의 가장 작은 단위인 분자를 생산하는 과정을 탐구하는 학문이다.

이를 위해 학부 2학년때에는 주로 생산 대상에 대한 이해인 화학을 주로 배웠는데, 내가 2학년 전공을 학습하며 배운 기본적인 화학의 원리는 다음과 같았다.

1. 분자 또한 물리학적 법칙에 따라 동작한다.
2. 분자의 물리학적 거동은 셀 수 없이 많은 분자들의 상호작용으로 하나의 현상(열역학)을 이룬다.
3. 그 열역학적 현상으로 인해 반응이 발생한다.

그래서 나는, 실제 분자의 동작 원리로부터 생겨난 현상을 직접 눈으로 확인해보고 싶었다. 컴퓨터 시뮬레이션이라는 방법으로.

2. 문제의 정의

어떠한 계의 열역학적 특성과 그로인해 나타나는 결과들을 파악하기 위해선 많은 양의 분자들을 필요로 한다.

즉, 이 문제는 컴퓨터를 사용해 많은 양의 연산을 처리하여 해결할 수 있는 문제라고 판단했다.

그리고, 필요한 연산은 복잡한 처리가 아닌 단순한 사칙연산의 반복(수만개 입자에 대한 물리법칙 연산)이므로, 병렬 연산을 통해 처리 속도를 극대화 시킬 수 있을 것이라고 판단하였다.

즉, GPU를 활용한 분자동력학 시뮬레이션을 만들어

1. 내가 만든 시뮬레이션의 물리학적 정합성을 확인하고
2. 실제 화학 반응 법칙이 기본적 물리학적 법칙으로부터 시스템 수준에서 구현되는지 확인하는것

을 목표로 하였다.

3. 구현해야 할 기능들

위의 문제를 해결하기 위해선 다음과 같은 기술적 과제를 달성해야 했다.

1. 입자의 구현 : 위치, 속도, 질량, 입자의 종류 정의
2. 속도 부여 : 초기 온도에 따른 속도 부여 및 운동량 보정
3. 입자의 운동 : 매 시간(dt)마다 위치 갱신
4. 충돌 처리
5. 1. 입자-입자 간 탄성 충돌 및 운동량 보존
6. 2. 입자-벽 간 경계 조건 처리
7. 화학 반응의 구현 : 활성화 에너지 조건에 따른 입자의 Type 변환

4. 사양

본 시뮬레이션 구동에 사용산 컴퓨터 사양은 다음과 같다.

5. 프로젝트 결과

프로젝트 github 링크

앞서 제시한 두 가지 목표

1. 내가 만든 시뮬레이션의 물리학적 정합성을 확인하고
2. 화학 반응 법칙이 기본적 물리학적 법칙으로부터 시스템 수준에서 구현되는지 확인하는것

에 대한 결과이다.

5.1. 시뮬레이션의 물리적 정합성 확인

다음의 조건에서 입자 기반 시뮬레이션을 수행하였다.

• 입자 수: 50,653개 (각 축 37개 격자 배치)
• 공간 크기: (6.7 × 10⁻⁸ m)³
• 시간 간격: 5e-16s
• 총 시뮬레이션 시간: 3e-10s (0.3ns, total 60만 step)

아래의 그림과 같은 초기 무작위 운동량 분포의 상황에서 충분한 시간을 지난 시점에서 다음의 볼츠만 분포(이론값, 붉은 선)에 피팅되는 것을 확인할 수 있었다.

또한, 두 시점에서의 온도는 각각 297.31K와 297.51K로, 계의 전체 내부에너지 또한 오차범위 내에서 보존되었음을 확인할 수 있었다.

5.2. 화학 반응 법칙의 구현

A + A -> B + B 이차반응의 구현

하단 이미지 클릭시 시각화 유튜브 영상으로 이동

다음의 조건에서 입자 기반 시뮬레이션을 수행하였다.

• 입자 수: 50,653개 (각 축 37개 격자 배치)
• 공간 크기: (6.7 × 10⁻⁸ m)³
• 시간 간격: 5e-16s
• 총 시뮬레이션 시간: 3e-09s (3ns, total 600만 step)
• 데이터 저장 간격: 1000step

1000 step당 평균 연산 시간은 약 1.6초로 측정되었다.

A + A <-> B + B 평형반응의 구현

하단 이미지 클릭시 시각화 유튜브 영상으로 이동

시뮬레이션 초기값은 다음과 같다.

• 입자 수: 50,653개 (각 축 37개 격자 배치)
• 공간 크기: (6.7 × 10⁻⁸ m)³
• 시간 간격: 5e-16s
• 총 시뮬레이션 시간: 5e-09s (5ns, total 1000만 step)
• 데이터 저장 간격: 5000step

A + B <-> C -> P 사전 평형 반응의 구현

하단 이미지 클릭시 시각화 유튜브 영상으로 이동

시뮬레이션 초기값은 다음과 같다.

• 입자 수: 50,653개 (각 축 37개 격자 배치)
• 공간 크기: (6.7 × 10⁻⁸ m)³
• 시간 간격: 5.0 × 10⁻¹⁶ s
• 총 시뮬레이션 시간: 5e-09 s (5ns, total 1000만 step)
• 데이터 저장 간격: 5000step

반응 초기 고농도 구간에서는 충돌 빈도가 높아 1000 step당 약 1.8초의 연산 시간이 소요되었으나, 반응이 진행되어 농도가 감소함에 따라 약 1.4초로 연산 시간이 감소하였다.

또한 위의 시뮬레이션 결과들과 실제 교과서 내의 이론적 그림 및 수식과 일치하는것을 확인하였다. 아래의 그래프가 사전 평형 반응의 실제 이론적인 그래프.

즉, 기본적 물리법칙(분자의 단순 직선운동)으로부터 계의 열역학적 특성과, 그에 의한 결과인 반응법칙또한 구현할 수 있음을 확인하였다.

5.3 프로젝트 파일 구조

Boltzman/
├── config/ # 실험 조건별 설정 파일 (JSON)
├── script/ # 소스 코드 모듈
│ ├── core/ # 시뮬레이션 핵심 로직 (Engine)
│ │ ├── boltzman.py # 시뮬레이션 실행
│ │ └── particle_setup.py # 입자 설정 코드
│ ├── cuda_kernel/ # GPU 가속을 위한 CUDA 커널 (.cu)
│ │ ├── build_list.cu # linked list 생성
│ │ └── collision.cu # 입자 충돌 이벤트 구현 (충돌 및 반응 제어)
│ ├── analysis/ # 데이터 시각화 및 분석 도구
│ │ ├── check_reaction.py # type별 입자의 농도변화 시각화
│ │ ├── plot_momentum.py # 입자의 운동량 분포 시각화
│ │ ├── visualize.py # 입자의 위치 시각화
│ │ └── visualize_reaction.py # type별 입자의 위치 시각화
│ ├── utility/ # 초기화 및 보조 유틸리티
│ │ └── project_setup.py # config 읽기 및 시뮬레이션 보조기능 수행
│ └── main.py # 프로그램 실행 진입점 (Entry Point)
└── result/ # 실험 결과 자동 저장소
  └── run_n/ # 실행 회차별 아카이빙
    ├── config.json # 해당 회차의 설정 백업
    ├── simulation.log # 실행 로그
    ├── distribution/ # 운동량/반응 속도 그래프
    └── step/ # 입자 데이터 (npz)

5.4 시뮬레이션 순서도

전체 프로세스

시뮬레이션 부분의 세부 프로세스

5.5. 실제 개발 일정

1일차 : 충돌 알고리즘 설계 및 기본 시뮬레이션 코드 구현 2일차 : 디버깅(충돌 문제, 온도의 발산 문제) 3일차 : 디버깅 완료 후 결과 검증, 코드 구조 개선 및 문서화 4일차 : 활성화 에너지 기반 화학 반응 모델 구현

6. 프로젝트의 초기 구상

이걸 만들기 전에 내가 알고 있던 부분은 이정도.

1. Python 기본 문법 (반복문, 조건문, 함수 이런거 할줄 암)
2. C++ 기본 문법 (코딩 맨 처음 할 때 dev C++ 써서 해봤음. 포인터도 학교에서 배움)
3. 자료구조 및 알고리즘의 존재 (그렇다고 잘 알지는 못함)
4. numpy와 matplotlib 사용법 (예전에 전기장 모델링한다고도 써봤고 아무튼 종종 자주 씀)
5. GPU가 병렬 연산에 특화되어있다는 사실 (어떻게 하는지는 모름)

그리하여, 초기 구상은 이랬다.

1. 입자를 정의하는데 필요한 정보는 총 8개 (위치 3개, 속도 3개, 질량, 반지름)
2. 충돌은 벽-입자, 입자-입자 충돌이 있음. 벽-입자는 범위 넘어가면 보정, 입자-입자는 반지름 합보다 위치 가까우면 충돌로 생각
3. 입자-입자 충돌은 nC2번, 즉 N^2번 하면 되지 않을까....

이렇게 해서 CPU에 넣고 돌리니까 8000개 입자만으로도 1step에 1분 30초가 걸렸다. 그래서 코드 들고 AI에게 총총총...

7. 입자 생성 및 속도 부여

본 프로젝트의 목적이 많은 양의 분자에 무작위 속도를 부여하여 실제 열역학적 성질을 구성하는지, 반응속도론이 적용되는지 확인하기 위함이므로

1. 입자가 존재햐야 하며
2. 입자에 속도가 들어가 있어야 하고
3. 입자가 충돌해야 한다.
   

아래와 같이 입자의 생성 및 배치를 하는 함수를 작성할 수 있었다.

# particle_setup.py

def setParticle(sx, sy, sz, num, r, m): 
    N = num**3 

    # 입자의 배열 간격을 공간의 크기 / 개수로 설정
    dx = sx / num
    dy = sy / num
    dz = sz / num

    # 격자 인덱스 생성 (0, 1, 2 ... num-1)
    gx, gy, gz = np.meshgrid(np.arange(num), np.arange(num), np.arange(num), indexing='ij') 
    gx, gy, gz = gx.ravel(), gy.ravel(), gz.ravel() 

    # 위치 계산 : 격자 인덱스 * 간격 + 반칸 오프셋
    # dx/2.0을 더해주는 이유는 0번 입자가 벽에 딱 붙지 않고, 자기 칸의 중앙에 오게 하기 위함.

    # + 약간의 무작위성(Jitter) 추가. 기체 상태에서는 결정성이 없으므로
    rng = np.random.default_rng()
    jitter_scale = 0.2 

    jitter_x = (rng.random(N) - 0.5) * (dx * jitter_scale)
    jitter_y = (rng.random(N) - 0.5) * (dy * jitter_scale)
    jitter_z = (rng.random(N) - 0.5) * (dz * jitter_scale)

    pos_cpu = np.stack([
        gx * dx + dx/2.0 + jitter_x,
        gy * dy + dy/2.0 + jitter_y,
        gz * dz + dz/2.0 + jitter_z
    ], axis=1).astype(np.float64)

    # 속도와 물성 배열 생성(값은 0)
    vel_cpu = np.zeros((N, 3), dtype=np.float64)
    prop_cpu = np.zeros((N, 2), dtype=np.float64)

    # 입자의 종류 배열 생성
    type_cpu = np.zeros((N,), dtype=np.int32)

    rng = np.random.default_rng()
    rand_vals = rng.random(N)

    # 입자 종류 배정
    type_cpu[rand_vals < 0.5] = 1 # 입자 A
    type_cpu[rand_vals >= 0.5] = 2 # 입자 B

    # 물성치 r,m 넣음
    prop_cpu[:, 0] = r
    prop_cpu[:, 1] = m

    return pos_cpu, vel_cpu, prop_cpu, type_cpu

아무튼, 이렇게 되면 입자의 위치 데이터가 이쁘게 생긴다.

속도값은 모든 값이 0인 배열이 되고, 속도와 물성은 config의 r,m값으로 저장된다.

사실 지금 하는 시뮬레이션의 수준에서는 r,m을 따로 어딘가 배열에 저장할 필요가 없긴 하지만, 나중에 입자가 충돌했을 때 활성화 에너지 이상을 가지면 반응하는 식으로도 만들어 보고 싶어서 일단은 r,m도 따로 두었음.
사실 이거는 한번 수정을 거친 버전인데, 구버전의 코드는 다음과 같았다.

# particle_setup.py (구버전)

#시뮬레이션의 초기 조건 설정

#1. 전체 공간 정의 및 입자 생성.
#sx, sy, sz : 공간의 크기, num : 한쪽 방향으로 배치되는 입자의 개수, r : 입자의 반지름, m : 입자의 질량
def setParticle(sx, sy, sz, num, r, m):

    #전체 공간의 크기 : 가로 * 세로 * 높이
    size_space = sx*sy*sz
    size_num = num*num*num

    if (r+0.5)*2 > sx/num or (r+0.5)*2 > sy/num or (r+0.5)*2 > sz/num:

        print(f"오류 : 입자의 격자가 공간보다 큽니다.")
        return None, None, None

    #num^3개의 파티클 생성

    #1) 위치 정보 : x,y,z
    positions = np.zeros((size_num, 3), dtype = float)
    #2) 속도 정보 : vx,vy,vz
    velocities = np.zeros((size_num, 3), dtype = float)
    #3) 물질 정보 : r,m
    properties = np.zeros((size_num, 2), dtype = float)

    #격자의 각 위치에 각각의 파티클 배치

    #배치의 시작지점 결정. x,y,z가 가장 작은 모서리부터 x->y->z순으로 배열.
    x_start = (sx / 2) - ((num - 1) * ((r+0.5) * 2)) / 2
    y_start = (sy / 2) - ((num - 1) * ((r+0.5) * 2)) / 2
    z_start = (sz / 2) - ((num - 1) * ((r+0.5) * 2)) / 2

    count = 0

    #x>y>z 순으로 격자에 입자 배치
    for i in range(num):
        current_z = z_start + i * (r+0.5)*2

        for j in range(num):
            current_y = y_start + j * (r+0.5)*2

            for k in range(num):
                current_x = x_start + k * (r+0.5)*2

                positions[count] = [current_x, current_y, current_z]
                velocities[count] = [0.0, 0.0, 0.0]
                properties[count] = [r,m]


                count = count + 1

    return positions, velocities, properties

나는 원시적으로 반복문 3번 돌려서 배열하는걸 생각했었는데, AI가 이게 더 빠를거에염 하고 단순히 행렬계산 몇번으로 딱딱딱 해버리는 코드를 만들어서 줬다. 엄청 신기했음.

아무튼, 이렇게 하면 공간에 격자 형태로 입자가 배치되지만 아직 속도는 들어가있지 않다.

속도는 물리화학 시간에 했던 기체분자운동론 수식으로 계산.

무작위하게 뽑았을지라도 계의 운동량이 0이 되지 않을 수 있어서 보정하는 부분도 넣었다.

# particle_setup.py

# particle의 초기 속도 설정 함수
def setParticleSpeed(velocities, properties, T, m):

    k_B = 1.380649e-23

    v = math.sqrt((k_B * 3 * T) / m)

    # velocities 배열에 v를 매개로 한 무작위 값 넣음. velocities의 dtype을 가지며, velocities와 동일한 형태를 가지는 새 배열 생성.
    # 2를 곱하고 1을 빼는 것으로 값을 [0, 1]의 값을 [-1, 1]의 배열들로 변경. 이후 상수 v 곱해 스케일을 크게 만듬.
    # velocities[:] : 전체 배열을 의미. 새로 생성한 배열을 덮어쓰기함
    rng = cp.random.default_rng()
    velocities[:] = (rng.random(velocities.shape, dtype=velocities.dtype) * 2 - 1) * v # 속도 난수 형성. 

    # 계의 전체 운동량 계산. Total 운동량을 0으로 만들기 위한 보정작업. 이론상으로는 무작위이므로 전체 속도의 합이 0일 것이지만, 큰 수의 법칙이 적용되지 않는 영역에서는 보정이 필요.
    mass = properties[:, 1].reshape(-1, 1) # 입자의 모든 질량값을 가져옴. velocities와 계산하기 위해 [N*1] 벡터로 저장. 
    total_mom = (velocities * mass).sum(axis=0, keepdims=True) # 각 입자의 속도 벡터에 질량을 모두 곱해 운동량 계산. .sum(axis=0, keepdims=True)를 통해 열 방향(각 차원의 운동량 방향)의 값을 모두 더하고 [1*3] 행렬의 형태로 만듬
    total_mass = mass.sum()

    # 전체 운동량의 합계가 0보다 크다면 편차만큼을 전체 속도에서 빼줌
    if total_mass > 0:
        velocities -= total_mom / total_mass

    return velocities

AI가 짜주긴 했는데 이게 대체 뭔가 싶어서 한줄한줄 읽으며 주석을 주렁주렁 달아도보았다.

0,1의 값을 2 곱하고 1빼는걸로 -1에서 1 만드는것도 신기하고 무엇보다 솔직히 행렬연산이 익숙하지 않다. 선형대수학을 따로 배운것도 아니고, 저번학기에 공학수학 1에서 상미분방정식 푸려고 행렬 조금 배운것밖에 없어서...

왜 컴공에서 선형대수학을 배우는지 이거 하면서 뼈저리게 느꼈음. 그냥 코드 몇줄 딸깍하니까 계산이 다 되어있더라... 그냥 궁금해서라도 시간 내서 선형대수학 조금 공부해볼 생각.

위의 두 코드에서 setparticle은 np, setParticleSpeed는 cp로 되어있는데 속도 설정은 그래픽카드 사용해서 처리했다. 그래서 중간에 이런 과정을 거쳐야한다.

    # 2. 입자 생성 및 CPU 데이터를 GPU로 전송
    pos_cpu, vel_cpu, prop_cpu, type_cpu = setParticle(sx, sy, sz, num_side, r, m) # 입자의 Position 생성
    positions = cp.asarray(pos_cpu)
    velocities = cp.asarray(vel_cpu)
    properties = cp.asarray(prop_cpu)
    types = cp.asarray(type_cpu)
    N = positions.shape[0] # 생성된 입자의 전체 개수 정의
    velocities = setParticleSpeed(velocities, properties, T, m) #입자에 속도 부여

8. 입자의 충돌 구현

충돌을 하기 이전에 입자의 운동부터 생각해보자.

사실 운동은 매우 간단하다. 맨 처음 들었던 생각은 그냥 dt 간격으로 속도에 곱해 더해주면 되는것. 지금의 가정에서는 입자가 전혀 가속을 받고 있지 않으므로(떠돌아다닐 때 등속운동) 크게 부담없이 적용할 수 있다고 생각했다.

그러면 다음으로 간단한거는 벽과의 충돌인데, 맨 처음 짰던 코드는 이거였다.

# boltzman.py (구버전)

#3. 시간 간격에 따른 입자의 운동 시뮬레이션
def moveOnParticle(positions, velocities, properties, dt, sx, sy, sz, grid_params):

    num = positions.shape[0]

    #입자의 이동 연산
    positions = positions + velocities * dt


    #벽과의 충돌 연산
    collision_x = (positions[:, 0] < 0) | (positions[:, 0] > sx)
    collision_y = (positions[:, 1] < 0) | (positions[:, 1] > sy)
    collision_z = (positions[:, 2] < 0) | (positions[:, 2] > sz)

    #충돌횟수 count
    wallCollision = collision_x | collision_y | collision_z
    wallCollision_count = wallCollision.sum()

    #운동량변화 (벽에 대한 충격량)
    Pressure_x = 2 * velocities[collision_x, 0] * properties[collision_x, 1]
    Pressure_y = 2 * velocities[collision_y, 1] * properties[collision_y, 1]
    Pressure_z = 2 * velocities[collision_z, 2] * properties[collision_z, 1]

    Sum_Pressure_x = Pressure_x.sum()
    Sum_Pressure_y = Pressure_y.sum()
    Sum_Pressure_z = Pressure_z.sum()

    Pressure = Sum_Pressure_x + Sum_Pressure_y + Sum_Pressure_z

    #속도 vector 방향 바꾸기
    velocities[collision_x, 0] *= -1.0
    velocities[collision_y, 1] *= -1.0
    velocities[collision_z, 2] *= -1.0

    #cliping
    positions[:, 0] = cp.clip(positions[:, 0], 0, sx)
    positions[:, 1] = cp.clip(positions[:, 1], 0, sy)
    positions[:, 2] = cp.clip(positions[:, 2], 0, sz)


    #입자간의 충돌 연산

    Total_Impulse_PP = cellCollision(positions, velocities, properties, grid_params, sx, sy, sz)


    return positions, velocities, wallCollision_count, Pressure, Total_Impulse_PP

진짜 아주아주 직관적이고 간단하죠...? 진짜 그냥 움직이고, 경계 넘으면 운동량 바꿔서 반대로 하고, 위치도 바꿔주고 입자간의 충돌연산 시켜주고.

근데 AI한테 맡기니까 뭔가 많이 달라졌다.

# boltzman.py

    # 시뮬레이션은 step 단위로 진행. 하나의 step 루프는 CPU에서 진행되므로, 이전 단계가 완료되어야 다음 단계로 진행됨.
    for step in range(steps): # step의 개수만큼 반복
        positions += velocities * dt # Particle를 움직이는 연산

        # 벽과의 충돌
        for axis, bound in enumerate([sx, sy, sz]): # 총 3번의 루프로 x축, y축, z축을 검사.
            mask = (positions[:, axis] < 0) | (positions[:, axis] > bound) # 입자의 위치가 경계 밖에 있는지를 판단
            if mask.any(): 
                velocities[mask, axis] *= -1.0 # 해당 방향의 속도를 반대로 바꾸어 튕겨내는 연산 수행
            positions[:, axis] = cp.clip(positions[:, axis], 0.0, bound) # position을 조정(경계 밖에서 안으로 데려오기)

움직이는것까지야 워낙 간단해서 똑같은데, 벽과의 충돌 부분이 잉? 싶을 정도로 뭔가 천지개벽했다...

이게 벽과의 충돌 부분에서 [sx, sy, sz] 이걸 돌리는거니까... 루프가 이렇게 돌아간다.

1번째 : axis = 0(sx), bound = sx(x축 경곗값) 2번째 : axis = 1(sy), bound = sy(y축 경곗값) 3번째 : axis = 2(sz), bound = sz(z축 경곗값)
그러면 이 값을 가지고 비교 연산을 수행한다.

mask = (positions[:, axis] < 0) | (posions[:, axis] > bound)
# positions[:, axis] : 모든 입자의 현재 축(axis = 0 : x축, 1 : y축, 2 : z축)의 위치를 가져온다.
# 0, bound랑 비교하며 결과는 [True, False, False, ...] 이런 식으로 반환되어 mask에 저장된다. (조건 둘중하나 맞음 True)
# ! 0보다 작거나 bound보다 크면 당연히 해당 입자는 계 밖에 있는 것이다. (존재할 수 없는 값)

# 결론 : 벽 밖에 있는 모든 입자에 True가 표시되어 있다.

하나라도 True가 있는지 확인 (if문)하고, True가 표시된 입자의 해당 방향 운동량을 바꾼다. 그리고 경계 넘어간 입자의 위치를 보정해서 안으로 끌어온다

사실 이게 입자가 반지름이 가정되어있는지라 r만큼 더 땡겨서 가져오긴 해야하는데... 바꾸어 생각하면 계산 대상이 입자 중심이니까, 중심 밖으로 나가있는거니까 반지름만큼 더 땡겨오면 오차가 더 커지는거다. 아무튼 이게 맞는것.

진짜 큰 문제는 입자-입자 충돌이었다.

진짜진짜진짜 큰 문제는 입자간 충돌이었다.

맨 처음 생각했던 알고리즘은 100개의 입자가 있으면 1번째 입자를 나머지 99개의 입자와, 2번쨰 입자를 나머지 98개의 입자와 충돌할 수 있다고 생각하고 계산하는거였다.

근데 이 경우는 연산량이 말도안되게 많아서... 입자 8000개로 1step 돌리는데 1분 30초가 걸렸다. 왜이렇게 오래 걸리냐고 AI한테 물어보니까 돌아온 답변은
너가 만든 코드는 너무 구려요.

그래서 linked list를 구현해서 격자 알고리즘을 짜고, 지금 사용하는 cupy가 아닌 직접 CUDA 커널을 짜서 연산을 돌리면 더 빨라질거라는 AI의 답변을 받고

아니 나 잘 모르겠으니까 그냥 해줘. 라고 했더니 애가 뭔가 결과물을 주긴 했는데 잘 돌아가지는 않았다. 음. 완전 안돌아갔다.

그래서 이것저것 찾아보고 고생을 조금 하고 나서 나온 결과물은 이거였다. 사실 초기 버전은 이게 아닌데, 나중에 반응을 추가한 최종 결과물이다.

//collision.cu 일부 발췌

extern "C" __global__
void collide_particles(...)
{
    // 1. 연산의 시행 여부 결정

    int i = blockDim.x * blockIdx.x + threadIdx.x;
    if (i >= N) return;

    // ... (입자 정보 수집 및 격자 계산 생략)

    // 3. 입자가 존재하는 격자와 충돌 연산을 수행할 격자의 판별 (x,y,z로 -1, +1에 해당하는 입자와만 연산. 총 27개)

    for (int dz = -1; dz <= 1; ++dz) {
        int nz_c = iz + dz;
        if (nz_c < 0 || nz_c >= nz) continue; 
        for (int dy = -1; dy <= 1; ++dy) {
            int ny_c = iy + dy;
            if (ny_c < 0 || ny_c >= ny) continue;
            for (int dx = -1; dx <= 1; ++dx) {
                int nx_c = ix + dx;
                if (nx_c < 0 || nx_c >= nx) continue;

                // 4. 판별한 격자와의 충돌 연산을 시행할 격자 내 입자 가져오기

                // 양 옆 격자에 대한 판별 종료 및 격자의 index 구하기
                int cell = nx_c + ny_c * nx + nz_c * (nx * ny);

                int j = head[cell];

                // 해당 격자 내 모든 입자에 대한 연산을 수행
                while (j != -1) {

                    // 5. 입자와의 충돌 여부 판단
                    if (j > i) {
                        double xj = pos[3*j+0], yj = pos[3*j+1], zj = pos[3*j+2]; //입자 j의 위치 데이터 가져옴

                        // 입자 i와 j 사이의 입체적 거리(정확히는 거리의 제곱) 계산
                        double dxij = xi - xj;
                        double dyij = yi - yj;
                        double dzij = zi - zj;
                        double dist2 = dxij*dxij + dyij*dyij + dzij*dzij;

                        // 충돌의 판별. 두 입자의 반지름의 합보다 입자 사이의 거리가 작으면 충돌으로 판단
                        double r_j = props[2*j + 0];
                        double rsum = r_i + r_j;

                        if (dist2 > 1e-20 && dist2 < rsum*rsum) {

                            // ... (운동량 계산 연산 수행) ...

                            // 입자가 서로 가까워지고 있을 때. 즉, 상대속도가 음수일 때만 연산을 시행
                            if (v_rel_n < 0.0) {

                                // 화학 반응의 여부 판별
                                int type_j = types[j]; // 상대 입자의 종류 가져옴

                                if ((type_i == 1 && type_j == 2) || (type_i == 2 && type_j == 1)) {
                                    if (fabs(v_rel_n) > 500.0) {

                                        atomicCAS(&types[i], type_i, 3); //중간체로 변환
                                        atomicCAS(&types[j], type_j, 3);

                                    }
                                }else if (type_i == 3 && type_j == 3) {

                                    if (fabs(v_rel_n) > 1500.0) {
                                        atomicCAS(&types[i], 3, 4); //생성물로 변환
                                        atomicCAS(&types[j], 3, 4);
                                    }
                                }

                                // ...(충격량 연산 수행)...

                                // 입자의 겹침에 대한 보정 수행
                                double overlap = rsum - dist;

                                // 작은 겹침(부동소숫점 오차)은 무시
                                if (overlap > 1e-15){
                                    double correction = overlap * 0.5; // 각 입자를 절반씩 나누어서 밀어냄

                                    // 밀어낼 길이 계산 (x,y,z방향)
                                    double cx = nxu * correction;
                                    double cy = nyu * correction;
                                    double cz = nzu * correction;

                                    // 보정 수행
                                    atomicAdd(&pos[3*i+0], cx);
                                    atomicAdd(&pos[3*i+1], cy);
                                    atomicAdd(&pos[3*i+2], cz);

                                    atomicAdd(&pos[3*j+0], -cx);
                                    atomicAdd(&pos[3*j+1], -cy);
                                    atomicAdd(&pos[3*j+2], -cz);
                                }

                            }
                        }
                    }

                    // linked list에서 다음 입자 불러옴
                    j = nxt[j];
                }
            }
        }
    }
}

내용이 길긴 한데, 생각보다 간단하다. 전체적인 흐름을 정리하면

1. GPU에 연산(입자)을 할당한다.
2. 입자의 위치와 물성 정보를 가져온다. (충돌 여부 판단에 필요한 정보)
3. 입자가 어떤 격자에 속해있는지 계산한다. (자기 자신을 포함한 27개의 격자 내 입자만 계산할것이므로)
4. 계산을 수행할 격자를 가져오고, 격자 내 입자를 가져온다. (잠재적 충돌 대상)
5. 입자와의 충돌 여부를 판단한다.
6. 충돌했다면 입자의 속도를 가져오고 충돌 연산을 수행한다.
7. 만약 상대속도가 활성화에너지를 넘으면 입자의 Type을 변경하는 반응 연산을 수행한다.
8. 입자 겹침을 보정한다.
9. 다음 잠재적 충돌 대상 입자를 가져온다.
   

와 이걸 처음 생각한 사람은 대체 누구일까...?

아니 물론, 하나하나 생각해보면 아주아주 타당한 접근이다. 굳이 모든 입자에 대해 충돌 여부를 판단할 필요가 없고, 가장 충돌할 가능성이 높은 대상들을 뽑아 판별하면 되는건데...

이걸 구현한 방식이 신기했다. 아무튼, 격자에 입자를 배치하는 연산은 이렇게 구성되어 있다.

build_list.cu

// 입자에 격자를 배정하는 커널
extern "C" __global__
void build_linked_list(...)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x; // 코어에 할당된 스레드의 ID. 1개의 thread당 1개의 입자 연산을 수행
    if (i >= N) return; // 입자의 개수를 초과하면 연산을 할 필요가 없음 (오류만 발생)

    // 할당된 입자의 위치 정보 수집
    // pos 배열은 pos[0]-[1]-[2]... 이며 각각의 [0]-[1]-[2]는 [x,y,z]의 쌍으로 구성되어 있음.
    // 따라서, 3*i로 입자의 ID를 얻고, 0을 더하면 x, 1을 더하면 y, 2를 더하면 z 좌표를 얻을 수 있음. 
    double x = pos[3*i + 0];
    double y = pos[3*i + 1];
    double z = pos[3*i + 2];

    // 현재 입자가 몇번째 격자에 속해있는지 계산
    int ix = floor(x / cell_size);
    int iy = floor(y / cell_size);
    int iz = floor(z / cell_size);

    // 입자가 불가능한 위치(0번 격자와 마지막 격자의 밖)에 있다면 보정
    if (ix < 0) ix = 0; if (iy < 0) iy = 0; if (iz < 0) iz = 0;
    if (ix >= nx) ix = nx - 1; if (iy >= ny) iy = ny - 1; if (iz >= nz) iz = nz - 1;

    // linked list의 형성.

    // cell의 3차원 좌표를 1차원 ID로 변환. 변환된 모든 1차원 index는 고유함
    int cell = ix + iy * nx + iz * (nx * ny);

    int old = atomicExch(&head[cell], i);
    nxt[i] = old; 
}

9. AI 시대에 내가 뭘 할 수 있을까?

사실 이거 처음 할때는 그냥 이전에 하던대로 슥슥 하면 될줄 알았는데 생각보다 깊다고 생각하는 부분들이 많이 튀어나와서 놀랐다.

이게 생성형 AI의 장점이 아닐까...?

어떤 내용을 학습하는 것은 대강 몇가지 정도로 나뉠 수 있다고 생각하는데,

1. 지금 가지고 있는 배경지식으로도 충분히 쉽게 알 수 있는 내용
2. 지금 가지고 있는 배경지식 외에도 추가적인 지식이 필요한 내용
3. 완전히 배경지식조차 존재하지 않는 처음 보는 내용
4. 내 뇌의 성능을 완전히 벗어난 내용

1번은 쉽게 공부하고, 2번과 3번의 경우 대체 내가 뭘 더 알아야 하는지 모를 것 같은 때가 많은 것 같다. 이럴 때 AI가 참 너무나도 유용하다.

인터넷에서 검색을 할 때에는 검색어 키워드를 잘 생각해가면서, 그 아무짝에도 쓸모없는, 상관없는 내용들을 걸러가면서 하나하나 수고스럽게 찾아야 하는데 AI와 대화하다 보면 필요한 정보를 쉽게 얻을 수 있을 뿐만 아니라 내가 모르던 영역까지도 애가 알아서 가져와서 '이런것도 있어요' 하고 보여주니 새로운 영역으로의 확장이 조금은 쉬워지는 기분이다.

그래서인지 맨 처음에 쉽다고 생각하면서 접근했던 부분들이 뭔가 새로운 알고리즘이나 자료구조를 알게되고, 그걸 도입했을 때(물론 내가 한건 아니고 AI가 추천해주고, AI가 코드를 짜줬지만) 성능이 향상되는것을 눈으로 보니까 AI의 파괴적인 힘을 느낄 수 있었다.

그럼에도, 정말 그럼에도 아직까지 AI가 약간 삐걱거리는 케이스도 확실히 많이 있었다.

초기 버전에서 충돌처리가 되지 않는 문제가 있어서 수정하기도 했고,

시뮬레이션이 말도안되는 온도값으로 발산해버리거나 하는 부분은 AI가 버그 수정하는데 심각한 어려움을 겪길래 직접 케이스별로 하나하나 시도해보며 수정했으며,

무엇보다 일단 이걸 해보고싶다는 생각이 있어야 AI가 해준다는 것 정도...

그럼에도...

진심으로... AI시대에 내가 뭘 할 수 있을지... 살짝 고민하게되는 프로젝트였던 것 같다...

fin.