m-a-d-n-e-s-s / madness

Multiresolution Adaptive Numerical Environment for Scientific Simulation
GNU General Public License v2.0
180 stars 61 forks source link

Periodic `CoulombOperator` Gives Garbage #535

Open JonathonMisiewicz opened 5 months ago

JonathonMisiewicz commented 5 months ago

@evaleev and I are working on periodic code, and I find that values of the electric potential changed after rebasing onto a more recent version of MADNESS. The current results appear to have an additional "contribution" that is pure garbage (e.g., not even monotonic as you get away from the center of the sum-of-Gaussian charge distribution).

The problem is this line of code added during #498. The CoulombOperator required the last argument to truncate_periodic_expansion to be true, not false, prior to the PR. Manually changing that argument reproduces our electric potential prior to rebasing onto #498, but doing that would break many other operators.

@fbischoff, I noticed that you generalized some operators in operator.h to use your new function but left others the same. Is there an underlying reason? I want to make sure that my fixes to the CoulombOperator don't interfere with anything you were trying to do.

fbischoff commented 5 months ago

Hi all,

the line in question I have just copy-pasted, without knowing what it actually does, esp. as the periodic code seemed disfunctional.

The idea behind the generalization of the operator is to handle all operators of the following forms:

Coulomb: 1/r12 BSH: exp(-r12)/r12 Slater: exp(-r12) F12: 1-exp(-r12) F12^2: (1-exp(-r12))^2 = 1-2 exp(-r12) + exp(-2 r12) FG: F12 *Coulomb

by using the same expansion exponents (they depend on eps only, but not on the exponents of the operator), and add their coefficients (they depend on the actual form of the operator). That way all these operators are simply represented as a single operator.

So if I understand correctly, the parameter that breaks your code either does nothing or truncates the expansion, based on the assumption of decreasing exponents of the expansion. That should still be the case with the exception of F12 (and related operators), where I represent the “1” as a single term in the operator expansion with a tiny exponent: see gfit.h:660

Hope this helps.

Best, Florian

Florian Bischoff Humboldt-Universitaet zu Berlin, Institut fuer Chemie Unter den Linden 6, 10099 Berlin phone: +49-30-2093-82708 office: Brook-Taylor-Str 2, room 3'322 @.*** https://tinyurl.com/bischofflab

On May 22, 2024, at 2:42 PM, Jonathon Misiewicz @.***> wrote:

@evaleev and I are working on periodic code, and I find that values of the electric potential changed after rebasing onto a more recent version of MADNESS. The current results appear to have an additional "contribution" that is pure garbage (e.g., not even monotonic as you get away from the center of the sum-of-Gaussian charge distribution). The problem is this line of code added during #498. The CoulombOperator required the last argument to truncate_periodic_expansion to be true, not false, prior to the PR. Manually changing that argument reproduces our electric potential prior to rebasing onto #498, but doing that would break many other operators. @fbischoff, I noticed that you generalized some operators in operator.h to use your new function but left others the same. Is there an underlying reason? I want to make sure that my fixes to the CoulombOperator don't interfere with anything you were trying to do. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>

JonathonMisiewicz commented 5 months ago

The periodic code is certainly dysfunctional for our purposes, and we're hoping to make it less so.

That helps, thank you. As a follow-up, why did BSHOperator3D get updated but not BSHOperatorPtr3D? Is there a reason for this, or was it simply a matter of convenience?

The trouble in our case is that Gaussians which should be discarded are not, and I'd also need to look into why they need to be discarded in some cases but not others... It could very well be a defect in the original code.

fbischoff commented 5 months ago

The periodic code is certainly dysfunctional for our purposes, and we're hoping to make it less so. I’m looking forward to that!

That helps. As a follow-up, why did BSHOperator3D get updated but not BSHOperatorPtr3D? Is there a reason for this, or was it simply a matter of convenience?

no real reason, I tried to keep legacy code functional.

Best, Florian

Florian Bischoff Humboldt-Universitaet zu Berlin, Institut fuer Chemie Unter den Linden 6, 10099 Berlin phone: +49-30-2093-82708 office: Brook-Taylor-Str 2, room 3'322 @.*** https://tinyurl.com/bischofflab

On May 22, 2024, at 9:16 PM, Jonathon Misiewicz @.***> wrote:

The periodic code is certainly dysfunctional for our purposes, and we're hoping to make it less so. That helps. As a follow-up, why did BSHOperator3D get updated but not BSHOperatorPtr3D? Is there a reason for this, or was it simply a matter of convenience? The trouble in our case is that Gaussians which should be discarded are not, and I'd also need to look into why they need to be discarded in some cases but not others... It could very well be a defect in the original code. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>

gifornl commented 5 months ago

Did you set the function and operator to be periodic in the heading?

George

On Thu, May 23, 2024 at 3:18 AM Florian Bischoff @.***> wrote:

The periodic code is certainly dysfunctional for our purposes, and we're hoping to make it less so. I’m looking forward to that!

That helps. As a follow-up, why did BSHOperator3D get updated but not BSHOperatorPtr3D? Is there a reason for this, or was it simply a matter of convenience?

no real reason, I tried to keep legacy code functional.

Best, Florian

Florian Bischoff Humboldt-Universitaet zu Berlin, Institut fuer Chemie Unter den Linden 6, 10099 Berlin phone: +49-30-2093-82708 office: Brook-Taylor-Str 2, room 3'322 @.*** https://tinyurl.com/bischofflab

On May 22, 2024, at 9:16 PM, Jonathon Misiewicz @.***> wrote:

The periodic code is certainly dysfunctional for our purposes, and we're hoping to make it less so. That helps. As a follow-up, why did BSHOperator3D get updated but not BSHOperatorPtr3D? Is there a reason for this, or was it simply a matter of convenience? The trouble in our case is that Gaussians which should be discarded are not, and I'd also need to look into why they need to be discarded in some cases but not others... It could very well be a defect in the original code. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/m-a-d-n-e-s-s/madness/issues/535#issuecomment-2126743324, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACQQMOIBILPTCSTNTAFYYG3ZDW663AVCNFSM6AAAAABIDSAATCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCMRWG42DGMZSGQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

JonathonMisiewicz commented 5 months ago

Yes. The same application behaved differently before/after #498, and I can restore the old behavior by restoring an argument flag to how it was pre-#498, so I doubt my application code is the issue.

The relevant part of our application code looks like

    BoundaryConditions<3> bc_periodic(BC_PERIODIC);
    FunctionDefaults<3>::set_bc(bc_periodic);

    auto r000 = std::vector<coord_3d>{coord_3d(std::array<double,3>{0.0, 0.0, 0.0})};
    // g*_000_func3d are standard functions to evaluate a gaussian
    Function<double, 3> gdiffuse_000 = FunctionFactory<double, 3>(world).initial_level(4).special_points(r000).f(gdiffuse_000_func3d);
    printf("projected gdiffuse_000\n");
    Function<double, 3> gtight_000 = FunctionFactory<double, 3>(world).initial_level(4).special_points(r000).f(gtight_000_func3d);
    auto gdiffuse = gdiffuse_000;
    auto gtight = gtight_000;
    auto rho = gdiffuse - gtight;

    SeparatedConvolution<double, 3> pop = CoulombOperator(world, 1e-10, eps, bc_periodic);
    Function<double, 3> V_periodic = apply(pop, rho);
JonathonMisiewicz commented 4 weeks ago

The picture is a little clearer now, so I'll try to clean this issue up soon.

JonathonMisiewicz commented 4 weeks ago

@fbischoff

Another clarification regarding the recent refactor. Was the new constructor was added for backwards-compatibility, so it therefore should be moved away from, or is that just an unfortunate copy-paste artifact?

I need to know what constructor to use when I'm repairing the operator construction.

fbischoff commented 4 weeks ago

Hi Jonathon,

I’m not quite sure which constructor which copy-paste artifact you mean. My commit you’re referring to was quite (too?) large.. Maybe you can give me some specific line or so?

Best, Florian

Florian Bischoff Humboldt-Universitaet zu Berlin, Institut fuer Chemie Unter den Linden 6, 10099 Berlin phone: +49-30-2093-82708 office: Brook-Taylor-Str 2, room 3'322 @.*** https://tinyurl.com/bischofflab

On Oct 9, 2024, at 5:42 PM, Jonathon Misiewicz @.***> wrote:

@fbischoff https://github.com/fbischoff Another clarification regarding the recent refactor. Was the new constructor was added, so it therefore should be moved away from, or is that just an unfortunate copy-paste artifact?

https://github.com/m-a-d-n-e-s-s/madness/pull/498/files#diff-22333d357a4becfe78caa180594a5472c9adc7f780585e267719f78261b0c497R1007

