p5.js Web Editor | Na 2+4+4 +1 100

xxxxxxxxxx
let u = [];
let K = [];
let P = [];
let w = [];
let charge = [];
let r = 0;
let E =0;
let c=[];
let norm = [];
let N=100;
let N2=50;
let dt;
let dt1;
let h;
let t=1;
let d=0.1;
let M=12
let M2=12;
​
​
function setup() {
  createCanvas(400, 400);
  h=8/N;
  dt=d*pow(h,2);
  dt1=dt/10;
  for (var m=1;m<M2;m++){
  w[m]=[];
  u[m]=[];
  P[m]=[];
  c[m]=1;
  charge[m]=1;
  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)*h,2)+pow((j-N2)*h,2)+pow((k-N2)*h,2)+0.18*h*h;
   r=sqrt(r);
   if (r<0.7 & i<N2){
   u[1][i][j][k]=exp(-5*r);
   w[1][i][j][k]=1;
   }
   if (r<0.7 & i>N2){
   u[2][i][j][k]=exp(-5*r);
   w[2][i][j][k]=1;
   }
   if (r>0.7 & r<1.5 & i<N2 & j<N2){
   u[3][i][j][k]=exp(-2*r);
   w[3][i][j][k]=1;
   }
   if (r>0.7 & r<1.5 & i<N2 & j>N2){
   u[4][i][j][k]=exp(-2*r);
   w[4][i][j][k]=1;
   }
    if (r>0.7 & r<1.5 & i>N2 & j<N2){
   u[5][i][j][k]=exp(-2*r);
   w[5][i][j][k]=1;
   }
   if (r>0.7 & r<1.5 & i>N2 & j>N2){
   u[6][i][j][k]=exp(-2*r);
   w[6][i][j][k]=1;
   }
   if (r>1.5 & r<2.3 & i<N2 & j<N2){
   u[7][i][j][k]=exp(-2*r);
   w[7][i][j][k]=1;
   }
   if (r>1.5 & r<2.3 & i<N2 & j>N2){
   u[8][i][j][k]=exp(-2*r);
   w[8][i][j][k]=1;
   }
   if (r>1.5 & r<2.3 & i>N2 & j<N2){
   u[9][i][j][k]=exp(-2*r);
   w[9][i][j][k]=1;
   }
   if (r>1.5 & r<2.3 & i>N2 & j>N2){
   u[10][i][j][k]=exp(-2*r);
   w[10][i][j][k]=1;
   }
   if (r>2.3){
   u[11][i][j][k]=exp(-r);
   w[11][i][j][k]=1;
   }
  K[i][j][k]=11/r;
  
  for (var m=3;m<M;m++){
  P[m][i][j][k]=1/r;
  }
  for (var m=7;m<M;m++){
  P[m][i][j][k]=P[m][i][j][k]+2/r;
  }
  for (var m=1;m<3;m++){
  if (r<0.9){
  P[m][i][j][k]=2/0.9;
  }
  if (r>0.9){
  P[m][i][j][k]=2/r;
  }
  }
  for (var m=1;m<7;m++){
  if (r<1.6){
  P[m][i][j][k]=P[m][i][j][k]+2.5/1.6;
  }
  if (r>1.6){
  P[m][i][j][k]=P[m][i][j][k]+2.5/r;
  }
  }
  }
  }
  }
}
​
function draw() {
  
  background(220);
  noStroke();
​
  if (t<100){
  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 functions
​
w[m][i][j][k]=w[m][i][j][k]+dt1*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*dt1*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.6){
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,200*u[m][i][j][N2]);
  square(4*i,4*j,4);
  }
  }
  }
  }
  
  fill(0);
  circle(200,200,10);
  
  for (var i=0;i<N;i++){
  for (var m=1;m<M;m++){
  fill(0);
  ellipse(4*i,300-100*w[m][i][N2-2][N2],4);
  fill(255,255,0,255);
  ellipse(4*i,300-40*u[m][i][N2][N2],4);  
  fill(0,255,255,255);
  ellipse(4*i,300-10*P[m][i][N2][N2],4);
  }
//  fill(255,0,255,255);
//  ellipse(4*i,184,4);
  fill(0,0,255,255);
  ellipse(4*i,300-10*K[i][N2][N2],4);
  }
  
  fill(0);
  text(t,100,378);
  text(E,130,395);
  text("time step =",10,378);
  text("Total Energy =",10,395);
  text("(ref -161.7 Hartree)",260,395);
  text("Na 2+4+4+1 100 3d",5,20);
  text("u = sum_i u_i, u_i non-overlap elect density i=1,..,8",5,40);
  text("w_i characteristic function for support of u_i",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_i(x) = sum_jnoti integral dy*u_j(y)^2/(2*|x-y|) acting on u_i",5,180);
 text("Compute E by time stepping du_i/dt = -H_iu_i + charge normalisation",5,160);
 text("Compute P_i by solving -Laplacian P_i= sum_jnoti 2*PI*(uj^2) et cet",5,220);
text("yellow/red: u, blue:kernel pot, lightblue: elec pot, black: w, cut magenta,",10,360);
text("Update w_i by front track: dw_i/dt+jump*abs(grad w_i)",5,240);
}
​
​
Cookies

Sketch Files

Console

Preview
