io_mesh_atomic/pdb_import.py

   1 # SPDX-License-Identifier: GPL-2.0-or-later
   2
   3 import os
   4 import bpy
   5 import bmesh
   6 from math import pi, cos, sin, sqrt, ceil
   7 from mathutils import Vector, Matrix
   8 from copy import copy
   9
  10 # -----------------------------------------------------------------------------
  11 #                                                  Atom, stick and element data
  12
  13
  14 # This is a list that contains some data of all possible elements. The structure
  15 # is as follows:
  16 #
  17 # 1, "Hydrogen", "H", [0.0,0.0,1.0], 0.32, 0.32, 0.32 , -1 , 1.54   means
  18 #
  19 # No., name, short name, color, radius (used), radius (covalent), radius (atomic),
  20 #
  21 # charge state 1, radius (ionic) 1, charge state 2, radius (ionic) 2, ... all
  22 # charge states for any atom are listed, if existing.
  23 # The list is fixed and cannot be changed ... (see below)
  24
  25 ELEMENTS_DEFAULT = (
  26 ( 1,      "Hydrogen",        "H",  (  1.0,   1.0,   1.0, 1.0), 0.32, 0.32, 0.79 , -1 , 1.54 ),
  27 ( 2,        "Helium",       "He",  ( 0.85,   1.0,   1.0, 1.0), 0.93, 0.93, 0.49 ),
  28 ( 3,       "Lithium",       "Li",  (  0.8,  0.50,   1.0, 1.0), 1.23, 1.23, 2.05 ,  1 , 0.68 ),
  29 ( 4,     "Beryllium",       "Be",  ( 0.76,   1.0,   0.0, 1.0), 0.90, 0.90, 1.40 ,  1 , 0.44 ,  2 , 0.35 ),
  30 ( 5,         "Boron",        "B",  (  1.0,  0.70,  0.70, 1.0), 0.82, 0.82, 1.17 ,  1 , 0.35 ,  3 , 0.23 ),
  31 ( 6,        "Carbon",        "C",  ( 0.56,  0.56,  0.56, 1.0), 0.77, 0.77, 0.91 , -4 , 2.60 ,  4 , 0.16 ),
  32 ( 7,      "Nitrogen",        "N",  ( 0.18,  0.31,  0.97, 1.0), 0.75, 0.75, 0.75 , -3 , 1.71 ,  1 , 0.25 ,  3 , 0.16 ,  5 , 0.13 ),
  33 ( 8,        "Oxygen",        "O",  (  1.0,  0.05,  0.05, 1.0), 0.73, 0.73, 0.65 , -2 , 1.32 , -1 , 1.76 ,  1 , 0.22 ,  6 , 0.09 ),
  34 ( 9,      "Fluorine",        "F",  ( 0.56,  0.87,  0.31, 1.0), 0.72, 0.72, 0.57 , -1 , 1.33 ,  7 , 0.08 ),
  35 (10,          "Neon",       "Ne",  ( 0.70,  0.89,  0.96, 1.0), 0.71, 0.71, 0.51 ,  1 , 1.12 ),
  36 (11,        "Sodium",       "Na",  ( 0.67,  0.36,  0.94, 1.0), 1.54, 1.54, 2.23 ,  1 , 0.97 ),
  37 (12,     "Magnesium",       "Mg",  ( 0.54,   1.0,   0.0, 1.0), 1.36, 1.36, 1.72 ,  1 , 0.82 ,  2 , 0.66 ),
  38 (13,     "Aluminium",       "Al",  ( 0.74,  0.65,  0.65, 1.0), 1.18, 1.18, 1.82 ,  3 , 0.51 ),
  39 (14,       "Silicon",       "Si",  ( 0.94,  0.78,  0.62, 1.0), 1.11, 1.11, 1.46 , -4 , 2.71 , -1 , 3.84 ,  1 , 0.65 ,  4 , 0.42 ),
  40 (15,    "Phosphorus",        "P",  (  1.0,  0.50,   0.0, 1.0), 1.06, 1.06, 1.23 , -3 , 2.12 ,  3 , 0.44 ,  5 , 0.35 ),
  41 (16,        "Sulfur",        "S",  (  1.0,   1.0,  0.18, 1.0), 1.02, 1.02, 1.09 , -2 , 1.84 ,  2 , 2.19 ,  4 , 0.37 ,  6 , 0.30 ),
  42 (17,      "Chlorine",       "Cl",  ( 0.12,  0.94,  0.12, 1.0), 0.99, 0.99, 0.97 , -1 , 1.81 ,  5 , 0.34 ,  7 , 0.27 ),
  43 (18,         "Argon",       "Ar",  ( 0.50,  0.81,  0.89, 1.0), 0.98, 0.98, 0.88 ,  1 , 1.54 ),
  44 (19,     "Potassium",        "K",  ( 0.56,  0.25,  0.83, 1.0), 2.03, 2.03, 2.77 ,  1 , 0.81 ),
  45 (20,       "Calcium",       "Ca",  ( 0.23,   1.0,   0.0, 1.0), 1.74, 1.74, 2.23 ,  1 , 1.18 ,  2 , 0.99 ),
  46 (21,      "Scandium",       "Sc",  ( 0.90,  0.90,  0.90, 1.0), 1.44, 1.44, 2.09 ,  3 , 0.73 ),
  47 (22,      "Titanium",       "Ti",  ( 0.74,  0.76,  0.78, 1.0), 1.32, 1.32, 2.00 ,  1 , 0.96 ,  2 , 0.94 ,  3 , 0.76 ,  4 , 0.68 ),
  48 (23,      "Vanadium",        "V",  ( 0.65,  0.65,  0.67, 1.0), 1.22, 1.22, 1.92 ,  2 , 0.88 ,  3 , 0.74 ,  4 , 0.63 ,  5 , 0.59 ),
  49 (24,      "Chromium",       "Cr",  ( 0.54,   0.6,  0.78, 1.0), 1.18, 1.18, 1.85 ,  1 , 0.81 ,  2 , 0.89 ,  3 , 0.63 ,  6 , 0.52 ),
  50 (25,     "Manganese",       "Mn",  ( 0.61,  0.47,  0.78, 1.0), 1.17, 1.17, 1.79 ,  2 , 0.80 ,  3 , 0.66 ,  4 , 0.60 ,  7 , 0.46 ),
  51 (26,          "Iron",       "Fe",  ( 0.87,   0.4,   0.2, 1.0), 1.17, 1.17, 1.72 ,  2 , 0.74 ,  3 , 0.64 ),
  52 (27,        "Cobalt",       "Co",  ( 0.94,  0.56,  0.62, 1.0), 1.16, 1.16, 1.67 ,  2 , 0.72 ,  3 , 0.63 ),
  53 (28,        "Nickel",       "Ni",  ( 0.31,  0.81,  0.31, 1.0), 1.15, 1.15, 1.62 ,  2 , 0.69 ),
  54 (29,        "Copper",       "Cu",  ( 0.78,  0.50,   0.2, 1.0), 1.17, 1.17, 1.57 ,  1 , 0.96 ,  2 , 0.72 ),
  55 (30,          "Zinc",       "Zn",  ( 0.49,  0.50,  0.69, 1.0), 1.25, 1.25, 1.53 ,  1 , 0.88 ,  2 , 0.74 ),
  56 (31,       "Gallium",       "Ga",  ( 0.76,  0.56,  0.56, 1.0), 1.26, 1.26, 1.81 ,  1 , 0.81 ,  3 , 0.62 ),
  57 (32,     "Germanium",       "Ge",  (  0.4,  0.56,  0.56, 1.0), 1.22, 1.22, 1.52 , -4 , 2.72 ,  2 , 0.73 ,  4 , 0.53 ),
  58 (33,       "Arsenic",       "As",  ( 0.74,  0.50,  0.89, 1.0), 1.20, 1.20, 1.33 , -3 , 2.22 ,  3 , 0.58 ,  5 , 0.46 ),
  59 (34,      "Selenium",       "Se",  (  1.0,  0.63,   0.0, 1.0), 1.16, 1.16, 1.22 , -2 , 1.91 , -1 , 2.32 ,  1 , 0.66 ,  4 , 0.50 ,  6 , 0.42 ),
  60 (35,       "Bromine",       "Br",  ( 0.65,  0.16,  0.16, 1.0), 1.14, 1.14, 1.12 , -1 , 1.96 ,  5 , 0.47 ,  7 , 0.39 ),
  61 (36,       "Krypton",       "Kr",  ( 0.36,  0.72,  0.81, 1.0), 1.31, 1.31, 1.24 ),
  62 (37,      "Rubidium",       "Rb",  ( 0.43,  0.18,  0.69, 1.0), 2.16, 2.16, 2.98 ,  1 , 1.47 ),
  63 (38,     "Strontium",       "Sr",  (  0.0,   1.0,   0.0, 1.0), 1.91, 1.91, 2.45 ,  2 , 1.12 ),
  64 (39,       "Yttrium",        "Y",  ( 0.58,   1.0,   1.0, 1.0), 1.62, 1.62, 2.27 ,  3 , 0.89 ),
  65 (40,     "Zirconium",       "Zr",  ( 0.58,  0.87,  0.87, 1.0), 1.45, 1.45, 2.16 ,  1 , 1.09 ,  4 , 0.79 ),
  66 (41,       "Niobium",       "Nb",  ( 0.45,  0.76,  0.78, 1.0), 1.34, 1.34, 2.08 ,  1 , 1.00 ,  4 , 0.74 ,  5 , 0.69 ),
  67 (42,    "Molybdenum",       "Mo",  ( 0.32,  0.70,  0.70, 1.0), 1.30, 1.30, 2.01 ,  1 , 0.93 ,  4 , 0.70 ,  6 , 0.62 ),
  68 (43,    "Technetium",       "Tc",  ( 0.23,  0.61,  0.61, 1.0), 1.27, 1.27, 1.95 ,  7 , 0.97 ),
  69 (44,     "Ruthenium",       "Ru",  ( 0.14,  0.56,  0.56, 1.0), 1.25, 1.25, 1.89 ,  4 , 0.67 ),
  70 (45,       "Rhodium",       "Rh",  ( 0.03,  0.49,  0.54, 1.0), 1.25, 1.25, 1.83 ,  3 , 0.68 ),
  71 (46,     "Palladium",       "Pd",  (  0.0,  0.41,  0.52, 1.0), 1.28, 1.28, 1.79 ,  2 , 0.80 ,  4 , 0.65 ),
  72 (47,        "Silver",       "Ag",  ( 0.75,  0.75,  0.75, 1.0), 1.34, 1.34, 1.75 ,  1 , 1.26 ,  2 , 0.89 ),
  73 (48,       "Cadmium",       "Cd",  (  1.0,  0.85,  0.56, 1.0), 1.48, 1.48, 1.71 ,  1 , 1.14 ,  2 , 0.97 ),
  74 (49,        "Indium",       "In",  ( 0.65,  0.45,  0.45, 1.0), 1.44, 1.44, 2.00 ,  3 , 0.81 ),
  75 (50,           "Tin",       "Sn",  (  0.4,  0.50,  0.50, 1.0), 1.41, 1.41, 1.72 , -4 , 2.94 , -1 , 3.70 ,  2 , 0.93 ,  4 , 0.71 ),
  76 (51,      "Antimony",       "Sb",  ( 0.61,  0.38,  0.70, 1.0), 1.40, 1.40, 1.53 , -3 , 2.45 ,  3 , 0.76 ,  5 , 0.62 ),
  77 (52,     "Tellurium",       "Te",  ( 0.83,  0.47,   0.0, 1.0), 1.36, 1.36, 1.42 , -2 , 2.11 , -1 , 2.50 ,  1 , 0.82 ,  4 , 0.70 ,  6 , 0.56 ),
  78 (53,        "Iodine",        "I",  ( 0.58,   0.0,  0.58, 1.0), 1.33, 1.33, 1.32 , -1 , 2.20 ,  5 , 0.62 ,  7 , 0.50 ),
  79 (54,         "Xenon",       "Xe",  ( 0.25,  0.61,  0.69, 1.0), 1.31, 1.31, 1.24 ),
  80 (55,       "Caesium",       "Cs",  ( 0.34,  0.09,  0.56, 1.0), 2.35, 2.35, 3.35 ,  1 , 1.67 ),
  81 (56,        "Barium",       "Ba",  (  0.0,  0.78,   0.0, 1.0), 1.98, 1.98, 2.78 ,  1 , 1.53 ,  2 , 1.34 ),
  82 (57,     "Lanthanum",       "La",  ( 0.43,  0.83,   1.0, 1.0), 1.69, 1.69, 2.74 ,  1 , 1.39 ,  3 , 1.06 ),
  83 (58,        "Cerium",       "Ce",  (  1.0,   1.0,  0.78, 1.0), 1.65, 1.65, 2.70 ,  1 , 1.27 ,  3 , 1.03 ,  4 , 0.92 ),
  84 (59,  "Praseodymium",       "Pr",  ( 0.85,   1.0,  0.78, 1.0), 1.65, 1.65, 2.67 ,  3 , 1.01 ,  4 , 0.90 ),
  85 (60,     "Neodymium",       "Nd",  ( 0.78,   1.0,  0.78, 1.0), 1.64, 1.64, 2.64 ,  3 , 0.99 ),
  86 (61,    "Promethium",       "Pm",  ( 0.63,   1.0,  0.78, 1.0), 1.63, 1.63, 2.62 ,  3 , 0.97 ),
  87 (62,      "Samarium",       "Sm",  ( 0.56,   1.0,  0.78, 1.0), 1.62, 1.62, 2.59 ,  3 , 0.96 ),
  88 (63,      "Europium",       "Eu",  ( 0.38,   1.0,  0.78, 1.0), 1.85, 1.85, 2.56 ,  2 , 1.09 ,  3 , 0.95 ),
  89 (64,    "Gadolinium",       "Gd",  ( 0.27,   1.0,  0.78, 1.0), 1.61, 1.61, 2.54 ,  3 , 0.93 ),
  90 (65,       "Terbium",       "Tb",  ( 0.18,   1.0,  0.78, 1.0), 1.59, 1.59, 2.51 ,  3 , 0.92 ,  4 , 0.84 ),
  91 (66,    "Dysprosium",       "Dy",  ( 0.12,   1.0,  0.78, 1.0), 1.59, 1.59, 2.49 ,  3 , 0.90 ),
  92 (67,       "Holmium",       "Ho",  (  0.0,   1.0,  0.61, 1.0), 1.58, 1.58, 2.47 ,  3 , 0.89 ),
  93 (68,        "Erbium",       "Er",  (  0.0,  0.90,  0.45, 1.0), 1.57, 1.57, 2.45 ,  3 , 0.88 ),
  94 (69,       "Thulium",       "Tm",  (  0.0,  0.83,  0.32, 1.0), 1.56, 1.56, 2.42 ,  3 , 0.87 ),
  95 (70,     "Ytterbium",       "Yb",  (  0.0,  0.74,  0.21, 1.0), 1.74, 1.74, 2.40 ,  2 , 0.93 ,  3 , 0.85 ),
  96 (71,      "Lutetium",       "Lu",  (  0.0,  0.67,  0.14, 1.0), 1.56, 1.56, 2.25 ,  3 , 0.85 ),
  97 (72,       "Hafnium",       "Hf",  ( 0.30,  0.76,   1.0, 1.0), 1.44, 1.44, 2.16 ,  4 , 0.78 ),
  98 (73,      "Tantalum",       "Ta",  ( 0.30,  0.65,   1.0, 1.0), 1.34, 1.34, 2.09 ,  5 , 0.68 ),
  99 (74,      "Tungsten",        "W",  ( 0.12,  0.58,  0.83, 1.0), 1.30, 1.30, 2.02 ,  4 , 0.70 ,  6 , 0.62 ),
 100 (75,       "Rhenium",       "Re",  ( 0.14,  0.49,  0.67, 1.0), 1.28, 1.28, 1.97 ,  4 , 0.72 ,  7 , 0.56 ),
 101 (76,        "Osmium",       "Os",  ( 0.14,   0.4,  0.58, 1.0), 1.26, 1.26, 1.92 ,  4 , 0.88 ,  6 , 0.69 ),
 102 (77,       "Iridium",       "Ir",  ( 0.09,  0.32,  0.52, 1.0), 1.27, 1.27, 1.87 ,  4 , 0.68 ),
 103 (78,      "Platinum",       "Pt",  ( 0.81,  0.81,  0.87, 1.0), 1.30, 1.30, 1.83 ,  2 , 0.80 ,  4 , 0.65 ),
 104 (79,          "Gold",       "Au",  (  1.0,  0.81,  0.13, 1.0), 1.34, 1.34, 1.79 ,  1 , 1.37 ,  3 , 0.85 ),
 105 (80,       "Mercury",       "Hg",  ( 0.72,  0.72,  0.81, 1.0), 1.49, 1.49, 1.76 ,  1 , 1.27 ,  2 , 1.10 ),
 106 (81,      "Thallium",       "Tl",  ( 0.65,  0.32,  0.30, 1.0), 1.48, 1.48, 2.08 ,  1 , 1.47 ,  3 , 0.95 ),
 107 (82,          "Lead",       "Pb",  ( 0.34,  0.34,  0.38, 1.0), 1.47, 1.47, 1.81 ,  2 , 1.20 ,  4 , 0.84 ),
 108 (83,       "Bismuth",       "Bi",  ( 0.61,  0.30,  0.70, 1.0), 1.46, 1.46, 1.63 ,  1 , 0.98 ,  3 , 0.96 ,  5 , 0.74 ),
 109 (84,      "Polonium",       "Po",  ( 0.67,  0.36,   0.0, 1.0), 1.46, 1.46, 1.53 ,  6 , 0.67 ),
 110 (85,      "Astatine",       "At",  ( 0.45,  0.30,  0.27, 1.0), 1.45, 1.45, 1.43 , -3 , 2.22 ,  3 , 0.85 ,  5 , 0.46 ),
 111 (86,         "Radon",       "Rn",  ( 0.25,  0.50,  0.58, 1.0), 1.00, 1.00, 1.34 ),
 112 (87,      "Francium",       "Fr",  ( 0.25,   0.0,   0.4, 1.0), 1.00, 1.00, 1.00 ,  1 , 1.80 ),
 113 (88,        "Radium",       "Ra",  (  0.0,  0.49,   0.0, 1.0), 1.00, 1.00, 1.00 ,  2 , 1.43 ),
 114 (89,      "Actinium",       "Ac",  ( 0.43,  0.67,  0.98, 1.0), 1.00, 1.00, 1.00 ,  3 , 1.18 ),
 115 (90,       "Thorium",       "Th",  (  0.0,  0.72,   1.0, 1.0), 1.65, 1.65, 1.00 ,  4 , 1.02 ),
 116 (91,  "Protactinium",       "Pa",  (  0.0,  0.63,   1.0, 1.0), 1.00, 1.00, 1.00 ,  3 , 1.13 ,  4 , 0.98 ,  5 , 0.89 ),
 117 (92,       "Uranium",        "U",  (  0.0,  0.56,   1.0, 1.0), 1.42, 1.42, 1.00 ,  4 , 0.97 ,  6 , 0.80 ),
 118 (93,     "Neptunium",       "Np",  (  0.0,  0.50,   1.0, 1.0), 1.00, 1.00, 1.00 ,  3 , 1.10 ,  4 , 0.95 ,  7 , 0.71 ),
 119 (94,     "Plutonium",       "Pu",  (  0.0,  0.41,   1.0, 1.0), 1.00, 1.00, 1.00 ,  3 , 1.08 ,  4 , 0.93 ),
 120 (95,     "Americium",       "Am",  ( 0.32,  0.36,  0.94, 1.0), 1.00, 1.00, 1.00 ,  3 , 1.07 ,  4 , 0.92 ),
 121 (96,        "Curium",       "Cm",  ( 0.47,  0.36,  0.89, 1.0), 1.00, 1.00, 1.00 ),
 122 (97,     "Berkelium",       "Bk",  ( 0.54,  0.30,  0.89, 1.0), 1.00, 1.00, 1.00 ),
 123 (98,   "Californium",       "Cf",  ( 0.63,  0.21,  0.83, 1.0), 1.00, 1.00, 1.00 ),
 124 (99,   "Einsteinium",       "Es",  ( 0.70,  0.12,  0.83, 1.0), 1.00, 1.00, 1.00 ),
 125 (100,       "Fermium",       "Fm", ( 0.70,  0.12,  0.72, 1.0), 1.00, 1.00, 1.00 ),
 126 (101,   "Mendelevium",       "Md", ( 0.70,  0.05,  0.65, 1.0), 1.00, 1.00, 1.00 ),
 127 (102,      "Nobelium",       "No", ( 0.74,  0.05,  0.52, 1.0), 1.00, 1.00, 1.00 ),
 128 (103,    "Lawrencium",       "Lr", ( 0.78,   0.0,   0.4, 1.0), 1.00, 1.00, 1.00 ),
 129 (104,       "Vacancy",      "Vac", (  0.5,   0.5,   0.5, 1.0), 1.00, 1.00, 1.00),
 130 (105,       "Default",  "Default", (  1.0,   1.0,   1.0, 1.0), 1.00, 1.00, 1.00),
 131 (106,         "Stick",    "Stick", (  0.5,   0.5,   0.5, 1.0), 1.00, 1.00, 1.00),
 132 )
 133
 134 # This list here contains all data of the elements and will be used during
 135 # runtime. It is a list of classes.
 136 # During executing Atomic Blender, the list will be initialized with the fixed
 137 # data from above via the class structure below (ElementProp). We
 138 # have then one fixed list (above), which will never be changed, and a list of
 139 # classes with same data. The latter can be modified via loading a separate
 140 # custom data file.
 141 ELEMENTS = []
 142
 143 # This is the class, which stores the properties for one element.
 144 class ElementProp(object):
 145     __slots__ = ('number', 'name', 'short_name', 'color', 'radii', 'radii_ionic')
 146     def __init__(self, number, name, short_name, color, radii, radii_ionic):
 147         self.number = number
 148         self.name = name
 149         self.short_name = short_name
 150         self.color = color
 151         self.radii = radii
 152         self.radii_ionic = radii_ionic
 153
 154 # This is the class, which stores the properties of one atom.
 155 class AtomProp(object):
 156     __slots__ = ('element', 'name', 'location', 'radius', 'color', 'material')
 157     def __init__(self, element, name, location, radius, color, material):
 158         self.element = element
 159         self.name = name
 160         self.location = location
 161         self.radius = radius
 162         self.color = color
 163         self.material = material
 164
 165 # This is the class, which stores the two atoms of one stick.
 166 class StickProp(object):
 167     __slots__ = ('atom1', 'atom2', 'number', 'dist')
 168     def __init__(self, atom1, atom2, number, dist):
 169         self.atom1 = atom1
 170         self.atom2 = atom2
 171         self.number = number
 172         self.dist = dist
 173
 174 # -----------------------------------------------------------------------------
 175 #                                                           Some basic routines
 176
 177
 178 # The function, which reads all necessary properties of the elements.
 179 def read_elements():
 180
 181     del ELEMENTS[:]
 182
 183     for item in ELEMENTS_DEFAULT:
 184
 185         # All three radii into a list
 186         radii = [item[4],item[5],item[6]]
 187         # The handling of the ionic radii will be done later. So far, it is an
 188         # empty list.
 189         radii_ionic = []
 190
 191         li = ElementProp(item[0],item[1],item[2],item[3],
 192                                      radii,radii_ionic)
 193         ELEMENTS.append(li)
 194
 195
 196 # The function, which reads the x,y,z positions of all atoms in a PDB
 197 # file.
 198 #
 199 # filepath_pdb: path to pdb file
 200 # radiustype  : '0' default
 201 #               '1' atomic radii
 202 #               '2' van der Waals
 203 def read_pdb_file(filepath_pdb, radiustype):
 204
 205     # The list of all atoms as read from the PDB file.
 206     all_atoms  = []
 207
 208     # Open the pdb file ...
 209     filepath_pdb_p = open(filepath_pdb, "r")
 210
 211     #Go to the line, in which "ATOM" or "HETATM" appears.
 212     for line in filepath_pdb_p:
 213         split_list = line.split(' ')
 214         if "ATOM" in split_list[0]:
 215             break
 216         if "HETATM" in split_list[0]:
 217             break
 218
 219     j = 0
 220     # This is in fact an endless 'while loop', ...
 221     while j > -1:
 222
 223         # ... the loop is broken here (EOF) ...
 224         if line == "":
 225             break
 226
 227         # If there is a "TER" we need to put empty entries into the lists
 228         # in order to not destroy the order of atom numbers and same numbers
 229         # used for sticks. "TER? What is that?" TER indicates the end of a
 230         # list of ATOM/HETATM records for a chain.
 231         if "TER" in line:
 232             short_name = "TER"
 233             name = "TER"
 234             radius = 0.0
 235             # 2019-03-14, New
 236             color = [0,0,0, 0]
 237             location = Vector((0,0,0))
 238             # Append the TER into the list. Material remains empty so far.
 239             all_atoms.append(AtomProp(short_name,
 240                                       name,
 241                                       location,
 242                                       radius,
 243                                       color,[]))
 244
 245         # If 'ATOM or 'HETATM' appears in the line then do ...
 246         elif "ATOM" in line or "HETATM" in line:
 247
 248             # What follows is due to deviations which appear from PDB to
 249             # PDB file. It is very special!
 250             #
 251             # PLEASE, DO NOT CHANGE! ............................... from here
 252             if line[12:13] == " " or line[12:13].isdigit() == True:
 253                 short_name = line[13:14]
 254                 if line[14:15].islower() == True:
 255                     short_name = short_name + line[14:15]
 256             elif line[12:13].isupper() == True:
 257                 short_name = line[12:13]
 258                 if line[13:14].isalpha() == True:
 259                     short_name = short_name + line[13:14]
 260             else:
 261                 print("Atomic Blender: Strange error in PDB file.\n"
 262                       "Look for element names at positions 13-16 and 78-79.\n")
 263                 return -1
 264
 265             if len(line) >= 78:
 266
 267                 if line[76:77] == " ":
 268                     short_name2 = line[76:77]
 269                 else:
 270                     short_name2 = line[76:78]
 271
 272                 if short_name2.isalpha() == True:
 273                     FOUND = False
 274                     for element in ELEMENTS:
 275                         if str.upper(short_name2) == str.upper(element.short_name):
 276                             FOUND = True
 277                             break
 278                     if FOUND == False:
 279                         short_name = short_name2
 280
 281             # ....................................................... to here.
 282
 283             # Go through all elements and find the element of the current atom.
 284             FLAG_FOUND = False
 285             for element in ELEMENTS:
 286                 if str.upper(short_name) == str.upper(element.short_name):
 287                     # Give the atom its proper names, color and radius:
 288                     short_name = str.upper(element.short_name)
 289                     name = element.name
 290                     # int(radiustype) => type of radius:
 291                     # pre-defined (0), atomic (1) or van der Waals (2)
 292                     radius = float(element.radii[int(radiustype)])
 293                     color = element.color
 294                     FLAG_FOUND = True
 295                     break
 296
 297             # Is it a vacancy or an 'unknown atom' ?
 298             if FLAG_FOUND == False:
 299                 # Give this atom also a name. If it is an 'X' then it is a
 300                 # vacancy. Otherwise ...
 301                 if "X" in short_name:
 302                     short_name = "VAC"
 303                     name = "Vacancy"
 304                     radius = float(ELEMENTS[-3].radii[int(radiustype)])
 305                     color = ELEMENTS[-3].color
 306                 # ... take what is written in the PDB file. These are somewhat
 307                 # unknown atoms. This should never happen, the element list is
 308                 # almost complete. However, we do this due to security reasons.
 309                 else:
 310                     short_name = str.upper(short_name)
 311                     name = str.upper(short_name)
 312                     radius = float(ELEMENTS[-2].radii[int(radiustype)])
 313                     color = ELEMENTS[-2].color
 314
 315             # x,y and z are at fixed positions in the PDB file.
 316             x = float(line[30:38].rsplit()[0])
 317             y = float(line[38:46].rsplit()[0])
 318             z = float(line[46:55].rsplit()[0])
 319
 320             location = Vector((x,y,z))
 321
 322             j += 1
 323
 324             # Append the atom to the list. Material remains empty so far.
 325             all_atoms.append(AtomProp(short_name,
 326                                       name,
 327                                       location,
 328                                       radius,
 329                                       color,[]))
 330
 331         line = filepath_pdb_p.readline()
 332         line = line[:-1]
 333
 334     filepath_pdb_p.close()
 335     # From above it can be clearly seen that j is now the number of all atoms.
 336     Number_of_total_atoms = j
 337
 338     return (Number_of_total_atoms, all_atoms)
 339
 340
 341 # The function, which reads the sticks in a PDB file.
 342 def read_pdb_file_sticks(filepath_pdb, use_sticks_bonds, all_atoms):
 343
 344     # The list of all sticks.
 345     all_sticks = []
 346
 347     # Open the PDB file.
 348     filepath_pdb_p = open(filepath_pdb, "r")
 349
 350     line = filepath_pdb_p.readline()
 351     split_list = line.split(' ')
 352
 353     # Go to the first entry
 354     # DO NOT CHANGE 'CONECT', read below.
 355     if "CONECT" not in split_list[0]:
 356         for line in filepath_pdb_p:
 357             split_list = line.split(' ')
 358             if "CONECT" in split_list[0]:
 359                 break
 360
 361     Number_of_sticks = 0
 362     sticks_double = 0
 363     j = 0
 364     # This is in fact an endless while loop, ...
 365     while j > -1:
 366
 367         # ... which is broken here (EOF) ...
 368         if line == "":
 369             break
 370         # ... or here, when no 'CONECT' appears anymore.
 371         if "CONECT" not in line:
 372             break
 373
 374         # Note 2019-03-16: in a PDB file the identifier for sticks is called
 375         # 'CONECT' and NOT 'CONNECT'! Please leave this as is, otherwise the
 376         # sticks are NOT correctly imported.
 377
 378         # The strings of the atom numbers do have a clear position in the file
 379         # (From 7 to 12, from 13 to 18 and so on.) and one needs to consider
 380         # this. One could also use the split function but then one gets into
 381         # trouble if there are lots of atoms: For instance, it may happen that
 382         # one has
 383         #                   CONECT 11111  22244444
 384         #
 385         # In Fact it means that atom No. 11111 has a connection with atom
 386         # No. 222 but also with atom No. 44444. The split function would give
 387         # me only two numbers (11111 and 22244444), which is wrong.
 388
 389         # Cut spaces from the right and 'CONECT' at the beginning
 390         line = line.rstrip()
 391         line = line[6:]
 392         # Amount of loops
 393         length = len(line)
 394         loops  = int(length/5)
 395
 396         # List of atoms
 397         atom_list = []
 398         for i in range(loops):
 399             number = line[5*i:5*(i+1)].rsplit()
 400             if number != []:
 401                 if number[0].isdigit() == True:
 402                     atom_number = int(number[0])
 403                     atom_list.append(atom_number)
 404
 405         # The first atom is connected with all the others in the list.
 406         atom1 = atom_list[0]
 407
 408         # For all the other atoms in the list do:
 409         for atom2 in atom_list[1:]:
 410
 411             if use_sticks_bonds == True:
 412                 number = atom_list[1:].count(atom2)
 413
 414                 if number == 2 or number == 3:
 415                     basis_list = list(set(atom_list[1:]))
 416
 417                     if len(basis_list) > 1:
 418                         basis1 = (all_atoms[atom1-1].location
 419                                 - all_atoms[basis_list[0]-1].location)
 420                         basis2 = (all_atoms[atom1-1].location
 421                                 - all_atoms[basis_list[1]-1].location)
 422                         plane_n = basis1.cross(basis2)
 423
 424                         dist_n = (all_atoms[atom1-1].location
 425                                 - all_atoms[atom2-1].location)
 426                         dist_n = dist_n.cross(plane_n)
 427                         dist_n = dist_n / dist_n.length
 428                     else:
 429                         dist_n = (all_atoms[atom1-1].location
 430                                 - all_atoms[atom2-1].location)
 431                         dist_n = Vector((dist_n[1],-dist_n[0],0))
 432                         dist_n = dist_n / dist_n.length
 433
 434                 elif number > 3:
 435                     number = 1
 436                     dist_n = None
 437                 else:
 438                     dist_n = None
 439             else:
 440                 number = 1
 441                 dist_n = None
 442
 443             # Note that in a PDB file, sticks of one atom pair can appear a
 444             # couple of times. (Only god knows why ...)
 445             # So, does a stick between the considered atoms already exist?
 446             FLAG_BAR = False
 447             for k in range(Number_of_sticks):
 448                 if ((all_sticks[k].atom1 == atom1 and all_sticks[k].atom2 == atom2) or
 449                     (all_sticks[k].atom2 == atom1 and all_sticks[k].atom1 == atom2)):
 450                     sticks_double += 1
 451                     # If yes, then FLAG on 'True'.
 452                     FLAG_BAR       = True
 453                     break
 454
 455             # If the stick is not yet registered (FLAG_BAR == False), then
 456             # register it!
 457             if FLAG_BAR == False:
 458                 all_sticks.append(StickProp(atom1,atom2,number,dist_n))
 459                 Number_of_sticks += 1
 460                 j += 1
 461
 462         line = filepath_pdb_p.readline()
 463         line = line.rstrip()
 464
 465     filepath_pdb_p.close()
 466
 467     return all_sticks
 468
 469
 470 # Function, which produces a cylinder. All is somewhat easy to understand.
 471 def build_stick(radius, length, sectors, element_name):
 472
 473     dphi = 2.0 * pi/(float(sectors)-1)
 474
 475     # Vertices
 476     vertices_top    = [Vector((0,0,length / 2.0))]
 477     vertices_bottom = [Vector((0,0,-length / 2.0))]
 478     vertices = []
 479     for i in range(sectors-1):
 480         x = radius * cos( dphi * i )
 481         y = radius * sin( dphi * i )
 482         z =  length / 2.0
 483         vertex = Vector((x,y,z))
 484         vertices_top.append(vertex)
 485         z = -length / 2.0
 486         vertex = Vector((x,y,z))
 487         vertices_bottom.append(vertex)
 488     vertices = vertices_top + vertices_bottom
 489
 490     # Side facets (Cylinder)
 491     faces1 = []
 492     for i in range(sectors-1):
 493         if i == sectors-2:
 494             faces1.append(  [i+1, 1, 1+sectors, i+1+sectors] )
 495         else:
 496             faces1.append(  [i+1, i+2, i+2+sectors, i+1+sectors] )
 497
 498     # Top facets
 499     faces2 = []
 500     for i in range(sectors-1):
 501         if i == sectors-2:
 502             face_top = [0,sectors-1,1]
 503             face_bottom = [sectors,2*sectors-1,sectors+1]
 504         else:
 505             face_top    = [0]
 506             face_bottom = [sectors]
 507             for j in range(2):
 508                 face_top.append(i+j+1)
 509                 face_bottom.append(i+j+1+sectors)
 510         faces2.append(face_top)
 511         faces2.append(face_bottom)
 512
 513     # Build the mesh, Cylinder
 514     cylinder = bpy.data.meshes.new(element_name+"_sticks_cylinder")
 515     cylinder.from_pydata(vertices, [], faces1)
 516     cylinder.update()
 517     new_cylinder = bpy.data.objects.new(element_name+"_sticks_cylinder", cylinder)
 518     # Attention: the linking will be done a few moments later, after this
 519     # is done definition.
 520
 521     # Build the mesh, Cups
 522     cups = bpy.data.meshes.new(element_name+"_sticks_cup")
 523     cups.from_pydata(vertices, [], faces2)
 524     cups.update()
 525     new_cups = bpy.data.objects.new(element_name+"_sticks_cup", cups)
 526     # Attention: the linking will be done a few moments later, after this
 527     # is done definition.
 528
 529     return new_cylinder, new_cups
 530
 531
 532 # Rotate an object.
 533 def rotate_object(rot_mat, obj):
 534
 535     bpy.ops.object.select_all(action='DESELECT')
 536     obj.select_set(True)
 537
 538     # Decompose world_matrix's components, and from them assemble 4x4 matrices.
 539     orig_loc, orig_rot, orig_scale = obj.matrix_world.decompose()
 540
 541     orig_loc_mat   = Matrix.Translation(orig_loc)
 542     orig_rot_mat   = orig_rot.to_matrix().to_4x4()
 543     orig_scale_mat = (Matrix.Scale(orig_scale[0],4,(1,0,0)) @
 544                       Matrix.Scale(orig_scale[1],4,(0,1,0)) @
 545                       Matrix.Scale(orig_scale[2],4,(0,0,1)))
 546
 547     # Assemble the new matrix.
 548     obj.matrix_world = orig_loc_mat @ rot_mat @ orig_rot_mat @ orig_scale_mat
 549
 550
 551 # Function, which puts a camera and light source into the 3D scene
 552 def camera_light_source(use_camera,
 553                         use_light,
 554                         object_center_vec,
 555                         object_size):
 556
 557     camera_factor = 15.0
 558
 559     # If chosen a camera is put into the scene.
 560     if use_camera == True:
 561
 562         # Assume that the object is put into the global origin. Then, the
 563         # camera is moved in x and z direction, not in y. The object has its
 564         # size at distance sqrt(object_size) from the origin. So, move the
 565         # camera by this distance times a factor of camera_factor in x and z.
 566         # Then add x, y and z of the origin of the object.
 567         object_camera_vec = Vector((sqrt(object_size) * camera_factor,
 568                                     0.0,
 569                                     sqrt(object_size) * camera_factor))
 570         camera_xyz_vec = object_center_vec + object_camera_vec
 571
 572         # Create the camera
 573         camera_data = bpy.data.cameras.new("A_camera")
 574         camera_data.lens = 45
 575         camera_data.clip_end = 500.0
 576         camera = bpy.data.objects.new("A_camera", camera_data)
 577         camera.location = camera_xyz_vec
 578         bpy.context.collection.objects.link(camera)
 579
 580         # Here the camera is rotated such it looks towards the center of
 581         # the object. The [0.0, 0.0, 1.0] vector along the z axis
 582         z_axis_vec             = Vector((0.0, 0.0, 1.0))
 583         # The angle between the last two vectors
 584         angle                  = object_camera_vec.angle(z_axis_vec, 0)
 585         # The cross-product of z_axis_vec and object_camera_vec
 586         axis_vec               = z_axis_vec.cross(object_camera_vec)
 587         # Rotate 'axis_vec' by 'angle' and convert this to euler parameters.
 588         # 4 is the size of the matrix.
 589         camera.rotation_euler  = Matrix.Rotation(angle, 4, axis_vec).to_euler()
 590
 591         # Rotate the camera around its axis by 90° such that we have a nice
 592         # camera position and view onto the object.
 593         bpy.ops.object.select_all(action='DESELECT')
 594         camera.select_set(True)
 595
 596         # Rotate the camera around its axis 'object_camera_vec' by 90° such
 597         # that we have a nice camera view onto the object.
 598         matrix_rotation = Matrix.Rotation(90/360*2*pi, 4, object_camera_vec)
 599         rotate_object(matrix_rotation, camera)
 600
 601     # Here a lamp is put into the scene, if chosen.
 602     if use_light == True:
 603
 604         # This is the distance from the object measured in terms of %
 605         # of the camera distance. It is set onto 50% (1/2) distance.
 606         lamp_dl = sqrt(object_size) * 15 * 0.5
 607         # This is a factor to which extend the lamp shall go to the right
 608         # (from the camera  point of view).
 609         lamp_dy_right = lamp_dl * (3.0/4.0)
 610
 611         # Create x, y and z for the lamp.
 612         object_lamp_vec = Vector((lamp_dl,lamp_dy_right,lamp_dl))
 613         lamp_xyz_vec = object_center_vec + object_lamp_vec
 614         length = lamp_xyz_vec.length
 615
 616         # As a lamp we use a point source.
 617         lamp_data = bpy.data.lights.new(name="A_lamp", type="POINT")
 618         # We now determine the emission strength of the lamp. Note that the
 619         # intensity depends on 1/r^2. For this we use a value of 100000.0 at a
 620         # distance of 58. This value was determined manually inside Blender.
 621         lamp_data.energy = 500000.0 * ( (length * length) / (58.0 * 58.0) )
 622         lamp = bpy.data.objects.new("A_lamp", lamp_data)
 623         lamp.location = lamp_xyz_vec
 624         bpy.context.collection.objects.link(lamp)
 625
 626         # Some settings for the World: a bit ambient occlusion
 627         bpy.context.scene.world.light_settings.use_ambient_occlusion = True
 628         bpy.context.scene.world.light_settings.ao_factor = 0.1
 629
 630
 631
 632 # Function, which draws the atoms of one type (balls). This is one
 633 # dupliverts structure then.
 634 # Return: the dupliverts structure
 635 def draw_atoms_one_type(draw_all_atoms_type,
 636                         Ball_type,
 637                         Ball_azimuth,
 638                         Ball_zenith,
 639                         Ball_radius_factor,
 640                         object_center_vec,
 641                         collection_molecule):
 642
 643     # Create the vertices composed of the coordinates of all atoms of one type
 644     atom_vertices = []
 645     for atom in draw_all_atoms_type:
 646         # In fact, the object is created in the World's origin.
 647         # This is why 'object_center_vec' is subtracted. At the end
 648         # the whole object is translated back to 'object_center_vec'.
 649         atom_vertices.append(atom[2] - object_center_vec)
 650
 651     # IMPORTANT: First, we create a collection of the element, which contains
 652     # the atoms (balls + mesh) AND the sticks! The definition dealing with the
 653     # sticks will put the sticks inside this collection later on.
 654     coll_element_name = atom[0] # the element name
 655     # Create the new collection and ...
 656     coll_element = bpy.data.collections.new(coll_element_name)
 657     # ... link it to the collection, which contains all parts of the
 658     # molecule.
 659     collection_molecule.children.link(coll_element)
 660
 661     # Now, create a collection for the atoms, which includes the representative
 662     # ball and the mesh.
 663     coll_atom_name = atom[0] + "_atom"
 664     # Create the new collection and ...
 665     coll_atom = bpy.data.collections.new(coll_atom_name)
 666     # ... link it to the collection, which contains all parts of the
 667     # element (ball and mesh).
 668     coll_element.children.link(coll_atom)
 669
 670     # Build the mesh
 671     atom_mesh = bpy.data.meshes.new("Mesh_"+atom[0])
 672     atom_mesh.from_pydata(atom_vertices, [], [])
 673     atom_mesh.update()
 674     new_atom_mesh = bpy.data.objects.new(atom[0] + "_mesh", atom_mesh)
 675
 676     # Link active object to the new collection
 677     coll_atom.objects.link(new_atom_mesh)
 678
 679     # Now, build a representative sphere (atom).
 680     if atom[0] == "Vacancy":
 681         bpy.ops.mesh.primitive_cube_add(
 682                         align='WORLD', enter_editmode=False,
 683                         location=(0.0, 0.0, 0.0),
 684                         rotation=(0.0, 0.0, 0.0))
 685     else:
 686         # NURBS balls
 687         if Ball_type == "0":
 688             bpy.ops.surface.primitive_nurbs_surface_sphere_add(
 689                         align='WORLD', enter_editmode=False,
 690                         location=(0,0,0), rotation=(0.0, 0.0, 0.0))
 691         # UV balls
 692         elif Ball_type == "1":
 693             bpy.ops.mesh.primitive_uv_sphere_add(
 694                         segments=Ball_azimuth, ring_count=Ball_zenith,
 695                         align='WORLD', enter_editmode=False,
 696                         location=(0,0,0), rotation=(0, 0, 0))
 697         # Meta balls
 698         elif Ball_type == "2":
 699             bpy.ops.object.metaball_add(type='BALL', align='WORLD',
 700                         enter_editmode=False, location=(0, 0, 0),
 701                         rotation=(0, 0, 0))
 702
 703     ball = bpy.context.view_layer.objects.active
 704     # Hide this ball because its appearance has no meaning. It is just the
 705     # representative ball. The ball is visible at the vertices of the mesh.
 706     # Rememmber, this is a dupliverts construct!
 707     # However, hiding does not work with meta balls!
 708     if Ball_type == "0" or Ball_type == "1":
 709         ball.hide_set(True)
 710     # Scale up/down the ball radius.
 711     ball.scale  = (atom[3]*Ball_radius_factor,) * 3
 712
 713     if atom[0] == "Vacancy":
 714         ball.name = atom[0] + "_cube"
 715     else:
 716         ball.name = atom[0] + "_ball"
 717
 718     ball.active_material = atom[1]
 719     ball.parent = new_atom_mesh
 720     new_atom_mesh.instance_type = 'VERTS'
 721     # The object is back translated to 'object_center_vec'.
 722     new_atom_mesh.location = object_center_vec
 723
 724     # Note the collection where the ball was placed into.
 725     coll_all = ball.users_collection
 726     if len(coll_all) > 0:
 727         coll_past = coll_all[0]
 728     else:
 729         coll_past = bpy.context.scene.collection
 730
 731     # Put the atom into the new collection 'atom' and ...
 732     coll_atom.objects.link(ball)
 733     # ... unlink the atom from the other collection.
 734     coll_past.objects.unlink(ball)
 735
 736     return new_atom_mesh, coll_element
 737
 738
 739 # Function, which draws the sticks with help of the dupliverts technique.
 740 # Return: list of dupliverts structures.
 741 def draw_sticks_dupliverts(all_atoms,
 742                            atom_all_types_list,
 743                            center,
 744                            all_sticks,
 745                            Stick_diameter,
 746                            Stick_sectors,
 747                            Stick_unit,
 748                            Stick_dist,
 749                            use_sticks_smooth,
 750                            use_sticks_color,
 751                            list_coll_elements):
 752
 753     dl = Stick_unit
 754
 755     if use_sticks_color == False:
 756         stick_material = bpy.data.materials.new(ELEMENTS[-1].name)
 757         stick_material.use_nodes = True
 758         mat_P_BSDF = next(n for n in stick_material.node_tree.nodes
 759                           if n.type == "BSDF_PRINCIPLED")
 760         mat_P_BSDF.inputs['Base Color'].default_value = ELEMENTS[-1].color
 761
 762     # Sort the sticks and put them into a new list such that ...
 763     sticks_all_lists = []
 764     if use_sticks_color == True:
 765         for atom_type in atom_all_types_list:
 766             if atom_type[0] == "TER":
 767                 continue
 768             sticks_list = []
 769             for stick in all_sticks:
 770                 for repeat in range(stick.number):
 771
 772                     atom1 = copy(all_atoms[stick.atom1-1].location)-center
 773                     atom2 = copy(all_atoms[stick.atom2-1].location)-center
 774
 775                     dist =  Stick_diameter * Stick_dist
 776
 777                     if stick.number == 2:
 778                         if repeat == 0:
 779                             atom1 += (stick.dist * dist)
 780                             atom2 += (stick.dist * dist)
 781                         if repeat == 1:
 782                             atom1 -= (stick.dist * dist)
 783                             atom2 -= (stick.dist * dist)
 784
 785                     if stick.number == 3:
 786                         if repeat == 0:
 787                             atom1 += (stick.dist * dist)
 788                             atom2 += (stick.dist * dist)
 789                         if repeat == 2:
 790                             atom1 -= (stick.dist * dist)
 791                             atom2 -= (stick.dist * dist)
 792
 793                     dv = atom1 - atom2
 794                     n  = dv / dv.length
 795                     if atom_type[0] == all_atoms[stick.atom1-1].name:
 796                         location = atom1
 797                         name     = "_" + all_atoms[stick.atom1-1].name
 798                         material = all_atoms[stick.atom1-1].material
 799                         sticks_list.append([name, location, dv, material])
 800                     if atom_type[0] == all_atoms[stick.atom2-1].name:
 801                         location = atom1 - n * dl * int(ceil(dv.length / (2.0 * dl)))
 802                         name     = "_" + all_atoms[stick.atom2-1].name
 803                         material = all_atoms[stick.atom2-1].material
 804                         sticks_list.append([name, location, dv, material])
 805
 806             if sticks_list != []:
 807                 sticks_all_lists.append(sticks_list)
 808     else:
 809         sticks_list = []
 810         for stick in all_sticks:
 811
 812             if stick.number > 3:
 813                 stick.number = 1
 814
 815             for repeat in range(stick.number):
 816
 817                 atom1 = copy(all_atoms[stick.atom1-1].location)-center
 818                 atom2 = copy(all_atoms[stick.atom2-1].location)-center
 819
 820                 dist =  Stick_diameter * Stick_dist
 821
 822                 if stick.number == 2:
 823                     if repeat == 0:
 824                         atom1 += (stick.dist * dist)
 825                         atom2 += (stick.dist * dist)
 826                     if repeat == 1:
 827                         atom1 -= (stick.dist * dist)
 828                         atom2 -= (stick.dist * dist)
 829                 if stick.number == 3:
 830                     if repeat == 0:
 831                         atom1 += (stick.dist * dist)
 832                         atom2 += (stick.dist * dist)
 833                     if repeat == 2:
 834                         atom1 -= (stick.dist * dist)
 835                         atom2 -= (stick.dist * dist)
 836
 837                 dv = atom1 - atom2
 838                 n  = dv / dv.length
 839                 location = atom1
 840                 material = stick_material
 841                 sticks_list.append(["", location, dv, material])
 842
 843         sticks_all_lists.append(sticks_list)
 844
 845     atom_object_list = []
 846     # ... the sticks in the list can be drawn:
 847     for stick_list in sticks_all_lists:
 848         vertices = []
 849         faces    = []
 850         i = 0
 851
 852         # What follows is school mathematics! :-) We construct equidistant
 853         # planes, on which the stcik sections (cylinders) are perpendicular on.
 854         for stick in stick_list:
 855
 856             dv = stick[2]
 857             v1 = stick[1]
 858             n  = dv / dv.length
 859             gamma = -n.dot(v1)
 860             b     = v1 + gamma * n
 861             n_b   = b / b.length
 862
 863             if use_sticks_color == True:
 864                 loops = int(ceil(dv.length / (2.0 * dl)))
 865             else:
 866                 loops = int(ceil(dv.length / dl))
 867
 868             for j in range(loops):
 869
 870                 # The plane, which is normal to the length of the cylinder,
 871                 # will have a 1/100 of the stick diameter. => When decreasing
 872                 # the size of the stick diameter, the plane will not be visible.
 873                 f = 0.01
 874                 g  = v1 - n * dl / 2.0 - n * dl * j
 875                 p1 = g + n_b * Stick_diameter * f
 876                 p2 = g - n_b * Stick_diameter * f
 877                 p3 = g - n_b.cross(n) * Stick_diameter * f
 878                 p4 = g + n_b.cross(n) * Stick_diameter * f
 879
 880                 vertices.append(p1)
 881                 vertices.append(p2)
 882                 vertices.append(p3)
 883                 vertices.append(p4)
 884                 faces.append((i*4+0,i*4+2,i*4+1,i*4+3))
 885                 i += 1
 886
 887         # Create a collection for the sticks, which includes the representative
 888         # cylinders, cups and the mesh.
 889         coll_name = stick[0][1:] + "_sticks"
 890         # Create the collection and ...
 891         coll = bpy.data.collections.new(coll_name)
 892         # ... link it to the collection, which contains all parts of the
 893         # element. 'stick[0][1:]' contains the name of the element!
 894         for coll_element_from_list in list_coll_elements:
 895             if stick[0][1:] in coll_element_from_list.name:
 896                 break
 897         coll_element_from_list.children.link(coll)
 898
 899         # Build the mesh.
 900         mesh = bpy.data.meshes.new("Sticks_"+stick[0][1:])
 901         mesh.from_pydata(vertices, [], faces)
 902         mesh.update()
 903         new_mesh = bpy.data.objects.new(stick[0][1:]+"_sticks_mesh", mesh)
 904         # Link active object to the new collection
 905         coll.objects.link(new_mesh)
 906
 907         # Build the object. Get the cylinder from the 'build_stick' function.
 908         stick_cylinder, stick_cups = build_stick(Stick_diameter,
 909                                                  dl,
 910                                                  Stick_sectors,
 911                                                  stick[0][1:])
 912         # Link active object to the new collection.
 913         coll.objects.link(stick_cylinder)
 914         coll.objects.link(stick_cups)
 915
 916         # Assign the material.
 917         stick_cylinder.active_material = stick[3]
 918         stick_cups.active_material = stick[3]
 919
 920         # Smooth the cylinders.
 921         if use_sticks_smooth == True:
 922             bpy.ops.object.select_all(action='DESELECT')
 923             stick_cylinder.select_set(True)
 924             stick_cups.select_set(True)
 925             bpy.ops.object.shade_smooth()
 926
 927         # Hide these objects because their appearance has no meaning. They are
 928         # just the representative objects. The cylinder and cups are visible at
 929         # the vertices of the mesh. Rememmber, this is a dupliverts construct!
 930         stick_cylinder.hide_set(True)
 931         stick_cups.hide_set(True)
 932
 933         # Parenting the mesh to the cylinder.
 934         stick_cylinder.parent = new_mesh
 935         stick_cups.parent = new_mesh
 936         new_mesh.instance_type = 'FACES'
 937         new_mesh.location = center
 938         atom_object_list.append(new_mesh)
 939
 940     # Return the list of dupliverts structures.
 941     return atom_object_list
 942
 943
 944 # Function, which draws the sticks with help of the skin and subdivision
 945 # modifiers.
 946 def draw_sticks_skin(all_atoms,
 947                      all_sticks,
 948                      Stick_diameter,
 949                      use_sticks_smooth,
 950                      sticks_subdiv_view,
 951                      sticks_subdiv_render,
 952                      coll_molecule):
 953
 954     # These counters are for the edges, in the shape [i,i+1].
 955     i = 0
 956
 957     # This is the list of vertices, containing the atom position
 958     # (vectors)).
 959     stick_vertices = []
 960     # This is the 'same' list, which contains not vector position of
 961     # the atoms but their numbers. It is used to handle the edges.
 962     stick_vertices_nr = []
 963     # This is the list of edges.
 964     stick_edges = []
 965
 966     # Go through the list of all sticks. For each stick do:
 967     for stick in all_sticks:
 968
 969         # Each stick has two atoms = two vertices.
 970
 971         """
 972         [ 0,1 ,  3,4 ,  0,8 ,  7,3]
 973         [[0,1], [2,3], [4,5], [6,7]]
 974
 975         [ 0,1 ,  3,4 ,  x,8 ,   7,x]    x:deleted
 976         [[0,1], [2,3], [0,5], [6,2]]
 977         """
 978
 979         # Check, if the vertex (atom) is already in the vertex list.
 980         # edge: [s1,s2]
 981         FLAG_s1 = False
 982         s1 = 0
 983         for stick2 in stick_vertices_nr:
 984             if stick2 == stick.atom1-1:
 985                 FLAG_s1 = True
 986                 break
 987             s1 += 1
 988         FLAG_s2 = False
 989         s2 = 0
 990         for stick2 in stick_vertices_nr:
 991             if stick2 == stick.atom2-1:
 992                 FLAG_s2 = True
 993                 break
 994             s2 += 1
 995
 996         # If the vertex (atom) is not yet in the vertex list:
 997         # append the number of atom and the vertex to the two lists.
 998         # For the first atom:
 999         if FLAG_s1 == False:
