diff --git a/README.md b/README.md index be905e1..0cd0322 100644 --- a/README.md +++ b/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 diff --git a/docs/_static/img/ch4_valence_orbitals.png b/docs/_static/img/ch4_valence_orbitals.png new file mode 100644 index 0000000..0e0b1a6 Binary files /dev/null and b/docs/_static/img/ch4_valence_orbitals.png differ diff --git a/docs/_static/img/co_adsorption.png b/docs/_static/img/co_adsorption.png new file mode 100644 index 0000000..6c6fcd0 Binary files /dev/null and b/docs/_static/img/co_adsorption.png differ diff --git a/docs/_static/img/tutorials/benzene_homo_blender_01.jpg b/docs/_static/img/tutorials/benzene_homo_blender_01.jpg new file mode 100644 index 0000000..497ef55 Binary files /dev/null and b/docs/_static/img/tutorials/benzene_homo_blender_01.jpg differ diff --git a/docs/_static/img/tutorials/benzene_homo_blender_02.jpg b/docs/_static/img/tutorials/benzene_homo_blender_02.jpg new file mode 100644 index 0000000..50f8b81 Binary files /dev/null and b/docs/_static/img/tutorials/benzene_homo_blender_02.jpg differ diff --git a/docs/_static/img/tutorials/benzene_homo_result.png b/docs/_static/img/tutorials/benzene_homo_result.png new file mode 100644 index 0000000..42935a4 Binary files /dev/null and b/docs/_static/img/tutorials/benzene_homo_result.png differ diff --git a/docs/_static/img/tutorials/genus2_tetrahedra_versus_cubes.png b/docs/_static/img/tutorials/genus2_tetrahedra_versus_cubes.png new file mode 100644 index 0000000..e52afc9 Binary files /dev/null and b/docs/_static/img/tutorials/genus2_tetrahedra_versus_cubes.png differ diff --git a/docs/background.rst b/docs/background.rst index 39bddf6..2e1deda 100644 --- a/docs/background.rst +++ b/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. \ No newline at end of file diff --git a/docs/community_guidelines.rst b/docs/community_guidelines.rst new file mode 100644 index 0000000..491d6ee --- /dev/null +++ b/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 `_ + 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 `_. +* If you seek support in using :program:`Den2Obj`, please + `open an issue with the question tag `_. +* If you wish to contact the developers, please send an e-mail to i.a.w.filot@tue.nl. diff --git a/docs/d2o_fileformat.rst b/docs/d2o_fileformat.rst index 4c59945..65cf8c6 100644 --- a/docs/d2o_fileformat.rst +++ b/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 `_. Protocol tokens -=============== +--------------- 1. GZIP compression 2. LZMA compression diff --git a/docs/index.rst b/docs/index.rst index 49c6986..bac7080 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -25,6 +25,8 @@ such as `Blender `_ to produce stunning visuals of your research. Below, a few examples are shown. See also the :ref:`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 diff --git a/docs/tutorial.rst b/docs/tutorial.rst index ec50317..735a759 100644 --- a/docs/tutorial.rst +++ b/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 + 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 \ No newline at end of file diff --git a/docs/user_interface.rst b/docs/user_interface.rst index 57838a6..7ab2703 100644 --- a/docs/user_interface.rst +++ b/docs/user_interface.rst @@ -27,7 +27,7 @@ Isosurface generation To perform an isosurface creation, one simply runs:: - ./den2obj -i -o -v [-c] [-d] + ./den2obj -i -o -v [-c] [-d] [-a ] 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. \ No newline at end of file +checking the inflation ratio of all possible compression algorithms. + +The following datasets are available: +* ``genus2`` +* ``benzene_homo`` +* ``benzene_lumo`` \ No newline at end of file diff --git a/src/den2obj.cpp b/src/den2obj.cpp index 50c9ed7..aa15646 100644 --- a/src/den2obj.cpp +++ b/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);