# Finding a cozy n comfy enough a spot…

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 $3D$ box.

As indicated in my last post here, I had found that it’s more useful to focus on the cell-side 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 $-\infty$ 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 cell-side, 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 always-finite 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 $-\infty$ 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 $2.0$ turns out to be exactly the same as using the half cell-side 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 mesh-refinement 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 $1.8$, then it turns out that the PE well does not go down enough in the neighbourhood of the singularity, and therefore, the reported ground-state eigenvalues can easily go to about $-0.45$, $-0.2$, 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: $-0.5$, exactly.

Conclusion: Using even a slightly lower CFSF spoils the results.

OTOH, if you use a CFSF of $2.2$ to $2.5$, then the ground-state energy can go lower than the exact value of $-0.5$. Now, this is a sure-shot 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 one-sided, 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 $-0.5$.

So, clearly, using a CFSF value even slightly greater than $2.0$ 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 $2.0$ 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 $2.0$ indeed happens to be quite the best value to choose, given the rest of the parameters like the cell-side, 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. Simple-minded mesh refinement doesn’t produce consistent results:

Suppose you keep the domain size fixed at, say, $20.0$, 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 Poisson-Laplace 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 $-20$ or so. In fact, I even found that a highly refined mesh might actually report a positive value for the ground-state 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 ground-state energy does monotonically approach the exact figure of $-0.5$ 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 H-atom, for the ground-state 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 $\psi(x)$-max goes like: $0.564189584\dots$. (That’s because in the atomic units, the Bohr radius $a_0$ is chosen to be exactly equal to $1$, and so, at $r = 0$, the ground-state wavefunction for the H-atom becomes $\psi(x_{\text{at nucleus}}) = 1/\sqrt{\pi}$.)

With mesh refinement, even as the energy is nicely converging to something like $-0.4938884$ (against $-0.5$), the $\psi$-max might still be lingering around a lower figure like $0.516189$. The $\psi$-max values converge more slowly, and their convergence shows opposite trend!

For relatively coarse meshes (i.e. high $\Delta x$ of the FDM mesh), the $\psi$-max value is actually way much higher than the analytical solution; it even becomes as bad as $3.276834$ or $1.393707$. As you refine the mesh, they do begin to fall down and approach the analytical solution.

However, with further mesh refinement, the $\psi$-max values continue to fall down! They cross the analytical solution level of $0.564189584$ too, and still continue to fall further! And, this behaviour occurs even as energy result is still approaching the exact solution in a nice-and-expected monotonic manner.

So, the trouble is: Using the right mesh size is actually a trade-off! You have to sacrifice some convergence on the energy number, so as to have a good (reliable) value for the $\psi$-max measure.

The trouble doesn’t stop there; see the next section.

4. Energies for the excited-states don’t always come out very well:

With appropriately high levels of mesh-refinement, the ground-state energy might be showing good convergence trends. Even the $\psi$-max values might be good enough (like $0.52$ or so). But the energy and/or $\psi$-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, $-0.125$, exactly.

The numerical values, when the simulation is working right, could be like $-0.11$, or even better, say $-0.123$, or thereabout. But that happens only when the mesh is of the intermediate refinement (the cell-side is neither too small nor too large).

However, with a more refined mesh (smaller cell-sides), the PE well can remain more or less rightly shaped for the ground-state energy, but it can still become too deep for the first-excited state energy! The first excited state energy can suddenly get degraded to a value like $-0.04471003$.

Indeed, there seems to be some kind of a numerical compensation going on in between the $\psi$-max values and the energy values, especially for the first-excited state energies. The ground-state energies remain much better, in relative terms. (If the mesh refinement is very high, even the ground-state energy goes off the track to something like $-0.2692952$ 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 $\psi$-max values for second-excited 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 cell-side), I conducted a series of numerical trials using different domain sizes (from $7.5$ through $40.0$ ), and implemented a rough-and-ready code to estimate the following measure:

