Stochastic GW Manual V1.0

User Manual: Pdf

Open the PDF directly: View PDF .
Page Count: 26

Download
Open PDF In Browser	View PDF

stochasticGW – User manual and Tutorial (version 1.0)

Vojtěch Vlček,1,2 Christopher Arntsen,3 Wenfei Li,1 Roi Baer,4
Eran Rabani,5,6 and Daniel Neuhauser1
1

Department of Chemistry and Biochemistry,
University of California, Los Angeles California 90095, USA.
2

After July 1st, 2018: Department of Chemistry and Biochemistry, University
of California, Santa Barbara California 93106, USA
3

Department of Chemistry,
Youngstown State University, Youngstown, Ohio 44555
4

Fritz Haber Center for Molecular Dynamics, Institute of Chemistry,
The Hebrew University of Jerusalem, Jerusalem 91904, Israel
5

Department of Chemistry,
University of California, Berkeley, California 94720, USA
6

The Sackler Center for Computational Molecular Science,
Tel Aviv University, Tel Aviv 69978, Israel

Contents
1 Introduction

3

2 License

4

3 Overview

5

4 Code status

7

5 Downloading, compilation and installation

8

5.1

Downloading and installation . . . . . . . . . . . . . . . . . . . . .

8

5.2

Compilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

5.3

Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

6 Input files
6.1

Brief overview of the input files . . . . . . . . . . . . . . . . . . . . 10

6.2

Detailed description of the files and flags . . . . . . . . . . . . . . . 10

7 Output

18

8 Tutorial

20

8.1

H2 : Combined Stochastic and Deterministic treatment . . . . . . . 20

8.2

C60 : Fully stochastic calculation . . . . . . . . . . . . . . . . . . . 23

9 Acknowledgements

2

10

26

stochasticGW - User manual and Tutorial

1

Introduction
This manual details stochasticGW , a software for large scale linear-scaling
G0 W0 stochastic calculations. The code combines the success of the G0 W0 approximation, a quantitative approach for vertical ionization energies and electron
affinities, with the ability of stochastic methods to tackle very large systems. The
methodology was introduced and detailed in several recent references, especially:
Ref. I :

→ Sections: 3, 5 & 8

D. Neuhauser, Y. Gao, C. Arntsen, C. Karshenas, E. Rabani and R. Baer,
Breaking the theoretical scaling limit for predicting quasiparticle energies:
The stochastic GW approach, Phys. Rev. Lett., 113, 076402 (2014).
https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.113.076402

Ref. II :

V. Vlcek, E. Rabani, D. Neuhauser and R. Baer, Stochastic GW calculations
for molecules, J. Chem. Theo. Comput. 13, 4997 (2017).
https://arxiv.org/abs/1612.08999

Ref. III :

V. Vlcek, W. Li, R. Baer, E. Rabani and D. Neuhauser, Turbo-boosting
GW beyond > 10,000 electrons using fractured stochastic orbitals, to be
placed on the Arxiv (2018)

If this is your first time using this manual, read the Overview and Downloading
and compilation section, and then jump to the Tutorial.

stochasticGW - User manual and Tutorial

3

2

License
Unless otherwise stated, all files distributed in this package are licensed under
the following terms:
stochasticGW
Developed by:
stochasticGW, University of California, Los Angeles
http://www.stochasticgw.com/
Permission is hereby granted, free of charge, to any person obtaining a copy of the
stochasticGW software and associated documentation files (the “Software”), to
deal with the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimers.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimers in the documentation
and/or other materials provided with the distribution.
3. Neither the names of stochasticGW , stochasticGW, University of California, Los Angeles, nor the names of its contributors may be used to
endorse or promote products derived from this Software without specific
prior written permission.
THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY
KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
WITH THE SOFTWARE.
You are under no obligation whatsoever to provide any bug fixes, patches, or
upgrades to the features, functionality or performance of the source code (“Enhancements”) to anyone; however, if you choose to make your Enhancements
available either publicly, or directly to stochasticGW University of California,
Los Angeles, without imposing a separate written license agreement for such Enhancements, then you hereby grant the following license: a non-exclusive, royaltyfree perpetual license to install, use, modify, prepare derivative works, incorporate
into other computer software, distribute, and sublicense such enhancements or
derivative works thereof, in binary and source code form.

4

stochasticGW - User manual and Tutorial

3

Overview
stochasticGW calculates a matrix element of the self-energy in the G0 W0 approximation
E
D
(1)
Σ (ω) = φ Σ̂ (ω) φ ,
where φ is a single particle orbital of the molecule studied (typically a HOMO or
LUMO), associated with a DFT Hamiltonian Ĥ; Σ̂ is the self-energy operator,
composed of a static exchange component (X) and a polarization part (P )
E
E D
E D
D
(2)
φ Σ̂ (ω) φ = φ Σ̂X φ + φ Σ̂P (ω) φ ,
where the latter is a Fourier transform
E
D
E Z
D
φ Σ̂P (ω) φ = eiωt φ Σ̂P (t) φ dt.

(3)

In the G0 W0 approximation the coordinate representation of Σ̂P (t) is extremely
simple:


ΣP (r, r0 , t) = G0 (r, r0 , t) WP r, r0 , t+ ,
(4)
where we introduced the non-interacting Green’s function and the dynamical part
of the effective interaction.
The starting point (used for G0 and W ), is a Kohn-Sham Hamiltonian which
is constructed from ground state orbitals φ. At the moment (version 1.0) the
stochasticGW code is limited to an LDA starting point – this restriction will
be lifted in the upcoming versions of the code. Orbitals φ are represented on a
real-space grid and have to be supplied by the user.
The matrix element of the self-energy is evaluated as a statistical average over
multiple stochastic realizations - the number of Monte Carlo samples. The
stochastic approach uses multiple resolutions of the identity with random functions for several purposes:
• Separating the product of G0 and WP , converting the calculation of ΣP
into evaluating a matrix element of WP .
• Evaluating a matrix element of WP R (where the R subscript indicates a
causal, i.e., retarded effective interaction) by stochastic TDH (time-dependent
Hartree, yielding RPA) or stochastic TDDFT .
• Converting the matrix elements of WP R to matrix elements of WP (the timeordered effective interaction) through a stochastic fragment basis approach
(Ref. III).
The eventual method scales practically linearly (or better) with system size; for
details of the method and scaling read Refs. I & II and, for the closest description
to the present version, Ref. III
The
D programEoutputs the matrix element of the polarization self-energy,
φ Σ̂P (ω) φ , over a wide range of frequencies, as well as the exchange selfD
E
energy φ Σ̂X φ . The code also prints the estimate of the quasiparticle energy
stochasticGW - User manual and Tutorial

5

of state φ for the given number of Monte Carlo samples, obtained from solving
the quasiparticle equation:
E
D
(5)
ε = εKS + φ Σ̂X + Σ̂P (ε) − v̂xc φ ,
where εKS is Kohn-Sham eigenvalue of the eigenstate φ and v̂xc is the (meanfield) exchange-correlation potential. The input and output are detailed in the
sections below.

6

stochasticGW - User manual and Tutorial

4

Code status
Several items in this 1.0 version have not been implemented or tested. Specifically,
note that:
• The code has not been tested for open-shell systems.
• The core-shell exchange correction is implemented but not tested.
• Only a local LDA exchange-correlation potential is implemented.
• The numbers of points in each dimension, nx, ny, nz need all to be even.
• The maximum angular momentum in the nonlocal pseudopotentials needs
to be 1 or 2 or 3. If it is 1 then it is necessary to have lloc=2.

stochasticGW - User manual and Tutorial

7

5
5.1

Downloading, compilation and installation
Downloading and installation
stochasticGW is a Fortran code requiring Open MPI parallelization; The tgz
tar file is available from:
https://github.com/stochasticGW/stochasticGW
For installation, first unpack the archive:
tar -xzvf sgw-v1.0.tgz
Upon unpacking, you will get a directory sgw-v1.0 with a README file and 3
subdirectories:
examples/

PP/

README

src/

The installation of the package is detailed then in README and is repeated here.

! → An installation of the FFTW library (version 3.1 or higher) is a prerequisite; it can
be downloaded from http://www.fftw.org.
5.2

Compilation
In the src directory you will find a makefile that needs a little modification
depending on your system settings, as described below.
Edit makefile and specify your parallel compiler (FCMPI) and FFTW library (FFTFLG).
The flags for the parallel compiler (MPIFLG) can usually be left unaltered.
The supplied header of makefile reads:
FCMPI
MPIFLG
FFTFLG

= mpif90
= -DMPI -O3
= -lfftw3

The stochasticGW code is built simply by writing:
make
in the src directory.
After successful compilation, the stochasticGW executable, sgw.x, will be found
in the src directory.
Make sure sgw.x has been created
ls -l sgw.x
5.3

