ic_basic.hpp

//////////////////////////
// ic_basic.hpp
//////////////////////////
// 
// basic initial condition generator for gevolution
//
// Author: Julian Adamek (Université de Genève & Observatoire de Paris & Queen Mary University of London & Universität Zürich & ETH Zürich)
//
// Last modified: May 2026
//
//////////////////////////

#ifndef IC_BASIC_HEADER
#define IC_BASIC_HEADER

#include "prng_engine.hpp"
#include <gsl/gsl_spline.h>
#include "parser.hpp"
#include "tools.hpp"
#include <thrust/device_vector.h>
#include <nvtx3/nvToolsExt.h>

#ifndef Cplx
#define Cplx Imag
#endif  

// maximum number of characters per line in a Tk file
#ifndef MAX_LINESIZE
#define MAX_LINESIZE 2048
#endif

using namespace std;
using namespace LATfield2;

// should be larger than maximum Ngrid
#ifndef HUGE_SKIP
#define HUGE_SKIP   65536
#endif


//////////////////////////
// displace_pcls_ic_basic
//////////////////////////
// Description:
//   displaces particles according to gradient of displacement field
//   (accepts two displacement fields for baryon treatment = hybrid)
// 
// Arguments:
//   coeff             coefficient to be applied to displacement
//   lat_resolution    1 / Ngrid
//   part              particle to be displaced
//   ref_dist          distance vector (in lattice units) to reference lattice point
//   partInfo          particle metadata
//   fields            array of pointers to displacement fields
//   sites             array of respective sites
//   nfields           number of fields passed
//   params            additional parameters (ignored)
//   outputs           array for reduction variables; first entry will contain the displacement
//   noutputs          number of reduction variables
//
// Returns:
// 
//////////////////////////

__host__ __device__ void displace_pcls_ic_basic(double coeff, double lat_resolution, part_simple * part, double * ref_dist, part_simple_info partInfo, Field<Real> ** fields, Site * sites, int nfield, double * params, double * outputs, int noutputs)
{
	int i;
	Real gradxi[3] = {0, 0, 0};
	
	if (nfield > 1 && (*part).ID % 8 == 0)
		i = 1;
	else
		i = 0;
	
	gradxi[0] = (1.-ref_dist[1]) * (1.-ref_dist[2]) * ((*fields[i])(sites[i]+0) - (*fields[i])(sites[i]));
	gradxi[1] = (1.-ref_dist[0]) * (1.-ref_dist[2]) * ((*fields[i])(sites[i]+1) - (*fields[i])(sites[i]));
	gradxi[2] = (1.-ref_dist[0]) * (1.-ref_dist[1]) * ((*fields[i])(sites[i]+2) - (*fields[i])(sites[i]));
	gradxi[0] += ref_dist[1] * (1.-ref_dist[2]) * ((*fields[i])(sites[i]+1+0) - (*fields[i])(sites[i]+1));
	gradxi[1] += ref_dist[0] * (1.-ref_dist[2]) * ((*fields[i])(sites[i]+1+0) - (*fields[i])(sites[i]+0));
	gradxi[2] += ref_dist[0] * (1.-ref_dist[1]) * ((*fields[i])(sites[i]+2+0) - (*fields[i])(sites[i]+0));
	gradxi[0] += (1.-ref_dist[1]) * ref_dist[2] * ((*fields[i])(sites[i]+2+0) - (*fields[i])(sites[i]+2));
	gradxi[1] += (1.-ref_dist[0]) * ref_dist[2] * ((*fields[i])(sites[i]+2+1) - (*fields[i])(sites[i]+2));
	gradxi[2] += (1.-ref_dist[0]) * ref_dist[1] * ((*fields[i])(sites[i]+2+1) - (*fields[i])(sites[i]+1));
	gradxi[0] += ref_dist[1] * ref_dist[2] * ((*fields[i])(sites[i]+2+1+0) - (*fields[i])(sites[i]+2+1));
	gradxi[1] += ref_dist[0] * ref_dist[2] * ((*fields[i])(sites[i]+2+1+0) - (*fields[i])(sites[i]+2+0));
	gradxi[2] += ref_dist[0] * ref_dist[1] * ((*fields[i])(sites[i]+2+1+0) - (*fields[i])(sites[i]+1+0));
	
	gradxi[0] /= lat_resolution;
	gradxi[1] /= lat_resolution;
	gradxi[2] /= lat_resolution;
	
	if (noutputs > 0)
		*outputs = coeff * sqrt(gradxi[0]*gradxi[0] + gradxi[1]*gradxi[1] + gradxi[2]*gradxi[2]);

	for (i = 0; i < 3; i++) (*part).pos[i] += coeff*gradxi[i];
}

// callable struct for displace_pcls_ic_basic

struct displace_pcls_ic_basic_functor
{	
	__host__ __device__ void operator()(double dtau, double dx, part_simple * part, double * ref_dist, part_simple_info partInfo, Field<Real> * fields[], Site * sites, int nfield, double * params, double * outputs, int noutputs)
	{
		displace_pcls_ic_basic(dtau, dx, part, ref_dist, partInfo, fields, sites, nfield, params, outputs, noutputs);
	}
};


//////////////////////////
// initialize_q_ic_basic
//////////////////////////
// Description:
//   initializes velocities proportional to gradient of potential
//   (accepts two potentials for baryon treatment = hybrid)
// 
// Arguments:
//   coeff             coefficient to be applied to velocity
//   lat_resolution    1 / Ngrid
//   part              particle to be manipulated
//   ref_dist          distance vector (in lattice units) to reference lattice point
//   partInfo          particle metadata
//   fields            array of pointers to velocity potentials
//   sites             array of respective sites
//   nfields           number of fields passed
//   params            additional parameters (ignored)
//   outputs           array for reduction variables (ignored)
//   noutputs          number of reduction variables
//
// Returns: square of the velocity, (q/m)^2
// 
//////////////////////////

__host__ __device__ Real initialize_q_ic_basic(double coeff, double lat_resolution, part_simple * part, double * ref_dist, part_simple_info partInfo, Field<Real> ** fields, Site * sites, int nfield, double * params, double * outputs, int noutputs)
{
	int i;
	Real gradPhi[3] = {0, 0, 0};
	Real v2 = 0.;
	
	if (nfield > 1 && (*part).ID % 8 == 0)
		i = 1;
	else
		i = 0;
	
	gradPhi[0] = (1.-ref_dist[1]) * (1.-ref_dist[2]) * ((*fields[i])(sites[i]+0) - (*fields[i])(sites[i]));
	gradPhi[1] = (1.-ref_dist[0]) * (1.-ref_dist[2]) * ((*fields[i])(sites[i]+1) - (*fields[i])(sites[i]));
	gradPhi[2] = (1.-ref_dist[0]) * (1.-ref_dist[1]) * ((*fields[i])(sites[i]+2) - (*fields[i])(sites[i]));
	gradPhi[0] += ref_dist[1] * (1.-ref_dist[2]) * ((*fields[i])(sites[i]+1+0) - (*fields[i])(sites[i]+1));
	gradPhi[1] += ref_dist[0] * (1.-ref_dist[2]) * ((*fields[i])(sites[i]+1+0) - (*fields[i])(sites[i]+0));
	gradPhi[2] += ref_dist[0] * (1.-ref_dist[1]) * ((*fields[i])(sites[i]+2+0) - (*fields[i])(sites[i]+0));
	gradPhi[0] += (1.-ref_dist[1]) * ref_dist[2] * ((*fields[i])(sites[i]+2+0) - (*fields[i])(sites[i]+2));
	gradPhi[1] += (1.-ref_dist[0]) * ref_dist[2] * ((*fields[i])(sites[i]+2+1) - (*fields[i])(sites[i]+2));
	gradPhi[2] += (1.-ref_dist[0]) * ref_dist[1] * ((*fields[i])(sites[i]+2+1) - (*fields[i])(sites[i]+1));
	gradPhi[0] += ref_dist[1] * ref_dist[2] * ((*fields[i])(sites[i]+2+1+0) - (*fields[i])(sites[i]+2+1));
	gradPhi[1] += ref_dist[0] * ref_dist[2] * ((*fields[i])(sites[i]+2+1+0) - (*fields[i])(sites[i]+2+0));
	gradPhi[2] += ref_dist[0] * ref_dist[1] * ((*fields[i])(sites[i]+2+1+0) - (*fields[i])(sites[i]+1+0));
	
	gradPhi[0] /= lat_resolution;
	gradPhi[1] /= lat_resolution;
	gradPhi[2] /= lat_resolution;  
	
	for (i = 0 ; i < 3; i++)
	{
		(*part).vel[i] = -gradPhi[i] * coeff;
		v2 += (*part).vel[i] * (*part).vel[i];
	}
	
	return v2;
}

// callable struct for initialize_q_ic_basic

struct initialize_q_ic_basic_functor
{	
	__host__ __device__ Real operator()(double dtau, double dx, part_simple * part, double * ref_dist, part_simple_info partInfo, Field<Real> * fields[], Site * sites, int nfield, double * params, double * outputs, int noutputs)
	{
		return initialize_q_ic_basic(dtau, dx, part, ref_dist, partInfo, fields, sites, nfield, params, outputs, noutputs);
	}
};


//////////////////////////
// Pk_primordial
//////////////////////////
// Description:
//   power spectrum of the primordial curvature perturbation
// 
// Arguments:
//   k          wavenumber in units of inverse Mpc (!)
//   ic         settings for IC generation (datastructure containing parameters)
//
// Returns: amplitude of primordial power spectrum at given wavenumber
// 
//////////////////////////

inline double Pk_primordial(const double k, const icsettings & ic)
{
	return ic.A_s * pow(k / ic.k_pivot, ic.n_s - 1.);  // note that k_pivot is in units of inverse Mpc!
}


//////////////////////////
// loadHomogeneousTemplate
//////////////////////////
// Description:
//   loads a homogeneous template from a GADGET-2 file
// 
// Arguments:
//   filename   string containing the path to the template file
//   numpart    will contain the number of particles of the template
//   partdata   will contain the particle positions (memory will be allocated)
//
// Returns:
// 
//////////////////////////

void loadHomogeneousTemplate(const char * filename, long & numpart, float * & partdata)
{
	int i;

	if (parallel.grid_rank()[0] == 0) // read file
	{
		FILE * templatefile = NULL;
		int blocksize1 = 0, blocksize2 = 0, num_read = 0;
		gadget2_header filehdr;
		
		templatefile = fopen(filename, "r");
		
		if (templatefile == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Unable to open template file " << filename << "." << endl;
			parallel.abortForce();
		}
		
		if (fread(&blocksize1, sizeof(int), 1, templatefile) != 1 || blocksize1 != sizeof(filehdr))
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Unknown template file format - header not recognized." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		if (fread(&filehdr, sizeof(filehdr), 1, templatefile) != 1)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Unable to read header." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		if (fread(&blocksize2, sizeof(int), 1, templatefile) != 1 || blocksize1 != blocksize2)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Unknown template file format - block size mismatch while reading header." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		
		// analyze header for compatibility
		if (filehdr.num_files != 1)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Multiple input files (" << filehdr.num_files << ") currently not supported." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		if (filehdr.BoxSize <= 0.)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! BoxSize = " << filehdr.BoxSize << " not allowed." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		if (filehdr.npart[1] <= 0)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! No particles declared." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		
		partdata = (float *) malloc(3 * sizeof(float) * filehdr.npart[1]);
		if (partdata == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Memory error." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		
		num_read = fread(&blocksize1, sizeof(int), 1, templatefile);
		if (filehdr.npart[0] > 0)
		{
			if (fseek(templatefile, 3 * sizeof(float) * filehdr.npart[0], SEEK_CUR))
			{
				cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Unsuccesful attempt to skip gas particles (" << filehdr.npart[0] << ")." << endl;
				fclose(templatefile);
				parallel.abortForce();
			}
		}
		num_read = fread(partdata, sizeof(float), 3 * filehdr.npart[1], templatefile);
		if (num_read != (int) (3 * filehdr.npart[1]))
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Unable to read particle data." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		for (i = 2; i < 6; i++)
		{
			if (filehdr.npart[i] > 0)
			{
				if (fseek(templatefile, 3 * sizeof(float) * filehdr.npart[i], SEEK_CUR))
				{
					cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Unsuccesful attempt to skip particles (" << filehdr.npart[i] << ")." << endl;
					parallel.abortForce();
				}
			}
		}
		if (fread(&blocksize2, sizeof(int), 1, templatefile) != 1 || blocksize1 != blocksize2)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Unknown template file format - block size mismatch while reading particles." << endl;
			fclose(templatefile);
			parallel.abortForce();
		}
		
		fclose(templatefile);
		
		// reformat and check particle data
		for (i = 0; i < (int) (3 * filehdr.npart[1]); i++)
		{
			partdata[i] /= filehdr.BoxSize;
			if (partdata[i] < 0. || partdata[i] > 1.)
			{
				cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Particle data corrupted." << endl;
				parallel.abortForce();
			}
		}
		numpart = (long) filehdr.npart[1];
		
		parallel.broadcast_dim0<long>(numpart, 0);
	}
	else
	{
		parallel.broadcast_dim0<long>(numpart, 0);
		
		if (numpart <= 0)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Communication error." << endl;
			parallel.abortForce();
		}
		
		partdata = (float *) malloc(3 * sizeof(float) * numpart);
		if (partdata == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadHomogeneousTemplate! Memory error." << endl;
			parallel.abortForce();
		}
	}
	
	parallel.broadcast_dim0<float>(partdata, 3 * numpart, 0);
}


