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 D=15;

let D1=5;

let C=1;

let dt;

let dt1;

let h;

let t=1;

let d=0.05;

let M=7;

let M2=8;

let R=1.05;

​

function setup() {

  createCanvas(400, 400);

  h=20/N;

  dt=d*pow(h,2);

  dt1=dt/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)*h,2)+pow((j-N2)*h,2)+pow((k-N2+D)*h,2)+0.25*h*h;

   r=sqrt(r);

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

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

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

   }

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

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

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

   }

  if (r>R & i>N2 & j>N2){ 

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

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

   }

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

   P[4][i][j][k]=1.5/r;

   P[5][i][j][k]=1.5/r;

   P[6][i][j][k]=1.5/r;

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

   r=sqrt(r);

   if (i<N2-D/2){ 

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

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

   }

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

   P[1][i][j][k]=0.5/r;

   P[2][i][j][k]=0.5/r; 

   P[3][i][j][k]=0.5/r;

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

   r=sqrt(r);

  if (i>N2+D/2){ 

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

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

   }

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

   P[1][i][j][k]=0.5/r;

   P[2][i][j][k]=0.5/r; 

   P[3][i][j][k]=0.5/r;

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

   r=sqrt(r);

  if (j<N2-D/2){ 

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

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

   }

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

   P[1][i][j][k]=0.5/r;

   P[2][i][j][k]=0.5/r; 

   P[3][i][j][k]=0.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++){

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

if (r>R){

​

//Front tracking of electron characteristic functions    

w[m][i][j][k]=w[m][i][j][k]+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)+ 5*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];

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

​

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

​

​

 }

 }

 }

}

​

​

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

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

  }

  }

  }

//  }

  fill(255,255,255,255);

  fill(0);

  circle(200-D*4,200,10);

  circle(200+D*4,200,10);

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

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

  fill(0);

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

  fill(255,255,0,255);

  ellipse(4*i,300-400*u[m][i][N2][N2],4);  

  fill(0,255,255,255);

  ellipse(4*i,300-100*P[m][i][N2][N2],3);

  }

  fill(0,0,255,255);  

  ellipse(4*i,300-30*K[i][N2][N2],4);

  }

  fill(0);

  text(t,100,378);

  text(E+9/(sqrt(2)*D*h)+1/(2*D*h)+sqrt(2)/(D*h)+5.61,140,395);

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

  text("Dissociation Energy =",10,395);

  text("(ref -0.446 Hartree)",270,395);

  text("NH3 N kernel +3 R=1.05",5,20);

  text("u = sum_i u_i, w_i char function for supp u_i non-overlap",5,40);

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

  text("over u_i with continuity of u across free bdy updated by front track w_i",5,80); 

  text("H_i = - 0.5*Laplacian - K(x) + 2*P_i(x)",5,120);

  text("K(x) kernel potential, P_i(x) electron potentials",5,100);

  text("P_1(x) = integral dy*u_2(y)^2/(2*|x-y|) acting on u_1 et cet",5,160);

 text("Compute E by time stepping du_i/dt = -H_iu_i + normalisation",5,140);

 text("Compute P_1 by solving -Laplacian P_1 =2*PI*(u_2^2) et cet",5,180);

text("green: u_1,u_2, blue: potentials, black/red: charac func, mid cross cut,",10,360);

text("Update charac func w_i by front track: dw_i/dt+jump*abs(grad w_i)",5,320);

}

​

​