分子動力学(MD)計算で得られた座標データに対して、SciPyで主成分分析をしてみた。通常、座標のトラジェクトリーに対して最小二乗法によるフィッティングを行ってからPCAを行うのだが、ここで示すPCAの仮引数trjに渡すトラジェクトリーはフィッティングを前提にしてある。
実際のコードではフィッティングやら、ARPACKの利用やら、原子の重み付けの方法やらをオプションで引き渡して処理するようになっているのだが、ここでは内容を明確にするために、MDの一般的なPCAに絞ったコードを示す。
MD計算により10個の構造が出力される場合、N原子を持つ一構造のXYZ座標のリスト([N][3])をトラジェクトリー([10][N][3])として渡す(trj)。原子量は要素がN個のリストとなる(mass)。原子量のリストを渡さなければ原子の重み付けは行われず、Essential Dynamics (ED)の計算となる。
以下、コード。
from scipy import * def pca(trj, mass=[], num=2): """trj: 座標データ. / mass: 原子量のリスト(重み付けが必要であれば) / num: 求める主成分の数""" # 座標の平滑化. R = [mat(crd).A1 for crd in trj] NDIM, N = len(R), len(R[0]) # 平均構造の作成. R_mean = zeros(N) for i in range(NDIM): R_mean += R[i] R_mean /= NDIM # 各構造と平均構造の差. for i in range(NDIM): R[i] -= R_mean # 必要なら原子の重み付けを行う. if mass: sq_mass = [sqrt(m) for m in mass] for i in range(NDIM): for j in range(N): R[i][j] *= sq_mass[j // 3] # 分散共分散行列の作成. covR = zeros((N, N)) for R_i in R: covR += mat(R_i).T * mat(R_i) covR /= NDIM # 固有値・固有ベクトルの算出. (eigenval, eigenvec) = linalg.eig(array(covR)) axis = zip(eigenval, eigenvec.T) axis.sort(reverse=True) # 第(num)主成分まで算出. Q = [[(mat(axis[i][1]) * mat(R_i).T)[0, 0].real for i in range(num)] for R_i in R] return Q
追記(2009/7/13):
programming/python/scipy - PukiWikiに、より洗練されたPCAのコードが示されていたので、それを参考に修正してみた。ついでに気に入らない部分も変更した。やっていることは同じ。
from numpy import * def pca(trj, mass=[], num=2): """trj: 座標データ. / mass: 原子量のリスト(重み付けが必要であれば) / num: 求める主成分の数""" # 座標の平滑化. R = array([mat(crd).A1 for crd in trj]) # 平均構造の作成. R_mean = array(map(sum, R.T)) / len(R) # 各構造と平均構造の差. R -= R_mean # 必要なら原子の重み付けを行う. if mass: mass = array([sqrt(m) for m in mass for i in range(3)]) R *= mass # 分散共分散行列の作成. covR = dot(R.T, R) / len(R) # 固有値・固有ベクトルの算出. (eigenval, eigenvec) = linalg.eig(covR) axis = zip(eigenval, eigenvec.T) axis.sort(key=lambda x: x[0], reverse=True) # 第(num)主成分まで算出. Q = [[(mat(axis[i][1]) * mat(R_i).T)[0, 0].real for i in range(num)] for R_i in R] return Q
実際のコードではフィッティングやら、ARPACKの利用やら、原子の重み付けの方法やらをオプションで引き渡して処理するようになっているのだが、ここでは内容を明確にするために、MDの一般的なPCAに絞ったコードを示す。
MD計算により10個の構造が出力される場合、N原子を持つ一構造のXYZ座標のリスト([N][3])をトラジェクトリー([10][N][3])として渡す(trj)。原子量は要素がN個のリストとなる(mass)。原子量のリストを渡さなければ原子の重み付けは行われず、Essential Dynamics (ED)の計算となる。
以下、コード。
from scipy import * def pca(trj, mass=[], num=2): """trj: 座標データ. / mass: 原子量のリスト(重み付けが必要であれば) / num: 求める主成分の数""" # 座標の平滑化. R = [mat(crd).A1 for crd in trj] NDIM, N = len(R), len(R[0]) # 平均構造の作成. R_mean = zeros(N) for i in range(NDIM): R_mean += R[i] R_mean /= NDIM # 各構造と平均構造の差. for i in range(NDIM): R[i] -= R_mean # 必要なら原子の重み付けを行う. if mass: sq_mass = [sqrt(m) for m in mass] for i in range(NDIM): for j in range(N): R[i][j] *= sq_mass[j // 3] # 分散共分散行列の作成. covR = zeros((N, N)) for R_i in R: covR += mat(R_i).T * mat(R_i) covR /= NDIM # 固有値・固有ベクトルの算出. (eigenval, eigenvec) = linalg.eig(array(covR)) axis = zip(eigenval, eigenvec.T) axis.sort(reverse=True) # 第(num)主成分まで算出. Q = [[(mat(axis[i][1]) * mat(R_i).T)[0, 0].real for i in range(num)] for R_i in R] return Q
追記(2009/7/13):
programming/python/scipy - PukiWikiに、より洗練されたPCAのコードが示されていたので、それを参考に修正してみた。ついでに気に入らない部分も変更した。やっていることは同じ。
from numpy import * def pca(trj, mass=[], num=2): """trj: 座標データ. / mass: 原子量のリスト(重み付けが必要であれば) / num: 求める主成分の数""" # 座標の平滑化. R = array([mat(crd).A1 for crd in trj]) # 平均構造の作成. R_mean = array(map(sum, R.T)) / len(R) # 各構造と平均構造の差. R -= R_mean # 必要なら原子の重み付けを行う. if mass: mass = array([sqrt(m) for m in mass for i in range(3)]) R *= mass # 分散共分散行列の作成. covR = dot(R.T, R) / len(R) # 固有値・固有ベクトルの算出. (eigenval, eigenvec) = linalg.eig(covR) axis = zip(eigenval, eigenvec.T) axis.sort(key=lambda x: x[0], reverse=True) # 第(num)主成分まで算出. Q = [[(mat(axis[i][1]) * mat(R_i).T)[0, 0].real for i in range(num)] for R_i in R] return Q
コメント