search

QMCPACK / qmcpack

Main repository for QMCPACK, an open-source production level many-body ab initio Quantum Monte Carlo code for computing the electronic structure of atoms, molecules, and solids with full performance portable GPU support
http://www.qmcpack.org
Other
295 stars 138 forks source link

Fatal Error with Twists for graphene example in nexus. #2411

Closed kryczko closed 3 years ago

kryczko commented 4 years ago

Describe the bug When trying to use a 2x2x1 k-point grid with the graphene example in nexus I receive the following error during the qmcpack optimization: 

Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2

To Reproduce Run script with PWSCF installed with qmcpack patch: 

#! /usr/bin/env python3

from nexus import settings,job,run_project
from nexus import generate_physical_system
from nexus import generate_pwscf
from nexus import generate_pw2qmcpack
from nexus import generate_qmcpack
from nexus import loop,linear,vmc,dmc

import sys

n_tasks = sys.argv[1]

# general settings for nexus
settings(
    pseudo_dir    = '../../pseudopotentials',# directory with all pseudopotentials
    sleep         = 3,                       # check on runs every 'sleep' seconds
    generate_only = 0,                       # only make input files
    status_only   = 0,                       # only show status of runs
    machine       = 'ws40',                  # local machine is 16 core workstation
    )

# generate the graphene physical system
graphene = generate_physical_system(
    lattice   = 'hexagonal',      # hexagonal cell shape
    cell      = 'primitive',      # primitive cell
    centering = 'P',              # primitive basis centering
    constants = (2.462,10.0),     # a,c constants
    units     = 'A',              # in Angstrom
    atoms     = ('C','C'),        # C primitive atoms
    basis     = [[ 0  , 0  , 0],  # basis vectors
                 [2./3,1./3, 0]],
    tiling    = (5,5,1),          # tiling of primitive cell
    kgrid     = (2,2,1),          # Monkhorst-Pack grid
    kshift    = (.5,.5,.5),       # and shift
    C         = 4                 # C has 4 valence electrons
    ) 

# scf run produces charge density
scf = generate_pwscf(
    # nexus inputs
    identifier   = 'scf',           # identifier/file prefix
    path         = 'graphene/scf',  # directory for scf run
    job          = job(cores=n_tasks),   # run on 16 cores
    pseudos      = ['C.BFD.upf'],   # pwscf PP file
    system       = graphene,        # run graphene
    # input format selector
    input_type   = 'scf',           # scf, nscf, relax, or generic
    # pwscf input parameters
    input_dft    = 'lda',           # dft functional
    ecut         =  200 ,           # planewave energy cutoff (Ry)
    conv_thr     =  1e-6,           # scf convergence threshold (Ry)
    mixing_beta  =    .7,           # charge mixing factor
    kgrid        = (2,2,1),         # MP grid of primitive cell
    kshift       = (.5,.5,.5),         #  to converge charge density
    wf_collect   = False,           # don't collect orbitals
    use_folded   = True,            # use primitive rep of graphene
    )

# nscf run to produce orbitals for jastrow optimization
nscf_opt = generate_pwscf(
    # nexus inputs
    identifier   = 'nscf',          # identifier/file prefix      
    path         = 'graphene/nscf_opt', # directory for nscf run       
    job          = job(cores=n_tasks),   # run on 16 cores             
    pseudos      = ['C.BFD.upf'],   # pwscf PP file               
    system       = graphene,        # run graphene                
    # input format selector                                      
    input_type   = 'nscf',          # scf, nscf, relax, or generic
    # pwscf input parameters
    input_dft    = 'lda',           # dft functional
    ecut         =  200 ,           # planewave energy cutoff (Ry)
    conv_thr     =  1e-6,           # scf convergence threshold (Ry)
    mixing_beta  =    .7,           # charge mixing factor
    nosym        = True,            # don't symmetrize k-points
    use_folded   = True,            # use primitive rep of graphene
    wf_collect   = True,            # write out orbitals
    kgrid        = (2,2,1),         # single k-point for opt
    kshift       = (0.5,0.5,0.5),         # gamma point
    # workflow dependencies
    dependencies = (scf,'charge_density'),
    )

