diff --git a/Sources/python/SimEx/Analysis/XMDYNPhotonMatterAnalysis.py b/Sources/python/SimEx/Analysis/XMDYNPhotonMatterAnalysis.py index 784d7be9..a1d3f209 100644 --- a/Sources/python/SimEx/Analysis/XMDYNPhotonMatterAnalysis.py +++ b/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,14 +131,25 @@ 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 @@ -142,8 +157,11 @@ def load_snapshot(self, snapshot_index): 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,30 +171,66 @@ 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. """ @@ -184,30 +238,37 @@ def load_trajectory(self): 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) \ No newline at end of file diff --git a/Sources/python/SimEx/Parameters/DetectorGeometry.py b/Sources/python/SimEx/Parameters/DetectorGeometry.py index 629052b3..becfc986 100644 --- a/Sources/python/SimEx/Parameters/DetectorGeometry.py +++ b/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. diff --git a/Sources/python/SimEx/Utilities/IOUtilities.py b/Sources/python/SimEx/Utilities/IOUtilities.py index 52209815..e57a20de 100644 --- a/Sources/python/SimEx/Utilities/IOUtilities.py +++ b/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 \ No newline at end of file