Homec4science
Diffusion Computational physics (master) Browse

compphys/exo4/molecule_problem.m
35dbb7a4f6d5master

Last Change
Raw File

molecule_problem.m
View Options

nTwo = false;
nThree = false;
nFour = false;
nFive = false;
nSix = false;
nSeven = false;
nEight = true;
k_lower = 0.5;
k_high = 3.0;
k_highest = 7.0;
global kz
kz = 0;
%%%%%%%%%%%%%
% test 2,3,4, random configuration
if nTwo
E = [];
for N = 2:4
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, []);
%X = solve_nonl_CG_num(@(x) f(x), 3*N, []);
E = [E f(X)]
% remap solution
X =reshape(X, 3, N)'
visualize_molecule(X, 3, N-1, sprintf('output/molecule%d',N))
end
figure
plot(2:4,E)
title("Case 2,3,4")
xlabel("N")
ylabel("E")
end
%%%%%%%%%%
if nThree
N = 5 % , locate minimum
aptempt = 100;
E = zeros(aptempt,1);
%N = 5;
norm_range = [0.1,2.0]
X = cell(aptempt,1);
for k = 1:aptempt
k
hint = rand(3*5,1);
hnorm = rand() * (norm_range(2) - norm_range(1)) + norm_range(1)
hint = hint * hnorm / norm(hint);
X{k} = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, hint);
U = f(X{k})
X{k} = reshape(X{k}, 3, N)';
E(k) = U;
end
cases = input('Case to visualise? ');
for i = 1:length(cases)
visualize_molecule(X{i}, 3, 17 + i, sprintf('output/molecule5_statcase_%d',i))
end
outcols('output/molecule_stat.out', {'case', 'E'}, [(1:aptempt)', E]);
figure
plot(1:aptempt,E)
title("Case N = 5, energy statistics")
xlabel("case number")
ylabel("E")
end
%%%%%%%%%%%%%%%%%
N = 5;
kz = 0.5;
if nFour || nFive || nSix
hint = [];
if nSix
kz = 0.5
hint = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, []);
% start by tetrahedron
end
if ~nFour
kz = 2.0
end
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, hint);
E = f(X)
% remap solution
X =reshape(X, 3, N)'
visualize_molecule(X, 3, 6, 'output/molecule')
end
if nSeven
disp("compute hint (2)")
kz = 0.5;
hint2 = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, []);
fprintf("hint2: E = %g\n", f(hint2))
disp("compute hint (3)")
kz = 2.0;
hint3 = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, []);
fprintf("hint3: E = %g\n", f(hint3))
disp("compute hint (4)")
hint4 = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, hint2);
fprintf("hint4: E = %g\n", f(hint4))
kz = 0;
disp("Set kz = 0, compute results using the previous hints")
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, hint2);
fprintf("hint2: E = %g\n", f(X))
X =reshape(X, 3, N)'
visualize_molecule(X, 3, 7, 'output/molecule_hint2')
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, hint3);
fprintf("hint3: E = %g\n", f(X))
X =reshape(X, 3, N)'
visualize_molecule(X, 3, 8, 'output/molecule_hint3')
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, hint4);
fprintf("hint4: E = %g\n", f(X))
X =reshape(X, 3, N)'
visualize_molecule(X, 3, 9, 'output/molecule_hint4')
end
if nEight
N = 6;
kz = 0.0;
fprintf("N = 6\n\n")
fprintf("Random configuration\n")
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, []);
E = f(X)
X =reshape(X, 3, N)'
visualize_molecule(X, 3, 10, 'output/molecule_rand6')
fprintf("Softly increase energy, kz = %g\n", k_lower)
kz = k_lower;
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, []);
E = f(X)
Y =reshape(X, 3, N)'
visualize_molecule(Y, 3, 11, 'output/molecule_inc6')
fprintf("Decrease energy, kz = 0, set previous as hint\n")
kz = 0;
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, X);
E = f(X)
X =reshape(X, 3, N)'
visualize_molecule(X, 3, 12, 'output/molecule_incdec6')
fprintf("Strongly increase energy, kz = %g\n", k_highest)
kz = k_highest;
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, []);
E = f(X)
Y =reshape(X, 3, N)'
visualize_molecule(Y, 3, 13, 'output/molecule_incstrong6')
Z = X(:); % copy it
fprintf("Decrease energy, kz = %g, set previous as hint\n", k_high)
kz = k_high;
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, X);
E = f(X)
Y =reshape(X, 3, N)'
visualize_molecule(Y, 3, 14, 'output/molecule_incstrongdec6')
fprintf("Decrease energy, kz = %g, set previous as hint\n", k_lower)
kz = k_lower;
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, X);
E = f(X)
Y =reshape(X, 3, N)'
visualize_molecule(Y, 3, 15, 'output/molecule_incstrongdecdec6')
fprintf("Decrease energy level, kz = 0, set previous as hint\n")
kz = 0.0;
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, X);
E = f(X)
Y =reshape(X, 3, N)'
visualize_molecule(Y, 3, 16, 'output/molecule_softfall6')
fprintf("Jump to zero from highest energy level, kz = %g -> 0\n", k_highest)
kz = 0.0;
X = solve_nonl_CG(@(x) f(x), @(i,x) grad_f(i,x), @(i,j,x) hess_f(i,j,x), 3*N, Z);
E = f(X)
Y =reshape(X, 3, N)'
visualize_molecule(Y, 3, 17, 'output/molecule_strongfall6')
end
% problem specific function definition
% scalar field, nothing to say
function y = f(X)
global kz
N = fix(length(X)/ 3);
y = X' * X + kz * X(1:3:N)' * X(1:3:N);
for i = 1:N
for j = (i+1):N
y = y + 1.0 / norm(X(threeV(i)) - X(threeV(j)));
end
end
end
% i is the index w.r.t the gradient of f
% i can be an array
% y is returned as a column vector
function y = grad_f(i, X)
global kz
N = length(X)/ 3;
% compute result as a column vector
y = zeros(length(i),1);
for arg = 1:length(i)
[i_i, l] = threeD(i(arg), 3);
y_y = 2 * X(three(i_i, l));
if l == 3
y_y = 2 * (1.0 + kz) * X(three(i_i, l));
end
for k = 1:N
if k ~= i_i
r = X(threeV(i_i)) - X(threeV(k));
y_y = y_y - (r(l) / norm(r)^3);
end
end
y(arg) = y_y;
end
end
% i,n are the row and resp. column indexes referring
% to the hessian matrix
% i,n can be arrays
function y = hess_f(i, n, X)
global kz
N = length(X) / 3;
y = zeros(length(i),length(n));
for arg_i = 1:length(i)
for arg_n = 1:length(n)
[i_i, l] = threeD(i(arg_i), 3);
[n_n, m] = threeD(n(arg_n), 3);
if i_i == n_n
y_y = 0;
for k = 1:N
if k ~= i_i
r = X(threeV(i_i)) - X(threeV(k));
y_y = y_y + 3 * r(m) * r(l) / norm(r)^5;
if l == m
y_y = y_y - 1.0 / norm(r)^3;
end
end
end
if l == m
y_y = y_y + 2 * (1.0 + (l == 3) * kz);
end
else
r = X(threeV(i_i)) - X(threeV(n_n));
y_y = -3 * r(m) * r(l) / norm(r)^5;
if l == m
y_y = y_y + 1.0 / norm(r)^3;
end
end % if
y(arg_i, arg_n) = y_y;
end % for
end % for
end
function [i,j] = threeD(k, m)
i = fix((k-1) / m) + 1;
j = mod((k-1), m) + 1;
% idea: k = m*(i-1) + j
end
function I = threeV(k)
I = 3*(k-1)+1:3*k;
end
% j = 0,1,2
function k = three(i,j)
k = 3 * (i-1) + j;
end
c4science · Help