Expanding documentation and adding algo selection

README.md

@@ -12,16 +12,7 @@ scalar fields. Den2Obj supports VASP charge files such as CHGCAR and PARCHG,
 Gaussian .cube files as well as its own .d2o format.
 
 ## Example images
-![3D Reaction-Diffusion system](img/reac_diff_3d_network_small.png)
-
-*The isosurface above represents the concentration profile of a
- reaction-diffusion system in 3D using Gray-Scott kinetics. The isosurface has
- been generated using den2obj and rendered using [Blender]
- (https://www.blender.org/).*
-
-![Molecular orbitals of CO](img/CO_mos.jpg)
-
-*The isosurfaces of the first 10 molecular orbitals of the CO molecule.*
+![Canonical valence orbitals of CH4](docs/_static/img/ch4_valence_orbitals.png)
 
 ## Compilation instructions
 

docs/background.rst

@@ -53,19 +53,26 @@ isosurface intersects each edge is determined using trilinear interpolation.
    **Figure 2**: The 256 possible ways that a scalar field can interact with
    each cube.
 
-To repeat on this. First the space wherein the function is evaluated is
-discretized into small cubes. Next, for each vertex on these cubes the function
-is evaluated and it is determined whether the value for the function at each of
-the vertices is larger or smaller than the isovalue. This leads to a number of
-edge intersections for each cube from which the polygons (triangles) for each
-cube can be established. As a final step, all polygons are gathered to form the
-threedimensional isosurface. It should be noted that this algorithm can be
-executed in a highly efficient fashion using trivial parallellization as the
-result for each cube is completely independent from all the other cubes. It
-turns out that the generation of the scalar field, which by itself is also a
-highly parallellizable step, is typically the most time-consuming.
+Thus the procedure is as follows:
+
+1. The space wherein the function is evaluated is discretized into small
+   cubes. 
+2. For each vertex on these cubes the function is evaluated and it is
+   determined whether the value for the function at each of the vertices is
+   larger or smaller than the isovalue. This leads to a number of edge
+   intersections for each cube from which the polygons (triangles) for each cube
+   can be established. 
+3. All polygons are gathered to form the threedimensional isosurface. 
+
+It should be noted that this algorithm can be executed in a highly efficient
+fashion using trivial parallellization as the result for each cube is
+completely independent from all the other cubes. It turns out that the
+generation of the scalar field, which by itself is also a highly
+parallellizable step, is typically the most time-consuming.
 
 The quality of the isosurface depends on the resolution of the grid used to
 sample the scalar field. The finer the grid (or mesh) used, the better the
-quality of the isosurface. Of course, a finer mesh comes at the expense of more
-values that need to evaluated in the algorithm.
+quality of the isosurface. Thus, rather than using cubes to sample the space,
+one can use tetrahedra. The algorithm is essentially the same as the one
+as described above, yet each cube is once more divided into six tetrahedra
+for which the above algorithm is executed.

docs/community_guidelines.rst

@@ -0,0 +1,14 @@
+.. _community_guidelines:
+.. index:: Community Guidelines
+
+Community guidelines
+********************
+
+* Contributions to :program:`Den2Obj` are always welcome and appreciated. Before doing
+  so, please first read the `CONTRIBUTING <https://github.com/ifilot/den2obj/blob/master/CONTRIBUTING.md>`_
+  guide.
+* For reporting issues or problems with the software, you are kindly invited to
+  `to open a new issue on Github with the bug tag <https://github.com/ifilot/den2obj/issues/new?labels=bug>`_.
+* If you seek support in using :program:`Den2Obj`, please
+  `open an issue with the question tag <https://github.com/ifilot/edp/issues/new?labels=question>`_.
+* If you wish to contact the developers, please send an e-mail to i.a.w.filot@tue.nl.

docs/d2o_fileformat.rst

@@ -60,7 +60,7 @@ To convert a ``CHGCAR`` file to ``.d2o`` file format, run the following command:
    cases, this corresponds to the `LZMA type of compression <https://en.wikipedia.org/wiki/Lempel%E2%80%93Ziv%E2%80%93Markov_chain_algorithm>`_.
 
 Protocol tokens
-===============
+---------------
 
 1. GZIP compression
 2. LZMA compression

docs/index.rst

@@ -25,6 +25,8 @@ such as `Blender <https://www.blender.org/>`_ to produce stunning visuals
 of your research. Below, a few examples are shown. See also the
 :ref:`Examples<Examples>` for a nice overview.
 
+.. figure:: _static/img/ch4_valence_orbitals.png
+   :alt: Canonical valence orbitals of CH4
 
 .. toctree::
    :maxdepth: 2
@@ -36,6 +38,7 @@ of your research. Below, a few examples are shown. See also the
    examples
    user_interface
    d2o_fileformat
+   community_guidelines
    faq
 
 Indices and tables

docs/tutorial.rst

@@ -116,3 +116,147 @@ we have used the settings as can be seen in the figure below.
 
 The only step that remains is to render the image, which will give the image
 as shown at the start of this section.
+
+Benzene highest occupied molecular orbital
+------------------------------------------
+
+In this tutorial, we will be reproducing the Figure as shown below.
+
+.. figure:: _static/img/tutorials/benzene_homo_result.png
+   :alt: Benzene highest occupied molecular orbital
+
+To generate the scalar field, run::
+
+    ./den2obj -g benzene_homo -o benzene_homo.d2o
+
+Next, the isosurface is generated. An isovalue of 0.03 is chosen. Because
+molecular orbital have positive and negative lobes, we use the ``-d`` tag
+to create both isosurfaces. Furthermore, we center the object by defining
+``-c`` and we explicitly ask to use the marching tetrahedra algorithm
+via ``--algo marching-tetrahedra``::
+
+    ./den2obj -i benzene_homo.d2o -o benzene_homo.ply -v 0.03 -c -d --algo marching-tetrahedra
+
+The following output (or similar) is generated::
+
+    --------------------------------------------------------------
+    Executing DEN2OBJ v.1.1.0
+    Author:  Ivo Filot <i.a.w.filot@tue.nl>
+    Website: https://den2obj.imc-tue.nl
+    Github:  https://github.com/ifilot/den2obj
+    --------------------------------------------------------------
+    Opening benzene_homo.d2o as D2O binary file
+    Recognizing floating point size: 4 bytes.
+    Reading 11415368 bytes from file.
+    Building LZMA decompressor
+    Decompressed data
+    Done reading D2O binary file
+    Read 3375000 values.
+    Using isovalue: 0.03
+    Lowest value in scalar field: -0.25383
+    Highest value in scalar field: 0.25383
+    Calculating normal vectors using two-point stencil
+    Centering structure
+    Writing mesh as Standford Triangle Format file (.ply).
+    Writing as Stanford (.ply) file: benzene_homo_pos.ply (4560.3kb).
+    Identified 59512 faces.
+    Calculating normal vectors using two-point stencil
+    Centering structure
+    Writing mesh as Standford Triangle Format file (.ply).
+    Writing as Stanford (.ply) file: benzene_homo_neg.ply (1454.4kb).
+    -------------------------------------------------------------------------------
+    Done in 2.17096 seconds.
+
+Observe that two isosurfaces are created and stored as ``.ply`` files:
+
+* benzene_homo_pos.ply
+* benzene_homo_neg.ply
+
+Importing these two files in Blender gives us the following result
+
+.. figure:: _static/img/tutorials/benzene_homo_blender_01.JPG
+   :alt: HOMO orbital of benzene imported into Blender
+
+Of course, this result is rather blend and we would like to add
+the positions of the carbon and hydrogen atoms and the bonds between
+them. For this, we are going to use the hardcoded Python script as shown
+below which we can readily execute in Blender
+
+.. code-block:: python
+
+    import bpy
+    import numpy as np
+
+    def main():
+        # define molecule
+        mol = []
+        
+        mol.append(['C',  0.0000000015, -1.3868467444, 0.0000000000])
+        mol.append(['C',  1.2010445126, -0.6934233709, 0.0000000000])
+        mol.append(['C',  1.2010445111,  0.6934233735, 0.0000000000])
+        mol.append(['C', -0.0000000015,  1.3868467444, 0.0000000000])
+        mol.append(['C', -1.2010445126,  0.6934233709, 0.0000000000])
+        mol.append(['C', -1.2010445111, -0.6934233735, 0.0000000000])
+        mol.append(['H',  0.0000000027, -2.4694205285, 0.0000000000])
+        mol.append(['H',  2.1385809117, -1.2347102619, 0.0000000000])
+        mol.append(['H',  2.1385809090,  1.2347102666, 0.0000000000])
+        mol.append(['H', -0.0000000027,  2.4694205285, 0.0000000000])
+        mol.append(['H', -2.1385809117,  1.2347102619, 0.0000000000])
+        mol.append(['H', -2.1385809090, -1.2347102666, 0.0000000000])
+
+        build_atoms(mol)
+        build_bonds(mol)
+
+    def build_atoms(molecule):
+        for i,at in enumerate(molecule):
+            sc = 0.4 if at[0] == 'C' else 0.3
+            obj = bpy.ops.mesh.primitive_ico_sphere_add(radius=sc, 
+                                                        enter_editmode=False, 
+                                                        align='WORLD', 
+                                                        location=(at[1], at[2], at[3]), 
+                                                        scale=(1, 1, 1),
+                                                        subdivisions=3)
+            bpy.context.selected_objects[0].name = at[0] + str(i+1)
+            bpy.ops.object.shade_smooth()
+
+    def build_bonds(molecule):
+        z = np.array([0,0,1])
+        for i,at1 in enumerate(molecule):
+            p1 = np.array([at1[1], at1[2], at1[3]])
+            for j,at2 in enumerate(molecule[i+1:]):
+                p2 = np.array([at2[1], at2[2], at2[3]])
+                d = np.linalg.norm(p1-p2)
+                
+                if d < 1.5:
+                    ctr = (p1 + p2) / 2
+                    angle = np.arccos(np.dot(z, (p2-p1) / d))
+                    axis = np.cross(z, (p2-p1) / d)
+                    axis /= np.linalg.norm(axis)
+                    
+                    bpy.ops.mesh.primitive_cylinder_add(radius=0.2, 
+                                                        depth=1, 
+                                                        enter_editmode=False, 
+                                                        align='WORLD', 
+                                                        location=(ctr[0], ctr[1], ctr[2]), 
+                                                        scale=(1, 1, 1))
+                    bpy.context.object.rotation_mode = 'AXIS_ANGLE'
+                    bpy.context.object.rotation_axis_angle[0] = angle
+                    bpy.context.object.rotation_axis_angle[1:4] = axis
+                    bpy.context.selected_objects[0].name = 'bond' + str(i+1) + '-' + str(j+1)
+                    bpy.ops.object.shade_smooth()
+
+
+    if __name__ == '__main__':
+        main()
+
+This will generate all the atoms and bonds between them. Next,
+materials are assigned to all atoms and bonds and the final result
+looks as seen in the image below.
+
+.. figure:: _static/img/tutorials/benzene_homo_blender_02.JPG
+   :alt: HOMO orbital of benzene together with the atoms and bonds
+
+Finally, we can render the scene to create a nice picture of the molecular orbital.
+
+.. figure:: _static/img/tutorials/benzene_homo_result.png
+   :alt: Benzene highest occupied molecular orbital

docs/user_interface.rst

@@ -27,7 +27,7 @@ Isosurface generation
 
 To perform an isosurface creation, one simply runs::
 
-    ./den2obj -i <path-to-scalarfield> -o <mesh-file> -v <isovalue> [-c] [-d]
+    ./den2obj -i <path-to-scalarfield> -o <mesh-file> -v <isovalue> [-c] [-d] [-a <algo>]
 
 wherein ``path-to-scalarfield`` is a scalar field file format of either
 of the following types
@@ -55,6 +55,16 @@ and ``isovalue`` the isovalue for the isosurface.
   that the center of the unit cell is located at the origin.
 * The argument ``-d`` is optional. If used, both signs of the isosurface will be
   rendered. This is useful when rendering wave functions rather than densities.
+* The argument ``--algo`` is optional. It can be set to either ``marching-cubes``
+  or to ``marching-tetrahedra``. The default option is ``marching-cubes``.
+
+.. figure:: _static/img/tutorials/genus2_tetrahedra_versus_cubes.png
+   :alt: Marching tetrahedra versus marching cubes
+
+   Marching tetrahedra (left) versus marching cubes (right). Edges are explicitly
+   visualized. Observe that whilst using the same resolution of the dataset,
+   when using marching tetrahedra more facets are created and thus a higher
+   resolution isosurface.
 
 Filetype conversion
 -------------------
@@ -100,4 +110,9 @@ directive. For example, to force BZIP2 type of compression, one runs::
     ./den2obj -g genus2 -o genus2.d2o -a bzip2
 
 When no ``-a`` is provided, automatically the best compression algorithm is used by
-checking the inflation ratio of all possible compression algorithms.
+checking the inflation ratio of all possible compression algorithms.
+
+The following datasets are available:
+* ``genus2``
+* ``benzene_homo``
+* ``benzene_lumo``

src/den2obj.cpp

@@ -103,7 +103,11 @@ int main(int argc, char* argv[]) {
             }
         }
 
