등온등압 앙상블

편집 역사 토론
계산물리 분자동역학 물리 마지막 수정: 2026-08-04 04:41:20

1. 개요[편집]

등온등압 앙상블
Isothermal-Isobaric Ensemble
다른 이름NPT 앙상블
고정량입자수 N, 압력 P, 온도 T
요동하는 양에너지 E, 부피 V
확률 가중exp[−β(E + PV)]
특성 함수깁스 자유에너지 G = −kBT ln Δ
주 구현바로스탯(MTK·파리넬로-라만), MC 부피 이동

실험은 대개 뚜껑 열린 플라스크에서 한다. 부피는 우리가 정하는 게 아니라 계가 정한다.

등온등압 앙상블(isothermal-isobaric ensemble)은 입자수 NN·압력 PP·온도 TT가 고정되고 에너지와 부피가 함께 요동치는 미시상태 분포로, 부피 VV인 배치가 나타날 확률이 e−β(E+PV)e^{-\beta(E + PV)}에 비례하는 앙상블이다. 줄여서 NPT라 부른다.

실험실 조건에 가장 가까운 앙상블이라서, 밀도·열팽창계수·용해 자유에너지처럼 “1 atm, 300 K에서의 값”을 재는 계산은 거의 전부 여기서 돌아간다. 정준 앙상블에서 상자 크기를 잘못 잡으면 액체가 은근히 과압축·과팽창된 상태로 굳어 버리는데, NPT는 그 부담을 계에 떠넘긴다.

범위 정리: 앙상블의 기본 개념과 요동-응답 관계는 정준 앙상블·미시정준 앙상블, 입자수를 여는 쪽은 대정준 앙상블, 온도를 붙잡는 알고리즘 비교는 서모스탯에 있다. 이 문서는 부피 자유도만 책임진다.

2. 부피 적분이 붙은 분배함수[편집]

정준 분배함수 Z(N,V,T)Z(N,V,T) 를 부피에 대해 한 번 더 적분하면 NPT 분배함수가 나온다.

Δ(N,P,T)=1V0∫0∞dV  e−βPVZ(N,V,T)\Delta(N,P,T) = \frac{1}{V_0}\int_0^\infty dV\; e^{-\beta PV} Z(N,V,T)

V0V_0 는 차원을 맞추는 임의의 기준 부피이고, 이 선택은 GG 에 상수만 더하므로 열역학적 극한에서 무해하다.1 좌표를 상자 한 변으로 나눈 환산좌표 s=r/L\mathbf{s} = \mathbf{r}/L 로 바꾸면 적분 영역이 부피와 무관한 단위 정육면체가 되고, 야코비안이 VNV^N 으로 떨어져 나온다.

Δ∝∫0∞dV VNe−βPV∫dsNe−βU(s;V)\Delta \propto \int_0^\infty dV\, V^N e^{-\beta PV} \int d\mathbf{s}^N e^{-\beta U(\mathbf{s};V)}

