% Load data
clear all;
close all;
load Tests/cubeions.mat

% Compute spline coefficients
Lcut = 8.0;
Dtol = 1.0e-5;
ewcoeff = getalpha(Lcut, Dtol);
ng = [64 64 64];
[SPx, dSPx, qmapx] = splcof(crdq, ng, gdim, 4, 1);
[SPy, dSPy, qmapy] = splcof(crdq, ng, gdim, 4, 2);
[SPz, dSPz, qmapz] = splcof(crdq, ng, gdim, 4, 3);

% Compute the charge density mesh
Q = MapCharges(ng, SPx, SPy, SPz, qmapx, qmapy, qmapz);
colmap = zeros(256, 3);
for i = 1:1:128
  colmap(i,:) = [ 1.0 i/128.0 i/128.0 ];
  colmap(i+128,:) = [ (128-i)/128.0 (128-i)/128.0 1.0 ];
end
figure;
hold on;
colormap(colmap);
surf(Q(:,:,32), 'edgecolor', 'none');
daspect([1 1 1]);
caxis([-332.0636/4.0, 332.0636/4.0]);

% Compute the long-ranged, reciprocal-space energy
[U,Enrg,B,C] = pmeOrthoConvBC(Q, gdim, [4 4 4], 0.5/ewcoeff);
colmap = zeros(256, 3);
for i = 1:1:128
  colmap(i,:) = [ 1.0 (64.0 + 0.5*i)/128.0 i/128.0 ];
  colmap(i+128,:) = [ (128-0.4*i)/128.0 (128-i)/128.0 1.0 ];
end
figure;
hold on;
colormap(colmap);
surf(real(U(:,:,32)), 'edgecolor', 'none');
daspect([1 1 1]);
caxis([-400.0, 700.0]);