//////////////////////////
// loadPowerSpectrum
//////////////////////////
// Description:
//   loads a tabulated matter power spectrum from a file
// 
// Arguments:
//   filename   string containing the path to the template file
//   pkspline   will point to the gsl_spline which holds the tabulated
//              power spectrum (memory will be allocated)
//   boxsize    comoving box size (in the same units as used in the file)
//
// Returns:
// 
//////////////////////////

void loadPowerSpectrum(const char * filename, gsl_spline * & pkspline, const double boxsize)
{
	int i = 0, numpoints = 0;
	double * k;
	double * pk;

	if (parallel.grid_rank()[0] == 0) // read file
	{
		FILE * pkfile;
		char line[MAX_LINESIZE];
		double dummy1, dummy2;
		
		line[MAX_LINESIZE-1] = 0;
		
		pkfile = fopen(filename, "r");
		
		if (pkfile == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! Unable to open file " << filename << "." << endl;
			parallel.abortForce();
		}
		
		while (!feof(pkfile) && !ferror(pkfile))
		{
			if (fgets(line, MAX_LINESIZE, pkfile) != NULL && line[MAX_LINESIZE-1] != 0)
			{
				cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! Character limit (" << (MAX_LINESIZE-1) << "/line) exceeded in file " << filename << "." << endl;
				fclose(pkfile);
				parallel.abortForce();
			}
			
			if (sscanf(line, " %lf %lf", &dummy1, &dummy2) == 2 && !feof(pkfile) && !ferror(pkfile)) numpoints++;
		}
		
		if (numpoints < 2)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! No valid data found in file " << filename << "." << endl;
			fclose(pkfile);
			parallel.abortForce();
		}
		
		k = (double *) malloc(sizeof(double) * numpoints);
		pk = (double *) malloc(sizeof(double) * numpoints);
		
		if (k == NULL || pk == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! Memory error." << endl;
			fclose(pkfile);
			parallel.abortForce();
		}
		
		rewind(pkfile);
		
		while (!feof(pkfile) && !ferror(pkfile))
		{
			if (fgets(line, MAX_LINESIZE, pkfile) != NULL && sscanf(line, " %lf %lf", &dummy1, &dummy2) == 2 && !feof(pkfile) && !ferror(pkfile))
			{
				if (dummy1 < 0. || dummy2 < 0.)
				{
					cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! Negative entry encountered." << endl;
					free(k);
					free(pk);
					fclose(pkfile);
					parallel.abortForce();
				}
				
				if (i > 0)
				{
					if (k[i-1] >= dummy1 * boxsize)
					{
						cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! k-values are not strictly ordered." << endl;
						free(k);
						free(pk);
						fclose(pkfile);
						parallel.abortForce();
					}
				}
				
				k[i] = dummy1 * boxsize;
				pk[i] = sqrt(0.5 * dummy2 * boxsize);
				i++;
			}
		}
		
		fclose(pkfile);
		
		if (i != numpoints)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! File may have changed or file pointer corrupted." << endl;
			free(k);
			free(pk);
			parallel.abortForce();
		}
		
		parallel.broadcast_dim0<int>(numpoints, 0);
	}
	else
	{
		parallel.broadcast_dim0<int>(numpoints, 0);
		
		if (numpoints < 2)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! Communication error." << endl;
			parallel.abortForce();
		}
		
		k = (double *) malloc(sizeof(double) * numpoints);
		pk = (double *) malloc(sizeof(double) * numpoints);
		
		if (k == NULL || pk == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadPowerSpectrum! Memory error." << endl;
			parallel.abortForce();
		}
	}
	
	parallel.broadcast_dim0<double>(k, numpoints, 0);
	parallel.broadcast_dim0<double>(pk, numpoints, 0);
	
	pkspline = gsl_spline_alloc(gsl_interp_cspline, numpoints);
	
	gsl_spline_init(pkspline, k, pk, numpoints);
	
	free(k);
	free(pk);
}


//////////////////////////
// loadTransferFunctions (1)
//////////////////////////
// Description:
//   loads a set of tabulated transfer functions from a file
// 
// Arguments:
//   filename   string containing the path to the template file
//   tk_delta   will point to the gsl_spline which holds the tabulated
//              transfer function for delta (memory will be allocated)
//   tk_theta   will point to the gsl_spline which holds the tabulated
//              transfer function for theta (memory will be allocated)
//   qname      string containing the name of the component (e.g. "cdm")
//   boxsize    comoving box size (in the same units as used in the file)
//   h          conversion factor between 1/Mpc and h/Mpc (theta is in units of 1/Mpc)
//
// Returns:
// 
//////////////////////////

void loadTransferFunctions(const char * filename, gsl_spline * & tk_delta, gsl_spline * & tk_theta, const char * qname, const double boxsize, double h)
{
	int i = 0, numpoints = 0;
	double * k;
	double * tk_d;
	double * tk_t;

	if (parallel.grid_rank()[0] == 0) // read file
	{
		FILE * tkfile;
		char line[MAX_LINESIZE];
		char format[MAX_LINESIZE];
		char * ptr;
		double dummy[3];
		char dname[16];
		char tname[16];
		char kname[16];
		int kcol = -1, dcol = -1, tcol = -1, colmax;
		
		if (qname != NULL)
		{
			if (strcmp(qname, "vx") == 0)
			{
				sprintf(dname, "vx_smg");
				sprintf(tname, "vx_prime_smg");
				h = 1.;
			}
			else if (strcmp(qname, "eta") == 0)
			{
				sprintf(dname, "eta");
				sprintf(tname, "eta_prime");
				h /= boxsize;
			}
			else if (strcmp(qname, "h") == 0)
			{
				sprintf(dname, "h");
				sprintf(tname, "h_prime");
				h /= boxsize;
			}
			else if (strcmp(qname, "Nb") == 0)
			{
				sprintf(dname, "H_T_Nb_prime");
				sprintf(tname, "k2gamma_Nb");
				h = 1.;
			}
			else if (strcmp(qname, "d_m") == 0)
			{
				sprintf(dname, "d_cdm");
				sprintf(tname, "d_b");
				h = 1.;
			}
			else
			{
				sprintf(dname, "d_%s", qname);
				sprintf(tname, "t_%s", qname);
				h /= boxsize;
			}
    	}
		else
		{
			sprintf(dname, "phi");
			sprintf(tname, "psi");
			h = 1.;
		}
		sprintf(kname, "k (h/Mpc)");
		
		line[MAX_LINESIZE-1] = 0;
		
		tkfile = fopen(filename, "r");
		
		if (tkfile == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! Unable to open file " << filename << "." << endl;
			parallel.abortForce();
		}
		
		while (!feof(tkfile) && !ferror(tkfile))
		{
			if (fgets(line, MAX_LINESIZE, tkfile) != NULL && line[MAX_LINESIZE-1] != 0)
			{
				cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! Character limit (" << (MAX_LINESIZE-1) << "/line) exceeded in file " << filename << "." << endl;
				fclose(tkfile);
				parallel.abortForce();
			}
			
			if (line[0] != '#' && !feof(tkfile) && !ferror(tkfile)) numpoints++;
		}
		
		if (numpoints < 2)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! No valid data found in file " << filename << "." << endl;
			fclose(tkfile);
			parallel.abortForce();
		}
		
		k = (double *) malloc(sizeof(double) * numpoints);
		tk_d = (double *) malloc(sizeof(double) * numpoints);
		tk_t = (double *) malloc(sizeof(double) * numpoints);
		
		if (k == NULL || tk_d == NULL || tk_t == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! Memory error." << endl;
			fclose(tkfile);
			parallel.abortForce();
		}
		
		rewind(tkfile);
		
		while (!feof(tkfile) && !ferror(tkfile))
		{
			if (fgets(line, MAX_LINESIZE, tkfile) == NULL) break;

			for (ptr = line, i = 0; (ptr = strchr(ptr, ':')) != NULL; i++)
			{
				ptr++;
				if (strncmp(ptr, kname, strlen(kname)) == 0) kcol = i;
				else if (strncmp(ptr, tname, strlen(tname)) == 0) tcol = i;
				else if (strncmp(ptr, dname, strlen(dname)) == 0) dcol = i;
			}
			
			if (kcol >= 0 && dcol >= 0 && tcol >= 0) break;
		}
		
		if (kcol < 0 || dcol < 0 || tcol < 0)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! Unable to identify requested columns! (" << kname << ": " << kcol << ", " << dname << ": " << dcol << ", " << tname << ": " << tcol << ")" << endl;
			fclose(tkfile);
			free(k);
			free(tk_d);
			free(tk_t);
			parallel.abortForce();
		}
		
		colmax = i;
		for (i = 0, ptr=format; i < colmax; i++)
		{
		    if (i == kcol || i == dcol || i == tcol)
		    {
		        strncpy(ptr, " %lf", 4);
		        ptr += 4;
		    }
		    else
		    {
		        strncpy(ptr, " %*lf", 5);
		        ptr += 5;
		    }
		}
		*ptr = '\0';
		
		if (kcol < dcol && dcol < tcol)
		{
			kcol = 0; dcol = 1; tcol = 2;
		}
		else if (kcol < tcol && tcol < dcol)
		{
			kcol = 0; dcol = 2; tcol = 1;
		}
		else if (dcol < kcol && kcol < tcol)
		{
			kcol = 1; dcol = 0; tcol = 2;
		}
		else if (dcol < tcol && tcol < kcol)
		{
			kcol = 2; dcol = 0; tcol = 1;
		}
		else if (tcol < kcol && kcol < dcol)
		{
			kcol = 1; dcol = 2; tcol = 0;
		}
		else if (tcol < dcol && dcol < kcol)
		{
			kcol = 2; dcol = 1; tcol = 0;
		}
		else
		{
			cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! Inconsistent columns!" << endl;
			fclose(tkfile);
			free(k);
			free(tk_d);
			free(tk_t);
			parallel.abortForce();
		}
		
		i = 0;
		while (!feof(tkfile) && !ferror(tkfile))
		{
			if (fgets(line, MAX_LINESIZE, tkfile) == NULL) break;
			
			if (sscanf(line, format, dummy, dummy+1, dummy+2) == 3 && !feof(tkfile) && !ferror(tkfile))
			{
				if (dummy[kcol] < 0.)
				{
					cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! Negative k-value encountered." << endl;
					free(k);
					free(tk_d);
					free(tk_t);
					fclose(tkfile);
					parallel.abortForce();
				}
				
				if (i > 0)
				{
					if (k[i-1] >= dummy[kcol] * boxsize)
					{
						cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! k-values are not strictly ordered." << endl;
						free(k);
						free(tk_d);
						free(tk_t);
						fclose(tkfile);
						parallel.abortForce();
					}
				}
				
				k[i] = dummy[kcol] * boxsize;
				tk_d[i] = dummy[dcol];
				tk_t[i] = dummy[tcol] / h;
				i++;
			}
		}
		
		fclose(tkfile);
		
		if (i != numpoints)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! File may have changed or file pointer corrupted." << endl;
			free(k);
			free(tk_d);
			free(tk_t);
			parallel.abortForce();
		}
		
		parallel.broadcast<int>(numpoints, 0);
	}
	else
	{
		parallel.broadcast<int>(numpoints, 0);
		
		if (numpoints < 2)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! Communication error." << endl;
			parallel.abortForce();
		}
		
		k = (double *) malloc(sizeof(double) * numpoints);
		tk_d = (double *) malloc(sizeof(double) * numpoints);
		tk_t = (double *) malloc(sizeof(double) * numpoints);
		
		if (k == NULL || tk_d == NULL || tk_t == NULL)
		{
			cerr << " proc#" << parallel.rank() << ": error in loadTransferFunctions! Memory error." << endl;
			parallel.abortForce();
		}
	}
	
	parallel.broadcast_dim0<double>(k, numpoints, 0);
	parallel.broadcast_dim0<double>(tk_d, numpoints, 0);
	parallel.broadcast_dim0<double>(tk_t, numpoints, 0);
	
	tk_delta = gsl_spline_alloc(gsl_interp_cspline, numpoints);
	tk_theta = gsl_spline_alloc(gsl_interp_cspline, numpoints);
	
	gsl_spline_init(tk_delta, k, tk_d, numpoints);
	gsl_spline_init(tk_theta, k, tk_t, numpoints);
	
	free(k);
	free(tk_d);
	free(tk_t);
}


