Skip to content

Initial Attempt at PR for issue #3603 #5013

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 9 commits into
base: develop
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion package/AUTHORS
Original file line number Diff line number Diff line change
@@ -254,6 +254,7 @@ Chronological list of authors
- Yu-Yuan (Stuart) Yang
- James Rowe
- Debasish Mohanty
- Echo Fia


External code
@@ -302,4 +303,4 @@ Logo

The MDAnalysis 'Atom' logo was designed by Christian Beckstein; it is
Copyright (c) 2011 Christian Beckstein and made available under a
Creative Commons Attribution-NoDerivs 3.0 Unported License.
Creative Commons Attribution-NoDerivs 3.0 Unported License.
3 changes: 2 additions & 1 deletion package/CHANGELOG
Original file line number Diff line number Diff line change
@@ -59,7 +59,8 @@ Enhancements
* Enable parallelization for analysis.nucleicacids.NucPairDist (Issue #4670)
* Add check and warning for empty (all zero) coordinates in RDKit converter
(PR #4824)
* Added `precision` for XYZWriter (Issue #4775, PR #4771)
* Added `precision` for XYZWriter (Issue #4775, PR #4771)
* Improved tests for lineardensity.py (Issue #3603, PR #5013)

Changes
* MDAnalysis.analysis.psa, MDAnalysis.analysis.waterdynamics and
387 changes: 278 additions & 109 deletions testsuite/MDAnalysisTests/analysis/test_lineardensity.py
Original file line number Diff line number Diff line change
@@ -31,6 +31,8 @@
from MDAnalysis.core._get_readers import get_reader_for
from MDAnalysisTests.util import no_deprecated_call

from MDAnalysis.units import constants


def test_invalid_grouping():
"""Invalid groupings raise AttributeError"""
@@ -43,131 +45,298 @@ def test_invalid_grouping():
ld.run()


# test data for grouping='atoms'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't delete existing tests.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi! So, my understanding of the problem was that the test for lineardensity.py was using itself to check itself (i.e. the saved 'answers' to the test were generated using the function) and that the test cases didn't cover the idea of charged residues (or larger) as water is neutral.

As such, I used the code given that creates a universe of 100 atoms, modified it such that the charges + positions matched the 3 different types of systems (neutral particles, charged particles & charged dimers), and then calculated the values of (total masses, total charges, mass densities and charge densities) for atoms, residues and segments. I then edited the pytest function to compare these to the values generated by LinearDensity.

Let me know if you want me to explain further, I feel that the code itself is simple enough to read through (though maybe that's because I wrote it, so I might be tricking myself), so I assume the uncertainty comes in what I'm trying to achieve, which, I've hopefully explained. I had a lot of questions at the start which I put in issue #3603, but those haven't been seen yet so I interpreted it as best I could.

Apologies for removing the tests with the water system, I can add it back, updating it for how I've restructured the test (essentially, just adding in the universe for the test, which won't change anything for the water system). I can also go through and change the spellings to American English.

In regards to the comments about naming things 'make_universe' as opposed to 'make_Universe', should I change that or leave it? From looking at other code, it appeared that the standard thing done in the module was capitalised after the first word with underscores.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, whilst I have your attention, can I clarify if my GSoC proposal is still valid or not?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, what do you mean by "valid"?

(I also commented on your question #5015 )

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was under the impression that for a valid GSoC application to MDAnalysis, you needed:

  • Pre-proposal accepted
  • Proposal on GSoC website
  • PR merged

and I had assumed that the PR had to be merged by the 8th as well, with the proposal, so I was worried that I'd missed that deadline. Thanks for answering question #5015, that has clarified things, and hopefully I'll be able to get this PR merged before you make assessments.

expected_masses_atoms = np.array(
[
15.9994,
1.008,
1.008,
15.9994,
1.008,
1.008,
15.9994,
1.008,
1.008,
15.9994,
1.008,
1.008,
15.9994,
1.008,
1.008,
]
)
expected_charges_atoms = np.array(
[
-0.834,
0.417,
0.417,
-0.834,
0.417,
0.417,
-0.834,
0.417,
0.417,
-0.834,
0.417,
0.417,
-0.834,
0.417,
0.417,
]
)
expected_xmass_atoms = np.array(
[0.0, 0.0, 0.0, 0.00723323, 0.00473288, 0.0, 0.0, 0.0, 0.0, 0.0]
)
expected_xcharge_atoms = np.array(
[0.0, 0.0, 0.0, 2.21582311e-05, -2.21582311e-05, 0.0, 0.0, 0.0, 0.0, 0.0]
)
### Creating the Test Universes ###

# test data for grouping='residues'
expected_masses_residues = np.array(
[18.0154, 18.0154, 18.0154, 18.0154, 18.0154]
)
expected_charges_residues = np.array([0, 0, 0, 0, 0])
expected_xmass_residues = np.array(
[0.0, 0.0, 0.0, 0.00717967, 0.00478644, 0.0, 0.0, 0.0, 0.0, 0.0]
)
expected_xcharge_residues = np.array(
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
)

# test data for grouping='segments'
expected_masses_segments = np.array([90.0770])
expected_charges_segments = np.array([0])
expected_xmass_segments = np.array(
[0.0, 0.0, 0.0, 0.01196611, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
)
expected_xcharge_segments = np.array(
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
)
def make_Universe(
coords, charges, masses, n_atoms, n_frames, atomsPerRes, resPerSegs
):
"""Generate a reference universe of 100 atoms with uniformly drawn random positions."""
n_residues = n_atoms // atomsPerRes # Arbitrarily 5 atoms per residue
n_segments = n_residues // resPerSegs # Arbitrarily 4 residues per segment

# test data for grouping='fragments'
expected_masses_fragments = np.array(
[18.0154, 18.0154, 18.0154, 18.0154, 18.0154]
)
expected_charges_fragments = np.array([0, 0, 0, 0, 0])
expected_xmass_fragments = np.array(
[0.0, 0.0, 0.0, 0.00717967, 0.00478644, 0.0, 0.0, 0.0, 0.0, 0.0]
)
expected_xcharge_fragments = np.array(
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
)
# Indexing atoms into residues & residues into segments
atom_resindex = np.repeat(np.arange(n_residues), atomsPerRes)
residue_segindex = np.repeat(np.arange(n_segments), resPerSegs)

# Creating the universe
u = mda.Universe.empty(
n_atoms=n_atoms,
n_residues=n_residues,
n_segments=n_segments,
atom_resindex=atom_resindex,
residue_segindex=residue_segindex,
)

# Assigning the Charges & Masses
u.add_TopologyAttr("charges", values=charges)
u.add_TopologyAttr("masses", values=masses)

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


### Defining the Systems ###

# The masses are all set to 1
# The charges are dependent on each system
# The positions are randomly taken from a uniform distribution between 0 and 1
# NOTE: The positions at each time frame are independent on the previous time step


def neutral_Particles(n_atoms, n_frames, atomsPerRes, resPerSegs):
charges = np.zeros(n_atoms)
masses = np.ones(n_atoms)

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

return make_Universe(
coords, charges, masses, n_atoms, n_frames, atomsPerRes, resPerSegs
)


def charged_Particles(n_atoms, n_frames, atomsPerRes, resPerSegs):
charges = np.random.random(n_atoms) * 2 - np.ones(
n_atoms
) # Between -1 and 1
masses = np.ones(n_atoms)

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

return make_Universe(
coords, charges, masses, n_atoms, n_frames, atomsPerRes, resPerSegs
)


def charged_Dimers(
n_dimers=100, n_frames=1, dimersPerRes=5, resPerSegs=4, dimerLength=0.05
):
n_atoms = 2 * n_dimers

charges = np.random.random(n_atoms) * 2 - np.ones(
n_atoms
) # Between -1 and 1
masses = np.ones(n_atoms)

shape = (n_frames, n_dimers, 3)
coords = np.random.random(shape) * 0.9 + np.ones(shape) * 0.05
# Puts in the same coordinate twice per dimer
coords = np.repeat(coords, 2, axis=1)

# Shifts one of the atoms of each dimer by their bondLength
# in a random direction (defined to be in the box)
for time in coords:
for coord in time[::2, :]:
phi = np.random.random() * 2 * np.pi
theta = np.random.random() * np.pi
x = (
np.array(
[
np.sin(theta) * np.cos(phi),
np.sin(theta) * np.sin(phi),
np.cos(theta),
]
)
* dimerLength
)
coord += (
np.array(
[
np.sin(theta) * np.cos(phi),
np.sin(theta) * np.sin(phi),
np.cos(theta),
]
)
* dimerLength
)

return make_Universe(
coords,
charges,
masses,
n_atoms,
n_frames,
dimersPerRes * 2,
resPerSegs,
)


### Calculating the Expected Values ###


def calc_Prop(u, prop="masses"): # Property can be 'masses' or 'charges'
expected_atoms = getattr(u.atoms, prop)
expected_residues = np.array(
[sum(getattr(res.atoms, prop)) for res in u.residues]
)
expected_segments = np.array(
[sum(getattr(seg.atoms, prop)) for seg in u.segments]
)

return expected_atoms, expected_residues, expected_segments


def find_Centres(
groups,
): # Always based on Centre of Mass (NOT Centre of Charge)
centres = []
for group in groups:
total_Mass = sum(group.atoms.masses)
if total_Mass != 0:
centres.append(
np.sum(
group.atoms.positions.transpose()
* abs(group.atoms.masses),
axis=1,
)
/ total_Mass
)
# Just in case there are particles with 0 mass (somehow)
elif total_Mass == 0:
centres.append(
np.sum(
group.atoms.positions.transpose()
* abs(group.atoms.masses),
axis=1,
)
/ len(group.atoms)
)

return np.array(centres)


def calc_Densities(
u, prop="masses", spliceLen=0.25
): # Property can be 'masses' or 'charges'

propShort = "mass"
if prop == "charges":
propShort = "charge"

### Atoms
expected_atoms = np.zeros((3, int(u.dimensions[0] // spliceLen))).astype(
float
) # Works for cubic Universe
for atom in u.atoms:
for i in range(3):
expected_atoms[i][int(atom.position[i] // spliceLen)] += getattr(
atom, propShort
)

_, residue_Totals, segment_Totals = calc_Prop(
u, prop
) # Total of Charge OR Mass
### Residues
expected_residues = np.zeros(
(3, int(u.dimensions[0] // spliceLen))
).astype(float)
residue_Centres = find_Centres(u.residues)

for i in range(len(residue_Centres)):
for j in range(3):
expected_residues[j][
int(residue_Centres[i][j] // spliceLen)
] += residue_Totals[i]

### Segments
expected_segments = np.zeros(
(3, int(u.dimensions[0] // spliceLen))
).astype(float)
segment_Centres = find_Centres(u.segments)

for i in range(len(segment_Centres)):
for j in range(3):
expected_segments[j][
int(segment_Centres[i][j] // spliceLen)
] += segment_Totals[i]

# Scaling based on splice volumes & converting units
for i in range(3):
expected_atoms[i, :] /= (
spliceLen * u.dimensions[(i + 1) % 3] * u.dimensions[(i + 2) % 3]
)
expected_residues[i, :] /= (
spliceLen * u.dimensions[(i + 1) % 3] * u.dimensions[(i + 2) % 3]
)
expected_segments[i, :] /= (
spliceLen * u.dimensions[(i + 1) % 3] * u.dimensions[(i + 2) % 3]
)
expected_atoms /= (
constants["N_Avogadro"] * 1e-24
) # To be consistent with lineardensity.py
expected_residues /= (
constants["N_Avogadro"] * 1e-24
) # To be consistent with lineardensity.py
expected_segments /= (
constants["N_Avogadro"] * 1e-24
) # To be consistent with lineardensity.py

return expected_atoms, expected_residues, expected_segments


### Creating the Expected Values ###

spliceLen = 0.25

universes = []
groupings = ["atoms", "residues", "segments"]
# Will be [[Atoms, Residues, Segments (of Universe 1)], [... of Universe 2], ...]
expected_masses = []
expected_charges = []
expected_xmass = []
expected_xcharge = []

test_Universes = [neutral_Particles, charged_Particles, charged_Dimers]
test_Params = [[100, 1, 5, 4], [100, 1, 5, 4], [100, 1, 5, 4, 0.05]]

for i in range(len(test_Universes)):

universe = test_Universes[i](*test_Params[i])
universes.append(universe)

expected_masses.append(calc_Prop(universe, "masses"))
expected_charges.append(calc_Prop(universe, "charges"))
expected_xmass.append(calc_Densities(universe, "masses", spliceLen))
expected_xcharge.append(calc_Densities(universe, "charges", spliceLen))

### Parametrizing for Pytest ###

# NOTE: There is an extra [0] after the expected_xmass and expected_xcharge as the original data has all 3 co-ords, but only comparing to x here
pytest_Params = [
(
universes[i],
groupings[j],
expected_masses[i][j],
expected_charges[i][j],
expected_xmass[i][j][0],
expected_xcharge[i][j][0],
)
for j in range(len(groupings))
for i in range(len(universes))
]


@pytest.mark.parametrize(
"grouping, expected_masses, expected_charges, expected_xmass, expected_xcharge",
[
(
"atoms",
expected_masses_atoms,
expected_charges_atoms,
expected_xmass_atoms,
expected_xcharge_atoms,
),
(
"residues",
expected_masses_residues,
expected_charges_residues,
expected_xmass_residues,
expected_xcharge_residues,
),
(
"segments",
expected_masses_segments,
expected_charges_segments,
expected_xmass_segments,
expected_xcharge_segments,
),
(
"fragments",
expected_masses_fragments,
expected_charges_fragments,
expected_xmass_fragments,
expected_xcharge_fragments,
),
],
"universe, grouping, expected_masses, expected_charges, expected_xmass, expected_xcharge",
pytest_Params,
)
def test_lineardensity(
universe,
grouping,
expected_masses,
expected_charges,
expected_xmass,
expected_xcharge,
):
universe = mda.Universe(waterPSF, waterDCD)
sel_string = "all"
selection = universe.select_atoms(sel_string)
ld = LinearDensity(selection, grouping, binsize=5).run()
ld = LinearDensity(selection, grouping, binsize=0.25).run()
assert_allclose(ld.masses, expected_masses)
assert_allclose(ld.charges, expected_charges)
# rtol changed here due to floating point imprecision