# orbital conversion job for jastrow optimization
p2q_opt = generate_pw2qmcpack(
    # nexus inputs
    identifier   = 'p2q',
    path         = 'graphene/nscf_opt',
    job          = job(cores=1),
    # pw2qmcpack input parameters
    write_psir   = False,
    # workflow dependencies
    dependencies = (nscf_opt,'orbitals'),
    )

# Jastrow optimization
opt = generate_qmcpack(
    # nexus inputs
    identifier   = 'opt',           # identifier/file prefix
    path         = 'graphene/opt',  # directory for opt run
    job          = job(cores=40,app='qmcpack'),
    pseudos      = ['C.BFD.xml'],   # qmcpack PP file
    system       = graphene,        # run graphene
    # input format selector   
    input_type   = 'basic',
    # qmcpack input parameters
    corrections  = [], 
    jastrows     = [('J1','bspline',8),   # 1 body bspline jastrow
                    ('J2','bspline',8)],  # 2 body bspline jastrow
    calculations = [
        loop(max = 6,                        # No. of loop iterations
             qmc = linear(                   # linearized optimization method
                energy               =  0.0, # cost function
                unreweightedvariance =  1.0, #   is all unreweighted variance
                reweightedvariance   =  0.0, #   no energy or r.w. var. 
                timestep             =  0.5, # vmc timestep (1/Ha)
                warmupsteps          =  100, # MC steps before data collected 
                samples              = 16000,# samples used for cost function 
                stepsbetweensamples  =   10, # steps between uncorr. samples
                blocks               =   10, # ignore this  
                minwalkers           =   0.1,#  and this
                bigchange            =  15.0,#  and this
                alloweddifference    =  1e-4 #  and this, for now
                )
             )        
        ],
    # workflow dependencies
    dependencies = (p2q_opt,'orbitals'),
    )

# nscf run to produce orbitals for final dmc
nscf = generate_pwscf(
    # nexus inputs
    identifier   = 'nscf',          # identifier/file prefix      
    path         = 'graphene/nscf', # directory for nscf run       
    job          = job(cores=n_tasks),   # run on 16 cores             
    pseudos      = ['C.BFD.upf'],   # pwscf PP file               
    system       = graphene,        # run graphene                
    # input format selector                                      
    input_type   = 'nscf',          # scf, nscf, relax, or generic
    # pwscf input parameters
    input_dft    = 'lda',           # dft functional
    ecut         =  200 ,           # planewave energy cutoff (Ry)
    conv_thr     =  1e-6,           # scf convergence threshold (Ry)
    mixing_beta  =    .7,           # charge mixing factor
    nosym        = True,            # don't symmetrize k-points
    use_folded   = True,            # use primitive rep of graphene
    wf_collect   = True,            # write out orbitals
    # workflow dependencies
    dependencies = (scf,'charge_density'),
    )

# orbital conversion job for final dmc
p2q = generate_pw2qmcpack(
    # nexus inputs
    identifier   = 'p2q',
    path         = 'graphene/nscf',
    job          = job(cores=1),
    # pw2qmcpack input parameters
    write_psir   = False,
    # workflow dependencies
    dependencies = (nscf,'orbitals'),
    )

# final dmc run
qmc = generate_qmcpack( 
    # nexus inputs
    identifier   = 'qmc',           # identifier/file prefix       
    path         = 'graphene/qmc',  # directory for dmc run       
    job          = job(cores=n_tasks,app='qmcpack'),
    pseudos      = ['C.BFD.xml'],   # qmcpack PP file
    system       = graphene,        # run graphene
    # input format selector                                      
    input_type   = 'basic',
    # qmcpack input parameters
    corrections  = [],              # no finite size corrections
    jastrows     = [],              # overwritten from opt
    calculations = [                # qmcpack input parameters for qmc
        vmc(                        # vmc parameters 
            timestep      = 0.5,    # vmc timestep (1/Ha)
            warmupsteps   = 100,    # No. of MC steps before data is collected
            blocks        = 200,    # No. of data blocks recorded in scalar.dat
            steps         =  10,    # No. of steps per block
            substeps      =   3,    # MC steps taken w/o computing E_local
            samplesperthread = 40   # No. of dmc walkers per thread
            ),                      
        dmc(                        # dmc parameters
            timestep      = 0.01,   # dmc timestep (1/Ha)
            warmupsteps   =  50,    # No. of MC steps before data is collected
            blocks        = 400,    # No. of data blocks recorded in scalar.dat
            steps         =   5,    # No. of steps per block
            nonlocalmoves = True    # use Casula's T-moves
            ),                      #  (retains variational principle for NLPP's)
        ],
    # workflow dependencies
    dependencies = [(p2q,'orbitals'),
                    (opt,'jastrow')],
    )

