let f = [];

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 D1=10;

let D2=10;

let C=1;

let dt;

let h;

let t=1;

let d=0.02;

let M=6;

let M2=7;

let R=1.2;

let Z1=8;

let Z2=11;

​

function setup() {

  createCanvas(400, 400);

  h=20/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]=[];

    P[m][i]=[];

  for (var j=0;j<N+1;j++){

    w[m][i][j]=[];

    u[m][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;

    P[m][i][j][k]=0;

    }

    }

    }

  }

  for (var i=0;i<N+1;i++){

    K[i]=[];

  for (var j=0;j<N+1;j++){

    K[i][j]=[];

  for (var k=0;k<N+1;k++){

    K[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+D1)*h,2)+pow((j-N2)*h,2)+pow((k-N2)*h,2)+0.5*h*h;

    r=sqrt(r);

   K[i][j][k]=Z1/r;

   if (r>R & i<N2-D1){

   u[1][i][j][k]=exp(-r);

   w[1][i][j][k]=1;

   }

   if (r>R & i>N2-D1 & i<N2){

   u[2][i][j][k]= exp(-r);

   w[2][i][j][k]=1;

   }

   if (r<R){

   u[3][i][j][k]=exp(-4*r);

   w[3][i][j][k]=1;

   }

   charge[3]=Z1-2;

   P[3][i][j][k]=1/R;

   P[1][i][j][k]=P[1][i][j][k]+(Z1-1)/(2*r);

   P[2][i][j][k]=P[2][i][j][k]+(Z1-1)/(2*r);

   P[3][i][j][k]=P[3][i][j][k]+(Z1-3)/(2*r);

   P[4][i][j][k]=P[4][i][j][k]+Z1/(2*r);

   P[5][i][j][k]=P[5][i][j][k]+Z1/(2*r);

   P[6][i][j][k]=P[6][i][j][k]+Z1/(2*r);

  r=pow((i-N2-D2)*h,2)+pow((j-N2)*h,2)+pow((k-N2)*h,2)+0.5*h*h;

   r=sqrt(r);

   K[i][j][k]=K[i][j][k]+Z2/r;

   if (r>R & i>N2){

   u[4][i][j][k]=exp(-r);

   w[4][i][j][k]=1;

   }

   if (r<R & i>N2){

   u[5][i][j][k]=exp(-r);

   w[5][i][j][k]=1;

   }

   charge[5]=Z2-1;

   P[1][i][j][k]=P[1][i][j][k]+ Z2/(2*r);

   P[2][i][j][k]=P[2][i][j][k]+ Z2/(2*r);

   P[3][i][j][k]=P[3][i][j][k]+ Z2/(2*r);

   P[5][i][j][k]=P[5][i][j][k]+0.5/R;

   P[4][i][j][k]=P[4][i][j][k]+(Z2-3)/(2*r);

//   P[5][i][j][k]=P[5][i][j][k]+(Z2-3)/(2*r);

  }

  }

  }

}

​

​

function draw() {

  background(220);

  noStroke();

​

  if (t<150){

  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++){

​

//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];

}

//Solve Poisson equation for electron potentials//

for (var m=1;m<M;m++){    

  c[m]=0;

  for (var n=1;n<M;n++){

  if (n!=m){

  c[m]=c[m]+pow(u[n][i][j][k],2); 

  }

  }

 } 

 c[3]=c[3]+pow(u[3][i][j][k],2)*(Z1-3)/(Z1-2);

 c[5]=c[5]+pow(u[5][i][j][k],2)*(Z2-3)/(Z2-2);

for (var m=1;m<M;m++){  

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];

}

​

//Compute energy//

for (var m=1;m<M;m++){  

if (m<7){

if (w[m][i][j][k]>0.9){

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);

}

}

}

}

}

​

//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]=sqrt(charge[m])*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,1000*u[m][i][j][N2]);

  square(4*i,4*j,4);

  }

  }

  }

  }

  fill(255,255,255,255);

  fill(0);

  circle(200,200,10);

  for (var m=1;m<M;m++){

  for (var i=0;i<N;i++){

  fill(255,255,0,255);

  ellipse(4*i,200-100*u[m][i][N2][N2],4);  

  fill(0,255,255,255);

  ellipse(4*i,200-10*P[m][i][N2][N2],4);

  fill(0,0,255,255);

  ellipse(4*i,200-10*K[i][N2][N2],4);

  }

}

  fill(0);

  text(t,100,378);

  text(E+Z1*Z2/((D1+D2)*h)+290.42,150,395);

  text("Ground state energy",5,395);

  text("Valence+2+1 kernels Z1, Z2",5,20);

  text("time step =",10,378);

  text("LiO: Z=3,8, 0.5*h*h, R=1.2, (ref 0.13)",5,260);

  text("NaO: Z=11,8, 0.4*h*h, R=1, -74.8, (ref 0.10)",5,280);

  text("Mg: Z=12, 0.5*h*h, R=1.2, -200.3, 0.005 (ref 0.003)",5,300);

  text("S: Z=16, 0.5*h*h, R=1.2, -399, 0.219 (ref 0.17)",5,320);

}

​

​

​