Add Einstein-Helfand viscosity implementation #25

xhgchen · 2023-07-14T23:42:12Z

Fixes #17 and #9

Changes made in this Pull Request:

Use class VelocityAutocorr as template
Add items to __init__()
Write _prepare() and _single_frame()
Early work on _conclude()

Areas to work on:

Complete _conclude()
Units
Boltzmann constant usage in Python

For later:

Add FFT calculation? Check tidynamics source code to see how feasible this is (raised in Consider adding FFT support to Einstein-Helfand viscosity #27)

codecov · 2023-07-14T23:46:18Z

Codecov Report

Merging #25 (b862476) into main (1211874) will not change coverage.
The diff coverage is 100.00%.

Additional details and impacted files

orionarcher · 2023-07-17T18:37:37Z

transport_analysis/viscosity.py

+
+        # args
+        self.dim_type = dim_type
+        self._parse_dim_type()


I think this should be a static method.

self.dim_type = dim_type self._parse_dim_type()

is much less clear than

self.dim_type, self.dim_fac = self._parse_dim_type(dim_type)

The latter makes it clearer what _parse_dim_type is doing. In the former, you don't even know that a self.dim_fac variable has been created.

orionarcher · 2023-07-17T18:37:43Z

transport_analysis/viscosity.py

+    def _parse_dim_type(self):
+        r"""Sets up the desired dimensionality of the calculation."""
+        keys = {
+            "x": [0],
+            "y": [1],
+            "z": [2],
+            "xy": [0, 1],
+            "xz": [0, 2],
+            "yz": [1, 2],
+            "xyz": [0, 1, 2],
+        }
+
+        self.dim_type = self.dim_type.lower()
+
+        try:
+            self._dim = keys[self.dim_type]
+        except KeyError:
+            raise ValueError(
+                "invalid dim_type: {} specified, please specify one of xyz, "
+                "xy, xz, yz, x, y, z".format(self.dim_type)
+            )
+
+        self.dim_fac = len(self._dim)


Could we make this a static method that returns dim_type and dim_fac?

Generally, I am of the opinion that it's always best to use static methods when possible. It makes it much easier to inspect the function call and understand what's happening.

So do you mean we would set the instance variables self._dim and self.dim_fac in __init__()?

Oh, I see what you mean now. Missed the comment above, my bad!

orionarcher · 2023-07-17T18:41:34Z

transport_analysis/viscosity.py

+        masses = self._mass_reshape.astype(np.float64)
+        velocities = self._velocity_array.astype(np.float64)
+        positions = self._position_array.astype(np.float64)


I believe that the default float in python is already 64 bits, what will this do?

orionarcher · 2023-07-17T18:51:28Z

transport_analysis/viscosity.py

+class ViscosityHelfand(AnalysisBase):
+    """


What's our plan for benchmarking the viscosity analysis?

@xhgchen you can benmark as part of your reproduction.

xhgchen · 2023-07-21T00:25:28Z

With the most recent commit, the calculation now actually has all the necessary elements! Still very early though and I definitely need to do some proper testing.

orionarcher · 2023-07-27T17:24:08Z

transport_analysis/viscosity.py

+            diff = (
+                self._masses_rs
+                * self._velocities[:-lag, :, :]
+                * self._positions[:-lag, :, :]
+                - self._masses_rs
+                * self._velocities[lag:, :, :]
+                * self._positions[lag:, :, :]
+            )


Can we break down this math to make it more understandable? something like

something = self._masses_rs * self._velocities[:-lag, :, :] * self._positions[:-lag, :, :] something_else = self._masses_rs * self._velocities[:-lag, :, :] * self._positions[:-lag, :, :] diff = something - something_else

orionarcher · 2023-07-27T17:25:23Z

transport_analysis/viscosity.py

+            self.results.visc_by_particle = self.results.visc_by_particle / (
+                2
+                * constants["Boltzman_constant"]
+                * self._vol_avg
+                * self.temp_avg
+            )


same as above here, I'd advocate for adding intermediate variables over putting a single term on each line. Ultimately a matter of taste though. I sense this will be a reoccuring thing so curious what @hmacdope thinks.

Wouldn't adding intermediate variables take more memory and make the calculation scale worse with larger datasets? If the difference doesn't really matter, I would be good with adding them. I'm also curious what @hmacdope thinks.

I'm actually not sure. I suspect the difference is negligable if it exists at all. This could be something to try out with %timeit in a jupyter notebook.

transport_analysis/viscosity.py

* Use class `VelocityAutocorr` as template * Add items to `__init__()` * Write `_prepare()` and `_single_frame()` * Early work on `_conclude()`

* Change `self.*_array` to simply `self.*` * Fix `self._masses_rs`

* Set up `self._volumes` in `_prepare()` and `_single_frame()`

* Complete `_conclude()` with 2, Boltzman constant, average volume, and average temp

* Check that AtomGroup is accepted * Test some exceptions * Functionality not yet tested here

* Move `visc_by_particle` division step outside of for loop

* Contains mass, velocities, positions, and volume

* Format `viscosity.py` with Black

* Simplify `characteristic_poly_helfand()`

* CI will not work otherwise

* use try except and update with mda 2.6.0 release

* Improve accuracy and agreement with `ViscosityHelfand` class * Improve performance

xhgchen · 2023-08-11T22:19:44Z

characteristic_poly_helfand() is now vectorized, so that is half of #26 taken care of.

hmacdope

Quick review, looking good!

hmacdope · 2023-08-12T04:23:58Z

transport_analysis/tests/test_viscosity.py

+
+
+@pytest.fixture(scope="module")
+def NSTEP():


No need to back this in as a fixtue, just use 5001 in place.

Sure, just wondering why 5001 in place is better? I thought setting NSTEP as a fixture would be more maintainable in case we ever decide to test it with another value (and your MSD module uses the fixture too)?

fair enough I agree.

hmacdope · 2023-08-12T04:24:36Z

transport_analysis/tests/test_viscosity.py

+
+@pytest.fixture(scope="module")
+def visc_helfand(ag):
+    vh_t = VH(ag, fft=False)


Parametrize over fft=True and fft=False

I'll remove FFT for now then and raise an issue.

hmacdope · 2023-08-12T04:25:24Z

transport_analysis/tests/test_viscosity.py

+    start=0,
+    step=1,
+):
+    # update when mda 2.6.0 releases with typo fix


Link to the issue number.

hmacdope · 2023-08-12T04:25:54Z

transport_analysis/tests/test_viscosity.py

+    step=1,
+):
+    # update when mda 2.6.0 releases with typo fix
+    try:


Grab the actual value here.

Was wondering also - why is grabbing the actual value better here?

Just to reduce code volume, you still have to index the dict later with this approach.

hmacdope · 2023-08-12T04:26:27Z

transport_analysis/tests/test_viscosity.py

+
+
+class TestViscosityHelfand:
+    # fixtures are helpful functions that set up a test


Can remove this comment now.

hmacdope · 2023-08-12T04:27:19Z

transport_analysis/viscosity.py

+class ViscosityHelfand(AnalysisBase):
+    """


@xhgchen you can benmark as part of your reproduction.

hmacdope · 2023-08-12T04:30:58Z

transport_analysis/viscosity.py

+        self.temp_avg = temp_avg
+        self.dim_type = dim_type.lower()
+        self._dim, self.dim_fac = self._parse_dim_type(self.dim_type)
+        # self.fft = fft # consider whether fft is possible later


If its not possibe raise exception or remove.

hmacdope · 2023-08-12T04:31:49Z

transport_analysis/viscosity.py

+            2 * constants[self.boltzmann] * self._vol_avg * self.temp_avg
+        )
+        # average over # particles and update results array
+        self.results.timeseries = self.results.visc_by_particle.mean(axis=1)