The side of the nucleus-centered sub-cube in which roughly $95 \%$ (or $99 \%$) 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 $95 \%$, because I didn’t care to interpolate the wavefunction fields in between their nodal values. What this means is that the side of the sub-cube effectively changes only in the integer steps, and therefore, the percentage of the sub-cube contained may not be exactly $95 \%$; it could be $97.5 \%$ for one domain size, and $95.3 \%$ 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-/edge-effects, 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 $10.0$, and varied the mesh refinement (from roughly $41$ nodes per side to $121$ and $131$). 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 $5.60$ while using a coarse mesh (requiring eigenvalue computation time of about $10$ seconds) to a value like $5.25$ for an intermediate refinement of the mesh (exe. time 2 min. 32 seconds i.e. 152 seconds), to $5.23$ for as fine a mesh as my machine can handle ($131 \times 131 \times 131$, 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 $10.0$.

Next, I changed the domain cube side to $15.0$, and repeated the trials, for various levels of mesh refinements. Then, ditto, for the domain side of $20.0$ and $40.0$.

Collecting the results together:

• $10.0$
• coarse: $5.60$
• intermediate: $5.25$
• fine: $5.23$
• $15.0$
• coarse: $6.0$
• intermediate: $5.62$
• $20.0$
• coarse: $6.40$
• intermediate: $5.50$
• fine: $5.67$
• $40.0$
• coarse: $4.0$
• intermediate: $5.5$
• fine: $6.0$
• very fine: $6.15$

You might be expecting very clear-cut 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 size-estimation procedure (only integer-sized sub-cubes considered, and no interpolations of $\psi$ 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 mid-level mesh refinement results didn’t vary much (they were in the range of $5.25$ to $5.6$) 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) $0.94486306$ using some “empirical” curve-fitting to some indirect properties of gases; (2) $4.5353427$ while using the van der Waal criterion, and (3) $2.0031097$ 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 ground-state energy values should come out algebraically greater (they do, say $-0.49$ vs. the exact datum of $-0.5$), the FDM $\psi$-max measure should come out smaller (it does, say $0.52$ vs. the analytical solution of $\approx 0.56$), similarly, for the same reasons, the rough-and-ready estimated atomic size should come out as greater (it does, say $5.25$ to $5.67$ as the domain size increases from $10.0$ to $40.0$, the. van der Waal value being $4.54$ ).

Inline update on 2021.02.02 19:26 IST: After the publication of this post, I compared the above-mentioned results with the analytical solution. I now find that the sizes of the sub-cubes 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 sub-cubes method) with the best datum, which is, the analytical solution! To put this finding in some perspective, realize that the FDM-computed wavefunctions still do differ a good deal from the analytical solution, but the volume integral for an easy measure like $95 \%$ 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, $10$ to $40$), 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 sub-cubes go on increasing, the Simpson procedure can give volume integrals in excess of $1.0$ 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 $10.0$-side cube, which implies that I can use a relatively crude mesh of just $51$ nodes per side too—which means a sufficiently small run-time 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 sweet-zone 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 (eigenvalue-eigenvector with singular potentials) there are no easy answers. Not even for just the hydrogen atom. There are trade-offs to be made.

However, for computation of bonding energy using my new approach, it’s OK even if a good trade-off could be reached only for the ground-state.

On this count, my recent numerical experimentation seems to suggest that using a mesh cell-side of $0.2$ or $0.25$ should give the most consistent results across a range of physical domain sizes (from $7.5$ through $30.0$ ). 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 all-important mesh cell-side parameter near the sweet spot of about $0.20$ to $0.25$. 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 cell-sides that are smaller than the sweet-spot of $0.20$ to $0.25$. 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 sweet-spot at smaller cell-sizes could be reached too. There is nothing physical about the $0.20$ to $0.25$ 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 $131$ 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 $41$ to $81$ nodes per side of the cube, not even $131$ 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 “client-side” 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 code-snippets 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 helium-atom code will take much longer, may be 3–4 weeks. But notice, the helium-atom 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 $3D$ 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 user-script 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 day-job.))

A song I like:

(Hindi) दिल कहे रुक जा रे रुक जा (“dil kahe ruk jaa re ruk jaa”)
Lyrics: Sahir Ludhiyanvi
Music: Laxmikant-Pyarelal
Singer: Mohammed Rafi

