path: root/io_mesh_pdb/import_pdb.py

# ##### BEGIN GPL LICENSE BLOCK #####
#
#  This program is free software; you can redistribute it and/or
#  modify it under the terms of the GNU General Public License
#  as published by the Free Software Foundation; either version 2
#  of the License, or (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software Foundation,
#  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# ##### END GPL LICENSE BLOCK #####

# <pep8 compliant>

# TODO, currently imported names are assumed to make valid blender names
#       this is _not_ assured, so we should use a reliable dict mapping.

import bpy


class Element:
    """Element class with properties ([R, G, B], cov_radius, vdw_radius, name)"""
    def __init__(self, color, cov_radius, vdw_radius, name):
        self.color = color
        self.cov_radius = cov_radius
        self.vdw_radius = vdw_radius
        self.name = name


class Atom:
    """Atom class with properties (serial, name, altloc, resname,chainid,
        resseq, icode, x, y, z, occupancy, tempfactor, element, charge)"""

    def __init__(self, serial, name, altloc, resname, chainid, resseq, icode,
                 x, y, z, occupancy, tempfactor, element, charge):
        self.serial = serial
        self.name = name
        self.altloc = altloc
        self.resname = resname
        self.chainid = chainid
        self.resseq = resseq
        self.icode = icode
        self.x = x
        self.y = y
        self.z = z
        self.occupancy = occupancy
        self.tempfactor = tempfactor
        self.element = element
        self.charge = charge


# collection of biomolecules based on model
# all chains in model stored here
class Model:
    '''Model class'''
    def __init__(self, model_id):
        self.model_id = model_id
        self.atoms = {}
        self.atom_count = 0
        self.vert_list = []
        # Dictionary of {vert index: [list of vertex groups it belongs to]}
        # Now element only
        self.vert_group_index = {}
        # Dictionary of {vertex group: number of group members}
        self.vert_group_counts = {}
        self.chains = {}


# new object level class
class Biomolecule:
    '''Biomolecule'''
    def __init__(self, serial):
        self.serial = serial
        self.atom_count = 0
        self.vert_list = []
        self.vert_group_index = {}
        self.vert_group_counts = {}
        self.chain_transforms = {}


# Atom collection
class Chain:
    '''Chain'''
    def __init__(self, chain_id):
        self.chain_id = chain_id
        self.atoms = {}

# Atomic data from http://www.ccdc.cam.ac.uk/products/csd/radii/
# Color palatte adapted from Jmol
# "Element symbol":[[Red, Green, Blue], Covalent radius, van der Waals radius,
#                   Element name]
# Atomic radii are in angstroms (100 pm)
# Unknown covalent radii are assigned 1.5 A, unknown van der Waals radiii are
# assigned 2 A,