Add plotting?

I wanted to do that in a separate PR to avoid making this PR too bloated. Would it be better to add it here?

thats fine to do next PR

hmacdope

I am going to approve so we can move forward here. If any issues come up when testing, lets address them in a new PR. I would however love an additional look over by @orionarcher :)

xhgchen added the enhancement New feature or request label Jul 14, 2023

orionarcher reviewed Jul 17, 2023

View reviewed changes

xhgchen mentioned this pull request Jul 20, 2023

Add self-diffusivity calculation #24

Merged

orionarcher reviewed Jul 27, 2023

View reviewed changes

xhgchen added 12 commits August 8, 2023 11:48

Add Einstein-Helfand viscosity class skeleton

4537d45

* Use class `VelocityAutocorr` as template * Add items to `__init__()` * Write `_prepare()` and `_single_frame()` * Early work on `_conclude()`

Format viscosity.py with Black

f5d4371

Update docstrings in viscosity.py

c9187c4

Make variable names more readable

329e17b

* Change `self.*_array` to simply `self.*` * Fix `self._masses_rs`

Remove unnecessary type casting

0405ff3

Add initial volume handling

e00bd8a

* Set up `self._volumes` in `_prepare()` and `_single_frame()`

Add remaining items to viscosity calculation

56e9e14

* Complete `_conclude()` with 2, Boltzman constant, average volume, and average temp

Refactor _parse_dim_type() to static method

7416979

