Skip to content
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

Transport Properties using Onsager Formalism #56

Open
wants to merge 6 commits into
base: main
Choose a base branch
from

Conversation

rohithsrinivaas
Copy link

Fixes #44 #45 #18

Changes made in this Pull Request:

  • Added code to compute diffusion coefficient(Einstein Approach aka Slope fitting) based on the Onsager formalism with the FFT approach. The standard deviation of the Onsager coefficients are also computed with respect to time. If the linear fitting regime is not provided, weighted linear regression is performed using the inverse of the standard deviation as the weights. This solves the manual intervention for the identifying the suitable fitting region.
  • Added code to compute the conductivity, true transference number(Einstein Approach aka Slope fitting) based on the Onsager formalism with the FFT approach. The standard deviation of the Onsager coefficients are also computed with respect to time. If the linear fitting regime is not provided, weighted linear regression is performed using the inverse of the standard deviation as the weights. This solves the manual intervention for the identifying the suitable fitting region.
  • Added code to compute the conductivity, true transference number(Green Kubo Approach aka Integration) based on the Onsager formalism. The velocity autocorrelations are computed with FFT

References:
https://github.com/kdfong/transport-coefficients-MSD
https://github.com/kdfong/transport-coefficients
https://pubs.acs.org/doi/10.1021/acs.macromol.0c02001
https://arxiv.org/abs/2006.16164

PR Checklist

  • Tests?
  • Docs?
  • CHANGELOG updated?
  • Issue raised/referenced?

@pep8speaks
Copy link

pep8speaks commented Nov 25, 2024

Hello @rohithsrinivaas! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

Line 291:81: E501 line too long (87 > 80 characters)
Line 326:81: E501 line too long (83 > 80 characters)
Line 375:81: E501 line too long (88 > 80 characters)
Line 404:81: E501 line too long (84 > 80 characters)
Line 681:81: E501 line too long (82 > 80 characters)

Line 16:1: E266 too many leading '#' for block comment
Line 16:81: E501 line too long (89 > 80 characters)
Line 21:81: E501 line too long (86 > 80 characters)
Line 22:81: E501 line too long (127 > 80 characters)
Line 23:81: E501 line too long (103 > 80 characters)
Line 26:81: E501 line too long (82 > 80 characters)
Line 27:81: E501 line too long (83 > 80 characters)
Line 177:81: E501 line too long (84 > 80 characters)
Line 196:81: E501 line too long (101 > 80 characters)
Line 254:81: E501 line too long (88 > 80 characters)
Line 255:81: E501 line too long (86 > 80 characters)
Line 280:81: E501 line too long (87 > 80 characters)
Line 283:81: E501 line too long (87 > 80 characters)
Line 290:81: E501 line too long (83 > 80 characters)
Line 292:81: E501 line too long (81 > 80 characters)
Line 320:81: E501 line too long (89 > 80 characters)

Line 20:1: E266 too many leading '#' for block comment
Line 20:81: E501 line too long (89 > 80 characters)
Line 25:81: E501 line too long (86 > 80 characters)
Line 26:81: E501 line too long (127 > 80 characters)
Line 27:81: E501 line too long (103 > 80 characters)
Line 30:81: E501 line too long (82 > 80 characters)
Line 31:81: E501 line too long (83 > 80 characters)
Line 185:81: E501 line too long (84 > 80 characters)
Line 204:81: E501 line too long (101 > 80 characters)
Line 214:81: E501 line too long (84 > 80 characters)
Line 256:81: E501 line too long (86 > 80 characters)
Line 257:81: E501 line too long (84 > 80 characters)
Line 323:81: E501 line too long (183 > 80 characters)
Line 324:81: E501 line too long (218 > 80 characters)
Line 327:81: E501 line too long (83 > 80 characters)
Line 329:81: E501 line too long (82 > 80 characters)
Line 332:81: E501 line too long (82 > 80 characters)
Line 337:81: E501 line too long (81 > 80 characters)
Line 348:54: E203 whitespace before ':'
Line 348:81: E501 line too long (83 > 80 characters)
Line 350:46: E203 whitespace before ':'
Line 354:54: E203 whitespace before ':'
Line 354:81: E501 line too long (83 > 80 characters)
Line 356:46: E203 whitespace before ':'
Line 360:54: E203 whitespace before ':'
Line 360:81: E501 line too long (83 > 80 characters)
Line 362:46: E203 whitespace before ':'
Line 368:81: E501 line too long (81 > 80 characters)
Line 374:81: E501 line too long (87 > 80 characters)
Line 378:57: E203 whitespace before ':'
Line 378:81: E501 line too long (85 > 80 characters)
Line 381:58: E203 whitespace before ':'
Line 381:81: E501 line too long (86 > 80 characters)
Line 409:81: E501 line too long (81 > 80 characters)
Line 412:81: E501 line too long (103 > 80 characters)
Line 421:81: E501 line too long (87 > 80 characters)
Line 425:81: E501 line too long (88 > 80 characters)
Line 444:81: E501 line too long (89 > 80 characters)
Line 469:81: E501 line too long (84 > 80 characters)
Line 471:81: E501 line too long (91 > 80 characters)
Line 472:81: E501 line too long (186 > 80 characters)

