1
2
3
4
5
6
7
8
9
10
11 from rdkit import Chem
12 from rdkit import RDConfig
13 import numpy
14 import math
15 import sys
16 import copy
17 import pprint
18 from rdkit.six import cmp
19
20 periodicTable=Chem.GetPeriodicTable()
21
32
34 dotsPerAngstrom= 30
35 useFraction= 0.85
36
37 atomLabelFontFace= "sans"
38 atomLabelFontSize= 12
39 atomLabelMinFontSize= 7
40
41 bondLineWidth= 1.2
42 dblBondOffset= .25
43 dblBondLengthFrac= .8
44
45 defaultColor= (1,0,0)
46 selectColor= (1,0,0)
47 bgColor = (1,1,1)
48
49 colorBonds= True
50 noCarbonSymbols= True
51 includeAtomNumbers= False
52 atomNumberOffset= 0
53 radicalSymbol= u'\u2219'
54
55 dash= (4,4)
56
57 wedgeDashedBonds= True
58 showUnknownDoubleBonds= True
59
60
61
62 coordScale= 1.0
63
64 elemDict={
65 1:(0.55,0.55,0.55),
66 7:(0,0,1),
67 8:(1,0,0),
68 9:(.2,.8,.8),
69 15:(1,.5,0),
70 16:(.8,.8,0),
71 17:(0,.8,0),
72 35:(.5,.3,.1),
73 53:(.63,.12,.94),
74 0:(.5,.5,.5),
75 }
76
78 atomPs = None
79 canvas = None
80 canvasSize = None
81
82 - def __init__(self,canvas=None, drawingOptions=None):
83 self.canvas = canvas
84 if canvas:
85 self.canvasSize=canvas.size
86 self.atomPs = {}
87 if drawingOptions is None:
88 self.drawingOptions=DrawingOptions()
89 else:
90 self.drawingOptions=drawingOptions
91 self.boundingBoxes = {}
92
93 if self.drawingOptions.bgColor is not None:
94 self.canvas.addCanvasPolygon(((0,0),
95 (canvas.size[0],0),
96 (canvas.size[0],canvas.size[1]),
97 (0,canvas.size[1])),
98 color=self.drawingOptions.bgColor,
99 fill=True,stroke=False)
100
101
108
109
111
112 dx = p2[0]-p1[0]
113 dy = p2[1]-p1[1]
114
115
116 ang = math.atan2(dy,dx)
117 perp = ang + math.pi/2.
118
119
120 offsetX = math.cos(perp)*self.drawingOptions.dblBondOffset*self.currDotsPerAngstrom
121 offsetY = math.sin(perp)*self.drawingOptions.dblBondOffset*self.currDotsPerAngstrom
122
123 return perp,offsetX,offsetY
124
128 if not lenFrac:
129 lenFrac = self.drawingOptions.dblBondLengthFrac
130
131 dx = p2[0]-p1[0]
132 dy = p2[1]-p1[1]
133
134
135
136
137 fracP1 = p1[0]+offsetX,p1[1]+offsetY
138
139
140 frac = (1.-lenFrac)/2
141 fracP1 = fracP1[0]+dx*frac,\
142 fracP1[1]+dy*frac
143
144 fracP2 = fracP1[0]+dx*lenFrac,\
145 fracP1[1]+dy*lenFrac
146 return fracP1,fracP2
147
148 - def _offsetDblBond(self,p1,p2,bond,a1,a2,conf,dir=1,
149 lenFrac=None):
150 perp,offsetX,offsetY = self._getBondOffset(p1,p2)
151 offsetX = offsetX*dir
152 offsetY = offsetY*dir
153
154
155 if bond.IsInRing():
156 bondIdx = bond.GetIdx()
157 a1Idx = a1.GetIdx()
158 a2Idx = a2.GetIdx()
159
160 for otherBond in a1.GetBonds():
161 if otherBond.GetIdx()!=bondIdx and \
162 otherBond.IsInRing():
163 sharedRing=False
164 for ring in self.bondRings:
165 if bondIdx in ring and otherBond.GetIdx() in ring:
166 sharedRing=True
167 break
168 if not sharedRing:
169 continue
170 a3 = otherBond.GetOtherAtom(a1)
171 if a3.GetIdx() != a2Idx:
172 p3 = self.transformPoint(conf.GetAtomPosition(a3.GetIdx())*self.drawingOptions.coordScale)
173 dx2 = p3[0] - p1[0]
174 dy2 = p3[1] - p1[1]
175 dotP = dx2*offsetX + dy2*offsetY
176 if dotP < 0:
177 perp += math.pi
178 offsetX = math.cos(perp)*self.drawingOptions.dblBondOffset*self.currDotsPerAngstrom
179 offsetY = math.sin(perp)*self.drawingOptions.dblBondOffset*self.currDotsPerAngstrom
180
181 fracP1,fracP2 = self._getOffsetBondPts(p1,p2,
182 offsetX,offsetY,
183 lenFrac=lenFrac)
184 return fracP1,fracP2
185
187 newpos = [None, None]
188 if labelSize != None:
189 labelSizeOffset = [labelSize[0][0]/2 + (cmp(p2[0], p1[0]) * labelSize[0][2]), labelSize[0][1]/2]
190 if p1[1] == p2[1]:
191 newpos[0] = p1[0] + cmp(p2[0], p1[0]) * labelSizeOffset[0]
192 else:
193 if abs(labelSizeOffset[1] * (p2[0] - p1[0]) / (p2[1] - p1[1])) < labelSizeOffset[0]:
194 newpos[0] = p1[0] + cmp(p2[0], p1[0]) * abs(labelSizeOffset[1] * (p2[0] - p1[0]) / (p2[1] - p1[1]))
195 else:
196 newpos[0] = p1[0] + cmp(p2[0], p1[0]) * labelSizeOffset[0]
197 if p1[0] == p2[0]:
198 newpos[1] = p1[1] + cmp(p2[1], p1[1]) * labelSizeOffset[1]
199 else:
200 if abs(labelSizeOffset[0] * (p1[1] - p2[1]) / (p2[0] - p1[0])) < labelSizeOffset[1]:
201 newpos[1] = p1[1] + cmp(p2[1], p1[1]) * abs(labelSizeOffset[0] * (p1[1] - p2[1]) / (p2[0] - p1[0]))
202 else:
203 newpos[1] = p1[1] + cmp(p2[1], p1[1]) * labelSizeOffset[1]
204 else:
205 newpos = copy.deepcopy(p1)
206 return newpos
207
208 - def _drawWedgedBond(self,bond,pos,nbrPos,
209 width=None,color=None,
210 dash=None):
211 if width is None:
212 width = self.drawingOptions.bondLineWidth
213 if color is None:
214 color = self.drawingOptions.defaultColor
215 perp,offsetX,offsetY = self._getBondOffset(pos,nbrPos)
216 offsetX *=.75
217 offsetY *=.75
218 poly = ((pos[0],pos[1]),
219 (nbrPos[0]+offsetX,nbrPos[1]+offsetY),
220 (nbrPos[0]-offsetX,nbrPos[1]-offsetY))
221
222 if not dash:
223 self.canvas.addCanvasPolygon(poly,color=color)
224 elif self.drawingOptions.wedgeDashedBonds and self.canvas.addCanvasDashedWedge:
225 self.canvas.addCanvasDashedWedge(poly[0],poly[1],poly[2],color=color)
226 else:
227 self.canvas.addCanvasLine(pos,nbrPos,linewidth=width*2,color=color,
228 dashes=dash)
229
230 - def _drawBond(self,bond,atom,nbr,pos,nbrPos,conf,
231 width=None,color=None,color2=None,labelSize1=None,labelSize2=None):
232 if width is None:
233 width = self.drawingOptions.bondLineWidth
234 if color is None:
235 color = self.drawingOptions.defaultColor
236 p1_raw = copy.deepcopy(pos)
237 p2_raw = copy.deepcopy(nbrPos)
238 newpos = self._getBondAttachmentCoordinates(p1_raw, p2_raw, labelSize1)
239 newnbrPos = self._getBondAttachmentCoordinates(p2_raw, p1_raw, labelSize2)
240 bType=bond.GetBondType()
241 if bType == Chem.BondType.SINGLE:
242 bDir = bond.GetBondDir()
243 if bDir in (Chem.BondDir.BEGINWEDGE,Chem.BondDir.BEGINDASH):
244
245 if bond.GetBeginAtom().GetChiralTag() in (Chem.ChiralType.CHI_TETRAHEDRAL_CW,
246 Chem.ChiralType.CHI_TETRAHEDRAL_CCW):
247 p1,p2 = newpos,newnbrPos
248 wcolor=color
249 else:
250 p2,p1 = newpos,newnbrPos
251 if color2 is not None:
252 wcolor=color2
253 else:
254 wcolor=self.drawingOptions.defaultColor
255 if bDir==Chem.BondDir.BEGINWEDGE:
256 self._drawWedgedBond(bond,p1,p2,color=wcolor,width=width)
257 elif bDir==Chem.BondDir.BEGINDASH:
258 self._drawWedgedBond(bond,p1,p2,color=wcolor,width=width,
259 dash=self.drawingOptions.dash)
260
261 else:
262 self.canvas.addCanvasLine(newpos, newnbrPos, linewidth=width, color=color, color2=color2)
263 elif bType == Chem.BondType.DOUBLE:
264 crossBond = (self.drawingOptions.showUnknownDoubleBonds and \
265 bond.GetStereo() == Chem.BondStereo.STEREOANY)
266 if not crossBond and \
267 ( bond.IsInRing() or (atom.GetDegree()!=1 and bond.GetOtherAtom(atom).GetDegree()!=1) ):
268 self.canvas.addCanvasLine(newpos,newnbrPos,linewidth=width,color=color,color2=color2)
269 fp1,fp2 = self._offsetDblBond(newpos,newnbrPos,bond,atom,nbr,conf)
270 self.canvas.addCanvasLine(fp1,fp2,linewidth=width,color=color,color2=color2)
271 else:
272 fp1,fp2 = self._offsetDblBond(newpos,newnbrPos,bond,atom,nbr,conf,dir=.5,
273 lenFrac=1.0)
274 fp3,fp4 = self._offsetDblBond(newpos,newnbrPos,bond,atom,nbr,conf,dir=-.5,
275 lenFrac=1.0)
276 if crossBond:
277 fp2,fp4=fp4,fp2
278 self.canvas.addCanvasLine(fp1,fp2,linewidth=width,color=color,color2=color2)
279 self.canvas.addCanvasLine(fp3,fp4,linewidth=width,color=color,color2=color2)
280
281 elif bType == Chem.BondType.AROMATIC:
282 self.canvas.addCanvasLine(newpos,newnbrPos,linewidth=width,color=color,color2=color2)
283 fp1,fp2 = self._offsetDblBond(newpos,newnbrPos,bond,atom,nbr,conf)
284 self.canvas.addCanvasLine(fp1,fp2,linewidth=width,color=color,color2=color2,
285 dash=self.drawingOptions.dash)
286 elif bType == Chem.BondType.TRIPLE:
287 self.canvas.addCanvasLine(newpos,newnbrPos,linewidth=width,color=color,color2=color2)
288 fp1,fp2 = self._offsetDblBond(newpos,newnbrPos,bond,atom,nbr,conf)
289 self.canvas.addCanvasLine(fp1,fp2,linewidth=width,color=color,color2=color2)
290 fp1,fp2 = self._offsetDblBond(newpos,newnbrPos,bond,atom,nbr,conf,dir=-1)
291 self.canvas.addCanvasLine(fp1,fp2,linewidth=width,color=color,color2=color2)
292 else:
293 self.canvas.addCanvasLine(newpos, newnbrPos, linewidth=width, color=color, color2=color2,
294 dash=(1,2))
295
296 - def scaleAndCenter(self,mol,conf,coordCenter=False,canvasSize=None,ignoreHs=False):
297 if canvasSize is None:
298 canvasSize=self.canvasSize
299 xAccum = 0
300 yAccum = 0
301 minX = 1e8
302 minY = 1e8
303 maxX = -1e8
304 maxY = -1e8
305
306 nAts = mol.GetNumAtoms()
307 for i in range(nAts):
308 if ignoreHs and mol.GetAtomWithIdx(i).GetAtomicNum()==1: continue
309 pos = conf.GetAtomPosition(i)*self.drawingOptions.coordScale
310 xAccum += pos[0]
311 yAccum += pos[1]
312 minX = min(minX,pos[0])
313 minY = min(minY,pos[1])
314 maxX = max(maxX,pos[0])
315 maxY = max(maxY,pos[1])
316
317 dx = abs(maxX-minX)
318 dy = abs(maxY-minY)
319 xSize = dx*self.currDotsPerAngstrom
320 ySize = dy*self.currDotsPerAngstrom
321
322 if coordCenter:
323 molTrans = -xAccum/nAts,-yAccum/nAts
324 else:
325 molTrans = -(minX+(maxX-minX)/2),-(minY+(maxY-minY)/2)
326 self.molTrans = molTrans
327
328 if xSize>=.95*canvasSize[0]:
329 scale = .9*canvasSize[0]/xSize
330 xSize*=scale
331 ySize*=scale
332 self.currDotsPerAngstrom*=scale
333 self.currAtomLabelFontSize = max(self.currAtomLabelFontSize*scale,
334 self.drawingOptions.atomLabelMinFontSize)
335 if ySize>=.95*canvasSize[1]:
336 scale = .9*canvasSize[1]/ySize
337 xSize*=scale
338 ySize*=scale
339 self.currDotsPerAngstrom*=scale
340 self.currAtomLabelFontSize = max(self.currAtomLabelFontSize*scale,
341 self.drawingOptions.atomLabelMinFontSize)
342 drawingTrans = canvasSize[0]/2,canvasSize[1]/2
343 self.drawingTrans = drawingTrans
344
345 - def _drawLabel(self,label,pos,baseOffset,font,color=None,**kwargs):
346 if color is None:
347 color = self.drawingOptions.defaultColor
348 x1 = pos[0]
349 y1 = pos[1]
350 labelP = x1,y1
351 labelSize = self.canvas.addCanvasText(label,(x1,y1,baseOffset),font,color,**kwargs)
352 return labelSize
353
354 - def AddMol(self,mol,centerIt=True,molTrans=None,drawingTrans=None,
355 highlightAtoms=[],confId=-1,flagCloseContactsDist=2,
356 highlightMap=None, ignoreHs=False,highlightBonds=[],**kwargs):
357 """Set the molecule to be drawn.
358
359 Parameters:
360 hightlightAtoms -- list of atoms to highlight (default [])
361 highlightMap -- dictionary of (atom, color) pairs (default None)
362
363 Notes:
364 - specifying centerIt will cause molTrans and drawingTrans to be ignored
365 """
366 conf = mol.GetConformer(confId)
367 if 'coordScale' in kwargs:
368 self.drawingOptions.coordScale=kwargs['coordScale']
369
370 self.currDotsPerAngstrom=self.drawingOptions.dotsPerAngstrom
371 self.currAtomLabelFontSize=self.drawingOptions.atomLabelFontSize
372 if centerIt:
373 self.scaleAndCenter(mol,conf,ignoreHs=ignoreHs)
374 else:
375 if molTrans is None:
376 molTrans = (0,0)
377 self.molTrans = molTrans
378 if drawingTrans is None:
379 drawingTrans = (0,0)
380 self.drawingTrans = drawingTrans
381
382 font = Font(face=self.drawingOptions.atomLabelFontFace,size=self.currAtomLabelFontSize)
383
384 obds=None
385 if not mol.HasProp('_drawingBondsWedged'):
386
387 obds=[x.GetBondDir() for x in mol.GetBonds()]
388 Chem.WedgeMolBonds(mol,conf)
389
390 includeAtomNumbers = kwargs.get('includeAtomNumbers',self.drawingOptions.includeAtomNumbers)
391 self.atomPs[mol] = {}
392 self.boundingBoxes[mol] = [0]*4
393 self.activeMol = mol
394 self.bondRings = mol.GetRingInfo().BondRings()
395 labelSizes = {}
396 for atom in mol.GetAtoms():
397 labelSizes[atom.GetIdx()] = None
398 if ignoreHs and atom.GetAtomicNum()==1:
399 drawAtom=False
400 else:
401 drawAtom=True
402 idx = atom.GetIdx()
403 pos = self.atomPs[mol].get(idx,None)
404 if pos is None:
405 pos = self.transformPoint(conf.GetAtomPosition(idx)*self.drawingOptions.coordScale)
406 self.atomPs[mol][idx] = pos
407 if drawAtom:
408 self.boundingBoxes[mol][0]=min(self.boundingBoxes[mol][0],pos[0])
409 self.boundingBoxes[mol][1]=min(self.boundingBoxes[mol][1],pos[1])
410 self.boundingBoxes[mol][2]=max(self.boundingBoxes[mol][2],pos[0])
411 self.boundingBoxes[mol][3]=max(self.boundingBoxes[mol][3],pos[1])
412
413 if not drawAtom: continue
414 nbrSum = [0,0]
415 for bond in atom.GetBonds():
416 nbr = bond.GetOtherAtom(atom)
417 if ignoreHs and nbr.GetAtomicNum()==1: continue
418 nbrIdx = nbr.GetIdx()
419 if nbrIdx > idx:
420 nbrPos = self.atomPs[mol].get(nbrIdx,None)
421 if nbrPos is None:
422 nbrPos = self.transformPoint(conf.GetAtomPosition(nbrIdx)*self.drawingOptions.coordScale)
423 self.atomPs[mol][nbrIdx] = nbrPos
424 self.boundingBoxes[mol][0]=min(self.boundingBoxes[mol][0],nbrPos[0])
425 self.boundingBoxes[mol][1]=min(self.boundingBoxes[mol][1],nbrPos[1])
426 self.boundingBoxes[mol][2]=max(self.boundingBoxes[mol][2],nbrPos[0])
427 self.boundingBoxes[mol][3]=max(self.boundingBoxes[mol][3],nbrPos[1])
428
429 else:
430 nbrPos = self.atomPs[mol][nbrIdx]
431 nbrSum[0] += nbrPos[0]-pos[0]
432 nbrSum[1] += nbrPos[1]-pos[1]
433
434 iso = atom.GetIsotope()
435
436 labelIt= not self.drawingOptions.noCarbonSymbols or \
437 atom.GetAtomicNum()!=6 or \
438 atom.GetFormalCharge()!=0 or \
439 atom.GetNumRadicalElectrons() or \
440 includeAtomNumbers or \
441 iso or \
442 atom.HasProp('molAtomMapNumber') or \
443 atom.GetDegree()==0
444 orient=''
445 if labelIt:
446 baseOffset = 0
447 if includeAtomNumbers:
448 symbol = str(atom.GetIdx())
449 symbolLength = len(symbol)
450 else:
451 base = atom.GetSymbol()
452 symbolLength = len(base)
453 if not atom.HasQuery():
454 nHs = atom.GetTotalNumHs()
455 else:
456 nHs = 0
457 if nHs>0:
458 if nHs>1:
459 hs='H<sub>%d</sub>'%nHs
460 symbolLength += 1 + len(str(nHs))
461 else:
462 hs ='H'
463 symbolLength += 1
464 else:
465 hs = ''
466 chg = atom.GetFormalCharge()
467 if chg!=0:
468 if chg==1:
469 chg = '+'
470 elif chg==-1:
471 chg = '-'
472 elif chg>1:
473 chg = '+%d'%chg
474 elif chg<-1:
475 chg = '-%d'%chg
476 symbolLength += len(chg)
477 else:
478 chg = ''
479 if chg:
480 chg = '<sup>%s</sup>'%chg
481
482 if atom.GetNumRadicalElectrons():
483 rad = self.drawingOptions.radicalSymbol*atom.GetNumRadicalElectrons()
484 rad = '<sup>%s</sup>'%rad
485 symbolLength += atom.GetNumRadicalElectrons()
486 else:
487 rad = ''
488
489 isotope=''
490 isotopeLength = 0
491 if iso:
492 isotope='<sup>%d</sup>'%atom.GetIsotope()
493 isotopeLength = len(str(atom.GetIsotope()))
494 symbolLength += isotopeLength
495 mapNum=''
496 mapNumLength = 0
497 if atom.HasProp('molAtomMapNumber'):
498 mapNum=':'+atom.GetProp('molAtomMapNumber')
499 mapNumLength = 1 + len(str(atom.GetProp('molAtomMapNumber')))
500 symbolLength += mapNumLength
501 deg = atom.GetDegree()
502
503
504
505 if deg==0:
506 if periodicTable.GetElementSymbol(atom.GetAtomicNum()) in ('O','S','Se','Te','F','Cl','Br','I','At'):
507 symbol = '%s%s%s%s%s%s'%(hs,isotope,base,chg,rad,mapNum)
508 else:
509 symbol = '%s%s%s%s%s%s'%(isotope,base,hs,chg,rad,mapNum)
510 elif deg>1 or nbrSum[0]<1:
511 symbol = '%s%s%s%s%s%s'%(isotope,base,hs,chg,rad,mapNum)
512 baseOffset = 0.5 - (isotopeLength + len(base) / 2.) / symbolLength
513 else:
514 symbol = '%s%s%s%s%s%s'%(rad,chg,hs,isotope,base,mapNum)
515 baseOffset = -0.5 + (mapNumLength + len(base) / 2.) / symbolLength
516 if deg==1:
517 if abs(nbrSum[1])>1:
518 islope=nbrSum[0]/abs(nbrSum[1])
519 else:
520 islope=nbrSum[0]
521 if abs(islope)>.3:
522 if islope>0:
523 orient='W'
524 else:
525 orient='E'
526 elif abs(nbrSum[1])>10:
527 if nbrSum[1]>0:
528 orient='N'
529 else :
530 orient='S'
531 else:
532 orient = 'C'
533 if highlightMap and idx in highlightMap:
534 color = highlightMap[idx]
535 elif highlightAtoms and idx in highlightAtoms:
536 color = self.drawingOptions.selectColor
537 else:
538 color = self.drawingOptions.elemDict.get(atom.GetAtomicNum(),(0,0,0))
539 labelSize = self._drawLabel(symbol, pos, baseOffset, font, color=color,orientation=orient)
540 labelSizes[atom.GetIdx()] = [labelSize, orient]
541
542 for bond in mol.GetBonds():
543 atom, idx = bond.GetBeginAtom(), bond.GetBeginAtomIdx()
544 nbr, nbrIdx = bond.GetEndAtom(), bond.GetEndAtomIdx()
545 pos = self.atomPs[mol].get(idx,None)
546 nbrPos = self.atomPs[mol].get(nbrIdx,None)
547 if highlightBonds and bond.GetIdx() in highlightBonds:
548 width=2.0*self.drawingOptions.bondLineWidth
549 color = self.drawingOptions.selectColor
550 color2 = self.drawingOptions.selectColor
551 elif highlightAtoms and idx in highlightAtoms and nbrIdx in highlightAtoms:
552 width=2.0*self.drawingOptions.bondLineWidth
553 color = self.drawingOptions.selectColor
554 color2 = self.drawingOptions.selectColor
555 elif highlightMap is not None and idx in highlightMap and nbrIdx in highlightMap:
556 width=2.0*self.drawingOptions.bondLineWidth
557 color = highlightMap[idx]
558 color2 = highlightMap[nbrIdx]
559 else:
560 width=self.drawingOptions.bondLineWidth
561 if self.drawingOptions.colorBonds:
562 color = self.drawingOptions.elemDict.get(atom.GetAtomicNum(),(0,0,0))
563 color2 = self.drawingOptions.elemDict.get(nbr.GetAtomicNum(),(0,0,0))
564 else:
565 color = self.drawingOptions.defaultColor
566 color2= color
567 self._drawBond(bond,atom,nbr,pos,nbrPos,conf,
568 color=color,width=width,color2=color2,labelSize1=labelSizes[idx],labelSize2=labelSizes[nbrIdx])
569
570
571 if obds:
572 for i,d in enumerate(obds):
573 mol.GetBondWithIdx(i).SetBondDir(d)
574
575
576 if flagCloseContactsDist>0:
577 tol = flagCloseContactsDist*flagCloseContactsDist
578 for i,atomi in enumerate(mol.GetAtoms()):
579 pi = numpy.array(self.atomPs[mol][i])
580 for j in range(i+1,mol.GetNumAtoms()):
581 pj = numpy.array(self.atomPs[mol][j])
582 d = pj-pi
583 dist2 = d[0]*d[0]+d[1]*d[1]
584 if dist2<=tol:
585 self.canvas.addCanvasPolygon(((pi[0]-2*flagCloseContactsDist,
586 pi[1]-2*flagCloseContactsDist),
587 (pi[0]+2*flagCloseContactsDist,
588 pi[1]-2*flagCloseContactsDist),
589 (pi[0]+2*flagCloseContactsDist,
590 pi[1]+2*flagCloseContactsDist),
591 (pi[0]-2*flagCloseContactsDist,
592 pi[1]+2*flagCloseContactsDist)),
593 color=(1.,0,0),
594 fill=False,stroke=True)
595