Implementation
Now you are ready to implement the code. Then, say you want to run it in a
directory you will call ∼ /mygwrun. Go there, and copy the H2 example input,
the executable, and pseudopotentials to this directory:
mkdir ~/mygwrun
cd ~/mygwrun
cp ~/sGW_v1/examples/H2/* .
cp ~/sGW_v1/src/sgw.x .
cp -r ~/sGW_v1/PP/ .

8

stochasticGW - User manual and Tutorial

At this point, your mygwrun directory should contain the following:
cnt.ini counter.inp INPUT log_h2example PP random.inp sgw.x wf.txt
where PP is a directory (with pseudopotentials). All these files are inputs except
for log_h2example.
To run the H2 example (from the Tutorial) on, say, 30 cores, write
mpirun -n 30 ./sgw.x >& log &
The output file log should match log_h2example.
After the H2 test is finished, you can run any other molecule/material you want.
You’ll need to change a few items, detailed below and in the Tutorial.
• Change INPUT as needed – especially make sure that you change ntddft
to be between 8 and 30.
• Change the atomic coordinates (cnt.ini).
• Bring your own wf.txt (or better yet, wf.bin ) file that contains the occupied DFT energies and functions.
• And remember to include in PP/ any pseudopotential you need.

stochasticGW - User manual and Tutorial

9

6
6.1

Input files
Brief overview of the input files
Two groups of input files are required in the working directory. One set includes
general parameters, and the other system-specific files.

Seed, labels and flags
random.inp :
counter.inp :

INPUT :

seed file for the random number generator.
a label file with a counter for modifying the random number input and
for labeling the input. Using this file simplifies working simultaneously on
multiple runs.
general input file with flags and parameters.

System specific files
cnt.ini :

6.2

atomic position file (in Angstroms or atomic units).

*UPF/*fhi :

pseudopotential files (only two formats are currently supported). These
files should be placed in a subdirectory of the working directory, PP/

wf.txt or wf.bin :

text or binary file with all occupied (and possibly a few unoccupied) wavefunctions.

orb.txt (optional) :

an orbital wavefunction (“φ” in the previous discussion) for which the selfenergy matrix element is needed. Usually this file will not be given and
instead φ will be extracted from wf.txt or wf.bin.

Detailed description of the files and flags

FILE: random.inp
A file with 6 lines, each with 4 integer numbers (Fortran format integer*8).
These are the seeds for the KISS random-number algorithm. Each line in is a
seed for a different random variable used in the StochasticGW algorithm.

FILE: counter.inp
A file with a single integer that labels the output directories and files and also
modifies the seed. It alleviates the need for different runs to use different random.inp files when the calculation is repeated to reduce the statistical error (see
example below)

10

stochasticGW - User manual and Tutorial

Example
Let’s say we want nmctot=2000 stochastic samples.a This can be achieved
by a single run with 2000 stochastic samples, or perhaps, due to computer
jobs scheduling issues, it may be easier to submit 5 different jobs, each with
400 stochastic samples.
For these 5 jobs we would need, in principle, 5 different random.inp files with
totally different values, which is cumbersome. We therefore use counter.inp.
Each counter modifies the values of the variables from random.inp. So in
our example we would only need to submit 5 different jobs, each with the
same random.inp file, and in each run the counter.inp file contains a different integer from, say, 0 to 4 (or any other 5 different integers).
Then, at the end of the 5 runs, average the 5 output quasiparticle energies
or 5 output self-energies files.
a

see also the buffer_size variable, which controls the parallelization of multiple stochastic states.

FILE: INPUT
Main file for flags and parameters. For most variables the default values are
appropriate. Examples are shown in the Tutorial.
Specifically, only one set of variables has to be specified in INPUT: the names of
the pseudopotentials (variable: pps). Other parameters need to be included in
the INPUT file only if they differ from their default values. The format for most
variables (except for the pseudopotential names) is:
variable name


The order is not important, but no variable should be repeated; lines starting
with # are ignored. The input variables are:

pps :

a label indicating the beginning of the list of pseudopotentials. After the
last pseudopotential put a line starting with #. The pseudopotentials listed
must be present in a subdirectory of the working directory (where the INPUT
file is) called PP/. It is recommended to use well-tested pseudopotentials.
At present, the supported formats are *UPF and *fhi. In addition, the
code currently supports only the LDA functional. In the next versions,
other functionals will be included too.

! →

Note that the program will issue a warning if the pseudopotential wavefunctions are not properly normalized, but it will then normalize them properly.

scratch :

(default: GW_SCRATCH) the prefix of the name of the scratch directory; the
scratch directory is then appended by the counter (from the counter.inp
file).

stochasticGW - User manual and Tutorial

11

Example
If counter.inp contains the number 0, and the relevant segment from
INPUT is
scratch
/scratch/john/GW_SCRATCH

then the scratch directory will be /scratch/john/GW_SCRATCH.0
The scratch directory contains several intermediate files. It should be
ideally not placed on the main shared disk but on fast disks connected
to each core, if such disks are available.
Each core will try to create this scratch directory, unless it was created
earlier by another core accessing the same scratch disk.
nmctot :

(default: 1000) the number of stochastic realizations used for statistical averaging (Monte Carlo iterations). Typically 300-2000 iterations are needed
for acceptable convergence. The error in the output self-energy and quasi1
particle energies decreases as √
.
nmctot
If nmctot is larger than the number of MPI processes Nproc the calculation will be repeated (nmctot/Nproc ) times. For instance: if you set
nmctot=1000 and distribute the job over Nproc =100 cores, the calculation
will be repeated 10 times on each core.
If needed, the program increases nmctot so that mod(nmctot, Nproc ) = 0.

binary :

(default: .TRUE.) a logical variable; designating whether the occupied wavefunctions are stored in a binary file wf.bin or a text one wf.txt.

periodic :

(default: .FALSE.) a logical variable designating whether the calculation is
periodic (restricted to Γ-point sampling).

box :

(default: pos) a character*3 variable designating the way the calculation
box is constructed. Two values are allowed, pos and sym .
The first choice, pos, implies that the calculation box is (in Bohr!)
[0,nx*dx][0,ny*dy][0,nz*dz]. The wavefunctions in wf.txt are then
interpreted to start at the corner of this box, (0,0,0). Such box convention
is used in Quantum Espresso and most planewave codes. Note that in our
code the coordinates in the cnt.ini file need to fit in the box (or more
preceisely, if they are given in Angstrom, their values in Bohr need to fit in
the box). If they do not fit in the box the program will stop.
The second, sym, assumes that the box of the calculation is
[-nx*dx/2,nx*dx/2][-ny*dy/2,ny*dy/2][-nz*dz/2,nz*dz/2]
(again in Bohr), i.e., its center is the origin. The wavefunctions in wf.txt
are again interpreted to start at the corner of the box, which is now
(-nx*dx/2,-ny*dy/2,-nz*dz/2). Again, the coordinates in cnt.ini need
to fit in the box (once expressed in Bohr). This symmetric box convention
is more natural for non-periodic molecules.

multiplicity :
12

(default: 1) either 1 or 2; 1 indicates a closed-shell calculation, 2 an open
stochasticGW - User manual and Tutorial

shell. Note: the program was not tested for open-shell systems,
and at this stage should be trusted only for closed-shell calculations.
ekcut :

(default: 25.0) a cutoff (in Hartree) on the kinetic energy. Even though
there is a formal maximum on the kinetic energy on the grid (and the grid
is supplied elsewhere, see below) it is numerically better to use a lower
cutoff. In practice one should use a value somewhat higher than the depth
of the pseudopotentials. For calculations with a new element check the
convergence of the calculations with ekcut=15 to 30 Hartree.

gamma :

(default: 0.06) energy broadening parameter for ΣP (Section 3 in Ref. II and
the supplementary material in Ref. I). Generally 0.06 should be fine. An
increased gamma reduces the computational effort and improves somewhat
the statistics, but if gamma is too large the broadening will be too severe
leading to wrong quasiparticle predictions.

orb_indx :

(default: -1). If orb_indx>0 its value designates the orbital φ which will
be used in the calculation of the self-energy (e.g., if orb_indx=5 then φ is
the 5th eigenstate). If orb_indx is -1 then φ should be read from orb.txt .

ntddft :

(default: 15) the number of stochastic states used in the time-dependent
RPA or TDDFT calculation (in the W part of GW ). For very small systems,
use ∼ 30. For large systems with hundreds of electrons or more, it is enough
to use ∼ 8 orbitals.
For open-shell systems (multiplicity=2) ntddft denotes the number of
stochastic TDRPA/TDDFT states for each spin.

! →

For a non-stochastic TDRPA/TDDFT propagation, use ntddft = -1. Then
a deterministic TDRPA/TDDFT calculation is used for calculating the action of W . In this case all occupied orbitals are excited and propagated.
This choice should only be done for small systems, since it increases the
computational time considerably. A non-stochastic calculation of the action of W is neither feasible nor needed for large systems, since the stochastic calculation of the action of W is getting progressively better for bigger
systems.

units_nuc :

(default: Bohr) either Angstrom, or Bohr. Controls the units for the nuclear
positions input (in the cnt.ini file).

flg_dyn :

(default: .FALSE.) logical variable; .FALSE. implies an RPA calculation;
.TRUE. yields a TDDFT calculation of the action of W where the exchangecorrelation potential is updated at each time-step.

The following variables could usually be safely kept at their default values: :
dt :

(default: 0.05) time step in atomic units for the TD propagation of W (note:
0.05 a.u. ' 1.2 attosec). If desired, dt could be reduced to 0.03, although
0.05 generally works well.

nxi :

(default: 10,000) Number of random spanning vectors (denoted Nξ in Ref.III). Usually there’s no need to change the default.

segment_fraction :

(default: 0.003) The ratio of the size of each stochastic vector ξ to the total
grid length. See Ref. III. Usually there’s no need to change the default.

stochasticGW - User manual and Tutorial

13

sm :

(default: 0.0001) The perturbation parameter (denoted τ in Eq. 29 in
Ref. II). The results should not vary for sm in the range 10−5 –10−3 .

scale_vh :

(default: 2; but will be set to 1 by the program if periodic) allowed values
1 or 2. Controls whether a Martyna-Tuckerman1 grid-doubling is used in
the TDRPA/TDDFT calculations. If the calculation is periodic, scale_vh
would be set in the program to 1 regardless of the input value. For nonperiodic molecules, scale_vh formally needs to be set at 2. So since the
program automatically sets scale_vh, there is usually no need to touch this
variable.

nrppmx :

(default: 1200) an upper bound on the number of radial grid points for
the pseudopotential. No need to change unless a new pseudopotential with
more than 1200 grid points is used.

buffer_size :

(default: 1) the number of cores used in each independent stochastic runs.
For intermediate or large size systems (hundreds or thousands of electrons)
choose the default (1) which uses the least total CPU resources. But for
very large systems, the wall-time of the calculation could be too long when
using a buffer_size of 1. Generally, it is a good idea to use a buffer_size
which either equals to the number of cores on each CPU, or divides it (e.g.,
if there are 12 cores on each CPU, use a buffer_size of 1,2,3,4,6 or 12).
Also, if buffer_size is >1, aim at having mod(ntddft+1,buffer_size)=0.
For instance: if buffer_size is 12, use ntddft=11, 23 or 35. If ntddft is
bigger the convergence is somewhat faster, but the effort grows, so values
in the range 8-30 are usually the most efficient.
Example
Let’s say that a giant system requires about 1000 Monte Carlo iterations. One could use 1000 cores and then each core will be used for
a single Monte Carlo iteration. But say that such a single calculation
took 2 days (i.e., each core takes 2 days), while the computer-center allocation limits each run to be 1 day long. In that case, it is necessary
to spread each single calculation to several cores. For example, use
then a buffer_size of 4 and 4000 mpi cores to do the 1000 MonteCarlo iterations, and each group of 4 cores will do a single stochastic
run. The calculation will then take a shorter wall-time.
Note: to spread the work in each separate run to several cores, we
use the mpi_comm_split routine, which divides the cores to subsets,
usually labeled as colors. Further details can be found on the web.
a
a

1

14

a good introduction is in http://www.bu.edu/tech/support/research/
training-consulting/online-tutorials/mpi/more/mpi_comm_split/

G.J. Martyna, and M.E. Tuckerman, “A reciprocal space based method for treating long range
interactions in ab-initio and force-field-based calculation in clusters”, J. Chem. Phys. 110,
2810 (1999),doi:10.1063/1.477923

stochasticGW - User manual and Tutorial

FILE: cnt.ini
Nuclear positions in a Cartesian format:








The element name is in a 1- or 2-character format (e.g., H or Si). The nuclear
units could be in Angstrom or Bohr (the units are specified by the units_nuc
variable, the default of which is Bohr), but recall that all other quantities in
the code are calculated in atomic units, i.e., Hartree for energies and Bohr for
distances.
FILE: wf.txt
A wavefunction input file (required if binary = .FALSE.) with all occupied (and
potentially a few unoccupied) eigenstates. The first 6 lines detail the grid used
(lengths in Bohr), followed by the multiplicity (which needs to match the default
value or, if given, the value in INPUT), and the Kohn-Sham eigenvalues. Then
each wavefunction is given, preceded by the state-indicies (and possibly the spinindex, if multiplicity=2). The file contains the state index on the 12th , 14th , . . .
lines:
nx
ny
nz
dx
dy
dz
nsp
nstates
evls

orbitals
1
<1st orbital>
2
<2nd orbital>
..
.













Note that the list of Kohn-Sham eigenvalues should include all "nstates" eigenvalues.
The program determines the number of states to be read based on the ionic
charges provided in the pseudopotential files. If the normalization of any orbital
is incorrect a warning is issued but the program continues.
The orbital is given as usual in Fortran with the left-most index (x) varying the
fastest. Each orbital should contain nx · ny · nz points.
If the multiplicity is 2 (open-shell), the spin-index of each orbital should be 1 or
2 (indicating an α or β spin). The spin-index is not required if the multiplicity
is 1 (closed-shell).

stochasticGW - User manual and Tutorial

15

FILE: wf.bin
The wavefunction input file (required if binary = .TRUE.). Analogous to wf.txt.
When preparing this binary file from your favorite DFT program output, you
should follow the example of the following code segment:
real*8 hw(1:nstates)
real*8 orbitals(1:ngrid,1:nstates)
...
open(9,file=’wf.bin’,status=’old’,form=’unformatted’)
write(9)’nx
’,nx
write(9)’ny
’,ny
write(9)’nz
’,nz
write(9)’dx
’,dx
write(9)’dy
’,dy
write(9)’dz
’,dz
write(9)’nsp
’,nsp
write(9)’nstates ’,nstates
write(9)’evls
’
write(9)(hw(i),i=1,nstates)
write(9)’orbitals ’
do i=1,nstates
! ispin: 1 for closed shell; 1 or 2 for open shell.
write(9)i, ispin(i)
write(9)u(1:ngrid,i)
enddo
close(9)
Each text string in lines 1-9 and 11 is of length 9.
At present, the user needs to prepare wf.bin, but in future versions we would
supply preparation routines that process the output of Quantum Espresso, Abinit,
and other DFT programs and produce the proper wf.bin file.

FILE: orb.txt
Orbital φ for which the quasiparticle energy is sought. This orbital is usually
taken from a DFT code. Typically it will be the HOMO or the LUMO. The file
is organized similarly to the wavefunction file, and the heading values need to
match those in wf.txt or wf.bin
nx
ny
nz
dx
dy
dz
nsp
orb

16











stochasticGW - User manual and Tutorial

The program issues an error warning if the normalization of φ is wrong (but will
not stop).
If the orb_indx variable is specified, with a value bigger than 0, then orb.txt is
not needed and φ will be taken directly from the wavefunction file.

FILE: *UPF/*fhi
A pseudopotential file is needed for each of the elements in the molecule (i.e.,
each of the different element in cnt.ini). The names of the pseudopotential
files are specified by variable pps. The pseudopotentials should be placed in a
subdirectory, PP/
At present the program only handles pseudopotentials with a single KleinmanBylander2 state for each angular momentum. Note that such pseudopotentials
are known to have ghost states for heavier elements and some choices of specific
pseudopotentials. 3 The user should verify that there are no ghost states for
each of the pseudopotentials. This can be verified, for example, with a DFT code
that uses the Kleinman-Bylander decomposition (e.g., Quantum Espresso). Such
a DFT code should be used for checking for ghost states in very small systems
(atoms or small molecules) each of which contains one or more of the atoms in
the system with its associated pseudopotentials.
Also note that we have not extensively checked/verified pseudopotentials with
core-charge exchange corrections.
Future versions of the code would incorporate more general pseudopotentials.

2

L. Kleinman and D. M. Bylander, “Efficacious form for model pseudopotentials.",
Phys. Rev. Lett. 48, 1425 (1982), https://doi.org/10.1103/PhysRevLett.48.1425
3 A specific example for a problematic pseudopotential is the usual one for hydrogen when
lmax=3 is used. Therefore, for H use lmax=1 and lloc=2, see the H.pw-mt_fhi.UPF pseudopotential in the supplied PP directory

stochasticGW - User manual and Tutorial

17

7

Output
During the run a standard output (log file) details the input variables, possible errors, warnings, timing and a few additional details, as well as the output
energies. We separately supply the self-energy as a function of frequency.
If the calculation is repeated for several Monte-Carlo steps on each core, the
variable designating the current step has the range:
MC_STEP = 1, ...,

nmctot × buffer_size
ncores

(6)

(see nmctot and buffer_size for details). The quasiparticle energy is estimated
at each cycle together with its stochastic error (see the Tutorial).
Here is an example of the tail of the log_h2example file from the Tutorial, for
an MPI run using n_cores = 30 cores and with buffer_size = 1.
######## Accumulative # of Monte Carlo steps : MC=
Completed Monte Carlo steps in each core : MC_STEP=

120
120
120
120
120
120
120
120
Time:
Time:

MC_STEP:
program end

-0.378696
-0.424763
-0.646902
0.018595
-0.018666
-0.000045
0.894364
-0.594063
-0.582239
0.003750

0.000000
0.012778
0.008620
0.000009

0.012778
0.000543

4 => time:
306.199127

120 ########
4

E

MC,
MC,
MC,
MC,
MC,
MC,
MC,
MC,


Sig_R(e_qp)
Sig_R(e_ks)
Sig_I(e_ks)
Z
Elin
E_qp
Ei_qp

+-stat.err
+-stat.err
+-stat.err
+-stat.err

+-stat.err
+-stat.err

306.175140

The tail details the MC_STEP = 4th calculation step, so
nmctot = 120 = MC_STEP × n_cores total calculations were done.
The DFT orbital energy is -0.378696, followed by  ≡ hφ|v̂xc |φi.
The following lines are all in the same format: number of Monte Carlo samples (
120 – denoted MC; since this is the final step of the calculation, it coincides with
the value of nmctot); variable name, its statistical error and labels. Specifically:
D
E
• Expectation value of the non-local Hartree-Fock exchange, < X > = φ Σ̂X φ ,
calculated deterministically, without statistical error.
• Sig_R is the real part of the self-energy evaluated either at the quasiparticle
energy (e_qp) or at the starting point KS eigenvalue energy (e_ks).
• Sig_I is the imaginary part of the self-energy.
• Z is the renormalization factor.
• Elin is the often-used but crude linear extrapolation quasiparticle energy
estimate based on a renormalization factor Z.
• Finally, E_qp and Ei_qp are the real and imaginary parts of the quasiparticle
energy obtained by properly solving Eq. 5.
18

stochasticGW - User manual and Tutorial

The last two lines provide CPU wall-timing in seconds.
The quasiparticle energy can also be read directly from the header of the selfenergy curve which is provided in the file SigW.txt in the GW_OUTPUT.X directory
(where X is the number in counter.inp file). If the run crashes for some reason
prior to completion, SigW.txt will contain the most up-to-date information.
Example
Here is a sample SigW.txt file header.
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#

SigW output from stochastic GW
8192 plotting frequencies
MC_STEP (Monte Carlo steps on each core)
4
N_cores
30
Buffer_size (# cores per indep. run)
1
MC (accumulative # Monte Carlo runs) =MC_STEP*(N_cores/buffer_size)
120
KS energy
-0.378696352

-0.646901488

d
0.00000000


-0.424763411
e_qp_real,
-0.582239330

de_qp_real
1.27781946E-02

e_qp_imag,
de_qp_imag
3.74985230E-03
5.42658614E-04
w

SigW_real

SigW_imag

dSigW_real

dSigw_imag

The header first reports that the file contains the self-energy at 8192 frequency points. It then repeats the information from the output log file on
the number of MC steps, the Kohn-Sham energy, the expectation values, the
polarization self-energies and the quasiparticle energies. The rest of the file,
after the header, contains the frequency-dependent polarization self-energies
and their statistical deviation.
Each line of the header is introduced with #, so that it is ignored by many
plotting programs. Therefore, the polarization part of the self-energy contained in the SigW.txt file can easily be visualized (see the example in the
Tutorial).
Further, intermediate information is placed in the directory GW_WORK.X (where X
stands for number specified in counter.inp); it contains files similar to SigW.txt
taken at intermediate number of Monte Carlo iterations. It also contains a file
(details_output.txt) with detailed information on the run that in case of trouble could be useful for debugging. The scratch directory (specified by the scratch
variable) contains intermediate files that can be safely discarded after the run.

stochasticGW - User manual and Tutorial

19

8

Tutorial
In this section we illustrate using two examples how to execute the program, and
discuss the meaning of the most important input variables.

! → The examples presented in this section were selected to quickly illustrate how
the stochasticGW code works. For production runs the user should carefully
investigate the convergence with respect to the grid size, spacing and other parameters.

8.1

H2 : Combined Stochastic and Deterministic treatment
In this example, we calculate the ionization potential of an H2 molecule. Note that
calculations for small molecules require relatively much more effort vs.
the extremely large systems for which stochasticGW was developed.
The DFT estimate of the H2 ionization potential is −10.3 eV, which is too
low compared with the experimental result4 (−15.4 eV) and is corrected by
the GW calculation in the next step. Before we can proceed with running the
stochasticGW code, we need the occupied wavefunctions, wf.txt, which is
provided in examples directory. User can obtain it from other codes, such as
Quantum Espresso5 or others.

(i) running the stochasticGW code
The example INPUT file is:
pps
H.pw-mt_fhi.UPF
#
scratch_path
.
ekcut
15d0
nmctot
120
ntddft
-1
gamma
0.06
binary
F
4
5

20

http://webbook.nist.gov/cgi/cbook.cgi?ID=C1333740&Mask=20
http://www.quantum-espresso.org/wp-content/uploads/Doc/INPUT_PP.html

stochasticGW - User manual and Tutorial

orb_indx
1
scale_vh
1
box
sym

! → The user should be familiar with all the input variables listed here. Note specifically:
• We used scale_vh=1 to expedite the calculation, although formally it needs
to be 2.
• We set ntddft=-1 to indicate that the time-dependent Hartree (RPA) calculation of the effect of WP should be deterministic, using here the single
occupied state.
As explained in Section 5, the working directory should contain and INPUT
and several other input files: wf.txt, cnt.ini, counter.inp, random.inp as
well as the PP/ pseudopotential directory that needs to contain (at least) the
H.pw-mt_fhi.UPF hydrogen pseudopotential file.
For simplicity we assume that the working directory contains also the sgw.x
executable.
The stochasticGW calculation is then run; for example, to run on 120 MPI
cores, write
mpirun -n 120

./sgw.x

>& log

The log file prints all parameters and information about all files as well as the
starting point energy in Hartree units :
 from orbital in orb.txt is:
vs. eorb as read:

-0.378696352
-0.378631353

This orbital energy slightly differs from the DFT input value; this is because the
Hamiltonian in the underlying DFT calculation slightly differs from the Hamiltonian constructed in stochasticGW .
After the calculation finishes, we look at the end of the log file where the quasiparticle energy is provided :
120

-0.582239

0.012778

MC, E_qp

+-stat.err

This value, which is the solution of Eq. 5, is in the right ballpark (−15.8 ± 0.3 eV
compared to the experimental value of −15.4 eV), but shows significant
√ statistical
error. Since the statistical error is generally proportional to 1/ nmctot, it is
necessary to use nmctot ∼ 120 ∗ (0.3/0.05)2 ∼ 4000 Monte Carlo iterations for a
statistical accuracy of 0.05 eV.
The program finds the ε that solves Eq. 5 self-consistently, automatically, i.e.,
since Σ (ω) is given at all frequencies the program searches for the energy (ε)
stochasticGW - User manual and Tutorial

21

where Eq. 5 is fulfilled (or more precisely, the it searches for the closest solution
to the Kohn-Sham energy).

It is instructive however to graphically do this procedure, by plotting the real-part
of the self-energy from the SigW.txt located here in the GW_OUTPUT.0 directory
(where 0 is the number in the counter.inp file).

The SigW.txt file has a header that, for completeness, essentially duplicates the
output in the log file, as presented in Section 7. The header was presented in the
8.

After the header the frequency dependent polarization part of the self-energy is
given. We concentrate on the first and second columns showing to the frequency
(w) and the real part of the self-energy (SigW_real, i.e., Σp (ω)).
Specifically, we use gnuplot6 to find the graphical solution to Eq. 5. As mentioned,
lines starting with # are ignored. In the GW_OUTPUT.X directory start gnuplot by
typing

gnuplot

and (inside gnuplot) write:

gnuplot> set xrange [-1.0:0.0]
gnuplot> set yrange [-1.0:0.0]
gnuplot> p ’SigW.txt’ u 1:(-0.37870+0.64690+$2+0.42476) w l, x w l

The first two lines set the range of the horizontal ω axis (x) and vertical axis to
be in a large region in which the quasiparticle energy is expected.
E
D
The third line plots the right hand side of Eq. 5, εKS + φ Σ̂X + Σ̂P (ε) − v̂xc φ =
(-0.37870+0.64690+$2+0.42476). Here, $2 stands for the second column of
the file, SigW_real= ΣP (ω). The remaining values in the brackets are the
Kohn-Sham orbital energy, exchange part of the self-energy, and negative of the
exchange-correlation potential energy, all given in the header.

The second plotted line is the frequency (its values are identical to the horizontal
x axis). The crossing of these two curves occurs at the quasiparticle energy, as
shown below.

6

22

http://www.gnuplot.info/

stochasticGW - User manual and Tutorial

By zooming in we see that the intersection occurs at ε = −0.582 Hartree =
−15.8 eV, i.e., close to the value found in the SigW.txt file.

! → If the Kohn-Sham eigenvalue is very far from the predicted quasiparticle energy
or the self-energy is very oscillatory, then check the results graphically, since there
could be several graphical solutions to Eq. 5.
The calculation should now be repeated with a higher value of nmctot to see how
the statistical error decreases with the number of MC samples.
In the production runs, carefully check the value of the energy broadening parameter gamma.
8.2

C60 : Fully stochastic calculation
Here we discuss a fully stochastic calculation of the electron affinity of C60 ; this
system is already big enough to profit from the stochastic treatment, since using
all the 120 occupied states for the time propagation is tedious and won’t bring
a large speedup over conventional approaches. We also illustrate how the timing and accuracy depends on the number of stochastic states (ntddft), on the
energy broadening parameter (gamma), and on the nxi and segment_fraction
parameters,
The user should first finish the previous example and get familiar with preparation
of the input files and the stochasticGW output.
We use a grid of 803 points, with a spacing of 0.4 Bohr. The LDA estimate of
electron affinity is 0.1645 Hartree, i.e., 4.48 eV, which is too high compared with
the experimental value of 2.69 eV.7
You should copy the files in the example directory ∼/sGW_v1/examples/C60
to your desired directory
mkdir ~/mygwrun
cd ~/mygwrun
cp ~/sGW_v1/examples/C60/* .
cp ~/sGW_v1/src/sgw.x .
7

https://aip.scitation.org/doi/abs/10.1063/1.4881421?journalCode=jcp

stochasticGW - User manual and Tutorial

23

cp -r ~/sGW_v1/PP/ .
The INPUT file is now:
pps
06-C.LDA.fhi
#
scratch_path
.
ekcut
20d0
nmctot
840
ntddft
8
units_nuc
bohr
gamma
0.06
binary
T
orb_indx
121
flg_dyn
T
box
sym

Note that for the LUMO energy we use a TDDFT effective potential instead of
an RPA (i.e., we set flg_dyn to be T ). As explained in Ref. 1, this yields a
much better LUMO energy.
In addition, the wf.bin waveufunction file for C60 is quite large (almost 1 gigabyte) so it is stored separately and you should download a zipped version of it
from
http://www.chem.ucla.edu/dept/Faculty/dxn/pdf/wf.bin.gz
and then write
gunzip wf.bin.gz
The calculation should be again done via mpi, e.g.,
mpirun -n 840 ./sgw.x >& log &
or, if one wants to use a lower number of cores and get the same exact numbers,
use, e.g.,
mpirun -n 120 ./sgw.x >& log &
albeit the calculation will be then approximately 840/120=7 times slower.
At the end of the calculation, the output file shows a quasiparticle energy (real
and imaginary parts) of
840

-0.096131

0.001864

MC, E_qp

+-stat.err

so the real part of the GW quasi-partilce energy is −0.096131 a.u.=-2.62eV with
a statistical error of 0.00186 a.u.=0.05 eV. This value is very close to the 2.69
eV experimental electron affinity of a single C60 .
24

stochasticGW - User manual and Tutorial

When run on a modern supercomputer, the wall time of the calculation (read in
the last few lines of the log file) was only 2994 sec=0.83 hr per each iteration.
For verification, we also ran, beyond this first run, several runs with different parameters: an increased number of stochastic TDDFT states ntddft=16, a reduced
energy width γ = 0.04 and different combination of the stochastic fragment parameters, nxi=5000-20,000 and segment_fraction=0.001, 0.003. There was
essentially no change in the results, and the real-part of the quasiparticle energy
changed by only 0.03 eV, i.e., less than the statistical error.
Finally, we also ran the RPA simulations (by removing the flg_dyn parameter,
which defaults to .false.), leading to a much deeper LUMO,
840

-0.124621

0.001648

MC, E_qp

+-stat.err

i.e., −0.124621 a.u.=-3.39eV.

stochasticGW - User manual and Tutorial

25

9

Acknowledgements
We are grateful for support by multiple sources.
This work was supported by the Center for Computational Study of Excited-State
Phenomena in Energy Materials at the Lawrence Berkeley National Laboratory,
which is funded by the U.S. Department of Energy, Office of Science, Basic Energy
Sciences, Materials Sciences and Engineering Division under Contract No. DEAC02-05CH11231, as part of the Computational Materials Sciences Program.
The work was further supported by the National Science Foundation, grants
DMR-1611382 and CHE-1112500 of Daniel Neuhauser, and CHE-1465064 of Eran
Rabani.
The Israel Science Foundation supported the work through grant 189/14 to Roi
Baer and through a FIRST Program grant (1700/14) to Roi Baer and Eran Rabani. Roi Baer also acknowledges support by the Binational Science Foundation,
Grant 2015687.
The code was tested and optimized within the XSEDE8 computational project
TG-CHE170058.
The stochasticGW code was mostly self-written, but it has several publicly
available routines, primarily the KISS random number generator of George Marsaglia;
we use, and slightly modified, the implementation of Jean-Michel Brankart. The
stochasticGW package also includes several spline routines from Numerical
Recipes.9

8

J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, et al. Xsede: accelerating scientific discovery. Comput. Sci.
Eng. 16,62 (2014). https://aip.scitation.org/doi/abs/10.1109/MCSE.2014.80
9 http://www.nr.com/

26

stochasticGW - User manual and Tutorial



---

Source Exif Data:

File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.5
Linearized                      : No
Page Count                      : 26
Page Mode                       : UseOutlines
Author                          : 
Title                           : 
Subject                         : 
Creator                         : LaTeX with hyperref package
Producer                        : pdfTeX-1.40.17
Create Date                     : 2018:05:16 19:21:23Z
Modify Date                     : 2018:05:16 19:21:23Z
Trapped                         : False
PTEX Fullbanner                 : This is pdfTeX, Version 3.14159265-2.6-1.40.17 (TeX Live 2016/Debian) kpathsea version 6.2.2

EXIF Metadata provided by EXIF.tools