Line 10:81: E501 line too long (84 > 80 characters)
Line 16:1: E266 too many leading '#' for block comment
Line 16:81: E501 line too long (89 > 80 characters)
Line 21:81: E501 line too long (86 > 80 characters)
Line 22:81: E501 line too long (127 > 80 characters)
Line 23:81: E501 line too long (103 > 80 characters)
Line 26:81: E501 line too long (82 > 80 characters)
Line 27:81: E501 line too long (83 > 80 characters)
Line 128:81: E501 line too long (85 > 80 characters)
Line 138:81: E501 line too long (84 > 80 characters)
Line 156:81: E501 line too long (101 > 80 characters)
Line 165:81: E501 line too long (84 > 80 characters)
Line 169:81: E501 line too long (86 > 80 characters)
Line 186:81: E501 line too long (83 > 80 characters)
Line 189:81: E501 line too long (82 > 80 characters)
Line 193:81: E501 line too long (83 > 80 characters)
Line 210:54: E203 whitespace before ':'
Line 210:81: E501 line too long (83 > 80 characters)
Line 211:57: E203 whitespace before ':'
Line 211:81: E501 line too long (86 > 80 characters)
Line 216:50: E203 whitespace before ':'
Line 220:58: E203 whitespace before ':'
Line 220:81: E501 line too long (86 > 80 characters)
Line 230:81: E501 line too long (107 > 80 characters)
Line 248:81: E501 line too long (89 > 80 characters)
Line 273:81: E501 line too long (84 > 80 characters)
Line 275:81: E501 line too long (91 > 80 characters)
Line 276:81: E501 line too long (186 > 80 characters)

Line 68:81: E501 line too long (81 > 80 characters)

Line 96:81: E501 line too long (84 > 80 characters)
Line 263:81: E501 line too long (88 > 80 characters)
Line 327:81: E501 line too long (84 > 80 characters)
Line 365:81: E501 line too long (88 > 80 characters)
Line 386:81: E501 line too long (86 > 80 characters)
Line 406:81: E501 line too long (86 > 80 characters)
Line 440:81: E501 line too long (81 > 80 characters)
Line 478:81: E501 line too long (88 > 80 characters)
Line 499:81: E501 line too long (83 > 80 characters)
Line 542:81: E501 line too long (87 > 80 characters)

Line 180:81: E501 line too long (82 > 80 characters)

Line 53:81: E501 line too long (84 > 80 characters)

Line 8:81: E501 line too long (81 > 80 characters)
Line 22:81: E501 line too long (94 > 80 characters)
Line 23:81: E501 line too long (125 > 80 characters)
Line 77:81: E501 line too long (94 > 80 characters)

Line 129:81: E501 line too long (84 > 80 characters)
Line 144:81: E501 line too long (83 > 80 characters)
Line 147:81: E501 line too long (84 > 80 characters)
Line 182:81: E501 line too long (87 > 80 characters)
Line 185:81: E501 line too long (85 > 80 characters)
Line 216:81: E501 line too long (84 > 80 characters)

Line 117:81: E501 line too long (83 > 80 characters)
Line 126:81: E501 line too long (84 > 80 characters)
Line 128:81: E501 line too long (83 > 80 characters)
Line 173:81: E501 line too long (87 > 80 characters)
Line 184:81: E501 line too long (85 > 80 characters)
Line 187:81: E501 line too long (83 > 80 characters)
Line 251:81: E501 line too long (82 > 80 characters)

Comment last updated at 2024-12-02 05:30:19 UTC

@rohithsrinivaas rohithsrinivaas changed the title Diffusion coeff Transport Properties using Onsager Formalism Nov 25, 2024
@orionarcher
Copy link
Collaborator

Would you mind linting this with something like black before I review?

@rohithsrinivaas
Copy link
Author

@orionarcher I have used black to lint as suggested. Kindly review the code and provide your comments

Copy link
Collaborator

@orionarcher orionarcher left a comment

Choose a reason for hiding this comment

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

Thank you @rohithsrinivaas! A of really awesome work here and I'm excited to have it in the package.

A few suggestions

  • Put both conductivity implementations in the same file and just call it conductivity.py
  • change diffusion_msd.py to diffusion.py
  • The function signatures need to be cleaned up in several places and the results object documented. I've mentioned this in my comments
  • Write tests for the new functionality. For that, I'd suggest adding a small ionic system to the test data . Just need to make sure everything can run end to end, it doesn't need to be scientifically rigorous, I trust you've done that already. solvation_analysis could be a good template for this. See the datafile and the contest.
  • There are more suggestions in the comments that I won't repeat here.