atom_data = {
    "H" : Element((1.00000, 1.00000, 1.00000), 0.23, 1.09, "Hydrogen"     ),
    "HE": Element((0.85098, 1.00000, 1.00000), 1.5 , 1.4 , "Helium"       ),
    "LI": Element((0.80000, 0.50196, 1.00000), 1.28, 1.82, "Lithium"      ),
    "BE": Element((0.76078, 1.00000, 0.00000), 0.96, 2   , "Beryllium"    ),
    "B" : Element((1.00000, 0.70980, 0.70980), 0.83, 2   , "Boron"        ),
    "C" : Element((0.56471, 0.56471, 0.56471), 0.68, 1.7 , "Carbon"       ),
    "N" : Element((0.18824, 0.31373, 0.97255), 0.68, 1.55, "Nitrogen"     ),
    "O" : Element((1.00000, 0.05098, 0.05098), 0.68, 1.52, "Oxygen"       ),
    "F" : Element((0.56471, 0.87843, 0.31373), 0.64, 1.47, "Fluorine"     ),
    "NE": Element((0.70196, 0.89020, 0.96078), 1.5 , 1.54, "Neon"         ),
    "NA": Element((0.67059, 0.36078, 0.94902), 1.66, 2.27, "Sodium"       ),
    "MG": Element((0.54118, 1.00000, 0.00000), 1.41, 1.73, "Magnesium"    ),
    "AL": Element((0.74902, 0.65098, 0.65098), 1.21, 2   , "Aluminum"     ),
    "SI": Element((0.94118, 0.78431, 0.62745), 1.2 , 2.1 , "Silicon"      ),
    "P" : Element((1.00000, 0.50196, 0.00000), 1.05, 1.8 , "Phosphorus"   ),
    "S" : Element((1.00000, 1.00000, 0.18824), 1.02, 1.8 , "Sulfur"       ),
    "CL": Element((0.12157, 0.94118, 0.12157), 0.99, 1.75, "Chlorine"     ),
    "AR": Element((0.50196, 0.81961, 0.89020), 1.51, 1.88, "Argon"        ),
    "K" : Element((0.56078, 0.25098, 0.83137), 2.03, 2.75, "Potassium"    ),
    "CA": Element((0.23922, 1.00000, 0.00000), 1.76, 2   , "Calcium"      ),
    "SC": Element((0.90196, 0.90196, 0.90196), 1.7 , 2   , "Scandium"     ),
    "TI": Element((0.74902, 0.76078, 0.78039), 1.6 , 2   , "Titanium"     ),
    "V" : Element((0.65098, 0.65098, 0.67059), 1.53, 2   , "Vanadium"     ),
    "CR": Element((0.54118, 0.60000, 0.78039), 1.39, 2   , "Chromium"     ),
    "MN": Element((0.61176, 0.47843, 0.78039), 1.61, 2   , "Manganese"    ),
    "FE": Element((0.87843, 0.40000, 0.20000), 1.52, 2   , "Iron"         ),
    "CO": Element((0.94118, 0.56471, 0.62745), 1.26, 2   , "Cobalt"       ),
    "NI": Element((0.31373, 0.81569, 0.31373), 1.24, 1.63, "Nickel"       ),
    "CU": Element((0.78431, 0.50196, 0.20000), 1.32, 1.4 , "Copper"       ),
    "ZN": Element((0.49020, 0.50196, 0.69020), 1.22, 1.39, "Zinc"         ),
    "GA": Element((0.76078, 0.56078, 0.56078), 1.22, 1.87, "Gallium"      ),
    "GE": Element((0.40000, 0.56078, 0.56078), 1.17, 2   , "Germanium"    ),
    "AS": Element((0.74118, 0.50196, 0.89020), 1.21, 1.85, "Arsenic"      ),
    "SE": Element((1.00000, 0.63137, 0.00000), 1.22, 1.9 , "Selenium"     ),
    "BR": Element((0.65098, 0.16078, 0.16078), 1.21, 1.85, "Bromine"      ),
    "KR": Element((0.36078, 0.72157, 0.81961), 1.5 , 2.02, "Krypton"      ),
    "RB": Element((0.43922, 0.18039, 0.69020), 2.2 , 2   , "Rubidium"     ),
    "SR": Element((0.00000, 1.00000, 0.00000), 1.95, 2   , "Strontium"    ),
    "Y" : Element((0.58039, 1.00000, 1.00000), 1.9 , 2   , "Yttrium"      ),
    "ZR": Element((0.58039, 0.87843, 0.87843), 1.75, 2   , "Zirconium"    ),
    "NB": Element((0.45098, 0.76078, 0.78824), 1.64, 2   , "Niobium"      ),
    "MO": Element((0.32941, 0.70980, 0.70980), 1.54, 2   , "Molybdenum"   ),
    "TC": Element((0.23137, 0.61961, 0.61961), 1.47, 2   , "Technetium"   ),
    "RU": Element((0.14118, 0.56078, 0.56078), 1.46, 2   , "Ruthenium"    ),
    "RH": Element((0.03922, 0.49020, 0.54902), 1.45, 2   , "Rhodium"      ),
    "PD": Element((0.00000, 0.41176, 0.52157), 1.39, 1.63, "Palladium"    ),
    "AG": Element((0.75294, 0.75294, 0.75294), 1.45, 1.72, "Silver"       ),
    "CD": Element((1.00000, 0.85098, 0.56078), 1.44, 1.58, "Cadmium"      ),
    "IN": Element((0.65098, 0.45882, 0.45098), 1.42, 1.93, "Indium"       ),
    "SN": Element((0.40000, 0.50196, 0.50196), 1.39, 2.17, "Tin"          ),
    "SB": Element((0.61961, 0.38824, 0.70980), 1.39, 2   , "Antimony"     ),
    "TE": Element((0.83137, 0.47843, 0.00000), 1.47, 2.06, "Tellurium"    ),
    "I" : Element((0.58039, 0.00000, 0.58039), 1.4 , 1.98, "Iodine"       ),
    "XE": Element((0.25882, 0.61961, 0.69020), 1.5 , 2.16, "Xenon"        ),
    "CS": Element((0.34118, 0.09020, 0.56078), 2.44, 2   , "Cesium"       ),
    "BA": Element((0.00000, 0.78824, 0.00000), 2.15, 2   , "Barium"       ),
    "LA": Element((0.43922, 0.83137, 1.00000), 2.07, 2   , "Lanthanum"    ),
    "CE": Element((1.00000, 1.00000, 0.78039), 2.04, 2   , "Cerium"       ),
    "PR": Element((0.85098, 1.00000, 0.78039), 2.03, 2   , "Praseodymium" ),
    "ND": Element((0.78039, 1.00000, 0.78039), 2.01, 2   , "Neodymium"    ),
    "PM": Element((0.63922, 1.00000, 0.78039), 1.99, 2   , "Promethium"   ),
    "SM": Element((0.56078, 1.00000, 0.78039), 1.98, 2   , "Samarium"     ),
    "EE": Element((0.38039, 1.00000, 0.78039), 1.98, 2   , "Europium"     ),
    "GD": Element((0.27059, 1.00000, 0.78039), 1.96, 2   , "Gadolinium"   ),
    "TB": Element((0.18824, 1.00000, 0.78039), 1.94, 2   , "Terbium"      ),
    "DY": Element((0.12157, 1.00000, 0.78039), 1.92, 2   , "Dysprosium"   ),
    "HO": Element((0.00000, 1.00000, 0.61176), 1.92, 2   , "Holmium"      ),
    "ER": Element((0.00000, 0.90196, 0.45882), 1.89, 2   , "Erbium"       ),
    "TM": Element((0.00000, 0.83137, 0.32157), 1.9 , 2   , "Thulium"      ),
    "YB": Element((0.00000, 0.74902, 0.21961), 1.87, 2   , "Ytterbium"    ),
    "LU": Element((0.00000, 0.67059, 0.14118), 1.87, 2   , "Lutetium"     ),
    "HF": Element((0.30196, 0.76078, 1.00000), 1.75, 2   , "Hafnium"      ),
    "TA": Element((0.30196, 0.65098, 1.00000), 1.7 , 2   , "Tantalum"     ),
    "W" : Element((0.12941, 0.58039, 0.83922), 1.62, 2   , "Tungsten"     ),
    "RE": Element((0.14902, 0.49020, 0.67059), 1.51, 2   , "Rhenium"      ),
    "OS": Element((0.14902, 0.40000, 0.58824), 1.44, 2   , "Osmium"       ),
    "IR": Element((0.09020, 0.32941, 0.52941), 1.41, 2   , "Iridium"      ),
    "PT": Element((0.81569, 0.81569, 0.87843), 1.36, 1.72, "Platinum"     ),
    "AU": Element((1.00000, 0.81961, 0.13725), 1.5 , 1.66, "Gold"         ),
    "HG": Element((0.72157, 0.72157, 0.81569), 1.32, 1.55, "Mercury"      ),
    "TL": Element((0.65098, 0.32941, 0.30196), 1.45, 1.96, "Thallium"     ),
    "PB": Element((0.34118, 0.34902, 0.38039), 1.46, 2.02, "Lead"         ),
    "BI": Element((0.61961, 0.30980, 0.70980), 1.48, 2   , "Bismuth"      ),
    "PO": Element((0.67059, 0.36078, 0.00000), 1.4 , 2   , "Polonium"     ),
    "AT": Element((0.45882, 0.30980, 0.27059), 1.21, 2   , "Astatine"     ),
    "RN": Element((0.25882, 0.50980, 0.58824), 1.5 , 2   , "Radon"        ),
    "FR": Element((0.25882, 0.00000, 0.40000), 2.6 , 2   , "Francium"     ),
    "RA": Element((0.00000, 0.49020, 0.00000), 2.21, 2   , "Radium"       ),
    "AC": Element((0.43922, 0.67059, 0.98039), 2.15, 2   , "Actinium"     ),
    "TH": Element((0.00000, 0.72941, 1.00000), 2.06, 2   , "Thorium"      ),
    "PA": Element((0.00000, 0.63137, 1.00000), 2   , 2   , "Protactinium" ),
    "U" : Element((0.00000, 0.56078, 1.00000), 1.96, 1.86, "Uranium"      ),
    "NP": Element((0.00000, 0.50196, 1.00000), 1.9 , 2   , "Neptunium"    ),
    "PU": Element((0.00000, 0.41961, 1.00000), 1.87, 2   , "Plutonium"    ),
    "AM": Element((0.32941, 0.36078, 0.94902), 1.8 , 2   , "Americium"    ),
    "CM": Element((0.47059, 0.36078, 0.89020), 1.69, 2   , "Curium"       ),
    "BK": Element((0.54118, 0.30980, 0.89020), 1.54, 2   , "Berkelium"    ),
    "CF": Element((0.63137, 0.21176, 0.83137), 1.83, 2   , "Californium"  ),
    "ES": Element((0.70196, 0.12157, 0.83137), 1.5 , 2   , "Einsteinium"  ),
    "FM": Element((0.70196, 0.12157, 0.72941), 1.5 , 2   , "Fermium"      ),
    "MD": Element((0.70196, 0.05098, 0.65098), 1.5 , 2   , "Mendelevium"  ),
    "NO": Element((0.74118, 0.05098, 0.52941), 1.5 , 2   , "Nobelium"     ),
    "LR": Element((0.78039, 0.00000, 0.40000), 1.5 , 2   , "Lawrencium"   ),
    "RF": Element((0.80000, 0.00000, 0.34902), 1.5 , 2   , "Rutherfordium"),
    "DB": Element((0.81961, 0.00000, 0.30980), 1.5 , 2   , "Dubnium"      ),
    "SG": Element((0.85098, 0.00000, 0.27059), 1.5 , 2   , "Seaborgium"   ),
    "BH": Element((0.87843, 0.00000, 0.21961), 1.5 , 2   , "Bohrium"      ),
    "HS": Element((0.90196, 0.00000, 0.18039), 1.5 , 2   , "Hassium"      ),
    "MT": Element((0.92157, 0.00000, 0.14902), 1.5 , 2   , "Meitnerium"   ),
    "DS": Element((0.93725, 0.00000, 0.12157), 1.5 , 2   , "Darmstadtium" )
}