-        // check if a generation is requested
+        /************************************************************************************************************
+         *
+         * SCALAR FIELD GENERATION
+         *
+         ***********************************************************************************************************/
         if(!arg_generator.getValue().empty()) {
             Generator gen;
 
@@ -143,6 +147,12 @@ int main(int argc, char* argv[]) {
                 throw std::runtime_error("Cannot interpret input file format. Please check the filename.");
             }
 
+            /************************************************************************************************************
+             *
+             * DATASET CONVERSION (TRANSFORMATION)
+             *
+             ***********************************************************************************************************/
+
             // check whether only a conversion is requested, if not, continue
             if(arg_t.getValue()) {
                 if(output_filename.substr(output_filename.size()-4) == ".d2o") {
@@ -162,6 +172,12 @@ int main(int argc, char* argv[]) {
                 }
 
             } else {
+
+                /************************************************************************************************************
+                *
+                * ISOSURFACE GENERATION
+                *
+                ***********************************************************************************************************/
                 fpt isovalue = arg_isovalue.getValue();
 
                 // automatically convert the isovalue to a positive number if d is set
@@ -177,7 +193,16 @@ int main(int argc, char* argv[]) {
 
                 // construct isosurface generator object
                 IsoSurface is(sf.get());
-                is.marching_cubes(isovalue);
+
+                if(arg_algo.getValue().empty() || arg_algo.getValue() == "marching-cubes") {
+                    std::cout << "Using marching cubes algorithm for isosurface generation." << std::endl;
+                    is.marching_cubes(isovalue);    
+                } else if(arg_algo.getValue() == "marching-tetrahedra") {
+                    std::cout << "Using marching tetrahedra algorithm for isosurface generation." << std::endl;
+                    is.marching_tetrahedra(isovalue);
+                } else {
+                    throw std::runtime_error("Invalid isosurface generation algo: " + arg_algo.getValue());
+                }
 
                 // store path to extract filename
                 boost::filesystem::path path(output_filename);