//gfactor3 is the random force \sqrt{\frac{2\gamma{}m_{i}}{\alpha*\delta{}t}}, \sqrt{12} makes the random array variance equal to unit.
gfactor3[i] = sqrt(2*fric_coef*atom->mass[i])*sqrt(force->mvv2e)*sqrt(12/h_timestep); //this still leaves a square energy term from the power spectrum H.
}
}
// generate random number array with zero mean and variance equal 1/12.
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
for (int m = 0; m < 2*N_f; m++) {
random_array_0[i][m] = random->uniform()-0.5;
random_array_1[i][m] = random->uniform()-0.5;
random_array_2[i][m] = random->uniform()-0.5;
}
}
// load omega_H with calculated spectrum at a specific temperature (corrected spectrum), omega_H is the Fourier transformation of time_H
for (int k = 0; k < 2*N_f; k++) {
double f_k=(k-N_f)/(2*N_f*h_timestep); //\omega_k=\frac{2\pi}{\delta{}h}\frac{k}{2N_f} for k from -N_f to N_f-1