io_mesh_atomic/pdb_import.py

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