It is a 1D MSM demo code to check the force error calculations with the expected values. In this minimal code, there is no restriction and prolongation steps as it is one level MSM.
I first get rid of old variables and figures.
clear variables clc close all
Setting random seed to recreate the results.
s = RandStream('mt19937ar','Seed',7); RandStream.setGlobalStream(s);
Grid spacing is set to h=2.5 to mimic large water sphere.
h = 2.5;
The number of particles.
The distribution of the particles on a line. I choose the boundaries -50 and 50 and the n=41 to mimic the large water sphere data in the NAMD-lite. In that data, there are 30006 atoms in a sphere with radius r=44A. The number of atoms in each direction is approximately 2*r/((4/3*pi*r^3/30006)^(1/3)) = 38. I set n=41 in order to locate each particles on the grid with h=2.5. In this way, I am sure that there is no error in anterpolation.
x0 = linspace(-50,50,n); %x0 = sort(100*rand(n,1)-50);
I take all the masses one. I don't know if this effects my calculations.
q0 = ones(n,1);
Since I know h=2.5 and the boundaries of [-50 50], I can preallocate the finest mesh beforehand. I enlarge the boundary by 2*h on each side.
x1 = linspace(-55,55,45);
The cutoff values that I try. In the dissertation, a=6:20 is used almost always. But here I want to see what happens outside of that safe region.
arange = [1:10 20:10:100 200:100:1000];
This is the C3 Taylor smoothing and its derivative.
s = @(p) 1 + (p*p-1)*(-1/2 + (p*p-1)*(3/8 + (p*p-1)*(-5 /16))); sd = @(p) (2*p) *(-1/2 + (p*p-1)*(3/4 + (p*p-1)*(-15/16)));
The cubic "numerical" Hermite nodal base and its derivative. The enumeration of the pieces of the functions is as follows:
<--------|--------|--------|--------|--------|-------> -2 -1 0 1 2 c1 c2 c3 c4 c1d c2d c3d c4d --> derivatives
To make sure, I don't make any mistake, I took the functions definitions directlt from NAMD-Lite. However, in NAMD-Lite the enumeration of the pieces are little different.
Enumeration of the cubic pieces in the NAMD-Lite:
<--------|--------|--------|--------|--------|-------> -2 -1 0 1 2 phi3 phi2 phi1 phi0 dphi3 dphi2 dphi1 dphi0 --> derivatives
c1 = @(t) 0.5 * (1 + t) * (2 + t) * (2 + t); % phi3 in NAMD-Lite c1d = @(t) (1.5 * t + 2) * (2 + t); % dphi3 in NAMD-Lite c2 = @(t) (1 + t) * (1 - t - 1.5 * t * t); % phi2 in NAMD-Lite c2d = @(t) (-5 - 4.5 * t) * t; % dphi2 in NAMD-Lite c3 = @(t) (1 - t) * (1 + t - 1.5 * t * t); % phi1 in NAMD-Lite c3d = @(t) (-5 + 4.5 * t) * t ; % dphi1 in NAMD-Lite c4 = @(t) 0.5 * (1 - t) * (2 - t) * (2 - t); % phi0 in NAMD-Lite c4d = @(t) (1.5 * t - 2) * (2 - t) ; % dphi0 in NAMD-Lite
In the following double loop, I create the original kernel matrix and calculate the exact forces.
G = zeros(n,n); freal = zeros(n,1); for i=1:n-1 for j=i+1:n rij = x0(i)-x0(j); r = abs(rij); r_1 = 1/r; rij_r3 = q0(i)*q0(j)*rij/r^3; G(i,j) = r_1; G(j,i) = G(i,j); freal(i) = freal(i) - rij_r3; freal(j) = freal(j) + rij_r3; end end
Here is the exact potential energy.
ureal = 0.5*q0'*G*q0;
For each cutoff value, MSM is run to calculate the approximate force and potentials.
Here I prefer to create the matrices just to check the linear algebra side of the MSM match with the loop calculations.
We split the G matrix into two parts:
Ghat = zeros(n,n); Gbar = zeros(n,n); fshort = zeros(n,1);
I calculate the short-range potential in two ways: using the matrix algebra and accumulating the potentials in the loop. That's why I use two different variables: ushort and ushortloop. In the end, I will compare them for sanity check.
ushortloop = 0; for i=1:n qi = q0(i);
Gbar(i,i) = s(0)/a; Ghat(i,i) = -Gbar(i,i);
ushortloop = ushortloop + qi*qi*0.5*Ghat(i,i); for j=i+1:n qj = q0(j); rij = x0(i)-x0(j); r = abs(rij); r_1 = 1/r; rij_r3 = rij/r^3; % Long-range if r >= a Gbar(i,j) = r_1; Gbar(j,i) = r_1; Ghat(i,j) = 0; Ghat(j,i) = 0; else % Short range g = s(r/a)/a; Gbar(i,j) = g; Gbar(j,i) = g; Ghat(i,j) = r_1-g; Ghat(j,i) = r_1-g; force = qi*qj*(1/r^2 + 1/a^2*sd(r/a))*rij/r; fshort(i) = fshort(i) - force; fshort(j) = fshort(j) + force; ushortloop = ushortloop + qi*qj*(r_1 - g); end end end
short-range potential energy is directly calculated by using
ushort = 0.5*q0'*Ghat*q0;
Here I use similar strategy. I calculate the anterpolation matrix explicitly to check if the grid masses (q1) calculated by q1=A*q0 is equal to the grid masses (q1loop) calculated in the loop. Later, I will use A' in the interpolation step.
nx1 = numel(x1); A = zeros(nx1,n); q1loop = zeros(nx1,1); for i=1:n i0 = floor((x0(i)-x1(1))/h); t = (x0(i) - x1(i0))/h; A(i0 ,i) = c4(t); A(i0+1,i) = c3(t-1); A(i0+2,i) = c2(t-2); A(i0+3,i) = c1(t-3); q1loop(i0 ) = q1loop(i0 ) + c4(t )*q0(i); q1loop(i0+1) = q1loop(i0+1) + c3(t-1)*q0(i); q1loop(i0+2) = q1loop(i0+2) + c2(t-2)*q0(i); q1loop(i0+3) = q1loop(i0+3) + c1(t-3)*q0(i); end q1 = A*q0;
The matrix is used to calculate the pairwise interactions between the grid masses q1.
G1 = zeros(nx1,nx1); for i=1:nx1 G1(i,i) = s(0)/a; for j=i+1:nx1 r = abs(x1(j)-x1(i)); if r >= a G1(i,j) = 1/r; else G1(i,j) = s(r/a)/a; end G1(j,i)=G1(i,j); end end
Grid potential e1 is just the matrix-vector product of G1 and q1.
e1 = G1*q1;
Since we already have the transformation matrix A in the anterpolation, I can use its transpose in the interpolation to get the long-range potentials in the particle level.
e0long = A'*e1;
I'll calculate the same e0long by using the loops to make sure I do the interpolation right.
e0longloop = zeros(n,1);
In the loop, I will calculate the long-range forces.
flong = zeros(n,1);
For each particle ...
for i=1:n % Location of the particle x = x0(i); % lowest grid location having non-zero contribution. i0 = floor((x-x1(1))/h); % Coordinate of the particle in terms of h. t = (x - x1(i0))/h;
Sanity check: Since I'm using cubic bases, t should be between [1 2].
if t < 1 || t > 2 error('t should be between [1 2]'); end
flong(i) = e1(i0 ) * c4d(t )/h + ... e1(i0+1) * c3d(t-1)/h + ... e1(i0+2) * c2d(t-2)/h + ... e1(i0+3) * c1d(t-3)/h; e0longloop(i) = e1(i0 ) * c4(t ) + ... e1(i0+1) * c3(t-1) + ... e1(i0+2) * c2(t-2) + ... e1(i0+3) * c1(t-3); end
Long-range potential energy is
ulong = 0.5*q1'*e1;
Finally, I add the short and long-range energy and forces to obtain the final approximations.
umsm = ushort+ulong; fmsm = fshort+flong;
Relative error in average mass-weighted force.
errfunc = @(f,m) sqrt(mean((f.^2)./m)); ferror(ia) = errfunc(freal-fmsm,q0)/errfunc(freal,q0);
Relative error in potential energy.
uerror(ia) = abs(ureal-umsm)/ureal;
Short-range potential energies calculated two different ways should be equal to each other.
if abs(ushort-ushortloop)/abs(ushort) > 1E-10 error('ushort must be equal to ushortloop\n'); end
Long-range potentials calculated by using both ways should be equal to each other.
if norm(e0long-e0longloop)/norm(e0long) > 1E-10 error('e0 should be equal to e0loop'); end
The sum of grid masses should be equal to the total mass.
if abs(sum(q0)-sum(q1))/sum(q0) > 1E-10 error('sum(q0) should be equal to sum(q1)'); end
The grid masses calculated both ways should be equal
if norm(q1-q1loop)/norm(q1) > 1E-10 error('q1 should be equal to q1loop'); end
The long-range interaction kernel is approaximated with the following relation.
If the particles are located on the grid nodes, then the approximation should be exact.
Gtilde = A'*G1*A; if norm(Gbar-Gtilde)/norm(Gbar) > 1E-10 error('The error in the approximation is high'); end
Long-range potential energy can be calculated in two different ways. They should be equal.
if abs(ulong-0.5*q0'*e0long) > 1E-10 error('0.5*q0T*e0 shoulb be equal to ulong'); end
The percent relative error in average force is displayed in loglog plot. Here, I also plot two lines of slopes, -2.5 and -5, the later is for asymptotic behaviour. The expected slopes are -4 and -5 since I'm using cubic (p=3) interpolating polynomials. But the observations in the dissertation and my experiments with NAMD-Lite indicates that we should observe -3 and -5. However, my code produces a slope -2.5. I'm still not sure why I can't get a slope of -3.
loglog(arange,ferror,'b-*'); hold on loglog(arange,(ferror(10)*arange(10)^2.5)*arange.^-2.5,'k-'); loglog(arange,(ferror(20)*arange(20)^5)*arange.^-5,'r-'); title('Cubic Numerical Hermite & C3 Taylor','FontSize',14); xlabel('cutoff (a)','FontSize',14); ylabel('percent relative error in average force','FontSize',14); h=legend('Force','Slope = -2.5','Slope = -5'); set(h,'FontSize',14); grid on set(gcf,'color','w'); export_fig('matlab_grid_Cubic_C3Taylor.pdf');