//////////////////////////
// generateCICKernel
//////////////////////////
// Description:
//   generates convolution kernel for CIC projection
// 
// Arguments:
//   ker        reference to allocated field that will contain the convolution kernel
//   numpcl     number of particles in the pcldata pointer
//   pcldata    raw particle data from which the kernel will be constructed;
//              a standard kernel will be provided in case no particles are specified
//   numtile    tiling factor used for the particle template
//
// Returns:
// 
//////////////////////////

void generateCICKernel(Field<Real> & ker, const long numpcl = 0, float * pcldata = NULL, const int numtile = 1)
{
	const long linesize = ker.lattice().sizeLocal(0);
	Real renorm = linesize * linesize;
	long sx, sy, sz;
	float wx, wy, wz, q1, ww;
	Site x(ker.lattice());

	thrust::fill_n(thrust::device, ker.data(), ker.lattice().sitesLocalGross(), Real(0));
	
	if (numpcl == 0 || pcldata == NULL) // standard kernel
	{
		if (x.setCoord(0, 0, 0))
		{
			ker(x) = Real(6) * renorm;
			ker(x+0) = -renorm;
			x.setCoord(linesize-1, 0, 0);
			ker(x) = -renorm;
		}
	
		if (x.setCoord(0, 1, 0))
			ker(x) = -renorm;
	
		if (x.setCoord(0, 0, 1))
			ker(x) = -renorm;
	
		if (x.setCoord(0, linesize-1, 0))
			ker(x) = -renorm;
	
		if (x.setCoord(0, 0, linesize-1))
			ker(x) = -renorm;
		
		return;
	}
	
	// compute kernel explicitly
	
	renorm /= (Real) (numpcl * (long) numtile * (long) numtile * (long) numtile) / (Real) (linesize * linesize * linesize);
	
	for (long i = 0; i < numpcl; i++)
	{
		for (int oct = 0; oct < 8; oct++)
		{
			if (oct == 0)
			{
				if (linesize * pcldata[3*i] / numtile >= 1. || linesize * pcldata[3*i+1] / numtile >= 1. || linesize * pcldata[3*i+2] / numtile >= 1.)
					continue;
				
				// particle is in the first octant   
				wx = 1.f - linesize * pcldata[3*i] / numtile;
				wy = 1.f - linesize * pcldata[3*i+1] / numtile;
				wz = 1.f - linesize * pcldata[3*i+2] / numtile;
				
				sx = 1;
				sy = 1;
				sz = 1;
			}
			else if (oct == 1)
			{
				if (linesize * (1.f-pcldata[3*i]) / numtile >= 1.f || linesize * pcldata[3*i+1] / numtile >= 1.f || linesize * pcldata[3*i+2] / numtile >= 1.f)
					continue;
					
				// particle is in the second octant   
				wx = 1.f - linesize * (1.f-pcldata[3*i]) / numtile;
				wy = 1.f - linesize * pcldata[3*i+1] / numtile;
				wz = 1.f - linesize * pcldata[3*i+2] / numtile;
				
				sx = -1;
				sy = 1;
				sz = 1;
			}
			else if (oct == 2)
			{
				if (linesize * (1.f-pcldata[3*i]) / numtile >= 1.f || linesize * (1.f-pcldata[3*i+1]) / numtile >= 1.f || linesize * pcldata[3*i+2] / numtile >= 1.f)
					continue;
					
				// particle is in the third octant   
				wx = 1.f - linesize * (1.f-pcldata[3*i]) / numtile;
				wy = 1.f - linesize * (1.f-pcldata[3*i+1]) / numtile;
				wz = 1.f - linesize * pcldata[3*i+2] / numtile;

				sx = -1;
				sy = -1;
				sz = 1;
			}
			else if (oct == 3)
			{
				if (linesize * pcldata[3*i] / numtile >= 1.f || linesize * (1.f-pcldata[3*i+1]) / numtile >= 1.f || linesize * pcldata[3*i+2] / numtile >= 1.f)
					continue;
					
				// particle is in the fourth octant   
				wx = 1.f - linesize * pcldata[3*i] / numtile;
				wy = 1.f - linesize * (1.f-pcldata[3*i+1]) / numtile;
				wz = 1.f - linesize * pcldata[3*i+2] / numtile;

				sx = 1;
				sy = -1;
				sz = 1;
			}
			else if (oct == 4)
			{
				if (linesize * pcldata[3*i] / numtile >= 1.f || linesize * pcldata[3*i+1] / numtile >= 1.f || linesize * (1.f-pcldata[3*i+2]) / numtile >= 1.f)
					continue;

				// particle is in the fifth octant
				wx = 1.f - linesize * pcldata[3*i] / numtile;
				wy = 1.f - linesize * pcldata[3*i+1] / numtile;
				wz = 1.f - linesize * (1.f-pcldata[3*i+2]) / numtile;

				sx = 1;
				sy = 1;
				sz = -1;
			}
			else if (oct == 5)
			{
				if (linesize * (1.f-pcldata[3*i]) / numtile >= 1.f || linesize * pcldata[3*i+1] / numtile >= 1.f || linesize * (1.f-pcldata[3*i+2]) / numtile >= 1.f)
					continue;

				// particle is in the sixth octant
				wx = 1.f - linesize * (1.f-pcldata[3*i]) / numtile;
				wy = 1.f - linesize * pcldata[3*i+1] / numtile;
				wz = 1.f - linesize * (1.f-pcldata[3*i+2]) / numtile;

				sx = -1;
				sy = 1;
				sz = -1;
			}
			else if (oct == 6)
			{
				if (linesize * (1.f-pcldata[3*i]) / numtile >= 1.f || linesize * (1.f-pcldata[3*i+1]) / numtile >= 1.f || linesize * (1.f-pcldata[3*i+2]) / numtile >= 1.f)
					continue;

				// particle is in the seventh octant
				wx = 1.f - linesize * (1.f-pcldata[3*i]) / numtile;
				wy = 1.f - linesize * (1.f-pcldata[3*i+1]) / numtile;
				wz = 1.f - linesize * (1.f-pcldata[3*i+2]) / numtile;

				sx = -1;
				sy = -1;
				sz = -1;
			}
			else if (oct == 7)
			{
				if (linesize * pcldata[3*i] / numtile >= 1.f || linesize * (1.f-pcldata[3*i+1]) / numtile >= 1.f || linesize * (1.f-pcldata[3*i+2]) / numtile >= 1.f)
					continue;

				// particle is in the eight-th octant
				wx = 1.f - linesize * pcldata[3*i] / numtile;
				wy = 1.f - linesize * (1.f-pcldata[3*i+1]) / numtile;
				wz = 1.f - linesize * (1.f-pcldata[3*i+2]) / numtile;

				sx = 1;
				sy = -1;
				sz = -1;
			}
			else
				continue;
			
			// 0-direction
			
			ww = wy*wz*renorm;		   
			q1 = (wx > 0.9f) ? 2.f : 1.f;
			
			if (x.setCoord(0, 0, 0))
			{
				ker(x) += ww * wy * wz * q1;
				x.setCoord((linesize+sx)%linesize, 0, 0);
				ker(x) -= ww * wy * wz;
				if (wx > 0.9f)
				{
					x.setCoord((linesize-sx)%linesize, 0, 0);
					ker(x) -= ww * wy * wz;
				}
			}
			
			if (x.setCoord(0, 0, (linesize+sz)%linesize))
			{
				ker(x) += ww * wy * (1.f-wz) * q1;
				x.setCoord((linesize+sx)%linesize, 0, (linesize+sz)%linesize);
				ker(x) -= ww * wy * (1.f-wz);
				if (wx > 0.9f)
				{
					x.setCoord((linesize-sx)%linesize, 0, (linesize+sz)%linesize);
					ker(x) -= ww * wy * (1.f-wz);
				}
			}
			
			if (x.setCoord(0, (linesize+sy)%linesize, 0))
			{
				ker(x) += ww * (1.f-wy) * wz * q1;
				x.setCoord((linesize+sx)%linesize, (linesize+sy)%linesize, 0);
				ker(x) -= ww * (1.f-wy) * wz;
				if (wx > 0.9)
				{
					x.setCoord((linesize-sx)%linesize, (linesize+sy)%linesize, 0);
					ker(x) -= ww * (1.f-wy) * wz;
				}
			}
			
			if (x.setCoord(0, (linesize+sy)%linesize, (linesize+sz)%linesize))
			{
				ker(x) += ww * (1.f-wy) * (1.f-wz) * q1;
				x.setCoord((linesize+sx)%linesize, (linesize+sy)%linesize, (linesize+sz)%linesize);
				ker(x) -= ww * (1.f-wy) * (1.f-wz);
				if (wx > 0.9f)
				{
					x.setCoord((linesize-sx)%linesize, (linesize+sy)%linesize, (linesize+sz)%linesize);
					ker(x) -= ww * (1.f-wy) * (1.f-wz);
				}
			}
			
			// 1-direction
			
			ww = wx*wz*renorm;
			q1 = (wy > 0.9f) ? 2.f : 1.f;
			
			if (x.setCoord(0, 0, 0))
			{
				ker(x) += ww * wx * wz * q1;
				x.setCoord((linesize+sx)%linesize, 0, 0);
				ker(x) += ww * (1.f-wx) * wz * q1;
			}
			
			if (x.setCoord(0, 0, (linesize+sz)%linesize))
			{
				ker(x) += ww * wx * (1.f-wz) * q1;
				x.setCoord((linesize+sx)%linesize, 0, (linesize+sz)%linesize);
				ker(x) += ww * (1.f-wx) * (1.f-wz) * q1;
			}
			
			if (x.setCoord(0, (linesize+sy)%linesize, 0))
			{
				ker(x) -= ww * wx * wz;
				x.setCoord((linesize+sx)%linesize, (linesize+sy)%linesize, 0);
				ker(x) -= ww * (1.f-wx) * wz;
			}
			
			if (x.setCoord(0, (linesize+sy)%linesize, (linesize+sz)%linesize))
			{
				ker(x) -= ww * wx * (1.f-wz);
				x.setCoord((linesize+sx)%linesize, (linesize+sy)%linesize, (linesize+sz)%linesize);
				ker(x) -= ww * (1.f-wx) * (1.f-wz);
			}
			
			if (wy > 0.9f)
			{
				if (x.setCoord(0, (linesize-sy)%linesize, 0))
				{
					ker(x) -= ww * wx * wz;
					x.setCoord((linesize+sx)%linesize, (linesize-sy)%linesize, 0);
					ker(x) -= ww * (1.f-wx) * wz;
				}
				
				if (x.setCoord(0, (linesize-sy)%linesize, (linesize+sz)%linesize))
				{
					ker(x) -= ww * wx * (1.f-wz);
					x.setCoord((linesize+sx)%linesize, (linesize-sy)%linesize, (linesize+sz)%linesize);
					ker(x) -= ww * (1.f-wx) * (1.f-wz);
				}
			}
						
			// 2-direction
			
			ww = wx*wy*renorm;
			q1 = (wz > 0.9f) ? 2.f : 1.f;
			
			if (x.setCoord(0, 0, 0))
			{
				ker(x) += ww * wx * wy * q1;
				x.setCoord((linesize+sx)%linesize, 0, 0);
				ker(x) += ww * (1.f-wx) * wy * q1;
			}
			
			if (x.setCoord(0, (linesize+sy)%linesize, 0))
			{
				ker(x) += ww * wx * (1.f-wy) * q1;
				x.setCoord((linesize+sx)%linesize, (linesize+sy)%linesize, 0);
				ker(x) += ww * (1.f-wx) * (1.f-wy) * q1;
			}
			
			if (x.setCoord(0, 0, (linesize+sz)%linesize))
			{
				ker(x) -= ww * wx * wy;
				x.setCoord((linesize+sx)%linesize, 0, (linesize+sz)%linesize);
				ker(x) -= ww * (1.f-wx) * wy;
			}
			
			if (x.setCoord(0, (linesize+sy)%linesize, (linesize+sz)%linesize))
			{
				ker(x) -= ww * wx * (1.f-wy);
				x.setCoord((linesize+sx)%linesize, (linesize+sy)%linesize, (linesize+sz)%linesize);
				ker(x) -= ww * (1.f-wx) * (1.f-wy);
			}
			
			if (wz > 0.9f)
			{
				if (x.setCoord(0, 0, (linesize-sz)%linesize))
				{
					ker(x) -= ww * wx * wy;
					x.setCoord((linesize+sx)%linesize, 0, (linesize-sz)%linesize);
					ker(x) -= ww * (1.f-wx) * wy;
				}
				
				if (x.setCoord(0, (linesize+sy)%linesize, (linesize-sz)%linesize))
				{
					ker(x) -= ww * wx * (1.f-wy);
					x.setCoord((linesize+sx)%linesize, (linesize+sy)%linesize, (linesize-sz)%linesize);
					ker(x) -= ww * (1.f-wx) * (1.f-wy);
				}
			}
		}
	}
}


