« back to all changes in this revision
Viewing changes to Src/basis_specs.f
- browse files at revision 591
- compare with another revision
- download diff
- download tarball
- view history from revision 591
- Committer: Alberto Garcia
- Date: 2018-06-28 12:07:22 UTC
- mfrom: (589.1.10 psml-nopols)
- Revision ID: albertog@icmab.es-20180628120722-p4tda7m283wl6xc2
Non-perturbative polarization option. Assorted clarifications
The default perturbative polarization scheme can be changed to another
in which polarization shells are promoted to 'first-class'
shells. This could be done explicitly in the PAO.Basis block, but now
it can be applied in cases where a compact (i.e. 'DZP') basis
specification is used, or even to cases in which a PAO.Basis block
with 'P' options is used.
While simple perturbative polarization seems to give better results
(in terms of final energy), one motivation for the new option has been
to cover cases in which the code mis-behaves. This happens
for (Hamann) pseudos with only s and p pseudopotential channels, and
with (at least) an 's' semicore shell.
This work has led to other clarifications throughout the basis-set
specification code. Among them:
- Set lmxo always to the maximum l of the basis orbitals
When perturbative polarization orbitals are present, the maximum l
reported for the basis set (lmxo) could sometimes be underestimated
(only the 'working' shells were taken into account). Now lmxo is set
always to the real maximum.
In case a highest-l shell is 'empty', it is reported as a
polarization orbital.
- Improve the reporting of polarization relationships.
- Improve the defaults for the 'nonodes' basis-set generation option.
- Fix the Util/Genbasis code.
The default perturbative polarization scheme can be changed to another
in which polarization shells are promoted to 'first-class'
shells. This could be done explicitly in the PAO.Basis block, but now
it can be applied in cases where a compact (i.e. 'DZP') basis
specification is used, or even to cases in which a PAO.Basis block
with 'P' options is used.
While simple perturbative polarization seems to give better results
(in terms of final energy), one motivation for the new option has been
to cover cases in which the code mis-behaves. This happens
for (Hamann) pseudos with only s and p pseudopotential channels, and
with (at least) an 's' semicore shell.
This work has led to other clarifications throughout the basis-set
specification code. Among them:
- Set lmxo always to the maximum l of the basis orbitals
When perturbative polarization orbitals are present, the maximum l
reported for the basis set (lmxo) could sometimes be underestimated
(only the 'working' shells were taken into account). Now lmxo is set
always to the real maximum.
In case a highest-l shell is 'empty', it is reported as a
polarization orbital.
- Improve the reporting of polarization relationships.
- Improve the defaults for the 'nonodes' basis-set generation option.
- Fix the Util/Genbasis code.
added
removed
Src/basis_specs.f
28
28
! * Read information about the valence charge density in the
29
29
! pseudopotential file and determine whether there are semicore
30
30
! states.
31
! (A further semicore analysis is 'autobasis')
31
! (A further semicore analysis is in 'autobasis')
32
32
! * Read the optional fdf blocks:
33
33
! AtomicMass (routine remass)
34
34
! PAO.BasisSize - Overrides the default 'basis_size' on a per-species basis.
35
! PAO.PolarizationScheme - Whether to use the standard perturbative approach
36
! for polarization orbitals, or to promote the shell
37
! to a stand-alone status.
35
38
! PAO.Basis - This is the most complex block, very flexible but in
36
39
! need of spelling-out the specific restrictions. Line by
37
40
! line, these are:
80
83
! There are no 'per l-shell' polarization orbitals, except if the
81
84
! third character of 'basis_size' is 'p' (as in 'dzp'), in which
82
85
! case polarization orbitals are defined so they have the minimum
83
! angular momentum l such that there are not occupied orbitals
86
! angular momentum l such that there are no occupied orbitals
84
87
! with the same l in the valence shell of the ground-state
85
88
! atomic configuration. They polarize the corresponding l-1 shell.
86
!
89
! (See above for generation scheme)
90
!
87
91
! The Soft-Confinement parameters 'rinn' and 'vcte' are set to 0.0
88
92
! The Charge-Confinement parameters 'qcoe', 'qyuk' and 'qwid'
89
93
! are set to 0.0, 0.0 and 0.01
142
146
143
147
type(basis_def_t), pointer :: basp => null()
144
148
type(shell_t), pointer :: s => null()
149
type(shell_t), pointer :: polarized_shell_ptr
145
150
type(lshell_t), pointer :: ls => null()
146
151
type(kbshell_t), pointer :: k => null()
147
152
388
393
389
394
call remass()
390
395
call resizes()
396
call read_polarization_scheme()
391
397
call repaobasis()
392
398
call autobasis()
393
399
call relmxkb()
666
672
. basp%ionic_charge = fdf_bvalues(pline,2)
667
673
allocate(basp%tmp_shell(basp%nshells_tmp))
668
674
675
! These are (l,n) shells
669
676
shells: do ish= 1, basp%nshells_tmp
670
677
s => basp%tmp_shell(ish)
671
678
call initialize(s)
728
735
else
729
736
s%nzeta_pol = 1
730
737
endif
738
basp%lmxo = max(basp%lmxo,s%l+1) ! NOTE new behavior
731
739
endif
732
740
!
733
741
! Soft-confinement
828
836
! Clean up for this species
829
837
enddo
830
838
!
831
! OK, now classify the states by l-shells
839
! OK, now classify the (l,n) states by l-shells
832
840
!
833
841
do isp = 1, nsp
834
842
basp => basis_parameters(isp)
835
843
if (basp%lmxo .eq. -1) cycle !! Species not in block
836
844
!!
845
polarized_shell_ptr => null()
837
846
allocate (basp%lshell(0:basp%lmxo))
847
838
848
loop_l: do l= 0, basp%lmxo
839
849
ls => basp%lshell(l)
840
850
call initialize(ls)
844
854
do ish= 1, basp%nshells_tmp
845
855
s => basp%tmp_shell(ish)
846
856
if (s%l .eq. l) nn=nn+1
857
858
if (basp%non_perturbative_polorbs) then
859
! Make room for a new one (explicit polarization orbital)
860
if (s%polarized .and. (s%l == (l-1))) nn=nn+1
861
endif
847
862
enddo
848
863
ls%nn = nn
849
864
if (nn.eq.0) then
856
871
ind = 0
857
872
do ish=1, basp%nshells_tmp
858
873
s => basp%tmp_shell(ish)
874
859
875
if (s%l .eq. l) then
860
876
ind = ind+1
861
877
call copy_shell(source=s,target=ls%shell(ind))
862
endif
863
enddo
864
if (nn.eq.1) then
878
if (basp%non_perturbative_polorbs) then
879
if (s%polarized) then
880
! Save for later. There should be just one 'polarized' shell
881
polarized_shell_ptr => ls%shell(ind)
882
! Remove markers (in new objects)
883
ls%shell(ind)%polarized = .false.
884
ls%shell(ind)%was_polarized = .true.
885
ls%shell(ind)%nzeta_pol = 0
886
endif
887
endif
888
endif
889
890
if (basp%non_perturbative_polorbs) then
891
if (s%polarized .and. (s%l == (l-1))) then
892
! Note that we have already seen this (parent) shell
893
! in the previous iteration of loop_l
894
ind = ind + 1
895
! Copy again to inherit data
896
call copy_shell(source=s,target=ls%shell(ind))
897
ls%shell(ind)%n = -1 ! reset to find later
898
ls%shell(ind)%l = l
899
ls%shell(ind)%nzeta = s%nzeta_pol
900
ls%shell(ind)%nzeta_pol = 0
901
ls%shell(ind)%polarized = .false.
902
ls%shell(ind)%polarization_shell = .true.
903
ls%shell(ind)%shell_being_polarized =>
904
$ polarized_shell_ptr
905
polarized_shell_ptr => null()
906
! rc's, lambdas, etc, will remain as in s.
907
endif
908
endif
909
910
enddo
911
912
if (nn.eq.1) then
865
913
! If n was not specified, set it to ground state n
866
914
if (ls%shell(1)%n.eq.-1)
867
915
. ls%shell(1)%n=basp%ground_state%n(l)
882
930
integer lmxkb
883
931
integer, intent(in) :: is
884
932
885
integer lpol, l, i
886
887
933
lmxkb = -1
888
934
basp=>basis_parameters(is)
889
935
if (basp%z .le. 0) return ! Leave it at -1 for floating orbs.
890
936
891
lmxkb = basp%lmxo + 1
892
893
! But watch out for polarization orbitals...
894
! We ASSUME that these do not count towards lshell%nn...
895
896
lpol = 0
897
do l = 0, basp%lmxo
898
ls=>basp%lshell(l)
899
do i = 1, ls%nn
900
s=>ls%shell(i)
901
if (s%polarized) lpol = s%l + 1
902
enddo
903
enddo
904
if (lpol .gt. basp%lmxo) lmxkb = lpol + 1
937
lmxkb = basp%lmxo + 1 ! This now includes any polarization orbitals
905
938
906
939
write(6,'(3a,i1,/,2a,/,a)') 'For ', trim(basp%label),
907
940
. ', standard SIESTA heuristics set lmxkb to ',
961
994
end subroutine resizes
962
995
963
996
!-----------------------------------------------------------------------
997
subroutine read_polarization_scheme()
998
999
! Optionally read polarization orbital scheme for different species.
1000
1001
type(block_fdf) :: bfdf
1002
type(parsed_line), pointer :: pline
1003
1004
character(len=40) :: string
1005
logical :: non_perturbative_pols
1006
integer isp
1007
1008
non_perturbative_pols =
1009
$ fdf_boolean('PAO.NonPerturbative.Polarization.Orbitals',
1010
$ .false.)
1011
1012
loop: do isp=1, nsp
1013
basp=>basis_parameters(isp)
1014
basp%non_perturbative_polorbs = non_perturbative_pols
1015
end do loop
1016
1017
1018
if (fdf_block('PAO.PolarizationScheme',bfdf)) then
1019
do while(fdf_bline(bfdf,pline))
1020
if (.not. fdf_bmatch(pline,'nn'))
1021
. call die("Wrong format in PAO.PolarizationScheme")
1022
isp = label2species(fdf_bnames(pline,1))
1023
if (isp .eq. 0) then
1024
write(6,'(a,1x,a)')
1025
. "WRONG species symbol in PAO.PolarizationScheme:",
1026
. trim(fdf_bnames(pline,1))
1027
call die("Wrong species in PAO.PolarizationScheme")
1028
else
1029
1030
basp => basis_parameters(isp)
1031
string = fdf_bnames(pline,2)
1032
1033
select case (trim(string))
1034
case ("non-perturbative")
1035
basp%non_perturbative_polorbs = .true.
1036
case ("perturbative")
1037
basp%non_perturbative_polorbs = .false.
1038
case default
1039
call die("Bad keyword in PAO.PolarizationScheme")
1040
end select
1041
1042
endif
1043
enddo
1044
endif
1045
1046
end subroutine read_polarization_scheme
1047
1048
!-----------------------------------------------------------------------
964
1049
965
1050
subroutine relmxkb()
966
1051
1157
1242
!
1158
1243
integer l, nzeta, nzeta_pol
1159
1244
integer :: nsemic_shells(0:3), nsh, i
1245
type(lshell_t), pointer :: ls_parent
1160
1246
1161
1247
loop: do isp=1, nsp
1162
1248
basp=>basis_parameters(isp)
1181
1267
if (nsemic_shells(l) > 0) basp%lmxo = max(l,basp%lmxo)
1182
1268
enddo
1183
1269
1184
!
1270
! Check whether we need to consider larger l's due to
1271
! polarization orbitals (this is a new behavior -- in the
1272
! case of perturbative pol orbs, the highest-l shell will remain empty)
1273
! But lmxo in the 'specs' output will be correct, and not confusing
1185
1274
1275
if (basp%basis_size(3:3) .eq. 'p') then
1276
loop_pol: do l=1,4 ! Note that l starts at 1
1277
if ( (.not. basp%ground_state%occupied(l)) ) then
1278
basp%lmxo = max(l,basp%lmxo)
1279
EXIT loop_pol
1280
endif
1281
enddo loop_pol
1282
endif
1283
1186
1284
allocate (basp%lshell(0:basp%lmxo))
1187
1285
1188
1286
if (basp%basis_size(1:2) .eq. 'sz') nzeta = 1
1209
1307
s%n = basp%ground_state%n(l) - (nsh-i)
1210
1308
s%nzeta = 1
1211
1309
s%polarized = .false.
1310
s%polarization_shell = .false.
1212
1311
s%split_norm = global_splnorm
1213
1312
s%nzeta_pol = 0
1214
1313
1242
1341
s%nzeta = 0
1243
1342
endif
1244
1343
s%polarized = .false.
1344
s%polarization_shell = .false.
1245
1345
if (abs(basp%z).eq.1) then
1246
1346
s%split_norm = global_splnorm_H
1247
1347
else
1265
1365
allocate(s%lambda(1:s%nzeta))
1266
1366
s%rc(1:s%nzeta) = 0.0d0
1267
1367
s%lambda(1:s%nzeta) = 1.0d0
1368
!
1369
! This option needs scale factors for multiple-z
1370
if (basp%basis_type.eq.'nonodes') then
1371
s%lambda(1) = 1.0d0
1372
if (s%nzeta > 1) then
1373
write(6,"(a)") "WARNING: NONODES: " //
1374
$ "Default scale factors set for shell (x0.8)"
1375
endif
1376
do i=2,s%nzeta
1377
s%lambda(i) = 0.8d0 * s%lambda(i-1)
1378
enddo
1379
endif
1268
1380
endif
1269
1381
enddo loop_l
1270
1382
1271
1383
if (basp%basis_size(3:3) .eq. 'p') then
1272
1384
1273
1385
! Polarization orbitals are defined so they have the minimum
1274
! angular momentum l such that there are not occupied orbitals
1386
! angular momentum l such that there are no occupied orbitals
1275
1387
! with the same l in the valence shell of the ground-state
1276
1388
! atomic configuration. They polarize the corresponding l-1
1277
1389
! shell.
1292
1404
! For example, in Ti, p is not occupied in the GS
1293
1405
! The following code will mark the top-most s shell
1294
1406
! as polarizable
1295
1296
if (.not. basp%ground_state%occupied(l)) then
1407
1408
! Note that l goes from 1 to 4
1409
if ( (.not. basp%ground_state%occupied(l)) ) then
1410
$
1411
if (basp%non_perturbative_polorbs) then
1412
ls=>basp%lshell(l)
1413
nsh = size(ls%shell)
1414
! This can happen, for example, for Ge (3d is semicore,
1415
! but the polarization shell is 4d (from 4p)
1416
!if (nsh /= 1) call die("Empty l-shell with nsh>1?")
1417
s => ls%shell(nsh)
1418
s%nzeta = nzeta_pol
1419
allocate(s%rc(1:s%nzeta))
1420
allocate(s%lambda(1:s%nzeta))
1421
s%rc(1:s%nzeta) = 0.0d0
1422
s%lambda(1:s%nzeta) = 1.0d0
1423
1424
s%polarization_shell = .true.
1425
1426
ls_parent=> basp%lshell(l-1)
1427
nsh = size(ls_parent%shell) ! Parent is the last (l-1) shell
1428
s%shell_being_polarized => ls_parent%shell(nsh)
1429
! Now, reset the polarization flags of the parent
1430
ls_parent%shell(nsh)%polarized = .false.
1431
ls_parent%shell(nsh)%was_polarized = .true.
1432
ls_parent%shell(nsh)%nzeta_pol = 0
1433
1434
EXIT loop_angmom ! Polarize only one shell!
1435
endif
1436
1437
! Normal treatment
1297
1438
ls=>basp%lshell(l-1)
1298
nsh = size(ls%shell) ! Use only the last shell
1439
nsh = size(ls%shell) ! Use only the last shell
1299
1440
s => ls%shell(nsh)
1300
1441
1301
1442
! Check whether shell to be polarized is occupied in the gs.
1302
1443
! (i.e., whether PAOs are going to be generated for it)
1303
1444
! If it is not, mark it for PAO generation anyway.
1315
1456
s%rc(1:s%nzeta) = 0.0d0
1316
1457
s%lambda(1:s%nzeta) = 1.0d0
1317
1458
endif
1459
! Put polarization flags
1318
1460
s%polarized = .true.
1319
1461
s%nzeta_pol = nzeta_pol
1320
1462
1321
exit loop_angmom ! Polarize only one shell!
1463
EXIT loop_angmom ! Polarize only one shell!
1322
1464
endif
1323
1465
enddo loop_angmom
1324
1466
Loggerhead is a web-based interface for Breezy
Version: 2.0.1