Simulating crystals

Dealing with crystals

From CIF to PDB

Well, well, well, I tried to generate a crystal to feed GROMACS only using free software... and that was not easy at the beginning, but then I found mercury and the clouds opened, and I could see the light. Here is the solution I found

what we will dosoftware
First we will get a CIF file. You can download it, or generate it) mercury
Then we will translate it to PDB for the assymetric unit mercury
We will then, eventually, generate the molecules by symmetry operations mercury or SPDB

Generating CIF

CIF file is a crystallographic file with the description of the unit cell and the fractional coordinates writen in the axis set of the cell, i.e. these coordinates are in a non-orthonormal axis set!!!! You can find an example here (is quite self explanatory). You can translate a fractional coordinates set by very simple mathematics operations... but why to do that if a program like MERCURY can do it? So you load the CIF file and say save as PDB. The problem is that you get the minimum number of molecules necessary to construct the whole cell (the assymetric unit) with all molecules that should be generated by symmetry operations (the assymtric unit). To build a whole crystal you must  go to next step.

Applying symmetry operations

To generate all molecules in a crystal go to calculate>packing/slicing> and then select the number of cells you want to build. Please, be aware that in the cif file you MUST have information about the symmetry of your crystal, otherwise you will only "replicate" again the assymetric unit...  to be more clear you must have a line like this one:

_symmetry_space_group_name_H-M     'Cc'

in your cif file. Otherwise you will not generate all molecules in your unit cell!!!

You can alternativelly load the PDB file with Deep View and then you generate the rest of the molecules using the tools "Build crystallographic symmetry". You then select your space group et voilà... you can relax and listen to some music. You can now save your unit cell as PDB. Be aware that you must select "save whole project" option to save all generated molecules.

Repeating the crystallographic cell

You might have done that with MERCURY. Otherwise you can also "repeat" the unit cell as many times as you want using:

  •  gmx genconf -f crystal_one.pdb -o crystal_box.gro -nbox 5 5 5

This will repeat the cell five times in each crystalographic direction... which might not be a good idea. Want to know why? Keep reading.

What if the atom numbering of PDB from CIF and .ITP is different???

Well, the problem here is that normally you will get the forcefield together with a pdb file. That is nice. The numbering of atoms is the same in both files and you just have to careless start your simulations. But, if your starting PDB file is coming from a CIF is probable (indeed usual!) that atom numbering coming from PDB and ForceField is not the same. No worries. you can use the CONTROL file from ANGULA to do the renumbering... please use carefully it is not fully checked!



Gromacs likes Spanish omelette! That means that the box you're building must respect some rules to be able to apply Periodic Boundary Conditions, if the cell is monoclinic or triclinic. In a few words: do not make the cell to high (big z values), is better to make it flat (i.e. big x and y and small z). To be more specific, the rules you can read in the manual to deal with a triclinic box are in page 12:

ay=az=bz=0      (1)

ax>0 , by>0 and cz>0      (2)

|bx|<0.5ax  ,  |cx|<0.5ax     and |cy|<0.5by      (3)

You can find the vectors of your box at the end of the gro file. Be carefull they are written in a strange way!! (to be able to write in an easy way an ortorrombic box)

ax by cz ay az bx bz cx cy

If you do not build a "flatish" crystal you will get a warning

Triclinic box is too skewed

And then GROMACS will write the vectors of the box... if you look at them you will realize that the warning is usually because the rules above are not fullfilled. Specially rule number 3

If you are doing an NPT simulation do not be slacker: do z small and x and y really big, otherwise while adapting the cell shape to get the right pressure at some point you will get the error since the condition 3 above is not fullfilled.


Viewing the molecules

I tried to view the molecules with VMD using pbc cell... however someting weird happened since the box seems not to be the right one. In doing that MERCURY performs better and you can switch more easily between crystallographic views... Also, it is easier to add the labels to the atoms.

Again this is what I have find out testing and maybe is not the best and easiest way to do tings... but it works

Calculating the Bragg diffraction pattern

Well, if you are dealing with crystals, you will probably compare your results with some diffraction pattern. The problem arises when you try to use the standard way to generate a f(q) (like using Sassena, or my fortran code): In order to get nice and sharp Bragg peaks you must simulate a giant crystal... and if you do that calculating the f(q) will take forever...

In order to get the Bragg peaks you can simply simulate a reasonable big crystal and then load it in the Mercury software. Then you click on calculations >> powder pattern et voilà. However, two warnings:

  • You are making one assumption: since waht the program is doing is taking your simulation as the unit cell and repeating it to infinity... that means that you are generating periodicities that are not real. However, for big enough cells, you will eventually see peaks only at vey small q's where no Bragg peaks are expected
  • GROMACS changes names: Maybe I am doing something wrong... but when using trjconv to generate the PDB files, the program misteriously changes the atom names at the end of the PDB file... and therefore the diffractogram is meaningless. Of course I wrote a small fortran code to "repair" the pdb files, just cutting the final name of the atoms.

 Dealing with disorder in crystals

If you have disorder in the crystal, you are going to have in the PDB generated from the CIF file a "pseuo-molecule" with all possibilities. For example, if a site is 25% occupied by an H, and 75% by an O, you will obtain a (fake) molecule with both possibilities. 

We could think: "ahá, i choose an option and I generate the crystal with a "correct" molecule"... yes, but if you do it before constructing the crystal, you will have no disorder... and this will probably cause an anysotropy that will destroy your crystal. And we don't want that...

An alternative is to generate the crystal and randomly delete the atoms so that you get a right molecule... but with random orientations... and you can do that with my software "anangle" with the option "rndstructure". Just follow the dialog and you will (eventually) generate a structure capturing the randomness of the molecular orientation (what a cool phrase to finish the paragraph!).




(many thanks to Dr. Arul Murugan that helped me in this part when I was lost...)