[cctbxbb] Problems with conversion from conventional to primitive cell

Fri Nov 25 10:41:15 PST 2005

--- Jörg-Rüdiger Hill <Joerg-Ruediger.Hill at scienomics.com> wrote:

> Thanks for your script. I am doing this in C++. I can't find a
> xray::structure 
> class in the C++ documentation. Does this only exist in Python ?

Yes. C++ is used only for low-level algorithms. It is too expensive to write
everything in C++, too frustrating, too unsafe, and fortunately unnecessary.

> Nevertheless, I run your script in Python and the result is somewhat
> different 
> from what I get, but I still think it is wrong. If I create and visualize a 
> supercell from the coordinates given by the script (e. g., with gdis) I no 
> longer see the hexagonal arrangement of the atoms typical for the 
> conventional cell. Since the conversion from conventional to primitive cell 
> should only redefine the cell but not the atom positions (in Cartesian 
> coordinates) I believe the output is incorrect.

Attached is an expanded script that shows

1. The coordination sequences for the original icsd_structure and the
transformed p1_structure are identical. I.e. the local environment of each atom
in the p1_structure is identical to that of the corresponding atom in the
icsd_structure.

2. f_calc for both structures are identical; note that there is a factor of 3
in the amplitudes resulting from the transformation from the R cell to the
primitive cell.

About 1/3 of the cctbx code is just for internal testing. Here is an example
related to your doubts:

http://phenix-online.org/cctbx_sources/cctbx/cctbx/regression/tst_expand_to_p1.py

It is highly unlikely that there are gross bugs like the one you suspect.

> Since the conversion from conventional to primitive cell 
> should only redefine the cell but not the atom positions (in Cartesian 
> coordinates)

You can convince yourself that this is not the case by working through your
Cr2O3 example. You have to consider three matrices:

  - orthogonalization matrix of the R-centered cell
  - change-of-basis matrix R cell -> primitive cell
  - orthogonalization matrix of the primitive cell

Here is a script that shows the three matrices (in Mathematica format):

from cctbx import crystal
from scitbx import matrix

def demo():
  icsd_symmetry = crystal.symmetry(
    unit_cell="4.961950 4.961950 13.597400 90.000000 90.000000 120.000000",
    space_group_symbol="R -3 C")
  icsd_symmetry.show_summary()
  print
  icsd_ortho_matrix = matrix.sqr(
    icsd_symmetry.unit_cell().orthogonalization_matrix())
  print icsd_ortho_matrix.mathematica_form(
    label="icsd_ortho_matrix:", one_row_per_line=True)
  print
  z2p_op = icsd_symmetry.space_group().z2p_op()
  primitive_symmetry = icsd_symmetry.change_basis(cb_op=z2p_op)
  primitive_symmetry.show_summary()
  print
  primitive_ortho_matrix = matrix.sqr(
    primitive_symmetry.unit_cell().orthogonalization_matrix())
  print primitive_ortho_matrix.mathematica_form(
    label="primitive_ortho_matrix:", one_row_per_line=True)
  print
  assert z2p_op.c().t().num() == (0,0,0)
  z2p_matrix = z2p_op.c().r().as_rational().as_float()
  print z2p_matrix.mathematica_form(
    label="z2p_matrix:", one_row_per_line=True)
  print
  print (primitive_ortho_matrix * z2p_matrix).mathematica_form(
    label="product:", one_row_per_line=True)
  print

if (__name__ == "__main__"):
  demo()

If you run the demo you will see that the "product" is not identical to the
icsd_ortho_matrix. This means the Cartesian coordinates will be different.

You can use Mathematica to recompute the matrix product, or just do it by hand.
The convention for the orthogonalization matrix is shown in the uctbx.h C++
documentation. I believe you can find the z2p_matrix in Int. Vol. A, but I
don't have the volume at hand right now to verify.

I am certain the output of the script in my previous message is correct. Could
you please double-check your input for gdis? Does gdis display your structure
correctly if you simply transform to the rhombohedral basis? You could do the
transformation of the two coordinates and the unit cell parameters by hand. The
space group is simply the rhombohedral setting. I just saw gdis uses sginfo,
therefore I guess it should understand r3c:r.
Once you have verified gids deals correctly with r3c:r, you could expand the
fractional coordinates manually. Keep the unit cell parameters and change the
space group to P1. Once you know this is OK, try to find out what went wrong
converting the output of my script to gdis input.

Cheers,
        Ralf

Here is the expanded version of my previous script, as explained above:

from cctbx import xray
from cctbx import miller
from cctbx import crystal
import cctbx.crystal.coordination_sequences
from cctbx.array_family import flex

def demo():
  """
  Result of ICSD query:
    N * -Cr2O3-[R3-CH] Baster, M.;Bouree, F.;Kowalska, A.;Latacz, Z(2000)
    C 4.961950 4.961950 13.597400 90.000000 90.000000 120.000000
    S GRUP R -3 C
    A Cr1    0.000000 0.000000 0.347570 0.000000
    A O1    0.305830 0.000000 0.250000
  """
  crystal_symmetry = crystal.symmetry(
    unit_cell="4.961950 4.961950 13.597400 90.000000 90.000000 120.000000",
    space_group_symbol="R -3 C")
  scatterers = flex.xray_scatterer()
  scatterers.append(xray.scatterer(
    label="Cr1", site=(0.000000,0.000000,0.347570)))
  scatterers.append(xray.scatterer(
    label="O1", site=(0.305830,0.000000,0.250000)))
  icsd_structure = xray.structure(
    crystal_symmetry=crystal_symmetry,
    scatterers=scatterers)
  icsd_structure.show_summary().show_scatterers()
  print
  icsd_pairs = icsd_structure.show_distances(
    distance_cutoff=2.5, keep_pair_asu_table=True)
  print
  primitive_structure = icsd_structure.primitive_setting()
  primitive_structure.show_summary().show_scatterers()
  print
  p1_structure = primitive_structure.expand_to_p1()
  p1_structure.show_summary().show_scatterers()
  print
  p1_pairs = p1_structure.show_distances(
    distance_cutoff=2.5, keep_pair_asu_table=True)
  print
  for label,structure,pairs in [("ICSD", icsd_structure,icsd_pairs),
                                ("P1", p1_structure,p1_pairs)]:
    print "Coordination sequences for", label, "structure"
    term_table = crystal.coordination_sequences.simple(
      pair_asu_table=pairs.pair_asu_table,
      max_shell=10)
    crystal.coordination_sequences.show_terms(
      structure=structure,
      term_table=term_table)
    print
  icsd_f_calc = icsd_structure.structure_factors(
    d_min=1, algorithm="direct").f_calc()
  icsd_f_calc_in_p1 = icsd_f_calc.primitive_setting().expand_to_p1()
  p1_f_calc = icsd_f_calc_in_p1.structure_factors_from_scatterers(
    xray_structure=p1_structure, algorithm="direct").f_calc()
  for h,i,p in zip(icsd_f_calc_in_p1.indices(),
                   icsd_f_calc_in_p1.data(),
                   p1_f_calc.data()):
    print h, abs(i), abs(p)*3
  print "OK"

if (__name__ == "__main__"):
  demo()

__________________________________ 
Yahoo! Music Unlimited 
Access over 1 million songs. Try it free. 
http://music.yahoo.com/unlimited/