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

let C=1;

let dt;

let h;

let t=1;

let d=0.2;

let R=1;

let M=3;

​

function setup() {

  createCanvas(400, 400);

  slider1 = createSlider(0,20,10);

//  slider2 = createSlider(0,2,1);

  h=10/N;

  dt=d*pow(h,2);

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

  w[m]=[];

  u[m]=[];

  c[m]=1;

  charge[m]=1;

  P[m]=[];

  norm[m]=0;

  }

  for (var m=1;m<M;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)+0.4*h*h;

   r=sqrt(r);

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

   if(i<N2 &r<R){

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

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

   }

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

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

   r=sqrt(r);

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

   if(i>N2 &r<R){

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

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

   }

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

   }

   }

  }

​

​

}

​

function draw() {

  background(220);

  noStroke();

    D=slider1.value();

//  R=slider2.value();

​

  if (t<2000){

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

  c[1]=0.5*(u[1][i][j][k]-u[2][i][j][k]);

  c[2]=0.5*(u[2][i][j][k]-u[1][i][j][k]);  

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

r=sqrt(r);

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

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

r=sqrt(r);

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

​

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

​

//Front tracking of electron characteristic functions    

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

//Solve Poisson equation for electron potentials//

P[1][i][j][k]=P[1][i][j][k]+dt*(P[1][i+1][j][k]-6*P[1][i][j][k]+P[1][i-1][j][k]+P[1][i][j+1][k]+P[1][i][j-1][k]+P[1][i][j][k+1]+P[1][i][j][k-1])/pow(h,2)+2*PI*dt*pow(u[2][i][j][k],2);

P[2][i][j][k]=P[2][i][j][k]+dt*(P[2][i+1][j][k]-6*P[2][i][j][k]+P[2][i-1][j][k]+P[2][i][j+1][k]+P[2][i][j-1][k]+P[2][i][j][k+1]+P[2][i][j][k-1])/pow(h,2)+2*PI*dt*pow(u[1][i][j][k],2);

//Compute electron energy//

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

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

 }

 }

 }

 }

//Add kernel repulsion energy

E=E+1/(2*D*h);

​

​

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

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

  fill(255,0,0,1000*u[1][i][j][25]);

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

  fill(255,0,0,1000*u[2][i][j][25]);

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

  }

  }

//  }

  fill(255,255,255,255);

  fill(0);

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

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

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

  fill(0);

  ellipse(8*i,200-100*w[1][i][25][25],4);

  ellipse(8*i,200-100*w[2][i][25][25],4);

  fill(255,255,0,255);

  ellipse(8*i,300-100*u[1][i][25][25],6);  

  ellipse(8*i,300-100*u[2][i][25][25],4);

  fill(0,0,255,255);

  ellipse(8*i,300-100*P[1][i][25][25],3);

  ellipse(8*i,300-100*P[2][i][25][25],3);

  ellipse(8*i,300-30*K[i][25][25],6);

//  ellipse(8*i,300-30*(K[i][25][25]-2*P[1][i][25][25]),6);

//  fill(0,255,255,255);

//  ellipse(8*i,300-30*(K[i][25][25]-2*P[2][i][25][25]),6);

  }

​

​

  fill(0);

  text(t,100,378);

​

  text(E,140,395);

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

  text("Ground State Energy =",5,395);

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

  text("H+H Interactive: u = sum_i u_i, u_i non-overlap wavef elect i=1,2.",5,20);

  text("w_i charac func for supp of u_i updated by front track on fixed grid.",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("Vary distance between kernels by Slider to minimize total energy.",5,320);

  text("Kernel distance =  ",5,340);

  text(2*D*h,100,340);

text("red/yellow: u_1,u_2, blue:potentials, black: charac func, mid cut,",5,360);

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

}

​

​