Évaluer le MSD de ma simulation

Question

J'exécute une simulation de dynamique moléculaire de l'eau dans TIP3P, et je stocke les positions de mes particules dans un tableau 2D appelé relevant_positions. Le nombre de particules dans ma simulation est de numPart. J'exécute une simulation de t=0à t=n_time_points-1. En effet, j'ai des instantanés des positions des particules à des n_time_pointsmoments précis.

J'essaye d'évaluer le MSD de ma simulation, et voici le code que j'utilise:

for d in range(1, n_time_points):
    for i in range(0, n_time_points-d):
        msd[d] += np.sum(np.square(\
        relevant_positions[numPart*(d+i):numPart*(d+i+1),:] -\
        relevant_positions[numPart*i:numPart*(i+1),:]))
    msd[d] = msd[d]/(n_time_points-d) 

msd = msd/numPart

Le résultat que j'obtiens est:

Je m'attends à ce que ce soit une ligne droite, mais ce n'est clairement pas le cas. Qu'est-ce que je fais mal ici?

Answer 1

Je n'ai pas vraiment d'expérience préalable avec la dynamique moléculaire, mais pour commencer, j'ai trouvé une ressource avec une méthode assez détaillée de calcul du MSD (déplacement carré moyen).

En regardant votre intrigue, je soupçonne que les conditions aux limites périodiques vous causent en quelque sorte un problème. Il semble qu'il converge vers une valeur qui se produit lorsque tout se disperse aussi loin que possible en moyenne à partir de sa position de départ. Lorsqu'un atome s'enroule autour de la cellule, sa position passera d'un grand nombre (tel que 10,5 angströms) à un petit nombre (0,1 angströms). Cela sera considéré comme un déplacement par votre code, je crois.

Si vous voulez une très mauvaise solution pour corriger votre code, vous pouvez supprimer tout déplacement supérieur à une valeur spécifique correspondant à quelque chose comme la moitié de la cellule unitaire (une distance impossible à déplacer dans votre simulation MD). Cela supprimerait les points de données problématiques.

Utiliser ASE pour les calculs de distance entre les atomes pour chaque pas de temps vous permettrait de prendre en compte les conditions aux limites périodiques car il prend en charge la prise de distances par rapport à lui. Ce serait une bien meilleure solution.

Answer 2

D'accord, il s'avère donc que @TristanMaxson avait raison - ce sont les conditions aux limites périodiques qui perturbaient mon calcul. La solution était de trouver les coordonnées non emballées de mon système.

Voici comment j'ai déballé mes coordonnées:

unwrapped_positions = relevant_positions.copy()
for ts in range(0,n_time_points-1):
    periodic_displacement = \
    relevant_positions[numPart*(ts+1):numPart*(ts+2),:]-\
    relevant_positions[numPart*(ts):numPart*(ts+1),:]

    boundary_crossing = (np.abs(periodic_displacement) > (L/2))*1
    boundary_crossing_sign = np.sign(periodic_displacement)

    absolute_displacement = periodic_displacement\
    - boundary_crossing*boundary_crossing_sign*L

    unwrapped_positions[numPart*(ts+1):numPart*(ts+2),:] = \
    unwrapped_positions[numPart*(ts):numPart*(ts+1),:] + \
    absolute_displacement

suivie par:

msd = np.zeros(np.shape(t))
msd[0] = 0

for d in range(1, n_time_points):
    for i in range(0, n_time_points-d):
        msd[d] += np.sum(np.square(\
        unwrapped_positions[numPart*(d+i):numPart*(d+i+1),:] -\
        unwrapped_positions[numPart*i:numPart*(i+1),:]))
    msd[d] = msd[d]/(n_time_points-d)

msd = msd/numPart

Comme vous pouvez le voir, il y avait une faille inhérente à mon code - je ne faisais pas la moyenne de mon MSD, mis à part le fait que j'avais des problèmes avec les coordonnées encapsulées.