Resurrect Ewald2D

[Paul-St-Young]

Please review the developer documentation
on the wiki of this project that contains help and requirements.

Proposed changes

Resurrect true 2D Ewald sum without using OHMMS_DIM.

What type(s) of changes does this code introduce?

• New feature

Does this introduce a breaking change?

• No

What systems has this change been tested on?

Workstation

Checklist

Update the following with a yes where the items apply. If you're unsure about any of them, don't hesitate to ask. This is
simply a reminder of what we are going to look for before merging your code.

• Yes. This PR is up to date with current the current state of 'develop'
• No. Code added or changed in the PR has been clang-formatted
• Yes. This PR adds tests to cover any new code, or to catch a bug that is being fixed
• No. Documentation has been added (if appropriate)

[prckent]

Paul -

Thanks for looking at this. It would be good to have a proper 2D evaluation... if this is a proper implementation. How thoroughly have you checked this code and its outputs? I think it is safe to say that shortcuts were made in the historical Coulomb code, so before we can advertise 2D Coulomb as working and add examples etc. I think we need to be very thorough with testing, convergence checks, validation etc.

[Paul-St-Young]

@prckent I looked at Madelung energy of square and triangular lattices using the unit cell and select supercells.
The potential energy of the perfect crystal checks out to 8 digits.
I also see reasonable numbers in VMC and DMC runs (which will require a few more PRs to enable).

That said, I'm not comfortable with this becoming a production-level feature right now.
I do want to get to that point in the near future (after more integration tests).

For now, this feature can be enable only via an undocumented flag <parameter name="ndim"> 2 </parameter> in <simulationcell>. Otherwise, the usual 3D code path is taken, so current users should be unaffected.

[prckent]

I looked at Madelung energy of square and triangular lattices using the unit cell and select supercells. The potential energy of the perfect crystal checks out to 8 digits. That is a relief to hear. Can you add some of these to test2D? It is good to capture these checks and establish tests that anyone fiddling with the code in future will have to keep working. An anisotropic case would be good as well.

[Paul-St-Young]

I'll need to figure out how to add these to the unit test.
Attached are the inputs and outputs in case I lose them ...

tb11f2_vmad.zip

[Paul-St-Young]

@prckent I added all my checks to the unit test.
This PR is ready for review.

I'm not sure why the CI cannot build the code.
I am unable to reproduce the CI failures using Clang 11, GCC 11.2.0, or GCC 10.3.0.

[jtkrogel]

I would like to discuss a little about the use of ndim in the input. As I understand, the physical systems we can simulate right now are fully 3D, while Ewald2D is meant to apply to systems that have two periodic dimensions. Is this correct?

If so, I don't think we should use ndim in the <simulationcell/> input. Rather, the use of Ewald2D should be triggered by the boundary conditions, i.e. bconds==ppp => Ewald3D while bconds==ppn/pnp/npp =>Ewald2D, and bconds==pnn/npn/nnp => abort.

As a follow-on, does the current code support ppn/pnp/npp or just ppn?

[prckent]

Rather, the use of Ewald2D should be triggered by the boundary conditions, i.e. bconds==ppp => Ewald3D while bconds==ppn/pnp/npp =>Ewald2D, and bconds==pnn/npn/nnp => abort. Agree. Note that there is no need to expose an option here - if someone wants to reproduce a 3D ion-ion energy they can use ppp. And we should definitely abort if someone tries pnn. This would be consistent with the general plan of having sensible default options for users and avoid/reduce the number of input parameters. Note that some ion-ion Coulomb checks based on the input wavefunctions might need bypassing for 2D, e.g. QE inputs. I think PySCF has proper low dimensional Coulomb handling now.

[Paul-St-Young]

@jtkrogel Yes, Ewald2D is meant for a system that's periodic only in 2 directions.
I hard-coded for ppn, so this code is not correct for pnp and npp conditions.

I don't think we can have ppn directly trigger Ewald2D, because one could be simulating a quasi-2D slab.
In that case, the EwaldHandlerQuasi2D should be used instead.

[jtkrogel]

Could you explain more? I don't understand the distinction (2D/quasi-2D) if all systems in question exist in a 3D space.

[Paul-St-Young]

By true 2D, I mean that all particles stay exactly on the same plane.
By quasi 2D, I mean that particles are much more confined in one of the spatial direction than the others, but can still move around a little bit along that direction. This out-of-plane movement is completely forbidden in the true 2D case.

[Paul-St-Young]

Do you have data to show 3D expanding converge to quasi2D value?

No, I do not have a 3D input with expanding z direction.

[Paul-St-Young]

On second thought, 3D ppn with large z will not match the quasi 2D value because they have different neutralizing backgrounds.
The 3D positive background charge density for an N-electron simulation cell is +e*N/V, i.e. \sum_i +e/V uniformly spread out over the entire cell volume V, whereas the quasi 2D background charge density is \sum_i +e/A * \delta(z-z_i), i.e. one plane of 2D positive charge at the z height of each electron. A is the 2D cell area.

[PDoakORNL]

It looks like my previous review was addressed where possible for this PR except for this.

[PDoakORNL]

If you want behavior to be preserved it should be tested, sometimes you need to be as clever in test code as you are in the actual implementation. You should be able to devise a synthetic case where this is true.

From the standpoint of the exception, putting the failing value in the app_log and it not actually being the alpha value means digging back into the implementation of the handlers for anyone that runs into this. Use a string stream or some other method and get the bad value into the exception's what string. You might need to indicate what the relationship between the short range part of the ewald being > 1e-6 at the WSR and apha being too small. Is this message valid for other choices of LRHANDLER? That is what the terribly name AA is right?

[ye-luo]

Could you also document the new option value in the manual?

[Paul-St-Young]

@PDoakORNL I clarified the exception message and added a test with CHECK_THROWS.
Turns out I just had to set a very small LR_dim_cutoff to get a violation of my condition.
The short range part being < 1e-6 at the WSR should apply to any LRHandler.
Although the message about alpha is specific to the Ewald handler.

[Paul-St-Young]

@ye-luo I added documentation for LR_handler in the manual.

[ye-luo]

Only minor things left.

Not directly related to this PR. I noticed quasi2D doesn't seem to be connected in PBCAB.

[ye-luo]

LGTM

[ye-luo]

Test this please

[ye-luo]

@PDoakORNL OK to merge?

[prckent]

It think that this meets easily meet the merge criteria now, and it can clearly be refined further as/if needed. Thanks for enduring the review Paul. I think that the code is better for it.

Keep things moving. Can be revisited in future if not considered sufficiently addressed.

[PDoakORNL]

@Paul-St-Young satsified my important concerns.