minor fixes

Author: SimoneMartino98

10_onion_plots.py

@@ -1,4 +1,6 @@
-#%%
+# COMPUTE ONION
+# THIS CODE USES A DEPRECATED VERSION OF ONION
+# CONSIDER TO CHANGE THIS CODE WITH THE NEW VERSION
 import os
 import shutil
 from pathlib import Path

11_onion_results.py

@@ -1,4 +1,3 @@
-#%%
 import numpy as np 
 import matplotlib.pyplot as plt
 from matplotlib.ticker import MaxNLocator

1_compute_hdf5.py

@@ -1,9 +1,13 @@
-#%% Libraries
+#THIS CODE IS USED TO CONVERT .GRO AND .XTC TRAJECTORY INTO
+#HDF5 DATABASE
+#.GRO AND .XTC ZENODO LINK CAN BE FOUND IN THE ARTICLE
+
+#Libraries
 import MDAnalysis as mda
 import dynsight 
 import h5py
 import numpy as np
-#%% HDF5 file build
+# HDF5 file build
 print(f"{'-'*10}\nSTART INITILIAZIATION PART\n{'-'*10}")
 traj_name = "ice_water_O"
 simulation_folder = "simulation"
@@ -21,6 +25,7 @@
 u = mda.Universe(topo_file,traj_file)
 dynsight.hdf5er.mda_to_hdf5(u, out_file, "ice_water")
 in_file = out_file
+#O atoms will be used as reference
 with h5py.File(in_file,"r") as file:
     dataset_box = file["Trajectories"]["ice_water"]["Box"]
     dataset_traj = file["Trajectories"]["ice_water"]["Trajectory"]
@@ -36,10 +41,11 @@
     file["Trajectories"][traj_name].create_dataset("Box", data=box)
     file["Trajectories"][traj_name].create_dataset("Trajectory", data = Ox_traj, chunks=(100,len(Ox_types),3))
     file["Trajectories"][traj_name].create_dataset("Types", data = Ox_types)
-    
+    #Prepare for descriptors computation...
     file.create_group("LENS")
     file.create_group("SOAP")
 steps = Ox_traj.shape[0]
 atoms = Ox_traj.shape[1]
+#Just check
 print("\nReading HDF5 file:\n")
 print(f"Number of steps: {steps}\nNumber of particles: {atoms}\n")

2_compute_LENS.py

@@ -1,9 +1,10 @@
-#%%
+#THIS CODE IS USED TO COMPUTE LENS DESCRIPTOR
 import h5py
 import numpy as np
 import dynsight 
+#LENS cutoff
 LENS_CUTOFF = 10
-#%%
+#Load HDF5 file (see code 1_*)
 in_file = "ice_water_O.hdf5"
 traj_name = "ice_water_O"
 frames_range = slice(0,500)
@@ -23,7 +24,8 @@
 with h5py.File(in_file, "r") as file:
     LENS = np.array(file["LENS"][f"LENS_{int(LENS_CUTOFF)}"][0,:,:])
     np.save(f"arrays/LENS_{int(LENS_CUTOFF)}",LENS)
-# %% Spatial smoothing
+
+# Local denoising (Spatial smoothing)
 input_file = "ice_water_O.hdf5"
 with h5py.File(input_file, "r") as file:
     traj_array = np.array(file["Trajectories/ice_water_O/Trajectory"])
@@ -44,5 +46,3 @@
                                                     cutoff=cutoff, 
                                                     volume_shape = volume_shape)
     np.save(res_array,averaged.T)
-
-# %%

3_compute_SOAP.py

@@ -1,10 +1,10 @@
-#%%
+#THIS CODE IS USED TO COMPUTE SOAP DESCRIPTOR
 import h5py
 import dynsight
 import numpy as np
+#SOAP cutoff
 SOAP_CUTOFF = 10
 
-# %%
 traj_name = "ice_water_O"
 in_file = f"{traj_name}.hdf5"
 
@@ -26,10 +26,9 @@
 filled_soap = dynsight.soapify.fill_soap_vector_from_dscribe(soap,lmax=8,nmax=8)
 with h5py.File(in_file, "a") as file:
     file["SOAP"][f"SOAP_{int(SOAP_CUTOFF)}"].create_dataset("fill_SOAP", data=filled_soap)
-#np.save(f"arrays/SOAP_{int(SOAP_CUTOFF)}", soap)
 np.save(f"arrays/fullvect_SOAP_{int(SOAP_CUTOFF)}", filled_soap)
 
-# %% Spatial smoothing
+# Local denoising (Spatial smoothing)
 input_file = "ice_water_O.hdf5"
 with h5py.File(input_file, "r") as file:
     traj_array = np.array(file["Trajectories/ice_water_O/Trajectory"])

4_compute_PCA.py

@@ -1,11 +1,12 @@
-#%%
+#THIS CODE IS USED TO COMPUTE THE PCA ON SOAP
 import numpy as np
 from sklearn.decomposition import PCA
 import matplotlib.pyplot as plt
 import pandas as pd
+#SOAP cutoff
 SOAP_CUTOFF = 10
+#number of PCA components
 pc_components = 3
-#%%
 print(f"{'-'*10}\nPCA SOAP\n{'-'*10}")
 soap_av = np.load(f"arrays/fullvect_SOAP_{SOAP_CUTOFF}.npy")
 print(soap_av.shape)
@@ -14,9 +15,10 @@
 pca = PCA(n_components=pc_components)
 pc_soap = pca.fit_transform(SOAP)
 np.save(f"arrays/SOAP_{SOAP_CUTOFF}_PC1.npy",np.transpose(pc_soap[:,0].reshape(500,2048)))
+#To compute other components....
 #np.save(f"arrays/SOAP_{SOAP_CUTOFF}_PC2.npy",np.transpose(pc_soap[:,1].reshape(500,2048)))
 #np.save(f"arrays/SOAP_{SOAP_CUTOFF}_PC3.npy",np.transpose(pc_soap[:,2].reshape(500,2048)))
-# %%
+
 print(f"{'-'*10}\nPCA SP SOAP\n{'-'*10}")
 sp_cutoff = [10]
 for cutoff in sp_cutoff:
@@ -27,6 +29,7 @@
 
     pca = PCA(n_components=pc_components)
     pc_soap = pca.fit_transform(SOAP)
+    #To save other components....
     np.save(f"arrays/sp_{cutoff}_SOAP_{SOAP_CUTOFF}_PC1.npy",np.transpose(pc_soap[:,0].reshape(500,2048)))
     #np.save(f"arrays/sp_{cutoff}_SOAP_{SOAP_CUTOFF}_PC2.npy",np.transpose(pc_soap[:,1].reshape(500,2048)))
     #np.save(f"arrays/sp_{cutoff}_SOAP_{SOAP_CUTOFF}_PC3.npy",np.transpose(pc_soap[:,2].reshape(500,2048)))

5_compute_disp.py renamed to 5_compute_vel.py

@@ -1,8 +1,10 @@
-#%%
+#THIS CODE IS USED TO COMPUTE VELOCITIES
 import h5py
 import numpy as np
 import dynsight
-#%%
+
+#100 ps
+sampling = 100
 def read_from_xyz(filename):
     with open(filename, "r") as file:
         lines = file.readlines()
@@ -66,12 +68,12 @@ def compute_displacement(absolute,box):
     traj_array = traj_array.transpose(1,0,2)
     box_array = np.array(file[f"Trajectories/{traj_name}/Box"])
 trajectory = traj_array
-#disp_abs = compute_displacement(True)
+
 disp_rel = compute_displacement(False,box_array)
-disp_rel = disp_rel / 100 
+disp_rel = disp_rel / sampling 
 np.save("arrays/vel.npy", disp_rel)
-#np.save("arrays/disp_abs.npy", disp_abs)
-# %% Spatial smoothing
+
+# Local denoising (Spatial smoothing)
 input_file = "ice_water_O.hdf5"
 with h5py.File(input_file, "r") as file:
     traj_array = np.array(file["Trajectories/ice_water_O/Trajectory"])

6_compute_CN.py

@@ -1,12 +1,12 @@
-#%% 
+#THIS CODE IS USED TO COMPUTE Nneigh
 import h5py
 import MDAnalysis as mda
 import numpy as np
 from MDAnalysis.analysis.distances import distance_array
 import dynsight
 
+#Nneigh cutoff..
 nn_cutoff = 10
-#%%
 def compute_distance(p1, p2, Lx, Ly, Lz):
     dp = np.abs(p2-p1)
     dp = np.minimum(dp, np.array([Lx,Ly,Lz]) - dp)
@@ -16,7 +16,6 @@ def compute_distance(p1, p2, Lx, Ly, Lz):
 topo_file = f"{simulation_folder}/ice_water.gro"
 traj_file = f"{simulation_folder}/ice_water_500.xtc"
 u = mda.Universe(topo_file, traj_file)
-#selected_atoms = u.select_atoms("type 1 or type 2")
 print(u.dimensions)
 
 
@@ -32,7 +31,7 @@ def compute_distance(p1, p2, Lx, Ly, Lz):
 nn = nn.T
 np.save(f"arrays/nn_{nn_cutoff}.npy", nn[::3,:])
 
-# %% Spatial smoothing
+# Local denoising (Spatial smoothing)
 input_file = "ice_water_O.hdf5"
 with h5py.File(input_file, "r") as file:
     traj_array = np.array(file["Trajectories/ice_water_O/Trajectory"])

7_compute_dist5.py

@@ -1,11 +1,10 @@
-#%%
+#THIS CODE IS USED TO COMPUTE d5
 import h5py
 import MDAnalysis as mda
 import numpy as np
 from MDAnalysis.analysis.distances import distance_array
 import dynsight
 