Write starter tests for Helfand viscosity class

512d57f

* Check that AtomGroup is accepted * Test some exceptions * Functionality not yet tested here

Fix visc_by_particle division step

98ad724

* Move `visc_by_particle` division step outside of for loop

Create full unit velocity trajectory for tests

577e131

* Contains mass, velocities, positions, and volume

Add characteristic polynomial for viscosity class

49216aa

* Format `viscosity.py` with Black

xhgchen force-pushed the viscosity-helfand branch from b2f9366 to 5000af0 Compare August 9, 2023 00:10

Add basic test for Helfand viscosity

b0e3033

* Simplify `characteristic_poly_helfand()`

xhgchen force-pushed the viscosity-helfand branch from 5000af0 to b0e3033 Compare August 9, 2023 00:15

xhgchen added 2 commits August 8, 2023 17:23

Use key "Boltzman_constant" until mda 2.6.0

06539e7

* CI will not work otherwise

Fix Boltzmann KeyError

760bf9d

* use try except and update with mda 2.6.0 release

xhgchen force-pushed the viscosity-helfand branch from a9857be to 760bf9d Compare August 9, 2023 03:25

xhgchen added 2 commits August 11, 2023 14:45

Vectorize basic Helfand viscosity test

248a5c1

* Improve accuracy and agreement with `ViscosityHelfand` class * Improve performance

Add tests for start, stop, step

4759a63

xhgchen added 2 commits August 11, 2023 15:25

Store Boltzmann constant key instead of value

016effe

Update docstrings for clarity

53e9cf3

xhgchen changed the title ~~Add Einstein-Helfand viscosity class skeleton~~ Add Einstein-Helfand viscosity implementation Aug 11, 2023

xhgchen marked this pull request as ready for review August 11, 2023 22:41

xhgchen requested review from orionarcher and hmacdope August 11, 2023 22:41

hmacdope requested changes Aug 12, 2023

View reviewed changes

xhgchen mentioned this pull request Aug 13, 2023

Consider adding FFT support to Einstein-Helfand viscosity #27

Open

xhgchen added 4 commits August 13, 2023 12:33

Remove leftover comments/unused code

17a3593

Get Boltzmann value instead of key

0e2a0c5

Use 5001 in place instead of fixture for tests

e482765

Comment Boltzmann issue from mdanalysis repo

9c7dee8

xhgchen requested a review from hmacdope August 13, 2023 20:10

hmacdope approved these changes Aug 13, 2023

View reviewed changes

Restore NSTEP and link issue number in tests

b862476

xhgchen merged commit 7c9b10d into MDAnalysis:main Aug 15, 2023
24 checks passed

xhgchen deleted the viscosity-helfand branch August 15, 2023 01:18

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add Einstein-Helfand viscosity implementation #25

Add Einstein-Helfand viscosity implementation #25

xhgchen commented Jul 14, 2023 •

edited

codecov bot commented Jul 14, 2023 •

edited

orionarcher Jul 17, 2023

orionarcher Jul 17, 2023

xhgchen Jul 18, 2023

xhgchen Jul 18, 2023

orionarcher Jul 17, 2023

orionarcher Jul 17, 2023

hmacdope Aug 12, 2023

xhgchen commented Jul 21, 2023

orionarcher Jul 27, 2023

orionarcher Jul 27, 2023

xhgchen Jul 27, 2023

orionarcher Aug 3, 2023

xhgchen commented Aug 11, 2023

hmacdope left a comment

hmacdope Aug 12, 2023

xhgchen Aug 13, 2023

hmacdope Aug 13, 2023

hmacdope Aug 12, 2023

xhgchen Aug 13, 2023

hmacdope Aug 13, 2023

hmacdope Aug 12, 2023

hmacdope Aug 12, 2023

xhgchen Aug 13, 2023

hmacdope Aug 13, 2023

hmacdope Aug 12, 2023

hmacdope Aug 12, 2023

hmacdope Aug 12, 2023

hmacdope Aug 12, 2023

xhgchen Aug 12, 2023

hmacdope Aug 12, 2023

hmacdope left a comment



		class TestViscosityHelfand:
		# fixtures are helpful functions that set up a test

Add Einstein-Helfand viscosity implementation #25

Add Einstein-Helfand viscosity implementation #25

Conversation

xhgchen commented Jul 14, 2023 • edited

codecov bot commented Jul 14, 2023 • edited

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

xhgchen commented Jul 21, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

xhgchen commented Aug 11, 2023

hmacdope left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hmacdope left a comment

Choose a reason for hiding this comment

xhgchen commented Jul 14, 2023 •

edited

codecov bot commented Jul 14, 2023 •

edited