mdanalysis icon indicating copy to clipboard operation
mdanalysis copied to clipboard

LinearDensity should have more reliable tests

Open BFedder opened this issue 2 years ago • 8 comments

Is your feature request related to a problem?

Currently, some of the tests in test_lineardensity.py are basically tautologies because the expected values the tests check against were obtained by running lineardensity.py. This is the case for the following expected values:

expected_xchar_atoms expected_xpos_residues expected_xchar_residues expected_xpos_segments expected_xchar_segments expected_xpos_fragments expected_xchar_fragments

The current test Universe u = mda.Universe(waterPSF, waterDCD) is not ideal for determining charge densities by residues, segments, or fragments, because in these groupings the overall charge is zero.

Describe the solution you'd like

The expected values should be obtained independently from LinearDensity so that the tests are actually meaningful.

@PicoCentauri had the following suggestion for how this might be achieved:

Generate a trajectory based on random positions of particles in a cubic box. Draw the positions from a uniform distribution. Run the analysis and compare the densities. You could do three kinds of system: uncharged particles, charged particles and maybe charged dimers.

Describe alternatives you've considered

Currently, expected_xpos_atoms came from a different source. One could look into how these values were originally determined and replicate this for the other test cases.

Additional context

This came up in a discussion on PR #3572

BFedder avatar Apr 05 '22 11:04 BFedder

Thanks for this issue. If somebody wants to start on this. The following lines can be used to generate a Universe of individual atoms.

import MDAnalysis as mda
import numpy as np
from MDAnalysis.core._get_readers import get_reader_for

def make_Universe():
    """Generate a reference universe of 100 atoms with uniformly drawn random positions."""
    n_atoms = 100
    n_frames = 100

    u = mda.Universe.empty(n_atoms=n_atoms,
                           n_residues=n_atoms,
                           n_segments=n_atoms,
                           atom_resindex=np.arange(n_atoms),
                           residue_segindex=np.arange(n_atoms))

    for attr in ["charges", "masses"]:
        u.add_TopologyAttr(attr, values=np.ones(n_atoms))

    shape = (n_frames, n_atoms, 3)
    coords = np.random.random(shape)

    u.trajectory = get_reader_for(coords)(coords,
                                          order='fac',
                                          n_atoms=n_atoms)

    for ts in u.trajectory:
        ts.dimensions = np.array([1, 1, 1, 90, 90, 90])

    return u

PicoCentauri avatar Apr 06 '22 16:04 PicoCentauri

i want to start on this

olivia632 avatar Apr 09 '22 02:04 olivia632

@olivia632 cool!

I suggest that you calculate the density using the make_universe function and compare it to the values you would calculate based on the system parameters. The number of particles is 100, the volume is 1 Å^3, the mass per particle is 1 u and the charge is 1 e. If I am not mistaken the mass density in each bin and direction should be 100 u/Å^3 and the charge 100 e/Å^3.

PicoCentauri avatar Apr 10 '22 06:04 PicoCentauri

k I will work on that and then reach you back

olivia632 avatar Apr 10 '22 09:04 olivia632

You can also reach me on discord if you have any general questions.

PicoCentauri avatar Apr 10 '22 09:04 PicoCentauri

there are 3 test cases failing in the original code can you clarify on that

olivia632 avatar Apr 15 '22 14:04 olivia632

==================================== short test summary info ==================================== FAILED analysis/test_lineardensity.py::test_invalid_grouping - TypeError: 'method' object is n... FAILED analysis/test_lineardensity.py::test_lineardensity[residues-expected_masses1-expected_charges1-expected_xpos1-expected_xchar1] FAILED analysis/test_lineardensity.py::test_lineardensity[segments-expected_masses2-expected_charges2-expected_xpos2-expected_xchar2] =========================== 3 failed, 2 passed, 15 warnings in 0.68s ============================

olivia632 avatar Apr 15 '22 14:04 olivia632

@olivia632 did you figures out what is going wrong (Pull the to latest version from github etc?). The tests seem to work for the CI...

PicoCentauri avatar Apr 21 '22 12:04 PicoCentauri