#ifdef FFT3D

inline float boxMullerRadialAmplitude(const float uniform)
{
#ifdef FIXED_ICS
	(void) uniform;
	return sqrt(2.0f);
#else
	return sqrt(-2.0f * log(uniform));
#endif
}

//////////////////////////
// generateDisplacementField (generateRealization)
//////////////////////////
// Description:
//   generates particle displacement field
//
// Non-type template parameters:
//   ignorekernel  this is effectively an optimization flag defaulted to 0; instantiating with 1 instead will cause
//                 the function to ignore the convolution kernel, allowing the function to be used for generating
//                 realizations (generateRealization is simply an alias for generateDisplacementField<1>)
// 
// Arguments:
//   potFT         reference to allocated field that contains the convolution kernel relating the potential
//                 (generating the displacement field) with the bare density perturbation; will contain the
//                 Fourier image of the potential generating the displacement field
//   coeff         gauge correction coefficient "H_conformal^2"
//   pkspline      pointer to a gsl_spline which holds a tabulated power spectrum
//   seed          initial seed for random number generator
//   ksphere       flag to indicate that only a sphere in k-space should be initialized
//                 (default = 0: full k-space cube is initialized)
//   deconvolve_f  flag to indicate deconvolution function
//                 0: no deconvolution
//                 1: sinc (default)
//
// Returns:
// 
//////////////////////////

#ifndef generateRealization
#define generateRealization generateDisplacementField<1>
#endif

template<int ignorekernel = 0>
void generateDisplacementField(Field<Cplx> & potFT, const Real coeff, const gsl_spline * pkspline, const unsigned int seed, const int ksphere = 0, const int deconvolve_f = 1)
{
	const int linesize = potFT.lattice().size(1);
	const int kmax = (linesize / 2) - 1;
	const int kmax2 = kmax * kmax;
	rKSite k(potFT.lattice());
	int kx;
	int kymin, kymax, kzmin, kzmax, kyend, kzend;
	float r1, r2, k2, s;
	float * sinc;
	sitmo::prng_engine prng;
	uint64_t huge_skip = HUGE_SKIP;
	gsl_interp_accel * acc;

	if constexpr (ignorekernel == 0)
	{
		nvtxRangePushA("generateDisplacementField");
	}
	else
	{
		nvtxRangePushA("generateRealization");
	}
	
	if (linesize <= STACK_ALLOCATION_LIMIT)
		sinc = (float *) alloca(linesize * sizeof(float));
	else
		sinc = (float *) malloc(linesize * sizeof(float));
	
	sinc[0] = 1.;
	if (deconvolve_f == 1)
	{
		for (int i = 1; i < linesize; i++)
			sinc[i] = sin(M_PI * (float) i / (float) linesize) * (float) linesize / (M_PI * (float) i);
	}
	else
	{
		for (int i = 1; i < linesize; i++)
			sinc[i] = 1.0f;
	}
	
	k.initialize(potFT.lattice(), potFT.lattice().siteLast());
	kymax = k.coord(1);
	kzmax = k.coord(2);
	k.initialize(potFT.lattice(), potFT.lattice().siteFirst());
	kymin = k.coord(1);
	kzmin = k.coord(2);

	kzend = (linesize / 2 < kzmax) ? linesize / 2 : kzmax;
		
	if (kymin < (linesize / 2) + 1 && kzmin < (linesize / 2) + 1)
	{
		kyend = (linesize / 2 < kymax) ? linesize / 2 : kymax;

		#pragma omp parallel default(shared) private(acc, kx, k2, r1, r2, s) firstprivate(prng, k)
		{
		acc = gsl_interp_accel_alloc();
		
		#pragma omp for collapse(2)
		for (int kz = kzmin; kz <= kzend; kz++)
		{
			for (int ky = kymin; ky <= kyend; ky++)
			{
				prng.seed(seed);
		   
				if (ky > 0 || kz > 0)
				{
					prng.discard(((uint64_t) kz * huge_skip + (uint64_t) ky) * huge_skip);
				}
				//i = 0;
				for (kx = 0; kx < (linesize / 2) + 1; kx++)
				{
					k.setCoord(kx, ky, kz);
					
					k2 = (float) (kx * kx) + (float) (ky * ky) + (float) (kz * kz);
					
					if (k2 == 0 || kx >= kmax || ky >= kmax || kz >= kmax || (k2 >= kmax2 && ksphere > 0))
					{
						potFT(k) = Cplx(0, 0);
					}
					else
					{
						s = sinc[kx] * sinc[ky] * sinc[kz];
						k2 *= 4.0f * M_PI * M_PI;
						do
						{
							r1 = (float) prng() / (float) sitmo::prng_engine::max();
							//i++;
						}
						while (r1 == 0);

						r2 = (float) prng() / (float) sitmo::prng_engine::max();
						r1 = boxMullerRadialAmplitude(r1);
						//i++;
					
						if constexpr (ignorekernel == 0)
						{
							potFT(k) = Cplx(cos(2.0f * M_PI * r2), sin(2.0f * M_PI * r2)) * (1.0f + 7.5f * static_cast<float>(coeff) / k2) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s / potFT(k);
						}
						else
						{
							potFT(k) = Cplx(cos(2.0f * M_PI * r2), sin(2.0f * M_PI * r2)) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s;
						}
					}
				}
				//prng.discard(huge_skip - (uint64_t) i);
			}
			//prng.discard(huge_skip * (huge_skip - (uint64_t) (kyend - kymin + 1)));
		}

		gsl_interp_accel_free(acc);
		}
	}
	
	if (kymax >= (linesize / 2) + 1 && kzmin < (linesize / 2) + 1)
	{
		//prng.seed(seed);
		//prng.discard(((huge_skip + (uint64_t) kzmin) * huge_skip + (uint64_t) (linesize - kymax)) * huge_skip);

		kyend = (linesize / 2 >= kymin) ? (linesize / 2) + 1 : kymin;

		#pragma omp parallel default(shared) private(acc, kx, k2, r1, r2, s) firstprivate(prng, k)
		{
		acc = gsl_interp_accel_alloc();
		
		#pragma omp for collapse(2)
		for (int kz = kzmin; kz <= kzend; kz++)
		{
			for (int ky = kymax; ky >= kyend; ky--)
			{
				prng.seed(seed);
				prng.discard(((huge_skip + (uint64_t) kz) * huge_skip + (uint64_t) (linesize - ky)) * huge_skip);

				//i = 0;
				for (kx = 0; kx < (linesize / 2) + 1; kx++)
				{
					k.setCoord(kx, ky, kz);
										
					k2 = (float) (kx * kx) + (float) ((linesize-ky) * (linesize-ky)) + (float) (kz * kz);
					
					if (kx >= kmax || (linesize-ky) >= kmax || kz >= kmax || (k2 >= kmax2 && ksphere > 0))
					{
						potFT(k) = Cplx(0, 0);
					}
					else
					{
						s = sinc[kx] * sinc[linesize-ky] * sinc[kz];
						k2 *= 4.0f * M_PI * M_PI;
						do
						{
							r1 = (float) prng() / (float) sitmo::prng_engine::max();
							//i++;
						}
						while (r1 == 0);
						r2 = (float) prng() / (float) sitmo::prng_engine::max();
						r1 = boxMullerRadialAmplitude(r1);
						if (kx == 0 && kz == 0)
						{
							// Conjugate (0, -ky, 0), retaining the original stream
							// advancement above so kx > 0 modes remain unchanged.
							sitmo::prng_engine conjugate_prng;
							conjugate_prng.seed(seed);
							conjugate_prng.discard((uint64_t) (linesize - ky) * huge_skip);
							do
							{
								r1 = (float) conjugate_prng() / (float) sitmo::prng_engine::max();
							}
							while (r1 == 0);
							r2 = -(float) conjugate_prng() / (float) sitmo::prng_engine::max();
							r1 = boxMullerRadialAmplitude(r1);
						}

						//i++;
						
						if constexpr (ignorekernel == 0)
						{
							potFT(k) = Cplx(cos(2.0f * M_PI * r2), sin(2.0f * M_PI * r2)) * (1.0f + 7.5f * static_cast<float>(coeff) / k2) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s / potFT(k);
						}
						else
						{
							potFT(k) = Cplx(cos(2.0f * M_PI * r2), sin(2.0f * M_PI * r2)) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s;
						}
					}
				}
				//prng.discard(huge_skip - (uint64_t) i);
			}
			//prng.discard(huge_skip * (huge_skip - (uint64_t) (kymax - kyend + 1)));
		}

		gsl_interp_accel_free(acc);
		}
	}

	kzend = (linesize / 2 >= kzmin) ? (linesize / 2) + 1 : kzmin;
	
	if (kymin < (linesize / 2) + 1 && kzmax >= (linesize / 2) + 1)
	{
		//prng.seed(seed);
		//prng.discard(((huge_skip + huge_skip + (uint64_t) (linesize - kzmax)) * huge_skip + (uint64_t) kymin) * huge_skip);

		kyend = (linesize / 2 < kymax) ? linesize / 2 : kymax;

		#pragma omp parallel default(shared) private(acc, kx, k2, r1, r2, s) firstprivate(prng, k)
		{
		acc = gsl_interp_accel_alloc();

		#pragma omp for collapse(2)
		for (int kz = kzmax; kz >= kzend; kz--)
		{
			for (int ky = kymin; ky <= kyend; ky++)
			{
				prng.seed(seed);
				prng.discard(((huge_skip + huge_skip + (uint64_t) (linesize - kz)) * huge_skip + (uint64_t) ky) * huge_skip);

				//i = 0;
				for (kx = 1; kx < (linesize / 2) + 1; kx++)
				{
					k.setCoord(kx, ky, kz);
					
					k2 = (float) (kx * kx) + (float) (ky * ky) + (float) ((linesize-kz) * (linesize-kz));
					
					if (kx >= kmax || ky >= kmax || (linesize-kz) >= kmax || (k2 >= kmax2 && ksphere > 0))
					{
						potFT(k) = Cplx(0, 0);
					}
					else
					{
						s = sinc[kx] * sinc[ky] * sinc[linesize-kz];
						k2 *= 4.0f * M_PI * M_PI;
						do
						{
							r1 = (float) prng() / (float) sitmo::prng_engine::max();
							//i++;
						}
						while (r1 == 0);
						r2 = (float) prng() / (float) sitmo::prng_engine::max();
						r1 = boxMullerRadialAmplitude(r1);
						//i++;
					
						if constexpr (ignorekernel == 0)
						{
							potFT(k) = Cplx(cos(2.0f * M_PI * r2), sin(2.0f * M_PI * r2)) * (1.0f + 7.5f * static_cast<float>(coeff) / k2) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s / potFT(k);
						}
						else
						{
							potFT(k) = Cplx(cos(2.0f * M_PI * r2), sin(2.0f * M_PI * r2)) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s;
						}
					}
				}
				//prng.discard(huge_skip - (uint64_t) i);
			}
			//prng.discard(huge_skip * (huge_skip - (uint64_t) (kyend - kymin + 1)));
		}
		
		//prng.seed(seed);
		//prng.discard(((uint64_t) (linesize - kzmax) * huge_skip + (uint64_t) kymin) * huge_skip);
		kx = 0;
		
		#pragma omp for
		for (int kz = kzmax; kz >= kzend; kz--)
		{
			for (int ky = kymin; ky <= kyend; ky++)
			{
				prng.seed(seed);
				// Match (0, -ky, -kz); ky == 0 belongs to the first quadrant stream.
				prng.discard((((ky == 0 ? 0ULL : huge_skip) + (uint64_t) (linesize - kz)) * huge_skip + (uint64_t) ky) * huge_skip);

				k.setCoord(kx, ky, kz);
					
				k2 = (float) (ky * ky) + (float) ((linesize-kz) * (linesize-kz));
				//i = 0;
				
				if (ky >= kmax || (linesize-kz) >= kmax || (k2 >= kmax2 && ksphere > 0))
				{
					potFT(k) = Cplx(0, 0);
				}
				else
				{
					s = sinc[ky] * sinc[linesize-kz];
					k2 *= 4.0f * M_PI * M_PI;
					do
					{
						r1 = (float) prng() / (float) sitmo::prng_engine::max();
						//i++;
					}
					while (r1 == 0);
					r2 = (float) prng() / (float) sitmo::prng_engine::max();
					r1 = boxMullerRadialAmplitude(r1);
					//i++;
				
					if constexpr (ignorekernel == 0)
					{
						potFT(k) = Cplx(cos(2.0f * M_PI * r2), -sin(2.0f * M_PI * r2)) * (1.0f + 7.5f * static_cast<float>(coeff) / k2) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s / potFT(k);
					}
					else
					{
						potFT(k) = Cplx(cos(2.0f * M_PI * r2), -sin(2.0f * M_PI * r2)) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s;
					}
				}
								
				//prng.discard(huge_skip - (uint64_t) i);
			}
			//prng.discard(huge_skip * (huge_skip - (uint64_t) (kyend - kymin + 1)));
		}

		gsl_interp_accel_free(acc);
		}
	}
	
	if (kymax >= (linesize / 2) + 1 && kzmax >= (linesize / 2) + 1)
	{
		//prng.seed(seed);
		//prng.discard(((huge_skip + huge_skip + huge_skip + (uint64_t) (linesize - kzmax)) * huge_skip + (uint64_t) (linesize - kymax)) * huge_skip);

		kyend = (linesize / 2 >= kymin) ? (linesize / 2) + 1 : kymin;
		
		#pragma omp parallel default(shared) private(acc, kx, k2, r1, r2, s) firstprivate(prng, k)
		{
		acc = gsl_interp_accel_alloc();

		#pragma omp for collapse(2)
		for (int kz = kzmax; kz >= kzend; kz--)
		{
			for (int ky = kymax; ky >= kyend; ky--)
			{
				prng.seed(seed);
				prng.discard(((huge_skip + huge_skip + huge_skip + (uint64_t) (linesize - kz)) * huge_skip + (uint64_t) (linesize - ky)) * huge_skip);

				//i = 0;
				for (kx = 1; kx < (linesize / 2) + 1; kx++)
				{
					k.setCoord(kx, ky, kz);
					
					k2 = (float) (kx * kx) + (float) ((linesize-ky) * (linesize-ky)) + (float) ((linesize-kz) * (linesize-kz));
					
					if (kx >= kmax || (linesize-ky) >= kmax || (linesize-kz) >= kmax || (k2 >= kmax2 && ksphere > 0))
					{
						potFT(k) = Cplx(0, 0);
					}
					else
					{
						s = sinc[kx] * sinc[linesize-ky] * sinc[linesize-kz];
						k2 *= 4.0f * M_PI * M_PI;
						do
						{
							r1 = (float) prng() / (float) sitmo::prng_engine::max();
							//i++;
						}
						while (r1 == 0);
						r2 = (float) prng() / (float) sitmo::prng_engine::max();
						r1 = boxMullerRadialAmplitude(r1);
						//i++;
					
						if constexpr (ignorekernel == 0)
						{
							potFT(k) = Cplx(cos(2.0f * M_PI * r2), sin(2.0f * M_PI * r2)) * (1.0f + 7.5f * static_cast<float>(coeff) / k2) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s / potFT(k);
						}
						else
						{
							potFT(k) = Cplx(cos(2.0f * M_PI * r2), sin(2.0f * M_PI * r2)) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s;
						}
					}
				}
				//prng.discard(huge_skip - (uint64_t) i);
			}
			//prng.discard(huge_skip * (huge_skip - (uint64_t) (kymax - kyend + 1)));
		}
		
		//prng.seed(seed);
		//prng.discard(((huge_skip + huge_skip + (uint64_t) (linesize - kzmax)) * huge_skip + (uint64_t) (linesize - kymax)) * huge_skip);
		kx = 0;
		
		#pragma omp for collapse(2)
		for (int kz = kzmax; kz >= kzend; kz--)
		{
			for (int ky = kymax; ky >= kyend; ky--)
			{
				prng.seed(seed);
				// Match the first quadrant mode (0, -ky, -kz).
				prng.discard(((uint64_t) (linesize - kz) * huge_skip + (uint64_t) (linesize - ky)) * huge_skip);

				k.setCoord(kx, ky, kz);
					
				k2 = (float) ((linesize-ky) * (linesize-ky)) + (float) ((linesize-kz) * (linesize-kz));
				//i = 0;
				
				if ((linesize-ky) >= kmax || (linesize-kz) >= kmax || (k2 >= kmax2 && ksphere > 0))
				{
					potFT(k) = Cplx(0, 0);
				}
				else
				{
					s = sinc[linesize-ky] * sinc[linesize-kz];
					k2 *= 4.0f * M_PI * M_PI;
					do
					{
						r1 = (float) prng() / (float) sitmo::prng_engine::max();
						//i++;
					}
					while (r1 == 0);
					r2 = (float) prng() / (float) sitmo::prng_engine::max();
					r1 = boxMullerRadialAmplitude(r1);
					//i++;
				
					if constexpr (ignorekernel == 0)
					{
						potFT(k) = Cplx(cos(2.0f * M_PI * r2), -sin(2.0f * M_PI * r2)) * (1.0f + 7.5f * static_cast<float>(coeff) / k2) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s / potFT(k);
					}
					else
					{
						potFT(k) = Cplx(cos(2.0f * M_PI * r2), -sin(2.0f * M_PI * r2)) * r1 * gsl_spline_eval(pkspline, sqrt(k2), acc) * s;
					}
				}
				
				//prng.discard(huge_skip - (uint64_t) i);
			}
			//prng.discard(huge_skip * (huge_skip - (uint64_t) (kymax - kyend + 1)));
		}

		gsl_interp_accel_free(acc);
		}
	}

	if (linesize > STACK_ALLOCATION_LIMIT)
		free(sinc);

	nvtxRangePop();
}
#endif


