DIALS で積分してみる

2016/05/04: 本記事は歴史的な意味で残してある。現在では、チュートリアル の通りに操作すれば処理できるし、SPring-8 の reverse phi も自動で認識される。

DIALS の積分が、ついに SPring-8 BL41-XU のデータセットで完走できたのでメモ。

データブロックの準備

dials.import コマンドに画像を与えると、ヘッダ情報を分析して datablock.json を作ってくれる。

$ dials.import ../collect10_*.img

--------------------------------------------------------------------------------
DataBlock 0
  format: <class 'dxtbx.format.FormatTIFFRayonix.FormatTIFFRayonix'>
  num images: 180
  num sweeps: 1
  num stills: 0
--------------------------------------------------------------------------------
Writing datablocks to datablock.json

SPring-8 は reverse phi なので、作成された datablock.json を編集。

    "goniometer": [
      {
        "fixed_rotation": [
          1.0, 
          0.0, 
          0.0, 
          0.0, 
          1.0, 
          0.0, 
          0.0, 
          0.0, 
          1.0
        ], 
        "rotation_axis": [
          -1.0, // マイナスをつける
          0.0, 
          0.0
        ]
      }
    ], 

ビームセンタや波長なども確認して、おかしなところがあったら修正する。

スポット検出

全てのフレームを使ってスポット検出が行われる。

$ dials.find_spots datablock.json 

Configuring spot finder from input parameters
--------------------------------------------------------------------------------
Finding strong spots in imageset 0
--------------------------------------------------------------------------------

Finding spots in image 0 to 180...
Extracted strong pixels from images.......................................61.25s
Merged 8 pixel lists with 5645051 pixels...................................0.12s
Extracted 173899 spots.....................................................7.12s
Calculated 173899 spot centroids...........................................0.63s
Calculated 173899 spot intensities.........................................0.08s
Filtered 122558 spots by number of pixels..................................0.16s
Filtered 120609 spots by peak-centroid distance............................0.08s

--------------------------------------------------------------------------------
Saved 120609 reflections to strong.pickle..................................0.44s
Total time:  72.0769660473

指数付け

デフォルトでは、三次元FFT を利用した指数付けが行われる。method= 引数によって fft1d (DPS アルゴリズム)も指定可能。

$ dials.index datablock.json strong.pickle 

reference {
  detector = None
  beam = None
}
...中略...
method = *fft3d fft1d real_space_grid_search
multiple_lattice_search {
  cluster_analysis_search = False
  recycle_unindexed_reflections = False
  recycle_unindexed_reflections_cutoff = 0.1
  minimum_angular_separation = 5
  max_lattices = None
  cluster_analysis {
    method = *dbscan hcluster
    hcluster {
      linkage {
        method = *ward
        metric = *euclidean
      }
      cutoff = 15
      cutoff_criterion = *distance inconsistent
    }
    dbscan {
      eps = 0.05
      min_samples = 30
    }
    min_cluster_size = 20
    intersection_union_ratio_cutoff = 0.4
  }
}
output {
  experiments = experiments.json
  reflections = indexed.pickle
}
Found max_cell: 217.0 Angstrom
FFT gridding: (256,256,256)
Number of centroids used: 52423
model 1 (13613 reflections):
Crystal:
    Unit cell: (秘密。まあ、下から計算できるけどw)
    Space group: P 1
    U matrix:  {{-0.9496, -0.2988, -0.0946},
                {-0.3135,  0.9035,  0.2922},
                {-0.0018,  0.3072, -0.9517}}
    B matrix:  {{ 0.0157,  0.0000,  0.0000},
                {-0.0042,  0.0084,  0.0000},
                {-0.0039, -0.0000,  0.0083}}
    A = UB:    {{-0.0133, -0.0025, -0.0008},
                {-0.0098,  0.0076,  0.0024},
                { 0.0024,  0.0026, -0.0079}}


38817 unindexed reflections

################################################################################
Starting refinement (macro-cycle 1)
################################################################################

...中略...

################################################################################
Starting refinement (macro-cycle 3)
################################################################################


Running refinement
------------------
0 1 2 3 4 5 6 7 8

