let f = [];

let K = [];

let P = [];

let w = [];

let charge = [];

let r = 0;

let r1=0;

let r4=0;

let E =0;

let c=[];

let N=100;

let N2=50;

let N3=50;

let N4=45;

let D=15;

let C=1;

let dt;

let dt1;

let h;

let t=1;

let d=0.2;

let M=6;

let M2=6;

let R1=2;

let R2=0.9;

let R3=2;

//  R2 = 0.5 D = 18 E = -2.83

//  R2 = 0.5 D = 35 E = =-2.48

// Dissociation energy = 0.35 Hartree

//  R2 = 1  D = 18  E = -1.927

//  R2 = 1  D = 35  E = - 1.836

// Dissociation energy = 0.091

//  R2 = 0.75 D = 18 E = -2.105

//  R2 = 0.75 D = 35 E = - 2.055

//  R2 = 0.6 D = 35 E =

​

function setup() {

  createCanvas(400, 400);

  h=10/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-N3)*h,2)+pow((j-N2)*h,2)+pow((k-N2)*h,2)+0.16*h*h;

   r=sqrt(r);

   if (r>R2 & r<R3 & i<N4){ 

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

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

   }

   if (r>R2 & r<R3 & i>N4 & j<N2 ){ 

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

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

   }

   if (r>R2 & r<R3 & i>N4 & j>N2){ 

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

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

   }

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

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

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

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

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

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

  r1=pow((i-N3-D)*h,2)+pow((j-N2+D)*h,2)+pow((k-N2)*h,2)+0.16*h*h;

   r1=sqrt(r1);

   if (r1<R1 & r>R3) { 

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

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

   }

//  charge[1]=0.01;

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

   P[1][i][j][k]=P[1][i][j][k]+0.5/r1

   P[2][i][j][k]=P[2][i][j][k]+0.5/r1;

   P[3][i][j][k]=P[3][i][j][k]+0.5/r1;

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

   r4=pow((i-N3-D)*h,2)+pow((j-N2-D)*h,2)+pow((k-N2)*h,2)+0.16*h*h;

   r4=sqrt(r4);

   if (r4<R1 & r>R3){ 

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

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

   }

//   charge[4]=0.01; 

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

   P[1][i][j][k]=P[1][i][j][k]+0.5/r4;

   P[2][i][j][k]=P[2][i][j][k]+0.5/r4;

   P[3][i][j][k]=P[3][i][j][k]+0.5/r4;

   P[4][i][j][k]=P[4][i][j][k]+0.5/r4;

  }

  }

  }

​

}

​

function draw() {

  background(220);

  noStroke();

​

  if (t<99){

  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 

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

​

//for (var q=1;q<6;q++){

if (r>R2){

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.44){

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 i=1;i<N;i++){

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

  fill(0,0,255,500*u[1][i][j][N2]);

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

  fill(255,255,0,500*u[2][i][j][N2]);

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

  fill(255,0,0,500*u[3][i][j][N2]);

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

  fill(0,255,0,500*u[4][i][j][N2]);

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

  fill(0,255,0,500*u[5][i][j][N2]);

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

  }

  }

  }

  fill(255,255,255,255);

  fill(0);

  circle(N3*4,200,10);

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

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

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

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

  fill(0);

//  ellipse(4*i,400-4*i-32,1);

  ellipse(4*i,300-100*w[m][i][N2+8][N2],4);

  fill(255,255,0,255);

  ellipse(4*i,300-100*u[m][i][N2+5][N2],4);  

  fill(0,255,255,255);

  ellipse(4*i,300-30*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-1/(2*D*h)-4/(sqrt(2)*D*h)-3.328,140,395);

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

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

  text("(ref 0.353)",270,395);

  text("HOH 2+1 valence electrons 3d 100^3 mesh",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_i(x) = sum_jnoti integral dy*u_j(y)^2/(2*|x-y|) acting on u_i",5,160);

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

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

text("yellow u black characteristic function mid cross cut ,",10,340);

text("blue: potentials, black: charac func dotted line",10,360);

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

}

​

​