Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Anisotropic Ewald code #2182

Open
wants to merge 13 commits into
base: develop
Choose a base branch
from

Conversation

jtkrogel
Copy link
Contributor

@jtkrogel jtkrogel commented Jan 3, 2020

This PR contains an adaptive anisotropic Ewald code, similar to the reference code previously developed.

The anisotropicEwald class can operate in two modes:

  1. Perform Ewald sums for all particle pairs while building up adaptive anisotropic grids until the sum converges. The grid stops growing in directions that do not contribute. The maximum grid size in each direction is recorded across all pairs in the sum (and subsequent evaluations if desired).

  2. Perform Ewald sums over a fixed rectangular grid that was determined during the adaptive phase. Grid sums of this type are much faster, and are at least as accurate as those determined during the adaptive phase.

This code was used to test the anisotropic nature of the Ewald sum in the quasi-2D graphene case for cells with vacuum ranging from 10-200 Bohr. For fixed accuracy, more points are needed in the z direction of the Ewald sum, while currently QMCPACK sums over a sphere (3x more points for z=40, 5x for z=200).

The resulting adaptive anisotropic grids contained 0.5 to 3 times the number of points used by QMCPACK by default in the k-space sum. In all cases tested (graphene, NiO) the adaptive code significantly outperformed QMCPACK's optimized breakup and Ewald3D implementations while requiring a similar number of points.

For this reason, I think we should consider replacing QMCPACK's LR code with an optimized version of the code in this PR.

@qmc-robot
Copy link

Can one of the admins verify this patch?

@jtkrogel jtkrogel changed the title Anisotropic ewald code Anisotropic Ewald code Jan 3, 2020
Copy link
Contributor

@ye-luo ye-luo left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A bit more documentation will be helpful for readers.

namespace qmcplusplus
{

using namespace ewaldref;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"using namespace" in header files may propagate very far away. I think all the functions and classes can be put inside the ewaldref name space.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When adding new code to the source, how do you determine what needs separate namspace(s) and what does not? Most of the code simply falls within the qmcplusplus namespace (which really should be qmcpack since the name of the code changed from QMC++ long ago...).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are two topics.

