atom.rst 35.3 KB
Newer Older
1
2
ATOM User manual
****************
Alberto Garcia's avatar
Alberto Garcia committed
3

4
5
:author: Alberto Garcı́a, ICMAB-CSIC, Barcelona
         albertog@icmab.es
Alberto Garcia's avatar
Alberto Garcia committed
6
7
8
9

PREFACE
=======

Alberto Garcia's avatar
Alberto Garcia committed
10
11
12
13
14
15
ATOM is the name of a program originally written (circa 1982) by
Sverre Froyen at the University of California at Berkeley, modified
starting in 1990 by Norman Troullier and Jose Luis Martins at the
University of Minnesota, and currently maintained by Alberto Garcia,
who added some features and made substantial structural changes to the
April 1990 (5.0) Minnesota version while at Berkeley and elsewhere.
Alberto Garcia's avatar
Alberto Garcia committed
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

Jose Luis Martins is maintaining his own version of the code:
::

   http://bohr.inesc-mn.pt/~jlm/pseudo.html

The program’s basic capabilities are:

-  All-electron DFT atomic calculations for arbitrary electronic
   configurations.

-  Generation of ab-initio pseudopotentials (several flavors).

-  Atomic calculations in which the effect of the core is represented by
   a previously generated pseudopotential. These are useful to make sure
   that the pseudopotential correctly reproduces the all-electron
   results for the valence complex.

A PRIMER ON AB-INITIO PSEUDOPOTENTIALS
======================================

37
38
In this case more than ever, there is a lot to be gained from reading
the original literature... Here are some basic references:
Alberto Garcia's avatar
Alberto Garcia committed
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995

-  Original idea of the ab-initio pseudopotential:

   | Kerker, J. Phys. C 13, L189-94 (1980)
   | Hamann, Schluter, Chiang, Phys. Rev. Lett. 43, 1494 (1979)

-  More on HSC scheme:

   | Bachelet, Schluter, Phys. Rev. B 25, 2103 (1982)
   | Bachelet, Hamann, Schluter, Phys. Rev. B 26, 4199 (1982)

-  Troullier-Martins elaboration of Kerker method:

   | Troullier, Martins, Phys. Rev. B 43, 1993 (1991)
   | Troullier, Martins, Phys. Rev. B 43, 8861 (1991)

-  Core corrections:

   Louie, Froyen, Cohen, Phys. Rev. B 26, 1738 (1982)

-  The full picture of plane-wave pseudopotential ab-initio
   calculations:

   W. E. Pickett, “Pseudopotential Methods in Condensed Matter
   Applications”, Computer Physics Reports 9, 115 (1989)

   M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D.
   Joannopoulos, “Iterative minimization techniques for ab initio
   total-energy calculations: molecular dynamics and conjugate
   gradients”, Rev. Mod. Phys. 64, 1045, (1992)

   Also, the book by Richard Martin “Electronic Structure: Basic Theory
   and Practical Methods” (Cambridge University Press) has a chapter on
   pseudopotentials.

-  Use in SIESTA:

   J.M. Soler, E. Artacho, J.D. Gale, A. Garcia, J. Junquera, P.
   Ordejon, D. Sanchez-Portal, “The SIESTA method for ab initio O(N)
   materials simulation”, Jour. Phys.: Condens. Matter, 14, 2745-2779
   (2002).

COMPILING THE PROGRAM
=====================

(Note that ATOM now depends on the GridXC and xmlf90 libraries for its
correct compilation. Please follow the instructions in the
``000_INSTALL`` file.)

Directory ``Tutorial`` in the source distribution contains a set of
scripts to automate the process of running ATOM and to analyze the
results. The file-manipulation details involved in each of the basic
functions of all-electron calculations, generation of pseudopotentials,
and testing of the pseudopotentias, are taken care of by ``ae.sh``,
``pg.sh``, and ``pt.sh``, respectively, all in the ``Tutorial/Utils``
directory. These scripts need to know where the ATOM executable ``atm``
is. If you have moved the ``Tutorial`` directory around, or you do not
have the source, the default location might not be right for you. The
easiest way to fix it is to define an environmental variable
``ATOM_PROGRAM``. Assuming ``atm`` is in ``/somedir/somewhere``, you
would do:

::

   ATOM_PROGRAM=/somedir/somewhere/atm ; export ATOM_PROGRAM  # sh-derived shells
   setenv ATOM_PROGRAM /somedir/somewhere/atm                 # csh-derived shells

Due to the shortcommings of the basic (GNUplot) plotting package used in
the Tutorial section, it is also necessary to copy some scripts from a
central repository. Again, if the default does not work for you, define
the ``ATOM_UTILS_DIR`` variable:

::

   ATOM_UTILS_DIR=/somewhere ; export ATOM_UTILS_DIR  # sh-derived shells
   setenv ATOM_UTILS_DIR /somewhere                   # csh-derived shells

USING THE ATOM PROGRAM
======================

All-electron calculations
-------------------------

Assume we want to find the orbital eigenvalues, total energy, and/or
charge density of Si in its ground state. (You should now go to the
``Tutorial/All_electron`` directory and try the following.) Our input
file is named ``si.ae.inp`` and contains the lines (see `Input File description`_
for more details):

::

      ae Si ground state all-electron
      Si   ca
          0.0
       3    2
       3    0      2.00      0.00
       3    1      2.00      0.00

   #2345678901234567890123456789012345678901234567890      Ruler

We can run the calculation by using the ``ae.sh`` script. Following the
layout of the ``Tutorial`` directory, we will assume that the script is
in the ``Tutorial/Utils`` directory. We run the script and go into the
directory created for the calculation (named as the input file without
the extension ``.inp``):

::

   $ sh ../Utils/ae.sh si.ae.inp
   ==> Output data in directory si.ae
   $ cd si.ae
   $ ls
   AECHARGE  CHARGE  RHO       charge.gplot   vcharge.gps
   AEWFNR0   INP     ae.gplot  charge.gps     vspin.gplot
   AEWFNR1   OUT     ae.gps    vcharge.gplot  vspin.gps
   $

We see some data files (those in all caps) and a few GNUPLOT plotting
scripts [1]_ .

The files are:

-  ``INP:`` A copy of the input file for the calculation.

-  ``OUT:`` Contains detailed information about the run.

-  ``AECHARGE:`` Contains in four columns values of :math:`r`, the “up”
   and “down” parts of the *total* charge density, and the total core
   charge density (the charges multiplied by :math:`4\pi r^2`).
   ``CHARGE`` is exactly identical and is generated for backwards
   compatibility.

-  ``RHO:`` Like ``CHARGE``, but without the :math:`4\pi r^2` factor.

-  ``AEWFNR0...AEWFNR3:`` All-electron valence wavefunctions as function
   of radius, for :math:`s`, :math:`p`, :math:`d`, and :math:`f` valence
   orbitals (0, 1, 2, 3, respectively — some channels might not be
   available). They include a factor of :math:`r`, the :math:`s`
   orbitals also going to zero at the origin.

It is interesting to peruse the ``OUT`` file. In particular, it lists
the orbital eigenvalues (in Rydbergs, as every other energy in the
program):

::

    nl    s      occ         eigenvalue    kinetic energy      pot energy

    1s   0.0    2.0000    -130.36911241     183.01377616    -378.73491463
    2s   0.0    2.0000     -10.14892694      25.89954259     -71.62102169
    2p   0.0    6.0000      -7.02876268      24.42537874     -68.74331203
    3s   0.0    2.0000      -0.79662742       3.23745215     -17.68692611
    3p   0.0    2.0000      -0.30705179       2.06135782     -13.62572515

(For a relativistic or spin-polarized calculation, there would be “up”
and “down” flags in the ``s`` column above.)

The **plotting** scripts come in two flavors: ``.gplot`` for terminal
use (default X11, use ``gnuplot -persist``), and ``.gps`` for postscript
output.

For all-electron calculations, the relevant scripts (without ``.gplot``
or ``.gps`` extensions) are:

-  ``charge:`` Charge density (separated core and valence contributions,
   multiplied by :math:`4\pi r^2`).

-  ``vcharge:`` Valence charge density (same normalization).

-  ``ae:`` Orbital valence wavefunctions (radial part multiplied by
   :math:`r`)

Pseudopotential generation
--------------------------

(You should now go to the ``Tutorial/Si`` directory and try the
following.) We are going to generate a pseudopotential for Si, using the
Troullier-Martins scheme. The calculation is relativistic and we use the
LDA (Ceperley-Alder flavor). The input file is named ``Si.tm2.inp`` and
contains the lines (see `Input File description`_ for more details):

::

   #
   #  Pseudopotential generation for Silicon
   #  pg: simple generation
   #
      pg      Silicon
           tm2      3.0             # PS flavor, logder R
    n=Si c=car                      # Symbol, XC flavor,{ |r|s}
          0.0       0.0       0.0       0.0       0.0       0.0
       3    4                       # norbs_core, norbs_valence
       3    0      2.00      0.00   # 3s2
       3    1      2.00      0.00   # 3p2
       3    2      0.00      0.00   # 3d0
       4    3      0.00      0.00   # 4f0
         1.90      1.90      1.90      1.90      0.00      0.00
   #
   # Last line (above): 
   #    rc(s)     rc(p)     rc(d)     rc(f)   rcore_flag  rcore
   #
   #23456789012345678901234567890123456789012345678901234567890

