--- Jörg-Rüdiger Hill 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_... 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/