# nexus monitors all runs
run_project()

# print out the total energy
performed_runs = not settings.generate_only and not settings.status_only
if performed_runs:
    # get the qmcpack analyzer object
    # it contains all of the statistically analyzed data from the run
    qa = qmc.load_analyzer_image()
    # get the local energy from dmc.dat
    le = qa.dmc[1].dmc.LocalEnergy  # dmc series 1, dmc.dat, local energy
    #  print the total energy for the 8 atom system
    print('The DMC ground state energy for graphene is:')
    print('    {0} +/- {1} Ha'.format(le.mean,le.error))
#end if

Expected behavior See Nexus documentation for expected output.

System: Niagara system on Compute Canada. Modules used are intel, openmpi, mkl, and hdf5-mpi.

Additional context I am trying to run the tutorial with twist averaging.

jtkrogel commented 4 years ago

Would you mind posting the QMCPACK input xml and log output for the failed case?

kryczko commented 4 years ago

Actually, I tried running the original example for graphene given with qmcpack and got another twist error. The input is: 

<?xml version="1.0"?>
<simulation>
   <project id="opt" series="0">
      <application name="qmcpack" role="molecu" class="serial" version="1.0"/>
   </project>
   <qmcsystem>
      <simulationcell>
         <parameter name="lattice" units="bohr">
                  9.30501148        0.00000000        0.00000000
                 -4.65250574        8.05837632        0.00000000
                  0.00000000        0.00000000       18.89726133
         </parameter>
         <parameter name="bconds">
            p p p
         </parameter>
         <parameter name="LR_dim_cutoff"       >    15                 </parameter>
      </simulationcell>
      <particleset name="e" random="yes">
         <group name="u" size="16" mass="1.0">
            <parameter name="charge"              >    -1                    </parameter>
            <parameter name="mass"                >    1.0                   </parameter>
         </group>
         <group name="d" size="16" mass="1.0">
            <parameter name="charge"              >    -1                    </parameter>
            <parameter name="mass"                >    1.0                   </parameter>
         </group>
      </particleset>
      <particleset name="ion0">
         <group name="C" size="8" mass="21894.71359057295">
            <parameter name="charge"              >    4                     </parameter>
            <parameter name="valence"             >    4                     </parameter>
            <parameter name="atomicnumber"        >    6                     </parameter>
            <parameter name="mass"                >    21894.71359057295            </parameter
>
            <attrib name="position" datatype="posArray" condition="0">
                     0.00000000        0.00000000        0.00000000
                     2.32625287        1.34306272        0.00000000
                     4.65250574        0.00000000        0.00000000
                     6.97875861        1.34306272        0.00000000
                    -2.32625287        4.02918816        0.00000000
                    -0.00000000        5.37225088        0.00000000
                     2.32625287        4.02918816        0.00000000
                     4.65250574        5.37225088        0.00000000
            </attrib>
         </group>
      </particleset>
      <wavefunction name="psi0" target="e">
         <sposet_builder type="bspline" href="../nscf_opt/pwscf_output/pwscf.pwscf.h5" tilematr
ix="2 0 0 0 2 0 0 0 1" twistnum="0" source="ion0" version="0.10" meshfactor="1.0" precision="fl
oat" truncate="no">
            <sposet type="bspline" name="spo_ud" size="16" spindataset="0"/>
         </sposet_builder>
         <determinantset>
            <slaterdeterminant>
               <determinant id="updet" group="u" sposet="spo_ud" size="16"/>
               <determinant id="downdet" group="d" sposet="spo_ud" size="16"/>
            </slaterdeterminant>
         </determinantset>
         <jastrow type="One-Body" name="J1" function="bspline" source="ion0" print="yes">
            <correlation elementType="C" size="8" rcut="4.65250573818748" cusp="0.0">
               <coefficients id="eC" type="Array">                  