Note the two extra lines with respect to an all-electron calculation.
The pseudopotential core radii for all channels are 1.90 bohrs. Even
though they are nominally empty in the ground state, we include the
:math:`3d` and :math:`4f` states in order to generate the corresponding
pseudopotentials.

We can run the calculation by using the ``pg.sh`` script. Following the
layout of the ``Tutorial`` directory, we will assume that the script is
in the ``Tutorial/Utils`` directory. We run the script and go into the
directory created for the calculation (named as the input file without
the extension ``.inp``):

::

   $ sh ../../Utils/pg.sh Si.tm2.inp
   ==> Output data in directory Si.tm2
   ==> Pseudopotential in Si.tm2.vps and Si.tm2.psf
   $ cd Si.tm2
   $ ls [A-Z]*    # show only the data filesAE
   CHARGE  AEWFNR3   PSLOGD3  PSPOTR3  PSWFNR3     
   AELOGD0   CHARGE    PSPOTQ0  PSWFNQ0  RHO         
   AELOGD1   INP       PSPOTQ1  PSWFNQ1  SCRPSPOTR0  
   AELOGD2   OUT       PSPOTQ2  PSWFNQ2  SCRPSPOTR1  
   AELOGD3   PSCHARGE  PSPOTQ3  PSWFNQ3  SCRPSPOTR2  
   AEWFNR0   PSLOGD0   PSPOTR0  PSWFNR0  SCRPSPOTR3  
   AEWFNR1   PSLOGD1   PSPOTR1  PSWFNR1  VPSFMT      
   AEWFNR2   PSLOGD2   PSPOTR2  PSWFNR2  VPSOUT      

There are quite a few data files now. The new ones are:

-  ``PSCHARGE:`` Contains in four columns values of :math:`r`, the “up”
   and “down” parts of the *pseudo valence* charge density, and the
   pseudo core charge density (see Sect. `4.2.1 <#sec:cc>`__) (the
   charges multiplied by :math:`4\pi r^2`).

-  ``PSWFNR0...PSWFNR3:`` Valence pseudowavefunctions as function of
   radius, for :math:`s`, :math:`p`, :math:`d`, and :math:`f` valence
   orbitals (0, 1, 2, 3, respectively — some channels might not be
   available). They include a factor of :math:`r`, the :math:`s`
   orbitals also going to zero at the origin.

-  ``PSPOTR0...PSPOTR3:`` Ionic pseudopotentials (i.e. unscreened) as a
   function of :math:`r`, for :math:`s`, :math:`p`, :math:`d`, and
   :math:`f` channels (0, 1, 2, 3, respectively — some channels might
   not be available). The last column is :math:`-2Z_{ps}/r`, that is,
   the Coulomb potential of the pseudo atom. All the ionic
   pseudopotentials tend to this Coulomb tail for :math:`r` beyond the
   range of the core electrons.

-  ``PSPOTQ0...PSPOTQ3:`` Fourier transform :math:`V(q)` (times
   :math:`q^2/Z_{ps}`) of the ionic pseudopotentials as a function of
   :math:`q` (in bohr\ :math:`^{-1}`), for :math:`s`, :math:`p`,
   :math:`d`, and :math:`f` channels (0, 1, 2, 3, respectively — some
   channels might not be available).

-  ``PSWFNQ0...PSWFNQ3:`` Fourier transform :math:`\Psi(q)` of the
   valence pseudowavefunctions as a function of :math:`q` (in
   bohr\ :math:`^{-1}`), for :math:`s`, :math:`p`, :math:`d`, and
   :math:`f` channels (0, 1, 2, 3, respectively — some channels might
   not be available).

-  ``VPSOUT, VPSFMT:`` Files (formatted and unformatted) containing
   pseudopotential information. They are used for *ab-initio* codes such
   as SIESTA and PW. Copies of these files are deposited in the top
   directory after the run.

The ``OUT`` file has two sections, one for the all-electron (AE) run,
and another for the pseudopotential (PS) generation itself. It is
instructive to compare the AE and PS eigenvalues. Simply do

::

   $ grep '&v' OUT
    ATM3      12-JUL-02        Silicon
    3s   0.5    2.0000      -0.79937161       0.00000000     -17.74263363
    3p  -0.5    0.6667      -0.30807129       0.00000000     -13.66178958
    3p   0.5    1.3333      -0.30567134       0.00000000     -13.60785822
    3d  -0.5    0.0000       0.00000000       0.00000000      -0.27407047
    3d   0.5    0.0000       0.00000000       0.00000000      -0.27407047
    4f  -0.5    0.0000       0.00000000       0.00000000      -0.26482365
    4f   0.5    0.0000       0.00000000       0.00000000      -0.26482365
   ---------------------------- &v
    3s   0.5    2.0000      -0.79936061       0.50555315      -3.74113059
    3p  -0.5    0.6667      -0.30804995       0.77243805      -3.26356669
    3p   0.5    1.3333      -0.30565760       0.76702460      -3.25197500
    3d  -0.5    0.0000       0.00000000       0.00140505      -0.07847269
    3d   0.5    0.0000       0.00000000       0.00140505      -0.07847269
    4f  -0.5    0.0000       0.00000000       0.00243411      -0.07586534
    4f   0.5    0.0000       0.00000000       0.00243411      -0.07586534
   ---------------------------- &v