//////////////////////////
// initializeParticlePositions
//////////////////////////
// Description:
//   initializes particle positions using a homogeneous template
// 
// Arguments:
//   numpart    number of particles of the template
//   partdata   particle positions in the template
//   numtile    tiling factor for homogeneous template - total particle number will be
//              numpart * numtile^3
//   pcls       reference to (empty) particle object which will contain the new particle ensemble
//
// Returns:
// 
//////////////////////////

template <typename part_handler>
void initializeParticlePositions(const long numpart, const float * partdata, const int numtile, part_handler & pcls)
{	
	part_simple part;
	
	part.vel[0] = Real(0);
	part.vel[1] = Real(0);
	part.vel[2] = Real(0);
	
	for (int ztile = (pcls.lattice().coordSkip()[0] * numtile) / pcls.lattice().size(2); ztile <= ((pcls.lattice().coordSkip()[0] + pcls.lattice().sizeLocal(2)) * numtile) / pcls.lattice().size(2); ztile++)
	{
		if (ztile >= numtile) break;
		
		for (int ytile = (pcls.lattice().coordSkip()[1] * numtile) / pcls.lattice().size(1); ytile <= ((pcls.lattice().coordSkip()[1] + pcls.lattice().sizeLocal(1)) * numtile) / pcls.lattice().size(1); ytile++)
		{
			if (ytile >= numtile) break;
			
			for (int xtile = 0; xtile < numtile; xtile++)
			{
				for (int i = 0; i < numpart; i++)
				{
					part.pos[0] = ((Real) xtile + partdata[3*i]) / (Real) numtile;
					part.pos[1] = ((Real) ytile + partdata[3*i+1]) / (Real) numtile;
					part.pos[2] = ((Real) ztile + partdata[3*i+2]) / (Real) numtile;
					
					part.ID = (long) i + (long) numpart * ((long) xtile + (long) numtile * ((long) ytile + (long) numtile * (long) ztile));
					
					pcls.addParticle_global(part);
				}
			}
		}
	}
}


//////////////////////////
// applyMomentumDistribution
//////////////////////////
// Description:
//   adds a random momentum vector drawn from a Fermi-Dirac distribution to
//   each particle. The current implementation uses a simple rejection-sampling
//   method based on the ziggurat algorithm [G. Marsaglia and W.W. Tsang,
//   J. Stat. Softw. 5 (2000) 1]. The "ziggurat" is hard-coded and was precomputed
//   for the ultrarelativistic limit of a Fermi-Dirac distribution with zero
//   chemical potential. A method to construct the ziggurat on-the-fly for more
//   general distribution functions could be implemented in the future.
// 
// Arguments:
//   pcls   pointer to particle handler
//   seed   seed for random number generator
//   T_m    dimensionless parameter in the distribution function
//          in most cases the ratio of temperature and fundamental mass
//   delta  Field containing a local dT/T (optional)
//
// Returns: sum of momenta over all particles (for reduction)
// 
//////////////////////////

double applyMomentumDistribution(Particles<part_simple,part_simple_info,part_simple_dataType> * pcls, unsigned int seed, float T_m = 0., Field<Real> * delta = NULL)
{	
	Site xPart(pcls->lattice());
	Site x;
	sitmo::prng_engine prng;
	float r1, r2, q, dT;
	double l, dummy, d[3];
	uint32_t i, r;
	double sum_q = 0.0;
	
	const float ql[] = {0.0f,       0.0453329523f, 0.0851601009f, 0.115766097f,
	                    0.142169202f, 0.166069623f, 0.188283033f, 0.209275639f,
	                    0.229344099f, 0.248691555f, 0.267464827f, 0.285774552f,
	                    0.303706948f, 0.321331123f, 0.338703809f, 0.355872580f,
	                    0.372878086f, 0.389755674f, 0.406536572f, 0.423248791f,
	                    0.439917819f, 0.456567164f, 0.473218795f, 0.489893502f,
	                    0.506611198f, 0.523391180f, 0.540252353f, 0.557213439f,
	                    0.574293166f, 0.591510448f, 0.608884565f, 0.626435339f,
	                    0.644183319f, 0.662149971f, 0.680357893f, 0.698831041f,
	                    0.717594982f, 0.736677192f, 0.756107381f, 0.775917886f,
	                    0.796144127f, 0.816825143f, 0.838004248f, 0.859729814f,
	                    0.882056241f, 0.905045149f, 0.928766878f, 0.953302387f,
	                    0.978745698f, 1.00520708f,  1.03281729f,  1.06173322f,
	                    1.09214584f,  1.12429121f,  1.15846661f,  1.19505456f,
	                    1.23456031f,  1.27767280f,  1.32536981f,  1.37911302f,
	                    1.44124650f,  1.51592808f,  1.61180199f,  1.75307820f, 29.0f};
	                    
	const float qr[] = {29.0f,       11.8477879f, 10.3339062f, 9.58550750f,
	                    9.08034038f, 8.69584518f, 8.38367575f, 8.11975908f,
	                    7.89035001f, 7.68685171f, 7.50352431f, 7.33634187f,
	                    7.18236952f, 7.03940031f, 6.90573157f, 6.78002122f,
	                    6.66119176f, 6.54836430f, 6.44081188f, 6.33792577f,
	                    6.23919052f, 6.14416529f, 6.05246957f, 5.96377208f,
	                    5.87778212f, 5.79424263f, 5.71292467f, 5.63362289f,
	                    5.55615183f, 5.48034280f, 5.40604138f, 5.33310512f,
	                    5.26140176f, 5.19080749f, 5.12120558f, 5.05248501f,
	                    4.98453932f, 4.91726547f, 4.85056276f, 4.78433177f,
	                    4.71847331f, 4.65288728f, 4.58747144f, 4.52212011f,
	                    4.45672261f, 4.39116154f, 4.32531058f, 4.25903193f,
	                    4.19217309f, 4.12456265f, 4.05600493f, 3.98627278f,
	                    3.91509779f, 3.84215661f, 3.76705129f, 3.68928014f,
	                    3.60819270f, 3.52291720f, 3.43223655f, 3.33436084f,
	                    3.22646839f, 3.10364876f, 2.95592669f, 2.75624893f, 0.0f};
	                    
	const float f[] = {0.0f,          0.0010042516f, 0.0034718142f, 0.00631345904f,
	                  0.00938886471f, 0.0126471707f, 0.0160614805f, 0.0196150986f,
	                   0.0232967063f, 0.0270982040f, 0.0310135938f, 0.0350383391f,
	                   0.0391689707f, 0.0434028310f, 0.0477379005f, 0.0521726764f,
	                   0.0567060859f, 0.0613374223f, 0.0660662983f, 0.0708926105f,
	                   0.0758165136f, 0.0808384013f, 0.0859588926f, 0.0911788222f,
	                   0.0964992355f, 0.101921386f,  0.107446735f,  0.113076958f,
	                   0.118813945f,  0.124659815f,  0.130616922f,  0.136687871f,
	                   0.142875536f,  0.149183078f,  0.155613967f,  0.162172016f,
	                   0.168861408f,  0.175686736f,  0.182653051f,  0.189765914f,
	                   0.197031455f,  0.204456455f,  0.212048433f,  0.219815749f,
	                   0.227767740f,  0.235914877f,  0.244268957f,  0.252843349f,
	                   0.261653294f,  0.270716296f,  0.280052614f,  0.289685921f,
	                   0.299644172f,  0.309960782f,  0.320676286f,  0.331840691f,
	                   0.343516979f,  0.355786485f,  0.368757588f,  0.382580625f,
	                   0.397475563f,  0.413789108f,  0.432131941f,  0.453799050f, 0.482830296f};
	
	if (delta != NULL)
	{
		l = (double) delta->lattice().size(0);
		x.initialize(delta->lattice());
		x.first();
	}
	                    
	for (xPart.first(); xPart.test(); xPart.next())
	{
		for (auto it=(pcls->field())(xPart).parts.begin(); it != (pcls->field())(xPart).parts.end(); ++it)
		{
			prng.seed(seed);
			prng.discard((uint64_t) (7l * (*it).ID));
				
			while (true)
			{
				r = prng();
				i = r % 64;
				r /= 64;
					
				q = ql[i] + 64.0f * (qr[i]-ql[i]) * ((float) r / (float) sitmo::prng_engine::max());
					
				if (q > ql[i+1] && q < qr[i+1]) break;
					
				if (f[i] + (f[i+1]-f[i]) * ((float) prng() / (float) sitmo::prng_engine::max()) < q * q / (exp(q) + 1.0f)) break;
			}
				
			r1 = 2 * ((float) prng() / (float) sitmo::prng_engine::max()) - 1;
			r2 = 2 * M_PI * ((float) prng() / (float) sitmo::prng_engine::max());
				
			if (delta != NULL)
			{
				for (i = 0; i < 3; i++)
					d[i] = modf((*it).pos[i] * l, &dummy);
					
				dT = (1. - d[0]) * (1. - d[1]) * (1. - d[2]) * (*delta)(x);
				dT += d[0] * (1. - d[1]) * (1. - d[2]) * (*delta)(x+0);
				dT += (1. - d[0]) * d[1] * (1. - d[2]) * (*delta)(x+1);
				dT += d[0] * d[1] * (1. - d[2]) * (*delta)(x+0+1);
				dT += (1. - d[0]) * (1. - d[1]) * d[2] * (*delta)(x+2);
				dT += d[0] * (1. - d[1]) * d[2] * (*delta)(x+0+2);
				dT += (1. - d[0]) * d[1] * d[2] * (*delta)(x+1+2);
				dT += d[0] * d[1] * d[2] * (*delta)(x+2);
					
				q *= T_m * (1. + dT);
			}
			else
				q *= T_m;
				
			(*it).vel[0] += cos(r2) * sqrt(1 - r1*r1) * q;
			(*it).vel[1] += sin(r2) * sqrt(1 - r1*r1) * q;
			(*it).vel[2] += r1 * q;
			
			sum_q += q;
		}
		
		if (delta != NULL) x.next();
	}
	
	return sum_q;
}


