xxxxxxxxxx268
let u = [];let f = [];let K = [];let P = [];let w = [];let charge = [];let r = 0;let E =0;let c=[];let norm = [];let N=50;let N2=25;let D=12;let dt;let h;let t=1;let d=0.05;let M=11;let M2=16;function setup() { createCanvas(400, 400); h=8/N; dt=d*pow(h,2); for (var m=1;m<M2;m++){ w[m]=[]; u[m]=[]; c[m]=1; charge[m]=1; P[m]=[]; norm[m]=0; } for (var m=1;m<M2;m++){ for (var i=0;i<N+1;i++){ w[m][i]=[]; u[m][i]=[]; K[i]=[]; P[m][i]=[]; for (var j=0;j<N+1;j++){ w[m][i][j]=[]; u[m][i][j]=[]; K[i][j]=[]; P[m][i][j]=[]; for (var k=0;k<N+1;k++){ w[m][i][j][k]=0; u[m][i][j][k]=0; K[i][j][k]=0; P[m][i][j][k]=0; } } } } for (var i=0;i<N+1;i++){ for (var j=0;j<N+1;j++){ for (var k=0;k<N+1;k++){ r=pow((i-N2+D)*h,2)+pow((j-N2)*h,2)+pow((k-N2)*h,2); r=sqrt(r); if (r<1 & i<26-D){ u[1][i][j][k]=exp(-5*r); w[1][i][j][k]=1; } if (r<1 & i>24-D){ u[2][i][j][k]=exp(-5*r); w[2][i][j][k]=1; } if (r>1.5 & r<3.5 & i>(30-D) & j<24){ u[3][i][j][k]=exp(-3*r); w[3][i][j][k]=1; } if (r>1.5 & r<3.5 & i<(28-D)){ u[5][i][j][k]=exp(-3*r); w[5][i][j][k]=1; } if (r>1.5 & r<3.5 & i>(30-D) & j>26){ u[4][i][j][k]=exp(-3*r); w[4][i][j][k]=1; } if (r>1.5*h){ K[i][j][k]=5/r; } else{ K[i][j][k]=5/(1.5*h); } r=pow((i-N2-D)*h,2)+pow((j-N2)*h,2)+pow((k-N2)*h,2); r=sqrt(r); if (r<1 & i<26+D){ u[6][i][j][k]=exp(-5*r); w[6][i][j][k]=1; } if (r<1 & i>24+D){ u[7][i][j][k]=exp(-5*r); w[7][i][j][k]=1; } if (r>1.5 & r<3.5 & i>(22+D)){ u[8][i][j][k]=exp(-3*r); w[8][i][j][k]=1; } if (r>1.5 & r<3.5 & i<(20+D) & j<24){ u[9][i][j][k]=exp(-3*r); w[9][i][j][k]=1; } if (r>1.5 & r<3.5 & i<(20+D) & j>26){ u[10][i][j][k]=exp(-3*r); w[10][i][j][k]=1; } if (r>1.5*h){ K[i][j][k]=K[i][j][k]+5/r; } else{ K[i][j][k]=K[i][j][k]+5/(1.5*h); } } } }}function draw() { background(220); noStroke(); if (t<200){ t=t+1; E=0; for (var i=1;i<N;i++){ for (var j=1;j<N;j++){ for (var k=1;k<N;k++){ for (var m=1;m<M;m++){ c[m]=0; for (var n=1;n<M;n++){ if (n!=m){ c[m]=c[m]-u[n][i][j][k]; } c[m]=0.5*(c[m]+u[m][i][j][k]); } }for (var m=1;m<M;m++){//Front tracking of electron characteristic functionsif (m<6){ r=pow((i-N2+D)*h,2)+pow((j-N2)*h,2)+pow((k-N2)*h,2)+h*h;r=sqrt(r);if (r>0.5){w[m][i][j][k]=w[m][i][j][k]+2*dt*abs(c[m])*(w[m][i+1][j][k]+w[m][i-1][j][k]-6*w[m][i][j][k]+w[m][i][j+1][k]+w[m][i][j-1][k]+w[m][i][j][k+1]+w[m][i][j][k-1])/pow(h,2)+ 10*dt*c[m]*sqrt(pow((w[m][i+1][j][k]-w[m][i-1][j][k])/h,2)+pow((w[m][i][j+1][k]-w[m][i][j-1][k])/h,2)+pow((w[m][i][j][k+1]-w[m][i][j][k-1])/h,2));}}if (m>5){ r=pow((i-N2-D)*h,2)+pow((j-N2)*h,2)+pow((k-N2)*h,2)+h*h;r=sqrt(r);if (r>0.5){w[m][i][j][k]=w[m][i][j][k]+2*dt*abs(c[m])*(w[m][i+1][j][k]+w[m][i-1][j][k]-6*w[m][i][j][k]+w[m][i][j+1][k]+w[m][i][j-1][k]+w[m][i][j][k+1]+w[m][i][j][k-1])/pow(h,2)+ 10*dt*c[m]*sqrt(pow((w[m][i+1][j][k]-w[m][i-1][j][k])/h,2)+pow((w[m][i][j+1][k]-w[m][i][j-1][k])/h,2)+pow((w[m][i][j][k+1]-w[m][i][j][k-1])/h,2));}}//Update of electron density functions////Hom Neumann problem for each electron// u[m][i][j][k]=u[m][i][j][k]+0.5*d*((u[m][i+1][j][k]-u[m][i][j][k])*(w[m][i+1][j][k]+w[m][i][j][k])/2-(u[m][i][j][k]-u[m][i-1][j][k])*(w[m][i][j][k]+w[m][i-1][j][k])/2)+0.5*d*((u[m][i][j+1][k]-u[m][i][j][k])*(w[m][i][j+1][k]+w[m][i][j][k])/2-(u[m][i][j][k]-u[m][i][j-1][k])*(w[m][i][j][k]+w[m][i][j-1][k])/2)+0.5*d*((u[m][i][j][k+1]-u[m][i][j][k])*(w[m][i][j][k+1]+w[m][i][j][k])/2-(u[m][i][j][k]-u[m][i][j][k-1])*(w[m][i][j][k]+w[m][i][j][k-1])/2)+dt*(K[i][j][k]-2*P[m][i][j][k])*u[m][i][j][k]*w[m][i][j][k]; //Compute energy// if (w[m][i][j][k]>0.5){E = E+0.5*(pow(u[m][i+1][j][k]-u[m][i][j][k],2)+pow(u[m][i][j+1][k]-u[m][i][j][k],2)+pow(u[m][i][j][k+1]-u[m][i][j][k],2))*h;}E=E+(P[m][i][j][k]-K[i][j][k])*pow(u[m][i][j][k],2)*pow(h,3); //Solve Poisson equation for electron potentials// c[m]=0; for (var n=1;n<M;n++){ if (n!=m){ c[m]=c[m]+pow(u[n][i][j][k],2); } } P[m][i][j][k]=P[m][i][j][k]+dt*(P[m][i+1][j][k]-6*P[m][i][j][k]+P[m][i-1][j][k]+P[m][i][j+1][k]+P[m][i][j-1][k]+P[m][i][j][k+1]+P[m][i][j][k-1])/pow(h,2)+2*PI*dt*c[m];}}}}} //Normalisation of charge densities// for (var m=1;m<M;m++){ norm[m] = 0; for (var i=1;i<N;i++){ for (var j=1;j<N;j++){ for (var k=1;k<N;k++){ norm[m] = norm[m] + pow(u[m][i][j][k],2)*pow(h,3); } } } for (var i=1;i<N;i++){ for (var j=1;j<N;j++){ for (var k=1;k<N;k++){ u[m][i][j][k]=u[m][i][j][k]/sqrt(norm[m]); } } } } //Plots// // if (t/10 - floor(t/10)>0.8){ for (var m=1;m<M;m++){ for (var i=1;i<N;i++){ for (var j=1;j<N;j++){ fill(255,0,0,500*u[m][i][j][25]); square(8*i,8*j,8); } } }// } fill(255,255,255,255); fill(0); circle(200-D*8,200,10); circle(200+D*8,200,10); for (var i=0;i<N;i++){ for (var m=1;m<M;m++){ fill(0); ellipse(8*i,340-100*w[m][i][30][25],6); fill(255,255,0,255); ellipse(8*i,340-100*u[m][i][30][25],6); fill(0,255,255,255); ellipse(8*i,340-400*P[m][i][30][25],4); // ellipse(8*i,300-100*P[2][i][25][25],3);// ellipse(8*i,300-100*P[3][i][25][25],6);// ellipse(8*i,300-100*P[4][i][25][25],3); } fill(255,0,255,255); ellipse(8*i,160,4); fill(0,0,255,255); ellipse(8*i,300-10*K[i][25][25],6); } fill(0); text(t,100,378); text(E+25/(2*D*h)-0*75.4,130,395); text("time step =",10,378); text("Total Energy =",10,395); text("(ref -49.2 Hartree)",260,395); text("Boron2 Molecule 3d 50^3 B: 1st shell 2 elect, 2nd 3 elect",5,20); text("u = sum_i u_i, u_i non-overlap elect density i=1,..,12",5,40); text("w_i characteristic function for support of u_i, initiate with gap 2nd shell",5,60); text("Minimise E = sum_i integral (0.5*(w_i*|grad(u_i)|^2-K*u_i^2+P_i*u_i^2)dx",5,80); text("over u_i with continuity of u across free bdy updated by front track w_i",5,100); text("H_i = - 0.5*Laplacian - K(x) + 2*P_i(x)",5,140); text("K(x) kernel potential, P_i(x) electron potentials",5,120); text("P_1(x) = integral dy*u_2(y)^2/(2*|x-y|) acting on u_1 et cet",5,180); text("Compute E by time stepping du_i/dt = -H_iu_i + normalisation",5,160); text("Compute P_1 by solving -Laplacian P_1 =2*PI*(u_2^2) et cet",5,220);text("yellow: u, blue:kernel pot, lightblue: elec pot, black: w, cross cut magenta,",10,360);text("Update w_i by front track: dw_i/dt+jump*abs(grad w_i)",5,240);}