Update on 2021.02.02:
I have made a couple of inline updates in the sections 5. and 7. below.
I have already begun cleaning up and reorganizing the code, and writing the boring document. The work so far has come along pretty well. I have by now achieved a satisfactory level of consistency in the numerical results for the hydrogen atom in a box.
As indicated in my last post here, I had found that it’s more useful to focus on the cellside of the mesh rather than on the physical size of the box or the number of nodes per side of the cube.
Now, given below are some of the details of certain further, systematic, trials which I conducted, in order to arrive at optimum ranges for numerical parameters.
Since the analytical solution is available only for the hydrogenic atoms (i.e. systems with a positively charged nucleus and just one electron, e.g. the hydrogen atom and the He+ ion), these systematic studies were conducted only for them.
If you came here expecting that I might have something to say about the reproducibility for the helium atom, then well, you will have to wait for some 2–3 weeks. The nature of the issues themselves is like that. You can’t hurry things like these too much, as the studies below sure bring out.
So, anyway, here are the highlights of the systematic studies I conducted and some of the representative results.
All the results reported in this post are in the atomic units.
1. Finitization policy to replace the singularity of the Coulomb field:
In my last post here, I had mentioned that in FDM, we have to use a finite value in place of the at the nucleus. As a necessary consequence, we have to adopt some policy for this finitization.
While conducting my studies, now I found that it is better to frame this policy not in terms of a chosen fraction of the cellside, but in terms of a certain relevant datum and a multiplying factor. The better procedure is this:
Whatever be the degree of mesh refinement, having constructed the mesh, calculate the alwaysfinite PE value at the FDM node right next to the singularity. Then, multiply this PE value by a certain factor, and use it in place of the theoretically value at the singular node. Let’s give a name to this multiplying factor; let’s call it the Coulomb Field’s Singularity Finitization Factor (“CFSF” for short).
Notice that using this terminology, a CFSF value of turns out to be exactly the same as using the half cellside rule. However, framing the finitization policy in terms of the CFSF factor has the advantage that it makes it easier to compare the differences in the relative sharpness of the potential well at different meshrefinement levels.
OK.
I then found that while using FDM, the eigenvalue solver is very sensitive to even small variations to the value of CFSF.
If you use a CFSF value of , then it turns out that the PE well does not go down enough in the neighbourhood of the singularity, and therefore, the reported groundstate eigenvalues can easily go to about , , or even worse. Refining the mesh doesn’t help—within the domain and mesh sizes which are both practicable on my laptop and relevant for the H and He atom modelling. Note, the analytical solution is: , exactly.
Conclusion: Using even a slightly lower CFSF spoils the results.
OTOH, if you use a CFSF of to , then the groundstate energy can go lower than the exact value of . Now, this is a sureshot indication that your numerical modelling has gone wrong.
In general, with FDM, you would expect that with mesh refinement, the convergence in the energy values would be not just monotonic but also onesided, and also that the convergence would occur from “above” (because the energy values here are negative). In other words, if my understanding of the theory of numerical analysis is correct, then a properly done (meaningful) numerical simulation cannot produce energies below the analytical solution of .
So, clearly, using a CFSF value even slightly greater than is bad for the health of the numerical simulations.
In the earlier trials reported in the last post, I had simply guessed that the value of might be good enough for my initial trials. Now it turns out that my computational modeller’s intuition was pretty much on target—or at least, that I was plain lucky! The CFSF value of indeed happens to be quite the best value to choose, given the rest of the parameters like the cellside, the domain size, and the rest of the details of this problem (viz., the strength of the Coulomb singularity, the nature of the Schrodinger equation, the use of uniform and structured meshes, the FDM discretization, etc.).
2. Simpleminded mesh refinement doesn’t produce consistent results:
Suppose you keep the domain size fixed at, say, , and vary the mesh refinement levels.
Now, your naive expectation might be that as you refine the mesh by increasing the number of nodes per side of the cube, you should get more and more accurate results. That’s our usual experience with problems like diffusion in continua, and even for problems like convection in fluids.
However, the category of the QM problems is different! Here, we have an eigenvalue problem that must be solved with a singular potential field. The naive expectations built on simple problems like the PoissonLaplace equation or the diffusion equation, go for a toss. Harmonic analysis might still apply in some form (frankly I don’t even know if it does!), but the singularity sure plays tricks!
This is an ugly fact, and frankly, I had not foreseen it. But it’s there. I had to keep myself reminding of the different nature of the eigenvalue problem together with the singular fields.
As you refine the mesh too much, then the absolute value of the PE at a node right next to the point of singularity increases without bound! This fact mandates that the finite value we (more or less) arbitrarily chose to use in place of the actually infinite value for the singular point, has itself to increase further too.
But, for some reasons not known to me (but which by now do feel vaguely reasonable!) the eigenvalue solver begins to experience difficulties with such increases in the absolute value of the PE value at the singularity. Roughly, the trouble begins to happen as the minimum potential energy (at the singular node) goes below or so. In fact, I even found that a highly refined mesh might actually report a positive value for the groundstate energy—no bonding but, on the contrary, a repulsion of the electron!
3. Wavefunction fields are important in my new approach, but they don’t always converge to the analytical solution very well!:
With a reasonable level of mesh refinement, the groundstate energy does monotonically approach the exact figure of However, I’ve found that a convergence in energy is not necessarily accompanied also by a good convergence trend in the absolute values of the wavefunction!
In the Hatom, for the groundstate analytical solution, the absolute value of the wavefunction has its maximum right at the nucleus; the wavefunction field forms a cusp at the nucleus, in fact. The analytical value for max goes like: . (That’s because in the atomic units, the Bohr radius is chosen to be exactly equal to , and so, at , the groundstate wavefunction for the Hatom becomes .)
Now the trouble I found was this:
With mesh refinement, even as the energy is nicely converging to something like (against ), the max might still be lingering around a lower figure like . The max values converge more slowly, and their convergence shows opposite trend!
For relatively coarse meshes (i.e. high of the FDM mesh), the max value is actually way much higher than the analytical solution; it even becomes as bad as or . As you refine the mesh, they do begin to fall down and approach the analytical solution.
However, with further mesh refinement, the max values continue to fall down! They cross the analytical solution level of too, and still continue to fall further! And, this behaviour occurs even as energy result is still approaching the exact solution in a niceandexpected monotonic manner.
So, the trouble is: Using the right mesh size is actually a tradeoff! You have to sacrifice some convergence on the energy number, so as to have a good (reliable) value for the max measure.
The trouble doesn’t stop there; see the next section.
4. Energies for the excitedstates don’t always come out very well:
With appropriately high levels of meshrefinement, the groundstate energy might be showing good convergence trends. Even the max values might be good enough (like or so). But the energy and/or max for the first excited state still easily give trouble.
The energy for the first excited state for the hydrogen atom is, by analytical solution, , exactly.
The numerical values, when the simulation is working right, could be like , or even better, say , or thereabout. But that happens only when the mesh is of the intermediate refinement (the cellside is neither too small nor too large).
However, with a more refined mesh (smaller cellsides), the PE well can remain more or less rightly shaped for the groundstate energy, but it can still become too deep for the firstexcited state energy! The first excited state energy can suddenly get degraded to a value like .
Indeed, there seems to be some kind of a numerical compensation going on in between the max values and the energy values, especially for the firstexcited state energies. The groundstate energies remain much better, in relative terms. (If the mesh refinement is very high, even the groundstate energy goes off the track to something like or even positive values. That’s what I meant by “appropriately” high levels of mesh refinement.)
I didn’t compare the numerical results with the analytical solutions for energies or max values for secondexcited states or higher. Computation of the bonding energy makes reference only to the ground state, and so, I stopped my exploration of this side of the FDM + eigenvalue solver behaviour at this stage.
5. Atomic sizes reported by the numerical modeling show very good trends:
Another important consideration in my new approach has to do with the atomic radius of the atoms being modelled (hydrogen and helium).
After optimizing the mesh refinement (i.e., effectively, the cellside), I conducted a series of numerical trials using different domain sizes (from through ), and implemented a roughandready code to estimate the following measure:
The side of the nucleuscentered subcube in which roughly (or ) of the probability cloud is contained.
This size can be taken as a good measure for the atomic size.
In the above working definition, I say roughly , because I didn’t care to interpolate the wavefunction fields in between their nodal values. What this means is that the side of the subcube effectively changes only in the integer steps, and therefore, the percentage of the subcube contained may not be exactly ; it could be for one domain size, and for another domain size, just to pick up some numbers.
But even while using this rough and ready measure (and implementation), I found that the results were quite meaningfully consistent.
But why conduct these trials?
Well, realize that (1) the simulation box has a finite size, and (2) the Dirichlet conditions are being imposed at all the boundary nodes. Given these two constraints, the solution is going to show boundary/edgeeffects, i.e., the solution is going to depend on the domain size.
Now, in my approach, the spread of the probability cloud enters the calculations in a crucial manner. Numerically “extracting” the size of the simulated atom was, therefore, an important part of optimizing the simulations.
The expected behaviour of the above mentioned “size effect” was that as the domain size increases, the calculated atomic size, too, should increase. The question was: Were these differences in the numerically determined sizes important enough? did they vary too much? if yes, how much? The following is what I found:
First, I fixed the domain size (cube side) at , and varied the mesh refinement (from roughly nodes per side to and ). I found that the calculated atomic sizes for the hydrogen atom varied but in a relatively small range—which was a big and happy surprise to me. The calculated size went from while using a coarse mesh (requiring eigenvalue computation time of about seconds) to a value like for an intermediate refinement of the mesh (exe. time 2 min. 32 seconds i.e. 152 seconds), to for as fine a mesh as my machine can handle (, which required an exe. time of about 20 minutes i.e. 1200 seconds, for each eigenvalue computation call). Remember, all these results were for a domain size of .
Next, I changed the domain cube side to , and repeated the trials, for various levels of mesh refinements. Then, ditto, for the domain side of and .
Collecting the results together:

 coarse:
 intermediate:
 fine:

 coarse:
 intermediate:

 coarse:
 intermediate:
 fine:

 coarse:
 intermediate:
 fine:
 very fine:
