Skip to content

Commit

Permalink
Added new data for unit tests; Tweaked stitch_time_series and added i…
Browse files Browse the repository at this point in the history
…ts unit test
  • Loading branch information
wehs7661 committed Apr 8, 2024
1 parent dad30e3 commit 8376dcd
Show file tree
Hide file tree
Showing 7 changed files with 271 additions and 42 deletions.
21 changes: 5 additions & 16 deletions ensemble_md/analysis/analyze_traj.py
Original file line number Diff line number Diff line change
Expand Up @@ -94,32 +94,21 @@ def stitch_time_series(files, rep_trajs, shifts=None, dhdl=True, col_idx=-1, sav
files_sorted[i].append(files[rep_trajs[i][j]][j])

# Then, stitch the trajectories for each starting configuration
# Unlike stitch_time_series_for_sim, there is no way to check the continuity.
trajs = [[] for i in range(n_configs)] # for each starting configuration
t_last, val_last = None, None # just for checking the continuity of the trajectory
for i in range(n_configs):
for j in range(n_iter):
if dhdl:
traj, t = extract_state_traj(files_sorted[i][j])
traj, _ = extract_state_traj(files_sorted[i][j])
else:
traj = np.loadtxt(files_sorted[i][j], comments=['#', '@'])[:, col_idx]
t = np.loadtxt(files_sorted[i][j], comments=['#', '@'])[:, 0]

# Shift the indices so that global indices are used.
shift_idx = rep_trajs[i][j]
traj = list(np.array(traj) + shifts[shift_idx])

if j != 0:
# Check the continuity of the trajectory
if traj[0] != val_last or t[0] != t_last:
err_str = f'The first frame of iteration {j} of starting configuration {i} is not continuous with the last frame of the previous iteration. ' # noqa: E501
err_str += f'Please check files {files_sorted[i][j - 1]} and {files_sorted[i][j]}.'
raise ValueError(err_str)

t_last = t[-1]
val_last = traj[-1]

if j != 0:
traj = traj[:-1] # remove the last frame, which is the same as the first of the next time series.
if j != n_iter - 1:
traj = traj[:-1]

trajs[i].extend(traj)

Expand Down Expand Up @@ -200,7 +189,7 @@ def stitch_time_series_for_sim(files, dhdl=True, col_idx=-1, save=True):
return trajs


def stitch_trajs(gmx_executable, files, rep_trajs):
def stitch_xtc_trajs(gmx_executable, files, rep_trajs):
"""
Demuxes GROMACS trajectories from different replicas into individual continuous trajectories.
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
# This file was created Mon Jun 19 04:08:44 2023
# Created by:
# :-) GROMACS - gmx mdrun, 2022.5-dev-20230428-fdf57150ad (-:
#
# Executable: /jet/home/wehs7661/pkgs/gromacs/2022.5/bin/gmx
# Data prefix: /jet/home/wehs7661/pkgs/gromacs/2022.5
# Working dir: /ocean/projects/cts160011p/wehs7661/EEXE_experiments/Anthracene/EEXE/fixed_weight/Redo/test_7/sim_0/iteration_2
# Command line:
# gmx mdrun -s sys_EE.tpr -nt 16 -ntmpi 1
# gmx mdrun is part of G R O M A C S:
#
# GROup of MAchos and Cynical Suckers
#
@ title "dH/d\xl\f{} and \xD\f{}H"
@ xaxis label "Time (ps)"
@ yaxis label "dH/d\xl\f{} and \xD\f{}H (kJ/mol [\xl\f{}]\S-1\N)"
@TYPE xy
@ subtitle "T = 300 (K) "
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "Thermodynamic state"
@ s1 legend "Total Energy (kJ/mol)"
@ s2 legend "dH/d\xl\f{} vdw-lambda = 0.1800"
@ s3 legend "\xD\f{}H \xl\f{} to 0.0000"
@ s4 legend "\xD\f{}H \xl\f{} to 0.1800"
@ s5 legend "\xD\f{}H \xl\f{} to 0.4200"
@ s6 legend "\xD\f{}H \xl\f{} to 0.5700"
@ s7 legend "\xD\f{}H \xl\f{} to 0.6800"
8.0000 1 -34629.660 76.294617 -12.171330 0.0000000 20.012351 33.639191 43.923604
8.2000 1 -34651.223 53.820854 -5.7918213 0.0000000 16.636073 29.299222 39.133190
8.4000 1 -34551.805 76.754662 -12.062535 0.0000000 20.355419 34.357167 44.960779
8.6000 0 -34470.441 60.429405 0.0000000 12.838254 33.812664 48.070257 58.823763
8.8000 1 -34470.359 50.146221 -5.6069839 0.0000000 15.713153 27.967770 37.615367
9.0000 1 -34498.391 39.090370 -2.5533874 0.0000000 13.934211 25.543530 34.808979
9.2000 1 -34556.098 84.098511 -14.461316 0.0000000 20.903842 34.415566 44.420412
9.4000 0 -34616.480 70.126610 0.0000000 13.909344 35.086697 49.137125 59.653158
9.6000 1 -34534.547 26.037083 1.3580391 0.0000000 12.167194 23.478886 32.722436
9.8000 2 -34312.969 71.234009 -19.139816 -14.436126 0.0000000 11.326751 20.238176
10.0000 0 -34378.742 76.646988 0.0000000 14.550442 35.442341 48.924097 58.903971
10.2000 2 -34285.172 84.664841 -27.865017 -18.636821 0.0000000 13.043783 22.901348
10.4000 1 -34497.484 63.202911 -8.3794296 0.0000000 18.313194 31.802937 42.233846
10.6000 1 -34578.695 87.972961 -15.083676 0.0000000 21.971679 36.270589 46.895210
10.8000 0 -34541.055 57.429951 0.0000000 12.311475 32.571461 46.366842 56.775666
11.0000 0 -34542.328 54.882725 0.0000000 12.163686 32.820439 47.037683 57.802180
11.2000 1 -34556.328 73.172264 -10.645428 0.0000000 20.349487 34.931084 46.128175
11.4000 0 -34620.914 66.442963 0.0000000 13.768872 35.857712 50.844045 62.154873
11.6000 1 -34630.789 87.372787 -14.095708 0.0000000 22.653880 37.868951 49.277689
11.8000 0 -34617.402 84.743134 0.0000000 16.467885 40.928452 56.989053 68.961985
12.0000 1 -34608.320 -64.654449 36.570163 0.0000000 2.2114481 12.925400 22.739461
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
# This file was created Mon Jun 19 04:08:44 2023
# Created by:
# :-) GROMACS - gmx mdrun, 2022.5-dev-20230428-fdf57150ad (-:
#
# Executable: /jet/home/wehs7661/pkgs/gromacs/2022.5/bin/gmx
# Data prefix: /jet/home/wehs7661/pkgs/gromacs/2022.5
# Working dir: /ocean/projects/cts160011p/wehs7661/EEXE_experiments/Anthracene/EEXE/fixed_weight/Redo/test_7/sim_1/iteration_2
# Command line:
# gmx mdrun -s sys_EE.tpr -nt 16 -ntmpi 1
# gmx mdrun is part of G R O M A C S:
#
# GROup of MAchos and Cynical Suckers
#
@ title "dH/d\xl\f{} and \xD\f{}H"
@ xaxis label "Time (ps)"
@ yaxis label "dH/d\xl\f{} and \xD\f{}H (kJ/mol [\xl\f{}]\S-1\N)"
@TYPE xy
@ subtitle "T = 300 (K) "
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "Thermodynamic state"
@ s1 legend "Total Energy (kJ/mol)"
@ s2 legend "dH/d\xl\f{} vdw-lambda = 0.1800"
@ s3 legend "\xD\f{}H \xl\f{} to 0.1800"
@ s4 legend "\xD\f{}H \xl\f{} to 0.4200"
@ s5 legend "\xD\f{}H \xl\f{} to 0.5700"
@ s6 legend "\xD\f{}H \xl\f{} to 0.6800"
@ s7 legend "\xD\f{}H \xl\f{} to 0.7600"
8.0000 0 -34755.293 46.771725 0.0000000 15.872268 28.695901 38.829498 46.447546
8.2000 0 -34785.953 74.836411 0.0000000 20.020610 33.922478 44.493968 52.312291
8.4000 0 -34647.559 56.015610 0.0000000 17.202335 30.347924 40.604961 48.276854
8.6000 0 -34641.148 -4.2277699 0.0000000 8.7699007 20.164078 29.975124 37.611529
8.8000 1 -34616.250 69.080246 -10.322360 0.0000000 11.894403 22.082836 30.019714
9.0000 1 -34614.609 42.226288 0.41087019 0.0000000 8.9144948 17.931109 25.390445
9.2000 1 -34468.852 47.732052 -2.6941484 0.0000000 9.2777303 18.100653 25.243219
9.4000 0 -34504.332 55.423943 0.0000000 16.886680 29.825406 39.961572 47.562781
9.6000 1 -34515.789 59.769459 -6.4545549 0.0000000 10.655896 19.998678 27.305587
9.8000 1 -34515.922 -3.8737392 21.757142 0.0000000 3.6053060 9.9567376 15.796676
10.0000 2 -34563.215 -34.931252 159.09013 15.563947 0.0000000 -0.40763683 2.4860333
10.2000 3 -34533.957 -2.8320389 295.73599 31.136484 4.1779768 0.0000000 1.2499733
10.4000 3 -34558.859 -2.1180019 280.30180 30.889294 4.1170456 0.0000000 1.3188639
10.6000 4 -34530.309 -30.518177 751.66180 96.462129 24.109847 5.0418812 0.0000000
10.8000 4 -34537.238 -71.261162 930.65434 133.39602 36.839697 9.0098523 0.0000000
11.0000 4 -34461.461 -87.927658 1239.1910 157.80551 43.247558 10.757349 0.0000000
11.2000 4 -34466.941 -82.624626 889.63010 142.10997 40.461601 10.151953 0.0000000
11.4000 4 -34466.344 -123.39691 1340.7995 193.15680 55.609261 14.434419 0.0000000
11.6000 4 -34515.422 -127.27545 1497.2377 203.23258 57.821370 14.932713 0.0000000
11.8000 4 -34671.367 -175.34641 1967.8005 259.49930 74.645520 19.813017 0.0000000
12.0000 4 -34669.914 -206.61479 2211.1280 291.95570 84.940968 22.912667 0.0000000
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
# This file was created Mon Jun 19 04:08:44 2023
# Created by:
# :-) GROMACS - gmx mdrun, 2022.5-dev-20230428-fdf57150ad (-:
#
# Executable: /jet/home/wehs7661/pkgs/gromacs/2022.5/bin/gmx
# Data prefix: /jet/home/wehs7661/pkgs/gromacs/2022.5
# Working dir: /ocean/projects/cts160011p/wehs7661/EEXE_experiments/Anthracene/EEXE/fixed_weight/Redo/test_7/sim_2/iteration_2
# Command line:
# gmx mdrun -s sys_EE.tpr -nt 16 -ntmpi 1
# gmx mdrun is part of G R O M A C S:
#
# GROup of MAchos and Cynical Suckers
#
@ title "dH/d\xl\f{} and \xD\f{}H"
@ xaxis label "Time (ps)"
@ yaxis label "dH/d\xl\f{} and \xD\f{}H (kJ/mol [\xl\f{}]\S-1\N)"
@TYPE xy
@ subtitle "T = 300 (K) "
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "Thermodynamic state"
@ s1 legend "Total Energy (kJ/mol)"
@ s2 legend "dH/d\xl\f{} vdw-lambda = 0.5700"
@ s3 legend "\xD\f{}H \xl\f{} to 0.4200"
@ s4 legend "\xD\f{}H \xl\f{} to 0.5700"
@ s5 legend "\xD\f{}H \xl\f{} to 0.6800"
@ s6 legend "\xD\f{}H \xl\f{} to 0.7600"
@ s7 legend "\xD\f{}H \xl\f{} to 0.8600"
8.0000 1 -34313.500 53.922421 -4.9744428 0.0000000 7.0177010 13.155697 21.633762
8.2000 1 -34317.035 54.985249 -5.1219533 0.0000000 7.1659649 13.451778 22.153823
8.4000 0 -34343.859 64.260559 0.0000000 11.269674 20.979761 28.515765 38.314333
8.6000 0 -34327.445 42.952110 0.0000000 8.9256910 17.692659 24.808122 34.303222
8.8000 0 -34193.602 30.741173 0.0000000 7.3090871 15.119509 21.632564 30.461708
9.0000 1 -34165.562 17.672602 3.5807292 0.0000000 4.0475147 8.9581755 16.590528
9.2000 2 -34196.078 53.809811 -1.9005395 -4.2808140 0.0000000 4.9581900 12.597696
9.4000 1 -34121.391 30.466057 -0.028901401 0.0000000 4.9883106 10.199266 17.972795
9.6000 2 -34181.188 12.930644 22.614589 1.8607628 0.0000000 2.2793522 7.7807644
9.8000 2 -34159.910 -83.369080 82.336740 16.540114 0.0000000 -3.9512120 -3.2089027
10.0000 3 -34253.008 -79.733406 147.43616 40.507554 9.9642161 0.0000000 -3.8256816
10.2000 3 -34316.266 -100.41053 169.57082 47.693830 12.082289 0.0000000 -5.3537159
10.4000 3 -34287.016 -38.847527 117.40750 29.041785 6.1480040 0.0000000 -0.39827370
10.6000 3 -34357.711 -7.0617037 73.913992 16.967344 2.8104849 0.0000000 1.9357874
10.8000 2 -34354.020 34.890873 8.7370223 -1.3708427 0.0000000 3.7911577 10.761747
11.0000 1 -34391.375 -18.832100 11.990726 0.0000000 1.0368603 4.7231291 11.600601
11.2000 2 -34355.613 53.997337 -1.6474591 -4.2284115 0.0000000 5.0100889 12.806016
11.4000 1 -34464.656 -3.4788733 8.4742171 0.0000000 2.2977398 6.4826487 13.639728
11.6000 1 -34486.930 -16.786358 10.481138 0.0000000 1.0813492 4.7300642 11.606537
11.8000 2 -34525.945 -17.548145 39.956987 6.5848033 0.0000000 0.41567518 4.8993547
12.0000 2 -34481.375 -113.31747 102.25158 21.396274 0.0000000 -5.7530722 -6.0001237
Binary file not shown.
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
# This file was created Mon Jun 19 04:08:44 2023
# Created by:
# :-) GROMACS - gmx mdrun, 2022.5-dev-20230428-fdf57150ad (-:
#
# Executable: /jet/home/wehs7661/pkgs/gromacs/2022.5/bin/gmx
# Data prefix: /jet/home/wehs7661/pkgs/gromacs/2022.5
# Working dir: /ocean/projects/cts160011p/wehs7661/EEXE_experiments/Anthracene/EEXE/fixed_weight/Redo/test_7/sim_3/iteration_2
# Command line:
# gmx mdrun -s sys_EE.tpr -nt 16 -ntmpi 1
# gmx mdrun is part of G R O M A C S:
#
# GROup of MAchos and Cynical Suckers
#
@ title "dH/d\xl\f{} and \xD\f{}H"
@ xaxis label "Time (ps)"
@ yaxis label "dH/d\xl\f{} and \xD\f{}H (kJ/mol [\xl\f{}]\S-1\N)"
@TYPE xy
@ subtitle "T = 300 (K) "
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "Thermodynamic state"
@ s1 legend "Total Energy (kJ/mol)"
@ s2 legend "dH/d\xl\f{} vdw-lambda = 0.8600"
@ s3 legend "\xD\f{}H \xl\f{} to 0.5700"
@ s4 legend "\xD\f{}H \xl\f{} to 0.6800"
@ s5 legend "\xD\f{}H \xl\f{} to 0.7600"
@ s6 legend "\xD\f{}H \xl\f{} to 0.8600"
@ s7 legend "\xD\f{}H \xl\f{} to 1.0000"
8.0000 3 -34784.020 -77.079910 103.68410 38.166362 14.101582 0.0000000 -2.3011742
8.2000 3 -34860.699 -58.251957 86.461296 31.417711 11.356713 0.0000000 -0.75197034
8.4000 3 -34859.223 -16.553837 54.195003 18.130927 5.7016413 0.0000000 3.2131857
8.6000 2 -34805.293 -131.15421 58.714853 15.259576 0.0000000 -7.6177195 -5.4159167
8.8000 3 -34891.207 -74.783188 98.909443 36.761498 13.646595 0.0000000 -2.1796022
9.0000 3 -34795.844 -85.563187 106.10056 39.968034 15.062222 0.0000000 -3.2385274
9.2000 3 -34817.168 -62.806141 89.513820 32.832173 11.981368 0.0000000 -1.1259146
9.4000 2 -34871.062 -132.36244 57.814263 15.243424 0.0000000 -7.8148593 -5.9443458
9.6000 2 -34903.781 -51.724876 31.809730 7.2676394 0.0000000 -1.5270443 4.9536180
9.8000 2 -34873.547 -65.990837 37.090620 8.7470559 0.0000000 -2.6001132 3.1049916
10.0000 2 -34842.301 -68.331505 36.699844 8.8441776 0.0000000 -2.8967041 2.3632594
10.2000 2 -34814.570 -28.212650 23.207305 4.8157071 0.0000000 0.18411137 7.6140201
10.4000 1 -34842.922 -17.237011 6.6958868 0.0000000 0.43287851 4.8234432 15.423047
10.6000 1 -34889.977 -50.003258 11.445712 0.0000000 -1.7323375 0.94870832 10.348200
10.8000 2 -34907.156 -159.88383 67.373923 17.997990 0.0000000 -9.9937868 -9.9714539
11.0000 3 -34936.309 -77.132729 98.832614 36.886646 13.790942 0.0000000 -2.6671138
11.2000 3 -34800.629 -68.980125 94.532006 34.906907 12.851202 0.0000000 -1.6614399
11.4000 3 -34740.562 -113.63854 129.12077 49.226807 18.943248 0.0000000 -5.8271056
11.6000 4 -34780.863 11.287065 139.94837 58.477653 27.002308 6.8633881 0.0000000
11.8000 3 -34739.133 -155.06459 154.36888 60.818883 24.201765 0.0000000 -10.067938
12.0000 4 -34730.758 -16.987450 183.54891 80.260639 39.675461 12.434976 0.0000000
80 changes: 54 additions & 26 deletions ensemble_md/tests/test_analyze_traj.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,60 @@ def test_extract_state_traj():


def test_stitch_time_series():
folder = os.path.join(input_path, 'dhdl/simulation_example')
files = [[f'{folder}/sim_{i}/iteration_{j}/dhdl.xvg' for j in range(3)] for i in range(4)]
rep_trajs = np.array([[0, 0, 1], [1, 1, 0], [2, 2, 2], [3, 3, 3]])
shifts = [0, 1, 2, 3]

trajs = analyze_traj.stitch_time_series(files, rep_trajs, shifts)
assert trajs[0] == [
0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1,
1, 1, 1, 2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 2, 2, 1, 0, 1, 1,
1, 1, 1, 1, 2, 2, 2, 1, 2, 2, 3, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5
]
assert trajs[1] == [
1, 1, 2, 3, 3, 3, 2, 2, 1, 1, 1, 1, 1, 2, 3, 2, 1, 1, 1, 1,
2, 2, 1, 1, 1, 1, 1, 2, 3, 2, 1, 1, 1, 1, 2, 3, 3, 3, 2, 2,
1, 1, 1, 0, 1, 1, 1, 0, 1, 2, 0, 2, 1, 1, 0, 0, 1, 0, 1, 0, 1
]
assert trajs[2] == [
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 2, 2, 2, 2, 3, 3,
3, 3, 3, 3, 3, 3, 3, 2, 3, 2, 3, 3, 3, 2, 2, 3, 4, 3, 3, 2,
3, 3, 2, 2, 2, 3, 4, 3, 4, 4, 5, 5, 5, 5, 4, 3, 4, 3, 3, 4, 4
]
assert trajs[3] == [
3, 3, 3, 3, 3, 3, 3, 5, 4, 4, 5, 4, 4, 5, 4, 5, 5, 5, 4, 5,
4, 4, 5, 4, 5, 5, 4, 5, 5, 5, 4, 5, 5, 4, 5, 4, 5, 4, 5, 5,
6, 6, 6, 5, 6, 6, 6, 5, 5, 5, 5, 5, 4, 4, 5, 6, 6, 6, 7, 6, 7
]

assert os.path.exists('state_trajs.npy')
os.remove('state_trajs.npy')


def test_stitch_time_series_for_sim():
# Set up files for testing
for sim in range(2):
for iteration in range(2):
target_dir = f'ensemble_md/tests/data/stitch_test/sim_{sim}/iteration_{iteration}'
os.makedirs(target_dir)
shutil.copy(f'ensemble_md/tests/data/dhdl/dhdl_{sim * 2 + iteration}.xvg', f'{target_dir}/dhdl.xvg')
save_and_exclude(f'{target_dir}/dhdl.xvg', 40) # just keep the first 10 frames

# files = [[f'ensemble_md/tests/data/stitch_test/sim_{i}/iteration_{j}/dhdl_short.xvg' for j in range(2)] for i in range(2)] # noqa: E501
# shifts = [1, 1]

# More to come ...
# trajs_test = analyze_traj.stitch_time_series_for_sim(files, shifts, save=True)
# trajs_expected = [
# [0, 0, 3, 1, 4, 4, 5, 4, 5, 5, 4]
# ]

# Clean up
shutil.rmtree('ensemble_md/tests/data/stitch_test')


def test_stitch_xtc_trajs():
pass


Expand Down Expand Up @@ -115,32 +169,6 @@ def test_convert_npy2xvg():
os.chdir('../../../')


def test_stitch_time_series_for_sim():
# Set up files for testing
for sim in range(2):
for iteration in range(2):
target_dir = f'ensemble_md/tests/data/stitch_test/sim_{sim}/iteration_{iteration}'
os.makedirs(target_dir)
shutil.copy(f'ensemble_md/tests/data/dhdl/dhdl_{sim * 2 + iteration}.xvg', f'{target_dir}/dhdl.xvg')
save_and_exclude(f'{target_dir}/dhdl.xvg', 40) # just keep the first 10 frames

# files = [[f'ensemble_md/tests/data/stitch_test/sim_{i}/iteration_{j}/dhdl_short.xvg' for j in range(2)] for i in range(2)] # noqa: E501
# shifts = [1, 1]

# More to come ...
# trajs_test = analyze_traj.stitch_time_series_for_sim(files, shifts, save=True)
# trajs_expected = [
# [0, 0, 3, 1, 4, 4, 5, 4, 5, 5, 4]
# ]

# Clean up
shutil.rmtree('ensemble_md/tests/data/stitch_test')


def test_stitch_trajs():
pass


def test_traj2transmtx():
traj = [0, 1, 2, 1, 0, 3]
N = 4 # matrix size
Expand Down

0 comments on commit 8376dcd

Please sign in to comment.