Processing面白いなぁ。以前はhttp://proce55ing.orgが公式サイトだったことから、Proce55ingとしても知られている。Java仕様のビジュアルデザイン用プログラミング言語なのだが、とても簡単に扱うことができる。プログラミングでアートする感じ。
そこで今回はなんちゃって分子動力学計算プログラムを作ってみた。もちろん原子をリアルタイムに計算しながら表示させている。ただし、コードを短くするために分子ではなく疑似原子単体同士の相互作用にしている。そのため計算している相互作用はvdW力とクーロン力のみだ。また特に工夫をしているわけでもなく全原子同士のすべての相互作用をそのまま計算している。分子動力学計算の詳細についてはPEACHによる生体分子の分子動力学シミュレーション(1)原理と方法あたりを参考にすると良いと思う。
実際のプログラムを以下に示す。num_atomsで原子数を指定しており、初期配置を比較的密集させてその状態で力を計算させている。また、左上にエネルギーを表示している。初期の高エネルギー状態から低エネルギー状態に安定化するのが見て取れる。その他の詳しいことは実際にソースコードを読めば分かると思う。短いしね。
以下、ソースコード。
simple_md.pde
/** * Simple molecular dynamics by nox, 2008.10.21. * * Perform molecular dynamics for pseudo-atoms without internal * interactions, such as bonds, angles, dihedrals. Only vdW and * Coulomb interactions. */ final int num_atoms = 20; Atom[] atoms = new Atom[num_atoms]; void setup() { size(360, 360, P3D); noStroke(); // initialize positions of atoms int cnt = 0; int num_side = ceil(pow(num_atoms, 1.0 / 3.0)); for (int i = 0; i < num_side && cnt < num_atoms; i++) for (int j = 0; j < num_side && cnt < num_atoms; j++) for (int k = 0; k < num_side && cnt < num_atoms; k++) atoms[cnt++] = new Atom((i - num_side / 2) * 20, (j - num_side / 2) * 20, (k - num_side / 2) * 20); PFont font = createFont("Times New Roman", 32); textFont(font); } float calcDistance(int atom1, int atom2) { float distance, sum = 0.0; for (int i = 0; i < 3; i++) { distance = atoms[atom1].crd[i] - atoms[atom2].crd[i]; sum += distance * distance; } return sqrt(sum); } float calcDistanceXYZ(int atom1, int atom2, int xyz) { return atoms[atom1].crd[xyz] - atoms[atom2].crd[xyz]; } float calcEnergy(Atom[] atoms, int num_atoms) { float sum = 0.0; for (int i = 0; i < num_atoms; i++) sum += atoms[i].energy(); return sum / 2.0; } void draw() { background(0); lights(); float[] r = new float[num_atoms]; float[][] r_xyz = new float[num_atoms][3]; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < num_atoms; j++) { r[j] = calcDistance(i, j); for (int k = 0; k < 3; k++) r_xyz[j][k] = calcDistanceXYZ(i, j, k); } atoms[i].calculate(r, r_xyz, num_atoms); } for (int i = 0; i < num_atoms; i++) atoms[i].move(); text("Energy: " + calcEnergy(atoms, num_atoms), 10, 40); }
simple_atom.pde
class Atom { float[] crd, f, max_range; float ene; Atom() { this(0, 0, 0); } Atom(float x, float y, float z) { this.crd = new float[3]; this.f = new float[3]; this.crd[0] = x; this.crd[1] = y; this.crd[2] = z; for (int i = 0; i < 3; i++) this.f[i] = 0.0; this.max_range = new float[3]; this.max_range[0] = width / 2; this.max_range[1] = height / 2; this.max_range[2] = 200; } // Calculate energies and forces // Total energy = vdW + Coulomb // vdW // U = eps * [(Rij / r[i])^12 - 2 * (Rij / r[i]^6)] // F = -12 * eps / Rij * [(Rij / r[i])^13 - (Rij / r[i])^7] * r_xyz / r[i] // Coulomb // U = SUM_i>j qiqj / r[i] // F = SUM_j qiqj / r[i]^3 * r_xyz void calculate(float[] r, float[][] r_xyz, int num_atoms) { float[] delta_f = new float[3]; final float eps = 10.0; final float Rij = 10.0; final float q = 10.0; // qi*qj this.ene = 0.0; for (int i = 0; i < num_atoms; i++) { if (r[i] == 0.0) continue; final float R12 = pow(Rij / r[i], 12); final float R6 = pow(Rij / r[i], 6); this.ene += eps * (R12 - 2.0 * R6); this.ene += q / r[i]; for (int j = 0; j < 3; j++) { delta_f[j] += -(12 * eps / Rij) * (R12 * Rij / r[i] - R6 * Rij / r[i]) * r_xyz[i][j] / r[i]; delta_f[j] += q / pow(r[i], 3) * r_xyz[i][j]; } } for (int i = 0; i < 3; i++) this.f[i] += delta_f[i]; } float energy() { return this.ene; } void move() { pushMatrix(); for (int i = 0; i < 3; i++) { this.crd[i] += this.f[i]; if (this.crd[i] > max_range[i] || this.crd[i] < -max_range[i]) this.f[i] *= -1; } translate(max_range[0] + this.crd[0], max_range[1] + this.crd[1], -max_range[2] + this.crd[2]); sphere(20); popMatrix(); } }
そこで今回はなんちゃって分子動力学計算プログラムを作ってみた。もちろん原子をリアルタイムに計算しながら表示させている。ただし、コードを短くするために分子ではなく疑似原子単体同士の相互作用にしている。そのため計算している相互作用はvdW力とクーロン力のみだ。また特に工夫をしているわけでもなく全原子同士のすべての相互作用をそのまま計算している。分子動力学計算の詳細についてはPEACHによる生体分子の分子動力学シミュレーション(1)原理と方法あたりを参考にすると良いと思う。
実際のプログラムを以下に示す。num_atomsで原子数を指定しており、初期配置を比較的密集させてその状態で力を計算させている。また、左上にエネルギーを表示している。初期の高エネルギー状態から低エネルギー状態に安定化するのが見て取れる。その他の詳しいことは実際にソースコードを読めば分かると思う。短いしね。
以下、ソースコード。
simple_md.pde
/** * Simple molecular dynamics by nox, 2008.10.21. * * Perform molecular dynamics for pseudo-atoms without internal * interactions, such as bonds, angles, dihedrals. Only vdW and * Coulomb interactions. */ final int num_atoms = 20; Atom[] atoms = new Atom[num_atoms]; void setup() { size(360, 360, P3D); noStroke(); // initialize positions of atoms int cnt = 0; int num_side = ceil(pow(num_atoms, 1.0 / 3.0)); for (int i = 0; i < num_side && cnt < num_atoms; i++) for (int j = 0; j < num_side && cnt < num_atoms; j++) for (int k = 0; k < num_side && cnt < num_atoms; k++) atoms[cnt++] = new Atom((i - num_side / 2) * 20, (j - num_side / 2) * 20, (k - num_side / 2) * 20); PFont font = createFont("Times New Roman", 32); textFont(font); } float calcDistance(int atom1, int atom2) { float distance, sum = 0.0; for (int i = 0; i < 3; i++) { distance = atoms[atom1].crd[i] - atoms[atom2].crd[i]; sum += distance * distance; } return sqrt(sum); } float calcDistanceXYZ(int atom1, int atom2, int xyz) { return atoms[atom1].crd[xyz] - atoms[atom2].crd[xyz]; } float calcEnergy(Atom[] atoms, int num_atoms) { float sum = 0.0; for (int i = 0; i < num_atoms; i++) sum += atoms[i].energy(); return sum / 2.0; } void draw() { background(0); lights(); float[] r = new float[num_atoms]; float[][] r_xyz = new float[num_atoms][3]; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < num_atoms; j++) { r[j] = calcDistance(i, j); for (int k = 0; k < 3; k++) r_xyz[j][k] = calcDistanceXYZ(i, j, k); } atoms[i].calculate(r, r_xyz, num_atoms); } for (int i = 0; i < num_atoms; i++) atoms[i].move(); text("Energy: " + calcEnergy(atoms, num_atoms), 10, 40); }
simple_atom.pde
class Atom { float[] crd, f, max_range; float ene; Atom() { this(0, 0, 0); } Atom(float x, float y, float z) { this.crd = new float[3]; this.f = new float[3]; this.crd[0] = x; this.crd[1] = y; this.crd[2] = z; for (int i = 0; i < 3; i++) this.f[i] = 0.0; this.max_range = new float[3]; this.max_range[0] = width / 2; this.max_range[1] = height / 2; this.max_range[2] = 200; } // Calculate energies and forces // Total energy = vdW + Coulomb // vdW // U = eps * [(Rij / r[i])^12 - 2 * (Rij / r[i]^6)] // F = -12 * eps / Rij * [(Rij / r[i])^13 - (Rij / r[i])^7] * r_xyz / r[i] // Coulomb // U = SUM_i>j qiqj / r[i] // F = SUM_j qiqj / r[i]^3 * r_xyz void calculate(float[] r, float[][] r_xyz, int num_atoms) { float[] delta_f = new float[3]; final float eps = 10.0; final float Rij = 10.0; final float q = 10.0; // qi*qj this.ene = 0.0; for (int i = 0; i < num_atoms; i++) { if (r[i] == 0.0) continue; final float R12 = pow(Rij / r[i], 12); final float R6 = pow(Rij / r[i], 6); this.ene += eps * (R12 - 2.0 * R6); this.ene += q / r[i]; for (int j = 0; j < 3; j++) { delta_f[j] += -(12 * eps / Rij) * (R12 * Rij / r[i] - R6 * Rij / r[i]) * r_xyz[i][j] / r[i]; delta_f[j] += q / pow(r[i], 3) * r_xyz[i][j]; } } for (int i = 0; i < 3; i++) this.f[i] += delta_f[i]; } float energy() { return this.ene; } void move() { pushMatrix(); for (int i = 0; i < 3; i++) { this.crd[i] += this.f[i]; if (this.crd[i] > max_range[i] || this.crd[i] < -max_range[i]) this.f[i] *= -1; } translate(max_range[0] + this.crd[0], max_range[1] + this.crd[1], -max_range[2] + this.crd[2]); sphere(20); popMatrix(); } }
コメント