I need to know what constructor to use when I'm repairing the operator construction.

— Reply to this email directly, view it on GitHub https://github.com/m-a-d-n-e-s-s/madness/issues/535#issuecomment-2402689276, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABZXPPGCI75ABVJQZRNWNDTZ2VFG7AVCNFSM6AAAAABIDSAATCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBSGY4DSMRXGY. You are receiving this because you were mentioned.

JonathonMisiewicz commented 3 weeks ago

Ah, I didn't realize that GitHub's link-to-specific-line feature doesn't work for files with diffs too large to show...

operator.h:1007

Note that the same comment is on line 1020.

Hi Jonathon, I’m not quite sure which constructor which copy-paste artifact you mean. My commit you’re referring to was quite (too?) large.. Maybe you can give me some specific line or so? Best, Florian Florian Bischoff Humboldt-Universitaet zu Berlin, Institut fuer Chemie Unter den Linden 6, 10099 Berlin phone: +49-30-2093-82708 office: Brook-Taylor-Str 2, room 3'322 @. https://tinyurl.com/bischofflab On Oct 9, 2024, at 5:42 PM, Jonathon Misiewicz @.> wrote: @fbischoff https://github.com/fbischoff Another clarification regarding the recent refactor. Was the new constructor was added, so it therefore should be moved away from, or is that just an unfortunate copy-paste artifact? https://github.com/m-a-d-n-e-s-s/madness/pull/498/files#diff-22333d357a4becfe78caa180594a5472c9adc7f780585e267719f78261b0c497R1007 I need to know what constructor to use when I'm repairing the operator construction. — Reply to this email directly, view it on GitHub <#535 (comment)>, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABZXPPGCI75ABVJQZRNWNDTZ2VFG7AVCNFSM6AAAAABIDSAATCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBSGY4DSMRXGY. You are receiving this because you were mentioned.

fbischoff commented 3 weeks ago

Hi Jonathon,

I see. Comment line 1007 is a copy-paste artifact, sorry for that..

Best, Florian


Florian Bischoff Humboldt-Universitaet zu Berlin, Institut fuer Chemie Unter den Linden 6, 10099 Berlin phone: +49-30-2093-82708 office: Brook-Taylor-Str 2, room 3'322 @.*** https://tinyurl.com/bischofflab

On Oct 10, 2024, at 3:22 PM, Jonathon Misiewicz @.***> wrote:

Ah, I didn't realize that GitHub's link-to-specific-line feature doesn't work for files with diffs too large to show...

operator.h:1007 https://github.com/m-a-d-n-e-s-s/madness/pull/498/files#diff-22333d357a4becfe78caa180594a5472c9adc7f780585e267719f78261b0c497 Note that the same comment is on line 1020.

Hi Jonathon, I’m not quite sure which constructor which copy-paste artifact you mean. My commit you’re referring to was quite (too?) large.. Maybe you can give me some specific line or so? Best, Florian Florian Bischoff Humboldt-Universitaet zu Berlin, Institut fuer Chemie Unter den Linden 6, 10099 Berlin phone: +49-30-2093-82708 office: Brook-Taylor-Str 2, room 3'322 @. https://tinyurl.com/bischofflab On Oct 9, 2024, at 5:42 PM, Jonathon Misiewicz @.> wrote: @fbischoff https://github.com/fbischoff https://github.com/fbischoff Another clarification regarding the recent refactor. Was the new constructor was added, so it therefore should be moved away from, or is that just an unfortunate copy-paste artifact? https://github.com/m-a-d-n-e-s-s/madness/pull/498/files#diff-22333d357a4becfe78caa180594a5472c9adc7f780585e267719f78261b0c497R1007 … <x-msg://58/#> I need to know what constructor to use when I'm repairing the operator construction. — Reply to this email directly, view it on GitHub <#535 (comment) https://github.com/m-a-d-n-e-s-s/madness/issues/535#issuecomment-2402689276>, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABZXPPGCI75ABVJQZRNWNDTZ2VFG7AVCNFSM6AAAAABIDSAATCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBSGY4DSMRXGY. You are receiving this because you were mentioned.

— Reply to this email directly, view it on GitHub https://github.com/m-a-d-n-e-s-s/madness/issues/535#issuecomment-2405083806, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABZXPPD33ZE64KHSD54EX43Z2Z5THAVCNFSM6AAAAABIDSAATCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBVGA4DGOBQGY. You are receiving this because you were mentioned.