Calculation of forces in QMC

[Vladimir_Konjkov]

It seems that I forgot to change the valence charge of C atom from 6 to 4 (my script molden2qmc does not perform it), but despite this the presence of c_pp.data and h_pp.data files initiated casino to perform calculations with a pseudopotential.
Energy obtained agrees with that from ORCA with same pseudo. Whether the valence charge has any effect on QMC calculation?

[Mike Towler]

Hi Vladimir,

The way that CASINO detects whether atoms defined by the external code are pseudoatoms or not is more difficult than you might think (because it relies on the external code and/or the converter behaving nicely, which can't be assumed).

If I can quote the DIARY entry 2.13.379 from one of my changes last summer:

Code: Select all

---[v2.13.379]---
* Revamped method for distinguishing between all-electron and pseudoatoms
  (which is harder than it sounds).
  -- Mike Towler, 2014-06-19

  Seemingly simple task of detecting whether an atom is a pseudoatom was not
  always done correctly (only to the level of e.g. error traps being
  incorrectly triggered, or wasting a bit of time/memory, rather than anything
  serious - like getting the wrong answer). Also CASINO requires explicit
  setting of input flags if you want to do anything 'weird' like mixing ae
  and pp atoms, or using pps with numerical orbitals etc. (which should not
  be necessary)

  Exact behaviour previously dependent on confusing mishmash of internal
  keywords HAVE_AE, ALLOW_AE_PPOTS, CUSP_CORRECTION, HAVE_AE_DEFINED, the first
  three of which can be redefined in input but hardly ever are, and if they
  are, can have the 'wrong' value if the user doesn't specifically bother
  changing them.

  Problems with current implementation:

  - Definition of pseudoatom should come from external code, but current
  implementation makes it depend on 'more likely defaults' and hence on
  stuff in CASINO input file like the CUSP_CORRECTION flag.

  - Internal assumption in some places (e.g. definition of HAVE_AE keyword)
  that all external Gaussian codes flag pseudos with atno>200 (true
  historically only for CRYSTAL and Badinski GAMESS but not GAUSSIAN, for some
  reason. Recently implemented interfaces MOLDEN and Defusco internal GAMESS
  also do not do this). Can lead to incorrect errstop if e.g. pseudopotential
  file exists AND gwfn.data produced by code which doesn't obey +200
  convention AND neither HAVE_AE, CUSP_CORRECTION, or ALLOW_AE_PPOTS is
  defined in input (hence taking on their defaults)

  Changes made :

  (1) HAVE_AE keyword flagged as redundant and is now only used internally.
  Default was previously set as 'the most likely value' for a given basis set
  i.e. F for plane-wave/blip, T for numerical, 'cusp_correction' for Gaussian
  etc., which is a bit stupid. Now no default is assumed except for Slater
  functions (since the only Slater code -- ADF -- does not support
  pseudopotentials), and special wave function types like nonint_he,h2,h3plus
  etc.. since they don't use pseudos by construction.  For other
  ATOM_BASIS_TYPEs, see (3).

  (2) ALLOW_AE_PPOTS keyword is flagged as redundant and removed internally.
  Only reason for its existence is to prevent the error where the user
  accidentally omits one or more pseudopotential files in systems with more
  than one kind of atom, but (a) for Gaussians there are other essentially
  infallible methods of checking that, and (b) on the very rare occasions you
  want to do this with some other basis, well, sometimes you just have to let
  the user screw things up.

  (3) ATOM_BASIS_TYPEs where the user can in principle make a choice of whether
  to have pseudopotentials or not now behave as follows:

  - Gaussians : in addition to checking whether the x_pp.data files actually
  exist, it uses two other pieces of information to check whether they ought to
  exist: (a) if atomic_number has 200 added to it, this definitely flags a
  pseudoatom (but < 200 doesn't imply anything, as some external codes don't
  support this convention), and (b) if valence charges defined in gwfn.data /=
  atomic number that definitely flags a pseudoatom (but if equal, that doesn't
  necessarily mean all-electron, e.g. H pseudoatom retaining its single
  electron. Code will now stop if it detects any blatant contradictions. Note
  that the +200 convention is not essential; if present it offers complete
  protection against the missing pseudo file error. The only error that would
  not also be detected by the valence_charge check is the missing H pseudo
  where no core electrons are explicitly removed.

  - plane waves/blips : default HAVE_AE=F so that it expects to find a full set
  of x_pp.data files (which it should in almost all cases). If it doesn't, it
  searches for the gzipped and bzipped equivalents x_pp.data.gz and
  x_pp.data.bz2, stopping if it finds them. If it doesn't find the compressed
  files either then it assumes all-electron and carries on as normal. [EDIT:
  and from 2.13.381 runqmc will issue a warning about this.]

  - numerical/dimer: orbitals tabulated on grids are normally all-electron
  (an exception is examples/atom/nickel_pp_n). No explicit checking is done; if
  CASINO finds a pseudopotential file, it will attempt to use it. If you don't
  want it to do that, don't give it a pseudopotential. This is an improvement
  on the previous behaviour, where the default HAVE_AE was T, and the user had
  to remember to set HAVE_AE=F explicitly in input if he wanted to use a
  pseudopotential.

  (4) Made sure gwfn.data specification in the manual specifically states the
  +200 convention (previously it was only implied through an example).

  (5) Had the MOLDEN converter changed to follow the +200 invention (thanks to
  Mike Deible)

  (6) Had the forthcoming internal GAMESS converter adopt the +200 convention
  (thanks to Albert Defusco). The old Badinski one distributed with CASINO
  already does.

  (7) Broke link between the value of CUSP_CORRECTION or USE_GPCC in input and
  the definition of a pseudoatom and HAVE_AE (e.g. since CUSP_CORRECTION just
  might accidentally be T in input because the user hasn't bothered to change
  it - should just be ignored if all atoms are pseudoatoms). Made sure that
  having CUSP_CORRECTION=T in input does not lead to running of Gaussian
  cusp_setup in the case where every atom is a pseudoatom, which takes time
  and whose data occupies memory and is never used (this could happen in
  certain weird circumstances before).

  TODO: The GAUSSIANXX converter still does not implement the +200 convention,
  as it's not immediately obvious how to extract that information from the
  hideous output (no offence, GAUSSIAN fans). Feel free to persuade it to do
  so..

The important points being:

(1) the main thing used by CASINO is the existence or not of the x_pp.data file. If it exists, that means it's a pseudoatom. If it doesn't exist, then it's all-electron. However, the user could simply have got that wrong (forgetting the pseudo file, for instance) so CASINO uses other checks (valence_change/=atomic_charge, or for the existence of the +200 flag) so that it can stop or warn about inconsistencies.

(2) in your MOLDEN converter you need to add +200 to the atomic number of atoms that are represented by a pseudopotential (the only guaranteed way
to always get it right) - sorry should have mentioned that before. Not all external codes/converters use this convention but when you're writing a new converter it's obviously best to do so.

(3) valence charge equalling the atomic charge is not a guarantee of all-electron (since e.g. for H these are often the same..)

The answer to your question is that the valence charge is used explicitly in a CASINO total energy calculation to compute the nuclear-repulsion term, but that this is obviously zero for a C atom - so you're not seeing a difference. It isn't used explicitly in the e-n term because the valence charge is implicit in the pseudopotential data. However it will be used to define the local cutoff radius, and having the wrong value of that might lead to small differences in the local electron-ion term?

Mike

PS: I haven't forgotten I said I would take a look at the post-calculation force analysis.. I'll try to get onto that tomorrow. Sorry for the delay..

[Mike Towler]

Hi Vladimir,

So yes, the reblock utility is used to calculate the forces from the info in the vmc.hist file.

I've just checked it over, and it seems to work. However, I have cleaned up the output a bit (the forces guy Alex had a very interesting notion of the English language..) and I've fixed a small error which occurs when it's asking you how you want to reblock the force error bars (you probably wouldn't have noticed this). The modified reblock utility is now in the latest CASINO current beta distribution (2.13.501).

Fixing plot_hist to plot forces is a little more involved than I thought (the accumulation of 22 separate quantities is flagged when you set forces = T in input). You don't need it, and I'll make sure it's done for the forthcoming release (there are more important things going on right now)..

Cheers,
Mike

[Vladimir_Konjkov]

It seems to me that the utility reblock does not change REBLOCKED ERROR BARS when increasing the number of steps.
For example, two calculations of forces in cubane molecule with basis aug-cc-pVDZ-CDF using "DF AREP Trail Needs ECP pseudopotential"
for 10,000 and 100,000 steps.

P.S. I did not attach files vmc.hist because they have a large volume (50mb и 500mb).

[Vladimir_Konjkov]

I regret the error was in my script that performs an additional averaging based on cubane molecule symmetry.

[Vladimir_Konjkov]

I discovered that local pseudopotentials without core electrons have been developed for H-Be, they also contain only local channel.
(https://pubs.aip.org/aip/jcp/article/15 ... tials-from).
They are available on the website https://pseudopotentiallibrary.org
They will have less problems when calculating VMC and DMC forces.