strange results when changing lattice paramaters

[bxm...@gmail.com]

Hello,

We have been using the qmcpack recently for calculations (ground-state) of phosphorene under strain. We did calculations along a general line (i.e. we change all structural parameters, 2 lattice vectors (a, b) + 2 internal atomic coordinates (y, z), please see the respective projections below). We get a very weired behaviour of total energy comparing to the DFT calculated in quantum espresso used subsequently as inputs for the nodal surfaces in the DMC; while the structure changes smoothly in infinitesimal increments, we find a huge discontinuity in the total energy, see below. This is, of course totally unreasonable and there must be something going wrong in these calculations. We wonder if you could suggest what could go wrong. We append the basic inputs/outputs to/from these calculations. 4 points altogether done in a 5x5 supercell of phosphorene. We note that we have found a similar behavior also elsewhere, however, only if we change the  the lattice parameters (a and b). If we instead change only the internal atomic coordinates (y, z) everything looks fine. 

Your suggestion as to what could be the reason would be very helpful indeed.

Thank you. Best regards,

Jan Brndiar.

link for data http://147.213.112.235/pub/jano/qmcpack_forum.tgz   

[Paul R. C. Kent]

This looks like the kind of change that occurs when different states start to be occupied, for example a metal-insulator transition. It might be useful to check if the energies from every twist change similarly or if the energy from a single twist is dramatically different. I have not checked your tar file to see how many twists you have, but potentially you may need to do more twist averaging. Eyeing the DFT band structures can be informative. It is also possible that you have a failed/misbehaving/underconverged DMC run, so plotting the trace of energies of the individual DMC runs might point to the issue.

[Luke]

I heard that you are looking at the LR_DIM_CUTOFF parameter.  There is some description of this in the manual (see the simulationcell section).  Basically this controls how many terms are kept in the ewald summation.  Larger values are better and the default value of 15 tends to be rather conservative for most applications.

Luke

[Ye Luo]

In runDMC.xml,   <include href="P.xml"/> is wrong.

P.xml is a QE xml file. pw2qmcpack spits out P.ptcl.xml for particle and cell information.

However, pw2qmcpack can not tell your intention of a 5x5x1 supercell run with 1 1 1 0 0 0 supertwists or 1x1x1 supercell run with 5 5 1 0 0 0 supertwists.

Currently the P.ptcl.xml is written for the later case. Then in the wfs file, you can only use tilematrix="1 0 0 0 1 0 0 0 1".

If your intention is to run 5x5x1 cell, a P.ptcl.xml file needs to be prepared with proper cell and particle information.

Ye

在 2019年9月30日星期一 UTC-5上午10:06:12，bxm...@gmail.com写道：

[bxm...@gmail.com]

I'm posting what I have sent to Luke and Paul

Luke:

Hi

 Thanks for reply. I don't no if you have read also my post in QMCPACK google groups, so please see it https://groups.google.com/forum/#!topic/qmcpack/3Hkt5sWY_e4

I know that you are the first author of https://pubs.acs.org/doi/10.1021/acs.nanolett.5b03615 so It looks like you are the best person to help us to resolve our problems with interpretation of results of phosphorene obtained with QMCPACK. You have calculate DMC energies as function of the longer in-layer lattice constant and we want to do same somehow for whole parameters, so 4D instead of 1D. But we've got in troubles to fit our data even that our approach works fine for pure DFT data obtained with Espresso. Such a drop in energy in my post in groups is cut throw some line in 4D so changing everything (2 lattice, 2 for atom) in phosphorene. The strange is that we are able to successfully fit our data for fix a b, for energy dependence in y a z we got nice 2D paraboloid. That's I don't believe we have physical problem in inputs from Espresso as Paul R.C. Kent proposed or maybe wrongly compiled executable. So we have focused on changing of lattice constants a b or volume changes as a potential source of troubles. I have found that some part of calculation in QMCPACK for periodic system is evaluated with the help of LR Breakup where LR-dim = r_c* k_c and kcut is fixed to 60Pi times cube root of cell volume. So for test purpose we have calculated 100 VMC (10 steps with 0.7 timestep) blocks without warmup with different LR-dim and for large LR-dim we also enlarge kcut to 120Pi times cube root of cell volume (I have added for this test purpose code to change LR-dim and k_cut to QMCPACK because with inputs from espresso LR-dim is fixed in method ParticleSet* ParticleSetPool::createESParticleSet and k_cut is const in whole program) with same random number generator seed, to be able to directly compare  such VMC runs and we obtain for some points on above mentioned line cut massive changes in VMC block energies and for others it were nicely rapidly converged in LR-dim. And it was perfectly correlated with the drop in energy on line cut, so points before drop were converged, points after not. Did you observe such odd behavior (discontinuity in energy) for your calculations because you also have changed volume or you do not test such dependence. And the main think is that after convergence in the VMC energies in LR-dim we have suppose to obtain smooth dependence of energy on line cut, we again finished with different value discontinuity. We have used 3.6.0 version but also 3.8.0 version produce same convergence rates. So I have attached all inputs/outputs and also convergence VMC figures. Could you look into our input xml maybe I have some stupid error or misunderstand of manual which produce such strange results.

