우산 표본추출

1. 개요[편집]

우산 표본추출 Umbrella Sampling
제안	1977년, G. M. Torrie & J. P. Valleau
편향 퍼텐셜	보통 조화형 w(ξ) = ½k(ξ − ξ₀)²
산출물	평균력 퍼텐셜(PMF) A(ξ)
복원 알고리즘	WHAM · MBAR · 우산 적분
필수 진단	이웃 창 히스토그램 겹침
대표 구현	PLUMED, Colvars, GROMACS, LAMMPS

우산 표본추출(umbrella sampling)은 자유에너지 장벽 때문에 평범한 시뮬레이션으로는 평생 가 보지 못할 영역을, 미리 고른 반응좌표 $\xi$ 에 인위적인 편향 퍼텐셜을 걸어 강제로 방문시킨 뒤, 나중에 그 편향의 효과를 통계적으로 정확히 벗겨 내 원래 앙상블의 자유에너지 프로파일을 복원하는 향상된 표본추출 기법이다. 1977년 토리와 발로가 제안했다.¹

동기는 단순하고 잔인하다. 정준 앙상블에서 상태를 방문할 확률은 $e^{-\beta A}$ 에 비례하므로, 높이 $20\,k_BT$ 짜리 장벽 너머는 사실상 확률 $10^{-9}$ 의 세계다. 분자동역학을 마이크로초 돌려도 전이가 한 번도 안 일어나고, 그러면 우리가 얻는 것은 “전이가 없다”가 아니라 **“아무것도 모른다”**이다. 표본이 0개인 구간의 히스토그램은 자유에너지가 무한대라고 말하지만 그건 통계가 아니라 침묵이다.

처방은 이거다. 장벽 꼭대기가 안 나오면 거기 있으라고 스프링으로 묶어 놓고 재라. 물론 그렇게 재면 값이 틀리는데, 얼마나 틀리는지는 정확히 안다(스프링을 우리가 직접 걸었으므로). 이 “일부러 틀리게 재고 정확히 되돌린다”는 발상은 중요도 표본추출 그 자체이며, 우산 표본추출은 그것을 반응좌표 공간에서 수행하는 것이다.

2. 편향을 걸고 벗기기[편집]

원래 퍼텐셜 $U(\mathbf{x})$ 에 반응좌표만의 함수인 편향 $w(\xi(\mathbf{x}))$ 를 더한다.

$$U^b(\mathbf{x}) = U(\mathbf{x}) + w(\xi(\mathbf{x})), \qquad w_i(\xi) = \tfrac{1}{2}k(\xi - \xi_i^0)^2$$

편향된 앙상블에서 잰 $\xi$ 의 분포 $P_i^b(\xi)$ 와 원래 분포 $P^0(\xi)$ 의 관계는 볼츠만 인자 한 줄에서 바로 나온다.

$$P_i^b(\xi) = \frac{P^0(\xi)\, e^{-\beta w_i(\xi)}}{\left\langle e^{-\beta w_i} \right\rangle_0} \quad\Longrightarrow\quad P^0(\xi) \;\propto\; P_i^b(\xi)\, e^{+\beta w_i(\xi)}$$

즉 히스토그램에 $e^{+\beta w}$ 를 곱해 주면 끝이다. 편향이 좌표 $\xi$ 에만 의존한다는 것이 핵심 — 그래야 직교 자유도에 대한 적분에서 편향 인자가 밖으로 빠져나온다. 그리고 자유에너지 프로파일(평균력 퍼텐셜, PMF)은

$$A(\xi) = -k_B T \ln P^0(\xi) + \text{const}$$

문제는 저 상수항이다. 창 $i$ 하나만 보면 규격화 상수 $\langle e^{-\beta w_i}\rangle_0$ 를 모르므로 PMF 조각의 세로 위치를 알 수 없다. 창마다 상수 $F_i = -k_BT \ln \langle e^{-\beta w_i}\rangle_0$ 가 미지수로 하나씩 붙고, 조각들을 하나의 곡선으로 이어 붙이는 것이 이 기법의 진짜 계산 문제다. 창이 30개면 미지수도 30개.

여기서 조용히 사람을 잡는 함정 하나. $\xi$ 가 데카르트 좌표가 아니면 야코비안 보정이 필요하다. 예컨대 3차원에서 두 입자 간 거리 $r$ 을 좌표로 쓰면 편향 없는 계에서도 $P(r) \propto 4\pi r^2$ 라 히스토그램이 저절로 오른쪽으로 쏠린다. 이걸 그대로 로그 취하면 존재하지도 않는 $-2k_BT\ln r$ 짜리 인력이 PMF에 얹힌다.²

3. 창 배치와 힘상수[편집]

