← Back to branch summary

~ubuntu-branches/ubuntu/trusty/travis/trusty-proposed

« back to all changes in this revision

Viewing changes to src/pdf.cpp

  • Committer: Package Import Robot
  • Author(s): Daniel Leidert
  • Date: 2014-01-18 20:07:16 UTC
  • mfrom: (1.1.8)
  • Revision ID: package-import@ubuntu.com-20140118200716-whsmcg7fa1eyqecq
Tags: 140117-1
New upstream release.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/pdf.cpp
1
 
/*****************************************************************************
2
 
    TRAVIS - Trajectory Analyzer and Visualizer
3
 
    http://www.travis-analyzer.de/
4
 
 
5
 
    Copyright (c) 2009-2013 Martin Brehm
6
 
                  2012-2013 Martin Thomas
7
 
 
8
 
    This file written by Martin Thomas.
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
 
#include "pdf.h"
25
 
 
26
 
#include "atomgroup.h"
27
 
#include "df.h"
28
 
#include "globalvar.h"
29
 
#include "maintools.h"
30
 
#include "xobarray.h"
31
 
#include "xvector3.h"
32
 
 
33
 
#define BUF_SIZE 1024
34
 
 
35
 
static CxVector3 g_normalVector;
36
 
static CxVector3 g_fixPoint;
37
 
 
38
 
static CxObArray g_pdfObserv;
39
 
 
40
 
CPDF::CPDF(int showMol) {
41
 
        m_iShowMol = showMol;
42
 
        _showAtomGes = 0;
43
 
        
44
 
        try { _df = new CDF(); } catch(...) { _df = NULL; }
45
 
        if(_df == NULL) NewException((double)sizeof(CDF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
46
 
        
47
 
        mprintf(WHITE, "\n>>> Fixed Plane Density Profile >>>\n\n");
48
 
        
49
 
        try { _ag = new CAtomGroup(); } catch(...) { _ag = NULL; }
50
 
        if(_ag == NULL) NewException((double)sizeof(CAtomGroup), __FILE__, __LINE__, __PRETTY_FUNCTION__);
51
 
        
52
 
        while(true) {
53
 
                if(((CMolecule *)g_oaMolecules[m_iShowMol])->m_iAtomGes == 3) {
54
 
                        mprintf("    %s is only one atom, there is no choice.\n", ((CMolecule *)g_oaMolecules[m_iShowMol])->m_sName);
55
 
                        _ag->Reset();
56
 
                        _ag->m_pMolecule = (CMolecule *)g_oaMolecules[m_iShowMol];
57
 
                        _ag->AddAtom(0, 0, false);
58
 
                        _ag->SortAtoms();
59
 
                        _ag->BuildName();
60
 
                } else {
61
 
                        mprintf("    Which atom(s) to take from OM %s (e.g. \"C1,C3-5,H\", \"*\"=all)? [#2] ",((CMolecule*)g_oaMolecules[m_iShowMol])->m_sName);
62
 
                        inpprintf("! Which atom(s) to take from OM %s (e.g. \"C1,C3-5,H\", \"*\"=all)? [#2]\n",((CMolecule*)g_oaMolecules[m_iShowMol])->m_sName);
63
 
                        
64
 
                        char buf[256];
65
 
                        myget(buf);
66
 
                        if(strlen(buf) == 0) {
67
 
                                if(!_ag->ParseAtoms((CMolecule *)g_oaMolecules[m_iShowMol], "#2")) {
68
 
                                        eprintf("Weird error.\n");
69
 
                                        continue;
70
 
                                }
71
 
                        } else if(!_ag->ParseAtoms((CMolecule *)g_oaMolecules[m_iShowMol], buf)) {
72
 
                                continue;
73
 
                        }
74
 
                }
75
 
                break;
76
 
        }
77
 
        _showAtomGes += _ag->m_iAtomGes;
78
 
        
79
 
        ParseDeriv();
80
 
        switch(m_iDeriv) {
81
 
                case 0:
82
 
                        _minDist = AskFloat("    Enter the minimal distance of this Density Profile in pm: [%d.0] ", -(float)HalfBox(), -HalfBox());
83
 
                        _maxDist = AskFloat("    Enter the maximal distance of this Density Profile in pm: [%d.0] ", (float)HalfBox(), HalfBox());
84
 
                        break;
85
 
/*              case 1:
86
 
                        if(m_bDerivAbs)
87
 
                                _minDist = AskFloat("    Enter the minimal value of this d1-PDF in pm/ps: [0] ", 0.0f);
88
 
                        else
89
 
                                _minDist = AskFloat("    Enter the minimal value of this d1-PDF in pm/ps: [-10.0] ", -10.0f);
90
 
                        _maxDist = AskFloat("    Enter the maximal value of this d1-RDF in pm/ps: [10.0] ", 10.0f);
91
 
                        break;
92
 
                case 2:
93
 
                        if(m_bDerivAbs)
94
 
                                _minDist = AskFloat("    Enter the minimal value of this d2-RDF in pm/ps^2: [0] ", 0.0f);
95
 
                        else
96
 
                                _minDist = AskFloat("    Enter the minimal value of this d2-RDF in pm/ps^2: [-10.0] ", -10.0f);
97
 
                        _maxDist = AskFloat("    Enter the maximal value of this d2-RDF in pm/ps^2: [10.0] ",10.0f);
98
 
                        break;*/
