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

make two-body Jastrow work in 2D #4289

Merged
merged 16 commits into from
Nov 2, 2022
Merged

Conversation

Paul-St-Young
Copy link
Contributor

Proposed changes

Follow up on #3868, enable Jastrow in 2D.
The WavefunctionTester was modified to run in 2D (in addition to 3D) to test these changes.

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

@Paul-St-Young Paul-St-Young self-assigned this Oct 21, 2022
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.

Could you please add a unit test. WaveFunctionTest is not the optimal way of making a test because it cannot provide enough safeguard. If the value, gradients are set zero, WaveFunctionTest cannot catch it.

@prckent prckent self-requested a review October 21, 2022 14:28
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.

Thanks Paul. Echoing Ye's comments, please stay away from wftester for formal testing. We have established routes for that via unit testing (by far preferred), followed by integration testing.

Background: wftester is a historical clearly not looked-after part of the code that we should not be increasing dependencies on. Other than this PR from you and a PR from me, there are no existing uses or tests. Perhaps we'll rewrite it, remove it, or trim it. Some days I have been tempted to just delete it, since it would be faster to start afresh with clean code. If there is something that is unclear how to write from a unit test please ask.

@Paul-St-Young
Copy link
Contributor Author

Paul-St-Young commented Oct 21, 2022

@ye-luo @prckent The WavefunctionTester is by far the easiest route to test particle and parameter derivatives of a complicated wavefunction. When I doubt the correctness of a run, it's easy replace the <qmc method="vmc/dmc"> section with <qmc method="wftest"> for a quick test. I know that @rcclay and my self have made many uses of this feature, for example when debugging backflow. @markdewing also made the FiniteDifference class fairly clean.

WavefunctionTester is also my go-to when developing a new wavefunction component, for example when resurrecting 2D electron gas orbitals #4084.
Since we have no other facility for generally checking the correctness of analytical derivatives (either for single-particle orbital, or for the many-body wavefunction), I don't see why we should remove this feature.
Having to re-code finite-difference for every unit test is not productive IMO.

I am happy to put in a simple regression unit test, but I would like to keep the wftest to check for correctness.