-#%%
 print(f"{'-'*10}\nDIST 5\n{'-'*10}")
 simulation_folder = "simulation"
 topo_file = f"{simulation_folder}/ice_water.gro"
@@ -20,14 +19,12 @@
     id = np.argsort(distances)
     for i in range(0,6):
         print(f"atom {i}) coord: {selection.positions[i]}")
-    #print(f"distances: {sort_distances[0,1:5]}")
-    #print(f"id: {id[0,1:5]}")
     for i in range(sort_distances.shape[0]):
         dist_5[i,a] = sort_distances[i,5]
     a += 1
     dist_5 = np.array(dist_5)
 np.save("arrays/dist_5.npy",dist_5)
-# %% Spatial smoothing
+# Local denoising (Spatial smoothing)
 input_file = "ice_water_O.hdf5"
 with h5py.File(input_file, "r") as file:
     traj_array = np.array(file["Trajectories/ice_water_O/Trajectory"])

8_compute_OTO.py

@@ -1,8 +1,8 @@
-#%%
+#THIS CODE IS USED TO COMPUTE q_tet
 import numpy as np
 import h5py
 import dynsight
-#%%
+
 def read_from_xyz(filename):
     with open(filename, "r") as file:
         lines = file.readlines()
@@ -146,7 +146,7 @@ def compute_oto(filename, box):
 # write_xyz(trajectory,q,"oto.xyz", "Properties=pos:R:3:color:S:1")
 np.save("arrays/OTO.npy",q)
 
-# %% Spatial smoothing
+# Local denoising (Spatial smoothing)
 input_file = "ice_water_O.hdf5"
 with h5py.File(input_file, "r") as file:
     traj_array = np.array(file["Trajectories/ice_water_O/Trajectory"])

9_compute_checks.py

@@ -1,4 +1,4 @@
-#%%
+# JUST CHECK
 import numpy as np
 import os 
 
@@ -15,4 +15,3 @@
 for d in descriptors:
     arr = np.load(f"arrays/{d}")
     print(f"descriptor:     {d[:-4].ljust(20)} -      shape: {arr.shape}")
-# %%

D_plotter.ipynb

@@ -43,7 +43,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -62,11 +62,9 @@
     "import seaborn as sns\n",
     "import numpy as np\n",
     "\n",
-    "# Esempio di dati per pc_soap\n",
     "# pc_soap = np.random.rand(1000, 2) # Sostituire con i propri dati\n",
     "plt.scatter(pc_soap[::skip, 0], pc_soap[::skip, 1], s=0.5, color=\"black\", alpha=0.5)\n",
     "\n",
-    "# Aggiungere linee di densità\n",
     "sns.kdeplot(x=pc_soap[::skip, 0], y=pc_soap[::skip, 1], levels=10, color=\"red\", linewidths=1)\n",
     "\n",
     "plt.title('SOAP PCA components')\n",
@@ -171,7 +169,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -189,36 +187,28 @@
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "\n",
-    "# Esempio di due array con 3 valori ciascuno\n",
     "array1 = ev * 100\n",
     "array2 = ev_10 * 100\n",
     "\n",
-    "# Le etichette per le barre\n",
     "labels = ['PC1', 'PC2', 'PC3']\n",
     "\n",
-    "# Posizione delle barre sul grafico\n",
     "x = np.arange(len(labels))\n",
     "\n",
-    "# Larghezza delle barre\n",
     "width = 0.35\n",
     "\n",
     "fig, ax = plt.subplots(figsize=(10, 6))\n",
     "\n",
-    "# Tracciare le barre per ciascun array\n",
     "bars1 = ax.bar(x - width/2, array1, width, label='SOAP', color='blue')\n",
     "bars2 = ax.bar(x + width/2, array2, width, label='⟨SOAP⟩', color='orange')\n",
     "\n",
-    "# Aggiungere annotazioni sopra le barre per il primo array\n",
     "for bar in bars1:\n",
     "    yval = bar.get_height()\n",
     "    ax.text(bar.get_x() + bar.get_width() / 2.0, yval, f'{yval:.2f}', va='bottom', ha='center', fontweight='bold',fontsize=fontsize)\n",
     "\n",
-    "# Aggiungere annotazioni sopra le barre per il secondo array\n",
     "for bar in bars2:\n",
     "    yval = bar.get_height()\n",
     "    ax.text(bar.get_x() + bar.get_width() / 2.0, yval, f'{yval:.2f}', va='bottom', ha='center', fontweight='bold',fontsize=fontsize)\n",
     "\n",
-    "# Aggiungere etichette, titolo e legenda\n",
     "ax.set_xlabel('PCA Components', fontweight='bold',fontsize=fontsize)\n",
     "ax.set_ylabel('Explained variance (%)', fontweight='bold',fontsize=fontsize)\n",
     "ax.set_title('Explained variance', fontsize=fontsize+2)\n",