« back to all changes in this revision
Viewing changes to src/cluster.cpp
- browse files at revision 13
- compare with another revision
- download diff
- download tarball
- view history from revision 13
- Committer: Package Import Robot
- Author(s): Daniel Leidert, Daniel Leidert, Michael Banck
- Date: 2014-09-22 10:00:17 UTC
- mfrom: (1.1.9)
- Revision ID: package-import@ubuntu.com-20140922100017-voentkpft033vxa6
Tags: 140902-1
* New upstream release.
[ Daniel Leidert ]
* debian/copyright: Updated.
* debian/upstream: Renamed to debian/upstream/metadata.
[ Michael Banck ]
* debian/travis.1: Add LAMMPS and DLPOLY to trajectory filetypes and
document -stream option.
* debian/control (Description): List supported fileteypes and
available analysis features.
[ Daniel Leidert ]
* debian/copyright: Updated.
* debian/upstream: Renamed to debian/upstream/metadata.
[ Michael Banck ]
* debian/travis.1: Add LAMMPS and DLPOLY to trajectory filetypes and
document -stream option.
* debian/control (Description): List supported fileteypes and
available analysis features.
- files added:
- debian/upstream/metadata
- files removed:
- src/cluster.cpp
- files modified:
- Makefile
- debian/upstream
added
removed
src/cluster.cpp
1
/*****************************************************************************
2
TRAVIS - Trajectory Analyzer and Visualizer
3
http://www.travis-analyzer.de/
4
5
Copyright (c) 2009-2014 Martin Brehm
6
2012-2014 Martin Thomas
7
8
This file written by Martin Brehm.
9
10
This program is free software: you can redistribute it and/or modify
11
it under the terms of the GNU General Public License as published by
12
the Free Software Foundation, either version 3 of the License, or
13
(at your option) any later version.
14
15
This program is distributed in the hope that it will be useful,
16
but WITHOUT ANY WARRANTY; without even the implied warranty of
17
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18
GNU General Public License for more details.
19
20
You should have received a copy of the GNU General Public License
21
along with this program. If not, see <http://www.gnu.org/licenses/>.
22
*****************************************************************************/
23
24
25
#include "cluster.h"
26
#include "grace.h"
27
#include "maintools.h"
28
#include "globalvar.h"
29
30
31
CClusterNode::CClusterNode()
32
{
33
m_vaAtoms.SetName("CClusterNode::m_vaAtoms");
34
m_iaAtoms.SetName("CClusterNode::m_iaAtoms");
35
}
36
37
38
CClusterNode::~CClusterNode()
39
{
40
}
41
42
43
CClusterAnalysis::CClusterAnalysis()
44
{
45
m_poaClusterPoly = NULL;
46
m_pDistCache = NULL;
47
m_pClusterDistanceDF = NULL;
48
m_pPolymerDF = NULL;
49
m_pClusterCountDF = NULL;
50
m_pClusterDistributionDF = NULL;
51
m_pClusterSizeDF = NULL;
52
m_pClusterDistance2D = NULL;
53
m_pClusterDistribution2D = NULL;
54
m_pClusterSize2D = NULL;
55
m_pClusterCount2D = NULL;
56
m_pRandom = NULL;
57
58
m_oaAtomGroups.SetName("CClusterAnalysis::m_oaAtomGroups");
59
m_iaAtomList.SetName("CClusterAnalysis::m_iaAtomList");
60
m_oaClusters.SetName("CClusterAnalysis::m_oaClusters");
61
m_oaBaseClusters.SetName("CClusterAnalysis::m_oaBaseClusters");
62
m_oaTopClusters.SetName("CClusterAnalysis::m_oaTopClusters");
63
m_faMinDist.SetName("CClusterAnalysis::m_faMinDist");
64
m_faMaxDist.SetName("CClusterAnalysis::m_faMaxDist");
65
m_faHetNumber.SetName("CClusterAnalysis::m_faHetNumber");
66
m_faCenter.SetName("CClusterAnalysis::m_faCenter");
67
m_faIntHetNumber.SetName("CClusterAnalysis::m_faIntHetNumber");
68
}
69
70
71
CClusterAnalysis::~CClusterAnalysis()
72
{
73
if (m_poaClusterPoly != NULL)
74
delete[] m_poaClusterPoly;
75
if (m_pDistCache != NULL)
76
delete[] m_pDistCache;
77
}
78
79
80
void CClusterAnalysis::AddCluster(int i)
81
{
82
CClusterNode *n;
83
84
try { n = new CClusterNode(); } catch(...) { n = NULL; }
85
if (n == NULL) NewException((double)sizeof(CClusterNode),__FILE__,__LINE__,__PRETTY_FUNCTION__);
86
87
// n->m_iParent = -1;
88
n->m_iIndex = m_oaClusters.GetSize();
89
n->m_fDistance = 1.0E30f;
90
n->m_iMonomers = 1;
91
n->m_pParent = NULL;
92
if (m_bDistCache)
93
n->m_iaAtoms.SetMaxSize(i);
94
else n->m_vaAtoms.SetMaxSize(i);
95
n->m_iChildren[0] = -1;
96
97
m_oaClusters.Add(n);
98
m_oaBaseClusters.Add(n);
99
m_poaClusterPoly[0].Add(n);
100
}
101
102
103
void CClusterAnalysis::AddParticle(float x, float y, float z)
104
{
105
CClusterNode *n;
106
107
n = (CClusterNode*)m_oaClusters[m_oaClusters.GetSize()-1];
108
n->m_vaAtoms.Add(CxVector3(x,y,z));
109
}
110
111
112
void CClusterAnalysis::AddParticle(int i)
113
{
114
CClusterNode *n;
115
116
n = (CClusterNode*)m_oaClusters[m_oaClusters.GetSize()-1];
117
n->m_iaAtoms.Add(i);
118
}
119
120
121
void CClusterAnalysis::BuildTree()
122
{
123
int z, i, j, st/*, m*/;
124
float d;
125
CClusterNode *n, *c1, *c2;
126
127
n = NULL;
128
129
// mprintf("BuildTree online.\n");
130
for (z=0;z<m_oaClusters.GetSize();z++)
131
m_oaTopClusters.Add(m_oaClusters[z]);
132
133
for (z=0;z<m_oaClusters.GetSize();z++)
134
{
135
n = (CClusterNode*)m_oaClusters[z];
136
FindNearestNeighbor(n);
137
}
138
139
// mprintf("Initialization done.\n");
140
141
st = 0;
142
while (m_oaTopClusters.GetSize() > 1)
143
{
144
st++;
145
// mprintf("Step %d, %d Top-Clusters remaining.\n",st,m_oaTopClusters.GetSize());
146
d = 1.0E30f;
147
i = -1;
148
for (z=0;z<m_oaTopClusters.GetSize();z++)
149
{
150
if (((CClusterNode*)m_oaTopClusters[z])->m_fNNDistance < d)
151
{
152
d = ((CClusterNode*)m_oaTopClusters[z])->m_fNNDistance;
153
i = z;
154
}
155
}
156
if (i == -1)
157
abort();
158
j = ((CClusterNode*)m_oaTopClusters[i])->m_iNextNeighbor;
159
160
c1 = (CClusterNode*)m_oaTopClusters[i];
161
c2 = (CClusterNode*)m_oaClusters[j];
162
163
try { n = new CClusterNode(); } catch(...) { n = NULL; }
164
if (n == NULL) NewException((double)sizeof(CClusterNode),__FILE__,__LINE__,__PRETTY_FUNCTION__);
165
166
if (m_bDistCache)
167
n->m_iaAtoms.SetMaxSize(c1->m_iaAtoms.GetSize()+c2->m_iaAtoms.GetSize());
168
else n->m_vaAtoms.SetMaxSize(c1->m_vaAtoms.GetSize()+c2->m_vaAtoms.GetSize());
169
170
n->m_pParent = NULL;
171
// n->m_iParent = -1;
172
n->m_fDistance = d;
173
// n->m_laChildren.Add(((CClusterNode*)m_oaTopClusters[i])->m_iIndex);
174
// n->m_laChildren.Add(j);
175
n->m_iChildren[0] = c1->m_iIndex;
176
n->m_iChildren[1] = j;
177
n->m_iMonomers = c1->m_iMonomers + c2->m_iMonomers;
178
// m = n->m_iMonomers - 1;
179
m_poaClusterPoly[n->m_iMonomers-1].Add(n);
180
// mprintf(" Combining %d and %d, d=%f...\n",((CClusterNode*)m_oaTopClusters[i])->m_iIndex+1,j+1,d);
181
n->m_iIndex = m_oaClusters.GetSize();
182
183
/* m_pClusterSizeDF->AddToBinInt_Fast(m);
184
185
m_pClusterDistribution2DF->AddToBinInt_Fast(m,d-((CClusterNode*)m_oaTopClusters[i])->m_fDistance);
186
187
if (((CClusterNode*)m_oaTopClusters[i])->m_iMonomers == 1)
188
m_pClusterDistribution2DF->AddToBinInt_Fast(0,d);
189
else m_pClusterDistribution2DF->AddToBinInt_Fast(((CClusterNode*)m_oaTopClusters[i])->m_iMonomers-1,d-((CClusterNode*)m_oaTopClusters[i])->m_fDistance);
190
191
if (((CClusterNode*)m_oaClusters[j])->m_iMonomers == 1)
192
m_pClusterDistribution2DF->AddToBinInt_Fast(0,d);
193
else m_pClusterDistribution2DF->AddToBinInt_Fast(((CClusterNode*)m_oaClusters[j])->m_iMonomers-1,d-((CClusterNode*)m_oaClusters[j])->m_fDistance);
194
*/
195
// ((CClusterNode*)m_oaTopClusters[i])->m_iParent = n->m_iIndex;
196
((CClusterNode*)m_oaTopClusters[i])->m_pParent = n;
197
198
// ((CClusterNode*)m_oaClusters[j])->m_iParent = n->m_iIndex;
199
((CClusterNode*)m_oaClusters[j])->m_pParent = n;
200
201
if (m_bDistCache)
202
{
203
for (z=0;z<c1->m_iaAtoms.GetSize();z++)
204
n->m_iaAtoms.Add(c1->m_iaAtoms[z]);
205
206
for (z=0;z<c2->m_iaAtoms.GetSize();z++)
207
n->m_iaAtoms.Add(c2->m_iaAtoms[z]);
208
} else
209
{
210
for (z=0;z<c1->m_vaAtoms.GetSize();z++)
211
n->m_vaAtoms.Add(c1->m_vaAtoms[z]);
212
213
for (z=0;z<c2->m_vaAtoms.GetSize();z++)
214
n->m_vaAtoms.Add(c2->m_vaAtoms[z]);
215
}
216
217
m_oaClusters.Add(n);
218
m_oaTopClusters.Add(n);
219
m_oaTopClusters.RemoveAt_NoShrink(i,1);
220
for (z=0;z<m_oaTopClusters.GetSize();z++)
221
{
222
if (((CClusterNode*)m_oaTopClusters[z])->m_iIndex == j)
223
{
224
m_oaTopClusters.RemoveAt_NoShrink(z,1);
225
break;
226
}
227
}
228
// for (z=0;z<m_oaTopClusters.GetSize();z++)
229
// FindNearestNeighbor((CClusterNode*)m_oaTopClusters[z]);
230
231
FindNearestNeighbor(n);
232
233
TraceNeighbors();
234
}
235
m_pTop = n;
236
// m_pClusterDistribution2DF->AddToBinInt_Fast(m_iMonomers-1,m_fMaxDist-n->m_fDistance);
237
// mprintf("BuildTree done.\n");
238
}
239
240
241
void CClusterAnalysis::FindNearestNeighbor(CClusterNode *n)
242
{
243
int z2, z3, z4;
244
float t;
245
CClusterNode *n2;
246
247
n->m_fNNDistance = 1.0E30f;
248
n->m_iNextNeighbor = -1;
249
250
if (m_bDistCache)
251
{
252
for (z2=0;z2<m_oaTopClusters.GetSize();z2++)
253
{
254
n2 = (CClusterNode*)m_oaTopClusters[z2];
255
if (n->m_iIndex == n2->m_iIndex)
256
continue;
257
for (z3=0;z3<n->m_iaAtoms.GetSize();z3++)
258
{
259
for (z4=0;z4<n2->m_iaAtoms.GetSize();z4++)
260
{
261
// mprintf(" z2=%d, z3=%d, z4=%d, dist=%f.\n",z2,z3,z4,n->m_fNNDistance);
262
t = Distance_Cache(n->m_iaAtoms[z3],n2->m_iaAtoms[z4]);
263
if (t < n->m_fNNDistance)
264
{
265
n->m_fNNDistance = t;
266
n->m_iNextNeighbor = n2->m_iIndex;
267
}
268
}
269
}
270
}
271
} else
272
{
273
for (z2=0;z2<m_oaTopClusters.GetSize();z2++)
274
{
275
n2 = (CClusterNode*)m_oaTopClusters[z2];
276
if (n->m_iIndex == n2->m_iIndex)
277
continue;
278
for (z3=0;z3<n->m_vaAtoms.GetSize();z3++)
279
{
280
for (z4=0;z4<n2->m_vaAtoms.GetSize();z4++)
281
{
282
// mprintf(" z2=%d, z3=%d, z4=%d, dist=%f.\n",z2,z3,z4,n->m_fNNDistance);
283
t = Distance(&n->m_vaAtoms[z3],&n2->m_vaAtoms[z4]);
284
if (t < n->m_fNNDistance)
285
{
286
n->m_fNNDistance = t;
287
n->m_iNextNeighbor = n2->m_iIndex;
288
}
289
}
290
}
291
}
292
}
293
}
294
295
296
void CClusterAnalysis::TraceNeighbors()
297
{
298
int z;
299
CClusterNode *n, *n2;
300
301
for (z=0;z<m_oaTopClusters.GetSize();z++)
302
{
303
n = (CClusterNode*)m_oaTopClusters[z];
304
305
if (n->m_iNextNeighbor == -1)
306
continue;
307
308
n2 = (CClusterNode*)m_oaClusters[n->m_iNextNeighbor];
309
310
// while (n2->m_iParent != -1)
311
// n2 = (CClusterNode*)m_oaClusters[n2->m_iParent];
312
313
while (n2->m_pParent != NULL)
314
n2 = n2->m_pParent;
315
316
n->m_iNextNeighbor = n2->m_iIndex;
317
}
318
}
319
320
321
void CClusterAnalysis::DumpDot(const char *s)
322
{
323
FILE *a;
324
int z;
325
CClusterNode *n;
326
327
a = OpenFileWrite(s,true);
328
329
mfprintf(a,"digraph test {\n");
330
for (z=0;z<m_oaClusters.GetSize();z++)
331
{
332
n = (CClusterNode*)m_oaClusters[z];
333
// if (n->m_laChildren.GetSize() == 0)
334
if (n->m_iChildren[0] == -1)
335
mfprintf(a," n%d [label=\"%d X%f\"];\n",z+1,n->m_iIndex+1,n->m_fPosX);
336
else mfprintf(a," n%d [label=\"X%f Y%.2f\"];\n",z+1,n->m_fPosX,n->m_fDistance);
337
}
338
339
mfprintf(a,"\n");
340
341
REC_DumpDot(a,m_pTop);
342
343
mfprintf(a,"}\n");
344
345
fclose(a);
346
}
347
348
349
void CClusterAnalysis::REC_DumpDot(FILE *a, CClusterNode *n)
350
{
351
/* int z;
352
353
for (z=0;z<n->m_laChildren.GetSize();z++)
354
{
355
mfprintf(a," n%d -> n%d;\n",n->m_iIndex+1,(int)(n->m_laChildren[z]+1));
356
REC_DumpDot(a,(CClusterNode*)m_oaClusters[n->m_laChildren[z]]);
357
}*/
358
359
if (n->m_iChildren[0] != -1)
360
{
361
mfprintf(a," n%d -> n%d;\n",n->m_iIndex+1,(int)(n->m_iChildren[0]+1));
362
REC_DumpDot(a,(CClusterNode*)m_oaClusters[n->m_iChildren[0]]);
363
mfprintf(a," n%d -> n%d;\n",n->m_iIndex+1,(int)(n->m_iChildren[1]+1));
364
REC_DumpDot(a,(CClusterNode*)m_oaClusters[n->m_iChildren[1]]);
365
}
366
}
367
368
369
int CA_compareX(const void *arg1, const void *arg2 )
370
{
371
if ((*((CClusterNode**)arg1))->m_fPosX > (*((CClusterNode**)arg2))->m_fPosX)
372
return 1;
373
else return -1;
374
}
375
376
377
void CClusterAnalysis::DumpAgr(const char *s)
378
{
379
int z;
380
// FILE *a;
381
CxFloatArray **fa;
382
// CClusterNode *n, *n2;
383
int *p;
384
float /*d, */d2/*, d2l*/;
385
bool b, c;
386
int i;
387
CGrace *g;
388
char buf[256];
389
390
REC_SortX(1,0,m_pTop);
391
392
qsort(&m_oaBaseClusters[0],m_oaBaseClusters.GetSize(),sizeof(CxObject*),CA_compareX);
393
394
for (z=0;z<m_oaBaseClusters.GetSize();z++)
395
{
396
// mprintf("%d: %f\n",z+1,((CClusterNode*)m_oaBaseClusters[z])->m_fPosX);
397
((CClusterNode*)m_oaBaseClusters[z])->m_fPosX = z+1;
398
}
399
400
REC_AvgX(m_pTop);
401
402
try { g = new CGrace(); } catch(...) { g = NULL; }
403
if (g == NULL) NewException((double)sizeof(CGrace),__FILE__,__LINE__,__PRETTY_FUNCTION__);
404
405
// Revised Version, better and clearer
406
if (m_bPOVTrajectory && m_bPOVDiagram)
407
{
408
g->CurrentGraph()->m_bShowXAxis = false;
409
g->CurrentGraph()->m_bShowFrame = false;
410
g->CurrentGraph()->m_fViewportX1 = 0.10f;
411
g->CurrentGraph()->m_fViewportX2 = 1.30f;
412
g->CurrentGraph()->m_fViewportY1 = 0.6f;
413
g->CurrentGraph()->m_fViewportY2 = 0.95f;
414
g->CurrentGraph()->m_fYAxisBarWidth = 2.0f;
415
g->SetRangeX(-5.0f,m_oaBaseClusters.GetSize()+1);
416
g->SetRangeY(250.0f,300.0f);
417
418
REC_DumpAgr(m_pTop,g,-1,0);
419
} else
420
{
421
try { fa = new CxFloatArray*[m_oaBaseClusters.GetSize()]; } catch(...) { fa = NULL; }
422
if (fa == NULL) NewException((double)m_oaBaseClusters.GetSize()*sizeof(CxFloatArray*),__FILE__,__LINE__,__PRETTY_FUNCTION__);
423
424
try { p = new int[m_oaBaseClusters.GetSize()]; } catch(...) { p = NULL; }
425
if (p == NULL) NewException((double)m_oaBaseClusters.GetSize()*sizeof(int),__FILE__,__LINE__,__PRETTY_FUNCTION__);
426
427
for (z=0;z<m_oaBaseClusters.GetSize();z++)
428
{
429
try { fa[z] = new CxFloatArray("CClusterAnalysis::DumpAgr():fa[z]"); } catch(...) { fa[z] = NULL; }
430
if (fa[z] == NULL) NewException((double)sizeof(CxFloatArray),__FILE__,__LINE__,__PRETTY_FUNCTION__);
431
432
p[z] = 0;
433
}
434
// d = 0;
435
436
// DumpDot("C:\\Software\\test.dot");
437
438
// a = OpenFileWrite(s,true);
439
440
441
for (z=0;z<m_oaBaseClusters.GetSize();z++)
442
{
443
g->AddDataset();
444
g->SetSetLineColor(z,0,0,0);
445
sprintf(buf,"Molecule %d",((CClusterNode*)m_oaBaseClusters[z])->m_iIndex+1);
446
g->SetDatasetName(z,buf);
447
}
448
449
// mfprintf(a,"0");
450
for (z=0;z<m_oaBaseClusters.GetSize();z++)
451
{
452
fa[z]->Add(((CClusterNode*)m_oaBaseClusters[z])->m_fPosX);
453
fa[z]->Add(0);
454
REC_DumpPoints(fa[z],(CClusterNode*)m_oaBaseClusters[z]);
455
// mfprintf(a,"; %f",((CClusterNode*)m_oaBaseClusters[z])->m_fPosX);
456
g->AddXYTupel(z,((CClusterNode*)m_oaBaseClusters[z])->m_fPosX,0);
457
}
458
// mfprintf(a,"\n");
459
460
/* for (z=0;z<fa[0]->GetSize()/2;z++)
461
mfprintf(a,"%f; %f\n",fa[0]->GetAt(z*2),fa[0]->GetAt(z*2+1));
462
fclose(a);
463
return;*/
464
465
// d2 = 0; // ?????
466
467
i = 0;
468
while (true)
469
{
470
i++;
471
// mprintf("Stage %d starting.\n",i);
472
// d2l = d2;
473
d2 = 1.0e35f;
474
for (z=0;z<m_oaBaseClusters.GetSize();z++)
475
{
476
// mprintf(" Cluster %d: %d Pairs, Pos=%d.\n",z+1,fa[z]->GetSize()/2,p[z]);
477
if ((p[z]+1)*2+1 < fa[z]->GetSize())
478
{
479
// mprintf(" %f < %f?\n",fa[z]->GetAt((p[z]+1)*2+1),d2);
480
if (fa[z]->GetAt((p[z]+1)*2+1) < d2)
481
d2 = fa[z]->GetAt((p[z]+1)*2+1);
482
}
483
}
484
// mprintf("d2=%f.\n",d2);
485
if (d2 > 1.0e30)
486
break;
487
_next:
488
b = false;
489
c = false;
490
for (z=0;z<m_oaBaseClusters.GetSize();z++)
491
{
492
if ((p[z]+1)*2+1 < fa[z]->GetSize())
493
{
494
if (fa[z]->GetAt((p[z]+1)*2+1) <= d2)
495
{
496
b = true;
497
p[z]++;
498
}
499
if (!c)
500
{
501
c = true;
502
// mfprintf(a,"%f",d2);
503
}
504
// mfprintf(a,"; %f",fa[z]->GetAt(p[z]*2));
505
g->AddXYTupel(z,fa[z]->GetAt(p[z]*2),d2);
506
}
507
}
508
if (!c)
509
{
510
// mfprintf(a,"%f",d2+1.0);
511
for (z=0;z<m_oaBaseClusters.GetSize();z++)
512
{
513
g->AddXYTupel(z,fa[0]->GetAt(p[0]*2),d2+100.0);
514
// mfprintf(a,"; %f",fa[0]->GetAt(p[0]*2));
515
}
516
// mfprintf(a,"\n");
517
break;
518
}
519
// mfprintf(a,"\n");
520
if (b)
521
goto _next;
522
}
523
524
g->SetRangeX(0,m_oaBaseClusters.GetSize()+1);
525
g->SetRangeY(0,m_fMaxDist);
526
}
527
528
g->MakeTicks();
529
g->SetLabelX("Molecule");
530
531
if (m_bVoroMetric)
532
g->SetLabelY("Cutoff Metric value [nm^-2]");
533
else g->SetLabelY("Cutoff Distance [pm]");
534
535
g->WriteAgr(s,true);
536
537
// fclose(a);
538
}
539
540
541
void CClusterAnalysis::REC_DumpAgr(CClusterNode *n, CGrace *g, int ec, unsigned long color)
542
{
543
int z;
544
CClusterNode *n2;
545
546
if (n->m_iMonomers == 1)
547
return;
548
549
if (ec == -1)
550
{
551
for (z=0;z<m_oaPOVExtendedCluster.GetSize();z++)
552
{
553
if (((CExtendedCluster*)m_oaPOVExtendedCluster[z])->m_iClusterIndex == n->m_iIndex)
554
{
555
ec = z;
556
break;
557
}
558
}
559
}
560
561
if (ec != -1)
562
{
563
for (z=0;z<m_iaPOVSMCluster.GetSize();z++)
564
{
565
if (m_iaPOVSMCluster[z] == ec)
566
{
567
// mprintf("z=%d, ec=%d, color=%d.\n",z,ec,m_iaPOVSMColor[z]);
568
color = m_iaPOVSMColor[z];
569
break;
570
}
571
}
572
} else color = 0x808080;
573
574
if (color == 0xFFFF00)
575
color = 0xD0D000;
576
577
g->AddDataset();
578
g->SetSetLineColor((color&0xFF0000)>>16,(color&0xFF00)>>8,color&0xFF);
579
580
if (ec == -1)
581
g->SetSetLineWidth(g->CurrentGraph()->m_oaDatasets.GetSize()-1,2.0);
582
else g->SetSetLineWidth(g->CurrentGraph()->m_oaDatasets.GetSize()-1,2.5);
583
584
// sprintf(buf,"Cluster %d",z+1);
585
// g->SetDatasetName(buf);
586
587
n2 = (CClusterNode*)m_oaClusters[n->m_iChildren[0]];
588
if (n2->m_iMonomers == 1)
589
g->AddXYTupel(n2->m_fPosX,0);
590
else g->AddXYTupel(n2->m_fPosX,n2->m_fDistance);
591
g->AddXYTupel(n2->m_fPosX,n->m_fDistance);
592
593
n2 = (CClusterNode*)m_oaClusters[n->m_iChildren[1]];
594
g->AddXYTupel(n2->m_fPosX,n->m_fDistance);
595
if (n2->m_iMonomers == 1)
596
g->AddXYTupel(n2->m_fPosX,0);
597
else g->AddXYTupel(n2->m_fPosX,n2->m_fDistance);
598
599
if (((CClusterNode*)m_oaClusters[n->m_iChildren[0]])->m_iMonomers > 1)
600
REC_DumpAgr((CClusterNode*)m_oaClusters[n->m_iChildren[0]],g,ec,color);
601
602
if (((CClusterNode*)m_oaClusters[n->m_iChildren[1]])->m_iMonomers > 1)
603
REC_DumpAgr((CClusterNode*)m_oaClusters[n->m_iChildren[1]],g,ec,color);
604
}
605
606
607
void CClusterAnalysis::REC_SortX(double pos, int depth, CClusterNode *n)
608
{
609
n->m_fPosX = pos;
610
611
if (n->m_iChildren[0] == -1)
612
return;
613
614
if (((CClusterNode*)m_oaClusters[n->m_iChildren[0]])->m_fDistance > ((CClusterNode*)m_oaClusters[n->m_iChildren[1]])->m_fDistance)
615
{
616
REC_SortX(pos+pow(0.5,depth+1),depth+1,(CClusterNode*)m_oaClusters[n->m_iChildren[0]]);
617
REC_SortX(pos-pow(0.5,depth+1),depth+1,(CClusterNode*)m_oaClusters[n->m_iChildren[1]]);
618
} else
619
{
620
REC_SortX(pos-pow(0.5,depth+1),depth+1,(CClusterNode*)m_oaClusters[n->m_iChildren[0]]);
621
REC_SortX(pos+pow(0.5,depth+1),depth+1,(CClusterNode*)m_oaClusters[n->m_iChildren[1]]);
622
}
623
}
624
625
626
void CClusterAnalysis::REC_AvgX(CClusterNode *n)
627
{
628
double d;
629
// int z;
630
631
if (n->m_iChildren[0] == -1)
632
return;
633
634
d = 0;
635
636
/* for (z=0;z<n->m_laChildren.GetSize();z++)
637
{
638
REC_AvgX((CClusterNode*)m_oaClusters[n->m_laChildren[z]]);
639
d += ((CClusterNode*)m_oaClusters[n->m_laChildren[z]])->m_fPosX;
640
}*/
641
642
REC_AvgX((CClusterNode*)m_oaClusters[n->m_iChildren[0]]);
643
d += ((CClusterNode*)m_oaClusters[n->m_iChildren[0]])->m_fPosX;
644
REC_AvgX((CClusterNode*)m_oaClusters[n->m_iChildren[1]]);
645
d += ((CClusterNode*)m_oaClusters[n->m_iChildren[1]])->m_fPosX;
646
647
d /= 2.0; //n->m_laChildren.GetSize();
648
// mprintf("PosX=%f.\n",d);
649
n->m_fPosX = d;
650
}
651
652
653
void CClusterAnalysis::REC_DumpPoints(CxFloatArray *fa, CClusterNode *n)
654
{
655
CClusterNode *n2;
656
657
// if (n->m_iParent != -1)
658
if (n->m_pParent != NULL)
659
{
660
// n2 = (CClusterNode*)m_oaClusters[n->m_iParent];
661
n2 = n->m_pParent;
662
fa->Add(n->m_fPosX);
663
fa->Add(n2->m_fDistance);
664
fa->Add(n2->m_fPosX);
665
fa->Add(n2->m_fDistance);
666
REC_DumpPoints(fa,n2);
667
}
668
}
669
670
671
void CClusterAnalysis::BinDistances(CTimeStep *ts)
672
{
673
int z, m, c, z2, i, z3, z0;
674
double tf, mi, ma, ce, tm;
675
CClusterNode *n;
676
CExtendedCluster *ec;
677
bool tb;
678
char buf[256];
679
FILE *a;
680
681
tf = (double)m_iCounter * g_iStride * g_fTimestepLength / 1000.0;
682
683
ce = 0;
684
c = 0;
685
mi = 1e30;
686
ma = 0;
687
688
if (m_bPOVTrajectory)
689
{
690
for (z=0;z<g_oaSingleMolecules.GetSize();z++)
691
{
692
m_iaPOVSMCluster[z] = -1;
693
m_iaPOVSMColor[z] = -1;
694
}
695
}
696
697
for (z=0;z<m_oaClusters.GetSize();z++)
698
{
699
n = (CClusterNode*)m_oaClusters[z];
700
701
if (n->m_fDistance > 1.0E10)
702
{
703
m_pClusterSizeDF->AddToBinInt_Fast(0);
704
m_pClusterDistributionDF->AddToBinInt_Fast(0,n->m_pParent->m_fDistance);
705
if (m_b2DPlots)
706
{
707
m_pClusterSize2D->AddToBin(tf,1.0);
708
m_pClusterDistribution2D->AddToBin(tf,0.0,n->m_pParent->m_fDistance);
709
}
710
continue;
711
}
712
713
m = n->m_iMonomers-1;
714
m_pClusterSizeDF->AddToBinInt_Fast(m);
715
716
if (n->m_pParent != NULL)
717
m_pClusterDistributionDF->AddToBinInt_Fast(m,n->m_pParent->m_fDistance-n->m_fDistance);
718
else m_pClusterDistributionDF->AddToBinInt_Fast(m,m_fMaxDist-n->m_fDistance);
719
720
m_pClusterDistanceDF->AddToBin_Fast(n->m_fDistance);
721
722
if (m_bHetMeasure)
723
{
724
if (n->m_fDistance < mi)
725
mi = n->m_fDistance;
726
if (n->m_fDistance > ma)
727
ma = n->m_fDistance;
728
729
ce += n->m_fDistance;
730
c++;
731
}
732
733
if (m_b2DPlots)
734
{
735
m_pClusterSize2D->AddToBin(tf,sqrt(m+1)*sqrt(m_iMonomers));
736
if (n->m_pParent != NULL)
737
m_pClusterDistribution2D->AddToBin(tf,sqrt(m+1)*sqrt(m_iMonomers),n->m_pParent->m_fDistance-n->m_fDistance);
738
else m_pClusterDistribution2D->AddToBin(tf,sqrt(m+1)*sqrt(m_iMonomers),m_fMaxDist-n->m_fDistance);
739
m_pClusterDistance2D->AddToBin(tf,n->m_fDistance);
740
}
741
742
if (m_bDiffPlot)
743
{
744
if (n->m_pParent != NULL)
745
m_pClusterDiff2D->AddToBin_IntX(n->m_iMonomers-1,n->m_pParent->m_fDistance - n->m_fDistance,1.0);
746
}
747
748
if (m_bClusterTopo)
749
{
750
if ((n->m_iMonomers > 1) && (n->m_iMonomers <= m_iClusterTopoMax))
751
{
752
ec = new CExtendedCluster();
753
ec->m_iCount = 1;
754
ec->m_iCountUndir = 1;
755
ec->m_fSignificance = n->m_pParent->m_fDistance - n->m_fDistance;
756
ec->m_fSignificanceUndir = n->m_pParent->m_fDistance - n->m_fDistance;
757
ec->CreateFrom(n,ts,this,true);
758
759
tb = false;
760
761
if (m_bClusterTopoDrawDirected || m_bClusterTopoDrawAtom)
762
{
763
for (z2=0;z2<ec->m_oaBonds.GetSize();z2++)
764
if (((CxIntArray*)ec->m_oaBonds[z2])->GetSize() == 0)
765
goto _skip1;
766
ec->BuildAtomCodes();
767
if (AddToClusterTopo(ec,&i))
768
tb = true;
769
770
if (m_bClusterTopoTrajectories)
771
{
772
#ifdef TARGET_WINDOWS
773
sprintf(buf,"ClusterTopo\\xtraj_directed_%02d_%04d.xyz",n->m_iMonomers,i+1);
774
#else
775
sprintf(buf,"ClusterTopo/xtraj_directed_%02d_%04d.xyz",n->m_iMonomers,i+1);
776
#endif
777
a = fopen(buf,"a+t");
778
WriteClusterXYZ(n,ts,a);
779
fclose(a);
780
}
781
_skip1:;
782
}
783
784
if (m_bClusterTopoDrawUndirected)
785
{
786
for (z2=0;z2<ec->m_oaRelBondsUndir.GetSize();z2++)
787
{
788
if (((CxIntArray*)ec->m_oaRelBondsUndir[z2])->GetSize() == 0)
789
{
790
// mprintf("Skip A.\n");
791
goto _skip2;
792
}
793
}
794
ec->BuildAtomCodesUndir();
795
if (AddToClusterTopoUndir(ec,&i))
796
tb = true;
797
798
if (m_bClusterTopoTrajectories)
799
{
800
#ifdef TARGET_WINDOWS
801
sprintf(buf,"ClusterTopo\\xtraj_undir_%02d_%04d.xyz",n->m_iMonomers,i+1);
802
#else
803
sprintf(buf,"ClusterTopo/xtraj_undir_%02d_%04d.xyz",n->m_iMonomers,i+1);
804
#endif
805
a = fopen(buf,"a+t");
806
WriteClusterXYZ(n,ts,a);
807
fclose(a);
808
}
809
_skip2:;
810
}
811
812
if (!tb)
813
delete ec;
814
}
815
}
816
817
if (m_bCutClusters)
818
{
819
if (n->m_iMonomers <= m_iCutClusterMaxSize)
820
{
821
if (n->m_pParent != NULL)
822
{
823
if (n->m_pParent->m_fDistance - n->m_fDistance >= m_pCutClusterDifference[n->m_iMonomers-2])
824
{
825
m_pCutClusterCounter[n->m_iMonomers-2]++;
826
827
WriteClusterXYZ(n,ts,m_pCutClusterFiles[n->m_iMonomers-2]);
828
}
829
}
830
}
831
}
832
}
833
834
if (m_bPOVTrajectory)
835
{
836
for (z0=m_iPOVMonomersMax;z0>=m_iPOVMonomersMin;z0--)
837
{
838
for (z=0;z<m_oaClusters.GetSize();z++)
839
{
840
_beg:
841
n = (CClusterNode*)m_oaClusters[z];
842
843
if (n->m_fDistance > 1.0E10)
844
continue;
845
846
if (n->m_iMonomers == z0)
847
{
848
ec = new CExtendedCluster();
849
ec->CreateFrom(n,ts,this,false);
850
851
for (z2=0;z2<ec->m_iaMonomerSM.GetSize();z2++)
852
{
853
if (m_iaPOVSMCluster[ec->m_iaMonomerSM[z2]] != -1)
854
{
855
z++;
856
goto _beg;
857
}
858
}
859
860
if (IsConnected(ec))
861
{
862
/* int colstat[10];
863
864
for (z2=0;z2<10;z2++)
865
colstat[z2] = 0;
866
867
for (z2=0;z2<ec->m_iaMonomerSM.GetSize();z2++)
868
{
869
for (z3=1;z3<100;z3++)
870
{
871
if (((CxIntArray*)m_oaPOVSMColorHistory[ec->m_iaMonomerSM[z2]])->GetAt(z3) != -1)
872
colstat[((CxIntArray*)m_oaPOVSMColorHistory[ec->m_iaMonomerSM[z2]])->GetAt(z3)]++;
873
}
874
}
875
876
z3 = 0;
877
i = -1;
878
879
for (z2=0;z2<10;z2++)
880
{
881
if (colstat[z2] > z3)
882
{
883
z3 = colstat[z2];
884
i = z2;
885
}
886
}*/
887
888
i = -1;
889
for (z3=1;z3<100;z3++)
890
{
891
if (((CxIntArray*)m_oaPOVSMColorHistory[ec->m_iaMonomerSM[0]])->GetAt(z3) == -1)
892
continue;
893
894
for (z2=1;z2<ec->m_iaMonomerSM.GetSize();z2++)
895
if (((CxIntArray*)m_oaPOVSMColorHistory[ec->m_iaMonomerSM[z2]])->GetAt(z3) != ((CxIntArray*)m_oaPOVSMColorHistory[ec->m_iaMonomerSM[0]])->GetAt(z3))
896
goto _not;
897
i = ((CxIntArray*)m_oaPOVSMColorHistory[ec->m_iaMonomerSM[0]])->GetAt(z3);
898
break;
899
_not:;
900
}
901
902
if (i == -1)
903
{
904
i = m_iPOVColCount;
905
m_iPOVColCount++;
906
if (m_iPOVColCount > 9)
907
m_iPOVColCount = 0;
908
}
909
910
for (z2=0;z2<ec->m_iaMonomerSM.GetSize();z2++)
911
((CxIntArray*)m_oaPOVSMColorHistory[ec->m_iaMonomerSM[z2]])->SetAt(0,i);
912
913
for (z2=0;z2<n->m_iaAtoms.GetSize();z2++)
914
{
915
m_iaPOVSMCluster[g_laAtomSMIndex[m_iaAtomList[n->m_iaAtoms[z2]]]] = m_oaPOVExtendedCluster.GetSize();
916
m_iaPOVSMColor[g_laAtomSMIndex[m_iaAtomList[n->m_iaAtoms[z2]]]] = POVColorFunction(i);
917
}
918
919
m_oaPOVExtendedCluster.Add(ec);
920
921
// mprintf("Put");
922
} else delete ec;
923
}
924
}
925
}
926
}
927
928
929
if (m_bHetMeasure)
930
{
931
ce /= c;
932
933
tm = 0;
934
for (z=0;z<m_oaClusters.GetSize();z++)
935
{
936
n = (CClusterNode*)m_oaClusters[z];
937
938
if (n->m_fDistance > 1.0E10)
939
continue;
940
941
tm += pow(n->m_fDistance - ce, 2);
942
}
943
tm /= c;
944
tm = sqrt(tm) / ce;
945
946
m_faMinDist.Add(mi);
947
m_faMaxDist.Add(ma);
948
m_faHetNumber.Add((ma-mi)/(ma+mi));
949
m_faCenter.Add(ce);
950
m_faIntHetNumber.Add(tm);
951
952
tm = 0;
953
for (z=0;z<m_oaClusters.GetSize();z++)
954
{
955
n = (CClusterNode*)m_oaClusters[z];
956
957
if (n->m_fDistance > 1.0E10)
958
continue;
959
960
tm += fabs(n->m_fDistance - ce);
961
}
962
tm /= c;
963
tm = tm / ce;
964
965
m_faIntHetNumber2.Add(tm);
966
}
967
}
968
969
970
void CClusterAnalysis::Parse()
971
{
972
int ti, ti2, z, z2, z3, z4, z5;
973
char buf[256], buf2[256];
974
CAtomGroup *ag;
975
CSingleMolecule *sm;
976
CMolecule *m;
977
978
mprintf(WHITE,">>> Cluster Analysis >>>\n");
979
980
_canextmol:
981
mprintf("\n");
982
sprintf(buf," Which molecule type to use for cluster analysis (");
983
for (z=0;z<g_oaMolecules.GetSize();z++)
984
{
985
sprintf(buf2,"%s=%d",((CMolecule*)g_oaMolecules[z])->m_sName,z+1);
986
strcat(buf,buf2);
987
if (z < g_oaMolecules.GetSize()-1)
988
strcat(buf,", ");
989
}
990
strcat(buf,")? ");
991
992
ti = AskRangeInteger_ND(buf,1,g_oaMolecules.GetSize()) - 1;
993
994
for (z=0;z<m_oaAtomGroups.GetSize();z++)
995
{
996
if (((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_iIndex == ti)
997
{
998
eprintf("This molecule type is already used for cluster analysis.\n");
999
goto _canextmol;
1000
}
1001
}
1002
1003
try { ag = new CAtomGroup(); } catch(...) { ag = NULL; }
1004
if (ag == NULL) NewException((double)sizeof(CAtomGroup),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1005
1006
_caagain:
1007
AskString(" Which atoms to take into account from %s (*=all)? [*] ",buf,"*",((CMolecule*)g_oaMolecules[ti])->m_sName);
1008
if (!ag->ParseAtoms((CMolecule*)g_oaMolecules[ti],buf))
1009
goto _caagain;
1010
1011
m_oaAtomGroups.Add(ag);
1012
1013
mprintf("\n");
1014
if (AskYesNo(" Add another molecule type to the cluster analysis (y/n)? [no] ",false))
1015
goto _canextmol;
1016
1017
for (z=0;z<m_oaAtomGroups.GetSize();z++)
1018
{
1019
ag = (CAtomGroup*)m_oaAtomGroups[z];
1020
m = ag->m_pMolecule;
1021
for (z2=0;z2<m->m_laSingleMolIndex.GetSize();z2++)
1022
{
1023
sm = (CSingleMolecule*)g_oaSingleMolecules[m->m_laSingleMolIndex[z2]];
1024
for (z3=0;z3<ag->m_oaAtoms.GetSize();z3++)
1025
{
1026
for (z4=0;z4<((CxIntArray*)ag->m_oaAtoms[z3])->GetSize();z4++)
1027
{
1028
m_iaAtomList.Add(((CxIntArray*)sm->m_oaAtomOffset[ag->m_baAtomType[z3]])->GetAt(((CxIntArray*)ag->m_oaAtoms[z3])->GetAt(z4)));
1029
}
1030
}
1031
}
1032
}
1033
m_iAtomListSize = m_iaAtomList.GetSize();
1034
1035
mprintf("\n");
1036
m_bCOMTrick = AskYesNo(" Use \"CoM Trick\" (y/n)? [no] ",false);
1037
1038
if (g_bVoro)
1039
m_bVoroMetric = AskYesNo(" Use Voronoi metric instead of distances (y/n)? [no] ",false);
1040
else m_bVoroMetric = false;
1041
1042
m_bSelectBonds = AskYesNo(" Allow only selected cluster bonds (y) or all bonds (n)? [no] ",false);
1043
1044
if (m_bSelectBonds)
1045
{
1046
m_iaAtomBondFrom.SetSize(m_iAtomListSize);
1047
m_iaAtomBondTo.SetSize(m_iAtomListSize);
1048
for (z=0;z<m_iAtomListSize;z++)
1049
{
1050
m_iaAtomBondFrom[z] = 0;
1051
m_iaAtomBondTo[z] = 0;
1052
}
1053
1054
mprintf("\n * Definition of bond donor atoms\n");
1055
for (z=0;z<m_oaAtomGroups.GetSize();z++)
1056
{
1057
_againx1:
1058
AskString(" Which atoms from %s to allow as bond donor (e.g. C1,C5-7, *=all)? [no atoms] ",buf,"",((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_sName);
1059
if (strlen(buf) == 0)
1060
continue;
1061
try { ag = new CAtomGroup(); } catch(...) { ag = NULL; }
1062
if (ag == NULL) NewException((double)sizeof(CAtomGroup),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1063
if (!ag->ParseAtoms(((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule,buf))
1064
goto _againx1;
1065
for (z2=0;z2<ag->m_baAtomType.GetSize();z2++)
1066
{
1067
// mprintf("z2=%d\n",z2);
1068
for (z3=0;z3<((CxIntArray*)ag->m_oaAtoms[z2])->GetSize();z3++)
1069
{
1070
// mprintf("z3=%d\n",z3);
1071
for (z4=0;z4<((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_laSingleMolIndex.GetSize();z4++)
1072
{
1073
// mprintf("z4=%d\n",z4);
1074
ti = ((CxIntArray*)((CSingleMolecule*)g_oaSingleMolecules[((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_laSingleMolIndex[z4]])->m_oaAtomOffset[ag->m_baAtomType[z2]])->GetAt(((CxIntArray*)ag->m_oaAtoms[z2])->GetAt(z3));
1075
// mprintf("ti=%d\n",ti);
1076
for (z5=0;z5<m_iAtomListSize;z5++)
1077
{
1078
// mprintf("z5=%d: %d vs %d\n",z5,m_iaAtomList[z5],ti);
1079
if (m_iaAtomList[z5] == ti)
1080
{
1081
ti = z5;
1082
goto _found1;
1083
}
1084
}
1085
eprintf("Weird Error A.\n");
1086
_found1:
1087
m_iaAtomBondFrom[ti] = 1;
1088
}
1089
}
1090
}
1091
delete ag;
1092
}
1093
1094
mprintf("\n * Definition of bond acceptor atoms\n");
1095
for (z=0;z<m_oaAtomGroups.GetSize();z++)
1096
{
1097
_againx2:
1098
AskString(" Which atoms from %s to allow as bond acceptor (e.g. C1,C5-7, *=all)? [no atoms] ",buf,"",((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_sName);
1099
if (strlen(buf) == 0)
1100
continue;
1101
try { ag = new CAtomGroup(); } catch(...) { ag = NULL; }
1102
if (ag == NULL) NewException((double)sizeof(CAtomGroup),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1103
if (!ag->ParseAtoms(((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule,buf))
1104
goto _againx2;
1105
for (z2=0;z2<ag->m_baAtomType.GetSize();z2++)
1106
{
1107
for (z3=0;z3<((CxIntArray*)ag->m_oaAtoms[z2])->GetSize();z3++)
1108
{
1109
for (z4=0;z4<((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_laSingleMolIndex.GetSize();z4++)
1110
{
1111
ti = ((CxIntArray*)((CSingleMolecule*)g_oaSingleMolecules[((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_laSingleMolIndex[z4]])->m_oaAtomOffset[ag->m_baAtomType[z2]])->GetAt(((CxIntArray*)ag->m_oaAtoms[z2])->GetAt(z3));
1112
for (z5=0;z5<m_iAtomListSize;z5++)
1113
{
1114
if (m_iaAtomList[z5] == ti)
1115
{
1116
ti = z5;
1117
goto _found2;
1118
}
1119
}
1120
eprintf("Weird Error A.\n");
1121
_found2:
1122
m_iaAtomBondTo[ti] = 1;
1123
}
1124
}
1125
}
1126
delete ag;
1127
}
1128
1129
/* mprintf("\nBond donors:\n");
1130
for (z=0;z<m_iAtomListSize;z++)
1131
mprintf(" %d: %d\n",z,m_iaAtomBondFrom[z]);
1132
mprintf("\nBond acceptors:\n");
1133
for (z=0;z<m_iAtomListSize;z++)
1134
mprintf(" %d: %d\n",z,m_iaAtomBondTo[z]);*/
1135
1136
mprintf("\n");
1137
}
1138
1139
m_bIdealGas = AskYesNo(" Use \"ideal gas mode\" (replace all coords by random numbers) (y/n)? [no] ",false);
1140
1141
if (m_bIdealGas)
1142
m_pRandom = new CRandom();
1143
1144
m_fMaxDist = AskFloat(" Enter max. distance for cluster analysis distribution (in pm): [%d] ",(float)HalfBoxSq3(),HalfBoxSq3());
1145
m_iRes = AskUnsignedInteger(" Enter resolution of the cluster distance distribution: [300] ",300);
1146
1147
m_bHetMeasure = AskYesNo(" Write out temporal development of heterogeneity measures (y/n)? [no] ",false);
1148
if (m_bHetMeasure)
1149
{
1150
if (g_fTimestepLength == 0)
1151
{
1152
g_fTimestepLength = AskFloat(" Enter the length of one trajectory time step in fs: [0.5] ",0.5f);
1153
mprintf("\n");
1154
}
1155
}
1156
1157
m_b2DPlots = AskYesNo(" Create 2D plots with temporal development (y/n)? [no] ",false);
1158
if (m_b2DPlots)
1159
{
1160
if (g_fTimestepLength == 0)
1161
g_fTimestepLength = AskFloat(" Enter the length of one trajectory time step in fs: [0.5] ",0.5f);
1162
m_i2DResX = AskUnsignedInteger(" Which resolution to use for the temporal axis of the 2D plots? [150] ",150);
1163
m_i2DResY = AskUnsignedInteger(" Which resolution to use for the ordinate axis of the 2D plots? [150] ",150);
1164
m_i2DStride = AskUnsignedInteger(" Use each n-th time step for temporal axis (stride)? [1] ",1);
1165
mprintf("\n");
1166
}
1167
1168
m_bDiffPlot = AskYesNo(" Create Cluster Distance Difference plot (y/n)? [no] ",false);
1169
if (m_bDiffPlot)
1170
{
1171
m_fDiffInterval = AskFloat(" Enter max. value of difference inverval (in pm): [200] ",200.0);
1172
mprintf("\n");
1173
}
1174
1175
m_bCutClusters = AskYesNo(" Export cluster structures from cluster analysis (y/n)? [no] ",false);
1176
if (m_bCutClusters)
1177
{
1178
m_iCutClusterMaxSize = AskUnsignedInteger(" Enter the largest cluster size to export: [3] ",3);
1179
1180
try { m_pCutClusterDifference = new double[m_iCutClusterMaxSize-1]; } catch(...) { m_pCutClusterDifference = NULL; }
1181
if (m_pCutClusterDifference == NULL) NewException((double)(m_iCutClusterMaxSize-1)*sizeof(double),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1182
1183
try { m_pCutClusterFiles = new FILE*[m_iCutClusterMaxSize-1]; } catch(...) { m_pCutClusterFiles = NULL; }
1184
if (m_pCutClusterFiles == NULL) NewException((double)(m_iCutClusterMaxSize-1)*sizeof(FILE*),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1185
1186
try { m_pCutClusterCounter = new int[m_iCutClusterMaxSize-1]; } catch(...) { m_pCutClusterCounter = NULL; }
1187
if (m_pCutClusterCounter == NULL) NewException((double)(m_iCutClusterMaxSize-1)*sizeof(int),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1188
1189
for (z=0;z<m_iCutClusterMaxSize-1;z++)
1190
{
1191
sprintf(buf,"cluster_%dmer.xyz",z+2);
1192
m_pCutClusterFiles[z] = OpenFileWrite(buf,true);
1193
m_pCutClusterCounter[z] = 0;
1194
}
1195
1196
for (z=0;z<m_iCutClusterMaxSize-1;z++)
1197
m_pCutClusterDifference[z] = AskFloat(" Enter required distance difference for %d-mer (in pm): [50] ",50.0,z+2);
1198
1199
mprintf("\n");
1200
}
1201
1202
m_bAnim = AskYesNo(" Create Cluster Diagram temporal animation (y/n)? [no] ",false);
1203
if (m_bAnim)
1204
{
1205
m_iResX = AskUnsignedInteger(" Enter width (in pixel) of the animation images: [640] ",640);
1206
m_iResY = AskUnsignedInteger(" Enter height (in pixel) of the animation images: [480] ",480);
1207
mprintf("\n");
1208
}
1209
1210
m_bClusterTopo = AskYesNo(" Perform Cluster Topology Analysis (y/n)? [no] ",false);
1211
1212
if (m_bClusterTopo)
1213
{
1214
mprintf("\n");
1215
m_iClusterTopoMax = AskUnsignedInteger(" Analyze Topology of clusters up to how many monomers? [8] ",8);
1216
1217
if (m_bVoroMetric)
1218
m_bTopoVoroMetric = AskYesNo(" Use Voronoi metric also for bond definition (y/n)? [no] ",false);
1219
1220
m_bClusterTopoAllBonds = AskYesNo(" Consider only hydrogen bonds (n), or allow also dispersive bonds (y)? [no] ",false);
1221
1222
m_bUseBondCutoffs = AskYesNo(" Use bond cutoff values (allow for ring clusters, etc.) (y/n)? [yes] ",true);
1223
1224
if (m_bUseBondCutoffs)
1225
{
1226
m_fClusterTopoHBThres = AskFloat(" Enter hydrogen bond threshold distance (pm): [200] ",200.0f);
1227
1228
if (m_bClusterTopoAllBonds)
1229
m_fClusterTopoDispThres = AskFloat(" Enter dispersive bond threshold distance (pm): [200] ",200.0f);
1230
} else
1231
{
1232
m_fClusterTopoHBThres = 0;
1233
m_fClusterTopoDispThres = 0;
1234
}
1235
1236
m_bClusterTopoDrawAtom = AskYesNo(" Draw atom-level graphs (y/n)? [no] ",false);
1237
m_bClusterTopoDrawDirected = AskYesNo(" Draw directed graphs (y/n)? [yes] ",true);
1238
m_bClusterTopoDrawUndirected = AskYesNo(" Draw undirected graphs (y/n)? [yes] ",true);
1239
1240
if (m_oaAtomGroups.GetSize() > 1)
1241
{
1242
mprintf("\n More than 1 type of molecule used for cluster analysis.\n\n");
1243
for (z=0;z<m_oaAtomGroups.GetSize();z++)
1244
{
1245
mprintf(" * Color for molecule %s\n",((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_sName);
1246
1247
if (m_oaAtomGroups.GetSize() > 1)
1248
{
1249
switch(z)
1250
{
1251
case 0: ti2 = 255; break; // Red
1252
case 1: ti2 = 128; break; // Blue
1253
case 2: ti2 = 128; break; // Green
1254
case 3: ti2 = 255; break; // Yellow
1255
default: ti2 = 211;
1256
}
1257
} else ti2 = 211;
1258
ti = AskUnsignedInteger(" Enter red component (0-255): [%d] ",ti2,ti2);
1259
m_iaClusterTopoMolColor.Add(ti);
1260
1261
if (m_oaAtomGroups.GetSize() > 1)
1262
{
1263
switch(z)
1264
{
1265
case 0: ti2 = 128; break; // Red
1266
case 1: ti2 = 128; break; // Blue
1267
case 2: ti2 = 255; break; // Green
1268
case 3: ti2 = 255; break; // Yellow
1269
default: ti2 = 211;
1270
}
1271
} else ti2 = 211;
1272
ti = AskUnsignedInteger(" Enter green component (0-255): [%d] ",ti2,ti2);
1273
m_iaClusterTopoMolColor.Add(ti);
1274
1275
if (m_oaAtomGroups.GetSize() > 1)
1276
{
1277
switch(z)
1278
{
1279
case 0: ti2 = 128; break; // Red
1280
case 1: ti2 = 255; break; // Blue
1281
case 2: ti2 = 128; break; // Green
1282
case 3: ti2 = 0; break; // Yellow
1283
default: ti2 = 211;
1284
}
1285
} else ti2 = 211;
1286
ti = AskUnsignedInteger(" Enter blue component (0-255): [%d] ",ti2,ti2);
1287
m_iaClusterTopoMolColor.Add(ti);
1288
1289
// if (z < m_oaAtomGroups.GetSize() -1)
1290
mprintf("\n");
1291
}
1292
} else
1293
{
1294
m_iaClusterTopoMolColor.Add(211);
1295
m_iaClusterTopoMolColor.Add(211);
1296
m_iaClusterTopoMolColor.Add(211);
1297
}
1298
m_iClusterTopoTries = AskUnsignedInteger(" How many tries for graph drawing? [10] ",10);
1299
m_bClusterTopoTrajectories = AskYesNo(" Write trajectory for each cluster topology (y/n)? [no] ",false);
1300
mprintf("\n");
1301
}
1302
1303
m_bPOVTrajectory = AskYesNo(" Render POV-Ray trajectory of whole box indicating clusters (y/n)? [no] ",false);
1304
if (m_bPOVTrajectory)
1305
{
1306
mprintf("\n");
1307
m_bClusterTopoDrawDirected = true;
1308
m_bClusterTopoDrawUndirected = true;
1309
m_bClusterTopoAllBonds = AskYesNo(" Consider only hydrogen bonds (n), or allow also dispersive bonds (y)? [no] ",false);
1310
m_fClusterTopoHBThres = AskFloat(" Enter hydrogen bond threshold distance (pm): [200] ",200.0f);
1311
if (m_bClusterTopoAllBonds)
1312
m_fClusterTopoDispThres = AskFloat(" Enter dispersive bond threshold distance (pm): [200] ",200.0f);
1313
m_fPOVAngleIncrement = AskFloat(" Rotate camera (degree per time step)? [0.3] ",0.3f) * Pi / 180.0f;
1314
m_iPOVMonomersMin = AskUnsignedInteger(" Minimal monomer number for emphasized clusters? [5] ",5);
1315
m_iPOVMonomersMax = AskUnsignedInteger(" Maximal monomer number for emphasized clusters? [5] ",5);
1316
if (m_bAnim)
1317
m_bPOVDiagram = AskYesNo(" Colorize Cluster Diagram according to POV files (y/n)? [yes] ",true);
1318
mprintf("\n");
1319
}
1320
1321
m_bDistCache = AskYesNo(" Use distance caching (y/n)? [yes] ",true);
1322
1323
mprintf(WHITE,"\n<<< End of Cluster Analysis <<<\n\n");
1324
}
1325
1326
1327
void CClusterAnalysis::Create()
1328
{
1329
int z, z2;//, z3, z4;
1330
char buf[64];
1331
// CAtomGroup *ag;
1332
// CSingleMolecule *sm;
1333
// CMolecule *m;
1334
1335
if (m_bAnim)
1336
{
1337
system("mkdir cluster_agr");
1338
1339
#ifdef TARGET_WINDOWS
1340
m_fAnim = OpenFileWrite("cluster_agr\\render_cluster_anim",true);
1341
#else
1342
m_fAnim = OpenFileWrite("cluster_agr/render_cluster_anim",true);
1343
#endif
1344
}
1345
1346
m_iCounter = 0;
1347
m_iMonomers = 0;
1348
1349
for (z=0;z<m_oaAtomGroups.GetSize();z++)
1350
m_iMonomers += ((CAtomGroup*)m_oaAtomGroups[z])->m_pMolecule->m_laSingleMolIndex.GetSize();
1351
1352
try { m_poaClusterPoly = new CxObArray[m_iMonomers]; } catch(...) { m_poaClusterPoly = NULL; }
1353
if (m_poaClusterPoly == NULL) NewException((double)m_iMonomers*sizeof(CxObArray),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1354
1355
if (m_bDistCache)
1356
{
1357
try { m_pAtomIndex = new int[g_iGesVirtAtomCount]; } catch(...) { m_pAtomIndex = NULL; }
1358
if (m_pAtomIndex == NULL) NewException((double)g_iGesVirtAtomCount*sizeof(int),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1359
1360
for (z=0;z<g_iGesVirtAtomCount;z++)
1361
m_pAtomIndex[z] = -1;
1362
1363
for (z=0;z<m_iaAtomList.GetSize();z++)
1364
m_pAtomIndex[m_iaAtomList[z]] = z;
1365
1366
try { m_pDistCache = new float[m_iaAtomList.GetSize()*m_iaAtomList.GetSize()]; } catch(...) { m_pDistCache = NULL; }
1367
if (m_pDistCache == NULL) NewException((double)m_iaAtomList.GetSize()*m_iaAtomList.GetSize()*sizeof(float),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1368
}
1369
1370
try { m_pClusterDistributionDF = new CDF(); } catch(...) { m_pClusterDistributionDF = NULL; }
1371
if (m_pClusterDistributionDF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1372
1373
m_pClusterDistributionDF->m_iResolution = m_iMonomers;
1374
m_pClusterDistributionDF->m_fMinVal = 1;
1375
m_pClusterDistributionDF->m_fMaxVal = m_iMonomers+1.0;
1376
m_pClusterDistributionDF->SetLabelX("Cluster Size");
1377
m_pClusterDistributionDF->SetLabelY("Percentage");
1378
m_pClusterDistributionDF->m_bLeft = true;
1379
m_pClusterDistributionDF->Create();
1380
1381
if (m_b2DPlots)
1382
{
1383
try { m_pClusterDistribution2D = new C2DF(); } catch(...) { m_pClusterDistribution2D = NULL; }
1384
if (m_pClusterDistribution2D == NULL) NewException((double)sizeof(C2DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1385
1386
m_pClusterDistribution2D->m_iRes[0] = m_i2DResX;
1387
m_pClusterDistribution2D->m_iRes[1] = m_i2DResY;
1388
m_pClusterDistribution2D->m_fMinVal[0] = 0;
1389
m_pClusterDistribution2D->m_fMaxVal[0] = (double)m_i2DResX*m_i2DStride*g_iStride*g_fTimestepLength/1000.0;
1390
m_pClusterDistribution2D->m_fMinVal[1] = 0;
1391
m_pClusterDistribution2D->m_fMaxVal[1] = m_iMonomers;
1392
m_pClusterDistribution2D->SetLabelX("Simulation Time [ps]");
1393
m_pClusterDistribution2D->SetLabelY("Cluster Size");
1394
m_pClusterDistribution2D->m_bContourLines = false;
1395
m_pClusterDistribution2D->m_fPlotExp = 0.1;
1396
m_pClusterDistribution2D->Create();
1397
}
1398
1399
if (m_bDiffPlot)
1400
{
1401
try { m_pClusterDiff2D = new C2DF(); } catch(...) { m_pClusterDiff2D = NULL; }
1402
if (m_pClusterDiff2D == NULL) NewException((double)sizeof(C2DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1403
1404
m_pClusterDiff2D->m_iRes[0] = m_iMonomers;
1405
m_pClusterDiff2D->m_iRes[1] = 150;
1406
m_pClusterDiff2D->m_fMinVal[0] = 0;
1407
m_pClusterDiff2D->m_fMaxVal[0] = m_iMonomers;
1408
m_pClusterDiff2D->m_fMinVal[1] = 0;
1409
m_pClusterDiff2D->m_fMaxVal[1] = m_fDiffInterval;
1410
m_pClusterDiff2D->SetLabelX("Cluster Size");
1411
m_pClusterDiff2D->SetLabelY("Cluster Distance Difference [pm]");
1412
m_pClusterDiff2D->m_bContourLines = false;
1413
m_pClusterDiff2D->m_fPlotExp = 0.1;
1414
m_pClusterDiff2D->Create();
1415
}
1416
1417
/* try { m_pClusterDistribution2DF = new CDF(); } catch(...) { m_pClusterDistribution2DF = NULL; }
1418
if (m_pClusterDistribution2DF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1419
1420
m_pClusterDistribution2DF->m_iResolution = m_iMonomers;
1421
m_pClusterDistribution2DF->m_fMinVal = 1;
1422
m_pClusterDistribution2DF->m_fMaxVal = m_iMonomers*((m_iMonomers+1.0)/m_iMonomers);
1423
m_pClusterDistribution2DF->SetLabelX("Cluster size");
1424
m_pClusterDistribution2DF->SetLabelY("Percentage");
1425
m_pClusterDistribution2DF->m_bLeft = true;
1426
m_pClusterDistribution2DF->Create();*/
1427
1428
try { m_pClusterSizeDF = new CDF(); } catch(...) { m_pClusterSizeDF = NULL; }
1429
if (m_pClusterSizeDF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1430
1431
m_pClusterSizeDF->m_iResolution = m_iMonomers;
1432
m_pClusterSizeDF->m_fMinVal = 1;
1433
m_pClusterSizeDF->m_fMaxVal = m_iMonomers+1.0;
1434
m_pClusterSizeDF->SetLabelX("Cluster Size");
1435
m_pClusterSizeDF->SetLabelY("Percentage");
1436
m_pClusterSizeDF->m_bLeft = true;
1437
m_pClusterSizeDF->Create();
1438
1439
if (m_b2DPlots)
1440
{
1441
try { m_pClusterSize2D = new C2DF(); } catch(...) { m_pClusterSize2D = NULL; }
1442
if (m_pClusterSize2D == NULL) NewException((double)sizeof(C2DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1443
1444
m_pClusterSize2D->m_iRes[0] = m_i2DResX;
1445
m_pClusterSize2D->m_iRes[1] = m_i2DResY;
1446
m_pClusterSize2D->m_fMinVal[0] = 0;
1447
m_pClusterSize2D->m_fMaxVal[0] = (double)m_i2DResX*m_i2DStride*g_iStride*g_fTimestepLength/1000.0;
1448
m_pClusterSize2D->m_fMinVal[1] = 0;
1449
m_pClusterSize2D->m_fMaxVal[1] = m_iMonomers;
1450
m_pClusterSize2D->SetLabelX("Simulation Time [ps]");
1451
m_pClusterSize2D->SetLabelY("Cluster Size");
1452
m_pClusterSize2D->m_bContourLines = false;
1453
m_pClusterSize2D->m_fPlotExp = 0.3;
1454
m_pClusterSize2D->Create();
1455
}
1456
1457
try { m_pClusterDistanceDF = new CDF(); } catch(...) { m_pClusterDistanceDF = NULL; }
1458
if (m_pClusterDistanceDF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1459
1460
m_pClusterDistanceDF->m_iResolution = m_iRes;
1461
m_pClusterDistanceDF->m_fMinVal = 0;
1462
m_pClusterDistanceDF->m_fMaxVal = m_fMaxDist;
1463
m_pClusterDistanceDF->SetLabelX("Cutoff Distance [pm]");
1464
m_pClusterDistanceDF->SetLabelY("Occurrence");
1465
m_pClusterDistanceDF->Create();
1466
1467
if (m_b2DPlots)
1468
{
1469
try { m_pClusterDistance2D = new C2DF(); } catch(...) { m_pClusterDistance2D = NULL; }
1470
if (m_pClusterDistance2D == NULL) NewException((double)sizeof(C2DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1471
1472
m_pClusterDistance2D->m_iRes[0] = m_i2DResX;
1473
m_pClusterDistance2D->m_iRes[1] = m_i2DResY;
1474
m_pClusterDistance2D->m_fMinVal[0] = 0;
1475
m_pClusterDistance2D->m_fMaxVal[0] = (double)m_i2DResX*m_i2DStride*g_iStride*g_fTimestepLength/1000.0;
1476
m_pClusterDistance2D->m_fMinVal[1] = 0;
1477
m_pClusterDistance2D->m_fMaxVal[1] = m_fMaxDist;
1478
m_pClusterDistance2D->SetLabelX("Simulation Time [ps]");
1479
m_pClusterDistance2D->SetLabelY("Cutoff Distance [pm]");
1480
m_pClusterDistance2D->m_bContourLines = false;
1481
m_pClusterDistance2D->m_fPlotExp = 0.15;
1482
m_pClusterDistance2D->Create();
1483
}
1484
1485
try { m_pPolymerDF = new CDF(); } catch(...) { m_pPolymerDF = NULL; }
1486
if (m_pPolymerDF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1487
1488
m_pPolymerDF->m_iResolution = m_iRes;
1489
m_pPolymerDF->m_fMinVal = 0;
1490
m_pPolymerDF->m_fMaxVal = m_fMaxDist;
1491
m_pPolymerDF->SetLabelX("Cutoff Distance [pm]");
1492
m_pPolymerDF->SetLabelY("Occurrence");
1493
m_pPolymerDF->CreateMulti(m_iMonomers);
1494
m_pPolymerDF->SetLabelMulti(0,"Monomer");
1495
1496
for (z=1;z<m_iMonomers;z++)
1497
{
1498
sprintf(buf,"%d-mer",z+1);
1499
m_pPolymerDF->SetLabelMulti(z,buf);
1500
}
1501
1502
try { m_pClusterCountDF = new CDF(); } catch(...) { m_pClusterCountDF = NULL; }
1503
if (m_pClusterCountDF == NULL) NewException((double)sizeof(CDF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1504
1505
m_pClusterCountDF->m_iResolution = m_iRes;
1506
m_pClusterCountDF->m_fMinVal = 0;
1507
m_pClusterCountDF->m_fMaxVal = m_fMaxDist;
1508
m_pClusterCountDF->SetLabelX("Cutoff Distance [pm]");
1509
m_pClusterCountDF->SetLabelY("Number of Clusters");
1510
m_pClusterCountDF->Create();
1511
1512
if (m_b2DPlots)
1513
{
1514
try { m_pClusterCount2D = new C2DF(); } catch(...) { m_pClusterCount2D = NULL; }
1515
if (m_pClusterCount2D == NULL) NewException((double)sizeof(C2DF),__FILE__,__LINE__,__PRETTY_FUNCTION__);
1516
1517
m_pClusterCount2D->m_iRes[0] = m_i2DResX;
1518
m_pClusterCount2D->m_iRes[1] = m_i2DResY;
1519
m_pClusterCount2D->m_fMinVal[0] = 0;
1520
m_pClusterCount2D->m_fMaxVal[0] = (double)m_i2DResX*m_i2DStride*g_iStride*g_fTimestepLength/1000.0;
1521
m_pClusterCount2D->m_fMinVal[1] = 0;
1522
m_pClusterCount2D->m_fMaxVal[1] = m_fMaxDist;
1523
m_pClusterCount2D->SetLabelX("Simulation Time [ps]");
1524
m_pClusterCount2D->SetLabelY("Cutoff Distance [pm]");
1525
m_pClusterCount2D->m_bContourLines = false;
1526
m_pClusterCount2D->m_fPlotExp = 0.5;
1527
m_pClusterCount2D->Create();
1528
}
1529
1530
if (m_bClusterTopo)
1531
{
1532
m_iaClusterTopoSMTemp.SetSize(g_oaSingleMolecules.GetSize());
1533
1534
if (m_bClusterTopoDrawDirected || m_bClusterTopoDrawAtom)
1535
{
1536
m_oaClusterTopo.SetSize(m_iMonomers);
1537
for (z=0;z<m_iMonomers;z++)
1538
m_oaClusterTopo[z] = new CxObArray();
1539
}
1540
if (m_bClusterTopoDrawUndirected)
1541
{
1542
m_oaClusterTopoUndir.SetSize(m_iMonomers);
1543
for (z=0;z<m_iMonomers;z++)
1544
m_oaClusterTopoUndir[z] = new CxObArray();
1545
}
1546
}
1547
1548
if (m_bPOVTrajectory)
1549
{
1550
system("mkdir POV");
1551
m_iaPOVSMCluster.SetSize(g_oaSingleMolecules.GetSize());
1552
m_iaPOVSMColor.SetSize(g_oaSingleMolecules.GetSize());
1553
for (z=0;z<g_oaSingleMolecules.GetSize();z++)
1554
{
1555
m_oaPOVSMColorHistory.Add(new CxIntArray());
1556
((CxIntArray*)m_oaPOVSMColorHistory[z])->SetSize(100);
1557
for (z2=0;z2<100;z2++)
1558
((CxIntArray*)m_oaPOVSMColorHistory[z])->SetAt(z2,-1);
1559
}
1560
1561
#ifdef TARGET_WINDOWS
1562
m_fPOVScript = OpenFileWrite("POV\\pov_render_script.bat",true);
1563
#else
1564
m_fPOVScript = OpenFileWrite("POV/pov_render_script.bat",true);
1565
#endif
1566
1567
m_iPOVFrameCounter = 1;
1568
m_iPOVColCount = 0;
1569
m_fPOVAngle = 0;
1570
}
1571
}
1572
1573
1574
void CClusterAnalysis::Process(CTimeStep *ts)
1575
{
1576
int z0, z, z2, z3, ti;
1577
CMolecule *m;
1578
CSingleMolecule *sm;
1579
CAtomGroup *ag;
1580
char buf[256];
1581
1582
CleanUp();
1583
1584
if (m_bDistCache)
1585
BuildDistCache(ts);
1586
1587
if (m_bPOVTrajectory)
1588
{
1589
ts->FoldMolecules();
1590
for (z=0;z<g_oaSingleMolecules.GetSize();z++)
1591
for (z2=99;z2>0;z2--)
1592
((CxIntArray*)m_oaPOVSMColorHistory[z])->SetAt(z2,((CxIntArray*)m_oaPOVSMColorHistory[z])->GetAt(z2-1));
1593
}
1594
1595
for (z0=0;z0<m_oaAtomGroups.GetSize();z0++)
1596
{
1597
ag = (CAtomGroup*)m_oaAtomGroups[z0];
1598
m = (CMolecule*)g_oaMolecules[ag->m_pMolecule->m_iIndex];
1599
for (z=0;z<m->m_laSingleMolIndex.GetSize();z++)
1600
{
1601
sm = (CSingleMolecule*)g_oaSingleMolecules[m->m_laSingleMolIndex[z]];
1602
AddCluster(ag->m_iAtomGes);
1603
for (z2=0;z2<ag->m_baAtomType.GetSize();z2++)
1604
{
1605
for (z3=0;z3<((CxIntArray*)ag->m_oaAtoms[z2])->GetSize();z3++)
1606
{
1607
ti = ((CxIntArray*)sm->m_oaAtomOffset[ag->m_baAtomType[z2]])->GetAt(((CxIntArray*)ag->m_oaAtoms[z2])->GetAt(z3));
1608
if (m_bDistCache)
1609
AddParticle(m_pAtomIndex[ti]);
1610
else AddParticle(ts->m_vaCoords[ti][0],ts->m_vaCoords[ti][1],ts->m_vaCoords[ti][2]);
1611
}
1612
}
1613
}
1614
}
1615
1616
BuildTree();
1617
BinDistances(ts);
1618
BinPoly();
1619
1620
if (m_iCounter == 0)
1621
{
1622
DumpAgr("cluster_diagram.agr");
1623
}
1624
1625
if (m_bAnim)
1626
{
1627
1628
#ifdef TARGET_WINDOWS
1629
sprintf(buf,"cluster_agr\\cluster_anim_%05lu.agr",g_iSteps/g_iStride);
1630
#else
1631
sprintf(buf,"cluster_agr/cluster_anim_%05lu.agr",g_iSteps/g_iStride);
1632
#endif
1633
1634
DumpAgr(buf);
1635
mfprintf(m_fAnim,"echo 'Printing frame %lu' ; \n",g_iSteps/g_iStride);
1636
sprintf(buf,"cluster_anim_%05lu.agr",g_iSteps/g_iStride);
1637
mfprintf(m_fAnim,"xmgrace %s -batch gracebatch -nosafe -hardcopy ; \n",buf);
1638
mfprintf(m_fAnim,"mv output.png frame%05lu.png ; \n",g_iSteps/g_iStride);
1639
}
1640
1641
if (m_bPOVTrajectory)
1642
{
1643
1644
#ifdef TARGET_WINDOWS
1645
sprintf(buf,"POV\\frame_%06d.pov",(int)g_iSteps);
1646
#else
1647
sprintf(buf,"POV/frame_%06d.pov",(int)g_iSteps);
1648
#endif
1649
1650
RenderStepPOV(ts,buf);
1651
1652
for (z=0;z<m_oaPOVExtendedCluster.GetSize();z++)
1653
delete (CExtendedCluster*)m_oaPOVExtendedCluster[z];
1654
m_oaPOVExtendedCluster.RemoveAll_KeepSize();
1655
1656
#ifdef TARGET_WINDOWS
1657
mfprintf(m_fPOVScript,"pvengine64 /NR /EXIT +a0.3 +w600 +h600 +fn -j0.0 frame_%06d.pov\n",m_iPOVFrameCounter);
1658
#else
1659
mfprintf(m_fPOVScript,"povray +a0.3 +w600 +h600 +fn -j0.0 frame_%06d.pov\n",m_iPOVFrameCounter);
1660
#endif
1661
1662
fflush(m_fPOVScript);
1663
1664
m_iPOVFrameCounter++;
1665
m_fPOVAngle += m_fPOVAngleIncrement;
1666
}
1667
1668
m_iCounter++;
1669
}
1670
1671
1672
void CClusterAnalysis::CleanUp()
1673
{
1674
int z;
1675
1676
for (z=0;z<m_oaClusters.GetSize();z++)
1677
delete (CClusterNode*)m_oaClusters[z];
1678
m_oaClusters.RemoveAll_KeepSize();
1679
1680
m_oaTopClusters.RemoveAll_KeepSize();
1681
1682
m_oaBaseClusters.RemoveAll_KeepSize();
1683
1684
for (z=0;z<m_iMonomers;z++)
1685
m_poaClusterPoly[z].RemoveAll_KeepSize();
1686
}
1687
1688
1689
void CClusterAnalysis::BinPoly()
1690
{
1691
int z, z2, z3, l, u;
1692
double tf;
1693
CClusterNode *n;
1694
1695
tf = (double)m_iCounter * g_iStride * g_fTimestepLength / 1000.0;
1696
1697
for (z=0;z<m_iMonomers;z++)
1698
{
1699
for (z2=0;z2<m_poaClusterPoly[z].GetSize();z2++)
1700
{
1701
n = (CClusterNode*)(m_poaClusterPoly[z])[z2];
1702
1703
if (n->m_iChildren[0] == -1)
1704
l = 0;
1705
else l = (int)(n->m_fDistance / m_fMaxDist * m_iRes);
1706
1707
if (l >= m_iRes)
1708
l = m_iRes-1;
1709
1710
// if (n->m_iParent == -1)
1711
// u = m_iRes-1;
1712
// else u = ((CClusterNode*)m_oaClusters[n->m_iParent])->m_fDistance / m_fMaxDist * m_iRes - 1;
1713
1714
if (n->m_pParent == NULL)
1715
u = m_iRes-1;
1716
else u = (int)(n->m_pParent->m_fDistance / m_fMaxDist * m_iRes - 1);
1717
1718
if (u >= m_iRes)
1719
u = m_iRes-1;
1720
1721
for (z3=l;z3<=u;z3++)
1722
{
1723
m_pPolymerDF->AddToBin_Multi_Int_Fast(z,z3,z+1.0);
1724
m_pClusterCountDF->AddToBinInt_Fast(z3);
1725
}
1726
1727
if (m_b2DPlots)
1728
{
1729
if (n->m_iChildren[0] == -1)
1730
l = 0;
1731
else l = (int)(n->m_fDistance / m_fMaxDist * m_pClusterCount2D->m_iRes[1]);
1732
1733
if (l >= m_pClusterCount2D->m_iRes[1])
1734
l = m_pClusterCount2D->m_iRes[1]-1;
1735
1736
if (n->m_pParent == NULL)
1737
u = m_pClusterCount2D->m_iRes[1]-1;
1738
else u = (int)(n->m_pParent->m_fDistance / m_fMaxDist * m_pClusterCount2D->m_iRes[1] - 1);
1739
1740
if (u >= m_pClusterCount2D->m_iRes[1])
1741
u = m_pClusterCount2D->m_iRes[1]-1;
1742
1743
for (z3=l;z3<=u;z3++)
1744
m_pClusterCount2D->AddToBin(tf,z3*m_fMaxDist/m_pClusterCount2D->m_iRes[1]);
1745
}
1746
}
1747
}
1748
}
1749
1750
1751
void CClusterAnalysis::BuildDistCache(CTimeStep *ts)
1752
{
1753
int z, z2;
1754
CxVector3 *a, *b;
1755
float t;
1756
1757
if (m_bIdealGas)
1758
{
1759
for (z=0;z<m_iAtomListSize;z++)
1760
{
1761
t = m_pRandom->RandomUniform();
1762
ts->m_vaCoords[m_iaAtomList[z]][0] = t*g_fBoxX;
1763
ts->m_vaCoords[m_iaAtomList[z]][1] = m_pRandom->RandomUniform()*g_fBoxY;
1764
ts->m_vaCoords[m_iaAtomList[z]][2] = m_pRandom->RandomUniform()*g_fBoxZ;
1765
}
1766
}
1767
1768
if (m_bCOMTrick)
1769
{
1770
for (z=0;z<m_iAtomListSize;z++)
1771
{
1772
a = &ts->m_vaCoords[MassCenter(m_iaAtomList[z])];
1773
for (z2=z+1;z2<m_iAtomListSize;z2++)
1774
{
1775
if (m_bSelectBonds)
1776
{
1777
if (((m_iaAtomBondFrom[z] != 0) && (m_iaAtomBondTo[z2] != 0)) || ((m_iaAtomBondFrom[z2] != 0) && (m_iaAtomBondTo[z] != 0)))
1778
{
1779
b = &ts->m_vaCoords[m_iaAtomList[z2]];
1780
t = Distance(a,b);
1781
} else t = 10000000.0f;
1782
} else
1783
{
1784
b = &ts->m_vaCoords[MassCenter(m_iaAtomList[z2])];
1785
t = Distance(a,b);
1786
}
1787
m_pDistCache[z*m_iAtomListSize+z2] = t;
1788
m_pDistCache[z2*m_iAtomListSize+z] = t;
1789
}
1790
}
1791
} else if (m_bVoroMetric)
1792
{
1793
g_pVoroWrapper->BuildVoroMetric(ts);
1794
1795
for (z=0;z<m_iAtomListSize;z++)
1796
{
1797
for (z2=z+1;z2<m_iAtomListSize;z2++)
1798
{
1799
if (m_bSelectBonds)
1800
{
1801
if (((m_iaAtomBondFrom[z] != 0) && (m_iaAtomBondTo[z2] != 0)) || ((m_iaAtomBondFrom[z2] != 0) && (m_iaAtomBondTo[z] != 0)))
1802
{
1803
t = g_pVoroWrapper->m_faVoroMetric[m_iaAtomList[z]*g_iGesAtomCount+m_iaAtomList[z2]];
1804
if (t == 0)
1805
t = 1e20f;
1806
} else t = 1e20f;
1807
} else
1808
{
1809
t = g_pVoroWrapper->m_faVoroMetric[m_iaAtomList[z]*g_iGesAtomCount+m_iaAtomList[z2]];
1810
if (t == 0)
1811
t = 1e20f;
1812
}
1813
m_pDistCache[z*m_iAtomListSize+z2] = t;
1814
m_pDistCache[z2*m_iAtomListSize+z] = t;
1815
}
1816
}
1817
} else // Normal distances
1818
{
1819
for (z=0;z<m_iAtomListSize;z++)
1820
{
1821
a = &ts->m_vaCoords[m_iaAtomList[z]];
1822
for (z2=z+1;z2<m_iAtomListSize;z2++)
1823
{
1824
if (m_bSelectBonds)
1825
{
1826
if (((m_iaAtomBondFrom[z] != 0) && (m_iaAtomBondTo[z2] != 0)) || ((m_iaAtomBondFrom[z2] != 0) && (m_iaAtomBondTo[z] != 0)))
1827
{
1828
b = &ts->m_vaCoords[m_iaAtomList[z2]];
1829
t = Distance(a,b);
1830
} else t = 10000000.0f;
1831
} else
1832
{
1833
b = &ts->m_vaCoords[m_iaAtomList[z2]];
1834
t = Distance(a,b);
1835
}
1836
m_pDistCache[z*m_iAtomListSize+z2] = t;
1837
m_pDistCache[z2*m_iAtomListSize+z] = t;
1838
}
1839
}
1840
}
1841
}
1842
1843
1844
void CClusterAnalysis::BuildClusterDistribution()
1845
{
1846
int z, z2;
1847
1848
for (z=0;z<m_iMonomers;z++)
1849
{
1850
// for (z2=0;z2<m_iRes;z2++)
1851
// m_pClusterDistributionDF->m_pBin[z] += m_pPolymerDF->m_pMultiBin[z][z2];
1852
1853
// m_pClusterDistributionDF->m_pBin[z] *= 100.0 / m_iRes / m_iMonomers / m_iCounter;
1854
1855
m_pClusterDistributionDF->m_pBin[z] *= 100.0 * (z+1.0) / m_iMonomers / m_iCounter / m_fMaxDist;
1856
1857
m_pClusterSizeDF->m_pBin[z] *= 100.0 * (z+1.0) / m_iMonomers / m_iCounter;
1858
}
1859
1860
if (m_b2DPlots)
1861
{
1862
for (z=0;z<m_i2DResY;z++)
1863
{
1864
for (z2=0;z2<m_i2DResX;z2++)
1865
{
1866
m_pClusterDistribution2D->m_pBin[z*m_i2DResY+z2] *= 100.0 * (z+1.0)/m_i2DResY / m_fMaxDist;
1867
m_pClusterSize2D->m_pBin[z*m_i2DResY+z2] *= 100.0 * (z+1.0)/m_i2DResY;
1868
}
1869
}
1870
m_pClusterCount2D->MultiplyBin(1.0/m_i2DStride);
1871
}
1872
}
1873
1874
1875
void CClusterAnalysis::WriteOutput(const char *multibuf)
1876
{
1877
int z, z2, z3, i, j;
1878
CExtendedCluster *ec;
1879
char buf[256], buf2[256];
1880
FILE *a;
1881
1882
if (m_bPOVTrajectory)
1883
{
1884
fclose(m_fPOVScript);
1885
1886
#ifdef TARGET_WINDOWS
1887
a = OpenFileWrite("POV\\video.avs",true);
1888
#else
1889
a = OpenFileWrite("POV/video.avs",true);
1890
#endif
1891
1892
mfprintf(a,"video = ImageSource(\"frame_%c06d.png\", start=1, end=%d, fps=15, pixel_type = \"rgb24\" )\n",'%',m_iPOVFrameCounter-1);
1893
mfprintf(a,"return video\n");
1894
fclose(a);
1895
}
1896
1897
if (m_bClusterTopo)
1898
{
1899
mprintf(WHITE," * Cluster Topology Analysis\n");
1900
system("mkdir ClusterTopo");
1901
1902
if (m_bClusterTopoDrawDirected || m_bClusterTopoDrawAtom)
1903
{
1904
mprintf(" * Directed Graphs\n");
1905
a = OpenFileWrite("cluster_topologies_directed.csv",true);
1906
mfprintf(a,"n-mer; id; count; countperc; avgsigni [pm]; signisum [pm]\n");
1907
for (z=0;z<m_oaClusterTopo.GetSize();z++)
1908
{
1909
if (((CxObArray*)m_oaClusterTopo[z])->GetSize() == 0)
1910
continue;
1911
1912
for (z2=0;z2<((CxObArray*)m_oaClusterTopo[z])->GetSize();z2++)
1913
((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_iIndex = z2;
1914
1915
// Stack Sort
1916
for (z2=0;z2<((CxObArray*)m_oaClusterTopo[z])->GetSize();z2++)
1917
{
1918
j = 0;
1919
i = -1;
1920
for (z3=z2;z3<((CxObArray*)m_oaClusterTopo[z])->GetSize();z3++)
1921
{
1922
if (((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z3))->m_iCount > j)
1923
{
1924
j = ((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z3))->m_iCount;
1925
i = z3;
1926
}
1927
}
1928
ec = (CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2);
1929
((CxObArray*)m_oaClusterTopo[z])->SetAt(z2,((CxObArray*)m_oaClusterTopo[z])->GetAt(i));
1930
((CxObArray*)m_oaClusterTopo[z])->SetAt(i,ec);
1931
}
1932
1933
i = 0;
1934
for (z2=0;z2<((CxObArray*)m_oaClusterTopo[z])->GetSize();z2++)
1935
i += ((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_iCount;
1936
1937
if (m_bClusterTopoTrajectories)
1938
{
1939
for (z2=0;z2<((CxObArray*)m_oaClusterTopo[z])->GetSize();z2++)
1940
{
1941
1942
#ifdef TARGET_WINDOWS
1943
sprintf(buf,"ClusterTopo\\xtraj_directed_%02d_%04d.xyz",z+1,((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_iIndex+1);
1944
sprintf(buf2,"ClusterTopo\\traj_directed_%02d_%04d.xyz",z+1,z2+1);
1945
#else
1946
sprintf(buf,"ClusterTopo/xtraj_directed_%02d_%04d.xyz",z+1,((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_iIndex+1);
1947
sprintf(buf2,"ClusterTopo/traj_directed_%02d_%04d.xyz",z+1,z2+1);
1948
#endif
1949
1950
if (rename(buf,buf2) != 0)
1951
eprintf("Could not rename %s to %s.\n",buf,buf2);
1952
}
1953
}
1954
1955
mprintf(" %d %d-mer topologies found. Rendering:",((CxObArray*)m_oaClusterTopo[z])->GetSize(),z+1);
1956
for (z2=0;z2<((CxObArray*)m_oaClusterTopo[z])->GetSize();z2++)
1957
{
1958
mprintf(".");
1959
if (m_bClusterTopoDrawDirected)
1960
{
1961
1962
#ifdef TARGET_WINDOWS
1963
sprintf(buf,"ClusterTopo\\cluster_directed_%02d_%04d",z+1,z2+1);
1964
#else
1965
sprintf(buf,"ClusterTopo/cluster_directed_%02d_%04d",z+1,z2+1);
1966
#endif
1967
1968
DumpSchematicDirectedDOT((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2),buf);
1969
}
1970
if (m_bClusterTopoDrawAtom)
1971
{
1972
1973
#ifdef TARGET_WINDOWS
1974
sprintf(buf,"ClusterTopo\\cluster_atom_%02d_%04d",z+1,z2+1);
1975
#else
1976
sprintf(buf,"ClusterTopo/cluster_atom_%02d_%04d",z+1,z2+1);
1977
#endif
1978
1979
DumpAtomDOT((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2),buf);
1980
}
1981
mfprintf(a,"%d; %d; %d; %f; %f; %f\n",z+1,z2+1,((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_iCount,100.0*((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_iCount/i,((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_fSignificance/((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_iCount,((CExtendedCluster*)((CxObArray*)m_oaClusterTopo[z])->GetAt(z2))->m_fSignificance);
1982
fflush(a);
1983
// sprintf(buf2,"del %s",buf);
1984
// system(buf2);
1985
}
1986
mfprintf(a,"\n");
1987
mprintf("\n");
1988
}
1989
fclose(a);
1990
}
1991
1992
if (m_bClusterTopoDrawUndirected)
1993
{
1994
mprintf(" * Unirected Graphs\n");
1995
a = OpenFileWrite("cluster_topologies_undirected.csv",true);
1996
mfprintf(a,"n-mer; id; count; countperc; avgsigni [pm]; signisum [pm]\n");
1997
for (z=0;z<m_oaClusterTopoUndir.GetSize();z++)
1998
{
1999
if (((CxObArray*)m_oaClusterTopoUndir[z])->GetSize() == 0)
2000
continue;
2001
2002
for (z2=0;z2<((CxObArray*)m_oaClusterTopoUndir[z])->GetSize();z2++)
2003
((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_iIndex = z2;
2004
2005
// Stack Sort
2006
for (z2=0;z2<((CxObArray*)m_oaClusterTopoUndir[z])->GetSize();z2++)
2007
{
2008
j = 0;
2009
i = -1;
2010
for (z3=z2;z3<((CxObArray*)m_oaClusterTopoUndir[z])->GetSize();z3++)
2011
{
2012
if (((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z3))->m_iCountUndir > j)
2013
{
2014
j = ((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z3))->m_iCountUndir;
2015
i = z3;
2016
}
2017
}
2018
ec = (CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2);
2019
((CxObArray*)m_oaClusterTopoUndir[z])->SetAt(z2,((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(i));
2020
((CxObArray*)m_oaClusterTopoUndir[z])->SetAt(i,ec);
2021
}
2022
2023
i = 0;
2024
for (z2=0;z2<((CxObArray*)m_oaClusterTopoUndir[z])->GetSize();z2++)
2025
i += ((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_iCountUndir;
2026
2027
if (m_bClusterTopoTrajectories)
2028
{
2029
for (z2=0;z2<((CxObArray*)m_oaClusterTopoUndir[z])->GetSize();z2++)
2030
{
2031
#ifdef TARGET_WINDOWS
2032
sprintf(buf,"ClusterTopo\\xtraj_undir_%02d_%04d.xyz",z+1,((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_iIndex+1);
2033
sprintf(buf2,"ClusterTopo\\traj_undir_%02d_%04d.xyz",z+1,z2+1);
2034
#else
2035
sprintf(buf,"ClusterTopo/xtraj_undir_%02d_%04d.xyz",z+1,((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_iIndex+1);
2036
sprintf(buf2,"ClusterTopo/traj_undir_%02d_%04d.xyz",z+1,z2+1);
2037
#endif
2038
2039
if (rename(buf,buf2) != 0)
2040
eprintf("Could not rename %s to %s.\n",buf,buf2);
2041
}
2042
}
2043
2044
mprintf(" %d %d-mer topologies found. Rendering:",((CxObArray*)m_oaClusterTopoUndir[z])->GetSize(),z+1);
2045
for (z2=0;z2<((CxObArray*)m_oaClusterTopoUndir[z])->GetSize();z2++)
2046
{
2047
mprintf(".");
2048
2049
#ifdef TARGET_WINDOWS
2050
sprintf(buf,"ClusterTopo\\cluster_undir_%02d_%04d",z+1,z2+1);
2051
#else
2052
sprintf(buf,"ClusterTopo/cluster_undir_%02d_%04d",z+1,z2+1);
2053
#endif
2054
2055
DumpSchematicUndirectedDOT((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2),buf);
2056
mfprintf(a,"%d; %d; %d; %f; %f; %f\n",z+1,z2+1,((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_iCountUndir,100.0*((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_iCountUndir/i,((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_fSignificanceUndir/((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_iCountUndir,((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_fSignificanceUndir);
2057
// mfprintf(a," %8.4f%c: %d\n",100.0*((CExtendedCluster*)((CxObArray*)m_oaClusterTopoUndir[z])->GetAt(z2))->m_iCountUndir/i,'%',z2+1);
2058
fflush(a);
2059
// sprintf(buf2,"del %s",buf);
2060
// system(buf2);
2061
}
2062
mfprintf(a,"\n");
2063
mprintf("\n");
2064
}
2065
fclose(a);
2066
}
2067
}
2068
2069
mprintf(WHITE," * Cluster Distance Distribution\n");
2070
// mprintf(" %.0f bin entries, %.0f out of bin range (%.2f percent).\n",m_pClusterDistanceDF->m_fBinEntries,m_pClusterDistanceDF->m_fSkipEntries,ZeroDivide(m_pClusterDistanceDF->m_fSkipEntries,m_pClusterDistanceDF->m_fBinEntries+m_pClusterDistanceDF->m_fSkipEntries)*100.0);
2071
// m_pClusterDistanceDF->CalcMeanSD();
2072
m_pClusterDistanceDF->MultiplyBin(100.0/m_iCounter);
2073
// mprintf(" Mean value: %.10G pm Standard deviation: %.10G pm\n",m_pClusterDistanceDF->m_fMean,m_pClusterDistanceDF->m_fSD);
2074
// mprintf(" Min. value: %.10G pm Max. value: %.10G pm\n",m_pClusterDistanceDF->m_fMinInput,m_pClusterDistanceDF->m_fMaxInput);
2075
sprintf(buf,"cluster_distance_df%s.csv",multibuf);
2076
mprintf(" Saving Cluster Distance distribution as %s ...\n",buf);
2077
m_pClusterDistanceDF->Write("",buf,"",false);
2078
sprintf(buf,"cluster_distance_df%s.agr",multibuf);
2079
mprintf(" Saving Cluster Distance distribution AGR file as \"%s\"...\n",buf);
2080
m_pClusterDistanceDF->WriteAgr("",buf,"","Cluster Distance distribution",false);
2081
2082
mprintf(WHITE,"\n * Cluster Significance Distribution\n");
2083
sprintf(buf,"cluster_significance_df%s.csv",multibuf);
2084
mprintf(" Saving Cluster Significance distribution as %s ...\n",buf);
2085
m_pClusterDistributionDF->WriteMulti("",buf,"");
2086
sprintf(buf,"cluster_significance_df%s.agr",multibuf);
2087
mprintf(" Saving Cluster Significance distribution AGR file as \"%s\"...\n",buf);
2088
m_pClusterDistributionDF->WriteMultiAgr("",buf,"","Cluster Significance distribution",false);
2089
2090
/* sprintf(buf,"cluster_sizeX_df%s.csv",multibuf);
2091
mprintf(" Saving Cluster size X distribution as %s ...\n",buf);
2092
m_pClusterDistribution2DF->WriteMulti("",buf,"");
2093
sprintf(buf,"cluster_sizeX_df%s.agr",multibuf);
2094
mprintf(" Saving Cluster size X distribution AGR file as \"%s\"...\n",buf);
2095
m_pClusterDistribution2DF->WriteMultiAgr("",buf,"","Cluster size X distribution",false);*/
2096
2097
mprintf(WHITE,"\n * Cluster Size Distribution\n");
2098
sprintf(buf,"cluster_size_df%s.csv",multibuf);
2099
mprintf(" Saving Cluster Size distribution as %s ...\n",buf);
2100
m_pClusterSizeDF->WriteMulti("",buf,"");
2101
sprintf(buf,"cluster_size_df%s.agr",multibuf);
2102
mprintf(" Saving Cluster Size distribution AGR file as \"%s\"...\n",buf);
2103
m_pClusterSizeDF->WriteMultiAgr("",buf,"","Cluster Size distribution",false);
2104
2105
mprintf(WHITE,"\n * Cluster Count Distribution\n");
2106
// mprintf(" %.0f bin entries, %.0f out of bin range (%.2f percent).\n",m_pClusterCountDF->m_fBinEntries,m_pClusterCountDF->m_fSkipEntries,ZeroDivide(m_pClusterCountDF->m_fSkipEntries,m_pClusterCountDF->m_fBinEntries+m_pClusterCountDF->m_fSkipEntries)*100.0);
2107
// m_pClusterCountDF->CalcMeanSD();
2108
m_pClusterCountDF->MultiplyBin(1.0/m_iCounter);
2109
// mprintf(" Mean value: %.10G pm Standard deviation: %.10G pm\n",m_pClusterCountDF->m_fMean,m_pClusterCountDF->m_fSD);
2110
// mprintf(" Min. value: %.10G pm Max. value: %.10G pm\n",m_pClusterCountDF->m_fMinInput,m_pClusterCountDF->m_fMaxInput);
2111
sprintf(buf,"cluster_count_df%s.csv",multibuf);
2112
mprintf(" Saving Cluster Count distribution as %s ...\n",buf);
2113
m_pClusterCountDF->Write("",buf,"",false);
2114
sprintf(buf,"cluster_count_df%s.agr",multibuf);
2115
mprintf(" Saving Cluster Count distribution AGR file as \"%s\"...\n",buf);
2116
m_pClusterCountDF->WriteAgr("",buf,"","Cluster Count distribution",false);
2117
2118
mprintf(WHITE,"\n * Polymer Distribution\n");
2119
// mprintf(" %.0f bin entries, %.0f out of bin range (%.2f percent).\n",m_pPolymerDF->m_fBinEntries,m_pPolymerDF->m_fSkipEntries,ZeroDivide(m_pPolymerDF->m_fSkipEntries,m_pPolymerDF->m_fBinEntries+m_pPolymerDF->m_fSkipEntries)*100.0);
2120
// m_pPolymerDF->CalcMeanSD();
2121
m_pPolymerDF->MultiplyBin(1.0/m_iCounter);
2122
// mprintf(" Mean value: %.10G pm Standard deviation: %.10G pm\n",m_pPolymerDF->m_fMean,m_pPolymerDF->m_fSD);
2123
// mprintf(" Min. value: %.10G pm Max. value: %.10G pm\n",m_pPolymerDF->m_fMinInput,m_pPolymerDF->m_fMaxInput);
2124
sprintf(buf,"cluster_polymer_df%s.csv",multibuf);
2125
mprintf(" Saving Polymer distribution as %s ...\n",buf);
2126
m_pPolymerDF->WriteMulti("",buf,"");
2127
sprintf(buf,"cluster_polymer_df_cumulative%s.csv",multibuf);
2128
mprintf(" Saving cumulative Polymer distribution as %s ...\n",buf);
2129
m_pPolymerDF->WriteMulti_Cumulative("",buf,"");
2130
sprintf(buf,"cluster_polymer_df%s.agr",multibuf);
2131
mprintf(" Saving Polymer distribution AGR file as \"%s\"...\n",buf);
2132
m_pPolymerDF->WriteMultiAgr("",buf,"","Polymer distribution",false);
2133
sprintf(buf,"cluster_polymer_df_cumulative%s.agr",multibuf);
2134
mprintf(" Saving cumulative Polymer distribution AGR file as \"%s\"...\n",buf);
2135
m_pPolymerDF->WriteMultiAgr_Cumulative("",buf,"","Cumulative Polymer distribution",false);
2136
2137
if (m_bHetMeasure)
2138
{
2139
mprintf(WHITE,"\n * Heterogeneity Time Development\n");
2140
sprintf(buf,"cluster_hetero_td%s.csv",multibuf);
2141
mprintf(" Saving Heterogeneity Time Development as %s ...\n",buf);
2142
a = OpenFileWrite(buf,true);
2143
mfprintf(a,"# Time [ps]; Min. Distance; Max. Distance; Center [pm]; Het. Measure; Int. Het. Measure; Int. Het. Measure NoSqr\n");
2144
for (z=0;z<m_faMinDist.GetSize();z++)
2145
mfprintf(a,"%G; %G; %G; %G; %G; %G; %G\n",z*g_iStride*g_fTimestepLength/1000.0,m_faMinDist[z],m_faMaxDist[z],m_faCenter[z],m_faHetNumber[z],m_faIntHetNumber[z],m_faIntHetNumber2[z]);
2146
fclose(a);
2147
}
2148
2149
if (m_b2DPlots)
2150
{
2151
mprintf(WHITE,"\n * 2D Cluster Distance Distribution\n");
2152
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_TRIPLES"))
2153
{
2154
sprintf(buf,"cluster_distance_df%s_triples.csv",multibuf);
2155
mprintf(" Saving 2D Cluster Distance Distribution triples as %s ...\n",buf);
2156
m_pClusterDistance2D->Write("",buf,"");
2157
}
2158
2159
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATRIX"))
2160
{
2161
sprintf(buf,"cluster_distance_df%s_matrix.csv",multibuf);
2162
mprintf(" Saving 2D Cluster Distance Distribution matrix as %s ...\n",buf);
2163
m_pClusterDistance2D->WriteCSV("",buf,"");
2164
}
2165
2166
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATHEMATICA"))
2167
{
2168
sprintf(buf,"cluster_distance_df%s.nb",multibuf);
2169
mprintf(" Saving 2D Cluster Distance Distribution Mathematica notebook as %s ...\n",buf);
2170
m_pClusterDistance2D->WriteMathematicaNb("",buf,"",false);
2171
}
2172
2173
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_GNUPLOT"))
2174
{
2175
sprintf(buf,"cluster_distance_df%s",multibuf);
2176
mprintf(" Saving 2D Cluster Distance Distribution Gnuplot Input as %s.gp ...\n",buf);
2177
m_pClusterDistance2D->WriteGnuplotInput("",buf,"",false);
2178
}
2179
2180
mprintf(WHITE,"\n * 2D Cluster Size Distribution\n");
2181
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_TRIPLES"))
2182
{
2183
sprintf(buf,"cluster_size_df%s_triples.csv",multibuf);
2184
mprintf(" Saving 2D Cluster Size Distribution triples as %s ...\n",buf);
2185
m_pClusterSize2D->Write("",buf,"");
2186
}
2187
2188
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATRIX"))
2189
{
2190
sprintf(buf,"cluster_size_df%s_matrix.csv",multibuf);
2191
mprintf(" Saving 2D Cluster Size Distribution matrix as %s ...\n",buf);
2192
m_pClusterSize2D->WriteCSV("",buf,"");
2193
}
2194
2195
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATHEMATICA"))
2196
{
2197
sprintf(buf,"cluster_size_df%s.nb",multibuf);
2198
mprintf(" Saving 2D Cluster Size Distribution Mathematica notebook as %s ...\n",buf);
2199
m_pClusterSize2D->WriteMathematicaNb("",buf,"",false);
2200
}
2201
2202
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_GNUPLOT"))
2203
{
2204
sprintf(buf,"cluster_size_df%s",multibuf);
2205
mprintf(" Saving 2D Cluster Size Distribution Gnuplot Input as %s.gp ...\n",buf);
2206
m_pClusterSize2D->WriteGnuplotInput("",buf,"",false);
2207
}
2208
2209
mprintf(WHITE,"\n * 2D Cluster Significance Distribution\n");
2210
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_TRIPLES"))
2211
{
2212
sprintf(buf,"cluster_significance_df%s_triples.csv",multibuf);
2213
mprintf(" Saving 2D Cluster Significance Distribution triples as %s ...\n",buf);
2214
m_pClusterDistribution2D->Write("",buf,"");
2215
}
2216
2217
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATRIX"))
2218
{
2219
sprintf(buf,"cluster_significance_df%s_matrix.csv",multibuf);
2220
mprintf(" Saving 2D Cluster Significance Distribution matrix as %s ...\n",buf);
2221
m_pClusterDistribution2D->WriteCSV("",buf,"");
2222
}
2223
2224
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATHEMATICA"))
2225
{
2226
sprintf(buf,"cluster_significance_df%s.nb",multibuf);
2227
mprintf(" Saving 2D Cluster Significance Distribution Mathematica notebook as %s ...\n",buf);
2228
m_pClusterDistribution2D->WriteMathematicaNb("",buf,"",false);
2229
}
2230
2231
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_GNUPLOT"))
2232
{
2233
sprintf(buf,"cluster_significance_df%s",multibuf);
2234
mprintf(" Saving 2D Cluster Significance Distribution Gnuplot Input as %s.gp ...\n",buf);
2235
m_pClusterDistribution2D->WriteGnuplotInput("",buf,"",false);
2236
}
2237
2238
mprintf(WHITE,"\n * 2D Cluster Count Distribution\n");
2239
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_TRIPLES"))
2240
{
2241
sprintf(buf,"cluster_count_df%s_triples.csv",multibuf);
2242
mprintf(" Saving 2D Cluster Count Distribution triples as %s ...\n",buf);
2243
m_pClusterCount2D->Write("",buf,"");
2244
}
2245
2246
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATRIX"))
2247
{
2248
sprintf(buf,"cluster_count_df%s_matrix.csv",multibuf);
2249
mprintf(" Saving 2D Cluster Count Distribution matrix as %s ...\n",buf);
2250
m_pClusterCount2D->WriteCSV("",buf,"");
2251
}
2252
2253
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATHEMATICA"))
2254
{
2255
sprintf(buf,"cluster_count_df%s.nb",multibuf);
2256
mprintf(" Saving 2D Cluster Count Distribution Mathematica notebook as %s ...\n",buf);
2257
m_pClusterCount2D->WriteMathematicaNb("",buf,"",false);
2258
}
2259
2260
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_GNUPLOT"))
2261
{
2262
sprintf(buf,"cluster_count_df%s",multibuf);
2263
mprintf(" Saving 2D Cluster Count Distribution Gnuplot Input as %s.gp ...\n",buf);
2264
m_pClusterCount2D->WriteGnuplotInput("",buf,"",false);
2265
}
2266
}
2267
2268
if (m_bDiffPlot)
2269
{
2270
mprintf(WHITE,"\n * 2D Cluster Distance Difference Plot\n");
2271
2272
m_pClusterDiff2D->NormalizeXCount();
2273
m_pClusterDiff2D->NormalizeBinIntegral(1000000.0);
2274
2275
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_TRIPLES"))
2276
{
2277
sprintf(buf,"cluster_distdiff_df%s_triples.csv",multibuf);
2278
mprintf(" Saving 2D Cluster Distance Difference triples as %s ...\n",buf);
2279
m_pClusterDiff2D->Write("",buf,"");
2280
}
2281
2282
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATRIX"))
2283
{
2284
sprintf(buf,"cluster_distdiff_df%s_matrix.csv",multibuf);
2285
mprintf(" Saving 2D Cluster Distance Difference matrix as %s ...\n",buf);
2286
m_pClusterDiff2D->WriteCSV("",buf,"");
2287
}
2288
2289
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_MATHEMATICA"))
2290
{
2291
sprintf(buf,"cluster_distdiff_df%s.nb",multibuf);
2292
mprintf(" Saving 2D Cluster Distance Difference Mathematica notebook as %s ...\n",buf);
2293
m_pClusterDiff2D->WriteMathematicaNb("",buf,"",false);
2294
}
2295
2296
if (g_pDatabase->GetBool("/PLOT2D/FORMATS/WRITE_GNUPLOT"))
2297
{
2298
sprintf(buf,"cluster_distdiff_df%s",multibuf);
2299
mprintf(" Saving 2D Cluster Distance Difference Gnuplot Input as %s.gp ...\n",buf);
2300
m_pClusterDiff2D->WriteGnuplotInput("",buf,"",false);
2301
}
2302
}
2303
2304
if (m_bCutClusters)
2305
{
2306
mprintf(WHITE,"\n * 2D Cluster Distance Distribution\n\n");
2307
for (z=0;z<m_iCutClusterMaxSize-1;z++)
2308
{
2309
mprintf(" - %d %d-mer clusters written to file cluster_%dmer.xyz\n",m_pCutClusterCounter[z],z+2,z+2);
2310
fclose(m_pCutClusterFiles[z]);
2311
}
2312
}
2313
2314
if (m_bAnim)
2315
{
2316
mprintf("\n");
2317
2318
#ifdef TARGET_WINDOWS
2319
a = OpenFileWrite("cluster_agr\\gracebatch",true);
2320
#else
2321
a = OpenFileWrite("cluster_agr/gracebatch",true);
2322
#endif
2323
2324
mfprintf(a,"PRINT TO \"output.png\"\n");
2325
mfprintf(a,"HARDCOPY DEVICE \"PNG\"\n");
2326
mfprintf(a,"PAGE SIZE %d, %d\n",m_iResX,m_iResY);
2327
mfprintf(a,"DEVICE \"PNG\" FONT ANTIALIASING on\n");
2328
mfprintf(a,"DEVICE \"PNG\" OP \"compression:9\"\n");
2329
mfprintf(a,"PRINT\n");
2330
fclose(a);
2331
mprintf(" Saved batch script as \"render_cluster_anim\".\n");
2332
fclose(m_fAnim);
2333
}
2334
mprintf("\n");
2335
}
2336
2337
2338
CExtendedCluster::CExtendedCluster()
2339
{
2340
}
2341
2342
2343
CExtendedCluster::~CExtendedCluster()
2344
{
2345
int z;
2346
2347
for (z=0;z<m_oaBonds.GetSize();z++)
2348
delete (CxIntArray*)m_oaBonds[z];
2349
2350
for (z=0;z<m_oaRelBonds.GetSize();z++)
2351
delete (CxIntArray*)m_oaRelBonds[z];
2352
2353
for (z=0;z<m_oaRelBondsUndir.GetSize();z++)
2354
delete (CxIntArray*)m_oaRelBondsUndir[z];
2355
}
2356
2357
2358
void CExtendedCluster::CreateFrom(CClusterNode *n, CTimeStep *ts, CClusterAnalysis *ca, bool wrap)
2359
{
2360
int z, z2, z3, ti, ti2;
2361
CxIntArray *ia;
2362
double dist, hbthres, dispthres;
2363
CSingleMolecule *sm;
2364
2365
m_iaAtoms.SetSize(n->m_iaAtoms.GetSize());
2366
for (z=0;z<n->m_iaAtoms.GetSize();z++)
2367
m_iaAtoms[z] = ca->m_iaAtomList[n->m_iaAtoms[z]];
2368
// m_iaAtoms.CopyFrom(&n->m_iaAtoms);
2369
m_oaBonds.SetSize(m_iaAtoms.GetSize());
2370
m_oaRelBonds.SetSize(m_iaAtoms.GetSize());
2371
m_iClusterIndex = n->m_iIndex;
2372
2373
if (ca->m_bUseBondCutoffs)
2374
{
2375
hbthres = ca->m_fClusterTopoHBThres;
2376
if (!ca->m_bCOMTrick)
2377
if (hbthres < n->m_fDistance)
2378
hbthres = n->m_fDistance;
2379
} else hbthres = n->m_fDistance;
2380
2381
if (ca->m_bUseBondCutoffs)
2382
{
2383
dispthres = ca->m_fClusterTopoDispThres;
2384
if (!ca->m_bCOMTrick)
2385
if (dispthres < n->m_fDistance)
2386
dispthres = n->m_fDistance;
2387
} else dispthres = n->m_fDistance;
2388
2389
// mprintf("cluster = %f pm, hbthres = %f pm, dispthres = %f pm.\n",n->m_fDistance,hbthres,dispthres);
2390
2391
m_iMonomers = n->m_iMonomers;
2392
m_iInterBonds = 0;
2393
2394
for (z=0;z<m_iaAtoms.GetSize();z++)
2395
{
2396
ia = new CxIntArray();
2397
m_oaBonds[z] = ia;
2398
ia = new CxIntArray();
2399
m_oaRelBonds[z] = ia;
2400
}
2401
2402
for (z=0;z<m_iaAtoms.GetSize();z++)
2403
{
2404
for (z2=z+1;z2<m_iaAtoms.GetSize();z2++)
2405
{
2406
if (ca->m_bTopoVoroMetric)
2407
{
2408
dist = ca->Distance_Cache(n->m_iaAtoms[z],n->m_iaAtoms[z2]);
2409
} else
2410
{
2411
if (wrap)
2412
dist = FoldedLength(ts->m_vaCoords[m_iaAtoms[z]] - ts->m_vaCoords[m_iaAtoms[z2]]);
2413
else dist = (ts->m_vaCoords[m_iaAtoms[z]] - ts->m_vaCoords[m_iaAtoms[z2]]).GetLength();
2414
}
2415
2416
if (g_laAtomSMIndex[m_iaAtoms[z]] == g_laAtomSMIndex[m_iaAtoms[z2]])
2417
{
2418
sm = (CSingleMolecule*)g_oaSingleMolecules[g_laAtomSMIndex[m_iaAtoms[z]]];
2419
for (z3=0;z3<sm->m_laBonds.GetSize()/2;z3++)
2420
{
2421
if (((sm->m_laBonds[z3*2] == m_iaAtoms[z]) && (sm->m_laBonds[z3*2+1] == m_iaAtoms[z2])) ||
2422
((sm->m_laBonds[z3*2] == m_iaAtoms[z2]) && (sm->m_laBonds[z3*2+1] == m_iaAtoms[z])))
2423
{
2424
((CxIntArray*)m_oaBonds[z])->Add(m_iaAtoms[z2]);
2425
((CxIntArray*)m_oaBonds[z2])->Add(m_iaAtoms[z]);
2426
2427
((CxIntArray*)m_oaRelBonds[z])->Add(z2);
2428
((CxIntArray*)m_oaRelBonds[z2])->Add(z);
2429
}
2430
}
2431
} else
2432
{
2433
if (ca->m_bSelectBonds)
2434
{
2435
// mprintf("%d -- %d, %d -- %d, (%d,%d), (%d,%d)\n",z,z2,n->m_iaAtoms[z],n->m_iaAtoms[z2],ca->m_iaAtomBondFrom[n->m_iaAtoms[z]],ca->m_iaAtomBondTo[n->m_iaAtoms[z2]],ca->m_iaAtomBondFrom[n->m_iaAtoms[z2]],ca->m_iaAtomBondTo[n->m_iaAtoms[z]]);
2436
if (((ca->m_iaAtomBondFrom[n->m_iaAtoms[z]] != 0) && (ca->m_iaAtomBondTo[n->m_iaAtoms[z2]] != 0)) || ((ca->m_iaAtomBondFrom[n->m_iaAtoms[z2]] != 0) && (ca->m_iaAtomBondTo[n->m_iaAtoms[z]] != 0)))
2437
{
2438
// mprintf(" # accepted, dist = %f vs %f\n",dist,hbthres);
2439
if (dist <= hbthres+0.1)
2440
{
2441
m_iInterBonds++;
2442
2443
((CxIntArray*)m_oaBonds[z])->Add(m_iaAtoms[z2]);
2444
((CxIntArray*)m_oaBonds[z2])->Add(m_iaAtoms[z]);
2445
2446
((CxIntArray*)m_oaRelBonds[z])->Add(z2);
2447
((CxIntArray*)m_oaRelBonds[z2])->Add(z);
2448
}
2449
}
2450
} else
2451
{
2452
if (IsHydrogenBond(g_waAtomRealElement[m_iaAtoms[z]],g_waAtomRealElement[m_iaAtoms[z2]]))
2453
{
2454
if (dist <= hbthres+0.1)
2455
{
2456
m_iInterBonds++;
2457
2458
// mprintf("Adding bond A: %f\n",dist);
2459
2460
((CxIntArray*)m_oaBonds[z])->Add(m_iaAtoms[z2]);
2461
((CxIntArray*)m_oaBonds[z2])->Add(m_iaAtoms[z]);
2462
2463
((CxIntArray*)m_oaRelBonds[z])->Add(z2);
2464
((CxIntArray*)m_oaRelBonds[z2])->Add(z);
2465
}
2466
} else if (ca->m_bClusterTopoAllBonds)
2467
{
2468
if (dist <= dispthres+0.1)
2469
{
2470
m_iInterBonds++;
2471
2472
// mprintf("Adding bond B: %f\n",dist);
2473
2474
((CxIntArray*)m_oaBonds[z])->Add(m_iaAtoms[z2]);
2475
((CxIntArray*)m_oaBonds[z2])->Add(m_iaAtoms[z]);
2476
2477
((CxIntArray*)m_oaRelBonds[z])->Add(z2);
2478
((CxIntArray*)m_oaRelBonds[z2])->Add(z);
2479
}
2480
}
2481
}
2482
}
2483
}
2484
}
2485
2486
// mprintf("%d Monomers, %d Bonds.\n",m_iMonomers,m_iInterBonds);
2487
2488
if (ca->m_bClusterTopoDrawUndirected)
2489
{
2490
m_oaRelBondsUndir.SetSize(m_iMonomers);
2491
2492
for (z=0;z<m_iMonomers;z++)
2493
{
2494
ia = new CxIntArray();
2495
m_oaRelBondsUndir[z] = ia;
2496
}
2497
2498
for (z=0;z<m_iaAtoms.GetSize();z++)
2499
{
2500
for (z2=0;z2<m_iaMonomerSM.GetSize();z2++)
2501
if (g_laAtomSMIndex[m_iaAtoms[z]] == m_iaMonomerSM[z2])
2502
goto _found;
2503
m_iaMonomerSM.Add(g_laAtomSMIndex[m_iaAtoms[z]]);
2504
_found:;
2505
}
2506
2507
// mprintf("Building up %d-mer (%d), d=%fpm, %d InterBonds...\n",m_iMonomers,m_iaMonomerSM.GetSize(),n->m_fDistance,m_iInterBonds);
2508
2509
for (z=0;z<m_iaAtoms.GetSize();z++)
2510
{
2511
ti = g_laAtomSMIndex[m_iaAtoms[z]];
2512
for (z2=0;z2<m_iaMonomerSM.GetSize();z2++)
2513
{
2514
if (ti == m_iaMonomerSM[z2])
2515
{
2516
ti = z2;
2517
goto _found2;
2518
}
2519
}
2520
eprintf("Weird error 1.\n");
2521
_found2:
2522
for (z2=0;z2<((CxIntArray*)m_oaRelBonds[z])->GetSize();z2++)
2523
{
2524
/* if (ca->m_bSelectBonds)
2525
{
2526
if (((ca->m_iaAtomBondFrom[n->m_iaAtoms[z]] == 0) || (ca->m_iaAtomBondTo[n->m_iaAtoms[((CxIntArray*)m_oaRelBonds[z])->GetAt(z2)]] == 0)) && ((ca->m_iaAtomBondFrom[n->m_iaAtoms[((CxIntArray*)m_oaRelBonds[z])->GetAt(z2)]] == 0) || (ca->m_iaAtomBondTo[n->m_iaAtoms[z]] == 0)))
2527
goto _have;
2528
}*/
2529
ti2 = g_laAtomSMIndex[((CxIntArray*)m_oaBonds[z])->GetAt(z2)];
2530
if (ti2 != g_laAtomSMIndex[m_iaAtoms[z]])
2531
{
2532
for (z3=0;z3<m_iaMonomerSM.GetSize();z3++)
2533
{
2534
if (ti2 == m_iaMonomerSM[z3])
2535
{
2536
ti2 = z3;
2537
goto _found3;
2538
}
2539
}
2540
eprintf("Weird error 2.\n");
2541
_found3:
2542
for (z3=0;z3<((CxIntArray*)m_oaRelBondsUndir[ti])->GetSize();z3++)
2543
{
2544
if (((CxIntArray*)m_oaRelBondsUndir[ti])->GetAt(z3) == ti2)
2545
goto _have;
2546
}
2547
((CxIntArray*)m_oaRelBondsUndir[ti])->Add(ti2);
2548
((CxIntArray*)m_oaRelBondsUndir[ti2])->Add(ti);
2549
// mprintf(" Adding Bond %d - %d.\n",ti,ti2);
2550
_have:;
2551
}
2552
}
2553
}
2554
}
2555
}
2556
2557
2558
void CClusterAnalysis::DumpAtomDOT(CExtendedCluster *ec, const char *s)
2559
{
2560
int z, z2;
2561
float tf, tf2;
2562
FILE *a;
2563
char buf[256];
2564
2565
sprintf(buf,"%s.dot",s);
2566
a = OpenFileWrite(buf,true);
2567
2568
mfprintf(a,"graph molecule {\n");
2569
mfprintf(a," graph [pack=true,splines=true,overlap=false];\n");
2570
mfprintf(a," node [shape=none,fontsize=16,fontname=\"Arial\",margin=0,fixedsize,height=0.28];\n");
2571
mfprintf(a," edge [style=bold,len=0.70];\n");
2572
2573
for (z=0;z<ec->m_iaAtoms.GetSize();z++)
2574
mfprintf(a," _%d [ label=\"%s\", width=%.2f ];\n",ec->m_iaAtoms[z]+1,((CAtom*)g_oaAtoms[g_waAtomRealElement[ec->m_iaAtoms[z]]])->m_sName,strlen(((CAtom*)g_oaAtoms[g_waAtomRealElement[ec->m_iaAtoms[z]]])->m_sName)*0.17f);
2575
2576
for (z=0;z<ec->m_iaAtoms.GetSize();z++)
2577
{
2578
for (z2=0;z2<((CxIntArray*)ec->m_oaBonds[z])->GetSize();z2++)
2579
{
2580
if (((CxIntArray*)ec->m_oaRelBonds[z])->GetAt(z2) < z)
2581
continue;
2582
2583
if (g_laAtomSMIndex[ec->m_iaAtoms[z]] == g_laAtomSMIndex[((CxIntArray*)ec->m_oaBonds[z])->GetAt(z2)])
2584
{
2585
tf = 2.5f;
2586
tf2 = 0.4f;
2587
} else
2588
{
2589
tf = 1.0f;
2590
tf2 = 0.8f;
2591
}
2592
2593
mfprintf(a," _%d -- _%d [ len=%.3f, weight=%d, penwidth=%f ];\n",ec->m_iaAtoms[z]+1,((CxIntArray*)ec->m_oaBonds[z])->GetAt(z2)+1,tf2,100,tf);
2594
}
2595
}
2596
mfprintf(a,"}\n");
2597
2598
fclose(a);
2599
2600
// sprintf(buf,"dot %s -Tpng -o%s.png -Kneato -Gepsilon=0.000001 -Gns0limit=5000 -Gmclimit=5000 -v -Gstart=rand",s,s);
2601
2602
// mprintf("Command: \"%s\".\n",buf);
2603
2604
RenderFormula(s);
2605
2606
// system(buf);
2607
}
2608
2609
2610
void CClusterAnalysis::DumpSchematicDirectedDOT(CExtendedCluster *ec, const char *s)
2611
{
2612
int z, z2, cr, cg, cb;
2613
FILE *a;
2614
char buf[256];
2615
2616
sprintf(buf,"%s.dot",s);
2617
a = OpenFileWrite(buf,true);
2618
2619
mfprintf(a,"digraph molecule {\n");
2620
2621
mfprintf(a," graph [ pack=true, overlap=false ];\n");
2622
mfprintf(a," node [ label = \"\", penwidth=3.0, shape = circle, height=0.3, style=\"filled\" ];\n");
2623
mfprintf(a," edge [ penwidth=3.0, arrowsize = 1.5];\n\n");
2624
2625
for (z=0;z<m_iaClusterTopoSMTemp.GetSize();z++)
2626
m_iaClusterTopoSMTemp[z] = 0;
2627
2628
for (z=0;z<ec->m_iaAtoms.GetSize();z++)
2629
{
2630
if (m_iaClusterTopoSMTemp[g_laAtomSMIndex[ec->m_iaAtoms[z]]] != 0)
2631
continue;
2632
m_iaClusterTopoSMTemp[g_laAtomSMIndex[ec->m_iaAtoms[z]]] = 1;
2633
for (z2=0;z2<m_oaAtomGroups.GetSize();z2++)
2634
{
2635
if (((CAtomGroup*)m_oaAtomGroups[z2])->m_pMolecule == (CMolecule*)g_oaMolecules[g_waAtomMolIndex[ec->m_iaAtoms[z]]])
2636
{
2637
cr = m_iaClusterTopoMolColor[z2*3];
2638
cg = m_iaClusterTopoMolColor[z2*3+1];
2639
cb = m_iaClusterTopoMolColor[z2*3+2];
2640
goto _found;
2641
}
2642
}
2643
cr = 211;
2644
cg = 211;
2645
cb = 211;
2646
_found:
2647
mfprintf(a," _%d [ fillcolor=\"#%02X%02X%02X\" ];\n",g_laAtomSMIndex[ec->m_iaAtoms[z]],cr,cg,cb);
2648
}
2649
2650
mfprintf(a,"\n");
2651
2652
for (z=0;z<ec->m_iaAtoms.GetSize();z++)
2653
{
2654
for (z2=0;z2<((CxIntArray*)ec->m_oaBonds[z])->GetSize();z2++)
2655
{
2656
if (((CxIntArray*)ec->m_oaRelBonds[z])->GetAt(z2) < z)
2657
continue;
2658
2659
if (g_laAtomSMIndex[ec->m_iaAtoms[z]] == g_laAtomSMIndex[((CxIntArray*)ec->m_oaBonds[z])->GetAt(z2)])
2660
continue;
2661
2662
if (ec->IsHydrogenBond(g_waAtomRealElement[ec->m_iaAtoms[z]],g_waAtomRealElement[((CxIntArray*)ec->m_oaBonds[z])->GetAt(z2)]))
2663
{
2664
if (((CAtom*)g_oaAtoms[g_waAtomRealElement[ec->m_iaAtoms[z]]])->m_pElement->m_fMass > ((CAtom*)g_oaAtoms[g_waAtomRealElement[((CxIntArray*)ec->m_oaBonds[z])->GetAt(z2)]])->m_pElement->m_fMass)
2665
mfprintf(a," _%d -> _%d;\n",g_laAtomSMIndex[((CxIntArray*)ec->m_oaBonds[z])->GetAt(z2)],g_laAtomSMIndex[ec->m_iaAtoms[z]]);
2666
else mfprintf(a," _%d -> _%d;\n",g_laAtomSMIndex[ec->m_iaAtoms[z]],g_laAtomSMIndex[((CxIntArray*)ec->m_oaBonds[z])->GetAt(z2)]);
2667
} else mfprintf(a," _%d -> _%d [penwidth=1.5,arrowsize = 0];\n",g_laAtomSMIndex[((CxIntArray*)ec->m_oaBonds[z])->GetAt(z2)],g_laAtomSMIndex[ec->m_iaAtoms[z]]);
2668
}
2669
}
2670
mfprintf(a,"}\n");
2671
2672
fclose(a);
2673
2674
RenderFormula(s);
2675
2676
// sprintf(buf,"dot %s -Tpng -o%s.png -Kneato -Gepsilon=0.000001 -Gns0limit=5000 -Gmclimit=5000 -v -Gstart=rand",s,s);
2677
2678
// mprintf("Command: \"%s\".\n",buf);
2679
2680
// system(buf);
2681
}
2682
2683
2684
void CClusterAnalysis::DumpSchematicUndirectedDOT(CExtendedCluster *ec, const char *s)
2685
{
2686
int z, z2, cr, cg, cb;
2687
FILE *a;
2688
char buf[256];
2689
2690
sprintf(buf,"%s.dot",s);
2691
a = OpenFileWrite(buf,true);
2692
2693
mfprintf(a,"graph molecule {\n");
2694
2695
mfprintf(a," graph [ pack=true, overlap=false ];\n");
2696
mfprintf(a," node [ label = \"\", penwidth=3.0, shape = circle, height=0.3, style=\"filled\" ];\n");
2697
mfprintf(a," edge [ penwidth=3.0 ];\n\n");
2698
2699
for (z=0;z<ec->m_iaMonomerSM.GetSize();z++)
2700
{
2701
for (z2=0;z2<m_oaAtomGroups.GetSize();z2++)
2702
{
2703
if (((CAtomGroup*)m_oaAtomGroups[z2])->m_pMolecule == (CMolecule*)g_oaMolecules[((CSingleMolecule*)g_oaSingleMolecules[ec->m_iaMonomerSM[z]])->m_iMolType])
2704
{
2705
cr = m_iaClusterTopoMolColor[z2*3];
2706
cg = m_iaClusterTopoMolColor[z2*3+1];
2707
cb = m_iaClusterTopoMolColor[z2*3+2];
2708
goto _found;
2709
}
2710
}
2711
cr = 211;
2712
cg = 211;
2713
cb = 211;
2714
_found:
2715
mfprintf(a," _%d [ fillcolor=\"#%02X%02X%02X\" ];\n",z,cr,cg,cb);
2716
}
2717
2718
mfprintf(a,"\n");
2719
2720
for (z=0;z<ec->m_iaMonomerSM.GetSize();z++)
2721
{
2722
for (z2=0;z2<((CxIntArray*)ec->m_oaRelBondsUndir[z])->GetSize();z2++)
2723
{
2724
if (((CxIntArray*)ec->m_oaRelBondsUndir[z])->GetAt(z2) < z)
2725
continue;
2726
2727
mfprintf(a," _%d -- _%d;\n",z,((CxIntArray*)ec->m_oaRelBondsUndir[z])->GetAt(z2));
2728
}
2729
}
2730
mfprintf(a,"}\n");
2731
2732
fclose(a);
2733
2734
RenderFormula(s);
2735
2736
// sprintf(buf,"dot %s -Tpng -o%s.png -Kneato -Gepsilon=0.000001 -Gns0limit=5000 -Gmclimit=5000 -v -Gstart=rand",s,s);
2737
2738
// mprintf("Command: \"%s\".\n",buf);
2739
2740
// system(buf);
2741
}
2742
2743
2744
void CExtendedCluster::BuildAtomCodes()
2745
{
2746
int z, z2, c1, c2, i;
2747
double ac;
2748
2749
m_faAtomCodes.SetSize(m_iaAtoms.GetSize());
2750
m_faTempAtomCodes.SetSize(m_iaAtoms.GetSize());
2751
m_faCountAC.SetSize(m_iaAtoms.GetSize());
2752
2753
// Die Anfangswerte der AtomCodes: [Ordnungszahl] * 10.0 + [Zahl der Nicht-Wasserstoff-Bindungen]
2754
for (z=0;z<m_iaAtoms.GetSize();z++)
2755
{
2756
m_faAtomCodes[z] = 10.0 * ((CAtom*)g_oaAtoms[g_waAtomRealElement[m_iaAtoms[z]]])->m_pElement->m_fMass;
2757
2758
i = 0;
2759
for (z2=0;z2<((CxIntArray*)m_oaBonds[z])->GetSize();z2++)
2760
{
2761
// Alle Wasserstoff-Atome ueberspringen
2762
if (((CAtom*)g_oaAtoms[g_waAtomRealElement[((CxIntArray*)m_oaBonds[z])->GetAt(z2)]])->m_pElement->m_fMass < 4.5)
2763
continue;
2764
m_faAtomCodes[z]++;
2765
i++;
2766
}
2767
}
2768
2769
i = 0;
2770
do {
2771
c1 = CountDifferentAtomCodes();
2772
2773
// if (g_bVerbose)
2774
// mprintf(WHITE,"\n Cycle %d: %d different atom codes exist.\n\n",i+1,c1);
2775
2776
for (z=0;z<m_iaAtoms.GetSize();z++)
2777
{
2778
m_faTempAtomCodes[z] = m_faAtomCodes[z] * 5.0;
2779
2780
for (z2=0;z2<((CxIntArray*)m_oaRelBonds[z])->GetSize();z2++)
2781
m_faTempAtomCodes[z] += m_faAtomCodes[((CxIntArray*)m_oaRelBonds[z])->GetAt(z2)];
2782
}
2783
for (z=0;z<m_iaAtoms.GetSize();z++)
2784
m_faAtomCodes[z] = m_faTempAtomCodes[z];
2785
2786
c2 = CountDifferentAtomCodes();
2787
i++;
2788
2789
} while (c1 != c2);
2790
2791
// Sortieren mittels StackSort
2792
for (z=0;z<m_iaAtoms.GetSize();z++)
2793
{
2794
ac = -1;
2795
i = -1;
2796
for (z2=z;z2<m_iaAtoms.GetSize();z2++)
2797
{
2798
if (m_faAtomCodes[z2] > ac)
2799
{
2800
ac = m_faAtomCodes[z2];
2801
i = z2;
2802
}
2803
}
2804
if (i != -1)
2805
{
2806
m_faAtomCodes[i] = m_faAtomCodes[z];
2807
m_faAtomCodes[z] = ac;
2808
} else
2809
{
2810
eprintf("CExtendedCluster::BuildAtomCodes(): Weird error.\n");
2811
return;
2812
}
2813
}
2814
2815
/* mprintf("AtomCodes:\n");
2816
for (z=0;z<m_iaAtoms.GetSize();z++)
2817
mprintf(" - %.0f\n",m_faAtomCodes[z]);*/
2818
}
2819
2820
2821
int CExtendedCluster::CountDifferentAtomCodes()
2822
{
2823
int z, z2, i;
2824
2825
i = 0;
2826
for (z=0;z<m_faAtomCodes.GetSize();z++)
2827
{
2828
for (z2=0;z2<i;z2++)
2829
if (m_faCountAC[z2] == m_faAtomCodes[z])
2830
goto _next;
2831
m_faCountAC[i] = m_faAtomCodes[z];
2832
i++;
2833
_next:;
2834
}
2835
return i;
2836
}
2837
2838
2839
void CExtendedCluster::BuildAtomCodesUndir()
2840
{
2841
int z, z2, c1, c2, i;
2842
double ac;
2843
2844
m_faAtomCodesUndir.SetSize(m_iaMonomerSM.GetSize());
2845
m_faTempAtomCodesUndir.SetSize(m_iaMonomerSM.GetSize());
2846
if (m_faCountAC.GetSize() < m_iaMonomerSM.GetSize())
2847
m_faCountAC.SetSize(m_iaMonomerSM.GetSize());
2848
2849
// Die Anfangswerte der AtomCodes: [Ordnungszahl] * 10.0 + [Zahl der Nicht-Wasserstoff-Bindungen]
2850
for (z=0;z<m_iaMonomerSM.GetSize();z++)
2851
{
2852
m_faAtomCodesUndir[z] = 10.0 * ((CMolecule*)g_oaMolecules[((CSingleMolecule*)g_oaSingleMolecules[m_iaMonomerSM[z]])->m_iMolType])->m_fMass;
2853
2854
i = 0;
2855
for (z2=0;z2<((CxIntArray*)m_oaRelBondsUndir[z])->GetSize();z2++)
2856
{
2857
m_faAtomCodesUndir[z]++;
2858
i++;
2859
}
2860
}
2861
2862
i = 0;
2863
do {
2864
c1 = CountDifferentAtomCodesUndir();
2865
2866
// if (g_bVerbose)
2867
// mprintf(WHITE,"\n Cycle %d: %d different atom codes exist.\n\n",i+1,c1);
2868
2869
for (z=0;z<m_iaMonomerSM.GetSize();z++)
2870
{
2871
m_faTempAtomCodesUndir[z] = m_faAtomCodesUndir[z] * 5.0;
2872
2873
for (z2=0;z2<((CxIntArray*)m_oaRelBondsUndir[z])->GetSize();z2++)
2874
m_faTempAtomCodesUndir[z] += m_faAtomCodesUndir[((CxIntArray*)m_oaRelBondsUndir[z])->GetAt(z2)];
2875
}
2876
for (z=0;z<m_iaMonomerSM.GetSize();z++)
2877
m_faAtomCodesUndir[z] = m_faTempAtomCodesUndir[z];
2878
2879
c2 = CountDifferentAtomCodesUndir();
2880
i++;
2881
2882
} while (c1 != c2);
2883
2884
// Sortieren mittels StackSort
2885
for (z=0;z<m_iaMonomerSM.GetSize();z++)
2886
{
2887
ac = -1;
2888
i = -1;
2889
for (z2=z;z2<m_iaMonomerSM.GetSize();z2++)
2890
{
2891
if (m_faAtomCodesUndir[z2] > ac)
2892
{
2893
ac = m_faAtomCodesUndir[z2];
2894
i = z2;
2895
}
2896
}
2897
if (i != -1)
2898
{
2899
m_faAtomCodesUndir[i] = m_faAtomCodesUndir[z];
2900
m_faAtomCodesUndir[z] = ac;
2901
} else
2902
{
2903
eprintf("CExtendedCluster::BuildAtomCodesUndir(): Weird error.\n");
2904
return;
2905
}
2906
}
2907
2908
/* mprintf("AtomCodes:\n");
2909
for (z=0;z<m_iaAtoms.GetSize();z++)
2910
mprintf(" - %.0f\n",m_faAtomCodes[z]);*/
2911
}
2912
2913
2914
int CExtendedCluster::CountDifferentAtomCodesUndir()
2915
{
2916
int z, z2, i;
2917
2918
i = 0;
2919
for (z=0;z<m_faAtomCodesUndir.GetSize();z++)
2920
{
2921
for (z2=0;z2<i;z2++)
2922
if (m_faCountAC[z2] == m_faAtomCodesUndir[z])
2923
goto _next;
2924
m_faCountAC[i] = m_faAtomCodesUndir[z];
2925
i++;
2926
_next:;
2927
}
2928
return i;
2929
}
2930
2931
2932
bool CClusterAnalysis::AddToClusterTopo(CExtendedCluster *ec, int *i)
2933
{
2934
int z, z2;
2935
CxObArray *oa;
2936
CExtendedCluster *ec2;
2937
2938
oa = ((CxObArray*)m_oaClusterTopo[ec->m_iMonomers-1]);
2939
2940
for (z=0;z<oa->GetSize();z++)
2941
{
2942
ec2 = (CExtendedCluster*)oa->GetAt(z);
2943
for (z2=0;z2<ec->m_faAtomCodes.GetSize();z2++)
2944
if (ec->m_faAtomCodes[z2] != ec2->m_faAtomCodes[z2])
2945
goto _diff;
2946
ec2->m_iCount++;
2947
ec2->m_fSignificance += ec->m_fSignificance;
2948
ec2->m_fSignificanceUndir += ec->m_fSignificanceUndir;
2949
*i = z;
2950
goto _found;
2951
_diff:;
2952
}
2953
// mprintf("Add\n");
2954
if (IsConnected(ec))
2955
{
2956
*i = oa->GetSize();
2957
oa->Add(ec);
2958
return true;
2959
}
2960
2961
_found:
2962
return false;
2963
}
2964
2965
2966
bool CClusterAnalysis::AddToClusterTopoUndir(CExtendedCluster *ec, int *i)
2967
{
2968
int z, z2;
2969
CxObArray *oa;
2970
CExtendedCluster *ec2;
2971
2972
oa = ((CxObArray*)m_oaClusterTopoUndir[ec->m_iMonomers-1]);
2973
2974
for (z=0;z<oa->GetSize();z++)
2975
{
2976
ec2 = (CExtendedCluster*)oa->GetAt(z);
2977
for (z2=0;z2<ec->m_faAtomCodesUndir.GetSize();z2++)
2978
if (ec->m_faAtomCodesUndir[z2] != ec2->m_faAtomCodesUndir[z2])
2979
goto _diff;
2980
ec2->m_iCountUndir++;
2981
*i = z;
2982
goto _found;
2983
_diff:;
2984
}
2985
// mprintf("Add\n");
2986
if (IsConnected(ec))
2987
{
2988
*i = oa->GetSize();
2989
oa->Add(ec);
2990
// mprintf("Take.\n");
2991
return true;
2992
}// else mprintf("Skip B.\n");
2993
2994
_found:
2995
// mprintf("Skip C.\n");
2996
return false;
2997
}
2998
2999
3000
void CClusterAnalysis::RenderFormula(const char *s)
3001
{
3002
int z, zm, tries;
3003
char buf[256], buf2[256], *p, *q;
3004
double tf, mi, ma, av;
3005
FILE *a;
3006
3007
tries = m_iClusterTopoTries;
3008
3009
// mprintf("%d Tries:\n",tries);
3010
zm = -1;
3011
av = 0;
3012
ma = 0;
3013
mi = 1e30;
3014
// mprintf(" Command: dot %s.dot -Tpng -o%s.png -Kneato -Gepsilon=0.000001 -Gnslimit=5000 -Gmclimit=5000 -v -Gstart=rand\n",s,s);
3015
// mprintf(" Optimizing (%d tries): ",tries);
3016
// mprintf(WHITE,"[");
3017
for (z=0;z<tries;z++)
3018
{
3019
// mprintf(WHITE,"#");
3020
3021
#ifdef TARGET_WINDOWS
3022
sprintf(buf,"dot %s.dot -Tpng -o%s.%d.png -Kneato -Gepsilon=0.000001 -Gnslimit=5000 -Gmclimit=5000 -v -Gstart=%d > ClusterTopo\\dot%d.log 2>&1",s,s,z,rand(),z);
3023
#else
3024
sprintf(buf,"dot %s.dot -Tpng -o%s.%d.png -Kneato -Gepsilon=0.000001 -Gnslimit=5000 -Gmclimit=5000 -v -Gstart=%d > ClusterTopo/dot%d.log 2>&1",s,s,z,rand(),z);
3025
#endif
3026
3027
// mprintf(" (%s)\n",buf);
3028
3029
system(buf);
3030
3031
#ifdef TARGET_WINDOWS
3032
sprintf(buf,"ClusterTopo\\dot%d.log",z);
3033
#else
3034
sprintf(buf,"ClusterTopo/dot%d.log",z);
3035
#endif
3036
3037
a = fopen(buf,"rt");
3038
if (a == NULL)
3039
{
3040
eprintf("\nRenderFormula(): Error opening %s for reading.\n",buf);
3041
eprintf("It seems that GraphViz is not working (try command \"dot\").\n\n");
3042
return;
3043
}
3044
while (!feof(a))
3045
{
3046
fgets(buf,256,a);
3047
if (strstr(buf,"final") != NULL)
3048
{
3049
p = buf;
3050
// mprintf(" buf=\"%s\".\n",buf);
3051
while (((*p < '0') || (*p > '9')) && (*p != 0))
3052
p++;
3053
// mprintf(" p=\"%s\".\n",p);
3054
if (*p == 0)
3055
continue;
3056
q = p;
3057
// mprintf(" *p = %d.\n",*p);
3058
while (((*p >= '0') && (*p <= '9')) || (*p == '.'))
3059
{
3060
// mprintf(" *p = %d --> p++;\n",*p);
3061
p++;
3062
}
3063
// mprintf(" p2=\"%s\".\n",p);
3064
*p = 0;
3065
// mprintf(" q=\"%s\".\n",q);
3066
tf = atof(q);
3067
// mprintf(" tf = %g.\n",tf);
3068
if (tf > ma)
3069
ma = tf;
3070
if (tf < mi)
3071
{
3072
mi = tf;
3073
zm = z;
3074
}
3075
av += tf;
3076
goto _done;
3077
}
3078
}
3079
eprintf("\nError with GraphViz output. See dot%d.log.\n",z);
3080
eprintf("It seems like GraphViz is not working (try command \"dot\").\n\n");
3081
return;
3082
_done:
3083
fclose(a);
3084
3085
#ifdef TARGET_WINDOWS
3086
sprintf(buf,"ClusterTopo\\dot%d.log",z);
3087
#else
3088
sprintf(buf,"ClusterTopo/dot%d.log",z);
3089
#endif
3090
3091
remove(buf);
3092
}
3093
// mprintf(WHITE,"]\n");
3094
av /= tries;
3095
for (z=0;z<tries;z++)
3096
{
3097
if (z == zm)
3098
continue;
3099
sprintf(buf,"%s.%d.png",s,z);
3100
if (remove(buf) != 0)
3101
eprintf(" Error removing \"%s\".\n",buf);
3102
}
3103
sprintf(buf,"%s.%d.png",s,zm);
3104
sprintf(buf2,"%s.png",s);
3105
remove(buf2);
3106
if (rename(buf,buf2) != 0)
3107
eprintf(" Error renaming \"%s\" to \"%s\".\n",buf,buf2);
3108
sprintf(buf,"%s.dot",s);
3109
remove(buf);
3110
}
3111
3112
3113
bool CExtendedCluster::IsHydrogenBond(int i1, int i2)
3114
{
3115
CAtom *a1, *a2;
3116
3117
a1 = (CAtom*)g_oaAtoms[i1];
3118
a2 = (CAtom*)g_oaAtoms[i2];
3119
3120
if (mystricmp(a1->m_sName,"H") == 0)
3121
{
3122
if (mystricmp(a2->m_sName,"O") == 0)
3123
return true;
3124
if (mystricmp(a2->m_sName,"N") == 0)
3125
return true;
3126
}
3127
3128
if (mystricmp(a2->m_sName,"H") == 0)
3129
{
3130
if (mystricmp(a1->m_sName,"O") == 0)
3131
return true;
3132
if (mystricmp(a1->m_sName,"N") == 0)
3133
return true;
3134
}
3135
3136
return false;
3137
}
3138
3139
3140
bool CClusterAnalysis::IsConnected(CExtendedCluster *ec)
3141
{
3142
int z;
3143
3144
if (m_iaConnectedBuffer.GetSize() < ec->m_iMonomers)
3145
m_iaConnectedBuffer.SetSize(ec->m_iMonomers);
3146
3147
for (z=0;z<ec->m_iMonomers;z++)
3148
m_iaConnectedBuffer[z] = 0;
3149
3150
REC_IsConnected(ec,0);
3151
3152
for (z=0;z<ec->m_iMonomers;z++)
3153
if (m_iaConnectedBuffer[z] == 0)
3154
return false;
3155
3156
return true;
3157
}
3158
3159
3160
void CClusterAnalysis::REC_IsConnected(CExtendedCluster *ec, int i)
3161
{
3162
int z;
3163
3164
m_iaConnectedBuffer[i] = 1;
3165
3166
for (z=0;z<((CxIntArray*)ec->m_oaRelBondsUndir[i])->GetSize();z++)
3167
{
3168
if (m_iaConnectedBuffer[((CxIntArray*)ec->m_oaRelBondsUndir[i])->GetAt(z)] != 0)
3169
continue;
3170
REC_IsConnected(ec,((CxIntArray*)ec->m_oaRelBondsUndir[i])->GetAt(z));
3171
}
3172
}
3173
3174
3175
void CClusterAnalysis::RenderStepPOV(CTimeStep *ts, const char *s)
3176
{
3177
FILE *b;
3178
CMolecule *m;
3179
CSingleMolecule *sm;
3180
CElement *el, *el2;
3181
int z, z2, z3, z4, o, o2;
3182
CxVector3 vec1, vec2, vec3, vec1b, vec2b, vec3b, vecA, vecB, vecC, vecD;
3183
CMolBond *mb;
3184
CxVector3 cam;
3185
CExtendedCluster *ec;
3186
float lb, bleach;
3187
// float cr, cg, cb;
3188
3189
b = OpenFileWrite(s,true);
3190
3191
mfprintf(b,"// Written by TRAVIS\n");
3192
mfprintf(b,"// See http://www.travis-analyzer.de\n\n");
3193
mfprintf(b,"#version 3.6;\n");
3194
mfprintf(b,"\n");
3195
3196
mfprintf(b,"/**** Atoms ****/\n");
3197
3198
mfprintf(b,"#declare atom_draw = true;\n");
3199
mfprintf(b,"#declare atom_r = 0.65;\n");
3200
mfprintf(b,"#declare atom_specular = 0.7;\n");
3201
mfprintf(b,"#declare atom_reflect = 0;\n");
3202
mfprintf(b,"#declare atom_ambient = 0.2;\n");
3203
mfprintf(b,"#declare atom_diffuse = 0.7;\n");
3204
mfprintf(b,"//#declare atom_color = < 1.0, 1.0, 1.0, 0, 0 >;\n");
3205
mfprintf(b,"//#declare atom_trans = 0.7;\n");
3206
3207
mfprintf(b,"\n");
3208
3209
mfprintf(b,"#declare atom_draw_halo1 = true;\n");
3210
mfprintf(b,"#declare atom_r_halo1 = 0.007;\n");
3211
mfprintf(b,"#declare atom_d_halo1 = 0.0125;\n");
3212
mfprintf(b,"#declare atom_color_halo1 = < 0, 0, 0, 0, 0 >;\n");
3213
mfprintf(b,"#declare atom_specular_halo1 = 0;\n");
3214
mfprintf(b,"#declare atom_reflect_halo1 = 0;\n");
3215
mfprintf(b,"#declare atom_ambient_halo1 = 1.0;\n");
3216
mfprintf(b,"#declare atom_diffuse_halo1 = 0;\n");
3217
3218
mfprintf(b,"\n");
3219
3220
mfprintf(b,"/**** Bonds ****/\n");
3221
3222
mfprintf(b,"#declare bond_draw = true;\n");
3223
mfprintf(b,"#declare bond_r = 0.015 * 0.8;\n");
3224
mfprintf(b,"#declare bond_cluster_r = 0.015 * 0.3;\n");
3225
mfprintf(b,"#declare bond_specular = 0.7;\n");
3226
mfprintf(b,"#declare bond_reflect = 0;\n");
3227
mfprintf(b,"#declare bond_ambient = 0.2;\n");
3228
mfprintf(b,"#declare bond_diffuse = 0.7;\n");
3229
3230
mfprintf(b,"\n");
3231
3232
mfprintf(b,"#declare bond_draw_halo1 = true;\n");
3233
mfprintf(b,"#declare bond_r_halo1 = 0.007;\n");
3234
mfprintf(b,"#declare bond_d_halo1 = 0;\n");
3235
mfprintf(b,"#declare bond_color_halo1 = < 0, 0, 0, 0, 0 >;\n");
3236
mfprintf(b,"#declare bond_specular_halo1 = 0;\n");
3237
mfprintf(b,"#declare bond_reflect_halo1 = 0;\n");
3238
mfprintf(b,"#declare bond_ambient_halo1 = 1.0;\n");
3239
mfprintf(b,"#declare bond_diffuse_halo1 = 0;\n");
3240
3241
mfprintf(b,"\n");
3242
3243
mfprintf(b,"#declare bond_draw_stub = true;\n");
3244
mfprintf(b,"#declare bond_r_stub = 0.004;\n");
3245
mfprintf(b,"#declare bond_l_stub = 0.004;\n");
3246
mfprintf(b,"#declare bond_color_stub = < 0, 0, 0, 0, 0 >;\n");
3247
mfprintf(b,"#declare bond_specular_stub = 0;\n");
3248
mfprintf(b,"#declare bond_reflect_stub = 0;\n");
3249
mfprintf(b,"#declare bond_ambient_stub = 1.0;\n");
3250
mfprintf(b,"#declare bond_diffuse_stub = 0;\n");
3251
3252
mfprintf(b,"\n");
3253
3254
bleach = 0.75f;
3255
3256
mfprintf(b,"/**** Element Colors ****/\n");
3257
for (z=0;z<g_oaAtoms.GetSize();z++)
3258
{
3259
if (z == g_iVirtAtomType)
3260
continue;
3261
el = ((CAtom*)g_oaAtoms[z])->m_pElement;
3262
mfprintf(b,"#declare elem_%s_col = < %f, %f, %f, 0, 0 >;\n",el->m_sLabel,el->m_iColorR/255.0*(1.0-bleach)+bleach,el->m_iColorG/255.0*(1.0-bleach)+bleach,el->m_iColorB/255.0*(1.0-bleach)+bleach);
3263
}
3264
3265
mfprintf(b,"\n/**** Element Radii ****/\n");
3266
for (z=0;z<g_oaAtoms.GetSize();z++)
3267
{
3268
if (z == g_iVirtAtomType)
3269
continue;
3270
el = ((CAtom*)g_oaAtoms[z])->m_pElement;
3271
mfprintf(b,"#declare elem_%s_r = %f * 0.8;\n",el->m_sLabel,el->m_fRadius/1000.0);
3272
}
3273
mfprintf(b,"\n");
3274
3275
// cam[0] = 3.0;
3276
// cam[1] = 2.0;
3277
// cam[2] = 8.0;
3278
3279
lb = g_fBoxX;
3280
if (g_fBoxY > lb)
3281
lb = g_fBoxY;
3282
if (g_fBoxZ > lb)
3283
lb = g_fBoxZ;
3284
3285
cam[0] = 8.544 * lb/1566.0 * sin(m_fPOVAngle);
3286
cam[1] = 2.0 * lb/1566.0;
3287
cam[2] = 8.544 * lb/1566.0 * cos(m_fPOVAngle);
3288
3289
mfprintf(b,"global_settings {\n");
3290
mfprintf(b," assumed_gamma 1\n");
3291
mfprintf(b,"/* radiosity {\n");
3292
mfprintf(b," pretrace_start 0.08\n");
3293
mfprintf(b," pretrace_end 0.04\n");
3294
mfprintf(b," count 100\n\n");
3295
mfprintf(b," nearest_count 5\n");
3296
mfprintf(b," error_bound 0.4\n");
3297
mfprintf(b," recursion_limit 1\n\n");
3298
mfprintf(b," low_error_factor .5\n");
3299
mfprintf(b," gray_threshold 0.0\n");
3300
mfprintf(b," minimum_reuse 0.015\n");
3301
mfprintf(b," brightness 1\n\n");
3302
mfprintf(b," adc_bailout 0.01/2\n");
3303
mfprintf(b," }*/\n");
3304
mfprintf(b,"}\n");
3305
3306
mfprintf(b,"\ncamera {\n");
3307
mfprintf(b," location <%f, %f, %f>\n",cam[0],cam[1],cam[2]);
3308
mfprintf(b," sky y\n");
3309
mfprintf(b," right -0.3*x*image_width/image_height\n");
3310
mfprintf(b," up 0.3*y\n");
3311
mfprintf(b," look_at <0, 0, 0>\n");
3312
mfprintf(b,"}\n");
3313
mfprintf(b,"\n");
3314
3315
mfprintf(b,"// Solid background\n");
3316
mfprintf(b,"background { rgb < 1.0, 1.0, 1.0 > }\n");
3317
mfprintf(b,"\n");
3318
3319
mfprintf(b,"// Gradient background\n");
3320
mfprintf(b,"/*sky_sphere {\n");
3321
mfprintf(b," pigment {\n");
3322
mfprintf(b," gradient y\n");
3323
mfprintf(b," color_map {\n");
3324
mfprintf(b," [ 0 color rgb < 0.05, 0.05, 0.05 > ]\n");
3325
mfprintf(b," [ 1 color rgb < 0.20, 0.16, 0.50 > ]\n");
3326
mfprintf(b," }\n");
3327
mfprintf(b," scale 0.1\n");
3328
mfprintf(b," translate -0.05\n");
3329
mfprintf(b," }\n");
3330
mfprintf(b,"}*/\n\n");
3331
3332
mfprintf(b,"/**** Invisible, only for Radiosity ****/\n");
3333
mfprintf(b,"sphere {\n");
3334
mfprintf(b," <0, 0, 0>, 1\n");
3335
mfprintf(b," texture {\n");
3336
mfprintf(b," pigment {color rgb < 1.0, 1.0, 1.0 > }\n");
3337
mfprintf(b," finish { diffuse 0 ambient 1 }\n");
3338
mfprintf(b," }\n");
3339
mfprintf(b," hollow on\n");
3340
mfprintf(b," no_shadow\n");
3341
mfprintf(b," no_image\n");
3342
mfprintf(b," scale 30000\n");
3343
mfprintf(b,"}\n\n");
3344
3345
// mfprintf(b,"light_source { < -8, 20, 20 > color rgb 0.8 }\n");
3346
mfprintf(b,"light_source { < %f, 20, %f > color rgb 0.8 }\n",21.54*sin(m_fPOVAngle-0.3805),21.54*cos(m_fPOVAngle-0.3805));
3347
mfprintf(b,"//light_source { < 25, 12, 20 > color rgb 0.5 }\n\n");
3348
3349
mfprintf(b,"union {\n");
3350
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n",-g_fBoxX/2000.0,-g_fBoxY/2000.0,-g_fBoxZ/2000.0, g_fBoxX/2000.0,-g_fBoxY/2000.0,-g_fBoxZ/2000.0);
3351
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n", g_fBoxX/2000.0,-g_fBoxY/2000.0,-g_fBoxZ/2000.0, g_fBoxX/2000.0, g_fBoxY/2000.0,-g_fBoxZ/2000.0);
3352
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n", g_fBoxX/2000.0, g_fBoxY/2000.0,-g_fBoxZ/2000.0,-g_fBoxX/2000.0, g_fBoxY/2000.0,-g_fBoxZ/2000.0);
3353
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n",-g_fBoxX/2000.0, g_fBoxY/2000.0,-g_fBoxZ/2000.0,-g_fBoxX/2000.0,-g_fBoxY/2000.0,-g_fBoxZ/2000.0);
3354
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n",-g_fBoxX/2000.0,-g_fBoxY/2000.0, g_fBoxZ/2000.0, g_fBoxX/2000.0,-g_fBoxY/2000.0, g_fBoxZ/2000.0);
3355
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n", g_fBoxX/2000.0,-g_fBoxY/2000.0, g_fBoxZ/2000.0, g_fBoxX/2000.0, g_fBoxY/2000.0, g_fBoxZ/2000.0);
3356
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n", g_fBoxX/2000.0, g_fBoxY/2000.0, g_fBoxZ/2000.0,-g_fBoxX/2000.0, g_fBoxY/2000.0, g_fBoxZ/2000.0);
3357
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n",-g_fBoxX/2000.0, g_fBoxY/2000.0, g_fBoxZ/2000.0,-g_fBoxX/2000.0,-g_fBoxY/2000.0, g_fBoxZ/2000.0);
3358
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n",-g_fBoxX/2000.0,-g_fBoxY/2000.0,-g_fBoxZ/2000.0,-g_fBoxX/2000.0,-g_fBoxY/2000.0, g_fBoxZ/2000.0);
3359
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n", g_fBoxX/2000.0,-g_fBoxY/2000.0,-g_fBoxZ/2000.0, g_fBoxX/2000.0,-g_fBoxY/2000.0, g_fBoxZ/2000.0);
3360
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n", g_fBoxX/2000.0, g_fBoxY/2000.0,-g_fBoxZ/2000.0, g_fBoxX/2000.0, g_fBoxY/2000.0, g_fBoxZ/2000.0);
3361
mfprintf(b," cylinder { < %f, %f, %f >, < %f, %f, %f >, 0.01 open }\n",-g_fBoxX/2000.0, g_fBoxY/2000.0,-g_fBoxZ/2000.0,-g_fBoxX/2000.0, g_fBoxY/2000.0, g_fBoxZ/2000.0);
3362
mfprintf(b," no_shadow pigment { rgbft < 0, 0, 1, 0, 0> }\n finish { reflection 0 specular 0.7 ambient 0.2 diffuse 0.7 }\n}\n\n");
3363
3364
mfprintf(b,"#macro m_atom_color(col)\n");
3365
mfprintf(b," #ifdef(atom_color)\n");
3366
mfprintf(b," atom_color\n");
3367
mfprintf(b," #else #if (defined(atom_trans))\n");
3368
mfprintf(b," col + < 0, 0, 0, 0, atom_trans >\n");
3369
mfprintf(b," #else\n");
3370
mfprintf(b," col\n");
3371
mfprintf(b," #end #end\n");
3372
mfprintf(b,"#end\n");
3373
3374
mfprintf(b,"\nunion {\n");
3375
3376
for (z=0;z<g_oaMolecules.GetSize();z++)
3377
{
3378
m = (CMolecule*)g_oaMolecules[z];
3379
3380
for (z2=0;z2<m->m_laSingleMolIndex.GetSize();z2++)
3381
{
3382
sm = (CSingleMolecule*)g_oaSingleMolecules[m->m_laSingleMolIndex[z2]];
3383
3384
for (z3=0;z3<m->m_baAtomIndex.GetSize();z3++)
3385
{
3386
if (m->m_baAtomIndex[z3] == g_iVirtAtomType)
3387
continue;
3388
3389
el = ((CAtom*)g_oaAtoms[m->m_baAtomIndex[z3]])->m_pElement;
3390
3391
for (z4=0;z4<((CxIntArray*)sm->m_oaAtomOffset[z3])->GetSize();z4++)
3392
{
3393
o = ((CxIntArray*)sm->m_oaAtomOffset[z3])->GetAt(z4);
3394
vec1 = ts->m_vaCoords[o];
3395
vec1 /= 1000.0;
3396
3397
// mfprintf(b,"#if (atom_draw)\n");
3398
mfprintf(b," sphere { <%g, %g, %g>, elem_%s_r*atom_r\n",vec1[0],vec1[1],vec1[2],el->m_sLabel);
3399
3400
if (m_iaPOVSMCluster[m->m_laSingleMolIndex[z2]] != -1)
3401
{
3402
mfprintf(b," pigment { rgbft < %f, %f, %f, 0, 0 > } ",((m_iaPOVSMColor[m->m_laSingleMolIndex[z2]]&0x00FF0000)>>16)/255.0,((m_iaPOVSMColor[m->m_laSingleMolIndex[z2]]&0x0000FF00)>>8)/255.0,(m_iaPOVSMColor[m->m_laSingleMolIndex[z2]]&0x000000FF)/255.0);
3403
} else mfprintf(b," pigment { rgbft m_atom_color(elem_%s_col) } ",el->m_sLabel);
3404
3405
mfprintf(b,"finish { reflection atom_reflect specular atom_specular ambient atom_ambient diffuse atom_diffuse } }\n");
3406
// mfprintf(b,"#end\n");
3407
3408
vec2 = cam - vec1;
3409
vec2.Normalize();
3410
3411
// if (shadow)
3412
{
3413
// mfprintf(b,"#if (atom_draw_halo1)\n");
3414
mfprintf(b," disc { < %g - (%g) * atom_d_halo1, %g - (%g) * atom_d_halo1, %g - (%g) * atom_d_halo1 >,\n",vec1[0],vec2[0],vec1[1],vec2[1],vec1[2],vec2[2]);
3415
// mfprintf(b," < %g, %g, %g >, (elem_%s_r * atom_r) + atom_r_halo1, elem_%s_r * atom_r\n",vec2[0],vec2[1],vec2[2],el->m_sLabel,el->m_sLabel);
3416
mfprintf(b," < %g, %g, %g >, (elem_%s_r * atom_r) + atom_r_halo1\n",vec2[0],vec2[1],vec2[2],el->m_sLabel);
3417
mfprintf(b," pigment { rgbft atom_color_halo1 } finish { reflection atom_reflect_halo1 specular atom_specular_halo1 ambient atom_ambient_halo1 diffuse atom_diffuse_halo1 } no_reflection no_radiosity }\n");
3418
// mfprintf(b,"#end\n");
3419
}
3420
}
3421
}
3422
3423
for (z3=0;z3<sm->m_oaBonds.GetSize();z3++)
3424
{
3425
mb = (CMolBond*)sm->m_oaBonds[z3];
3426
o = mb->m_iAtomOffset[0];
3427
o2 = mb->m_iAtomOffset[1];
3428
3429
el = ((CAtom*)g_oaAtoms[g_waAtomRealElement[o]])->m_pElement;
3430
el2 = ((CAtom*)g_oaAtoms[g_waAtomRealElement[o2]])->m_pElement;
3431
3432
vec1 = ts->m_vaCoords[o];
3433
vec1 /= 1000.0;
3434
3435
vec2 = ts->m_vaCoords[o2];
3436
vec2 /= 1000.0;
3437
3438
if ((vec1-vec2).GetLength() > 0.3)
3439
continue;
3440
3441
vec3 = (vec1/el->m_fRadius + vec2/el2->m_fRadius) / (1.0/el->m_fRadius+1.0/el2->m_fRadius);
3442
3443
vec3b = vec2 - vec1;
3444
vec3b.Normalize();
3445
3446
mfprintf(b,"#local vec1c = < %g + (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_r*bond_r ) ) ),\n",vec1[0],vec3b[0],el->m_sLabel,el->m_sLabel,el->m_sLabel,el->m_sLabel);
3447
mfprintf(b," %g + (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_r*bond_r ) ) ),\n",vec1[1],vec3b[1],el->m_sLabel,el->m_sLabel,el->m_sLabel,el->m_sLabel);
3448
mfprintf(b," %g + (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_r*bond_r ) ) ) >;\n",vec1[2],vec3b[2],el->m_sLabel,el->m_sLabel,el->m_sLabel,el->m_sLabel);
3449
mfprintf(b,"#local vec2c = < %g - (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_r*bond_r ) ) ),\n",vec2[0],vec3b[0],el2->m_sLabel,el2->m_sLabel,el2->m_sLabel,el2->m_sLabel);
3450
mfprintf(b," %g - (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_r*bond_r ) ) ),\n",vec2[1],vec3b[1],el2->m_sLabel,el2->m_sLabel,el2->m_sLabel,el2->m_sLabel);
3451
mfprintf(b," %g - (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_r*bond_r ) ) ) >;\n",vec2[2],vec3b[2],el2->m_sLabel,el2->m_sLabel,el2->m_sLabel,el2->m_sLabel);
3452
3453
// mfprintf(b,"#if (bond_draw)\n");
3454
if (m_iaPOVSMCluster[m->m_laSingleMolIndex[z2]] != -1)
3455
{
3456
mfprintf(b," cylinder { vec1c, vec2c, bond_r open\n");
3457
mfprintf(b," pigment { rgbft < %f, %f, %f, 0, 0 > } finish { reflection bond_reflect specular bond_specular ambient bond_ambient diffuse bond_diffuse } }\n",((m_iaPOVSMColor[m->m_laSingleMolIndex[z2]]&0x00FF0000)>>16)/255.0,((m_iaPOVSMColor[m->m_laSingleMolIndex[z2]]&0x0000FF00)>>8)/255.0,(m_iaPOVSMColor[m->m_laSingleMolIndex[z2]]&0x000000FF)/255.0);
3458
} else
3459
{
3460
mfprintf(b," cylinder { < %g, %g, %g >, vec1c, bond_r open\n",vec3[0],vec3[1],vec3[2]);
3461
mfprintf(b," pigment { rgbft m_atom_color(elem_%s_col) } finish { reflection bond_reflect specular bond_specular ambient bond_ambient diffuse bond_diffuse } }\n",el->m_sLabel);
3462
mfprintf(b," cylinder { < %g, %g, %g >, vec2c, bond_r open\n",vec3[0],vec3[1],vec3[2]);
3463
mfprintf(b," pigment { rgbft m_atom_color(elem_%s_col) } finish { reflection bond_reflect specular bond_specular ambient bond_ambient diffuse bond_diffuse } }\n",el2->m_sLabel);
3464
}
3465
// mfprintf(b,"#end\n");
3466
3467
vec3 = cam - (vec1 + vec2) / 2.0;
3468
vec3.Normalize();
3469
3470
vec2b = vec2-vec1;
3471
vec1b = CrossP(vec3,vec2b);
3472
vec1b.Normalize();
3473
3474
mfprintf(b,"#local vec3 = < %g, %g, %g >;\n",vec3[0],vec3[1],vec3[2]);
3475
mfprintf(b,"#local vec1b = < %g, %g, %g >;\n",vec1b[0],vec1b[1],vec1b[2]);
3476
3477
// if (shadow)
3478
{
3479
// mfprintf(b,"#if (bond_draw_halo1)\n");
3480
mfprintf(b," #local vecA = vec1c + vec1b * (bond_r + bond_r_halo1) - vec3 * bond_d_halo1;\n");
3481
mfprintf(b," #local vecB = vec2c + vec1b * (bond_r + bond_r_halo1) - vec3 * bond_d_halo1;\n");
3482
mfprintf(b," #local vecC = vec2c + vec1b * (bond_r) - vec3 * bond_d_halo1;\n");
3483
mfprintf(b," #local vecD = vec1c + vec1b * (bond_r) - vec3 * bond_d_halo1;\n");
3484
mfprintf(b," triangle { vecA, vecB, vecC\n");
3485
mfprintf(b," pigment { rgbft bond_color_halo1 } finish { reflection bond_reflect_halo1 specular bond_specular_halo1 ambient bond_ambient_halo1 diffuse bond_diffuse_halo1 } no_reflection no_radiosity }\n");
3486
mfprintf(b," triangle { vecA, vecD, vecC\n");
3487
mfprintf(b," pigment { rgbft bond_color_halo1 } finish { reflection bond_reflect_halo1 specular bond_specular_halo1 ambient bond_ambient_halo1 diffuse bond_diffuse_halo1 } no_reflection no_radiosity }\n");
3488
mfprintf(b," #local vecA = vec1c - vec1b * (bond_r + bond_r_halo1) - vec3 * bond_d_halo1;\n");
3489
mfprintf(b," #local vecB = vec2c - vec1b * (bond_r + bond_r_halo1) - vec3 * bond_d_halo1;\n");
3490
mfprintf(b," #local vecC = vec2c - vec1b * (bond_r) - vec3 * bond_d_halo1;\n");
3491
mfprintf(b," #local vecD = vec1c - vec1b * (bond_r) - vec3 * bond_d_halo1;\n");
3492
mfprintf(b," triangle { vecA, vecB, vecC\n");
3493
mfprintf(b," pigment { rgbft bond_color_halo1 } finish { reflection bond_reflect_halo1 specular bond_specular_halo1 ambient bond_ambient_halo1 diffuse bond_diffuse_halo1 } no_reflection no_radiosity }\n");
3494
mfprintf(b," triangle { vecA, vecD, vecC\n");
3495
mfprintf(b," pigment { rgbft bond_color_halo1 } finish { reflection bond_reflect_halo1 specular bond_specular_halo1 ambient bond_ambient_halo1 diffuse bond_diffuse_halo1 } no_reflection no_radiosity }\n");
3496
// mfprintf(b,"#end\n");
3497
}
3498
3499
// if (shadow)
3500
{
3501
vec3 = vec2 - vec1;
3502
vec3.Normalize();
3503
3504
// mfprintf(b,"#if (bond_draw_stub)\n");
3505
mfprintf(b," #local vec3 = < %g, %g, %g >;\n",vec3[0],vec3[1],vec3[2]);
3506
mfprintf(b," difference {\n");
3507
mfprintf(b," cylinder { vec1c, vec1c + vec3 * bond_l_stub, bond_r + bond_r_stub }\n");
3508
mfprintf(b," cylinder { vec1c - vec3 * (bond_l_stub+0.001), vec1c + vec3 * (bond_l_stub+0.001), bond_r }\n");
3509
mfprintf(b," pigment { rgbft bond_color_stub } finish { reflection bond_reflect_stub specular bond_specular_stub ambient bond_ambient_stub diffuse bond_diffuse_stub } }\n");
3510
mfprintf(b," difference {\n");
3511
mfprintf(b," cylinder { vec2c, vec2c - vec3 * bond_l_stub, bond_r + bond_r_stub }\n");
3512
mfprintf(b," cylinder { vec2c + vec3 * (bond_l_stub+0.001), vec2c - vec3 * (bond_l_stub+0.001), bond_r }\n");
3513
mfprintf(b," pigment { rgbft bond_color_stub } finish { reflection bond_reflect_stub specular bond_specular_stub ambient bond_ambient_stub diffuse bond_diffuse_stub } }\n");
3514
// mfprintf(b,"#end\n");
3515
}
3516
}
3517
3518
}
3519
}
3520
3521
for (z=0;z<m_oaPOVExtendedCluster.GetSize();z++)
3522
{
3523
ec = (CExtendedCluster*)m_oaPOVExtendedCluster[z];
3524
3525
for (z2=0;z2<ec->m_iaAtoms.GetSize();z2++)
3526
{
3527
for (z3=0;z3<((CxIntArray*)ec->m_oaBonds[z2])->GetSize();z3++)
3528
{
3529
if (g_laAtomSMIndex[ec->m_iaAtoms[z2]] == g_laAtomSMIndex[((CxIntArray*)ec->m_oaBonds[z2])->GetAt(z3)])
3530
continue;
3531
3532
if (ec->m_iaAtoms[z2] > ((CxIntArray*)ec->m_oaBonds[z2])->GetAt(z3))
3533
continue;
3534
3535
o = ec->m_iaAtoms[z2];
3536
o2 = ((CxIntArray*)ec->m_oaBonds[z2])->GetAt(z3);
3537
3538
el = ((CAtom*)g_oaAtoms[g_waAtomRealElement[o]])->m_pElement;
3539
el2 = ((CAtom*)g_oaAtoms[g_waAtomRealElement[o2]])->m_pElement;
3540
3541
vec1 = ts->m_vaCoords[o];
3542
vec1 /= 1000.0;
3543
3544
vec2 = ts->m_vaCoords[o2];
3545
vec2 /= 1000.0;
3546
3547
// if ((vec1-vec2).GetLength() > 0.3)
3548
// continue;
3549
3550
vec3 = (vec1/el->m_fRadius + vec2/el2->m_fRadius) / (1.0/el->m_fRadius+1.0/el2->m_fRadius);
3551
3552
vec3b = vec2 - vec1;
3553
vec3b.Normalize();
3554
3555
mfprintf(b,"#local vec1c = < %g + (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_cluster_r*bond_cluster_r ) ) ),\n",vec1[0],vec3b[0],el->m_sLabel,el->m_sLabel,el->m_sLabel,el->m_sLabel);
3556
mfprintf(b," %g + (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_cluster_r*bond_cluster_r ) ) ),\n",vec1[1],vec3b[1],el->m_sLabel,el->m_sLabel,el->m_sLabel,el->m_sLabel);
3557
mfprintf(b," %g + (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_cluster_r*bond_cluster_r ) ) ) >;\n",vec1[2],vec3b[2],el->m_sLabel,el->m_sLabel,el->m_sLabel,el->m_sLabel);
3558
mfprintf(b,"#local vec2c = < %g - (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_cluster_r*bond_cluster_r ) ) ),\n",vec2[0],vec3b[0],el2->m_sLabel,el2->m_sLabel,el2->m_sLabel,el2->m_sLabel);
3559
mfprintf(b," %g - (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_cluster_r*bond_cluster_r ) ) ),\n",vec2[1],vec3b[1],el2->m_sLabel,el2->m_sLabel,el2->m_sLabel,el2->m_sLabel);
3560
mfprintf(b," %g - (%g) * ( elem_%s_r*atom_r - ( elem_%s_r*atom_r - 0.5*sqrt(4*elem_%s_r*elem_%s_r*atom_r*atom_r - 4*bond_cluster_r*bond_cluster_r ) ) ) >;\n",vec2[2],vec3b[2],el2->m_sLabel,el2->m_sLabel,el2->m_sLabel,el2->m_sLabel);
3561
3562
// mfprintf(b,"#if (bond_draw)\n");
3563
mfprintf(b," cylinder { vec1c, vec2c, bond_cluster_r open\n");
3564
mfprintf(b," pigment { rgbft < %f, %f, %f, 0, 0 > } finish { reflection bond_reflect specular bond_specular ambient bond_ambient diffuse bond_diffuse } }\n",((m_iaPOVSMColor[g_laAtomSMIndex[ec->m_iaAtoms[z2]]]&0x00FF0000)>>16)/255.0,((m_iaPOVSMColor[g_laAtomSMIndex[ec->m_iaAtoms[z2]]]&0x0000FF00)>>8)/255.0,(m_iaPOVSMColor[g_laAtomSMIndex[ec->m_iaAtoms[z2]]]&0x000000FF)/255.0);
3565
// mfprintf(b,"#end\n");
3566
3567
vec3 = cam - (vec1 + vec2) / 2.0;
3568
vec3.Normalize();
3569
3570
vec2b = vec2-vec1;
3571
vec1b = CrossP(vec3,vec2b);
3572
vec1b.Normalize();
3573
3574
mfprintf(b,"#local vec3 = < %g, %g, %g >;\n",vec3[0],vec3[1],vec3[2]);
3575
mfprintf(b,"#local vec1b = < %g, %g, %g >;\n",vec1b[0],vec1b[1],vec1b[2]);
3576
3577
// if (shadow)
3578
{
3579
// mfprintf(b,"#if (bond_draw_halo1)\n");
3580
mfprintf(b," #local vecA = vec1c + vec1b * (bond_cluster_r + bond_r_halo1) - vec3 * bond_d_halo1;\n");
3581
mfprintf(b," #local vecB = vec2c + vec1b * (bond_cluster_r + bond_r_halo1) - vec3 * bond_d_halo1;\n");
3582
mfprintf(b," #local vecC = vec2c + vec1b * (bond_cluster_r) - vec3 * bond_d_halo1;\n");
3583
mfprintf(b," #local vecD = vec1c + vec1b * (bond_cluster_r) - vec3 * bond_d_halo1;\n");
3584
mfprintf(b," triangle { vecA, vecB, vecC\n");
3585
mfprintf(b," pigment { rgbft bond_color_halo1 } finish { reflection bond_reflect_halo1 specular bond_specular_halo1 ambient bond_ambient_halo1 diffuse bond_diffuse_halo1 } no_reflection no_radiosity }\n");
3586
mfprintf(b," triangle { vecA, vecD, vecC\n");
3587
mfprintf(b," pigment { rgbft bond_color_halo1 } finish { reflection bond_reflect_halo1 specular bond_specular_halo1 ambient bond_ambient_halo1 diffuse bond_diffuse_halo1 } no_reflection no_radiosity }\n");
3588
mfprintf(b," #local vecA = vec1c - vec1b * (bond_cluster_r + bond_r_halo1) - vec3 * bond_d_halo1;\n");
3589
mfprintf(b," #local vecB = vec2c - vec1b * (bond_cluster_r + bond_r_halo1) - vec3 * bond_d_halo1;\n");
3590
mfprintf(b," #local vecC = vec2c - vec1b * (bond_cluster_r) - vec3 * bond_d_halo1;\n");
3591
mfprintf(b," #local vecD = vec1c - vec1b * (bond_cluster_r) - vec3 * bond_d_halo1;\n");
3592
mfprintf(b," triangle { vecA, vecB, vecC\n");
3593
mfprintf(b," pigment { rgbft bond_color_halo1 } finish { reflection bond_reflect_halo1 specular bond_specular_halo1 ambient bond_ambient_halo1 diffuse bond_diffuse_halo1 } no_reflection no_radiosity }\n");
3594
mfprintf(b," triangle { vecA, vecD, vecC\n");
3595
mfprintf(b," pigment { rgbft bond_color_halo1 } finish { reflection bond_reflect_halo1 specular bond_specular_halo1 ambient bond_ambient_halo1 diffuse bond_diffuse_halo1 } no_reflection no_radiosity }\n");
3596
// mfprintf(b,"#end\n");
3597
}
3598
3599
// if (shadow)
3600
{
3601
vec3 = vec2 - vec1;
3602
vec3.Normalize();
3603
3604
// mfprintf(b,"#if (bond_draw_stub)\n");
3605
mfprintf(b," #local vec3 = < %g, %g, %g >;\n",vec3[0],vec3[1],vec3[2]);
3606
mfprintf(b," difference {\n");
3607
mfprintf(b," cylinder { vec1c, vec1c + vec3 * bond_l_stub, bond_cluster_r + bond_r_stub }\n");
3608
mfprintf(b," cylinder { vec1c - vec3 * (bond_l_stub+0.001), vec1c + vec3 * (bond_l_stub+0.001), bond_cluster_r }\n");
3609
mfprintf(b," pigment { rgbft bond_color_stub } finish { reflection bond_reflect_stub specular bond_specular_stub ambient bond_ambient_stub diffuse bond_diffuse_stub } }\n");
3610
mfprintf(b," difference {\n");
3611
mfprintf(b," cylinder { vec2c, vec2c - vec3 * bond_l_stub, bond_cluster_r + bond_r_stub }\n");
3612
mfprintf(b," cylinder { vec2c + vec3 * (bond_l_stub+0.001), vec2c - vec3 * (bond_l_stub+0.001), bond_cluster_r }\n");
3613
mfprintf(b," pigment { rgbft bond_color_stub } finish { reflection bond_reflect_stub specular bond_specular_stub ambient bond_ambient_stub diffuse bond_diffuse_stub } }\n");
3614
// mfprintf(b,"#end\n");
3615
}
3616
}
3617
}
3618
}
3619
3620
mfprintf(b,"\n no_shadow\n}\n\n");
3621
3622
fclose(b);
3623
}
3624
3625
3626
void CClusterAnalysis::WriteClusterXYZ(CClusterNode *n, CTimeStep *ts, FILE *a)
3627
{
3628
int z2;
3629
float avgx, avgy, avgz;
3630
3631
mfprintf(a,"%d\nStep %d\n",n->m_iaAtoms.GetSize(),g_iSteps);
3632
3633
avgx = 0;
3634
avgy = 0;
3635
avgz = 0;
3636
for (z2=0;z2<n->m_iaAtoms.GetSize();z2++)
3637
{
3638
if (z2 != 0)
3639
{
3640
// mprintf("\nOrig: (%f|%f|%f), Avg: (%f|%f|%f)",ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][0],ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][1],ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][2],avgx/z2, avgy/z2, avgz/z2);
3641
while (ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][0]-(avgx/z2) > g_fBoxX/2)
3642
ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][0] -= g_fBoxX;
3643
while (ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][0]-(avgx/z2) < -g_fBoxX/2)
3644
ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][0] += g_fBoxX;
3645
while (ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][1]-(avgy/z2) > g_fBoxY/2)
3646
ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][1] -= g_fBoxY;
3647
while (ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][1]-(avgy/z2) < -g_fBoxY/2)
3648
ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][1] += g_fBoxY;
3649
while (ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][2]-(avgz/z2) > g_fBoxZ/2)
3650
ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][2] -= g_fBoxZ;
3651
while (ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][2]-(avgz/z2) < -g_fBoxZ/2)
3652
ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][2] += g_fBoxZ;
3653
// mprintf("\n New: (%f|%f|%f)",ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][0],ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][1],ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][2]);
3654
}
3655
avgx += ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][0];
3656
avgy += ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][1];
3657
avgz += ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][2];
3658
}
3659
avgx /= n->m_iaAtoms.GetSize();
3660
avgy /= n->m_iaAtoms.GetSize();
3661
avgz /= n->m_iaAtoms.GetSize();
3662
3663
for (z2=0;z2<n->m_iaAtoms.GetSize();z2++)
3664
mfprintf(a,"%s %f %f %f\n",((CAtom*)g_oaAtoms[g_waAtomRealElement[m_iaAtomList[n->m_iaAtoms[z2]]]])->m_sName,(ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][0]-avgx)/100.0,(ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][1]-avgy)/100.0,(ts->m_vaCoords[m_iaAtomList[n->m_iaAtoms[z2]]][2]-avgz)/100.0);
3665
3666
fflush(a);
3667
}
3668
3669
3670
unsigned int CClusterAnalysis::POVColorFunction(int i)
3671
{
3672
switch(i % 10)
3673
{
3674
case 0: return 0x000000E0; // Dunkelblau
3675
case 1: return 0x00009000; // Dunkelgruen
3676
case 2: return 0x00C00000; // Dunkelrot
3677
case 3: return 0x00FFFF00; // Gelb
3678
case 4: return 0x00FF00FF; // Violett
3679
case 5: return 0x0000FFFF; // Cyan
3680
case 6: return 0x005050FF; // Hellblau
3681
case 7: return 0x0060FF60; // Hellgruen
3682
case 8: return 0x00FF3060; // Hellrot
3683
case 9: return 0x00E07000; // Braun
3684
}
3685
3686
return 0;
3687
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1