Bugfixes for XMDYNPhotonMatterAnalysis

Sources/python/SimEx/Analysis/XMDYNPhotonMatterAnalysis.py

@@ -1,15 +1,15 @@
 #!/usr/bin/env python
 """:module XMDYNPhotonMatterAnalysis: Hosting utilities to analyse and visualize photon-matter trajectories generated by XMDYN."""
 
+from random import sample
 import h5py
 import numpy
 import os
-import pylab
 import periodictable as pte
 
 ELEMENT_SYMBOL = ['All'] + [e.symbol for e in pte.elements]
 from SimEx.Analysis.AbstractAnalysis import AbstractAnalysis, plt
-from SimEx.Utilities.IOUtilities import loadPDB
+from SimEx.Utilities.IOUtilities import loadPDB, loadXYZ
 
 
 class XMDYNPhotonMatterAnalysis(AbstractAnalysis):
@@ -38,6 +38,10 @@ def __init__( self, input_path=None, snapshot_indices=None, elements=None, sampl
 
         self.load_trajectory()
 
+    @property
+    def t(self):
+        return self.__trajectory["time"]
+
     @property
     def input_path(self):
         """ Query the input path. """
@@ -127,23 +131,37 @@ def elements(self, value):
 
         self.__elements = value
 
+    def get_element_symbols(self):
+        """return the symbol of the amalyzed elements as a list of str"""
+        sample = self.get_sample()
+        z = sample['selZ'].keys()
+        if self.elements != ["All"]:
+            z = [a for a in z if a in self.elements]
+        return [pte.elements[a] for a in z]
+
+
     def load_snapshot(self, snapshot_index):
         """ Load snapshot data from hdf5 file into memory. """
 
         snp = snapshot_index
         dbase_root = "/data/snp_" + str( snp ).zfill(self.__num_digits) + "/"
+        time_path = "/misc/time/snp_" + str( snp ).zfill(self.__num_digits) + "/"
         xsnp = dict()
 
         with h5py.File(self.input_path) as fp:
+            xsnp['t'] = fp.get(time_path).value.item()
             xsnp['Z'] = fp.get(dbase_root+ 'Z').value
             xsnp['T'] = fp.get(dbase_root + 'T').value
             xsnp['ff'] = fp.get(dbase_root + 'ff').value
             xsnp['xyz'] = fp.get(dbase_root + 'xyz').value
             xsnp['r'] = fp.get(dbase_root + 'r').value
             xsnp['Nph'] = fp.get(dbase_root + 'Nph').value
             N = xsnp['Z'].size
-            xsnp['q'] = numpy.array([xsnp['ff'][pylab.find(xsnp['T']==x), 0] for x in xsnp['xyz']]).reshape(N,)
-            xsnp['snp'] = snp ;
+            xsnp['q'] = numpy.array([xsnp['ff'][numpy.nonzero(numpy.ravel(xsnp['T']==x))[0], 0] for x in xsnp['xyz']]).reshape(N,)
+            # xsnp['q'] = numpy.array([xsnp['ff'][(numpy.array(xsnp['T'])==numpy.array(x)), 0] for x in xsnp['xyz']]).reshape(N,)
+            # xsnp['q'] = numpy.array([xsnp['ff'][(numpy.equals(xsnp['T'],x), 0] for x in xsnp['xyz']]).reshape(N,)
+            #xsnp['q'] = numpy.array([numpy.array(xsnp['ff'][xsnp['T']==x][0] for x in xsnp['xyz']]).reshape(N,)
+            xsnp['snp'] = snp
 
         return xsnp
 
@@ -153,61 +171,104 @@ def number_of_snapshots(self) :
             count = len([k for k in xfp['data'].keys() if "snp_" in k])
         return count
 
+    def get_avg_displacement(self):
+        """ Get the average displacement per atomic species as function of time. """
+        return self.__trajectory['displacement'].T
+
+    def get_max_displacement(self):
+        """ Get the average displacement per atomic species as function of time. """
+        return self.__trajectory['max_displacement'].T
+
     def plot_displacement(self):
         """ Plot the average displacement per atomic species as function of time. """
 
-        #t = self.trajectory['time']
-
-        for d in self.__trajectory['displacement'].T:
-            plt.plot(d) # self.trajectory['disp'][ : , pylab.find( sel_Z == pylab.array( list(data['sample']['selZ'].keys()) ) ) ] , xcolor  )
+        for d, sym in zip(self.get_avg_displacement(), self.get_element_symbols()):
+            plt.plot(self.t*1e15, d, label=sym) # self.trajectory['disp'][ : , pylab.find( sel_Z == pylab.array( list(data['sample']['selZ'].keys()) ) ) ] , xcolor  )
         plt.xlabel( 'Time [fs]' )
         plt.ylabel( 'Average displacement [$\AA$]' )
+        plt.legend()
+
+    def get_sample(self):
+        """Load the sample and return it."""
+        try:
+            sample = load_sample(self.sample_path)
+        except:
+            if os.path.splitext(self.sample_path)[1] in ['.xy', '.xyz']:
+                try:
+                    sample = loadXYZ(self.sample_path)
+                    print("Using experimental XYZ loader")
+                except:
+                    sample = loadPDB(self.sample_path)
+            else:
+                sample = loadPDB(self.sample_path)
+
+        return sample
+
+    def get_avg_charge(self):
+        """ Plot the average number of electrons per atom per atomic species as function of time. """
+        return self.__trajectory['charge'].T
 
     def plot_charge(self):
         """ Plot the average number of electrons per atom per atomic species as function of time. """
 
-        for d in self.__trajectory['charge'].T:
-            plt.plot(d)
+        # for d in self.__trajectory['charge'].T:
+        for d, sym in zip(self.get_avg_charge(), self.get_element_symbols()):
+            plt.plot(self.t*1e15, d, label=sym)
 
         ### TODO: time axis, labels, legend
         plt.xlabel( 'Time [fs]' )
         plt.ylabel( 'Number of bound electrons per atom' )
+        plt.legend()
 
     def plot_energies(self):
         """ Plot the evolution of MD energies over the simulation time. """
         raise RuntimeError("Not implemented yet.")
 
+    def get_angle(self):
+        """ Get the angle quaternion of the sample. """
+
+        with h5py.File(self.input_path) as fp:
+            qt = fp.get("/params/angle/").value
+        return qt
+
     def load_trajectory(self):
         """ Load the selected snapshots and extract data to analyze. """
 
         trajectory = dict()
 
         sample = None
         disp = []
+        max_disp = []
         charge = []
         time = []
+
+        sample = self.get_sample()
+
         # Read sample data.
-        try:
-            sample = load_sample(self.sample_path)
-        except:
-            sample = loadPDB(self.sample_path)
+        # try:
+        #     sample = load_sample(self.sample_path)
+        # except:
+        #    sample = loadPDB(self.sample_path)
 
         snapshot_indices = self.snapshot_indices
-        if snapshot_indices == "All":
+        if snapshot_indices == "All" or snapshot_indices == ["All"]:
             snapshot_indices = range(1, self.number_of_snapshots() + 1)
 
         for si in snapshot_indices:
             snapshot = self.load_snapshot(si)
-
+            time.append(snapshot['t'])
             disp.append(calculate_displacement(snapshot, r0=sample['r'], sample=sample))
+            max_disp.append(calculate_max_displacement(snapshot, r0=sample['r'], sample=sample))
             charge.append(calculate_ion_charge(snapshot, sample))
 
         trajectory['displacement'] = numpy.array(disp)
+        trajectory['max_displacement'] = numpy.array(max_disp)
         trajectory['charge'] = numpy.array(charge)
         trajectory['time'] = numpy.array(time)
 
         self.__trajectory = trajectory
 
+
     def animate(self):
         """ Generate an animation of atom trajectories and their ionization. """
         pass
@@ -224,7 +285,7 @@ def load_sample(sample_path) :
     sample['selZ'] = dict()
 
     for sel_Z in numpy.unique(sample['Z']) :
-        sample['selZ'][sel_Z] = pylab.find(sel_Z == sample['Z'])
+        sample['selZ'][sel_Z] = (numpy.array(sel_Z) == numpy.array(sample['Z']))
 
     return sample
 
@@ -234,6 +295,7 @@ def read_h5_dataset( path , dataset ) :
         data = xfp.get( dataset ).value
     return data
 
+
 def calculate_displacement(snapshot, r0, sample) :
     """ Calculate the average displacement per atomic species in a snapshot.
 
@@ -261,6 +323,35 @@ def calculate_displacement(snapshot, r0, sample) :
 
     return numpy.array(all_disp)
 
+
+def calculate_max_displacement(snapshot, r0, sample) :
+    """ Find the highest displacement value per atomic species in a snapshot.
+
+    :param snapshot: The snapshot to analyze
+    :type snapshot: dict
+
+    :param r0: Unperturbed positions of the sample atoms.
+    :type r0: numpy.array (shape=(Natoms, 3))
+    ### CHECKME: Can't we read r0 from the sample dict?
+
+    :param sample: Sample data
+    :type sample: dict
+
+    """
+
+    num_Z = len( list(sample['selZ'].keys()) )
+    all_disp = numpy.zeros( ( num_Z , ) )
+
+    count = 0
+
+    for sel_Z in list(sample['selZ'].keys()) :
+        dr = snapshot['r'][sample['selZ'][sel_Z],:] - r0[sample['selZ'][sel_Z],:]
+        all_disp[count] = numpy.amax(numpy.sqrt(numpy.sum(dr*dr, axis=1))) / 1e-10
+        count = count + 1
+
+    return numpy.array(all_disp)
+
+
 def calculate_ion_charge(snapshot, sample):
     """ Calculate the remaining electric charge per atomic species of a given snapshot.
 
@@ -280,4 +371,4 @@ def calculate_ion_charge(snapshot, sample):
         all_numE[count] = numpy.mean( snapshot['q'][sample['selZ'][sel_Z]] )
         count = count + 1
 
-    return numpy.array(all_numE)
+    return numpy.array(all_numE)

Sources/python/SimEx/Parameters/DetectorGeometry.py

@@ -53,7 +53,7 @@ def __init__(self,
         """
         :param ranges: The minimum and maximum values pixel numbers on the respective transverse axis.
         :type  ranges: Dictionary
-        ":example ranges: {"fast_scan_min : 11, "fast_scan_max" : 20, "slow_scan_min" : 1, "fast_scan_max" : 20} # First axis from 11 to 20 and second axis from 1 to 20."
+        ":example ranges: {"fast_scan_min : 11, "fast_scan_max" : 20, "slow_scan_min" : 1, "slow_scan_max" : 20} # First axis from 11 to 20 and second axis from 1 to 20."
 
         :param pixel_size: The physical size of the pixel (assuming quadratic shape) (SI units).
         :type  pixel_size: PhysicalQuantity with unit meter.

Sources/python/SimEx/Utilities/IOUtilities.py

@@ -31,6 +31,9 @@
 import requests
 import uuid
 
+import warnings
+warnings.simplefilter("default")
+
 def getTmpFileName():
     """ Create a unique filename
     :return: unique filename for temporary storage
@@ -400,6 +403,4 @@ def get_dict_from_lines(reader):
         ret[key] = val
 
     # Return finished dict.
-    return ret
-
-
+    return ret