DensityAnalysis Normalization

[nawtrey]

The current density function may have an issue with normalization. Since DensityAnalysis uses the MDAnalysis density function, it is currently unclear whether MDAnalysis.density.density_from_universe() is incorrect or the implementation of density_from_universe() into pmda is incorrect. We need to compare the results from a run of MDAnalysis.analysis.density.density_from_universe() to pmda.density.DensityAnalysis for other trajectories to see if they agree, as well as find a 3rd party method of calculating densities to compare values with.

[orbeckst]

I doubt that packaging the histogram from PMDA into Density instance in

pmda/pmda/density.py

Lines 314 to 317 in 1e99541

	density = Density(grid=self._grid, edges=self._edges,
	units={'length': "Angstrom"},
	parameters=parameters,
	metadata=metadata)

is the problem – but do check!

As an alternative density calculator you can use VMD's VolMap.

[orbeckst]

Are you sure that

pmda/pmda/density.py

Line 303 in 1e99541

self._grid /= float(self._n_frames)

is correct, especially if one changes the run(start, stop, step) parameters?

There's some "we're trying to be smart"-code in Density https://github.com/MDAnalysis/mdanalysis/blob/d1b0e659329b1d98fc2422c88614a196976f3497/package/MDAnalysis/analysis/density.py#L331-L336 (which you can blame on me...) – can you check that Density is not automatically setting parameters['isDensity'] to True, which would mess up everything?

[orbeckst]

@nawtrey can you fix this one so that we can make a 0.3.0 release #106 ?

[orbeckst]

Can you please also look into why DensityAnalysis failed after merging into master? See https://travis-ci.org/MDAnalysis/pmda/jobs/592141990

=================================== FAILURES ===================================

_____________________________ test_n_frames[0-5-2] _____________________________

u = <Universe with 43480 atoms>, start = 0, stop = 5, step = 2

    @pytest.mark.parametrize("step", [1, 2])

    @pytest.mark.parametrize("stop", [5, 4])

    @pytest.mark.parametrize("start", [0, 1])

    def test_n_frames(u, start, stop, step):

        D = pmda.density.DensityAnalysis(u.atoms)

        D.run(start, stop, step)

>       assert D._n_frames == len(range(start, stop, step))

E       assert 5 == 3

E        +  where 5 = <pmda.density.DensityAnalysis object at 0x7f1d0ac42e48>._n_frames

E        +  and   3 = len(range(0, 5, 2))

E        +    where range(0, 5, 2) = range(0, 5, 2)

/home/travis/build/MDAnalysis/pmda/pmda/test/test_density.py:85: AssertionError

_____________________________ test_n_frames[0-4-1] _____________________________

u = <Universe with 43480 atoms>, start = 0, stop = 4, step = 1

    @pytest.mark.parametrize("step", [1, 2])

    @pytest.mark.parametrize("stop", [5, 4])

    @pytest.mark.parametrize("start", [0, 1])

    def test_n_frames(u, start, stop, step):

        D = pmda.density.DensityAnalysis(u.atoms)

        D.run(start, stop, step)

>       assert D._n_frames == len(range(start, stop, step))

E       assert 5 == 4

E        +  where 5 = <pmda.density.DensityAnalysis object at 0x7f1d4494a160>._n_frames

E        +  and   4 = len(range(0, 4))

E        +    where range(0, 4) = range(0, 4, 1)

/home/travis/build/MDAnalysis/pmda/pmda/test/test_density.py:85: AssertionError

_____________________________ test_n_frames[0-4-2] _____________________________

u = <Universe with 43480 atoms>, start = 0, stop = 4, step = 2

    @pytest.mark.parametrize("step", [1, 2])

    @pytest.mark.parametrize("stop", [5, 4])

    @pytest.mark.parametrize("start", [0, 1])

    def test_n_frames(u, start, stop, step):

        D = pmda.density.DensityAnalysis(u.atoms)

        D.run(start, stop, step)

>       assert D._n_frames == len(range(start, stop, step))

E       assert 5 == 2

E        +  where 5 = <pmda.density.DensityAnalysis object at 0x7f1d2f8eca90>._n_frames

E        +  and   2 = len(range(0, 4, 2))

E        +    where range(0, 4, 2) = range(0, 4, 2)