def plot_acf(self):
"""
Plot the mead square displacement vs time and the fit region in the log-log scale
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
Plot the mead square displacement vs time and the fit region in the log-log scale
Plot the mean square displacement vs time and the fit region in the log-log scale

Comment on lines +76 to +78
allatomgroup: "AtomGroup",
cationgroup_query: str,
aniongroup_query: str,
Copy link
Collaborator

Choose a reason for hiding this comment

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

I would not use queries for selecting anions and cations, just provide them directly

Suggested change
allatomgroup: "AtomGroup",
cationgroup_query: str,
aniongroup_query: str,
cations: "AtomGroup",
anions: "AtomGroup",

Comment on lines +81 to +84
cation_num_atoms_per_species=int,
anion_num_atoms_per_species=int,
cation_mass_per_species=float,
anion_mass_per_species=float,
Copy link
Collaborator

Choose a reason for hiding this comment

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

X_num_atoms_per_species can be calculated internally and masses are available at atoms.masses, no need to include any of these as kwargs

cations = # define cation atom group
cation_num_atoms_per_species = int(len(cations) / len(cations.residues))
# same for anions
cation_mass_per_species = 

Comment on lines +19 to +44
Class to calculate the diffusion_coefficient of a species.
Note that the slope of the mean square displacement provides
the diffusion coeffiecient. This implementation uses FFT to compute MSD faster
as explained here: https://stackoverflow.com/questions/34222272/computing-mean-square-displacement-using-python-and-fft
Since we are using FFT, the sampling frequency (or) dump frequency of the simulation trajectory
plays a critical role in the accuracy of the result.

Particle velocities are required to calculate the viscosity function. Thus
you are limited to MDA trajectories that contain velocity information, e.g.
GROMACS .trr files, H5MD files, etc. See the MDAnalysis documentation:
https://docs.mdanalysis.org/stable/documentation_pages/coordinates/init.html#writers.

Parameters
----------
atomgroup : AtomGroup
An MDAnalysis :class:`~MDAnalysis.core.groups.AtomGroup`.
Note that :class:`UpdatingAtomGroup` instances are not accepted.
temp_avg : float (optional, default 300)
Average temperature over the course of the simulation, in Kelvin.
dim_type : {'xyz', 'xy', 'yz', 'xz', 'x', 'y', 'z'}
Desired dimensions to be included in the viscosity calculation.
Defaults to 'xyz'.
linear_fit_window : tuple of int (optional)
A tuple of two integers specifying the start and end lag-time for
the linear fit of the viscosity function. If not provided, the
linear fit is not performed and viscosity is not calculated.
Copy link
Collaborator

Choose a reason for hiding this comment

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

Indentation is funky here, look at another class for an example

Comment on lines +85 to +86
cation_charge=int,
anion_charge=int,
Copy link
Collaborator

Choose a reason for hiding this comment

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

Include defaults

Suggested change
cation_charge=int,
anion_charge=int,
cation_charge: int=1,
anion_charge: int=-1,

Comment on lines +52 to +61
results.timeseries : :class:`numpy.ndarray`
The averaged viscosity function over all the particles with respect
to lag-time. Obtained after calling :meth:`ViscosityHelfand.run`
results.visc_by_particle : :class:`numpy.ndarray`
The viscosity function of each individual particle with respect to
lag-time.
results.viscosity : float
The viscosity coefficient of the solvent. The argument
`linear_fit_window` must be provided to calculate this to
avoid misinterpretation of the viscosity function.
Copy link
Collaborator

Choose a reason for hiding this comment

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

I don't think there are correct for this class

Comment on lines +50 to +51
dim_fac : int
Dimensionality :math:`d` of the viscosity computation.
Copy link
Collaborator

Choose a reason for hiding this comment

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

Not viscosity here

Comment on lines +427 to +440
self.results.cation_transference_com = (
(self.cation_charge**2) * lij_cation_cation
+ (self.cation_charge * self.anion_charge) * lij_cation_anion
) / cond
self.results.anion_transference_com = (
(self.anion_charge**2) * lij_anion_anion
+ (self.cation_charge * self.anion_charge) * lij_cation_anion
) / cond
self.results.cation_transference_solvent = (
self.results.cation_transference_com - anion_mass_fraction
) / solvent_fraction
self.results.anion_transference_solvent = (
self.results.anion_transference_com - cation_mass_fraction
) / solvent_fraction
Copy link
Collaborator

Choose a reason for hiding this comment

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

The contents of results should be documented here. It's not clear to me what the difference between the _solvent and _com quantities is?

Comment on lines +52 to +61
results.timeseries : :class:`numpy.ndarray`
The averaged viscosity function over all the particles with respect
to lag-time. Obtained after calling :meth:`ViscosityHelfand.run`
results.visc_by_particle : :class:`numpy.ndarray`
The viscosity function of each individual particle with respect to
lag-time.
results.viscosity : float
The viscosity coefficient of the solvent. The argument
`linear_fit_window` must be provided to calculate this to
avoid misinterpretation of the viscosity function.
Copy link
Collaborator

Choose a reason for hiding this comment

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

Document results signature

Comment on lines +76 to +77
allatomgroup: "AtomGroup",
atomgroup_query: str,
Copy link
Collaborator

Choose a reason for hiding this comment

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

Not clear why we need a query here?

Suggested change
allatomgroup: "AtomGroup",
atomgroup_query: str,
atomgroup: "AtomGroup",

Would need to be changed below too.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Add conductivity impementation
3 participants