/******************************************************************************* * * McXtrace, X-ray tracing package * Copyright, All rights reserved * DTU Physics, Kgs. Lyngby, Denmark * Synchrotron SOLEIL, Saint-Aubin, France * * Component: Single_crystal * * %Identification * Written by: Kristian Nielsen * Date: December 1999 * Origin: Risoe * Modified by: EF, 22nd Apr 2003 : now uses Read_Table library * Modified by: EBK 2010: adapted for x-rays: * Modified by: M. Schulz, March 2012 : allow to curve the crystal planes * Modified by: EF, PW, May 2014: code efficiency improvement when SPLIT is used * Modified by: EF, PW, 2015: powder and texture mode * Modified by: PW, May 2017: Remove statement about being under validation * Modified by: PW, June 2017: Doc updates * Modified by: T. Weber, October 2017: new focus code * Modified by: PW, Feb 2018: GPU edits * Modified by: EF, May 2023: (re)port from n to X * Modified by: EF, June 2023: can now read CIF files, via cif2hkl * * Mosaic single crystal with multiple scattering vectors, optimised for speed * with large crystals and many reflections. * * %Description * Single crystal with mosaic. Delta-D/D option for finite-size effects. * Rectangular geometry. Multiple scattering and secondary extinction included. * The mosaic may EITHER be specified isotropic by setting the mosaic input * parameter, OR anisotropic by setting the mosaic_a, mosaic_b, and mosaic_c * parameters. * If you get strange results, check the mosaicity and delta(d)/d parameters, as this * component is not suited for ideal/perfect mosaic crystals. * The crystal lattice can be bent locally, keeping the external geometry unchanged. * Curvature is spherical along vertical and horizontal axes. * * Speed/stat optimisation using SPLIT * In order to dramatically improve the simulation efficiency, we recommend to * use a SPLIT keyword on this component (or prior to it), as well as to disable * the multiple scattering handling by setting order=1. This is especially powerful * for large reflection lists such as with macromolecular proteins. When an incoming * particle is identical to the preceeding, reciprocal space initialisation is * skipped, and a Monte Carlo choice is done on available reflections from the last * repciprocal space calculation! To assist the user in choosing a "relevant" value * of the SPLIT, a rolling average of the number of available reflections is * calculated and presented in the component output. * * Mosacitiy modes * The component features three independent ways of parametrising mosaicity: * a) The original algorithm where mosaicity is implemented by extending each * reflection by a Gaussian "cigar" in reciprocal space, characterised by * the parameters mosaic and delta_d_d. * (Also known as "isotropic mosaicity"). * b) A similar mode where mosaicities can be non-isotropic and given as the * parameters mosaic_a, mosaic_b and mosaic_c, around the unit cell axes. * (Also known as "anisotropic mosaicity"). * c) Given two "macroscopically"/experimentally measured width/mosaicities * of two independent reflections, parametrised by the list * mosaic_AB = {mos_a, mos_b, a_h, a_k, a_l, b_h, b_k, b_l}, a set of * microscopic mosaicities as in b) are estimated (internally) and applied. * (Also known as "phenomenological mosaicity"). * * Powder-mode * When the powder mode is used (powder=0-1), a randomised transformation * of the particle direction is made before and after scattering, thereby letting * the single crystal behave as a crystallite of either a powder (crystallite * orientation fully randomised). * * Curved crystal mode * The component features a method to curve the lattice planes slightly with respect * to the outer geometry of the crystal. The method is implemented as a transformation * on the particle direction vector, and should be used only in cases where: * a) The reflection lattice vector is ~ orthogonal to the crystal surface. * b) The modelled curvarture is "small" with respect to the crystal surface. * * Sample shape * Sample shape may be a cylinder, a sphere, a box or any other shape: * box/plate: xwidth x yheight x zdepth * cylinder: radius x yheight * sphere: radius (yheight=0) * any shape: geometry=OFF/PLY file * * The complex geometry option handles any closed non-convex polyhedra. * It computes the intersection points of the photon ray with the object * transparently, so that it can be used like a regular sample object. * It supports the PLY, OFF and NOFF file format but not COFF (colored faces). * Such files may be generated from XYZ data using: * qhull < coordinates.xyz Qx Qv Tv o > geomview.off * or * powercrust coordinates.xyz * and viewed with geomview or java -jar jroff.jar (see below). * The default size of the object depends on the OFF/PLY file data, but its * bounding box may be resized using xwidth,yheight and zdepth. * * Crystal definition file format * Crystal structure is specified with an ascii data file. Each line contains * 4 or more numbers, separated by white spaces: * * h k l ... F2 * * The first three numbers are the (h,k,l) indices of the reciprocal lattice * point, and the 7-th number is the value of the structure factor |F|**2, in * barns. The rest of the numbers are not used; the file is in the format * output by the Crystallographica program. * The reflection list should be ordered by decreasing d-spacing values. * Lines begining by '#' are read as comments (ignored). Most sample parameters * may be defined from the data file header, following the same mechanism as * PowderN. * * Current data file header keywords include, for data format specification: * #column_h * #column_k * #column_l * #column_F2 * #column_F * and for material specification: * #sigma_inc * #Delta_d/d * #lattice_a * #lattice_b * #lattice_c * #lattice_aa * #lattice_bb * #lattice_cc * * Last, CIF, FullProf and ShelX files can be read, and converted to F2(hkl) lists * when 'cif2hkl' is installed. The CIF2HKL env variable can be used to point to a * proper executable, else the McCode or the system installed versions are used. * * Satellite Bragg peaks - surface crystal truncation rods (CTR) * It is known that scattering from a finite crystal introduces a broadening of * Bragg peaks, seen in surface diffraction at grazing angle [Robinson and Tweet, * Rep. Prog. Phys. 55 (1992) 599)]. The CTR is specified as two vectors, which * hold the squared Fourier transform of the crystal geometry, for instance: * bulk: Dirac peak (FT of infinity) * half-bulk: Dirac+1/k^2 (FT of half plane, k in rlu) * layer: sinc(PI*k*d)^2 (FT of a top-hat, thickness 'd') * These vectors are given as 'surf_k' [in 1/Angs] and 'surf_FT2', both of length * 'surf_size'. The CTR is to be applied along vector 'surf_dir' which indicates * the surface normal {nx,ny,nz} in real space. When surf_size=-1, the component * sets the proper truncation function (only for thin box and disk shapes). * This feature is an approximation of the real surface scattering, and does not * handle complex geometries (clusters, wetting, multi-layers, ...). It may be * used as well to describe over-structure satellite peaks. * * Example: Single_crystal(xwidth=0.01, yheight=0.01, zdepth=0.01, mosaic = 5, * reflections="Si.lau") * * A diamond crystal plate, cut for (002) reflections * Single_crystal(xwidth = 0.002, yheight = 0.1, zdepth = 0.1, * mosaic = 5, delta_d_d=3e-4, reflections = "C-diamond.lau", * ax=0, ay=2.14, az=-1.24, * bx = 0, by = 0, bz = 2.47, * cx = 6.71, cy = 0, cz = 0) * * A adrenaline protein * Single_crystal(xwidth=0.005, yheight=0.005, zdepth=0.005, * mosaic = 5, reflections="adrenaline.lau") * * Also, always use a non-zero value of delta_d_d. * * This sample component can advantageously benefit from the SPLIT feature, e.g. * SPLIT COMPONENT sx = Single_crystal(...) * * %Parameters * INPUT PARAMETERS * radius: [m] Outer radius of sample in (x,z) plane * xwidth: [m] Width of crystal * yheight: [m] Height of crystal * zdepth: [m] Depth of crystal (no extinction simulated) * geometry: [str] Name of an Object File Format (OFF) or PLY file for complex geometry. The OFF/PLY file may be generated from XYZ coordinates using qhull/powercrust * delta_d_d: [1] Lattice spacing variance, gaussian RMS (longitudinal mosaic) e.g. 1e-4 to 1e-3. * mosaic: [arc minutes] Crystal mosaic (isotropic), gaussian RMS. Puts the crystal in the isotropic mosaic model state, thus disregarding other mosaicity parameters, e.g. 1-10. * mosaic_a: [arc minutes] Horizontal (rotation around lattice vector a) mosaic (anisotropic), gaussian RMS. Put the crystal in the anisotropic crystal vector state. i.e. model mosaicity through rotation around the crystal lattice vectors. Has precedence over in-plane mosaic model. * mosaic_b: [arc minutes] Vertical (rotation around lattice vector b) mosaic (anisotropic), gaussian RMS. * mosaic_c: [arc minutes] Out-of-plane (Rotation around lattice vector c) mosaic (anisotropic), gaussian RMS * mosaic_AB: [arc_minutes, arc_minutes,1, 1, 1, 1, 1, 1] In-plane mosaic rotation and plane vectors (anisotropic), mosaic_A, mosaic_B, A_h,A_k,A_l, B_h,B_k,B_l. Puts the crystal in the in-plane mosaic state. Vectors A and B define plane in which the crystal roation is defined, and mosaic_A, mosaic_B, denotes the resp. mosaicities (gaussian RMS) with respect to the two reflections chosen by A and B (Miller indices). * * recip_cell: [1] Choice of direct/reciprocal (0/1) unit cell definition * ax: [AA or AA^-1] Coordinates of first (direct/recip) unit cell vector * ay: [AA or AA^-1] a on y axis * az: [AA or AA^-1] a on z axis * bx: [AA or AA^-1] Coordinates of second (direct/recip) unit cell vector * by: [AA or AA^-1] b on y axis * bz: [AA or AA^-1] b on z axis * cx: [AA or AA^-1] Coordinates of third (direct/recip) unit cell vector * cy: [AA or AA^-1] c on y axis * cz: [AA or AA^-1] c on z axis * reflections: [string] File name containing structure factors of reflections (LAZ LAU CIF, FullProf, ShelX). Use empty ("") or NULL for incoherent scattering only * order: [1] Limit multiple scattering up to given order (0: all, 1: first, 2: second, ...) * extra_order: [1] When using order, allow additional multiple scattering without coherent scattering, sensible with very large unit cells (0: disable, 1: one extra, 2: two extra, ...) * * Optional input parameters * * p_transmit: [1] Monte Carlo probability for photons to be transmitted without any scattering. Used to improve statistics from weak reflections * sigma_inc: [barns] Incoherent scattering cross-section per unit cell (uniform). Fully isotropic and constant. Use -1 to inactivate * aa: [deg] Unit cell angles alpha, beta and gamma. Then uses norms of vectors a,b and c as lattice parameters * bb: [deg] Beta angle * cc: [deg] Gamma angle * barns: [1] Flag to indicate if |F|^2 from 'reflections' is in barns or fm^2. barns=1 for laz/cif and isotropic constant elastic scattering (reflections=NULL), barns=0 for lau type files * RX: [m] Radius of horizontal along X lattice curvature. flat for 0 * RY: [m] Radius of vertical along Y lattice curvature. flat for 0 * powder: [1] Flag to indicate powder mode, for simulation of Debye-Scherrer cones via random crystallite orientation. A powder texture can be approximated with powder within 0-1 * deltak: [AA-1] Equality-threshold for use in SPLIT settings. If difference between all ki_{x,y,z} are less than deltak from previous particle, the two are considered alike enough to jump directly to the MC choice between 'active' reflections * material_datafile: [Be.txt] File where the material parameters for the absorption may be found. Format is similar to what may be found off the NIST website. * surf_size: [1] Length of the surf_k and surf_FT vectors. When set as -1, CTR is automatically set for the box/thin disk geometry. * surf_k: [1/Angs] Momentum 'k' distribution around 0 for the CTR, length 'surf_size'. * surf_FT2: [1] Intensity |FT(real space)|^2 distribution as CTR of a single Bragg peak, length 'surf_size'. * surf_dir: [nx,ny,nz] Normal vector in real space for the truncation rod spread direction. {1,0,0} is for an YZ surface, with normal along X. * * CALCULATED PARAMETERS: * * hkl_info: [structure] Internal * hkl_info.type: interaction type of event 't'=Transmit, 'i'=Incoherent, 'c'=Coherent [char] * hkl_info.h: * hkl_info.k: wavevector indices of last coherent scattering event [Angs-1] * hkl_info.l: * * %Link * See ICSD Inorganic Crystal Structure Database * %Link * Web Elements * %Link * Fullprof powder refinement * %Link * Crystallographica software * %Link * Geomview and Object File Format (OFF) * %Link * Java version of Geomview (display only) jroff.jar * %Link * qhull * %Link * powercrust * %Link * material datafile obtained from http://physics.nist.gov/cgi-bin/ffast/ffast.pl * %Link * cif2hkl https://gitlab.com/soleil-data-treatment/soleil-software-projects/cif2hkl * * %End ****************************************************************************/ /* Overview of algorithm: (1). The photon intersects the crystal at (x,y,z) with given incoming wavevector ki=(kix,kiy,kiz). (2). Every reciprocal lattice point tau of magnitude less than 2*ki is considered for scattering. The scattering probability is the area of the intersection of the Ewald sphere (approximated by the tangential plane) with the 3-D Gaussian mosaic of the point tau. (3). The total coherent scattering cross section is computed as the sum over all tau. Together with the absorption and incoherent scattering cross sections and known potential flight-length l_full through the sample, we can compute the probability of the four events absorption, coherent scattering, incoherent scattering, and transmission. (4). absorption is never simulated explicitly, just incorporated in the photon weight. (5). Transmission in the first event is selected with the Monte Carlo probability p_transmit, which defaults to the actual transmission probability. After the first event, transmission is selected with the correct Monte Carlo probability. (6). Incoherent scattering is done simply by selecting a random direction for the outgoing wave vector kf. (7). For coherent scattering, a reciprocal lattice point is selected using the relative probabilities computed in (2), and the weight is adjusted with the contribution from the structure factors (this way all reflections will get equally good statistics in the detector). (8). The outgoing wave vector direction is picked at random using the intersecting 2-D Gauss computed in (2). The vector is normalized to the length of ki (elastic scattering) to account for the error caused by the planar approximation of the Ewald sphere. (9). The process is repeated from (2) with kf as new initial wave vector ki. */ DEFINE COMPONENT Single_crystal SETTING PARAMETERS(string reflections=0, string geometry=0, vector mosaic_AB={0,0, 0,0,0, 0,0,0}, xwidth=0, yheight=0, zdepth=0, radius=0, delta_d_d=1e-3, mosaic = -1, mosaic_a = -1, mosaic_b = -1, mosaic_c = -1, recip_cell=0, barns=0, ax = 0, ay = 0, az = 0, bx = 0, by = 0, bz = 0, cx = 0, cy = 0, cz = 0, p_transmit = 0.001, sigma_inc = 0, aa=0, bb=0, cc=0, order=1, extra_order=0, RX=0, RY=0, powder=0, deltak=1e-6, string material_datafile="Si.txt", int surf_size=0, vector surf_k=NULL, vector surf_FT2=NULL, vector surf_dir={0,0,0}) /* photon parameters: (x,y,z,kx,ky,kz,phi,t,Ex,Ey,Ez,p) */ SHARE %{ /* used for reading data table from file */ %include "read_table-lib" %include "interoff-lib" /* Declare structures and functions only once in each instrument. */ #ifndef SINGLE_CRYSTAL_DECL #define SINGLE_CRYSTAL_DECL #ifndef Mosaic_AB_Undefined #define Mosaic_AB_Undefined {0,0, 0,0,0, 0,0,0} #endif #ifndef MCSX_REFL_SLIST_SIZE #define MCSX_REFL_SLIST_SIZE 128 #endif struct hkl_data { int h, k, l; /* Indices for this reflection */ double F2; /* Value of structure factor */ double tau_x, tau_y, tau_z; /* Coordinates in reciprocal space */ double tau; /* Length of (tau_x, tau_y, tau_z) */ double u1x, u1y, u1z; /* First axis of local coordinate system */ double u2x, u2y, u2z; /* Second axis of local coordinate system */ double u3x, u3y, u3z; /* Third axis of local coordinate system */ double sig123; /* The product sig1*sig2*sig3 = volume of spot */ double m1, m2, m3; /* Diagonal matrix representation of Gauss */ double cutoff; /* Cutoff value for Gaussian tails */ }; struct tau_data { int index; /* Index into reflection table */ double refl; double xsect; /* The following vectors are in local koordinates. */ double rho_x, rho_y, rho_z; /* The vector ki - tau */ double rho; /* Length of rho vector */ double ox, oy, oz; /* Origin of Ewald sphere tangent plane */ double b1x, b1y, b1z; /* Spanning vectors of Ewald sphere tangent */ double b2x, b2y, b2z; double l11, l12, l22; /* Cholesky decomposition L of 2D Gauss */ double y0x, y0y; /* 2D Gauss center in tangent plane */ }; struct hkl_info_struct { int count; /* Number of reflections */ double m_delta_d_d; /* Delta-d/d FWHM */ double m_ax, m_ay, m_az; /* First unit cell axis (direct space, AA) */ double m_bx, m_by, m_bz; /* Second unit cell axis */ double m_cx, m_cy, m_cz; /* Third unit cell axis */ double asx, asy, asz; /* First reciprocal lattice axis (1/AA) */ double bsx, bsy, bsz; /* Second reciprocal lattice axis */ double csx, csy, csz; /* Third reciprocal lattice axis */ double m_a, m_b, m_c; /* length of lattice parameter lengths */ double m_aa, m_bb, m_cc; /* lattice angles */ double sigma_i; /* inc X sect */ double rho; /* density [at/Angs^3] */ double mat_weight; /* material atomic weight [g/mol] */ double mat_density; /* material atomic density [g/cm3] */ double V0; /* Unit cell volume (AA**3) */ int column_order[5]; /* column signification [h,k,l,F,F2] */ int recip; /* Flag to indicate if recip or direct cell axes given */ int shape; /* 0:cylinder, 1:box, 2:sphere 3:any shape*/ int flag_warning; /* number of warnings */ int flag_barns; /* 1: F2 in barns, 0: in fm^2 */ char type; /* type of last event: t=transmit,c=coherent or i=incoherent */ int h, k, l; /* last coherent scattering momentum transfer indices */ int tau_count; /* Number of reflections within cutoff */ double coh_refl, coh_xsect; /* cross section computed with last tau_list */ double kix, kiy, kiz; /* last incoming photon ki */ int nb_reuses, nb_refl, nb_refl_count; int max_tau_count; int ctr_size; double* ctr_k; double* ctr_FT2; double* ctr_dir; char compname[256]; }; #pragma acc routine // sign = SX_list_compare(a,b) // used for qsort, to sort reflections // input: a,b: two 'line_data' reflection pointers. // output: -1,0,1 for ab int SX_list_compare (void const* a, void const* b) { struct hkl_data const* pa = a; struct hkl_data const* pb = b; /* Sort by tau */ if (pa->tau < pb->tau) return -1; if (pa->tau > pb->tau) return 1; /* Sort by tau_x */ if (pa->tau_x < pb->tau_x) return -1; if (pa->tau_x > pb->tau_x) return 1; /* Sort by tau_y */ if (pa->tau_y < pb->tau_y) return -1; if (pa->tau_y > pb->tau_y) return 1; /* Sort by tau_z */ if (pa->tau_z < pb->tau_z) return -1; if (pa->tau_z > pb->tau_z) return 1; /* In case of tie, sort by F2 also */ if (pa->F2 < pb->F2) return -1; if (pa->F2 > pb->F2) return 1; return 0; } /* SX_list_compare */ #ifndef CIF2HKL #define CIF2HKL // hkl_filename = cif2hkl(file, options) // used to convert CIF/CFL/INS file into F2(hkl) // the CIF2HKL env var can point to a cif2hkl executable // else the McCode binary is attempted, then the system. char* cif2hkl (char* infile, char* options) { char cmd[1024]; int ret = 0; int found = 0; char* OUTFILE; char* inpath; // get filename extension const char* ext = strrchr (infile, '.'); if (!ext || ext == infile) return infile; else ext++; // return input when no extension or not a CIF/FullProf/ShelX file if (strcasecmp (ext, "cif") && strcasecmp (ext, "pcr") && strcasecmp (ext, "cfl") && strcasecmp (ext, "shx") && strcasecmp (ext, "ins") && strcasecmp (ext, "res")) return infile; OUTFILE = malloc (1024); if (!OUTFILE) { free (OUTFILE); return infile; } inpath = malloc (1024); if (!inpath) { free (OUTFILE); free (inpath); return infile; } // get input file path from read-table:Open_File FILE* f_infile = Open_File (infile, "r", inpath); if (!f_infile) { free (OUTFILE); free (inpath); free (f_infile); return infile; } fclose (f_infile); strncpy (OUTFILE, tmpnam (NULL), 1024); // create an output temporary file name // try in order the CIF2HKL env var, then the system cif2hkl, then the McCode one if (!found && getenv ("CIF2HKL")) { snprintf (cmd, 1024, "%s -o %s %s %s", getenv ("CIF2HKL"), OUTFILE, options, inpath); ret = system (cmd); if (ret != -1 && ret != 127) found = 1; } if (!found) { // try with cif2hkl command from the system PATH snprintf (cmd, 1024, "%s -o %s %s %s", "cif2hkl", OUTFILE, options, inpath); ret = system (cmd); if (ret != -1 && ret != 127) found = 1; } if (!found) { // As a last resort, attempt with cif2hkl from $MCXTRACE/bin snprintf (cmd, 1024, "%s%c%s%c%s -o %s %s %s", getenv (FLAVOR_UPPER) ? getenv (FLAVOR_UPPER) : MCXTRACE, MC_PATHSEP_C, "bin", MC_PATHSEP_C, "cif2hkl", OUTFILE, options, inpath); ret = system (cmd); } // ret = -1: child process could not be created // ret = 127: shell could not be executed in the child process if (ret == -1 || ret == 127) { free (OUTFILE); return (NULL); } // test if the result file has been created FILE* file = fopen (OUTFILE, "r"); if (!file) { free (OUTFILE); return (NULL); } MPI_MASTER (printf ("%s: INFO: Converting %s into F2(HKL) list %s\n", __FILE__, infile, OUTFILE); printf ("%s\n", cmd);); fflush (NULL); fclose (file); return (OUTFILE); } // cif2hkl #endif /* ------------------------------------------------------------------------ */ // count = read_hkl_data(filename, &info, ...) // read the given file to store F2(HKL) reflections and metadata. int read_hkl_data (char* SC_file, struct hkl_info_struct* info, struct hkl_data** hkl_list, double SC_mosaic, double SC_mosaic_a, double SC_mosaic_b, double SC_mosaic_c, double* SC_mosaic_AB, int surf_size, double* surf_k, double* surf_FT2, double* surf_dir) { struct hkl_data* list = NULL; int size = 0; t_Table sTable; /* sample data table structure from SC_file */ int i = 0, j = 0, J = 0; double tmp_x, tmp_y, tmp_z; char** parsing; char flag = 0; double nb_atoms = 1; double sum_F2 = 0; char* filename = NULL; if (!SC_file || !strlen (SC_file) || !strcmp (SC_file, "NULL") || !strcmp (SC_file, "0")) return (0); filename = cif2hkl (SC_file, "--xtal --mode XRA"); if (filename && filename != SC_file) info->flag_barns = 1; // cif2hkl returns barns if (!filename || Table_Read (&sTable, filename, 1) < 0) exit (fprintf (stderr, "Single_crystal: %s: ERROR: Could not open file %s - exiting!\n", info->compname, filename)); // remove temporary F2(hkl) file when giving CFL/CIF/ShelX file if (filename && filename != SC_file) unlink (filename); if (sTable.columns < 4) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: The number of columns in %s should be at least %d for [h,k,l,F2]\n", info->compname, SC_file, 4);); return (0); } if (!sTable.rows) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: The number of rows in %s should be at least %d\n", info->compname, SC_file, 1);); return (0); } else size = sTable.rows; /* parsing of header */ parsing = Table_ParseHeader (sTable.header, "sigma_abs", "sigma_a ", // 0-1 "sigma_inc", "sigma_i ", "column_h", "column_k", "column_l", "column_F ", "column_F2", "Delta_d/d", // 9 "lattice_a ", "lattice_b ", "lattice_c ", "lattice_aa", "lattice_bb", "lattice_cc", "nb_atoms", "multiplicity", // 17 "Vc", "V_0", "V_rho", "density", "weight", NULL); if (parsing) { // if (parsing[0] && !info->sigma_a) info->sigma_a=atof(parsing[0]); // if (parsing[1] && !info->sigma_a) info->sigma_a=atof(parsing[1]); if (parsing[2] && !info->sigma_i) info->sigma_i = atof (parsing[2]); if (parsing[3] && !info->sigma_i) info->sigma_i = atof (parsing[3]); if (parsing[4]) info->column_order[0] = atoi (parsing[4]); if (parsing[5]) info->column_order[1] = atoi (parsing[5]); if (parsing[6]) info->column_order[2] = atoi (parsing[6]); if (parsing[7]) info->column_order[3] = atoi (parsing[7]); if (parsing[8]) info->column_order[4] = atoi (parsing[8]); if (parsing[9] && info->m_delta_d_d <= 0) info->m_delta_d_d = atof (parsing[9]); if (parsing[10] && !info->m_a) info->m_a = atof (parsing[10]); if (parsing[11] && !info->m_b) info->m_b = atof (parsing[11]); if (parsing[12] && !info->m_c) info->m_c = atof (parsing[12]); if (parsing[13] && !info->m_aa) info->m_aa = atof (parsing[13]); if (parsing[14] && !info->m_bb) info->m_bb = atof (parsing[14]); if (parsing[15] && !info->m_cc) info->m_cc = atof (parsing[15]); if (parsing[16]) nb_atoms = atof (parsing[16]); if (parsing[17]) nb_atoms = atof (parsing[17]); if (parsing[18] && !info->rho) info->rho = 1 / atof (parsing[18]); if (parsing[19] && !info->rho) info->rho = 1 / atof (parsing[19]); if (parsing[20] && !info->rho) info->rho = atof (parsing[20]); if (parsing[21] && !info->mat_density) info->mat_density = atof (parsing[21]); if (parsing[22] && !info->mat_weight) info->mat_weight = atof (parsing[22]); for (i = 0; i <= 22; i++) if (parsing[i]) free (parsing[i]); free (parsing); } if (!info->column_order[0] || !info->column_order[1] || !info->column_order[2]) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Wrong h,k,l column definition\n", SC_file);); return (0); } if (!info->column_order[3] && !info->column_order[4]) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Wrong F,F2 column definition\n", SC_file);); return (0); } if (!info->rho && info->mat_density > 0 && info->mat_weight > 0 && nb_atoms > 0) { /* molar volume [cm^3/mol] = weight [g/mol] / density [g/cm^3] */ /* atom density per Angs^3 = [mol/cm^3] * N_Avogadro *(1e-8)^3 */ info->rho = info->mat_density / (info->mat_weight * nb_atoms) / 1e24 * NA; } if (nb_atoms > 1) { info->sigma_i *= nb_atoms; } if (info->m_delta_d_d <= 0) { info->m_delta_d_d = SC_mosaic > 0 ? SC_mosaic * MIN2RAD : MIN2RAD; // MIN2RAD = 3e-4 MPI_MASTER ( fprintf (stderr, "Single_crystal: %s: WARNING: Delta_d/d is <=0 (longitudinal mosaic). Setting delta_d_d=%g\n", info->compname, info->m_delta_d_d);); } else if (SC_mosaic < 0 && SC_mosaic_a < 0 && SC_mosaic_b < 0 && SC_mosaic_c < 0) { if (info->m_delta_d_d > 0) SC_mosaic = info->m_delta_d_d / MIN2RAD; } /* special cases for the structure definition */ if (info->m_ax || info->m_ay || info->m_az) { info->m_a = 0; info->m_aa = 0; } /* means we specify by hand the vectors */ if (info->m_bx || info->m_by || info->m_bz) { info->m_b = 0; info->m_bb = 0; } if (info->m_cx || info->m_cy || info->m_cz) { info->m_c = 0; info->m_cc = 0; }; /* compute the norm from vector a if missing */ if (info->m_ax || info->m_ay || info->m_az) { double as = sqrt (info->m_ax * info->m_ax + info->m_ay * info->m_ay + info->m_az * info->m_az); if (!info->m_bx && !info->m_by && !info->m_bz) info->m_a = info->m_b = as; if (!info->m_cx && !info->m_cy && !info->m_cz) info->m_a = info->m_c = as; } if (info->m_a && !info->m_b) info->m_b = info->m_a; if (info->m_b && !info->m_c) info->m_c = info->m_b; /* compute the lattive angles if not set from data file. Not used when in vector mode. */ if (info->m_a && !info->m_aa) info->m_aa = 90; if (info->m_aa && !info->m_bb) info->m_bb = info->m_aa; if (info->m_bb && !info->m_cc) info->m_cc = info->m_bb; /* parameters consistency checks */ if (!info->m_ax && !info->m_ay && !info->m_az && !info->m_a) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Wrong a lattice vector definition\n", info->compname);); return (0); } if (!info->m_bx && !info->m_by && !info->m_bz && !info->m_b) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Wrong b lattice vector definition\n", info->compname);); return (0); } if (!info->m_cx && !info->m_cy && !info->m_cz && !info->m_c) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Wrong c lattice vector definition\n", info->compname);); return (0); } if (info->m_aa && info->m_bb && info->m_cc && info->recip) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Selecting reciprocal cell and angles is unmeaningful.\n", info->compname);); return (0); } /* when lengths a,b,c + angles are given (instead of vectors a,b,c) */ if (info->m_aa && info->m_bb && info->m_cc) { double as, bs, cs; if (info->m_a) as = info->m_a; else as = sqrt (info->m_ax * info->m_ax + info->m_ay * info->m_ay + info->m_az * info->m_az); if (info->m_b) bs = info->m_b; else bs = sqrt (info->m_bx * info->m_bx + info->m_by * info->m_by + info->m_bz * info->m_bz); if (info->m_c) cs = info->m_c; else cs = sqrt (info->m_cx * info->m_cx + info->m_cy * info->m_cy + info->m_cz * info->m_cz); /* Real-space vectors with b along z and a in the yz-plane at the angle gamma to b; c then has cos(alpha) along b, cos(beta) against a, and a positive x component, which makes (a, b, c) right-handed. Each vector carries its own length: a = as, b = bs, c = cs. (Before, a was given the length bs and b the length as, and alpha/beta were exchanged in c, so for a != b every |F(h,k,l)|^2 ended up at the position of (k,h,l).) */ info->m_bx = 0; info->m_by = 0; info->m_bz = bs; info->m_ax = 0; info->m_ay = as * sin (info->m_cc * DEG2RAD); info->m_az = as * cos (info->m_cc * DEG2RAD); info->m_cz = cs * cos (info->m_aa * DEG2RAD); info->m_cy = cs * (cos (info->m_bb * DEG2RAD) - cos (info->m_cc * DEG2RAD) * cos (info->m_aa * DEG2RAD)) / sin (info->m_cc * DEG2RAD); info->m_cx = sqrt (cs * cs - info->m_cz * info->m_cz - info->m_cy * info->m_cy); MPI_MASTER (printf ("Single_crystal: %s: %s structure a=%g b=%g c=%g aa=%g bb=%g cc=%g ", info->compname, SC_file, as, bs, cs, info->m_aa, info->m_bb, info->m_cc);); } else { MPI_MASTER ( if (!info->recip) { printf ("Single_crystal: %s: %s structure a=[%g,%g,%g] b=[%g,%g,%g] c=[%g,%g,%g] ", info->compname, SC_file, info->m_ax, info->m_ay, info->m_az, info->m_bx, info->m_by, info->m_bz, info->m_cx, info->m_cy, info->m_cz); } else { printf ("Single_crystal: %s: %s structure a*=[%g,%g,%g] b*=[%g,%g,%g] c*=[%g,%g,%g] ", info->compname, SC_file, info->m_ax, info->m_ay, info->m_az, info->m_bx, info->m_by, info->m_bz, info->m_cx, info->m_cy, info->m_cz); }); } /* Compute reciprocal or direct lattice vectors. */ if (!info->recip) { vec_prod (tmp_x, tmp_y, tmp_z, info->m_bx, info->m_by, info->m_bz, info->m_cx, info->m_cy, info->m_cz); info->V0 = fabs (scalar_prod (info->m_ax, info->m_ay, info->m_az, tmp_x, tmp_y, tmp_z)); if (!info->V0 || isnan (info->V0)) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Lattice cell volume is found invalid V0=%g. Check crystal cell parameters.\n", info->compname, info->V0)); return (0); } info->asx = 2 * PI / info->V0 * tmp_x; info->asy = 2 * PI / info->V0 * tmp_y; info->asz = 2 * PI / info->V0 * tmp_z; vec_prod (tmp_x, tmp_y, tmp_z, info->m_cx, info->m_cy, info->m_cz, info->m_ax, info->m_ay, info->m_az); info->bsx = 2 * PI / info->V0 * tmp_x; info->bsy = 2 * PI / info->V0 * tmp_y; info->bsz = 2 * PI / info->V0 * tmp_z; vec_prod (tmp_x, tmp_y, tmp_z, info->m_ax, info->m_ay, info->m_az, info->m_bx, info->m_by, info->m_bz); info->csx = 2 * PI / info->V0 * tmp_x; info->csy = 2 * PI / info->V0 * tmp_y; info->csz = 2 * PI / info->V0 * tmp_z; } else { info->asx = info->m_ax; info->asy = info->m_ay; info->asz = info->m_az; info->bsx = info->m_bx; info->bsy = info->m_by; info->bsz = info->m_bz; info->csx = info->m_cx; info->csy = info->m_cy; info->csz = info->m_cz; vec_prod (tmp_x, tmp_y, tmp_z, info->bsx / (2 * PI), info->bsy / (2 * PI), info->bsz / (2 * PI), info->csx / (2 * PI), info->csy / (2 * PI), info->csz / (2 * PI)); info->V0 = scalar_prod (info->asx / (2 * PI), info->asy / (2 * PI), info->asz / (2 * PI), tmp_x, tmp_y, tmp_z); if (!info->V0 || isnan (info->V0)) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Lattice reciprocal cell volume is found invalid V0=%g. Check crystal cell parameters.\n", info->compname, info->V0)); return (0); } info->V0 = 1 / fabs (info->V0); /*compute the direct cell parameters, for completeness*/ info->m_ax = tmp_x * info->V0; info->m_ay = tmp_y * info->V0; info->m_az = tmp_z * info->V0; vec_prod (tmp_x, tmp_y, tmp_z, info->csx / (2 * PI), info->csy / (2 * PI), info->csz / (2 * PI), info->asx / (2 * PI), info->asy / (2 * PI), info->asz / (2 * PI)); info->m_bx = tmp_x * info->V0; info->m_by = tmp_y * info->V0; info->m_bz = tmp_z * info->V0; vec_prod (tmp_x, tmp_y, tmp_z, info->asx / (2 * PI), info->asy / (2 * PI), info->asz / (2 * PI), info->bsx / (2 * PI), info->bsy / (2 * PI), info->bsz / (2 * PI)); info->m_cx = tmp_x * info->V0; info->m_cy = tmp_y * info->V0; info->m_cz = tmp_z * info->V0; } MPI_MASTER (printf ("V0=%g\n", info->V0);); /* allocate hkl_data array */ if (surf_size <= 0 || !surf_k || !surf_FT2 || !surf_dir) surf_size = 1; list = (struct hkl_data*)malloc (size * surf_size * sizeof (struct hkl_data)); if (!list) { fprintf (stderr, "Single_crystal: Error allocating reflection list\n"); return (0); } for (i = 0; i < size; i++) for (j = 0; j < surf_size; j++) { double h = 0, k = 0, l = 0, F2 = 0; double b1[3], b2[3]; double sig1, sig2, sig3; /* get data from table */ h = Table_Index (sTable, i, info->column_order[0] - 1); k = Table_Index (sTable, i, info->column_order[1] - 1); l = Table_Index (sTable, i, info->column_order[2] - 1); if (info->column_order[3]) { F2 = Table_Index (sTable, i, info->column_order[3] - 1); F2 *= F2; } else if (info->column_order[4]) F2 = Table_Index (sTable, i, info->column_order[4] - 1); list[J].h = h; list[J].k = k; list[J].l = l; list[J].F2 = F2; sum_F2 += F2; /* Precompute some values */ list[J].tau_x = h * info->asx + k * info->bsx + l * info->csx; list[J].tau_y = h * info->asy + k * info->bsy + l * info->csy; list[J].tau_z = h * info->asz + k * info->bsz + l * info->csz; // create satellite peaks from e.g. surface truncation (CTR) if (surf_size > 1 && surf_FT2 && surf_k && surf_FT2[j] > 0 && surf_k[j]) { list[J].tau_x += surf_dir[0] * surf_k[j]; list[J].tau_y += surf_dir[1] * surf_k[j]; list[J].tau_z += surf_dir[2] * surf_k[j]; list[J].F2 *= surf_FT2[j]; } list[J].tau = sqrt (list[J].tau_x * list[J].tau_x + list[J].tau_y * list[J].tau_y + list[J].tau_z * list[J].tau_z); list[J].u1x = list[J].tau_x / list[J].tau; list[J].u1y = list[J].tau_y / list[J].tau; list[J].u1z = list[J].tau_z / list[J].tau; sig1 = FWHM2RMS * info->m_delta_d_d * list[J].tau; /* Find two arbitrary axes perpendicular to tau and each other. */ normal_vec (&b1[0], &b1[1], &b1[2], list[J].u1x, list[J].u1y, list[J].u1z); vec_prod (b2[0], b2[1], b2[2], list[J].u1x, list[J].u1y, list[J].u1z, b1[0], b1[1], b1[2]); /* Find the two mosaic axes perpendicular to tau. */ if (SC_mosaic > 0) { /* Use isotropic mosaic. */ list[J].u2x = b1[0]; list[J].u2y = b1[1]; list[J].u2z = b1[2]; sig2 = FWHM2RMS * list[J].tau * MIN2RAD * SC_mosaic; list[J].u3x = b2[0]; list[J].u3y = b2[1]; list[J].u3z = b2[2]; sig3 = FWHM2RMS * list[J].tau * MIN2RAD * SC_mosaic; } else if (SC_mosaic_a > 0 && SC_mosaic_b > 0 && SC_mosaic_c > 0) { /* Use anisotropic mosaic. */ MPI_MASTER (fprintf (stderr, "Single_crystal: WARNING: you are using an experimental feature:\n" " anistropic mosaicity. Please examine your data carefully.\n");); /* compute the jacobian of (tau_v,tau_n) from rotations around the unit cell vectors. */ struct hkl_data* l = &(list[J]); double xia_x, xia_y, xia_z, xib_x, xib_y, xib_z, xic_x, xic_y, xic_z; /*input parameters are in arc minutes*/ double sig_fi_a = SC_mosaic_a * MIN2RAD; double sig_fi_b = SC_mosaic_b * MIN2RAD; double sig_fi_c = SC_mosaic_c * MIN2RAD; if (info->m_a == 0) info->m_a = sqrt (scalar_prod (info->m_ax, info->m_ay, info->m_az, info->m_ax, info->m_ay, info->m_az)); if (info->m_b == 0) info->m_b = sqrt (scalar_prod (info->m_bx, info->m_by, info->m_bz, info->m_bx, info->m_by, info->m_bz)); if (info->m_c == 0) info->m_c = sqrt (scalar_prod (info->m_cx, info->m_cy, info->m_cz, info->m_cx, info->m_cy, info->m_cz)); l->u2x = b1[0]; l->u2y = b1[1]; l->u2z = b1[2]; l->u3x = b2[0]; l->u3y = b2[1]; l->u3z = b2[2]; xia_x = l->tau_x - (M_2_PI * h / info->m_a) * info->asx; xia_y = l->tau_y - (M_2_PI * h / info->m_a) * info->asy; xia_z = l->tau_z - (M_2_PI * h / info->m_a) * info->asz; xib_x = l->tau_x - (M_2_PI * h / info->m_b) * info->bsx; xib_y = l->tau_y - (M_2_PI * h / info->m_b) * info->bsy; xib_z = l->tau_z - (M_2_PI * h / info->m_b) * info->bsz; xic_x = l->tau_x - (M_2_PI * h / info->m_c) * info->csx; xic_y = l->tau_y - (M_2_PI * h / info->m_c) * info->csy; xic_z = l->tau_z - (M_2_PI * h / info->m_c) * info->csz; double xia = sqrt (xia_x * xia_x + xia_y * xia_y + xia_z * xia_z); double xib = sqrt (xib_x * xib_x + xib_y * xib_y + xib_z * xib_z); double xic = sqrt (xic_x * xic_x + xic_y * xic_y + xic_z * xic_z); vec_prod (tmp_x, tmp_y, tmp_z, l->tau_x, l->tau_y, l->tau_z, l->u2x, l->u2y, l->u2z); double J_n_fia = xia / info->m_a / l->tau * scalar_prod (info->asx, info->asy, info->asz, tmp_x, tmp_y, tmp_z); vec_prod (tmp_x, tmp_y, tmp_z, l->tau_x, l->tau_y, l->tau_z, l->u2x, l->u2y, l->u2z); double J_n_fib = xib / info->m_b / l->tau * scalar_prod (info->bsx, info->bsy, info->bsz, tmp_x, tmp_y, tmp_z); vec_prod (tmp_x, tmp_y, tmp_z, l->tau_x, l->tau_y, l->tau_z, l->u2x, l->u2y, l->u2z); double J_n_fic = xic / info->m_c / l->tau * scalar_prod (info->csx, info->csy, info->csz, tmp_x, tmp_y, tmp_z); vec_prod (tmp_x, tmp_y, tmp_z, l->tau_x, l->tau_y, l->tau_z, l->u3x, l->u3y, l->u3z); double J_v_fia = xia / info->m_a / l->tau * scalar_prod (info->asx, info->asy, info->asz, tmp_x, tmp_y, tmp_z); vec_prod (tmp_x, tmp_y, tmp_z, l->tau_x, l->tau_y, l->tau_z, l->u3x, l->u3y, l->u3z); double J_v_fib = xib / info->m_b / l->tau * scalar_prod (info->bsx, info->bsy, info->bsz, tmp_x, tmp_y, tmp_z); vec_prod (tmp_x, tmp_y, tmp_z, l->tau_x, l->tau_y, l->tau_z, l->u3x, l->u3y, l->u3z); double J_v_fic = xic / info->m_c / l->tau * scalar_prod (info->csx, info->csy, info->csz, tmp_x, tmp_y, tmp_z); /*with the jacobian we can compute the sigmas in terms of the orthogonal vectors u2 and u3*/ sig2 = sig_fi_a * fabs (J_v_fia) + sig_fi_b * fabs (J_v_fib) + sig_fi_c * fabs (J_v_fic); sig3 = sig_fi_a * fabs (J_n_fia) + sig_fi_b * fabs (J_n_fib) + sig_fi_c * fabs (J_n_fic); } else if (SC_mosaic_AB[0] != 0 && SC_mosaic_AB[1] != 0) { /* Use isotropic mosaic via mosaic matrix */ if ((SC_mosaic_AB[2] == 0 && SC_mosaic_AB[3] == 0 && SC_mosaic_AB[4] == 0) || (SC_mosaic_AB[5] == 0 && SC_mosaic_AB[6] == 0 && SC_mosaic_AB[7] == 0)) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: in-plane mosaics are specified but one (or both)\n" " in-plane reciprocal vector is the zero vector.\n", info->compname);); return (0); } MPI_MASTER (fprintf (stderr, "Single_crystal: %s: WARNING: you are using an experimental feature: \n" " \"in-plane\" anistropic mosaicity. Please examine your data carefully.\n", info->compname);); /*for given reflection in list - compute linear comb of tau_a and tau_b*/ /*check for not in plane - f.i. check if (tau_a X tau_b).tau_i)==0*/ struct hkl_data* l = &(list[J]); double det, c1, c2, sig_tau_c; double em_x, em_y, em_z, tmp_x, tmp_y, tmp_z; double tau_a[3], tau_b[3]; /*convert Miller indices to taus*/ if (info->m_a == 0) info->m_a = sqrt (scalar_prod (info->m_ax, info->m_ay, info->m_az, info->m_ax, info->m_ay, info->m_az)); if (info->m_b == 0) info->m_b = sqrt (scalar_prod (info->m_bx, info->m_by, info->m_bz, info->m_bx, info->m_by, info->m_bz)); if (info->m_c == 0) info->m_c = sqrt (scalar_prod (info->m_cx, info->m_cy, info->m_cz, info->m_cx, info->m_cy, info->m_cz)); tau_a[0] = M_2_PI * ((SC_mosaic_AB[2] / info->m_a) * info->asx + (SC_mosaic_AB[3] / info->m_b) * info->bsx + (SC_mosaic_AB[4] / info->m_c) * info->csx); tau_a[1] = M_2_PI * ((SC_mosaic_AB[2] / info->m_a) * info->asy + (SC_mosaic_AB[3] / info->m_b) * info->bsy + (SC_mosaic_AB[4] / info->m_c) * info->csy); tau_a[2] = M_2_PI * ((SC_mosaic_AB[2] / info->m_a) * info->asz + (SC_mosaic_AB[3] / info->m_b) * info->bsz + (SC_mosaic_AB[4] / info->m_c) * info->csz); tau_b[0] = M_2_PI * ((SC_mosaic_AB[5] / info->m_a) * info->asx + (SC_mosaic_AB[6] / info->m_b) * info->bsx + (SC_mosaic_AB[7] / info->m_c) * info->csx); tau_b[1] = M_2_PI * ((SC_mosaic_AB[5] / info->m_a) * info->asy + (SC_mosaic_AB[6] / info->m_b) * info->bsy + (SC_mosaic_AB[7] / info->m_c) * info->csy); tau_b[2] = M_2_PI * ((SC_mosaic_AB[5] / info->m_a) * info->asz + (SC_mosaic_AB[6] / info->m_b) * info->bsz + (SC_mosaic_AB[7] / info->m_c) * info->csz); /*check determinants to see how we should compute the linear combination of a and b (to match c)*/ c1 = c2 = 0; if ((det = tau_a[0] * tau_b[1] - tau_a[1] * tau_b[0]) != 0) { c1 = (l->tau_x * tau_b[1] - l->tau_y * tau_b[0]) / det; c2 = (tau_a[0] * l->tau_y - tau_a[1] * l->tau_x) / det; } else if ((det = tau_a[1] * tau_b[2] - tau_a[2] * tau_b[1]) != 0) { c1 = (l->tau_y * tau_b[2] - l->tau_z * tau_b[1]) / det; c2 = (tau_a[1] * l->tau_z - tau_a[2] * l->tau_y) / det; } else if ((det = tau_a[0] * tau_b[2] - tau_a[2] * tau_b[0]) != 0) { c1 = (l->tau_x * tau_b[2] - l->tau_z * tau_b[0]) / det; c2 = (tau_a[0] * l->tau_z - tau_a[2] * l->tau_x) / det; } if ((c1 == 0) && (c2 == 0)) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: WARNING: reflection tau[%i]=(%g %g %g) " "has no component in defined mosaic plane\n", info->compname, i, l->tau_x, l->tau_y, l->tau_z);); } /*compute linear combination => sig_tau_i = | c1*sig_tau_a + c2*sig_tau_b | - also add in the minute to radian scaling factor*/; sig_tau_c = MIN2RAD * sqrt (c1 * SC_mosaic_AB[0] * c1 * SC_mosaic_AB[0] + c2 * SC_mosaic_AB[1] * c2 * SC_mosaic_AB[1]); l->u2x = b1[0]; l->u2y = b1[1]; l->u2z = b1[2]; l->u3x = b2[0]; l->u3y = b2[1]; l->u3z = b2[2]; /*so now let's compute the rotation around planenormal tau_a X tau_b*/ /*g_bar (unit normal of rotation plane) = tau_a X tau_b / norm(tau_a X tau_b)*/ vec_prod (tmp_x, tmp_y, tmp_z, tau_a[0], tau_a[1], tau_a[2], tau_b[0], tau_b[1], tau_b[2]); vec_prod (em_x, em_y, em_z, l->tau_x, l->tau_y, l->tau_z, tmp_x, tmp_y, tmp_z); NORM (em_x, em_y, em_z); sig2 = l->tau * sig_tau_c * fabs (scalar_prod (em_x, em_y, em_z, l->u2x, l->u2y, l->u2z)); sig3 = l->tau * sig_tau_c * fabs (scalar_prod (em_x, em_y, em_z, l->u3x, l->u3y, l->u3z)); /*protect against collapsing gaussians. These seem to be sensible values.*/ if (sig2 < 1e-5) sig2 = 1e-5; if (sig3 < 1e-5) sig3 = 1e-5; } /* end moisaic stuff */ else { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: EITHER mosaic OR (mosaic_a, mosaic_b, mosaic_c)\n" " must be given and be >0.\n", info->compname);); return (0); } list[J].sig123 = sig1 * sig2 * sig3; list[J].m1 = 1 / (2 * sig1 * sig1); list[J].m2 = 1 / (2 * sig2 * sig2); list[J].m3 = 1 / (2 * sig3 * sig3); /* Set Gauss cutoff to 5 times the maximal sigma. */ if (sig1 > sig2) list[J].cutoff = 5 * (sig1 > sig3 ? sig1 : sig3); else list[J].cutoff = 5 * (sig2 > sig3 ? sig2 : sig3); if (!isnan (list[J].F2) && !isnan (list[J].cutoff) && !isnan (list[J].tau)) J++; // switch to next HKL Bragg element in hkl_data } // for i, j Table_Free (&sTable); if (!sum_F2) { MPI_MASTER (exit (fprintf (stderr, "Single_crystal: %s: ERROR: all %i structure factors in file '%s' are null. Check the reflection list.\n", info->compname, i, SC_file));); } /* sort the list with increasing tau */ qsort (list, J, sizeof (struct hkl_data), SX_list_compare); *hkl_list = list; info->count = J; return (info->count); } /* read_hkl_data */ /* ------------------------------------------------------------------------ */ /* hkl_search search the HKL reflections which are on the Ewald sphere input: L,T,count,V0: constants for all calls kix,kiy,kiz: may be different for each call this function returns: tau_count (return), coh_refl, coh_xsect, T (updated elements in the array up to [j]) */ #pragma acc routine int hkl_search (struct hkl_data* L, void* TT, int count, double V0, double kix, double kiy, double kiz, double tau_max, double* coh_refl, double* coh_xsect) { int i, j = 0; double ki = sqrt (kix * kix + kiy * kiy + kiz * kiz); struct tau_data* T = (struct tau_data*)TT; /* Common factor in coherent cross-section */ double xsect_factor = pow (2 * PI, 5.0 / 2.0) / (V0 * ki * ki); for (i = 0; i < count; i++) { /* Assuming reflections are sorted, stop search when max tau exceeded. */ double rho, rho_x, rho_y, rho_z; double diff; if (L[i].tau > tau_max) break; // passed max tau found in crystal reflections /* Check if this reciprocal lattice point is close enough to the Ewald sphere to make scattering possible. */ rho_x = kix - L[i].tau_x; /* rho = |ki - tau| = |kf| = |ki| within cutoff on the Ewald sphere */ rho_y = kiy - L[i].tau_y; rho_z = kiz - L[i].tau_z; rho = sqrt (rho_x * rho_x + rho_y * rho_y + rho_z * rho_z); diff = fabs (rho - ki); /* Check if scattering is possible (cutoff of Gaussian tails). */ if (diff <= L[i].cutoff) { double ox, oy, oz; double b1x, b1y, b1z, b2x, b2y, b2z, kx, ky, kz, nx, ny, nz; double n11, n22, n12, det_N, inv_n11, inv_n22, inv_n12, l11, l22, l12, det_L; double Bt_D_O_x, Bt_D_O_y, y0x, y0y, alpha; /* Store reflection. */ T[j].index = i; /* Get ki vector in local coordinates. */ kx = kix * L[i].u1x + kiy * L[i].u1y + kiz * L[i].u1z; ky = kix * L[i].u2x + kiy * L[i].u2y + kiz * L[i].u2z; kz = kix * L[i].u3x + kiy * L[i].u3y + kiz * L[i].u3z; T[j].rho_x = kx - L[i].tau; T[j].rho_y = ky; T[j].rho_z = kz; T[j].rho = rho; /* Compute the tangent plane of the Ewald sphere. */ nx = T[j].rho_x / T[j].rho; ny = T[j].rho_y / T[j].rho; nz = T[j].rho_z / T[j].rho; ox = (ki - T[j].rho) * nx; oy = (ki - T[j].rho) * ny; oz = (ki - T[j].rho) * nz; T[j].ox = ox; T[j].oy = oy; T[j].oz = oz; /* Compute unit vectors b1 and b2 that span the tangent plane. */ normal_vec (&b1x, &b1y, &b1z, nx, ny, nz); vec_prod (b2x, b2y, b2z, nx, ny, nz, b1x, b1y, b1z); T[j].b1x = b1x; T[j].b1y = b1y; T[j].b1z = b1z; T[j].b2x = b2x; T[j].b2y = b2y; T[j].b2z = b2z; /* Compute the 2D projection of the 3D Gauss of the reflection. */ /* The symmetric 2x2 matrix N describing the 2D gauss. */ n11 = L[i].m1 * b1x * b1x + L[i].m2 * b1y * b1y + L[i].m3 * b1z * b1z; n12 = L[i].m1 * b1x * b2x + L[i].m2 * b1y * b2y + L[i].m3 * b1z * b2z; n22 = L[i].m1 * b2x * b2x + L[i].m2 * b2y * b2y + L[i].m3 * b2z * b2z; /* The (symmetric) inverse matrix of N. */ det_N = n11 * n22 - n12 * n12; inv_n11 = n22 / det_N; inv_n12 = -n12 / det_N; inv_n22 = n11 / det_N; /* The Cholesky decomposition of 1/2*inv_n (lower triangular L). */ l11 = sqrt (inv_n11 / 2); l12 = inv_n12 / (2 * l11); l22 = sqrt (inv_n22 / 2 - l12 * l12); T[j].l11 = l11; T[j].l12 = l12; T[j].l22 = l22; det_L = l11 * l22; /* The product B^T D o. */ Bt_D_O_x = b1x * L[i].m1 * ox + b1y * L[i].m2 * oy + b1z * L[i].m3 * oz; Bt_D_O_y = b2x * L[i].m1 * ox + b2y * L[i].m2 * oy + b2z * L[i].m3 * oz; /* Center of 2D Gauss in plane coordinates. */ y0x = -(Bt_D_O_x * inv_n11 + Bt_D_O_y * inv_n12); y0y = -(Bt_D_O_x * inv_n12 + Bt_D_O_y * inv_n22); T[j].y0x = y0x; T[j].y0y = y0y; /* Factor alpha for the distance of the 2D Gauss from the origin. */ alpha = L[i].m1 * ox * ox + L[i].m2 * oy * oy + L[i].m3 * oz * oz - (y0x * y0x * n11 + y0y * y0y * n22 + 2 * y0x * y0y * n12); T[j].refl = xsect_factor * det_L * exp (-alpha) / L[i].sig123; /* intensity of that Bragg */ *coh_refl += T[j].refl; /* total scatterable intensity */ T[j].xsect = T[j].refl * L[i].F2; *coh_xsect += T[j].xsect; j++; } /*protect against tau shortlist buffer overrrun*/ if (j == MCSX_REFL_SLIST_SIZE) { break; } } /* end for */ return (j); // this is 'tau_count', i.e. number of reachable reflections } /* end hkl_search */ #pragma acc routine int hkl_select (struct tau_data* T, int tau_count, double coh_refl, double* sum, _class_particle* _particle) { int j; double r = rand0max (coh_refl); *sum = 0; for (j = 0; j < tau_count; j++) { *sum += T[j].refl; if (*sum > r) break; } return j; } /* Functions for "reorientation", powder mode */ /* Powder, forward */ #pragma acc routine void randrotate (double* nx, double* ny, double* nz, double a, double b, double c) { double x1, y1, z1, x2, y2, z2; rotate (x1, y1, z1, *nx, *ny, *nz, a, 1, 0, 0); /* <1> = rot(,a) */ rotate (x2, y2, z2, x1, y1, z1, b, 0, 1, 0); /* <2> = rot(<1>,b) */ rotate (*nx, *ny, *nz, x2, y2, z2, c, 0, 0, 1); /* = rot(<2>,c) */ } /* Powder, back */ #pragma acc routine void randderotate (double* nx, double* ny, double* nz, double a, double b, double c) { double x1, y1, z1, x2, y2, z2; rotate (x1, y1, z1, *nx, *ny, *nz, -c, 0, 0, 1); rotate (x2, y2, z2, x1, y1, z1, -b, 0, 1, 0); rotate (*nx, *ny, *nz, x2, y2, z2, -a, 1, 0, 0); } #pragma acc routine /* rotate vector counterclockwise */ void vec_rotate_2d (double* x, double* y, double angle) { double c, s; double newx, newy; c = cos (angle); s = sin (angle); newx = *x * c - *y * s; newy = *x * s + *y * c; *x = newx; *y = newy; } struct sx_abs_data { int muc; // column where mu is to be found double atomic_weight; double density; // mass density [g/cm^3] double delta_prefactor; t_Table table; }; int sx_read_abs_data (char* ABS_file, struct sx_abs_data* abs, struct hkl_info_struct info) { int status; char** parsing; double Z; if (ABS_file && strlen (ABS_file) && strcmp (ABS_file, "NULL") && strcmp (ABS_file, "0")) { if ((status = Table_Read (&(abs->table), ABS_file, 0)) == -1) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: ERROR: Could not parse file \"%s\"\n", info.compname, ABS_file);); exit (-1); } abs->muc = (abs->table.columns == 3 ? 1 : 5); abs->density = 0; parsing = Table_ParseHeader (abs->table.header, "Z", "A[r]", "rho", NULL); if (parsing) { if (parsing[0]) Z = strtol (parsing[0], NULL, 10); if (parsing[1] && !info.mat_weight) info.mat_weight = strtod (parsing[1], NULL); // aka A[r] in g/mol if (parsing[2] && !info.rho) info.rho = strtod (parsing[2], NULL); if (parsing[2]) abs->density = strtod (parsing[2], NULL); // nominal density of the absorption material } if (info.mat_density > 0) abs->density = info.mat_density; // crystal density from the reflections file wins if (info.rho && info.mat_weight) abs->delta_prefactor = NA * (info.rho * 1e-24) / info.mat_weight * 2.0 * M_PI * RE; abs->atomic_weight = info.mat_weight; return 1; } else { // init an empty table with absorption 0 Table_Init (&(abs->table), 2, 2); abs->table.data[0] = 0; abs->table.data[1] = 0; abs->table.data[2] = FLT_MAX; abs->table.data[3] = 0; MPI_MASTER (fprintf (stderr, "Single_crystal: %s: WARNING: material file %s not found. Absorption set to 0.\n", info.compname, ABS_file)); return 0; } } // sx_read_abs_data #endif /* !SINGLE_CRYSTAL_DECL */ %} DECLARE %{ struct hkl_info_struct hkl_info; off_struct offdata; struct hkl_data* hkl_list; struct tau_data tau_list[MCSX_REFL_SLIST_SIZE]; struct sx_abs_data abs_info; %} INITIALIZE %{ double as, bs, cs; int i = 0; /* transfer input parameters */ strncpy (hkl_info.compname, NAME_CURRENT_COMP, 256); hkl_info.m_delta_d_d = delta_d_d; hkl_info.m_a = 0; hkl_info.m_b = 0; hkl_info.m_c = 0; hkl_info.m_aa = aa; hkl_info.m_bb = bb; hkl_info.m_cc = cc; hkl_info.m_ax = ax; hkl_info.m_ay = ay; hkl_info.m_az = az; hkl_info.m_bx = bx; hkl_info.m_by = by; hkl_info.m_bz = bz; hkl_info.m_cx = cx; hkl_info.m_cy = cy; hkl_info.m_cz = cz; hkl_info.sigma_i = sigma_inc; hkl_info.recip = recip_cell; /* default format h,k,l,F,F2 */ hkl_info.column_order[0] = 1; hkl_info.column_order[1] = 2; hkl_info.column_order[2] = 3; hkl_info.column_order[3] = 0; hkl_info.column_order[4] = 7; hkl_info.kix = hkl_info.kiy = hkl_info.kiz = 0; hkl_info.nb_reuses = hkl_info.nb_refl = hkl_info.nb_refl_count = 0; hkl_info.tau_count = 0; hkl_info.flag_barns = barns; // initial tests/warnings if (powder && !(order == 1)) { MPI_MASTER (fprintf (stderr, "Single_crystal: %s: powder mode means implicit choice of no multiple scattering!\n" "WARNING setting order=1\n", NAME_CURRENT_COMP);); order = 1; } if (powder && (RX || RY)) { exit (fprintf (stderr, "Single_crystal: %s: powder mode can not be used together with crystal curvature.\n", NAME_CURRENT_COMP)); } /*this should not be in hkl_info*/ hkl_info.shape = -1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape */ if (geometry && strlen (geometry) && strcmp (geometry, "NULL") && strcmp (geometry, "0")) { #ifndef USE_OFF MPI_MASTER ( fprintf (stderr, "Single_crystal: %s: ERROR: You are attempting to use an OFF/PLY geometry without -DUSE_OFF. You will need to recompile with that define set!\n", NAME_CURRENT_COMP);); exit (-1); #else if (off_init (geometry, xwidth, yheight, zdepth, 0, &offdata)) { hkl_info.shape = 3; } #endif } else if (xwidth && yheight && zdepth) hkl_info.shape = 1; /* box */ else if (radius > 0 && yheight) hkl_info.shape = 0; /* cylinder */ else if (radius > 0 && !yheight) hkl_info.shape = 2; /* sphere */ if (hkl_info.shape < 0) MPI_MASTER (fprintf (stderr, "Single_crystal: %s: sample has invalid dimensions.\n" "WARNING Please check parameter values (xwidth, yheight, zdepth, radius).\n" " Component will be ignored.\n", NAME_CURRENT_COMP);); /* determine if we have a surface broadening/truncation rods (CTR) */ hkl_info.ctr_size = surf_size; hkl_info.ctr_k = surf_k; hkl_info.ctr_FT2 = surf_FT2; hkl_info.ctr_dir = surf_dir; if (hkl_info.ctr_size < 0) { /* CTR automatic mode */ double surf_dir_sum = 0; double surf_thickness = 0; for (i = 0; i < 3; hkl_info.ctr_dir[i++] = 0) ; if (hkl_info.shape == 1) { if (xwidth <= yheight && xwidth <= zdepth) { hkl_info.ctr_dir[0] = 1; surf_thickness = xwidth; } else if (yheight <= xwidth && yheight <= zdepth) { hkl_info.ctr_dir[1] = 1; surf_thickness = yheight; } else { hkl_info.ctr_dir[2] = 1; surf_thickness = zdepth; } } else if (hkl_info.shape == 0 && yheight <= radius) { hkl_info.ctr_dir[1] = 1; surf_thickness = yheight; } // surf_size=-1: set automatically the k and FT2 vectors if (0 < surf_thickness) { hkl_info.ctr_size = 21; // nb of satellite peaks double range_k = 1.0; // k range [1/Angs], k=[-range/2 : +range/2] double dk = range_k / (hkl_info.ctr_size - 1); // delta-k for btw satellites DArray1d s_k = create_darr1d (hkl_info.ctr_size); DArray1d s_F = create_darr1d (hkl_info.ctr_size); hkl_info.ctr_k = (double*)s_k; hkl_info.ctr_FT2 = (double*)s_F; for (i = 0; i < hkl_info.ctr_size; i++) { hkl_info.ctr_k[i] = -range_k / 2 + i * dk; // 1/Angs, symmetric around 0 [1/Angs] if (surf_thickness < 1e-7) { // thin layer thickness < 0.1 um -> sinc^2 double x = PI * hkl_info.ctr_k[i] * surf_thickness * 1e10; // 1e10 [m -> Angs] hkl_info.ctr_FT2[i] = hkl_info.ctr_k[i] ? sin (x) / x : 1; } else { // assume a half-bulk -> 1/k^2 [in rlu] hkl_info.ctr_FT2[i] = hkl_info.ctr_k[i] ? 1 / (hkl_info.ctr_k[i] * 2 * PI) : 4 * PI; // we assume a=2PI } hkl_info.ctr_FT2[i] *= hkl_info.ctr_FT2[i]; // squared } MPI_MASTER (printf ("Single_crystal: %s: Info: Setting CTR surface normal XYZ=[%g,%g,%g], as a %s.\n", NAME_CURRENT_COMP, hkl_info.ctr_dir[0], hkl_info.ctr_dir[1], hkl_info.ctr_dir[2], surf_thickness < 1e-7 ? "thin layer (sinc^2)" : "half-bulk (1/k^2)");); } } if (hkl_info.ctr_size > 1 && (!hkl_info.ctr_dir || !hkl_info.ctr_dir[0] && !hkl_info.ctr_dir[1] && !hkl_info.ctr_dir[2])) { // no CTR normal ? MPI_MASTER (printf ("Single_crystal: %s: Warning: No surf_dir surface normal set. Ignoring CTR.\n", i, NAME_CURRENT_COMP);); hkl_info.ctr_size = 1; } if (hkl_info.ctr_size > 1 && hkl_info.ctr_FT2 && hkl_info.ctr_k) { // do some CTR tests double sum_FT2 = 0, min_k = FLT_MAX, max_k = -FLT_MAX; for (i = 0; i < hkl_info.ctr_size; i++) { if (hkl_info.ctr_FT2[i] < 0) { MPI_MASTER (printf ("Single_crystal: %s: Warning: surf_FT2[%i] < 0. Must all be >= 0. Ignoring CTR.\n", i, NAME_CURRENT_COMP);); hkl_info.ctr_size = 1; break; } sum_FT2 += hkl_info.ctr_FT2[i]; if (hkl_info.ctr_k[i] < min_k) min_k = hkl_info.ctr_k[i]; if (max_k < hkl_info.ctr_k[i]) max_k = hkl_info.ctr_k[i]; } if (hkl_info.ctr_size > 1 && min_k * max_k >= 0) { // k vector not around 0 ? MPI_MASTER (printf ("Single_crystal: %s: Warning: surf_k = [%g:%g]. Must be around 0. Ignoring CTR.\n", min_k, max_k, NAME_CURRENT_COMP);); hkl_info.ctr_size = 1; } if (hkl_info.ctr_size > 1 && sum_FT2 <= 0) { // no positive FT2 value ? MPI_MASTER (printf ("Single_crystal: %s: Warning: sum(surf_FT2) <= 0. Must be > 0. Ignoring CTR.\n", min_k, max_k, NAME_CURRENT_COMP);); hkl_info.ctr_size = 1; } if (hkl_info.ctr_size > 1 && sum_FT2) { // all is fine. Normalize CTR. NORM (hkl_info.ctr_dir[0], hkl_info.ctr_dir[1], hkl_info.ctr_dir[2]); // normal vector |u|=1 // normalize surf_FT2 (distribute Bragg scattering along all satellites) for (i = 0; i < hkl_info.ctr_size; hkl_info.ctr_FT2[i++] /= sum_FT2) ; } } // if (CTR) /* ought to be cleaned up as mosaic_AB now is a proper vector/array and not a define */ double* mosaic_ABin = mosaic_AB; /* Read in structure factors, and do some pre-calculations. */ if (!read_hkl_data (reflections, &hkl_info, &hkl_list, mosaic, mosaic_a, mosaic_b, mosaic_c, mosaic_ABin, hkl_info.ctr_size, hkl_info.ctr_k, hkl_info.ctr_FT2, hkl_info.ctr_dir)) { MPI_MASTER (printf ("Single_crystal: %s: Can not set crystal structure.\n" "WARNING Check file and mosaic. Ignoring.\n", NAME_CURRENT_COMP);); hkl_info.shape = -1; // ignore component. The failed read_hkl_data has set hkl_info.count = 0 } if (!sx_read_abs_data (material_datafile, &abs_info, hkl_info)) abs_info.muc = 0; if (hkl_info.sigma_i < 0) hkl_info.sigma_i = 0; // final information MPI_MASTER ( if (hkl_info.count) { printf ("Single_crystal: %s: Read %d reflections from file '%s'\n", NAME_CURRENT_COMP, hkl_info.count, reflections); } printf ("Single_crystal: %s: Vc=%g [Angs^3] sigma_inc=%g [barn] reflections=%s\n", NAME_CURRENT_COMP, hkl_info.V0, hkl_info.sigma_i, reflections&& strlen (reflections) ? reflections : "NULL"); printf ("Direct space lattice orientation:\n"); printf (" a = [%g %g %g]\n", hkl_info.m_ax, hkl_info.m_ay, hkl_info.m_az); printf (" b = [%g %g %g]\n", hkl_info.m_bx, hkl_info.m_by, hkl_info.m_bz); printf (" c = [%g %g %g]\n", hkl_info.m_cx, hkl_info.m_cy, hkl_info.m_cz); printf ("Reciprocal space lattice orientation:\n"); printf (" a* = [%g %g %g]\n", hkl_info.asx, hkl_info.asy, hkl_info.asz); printf (" b* = [%g %g %g]\n", hkl_info.bsx, hkl_info.bsy, hkl_info.bsz); printf (" c* = [%g %g %g]\n", hkl_info.csx, hkl_info.csy, hkl_info.csz); ); %} TRACE %{ double l1, l2 = -1; /* Entry and exit lengths in sample */ struct hkl_data* L; /* Structure factor list */ int i; /* Index into structure factor list */ #ifndef OPENACC struct tau_data* T; /* List of reflections close to Ewald sphere */ #else struct tau_data T[MCSX_REFL_SLIST_SIZE]; #endif int tau_count; /* Number of reflections close to Ewald sphere*/ int j; /* Index into reflection list */ int event_counter; /* scattering event counter */ double kix, kiy, kiz, ki; /* Initial wave vector [1/AA] */ double kfx, kfy, kfz; /* Final wave vector */ double k; /* photon wave vector */ double rho_x, rho_y, rho_z; /* the vector ki - tau */ double rho; double diff; /* Deviation from Bragg condition */ double ox, oy, oz; /* Origin of Ewald sphere tangent plane */ double b1x, b1y, b1z; /* First vector spanning tangent plane */ double b2x, b2y, b2z; /* Second vector spanning tangent plane */ double n11, n12, n22; /* 2D Gauss description matrix N */ double det_N; /* Determinant of N */ double inv_n11, inv_n12, inv_n22; /* Inverse of N */ double l11, l12, l22; /* Cholesky decomposition L of 1/2*inv(N) */ double det_L; /* Determinant of L */ double Bt_D_O_x, Bt_D_O_y; /* Temporaries */ double y0x, y0y; /* Center of 2D Gauss in plane coordinates */ double alpha; /* Offset of 2D Gauss center from 3D center */ double V0; /* Volume of unit cell */ double l_full; /* photon path length for transmission */ double l; /* Path length to scattering event */ double abs_xlen; /* Absorption length */ double inc_xsect, inc_xlen; /* Incoherent scattering cross section and length */ double coh_xsect, coh_xlen; /* Coherent cross section and length */ double tot_xsect, tot_xlen; /* Total cross section and length */ double z1, z2, y1, y2; /* Temporaries to choose kf from 2D Gauss */ double adjust, sum; /* Temporaries */ double p_trans; /* Transmission probability */ double mc_trans, mc_interact; /* Transmission, interaction MC choices */ int intersect = 0; double curv_xangle; double curv_yangle; double _kx; double _ky; double _kz; char type; /* type of last event: t=transmit,c=coherent or i=incoherent */ int itype; /* type of last event: t=1,c=2 or i=3 */ #ifdef OPENACC #ifdef USE_OFF off_struct thread_offdata = offdata; #endif #else #define thread_offdata offdata #endif /* Intersection photon trajectory / sample (sample surface) */ if (hkl_info.shape == 0) intersect = cylinder_intersect (&l1, &l2, x, y, z, kx, ky, kz, radius, yheight); else if (hkl_info.shape == 1) intersect = box_intersect (&l1, &l2, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); else if (hkl_info.shape == 2) intersect = sphere_intersect (&l1, &l2, x, y, z, kx, ky, kz, radius); #ifdef USE_OFF else if (hkl_info.shape == 3) intersect = off_x_intersect (&l1, &l2, NULL, NULL, x, y, z, kx, ky, kz, thread_offdata); #endif if (l2 < 0) intersect = 0; /* we passed sample volume already */ if (intersect) { /* photon intersects crystal */ if (l1 > 0) PROP_DL (l1); /* Move to crystal surface if not inside */ ki = sqrt (kx * kx + ky * ky + kz * kz); event_counter = 0; V0 = hkl_info.V0; /* absorption: linear attenuation = mu [cm^2/g] * density [g/cm^3] = [1/cm], * 100 -> [1/m]. Without a density, fall back to one atom of atomic_weight per unit cell V0. */ abs_xlen = 0; if (abs_info.muc) { double mu = Table_Value (abs_info.table, ki * K2E, abs_info.muc); if (abs_info.density > 0) abs_xlen = mu * abs_info.density * 100; else abs_xlen = abs_info.atomic_weight / NA * mu * 1e26 / V0; // ponytail: per-atom weight over full cell, too weak for multi-atom cells; give a density } /*look into Compton scattering*/ inc_xsect = hkl_info.sigma_i; inc_xlen = inc_xsect / V0; if (hkl_info.flag_barns) { /*If cross sections are given in barns, we need a scaling factor of 100 to get scattering lengths in m, since V0 is assumed to be in AA*/ inc_xlen *= 100; } /* else assume fm^2 */ L = hkl_list; type = '\0'; itype = 0; #ifndef OPENACC T = tau_list; hkl_info.type = type; #endif do { /* Loop over multiple scattering events */ /* Angles for powder randomization */ double Alpha, Beta, Gamma; if (hkl_info.shape == 0) intersect = cylinder_intersect (&l1, &l2, x, y, z, kx, ky, kz, radius, yheight); else if (hkl_info.shape == 1) intersect = box_intersect (&l1, &l2, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); else if (hkl_info.shape == 2) intersect = sphere_intersect (&l1, &l2, x, y, z, kx, ky, kz, radius); #ifdef USE_OFF else if (hkl_info.shape == 3) intersect = off_x_intersect (&l1, &l2, NULL, NULL, x, y, z, kx, ky, kz, thread_offdata); #endif if (!intersect || l2 < -1e-7 || l1 > 1e-7) { /* photon is leaving the sample */ hkl_info.flag_warning++; if (hkl_info.flag_warning < 10) #ifndef OPENACC MPI_MASTER (fprintf (stderr, "Single_crystal: %s: WARNING: photon has unexpectedly left the crystal!\n" " l1=%g l2=%g x=%g y=%g z=%g kx=%g ky=%g kz=%g order=%i\n", NAME_CURRENT_COMP, l1, l2, x, y, z, kx, ky, kz, event_counter);); #endif break; } l_full = l2; if ((order && !(extra_order) && event_counter >= order) || (order && extra_order && event_counter >= order + extra_order)) { // Exit due to truncated order, weight with relevant cross-sections to distance l_full p *= exp (-abs_xlen * l_full); intersect = 0; break; } /* (1). Compute incoming wave vector ki */ if (powder) { /* orientation of crystallite is random */ Alpha = randpm1 () * PI * powder; Beta = randpm1 () * PI / 2; Gamma = randpm1 () * PI; randrotate (&kx, &ky, &kz, Alpha, Beta, Gamma); } /* ------------------------------------------------------------------------- */ /* lattice curvature option: rotate photon velocity */ curv_xangle = 0; curv_yangle = 0; _kx = kx; _ky = ky; _kz = kz; if (RY) { /* rotate v around x axis based on y pos, for vertical focus */ curv_yangle = atan2 (y, RY); vec_rotate_2d (&ky, &kz, curv_yangle); vec_rotate_2d (&Ey, &Ez, curv_yangle); /*changing y,z actually curves the crystal, not only the planes*/ /*comment out if only curvature of the lattice planes is needed*/ vec_rotate_2d (&y, &z, curv_yangle); } if (RX) { /* rotate v around y axis based on x pos, for horizontal focus */ curv_xangle = atan2 (x, RX); vec_rotate_2d (&kx, &kz, curv_xangle); vec_rotate_2d (&Ex, &Ez, curv_xangle); /*changing x,z actually curves the crystal, not only the planes*/ /*comment out if only curvature of the lattice planes is needed*/ vec_rotate_2d (&x, &z, curv_xangle); } kix = kx; kiy = ky; kiz = kz; kx = _kx; ky = _ky; kz = _kz; /* ------------------------------------------------------------------------- */ /* (2). Intersection of Ewald sphere with reciprocal lattice points */ double coh_xsect = 0, coh_refl = 0; // Condition to skip calculation of coherent cross section when, needed for extra_order feature if (order == 0 || extra_order == 0 || event_counter < order) { /* in case we use 'SPLIT' then consecutive photons can be identical when entering here and we may skip the hkl_search call */ #ifndef OPENACC if (order == 1 && fabs (kix - hkl_info.kix) < deltak && fabs (kiy - hkl_info.kiy) < deltak && fabs (kiz - hkl_info.kiz) < deltak) { hkl_info.nb_reuses++; /* Restore in case of matching event (e.g. SPLIT) */ coh_refl = hkl_info.coh_refl; coh_xsect = hkl_info.coh_xsect; tau_count = hkl_info.tau_count; } else { #endif /* Max possible tau for this ki with 5*sigma delta-d/d cutoff. */ double tau_max = 2 * ki / (1 - 5 * hkl_info.m_delta_d_d); /* call hkl_search */ tau_count = hkl_search (L, T, hkl_info.count, hkl_info.V0, kix, kiy, kiz, tau_max, &coh_refl, &coh_xsect); /* store ki so that we can check for further SPLIT iterations */ #ifndef OPENACC if (tau_count > hkl_info.max_tau_count) { hkl_info.max_tau_count = tau_count; } if (event_counter == 0) { /* only for incoming photon */ hkl_info.kix = kix; hkl_info.kiy = kiy; hkl_info.kiz = kiz; /* Store for potential re-use (e.g. SPLIT) */ hkl_info.coh_refl = coh_refl; hkl_info.coh_xsect = coh_xsect; hkl_info.tau_count = tau_count; hkl_info.nb_refl += tau_count; hkl_info.nb_refl_count++; } } #endif } else { // When extra_order used, disable coherent scattering after order reached, but continue // Set coherent cross section to zero to ignore coherent part coh_refl = 0; coh_xsect = 0; tau_count = 0; } /* (3). Probabilities of the different possible interactions. */ /* Cross-sections are in barns = 10**-28 m**2, and unit cell volumes are in AA**3 = 10**-30 m**2. Hence a factor of 100 is used to convert scattering lengths to m**-1 */ coh_xlen = coh_xsect / V0; if (hkl_info.flag_barns) { coh_xlen *= 100; } /* else assume fm^2 */ tot_xlen = abs_xlen + inc_xlen + coh_xlen; if (tot_xlen <= 0) { ABSORB; // Should we really absorb here? If "nothing" can happen we perhaps ought to "pass" instead? } /* (5). Transmission */ p_trans = exp (-tot_xlen * l_full); if (!event_counter && p_transmit >= 0 && p_transmit <= 1) { mc_trans = p_transmit; /* first event */ } else { mc_trans = p_trans; } mc_interact = 1 - mc_trans; if (mc_trans > 0 && (mc_trans >= 1 || rand01 () < mc_trans)) /* Transmit */ { p *= p_trans / mc_trans; intersect = 0; if (powder) { /* orientation of crystallite is no longer random */ randderotate (&kx, &ky, &kz, Alpha, Beta, Gamma); } type = 't'; if (!itype) itype = 1; #ifndef OPENACC hkl_info.type = type; #endif break; /* This break means that we are leaving the while-loop, exiting the crystal by "tunneling". */ } /* Scattering "proper", i.e. coh or incoh */ if (mc_interact <= 0) /* Protect against rounding errors */ { intersect = 0; if (powder) { /* orientation of crystallite is no longer random */ randderotate (&kx, &ky, &kz, Alpha, Beta, Gamma); } break; } /* First-pass considerations: */ if (!event_counter) p *= fabs (1 - p_trans) / mc_interact; /* Select a point at which to scatter the photon, taking secondary extinction into account. */ /* dP(l) = exp(-tot_xlen*l)dl P(l= tau_count) { #ifndef OPENACC if (hkl_info.flag_warning < 10) MPI_MASTER (fprintf (stderr, "Single_crystal: %s: WARNING: failed tau search " "(sum=%g, coh_refl=%g, j=%i, tau_count=%i). Using last reflection.\n", NAME_CURRENT_COMP, sum, coh_refl, j, tau_count);); #endif hkl_info.flag_warning++; j = tau_count - 1; } i = T[j].index; /* (8). Pick scattered wavevector kf from 2D Gauss distribution. */ z1 = randnorm (); z2 = randnorm (); y1 = T[j].l11 * z1 + T[j].y0x; y2 = T[j].l12 * z1 + T[j].l22 * z2 + T[j].y0y; kfx = T[j].rho_x + T[j].ox + T[j].b1x * y1 + T[j].b2x * y2; kfy = T[j].rho_y + T[j].oy + T[j].b1y * y1 + T[j].b2y * y2; kfz = T[j].rho_z + T[j].oz + T[j].b1z * y1 + T[j].b2z * y2; /* Normalize kf to length of ki, to account for planer approximation of the Ewald sphere. */ adjust = ki / sqrt (kfx * kfx + kfy * kfy + kfz * kfz); kfx *= adjust; kfy *= adjust; kfz *= adjust; /* Adjust photon weight (see manual for explanation). */ double pmul = T[j].xsect * coh_refl / (coh_xsect * T[j].refl); if (!isnan (pmul)) p *= pmul; kx = L[i].u1x * kfx + L[i].u2x * kfy + L[i].u3x * kfz; ky = L[i].u1y * kfx + L[i].u2y * kfy + L[i].u3y * kfz; kz = L[i].u1z * kfx + L[i].u2z * kfy + L[i].u3z * kfz; /* add thin layer k-broadening, only with box/cylinder shape */ if (hkl_info.shape == 1) { /* box */ double thickness_threshold = 1000 * 2 * PI / ki / 1e10; // 1000*lambda in Angs, /1e10 -> m if (xwidth < thickness_threshold) kx += 2 * PI / xwidth / 1e10 * randnorm (); if (yheight < thickness_threshold) ky += 2 * PI / yheight / 1e10 * randnorm (); if (zdepth < thickness_threshold) kz += 2 * PI / zdepth / 1e10 * randnorm (); } else if (hkl_info.shape == 0) { /* cylinder */ double thickness_threshold = 1000 * 2 * PI / ki / 1e10; if (yheight < thickness_threshold) ky += 2 * PI / yheight / 1e10 * randnorm (); } type = 'c'; if (!itype) itype = 3; #ifndef OPENACC hkl_info.type = type; hkl_info.h = L[i].h; hkl_info.k = L[i].k; hkl_info.l = L[i].l; #endif } /* ------------------------------------------------------------------------- */ /* lattice curvature option: rotate back photon velocity */ if (RX) { vec_rotate_2d (&kx, &kz, -curv_xangle); vec_rotate_2d (&Ex, &Ez, -curv_xangle); /*changing x,z actually curves the crystal, not only the planes*/ /*comment out if only curvature of the lattice planes is needed*/ vec_rotate_2d (&x, &z, -curv_xangle); } if (RY) { vec_rotate_2d (&ky, &kz, -curv_yangle); vec_rotate_2d (&Ey, &Ez, -curv_yangle); /*changing y,z actually curves the crystal, not only the planes*/ /*comment out if only curvature of the lattice planes is needed*/ vec_rotate_2d (&y, &z, -curv_yangle); } /* ------------------------------------------------------------------------- */ SCATTER; if (powder) { /* orientation of crystallite is no longer random */ randderotate (&kx, &ky, &kz, Alpha, Beta, Gamma); } /* Repeat loop for next scattering event. */ } while (intersect); /* end do (intersect) (multiple scattering loop) */ } /* if intersect */ %} FINALLY %{ MPI_MASTER ( if (hkl_info.flag_warning > 10) fprintf (stderr, "Single_crystal: %s: Warning messages were repeated %i times with absorbed photons/illegal tau search.\n", NAME_CURRENT_COMP, hkl_info.flag_warning); /* in case this instance is used in a SPLIT, we can recommend the optimal iteration value */ if (hkl_info.max_tau_count >= MCSX_REFL_SLIST_SIZE) { fprintf ( stderr, "Single_crystal: %s: WARNING: The reflection short list buffer was exhausted at least once. Please consider redefining MCSX_REFL_SLIST_SIZE > %d\n", NAME_CURRENT_COMP, MCSX_REFL_SLIST_SIZE); } if (hkl_info.nb_refl_count) { double split_iterations = (double)hkl_info.nb_reuses / hkl_info.nb_refl_count + 1; double split_optimal = (double)hkl_info.nb_refl / hkl_info.nb_refl_count; if (split_optimal > split_iterations + 5) { printf ("Single_crystal: %s: Info: you may improve the computation efficiency by using\n" " SPLIT %i COMPONENT %s=Single_crystal(order=1, ...)\n" " in the instrument description %s.\n", NAME_CURRENT_COMP, (int)split_optimal, NAME_CURRENT_COMP, instrument_source); } }); %} MCDISPLAY %{ if (hkl_info.shape == 0) { /* cylinder */ cylinder (0, 0, 0, radius, yheight, 0, 0, 1, 0); } else if (hkl_info.shape == 1) { /* box */ box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); } else if (hkl_info.shape == 2) { /* sphere */ sphere (0, 0, 0, radius); } else if (hkl_info.shape == 3) { /* OFF file */ off_display (offdata); } %} END