たとえば、水：
1つの分子内の水素と酸素が、同じ酸素と隣接する分子の水素とはまったく異なる方法で相互作用することは明らかです。したがって、分子内相互作用と分子間相互作用は区別されます。分子間相互作用は、以前の記事で説明した短距離およびクーロンペアポテンシャルによって指定できます。ここでは、分子内に焦点を当てます。
最も一般的なタイプの分子内相互作用は、化学的（原子価）結合です。化学結合は、結合した原子間の距離に対する電位エネルギーの機能依存性によって、つまり、実際には同じペア電位によって設定されます。ただし、通常のペアポテンシャルとは異なり、この相互作用は特定のタイプの粒子ではなく、特定のペアの粒子に対して（インデックスによって）指定されます。化学結合電位の最も一般的な関数形式は、高調波電位です。
rは粒子間の距離であり、kは結合スティフネス定数であり、r 0は、平衡結合長です。とモールスの可能性：
ここで、Dはポテンシャルウェルの深さであり、パラメーターαはポテンシャルウェルの幅を特徴づけます。
次のタイプの分子内相互作用は結合角です。タイトルの写真をもう一度見てみましょう。静電気の力が水素イオン間の最大距離を提供すると想定されていたので、なぜ分子が角で描かれているのですか？これは180°に等しい角度HOHに対応しますか？事実、すべてが図に描かれているわけではありません。学校の化学コースから、酸素にはさらに2つの孤立した電子ペアがあり、その相互作用によって角度が歪むことを思い出すことができます。
古典的な分子ダイナミクスでは、電子や電子雲などのオブジェクトは通常導入されないため、「正しい」角度をシミュレートするために原子価角ポテンシャルが使用されます。3つの粒子の座標に対する潜在的なエネルギーの機能的依存性。このような最も便利なポテンシャルの1つは、ハーモニックコサインです。
θは、粒子のトリプレットによって形成される角度であり、kは剛性定数であり、θ 0は平衡角です。ねじれ角
など、より高次の分子内電位がありますが、それらは結合角よりもさらに人工的です。
所定のインデックスを持つ粒子間に相互作用を追加するのは簡単です。リンクの場合、リンクされた粒子のインデックスと相互作用のタイプを含む配列を格納します。各スレッドに処理用の独自のリンク範囲を提供します。これで完了です。同様に結合角も。したがって、すぐにタスクを複雑にします。化学結合と実行時結合角度を作成/削除する機能を追加します。これはすぐに私たちを古典的な分子ダイナミクスの面から外し、可能性の新しい地平を開きます。それ以外の場合は、特に無料で配布されるため、LAMMPS、DL_POLY、GROMACSなどの既存のパッケージから何かをダウンロードするだけで済みます。
さて、いくつかのコードについて。メイン構造に適切なフィールドを追加しましょう：
//bonds:
int nBond; // number of bonds
int mxBond; // maximal number of bonds
int4* bonds; // array of bonds
int* nbonds; // count of bond for a given atom
int* neighToBind; // a neighbor of a given atom for binding
int* canBind; // flags that atom[iat] can be bind
int* r2Min; // distances for the nearest neighbor (used for binding)
int* parents; // indexes of one of the atom bonded with a given
cudaBond* bondTypes;
int** def_bonds; // array[nSpec][nSpec] of default bond types
int** bindBonds; // array[nSpec][nSpec] bond types created by binding
float** bindR2; // square of binding distance [nSpec][nSpec]
//angles:
int nAngle; // number of angles
int mxAngle;
int4* angles; // array of angles
int* nangles; // number of angles for given atom
int* oldTypes;
cudaAngle* angleTypes;
int* specAngles; // [nSp] angle type formed by given species
リンクと角度の数は可変ですが、ほとんどの場合、メモリを過剰に割り当てないように、可能な最大数を見積もり、メモリをすぐに最大に割り当てることができます。フィールドnBondとmxBondは、それぞれ現在のリンク数と最大値を意味します。結合配列には、結合する原子のインデックス、結合のタイプ、および結合の作成時刻が含まれます（平均結合寿命などの統計に突然関心がある場合）。角度配列には、結合角度と結合角度のタイプを構成する3つの原子のインデックスが格納されます。bondTypesとangleTypesアレイが可能結合電位と角度の特性を含むことになります。それらの構造は次のとおりです。
struct cudaBond
{
int type; // potential type
int spec1, spec2; // type of atoms that connected by this bond type
int new_type[2]; // bond type after mutation
int new_spec1[2], new_spec2[2];
int mxEx, mnEx; // flags: maximum or minimum of bond length exists
float p0, p1, p2, p3, p4; // potential parameters
float r2min, r2max; // square of minimal and maximal bond length
float (*force_eng)(float r2, float r, float &eng, cudaBond *bond); // return energy
int count; // quantity of such bonds
float rSumm; // summ of lentghs (for mean length calculation)
int rCount; // number of measured lengths (for mean length calculation)
int ltSumm, ltCount; // for calculation of lifetime
};
struct cudaAngle
{
int type; // potential type
float p0, p1, p2; // potential parameters
void (*force_eng)(int4* angle, cudaAngle* type, cudaMD* md, float& eng);
};
タイプ フィールドは、ポテンシャルの関数形式を定義します。new_type、new_spec1、およびnew_specは、結合タイプのインデックスであり、結合が変更された後（切断または別のタイプの結合に変わる）に結合される原子のタイプです。これらのフィールドは、2つの要素を持つ配列として表されます。 1つ目は長さがr2min1 / 2より短くなる状況に対応し、2つ目はr2max1 / 2を超える状況に対応します。..。アルゴリズムの最も難しい部分は、すべての結合のプロパティを適用することです。これは、結合の切断と変換の可能性、および他のフローが隣接する結合を切断して結合原子のタイプを変更する可能性があるという事実を考慮したものです。同じ水の例を使って説明しましょう。最初、分子は電気的に中性であり、化学結合は水素と酸素に共通の電子によって形成されます。大まかに言えば、水素原子と酸素原子の電荷はゼロであると言えます（実際、電子密度は酸素にシフトしているため、水素、δ+、および酸素-2δ-には小さなプラスがあります）。私たちが結合を断ち切ると、酸素は最終的にそれ自体のために電子を取り、水素はそれを与えます。得られた粒子は、Hであり+及びO - 。全部で5種類の粒子が得られますが、慣例的にH、H +、O、Oと呼びましょう。-、O 2 - 。後者は、水分子から両方の水素を分離すると形成されます。従って、反応：
H 2 O - > H + + OH -
及び
OH - - > H + + O 2 - 。
化学の専門家は、水のための標準的な条件下で、分解の第一段階は、実質的に10の平衡で、唯一の分子アウト（実装されていないことを私に修正します7イオンに解離し、それでも書かれたとおりではありません）。しかし、アルゴリズムの説明のために、そのようなスキームは実例となるでしょう。 1つのストリームが水分子の1つの結合を処理し、別のストリームが同じ分子の2番目の結合を処理するとします。そして、たまたま両方の結びつきを断ち切る必要がありました。次いで、一つのストリームは、H原子に変換する必要があり+及びO - 、及びHに第二+及びO 2 - 。フローが同時にこれを行う場合は、手順の初めの時点で、酸素はO状態であり、両方のフローがOに変換- 、間違っていました。どういうわけかそのような状況を防ぐ必要があります。化学結合を処理する関数のブロック図：
現在のアトムのタイプが接続のタイプに対応しているかどうかを確認し、対応していない場合は、デフォルトのタイプのテーブルから取得し（事前にコンパイルする必要があります）、アトム間の距離の2乗（r 2）を決定し、接続が最大長または最小長を意味する場合は、それが出ていないかどうかを確認します私たちがこれらの境界を超えているかどうか。そうした場合、接続のタイプを変更するか削除する必要があり、どちらの場合もアトムのタイプを変更します。このために、atomicCAS関数が使用されます-現在のタイプのアトムを本来あるべきタイプと比較し、この場合は新しいタイプに置き換えます。アトムのタイプがすでに別のスレッドによって変更されている場合は、最初に戻ってリンクタイプをオーバーライドします。最悪のシナリオは、最初のアトムのタイプを変更できたが、2番目のアトムは変更できなかった場合です。最初のアトムを変更した後、他のスレッドがすでに何かを実行している可能性があるため、ロールバックするには遅すぎます。抜け道は何ですか？最初に取った接続ではなく、別のタイプの接続を切断/変更しているふりをすることをお勧めします。最初の最初のアトムと変更された2番目のアトムの間にどのタイプの接続が必要であったかを見つけ、当初の予想と同じルールに従って処理します。この場合、原子のタイプが再び変更された場合は、同じスキームを再度使用します。ここでは暗黙的に暗示されていますが、新しいタイプの結合は同じ特性を持っていること-それは同じ長さで切断するなど、そして切断中に形成される粒子は必要に応じてあります。この情報はユーザーが触れているので、私たちは責任を私たちのプログラムから彼に移します。彼はすべてを正しく設定する必要があります。コード：
__global__ void apply_bonds(int iStep, int bndPerBlock, int bndPerThread, cudaMD* md)
{
int def;
int id1, id2; // atom indexes
int old, old_spec2, spec1, spec2, new_spec1, new_spec2; // atom types
int new_bond_type;
int save_lt, need_r, loop; // flags to save lifetime, to need to calculate r^2 and to be in ‘while’ loop
int mnmx; // flag minimum or maximum
int action; // flag: 0 - do nothing, 1 - delete bond, 2 - transform bond
cudaBond *old_bnd, *cur_bnd; // old bond type, current bond type
float dx, dy, dz, r2, r;
float f, eng = 0.0f;
__shared__ float shEng;
#ifdef DEBUG_MODE
int cnt; // count of change spec2 loops
#endif
if (threadIdx.x == 0)
{
shEng = 0.0f;
}
__syncthreads();
int id0 = blockIdx.x * bndPerBlock + threadIdx.x * bndPerThread;
int N = min(id0 + bndPerThread, md->nBond);
int iBnd;
for (iBnd = id0; iBnd < N; iBnd++)
if (md->bonds[iBnd].z) // the bond is not broken
{
// atom indexes
id1 = md->bonds[iBnd].x;
id2 = md->bonds[iBnd].y;
// atom types
spec1 = md->types[id1];
spec2 = md->types[id2];
old_bnd = &(md->bondTypes[md->bonds[iBnd].z]);
cur_bnd = old_bnd;
save_lt = 0;
need_r = 1;
loop = 1;
#ifdef DEBUG_MODE
cnt = 0;
#endif
if ((cur_bnd->spec1 == spec1)&&(cur_bnd->spec2 == spec2))
{
//ok
}
else
if ((cur_bnd->spec1 == spec2) && (cur_bnd->spec2 == spec1))
{
invert_bond(id1, id2, spec1, spec2, &(md->bonds[iBnd]));
//... then ok
}
else // atom types do not correspond to bond types
{
save_lt = 1;
}
// end initial stage
while (loop)
{
if (save_lt)
{
def = md->def_bonds[spec1][spec2];
if (def == 0) // these atom types do not form a bond
{
#ifdef DEBUG_MODE
printf("probably, something goes wrong\n");
#endif
action = 1; // delete
break;
}
else
{
//! change bond type and go on
if (def < 0)
{
invert_bond(id1, id2, spec1, spec2, &(md->bonds[iBnd]));
def = -def;
}
md->bonds[iBnd].z = def;
cur_bnd = &(md->bondTypes[def]);
}
} // end if (save_lt)
// calculate distance (only once)
if (need_r)
{
dx = md->xyz[id1].x - md->xyz[id2].x;
dy = md->xyz[id1].y - md->xyz[id2].y;
dz = md->xyz[id1].z - md->xyz[id2].z;
delta_periodic(dx, dy, dz, md);
r2 = dx * dx + dy * dy + dz * dz;
need_r = 0;
}
action = 0; // 0 - just cultivate bond 1 - delete bond 2 - transform bond
if ((cur_bnd->mxEx) && (r2 > cur_bnd->r2max))
{
mnmx = 1;
if (cur_bnd->new_type[mnmx] == 0) // delete bond
action = 1;
else
action = 2; // modify bond
}
else if ((cur_bnd->mnEx) && (r2 < cur_bnd->r2min))
{
mnmx = 0;
action = 2; // at minimum only bond modification possible
}
// end select action
// try to change atom types (if needed)
if (action)
{
save_lt = 1;
new_spec1 = cur_bnd->new_spec1[mnmx];
new_spec2 = cur_bnd->new_spec2[mnmx];
//the first atom
old = atomicCAS(&(md->types[id1]), spec1, new_spec1);
if (old != spec1)
{
spec1 = old;
spec2 = md->types[id2]; // refresh type of the 2nd atom
// return to begin of the ‘while’ loop
}
else // types[id1] have been changed
{
#ifdef USE_NEWANG // save changes in atom type
atomicCAS(&(md->oldTypes[id1]), -1, spec1);
#endif
old_spec2 = spec2;
while ((old = atomicCAS(&(md->types[id2]), old_spec2, new_spec2)) != old_spec2)
{
//! the worst variant: this thread changes atom 1, other thread changes atom 2
// imagine that we had A-old bond with the same behavior
def = md->def_bonds[spec1][old];
#ifdef DEBUG_MODE
if (def == 0)
{
printf("UBEH[001]: in apply_bonds, change atom types. There are no bond types between Species[%d] and [%d]\n", spec1, old);
break;
}
#endif
if (def < 0) // spec1 -> new_spec2 spec2 -> newSpec1
{
cur_bnd = &(md->bondTypes[-def]);
new_spec2 = cur_bnd->new_spec1[mnmx];
}
else // direct order
{
cur_bnd = &(md->bondTypes[def]);
new_spec2 = cur_bnd->new_spec2[mnmx];
}
old_spec2 = old;
#ifdef DEBUG_MODE
cnt++;
if (cnt > 10)
{
printf("UBEH[002]: too many atempst to change spec2 = %d\n", spec2);
break;
}
#endif
}
#ifdef USE_NEWANG // save changes in atom type
atomicCAS(&(md->oldTypes[id2]), -1, spec2);
#endif
loop = 0;
}
//end change types
} // end if (action)
else
loop = 0; // action == 0, out of cycle
} // end while(loop)
if (action == 2)
{
new_bond_type = cur_bnd->new_type[mnmx];
if (new_bond_type < 0)
{
md->bonds[iBnd].x = id2;
md->bonds[iBnd].y = id1;
new_bond_type = -new_bond_type;
}
md->bonds[iBnd].z = new_bond_type;
cur_bnd = &(md->bondTypes[new_bond_type]);
}
// perform calculations and save mean bond length
if (action != 1) // not delete
{
r = sqrt(r2);
f = cur_bnd->force_eng(r2, r, eng, cur_bnd);
atomicAdd(&(md->frs[id1].x), f * dx);
atomicAdd(&(md->frs[id2].x), -f * dx);
atomicAdd(&(md->frs[id1].y), f * dy);
atomicAdd(&(md->frs[id2].y), -f * dy);
atomicAdd(&(md->frs[id1].z), f * dz);
atomicAdd(&(md->frs[id2].z), -f * dz);
atomicAdd(&(cur_bnd->rSumm), r);
atomicAdd(&(cur_bnd->rCount), 1);
}
else //delete bond
{
// decrease the number of bonds for atoms
atomicSub(&(md->nbonds[id1]), 1);
atomicSub(&(md->nbonds[id2]), 1);
md->bonds[iBnd].z = 0;
// change parents
exclude_parents(id1, id2, md);
}
if (save_lt)
{
keep_bndlifetime(iStep, &(md->bonds[iBnd]), old_bnd);
if (action != 1) // not delete
atomicAdd(&(cur_bnd->count), 1);
atomicSub(&(old_bnd->count), 1);
}
} // end main loop
// split energy to shared and then to global memory
atomicAdd(&shEng, eng);
__syncthreads();
if (threadIdx.x == 0)
atomicAdd(&(md->engBond), shEng);
}
コードでは、プリプロセッサディレクティブを使用して、ユーザーの監視が原因で発生する可能性のある状況をチェックできるようにしました。それらをオフにして、パフォーマンスを高速化できます。この関数は上記のスキームを実装しますが、このスレッドが担当するリンクの範囲を実行するもう1つのループにラップされます。以下、結合タイプの識別子は負にすることができます。これは、結合内の原子の順序を逆にする必要があることを意味します（たとえば、OH結合とHO結合は同じ結合ですが、アルゴリズムでは順序が重要です。これを示すには、反対のインデックスを使用します。 sign）、invert_bond関数は、説明するのが簡単すぎます。Delta_periodic関数周期的な境界条件を適用して、差異を調整します。結合だけでなく結合角度も変更する必要がある場合（USE_NEWANGディレクティブ）、タイプを変更した原子をマークする必要があります（これについては後で詳しく説明します）。同じ原子の結合による再結合を除外するために、親配列はデータに関連付けられた原子の1つのインデックスを格納します（このセーフティネットはすべての場合に機能するわけではありませんが、私の場合はそれで十分です）。ある種の接続を切断した場合、対応するアトミックインデックスを親配列から削除する必要があります。これは、exclude_parents関数によって行われます。
__device__ void exclude_parents(int id1, int id2, cudaMD* md)
// exclude id1 and id2 from parents of each other (if they are)
// and seek other parents if able
{
// flags to clear parent
int clear_1 = 0;
int clear_2 = 0;
int i, flag;
if (md->parents[id1] == id2)
clear_1 = 1;
if (md->parents[id2] == id1)
clear_2 = 1;
i = 0;
while ((i < md->nBond) && (clear_1 || clear_2))
{
if (md->bonds[i].z != 0)
{
flag = 0;
if (clear_1)
{
if (md->bonds[i].x == id1)
{
md->parents[id1] = md->bonds[i].y;
flag = 1;
}
else if (md->bonds[i].y == id1)
{
md->parents[id1] = md->bonds[i].y;
flag = 1;
}
if (flag)
{
clear_1 = 0;
i++;
continue;
}
}
if (clear_2)
{
if (md->bonds[i].x == id2)
{
md->parents[id2] = md->bonds[i].y;
flag = 1;
}
else if (md->bonds[i].y == id2)
{
md->parents[id2] = md->bonds[i].y;
flag = 1;
}
if (flag)
{
clear_2 = 0;
i++;
continue;
}
}
}
i++;
}
// be on the safe side
if (clear_1)
md->parents[id1] = -1;
if (clear_2)
md->parents[id2] = -1;
}
残念ながら、この関数はリンクの配列全体を実行します。リンクを処理および削除する方法を学びました。次に、リンクを作成する方法を学ぶ必要があります。次の関数は、化学結合を形成するのに適した原子をマークします。
__device__ void try_to_bind(float r2, int id1, int id2, int spec1, int spec2, cudaMD *md)
{
int r2Int; // (int)r2 * const
// save parents to exclude re-linking
if (md->parents[id1] == id2)
return;
if (md->parents[id2] == id1)
return;
if (md->bindBonds[spec1][spec2] != 0)
{
if (r2 < md->bindR2[spec1][spec2])
{
r2Int = (int)(r2 * 100);
if (atomicMin(&(md->r2Min[id1]), r2Int) > r2Int) // replace was sucessfull
{
md->neighToBind[id1] = id2 + 1; // as 0 is reserved for no neighbour
md->canBind[id1] = 1;
}
// similar for the second atom
if (atomicMin(&(md->r2Min[id2]), r2Int) > r2Int) // replace was sucessfull
{
md->neighToBind[id2] = id1 + 1; // as 0 is reserved for no bind
md->canBind[id2] = 1;
}
}
}
}
bindBonds マトリックスは、これらのタイプの原子が結合を形成できるかどうか、および形成できる場合はどれかに関する情報を格納します。bindR2マトリックスは、結合に必要な原子間の最大距離を格納します。すべての条件が良好である場合、他の隣接する原子が結合に適しているかどうかを確認しますが、より近くなります。
ネイバーへの最も近い距離に関する情報はr2Min配列に格納されます（便宜上、配列はint型であり、値は定数100を掛けて変換されます）。検出されたネイバーが最も近いネイバーである場合、neighToBind配列にそのインデックスを記憶し、canBindフラグを設定します..。インデックスの更新に移っている間に、別のスレッドが最小値を上書きしてしまうという本当の危険がありますが、これはそれほど重要ではありません。この関数は、最初の部分で説明した、cell_listやall_pairなど、原子のペアをトラバースする関数で呼び出すことをお勧めします。バインディング自体：
__global__ void create_bonds(int iStep, int atPerBlock, int atPerThread, cudaMD* md)
// connect atoms which are selected to form bonds
{
int id1, id2, nei; // neighbour index
int btype, bind; // bond type index and bond index
cudaBond* bnd;
int spec1, spec2; // species indexes
int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;
int N = min(id0 + atPerThread, md->nAt);
int iat;
for (iat = id0; iat < N; iat++)
{
nei = md->neighToBind[iat];
if (nei) // neighbour exists
{
nei--; // (nei = spec_index + 1)
if (iat < nei)
{
id1 = iat;
id2 = nei;
}
else
{
id1 = nei;
id2 = iat;
}
// try to lock the first atom
if (atomicCAS(&(md->canBind[id1]), 1, 0) == 0) // the atom is already used
continue;
// try to lock the second atom
if (atomicCAS(&(md->canBind[id2]), 1, 0) == 0) // the atom is already used
{
// unlock the first one back
atomicExch(&(md->canBind[id1]), 1);
continue;
}
// create bond iat-nei
bind = atomicAdd(&(md->nBond), 1);
#ifdef DEBUG_MODE
if (bind >= md->mxBond)
{
printf("UBEH[003]: Exceed maximal number of bonds, %d\n", md->mxBond);
}
#endif
spec1 = md->types[id1];
spec2 = md->types[id2];
#ifdef USE_NEWANG // save that we have changed atom type
atomicCAS(&(md->oldTypes[id1]), -1, spec1);
atomicCAS(&(md->oldTypes[id2]), -1, spec2);
#endif
btype = md->bindBonds[spec1][spec2];
if (btype < 0)
{
// invert atoms order
md->bonds[bind].x = id2;
md->bonds[bind].y = id1;
md->bonds[bind].z = -btype;
bnd = &(md->bondTypes[-btype]);
// change atom types according the formed bond
md->types[id1] = bnd->spec2;
md->types[id2] = bnd->spec1;
}
else
{
md->bonds[bind].x = id1;
md->bonds[bind].y = id2;
md->bonds[bind].z = btype;
bnd = &(md->bondTypes[btype]);
// change atom types according the formed bond
md->types[id1] = bnd->spec1;
md->types[id2] = bnd->spec2;
}
atomicAdd((&bnd->count), 1);
md->bonds[bind].w = iStep; // keep time of the bond creation for lifetime calculation
atomicAdd(&(md->nbonds[id1]), 1);
atomicAdd(&(md->nbonds[id2]), 1);
// replace parents if none:
atomicCAS(&(md->parents[id1]), -1, id2);
atomicCAS(&(md->parents[id2]), -1, id1);
}
}
// end loop by atoms
}
この関数は、1つのアトムをブロックし、次に2番目のアトムをブロックし、両方のブロックに成功した場合は、それらの間に接続を作成します。関数の最初に、一方のスレッドがペアの最初のアトムをブロックし、もう一方のスレッドが同じペアの2番目のアトムをブロックする状況を除外するために、アトミックインデックスが並べ替えられます。両方のスレッドが最初のチェックに合格し、2番目のチェックに失敗します。その結果、どちらも接続はそれらを作成しません。そして最後に、apply_bonds関数で削除のマークを付けたリンクを削除する必要があります。
__global__ void clear_bonds(cudaMD* md)
// clear bonds with .z == 0
{
int i = 0;
int j = md->nBond - 1;
while (i < j)
{
while ((md->bonds[j].z == 0) && (j > i))
j--;
while ((md->bonds[i].z != 0) && (i < j))
i++;
if (i < j)
{
md->bonds[i] = md->bonds[j];
md->bonds[j].z = 0;
i++;
j--;
}
}
if ((i == j) && (md->bonds[i].z == 0))
md->nBond = j;
else
md->nBond = j + 1;
}
「無効化された」リンクを配列の最後に移動し、実際のリンク数を減らすだけです。残念ながら、コードはシリアルですが、並列化することで具体的な効果が得られるかどうかはわかりません。cudaBond構造のforce_engフィールドで示される、実際の結合エネルギーと原子への力を計算する関数はまだ省略されていますが、最初のセクションで説明したペアポテンシャルの関数と完全に類似しています。
原子価角
原子価角を使用して、アルゴリズムと関数を簡単にするためのいくつかの仮定を紹介します。その結果、それらは原子価結合の場合よりもさらに単純になります。まず、結合角のパラメータは3つの原子すべてに依存する必要がありますが、ここでは、結合角のタイプがその頂点の原子のみを決定すると仮定します。コーナーを形成/削除するための最も単純なアルゴリズムを提案します。原子のタイプを変更するたびに、対応する配列oldTypes []にこの事実を覚えています。配列のサイズは原子の数と同じで、最初は-1で埋められます。関数がアトムのタイプを変更すると、-1が元のタイプのインデックスに置き換えられます。タイプを変更したすべての原子について、すべての結合角度を削除し、この原子のすべての結合を実行して、対応する角度を追加します。
__global__ void refresh_angles(int iStep, int atPerBlock, int atPerThread, cudaMD *md)
// delete old angles and create new ones for atoms which change their type
{
int i, j, n, t, ang;
int nei[8]; // bonded neighbors of given atom
int cnt;
int id0 = blockIdx.x * atPerBlock + threadIdx.x * atPerThread;
int N = min(id0 + atPerThread, md->nAt);
int iat;
for (iat = id0; iat < N; iat++)
if (md->oldTypes[iat] != -1)
{
i = 0;
n = md->nangles[iat];
while (n && (i < md->nAngle))
{
if (md->angles[i].w)
if (md->angles[i].x == iat)
{
n--;
md->angles[i].w = 0;
}
i++;
}
// create new angles
t = md->specAngles[md->types[iat]]; // get type of angle, which formed by current atom type
if (t && (md->nbonds[iat] > 1)) // atom type supports angle creating and number of bonds is enough
{
// search of neighbors by bonds
i = 0; cnt = 0;
n = md->nbonds[iat];
while (n && (i < md->nBond))
{
if (md->bonds[i].z) // if bond isn't deleted
{
if (md->bonds[i].x == iat)
{
nei[cnt] = md->bonds[i].y;
cnt++;
n--;
}
else if (md->bonds[i].y == iat)
{
nei[cnt] = md->bonds[i].x;
cnt++;
n--;
}
}
i++;
}
// add new angles based on found neighbors:
for (i = 0; i < cnt-1; i++)
for (j = i + 1; j < cnt; j++)
{
ang = atomicAdd(&(md->nAngle), 1);
md->angles[ang].x = iat;
md->angles[ang].y = nei[i];
md->angles[ang].z = nei[j];
md->angles[ang].w = t;
}
n = (cnt * (cnt - 1)) / 2;
}
md->nangles[iat] = n;
// reset flag
md->oldTypes[iat] = -1;
}
}
specAngles配列には、指定された原子タイプに対応する結合角識別子が含まれています。次の関数は、すべての角度のエネルギーと力の計算を呼び出します。
__global__ void apply_angles(int iStep, int angPerBlock, int angPerThread, cudaMD* md)
// apply valence angle potentials
{
cudaAngle* ang;
// energies of angle potential
float eng;
__shared__ float shEng;
if (threadIdx.x == 0)
shEng = 0.0f;
__syncthreads();
int id0 = blockIdx.x * angPerBlock + threadIdx.x * angPerThread;
int N = min(id0 + angPerThread, md->nAngle);
int i;
for (i = id0; i < N; i++)
if (md->angles[i].w)
{
ang = &(md->angleTypes[md->angles[i].w]);
ang->force_eng(&(md->angles[i]), ang, md, eng);
}
// split energy to shared and then to global memory
atomicAdd(&shEng, eng);
__syncthreads();
if (threadIdx.x == 0)
atomicAdd(&(md->engAngl), shEng);
}さて、例えば、そのような角度の可能性は、調和コサイン関数を与えます。これは、フィールドforce_eng構造cudaAngleを示す可能性があります。
__device__ void angle_hcos(int4* angle, cudaAngle* type, cudaMD* md, float& eng)
// harmonic cosine valent angle potential:
// U = k / 2 * (cos(th)-cos(th0))^
{
float k = type->p0;
float cos0 = type->p1;
// indexes of central atom and ligands:
int c = angle->x;
int l1 = angle->y;
int l2 = angle->z;
// vector ij
float xij = md->xyz[l1].x - md->xyz[c].x;
float yij = md->xyz[l1].y - md->xyz[c].y;
float zij = md->xyz[l1].z - md->xyz[c].z;
delta_periodic(xij, yij, zij, md);
float r2ij = xij * xij + yij * yij + zij * zij;
float rij = sqrt(r2ij);
// vector ik
float xik = md->xyz[l2].x - md->xyz[c].x;
float yik = md->xyz[l2].y - md->xyz[c].y;
float zik = md->xyz[l2].z - md->xyz[c].z;
delta_periodic(xik, yik, zik, md);
float r2ik = xik * xik + yik * yik + zik * zik;
float rik = sqrt(r2ik);
float cos_th = (xij * xik + yij * yik + zij * zik) / rij / rik;
float dCos = cos_th - cos0; // delta cosinus
float c1 = -k * dCos;
float c2 = 1.0 / rij / rik;
atomicAdd(&(md->frs[c].x), -c1 * (xik * c2 + xij * c2 - cos_th * (xij / r2ij + xik / r2ik)));
atomicAdd(&(md->frs[c].y), -c1 * (yik * c2 + yij * c2 - cos_th * (yij / r2ij + yik / r2ik)));
atomicAdd(&(md->frs[c].z), -c1 * (zik * c2 + zij * c2 - cos_th * (zij / r2ij + zik / r2ik)));
atomicAdd(&(md->frs[l1].x), c1 * (xik * c2 - cos_th * xij / r2ij));
atomicAdd(&(md->frs[l1].y), c1 * (yik * c2 - cos_th * yij / r2ij));
atomicAdd(&(md->frs[l1].z), c1 * (zik * c2 - cos_th * zij / r2ij));
atomicAdd(&(md->frs[l2].x), c1 * (xij * c2 - cos_th * xik / r2ik));
atomicAdd(&(md->frs[l2].y), c1 * (yij * c2 - cos_th * yik / r2ik));
atomicAdd(&(md->frs[l2].z), c1 * (zij * c2 - cos_th * zik / r2ik));
eng += 0.5 * k * dCos * dCos;
}
「無効化された」コーナーを削除する機能は提供しません。基本的にclear_bondsと違いはありません。
の例
正確なふりをせずに、個々のイオンからの水分子の集合を描写しようとしました。対の電位は、バッキンガムの電位の形で任意に設定され、次に、水中のHO結合の長さに等しい平衡距離0.96Åで、調和電位の形で結合を作成する機能が追加されました。さらに、2番目のプロトンが酸素と結合したときに、酸素の頂点との結合角が追加されました。 100,000ステップ後、最初の分子はランダムに散乱したイオンから現れました。この図は、初期（左）と最終（右）の構成を示しています。
このような実験を設定できます。最初はすべての原子が同じであるとしますが、それらが互いに隣接している場合、それらは結合を形成します。結合した原子が、自由な原子または別の同様の結合した分子と別の結合を形成するようにします。その結果、原子が鎖状に並ぶ、ある種の自己組織化が得られます。
最終コメント
- ここでは、バインディングの基準を1つだけ使用しました。原則として、システムのエネルギーなど、他の基準もあります。実際には、化学結合が形成されると、原則として、エネルギーは熱の形で放出されます。これはここでは考慮されていませんが、たとえば粒子速度を変更するなど、実装を試みることができます。
- 化学結合ポテンシャルを介した粒子間の相互作用は、粒子が分子間ペア電位およびクーロン相互作用を介して相互作用できるという事実を否定するものではありません。もちろん、結合した原子の分子間相互作用を計算しないことは可能ですが、これは一般的な場合、長いチェックを必要とします。したがって、他の電位との合計が目的の機能を与えるように、化学結合電位を設定する方が簡単です。
- 粒子結合を並列に実装すると、速度が向上するだけでなく、粒子が互いに競合するため、よりリアルに見えます。
さて、これに非常に近いハブレに関するいくつかのプロジェクトがあります：