PS: here are just VMC (runVMC.xml). Our DMC runs in/out (runDMC.xml) are linked in google groups post, please check

to all:

My apologies for my insistence. I am struggling  to make sense of my QMCPACK energies. As you know, I am calculating strained phosphorene. I see huge discontinuities in total energies when I change the lattice parameters (and ionic positions) by very small amounts. What I empirically observe is that along with ETOT also the ion-ion energy and all the other terms calculated on the periodic lattice exhibit discontinuous behavior. It looks like they all are offset by a different but fixed amount. Noticing this, I have tried to enlarge LR-DIM in the effort to bring the energy to convergence. Energies of both structures appear to converge (at very different LR-DIM) but to two very different numbers. This does not make any sense to me. 

Could you, please, have a quick look at the calculation (I sent more in my previous email to Luke) and could you tell me what's going wrong? This is very urgent for me as we have already burnt a fair amount of corehours trying to map out ground-state energies as a function of the lattice parameters.

Your assistance is highly appreciated as I have already exhausted all I can do on my own.

Paul:

Dear Paul,

Thank you for your time, I do appreciate.

Here are a quick answers to your questions:

Q. Do the VMC energies show the same issue?

A. Yes, please see figure VMC.

Q. Do the no jastrow VMC energies track the DFT results as expected? (would catch input errors/bugs)

A. No, the runs without Jastrow are qualitatively the same as those with the Jastrow and different from the DFT results, please see figure VMC_no_jastrow.

We have ourselves looked more in detail on what goes wrong and we find that the energy components which show the jumps are:

- local potential

- ion-ion

- local ECP,

i.e. there is a discontinuity in all energies where periodic arrays of ions enters. I note that I use MPC on electrons.

As already said, I believe I have exhausted all I can do on my own. Your assistance is most appreciated.

I'm also attaching images:

Answer to Ye Luo:

  Actually I have realized that QMCPACK ignore P.xml and read everything from h5 file from PW where I have also cell coordinates and orbitals tilling and also set information about periodicity to 3D and cutoffs to constant 15.0 see answer to Luke

[bxm...@gmail.com]

Hello again,

I've been pursuing the problem I recently encountered further and found that the mentioned problem manifests itself already in the ion-ion energy (I do believe this is what is in the QMCPACK logfile called VTOT). I give a few examples in comparison with VASP, as an example, using exactly the same cell geometry and ionic positions:

size	nel	vasp [eV]	vasp per cell [eV]	qmcpack [eV]	qmcpack per cell [eV]	qmcpack-vasp [eV]
1	1	1903.65017598	1903.65017598	1903.33834497932	1903.33834497932	-0.311831000679376
2	4	7614.60044955	1903.6501123875	7614.50672316124	1903.62668079031	-0.093726388760842
3	9	17132.85096064	1903.65010673778	17144.1202585577	1904.90225095085	11.2692979176645
4	16	30458.40133538	1903.65008346125	30458.1200712728	1903.63250445455	-0.281264107215975
5	25	47591.25103742	1903.6500414968	47592.9109545248	1903.71643818099	1.65991710482922

The table above demonstrates that:

1) the QMCPACK Vtot energies are diiferent from those of VASP

2) more importantly, unlike the VASP energies, they are NOT constant per cell.

I presume that that also indicates that the other energy terms (which use the same long-range scheme) may be similarly biased which is indeed what I observe.

For your reference I enclose the 1x1 (first row in above table) geometry of the system used in the above test and also attached all i/o from qmcpack from this tests.

ATOMIC_POSITIONS angstroms
P       0       0.371128        2.12622
P       0       4.078822        0
P       1.66918 1.853847        2.12622
P       1.66918 2.596103        0

CELL_PARAMETERS angstroms
3.33836 0       0
0       4.44995 0
0       0       20

I wish to say again that at this point I feel I can't get any further on my own. Hence, I do rely on your kind support.

[Jaron Krogel]

Hi,

The default Ewald implementation in QMCPACK (the "optimized breakup") is both accurate and efficient for most condensed 3D systems.  However, this implementation has a bias in the energy for quasi-2D systems which grows as the vacuum region is increased.

An alternative implementation with controllable errors (based on standard Gaussian/error function Ewald formulae) has been made available to users in the development branch of QMCPACK ahead of the v3.9.0 release, albeit at the price of summing over additional k-vectors which adds to the computational cost.

