From 597e6778b2c3525889d30976fd9a5d44b6e3da5a Mon Sep 17 00:00:00 2001 From: SanderBorgmans Date: Thu, 3 Aug 2023 14:29:38 +0200 Subject: [PATCH] Added RED-COFFEE notebook as an example --- README.md | 3 +- mkdocs.yml | 3 - notebooks/PXRD.ipynb | 1023 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 1024 insertions(+), 5 deletions(-) create mode 100644 notebooks/PXRD.ipynb diff --git a/README.md b/README.md index 49e505f..7b3ce82 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,6 @@ ![gpxrdpy](./docs/gpxrdpy_banner_light.svg#gh-light-mode-only) ![gpxrdpy](./docs/gpxrdpy_banner_dark.svg#gh-dark-mode-only) -# gpxrdpy Python wrapper for PXRD pattern calculation based on pyobjcryst. There are four main usages: * calculate PXRD pattern based on cif file - `calc` * compare two PXRD patterns based on .tsv files - `comp` @@ -20,7 +19,7 @@ Usage: The default values for the `run_up_time` and `exp_filename` should be altered within this bash script before execution. -A version of gpxrdpy is also implemented in pyiron at https://github.com/SanderBorgmans/pyiron/tree/hpc_ugent_2020. +A version of gpxrdpy is also implemented in pyiron at https://github.com/SanderBorgmans/pyiron/tree/hpc_ugent_2020. An example notebook, created as supporting information for [this paper](https://doi.org/10.1039/D3TA00470H), can be found in the notebooks folder. In this notebook, static and dynamic simulations are performed, and provided as input to the gpxrdpy module of pyiron. ## Requirements * `pyobjcryst` - Python bindings to ObjCryst++, the Object-Oriented Crystallographic Library (see https://github.com/diffpy/pyobjcryst) diff --git a/mkdocs.yml b/mkdocs.yml index 5aa242d..0e73195 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -17,9 +17,6 @@ nav: - Introduction: index.md - Installation: installation.md -hide: -- toc - repo_url: https://github.com/SanderBorgmans/gpxrdpy markdown_extensions: - md_in_html diff --git a/notebooks/PXRD.ipynb b/notebooks/PXRD.ipynb new file mode 100644 index 0000000..88ecfab --- /dev/null +++ b/notebooks/PXRD.ipynb @@ -0,0 +1,1023 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Import statements\n", + "import numpy as np\n", + "import matplotlib.pyplot as pt\n", + "import glob\n", + "import os\n", + "\n", + "from ase.io import read\n", + "from pyiron import Project, ase_to_pyiron\n", + "\n", + "from molmod.units import *\n", + "from molmod.constants import *\n", + "\n", + "from collections import namedtuple\n", + "from dataclasses import dataclass # replaces namedtuple with mutable attributes\n", + "%matplotlib inline\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Gather structures \n", + "\n", + "'''@dataclass\n", + "class s_object:\n", + " structure: object\n", + " ffatypes: object\n", + " ffatype_ids: object\n", + " bonds : object'''\n", + " \n", + " \n", + "class Sobject(object):\n", + " def __init__(self,structure,ffatypes,ffatype_ids,bonds):\n", + " self.structure = structure\n", + " self.ffatypes = ffatypes\n", + " self.ffatype_ids = ffatype_ids\n", + " self.bonds = bonds\n", + " \n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Execution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "structures_database = {}\n", + "\n", + "for block in glob.glob('./input_files/PXRD/Structures_Database/*/*.chk'):\n", + " block_name = block.split('/')[-2]\n", + " print(block_name)\n", + " \n", + " tmp = pr.create_job(pr.job_type.Yaff,'tmp',delete_existing_job=True)\n", + " tmp.load_chk(block)\n", + " structures_database[block_name] = Sobject(tmp.structure,tmp.ffatypes,tmp.ffatype_ids,tmp.bonds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pr_database = Project('PXRD/Juul_database')\n", + "pr_database_uff = Project('PXRD/Juul_database_uff')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Optimizations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# yaff job function\n", + "def yaff_opt_job(pr, name, sobject, ffpars):\n", + " job = pr.create_job(pr.job_type.Yaff, name, delete_existing_job=True)\n", + " \n", + " job.calc_minimize(max_iter=10000, cell=True)\n", + "\n", + " job.set_ffpars(fnames=ffpars)\n", + " job.structure = sobject.structure\n", + " job.ffatypes = sobject.ffatypes\n", + " job.ffatype_ids = sobject.ffatype_ids\n", + " job.bonds = sobject.bonds\n", + " \n", + " job.input['rcut'] = 11.0*angstrom\n", + " job.input['alpha_scale'] = 2.86\n", + " job.input['gcut_scale'] = 1.0\n", + " job.input['tailcorrections'] = True\n", + "\n", + " job.executable.version = '2020'\n", + " job.server.queue = 'slaking'\n", + " job.server.cores = 1\n", + " job.server.run_time = 5*60*60 # in seconds\n", + "\n", + " job.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for bn,s in structures_database.items():\n", + " ffpars = os.path.join('./input_files/PXRD/Structures_Database', bn, 'pars_cluster.txt')\n", + " bn = bn.replace('-','_')\n", + " yaff_opt_job(pr_database, bn, s, ffpars)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for bn,s in structures_database.items():\n", + " ffpars = os.path.join('./input_files/PXRD/Structures_Database', bn, 'pars_uff.txt')\n", + " bn = bn.replace('-','_')\n", + " yaff_opt_job(pr_database_uff, bn, s, ffpars)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Molecular Dynamics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pr_database_md = Project('PXRD/Juul_database_md')\n", + "pr_database_uff_md = Project('PXRD/Juul_database_uff_md')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def md_job_v2(pr,name,sobject,structure,ffpars,temp=300*kelvin,press=1e5*pascal,nsteps=400000, \n", + " time_step=0.5*femtosecond,n_print=200,timecon_thermo=100.0*femtosecond,timecon_baro=1000.0*femtosecond, repeat=(1,1,1), log_lammps=False):\n", + " # MD code\n", + " job = pr.create_job(pr.job_type.Yaff, name, delete_existing_job=True)\n", + " job.calc_md(temperature=temp, pressure=press, nsteps=nsteps, time_step=time_step, n_print=n_print,\n", + " timecon_thermo=timecon_thermo, timecon_baro=timecon_baro)\n", + " \n", + " job.set_ffpars(ffpars)\n", + " \n", + " job.input['rcut'] = 11.0*angstrom\n", + " job.input['alpha_scale'] = 2.86\n", + " job.input['gcut_scale'] = 1.0\n", + " job.input['tailcorrections'] = True\n", + " \n", + " if log_lammps:\n", + " job.input['log_lammps'] = 'lammps.log'\n", + " \n", + " # Load structure and atomtypes\n", + " job.structure = structure.repeat(repeat)\n", + " job.ffatypes = sobject.ffatypes\n", + " job.ffatype_ids = np.tile(sobject.ffatype_ids,np.prod(repeat))\n", + " \n", + " # bonds will automatically be detected in the supercell\n", + " # assume that optimization has taken care of structural irregularities\n", + " \n", + " job.enable_lammps(executable='2020_lammps')\n", + " \n", + " job.server.queue = 'slaking'\n", + " job.server.cores = 8 \n", + " job.server.run_time = 72*60*60 # in seconds\n", + "\n", + " job.run()\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for bn,s in structures_database.items():\n", + " ffpars = os.path.join('./input_files/PXRD/Structures_Database', bn, 'pars_cluster.txt')\n", + " bn = bn.replace('-','_')\n", + " opt_structure = pr_database.load(bn).get_structure()\n", + "\n", + " # if the structure is 2D make a 1x1x5 supercell, if it is very small, use a 2x2x5 supercell\n", + " if any(bn.startswith(k) for k in ['hcb','kgm','sql']):\n", + " if len(opt_structure)<280: # 140 atoms per layer\n", + " repeat = (2,2,5)\n", + " else:\n", + " repeat = (1,1,5)\n", + " else:\n", + " repeat = (1,1,1)\n", + "\n", + " md_job_v2(pr_database_md,bn,s,opt_structure,ffpars,repeat=repeat)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for bn,s in structures_database.items():\n", + " ffpars = os.path.join('./input_files/PXRD/Structures_Database', bn, 'pars_uff.txt')\n", + " bn = bn.replace('-','_')\n", + " opt_structure = pr_database_uff.load(bn).get_structure()\n", + "\n", + " # if the structure is 2D make a 1x1x5 supercell, if it is very small, use a 2x2x5 supercell\n", + " if any(bn.startswith(k) for k in ['hcb','kgm','sql']):\n", + " if len(opt_structure)<280: # 140 atoms per layer\n", + " repeat = (2,2,5)\n", + " else:\n", + " repeat = (1,1,5)\n", + " else:\n", + " repeat = (1,1,1)\n", + "\n", + " md_job_v2(pr_database_uff_md,bn,s,opt_structure,ffpars,repeat=repeat)\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### PXRD" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Static" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pr_database_pxrd = pr_database.create_group('PXRD')\n", + "pr_database_uff_pxrd = pr_database_uff.create_group('PXRD')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def static_PXRD_job(pr,name,structure,refpattern=None):\n", + " job = pr.create_job(pr.job_type.GPXRD, name, delete_existing_job=True)\n", + " job.structure = structure\n", + " job.input['jobtype'] = 'static'\n", + " \n", + " if refpattern is not None:\n", + " job.set_reference_pattern(refpattern)\n", + "\n", + " job.server.queue = 'slaking'\n", + " job.server.cores = 1 \n", + " job.server.run_time = 15*60 # in seconds\n", + "\n", + " job.run()\n", + " #return job" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for bn in structures_database.keys():\n", + " name = bn.replace('-','_')\n", + " refpattern = glob.glob('./input_files/PXRD/Structures_Database/'+bn+'/*.dat')[0]\n", + " job = pr_database.load(name)\n", + " job_uff = pr_database_uff.load(name)\n", + "\n", + " static_PXRD_job(pr_database_pxrd,name,job.get_structure(),refpattern=refpattern)\n", + " static_PXRD_job(pr_database_uff_pxrd,name,job_uff.get_structure(),refpattern=refpattern)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.ticker import MaxNLocator\n", + "\n", + "for bn in structures_database.keys():\n", + " name = bn.replace('-','_')\n", + " \n", + " job = pr_database_pxrd.load(name)\n", + " job_uff = pr_database_uff_pxrd.load(name)\n", + " \n", + " ttheta_calc = job.get(\"output/ttheta_calc\")\n", + " int_calc = job.get(\"output/int_calc\")\n", + " \n", + " ttheta_uff_calc = job_uff.get(\"output/ttheta_calc\")\n", + " int_uff_calc = job_uff.get(\"output/int_calc\")\n", + "\n", + " ttheta_ref = job.get(\"output/ttheta_ref\")\n", + " int_ref = job.get(\"output/int_ref\")\n", + " \n", + " stat_res = job.compare(np.array([ttheta_calc,int_calc]).T, np.array([ttheta_ref,int_ref]).T,scale='optimal',verbose=False,plot=False)\n", + " stat_res_uff = job.compare(np.array([ttheta_uff_calc,int_uff_calc]).T, np.array([ttheta_ref,int_ref]).T,scale='optimal',verbose=False,plot=False)\n", + " \n", + " print(name)\n", + " pt.figure()\n", + " pt.plot(ttheta_calc, (int_calc-int_calc.min())*stat_res['scalefactor'],lw=1,label='QuickFF, {:3.2f}'.format(stat_res['R_wp']))\n", + " pt.plot(ttheta_uff_calc,(int_uff_calc-int_uff_calc.min())*stat_res_uff['scalefactor'],lw=1,label='UFF, {:3.2f}'.format(stat_res_uff['R_wp']))\n", + " pt.plot(ttheta_ref, (int_ref-int_ref.min()),lw=1,label='reference')\n", + " \n", + " ax1 = pt.gca()\n", + " ax1.set_xlabel('2θ (°)')\n", + " ax1.set_ylabel('Intensity (a.u.)')\n", + "\n", + " ax1.tick_params(\n", + " axis='y', # changes apply to the y-axis\n", + " which='both', # both major and minor ticks are affected\n", + " left=False, # ticks along the left edge are off\n", + " right=False, # ticks along the right edge are off\n", + " labelleft=False) # labels along the left edge are off\n", + "\n", + " ax1.ticklabel_format(axis='y', style='plain', useOffset=False) # avoid scientific notation and offsets\n", + " ax1.legend(bbox_to_anchor=(1.1,.5), loc='center left',frameon=False)\n", + " ax1.xaxis.set_major_locator(MaxNLocator(integer=True))\n", + " \n", + " pt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Background filtering" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.ticker import MaxNLocator\n", + "def remove_background(reference,locs=None,bkg_range=None,bkg_points=10,uniform=True,plot=False,fname=None,skiprows=0):\n", + " \"\"\"\n", + " Remove the background the provided reference pattern\n", + "\n", + " ***Args***\n", + " reference\n", + " file name of reference data (file with two columns, no header)\n", + " or a 2D array with (ttheta,intensity)\n", + " skiprows\n", + " number of rows to skip in reading file (when reference is a filename)\n", + " locs\n", + " node locations in degrees, default linspace in range with bkg_points\n", + " bkg_range\n", + " range for background points, default full ttheta range of reference pattern\n", + " bkg_points\n", + " set number of background points in fit\n", + " uniform\n", + " if false, the grid points are locally optimized towards the minima for a better interpolation\n", + " plot\n", + " plot the analysis and difference plot\n", + " fname\n", + " file location for the plot\n", + "\n", + " ***Returns***\n", + " a new reference pattern (ttheta,intensity) array\n", + " \"\"\"\n", + "\n", + " # Interactively use the pyobjcryst code here, it should not be necessary to run a separate job for this\n", + " try:\n", + " import pyobjcryst\n", + " except ImportError:\n", + " raise ImportError(\"The pyobjcryst package is required for this\")\n", + "\n", + " pp = pyobjcryst.powderpattern.PowderPattern()\n", + "\n", + " # Add the reference data\n", + " if isinstance(reference,str) and os.path.exists(reference):\n", + " pp.ImportPowderPattern2ThetaObs(reference,skiprows) # skip no lines\n", + " elif isinstance(reference,np.ndarray):\n", + " pp.SetPowderPatternX(reference[:,0]*deg)\n", + " pp.SetPowderPatternObs(reference[:,1])\n", + " ttheta = pp.GetPowderPatternX()/deg\n", + " reference = pp.GetPowderPatternObs()\n", + "\n", + " if ttheta.shape == (0,):\n", + " raise ValueError('Did not succeed in loading reference data. Try explicitely loading your data as an array instead.')\n", + "\n", + "\n", + " # Background\n", + " if locs is not None:\n", + " bx = np.array(locs)\n", + " # Keep first location index for clear plot\n", + " else:\n", + " if bkg_range is not None:\n", + " assert len(bkg_range)==2\n", + " bx=np.linspace(bkg_range[0],bkg_range[1],bkg_points)\n", + " else:\n", + " bx=np.linspace(ttheta.min(),ttheta.max(),bkg_points)\n", + " \n", + " # Keep first location index for clear plot\n", + " idx = [np.argmin(np.abs(ttheta-bxi)) for bxi in bx]\n", + " reference_idx = idx[0]\n", + " \n", + " # Adapt bx to minima of reference pattern in each neighbourhood (optional)\n", + " if locs is None and not uniform:\n", + " idx = [np.argmin(np.abs(ttheta-bxi)) for bxi in bx]\n", + " step = (idx[1] - idx[0])//4\n", + " for n in range(len(bx)):\n", + " mn = -step if idx[n]>step else 0\n", + " mx = step if n<(len(idx)-1) else 0\n", + " bx[n] = ttheta[idx[n]+ mn + np.argmin(reference[idx[n]+mn:idx[n]+mx])]\n", + " if n==0: reference_idx = idx[n]+ mn + np.argmin(reference[idx[n]+mn:idx[n]+mx])\n", + "\n", + " \n", + " bx*=deg\n", + " by=np.zeros(bx.shape)\n", + " b=pp.AddPowderPatternBackground()\n", + " b.SetInterpPoints(bx,by)\n", + " b.UnFixAllPar()\n", + " b.OptimizeBayesianBackground()\n", + "\n", + " no_bg = pp.GetPowderPatternObs()-pp.GetPowderPatternCalc()\n", + " no_bg -= np.min(no_bg)\n", + " \n", + "\n", + " # Plot the difference\n", + " if plot:\n", + " # Consider the difference for the delta plot\n", + " height = no_bg.max()-no_bg.min()\n", + "\n", + " fig = pt.figure(figsize=(10,6))\n", + " ax1 = pt.gca()\n", + " background = pp.GetPowderPatternCalc() - pp.GetPowderPatternCalc().min()\n", + " reference_pattern = reference - reference.min()\n", + " ax1.plot(ttheta,background - (background[reference_idx]- reference_pattern[reference_idx]),lw=1,label='background')\n", + " ax1.plot(ttheta,reference_pattern,lw=1,label='reference')\n", + " ax1.plot(ttheta,no_bg-height*0.1,color='g',lw=1,label=r'$\\Delta$')\n", + " ax1.set_xlabel('2θ (°)')\n", + " ax1.set_ylabel('Intensity (a.u.)')\n", + "\n", + " # Format the plot\n", + " lims = pt.xlim()\n", + " ax1.hlines(0,lims[0],lims[1],lw=0.1)\n", + " ax1.set_xlim(lims)\n", + "\n", + " lims = pt.ylim()\n", + " ax1.vlines(bx/deg,lims[0],lims[1],lw=0.1)\n", + "\n", + " ax1.tick_params(\n", + " axis='y', # changes apply to the y-axis\n", + " which='both', # both major and minor ticks are affected\n", + " left=False, # ticks along the left edge are off\n", + " right=False, # ticks along the right edge are off\n", + " labelleft=False) # labels along the left edge are off\n", + "\n", + " ax1.ticklabel_format(axis='y', style='plain', useOffset=False) # avoid scientific notation and offsets\n", + "\n", + " ax1.legend(bbox_to_anchor=(1.1,.5), loc='center left',frameon=False)\n", + " ax1.xaxis.set_major_locator(MaxNLocator(integer=True))\n", + "\n", + " if fname is None:\n", + " fig.tight_layout()\n", + " pt.show()\n", + " else:\n", + " pt.savefig(fname+'.pdf',bbox_inches='tight')\n", + " pt.close()\n", + "\n", + " if fname is not None:\n", + " with open(fname + '.tsv','w') as f:\n", + " for i in range(len(ttheta)):\n", + " f.write(\"{:4.3f}\\t{:10.8f}\\n\".format(ttheta[i],no_bg[i]))\n", + "\n", + " return np.array([ttheta, no_bg]).T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# SEQUENCE BELOW\n", + "'dia_30-03-12_06-12-12',\n", + "'kgm_31-03-04_01-02-04',\n", + "'hcb_34-11-10_01-06-10',\n", + "'ctn_33-11-02_29-01-02_None',\n", + "'bor_18-08-01_28-01-01_None',\n", + "'hcb_11-02-06_None',\n", + "'sql_24-01-01_10-08-01'," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 'dia_30-03-12_06-12-12'\n", + "bn = list(structures_database.keys())[0]\n", + "name = bn.replace('-','_')\n", + "job = pr_database_pxrd.load(name)\n", + "refpattern = np.array([job['output/ttheta_ref'],job['output/int_ref']]).T\n", + "\n", + "refpattern_1 = remove_background(refpattern, bkg_range=(2,10), bkg_points=2, uniform=True,plot=False)\n", + "refpattern_2 = remove_background(refpattern_1, bkg_points=3, uniform=True,plot=False)\n", + "refpattern_3 = remove_background(refpattern_2, bkg_range=(20,40), bkg_points=3, uniform=True,plot=False,fname='./input_files/PXRD/Structures_Database/'+bn+'/background_filtered')\n", + "#job.compare(refpattern,refpattern_3,scale=False,plot=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 'kgm_31-03-04_01-02-04'\n", + "bn = list(structures_database.keys())[1]\n", + "name = bn.replace('-','_')\n", + "print(name)\n", + "job = pr_database_pxrd.load(name)\n", + "#refpattern = np.array([job['output/ttheta_ref'],job['output/int_ref']]).T\n", + "\n", + "#refpattern_1 = remove_background(refpattern, bkg_range=(1,8), bkg_points=4, uniform=False,plot=True)\n", + "#refpattern_2 = remove_background(refpattern_1, bkg_points=3, uniform=True,plot=False)\n", + "#refpattern_3 = remove_background(refpattern_2, bkg_range=(20,40), bkg_points=3, uniform=True,plot=False,fname='./input_files/PXRD/Structures_Database/'+bn+'/background_filtered')\n", + "#job.compare(refpattern,refpattern_3,scale=False,plot=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 'hcb_34-11-10_01-06-10'\n", + "bn = list(structures_database.keys())[2]\n", + "name = bn.replace('-','_')\n", + "job = pr_database_pxrd.load(name)\n", + "refpattern = np.array([job['output/ttheta_ref'],job['output/int_ref']]).T\n", + "plot=False\n", + "\n", + "refpattern_1 = remove_background(refpattern, bkg_points=5, uniform=True,plot=plot)\n", + "refpattern_2 = remove_background(refpattern_1, bkg_range=(10,40), bkg_points=2, uniform=True,plot=plot)\n", + "refpattern_3 = remove_background(refpattern_2, bkg_points=5, uniform=True,plot=plot)\n", + "refpattern_4 = remove_background(refpattern_3, bkg_range=(10,40), bkg_points=3, uniform=True,plot=plot)\n", + "refpattern_5 = remove_background(refpattern_4, bkg_range=(35,40), bkg_points=3, uniform=True,plot=plot,fname='./input_files/PXRD/Structures_Database/'+bn+'/background_filtered')\n", + "#job.compare(refpattern,refpattern_5,scale=False,plot=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 'ctn_33-11-02_29-01-02_None'\n", + "bn = list(structures_database.keys())[3]\n", + "name = bn.replace('-','_')\n", + "job = pr_database_pxrd.load(name)\n", + "refpattern = np.array([job['output/ttheta_ref'],job['output/int_ref']]).T\n", + "\n", + "plot=True\n", + "# no fit necessary\n", + "#job.compare(refpattern,refpattern_5,scale=False,plot=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 'bor_18-08-01_28-01-01_None'\n", + "bn = list(structures_database.keys())[4]\n", + "name = bn.replace('-','_')\n", + "job = pr_database_pxrd.load(name)\n", + "refpattern = np.array([job['output/ttheta_ref'],job['output/int_ref']]).T\n", + "\n", + "plot=False\n", + "\n", + "refpattern_1 = remove_background(refpattern, bkg_points=3, uniform=True,plot=plot)\n", + "refpattern_2 = remove_background(refpattern_1, bkg_range=(9,35), bkg_points=4, uniform=True,plot=plot)\n", + "refpattern_3 = remove_background(refpattern_2, bkg_points=4, uniform=True,plot=plot,fname='./input_files/PXRD/Structures_Database/'+bn+'/background_filtered')\n", + "#job.compare(refpattern,refpattern_3,scale=False,plot=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 'hcb_11-02-06_None'\n", + "bn = list(structures_database.keys())[5]\n", + "name = bn.replace('-','_')\n", + "job = pr_database_pxrd.load(name)\n", + "refpattern = np.array([job['output/ttheta_ref'],job['output/int_ref']]).T\n", + "\n", + "plot=False\n", + "\n", + "refpattern_1 = remove_background(refpattern, bkg_points=4, uniform=True,plot=plot)\n", + "refpattern_2 = remove_background(refpattern_1, bkg_points=2, uniform=True,plot=plot)\n", + "refpattern_3 = remove_background(refpattern_2, bkg_points=4, uniform=True,plot=plot)\n", + "refpattern_4 = remove_background(refpattern_3, bkg_range=(16,60), bkg_points=2, uniform=True,plot=plot)\n", + "refpattern_5 = remove_background(refpattern_4, bkg_points=4, uniform=False,plot=plot)\n", + "refpattern_6 = remove_background(refpattern_5, bkg_range=(38,60), bkg_points=3, uniform=False,plot=plot,fname='./input_files/PXRD/Structures_Database/'+bn+'/background_filtered')\n", + "\n", + "#job.compare(refpattern,refpattern_6,scale=False,plot=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#\n", + "bn = list(structures_database.keys())[6]\n", + "name = bn.replace('-','_')\n", + "job = pr_database_pxrd.load(name)\n", + "refpattern = np.array([job['output/ttheta_ref'],job['output/int_ref']]).T\n", + "\n", + "plot=False\n", + "\n", + "refpattern_1 = remove_background(refpattern, bkg_range=(2,13), bkg_points=3, uniform=False,plot=plot)\n", + "refpattern_2 = remove_background(refpattern_1,bkg_range=(2,50), bkg_points=2, uniform=True,plot=plot)\n", + "refpattern_3 = remove_background(refpattern_2, bkg_points=4, uniform=True,plot=plot)\n", + "refpattern_4 = remove_background(refpattern_3, bkg_points=7, uniform=False,plot=plot)\n", + "refpattern_5 = remove_background(refpattern_4, bkg_range=(15,50), bkg_points=7, uniform=False,plot=plot)\n", + "refpattern_6 = remove_background(refpattern_5, bkg_points=4, uniform=False,plot=plot,fname='./input_files/PXRD/Structures_Database/'+bn+'/background_filtered')\n", + "\n", + "#job.compare(refpattern,refpattern_6,scale=False,plot=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Static w background filtering" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pr_database_pxrd_bg = pr_database.create_group('PXRD_bg')\n", + "pr_database_uff_pxrd_bg = pr_database_uff.create_group('PXRD_bg')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for bn in structures_database.keys():\n", + " name = bn.replace('-','_')\n", + " background_filtering = glob.glob('./input_files/PXRD/Structures_Database/'+bn+'/background_filtered.tsv')\n", + " if len(background_filtering)==0:\n", + " refpattern = glob.glob('./input_files/PXRD/Structures_Database/'+bn+'/*.dat')[0]\n", + " else:\n", + " refpattern = background_filtering[0]\n", + " job = pr_database.load(name)\n", + " job_uff = pr_database_uff.load(name)\n", + "\n", + " static_PXRD_job(pr_database_pxrd_bg,name,job.get_structure(),refpattern=refpattern)\n", + " static_PXRD_job(pr_database_uff_pxrd_bg,name,job_uff.get_structure(),refpattern=refpattern)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.ticker import MaxNLocator\n", + "\n", + "for bn in structures_database.keys():\n", + " name = bn.replace('-','_')\n", + " \n", + " job = pr_database_pxrd_bg.load(name)\n", + " job_uff = pr_database_uff_pxrd_bg.load(name)\n", + " \n", + " ttheta_calc = job.get(\"output/ttheta_calc\")\n", + " int_calc = job.get(\"output/int_calc\")\n", + " \n", + " ttheta_uff_calc = job_uff.get(\"output/ttheta_calc\")\n", + " int_uff_calc = job_uff.get(\"output/int_calc\")\n", + "\n", + " ttheta_ref = job.get(\"output/ttheta_ref\")\n", + " int_ref = job.get(\"output/int_ref\")\n", + " \n", + " stat_res = job.compare(np.array([ttheta_calc,int_calc]).T, np.array([ttheta_ref,int_ref]).T,scale='optimal',verbose=False,plot=False)\n", + " stat_res_uff = job.compare(np.array([ttheta_uff_calc,int_uff_calc]).T, np.array([ttheta_ref,int_ref]).T,scale='optimal',verbose=False,plot=False)\n", + " \n", + " print(name)\n", + " pt.figure()\n", + " pt.plot(ttheta_calc, (int_calc-int_calc.min())*stat_res['scalefactor'],lw=1,label='QuickFF, {:3.2f}'.format(stat_res['R_wp']))\n", + " pt.plot(ttheta_uff_calc,(int_uff_calc-int_uff_calc.min())*stat_res_uff['scalefactor'],lw=1,label='UFF, {:3.2f}'.format(stat_res_uff['R_wp']))\n", + " pt.plot(ttheta_ref, (int_ref-int_ref.min()),lw=1,label='reference')\n", + " \n", + " ax1 = pt.gca()\n", + " ax1.set_xlabel('2θ (°)')\n", + " ax1.set_ylabel('Intensity (a.u.)')\n", + "\n", + " ax1.tick_params(\n", + " axis='y', # changes apply to the y-axis\n", + " which='both', # both major and minor ticks are affected\n", + " left=False, # ticks along the left edge are off\n", + " right=False, # ticks along the right edge are off\n", + " labelleft=False) # labels along the left edge are off\n", + "\n", + " ax1.ticklabel_format(axis='y', style='plain', useOffset=False) # avoid scientific notation and offsets\n", + " ax1.legend(bbox_to_anchor=(1.1,.5), loc='center left',frameon=False)\n", + " ax1.xaxis.set_major_locator(MaxNLocator(integer=True))\n", + " \n", + " pt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Dynamic" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pr_database_md_pxrd = pr_database_md.create_group('PXRD')\n", + "pr_database_uff_md_pxrd = pr_database_uff_md.create_group('PXRD')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def dynamic_PXRD_job_v2(pr,name,md_job,refpattern=None,start=1000,num_frames=50,cores=1,cluster='victini',reservation_tag=None):\n", + " job = pr.create_job(pr.job_type.GPXRD, name, delete_existing_job=True)\n", + " job.load_trajectory(md_job,start=start,num_frames=num_frames)\n", + " \n", + " if refpattern is not None:\n", + " job.set_reference_pattern(refpattern)\n", + " \n", + " # We dont want individual fhkl indices\n", + " job.input['save_fhkl'] = False\n", + " \n", + " job.server.queue = cluster\n", + " job.server.cores = cores\n", + " job.server.run_time = 5*60*60 # in seconds\n", + " job.server.reservation_tag = reservation_tag\n", + "\n", + " job.run()\n", + " #return job" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for bn in structures_database.keys():\n", + " name = bn.replace('-','_')\n", + " refpattern = glob.glob('./input_files/PXRD/Structures_Database/'+bn+'/*.dat')[0]\n", + " job = pr_database_md.load(name)\n", + " job_uff = pr_database_uff_md.load(name)\n", + "\n", + " if job.status.finished:\n", + " dynamic_PXRD_job_v2(pr_database_md_pxrd,name,job,refpattern=refpattern,cluster='slaking',cores=1)\n", + " if job_uff.status.finished:\n", + " dynamic_PXRD_job_v2(pr_database_uff_md_pxrd,name,job_uff,refpattern=refpattern,cluster='slaking',cores=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.ticker import MaxNLocator\n", + "\n", + "for bn in structures_database.keys():\n", + " name = bn.replace('-','_')\n", + " \n", + " try:\n", + " job = pr_database_md_pxrd.load(name)\n", + " job_uff = pr_database_uff_md_pxrd.load(name)\n", + "\n", + " ttheta_calc = job.get(\"output/ttheta_calc\")\n", + " int_calc = job.get(\"output/int_calc\")\n", + "\n", + " ttheta_uff_calc = job_uff.get(\"output/ttheta_calc\")\n", + " int_uff_calc = job_uff.get(\"output/int_calc\")\n", + "\n", + " ttheta_ref = job.get(\"output/ttheta_ref\")\n", + " int_ref = job.get(\"output/int_ref\")\n", + " \n", + " except AttributeError:\n", + " continue\n", + " \n", + " stat_res = job.compare(np.array([ttheta_calc,int_calc]).T, np.array([ttheta_ref,int_ref]).T,scale='optimal',verbose=False,plot=False)\n", + " stat_res_uff = job.compare(np.array([ttheta_uff_calc,int_uff_calc]).T, np.array([ttheta_ref,int_ref]).T,scale='optimal',verbose=False,plot=False)\n", + " \n", + " print(name)\n", + " pt.figure(figsize=(20,10))\n", + " pt.plot(ttheta_calc, (int_calc-int_calc.min())*stat_res['scalefactor'],lw=1,label='QuickFF, {:3.2f}'.format(stat_res['R_wp']))\n", + " pt.plot(ttheta_uff_calc,(int_uff_calc-int_uff_calc.min())*stat_res_uff['scalefactor'],lw=1,label='UFF, {:3.2f}'.format(stat_res_uff['R_wp']))\n", + " pt.plot(ttheta_ref, (int_ref-int_ref.min()),lw=1,label='reference')\n", + " \n", + " ax1 = pt.gca()\n", + " ax1.set_xlabel('2θ (°)')\n", + " ax1.set_ylabel('Intensity (a.u.)')\n", + "\n", + " ax1.tick_params(\n", + " axis='y', # changes apply to the y-axis\n", + " which='both', # both major and minor ticks are affected\n", + " left=False, # ticks along the left edge are off\n", + " right=False, # ticks along the right edge are off\n", + " labelleft=False) # labels along the left edge are off\n", + "\n", + " ax1.ticklabel_format(axis='y', style='plain', useOffset=False) # avoid scientific notation and offsets\n", + " ax1.legend(bbox_to_anchor=(1.1,.5), loc='center left',frameon=False)\n", + " ax1.xaxis.set_major_locator(MaxNLocator(integer=True))\n", + " \n", + " pt.savefig(name+'.png',bbox_inches='tight')\n", + " pt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Dynamic w background filtering" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pr_database_md_pxrd_bg = pr_database_md.create_group('PXRD_bg')\n", + "pr_database_uff_md_pxrd_bg = pr_database_uff_md.create_group('PXRD_bg')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for bn in structures_database.keys():\n", + " name = bn.replace('-','_')\n", + " background_filtering = glob.glob('./input_files/PXRD/Structures_Database/'+bn+'/background_filtered.tsv')\n", + " if len(background_filtering)==0:\n", + " refpattern = glob.glob('./input_files/PXRD/Structures_Database/'+bn+'/*.dat')[0]\n", + " else:\n", + " refpattern = background_filtering[0]\n", + "\n", + " job = pr_database_md.load(name)\n", + " job_uff = pr_database_uff_md.load(name)\n", + "\n", + " if job.status.finished:\n", + " if pr_database_md_pxrd_bg.load(name) is None:\n", + " dynamic_PXRD_job_v2(pr_database_md_pxrd_bg,name,job,refpattern=refpattern,cluster='slaking',cores=1)\n", + " elif pr_database_md_pxrd_bg.load(name).status.aborted:\n", + " dynamic_PXRD_job_v2(pr_database_md_pxrd_bg,name,job,refpattern=refpattern,cluster='slaking',cores=8)\n", + " if job_uff.status.finished:\n", + " if pr_database_uff_md_pxrd_bg.load(name) is None:\n", + " dynamic_PXRD_job_v2(pr_database_uff_md_pxrd_bg,name,job_uff,refpattern=refpattern,cluster='slaking',cores=1)\n", + " elif pr_database_uff_md_pxrd_bg.load(name).status.aborted:\n", + " dynamic_PXRD_job_v2(pr_database_uff_md_pxrd_bg,name,job_uff,refpattern=refpattern,cluster='slaking',cores=8)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.ticker import MaxNLocator\n", + "\n", + "for bn in structures_database.keys():\n", + " name = bn.replace('-','_')\n", + " \n", + " try:\n", + " job = pr_database_md_pxrd_bg.load(name)\n", + " job_uff = pr_database_uff_md_pxrd_bg.load(name)\n", + "\n", + " ttheta_calc = job.get(\"output/ttheta_calc\")\n", + " int_calc = job.get(\"output/int_calc\")\n", + "\n", + " ttheta_uff_calc = job_uff.get(\"output/ttheta_calc\")\n", + " int_uff_calc = job_uff.get(\"output/int_calc\")\n", + "\n", + " ttheta_ref = job.get(\"output/ttheta_ref\")\n", + " int_ref = job.get(\"output/int_ref\")\n", + " \n", + " except AttributeError:\n", + " continue\n", + " \n", + " stat_res = job.compare(np.array([ttheta_calc,int_calc]).T, np.array([ttheta_ref,int_ref]).T,scale='optimal',verbose=False,plot=False)\n", + " stat_res_uff = job.compare(np.array([ttheta_uff_calc,int_uff_calc]).T, np.array([ttheta_ref,int_ref]).T,scale='optimal',verbose=False,plot=False)\n", + " \n", + " print(name)\n", + " pt.figure(figsize=(20,10))\n", + " pt.plot(ttheta_calc, (int_calc-int_calc.min())*stat_res['scalefactor'],lw=1,label='QuickFF, {:3.2f}'.format(stat_res['R_wp']))\n", + " pt.plot(ttheta_uff_calc,(int_uff_calc-int_uff_calc.min())*stat_res_uff['scalefactor'],lw=1,label='UFF, {:3.2f}'.format(stat_res_uff['R_wp']))\n", + " pt.plot(ttheta_ref, (int_ref-int_ref.min()),lw=1,label='reference')\n", + " \n", + " ax1 = pt.gca()\n", + " ax1.set_xlabel('2θ (°)')\n", + " ax1.set_ylabel('Intensity (a.u.)')\n", + "\n", + " ax1.tick_params(\n", + " axis='y', # changes apply to the y-axis\n", + " which='both', # both major and minor ticks are affected\n", + " left=False, # ticks along the left edge are off\n", + " right=False, # ticks along the right edge are off\n", + " labelleft=False) # labels along the left edge are off\n", + "\n", + " ax1.ticklabel_format(axis='y', style='plain', useOffset=False) # avoid scientific notation and offsets\n", + " ax1.legend(bbox_to_anchor=(1.1,.5), loc='center left',frameon=False)\n", + " ax1.xaxis.set_major_locator(MaxNLocator(integer=True))\n", + " \n", + " #pt.savefig(name+'.png',bbox_inches='tight')\n", + " pt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.10.12" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": { + "height": "calc(100% - 180px)", + "left": "10px", + "top": "150px", + "width": "276px" + }, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}