#ifdef FFT3D

//////////////////////////
// generateIC_basic
//////////////////////////
// Description:
//   basic initial condition generator
// 
// Arguments:
//   sim            simulation metadata structure
//   ic             settings for IC generation
//   cosmo          cosmological parameter structure
//   fourpiG        4 pi G (in code units)
//   pcls_cdm       pointer to (uninitialized) particle handler for CDM
//   pcls_b         pointer to (uninitialized) particle handler for baryons
//   pcls_ncdm      array of (uninitialized) particle handlers for
//                  non-cold DM (may be set to NULL)
//   maxvel         array that will contain the maximum q/m/a (max. velocity)
//   phi            pointer to allocated field
//   chi            pointer to allocated field
//   Bi             pointer to allocated field
//   source         pointer to allocated field
//   Sij            pointer to allocated field
//   scalarFT       pointer to allocated field
//   BiFT           pointer to allocated field
//   SijFT          pointer to allocated field
//   plan_phi       pointer to FFT planner
//   plan_chi       pointer to FFT planner
//   plan_Bi        pointer to FFT planner
//   plan_source    pointer to FFT planner
//   plan_Sij       pointer to FFT planner
//   class_background  reference to CLASS background structure
//   class_perturbs    reference to CLASS perturbation structure
//   params         pointer to array of precision settings for CLASS (can be NULL)
//   numparam       number of precision settings for CLASS (can be 0)
//
// Returns:
// 
//////////////////////////