  1. header file "using namespace" is a bad practice. cpp file is fine. "using namespace" propagates in header files but not in cpp files.
  2. There is no strict rule whether extra namespace should be used or not. The original EwaldRef.h defines a few datatypes directly in the qmcplusplus name space. They can propagate and can be misused or cause potential name collisions. So I put most of the content of EwaldRef.h in ewaldref namespace within qmcplusplus namespace as a protection. Here the added code is quite standalone and closely related to EwaldRef then putting it in the same namespace as EwaldRef seems appropriate.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm fine with this, though the new code is not going to be used as reference. If it has a future, it will be in the main application, otherwise it will be deleted.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you intend to implement a class derived from OperatorBase, then don't put the class in ewaldref namespace. Instead, you can either put ewaldref:: for using any piece from EwaldRef or "using ewaldref::ABC" in the class to bring in ABC to the new class.

@@ -416,7 +416,7 @@ real_t ewaldEnergy(const RealMat& a, const PosArray& R, const ChargeArray& Q, re

#pragma omp parallel for reduction(+ : ve)
for(size_t n=0; n<Npairs; ++n)
ve += qq[n]*ewaldSum(rr[n],a,tol/qq[n]);
ve += qq[n]*ewaldSum(rr[n],a,tol/qq[n]/Npairs);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If possible, please run clang-format on both files.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not time efficient. I'm waiting for the automatic part of #794 to be implemented.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The automatic part has been implemented; see https://github.com/QMCPACK/qmcpack/wiki/Source-formatting

This was a decision made as a project to not have a robot do formatting without supervision, but to provide various methods for different workflow preferences that everyone can choose for themselves, from assisted manual to full automatic.

}


real_t madelungConstant()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please revise function names in this file following our function naming rules. For example here, calculateMadelungConstant(). This function also doesn't modify any class member, so please mark this function const.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The naming already follows the coding conventions, as in the manual:

27.3.5 (Member) function names
Function names should start with a lowercase character and have a capital letter for each new word.

Adding calculate to every function that calculates something is IMO trite and lexically wasteful.

Will add const.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

madelungConstant() is not clear whether it to access a stored value in the class or run a calculation. When I name functions, evaluate/calculateXXX() means it can be costly. getXXX() means cheap. Probably in your case, it doesn't matter much.

int_t km = 0;
IntVec iv;
// Add the value at the origin, if requested.
if (zero)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if we define a member function addZeroContribution() to both T& function, there is no need of the zero switch.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The term should not need to know how it is being summed.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ignore my last comment. zero and function need to be controlled independently.

@ye-luo
Copy link
Contributor

ye-luo commented Jan 6, 2020

Leave some notes of the understanding from my last optimization attempt.
The k-space term is expensive when calculation sincos in exp(G (Ri-Rj)), a direct computation cost Ng x Natom x Natom. We can separate this term into exp(Gs Ri) x exp(-Gs Rj) and the cost reduces to 2 x Ng x Natom.
The r-space term is expensive in erfc(|Rs -(Ri - Rj)|). It costs Nr x Natom x Natom. I don't have a way to reduce its scaling.
The kappa choice we have in the current code assumes the same scaling for both K and R space computation. I saw more or less the same time spent in both parts.
To have a fully balanced K and R space computation, we need to take into account
1. the single function call cost of sincos() and erfc()
2. the scaling difference between K and R space computation if we use the separate form to reduce computation.
moved to #2176 which is more relavent.

@jtkrogel
Q. How was the performance comparison made? just Ion-ion coulomb calculation or electron-electron term as well?

I misunderstood "outperform" for computational cost.

@jtkrogel
Copy link
Contributor Author

jtkrogel commented Jan 6, 2020

I'm not sure what you are asking for in terms of my documenting your last optimization attempt. This code is not presently visible, and optimizations of that type are already well documented in the manual, e.g. in equations 26.11-26.13.

The code in this PR is unoptimized as yet. When I mentioned "outperformed" above, I meant in terms of accuracy. I measured the potential for similar efficiency based on the number of points that need to be summed (which is about equivalent) for the ion-ion sum.

In the next step (after this PR), I will attempt to replace QMCPACK's use of the opt breakup with this code for the e-e coulomb part and compare speed. Optimization may follow at that point (using OMP, inlining and reordering sums, etc).

@jtkrogel
Copy link
Contributor Author

jtkrogel commented Jan 6, 2020

I've added new namespace scoping, the requested const function and light function documentation.

I'm willing to improve these further in the later event that this code is adapted and used in the main application. For now, the code is exploratory. Good enough for now?

ye-luo
ye-luo previously approved these changes Jan 6, 2020
Copy link
Contributor

@ye-luo ye-luo left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM

@grahamlopez
Copy link
Contributor

I'm willing to improve these further in the later event that this code is adapted and used in the main application. For now, the code is exploratory. Good enough for now?

Please be careful with this kind of thinking. These PR's are the main gate to code getting merged into develop. Once in develop, intentional steps have to be taken by a human to keep it from eventually going into production. Of course intentions can be good, but we know how plans/priorities change, tasks shift to different people, code gets forgotten/neglected, etc. If this code works well enough to be merged now, why shouldn't it deserve the same kind of care as any other? The amount of work is the same (or probably less) if it happens now especially since it currently has attention vs. when it possibly gets revisited later for separate reasons.

@prckent prckent self-requested a review January 6, 2020 16:27
Copy link
Contributor

@prckent prckent left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This code should not be merged until we have had a full discussion around the criteria for updating the current Coulomb code.

Adding some comments here is on my to do list for today.

@jtkrogel
Copy link
Contributor Author

jtkrogel commented Jan 6, 2020

N.B. This code is not used anywhere outside of its own header file right now.

I'm happy to retract this PR and keep on exploring, but I thought making a waypoint for everyone else to see would be useful.

@jtkrogel
Copy link
Contributor Author

jtkrogel commented Jan 6, 2020

I'll make an issue elaborating where I am in the exploration process and where it might lead in later steps.

@jtkrogel
Copy link
Contributor Author

jtkrogel commented Jan 6, 2020

See #2185.

@jtkrogel jtkrogel changed the title Anisotropic Ewald code [WIP] Anisotropic Ewald code Jan 9, 2020
@jtkrogel
Copy link
Contributor Author

jtkrogel commented Jan 9, 2020

Will update the code here occasionally as #2185 unfolds.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants