createCanvas(400, 400);

}

​

function draw() {

  background(220);

}

let N=800;

let Q=5;

let Q0=0;

let D=10;

let A=2;

let u = [];

let P = [];

let K = [];

let K1=[];

let M =[];

let Z;

let E=[];

let normu=[];

let energy=[];

let energytot;

let slider=[];

let dt;

let t=0;

let h;

let start;

let stop;

let fix=0;

let S=0;

let diff=0.1;

​

function setup() {

  createCanvas(400, 400);

  h=D/N;

  dt=0.5*pow(h,2);

  M[0]=0;

  M[1]=5;

  M[2]=20;

  M[3]=100;

  M[4]=220;

  M[5]=800;

  M[6]=800;

  slider[0]=createSlider(2,79,37,1);

  slider[1]=createSlider(2,2,2,1);

  slider[2]=createSlider(0,8,8,1);

  slider[3]=createSlider(0,18,18,1);

  slider[4]=createSlider(0,32,8,1);

  slider[5]=createSlider(0,18,1,1);

  slider[6]=createSlider(0,8,1,1);

  start=createButton('Start');

  stop=createButton('Stop');

  start.position(350,365);

  start.mousePressed(startsimulation);

  stop.mousePressed(stopsimulation);

  stop.position(350,385);

  Z=slider[0].value();

  for (var q=1;q<7;q++){

  E[q]=slider[q].value();

  normu[q]=0;

  energy[q]=0;

  }

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

  u[i]=100*exp(-2*(i*h));  

  P[i]=0;

  }

​

  u[M[Q]]=0;

​

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

  if (i>A){

  K[i]=Z/(i*h);

  }

  if (i<A+1){

  K[i]=Z/(A*h); 

  }

  }

  for (var q=Q0;q<Q;q++){

  for (var i=M[q];i<M[q+1];i++){

  P[i]=0;

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

P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h*min(1/(i*h),1/(j*h));

  if (j>M[q]-1 && j<M[q+1]){

P[i]=P[i]-0.5*(1/E[q+1])*pow(u[j],2)*pow(j*h,2)*h*min(1/(i*h),1/(j*h));

  }

  }

  }

  }

}

​

​

function draw() {

​

//if (t>2000){

//Q0=2;

//}

  background(255);

if (S==1){

  t=t+1;

//  M[1]=13;

  for (var q=Q0+1;q<Q;q++){

  if (u[M[q]]>u[M[q]-1]+diff){  

  M[q]=M[q]-2;

  }

  if (u[M[q]]<u[M[q]-1]-diff){

  M[q]=M[q]+1;

  }

  }

  u[0]=u[1]/(1-Z*h); 

  u[M[Q]]=0;

  for (var q=1;q<Q+1;q++){

  energy[q]=0;

  }

  for (var q=Q0+1;q<Q;q++){

  u[M[q]]=u[M[q]+1];

//  if (q>Q0+2){ 

  u[M[q]-1]=u[M[q]-2];

//  }

  }

  for (var q=Q0+1;q<Q+1;q++){

  for (var i=M[q-1]+1;i<M[q]-1;i++){

  for (var k=0;k<2;k++){

 u[i]=u[i]+0.5*dt*(u[i+1]-2*u[i]+u[i-1])/pow(h,2)+0.5*dt*(u[i+1]-u[i-1])/(h*i*h) + dt*K[i]*u[i] - dt*2*P[i]*u[i];

  }

  }

  normu[q]=0;

  for (var i=M[q-1];i<M[q];i++){ 

  normu[q]= normu[q] + pow(u[i]*(i*h),2)*h;

  }

  for (var i=M[q-1];i<M[q];i++){ 

  u[i]=sqrt(E[q])*u[i]/sqrt(normu[q]);

  }

  energy[q]=0;

  for (var i=M[q-1]+1;i<M[q];i++){   

  energy[q]= energy[q] + 0.5*pow((u[i+1]-u[i])/h,2)*pow(i*h,2)*h - pow(u[i],2)*K[i]*pow((i*h),2)*h + P[i]*pow(u[i],2)*pow(i*h,2)*h;

  }

  }

  energytot=0;

  for (var q=Q0+1;q<Q+1;q++){

  energytot = energytot +energy[q];

  }

​

  for (var q=Q0;q<Q+1;q++){

  for (var i=M[q]+1;i<M[q+1]+1;i++){

  P[i]=0;

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

  P[i]=P[i]+0.5*pow(u[j],2)*pow(j*h,2)*h*min(1/(i*h),1/(j*h));

  if (j>M[q] && j<M[q+1]){

  P[i]=P[i]-0.25*pow(u[j],2)*pow(j*h,2)*h*min(1/(i*h),1/(j*h));

  }

  }

  }

  }

 noStroke();

  for (var q=1;q<Q+1;q++){

  text(floor(1000*M[q]*h)/1000,150,280+q*10);

  text(floor(100*normu[q])/100,250,280+q*10);

  text(floor(100*energy[q])/100,20,280+q*10);

  text(M[q],200,280+q*10);

  }

  text(floor(10*energytot)/10,20,300+70);

  text("Shell Energies",20,250);

  text("Total Energy",20,350);

  text("Shell radii",150,250);

  text("Shell charges",250,250); 

  text("red: electron wave function, blue: kernel pot, green: electron pot",10,220);

//  text(M[Q-1]*h,200,200);

  text(Z,20,395);

  text(E[1],80,395);

  text(E[2],130,395);

  text(E[3],180,395);

  text(E[4],230,395);

  text(E[5],280,395);

  text(E[6],330,395);

//  ellipse(M[Q-1],200,10);

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

  fill(255,0,0,255);

  ellipse(0.5*i,200-10*u[i],3);

  fill(0,0,255,255);

  ellipse(0.5*i,200-2*K[i],3);

  fill(0,255,0,255);

  ellipse(0.5*i,200-2*P[i],3);

  fill(0,255,255,255);

//  ellipse(i,200-10*(K[i]-2*P[i]),4);

  fill(255,255,0,255);

//  ellipse(i,200-200*K1[i],4);

  fill(0);

  ellipse(i,200,2);

  }