Refinement steps
----------------
------------------------------------------------------------
| Step | Nref | Objective | RMSD_X   | RMSD_Y   | RMSD_Phi |
|      |      |           | (mm)     | (mm)     | (deg)    |
------------------------------------------------------------
| 0    | 9000 | 8342.2    | 0.03331  | 0.034077 | 0.20053  |
| 1    | 9000 | 8192.3    | 0.032854 | 0.033841 | 0.20122  |
| 2    | 9000 | 8181.2    | 0.032869 | 0.033777 | 0.20122  |
| 3    | 9000 | 8165.4    | 0.032874 | 0.033704 | 0.20118  |
| 4    | 9000 | 8141.7    | 0.032852 | 0.033624 | 0.20109  |
| 5    | 9000 | 8117      | 0.032838 | 0.033534 | 0.2009   |
| 6    | 9000 | 8105.2    | 0.032843 | 0.033483 | 0.20069  |
| 7    | 9000 | 8103.6    | 0.032847 | 0.033475 | 0.2006   |
| 8    | 9000 | 8103.6    | 0.032848 | 0.033474 | 0.20058  |
------------------------------------------------------------
RMSD no longer decreasing
Final refined crystal models:
model 1 (118188 reflections):
Crystal:
    Unit cell: (秘密。まあ、下から計算できるけどw)
    Space group: P 1
    U matrix:  {{-0.9510, -0.2965, -0.0881},
                {-0.3093,  0.9079,  0.2831},
                {-0.0039,  0.2965, -0.9550}}
    B matrix:  {{ 0.0157,  0.0000,  0.0000},
                {-0.0042,  0.0084,  0.0000},
                {-0.0042,  0.0000,  0.0084}}
    A = UB:    {{-0.0133, -0.0025, -0.0007},
                {-0.0098,  0.0076,  0.0024},
                { 0.0027,  0.0025, -0.0080}}

結果が、experiments.json と indexed.pickle に保存される。

Bravais lattice へのあてはめ

ここまでは triclinic な primitive lattice を想定しているので、Bravais 格子へのあてはめを行う。

$ dials.refine_bravais_settings experiments.json indexed.pickle 

refinement {
  verbosity = 1
  parameterisation {
    beam {
      fix = all *in_spindle_plane out_spindle_plane *wavelength
      fix_list = None
    }
    crystal {
      fix = all cell orientation
      cell_fix_list = None
      orientation_fix_list = None
      scan_varying = False
      num_intervals = *fixed_width absolute
      interval_width_degrees = 36.0
      absolute_num_intervals = 5
      UB_model_per = reflection *image
    }
    detector {
      panels = *automatic single multiple hierarchical
      hierarchy_level = 0
      fix = all position orientation
      fix_list = None
    }
    sparse = False
  }
  refinery {
    engine = SimpleLBFGS LBFGScurvs GaussNewton *LevMar
    track_step = False
    track_gradient = False
    track_parameter_correlation = False
    log = None
    max_iterations = None
  }
  target {
    rmsd_cutoff = *fraction_of_bin_size absolute
    bin_size_fraction = 0.33333
    absolute_cutoffs = None
    jacobian_max_nref = None
  }
  reflections {
    reflections_per_degree = 50
    minimum_sample_size = 1000
    maximum_sample_size = None
    use_all_reflections = False
    random_seed = 42
    minimum_number_of_reflections = 20
    close_to_spindle_cutoff = 0.1
    do_outlier_rejection = False
    iqr_multiplier = 1.5
  }
}
lepage_max_delta = 5
verbosity = 0
nproc = Auto
--------------------------------------------------------------------------------------------------
Solution Metric fit  rmsd #spots  crystal_system                                 unit_cell  volume
--------------------------------------------------------------------------------------------------
       9  0.1894 dg 0.069   9000   tetragonal tI 秘密
       8  0.1894 dg 0.067   9000 orthorhombic oI 秘密
       7  0.1894 dg 0.067   9000   monoclinic mC 秘密
       6  0.1750 dg 0.067   9000 orthorhombic oF 秘密
       5  0.1750 dg 0.064   9000   monoclinic mC 秘密
       4  0.1505 dg 0.063   9000   monoclinic mC 秘密
       3  0.1279 dg 0.056   9000   monoclinic mC 秘密
       2  0.0768 dg 0.050   9000   monoclinic mC 秘密
       1  0.0000 dg 0.047   9000    triclinic aP 秘密
