% Computes PCS using different models in basic Cu(II) porphyrin complex. See 
% the "getting started" manual at
%
% http://spindynamics.org/wiki/index.php?title=Pseudocontact_shift_analysis
%
% The paper describing the distributed PCS model used below is available at
%
%                   http://dx.doi.org/10.1039/c6cp05437d
%
% e.suturina@soton.ac.uk

function Eu_III_cube_export()

% proton coordinates (nuclei at which PCS is to be evaluated)
nxyz=[  5.272196	17.426539	-8.878158
2.985248	18.364758	-8.936228
1.16006	17.255536	-7.684605
1.6354	15.217199	-6.362585
3.930697	14.300749	-6.285785
4.830619	14.569099	-9.638303
6.695498	12.159419	-9.660988
5.13315	12.245408	-10.501389
5.201434	12.305951	-8.729255
6.5604	15.673716	-11.035515
5.884525	14.253886	-11.856142
7.477453	14.16297	-11.091171
5.578765	16.659819	-5.489822
6.01112	14.701326	-4.027065
6.479157	16.213977	-3.22715
7.71599	15.157628	-3.922625
8.590744	17.107785	-5.385064
7.382922	18.143029	-4.59585
7.436444	18.019294	-6.364837
7.843654	9.570401	-8.993885
8.188288	7.122235	-9.022411
9.938082	6.101716	-7.599405
11.335113	7.5378	-6.145315
10.971346	9.98306	-6.106039
8.04549	10.204739	-5.567171
8.18618	12.793981	-3.968163
9.426493	11.534247	-3.99241
7.827289	11.193563	-3.304466
6.159965	12.598487	-5.570391
5.805446	11.027263	-4.821112
6.009836	11.15554	-6.579051
10.576409	10.636895	-9.513224
9.702744	13.130919	-11.029707
8.863657	11.574687	-11.047022
10.489798	11.71495	-11.741673
11.706785	13.465944	-9.43338
12.488598	12.079154	-10.221801
12.262452	12.098594	-8.461922
11.160212	18.221664	-5.893204
13.057208	19.804282	-5.790816
15.101867	19.373175	-7.118154
15.233808	17.364359	-8.558924
13.321394	15.801382	-8.678987
12.424682	15.613232	-5.284285
11.330942	12.770227	-5.359897
12.803973	13.270321	-4.502723
12.718897	13.373115	-6.272211
10.458585	16.1245	-3.857633
11.514061	14.926955	-3.084297
10.010844	14.414102	-3.864163
11.055965	17.561877	-9.351343
11.210235	15.628737	-10.900835
10.301822	16.965295	-11.632888
9.445889	15.558222	-10.986828
8.046244	17.096152	-9.446737
8.884769	18.481881	-10.175077
8.885751	18.259662	-8.414927];
   