0 0 0 0 0 0 0 0
               </coefficients>
            </correlation>
         </jastrow>
         <jastrow type="Two-Body" name="J2" function="bspline" print="yes">
            <correlation speciesA="u" speciesB="u" size="8" rcut="4.65250573818748">
               <coefficients id="uu" type="Array">                  
0.48187490362234775 0.3894836478661487 0.2798979091850247 0.1788222270584585 
0.10155095246178476 0.05124984473968947 0.02297906259440177 0.009151077742247856 
               </coefficients>
            </correlation>
            <correlation speciesA="u" speciesB="d" size="8" rcut="4.65250573818748">
               <coefficients id="ud" type="Array">                  
0.6800217727257846 0.5185841205635617 0.3542238636562356 0.21663372481070922 
0.11854809110839902 0.05800123631539221 0.025348648741029732 0.009885944877629274 
               </coefficients>
            </correlation>
         </jastrow>
      </wavefunction>
      <hamiltonian name="h0" type="generic" target="e">
         <pairpot type="coulomb" name="ElecElec" source="e" target="e"/>
         <pairpot type="coulomb" name="IonIon" source="ion0" target="ion0"/>
         <pairpot type="pseudo" name="PseudoPot" source="ion0" wavefunction="psi0" format="xml"
>
            <pseudo elementType="C" href="C.BFD.xml"/>
         </pairpot>
      </hamiltonian>
   </qmcsystem>
   <loop max="6">
      <qmc method="linear" move="pbyp" checkpoint="-1">
         <cost name="energy"              >    0.0                </cost>
         <cost name="unreweightedvariance">    1.0                </cost>
         <cost name="reweightedvariance"  >    0.0                </cost>
         <parameter name="warmupSteps"         >    100                </parameter>
         <parameter name="blocks"              >    10                 </parameter>
         <parameter name="timestep"            >    0.5                </parameter>
         <parameter name="stepsbetweensamples" >    10                 </parameter>
         <parameter name="samples"             >    16000              </parameter>
         <parameter name="minwalkers"          >    0.1                </parameter>
         <parameter name="nonlocalpp"          >    yes                </parameter>
         <parameter name="use_nonlocalpp_deriv">    yes                </parameter>
         <parameter name="alloweddifference"   >    0.0001             </parameter>
         <parameter name="bigchange"           >    15.0               </parameter>
      </qmc>
   </loop>
</simulation>

Output is:

  Input file(s): opt.in.xml 

=====================================================
                    QMCPACK 3.9.9

       (c) Copyright 2003-  QMCPACK developers

                    Please cite:
 J. Kim et al. J. Phys. Cond. Mat. 30 195901 (2018)
      https://doi.org/10.1088/1361-648X/aab9c3

  Git branch: develop
  Last git commit: 5a38ba78887b391917214141b6d4d4510b212703-dirty
  Last git commit date: Thu Apr 23 10:25:12 2020 -0500
  Last git commit subject: Merge pull request #2399 from mcbennet/nxs_xmlreader
