forked from rikigigi/analisi
-
Notifications
You must be signed in to change notification settings - Fork 0
/
correlatorespaziale.cpp_copia
177 lines (144 loc) · 7.38 KB
/
correlatorespaziale.cpp_copia
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
#include "correlatorespaziale.h"
#include "traiettoria.h"
#include "operazionisulista.h"
#include <cmath>
#include "config.h"
#ifdef HAVEfftw3
#include <fftw3.h>
#else
#include <fftw.h>
#endif
#define PI 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628
CorrelatoreSpaziale::CorrelatoreSpaziale(Traiettoria *t,
unsigned int nk,
double sigma2,
unsigned int nthreads,
unsigned int skip,
bool debug) :
sfac(0), nk(nk), sigma2(sigma2),nthreads(nthreads),debug(debug),skip(skip),tipi_atomi(0),t(t)
{
}
fftw_plan CorrelatoreSpaziale::fftw3;
void CorrelatoreSpaziale::reset(const unsigned int numeroTimestepsPerBlocco) {
tipi_atomi=t->get_ntypes();
if (sfac==0){
sfac = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*3*nk*nk*(nk/2+1)*tipi_atomi*(tipi_atomi+1)/2);
//lista = new double[nk*nk*(nk/2+1)*tipi_atomi*(tipi_atomi+1)/2*3];
lunghezza_lista=nk*nk*nk*tipi_atomi*(tipi_atomi+1)/2*3;
lista = new double[lunghezza_lista];
// lista = (double *) fftw_malloc(sizeof(double)*nk*nk*nk*tipi_atomi*(tipi_atomi+1)/2*3);
}
ntimesteps=numeroTimestepsPerBlocco;
// std::cerr<< "Sto usando il CorrSpaziale Davide"<<std::endl;
}
void CorrelatoreSpaziale::s_fac_k(double k[3], unsigned int i_t, fftw_complex * out1 ) {
for (unsigned int i_at=0;i_at<t->get_natoms();i_at++){
double * pos = t->posizioni(i_t,i_at);
double * vel = t->velocita(i_t,i_at);
unsigned int type=t->get_type(i_at);
//std::cerr << type << "\n";
double arg=k[0]*pos[0]+k[1]*pos[1]+k[2]*pos[2];
double s=sin(arg),c=cos(arg);
for (unsigned int icoord=0;icoord<3;icoord++){
out1[3*type+icoord][0]+=c*vel[icoord];
out1[3*type+icoord][1]+=s*vel[icoord];
}
}
}
void CorrelatoreSpaziale::calcola(unsigned int primo) {
//size = nk*nk*(nk/2+1)*tipi_atomi*(tipi_atomi+1)/2;
size = nk*nk*nk*tipi_atomi*(tipi_atomi+1)/2;
//int size_half=nk*nk*(nk/2+1)*tipi_atomi*(tipi_atomi+1)/2;
int size_half=nk*nk*(nk/2+1)*tipi_atomi*(tipi_atomi+1)/2;
int itime = 0;
fftw_complex * sfac_t=(fftw_complex *) fftw_malloc(sizeof(fftw_complex) * tipi_atomi*3);
int n[]={nk,nk,nk};
//int n[]={1,1,1};
//int n[]={nk,nk,nk/2+1};
fftw3 = fftw_plan_many_dft_c2r(3, // rango della trasformata (3D)
n, // lunghezza di ciascuna trasformata
tipi_atomi*(tipi_atomi+1)/2, // numero di trasformate
/* ****** input ****** */
sfac,//traiettoria->velocita_inizio(), // puntatore ai dati
NULL, // i dati sono tutti compatti, non fanno parte di array più grandi
tipi_atomi*(tipi_atomi+1)/2, // la distanza fra due dati successivi
1, // la distanza fra due serie di dati adiacenti
/* ****** output ****** */
lista, // puntatore all'array di output
NULL,
tipi_atomi*(tipi_atomi+1)/2, // la distanza fra due dati successivi
1, // la distanza fra due serie di dati adiacenti
FFTW_DESTROY_INPUT
);
for (unsigned int i=0;i<3*nk*nk*(nk/2+1)*tipi_atomi*(tipi_atomi+1)/2;i++) {
sfac[i][0]=0.0;
sfac[i][1]=0.0;
}
for (unsigned int itimestep=primo;itimestep<ntimesteps+primo;itimestep+=skip){
itime += 1;
double * box=t->scatola(itimestep);
//ciclo su kx,ky e kz. kz va da 0 a kz/2+1
double k[3]={0.0};
for (int kx=0;kx<nk;kx++){
unsigned int ikx=kx*nk*(nk/2+1)*(tipi_atomi+1)*tipi_atomi/2;
k[0]=PI/(box[1]-box[0])*kx;
for (int ky=0;ky<nk;ky++) {
unsigned int iky=ky*(nk/2+1)*(tipi_atomi+1)*tipi_atomi/2;
k[1]=PI/(box[1]-box[0])*ky;
for (int kz=0;kz<nk/2+1;kz++) {
for (unsigned int i=0;i<tipi_atomi*3;i++) {
sfac_t[i][0]=0.0;
sfac_t[i][1]=0.0;
}
unsigned int ik=ikx+iky+kz*(tipi_atomi+1)*tipi_atomi/2;
//k[2]=PI/(box[1]-box[0])*ky; ///// MERDA IL COPIA/INCOLLA
k[2]=PI/(box[1]-box[0])*kz;
// sfac_t[0][0]=0.1 ;/////// TOGLIERE
s_fac_k(k,itimestep,sfac_t);
unsigned int m=0;
for (unsigned int i=0;i<tipi_atomi;i++)
for (unsigned int j=0;j<=i;j++){
//unsigned int m = (tipi_atomi+1)*tipi_atomi/2 - (i+1)*(i+2)/2 + j;
sfac[ik+m][0]+= ( sfac_t[i*3][0]* sfac_t[j*3][0] +sfac_t[i*3][1]* sfac_t[j*3][1] -sfac[ik+m][0])/(itime);
sfac[ik+m][1]+= (-sfac_t[i*3][0]* sfac_t[j*3][1] + sfac_t[i*3][1]* sfac_t[j*3][0] -sfac[ik+m][1])/(itime);
sfac[size_half +ik+m][0]+= ( sfac_t[i*3+1][0]*sfac_t[j*3+1][0]+sfac_t[i*3+1][1]*sfac_t[j*3+1][1] - sfac[ size_half+ik+m][0])/(itime);
sfac[size_half +ik+m][1]+= (-sfac_t[i*3+1][0]*sfac_t[j*3+1][1]+sfac_t[i*3+1][1]*sfac_t[j*3+1][0] - sfac[ size_half+ik+m][1])/(itime);
sfac[2*size_half+ik+m][0]+=( sfac_t[i*3+2][0]*sfac_t[j*3+2][0]+ sfac_t[i*3+2][1]*sfac_t[j*3+2][1]- sfac[2*size_half+ik+m][0])/(itime);
sfac[2*size_half+ik+m][1]+=(-sfac_t[i*3+2][0]*sfac_t[j*3+2][1]+sfac_t[i*3+2][1]*sfac_t[j*3+2][0]- sfac[2*size_half+ik+m][1])/(itime);
m++;
}
//std::cout<<" m= "<<m<<std::endl ;
}
}
}
}
//fare trasformata di fourier discreta (copia da spettrovibrazionale.cpp:84)
fftw_free( sfac_t);
fftw_execute_dft_c2r(fftw3,sfac,lista);
fftw_execute_dft_c2r(fftw3,&sfac[ size_half], &lista[size]);
fftw_execute_dft_c2r(fftw3,&sfac[2*size_half], &lista[2*size]);
}
CorrelatoreSpaziale::~CorrelatoreSpaziale(){
fftw_free(sfac);
}
CorrelatoreSpaziale & CorrelatoreSpaziale::operator =(const CorrelatoreSpaziale & destra) {
OperazioniSuLista<CorrelatoreSpaziale>::operator =(destra);
}
double CorrelatoreSpaziale::corr(unsigned int rx, unsigned int ry, unsigned int rz, unsigned int itype,unsigned int idim) {
if (itype >=tipi_atomi*(tipi_atomi+1)/2) {
std::cerr << " Errore: indice dei tipi fuori dal range in " __FILE__ "\n";
abort();
}
if (rx>=nk || ry>=nk || rz>=nk) {
std::cerr << " Errore: indice spaziale fuori dal range in " __FILE__ "\n";
abort();
}
if (idim >= 3) {
std::cerr << " Errore: indice della coordinata fuori dal range in " __FILE__ "\n";
}
return lista[idim*nk*nk*nk*tipi_atomi*(tipi_atomi+1)/2+ rx*nk*nk*tipi_atomi*(tipi_atomi+1)/2+
ry*nk*tipi_atomi*(tipi_atomi+1)/2+
rz*tipi_atomi*(tipi_atomi+1)/2+
itype
];
}