--------------------------------------------------------------------------------------------------
usr+sys time: 1.20 seconds
wall clock time: 5.82 seconds

各 lattice type の制約を適用した方位行列が bravais_settings_N.json というファイルに保存される。RMSD と Metric を参考にそれらしい番号を選んで、これを元に精密化する。

$ dials.refine bravais_setting_9.json indexed.pickle 
Configuring refiner
Performing refinement

Running refinement
------------------
0 1 2

しかし、精密化が不安定なようで、assertion を吐いて死ぬので諦める。

積分

いよいよ積分だ。デフォルトでは、XDS 的な三次元プロファイル・フィッティングが行われる。

$ dials.integrate experiments.json -r indexed.pickle -o integrated.pickle

Loaded experiments from experiments.json...................................0.02s
Loaded reference spots from indexed.pickle.................................0.31s
Removed reference spots with invalid coordinates, 9000 remaining...........0.18s
Filtered 9000 reflections with zeta > 0.050000.............................0.00s
Calculated E.S.D Beam Divergence...........................................0.29s
Calculated E.S.D Reflecting Range..........................................0.12s
Sigma B: 0.037465
Sigma M: 0.280245
 Prediction type: scan static prediction
Predicted 1201272 reflections..............................................5.93s
Calculated 1201272 bounding boxes..........................................0.37s
Filtered 1161015 reflections by detector mask..............................0.72s
Filtered 1160052 reflections by zeta >= 0.050000...........................0.09s
Matched 8068 reference spots with predicted reflections....................0.97s
Extracted profiles from 180 frames........................................68.89s
Extracting reflections from the following blocks of images:
   0 ->  10
  10 ->  20
  20 ->  30
  30 ->  40
  40 ->  50
  50 ->  60
  60 ->  70
  70 ->  80
  80 ->  90
  90 -> 100
 100 -> 110
 110 -> 120
 120 -> 130
 130 -> 140
 140 -> 150
 150 -> 160
 160 -> 170
 170 -> 180
Extracted 65330 reflections from block 0...................................3.46s
Masked foreground pixels for 65330 reflections.............................5.48s
Calculated 65330 background values.........................................3.02s
Calculated 65330 reflection centroids......................................3.30s
Integrated 65330 reflections...............................................0.36s
Initialised reciprocal space transform.....................................0.24s
Transformed 65330 reflections.............................................33.76s
Computed correlations between profile and ideal............................0.30s
Number of reflections contributing to each profile:
0: (512, 512, 90); 154
1: (1536, 512, 90); 407
2: (2560, 512, 90); 160
3: (512, 1536, 90); 363
4: (1536, 1536, 90); 580
5: (2560, 1536, 90); 360
6: (512, 2560, 90); 135
7: (1536, 2560, 90); 382
8: (2560, 2560, 90); 127
Average Profile:

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 1 1 1 0 0 0 0 
0 0 0 0 1 2 1 0 0 0 0 
0 0 0 0 1 1 1 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 1 2 1 0 0 0 0 
0 0 0 3 7 8 7 3 0 0 0 
0 0 1 7 B B A 7 2 0 0 
0 0 2 8 B C B 8 2 0 0 
0 0 2 7 B B B 7 2 0 0 
0 0 0 3 7 8 7 3 0 0 0 
0 0 0 0 1 2 1 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 1 0 0 0 0 0 
0 0 0 4 9 B 9 4 0 0 0 
0 0 4 C E E E C 4 0 0 
0 0 9 E E F E E A 0 0 
0 1 B E F G F E B 1 0 
0 0 A E E F E E A 0 0 
0 0 4 C E E E C 4 0 0 
0 0 0 4 9 B 9 4 0 0 0 
0 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 1 3 1 0 0 0 0 
0 0 0 9 F G F 8 0 0 0 
0 0 9 H H I H H 9 0 0 
0 2 F H J L J H F 2 0 
0 3 G I L P L I G 3 0 
0 2 F H K M J H F 2 0 
0 0 9 H H I H H 9 0 0 
0 0 1 9 F G F 9 0 0 0 
0 0 0 0 1 3 1 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 2 3 2 0 0 0 0 
0 0 1 A G I G 9 1 0 0 
0 0 A I J K J I A 0 0 
0 2 H K N Q M J H 2 0 
0 4 I K R Z Q K I 4 0 
0 2 H K O R N J H 2 0 
0 0 A J K K J I A 0 0 
0 0 1 A G I G A 1 0 0 
0 0 0 0 2 3 2 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 1 3 1 0 0 0 0 
0 0 1 9 F G F 8 0 0 0 
0 0 9 H H I H H 9 0 0 
0 2 F H J L J H F 2 0 
0 3 G I L P L I G 3 0 
0 2 F I K L J H F 2 0 
0 0 9 H I I H H 9 0 0 
0 0 1 9 F G F 9 0 0 0 
0 0 0 0 1 3 1 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 1 0 0 0 0 0 
0 0 0 4 9 B 9 4 0 0 0 
0 0 4 C E E E C 4 0 0 
0 0 A E F F F E A 0 0 
0 1 B F F G F E B 1 0 
0 0 A E F F F E A 0 0 
0 0 4 C E F E C 4 0 0 
0 0 0 4 9 B 9 4 0 0 0 
0 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 1 2 2 0 0 0 0 
0 0 0 3 7 8 7 3 0 0 0 
0 0 2 7 B B B 7 2 0 0 
0 0 2 9 B C B 8 3 0 0 
0 0 2 7 A B A 7 2 0 0 
0 0 0 3 7 8 7 3 0 0 0 
0 0 0 0 2 2 2 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 1 1 1 0 0 0 0 
0 0 0 0 1 1 1 0 0 0 0 
0 0 0 0 1 1 1 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 

