Example: CO2 Crystal ------------------------ .. tip:: The sample input and output files can be found in ``testfiles/1-co2``. In this first example, we will try to predict the crystal structure of :math:`\text{CO}_2`. Run the Search =============== Before proceeding, we should have a **parameter file** for :math:`\text{CO}_2` molecule and an initial guess of crystal shape. - For the parameter file, you can check if it has already be provided in ``mics/charmm36``. In this case, yes: you can find a file ``mics/charmm36/co2.xyz``. That is. In not, you have to build the paramter by yourself. This will be introduced :doc:`eg-nhch32co`. - In addition, you need an initial guess for crystal shape: the number of molecules in a cell (:math:`Z`) and the crystal system. If experimental values are known, we can just use them as input. Otherwise, we have to try different values. Step 1: Prepare an input file named ``sys.inp`` with the following content (other file names are also acceptable): .. code-block:: bash :linenos: :caption: sys.inp sys # Result file name sys.crystal # Crystal components file name 34 # Crystal family: 31-37 20 # ABC argument: Population 5 # ABC argument: Maximum generation 5 # ABC argument: Scout limit 1000 # Number of local minima to be saved 3 # Cutoff of the neighboring cells. For non-polar: 3; polar or charged: 4-6 1.0 # Structures with intramolecular atoms smaller than this distance will be discarded. 0.1 # Structures with density smaller than this value will be discarded. 0.02 # When two energies differ smaller than this value, they (of the same symmetry) are treated as identical. The texts after ``#`` are comments and can be arbitrary. The input will be explained line by line: #. The folder name to save predicted crystal structures. #. The name of the file that contains crystal components, cell initial guess and optionally geometrical information. It will be explained in the next step. #. An integer representing the crystal system that will be kept during the search. The following input are acceptable: +-------------------+------------------+ | Input | Unit | +===================+==================+ | ``31`` | triclinic | +-------------------+------------------+ | ``32`` | monoclinic | +-------------------+------------------+ | ``33`` | orthorhombic | +-------------------+------------------+ | ``34`` | tetragonal | +-------------------+------------------+ | ``35`` | trigonal | +-------------------+------------------+ | ``36`` | hexagonal | +-------------------+------------------+ | ``37`` | cubic | +-------------------+------------------+ One can consider more than cyrstal systems, for example ``35 37`` stands for trigonal and cubic lattice. Note that, during optimization, the symmetry may **raise** (say, triclinic cell may become a cubic one), but will rarely be lowered. #. The population size :math:`SN`. See :doc:`theory` for details. #. The maximum cycle number :math:`g_{\mathrm{max}}`. See :doc:`theory` for details. #. The scout limit :math:`g_{\mathrm{limit}}`. See :doc:`theory` for details. #. The number of crystal structures to be saved. During the calculations, there may be a large number of generated structures. This is a upper limit of the structures to be saved. #. This is a cutoff parameters for short-ranged force evaluations. Usually, ``3`` is suffucient. If higher accuracy is needed, a larger value like ``4`` can be set. However, larger values will lead to more expensive energy evaluation. #. This is a distance in Angstrom. A structure with intramolecular atoms **smaller** than this distance will be discarded. Usually, ``1.0`` is a good choice. #. This is a density in g/cm :sup:`3`. A structure with density **samller** than this value will be discarded. For example, if you know (maybe from experiments) that the crystal density is around 2 g/cm:sup:`3`, then you can set this value to ``1.7``. If you are not sure, simply set it to ``0.1``. #. This is an energy difference in kJ/mol. Two structures with energy difference **smaller** than this value will be discarded. Usually, ``0.02`` is a good choice. Step 2: **Copy** ``misc/charmm36/co2.xyz`` **to the current path.** Then prepare the cluster file named ``sys.crystal`` with the following content: .. code-block:: bash :linenos: :caption: sys.crystal 1 1.01325 co2.xyz 4 6 6 6 90 90 90 * #. The first number is kind of molecules in the crystal. The second number is the pressure in bar. The ambient pressure is ``1.01325`` bar. #. The molecule parameter file name, and its number (:math:`Z`). .. tip:: In the crystal file, the parameter files can be given with **relative or absolute path.** For example, if the water parameter file is at ``/home/you/abcrystal/misc/tip4p.xyz`` and you do not want to copy it, you can simply give the absolute path of this file in ``sys.cluster`` like ``/home/you/abcrystal/misc/tip4p.xyz``. 3. The initial guess of crystal shape: :math:`a,b,c,\alpha,\beta,\gamma` (see the Figure below). Note that, the number does **NOT** have to meet the crystal system constraint. The symmetrization will be treated by ABCrystal. .. image:: pics/cell.png :scale: 50 % :alt: alternate text :align: center 4. A ``*`` is needed. Note that, ``abcrystal`` will also generate ``*.crystal`` files, but it does not end with ``*`` but with some real numbers, for example: .. code-block:: bash :linenos: :caption: A crystal file generated by ABCrystal 1 1.01325000 co2.xyz 4 5.53512686 5.53522121 6.00000000 90.00000000 90.00000000 90.00000000 Energy: -107.74033984 3.93647424 0.95531967 2.35617388 0.00000000 0.00000000 0.00000000 4.68931899 0.95532737 5.49777255 0.49999755 0.00000000 0.49999760 1.68906741 5.32785973 0.78540539 0.00000052 0.49999714 0.49999743 3.38675967 0.95532317 0.78541716 0.49999807 0.49999755 0.00000019 These numbers are rigid coordinates of each molecule in the crystal. You do not need to care about them. Step 3: Run the crystal structure prediction: .. code-block:: bash $ abcrystal sys.inp > sys.out After a few seconds, you will find several new files: - ``sys.out`` The main output file. - ``sys-OPT.crystal`` The most stable crystal structure in ABCrystal format. - ``sys-OPT.xyz`` The most stable crystal structure in XYZ format. - ``sys-OPT.gjf`` The most stable crystal structure in GJF format. - ``sys-OPT.cif`` The most stable crystal structure in CIF format. - ``sys`` A folder containing all crystal strucures, each one having 4 files in ABCrystal, XYZ, GJF, and CIF format, respectively. They are sorted in energy-increasing order, e.g. ``0.xyz`` is lower in energy than ``13.xyz``. - ``abcrystal*.crystal/xyz/gjf/cif`` The file containing the currently found most stable sturcture during the running of ``abcrystal``. You can check the current stable structure before ``abcrystal`` terminates. If it crashes, one can use this ``abcrystal*.crystal`` to start a new search. Check the Output ======================= We can now open ``sys.out``. You can go to ``Optimization Results`` to check the may results: .. code-block:: bash :linenos: :caption: sys.out # Energy Vcell Density Space ID Bravais Good? (kJ/mol) (AA^3) (g/cm^3) Group Lattice 0 -107.7403 169.5847 1.7237 Pa-3 205 cP Y 1 -107.6989 169.5989 1.7235 Pa-3 205 cP Y 2 -107.6493 169.6153 1.7233 Pa-3 205 cP Y 3 -107.4886 169.6683 1.7228 Pa-3 205 cP Y 4 -106.8056 169.9081 1.7204 Pa-3 205 cP Y 5 -106.3587 170.0616 1.7188 Pa-3 205 cP Y 6 -104.1541 168.7669 1.7320 Pca2_1 29 oP Y 7 -103.3376 170.5019 1.7144 Pnnm 58 oP Y 8 -103.3130 170.5116 1.7143 Pnnm 58 oP Y 9 -102.5955 170.8538 1.7109 Pnnm 58 oP Y 10 -96.5469 168.6160 1.7336 Pbca 61 oP Y 11 -94.2686 177.0464 1.6510 P2_1/c 14 mP Y 12 -94.2302 177.0638 1.6508 P2_1/c 14 mP Y 13 -94.0207 177.1590 1.6500 P2_1/c 14 mP Y 14 -92.9351 167.4234 1.7459 P-1 2 aP Y 15 -92.4825 175.7740 1.6630 P2_1/c 14 mP Y Let's see line 3: .. code-block:: bash 0 -107.7403 169.5847 1.7237 Pa-3 205 cP Y - The 1st index means that this structure corresponds to file ``sys/0.crystal/xyz/gjf/cif``. - The 2nd number is the enthalpy: ``-107.7403`` kJ/mol. - The 3rd number is the cell volume: ``169.5847`` Angstrom :sup:`3`. - The 4th number is the density: ``1.7237`` g/cm :sup:`3`. - The 5th string is the space group name international symbol: ``Pa-3``, i.e., :math:`\text{Pa}\bar{3}`. - The 6th number is the space group ID: ``205``. - The 7th string is the Bravais lattice symbol: ``cP`` (primitive cubic lattice). Visualization of Structures ================================== The structures can be visualized in different ways. XYZ file does not contain cell information. It only contains Cartesian coordinates and can be used in, e.g., molecular dynamics. GIF file contains both atomic and cell information and can be visualized with GaussView. All molecules can be shown in a wrapped way. CIF file is the standard crystal file. If can be viewed with Vesta and GaussView, but molecules may not be wrapped. Below is a summary of visualization of structure ``sys/0.xyz/gjf/cif``. .. image:: pics/co2.png :scale: 50 % :alt: alternate text :align: center DFT Calculations ================================== The energy in ABCrystal is caluclated with CHARMM force field. A better energy estimation way is to use density functional theory. Here, we will use CP2K to do a PBE-D3 calculation for 2 crystal structures. 1. Copy ``misc/crystal-cp2k-dft.inp`` to a suitable path, say, ``cp2k-calcs``, and change the name to ``0.inp`` and ``11.inp`` 2. We want to see structure ``0`` and ``11`` energy difference. Copy ``sys/0.cif`` and ``sys/11.cif`` to ``cp2k-calcs``. 3. Open ``0.inp``, and change ``@SET sysname`` line to .. code-block:: bash @SET sysname 0 4. In ``&SUBSYS ... &END SUBSYS``, you cna add corresponding basis set and pseudopotential. In ``&DFT ... &END DFT``, you can change the calculation method. Here, nothing needs to be changed. 5. In ``&CELL_OPT ... &END CELL_OPT``, you can change pressure, symmtery constraint, etc. according to your need. 6. After editting ``0.inp``, run the calculation: .. code-block:: bash cp2k -i 0.inp -o 0.out 7. Run the calculation for ``11`` in a similar way. When the calculations are finished, the optimized atomic position is in ``0/11-pos-1.xyz``, and cell information is in ``0/11.out``. We can compare their results: +------------------------------------+------------------+----------------------+ | Method | CHARMM | PBE-D3 | +====================================+==================+======================+ | E(11)-E(0) (kJ/mol) | 13.47 | 15.91 | +------------------------------------+------------------+----------------------+ | Cell volume 0 (Angstrom :sup:`3`) | 169.58 | 164.33 | +------------------------------------+------------------+----------------------+ | Cell volume 11 (Angstrom :sup:`3`) | 177.05 | 178.39 | +------------------------------------+------------------+----------------------+ The structures are shown below: .. image:: pics/co2011.png :scale: 50 % :alt: alternate text :align: center