This method is selected by providing the following parameter in the <simulationcell> element in QMCPACK's input file:

<parameter name="LR_handler"   >  ewald  </parameter>

The accuracy of the Ewald sum is increased by increasing the "LR_dim_cutoff" parameter (in <simulationcell>), which defaults to a value of 15.0 (unitless):

<parameter name="LR_dim_cutoff">  15     </parameter>

A high accuracy check of the ion-ion total energy is now performed during the initialization phase of every QMCPACK run.  If the error in the per-atom ion-ion energy exceeds a user-provided tolerance, the run is halted.  The tolerance is provided in <simulationcell> via the parameter "LR_tol", which is in units of Hatree/atom and defaults to 3e-4 Ha/atom:

<parameter name="LR_tol"       >  3e-4   </parameter>

Prior to any production runs, the accuracy of the Ewald sum can be checked on a laptop or workstation in the following way.  For concreteness, an example for 1x1 graphene is provided (see attached vmc.in.xml, no other files are needed to run "qmcpack vmc.in.xml").

First, make a QMCPACK input file for the system of interest.   Comment out any sections requiring external files such as <sposet_builder>, <determinantset>, or the pseudopotential <pairpot> and set the VMC section to take only a single step as in the attached vmc.in.xml file.

Next, decide on a target accuracy.  For quasi-2D exfoliation energies, one might wish to target 1 meV/A^2 or perhaps 0.1 meV/A^2.  For the graphene example, these targets correspond to about 1e-4 Ha/atom and 1e-5 Ha/atom, respectively, unless I've made a trivial mistake.

Search for the most compact (efficient) sum that meets your accuracy requirements.  Start by setting the Ewald-related values in the <simulationcell> element to the following (default LR_dim_cutoff):

<simulationcell>

   ... other parameters are excluded for brevity ...

   <parameter name="LR_handler"   >  ewald  </parameter>

   <parameter name="LR_dim_cutoff">  15     </parameter>

   <parameter name="LR_tol"       >  1e-12  </parameter>

</simulationcell>

Run QMCPACK once and observe the level of error in the ion-ion energy per atom.  For reference, also pay attention to the number of k-points in the sum.  Abbreviated output from running QMCPACK with the attached vmc.in.xml for graphene are as follows:

  Creating CoulombHandler with the 3D Ewald Breakup.

  Sigma=1.35939

   PBCAA self-interaction term -6.13562

   PBCAA total constant -6.20817

  Maximum K shell 435

  Number of k vectors 5338

  Fixed Coulomb potential for e

    e-e Madelung Const. =0.665507

    Vtot     =0

  QMCHamiltonian::addOperator ElecElec to H, physical Hamiltonian

QMCHamiltonian::addOperatorType added type coulomb named ElecElec

  Clone CoulombHandler.

   PBCAA self-interaction term -24.5425

   PBCAA total constant -24.615

Checking ion-ion Ewald energy against reference...

Error in ion-ion Ewald energy exceeds 1e-12 Ha/atom tolerance.

  Reference ion-ion energy: 51.710928350571

  QMCPACK   ion-ion energy: 51.708422776574

            ion-ion diff  : -0.0025055739966575

            diff/atom     : -0.0012527869983288

            tolerance     : 1e-12

Increase LR_dim_cutoff and rerun QMCPACK until the error per atom ("diff/atom") is less than or equal to your desired target accuracy.  For the graphene example, the behavior of the error with increasing cutoff is as follows:

    cutoff   kpoints  error

    15        5338    1.25e-3  Ha/atom

    20       12528    1.07e-4  Ha/atom

    25       24140    9.80e-6  Ha/atom

We see that the accuracy thresholds of 1 meV/A^2 (1e-4 Ha/atom) and 0.1 meV/A^2 (1e-5 Ha/atom) are met with LR_dim_cutoffs of about 20.0 and 25.0, respectively.  As noted above , the number of k-points required in the long-ranged portion of the Ewald sum grows rapidly, increasing the cost of the run.   The cost is linear in the number of k-points.

Finally, set LR_dim_cutoff to the converged value (e.g. 25.0 above) and LR_tol to the desired tolerance (e.g. 1e-5) and proceed to production runs.

Best,

Jaron

NB: With future algorithmic improvements, the number of k-points needed to reach the target accuracy might be reduced by as much as 5x relative to the current implementation.  If the speed of your run is much degraded relative to the default LR_dim_cutoff of 15.0, please reach out to the developers to communicate your need for the improved algorithm.

--
You received this message because you are subscribed to the Google Groups "qmcpack" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qmcpack+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/qmcpack/1d1a01fc-5ac2-4d82-8eea-bfd3a3fbee67%40googlegroups.com.