Skip to content

Commit

Permalink
rescructured data
Browse files Browse the repository at this point in the history
  • Loading branch information
@simongravelle
simongravelle committed Jan 4, 2024
1 parent 1517632 commit c7195dd
Show file tree
Hide file tree
Showing 20 changed files with 170 additions and 0 deletions.
170 changes: 170 additions & 0 deletions examples/illustrations/bulk-water-high-dumping-date/analysis.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,170 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import MDAnalysis as mda\n",
"import nmrformd"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"import sys, os, git\n",
"current_path = os.getcwd()\n",
"git_repo = git.Repo(current_path, search_parent_directories=True)\n",
"git_path = git_repo.git.rev_parse(\"--show-toplevel\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def save_result(data, name, aniso = False):\n",
" \"\"\"Save the correlation functions in dictionary\"\"\"\n",
" if not os.path.exists(\"raw_data/\"):\n",
" os.makedirs(\"raw_data/\")\n",
" saving_file = \"raw_data/\" + name + \".npy\"\n",
" t = data.t\n",
" f = data.f\n",
" if aniso:\n",
" C1 = data.gij[0]\n",
" C2 = data.gij[1]\n",
" C3 = data.gij[2]\n",
" else:\n",
" C = data.gij[0]\n",
" R1 = data.R1\n",
" R2 = data.R2\n",
" N = data.group_j.atoms.n_atoms\n",
" try:\n",
" previous_dictionary = np.load(saving_file, allow_pickle=True)\n",
" t_prev = np.real(previous_dictionary.item()[\"t\"])\n",
" assert len(t_prev) == len(t)\n",
" if aniso:\n",
" C1_prev = np.real(previous_dictionary.item()[\"C1\"])\n",
" C2_prev = np.real(previous_dictionary.item()[\"C2\"])\n",
" C3_prev = np.real(previous_dictionary.item()[\"C3\"])\n",
" else:\n",
" C_prev = np.real(previous_dictionary.item()[\"C\"])\n",
" R1_prev = np.real(previous_dictionary.item()[\"R1\"])\n",
" R2_prev = np.real(previous_dictionary.item()[\"R2\"])\n",
" N_prev = np.real(previous_dictionary.item()[\"N\"])\n",
" if aniso:\n",
" C1 = (C1*N + C1_prev*N_prev) / (N_prev + N)\n",
" C2 = (C2*N + C2_prev*N_prev) / (N_prev + N)\n",
" C3 = (C3*N + C3_prev*N_prev) / (N_prev + N)\n",
" else:\n",
" C = (C*N + C_prev*N_prev) / (N_prev + N)\n",
" R1 = (R1*N + R1_prev*N_prev) / (N_prev + N)\n",
" R2 = (R2*N + R2_prev*N_prev) / (N_prev + N)\n",
" N += N_prev\n",
" except:\n",
" pass\n",
" if aniso:\n",
" dictionary = {\n",
" \"t\": t,\n",
" \"f\": f,\n",
" \"C1\": C1,\n",
" \"C2\": C2,\n",
" \"C3\": C3,\n",
" \"N\": N,\n",
" \"R1\": R1,\n",
" \"R2\": R2,\n",
" }\n",
" else:\n",
" dictionary = {\n",
" \"t\": t,\n",
" \"f\": f,\n",
" \"C\": C,\n",
" \"N\": N,\n",
" \"R1\": R1,\n",
" \"R2\": R2,\n",
" }\n",
" np.save(saving_file, dictionary)\n",
" return N"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8000 1 137 0.2\n",
"8000 1 138 0.2\n"
]
},
{
"ename": "KeyboardInterrupt",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[8], line 17\u001b[0m\n\u001b[1;32m 15\u001b[0m n_hydrogen \u001b[38;5;241m=\u001b[39m hydrogen\u001b[38;5;241m.\u001b[39mn_atoms\n\u001b[1;32m 16\u001b[0m intra \u001b[38;5;241m=\u001b[39m nmrformd\u001b[38;5;241m.\u001b[39mNMR(u, atom_group \u001b[38;5;241m=\u001b[39m hydrogen, neighbor_group \u001b[38;5;241m=\u001b[39m hydrogen, type_analysis \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mintra_molecular\u001b[39m\u001b[38;5;124m\"\u001b[39m, number_i\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m5\u001b[39m)\n\u001b[0;32m---> 17\u001b[0m inter \u001b[38;5;241m=\u001b[39m \u001b[43mnmrformd\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mNMR\u001b[49m\u001b[43m(\u001b[49m\u001b[43mu\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43matom_group\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mhydrogen\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mneighbor_group\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mhydrogen\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtype_analysis\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43minter_molecular\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnumber_i\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 18\u001b[0m n_intra \u001b[38;5;241m=\u001b[39m save_result(intra, name \u001b[38;5;241m=\u001b[39m N \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m_intra_dt\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mstr\u001b[39m(dt))\n\u001b[1;32m 19\u001b[0m n_inter \u001b[38;5;241m=\u001b[39m save_result(inter, name \u001b[38;5;241m=\u001b[39m N \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m_inter_dt\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mstr\u001b[39m(dt))\n",
"File \u001b[0;32m~/.local/lib/python3.11/site-packages/nmrformd/NMR.py:100\u001b[0m, in \u001b[0;36mNMR.__init__\u001b[0;34m(self, u, atom_group, neighbor_group, type_analysis, number_i, isotropic, f0, actual_dt, hydrogen_per_atom, spin, pbc)\u001b[0m\n\u001b[1;32m 97\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mT2 \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[1;32m 99\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39minitialize()\n\u001b[0;32m--> 100\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcollect_data\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 101\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfinalize()\n",
"File \u001b[0;32m~/.local/lib/python3.11/site-packages/nmrformd/NMR.py:161\u001b[0m, in \u001b[0;36mNMR.collect_data\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 159\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m cpt_i \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[1;32m 160\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39minitialise_data()\n\u001b[0;32m--> 161\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mloop_over_trajectory\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 162\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcalculate_correlation_ij()\n",
"File \u001b[0;32m~/.local/lib/python3.11/site-packages/nmrformd/NMR.py:236\u001b[0m, in \u001b[0;36mNMR.loop_over_trajectory\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 232\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m np\u001b[38;5;241m.\u001b[39mall(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbox[\u001b[38;5;241m3\u001b[39m:] \u001b[38;5;241m==\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbox[\u001b[38;5;241m3\u001b[39m:][\u001b[38;5;241m0\u001b[39m]) \u001b[38;5;241m&\u001b[39m np\u001b[38;5;241m.\u001b[39mall(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbox[\u001b[38;5;241m3\u001b[39m:][\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m90.0\u001b[39m):\n\u001b[1;32m 233\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNMRforMD does not accept non-orthogonal box\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 234\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mYou can use triclinic_to_orthorhombic from the package\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 235\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mlipyphilic to convert the trajectory file.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m--> 236\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvector_rij\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 237\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcartesian_to_spherical()\n\u001b[1;32m 238\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mevaluate_function_F()\n",
"File \u001b[0;32m~/.local/lib/python3.11/site-packages/nmrformd/NMR.py:248\u001b[0m, in \u001b[0;36mNMR.vector_rij\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 243\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Calculate distance between position_i and position_j.\u001b[39;00m\n\u001b[1;32m 244\u001b[0m \u001b[38;5;124;03m\u001b[39;00m\n\u001b[1;32m 245\u001b[0m \u001b[38;5;124;03mBy defaults, periodic boundary conditions are assumed, but can also be turned off.\u001b[39;00m\n\u001b[1;32m 246\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 247\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpbc:\n\u001b[0;32m--> 248\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mrij \u001b[38;5;241m=\u001b[39m (\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mremainder\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mposition_i\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mposition_j\u001b[49m\n\u001b[1;32m 249\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbox\u001b[49m\u001b[43m[\u001b[49m\u001b[43m:\u001b[49m\u001b[38;5;241;43m3\u001b[39;49m\u001b[43m]\u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[38;5;241;43m2.\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbox\u001b[49m\u001b[43m[\u001b[49m\u001b[43m:\u001b[49m\u001b[38;5;241;43m3\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;241m-\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbox[:\u001b[38;5;241m3\u001b[39m]\u001b[38;5;241m/\u001b[39m\u001b[38;5;241m2.\u001b[39m)\u001b[38;5;241m.\u001b[39mT\n\u001b[1;32m 250\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 251\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mrij \u001b[38;5;241m=\u001b[39m (\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mposition_i \u001b[38;5;241m-\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mposition_j)\u001b[38;5;241m.\u001b[39mT\n",
"\u001b[0;31mKeyboardInterrupt\u001b[0m: "
]
}
],
"source": [
"while 1<2: \n",
" # high res\n",
" N = \"N4000_HR\"\n",
" for mystep in [1, 2, 4, 8, 16, 32]:\n",
" datapath = git_path + \"/data/bulk-water-high-dumping-rate/bulk-water/raw-data/N4000_HR/\"\n",
" xtc = [datapath+\"prod1.xtc\", datapath+\"prod2.xtc\",\n",
" datapath+\"prod3.xtc\", datapath+\"prod4.xtc\",\n",
" datapath+\"prod5.xtc\", datapath+\"prod6.xtc\",\n",
" datapath+\"prod7.xtc\", datapath+\"prod8.xtc\",\n",
" datapath+\"prod9.xtc\", datapath+\"prod10.xtc\"]\n",
" u = mda.Universe(datapath+\"prod1.tpr\", xtc)\n",
" u.transfer_to_memory(step = mystep)\n",
" dt = np.round(u.trajectory.dt,2)\n",
" hydrogen = u.select_atoms(\"name HW1 HW2\")\n",
" n_hydrogen = hydrogen.n_atoms\n",
" intra = nmrformd.NMR(u, atom_group = hydrogen, neighbor_group = hydrogen, type_analysis = \"intra_molecular\", number_i=5)\n",
" inter = nmrformd.NMR(u, atom_group = hydrogen, neighbor_group = hydrogen, type_analysis = \"inter_molecular\", number_i=1)\n",
" n_intra = save_result(intra, name = \"N4000_HR_intra_dt\" + str(dt))\n",
" n_inter = save_result(inter, name = \"N4000_HR_inter_dt\" + str(dt))\n",
" print(n_hydrogen, mystep, n_intra, dt)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Binary file added examples/illustrations/bulk-water-high-dumping-date/raw_data/N4000_HR_inter_dt0.2.npy
View file
Binary file not shown.
Binary file added examples/illustrations/bulk-water-high-dumping-date/raw_data/N4000_HR_inter_dt0.4.npy
View file
Binary file not shown.
Binary file added examples/illustrations/bulk-water-high-dumping-date/raw_data/N4000_HR_inter_dt0.8.npy
View file
Binary file not shown.
Binary file added examples/illustrations/bulk-water-high-dumping-date/raw_data/N4000_HR_inter_dt1.6.npy
View file
Binary file not shown.
Binary file added examples/illustrations/bulk-water-high-dumping-date/raw_data/N4000_HR_inter_dt3.2.npy
View file
Binary file not shown.
Binary file added examples/illustrations/bulk-water-high-dumping-date/raw_data/N4000_HR_inter_dt6.4.npy
View file
Binary file not shown.
Binary file renamed ...k-water/raw_data/N4000_HR_intra_dt0.2.npy → ...ng-date/raw_data/N4000_HR_intra_dt0.2.npy
View file
Binary file not shown.
Binary file renamed ...k-water/raw_data/N4000_HR_intra_dt0.4.npy → ...ng-date/raw_data/N4000_HR_intra_dt0.4.npy
View file
Binary file not shown.
Binary file renamed ...k-water/raw_data/N4000_HR_intra_dt0.8.npy → ...ng-date/raw_data/N4000_HR_intra_dt0.8.npy
View file
Binary file not shown.
Binary file renamed ...k-water/raw_data/N4000_HR_intra_dt1.6.npy → ...ng-date/raw_data/N4000_HR_intra_dt1.6.npy
View file
Binary file not shown.
Binary file renamed ...k-water/raw_data/N4000_HR_intra_dt3.2.npy → ...ng-date/raw_data/N4000_HR_intra_dt3.2.npy
View file
Binary file not shown.
Binary file renamed ...k-water/raw_data/N4000_HR_intra_dt6.4.npy → ...ng-date/raw_data/N4000_HR_intra_dt6.4.npy
View file
Binary file not shown.
Binary file removed examples/illustrations/bulk-water/raw_data/N4000_HR_inter_dt0.2.npy
View file
Binary file not shown.
Binary file removed examples/illustrations/bulk-water/raw_data/N4000_HR_inter_dt0.4.npy
View file
Binary file not shown.
Binary file removed examples/illustrations/bulk-water/raw_data/N4000_HR_inter_dt0.8.npy
View file
Binary file not shown.
Binary file removed examples/illustrations/bulk-water/raw_data/N4000_HR_inter_dt1.6.npy
View file
Binary file not shown.
Binary file removed examples/illustrations/bulk-water/raw_data/N4000_HR_inter_dt3.2.npy
View file
Binary file not shown.
Binary file removed examples/illustrations/bulk-water/raw_data/N4000_HR_inter_dt6.4.npy
View file
Binary file not shown.
0 tests/test_bulk_water.py → tests/old_test
View file
File renamed without changes.

0 comments on commit c7195dd

Please sign in to comment.