Calculation of reflections : update code to get list of permutations and...

Calculation of reflections : update code to get list of permutations and rotations for each crystal system; create function to read text file output from python

-  Five text files are created running gtCrystCalculateReflections:
      + [filename][type].txt with the list of reflections with the specific unit cell, lambda, sin(theta)/lambda limits
      + rot_[crystal_system].txt with the set of unitary rotation matrices corresponding to the indistinguasible lattice permutations
      + perm_[crystal_system].txt with the set of indistinguasible lattice permutations
      + sg_rot_[crystal_system].txt with rotations (complete list using symmetry operators)
      + sg_trans_[crystal_system].txt with transition values from symmetry operators

-  Renamed zUtilCryst/gtCrystGetReflections.m into zUtil_Help/gtReadTextFile.m

Signed-off-by: Laura Nervo <laura.nervo@esrf.fr>

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@953 4c865b51-4357-4376-afb4-474e03ccb993

zUtil_Cryst/gtCrystCalculateReflections.m

-function [hkl, theta, output_file] = gtCrystCalculateReflections(latticepar, mintheta, maxtheta, sg, typelist, filename)
+function [hkl, theta, hkil, sintheta] = gtCrystCalculateReflections(latticepar, mintheta, maxtheta, sg, typelist, filename)
 % GTCRYSTCALCULATEREFLECTIONS  Calculate the reflections list using python 
 %                              from fable library (genhkl_all or genhkl_unique
 %                              scripts)
-%     [hkl, theta, output_file] = gtCrystCalculateReflections(latticepar, mintheta, maxtheta, sg, typelist, filename)
-%     ---------------------------------------------------------------------------------------------------------------
+%     [hkl, theta, hkil, sintheta] = gtCrystCalculateReflections(latticepar, mintheta, maxtheta, sg, typelist, filename)
+%     ------------------------------------------------------------------------------------------------------------------
 %
 %     INPUT: 
 %       latticepar      = lattice parameters <double 1x6>