99
 
                default:
100
 
                        eprintf("Higher derivatives are not implemented.\n");
101
 
                        abort();
102
 
                        break;
103
 
        }
104
 
        m_iResolution = AskUnsignedInteger("    Enter the resolution (bin count) for this Density Profile: [300] ", 300);
105
 
        m_iHistogramRes = 0;
106
 
        
107
 
        buildName();
108
 
        m_bSelf = false;
109
 
        
110
 
        mprintf(WHITE, "\n<<< Enf of Fixed Plane Density Profile <<<\n\n");
111
 
}
112
 
 
113
 
CPDF::~CPDF() {
114
 
        delete _ag;
115
 
        delete _df;
116
 
        delete[] m_faData;
117
 
        delete[] m_sName;
118
 
        delete[] m_sShortName;
119
 
}
120
 
 
121
 
void CPDF::initialize(int showMolCount) {
122
 
        mprintf("    Creating Density Profile...\n");
123
 
        _df->m_fMinVal = _minDist;
124
 
        _df->m_fMaxVal = _maxDist;
125
 
        _df->m_iResolution = m_iResolution;
126
 
        _df->m_iHistogramRes = m_iHistogramRes;
127
 
        _df->SetLabelX("Distance / pm");
128
 
        _df->SetLabelY("g(d)");
129
 
        _df->Create();
130
 
        
131
 
        try { m_faData = new CxDoubleArray[showMolCount]; } catch(...) { m_faData = NULL; }
132
 
        if(m_faData == NULL) NewException((double)showMolCount*sizeof(CxDoubleArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
133
 
        
134
 
        // Add timedev here
135
 
}
136
 
 
137
 
void CPDF::finalize() {
138
 
        mprintf("    %.0f bin entries, %.0f out of bin range (%.2f percent).\n", _df->m_fBinEntries, _df->m_fSkipEntries, ZeroDivide(_df->m_fSkipEntries, _df->m_fBinEntries + _df->m_fSkipEntries) * 100.0f);
139
 
        _df->CalcMeanSD();
140
 
        mprintf("    Mean value: %10G pm    Standard deviation: %10G pm\n", _df->m_fMean, _df->m_fSD);
141
 
        mprintf("    Min. value: %10G pm    Max.value:          %10G pm\n", _df->m_fMinInput, _df->m_fMaxInput);
142
 
        
143
 
        mprintf("    Scaling Density Profile to uniform density...\n");
144
 
        _df->MultiplyBin((g_normalVector[0] * g_normalVector[0] * g_fBoxX + g_normalVector[1] * g_normalVector[1] * g_fBoxY + g_normalVector[2] * g_normalVector[2] * g_fBoxZ) * m_iResolution / (_maxDist-_minDist) / g_iSteps / ((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex.GetSize() / _showAtomGes);
145
 
        
146
 
        char filename[BUF_SIZE];
147
 
 
148
 
#ifdef TARGET_WINDOWS
149
 
        _snprintf(filename, BUF_SIZE, "dprof_%s.csv", m_sName);
150
 
#elif defined TARGET_LINUX
151
 
        snprintf(filename, BUF_SIZE, "dprof_%s.csv", m_sName);
152
 
#else
153
 
        sprintf(filename, "dprof_%s.csv", m_sName);
154
 
#endif
155
 
 
156
 
        mprintf("    Saving density profile as \"%s\"...\n", filename);
157
 
        _df->Write("", filename, "", true);
158
 
 
159
 
#ifdef TARGET_WINDOWS
160
 
        _snprintf(filename, BUF_SIZE, "dprof_%s.agr", m_sName);
161
 
#elif defined TARGET_LINUX
162
 
        snprintf(filename, BUF_SIZE, "dprof_%s.agr", m_sName);
163
 
#else
164
 
        sprintf(filename, "dprof_%s.agr", m_sName);
165
 
#endif
166
 
 
167
 
        mprintf("    Saving density profile AGR as \"%s\"...\n", filename);
168
 
        _df->WriteAgr("", filename, "", m_sName, true);
169
 
}
170
 
 
171
 
void CPDF::buildAtomList(CSingleMolecule *sm, CxIntArray *array) {
172
 
        int i, j;
173
 
        
174
 
        array->RemoveAll_KeepSize();
175
 
        
176
 
        for(i = 0; i < _ag->m_baAtomType.GetSize(); i++) {
177
 
                CxIntArray *a = (CxIntArray *)_ag->m_oaAtoms[i];
178
 
                for(j = 0; j < a->GetSize(); j++) {
179
 
                        array->Add(((CxIntArray *)sm->m_oaAtomOffset[_ag->m_baAtomType[i]])->GetAt(a->GetAt(j)));
180
 
                }
181
 
        }
182
 
}
183
 
 
184
 
void CPDF::addToDF(float value) {
185
 
        _df->AddToBin(value);
186
 
}
187
 
 
188
 
void CPDF::buildName() {
189
 
        char tmp[BUF_SIZE];
190
 
        
191
 
        // Check for overflow is missing!!!
192
 
        tmp[0] = 0;
193
 
        if(m_iDeriv != 0) {
194
 
                sprintf(tmp, "deriv%d_", m_iDeriv);
195
 
        }
196
 
        strcat(tmp, "[");
197
 
        strcat(tmp, _ag->m_sName);
198
 
        strcat(tmp, "]");
199
 
        
200
 
        try { m_sShortName = new char[1]; } catch(...) { m_sShortName = NULL; }
201
 
        if(m_sShortName == NULL) NewException((double)sizeof(char), __FILE__, __LINE__, __PRETTY_FUNCTION__);
202
 
        
203
 
        m_sShortName[0] = 0;
204
 
        
205
 
        try { m_sName = new char[strlen(tmp)+1]; } catch(...) { m_sName = NULL; }
206
 
        if(m_sName == NULL) NewException((double)(strlen(tmp)+1)*sizeof(char), __FILE__, __LINE__, __PRETTY_FUNCTION__);
207
 
        
208
 
        strcpy(m_sName, tmp);
209
 
}
210
 
 
211
 
CPDFObservation::CPDFObservation() {
212
 
        int i;
213
 
        
214
 
        m_pConditions = NULL;
215
 
        m_bTimeDev = false;
216
 
        m_bSelf = false;
217
 
        m_bOthers = true;
218
 
        
219
 
        if(g_oaMolecules.GetSize() > 1) {
220
 
                char buf[BUF_SIZE], buf2[BUF_SIZE];
221
 
                size_t remaining = BUF_SIZE;
222
 
                
223
 
#ifdef TARGET_WINDOWS
224
 
                remaining -= _snprintf(buf, remaining, "    Which molecule should be observed (");
225
 
#elif defined TARGET_LINUX
226
 
                remaining -= snprintf(buf, remaining, "    Which molecule should be observed (");
227
 
#else
228
 
                remaining -= sprintf(buf, "    Which molecule should be observed (");
229
 
#endif
230
 
                
231
 
                for(i = 0; i < g_oaMolecules.GetSize(); i++) {
232
 
                        if(remaining < 1)
233
 
                                break;
234
 
 
235
 
#ifdef TARGET_WINDOWS
236
 
                        size_t length = _snprintf(buf2, remaining, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
237
 
#elif defined TARGET_LINUX
238
 
                        size_t length = snprintf(buf2, remaining, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
239
 
#else
240
 
                        size_t length = sprintf(buf2, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
241
 
#endif
242
 
 
243
 
                        strncat(buf, buf2, remaining - 1);
244
 
                        remaining -= length;
245
 
                        if(i < g_oaMolecules.GetSize() - 1) {
246
 
 
247
 
#ifdef TARGET_WINDOWS
248
 
                                length = _snprintf(buf2, remaining, ", ");
249
 
#elif defined TARGET_LINUX
250
 
                                length = snprintf(buf2, remaining, ", ");
251
 
#else
252
 
                                length = sprintf(buf2, ", ");
253
 
#endif
254
 
 
255
 
                                strncat(buf, buf2, remaining - 1);
256
 
                                remaining -= length;
257
 
                        }
258
 
                }
259
 
                strncat(buf, ")? ", remaining - 1);
260
 
                m_iShowMol = AskRangeInteger_ND(buf, 1, g_oaMolecules.GetSize()) - 1;
261
 
        } else {
262
 
                m_iShowMol = 0;
263
 
        }
264
 
        m_iShowMolCount = ((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex.GetSize();
265
 
        m_bObsCertain = false;
266
 
        m_bDecompDist = false;
267
 
        m_bDecompType = false;
268
 
        
269
 
        try { _pdf = new CPDF(m_iShowMol); } catch(...) { _pdf = NULL; }
270
 
        if(_pdf == NULL) NewException((double)sizeof(CPDF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
271
 
        
272
 
        // Ask for conditions here
273
 
        
274
 
        // Add observation list here
275
 
}
276
 
 
277
 
CPDFObservation::~CPDFObservation() {
278
 
        delete _pdf;
279
 
}
280
 
 
281
 
void CPDFObservation::initialize() {
282
 
        _pdf->initialize(m_iShowMolCount);
283
 
}
284
 
 
285
 
void CPDFObservation::process(CTimeStep *ts) {
286
 
        int i, j;
287
 
        // Process conditions here
288
 
        
289
 
        for(i = 0; i < m_iShowMolCount; i++) {
290
 
                _pdf->m_faData[i].RemoveAll_KeepSize();
291
 
        }
292
 
        
293
 
        for(i = 0; i < m_iShowMolCount; i++) {
294
 
                CSingleMolecule *sm = (CSingleMolecule *)g_oaSingleMolecules[((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex[i]];
295
 
                CxIntArray temp;
296
 
                _pdf->buildAtomList(sm, &temp);
297
 
                for(j = 0; j < temp.GetSize(); j++) {
298
 
                        CxVector3 vec = FoldVector(ts->m_vaCoords[temp[j]] - g_fixPoint);
299
 
                        float val = DotP(vec, g_normalVector);
300
 
                        _pdf->m_faData[i].Add(val);
301
 
                }
302
 
        }
303
 
        
304
 
        for(i = 0; i < m_iShowMolCount; i++) {
305
 
                for(j = 0; j < _pdf->m_faData[i].GetSize(); j++) {
306
 
                        _pdf->addToDF(_pdf->m_faData[i][j]);
307
 
                }
308
 
        }
309
 
}
310
 
 
311
 
void CPDFObservation::finalize() {
312
 
        mprintf(WHITE, "* Planar Distribution Function\n");
313
 
        _pdf->finalize();
314
 
}
315
 
 
316
 
bool gatherPDF() {
317
 
        mprintf(YELLOW, "\n>>> Reference Plane >>>\n\n");
318
 
        
319
 
        mprintf("    Enter normal vector of reference plane:\n");
320
 
        float x, y, z;
321
 
        x = AskFloat("    x component [0.0] ", 0.0f);
322
 
        y = AskFloat("    y component [0.0] ", 0.0f);
323
 
        z = AskFloat("    z component [1.0] ", 1.0f);
324
 
        g_normalVector = CxVector3(x, y, z);
325
 
        g_normalVector.Normalize();
326
 
        mprintf("\n    Normalized normal vector is\n");
327
 
        mprintf("        ");
328
 
        g_normalVector.Dump();
329
 
        mprintf("\n");
330
 
        
331
 
        mprintf("\n    Enter coordinates of fix point in pm:\n");
332
 
        x = AskFloat("    x component [0.0] ", 0.0f);
333
 
        y = AskFloat("    y component [0.0] ", 0.0f);
334
 
        z = AskFloat("    z component [0.0] ", 0.0f);
335
 
        g_fixPoint = CxVector3(x, y, z);
336
 
        mprintf("\n    Fix point vector is\n");
337
 
        mprintf("        ");
338
 
        g_fixPoint.Dump();
339
 
        mprintf("\n");
340
 
        
341
 
        mprintf(YELLOW, "\n<<< End of Reference Plane <<<\n\n");
342
 
        
343
 
        while(true) {
344
 
                mprintf(YELLOW, "\n>>> Density Profile Observation %d >>>\n\n", g_pdfObserv.GetSize() + 1);
345
 
                
346
 
                CPDFObservation *obs;
347
 
                try { obs = new CPDFObservation(); } catch(...) { obs = NULL; }
348
 
                if(obs == NULL) NewException((double)sizeof(CPDFObservation), __FILE__, __LINE__, __PRETTY_FUNCTION__);
349
 
                g_pdfObserv.Add(obs);
350
 
                
351
 
                mprintf(YELLOW, "\n<<< End of Density Profile Observation %d <<<\n\n", g_pdfObserv.GetSize());
352
 
                        
353
 
                if(!AskYesNo("    Add another observation (y/n)? [no] ", false))
354
 
                        break;
355
 
                mprintf("\n");
356
 
        }
357
 
        
358
 
        return true;
359
 
}
360
 
 
361
 
bool initializePDF() {
362
 
        int i;
363
 
        for(i = 0; i < g_pdfObserv.GetSize(); i++) {
364
 
                ((CPDFObservation *)g_pdfObserv[i])->initialize();
365
 
        }
366
 
        return true;
367
 
}
368
 
 
369
 
void processPDF(CTimeStep *ts) {
370
 
        int i;
371
 
        for(i = 0; i < g_pdfObserv.GetSize(); i++) {
372
 
                ((CPDFObservation *)g_pdfObserv[i])->process(ts);
373
 
        }
374
 
}
375
 
 
376
 
void finalizePDF() {
377
 
        int i;
378
 
        for(i = 0; i < g_pdfObserv.GetSize(); i++) {
379
 
                mprintf(YELLOW, "\n>>> Density Profile Observation %d >>>\n\n", i+1);
380
 
                ((CPDFObservation *)g_pdfObserv[i])->finalize();
381
 
                mprintf(YELLOW, "\n<<< End of Density Profile Observation %d <<<\n\n", i+1);
382
 
        }
383
 
        
384
 
        for(i = 0; i < g_pdfObserv.GetSize(); i++) {
385
 
                delete (CPDFObservation *)g_pdfObserv[i];
386
 
        }
387
 
}
 
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 Thomas.
 
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
#include "pdf.h"
 
25
 
 
26
#include "atomgroup.h"
 
27
#include "df.h"
 
28
#include "globalvar.h"
 
29
#include "maintools.h"
 
30
#include "xobarray.h"
 
31
#include "xvector3.h"
 
32
 
 
33
#define BUF_SIZE 1024
 
34
 
 
35
static CxVector3 g_normalVector;
 
36
static CxVector3 g_fixPoint;
 
37
 
 
38
static CxObArray g_pdfObserv;
 
39
 
 
40
CPDF::CPDF(int showMol) {
 
41
        m_iShowMol = showMol;
 
42
        _showAtomGes = 0;
 
43
        
 
44
        try { _df = new CDF(); } catch(...) { _df = NULL; }
 
45
        if(_df == NULL) NewException((double)sizeof(CDF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
46
        
 
47
        mprintf(WHITE, "\n>>> Fixed Plane Density Profile >>>\n\n");
 
48
        
 
49
        try { _ag = new CAtomGroup(); } catch(...) { _ag = NULL; }
 
50
        if(_ag == NULL) NewException((double)sizeof(CAtomGroup), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
51
        
 
52
        while(true) {
 
53
                if(((CMolecule *)g_oaMolecules[m_iShowMol])->m_iAtomGes == 3) {
 
54
                        mprintf("    %s is only one atom, there is no choice.\n", ((CMolecule *)g_oaMolecules[m_iShowMol])->m_sName);
 
55
                        _ag->Reset();
 
56
                        _ag->m_pMolecule = (CMolecule *)g_oaMolecules[m_iShowMol];
 
57
                        _ag->AddAtom(0, 0, false);
 
58
                        _ag->SortAtoms();
 
59
                        _ag->BuildName();
 
60
                } else {
 
61
                        mprintf("    Which atom(s) to take from OM %s (e.g. \"C1,C3-5,H\", \"*\"=all)? [#2] ",((CMolecule*)g_oaMolecules[m_iShowMol])->m_sName);
 
62
                        inpprintf("! Which atom(s) to take from OM %s (e.g. \"C1,C3-5,H\", \"*\"=all)? [#2]\n",((CMolecule*)g_oaMolecules[m_iShowMol])->m_sName);
 
63
                        
 
64
                        char buf[256];
 
65
                        myget(buf);
 
66
                        if(strlen(buf) == 0) {
 
67
                                if(!_ag->ParseAtoms((CMolecule *)g_oaMolecules[m_iShowMol], "#2")) {
 
68
                                        eprintf("Weird error.\n");
 
69
                                        continue;
 
70
                                }
 
71
                        } else if(!_ag->ParseAtoms((CMolecule *)g_oaMolecules[m_iShowMol], buf)) {
 
72
                                continue;
 
73
                        }
 
74
                }
 
75
                break;
 
76
        }
 
77
        _showAtomGes += _ag->m_iAtomGes;
 
78
        
 
79
        ParseDeriv();
 
80
        switch(m_iDeriv) {
 
81
                case 0:
 
82
                        _minDist = AskFloat("    Enter the minimal distance of this Density Profile in pm: [%d.0] ", -(float)HalfBox(), -HalfBox());
 
83
                        _maxDist = AskFloat("    Enter the maximal distance of this Density Profile in pm: [%d.0] ", (float)HalfBox(), HalfBox());
 
84
                        break;
 
85
/*              case 1:
 
86
                        if(m_bDerivAbs)
 
87
                                _minDist = AskFloat("    Enter the minimal value of this d1-PDF in pm/ps: [0] ", 0.0f);
 
88
                        else
 
89
                                _minDist = AskFloat("    Enter the minimal value of this d1-PDF in pm/ps: [-10.0] ", -10.0f);
 
90
                        _maxDist = AskFloat("    Enter the maximal value of this d1-RDF in pm/ps: [10.0] ", 10.0f);
 
91
                        break;
 
92
                case 2:
 
93
                        if(m_bDerivAbs)
 
94
                                _minDist = AskFloat("    Enter the minimal value of this d2-RDF in pm/ps^2: [0] ", 0.0f);
 
95
                        else
 
96
                                _minDist = AskFloat("    Enter the minimal value of this d2-RDF in pm/ps^2: [-10.0] ", -10.0f);
 
97
                        _maxDist = AskFloat("    Enter the maximal value of this d2-RDF in pm/ps^2: [10.0] ",10.0f);
 
98
                        break;*/
 
99
                default:
 
100
                        eprintf("Higher derivatives are not implemented.\n");
 
101
                        abort();
 
102
                        break;
 
103
        }
 
104
        m_iResolution = AskUnsignedInteger("    Enter the resolution (bin count) for this Density Profile: [300] ", 300);
 
105
        m_iHistogramRes = 0;
 
106
        
 
107
        buildName();
 
108
        m_bSelf = false;
 
109
        
 
110
        mprintf(WHITE, "\n<<< Enf of Fixed Plane Density Profile <<<\n\n");
 
111
}
 
112
 
 
113
CPDF::~CPDF() {
 
114
        delete _ag;
 
115
        delete _df;
 
116
        delete[] m_faData;
 
117
        delete[] m_sName;
 
118
        delete[] m_sShortName;
 
119
}
 
120
 
 
121
void CPDF::initialize(int showMolCount) {
 
122
        mprintf("    Creating Density Profile...\n");
 
123
        _df->m_fMinVal = _minDist;
 
124
        _df->m_fMaxVal = _maxDist;
 
125
        _df->m_iResolution = m_iResolution;
 
126
        _df->m_iHistogramRes = m_iHistogramRes;
 
127
        _df->SetLabelX("Distance / pm");
 
128
        _df->SetLabelY("g(d)");
 
129
        _df->Create();
 
130
        
 
131
        try { m_faData = new CxDoubleArray[showMolCount]; } catch(...) { m_faData = NULL; }
 
132
        if(m_faData == NULL) NewException((double)showMolCount*sizeof(CxDoubleArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
133
        
 
134
        // Add timedev here
 
135
}
 
136
 
 
137
void CPDF::finalize() {
 
138
        mprintf("    %.0f bin entries, %.0f out of bin range (%.2f percent).\n", _df->m_fBinEntries, _df->m_fSkipEntries, ZeroDivide(_df->m_fSkipEntries, _df->m_fBinEntries + _df->m_fSkipEntries) * 100.0f);
 
139
        _df->CalcMeanSD();
 
140
        mprintf("    Mean value: %10G pm    Standard deviation: %10G pm\n", _df->m_fMean, _df->m_fSD);
 
141
        mprintf("    Min. value: %10G pm    Max.value:          %10G pm\n", _df->m_fMinInput, _df->m_fMaxInput);
 
142
        
 
143
        mprintf("    Scaling Density Profile to uniform density...\n");
 
144
        _df->MultiplyBin((g_normalVector[0] * g_normalVector[0] * g_fBoxX + g_normalVector[1] * g_normalVector[1] * g_fBoxY + g_normalVector[2] * g_normalVector[2] * g_fBoxZ) * m_iResolution / (_maxDist-_minDist) / g_iSteps / ((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex.GetSize() / _showAtomGes);
 
145
        
 
146
        char filename[BUF_SIZE];
 
147
 
 
148
#ifdef TARGET_WINDOWS
 
149
        _snprintf(filename, BUF_SIZE, "dprof_%s.csv", m_sName);
 
150
#elif defined TARGET_LINUX
 
151
        snprintf(filename, BUF_SIZE, "dprof_%s.csv", m_sName);
 
152
#else
 
153
        sprintf(filename, "dprof_%s.csv", m_sName);
 
154
#endif
 
155
 
 
156
        mprintf("    Saving density profile as \"%s\"...\n", filename);
 
157
        _df->Write("", filename, "", true);
 
158
 
 
159
#ifdef TARGET_WINDOWS
 
160
        _snprintf(filename, BUF_SIZE, "dprof_%s.agr", m_sName);
 
161
#elif defined TARGET_LINUX
 
162
        snprintf(filename, BUF_SIZE, "dprof_%s.agr", m_sName);
 
163
#else
 
164
        sprintf(filename, "dprof_%s.agr", m_sName);
 
165
#endif
 
166
 
 
167
        mprintf("    Saving density profile AGR as \"%s\"...\n", filename);
 
168
        _df->WriteAgr("", filename, "", m_sName, true);
 
169
}
 
170
 
 
171
void CPDF::buildAtomList(CSingleMolecule *sm, CxIntArray *array) {
 
172
        int i, j;
 
173
        
 
174
        array->RemoveAll_KeepSize();
 
175
        
 
176
        for(i = 0; i < _ag->m_baAtomType.GetSize(); i++) {
 
177
                CxIntArray *a = (CxIntArray *)_ag->m_oaAtoms[i];
 
178
                for(j = 0; j < a->GetSize(); j++) {
 
179
                        array->Add(((CxIntArray *)sm->m_oaAtomOffset[_ag->m_baAtomType[i]])->GetAt(a->GetAt(j)));
 
180
                }
 
181
        }
 
182
}
 
183
 
 
184
void CPDF::addToDF(float value) {
 
185
        _df->AddToBin(value);
 
186
}
 
187
 
 
188
void CPDF::buildName() {
 
189
        char tmp[BUF_SIZE];
 
190
        
 
191
        // Check for overflow is missing!!!
 
192
        tmp[0] = 0;
 
193
        if(m_iDeriv != 0) {
 
194
                sprintf(tmp, "deriv%d_", m_iDeriv);
 
195
        }
 
196
        strcat(tmp, "[");
 
197
        strcat(tmp, _ag->m_sName);
 
198
        strcat(tmp, "]");
 
199
        
 
200
        try { m_sShortName = new char[1]; } catch(...) { m_sShortName = NULL; }
 
201
        if(m_sShortName == NULL) NewException((double)sizeof(char), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
202
        
 
203
        m_sShortName[0] = 0;
 
204
        
 
205
        try { m_sName = new char[strlen(tmp)+1]; } catch(...) { m_sName = NULL; }
 
206
        if(m_sName == NULL) NewException((double)(strlen(tmp)+1)*sizeof(char), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
207
        
 
208
        strcpy(m_sName, tmp);
 
209
}
 
210
 
 
211
CPDFObservation::CPDFObservation() {
 
212
        int i;
 
213
        
 
214
        m_pConditions = NULL;
 
215
        m_bTimeDev = false;
 
216
        m_bSelf = false;
 
217
        m_bOthers = true;
 
218
        
 
219
        if(g_oaMolecules.GetSize() > 1) {
 
220
                char buf[BUF_SIZE], buf2[BUF_SIZE];
 
221
                size_t remaining = BUF_SIZE;
 
222
                
 
223
#ifdef TARGET_WINDOWS
 
224
                remaining -= _snprintf(buf, remaining, "    Which molecule should be observed (");
 
225
#elif defined TARGET_LINUX
 
226
                remaining -= snprintf(buf, remaining, "    Which molecule should be observed (");
 
227
#else
 
228
                remaining -= sprintf(buf, "    Which molecule should be observed (");
 
229
#endif
 
230
                
 
231
                for(i = 0; i < g_oaMolecules.GetSize(); i++) {
 
232
                        if(remaining < 1)
 
233
                                break;
 
234
 
 
235
#ifdef TARGET_WINDOWS
 
236
                        size_t length = _snprintf(buf2, remaining, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
 
237
#elif defined TARGET_LINUX
 
238
                        size_t length = snprintf(buf2, remaining, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
 
239
#else
 
240
                        size_t length = sprintf(buf2, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
 
241
#endif
 
242
 
 
243
                        strncat(buf, buf2, remaining - 1);
 
244
                        remaining -= length;
 
245
                        if(i < g_oaMolecules.GetSize() - 1) {
 
246
 
 
247
#ifdef TARGET_WINDOWS
 
248
                                length = _snprintf(buf2, remaining, ", ");
 
249
#elif defined TARGET_LINUX
 
250
                                length = snprintf(buf2, remaining, ", ");
 
251
#else
 
252
                                length = sprintf(buf2, ", ");
 
253
#endif
 
254
 
 
255
                                strncat(buf, buf2, remaining - 1);
 
256
                                remaining -= length;
 
257
                        }
 
258
                }
 
259
                strncat(buf, ")? ", remaining - 1);
 
260
                m_iShowMol = AskRangeInteger_ND(buf, 1, g_oaMolecules.GetSize()) - 1;
 
261
        } else {
 
262
                m_iShowMol = 0;
 
263
        }
 
264
        m_iShowMolCount = ((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex.GetSize();
 
265
        m_bObsCertain = false;
 
266
        m_bDecompDist = false;
 
267
        m_bDecompType = false;
 
268
        
 
269
        try { _pdf = new CPDF(m_iShowMol); } catch(...) { _pdf = NULL; }
 
270
        if(_pdf == NULL) NewException((double)sizeof(CPDF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
271
        
 
272
        // Ask for conditions here
 
273
        
 
274
        // Add observation list here
 
275
}
 
276
 
 
277
CPDFObservation::~CPDFObservation() {
 
278
        delete _pdf;
 
279
}
 
280
 
 
281
void CPDFObservation::initialize() {
 
282
        _pdf->initialize(m_iShowMolCount);
 
283
}
 
284
 
 
285
void CPDFObservation::process(CTimeStep *ts) {
 
286
        int i, j;
 
287
        // Process conditions here
 
288
        
 
289
        for(i = 0; i < m_iShowMolCount; i++) {
 
290
                _pdf->m_faData[i].RemoveAll_KeepSize();
 
291
        }
 
292
        
 
293
        for(i = 0; i < m_iShowMolCount; i++) {
 
294
                CSingleMolecule *sm = (CSingleMolecule *)g_oaSingleMolecules[((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex[i]];
 
295
                CxIntArray temp;
 
296
                _pdf->buildAtomList(sm, &temp);
 
297
                for(j = 0; j < temp.GetSize(); j++) {
 
298
                        CxVector3 vec = FoldVector(ts->m_vaCoords[temp[j]] - g_fixPoint);
 
299
                        float val = DotP(vec, g_normalVector);
 
300
                        _pdf->m_faData[i].Add(val);
 
301
                }
 
302
        }
 
303
        
 
304
        for(i = 0; i < m_iShowMolCount; i++) {
 
305
                for(j = 0; j < _pdf->m_faData[i].GetSize(); j++) {
 
306
                        _pdf->addToDF(_pdf->m_faData[i][j]);
 
307
                }
 
308
        }
 
309
}
 
310
 
 
311
void CPDFObservation::finalize() {
 
312
        mprintf(WHITE, "* Planar Distribution Function\n");
 
313
        _pdf->finalize();
 
314
}
 
315
 
 
316
bool gatherPDF() {
 
317
        mprintf(YELLOW, "\n>>> Reference Plane >>>\n\n");
 
318
        
 
319
        mprintf("    Enter normal vector of reference plane:\n");
 
320
        float x, y, z;
 
321
        x = AskFloat("    x component [0.0] ", 0.0f);
 
322
        y = AskFloat("    y component [0.0] ", 0.0f);
 
323
        z = AskFloat("    z component [1.0] ", 1.0f);
 
324
        g_normalVector = CxVector3(x, y, z);
 
325
        g_normalVector.Normalize();
 
326
        mprintf("\n    Normalized normal vector is\n");
 
327
        mprintf("        ");
 
328
        g_normalVector.Dump();
 
329
        mprintf("\n");
 
330
        
 
331
        mprintf("\n    Enter coordinates of fix point in pm:\n");
 
332
        x = AskFloat("    x component [0.0] ", 0.0f);
 
333
        y = AskFloat("    y component [0.0] ", 0.0f);
 
334
        z = AskFloat("    z component [0.0] ", 0.0f);
 
335
        g_fixPoint = CxVector3(x, y, z);
 
336
        mprintf("\n    Fix point vector is\n");
 
337
        mprintf("        ");
 
338
        g_fixPoint.Dump();
 
339
        mprintf("\n");
 
340
        
 
341
        mprintf(YELLOW, "\n<<< End of Reference Plane <<<\n\n");
 
342
        
 
343
        while(true) {
 
344
                mprintf(YELLOW, "\n>>> Density Profile Observation %d >>>\n\n", g_pdfObserv.GetSize() + 1);
 
345
                
 
346
                CPDFObservation *obs;
 
347
                try { obs = new CPDFObservation(); } catch(...) { obs = NULL; }
 
348
                if(obs == NULL) NewException((double)sizeof(CPDFObservation), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
349
                g_pdfObserv.Add(obs);
 
350
                
 
351
                mprintf(YELLOW, "\n<<< End of Density Profile Observation %d <<<\n\n", g_pdfObserv.GetSize());
 
352
                        
 
353
                if(!AskYesNo("    Add another observation (y/n)? [no] ", false))
 
354
                        break;
 
355
                mprintf("\n");
 
356
        }
 
357
        
 
358
        return true;
 
359
}
 
360
 
 
361
bool initializePDF() {
 
362
        int i;
 
363
        for(i = 0; i < g_pdfObserv.GetSize(); i++) {
 
364
                ((CPDFObservation *)g_pdfObserv[i])->initialize();
 
365
        }
 
366
        return true;
 
367
}
 
368
 
 
369
void processPDF(CTimeStep *ts) {
 
370
        int i;
 
371
        for(i = 0; i < g_pdfObserv.GetSize(); i++) {
 
372
                ((CPDFObservation *)g_pdfObserv[i])->process(ts);
 
373
        }
 
374
}
 
375
 
 
376
void finalizePDF() {
 
377
        int i;
 
378
        for(i = 0; i < g_pdfObserv.GetSize(); i++) {
 
379
                mprintf(YELLOW, "\n>>> Density Profile Observation %d >>>\n\n", i+1);
 
380
                ((CPDFObservation *)g_pdfObserv[i])->finalize();
 
381
                mprintf(YELLOW, "\n<<< End of Density Profile Observation %d <<<\n\n", i+1);
 
382
        }
 
383
        
 
384
        for(i = 0; i < g_pdfObserv.GetSize(); i++) {
 
385
                delete (CPDFObservation *)g_pdfObserv[i];
 
386
        }
 
387
}