Siesta

Merge lp:~albertog/siesta/4.0-semic-pol-fix into lp:siesta/4.0

  1. 4.0-semic-pol-fix
  2. Merge into rel-4.0
Proposed by Alberto Garcia
Status: Merged
Merged at revision: 590
Proposed branch: lp:~albertog/siesta/4.0-semic-pol-fix
Merge into: lp:siesta/4.0
Diff against target: 1478 lines (+1092/-101)
11 files modified
Docs/siesta.tex (+13/-2)
Src/atom.F (+82/-50)
Src/basis_specs.f (+72/-16)
Src/basis_types.f (+27/-3)
Src/old_atmfuncs.f (+24/-28)
Tests/Makefile (+1/-1)
Tests/Reference/bessel-rich.out (+817/-0)
Tests/bessel-rich/bessel-rich.fdf (+45/-0)
Tests/bessel-rich/bessel-rich.pseudos (+1/-0)
Tests/bessel-rich/makefile (+6/-0)
version.info (+4/-1)
To merge this branch: bzr merge lp:~albertog/siesta/4.0-semic-pol-fix
Reviewer Review Type Date Requested Status
Nick Papior Pending
Javier Junquera Pending
Review via email: mp+362250@code.launchpad.net

Commit message

Fix some issues with polarization orbitals. Bessel orbs improvements.

A number of issues (mostly cosmetic, or appearing in rare occasions) with
polarization orbitals have been fixed.