(The AE and PS eigenvalues are not *exactly* identical because the
pseudopotentials are changed slightly to make them approach their limit
tails faster).

The relevant plotting scripts (without ``.gplot`` or ``.gps``
extensions) are:

-  ``charge:`` It compares the AE and PS charge densities.

-  ``pseudo:`` A multi-page plot showing, on one page/window per
   channel:

   -  The AE and PS wavefunctions

   -  The AE and PS logarithmic derivatives.

   -  The real-space pseudopotential

   -  The Fourier-transformed pseudopotential (times :math:`q^2/Z_{ps}`)

-  ``pots:`` All the real-space pseudopotentials.

.. _sec:cc:

Core Corrections
~~~~~~~~~~~~~~~~

The program can generate pseudopotentials with the non-linear
exchange-correlation correction proposed in S. G. Louie, S. Froyen, and
M. L. Cohen, Phys. Rev. B 26, 1738 (1982).

In the traditional approach (which is the default for LDA calculations),
the pseudocore charge density equals the charge density outside a given
radius :math:`r_{pc}`, and has the smooth form

.. math:: \rho_{pc}(r) = A r   \sin(b r)

inside that radius. A smooth matching is provided with suitable
:math:`A` and :math:`b` parameters calculated by the program.

A new scheme has been implemented to fix some problems in the generation
of GGA pseudopotentials. The smooth function is now

.. math:: \rho_{pc}(r) =  r^2  \exp{(a + b r^2 +c r^4)}

and derivatives up to the second are continuous at :math:`r_{pc}`.

To use core corrections in the pseudopotential generation the jobcode in
the first line should be ``pe`` instead of ``pg``.

The radius :math:`r_{pc}` should be given in the sixth slot in the last
input line (see above). If it is negative or zero (or blank), the radius
is then computed using the fifth number in that line (``rcore_flag``,
see the example input file above) and the following criterion: at
:math:`r_{pc}` the core charge density equals ``rcore_flag``\ \*(valence
charge density). It is *highly recommended* to set an explicit value for
the pseudocore radius :math:`r_{pc}`, rather than letting the program
provide a default.

If ``rcore_flag`` is input as negative, the full core charge is used. If
``rcore_flag`` is input as zero, it is set equal to one, which will be
thus the default if ``pe`` is given but no numbers are given for these
two variables.

The output file contains the radius used and the :math:`A` (:math:`a`)
and :math:`b` (and :math:`c`) parameters used for the matching. The
``VPSOUT`` and ``VPSFMT`` files will contain the pseudocore charge in
addition to the pseudopotential.

It is possible to override the default (new scheme for GGA calculations,
old scheme for LDA calculations) by using the directives

::

   %define NEW_CC
   %define OLD_CC

The program will issue the appropriate warnings. (See
the  `Input File description`_)

Relevant files:

-  ``PSCHARGE:`` Contains the pseudocore charge in column four.
   (multiplied by :math:`4\pi r^2`).

-  ``COREQ:`` Fourier transform of the pseudocore charge density
   :math:`\rho_{pc}(q)` in units of electrons, with :math:`q` in
   bohr\ :math:`^{-1}`.

Useful plotting scripts (without ``.gplot`` or ``.gps`` extensions) are:

-  ``charge:`` Shows also the pseudocore charge.

-  ``coreq:`` Shows the Fourier transform of the pseudocore charge.

Pseudopotential test
--------------------

While it is helpful to “have a look” at the plots of the pseudopotential
generation to get a feeling for its quality, there is no substitute for
a proper **transferability testing**. A pseudopotential with good
transferability will reproduce the all-electron energy levels and
wavefunctions in arbitrary environments, (i.e., in the presence of
bonding, which always takes place when forming solids and molecules). We
know that norm conservation guarantees a certain degree of
transferability (usually seen clearly in the plot of the logarithmic
derivative), but we can get a better assessment by performing
all-electron and “pseudo” calculations on the same series of atomic
configurations, and comparing the eigenvalues and excitation energies.

In the same ``Tutorial/Si`` directory we can find file ``Si.test.inp``,
containing the concatenation of ten jobs. The first five are
all-electron (``ae``) calculations, and the last five, pseudopotential
test (``pt``) runs for the same configurations:

::

   #
   # All-electron calculations for a series of Si configurations
   #
      ae Si Test -- GS 3s2 3p2
      Si   ca
          0.0
       3    2
       3    0      2.00
       3    1      2.00
      ae Si Test -- 3s2 3p1 3d1
      Si   ca
          0.0
       3    3
       3    0      2.00
       3    1      1.00
       3    2      1.00
      ae Si Test -- 3s1 3p3
      Si   ca
          0.0
       3    2
       3    0      1.00
       3    1      3.00
      ae Si Test -- 3s1 3p2 3d1
      Si   ca
          0.0
       3    3
       3    0      1.00
       3    1      2.00
       3    2      1.00
      ae Si Test -- 3s0 3p3 3d1
      Si   ca
          0.0
       3    3
       3    0      0.00
       3    1      3.00
       3    2      1.00
   #
   # Pseudopotential test calculations
   #
      pt Si Test -- GS 3s2 3p2
      Si   ca
          0.0
       3    2
       3    0      2.00
       3    1      2.00
      pt Si Test -- 3s2 3p1 3d1
      Si   ca
          0.0
       3    3
       3    0      2.00
       3    1      1.00
       3    2      1.00
      pt Si Test -- 3s1 3p3
      Si   ca
          0.0
       3    2
       3    0      1.00
       3    1      3.00
      pt Si Test -- 3s1 3p2 3d1
      Si   ca
          0.0
       3    3
       3    0      1.00
       3    1      2.00
       3    2      1.00
      pt Si Test -- 3s0 3p3 3d1
      Si   ca
          0.0
       3    3
       3    0      0.00
       3    1      3.00
       3    2      1.00

The configurations differ in the promotion of electrons from one level
to another (it is also possible to transfer *fractions* of an electron).

We can run the file by using the ``pt.sh`` script. Following the layout
of the ``Tutorial`` directory, we will assume that the script is in the
directory directly above the current one. We need to give it **two**
arguments: the calculation input file, and the file containing the
pseudopotential we want to test. Let’s make the latter ``Si.tm2.vps``:

::

   $ sh ../../Utils/pt.sh Si.test.inp Si.tm2.vps
   ==> Output data in directory Si.test-Si.tm2
   $ cd Si.test-Si.tm2/
   $ ls [A-Z]*
   AECHARGE  AEWFNR1  CHARGE  OUT       PTWFNR0  PTWFNR2  VPSIN
   AEWFNR0   AEWFNR2  INP     PTCHARGE  PTWFNR1  RHO

The working directory is named after *both* the test and pseudopotential
files. It contains several new files:

-  ``VPSIN:`` A copy of the pseudopotential file to be tested.

-  ``PTCHARGE:`` Contains in four columns values of :math:`r`, the “up”
   and “down” parts of the *pseudo valence* charge density, and the
   pseudo core charge density (see Sect. `4.2.1 <#sec:cc>`__) (the
   charges multiplied by :math:`4\pi r^2`).

-  ``PTWFNR0...PTWFNR3:`` Valence pseudowavefunctions as function of
   radius, for :math:`s`, :math:`p`, :math:`d`, and :math:`f` valence
   orbitals (0, 1, 2, 3, respectively — some channels might not be
   available). They include a factor of :math:`r`, the :math:`s`
   orbitals also going to zero at the origin.

The ``OUT`` file has two sections, one for the all-electron (AE) runs,
and another for the pseudopotential tests (PT). At the end of each
series of runs there is a table showing the excitation energies. A handy
way to compare the AE and PT energies is:

::

   $ grep '&d' OUT
   [...elided...]
    &d total energy differences in series
    &d          1        2        3        4        5
    &d  1    0.0000
    &d  2    0.4308   0.0000
    &d  3    0.4961   0.0653   0.0000
    &d  4    0.9613   0.5305   0.4652   0.0000
    &d  5    1.4997   1.0689   1.0036   0.5384   0.0000
   *----- End of series ----* spdfg &d&v
    ATM3      12-JUL-02   Si Test -- GS 3s2 3p2
    ATM3      12-JUL-02   Si Test -- 3s2 3p1 3d1
    ATM3      12-JUL-02   Si Test -- 3s1 3p3
    ATM3      12-JUL-02   Si Test -- 3s1 3p2 3d1
    ATM3      12-JUL-02   Si Test -- 3s0 3p3 3d1
    &d total energy differences in series
    &d          1        2        3        4        5
    &d  1    0.0000
    &d  2    0.4299   0.0000
    &d  3    0.4993   0.0694   0.0000
    &d  4    0.9635   0.5336   0.4642   0.0000
    &d  5    1.5044   1.0745   1.0051   0.5409   0.0000
   *----- End of series ----* spdfg &d&v

The tables (top AE, bottom PT) give the cross-excitations among all
configurations. Typically, one should be all right if the AE-PT
differences are not much larger than 1 mRy.

You can also compare the AE and PT eigenvalues. Simply do

::

   $ grep '&v' OUT | grep s
    ATM3      12-JUL-02   Si Test -- GS 3s2 3p2
    3s   0.0    2.0000      -0.79662742       3.23745215     -17.68692611
    ATM3      12-JUL-02   Si Test -- 3s2 3p1 3d1
    3s   0.0    2.0000      -1.08185979       3.53885995     -18.40569836
    ATM3      12-JUL-02   Si Test -- 3s1 3p3
    3s   0.0    1.0000      -0.85138783       3.35438895     -17.96219240
    ATM3      12-JUL-02   Si Test -- 3s1 3p2 3d1
    3s   0.0    1.0000      -1.11431855       3.62997498     -18.60814708
    ATM3      12-JUL-02   Si Test -- 3s0 3p3 3d1
    3s   0.0    0.0000      -1.14358268       3.71462770     -18.79448684
   *----- End of series ----* spdfg &d&v
    ATM3      12-JUL-02   Si Test -- GS 3s2 3p2
    1s   0.0    2.0000      -0.79938037       0.50556261      -3.74114712
    ATM3      12-JUL-02   Si Test -- 3s2 3p1 3d1
    1s   0.0    2.0000      -1.08384468       0.55070398      -3.81988817
    ATM3      12-JUL-02   Si Test -- 3s1 3p3
    1s   0.0    1.0000      -0.85392666       0.52020429      -3.76852577
    ATM3      12-JUL-02   Si Test -- 3s1 3p2 3d1
    1s   0.0    1.0000      -1.11546463       0.56048425      -3.83646615
    ATM3      12-JUL-02   Si Test -- 3s0 3p3 3d1
    1s   0.0    0.0000      -1.14353959       0.56945741      -3.85106049
   *----- End of series ----* spdfg &d&v

(and similarly for :math:`p`, :math:`d`, and :math:`f`, if desired).
Again, the typical difference should be of around 1 mRyd for a “good”
pseudopotential. (The *real* proof of good transferability, remember,
can only come from a molecular or solid-state calculation). Note that
the PT levels (in older versions of ATOM) are labeled starting from
principal quantum number 1.

The relevant plotting scripts (without ``.gplot`` or ``.gps``
extensions) are:

-  ``charge:`` It compares the AE and PT charge densities.

-  ``pt:`` Compares the valence all-electron and pseudo-wavefunctions.

.. _Input File description:

APPENDIX: THE INPUT FILE
========================

For historical reasons, the input file is in a rigid column format.
Fortunately, most of the column fields line up, so the possibility of
errors is reduced. We will begin by describing in detail a very simple
input file for an **all-electron calculation** for the ground state of
Si. More examples can be found in the ``Tutorial`` directory.

The file itself is:

::

   #
   #  Comments allowed here
   #
      ae Si ground state all-electron
      Si   car  
          0.0
       3    2
       3    0      2.00      0.00
       3    1      2.00      0.00
   #
   # Comments allowed here
   #
   #2345678901234567890123456789012345678901234567890      Ruler

-  The first line specifies:

   -  The calculation code (``ae`` here stands for “all-electron”).

   -  A title for the job (here ``Si ground state all-electron``).

   (format 3x,a2,a50)