@@ -223,7 +224,7 @@ std::unique_ptr<WaveFunctionComponent> RadialJastrowBuilder::createJ2(xmlNodePtr
#if OHMMS_DIM == 1
RealType dim_factor = 1.0 / (OHMMS_DIM + 1);
#else
RealType dim_factor = (ia == ib) ? 1.0 / (OHMMS_DIM + 1) : 1.0 / (OHMMS_DIM - 1);
RealType dim_factor = (ia == ib) ? 0.5 / (ndim - 1) : 1.0 / (ndim - 1);
Copy link
Contributor

Choose a reason for hiding this comment

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

For 2 dimensions, the factor for the case ia==ib changes from 1/3 to 1/2. Is that change intended?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, this is an intentional change.
The same-spin cusp should always be half of that of the opposite-spin cusp.
Did I miss something?

Copy link
Contributor

Choose a reason for hiding this comment

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

The previous code was added in #2806 , and it references equation 9 of https://arxiv.org/abs/2003.06506 (and footnote 45 for the 1-D case). According to that formula, it looks like the same-spin cusp is half of the opposite spin cusp for d=3, but not for d=2.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yes, you are right. I must have made a mistake in my own derivation.
Thanks for catching this!

@prckent
Copy link
Contributor

prckent commented Oct 21, 2022

Fair point on not recoding the finite derivatives.

@prckent
Copy link
Contributor

prckent commented Oct 21, 2022

@Paul-St-Young You use <parameter name="ratio">deriv</parameter> or perhaps <parameter name="basic">yes</parameter>?

@@ -69,6 +70,12 @@ WaveFunctionTester::WaveFunctionTester(const ProjectData& project_data,

deltaR.resize(w.getTotalNum());
makeGaussRandom(deltaR);
Copy link
Contributor

Choose a reason for hiding this comment

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

Any reason not to make makeGaussRandom dimension aware and avoid the changes and logging on later lines?

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 reason is that I wanted the change to be localized in the WaveFunctionTester.

Eventually, what you suggested is nicer. Although, I would have to touch a lot of files to make this API change.
I wonder if it's worth making another PR for that.

Copy link
Contributor

Choose a reason for hiding this comment

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

Hi Paul - In that case I agree that making another PR would be a good way of splitting this up

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.

Just noticed: the tests are not being run in CI. Please can you rename them with the prefix deterministic- so that "ctest -L deterministic" will pick them up. I think the name is sufficient. Thanks.

edit: for consistency with the existing heg tests, how about "deterministic-heg2d-...". That avoids introducing a new abbreviation that people will have to learn the meaning of.

@prckent
Copy link
Contributor

prckent commented Oct 21, 2022

Test this please

@ye-luo
Copy link
Contributor

ye-luo commented Oct 22, 2022

@ye-luo @prckent The WavefunctionTester is by far the easiest route to test particle and parameter derivatives of a complicated wavefunction. When I doubt the correctness of a run, it's easy replace the <qmc method="vmc/dmc"> section with <qmc method="wftest"> for a quick test. I know that @rcclay and my self have made many uses of this feature, for example when debugging backflow. @markdewing also made the FiniteDifference class fairly clean.

WavefunctionTester is also my go-to when developing a new wavefunction component, for example when resurrecting 2D electron gas orbitals #4084. Since we have no other facility for generally checking the correctness of analytical derivatives (either for single-particle orbital, or for the many-body wavefunction), I don't see why we should remove this feature. Having to re-code finite-difference for every unit test is not productive IMO.

I am happy to put in a simple regression unit test, but I would like to keep the wftest to check for correctness.

You can use wftest for your development but the final test being added needs to be a unit test. wftest does random sampling which is not desired in unit tests. I don't see where you get the impression of the need of finite-difference in unit test. Unit test should check against stored values. See https://github.com/QMCPACK/qmcpack/blob/develop/src/QMCWaveFunctions/tests/test_polynomial_eeI_jastrow.cpp

@Paul-St-Young
Copy link
Contributor Author

@Paul-St-Young You use <parameter name="ratio">deriv</parameter> or perhaps <parameter name="basic">yes</parameter>?

Yes, I use only two features of the WavefunctionTester: the "basic" test for particle derivatives and move ratios and the "derivative" test for parameter derivatives.

@Paul-St-Young
Copy link
Contributor Author

Just noticed: the tests are not being run in CI. Please can you rename them with the prefix deterministic- so that "ctest -L deterministic" will pick them up. I think the name is sufficient. Thanks.

edit: for consistency with the existing heg tests, how about "deterministic-heg2d-...". That avoids introducing a new abbreviation that people will have to learn the meaning of.

Thanks fo advice! I have changed the name of the test.

@prckent
Copy link
Contributor

prckent commented Oct 24, 2022

Test this please

@ye-luo
Copy link
Contributor

ye-luo commented Oct 27, 2022

@prckent
Copy link
Contributor

prckent commented Oct 31, 2022

Paul - Are you able to reproduce any of these problems locally? e.g. Our gcc9 " no mpi and no omp " build is crashing with a segfault.

@Paul-St-Young
Copy link
Contributor Author

Yes, I was able to reproduce the crash.
The problem is only in the real build.

The input has 2 up electrons in 2 PW orbitals.
The 2 complex PW orbitals 0 and e^ikr are fine, but there are 3 real PW orbitals: 0, sin(kr), and cos(kr).
The extra memory was not allocated in DiracDeterminant.

How may I resize DiracDeterminant from within an SPOSet builder?

@ye-luo
Copy link
Contributor

ye-luo commented Nov 1, 2022

DiracDeterminant query the size from SPOSet. Make sure it is correct at the first place.

@Paul-St-Young
Copy link
Contributor Author

DiracDeterminant query the size from SPOSet. Make sure it is correct at the first place.

@ye-luo I'm confused by this.
In DiractDeterminant.cpp, resize is called with NumPtcls in both arguments.
If I make the following change, then my error goes way

-  resize(NumPtcls, NumPtcls);
+  const int norb = Phi->getOrbitalSetSize();
+  resize(NumPtcls, norb);

However, it breaks other tests.
Where does DiracDeterminant query the size from SPOSet to resize the determinant in the current version of the code?

@prckent
Copy link
Contributor

prckent commented Nov 1, 2022

I agree with Ye on this one. We should create diracdeterminants with the correct size to begin with. Cleaner and clearer. The old trope of changing the size of things after initial construction hurts the accessibility of the code and is not good for maintenance.

@PDoakORNL
Copy link
Contributor

I feel like this is in the end a multiple sources of truth issue. The fix is not resizing a determinant on the fly.

@Paul-St-Young
Copy link
Contributor Author

@prckent @PDoakORNL Yes, I agree with your view.
That that resize in the snippet above is from the constructor of DiracDeterminant, so it's not changing the size of things after initial construction.
The question is whether to initialize the determinants to (NumPtcls, NumPtcls) or (NumPtcls, norb).

@ye-luo
Copy link
Contributor

ye-luo commented Nov 1, 2022

The determinant matrices must be strictly square. (NumPtcls, NumPtcls). In the FreeOrbital input, you specified size 2 in the input. So the object should work for exactly or up to 2 orbitals. Based on what you mentioned that FreeOrbital creates 3 orbitals in the real case, then a truncation must be done at the construction of FreeOrbital object.

add some checks in FreeOrbitalBuilder
@Paul-St-Young
Copy link
Contributor Author

Thanks @ye-luo! Your explanation was very helpful.

I think the fundamental problem was that I set up a system that truly needs a complex wavefunction.
The 2 up 2 dn system at Gamma should have a complex single determinant wavefunction, because there is no -k PW to pair up with the +k PW to make sin and cos.

I added some checks in FreeOrbitalBuilder to abort in a situation as described above.

I also edited the example to contain 5 up 5 dn (closed-shell in 2D) and now things run through.

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.

Minor changes requested.

FiniteDifference fd(FiniteDifference::FiniteDiff_LowOrder);
//FiniteDifference fd(FiniteDifference::FiniteDiff_Richardson);
FiniteDifference fd(ndim, FiniteDifference::FiniteDiff_LowOrder);
//FiniteDifference fd(ndim, FiniteDifference::FiniteDiff_Richardson);
Copy link
Contributor

Choose a reason for hiding this comment

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

Remove this commented line

//RealType delta = 1.0;
//FiniteDifference fd(FiniteDifference::FiniteDiff_Richardson);
//FiniteDifference fd(ndim, FiniteDifference::FiniteDiff_Richardson);
Copy link
Contributor

Choose a reason for hiding this comment

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

Remove commented lines.

std::ostringstream msg;
msg << "norb = " << norb << " npw = " << npw;
msg << " cannot be ran in real PWs (sin, cos)" << std::endl;
msg << "either use complex build or change the size of SPO set" << std::endl;
Copy link
Contributor

Choose a reason for hiding this comment

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

Do you have recommendation to print about how to "change the size of SPO set"?

@@ -31,6 +34,14 @@ std::unique_ptr<SPOSet> FreeOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur)
app_log() << "twist cartesian = " << tvec << std::endl;
#else
npw = std::ceil((norb + 1.0) / 2);
Copy link
Contributor

Choose a reason for hiding this comment

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

Remove int npw in line 15 and put
int npw = here and line 31
Also please do a clang-formatting.

RealType qq = species(chargeInd, ia) * species(chargeInd, ib);
RealType red_mass = species(massInd, ia) * species(massInd, ib) / (species(massInd, ia) + species(massInd, ib));
#if OHMMS_DIM == 1
RealType dim_factor = 1.0 / (OHMMS_DIM + 1);
#else
RealType dim_factor = (ia == ib) ? 1.0 / (OHMMS_DIM + 1) : 1.0 / (OHMMS_DIM - 1);
RealType dim_factor = (ia == ib) ? 1.0 / (ndim + 1) : 1.0 / (ndim - 1);
Copy link
Contributor

Choose a reason for hiding this comment

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

Remove #if OHMMS_DIM == 1

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.

Please add random seeds in both determinsitic test xml.

@Paul-St-Young
Copy link
Contributor Author

@ye-luo I think I addressed all of your comments. Please see if this is satisfactory now.
The macos error seems to be an h5py install error, so not related to this PR

@ye-luo
Copy link
Contributor

ye-luo commented Nov 2, 2022

Test this please

@prckent
Copy link
Contributor

prckent commented Nov 2, 2022

I reran the failed macos test and it passed. Seems to have failed due to a transient infrastructure problem, nothing to do with this PR.

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
@prckent do you still have remaining concerns?

@ye-luo ye-luo enabled auto-merge November 2, 2022 14:59
@ye-luo ye-luo merged commit e6169b5 into QMCPACK:develop Nov 2, 2022
@prckent
Copy link
Contributor

prckent commented Nov 2, 2022

Thanks @Paul-St-Young !

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