실무의 절반은 창(window)을 어디에 몇 개 놓고 스프링을 얼마나 세게 걸지 정하는 일이다.

• 겹침이 전부다. 이웃 창의 히스토그램이 겹치지 않으면 두 조각의 상대 높이 $F_i - F_{i+1}$ 를 결정할 정보가 아예 없다. 곡선은 그 지점에서 임의의 값만큼 위아래로 끊긴다. 그래서 계산 후 제일 먼저 그려 봐야 할 그림은 PMF가 아니라 창별 히스토그램을 한 판에 겹쳐 그린 것이다. 낙타 등처럼 봉우리가 줄줄이 이어지되 서로 발치가 겹쳐야 정상.
• 힘상수의 절충. 조화 편향에서 표본 폭은 대략 $\sigma \approx 1/\sqrt{\beta k}$ 다. $k$ 가 크면 창이 좁아져 겹침을 확보하려 창을 많이 박아야 하고, $k$ 가 작으면 계가 장벽에서 미끄러져 골짜기로 도망가 정작 꼭대기를 안 잰다. 하한의 감은 편향이 PMF의 곡률을 이겨야 한다는 것 — 장벽 마루에서 $A'' < 0$ 이므로 $k > \lvert A''\rvert_{\max}$ 여야 유효 퍼텐셜 $A + w$ 가 그 창 안에서 볼록해진다. 물론 $A''$ 를 알면 이미 답을 아는 셈이라, 현실에서는 짧은 시험 계산으로 눈대중한다. 창 간격도 등간격일 이유가 없어서 장벽 근처는 촘촘하게, 평평한 골짜기는 성기게 까는 것이 정석이다.
• 초기 구조를 어떻게 만드느냐가 조용한 사고의 원천이다. 흔한 방식은 조종 MD(steered MD)로 한 번 쭉 끌어당기며 각 창의 시작 구조를 뽑는 것인데, 이러면 모든 창이 같은 궤적의 후손이라 직교 자유도가 전부 같은 값에 갇힌 채 “수렴”한다. 반대 방향으로도 끌어 두 세트를 만들어 히스테리시스를 보는 것이 최소한의 방어다.

4. WHAM으로 PMF 복원하기[편집]

창별 조각과 미지 상수 $\{F_i\}$ 를 동시에 푸는 표준 도구가 WHAM(Weighted Histogram Analysis Method)이다. 페렌버그-스벤센의 다중 히스토그램 재가중을 우산 표본추출에 맞춰 정리한 것으로,³ 각 구간의 추정치를 분산이 최소가 되도록 가중 평균한다는 원리에서 두 개의 방정식이 나온다.

$$P^0(\xi) = \frac{\displaystyle\sum_{i=1}^{S} n_i(\xi)}{\displaystyle\sum_{i=1}^{S} N_i\, e^{-\beta\left(w_i(\xi) - F_i\right)}}, \qquad e^{-\beta F_i} = \int d\xi\; P^0(\xi)\, e^{-\beta w_i(\xi)}$$

여기서 $n_i(\xi)$ 는 창 $i$ 가 그 구간에 넣은 표본 수, $N_i$ 는 창 $i$ 의 총 표본 수다. 구조를 읽으면 이렇다 — $\xi$ 근처의 표본을 낸 모든 창의 기여를 분자에서 합치되, 분모에서 각 창이 그 지점을 얼마나 편향해서 봤는지로 나눠 준다. 표본이 많고 편향이 유리한 창일수록 그 구간의 추정에 크게 기여한다.

두 식은 $P^0$ 과 $\{F_i\}$ 가 서로를 참조하는 자기무모순 방정식이라 닫힌 해가 없고, $F_i = 0$ 에서 출발해 번갈아 대입하며 수렴할 때까지 돌린다. 창이 수십 개면 보통 수백~수천 회면 붙지만, 겹침이 부실하면 반복이 끝없이 기어간다 — “WHAM이 수렴 안 한다”는 증상은 대개 알고리즘 문제가 아니라 데이터가 정보를 안 갖고 있다는 신호다.

5. MBAR — 히스토그램을 버리다[편집]

WHAM은 히스토그램을 만드는 순간 구간 폭이라는 임의의 손잡이를 도입한다. 좁게 자르면 잡음이, 넓게 자르면 편향이 생긴다. 그 극한(구간 폭 → 0)을 정면으로 취한 것이 MBAR(Multistate Bennett Acceptance Ratio)다. 표본 하나하나를 그대로 쓰며, 결합 자기무모순 방정식은

