« back to all changes in this revision
Viewing changes to src/raman.cpp
- browse files at revision 12
- compare with another revision
- download diff
- download tarball
- view history from revision 12
- 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.
- files added:
- src/plproj.cpp
- files modified:
- Makefile
added
removed
src/raman.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 "raman.h"
25
26
#include "acf.h"
27
#include "df.h"
28
#include "globalvar.h"
29
#include "maintools.h"
30
#include "moltools.h"
31
#include "timestep.h"
32
#include "tools.h"
33
#include "xfloatarray.h"
34
#include "xobarray.h"
35
36
#include <errno.h>
37
#include <stdio.h>
38
#include <string.h>
39
40
#ifdef TARGET_LINUX
41
#include <sys/stat.h>
42
#include <sys/types.h>
43
#endif
44
45
#define BUF_SIZE 4096
46
47
static bool g_newRaman;
48
static char *g_ramanDir;
49
static bool g_orientAvg;
50
static float g_fieldStrength;
51
static int g_step = 0;
52
static int g_stride;
53
static float g_laser;
54
static float g_temp;
55
static char *g_inputTemplate;
56
static char *g_templateFieldPos;
57
static char *g_templatePolPos;
58
static char *g_templateCoordPos;
59
static char *g_templateStepsPos;
60
61
static CxObArray g_ramObserv;
62
63
static FILE *g_polFile[3];
64
static CTimeStep *g_timestep[3];
65
66
static FILE *g_reftrajFile;
67
static int g_steps = 0;
68
69
CRamanDyn::CRamanDyn(int showMol, bool global) {
70
_global = global;
71
if(_global) {
72
m_iShowMol = -1;
73
m_iMolecules = g_oaSingleMolecules.GetSize();
74
} else {
75
m_iShowMol = showMol;
76
m_iMolecules = ((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex.GetSize();
77
}
78
79
/*** Set Names of dynamic arrays. By M. Brehm */
80
_dipole0.SetName("CRamanDyn::_dipole0");
81
char tbuf[256];
82
for (int tz1=0;tz1<3;tz1++)
83
{
84
for (int tz2=0;tz2<3;tz2++)
85
{
86
sprintf(tbuf,"CRamanDyn::_polarizability[%d][%d]",tz1,tz2);
87
_polarizability[tz1][tz2].SetName(tbuf);
88
}
89
}
90
/* End of set names */
91
92
mprintf(YELLOW, ">>> Raman Spectrum >>>\n\n");
93
94
if(_global)
95
mprintf(" All atoms will be taken.\n\n");
96
else
97
mprintf(" All atoms will be taken from the OM %s.\n\n", ((CMolecule *)g_oaMolecules[m_iShowMol])->m_sName);
98
99
m_iVecType = 1;
100
m_iCombinations = 1;
101
g_bDipole = true;
102
ParseDipole();
103
104
if(g_iTrajSteps != -1) {
105
int depth = (int)(g_iTrajSteps / g_stride * 0.75);
106
if(depth > 4096)
107
depth = 4096;
108
m_iDepth = AskUnsignedInteger(" Enter the resolution (=depth) of the ACF (in time steps): [%d] ", depth, depth);
109
} else {
110
m_iDepth = AskUnsignedInteger(" Enter the resolution (=depth) of the ACF (in time steps): [2000] ", 2000);
111
}
112
113
int size = CalcFFTSize(m_iDepth, false);
114
if(m_iDepth != size) {
115
mprintf(WHITE," The next \"fast\" size for FFT is %d. Using this instead of %d as depth.\n", size, m_iDepth);
116
m_iDepth = size;
117
}
118
119
m_iStride = 1;
120
m_bSpectrum = true;
121
122
_derivativeOrder = 1;
123
bool window = true;
124
bool derive = true;
125
if(g_bAdvanced2) {
126
derive = AskYesNo(" Derive the vectors before autocorrelation (y/n)? [yes] ", true);
127
if(derive)
128
_derivativeOrder = AskRangeInteger(" Please enter degree of vector derivative (1-2): [1] ", 1, 2, 1);
129
window = AskYesNo(" Apply window function (Cos^2) to autocorrelation function (y/n)? [yes] ", true);
130
}
131
float possibleRange = 33356.41f / g_fTimestepLength / 2.0f / g_stride;
132
mprintf(" A time step length of %.2f fs with a stride of %d allows a spectral range up to %.2f cm^-1.\n", g_fTimestepLength, g_stride, possibleRange);
133
float specWaveNumber = AskRangeFloat("\n Calculate spectrum up to which wave number (cm^-1)? [%.2f cm^-1] ", 0, possibleRange, (possibleRange < 5000.0f) ? possibleRange : 5000.0f, (possibleRange < 5000.0f) ? possibleRange : 5000.0f);
134
int mirror = 1;
135
int zeroPadding = m_iDepth * 3;
136
if(g_bAdvanced2) {
137
mirror = AskRangeInteger(" No mirroring (0) or short-time enhancing (1)? [1] ", 0, 1, 1);
138
zeroPadding = AskUnsignedInteger(" Zero Padding: How many zeros to insert? [%d] ", m_iDepth * 3, m_iDepth * 3);
139
}
140
141
size = CalcFFTSize(m_iDepth + zeroPadding, false);
142
if(m_iDepth + zeroPadding != size) {
143
mprintf(WHITE, " The next \"fast\" size for FFT is %d. Using %d zeros for zero padding.\n", size, size-m_iDepth);
144
zeroPadding = size-m_iDepth;
145
}
146
147
mprintf(" This results in a spectral resolution to %.2f cm^-1.\n\n", 33356.41 / g_fTimestepLength / 2.0f / size);
148
149
try { _isoACF = new CACF(); } catch(...) { _isoACF = NULL; }
150
if(_isoACF == NULL) NewException((double)sizeof(CACF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
151
_isoACF->m_iSize = m_iDepth;
152
_isoACF->m_bSpectrum = true;
153
_isoACF->m_bDerivative = derive;
154
_isoACF->m_iDerivative = _derivativeOrder;
155
_isoACF->m_bWindowFunction = window;
156
_isoACF->m_fSpecWaveNumber = specWaveNumber;
157
_isoACF->m_iMirror = mirror;
158
_isoACF->m_iZeroPadding = zeroPadding;
159
_isoACF->m_bACF_DB = false;
160
161
try { _anisoACF = new CACF(); } catch(...) { _anisoACF = NULL; }
162
if(_anisoACF == NULL) NewException((double)sizeof(CACF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
163
_anisoACF->m_iSize = m_iDepth;
164
_anisoACF->m_bSpectrum = true;
165
_anisoACF->m_bDerivative = derive;
166
_anisoACF->m_iDerivative = _derivativeOrder;
167
_anisoACF->m_bWindowFunction = window;
168
_anisoACF->m_fSpecWaveNumber = specWaveNumber;
169
_anisoACF->m_iMirror = mirror;
170
_anisoACF->m_iZeroPadding = zeroPadding;
171
_anisoACF->m_bACF_DB = false;
172
173
BuildName();
174
175
mprintf(YELLOW, "<<< End of Raman Spectrum <<<\n\n");
176
}
177
178
CRamanDyn::~CRamanDyn() {
179
int i, j, k;
180
181
for(i = 0; i < m_iMolecules; i++)
182
delete (CxFloatArray *)_dipole0[i];
183
for(i = 0; i < m_iMolecules; i++) {
184
for(j = 0; j < 3; j++) {
185
for(k = 0; k < (g_orientAvg ? 3 : 1); k++) {
186
delete (CxFloatArray *)_polarizability[j][k][i];
187
}
188
}
189
}
190
delete _isoACF;
191
delete _anisoACF;
192
}
193
194
void CRamanDyn::initialize() {
195
int i, j, k;
196
197
_isoACF->Create();
198
_anisoACF->Create();
199
200
if (g_iTrajSteps != -1)
201
mprintf(" Raman Cache: Trying to allocate %s of memory...\n", FormatBytes((double)m_iMolecules * g_iTrajSteps / g_iStride / g_stride * (g_orientAvg ? 9.9 : 3.3) * sizeof(float)));
202
else
203
mprintf(" Raman Cache: Trying to allocate %s of memory...\n", FormatBytes((double)m_iMolecules * 2000 / g_iStride /g_stride * (g_orientAvg ? 9.9 : 3.3) * sizeof(float)));
204
for(i = 0; i < m_iMolecules; i++) {
205
CxVector3 *dipoleVector;
206
try { dipoleVector = new CxVector3(); } catch(...) { dipoleVector = NULL; }
207
if(dipoleVector == NULL) NewException((double)sizeof(CxVector3), __FILE__, __LINE__, __PRETTY_FUNCTION__);
208
_dipole0.Add(dipoleVector);
209
}
210
for(i = 0; i < m_iMolecules; i++) {
211
for(j = 0; j < 3; j++) {
212
for(k = 0; k < (g_orientAvg ? 3 : 1); k++) {
213
CxFloatArray *floatArray;
214
try { floatArray = new CxFloatArray("CRamanDyn::initialize():floatArray"); } catch(...) { floatArray = NULL; }
215
if(floatArray == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
216
if(g_iTrajSteps != -1) {
217
floatArray->SetMaxSize((long)(g_iTrajSteps / g_iStride / g_stride * 1.1));
218
floatArray->SetGrow((long)(g_iTrajSteps / g_iStride / g_stride * 0.1));
219
} else {
220
floatArray->SetGrow(10000);
221
}
222
_polarizability[j][k].Add(floatArray);
223
}
224
}
225
}
226
}
227
228
void CRamanDyn::getDipole0() {
229
int i;
230
231
for(i = 0; i < m_iMolecules; i++) {
232
CSingleMolecule *sm;
233
if(_global)
234
sm = (CSingleMolecule *)g_oaSingleMolecules[i];
235
else
236
sm = (CSingleMolecule *)g_oaSingleMolecules[((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex[i]];
237
*((CxVector3 *)_dipole0[i]) = sm->m_vDipole;
238
}
239
}
240
241
void CRamanDyn::calcPolarizability(int fieldDirection) {
242
int i;
243
244
for(i = 0; i < m_iMolecules; i++) {
245
CSingleMolecule *sm;
246
if(_global)
247
sm = (CSingleMolecule *)g_oaSingleMolecules[i];
248
else
249
sm = (CSingleMolecule *)g_oaSingleMolecules[((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex[i]];
250
for(int j = 0; j < 3; j++) {
251
((CxFloatArray *)_polarizability[j][fieldDirection][i])->Add(((*((CxVector3 *)_dipole0[i]))[j] - sm->m_vDipole[j]) / g_fieldStrength * 0.393430f); // 0.393430 - Umrechnung von Debye in a.u.
252
}
253
}
254
}
255
256
void CRamanDyn::finalize() {
257
int i, j, k, l;
258
259
char filename[BUF_SIZE];
260
#ifdef TARGET_LINUX
261
snprintf(filename, BUF_SIZE, "%s/polarizability_%s.dat", g_ramanDir, m_sName);
262
#else
263
sprintf(filename, "%s/polarizability_%s.dat", g_ramanDir, m_sName);
264
#endif
265
mprintf(" Writing polarizabilities for first molecule to \"%s\"...\n", filename);
266
FILE *pol_file;
267
pol_file = OpenFileWrite(filename, false);
268
for(i = 0; i < ((CxFloatArray *)_polarizability[0][0][0])->GetSize(); i++) {
269
fprintf(pol_file, "%10.2f", i * g_fTimestepLength * g_iStride * g_stride);
270
for(j = 0; j < (g_orientAvg ? 3 : 1); j++) {
271
for(k = 0; k < 3; k++) {
272
fprintf(pol_file, " %14.8f", (*((CxFloatArray *)_polarizability[k][j][0]))[i]);
273
}
274
}
275
fprintf(pol_file, "\n");
276
}
277
fclose(pol_file);
278
279
float step;
280
switch(_derivativeOrder) {
281
case 0:
282
mprintf(" Not deriving polarizabilities.\n");
283
break;
284
case 1:
285
mprintf(" Deriving polarizabilities (1st derivative)...\n");
286
mprintf(WHITE, " [");
287
step = m_iMolecules / 20.0f;
288
for(i = 0; i < m_iMolecules; i++) {
289
if(fmodf((float)i, step) < 1.0f)
290
mprintf(WHITE, "#");
291
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - 1; j++) {
292
for(k = 0; k < 3; k++) {
293
for(l = 0; l < (g_orientAvg ? 3 : 1); l++) {
294
(*((CxFloatArray *)_polarizability[k][l][i]))[j] = 0.5f * ((*((CxFloatArray *)_polarizability[k][l][i]))[j+1] - (*((CxFloatArray *)_polarizability[k][l][i]))[j]);
295
}
296
}
297
}
298
}
299
mprintf(WHITE, "]\n");
300
break;
301
case 2:
302
mprintf(" Deriving polarizabilities (2nd derivative)...\n");
303
mprintf(WHITE, " [");
304
step = m_iMolecules / 20.0f;
305
for(i = 0; i < m_iMolecules; i++) {
306
if(fmodf((float)i, step) < 1.0f)
307
mprintf(WHITE, "#");
308
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - 2; j++) {
309
for(k = 0; k < 3; k++) {
310
for(l = 0; l < (g_orientAvg ? 3 : 1); l++) {
311
(*((CxFloatArray *)_polarizability[k][l][i]))[j] = (*((CxFloatArray *)_polarizability[k][l][i]))[j+2] - 2.0f * (*((CxFloatArray *)_polarizability[k][l][i]))[j+1] + (*((CxFloatArray *)_polarizability[k][l][i]))[j];
312
}
313
}
314
}
315
}
316
mprintf(WHITE, "]\n");
317
break;
318
default:
319
mprintf(RED, "Higher derivatives not implemented.\n");
320
abort();
321
break;
322
}
323
324
mprintf(" Processing polarizability tensor components...\n");
325
step = m_iMolecules / 20.0f;
326
327
CAutoCorrelation *autoCorrelation;
328
try { autoCorrelation = new CAutoCorrelation(); } catch(...) { autoCorrelation = NULL; }
329
if(autoCorrelation == NULL) NewException((double)sizeof(CAutoCorrelation), __FILE__, __LINE__, __PRETTY_FUNCTION__);
330
autoCorrelation->Init(((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder, m_iDepth, g_bACFFFT);
331
332
if(g_orientAvg) {
333
CxFloatArray *isotropyPol, *isotropyACF;
334
try { isotropyACF = new CxFloatArray("CRamanDyn::finalize():isotropyACF"); } catch(...) { isotropyACF = NULL; }
335
if(isotropyACF == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
336
isotropyACF->SetSize(m_iDepth);
337
try { isotropyPol = new CxFloatArray("CRamanDyn::finalize():isotropyPol"); } catch(...) { isotropyPol = NULL; }
338
if(isotropyPol == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
339
isotropyPol->SetSize(((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder);
340
341
mprintf(" Isotropic part...\n");
342
mprintf(WHITE, " [");
343
for(i = 0; i < m_iDepth; i++) {
344
_isoACF->m_pData[i] = 0.0f;
345
}
346
for(i = 0; i < m_iMolecules; i++) {
347
if(fmodf((float)i, step) < 1.0f)
348
mprintf(WHITE, "#");
349
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
350
(*isotropyPol)[j] = 0.0f;
351
for(k = 0; k < 3; k++) {
352
(*isotropyPol)[j] += (*((CxFloatArray *)_polarizability[k][k][i]))[j];
353
}
354
(*isotropyPol)[j] /= 3.0f;
355
}
356
autoCorrelation->AutoCorrelate(isotropyPol, isotropyACF);
357
for(j = 0; j < m_iDepth; j++) {
358
_isoACF->m_pData[j] += (*isotropyACF)[j];
359
}
360
}
361
for(i = 0; i < m_iDepth; i++) {
362
if(_global)
363
_isoACF->m_pData[i] /= g_iSteps / g_stride;
364
else
365
_isoACF->m_pData[i] /= g_iSteps / g_stride * m_iMolecules;
366
}
367
mprintf(WHITE, "]\n");
368
369
// #ifdef TARGET_LINUX
370
// snprintf(filename, BUF_SIZE, "%s/acf_iso_%s.dat", g_ramanDir, m_sName);
371
// #else
372
// sprintf(filename, "%s/acf_iso_%s.dat", g_ramanDir, m_sName);
373
// #endif
374
// mprintf(" Saving ACF as %s...\n", filename);
375
// FILE *acf_file = OpenFileWrite(filename, false);
376
// for(i = 0; i < m_iDepth; i++)
377
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _isoACF->m_pData[i]);
378
// fclose(acf_file);
379
380
if(_isoACF->m_iMirror != 0) {
381
mprintf(" Mirroring ACF...\n");
382
_isoACF->Mirror(_isoACF->m_iMirror);
383
}
384
if(_isoACF->m_bWindowFunction != 0) {
385
mprintf(" Applying window function to ACF...\n");
386
_isoACF->Window();
387
}
388
389
// #ifdef TARGET_LINUX
390
// snprintf(filename, BUF_SIZE, "%s/acf_iso_%s.mw.dat", g_ramanDir, m_sName);
391
// #else
392
// sprintf(filename, "%s/acf_iso_%s.mw.dat", g_ramanDir, m_sName);
393
// #endif
394
// mprintf(" Saving ACF as %s...\n", filename);
395
// acf_file = OpenFileWrite(filename, false);
396
// for(i = 0; i < _isoACF->m_iSize; i++)
397
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _isoACF->m_pData[i]);
398
// fclose(acf_file);
399
400
mprintf(" Performing Fourier transformation...\n");
401
CFFT *fft;
402
try { fft = new CFFT(); } catch(...) { fft = NULL; }
403
if(fft == NULL) NewException((double)sizeof(CFFT), __FILE__, __LINE__, __PRETTY_FUNCTION__);
404
fft->PrepareFFT_C2C(_isoACF->m_iSize + _isoACF->m_iZeroPadding);
405
_isoACF->Transform(fft);
406
delete fft;
407
_isoACF->m_pSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
408
409
// #ifdef TARGET_LINUX
410
// snprintf(filename, BUF_SIZE, "%s/spectrum_iso_%s.dat", g_ramanDir, m_sName);
411
// #else
412
// sprintf(filename, "%s/spectrum_iso_%s.dat", g_ramanDir, m_sName);
413
// #endif
414
// mprintf(" Saving spectrum as %s...\n", filename);
415
// _isoACF->m_pSpectrum->Write("", filename, "");
416
417
delete isotropyACF;
418
delete isotropyPol;
419
420
CxFloatArray *anisotropyPol, *anisotropyACF, *anisotropyACFSum;
421
try { anisotropyPol = new CxFloatArray("CRamanDyn::finalize():anisotropyPol"); } catch(...) { anisotropyPol = NULL; }
422
if(anisotropyPol == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
423
anisotropyPol->SetSize(((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder);
424
try { anisotropyACF = new CxFloatArray("CRamanDyn::finalize():anisotropyACF"); } catch(...) { anisotropyACF = NULL; }
425
if(anisotropyACF == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
426
anisotropyACF->SetSize(m_iDepth);
427
try { anisotropyACFSum = new CxFloatArray("CRamanDyn::finalize():anisotropyACFSum"); } catch(...) { anisotropyACFSum = NULL; }
428
if(anisotropyACFSum == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
429
anisotropyACFSum->SetSize(m_iDepth);
430
431
mprintf(" First anisotropic part...\n");
432
mprintf(WHITE, " [");
433
for(i = 0; i < m_iDepth; i++) {
434
_anisoACF->m_pData[i] = 0.0f;
435
}
436
for(i = 0; i < m_iMolecules; i++) {
437
if(fmodf((float)i, step) < 1.0f)
438
mprintf(WHITE, "#");
439
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
440
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[0][0][i]))[j];
441
(*anisotropyPol)[j] -= (*((CxFloatArray *)_polarizability[1][1][i]))[j];
442
}
443
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
444
for(int j = 0; j < m_iDepth; j++) {
445
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
446
}
447
}
448
for(i = 0; i < m_iDepth; i++) {
449
(*anisotropyACFSum)[i] = 0.5f * _anisoACF->m_pData[i];
450
}
451
mprintf(WHITE, "]\n");
452
453
mprintf(" Second anisotropic part...\n");
454
mprintf(WHITE, " [");
455
for(i = 0; i < m_iDepth; i++) {
456
_anisoACF->m_pData[i] = 0.0f;
457
}
458
for(i = 0; i < m_iMolecules; i++) {
459
if(fmodf((float)i, step) < 1.0f)
460
mprintf(WHITE, "#");
461
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
462
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[1][1][i]))[j];
463
(*anisotropyPol)[j] -= (*((CxFloatArray *)_polarizability[2][2][i]))[j];
464
}
465
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
466
for(j = 0; j < m_iDepth; j++) {
467
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
468
}
469
}
470
for(i = 0; i < m_iDepth; i++) {
471
(*anisotropyACFSum)[i] += 0.5f * _anisoACF->m_pData[i];
472
}
473
mprintf(WHITE, "]\n");
474
475
mprintf(" Third anisotropic part...\n");
476
mprintf(WHITE, " [");
477
for(i = 0; i < m_iDepth; i++) {
478
_anisoACF->m_pData[i] = 0.0f;
479
}
480
for(i = 0; i < m_iMolecules; i++) {
481
if(fmodf((float)i, step) < 1.0f)
482
mprintf(WHITE, "#");
483
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
484
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[2][2][i]))[j];
485
(*anisotropyPol)[j] -= (*((CxFloatArray *)_polarizability[0][0][i]))[j];
486
}
487
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
488
for(j = 0; j < m_iDepth; j++) {
489
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
490
}
491
}
492
for(i = 0; i < m_iDepth; i++) {
493
(*anisotropyACFSum)[i] += 0.5f * _anisoACF->m_pData[i];
494
}
495
mprintf(WHITE, "]\n");
496
497
mprintf(" Fourth anisotropic part...\n");
498
mprintf(WHITE, " [");
499
for(i = 0; i < m_iDepth; i++) {
500
_anisoACF->m_pData[i] = 0.0f;
501
}
502
for(i = 0; i < m_iMolecules; i++) {
503
if(fmodf((float)i, step) < 1.0f)
504
mprintf(WHITE, "#");
505
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
506
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[0][1][i]))[j];
507
(*anisotropyPol)[j] += (*((CxFloatArray *)_polarizability[1][0][i]))[j];
508
(*anisotropyPol)[j] *= 0.5f;
509
}
510
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
511
for(j = 0; j < m_iDepth; j++) {
512
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
513
}
514
}
515
for(i = 0; i < m_iDepth; i++) {
516
(*anisotropyACFSum)[i] += 3.0f * _anisoACF->m_pData[i];
517
}
518
mprintf(WHITE, "]\n");
519
520
mprintf(" Fifth anisotropic part...\n");
521
mprintf(WHITE, " [");
522
for(i = 0; i < m_iDepth; i++) {
523
_anisoACF->m_pData[i] = 0.0f;
524
}
525
for(i = 0; i < m_iMolecules; i++) {
526
if(fmodf((float)i, step) < 1.0f)
527
mprintf(WHITE, "#");
528
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
529
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[1][2][i]))[j];
530
(*anisotropyPol)[j] += (*((CxFloatArray *)_polarizability[2][1][i]))[j];
531
(*anisotropyPol)[j] *= 0.5f;
532
}
533
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
534
for(j = 0; j < m_iDepth; j++) {
535
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
536
}
537
}
538
for(i = 0; i < m_iDepth; i++) {
539
(*anisotropyACFSum)[i] += 3.0f * _anisoACF->m_pData[i];
540
}
541
mprintf(WHITE, "]\n");
542
543
mprintf(" Sixth anisotropic part...\n");
544
mprintf(WHITE, " [");
545
for(i = 0; i < m_iDepth; i++) {
546
_anisoACF->m_pData[i] = 0.0f;
547
}
548
for(i = 0; i < m_iMolecules; i++) {
549
if(fmodf((float)i, step) < 1.0f)
550
mprintf(WHITE, "#");
551
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
552
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[2][0][i]))[j];
553
(*anisotropyPol)[j] += (*((CxFloatArray *)_polarizability[0][2][i]))[j];
554
(*anisotropyPol)[j] *= 0.5f;
555
}
556
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
557
for(j = 0; j < m_iDepth; j++) {
558
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
559
}
560
}
561
for(i = 0; i < m_iDepth; i++) {
562
(*anisotropyACFSum)[i] += 3.0f * _anisoACF->m_pData[i];
563
}
564
mprintf(WHITE, "]\n");
565
566
for(i = 0; i < m_iDepth; i++) {
567
if(_global)
568
_anisoACF->m_pData[i] = (*anisotropyACFSum)[i] / g_iSteps * g_stride;
569
else
570
_anisoACF->m_pData[i] = (*anisotropyACFSum)[i] / g_iSteps * g_stride / m_iMolecules;
571
}
572
573
// #ifdef TARGET_LINUX
574
// snprintf(filename, BUF_SIZE, "%s/acf_aniso_%s.dat", g_ramanDir, m_sName);
575
// #else
576
// sprintf(filename, "%s/acf_aniso_%s.dat", g_ramanDirm_sName);
577
// #endif
578
// mprintf(" Saving ACF as %s...\n", filename);
579
// acf_file = OpenFileWrite(filename, false);
580
// for(i = 0; i < m_iDepth; i++)
581
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _anisoACF->m_pData[i]);
582
// fclose(acf_file);
583
584
if(_anisoACF->m_iMirror != 0) {
585
mprintf(" Mirroring ACF...\n");
586
_anisoACF->Mirror(_anisoACF->m_iMirror);
587
}
588
if(_anisoACF->m_bWindowFunction != 0) {
589
mprintf(" Applying window function to ACF...\n");
590
_anisoACF->Window();
591
}
592
593
// #ifdef TARGET_LINUX
594
// snprintf(filename, BUF_SIZE, "%s/acf_aniso_%s.mw.dat", g_ramanDir, m_sName);
595
// #else
596
// sprintf(filename, "%s/acf_aniso_%s.mw.dat", g_ramanDir, m_sName);
597
// #endif
598
// mprintf(" Saving ACF as %s...\n", filename);
599
// acf_file = OpenFileWrite(filename, false);
600
// for(i = 0; i < _anisoACF->m_iSize; i++)
601
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _anisoACF->m_pData[i]);
602
// fclose(acf_file);
603
604
mprintf(" Performing Fourier transformation...\n");
605
try { fft = new CFFT(); } catch(...) { fft = NULL; }
606
if(fft == NULL) NewException((double)sizeof(CFFT), __FILE__, __LINE__, __PRETTY_FUNCTION__);
607
fft->PrepareFFT_C2C(_anisoACF->m_iSize + _anisoACF->m_iZeroPadding);
608
_anisoACF->Transform(fft);
609
delete fft;
610
_anisoACF->m_pSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
611
612
// #ifdef TARGET_LINUX
613
// snprintf(filename, BUF_SIZE, "%s/spectrum_aniso_%s.dat", g_ramanDir, m_sName);
614
// #else
615
// sprintf(filename, "%s/spectrum_aniso_%s.dat", g_ramanDir, m_sName);
616
// #endif
617
// mprintf(" Saving spectrum as %s...\n", filename);
618
// _anisoACF->m_pSpectrum->Write("", filename, "");
619
620
int specSize = _isoACF->m_pSpectrum->m_iSize;
621
CSpectrum *paraSpectrum, *orthoSpectrum, *sumSpectrum, *depolSpectrum;
622
try { paraSpectrum = new CSpectrum(); } catch(...) { paraSpectrum = NULL; }
623
if(paraSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
624
paraSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
625
paraSpectrum->Create(specSize);
626
paraSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
627
try { orthoSpectrum = new CSpectrum(); } catch(...) { orthoSpectrum = NULL; }
628
if(orthoSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
629
orthoSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
630
orthoSpectrum->Create(specSize);
631
orthoSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
632
try { sumSpectrum = new CSpectrum(); } catch(...) { sumSpectrum = NULL; }
633
if(sumSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
634
sumSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
635
sumSpectrum->Create(specSize);
636
sumSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
637
try { depolSpectrum = new CSpectrum(); } catch(...) { depolSpectrum = NULL; }
638
if(depolSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
639
depolSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
640
depolSpectrum->Create(specSize);
641
depolSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
642
643
for(i = 0; i < specSize; i++) {
644
paraSpectrum->m_pData[i] = _isoACF->m_pSpectrum->m_pData[i] + 4.0f / 45.0f * _anisoACF->m_pSpectrum->m_pData[i];
645
}
646
647
// #ifdef TARGET_LINUX
648
// snprintf(filename, BUF_SIZE, "%s/spectrum_para_%s.dat", g_ramanDir, m_sName);
649
// #else
650
// sprintf(filename, "%s/spectrum_para_%s.dat", g_ramanDir, m_sName);
651
// #endif
652
// mprintf(" Saving para spectrum as %s...\n", filename);
653
// paraSpectrum->Write("", filename, "");
654
655
for(i = 0; i < specSize; i++) {
656
orthoSpectrum->m_pData[i] = _anisoACF->m_pSpectrum->m_pData[i] / 15.0f;
657
}
658
659
// #ifdef TARGET_LINUX
660
// snprintf(filename, BUF_SIZE, "%s/spectrum_ortho_%s.dat", g_ramanDir, m_sName);
661
// #else
662
// sprintf(filename, "%s/spectrum_ortho_%s.dat", g_ramanDir, m_sName);
663
// #endif
664
// mprintf(" Saving ortho spectrum as %s...\n", filename);
665
// orthoSpectrum->Write("", filename, "");
666
667
for(i = 0; i < specSize; i++) {
668
sumSpectrum->m_pData[i] = paraSpectrum->m_pData[i] + orthoSpectrum->m_pData[i];
669
}
670
671
// #ifdef TARGET_LINUX
672
// snprintf(filename, BUF_SIZE, "%s/spectrum_unpol_%s.dat", g_ramanDir, m_sName);
673
// #else
674
// sprintf(filename, "%s/spectrum_unpol_%s.dat", g_ramanDir, m_sName);
675
// #endif
676
// mprintf(" Saving unpolarized spectrum as %s...\n", filename);
677
// sumSpectrum->Write("", filename, "");
678
679
for(i = 0; i < specSize; i++) {
680
float freq = i * paraSpectrum->m_fMaxRWL / paraSpectrum->m_iSize;
681
float crossSecFactor = powf(g_laser - freq, 4) / freq / (1 - expf(-1.438777f * freq / g_temp));
682
paraSpectrum->m_pData[i] *= crossSecFactor;
683
orthoSpectrum->m_pData[i] *= crossSecFactor;
684
sumSpectrum->m_pData[i] *= crossSecFactor;
685
}
686
687
for(i = 0; i < specSize; i++) {
688
depolSpectrum->m_pData[i] = orthoSpectrum->m_pData[i] / paraSpectrum->m_pData[i];
689
}
690
691
#ifdef TARGET_LINUX
692
snprintf(filename, BUF_SIZE, "%s/spectrum_para_%s.csv", g_ramanDir, m_sName);
693
#else
694
sprintf(filename, "%s/spectrum_para_%s.csv", g_ramanDir, m_sName);
695
#endif
696
mprintf(" Saving parallel spectrum as %s...\n", filename);
697
paraSpectrum->Write("", filename, "");
698
699
#ifdef TARGET_LINUX
700
snprintf(filename, BUF_SIZE, "%s/spectrum_ortho_%s.csv", g_ramanDir, m_sName);
701
#else
702
sprintf(filename, "%s/spectrum_ortho_%s.csv", g_ramanDir, m_sName);
703
#endif
704
mprintf(" Saving orthogonal spectrum as %s...\n", filename);
705
orthoSpectrum->Write("", filename, "");
706
707
#ifdef TARGET_LINUX
708
snprintf(filename, BUF_SIZE, "%s/spectrum_unpol_%s.csv", g_ramanDir, m_sName);
709
#else
710
sprintf(filename, "%s/spectrum_unpol_%s.csv", g_ramanDir, m_sName);
711
#endif
712
mprintf(" Saving unpolarized spectrum as %s...\n", filename);
713
sumSpectrum->Write("", filename, "");
714
715
#ifdef TARGET_LINUX
716
snprintf(filename, BUF_SIZE, "%s/depol_ratio_%s.csv", g_ramanDir, m_sName);
717
#else
718
sprintf(filename, "%s/depol_ratio_%s.csv", g_ramanDir, m_sName);
719
#endif
720
mprintf(" Saving depolarization ratio as %s...\n", filename);
721
depolSpectrum->Write("", filename, "");
722
723
delete anisotropyACFSum;
724
delete anisotropyACF;
725
delete anisotropyPol;
726
727
delete paraSpectrum;
728
delete orthoSpectrum;
729
delete sumSpectrum;
730
delete depolSpectrum;
731
} else {
732
CxFloatArray *ACF;
733
try { ACF = new CxFloatArray("CRamanDyn::finalize():ACF"); } catch(...) { ACF = NULL; }
734
if(ACF == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
735
ACF->SetSize(m_iDepth);
736
737
mprintf(" Parallel part...\n");
738
mprintf(WHITE, " [");
739
for(i = 0; i < m_iDepth; i++) {
740
_isoACF->m_pData[i] = 0.0f;
741
}
742
for(i = 0; i < m_iMolecules; i++) {
743
if(fmodf((float)i, step) < 1.0f)
744
mprintf(WHITE, "#");
745
autoCorrelation->AutoCorrelate((CxFloatArray *)_polarizability[0][0][i], ACF);
746
for(int j = 0; j < m_iDepth; j++) {
747
_isoACF->m_pData[j] += (*ACF)[j];
748
}
749
}
750
for(i = 0; i < m_iDepth; i++) {
751
if(_global)
752
_isoACF->m_pData[i] /= g_iSteps / g_stride;
753
else
754
_isoACF->m_pData[i] /= g_iSteps / g_stride * m_iMolecules;
755
}
756
mprintf(WHITE, "]\n");
757
758
// #ifdef TARGET_LINUX
759
// snprintf(filename, BUF_SIZE, "%s/acf_para_%s.dat", g_ramanDir, m_sName);
760
// #else
761
// sprintf(filename, "%s/acf_para_%s.dat", g_ramanDir, m_sName);
762
// #endif
763
// mprintf(" Saving ACF as %s...\n", filename);
764
// FILE *acf_file = OpenFileWrite(filename, false);
765
// for(i = 0; i < m_iDepth; i++)
766
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _isoACF->m_pData[i]);
767
// fclose(acf_file);
768
769
if(_isoACF->m_iMirror != 0) {
770
mprintf(" Mirroring ACF...\n");
771
_isoACF->Mirror(_isoACF->m_iMirror);
772
}
773
if(_isoACF->m_bWindowFunction != 0) {
774
mprintf(" Applying window function to ACF...\n");
775
_isoACF->Window();
776
}
777
778
// #ifdef TARGET_LINUX
779
// snprintf(filename, BUF_SIZE, "%s/acf_para_%s.mw.dat", g_ramanDir, m_sName);
780
// #else
781
// sprintf(filename, "%s/acf_para_%s.mw.dat", g_ramanDir, m_sName);
782
// #endif
783
// mprintf(" Saving ACF as %s...\n", filename);
784
// acf_file = OpenFileWrite(filename, false);
785
// for(i = 0; i < _isoACF->m_iSize; i++)
786
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _isoACF->m_pData[i]);
787
// fclose(acf_file);
788
789
mprintf(" Performing Fourier transformation...\n");
790
CFFT *fft;
791
try { fft = new CFFT(); } catch(...) { fft = NULL; }
792
if(fft == NULL) NewException((double)sizeof(CFFT), __FILE__, __LINE__, __PRETTY_FUNCTION__);
793
fft->PrepareFFT_C2C(_isoACF->m_iSize + _isoACF->m_iZeroPadding);
794
_isoACF->Transform(fft);
795
delete fft;
796
_isoACF->m_pSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
797
798
// #ifdef TARGET_LINUX
799
// snprintf(filename, BUF_SIZE, "%s/spectrum_para_%s.dat", g_ramanDir, m_sName);
800
// #else
801
// sprintf(filename, "%s/spectrum_para_%s.dat", g_ramanDir, m_sName);
802
// #endif
803
// mprintf(" Saving spectrum as %s...\n", filename);
804
// _isoACF->m_pSpectrum->Write("", filename, "");
805
806
mprintf(" Orthogonal part...\n");
807
mprintf(WHITE, " [");
808
for(i = 0; i < m_iDepth; i++) {
809
_anisoACF->m_pData[i] = 0.0f;
810
}
811
for(i = 0; i < m_iMolecules; i++) {
812
if(fmodf((float)i, step) < 1.0f)
813
mprintf(WHITE, "#");
814
autoCorrelation->AutoCorrelate((CxFloatArray *)_polarizability[1][0][i], ACF);
815
for(int j = 0; j < m_iDepth; j++) {
816
_anisoACF->m_pData[j] += (*ACF)[j];
817
}
818
}
819
for(i = 0; i < m_iDepth; i++) {
820
if(_global)
821
_anisoACF->m_pData[i] /= g_iSteps / g_stride;
822
else
823
_anisoACF->m_pData[i] /= g_iSteps / g_stride * m_iMolecules;
824
}
825
mprintf(WHITE, "]\n");
826
827
// #ifdef TARGET_LINUX
828
// snprintf(filename, BUF_SIZE, "%s/acf_ortho_%s.dat", g_ramanDir, m_sName);
829
// #else
830
// sprintf(filename, "%s/acf_ortho_%s.dat", g_ramanDir, m_sName);
831
// #endif
832
// mprintf(" Saving ACF as %s...\n", filename);
833
// acf_file = OpenFileWrite(filename, false);
834
// for(i = 0; i < m_iDepth; i++)
835
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _anisoACF->m_pData[i]);
836
// fclose(acf_file);
837
838
if(_anisoACF->m_iMirror != 0) {
839
mprintf(" Mirroring ACF...\n");
840
_anisoACF->Mirror(_isoACF->m_iMirror);
841
}
842
if(_anisoACF->m_bWindowFunction != 0) {
843
mprintf(" Applying window function to ACF...\n");
844
_anisoACF->Window();
845
}
846
847
// #ifdef TARGET_LINUX
848
// snprintf(filename, BUF_SIZE, "%s/acf_ortho_%s.mw.dat", g_ramanDir, m_sName);
849
// #else
850
// sprintf(filename, "%s/acf_ortho_%s.mw.dat", g_ramanDir, m_sName);
851
// #endif
852
// mprintf(" Saving ACF as %s...\n", filename);
853
// acf_file = OpenFileWrite(filename, false);
854
// for(i = 0; i < _isoACF->m_iSize; i++)
855
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _anisoACF->m_pData[i]);
856
// fclose(acf_file);
857
858
mprintf(" Performing Fourier transformation...\n");
859
try { fft = new CFFT(); } catch(...) { fft = NULL; }
860
if(fft == NULL) NewException((double)sizeof(CFFT), __FILE__, __LINE__, __PRETTY_FUNCTION__);
861
fft->PrepareFFT_C2C(_anisoACF->m_iSize + _anisoACF->m_iZeroPadding);
862
_anisoACF->Transform(fft);
863
delete fft;
864
_anisoACF->m_pSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
865
866
// #ifdef TARGET_LINUX
867
// snprintf(filename, BUF_SIZE, "%s/spectrum_ortho_%s.dat", g_ramanDir, m_sName);
868
// #else
869
// sprintf(filename, "%s/spectrum_ortho_%s.dat", g_ramanDir, m_sName);
870
// #endif
871
// mprintf(" Saving spectrum as %s...\n", filename);
872
// _anisoACF->m_pSpectrum->Write("", filename, "");
873
874
delete ACF;
875
876
int specSize = _isoACF->m_pSpectrum->m_iSize;
877
CSpectrum *sumSpectrum, *depolSpectrum;
878
try { sumSpectrum = new CSpectrum(); } catch(...) { sumSpectrum = NULL; }
879
if(sumSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
880
sumSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
881
sumSpectrum->Create(specSize);
882
sumSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
883
try { depolSpectrum = new CSpectrum(); } catch(...) { depolSpectrum = NULL; }
884
if(depolSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
885
depolSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
886
depolSpectrum->Create(specSize);
887
depolSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
888
889
for(i = 0; i < specSize; i++) {
890
sumSpectrum->m_pData[i] = _isoACF->m_pSpectrum->m_pData[i] + _anisoACF->m_pSpectrum->m_pData[i];
891
}
892
893
// #ifdef TARGET_LINUX
894
// snprintf(filename, BUF_SIZE, "%s/spectrum_unpol_%s.dat", g_ramanDir, m_sName);
895
// #else
896
// sprintf(filename, "%s/spectrum_unpol_%s.dat", g_ramanDir, m_sName);
897
// #endif
898
// mprintf(" Saving unpolarized spectrum as %s...\n", filename);
899
// sumSpectrum->Write("", filename, "");
900
901
for(i = 0; i < specSize; i++) {
902
float freq = i * _isoACF->m_pSpectrum->m_fMaxRWL / _isoACF->m_pSpectrum->m_iSize;
903
float crossSecFactor = powf(g_laser - freq, 4) / freq / (1 - expf(-1.438777f * freq / g_temp));
904
_isoACF->m_pSpectrum->m_pData[i] *= crossSecFactor;
905
_anisoACF->m_pSpectrum->m_pData[i] *= crossSecFactor;
906
sumSpectrum->m_pData[i] *= crossSecFactor;
907
}
908
909
for(i = 0; i < specSize; i++) {
910
depolSpectrum->m_pData[i] = _anisoACF->m_pSpectrum->m_pData[i] / _isoACF->m_pSpectrum->m_pData[i];
911
}
912
913
#ifdef TARGET_LINUX
914
snprintf(filename, BUF_SIZE, "%s/spectrum_para_%s.csv", g_ramanDir, m_sName);
915
#else
916
sprintf(filename, "%s/spectrum_para_%s.csv", g_ramanDir, m_sName);
917
#endif
918
mprintf(" Saving parallel spectrum as %s...\n", filename);
919
_isoACF->m_pSpectrum->Write("", filename, "");
920
921
#ifdef TARGET_LINUX
922
snprintf(filename, BUF_SIZE, "%s/spectrum_ortho_%s.csv", g_ramanDir, m_sName);
923
#else
924
sprintf(filename, "%s/spectrum_ortho_%s.csv", g_ramanDir, m_sName);
925
#endif
926
mprintf(" Saving orthogonal spectrum as %s...\n", filename);
927
_anisoACF->m_pSpectrum->Write("", filename, "");
928
929
#ifdef TARGET_LINUX
930
snprintf(filename, BUF_SIZE, "%s/spectrum_unpol_%s.csv", g_ramanDir, m_sName);
931
#else
932
sprintf(filename, "%s/spectrum_unpol_%s.csv", g_ramanDir, m_sName);
933
#endif
934
mprintf(" Saving unpolarized spectrum as %s...\n", filename);
935
sumSpectrum->Write("", filename, "");
936
937
#ifdef TARGET_LINUX
938
snprintf(filename, BUF_SIZE, "%s/depol_ratio_%s.csv", g_ramanDir, m_sName);
939
#else
940
sprintf(filename, "%s/depol_ratio_%s.csv", g_ramanDir, m_sName);
941
#endif
942
mprintf(" Saving depolarization ratio as %s...\n", filename);
943
depolSpectrum->Write("", filename, "");
944
945
delete sumSpectrum;
946
delete depolSpectrum;
947
}
948
949
delete autoCorrelation;
950
}
951
952
CRamanObservation::CRamanObservation(bool global) {
953
int i;
954
955
m_bTimeDev = false;
956
m_bTimeDiff = false;
957
m_bSelf = false;
958
m_bOthers = true;
959
960
if(global) {
961
m_iShowMol = -1;
962
m_iShowMolCount = g_oaSingleMolecules.GetSize();
963
m_bObsCertain = false;
964
m_bDecompDist = false;
965
m_bDecompType = false;
966
} else {
967
char buf[BUF_SIZE];
968
char buf2[BUF_SIZE];
969
size_t remaining = BUF_SIZE;
970
if(g_oaMolecules.GetSize() > 1) {
971
#ifdef TARGET_LINUX
972
remaining -= snprintf(buf, remaining, " Which molecule should be observed (");
973
#else
974
remaining -= sprintf(buf, " Which molecule should be observed (");
975
#endif
976
for(i = 0; i < g_oaMolecules.GetSize(); i++) {
977
if(remaining < 1)
978
break;
979
#ifdef TARGET_LINUX
980
size_t length = snprintf(buf2, remaining, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
981
#else
982
size_t length = sprintf(buf2, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
983
#endif
984
strncat(buf, buf2, remaining - 1);
985
remaining -= length;
986
if(i < g_oaMolecules.GetSize() - 1) {
987
#ifdef TARGET_LINUX
988
length = snprintf(buf2, remaining, ", ");
989
#else
990
length = sprintf(buf2, ", ");
991
#endif
992
strncat(buf, buf2, remaining - 1);
993
remaining -= length;
994
}
995
}
996
strncat(buf, ")? ", remaining - 1);
997
m_iShowMol = AskRangeInteger_ND(buf, 1, g_oaMolecules.GetSize()) - 1;
998
} else {
999
m_iShowMol = 0;
1000
}
1001
m_iShowMolCount = ((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex.GetSize();
1002
m_bObsCertain = false;
1003
m_bDecompDist = false;
1004
m_bDecompType = false;
1005
mprintf("\n");
1006
}
1007
1008
try { _ramanDyn = new CRamanDyn(m_iShowMol, global); } catch(...) { _ramanDyn = NULL; }
1009
if(_ramanDyn == NULL) NewException((double)sizeof(CRamanDyn), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1010
}
1011
1012
CRamanObservation::~CRamanObservation() {
1013
delete _ramanDyn;
1014
}
1015
1016
void CRamanObservation::initialize() {
1017
_ramanDyn->initialize();
1018
}
1019
1020
void CRamanObservation::getDipole0() {
1021
_ramanDyn->getDipole0();
1022
}
1023
1024
void CRamanObservation::calcPolarizability(int fieldDirection) {
1025
_ramanDyn->calcPolarizability(fieldDirection);
1026
}
1027
1028
void CRamanObservation::finalize() {
1029
_ramanDyn->finalize();
1030
}
1031
1032
static bool parseSettings(FILE *settingsFile) {
1033
char buf[BUF_SIZE];
1034
1035
if(fgets(buf, BUF_SIZE, settingsFile) == NULL)
1036
return false;
1037
int temp = -1;
1038
if(sscanf(buf, "%d", &temp) < 1)
1039
return false;
1040
if(temp == 0)
1041
g_orientAvg = false;
1042
else if(temp == 1)
1043
g_orientAvg = true;
1044
else
1045
return false;
1046
1047
if(fgets(buf, BUF_SIZE, settingsFile) == NULL)
1048
return false;
1049
if(sscanf(buf, "%f", &g_fieldStrength) < 1)
1050
return false;
1051
1052
if(fgets(buf, BUF_SIZE, settingsFile) == NULL)
1053
return false;
1054
if(sscanf(buf, "%d", &g_stride) < 1)
1055
return false;
1056
1057
return true;
1058
}
1059
1060
bool gatherRaman() {
1061
#ifndef TARGET_LINUX
1062
mprintf(RED, "Raman calculations are currently possible only under Linux.\n");
1063
return false;
1064
#else
1065
mprintf(" To calculate Raman spectra, polarizabilities have to be determined.\n");
1066
mprintf(" TRAVIS creates CP2K input files for numerical polarizabilities\n and will process the resulting data in a second run.\n\n");
1067
1068
g_newRaman = AskYesNo(" Do you wish to create new CP2K input files (y) or process existing results (n)? [y] ", true);
1069
char buf[BUF_SIZE];
1070
if(g_newRaman) {
1071
AskString(" Please enter a name for the directory to collect the data: [raman] ", buf, "raman");
1072
} else {
1073
AskString(" Please enter the name of the directory to take the data from: [raman] ", buf, "raman");
1074
}
1075
try { g_ramanDir = new char[strlen(buf)+1]; } catch(...) { g_ramanDir = NULL; }
1076
if(g_ramanDir == NULL) NewException((double)(strlen(buf)+1)*sizeof(char), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1077
strcpy(g_ramanDir, buf);
1078
1079
if(g_newRaman) {
1080
g_bKeepOriginalCoords = true;
1081
if(FileExist(g_ramanDir)) {
1082
mprintf(RED, "A file or a directory \"%s\" already exists. Please remove it first.\n", g_ramanDir);
1083
return false;
1084
}
1085
if(mkdir(g_ramanDir, S_IRWXU) != 0) {
1086
mprintf(RED, "Directory \"%s\" could not be created: %s\n", g_ramanDir, strerror(errno));
1087
return false;
1088
}
1089
1090
char filename[BUF_SIZE];
1091
snprintf(filename, BUF_SIZE, "%s/settings.dat", g_ramanDir);
1092
FILE *settingsFile = fopen(filename, "w");
1093
if(settingsFile == NULL) {
1094
mprintf(RED, "Could not open settings file \"%s\": %s\n", filename, strerror(errno));
1095
return false;
1096
}
1097
1098
g_orientAvg = AskYesNo("\n Use orientational averaging? [no] ", false);
1099
if(g_orientAvg)
1100
fprintf(settingsFile, "%d\n", 1);
1101
else
1102
fprintf(settingsFile, "%d\n", 0);
1103
1104
g_fieldStrength = AskFloat(" Field strength in atomic units [5.0e-4] ", 5.0e-4);
1105
fprintf(settingsFile, "%.6E\n", g_fieldStrength);
1106
1107
g_stride = AskInteger(" Calculate polarizability for every n-th timestep [1] ", 1);
1108
fprintf(settingsFile, "%d\n", g_stride);
1109
1110
fclose(settingsFile);
1111
1112
snprintf(filename, BUF_SIZE, "%s/1", g_ramanDir);
1113
if(mkdir(filename, S_IRWXU) != 0) {
1114
mprintf(RED, "Directory \"%s\" could not be created: %s\n", filename, strerror(errno));
1115
return false;
1116
}
1117
if(g_orientAvg) {
1118
snprintf(filename, BUF_SIZE, "%s/2", g_ramanDir);
1119
if(mkdir(filename, S_IRWXU) != 0) {
1120
mprintf(RED, "Directory \"%s\" could not be created: %s\n", filename, strerror(errno));
1121
return false;
1122
}
1123
snprintf(filename, BUF_SIZE, "%s/3", g_ramanDir);
1124
if(mkdir(filename, S_IRWXU) != 0) {
1125
mprintf(RED, "Directory \"%s\" could not be created: %s\n", filename, strerror(errno));
1126
return false;
1127
}
1128
}
1129
1130
FILE* templateFile;
1131
snprintf(filename, BUF_SIZE, "%s/template.inp", g_ramanDir);
1132
templateFile = fopen(filename, "w");
1133
if(templateFile == NULL) {
1134
mprintf(RED, "Could not open template file \"%s\": %s\n", filename, strerror(errno));
1135
return false;
1136
}
1137
fprintf(templateFile, "&GLOBAL\n");
1138
fprintf(templateFile, " PROJECT_NAME polarizability\n");
1139
fprintf(templateFile, " RUN_TYPE MD\n");
1140
fprintf(templateFile, " PRINT_LEVEL LOW\n");
1141
fprintf(templateFile, "&END\n");
1142
fprintf(templateFile, "&FORCE_EVAL\n");
1143
fprintf(templateFile, " &DFT\n");
1144
fprintf(templateFile, " BASIS_SET_FILE_NAME BASIS_MOLOPT\n");
1145
fprintf(templateFile, " POTENTIAL_FILE_NAME POTENTIAL\n");
1146
fprintf(templateFile, " &MGRID\n");
1147
fprintf(templateFile, " NGRIDS 5\n");
1148
fprintf(templateFile, " CUTOFF 280\n");
1149
fprintf(templateFile, " REL_CUTOFF 40\n");
1150
fprintf(templateFile, " &END\n");
1151
fprintf(templateFile, " &SCF\n");
1152
fprintf(templateFile, " SCF_GUESS ATOMIC\n");
1153
fprintf(templateFile, " MAX_SCF 200\n");
1154
fprintf(templateFile, " EPS_SCF 1.0E-5\n");
1155
fprintf(templateFile, " &OT\n");
1156
fprintf(templateFile, " MINIMIZER DIIS\n");
1157
fprintf(templateFile, " PRECONDITIONER FULL_SINGLE_INVERSE\n");
1158
fprintf(templateFile, " &END\n");
1159
fprintf(templateFile, " &END\n");
1160
fprintf(templateFile, " &XC\n");
1161
fprintf(templateFile, " &XC_FUNCTIONAL BLYP\n");
1162
fprintf(templateFile, " &END\n");
1163
fprintf(templateFile, " &XC_GRID\n");
1164
fprintf(templateFile, " XC_SMOOTH_RHO NN10\n");
1165
fprintf(templateFile, " XC_DERIV NN10_SMOOTH\n");
1166
fprintf(templateFile, " &END\n");
1167
fprintf(templateFile, " &VDW_POTENTIAL\n");
1168
fprintf(templateFile, " POTENTIAL_TYPE PAIR_POTENTIAL\n");
1169
fprintf(templateFile, " &PAIR_POTENTIAL\n");
1170
fprintf(templateFile, " TYPE DFTD3\n");
1171
fprintf(templateFile, " REFERENCE_FUNCTIONAL BLYP\n");
1172
fprintf(templateFile, " PARAMETER_FILE_NAME dftd3.dat\n");
1173
fprintf(templateFile, " &END\n");
1174
fprintf(templateFile, " &END\n");
1175
fprintf(templateFile, " &END\n");
1176
fprintf(templateFile, " &LOCALIZE\n");
1177
fprintf(templateFile, " METHOD CRAZY\n");
1178
fprintf(templateFile, " MAX_ITER 2000\n");
1179
fprintf(templateFile, " &PRINT\n");
1180
fprintf(templateFile, " &WANNIER_CENTERS\n");
1181
fprintf(templateFile, " IONS+CENTERS\n");
1182
fprintf(templateFile, " FILENAME =polarizability_wannier.xyz\n");
1183
fprintf(templateFile, " &END\n");
1184
fprintf(templateFile, " &END\n");
1185
fprintf(templateFile, " &END\n");
1186
fprintf(templateFile, " &PERIODIC_EFIELD\n");
1187
fprintf(templateFile, " INTENSITY ###!field strength will be put here###\n");
1188
fprintf(templateFile, " POLARISATION ###!polarisation vector will be put here###\n");
1189
fprintf(templateFile, " &END\n");
1190
fprintf(templateFile, " &END\n");
1191
fprintf(templateFile, " &SUBSYS\n");
1192
fprintf(templateFile, " &CELL\n");
1193
fprintf(templateFile, " ABC %.6f %.6f %.6f\n", g_fBoxX / 100.0f, g_fBoxY / 100.0f, g_fBoxZ / 100.0f);
1194
fprintf(templateFile, " &END\n");
1195
fprintf(templateFile, " &COORD\n");
1196
fprintf(templateFile, "###!coordinates will be put here###\n");
1197
fprintf(templateFile, " &END\n");
1198
fprintf(templateFile, " &KIND H\n");
1199
fprintf(templateFile, " BASIS_SET DZVP-MOLOPT-SR-GTH\n");
1200
fprintf(templateFile, " POTENTIAL GTH-BLYP-q1\n");
1201
fprintf(templateFile, " &END\n");
1202
fprintf(templateFile, " &KIND C\n");
1203
fprintf(templateFile, " BASIS_SET DZVP-MOLOPT-SR-GTH\n");
1204
fprintf(templateFile, " POTENTIAL GTH-BLYP-q4\n");
1205
fprintf(templateFile, " &END\n");
1206
fprintf(templateFile, " &KIND N\n");
1207
fprintf(templateFile, " BASIS_SET DZVP-MOLOPT-SR-GTH\n");
1208
fprintf(templateFile, " POTENTIAL GTH-BLYP-q5\n");
1209
fprintf(templateFile, " &END\n");
1210
fprintf(templateFile, " &KIND O\n");
1211
fprintf(templateFile, " BASIS_SET DZVP-MOLOPT-SR-GTH\n");
1212
fprintf(templateFile, " POTENTIAL GTH-BLYP-q6\n");
1213
fprintf(templateFile, " &END\n");
1214
fprintf(templateFile, " &END\n");
1215
fprintf(templateFile, "&END\n");
1216
fprintf(templateFile, "&MOTION\n");
1217
fprintf(templateFile, " &MD\n");
1218
fprintf(templateFile, " ENSEMBLE REFTRAJ\n");
1219
fprintf(templateFile, " STEPS ###!number of steps will be put here###\n");
1220
fprintf(templateFile, " TIMESTEP %f\n", g_fTimestepLength * g_stride);
1221
fprintf(templateFile, " &REFTRAJ\n");
1222
fprintf(templateFile, " EVAL_ENERGY_FORCES\n");
1223
fprintf(templateFile, " TRAJ_FILE_NAME ../polarizability_reftraj.xyz\n");
1224
fprintf(templateFile, " &END\n");
1225
fprintf(templateFile, " &END\n");
1226
fprintf(templateFile, " &PRINT\n");
1227
fprintf(templateFile, " &RESTART\n");
1228
fprintf(templateFile, " &EACH\n");
1229
fprintf(templateFile, " MD 1\n");
1230
fprintf(templateFile, " &END\n");
1231
fprintf(templateFile, " &END\n");
1232
fprintf(templateFile, " &END\n");
1233
fprintf(templateFile, "&END\n");
1234
fclose(templateFile);
1235
1236
mprintf("\n An input template for the polarizability calculations has been written to \"%s/template.inp\"\n Please modify it according to your needs.\n Press any key when you are finished.\n", g_ramanDir);
1237
getchar();
1238
1239
if(g_orientAvg) {
1240
mprintf(" CP2K input files will be created in \"%s/1\", \"%s/2\", and \"%s/3\".\n", g_ramanDir, g_ramanDir, g_ramanDir);
1241
mprintf(" After TRAVIS has finished, please run them and make sure that the resulting trajectories \"polarizability_wannier.xyz\"\n are placed in the same directories before you execute TRAVIS for evaluation.\n\n");
1242
} else {
1243
mprintf(" A CP2K input file will be created in \"%s/1\".\n", g_ramanDir);
1244
mprintf(" After TRAVIS has finished, please run it and make sure that the resulting trajectory \"polarizability_wannier.xyz\"\n is placed in the same directory before you run TRAVIS for evaluation.\n\n");
1245
}
1246
} else {
1247
if(!FileExist(g_ramanDir)) {
1248
mprintf(RED, "The directory \"%s\" was not found.\n", g_ramanDir);
1249
return false;
1250
}
1251
1252
char filename[BUF_SIZE];
1253
snprintf(filename, BUF_SIZE, "%s/settings.dat", g_ramanDir);
1254
FILE *settingsFile;
1255
settingsFile = fopen(filename, "r");
1256
if(settingsFile == NULL) {
1257
mprintf(RED, "Could not open settings file \"%s\": %s\n", filename, strerror(errno));
1258
return false;
1259
}
1260
if(!parseSettings(settingsFile)) {
1261
mprintf(RED, "Could not parse settings file.\n");
1262
return false;
1263
}
1264
1265
mprintf("\n");
1266
if(g_orientAvg)
1267
mprintf(" Using orientational averaging\n");
1268
else
1269
mprintf(" Not using orientational averaging\n");
1270
mprintf(" Field strength: %.6e a. u.\n", g_fieldStrength);
1271
mprintf(" Using every %d%s timestep\n\n", g_stride, (g_stride == 1) ? "st" : ((g_stride == 2) ? "nd" : ((g_stride == 3) ? "rd" : "th")));
1272
1273
g_laser = AskFloat(" Calculate scattering cross section for which laser frequency (cm^-1)? [20000.0] ", 20000.0f);
1274
g_temp = AskFloat(" Calculate scattering cross section for which temperature (K)? [300.0] ", 300.0f);
1275
1276
if(AskYesNo(" Compute Raman spectrum of whole system (y/n)? [yes] ", true)) {
1277
mprintf(YELLOW, "\n>>> Global Raman Observation >>>\n\n");
1278
1279
CRamanObservation *obs;
1280
try { obs = new CRamanObservation(true); } catch(...) { obs = NULL; }
1281
if(obs == NULL) NewException((double)sizeof(CRamanObservation), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1282
g_ramObserv.Add(obs);
1283
1284
mprintf(YELLOW, "<<< End of Global Raman Observation <<<\n\n");
1285
}
1286
1287
if(AskYesNo(" Compute Raman spectra for certain molecule types (y/n)? [no] ", false)) {
1288
while(true) {
1289
mprintf(YELLOW, "\n>>> Raman Observation %d >>>\n\n", g_ramObserv.GetSize() + 1);
1290
1291
CRamanObservation *obs;
1292
try { obs = new CRamanObservation(); } catch(...) { obs = NULL; }
1293
if(obs == NULL) NewException((double)sizeof(CRamanObservation), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1294
g_ramObserv.Add(obs);
1295
1296
mprintf(YELLOW, "<<< End of Raman Observation %d <<<\n\n", g_ramObserv.GetSize());
1297
1298
if(!AskYesNo(" Add another observation (y/n)? [no] ", false))
1299
break;
1300
mprintf("\n");
1301
}
1302
}
1303
mprintf("\n");
1304
}
1305
1306
return true;
1307
#endif
1308
}
1309
1310
bool initializeRaman() {
1311
int i;
1312
1313
if(g_newRaman) {
1314
char filename[BUF_SIZE];
1315
#ifdef TARGET_LINUX
1316
snprintf(filename, BUF_SIZE, "%s/template.inp", g_ramanDir);
1317
#else
1318
sprintf(filename, "%s/template.inp", g_ramanDir);
1319
#endif
1320
FILE *templateFile;
1321
templateFile = fopen(filename, "r");
1322
if(templateFile == NULL) {
1323
mprintf(RED, "Could not open template file \"%s\": %s\n", filename, strerror(errno));
1324
return false;
1325
}
1326
fseek(templateFile, 0L, SEEK_END);
1327
long length = ftell(templateFile);
1328
if(length < 0) {
1329
mprintf(RED, "Could not determine size of template file: %s\n", strerror(errno));
1330
return false;
1331
}
1332
1333
try { g_inputTemplate = new char[length + 1]; } catch(...) { g_inputTemplate = NULL; }
1334
if(g_inputTemplate == NULL) NewException((double)(length + 1)*sizeof(char), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1335
1336
rewind(templateFile);
1337
if((long)fread(g_inputTemplate, sizeof(char), length, templateFile) < length) {
1338
mprintf(RED, "Could not read template file: %s\n", strerror(errno));
1339
return false;
1340
}
1341
g_inputTemplate[length] = '\0';
1342
1343
fclose(templateFile);
1344
1345
g_templateFieldPos = strstr(g_inputTemplate, "###!field strength will be put here###");
1346
if(g_templateFieldPos == NULL) {
1347
mprintf(RED, "Position mark for field strength missing in template.\n");
1348
return false;
1349
}
1350
g_templatePolPos = strstr(g_inputTemplate, "###!polarisation vector will be put here###");
1351
if(g_templatePolPos == NULL) {
1352
mprintf(RED, "Position mark for polarisation vector missing in template.\n");
1353
return false;
1354
}
1355
g_templateCoordPos = strstr(g_inputTemplate, "###!coordinates will be put here###");
1356
if(g_templateCoordPos == NULL) {
1357
mprintf(RED, "Position mark for coordinates missing in template.\n");
1358
return false;
1359
}
1360
g_templateStepsPos = strstr(g_inputTemplate, "###!number of steps will be put here###");
1361
if(g_templateStepsPos == NULL) {
1362
mprintf(RED, "Position mark for number of steps missing in template.\n");
1363
return false;
1364
}
1365
g_templateFieldPos[0] = '\0';
1366
g_templatePolPos[0] = '\0';
1367
g_templateCoordPos[0] = '\0';
1368
g_templateStepsPos[0] = '\0';
1369
1370
#ifdef TARGET_LINUX
1371
snprintf(filename, BUF_SIZE, "%s/polarizability_reftraj.xyz", g_ramanDir);
1372
#else
1373
sprintf(filename, "%s/polarizability_reftraj.xyz", g_ramanDir);
1374
#endif
1375
g_reftrajFile = fopen(filename, "w");
1376
if(g_reftrajFile == NULL) {
1377
mprintf(RED, "Could not open reference trajectory \"%s\": %s\n", filename, strerror(errno));
1378
return false;
1379
}
1380
} else {
1381
for(i = 0; i < g_ramObserv.GetSize(); i++) {
1382
mprintf("Initializing Raman Observation %d...\n", i+1);
1383
CRamanObservation *obs = (CRamanObservation *)g_ramObserv[i];
1384
obs->initialize();
1385
}
1386
1387
char filename[BUF_SIZE];
1388
for(i = 0; i < (g_orientAvg ? 3 : 1); i++) {
1389
#ifdef TARGET_LINUX
1390
snprintf(filename, BUF_SIZE, "%s/%d/polarizability_wannier.xyz", g_ramanDir, i + 1);
1391
#else
1392
sprintf(filename, "%s/%d/polarizability_wannier.xyz", g_ramanDir, i + 1);
1393
#endif
1394
g_polFile[i] = fopen(filename, "r");
1395
if(g_polFile[i] == NULL) {
1396
mprintf(RED, "Could not open trajectory \"%s\": %s\n", filename, strerror(errno));
1397
return false;
1398
}
1399
}
1400
1401
for(i = 0; i < (g_orientAvg ? 3 : 1); i++) {
1402
try { g_timestep[i] = new CTimeStep(); } catch(...) { g_timestep[i] = NULL; }
1403
if(g_timestep[i] == NULL) NewException((double)sizeof(CTimeStep), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1404
}
1405
}
1406
1407
return true;
1408
}
1409
1410
void processRaman(CTimeStep *ts) {
1411
int i, j;
1412
g_step++;
1413
if(g_newRaman) {
1414
if(g_step % g_stride == 0) {
1415
g_steps++;
1416
int numAtoms = 0;
1417
for(i = 0; i < g_iGesAtomCount; i++)
1418
if(g_waAtomMolIndex[i] != 60000)
1419
numAtoms++;
1420
fprintf(g_reftrajFile, "%d\nStep %d\n", numAtoms, g_step);
1421
for(i = 0; i < g_iGesAtomCount; i++) {
1422
if(g_waAtomMolIndex[i] == 60000)
1423
continue;
1424
fprintf(g_reftrajFile, "%4s %14.10f %14.10f %14.10f\n", ((CAtom *)g_oaAtoms[g_waAtomRealElement[i]])->m_sName, ts->m_vaCoords_Original[i][0] / 100.0f, ts->m_vaCoords_Original[i][1] / 100.0f, ts->m_vaCoords_Original[i][2] / 100.0f);
1425
}
1426
}
1427
} else {
1428
if(g_step % g_stride == 0) {
1429
for(i = 0; i < g_ramObserv.GetSize(); i++)
1430
((CRamanObservation *)g_ramObserv[i])->getDipole0();
1431
for(i = 0; i < (g_orientAvg ? 3 : 1); i++) {
1432
if(!g_timestep[i]->ReadTimestep(g_polFile[i], false)) {
1433
mprintf(RED, "Error while reading trajectory for polarizabilities.\n");
1434
abort();
1435
}
1436
if(!g_bSaveCoordsUnchanged) {
1437
g_timestep[i]->UniteMolecules(false);
1438
if(g_bRemoveCOM)
1439
g_timestep[i]->CenterCOM();
1440
}
1441
g_timestep[i]->CalcCenters();
1442
if(g_bWannier)
1443
g_timestep[i]->ScanWannier(false);
1444
g_timestep[i]->CalcDipoles();
1445
for(j = 0; j < g_ramObserv.GetSize(); j++) {
1446
((CRamanObservation *)g_ramObserv[j])->calcPolarizability(i);
1447
}
1448
}
1449
}
1450
}
1451
}
1452
1453
void finalizeRaman() {
1454
int i, j, k, l;
1455
if(g_newRaman) {
1456
fclose(g_reftrajFile);
1457
1458
char filename[BUF_SIZE];
1459
FILE *inputFile;
1460
for(i = 1; i < (g_orientAvg ? 4 : 2); i++) {
1461
#ifdef TARGET_LINUX
1462
snprintf(filename, BUF_SIZE, "%s/%d/polarizability.inp", g_ramanDir, i);
1463
#else
1464
sprintf(filename, "%s/%d/polarizability.inp", g_ramanDir, i);
1465
#endif
1466
inputFile = fopen(filename, "w");
1467
if(inputFile == NULL) {
1468
mprintf(RED, "Could not open input file \"%s\": %s\n", filename, strerror(errno));
1469
abort();
1470
}
1471
1472
fprintf(inputFile, "%s", g_inputTemplate);
1473
for(j = 0; j < 4; j++) {
1474
char *tempPos = g_inputTemplate;
1475
for(k = 0; k <= j; k++)
1476
tempPos = strchr(&tempPos[1], '\0');
1477
if(tempPos == g_templateFieldPos) {
1478
fprintf(inputFile, "%.6E", g_fieldStrength);
1479
fprintf(inputFile, "%s", &g_templateFieldPos[38]);
1480
} else if(tempPos == g_templatePolPos) {
1481
if(i == 1)
1482
fprintf(inputFile, "1.0 0.0 0.0");
1483
if(i == 2)
1484
fprintf(inputFile, "0.0 1.0 0.0");
1485
if(i == 3)
1486
fprintf(inputFile, "0.0 0.0 1.0");
1487
fprintf(inputFile, "%s", &g_templatePolPos[43]);
1488
} else if(tempPos == g_templateCoordPos) {
1489
CTimeStep *ts = GetTimeStep(0);
1490
for(l = 0; l < g_iGesAtomCount; l++) {
1491
if(g_waAtomMolIndex[l] == 60000)
1492
continue;
1493
fprintf(inputFile, " %s %14.10f %14.10f %14.10f\n", ((CAtom *)g_oaAtoms[g_waAtomRealElement[l]])->m_sName, ts->m_vaCoords_Original[l][0] / 100.0f, ts->m_vaCoords_Original[l][1] / 100.0f, ts->m_vaCoords_Original[l][2] / 100.0f);
1494
}
1495
fprintf(inputFile, "%s", &g_templateCoordPos[36]);
1496
} else if(tempPos == g_templateStepsPos) {
1497
fprintf(inputFile, "%d", g_steps);
1498
fprintf(inputFile, "%s", &g_templateStepsPos[39]);
1499
} else {
1500
mprintf(RED, "Unexpected error while processing the input template.\n");
1501
}
1502
}
1503
fclose(inputFile);
1504
}
1505
delete[] g_inputTemplate;
1506
} else {
1507
for(i = 0; i < g_ramObserv.GetSize(); i++) {
1508
mprintf(YELLOW, ">>> Raman Observation %d >>>\n\n", i + 1);
1509
((CRamanObservation *)g_ramObserv[i])->finalize();
1510
mprintf(YELLOW, "\n<<< End of Raman Observation %d <<<\n\n", i + 1);
1511
}
1512
1513
for(i = 0; i < (g_orientAvg ? 3 : 1); i++)
1514
fclose(g_polFile[i]);
1515
for(i = 0; i < (g_orientAvg ? 3 : 1); i++)
1516
delete g_timestep[i];
1517
}
1518
delete[] g_ramanDir;
1519
}
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 "raman.h"
25
26
#include "acf.h"
27
#include "df.h"
28
#include "globalvar.h"
29
#include "maintools.h"
30
#include "moltools.h"
31
#include "timestep.h"
32
#include "tools.h"
33
#include "xfloatarray.h"
34
#include "xobarray.h"
35
36
#include <errno.h>
37
#include <stdio.h>
38
#include <string.h>
39
40
#ifdef TARGET_LINUX
41
#include <sys/stat.h>
42
#include <sys/types.h>
43
#endif
44
45
#define BUF_SIZE 4096
46
47
static bool g_newRaman;
48
static char *g_ramanDir;
49
static bool g_orientAvg;
50
static float g_fieldStrength;
51
static int g_step = 0;
52
static int g_stride;
53
static float g_laser;
54
static float g_temp;
55
static char *g_inputTemplate;
56
static char *g_templateFieldPos;
57
static char *g_templatePolPos;
58
static char *g_templateCoordPos;
59
static char *g_templateStepsPos;
60
61
static CxObArray g_ramObserv;
62
63
static FILE *g_polFile[3];
64
static CTimeStep *g_timestep[3];
65
66
static FILE *g_reftrajFile;
67
static int g_steps = 0;
68
69
CRamanDyn::CRamanDyn(int showMol, bool global) {
70
_global = global;
71
if(_global) {
72
m_iShowMol = -1;
73
m_iMolecules = g_oaSingleMolecules.GetSize();
74
} else {
75
m_iShowMol = showMol;
76
m_iMolecules = ((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex.GetSize();
77
}
78
79
/*** Set Names of dynamic arrays. By M. Brehm */
80
_dipole0.SetName("CRamanDyn::_dipole0");
81
char tbuf[256];
82
for (int tz1=0;tz1<3;tz1++)
83
{
84
for (int tz2=0;tz2<3;tz2++)
85
{
86
sprintf(tbuf,"CRamanDyn::_polarizability[%d][%d]",tz1,tz2);
87
_polarizability[tz1][tz2].SetName(tbuf);
88
}
89
}
90
/* End of set names */
91
92
mprintf(YELLOW, ">>> Raman Spectrum >>>\n\n");
93
94
if(_global)
95
mprintf(" All atoms will be taken.\n\n");
96
else
97
mprintf(" All atoms will be taken from the OM %s.\n\n", ((CMolecule *)g_oaMolecules[m_iShowMol])->m_sName);
98
99
m_iVecType = 1;
100
m_iCombinations = 1;
101
g_bDipole = true;
102
ParseDipole();
103
104
if(g_iTrajSteps != -1) {
105
int depth = (int)(g_iTrajSteps / g_stride * 0.75);
106
if(depth > 4096)
107
depth = 4096;
108
m_iDepth = AskUnsignedInteger(" Enter the resolution (=depth) of the ACF (in time steps): [%d] ", depth, depth);
109
} else {
110
m_iDepth = AskUnsignedInteger(" Enter the resolution (=depth) of the ACF (in time steps): [2000] ", 2000);
111
}
112
113
int size = CalcFFTSize(m_iDepth, false);
114
if(m_iDepth != size) {
115
mprintf(WHITE," The next \"fast\" size for FFT is %d. Using this instead of %d as depth.\n", size, m_iDepth);
116
m_iDepth = size;
117
}
118
119
m_iStride = 1;
120
m_bSpectrum = true;
121
122
_derivativeOrder = 1;
123
bool window = true;
124
bool derive = true;
125
if(g_bAdvanced2) {
126
derive = AskYesNo(" Derive the vectors before autocorrelation (y/n)? [yes] ", true);
127
if(derive)
128
_derivativeOrder = AskRangeInteger(" Please enter degree of vector derivative (1-2): [1] ", 1, 2, 1);
129
window = AskYesNo(" Apply window function (Cos^2) to autocorrelation function (y/n)? [yes] ", true);
130
}
131
float possibleRange = 33356.41f / g_fTimestepLength / 2.0f / g_stride;
132
mprintf(" A time step length of %.2f fs with a stride of %d allows a spectral range up to %.2f cm^-1.\n", g_fTimestepLength, g_stride, possibleRange);
133
float specWaveNumber = AskRangeFloat("\n Calculate spectrum up to which wave number (cm^-1)? [%.2f cm^-1] ", 0, possibleRange, (possibleRange < 5000.0f) ? possibleRange : 5000.0f, (possibleRange < 5000.0f) ? possibleRange : 5000.0f);
134
int mirror = 1;
135
int zeroPadding = m_iDepth * 3;
136
if(g_bAdvanced2) {
137
mirror = AskRangeInteger(" No mirroring (0) or short-time enhancing (1)? [1] ", 0, 1, 1);
138
zeroPadding = AskUnsignedInteger(" Zero Padding: How many zeros to insert? [%d] ", m_iDepth * 3, m_iDepth * 3);
139
}
140
141
size = CalcFFTSize(m_iDepth + zeroPadding, false);
142
if(m_iDepth + zeroPadding != size) {
143
mprintf(WHITE, " The next \"fast\" size for FFT is %d. Using %d zeros for zero padding.\n", size, size-m_iDepth);
144
zeroPadding = size-m_iDepth;
145
}
146
147
mprintf(" This results in a spectral resolution to %.2f cm^-1.\n\n", 33356.41 / g_fTimestepLength / 2.0f / size);
148
149
try { _isoACF = new CACF(); } catch(...) { _isoACF = NULL; }
150
if(_isoACF == NULL) NewException((double)sizeof(CACF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
151
_isoACF->m_iSize = m_iDepth;
152
_isoACF->m_bSpectrum = true;
153
_isoACF->m_bDerivative = derive;
154
_isoACF->m_iDerivative = _derivativeOrder;
155
_isoACF->m_bWindowFunction = window;
156
_isoACF->m_fSpecWaveNumber = specWaveNumber;
157
_isoACF->m_iMirror = mirror;
158
_isoACF->m_iZeroPadding = zeroPadding;
159
_isoACF->m_bACF_DB = false;
160
161
try { _anisoACF = new CACF(); } catch(...) { _anisoACF = NULL; }
162
if(_anisoACF == NULL) NewException((double)sizeof(CACF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
163
_anisoACF->m_iSize = m_iDepth;
164
_anisoACF->m_bSpectrum = true;
165
_anisoACF->m_bDerivative = derive;
166
_anisoACF->m_iDerivative = _derivativeOrder;
167
_anisoACF->m_bWindowFunction = window;
168
_anisoACF->m_fSpecWaveNumber = specWaveNumber;
169
_anisoACF->m_iMirror = mirror;
170
_anisoACF->m_iZeroPadding = zeroPadding;
171
_anisoACF->m_bACF_DB = false;
172
173
BuildName();
174
175
mprintf(YELLOW, "<<< End of Raman Spectrum <<<\n\n");
176
}
177
178
CRamanDyn::~CRamanDyn() {
179
int i, j, k;
180
181
for(i = 0; i < m_iMolecules; i++)
182
delete (CxFloatArray *)_dipole0[i];
183
for(i = 0; i < m_iMolecules; i++) {
184
for(j = 0; j < 3; j++) {
185
for(k = 0; k < (g_orientAvg ? 3 : 1); k++) {
186
delete (CxFloatArray *)_polarizability[j][k][i];
187
}
188
}
189
}
190
delete _isoACF;
191
delete _anisoACF;
192
}
193
194
void CRamanDyn::initialize() {
195
int i, j, k;
196
197
_isoACF->Create();
198
_anisoACF->Create();
199
200
if (g_iTrajSteps != -1)
201
mprintf(" Raman Cache: Trying to allocate %s of memory...\n", FormatBytes((double)m_iMolecules * g_iTrajSteps / g_iStride / g_stride * (g_orientAvg ? 9.9 : 3.3) * sizeof(float)));
202
else
203
mprintf(" Raman Cache: Trying to allocate %s of memory...\n", FormatBytes((double)m_iMolecules * 2000 / g_iStride /g_stride * (g_orientAvg ? 9.9 : 3.3) * sizeof(float)));
204
for(i = 0; i < m_iMolecules; i++) {
205
CxVector3 *dipoleVector;
206
try { dipoleVector = new CxVector3(); } catch(...) { dipoleVector = NULL; }
207
if(dipoleVector == NULL) NewException((double)sizeof(CxVector3), __FILE__, __LINE__, __PRETTY_FUNCTION__);
208
_dipole0.Add(dipoleVector);
209
}
210
for(i = 0; i < m_iMolecules; i++) {
211
for(j = 0; j < 3; j++) {
212
for(k = 0; k < (g_orientAvg ? 3 : 1); k++) {
213
CxFloatArray *floatArray;
214
try { floatArray = new CxFloatArray("CRamanDyn::initialize():floatArray"); } catch(...) { floatArray = NULL; }
215
if(floatArray == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
216
if(g_iTrajSteps != -1) {
217
floatArray->SetMaxSize((long)(g_iTrajSteps / g_iStride / g_stride * 1.1));
218
floatArray->SetGrow((long)(g_iTrajSteps / g_iStride / g_stride * 0.1));
219
} else {
220
floatArray->SetGrow(10000);
221
}
222
_polarizability[j][k].Add(floatArray);
223
}
224
}
225
}
226
}
227
228
void CRamanDyn::getDipole0() {
229
int i;
230
231
for(i = 0; i < m_iMolecules; i++) {
232
CSingleMolecule *sm;
233
if(_global)
234
sm = (CSingleMolecule *)g_oaSingleMolecules[i];
235
else
236
sm = (CSingleMolecule *)g_oaSingleMolecules[((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex[i]];
237
*((CxVector3 *)_dipole0[i]) = sm->m_vDipole;
238
}
239
}
240
241
void CRamanDyn::calcPolarizability(int fieldDirection) {
242
int i;
243
244
for(i = 0; i < m_iMolecules; i++) {
245
CSingleMolecule *sm;
246
if(_global)
247
sm = (CSingleMolecule *)g_oaSingleMolecules[i];
248
else
249
sm = (CSingleMolecule *)g_oaSingleMolecules[((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex[i]];
250
for(int j = 0; j < 3; j++) {
251
((CxFloatArray *)_polarizability[j][fieldDirection][i])->Add(((*((CxVector3 *)_dipole0[i]))[j] - sm->m_vDipole[j]) / g_fieldStrength * 0.393430f); // 0.393430 - Umrechnung von Debye in a.u.
252
}
253
}
254
}
255
256
void CRamanDyn::finalize() {
257
int i, j, k, l;
258
259
char filename[BUF_SIZE];
260
#ifdef TARGET_LINUX
261
snprintf(filename, BUF_SIZE, "%s/polarizability_%s.dat", g_ramanDir, m_sName);
262
#else
263
sprintf(filename, "%s/polarizability_%s.dat", g_ramanDir, m_sName);
264
#endif
265
mprintf(" Writing polarizabilities for first molecule to \"%s\"...\n", filename);
266
FILE *pol_file;
267
pol_file = OpenFileWrite(filename, false);
268
for(i = 0; i < ((CxFloatArray *)_polarizability[0][0][0])->GetSize(); i++) {
269
fprintf(pol_file, "%10.2f", i * g_fTimestepLength * g_iStride * g_stride);
270
for(j = 0; j < (g_orientAvg ? 3 : 1); j++) {
271
for(k = 0; k < 3; k++) {
272
fprintf(pol_file, " %14.8f", (*((CxFloatArray *)_polarizability[k][j][0]))[i]);
273
}
274
}
275
fprintf(pol_file, "\n");
276
}
277
fclose(pol_file);
278
279
float step;
280
switch(_derivativeOrder) {
281
case 0:
282
mprintf(" Not deriving polarizabilities.\n");
283
break;
284
case 1:
285
mprintf(" Deriving polarizabilities (1st derivative)...\n");
286
mprintf(WHITE, " [");
287
step = m_iMolecules / 20.0f;
288
for(i = 0; i < m_iMolecules; i++) {
289
if(fmodf((float)i, step) < 1.0f)
290
mprintf(WHITE, "#");
291
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - 1; j++) {
292
for(k = 0; k < 3; k++) {
293
for(l = 0; l < (g_orientAvg ? 3 : 1); l++) {
294
(*((CxFloatArray *)_polarizability[k][l][i]))[j] = 0.5f * ((*((CxFloatArray *)_polarizability[k][l][i]))[j+1] - (*((CxFloatArray *)_polarizability[k][l][i]))[j]);
295
}
296
}
297
}
298
}
299
mprintf(WHITE, "]\n");
300
break;
301
case 2:
302
mprintf(" Deriving polarizabilities (2nd derivative)...\n");
303
mprintf(WHITE, " [");
304
step = m_iMolecules / 20.0f;
305
for(i = 0; i < m_iMolecules; i++) {
306
if(fmodf((float)i, step) < 1.0f)
307
mprintf(WHITE, "#");
308
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - 2; j++) {
309
for(k = 0; k < 3; k++) {
310
for(l = 0; l < (g_orientAvg ? 3 : 1); l++) {
311
(*((CxFloatArray *)_polarizability[k][l][i]))[j] = (*((CxFloatArray *)_polarizability[k][l][i]))[j+2] - 2.0f * (*((CxFloatArray *)_polarizability[k][l][i]))[j+1] + (*((CxFloatArray *)_polarizability[k][l][i]))[j];
312
}
313
}
314
}
315
}
316
mprintf(WHITE, "]\n");
317
break;
318
default:
319
mprintf(RED, "Higher derivatives not implemented.\n");
320
abort();
321
break;
322
}
323
324
mprintf(" Processing polarizability tensor components...\n");
325
step = m_iMolecules / 20.0f;
326
327
CAutoCorrelation *autoCorrelation;
328
try { autoCorrelation = new CAutoCorrelation(); } catch(...) { autoCorrelation = NULL; }
329
if(autoCorrelation == NULL) NewException((double)sizeof(CAutoCorrelation), __FILE__, __LINE__, __PRETTY_FUNCTION__);
330
autoCorrelation->Init(((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder, m_iDepth, g_bACFFFT);
331
332
if(g_orientAvg) {
333
CxFloatArray *isotropyPol, *isotropyACF;
334
try { isotropyACF = new CxFloatArray("CRamanDyn::finalize():isotropyACF"); } catch(...) { isotropyACF = NULL; }
335
if(isotropyACF == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
336
isotropyACF->SetSize(m_iDepth);
337
try { isotropyPol = new CxFloatArray("CRamanDyn::finalize():isotropyPol"); } catch(...) { isotropyPol = NULL; }
338
if(isotropyPol == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
339
isotropyPol->SetSize(((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder);
340
341
mprintf(" Isotropic part...\n");
342
mprintf(WHITE, " [");
343
for(i = 0; i < m_iDepth; i++) {
344
_isoACF->m_pData[i] = 0.0f;
345
}
346
for(i = 0; i < m_iMolecules; i++) {
347
if(fmodf((float)i, step) < 1.0f)
348
mprintf(WHITE, "#");
349
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
350
(*isotropyPol)[j] = 0.0f;
351
for(k = 0; k < 3; k++) {
352
(*isotropyPol)[j] += (*((CxFloatArray *)_polarizability[k][k][i]))[j];
353
}
354
(*isotropyPol)[j] /= 3.0f;
355
}
356
autoCorrelation->AutoCorrelate(isotropyPol, isotropyACF);
357
for(j = 0; j < m_iDepth; j++) {
358
_isoACF->m_pData[j] += (*isotropyACF)[j];
359
}
360
}
361
for(i = 0; i < m_iDepth; i++) {
362
if(_global)
363
_isoACF->m_pData[i] /= g_iSteps / g_stride;
364
else
365
_isoACF->m_pData[i] /= g_iSteps / g_stride * m_iMolecules;
366
}
367
mprintf(WHITE, "]\n");
368
369
// #ifdef TARGET_LINUX
370
// snprintf(filename, BUF_SIZE, "%s/acf_iso_%s.dat", g_ramanDir, m_sName);
371
// #else
372
// sprintf(filename, "%s/acf_iso_%s.dat", g_ramanDir, m_sName);
373
// #endif
374
// mprintf(" Saving ACF as %s...\n", filename);
375
// FILE *acf_file = OpenFileWrite(filename, false);
376
// for(i = 0; i < m_iDepth; i++)
377
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _isoACF->m_pData[i]);
378
// fclose(acf_file);
379
380
if(_isoACF->m_iMirror != 0) {
381
mprintf(" Mirroring ACF...\n");
382
_isoACF->Mirror(_isoACF->m_iMirror);
383
}
384
if(_isoACF->m_bWindowFunction != 0) {
385
mprintf(" Applying window function to ACF...\n");
386
_isoACF->Window();
387
}
388
389
// #ifdef TARGET_LINUX
390
// snprintf(filename, BUF_SIZE, "%s/acf_iso_%s.mw.dat", g_ramanDir, m_sName);
391
// #else
392
// sprintf(filename, "%s/acf_iso_%s.mw.dat", g_ramanDir, m_sName);
393
// #endif
394
// mprintf(" Saving ACF as %s...\n", filename);
395
// acf_file = OpenFileWrite(filename, false);
396
// for(i = 0; i < _isoACF->m_iSize; i++)
397
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _isoACF->m_pData[i]);
398
// fclose(acf_file);
399
400
mprintf(" Performing Fourier transformation...\n");
401
CFFT *fft;
402
try { fft = new CFFT(); } catch(...) { fft = NULL; }
403
if(fft == NULL) NewException((double)sizeof(CFFT), __FILE__, __LINE__, __PRETTY_FUNCTION__);
404
fft->PrepareFFT_C2C(_isoACF->m_iSize + _isoACF->m_iZeroPadding);
405
_isoACF->Transform(fft);
406
delete fft;
407
_isoACF->m_pSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
408
409
// #ifdef TARGET_LINUX
410
// snprintf(filename, BUF_SIZE, "%s/spectrum_iso_%s.dat", g_ramanDir, m_sName);
411
// #else
412
// sprintf(filename, "%s/spectrum_iso_%s.dat", g_ramanDir, m_sName);
413
// #endif
414
// mprintf(" Saving spectrum as %s...\n", filename);
415
// _isoACF->m_pSpectrum->Write("", filename, "");
416
417
delete isotropyACF;
418
delete isotropyPol;
419
420
CxFloatArray *anisotropyPol, *anisotropyACF, *anisotropyACFSum;
421
try { anisotropyPol = new CxFloatArray("CRamanDyn::finalize():anisotropyPol"); } catch(...) { anisotropyPol = NULL; }
422
if(anisotropyPol == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
423
anisotropyPol->SetSize(((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder);
424
try { anisotropyACF = new CxFloatArray("CRamanDyn::finalize():anisotropyACF"); } catch(...) { anisotropyACF = NULL; }
425
if(anisotropyACF == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
426
anisotropyACF->SetSize(m_iDepth);
427
try { anisotropyACFSum = new CxFloatArray("CRamanDyn::finalize():anisotropyACFSum"); } catch(...) { anisotropyACFSum = NULL; }
428
if(anisotropyACFSum == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
429
anisotropyACFSum->SetSize(m_iDepth);
430
431
mprintf(" First anisotropic part...\n");
432
mprintf(WHITE, " [");
433
for(i = 0; i < m_iDepth; i++) {
434
_anisoACF->m_pData[i] = 0.0f;
435
}
436
for(i = 0; i < m_iMolecules; i++) {
437
if(fmodf((float)i, step) < 1.0f)
438
mprintf(WHITE, "#");
439
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
440
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[0][0][i]))[j];
441
(*anisotropyPol)[j] -= (*((CxFloatArray *)_polarizability[1][1][i]))[j];
442
}
443
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
444
for(int j = 0; j < m_iDepth; j++) {
445
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
446
}
447
}
448
for(i = 0; i < m_iDepth; i++) {
449
(*anisotropyACFSum)[i] = 0.5f * _anisoACF->m_pData[i];
450
}
451
mprintf(WHITE, "]\n");
452
453
mprintf(" Second anisotropic part...\n");
454
mprintf(WHITE, " [");
455
for(i = 0; i < m_iDepth; i++) {
456
_anisoACF->m_pData[i] = 0.0f;
457
}
458
for(i = 0; i < m_iMolecules; i++) {
459
if(fmodf((float)i, step) < 1.0f)
460
mprintf(WHITE, "#");
461
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
462
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[1][1][i]))[j];
463
(*anisotropyPol)[j] -= (*((CxFloatArray *)_polarizability[2][2][i]))[j];
464
}
465
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
466
for(j = 0; j < m_iDepth; j++) {
467
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
468
}
469
}
470
for(i = 0; i < m_iDepth; i++) {
471
(*anisotropyACFSum)[i] += 0.5f * _anisoACF->m_pData[i];
472
}
473
mprintf(WHITE, "]\n");
474
475
mprintf(" Third anisotropic part...\n");
476
mprintf(WHITE, " [");
477
for(i = 0; i < m_iDepth; i++) {
478
_anisoACF->m_pData[i] = 0.0f;
479
}
480
for(i = 0; i < m_iMolecules; i++) {
481
if(fmodf((float)i, step) < 1.0f)
482
mprintf(WHITE, "#");
483
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
484
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[2][2][i]))[j];
485
(*anisotropyPol)[j] -= (*((CxFloatArray *)_polarizability[0][0][i]))[j];
486
}
487
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
488
for(j = 0; j < m_iDepth; j++) {
489
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
490
}
491
}
492
for(i = 0; i < m_iDepth; i++) {
493
(*anisotropyACFSum)[i] += 0.5f * _anisoACF->m_pData[i];
494
}
495
mprintf(WHITE, "]\n");
496
497
mprintf(" Fourth anisotropic part...\n");
498
mprintf(WHITE, " [");
499
for(i = 0; i < m_iDepth; i++) {
500
_anisoACF->m_pData[i] = 0.0f;
501
}
502
for(i = 0; i < m_iMolecules; i++) {
503
if(fmodf((float)i, step) < 1.0f)
504
mprintf(WHITE, "#");
505
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
506
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[0][1][i]))[j];
507
(*anisotropyPol)[j] += (*((CxFloatArray *)_polarizability[1][0][i]))[j];
508
(*anisotropyPol)[j] *= 0.5f;
509
}
510
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
511
for(j = 0; j < m_iDepth; j++) {
512
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
513
}
514
}
515
for(i = 0; i < m_iDepth; i++) {
516
(*anisotropyACFSum)[i] += 3.0f * _anisoACF->m_pData[i];
517
}
518
mprintf(WHITE, "]\n");
519
520
mprintf(" Fifth anisotropic part...\n");
521
mprintf(WHITE, " [");
522
for(i = 0; i < m_iDepth; i++) {
523
_anisoACF->m_pData[i] = 0.0f;
524
}
525
for(i = 0; i < m_iMolecules; i++) {
526
if(fmodf((float)i, step) < 1.0f)
527
mprintf(WHITE, "#");
528
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
529
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[1][2][i]))[j];
530
(*anisotropyPol)[j] += (*((CxFloatArray *)_polarizability[2][1][i]))[j];
531
(*anisotropyPol)[j] *= 0.5f;
532
}
533
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
534
for(j = 0; j < m_iDepth; j++) {
535
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
536
}
537
}
538
for(i = 0; i < m_iDepth; i++) {
539
(*anisotropyACFSum)[i] += 3.0f * _anisoACF->m_pData[i];
540
}
541
mprintf(WHITE, "]\n");
542
543
mprintf(" Sixth anisotropic part...\n");
544
mprintf(WHITE, " [");
545
for(i = 0; i < m_iDepth; i++) {
546
_anisoACF->m_pData[i] = 0.0f;
547
}
548
for(i = 0; i < m_iMolecules; i++) {
549
if(fmodf((float)i, step) < 1.0f)
550
mprintf(WHITE, "#");
551
for(j = 0; j < ((CxFloatArray *)_polarizability[0][0][0])->GetSize() - _derivativeOrder; j++) {
552
(*anisotropyPol)[j] = (*((CxFloatArray *)_polarizability[2][0][i]))[j];
553
(*anisotropyPol)[j] += (*((CxFloatArray *)_polarizability[0][2][i]))[j];
554
(*anisotropyPol)[j] *= 0.5f;
555
}
556
autoCorrelation->AutoCorrelate(anisotropyPol, anisotropyACF);
557
for(j = 0; j < m_iDepth; j++) {
558
_anisoACF->m_pData[j] += (*anisotropyACF)[j];
559
}
560
}
561
for(i = 0; i < m_iDepth; i++) {
562
(*anisotropyACFSum)[i] += 3.0f * _anisoACF->m_pData[i];
563
}
564
mprintf(WHITE, "]\n");
565
566
for(i = 0; i < m_iDepth; i++) {
567
if(_global)
568
_anisoACF->m_pData[i] = (*anisotropyACFSum)[i] / g_iSteps * g_stride;
569
else
570
_anisoACF->m_pData[i] = (*anisotropyACFSum)[i] / g_iSteps * g_stride / m_iMolecules;
571
}
572
573
// #ifdef TARGET_LINUX
574
// snprintf(filename, BUF_SIZE, "%s/acf_aniso_%s.dat", g_ramanDir, m_sName);
575
// #else
576
// sprintf(filename, "%s/acf_aniso_%s.dat", g_ramanDirm_sName);
577
// #endif
578
// mprintf(" Saving ACF as %s...\n", filename);
579
// acf_file = OpenFileWrite(filename, false);
580
// for(i = 0; i < m_iDepth; i++)
581
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _anisoACF->m_pData[i]);
582
// fclose(acf_file);
583
584
if(_anisoACF->m_iMirror != 0) {
585
mprintf(" Mirroring ACF...\n");
586
_anisoACF->Mirror(_anisoACF->m_iMirror);
587
}
588
if(_anisoACF->m_bWindowFunction != 0) {
589
mprintf(" Applying window function to ACF...\n");
590
_anisoACF->Window();
591
}
592
593
// #ifdef TARGET_LINUX
594
// snprintf(filename, BUF_SIZE, "%s/acf_aniso_%s.mw.dat", g_ramanDir, m_sName);
595
// #else
596
// sprintf(filename, "%s/acf_aniso_%s.mw.dat", g_ramanDir, m_sName);
597
// #endif
598
// mprintf(" Saving ACF as %s...\n", filename);
599
// acf_file = OpenFileWrite(filename, false);
600
// for(i = 0; i < _anisoACF->m_iSize; i++)
601
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _anisoACF->m_pData[i]);
602
// fclose(acf_file);
603
604
mprintf(" Performing Fourier transformation...\n");
605
try { fft = new CFFT(); } catch(...) { fft = NULL; }
606
if(fft == NULL) NewException((double)sizeof(CFFT), __FILE__, __LINE__, __PRETTY_FUNCTION__);
607
fft->PrepareFFT_C2C(_anisoACF->m_iSize + _anisoACF->m_iZeroPadding);
608
_anisoACF->Transform(fft);
609
delete fft;
610
_anisoACF->m_pSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
611
612
// #ifdef TARGET_LINUX
613
// snprintf(filename, BUF_SIZE, "%s/spectrum_aniso_%s.dat", g_ramanDir, m_sName);
614
// #else
615
// sprintf(filename, "%s/spectrum_aniso_%s.dat", g_ramanDir, m_sName);
616
// #endif
617
// mprintf(" Saving spectrum as %s...\n", filename);
618
// _anisoACF->m_pSpectrum->Write("", filename, "");
619
620
int specSize = _isoACF->m_pSpectrum->m_iSize;
621
CSpectrum *paraSpectrum, *orthoSpectrum, *sumSpectrum, *depolSpectrum;
622
try { paraSpectrum = new CSpectrum(); } catch(...) { paraSpectrum = NULL; }
623
if(paraSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
624
paraSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
625
paraSpectrum->Create(specSize);
626
paraSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
627
try { orthoSpectrum = new CSpectrum(); } catch(...) { orthoSpectrum = NULL; }
628
if(orthoSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
629
orthoSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
630
orthoSpectrum->Create(specSize);
631
orthoSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
632
try { sumSpectrum = new CSpectrum(); } catch(...) { sumSpectrum = NULL; }
633
if(sumSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
634
sumSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
635
sumSpectrum->Create(specSize);
636
sumSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
637
try { depolSpectrum = new CSpectrum(); } catch(...) { depolSpectrum = NULL; }
638
if(depolSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
639
depolSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
640
depolSpectrum->Create(specSize);
641
depolSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
642
643
for(i = 0; i < specSize; i++) {
644
paraSpectrum->m_pData[i] = _isoACF->m_pSpectrum->m_pData[i] + 4.0f / 45.0f * _anisoACF->m_pSpectrum->m_pData[i];
645
}
646
647
// #ifdef TARGET_LINUX
648
// snprintf(filename, BUF_SIZE, "%s/spectrum_para_%s.dat", g_ramanDir, m_sName);
649
// #else
650
// sprintf(filename, "%s/spectrum_para_%s.dat", g_ramanDir, m_sName);
651
// #endif
652
// mprintf(" Saving para spectrum as %s...\n", filename);
653
// paraSpectrum->Write("", filename, "");
654
655
for(i = 0; i < specSize; i++) {
656
orthoSpectrum->m_pData[i] = _anisoACF->m_pSpectrum->m_pData[i] / 15.0f;
657
}
658
659
// #ifdef TARGET_LINUX
660
// snprintf(filename, BUF_SIZE, "%s/spectrum_ortho_%s.dat", g_ramanDir, m_sName);
661
// #else
662
// sprintf(filename, "%s/spectrum_ortho_%s.dat", g_ramanDir, m_sName);
663
// #endif
664
// mprintf(" Saving ortho spectrum as %s...\n", filename);
665
// orthoSpectrum->Write("", filename, "");
666
667
for(i = 0; i < specSize; i++) {
668
sumSpectrum->m_pData[i] = paraSpectrum->m_pData[i] + orthoSpectrum->m_pData[i];
669
}
670
671
// #ifdef TARGET_LINUX
672
// snprintf(filename, BUF_SIZE, "%s/spectrum_unpol_%s.dat", g_ramanDir, m_sName);
673
// #else
674
// sprintf(filename, "%s/spectrum_unpol_%s.dat", g_ramanDir, m_sName);
675
// #endif
676
// mprintf(" Saving unpolarized spectrum as %s...\n", filename);
677
// sumSpectrum->Write("", filename, "");
678
679
for(i = 0; i < specSize; i++) {
680
float freq = i * paraSpectrum->m_fMaxRWL / paraSpectrum->m_iSize;
681
float crossSecFactor = powf(g_laser - freq, 4) / freq / (1 - expf(-1.438777f * freq / g_temp));
682
paraSpectrum->m_pData[i] *= crossSecFactor;
683
orthoSpectrum->m_pData[i] *= crossSecFactor;
684
sumSpectrum->m_pData[i] *= crossSecFactor;
685
}
686
687
for(i = 0; i < specSize; i++) {
688
depolSpectrum->m_pData[i] = orthoSpectrum->m_pData[i] / paraSpectrum->m_pData[i];
689
}
690
691
#ifdef TARGET_LINUX
692
snprintf(filename, BUF_SIZE, "%s/spectrum_para_%s.csv", g_ramanDir, m_sName);
693
#else
694
sprintf(filename, "%s/spectrum_para_%s.csv", g_ramanDir, m_sName);
695
#endif
696
mprintf(" Saving parallel spectrum as %s...\n", filename);
697
paraSpectrum->Write("", filename, "");
698
699
#ifdef TARGET_LINUX
700
snprintf(filename, BUF_SIZE, "%s/spectrum_ortho_%s.csv", g_ramanDir, m_sName);
701
#else
702
sprintf(filename, "%s/spectrum_ortho_%s.csv", g_ramanDir, m_sName);
703
#endif
704
mprintf(" Saving orthogonal spectrum as %s...\n", filename);
705
orthoSpectrum->Write("", filename, "");
706
707
#ifdef TARGET_LINUX
708
snprintf(filename, BUF_SIZE, "%s/spectrum_unpol_%s.csv", g_ramanDir, m_sName);
709
#else
710
sprintf(filename, "%s/spectrum_unpol_%s.csv", g_ramanDir, m_sName);
711
#endif
712
mprintf(" Saving unpolarized spectrum as %s...\n", filename);
713
sumSpectrum->Write("", filename, "");
714
715
#ifdef TARGET_LINUX
716
snprintf(filename, BUF_SIZE, "%s/depol_ratio_%s.csv", g_ramanDir, m_sName);
717
#else
718
sprintf(filename, "%s/depol_ratio_%s.csv", g_ramanDir, m_sName);
719
#endif
720
mprintf(" Saving depolarization ratio as %s...\n", filename);
721
depolSpectrum->Write("", filename, "");
722
723
delete anisotropyACFSum;
724
delete anisotropyACF;
725
delete anisotropyPol;
726
727
delete paraSpectrum;
728
delete orthoSpectrum;
729
delete sumSpectrum;
730
delete depolSpectrum;
731
} else {
732
CxFloatArray *ACF;
733
try { ACF = new CxFloatArray("CRamanDyn::finalize():ACF"); } catch(...) { ACF = NULL; }
734
if(ACF == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
735
ACF->SetSize(m_iDepth);
736
737
mprintf(" Parallel part...\n");
738
mprintf(WHITE, " [");
739
for(i = 0; i < m_iDepth; i++) {
740
_isoACF->m_pData[i] = 0.0f;
741
}
742
for(i = 0; i < m_iMolecules; i++) {
743
if(fmodf((float)i, step) < 1.0f)
744
mprintf(WHITE, "#");
745
autoCorrelation->AutoCorrelate((CxFloatArray *)_polarizability[0][0][i], ACF);
746
for(int j = 0; j < m_iDepth; j++) {
747
_isoACF->m_pData[j] += (*ACF)[j];
748
}
749
}
750
for(i = 0; i < m_iDepth; i++) {
751
if(_global)
752
_isoACF->m_pData[i] /= g_iSteps / g_stride;
753
else
754
_isoACF->m_pData[i] /= g_iSteps / g_stride * m_iMolecules;
755
}
756
mprintf(WHITE, "]\n");
757
758
// #ifdef TARGET_LINUX
759
// snprintf(filename, BUF_SIZE, "%s/acf_para_%s.dat", g_ramanDir, m_sName);
760
// #else
761
// sprintf(filename, "%s/acf_para_%s.dat", g_ramanDir, m_sName);
762
// #endif
763
// mprintf(" Saving ACF as %s...\n", filename);
764
// FILE *acf_file = OpenFileWrite(filename, false);
765
// for(i = 0; i < m_iDepth; i++)
766
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _isoACF->m_pData[i]);
767
// fclose(acf_file);
768
769
if(_isoACF->m_iMirror != 0) {
770
mprintf(" Mirroring ACF...\n");
771
_isoACF->Mirror(_isoACF->m_iMirror);
772
}
773
if(_isoACF->m_bWindowFunction != 0) {
774
mprintf(" Applying window function to ACF...\n");
775
_isoACF->Window();
776
}
777
778
// #ifdef TARGET_LINUX
779
// snprintf(filename, BUF_SIZE, "%s/acf_para_%s.mw.dat", g_ramanDir, m_sName);
780
// #else
781
// sprintf(filename, "%s/acf_para_%s.mw.dat", g_ramanDir, m_sName);
782
// #endif
783
// mprintf(" Saving ACF as %s...\n", filename);
784
// acf_file = OpenFileWrite(filename, false);
785
// for(i = 0; i < _isoACF->m_iSize; i++)
786
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _isoACF->m_pData[i]);
787
// fclose(acf_file);
788
789
mprintf(" Performing Fourier transformation...\n");
790
CFFT *fft;
791
try { fft = new CFFT(); } catch(...) { fft = NULL; }
792
if(fft == NULL) NewException((double)sizeof(CFFT), __FILE__, __LINE__, __PRETTY_FUNCTION__);
793
fft->PrepareFFT_C2C(_isoACF->m_iSize + _isoACF->m_iZeroPadding);
794
_isoACF->Transform(fft);
795
delete fft;
796
_isoACF->m_pSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
797
798
// #ifdef TARGET_LINUX
799
// snprintf(filename, BUF_SIZE, "%s/spectrum_para_%s.dat", g_ramanDir, m_sName);
800
// #else
801
// sprintf(filename, "%s/spectrum_para_%s.dat", g_ramanDir, m_sName);
802
// #endif
803
// mprintf(" Saving spectrum as %s...\n", filename);
804
// _isoACF->m_pSpectrum->Write("", filename, "");
805
806
mprintf(" Orthogonal part...\n");
807
mprintf(WHITE, " [");
808
for(i = 0; i < m_iDepth; i++) {
809
_anisoACF->m_pData[i] = 0.0f;
810
}
811
for(i = 0; i < m_iMolecules; i++) {
812
if(fmodf((float)i, step) < 1.0f)
813
mprintf(WHITE, "#");
814
autoCorrelation->AutoCorrelate((CxFloatArray *)_polarizability[1][0][i], ACF);
815
for(int j = 0; j < m_iDepth; j++) {
816
_anisoACF->m_pData[j] += (*ACF)[j];
817
}
818
}
819
for(i = 0; i < m_iDepth; i++) {
820
if(_global)
821
_anisoACF->m_pData[i] /= g_iSteps / g_stride;
822
else
823
_anisoACF->m_pData[i] /= g_iSteps / g_stride * m_iMolecules;
824
}
825
mprintf(WHITE, "]\n");
826
827
// #ifdef TARGET_LINUX
828
// snprintf(filename, BUF_SIZE, "%s/acf_ortho_%s.dat", g_ramanDir, m_sName);
829
// #else
830
// sprintf(filename, "%s/acf_ortho_%s.dat", g_ramanDir, m_sName);
831
// #endif
832
// mprintf(" Saving ACF as %s...\n", filename);
833
// acf_file = OpenFileWrite(filename, false);
834
// for(i = 0; i < m_iDepth; i++)
835
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _anisoACF->m_pData[i]);
836
// fclose(acf_file);
837
838
if(_anisoACF->m_iMirror != 0) {
839
mprintf(" Mirroring ACF...\n");
840
_anisoACF->Mirror(_isoACF->m_iMirror);
841
}
842
if(_anisoACF->m_bWindowFunction != 0) {
843
mprintf(" Applying window function to ACF...\n");
844
_anisoACF->Window();
845
}
846
847
// #ifdef TARGET_LINUX
848
// snprintf(filename, BUF_SIZE, "%s/acf_ortho_%s.mw.dat", g_ramanDir, m_sName);
849
// #else
850
// sprintf(filename, "%s/acf_ortho_%s.mw.dat", g_ramanDir, m_sName);
851
// #endif
852
// mprintf(" Saving ACF as %s...\n", filename);
853
// acf_file = OpenFileWrite(filename, false);
854
// for(i = 0; i < _isoACF->m_iSize; i++)
855
// fprintf(acf_file, "%10.2f %12.8f\n", i * g_fTimestepLength * g_iStride * g_stride, _anisoACF->m_pData[i]);
856
// fclose(acf_file);
857
858
mprintf(" Performing Fourier transformation...\n");
859
try { fft = new CFFT(); } catch(...) { fft = NULL; }
860
if(fft == NULL) NewException((double)sizeof(CFFT), __FILE__, __LINE__, __PRETTY_FUNCTION__);
861
fft->PrepareFFT_C2C(_anisoACF->m_iSize + _anisoACF->m_iZeroPadding);
862
_anisoACF->Transform(fft);
863
delete fft;
864
_anisoACF->m_pSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
865
866
// #ifdef TARGET_LINUX
867
// snprintf(filename, BUF_SIZE, "%s/spectrum_ortho_%s.dat", g_ramanDir, m_sName);
868
// #else
869
// sprintf(filename, "%s/spectrum_ortho_%s.dat", g_ramanDir, m_sName);
870
// #endif
871
// mprintf(" Saving spectrum as %s...\n", filename);
872
// _anisoACF->m_pSpectrum->Write("", filename, "");
873
874
delete ACF;
875
876
int specSize = _isoACF->m_pSpectrum->m_iSize;
877
CSpectrum *sumSpectrum, *depolSpectrum;
878
try { sumSpectrum = new CSpectrum(); } catch(...) { sumSpectrum = NULL; }
879
if(sumSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
880
sumSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
881
sumSpectrum->Create(specSize);
882
sumSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
883
try { depolSpectrum = new CSpectrum(); } catch(...) { depolSpectrum = NULL; }
884
if(depolSpectrum == NULL) NewException((double)sizeof(CSpectrum), __FILE__, __LINE__, __PRETTY_FUNCTION__);
885
depolSpectrum->m_fWaveNumber = _isoACF->m_fSpecWaveNumber;
886
depolSpectrum->Create(specSize);
887
depolSpectrum->SetMaxRWL(1e15f/299792458.0f/100.0f/g_fTimestepLength/g_iStride/g_stride);
888
889
for(i = 0; i < specSize; i++) {
890
sumSpectrum->m_pData[i] = _isoACF->m_pSpectrum->m_pData[i] + _anisoACF->m_pSpectrum->m_pData[i];
891
}
892
893
// #ifdef TARGET_LINUX
894
// snprintf(filename, BUF_SIZE, "%s/spectrum_unpol_%s.dat", g_ramanDir, m_sName);
895
// #else
896
// sprintf(filename, "%s/spectrum_unpol_%s.dat", g_ramanDir, m_sName);
897
// #endif
898
// mprintf(" Saving unpolarized spectrum as %s...\n", filename);
899
// sumSpectrum->Write("", filename, "");
900
901
for(i = 0; i < specSize; i++) {
902
float freq = i * _isoACF->m_pSpectrum->m_fMaxRWL / _isoACF->m_pSpectrum->m_iSize;
903
float crossSecFactor = powf(g_laser - freq, 4) / freq / (1 - expf(-1.438777f * freq / g_temp));
904
_isoACF->m_pSpectrum->m_pData[i] *= crossSecFactor;
905
_anisoACF->m_pSpectrum->m_pData[i] *= crossSecFactor;
906
sumSpectrum->m_pData[i] *= crossSecFactor;
907
}
908
909
for(i = 0; i < specSize; i++) {
910
depolSpectrum->m_pData[i] = _anisoACF->m_pSpectrum->m_pData[i] / _isoACF->m_pSpectrum->m_pData[i];
911
}
912
913
#ifdef TARGET_LINUX
914
snprintf(filename, BUF_SIZE, "%s/spectrum_para_%s.csv", g_ramanDir, m_sName);
915
#else
916
sprintf(filename, "%s/spectrum_para_%s.csv", g_ramanDir, m_sName);
917
#endif
918
mprintf(" Saving parallel spectrum as %s...\n", filename);
919
_isoACF->m_pSpectrum->Write("", filename, "");
920
921
#ifdef TARGET_LINUX
922
snprintf(filename, BUF_SIZE, "%s/spectrum_ortho_%s.csv", g_ramanDir, m_sName);
923
#else
924
sprintf(filename, "%s/spectrum_ortho_%s.csv", g_ramanDir, m_sName);
925
#endif
926
mprintf(" Saving orthogonal spectrum as %s...\n", filename);
927
_anisoACF->m_pSpectrum->Write("", filename, "");
928
929
#ifdef TARGET_LINUX
930
snprintf(filename, BUF_SIZE, "%s/spectrum_unpol_%s.csv", g_ramanDir, m_sName);
931
#else
932
sprintf(filename, "%s/spectrum_unpol_%s.csv", g_ramanDir, m_sName);
933
#endif
934
mprintf(" Saving unpolarized spectrum as %s...\n", filename);
935
sumSpectrum->Write("", filename, "");
936
937
#ifdef TARGET_LINUX
938
snprintf(filename, BUF_SIZE, "%s/depol_ratio_%s.csv", g_ramanDir, m_sName);
939
#else
940
sprintf(filename, "%s/depol_ratio_%s.csv", g_ramanDir, m_sName);
941
#endif
942
mprintf(" Saving depolarization ratio as %s...\n", filename);
943
depolSpectrum->Write("", filename, "");
944
945
delete sumSpectrum;
946
delete depolSpectrum;
947
}
948
949
delete autoCorrelation;
950
}
951
952
CRamanObservation::CRamanObservation(bool global) {
953
int i;
954
955
m_bTimeDev = false;
956
m_bTimeDiff = false;
957
m_bSelf = false;
958
m_bOthers = true;
959
960
if(global) {
961
m_iShowMol = -1;
962
m_iShowMolCount = g_oaSingleMolecules.GetSize();
963
m_bObsCertain = false;
964
m_bDecompDist = false;
965
m_bDecompType = false;
966
} else {
967
char buf[BUF_SIZE];
968
char buf2[BUF_SIZE];
969
size_t remaining = BUF_SIZE;
970
if(g_oaMolecules.GetSize() > 1) {
971
#ifdef TARGET_LINUX
972
remaining -= snprintf(buf, remaining, " Which molecule should be observed (");
973
#else
974
remaining -= sprintf(buf, " Which molecule should be observed (");
975
#endif
976
for(i = 0; i < g_oaMolecules.GetSize(); i++) {
977
if(remaining < 1)
978
break;
979
#ifdef TARGET_LINUX
980
size_t length = snprintf(buf2, remaining, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
981
#else
982
size_t length = sprintf(buf2, "%s=%d", ((CMolecule *)g_oaMolecules[i])->m_sName, i+1);
983
#endif
984
strncat(buf, buf2, remaining - 1);
985
remaining -= length;
986
if(i < g_oaMolecules.GetSize() - 1) {
987
#ifdef TARGET_LINUX
988
length = snprintf(buf2, remaining, ", ");
989
#else
990
length = sprintf(buf2, ", ");
991
#endif
992
strncat(buf, buf2, remaining - 1);
993
remaining -= length;
994
}
995
}
996
strncat(buf, ")? ", remaining - 1);
997
m_iShowMol = AskRangeInteger_ND(buf, 1, g_oaMolecules.GetSize()) - 1;
998
} else {
999
m_iShowMol = 0;
1000
}
1001
m_iShowMolCount = ((CMolecule *)g_oaMolecules[m_iShowMol])->m_laSingleMolIndex.GetSize();
1002
m_bObsCertain = false;
1003
m_bDecompDist = false;
1004
m_bDecompType = false;
1005
mprintf("\n");
1006
}
1007
1008
try { _ramanDyn = new CRamanDyn(m_iShowMol, global); } catch(...) { _ramanDyn = NULL; }
1009
if(_ramanDyn == NULL) NewException((double)sizeof(CRamanDyn), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1010
}
1011
1012
CRamanObservation::~CRamanObservation() {
1013
delete _ramanDyn;
1014
}
1015
1016
void CRamanObservation::initialize() {
1017
_ramanDyn->initialize();
1018
}
1019
1020
void CRamanObservation::getDipole0() {
1021
_ramanDyn->getDipole0();
1022
}
1023
1024
void CRamanObservation::calcPolarizability(int fieldDirection) {
1025
_ramanDyn->calcPolarizability(fieldDirection);
1026
}
1027
1028
void CRamanObservation::finalize() {
1029
_ramanDyn->finalize();
1030
}
1031
1032
static bool parseSettings(FILE *settingsFile) {
1033
char buf[BUF_SIZE];
1034
1035
if(fgets(buf, BUF_SIZE, settingsFile) == NULL)
1036
return false;
1037
int temp = -1;
1038
if(sscanf(buf, "%d", &temp) < 1)
1039
return false;
1040
if(temp == 0)
1041
g_orientAvg = false;
1042
else if(temp == 1)
1043
g_orientAvg = true;
1044
else
1045
return false;
1046
1047
if(fgets(buf, BUF_SIZE, settingsFile) == NULL)
1048
return false;
1049
if(sscanf(buf, "%f", &g_fieldStrength) < 1)
1050
return false;
1051
1052
if(fgets(buf, BUF_SIZE, settingsFile) == NULL)
1053
return false;
1054
if(sscanf(buf, "%d", &g_stride) < 1)
1055
return false;
1056
1057
return true;
1058
}
1059
1060
bool gatherRaman() {
1061
#ifndef TARGET_LINUX
1062
mprintf(RED, "Raman calculations are currently possible only under Linux.\n");
1063
return false;
1064
#else
1065
mprintf(" To calculate Raman spectra, polarizabilities have to be determined.\n");
1066
mprintf(" TRAVIS creates CP2K input files for numerical polarizabilities\n and will process the resulting data in a second run.\n\n");
1067
1068
g_newRaman = AskYesNo(" Do you wish to create new CP2K input files (y) or process existing results (n)? [y] ", true);
1069
char buf[BUF_SIZE];
1070
if(g_newRaman) {
1071
AskString(" Please enter a name for the directory to collect the data: [raman] ", buf, "raman");
1072
} else {
1073
AskString(" Please enter the name of the directory to take the data from: [raman] ", buf, "raman");
1074
}
1075
try { g_ramanDir = new char[strlen(buf)+1]; } catch(...) { g_ramanDir = NULL; }
1076
if(g_ramanDir == NULL) NewException((double)(strlen(buf)+1)*sizeof(char), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1077
strcpy(g_ramanDir, buf);
1078
1079
if(g_newRaman) {
1080
g_bKeepOriginalCoords = true;
1081
if(FileExist(g_ramanDir)) {
1082
mprintf(RED, "A file or a directory \"%s\" already exists. Please remove it first.\n", g_ramanDir);
1083
return false;
1084
}
1085
if(mkdir(g_ramanDir, S_IRWXU) != 0) {
1086
mprintf(RED, "Directory \"%s\" could not be created: %s\n", g_ramanDir, strerror(errno));
1087
return false;
1088
}
1089
1090
char filename[BUF_SIZE];
1091
snprintf(filename, BUF_SIZE, "%s/settings.dat", g_ramanDir);
1092
FILE *settingsFile = fopen(filename, "w");
1093
if(settingsFile == NULL) {
1094
mprintf(RED, "Could not open settings file \"%s\": %s\n", filename, strerror(errno));
1095
return false;
1096
}
1097
1098
g_orientAvg = AskYesNo("\n Use orientational averaging? [no] ", false);
1099
if(g_orientAvg)
1100
fprintf(settingsFile, "%d\n", 1);
1101
else
1102
fprintf(settingsFile, "%d\n", 0);
1103
1104
g_fieldStrength = AskFloat(" Field strength in atomic units [5.0e-4] ", 5.0e-4);
1105
fprintf(settingsFile, "%.6E\n", g_fieldStrength);
1106
1107
g_stride = AskInteger(" Calculate polarizability for every n-th timestep [1] ", 1);
1108
fprintf(settingsFile, "%d\n", g_stride);
1109
1110
fclose(settingsFile);
1111
1112
snprintf(filename, BUF_SIZE, "%s/1", g_ramanDir);
1113
if(mkdir(filename, S_IRWXU) != 0) {
1114
mprintf(RED, "Directory \"%s\" could not be created: %s\n", filename, strerror(errno));
1115
return false;
1116
}
1117
if(g_orientAvg) {
1118
snprintf(filename, BUF_SIZE, "%s/2", g_ramanDir);
1119
if(mkdir(filename, S_IRWXU) != 0) {
1120
mprintf(RED, "Directory \"%s\" could not be created: %s\n", filename, strerror(errno));
1121
return false;
1122
}
1123
snprintf(filename, BUF_SIZE, "%s/3", g_ramanDir);
1124
if(mkdir(filename, S_IRWXU) != 0) {
1125
mprintf(RED, "Directory \"%s\" could not be created: %s\n", filename, strerror(errno));
1126
return false;
1127
}
1128
}
1129
1130
FILE* templateFile;
1131
snprintf(filename, BUF_SIZE, "%s/template.inp", g_ramanDir);
1132
templateFile = fopen(filename, "w");
1133
if(templateFile == NULL) {
1134
mprintf(RED, "Could not open template file \"%s\": %s\n", filename, strerror(errno));
1135
return false;
1136
}
1137
fprintf(templateFile, "&GLOBAL\n");
1138
fprintf(templateFile, " PROJECT_NAME polarizability\n");
1139
fprintf(templateFile, " RUN_TYPE MD\n");
1140
fprintf(templateFile, " PRINT_LEVEL LOW\n");
1141
fprintf(templateFile, "&END\n");
1142
fprintf(templateFile, "&FORCE_EVAL\n");
1143
fprintf(templateFile, " &DFT\n");
1144
fprintf(templateFile, " BASIS_SET_FILE_NAME BASIS_MOLOPT\n");
1145
fprintf(templateFile, " POTENTIAL_FILE_NAME POTENTIAL\n");
1146
fprintf(templateFile, " &MGRID\n");
1147
fprintf(templateFile, " NGRIDS 5\n");
1148
fprintf(templateFile, " CUTOFF 280\n");
1149
fprintf(templateFile, " REL_CUTOFF 40\n");
1150
fprintf(templateFile, " &END\n");
1151
fprintf(templateFile, " &SCF\n");
1152
fprintf(templateFile, " SCF_GUESS ATOMIC\n");
1153
fprintf(templateFile, " MAX_SCF 200\n");
1154
fprintf(templateFile, " EPS_SCF 1.0E-5\n");
1155
fprintf(templateFile, " &OT\n");
1156
fprintf(templateFile, " MINIMIZER DIIS\n");
1157
fprintf(templateFile, " PRECONDITIONER FULL_SINGLE_INVERSE\n");
1158
fprintf(templateFile, " &END\n");
1159
fprintf(templateFile, " &END\n");
1160
fprintf(templateFile, " &XC\n");
1161
fprintf(templateFile, " &XC_FUNCTIONAL BLYP\n");
1162
fprintf(templateFile, " &END\n");
1163
fprintf(templateFile, " &XC_GRID\n");
1164
fprintf(templateFile, " XC_SMOOTH_RHO NN10\n");
1165
fprintf(templateFile, " XC_DERIV NN10_SMOOTH\n");
1166
fprintf(templateFile, " &END\n");
1167
fprintf(templateFile, " &VDW_POTENTIAL\n");
1168
fprintf(templateFile, " POTENTIAL_TYPE PAIR_POTENTIAL\n");
1169
fprintf(templateFile, " &PAIR_POTENTIAL\n");
1170
fprintf(templateFile, " TYPE DFTD3\n");
1171
fprintf(templateFile, " REFERENCE_FUNCTIONAL BLYP\n");
1172
fprintf(templateFile, " PARAMETER_FILE_NAME dftd3.dat\n");
1173
fprintf(templateFile, " &END\n");
1174
fprintf(templateFile, " &END\n");
1175
fprintf(templateFile, " &END\n");
1176
fprintf(templateFile, " &LOCALIZE\n");
1177
fprintf(templateFile, " METHOD CRAZY\n");
1178
fprintf(templateFile, " MAX_ITER 2000\n");
1179
fprintf(templateFile, " &PRINT\n");
1180
fprintf(templateFile, " &WANNIER_CENTERS\n");
1181
fprintf(templateFile, " IONS+CENTERS\n");
1182
fprintf(templateFile, " FILENAME =polarizability_wannier.xyz\n");
1183
fprintf(templateFile, " &END\n");
1184
fprintf(templateFile, " &END\n");
1185
fprintf(templateFile, " &END\n");
1186
fprintf(templateFile, " &PERIODIC_EFIELD\n");
1187
fprintf(templateFile, " INTENSITY ###!field strength will be put here###\n");
1188
fprintf(templateFile, " POLARISATION ###!polarisation vector will be put here###\n");
1189
fprintf(templateFile, " &END\n");
1190
fprintf(templateFile, " &END\n");
1191
fprintf(templateFile, " &SUBSYS\n");
1192
fprintf(templateFile, " &CELL\n");
1193
fprintf(templateFile, " ABC %.6f %.6f %.6f\n", g_fBoxX / 100.0f, g_fBoxY / 100.0f, g_fBoxZ / 100.0f);
1194
fprintf(templateFile, " &END\n");
1195
fprintf(templateFile, " &COORD\n");
1196
fprintf(templateFile, "###!coordinates will be put here###\n");
1197
fprintf(templateFile, " &END\n");
1198
fprintf(templateFile, " &KIND H\n");
1199
fprintf(templateFile, " BASIS_SET DZVP-MOLOPT-SR-GTH\n");
1200
fprintf(templateFile, " POTENTIAL GTH-BLYP-q1\n");
1201
fprintf(templateFile, " &END\n");
1202
fprintf(templateFile, " &KIND C\n");
1203
fprintf(templateFile, " BASIS_SET DZVP-MOLOPT-SR-GTH\n");
1204
fprintf(templateFile, " POTENTIAL GTH-BLYP-q4\n");
1205
fprintf(templateFile, " &END\n");
1206
fprintf(templateFile, " &KIND N\n");
1207
fprintf(templateFile, " BASIS_SET DZVP-MOLOPT-SR-GTH\n");
1208
fprintf(templateFile, " POTENTIAL GTH-BLYP-q5\n");
1209
fprintf(templateFile, " &END\n");
1210
fprintf(templateFile, " &KIND O\n");
1211
fprintf(templateFile, " BASIS_SET DZVP-MOLOPT-SR-GTH\n");
1212
fprintf(templateFile, " POTENTIAL GTH-BLYP-q6\n");
1213
fprintf(templateFile, " &END\n");
1214
fprintf(templateFile, " &END\n");
1215
fprintf(templateFile, "&END\n");
1216
fprintf(templateFile, "&MOTION\n");
1217
fprintf(templateFile, " &MD\n");
1218
fprintf(templateFile, " ENSEMBLE REFTRAJ\n");
1219
fprintf(templateFile, " STEPS ###!number of steps will be put here###\n");
1220
fprintf(templateFile, " TIMESTEP %f\n", g_fTimestepLength * g_stride);
1221
fprintf(templateFile, " &REFTRAJ\n");
1222
fprintf(templateFile, " EVAL_ENERGY_FORCES\n");
1223
fprintf(templateFile, " TRAJ_FILE_NAME ../polarizability_reftraj.xyz\n");
1224
fprintf(templateFile, " &END\n");
1225
fprintf(templateFile, " &END\n");
1226
fprintf(templateFile, " &PRINT\n");
1227
fprintf(templateFile, " &RESTART\n");
1228
fprintf(templateFile, " &EACH\n");
1229
fprintf(templateFile, " MD 1\n");
1230
fprintf(templateFile, " &END\n");
1231
fprintf(templateFile, " &END\n");
1232
fprintf(templateFile, " &END\n");
1233
fprintf(templateFile, "&END\n");
1234
fclose(templateFile);
1235
1236
mprintf("\n An input template for the polarizability calculations has been written to \"%s/template.inp\"\n Please modify it according to your needs.\n Press any key when you are finished.\n", g_ramanDir);
1237
getchar();
1238
1239
if(g_orientAvg) {
1240
mprintf(" CP2K input files will be created in \"%s/1\", \"%s/2\", and \"%s/3\".\n", g_ramanDir, g_ramanDir, g_ramanDir);
1241
mprintf(" After TRAVIS has finished, please run them and make sure that the resulting trajectories \"polarizability_wannier.xyz\"\n are placed in the same directories before you execute TRAVIS for evaluation.\n\n");
1242
} else {
1243
mprintf(" A CP2K input file will be created in \"%s/1\".\n", g_ramanDir);
1244
mprintf(" After TRAVIS has finished, please run it and make sure that the resulting trajectory \"polarizability_wannier.xyz\"\n is placed in the same directory before you run TRAVIS for evaluation.\n\n");
1245
}
1246
} else {
1247
if(!FileExist(g_ramanDir)) {
1248
mprintf(RED, "The directory \"%s\" was not found.\n", g_ramanDir);
1249
return false;
1250
}
1251
1252
char filename[BUF_SIZE];
1253
snprintf(filename, BUF_SIZE, "%s/settings.dat", g_ramanDir);
1254
FILE *settingsFile;
1255
settingsFile = fopen(filename, "r");
1256
if(settingsFile == NULL) {
1257
mprintf(RED, "Could not open settings file \"%s\": %s\n", filename, strerror(errno));
1258
return false;
1259
}
1260
if(!parseSettings(settingsFile)) {
1261
mprintf(RED, "Could not parse settings file.\n");
1262
return false;
1263
}
1264
1265
mprintf("\n");
1266
if(g_orientAvg)
1267
mprintf(" Using orientational averaging\n");
1268
else
1269
mprintf(" Not using orientational averaging\n");
1270
mprintf(" Field strength: %.6e a. u.\n", g_fieldStrength);
1271
mprintf(" Using every %d%s timestep\n\n", g_stride, (g_stride == 1) ? "st" : ((g_stride == 2) ? "nd" : ((g_stride == 3) ? "rd" : "th")));
1272
1273
g_laser = AskFloat(" Calculate scattering cross section for which laser frequency (cm^-1)? [20000.0] ", 20000.0f);
1274
g_temp = AskFloat(" Calculate scattering cross section for which temperature (K)? [300.0] ", 300.0f);
1275
1276
if(AskYesNo(" Compute Raman spectrum of whole system (y/n)? [yes] ", true)) {
1277
mprintf(YELLOW, "\n>>> Global Raman Observation >>>\n\n");
1278
1279
CRamanObservation *obs;
1280
try { obs = new CRamanObservation(true); } catch(...) { obs = NULL; }
1281
if(obs == NULL) NewException((double)sizeof(CRamanObservation), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1282
g_ramObserv.Add(obs);
1283
1284
mprintf(YELLOW, "<<< End of Global Raman Observation <<<\n\n");
1285
}
1286
1287
if(AskYesNo(" Compute Raman spectra for certain molecule types (y/n)? [no] ", false)) {
1288
while(true) {
1289
mprintf(YELLOW, "\n>>> Raman Observation %d >>>\n\n", g_ramObserv.GetSize() + 1);
1290
1291
CRamanObservation *obs;
1292
try { obs = new CRamanObservation(); } catch(...) { obs = NULL; }
1293
if(obs == NULL) NewException((double)sizeof(CRamanObservation), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1294
g_ramObserv.Add(obs);
1295
1296
mprintf(YELLOW, "<<< End of Raman Observation %d <<<\n\n", g_ramObserv.GetSize());
1297
1298
if(!AskYesNo(" Add another observation (y/n)? [no] ", false))
1299
break;
1300
mprintf("\n");
1301
}
1302
}
1303
mprintf("\n");
1304
}
1305
1306
return true;
1307
#endif
1308
}
1309
1310
bool initializeRaman() {
1311
int i;
1312
1313
if(g_newRaman) {
1314
char filename[BUF_SIZE];
1315
#ifdef TARGET_LINUX
1316
snprintf(filename, BUF_SIZE, "%s/template.inp", g_ramanDir);
1317
#else
1318
sprintf(filename, "%s/template.inp", g_ramanDir);
1319
#endif
1320
FILE *templateFile;
1321
templateFile = fopen(filename, "r");
1322
if(templateFile == NULL) {
1323
mprintf(RED, "Could not open template file \"%s\": %s\n", filename, strerror(errno));
1324
return false;
1325
}
1326
fseek(templateFile, 0L, SEEK_END);
1327
long length = ftell(templateFile);
1328
if(length < 0) {
1329
mprintf(RED, "Could not determine size of template file: %s\n", strerror(errno));
1330
return false;
1331
}
1332
1333
try { g_inputTemplate = new char[length + 1]; } catch(...) { g_inputTemplate = NULL; }
1334
if(g_inputTemplate == NULL) NewException((double)(length + 1)*sizeof(char), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1335
1336
rewind(templateFile);
1337
if((long)fread(g_inputTemplate, sizeof(char), length, templateFile) < length) {
1338
mprintf(RED, "Could not read template file: %s\n", strerror(errno));
1339
return false;
1340
}
1341
g_inputTemplate[length] = '\0';
1342
1343
fclose(templateFile);
1344
1345
g_templateFieldPos = strstr(g_inputTemplate, "###!field strength will be put here###");
1346
if(g_templateFieldPos == NULL) {
1347
mprintf(RED, "Position mark for field strength missing in template.\n");
1348
return false;
1349
}
1350
g_templatePolPos = strstr(g_inputTemplate, "###!polarisation vector will be put here###");
1351
if(g_templatePolPos == NULL) {
1352
mprintf(RED, "Position mark for polarisation vector missing in template.\n");
1353
return false;
1354
}
1355
g_templateCoordPos = strstr(g_inputTemplate, "###!coordinates will be put here###");
1356
if(g_templateCoordPos == NULL) {
1357
mprintf(RED, "Position mark for coordinates missing in template.\n");
1358
return false;
1359
}
1360
g_templateStepsPos = strstr(g_inputTemplate, "###!number of steps will be put here###");
1361
if(g_templateStepsPos == NULL) {
1362
mprintf(RED, "Position mark for number of steps missing in template.\n");
1363
return false;
1364
}
1365
g_templateFieldPos[0] = '\0';
1366
g_templatePolPos[0] = '\0';
1367
g_templateCoordPos[0] = '\0';
1368
g_templateStepsPos[0] = '\0';
1369
1370
#ifdef TARGET_LINUX
1371
snprintf(filename, BUF_SIZE, "%s/polarizability_reftraj.xyz", g_ramanDir);
1372
#else
1373
sprintf(filename, "%s/polarizability_reftraj.xyz", g_ramanDir);
1374
#endif
1375
g_reftrajFile = fopen(filename, "w");
1376
if(g_reftrajFile == NULL) {
1377
mprintf(RED, "Could not open reference trajectory \"%s\": %s\n", filename, strerror(errno));
1378
return false;
1379
}
1380
} else {
1381
for(i = 0; i < g_ramObserv.GetSize(); i++) {
1382
mprintf("Initializing Raman Observation %d...\n", i+1);
1383
CRamanObservation *obs = (CRamanObservation *)g_ramObserv[i];
1384
obs->initialize();
1385
}
1386
1387
char filename[BUF_SIZE];
1388
for(i = 0; i < (g_orientAvg ? 3 : 1); i++) {
1389
#ifdef TARGET_LINUX
1390
snprintf(filename, BUF_SIZE, "%s/%d/polarizability_wannier.xyz", g_ramanDir, i + 1);
1391
#else
1392
sprintf(filename, "%s/%d/polarizability_wannier.xyz", g_ramanDir, i + 1);
1393
#endif
1394
g_polFile[i] = fopen(filename, "r");
1395
if(g_polFile[i] == NULL) {
1396
mprintf(RED, "Could not open trajectory \"%s\": %s\n", filename, strerror(errno));
1397
return false;
1398
}
1399
}
1400
1401
for(i = 0; i < (g_orientAvg ? 3 : 1); i++) {
1402
try { g_timestep[i] = new CTimeStep(); } catch(...) { g_timestep[i] = NULL; }
1403
if(g_timestep[i] == NULL) NewException((double)sizeof(CTimeStep), __FILE__, __LINE__, __PRETTY_FUNCTION__);
1404
}
1405
}
1406
1407
return true;
1408
}
1409
1410
void processRaman(CTimeStep *ts) {
1411
int i, j;
1412
g_step++;
1413
if(g_newRaman) {
1414
if(g_step % g_stride == 0) {
1415
g_steps++;
1416
int numAtoms = 0;
1417
for(i = 0; i < g_iGesAtomCount; i++)
1418
if(g_waAtomMolIndex[i] != 60000)
1419
numAtoms++;
1420
fprintf(g_reftrajFile, "%d\nStep %d\n", numAtoms, g_step);
1421
for(i = 0; i < g_iGesAtomCount; i++) {
1422
if(g_waAtomMolIndex[i] == 60000)
1423
continue;
1424
fprintf(g_reftrajFile, "%4s %14.10f %14.10f %14.10f\n", ((CAtom *)g_oaAtoms[g_waAtomRealElement[i]])->m_sName, ts->m_vaCoords_Original[i][0] / 100.0f, ts->m_vaCoords_Original[i][1] / 100.0f, ts->m_vaCoords_Original[i][2] / 100.0f);
1425
}
1426
}
1427
} else {
1428
if(g_step % g_stride == 0) {
1429
for(i = 0; i < g_ramObserv.GetSize(); i++)
1430
((CRamanObservation *)g_ramObserv[i])->getDipole0();
1431
for(i = 0; i < (g_orientAvg ? 3 : 1); i++) {
1432
if(!g_timestep[i]->ReadTimestep(g_polFile[i], false)) {
1433
mprintf(RED, "Error while reading trajectory for polarizabilities.\n");
1434
abort();
1435
}
1436
if(!g_bSaveCoordsUnchanged) {
1437
g_timestep[i]->UniteMolecules(false);
1438
if(g_bRemoveCOM)
1439
g_timestep[i]->CenterCOM();
1440
}
1441
g_timestep[i]->CalcCenters();
1442
if(g_bWannier)
1443
g_timestep[i]->ScanWannier(false);
1444
g_timestep[i]->CalcDipoles();
1445
for(j = 0; j < g_ramObserv.GetSize(); j++) {
1446
((CRamanObservation *)g_ramObserv[j])->calcPolarizability(i);
1447
}
1448
}
1449
}
1450
}
1451
}
1452
1453
void finalizeRaman() {
1454
int i, j, k, l;
1455
if(g_newRaman) {
1456
fclose(g_reftrajFile);
1457
1458
char filename[BUF_SIZE];
1459
FILE *inputFile;
1460
for(i = 1; i < (g_orientAvg ? 4 : 2); i++) {
1461
#ifdef TARGET_LINUX
1462
snprintf(filename, BUF_SIZE, "%s/%d/polarizability.inp", g_ramanDir, i);
1463
#else
1464
sprintf(filename, "%s/%d/polarizability.inp", g_ramanDir, i);
1465
#endif
1466
inputFile = fopen(filename, "w");
1467
if(inputFile == NULL) {
1468
mprintf(RED, "Could not open input file \"%s\": %s\n", filename, strerror(errno));
1469
abort();
1470
}
1471
1472
fprintf(inputFile, "%s", g_inputTemplate);
1473
for(j = 0; j < 4; j++) {
1474
char *tempPos = g_inputTemplate;
1475
for(k = 0; k <= j; k++)
1476
tempPos = strchr(&tempPos[1], '\0');
1477
if(tempPos == g_templateFieldPos) {
1478
fprintf(inputFile, "%.6E", g_fieldStrength);
1479
fprintf(inputFile, "%s", &g_templateFieldPos[38]);
1480
} else if(tempPos == g_templatePolPos) {
1481
if(i == 1)
1482
fprintf(inputFile, "1.0 0.0 0.0");
1483
if(i == 2)
1484
fprintf(inputFile, "0.0 1.0 0.0");
1485
if(i == 3)
1486
fprintf(inputFile, "0.0 0.0 1.0");
1487
fprintf(inputFile, "%s", &g_templatePolPos[43]);
1488
} else if(tempPos == g_templateCoordPos) {
1489
CTimeStep *ts = GetTimeStep(0);
1490
for(l = 0; l < g_iGesAtomCount; l++) {
1491
if(g_waAtomMolIndex[l] == 60000)
1492
continue;
1493
fprintf(inputFile, " %s %14.10f %14.10f %14.10f\n", ((CAtom *)g_oaAtoms[g_waAtomRealElement[l]])->m_sName, ts->m_vaCoords_Original[l][0] / 100.0f, ts->m_vaCoords_Original[l][1] / 100.0f, ts->m_vaCoords_Original[l][2] / 100.0f);
1494
}
1495
fprintf(inputFile, "%s", &g_templateCoordPos[36]);
1496
} else if(tempPos == g_templateStepsPos) {
1497
fprintf(inputFile, "%d", g_steps);
1498
fprintf(inputFile, "%s", &g_templateStepsPos[39]);
1499
} else {
1500
mprintf(RED, "Unexpected error while processing the input template.\n");
1501
}
1502
}
1503
fclose(inputFile);
1504
}
1505
delete[] g_inputTemplate;
1506
} else {
1507
for(i = 0; i < g_ramObserv.GetSize(); i++) {
1508
mprintf(YELLOW, ">>> Raman Observation %d >>>\n\n", i + 1);
1509
((CRamanObservation *)g_ramObserv[i])->finalize();
1510
mprintf(YELLOW, "\n<<< End of Raman Observation %d <<<\n\n", i + 1);
1511
}
1512
1513
for(i = 0; i < (g_orientAvg ? 3 : 1); i++)
1514
fclose(g_polFile[i]);
1515
for(i = 0; i < (g_orientAvg ? 3 : 1); i++)
1516
delete g_timestep[i];
1517
}
1518
delete[] g_ramanDir;
1519
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1