이 VNV^N 항이 나중에 몬테카를로 수용확률에 Nln⁡(V′/V)N\ln(V'/V) 로 등장하는 그 항이다. 특성 함수는 깁스 자유에너지다.

G(N,P,T)=−kBTln⁡Δ,⟨V⟩=(∂G∂P)N,T,S=−(∂G∂T)N,PG(N,P,T) = -k_B T \ln \Delta, \qquad \langle V \rangle = \left(\frac{\partial G}{\partial P}\right)_{N,T}, \qquad S = -\left(\frac{\partial G}{\partial T}\right)_{N,P}

부피 요동은 등온압축률과 직결된다. 정준 앙상블의 ⟨ΔE2⟩=kBT2CV\langle \Delta E^2\rangle = k_BT^2C_V 와 완전히 같은 구조다.

⟨ΔV2⟩=⟨V2⟩−⟨V⟩2=kBT ⟨V⟩ κT\langle \Delta V^2 \rangle = \langle V^2\rangle - \langle V\rangle^2 = k_B T\, \langle V\rangle\, \kappa_T

그래서 NPT 궤적 하나에서 압축률을 공짜로 뽑을 수 있다. 반대로 말하면, 부피 요동을 죽이는 바로스탯은 압축률을 계산할 수 없다.

3. 순간 압력은 왜 그렇게 요동치는가[편집]

압력은 온도처럼 미시상태에서 읽어내는 추정량이고, 그 표현이 비리얼 정리에서 나온다.

P=13V(∑imivi2+∑i<jrij⋅fij)\mathcal{P} = \frac{1}{3V}\left( \sum_i m_i v_i^2 + \sum_{i<j} \mathbf{r}_{ij}\cdot\mathbf{f}_{ij} \right)

앞이 이상기체(운동) 항, 뒤가 상호작용 비리얼 항이다. 문제는 두 항이 거의 상쇄된다는 것. 상온 액체 물의 밀도에서 운동 항만 계산하면 ρkBT≈1400\rho k_B T \approx 1400 bar 수준이고, 비리얼 항이 거의 같은 크기의 음수로 이를 깎아 남는 잔돈이 1 bar다. 즉 우리가 재려는 값은 큰 수 두 개의 차이이며, 각 항의 상대오차가 조금만 흔들려도 결과는 수백 bar씩 널뛴다.

그래서 원자 수천 개짜리 상자의 순간 압력이 ±수백~수천 bar로 진동하는 것은 지극히 정상이고, 이걸 보고 “발산했나” 싶어 바로스탯 시간상수를 조이면 부피 요동만 망가진다. 평균 압력이 목표값에 수렴하는 데 걸리는 시간도 온도보다 훨씬 길다.2 참고로 압력의 요동 크기에는 열용량 같은 깔끔한 상한이 없다 — 그 분산이 하이퍼비리얼이라는 또 다른 상관량에 걸려 있기 때문이다.

4. 바로스탯 — 무엇이 진짜 NPT인가[편집]

온도조절기와 마찬가지로, 부피조절기도 “평균 압력을 맞추는 것”과 “올바른 부피 분포를 만드는 것”은 별개 문제다.

  • 베렌센(약한 결합). 매 스텝 상자와 좌표에 스케일 인자를 곱해 압력을 목표값 쪽으로 지수적으로 끌어당긴다. μ=[1−(Δt/τP)κT(P0−P)]1/3\mu = [1 - (\Delta t/\tau_P)\kappa_T(P_0 - \mathcal{P})]^{1/3}. 평형화용으로는 무적이지만 부피 요동을 억눌러 올바른 NPT 분포가 아니다. 여기서 나온 압축률은 신뢰할 수 없다.
  • 안데르센 피스톤(1980). 부피를 질량 WW 를 가진 동역학 변수로 승격시킨 확장 라그랑지안 방식. 등방 스케일링만 가능하다. 제대로 된 서모스탯과 결합하면 정확한 NPT를 준다.
  • 파리넬로-라만(1981). 상자를 벡터 세 개를 담은 행렬 h\mathbf{h} 로 두고 그 성분 전체를 동역학 변수로 만든다. 상자 모양까지 바뀌므로 결정의 구조 상전이를 시뮬레이션에서 직접 재현할 수 있다 — 원 논문의 목적이 그거였다.
  • MTK(마르티나-토비아스-클라인, 1994). 노제-후버 체인과 셀 동역학을 함께 정리해, 유한 NN 에서도 정확한 NPT 위상공간 분포를 재현함이 증명된 운동방정식. 오늘날 결정론적 NPT의 표준이며 이전 형식에 빠져 있던 자유도 보정항이 여기서 채워졌다.
  • MC 바로스탯 / 확률적 셀 재조정. 부피 변경을 메트로폴리스 수용판정으로 처리하면(다음 절) 구성상 정확하다. 베렌센에 잡음항을 정확히 얹은 확률적 셀 재조정도 올바른 부피 분포를 준다.

정리하면 베렌센은 평형화 전용, 생산 구간은 MTK·파리넬로-라만·MC 계열이 국룰이다.

5. 셀 모양 자유도와 MC 부피 이동[편집]

상자를 얼마나 자유롭게 풀어줄지는 계에 따라 갈린다.

  • 등방(isotropic) — 세 변을 같은 비율로 스케일. 액체의 기본값이다. 액체는 전단응력을 못 버티므로 이방성 바로스탯을 걸면 상자가 아무 저항 없이 납작해지거나 길쭉해져 최소 이미지 조건을 위협한다.
  • 반등방(semi-isotropic) — xyxy 는 묶고 zz 는 따로. 지질 이중층 막처럼 면내와 법선 방향의 압력이 다른 계에 쓴다. 표면장력을 고정하는 변형도 여기 속한다.
  • 완전 유연(flexible / triclinic) — h\mathbf{h} 의 6개 독립 성분을 모두 푼다. 결정, 응력 하 고체, 압력 유도 상전이에 필수. 대신 상자가 심하게 기울면 정기적으로 격자를 다시 정규화해 줘야 한다.

몬테카를로 쪽은 오히려 단순명쾌하다. 부피를 V→V′V \to V' 로 바꾸고 모든 분자의 환산좌표를 그대로 유지한 뒤, 다음 확률로 수락한다.

acc=min⁡{1,  exp⁡[−β(ΔU+P ΔV)+Nln⁡V′V]}\mathrm{acc} = \min\left\{ 1,\; \exp\left[ -\beta\left(\Delta U + P\,\Delta V\right) + N \ln\frac{V'}{V} \right] \right\}

PΔVP\Delta V 는 압력이 하는 일이고, Nln⁡(V′/V)N\ln(V'/V) 가 앞서 나온 VNV^N 야코비안이다. VV 대신 ln⁡V\ln V 에서 균일하게 걷는 변형이 흔히 쓰이며, 그때는 계수가 N+1N+1 이 된다.3 여기엔 가상의 피스톤 질량도, 시간상수도 없다 — 상세균형만 지키면 분포가 보장된다.

마지막으로 실무 함정 하나. 구속 동역학이 걸린 계에서 원자 좌표를 그대로 스케일하면 안 된다. 고정된 결합 길이까지 늘어나 버리기 때문이다. 그래서 분자 무게중심만 스케일하고 내부 좌표는 강체로 옮기며, 압력도 그에 맞춰 분자 비리얼로 계산한다. 이걸 섞어 쓰면 밀도가 조용히 몇 퍼센트 틀어진다.4

6. 관련 문서[편집]

7. Footnotes[편집]

  1. 이 기준 부피를 뭘로 잡느냐로 논문 한 편 분량의 논쟁이 있었다. 결론은 ”GG 에 kBTln⁡(⋅)k_BT\ln(\cdot) 짜리 상수가 붙을 뿐이고, NN 이 커지면 시성량 대비 0으로 간다”는 것. 유한 클러스터에서는 얘기가 또 달라지지만, 그건 그때 가서 걱정하기로 다들 합의했다. ↩

  2. 그래서 NPT 로그를 보고 “압력이 -300 bar인데요?”라고 놀란 신입에게 선배가 하는 말이 정해져 있다. “1 ns 평균 내고 다시 와.” 실제로 밀도는 수십 ps만에 자리를 잡는데 순간 압력은 끝까지 널뛴다 — 밀도가 잘 맞으면 압력은 맞은 것이다. ↩

  3. ln⁡V\ln V 에서 걷는 이유는 부피가 몇 배씩 변하는 기체-액체 공존 계산에서 VV 균일 이동이 큰 부피 쪽에만 시간을 쓰기 때문이다. 계수가 NN 에서 N+1N+1 로 바뀌는 이 한 글자를 빠뜨린 코드가 세상에 은근히 많다. ↩

  4. 물 상자의 밀도가 실험값보다 2~3% 낮게 나오는 버그의 단골 원인. 힘장 탓을 하며 파라미터를 만지기 전에 바로스탯의 스케일링 단위가 원자인지 분자인지부터 확인하는 게 훨씬 싸게 먹힌다. ↩