/home/travis/build/MDAnalysis/pmda/pmda/test/test_density.py:85: AssertionError

_____________________________ test_n_frames[1-5-1] _____________________________

u = <Universe with 43480 atoms>, start = 1, stop = 5, step = 1

    @pytest.mark.parametrize("step", [1, 2])

    @pytest.mark.parametrize("stop", [5, 4])

    @pytest.mark.parametrize("start", [0, 1])

    def test_n_frames(u, start, stop, step):

        D = pmda.density.DensityAnalysis(u.atoms)

        D.run(start, stop, step)

>       assert D._n_frames == len(range(start, stop, step))

E       assert 5 == 4

E        +  where 5 = <pmda.density.DensityAnalysis object at 0x7f1d1ba772b0>._n_frames

E        +  and   4 = len(range(1, 5))

E        +    where range(1, 5) = range(1, 5, 1)

/home/travis/build/MDAnalysis/pmda/pmda/test/test_density.py:85: AssertionError

_____________________________ test_n_frames[1-5-2] _____________________________

u = <Universe with 43480 atoms>, start = 1, stop = 5, step = 2

    @pytest.mark.parametrize("step", [1, 2])

    @pytest.mark.parametrize("stop", [5, 4])

    @pytest.mark.parametrize("start", [0, 1])

    def test_n_frames(u, start, stop, step):

        D = pmda.density.DensityAnalysis(u.atoms)

        D.run(start, stop, step)

>       assert D._n_frames == len(range(start, stop, step))

E       assert 5 == 2

E        +  where 5 = <pmda.density.DensityAnalysis object at 0x7f1d4f8e1940>._n_frames

E        +  and   2 = len(range(1, 5, 2))

E        +    where range(1, 5, 2) = range(1, 5, 2)

/home/travis/build/MDAnalysis/pmda/pmda/test/test_density.py:85: AssertionError

_____________________________ test_n_frames[1-4-1] _____________________________

u = <Universe with 43480 atoms>, start = 1, stop = 4, step = 1

    @pytest.mark.parametrize("step", [1, 2])

    @pytest.mark.parametrize("stop", [5, 4])

    @pytest.mark.parametrize("start", [0, 1])

    def test_n_frames(u, start, stop, step):

        D = pmda.density.DensityAnalysis(u.atoms)

        D.run(start, stop, step)

>       assert D._n_frames == len(range(start, stop, step))

E       assert 5 == 3

E        +  where 5 = <pmda.density.DensityAnalysis object at 0x7f1d090fd6a0>._n_frames

E        +  and   3 = len(range(1, 4))

E        +    where range(1, 4) = range(1, 4, 1)

/home/travis/build/MDAnalysis/pmda/pmda/test/test_density.py:85: AssertionError

_____________________________ test_n_frames[1-4-2] _____________________________

u = <Universe with 43480 atoms>, start = 1, stop = 4, step = 2

    @pytest.mark.parametrize("step", [1, 2])

    @pytest.mark.parametrize("stop", [5, 4])

    @pytest.mark.parametrize("start", [0, 1])

    def test_n_frames(u, start, stop, step):

        D = pmda.density.DensityAnalysis(u.atoms)

        D.run(start, stop, step)

>       assert D._n_frames == len(range(start, stop, step))

E       assert 5 == 2

E        +  where 5 = <pmda.density.DensityAnalysis object at 0x7f1d1bb6f9e8>._n_frames

E        +  and   2 = len(range(1, 4, 2))

E        +    where range(1, 4, 2) = range(1, 4, 2)

/home/travis/build/MDAnalysis/pmda/pmda/test/test_density.py:85: AssertionError

[nawtrey]

It failed because the number of frames is grabbed from the length of the trajectory in __init__(), not from the user-input start, stop and step. So I wrote a test that runs through different pieces of the trajectory and checked it against len(range(start, stop, step)) to verify that this was the problem with the normalization, which it is.

I'm currently working on a solution to this since the number of frames is never saved in the run method or dask helper. I looked at the RMSF function for inspiration since it has to divide by a similar value in _conclude() for its final RMSF value calculation, and as it turns out this 'totalts' value is not 100% accurate either. So I need to find a solution for both functions.

[orbeckst]

Ok, my fault – I looked at the wrong travis run. Thanks for the explanation.