void generateIC_basic(metadata & sim, icsettings & ic, cosmology & cosmo, const double fourpiG, perfParticles<part_simple,part_simple_info> * pcls_cdm, perfParticles<part_simple,part_simple_info> * pcls_b, Particles<part_simple,part_simple_info,part_simple_dataType> * pcls_ncdm, double * maxvel, Field<Real> * phi, Field<Real> * chi, Field<Real> * Bi, Field<Real> * source, Field<Real> * Sij, Field<Cplx> * scalarFT, Field<Cplx> * BiFT, Field<Cplx> * SijFT, PlanFFT<Cplx> * plan_phi, PlanFFT<Cplx> * plan_chi, PlanFFT<Cplx> * plan_Bi, PlanFFT<Cplx> * plan_source, PlanFFT<Cplx> * plan_Sij,
#ifdef HAVE_CLASS
background & class_background, perturbs & class_perturbs,
#endif
parameter * params, int & numparam)
{
	double a = 1. / (1. + sim.z_in);
	float * pcldata = NULL;
	gsl_spline * pkspline = NULL;
	gsl_spline * dgaugespline = NULL;
	gsl_spline * vgaugespline = NULL;
	gsl_spline * tk_d1 = NULL;
	gsl_spline * tk_d2 = NULL;
	gsl_spline * tk_t1 = NULL;
	gsl_spline * tk_t2 = NULL;
	double * temp1 = NULL;
	double * temp2 = NULL;
	Site x(phi->lattice());
	rKSite kFT(scalarFT->lattice());
	double max_displacement;
	double rescale;
	double mean_q;
	part_simple_info pcls_cdm_info;
	part_simple_dataType pcls_cdm_dataType;
	part_simple_info pcls_b_info;
	part_simple_dataType pcls_b_dataType;
	part_simple_info pcls_ncdm_info[MAX_PCL_SPECIES];
	part_simple_dataType pcls_ncdm_dataType;
	Real boxSize[3] = {1.,1.,1.};
	char ncdm_name[24];
	Field<Real> * ic_fields[2];
	
	ic_fields[0] = chi;
	ic_fields[1] = phi;

#ifdef HAVE_CLASS
	gsl_interp_accel * acc = gsl_interp_accel_alloc();
#endif
	
	loadHomogeneousTemplate(ic.pclfile[0], sim.numpcl[0], pcldata);
	
	if (pcldata == NULL)
	{
		COUT << " error: particle data was empty!" << endl;
		parallel.abortForce();
	}

	uint64_t capacity = (16L * sim.numpcl[0] * (long) ic.numtile[0] * (long) ic.numtile[0] * (long) ic.numtile[0]) / (parallel.size() * 15L);

	plan_source->preallocate(device_workspace_preallocation_bytes(
		(capacity + PCL_EXTRA_CAPACITY) * 3L * sizeof(Real), PCLBUFFER));
	
	nvtxRangePushA("generateCICKernel");
	if (ic.flags & ICFLAG_CORRECT_DISPLACEMENT)
		generateCICKernel(*source, sim.numpcl[0], pcldata, ic.numtile[0]);
	else
		generateCICKernel(*source);	
	nvtxRangePop();
	
	nvtxRangePushA("generate displacement field(s)");
	plan_source->execute(FFT_FORWARD);
	
	if (ic.pkfile[0] != '\0')	// initial displacements & velocities are derived from a single power spectrum
	{
		loadPowerSpectrum(ic.pkfile, pkspline, sim.boxsize);
	
		if (pkspline == NULL)
		{
			COUT << " error: power spectrum was empty!" << endl;
			parallel.abortForce();
		}
		
		temp1 = (double *) malloc(pkspline->size * sizeof(double));
		temp2 = (double *) malloc(pkspline->size * sizeof(double));

#pragma omp parallel for		
		for (int i = 0; i < pkspline->size; i++)
		{
			temp1[i] = pkspline->x[i];
			temp2[i] = pkspline->y[i] / sim.boxsize / sim.boxsize;
		}

		int size = pkspline->size;
		gsl_spline_free(pkspline);
		pkspline = gsl_spline_alloc(gsl_interp_cspline, size);
		gsl_spline_init(pkspline, temp1, temp2, size);
	
		generateDisplacementField(*scalarFT, sim.gr_flag * Hconf(a, fourpiG, cosmo) * Hconf(a, fourpiG, cosmo), pkspline, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE);
	}
	else					// initial displacements and velocities are set by individual transfer functions
	{
#ifdef HAVE_CLASS
		if (ic.tkfile[0] == '\0')
		{
			nvtxRangePushA("CLASS");
			initializeCLASSstructures(sim, ic, cosmo, class_background, class_perturbs, params, numparam);
			nvtxRangePop();

			loadBGFunctions(class_background, cosmo.Hspline, "H [1/Mpc]", 1.05 * sim.z_in + 0.05, sim.boxsize/cosmo.h);
			cosmo.acc_H = gsl_interp_accel_alloc();
			
			loadTransferFunctions(class_background, class_perturbs, tk_d1, tk_t1, NULL, sim.boxsize, sim.z_in, cosmo.h);
		}
		else
#endif
		loadTransferFunctions(ic.tkfile, tk_d1, tk_t1, NULL, sim.boxsize, cosmo.h);
		
		if (tk_d1 == NULL || tk_t1 == NULL)
		{
			COUT << " error: total transfer function was empty!" << endl;
			parallel.abortForce();
		}
		
		temp1 = (double *) malloc(tk_d1->size * sizeof(double));
		temp2 = (double *) malloc(tk_d1->size * sizeof(double));

#pragma omp parallel for		
		for (int i = 0; i < tk_d1->size; i++) // construct phi
			temp1[i] = -M_PI * tk_d1->y[i] * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) * tk_d1->x[i];

		pkspline = gsl_spline_alloc(gsl_interp_cspline, tk_d1->size);
		gsl_spline_init(pkspline, tk_d1->x, temp1, tk_d1->size);
				
#ifdef HAVE_CLASS
		if (ic.tkfile[0] == '\0')
		{
			loadTransferFunctions(class_background, class_perturbs, tk_d2, tk_t2, "eta", sim.boxsize, sim.z_in, cosmo.h);
			rescale = -gsl_spline_eval_deriv(cosmo.Hspline, a, cosmo.acc_H) * a * a;
		}
		else
#endif
		{
			loadTransferFunctions(ic.tkfile, tk_d2, tk_t2, "eta", sim.boxsize, cosmo.h);
				
			if (tk_d2 == NULL || tk_t2 == NULL)
			{
				COUT << " error: transfer functions were empty!" << endl;
				parallel.abortForce();
			}
				
			rescale = (Hconf(a, fourpiG, cosmo) - 100. * (Hconf(1.005 * a, fourpiG, cosmo) - Hconf(0.995 * a, fourpiG, cosmo)));
		}

		if (sim.gr_flag == 0)
		{	
#pragma omp parallel for
			for (int i = 0; i < tk_t2->size; i++) // construct gauge correction for N-body gauge (3 Hconf theta_tot / k^2 = 3 Hconf eta' / (Hconf^2 - Hconf'))
				temp2[i] = 3. * tk_t2->y[i] / rescale;
		}
		else
		{
#pragma omp parallel for
			for (int i = 0; i < tk_t2->size; i++) // construct gauge correction for Poisson gauge (3 phi - 3 eta)
				temp2[i] = 3. * tk_d1->y[i] - 3. * tk_d2->y[i];
		}

		dgaugespline = gsl_spline_alloc(gsl_interp_cspline, tk_t2->size);
		gsl_spline_init(dgaugespline, tk_t2->x, temp2, tk_t2->size);

		gsl_spline_free(tk_d1);
		gsl_spline_free(tk_t1);

		if (sim.gr_flag > 0)
		{
#pragma omp parallel for
			for (int i = 0; i < tk_t2->size; i++)
				temp1[i] = 3. * tk_t2->y[i];
		}
		else
		{
			if (ic.tkfile[0] != '\0')
			{
				COUT << COLORTEXT_YELLOW << " /!\\ warning" << COLORTEXT_RESET << ": gauge correction for N-body gauge velocities cannot be computed accurately from the transfer functions at a single redshift!" << endl;
			}

#pragma omp parallel for
			for (int i = 0; i < tk_t2->size; i++)
				temp1[i] = 0.;
		}

#ifdef HAVE_CLASS
		if (ic.tkfile[0] == '\0')
		{	
			if (sim.gr_flag == 0)
			{
				loadTransferFunctions(class_background, class_perturbs, tk_d1, tk_t1, "eta", sim.boxsize, 1./1.01/a - 1., cosmo.h);
				rescale = -gsl_spline_eval_deriv(cosmo.Hspline, 1.01*a, cosmo.acc_H) * a * a * 1.0201;

#pragma omp parallel for
				for (int i = 0; i < tk_t1->size; i++)
					temp1[i] = (temp2[i] - 3. * tk_t1->y[i] / rescale) * 100. * Hconf(1.005 * a, fourpiG, cosmo);

				gsl_spline_free(tk_d1);
				gsl_spline_free(tk_t1);
			}
			
			loadTransferFunctions(class_background, class_perturbs, tk_d1, tk_t1, "h", sim.boxsize, sim.z_in, cosmo.h);
		}
		else
#endif
		{
			loadTransferFunctions(ic.tkfile, tk_d1, tk_t1, "h", sim.boxsize, cosmo.h);
			
			if (tk_d1 == NULL || tk_t1 == NULL)
			{
				COUT << " error: transfer functions were empty!" << endl;
				parallel.abortForce();
			}
		}

#pragma omp parallel for
		for (int i = 0; i < tk_t1->size; i++)
			temp1[i] += tk_t1->y[i] * 0.5;

		gsl_spline_free(tk_d1);
		gsl_spline_free(tk_t1);

		vgaugespline = gsl_spline_alloc(gsl_interp_cspline, tk_t2->size);
		gsl_spline_init(vgaugespline, tk_t2->x, temp1, tk_t2->size);

		gsl_spline_free(tk_d2);
		gsl_spline_free(tk_t2);

#ifdef HAVE_CLASS
		if (ic.tkfile[0] == '\0')
			loadTransferFunctions(class_background, class_perturbs, tk_d1, tk_t1, "d_m", sim.boxsize, sim.z_in, cosmo.h);
		else
#endif		
		loadTransferFunctions(ic.tkfile, tk_d1, tk_t1, "d_m", sim.boxsize, cosmo.h);	// get transfer functions for CDM
		
		if (tk_d1 == NULL || tk_t1 == NULL)
		{
			COUT << " error: matter transfer function was empty!" << endl;
			parallel.abortForce();
		}
		
		if (sim.baryon_flag > 0)
		{
#ifdef HAVE_CLASS
			if (ic.tkfile[0] == '\0')
				loadTransferFunctions(class_background, class_perturbs, tk_d2, tk_t2, "b", sim.boxsize, sim.z_in, cosmo.h);
			else
#endif
			loadTransferFunctions(ic.tkfile, tk_d2, tk_t2, "b", sim.boxsize, cosmo.h);	// get transfer functions for baryons
		
			if (tk_d2 == NULL || tk_t2 == NULL)
			{
				COUT << " error: baryon transfer function was empty!" << endl;
				parallel.abortForce();
			}
			if (tk_d2->size != tk_d1->size)
			{
				COUT << " error: baryon transfer function line number mismatch!" << endl;
				parallel.abortForce();
			}
		}
		
		if (sim.baryon_flag == 2)	// baryon treatment = blend; compute displacement & velocity from weighted average
		{
#pragma omp parallel for
			for (int i = 0; i < tk_d1->size; i++)
			{
					temp1[i] = (sim.gr_flag > 0 ? -3. * pkspline->y[i] / pkspline->x[i] / pkspline->x[i] : 0.) - ((cosmo.Omega_cdm * tk_d1->y[i] + cosmo.Omega_b * tk_d2->y[i]) / (cosmo.Omega_cdm + cosmo.Omega_b) + dgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
					temp2[i] = -a * ((cosmo.Omega_b * tk_t2->y[i]) / (cosmo.Omega_cdm + cosmo.Omega_b) + vgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
			}

			gsl_spline_free(tk_d1);
			gsl_spline_free(tk_t1);
			tk_d1 = gsl_spline_alloc(gsl_interp_cspline, tk_d2->size);
			tk_t1 = gsl_spline_alloc(gsl_interp_cspline, tk_d2->size);
			gsl_spline_init(tk_d1, tk_d2->x, temp1, tk_d2->size);
			gsl_spline_init(tk_t1, tk_d2->x, temp2, tk_d2->size);
			gsl_spline_free(tk_d2);
			gsl_spline_free(tk_t2);
		}
		
		if (sim.baryon_flag == 3)	// baryon treatment = hybrid; compute displacement & velocity from weighted average (sub-species)
		{
			if (8. * cosmo.Omega_b / (cosmo.Omega_cdm + cosmo.Omega_b) > 1.)
			{
#pragma omp parallel for
				for (int i = 0; i < tk_d1->size; i++)
				{
						temp1[i] = (sim.gr_flag > 0 ? -3. * pkspline->y[i] / pkspline->x[i] / pkspline->x[i] : 0.) - ((8. * cosmo.Omega_cdm * tk_d1->y[i] + (7. * cosmo.Omega_b - cosmo.Omega_cdm) * tk_d2->y[i]) / (cosmo.Omega_cdm + cosmo.Omega_b) / 7. + dgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
						temp2[i] = -a * (((7. * cosmo.Omega_b - cosmo.Omega_cdm) * tk_t2->y[i]) / (cosmo.Omega_cdm + cosmo.Omega_b) / 7. + vgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
				}
				
				gsl_spline_free(tk_d1);
				gsl_spline_free(tk_t1);
				tk_d1 = gsl_spline_alloc(gsl_interp_cspline, tk_d2->size);
				tk_t1 = gsl_spline_alloc(gsl_interp_cspline, tk_d2->size);
				gsl_spline_init(tk_d1, tk_d2->x, temp1, tk_d2->size);
				gsl_spline_init(tk_t1, tk_d2->x, temp2, tk_d2->size);
			}
			else
			{
#pragma omp parallel for
				for (int i = 0; i < tk_d1->size; i++)
				{
					temp1[i] = (sim.gr_flag > 0 ? -3. * pkspline->y[i] / pkspline->x[i] / pkspline->x[i] : 0.) - (((cosmo.Omega_cdm - 7. * cosmo.Omega_b) * tk_d1->y[i] + 8. * cosmo.Omega_b * tk_d2->y[i]) / (cosmo.Omega_cdm + cosmo.Omega_b) + dgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
					temp2[i] = -a * ((8. * cosmo.Omega_b * tk_t2->y[i]) / (cosmo.Omega_cdm + cosmo.Omega_b) + vgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
				}

				gsl_spline_free(tk_d2);
				gsl_spline_free(tk_t2);
				tk_d2 = gsl_spline_alloc(gsl_interp_cspline, tk_d1->size);
				tk_t2 = gsl_spline_alloc(gsl_interp_cspline, tk_d1->size);
				gsl_spline_init(tk_d2, tk_d1->x, temp1, tk_d1->size);
				gsl_spline_init(tk_t2, tk_d1->x, temp2, tk_d1->size);
			}
		}
		
		if (sim.baryon_flag == 1 || (sim.baryon_flag == 3 && 8. * cosmo.Omega_b / (cosmo.Omega_cdm + cosmo.Omega_b) > 1.)) // compute baryonic displacement & velocity
		{
#pragma omp parallel for
			for (int i = 0; i < tk_d1->size; i++)
			{
					temp1[i] = (sim.gr_flag > 0 ? -3. * pkspline->y[i] / pkspline->x[i] / pkspline->x[i] : 0.) - (tk_d2->y[i] + dgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d2->x[i] * cosmo.h / sim.boxsize, ic) / tk_d2->x[i]) / tk_d2->x[i];
					temp2[i] = -a * (tk_t2->y[i] + vgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d2->x[i] * cosmo.h / sim.boxsize, ic) / tk_d2->x[i]) / tk_d2->x[i];			
			}

			gsl_spline_free(tk_d2);
			gsl_spline_free(tk_t2);
			tk_d2 = gsl_spline_alloc(gsl_interp_cspline, tk_d1->size);
			tk_t2 = gsl_spline_alloc(gsl_interp_cspline, tk_d1->size);
			gsl_spline_init(tk_d2, tk_d1->x, temp1, tk_d1->size);
			gsl_spline_init(tk_t2, tk_d1->x, temp2, tk_d1->size);
		}
		
		if (sim.baryon_flag < 2 || (sim.baryon_flag == 3 && 8. * cosmo.Omega_b / (cosmo.Omega_cdm + cosmo.Omega_b) <= 1.))	// compute CDM displacement & velocity
		{
#pragma omp parallel for
			for (int i = 0; i < tk_d1->size; i++)
			{
					temp1[i] = (sim.gr_flag > 0 ? -3. * pkspline->y[i] / pkspline->x[i] / pkspline->x[i] : 0.) - (tk_d1->y[i] + dgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
					temp2[i] = -a * vgaugespline->y[i] * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
			}

			gsl_spline_free(tk_d1);
			gsl_spline_free(tk_t1);
			tk_d1 = gsl_spline_alloc(gsl_interp_cspline, pkspline->size);
			tk_t1 = gsl_spline_alloc(gsl_interp_cspline, pkspline->size);
			gsl_spline_init(tk_d1, pkspline->x, temp1, pkspline->size);
			gsl_spline_init(tk_t1, pkspline->x, temp2, pkspline->size);
		}
		
		if ((sim.baryon_flag == 1 && !(ic.flags & ICFLAG_CORRECT_DISPLACEMENT)) || sim.baryon_flag == 3)
		{
			generateDisplacementField(*scalarFT, 0., tk_d2, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE);
			gsl_spline_free(tk_d2);
			plan_phi->execute(FFT_BACKWARD);
			phi->updateHalo();	// phi now contains the baryonic displacement
			plan_source->execute(FFT_FORWARD);
		}
		
		generateDisplacementField(*scalarFT, 0., tk_d1, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE);
		gsl_spline_free(tk_d1);
	}

	plan_chi->execute(FFT_BACKWARD);
	chi->updateHalo();	// chi now contains the CDM displacement
	nvtxRangePop();
	
	strcpy(pcls_cdm_info.type_name, "part_simple");
	if (sim.baryon_flag == 1)
		pcls_cdm_info.mass = cosmo.Omega_cdm / (Real) (sim.numpcl[0]*(long)ic.numtile[0]*(long)ic.numtile[0]*(long)ic.numtile[0]);
	else
		pcls_cdm_info.mass = (cosmo.Omega_cdm + cosmo.Omega_b) / (Real) (sim.numpcl[0]*(long)ic.numtile[0]*(long)ic.numtile[0]*(long)ic.numtile[0]);
	pcls_cdm_info.relativistic = false;
	
	//pcls_cdm->initialize(pcls_cdm_info, pcls_cdm_dataType, &(phi->lattice()), boxSize);
	nvtxRangePushA("initialize CDM particles: memory allocation");
	pcls_cdm->initialize(pcls_cdm_info, &(phi->lattice()), boxSize, capacity+PCL_EXTRA_CAPACITY, PCL_EXTRA_CAPACITY);
	nvtxRangePop();
	
	nvtxRangePushA("initialize CDM particles: positions");
	initializeParticlePositions(sim.numpcl[0], pcldata, ic.numtile[0], *pcls_cdm);
	pcls_cdm->updateRowBuffers();
	nvtxRangePop();

	nvtxRangePushA("initialize CDM particles: displacements");
	int op = MAX;
	if (sim.baryon_flag == 3)	// baryon treatment = hybrid; displace particles using both displacement fields
		pcls_cdm->moveParticles(displace_pcls_ic_basic_functor(), 1., ic_fields, 2, NULL, &max_displacement, &op, 1);
	else
		pcls_cdm->moveParticles(displace_pcls_ic_basic_functor(), 1., &chi, 1, NULL, &max_displacement, &op, 1);	// displace CDM particles
	nvtxRangePop();
	
	sim.numpcl[0] *= (long) ic.numtile[0] * (long) ic.numtile[0] * (long) ic.numtile[0];
	
	COUT << " " << sim.numpcl[0] << " cdm particles initialized: maximum displacement = " << max_displacement * sim.numpts << " lattice units." << endl;
	
	free(pcldata);
	
	if (sim.baryon_flag == 1)
	{
		loadHomogeneousTemplate(ic.pclfile[1], sim.numpcl[1], pcldata);
	
		if (pcldata == NULL)
		{
			COUT << " error: particle data was empty!" << endl;
			parallel.abortForce();
		}
		
		if (ic.flags & ICFLAG_CORRECT_DISPLACEMENT)
		{
			generateCICKernel(*phi, sim.numpcl[1], pcldata, ic.numtile[1]);
			plan_phi->execute(FFT_FORWARD);
			generateDisplacementField(*scalarFT, 0., tk_d2, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE);
			gsl_spline_free(tk_d2);
			plan_phi->execute(FFT_BACKWARD);
			phi->updateHalo();
		}
		
		strcpy(pcls_b_info.type_name, "part_simple");
		pcls_b_info.mass = cosmo.Omega_b / (Real) (sim.numpcl[1]*(long)ic.numtile[1]*(long)ic.numtile[1]*(long)ic.numtile[1]);
		pcls_b_info.relativistic = false;

		capacity = (16L * sim.numpcl[1] * (long) ic.numtile[1] * (long) ic.numtile[1] * (long) ic.numtile[1]) / (parallel.size() * 15L);
	
		//pcls_b->initialize(pcls_b_info, pcls_b_dataType, &(phi->lattice()), boxSize);
		pcls_b->initialize(pcls_b_info, &(phi->lattice()), boxSize, capacity+PCL_EXTRA_CAPACITY, PCL_EXTRA_CAPACITY);
	
		initializeParticlePositions(sim.numpcl[1], pcldata, ic.numtile[1], *pcls_b);
		pcls_b->updateRowBuffers();

		pcls_b->moveParticles(displace_pcls_ic_basic_functor(), 1., &phi, 1, NULL, &max_displacement, &op, 1);	// displace baryon particles
	
		sim.numpcl[1] *= (long) ic.numtile[1] * (long) ic.numtile[1] * (long) ic.numtile[1];
	
		COUT << " " << sim.numpcl[1] << " baryon particles initialized: maximum displacement = " << max_displacement * sim.numpts << " lattice units." << endl;
	
		free(pcldata);
	}
	
	nvtxRangePushA("initialize CDM particles: velocities");
	if (ic.pkfile[0] == '\0')	// set velocities using transfer functions
	{
		if (ic.flags & ICFLAG_CORRECT_DISPLACEMENT)
			generateCICKernel(*source);
		
		plan_source->execute(FFT_FORWARD);
		
		if (sim.baryon_flag == 1 || sim.baryon_flag == 3)
		{
			generateDisplacementField(*scalarFT, 0., tk_t2, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE, 0);
			plan_phi->execute(FFT_BACKWARD);
			phi->updateHalo();	// phi now contains the baryonic velocity potential
			gsl_spline_free(tk_t2);
			plan_source->execute(FFT_FORWARD);
		}
		
		generateDisplacementField(*scalarFT, 0., tk_t1, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE, 0);
		plan_chi->execute(FFT_BACKWARD);
		chi->updateHalo();	// chi now contains the CDM velocity potential
		gsl_spline_free(tk_t1);			
		
		if (sim.baryon_flag == 3)	// baryon treatment = hybrid; set velocities using both velocity potentials
			maxvel[0] = pcls_cdm->updateVel(initialize_q_ic_basic_functor(), 1., ic_fields, 2) / a;
		else
			maxvel[0] = pcls_cdm->updateVel(initialize_q_ic_basic_functor(), 1., &chi, 1) / a;	// set CDM velocities
		
		if (sim.baryon_flag == 1)
			maxvel[1] = pcls_b->updateVel(initialize_q_ic_basic_functor(), 1., &phi, 1) / a;	// set baryon velocities
	}
	nvtxRangePop();
	
	if (sim.baryon_flag > 1) sim.baryon_flag = 0;
	
	for (int p = 0; p < cosmo.num_ncdm; p++)	// initialization of non-CDM species
	{
		if (ic.numtile[1+sim.baryon_flag+p] < 1) continue;

		loadHomogeneousTemplate(ic.pclfile[1+sim.baryon_flag+p], sim.numpcl[1+sim.baryon_flag+p], pcldata);
	
		if (pcldata == NULL)
		{
			COUT << " error: particle data was empty!" << endl;
			parallel.abortForce();
		}
		
		if (ic.pkfile[0] == '\0')
		{
			sprintf(ncdm_name, "ncdm[%d]", p);
#ifdef HAVE_CLASS
			if (ic.tkfile[0] == '\0')
				loadTransferFunctions(class_background, class_perturbs, tk_d1, tk_t1, ncdm_name, sim.boxsize, sim.z_in, cosmo.h);
			else
#endif
			loadTransferFunctions(ic.tkfile, tk_d1, tk_t1, ncdm_name, sim.boxsize, cosmo.h);
		
			if (tk_d1 == NULL || tk_t1 == NULL)
			{
				COUT << " error: ncdm transfer function was empty! (species " << p << ")" << endl;
				parallel.abortForce();
			}
			
#ifdef HAVE_CLASS
			sprintf(ncdm_name, "(.)rho_ncdm[%d]", p);
			loadBGFunctions(class_background, tk_d2, ncdm_name, 1.05 * sim.z_in + 0.05);
			sprintf(ncdm_name, "(.)p_ncdm[%d]", p);
			loadBGFunctions(class_background, tk_t2, ncdm_name, 1.05 * sim.z_in + 0.05);
			rescale = 1./gsl_spline_eval(tk_d2, a, acc);
			gsl_interp_accel_reset(acc);
			rescale *= gsl_spline_eval(tk_t2, a, acc);
			gsl_interp_accel_reset(acc);
			rescale += 1.;
			gsl_spline_free(tk_d2);
			gsl_spline_free(tk_t2);
#else
			rescale = (bg_ncdm(a * 1.005, cosmo, p) - bg_ncdm(a * 0.995, cosmo, p)) * 100. / bg_ncdm(a, cosmo, p);
			rescale = 1. - rescale / 3.;
#endif

#pragma omp parallel for		
			for (int i = 0; i < tk_d1->size; i++)
			{
				temp1[i] = (sim.gr_flag > 0 ? -3. * pkspline->y[i] / pkspline->x[i] / pkspline->x[i] : 0.) - (tk_d1->y[i] + dgaugespline->y[i] * rescale) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
				temp2[i] = -a * (tk_t1->y[i] + vgaugespline->y[i]) * M_PI * sqrt(Pk_primordial(tk_d1->x[i] * cosmo.h / sim.boxsize, ic) / tk_d1->x[i]) / tk_d1->x[i];
			}

			gsl_spline_free(tk_d1);
			gsl_spline_free(tk_t1);
			tk_d1 = gsl_spline_alloc(gsl_interp_cspline, pkspline->size);
			tk_t1 = gsl_spline_alloc(gsl_interp_cspline, pkspline->size);
			gsl_spline_init(tk_d1, pkspline->x, temp1, pkspline->size);
			gsl_spline_init(tk_t1, pkspline->x, temp2, pkspline->size);
			
			plan_source->execute(FFT_FORWARD);
			generateDisplacementField(*scalarFT, 0., tk_d1, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE);
			plan_chi->execute(FFT_BACKWARD);	// chi now contains the displacement for the non-CDM species
			chi->updateHalo();
			gsl_spline_free(tk_d1);
			
			plan_source->execute(FFT_FORWARD);
			generateDisplacementField(*scalarFT, 0., tk_t1, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE, 0);
			plan_phi->execute(FFT_BACKWARD);	// phi now contains the velocity potential for the non-CDM species
			phi->updateHalo();
			gsl_spline_free(tk_t1);
		}
		
		strcpy(pcls_ncdm_info[p].type_name, "part_simple");
		pcls_ncdm_info[p].mass = cosmo.Omega_ncdm[p] / (Real) (sim.numpcl[1+sim.baryon_flag+p]*(long)ic.numtile[1+sim.baryon_flag+p]*(long)ic.numtile[1+sim.baryon_flag+p]*(long)ic.numtile[1+sim.baryon_flag+p]);
		pcls_ncdm_info[p].relativistic = true;
		
		pcls_ncdm[p].initialize(pcls_ncdm_info[p], pcls_ncdm_dataType, &(phi->lattice()), boxSize);
		
		initializeParticlePositions(sim.numpcl[1+sim.baryon_flag+p], pcldata, ic.numtile[1+sim.baryon_flag+p], pcls_ncdm[p]);
		
		pcls_ncdm[p].moveParticles(displace_pcls_ic_basic, 1., &chi, 1, NULL, &max_displacement, &op, 1);	// displace non-CDM particles
		
		sim.numpcl[1+sim.baryon_flag+p] *= (long) ic.numtile[1+sim.baryon_flag+p] * (long) ic.numtile[1+sim.baryon_flag+p] * (long) ic.numtile[1+sim.baryon_flag+p];
	
		COUT << " " << sim.numpcl[1+sim.baryon_flag+p] << " ncdm particles initialized for species " << p+1 << ": maximum displacement = " << max_displacement * sim.numpts << " lattice units." << endl;
		
		free(pcldata);
		
		if (ic.pkfile[0] == '\0')	// set non-CDM velocities using transfer functions
			pcls_ncdm[p].updateVel(initialize_q_ic_basic, 1., &phi, 1);
	}

	free(temp1);
	free(temp2);
	
	nvtxRangePushA("initialize phi");
	if (ic.pkfile[0] == '\0')
	{
		plan_source->execute(FFT_FORWARD);
		generateDisplacementField(*scalarFT, 0., pkspline, (unsigned int) ic.seed, ic.flags & ICFLAG_KSPHERE, 0);
	}
	else
	{
		//projection_init(source);
		thrust::fill_n(thrust::device, source->data(), source->lattice().sitesLocalGross(), Real(0));
		scalarProjectionCIC_project(pcls_cdm, source);
		if (sim.baryon_flag)
			scalarProjectionCIC_project(pcls_b, source);
		for (int p = 0; p < cosmo.num_ncdm; p++)
		{
			if (ic.numtile[1+sim.baryon_flag+p] < 1) continue;
			scalarProjectionCIC_project(pcls_ncdm+p, source);
		}
		scalarProjectionCIC_comm(source);
	
		plan_source->execute(FFT_FORWARD);
	
		kFT.first();
		if (kFT.coord(0) == 0 && kFT.coord(1) == 0 && kFT.coord(2) == 0)
			(*scalarFT)(kFT) = Cplx(0.,0.);
				
		solveModifiedPoissonFT(*scalarFT, *scalarFT, fourpiG / a, 3. * sim.gr_flag * (Hconf(a, fourpiG, cosmo) * Hconf(a, fourpiG, cosmo) + fourpiG * cosmo.Omega_m / a));
	}
	
	plan_phi->execute(FFT_BACKWARD);
	phi->updateHalo();	// phi now finally contains phi
	nvtxRangePop();
	
	if (ic.pkfile[0] != '\0')	// if power spectrum is used instead of transfer functions, set velocities using linear approximation
	{	
		rescale = a / Hconf(a, fourpiG, cosmo) / (1.5 * Omega_m(a, cosmo) + 2. * Omega_rad(a, cosmo));
		maxvel[0] = pcls_cdm->updateVel(initialize_q_ic_basic_functor(), rescale, &phi, 1) / a;
		if (sim.baryon_flag)
			maxvel[1] = pcls_b->updateVel(initialize_q_ic_basic_functor(), rescale, &phi, 1) / a;
	}
			
	for (int p = 0; p < cosmo.num_ncdm; p++)
	{
		if (ic.numtile[1+sim.baryon_flag+p] < 1)
		{
			maxvel[1+sim.baryon_flag+p] = 0;
			continue;
		}

		if (ic.pkfile[0] != '\0') // if power spectrum is used instead of transfer functions, set bulk velocities using linear approximation
		{		
			rescale = a / Hconf(a, fourpiG, cosmo) / (1.5 * Omega_m(a, cosmo) + Omega_rad(a, cosmo));
			pcls_ncdm[p].updateVel(initialize_q_ic_basic, rescale, &phi, 1);
		}
		
		if (cosmo.m_ncdm[p] > 0.) // add velocity dispersion for non-CDM species
		{
			rescale = pow(cosmo.Omega_g * cosmo.h * cosmo.h / C_PLANCK_LAW, 0.25) * cosmo.T_ncdm[p] * C_BOLTZMANN_CST / cosmo.m_ncdm[p];
			mean_q = applyMomentumDistribution(pcls_ncdm+p, (unsigned int) (ic.seed + p), rescale);
			parallel.sum(mean_q);
			COUT << " species " << p+1 << " Fermi-Dirac distribution had mean q/m = " << mean_q / sim.numpcl[1+sim.baryon_flag+p] << endl;
		}
#ifdef ANISOTROPIC_EXPANSION
		double f_params[7] = {a, 0., 0., 0., 0., 0., 0.};
		maxvel[1+sim.baryon_flag+p] = pcls_ncdm[p].updateVel(update_q, 0., &phi, 1, f_params);
#else
		maxvel[1+sim.baryon_flag+p] = pcls_ncdm[p].updateVel(update_q, 0., &phi, 1, &a);
#endif
	}
	
	nvtxRangePushA("initialize B");
	//projection_init(Bi);
	thrust::fill_n(thrust::device, Bi->data(), 3*Bi->lattice().sitesLocalGross(), Real(0));
	projection_T0i_project(pcls_cdm, Bi, phi);
	if (sim.baryon_flag)
		projection_T0i_project(pcls_b, Bi, phi);
	projection_T0i_comm(Bi);
	plan_Bi->execute(FFT_FORWARD);
	projectFTvector(*BiFT, *BiFT, fourpiG / (double) sim.numpts / (double) sim.numpts);	
	plan_Bi->execute(FFT_BACKWARD);	
	Bi->updateHalo();	// B initialized
	nvtxRangePop();
	
	nvtxRangePushA("initialize chi");
	//projection_init(Sij);
	thrust::fill_n(thrust::device, Sij->data(), 6*Sij->lattice().sitesLocalGross(), Real(0));
	projection_Tij_project(pcls_cdm, Sij, a, phi);
	if (sim.baryon_flag)
		projection_Tij_project(pcls_b, Sij, a, phi);
	projection_Tij_comm(Sij);
	
	prepareFTsource(*phi, *Sij, *Sij, 2. * fourpiG / a / (double) sim.numpts / (double) sim.numpts);	
	plan_Sij->execute(FFT_FORWARD);	
	projectFTscalar(*SijFT, *scalarFT);
	plan_chi->execute(FFT_BACKWARD);		
	chi->updateHalo();	// chi now finally contains chi
	nvtxRangePop();

	gsl_spline_free(pkspline);
	if (dgaugespline != NULL)
		gsl_spline_free(dgaugespline);
	if (vgaugespline != NULL)
		gsl_spline_free(vgaugespline);

#ifdef HAVE_CLASS
	gsl_interp_accel_free(acc);
#endif
}

#endif

#endif