$$\hat{F}_i = -k_B T \ln \sum_{j=1}^{S} \sum_{n=1}^{N_j} \frac{e^{-\beta u_i(\mathbf{x}_{jn})}} {\displaystyle\sum_{k=1}^{S} N_k\, e^{\beta \left(\hat{F}_k - u_k(\mathbf{x}_{jn})\right)}}$$

이 되고, 여기서 $u_i$ 는 편향까지 포함한 상태 $i$ 의 퍼텐셜 에너지, $\mathbf{x}_{jn}$ 은 창 $j$ 가 낸 $n$ 번째 구조다. 창이 두 개일 때 정확히 BAR로 환원되며, 주어진 표본에서 점근적으로 분산이 최소인 추정량이라는 것이 증명되어 있다. 덤으로 표본별 가중치가 직접 나오므로 PMF뿐 아니라 임의의 관측량의 재가중 평균과 통계 오차를 같은 틀에서 뽑을 수 있다. 오늘날 새로 짜는 워크플로에서 WHAM 대신 MBAR을 쓰는 것이 사실상 국룰이고, WHAM은 “구간을 쓰는 MBAR의 근사”로 이해하면 정확하다.

단, 오차 막대를 곧이곧대로 믿으면 안 된다. 두 방법 모두 표본이 독립이라고 가정하는데, MD 궤적의 연속 프레임은 지독하게 상관되어 있다. 자기상관 시간으로 궤적을 솎아 내거나 유효표본크기로 $N_i$ 를 대체하지 않으면, 오차 막대만 예쁘게 작아진다.

6. 힘으로 우회하기[편집]

히스토그램을 아예 안 쓰는 길도 있다. PMF의 기울기(평균력)를 창마다 재고 적분하는 것이다.

• 우산 적분(umbrella integration). 각 창의 편향된 $\xi$ 분포가 대략 정규분포라고 보면, 그 창에서의 평균력이 표본 평균 $\langle\xi\rangle_i$ 와 분산 $\sigma_i^2$ 만으로 떨어진다.

$$\frac{dA}{d\xi} = \frac{1}{\beta}\frac{\xi - \langle \xi\rangle_i}{\sigma_i^2} - k\left(\xi - \xi_i^0\right)$$

창별 기울기를 표본 밀도로 가중 평균한 뒤 수치적분하면 PMF가 나온다. 창 간 상수 $F_i$ 를 푸는 반복 자체가 사라지므로 겹침 요구가 조금 느슨하고 잡음에도 강한 편인데, 대신 정규분포 가정이 깨지는 곳(장벽 마루의 비대칭, 이중 우물이 한 창에 들어온 경우)에서는 조용히 틀린다.

• 제약 MD / 블루문 앙상블. 스프링 대신 $\xi$ 를 아예 특정 값에 고정하고, 그 제약을 유지하는 데 든 라그랑주 승수(제약력)의 평균을 재면 그게 평균력이다. 다만 좌표를 고정하면 운동량 공간의 측도가 왜곡되므로 질량 계량 인자에 대한 보정(픽스먼 항)을 넣어야 하고, 빼먹으면 곡선이 계통적으로 기운다. 구속 동역학의 SHAKE/RATTLE 기계를 그대로 재활용할 수 있다는 것이 구현상의 장점.
• 적응적 바이어싱 힘(ABF)은 계산 도중 추정한 평균력을 실시간으로 상쇄해 $\xi$ 방향을 평평하게 만든다. 창을 나누지 않는다는 점에서 사상이 다르다.

7. 반응좌표라는 진짜 문제[편집]

우산 표본추출의 모든 알고리즘적 정교함은 $\xi$ 가 옳다는 전제 위에 서 있다. 이 전제가 깨지면 위의 어떤 추정량도 구해 주지 못한다.

$\xi$ 를 잘못 고르면 무슨 일이 벌어지는가. 계는 우리가 지정한 좌표를 따라 얌전히 이동하지만, 진짜 병목인 직교 자유도(단백질 곁사슬의 회전이성질체, 물 분자의 재배열, 고리 뒤집힘)는 각 창 안에서 여전히 갇혀 있다. 각 창은 자기 나름대로 잘 수렴한 것처럼 보이고, 히스토그램도 예쁘게 겹치고, WHAM도 잘 붙고, PMF도 매끄럽게 나온다. 그런데 값이 틀렸다.⁴ 이게 위험한 이유는 수렴 진단이 전부 통과한다는 점이다.

증상을 잡는 정직한 방법은 정해져 있다. 정방향·역방향으로 끌어 만든 두 창 세트가 다른 PMF를 주면(히스테리시스) 직교 자유도가 안 풀린 것이고, 각 창의 $\langle\xi\rangle$ 나 직교 관측량이 계속 표류하면 아직 평형이 아니다. 표준 처방은 이웃 창끼리 좌표를 주기적으로 교환하는 복제교환 우산 표본추출(REUS, 해밀토니안 복제교환의 일종)이다. 직교 자유도가 창을 오가며 섞이므로 갇힘이 크게 완화된다.

