diff --git a/GridDxFormat/grid_reader.c b/GridDxFormat/grid_reader.c
index bc210750ae18d7758493a8d98feb886b79a527ae..990ab146dfe9fed78eed81792ace7b93ce3de40b 100644
--- a/GridDxFormat/grid_reader.c
+++ b/GridDxFormat/grid_reader.c
@@ -9,6 +9,7 @@
#include <assert.h>
#include <unistd.h>
#include <ctype.h>
+#include <math.h>
#include "grid_reader.h"
typedef struct
{
@@ -386,6 +387,130 @@ void set_grid_data(void* src, compute_type compute, void* kernel)
}
}
+static void* create_new_grid(size_t nx, size_t ny, size_t nz, double* origin, double** basis)
+{
+ void* g = malloc(sizeof(grid));
+ ((grid*)g)->nx = nx;
+ ((grid*)g)->ny = ny;
+ ((grid*)g)->nz = nz;
+ ((grid*)g)->size = nx*ny*nz;
+ ((grid*)g)->data = (double*)malloc(sizeof(double)*(((grid*)g)->size));
+
+ ((grid*)g)->origin = (double*)malloc(sizeof(double)*3);
+ for(size_t i = 0; i < 3; ++i)
+ ((grid*)g)->origin[i] = origin[i];
+
+ ((grid*)g)->basis = (double**)malloc(sizeof(double*)*3);
+ for(size_t i = 0; i < 3; ++i)
+ ((grid*)g)->basis[i] = (double*)malloc(sizeof(double)*3);
+
+ for(size_t i = 0; i < 3; ++i)
+ for(size_t j = 0; j < 3; ++j)
+ ((grid*)g)->basis[i][j] = basis[i][j];
+ return g;
+
+}
+
+static size_t Index(size_t nx, size_t ny, size_t nz, int x, int y, int z)
+{
+ x = (x < 0) ? 0 : x;
+ y = (y < 0) ? 0 : y;
+ z = (z < 0) ? 0 : z;
+ x = ((size_t)x >= nx) ? (int)nx-1 : x;
+ y = ((size_t)y >= ny) ? (int)ny-1 : y;
+ z = ((size_t)z >= nz) ? (int)nz-1 : z;
+ return z + nz * (y + ny * x);
+}
+/*
+static void inverse_3by3(double** m, double inv[3][3])
+{
+ double det = m[0][0]*(m[1][1]*m[2][2]-m[1][2]*m[2][1])-
+ m[0][1]*(m[1][0]*m[2][2]-m[1][2]*m[2][0])+
+ m[0][2]*(m[1][0]*m[2][1]-m[1][1]*m[2][0]);
+
+ inv[0][0] = (m[1][1]*m[2][2]-m[1][2]*m[2][1])/det;
+ inv[0][1] =-(m[0][1]*m[2][2]-m[0][2]*m[2][1])/det;
+ inv[0][2] = (m[0][1]*m[1][2]-m[0][2]*m[1][1])/det;
+
+ inv[1][0] =-(m[1][0]*m[2][2]-m[1][2]*m[2][0])/det;
+ inv[1][1] = (m[0][0]*m[2][2]-m[0][2]*m[2][0])/det;
+ inv[1][2] =-(m[0][0]*m[1][2]-m[0][2]*m[1][0])/det;
+
+ inv[2][0] = (m[1][0]*m[2][1]-m[1][1]*m[2][0])/det;
+ inv[2][1] =-(m[0][0]*m[2][1]-m[0][1]*m[2][0])/det;
+ inv[2][2] = (m[0][0]*m[1][1]-m[0][1]*m[1][0])/det;
+
+}*/
+//resample the grid data such that each dimensition is even number by using linear interpolation
+//I currently implement for diagonal basis set
+void* resample_grid_data(void* src)
+{
+ grid* old_g = (grid*)src;
+ grid* new_g;
+ size_t nx, ny, nz;
+ double** b = (double**)malloc(sizeof(double*)*3);
+ for(size_t i = 0; i < 3; ++i)
+ {
+ b[i] = (double*)malloc(sizeof(double)*3);
+ memset(b[i],0,sizeof(double)*3);
+ }
+
+ nx = (old_g->nx/2)*2;
+ ny = (old_g->ny/2)*2;
+ nz = (old_g->nz/2)*2;
+
+ b[0][0] = old_g->basis[0][0]*(old_g->nx/(double)nx);
+ b[1][1] = old_g->basis[1][1]*(old_g->ny/(double)ny);
+ b[2][2] = old_g->basis[2][2]*(old_g->nz/(double)nz);
+
+ new_g = create_new_grid(nx,ny,nz,old_g->origin, b);
+
+ for(size_t i = 0; i < 3; ++i)
+ free(b[i]);
+ free(b);
+
+ for(size_t i = 0; i < new_g->size; ++i)
+ {
+ size_t ix = i / new_g->nz / new_g->ny;
+ size_t iy = i / new_g->nz % new_g->ny;
+ size_t iz = i % new_g->nz;
+
+ double posx = new_g->basis[0][0]*ix;
+ double posy = new_g->basis[1][1]*iy;
+ double posz = new_g->basis[2][2]*iz;
+
+ int x0 = floor(posx/old_g->basis[0][0]);
+ int y0 = floor(posy/old_g->basis[1][1]);
+ int z0 = floor(posz/old_g->basis[2][2]);
+
+ double posx0 = old_g->basis[0][0] * x0;
+ double posy0 = old_g->basis[1][1] * y0;
+ double posz0 = old_g->basis[2][2] * z0;
+
+ int x1 = x0 + 1;
+ int y1 = y0 + 1;
+ int z1 = z0 + 1;
+
+ double posx1 = old_g->basis[0][0] * x1;
+ double posy1 = old_g->basis[1][1] * y1;
+ double posz1 = old_g->basis[2][2] * z1;
+
+ double xd = (posx-posx0)/(posx1-posx0);
+ double yd = (posy-posy0)/(posy1-posy0);
+ double zd = (posz-posz0)/(posz1-posz0);
+ double c0 = ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y0, z0)] +
+ xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y0, z0)])*(1.-yd) +
+ ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y1, z0)] +
+ xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y1, z0)])*yd;
+
+ double c1 = ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y0, z1)] +
+ xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y0, z1)])*(1.-yd) +
+ ((1.-xd)* old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x0, y1, z1)] +
+ xd * old_g->data[Index(old_g->nx, old_g->ny, old_g->nz, x1, y1, z1)])*yd;
+ new_g->data[i] = c0*(1.-zd) + c1*zd;
+ }
+ return new_g;
+}
double* get_grid_origin(void* src)
{
return ((grid*)src)->origin;
diff --git a/Kernel/lenard_jones_kernel.c b/Kernel/lenard_jones_kernel.c
index cb3c81a0d3e146238f0c1f1131c56772eef4c078..f789212b695af2e5e25884aae60f9921da14908c 100644
--- a/Kernel/lenard_jones_kernel.c
+++ b/Kernel/lenard_jones_kernel.c
@@ -63,16 +63,16 @@ double lenard_jones_compute(size_t i, size_t j, size_t k, double* origin, double
double pos[3];
for(size_t dim = 0; dim < 3; ++dim)
pos[dim] = basis[dim][0]*i + basis[dim][1]*j + basis[dim][2]*k + origin[dim];
- double d = sqrt(pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2]);
- d = d*d*d*d*d*d;
+ double d = pos[0] * pos[0] + pos[1] * pos[1] + pos[2] * pos[2];
+ d = d*d*d;
if(d > 0.)
{
double e = -(kernel->c6)/d+(kernel->c12)/(d*d);
- return (e>1e4) ? 1e4 : e;
+ return (e>20) ? 20 : e;
}
else
{
- return 1e4;
+ return 20;
}
}
diff --git a/inc/grid_reader.h b/inc/grid_reader.h
index 25a9d2c055b92cc5351938c6c3eb0a6cbe6b5140..9b973947abf6fb3c261c079e3b21f6c092c344fb 100644
--- a/inc/grid_reader.h
+++ b/inc/grid_reader.h
@@ -16,4 +16,5 @@ size_t get_size(void* src);
void* copy_grid_format(void*);
typedef double (*compute_type)(size_t,size_t,size_t,double*,double**, void*);
void set_grid_data(void* src, compute_type compute, void* kernel);
+void* resample_grid_data(void* src);
#endif
diff --git a/main.c b/main.c
index e8c0d7f83ad210dd5124420a74cd2593ce0f9f98..57d593481d457f41df983af2a70e349b563bf2f3 100644
--- a/main.c
+++ b/main.c
@@ -12,8 +12,11 @@ static void inner_product_grid(fftw_complex* dest, fftw_complex* src, size_t nx,
{
re = dest[i][0];
im = dest[i][1];
- dest[i][0] = (re*src[i][0]-im*src[i][1])/V;
- dest[i][1] = (re*src[i][1]+im*src[i][0])/V;
+ size_t x = i / (nz/2+1) / ny;
+ size_t y = i / (nz/2+1) % ny;
+ size_t z = i % (nz/2+1);
+ dest[i][0] = (re*src[i][0]-im*src[i][1])/V*((0.5-(x+y+z)%2)*2);
+ dest[i][1] = (re*src[i][1]+im*src[i][0])/V*((0.5-(x+y+z)%2)*2);
}
}
static void shift_zero_frequency(size_t nx, size_t ny, size_t nz, double* data)
@@ -75,7 +78,11 @@ int main(int argc, char** argv)
char* kernel_name = argv[2];
char* kernel_infile = argv[3];
- void* density = gridformat_reader(density_file);
+ void* density = gridformat_reader(density_file);
+ void* resample_density = resample_grid_data(density);
+ grid_destroy(density);
+ density = resample_density;
+
void* potential = copy_grid_format(density);
create_type creator = find_list(create_map, kernel_name);