def load_pdb(context,
             filepath="",
             atom_size=1.0,
             scene_scale=1.0,
             atom_subdivisions=3,
             retain_alts=False,
             multi_models=False,
             multimers=False,
             ):

    scene = context.scene

    file = open(filepath, 'r')

    model_list = {}
    model_flag = False
    biomolecule_flag = False
    biomolecule_list = {}
    chain_list = []
    mat_list = []

    # -------------------------------------------------------------------------
    # Parse data

    for line in file:
        # print(line)
        if line[:6] == 'COMPND':
            if line[11:17] == 'CHAIN:':
                s = 17
                for i in range(1, (len(line[17:]) + 1)):
                    if line[16 + i:17 + i] == ',':
                        chain_id = line[s:16 + i].strip()
                        chain_list.append(chain_id)
                        s = 17 + i
                    elif line[16 + i:17 + i] == ';':
                        chain_id = line[s:16 + i].strip()
                        chain_list.append(chain_id)
                        break
                    elif i == len(line[17:]):
                        chain_id = line[s:].strip()
                        chain_list.append(chain_id)
        elif line[:10] == 'REMARK 300':
            if line[11:23] == 'BIOMOLECULE:':
                biomolecule_flag = True
                s = 23
                for i in range(1, (len(line[23:]) + 1)):
                    if line[22 + i:23 + i] == ',':
                        bm_serial = int(line[s:22 + i])
                        biomolecule_list[bm_serial] = Biomolecule(bm_serial)
                        s = 23 + i
                    elif i == len(line[23:]):
                        bm_serial = int(line[s:])
                        biomolecule_list[bm_serial] = Biomolecule(bm_serial)

        elif line[:10] == 'REMARK 350':
            if line[11:23] == 'BIOMOLECULE:':
                cur_biomolecule = biomolecule_list[int(line[23:].strip())]
            elif line[11:41] == 'APPLY THE FOLLOWING TO CHAINS:':
                s = 41
                cur_chain_list = []
                for i in range(1, (len(line[41:]) + 1)):
                    if line[40 + i:41 + i] == ',':
                        cur_chain_list.append(line[s:40 + i].strip())
                        s = 41 + i
                    elif i == len(line[41:]):
                        cur_chain_list.append(line[s:].strip())
                cur_biomolecule.chain_transforms[tuple(cur_chain_list)] = []
            elif line[13:18] == 'BIOMT':
                if line[18:19] == '1':
                    row1 = [float(line[24:33]), float(line[34:43]),
                            float(line[44:53]), float(line[60:68])]
                elif line[18:19] == '2':
                    row2 = [float(line[24:33]), float(line[34:43]),
                            float(line[44:53]), float(line[60:68])]
                elif line[18:19] == '3':
                    row3 = [float(line[24:33]), float(line[34:43]),
                            float(line[44:53]), float(line[60:68])]
                    cur_biomolecule.chain_transforms[tuple(cur_chain_list)].append([row1, row2, row3])

        elif line[:5] == 'MODEL':
            model_id = int(line[10:14])
            model_list[model_id] = Model(model_id)
            cur_model = model_list[model_id]
            model_flag = True
            for chain in chain_list:
                cur_model.chains[chain] = Chain(chain)
        elif line[:6] == 'ENDMDL':
            cur_model = None
        elif line[:4] == 'ATOM':
            if not model_flag:
                model_list[1] = Model(1)
                cur_model = model_list[1]
                model_flag = True
                for chain in chain_list:
                    cur_model.chains[chain] = Chain(chain)

            serial = int(line[6:11])
            name = line[12:16].strip()
            altloc = line[16:17]
            resname = line[17:20]
            chainid = line[21:22]
            resseq = int(line[22:26])
            icode = line[26:28].strip()
            x = float(line[30:38]) * scene_scale
            y = float(line[38:46]) * scene_scale
            z = float(line[46:54]) * scene_scale
            occupancy = float(line[54:60])
            tempfactor = float(line[60:66])
            element = line[76:78].strip()
            charge = line[78:80].strip()

            '''
            print('******')
            print('serial       : ' )
            print(serial)
            print('name         : ' )
            print(name)
            print('altloc       : ' )
            print(altloc)
            print('resname      : ' )
            print(resname)
            print('chainid      : ' )
            print(chainid)
            print('resseq       : ' )
            print(resseq)
            print('icode        : ' )
            print(icode)
            print('x            : ' )
            print(x)
            print('y            : ' )
            print(y)
            print('z            : ' )
            print(z)
            print('occupancy    : ' )
            print(occupancy)
            print('tempfactor   : ' )
            print(tempfactor)
            print('element      : ' )
            print(element)
            print('charge       : ' )
            print(charge)
            print('******')
            '''

            cur_model.chains[chainid].atoms[serial] = Atom(serial, name, altloc,
                                                           resname, chainid,
                                                           resseq, icode,
                                                           x, y, z, occupancy,
                                                           tempfactor, element,
                                                           charge)

    file.close()

    if (not multimers) or (not biomolecule_flag):
        # Create a default biomolecule w/ all chains and identity transform
        # Overwrites original biomolecule_list
        biomolecule_flag = True
        bm_serial = 1
        biomolecule = Biomolecule(bm_serial)
        biomolecule.chain_transforms[tuple(chain_list)] = [(((1.0, 0.0, 0.0, 0.0),
                                                             (0.0, 1.0, 0.0, 0.0),
                                                             (0.0, 0.0, 1.0, 0.0)))]
        biomolecule_list = {bm_serial: biomolecule}

    # -------------------------------------------------------------------------
    # Create atom mesh template
    bpy.ops.object.select_all(action='DESELECT')

    bpy.ops.mesh.primitive_ico_sphere_add(subdivisions=atom_subdivisions,
                                          size=(atom_size * scene_scale))
    bpy.ops.object.shade_smooth()
    atom_obj = scene.objects.active
    atom_mesh = atom_obj.data
    atom_mesh.name = 'Atom_Template'
    scene.objects.unlink(atom_obj)
    bpy.data.objects.remove(atom_obj)

    layers = scene.layers

    # After parsing and preparing the templates, generate the output
    for model_id, model in model_list.items():
        if (not multi_models) and model_id != 1:
            break

        for bm_serial, biomolecule in biomolecule_list.items():
            biom_mesh_name = "Biomolecule%d.%d" % (bm_serial, model_id)
            biom_mesh = bpy.data.meshes.new(biom_mesh_name)
            cur_biom = bpy.data.objects.new(biom_mesh_name, biom_mesh)
            scene.objects.link(cur_biom)

            for chain_clones, transforms in biomolecule.chain_transforms.items():
                for chain in chain_clones:

                    for transform in transforms:
                        # rotations
                        row1 = transform[0][0:3]
                        row2 = transform[1][0:3]
                        row3 = transform[2][0:3]
                        # translations
                        dx = transform[0][3]
                        dy = transform[1][3]
                        dz = transform[2][3]

                        for serial, atom in model.chains[chain].atoms.items():
                            # Prunes alternative locations for the atoms
                            # (should pick the one with highest occupancy but doesn't)
                            if (atom.altloc == ' ' or atom.altloc == 'A') or retain_alts:

                                element = atom_data[atom.element].name

                                # Generate master element models
                                if element not in mat_list:

                                    # Create a master atom
                                    mesh = atom_mesh.copy()
                                    mesh.name = element
                                    mat = bpy.data.materials.new(element)
                                    mat.diffuse_color = atom_data[atom.element].color
                                    mesh.materials.append(mat)
                                    master_atom = bpy.data.objects.new(element, mesh)
                                    master_atom.scale = [atom_data[atom.element].vdw_radius] * 3
                                    master_atom.layers = layers
                                    master_atom.name = element  # why set again?
                                    scene.objects.link(master_atom)

                                    mat_list.append(element)
                                else:
                                    pass

                                # Generate element vertex groups
                                if element not in cur_biom.vertex_groups:
                                    cur_vert_group = cur_biom.vertex_groups.new(element)
                                    #Adds a key in the vert_group_count
                                    biomolecule.vert_group_counts[cur_vert_group] = 0

                                else:
                                    cur_vert_group = cur_biom.vertex_groups[element]

                                # Generate particle systems
                                # (can be merged with vertex group generator)
                                if element not in cur_biom.particle_systems:
                                    scene.objects.active = cur_biom
                                    bpy.ops.object.particle_system_add()

                                    psys = cur_biom.particle_systems.active
                                    part = psys.settings

                                    psys.name = element
                                    part.name = "%s.%d.%d" % (element, bm_serial, model_id)
                                    part.frame_start = 0
                                    part.frame_end = 0
                                    part.lifetime = 10000
                                    part.emit_from = 'VERT'
                                    part.use_emit_random = False
                                    part.normal_factor = 0
                                    part.particle_size = 1
                                    part.render_type = 'OBJECT'
                                    part.dupli_object = bpy.data.objects[element]
                                    part.effector_weights.gravity = 0
                                    psys.vertex_group_density = element

                                biomolecule.vert_group_index[biomolecule.atom_count] = cur_vert_group
                                biomolecule.vert_group_counts[cur_vert_group] += 1
                                biomolecule.atom_count += 1
                                tx = atom.x * row1[0] + atom.y * row1[1] + atom.z * row1[2] + (dx * scene_scale)
                                ty = atom.x * row2[0] + atom.y * row2[1] + atom.z * row2[2] + (dy * scene_scale)
                                tz = atom.x * row3[0] + atom.y * row3[1] + atom.z * row3[2] + (dz * scene_scale)
                                biomolecule.vert_list.extend([tx, ty, tz])
                            else:
                                pass

            biom_mesh.vertices.add(biomolecule.atom_count)
            biom_mesh.vertices.foreach_set('co', biomolecule.vert_list)

            for vert_index, vert_group in biomolecule.vert_group_index.items():
                aa = []  # neeed array to vertex_groups.assign
                aa.append(vert_index)
                vert_group.add(aa, 1, "ADD")

            for vert_group, count in biomolecule.vert_group_counts.items():
                bpy.data.particles["%s.%d.%d" % (vert_group.name, bm_serial, model_id)].count = count