-  Second line:

   -  Chemical symbol of the nucleus (here ``Si``, obviously)

   -  Exchange-correlation type. Here, ``ca`` stands for Ceperley-Alder.
      The options are:

      Local density approximations (LDA):

      -  ``wi``: Wigner, PR 46, 1002 (1934)

      -  ``hl``: Hedin and Lundqvist, J. Phys. C 4, 2064 (1971)

      -  ``gl``: Gunnarson and Lundqvist, PRB 13, 4274 (1976)

      -  ``bh``: von Barth and Hedin, J. Phys. C 5, 1629 (1972)

      -  ``ca``: Ceperley-Alder, parametrized by Perdew and Zunger, PRB
         23, 5075 (1981)

      -  ``pw``: PW92 (Perdew and Wang, PRB, 45, 13244 (1992))
         (recommended LDA functional)

      Generalized gradient approximations (GGA), implemented as Balbás,
      Martins, and Soler, PRB 64, 165110 (2001):

      -  ``pb``: PBE (Perdew, Burke, and Ernzerhof, PRL 77, 3865 (1996))
         (recommended GGA)

      -  ``wp``: PW91 (Perdew and Wang, JCP, 100, 1290 (1994))

      -  ``rp``: RPBE (Hammer, Hansen, and Norskov, PRB 59, 7413 (1999))

      -  ``rv``: revPBE (Zhang and Yang, PRL 80, 890 (1998))

      -  ``ps``: PBEsol (Perdew et al, PRL 100, 136406 (2008))

      -  ``wc``: WC (Wu and Cohen, PRB 73, 235116 (2006))

      -  ``jo``: PBEJsJrLO. This and the next three functionals are
         reparametrizations of the PBE functional by Pedroza et al, PRB
         79, 201106 (2009) and Odashima et al, J. Chem. Theory Comp. 5,
         798 (2009). Js and Jr refer to jellium surface energy and
         jellium response, respectively. LO refers to Lieb-Oxford bound.
         Gx and Gc refer to gradient expansions for exchange and
         correlation, respectively. HGE refers to the low-density limit
         of the homogeneous electron gas.

      -  ``jh``: PBEJsJrHEG

      -  ``go``: PBEGcGxLO

      -  ``gh``: PBEGcGxHEG

      -  ``am``: AM05 (Armiento and Mattsson PRB 72, 085108 (2005), PRB
         79, 155101 (2009)),

      -  ``bl``: BLYP (Becke, PRA 38, 3098 (1988) and Lee, Yang, and
         Parr, PRB 37, 785 (1988))

      Van der Waals (VDW) density functionals, implemented as
      Román-Pérez and Soler, PRL 103, 096102 (2009):

      -  ``vw`` or ``vf``: DRSLL (Dion et al, PRL 92, 246401 (2004))

      -  ``vl``: LMKLL (Lee et al, PRB 82, 081101 (2010))

      -  ``vk``: KBM (Klimes, Bowler, and Michaelides, JPCM 22, 022201
         (2009))

      -  ``vc``: C09 (Cooper, PRB 81, 161104 (2010))

      -  ``vb``: BH (Berland and Hyldgaard, PRB 89, 035412 (2014))

      -  ``vv``: VV (Vydrov and VanVoorhis, JCP 133, 244103 (2010))

      Support for libxc functionals (if compiled in in libGridXC) is
      available through the use of the special code ``xc``, and a
      composite integer of the form 0XXX0YYY that follows later in the
      line. For example:

      ::

         #  Example with libxc functional
            pg      Silicon
                 tm2      3.0             # PS flavor, logder R
          n=Si c=xcr 01010130      # Symbol, XC flavor,{ |r|s} {Libxc code}
                0.0       0.0       0.0       0.0       0.0       0.0
             3    4                       # norbs_core, norbs_valence
             3    0      2.00      0.00   # 3s2
             3    1      2.00      0.00   # 3p2
             3    2      0.00      0.00   # 3d0
             4    3      0.00      0.00   # 4f0
               1.90      1.90      1.90      1.90      0.00      0.00
         #23456789012345678901234567890123456789012345678901234567890

      where the code “01010130” packs the two codes “101” and “130” for
      the PBE functionals.

   -  The character ``r`` next to ``ca`` is a flag to perform the
      calculation relativistically, that is, solving the Dirac equation
      instead of the Schrodinger equation. The full range of options is:

      -  ``s`` : Spin-polarized calculation, non-relativistic.

      -  ``r``: Relativistic calculation, obviously polarized.

      -  (blank) : Non-polarized (spin ignored), non-relativistic
         calculation.

   (format 3x,a2,3x,a2,a1,1x,i8)

-  Third line. Its use is somewhat esoteric and for most calculations it
   should contain just a 0.0 in the position shown, but that first field
   might be useful to generate pseudopotentials for “atoms” with a
   fractional atomic number (see the example for ON in the
   ``Tutorial/PS_Generation`` directory).

The rest of the file is devoted to the specification of the electronic
configuration:

-  | Fourth line:
   | Number of core and valence orbitals. For example, for Si, we have
     :math:`1s`, :math:`2s`, and :math:`2p` in the core (a total of 3
     orbitals), and :math:`3s` and :math:`3p` in the valence complex (2
     orbitals).

   (format 2i5)

-  Fifth, sixth... lines: (there is one line for each valence orbital)

   -  ``n`` (principal quantum number)

   -  ``l`` (angular momentum quantum number)

   -  Occupation of the orbital in electrons.

   (format 2i5,2f10.3)

   (There are two f input descriptors to allow the input of “up” and
   “down” occupations in spin-polarized calculations (see example
   below))

Comments or blank lines may appear in the file at the beginning and at
the end. It is possible to perform two or more calculations in
succession by simply concatenating blocks as the one described above.
For example, the following file is used to study the ground state of N
and an excited state with one electron promoted from the :math:`2s` to
the :math:`2p` orbital taking into account the spin polarization:

::

   #
      ae N ground state all-electron
      N    cas
          0.0
       1    2
       2    0      2.00      0.00
       2    1      3.00      0.00
   #
   #  Second calculation starts here
   #
      ae N 1s2 2s1 2p4  all-electron
      N    cas
          0.0
       1    2
       2    0      1.00      0.00
       2    1      3.00      1.00

   #2345678901234567890123456789012345678901234567890      Ruler