=====================================================
  Global options 

  Total number of MPI ranks = 16
  Number of MPI groups      = 1
  MPI group ID              = 0
  Number of ranks in group  = 16
  MPI ranks per node        = 16
  OMP 1st level threads     = 1
  OMP nested threading disabled or only 1 thread on the 2nd level

  Precision used in this calculation, see definitions in the manual:
  Base precision      = double
  Full precision      = double

  Structure-of-arrays (SoA) optimization enabled

  Input XML = opt.in.xml

  Project = opt
  date    = 2020-04-27 09:03:10 EDT
  host    = nia0128.scinet.local

 Random Number
 -------------
  Offset for the random number seeds based on time: 14

  Range of prime numbers to use as seeds over processors and threads = 53-131

 Lattice
 -------
  Simulation cell radius   = 4.029188 bohr
  Wigner-Seitz cell radius = 4.652506 bohr

  Lattice (bohr):      9.3050114800      0.0000000000      0.0000000000
                      -4.6525057400      8.0583763200      0.0000000000
                       0.0000000000      0.0000000000     18.8972613300

  Boundary Conditions:  p  p  p 

  Volume (bohr^3) = 1416.9787162998

  Reciprocal vectors without 2*pi.
    g_1 =       0.1074689700      0.0620472388     -0.0000000000
    g_2 =       0.0000000000      0.1240944776      0.0000000000
    g_3 =       0.0000000000     -0.0000000000      0.0529177209

  Metric tensor in real-space.
    h_1 = 86.5832386429 -43.2916193215 0.0000000000 
    h_2 = -43.2916193215 86.5832385755 0.0000000000 
    h_3 = 0.0000000000 0.0000000000 357.1064857743 

  Metric tensor in g-space.
    h_1 = 0.6079454982 0.3039727493 0.0000000000 
    h_2 = 0.3039727493 0.6079454987 0.0000000000 
    h_3 = 0.0000000000 0.0000000000 0.1105508278 
 Particle Set 
 ------------
  Name: e
  Initializing the lattice by the global supercell
  All the species have the same mass 1.0000000000
  Long-range breakup parameters:
    rc*kc = 15.0000000000; rc = 4.0291881600; kc = 3.7228343290; tol = 0.0003000000

  Creating Structure Factor for periodic systems 3.7228343290
  KContainer initialised with cutoff 3.7228343290
   # of K-shell  = 99
   # of K points = 1270

  Particle set size: 32

 Particle Set 
 ------------
  Name: ion0
  Initializing the lattice by the global supercell
  All the species have the same mass 21894.7135905729
  Long-range breakup parameters:
    rc*kc = 15.0000000000; rc = 4.0291881600; kc = 3.7228343290; tol = 0.0003000000

  Creating Structure Factor for periodic systems 3.7228343290
  KContainer initialised with cutoff 3.7228343290
   # of K-shell  = 99
   # of K points = 1270

  Particle set size: 8

 Wavefunction setup: 
 ------------------- 
  Name: psi0
building sposet collection of type bspline
  Created SPOSet builder named 'bspline' of type bspline
  Building SPOSet "spo_ud" with bspline SPOSetBuilder
  TileMatrix = 
 [  2  0  0
    0  2  0
    0  0  1 ]
  Reading 16 orbitals from HDF5 file.
  HDF5 orbital file version 2.1.0
  Reading orbital file in ESHDF format.
  ESHDF orbital file version 2.1.0
  Lattice = 
    [  4.652506  0.000000  0.000000
      -2.326253  4.029188  0.000000
       0.000000  0.000000 18.897261 ]
  SuperLattice = 
    [  9.305011  0.000000  0.000000
      -4.652506  8.058376  0.000000
       0.000000  0.000000 18.897261 ]
bands=8, elecs=32, spins=1, twists=1, muffin tins=0, core states=0
atomic orbital=0
Atom type(0) = 6
Atom type(1) = 6
   Skip initialization of the density
TIMER  EinsplineSetBuilder::ReadOrbitalInfo 0.0050780773
TIMER  EinsplineSetBuilder::BroadcastOrbitalInfo 0.0000219345
Found 1 distinct supercell twists.
number of things
1
1
Super twist #0:  [   0.00000   0.00000   0.00000 ]
  Using supercell twist 0:  [   0.00000   0.00000   0.00000]

and Error is:


Rank =    1  Free Memory = 171949 MB
Rank =    2  Free Memory = 171949 MB
Rank =    3  Free Memory = 171949 MB
Rank =    4  Free Memory = 171949 MB
Rank =    5  Free Memory = 171949 MB
Rank =    6  Free Memory = 171949 MB
Rank =    7  Free Memory = 171949 MB
Rank =    8  Free Memory = 171949 MB
Rank =    9  Free Memory = 171949 MB
Rank =   10  Free Memory = 171949 MB
Rank =   11  Free Memory = 171949 MB
Rank =   12  Free Memory = 171949 MB
Rank =   13  Free Memory = 171949 MB
Rank =   14  Free Memory = 171949 MB
Rank =   15  Free Memory = 171949 MB
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 15 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
ERROR Super twist 0 should own 4 k-points, but owns 1.
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
Fatal Error. Aborting at EinsplineSetBuilder::AnalyzeTwists2
[nia0128.scinet.local:305794] 15 more processes have sent help message help-mpi-api.txt / mpi-a
bort
[nia0128.scinet.local:305794] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help
 / error messages```
kryczko commented 4 years ago

I was getting similar errors when trying to use qbox wavefunctions as well, so it must be the converter that is not working as expected.

jtkrogel commented 4 years ago

What you are seeing could also be caused by a mismatch between the primitive cell kpoints present in the converted orbital file and the supercell you want to use in QMCPACK. Mismatches of this type can happen more easily if you try to specify the k-points to the nscf (for the orbitals) yourself. I would recommend deleting kgrid/kshift input from the nscf run and allowing Nexus to populate it instead.

kryczko commented 4 years ago

I did not alter anything myself. I just ran the script and let nexus handle everything.

jtkrogel commented 4 years ago

Is this the example you are running?

nexus/examples/qmcpack/rsqmc_misc/graphene/graphene.py

kryczko commented 4 years ago

yes

jtkrogel commented 4 years ago

This example is a bit dusty. I'm going to run through it and post any corrections I find.

kryczko commented 4 years ago

Okay, great! Thanks.

jtkrogel commented 4 years ago

Try this one:

#! /usr/bin/env python3

from nexus import settings,job,run_project
from nexus import generate_physical_system
from nexus import generate_pwscf
from nexus import generate_pw2qmcpack
from nexus import generate_qmcpack
from nexus import loop,linear,vmc,dmc

# general settings for nexus
settings(
    pseudo_dir    = '../../pseudopotentials',# directory with all pseudopotentials
    sleep         = 3,                       # check on runs every 'sleep' seconds
    generate_only = 0,                       # only make input files
    status_only   = 0,                       # only show status of runs
    machine       = 'ws16',                  # local machine is 16 core workstation
    )

# generate the graphene physical system
graphene = generate_physical_system(
    lattice   = 'hexagonal',      # hexagonal cell shape
    cell      = 'primitive',      # primitive cell
    centering = 'P',              # primitive basis centering
    constants = (2.462,10.0),     # a,c constants
    units     = 'A',              # in Angstrom
    atoms     = ('C','C'),        # C primitive atoms
    basis     = [[ 0  , 0  , 0],  # basis vectors
                 [2./3,1./3, 0]],
    tiling    = (2,2,1),          # tiling of primitive cell
    kgrid     = (1,1,1),          # Monkhorst-Pack grid
    kshift    = (.5,.5,.5),       # and shift
    C         = 4                 # C has 4 valence electrons
    ) 

# scf run produces charge density
scf = generate_pwscf(
    # nexus inputs
    identifier   = 'scf',           # identifier/file prefix
    path         = 'graphene/scf',  # directory for scf run
    job          = job(cores=16),   # run on 16 cores
    pseudos      = ['C.BFD.upf'],   # pwscf PP file
    system       = graphene,        # run graphene
    # input format selector
    input_type   = 'scf',           # scf, nscf, relax, or generic
    # pwscf input parameters
    input_dft    = 'lda',           # dft functional
    ecut         =  150 ,           # planewave energy cutoff (Ry)
    conv_thr     =  1e-6,           # scf convergence threshold (Ry)
    mixing_beta  =    .7,           # charge mixing factor
    kgrid        = (8,8,8),         # MP grid of primitive cell
    kshift       = (1,1,1),         #  to converge charge density
    wf_collect   = False,           # don't collect orbitals
    use_folded   = True,            # use primitive rep of graphene
    )

# nscf run to produce orbitals for qmc
nscf = generate_pwscf(
    # nexus inputs
    identifier   = 'nscf',          # identifier/file prefix      
    path         = 'graphene/nscf', # directory for nscf run       
    job          = job(cores=16),   # run on 16 cores             
    pseudos      = ['C.BFD.upf'],   # pwscf PP file               
    system       = graphene,        # run graphene                
    # input format selector                                      
    input_type   = 'nscf',          # scf, nscf, relax, or generic
    # pwscf input parameters
    input_dft    = 'lda',           # dft functional
    ecut         =  150 ,           # planewave energy cutoff (Ry)
    conv_thr     =  1e-6,           # scf convergence threshold (Ry)
    mixing_beta  =    .7,           # charge mixing factor
    nosym        = True,            # don't symmetrize k-points
    use_folded   = True,            # use primitive rep of graphene
    wf_collect   = True,            # write out orbitals
    # workflow dependencies
    dependencies = (scf,'charge_density'),
    )

# orbital conversion job for qmc
p2q = generate_pw2qmcpack(
    # nexus inputs
    identifier   = 'p2q',
    path         = 'graphene/nscf',
    job          = job(cores=1),
    # pw2qmcpack input parameters
    write_psir   = False,
    # workflow dependencies
    dependencies = (nscf,'orbitals'),
    )

# Jastrow optimization
opt = generate_qmcpack(
    # nexus inputs
    identifier   = 'opt',           # identifier/file prefix
    path         = 'graphene/opt',  # directory for opt run
    job          = job(cores=16,app='qmcpack'),
    pseudos      = ['C.BFD.xml'],   # qmcpack PP file
    system       = graphene,        # run graphene
    # input format selector   
    input_type   = 'basic',
    # qmcpack input parameters
    twistnum     = 0,               # run at a single twist
    corrections  = [], 
    jastrows     = [('J1','bspline',8),   # 1 body bspline jastrow
                    ('J2','bspline',8)],  # 2 body bspline jastrow
    calculations = [
        loop(max = 6,                        # No. of loop iterations
             qmc = linear(                   # linearized optimization method
                energy               =  0.0, # cost function
                unreweightedvariance =  1.0, #   is all unreweighted variance
                reweightedvariance   =  0.0, #   no energy or r.w. var. 
                timestep             =  0.5, # vmc timestep (1/Ha)
                warmupsteps          =  100, # MC steps before data collected 
                samples              = 16000,# samples used for cost function 
                stepsbetweensamples  =   10, # steps between uncorr. samples
                blocks               =   10, # ignore this  
                minwalkers           =   0.1,#  and this
                bigchange            =  15.0,#  and this
                alloweddifference    =  1e-4 #  and this, for now
                )
             )        
        ],
    # workflow dependencies
    dependencies = (p2q,'orbitals'),
    )

# final dmc run
qmc = generate_qmcpack( 
    # nexus inputs
    identifier   = 'qmc',           # identifier/file prefix       
    path         = 'graphene/qmc',  # directory for dmc run       
    job          = job(cores=16,app='qmcpack'),
    pseudos      = ['C.BFD.xml'],   # qmcpack PP file
    system       = graphene,        # run graphene
    # input format selector                                      
    input_type   = 'basic',
    # qmcpack input parameters
    corrections  = [],              # no finite size corrections
    jastrows     = [],              # overwritten from opt
    calculations = [                # qmcpack input parameters for qmc
        vmc(                        # vmc parameters 
            timestep      = 0.5,    # vmc timestep (1/Ha)
            warmupsteps   = 100,    # No. of MC steps before data is collected
            blocks        = 200,    # No. of data blocks recorded in scalar.dat
            steps         =  10,    # No. of steps per block
            substeps      =   3,    # MC steps taken w/o computing E_local
            samplesperthread = 40   # No. of dmc walkers per thread
            ),                      
        dmc(                        # dmc parameters
            timestep      = 0.01,   # dmc timestep (1/Ha)
            warmupsteps   =  50,    # No. of MC steps before data is collected
            blocks        = 400,    # No. of data blocks recorded in scalar.dat
            steps         =   5,    # No. of steps per block
            nonlocalmoves = True    # use Casula's T-moves
            ),                      #  (retains variational principle for NLPP's)
        ],
    # workflow dependencies
    dependencies = [(p2q,'orbitals'),
                    (opt,'jastrow')],
    )

# nexus monitors all runs
run_project()

# print out the total energy
performed_runs = not settings.generate_only and not settings.status_only
if performed_runs:
    # get the qmcpack analyzer object
    # it contains all of the statistically analyzed data from the run
    qa = qmc.load_analyzer_image()
    # get the local energy from dmc.dat
    le = qa.dmc[1].dmc.LocalEnergy  # dmc series 1, dmc.dat, local energy
    #  print the total energy for the 8 atom system
    print('The DMC ground state energy for graphene is:')
    print('    {0} +/- {1} Ha'.format(le.mean,le.error))
#end if
kryczko commented 4 years ago

This seems to be working fine now. What changes were made and why? If I wanted to turn this into a production calculation, I was thinking I should:

  1. Increase the number of tiles to have the number of atoms >= 50
  2. Increase the ECUT to 200

Do you have other recomendations?

jtkrogel commented 4 years ago

In both the original and revised example, the Jastrow factor optimization was performed at a single supercell twist.

The original example attempted to create orbitals for the optimization via a separate nscf calculation, but it erroneously generated orbitals for a single primitive cell twist (via the user provided kgrid/kshift to nscf), not enough to populate the supercell.

The revised example uses a single nscf to create orbitals for all supercell twists relying on the information provided to generate_physical_system. The Jastrow optimization still occurs at a single supercell twist (the twist indexed "0" in QMCPACK, requested via setting twistnum=0).

kryczko commented 4 years ago

how would you approach doing twist averages?

jtkrogel commented 4 years ago

For production calculations, you have identified important issues.

First you should converge the quality of the orbitals. You can do part of this this e.g. by increasing ECUT until the DFT energy changes by less than 1 meV per formula unit. Next, optimize a Jastrow factor using meshfactor=1.0 and then change the meshfactor (over a range between 0.8 and 1.6) until the VMC variance is roughly converged. The meshfactor corresponds to the fineness of the 3D B-spline grid used to represent the orbitals in real space within QMCPACK.

Next, converge the supercell twist grid for a supercell of interest. Do this by observing the VMC total energy per formula unit as a function of twist grid density. Use the converged twist grid density for all subsequent supercell calculations (which may be of different sizes/tilings).

Finally, account for finite size effects. Extrapolation over multiple cell sizes may be needed. This is easier in bulk solids (single variable extrap) that in quasi-2D systems (two variable extrap). Using finite size corrections schemes can accelerate the convergence here, and sometimes eliminate the need for extrapolation. Good finite size corrections to consider are the KZK DFT-based corrections (difference between LDA and KZK functional energies) and the many-body S(k) based corrections (see SkAll estimator in QMCPACK). I would explore finite size extrapolation/correction in VMC before moving on to more expensive DMC calculations.

kryczko commented 4 years ago

This is great information. Thank you very much. How do you suggest including estimators with nexus? I am digging through the source code but I am unsure how to import the estimators to be appended in a list for generate_qmcpack.

jtkrogel commented 4 years ago

If you want to include the density for example, you can do this:

from qmcpack_input import spindensity

qmc = generate_qmcpack(
    ...
    estimators = [spindensity(grid=(40,40,40))],
    ...
     )

It is sometimes nice to inspect the qmcpack xml input directly before running. You can do it this way:

print(qmc.input.write())
exit()
kryczko commented 4 years ago

great thanks so much!!

kryczko commented 4 years ago

How would I also include twist averaging? I see that there is energydensity as well. How would I calculate the energy density on an equivalent grid? I.e. a 40x40x40 uniformly spaced grid?

jtkrogel commented 4 years ago

For twist averaging, you can just increase the kgrid within generate_physical_system. This is the supercell twist grid following tiling. Nexus will automatically make one input file per twist, similar to this example:

https://github.com/QMCPACK/qmcpack_workshop_2019/tree/master/day2_nexus/quantum_espresso/03_qe_diamond_dft_vmc_twistavg

One point I should bring up in this context is that some care will be needed for twist averaged DMC calculations. For any DMC calculation, population control bias starts to become an issue with a small enough walker population. Provided you have 1000-2000 total walkers per twist in the DMC calculations, you shouldn't have to worry about this.

You could specify the energy density XML element w/o much difficulty, but the estimator is still not functional.

kryczko commented 4 years ago

Great, thank you for this information. Do you have an estimate for when the energy density estimator will be functional?

jtkrogel commented 4 years ago

No, unfortunately I do not.

kryczko commented 4 years ago

That is too bad, I would really like to have that part of the code functioning!!

prckent commented 4 years ago

Noted!

gokaran commented 3 years ago

Dear Krogel, How one can specify starting magnetization (let say LiH, where let say one want to put 0.5 instead of zero in scf or nscf section) of heterogeneous system in Nexus? .Writing simply starting_magnetization(1)=0.5 give error.

Many thanks in advance

Yours Gokaran Shukla

jtkrogel commented 3 years ago

@gokaran please start a new issue with your question. I will be happy to answer it there.

jtkrogel commented 3 years ago

The many issues discussed here have been resolved.