좌표를 2차원으로 늘리는 것도 방법이지만 창 수가 곱으로 늘어나고(30 × 30 = 900개) 3차원부터는 사실상 불가능하다. 그래서 여기서부터는 좌표를 늘리는 대신 더 나은 좌표 하나를 찾는 쪽이 이긴다. 커밋터 확률이나 확산 지도, 최근에는 학습된 집단변수가 그 자리를 노린다.

8. 이웃 기법과의 위치[편집]

기법	편향의 성격	강점	약점
우산 표본추출	고정된 조화 편향, 창 분할	창별로 완전 병렬, 진단이 명확	창·힘상수 설계 노동
메타다이내믹스	방문한 곳에 가우시안을 쌓는 적응 편향	미리 설계할 게 적음	수렴 판정이 미묘, 잘 채워야 평평
자유에너지 섭동·TI	물리 좌표가 아닌 $\lambda$ 를 따라 이동	연금술적 변화에 최적	물리적 경로의 PMF는 못 줌
복제교환(온도)	편향 없음, 온도를 올려 장벽을 낮춤	반응좌표를 안 골라도 됨	계가 커지면 복제 수가 폭발

정리하면 우산 표본추출과 메타다이내믹스는 물리적 좌표를 따라가는 자유에너지 곡면을 다루고, 자유에너지 섭동 계열은 비물리적 연금술 경로를 다룬다. 목적이 다르므로 경쟁이 아니고, 실제 논문에서는 리간드 결합 PMF를 우산 표본추출로 뽑고 그 위에 표준상태 보정을 얹은 뒤 FEP 결과와 교차 검증하는 식으로 함께 쓰인다. 온도 기반 복제교환은 반응좌표를 안 고른다는 점에서 매력적이지만, 필요한 복제 수가 자유도 수의 제곱근에 비례해 늘어나는 탓에 용매까지 넣은 큰 계에서는 값이 안 나온다. 그리고 이 모든 것의 천장은 결국 힘장이다 — 표본추출이 완벽해도 퍼텐셜이 틀렸으면 정확히 틀린 PMF로 수렴한다. 그건 검증 및 확인의 영역이다.

9. 관련 문서[편집]

• 분자동역학 · 랑주뱅 동역학 · 정준 앙상블
• 자유에너지 섭동 · 헬름홀츠 자유에너지
• 중요도 표본추출 · 마르코프 연쇄 몬테카를로 · 해밀토니안 몬테카를로
• 병렬 템퍼링 · 담금질 모사
• 유효표본크기 · 수치적분
• 힘장 · 구속 동역학 · 검증 및 확인
• GROMACS · LAMMPS

10. Footnotes[편집]

1. Torrie, G. M. & Valleau, J. P. (1977). “Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling”. J. Comput. Phys. 23, 187. 이름의 유래는 원래 제안한 편향 함수가 관심 영역 전체를 우산처럼 덮도록 넓게 퍼진 가중치였기 때문이다. 오늘날 표준이 된 “좁은 조화 창을 줄줄이 세우는” 방식은 원조 우산보다는 우산 여러 개를 이어 붙인 아케이드에 가까운데, 이름만 살아남았다. ↩

2. 리간드-수용체 결합 PMF에서 이 보정을 빼먹고 결합 자유에너지를 보고하는 사고가 지금도 주기적으로 일어난다. 게다가 결합 상수로 환산하려면 야코비안 위에 표준상태(1 M 기준) 보정까지 얹어야 해서, 보정 두 개를 다 챙기지 않으면 실험값과 몇 kcal/mol씩 어긋난다. PMF 곡선의 모양은 멀쩡한데 숫자만 틀리므로 그림만 봐서는 절대 안 걸린다. ↩

3. Kumar, S. 외 (1992). J. Comput. Chem. 13, 1011. 뿌리는 Ferrenberg & Swendsen (1989)의 다중 히스토그램 재가중이고, 그쪽은 원래 이징 모형 임계 현상용이었다. 스핀 격자에서 태어난 알고리즘이 신약 설계 파이프라인에서 굴러가고 있는 셈. ↩

4. 이 실패 양상에 붙은 별명이 “hidden barrier”다. 겉으로 드러난 좌표의 장벽은 우리가 스프링으로 눌러 없앴는데, 정작 시간을 잡아먹는 장벽은 우리가 보지도 않는 방향에 그대로 서 있다. 조명 아래에서만 열쇠를 찾는 술 취한 사람 비유가 이 분야만큼 잘 맞는 곳도 드물다. ↩