* Fix 'n' quantum number for high-l polarization orbitals

  When a polarization orbital had angular momentum l greater than the
  rest of the the basis, its associated principal quantum number n was
  off (it was set to the polarized shell's n).

  This error affected the metadata in the .ion files, the wavefunction
  files, and the PDOS file, so some analysis operations might have
  failed to take into account those orbitals if they were requested
  explicitly (e.g., a COOP calculation with 'O_3d' symbols might show
  the bug; but not one with 'O_d' or 'O').

* Fix for the case of polarization orbitals having the same l as a
  semicore state.

  An example: For Ti, the electronic structure is 3s2 3p6 3d2 4s2 4p0*

  4p0* is a polarization orbital, and 3p6 can be a semicore orbital.
  The outer orbitals for l=0..3 reported by the 'periodic table'
  routines for Ti are:

      4s2 3p6 3d2 4f0

  When using a pseudopotential with 3p6 in the valence (a "semicore"
  state), if 4p did not appear explicitly in a PAO.Basis block (i.e.,
  with the entry for 4s marked as "polarized"), the program did not
  identify 3p as a "semicore state", and thus it did not generate
  *two* KB operators for l=1.

  This has been fixed, together with the wrong 'n' assignment to the
  polarization orbital.

* Improve output of basis specification

  Change misleading 'n=' prefix by 'i=' (it is simply an enumeration)
  and add a proper symbolic nl identifier at the end of the line.

Some cosmetic and usability improvements have also been made to the
handling of Bessel floating orbitals:

* Extend and clarify the handling of Bessel floating orbitals

  For Bessel floating orbitals different 'zetas' correspond to
  progressively excited states (more nodes) for a given l.

  A user might want to reproduce the basis coverage of a given 'nl'
  shell with a suite of bessel functions with increasing number of
  nodes, using the 'zeta' mechanism. Up to now, it was only possible to
  do this for a single 'n' for each l, and the 'n' quantum number
  specified by the user in the (compulsory for Bessel orbitals)
  PAO.Basis block was ignored.

  This patch enables the specification of multiple 'nl' shells for the
  same l, and keeps the user's choice for the 'n' quantum number in
  the output files of the program (.ion, .PDOS, etc).

* Enable short-hand for entering rc information in PAO.Basis block

  If the number of rc's for a given shell is less than the number of
  'zetas', the program will assign the last rc value to the remaining
  zetas, rather than stopping with an error. This is particularly
  useful for Bessel suites of orbitals.

  A new test (Tests/bessel-rich) exemplifies this, as well as the
  new Bessel extensions.

Description of the change

I had a new inspiration and I think I have fixed all the issues in a simpler and less intrusive way.
'n' quantum numbers are now assigned correctly, and the KB issue for Ti is solved.

The underlying problem is documented in an "archaeological note" at the end of the header of basis_specs.f.

The "basis specs" are now more clearly output, with symbolic names for shells (similar to the
output of the psml branch).

Thank you for your patience!

To post a comment you must log in.
Revision history for this message
Nick Papior (nickpapior) wrote :

Looks good to me!

lp:~albertog/siesta/4.0-semic-pol-fix updated
599. By Alberto Garcia

Update manual. Extend 'assumed rc' feature to contraction factors

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex 2018-07-19 19:31:44 +0000
+++ Docs/siesta.tex 2019-01-30 13:27:17 +0000
@@ -1669,7 +1669,11 @@
1669\item[-] {\tt l}: Angular momentum of1669\item[-] {\tt l}: Angular momentum of
1670basis orbitals of this shell1670basis orbitals of this shell
1671\item[-] {\tt nzls(lsh,is)}: Number of 'zetas' for this shell. For a filteret 1671\item[-] {\tt nzls(lsh,is)}: Number of 'zetas' for this shell. For a filteret
1672basis this number is ignored since the number is controlled by the cutoff.1672basis this number is ignored since the number is controlled by the
1673cutoff.
1674For bessel-floating orbitals, the different 'zetas' map to
1675increasingly excited states with the same angular momentum (with
1676increasing number of nodes).
1673\item[-] {\tt PolOrb(l+1)}: {\it Optional input}. If set equal to {\tt P}, a1677\item[-] {\tt PolOrb(l+1)}: {\it Optional input}. If set equal to {\tt P}, a
1674shell of1678shell of
1675polarization functions (with angular momentum $l+1$) will be constructed1679polarization functions (with angular momentum $l+1$) will be constructed
@@ -1718,9 +1722,16 @@
1718each 'zeta' for this shell. For the second zeta onwards, if this value1722each 'zeta' for this shell. For the second zeta onwards, if this value
1719is negative, the actual rc used will be the given fraction of the1723is negative, the actual rc used will be the given fraction of the
1720first zeta's rc.1724first zeta's rc.
1725If the number of rc's for a given shell is less than the number of
1726'zetas', the program will assign the last rc value to the remaining
1727zetas, rather than stopping with an error. This is particularly
1728useful for Bessel suites of orbitals.
1721\item[-] {\tt contrf(izeta,l,is)}: {\it Optional input}.1729\item[-] {\tt contrf(izeta,l,is)}: {\it Optional input}.
1722Contraction factor\index{scale factor} of1730Contraction factor\index{scale factor} of
1723each 'zeta' for this shell.1731each 'zeta' for this shell.
1732If the number of entries for a given shell is less than the number of
1733'zetas', the program will assign the last contraction value to the remaining
1734zetas, rather than stopping with an error.
1724{\it Default value}: {\tt 1.0}1735{\it Default value}: {\tt 1.0}
1725\end{itemize}1736\end{itemize}
17261737
@@ -1751,7 +1762,7 @@
1751a constant-pseudopotential atom, {\it i.e.}, the basis functions1762a constant-pseudopotential atom, {\it i.e.}, the basis functions
1752generated will be spherical Bessel functions\index{Bessel functions}1763generated will be spherical Bessel functions\index{Bessel functions}
1753with the specified $r_c$. In this case, $r_c$ has to be given, as1764with the specified $r_c$. In this case, $r_c$ has to be given, as
1754{\bf EnergyShift} will not calculate it.\index{basis!Bessel functions}1765{\bf PAO.EnergyShift} will not calculate it.\index{basis!Bessel functions}
17551766
1756Other negative atomic numbers will be interpreted by {\sc Siesta} as1767Other negative atomic numbers will be interpreted by {\sc Siesta} as
1757{\it ghosts}\index{ghost atoms}\index{basis!ghost atoms}1768{\it ghosts}\index{ghost atoms}\index{basis!ghost atoms}
17581769
=== modified file 'Src/atom.F'
--- Src/atom.F 2018-04-26 13:04:08 +0000
+++ Src/atom.F 2019-01-30 13:27:17 +0000
@@ -289,8 +289,9 @@
289 . ' ( Floating basis ) ' 289 . ' ( Floating basis ) '
290290
291 elseif(iz.eq.-100) then291 elseif(iz.eq.-100) then
292 write(6,'(a,i4,a)')292 write(6,'(3a,i4,a,a)')
293 . 'atom: Called for Z=',iz,'( Floating Bessel functions)' 293 . 'atom: Called for ', atm_label, ' (Z =', iz,')',
294 . ' ( Floating Bessel functions ) '
294295
295 endif296 endif
296! 297!
@@ -702,7 +703,7 @@
702 slfe(is) = 0.0d0703 slfe(is) = 0.0d0
703! Calculate Bessel functions704! Calculate Bessel functions
704 call Bessel(is,a,b,rofi,drdi,s,705 call Bessel(is,a,b,rofi,drdi,s,
705 . lmxo,nzeta,rco,lambda,filtercut,no)706 . lmxo,nzeta,rco,lambda,filtercut,no,nsemic)
706 707
707 write(6,'(/a,i3)')708 write(6,'(/a,i3)')
708 . 'atom: Total number of floating Bessel orbitals:', no709 . 'atom: Total number of floating Bessel orbitals:', no
@@ -2578,7 +2579,8 @@
2578C2579C
2579C Internal and common variables2580C Internal and common variables
2580C2581C
2581 integer l, lmax, nzetamax, nkblmx,nsm, nsm_max2582 integer l, lmax, nzetamax, nkblmx,nsm, nsm_max, n_user
2583 integer npols_lm1
2582C2584C
2583C Adding a new species to the list2585C Adding a new species to the list
2584C 2586C
@@ -2635,27 +2637,42 @@
2635 2637
2636 cnfigtb(:,:,is) = 02638 cnfigtb(:,:,is) = 0
26372639
2638 if (iz.ne.-100) then 2640 ! Generate the array with 'n' quantum number information
2639 do l=0,lmxo2641 if (iz.ne.-100) then ! Normal case
2640 do nsm=1, nsemic(l)+12642
2641 cnfigtb(l,nsm,is)=cnfigmx(l)-(nsemic(l)+1)+nsm2643 ! Note that cnfigmx and nsemic are now properly extended
2642 enddo2644 ! all the way to lmaxd (in basis_types) to give a proper cnfigtb
2643 enddo 2645
2644 do l=min(lmaxd,lmxo+1),lmaxd
2645 cnfigtb(l,1,is)=l+1 ! AG*** Why??
2646 ! To deal with polarization
2647 ! orbitals beyond lmxo?
2648! BUT, what happens to "in-between" polarization orbitals?
2649!
2650!
2651 enddo
2652 else
2653!
2654! Arbitrary "n" for "semicore" bessel states...
2655!
2656 do l=0,lmaxd2646 do l=0,lmaxd
2647 do nsm=1, nsemic(l)+1
2648 cnfigtb(l,nsm,is)=cnfigmx(l)-(nsemic(l)+1)+nsm
2649 enddo
2650
2651 if (l==0) cycle
2652
2653 ! Special case where a polarization orbital has the same l
2654 ! as a lower-lying state (e.g. Ti: 3s2 3p6 4s2 3d2 4p0*) 4p0 is a polarization
2655 ! orbital, with the same l as 3p.
2656 ! (See 'archaeological note' in the header of 'basis_specs.f')
2657 npols_lm1 = sum(npolorb(l-1,:))
2658 if (npols_lm1 > 0) then
2659 nsm = nsemic(l)+2
2660 cnfigtb(l,nsm,is)= cnfigmx(l) + 1 ! in the above case, 4p
2661 endif
2662
2663 enddo
2664
2665 else
2666!
2667! For bessel-function states, the 'n' quantum number
2668! does not follow the heuristics for real 'atomic'
2669! states, but since the user has to specify them in the
2670! PAO.Basis block, we use them.
2671
2672 do l=0,lmxo
2657 do nsm=1,nsemic(l)+12673 do nsm=1,nsemic(l)+1
2658 cnfigtb(l,nsm,is)=nsm2674 n_user = cnfigmx(l)-(nsemic(l)+1)+nsm
2675 cnfigtb(l,nsm,is) = n_user
2659 enddo2676 enddo
2660 enddo 2677 enddo
2661 endif ! iz.ne.-1002678 endif ! iz.ne.-100
@@ -6727,7 +6744,7 @@
6727!6744!
6728 subroutine BESSEL(is,a,b,rofi,drdi,s,lmxo,6745 subroutine BESSEL(is,a,b,rofi,drdi,s,lmxo,
6729 . nzeta,rco,lambda,filtercut,6746 . nzeta,rco,lambda,filtercut,
6730 . norb) 6747 . norb,nsemic)
67316748
6732 use basis_specs, only: restricted_grid6749 use basis_specs, only: restricted_grid
6733C6750C
@@ -6744,14 +6761,13 @@
6744 . lambda(nzetmx, 0:lmaxd,nsemx),6761 . lambda(nzetmx, 0:lmaxd,nsemx),
6745 . filtercut(0:lmaxd,nsemx)6762 . filtercut(0:lmaxd,nsemx)
67466763
6747
6748 integer6764 integer
6749 . lmxo, is, nzeta(0:lmaxd,nsemx), norb6765 . lmxo, is, nzeta(0:lmaxd,nsemx), norb, nsemic(0:lmaxd)
6750C6766C
6751C Internal variables6767C Internal variables
6752C6768C
6753 integer6769 integer
6754 . l,nprin, nnodes, nodd, nrc, ir, indx, izeta6770 . l,nprin, nnodes, nodd, nrc, ir, indx, izeta, nsm
67556771
6756 real(dp)6772 real(dp)
6757 . rc, dnrm, phi,6773 . rc, dnrm, phi,
@@ -6763,17 +6779,28 @@
6763 indx=06779 indx=0
6764 do l=0,lmxo 6780 do l=0,lmxo
67656781
6766 if (nzeta(l,1).gt.0) then6782 do nsm=1,nsemic(l)+1
67676783 if (nzeta(l,nsm).gt.0) then
6768 write(6,'(/2A,I2)')6784 write(6,'(/2A,I2)')
6769 . 'Bessel: floating Bessel functions ',6785 . 'Bessel: floating Bessel functions ',
6770 . 'with angular momentum L=',l6786 . 'with angular momentum L=',l
67716787 exit
6772 do izeta=1, nzeta(l,1) 6788 endif
6789 enddo
6790
6791 do nsm=1,nsemic(l)+1 ! Note extra shell loop
6792
6793 if (nzeta(l,nsm).eq.0) CYCLE
6794
6795 write(6,'(/A,I1,A)')
6796 . 'Bessel: Basis orbitals for "state" ',
6797 . cnfigtb(l,nsm,is), sym(l)
6798
6799 do izeta=1, nzeta(l,nsm)
6773C6800C
6774C Cut-off radius for Bessel functions must be an explicit input6801C Cut-off radius for Bessel functions must be an explicit input
6775C6802C
6776 if (rco(izeta,l,1).lt.1.0d-5) then 6803 if (rco(izeta,l,nsm).lt.1.0d-5) then
6777 write(6,'(a)')6804 write(6,'(a)')
6778 . 'Bessel: ERROR Zero cut-off radius with Z=-100 option'6805 . 'Bessel: ERROR Zero cut-off radius with Z=-100 option'
6779 write(6,'(a)')6806 write(6,'(a)')
@@ -6783,15 +6810,15 @@
6783 call die6810 call die
6784 endif6811 endif
6785C6812C
6786 if (abs(lambda(izeta,l,1)).lt.1.0d-3) 6813 if (abs(lambda(izeta,l,nsm)).lt.1.0d-3)
6787 . lambda(izeta,l,1)=1.0d06814 . lambda(izeta,l,nsm)=1.0d0
6788 if (abs(lambda(izeta,l,1)-1.0d0).gt.1.0d-3) then6815 if (abs(lambda(izeta,l,nsm)-1.0d0).gt.1.0d-3) then
6789 write(6,'(/,a)')6816 write(6,'(/,a)')
6790 . 'Bessel: WARNING Scale factor is not active with Z=-100 option' 6817 . 'Bessel: WARNING Scale factor is not active with Z=-100 option'
6791 endif 6818 endif
6792 lambda(izeta,l,1)=1.0d06819 lambda(izeta,l,nsm)=1.0d0
67936820
6794 rc=rco(izeta,l,1)6821 rc=rco(izeta,l,nsm)
6795 nrc=nint(log(rc/b+1.0d0)/a)+16822 nrc=nint(log(rc/b+1.0d0)/a)+1
6796 if (restricted_grid) then6823 if (restricted_grid) then
6797 nodd=mod(nrc,2)6824 nodd=mod(nrc,2)
@@ -6800,10 +6827,17 @@
6800 endif6827 endif
6801 endif6828 endif
6802 rc=b*(exp(a*(nrc-1))-1.0d0)6829 rc=b*(exp(a*(nrc-1))-1.0d0)
6803 rco(izeta,l,1)=rc6830 rco(izeta,l,nsm)=rc
68046831
6805 nnodes=izeta6832 ! Instead of generating directly j_l(k*r), where
6806 nprin=l+16833 ! k is chosen to have a zero at rc, the effective Schrodinger
6834 ! equation is integrated with a zero potential.
6835
6836 ! Note that the the different zetas here mean successively
6837 ! excited states for this l
6838
6839 nnodes=izeta ! Including the one at r=rc...
6840 nprin=l+1
6807 call schro_eq(1.0d0,rofi,v,v,s,drdi,6841 call schro_eq(1.0d0,rofi,v,v,s,drdi,
6808 . nrc,l,a,b,nnodes,nprin,6842 . nrc,l,a,b,nnodes,nprin,
6809 . eorb,g) 6843 . eorb,g)
@@ -6826,19 +6860,17 @@
68266860
6827 write(6,'(/,(3x,a,i2),2(/,a25,f12.6))')6861 write(6,'(/,(3x,a,i2),2(/,a25,f12.6))')
6828 . 'izeta =',izeta,6862 . 'izeta =',izeta,
6829 . 'rc =',rco(izeta,l,1),6863 . 'rc =',rco(izeta,l,nsm),
6830 . 'energy =',eorb 6864 . 'energy =',eorb
68316865
6832 norb=norb+(2*l+1)6866 norb=norb+(2*l+1)
6833 indx=indx+16867 indx=indx+1
6834 call comBasis(is,a,b,rofi,g,l,6868 call comBasis(is,a,b,rofi,g,l,
6835 . rco(izeta,l,1),lambda(izeta,l,1),6869 . rco(izeta,l,nsm),lambda(izeta,l,nsm),
6836 . filtercut(l,1),izeta,1,nrc,indx)6870 . filtercut(l,nsm),izeta,nsm,nrc,indx)
68376871
6838 enddo 6872 enddo
6839 6873 enddo
6840 endif
6841
6842 enddo 6874 enddo
68436875
6844 end subroutine bessel6876 end subroutine bessel
68456877
=== modified file 'Src/basis_specs.f'
--- Src/basis_specs.f 2017-12-15 10:24:23 +0000
+++ Src/basis_specs.f 2019-01-30 13:27:17 +0000
@@ -47,9 +47,11 @@
47! 47!
48! rc_1 rc_2 .... rc_nzeta 48! rc_1 rc_2 .... rc_nzeta
49! 49!
50! are the cutoff radii in bohrs. This line is mandatory and there must50! are the cutoff radii in bohrs. This line is mandatory
51! be at least nzeta values (the extra ones are discarded)51! If the number of rc's for a given shell is less than the number of
52! 52! 'zetas', the program will assign the last rc value to the remaining
53! zetas, rather than stopping with an error. This is particularly
54! useful for Bessel suites of orbitals.
53! A line containing contraction (scale) factors is optional, but if it55! A line containing contraction (scale) factors is optional, but if it
54! appears, the values *must be* real numbers, and there must be at56! appears, the values *must be* real numbers, and there must be at
55! least nzeta of them.57! least nzeta of them.
@@ -76,7 +78,7 @@
76! There are no 'per l-shell' polarization orbitals, except if the78! There are no 'per l-shell' polarization orbitals, except if the
77! third character of 'basis_size' is 'p' (as in 'dzp'), in which79! third character of 'basis_size' is 'p' (as in 'dzp'), in which
78! case polarization orbitals are defined so they have the minimum 80! case polarization orbitals are defined so they have the minimum
79! angular momentum l such that there are not occupied orbitals 81! angular momentum l such that there are no occupied orbitals
80! with the same l in the valence shell of the ground-state 82! with the same l in the valence shell of the ground-state
81! atomic configuration. They polarize the corresponding l-1 shell.83! atomic configuration. They polarize the corresponding l-1 shell.
82! 84!
@@ -86,7 +88,7 @@
86! 88!
87! rc(1:nzeta) is set to 0.089! rc(1:nzeta) is set to 0.0
88! lambda(1:nzeta) is set to 1.0 (this is a change from old practice)90! lambda(1:nzeta) is set to 1.0 (this is a change from old practice)
89! 91!
90! ----------------------------------92! ----------------------------------
91! 93!
92! Next come the blocks associated to the KB projectors:94! Next come the blocks associated to the KB projectors:
@@ -115,7 +117,28 @@
115! n-shells in the corresponding PAO shell with the same l. For l117! n-shells in the corresponding PAO shell with the same l. For l
116! greater than lmxo, it is set to 1. The reference energies are118! greater than lmxo, it is set to 1. The reference energies are
117! in all cases set to huge(1.d0) to mark the default.119! in all cases set to huge(1.d0) to mark the default.
118! 120!
121! Archaeological note: The first implementation of the basis-set generation
122! module had only non-polarization orbitals. Polarization orbitals were added
123! later as "second-class" companions. This shows in details like "lmxo" (the
124! maximum l of the basis set) not taking into account polarization orbitals.
125! Polarization orbitals were tagged at the end, without maintaining l-shell
126! ordering.
127! Even later, support for "semicore" orbitals was added. A new ('nsm') index was
128! used to distinguish the different orbitals in a l-shell. Polarization orbitals
129! were not brought into this classification.
130! The code is strained when semicore and polarization orbitals coexist for the same l,
131! as in Ti, whose electronic structure is []3s2 3p6 3d2 4s2 (4p0*) with 4p as the
132! polarization orbital.
133! Here 'nsemic' (the number of semicore states) is kept at zero for l=1
134! (the polarization orbital is not counted), with the side-effect that only one KB
135! projector is generated for l=1.
136! Special checks have been implemented to cover these cases.
137!
138! Future work should probably remove the separate treatment of polarization orbitals.
139! (Note that if *all* orbitals are specified in a PAO.Basis block, in effect turning
140! perturbative polarization orbitals into 'normal' orbitals, this problem is not present.)
141!
119! =======================================================================142! =======================================================================
120!143!
121 use precision144 use precision
@@ -489,8 +512,21 @@
489 k%l = l512 k%l = l
490 if (l.gt.basp%lmxo) then513 if (l.gt.basp%lmxo) then
491 k%nkbl = 1514 k%nkbl = 1
492 else ! Set equal to the number of PAO shells515 else
493 k%nkbl = basp%lshell(l)%nn516 ! Set equal to the number of PAO shells with this l
517 k%nkbl = basp%lshell(l)%nn
518 ! Should include polarization orbs (as in Ti case: 3p..4p*)
519 ! (See 'archaeological note' in the header of this file)
520 if (l>0) then
521 do i = 1, basp%lshell(l-1)%nn
522 if (basp%lshell(l-1)%shell(i)%polarized) then
523 k%nkbl = k%nkbl + 1
524 write(6,"(a,i1,a)") trim(basp%label) //
525 $ ': nkbl increased for l=',l,
526 $ ' due to the presence of a polarization orbital'
527 endif
528 enddo
529 endif
494 if (k%nkbl.eq.0) then530 if (k%nkbl.eq.0) then
495 write(6,*) 'Warning: Empty PAO shell. l =', l531 write(6,*) 'Warning: Empty PAO shell. l =', l
496 write(6,*) 'Will have a KB projector anyway...'532 write(6,*) 'Will have a KB projector anyway...'
@@ -592,6 +628,7 @@
592 subroutine repaobasis()628 subroutine repaobasis()
593629
594 integer isp, ish, nn, i, ind, l, indexp, index_splnorm630 integer isp, ish, nn, i, ind, l, indexp, index_splnorm
631 integer nrcs_zetas
595632
596 type(block_fdf) :: bfdf633 type(block_fdf) :: bfdf
597 type(parsed_line), pointer :: pline634 type(parsed_line), pointer :: pline
@@ -742,11 +779,21 @@
742 s%rc(:) = 0.d0779 s%rc(:) = 0.d0
743 s%lambda(:) = 1.d0780 s%lambda(:) = 1.d0
744 if (.not. fdf_bline(bfdf,pline)) call die("No rc's")781 if (.not. fdf_bline(bfdf,pline)) call die("No rc's")
745 if (fdf_bnvalues(pline) .ne. s%nzeta)782
746 . call die("Wrong number of rc's")783 ! Use the last rc entered for the successive zetas
784 ! if there are not enough values (useful for Bessel)
785 nrcs_zetas = fdf_bnvalues(pline)
786 if (nrcs_zetas < 1) then
787 call die("Need at least one rc per shell in PAO.Basis block")
788 endif
747 do i= 1, s%nzeta789 do i= 1, s%nzeta
748 s%rc(i) = fdf_bvalues(pline,i)790 if (i <= nrcs_zetas) then
791 s%rc(i) = fdf_bvalues(pline,i)
792 else
793 s%rc(i) = s%rc(nrcs_zetas)
794 endif
749 enddo795 enddo
796
750 if (s%split_norm_specified) then797 if (s%split_norm_specified) then
751 do i = 2,s%nzeta798 do i = 2,s%nzeta
752 if (s%rc(i) /= 0.0_dp) then799 if (s%rc(i) /= 0.0_dp) then
@@ -770,11 +817,20 @@
770 . call die('repaobasis: ERROR in PAO.Basis block')817 . call die('repaobasis: ERROR in PAO.Basis block')
771 cycle shells818 cycle shells
772 else819 else
773 if (fdf_bnreals(pline) .ne. s%nzeta)820 ! Read scale factors
774 . call die("Wrong number of lambda's")821 ! Use the last scale factor entered for the successive zetas
775 do i=1,s%nzeta822 ! if there are not enough values
776 s%lambda(i) = fdf_breals(pline,i)823 nrcs_zetas = fdf_bnreals(pline)
777 enddo824 if (nrcs_zetas < 1) then
825 call die("Need at least one scale factor in PAO.Basis")
826 endif
827 do i= 1, s%nzeta
828 if (i <= nrcs_zetas) then
829 s%lambda(i) = fdf_breals(pline,i)
830 else
831 s%lambda(i) = s%lambda(nrcs_zetas)
832 endif
833 enddo
778 endif834 endif
779 endif835 endif
780836
781837
=== modified file 'Src/basis_types.f'
--- Src/basis_types.f 2017-12-15 10:24:23 +0000
+++ Src/basis_types.f 2019-01-30 13:27:17 +0000
@@ -578,6 +578,24 @@
578 $ cnfigmx(l,isp) = basp%ground_state%n(l)578 $ cnfigmx(l,isp) = basp%ground_state%n(l)
579579
580 enddo580 enddo
581
582 ! NOTE: cnfigmx and nsemic are only initialized for l up to lmxo
583 ! in the above loop
584 ! Extend them so that we can deal properly with outer polarization states
585
586 do l=basp%lmxo+1, lmaxd
587 ! gs is only setup up to l=3 (f)
588 if (l <= 3) then
589 cnfigmx(l,isp) = basp%ground_state%n(l)
590 else
591 ! g orbitals. Use the "l+1" heuristic
592 ! For example, 5g pol orb associated to a 4f orb.
593 cnfigmx(l,isp) = l + 1
594 endif
595 nsemic(l,isp) = 0
596 enddo
597
598
581 do l=0,basp%lmxkb599 do l=0,basp%lmxkb
582 k=>basp%kbshell(l)600 k=>basp%kbshell(l)
583 nkbl(l,isp) = k%nkbl601 nkbl(l,isp) = k%nkbl
@@ -594,6 +612,10 @@
594 integer, intent(in) :: is612 integer, intent(in) :: is
595613
596 integer l, n, i614 integer l, n, i
615 integer :: nprin
616 character(len=4) :: orb_id
617 character(len=1), parameter ::
618 $ sym(0:4) = (/ 's','p','d','f','g' /)
597619
598 write(lun,'(/a/79("="))') '<basis_specs>'620 write(lun,'(/a/79("="))') '<basis_specs>'
599 write(lun,'(a20,1x,a2,i4,4x,a5,g12.5,4x,a7,g12.5)')621 write(lun,'(a20,1x,a2,i4,4x,a5,g12.5,4x,a7,g12.5)')
@@ -609,9 +631,11 @@
609 $ 'Cnfigmx=', cnfigmx(l,is)631 $ 'Cnfigmx=', cnfigmx(l,is)
610 do n=1,nsemic(l,is)+1632 do n=1,nsemic(l,is)+1
611 if (nzeta(l,n,is) == 0) exit633 if (nzeta(l,n,is) == 0) exit
612 write(lun,'(10x,a2,i1,2x,a6,i1,2x,a7,i1)')634 nprin = cnfigmx(l,is) - nsemic(l,is) + n - 1
613 $ 'n=', n, 'nzeta=',nzeta(l,n,is),635 write(orb_id,"(a1,i1,a1,a1)") "(",nprin, sym(l), ")"
614 $ 'polorb=', polorb(l,n,is)636 write(lun,'(10x,a2,i1,2x,a6,i1,2x,a7,i1,2x,a4)')
637 $ 'i=', n, 'nzeta=',nzeta(l,n,is),
638 $ 'polorb=', polorb(l,n,is), orb_id
615 if (basistype(is).eq.'filteret') then639 if (basistype(is).eq.'filteret') then
616 write(lun,'(10x,a10,2x,g12.5)') 640 write(lun,'(10x,a10,2x,g12.5)')
617 $ 'fcutoff:', filtercut(l,n,is)641 $ 'fcutoff:', filtercut(l,n,is)
618642
=== modified file 'Src/old_atmfuncs.f'
--- Src/old_atmfuncs.f 2016-01-25 16:00:16 +0000
+++ Src/old_atmfuncs.f 2019-01-30 13:27:17 +0000
@@ -399,18 +399,17 @@
399 integer, intent(in) :: is ! Species index399 integer, intent(in) :: is ! Species index
400 integer, intent(in) :: io ! Orbital index (within atom)400 integer, intent(in) :: io ! Orbital index (within atom)
401401
402C Returns the valence-shell configuration in the atomic ground state402C INTEGER CNFIGFIO: Principal quantum number of the shell to which
403C (i.e. the principal quatum number for orbitals of angular momentum l)403C the orbital belongs
404404
405C INTEGER CNFIGFIO: Principal quantum number of the shell to what 405C (Formerly, for polarization orbitals
406C the orbital belongs ( for polarization orbitals406C the quantum number corresponded to the shell which
407C the quantum number corresponds to the shell which407C is polarized by the orbital io.
408C is polarized by the orbital io) 408C This behavior for polarization orbitals is meaningless,
409409C and has been replaced.)
410 integer l, norb, lorb, izeta, ipol,nsm410
411 integer indx, nsmorb411 integer l, norb, izeta, ipol, nsm
412412
413C
414 call check_is('cnfigfio',is)413 call check_is('cnfigfio',is)
415 if ((io.gt.nomax(is)).or.(io.lt.1)) then414 if ((io.gt.nomax(is)).or.(io.lt.1)) then
416 write(6,*) 'CNFIGFIO: THERE ARE NO DATA FOR IO=',IO415 write(6,*) 'CNFIGFIO: THERE ARE NO DATA FOR IO=',IO
@@ -419,25 +418,32 @@
419 call die()418 call die()
420 endif419 endif
421420
421 ! This routine assumes that polarization orbitals are the last ones
422 ! in the list indexed by 'io'
423
422 norb=0424 norb=0
423 indx=0
424 do 10 l=0,lmxosave(is)425 do 10 l=0,lmxosave(is)
425 do 8 nsm=1,nsemicsave(l,is)+1426 do 8 nsm=1,nsemicsave(l,is)+1
426 do 5 izeta=1,nzetasave(l,nsm,is)427 do 5 izeta=1,nzetasave(l,nsm,is)
427 norb=norb+(2*l+1)428 norb=norb+(2*l+1)
428 indx=indx+1429 if(norb.ge.io) then
429 if(norb.ge.io) goto 30430 cnfigfio=cnfigtb(l,nsm,is)
431 return
432 endif
430 5 continue433 5 continue
431 8 continue434 8 continue
43210 continue43510 continue
433436
434 indx=0
435 do 20 l=0,lmxosave(is)437 do 20 l=0,lmxosave(is)
436 do 18 nsm=1,nsemicsave(l,is)+1438 do 18 nsm=1,nsemicsave(l,is)+1
437 do 15 ipol=1, npolorbsave(l,nsm,is)439 do 15 ipol=1, npolorbsave(l,nsm,is)
438 norb=norb+(2*(l+1)+1)440 norb=norb+(2*(l+1)+1)
439 indx=indx+1441 if(norb.ge.io) then
440 if(norb.ge.io) goto 40442 ! Deal properly with polarization orbitals
443 ! Return the highest n for l+1
444 cnfigfio=maxval(cnfigtb(l+1,:,is))
445 return
446 endif
44115 continue44715 continue
44218 continue44818 continue
44320 continue44920 continue
@@ -445,16 +451,6 @@
445 write(6,*) 'CNFIGFIO: ERROR: NOT FOUND'451 write(6,*) 'CNFIGFIO: ERROR: NOT FOUND'
446 call die()452 call die()
447453
44830 lorb=l
449 nsmorb=nsm
450 cnfigfio=cnfigtb(lorb,nsmorb,is)
451 return
452
45340 lorb=l
454 nsmorb=nsm
455 cnfigfio=cnfigtb(lorb,nsmorb,is)
456 return
457
458 end function cnfigfio454 end function cnfigfio
459!455!
460!456!
461457
=== modified file 'Tests/Makefile'
--- Tests/Makefile 2018-04-25 11:14:24 +0000
+++ Tests/Makefile 2019-01-30 13:27:17 +0000
@@ -36,7 +36,7 @@
36#36#
37label = work37label = work
3838
39tests = 32_h2o ag anneal-cont ar2_vdw batio3 benzene bessel born \39tests = 32_h2o ag anneal-cont ar2_vdw batio3 benzene bessel bessel-rich born \
40 born_spin carbon_nanoscroll ch4 chargeconf-h2o constant_volume \40 born_spin carbon_nanoscroll ch4 chargeconf-h2o constant_volume \
41 dipole_correction fe fe_broyden fe_clust_noncollinear fe_cohp \41 dipole_correction fe fe_broyden fe_clust_noncollinear fe_cohp \
42 fen fe_noncol_kp fire_benzene floating force_2 force_constants \42 fen fe_noncol_kp fire_benzene floating force_2 force_constants \
4343
=== added file 'Tests/Reference/bessel-rich.out'
--- Tests/Reference/bessel-rich.out 1970-01-01 00:00:00 +0000
+++ Tests/Reference/bessel-rich.out 2019-01-30 13:27:17 +0000
@@ -0,0 +1,817 @@
1Siesta Version : siesta-4.0--589--n-fix-4
2Architecture : gfortran-macosx64-openmpi
3Compiler version: GNU Fortran (Homebrew GCC 7.2.0) 7.2.0
4Compiler flags : mpif90 -g -O2
5PP flags : -DCDF -DMPI -DF2003
6PARALLEL version
7NetCDF support
8
9* Running in serial mode with MPI
10>> Start of run: 17-JAN-2019 15:47:45
11
12 ***********************
13 * WELCOME TO SIESTA *
14 ***********************
15
16reinit: Reading from standard input
17************************** Dump of input data file ****************************
18#
19# Needs new feature: handling of fewer rc's than nzetas in PAO.Basis block
20#
21write-ion-plot-files T
22#
23SystemName Water molecule with various Bessel Orbitals
24SystemLabel bessel-rich
25NumberofAtoms 7
26NumberOfSpecies 4
27%block ChemicalSpeciesLabel
28 1 8 O # Species index, atomic number, species label
29 2 1 H
30 3 -100 Bessel
31 4 -100 J
32%endblock ChemicalSpeciesLabel
33AtomicCoordinatesFormat Ang
34%block AtomicCoordinatesAndAtomicSpecies
35 0.000 0.000 0.000 1
36 0.757 0.586 0.000 2
37-0.757 0.586 0.000 2
38 0.3785 0.293 0.000 3
39-0.3785 0.293 0.000 3
40 0.3785 0.293 0.000 4
41-0.3785 0.293 0.000 4
42%endblock AtomicCoordinatesAndAtomicSpecies
43%block PAO.Basis
44Bessel 3
45 n=1 0 1
46 2.0
47 1.0
48 n=2 0 1
49 2.5
50 1.0
51 n=3 1 1
52 3.5
53 1.0
54J 2 # l-shells
55n=2 0 7 # Note new feature: fewer rc's than zetas
56 4.5
57n=2 1 3
58 4.5 4.5 5.0
59%endblock PAO.Basis
60************************** End of input data file *****************************
61
62reinit: -----------------------------------------------------------------------
63reinit: System Name: Water molecule with various Bessel Orbitals
64reinit: -----------------------------------------------------------------------
65reinit: System Label: bessel-rich
66reinit: -----------------------------------------------------------------------
67
68initatom: Reading input for the pseudopotentials and atomic orbitals ----------
69 Species number: 1 Label: O Atomic number: 8
70 Species number: 2 Label: H Atomic number: 1
71 Species number: 3 Label: Bessel (floating Bessel functions)
72 Species number: 4 Label: J (floating Bessel functions)
73Ground state valence configuration: 2s02 2p04
74Reading pseudopotential information in formatted form from O.psf
75
76Valence configuration for pseudopotential generation:
772s( 2.00) rc: 1.14
782p( 4.00) rc: 1.14
793d( 0.00) rc: 1.14
804f( 0.00) rc: 1.14
81Dumping pseudopotential information in formatted form in O.psdump
82Ground state valence configuration: 1s01
83Reading pseudopotential information in formatted form from H.psf
84
85Valence configuration for pseudopotential generation:
861s( 1.00) rc: 1.25
872p( 0.00) rc: 1.25
883d( 0.00) rc: 1.25
894f( 0.00) rc: 1.25
90Dumping pseudopotential information in formatted form in H.psdump
91For O, standard SIESTA heuristics set lmxkb to 3
92 (one more than the basis l, including polarization orbitals).
93Use PS.lmax or PS.KBprojectors blocks to override.
94For H, standard SIESTA heuristics set lmxkb to 2
95 (one more than the basis l, including polarization orbitals).
96Use PS.lmax or PS.KBprojectors blocks to override.
97
98<basis_specs>
99===============================================================================
100O Z= 8 Mass= 16.000 Charge= 0.17977+309
101Lmxo=1 Lmxkb= 3 BasisType=split Semic=F
102L=0 Nsemic=0 Cnfigmx=2
103 n=1 nzeta=2 polorb=0
104 splnorm: 0.15000
105 vcte: 0.0000
106 rinn: 0.0000
107 qcoe: 0.0000
108 qyuk: 0.0000
109 qwid: 0.10000E-01
110 rcs: 0.0000 0.0000
111 lambdas: 1.0000 1.0000
112L=1 Nsemic=0 Cnfigmx=2
113 n=1 nzeta=2 polorb=1
114 splnorm: 0.15000
115 vcte: 0.0000
116 rinn: 0.0000
117 qcoe: 0.0000
118 qyuk: 0.0000
119 qwid: 0.10000E-01
120 rcs: 0.0000 0.0000
121 lambdas: 1.0000 1.0000
122-------------------------------------------------------------------------------
123L=0 Nkbl=1 erefs: 0.17977+309
124L=1 Nkbl=1 erefs: 0.17977+309
125L=2 Nkbl=1 erefs: 0.17977+309
126L=3 Nkbl=1 erefs: 0.17977+309
127===============================================================================
128</basis_specs>
129
130atom: Called for O (Z = 8)
131
132read_vps: Pseudopotential generation method:
133read_vps: ATM3 Troullier-Martins
134Total valence charge: 6.00000
135
136xc_check: Exchange-correlation functional:
137xc_check: Ceperley-Alder
138V l=0 = -2*Zval/r beyond r= 1.1278
139V l=1 = -2*Zval/r beyond r= 1.1278
140V l=2 = -2*Zval/r beyond r= 1.1278
141V l=3 = -2*Zval/r beyond r= 1.1138
142All V_l potentials equal beyond r= 1.1278
143This should be close to max(r_c) in ps generation
144All pots = -2*Zval/r beyond r= 1.1278
145
146VLOCAL1: 99.0% of the norm of Vloc inside 34.126 Ry
147VLOCAL1: 99.9% of the norm of Vloc inside 77.774 Ry
148atom: Maximum radius for 4*pi*r*r*local-pseudopot. charge 1.37759
149atom: Maximum radius for r*vlocal+2*Zval: 1.18566
150GHOST: No ghost state for L = 0
151GHOST: No ghost state for L = 1
152GHOST: No ghost state for L = 2
153GHOST: No ghost state for L = 3
154
155KBgen: Kleinman-Bylander projectors:
156 l= 0 rc= 1.294105 el= -1.742414 Ekb= 9.135903 kbcos= 0.326910
157 l= 1 rc= 1.294105 el= -0.676589 Ekb= -8.124878 kbcos= -0.395047
158 l= 2 rc= 1.448233 el= 0.002386 Ekb= -2.039267 kbcos= -0.003484
159 l= 3 rc= 1.561052 el= 0.003508 Ekb= -0.799141 kbcos= -0.000344
160
161KBgen: Total number of Kleinman-Bylander projectors: 16
162atom: -------------------------------------------------------------------------
163
164atom: SANKEY-TYPE ORBITALS:
165atom: Selected multiple-zeta basis: split
166
167SPLIT: Orbitals with angular momentum L= 0
168
169SPLIT: Basis orbitals for state 2s
170
171SPLIT: PAO cut-off radius determined from an
172SPLIT: energy shift= 0.020000 Ry
173
174 izeta = 1
175 lambda = 1.000000
176 rc = 3.305093
177 energy = -1.723766
178 kinetic = 1.614911
179 potential(screened) = -3.338677
180 potential(ionic) = -11.304675
181
182 izeta = 2
183 rmatch = 2.510382
184 splitnorm = 0.150000
185 energy = -1.471299
186 kinetic = 2.446434
187 potential(screened) = -3.917732
188 potential(ionic) = -12.476133
189
190SPLIT: Orbitals with angular momentum L= 1
191
192SPLIT: Basis orbitals for state 2p
193
194SPLIT: PAO cut-off radius determined from an
195SPLIT: energy shift= 0.020000 Ry
196
197 izeta = 1
198 lambda = 1.000000
199 rc = 3.937239
200 energy = -0.658841
201 kinetic = 5.005986
202 potential(screened) = -5.664827
203 potential(ionic) = -13.452360
204
205 izeta = 2
206 rmatch = 2.541963
207 splitnorm = 0.150000
208 energy = -0.367441
209 kinetic = 7.530509
210 potential(screened) = -7.897949
211 potential(ionic) = -16.611953
212
213POLgen: Perturbative polarization orbital with L= 2
214
215POLgen: Polarization orbital for state 2p
216
217 izeta = 1
218 rc = 3.937239
219 energy = 2.398520
220 kinetic = 4.716729
221 potential(screened) = -2.318209
222 potential(ionic) = -8.603170
223atom: Total number of Sankey-type orbitals: 13
224
225atm_pop: Valence configuration (for local Pseudopot. screening):
226 2s( 2.00)
227 2p( 4.00)
228Vna: chval, zval: 6.00000 6.00000
229
230Vna: Cut-off radius for the neutral-atom potential: 3.937239
231
232atom: _________________________________________________________________________
233
234<basis_specs>
235===============================================================================
236H Z= 1 Mass= 1.0100 Charge= 0.17977+309
237Lmxo=0 Lmxkb= 2 BasisType=split Semic=F
238L=0 Nsemic=0 Cnfigmx=1
239 n=1 nzeta=2 polorb=1
240 splnorm: 0.15000
241 vcte: 0.0000
242 rinn: 0.0000
243 qcoe: 0.0000
244 qyuk: 0.0000
245 qwid: 0.10000E-01
246 rcs: 0.0000 0.0000
247 lambdas: 1.0000 1.0000
248-------------------------------------------------------------------------------
249L=0 Nkbl=1 erefs: 0.17977+309
250L=1 Nkbl=1 erefs: 0.17977+309
251L=2 Nkbl=1 erefs: 0.17977+309
252===============================================================================
253</basis_specs>
254
255atom: Called for H (Z = 1)
256
257read_vps: Pseudopotential generation method:
258read_vps: ATM3 Troullier-Martins
259Total valence charge: 1.00000
260
261xc_check: Exchange-correlation functional:
262xc_check: Ceperley-Alder
263V l=0 = -2*Zval/r beyond r= 1.2343
264V l=1 = -2*Zval/r beyond r= 1.2189
265V l=2 = -2*Zval/r beyond r= 1.2189
266All V_l potentials equal beyond r= 1.2343
267This should be close to max(r_c) in ps generation
268All pots = -2*Zval/r beyond r= 1.2343
269
270VLOCAL1: 99.0% of the norm of Vloc inside 28.493 Ry
271VLOCAL1: 99.9% of the norm of Vloc inside 64.935 Ry
272atom: Maximum radius for 4*pi*r*r*local-pseudopot. charge 1.45251
273atom: Maximum radius for r*vlocal+2*Zval: 1.21892
274GHOST: No ghost state for L = 0
275GHOST: No ghost state for L = 1
276GHOST: No ghost state for L = 2
277
278KBgen: Kleinman-Bylander projectors:
279 l= 0 rc= 1.364359 el= -0.467325 Ekb= -2.005361 kbcos= -0.336422
280 l= 1 rc= 1.434438 el= 0.001430 Ekb= -0.501708 kbcos= -0.021697
281 l= 2 rc= 1.470814 el= 0.002365 Ekb= -0.190555 kbcos= -0.002281
282
283KBgen: Total number of Kleinman-Bylander projectors: 9
284atom: -------------------------------------------------------------------------
285
286atom: SANKEY-TYPE ORBITALS:
287atom: Selected multiple-zeta basis: split
288
289SPLIT: Orbitals with angular momentum L= 0
290
291SPLIT: Basis orbitals for state 1s
292
293SPLIT: PAO cut-off radius determined from an
294SPLIT: energy shift= 0.020000 Ry
295
296 izeta = 1
297 lambda = 1.000000
298 rc = 4.828263
299 energy = -0.449375
300 kinetic = 0.929372
301 potential(screened) = -1.378747
302 potential(ionic) = -1.915047
303
304 izeta = 2
305 rmatch = 3.854947
306 splitnorm = 0.150000
307 energy = -0.336153
308 kinetic = 1.505294
309 potential(screened) = -1.841447
310 potential(ionic) = -2.413582
311
312POLgen: Perturbative polarization orbital with L= 1
313
314POLgen: Polarization orbital for state 1s
315
316 izeta = 1
317 rc = 4.828263
318 energy = 0.706972
319 kinetic = 1.396397
320 potential(screened) = -0.689424
321 potential(ionic) = -1.169792
322atom: Total number of Sankey-type orbitals: 5
323
324atm_pop: Valence configuration (for local Pseudopot. screening):
325 1s( 1.00)
326Vna: chval, zval: 1.00000 1.00000
327
328Vna: Cut-off radius for the neutral-atom potential: 4.828263
329
330atom: _________________________________________________________________________
331
332<basis_specs>
333===============================================================================
334Bessel Z=-100 Mass= 0.10000E+41 Charge= 0.17977+309
335Lmxo=1 Lmxkb=-1 BasisType=split Semic=F
336L=0 Nsemic=1 Cnfigmx=2
337 n=1 nzeta=1 polorb=0
338 splnorm: 0.15000
339 vcte: 0.0000
340 rinn: 0.0000
341 qcoe: 0.0000
342 qyuk: 0.0000
343 qwid: 0.10000E-01
344 rcs: 2.0000
345 lambdas: 1.0000
346 n=2 nzeta=1 polorb=0
347 splnorm: 0.15000
348 vcte: 0.0000
349 rinn: 0.0000
350 qcoe: 0.0000
351 qyuk: 0.0000
352 qwid: 0.10000E-01
353 rcs: 2.5000
354 lambdas: 1.0000
355L=1 Nsemic=0 Cnfigmx=3
356 n=1 nzeta=1 polorb=0
357 splnorm: 0.15000
358 vcte: 0.0000
359 rinn: 0.0000
360 qcoe: 0.0000
361 qyuk: 0.0000
362 qwid: 0.10000E-01
363 rcs: 3.5000
364 lambdas: 1.0000
365-------------------------------------------------------------------------------
366===============================================================================
367</basis_specs>
368
369atom: Called for Bessel (Z =-100) ( Floating Bessel functions )
370
371Bessel: floating Bessel functions with angular momentum L= 0
372
373Bessel: Basis orbitals for "state" 1s
374
375 izeta = 1
376 rc = 2.011274
377 energy = 2.439817
378
379Bessel: Basis orbitals for "state" 2s
380
381 izeta = 1
382 rc = 2.519390
383 energy = 1.554924
384
385Bessel: floating Bessel functions with angular momentum L= 1
386
387Bessel: Basis orbitals for "state" 3p
388
389 izeta = 1
390 rc = 3.487864
391 energy = 1.659713
392
393atom: Total number of floating Bessel orbitals: 5
394
395atom: _________________________________________________________________________
396
397<basis_specs>
398===============================================================================
399J Z=-100 Mass= 0.10000E+41 Charge= 0.17977+309
400Lmxo=1 Lmxkb=-1 BasisType=split Semic=F
401L=0 Nsemic=0 Cnfigmx=2
402 n=1 nzeta=7 polorb=0
403 splnorm: 0.15000
404 vcte: 0.0000
405 rinn: 0.0000
406 qcoe: 0.0000
407 qyuk: 0.0000
408 qwid: 0.10000E-01
409 rcs: 4.5000 4.5000 4.5000 4.5000
410 lambdas: 1.0000 1.0000 1.0000 1.0000
411L=1 Nsemic=0 Cnfigmx=2
412 n=1 nzeta=3 polorb=0
413 splnorm: 0.15000
414 vcte: 0.0000
415 rinn: 0.0000
416 qcoe: 0.0000
417 qyuk: 0.0000
418 qwid: 0.10000E-01
419 rcs: 4.5000 4.5000 5.0000
420 lambdas: 1.0000 1.0000 1.0000
421-------------------------------------------------------------------------------
422===============================================================================
423</basis_specs>
424
425atom: Called for J (Z =-100) ( Floating Bessel functions )
426
427Bessel: floating Bessel functions with angular momentum L= 0
428
429Bessel: Basis orbitals for "state" 2s
430
431 izeta = 1
432 rc = 4.479210
433 energy = 0.491923
434
435 izeta = 2
436 rc = 4.479210
437 energy = 1.967691
438
439 izeta = 3
440 rc = 4.479210
441 energy = 4.427304
442
443 izeta = 4
444 rc = 4.479210
445 energy = 7.870760
446
447 izeta = 5
448 rc = 4.479210
449 energy = 12.298054
450
451 izeta = 6
452 rc = 4.479210
453 energy = 17.709175
454
455 izeta = 7
456 rc = 4.479210
457 energy = 24.104105
458
459Bessel: floating Bessel functions with angular momentum L= 1
460
461Bessel: Basis orbitals for "state" 2p
462
463 izeta = 1
464 rc = 4.479210
465 energy = 1.006350
466
467 izeta = 2
468 rc = 4.479210
469 energy = 2.974558
470
471 izeta = 3
472 rc = 5.075940
473 energy = 4.614751
474
475atom: Total number of floating Bessel orbitals: 16
476
477atom: _________________________________________________________________________
478
479prinput: Basis input ----------------------------------------------------------
480
481PAO.BasisType split
482
483%block ChemicalSpeciesLabel
484 1 8 O # Species index, atomic number, species label
485 2 1 H # Species index, atomic number, species label
486 3 -100 Bessel # Species index, atomic number, species label
487 4 -100 J # Species index, atomic number, species label
488%endblock ChemicalSpeciesLabel
489
490%block PAO.Basis # Define Basis set
491O 2 # Species label, number of l-shells
492 n=2 0 2 # n, l, Nzeta
493 3.305 2.510
494 1.000 1.000
495 n=2 1 2 P 1 # n, l, Nzeta, Polarization, NzetaPol
496 3.937 2.542
497 1.000 1.000
498H 1 # Species label, number of l-shells
499 n=1 0 2 P 1 # n, l, Nzeta, Polarization, NzetaPol
500 4.828 3.855
501 1.000 1.000
502Bessel 3 # Species label, number of l-shells
503 n=1 0 1 # n, l, Nzeta
504 2.011
505 1.000
506 n=2 0 1 # n, l, Nzeta
507 2.519
508 1.000
509 n=3 1 1 # n, l, Nzeta
510 3.488
511 1.000
512J 2 # Species label, number of l-shells
513 n=2 0 7 # n, l, Nzeta
514 4.479 4.479 4.479 4.479 4.479 4.479 4.479
515 1.000 1.000 1.000 1.000 1.000 1.000 1.000
516 n=2 1 3 # n, l, Nzeta
517 4.479 4.479 5.076
518 1.000 1.000 1.000
519%endblock PAO.Basis
520
521prinput: ----------------------------------------------------------------------
522
523Dumping basis to NetCDF file O.ion.nc
524Dumping basis to NetCDF file H.ion.nc
525Dumping basis to NetCDF file Bessel.ion.nc
526Dumping basis to NetCDF file J.ion.nc
527coor: Atomic-coordinates input format = Cartesian coordinates
528coor: (in Angstroms)
529
530siesta: Atomic coordinates (Bohr) and species
531siesta: 0.00000 0.00000 0.00000 1 1
532siesta: 1.43052 1.10738 0.00000 2 2
533siesta: -1.43052 1.10738 0.00000 2 3
534siesta: 0.71526 0.55369 0.00000 3 4
535siesta: -0.71526 0.55369 0.00000 3 5
536siesta: 0.71526 0.55369 0.00000 4 6
537siesta: -0.71526 0.55369 0.00000 4 7
538
539siesta: Automatic unit cell vectors (Ang):
540siesta: 7.286412 0.000000 0.000000
541siesta: 0.000000 6.087484 0.000000
542siesta: 0.000000 0.000000 5.909356
543
544siesta: System type = molecule
545
546initatomlists: Number of atoms, orbitals, and projectors: 7 65 34
547
548siesta: ******************** Simulation parameters ****************************
549siesta:
550siesta: The following are some of the parameters of the simulation.
551siesta: A complete list of the parameters used, including default values,
552siesta: can be found in file out.fdf
553siesta:
554redata: Non-Collinear-spin run = F
555redata: SpinPolarized (Up/Down) run = F
556redata: Number of spin components = 1
557redata: Long output = F
558redata: Number of Atomic Species = 4
559redata: Charge density info will appear in .RHO file
560redata: Write Mulliken Pop. = NO
561redata: Matel table size (NRTAB) = 1024
562redata: Mesh Cutoff = 100.0000 Ry
563redata: Net charge of the system = 0.0000 |e|
564redata: Min. number of SCF Iter = 0
565redata: Max. number of SCF Iter = 50
566redata: Mix DM or H after convergence = F
567redata: Recompute H after scf cycle = F
568redata: Mixing is linear
569redata: Mix DM in first SCF step ? = F
570redata: Write Pulay info on disk? = F
571redata: Discard 1st Pulay DM after kick = F
572redata: New DM Mixing Weight = 0.2500
573redata: New DM Occupancy tolerance = 0.000000000001
574redata: No kicks to SCF
575redata: DM Mixing Weight for Kicks = 0.5000
576redata: DM Tolerance for SCF = 0.000100
577redata: Require (free) Energy convergence in SCF = F
578redata: DM (free)Energy tolerance for SCF = 0.000010 eV
579redata: Require Harris convergence for SCF = F
580redata: DM Harris energy tolerance for SCF = 0.000010 eV
581redata: Using Saved Data (generic) = F
582redata: Use continuation files for DM = F
583redata: Neglect nonoverlap interactions = F
584redata: Method of Calculation = Diagonalization
585redata: Divide and Conquer = T
586redata: Electronic Temperature = 0.0019 Ry
587redata: Fix the spin of the system = F
588redata: Dynamics option = Single-point calculation
589redata: ***********************************************************************
590Total number of electrons: 8.000000
591Total ionic charge: 8.000000
592
593* ProcessorY, Blocksize: 1 24
594
595
596* Orbital distribution balance (max,min): 65 65
597
598 Kpoints in: 1 . Kpoints trimmed: 1
599
600siesta: k-grid: Number of k-points = 1
601siesta: k-grid: Cutoff (effective) = 2.955 Ang
602siesta: k-grid: Supercell and displacements
603siesta: k-grid: 1 0 0 0.000
604siesta: k-grid: 0 1 0 0.000
605siesta: k-grid: 0 0 1 0.000
606
607 ====================================
608 Single-point calculation
609 ====================================
610
611outcell: Unit cell vectors (Ang):
612 7.286412 0.000000 0.000000
613 0.000000 6.087484 0.000000
614 0.000000 0.000000 5.909356
615
616outcell: Cell vector modules (Ang) : 7.286412 6.087484 5.909356
617outcell: Cell angles (23,13,12) (deg): 90.0000 90.0000 90.0000
618outcell: Cell volume (Ang**3) : 262.1149
619New_DM. Step: 1
620Initializing Density Matrix...
621New grid distribution: 1
622 1 1: 24 1: 20 1: 18
623
624InitMesh: MESH = 48 x 40 x 36 = 69120
625InitMesh: (bp) = 24 x 20 x 18 = 8640
626InitMesh: Mesh cutoff (required, used) = 100.000 102.571 Ry
627ExtMesh (bp) on 0 = 60 x 56 x 54 = 181440
628PhiOnMesh: Number of (b)points on node 0 = 8640
629PhiOnMesh: nlist on node 0 = 119152
630
631stepf: Fermi-Dirac step function
632
633siesta: Program's energy decomposition (eV):
634siesta: Ebs = -124.892071
635siesta: Eions = 815.854478
636siesta: Ena = 175.155695
637siesta: Ekin = 341.667406
638siesta: Enl = -52.736859
639siesta: DEna = -0.000003
640siesta: DUscf = 0.000000
641siesta: DUext = 0.000000
642siesta: Exc = -109.898170
643siesta: eta*DQ = 0.000000
644siesta: Emadel = 0.000000
645siesta: Emeta = 0.000000
646siesta: Emolmec = 0.000000
647siesta: Ekinion = 0.000000
648siesta: Eharris = -467.297614
649siesta: Etot = -461.666410
650siesta: FreeEng = -461.666410
651
652 scf: iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV)
653 scf: 1 -467.2976 -461.6664 -461.6664 1.77421 -8.4794
654timer: Routine,Calls,Time,% = IterSCF 1 0.073 0.70
655 scf: 2 -467.9568 -465.9473 -465.9473 0.17421 0.7795
656 scf: 3 -466.7456 -466.2075 -466.2075 0.03742 -1.0131
657 scf: 4 -466.6648 -466.3262 -466.3262 0.02266 -1.4679
658 scf: 5 -466.6584 -466.4096 -466.4096 0.02027 -1.5893
659 scf: 6 -466.6578 -466.4717 -466.4717 0.01648 -1.6214
660 scf: 7 -466.6578 -466.5182 -466.5182 0.01273 -1.6293
661 scf: 8 -466.6577 -466.5531 -466.5531 0.00965 -1.6311
662 scf: 9 -466.6577 -466.5793 -466.5793 0.00727 -1.6314
663 scf: 10 -466.6577 -466.5989 -466.5989 0.00545 -1.6314
664 scf: 11 -466.6577 -466.6136 -466.6136 0.00409 -1.6315
665 scf: 12 -466.6577 -466.6246 -466.6246 0.00306 -1.6316
666 scf: 13 -466.6577 -466.6329 -466.6329 0.00229 -1.6317
667 scf: 14 -466.6577 -466.6391 -466.6391 0.00171 -1.6318
668 scf: 15 -466.6577 -466.6438 -466.6438 0.00128 -1.6319
669 scf: 16 -466.6577 -466.6472 -466.6472 0.00096 -1.6319
670 scf: 17 -466.6577 -466.6499 -466.6499 0.00072 -1.6320
671 scf: 18 -466.6577 -466.6518 -466.6518 0.00054 -1.6320
672 scf: 19 -466.6577 -466.6533 -466.6533 0.00040 -1.6321
673 scf: 20 -466.6577 -466.6544 -466.6544 0.00030 -1.6321
674 scf: 21 -466.6577 -466.6552 -466.6552 0.00022 -1.6321
675 scf: 22 -466.6577 -466.6559 -466.6559 0.00017 -1.6321
676 scf: 23 -466.6577 -466.6563 -466.6563 0.00012 -1.6322
677 scf: 24 -466.6577 -466.6567 -466.6567 0.00009 -1.6322
678
679SCF Convergence by dMax criterion
680max |DM_out - DM_in|: 0.00009343
681SCF cycle converged after 24 iterations
682
683Using DM_out to compute the final energy and forces
684
685siesta: E_KS(eV) = -466.6577
686
687siesta: E_KS - E_eggbox = -466.6577
688
689siesta: Atomic forces (eV/Ang):
690----------------------------------------
691 Tot -0.000000 -0.032764 -0.000000
692----------------------------------------
693 Max 0.376555
694 Res 0.167768 sqrt( Sum f_i^2 / 3N )
695----------------------------------------
696 Max 0.376555 constrained
697
698Stress-tensor-Voigt (kbar): -3.42 -0.89 -1.16 0.00 0.00 0.00
699(Free)E + p*V (eV/cell) -466.3595
700Target enthalpy (eV/cell) -466.6577
701
702siesta: Program's energy decomposition (eV):
703siesta: Ebs = -107.092904
704siesta: Eions = 815.854478
705siesta: Ena = 175.155695
706siesta: Ekin = 351.559238
707siesta: Enl = -63.040421
708siesta: DEna = -2.430290
709siesta: DUscf = 0.775212
710siesta: DUext = 0.000000
711siesta: Exc = -112.822678
712siesta: eta*DQ = 0.000000
713siesta: Emadel = 0.000000
714siesta: Emeta = 0.000000
715siesta: Emolmec = 0.000000
716siesta: Ekinion = 0.000000
717siesta: Eharris = -466.657724
718siesta: Etot = -466.657724
719siesta: FreeEng = -466.657724
720
721siesta: Final energy (eV):
722siesta: Band Struct. = -107.092904
723siesta: Kinetic = 351.559238
724siesta: Hartree = 388.214304
725siesta: Ext. field = 0.000000
726siesta: Exch.-corr. = -112.822678
727siesta: Ion-electron = -1087.255665
728siesta: Ion-ion = -6.352922
729siesta: Ekinion = 0.000000
730siesta: Total = -466.657724
731
732siesta: Atomic forces (eV/Ang):
733siesta: 1 0.000000 -0.018600 0.000000
734siesta: 2 0.376555 0.275107 -0.000000
735siesta: 3 -0.376555 0.275107 -0.000000
736siesta: 4 -0.001135 -0.003186 0.000000
737siesta: 5 0.001135 -0.003186 0.000000
738siesta: 6 -0.005249 -0.279003 0.000000
739siesta: 7 0.005249 -0.279003 0.000000
740siesta: ----------------------------------------
741siesta: Tot -0.000000 -0.032764 -0.000000
742
743siesta: Stress tensor (static) (eV/Ang**3):
744siesta: -0.002134 0.000000 -0.000000
745siesta: 0.000000 -0.000556 -0.000000
746siesta: 0.000000 0.000000 -0.000724
747
748siesta: Cell volume = 262.114919 Ang**3
749
750siesta: Pressure (static):
751siesta: Solid Molecule Units
752siesta: 0.00001239 0.00000224 Ry/Bohr**3
753siesta: 0.00113783 0.00020528 eV/Ang**3
754siesta: 1.82302327 0.32889750 kBar
755(Free)E+ p_basis*V_orbitals = -464.513364
756(Free)Eharris+ p_basis*V_orbitals = -464.513364
757
758siesta: Electric dipole (a.u.) = 0.000000 0.657893 0.000000
759siesta: Electric dipole (Debye) = 0.000000 1.672200 0.000000
760
761timer: Elapsed wall time (sec) = 11.961
762timer: CPU execution times (sec):
763
764Routine Calls Time/call Tot.time %
765siesta 1 11.909 11.909 100.00
766Setup 1 0.365 0.365 3.06
767bands 1 0.000 0.000 0.00
768KSV_init 1 0.000 0.000 0.00
769IterGeom 1 11.539 11.539 96.90
770state_init 1 2.127 2.127 17.86
771hsparse 1 0.002 0.002 0.02
772overlap 1 2.123 2.123 17.83
773Setup_H0 1 7.766 7.766 65.22
774naefs 2 0.000 0.001 0.00
775MolMec 2 0.000 0.000 0.00
776kinefsm 2 0.990 1.979 16.62
777nlefsm 2 2.848 5.695 47.82
778DHSCF_Init 1 0.097 0.097 0.81
779DHSCF1 1 0.007 0.007 0.06
780INITMESH 1 0.000 0.000 0.00
781DHSCF2 1 0.090 0.090 0.75
782REMESH 1 0.008 0.008 0.07
783REORD 58 0.000 0.006 0.05
784PHION 1 0.067 0.067 0.56
785COMM_BSC 53 0.000 0.006 0.05
786POISON 26 0.007 0.174 1.46
787fft 52 0.003 0.152 1.27
788IterSCF 24 0.061 1.464 12.29
789setup_H 24 0.059 1.420 11.92
790DHSCF 25 0.064 1.592 13.37
791DHSCF3 25 0.059 1.475 12.38
792rhoofd 25 0.029 0.714 6.00
793cellXC 25 0.004 0.108 0.90
794vmat 25 0.018 0.449 3.77
795compute_dm 24 0.001 0.028 0.24
796diagon 24 0.001 0.027 0.22
797r-eigvec 24 0.001 0.025 0.21
798r-buildHS 24 0.000 0.000 0.00
799rdiag 24 0.001 0.025 0.21
800rdiag1 24 0.000 0.001 0.01
801rdiag2 24 0.000 0.002 0.02
802rdiag3 24 0.001 0.018 0.15
803rdiag4 24 0.000 0.003 0.02
804r-buildD 24 0.000 0.001 0.01
805MIXER 23 0.000 0.000 0.00
806WriteDM 24 0.000 0.010 0.09
807PostSCF 1 0.180 0.180 1.51
808DHSCF4 1 0.118 0.118 0.99
809dfscf 1 0.112 0.112 0.94
810overfsm 1 0.001 0.001 0.01
811state_analysis 1 0.002 0.002 0.02
812siesta_move 1 0.000 0.000 0.00
813siesta_analysis 1 0.003 0.003 0.02
814optical 1 0.000 0.000 0.00
815
816>> End of run: 17-JAN-2019 15:47:57
817Job completed
0818
=== added directory 'Tests/bessel-rich'
=== added file 'Tests/bessel-rich/bessel-rich.fdf'
--- Tests/bessel-rich/bessel-rich.fdf 1970-01-01 00:00:00 +0000
+++ Tests/bessel-rich/bessel-rich.fdf 2019-01-30 13:27:17 +0000
@@ -0,0 +1,45 @@
1#
2# Needs new feature: handling of fewer rc's than nzetas in PAO.Basis block
3#
4write-ion-plot-files T
5#
6SystemName Water molecule with various Bessel Orbitals
7SystemLabel bessel-rich
8NumberofAtoms 7
9NumberOfSpecies 4
10%block ChemicalSpeciesLabel
11 1 8 O # Species index, atomic number, species label
12 2 1 H
13 3 -100 Bessel
14 4 -100 J
15%endblock ChemicalSpeciesLabel
16
17AtomicCoordinatesFormat Ang
18%block AtomicCoordinatesAndAtomicSpecies
19 0.000 0.000 0.000 1
20 0.757 0.586 0.000 2
21-0.757 0.586 0.000 2
22 0.3785 0.293 0.000 3
23-0.3785 0.293 0.000 3
24 0.3785 0.293 0.000 4
25-0.3785 0.293 0.000 4
26%endblock AtomicCoordinatesAndAtomicSpecies
27
28%block PAO.Basis
29Bessel 3
30 n=1 0 1
31 2.0
32 1.0
33 n=2 0 1
34 2.5
35 1.0
36 n=3 1 1
37 3.5
38 1.0
39J 2 # l-shells
40n=2 0 7 # Note new feature: fewer rc's than zetas
41 4.5
42n=2 1 3
43 4.5 4.5 5.0
44%endblock PAO.Basis
45
0\ No newline at end of file46\ No newline at end of file
147
=== added file 'Tests/bessel-rich/bessel-rich.pseudos'
--- Tests/bessel-rich/bessel-rich.pseudos 1970-01-01 00:00:00 +0000
+++ Tests/bessel-rich/bessel-rich.pseudos 2019-01-30 13:27:17 +0000
@@ -0,0 +1,1 @@
1H O
02
=== added file 'Tests/bessel-rich/makefile'
--- Tests/bessel-rich/makefile 1970-01-01 00:00:00 +0000
+++ Tests/bessel-rich/makefile 2019-01-30 13:27:17 +0000
@@ -0,0 +1,6 @@
1#
2# Single-test makefile
3#
4name=bessel-rich
5#
6include ../test.mk
07
=== modified file 'version.info'
--- version.info 2018-12-19 13:52:34 +0000
+++ version.info 2019-01-30 13:27:17 +0000
@@ -1,1 +1,4 @@
1siesta-4.0--5891siesta-4.0--589--n-fix-5--semic-5
2
3
4

Subscribers

People subscribed via source and target branches