//  text(floor(10000*(energy1+energy2+energy3))/10000,10,340);

stroke(0);

line(5,0,5,200);

}

text(t,10,190);

if (S==0){

  for (var q=1;q<7;q++){

  E[q]=slider[q].value();

  normu[q]=0;

  energy[q]=0;

  }

  Z=slider[0].value()

  fill(0,0,255,100);

  circle(200,200,300);

  fill(0,255,0,255);

  circle(200,200,150);

  fill(255,0,0,255);

  circle(200,200,50);

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

  ellipse(i,200-100*exp(-abs(200-i)*0.02),4);

  }

  fill(0);

  circle(200,200,5);

  noStroke();

  fill(0);

  text("2",205,200);

  text("E[2]",240,200);

  text("E[3]",280,200);

  text("E[4]",360,200);

  text("RealQM Atom Simulator: Ground States + 1st,2nd Ions",50,20);

  text("Kernel Charge 1 < Z < 86.",50,40);

  text("1st Shell: E[1] = 2, 2nd Shell: E[2], 3rd Shell: E[3].... Electrons.",50,100);

  text("Choose Z, E[2], E[3],..with sliders below and press Start!",50,120);

  text("RealQM = System of nonoverlapping electron densities",50,60);

  text("meeting at Bernoulli free boundary.",50,80);

  text("Electronic shell struture in accordance with observation:",50,230);

  text("Lithium Z=3: 2+1 (7.48),1st 2 (7.28)",50,250);

  text("Beryllium Z=4: 2+2 (14.6), 1st 2+1 (14.2), 2nd 2 (13.5)",50,265);

  text("Boron Z=5: 2+2+1 (24.5),1st 2+2 (24.2)",50,280);

  text("Carbon Z=6: 2+2+2 (37.5),1st 2+2+1 (37.1), 2nd 2+2 (36.5)",50,295);

  text("Nitrogen Z=7: 2+3+2 (54.6), 1st 2+2+2 (54.1)",50,310);

  text("Oxygen Z=8: 2+3+3 (74.8), 1st 2+3+2 (74.2)",50,325);

  text("Fluorine Z=9: 2+4+3 (99.8), 1st 2+4+2 (99.2)",50,340);

  text("Neon Z=10: 2+4+4 (128.5), 1st 2+4+3 (127.7), 2nd 2+4+2 (126.2)",50,355);

//  text("Test other electron distributions, e g Oxygen 2+6 or 2+4+2:",50,370);

//  text("Argue that (for Oxygen) 4 or 6 do not fit into 2nd shell.",50,385);

  text("Z =",0,395);

  text(Z,20,395);

  text(E[1],80,395);

  text(E[2],130,395);

  text(E[3],180,395);

  text(E[4],230,395);

  text(E[5],280,395);

  text(E[6],330,395);

  text("E[1] =",45,395);

  text("E[2] =",95,395);

  text("E[3] =",145,395);

  text("E[4] =",195,395);

  text("E[5] =",245,395);

  text("E[6] =",295,395);

}

if (S==2){

noStroke();

  fill(0);

  text("3(sub)Shell Atoms: Ground State+Ions SpherSym.",50,10);

  text("Kernel K(x)= Z/|x|, 2<Z<11, u(r) wave function,",50,30);

  text("u(r)^2 density of electrons in shell1, shell2,3 subshells, r=|x|.",20,50);

  text("Normalisation: 0<E2<9, 0<=E3<5.",50,70);

   text("integral_1 u(r)^2*r^2*dr = 2, integral_2,3 u(r)^2*r^2*dr =E2,3.",50,90);

  text("Minimise Total Energy:",50,110);

   text("integral (0.5*(du/dr)^2 - K*u^2 + P*u^2)*r^2*dr,",50,130);

  text("P(r)=0.5*E1red*integral_1 u(s)^2*min(1/r,1/s)*s^2*ds ",50,150);

  text("+ 0.5*integral_2 u(s)^2*(1/r)*s^2*ds for r in shell 1.",50,170);

  text("P(r) similar in shell 2,3.",50,190);

  text("Time step du/dt = 0.5*(d2u/dr2+2*r*du/dr)-K(r)*u+2*P*u in shell1,2,3", 10,220);

  text("+ homogeneous Neumann for shell1,2,3 at free boundary + normalisation.", 10,240);

  text("Free boundary between shell 1,2,3 established by continuity of u(r).",10,260);

  text("red: computed solution u1(r)",40,280);

  text("energy shell1 = ",40,320);

  text("energy shell2 = ",40,335);

  text("energy shell3 = ",40,350);

  text("total energy = ",40,365);

  text("r-axis",350,210);

  text("u-axis",10,10);

  text("free boundary1 =",40,380);

  text("picometer",230,380);

  text("au",180,380);

  text("picometer",230,395);

  text("au",180,395)

}

}

function startsimulation(){

  S=1;

}

function stopsimulation(){

  S=0;

}

​