周期的境界条件による接触/配位数の計算

Question

各原子の半径、各原子の位置、およびボックスの大きさを含むN原子のデータがあります。私は各原子が接触している原子の数を計算し、その結果を保存したいと思います。これは、指定されたカットオフ距離内の原子数を計算するのと同じです。

このボックスは周期的な境界条件を持っているため、各原子間で計算された距離が最小になるように原子を再マッピングすることが困難です。ここ

私は今まで持っているものである。

import numpy as np 
Lx = 1.8609087768899999 
Ly = 1.8609087768899999 
Lz = 1.8609087768899999 
natoms = 8 
#Extract all atom positions 
atomID = [1, 2, 3, 4, 5, 6, 7, 8] 
x = [0.305153, 0.324049, 0.14708, 0.871106, 1.35602, 1.09239, 1.84683, 1.23874] 
y = [1.59871, 1.40195, 0.637672, 0.859931, 1.68457, 0.651108, 0.694614, 1.57241] 
z = [1.11666, 0.0698172, 1.60292, 0.787037, 0.381979, -0.00528695, 0.392633, 1.37821] 
radius = np.array([0.5 for i in range(natoms)]) 

#Create matrix of particle distance 
dx = np.subtract.outer(x,x) 
dy = np.subtract.outer(y,y) 
dz = np.subtract.outer(z,z) 
radius_sum = np.add.outer(radius,radius) 

#Account for periodic boundary conditions. 
dx = dx - Lx*np.around(dx/Lx) 
dy = dy - Ly*np.around(dy/Ly) 
dz = dz - Lz*np.around(dz/Lz) 
distance = np.array(np.sqrt(np.square(dx) + np.square(dy) + np.square(dz))) 

touching = np.where(distance < radius_sum, 1.0, 0.0) #elementwise conditional operator (condition ? 1 : 0 in this case) 
coordination_number = np.sum(touching, axis=0) - 1 # -1 to account for all diagonal terms being 1 in touching. 

X、Y、Zおよび半径は長さN（原子数）のアレイです。

ここで計算される調整値のいくつかは、シミュレーションを実行するために使用されるソフトウェアであるLAMMPS出力とは異なります。

何か助けていただければ幸いです。ここで予期しない動作がありますか？何か簡単なものがありませんか？

私は見つけるこの例では

：

1: [2, 4, 8] 
2: [1, 3, 5, 7] 
3: [2, 6, 7] 
4: [1, 6, 7, 8] 
5: [2, 6, 7, 8] 
6: [3, 4, 5, 7] 
7: [2, 3, 4, 5, 6] 
8: [1, 4, 5] 

これは、コードを用いて計算した：私は、さらに私のプログラムが接触しており、それらがあることが判明思っ原子を計算した

calculated expected 
3.0   4.0 
4.0   4.0 
3.0   4.0 
4.0   4.0 
4.0   4.0 
4.0   5.0 
5.0   5.0 
3.0   4.0 

：

contact_ID = [[] for i in range(natoms)] 
for i in range(natoms): 
    for j in range(natoms): 
     if i==j: continue 
     if (abs(touching[i,j] - radius_sum[i,j])) < .00001: 
      contact_ID[atomID[i]-1].append(atomID[j]) 
     else: continue 

私はシミュレーションソフトウェアから同じことをする方法はありません。私は原子3を見ることができます私のプログラムは原子3が原子2,6と7と接触していると言います。しかし4つの接続があると思います。原子（周期的な境界を有する）は、同じ原子と2回以上接触している可能性がある。これを今すぐコーディングするためのエレガントな方法を探しています。私はこれがわずか8原子でこの問題を解決すると思う。しかし、完全なシミュレーションでは、2000個のパーティクルとはるかに大きなボックスを使用していますが、まだ矛盾があります。

編集：以下のプルーンに示されているように、これは事実であり、原子3-6および1-8は実際にはダブルコンタクトであった。すばらしいです！しかし、このケースはすでに非常に非物理的であり、私のシミュレーションではエラーが発生することはありません。

この問題が発生しないような大きなサイズのボックスでも、問題は残ります。問題を再現できる最も簡単なケースは、L = 2.1、natoms = 21です。

ID calculated expected 
1 12.0 12.0 
8 11.0 11.0 
11 13.0 13.0 
14 10.0 10.0 
15 11.0 11.0 
2 11.0 12.0 
4 10.0 10.0 
5 12.0 12.0 
6 11.0 11.0 
7 12.0 12.0 
10 10.0 10.0 
3 12.0 12.0 
12 11.0 11.0 
13 11.0 11.0 
20 12.0 12.0 
21 10.0 10.0 
9 9.0 9.0 
17 11.0 11.0 
16 10.0 10.0 
18 10.0 10.0 
19 11.0 12.0 

ここでは、私は12が両方とも短いものを落下（そして唯一の不正な値であること持つためにそれらを計算している必要があり、一方、atomID 2及び19の両方を持つ原子は、11の接点を持っていることを参照してください。ここで私はの配位数を見つけますリスト中の）は、原子2と19が接触していることを意味する。手動で境界条件を考慮し、これらの原子間の最短距離を完了：

import numpy as np 
L=2.1 
atomID = [2, 19] 
x = [2.00551, 1.64908] 
y = [0.146456, 1.90822] 
z = [1.28094, 0.0518928] 

#Manually checking these values we find: 
dx = x[1] - x[0] #(= -0.35643) 
dy = y[1] - y[0] #(= 1.761764) 
dz = z[1] - z[0] #(= -1.2290472) 

#Account for period BC's: 
dx = dx   #(= -0.35643) 
dy = dy - L  #(= -0.338236) 
dz = dz + L  #(= 0.8709528) 

distance = np.sqrt(dx**2 + dy**2 + dz**2) 
print distance #(1.000002358) 

これは私のコードでは、実際にはないエラーがある意味...、シミュレーション・ソフトウェアは、これらの粒子を数えるということが理由ファンタスティックです実際に彼らはお互いから0.000002358離れているときに触れる？補正係数を追加せずに値を修正するにはどうしたらよいですか？

radius_sumに修正「スキン」ファクタを追加して完全なシミュレーションを再実行すると、自分の値が実際の結果を下回ることがよくありますが、場合によってはオーバーシュートします。それでも基本的な違いがあることを暗示していますか？

Answer 1

はい、ダブルタッチの問題では、8アトムシミュレーションの問題が説明されています。ダブルタッチのペアは、1-8と3-6で、どちらもxディメンションで囲みます。一方の次元の差がLx/2に近く、他方の2つの距離が比較的小さい場合に注意してください。これは、原子があなたのスペースラップの周りに完全に到達し、反対側に触れることを意味します。現在のコードでは、タッチの1つのみがカウントされます。

これを修正するには、Lx-radius-radius（LyとLzの繰り返し）の範囲の距離値をタッチするペアを見てください。あなたは "反転"する必要があると思います：距離= Lx距離。次に、あなたが変更したタッチするペアを含む別のセットの接触を探します。

これは、ボックスの寸法が半径の2倍よりも小さい場合にのみ発生する可能性があることに注意してください。そのような場合にエラー出力を提供できますか？大きなボックスで問題が発生しているので、上記の問題は大きなセットの問題を解決しません。おそらく8つの原子を持つが、2.1のボックスサイズ？