Convert multisphere to use MSTM

[barkls]

@AnnieStephenson pointed out that there is currently an upper limit of 20 spheres (specified in scfodim.for file). Is there a physical or computational reason for this hard limit?

[vnmanoharan]

The Multisphere model currently relies on Fortran 77 code. Fortran 77 doesn't do dynamic memory allocation for arrays. So what one has to do is to specify the largest size for the array ahead of time, which can be done by setting a few parameters in scfodim.for. We could increase the number of spheres, but that has an impact on memory utilization even for small numbers of spheres (since the arrays are allocated at the maximum size).

The longer-term solution is to switch to Mackowski's F90 MSTM code. We have a python wrapper around this code, but it would need some work to make it compatible with holopy. In particular, we would need to:

• modify the Fortran code to output the amplitude scattering matrix rather than the scattering matrix
• modify the code to work with F2py, so that we can call functions directly from holopy rather than going through a pipe or input/output from a file
• get permission from Mackowski to distribute the modified code under the terms of the GPLv3

[barkls]

I guess the short term question then becomes: What is the memory load per sphere? It looks to me like it's at least around 10 MB just from skimming the code. We could do some memory usage tests to decide if 20 spheres is a reasonable limit for modern computers, or we can just leave it until the longer-term solution is implemented.

[vnmanoharan]

I think we can leave it up to the user. Anyone who needs a larger number of spheres can simply fork the repository, make the change to scsmfo.dim, and then pull new revisions from the main branch if they want to keep up with changes. There should probably be a note somewhere in the documentation, though. @tdimiduk and @AnnieStephenson, any comments?

[tdimiduk]

We set the 20 sphere number back when memories were smaller so pushing it to 40 or 50 might be reasonable. I agree with @vnmanoharan that it makes sense to include instructions on changing the value in the docs, or even provide a script to change the value and recompile (that might actually be easier than describing it in enough detail for a non expert to do the change).

[vnmanoharan]

OK, then I suggest doing a few tests and seeing how memory use changes as we increase the value. Let's set an upper bound of say 1 GB of memory (so that users with 32-bit machines can still operate the code). @barkls, can you look at this?

[tdimiduk]

@anna-wang did a conversion to amplitude scattering matrices of another t-matrix code for non spherical scatterers in PR #7 I don't know if that is useful for a template for converting the Mackowski F90 MSTM code?

Also, we could just teach holopy to use the scattering matrices output by the code without converting them. Not sure what the relative effort levels are.

@barkls can you get in touch with Mackowski to see about permission? You can look at the permission note jerome wrote for a template of how to ask.

[vnmanoharan]

The default output is the scattering matrix, which gives the intensity and polarization state. So we do need the amplitude scattering matrices, because we need the fields. It's not too hard to do that, since the code calculates them in the course of calculating the scattering matrix. We also need the near-field results, which shouldn't be too hard to get. But I don't know how much effort is involved in converting the code so that python can talk directly to the compiled Fortran functions.

[barkls]

I think updating to F90 MSTM code is a definite case of scope creep for the Holopy 3 release. However, if you are both confident that is what we should do, then I can start working on it.

Otherwise, I think it makes more sense to establish a new default max number of spheres and include instructions in the docs for how to increase the limit.

[vnmanoharan]

I agree that for 3.0 we should simply investigate increasing the number of spheres (see my comment above). But I also think it's a good idea to get in touch about permission to modify and redistribute MSTM, so that we can start working on another branch where MSTM replaces the F77 code. We are not prioritizing development of that branch, but without permission, we can't even start working on it.

[barkls]

Sounds good, I will do both of these tasks.

[barkls]

After playing around with this for a bit, I have determined that the obstacle isn't total memory usage. Rather, the fortran code becomes unstable if either npd (max number of spheres) or notd (max cluster expansion order) parameters are increased in scfodim.for. It will run properly for some sphere clusters occasionally, but usually throws one of several error types having to do with array indexing and memory allocation. Interestingly, subsequent runs of identical code does not always produce the same errors.

I will get in touch with Mackowski about the F90 code, but I'm closing this issue in the meantime.

[vnmanoharan]

Note that Mackowski has now released the MSTM codes under an MIT license. See here: https://github.com/dmckwski/MSTM