1000             atom1 = copy(all_atoms[stick.atom1-1].location)
1001             stick_vertices.append(atom1)
1002             stick_vertices_nr.append(stick.atom1-1)
1003         # For the second atom:
1004         if FLAG_s2 == False:
1005             atom2 = copy(all_atoms[stick.atom2-1].location)
1006             stick_vertices.append(atom2)
1007             stick_vertices_nr.append(stick.atom2-1)
1008
1009         # Build the edges:
1010
1011         # If both vertices (atoms) were not in the lists, then
1012         # the edge is simply [i,i+1]. These are two new vertices
1013         # (atoms), so increase i by 2.
1014         if FLAG_s1 == False and FLAG_s2 == False:
1015             stick_edges.append([i,i+1])
1016             i += 2
1017         # Both vertices (atoms) were already in the list, so then
1018         # use the vertices (atoms), which already exist. They are
1019         # at positions s1 and s2.
1020         if FLAG_s1 == True and FLAG_s2 == True:
1021             stick_edges.append([s1,s2])
1022         # The following two if cases describe the situation that
1023         # only one vertex (atom) was in the list. Since only ONE
1024         # new vertex was added, increase i by one.
1025         if FLAG_s1 == True and FLAG_s2 == False:
1026             stick_edges.append([s1,i])
1027             i += 1
1028         if FLAG_s1 == False and FLAG_s2 == True:
1029             stick_edges.append([i,s2])
1030             i += 1
1031
1032     # Build the mesh of the sticks
1033     stick_mesh = bpy.data.meshes.new("Mesh_sticks")
1034     stick_mesh.from_pydata(stick_vertices, stick_edges, [])
1035     stick_mesh.update()
1036     new_stick_mesh = bpy.data.objects.new("Sticks", stick_mesh)
1037     # Link the active mesh to the molecule collection
1038     coll_molecule.objects.link(new_stick_mesh)
1039
1040     # Apply the skin modifier.
1041     new_stick_mesh.modifiers.new(name="Sticks_skin", type='SKIN')
1042     # Smooth the skin surface if this option has been chosen.
1043     new_stick_mesh.modifiers[0].use_smooth_shade = use_sticks_smooth
1044     # Apply the Subdivision modifier.
1045     new_stick_mesh.modifiers.new(name="Sticks_subsurf", type='SUBSURF')
1046     # Options: choose the levels
1047     new_stick_mesh.modifiers[1].levels = sticks_subdiv_view
1048     new_stick_mesh.modifiers[1].render_levels = sticks_subdiv_render
1049
1050     stick_material = bpy.data.materials.new(ELEMENTS[-1].name)
1051     stick_material.use_nodes = True
1052     mat_P_BSDF = next(n for n in stick_material.node_tree.nodes
1053                       if n.type == "BSDF_PRINCIPLED")
1054     mat_P_BSDF.inputs['Base Color'].default_value = ELEMENTS[-1].color
1055     new_stick_mesh.active_material = stick_material
1056
1057     # This is for putting the radius of the sticks onto
1058     # the desired value 'Stick_diameter'
1059     bpy.context.view_layer.objects.active = new_stick_mesh
1060     # EDIT mode
1061     bpy.ops.object.mode_set(mode='EDIT', toggle=False)
1062     bm = bmesh.from_edit_mesh(new_stick_mesh.data)
1063     bpy.ops.mesh.select_all(action='DESELECT')
1064
1065     # Select all vertices
1066     for v in bm.verts:
1067         v.select = True
1068
1069     # This is somewhat a factor for the radius.
1070     r_f = 4.0
1071     # Apply operator 'skin_resize'.
1072     bpy.ops.transform.skin_resize(
1073         value=(
1074             Stick_diameter * r_f,
1075             Stick_diameter * r_f,
1076             Stick_diameter * r_f,
1077         ),
1078         constraint_axis=(False, False, False),
1079         orient_type='GLOBAL',
1080         mirror=False,
1081         use_proportional_edit=False,
1082         snap=False,
1083         snap_target='CLOSEST',
1084         snap_point=(0, 0, 0),
1085         snap_align=False,
1086         snap_normal=(0, 0, 0),
1087         release_confirm=False,
1088     )
1089     # Back to the OBJECT mode.
1090     bpy.ops.object.mode_set(mode='OBJECT', toggle=False)
1091
1092     return new_stick_mesh
1093
1094
1095 # Draw the sticks the normal way: connect the atoms by simple cylinders.
1096 # Two options: 1. single cylinders parented to an empty
1097 #              2. one single mesh object
1098 def draw_sticks_normal(all_atoms,
1099                        all_sticks,
1100                        center,
1101                        Stick_diameter,
1102                        Stick_sectors,
1103                        use_sticks_smooth,
1104                        use_sticks_one_object,
1105                        use_sticks_one_object_nr,
1106                        coll_molecule):
1107
1108     stick_material = bpy.data.materials.new(ELEMENTS[-1].name)
1109     stick_material.use_nodes = True
1110     mat_P_BSDF = next(n for n in stick_material.node_tree.nodes
1111                       if n.type == "BSDF_PRINCIPLED")
1112     mat_P_BSDF.inputs['Base Color'].default_value = ELEMENTS[-1].color
1113
1114     up_axis = Vector([0.0, 0.0, 1.0])
1115
1116     # For all sticks, do ...
1117     list_group = []
1118     list_group_sub = []
1119     counter = 0
1120     for stick in all_sticks:
1121
1122         # The vectors of the two atoms
1123         atom1 = all_atoms[stick.atom1-1].location-center
1124         atom2 = all_atoms[stick.atom2-1].location-center
1125         # Location
1126         location = (atom1 + atom2) * 0.5
1127         # The difference of both vectors
1128         v = (atom2 - atom1)
1129         # Angle with respect to the z-axis
1130         angle = v.angle(up_axis, 0)
1131         # Cross-product between v and the z-axis vector. It is the
1132         # vector of rotation.
1133         axis = up_axis.cross(v)
1134         # Calculate Euler angles
1135         euler = Matrix.Rotation(angle, 4, axis).to_euler()
1136         # Create stick
1137         stick = bpy.ops.mesh.primitive_cylinder_add(vertices=Stick_sectors,
1138                                                     radius=Stick_diameter,
1139                                                     depth=v.length,
1140                                                     end_fill_type='NGON',
1141                                                     align='WORLD',
1142                                                     enter_editmode=False,
1143                                                     location=location,
1144                                                     rotation=(0, 0, 0))
1145         # Put the stick into the scene ...
1146         stick = bpy.context.view_layer.objects.active
1147         # ... and rotate the stick.
1148         stick.rotation_euler = euler
1149         # ... and name
1150         stick.name = "Stick_Cylinder"
1151         counter += 1
1152
1153         # Smooth the cylinder.
1154         if use_sticks_smooth == True:
1155             bpy.ops.object.select_all(action='DESELECT')
1156             stick.select_set(True)
1157             bpy.ops.object.shade_smooth()
1158
1159         list_group_sub.append(stick)
1160
1161         if use_sticks_one_object == True:
1162             if counter == use_sticks_one_object_nr:
1163                 bpy.ops.object.select_all(action='DESELECT')
1164                 for stick in list_group_sub:
1165                     stick.select_set(True)
1166                 bpy.ops.object.join()
1167                 list_group.append(bpy.context.view_layer.objects.active)
1168                 bpy.ops.object.select_all(action='DESELECT')
1169                 list_group_sub = []
1170                 counter = 0
1171         else:
1172             # Material ...
1173             stick.active_material = stick_material
1174
1175     if use_sticks_one_object == True:
1176         bpy.ops.object.select_all(action='DESELECT')
1177         for stick in list_group_sub:
1178             stick.select_set(True)
1179         bpy.ops.object.join()
1180         list_group.append(bpy.context.view_layer.objects.active)
1181         bpy.ops.object.select_all(action='DESELECT')
1182
1183         for group in list_group:
1184             group.select_set(True)
1185         bpy.ops.object.join()
1186         bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY',
1187                                    center='MEDIAN')
1188         sticks = bpy.context.view_layer.objects.active
1189         sticks.active_material = stick_material
1190
1191         sticks.location += center
1192
1193         # Collections
1194         # ===========
1195         # Note the collection where the sticks were placed into.
1196         coll_all = sticks.users_collection
1197         if len(coll_all) > 0:
1198             coll_past = coll_all[0]
1199         else:
1200             coll_past = bpy.context.scene.collection
1201
1202         # Link the sticks with the collection of the molecule ...
1203         coll_molecule.objects.link(sticks)
1204         # ... and unlink them from the collection it has been before.
1205         coll_past.objects.unlink(sticks)
1206
1207         return sticks
1208     else:
1209         # Here we use an empty ...
1210         bpy.ops.object.empty_add(type='ARROWS',
1211                                   align='WORLD',
1212                                   location=(0, 0, 0),
1213                                   rotation=(0, 0, 0))
1214         sticks_empty = bpy.context.view_layer.objects.active
1215         sticks_empty.name = "A_sticks_empty"
1216         # ... that is parent to all sticks. With this, we can better move
1217         # all sticks if necessary.
1218         for stick in list_group_sub:
1219             stick.parent = sticks_empty
1220
1221         sticks_empty.location += center
1222
1223         # Collections
1224         # ===========
1225         # Create a collection that will contain all sticks + the empty and ...
1226         coll = bpy.data.collections.new("Sticks")
1227         # ... link it to the collection, which contains all parts of the
1228         # molecule.
1229         coll_molecule.children.link(coll)
1230         # Now, create a collection that only contains the sticks and ...
1231         coll_cylinder = bpy.data.collections.new("Sticks_cylinders")
1232         # ... link it to the collection, which contains the sticks and empty.
1233         coll.children.link(coll_cylinder)
1234
1235         # Note the collection where the empty was placed into, ...
1236         coll_all = sticks_empty.users_collection
1237         if len(coll_all) > 0:
1238             coll_past = coll_all[0]
1239         else:
1240             coll_past = bpy.context.scene.collection
1241         # ... link the empty with the new collection  ...
1242         coll.objects.link(sticks_empty)
1243         # ... and unlink it from the old collection where it has been before.
1244         coll_past.objects.unlink(sticks_empty)
1245
1246         # Note the collection where the cylinders were placed into, ...
1247         coll_all = list_group_sub[0].users_collection
1248         if len(coll_all) > 0:
1249             coll_past = coll_all[0]
1250         else:
1251             coll_past = bpy.context.scene.collection
1252
1253         for stick in list_group_sub:
1254             # ... link each stick with the new collection  ...
1255             coll_cylinder.objects.link(stick)
1256             # ... and unlink it from the old collection.
1257             coll_past.objects.unlink(stick)
1258
1259         return sticks_empty
1260
1261
1262 # -----------------------------------------------------------------------------
1263 #                                                            The main routine
1264
1265 def import_pdb(Ball_type,
1266                Ball_azimuth,
1267                Ball_zenith,
1268                Ball_radius_factor,
1269                radiustype,
1270                Ball_distance_factor,
1271                use_sticks,
1272                use_sticks_type,
1273                sticks_subdiv_view,
1274                sticks_subdiv_render,
1275                use_sticks_color,
1276                use_sticks_smooth,
1277                use_sticks_bonds,
1278                use_sticks_one_object,
1279                use_sticks_one_object_nr,
1280                Stick_unit, Stick_dist,
1281                Stick_sectors,
1282                Stick_diameter,
1283                put_to_center,
1284                use_camera,
1285                use_light,
1286                filepath_pdb):
1287
1288     # List of materials
1289     atom_material_list = []
1290
1291     # A list of ALL objects which are loaded (needed for selecting the loaded
1292     # structure.
1293     atom_object_list = []
1294
1295     # ------------------------------------------------------------------------
1296     # INITIALIZE THE ELEMENT LIST
1297
1298     read_elements()
1299
1300     # ------------------------------------------------------------------------
1301     # READING DATA OF ATOMS
1302
1303     (Number_of_total_atoms, all_atoms) = read_pdb_file(filepath_pdb, radiustype)
1304
1305     # ------------------------------------------------------------------------
1306     # MATERIAL PROPERTIES FOR ATOMS
1307
1308     # The list that contains info about all types of atoms is created
1309     # here. It is used for building the material properties for
1310     # instance (see below).
1311     atom_all_types_list = []
1312
1313     for atom in all_atoms:
1314         FLAG_FOUND = False
1315         for atom_type in atom_all_types_list:
1316             # If the atom name is already in the list, FLAG on 'True'.
1317             if atom_type[0] == atom.name:
1318                 FLAG_FOUND = True
1319                 break
1320         # No name in the current list has been found? => New entry.
1321         if FLAG_FOUND == False:
1322             # Stored are: Atom label (e.g. 'Na'), the corresponding atom
1323             # name (e.g. 'Sodium') and its color.
1324             atom_all_types_list.append([atom.name, atom.element, atom.color])
1325
1326     # The list of materials is built.
1327     # Note that all atoms of one type (e.g. all hydrogens) get only ONE
1328     # material! This is good because then, by activating one atom in the
1329     # Blender scene and changing the color of this atom, one changes the color
1330     # of ALL atoms of the same type at the same time.
1331
1332     # Create first a new list of materials for each type of atom
1333     # (e.g. hydrogen)
1334     for atom_type in atom_all_types_list:
1335         material = bpy.data.materials.new(atom_type[1])
1336         material.diffuse_color = atom_type[2]
1337         material.use_nodes = True
1338         mat_P_BSDF = next(n for n in material.node_tree.nodes
1339                           if n.type == "BSDF_PRINCIPLED")
1340         mat_P_BSDF.inputs['Base Color'].default_value = atom_type[2]
1341         material.name = atom_type[0]
1342         atom_material_list.append(material)
1343
1344     # Now, we go through all atoms and give them a material. For all atoms ...
1345     for atom in all_atoms:
1346         # ... and all materials ...
1347         for material in atom_material_list:
1348             # ... select the correct material for the current atom via
1349             # comparison of names ...
1350             if atom.name in material.name:
1351                 # ... and give the atom its material properties.
1352                 # However, before we check if it is a vacancy.
1353                 # The vacancy is represented by a transparent cube.
1354                 if atom.name == "Vacancy":
1355                     # For cycles and eevee.
1356                     material.use_nodes = True
1357                     mat_P_BSDF = next(n for n in material.node_tree.nodes
1358                                       if n.type == "BSDF_PRINCIPLED")
1359                     mat_P_BSDF.inputs['Metallic'].default_value = 0.1
1360                     mat_P_BSDF.inputs['Specular'].default_value = 0.15
1361                     mat_P_BSDF.inputs['Roughness'].default_value = 0.05
1362                     mat_P_BSDF.inputs['Clearcoat Roughness'].default_value = 0.37
1363                     mat_P_BSDF.inputs['IOR'].default_value = 0.8
1364                     mat_P_BSDF.inputs['Transmission'].default_value = 0.6
1365                     mat_P_BSDF.inputs['Transmission Roughness'].default_value = 0.0
1366                     mat_P_BSDF.inputs['Alpha'].default_value = 0.5
1367                     # Some additional stuff for eevee.
1368                     material.blend_method = 'HASHED'
1369                     material.shadow_method = 'HASHED'
1370                     material.use_backface_culling = False
1371                 # The atom gets its properties.
1372                 atom.material = material
1373
1374     # ------------------------------------------------------------------------
1375     # READING DATA OF STICKS
1376
1377     all_sticks = read_pdb_file_sticks(filepath_pdb,
1378                                       use_sticks_bonds,
1379                                       all_atoms)
1380     #
1381     # So far, all atoms, sticks and materials have been registered.
1382     #
1383
1384     # ------------------------------------------------------------------------
1385     # TRANSLATION OF THE STRUCTURE TO THE ORIGIN
1386
1387     # It may happen that the structure in a PDB file already has an offset
1388     # If chosen, the structure is first put into the center of the scene
1389     # (the offset is subtracted).
1390
1391     if put_to_center == True:
1392         sum_vec = Vector((0.0,0.0,0.0))
1393         # Sum of all atom coordinates
1394         sum_vec = sum([atom.location for atom in all_atoms], sum_vec)
1395         # Then the average is taken
1396         sum_vec = sum_vec / Number_of_total_atoms
1397         # After, for each atom the center of gravity is subtracted
1398         for atom in all_atoms:
1399             atom.location -= sum_vec
1400
1401     # ------------------------------------------------------------------------
1402     # SCALING
1403
1404     # Take all atoms and adjust their radii and scale the distances.
1405     for atom in all_atoms:
1406         atom.location *= Ball_distance_factor
1407
1408     # ------------------------------------------------------------------------
1409     # DETERMINATION OF SOME GEOMETRIC PROPERTIES
1410
1411     # In the following, some geometric properties of the whole object are
1412     # determined: center, size, etc.
1413     sum_vec = Vector((0.0,0.0,0.0))
1414
1415     # First the center is determined. All coordinates are summed up ...
1416     sum_vec = sum([atom.location for atom in all_atoms], sum_vec)
1417
1418     # ... and the average is taken. This gives the center of the object.
1419     object_center_vec = sum_vec / Number_of_total_atoms
1420
1421     # Now, we determine the size.The farthest atom from the object center is
1422     # taken as a measure. The size is used to place well the camera and light
1423     # into the scene.
1424     object_size_vec = [atom.location - object_center_vec for atom in all_atoms]
1425     object_size = max(object_size_vec).length
1426
1427     # ------------------------------------------------------------------------
1428     # SORTING THE ATOMS
1429
1430     # Lists of atoms of one type are created. Example:
1431     # draw_all_atoms = [ data_hydrogen,data_carbon,data_nitrogen ]
1432     # data_hydrogen = [["Hydrogen", Material_Hydrogen, Vector((x,y,z)), 109], ...]
1433
1434     # Go through the list which contains all types of atoms. It is the list,
1435     # which has been created on the top during reading the PDB file.
1436     # Example: atom_all_types_list = ["hydrogen", "carbon", ...]
1437     draw_all_atoms = []
1438     for atom_type in atom_all_types_list:
1439
1440         # Don't draw 'TER atoms'.
1441         if atom_type[0] == "TER":
1442             continue
1443
1444         # This is the draw list, which contains all atoms of one type (e.g.
1445         # all hydrogens) ...
1446         draw_all_atoms_type = []
1447
1448         # Go through all atoms ...
1449         for atom in all_atoms:
1450             # ... select the atoms of the considered type via comparison ...
1451             if atom.name == atom_type[0]:
1452                 # ... and append them to the list 'draw_all_atoms_type'.
1453                 draw_all_atoms_type.append([atom.name,
1454                                             atom.material,
1455                                             atom.location,
1456                                             atom.radius])
1457
1458         # Now append the atom list to the list of all types of atoms
1459         draw_all_atoms.append(draw_all_atoms_type)
1460
1461     # ------------------------------------------------------------------------
1462     # COLLECTION
1463
1464     # Before we start to draw the atoms and sticks, we first create a
1465     # collection for the molecule. All atoms (balls) and sticks (cylinders)
1466     # are put into this collection.
1467     coll_molecule_name = os.path.basename(filepath_pdb)
1468     scene = bpy.context.scene
1469     coll_molecule = bpy.data.collections.new(coll_molecule_name)
1470     scene.collection.children.link(coll_molecule)
1471
1472     # ------------------------------------------------------------------------
1473     # DRAWING THE ATOMS
1474
1475     bpy.ops.object.select_all(action='DESELECT')
1476
1477     list_coll_elements = []
1478     # For each list of atoms of ONE type (e.g. Hydrogen)
1479     for draw_all_atoms_type in draw_all_atoms:
1480
1481         atom_mesh, coll_element = draw_atoms_one_type(draw_all_atoms_type,
1482                                                       Ball_type,
1483                                                       Ball_azimuth,
1484                                                       Ball_zenith,
1485                                                       Ball_radius_factor,
1486                                                       object_center_vec,
1487                                                       coll_molecule)
1488         atom_object_list.append(atom_mesh)
1489         list_coll_elements.append(coll_element)
1490
1491     # ------------------------------------------------------------------------
1492     # DRAWING THE STICKS: cylinders in a dupliverts structure
1493
1494     if use_sticks == True and use_sticks_type == '0' and all_sticks != []:
1495
1496         sticks = draw_sticks_dupliverts(all_atoms,
1497                                         atom_all_types_list,
1498                                         object_center_vec,
1499                                         all_sticks,
1500                                         Stick_diameter,
1501                                         Stick_sectors,
1502                                         Stick_unit,
1503                                         Stick_dist,
1504                                         use_sticks_smooth,
1505                                         use_sticks_color,
1506                                         list_coll_elements)
1507         for stick in sticks:
1508             atom_object_list.append(stick)
1509
1510     # ------------------------------------------------------------------------
1511     # DRAWING THE STICKS: skin and subdivision modifier
1512
1513     if use_sticks == True and use_sticks_type == '1' and all_sticks != []:
1514
1515         sticks = draw_sticks_skin(all_atoms,
1516                                   all_sticks,
1517                                   Stick_diameter,
1518                                   use_sticks_smooth,
1519                                   sticks_subdiv_view,
1520                                   sticks_subdiv_render,
1521                                   coll_molecule)
1522         atom_object_list.append(sticks)
1523
1524     # ------------------------------------------------------------------------
1525     # DRAWING THE STICKS: normal cylinders
1526
1527     if use_sticks == True and use_sticks_type == '2' and all_sticks != []:
1528
1529         sticks = draw_sticks_normal(all_atoms,
1530                                     all_sticks,
1531                                     object_center_vec,
1532                                     Stick_diameter,
1533                                     Stick_sectors,
1534                                     use_sticks_smooth,
1535                                     use_sticks_one_object,
1536                                     use_sticks_one_object_nr,
1537                                     coll_molecule)
1538         atom_object_list.append(sticks)
1539
1540     # ------------------------------------------------------------------------
1541     # CAMERA and LIGHT SOURCES
1542
1543     camera_light_source(use_camera,
1544                         use_light,
1545                         object_center_vec,
1546                         object_size)
1547
1548     # ------------------------------------------------------------------------
1549     # SELECT ALL LOADED OBJECTS
1550     bpy.ops.object.select_all(action='DESELECT')
1551     obj = None
1552     for obj in atom_object_list:
1553         obj.select_set(True)
1554
1555     # activate the last selected object
1556     if obj:
1557         bpy.context.view_layer.objects.active = obj