diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 5589dea..07fe8d8 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -1,4 +1,4 @@ # The following owners will be the default owners for everything # in the repo unless a later match takes precedence. -sportran/* @lorisercole -sportran_gui/* @rikigigi +sportran/ @lorisercole +sportran_gui/ @rikigigi diff --git a/MANIFEST.in b/MANIFEST.in index 754f4f9..55dc835 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -3,9 +3,9 @@ include LICENSE.txt include setup.json include sportran/README.md include sportran/metadata.json -include sportran/utils/plot_style.mplstyle -exclude sportran/utils/blocks.py -exclude sportran/utils/blockanalysis.py +exclude sportran/utils/obsolete/ include sportran_gui/README_GUI.md include sportran_gui/assets/icon.gif include sportran_gui/assets/languages.json +include sportran/plotter/styles/api_style.mplstyle +include sportran/plotter/styles/cli_style.mplstyle diff --git a/README.md b/README.md index 3489a78..5e38ba4 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,7 @@ A code to estimate transport coefficients from the cepstral analysis of a multi- https://sportran.readthedocs.io ### References - - [Ercole L., Bertossa R., Bisacchi S., and Baroni S., "_SporTran: a code to estimate transport coefficients from the cepstral analysis of (multivariate) current time series_", *arXiv*:2202.11571 (2022)](https://arxiv.org/abs/2202.11571), submitted to *Comput. Phys. Commun.* + - [Ercole L., Bertossa R., Bisacchi S., and Baroni S., "_SporTran: a code to estimate transport coefficients from the cepstral analysis of (multivariate) current time series_", *Comput. Phys. Commun.*, 108470](https://doi.org/10.1016/j.cpc.2022.108470), [*arXiv*:2202.11571 (2022)](https://arxiv.org/abs/2202.11571) - (cepstral analysis) [Ercole, Marcolongo, Baroni, *Sci. Rep.* **7**, 15835 (2017)](https://doi.org/10.1038/s41598-017-15843-2) - (multicomponent systems) [Bertossa, Grasselli, Ercole, Baroni, *Phys. Rev. Lett.* **122**, 255901 (2019)](https://doi.org/10.1103/PhysRevLett.122.255901) ([arXiv](https://arxiv.org/abs/1808.03341)) - (review) [Baroni, Bertossa, Ercole, Grasselli, Marcolongo, *Handbook of Materials Modeling* (2018)](https://doi.org/10.1007/978-3-319-50257-1_12-1) ([arXiv](https://arxiv.org/abs/1802.08006)) diff --git a/setup.json b/setup.json index f3b1bf0..63076bc 100644 --- a/setup.json +++ b/setup.json @@ -1,6 +1,6 @@ { "name": "sportran", - "version": "1.0.0rc2", + "version": "1.0.0rc3", "author": "Loris Ercole, Riccardo Bertossa, Sebastiano Bisacchi", "author_email": "loris.ercole@epfl.ch", "description": "Cepstral Data Analysis of current time series for Green-Kubo transport coefficients", diff --git a/sportran/analysis.py b/sportran/analysis.py index 8c69506..b0c44d3 100755 --- a/sportran/analysis.py +++ b/sportran/analysis.py @@ -211,6 +211,18 @@ def main(): return 0 +def concatenate_if_not_none_with_labels(concat, labels=None): + out_arr = [] + out_label = '' + if labels is None: + labels = ['' for i in concat] + for arr, label in zip(concat, labels): + if arr is not None: + out_arr.append(arr) + out_label += f' {label}' + return np.concatenate([out_arr], axis=1).transpose(), f'{out_label}\n' + + def run_analysis(args): inputfile = args.inputfile @@ -472,8 +484,8 @@ def average(data, name, units=''): if not no_text_out: outfile_name = output + '.psd.dat' - outarray = np.c_[j.freqs_THz, j.psd, j.fpsd, j.logpsd, j.flogpsd] - outfile_header = 'freqs_THz psd fpsd logpsd flogpsd\n' + outarray, outfile_header = concatenate_if_not_none_with_labels( + [j.freqs_THz, j.psd, j.fpsd, j.logpsd, j.flogpsd], ['freqs_THz', 'psd', 'fpsd', 'logpsd', 'flogpsd']) np.savetxt(outfile_name, outarray, header=outfile_header, fmt=fmt) if j.MANY_CURRENTS: outfile_name = output + '.cospectrum.dat' @@ -483,15 +495,17 @@ def average(data, name, units=''): np.savetxt(outfile_name, np.column_stack([outarray.real, outarray.imag]), fmt=fmt) outfile_name = output + '.cospectrum.filt.dat' - outarray = np.c_[j.freqs_THz, - j.fcospectrum.reshape( - (j.fcospectrum.shape[0] * j.fcospectrum.shape[1], j.fcospectrum.shape[2])).transpose()] - np.savetxt(outfile_name, np.column_stack([outarray.real, outarray.imag]), fmt=fmt) + if j.fcospectrum is not None: + outarray = np.c_[j.freqs_THz, + j.fcospectrum.reshape((j.fcospectrum.shape[0] * j.fcospectrum.shape[1], + j.fcospectrum.shape[2])).transpose()] + np.savetxt(outfile_name, np.column_stack([outarray.real, outarray.imag]), fmt=fmt) if resample: outfile_name = output + '.resampled_psd.dat' - outarray = np.c_[jf.freqs_THz, jf.psd, jf.fpsd, jf.logpsd, jf.flogpsd] - outfile_header = 'freqs_THz psd fpsd logpsd flogpsd\n' + outarray, outfile_header = concatenate_if_not_none_with_labels( + [jf.freqs_THz, jf.psd, jf.fpsd, jf.logpsd, jf.flogpsd], + ['freqs_THz', 'psd', 'fpsd', 'logpsd', 'flogpsd']) np.savetxt(outfile_name, outarray, header=outfile_header, fmt=fmt) outfile_name = output + '.cepstral.dat' @@ -623,7 +637,8 @@ def write_old_binary(self, output): """Write old binary format.""" opts = {'allow_pickle': False} optsa = {'axis': 1} - outarray = np.c_[self.j_freqs_THz, self.j_fpsd, self.j_flogpsd, self.j_psd, self.j_logpsd] + outarray, _ = concatenate_if_not_none_with_labels( + [self.j_freqs_THz, self.j_fpsd, self.j_flogpsd, self.j_psd, self.j_logpsd]) np.save(output + '.psd.npy', outarray, **opts) if self.j_cospectrum is not None: @@ -634,7 +649,8 @@ def write_old_binary(self, output): outarray = np.c_[self.j_freqs_THz, self.j_fcospectrum.reshape(-1, self.j_fcospectrum.shape[-1]).transpose()] np.save(output + '.cospectrum.filt.npy', outarray, **opts) - outarray = np.c_[self.jf_freqs_THz, self.jf_psd, self.jf_fpsd, self.jf_logpsd, self.jf_flogpsd] + outarray, _ = concatenate_if_not_none_with_labels( + [self.jf_freqs_THz, self.jf_psd, self.jf_fpsd, self.jf_logpsd, self.jf_flogpsd]) np.save(output + '.resampled_psd.npy', outarray, **opts) outarray = np.c_[self.jf_cepf_logpsdK, self.jf_cepf_logpsdK_THEORY_std, self.jf_cepf_logtau, diff --git a/sportran/md/mdsample.py b/sportran/md/mdsample.py index a3c44a4..f8acaec 100644 --- a/sportran/md/mdsample.py +++ b/sportran/md/mdsample.py @@ -320,7 +320,7 @@ def filter_psd(self, PSD_FILTER_W=None, freq_units='THz', window_type='rectangul """ Filter the periodogram with the given PSD_FILTER_W [freq_units]. - PSD_FILTER_W PSD filter window [freq_units] - - freq_units frequency units ['THz', 'red' (default)] + - freq_units frequency units ['THz' (default), 'red'] - window_type filtering window type ['rectangular'] """ if self.psd is None: diff --git a/sportran/md/tools/filter.py b/sportran/md/tools/filter.py index 557b1ff..3806260 100644 --- a/sportran/md/tools/filter.py +++ b/sportran/md/tools/filter.py @@ -13,5 +13,10 @@ def runavefilter(X, WF): WF = WF + 1 W = int(WF / 2) + if W > X.shape[0] - 1: + W = X.shape[0] - 1 + WF = 2 * W + print('Warning: reducing filtering window') + Y = np.concatenate((X[W:0:-1], X, X[-2:-W - 2:-1])) return np.convolve(Y, np.array([1.0 / WF] * WF), 'valid') diff --git a/tests/test_cli.py b/tests/test_cli.py index 9b6dc5b..8a047d6 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -35,6 +35,27 @@ def test_cli_NaCl(tmpdir, run_cli, data_NaCl_path, num_regression, file_regressi num_regression.check(readed) +def test_cli_NaCl_no_w(tmpdir, run_cli, data_NaCl_path, num_regression, file_regression): + fileout_pref = tmpdir.join('output_no_w') + output = run_cli(data_NaCl_path, '-k', 'flux', '-j', 'vcm[1]', '-t', '5.0', '--VOLUME', '65013.301261', + '--param-from-input-file-column', 'Temp', 'TEMPERATURE', '--FSTAR', '14.0', '-r', + '--test-suite-run', '-o', fileout_pref) + file_list = [ + 'output_no_w.cepstral.dat', 'output_no_w.cepstrumfiltered_psd.dat', 'output_no_w.cospectrum.dat', + 'output_no_w.psd.dat', 'output_no_w.resampled_psd.dat' + ] + with open(str(tmpdir.join('output_no_w.log'))) as l: + file_regression.check(l.read(), basename='output_no_w.log') + file_regression.check(output.stdout.str(), basename='stdout_no_w') + file_regression.check(output.stderr.str(), basename='stderr_no_w') + readed = {} + for f in file_list: + file_out = tmpdir.join(f) + readed[f] = np.loadtxt(file_out).flatten() + + num_regression.check(readed) + + def test_cli_bin_output_NaCl(tmpdir, run_cli, data_NaCl_path, num_regression, file_regression): fileout_pref = tmpdir.join('output') output = run_cli(data_NaCl_path, '-k', 'flux', '-j', 'vcm[1]', '-t', '5.0', '--VOLUME', '65013.301261', diff --git a/tests/test_cli/output_no_w.log.txt b/tests/test_cli/output_no_w.log.txt new file mode 100644 index 0000000..7b7e98a --- /dev/null +++ b/tests/test_cli/output_no_w.log.txt @@ -0,0 +1,66 @@ + Units: metal + Time step: 5.0 fs +# Molten NaCl - Fumi-Tosi potential +# https://journals.aps.org/prl/supplemental/10.1103/PhysRevLett.122.255901/Bertossa_et_al_Supplemental_Material.pdf +# 1728 atoms, T~1400K +# NVE, dt = 1.0 fs, 100 ps, print_step = 5.0 fs +# Volume = 65013.301261 A^3 +# LAMMPS metal units +Temp c_flux[1] c_flux[2] c_flux[3] c_vcm[1][1] c_vcm[1][2] c_vcm[1][3] + ##################################### + all_ckeys = [('Temp', [0]), ('flux', array([1, 2, 3])), ('vcm[1]', array([4, 5, 6]))] + ##################################### +Data length = 20000 + ckey = [('Temp', [0]), ('flux', array([1, 2, 3])), ('vcm[1]', array([4, 5, 6]))] + ( 20000 ) steps read. +VOLUME (input): 65013.301261 +Mean TEMPERATURE (computed): 1399.3477811999999 +/- 19.318785820942594 + Time step (input): 5.0 fs + currents shape is (2, 20000, 3) +snippet: +[[[ 2.5086549e+02 2.0619423e+01 2.0011500e+02] + [ 1.9622265e+02 8.2667342e+01 2.8433250e+02] + [ 1.2639441e+02 1.6075472e+02 3.4036769e+02] + ... + [ 1.7991856e+02 1.8612706e+01 -1.3265623e+02] + [ 2.0471193e+02 -4.6643315e-01 -2.0401650e+02] + [ 2.4123318e+02 -1.8295461e+01 -2.5246475e+02]] + + [[-1.5991832e-01 -7.1370426e-02 2.0687917e-02] + [-1.3755206e-01 -7.1002931e-02 -1.1279876e-02] + [-1.0615044e-01 -6.2381243e-02 -4.1568120e-02] + ... + [-9.1939899e-02 -8.4778292e-02 6.0011385e-02] + [-1.3384949e-01 -1.1474530e-01 8.9323546e-02] + [-1.8385053e-01 -1.3693430e-01 1.1434060e-01]]] +Using multicomponent code. + Number of currents = 2 + Number of equivalent components = 3 + KAPPA_SCALE = 1.4604390788939313e-07 + Nyquist_f = 100.0 THz +Using multicomponent code. +----------------------------------------------------- + RESAMPLE TIME SERIES +----------------------------------------------------- + Original Nyquist freq f_Ny = 100.00000 THz + Resampling freq f* = 14.28571 THz + Sampling time TSKIP = 7 steps + = 35.000 fs + Original n. of frequencies = 10001 + Resampled n. of frequencies = 1429 + min(PSD) (pre-filter&sample) = 0.31536 + min(PSD) (post-filter&sample) = 22166.11934 + % of original PSD Power f