Integrated 65330 reflections...............................................2.39s
Performed LP-correction on 65330 reflections...............................2.03s

Extracted 65277 reflections from block 1...................................3.86s

...中略...
Saved 1160052 reflections to integrated.pickle.............................2.04s

Total time taken:  1121.54160213

かなり時間がかかるうえ、スポットを切り出した 13GB もある shoebox.dat なるファイルを作る。これは削除していい。積分結果は integrated.pickle なので mtz に変換する。

$dials.export_mtz integrated.pickle experiments.json hklout.mtz

Title: from dials.export_mtz
Space group symbol from file: P1
Space group number from file: 1
Space group from matrices: P 1 (No. 1)
Point group symbol from file: 1
Number of batches: 180
Number of crystals: 1
Number of Miller indices: 1160052
Resolution range: 119.442 1.27141
History:
Crystal 1:
  Name: XTAL
  Project: DIALS
  Id: 1
  Unit cell: (秘密)
  Number of datasets: 1
  Dataset 1:
    Name: FROMDIALS
    Id: 1
    Wavelength: 1
    Number of columns: 11
    label         #valid  %valid      min       max type
    H            1160052 100.00%   -42.00     43.00 H: index h,k,l
    K            1160052 100.00%   -94.00     95.00 H: index h,k,l
    L            1160052 100.00%     0.00     90.00 H: index h,k,l
    M_ISYM       1160052 100.00%     1.00      2.00 Y: M/ISYM, packed partial/reject flag and symmetry number
    BATCH        1160052 100.00%     1.00    179.00 B: BATCH number
    I            1160052 100.00% -9254.44 231876.97 J: intensity
    SIGI         1160052 100.00%     0.00    517.14 Q: standard deviation
    FRACTIONCALC 1160052 100.00%     1.00      1.00 R: real
    XDET         1160052 100.00%     0.07   3057.12 R: real
    YDET         1160052 100.00%     0.01   3058.35 R: real
    ROT          1160052 100.00%    29.00    207.16 R: real

あとは、POINTLESS と AIMLESS で処理すればいい。積分を P1 で行ったので警告されるが、今回はしょうがない。

pointless hklin hklout.mtz hklout pointless.mtz | tee pointless.log
aimless hklin pointless.mtz hklout aimless.mtz | tee aimless.log

お疲れさまでした。