MDAnalysis / mdanalysis

MDAnalysis is a Python library to analyze molecular dynamics simulations.
https://mdanalysis.org
Other
1.32k stars 652 forks source link

Augment_coordinates ghost atoms with incorrect Z axis positions #4205

Open m-bone opened 1 year ago

m-bone commented 1 year ago

Expected behavior

Augment_coordinates should create ghost atoms within a specified cutoff from atoms existing within a universe. Used to mimic the effects of periodic boundaries when doing other calculations like counting neighbours.

Actual behavior

The majority of ghost atoms get added at the expected positions . However, a small selection of ghost atoms get placed within the system with the wrong Z coordinates, always far more negative than expected. For example, if the box is 5 Ang high then these atoms get placed with Z coordinates of -4, instead of 6. This can be seen in the attached files and the image below for methane, where two of the methane molecules are placed wrong. Using Ovito to overlay the LAMMPS data with the augmented coordinates in the xyz file shows the issue clearly.

image

The smaller atoms are the ghosts. The atoms in the lower part of the image should be at the top so that they've close enough to interact with any other atoms that might be in the original unit cell (in blue).

Code to reproduce the behavior

This Python code will recreate the ghost atoms xyz file and leave it in the specified directory.

import os
import numpy as np
import MDAnalysis as mda
from MDAnalysis.lib.distances import augment_coordinates

def augment_coord_resid(universe, searchRadius):   
    augmentList = []
    for resid in range(universe.atoms.resids.min(), universe.atoms.resids.max() +1): # +1 need as resids start at 1 not 0
        atoms = universe.select_atoms('resid ' + str(resid))

        ghostAtomCoords, _ = augment_coordinates(atoms.positions, universe.dimensions, searchRadius)
        augmentList.append(ghostAtomCoords)

    augmentArray = np.concatenate(augmentList)

    return augmentArray

# Consts
DIR = ''
FILENAME = 'methane.data'

os.chdir(DIR)
u = mda.Universe(FILENAME)
atoms = u.select_atoms('all')
ghostAtomCoords = augment_coord_resid(atoms, 5)

filename = "ghostAtoms.xyz"
with open(filename, "w") as f:
    f.write(f"{ghostAtomCoords.shape[0]}\n")  # Write the number of atoms/coordinates
    f.write("Generated grid\n")  # Write a comment line
    for i in range(ghostAtomCoords.shape[0]):
        f.write(f"X{i} {ghostAtomCoords[i, 0]} {ghostAtomCoords[i, 1]} {ghostAtomCoords[i, 2]}\n")  # Write each coordinate

Current version of MDAnalysis

Tested on MDA version 2.5.0 with Python 3.10.4 on Windows 10. Same issue seen on Ubuntu 20.04

Test Files (LAMMPS and current ghost atom xyz data) methane.zip

orbeckst commented 1 year ago

@hmacdope @richardjgowers can you have a look at this issue with the augment functionality?

hmacdope commented 1 year ago

Can do.

hmacdope commented 1 year ago

Sorry for the delay @m-bone can you clarify what you mean by original unit cell, vs the white cell?

m-bone commented 1 year ago

No problem, the original unit cell is the one with blue boundaries and contains the larger atoms. The white cell is the larger surrounding area with the white containing the ghost atoms represented with the smaller white balls.