diff --git a/README.rst b/README.rst index c8e3a30..d57504c 100644 --- a/README.rst +++ b/README.rst @@ -38,7 +38,16 @@ diffpy.cmi package for doing refinements with structure objects -* LONGER DESCRIPTION HERE +``diffpy.cmistructure`` adapts crystal and molecular structures to the +``diffpy.srfit`` ParameterSet interface, so lattice parameters, atom positions +and displacement parameters can be refined. It supports +``diffpy.structure.Structure``, ``pyobjcryst`` crystals and molecules, and +``cctbx`` crystals. It also generates space group constraints and bond valence +sum restraints. + +This code was moved from ``diffpy.srfit.structure``. If you are updating a +script, see the `migration guide +`_. For more information about the diffpy.cmistructure library, please consult our `online documentation `_. @@ -82,7 +91,7 @@ and run the following :: This package also provides command-line utilities. To check the software has been installed correctly, type :: - diffpy.cmistructure --version + diffpy-cmistructure --version You can also type the following command to verify the installation. :: @@ -91,7 +100,7 @@ You can also type the following command to verify the installation. :: To view the basic usage and available commands, type :: - diffpy.cmistructure -h + diffpy-cmistructure -h Getting Started --------------- diff --git a/docs/examples/data/C60.stru b/docs/examples/data/C60.stru new file mode 100644 index 0000000..14e306e --- /dev/null +++ b/docs/examples/data/C60.stru @@ -0,0 +1,369 @@ +title C60 +format pdffit +scale 1.000000 +sharp 0.000000, 0.000000, 1.000000, 0.000000 +spcgr P1 +cell 100.000000, 100.000000, 100.000000, 90.000000, 90.000000, 90.000000 +dcell 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000 +ncell 1, 1, 1, 60 +atoms +C 0.51224937 0.50000000 0.53351511 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50378526 0.51164984 0.53351511 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49009005 0.50720000 0.53351511 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49009005 0.49280000 0.53351511 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50378526 0.48835016 0.53351511 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.52449874 0.50000000 0.52594458 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.52828401 0.48835016 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.53440869 0.49280000 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.53440869 0.50720000 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.52828401 0.51164984 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50757053 0.52329969 0.52594458 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.51981990 0.52329969 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.51748048 0.53049969 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50378526 0.53494953 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49766058 0.53049969 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48018010 0.51440000 0.52594458 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48396536 0.52604984 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.47639484 0.52604984 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.46793073 0.51440000 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.47027015 0.50720000 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48018010 0.48560000 0.52594458 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.47027015 0.49280000 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.46793073 0.48560000 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.47639484 0.47395016 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48396536 0.47395016 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50757053 0.47670031 0.52594458 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49766058 0.46950031 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50378526 0.46505047 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.51748048 0.46950031 0.50612469 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.51981990 0.47670031 0.51837406 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48018010 0.52329969 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48251952 0.53049969 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49621474 0.53494953 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50233942 0.53049969 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49242947 0.52329969 0.47405542 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.51603464 0.52604984 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.52360516 0.52604984 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.53206927 0.51440000 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.52972985 0.50720000 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.51981990 0.51440000 0.47405542 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.52972985 0.49280000 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.53206927 0.48560000 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.52360516 0.47395016 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.51603464 0.47395016 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.51981990 0.48560000 0.47405542 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50233942 0.46950031 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49621474 0.46505047 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48251952 0.46950031 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48018010 0.47670031 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49242947 0.47670031 0.47405542 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.47171599 0.48835016 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.46559131 0.49280000 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.46559131 0.50720000 0.49387531 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.47171599 0.51164984 0.48162594 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.47550126 0.50000000 0.47405542 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49621474 0.51164984 0.46648489 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50990995 0.50720000 0.46648489 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.50990995 0.49280000 0.46648489 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.49621474 0.48835016 0.46648489 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 +C 0.48775063 0.50000000 0.46648489 1.0000 + 0.00000000 0.00000000 0.00000000 0.0000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 + 0.00000000 0.00000000 0.00000000 diff --git a/docs/examples/npintensity.py b/docs/examples/npintensity.py new file mode 100644 index 0000000..64d73a0 --- /dev/null +++ b/docs/examples/npintensity.py @@ -0,0 +1,553 @@ +#!/usr/bin/env python +######################################################################## +# +# (c) 2009 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.rst for license information. +# +######################################################################## +"""Example of using ProfileGenerators in FitContributions. + +This is an example of building a ProfileGenerator and using it in a +FitContribution in order to fit theoretical intensity data. ProfileGenerators +are used to organize profile calculators that require more information than can +be conveniently passed into a function call. The IntensityGenerator class is +an example of a ProfileGenerator that can be used by a FitContribution to help +generate a the profile. + +Instructions + +Run the example and note the last line of the output. It will be described in +the code. Then read through 'IntensityGenerator' class. This will motivate the +need for the ProfileGenerator class. Next read the 'makeRecipe' code. This +will demonstrate how to use the generator, the structure container needed by +the generator and introduce new operations that can be used in a registered +equation. + +Extensions + +- The IntensityGenerator class uses the 'add_parameter_set' method to associate + the structure adapter (DiffpyStructureParSet) with the generator. Most SrFit + classes have an 'add_parameter_set' class and can store ParameterSet objects. + Grab the phase object from the IntensityGenerator and try to add it to other + objects used in the fit recipe. Create variables from the moved Parameters + rather than from the 'phase' that lives in the IntensityGenerator and see if + everything still refines. +""" + +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy + +from diffpy.cmistructure.diffpyparset import DiffpyStructureParSet +from diffpy.srfit.fitbase import ( + FitContribution, + FitRecipe, + FitResults, + Profile, + ProfileGenerator, +) + +###### +# Example Code + + +class IntensityGenerator(ProfileGenerator): + """A class for calculating intensity using the Debye equation. + + Calculating intensity from a structure is difficult in general. This + class takes a diffpy.structure.Structure instance and from that + generates a theoretical intensity signal. Unlike the example in + gaussianrecipe.py, the intensity generator is not simple. It must + take a structure object and some Parameters, and from that generate + a signal. At the same time, the structure itself (the lattice, atom + positions, thermal parameters, etc.) needs to be refinable. Thus we + define this ProfileGenerator to help us interface which exposes the + Parameters required by the calculation and provides a way for a + FitContribution to perform that calculation. + + The purpose of a ProfileGenerator is to 1) provide a function that + generates a profile signal 2) organize the Parameters required for + the calculation + + This generator wraps the 'iofq' function defined below. Knowledge of + this function is not required for this example. + """ + + def __init__(self, name): + """Define our generator. + + In this example we will keep count of how many times the + calculation gets performed. The 'count' attribute will be used + to store the count. + """ + ProfileGenerator.__init__(self, name) + # Count the calls + self.count = 0 + return + + def setStructure(self, strufile): + """Set the structure used in the calculation. + + strufile -- The name of a structure file. A + diffpy.structure.Structure object will be created from + the file, and that object will be passed to the 'iofq' + function whenever it is called. + + This will create the refinement Parameters using the + DiffpyStructureParSet adapter from diffpy.cmistructure.diffpyparset. + DiffpyStructureParSet is a ParameterSet object that organizes and gives + attribute access to Parameters and ParameterSets adapted from a diffpy + Structure object. The Parameters embedded in the DiffpyStructureParSet + are proxies for attributes of the diffpy.structure.Structure object + that is needed by the 'iofq' function. The Parameters will be + accessible by name under the 'phase' attribute of this generator, and + are organized hierarchically: + + phase + - lattice (retrieved with 'get_lattice') + - a + - b + - c + - alpha + - beta + - gamma + - scatterers (retrieved with 'get_scatterers') + - atom1 (the name depends on the element) + - x + - y + - z + - occ + - U11 + - U22 + - U33 + - U12 + - U13 + - U23 + - Uiso + - etc. + + The diffpy.structure.Structure instance is held within the + DiffpyStructureParSet as the 'structure' attribute. + """ + # Load the structure from file + from diffpy.structure import Structure + + stru = Structure() + stru.read(str(strufile)) + + # Create a ParameterSet designed to interface with + # diffpy.structure.Structure objects that organizes the Parameter + # hierarchy. Note that the DiffpyStructureParSet holds a handle to the + # loaded structure that we use in the __call__ method below. + # + # We pass the diffpy.structure.Structure instance, and give the + # DiffpyStructureParSet the name "phase". + parset = DiffpyStructureParSet("phase", stru) + + # Put this ParameterSet in the ProfileGenerator. + self.add_parameter_set(parset) + + return + + def __call__(self, q): + """Calculate the intensity. + + This ProfileGenerator will be used in a FitContribution that + will be optimized to fit some data. By the time this function + is evaluated, the diffpy.structure.Structure instance has been + updated by the optimizer via the DiffpyStructureParSet defined + in setStructure. Thus, we need only call iofq with the internal + structure object. + """ + self.count += 1 + print("iofq called", self.count) + return iofq(self.phase.structure, q) + + +# End class IntensityGenerator + + +def makeRecipe(strufile, datname): + """Create a recipe that uses the IntensityGenerator. + + This will create a FitContribution that uses the IntensityGenerator, + associate this with a Profile, and use this to define a FitRecipe. + """ + # The Profile + # Create a Profile. This will hold the experimental and calculated signal. + profile = Profile() + + # Load data and add it to the profile + x, y, u = profile.loadtxt(datname) + + # The ProfileGenerator + # Create an IntensityGenerator named "I". This will be the name we use to + # refer to the generator from within the FitContribution equation. We also + # need to load the model structure we're using. + generator = IntensityGenerator("I") + generator.setStructure(strufile) + + # The FitContribution + # Create a FitContribution, that will associate the Profile with the + # ProfileGenerator. The ProfileGenerator will be accessible as an + # attribute of the FitContribution by its name ("I"). We also want to tell + # the FitContribution to name the x-variable of the profile "q", so we can + # use it in equations with this name. + contribution = FitContribution("bucky") + contribution.add_profile_generator(generator) + contribution.set_profile(profile, xname="q") + + # Now we're ready to define the fitting equation for the FitContribution. + # We need to modify the intensity calculation, and we'll do that from + # within the fitting equation for the sake of instruction. We want to + # modify the calculation in three ways. We want to scale it, add a + # polynomial background, and broaden the peaks. + # + # There is added benefit for defining these operations outside of the + # IntensityGenerator. By combining the different parts of the calculation + # within the fitting equation, the time-consuming iofq calculation is only + # performed when a structural Parameter is changed. If only non-structural + # parameters are changed, such as the background and broadening Parameters, + # then then previously computed iofq value will be used to compute the + # contribution equation. The benefit in this is very apparent when + # refining the recipe with the LM optimizer, which only changes two + # variables at a time most of the time. Note in the refinement output how + # many times the residual is calculated, versus how many times iofq is + # called when using the scipyOptimize function. + + # We will define the background as a string. + + bkgdstr = "b0 + b1*q + b2*q**2 + b3*q**3 + b4*q**4 + b5*q**5 + b6*q**6 +\ + b7*q**7 + b8*q**8 + b9*q**9" + + # This creates a callable equation named "bkgd" within the FitContribution, + # and turns the polynomial coefficients into Parameters. + contribution.register_string_function(bkgdstr, "bkgd") + + # We will create the broadening function that we need by creating a python + # function and registering it with the FitContribution. + pi = numpy.pi + exp = numpy.exp + + def gaussian(q, q0, width): + return ( + 1 + / (2 * pi * width**2) ** 0.5 + * exp(-0.5 * ((q - q0) / width) ** 2) + ) + + # This registers the python function and extracts the name and creates + # Parameters from the arguments. + contribution.register_function(gaussian) + + # Center the Gaussian so it is not truncated. + contribution.q0.value = x[len(x) // 2] + + # Now we can incorporate the scale and bkgd into our calculation. We also + # convolve the signal with the Gaussian to broaden it. Recall that we don't + # need to supply arguments to the registered functions unless we want to + # make changes to their input values. + contribution.set_equation("scale * convolve(I, gaussian) + bkgd") + + # Make the FitRecipe and add the FitContribution. + recipe = FitRecipe() + recipe.add_contribution(contribution) + + # Specify which parameters we want to refine. + recipe.add_variable(contribution.b0, 0) + recipe.add_variable(contribution.b1, 0) + recipe.add_variable(contribution.b2, 0) + recipe.add_variable(contribution.b3, 0) + recipe.add_variable(contribution.b4, 0) + recipe.add_variable(contribution.b5, 0) + recipe.add_variable(contribution.b6, 0) + recipe.add_variable(contribution.b7, 0) + recipe.add_variable(contribution.b8, 0) + recipe.add_variable(contribution.b9, 0) + + # We also want to adjust the scale and the convolution width + recipe.add_variable(contribution.scale, 1) + recipe.add_variable(contribution.width, 0.1) + + # We can also refine structural parameters. Here we extract the + # DiffpyStructureParSet from the intensity generator and use the parameters + # like we would any others. + phase = generator.phase + + # We want to allow for isotropic expansion, so we'll constrain the lattice + # parameters to the same value (the lattice is cubic). Note that we + # constrain to the "a" Parameter directly. In previous examples, we + # constrained to a Variable by name. This has the same effect. + lattice = phase.get_lattice() + a = lattice.a + recipe.add_variable(a) + recipe.add_constraint(lattice.b, a) + recipe.add_constraint(lattice.c, a) + # We want to refine the thermal parameters as well. We will add a new + # Variable that we call "Uiso" and constrain the atomic Uiso values to + # this. Note that we don't give Uiso an initial value. The initial value + # will be inferred from the following constraints. + Uiso = recipe.create_new_variable("Uiso") + for atom in phase.get_scatterers(): + recipe.add_constraint(atom.Uiso, Uiso) + + # Give the recipe away so it can be used! + return recipe + + +def main(): + + # Make the data and the recipe + strufile = Path(__file__).parent / "data/C60.stru" + q = numpy.arange(1, 20, 0.05) + makeData(strufile, q, "C60.iq", 1.0, 100.68, 0.005, 0.13, 2) + + # Make the recipe + recipe = makeRecipe(strufile, "C60.iq") + + # Optimize + scipyOptimize(recipe) + + # Generate and print the FitResults + res = FitResults(recipe) + # We want to see how much speed-up we get from bringing the scale and + # background outside of the intensity generator. Get the number of calls + # to the residual function from the FitRecipe, and the number of calls to + # 'iofq' from the IntensityGenerator. + rescount = recipe.fithooks[0].count + calcount = recipe.bucky.I.count + footer = "iofq called %i%% of the time" % int(100.0 * calcount / rescount) + res.print_results(footer=footer) + + # Plot! + plot_results(recipe) + + return + + +def plot_results(recipe): + """Plot the results contained within a refined FitRecipe.""" + # The background is not part of the standard observed/fit/diff plot + # that plot_recipe produces, so we overlay it afterwards. + q = recipe.bucky.profile.x + bkgd = recipe.bucky.evaluate_equation("bkgd") + + fig, ax = recipe.plot_recipe( + show=False, + return_fig=True, + data_color="b", + fit_color="r", + diff_color="g", + data_label="I(Q) Data", + fit_label="I(Q) Fit", + diff_label="I(Q) diff", + xlabel=r"$Q (\AA^{-1})$", + ylabel="Intensity (arb. units)", + ) + ax.plot(q, bkgd, "c-", label="Bkgd. Fit") + ax.legend(loc=1) + + plt.show() + return + + +def iofq(S, q): + """Calculate I(Q) (X-ray) using the Debye Equation. + + I(Q) = 2 sum(i,j) f_i(Q) f_j(Q) sinc(rij Q) exp(-0.5 ssij Q**2) + (The exponential term is the Debye-Waller factor.) + + S -- A diffpy.structure.Structure instance. It is assumed that the + structure is that of an isolated scatterer. Periodic boundary + conditions are not applied. + q -- The q-points to calculate over. + + This uses cctbx for the calculation of the f_i if it is available, + otherwise f_i = 1. + """ + # The functions we need + sinc = numpy.sinc + exp = numpy.exp + pi = numpy.pi + + # The brute-force calculation is very slow. Thus we optimize a little bit. + + # The precision of distance measurements + deltad = 1e-6 + dmult = int(1 / deltad) + deltau = deltad**2 + umult = int(1 / deltau) + + pairdict = {} + elcount = {} + n = len(S) + for i in range(n): + + # count the number of each element + eli = S[i].element + m = elcount.get(eli, 0) + elcount[eli] = m + 1 + + for j in range(i + 1, n): + + elj = S[j].element + + # Get the pair + els = [eli, elj] + els.sort() + + # Get the distance to the desired precision + d = S.distance(i, j) + D = int(d * dmult) + + # Get the DW factor to the same precision + ss = S[i].Uisoequiv + S[j].Uisoequiv + SS = int(ss * umult) + + # Record the multiplicity of this pair + key = (els[0], els[1], D, SS) + mult = pairdict.get(key, 0) + pairdict[key] = mult + 1 + + # Now we can calculate IofQ from the pair dictionary. Making the dictionary + # first reduces the amount of calls to sinc and exp we have to make. + + # First we must cache the scattering factors + fdict = {} + for el in elcount: + fdict[el] = getXScatteringFactor(el, q) + + # Now we can compute I(Q) for the i != j pairs + y = 0 + x = q * deltad / pi + for key, mult in pairdict.items(): + eli = key[0] + elj = key[1] + fi = fdict[eli] + fj = fdict[elj] + D = key[2] + SS = key[3] + + # Note that numpy's sinc(x) = sin(x*pi)/(x*pi) + y += fi * fj * mult * sinc(x * D) * exp(-0.5 * SS * deltau * q**2) + + # We must multiply by 2 since we only counted j > i pairs. + y *= 2 + + # Now we must add in the i == j pairs. + for el, f in fdict.items(): + y += f**2 * elcount[el] + + # And that's it! + + return y + + +def getXScatteringFactor(el, q): + """Get the x-ray scattering factor for an element over the q range. + + If cctbx is not available, f(q) = 1 is used. + """ + try: + import cctbx.eltbx.xray_scattering as xray + + wk1995 = xray.wk1995(el) + g = wk1995.fetch() + # at_stol - at sin(theta)/lambda = Q/(4*pi) + f = numpy.asarray(map(g.at_stol, q / (4 * numpy.pi))) + return f + except ImportError: + return 1 + + +def makeData(strufile, q, datname, scale, a, Uiso, sig, bkgc, nl=1): + """Make some fake data and save it to file. + + Make some data to fit. This uses iofq to calculate an intensity curve, and + adds to it a background, broadens the peaks, and noise. + + strufile-- A filename holding the sample structure + q -- The q-range to calculate over. + datname -- The name of the file we're saving to. + scale -- The scale factor + a -- The lattice constant to use + Uiso -- The thermal factor for all atoms + sig -- The broadening factor + bkgc -- A parameter that gives minor control of the background. + nl -- Noise level (0, inf), default 1, larger -> less noise. + """ + from diffpy.structure import Structure + + S = Structure() + S.read(str(strufile)) + + # Set the lattice parameters + S.lattice.set_latt_parms(a, a, a) + + # Set a DW factor + for a in S: + a.Uisoequiv = Uiso + y = iofq(S, q) + + # We want to broaden the peaks as well. This simulates instrument effects. + q0 = q[len(q) // 2] + g = numpy.exp(-0.5 * ((q - q0) / sig) ** 2) + y = numpy.convolve(y, g, mode="same") / sum(g) + + # Add a polynomial background. + bkgd = (q + bkgc) ** 2 * (1.5 * max(q) - q) ** 5 + bkgd *= 0.2 * max(y) / max(bkgd) + + y += bkgd + + # Multiply by a scale factor + y *= scale + + # Calculate the uncertainty + u = (y / nl) ** 0.5 + + # And apply the noise + if nl > 0: + y = numpy.random.poisson(y * nl) / nl + + # Now save it + numpy.savetxt(datname, numpy.transpose([q, y, u])) + return + + +def scipyOptimize(recipe): + """Optimize the recipe created above using scipy. + + The FitRecipe we created in makeRecipe has a 'residual' method that + we can be minimized using a scipy optimizer. The details are + described in the source. + """ + # We're going to use the least-squares (Levenberg-Marquardt) optimizer from + # scipy. We simply have to give it the function to minimize + # (recipe.residual) and the starting values of the Variables + # (recipe.get_values()). + from scipy.optimize import leastsq + + print("Fit using scipy's LM optimizer") + leastsq(recipe.residual, recipe.get_values()) + + return + + +if __name__ == "__main__": + + main() + +# End of file diff --git a/docs/examples/npintensityII.py b/docs/examples/npintensityII.py new file mode 100644 index 0000000..2987402 --- /dev/null +++ b/docs/examples/npintensityII.py @@ -0,0 +1,273 @@ +#!/usr/bin/env python +######################################################################## +# +# (c) 2009 The Trustees of Columbia University in the City of New York. +# (c) 2026 Contributors to diffpy.cmistructure. +# All rights reserved. +# +# File coded by: Chris Farrow and members of the diffpy community. +# +# Originally developed in diffpy.srfit by the DANSE Diffraction group and +# Simon J. L. Billinge. +# +# See GitHub contributions for a more detailed list of contributors. +# https://github.com/diffpy/diffpy.cmistructure/graphs/contributors +# +# See LICENSE.rst and LICENSE_DANSE.rst for license information. +# +######################################################################## +"""Example of extracting information from multiple data sets +simultaneously. + +This example builds on npintensitygenerator.py, and uses IntensityGenerator +from that example to build a recipe that simultaneously refines two data sets +generated from the same structure. + +Instructions + +Run the example and then read through the 'makeRecipe' code. You will see how +to refine a single structure to two data sets. + +Extensions + +- In 'makeRecipe' the fit contributions are identically configured except for + the profile. Factor out that configuration code and apply it in the + 'makeRecipe' method. This will reduce the amount of code required to get the + job done, make it clearer what is being done and therefore reduce potential + mistakes in the code. This encapsulation of configuration workflow is the + first step towards writing a user interface. +""" + +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy +from npintensity import IntensityGenerator, makeData, scipyOptimize + +from diffpy.srfit.fitbase import ( + FitContribution, + FitRecipe, + FitResults, + Profile, +) + +###### +# Example Code + + +def makeRecipe(strufile, datname1, datname2): + """Create a recipe that uses the IntensityGenerator. + + We will create two FitContributions that use the IntensityGenerator + from npintensitygenerator.py and associate each of these with a + Profile, and use this to define a FitRecipe. + + Both simulated data sets come from the same structure. We're going + to make two FitContributions that are identical, except for the + profile that is held in each. We're going to assure that the + structures are identical by using the same DiffpyStructureParSet + (which is generated by the IntensityGenerator when we load the + structure) in both generators. + """ + # The Profiles + # Create two Profiles for the two FitContributions. + profile1 = Profile() + profile2 = Profile() + + # Load data into the Profiles + profile1.loadtxt(datname1) + x, y, u = profile2.loadtxt(datname2) + + # The ProfileGenerators + # Create two IntensityGenerators named "I". There will not be a name + # conflict, since the name is only meaningful within the FitContribution + # that holds the ProfileGenerator. Load the structure into one and make + # sure that the second ProfileGenerator is using the same + # DiffyStructureParSet. This will assure that both ProfileGenerators are + # using the exact same Parameters, and underlying Structure object in the + # calculation of the profile. + generator1 = IntensityGenerator("I") + generator1.setStructure(strufile) + generator2 = IntensityGenerator("I") + generator2.add_parameter_set(generator1.phase) + + # The FitContributions + # Create the FitContributions. + contribution1 = FitContribution("bucky1") + contribution1.add_profile_generator(generator1) + contribution1.set_profile(profile1, xname="q") + contribution2 = FitContribution("bucky2") + contribution2.add_profile_generator(generator2) + contribution2.set_profile(profile2, xname="q") + + # Now we're ready to define the fitting equation for each FitContribution. + # The functions registered below will be independent, even though they take + # the same form and use the same Parameter names. By default, Parameters + # in different contributions are different Parameters even if they have the + # same names. FitContributions are isolated namespaces than only share + # information if you tell them to by using addParameter or + # add_parameter_set. + bkgdstr = "b0 + b1*q + b2*q**2 + b3*q**3 + b4*q**4 + b5*q**5 + b6*q**6 +\ + b7*q**7 +b8*q**8 + b9*q**9" + + contribution1.register_string_function(bkgdstr, "bkgd") + contribution2.register_string_function(bkgdstr, "bkgd") + + # We will create the broadening function by registering a python function. + pi = numpy.pi + exp = numpy.exp + + def gaussian(q, q0, width): + return ( + 1 + / (2 * pi * width**2) ** 0.5 + * exp(-0.5 * ((q - q0) / width) ** 2) + ) + + contribution1.register_function(gaussian) + contribution2.register_function(gaussian) + # Center the gaussian + contribution1.q0.value = x[len(x) // 2] + contribution2.q0.value = x[len(x) // 2] + + # Now we can incorporate the scale and bkgd into our calculation. We also + # convolve the signal with the gaussian to broaden it. + contribution1.set_equation("scale * convolve(I, gaussian) + bkgd") + contribution2.set_equation("scale * convolve(I, gaussian) + bkgd") + + # Make a FitRecipe and associate the FitContributions. + recipe = FitRecipe() + recipe.add_contribution(contribution1) + recipe.add_contribution(contribution2) + + # Specify which Parameters we want to refine. We want to refine the + # background that we just defined in the FitContributions. We have to do + # this separately for each FitContribution. We tag the variables so it is + # easy to retrieve the background variables. + recipe.add_variable(contribution1.b0, 0, name="b1_0", tag="bcoeffs1") + recipe.add_variable(contribution1.b1, 0, name="b1_1", tag="bcoeffs1") + recipe.add_variable(contribution1.b2, 0, name="b1_2", tag="bcoeffs1") + recipe.add_variable(contribution1.b3, 0, name="b1_3", tag="bcoeffs1") + recipe.add_variable(contribution1.b4, 0, name="b1_4", tag="bcoeffs1") + recipe.add_variable(contribution1.b5, 0, name="b1_5", tag="bcoeffs1") + recipe.add_variable(contribution1.b6, 0, name="b1_6", tag="bcoeffs1") + recipe.add_variable(contribution1.b7, 0, name="b1_7", tag="bcoeffs1") + recipe.add_variable(contribution1.b8, 0, name="b1_8", tag="bcoeffs1") + recipe.add_variable(contribution1.b9, 0, name="b1_9", tag="bcoeffs1") + recipe.add_variable(contribution2.b0, 0, name="b2_0", tag="bcoeffs2") + recipe.add_variable(contribution2.b1, 0, name="b2_1", tag="bcoeffs2") + recipe.add_variable(contribution2.b2, 0, name="b2_2", tag="bcoeffs2") + recipe.add_variable(contribution2.b3, 0, name="b2_3", tag="bcoeffs2") + recipe.add_variable(contribution2.b4, 0, name="b2_4", tag="bcoeffs2") + recipe.add_variable(contribution2.b5, 0, name="b2_5", tag="bcoeffs2") + recipe.add_variable(contribution2.b6, 0, name="b2_6", tag="bcoeffs2") + recipe.add_variable(contribution2.b7, 0, name="b2_7", tag="bcoeffs2") + recipe.add_variable(contribution2.b8, 0, name="b2_8", tag="bcoeffs2") + recipe.add_variable(contribution2.b9, 0, name="b2_9", tag="bcoeffs2") + + # We also want to adjust the scale and the convolution width + recipe.add_variable(contribution1.scale, 1, name="scale1") + recipe.add_variable(contribution1.width, 0.1, name="width1") + recipe.add_variable(contribution2.scale, 1, name="scale2") + recipe.add_variable(contribution2.width, 0.1, name="width2") + + # We can also refine structural parameters. We only have to do this once, + # since each generator holds the same DiffpyStructureParSet. + phase = generator1.phase + lattice = phase.get_lattice() + a = recipe.add_variable(lattice.a) + # We want to allow for isotropic expansion, so we'll make constraints for + # that. + recipe.add_constraint(lattice.b, a) + recipe.add_constraint(lattice.c, a) + # We want to refine the thermal parameters as well. We will add a new + # variable that we call "Uiso" and constrain the atomic Uiso values to + # this. Note that we don't give Uiso an initial value. The initial value + # will be inferred from the subsequent constraints. + Uiso = recipe.create_new_variable("Uiso") + for atom in phase.get_scatterers(): + recipe.add_constraint(atom.Uiso, Uiso) + + # Give the recipe away so it can be used! + return recipe + + +def plot_results(recipe): + """Plot the results contained within a refined FitRecipe. + + The recipe has two contributions ("bucky1" and "bucky2"), so + plot_recipe produces one figure per contribution. The backgrounds + are not part of the standard observed/fit/diff plot, so they are + overlaid on each figure afterwards. + """ + q = recipe.bucky1.profile.x + bkgd1 = recipe.bucky1.evaluate_equation("bkgd") + bkgd2 = recipe.bucky2.evaluate_equation("bkgd") + + figs, axes = recipe.plot_recipe( + show=False, + return_fig=True, + data_label="I(Q) Data", + fit_label="I(Q) Fit", + diff_label="I(Q) diff", + xlabel=r"$Q (\AA^{-1})$", + ylabel="Intensity (arb. units)", + ) + # "bucky1" was added to the recipe first, so its axes come first. + axes[0].plot(q, bkgd1, "c-", label="Bkgd1 Fit") + axes[0].legend(loc=1) + axes[1].plot(q, bkgd2, "c-", label="Bkgd2 Fit") + axes[1].legend(loc=1) + + plt.show() + return + + +def main(): + + # Make two different data sets, each from the same structure, but with + # different scale, noise, broadening and background. + strufile = Path(__file__).parent / "data/C60.stru" + q = numpy.arange(1, 20, 0.05) + makeData(strufile, q, "C60_1.iq", 8.1, 101.68, 0.008, 0.12, 2, 0.01) + makeData(strufile, q, "C60_2.iq", 3.2, 101.68, 0.02, 0.003, 0, 1) + + # Make the recipe + recipe = makeRecipe(strufile, "C60_1.iq", "C60_2.iq") + + # Optimize + # Since the backgrounds have a large effect on the profile, we will refine + # them first, but do so separately. + # To refine the background from the first contribution, we will fix + # all other parameters and give the second contribution no weight in the + # fit. + recipe.fix("all") + recipe.free("bcoeffs1") + recipe.set_weight(recipe.bucky2, 0) + scipyOptimize(recipe) + # Now do the same for the second background + recipe.fix("all") + recipe.free("bcoeffs1") + recipe.set_weight(recipe.bucky2, 1) + recipe.set_weight(recipe.bucky1, 0) + scipyOptimize(recipe) + # Now refine everything with the structure parameters included + recipe.free("all") + recipe.set_weight(recipe.bucky1, 1) + scipyOptimize(recipe) + + # Generate and print the FitResults + res = FitResults(recipe) + res.print_results() + + # Plot! + plot_results(recipe) + + return + + +if __name__ == "__main__": + + main() + +# End of file diff --git a/docs/source/getting-started.rst b/docs/source/getting-started.rst deleted file mode 100644 index 6b62bf8..0000000 --- a/docs/source/getting-started.rst +++ /dev/null @@ -1,79 +0,0 @@ -:tocdepth: -1 - -.. index:: getting-started - -.. _getting-started: - -================ -Getting started -================ - -Here are some example templates provided to help you get started with writing your documentation. You can use these templates to create your own documentation. - -Reuse ``.rst`` files across multiple pages ------------------------------------------- - -Here is how you can reuse a reusable block of ``.rst`` files across multiple pages: - -.. include:: snippets/example-table.rst - -.. warning:: - - Ensure that the ``.rst`` file you are including is not too long. If it is too long, it may be better to split it into multiple files and include them separately. - -Refer to a specific section in the documentation ------------------------------------------------- - -You can use the ``ref`` tag to refer to a specific section in the documentation. For example, you can refer to the section below using the ``:ref:`` tag as shown :ref:`here `. - -.. note:: - - Please check the raw ``.rst`` file of this page to see the exact use of the ``:ref:`` tag. - -Embed your code snippets in the documentation ---------------------------------------------- - -Here is how you can write a block of code in the documentation. You can use the ``code-block`` directive to write a block of code in the documentation. For example, you can write a block of code as shown below: - -.. code-block:: bash - - # Create a new environment, without build dependencies (pure Python package) - conda create -n -env python=3.14 \ - --file requirements/tests.txt \ - --file requirements/conda.txt - - # Create a new environment, with build dependencies (non-pure Python package) - conda create -n -env python=3.14 \ - --file requirements/tests.txt \ - --file requirements/conda.txt \ - --file requirements/build.txt - - # Activate the environment - conda activate _env - - # Install your package locally - # `--no-deps` to NOT install packages again from `requirements.pip.txt` - pip install -e . --no-deps - - # Run pytest locally - pytest - - # ... run example tutorials - -.. _attach-image: - -Attach an image to the documentation ------------------------------------- - -Here is how you attach an image to the documentation. The ``/docs/source/img/scikit-package-logo-text.png`` example image is provided in the template. - -.. image:: ./img/scikit-package-logo-text.png - :alt: codecov-in-pr-comment - :width: 400px - :align: center - - -Other useful directives ------------------------ - -Here is how you can do menu selection :menuselection:`Admin --> Settings` and display labels for buttons like :guilabel:`Privacy level`. diff --git a/docs/source/img/.placeholder b/docs/source/img/.placeholder deleted file mode 100644 index e69de29..0000000 diff --git a/docs/source/img/scikit-package-logo-text.png b/docs/source/img/scikit-package-logo-text.png deleted file mode 100644 index 823178d..0000000 Binary files a/docs/source/img/scikit-package-logo-text.png and /dev/null differ diff --git a/docs/source/index.rst b/docs/source/index.rst index 441b80d..f61aed3 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -9,14 +9,6 @@ | Software version |release| | Last updated |today|. -=============== -Getting started -=============== - -Welcome to the ``diffpy.cmistructure`` documentation! - -To get started, please visit the :ref:`Getting started ` page. - ======= Authors ======= @@ -43,7 +35,8 @@ Table of contents .. toctree:: :maxdepth: 2 - getting-started + migrating-from-srfit + making-a-structure-adapter Package API release license diff --git a/docs/source/making-a-structure-adapter.rst b/docs/source/making-a-structure-adapter.rst new file mode 100644 index 0000000..79e055f --- /dev/null +++ b/docs/source/making-a-structure-adapter.rst @@ -0,0 +1,379 @@ +:tocdepth: -1 + +.. index:: making-a-structure-adapter + +.. _making-a-structure-adapter: + +========================== +Making a structure adapter +========================== + +A structure adapter is a ``ParameterSet`` that wraps a structure object, so +that the lattice parameters, atomic positions and displacement parameters of +that object become refinable srfit ``Parameter`` objects. When a ``Parameter`` +changes during a refinement, the adapter writes the new value straight into +the wrapped structure object. + +``diffpy.cmistructure`` ships adapters for these structure representations: + +.. list-table:: + :header-rows: 1 + + * - Structure object + - Adapter + * - ``diffpy.structure.Structure`` + - ``diffpy.cmistructure.diffpyparset.DiffpyStructureParSet`` + * - ``pyobjcryst.crystal.Crystal`` + - ``diffpy.cmistructure.objcrystparset.ObjCrystCrystalParSet`` + * - ``pyobjcryst.molecule.Molecule`` + - ``diffpy.cmistructure.objcrystparset.ObjCrystMoleculeParSet`` + * - ``cctbx.crystal.special_position_settings`` + - ``diffpy.cmistructure.cctbxparset.CCTBXCrystalParSet`` + +If your structures live in some other representation, for example a class +from another package or one written for your own project, write an adapter +for it. The adapter then works with ``constrain_as_space_group`` and, if it +derives from ``SrRealParSet``, with the PDF generators in ``diffpy.cmipdf``. + +What an adapter must provide +---------------------------- + +An adapter is a subclass of one of two base classes: + +``BaseStructureParSet`` + The minimal base class. Use it when the adapter only needs to work with + space group constraints. + +``SrRealParSet`` + A subclass of ``BaseStructureParSet`` that adds what the + ``diffpy.srreal`` calculators need: ``use_symmetry()``, + ``using_symmetry()``, ``restrain_bvs()`` and ``_get_srreal_structure()``. + Use it when the structure will be used to calculate a PDF or other + srreal quantity. This is the usual choice. + +The rest of the package relies on the following parts of an adapter: + +.. list-table:: + :header-rows: 1 + :widths: 30 70 + + * - Member + - Requirement + * - ``structure`` + - The wrapped structure object. + * - ``can_adapt(structure)`` + - A classmethod that returns True if ``structure`` is of the type this + adapter wraps, and False otherwise. It must not raise if an optional + package is missing; return False instead. + * - ``get_lattice()`` + - Return a ``ParameterSet`` holding ``Parameter`` objects named ``a``, + ``b``, ``c``, ``alpha``, ``beta`` and ``gamma``. That ``ParameterSet`` + must have an ``angle_units`` attribute set to ``"deg"`` or ``"rad"``. + * - ``get_scatterers()`` + - Return a list of ``ParameterSet`` objects, one per atom or scatterer. + Each holds the fractional coordinates as ``Parameter`` objects named + ``x``, ``y`` and ``z``, plus its displacement parameters (see + :ref:`adapter-adps`). + * - ``_get_srreal_structure()`` + - ``SrRealParSet`` only. Return an object that the ``diffpy.srreal`` + calculators can read (see :ref:`adapter-srreal`). + +Example: adapting a simple crystal class +---------------------------------------- + +The example below adapts a minimal crystal representation, ``SimpleCrystal``, +which stores the unit cell as a list of six numbers and the atoms as a list +of ``SimpleSite`` objects. Swap these two classes for the structure classes +you actually use. + +.. code-block:: python + + from operator import getitem, setitem + + from diffpy.cmistructure.srrealparset import SrRealParSet + from diffpy.srfit.fitbase.parameter import ParameterAdapter + from diffpy.srfit.fitbase.parameterset import ParameterSet + from diffpy.srfit.util.argbinders import bind2nd + + + # The structure representation being adapted. + class SimpleSite: + def __init__(self, element, xyz, uiso): + self.element = element + self.xyz = list(xyz) + self.uiso = uiso + + + class SimpleCrystal: + def __init__(self, cell, sites): + self.cell = list(cell) + self.sites = list(sites) + + + # The adapter. + def _add_list_item_parameters(parameter_set, names, values): + for index, name in enumerate(names): + parameter_set.add_parameter( + ParameterAdapter( + name, + values, + bind2nd(getitem, index), + bind2nd(setitem, index), + ) + ) + + + class SimpleCellParSet(ParameterSet): + def __init__(self, cell): + ParameterSet.__init__(self, "lattice") + self.angle_units = "deg" + _add_list_item_parameters( + self, ("a", "b", "c", "alpha", "beta", "gamma"), cell + ) + + + class SimpleSiteParSet(ParameterSet): + def __init__(self, name, site): + ParameterSet.__init__(self, name) + self.site = site + _add_list_item_parameters(self, ("x", "y", "z"), site.xyz) + self.add_parameter(ParameterAdapter("Uiso", site, attr="uiso")) + + @property + def element(self): + return self.site.element + + + class SimpleCrystalParSet(SrRealParSet): + def __init__(self, name, structure): + SrRealParSet.__init__(self, name) + self.structure = structure + self.add_parameter_set(SimpleCellParSet(structure.cell)) + self.sites = [] + element_counts = {} + for site in structure.sites: + index = element_counts.get(site.element, 0) + element_counts[site.element] = index + 1 + site_parameter_set = SimpleSiteParSet( + f"{site.element}{index}", site + ) + self.add_parameter_set(site_parameter_set) + self.sites.append(site_parameter_set) + + @classmethod + def can_adapt(cls, structure): + return isinstance(structure, SimpleCrystal) + + def get_lattice(self): + return self.lattice + + def get_scatterers(self): + return self.sites + + def _get_srreal_structure(self): + from diffpy.srreal.structureadapter import nosymmetry + from diffpy.structure import Atom, Lattice, Structure + + atoms = [ + Atom(site.element, site.xyz, Uisoequiv=site.uiso) + for site in self.structure.sites + ] + lattice = Lattice(*self.structure.cell) + diffpy_structure = Structure(atoms, lattice=lattice) + if self.using_symmetry(): + return diffpy_structure + return nosymmetry(diffpy_structure) + +The sections below explain each part. + +Wrap each value with ``ParameterAdapter`` +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +A ``ParameterAdapter`` is a ``Parameter`` that does not store its value. +Instead it reads and writes the value in the wrapped object every time, so +the structure object and the ``Parameter`` always agree. There are two ways +to create one: + +- If the value is an attribute of an object, pass the attribute name, as in + ``ParameterAdapter("Uiso", site, attr="uiso")``. +- Otherwise, pass a getter and a setter. The getter is called as + ``getter(obj)`` and returns the value. The setter is called as + ``setter(obj, value)``. In the example, ``bind2nd(getitem, index)`` and + ``bind2nd(setitem, index)`` read and write one item of a list. + +Keep a reference to the live object, never a copy of its values. If you +wrap a copy, the refinement changes the copy and the structure object never +sees the new values. + +The lattice +^^^^^^^^^^^ + +``get_lattice()`` returns a ``ParameterSet`` with ``Parameter`` objects +named ``a``, ``b``, ``c``, ``alpha``, ``beta`` and ``gamma``. +``constrain_as_space_group`` looks the lattice parameters up by these names +and fixes the angles to 90 or 120 for the higher-symmetry crystal systems. +It reads the ``angle_units`` attribute to know whether to use degrees +(``"deg"``) or radians (``"rad"``). + +The ``ParameterSet`` in the example is named ``"lattice"``, so the adapter +can return it as ``self.lattice``. + +The scatterers +^^^^^^^^^^^^^^ + +``get_scatterers()`` returns one ``ParameterSet`` per atom, in a fixed +order. ``constrain_as_space_group`` reads the positions from the ``x``, +``y`` and ``z`` parameters of each one, so these must be fractional +coordinates. + +Each scatterer is added to the adapter with ``add_parameter_set``, which +makes it reachable as an attribute of the adapter, for example +``phase.Ni0``. The names must therefore be valid Python identifiers and +unique within the structure. The built-in adapters use the element symbol +followed by a running index for that element, starting from zero: ``Ni0``, +``Ni1``, ``O0`` and so on. + +An ``element`` property is not required by ``constrain_as_space_group``, +but the built-in adapters provide one, and users expect it. + +.. _adapter-adps: + +Displacement parameters +^^^^^^^^^^^^^^^^^^^^^^^ + +``constrain_as_space_group`` constrains the atomic displacement parameters +(ADPs) only for scatterers that have all six anisotropic parameters. By +default it looks for parameters named ``U11``, ``U22``, ``U33``, ``U12``, +``U13`` and ``U23``. It then uses the isotropic parameter, by default +``Uiso``, for sites whose symmetry makes them isotropic. It leaves +scatterers that lack any of the six anisotropic parameters out of the ADP +constraints, without raising an error. + +The ``SimpleSite`` class in the example has only an isotropic ``uiso``. The +space group constraints therefore leave the ADPs alone, and you constrain +``Uiso`` yourself in the fit recipe. If your structure stores anisotropic +ADPs, add all six as parameters. ``DiffpyAtomParSet`` in +``diffpy.cmistructure.diffpyparset`` is a complete example. + +If your structure uses different names, for example B-factors, pass them +to ``constrain_as_space_group`` with the ``adpsymbols`` and ``isosymbol`` +arguments instead of renaming the parameters: + +.. code-block:: python + + constrain_as_space_group( + phase, + "Fm-3m", + adpsymbols=["B11", "B22", "B33", "B12", "B13", "B23"], + isosymbol="Biso", + ) + +.. _adapter-srreal: + +The structure seen by srreal calculators +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +PDF generators do not read the adapter's ``Parameter`` objects. Before every +calculation they call ``_get_srreal_structure()`` and pass the result to a +``diffpy.srreal`` calculator. The default in ``SrRealParSet`` returns +``self.structure``, which works only if srreal already understands that +type, as it does for ``diffpy.structure.Structure`` and the ``pyobjcryst`` +classes. + +For any other type you have two options: + +- Convert the structure to one that srreal understands, as the example does + by building a ``diffpy.structure.Structure`` from the current values. This + is the simplest option. The conversion runs on every calculation, which is + fast enough for small structures. +- Register a converter for your type with + ``diffpy.srreal.structureadapter.RegisterStructureAdapter``, so that + srreal can read the structure directly. See the ``diffpy.srreal`` + documentation for details. + +In both cases, wrap the result with +``diffpy.srreal.structureadapter.nosymmetry`` when ``self.using_symmetry()`` +is False, so that the generator's ``periodic`` setting is respected. + +Use the adapter +--------------- + +Create the adapter from a structure object and use it like a built-in +adapter. Here the example adapter is used to constrain a nickel unit cell, +given in P1 with four atoms, to the space group Fm-3m: + +.. code-block:: python + + from diffpy.cmistructure import constrain_as_space_group + + positions = [(0, 0, 0), (0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0)] + crystal = SimpleCrystal( + [3.52, 3.52, 3.52, 90, 90, 90], + [SimpleSite("Ni", xyz, 0.005) for xyz in positions], + ) + phase = SimpleCrystalParSet("nickel", crystal) + + space_group_parameters = constrain_as_space_group(phase, "Fm-3m") + print([parameter.name for parameter in space_group_parameters.latpars]) + +which prints: + +.. code-block:: text + + ['a'] + +For a cubic space group only ``a`` is left free. ``b`` and ``c`` are +constrained to it, and the angles are fixed to 90 degrees. + +To calculate a PDF from the adapter, pass it to the generator's +``set_structure_from_parset`` method instead of passing a structure object +to ``set_structure``: + +.. code-block:: python + + from diffpy.cmipdf.pdfgenerator import PDFGenerator + + generator = PDFGenerator("nickel_pdf") + generator.set_structure_from_parset(phase) + +Add the adapter to ``structure_to_parameter_set`` +------------------------------------------------- + +``diffpy.cmistructure.structure_to_parameter_set``, which is also what a +generator's ``set_structure`` method uses, only knows about the adapters +listed at the top of this page. For any other type it raises +``TypeError("Unadaptable structure format")``. Outside this package, create +your adapter directly as shown above. + +If you are adding an adapter to ``diffpy.cmistructure`` itself: + +#. Put the adapter in a new module, ``src/diffpy/cmistructure/parset.py``, + with the package header and NumPy-style docstrings used by the other + adapter modules. +#. Add a ``can_adapt`` branch to ``structure_to_parameter_set`` in + ``src/diffpy/cmistructure/__init__.py``. Import the adapter inside the + function, as the existing branches do, so that the optional package it + depends on is imported only when it is needed. +#. Add tests in ``tests/test_parset.py``. At a minimum, cover + ``can_adapt`` for supported and unsupported objects, that changing each + ``Parameter`` changes the wrapped structure, and that + ``constrain_as_space_group`` produces the expected free parameters. +#. Add the module to ``docs/source/api/diffpy.cmistructure.rst`` and a news + item to ``news/``. + +Checklist +--------- + +- The adapter subclasses ``SrRealParSet``, or ``BaseStructureParSet`` if it + is never used with srreal calculators. +- ``structure`` holds the wrapped object, and every ``ParameterAdapter`` + reads and writes that object, not a copy. +- ``can_adapt`` returns False, without raising, for unsupported objects. +- ``get_lattice()`` returns parameters named ``a``, ``b``, ``c``, + ``alpha``, ``beta`` and ``gamma``, and has ``angle_units`` set. +- ``get_scatterers()`` returns one ``ParameterSet`` per atom, with + fractional ``x``, ``y`` and ``z`` and uniquely named. +- The ADPs use the default ``U11``…``U23`` and ``Uiso`` names, or you pass + ``adpsymbols`` and ``isosymbol`` to ``constrain_as_space_group``. +- ``_get_srreal_structure()`` returns an object srreal can read and + respects ``using_symmetry()``. diff --git a/docs/source/migrating-from-srfit.rst b/docs/source/migrating-from-srfit.rst new file mode 100644 index 0000000..ba17c60 --- /dev/null +++ b/docs/source/migrating-from-srfit.rst @@ -0,0 +1,203 @@ +:tocdepth: -1 + +.. index:: migrating-from-srfit + +.. _migrating-from-srfit: + +========================================= +Migrating from ``diffpy.srfit.structure`` +========================================= + +``diffpy.cmistructure`` holds the code that used to live in +``diffpy.srfit.structure``. The behavior is the same, but the names follow +the Python style: camelCase methods are now snake_case, and abbreviated +attribute and argument names are spelled out. + +``diffpy.srfit`` 3.5.0 warns when you use an old method or attribute and names +its replacement. **Renamed function arguments can't be warned about**, so check +the :ref:`renamed-arguments` table if your script passes arguments by keyword. + +Update your imports +------------------- + +Replace ``diffpy.srfit.structure`` with ``diffpy.cmistructure``. The module +names are unchanged: + +.. code-block:: python + + # Before + from diffpy.srfit.structure import constrainAsSpaceGroup + from diffpy.srfit.structure.objcrystparset import ObjCrystCrystalParSet + + # After + from diffpy.cmistructure import constrain_as_space_group + from diffpy.cmistructure.objcrystparset import ObjCrystCrystalParSet + +``constrain_as_space_group`` accepts a generator +------------------------------------------------ + +``constrain_as_space_group`` now also accepts a generator, such as a +``PDFGenerator``, in place of its ``.phase``. Both of these calls work: + +.. code-block:: python + + constrain_as_space_group(generator.phase, "Fm-3m") + constrain_as_space_group(generator, "Fm-3m") + +Anything that is neither a structure parameter set nor an object with one as +its ``.phase`` raises a ``TypeError``. + +What has not changed +-------------------- + +- Class names, such as ``DiffpyStructureParSet`` and + ``ObjCrystCrystalParSet``. +- The space group containers ``.sgpars``, ``.latpars``, ``.xyzpars`` and + ``.adppars``. They also name Parameters in fit results, for example + ``sgpars.latpars.a``. +- The ``constrainadps`` and ``periodic`` keyword arguments. + +Renamed functions and methods +----------------------------- + +.. list-table:: + :header-rows: 1 + + * - diffpy.srfit.structure + - diffpy.cmistructure + * - ``struToParameterSet(name, stru)`` + - ``structure_to_parameter_set(name, structure)`` + * - ``constrainAsSpaceGroup`` + - ``constrain_as_space_group`` + * - ``canAdapt`` (all adapters) + - ``can_adapt`` + * - ``getLattice`` (all adapters) + - ``get_lattice`` + * - ``getScatterers`` (all adapters) + - ``get_scatterers`` + * - ``CCTBXCrystalParSet.getSpaceGroup`` + - ``get_space_group`` + * - ``SrRealParSet.restrainBVS`` + - ``restrain_bvs`` + * - ``useSymmetry`` + - ``use_symmetry`` + * - ``usingSymmetry`` + - ``using_symmetry`` + * - ``isDummy`` + - ``is_dummy`` + * - ``ObjCrystScattererParSet.hasScatterers`` + - ``has_scatterers`` + * - ``ObjCrystMoleculeParSet.wrapRestraints`` + - ``wrap_restraints`` + * - ``ObjCrystMoleculeParSet.wrapStretchModeParameters`` + - ``wrap_stretch_mode_parameters`` + * - ``ObjCrystMoleculeParSet.restrainBondLength`` + - ``restrain_bond_length`` + * - ``ObjCrystMoleculeParSet.restrainBondLengthParameter`` + - ``restrain_bond_length_parameter`` + * - ``ObjCrystMoleculeParSet.restrainBondAngle`` + - ``restrain_bond_angle`` + * - ``ObjCrystMoleculeParSet.restrainBondAngleParameter`` + - ``restrain_bond_angle_parameter`` + * - ``ObjCrystMoleculeParSet.restrainDihedralAngle`` + - ``restrain_dihedral_angle`` + * - ``ObjCrystMoleculeParSet.restrainDihedralAngleParameter`` + - ``restrain_dihedral_angle_parameter`` + * - ``ObjCrystMoleculeParSet.addBondLengthParameter`` + - ``add_bond_length_parameter`` + * - ``ObjCrystMoleculeParSet.addBondAngleParameter`` + - ``add_bond_angle_parameter`` + * - ``ObjCrystMoleculeParSet.addDihedralAngleParameter`` + - ``add_dihedral_angle_parameter`` + * - ``StretchModeParameter.addAtoms`` + - ``add_atoms`` + * - ``setConst`` on the bond length, bond angle and dihedral angle + Parameters + - ``set_constant`` + * - ``BaseSpaceGroupParameters.addParameter`` + - ``add_parameter`` + +Renamed attributes +------------------ + +.. list-table:: + :header-rows: 1 + + * - Class + - diffpy.srfit.structure + - diffpy.cmistructure + * - ``SrRealParSet``, ``DiffpyStructureParSet``, + ``ObjCrystCrystalParSet``, ``ObjCrystMoleculeParSet``, + ``CCTBXCrystalParSet`` + - ``.stru`` + - ``.structure`` + * - ``DiffpyLatticeParSet``, ``ObjCrystCrystalParSet`` + - ``.angunits`` + - ``.angle_units`` + * - ``ObjCrystScattererParSet`` + - ``.scat`` + - ``.scatterer`` + * - ``ObjCrystMoleculeRestraint`` + - ``.res`` + - ``.restraint`` + * - Bond length, bond angle and dihedral angle Parameters + - ``.matoms`` + - ``.mutated_atoms`` + * - ``SpaceGroupParameters`` + - ``.sg`` + - ``.space_group`` + * - ``CCTBXScattererParSet``, ``CCTBXUnitCellParSet`` + - ``.strups`` + - ``.structure_parameter_set`` + * - ``CCTBXScattererParSet`` + - ``.idx`` + - ``.index`` + +.. _renamed-arguments: + +Renamed arguments +----------------- + +These only matter if you pass the argument by keyword. Positional calls keep +working, because no argument has moved. + +.. list-table:: + :header-rows: 1 + + * - Function or class + - diffpy.srfit.structure + - diffpy.cmistructure + * - ``structure_to_parameter_set``, ``can_adapt``, + ``DiffpyStructureParSet``, ``CCTBXCrystalParSet`` + - ``stru`` + - ``structure`` + * - ``ObjCrystCrystalParSet`` + - ``cryst`` + - ``crystal`` + * - ``ObjCrystScattererParSet``, ``ObjCrystMolAtomParSet`` + - ``scat`` + - ``scatterer`` + * - ``ObjCrystMoleculeRestraint`` + - ``res`` + - ``restraint`` + * - ``BVSRestraint`` + - ``parset`` + - ``parameter_set`` + * - ``CCTBXScattererParSet``, ``CCTBXUnitCellParSet`` + - ``strups`` + - ``structure_parameter_set`` + * - ``CCTBXScattererParSet`` + - ``idx`` + - ``index`` + * - ``SpaceGroupParameters`` + - ``sg`` + - ``space_group`` + * - ``restrain_bond_length_parameter``, + ``restrain_bond_angle_parameter``, + ``restrain_dihedral_angle_parameter``, + ``BaseSpaceGroupParameters.add_parameter`` + - ``par`` + - ``parameter`` + * - ``StretchModeParameter.set_value`` + - ``val`` + - ``value`` diff --git a/docs/source/snippets/.placeholder b/docs/source/snippets/.placeholder deleted file mode 100644 index e69de29..0000000 diff --git a/docs/source/snippets/example-table.rst b/docs/source/snippets/example-table.rst deleted file mode 100644 index 7c4c11d..0000000 --- a/docs/source/snippets/example-table.rst +++ /dev/null @@ -1,28 +0,0 @@ -.. list-table:: 5 levels of reusing/sharing code - :widths: 5 15 40 40 - :header-rows: 1 - - * - Level - - Name - - Scope - - How to setup - * - 1 - - ``function`` - - Reuse code in the single file. - - See Level 1 tutorial - * - 2 - - ``module`` - - Reuse code across files. - - See Level 2 tutorial - * - 3 - - ``workspace`` - - Reuse code across project folders. - - ``package create workspace`` - * - 4 - - ``system`` - - Reuse code across any files in the computer. - - ``package create system`` - * - 5 - - ``public`` - - Share code as publicly installable package. - - ``package create public`` diff --git a/news/structure-tests-and-docs.rst b/news/structure-tests-and-docs.rst new file mode 100644 index 0000000..abd206e --- /dev/null +++ b/news/structure-tests-and-docs.rst @@ -0,0 +1,25 @@ +**Added:** + +* A "Migrating from diffpy.srfit.structure" documentation page that lists every renamed function, method, attribute and argument. +* The ``npintensity.py`` and ``npintensityII.py`` examples from diffpy.srfit, updated to the diffpy.cmistructure names. +* A "Making a structure adapter" documentation page that explains how to adapt another structure representation for use with diffpy.cmistructure. + +**Changed:** + +* + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/tests/conftest.py b/tests/conftest.py index 584fbbb..f322641 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -4,11 +4,79 @@ import sys from pathlib import Path +import numpy import pytest +from pyobjcryst.atom import Atom +from pyobjcryst.crystal import Crystal +from pyobjcryst.molecule import Molecule +from pyobjcryst.scatteringpower import ScatteringPowerAtom from diffpy.cmistructure.diffpyparset import DiffpyStructureParSet from diffpy.structure import Structure +C60_XYZ = """\ +3.451266498 0.685000000 0.000000000 +3.451266498 -0.685000000 0.000000000 +-3.451266498 0.685000000 0.000000000 +-3.451266498 -0.685000000 0.000000000 +0.685000000 0.000000000 3.451266498 +-0.685000000 0.000000000 3.451266498 +0.685000000 0.000000000 -3.451266498 +-0.685000000 0.000000000 -3.451266498 +0.000000000 3.451266498 0.685000000 +0.000000000 3.451266498 -0.685000000 +0.000000000 -3.451266498 0.685000000 +0.000000000 -3.451266498 -0.685000000 +3.003809890 1.409000000 1.171456608 +3.003809890 1.409000000 -1.171456608 +3.003809890 -1.409000000 1.171456608 +3.003809890 -1.409000000 -1.171456608 +-3.003809890 1.409000000 1.171456608 +-3.003809890 1.409000000 -1.171456608 +-3.003809890 -1.409000000 1.171456608 +-3.003809890 -1.409000000 -1.171456608 +1.409000000 1.171456608 3.003809890 +1.409000000 -1.171456608 3.003809890 +-1.409000000 1.171456608 3.003809890 +-1.409000000 -1.171456608 3.003809890 +1.409000000 1.171456608 -3.003809890 +1.409000000 -1.171456608 -3.003809890 +-1.409000000 1.171456608 -3.003809890 +-1.409000000 -1.171456608 -3.003809890 +1.171456608 3.003809890 1.409000000 +-1.171456608 3.003809890 1.409000000 +1.171456608 3.003809890 -1.409000000 +-1.171456608 3.003809890 -1.409000000 +1.171456608 -3.003809890 1.409000000 +-1.171456608 -3.003809890 1.409000000 +1.171456608 -3.003809890 -1.409000000 +-1.171456608 -3.003809890 -1.409000000 +2.580456608 0.724000000 2.279809890 +2.580456608 0.724000000 -2.279809890 +2.580456608 -0.724000000 2.279809890 +2.580456608 -0.724000000 -2.279809890 +-2.580456608 0.724000000 2.279809890 +-2.580456608 0.724000000 -2.279809890 +-2.580456608 -0.724000000 2.279809890 +-2.580456608 -0.724000000 -2.279809890 +0.724000000 2.279809890 2.580456608 +0.724000000 -2.279809890 2.580456608 +-0.724000000 2.279809890 2.580456608 +-0.724000000 -2.279809890 2.580456608 +0.724000000 2.279809890 -2.580456608 +0.724000000 -2.279809890 -2.580456608 +-0.724000000 2.279809890 -2.580456608 +-0.724000000 -2.279809890 -2.580456608 +2.279809890 2.580456608 0.724000000 +-2.279809890 2.580456608 0.724000000 +2.279809890 2.580456608 -0.724000000 +-2.279809890 2.580456608 -0.724000000 +2.279809890 -2.580456608 0.724000000 +-2.279809890 -2.580456608 0.724000000 +2.279809890 -2.580456608 -0.724000000 +-2.279809890 -2.580456608 -0.724000000 +""" + @pytest.fixture def user_filesystem(tmp_path): @@ -61,9 +129,73 @@ def __init__(self, phase): @pytest.fixture -def lamno3_generator(datafile): +def lamno3_p1_structure(datafile): + """Return the full LaMnO3 unit cell in P1 as a diffpy Structure.""" structure = Structure() structure.read(datafile("LaMnO3.stru")) - parameter_set = DiffpyStructureParSet("LaMnO3", structure) + return structure + + +@pytest.fixture +def lamno3_generator(lamno3_p1_structure): + parameter_set = DiffpyStructureParSet("LaMnO3", lamno3_p1_structure) generator = GeneratorStandIn(parameter_set) return generator + + +@pytest.fixture +def lamno3_crystal(): + """Return the LaMnO3 asymmetric unit in P b n m as a pyobjcryst + Crystal.""" + pi = numpy.pi + # It appears that ObjCryst only supports standard symbols + crystal = Crystal(5.486341, 5.619215, 7.628206, "P b n m") + crystal.SetName("LaMnO3") + # La1 + sp = ScatteringPowerAtom("La1", "La") + sp.SetBiso(8 * pi * pi * 0.003) + atom = Atom(0.996096, 0.0321494, 0.25, "La1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # Mn1 + sp = ScatteringPowerAtom("Mn1", "Mn") + sp.SetBiso(8 * pi * pi * 0.003) + atom = Atom(0, 0.5, 0, "Mn1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # O1 + sp = ScatteringPowerAtom("O1", "O") + sp.SetBiso(8 * pi * pi * 0.003) + atom = Atom(0.0595746, 0.496164, 0.25, "O1", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + # O2 + sp = ScatteringPowerAtom("O2", "O") + sp.SetBiso(8 * pi * pi * 0.003) + atom = Atom(0.720052, 0.289387, 0.0311126, "O2", sp) + crystal.AddScatteringPower(sp) + crystal.AddScatterer(atom) + return crystal + + +@pytest.fixture +def c60_crystal(): + """Return a P1 pyobjcryst Crystal holding the C60 molecule "c60".""" + pi = numpy.pi + c = Crystal(100, 100, 100, "P1") + c.SetName("c60frame") + m = Molecule(c, "c60") + c.AddScatterer(m) + sp = ScatteringPowerAtom("C", "C") + sp.SetBiso(8 * pi * pi * 0.003) + for i, l in enumerate(C60_XYZ.strip().splitlines()): + x, y, z = map(float, l.split()) + m.AddAtom(x, y, z, sp, "C%i" % i) + return c + + +@pytest.fixture +def c60_molecule(c60_crystal): + """Return the C60 Molecule held by the c60_crystal fixture.""" + molecule = c60_crystal.GetScatterer("c60") + return molecule diff --git a/tests/test_diffpyparset.py b/tests/test_diffpyparset.py index a8ef41f..b0f6a97 100644 --- a/tests/test_diffpyparset.py +++ b/tests/test_diffpyparset.py @@ -23,12 +23,11 @@ import numpy as np from diffpy.cmistructure.diffpyparset import DiffpyStructureParSet +from diffpy.structure import Atom, Lattice, Structure -def testDiffpyStructureParSet(): +def test_diffpy_structure_parset(): """Test the structure conversion.""" - from diffpy.structure import Atom, Lattice, Structure - a1 = Atom("Cu", xyz=np.array([0.0, 0.1, 0.2]), Uisoequiv=0.003) a2 = Atom("Ag", xyz=np.array([0.3, 0.4, 0.5]), Uisoequiv=0.002) lattice = Lattice(2.5, 2.5, 2.5, 90, 90, 90) @@ -138,8 +137,6 @@ def _testLattice(): def test___repr__(): """Test representation of DiffpyStructureParSet objects.""" - from diffpy.structure import Atom, Lattice, Structure - lat = Lattice(3, 3, 2, 90, 90, 90) atom = Atom("C", [0, 0.2, 0.5]) structure = Structure([atom], lattice=lat) @@ -154,8 +151,6 @@ def test___repr__(): def test_pickling(): """Test pickling of DiffpyStructureParSet.""" - from diffpy.structure import Atom, Structure - structure = Structure([Atom("C", [0, 0.2, 0.5])]) dsps = DiffpyStructureParSet("dsps", structure) data = pickle.dumps(dsps) @@ -176,8 +171,6 @@ def test_show_bvs_restraint(capturestdout): # nacl. # Expected: The Restraints section lists it as the BVS of nacl with its # sig and scaling. - from diffpy.structure import Atom, Lattice, Structure - lattice = Lattice(5.6, 5.6, 5.6, 90, 90, 90) atoms = [Atom("Na", [0, 0, 0]), Atom("Cl", [0.5, 0.5, 0.5])] parameter_set = DiffpyStructureParSet( @@ -188,6 +181,3 @@ def test_show_bvs_restraint(capturestdout): actual_restraint_lines = lines[lines.index("Restraints") + 2 :] expected_restraint_lines = ["BVS of nacl: sig = 0.500000, scaled = True"] assert actual_restraint_lines == expected_restraint_lines - - -# End of class TestParameterAdapter diff --git a/tests/test_objcrystparset.py b/tests/test_objcrystparset.py index 7226371..cfbcffe 100644 --- a/tests/test_objcrystparset.py +++ b/tests/test_objcrystparset.py @@ -20,836 +20,699 @@ import numpy import pytest +from pyobjcryst.spacegroup import SpaceGroup -# Global variables to be assigned in setUp -ObjCrystCrystalParSet = spacegroups = None -Crystal = Atom = Molecule = ScatteringPowerAtom = None - - -c60xyz = """\ -3.451266498 0.685000000 0.000000000 -3.451266498 -0.685000000 0.000000000 --3.451266498 0.685000000 0.000000000 --3.451266498 -0.685000000 0.000000000 -0.685000000 0.000000000 3.451266498 --0.685000000 0.000000000 3.451266498 -0.685000000 0.000000000 -3.451266498 --0.685000000 0.000000000 -3.451266498 -0.000000000 3.451266498 0.685000000 -0.000000000 3.451266498 -0.685000000 -0.000000000 -3.451266498 0.685000000 -0.000000000 -3.451266498 -0.685000000 -3.003809890 1.409000000 1.171456608 -3.003809890 1.409000000 -1.171456608 -3.003809890 -1.409000000 1.171456608 -3.003809890 -1.409000000 -1.171456608 --3.003809890 1.409000000 1.171456608 --3.003809890 1.409000000 -1.171456608 --3.003809890 -1.409000000 1.171456608 --3.003809890 -1.409000000 -1.171456608 -1.409000000 1.171456608 3.003809890 -1.409000000 -1.171456608 3.003809890 --1.409000000 1.171456608 3.003809890 --1.409000000 -1.171456608 3.003809890 -1.409000000 1.171456608 -3.003809890 -1.409000000 -1.171456608 -3.003809890 --1.409000000 1.171456608 -3.003809890 --1.409000000 -1.171456608 -3.003809890 -1.171456608 3.003809890 1.409000000 --1.171456608 3.003809890 1.409000000 -1.171456608 3.003809890 -1.409000000 --1.171456608 3.003809890 -1.409000000 -1.171456608 -3.003809890 1.409000000 --1.171456608 -3.003809890 1.409000000 -1.171456608 -3.003809890 -1.409000000 --1.171456608 -3.003809890 -1.409000000 -2.580456608 0.724000000 2.279809890 -2.580456608 0.724000000 -2.279809890 -2.580456608 -0.724000000 2.279809890 -2.580456608 -0.724000000 -2.279809890 --2.580456608 0.724000000 2.279809890 --2.580456608 0.724000000 -2.279809890 --2.580456608 -0.724000000 2.279809890 --2.580456608 -0.724000000 -2.279809890 -0.724000000 2.279809890 2.580456608 -0.724000000 -2.279809890 2.580456608 --0.724000000 2.279809890 2.580456608 --0.724000000 -2.279809890 2.580456608 -0.724000000 2.279809890 -2.580456608 -0.724000000 -2.279809890 -2.580456608 --0.724000000 2.279809890 -2.580456608 --0.724000000 -2.279809890 -2.580456608 -2.279809890 2.580456608 0.724000000 --2.279809890 2.580456608 0.724000000 -2.279809890 2.580456608 -0.724000000 --2.279809890 2.580456608 -0.724000000 -2.279809890 -2.580456608 0.724000000 --2.279809890 -2.580456608 0.724000000 -2.279809890 -2.580456608 -0.724000000 --2.279809890 -2.580456608 -0.724000000 -""" - - -def makeC60(): - """Make a crystal containing the C60 molecule using pyobjcryst.""" - pi = numpy.pi - c = Crystal(100, 100, 100, "P1") - c.SetName("c60frame") - m = Molecule(c, "c60") - - c.AddScatterer(m) - - sp = ScatteringPowerAtom("C", "C") - sp.SetBiso(8 * pi * pi * 0.003) - # c.AddScatteringPower(sp) - - for i, l in enumerate(c60xyz.strip().splitlines()): - x, y, z = map(float, l.split()) - m.AddAtom(x, y, z, sp, "C%i" % i) - - return c - +from diffpy.cmistructure.objcrystparset import ObjCrystCrystalParSet +from diffpy.structure import spacegroups # ---------------------------------------------------------------------------- -class TestParameterAdapter: - @pytest.fixture(autouse=True) - def setup(self): - # shared setup - global ObjCrystCrystalParSet, Crystal, Atom, Molecule - global ScatteringPowerAtom - from pyobjcryst.atom import Atom - from pyobjcryst.crystal import Crystal - from pyobjcryst.molecule import Molecule - from pyobjcryst.scatteringpower import ScatteringPowerAtom - - from diffpy.cmistructure.objcrystparset import ObjCrystCrystalParSet - - self.occryst = makeC60() - self.ocmol = self.occryst.GetScatterer("c60") - return - - def tearDown(self): - del self.occryst - del self.ocmol - return - - def testImplicitBondAngleRestraints(self): - """Test the structure with implicit bond angles.""" - occryst = self.occryst - ocmol = self.ocmol - - # Add some bond angles to the molecule - ocmol.AddBondAngle(ocmol[0], ocmol[5], ocmol[8], 1.1, 0.1, 0.1) - ocmol.AddBondAngle(ocmol[0], ocmol[7], ocmol[44], 1.3, 0.1, 0.1) - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - m.wrap_restraints() - - # C1: Two restraints are added to the molecule. - # Expected: The molecule holds both restraints. - actual_restraint_count = len(m._restraints) - expected_restraint_count = 2 - assert actual_restraint_count == expected_restraint_count - - # C2: The restraint penalties are evaluated. - # Expected: They equal the pyobjcryst log-likelihoods. - res0, res1 = m._restraints - actual_penalties = set([res0.penalty(), res1.penalty()]) - angles = ocmol.GetBondAngleList() - expected_penalties = set( - [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] - ) - assert actual_penalties == expected_penalties - +def test_implicit_bond_angle_restraints(c60_crystal, c60_molecule): + """Test the structure with implicit bond angles.""" + occryst = c60_crystal + ocmol = c60_molecule + + # Add some bond angles to the molecule + ocmol.AddBondAngle(ocmol[0], ocmol[5], ocmol[8], 1.1, 0.1, 0.1) + ocmol.AddBondAngle(ocmol[0], ocmol[7], ocmol[44], 1.3, 0.1, 0.1) + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + m.wrap_restraints() + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + res0, res1 = m._restraints + actual_penalties = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetBondAngleList() + expected_penalties = set( + [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + return + + +def test_objcryst_parset(c60_crystal, c60_molecule): + """Test the structure conversion.""" + occryst = c60_crystal + ocmol = c60_molecule + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + actual_name = crystal.name + expected_name = "bucky" + assert actual_name == expected_name + + def _testCrystal(): + # Test the lattice + actual_lattice = [ + crystal.a.value, + crystal.b.get_value(), + crystal.c.get_value(), + crystal.alpha.get_value(), + crystal.beta.get_value(), + crystal.gamma.get_value(), + ] + expected_lattice = [ + occryst.a, + occryst.b, + occryst.c, + occryst.alpha, + occryst.beta, + occryst.gamma, + ] + assert actual_lattice == pytest.approx(expected_lattice) return - def testObjCrystParSet(self): - """Test the structure conversion.""" - occryst = self.occryst - ocmol = self.ocmol - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - - actual_name = crystal.name - expected_name = "bucky" - assert actual_name == expected_name - - def _testCrystal(): - # Test the lattice - actual_lattice = [ - crystal.a.value, - crystal.b.get_value(), - crystal.c.get_value(), - crystal.alpha.get_value(), - crystal.beta.get_value(), - crystal.gamma.get_value(), - ] - expected_lattice = [ - occryst.a, - occryst.b, - occryst.c, - occryst.alpha, - occryst.beta, - occryst.gamma, + def _testMolecule(): + # Test position, occupancy and orientation + actual_molecule = [ + m.x.get_value(), + m.y.get_value(), + m.z.get_value(), + m.occ.get_value(), + m.q0.get_value(), + m.q1.get_value(), + m.q2.get_value(), + m.q3.get_value(), + ] + expected_molecule = [ + ocmol.X, + ocmol.Y, + ocmol.Z, + ocmol.Occupancy, + ocmol.Q0, + ocmol.Q1, + ocmol.Q2, + ocmol.Q3, + ] + assert actual_molecule == pytest.approx(expected_molecule) + + # Check the atoms thoroughly + actual_elements = [a.element for a in m.atoms] + expected_elements = [ + ocmol[i].GetScatteringPower().GetSymbol() + for i in range(len(ocmol)) + ] + assert actual_elements == expected_elements + actual_atoms = [ + [ + a.x.get_value(), + a.y.get_value(), + a.z.get_value(), + a.occ.get_value(), + a.Biso.get_value(), ] - assert actual_lattice == pytest.approx(expected_lattice) - return - - def _testMolecule(): - # Test position, occupancy and orientation - actual_molecule = [ - m.x.get_value(), - m.y.get_value(), - m.z.get_value(), - m.occ.get_value(), - m.q0.get_value(), - m.q1.get_value(), - m.q2.get_value(), - m.q3.get_value(), - ] - expected_molecule = [ - ocmol.X, - ocmol.Y, - ocmol.Z, - ocmol.Occupancy, - ocmol.Q0, - ocmol.Q1, - ocmol.Q2, - ocmol.Q3, - ] - assert actual_molecule == pytest.approx(expected_molecule) - - # Check the atoms thoroughly - actual_elements = [a.element for a in m.atoms] - expected_elements = [ - ocmol[i].GetScatteringPower().GetSymbol() - for i in range(len(ocmol)) - ] - assert actual_elements == expected_elements - actual_atoms = [ + for a in m.atoms + ] + expected_atoms = [ + pytest.approx( [ - a.x.get_value(), - a.y.get_value(), - a.z.get_value(), - a.occ.get_value(), - a.Biso.get_value(), + ocmol[i].X, + ocmol[i].Y, + ocmol[i].Z, + ocmol[i].Occupancy, + ocmol[i].GetScatteringPower().Biso, ] - for a in m.atoms - ] - expected_atoms = [ - pytest.approx( - [ - ocmol[i].X, - ocmol[i].Y, - ocmol[i].Z, - ocmol[i].Occupancy, - ocmol[i].GetScatteringPower().Biso, - ] - ) - for i in range(len(ocmol)) - ] - assert actual_atoms == expected_atoms - return - - # C1: The ParameterSet has just been created from the crystal. - # Expected: The Parameters match the pyobjcryst values. - _testCrystal() - _testMolecule() - - # C2: Values are changed through pyobjcryst. - # Expected: The Parameters follow the changes. - ocmol[0].X *= 1.1 - ocmol[0].Occupancy *= 1.1 - ocmol[0].GetScatteringPower().Biso *= 1.1 - ocmol.Q0 *= 1.1 - occryst.a *= 1.1 - - _testCrystal() - _testMolecule() - - # C3: Values are changed through the ParameterSet. - # Expected: The pyobjcryst objects follow the changes. - crystal.c60.C44.x.set_value(1.1) - crystal.c60.C44.occ.set_value(1.1) - crystal.c60.C44.Biso.set_value(1.1) - crystal.c60.q3.set_value(1.1) - crystal.a.set_value(1.1) - - _testCrystal() - _testMolecule() - return - - def testImplicitBondLengthRestraints(self): - """Test the structure with implicit bond lengths.""" - occryst = self.occryst - ocmol = self.ocmol - - # Add some bonds to the molecule - ocmol.AddBond(ocmol[0], ocmol[5], 3.3, 0.1, 0.1) - ocmol.AddBond(ocmol[0], ocmol[7], 3.3, 0.1, 0.1) - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - m.wrap_restraints() - - # C1: Two restraints are added to the molecule. - # Expected: The molecule holds both restraints. - actual_restraint_count = len(m._restraints) - expected_restraint_count = 2 - assert actual_restraint_count == expected_restraint_count - - # C2: The restraint penalties are evaluated. - # Expected: They equal the pyobjcryst log-likelihoods. - res0, res1 = m._restraints - actual_penalties = set([res0.penalty(), res1.penalty()]) - bonds = ocmol.GetBondList() - expected_penalties = set( - [bonds[0].GetLogLikelihood(), bonds[1].GetLogLikelihood()] - ) - assert actual_penalties == expected_penalties - + ) + for i in range(len(ocmol)) + ] + assert actual_atoms == expected_atoms return - def testImplicitDihedralAngleRestraints(self): - """Test the structure with implicit dihedral angles.""" - occryst = self.occryst - ocmol = self.ocmol - - # Add some bond angles to the molecule - ocmol.AddDihedralAngle( - ocmol[0], ocmol[5], ocmol[8], ocmol[41], 1.1, 0.1, 0.1 - ) - ocmol.AddDihedralAngle( - ocmol[0], ocmol[7], ocmol[44], ocmol[2], 1.3, 0.1, 0.1 - ) - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - m.wrap_restraints() - - # C1: Two restraints are added to the molecule. - # Expected: The molecule holds both restraints. - actual_restraint_count = len(m._restraints) - expected_restraint_count = 2 - assert actual_restraint_count == expected_restraint_count - - # C2: The restraint penalties are evaluated. - # Expected: They equal the pyobjcryst log-likelihoods. - res0, res1 = m._restraints - actual_penalties = set([res0.penalty(), res1.penalty()]) - angles = ocmol.GetDihedralAngleList() - expected_penalties = set( - [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] - ) - assert actual_penalties == expected_penalties - - return - - def testImplicitStretchModes(self): - """Test the molecule with implicit stretch modes.""" - # Not sure how to make this happen. - pass - - def testExplicitBondLengthRestraints(self): - """Test the structure with explicit bond lengths.""" - occryst = self.occryst - ocmol = self.ocmol - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - - # make some bond angle restraints - res0 = m.restrain_bond_length(m.atoms[0], m.atoms[5], 3.3, 0.1, 0.1) - res1 = m.restrain_bond_length(m.atoms[0], m.atoms[7], 3.3, 0.1, 0.1) - - # C1: Two restraints are added to the molecule. - # Expected: The molecule holds both restraints. - actual_restraint_count = len(m._restraints) - expected_restraint_count = 2 - assert actual_restraint_count == expected_restraint_count - - # C2: The restraint penalties are evaluated. - # Expected: They equal the pyobjcryst log-likelihoods. - bonds = ocmol.GetBondList() - actual_bond_count = len(bonds) - expected_bond_count = 2 - assert actual_bond_count == expected_bond_count - actual_penalties = [res0.penalty(), res1.penalty()] - expected_penalties = [b.GetLogLikelihood() for b in bonds] - assert actual_penalties == expected_penalties - - return - - def test_show_explicit_molecule_restraints(self, capturestdout): - # C1: A bond length, a bond angle and a dihedral angle restraint are - # added through the ObjCrystMoleculeParSet. - # Expected: The Restraints section lists each restraint by its atom - # names with its target value, sigma, delta and scaling. - crystal = ObjCrystCrystalParSet("bucky", self.occryst) - molecule = crystal.c60 - atoms = molecule.atoms - molecule.restrain_bond_length(atoms[0], atoms[5], 3.3, 0.1, 0.2) - molecule.restrain_bond_angle( - atoms[0], atoms[5], atoms[7], 1.0, 0.1, 0.2 - ) - molecule.restrain_dihedral_angle( - atoms[0], atoms[5], atoms[7], atoms[8], 2.0, 0.1, 0.2 - ) - lines = capturestdout(crystal.show, textwidth=0).splitlines() - actual_restraint_lines = set(lines[lines.index("Restraints") + 2 :]) - expected_restraint_lines = { - "C0-C5: length = 3.300000, sigma = 0.100000, delta = 0.200000, " - "scaled = False", - "C0-C5-C7: angle = 1.000000, sigma = 0.100000, delta = 0.200000, " - "scaled = False", - "C0-C5-C7-C8: angle = 2.000000, sigma = 0.100000, " - "delta = 0.200000, scaled = False", - } - assert actual_restraint_lines == expected_restraint_lines - - def test_show_implicit_molecule_restraints(self, capturestdout): - # C1: A bond length, a bond angle and a dihedral angle restraint are - # defined on the pyobjcryst Molecule and wrapped with - # wrap_restraints. - # Expected: The Restraints section lists each restraint by its atom - # names with its target value, sigma, delta and scaling. - ocmol = self.ocmol - ocmol.AddBond(ocmol[0], ocmol[5], 3.3, 0.1, 0.2) - ocmol.AddBondAngle(ocmol[0], ocmol[5], ocmol[7], 1.0, 0.1, 0.2) - ocmol.AddDihedralAngle( - ocmol[0], ocmol[5], ocmol[7], ocmol[8], 2.0, 0.1, 0.2 - ) - crystal = ObjCrystCrystalParSet("bucky", self.occryst) - crystal.c60.wrap_restraints() - lines = capturestdout(crystal.show, textwidth=0).splitlines() - actual_restraint_lines = set(lines[lines.index("Restraints") + 2 :]) - expected_restraint_lines = { - "C0-C5: length = 3.300000, sigma = 0.100000, delta = 0.200000, " - "scaled = False", - "C0-C5-C7: angle = 1.000000, sigma = 0.100000, delta = 0.200000, " - "scaled = False", - "C0-C5-C7-C8: angle = 2.000000, sigma = 0.100000, " - "delta = 0.200000, scaled = False", - } - assert actual_restraint_lines == expected_restraint_lines - - def testExplicitBondAngleRestraints(self): - """Test the structure with explicit bond angles. - - Note that this cannot work with co-linear points as the - direction of rotation cannot be defined in this case. - """ - occryst = self.occryst - ocmol = self.ocmol - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - - # restrain some bond angles - res0 = m.restrain_bond_angle( - m.atoms[0], m.atoms[5], m.atoms[8], 3.3, 0.1, 0.1 - ) - res1 = m.restrain_bond_angle( - m.atoms[0], m.atoms[7], m.atoms[44], 3.3, 0.1, 0.1 - ) - - # C1: Two restraints are added to the molecule. - # Expected: The molecule holds both restraints. - actual_restraint_count = len(m._restraints) - expected_restraint_count = 2 - assert actual_restraint_count == expected_restraint_count - - # C2: The restraint penalties are evaluated. - # Expected: They equal the pyobjcryst log-likelihoods. - actual_penalties = set([res0.penalty(), res1.penalty()]) - angles = ocmol.GetBondAngleList() - expected_penalties = set( - [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] - ) - assert actual_penalties == expected_penalties - - return - - def testExplicitDihedralAngleRestraints(self): - """Test the structure with explicit dihedral angles.""" - occryst = self.occryst - ocmol = self.ocmol - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - - # Restrain some dihedral angles. - res0 = m.restrain_dihedral_angle( - m.atoms[0], m.atoms[5], m.atoms[8], m.atoms[41], 1.1, 0.1, 0.1 - ) - res1 = m.restrain_dihedral_angle( - m.atoms[0], m.atoms[7], m.atoms[44], m.atoms[2], 1.1, 0.1, 0.1 - ) - - # C1: Two restraints are added to the molecule. - # Expected: The molecule holds both restraints. - actual_restraint_count = len(m._restraints) - expected_restraint_count = 2 - assert actual_restraint_count == expected_restraint_count - - # C2: The restraint penalties are evaluated. - # Expected: They equal the pyobjcryst log-likelihoods. - actual_penalties = set([res0.penalty(), res1.penalty()]) - angles = ocmol.GetDihedralAngleList() - expected_penalties = set( - [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] - ) - assert actual_penalties == expected_penalties - - return - - def testExplicitBondLengthParameter(self): - """Test adding bond length parameters to the molecule.""" - occryst = self.occryst - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - - a0 = m.atoms[0] - a7 = m.atoms[7] - a20 = m.atoms[20] - - # Add a parameter - p1 = m.add_bond_length_parameter("C07", a0, a7) - # Have another atom tag along for the ride - p1.add_atoms([a20]) - - xyz0 = numpy.array( - [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] - ) - xyz7 = numpy.array( - [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] - ) - xyz20 = numpy.array( - [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] - ) - - dd = xyz0 - xyz7 - d0 = numpy.dot(dd, dd) ** 0.5 - # C1: A bond length Parameter is added. - # Expected: Its value is the current bond length. - actual_length = p1.get_value() - expected_length = d0 - assert actual_length == pytest.approx(expected_length, abs=1e-6) - - # Record the unit direction of change for later - u = dd / d0 - - # C2: The bond length Parameter is stretched by 5%. - # Expected: The Parameter and the measured bond length both take - # the new value, the first atom stays put, and the second and - # tag-along atoms move along the bond. - scale = 1.05 - p1.set_value(scale * d0) - - actual_length = p1.get_value() - expected_length = scale * d0 - assert actual_length == pytest.approx(expected_length, abs=1e-6) - - xyz0a = numpy.array( - [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] - ) - xyz7a = numpy.array( - [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] - ) - xyz20a = numpy.array( - [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] - ) - - dda = xyz0a - xyz7a - d1 = numpy.dot(dda, dda) ** 0.5 - - actual_measured_length = d1 - expected_measured_length = scale * d0 - assert actual_measured_length == pytest.approx( - expected_measured_length, abs=1e-6 - ) - - actual_xyz0 = xyz0a.tolist() - expected_xyz0 = xyz0.tolist() - assert actual_xyz0 == expected_xyz0 - - actual_xyz7 = xyz7a - expected_xyz7 = xyz7 + (1 - scale) * d0 * u - assert actual_xyz7 == pytest.approx(expected_xyz7, abs=1e-5) - - actual_xyz20 = xyz20a - expected_xyz20 = xyz20 + (1 - scale) * d0 * u - assert actual_xyz20 == pytest.approx(expected_xyz20, abs=1e-6) - - return - - def testExplicitBondAngleParameter(self): - """Test adding bond angle parameters to the molecule.""" - occryst = self.occryst - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - - a0 = m.atoms[0] - a7 = m.atoms[7] - a20 = m.atoms[20] - a25 = m.atoms[25] - - xyz0 = numpy.array( - [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] - ) - xyz7 = numpy.array( - [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] - ) - xyz20 = numpy.array( - [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] - ) - xyz25 = numpy.array( - [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] - ) - - v1 = xyz7 - xyz0 - d1 = numpy.dot(v1, v1) ** 0.5 - v2 = xyz7 - xyz20 - d2 = numpy.dot(v2, v2) ** 0.5 - - angle0 = numpy.arccos(numpy.dot(v1, v2) / (d1 * d2)) - - # Add a parameter - p1 = m.add_bond_angle_parameter("C0720", a0, a7, a20) - # Have another atom tag along for the ride - p1.add_atoms([a25]) - - # C1: A bond angle Parameter is added. - # Expected: Its value is the current bond angle. - actual_angle = p1.get_value() - expected_angle = angle0 - assert actual_angle == pytest.approx(expected_angle, abs=1e-6) - - # C2: The bond angle Parameter is stretched by 5%. - # Expected: The Parameter and the measured angle both take the new - # value, and only the third and tag-along atoms move. - scale = 1.05 - p1.set_value(scale * angle0) - - actual_angle = p1.get_value() - expected_angle = scale * angle0 - assert actual_angle == pytest.approx(expected_angle, abs=1e-6) - - xyz0a = numpy.array( - [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] - ) - xyz7a = numpy.array( - [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] - ) - xyz20a = numpy.array( - [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] - ) - xyz25a = numpy.array( - [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] - ) - - v1a = xyz7a - xyz0a - d1a = numpy.dot(v1a, v1a) ** 0.5 - v2a = xyz7a - xyz20a - d2a = numpy.dot(v2a, v2a) ** 0.5 - - angle1 = numpy.arccos(numpy.dot(v1a, v2a) / (d1a * d2a)) - - actual_measured_angle = angle1 - expected_measured_angle = scale * angle0 - assert actual_measured_angle == pytest.approx( - expected_measured_angle, abs=1e-6 - ) - - actual_moved = { - "C0": not numpy.array_equal(xyz0, xyz0a), - "C7": not numpy.array_equal(xyz7, xyz7a), - "C20": not numpy.array_equal(xyz20, xyz20a), - "C25": not numpy.array_equal(xyz25, xyz25a), - } - expected_moved = {"C0": False, "C7": False, "C20": True, "C25": True} - assert actual_moved == expected_moved - - return - - def testExplicitDihedralAngleParameter(self): - """Test adding dihedral angle parameters to the molecule.""" - occryst = self.occryst - - # make our crystal - crystal = ObjCrystCrystalParSet("bucky", occryst) - m = crystal.c60 - - a0 = m.atoms[0] - a7 = m.atoms[7] - a20 = m.atoms[20] - a25 = m.atoms[25] - a33 = m.atoms[33] - - xyz0 = numpy.array( - [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] - ) - xyz7 = numpy.array( - [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] - ) - xyz20 = numpy.array( - [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] - ) - xyz25 = numpy.array( - [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] - ) - xyz33 = numpy.array( - [a33.x.get_value(), a33.y.get_value(), a33.z.get_value()] - ) - - v12 = xyz0 - xyz7 - v23 = xyz7 - xyz20 - v34 = xyz20 - xyz25 - v123 = numpy.cross(v12, v23) - v234 = numpy.cross(v23, v34) - - d123 = numpy.dot(v123, v123) ** 0.5 - d234 = numpy.dot(v234, v234) ** 0.5 - angle0 = -numpy.arccos(numpy.dot(v123, v234) / (d123 * d234)) - - # Add a parameter - p1 = m.add_dihedral_angle_parameter("C072025", a0, a7, a20, a25) - # Have another atom tag along for the ride - p1.add_atoms([a33]) - - # C1: A dihedral angle Parameter is added. - # Expected: Its value is the current dihedral angle. - actual_angle = p1.get_value() - expected_angle = angle0 - assert actual_angle == pytest.approx(expected_angle, abs=1e-6) - - # C2: The dihedral angle Parameter is stretched by 5%. - # Expected: The Parameter and the measured angle both take the new - # value, and only the fourth and tag-along atoms move. - scale = 1.05 - p1.set_value(scale * angle0) - - actual_angle = p1.get_value() - expected_angle = scale * angle0 - assert actual_angle == pytest.approx(expected_angle, abs=1e-6) - - xyz0a = numpy.array( - [a0.x.get_value(), a0.y.get_value(), a0.z.get_value()] - ) - xyz7a = numpy.array( - [a7.x.get_value(), a7.y.get_value(), a7.z.get_value()] - ) - xyz20a = numpy.array( - [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] - ) - xyz25a = numpy.array( - [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] - ) - xyz33a = numpy.array( - [a33.x.get_value(), a33.y.get_value(), a33.z.get_value()] - ) - - v12a = xyz0a - xyz7a - v23a = xyz7a - xyz20a - v34a = xyz20a - xyz25a - v123a = numpy.cross(v12a, v23a) - v234a = numpy.cross(v23a, v34a) - - d123a = numpy.dot(v123a, v123a) ** 0.5 - d234a = numpy.dot(v234a, v234a) ** 0.5 - angle1 = -numpy.arccos(numpy.dot(v123a, v234a) / (d123a * d234a)) - actual_measured_angle = angle1 - expected_measured_angle = scale * angle0 - assert actual_measured_angle == pytest.approx( - expected_measured_angle, abs=1e-6 - ) - - actual_moved = { - "C0": not numpy.array_equal(xyz0, xyz0a), - "C7": not numpy.array_equal(xyz7, xyz7a), - "C20": not numpy.array_equal(xyz20, xyz20a), - "C25": not numpy.array_equal(xyz25, xyz25a), - "C33": not numpy.array_equal(xyz33, xyz33a), - } - expected_moved = { - "C0": False, - "C7": False, - "C20": False, - "C25": True, - "C33": True, - } - assert actual_moved == expected_moved - - return - - -class TestCreateSpaceGroup: - """Test space group creation from pyobjcryst structures. - - This makes sure that the space groups created by the structure - parameter set are correct. + # C1: The ParameterSet has just been created from the crystal. + # Expected: The Parameters match the pyobjcryst values. + _testCrystal() + _testMolecule() + + # C2: Values are changed through pyobjcryst. + # Expected: The Parameters follow the changes. + ocmol[0].X *= 1.1 + ocmol[0].Occupancy *= 1.1 + ocmol[0].GetScatteringPower().Biso *= 1.1 + ocmol.Q0 *= 1.1 + occryst.a *= 1.1 + + _testCrystal() + _testMolecule() + + # C3: Values are changed through the ParameterSet. + # Expected: The pyobjcryst objects follow the changes. + crystal.c60.C44.x.set_value(1.1) + crystal.c60.C44.occ.set_value(1.1) + crystal.c60.C44.Biso.set_value(1.1) + crystal.c60.q3.set_value(1.1) + crystal.a.set_value(1.1) + + _testCrystal() + _testMolecule() + return + + +def test_implicit_bond_length_restraints(c60_crystal, c60_molecule): + """Test the structure with implicit bond lengths.""" + occryst = c60_crystal + ocmol = c60_molecule + + # Add some bonds to the molecule + ocmol.AddBond(ocmol[0], ocmol[5], 3.3, 0.1, 0.1) + ocmol.AddBond(ocmol[0], ocmol[7], 3.3, 0.1, 0.1) + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + m.wrap_restraints() + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + res0, res1 = m._restraints + actual_penalties = set([res0.penalty(), res1.penalty()]) + bonds = ocmol.GetBondList() + expected_penalties = set( + [bonds[0].GetLogLikelihood(), bonds[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + return + + +def test_implicit_dihedral_angle_restraints(c60_crystal, c60_molecule): + """Test the structure with implicit dihedral angles.""" + occryst = c60_crystal + ocmol = c60_molecule + + # Add some bond angles to the molecule + ocmol.AddDihedralAngle( + ocmol[0], ocmol[5], ocmol[8], ocmol[41], 1.1, 0.1, 0.1 + ) + ocmol.AddDihedralAngle( + ocmol[0], ocmol[7], ocmol[44], ocmol[2], 1.3, 0.1, 0.1 + ) + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + m.wrap_restraints() + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + res0, res1 = m._restraints + actual_penalties = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetDihedralAngleList() + expected_penalties = set( + [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + return + + +def test_implicit_stretch_modes(): + """Test the molecule with implicit stretch modes.""" + # Not sure how to make this happen. + pass + + +def test_explicit_bond_length_restraints(c60_crystal, c60_molecule): + """Test the structure with explicit bond lengths.""" + occryst = c60_crystal + ocmol = c60_molecule + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + # make some bond angle restraints + res0 = m.restrain_bond_length(m.atoms[0], m.atoms[5], 3.3, 0.1, 0.1) + res1 = m.restrain_bond_length(m.atoms[0], m.atoms[7], 3.3, 0.1, 0.1) + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + bonds = ocmol.GetBondList() + actual_bond_count = len(bonds) + expected_bond_count = 2 + assert actual_bond_count == expected_bond_count + actual_penalties = [res0.penalty(), res1.penalty()] + expected_penalties = [b.GetLogLikelihood() for b in bonds] + assert actual_penalties == expected_penalties + + return + + +def test_show_explicit_molecule_restraints(c60_crystal, capturestdout): + # C1: A bond length, a bond angle and a dihedral angle restraint are + # added through the ObjCrystMoleculeParSet. + # Expected: The Restraints section lists each restraint by its atom + # names with its target value, sigma, delta and scaling. + crystal = ObjCrystCrystalParSet("bucky", c60_crystal) + molecule = crystal.c60 + atoms = molecule.atoms + molecule.restrain_bond_length(atoms[0], atoms[5], 3.3, 0.1, 0.2) + molecule.restrain_bond_angle(atoms[0], atoms[5], atoms[7], 1.0, 0.1, 0.2) + molecule.restrain_dihedral_angle( + atoms[0], atoms[5], atoms[7], atoms[8], 2.0, 0.1, 0.2 + ) + lines = capturestdout(crystal.show, textwidth=0).splitlines() + actual_restraint_lines = set(lines[lines.index("Restraints") + 2 :]) + expected_restraint_lines = { + "C0-C5: length = 3.300000, sigma = 0.100000, delta = 0.200000, " + "scaled = False", + "C0-C5-C7: angle = 1.000000, sigma = 0.100000, delta = 0.200000, " + "scaled = False", + "C0-C5-C7-C8: angle = 2.000000, sigma = 0.100000, " + "delta = 0.200000, scaled = False", + } + assert actual_restraint_lines == expected_restraint_lines + + +def test_show_implicit_molecule_restraints( + c60_crystal, c60_molecule, capturestdout +): + # C1: A bond length, a bond angle and a dihedral angle restraint are + # defined on the pyobjcryst Molecule and wrapped with + # wrap_restraints. + # Expected: The Restraints section lists each restraint by its atom + # names with its target value, sigma, delta and scaling. + ocmol = c60_molecule + ocmol.AddBond(ocmol[0], ocmol[5], 3.3, 0.1, 0.2) + ocmol.AddBondAngle(ocmol[0], ocmol[5], ocmol[7], 1.0, 0.1, 0.2) + ocmol.AddDihedralAngle( + ocmol[0], ocmol[5], ocmol[7], ocmol[8], 2.0, 0.1, 0.2 + ) + crystal = ObjCrystCrystalParSet("bucky", c60_crystal) + crystal.c60.wrap_restraints() + lines = capturestdout(crystal.show, textwidth=0).splitlines() + actual_restraint_lines = set(lines[lines.index("Restraints") + 2 :]) + expected_restraint_lines = { + "C0-C5: length = 3.300000, sigma = 0.100000, delta = 0.200000, " + "scaled = False", + "C0-C5-C7: angle = 1.000000, sigma = 0.100000, delta = 0.200000, " + "scaled = False", + "C0-C5-C7-C8: angle = 2.000000, sigma = 0.100000, " + "delta = 0.200000, scaled = False", + } + assert actual_restraint_lines == expected_restraint_lines + + +def test_explicit_bond_angle_restraints(c60_crystal, c60_molecule): + """Test the structure with explicit bond angles. + + Note that this cannot work with co-linear points as the direction of + rotation cannot be defined in this case. """ - - @pytest.fixture(autouse=True) - def setup(self): - # shared setup - global ObjCrystCrystalParSet, spacegroups - from diffpy.cmistructure.objcrystparset import ObjCrystCrystalParSet - from diffpy.structure import spacegroups - - @staticmethod - def getObjCrystParSetSpaceGroup(space_group): - """Make an ObjCrystCrystalParSet with the proper space group.""" - from pyobjcryst.spacegroup import SpaceGroup - - sgobjcryst = SpaceGroup(space_group.short_name) - sgnew = ObjCrystCrystalParSet._create_space_group(sgobjcryst) - return sgnew - - @staticmethod - def hashDiffPySpaceGroup(space_group): - lines = [str(space_group.number % 1000)] + sorted( - map(str, space_group.iter_symops()) - ) - s = "\n".join(lines) - return s - - def sgsEquivalent(self, sg1, sg2): - """Check to see if two space group objects are the same.""" - hash1 = self.hashDiffPySpaceGroup(sg1) - hash2 = self.hashDiffPySpaceGroup(sg2) - return hash1 == hash2 - - # FIXME: only about 50% of the spacegroups pass the assertion - # test disabled even if cctbx is installed - def xtestCreateSpaceGroup(self): - """Check all sgtbx space groups for proper conversion to - SpaceGroup.""" - from cctbx import sgtbx - - for smbls in sgtbx.space_group_symbol_iterator(): - shn = smbls.hermann_mauguin() - short_name = shn.replace(" ", "") - if spacegroups.is_space_group_identifier(short_name): - space_group = spacegroups.get_space_group(shn) - sgnew = self.getObjCrystParSetSpaceGroup(space_group) - actual_equivalent = self.sgsEquivalent(space_group, sgnew) - expected_equivalent = True - assert actual_equivalent == expected_equivalent - return - - -# End of class TestCreateSpaceGroup + occryst = c60_crystal + ocmol = c60_molecule + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + # restrain some bond angles + res0 = m.restrain_bond_angle( + m.atoms[0], m.atoms[5], m.atoms[8], 3.3, 0.1, 0.1 + ) + res1 = m.restrain_bond_angle( + m.atoms[0], m.atoms[7], m.atoms[44], 3.3, 0.1, 0.1 + ) + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + actual_penalties = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetBondAngleList() + expected_penalties = set( + [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + return + + +def test_explicit_dihedral_angle_restraints(c60_crystal, c60_molecule): + """Test the structure with explicit dihedral angles.""" + occryst = c60_crystal + ocmol = c60_molecule + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + # Restrain some dihedral angles. + res0 = m.restrain_dihedral_angle( + m.atoms[0], m.atoms[5], m.atoms[8], m.atoms[41], 1.1, 0.1, 0.1 + ) + res1 = m.restrain_dihedral_angle( + m.atoms[0], m.atoms[7], m.atoms[44], m.atoms[2], 1.1, 0.1, 0.1 + ) + + # C1: Two restraints are added to the molecule. + # Expected: The molecule holds both restraints. + actual_restraint_count = len(m._restraints) + expected_restraint_count = 2 + assert actual_restraint_count == expected_restraint_count + + # C2: The restraint penalties are evaluated. + # Expected: They equal the pyobjcryst log-likelihoods. + actual_penalties = set([res0.penalty(), res1.penalty()]) + angles = ocmol.GetDihedralAngleList() + expected_penalties = set( + [angles[0].GetLogLikelihood(), angles[1].GetLogLikelihood()] + ) + assert actual_penalties == expected_penalties + + return + + +def test_explicit_bond_length_parameter(c60_crystal): + """Test adding bond length parameters to the molecule.""" + occryst = c60_crystal + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + a0 = m.atoms[0] + a7 = m.atoms[7] + a20 = m.atoms[20] + + # Add a parameter + p1 = m.add_bond_length_parameter("C07", a0, a7) + # Have another atom tag along for the ride + p1.add_atoms([a20]) + + xyz0 = numpy.array([a0.x.get_value(), a0.y.get_value(), a0.z.get_value()]) + xyz7 = numpy.array([a7.x.get_value(), a7.y.get_value(), a7.z.get_value()]) + xyz20 = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + + dd = xyz0 - xyz7 + d0 = numpy.dot(dd, dd) ** 0.5 + # C1: A bond length Parameter is added. + # Expected: Its value is the current bond length. + actual_length = p1.get_value() + expected_length = d0 + assert actual_length == pytest.approx(expected_length, abs=1e-6) + + # Record the unit direction of change for later + u = dd / d0 + + # C2: The bond length Parameter is stretched by 5%. + # Expected: The Parameter and the measured bond length both take + # the new value, the first atom stays put, and the second and + # tag-along atoms move along the bond. + scale = 1.05 + p1.set_value(scale * d0) + + actual_length = p1.get_value() + expected_length = scale * d0 + assert actual_length == pytest.approx(expected_length, abs=1e-6) + + xyz0a = numpy.array([a0.x.get_value(), a0.y.get_value(), a0.z.get_value()]) + xyz7a = numpy.array([a7.x.get_value(), a7.y.get_value(), a7.z.get_value()]) + xyz20a = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + + dda = xyz0a - xyz7a + d1 = numpy.dot(dda, dda) ** 0.5 + + actual_measured_length = d1 + expected_measured_length = scale * d0 + assert actual_measured_length == pytest.approx( + expected_measured_length, abs=1e-6 + ) + + actual_xyz0 = xyz0a.tolist() + expected_xyz0 = xyz0.tolist() + assert actual_xyz0 == expected_xyz0 + + actual_xyz7 = xyz7a + expected_xyz7 = xyz7 + (1 - scale) * d0 * u + assert actual_xyz7 == pytest.approx(expected_xyz7, abs=1e-5) + + actual_xyz20 = xyz20a + expected_xyz20 = xyz20 + (1 - scale) * d0 * u + assert actual_xyz20 == pytest.approx(expected_xyz20, abs=1e-6) + + return + + +def test_explicit_bond_angle_parameter(c60_crystal): + """Test adding bond angle parameters to the molecule.""" + occryst = c60_crystal + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + a0 = m.atoms[0] + a7 = m.atoms[7] + a20 = m.atoms[20] + a25 = m.atoms[25] + + xyz0 = numpy.array([a0.x.get_value(), a0.y.get_value(), a0.z.get_value()]) + xyz7 = numpy.array([a7.x.get_value(), a7.y.get_value(), a7.z.get_value()]) + xyz20 = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + xyz25 = numpy.array( + [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] + ) + + v1 = xyz7 - xyz0 + d1 = numpy.dot(v1, v1) ** 0.5 + v2 = xyz7 - xyz20 + d2 = numpy.dot(v2, v2) ** 0.5 + + angle0 = numpy.arccos(numpy.dot(v1, v2) / (d1 * d2)) + + # Add a parameter + p1 = m.add_bond_angle_parameter("C0720", a0, a7, a20) + # Have another atom tag along for the ride + p1.add_atoms([a25]) + + # C1: A bond angle Parameter is added. + # Expected: Its value is the current bond angle. + actual_angle = p1.get_value() + expected_angle = angle0 + assert actual_angle == pytest.approx(expected_angle, abs=1e-6) + + # C2: The bond angle Parameter is stretched by 5%. + # Expected: The Parameter and the measured angle both take the new + # value, and only the third and tag-along atoms move. + scale = 1.05 + p1.set_value(scale * angle0) + + actual_angle = p1.get_value() + expected_angle = scale * angle0 + assert actual_angle == pytest.approx(expected_angle, abs=1e-6) + + xyz0a = numpy.array([a0.x.get_value(), a0.y.get_value(), a0.z.get_value()]) + xyz7a = numpy.array([a7.x.get_value(), a7.y.get_value(), a7.z.get_value()]) + xyz20a = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + xyz25a = numpy.array( + [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] + ) + + v1a = xyz7a - xyz0a + d1a = numpy.dot(v1a, v1a) ** 0.5 + v2a = xyz7a - xyz20a + d2a = numpy.dot(v2a, v2a) ** 0.5 + + angle1 = numpy.arccos(numpy.dot(v1a, v2a) / (d1a * d2a)) + + actual_measured_angle = angle1 + expected_measured_angle = scale * angle0 + assert actual_measured_angle == pytest.approx( + expected_measured_angle, abs=1e-6 + ) + + actual_moved = { + "C0": not numpy.array_equal(xyz0, xyz0a), + "C7": not numpy.array_equal(xyz7, xyz7a), + "C20": not numpy.array_equal(xyz20, xyz20a), + "C25": not numpy.array_equal(xyz25, xyz25a), + } + expected_moved = {"C0": False, "C7": False, "C20": True, "C25": True} + assert actual_moved == expected_moved + + return + + +def test_explicit_dihedral_angle_parameter(c60_crystal): + """Test adding dihedral angle parameters to the molecule.""" + occryst = c60_crystal + + # make our crystal + crystal = ObjCrystCrystalParSet("bucky", occryst) + m = crystal.c60 + + a0 = m.atoms[0] + a7 = m.atoms[7] + a20 = m.atoms[20] + a25 = m.atoms[25] + a33 = m.atoms[33] + + xyz0 = numpy.array([a0.x.get_value(), a0.y.get_value(), a0.z.get_value()]) + xyz7 = numpy.array([a7.x.get_value(), a7.y.get_value(), a7.z.get_value()]) + xyz20 = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + xyz25 = numpy.array( + [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] + ) + xyz33 = numpy.array( + [a33.x.get_value(), a33.y.get_value(), a33.z.get_value()] + ) + + v12 = xyz0 - xyz7 + v23 = xyz7 - xyz20 + v34 = xyz20 - xyz25 + v123 = numpy.cross(v12, v23) + v234 = numpy.cross(v23, v34) + + d123 = numpy.dot(v123, v123) ** 0.5 + d234 = numpy.dot(v234, v234) ** 0.5 + angle0 = -numpy.arccos(numpy.dot(v123, v234) / (d123 * d234)) + + # Add a parameter + p1 = m.add_dihedral_angle_parameter("C072025", a0, a7, a20, a25) + # Have another atom tag along for the ride + p1.add_atoms([a33]) + + # C1: A dihedral angle Parameter is added. + # Expected: Its value is the current dihedral angle. + actual_angle = p1.get_value() + expected_angle = angle0 + assert actual_angle == pytest.approx(expected_angle, abs=1e-6) + + # C2: The dihedral angle Parameter is stretched by 5%. + # Expected: The Parameter and the measured angle both take the new + # value, and only the fourth and tag-along atoms move. + scale = 1.05 + p1.set_value(scale * angle0) + + actual_angle = p1.get_value() + expected_angle = scale * angle0 + assert actual_angle == pytest.approx(expected_angle, abs=1e-6) + + xyz0a = numpy.array([a0.x.get_value(), a0.y.get_value(), a0.z.get_value()]) + xyz7a = numpy.array([a7.x.get_value(), a7.y.get_value(), a7.z.get_value()]) + xyz20a = numpy.array( + [a20.x.get_value(), a20.y.get_value(), a20.z.get_value()] + ) + xyz25a = numpy.array( + [a25.x.get_value(), a25.y.get_value(), a25.z.get_value()] + ) + xyz33a = numpy.array( + [a33.x.get_value(), a33.y.get_value(), a33.z.get_value()] + ) + + v12a = xyz0a - xyz7a + v23a = xyz7a - xyz20a + v34a = xyz20a - xyz25a + v123a = numpy.cross(v12a, v23a) + v234a = numpy.cross(v23a, v34a) + + d123a = numpy.dot(v123a, v123a) ** 0.5 + d234a = numpy.dot(v234a, v234a) ** 0.5 + angle1 = -numpy.arccos(numpy.dot(v123a, v234a) / (d123a * d234a)) + actual_measured_angle = angle1 + expected_measured_angle = scale * angle0 + assert actual_measured_angle == pytest.approx( + expected_measured_angle, abs=1e-6 + ) + + actual_moved = { + "C0": not numpy.array_equal(xyz0, xyz0a), + "C7": not numpy.array_equal(xyz7, xyz7a), + "C20": not numpy.array_equal(xyz20, xyz20a), + "C25": not numpy.array_equal(xyz25, xyz25a), + "C33": not numpy.array_equal(xyz33, xyz33a), + } + expected_moved = { + "C0": False, + "C7": False, + "C20": False, + "C25": True, + "C33": True, + } + assert actual_moved == expected_moved + + return + + +def get_objcryst_parset_space_group(space_group): + """Make an ObjCrystCrystalParSet with the proper space group.""" + sgobjcryst = SpaceGroup(space_group.short_name) + sgnew = ObjCrystCrystalParSet._create_space_group(sgobjcryst) + return sgnew + + +def get_diffpy_space_group_hash(space_group): + lines = [str(space_group.number % 1000)] + sorted( + map(str, space_group.iter_symops()) + ) + s = "\n".join(lines) + return s + + +def are_space_groups_equivalent(sg1, sg2): + """Check to see if two space group objects are the same.""" + hash1 = get_diffpy_space_group_hash(sg1) + hash2 = get_diffpy_space_group_hash(sg2) + return hash1 == hash2 + + +# FIXME: only about 50% of the spacegroups pass the assertion +# test disabled even if cctbx is installed +def xtest_create_space_group(): + """Check all sgtbx space groups for proper conversion to + SpaceGroup.""" + from cctbx import sgtbx + + for smbls in sgtbx.space_group_symbol_iterator(): + shn = smbls.hermann_mauguin() + short_name = shn.replace(" ", "") + if spacegroups.is_space_group_identifier(short_name): + space_group = spacegroups.get_space_group(shn) + sgnew = get_objcryst_parset_space_group(space_group) + actual_equivalent = are_space_groups_equivalent(space_group, sgnew) + expected_equivalent = True + assert actual_equivalent == expected_equivalent + return diff --git a/tests/test_sgconstraints.py b/tests/test_sgconstraints.py index 3707d3e..9ce4296 100644 --- a/tests/test_sgconstraints.py +++ b/tests/test_sgconstraints.py @@ -22,9 +22,6 @@ import numpy import pytest -from pyobjcryst.atom import Atom -from pyobjcryst.crystal import Crystal -from pyobjcryst.scatteringpower import ScatteringPowerAtom from diffpy.cmistructure.diffpyparset import DiffpyStructureParSet from diffpy.cmistructure.objcrystparset import ObjCrystCrystalParSet @@ -36,7 +33,7 @@ # ---------------------------------------------------------------------------- -def test_ObjCryst_constrain_space_group(): +def test_objcryst_constrain_space_group(lamno3_crystal): """Make sure that all Parameters are constrained properly. This tests constrainSpaceGroup from @@ -45,7 +42,7 @@ def test_ObjCryst_constrain_space_group(): """ pi = numpy.pi - occryst = makeLaMnO3() + occryst = lamno3_crystal structure = ObjCrystCrystalParSet(occryst.GetName(), occryst) # Make sure we actually create the constraints structure._constrain_space_group() @@ -124,9 +121,9 @@ def test_ObjCryst_constrain_space_group(): return -def test_DiffPy_constrain_as_space_group(datafile): +def test_diffpy_constrain_as_space_group(lamno3_p1_structure): """Test the constrain_as_space_group function.""" - structure = makeLaMnO3_P1(datafile) + structure = lamno3_p1_structure parameter_set = DiffpyStructureParSet("LaMnO3", structure) space_group_parameters = constrain_as_space_group( @@ -226,17 +223,17 @@ def _alltests(parameter): return -def test_constrain_as_space_group_args(datafile): +def test_constrain_as_space_group_args(lamno3_p1_structure): """Test the arguments processing of constrain_as_space_group function.""" # C1: The space group is given as a symbol or as a SpaceGroup object. # Expected: Both create the same space group Parameters. - structure = makeLaMnO3_P1(datafile) + structure = lamno3_p1_structure parameter_set = DiffpyStructureParSet("LaMnO3", structure) symbol_parameters = constrain_as_space_group(parameter_set, "P b n m") space_group = get_space_group("P b n m") object_parameter_set = DiffpyStructureParSet( - "LMO", makeLaMnO3_P1(datafile) + "LMO", lamno3_p1_structure.copy() ) object_parameters = constrain_as_space_group( object_parameter_set, space_group @@ -283,7 +280,6 @@ def test_constrain_as_space_group_generator(lamno3_generator): ], ) def test_constrain_as_space_group_bad_phase(input_phase, expected_type_name): - expected_msg = ( "Expected a structure parameter set or a generator with a .phase " f"attribute, got {expected_type_name}. " @@ -308,67 +304,24 @@ def test_constrain_as_space_group_bad_phase(input_phase, expected_type_name): ], ) def test_space_group_parameter_container_names( - datafile, input_container_attribute, expected_container_name + lamno3_p1_structure, input_container_attribute, expected_container_name ): - - parameter_set = DiffpyStructureParSet("LaMnO3", makeLaMnO3_P1(datafile)) + parameter_set = DiffpyStructureParSet("LaMnO3", lamno3_p1_structure) space_group_parameters = constrain_as_space_group(parameter_set, "P b n m") container = getattr(space_group_parameters, input_container_attribute) actual_container_name = container.name assert actual_container_name == expected_container_name -def test_objcryst_space_group_parameters_name(): +def test_objcryst_space_group_parameters_name(lamno3_crystal): # C1: The space group Parameters are read from an ObjCrystCrystalParSet. # Expected: They are available as .sgpars, in a container named # "sgpars". - - occryst = makeLaMnO3() + occryst = lamno3_crystal structure = ObjCrystCrystalParSet(occryst.GetName(), occryst) actual_container_name = structure.sgpars.name expected_container_name = "sgpars" assert actual_container_name == expected_container_name -def makeLaMnO3_P1(datafile): - - structure = Structure() - structure.read(datafile("LaMnO3.stru")) - return structure - - -def makeLaMnO3(): - - pi = numpy.pi - # It appears that ObjCryst only supports standard symbols - crystal = Crystal(5.486341, 5.619215, 7.628206, "P b n m") - crystal.SetName("LaMnO3") - # La1 - sp = ScatteringPowerAtom("La1", "La") - sp.SetBiso(8 * pi * pi * 0.003) - atom = Atom(0.996096, 0.0321494, 0.25, "La1", sp) - crystal.AddScatteringPower(sp) - crystal.AddScatterer(atom) - # Mn1 - sp = ScatteringPowerAtom("Mn1", "Mn") - sp.SetBiso(8 * pi * pi * 0.003) - atom = Atom(0, 0.5, 0, "Mn1", sp) - crystal.AddScatteringPower(sp) - crystal.AddScatterer(atom) - # O1 - sp = ScatteringPowerAtom("O1", "O") - sp.SetBiso(8 * pi * pi * 0.003) - atom = Atom(0.0595746, 0.496164, 0.25, "O1", sp) - crystal.AddScatteringPower(sp) - crystal.AddScatterer(atom) - # O2 - sp = ScatteringPowerAtom("O2", "O") - sp.SetBiso(8 * pi * pi * 0.003) - atom = Atom(0.720052, 0.289387, 0.0311126, "O2", sp) - crystal.AddScatteringPower(sp) - crystal.AddScatterer(atom) - - return crystal - - # ----------------------------------------------------------------------------