-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathhklGen.py
executable file
·1632 lines (1500 loc) · 70.7 KB
/
hklGen.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# hklGen.py
# Generates predicted diffraction peak positions for a given unit cell and space
# group, using the Fortran CFML library.
# Also uses the program "bumps" to perform fitting of calculated diffraction
# patterns to observed data.
# Created 6/4/2013
# Last edited 7/25/2013
import sys
import os
from math import floor, sqrt, log, tan, radians
from string import rstrip, ljust, rjust, center
from copy import deepcopy
from ordereddict import OrderedDict
from ctypes import cdll, Structure, c_int, c_float, c_char, c_bool, c_char_p, \
c_void_p, c_ulong, addressof, sizeof, POINTER
import numpy as np
import pylab
try:
from bumps.names import Parameter, FitProblem
except(ImportError):
pass
# This should be the location of the CFML library
LIBFILE = "libcrysfml.so"
LIBPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), LIBFILE)
lib = cdll[LIBPATH]
# DVDim: contains array dimension information. Attributes:
# stride_mult - stride multiplier for the dimension
# lower_bound - first index for the dimension
# upper_bound - last index for the dimension
class DVDim(Structure):
_fields_ = [("stride_mult", c_ulong), ("lower_bound", c_ulong),
("upper_bound", c_ulong)]
# DV: dope vector for gfortran that passes array information. Attributes:
# base_addr - base address for the array
# base - base offset
# dtype - contains the element size, type (3 bits), and rank (3 bits)
# dim - DVDim object for the vector
class DV(Structure):
_fields_ = [("base_addr",c_void_p), ("base", c_void_p), ("dtype", c_ulong),
("dim", DVDim*7)]
def __len__(self):
return self.dim[0].upper_bound - self.dim[0].lower_bound + 1
def elemsize(self):
return self.dv.dtype >> 6
def rank(self):
return self.dv.dtype&7
def data(self, dataType):
return (dataType*len(self)).from_address(self.base_addr)
# dv_dtype: calculates the "dtype" attribute for a dope vector
def dv_dtype(size,type,rank): return size*64+type*8+rank
# build_struct_dv: constructs a dope vector for an array of derived types
def build_struct_dv(array):
dv = DV()
dv.base_addr = addressof(array)
dv.base = c_void_p()
dv.dtype = dv_dtype(sizeof(array[0]), 5, 1) # 5 = derived type
dv.dim[0].stride_mult = 1
dv.dim[0].lower_bound = 1
dv.dim[0].upper_bound = len(array)
return dv
# deconstruct_dv: converts a dope vector into an array of any type (currently
# only implemented for 1-D arrays)
def deconstruct_dv(dv, dataType):
# dims = dv.dtype & 7
elem_size = dv.dtype >> 6
size = dv.dim[0].upper_bound - dv.dim[0].lower_bound + 1
array = [None for i in xrange(size)]
for i in xrange(size):
offset = i * elem_size * dv.dim[0].stride_mult
array[i] = dataType.from_address(dv.base_addr + offset)
return array
# SymmetryOp Attributes:
# rot - rotational part of symmetry operator (3 by 3 matrix)
# trans - translational part of symmetry operator (vector)
class SymmetryOp(Structure):
_fields_ = [("rot", c_int*3*3), ("trans", c_float*3)]
# MagSymmetryOp attributes:
# rot - roatational part of symmetry operator
# phase - phase as a fraction of 2*pi
class MagSymmetryOp(Structure):
_fields_ = [("rot", c_int*3*3), ("phase", c_float)]
# WyckoffPos attributes:
# multip - multiplicity
# site - site symmetry
# numElements - number of elements in orbit
# origin - origin
# orbit - strings containing orbit information
class WyckoffPos(Structure):
_fields_ = [("multip", c_int), ("site", c_char*6),
("numElements", c_int), ("origin", c_char*40),
("orbit", c_char*40*48)]
# Wyckoff attributes:
# numOrbits - number of orbits
# orbits - list of Wyckoff position objects
class Wyckoff(Structure):
_fields_ = [("numOrbits", c_int), ("orbits", WyckoffPos*26)]
# SpaceGroup attributes:
# number - space group number
# symbol - Hermann-Mauguin symbol
# hallSymbol - Hall symbol
# xtalSystem - Crystal system
# laue - Laue class
# pointGroup - corresponding point group
# info - additional information
# setting - space group setting information (IT, KO, ML, ZA, Table,
# Standard, or UnConventional)
# hex - true if space group is hexagonal
# lattice - lattice type
# latticeSymbol - lattice type symbol
# latticeNum - number of lattice points in a cell
# latticeTrans - lattice translations
# bravais - Bravais symbol and translations
# centerInfo - information about symmetry center
# centerType - 0 = centric (-1 not at origin), 1 = acentric,
# 2 = centric (-1 at origin)
# centerCoords - fractional coordinates of inversion center
# numOps - number of symmetry operators in the reduced set
# multip - multiplicity of the general position
# numGens - minimum number of operators to generate the group
# symmetryOps - list of symmetry operators
# symmetryOpsSymb - string form of symmetry operator objects
# wyckoff - object containing Wyckoff information
# asymmetricUnit - direct space parameters for the asymmetric unit
class SpaceGroup(Structure):
_fields_ = [("number", c_int), ("symbol", c_char*20),
("hallSymbol", c_char*16), ("xtalSystem", c_char*12),
("laue", c_char*5), ("pointGroup", c_char*5),
("info", c_char*5), ("setting", c_char*80),
("hex", c_int), ("lattice", c_char),
("latticeSymbol", c_char*2), ("latticeNum", c_int),
("latticeTrans", c_float*3*12), ("bravais", c_char*51),
("centerInfo", c_char*80), ("centerType", c_int),
("centerCoords", c_float*3), ("numOps", c_int),
("multip", c_int), ("numGens", c_int),
("symmetryOps", SymmetryOp*192),
("symmetryOpsSymb", c_char*40*192),
("wyckoff", Wyckoff), ("asymmetricUnit", c_float*6)]
def __init__(self, groupName=None):
if (groupName != None):
fn = getattr(lib, "__cfml_crystallographic_symmetry_MOD_set_spacegroup")
fn.argtypes = [c_char_p, POINTER(SpaceGroup), c_void_p, POINTER(c_int),
c_char_p, c_char_p, c_int, c_int, c_int]
fn.restype = None
fn(groupName, self, None, None, None, None, len(groupName), 0, 0)
# CrystalCell attributes:
# length, angle - arrays of unit cell parameters
# lengthSD, angleSD - standard deviations of parameters
# rLength, rAngle - arrays of reciprocal cell parameters
# GD, GR - direct and reciprocal space metric tensors
# xtalToOrth, orthToXtal - matrices to convert between orthonormal and
# crystallographic bases
# BLB, invBLB - Busing-Levy B-matrix and its inverse
# volume, rVolume - direct and reciprocal cell volumes
# cartType - Cartesian reference frame type (cartType = 'A'
# designates x || a)
class CrystalCell(Structure):
_fields_ = [("length", c_float*3), ("angle", c_float*3),
("lengthSD", c_float*3), ("angleSD", c_float*3),
("rLength", c_float*3), ("rAngle", c_float*3),
("GD", c_float*9), ("GR", c_float*9),
("xtalToOrth", c_float*9), ("orthToXtal", c_float*9),
("BLB", c_float*9), ("invBLB", c_float*9),
("volume", c_float), ("rVolume", c_float),
("cartType", c_char)]
def __init__(self, length=None, angle=None):
if (length != None):
self.setCell(length, angle)
def setCell(self, length, angle):
fn = getattr(lib, "__cfml_crystal_metrics_MOD_set_crystal_cell")
float3 = c_float*3
fn.argtypes = [POINTER(float3), POINTER(float3), POINTER(CrystalCell),
POINTER(c_char), POINTER(float3)]
fn.restype = None
fn(float3(*length), float3(*angle), self, None, None)
# MagSymmetry attributes:
# [corresponds to CFML MagSymm_k_type]
# name - name describing magnetic symmetry
# SkType - "Spherical_Frame" designates input Fourier coefficients
# should be in spherical coordinates
# lattice - lattice type
# irrRepsNum - number of irreducible representations (max 4)
# magSymOpsNum - number of magnetic symmetry operators per crystallographic
# symmetry operator (max 8)
# centerType - 0 = centric (-1 not at origin), 1 = acentric,
# 2 = centric (-1 at origin)
# magCenterType - 1 = acentric magnetic symmetry, 2 = centric
# numk - number of independent propagation vectors
# k - propagation vectors
# numCentVec - number of centering lattice vectors
# centVec - centering lattice vectors
# numSymOps - number of crystallographic symmetry operators
# multip - multiplicity of the space group
# numBasisFunc - number of basis functions per representation (can be
# given a negative value to indicate a complex basis)
# coeffType - 0 = real basis function coefficients, 1 = pure imaginary
# symOpsSymb - symbolic form of symmetry operators
# symOps - crystallographic symmetry operators
# magSymOpsSymb - symbolic form of magnetic operators
# magSymOps - magnetic symmetry operators
# basis - coeffs of basis functions of irreducible representations
class MagSymmetry(Structure):
_fields_ = [("name", c_char*31), ("SkType", c_char*15), ("lattice", c_char),
("numIrreps", c_int), ("magSymOpsNum", c_int),
("centerType", c_int), ("magCenterType", c_int),
("numk", c_int), ("k", c_float*3*12), ("numCentVec", c_int),
("centVec", c_float*3*4), ("numSymOps", c_int),
("multip", c_int), ("numBasisFunc", c_int*4),
("coeffType", c_int*12*4), ("symOpsSymb", c_char*40*48),
("symOps", SymmetryOp*48), ("magSymOpsSymb", c_char*40*48*8),
("magSymOps", MagSymmetryOp*48*8), ("basis", c_float*2*3*12*48*4)]
def setBasis(self, irrRepNum, symOpNum, vectorNum, v):
c_array2 = c_float*2
self.basis[irrRepNum][symOpNum][vectorNum] = \
(c_array2*3)(c_array2(v[0].real, v[0].imag),
c_array2(v[1].real, v[1].imag),
c_array2(v[2].real, v[2].imag))
# Atom attributes:
# label - label for the atom
# element - chemical symbol of the element
# strFactSymb - chemical symbol used in the structure factor
# active - used for program control
# atomicNum - atomic number
# multip - site multiplicity
# coords - fractional coordinates
# occupancy - site occupancy
# BIso - isotropic temperature (Debye-Waller) factor
# uType - type of anisotropic thermal factor: "u_ij", "b_ij", "beta",
# or none
# thType - thermal factor type: "isotr", "aniso", or "other"
# U - matrix entries U11, U22, U33, U12, U13, and U23
# UEquiv - equivalent U
# charge - self-explanatory
# moment - magnetic moment
# index - index "for different purposes", whatever mysterious purposes
# those might be
# numVars - a variable completely lacking in documentation
# freeVars - free variables used "for different purposes" again
# atomInfo - string containing miscellaneous information
# *** mult, lsq prefixes indicate multipliers and positions in the list of
# LSQ parameters, respectively; SD indicates standard deviation
class Atom(Structure):
_fields_ = [("label", c_char*10), ("element", c_char*2),
("strFactSymb", c_char*4), ("active", c_int),
("atomicNum", c_int), ("multip", c_int), ("coords", c_float*3),
("coordsSD", c_float*3), ("multCoords", c_float*3),
("lsqCoords", c_int*3), ("occupancy", c_float),
("occupancySD", c_float), ("multOccupancy", c_float),
("lsqOccupancy", c_int), ("BIso", c_float), ("BIsoSD", c_float),
("multBIso", c_float), ("lsqBIso", c_float),
("uType", c_char*4), ("thType", c_char*5), ("U", c_float*6),
("USD", c_float*6), ("UEquiv", c_float), ("multU", c_float*6),
("lsqU", c_int*6), ("charge", c_float), ("moment", c_float),
("index", c_int*5), ("numVars", c_int),
("freeVars", c_float*10), ("atomInfo", c_char*40)]
def __init__(self, *args):
# construct an atom from a list of attributes
if (len(args) == 6):
init = getattr(lib, "__cfml_atom_typedef_MOD_init_atom_type")
init.argtypes = [POINTER(Atom)]
init.restype = None
init(self)
self.label = ljust(args[0], 10)
self.element = ljust(args[1], 2)
self.strFactSymb = ljust(self.element, 4)
self.coords = (c_float*3)(*args[2])
self.multip = args[3]
self.occupancy = c_float(args[4])
self.BIso = c_float(args[5])
# # copy over attributes from a magnetic atom
# for field in Atom._fields_:
# setattr(self, field[0], getattr(magAtom, field[0]))
def sameSite(self, other):
# TODO: make this work for equivalent sites, not just identical ones
# returns true if two atoms occupy the same position
# Warning: they must be specified with identical starting coordinates
eps = 0.001
return all([approxEq(self.coords[i], other.coords[i], eps)
for i in xrange(3)])
# MagAtom attributes: same as atom, plus:
# numkVectors - number of propagation vectors (excluding -k)
# irrepNum - index of the irreducible representation to be used
# SkReal - real part of Fourier coefficient
# SkRealSphere - real part of Fourier coefficient (spherical coordinates)
# SkIm - imaginary part of Fourier coefficient
# SkimSphere - imaginary part of Fourier coefficient (spherical coords)
# phase - magnetic phase (fraction of 2*pi)
# basis - coefficients of the basis functions
class MagAtom(Structure):
_fields_ = [("label", c_char*10), ("element", c_char*2),
("strFactSymb", c_char*4), ("active", c_int),
("atomicNum", c_int), ("multip", c_int), ("coords", c_float*3),
("coordsSD", c_float*3), ("multCoords", c_float*3),
("lsqCoords", c_int*3), ("occupancy", c_float),
("occupancySD", c_float), ("multOccupancy", c_float),
("lsqOccupancy", c_int), ("BIso", c_float), ("BIsoSD", c_float),
("multBIso", c_float), ("lsqBIso", c_float),
("uType", c_char*4), ("thType", c_char*5), ("U", c_float*6),
("USD", c_float*6), ("UEquiv", c_float), ("multU", c_float*6),
("lsqU", c_int*6), ("charge", c_float), ("moment", c_float),
("index", c_int*5), ("numVars", c_int),
("freeVars", c_float*10), ("atomInfo", c_char*40),
("numkVectors", c_int), ("irrepNum", c_int*12),
("SkReal", c_float*3*12), ("SkRealSphere", c_float*3*12),
("multSkReal", c_float*3*12), ("lsqSkReal", c_int*3*12),
("SkIm", c_float*3*12), ("SkImSphere", c_float*3*12),
("multSkIm", c_float*3*12), ("lsqSkIm", c_int*3*12),
("phase", c_float*12), ("multPhase", c_float*12),
("lsqPhase", c_int*12), ("basis", c_float*12*12),
("multBasis", c_float*12*12), ("lsqBasis", c_int*12*12)]
def sameSite(self, other):
# returns true if two atoms occupy the same position
# Warning: they must be specified with identical starting coordinates
eps = 0.001
return all([approxEq(self.coords[i], other.coords[i], eps)
for i in xrange(3)])
# AtomList attributes:
# numAtoms - the number of atoms
# atoms - a list of Atom objects
# magnetic - True if this is a list of MagAtoms
class AtomList(Structure):
_fields_ = [("numAtoms", c_int), ("atoms", DV)]
def __init__(self, atoms=None, magnetic=False):
self.magnetic = magnetic
if (atoms != None):
init = getattr(lib, "__cfml_atom_typedef_MOD_allocate_atom_list")
init.argtypes = [POINTER(c_int), POINTER(AtomList), POINTER(c_int)]
init.restype = None
numAtoms = len(atoms)
init(c_int(numAtoms), self, None)
self.numAtoms = c_int(numAtoms)
# copy information from provided atom list
for i, atom in enumerate(self):
for field in atom._fields_:
setattr(atom, field[0], getattr(atoms[i], field[0]))
# atoms = self.atoms.data(Atom)
# c_array3 = c_float*3
# if not isinstance(elements[0], Atom):
# # build atoms from scratch
# for i in xrange(numAtoms):
# atoms[i].element = elements[i]
# atoms[i].atomicNum = c_int(atomicNums[i])
# atoms[i].coords = c_array3(*coords[i])
# atoms[i].occupancy = c_float(occupancies[i])
# atoms[i].BIso = c_float(BIsos[i])
# atoms[i].label = elements[i]
def __len__(self):
return int(self.numAtoms)
def __getitem__(self, index):
if (index < 0): index += len(self)
if (self.magnetic): dtype = MagAtom
else: dtype = Atom
return self.atoms.data(dtype)[index]
# Reflection attributes:
# hkl - list containing hkl indices for the reflection
# multip - multiplicity
# FObs - observed structure factor
# FCalc - calculated structure factor
# FSD - standard deviation of structure factor
# s - s = sin(theta)/lambda = 1/(2d) [No 4*pi factor!]
# weight - weight of reflection
# phase - phase angle in degrees
# realPart - real part of the structure factor
# imPart - imaginary part of the structure factor
# aa - currently unused
# bb - currently unused
class Reflection(Structure):
_fields_ = [("hkl", c_int*3), ("multip", c_int), ("FObs", c_float),
("FCalc", c_float), ("FSD", c_float), ("s", c_float),
("weight", c_float), ("phase", c_float), ("realPart", c_float),
("imPart", c_float), ("aa", c_float), ("bb", c_float)]
# MagReflection attributes:
# [corresponds to CFML MagH_type]
# equalMinus - True if k is equivalent to -k
# multip - multiplicity
# knum - index of the propagation vector (k)
# signk - equal to +1 for -k and -1 for +k, because somebody
# thought that labeling system was logical
# s - sin(theta)/lambda
# magIntVecSq - norm squared of the magnetic interaction vector
# hkl - reciprocal scattering vector +/- k
# magStrFact - magnetic structure factor
# magIntVec - magnetic interaction vector
# magIntVecCart - magnetic interaction vector (Cartesian coordinates)
class MagReflection(Structure):
_fields_ = [("equalMinus", c_int), ("multip", c_int), ("knum", c_int),
("signk", c_float), ("s", c_float), ("magIntVecSq", c_float),
("hkl", c_float*3), ("magStrFact", c_float*2*3),
("magIntVec", c_float*2*3), ("magIntVecCart", c_float*2*3)]
# can initialize this from a regular (non-magnetic) reflection
def __init__(self, reflection=None):
if (reflection != None):
self.hkl = reflection.hkl
self.multip = reflection.multip
self.s = reflection.s
# ReflectionList attributes
# numReflections - the number of reflections
# reflections - list of Reflection objects
# magnetic - True if this is a list of MagReflections
class ReflectionList(Structure):
_fields_ = [("numReflections", c_int), ("reflections", DV)]
def __init__(self, magnetic=False):
self.magnetic = magnetic
def __len__(self):
return int(self.numReflections)
def __getitem__(self, index):
if (index < 0): index += len(self)
if (self.magnetic): dtype = MagReflection
else: dtype = Reflection
return self.reflections.data(dtype)[index]
# FileList: represents a Fortran file object
class FileList(Structure):
_fields_ = [("numLines", c_int), ("lines", DV)]
def __init__(self, filename):
makeList = getattr(lib, "__cfml_io_formats_MOD_file_to_filelist")
makeList.argtypes = [c_char_p, POINTER(FileList), c_int]
makeList.restype = None
makeList(filename, self, len(filename))
# Gaussian: represents a Gaussian function that can be evaluated at any
# 2*theta value. u, v, and w are fitting parameters.
class Gaussian(object):
# TODO: add option for psuedo-Voigt peak shape
scaleFactor = 1
def __init__(self, center, u, v, w, I, hkl=[0,0,0]):
self.C0 = 4*log(2)
self.center = center # 2*theta position
self.u = u
self.v = v
self.w = w
self.I = I
try:
self.H = sqrt(u*(tan(radians(center/2))**2)
+ v*tan(radians(center/2)) + w)
self.scale = self.I * sqrt(self.C0/np.pi)/self.H * Gaussian.scaleFactor
except ValueError:
self.H = 0
self.scale = 0
self.hkl = hkl
# __call__: returns the value of the Gaussian at some 2*theta positions
def __call__(self, x):
return self.scale * np.exp(-self.C0*(x-self.center)**2/self.H**2)
def add(self, v, x):
# only add to nearby 2*theta positions
idx = (x>self.center-self.H*3) & (x<self.center+self.H*3)
v[idx] += self.__call__(x[idx])
# LinSpline: represents a linear spline function to be used for the background
class LinSpline(object):
def __init__(self, arg1=None, arg2=None):
if (arg1 == None):
# create default uniform 0 background
self.x = [0,1]
self.y = [0,0]
elif isSequence(arg1):
# read in x and y coordinates from lists
self.x = np.copy(arg1)
self.y = np.copy(arg2)
elif (type(arg1) == str):
# read in x and y coordinates from a file
self.x, self.y = np.loadtxt(arg1, dtype=float, skiprows=5, unpack=True)
else:
# create a uniform background with the numeric value of arg1
self.x = [0,1]
self.y = [arg1, arg1]
# __call__: returns the interpolated y value at some x position
def __call__(self, x):
# locate the two points to interpolate between
return np.interp(x, self.x, self.y)
def __repr__(self):
return "LinSpline(" + str(self.x) + ", " + str(self.y) + ")"
# approxEq: returns True if two floats are equal to within some tolerance
def approxEq(num1, num2, eps):
return abs(num1-num2) <= eps
# isSequence: returns True if a value is a list, tuple, numpy array, etc. and
# False if it is a string or a non-sequence
def isSequence(x):
return ((not hasattr(x, "strip") and hasattr(x, "__getslice__")) or
hasattr(x, "__iter__"))
# twoTheta: converts a sin(theta)/lambda position to a 2*theta position
def twoTheta(s, wavelength):
if (s*wavelength >= 1): return 180.0
if (s*wavelength <= 0): return 0.0
return 2*np.degrees(np.arcsin(s*wavelength))
# getS: converts a 2*theta position to a sin(theta)/lambda position
def getS(tt, wavelength):
return np.sin(np.radians(tt/2))/wavelength
# hklString: converts an hkl list into a string
def hklString(hkl):
try:
for x in hkl:
assert(x == int(x))
return "%d %d %d" % tuple(hkl)
except(AssertionError):
return "%.1f %.1f %.1f" % tuple(hkl)
# calcS: calculates the sin(theta)/lambda value for a given list of planes
def calcS(cell, hkl):
if isSequence(hkl[0]):
# list of hkl positions
return [calcS(cell, h) for h in hkl]
else:
# single hkl
fn = lib.__cfml_reflections_utilities_MOD_hs_r
float3 = c_float*3
fn.argtypes = [POINTER(float3), POINTER(CrystalCell)]
fn.restype = c_float
return float(fn(float3(*hkl), cell))
# applySymOp: applies a symmetry operator to a given vector and normalizes the
# resulting vector to stay within one unit cell
def applySymOp(v, symOp):
rotMat = np.mat(symOp.rot)
transMat = np.mat(symOp.trans)
vMat = np.mat(v).T
newV = rotMat * vMat + transMat
return np.array(newV.T[0])%1
# dotProduct: returns the dot product of two complex vectors (stored as 3x2
# Numpy arrays) or a list of two complex vectors
def dotProduct(v1, v2):
if (not isSequence(v1[0][0]) or len(v1[0][0]) == 1):
u1 = np.matrix([complex(v1[0][0], -v1[0][1]),
complex(v1[1][0], -v1[1][1]),
complex(v1[2][0], -v1[2][1])])
u2 = np.matrix([[complex(v2[0][0], v2[0][1])],
[complex(v2[1][0], v2[1][1])],
[complex(v2[2][0], v2[2][1])]])
dot = np.dot(u1, u2)
return np.array(dot)[0][0]
else:
return np.array([dotProduct(u1,u2) for u1, u2 in zip(v1, v2)])
# inputInfo: requests space group and cell information, wavelength, and range
# of interest
def inputInfo():
# Input the space group name and create it
groupName = raw_input("Enter space group "
"(HM symbol, Hall symbol, or number): ")
spaceGroup = SpaceGroup(groupName)
# Remove excess spaces from Fortran
spaceGroup.xtalSystem = rstrip(spaceGroup.xtalSystem)
length = [0, 0, 0]
angle = [0, 0, 0]
# Ask for parameters according to the crystal system and create the cell
if (spaceGroup.xtalSystem.lower() == "triclinic"):
length[0], length[1], length[2], angle[0], angle[1], angle[2] = \
input("Enter cell parameters (a, b, c, alpha, beta, gamma): ")
elif (spaceGroup.xtalSystem.lower() == "monoclinic"):
angle[0] = angle[2] = 90
length[0], length[1], length[2], angle[1] = \
input("Enter cell parameters (a, b, c, beta): ")
elif (spaceGroup.xtalSystem.lower() == "orthorhombic"):
angle[0] = angle[1] = angle[2] = 90
length[0], length[1], length[2] = \
input("Enter cell parameters (a, b, c): ")
elif (spaceGroup.xtalSystem.lower() == "tetragonal"):
angle[0] = angle[1] = angle[2] = 90
length[0], length[2] = input("Enter cell parameters (a, c): ")
length[1] = length[0]
elif (spaceGroup.xtalSystem.lower() in ["rhombohedral", "hexagonal"]):
angle[0] = angle[1] = 90
angle[2] = 120
length[0], length[2] = input("Enter cell parameters (a, c): ")
length[1] = length[0]
elif (spaceGroup.xtalSystem.lower() == "cubic"):
angle[0] = angle[1] = angle[2] = 90
length[0] = input("Enter cell parameter (a): ")
length[1] = length[2] = length[0]
cell = CrystalCell(length, angle)
# Input the wavelength and range for hkl calculation (and adjust the range
# if necessary)
wavelength = input("Enter the wavelength: ")
sMin, sMax = input("Enter the sin(theta)/lambda interval: ")
adjusted = False
if (sMin < 0):
sMin = 0
adjusted = True
if (sMax > 1.0/wavelength):
sMax = 1.0/wavelength
adjusted = True
if (adjusted):
print "sin(theta)/lambda interval adjusted to [%f, %f]" % (sMin, sMax)
return (spaceGroup, cell, wavelength, sMin, sMax)
# readInfo: acquires cell, space group, and atomic information from a .cif,
# .cfl, .pcr, or .shx file
def readInfo(filename):
fn = lib.__cfml_io_formats_MOD_readn_set_xtal_structure_split
fn.argtypes = [c_char_p, POINTER(CrystalCell), POINTER(SpaceGroup),
POINTER(AtomList), c_char_p, POINTER(c_int),
c_void_p, c_void_p, c_char_p, c_int, c_int, c_int]
fn.restype = None
cell = CrystalCell()
spaceGroup = SpaceGroup()
atomList = AtomList()
ext = filename[-3:]
fn(filename, cell, spaceGroup, atomList, ext, None, None, None,
None, len(filename), len(ext), 0)
return (spaceGroup, cell, atomList)
# readMagInfo: acquires cell, space group, atomic, and magnetic information
# from a .cfl file
def readMagInfo(filename):
fn = lib.__cfml_magnetic_symmetry_MOD_readn_set_magnetic_structure
fn.argtypes = [POINTER(FileList), POINTER(c_int), POINTER(c_int),
POINTER(MagSymmetry), POINTER(AtomList), c_void_p,
c_void_p, POINTER(CrystalCell)]
fn.restype = None
info = readInfo(filename)
spaceGroup = info[0]
cell = info[1]
symmetry = MagSymmetry()
atomList = AtomList(magnetic=True)
fileList = FileList(filename)
print "reading information from " + filename
fn(fileList, c_int(0), c_int(0), symmetry, atomList, None, None, cell)
return (spaceGroup, cell, atomList, symmetry)
# readData: reads in a data file for the observed intensities and 2*theta
# values
# xy: 2*theta and intensity values in two-column format
# y: data file only contains intensities for a linear 2*theta range
def readData(filename, kind="xy", skiplines=0, skipcols=0, colstep=1,
start=None, stop=None, step=None, exclusions=None):
if (kind == "xy"):
tt, observed = np.loadtxt(filename, dtype=float, usecols=(0,1),
skiprows=skiplines, unpack=True)
elif (kind == "y"):
tt = np.linspace(start, stop, round((stop-start)/float(step) + 1))
data = np.loadtxt(filename, dtype=float, skiprows=skiplines)
observed = data.flatten()
observed = observed[skipcols:skipcols+colstep*len(tt):colstep]
tt, observed = removeRange(tt, exclusions, observed)
return (tt, observed)
# getMaxNumRef: returns the maximum number of reflections for a given cell
def getMaxNumRef(sMax, volume, sMin=0.0, multip=2):
fn = lib.__cfml_reflections_utilities_MOD_get_maxnumref
fn.argtypes = [POINTER(c_float), POINTER(c_float), POINTER(c_float),
POINTER(c_int)]
fn.restype = c_int
numref = fn(c_float(sMax), c_float(volume), c_float(sMin), c_int(multip))
return numref
# hklGen: generates a list of reflections in a specified range
# If getList is true, returns a ReflectionList object
def hklGen(spaceGroup, cell, sMin, sMax, getList=False, xtal=False):
# Calculate the reflection positions
maxReflections = getMaxNumRef(sMax+0.2, cell.volume, multip=spaceGroup.multip)
# Create a reference that will be modified by calling Fortran
reflectionCount = c_int(maxReflections)
if (not getList):
c_ReflectionArray = Reflection*max(maxReflections,1)
reflections = c_ReflectionArray()
if xtal:
# single crystal reflections (also used for magnetic structures)
# This option is currently non-functioning (use genList=True instead)
fn = lib.__cfml_reflections_utilities_MOD_hkl_gen_sxtal_reflection
fn.argtypes = [POINTER(CrystalCell), POINTER(SpaceGroup),
POINTER(c_float), POINTER(c_float), POINTER(c_int),
POINTER(DV), POINTER(c_int*3), POINTER(c_int*3*2)]
fn.restype = None
fn(cell, spaceGroup, c_float(sMin), c_float(sMax), reflectionCount,
build_struct_dv(reflections), None, None)
else:
# powder reflections
fn = lib.__cfml_reflections_utilities_MOD_hkl_uni_reflection
fn.argtypes = [POINTER(CrystalCell), POINTER(SpaceGroup), POINTER(c_bool),
POINTER(c_float), POINTER(c_float), c_char_p,
POINTER(c_int), POINTER(DV), POINTER(c_bool)]
fn.restype = None
fn(cell, spaceGroup, c_bool(True), c_float(sMin), c_float(sMax), 's',
reflectionCount, build_struct_dv(reflections), c_bool(False))
else:
if xtal:
fn = lib.__cfml_reflections_utilities_MOD_hkl_gen_sxtal_list
fn.argtypes = [POINTER(CrystalCell), POINTER(SpaceGroup),
POINTER(c_float), POINTER(c_float), POINTER(c_int),
POINTER(ReflectionList), POINTER(c_int*3),
POINTER(c_int*3*2)]
fn.restype = None
reflections = ReflectionList()
fn(cell, spaceGroup, c_float(sMin), c_float(sMax), reflectionCount,
reflections, None, None)
else:
fn = lib.__cfml_reflections_utilities_MOD_hkl_uni_refllist
fn.argtypes = [POINTER(CrystalCell), POINTER(SpaceGroup),
POINTER(c_bool), POINTER(c_float), POINTER(c_float),
c_char_p, POINTER(c_int), POINTER(ReflectionList),
POINTER(c_bool)]
fn.restype = None
reflections = ReflectionList()
fn(cell, spaceGroup, c_bool(True), c_float(sMin), c_float(sMax), 's',
reflectionCount, reflections, c_bool(False))
if (not isinstance(reflections, ReflectionList)):
reflections = reflections[:reflectionCount.value]
return reflections
# satelliteGen: generates a list of magnetic satellite reflections below a
# maximum sin(theta)/lambda value
def satelliteGen(cell, symmetry, sMax):
fn = lib.__cfml_magnetic_structure_factors_MOD_gen_satellites
fn.argtypes = [POINTER(CrystalCell), POINTER(MagSymmetry), POINTER(c_float),
POINTER(ReflectionList), POINTER(c_bool), POINTER(c_bool),
POINTER(ReflectionList)]
fn.restype = None
refList = ReflectionList(True)
fn(cell, symmetry, c_float(sMax), refList, c_bool(True), c_bool(True), None)
return refList
# calcStructFact: calculates the structure factor squared for a list of planes
# using provided atomic positions
def calcStructFact(refList, atomList, spaceGroup, wavelength):
init = lib.__cfml_structure_factors_MOD_init_calc_strfactors
init.argtypes = [POINTER(ReflectionList), POINTER(AtomList),
POINTER(SpaceGroup), c_char_p, POINTER(c_float),
POINTER(c_int), c_int]
init.restype = None
init(refList, atomList, spaceGroup, 'NUC', c_float(wavelength), None, 3)
calc = lib.__cfml_structure_factors_MOD_calc_strfactor
calc.argtypes = [c_char_p, c_char_p, POINTER(c_int), POINTER(c_float),
POINTER(AtomList), POINTER(SpaceGroup), POINTER(c_float),
POINTER(DV), POINTER(c_float*2), c_int, c_int]
calc.restype = None
structFacts = [c_float() for i in xrange(refList.numReflections)]
reflections = refList[:]
for i, reflection in enumerate(reflections):
# calculates the square of the structure factor
calc('P', 'NUC', c_int(i+1), c_float(reflection.s**2), atomList,
spaceGroup, structFacts[i], None, None, 1, 3)
# convert from c_float to float
structFacts = [sf.value for sf in structFacts]
return structFacts
# calcMagStructFact: calculates the magnetic structure factors around a list
# of lattice reflections
def calcMagStructFact(refList, atomList, symmetry, cell):
init = lib.__cfml_magnetic_structure_factors_MOD_init_mag_structure_factors
init.argtypes = [POINTER(ReflectionList), POINTER(AtomList),
POINTER(MagSymmetry), POINTER(c_int)]
init.restype = None
init(refList, atomList, symmetry, None)
calc = lib.__cfml_magnetic_structure_factors_MOD_mag_structure_factors
calc.argtypes = [POINTER(AtomList), POINTER(MagSymmetry),
POINTER(ReflectionList)]
calc.restype = None
calc(atomList, symmetry, refList)
# calculate the "magnetic interaction vector" (the square of which is
# proportional to the intensity)
calcMiv = lib.__cfml_magnetic_structure_factors_MOD_calc_mag_interaction_vector
calcMiv.argtypes = [POINTER(ReflectionList), POINTER(CrystalCell)]
calcMiv.restype = None
calcMiv(refList, cell)
# strFacts = np.array([ref.magStrFact for ref in refList])
mivs = np.array([ref.magIntVec for ref in refList])
return mivs
# Ye Olde Function:
# fn = lib.__cfml_magnetic_structure_factors_MOD_calc_magnetic_strf_miv
# fn.argtypes = [POINTER(CrystalCell), POINTER(MagSymmetry),
# POINTER(AtomList), POINTER(MagReflection)]
# fn.restype = None
#
# float3 = c_float*3
# equiv = lib.__cfml_propagation_vectors_MOD_k_equiv_minus_k
# equiv.argtypes = [POINTER(float3), c_char_p, c_int]
# equiv.restype = c_bool
# equalMinus = [equiv(float3(*symmetry.k[i]), symmetry.lattice,
# len(symmetry.lattice)) for i in xrange(symmetry.numk)]
#
# magRefs = [MagReflection() for i in xrange(len(reflections)*2*symmetry.numk)]
# for i in xrange(len(reflections)):
# # loop through all k and -k vectors
# for ii in xrange(-symmetry.numk, symmetry.numk+1):
# if (ii == 0): continue
# index = ii + symmetry.numk
# if (ii > 0): index -= 1
# ref = magRefs[i*2*symmetry.numk+index]
# ref.signk = -copysign(1,ii)
# ref.knum = abs(ii)
# newhkl = np.array(reflections[i].hkl) - \
# ref.signk * np.array(symmetry.k[ref.knum-1])
# ref.hkl = float3(*newhkl)
# ref.s = calcS(cell, ref.hkl)
# ref.multip = reflections[i].multip
# ref.equalMinus = equalMinus[ref.knum-1]
# fn(cell, symmetry, atomList, ref)
# return magRefs
# calcIntensity: calculates the intensity for a given set of reflections,
# based on the structure factor
def calcIntensity(refList, atomList, spaceGroup, wavelength, cell=None,
magnetic=False):
# TODO: be smarter about determining whether the structure is magnetic
# TODO: make sure magnetic phase factor is properly being taken into account
if (refList.magnetic):
sfs = calcMagStructFact(refList, atomList, spaceGroup, cell)
sfs2 = np.array([np.sum(np.array(sf)**2) for sf in sfs])
else:
sfs2 = np.array(calcStructFact(refList, atomList, spaceGroup, wavelength))
multips = np.array([ref.multip for ref in refList])
tt = np.radians(np.array([twoTheta(ref.s, wavelength) for ref in refList]))
lorentz = (1+np.cos(tt)**2) / (np.sin(tt)*np.sin(tt/2))
return sfs2 * multips * lorentz
# makeGaussians() creates a series of Gaussians to represent the powder
# diffractionn pattern
def makeGaussians(reflections, coeffs, I, scale, wavelength):
Gaussian.scaleFactor = scale
gaussians = [Gaussian(twoTheta(rk.s, wavelength),
coeffs[0], coeffs[1], coeffs[2], Ik, rk.hkl)
for rk,Ik in zip(reflections,I)]
return gaussians
# getIntensity: calculates the intensity at a given 2*theta position, or for an
# array of 2*theta positions
def getIntensity(gaussians, background, tt):
#return background(tt) + sum(g(tt) for g in gaussians)
v = background(tt)
for g in gaussians:
g.add(v,tt)
return v
# removeRange: takes in an array of 2*theta intervals and removes them from
# consideration for data analysis, with an optional argument for removing the
# corresponding intensities as well
def removeRange(tt, remove, intensity=None):
if (remove == None):
if (intensity != None): return (tt, intensity)
else: return tt
if (not isSequence(remove[0]) or len(remove[0]) == 1):
# single interval
keepEntries = (tt < remove[0]) | (tt > remove[1])
tt = tt[keepEntries]
if (intensity != None):
intensity = intensity[keepEntries]
return (tt, intensity)
else: return tt
else:
# array of intervals
if (intensity != None):
for interval in remove:
tt, intensity = removeRange(tt, interval, intensity)
return (tt, intensity)
else:
for interval in remove:
tt = removeRange(tt, interval)
return tt
# diffPattern: generates a neutron diffraction pattern from a file containing
# crystallographic information or from the same information generated
# elsewhere
def diffPattern(infoFile=None, backgroundFile=None, wavelength=1.5403,
ttMin=0, ttMax=180, ttStep=0.05, exclusions=None,
spaceGroup=None, cell=None, atomList=None,
symmetry=None, basisSymmetry=None, magAtomList=None,
magnetic=False, info=False, plot=False, saveFile=None):
background = LinSpline(backgroundFile)
sMin, sMax = getS(ttMin, wavelength), getS(ttMax, wavelength)
if magnetic:
if (infoFile != None):
info = readMagInfo(infoFile)
if (spaceGroup == None): spaceGroup = info[0]
if (cell == None): cell = info[1]
if (magAtomList == None): magAtomList = info[2]
if (symmetry == None): symmetry = info[3]
if (basisSymmetry == None): basisSymmetry = symmetry
# magnetic peaks
magRefList = satelliteGen(cell, symmetry, sMax)
print "length of reflection list " + str(len(magRefList))
magIntensities = calcIntensity(magRefList, magAtomList, basisSymmetry,
wavelength, cell, True)
# add in structural peaks
if (atomList == None): atomList = readInfo(infoFile)[2]
refList = hklGen(spaceGroup, cell, sMin, sMax, True)
intensities = calcIntensity(refList, atomList, spaceGroup, wavelength)
reflections = magRefList[:] + refList[:]
intensities = np.append(magIntensities, intensities)
else:
if (infoFile != None):
info = readInfo(infoFile)
if (spaceGroup == None): spaceGroup = info[0]
if (cell == None): cell = info[1]
if (atomList == None): atomList = info[2]
refList = hklGen(spaceGroup, cell, sMin, sMax, True)
reflections = refList[:]
intensities = calcIntensity(refList, atomList, spaceGroup, wavelength)
gaussians = makeGaussians(reflections,[0,0,1], intensities, 1, wavelength)
numPoints = int(floor((ttMax-ttMin)/ttStep)) + 1
tt = np.linspace(ttMin, ttMax, numPoints)
intensity = getIntensity(gaussians, background, tt)
if info:
if magnetic:
printInfo(cell, spaceGroup, (atomList, magAtomList), (refList, magRefList),
wavelength, basisSymmetry)
else:
printInfo(cell, spaceGroup, atomList, refList, wavelength)
if plot:
plotPattern(gaussians, background, None, None, ttMin, ttMax, ttStep,
exclusions, "hkl")
pylab.show()
if saveFile:
np.savetxt(saveFile, (tt, intensity), delimiter=" ")
return (tt, intensity)
# printInfo: prints out information about the provided space group and atoms,
# as well as the generated reflections
def printInfo(cell, spaceGroup, atomLists, refLists, wavelength, symmetry=None):
if (isinstance(refLists, ReflectionList)):
atomLists = (atomLists,)
refLists = (refLists,)
divider = "-" * 40
print "Cell information (%s cell)" % rstrip(spaceGroup.xtalSystem)
print divider
print " a = %.3f alpha = %.3f" % (cell.length[0], cell.angle[0])
print " b = %.3f beta = %.3f" % (cell.length[1], cell.angle[1])
print " c = %.3f gamma = %.3f" % (cell.length[2], cell.angle[2])
print divider
print
print "Space group information"
print divider
print " Number: ", spaceGroup.number
print " H-M Symbol: ", spaceGroup.symbol
print " Hall Symbol: ", spaceGroup.hallSymbol
print " Crystal System: ", spaceGroup.xtalSystem
print " Laue Class: ", spaceGroup.laue
print " Point Group: ", spaceGroup.pointGroup
print " General Multiplicity: ", spaceGroup.multip
print divider
print
print "Atom information (%d atoms)" % len(atomLists[0])
print divider
atomList = atomLists[0]
magnetic = atomList.magnetic
label = [rstrip(atom.label) for atom in atomList]
x, y, z = tuple(["%.3f" % atom.coords[i] for atom in atomList]