Random generation of molecular crystals

[qzhu2017]

@scottfredericks , I suggest you start with the code which we had for atomic crystals.

1, copy structure.py to a new file (structure_molecule.py?)
2, for a given molecule (like H2O), we still randomly pick a space group and wyckoff site, and then align the molecules with some random numbers.
3, In the structure class, we now have
numIons, space group, lattice, coordinates, wyckoffs
it would be good to add some parameters which represents the necessary information to reproduce the structure, like the followings,
molecular xyz
the alignment of molecules and random numbers.

I will come up with more suggestions later

[scottfredericks]

@qzhu2017
Hopefully most of it should be simple. However, I'm not sure how the lattice's volume factor should change based on the molecule being used, or how we should check inter-molecular distances. Perhaps we could use some kind of minimum sphere or ellipsoid fitting to determine an approximate volume for the molecule? Using a spherical model would be easy to check, and would not require any change in the algorithm. However, calculating a sphere for the molecule may not be straightforward, since the center of mass of the molecule is not necessarily the geometric center of the binding sphere.

For an ellipsoid or other shape (perhaps something like an atomic orbital expansion?) I am not sure if an efficient algorithm exists for checking overlap (especially if we want to consider periodic boundary conditions). And again, calculating the binding surface could be difficult. But, it could prove slightly more accurate, and hopefully faster than checking each individual atomic distance.

Of course, it would be simplest to just check inter-atomic distances directly. However, I am not sure if the minimum distances need to be altered, since molecular bonds might be stronger than crystallographic ones. I think the computation cost goes up by a factor of n^2.

Please let me know what you think is best. In the mean time, I will work on getting the rest of the code finished.

[qzhu2017]

@scottfredericks
For the estimation of volume, you can use the following scheme,
1, reorient the molecule based on the eigenvalue and eigenvector of the molecular inertia tensor;
2, calculate the volume by treating it as an orthogonal box with a certain thickness say 2 \AA, so the volume is (x_max-x_min+2)(y_max-y_min+2)(z_max-z_min+2)

For the minimum distance, we only check the intermolecular distance (i.e., distance between two atoms belonging to two molecules). We should set the minimum distance larger than the convenient bond length.

For the easy case, I suggest you use the example of H2O and CH4. Please upload your code asap. I have some matlab codes dealing with the volume estimation and distance check and I can easily modify your code to make it work.

[scottfredericks]

@qzhu2017
I have implemented the volume and box-finding methods. The relevant code is in molecular_crystal.py. However, looking at check_compatibility, I realized that we will need to call orientation_in_wyckoff_position many times per crystal, which would slow down the code too much.

To avoid this problem, I can build a class for storing orientations with 0, 1, or 2 contraint vectors. orientation_in_wyckoff_position can return these objects, so that we only need to call it once per Wyckoff position/molecule (while checking the compatibility). Then, once a list of orientations is obtained for each Wyckoff position/molecule pair, we can easily choose a random orientation without recalling the function.

I think I can code this today. Once done, I will put the orientation class in matrix.py.

[scottfredericks]

@qzhu2017
I think most of the functionality is programmed. However, I still need to edit merge_coordinate to check for valid orientations, and to make sure that the final coordinates are transformed correctly. Right now, there is an error in the process, so the molecules tend to overlap. I will continue work tomorrow.

[scottfredericks]

@qzhu2017
I am now obtaining reasonable-looking results. In order to obtain results in a reasonable timeframe, I increased the default volume factor to 100. By default, molecular_crystal.py will attempt to generate water ice 1h: spacegroup P6_3/mmc (194), with 12 H2O molecules per unit cell. After a few attempts, some of the generated crystals have a hexagonal motif similar to that of actual ice. I am not sure if the structure is actually correct though, so please take a look.

I think that once we have a better distance checking method, we can reduce the volume factor without slowing down so much.

The main difficulties were:

1. Storing and accessing valid molecular orientations for each Wyckoff position/molecule
2. checking that merged wyckoff positions also have valid molecular orientations
3. generating the correct relative coordinates for the molecules once the centers are known
4. finding and using the correct generating (x,y,z) point for merged Wyckoff positions, so that the corresponding rotations could be generated

[qzhu2017]

@scottfredericks Your code doesn't work at the moment. I fixed a couple of issue, but still cannot get the code run. It returns the following errors:

/Users/qiangzhu/miniconda3/lib/python3.6/site-packages/pymatgen/io/cif.py:44: UserWarning: Please install optional dependency pybtex if youwant to extract references from CIF files.
  warnings.warn("Please install optional dependency pybtex if you"
Calculating molecular orientations... Done.
Traceback (most recent call last):
  File "molecular_crystal.py", line 573, in <module>
    rand_crystal = molecular_crystal(options.sg, system, numMols0, options.factor, orientations=orientations)
  File "molecular_crystal.py", line 309, in __init__
    self.generate_crystal()
  File "molecular_crystal.py", line 496, in generate_crystal
    op1 = choose(self.valid_orientations[i][j][k]).get_op()
AttributeError: 'SymmOp' object has no attribute 'get_op'

[scottfredericks]

@qzhu2017
I have updated the code - I had run into the same error, but forgot to upload the needed change in molecule.py. I believe everything is now updated, and the program works very fast on my computer - changing the stoichiometry has sped things up dramatically. Let me know if you run into any other issues.

[qzhu2017]

@scottfredericks The code runs very fast. It is a bit strange that we need to use such big volume factor. I am going to check the distance calculation function soon. Could you please check another issue in the mean time? If you check the symmetry of the generated structures, sometimes, it ends up with the structure with lower symmetry (P2_1/m), instead of the desired one (Cmc2_1). Maybe, there exist some bugs in orientation.

qiangzhu@u-10-240-0-36:~/Desktop/github/crystallography$ grep space_group_name out/*.cif
out/1.cif:_symmetry_space_group_name_H-M   P2_1/m
out/10.cif:_symmetry_space_group_name_H-M   Cmc2_1
out/2.cif:_symmetry_space_group_name_H-M   Cmc2_1
out/3.cif:_symmetry_space_group_name_H-M   Cmc2_1
out/4.cif:_symmetry_space_group_name_H-M   P2_1/m
out/5.cif:_symmetry_space_group_name_H-M   P2_1/m
out/6.cif:_symmetry_space_group_name_H-M   P2_1/m
out/7.cif:_symmetry_space_group_name_H-M   P2_1/m
out/8.cif:_symmetry_space_group_name_H-M   Cmc2_1
out/9.cif:_symmetry_space_group_name_H-M   P2_1/m

[qzhu2017]

@scottfredericks I have fixed some issue to avoid the using of very large volume factor.
can you also add an option of outputting the information about the used wyckoff sites. In molecular_crystal.py, add an option of verbosity, by default, it output only
Space group requested: 63 generated 63 vol: 169.345755

If verbosity is 1, output the wyckoff information about the multiplicity, letter, and the position of molecular center
like
4a (0.000, 0.193, 0.116)

[scottfredericks]

@qzhu2017
I renamed -v (--volume) to -f (--factor), and added the -v (--verbosity) option to the script. It defaults to 0. If set to greater than 0, it will print the first molecule in each Wyckoff position, and the center of the molecule.

I have also added the -a (--attempts) option to generate more than one structure per run. It still defaults to 1.

I would like to print the molecule's chemical symbol along with the Wyckoff letter, but I am not sure how to obtain it. Are there any other verbosity options you would like?

[qzhu2017]

@scottfredericks Very good design now. I just noticed another error. It failed to generate space group 64 for some reason.

qiangzhu@u-10-240-0-36:~/Desktop/github/crystallography$ python molecular_crystal.py -v 1 -a 10 -f 0.8 -s 63 -n 2
/Users/qiangzhu/miniconda3/lib/python3.6/site-packages/pymatgen/io/cif.py:44: UserWarning: Please install optional dependency pybtex if youwant to extract references from CIF files.
  warnings.warn("Please install optional dependency pybtex if you"
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.29704853 0.25      ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.69965202 0.25      ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.        0.1669475 0.25     ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.28650741 0.25      ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.88176841 0.25      ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.30762532 0.25      ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.18704583 0.25      ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.25537717 0.25      ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.65369338 0.25      ]
Space group requested:  63 generated 36 vol:  130.3625664
4c [0.         0.79558808 0.25      ]

[scottfredericks]

@qzhu2017
It was a careless mistake: I accidentally used mol to generate the final coordinates, instead of mo (which has the orientation already applied). Changing this fixes the issue for spacegroup 63.

Also, with the changes to distance checking, you'll need a volume factor larger than 1, since 1 corresponds to near-perfect orthogonal packing of the molecules. A factor of 4 worked fine for me, with the other parameters left unchanged.

[scottfredericks]

@qzhu2017
I am now able to generate molecular 2d crystals for all layer groups, within a reasonable time (less than a minute for each). I have also used spglib to verify that the output crystals have the correct spacegroup. The script 2dmolecular.py will perform this test automatically for water. Next, I will try H2, NH3, CH4, and C60 (you can just run 'python 2dmolecular.py -e H2', and some other command line options which are used in molecular_crystal.py).

[scottfredericks]

@qzhu2017
There seems to be an issue with generating lattices for 3d molecular crystals. With a volume factor of 1, lattice generation fails for many space groups. With a larger volume factor this is not the case, which leads me to believe it is due to an improper treatment of the minimum vector. I will look into this.

[qzhu2017]

@scottfredericks
it is good to know why it fails. On the other hand, the volume factor for molecular crystals can be larger, I would say some number between 1-2.5 should be okay.