The different treatment of core and valence orbitals in the input for an
all-electron calculation is purely cosmetic. The program “knows” how to
fill the internal orbitals in the right order, so it is only necessary
to give their number. That is handy for heavy atoms... Overzealous users
might want to check the output to make sure that the core orbitals are
indeed correctly treated.

For a **pseudopotential test calculation**, the format is exactly the
same, except that the job code is ``pt`` instead of ``ae``.

For a **pseudopotential generation run**, in addition to the electronic
configuration chosen for the generation of the pseudopotentials (which
is input in the same manner as above), one has to specify the “flavor”
(generation scheme) and the set of core radii :math:`r_c` for the
construction of the pseudowavefunction. Here is an example for Si using
the Hamann-Schluter-Chiang scheme:

::

   # 
      pg Si Pseudopotencial
           hsc     2.00
      Si   ca
            0
       3    3
       3    0      2.00
       3    1      0.50
       3    2      0.50
         1.12      1.35      1.17       0.0       0.0       0.0
   #
   #23456789012345678901234567890123456789012345678901234567890   Ruler
   ---------------------------------------

Apart from the ``pg`` (pseudopotential generation) job code in the first
line, there are two extra lines:

-  Second line:

   Flavor and radius at which to compute logarithmic derivatives for
   test purposes.

   The flavor can be one of :

   === ==========================
   hsc Hamann-Schluter-Chiang
   ker Kerker
   tm2 Improved Troullier-Martins
   === ==========================

   The ``ker`` and ``tm2`` schemes can get away with larger :math:`r_c`,
   due to their wavefunction matching conditions.

   (format 8x, a3, f9.3)

-  The last line (before the blank line) specifies:

   -  The values of the :math:`r_c` in atomic units (bohrs) for the
      :math:`s`, :math:`p`, :math:`d`, and :math:`f` orbitals (it is a
      good practice to input the valence orbitals in the order of
      increasing angular momentum, so that there is no possible
      confusion).

      (format 4f10.5)

   -  Two extra fields (2f10.5) which are relevant only if non-local
      core corrections are used (see Sect `4.2.1 <#sec:cc>`__).

In the ``hsc`` example above, only :math:`s` ,\ :math:`p`, and :math:`d`
:math:`r_c`\ ’s are given. Here is an example for Silicon in which we
are only interested in the :math:`s` and :math:`p` channels for our
pseudopotential, and use the Kerker scheme:

::

   # 
      pg Si Kerker generation
           ker     2.00
      Si   ca
            0
       3    2
       3    0      2.00
       3    1      2.00
         1.80      1.80      0.00       0.0       0.0       0.0

   #23456789012345678901234567890123456789012345678901234567890   Ruler

This completes the discussion of the more common features of the input
file. See the Appendix `6 <#sec:directives>`__ for more advanced
options.

.. _sec:directives:

APPENDIX: INPUT FILE DIRECTIVES
===============================

The fixed format can be a source of desperation for the beginner, and
its rigidity means that it is not easy to add new items to the input.
For this purpose, the program takes another route: several variables can
be entered in a specially flexible format by means of *directives* at
the top of the file. For example

::

   %define NEW_CC
   .... rest of the input file

would signal that we want to use a new core-correction scheme.

There are two kinds of directives, with syntax:

::

   %VARIABLE=value
   %define NAME

In the first case we assign the value ``value`` to the variable
``VARIABLE``. The program can look at the value via a special subroutine
call.

The second form is a bit more abstract, but can be understood as
assigning a special “existence” value of ``1`` to the variable ``NAME``.
Again, the program can check for the existence of the variable via a
special subroutine call.

Currently, the program understands the following ``NAME``\ s:

-  ``COMPAT_UCB:`` Revert to the standard circa 1990 UCB values. Note
   that these correspond to the first released version of Jose Luis
   Martins code, not to the old Froyen version. (The defaults are: use a
   denser grid up to larger radii. Use a larger value for the
   pseudopotential cutoff point. Use the Soler-Balbas XC package.)

-  ``NEW_CC:`` New core-correction scheme

-  ``OLD_CC:`` Old core-correction scheme (see
   Sect. `4.2.1 <#sec:cc>`__)

-  ``USE_OLD_EXCORR:`` Use the old exchange-correlation package.

-  ``NO_PS_CUTOFFS:`` Avoid cutting off the tails of the
   pseudopotentials. Currently, a simple exponential tapering function
   is used, which introduces a discontinuity in the first derivative of
   the ionic pseudopotential.

-  ``FREE_FORMAT_RC_INPUT:`` Use free-format for the input of the cutoff
   radii and the specification of the core-correction parameters. This
   is useful for externally-driven runs of ATOM. In this case the user
   should make sure that all six values (four rc’s plus the two cc
   parameters) are present in the input line.

.. [1]
   GNUPLOT has its issues, but it is free, and installed almost
   everywhere. Hence we have chosen it as the lowest-common denominator
   for basic plotting