経緯
今年の10月にOculus Quest 2が発売されましたね。VR機器が（比較的）安価に手に入るということだったので、これはと思って試しに購入し、遊んでみたところ期待通りの良いものでした。せっかくなので自分で何かアプリケーションを書いてみようと思い、普段分子動力学(以下MD)のコードに触れる機会が多いので、それなら分子を触って動かせるMDを書いてみようと思った次第です。
Oculusで動くアプリ開発は基本的にUnity上で行われているようです。UnityのOculus Integration がコントローラーでオブジェクトを掴んだりなんだりをよしなにやってくれるようなので、基本的にはUnityの仕様に乗ることでVRの煩雑さを回避していく方針で進めます。
※UnityやC#に触れたのは今回が初めてなので、回りくどい実装になっている部分もあるかと思います。
Lennard-jones流体を動かす
Lennard-jonesポテンシャル
まずシンプルなケースとして、Lennard-jones(LJ)流体を実装していきます。
Unityには物理衝突を計算するためのコンポーネントであるColliderが存在するので、これを使います。ただ、Unityビルトインの衝突処理は基本的に剛体間の衝突のサポートになっているので、この部分をLJポテンシャルでの衝突に上書きしていきます。
ポテンシャルの表式は下記です。
U(r) = 4 \epsilon \left[ \left(\frac{\sigma}{r} \right)^{12} - \left(\frac{\sigma}{r} \right)^6 \right]
$r$が二粒子間の距離、$\sigma$は衝突直径、$\epsilon$はポテンシャルの深さに相当します。粒子のペアについて$\sigma$と$\epsilon$を決める方法はいくつかあるのですが、今回はLorentz–Berthelot則を用いて、例えば粒子$i$と$j$について、下記のように決めています。
\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2} \\
\epsilon_{ij} = \sqrt{\epsilon_i \epsilon_j}
計算の際に必要なのは力なので、実際には下記の、ポテンシャル式の微分を実装します。
\frac{dU(r)}{dr} = \frac{24 \epsilon}{r} \left[ \left( \frac{\sigma}{r} \right)^6 - 2\left( \frac{\sigma}{r} \right)^{12} \right]
Lennard-jones相互作用の実装
数値積分法は取り敢えずUnityに委ねます。Unity内部ではPhysXという物理エンジンが使われているのですが、商用ソフトウェアということもあって詳細な仕様は公開されていないみたいです（？）。
https://gamedev.stackexchange.com/questions/79938/unitys-default-integration-method
大まかな実装方法としては、Unityの粒子に相当するGameObjectにRigidBody、SphereCollider、LJポテンシャル用のスクリプトを付加し、実際の衝突処理は無効にして、Colliderでは衝突の検出のみを行います。実際の衝突処理はLJポテンシャル用スクリプトのOnTriggerStay関数内で力を計算し、直接RigidBodyにかかる力に足しこむことで行います。OnTriggerStayは二つのColliderが重なっている間のみゲーム内時間で一定間隔に呼び出されるメソッドです。
上記を実装すると、LJポテンシャルのスクリプトは下記のようになります。
using System;
using UnityEngine;
public class LennardJonesParticle : MonoBehaviour
{
public float sphere_radius; // r
public float scaled_epsilon; // ε
private Rigidbody m_Rigidbody;
private SphereCollider m_SphereCollider;
private void OnTriggerStay(Collider other)
{
LennardJonesParticle other_lj = other.GetComponent<LennardJonesParticle>();
// （省略）
// 以下力の計算
Vector3 dist_vec = other.attachedRigidbody.position - transform.position;
float sigma = sphere_radius + other_lj.sphere_radius;
float epsilon = Mathf.Sqrt(scaled_epsilon * other_lj.scaled_epsilon);
float rinv = 1.0f / dist_vec.magnitude;
float r1s1 = sigma * rinv;
float r3s3 = r1s1 * r1s1* r1s1;
float r6s6 = r3s3 * r3s3;
float r12s12 = r6s6 * r6s6;
float derivative = 24.0f * epsilon * (r6s6 - 2.0f * r12s12) * rinv;
m_Rigidbody.AddForce(derivative * rinv * dist_vec);
}
// Initは最初のフレームが計算される前に処理されるようにする
internal void Init(float radius, float epsilon, float timescale)
{
// 各メンバ変数の初期化
//（省略）
m_Rigidbody = GetComponent<Rigidbody>();
m_SphereCollider = GetComponent<SphereCollider>();
// （省略）
// This radius mean cutoff radius
m_SphereCollider.radius = 1.25f; // カットオフ距離の設定
m_SphereCollider.isTrigger = true; // Colliderの衝突処理を無効にする
}
};
また、上で作成したGameObjectにOculus IntegrationのOVR Grabbableを付加し、シーン中のコントローラーオブジェクトにOVR Grabberを付加することで、VR空間で粒子を掴んで動かすことができるようになります。
初期速度の生成
初期速度はMaxwell-Boltzmann分布に従って生成します。分布の式は各成分について、
P(\dot{x}) = \sqrt{\frac{m}{2\pi k_B T})}\exp\left( - \frac{m \dot{x}^2}{2 k_B T} \right)
になりますが、Unityの乱数は一様乱数の生成しかサポートしていないので、一様乱数をこちらで正規分布に成形してやる必要があります。乱数の変換にはMarsaglia polar methodを用います。
Marsaglia poloar methodは、$-1 < x < 1$、$-1 < y < 1$、で分布する一様乱数$x,y$が$$0 < x^2 + y^2 < 1$$の条件を満たしているとき、正規分布に従う互いに独立な確率変数$u,v$を
u = x \sqrt{\frac{-2 \ln{(s)}}{s}}, \quad v = y \sqrt{\frac{-2\ln{(s)}}{s}}
で得るもので、Box-muller法に比べ、三角関数を使わないので幾分高速に計算できます。また、この方法では一度に二つの乱数を生成することになるので、乱数の変換器の中で生成した乱数をストックするようにします。
実装は以下のようになります。
using UnityEngine;
public class NormalizedRandom
{
private bool has_stock;
private float stock;
public NormalizedRandom()
{
has_stock = false;
}
public float Generate()
{
if (!has_stock)
{ // ストックがない場合は乱数を生成する
float u, v, s;
do
{
u = Random.value * 2 - 1;
v = Random.value * 2 - 1;
s = u * u + v * v;
} while ( s >= 1 || s == 0);
s = Mathf.Sqrt(-2.0f * Mathf.Log(s) / s);
stock = v * s;
has_stock = true;
return u * s;
}
else
{ // ストックがある場合はストックを返す
has_stock = false;
return stock;
}
}
}
その他の部分
境界条件には反射境界条件を使いました。また、系の温度、粒子の重さ、各粒子のLJポテンシャルについての各パラメータなどは下記のようなTOMLフォーマットのファイルで指定できるようにしました。
[[systems]]
attributes.temperature = 300.0
boundary_shape.lower = [-5.0, -5.0, -5.0]
boundary_shape.upper = [ 5.0, 5.0, 5.0]
particles = [
{m = 1.00, pos = [ 0.9084463866571024, 2.667970365022592, 0.6623741650618591]}, # particle index 1
{m = 1.00, pos = [-0.39914657537482867, 2.940257103317942, 3.5414659037905025]}, # particle index 2
#...
]
[[forcefields]]
[[forcefields.global]]
potential = "LennardJones"
parameters = [
{index = 0, sigma = 0.5, epsilon = 0.05},
{index = 1, sigma = 0.5, epsilon = 0.05},
#...
]
このフォーマットは、OSSの汎用MDシミュレーションソフトMjolnirのインプットのサブセットになっています 。
動作画面
実際に動かしてみるとこんな感じになります。
楽しいですね。
LJ流体のコードは以下のリポジトリで公開しています。Oculus Riftの環境がある方は、Questを開発者モードにして、README手順に従ってセットアップすれば（多分）プレイ可能です。
https://github.com/yutakasi634/LennardJonesFluidOnVR
タンパク質を動かす
Langevin dynamics
さて、LJ流体のようなシンプルな系を触れるようになると、タンパク質のようなもっと複雑な生体分子も触りたくなってきます。しかし、VRアプリケーションはリアルタイムで描画処理を行うので、素直に全原子分子動力学を行うと数十秒から数分の時間ではタンパク質に大きな動きは見られなくなってしまいます。そこで、今回はLangevin dynamicsを用いて溶媒を陰に扱い、計算を軽くします。
Langevin dynamicsの式は下記になります。
m\frac{d^2{\bf r}}{dt^2} = {\bf f} - m \gamma \frac{d{\bf r}}{dt} + m \xi
$m$は対象の粒子の重さ、$\bf r$は座標、$\bf f$はポテンシャルからの力、$\gamma$は摩擦係数、$\xi$は周囲の溶媒から受けるランダム力を表しています。$\xi$は
\left<\xi(t)\right> = 0, \quad \left< \xi(t)^2 \right> = \frac{2\gamma k_B T}{m}
の正規分布ですが、実際は離散時間での計算になるので、
\left< \xi(t)^2 \right> = \frac{2\gamma k_B T}{m \Delta t}
になります。
Langevin dynamicsの実装
基本的にはUnityによるアップデートの間に速度の減衰とランダム力を挟み込めばよいはずです。そこで、この処理用のGameObjectをシーン中に配置し、そこにLangevin dynamicsのスクリプトを付加して実装します。減衰項とランダム項からの力は、フレームレートに応じて呼び出されるUpdateメソッドではなく、ゲーム内時間一定間隔で呼び出されるFixedUpdateメソッドを用いて実装します。
実装は以下のようになります。
using System.Collections.Generic;
using UnityEngine;
public class UnderdampedLangevinManager : MonoBehaviour
{
private List<float> m_NoiseCoefs;
private List<float> m_ScaledGammas;
private List<Rigidbody> m_LJRigidbodies;
private NormalizedRandom m_NormalizedRandom;
private void FixedUpdate()
{
for (int part_idx = 0; part_idx < m_LJRigidbodies.Count; part_idx++)
{
Rigidbody ljrigid = m_LJRigidbodies[part_idx];
Vector3 accelerate = -m_ScaledGammas[part_idx] * ljrigid.velocity; // 減衰項の計算
Vector3 random_force = new Vector3(m_NormalizedRandom.Generate(),
m_NormalizedRandom.Generate(),
m_NormalizedRandom.Generate());
accelerate += m_NoiseCoefs[part_idx] * random_force; // ランダム力の項の計算
ljrigid.AddForce(accelerate, ForceMode.Acceleration);
}
}
internal void Init(float kb_scaled, float temperature,
List<GameObject> general_particles, float[] gammas, float timescale)
{
// メンバ変数の初期化（省略）
int particles_num = general_particles.Count;
m_LJRigidbodies = new List<Rigidbody>();
m_NoiseCoefs = new List<float>();
float invdt = 1.0f / Time.fixedDeltaTime;
for (int part_idx = 0; part_idx < particles_num; part_idx++)
{
Rigidbody ljrigid = general_particles[part_idx].GetComponent<Rigidbody>();
// 各粒子にかかるランダム力の項の係数を計算
float noise_coef
= Mathf.Sqrt(2.0f * m_ScaledGammas[part_idx] * kb_scaled * temperature * invdt / ljrigid.mass);
m_LJRigidbodies.Add(ljrigid);
m_NoiseCoefs.Add(noise_coef);
}
}
}
Off-lattice Goモデル
これで溶媒分子を用に計算する必要はなくなりましたが、それでもまだ、タンパク質を構成するすべての原子をそのまま粒子にして計算をすすめるのは大変です。そこで、タンパク質を各アミノ酸のCα原子の集まりで表現することで、さらなる粒子数の削減を図ります。このスケールでタンパク質を扱えるモデルとしてシンプルなものに、Off-lattice Goモデル(Clementi et al., 2000)があるので、こちらの実装を試みます。このモデルでは、タンパク質のフォールディングが構造依拠ポテンシャルを使って安定化されています。
Off-lattice Goモデルは下記表式で表されます。
E(\Gamma, \Gamma_0) = \sum_{\rm bonds} K_r \left( r - r_0\right)^2 +
\sum_{\rm angles} K_{\theta} \left( \theta - \theta_0 \right)
+ \sum_{\rm dihedral} K_{\phi}^{(n)}\left[ 1 + \cos(n \times \left( \phi - \phi_0 \right)) \right] \\
+ \sum_{i<j-3} \left\{ \epsilon(i,j) \left[ 5\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - 6\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{10} \right] \right\}
+ \epsilon_2(i,j)\left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12}
右辺第1~3項はそれぞれ結合長$r$、結合角$\theta$、二面角$\phi$に対してかかるポテンシャルです。$r_0$、$\theta_0$、$\phi_0$はそれぞれ対象のタンパク質の安定と考えられる構造（参照構造と呼び、例えば、X線結晶構造がとられているのならその構造）が使われます。
第4項は参照構造で掲載されているアミノ酸残基間の非共有結合のコンタクトを安定化させるためのポテンシャルで、第5項はアミノ酸残基の排除体積を表す項です。第5項はLJ流体におけるLJポテンシャルに相当しています。Off-lattice Goモデルでは排除体積を表す項が結合とコンタクトを形成しているアミノ酸残基のペア間にはかからないようになっています。したがって、参照構造で残基間にコンタクトが形成されない場合は$\epsilon(i,j) = {\rm constanct} > 0$で$\epsilon_2(i,j) = 0$、コンタクトが形成される場合は$\epsilon(i,j) = 0$で$\epsilon_2(i,j) = {\rm constanct} > 0$となります。$\rm constant$はパラメータです。
各項の係数$K_r, K_\theta, K_\phi^{(1)}, K_\phi^{(3)}$はそれぞれ、$K_r = 100 \epsilon$, $K_{\theta} = 20\epsilon$、$K_\phi^{(1)} = \epsilon$、$K_\phi^{(2)} = 0.5\epsilon$で決まります。$n = 1 \quad {\rm or} \quad 3$以外での$K_\phi^{(n)}$は$0$です。
※12/22日の記事執筆時点では、時間の都合上結合角と二面角を除いた以下の部分のみの実装になっています。
E(\Gamma, \Gamma_0) = \sum_{\rm bonds} K_r \left( r - r_0\right)^2
+ \sum_{i<j-3} \left\{ \epsilon(i,j) \left[ 5\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - 6\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{10} \right] \right\}
+ \epsilon_2(i,j)\left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12}
結合長ポテンシャルの実装
結合長に対してかかるポテンシャルは粒子ペアに対して定義されるので、それぞれを表現するクラスが粒子のペアを持ち、各ステップごとに力を計算していく実装にします。また、粒子間の排除体積の無視は、UnityのPhysics.IgnoreCollisionメソッドを用いて衝突を無視することで実現します。
実装は以下のようになります。
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public class HarmonicBondManager : MonoBehaviour
{
private List<float> m_V0s; // r_0
private List<float> m_ScaledKs; // K_r
private List<List<Rigidbody>> m_RigidPairs;
private void FixedUpdate()
{
for (int pair_idx = 0; pair_idx < m_RigidPairs.Count; pair_idx++)
{
List<Rigidbody> rigid_pair = m_RigidPairs[pair_idx];
Rigidbody rigid_first = rigid_pair[0];
Rigidbody rigid_second = rigid_pair[1];
Vector3 dist_vec = rigid_second.position - rigid_first.position;
Vector3 norm_vec = dist_vec.normalized;
float coef = 2.0f * m_ScaledKs[pair_idx] * (dist_vec.magnitude - m_V0s[pair_idx]);
rigid_first.AddForce(coef * norm_vec);
rigid_second.AddForce(-coef * norm_vec);
}
}
internal void Init(List<float> v0s, List<float> ks, List<List<Rigidbody>> rigid_pairs, float timescale)
{
// 各メンバ変数の初期化(省略)
// コリジョンを無視するペアの設定
foreach (List<Rigidbody> rigid_pair in m_RigidPairs)
{
Collider first_collider = rigid_pair[0].GetComponent<Collider>();
Collider second_collider = rigid_pair[1].GetComponent<Collider>();
Physics.IgnoreCollision(first_collider, second_collider);
}
}
}
その他の部分
コンタクトに対してかかるポテンシャルは、力の計算部分を変更するだけでロジックはほぼ同じになります。また、排除体積を表す項は、前節のLJポテンシャルの力の計算を変更するだけで実現することができます。
各ポテンシャル中のパラメータに関しては、LJ流体のケースと同様にTOMLフォーマットのインプットファイルを読み込ませることで指定できます。後述するリポジトリの\input以下にサンプルが置いてあります。
動作画面
上記実装で実際に動かしてみるとこんな感じになります。
パラメータ次第ではリフォールディング過程も部分的に見れたりします。コンタクトの強さと時間刻みを弄ったら割合取り回しやすくなった pic.twitter.com/lM2QBQzdGK
— 紙を破る (@yutakashi634) December 23, 2020
今回実装したシミュレータは以下のリポジトリで公開しています。
https://github.com/yutakasi634/CoarseGrainedMDonVR
LJ流体と同様に、Oculus Riftの環境がある方は、Questを開発者モードにして、README手順に従ってセットアップすれば（多分）プレイ可能です。