@@ -14,9 +14,10 @@ function [hkl, theta, output_file] = gtCrystCalculateReflections(latticepar, min
 %       filename = filename for output <string>
 %
 %     OUTPUT:
-%       hkl             = reflection list <double>
-%       theta           = theta value for each reflection <double>
-%       output_file     = path to the output filename <string>
+%       hkl      = reflection list <double>
+%       theta    = theta value for each reflection <double>
+%       hkil     = hkil indexes if needed <double>
+%       sintheta = sind(theta)/lambda <double>
 %
 %
 %     Version 001 03-12-2012 by LNervo
@@ -43,7 +44,7 @@ if ~exist('typelist','var') || isempty(typelist)
     typelist = 'all';
 end
 if ~exist('filename','var') || isempty(filename)
-    filename = 'reflections.txt';
+    filename = 'reflections';
 end
 
 latticepar_str = sprintf('%f ', latticepar);
@@ -55,15 +56,31 @@ cmd = sprintf('python26 %s %s %f %f %d %s %s',...
 [~,msg]=unix(cmd);
 disp(msg)
 
-[hkl,sintheta] = gtCrystGetReflections(filename);
+[C,Cmat] = gtReadTextFile([filename typelist '.txt'],'%f %f %f %f','\n','#',[1 3]);
+hkl = Cmat(:,1:3);
+sintheta = Cmat(:,4);
 
 if nargout > 1
     load parameters;
     energy = parameters.acq.energy;
     theta = asind(sintheta)*gtConvEnergyToWavelength(energy);
-    
     if nargout == 3
-        output_file = fullfile(pwd(),filename);
+      % for planes with Miller Indexes hkl -> hkil
+      hkil = hkl;
+      hkil(:,4) = hkl(:,3);
+      hkil(:,3) = -hkil(:,1)-hkil(:,2);
+
+      disp('  h   k   i   l   theta')
+      disp('-------------------------')
+      for ii=1:size(hkl,1)
+          fprintf('%3d %3d %3d %3d %7.3f\n',hkil(ii,:),theta(ii))
+      end
+    end
+else
+    disp('  h   k   l   sin(theta)/lambda')
+    disp('---------------------------------')
+    for ii=1:size(hkl,1)
+        fprintf('%3d %3d %3d %7.3f\n',hkl(ii,:),sintheta(ii))
     end
 end
 

zUtil_Cryst/gtCrystGetReflections.m

-function [hkl, sintheta] = gtCrystGetReflections(filename)
-% GTCRYSTGETREFLECTIONS  Reads output file from gtCrystCalculateReflections
-%
-%     [hkl, sintheta] = gtCrystGetReflections(filename)
-%     -------------------------------------------------
-%
-%     INPUT:
-%       filename = output file from gtCrystCalculateReflections <string> 
-%
-%     OUTPUT:
-%       hkl      = generated reflections <double>
-%       sintheta = sin(theta)/lambda <double>
-%
-%
-%     Version 001 03-12-2012 by LNervo
-
-
-fid = fopen(filename,'r');
-C = textscan(fid,'%f %f %f %f',...
-                 'Delimiter','\n',...
-                 'CommentStyle','#');
-fclose(fid);
-
-hkl = cell2mat(C(:,1:3));
-
-if nargout == 2
-    sintheta = cell2mat(C(:,4));
-end
-
-end % end of function

zUtil_Help/gtReadTextFile.m

+function [Ccell,Cmat] = gtReadTextFile(filename, formatS, delimiter, comments, minsize)
+% GTREADTEXTFILE  Reads textfiles into matrices
+%
+%     [Ccell, Cmat] = gtReadTextFile(filename, format, delimiter, comments, minsize)
+%     ------------------------------------------------------------------------------
+%
+%     INPUT:
+%       filename  = text file to be read <string>
+%       formatS   = format pattern as in fprintf <string> i.e. '%f %f %f'
+%       delimiter = character as delimiter i.e. ' '
+%       comments  = character as comments i.e. '#'
+%       minsize   = minimum size of element to be consider as object (like
+%                   matrix 3x3 or vector 1x3)
+%
+%     OUTPUT:
+%       Ccell = content of the text file (read as columns) <cell> 
+%       Cmat  = read matrix <double>
+%
+%
+%     Version 001 03-12-2012 by LNervo
+
+
+fid = fopen(filename,'r');
+C = textscan(fid, formatS,...
+             'Delimiter', delimiter,...
+             'CommentStyle', comments);
+fclose(fid);
+
+Cmat = cell2mat(C);
+
+if minsize(1)>0 && minsize(2)>0
+    for ii = 1 : size(Cmat,1)/minsize(1)
+        C{ii} = Cmat((ii-1)*minsize(1)+1:minsize(1)*ii,1:minsize(2));
+    end
+end
+
+Ccell = C;
+
+end % end of function

zUtil_Python/create_reflections_list.py

@@ -3,6 +3,12 @@
 '''
 Created on Dec 03, 2012
 @author: lnervo
+
+Based on the algorithm described in 
+    Le Page and Gabe (1979) J. Appl. Cryst., 12, 464-466
+
+and implemented in python by
+    Henning Osholm Sorensen, University of Copenhagen, July 22, 2010.
 '''
 
 import sys, os, numpy as np
@@ -11,10 +17,10 @@ from xfab import *
 from xfab import symmetry
 from xfab.symmetry import *
 from xfab.symmetry import tools
+from xfab import sg
 
-index = 1;
-_a = float(sys.argv[index])
-_b = float(sys.argv[index+1])
+_a = float(sys.argv[1])
+_b = float(sys.argv[2])
 _c = float(sys.argv[3])
 _alpha = float(sys.argv[4])
 _beta  = float(sys.argv[5])
@@ -32,9 +38,9 @@ print("    lattice parameters: (%f,%f,%f,%f,%f,%f)"
        % (_a, _b, _c, _alpha, _beta, _gamma))
 print("    sin(theta)/lambda range: %f,%f" % (_mintheta, _maxtheta))
 print("    spacegroup: %s" % _sgno)
-print("    type: %s" % _type)
+print("    type:       %s" % _type)
 print("")
-print("Reflections list calculation...")
+print("    *** Reflections list calculation ***")
 
 if _type == 'all':
     list = symmetry.tools.genhkl_all([_a,_b,_c,_alpha,_beta,_gamma],
@@ -48,14 +54,108 @@ else:
     print("Type not known...Existing")
     list = None
 
+_outputfile = _output + _type
+
 if list is not None:
-    print("...done!")
-    print("Saving reflections list in:\n    %s..." % _output)
+    print("Saving %d reflections generated for spacegroup %d - from symmetry.tools.genhkl_%s...\n"
+          % (len(list), _sgno, _type))
 
-    with file(_output,'w') as outfile:
+    with file(_outputfile + '.txt','w') as outfile:
         outfile.write('# Reflection list: {0}\n'.format(list.shape))
-        for data_hkl in list:
-            np.savetxt(outfile, data_hkl)
-            outfile.write('# New reflection\n')
+        outfile.write('# h k l sin(theta)/lambda\n')
+        np.savetxt(outfile, list, fmt=('%d','%d','%d','%f'))
 
     print("...done!")
+
+# list of symmetry operators for crystal system
+
+if _sgno != None:
+    spg = sg.sg(sgno=_sgno)
+
+    # number of symmetry operators
+    nsymop = spg.nsymop
+    # rotations list
+    rot    = spg.rot
+    # crystal system
+    crystal_system = spg.crystal_system
+    # trans info
+    trans = spg.trans
+
+    if rot is not None:
+        print("Saving %d rotations for crystal system '%s' - from sg.sg(sgno=%d).rot\n"
+              % (len(rot), crystal_system, _sgno))
+
+        with file("sg_rot_%s.txt" % crystal_system,'w') as outfile:
+            outfile.write('# Rotations list: {0}\n'.format(rot.shape))
+
+            for data_rot in rot:
+                np.savetxt(outfile, data_rot, fmt='%5.3f')
+
+                outfile.write('# New rotation\n')
+
+        print("...done!")
+
+    if trans is not None:
+        print("Saving %d trans for crystal system '%s' - from sg.sg(sgno=%d).trans\n"
+              % (len(trans), crystal_system, _sgno))
+
+        with file("sg_trans_%s.txt" % crystal_system,'w') as outfile:
+            outfile.write('# Trans list: {0}\n'.format(trans.shape))
+            np.savetxt(outfile, trans, fmt=('%f','%f','%f'))
+
+            outfile.write('# New trans\n')
+
+        print("...done!")
+
+
+    if crystal_system == "triclinic":
+        nsystem = 1
+    elif crystal_system == "monoclinic":
+        nsystem = 2
+    elif crystal_system == "orthorhombic":
+        nsystem = 3
+    elif crystal_system == "tetragonal":
+        nsystem = 4
+    elif crystal_system == "trigonal":
+        nsystem = 5
+    elif crystal_system == "hexagonal":
+        nsystem = 6
+    elif crystal_system == "cubic":
+        nsystem = 7
+    else:
+        print("Crystal system is not know")
+
+
+    # calculation of rotations
+    calc_rot = symmetry.rotations(nsystem)
+    # calc_rot_nf = calc_rot.reshape( calc_rot.shape[0]*calc_rot.shape[1], calc_rot.shape[2])
+
+    if calc_rot is not None:
+        print("Saving %d rotations for crystal system '%s' - from symmetry.rotations(%d)\n"
+              % (len(calc_rot), crystal_system, nsystem))
+
+        with file("rot_%s.txt" % crystal_system,'w') as outfile:
+            outfile.write('# Rotations list: {0}\n'.format(calc_rot.shape))
+            for data_rot in calc_rot:
+                np.savetxt(outfile, data_rot, fmt='%5.3f')
+                outfile.write('# New rotation\n')
+
+        print("...done!")
+
+    # calculation of permutations
+    calc_perm = symmetry.permutations(nsystem)
+
+    if calc_perm is not None:
+        print("Saving %d permutations for crystal system '%s' - from symmetry.permutations(%d)\n"
+              % (len(calc_perm), crystal_system, nsystem))
+
+        with file("perm_%s.txt" % crystal_system,'w') as outfile:
+            outfile.write('# Permutations list: {0}\n'.format(calc_perm.shape))
+            for data_perm in calc_perm:
+                np.savetxt(outfile, data_perm, fmt='%5.3f')
+                outfile.write('# New permutation\n')
+
+        print("...done!")
+
+else:
+    raise ValueError, 'No space group information given'