% coordinates of all atoms in the molecule (latest Euler-angle corrected geometry)  
xyz=[ 8.759911	14.200605	-7.478074
6.637987	14.584815	-8.573322
6.903338	15.367166	-6.475679
8.573545	12.007282	-6.502507
6.127544	15.239771	-7.542068
4.745006	15.805334	-7.578013
4.472936	16.947779	-8.323391
5.272196	17.426539	-8.878158
3.188224	17.471319	-8.357317
2.985248	18.364758	-8.936228
2.164371	16.849731	-7.654153
1.16006	17.255536	-7.684605
2.431268	15.706563	-6.91185
1.6354	15.217199	-6.362585
3.718352	15.189536	-6.869629
3.930697	14.300749	-6.285785
5.837385	14.138148	-9.691206
4.830619	14.569099	-9.638303
5.702953	12.619481	-9.646932
6.695498	12.159419	-9.660988
5.13315	12.245408	-10.501389
5.201434	12.305951	-8.729255
6.473701	14.586093	-10.998197
6.5604	15.673716	-11.035515
5.884525	14.253886	-11.856142
7.477453	14.16297	-11.091171
6.600713	16.275087	-5.39152
5.578765	16.659819	-5.489822
6.703176	15.544855	-4.061158
6.01112	14.701326	-4.027065
6.479157	16.213977	-3.22715
7.71599	15.157628	-3.922625
7.55691	17.462671	-5.433373
8.590744	17.107785	-5.385064
7.382922	18.143029	-4.59585
7.436444	18.019294	-6.364837
9.169225	11.407889	-7.523364
9.384132	9.929836	-7.546484
8.605581	9.120431	-8.367187
7.843654	9.570401	-8.993885
8.801789	7.746778	-8.383593
8.188288	7.122235	-9.022411
9.783271	7.174121	-7.584987
9.938082	6.101716	-7.599405
10.566504	7.979755	-6.768563
11.335113	7.5378	-6.145315
10.364991	9.352507	-6.746738
10.971346	9.98306	-6.106039
7.865555	11.281849	-5.471214
8.04549	10.204739	-5.567171
8.356556	11.722481	-4.100467
8.18618	12.793981	-3.968163
9.426493	11.534247	-3.99241
7.827289	11.193563	-3.304466
6.366793	11.525138	-5.615743
6.159965	12.598487	-5.570391
5.805446	11.027263	-4.821112
6.009836	11.15554	-6.579051
9.565338	12.184639	-8.521114
10.44302	11.723367	-9.573482
10.576409	10.636895	-9.513224
9.839183	12.050814	-10.930873
9.702744	13.130919	-11.029707
8.863657	11.574687	-11.047022
10.489798	11.71495	-11.741673
11.812145	12.377012	-9.416723
11.706785	13.465944	-9.43338
12.488598	12.079154	-10.221801
12.262452	12.098594	-8.461922
10.683828	15.144191	-6.35799
10.189241	15.891031	-8.426696
10.972118	15.960613	-7.359422
12.121853	16.91226	-7.292191
12.053755	18.039406	-6.479848
11.160212	18.221664	-5.893204
13.120338	18.925203	-6.421513
13.057208	19.804282	-5.790816
14.267136	18.683916	-7.166985
15.101867	19.373175	-7.118154
14.340826	17.557257	-7.975889
15.233808	17.364359	-8.558924
13.269684	16.677593	-8.042262
13.321394	15.801382	-8.678987
11.584726	14.909324	-5.251983
12.424682	15.613232	-5.284285
12.149668	13.495857	-5.348707
11.330942	12.770227	-5.359897
12.803973	13.270321	-4.502723
12.718897	13.373115	-6.272211
10.851909	15.108935	-3.933982
10.458585	16.1245	-3.857633
11.514061	14.926955	-3.084297
10.010844	14.414102	-3.864163
10.194304	16.894269	-9.468325
11.055965	17.561877	-9.351343
10.298215	16.224503	-10.829943
11.210235	15.628737	-10.900835
10.301822	16.965295	-11.632888
9.445889	15.558222	-10.986828
8.928135	17.739634	-9.373957
8.046244	17.096152	-9.446737
8.884769	18.481881	-10.175077
8.885751	18.259662	-8.414927];

% mxyz=[ -0.01289	0.001904	0.348096];

% Put the metal at the origin
% xyz=xyz-kron(mxyz,ones(size(xyz,1),1));

%possible g-Tensor values:
%ground state (from effective Hamiltonian):                         0.809329     0.952159     0.978804
%first 3 SO-states after ground state (pseudospin Hamiltonian):     1.66609838   1.64718265   1.05881485

% Eu(III) g-tensor eigenvalues
g_eu=diag([0.809722     0.953096     0.977340]);

% Curie susceptibility tensor
chi_eu=g2chi(g_eu,273,3);

% CAS susceptibility (Orca) in cm³*K/mol at 273 K
% chi_eu=[ 1.409858    -0.004367     0.020120 
%         -0.004367     1.414043     0.001456 
%          0.020120     0.001456     1.796160 ];
% 
% chi_eu=(1.660539067*chi_eu)/273;   
% [~,chi_eu]=eig(chi_eu);

% Parse ORCA log
% props=oparse('U_III.out');

