From d27e6b355a1f526cef853422f9e7a355998bd96c Mon Sep 17 00:00:00 2001 From: Tammo van der Heide Date: Thu, 4 Jul 2024 15:13:33 +0200 Subject: [PATCH] Hand over WS compression to slateratom --- sktools/src/sktools/calculators/slateratom.py | 75 +++++++++++---- sktools/src/sktools/common.py | 2 + sktools/src/sktools/compressions.py | 73 +++++++++++++++ slateratom/lib/CMakeLists.txt | 1 + slateratom/lib/confinement.f90 | 91 +++++++++++++++++++ slateratom/lib/core_overlap.f90 | 61 +------------ slateratom/lib/globals.f90 | 29 +++++- slateratom/lib/input.f90 | 79 +++++++++++----- slateratom/prog/main.f90 | 22 +++-- 9 files changed, 319 insertions(+), 114 deletions(-) create mode 100644 slateratom/lib/confinement.f90 diff --git a/sktools/src/sktools/calculators/slateratom.py b/sktools/src/sktools/calculators/slateratom.py index 5659ddb6..f429d284 100644 --- a/sktools/src/sktools/calculators/slateratom.py +++ b/sktools/src/sktools/calculators/slateratom.py @@ -10,8 +10,9 @@ import sktools.hsd.converter as conv import sktools.common as sc from sktools.taggedfile import TaggedFile -import sktools.compressions +import sktools.compressions as skcomp import sktools.radial_grid as oc +import sktools.xcfunctionals as xc LOGGER = logging.getLogger('slateratom') @@ -134,7 +135,7 @@ class SlateratomInput: functional : str DFT functional type ('lda', 'pbe', 'blyp', 'lcpbe', 'lcbnl') compressions : list - List of PowerCompression objects. Either empty (no compression applied) + List of Compression objects. Either empty (no compression applied) or has a compression object for every angular momentum of the atom. settings : SlaterAtom Further detailed settings of the program. @@ -193,16 +194,26 @@ def __init__(self, settings, atomconfig, functional, compressions): if compressions is None: compressions = [] for comp in compressions: - if not isinstance(comp, sktools.compressions.PowerCompression): - msg = "Invalid compressiont type {} for slateratom".format( + if not isinstance(comp, (skcomp.PowerCompression, + skcomp.WoodsSaxonCompression)): + msg = "Invalid compression type {} for slateratom".format( comp.__class__.__name__) raise sc.SkgenException(msg) maxang = atomconfig.maxang ncompr = len(compressions) - if ncompr and ncompr != maxang + 1: + if ncompr > 0 and ncompr != maxang + 1: msg = "Invalid number of compressions" \ "(expected {:d}, got {:d})".format(maxang + 1, ncompr) raise sc.SkgenException(msg) + + # block different compression types of different shells for now + if ncompr > 0: + compids = [comp.compid for comp in compressions] + if not len(set(compids)) == 1: + msg = "At the moment, shells may only be compressed by the " \ + + "same type of potential." + raise sc.SkgenException(msg) + self._compressions = compressions myrelativistics = sc.RELATIVISTICS_NONE, sc.RELATIVISTICS_ZORA if atomconfig.relativistics not in myrelativistics: @@ -213,7 +224,7 @@ def __init__(self, settings, atomconfig, functional, compressions): def isXCFunctionalSupported(self, functional): '''Checks if the given xc-functional is supported by the calculator, in particular: checks if AVAILABLE_FUNCTIONALS intersect with - sktools.xcfunctionals.XCFUNCTIONALS + xc.XCFUNCTIONALS Args: functional: xc-functional, defined in xcfunctionals.py @@ -224,8 +235,8 @@ def isXCFunctionalSupported(self, functional): tmp = [] for xx in SUPPORTED_FUNCTIONALS: - if xx in sktools.xcfunctionals.XCFUNCTIONALS: - tmp.append(sktools.xcfunctionals.XCFUNCTIONALS[xx]) + if xx in xc.XCFUNCTIONALS: + tmp.append(xc.XCFUNCTIONALS[xx]) return bool(functional.__class__ in tmp) @@ -297,14 +308,42 @@ def write(self, workdir): # Compressions if len(self._compressions) == 0: - out += ["{:g} {:d} \t\t{:s} Compr. radius and power ({:s})".format( - 1e30, 0, self._COMMENT, sc.ANGMOM_TO_SHELL[ll]) + # no compression for all shells + out += ["{:d} \t\t\t{:s} Compression ID".format( + skcomp.SUPPORTED_COMPRESSIONS['nocompression'], + self._COMMENT) + " ({:s})".format(sc.ANGMOM_TO_SHELL[ll]) for ll in range(maxang + 1)] else: - out += ["{:g} {:g} \t\t{:s} Compr. radius and power ({:s})".format( - compr.radius, compr.power, self._COMMENT, - sc.ANGMOM_TO_SHELL[ll]) - for ll, compr in enumerate(self._compressions)] + # define the type of compression for each shell + for ll, compr in enumerate(self._compressions): + if compr.compid == \ + skcomp.SUPPORTED_COMPRESSIONS['powercompression']: + out += ["{:d} \t\t{:s} Compression ID".format( + skcomp.SUPPORTED_COMPRESSIONS['powercompression'], + self._COMMENT) \ + + " ({:s})".format(sc.ANGMOM_TO_SHELL[ll])] + elif compr.compid == \ + skcomp.SUPPORTED_COMPRESSIONS['woodssaxoncompression']: + out += ["{:d} \t\t{:s} Compression ID".format( + skcomp.SUPPORTED_COMPRESSIONS['woodssaxoncompression'], + self._COMMENT) \ + + " ({:s})".format(sc.ANGMOM_TO_SHELL[ll])] + else: + msg = 'Invalid compression type.' + raise sc.SkgenException(msg) + # provide the compression parametrization for each shell + for ll, compr in enumerate(self._compressions): + if compr.compid == \ + skcomp.SUPPORTED_COMPRESSIONS['powercompression']: + out += ["{:g} {:g} \t\t{:s} Compr. radius and power ({:s})" + .format(compr.radius, compr.power, self._COMMENT, + sc.ANGMOM_TO_SHELL[ll])] + elif compr.compid == \ + skcomp.SUPPORTED_COMPRESSIONS['woodssaxoncompression']: + out += ["{:g} {:g} {:g} \t\t{:s}".format( + compr.onset, compr.cutoff, compr.vmax, self._COMMENT) \ + + " Compr. onset, cutoff and vmax ({:s})" + .format(sc.ANGMOM_TO_SHELL[ll])] out += ["{:d} \t\t\t{:s} nr. of occupied shells ({:s})".format( len(occ), self._COMMENT, sc.ANGMOM_TO_SHELL[ll]) @@ -313,9 +352,9 @@ def write(self, workdir): # STO powers and exponents exponents = self._settings.exponents maxpowers = self._settings.maxpowers - out += ["{:d} {:d} \t\t\t{:s} nr. of exponents, max. power ({:s})".format( - len(exponents[ll]), maxpowers[ll], self._COMMENT, - sc.ANGMOM_TO_SHELL[ll]) + out += ["{:d} {:d} \t\t\t{:s} nr. of exponents, max. power ({:s})" + .format(len(exponents[ll]), maxpowers[ll], self._COMMENT, + sc.ANGMOM_TO_SHELL[ll]) for ll in range(maxang + 1)] out.append("{:s} \t\t{:s} automatic exponent generation".format( self._LOGICALSTRS[False], self._COMMENT)) @@ -326,7 +365,7 @@ def write(self, workdir): out.append("{:s} \t\t{:s} write eigenvectors".format( self._LOGICALSTRS[False], self._COMMENT)) - out.append("{} {:g} \t\t{:s} broyden mixer, mixing factor".format( + out.append("{} {:g} \t\t\t{:s} broyden mixer, mixing factor".format( 2, 0.1, self._COMMENT)) # Occupations diff --git a/sktools/src/sktools/common.py b/sktools/src/sktools/common.py index 7c97d78d..8e3aa313 100644 --- a/sktools/src/sktools/common.py +++ b/sktools/src/sktools/common.py @@ -528,6 +528,8 @@ def is_shelf_file_matching(shelf_file, mydict): db = shelve.open(shelf_file, 'r') except dbm.error: return False + if not type(db) is type(mydict): + return False match = True for key in mydict: match = key in db and db[key] == mydict[key] diff --git a/sktools/src/sktools/compressions.py b/sktools/src/sktools/compressions.py index 62a1eb05..ec54c48b 100644 --- a/sktools/src/sktools/compressions.py +++ b/sktools/src/sktools/compressions.py @@ -39,6 +39,8 @@ def fromhsd(cls, root, query): node=child) myself = cls() + myself.compid = SUPPORTED_COMPRESSIONS[ + myself.__class__.__name__.lower()] myself.power = power myself.radius = radius @@ -62,9 +64,80 @@ def __eq__(self, other): return power_ok and radius_ok +class WoodsSaxonCompression(sc.ClassDict): + '''Compression by the Woods Saxon potential. + + Attributes + ---------- + onset : float + Onset radius of the compression + cutoff : float + Cutoff radius of the compression + vmax : float + Potential well depth/height + ''' + + @classmethod + def fromhsd(cls, root, query): + '''Creates instance from a HSD-node and with given query object.''' + + onset, child = query.getvalue(root, 'onset', conv.float0, + returnchild=True) + if onset < 0.0: + msg = 'Invalid onset radius {:f}'.format(onset) + raise hsd.HSDInvalidTagValueException(msg=msg, node=child) + + cutoff, child = query.getvalue(root, 'cutoff', conv.float0, + returnchild=True) + if cutoff <= onset: + msg = 'Invalid cutoff radius {:f}'.format(cutoff) + raise hsd.HSDInvalidTagValueException(msg=msg, node=child) + + vmax, child = query.getvalue( + root, 'vmax', conv.float0, defvalue=100.0, returnchild=True) + if vmax <= 0.0: + msg = 'Invalid potential well height {:f}'.format(vmax) + raise hsd.HSDInvalidTagValueException(msg=msg, node=child) + + myself = cls() + myself.compid = SUPPORTED_COMPRESSIONS[ + myself.__class__.__name__.lower()] + myself.onset = onset + myself.cutoff = cutoff + myself.vmax = vmax + + return myself + + + def tohsd(self, root, query, parentname=None): + '''''' + + if parentname is None: + mynode = root + else: + mynode = query.setchild(root, 'WoodsSaxonCompression') + + query.setchildvalue(mynode, 'onset', conv.float0, self.onset) + query.setchildvalue(mynode, 'cutoff', conv.float0, self.cutoff) + query.setchildvalue(mynode, 'vmax', conv.float0, self.vmax) + + + def __eq__(self, other): + onset_ok = abs(self.onset - other.onset) < 1e-8 + cutoff_ok = abs(self.cutoff - other.cutoff) < 1e-8 + vmax_ok = abs(self.vmax - other.vmax) < 1e-8 + return onset_ok and cutoff_ok and vmax_ok + + # Registered compressions with corresponding hsd name as key COMPRESSIONS = { 'powercompression': PowerCompression, + 'woodssaxoncompression': WoodsSaxonCompression, +} +SUPPORTED_COMPRESSIONS = { + 'nocompression': 0, + 'powercompression': 1, + 'woodssaxoncompression': 2, } diff --git a/slateratom/lib/CMakeLists.txt b/slateratom/lib/CMakeLists.txt index 78399778..9dea55e3 100644 --- a/slateratom/lib/CMakeLists.txt +++ b/slateratom/lib/CMakeLists.txt @@ -9,6 +9,7 @@ list(APPEND fypp_flags -I${projectdir}/common/include -DRELEASE="'${RELEASE}'") set(sources-f90 avgpot.f90 blas.f90 + confinement.f90 core_overlap.f90 coulomb_hfex.f90 coulomb_potential.f90 diff --git a/slateratom/lib/confinement.f90 b/slateratom/lib/confinement.f90 new file mode 100644 index 00000000..7822b017 --- /dev/null +++ b/slateratom/lib/confinement.f90 @@ -0,0 +1,91 @@ +!> Module that builds up various supervectors. +module confinement + + use common_accuracy, only : dp + use utilities, only : fak + use globals, only : TConf + use core_overlap, only : v + + implicit none + private + + public :: confType + public :: getConf_power + + + !> Enumerator for type of confinement potential. + type :: TConfEnum + + !> no compression + integer :: none = 0 + + !> power compression + integer :: power = 1 + + !> Woods-Saxon compression + integer :: ws = 2 + + end type TConfEnum + + !> Container for enumerated types of confinement potentials. + type(TConfEnum), parameter :: confType = TConfEnum() + + +contains + + !> Calculates analytic matrix elements of confining potential. + !! No checking for power, e.g. power==0 or power<0 etc. ! + pure subroutine getConf_power(vconf, max_l, num_alpha, alpha, poly_order, conf) + + !> confinement supervector + real(dp), intent(out) :: vconf(0:,:,:) + + !> maximum angular momentum + integer, intent(in) :: max_l + + !> number of exponents in each shell + integer, intent(in) :: num_alpha(0:) + + !> basis exponents + real(dp), intent(in) :: alpha(0:,:) + + !> highest polynomial order + l in each shell + integer, intent(in) :: poly_order(0:) + + !> confinement radii (power compression) + type(TConf), intent(in) :: conf + + !! temporary storage + real(dp) :: alpha1 + + !! auxiliary variables + integer :: ii, jj, kk, ll, mm, nn, oo, nlp, nlq + + vconf(:,:,:) = 0.0_dp + + do ii = 0, max_l + if (conf%power(ii) > 1.0e-06_dp) then + nn = 0 + do jj = 1, num_alpha(ii) + do ll = 1, poly_order(ii) + nn = nn + 1 + oo = 0 + nlp = ll + ii + do kk = 1, num_alpha(ii) + alpha1 = 0.5_dp * (alpha(ii, jj) + alpha(ii, kk)) + do mm = 1, poly_order(ii) + oo = oo + 1 + nlq = mm + ii + vconf(ii, nn, oo) = 1.0_dp / sqrt(v(alpha(ii, jj), 2 * nlp)& + & * v(alpha(ii, kk), 2 * nlq)) / (conf%r0(ii) * 2.0_dp)**conf%power(ii)& + & * v(alpha1, nlp + nlq + conf%power(ii)) + end do + end do + end do + end do + end if + end do + + end subroutine getConf_power + +end module confinement diff --git a/slateratom/lib/core_overlap.f90 b/slateratom/lib/core_overlap.f90 index 258386c6..91e9d523 100644 --- a/slateratom/lib/core_overlap.f90 +++ b/slateratom/lib/core_overlap.f90 @@ -7,7 +7,7 @@ module core_overlap implicit none private - public :: overlap, kinetic, nuclear, moments, v, confinement + public :: overlap, kinetic, nuclear, moments, v interface v @@ -192,65 +192,6 @@ pure subroutine kinetic(tt, max_l, num_alpha, alpha, poly_order) end subroutine kinetic - !> Calculates analytic matrix elements of confining potential. - !! No checking for power, e.g. power==0 or power<0 etc. ! - pure subroutine confinement(vconf, max_l, num_alpha, alpha, poly_order, conf_r0, conf_power) - - !> confinement supervector - real(dp), intent(out) :: vconf(0:,:,:) - - !> maximum angular momentum - integer, intent(in) :: max_l - - !> number of exponents in each shell - integer, intent(in) :: num_alpha(0:) - - !> basis exponents - real(dp), intent(in) :: alpha(0:,:) - - !> highest polynomial order + l in each shell - integer, intent(in) :: poly_order(0:) - - !> confinement radii - real(dp), intent(in) :: conf_r0(0:) - - !> power of confinement - real(dp), intent(in) :: conf_power(0:) - - !> temporary storage - real(dp) :: alpha1 - - !> auxiliary variables - integer :: ii, jj, kk, ll, mm, nn, oo, nlp, nlq - - vconf(:,:,:) = 0.0_dp - - do ii = 0, max_l - if (conf_power(ii) > 1.0e-06_dp) then - nn = 0 - do jj = 1, num_alpha(ii) - do ll = 1, poly_order(ii) - nn = nn + 1 - oo = 0 - nlp = ll + ii - do kk = 1, num_alpha(ii) - alpha1 = 0.5_dp * (alpha(ii, jj) + alpha(ii, kk)) - do mm = 1, poly_order(ii) - oo = oo + 1 - nlq = mm + ii - vconf(ii, nn, oo) = 1.0_dp / sqrt(v(alpha(ii, jj), 2 * nlp)& - & * v(alpha(ii, kk), 2 * nlq)) / (conf_r0(ii) * 2.0_dp)**conf_power(ii)& - & * v(alpha1, nlp + nlq + conf_power(ii)) - end do - end do - end do - end do - end if - end do - - end subroutine confinement - - !> Calculates arbitrary moments of electron distribution, e.g. expectation values of , !! etc.; this is implemented analytically for arbitrary powers. pure subroutine moments(moment, max_l, num_alpha, alpha, poly_order, cof, power) diff --git a/slateratom/lib/globals.f90 b/slateratom/lib/globals.f90 index 6687874d..5700fb4f 100644 --- a/slateratom/lib/globals.f90 +++ b/slateratom/lib/globals.f90 @@ -8,11 +8,32 @@ module globals implicit none - !> confinement radii - real(dp) :: conf_r0(0:4) + !> Confinement potential structure. + type :: TConf - !> power of confinement - real(dp) :: conf_power(0:4) + !> type of confinement potential (0: none, 1: power, 2: Woods-Saxon) + !! at the moment different shells are compressed by the same type of potential + integer :: type = 0 + + !> confinement radii (power compression) + real(dp) :: r0(0:4) + + !> power of confinement (power compression) + real(dp) :: power(0:4) + + !> onset radii (Woods-Saxon compression) + real(dp) :: onset(0:4) + + !> cutoff of confinement (Woods-Saxon compression) + real(dp) :: cutoff(0:4) + + !> potential well height of confinement (Woods-Saxon compression) + real(dp) :: vmax(0:4) + + end type TConf + + !> confinement potential + type(TConf) :: conf !> basis exponents real(dp) :: alpha(0:4, 10) diff --git a/slateratom/lib/input.f90 b/slateratom/lib/input.f90 index 51104403..a588e387 100644 --- a/slateratom/lib/input.f90 +++ b/slateratom/lib/input.f90 @@ -4,6 +4,8 @@ module input use common_accuracy, only : dp use common_poisson, only : TBeckeGridParams + use globals, only : TConf + use confinement, only : confType use xcfunctionals, only : xcFunctional implicit none @@ -16,9 +18,9 @@ module input !> Reads in all properties, except for occupation numbers. subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min_alpha,& - & max_alpha, num_alpha, tAutoAlphas, alpha, conf_r0, conf_power, num_occ, num_power,& - & num_alphas, xcnr, tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega,& - & camAlpha, camBeta, grid_params) + & max_alpha, num_alpha, tAutoAlphas, alpha, conf, num_occ, num_power, num_alphas, xcnr,& + & tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega, camAlpha, camBeta,& + & grid_params) !> nuclear charge, i.e. atomic number integer, intent(out) :: nuc @@ -53,11 +55,8 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min !> basis exponents real(dp), intent(out) :: alpha(0:4, 10) - !> confinement radii - real(dp), intent(out) :: conf_r0(0:4) - - !> power of confinement - real(dp), intent(out) :: conf_power(0:4) + !> confinement potential + type(TConf), intent(out) :: conf !> maximal occupied shell integer, intent(out) :: num_occ @@ -98,6 +97,9 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min !> holds parameters, defining a Becke integration grid type(TBeckeGridParams), intent(out) :: grid_params + !! confinement type of last query + integer :: last_conf_type + !! auxiliary variables integer :: ii, jj @@ -153,12 +155,38 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min stop end if - write(*, '(A)') 'Enter Confinement: r_0 and integer power, power=0.0 -> off' + last_conf_type = 0 do ii = 0, max_l - write(*, '(A,I3)') 'l=', ii - read(*,*) conf_r0(ii), conf_power(ii) + write(*, '(A,I3)') 'Enter confinement ID (0: none, 1: power, 2: WS) for l=', ii + read(*,*) conf%type + if (ii > 0) then + if (conf%type /= last_conf_type) then + error stop 'At the moment different shells are supposed to be compressed with the same& + & type of potential' + end if + end if + last_conf_type = conf%type end do + select case (conf%type) + case(confType%none) + continue + case(confType%power) + write(*, '(A)') 'Enter parameters r_0 and power' + do ii = 0, max_l + write(*, '(A,I3)') 'l=', ii + read(*,*) conf%r0(ii), conf%power(ii) + end do + case(confType%ws) + write(*, '(A)') 'Enter parameters Ronset, Rcut and Vmax' + do ii = 0, max_l + write(*, '(A,I3)') 'l=', ii + read(*,*) conf%onset(ii), conf%cutoff(ii), conf%vmax(ii) + end do + case default + error stop 'Invalid confinement potential.' + end select + write(*, '(A)') 'Enter number of occupied shells for each angular momentum.' do ii = 0, max_l write(*, '(A,I3)') 'l=', ii @@ -188,7 +216,8 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min end do else do ii = 0, max_l - write(*, '(A,I3,A,I3,A)') 'Enter ', num_alpha(ii), 'exponents for l=', ii, ', one per line.' + write(*, '(A,I3,A,I3,A)') 'Enter ', num_alpha(ii),& + & ' exponents for l=', ii, ', one per line.' do jj = 1, num_alpha(ii) read(*,*) alpha(ii, jj) end do @@ -265,8 +294,7 @@ end subroutine read_input_2 !> Echos gathered input to stdout. subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_alpha, alpha,& - & conf_r0, conf_power, occ, num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points,& - & xalpha_const) + & conf, occ, num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points, xalpha_const) !> nuclear charge, i.e. atomic number integer, intent(in) :: nuc @@ -292,11 +320,8 @@ subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_a !> basis exponents real(dp), intent(in) :: alpha(0:,:) - !> confinement radii - real(dp), intent(in) :: conf_r0(0:) - - !> power of confinement - real(dp), intent(in) :: conf_power(0:) + !> confinement potential + type(TConf), intent(in) :: conf !> occupation numbers real(dp), intent(in) :: occ(:,0:,:) @@ -391,12 +416,16 @@ subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_a write(*, '(A)') ' ' do ii = 0, max_l - if (conf_power(ii) > 1.0e-06_dp) then - write(*, '(A,I3,A,E15.7,A,E15.7)') 'l= ', ii, ', r0= ', conf_r0(ii), ' power= ',& - & conf_power(ii) - else - write(*, '(A,I3,A)') 'l= ', ii, ' no confinement ' - end if + select case (conf%type) + case(confType%none) + write(*, '(A,I3,A)') 'l= ', ii, ' no confinement' + case(confType%power) + write(*, '(A,I3,A,E15.7,A,E15.7)') 'l= ', ii, ', r0= ', conf%r0(ii),& + & ' power= ', conf%power(ii) + case(confType%ws) + write(*, '(A,I3,A,E15.7,A,E15.7,A,E15.7)') 'l= ', ii, ', Ronset= ', conf%onset(ii),& + & ' Rcutoff= ', conf%cutoff(ii), ' Vmax= ', conf%vmax(ii) + end select end do write(*, '(A)') ' ' diff --git a/slateratom/prog/main.f90 b/slateratom/prog/main.f90 index a4143149..dd53f925 100644 --- a/slateratom/prog/main.f90 +++ b/slateratom/prog/main.f90 @@ -4,7 +4,8 @@ program HFAtom use common_message, only : error use integration, only : gauss_chebyshev_becke_mesh use input, only : read_input_1, read_input_2, echo_input - use core_overlap, only : overlap, nuclear, kinetic, confinement + use core_overlap, only : overlap, nuclear, kinetic + use confinement, only : confType, getConf_power use coulomb_hfex, only : coulomb, hfex, hfex_lr use densitymatrix, only : densmatrix use hamiltonian, only : build_hamiltonian @@ -54,9 +55,8 @@ program HFAtom call parse_command_arguments() call read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min_alpha, max_alpha,& - & num_alpha, tAutoAlphas, alpha, conf_r0, conf_power, num_occ, num_power, num_alphas, xcnr,& - & tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega, camAlpha, camBeta,& - & grid_params) + & num_alpha, tAutoAlphas, alpha, conf, num_occ, num_power, num_alphas, xcnr, tPrintEigvecs,& + & tZora, mixnr, mixing_factor, xalpha_const, omega, camAlpha, camBeta, grid_params) problemsize = num_power * num_alphas @@ -73,8 +73,8 @@ program HFAtom if (nuc > 36) num_mesh_points = 1250 if (nuc > 54) num_mesh_points = 1500 - call echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_alpha, alpha, conf_r0,& - & conf_power, occ, num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points, xalpha_const) + call echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_alpha, alpha, conf, occ,& + & num_occ, num_power, num_alphas, xcnr, tZora, num_mesh_points, xalpha_const) ! allocate global stuff and zero out call allocate_globals() @@ -90,7 +90,15 @@ program HFAtom call overlap(ss, max_l, num_alpha, alpha, poly_order) call nuclear(uu, max_l, num_alpha, alpha, poly_order) call kinetic(tt, max_l, num_alpha, alpha, poly_order) - call confinement(vconf, max_l, num_alpha, alpha, poly_order, conf_r0, conf_power) + + select case (conf%type) + case(confType%none) + vconf(:,:,:) = 0.0_dp + case(confType%power) + call getConf_power(vconf, max_l, num_alpha, alpha, poly_order, conf) + case(confType%ws) + error stop 'Woods-Saxon compression not yet supported.' + end select ! test for linear dependency call diagonalize_overlap(max_l, num_alpha, poly_order, ss)