diff --git a/sktools/src/sktools/calculators/slateratom.py b/sktools/src/sktools/calculators/slateratom.py index f429d284..4fd25de3 100644 --- a/sktools/src/sktools/calculators/slateratom.py +++ b/sktools/src/sktools/calculators/slateratom.py @@ -341,8 +341,8 @@ def write(self, workdir): 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})" + compr.ww, compr.aa, compr.r0, self._COMMENT) \ + + " Compr. height, slope, half-height radius ({:s})" .format(sc.ANGMOM_TO_SHELL[ll])] out += ["{:d} \t\t\t{:s} nr. of occupied shells ({:s})".format( diff --git a/sktools/src/sktools/compressions.py b/sktools/src/sktools/compressions.py index ec54c48b..c2c69b4d 100644 --- a/sktools/src/sktools/compressions.py +++ b/sktools/src/sktools/compressions.py @@ -69,42 +69,40 @@ class WoodsSaxonCompression(sc.ClassDict): Attributes ---------- - onset : float - Onset radius of the compression - cutoff : float - Cutoff radius of the compression - vmax : float - Potential well depth/height + ww : float + Height (W) of the compression + aa : float + Slope (a) of the compression + r0 : float + Half-height radius of the compression ''' @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) + ww, child = query.getvalue( + root, 'w', conv.float0, defvalue=100.0, returnchild=True) + if ww < 0.0: + msg = 'Invalid potential height (W) {:f}'.format(ww) 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) + aa, child = query.getvalue(root, 'a', conv.float0, returnchild=True) + if aa <= 0.0: + msg = 'Invalid potential slope (a) {:f}'.format(aa) 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) + r0, child = query.getvalue(root, 'r0', conv.float0, returnchild=True) + if r0 < 0.0: + msg = 'Invalid potential half-height radius {:f}'.format(r0) 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 + myself.ww = ww + myself.aa = aa + myself.r0 = r0 return myself @@ -117,16 +115,16 @@ def tohsd(self, root, query, parentname=None): 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) + query.setchildvalue(mynode, 'w', conv.float0, self.ww) + query.setchildvalue(mynode, 'a', conv.float0, self.aa) + query.setchildvalue(mynode, 'r0', conv.float0, self.r0) 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 + ww_ok = abs(self.ww - other.ww) < 1e-8 + aa_ok = abs(self.aa - other.aa) < 1e-8 + r0_ok = abs(self.r0 - other.r0) < 1e-8 + return ww_ok and aa_ok and r0_ok # Registered compressions with corresponding hsd name as key diff --git a/slateratom/lib/CMakeLists.txt b/slateratom/lib/CMakeLists.txt index c23d9ad2..1c47c778 100644 --- a/slateratom/lib/CMakeLists.txt +++ b/slateratom/lib/CMakeLists.txt @@ -18,7 +18,6 @@ set(sources-f90 diagonalizations.f90 globals.f90 hamiltonian.f90 - input.f90 integration.f90 lapack.f90 mixer.f90 @@ -33,6 +32,7 @@ set(sources-fpp blasroutines.F90 broydenmixer.F90 confinement.F90 + input.F90 lapackroutines.F90 simplemixer.F90) diff --git a/slateratom/lib/confinement.F90 b/slateratom/lib/confinement.F90 index 777dcfaf..1ae96498 100644 --- a/slateratom/lib/confinement.F90 +++ b/slateratom/lib/confinement.F90 @@ -107,14 +107,14 @@ end subroutine getSupervec !> Input for the Woods-Saxon confinement. type :: TWsConfInp - !> onset radii (Woods-Saxon compression) - real(dp) :: rOnset(0:4) + !> potential heights (W) + real(dp) :: ww(0:4) - !> cutoff of confinement (Woods-Saxon compression) - real(dp) :: rCut(0:4) + !> potential slopes (a) + real(dp) :: aa(0:4) - !> potential well height of confinement (Woods-Saxon compression) - real(dp) :: vMax(0:4) + !> half-height radii (r0) + real(dp) :: r0(0:4) end type TWsConfInp @@ -135,10 +135,10 @@ end subroutine getSupervec type, extends(TConf) :: TPowerConf private - !> confinement radii (power compression) + !> confinement radii real(dp) :: r0(0:4) - !> power of confinement (power compression) + !> power of confinement real(dp) :: power(0:4) contains @@ -153,14 +153,14 @@ end subroutine getSupervec type, extends(TConf) :: TWsConf private - !> onset radii (Woods-Saxon compression) - real(dp) :: rOnset(0:4) + !> potential heights (W) + real(dp) :: ww(0:4) - !> cutoff of confinement (Woods-Saxon compression) - real(dp) :: rCut(0:4) + !> potential slopes (a) + real(dp) :: aa(0:4) - !> potential well height of confinement (Woods-Saxon compression) - real(dp) :: vMax(0:4) + !> half-height radii (r0) + real(dp) :: r0(0:4) contains @@ -214,9 +214,9 @@ subroutine TWsConf_init(this, input) !> input data type(TWsConfInp), intent(inout) :: input - this%rOnset(0:) = input%rOnset - this%rCut(0:) = input%rCut - this%vMax(0:) = input%vMax + this%ww(0:) = input%ww + this%aa(0:) = input%aa + this%r0(0:) = input%r0 end subroutine TWsConf_init @@ -239,16 +239,14 @@ subroutine TPowerConf_getPotOnGrid(this, max_l, num_mesh_points, abcissa, vconf) !> Woods-Saxon potential on grid real(dp), intent(out) :: vconf(:, 0:) - !! auxiliary variables - integer :: ii, ll + !! angular momentum + integer :: ll @:ASSERT(size(vconf, dim=1) == num_mesh_points) @:ASSERT(size(vconf, dim=2) == max_l + 1) do ll = 0, max_l - do ii = 1, num_mesh_points - vconf(ii, ll) = getPowerPotential(abcissa(ii), this%r0(ll), this%power(ll)) - end do + vconf(:, ll) = getPowerPotential(abcissa, this%r0(ll), this%power(ll)) end do end subroutine TPowerConf_getPotOnGrid @@ -272,16 +270,14 @@ subroutine TWsConf_getPotOnGrid(this, max_l, num_mesh_points, abcissa, vconf) !> Woods-Saxon potential on grid real(dp), intent(out) :: vconf(:, 0:) - !! auxiliary variables - integer :: ii, ll + !! angular momentum + integer :: ll @:ASSERT(size(vconf, dim=1) == num_mesh_points) @:ASSERT(size(vconf, dim=2) == max_l + 1) do ll = 0, max_l - do ii = 1, num_mesh_points - vconf(ii, ll) = getWSPotential(abcissa(ii), this%rOnset(ll), this%rCut(ll), this%vMax(ll)) - end do + vconf(:, ll) = getWSPotential(abcissa, this%ww(ll), this%aa(ll), this%r0(ll)) end do end subroutine TWsConf_getPotOnGrid @@ -499,7 +495,7 @@ end subroutine dft_conf_matrixelement !> Returns the Power potential for a given radial distance. - pure function getPowerPotential(rr, r0, power) result(pot) + elemental function getPowerPotential(rr, r0, power) result(pot) !> radial distance real(dp), intent(in) :: rr @@ -513,37 +509,35 @@ pure function getPowerPotential(rr, r0, power) result(pot) !> Power potential at evaluation point real(dp) :: pot + @:ASSERT(rr >= 0.0_dp) + pot = (rr / r0)**power end function getPowerPotential !> Returns the Woods-Saxon potential for a given radial distance. - pure function getWSPotential(rr, rOnset, rCut, vMax) result(pot) + !! see J. Chem. Theory Comput. 12, 1, 53-64 (2016) eqn. 4. + elemental function getWSPotential(rr, ww, aa, r0) result(pot) !> radial distance real(dp), intent(in) :: rr - !> onset radius - real(dp), intent(in) :: rOnset + !> potential height (W) + real(dp), intent(in) :: ww - !> cutoff radius - real(dp), intent(in) :: rCut + !> potential slope (a) + real(dp), intent(in) :: aa - !> potential well height - real(dp), intent(in) :: vMax + !> half-height radius (r0) + real(dp), intent(in) :: r0 !> Woods-Saxon potential at evaluation point real(dp) :: pot - ! case: rr <= rOnset - pot = 0.0_dp + @:ASSERT(rr >= 0.0_dp) - if (rr > rOnset .and. rr < rCut) then - pot = vMax / (rCut - rOnset) * (rr - rOnset) - elseif (rr >= rCut) then - pot = vMax - end if + pot = ww / (1.0 + exp(-aa * (rr - r0))) end function getWSPotential diff --git a/slateratom/lib/input.f90 b/slateratom/lib/input.F90 similarity index 97% rename from slateratom/lib/input.f90 rename to slateratom/lib/input.F90 index 16e9b994..138e3b94 100644 --- a/slateratom/lib/input.f90 +++ b/slateratom/lib/input.F90 @@ -1,3 +1,5 @@ +#:include 'common.fypp' + !> Module that reads input values from stdin. module input @@ -180,10 +182,13 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min read(*,*) confInp%power%r0(ii), confInp%power%power(ii) end do case(confType%ws) - write(*, '(A)') 'Enter parameters Ronset, Rcut and Vmax' + write(*, '(A)') 'Enter parameters compr. height, slope and half-height radius' do ii = 0, max_l write(*, '(A,I3)') 'l=', ii - read(*,*) confInp%ws%rOnset(ii), confInp%ws%rCut(ii), confInp%ws%vMax(ii) + read(*,*) confInp%ws%ww(ii), confInp%ws%aa(ii), confInp%ws%r0(ii) + @:ASSERT(confInp%ws%ww(ii) >= 0.0) + @:ASSERT(confInp%ws%aa(ii) > 0.0) + @:ASSERT(confInp%ws%r0(ii) >= 0.0) end do case default error stop 'Invalid confinement potential.' @@ -429,8 +434,8 @@ subroutine echo_input(nuc, max_l, occ_shells, maxiter, scftol, poly_order, num_a write(*, '(A,I3,A,E15.7,A,E15.7)') 'l= ', ii, ', r0= ', confInp%power%r0(ii),& & ' power= ', confInp%power%power(ii) case(confType%ws) - write(*, '(A,I3,A,E15.7,A,E15.7,A,E15.7)') 'l= ', ii, ', Ronset= ', confInp%ws%rOnset(ii),& - & ' Rcutoff= ', confInp%ws%rCut(ii), ' Vmax= ', confInp%ws%vMax(ii) + write(*, '(A,I3,A,E15.7,A,E15.7,A,E15.7)') 'l= ', ii, ', W= ', confInp%ws%ww(ii),& + & ' a= ', confInp%ws%aa(ii), ' r0= ', confInp%ws%r0(ii) end select end do