% Extract hyperfine tensors
% hfcs=props.hfc.full.matrix(26:37);

% Compute HFC PCS
% hfc_pcs_cu=zeros(numel(hfcs),1);
% for n=1:numel(hfcs)
%     hfc_pcs_cu(n)=hfc2pcs(hfcs{n},chi_cu,'1H',1);
% end

% Parse ORCA cube and pad the density with zeros to avoid PBC effects 
pad_size=2; [density,ext]=ocparse('Eu_III.spindens.3d',pad_size);

% Compute extentens for original density without padding
zoom_vol=[pad_size/(2*pad_size+1) (pad_size+1)/(2*pad_size+1)...
          pad_size/(2*pad_size+1) (pad_size+1)/(2*pad_size+1)...
          pad_size/(2*pad_size+1) (pad_size+1)/(2*pad_size+1)];

% Draw the probability density schematic
[zoom_den,zoom_ext]=zoom_3d(density,ext,zoom_vol);
figure(); volplot(zoom_den,zoom_ext,[0.05 0.05]);
hold on; molplot(xyz,[]); drawnow();

% Solve Kuprov equation
[pcs_vals,pcs_cube]=kpcs(density,chi_eu,ext,nxyz,'fft');

% Compare HFC PCS with PDE PCS
% disp('Pseudocontact shifts [hfc, pde], ppm');
% disp([hfc_pcs_cu pde_pcs_cu]);

% Plot the PCS field schematic
[zoom_pcs,zoom_ext]=zoom_3d(pcs_cube,ext,zoom_vol);
figure(); volplot(zoom_pcs,zoom_ext,[0.02 0.02]); 
hold on; molplot(xyz,[]);

%atom numbers according to xyz order
atoms=[ 63
7
7
7
6
6
6
1
6
1
6
1
6
1
6
1
6
1
6
1
1
6
1
1
1
6
1
6
1
1
1
6
1
1
1
6
6
6
1
6
1
6
1
6
1
6
1
6
1
6
1
1
6
1
1
1
6
6
1
6
1
6
1
6
1
6
1
6
6
6
1
6
1
6
1
6
1
6
1
1
6
1
1
1
6
1
6
1
1
1
6
1
1
1 ];
     
num_voxels = size(pcs_cube);
dim=num_voxels(1);

% Put the metal at the origin
% xyz=xyz-kron(mxyz,ones(size(xyz,1),1));

% Bohr radius in Angstroms
  rB=1.889725989;

% % Grid
% ext=ext*rB;
% x=linspace(ext(1),ext(2),dim);
% y=linspace(ext(3),ext(4),dim);
% z=linspace(ext(5),ext(6),dim);
% % Step
% xstep=x(2)-x(1);
% ystep=y(2)-y(1);
% zstep=z(2)-z(1);
% % writing a file
% filename='pcs_cube.cub';
% fileID = fopen(filename,'w');
% fprintf(fileID,"PCS cube\n");
% fprintf(fileID,"%s\n",filename);
% fprintf(fileID,"%d   % -E   % -E   % -E   \n",numel(atoms),ext(1),ext(3),ext(5));
% fprintf(fileID,"%d   %6.2f   0   0   \n",dim,xstep);
% fprintf(fileID,"%d   0   %6.2f   0   \n",dim,ystep);
% fprintf(fileID,"%d   0   0   %6.2f   \n",dim,zstep);
% A=[atoms atoms xyz*rB];
% fprintf(fileID,'%d %d %12.8f %12.8f %12.8f\n',A.');
% for i=1:numel(x)
%       for j=1:numel(y)
%          for k=1:numel(z)
%              pcs=pcs_cube(i, j, k)*rB;
%              fprintf(fileID,"% -E   ",pcs);
%             if mod(k,6) == 0
%                fprintf(fileID,"\n");
%             end
%          end
%          fprintf(fileID,"\n");
%       end
% end
% fclose(fileID);

pcs_vals

format long
chi_eu

[iso,ax,rh]=mat2axrh(chi_eu);

ax_m = ax*1e-30
rh_m = rh*1e-30
iso


end