[Another favourite right from my high-school days… A good quality audio is here [^]. Many would like the video too. A good quality video is here [^], but the aspect-ratio 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.

# Now I am become Bohmianism

1. About the title of this post:

Just before this Diwali, I had tweeted that I had made a resolution. The tweets went like this:

Let me note the text portions of these tweets (just in case I delete these some time later or so).

3:29 PM 13 Nov. 2020:

This year, Pune directly went from the monsoon air to the Diwali air. We seem to have tunnelled through the October heat!

3:55 PM, 13 Nov. 2020:

#Deepavali #Diwali #deepavali2020 #Diwali2020

[Diya lamp emoji, 3 times]

This is the *third* straight Diwali that I go jobless.

3:56 PM, 13 Nov. 2020:

My Diwali Resolution:

[Yes, there are going to be the usual New Year’s Resolutions according to the Western calender too!]

Alright.

We will come to the “tunnelling” part later. Also, the tweet related to my jobless-ness. [If the Indian IT industry has any sense of shame left at all, they would have prevented this circumstance. But more on this, too, later.]

For the time being, I want to focus on the last tweet, and say that, accordingly:

Now I am become Bohmianism.

As to the quaint grammar used in the expression, first consult this Wired article [^], also the Q&A at the Quora [^].

As to why I use “Bohmianism” instead of “a Bohmian”: Well, to know that, you have to understand Sanskrit. If you do, then refer to the Gita, Chapter 11, verse 32, the compound phrase “कालोऽस्मि” (“kaalo smi”). I just tried to keep a similar grammatical form. … But let me hasten to add that I am not a Sanskrit expert, and so, going wrong is always a possibility. However, I also think that here I have not.

Hence the title of this post.

Now, going over to the Bohmianism i.e. the Bohmian mechanics proper…

2. Material on the Bohmian mechanics (BM):

The following is a partial list of papers and other material on BM that I have downloaded. I am giving you the list in a roughly chronological order. However, my reading isn’t going to be in any particular order. I have not read them all yet. In fact, I’ve just got going with some them, as of now.

Also note, I expect that

• Some of this material might have become outdated by now
• I may run into some other related topics as my studies progress

Alright. On to the list…

2.1 Student theses:

Antony Valentini (1992) “On the pilot-wave theory of classical, quantum and subquantum physics,” Ph.D. Thesis, International School for Advanced Studies, Trieste

Caroline Colijn (2003) “The de Broglie-Bohm causal interpretation of quantum mechanics and its application to some simple systems,” Ph.D. Thesis, University of Waterloo.

Paulo Machado (2007) “Computational approach to Bohm’s quantum mechanics,” Ph.D. Thesis, McMaster University

Jeff Timko (2007) “Bohmian trajectories of the two-electron helium atom,” Master’s Thesis, University of Waterloo

Leopold Kellers (2017) “Making use of quantum trajectories for numerical purposes,” Master’s Thesis, Technische Universität München

2.2. Code:

Dane Odekirk (2012) “Python calculations of Bohmian trajectories,” GitHub, 12 December 2012. https://github.com/daneodekirk/bohm

2.3. Papers:

C. Philippidis, C. Dewdney and B. J. Hiley (1978) “Quantum interference and the quantum potential,” https://www.researchgate.net/publication/225228072

Berthold-Georg Englert, Marlan O. Scully, Georg Sussmann and Herbert Walther (1992) “Surrealistic Bohm trajectories,” Z. Naturforsch. 47 a, 1175–1186.

Robert E. Wyatt and Eric R. Bittner (2003) “Quantum mechanics with trajectories: quantum trajectories and adaptive grids,” arXiv:quant-phy/0302088v1 11 Feb 2003

Roderich Tumulka (2004) “Understanding Bohmian mechanics: A dialogue,” Am. J. Phys., vol. 72, no. 9, September 2004, pp. 1220–1226.

D.-A. Deckert, D. Dürr, P. Pickl (2007) “Quantum dynamics with Bohmian trajectories,” arXiv:quant-phy/0701190v2 13 May 2007

Guido Bacciiagaluppi and Antony Valentini (2009) “Quantum theory at the crossroads: Reconsidering the 1927 Solvay conference,” Cambridge UP, ISBN: 9780521814218 arXiv:quant-ph/0609184v2 24 Oct 2009 [Note: This is actually a book.]

M. D. Towler and N. J. Russell (2011) “Timescales for dynamical relaxation to the Born rule,” arXiv:1103.1589v2 [quant-ph] 27 Sep 2011

Michael Esfeld, Dustin Lazarovici, Mario Hubert, Detlef Dürr (2012) “The ontology of Bohmian mechanics,” preprint, British Journal for the Philosophy of Science

Travis Norsen (2013) “The pilot-wave perspective on quantum scattering and tunneling,” m. J. Phys., vol. 81, no. 4, April 2013, pp. 258–266. arXiv:1210.7265v2 [quant-ph] 9 Jan 2013

Travis Norsen (2013) “The pilot-wave perspective on spin,” arXiv:1305.1280v2 [quant-ph] 10 Sep 2013

Kurt Jung (2013) “Is the de Broglie-Bohm interpretation of quantum mechanics really plausible?,” Journal of Physics: Conference Series 442 (2013) 012060 doi:10.1088/1742-6596/442/1/012060

Samuel Colin and Antony Valentini (2014) “Instability of quantum equilibrium in Bohm’s dynamics,” Proc. R. Soc. A 470: 20140288. http://dx.doi.org/10.1098/rspa.2014.0288

W. B. Hodge, S. V. Migirditch and W. C. Kerr (2014) “Electron spin and probability current density in quantum mechanics,” Am. J. Phys., vol. 82, no. 7, July 2014, pp. 681–690

B. Zwiebach (2016) “Lecture 6,” Course Notes for MIT 8.04 Quantum Physics, Spring 2016.

Basil J. Hiley and Peter Van Reeth (2018) “Quantum trajectories: real or surreal?,” Entropy vol. 20, pp. 353 doi:10.3390/e20050353

Oliver Passon (2018) “On a common misconception regarding the de Broglie-Bohm theory,” Entropy vol. 20, no. 440. doi:10.3390/e20060440

Asher Yahalom (2018) “The fluid dynamics of spin,” Molecular Physics, April 2018, doi: 10.1080/00268976.2018.1457808. https://www.researchgate.net/publication/324512014, arXiv:1802:09331v1 [physics.flu-dyn] 3 Feb 2018

Siddhant Das and Detlef Dürr (2019) “Arrival time distributions of spin-1/2 particles,” Scientific Reports, https://doi.org/10.1038/s41598-018-38261-4

Siddhant Das, Markus Nöth, and Detlef Dür (2019) “Exotic Bohmian arrival times of spin-1/2 particles I—An analytical treatment,” arXiv:1901.08672v1 [quant-ph] 24 Jan 2019

2.5. Nonlinearity in the Bohmian mechanics:

To my surprise, I found that a form of non-linearity has been found to come up in the Bohmian mechanics too. I am sure it must have come as a surprise to many others too. [I will comment on this aspect quite some time later. For the time being, let me list some of the papers/presentations I’ve found so far.]

Sheldon Goldstein (1999) “Absence of chaos in Bohmian dynamics,” arXiv:quant-ph/9901005v1 6 Jan 1999

S. Sengupta, A. Poddar and P. K. Chattaraj (2000) “Quantum manifestations of the classical chaos in an undamped Duffing oscillator in presence of an external field: A quantum theory of motion study,” Indian Journal of Chemistry, vol. 39A, Jan–March 2000, pp. 316–322

A. Benseny, G. Albareda, A. S. Sanz, J. Mompart, and X. Oriols (2014) “Applied Bohmian mechanics,” arXiv:1406.3151v1 [quant-ph] 12 Jun 2014

Athanasios C. Tzemos (2016) “The mechanism of chaos in 3-D Bohmian trajectories,” Poster Presentation, https://www.researchgate.net/publication/305317081

Athanasios C. Tzemos (2018) “3-d Bohmian chaos: a short review,” Presentation Slides, RCAAM, Academy Of Athens

Athanasios C. Tzemos (2019) “Quantum entanglement and Bohmian Mechanics,” Presentation Slides 17 July 2019, RCAAM of the Academy of Athens

Klaus von Bloh (2020) “Bohm trajectories for the noncentral Hartmann potential,” Wolfram demonstration projects, https://www.researchgate.net/publication/344171771 (August 2020)

G. Contopoulos and A. C. Tzemos (2020) “Chaos in Bohmian quantum mechanics: a short review,” arXiv:2009.05867v1 [quant-ph] 12 Sep 2020

3. What happens to my new approach?

It was only yesterday that a neat thing struck me. Pending verification via simulations, it has the potential to finally bring together almost all of my research on the spinless particles. I’ve noted this insight in the hand-written journal (i.e. research notebook) that I maintain. I will be developing this idea further too. After all, Bohmians do study mainstream quantum mechanics and other interpretations, don’t they?

Due to the RSI, the simulations, however, will have to wait further. (The status is more or less the same. If I type for 2–3 hours, it’s easily possible that I can’t do much anything for the next 2–3 days.)

OK. Take care and bye for now.

A song I like:

(Hindi) देखा ना हाय रे सोचा ना (“dekhaa naa haay re sochaa naa”)
Singer: Kishore Kumar
Music: R. D. Burman
Lyrics: Rajinder Krishan

[Another song I used to love in my high-school days—who wouldn’t? … And, of course, I still do! A good quality audio I found is here [^]. I had not watched this movie until about a decade ago, on a CD (or may be on the TV). I’ve forgotten the movie by now. I don’t mind giving you the link for the video of this song; see here [^]. (In any case, it’s at least 3 orders of magnitude better than any so-called Lyrical Video Saregama has released for any song. The very idea of the Lyrical is, IMO, moronic.)]

# Python scripts for simulating QM, part 2: Vectorized code for the H atom in a 1D/2D/3D box. [Also, un-lockdown and covid in India.]

A Special note for the Potential Employers from the Data Science field:

Recently, in April 2020, I achieved a World Rank # 5 on the MNIST problem. The initial announcement can be found here [^], and a further status update, here [^].

All my data science-related posts can always be found here [^].

1.The first cut for the H atom in a 3D box:

The last time [^], I spoke of an enjoyable activity, namely, how to make the tea (and also how to have it).

Talking of other, equally enjoyable things, I have completed the Python code for simulating the H atom in a $3D$ box.

In the first cut for the 3D code (as also in the previous code in this series [^]), I used NumPy’s dense matrices, and the Python for” loops. Running this preliminary code, I obtained the following colourful diagrams, and twitted them:

H atom in a 3D box of 1 angstrom sides. Ground state (i.e., 1s, eigenvalue index = 0). All contours taken together show a single stationary state. Contour surfaces plotted with wiremesh. Plotted with Mayavi’s mlab.contour3d().

H atom in a 3D box of 1 angstrom sides. A ‘p’ state (eigenvalue index = 2). All contours taken together show a single stationary state. Contour surfaces with the Gouraud interpolation. Plotted with Mayavi’s mlab.contour3d().

H atom in a 3D box of 1 angstrom sides. A ‘p’ state (eigenvalue index = 2). All contours taken together show a single stationary state. Contour surfaces with wiremesh. Plotted with Mayavi’s mlab.contour3d().

H atom in a 3D box of 1 angstrom sides. Another ‘p’ state (eigenvalue index = 3). All contours taken together show a single stationary state. Contour surfaces with the Gourauad interpolation. Plotted with Mayavi’s mlab.contour3d().

OK, as far as many (most?) of you are concerned, the enjoyable part of this post is over. So, go read something else on the ‘net.

Coming back to my enjoyment…

2. Sparse matrices. Vectorization of code:

After getting to the above plots with dense matrices and Python “for” loops, I then completely rewrote the whole code using SciPy‘s sparse matrices, and put a vectorized code in place of the Python “for” loops. (As a matter of fact, in the process of rewriting, I seem to have deleted the plots too. So, today, I took the above plots from my twitter account!)

2.1 My opinions about vectorizing code, other programmers, interviews, etc:

Vectorization was not really necessary in this problem (an eigenvalue problem), because even if you incrementally build the FD-discretized Hamiltonian matrix, it takes less than 1 percent of the total execution time. 99 % of the execution time is spent in the SciPy library calls.

Python programmers have a habit of always looking down on the simple “for” loops—and hence, on any one who writes them. So, I decided to write this special note.

The first thing you do about vectorization is not to begin wondering how best to implement it for a particular problem. The first thing you do is to ask yourself: Is vectorization really necessary here? Ditto, for lambda expressions. Ditto, for list comprehension. Ditto for itertools. Ditto for almost anything that is a favourite of the dumb interviewers (which means, most Indian Python programmers).

Vectorized codes might earn you brownie points among the Python programmers (including those who interview you for jobs). But such codes are more prone to bugs, harder to debug, and definitely much harder to understand even by you after a gap of time. Why?

That’s because practically speaking, while writing in Python, you hardly if ever define C-struct like things. Python does have classes. But these are rather primitive classes. You are not expected to write code around classes, and then put objects into containers. Technically, you can do that, but it’s not at all efficient. So, practically speaking, you are almost always into using the NumPy ndarrays, or similar things (like Pandas, xarrays, dasks, etc.).

Now, once you have these array-like thingies, indexing becomes important. Why? Because, in Python, it is the design of a number of arrays, and the relationships among their indexing scheme which together “defines” the operative data structures. Python, for all its glory, has this un-removable flaw: The design of the data structures is always implicit; never directly visible through a language construct.

So, in Python, it’s the indexing scheme which plays the same part as the classes, inheritance, genericity play in C++. But it’s implicit. So, how you implement the indexing features becomes of paramount importance.

And here, in my informed opinion, the Python syntax for the slicing and indexing operations has been made unnecessarily intricate. I, for one, could easily design an equally powerful semantics that comes with a syntax that’s much easier on the eye.

In case some professionally employed data scientist (especially a young Indian one) takes an offence to my above claim: Yes, I do mean what I say above. And, I also know what I am talking about.

Though I no longer put it on my CV, once, in the late 1990s, I had implemented a Yacc-like tool to output table-driven parser for the LALR-1 languages (like Java and C++). It would take a language specification in the EBNF (Extended Backus-Noor Form) as the input file, and produce the tables for table-driven parsing of that language. I had implemented this thing completely on my own, looking just at the Dragon Book (Aho, Sethi, Ullman). I haven’t had a CS university education. So, I taught myself the compilers theory, and then, began straight implementing it.

I looked at no previous code. And even if I were to look at something, it would have been horrible. These were ancient projects, written in C, not in C++, and written using arrays, no STL containers like “map”s. A lot of hard-coding, pre-proc macros, and all that. Eventually, I did take a look at the others’ code, but it was only in the verification stage. How did my code fare? Well, I didn’t have to change anything critical.

I had taken about 8 months for this exercise (done part time, on evenings, as a hobby). The closest effort was by some US mountain-time university group (consisting of a professor, one or two post-docs, and four-five graduate students). They had taken some 4 years to reach roughly the same place. To be fair, their code had many more features. But yes, both their code and mine addressed those languages which belonged to the same class of grammar specification, and hence, both would have had the same parsing complexity.

I mention it all it here mainly in order to “assure” the Indian / American programmers (you know, those BE/BS CS guys, the ones who are right now itching to fail me in any interview, should their HR arrange one for me in the first place) that I do know a bit about it when I was talking about the actual computing operations on one hand, and the mere syntax for those operations on the other. There are a lot of highly paid Indian IT professionals who never do learn this difference (but take care to point out that your degree isn’t from the CS field).

So, my conclusion is that despite all its greatness (and I do actually love Python), its syntax does have some serious weaknesses. Not just idiosyncrasies (which are fine) but actual weaknesses. The syntax for slicing and indexing is a prominent part of it.

Anyway, coming back to my present code (for the H atom in the 3D box, using finite difference method), if the execution time was so short, and if vectorization makes a code prone to bugs (and difficult to maintain), why did I bother implementing it?

Two reasons:

1. I wanted to have a compact-looking code. I was writing this code mainly for myself, so maintenance wasn’t an issue.
2. In case some programmer / manager interviewing me began acting over-smart, I wanted to have something which I could throw at his face. (Recently, I ran into a woman who easily let out: “What’s there in a PoC (proof of concept, in case you don’t know)? Any one can do a PoC…” She ranted on a bit, but it was obvious that though she has been a senior manager and all, and lists managing innovations and all, she doesn’t know. There are a lot of people like her in the Indian IT industry. People who act over-smart. An already implemented vectorized code, especially one they find difficult to read, would be a nice projectile to have handy.

2.2 Sparse SciPy matrices:

Coming back to the present code for the H atom: As I was saying, though vectorization was not necessary, I have anyway implemented the vectorization part.

I also started using sparse matrices.

In case you don’t know, SciPy‘s and NumPy‘s sparse matrix calls look identical, but they go through different underlying implementations.

From what I have gathered, it seems safe to conclude this much: As a general rule, if doing some serious work, use SciPy’s calls, not NumPy’s. (But note, I am still learning this part.)

With sparse matrices, now, I can easily go to a $50 \times 50 \times 50$ domain. I haven’t empirically tested the upper limit on my laptop, though an even bigger mesh should be easily possible. In contrast, earlier, with dense matrices, I was stuck at at most at a $25 \times 25 \times 25$ mesh. The execution time too reduced drastically.

In my code, I have used only the dok_matrix() to build the sparse matrices, and only the tocsr(), and tocoo() calls for faster matrix computations in the SciPy eigenvalue calls. These are the only functions I’ve used—I haven’t tried all the pathways that SciPy opens up. However, I think that I have a pretty fast running code; that the execution time wouldn’t improve to any significant degree by using some other combination of calls.

2.3 A notable curiosity:

I also tried, and succeeded to a great degree, in having an exactly identical code for all dimensions: 1D, 2D, 3D, and then even, in principle, ND. That is to say, no “if–else” statements that lead to different execution paths depending on the dimensionality.

If you understand what I just stated, then you sure would want to have a look at my code, because nothing similar exists anywhere on the ‘net (i.e., within the first 10 pages thrown up by Google during several differently phrased searches covering many different domains).

However, eventually, I abandoned this approach, because it made things too complicated, especially while dealing with computing the Coulomb fields. The part dealing with the discretized Laplacian was, in contrast, easier to implement, and it did begin working fully well, which was when I decided to abandon this entire approach. In case you know a bit about this territory: I had to liberally use numpy.newaxis.

Eventually, I came to abandon this insistence on having only a single set of code lines regardless of the dimensionality, because my programmer’s/engineer’s instincts cried against it. (Remember I don’t even like the slicing syntax of Python?) And so, I scrapped it. (But yes, I do have a copy, just in case someone wants to have a look.)

2.4 When to use the “for” loops and when to use slicing + vectorization: A good example:

I always try to lift code if a suitable one is available ready made. So, I did a lot of search for Python/MatLab code for such things.

As far as the FD implementations of the Laplacian go, IMO, the best piece of Python code I saw (for this kind of a project) was that by Prof. Christian Hill [^]. His code is available for free from the site for a book he wrote; see here [^] for an example involving the finite difference discretization of the Laplacian.

Yes, Prof. Hill has wisely chosen to use only the Python “for” loops when it comes to specifying the IC. Thus, he reserves the vectorization only for the time-stepping part of the code.

Of course, unlike Prof. Hill’s code (transient diffusion), my code involves only eigenvalue computations—no time-stepping. So, one would be even more amply justified in using only the “for” loops for building the Laplacian matrix. Yet, as I noted, I vectorized everything in my code, merely because I felt like doing so. It’s during vectorization that the problem of differing dimensionality came up, which I solved, and then abandoned.

2.5 Use of indexing matrices:

While writing my code, I figured out that a simple trick with using index matrices and arrays makes the vectorization part even more compact (and less susceptible to bugs). So, I implemented this approach—indexing matrices and arrays.

“Well, this is a very well known approach. What’s new?” you might ask. The new part is the use of matrices for indexing, not arrays. Very well known, sure. But very few people use it anyway.

Again, I was cautious. I wrote the code, saw it a couple of days later again, and made sure that using indices really made the code easier to understand—to me, of course. Only then I decided to retain it.

By using the indexing matrices, the code indeed becomes very clean-looking. It certainly looks far better (i.e. easier to grasp structure) than the first lines of code in Prof. Hill’s “do_timestep” function [^].

2.6 No code-drop:

During my numerous (if not exhaustive) searches, I found that no one posts a 3D code for quantum simulation that also uses finite differences (i.e. the simplest numerical technique).

Note, people do post codes for 3D, but these are only for more complicated approaches like: FDTD (finite difference time domain), FEM, (pseudo)spectral methods, etc. People also post code for FDM, when the domain is 1D. But none posts a code that is both FD and 2D/3D. People only post the maths for such a case. Some rare times, they also post the results of the simulations. But they don’t post the 3D FDM code. I don’t know the reason for this.

May be there is some money to be made if you keep some such tricks all to yourself?

Once this idea occurred to me, it was impossible for me not to take it seriously. … You know that I have been going jobless for two years by now. And, further, I did have to invest a definite amount of time and effort in getting those indexing matrices working right so that the vectorization part becomes intuitive.

So, I too have decided not to post my 3D code anywhere on the ‘net for free. Not immediately anyway. Let me think about it for a while before I go, post my code.

3. Covid in India:

The process of unlocking down has begun in India. However, the numbers simply aren’t right for any one to get relaxed (except for the entertainment sections of the Indian media like the Times of India, Yahoo!, etc.).

In India, we are nowhere near turning the corner. The data about India are such that even the time when the flattening might occur, is not just hard to predict, but with the current data, it is impossible.

Yes, I said impossible. I could forward reasoning grounded in sound logic and good mathematics (e.g., things like Shannon’s theorem, von Neumann’s errors analysis, etc.), if you want. But I think to any one who really knows a fair amount of maths, it’s not necessary. I think they will understand my point.

Let me repeat: The data about India are such that even the time when the flattening might occur, is not just hard to predict, but with the current data, it is impossible.

India’s data show a certain unique kind of a challenge for the data scientist—and it definitely calls for some serious apprehension by every one concerned. The data themselves are such that predictions have to be made very carefully.

If any one is telling you that India will cross (or has already crossed), say, more than 20 lakh cases, then know that he/she is not speaking from the data, the population size, the social structures, the particular diffusive dynamics of this country, etc. He/she is talking purely from imagination—or very poor maths.

Ditto, if someone tells you that there are going be so many cases in this city or that, by this date or that, if the date runs into, say, August.

Given the actual data, in India, projections about number of cases in the future are likely to remain very tentative (having very big error bands).

Of course,  you may still make some predictions, like those based on the doubling rate. You would be even justified in using this measure, but only for a very short time-span into the future. The reason is that India’s data carry these two peculiarities:

1. The growth rate has been, on a large enough scale, quite steady for a relatively longer period of time. In India, there has been no exponential growth with a very large log-factor, not even initially (which I attribute to an early enough a lock-down). There also has been no flattening (for whatever reasons, but see the next one).
2. The number of cases per million population still remains small.

Because of 1., the doubling rate can serve as a good short-term estimator when it comes to activities like large-scale resource planning (but it would be valid only for the short term). You will have to continuously monitor the data, and be willing to adjust your plans. Yet, the fact is also that the doubling rate has remained steady long enough that it can certainly be used for short-term planning (including by corporates).

However, because of 2., everyone will have to revise their estimates starting from the third week of June, when the effects of the un-locking down begin to become visible (not just in the hospitals or the quarantine centres, but also in terms of aggregated numbers).

Finally, realize that 1. matters only to the policy-makers (whether in government or in corporate sectors).

What matters to the general public at large is this one single question: Have we turned around the corner already? if not, when will we do that?

The short answers are: “No” and  “Can’t Tell As of Today.”

In India’s case the data themselves are such that no data scientist worth his salt would be able to predict the time of flattening with any good accuracy—as of today. Nothing clear has emerged, even after 2.5 months, in the data. Since this sentence is very likely to be misinterpreted, let me explain.

I am not underestimating the efforts of the Indian doctors, nurses, support staff, police, and even other government agencies. If they were not to be in this fight, the data would’ve been far simpler to analyse—and far more deadly.

Given India’s population size, its poverty, its meagre medical resources, the absence of civic discipline, the illiteracy (which makes using symbols for political parties indispensable at the time of elections)… Given all such factors, the very fact that India’s data even today (after 2.5 months) still manages to remain hard to analyse suggests, to my mind, this conclusion:

There has been a very hard tussle going on between man and the virus so that no definitive trend could emerge either way.

There weren’t enough resources so that flattening could occur by now. If you kept that expectation to begin with, you were ignoring reality.

However, in India, the fight has been such that it must have been very tough on the virus too—else, the exponential function is too bad for us, and it is too easy for the virus.

The inability to project the date by which the flattening might be reached, must be seen in such a light.

The picture will become much clearer starting from two weeks in the future, because it would then begin reflecting the actual effects that the unlocking is producing right now.

So, if you are in India, take care even if the government has now allowed you to step out, go to office, and all that. But remember, you have to take even more care than you did during the lock-down, at least for the next one month or so, until the time that even if faint, some definitely discernible trends do begin to emerge, objectively speaking.

I sincerely hope that every one takes precautions so that we begin to see even just an approach towards the flattening. Realize, number of cases and number deaths increase until the flattening occurs. So, take extra care, now that the diffusivity of people has increased.

Good luck!

A song I like:

(Western, instrumental): Mozart, Piano concerto 21, k. 467, second movement (andante in F major).

Listen, e.g., at this [^] YouTube viedo.

[ I am not too much into Western classical, though I have listened to a fair deal of it. I would spend hours in UAB’s excellent music library listening to all sorts of songs, though mostly Western classical. I would also sometimes make on-the-fly requests to the classical music channel of UAB’s radio station (or was it a local radio station? I no longer remember). I didn’t always like what I listened to, but I continuing listening a lot anyway.

Then, as I grew older, I began discovering that, as far as the Western classical music goes, very often, I actually don’t much appreciate even some pieces that are otherwise very highly regarded by others. Even with a great like Mozart, there often are places where I can’t continue to remain in the flow of the music. Unknowingly, I come out of the music, and begin wondering: Here, in this piece, was the composer overtaken by a concern to show off his technical virtuosity rather than being absorbed in the music? He does seem to have a very neat tune somewhere in the neighbourhood of what he is doing here. Why doesn’t he stop tinkling the piano or stretching the violin, stop, think, and resume? I mean, he was composing music, not just blogging, wasn’t he?

The greater the composer or the tune suggested by the piece, the greater is this kind of a disappointment on my part.

Then, at other times, these Western classical folks do the equivalent of analysis-paralysis. They get stuck into the same thing for seemingly forever. If composing music is difficult, composing good music in the Western classical style is, IMHO, exponentially more difficult. That’s the reason why despite showing a definite “cultured-ness,” purely numbers-wise, most Western classical music tends to be boring. … Most Indian classical music also tends to be very boring. But I will cover it on some other day. Actually, one day won’t be enough. But still, this post is already too big…

Coming to the Western classical, Mozart, and the song selected for this time: I think that if Mozart were to do something different with his piano concerto no. 20 (k. 466), then I might have actually liked it as much as k. 467, perhaps even better. (For a good YouTube video on k. 466, see here [^].)

But as things stand, it’s k. 467. It is one of the rarest Western (or Eastern) classical pieces that can auto ride on my mind at some unpredictable moments; also one of the rare pieces that never disappoint me when I play it. Maybe that’s because I don’t play it unless I am in the right mood. A mood that’s not at all bright; a mood that suggests as if someone were plaintively raising the question “why? But why?”. (Perhaps even: “But why me?”) It’s a question not asked to any one in particular. It’s a question raised in the midst of continuing to bear either some tragic something, or, may be, a question raised while in the midst of having to suffer the consequences of someone else’s stupidity or so. … In fact, it’s not even a question explicitly raised. It’s to do with some feeling which comes before you even become aware of it, let alone translate it into a verbal question.  I don’t know, the mood is something like that. … I don’t get in that kind of a mood very often. But sometimes, this kind of a mood is impossible to avoid. And then, if the outward expression of such a mood also is this great, sometimes, you even feel like listening to it… The thing here is, any ordinary composer can evoke pathos. But what Mozart does is in an entirely different class. He captures the process of forming that question clearly, alright. But he captures the whole process in such a subdued manner. Extraordinary clarity, and extraordinary subdued way of expressing it. That’s what appeals to me in this piece… How do I put it?… It’s the epistemological clarity or something like that—I don’t know. Whatever be the reason, I simply love this piece. Even if I play it only infrequently.

Coming back to the dynamic k. 466 vs. the quiet, sombre, even plaintive k. 467, I think, the makers of the “Elvira Madigan” movie were smart; they correctly picked up the k. 467, and only the second movement, not others. It’s the second movement that’s musically extraordinary. My opinion, anyway…

Bye for now.

]