You might be expecting very clearcut trends and it’s not the case here. However, remember, due to the trickiness of the eigenvalue solver in the presence of a “singular” PE well, not to mention the roughness of the sizeestimation procedure (only integersized subcubes considered, and no interpolations of to internodal values), a monotonic sort of behaviour is simply not to be expected here.
Indeed, if you ask me, these are pretty good trends, even if they are only for the hydrogen atom.
Note, for the helium atom, my new approach would require giving eigenvalue computation calls thousands of times. So, at least on this count of atomic radius computation, the fact that even the coarse or midlevel mesh refinement results didn’t vary much (they were in the range of to ) was very good. Meaning, I don’t have to sacrifice a lot of accuracy due to this one factor taken by itself.
For comparison, the atomic size (diameter) for the hydrogen atom is given in the literature (Wiki), when translated into atomic units, comes out variously as: (1) using some “empirical” curvefitting to some indirect properties of gases; (2) while using the van der Waal criterion, and (3) using “calculations” (whose basis or criteria I do not know in detail).
Realize, the van der Waal measure is closest to the criterion used by me above. Also, it is only expected that when using FDM, due to the numerical approximations, just the way the FDM groundstate energy values should come out algebraically greater (they do, say vs. the exact datum of ), the FDM max measure should come out smaller (it does, say vs. the analytical solution of ), similarly, for the same reasons, the roughandready estimated atomic size should come out as greater (it does, say to as the domain size increases from to , the. van der Waal value being ).
Inline update on 2021.02.02 19:26 IST: After the publication of this post, I compared the abovementioned results with the analytical solution. I now find that the sizes of the subcubes found using FDM, and using the analytical solution for the hydrogen atom, come out as identical! This is a very happy news. In other words, making comparisons with the van der Waal size and the other measure was not so relevant anyway; I should have compared the atomic sizes (found using the subcubes method) with the best datum, which is, the analytical solution! To put this finding in some perspective, realize that the FDMcomputed wavefunctions still do differ a good deal from the analytical solution, but the volume integral for an easy measure like does turn out be the same. The following proviso’s apply for this finding: The good match between the analytical solution and the FDM solution are valid only for (i) the range of the domain sizes considered here (roughly, to ), not for the smaller box sizes (though the two solution would match even better for bigger boxes), and (ii) only when using the Simpson procedure for numerically evaluating the volume integrals. I might as well also note that the Simpson procedure is, relatively, pretty crude. As the sizes of the subcubes go on increasing, the Simpson procedure can give volume integrals in excess of for both the FDM and the analytical solutions. Inline update over.
These results are important because now I can safely use even a small sized domain like a side cube, which implies that I can use a relatively crude mesh of just nodes per side too—which means a sufficiently small runtime for each eigenvalue function call. Even then, I would still remain within a fairly good range on all the important parameters.
Of course, it is already known with certainty that the accuracy for the bonding energy for the helium atom is thereby going to get affected adversely. The accuracy will suffer, but the numerical results would be on the basis of a sweetzone of all the numerical parameters of relevance—when validated against hydrogen atom. So, the numerical results, even for the helium atom, should have greater reliability.
Considerations like conformance to expected behaviour in convergence, stability, and reliability are far more important considerations in numerical work of this nature. As to sheer accuracy itself, see the next section too.
6. Putting the above results in perspective:
All in all, for the convergence behaviour for this problem (eigenvalueeigenvector with singular potentials) there are no easy answers. Not even for just the hydrogen atom. There are tradeoffs to be made.
However, for computation of bonding energy using my new approach, it’s OK even if a good tradeoff could be reached only for the groundstate.
On this count, my recent numerical experimentation seems to suggest that using a mesh cellside of or should give the most consistent results across a range of physical domain sizes (from through ). The atomic size extracted from the simulations also show good behaviour.
Yes, all these results are only for the hydrogen atom. But it was important that I understand the solver behaviour well enough. It’s this understanding which will come in handy while optimizing for the helium atom—which will be my next step on the simulation side.
The trends for the hydrogen atom would be used in judging the results for the the bonding energy for the helium atom.
7. The discussed “optimization” of the numerical parameters is strictly for my laptop:
Notice, if I were employed in a Western university or even at an IIT (or in an Indian government/private research lab), I would have easy access to supercomputers. In that case, much of this study wouldn’t be so relevant.
The studies regarding the atomic size determination, in particular, would still be necessary, but the results are quite stable there. And it is these results which tell me that, had I have access to powerful computational resources, I could have used larger boxes (which would minimize the edge/size effect due to the finite size of the box), and I could have used much, much bigger meshes, while still maintaining the allimportant mesh cellside parameter near the sweet spot of about to . So, yes, optimization would still be required. But I would be doing it at a different level, and much faster. And, with much better accuracy levels to report for the helium atom calculations.
Inline update on 2021.02.02 19:36 IST: Addendum: I didn’t write this part very well, and a misleading statement crept in. The point is this: If my computational resources allow me to use very big meshes, and then I would also explore cellsides that are smaller than the sweetspot of to . I’ve been having a hunch that the eigenvalue solver would still not show up the kind of degeneracy due to very deep PE well, provided that the physical domain size also were to be made much bigger. In short, if very big meshes are permissible, then there is a possibility that another sweetspot at smaller cellsizes could be reached too. There is nothing physical about the to range alone, that’s the point. Inline update over.
The specifics of the study mentioned in this post was largely chosen keeping in the mind the constraint of working within the limits of my laptop.
Whatever accuracy levels I do eventually end up getting for the helium atom using my laptop, I’ll be using it not just for my planned document but also for my very first arXiv/journal paper. The reader of the paper would, then, have to make a mental note that my machine could only support a mesh size of only nodes at its highest end. For FDM computations, that still is a very crude mesh.
And, indeed, for the reasons given above, I would in fact be reporting the helium atom results for meshes in between to nodes per side of the cube, not even nodes. All the rest of the choices of the parameters were made keeping in view this limitation.
8. “When do you plan to ship the code?”
I should be uploading the code eventually. It may not be possible to upload the “clientside” scripts for all the trials reported here (simply because once you upload some code, the responsibility to maintain it comes too!). However, exactly the same “server” or “backend” side code will sure be distributed, in its entirety. I will also be giving some indication of the kind of codesnippets I used in order to implement the above mentioned studies. So, all in all, it should be possible for you to conduct the same/similar trials and verify the above given trends.
I plan to clean up and revise the code for the hydrogen atom a bit further, finalize it, and upload it to my GitHub account within, say, a week’s time. The cleaned up and revised version of the heliumatom code will take much longer, may be 3–4 weeks. But notice, the heliumatom code would be giving calls to exactly the same library as that for the hydrogen atom.
All in all, you should have a fairly good amount of time to go through the code for the boxes (something which I have never uploaded so far), run it, run the above kind of studies on the solid grounds of the hydrogen atom, and perhaps even spot bugs or suggest better alternatives to me. The code for the helium atom would arrive by the time you run through this gamut of activities.
So, hold on just a while, may be just a week or even less, for the first code to be put on the GitHub.
On another note, I’ve almost completed compiling a document on the various set of statements for the postulates of QM. I should be uploading it soon too.
OK, so look for an announcement here and on my Twitter thread, regarding the shipping of the basic code library and the userscript for the hydrogen atom, say, within a week’s time. (And remember, this all comes to you without any charge to you! (For that matter, I am not even in any dayjob.))
A song I like:
(Hindi) दिल कहे रुक जा रे रुक जा (“dil kahe ruk jaa re ruk jaa”)
Lyrics: Sahir Ludhiyanvi
Music: LaxmikantPyarelal
Singer: Mohammed Rafi
[Another favourite right from my highschool days… A good quality audio is here [^]. Many would like the video too. A good quality video is here [^], but the aspectratio has gone awry, as usual! ]
History:
— 2020.01.30 17:09 IST: First published.
— 2021.02.02 20:04 IST: Inline updates to sections 5. and 7 completed. Also corrected a couple of typos and streamlined just a few sentences. Now leaving this post in whatever shape it is in.