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 stick 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                        Stick_dist,
1104                        use_sticks_smooth,
1105                        use_sticks_one_object,
1106                        use_sticks_one_object_nr,
1107                        coll_molecule):
1108
1109     stick_material = bpy.data.materials.new(ELEMENTS[-1].name)
1110     stick_material.use_nodes = True
1111     mat_P_BSDF = next(n for n in stick_material.node_tree.nodes
1112                       if n.type == "BSDF_PRINCIPLED")
1113     mat_P_BSDF.inputs['Base Color'].default_value = ELEMENTS[-1].color
1114
1115     up_axis = Vector([0.0, 0.0, 1.0])
1116
1117     # For all sticks, do ...
1118     list_group = []
1119     list_group_sub = []
1120     counter = 0
1121     for i, stick in enumerate(all_sticks):
1122
1123         # We treat here single, double and tripple bonds: stick.number <= 3
1124         for repeat in range(stick.number):
1125
1126             # The vectors of the two atoms
1127             atom1 = copy(all_atoms[stick.atom1-1].location)-center
1128             atom2 = copy(all_atoms[stick.atom2-1].location)-center
1129
1130             dist =  Stick_diameter * Stick_dist
1131
1132             # The two sticks are on the left and right of the middle connection.
1133             if stick.number == 2:
1134                 if repeat == 0:
1135                     atom1 += (stick.dist * dist)
1136                     atom2 += (stick.dist * dist)
1137                 if repeat == 1:
1138                     atom1 -= (stick.dist * dist)
1139                     atom2 -= (stick.dist * dist)
1140
1141             if stick.number == 3:
1142                 if repeat == 0:
1143                     atom1 += (stick.dist * dist)
1144                     atom2 += (stick.dist * dist)
1145                 if repeat == 2:
1146                     atom1 -= (stick.dist * dist)
1147                     atom2 -= (stick.dist * dist)
1148
1149             # Vector pointing along the stick direction
1150             dv = atom1 - atom2
1151             # The normalized vector of this, with lenght 1
1152             n  = dv / dv.length
1153             # Starting point of the stick
1154             location = (atom1 + atom2) * 0.5
1155             # Angle with respect to the z-axis
1156             angle = dv.angle(up_axis, 0)
1157             # Cross-product between v and the z-axis vector. It is the
1158             # vector of rotation.
1159             axis = up_axis.cross(dv)
1160             # Calculate Euler angles
1161             euler = Matrix.Rotation(angle, 4, axis).to_euler()
1162             # Create stick
1163             stick_obj = bpy.ops.mesh.primitive_cylinder_add(vertices=Stick_sectors,
1164                                                             radius=Stick_diameter,
1165                                                             depth=dv.length,
1166                                                             end_fill_type='NGON',
1167                                                             align='WORLD',
1168                                                             enter_editmode=False,
1169                                                             location=location,
1170                                                             rotation=(0, 0, 0))
1171             # Put the stick into the scene ...
1172             stick_obj = bpy.context.view_layer.objects.active
1173             # ... and rotate the stick.
1174             stick_obj.rotation_euler = euler
1175             # ... and name
1176             if stick.number == 1:
1177                 stick_obj.name = "Stick_Cylinder_%04d" %(i)
1178             elif stick.number == 2:
1179                 if repeat == 0:
1180                     stick_obj.name = "Stick_Cylinder_%04d" %(i) + "_left"
1181                 elif repeat == 1:
1182                     stick_obj.name = "Stick_Cylinder_%04d" %(i) + "_right"
1183             elif stick.number == 3:
1184                 if repeat == 0:
1185                     stick_obj.name = "Stick_Cylinder_%04d" %(i) + "_left"
1186                 elif repeat == 1:
1187                     stick_obj.name = "Stick_Cylinder_%04d" %(i) + "_middle"
1188                 elif repeat == 2:
1189                     stick_obj.name = "Stick_Cylinder_%04d" %(i) + "_right"
1190                 # Never occurs:
1191                 else:
1192                     stick_obj.name = "Stick_Cylinder"
1193             # Never occurs:
1194             else:
1195                 stick_obj.name = "Stick_Cylinder"
1196             counter += 1
1197
1198             # Smooth the cylinder.
1199             if use_sticks_smooth == True:
1200                 bpy.ops.object.select_all(action='DESELECT')
1201                 stick_obj.select_set(True)
1202                 bpy.ops.object.shade_smooth()
1203
1204             list_group_sub.append(stick_obj)
1205
1206             if use_sticks_one_object == True:
1207                 if counter == use_sticks_one_object_nr:
1208                     bpy.ops.object.select_all(action='DESELECT')
1209                     for stick_select in list_group_sub:
1210                         stick_select.select_set(True)
1211                     bpy.ops.object.join()
1212                     list_group.append(bpy.context.view_layer.objects.active)
1213                     bpy.ops.object.select_all(action='DESELECT')
1214                     list_group_sub = []
1215                     counter = 0
1216             else:
1217                 # Material ...
1218                 stick_obj.active_material = stick_material
1219
1220     if use_sticks_one_object == True:
1221         bpy.ops.object.select_all(action='DESELECT')
1222         for stick in list_group_sub:
1223             stick.select_set(True)
1224         bpy.ops.object.join()
1225         list_group.append(bpy.context.view_layer.objects.active)
1226         bpy.ops.object.select_all(action='DESELECT')
1227
1228         for group in list_group:
1229             group.select_set(True)
1230         bpy.ops.object.join()
1231         bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY',
1232                                   center='MEDIAN')
1233         sticks = bpy.context.view_layer.objects.active
1234         sticks.active_material = stick_material
1235
1236         sticks.location += center
1237
1238         # Collections
1239         # ===========
1240         # Note the collection where the sticks were placed into.
1241         coll_all = sticks.users_collection
1242         if len(coll_all) > 0:
1243             coll_past = coll_all[0]
1244         else:
1245             coll_past = bpy.context.scene.collection
1246
1247         # Link the sticks with the collection of the molecule ...
1248         coll_molecule.objects.link(sticks)
1249         # ... and unlink them from the collection it has been before.
1250         coll_past.objects.unlink(sticks)
1251
1252         return sticks
1253     else:
1254         # Here we use an empty ...
1255         bpy.ops.object.empty_add(type='ARROWS',
1256                                   align='WORLD',
1257                                   location=(0, 0, 0),
1258                                   rotation=(0, 0, 0))
1259         sticks_empty = bpy.context.view_layer.objects.active
1260         sticks_empty.name = "A_sticks_empty"
1261         # ... that is parent to all sticks. With this, we can better move
1262         # all sticks if necessary.
1263         for stick in list_group_sub:
1264             stick.parent = sticks_empty
1265
1266         sticks_empty.location += center
1267
1268         # Collections
1269         # ===========
1270         # Create a collection that will contain all sticks + the empty and ...
1271         coll = bpy.data.collections.new("Sticks")
1272         # ... link it to the collection, which contains all parts of the
1273         # molecule.
1274         coll_molecule.children.link(coll)
1275         # Now, create a collection that only contains the sticks and ...
1276         coll_cylinder = bpy.data.collections.new("Sticks_cylinders")
1277         # ... link it to the collection, which contains the sticks and empty.
1278         coll.children.link(coll_cylinder)
1279
1280         # Note the collection where the empty was placed into, ...
1281         coll_all = sticks_empty.users_collection
1282         if len(coll_all) > 0:
1283             coll_past = coll_all[0]
1284         else:
1285             coll_past = bpy.context.scene.collection
1286         # ... link the empty with the new collection  ...
1287         coll.objects.link(sticks_empty)
1288         # ... and unlink it from the old collection where it has been before.
1289         coll_past.objects.unlink(sticks_empty)
1290
1291         # Note the collection where the cylinders were placed into, ...
1292         coll_all = list_group_sub[0].users_collection
1293         if len(coll_all) > 0:
1294             coll_past = coll_all[0]
1295         else:
1296             coll_past = bpy.context.scene.collection
1297
1298         for stick in list_group_sub:
1299             # ... link each stick with the new collection  ...
1300             coll_cylinder.objects.link(stick)
1301             # ... and unlink it from the old collection.
1302             coll_past.objects.unlink(stick)
1303
1304         return sticks_empty
1305
1306
1307 # -----------------------------------------------------------------------------
1308 #                                                            The main routine
1309
1310 def import_pdb(Ball_type,
1311                Ball_azimuth,
1312                Ball_zenith,
1313                Ball_radius_factor,
1314                radiustype,
1315                Ball_distance_factor,
1316                use_sticks,
1317                use_sticks_type,
1318                sticks_subdiv_view,
1319                sticks_subdiv_render,
1320                use_sticks_color,
1321                use_sticks_smooth,
1322                use_sticks_bonds,
1323                use_sticks_one_object,
1324                use_sticks_one_object_nr,
1325                Stick_unit, Stick_dist,
1326                Stick_sectors,
1327                Stick_diameter,
1328                put_to_center,
1329                use_camera,
1330                use_light,
1331                filepath_pdb):
1332
1333     # List of materials
1334     atom_material_list = []
1335
1336     # A list of ALL objects which are loaded (needed for selecting the loaded
1337     # structure.
1338     atom_object_list = []
1339
1340     # ------------------------------------------------------------------------
1341     # INITIALIZE THE ELEMENT LIST
1342
1343     read_elements()
1344
1345     # ------------------------------------------------------------------------
1346     # READING DATA OF ATOMS
1347
1348     (Number_of_total_atoms, all_atoms) = read_pdb_file(filepath_pdb, radiustype)
1349
1350     # ------------------------------------------------------------------------
1351     # MATERIAL PROPERTIES FOR ATOMS
1352
1353     # The list that contains info about all types of atoms is created
1354     # here. It is used for building the material properties for
1355     # instance (see below).
1356     atom_all_types_list = []
1357
1358     for atom in all_atoms:
1359         FLAG_FOUND = False
1360         for atom_type in atom_all_types_list:
1361             # If the atom name is already in the list, FLAG on 'True'.
1362             if atom_type[0] == atom.name:
1363                 FLAG_FOUND = True
1364                 break
1365         # No name in the current list has been found? => New entry.
1366         if FLAG_FOUND == False:
1367             # Stored are: Atom label (e.g. 'Na'), the corresponding atom
1368             # name (e.g. 'Sodium') and its color.
1369             atom_all_types_list.append([atom.name, atom.element, atom.color])
1370
1371     # The list of materials is built.
1372     # Note that all atoms of one type (e.g. all hydrogens) get only ONE
1373     # material! This is good because then, by activating one atom in the
1374     # Blender scene and changing the color of this atom, one changes the color
1375     # of ALL atoms of the same type at the same time.
1376
1377     # Create first a new list of materials for each type of atom
1378     # (e.g. hydrogen)
1379     for atom_type in atom_all_types_list:
1380         material = bpy.data.materials.new(atom_type[1])
1381         material.diffuse_color = atom_type[2]
1382         material.use_nodes = True
1383         mat_P_BSDF = next(n for n in material.node_tree.nodes
1384                           if n.type == "BSDF_PRINCIPLED")
1385         mat_P_BSDF.inputs['Base Color'].default_value = atom_type[2]
1386         material.name = atom_type[0]
1387         atom_material_list.append(material)
1388
1389     # Now, we go through all atoms and give them a material. For all atoms ...
1390     for atom in all_atoms:
1391         # ... and all materials ...
1392         for material in atom_material_list:
1393             # ... select the correct material for the current atom via
1394             # comparison of names ...
1395             if atom.name in material.name:
1396                 # ... and give the atom its material properties.
1397                 # However, before we check if it is a vacancy.
1398                 # The vacancy is represented by a transparent cube.
1399                 if atom.name == "Vacancy":
1400                     # For cycles and eevee.
1401                     material.use_nodes = True
1402                     mat_P_BSDF = next(n for n in material.node_tree.nodes
1403                                       if n.type == "BSDF_PRINCIPLED")
1404                     mat_P_BSDF.inputs['Metallic'].default_value = 0.1
1405                     mat_P_BSDF.inputs['Specular'].default_value = 0.15
1406                     mat_P_BSDF.inputs['Roughness'].default_value = 0.05
1407                     mat_P_BSDF.inputs['Clearcoat Roughness'].default_value = 0.37
1408                     mat_P_BSDF.inputs['IOR'].default_value = 0.8
1409                     mat_P_BSDF.inputs['Transmission'].default_value = 0.6
1410                     mat_P_BSDF.inputs['Transmission Roughness'].default_value = 0.0
1411                     mat_P_BSDF.inputs['Alpha'].default_value = 0.5
1412                     # Some additional stuff for eevee.
1413                     material.blend_method = 'HASHED'
1414                     material.shadow_method = 'HASHED'
1415                     material.use_backface_culling = False
1416                 # The atom gets its properties.
1417                 atom.material = material
1418
1419     # ------------------------------------------------------------------------
1420     # READING DATA OF STICKS
1421
1422     all_sticks = read_pdb_file_sticks(filepath_pdb,
1423                                       use_sticks_bonds,
1424                                       all_atoms)
1425     #
1426     # So far, all atoms, sticks and materials have been registered.
1427     #
1428
1429     # ------------------------------------------------------------------------
1430     # TRANSLATION OF THE STRUCTURE TO THE ORIGIN
1431
1432     # It may happen that the structure in a PDB file already has an offset
1433     # If chosen, the structure is first put into the center of the scene
1434     # (the offset is subtracted).
1435
1436     if put_to_center == True:
1437         sum_vec = Vector((0.0,0.0,0.0))
1438         # Sum of all atom coordinates
1439         sum_vec = sum([atom.location for atom in all_atoms], sum_vec)
1440         # Then the average is taken
1441         sum_vec = sum_vec / Number_of_total_atoms
1442         # After, for each atom the center of gravity is subtracted
1443         for atom in all_atoms:
1444             atom.location -= sum_vec
1445
1446     # ------------------------------------------------------------------------
1447     # SCALING
1448
1449     # Take all atoms and adjust their radii and scale the distances.
1450     for atom in all_atoms:
1451         atom.location *= Ball_distance_factor
1452
1453     # ------------------------------------------------------------------------
1454     # DETERMINATION OF SOME GEOMETRIC PROPERTIES
1455
1456     # In the following, some geometric properties of the whole object are
1457     # determined: center, size, etc.
1458     sum_vec = Vector((0.0,0.0,0.0))
1459
1460     # First the center is determined. All coordinates are summed up ...
1461     sum_vec = sum([atom.location for atom in all_atoms], sum_vec)
1462
1463     # ... and the average is taken. This gives the center of the object.
1464     object_center_vec = sum_vec / Number_of_total_atoms
1465
1466     # Now, we determine the size.The farthest atom from the object center is
1467     # taken as a measure. The size is used to place well the camera and light
1468     # into the scene.
1469     object_size_vec = [atom.location - object_center_vec for atom in all_atoms]
1470     object_size = max(object_size_vec).length
1471
1472     # ------------------------------------------------------------------------
1473     # SORTING THE ATOMS
1474
1475     # Lists of atoms of one type are created. Example:
1476     # draw_all_atoms = [ data_hydrogen,data_carbon,data_nitrogen ]
1477     # data_hydrogen = [["Hydrogen", Material_Hydrogen, Vector((x,y,z)), 109], ...]
1478
1479     # Go through the list which contains all types of atoms. It is the list,
1480     # which has been created on the top during reading the PDB file.
1481     # Example: atom_all_types_list = ["hydrogen", "carbon", ...]
1482     draw_all_atoms = []
1483     for atom_type in atom_all_types_list:
1484
1485         # Don't draw 'TER atoms'.
1486         if atom_type[0] == "TER":
1487             continue
1488
1489         # This is the draw list, which contains all atoms of one type (e.g.
1490         # all hydrogens) ...
1491         draw_all_atoms_type = []
1492
1493         # Go through all atoms ...
1494         for atom in all_atoms:
1495             # ... select the atoms of the considered type via comparison ...
1496             if atom.name == atom_type[0]:
1497                 # ... and append them to the list 'draw_all_atoms_type'.
1498                 draw_all_atoms_type.append([atom.name,
1499                                             atom.material,
1500                                             atom.location,
1501                                             atom.radius])
1502
1503         # Now append the atom list to the list of all types of atoms
1504         draw_all_atoms.append(draw_all_atoms_type)
1505
1506     # ------------------------------------------------------------------------
1507     # COLLECTION
1508
1509     # Before we start to draw the atoms and sticks, we first create a
1510     # collection for the molecule. All atoms (balls) and sticks (cylinders)
1511     # are put into this collection.
1512     coll_molecule_name = os.path.basename(filepath_pdb)
1513     scene = bpy.context.scene
1514     coll_molecule = bpy.data.collections.new(coll_molecule_name)
1515     scene.collection.children.link(coll_molecule)
1516
1517     # ------------------------------------------------------------------------
1518     # DRAWING THE ATOMS
1519
1520     bpy.ops.object.select_all(action='DESELECT')
1521
1522     list_coll_elements = []
1523     # For each list of atoms of ONE type (e.g. Hydrogen)
1524     for draw_all_atoms_type in draw_all_atoms:
1525
1526         atom_mesh, coll_element = draw_atoms_one_type(draw_all_atoms_type,
1527                                                       Ball_type,
1528                                                       Ball_azimuth,
1529                                                       Ball_zenith,
1530                                                       Ball_radius_factor,
1531                                                       object_center_vec,
1532                                                       coll_molecule)
1533         atom_object_list.append(atom_mesh)
1534         list_coll_elements.append(coll_element)
1535
1536     # ------------------------------------------------------------------------
1537     # DRAWING THE STICKS: cylinders in a dupliverts structure
1538
1539     if use_sticks == True and use_sticks_type == '0' and all_sticks != []:
1540
1541         sticks = draw_sticks_dupliverts(all_atoms,
1542                                         atom_all_types_list,
1543                                         object_center_vec,
1544                                         all_sticks,
1545                                         Stick_diameter,
1546                                         Stick_sectors,
1547                                         Stick_unit,
1548                                         Stick_dist,
1549                                         use_sticks_smooth,
1550                                         use_sticks_color,
1551                                         list_coll_elements)
1552         for stick in sticks:
1553             atom_object_list.append(stick)
1554
1555     # ------------------------------------------------------------------------
1556     # DRAWING THE STICKS: skin and subdivision modifier
1557
1558     if use_sticks == True and use_sticks_type == '1' and all_sticks != []:
1559
1560         sticks = draw_sticks_skin(all_atoms,
1561                                   all_sticks,
1562                                   Stick_diameter,
1563                                   use_sticks_smooth,
1564                                   sticks_subdiv_view,
1565                                   sticks_subdiv_render,
1566                                   coll_molecule)
1567         atom_object_list.append(sticks)
1568
1569     # ------------------------------------------------------------------------
1570     # DRAWING THE STICKS: normal cylinders
1571
1572     if use_sticks == True and use_sticks_type == '2' and all_sticks != []:
1573
1574         sticks = draw_sticks_normal(all_atoms,
1575                                     all_sticks,
1576                                     object_center_vec,
1577                                     Stick_diameter,
1578                                     Stick_sectors,
1579                                     Stick_dist,
1580                                     use_sticks_smooth,
1581                                     use_sticks_one_object,
1582                                     use_sticks_one_object_nr,
1583                                     coll_molecule)
1584         atom_object_list.append(sticks)
1585
1586     # ------------------------------------------------------------------------
1587     # CAMERA and LIGHT SOURCES
1588
1589     camera_light_source(use_camera,
1590                         use_light,
1591                         object_center_vec,
1592                         object_size)
1593
1594     # ------------------------------------------------------------------------
1595     # SELECT ALL LOADED OBJECTS
1596     bpy.ops.object.select_all(action='DESELECT')
1597     obj = None
1598     for obj in atom_object_list:
1599         obj.select_set(True)
1600
1601     # activate the last selected object
1602     if obj:
1603         bpy.context.view_layer.objects.active = obj