18e46578097e6139204e60e84289782297bacf0f
[blender-addons-contrib.git] / io_import_pdb_atomic_blender.py
1 #
2 #
3 #  Authors           : Clemens Barth (Blendphys@root-1.de), ...
4 #
5 #  Homepage(Wiki)    : http://development.root-1.de/Atomic_Blender.php
6 #  Tracker           : http://projects.blender.org/tracker/index.php?func=detail&aid=29226&group_id=153&atid=467
7 #
8 #  Start of project              : 2011-08-31 by Clemens Barth
9 #  First publication in Blender  : 2011-11-11
10 #  Last modified                 : 2011-11-18
11 #
12 #
13 # ##### BEGIN GPL LICENSE BLOCK #####
14 #
15 #  This program is free software; you can redistribute it and/or
16 #  modify it under the terms of the GNU General Public License
17 #  as published by the Free Software Foundation; either version 2
18 #  of the License, or (at your option) any later version.
19 #
20 #  This program is distributed in the hope that it will be useful,
21 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
22 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
23 #  GNU General Public License for more details.
24 #
25 #  You should have received a copy of the GNU General Public License
26 #  along with this program; if not, write to the Free Software Foundation,
27 #  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 #
29 # ##### END GPL LICENSE BLOCK #####
30 #
31 #
32 #
33 #
34 #
35 #
36
37 bl_info = {
38   "name": "Atomic Blender",
39   "description": "Loading and manipulating atoms from PDB files",
40   "author": "Dr. Clemens Barth",
41   "version": (1,11),
42   "blender": (2,6),
43   "api": 31236,
44   "location": "Properties > Physics => Panel Atomic Blender",
45   "warning": "",
46   "wiki_url": "http://development.root-1.de/Atomic_Blender.php",
47   "tracker_url": "http://projects.blender.org/tracker/?func=detail&atid=467&aid=29226&group_id=153",
48   "category": "Import-Export"
49 }
50
51 import bpy
52 import io
53 import sys
54 import os
55 from math import *
56 import mathutils, math
57 from mathutils import Vector
58
59
60 # These are variables, which contain the name of the PDB file,
61 # the path of the PDB file and, finally, the path to the DATAFILEPATH.
62 # They are used almost everywhere, which is the reason why they 
63 # should stay global. First, they are empty and get 'filled' directly
64 # after having chosen the PDB file (see discussion at 'class LoadPDB'
65 # further below).
66 #
67 # Note that the treatment of the data filepath needs to be adjusted for 
68 # the case that the script shall be called during startup (see discussion 
69 # in the class LoadPDB further below).
70
71 PDBFILEPATH       = ""
72 PDBFILENAME       = ""
73 DATAFILEPATH      = ""
74
75 # The name of this script and the data file. This is used in the class 
76 # LoadPDB, for determining the path of the data file. For further details 
77 # see below.
78 SCRIPTNAME   = "io_import_pdb_atomic_blender.py"
79 DATAFILENAME = "io_import_pdb_atomic_blender.dat"
80
81
82 # Some string stuff for the console.
83 Atomic_Blender_string     = "Atomic Blender 1.11\n==================="
84 Atomic_Blender_panel_name = "PDB - Atomic Blender"
85
86
87
88
89
90 # The panel, which is loaded after the file has been
91 # chosen via the menu 'File -> Import'
92 class CLASS_PDB_Panel(bpy.types.Panel):
93     bl_label       = Atomic_Blender_panel_name
94     bl_space_type  = "PROPERTIES"
95     bl_region_type = "WINDOW"
96     bl_context     = "physics"
97     # This could be also an option ... :
98     #bl_space_type  = "VIEW_3D"
99     #bl_region_type = "TOOL_PROPS"
100
101     # This 'poll thing' has taken 3 hours of a hard search and understanding.
102     # I explain it in the following from my point of view:
103     #
104     # Before this class is entirely treaten (here: drawing the panel) the
105     # poll method is called first. Basically, some conditions are 
106     # checked before other things in the class are done afterwards. If a 
107     # condition is not valid, one returns 'False' such that nothing further 
108     # is done. 'True' means: 'Go on'
109     #
110     # In the case here, it is verified if the PDBFILEPATH variable contains
111     # a name. If not - and this is the case directly after having started the
112     # script - the panel does not appear because 'False' is returned. However,
113     # as soon as a file has been chosen, the panel appears because PDBFILEPATH
114     # contains a name.
115     #
116     # Please, correct me if I'm wrong. 
117     @classmethod
118     def poll(self, context):
119         if PDBFILEPATH == "":
120             return False
121         else:
122             return True
123
124     def draw(self, context):
125         layout = self.layout
126         scn    = bpy.context.scene
127
128         row = layout.row()
129         col = row.column(align=True)
130         col.prop(scn, "atom_pdb_PDB_filename") 
131         col.prop(scn, "atom_pdb_PDB_file")
132         row = layout.row()
133         row = layout.row() 
134   
135         row.prop(scn, "atom_pdb_group_atoms_yesno")
136         row.prop(scn, "atom_pdb_group_atoms_dn")
137
138         row = layout.row()
139         col = row.column(align=True)        
140         col.prop(scn, "atom_pdb_mesh_yesno")
141         col.prop(scn, "atom_pdb_mesh_azimuth")
142         col.prop(scn, "atom_pdb_mesh_zenith")        
143         col = row.column(align=True)        
144         col.label(text="Scaling factors")
145         col.prop(scn, "atom_pdb_scale_ballradius")
146         col.prop(scn, "atom_pdb_scale_distances")
147         col = row.column(align=True) 
148         col.prop(scn, "atom_pdb_sticks_yesno")
149         col.prop(scn, "atom_pdb_sticks_sectors")
150         col.prop(scn, "atom_pdb_sticks_radius")
151
152         row = layout.row()        
153         col = row.column(align=True)  
154         col.prop(scn, "atom_pdb_offset_x")
155         col.prop(scn, "atom_pdb_offset_y")
156         col.prop(scn, "atom_pdb_offset_z")
157         col = row.column()         
158         col.prop(scn, "atom_pdb_center_yesno")
159
160         layout.separator()  
161         row = layout.row(align=True)        
162         col = row.column()
163         col.prop(scn, "atom_pdb_cam_yesno")
164         col.prop(scn, "atom_pdb_lamp_yesno")        
165         col = row.column() 
166         col.operator( "atom_pdb.button_start" )
167         row2 = col.row()
168         row2.label(text="Number of atoms")
169         row2.prop(scn, "atom_pdb_start_number_atoms")
170         layout.separator()
171               
172         if scn.atom_pdb_group_atoms_yesno == False:        
173
174             row = layout.row()             
175             row.operator( "atom_pdb.button_distance")
176             row.prop(scn, "atom_pdb_distance") 
177             layout.separator()
178             
179             row = layout.row()                   
180             row.label(text="Modification of the radii of one type of atom")            
181             
182             row = layout.row()     
183             split = row.split(percentage=0.40)
184             col = split.column()
185             col.prop(scn, "atom_pdb_mod_atomname") 
186             split = split.split(percentage=0.50)
187             col = split.column()
188             col.prop(scn, "atom_pdb_mod_pm_radius")
189             split = split.split(percentage=1.0)
190             col = split.column()
191             col.operator("atom_pdb.button_modify_single")
192             
193             row = layout.row()     
194             split = row.split(percentage=0.40)
195             col = split.column()
196             split = split.split(percentage=0.50)
197             col = split.column()
198             col.prop(scn, "atom_pdb_mod_rel_radius")
199             split = split.split(percentage=1.0)
200             col = split.column(align=True)
201             col.operator( "atom_pdb.button_bigger_single" )            
202             col.operator( "atom_pdb.button_smaller_single" ) 
203                                              
204             row = layout.row()            
205             row.label(text="Modification of all atom radii")
206             row = layout.row()
207             col = row.column() 
208             col.prop(scn, "atom_pdb_mod_all_radii")
209             col = row.column(align=True) 
210             col.operator( "atom_pdb.button_modify_all" )
211             col.operator( "atom_pdb.button_invert_all" )
212
213
214 class CLASS_Input_Output(bpy.types.PropertyGroup):
215     bpy.types.Scene.atom_pdb_PDB_filename        = bpy.props.StringProperty(name = "File name", default="", description = "PDB file name")
216     bpy.types.Scene.atom_pdb_PDB_file            = bpy.props.StringProperty(name = "Path to file", default="", description = "Path of the PDB file")
217     bpy.types.Scene.atom_pdb_group_atoms_yesno   = bpy.props.BoolProperty  (name = "Group atoms", default=False, description = "Grouping same type of atoms speeds up the loading of large-atom PDB files")    
218     bpy.types.Scene.atom_pdb_group_atoms_dn      = bpy.props.IntProperty   (name = "Delta N", default=200, min=0, description = "The number of single atoms, which are grouped together after their loading")
219     bpy.types.Scene.atom_pdb_mesh_yesno          = bpy.props.BoolProperty  (name = "Mesh balls", default=False, description = "Do you want to use mesh balls instead of NURBS?")    
220     bpy.types.Scene.atom_pdb_mesh_azimuth     = bpy.props.IntProperty   (name = "Azimuth", default=32, min=0, description = "Number of sectors (azimuth)")
221     bpy.types.Scene.atom_pdb_mesh_zenith      = bpy.props.IntProperty   (name = "Zenith", default=32, min=0, description = "Number of sectors (zenith)")
222     bpy.types.Scene.atom_pdb_scale_ballradius         = bpy.props.FloatProperty (name = "Balls", default=1.0, min=0.0, description = "Scale factor for all atom radii")
223     bpy.types.Scene.atom_pdb_scale_distances           = bpy.props.FloatProperty (name = "Distances", default=1.0, min=0.0, description = "Scale factor for all distances")
224     bpy.types.Scene.atom_pdb_center_yesno        = bpy.props.BoolProperty  (name = "Object to origin", default=True, description = "Shall the object first put into the global origin before applying the offsets on the left?")    
225     bpy.types.Scene.atom_pdb_offset_x             = bpy.props.FloatProperty (name = "X", default=0.0, description = "Offset in X")
226     bpy.types.Scene.atom_pdb_offset_y             = bpy.props.FloatProperty (name = "Y", default=0.0, description = "Offset in Y")
227     bpy.types.Scene.atom_pdb_offset_z             = bpy.props.FloatProperty (name = "Z", default=0.0, description = "Offset in Z")
228     bpy.types.Scene.atom_pdb_sticks_yesno        = bpy.props.BoolProperty  (name = "Use sticks", default=False, description = "Do you want to display also the sticks?")    
229     bpy.types.Scene.atom_pdb_sticks_sectors      = bpy.props.IntProperty   (name = "Sector", default = 20, min=0,   description = "Number of sectors of a stick")        
230     bpy.types.Scene.atom_pdb_sticks_radius       = bpy.props.FloatProperty (name = "Radius", default =  0.1, min=0.0, description = "Radius of a stick")  
231     bpy.types.Scene.atom_pdb_cam_yesno           = bpy.props.BoolProperty  (name = "Camera", default=False, description = "Do you need a camera?")   
232     bpy.types.Scene.atom_pdb_lamp_yesno          = bpy.props.BoolProperty  (name = "Lamp", default=False, description = "Do you need a lamp?")
233     bpy.types.Scene.atom_pdb_start_number_atoms  = bpy.props.StringProperty(name = "", default="Number", description = "This output shows the number of atoms which have been loaded")
234     bpy.types.Scene.atom_pdb_distance            = bpy.props.StringProperty(name = "", default="Distance (Angstrom)", description = "Distance of 2 objects in Angstrom")  
235     bpy.types.Scene.atom_pdb_mod_atomname        = bpy.props.StringProperty(name = "", default = "Atom name", description="Put in the name of the atom (e.g. Hydrogen)")
236     bpy.types.Scene.atom_pdb_mod_pm_radius       = bpy.props.FloatProperty (name = "", default = 100.0, min=0.0, description="Put in the radius of the atom (in pm)")
237     bpy.types.Scene.atom_pdb_mod_rel_radius      = bpy.props.FloatProperty (name = "", default = 1.05, min=1.0, description="Put in the scale factor")
238     bpy.types.Scene.atom_pdb_mod_all_radii       = bpy.props.FloatProperty (name = "Scale", default = 1.05, min=1.0, description="Put in the scale factor")
239
240
241 # Button for measuring the distance of the active objects
242 class CLASS_Distance_Button(bpy.types.Operator):
243     bl_idname = "atom_pdb.button_distance"
244     bl_label = "Measure ..."
245     bl_description = "Measure the distance between two objects."
246
247     def execute(self, context):
248         dist   = Measure_distance_in_scene()
249
250         if dist != "-1.0":
251            # The string length is cut, 3 digits after the first 3 digits 
252            # after the '.'. Append also "Angstrom". 
253            # Remember: 1 Angstrom = 10^(-10) m 
254            pos    = str.find(dist, ".")
255            dist   = dist[:pos+4] 
256            dist   = dist + " Angstrom"
257
258         # Put the distance into the string of the output field.
259         scn                = bpy.context.scene
260         scn.atom_pdb_distance = dist
261         return {'FINISHED'}
262   
263
264 # Button for changing the radii (in pm) of atoms of one type
265 class CLASS_Modify_Single_Button(bpy.types.Operator):
266     bl_idname = "atom_pdb.button_modify_single"
267     bl_label = "Modify ..."
268     bl_description = "Change the radii of atoms of one type in pm."
269
270     def execute(self, context):
271         scn = bpy.context.scene
272         atomname   = scn.atom_pdb_mod_atomname
273         radius_pm  = scn.atom_pdb_mod_pm_radius
274         Modify_atom_radii_type_pm(atomname, radius_pm)
275         return {'FINISHED'}
276
277
278 # Button for increasing the radii of atoms of one type
279 class CLASS_Bigger_Single_Button(bpy.types.Operator):
280     bl_idname = "atom_pdb.button_bigger_single"
281     bl_label = "Bigger ..."
282     bl_description = "Increase the radii of atoms of one type."
283
284     def execute(self, context):
285         scn = bpy.context.scene
286         atomname   = scn.atom_pdb_mod_atomname
287         radius_rel = scn.atom_pdb_mod_rel_radius
288         Modify_atom_radii_type_scale(atomname, radius_rel)
289         return {'FINISHED'}
290
291
292 # Button for decreasing the radii of atoms of one type
293 class CLASS_Smaller_Single_Button(bpy.types.Operator):
294     bl_idname = "atom_pdb.button_smaller_single"
295     bl_label = "Smaller ..."
296     bl_description = "Decrease the radii of atoms of one type."
297
298     def execute(self, context):
299         scn = bpy.context.scene
300         atomname   = scn.atom_pdb_mod_atomname
301         radius_rel = 1.0/scn.atom_pdb_mod_rel_radius
302         Modify_atom_radii_type_scale(atomname, radius_rel)
303         return {'FINISHED'}
304
305
306 # Button for increasing the radii of all atoms
307 class CLASS_Bigger_All_Button(bpy.types.Operator):
308     bl_idname = "atom_pdb.button_modify_all"
309     bl_label = "Bigger ..."
310     bl_description = "Increase the radii of all atoms."
311
312     def execute(self, context):
313         scn     = bpy.context.scene
314         scale   = scn.atom_pdb_mod_all_radii
315         Modify_all_atom_radii(scale)
316         return {'FINISHED'}
317
318
319 # Button for decreasing the radii of all atoms
320 class CLASS_Smaller_All_Button(bpy.types.Operator):
321     bl_idname = "atom_pdb.button_invert_all"
322     bl_label = "Smaller ..."
323     bl_description = "Decrease the radii of all atoms."
324
325     def execute(self, context):
326         scn     = bpy.context.scene
327         scale   = 1.0/scn.atom_pdb_mod_all_radii
328         Modify_all_atom_radii(scale)
329         return {'FINISHED'}
330
331
332 # The button for loading the atoms and creating the scene
333 class CLASS_Start_Button(bpy.types.Operator):
334     bl_idname = "atom_pdb.button_start"
335     bl_label = "DRAW THE ATOMS ..."
336     bl_description = "Start to load and draw the atoms and sticks."
337     
338     def execute(self, context):
339         scn = bpy.context.scene
340
341         azimuth   = scn.atom_pdb_mesh_azimuth
342         zenith    = scn.atom_pdb_mesh_zenith 
343         bradius   = scn.atom_pdb_scale_ballradius
344         bdistance = scn.atom_pdb_scale_distances
345         center    = scn.atom_pdb_center_yesno 
346         x         = scn.atom_pdb_offset_x
347         y         = scn.atom_pdb_offset_y
348         z         = scn.atom_pdb_offset_z
349         yn        = scn.atom_pdb_sticks_yesno 
350         ssector   = scn.atom_pdb_sticks_sectors
351         sradius   = scn.atom_pdb_sticks_radius
352         pdb       = PDBFILEPATH
353         data      = DATAFILEPATH
354         cam       = scn.atom_pdb_cam_yesno 
355         lamp      = scn.atom_pdb_lamp_yesno
356         mesh      = scn.atom_pdb_mesh_yesno 
357         group     = scn.atom_pdb_group_atoms_yesno
358         dn        = scn.atom_pdb_group_atoms_dn
359         
360         atom_number                  = Draw_scene(group,dn,mesh,azimuth,zenith,bradius,bdistance,x,y,z,yn,ssector,sradius,center,cam,lamp)
361         scn.atom_pdb_start_number_atoms = atom_number
362
363         return {'FINISHED'}
364
365
366 # This is the class for the file dialog.
367 class CLASS_LoadPDB(bpy.types.Operator):
368     bl_idname = "import_image.brushset"
369     bl_label  = "Import PDB"
370
371     filename = bpy.props.StringProperty(name = "PDB  File", description="Path to the PDB file", maxlen = 256, default = PDBFILEPATH, subtype='FILE_PATH', options={'HIDDEN'})
372     filepath = bpy.props.StringProperty(name = "PDB  File", description="Path to the PDB file", maxlen = 256, default = PDBFILEPATH, subtype='FILE_PATH', options={'HIDDEN'})
373     
374     def execute(self, context):
375         global PDBFILEPATH
376         global PDBFILENAME
377         global DATAFILEPATH
378         
379         # In the following the name and path of the PDB file
380         # is stored into the global variables.
381         scn = bpy.context.scene
382         PDBFILENAME     = self.filename
383         PDBFILEPATH     = self.filepath
384         scn.atom_pdb_PDB_filename = PDBFILENAME
385         scn.atom_pdb_PDB_file     = PDBFILEPATH
386         
387         # This needs to be changed in future, if one decides to permanently
388         # include the PDB importer.
389         # So far, it takes the path of the loaded Atomic Blender script
390         # and replaces the name of the python script with the name of the
391         # data file.
392         # This works, so far, only if the script is loaded manually. 
393         if bpy.data.texts[0].filepath != "":
394             DATAFILEPATH = bpy.data.texts[0].filepath
395         else:
396             DATAFILEPATH = bpy.data.texts[-1].filepath
397         DATAFILEPATH = DATAFILEPATH.replace(SCRIPTNAME, DATAFILENAME)   
398         print(DATAFILEPATH)
399
400         return {'FINISHED'}
401
402     def invoke(self, context, event):
403         wm = context.window_manager
404         wm.fileselect_add(self)
405         return {'RUNNING_MODAL'}
406
407
408 # The entry into the menu 'file -> import'
409 def menu_func(self, context):
410     default_name = "" 
411     self.layout.operator(CLASS_LoadPDB.bl_idname, text="PDB (.pdb)").filename = default_name
412
413
414 def register():
415     bpy.utils.register_class(CLASS_PDB_Panel)
416     bpy.utils.register_class(CLASS_Input_Output)
417     bpy.utils.register_class(CLASS_Start_Button)
418     bpy.utils.register_class(CLASS_Modify_Single_Button)
419     bpy.utils.register_class(CLASS_Bigger_Single_Button)
420     bpy.utils.register_class(CLASS_Smaller_Single_Button)
421     bpy.utils.register_class(CLASS_Bigger_All_Button)
422     bpy.utils.register_class(CLASS_Smaller_All_Button)
423     bpy.utils.register_class(CLASS_Distance_Button)
424     bpy.utils.register_module(__name__)
425     bpy.types.INFO_MT_file_import.append(menu_func)
426
427
428 def unregister():
429     bpy.utils.unregister_class(CLASS_PDB_Panel)
430     bpy.utils.unregister_class(CLASS_Input_Output)
431     bpy.utils.unregister_class(CLASS_Start_Button)
432     bpy.utils.unregister_class(CLASS_Modify_Single_Button)
433     bpy.utils.unregister_class(CLASS_Bigger_Single_Button)
434     bpy.utils.unregister_class(CLASS_Smaller_Single_Button)    
435     bpy.utils.unregister_class(CLASS_Bigger_All_Button)
436     bpy.utils.unregister_class(CLASS_Smaller_All_Button)
437     bpy.utils.unregister_class(CLASS_Distance_Button)  
438     bpy.utils.unregister_module(__name__)  
439     bpy.types.INFO_MT_file_import.remove(menu_func)
440         
441
442 if __name__ == "__main__":
443
444     register()
445
446
447
448
449
450
451
452
453
454
455
456
457 ########################################################
458 #
459 #
460 #
461 #
462 #
463 #
464 #          Some small routines
465 #
466 #
467 #
468 #
469 #
470 ########################################################
471
472
473
474 # This function measures the distance between two objects (atoms), which are marked.
475 def Measure_distance_in_scene():
476
477     if len(bpy.context.selected_bases) > 1:
478         object_1 = bpy.context.selected_objects[0]
479         object_2 = bpy.context.selected_objects[1]
480     else:
481         return "-1.0"
482
483     v1     = object_1.location
484     v2     = object_2.location
485     dv     = (v2 - v1)
486     length = str(dv.length)
487     return length 
488
489
490 # Routine to modify the radii of a specific type of atom
491 def Modify_atom_radii_type_pm(atomname, radius_pm):
492
493     for obj in bpy.data.objects:
494
495         if atomname in obj.name:
496
497             bpy.data.objects[obj.name].scale = (radius_pm/100,radius_pm/100,radius_pm/100)
498                 
499
500 # Routine to modify the radii of a specific type of atom
501 def Modify_atom_radii_type_scale(atomname, radius_rel):
502
503     for obj in bpy.data.objects:
504
505         if atomname in obj.name:
506
507             value = obj.scale[0]
508             bpy.data.objects[obj.name].scale = (radius_rel * value,radius_rel * value,radius_rel * value)
509
510
511 # Routine to scale the radii of all atoms
512 def Modify_all_atom_radii(scale):
513
514     for obj in bpy.data.objects:
515
516         if obj.type == "SURFACE" or obj.type == "MESH":
517
518             if "Stick" not in obj.name:
519                 radius = obj.scale[0]
520                 obj.scale = (radius * scale,radius * scale,radius * scale)
521
522
523
524 ########################################################
525 #
526 #
527 #
528 #
529 #
530 #          For reading the sticks inside the PDB file
531 #
532 #
533 #
534 #
535 #
536 ########################################################
537
538
539
540 def Read_atom_for_stick(string):        
541
542     string_length = len(string)
543
544     j               = 0
545     string_reversed = ""
546     atoms           = []
547     space           = False
548     # An atom number can have max 5 letters! This automatically means that
549     # up to 99999 atoms can be loaded max. - Well, this should be sufficient.
550     counter_letters = 5 
551    
552     # I know that what follows is somewhat 'confusing'! However, the strings of
553     # the atom numbers do have a clear position in the file (From 1 to 5, 
554     # from 6 to 10 and so on.) and one needs to consider this. One could also use
555     # the split function but then one gets into trouble if there are many atoms:
556     # For instance, it may happen that one has
557     #
558     # CONECT 11111  22244444
559     #
560     # In Fact it means that atom No. 11111 has a stick with atom No. 222 but also
561     # with atom No. 44444. The split function would give me only two numbers (11111
562     # and 22244444), which is wrong. However, the following code supplies 
563     # the three correct numbers: 
564     for i in list(range(string_length)):
565    
566         # If the 'T' of 'CONECT' is read => exit
567         if string[string_length-i-1] == 'T':
568             break
569
570         # Continue, if a space is read but no letter is present in 'string_reversed'.
571         # This happens, when there are spaces behind the last atom number in the
572         # string. 
573         if string[string_length-i-1] == ' ' and string_reversed == "":
574             continue
575    
576         if string[string_length-i-1] == ' ' or counter_letters == 0:
577       
578             string_correct         = ""
579             string_reversed_length = len(string_reversed)
580             l                      = 0
581             for k in list(range(string_reversed_length)):
582                 string_correct = string_correct + string_reversed[string_reversed_length-k-1]
583                 l += 1
584       
585             # If the first 'space' is found, we found the number of an atom
586             # Transform the string into an integer and append this to the overall list
587             if space == False:
588                 atom            = int(string_correct)
589                 atoms.append(atom)
590                 # Initialization of the variables
591                 string_reversed = ""
592                 space           = True
593             
594             
595             # If it was only a 'space' then go up the 'for loop'.
596             if counter_letters != 0:
597                 counter_letters  = 5
598                 continue
599             
600             counter_letters = 5
601    
602    
603         space            = False
604         string_reversed  = string_reversed + string[string_length-i-1]
605         j               += 1
606         # One letter has been read, so one down with the counter. 
607         # Max is 5! 
608         counter_letters -= 1
609       
610     # Return the list of atoms   
611     return atoms
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642 ########################################################
643 #
644 #
645 #
646 #
647 #
648 #          The main routine
649 #
650 #
651 #
652 #
653 #
654 ########################################################
655
656
657
658 def Draw_scene(FLAG_group_atoms,group_atoms_dn,mesh_yn,Ball_azimuth,Ball_zenith,Ball_radius_factor,Ball_distance_factor,offset_x,offset_y,offset_z,Stick_yn,Stick_sectors,Stick_diameter, put_to_center, camera_yn, lamp_yn):
659
660     global PDBFILEPATH
661     global PDBFILENAME
662     global DATAFILEPATH
663     
664
665     # This is in order to solve this strange 'relative path' thing.
666     PDBFILEPATH  = bpy.path.abspath(PDBFILEPATH)
667     DATAFILEPATH = bpy.path.abspath(DATAFILEPATH)
668    
669    
670     # Lists for all atoms in the data file
671     Data_Number        = []
672     Data_Atomname      = []
673     Data_Shortname     = []
674     Data_Color         = []
675     Data_Radius        = []
676
677     # Lists for atoms in the PDB file
678     atom_object    = []
679     atom_element   = []
680     atom_name      = []
681     atom_charge    = []
682     atom_color     = []
683     atom_material  = []   
684     atom_x         = []
685     atom_y         = []
686     atom_z         = []
687     atom_R         = []
688     bar_atom1      = []
689     bar_atom2      = []
690    
691     # Materials
692     atom_material_list = []
693    
694
695     #
696     #
697     #
698     #
699     #          READING DATA OF ALL POSSIBLE ATOMS FROM THE DATA FILE
700     #         
701     #
702     #
703     #
704
705
706     # Read the data file, which contains all data (atom name, radii, colors, etc.)
707     DATAFILEPATH_p = io.open(DATAFILEPATH, "r")
708
709     i = 0
710     for line in DATAFILEPATH_p:
711
712       if "Atom" in line:
713
714           line              = DATAFILEPATH_p.readline() 
715           line              = DATAFILEPATH_p.readline()
716           pos               = str.find(line, ":")
717           Data_Number.append(line[pos+2:-1])
718
719           line              = DATAFILEPATH_p.readline()
720           pos               = str.find(line, ":")
721           Data_Atomname.append(line[pos+2:-1])
722
723           line              = DATAFILEPATH_p.readline()
724           pos               = str.find(line, ":")
725           Data_Shortname.append(line[pos+2:-1])
726
727           line              = DATAFILEPATH_p.readline()
728           pos               = str.find(line, ":")
729           color_value       = line[pos+2:-1].split(',')
730           Data_Color.append([float(color_value[0]),float(color_value[1]),float(color_value[2])]) 
731
732           line              = DATAFILEPATH_p.readline()
733           pos               = str.find(line, ":")
734           Data_Radius.append(line[pos+2:-1])
735
736           i += 1
737
738     DATAFILEPATH_p.close()
739     
740     all_existing_atoms = i
741
742
743     #
744     #
745     #
746     #
747     #          READING DATA OF ATOMS
748     #
749     #
750     #
751
752
753     # Open the file ...
754     PDBFILEPATH_p = io.open(PDBFILEPATH, "r")
755
756     #Go to the line, in which "ATOM" or "HETATM" appears.
757     for line in PDBFILEPATH_p:
758         split_list = line.split(' ')
759         if "ATOM" in split_list[0]:
760             break
761         if "HETATM" in split_list[0]:
762             break
763  
764     # This is the list, which contains the names of all type of atoms.
765     atom_all_types_list = []
766
767     j = 0
768     # This is in fact an endless 'while loop', ...
769     while j > -1:
770
771         # ... the loop is broken here (EOF) ...
772         if line == "":
773             break  
774
775         # If 'ATOM4 or 'HETATM' appears in the line then do ...
776         if "ATOM" in line or "HETATM" in line:
777         
778             # Split the line into its parts (devided by a ' ') and analyse it. The first line is read.
779             split_list = line.rsplit()
780          
781             for i in list(range(all_existing_atoms)):
782                 if str.upper(split_list[-1]) == str.upper(Data_Shortname[i]):
783                     # Give the atom its proper name and radius:
784                     atom_element.append(str.upper(Data_Shortname[i]))
785                     atom_name.append(Data_Atomname[i])
786                     atom_R.append(float(Data_Radius[i]))
787                     atom_color.append(Data_Color[i])
788                     break
789
790             # 1. case: These are 'unknown' atoms. In some cases, atoms are named with an additional label like H1 (hydrogen1)
791             # 2. case: The last column 'split_list[-1]' does not exist, we take then column 3 in the PDB file.
792             if i == all_existing_atoms-1:
793
794                 # Give this atom also a name. If it is an 'X' then it is a vacancy. Otherwise ...
795                 if "X" in str.upper(split_list[2]):
796                     atom_element.append("VAC")
797                     atom_name.append("Vacancy")
798                 # ... take what is written in the PDB file.
799                 else:
800                     atom_element.append(str.upper(split_list[2]))
801                     atom_name.append(str.upper(split_list[2]))
802
803                 # Default values for the atom.
804                 atom_R.append(float(Data_Radius[all_existing_atoms-2]))
805                 atom_color.append(Data_Color[all_existing_atoms-2])
806          
807                  
808          
809             # The list that contains info about all types of atoms is created here.
810             # It is used for building the material properties for instance. 
811             
812             # If the name of the atom is already in the list, FLAG on 'True'. 
813             FLAG_FOUND = False
814             for atom_type in atom_all_types_list:
815                 if atom_type[0] == atom_name[-1]:
816                     FLAG_FOUND = True
817                     break
818             # No name in the current list has been found? => New entry.
819             if FLAG_FOUND == False:
820                 # Stored are: Atom label (e.g. 'Na'), the corresponding atom name (e.g. 'Sodium') and its color.
821                 atom_all_types_list.append([atom_name[-1],atom_element[-1],atom_color[-1]])
822
823
824                  
825             # Now the coordinates x, y and z are read.
826             coor   = 1
827             number = 0
828          
829             # The coordinates x, y and z are identified as such by the dot (e.g., 5.678) - a dot
830             # in the line appears only in the coordinates. The coordinates are listed as follwos: 
831             # x, y and z.
832             # If the first coordinate (x) is read, increase coor (y).
833             for each_element in split_list:
834     
835                 # If there is a dot, it is an coordinate.
836                 if "." in each_element:
837                     if coor == 1:
838                         atom_x.append(float(each_element))
839                         coor     += 1
840                     elif coor == 2:
841                         atom_y.append(float(each_element))
842                         coor     += 1
843                     elif coor == 3:
844                         atom_z.append(float(each_element))
845                         coor     += 1        
846       
847             j += 1
848            
849                
850         line = PDBFILEPATH_p.readline()
851         line = line[:-1]
852
853     PDBFILEPATH_p.close()
854     # From above it can be clearly seen that j is now the number of all atoms.
855     Number_of_total_atoms = j
856
857
858     #
859     #
860     #
861     #
862     #          MATERIAL PROPERTIES FOR ATOMS
863     #
864     #
865     #
866
867
868
869     # Here, the atoms get already their material properties. Why already here? Because
870     # then it is done and the atoms can be drawn in a fast way (see drawing part at the end 
871     # of this script, further below). 
872     # Note that all atoms of one type (e.g. all hydrogens) get only ONE material! This 
873     # is good because then, by activating one atom in the Blender scene and changing 
874     # the color of this atom, one changes the color of ALL atoms of the same type at the 
875     # same time.
876    
877     # Create first a new list of materials for each type of atom (e.g. hydrogen)
878     for atom_type in atom_all_types_list:
879    
880         bpy.ops.object.material_slot_add()
881         material               = bpy.data.materials.new(atom_type[1])
882         material.name          = atom_type[0]
883         material.diffuse_color = atom_type[2]
884         atom_material_list.append(material)
885    
886     # Now, we go through all atoms and give them a material. For all atoms ...   
887     for i in range(0, Number_of_total_atoms):
888         # ... and all materials ...
889         for material in atom_material_list:
890             # ... select the correct material for the current atom via name-comparison ...
891             if atom_name[i] in material.name:
892                 # ... and give the atom its material properties. 
893                 # However, before we check, if it is a vacancy, because then it
894                 # gets some additional preparation. The vacancy is represented by
895                 # a transparent cube.
896                 if atom_name[i] == "Vacancy":
897                     material.transparency_method                  = 'Z_TRANSPARENCY'
898                     material.alpha                                = 1.3
899                     material.raytrace_transparency.fresnel        = 1.6
900                     material.raytrace_transparency.fresnel_factor = 1.6                   
901                     material.use_transparency                     = True      
902                 # The atom gets its properties.
903                 atom_material.append(material)   
904
905
906     #
907     #
908     #
909     #
910     #          READING DATA OF STICKS
911     #
912     #
913     #
914
915
916     # Open the PDB file again such that the file pointer is in the first line ...
917     # Stupid, I know ... ;-)
918     PDBFILEPATH_p = io.open(PDBFILEPATH, "r")
919
920     split_list = line.split(' ')
921
922     # Go to the first entry
923     if "CONECT" not in split_list[0]:
924         for line in PDBFILEPATH_p:
925             split_list = line.split(' ')
926             if "CONECT" in split_list[0]:
927                 break
928
929   
930     Number_of_sticks = 0
931     doppelte_bars  = 0
932     j              = 0
933     # This is in fact an endless while loop, ...    
934     while j > -1:
935  
936         # ... which is broken here (EOF) ...
937         if line == "":
938             break  
939         # ... or here, when no 'CONECT' appears anymore.
940         if "CONECT" not in line:
941             break
942          
943         if line[-1] == '\n':
944             line = line[:-1]
945         
946         # Read the sticks for the actual atom (sub routine). One gets a list of
947         # sticks.
948         atoms_list        = Read_atom_for_stick(line)
949         # Determine the length of the list
950         atoms_list_length = len(atoms_list)
951
952         # For all sticks in the list do:
953         q = 0
954         for each_element in atoms_list:
955       
956             # End == break
957             if q == atoms_list_length - 1:
958                 break
959       
960             # The first atom is connected with all the others in the list.
961             atom1 = atoms_list[-1]
962             # The second, third, ... partner atom
963             atom2 = each_element
964
965             FLAG_BAR = False
966  
967             # Note that in a PDB file, sticks of one atom pair can appear a couple
968             # of times. (Only god knows why ...) 
969             # So, does a stick between the considered atoms already exist?
970             for k in list(range(j)):
971                 if (bar_atom1[k] == atom1 and bar_atom2[k] == atom2) or (bar_atom1[k] == atom2 and bar_atom2[k] == atom1):
972                     doppelte_bars += 1
973                     # If yes, then FLAG on 'True'.
974                     FLAG_BAR       = True
975                     break
976
977             # If the stick is not yet registered (FLAG_BAR == False), then 
978             # register it!
979             if FLAG_BAR == False:
980                 bar_atom1.append(atom1)
981                 bar_atom2.append(atom2)      
982                 Number_of_sticks += 1   
983                 j += 1
984  
985             q += 1
986
987         line = PDBFILEPATH_p.readline()
988         line = line[:-1]
989
990     PDBFILEPATH_p.close()
991     # So far, all atoms and sticks have been registered.
992
993
994
995     #
996     #
997     #
998     #
999     #          TRANSLATION OF THE OBJECT TO THE ORIGIN
1000     #
1001     #
1002     #
1003
1004
1005
1006     # If chosen, the objects are first put into the center of the scene.
1007     if put_to_center == True:
1008
1009         sum_x = 0
1010         sum_y = 0
1011         sum_z = 0
1012
1013         # Sum of all atom coordinates
1014         for i in list(range(Number_of_total_atoms)):
1015
1016             sum_x = sum_x + atom_x[i]
1017             sum_y = sum_y + atom_y[i]
1018             sum_z = sum_z + atom_z[i]
1019
1020         # Then the average is taken
1021         sum_x = sum_x / Number_of_total_atoms
1022         sum_y = sum_y / Number_of_total_atoms
1023         sum_z = sum_z / Number_of_total_atoms
1024
1025         # After, for each atom the center of gravity is substracted
1026         for i in list(range(Number_of_total_atoms)):
1027
1028             atom_x[i] = atom_x[i] - sum_x
1029             atom_y[i] = atom_y[i] - sum_y
1030             atom_z[i] = atom_z[i] - sum_z
1031
1032
1033
1034     #
1035     #
1036     #
1037     #
1038     #          SCALING GEOMETRIC PROPERTIES
1039     #
1040     #
1041     #
1042
1043
1044     # Take all atoms and ...
1045     # - adjust their radii,
1046     # - scale the distances,
1047     # - and move the center of the whole ('+= offset_x', in Angstroem)
1048     for i in list(range(Number_of_total_atoms)):
1049
1050         atom_charge.append(1.0)  
1051         atom_x[i] += offset_x
1052         atom_y[i] += offset_y
1053         atom_z[i] += offset_z
1054         atom_x[i] *= Ball_distance_factor
1055         atom_y[i] *= Ball_distance_factor
1056         atom_z[i] *= Ball_distance_factor
1057
1058
1059
1060     #
1061     #
1062     #
1063     #
1064     #          DETERMINATION OF SOME GEOMETRIC PROPERTIES
1065     #
1066     #
1067     #
1068
1069
1070     # In the following, some geometric properties of the whole object are 
1071     # determined: center, size, etc. 
1072     sum_x = 0
1073     sum_y = 0
1074     sum_z = 0
1075
1076     # First the center is determined. All coordinates are summed up ...
1077     for i in list(range(Number_of_total_atoms)):
1078         sum_x = sum_x + atom_x[i]
1079         sum_y = sum_y + atom_y[i]
1080         sum_z = sum_z + atom_z[i]
1081     # ... and the average is taken. This gives the center of the object.
1082     object_center = [sum_x / Number_of_total_atoms, sum_y / Number_of_total_atoms, sum_z / Number_of_total_atoms]
1083
1084     # Now, we determine the size. All coordinates are analyzed ...
1085     object_size = 0.0
1086     for i in list(range(Number_of_total_atoms)):
1087
1088         diff_x = atom_x[i] - object_center[0]
1089         diff_y = atom_y[i] - object_center[1]
1090         diff_z = atom_z[i] - object_center[2]
1091
1092         # This is needed in order to estimate the size of the object.
1093         # The farest atom from the object center is taken as a measure.
1094         distance_to_object_center = math.sqrt(diff_x*diff_x + diff_y*diff_y + diff_z*diff_z)
1095         if distance_to_object_center > object_size:
1096             object_size = distance_to_object_center
1097
1098
1099     #
1100     #
1101     #
1102     #
1103     #          CAMERA AND LAMP
1104     #
1105     #
1106     #
1107
1108     camera_factor = 15.0
1109
1110     # Here a camera is put into the scene, if chosen.
1111     if camera_yn == True:
1112
1113         # Assume that the object is put into the global origin. Then, the camera
1114         # is moved in x and z direction, not in y. The object has its size at distance
1115         # math.sqrt(object_size) from the origin. So, move the camera by this distance
1116         # times a factor of camera_factor in x and z. Then add x, y and z of the origin of the
1117         # object.   
1118         camera_x = object_center[0] + math.sqrt(object_size) * camera_factor
1119         camera_y = object_center[1]
1120         camera_z = object_center[2] + math.sqrt(object_size) * camera_factor
1121         camera_pos    = [camera_x,camera_y,camera_z]
1122         # Create the camera
1123         bpy.ops.object.camera_add(view_align=False, enter_editmode=False, location=camera_pos, rotation=(0.0, 0.0, 0.0), layers=(True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False))
1124         # Some properties of the camera are changed.
1125         camera               = bpy.context.scene.objects.active
1126         camera.name          = "A_camera"
1127         camera.data.name     = "A_camera"
1128         camera.data.lens     = 45
1129         camera.data.clip_end = 500.0
1130
1131         # Here the camera is rotated such it looks towards the center of the object.
1132         
1133         # The vector between camera and origin of the object
1134         vec_cam_obj            = (mathutils.Vector(camera_pos) - mathutils.Vector(object_center))
1135         # The [0.0, 0.0, 1.0] vector along the z axis
1136         vec_up_axis            = mathutils.Vector([0.0, 0.0, 1.0])
1137         # The angle between the last two vectors
1138         angle                  = vec_cam_obj.angle(vec_up_axis, 0)
1139         # The cross-product of the [0.0, 0.0, 1.0] vector and vec_cam_obj
1140         # It is the resulting vector which stands up perpendicular on vec_up_axis and vec_cam_obj
1141         axis                   = vec_up_axis.cross(vec_cam_obj)
1142         # Rotate axis 'axis' by angle 'angle' and convert this to euler parameters. 4 is the size
1143         # of the matrix.
1144         euler                  = mathutils.Matrix.Rotation(angle, 4, axis).to_euler()
1145         camera.rotation_euler  = euler
1146
1147         # Rotate the camera around its axis by 90° such that we have a nice camera position
1148         # and view onto the object.
1149         bpy.ops.transform.rotate(value=(90.0*2*math.pi/360.0,), axis=vec_cam_obj, constraint_axis=(False, False, False), constraint_orientation='GLOBAL', mirror=False, proportional='DISABLED', proportional_edit_falloff='SMOOTH', proportional_size=1, snap=False, snap_target='CLOSEST', snap_point=(0, 0, 0), snap_align=False, snap_normal=(0, 0, 0), release_confirm=False)
1150
1151
1152         # This does not work, I don't know why. 
1153         #
1154         #for area in bpy.context.screen.areas:
1155         #    if area.type == 'VIEW_3D':
1156         #        area.spaces[0].region_3d.view_perspective = 'CAMERA'
1157
1158
1159     # Here a lamp is put into the scene, if chosen.
1160     if lamp_yn == True:
1161
1162
1163         # This is the distance from the object measured in terms of % of the camera 
1164         # distance. It is set onto 50% (1/2) distance.
1165         lamp_dl                        = math.sqrt(object_size) * 15 * 0.5
1166         # This is a factor to which extend the lamp shall go to the right (from the camera 
1167         # point of view).
1168         lamp_dy_right                  = lamp_dl * (3.0/4.0)
1169         
1170         # Create x, y and z for the lamp.
1171         lamp_x                         = object_center[0] + lamp_dl
1172         lamp_y                         = object_center[1] + lamp_dy_right
1173         lamp_z                         = object_center[2] + lamp_dl
1174         lamp_pos                       = [lamp_x, lamp_y, lamp_z]
1175         # Create the lamp
1176         bpy.ops.object.lamp_add  (type = 'POINT', view_align=False,         location=lamp_pos,   rotation=(0.0, 0.0, 0.0), layers=(True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False))
1177         # Some properties of the lamp are changed.
1178         lamp                           = bpy.context.scene.objects.active
1179         lamp.data.name                 = "A_lamp"
1180         lamp.name                      = "A_lamp"
1181         lamp.data.distance             = 500.0 
1182         lamp.data.energy               = 3.0 
1183         lamp.data.shadow_method        = 'RAY_SHADOW'
1184
1185         bpy.context.scene.world.light_settings.use_ambient_occlusion = True
1186         bpy.context.scene.world.light_settings.ao_factor = 0.2
1187
1188
1189
1190     #
1191     #
1192     #
1193     #
1194     #          SOME OUTPUT ON THE CONSOLE
1195     #
1196     #
1197     #
1198
1199    
1200     # The following two loops give a huge printout in the terminal. If needed one can uncomment these lines
1201    
1202     # Atoms
1203     # print("\nCoordinates of the atoms:")
1204     # for i in list(range(Number_of_total_atoms)):
1205     #   print(str(i+1) + "      " + str(atom_x[i]) + "  " + str(atom_y[i]) + "  " + str(atom_z[i]) + "  " + str(atom_R[i]) + "  " + atom_element[i])
1206
1207     # Sticks
1208     # print("\nSticks, which connect two atoms with indices:")
1209     # for i in list(range(Number_of_sticks)):
1210     #    print(str(bar_atom1[i]) + "   " + str(bar_atom2[i]))
1211    
1212     print()
1213     print()
1214     print()
1215     print(Atomic_Blender_string)
1216     print()
1217     print("Total number of atoms   : " + str(Number_of_total_atoms))
1218     print("Total number of sticks  : " + str(Number_of_sticks))
1219     print("Center of object        : ", object_center)
1220     print("Size of object          : ", object_size)
1221     print()
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236     #
1237     #
1238     #
1239     #
1240     #          DRAWING OF ATOMS
1241     #
1242     #
1243     #
1244
1245     # This part was the main issue in earlier versions: the loading of atoms has taken too much time!
1246     # Example: a surface can be easily composed of 5000 atoms => Loading 5000 NURBS needs quite a long time. 
1247     # There are two things that can be done: 
1248     #
1249     #                                1. Remove any 'if's in the drawing loop <= Increased a bit the speed
1250     #                                2. Group atoms of one type such that it is ONE object <= Huge speed increase
1251     #
1252     # I have done some measurements with a stopwatch: 'Grouping' is definitely much, much faster for
1253     # large scale situations. We talk about differences in seconds or even, for many-atom-scenes, in 
1254     # minutes and hours. Try out and compare yourself. To clearly see the difference you need, say, 2000 atoms.
1255     # Try with 8000 atoms => you need to group the atoms ... . 
1256
1257
1258     # Lists of atoms of one type are first created. If it is atoms, all theses lists are put into one 
1259     # single list called 'draw_atom_type_list'. The vacancies have their extra list 'draw_atom_type_list_vacancy' 
1260    
1261     # So far, the empty lists:
1262     # The list containing all lists, which each contains all atoms of one type
1263     draw_atom_type_list           = []
1264     # The list which contains all vacancies
1265     draw_atom_type_list_vacancy   = []
1266
1267
1268     # Go through the list which contains all types of atoms. It is the list, which has been
1269     # created on the top during reading the PDB file. It contains elements like [hydrogen, carbon, ...]
1270     for atom_type in atom_all_types_list:
1271    
1272         # This is the draw list, which contains all atoms of one type (e.g. all hydrogens) ...
1273         draw_atom_list = []  
1274       
1275         # Go through all atoms ...
1276         for i in range(0, Number_of_total_atoms):
1277             # ... select the atoms of the considered type via comparison ...
1278             if atom_type[0] == atom_name[i]:
1279                 # ... and append them to the list 'draw_atom_list'.
1280                 draw_atom_list.append([atom_name[i], atom_material[i], [atom_x[i], atom_y[i], atom_z[i]], atom_R[i]])
1281       
1282         if atom_type[0] == "Vacancy":
1283             draw_atom_type_list_vacancy.append(draw_atom_list)
1284         else:
1285             draw_atom_type_list.append(draw_atom_list)
1286    
1287
1288     #
1289     # DRAW ATOMS
1290     #
1291
1292     # The comments in the follwoing block of code (NURBS) are basically the same
1293     # for the code, which is used to draw meshes and the vacancies.
1294    
1295     # This is the number of all atoms which are put into the scene.
1296     i = 0 
1297     # Deselect all objects.
1298     bpy.ops.object.select_all(action='DESELECT')    
1299     # Draw NURBS or ...
1300     if mesh_yn == False:
1301         # For all atom lists in the list 'draw_atom_type_list' do ...
1302         # So, in other words: for all lists of atoms of ONE type (e.g. Hydrogen)
1303         for atom_list in draw_atom_type_list:
1304             group_counter = 0
1305             # This is a list of groups. Each group contains 'group_atoms_dn' 
1306             # atoms; let's say it is 200 atoms which shall be loaded first and then grouped
1307             # together.
1308             groups        = []
1309             # For each atom in a group do ...
1310             for atom in atom_list:
1311                 # First print in the terminal the number of atom that will be build right now.
1312                 sys.stdout.write("Atom No. %d has been built\r" % (i+1) )
1313                 sys.stdout.flush()
1314
1315                 # Build a sphere (atom), with coordinates 'atom[2]'. the index 2 is a list
1316                 # of the coordinates [x,y,z]
1317                 bpy.ops.surface.primitive_nurbs_surface_sphere_add(view_align=False, enter_editmode=False, location=atom[2], rotation=(0.0, 0.0, 0.0), layers=(True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False))
1318                 # Put the ball into the scene
1319                 ball                 = bpy.context.scene.objects.active
1320                 # Scale the radius (atom[3])
1321                 ball.scale           = (atom[3]*Ball_radius_factor,atom[3]*Ball_radius_factor,atom[3]*Ball_radius_factor)
1322                 # Name and ...
1323                 ball.name            = atom[0]
1324                 # Material
1325                 ball.active_material = atom[1]
1326                 # Very important for grouping: append the atom to the list 'atom_object'
1327                 atom_object.append(ball)
1328             
1329                 # We have a new atom, so increase the counter.
1330                 i += 1   
1331                 # Increase also the group counter      
1332                 group_counter += 1 
1333                 # If we shall group then do ... 
1334                 if FLAG_group_atoms == True:
1335                     # First check if we have loaded 200 atoms ('group_atoms_dn').
1336                     if group_counter == group_atoms_dn:
1337                         # If so, then activate all the last 200 atoms. i is the current number
1338                         # of all loaded atoms and group_counter is 200.
1339                         for z in range(i-group_counter,i):
1340                             atom_object[z].select = True 
1341                         # ... and join them. This here is the 'grouping'.
1342                         bpy.ops.object.join()
1343                         # Put this list of 200 atoms into the group list
1344                         groups.append(bpy.context.scene.objects.active)
1345                         # Deselect all objects and ...
1346                         bpy.ops.object.select_all(action='DESELECT')
1347                         # ... start from the beginning with the group counter. Come
1348                         # back if you have again 200 atoms loaded.
1349                         group_counter = 0
1350  
1351             # So, if we shall group then group the very last loaded atoms ...
1352             if FLAG_group_atoms == True:
1353                 for z in range(i-group_counter,i):
1354                     atom_object[z].select = True 
1355                     bpy.context.scene.objects.active = atom_object[z]
1356                 bpy.ops.object.join()
1357                 groups.append(bpy.context.scene.objects.active)
1358                 bpy.ops.object.select_all(action='DESELECT')
1359          
1360                 # Group all groups into one huge group! Select the before ...
1361                 for group in groups:
1362                     group.select = True 
1363                 bpy.ops.object.join()   
1364                 bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY')
1365             
1366     # ... draw Mesh balls  
1367     else: 
1368         for atom_list in draw_atom_type_list:
1369             group_counter = 0
1370             groups        = []
1371             for atom in atom_list:
1372                 sys.stdout.write("Atom No. %d has been built\r" % (i+1) )
1373                 sys.stdout.flush()
1374
1375                 bpy.ops.mesh.primitive_uv_sphere_add(segments=Ball_azimuth, ring_count=Ball_zenith, size=1, view_align=False, enter_editmode=False, location=atom[2], rotation=(0, 0, 0), layers=(True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False))
1376                 ball                 = bpy.context.scene.objects.active
1377                 ball.scale           = (atom[3]*Ball_radius_factor,atom[3]*Ball_radius_factor,atom[3]*Ball_radius_factor)
1378                 ball.name            = atom[0]
1379                 ball.active_material = atom[1]
1380                 atom_object.append(ball)
1381
1382                 i += 1         
1383                 group_counter += 1 
1384                 if FLAG_group_atoms == True:
1385                     if group_counter == group_atoms_dn:
1386                         for z in range(i-group_counter,i):
1387                             atom_object[z].select = True 
1388                         bpy.ops.object.join()
1389                         groups.append(bpy.context.scene.objects.active)
1390                         bpy.ops.object.select_all(action='DESELECT')
1391                         group_counter = 0
1392          
1393             if FLAG_group_atoms == True:
1394                 for z in range(i-group_counter,i):
1395                     atom_object[z].select = True 
1396                 bpy.ops.object.join()
1397                 groups.append(bpy.context.scene.objects.active)
1398                 bpy.ops.object.select_all(action='DESELECT')  
1399
1400                 for group in groups:
1401                     group.select = True 
1402                 bpy.ops.object.join()  
1403                 bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY')
1404
1405     #
1406     # DRAW VACANCIES
1407     #
1408
1409
1410     bpy.ops.object.select_all(action='DESELECT') 
1411     for atom_list in draw_atom_type_list_vacancy:
1412         group_counter = 0
1413         groups        = []
1414         for atom in atom_list:
1415             sys.stdout.write("Atom No. %d has been built\r" % (i+1) )
1416             sys.stdout.flush()
1417             bpy.ops.mesh.primitive_cube_add(view_align=False, enter_editmode=False, location=atom[2], rotation=(0.0, 0.0, 0.0), layers=(True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False))
1418             ball                 = bpy.context.scene.objects.active
1419             ball.scale           = (atom[3]*Ball_radius_factor,atom[3]*Ball_radius_factor,atom[3]*Ball_radius_factor)
1420             ball.name            = atom[0]
1421             ball.active_material = atom[1]
1422             atom_object.append(ball)
1423          
1424             i += 1         
1425             group_counter += 1 
1426             if FLAG_group_atoms == True:
1427                 if group_counter == group_atoms_dn:
1428                     for z in range(i-group_counter,i):
1429                         atom_object[z].select = True
1430                     bpy.ops.object.join()
1431                     groups.append(bpy.context.scene.objects.active)
1432                     bpy.ops.object.select_all(action='DESELECT')
1433                     group_counter = 0
1434                     
1435         if FLAG_group_atoms == True:
1436             for z in range(i-group_counter,i):
1437                 atom_object[z].select = True
1438             bpy.ops.object.join()
1439             groups.append(bpy.context.scene.objects.active)
1440             bpy.ops.object.select_all(action='DESELECT')  
1441          
1442             for group in groups:
1443                 group.select = True 
1444             bpy.ops.object.join()  
1445             bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY')
1446       
1447     print()    
1448       
1449       
1450       
1451       
1452       
1453       
1454       
1455       
1456     #
1457     #
1458     #
1459     #
1460     #          DRAWING OF STICKS
1461     #
1462     #
1463     #
1464
1465
1466
1467     if Stick_yn == True:
1468  
1469         # Create a new material with the corresponding color. The
1470         # color is taken from the all_atom list, it is the last entry
1471         # in the data file (index -1).
1472         bpy.ops.object.material_slot_add()
1473         stick_material               = bpy.data.materials.new(Data_Shortname[all_existing_atoms-1])
1474         stick_material.name          = Data_Atomname[all_existing_atoms-1]
1475         stick_material.diffuse_color = Data_Color [all_existing_atoms-1]
1476  
1477         # This is the unit vector of the z axis
1478         up_axis = mathutils.Vector([0.0, 0.0, 1.0])
1479  
1480         # For all sticks, do ...
1481         for i in range(0,Number_of_sticks):
1482             # Print on the terminal the actual number of the stick that is build
1483             sys.stdout.write("Stick No. %d has been built\r" % (i+1) )
1484             sys.stdout.flush()
1485             # The vectors of the two atoms are build 
1486             k1 = mathutils.Vector([atom_x[bar_atom1[i]-1],atom_y[bar_atom1[i]-1],atom_z[bar_atom1[i]-1]])
1487             k2 = mathutils.Vector([atom_x[bar_atom2[i]-1],atom_y[bar_atom2[i]-1],atom_z[bar_atom2[i]-1]])
1488             # This is the difference of both vectors
1489             v = (k2-k1)
1490             # Angle with respect to the z-axis
1491             angle   = v.angle(up_axis, 0)
1492             # Cross-product between v and the z-axis vector. It is the vector of
1493             # rotation.
1494             axis    = up_axis.cross(v)
1495             # Calculate Euler angles
1496             euler   = mathutils.Matrix.Rotation(angle, 4, axis).to_euler()
1497             # Create stick
1498             bpy.ops.mesh.primitive_cylinder_add(vertices=Stick_sectors, radius=Stick_diameter, depth= v.length, cap_ends=True, view_align=False, enter_editmode=False, location= ((k1+k2)*0.5), rotation=(0,0,0), layers=(True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False))
1499             # Put the stick into the scene ...
1500             stick                 = bpy.context.scene.objects.active
1501             # ... and rotate the stick.
1502             stick.rotation_euler  = euler
1503             # Material ... 
1504             stick.active_material = stick_material
1505             # ... and name
1506             stick.name            = Data_Atomname[all_existing_atoms-1]
1507
1508     print()
1509     print()
1510     print("All atoms and sticks have been drawn - finished.")
1511     print()
1512     print()
1513
1514
1515     return str(Number_of_total_atoms)
1516             
1517