From 4a9709acda057d102b175431a43e159c77b3a9b1 Mon Sep 17 00:00:00 2001 From: Tammo van der Heide Date: Sat, 22 Jun 2024 08:56:19 +0200 Subject: [PATCH] Mixer modernisation (#87) * Replace mixer implementations with more modern variants * Implement FYPP macros and re-enable assertions --- CMakeLists.txt | 4 +- common/include/common.fypp | 19 +- common/lib/CMakeLists.txt | 13 +- common/lib/assert.F90 | 41 ++ sktools/src/sktools/calculators/slateratom.py | 2 +- slateratom/lib/CMakeLists.txt | 22 +- slateratom/lib/blas.f90 | 48 ++ slateratom/lib/blasroutines.F90 | 51 ++ slateratom/lib/broyden.f90 | 532 ------------------ slateratom/lib/broydenmixer.F90 | 320 +++++++++++ slateratom/lib/globals.f90 | 54 +- slateratom/lib/hamiltonian.f90 | 28 +- slateratom/lib/input.f90 | 11 +- slateratom/lib/lapack.f90 | 38 ++ slateratom/lib/lapackroutines.F90 | 165 ++++++ slateratom/lib/mixer.f90 | 150 +++++ slateratom/lib/simplemixer.F90 | 74 +++ slateratom/prog/main.f90 | 12 +- test/prog/sktable/bin/testwithworkdir.py | 4 +- 19 files changed, 1016 insertions(+), 572 deletions(-) create mode 100644 common/lib/assert.F90 create mode 100644 slateratom/lib/blas.f90 create mode 100644 slateratom/lib/blasroutines.F90 delete mode 100644 slateratom/lib/broyden.f90 create mode 100644 slateratom/lib/broydenmixer.F90 create mode 100644 slateratom/lib/lapack.f90 create mode 100644 slateratom/lib/lapackroutines.F90 create mode 100644 slateratom/lib/mixer.f90 create mode 100644 slateratom/lib/simplemixer.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c43e8eb..8d8e9452 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -61,11 +61,13 @@ set(PYTHON_VERSION_MAJOR_MINOR "${Python3_VERSION_MAJOR}.${Python3_VERSION_MINOR set(FYPP "${PROJECT_SOURCE_DIR}/external/fypp/bin/fypp" CACHE FILEPATH "Fypp preprocessor") skprogs_add_fypp_defines(FYPP_FLAGS) +set(FYPP_FLAGS "${FYPP_FLAGS}") set(FYPP_CONFIG_FLAGS "${FYPP_FLAGS}") # Make sure, the line-marker option is not set list(REMOVE_ITEM FYPP_CONFIG_FLAGS "-n") -set(FYPP_BUILD_FLAGS "${FYPP_FLAGS}" "$,-DDEBUG=1,-DDEBUG=0>") +set(FYPP_BUILD_FLAGS "${FYPP_FLAGS}" "--file-var-root=${CMAKE_SOURCE_DIR}" + "$,-DDEBUG=1,-DDEBUG=0>") set(PYTHON_INTERPRETER "python3" CACHE STRING "Python interpreter to use for installing and test python components") diff --git a/common/include/common.fypp b/common/include/common.fypp index 7478c688..17ad7d6a 100644 --- a/common/include/common.fypp +++ b/common/include/common.fypp @@ -18,9 +18,26 @@ #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -#! DEBUG related macros +#! ASSERT and DEBUG related macros #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +#! Check a condition if WITH_ASSERT is True and call assertError if condition is False. +#! If an optional text string is included, print this in addition as an error +#:def ASSERT(cond, msg=None) + #:if WITH_ASSERT + if (.not. (${cond}$)) then + block + use common_assert, only : assertError + #:if msg + call assertError("${_FILE_}$", ${_LINE_}$, ${msg}$) + #:else + call assertError("${_FILE_}$", ${_LINE_}$) + #:endif + end block + end if + #:endif +#:enddef ASSERT + #! Insert code if DEBUG level is greater than zero. #:def DEBUG_CODE(code) #:if DEBUG > 0 diff --git a/common/lib/CMakeLists.txt b/common/lib/CMakeLists.txt index 931c9f05..7a3a653a 100644 --- a/common/lib/CMakeLists.txt +++ b/common/lib/CMakeLists.txt @@ -1,10 +1,10 @@ -file(TO_NATIVE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/../../" projectdir) +set(projectdir ${PROJECT_SOURCE_DIR}) # # General options for all targets # -set(fypp_flags ${FYPP_BUILD_FLAGS}) -list(APPEND fypp_flags -I${CMAKE_CURRENT_SOURCE_DIR}/../include -DRELEASE="'${RELEASE}'") +set(fypp_flags ${FYPP_BUILD_FLAGS} ${FYPP_CONFIG_FLAGS}) +list(APPEND fypp_flags -I${projectdir}/common/include -DRELEASE="'${RELEASE}'") set(sources-f90 accuracy.F90 @@ -30,7 +30,12 @@ set(sources-f90 taggedout.F90 utils.F90) -add_library(skprogs-common ${sources-f90}) +set(sources-fpp + assert.F90) + +skprogs_preprocess("${FYPP}" "${fypp_flags}" "F90" "f90" "${sources-fpp}" sources-f90-preproc) + +add_library(skprogs-common ${sources-f90} ${sources-f90-preproc}) set(moddir ${CMAKE_CURRENT_BINARY_DIR}/modfiles) set_target_properties(skprogs-common PROPERTIES Fortran_MODULE_DIRECTORY ${moddir}) diff --git a/common/lib/assert.F90 b/common/lib/assert.F90 new file mode 100644 index 00000000..abb65c12 --- /dev/null +++ b/common/lib/assert.F90 @@ -0,0 +1,41 @@ +#:include 'common.fypp' + +!> Auxiliary subroutines for the ASSERT command +module common_assert + + implicit none + + private +#:block DEBUG_CODE + public :: assertError +#:endblock DEBUG_CODE + +contains + +#:block DEBUG_CODE + + !> Prints assertion error and abort program execution. + subroutine assertError(fileName, lineNr, message) + + !> Name of the file in which the error occurred. + character(*), intent(in) :: fileName + + !> Nr. of the line at which the error occurred. + integer, intent(in) :: lineNr + + !> Additional message for error + character(*), intent(in), optional :: message + + write(*, '(A)') "!!! UNFULLFILLED ASSERTION" + write(*, '(A,A)') "!!! FILE: ", fileName + write(*, '(A,I0)') "!!! LINE NR.: ", lineNr + if (present(message)) then + write(*, '(A,A,A)') '!!! MESSAGE: "', trim(message), '"' + end if + error stop + + end subroutine assertError + +#:endblock DEBUG_CODE + +end module common_assert diff --git a/sktools/src/sktools/calculators/slateratom.py b/sktools/src/sktools/calculators/slateratom.py index c8e18530..5659ddb6 100644 --- a/sktools/src/sktools/calculators/slateratom.py +++ b/sktools/src/sktools/calculators/slateratom.py @@ -327,7 +327,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( - self._LOGICALSTRS[True], 0.1, self._COMMENT)) + 2, 0.1, self._COMMENT)) # Occupations for ll, occperl in enumerate(self._atomconfig.occupations): diff --git a/slateratom/lib/CMakeLists.txt b/slateratom/lib/CMakeLists.txt index df11f24f..78399778 100644 --- a/slateratom/lib/CMakeLists.txt +++ b/slateratom/lib/CMakeLists.txt @@ -1,6 +1,14 @@ +set(projectdir ${PROJECT_SOURCE_DIR}) + +# +# General options for all targets +# +set(fypp_flags ${FYPP_BUILD_FLAGS} ${FYPP_CONFIG_FLAGS}) +list(APPEND fypp_flags -I${projectdir}/common/include -DRELEASE="'${RELEASE}'") + set(sources-f90 avgpot.f90 - broyden.f90 + blas.f90 core_overlap.f90 coulomb_hfex.f90 coulomb_potential.f90 @@ -12,6 +20,8 @@ set(sources-f90 hamiltonian.f90 input.f90 integration.f90 + lapack.f90 + mixer.f90 numerical_differentiation.f90 output.f90 total_energy.f90 @@ -19,7 +29,15 @@ set(sources-f90 xcfunctionals.f90 zora_routines.f90) -add_library(skprogs-slateratom ${sources-f90}) +set(sources-fpp + blasroutines.F90 + broydenmixer.F90 + lapackroutines.F90 + simplemixer.F90) + +skprogs_preprocess("${FYPP}" "${fypp_flags}" "F90" "f90" "${sources-fpp}" sources-f90-preproc) + +add_library(skprogs-slateratom ${sources-f90} ${sources-f90-preproc}) target_link_libraries(skprogs-slateratom skprogs-common Libxc::xcf03 Libxc::xc) target_link_libraries(skprogs-slateratom skprogs-common LAPACK::LAPACK) diff --git a/slateratom/lib/blas.f90 b/slateratom/lib/blas.f90 new file mode 100644 index 00000000..ce701ec3 --- /dev/null +++ b/slateratom/lib/blas.f90 @@ -0,0 +1,48 @@ +!> Interface wrapper for the blas routines. +!! +!! ALL BLAS routines which are called from the main code must be included here. +module blas + + use common_accuracy, only : rdp + public + + interface + + + !> Performs the rank 1 operation + !! A := alpha*x*y**T + A, + subroutine dger(mm, nn, alpha, xx, incx, yy, incy, aa, lda) + import rdp + + !> Matrix sizing + integer, intent(in) :: mm + + !> Matrix size + integer, intent(in) :: nn + + !> Scale factor + real(rdp), intent(in) :: alpha + + !> Vector + real(rdp), intent(in) :: xx(*) + + !> Stride + integer, intent(in) :: incx + + !> Vector + real(rdp), intent(in) :: yy(*) + + !> Stride + integer, intent(in) :: incy + + !> Leading matrix dimension + integer, intent(in) :: lda + + !> Matrix A + real(rdp), intent(inout) :: aa(lda, *) + + end subroutine dger + + end interface + +end module blas diff --git a/slateratom/lib/blasroutines.F90 b/slateratom/lib/blasroutines.F90 new file mode 100644 index 00000000..55f24aaf --- /dev/null +++ b/slateratom/lib/blasroutines.F90 @@ -0,0 +1,51 @@ +#:include 'common.fypp' + +!> Contains F90 wrapper functions for some commonly used blas calls needed in the code. +!! The interface of all BLAS calls must be defined in the module blas. +module blasroutines + + use common_accuracy, only : rdp + use blas, only : dger + implicit none + + private + public :: ger + + + !> Rank 1 update of a matrix A := alpha*x*y' + A + !! Wrapper for the level 2 blas routine xger to perform the rank 1 update of a general matrix + interface ger + module procedure ger_dble + end interface ger + + +contains + + !> Double precision rank 1 update of a general matrix + subroutine ger_dble(a, alpha, x, y) + + !> Contains the matrix for the update + real(rdp), intent(inout) :: a(:,:) + + !> Scaling value for the update contribution + real(rdp), intent(in) :: alpha + + !> Vector of values for the update + real(rdp), intent(in) :: x(:) + + !> Vector of values for the update + real(rdp), intent(in) :: y(:) + + integer :: n, m + + @:ASSERT(size(a,dim=1) == size(x)) + @:ASSERT(size(a,dim=2) == size(y)) + + m = size(x) + n = size(y) + + call dger(m, n, alpha, x, 1, y, 1, a, m) + + end subroutine ger_dble + +end module blasroutines diff --git a/slateratom/lib/broyden.f90 b/slateratom/lib/broyden.f90 deleted file mode 100644 index 0ae18072..00000000 --- a/slateratom/lib/broyden.f90 +++ /dev/null @@ -1,532 +0,0 @@ -!> Module that provides the density mixing functionality (simple, Broyden). -module broyden - - use common_accuracy, only : dp - - implicit none - private - - public :: mixing_driver - - -contains - - !> This is the main driver for simple and broyden mixers, both mix one big one-dimensional array. - subroutine mixing_driver(pot_old, pot_new, max_l, num_alpha, poly_order, problemsize, iter,& - & tBroyden, mixing_factor) - - !> old potential - real(dp), intent(in) :: pot_old(:,0:,:,:) - - !> contains current potential on entry and mixed potential on exit - real(dp), intent(inout) :: pot_new(:,0:,:,:) - - !> maximum angular momentum - integer, intent(in) :: max_l - - !> number of exponents in each shell - integer, intent(in) :: num_alpha(0:) - - !> highest polynomial order + l in each shell - integer, intent(in) :: poly_order(0:) - - !> maximum size of the eigenproblem - integer, intent(in) :: problemsize - - !> current SCF iteration - integer, intent(in) :: iter - - !> true, if Broyden mixing is desired, otherwise simple mixing is applied - logical, intent(in) :: tBroyden - - !> mixing factor - real(dp), intent(in) :: mixing_factor - - !> serialized potentials - real(dp), allocatable :: vecin(:), vecout(:) - - !> equals current SCF iteration or next one if (iter == 0) - integer :: titer - - !> auxiliary variables - integer :: ii, jj, kk, ll, pp - - allocate(vecout(size(pot_old))) - allocate(vecin(size(pot_old))) - - ! serialize potentials - pp = 0 - do ii = 1, 2 - do jj = 0, max_l - do kk = 1, num_alpha(jj) * poly_order(jj) - do ll = 1, problemsize - pp = pp + 1 - vecin(pp) = pot_old(ii, jj, kk, ll) - vecout(pp) = pot_new(ii, jj, kk, ll) - end do - end do - end do - end do - - ! this check is still necessary, since max. 10000 entries is hardcoded in broyden_mixer - if (pp > 10000) then - write(*,*) 'Static dimensions in broyden_mixer too small: ', pp - stop - end if - - titer = iter - ! broyden returns if (iter == 0) - if (iter == 0) titer = 1 - - if (tBroyden) then - call broyden_mixer(titer, mixing_factor, size(vecin), vecin, vecout) - else - call simple_mix(vecin, vecout, mixing_factor) - end if - - ! deserialize obtained potential - pp = 0 - do ii = 1, 2 - do jj = 0, max_l - do kk = 1, num_alpha(jj) * poly_order(jj) - do ll = 1, problemsize - pp = pp + 1 - pot_new(ii, jj, kk, ll) = vecin(pp) - end do - end do - end do - end do - - end subroutine mixing_driver - - - !> Simple mixer for last and current density. - subroutine simple_mix(last, cur, factor) - - !> old vector, holds mixed values at exit - real(dp), intent(inout) :: last(:) - - !> new vector - real(dp), intent(in) :: cur(:) - - !> mixing factor - real(dp), intent(in) :: factor - - last(:) = factor * cur + (1.0_dp - factor) * last - - end subroutine simple_mix - - -! -!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - SUBROUTINE BROYDEN_mixer(NITER,ALPHA,JTOP,VECIN,VECOUT) - -! This is the Broyden routine as also implemented in the old DFTB code. - - IMPLICIT real(dp) (A-H,O-Z) - IMPLICIT INTEGER (I-N) -! -!************************************************************ -!* THE VECTORS UI(MAXSIZ) AND VTI(MAXSIZ) ARE JOHNSON'S * -!* U(OF I ) AND DF(TRANSPOSE), RESPECTIVELY. THESE ARE * -!* CONTINUALLY UPDATED. ALL ITERATIONS ARE S7ORED ON TAPE * -!* 32 . THIS IS DONE TO PREVENT THE PROHIBITIVE STORAGE * -!* COSTS ASSOCIATED WITH HOLDING ONTO THE ENTIRE JACOBIAN. * -!* VECTOR TL IS THE VT OF EARLIER ITERATIONS. VECTOR F IS: * -!* VECTOR(OUTPUT) - VECTOR(IN). VECTOR DF IS: F(M+1)-F(M) * -!* FINALLY,VECTOR DUMVI(MAXSIZ) IS THE PREDICTED VECTOR. * -!* ON RETURN, VECIN CONTAINS THE NEW TRIAL VECTOR. * -!************************************************************ -!* FOR THE CRAY2-CIVIC ENVIRONMENT , FILES 32 AND 31 * -!* SHOULD BE INTRODUCED IN THE LINK STATEMENT. * -!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - PARAMETER (ZERO=0.0D0,ONE=1.0D0,IMATSZ=40,maxsiz=10000) -! formerly IMATSZ=90 -! -! ADDED PARAMETER MAXITER. POREZAG, MAY 1995 -! - PARAMETER(MAXITER=15) -! -! replaced writing to disk by storing values in -! arrays UNIT31, UNIT32 hajnal@scientist.com 2000-10-4 -!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! -! CHARACTER*7 NAMES -! -! SCRATCH COMMON BLOCK FOR LOCAL VARIABLES -! - DIMENSION VECIN(*),VECOUT(*) - DIMENSION F(MAXSIZ),UI(MAXSIZ),VTI(MAXSIZ),T1(MAXSIZ),& - & VECTOR(MAXSIZ,2),DUMVI(MAXSIZ),DF(MAXSIZ) -! DIMENSION NAMES(3) - DIMENSION A(IMATSZ,IMATSZ),B(IMATSZ,IMATSZ),CM(IMATSZ) - DIMENSION D(IMATSZ,IMATSZ),W(IMATSZ) - DIMENSION UNIT31(MAXSIZ,2),UNIT32(MAXSIZ,2,MAXITER+15) -! DATA NAMES/'BROYD01','BROYD02','BROYD03'/ - real(dp) UAMIX,WTMP - INTEGER ILASTIT - common /broyd/ uamix, w, WTMP, unit31, unit32, ilastit - save -! -! PRINT *,'IN MIXING, WHERE ARE YOU?' -! -! NEW LINES JULY 1996 -! - IF (JTOP .GT. MAXSIZ) THEN - PRINT *,'BROYDEN: JTOP > MAXSIZ' - STOP - END IF -! -! NEW LINES POREZAG, MAY 1995 -! - ITER=NITER - IF(NITER.GT.MAXITER)ITER=MOD(ITER,MAXITER)+1 - IF(ITER.EQ.0)RETURN -! -! END NEW LINES -! -! OPEN(66,FILE=NAMES(1),STATUS='UNKNOWN',FORM='FORMATTED') -! REWIND(66) -! OPEN(31,FILE=NAMES(2),STATUS='UNKNOWN',FORM='UNFORMATTED') -! OPEN(32,FILE=NAMES(3),STATUS='UNKNOWN',FORM='UNFORMATTED') -! REWIND(31) -! REWIND(32) - -! IF(ITER.EQ.1)THEN -! ENDFILE 31 -! ENDFILE 32 -! END IF -! -! -!++++++ SET UP THE VECTOR OF THE CURRENT ITERATION FOR MIXING ++++++ -! -! FOR THIS METHOD WE HAVE ONLY SAVED INPUT/OUTPUT CHG. DENSITIES, - DO K=1,JTOP - VECTOR(K,1)= VECIN(K) - VECTOR(K,2)= VECOUT(K) - END DO -!++++++ END OF PROGRAM SPECIFIC LOADING OF VECTOR FROM MAIN ++++++++ -! -! IVSIZ IS THE LENGTH OF THE VECTOR - IVSIZ=JTOP -! IF(ITER.LT.3)WRITE( 6,1001)IVSIZ - IF(IVSIZ.GT.MAXSIZ)THEN - PRINT *,'MIXING: EXCEEDED MAXIMAL VECTOR LENGTH' - STOP - END IF -! -! -!******************* BEGIN BROYDEN'S METHOD ********************** -! -! WEIGHTING FACTOR FOR THE ZEROTH ITERATION - W0=0.01D0 -! -! F: THE DIFFERENCE OF PREVIOUS OUTPUT AND INPUT VECTORS -! DUMVI: A DUMMY VECTOR, HERE IT IS THE PREVIOUS INPUT VECTOR -! REWIND(31) -! READ(31,END=119,ERR=119)AMIX,LASTIT - IF (ITER .EQ. 1) THEN - GOTO 119 - ELSE - AMIX=UAMIX - LASTIT=ILASTIT - END IF -! READ(31)(F(K),K=1,IVSIZ) - DO k=1,IVSIZ - F(k)=UNIT31(K,1) - END DO -! READ(31)(DUMVI(K),K=1,IVSIZ) - DO k=1,IVSIZ - DUMVI(k)=UNIT31(K,2) - END DO -! IF(ITER.EQ.1 .AND. LASTIT.GT.1)THEN -! READ(31)LTMP,((A(I,J),I=1,LTMP),J=1,LTMP) -! READ(31)(W(I),I=1,LTMP) -! ENDIF -! -! ALPHA(OR AMIX)IS SIMPLE MIXING PARAMETERS -! WRITE(66,1002)AMIX,ITER+1 -! - DO K=1,IVSIZ - DUMVI(K)=VECTOR(K,1)-DUMVI(K) - DF(K)=VECTOR(K,2)-VECTOR(K,1)-F(K) - END DO - DO K=1,IVSIZ - F(K)=VECTOR(K,2)-VECTOR(K,1) - END DO -! -! FOR I-TH ITER.,DFNORM IS ( F(I) MINUS F(I-1) ), USED FOR NORMALIZATION -! - DFNORM=ZERO - FNORM=ZERO - DO K=1,IVSIZ - DFNORM=DFNORM + DF(K)*DF(K) - FNORM=FNORM + F(K)*F(K) - END DO - DFNORM=SQRT(DFNORM) - FNORM=SQRT(FNORM) -! WRITE(66,'('' DFNORM '',E12.6,'' FNORM '',E12.6)')DFNORM,FNORM -! - FAC2=ONE/DFNORM - FAC1=AMIX*FAC2 -! - DO K=1,IVSIZ - UI(K) = FAC1*DF(K) + FAC2*DUMVI(K) - VTI(K)= FAC2*DF(K) - END DO -! -!*********** CALCULATION OF COEFFICIENT MATRICES ************* -!*********** AND THE SUM FOR CORRECTIONS ************* -! -! RECALL: A(I,J) IS A SYMMETRIC MATRIX -! : B(I,J) IS THE INVERSE OF [ W0**2 I + A ] -! - LASTIT=LASTIT+1 - LASTM1=LASTIT-1 - LASTM2=LASTIT-2 -! -! DUMVI IS THE U(OF I) AND T1 IS THE VT(OF I) -! FROM THE PREVIOUS ITERATIONS -! REWIND(32) -! WRITE(66,1003)LASTIT,LASTM1 - IF(LASTIT.GT.2)THEN - DO J=1,LASTM2 -! READ(32)(DUMVI(K),K=1,IVSIZ) - DO k=1,IVSIZ - DUMVI(k)=UNIT32(k,1,J) - END DO -! READ(32)(T1(K),K=1,IVSIZ) - DO k=1,IVSIZ - T1(k)=UNIT32(k,2,J) - END DO -! - AIJ=ZERO - CMJ=ZERO - DO K=1,IVSIZ - CMJ=CMJ + T1(K)*F(K) - AIJ=AIJ + T1(K)*VTI(K) - END DO - A(LASTM1,J)=AIJ - A(J,LASTM1)=AIJ - CM(J)=CMJ - END DO - ENDIF -! - AIJ=ZERO - CMJ=ZERO - DO K=1,IVSIZ - CMJ= CMJ + VTI(K)*F(K) - AIJ= AIJ + VTI(K)*VTI(K) - END DO - A(LASTM1,LASTM1)=AIJ - CM(LASTM1)=CMJ -! -! WRITE(32)(UI(K),K=1,IVSIZ) - DO k=1,IVSIZ - UNIT32(k,1,LASTM1)=UI(k) - END DO -! WRITE(32)(VTI(K),K=1,IVSIZ) - DO k=1,IVSIZ - UNIT32(k,2,LASTM1)=VTI(k) - END DO -! REWIND(32) -! -! THE WEIGHTING FACTORS FOR EACH ITERATION HAVE BEEN CHOSEN -! EQUAL TO ONE OVER THE R.M.S. ERROR. THIS NEED NOT BE THE CASE. - IF(FNORM .GT. 1.0D-7)THEN - WTMP=0.010D0/FNORM - ELSE - WTMP=1.0D5 - END IF - IF(WTMP.LT. 1.00D0) then - WTMP=1.00D0 - end if -! print *,wtmp,lastm1,w(lastm1) - W(LASTM1)=WTMP -! WRITE(66,'('' WEIGHTING SET = '',E12.6)')WTMP -! -! -! WITH THE CURRENT ITERATIONS F AND VECTOR CALCULATED, -! WRITE THEM TO UNIT 31 FOR USE LATER. -! REWIND(31) -! WRITE(31)AMIX,LASTIT - UAMIX=AMIX - ILASTIT=LASTIT -! WRITE(31)(F(K),K=1,IVSIZ) - DO k=1,IVSIZ - UNIT31(K,1)=F(k) - END DO -! WRITE(31)(VECTOR(K,1),K=1,IVSIZ) - DO k=1,IVSIZ - UNIT31(K,2)=VECTOR(K,1) - END DO -! WRITE(31)LASTM1,((A(I,J),I=1,LASTM1),J=1,LASTM1) -! WRITE(31)(W(I),I=1,LASTM1) -! -! SET UP AND CALCULATE BETA MATRIX - DO LM=1,LASTM1 - DO LN=1,LASTM1 - D(LN,LM)= A(LN,LM)*W(LN)*W(LM) - B(LN,LM)= ZERO - END DO - B(LM,LM)= ONE - D(LM,LM)= W0**2 + A(LM,LM)*W(LM)*W(LM) - END DO -! - CALL INVERSE(D,B,LASTM1) -! -! CALCULATE THE VECTOR FOR THE NEW ITERATION - DO K=1,IVSIZ - DUMVI(K)= VECTOR(K,1) + AMIX*F(K) - END DO -! - DO I=1,LASTM1 -! READ(32)(UI(K),K=1,IVSIZ) - DO k=1,IVSIZ - UI(k)=UNIT32(k,1,I) - END DO -! READ(32)(VTI(K),K=1,IVSIZ) - DO k=1,IVSIZ - VTI(k)=UNIT32(k,2,I) - END DO - GMI=ZERO - DO IP=1,LASTM1 - GMI=GMI + CM(IP)*B(IP,I)*W(IP) - END DO - DO K=1,IVSIZ - DUMVI(K)=DUMVI(K)-GMI*UI(K)*W(I) - END DO - END DO -! END OF THE CALCULATION OF DUMVI, THE NEW VECTOR -! -! REWIND(31) -! REWIND(32) -! - GOTO 120 -! IF THIS IS THE FIRST ITERATION, THEN LOAD -! F=VECTOR(OUT)-VECTOR(IN) AND VECTOR(IN) - 119 CONTINUE -! PRINT*,'SIMPLE MIXING THIS ITERATION' -! REWIND(31) - LASTIT=1 - AMIX=ALPHA -! WRITE(31)AMIX,LASTIT - UAMIX=AMIX - ILASTIT=LASTIT - DO K=1,IVSIZ - F(K)=VECTOR(K,2)-VECTOR(K,1) - END DO -! WRITE(31)(F(K),K=1,IVSIZ) - DO k=1,IVSIZ - UNIT31(K,1)=F(k) - END DO -! WRITE(31)(VECTOR(K,1),K=1,IVSIZ) - DO k=1,IVSIZ - UNIT31(K,2)=VECTOR(K,1) - END DO -! -! SINCE WE ARE ON THE FIRST ITERATION, SIMPLE MIX THE VECTOR. - DO K=1,IVSIZ - DUMVI(K)= VECTOR(K,1) + AMIX*F(K) - END DO -! WRITE( 6,1000) - 120 CONTINUE -! -! CLOSE(31,STATUS='KEEP') -! CLOSE(32,STATUS='KEEP') -! -!************* THE END OF THE BROYDEN METHOD ************** -! -!+++++++ PROGRAM SPECIFIC CODE OF RELOADING ARRAYS +++++++++ -! -! NEED TO UNLOAD THE NEW VECTOR INTO THE APPROPRIATE ARRAYS. - DO K=1,JTOP - VECIN(K)=DUMVI(K) - END DO -! -!+++++++++ END OF PROGRAM SPECIFIC RELOADING OF ARRAYS +++++++++ -! -! WRITE(66,1004)ITER+1 -! CLOSE(66) - RETURN -! - 1000 FORMAT(' ----> STRAIGHT MIXING ON THIS ITERATION') - 1001 FORMAT(' IN MIXING: IVSIZ =',I7,/) - 1002 FORMAT(' IN MIXING: SIMPLE MIX PARAMETER',1(F10.6,',')& - & ,' FOR ITER=',I5) - 1003 FORMAT(' CURRENT ITER= ',I5,' INCLUDES VALUES FROM ITER=',I5) - 1004 FORMAT(10X,'DENSITY FOR ITERATION',I4,' PREPARED') - END subroutine broyden_mixer -! -! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - SUBROUTINE INVERSE(A,B,M) - IMPLICIT real(dp) (A-H,O-Z) - IMPLICIT INTEGER (I-N) - -! ============================================================= -! -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - PARAMETER (IMATSZ=40) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -! - DIMENSION A(IMATSZ,IMATSZ),B(IMATSZ,IMATSZ) - DIMENSION TD(IMATSZ),AD(IMATSZ),BD(IMATSZ) - SAVE -! -! SUBROUTINE TO PREFORM GAUSSIAN ELIMINATION -! NO ZEROS ALONG THE DIAGONAL -! - N=M - IF(N.GT.IMATSZ)THEN - PRINT *,'INVERT: MATRIX A TOO LARGE' - STOP - END IF -! - DO I=1,N - ATMP=A(I,I) - IF(ABS(ATMP) .LT. 1.0D-08)THEN -! WRITE(66,'('' INVERT: MATRIX HAS ZERO DIAGONAL'', -! & '' ELEMENT IN THE '',I4,'' ROW'')')I - STOP - ENDIF - END DO -! - IF(N.EQ.1) GO TO 605 -! - DO I=1,N -! - DO J=1,N - TD(J)=A(J,I)/A(I,I) - END DO -! -! TD(I)=(0.0E+00,0.0E+00) - TD(I)=0.0D0 -! - DO K=1,N - BD(K)=B(I,K) - AD(K)=A(I,K) - END DO -! - DO K=1,N - DO J=1,N - B(J,K)=B(J,K)-(TD(J)*BD(K)) - A(J,K)=A(J,K)-(TD(J)*AD(K)) - END DO - END DO -! - END DO -! - DO I=1,N - DO J=1,N - B(J,I)=B(J,I)/A(J,J) - END DO - END DO -! - RETURN -! - 605 B(1,1)=1.0D0/A(1,1) - RETURN - END subroutine inverse -! - -end module broyden diff --git a/slateratom/lib/broydenmixer.F90 b/slateratom/lib/broydenmixer.F90 new file mode 100644 index 00000000..b29a7e03 --- /dev/null +++ b/slateratom/lib/broydenmixer.F90 @@ -0,0 +1,320 @@ +#:include 'common.fypp' + +!> Contains a modified Broyden mixer. +!! The modified Broyden mixer implemented here is practically the same as the one in the old DFTB +!! code. A detailed description of the method can be found in Johnson's paper. +!! See: D.D. Johnson, PRB 38, 12807 (1988) +!! DOI: 10.1103/PhysRevB.38.12807 +!! In order to use the mixer you have to create and reset it. +module broydenmixer + + use common_accuracy, only : dp + use blasroutines, only : ger + use lapackroutines, only : getrf, getrs + implicit none + + private + public :: TBroydenMixer, TBroydenMixer_init, TBroydenMixer_reset, TBroydenMixer_mix + + + !> Contains the necessary data for a Broyden mixer. + type TBroydenMixer + private + + !> Actual iteration + integer :: iIter + + !> Nr. of maximal iterations + integer :: mIter + + !> Nr. of elements in the vectors + integer :: nElem + + !> Jacobi matrix differences + real(dp) :: omega0 + + !> Mixing parameter + real(dp) :: alpha + + !> Minimal weight + real(dp) :: minWeight + + !> Maximal weight + real(dp) :: maxWeight + + !> Weighting factor (numerator) + real(dp) :: weightFac + + !> Weights for prev. iterations + real(dp), allocatable :: ww(:) + + !> Charge difference in last iteration + real(dp), allocatable :: qDiffLast(:) + + !> Input charge in last iteration + real(dp), allocatable :: qInpLast(:) + + !> Storage for the "a" matrix + real(dp), allocatable :: aa(:,:) + + !> DF vectors + real(dp), allocatable :: dF(:,:) + + !> uu vectors + real(dp), allocatable :: uu(:,:) + + end type TBroydenMixer + + +contains + + !> Creates a Broyden mixer instance. + !! The weight associated with an iteration is calculated as weigthFac/ww where ww is the Euclidean + !! norm of the charge difference vector. If the calculated weigth is outside of the + !! [minWeight, maxWeight] region it is replaced with the appropriate boundary value. + subroutine TBroydenMixer_init(this, mIter, mixParam, omega0, minWeight, maxWeight, weightFac) + + !> An initialized Broyden mixer on exit + type(TBroydenMixer), intent(out) :: this + + !> Maximum nr. of iterations (max. nr. of vectors to store) + integer, intent(in) :: mIter + + !> Mixing parameter + real(dp), intent(in) :: mixParam + + !> Weight for the Jacobi matrix differences + real(dp), intent(in) :: omega0 + + !> Minimal weight allowed + real(dp), intent(in) :: minWeight + + !> Maximal weight allowed + real(dp), intent(in) :: maxWeight + + !> Numerator of the weight + real(dp), intent(in) :: weightFac + + @:ASSERT(mIter > 0) + @:ASSERT(mixParam > 0.0_dp) + @:ASSERT(omega0 > 0.0_dp) + + this%nElem = 0 + this%mIter = mIter + this%alpha = mixParam + this%omega0 = omega0 + this%minWeight = minWeight + this%maxWeight = maxWeight + this%weightFac = weightFac + allocate(this%ww(mIter-1)) + allocate(this%qInpLast(this%nElem)) + allocate(this%qDiffLast(this%nElem)) + allocate(this%aa(mIter-1, mIter-1)) + allocate(this%dF(this%nElem, mIter - 1)) + allocate(this%uu(this%nElem, mIter - 1)) + + end subroutine TBroydenMixer_init + + + !> Makes the mixer ready for a new SCC cycle. + subroutine TBroydenMixer_reset(this, nElem) + + !> Broyden mixer instance + type(TBroydenMixer), intent(inout) :: this + + !> Length of the vectors to mix + integer, intent(in) :: nElem + + @:ASSERT(nElem > 0) + + if (nElem /= this%nElem) then + this%nElem = nElem + deallocate(this%qInpLast) + deallocate(this%qDiffLast) + allocate(this%qInpLast(this%nElem)) + allocate(this%qDiffLast(this%nElem)) + deallocate(this%dF) + allocate(this%dF(this%nElem, this%mIter - 1)) + deallocate(this%uu) + allocate(this%uu(this%nElem, this%mIter - 1)) + end if + this%iIter = 0 + this%ww(:) = 0.0_dp + this%aa(:,:) = 0.0_dp + + end subroutine TBroydenMixer_reset + + + !> Mixes charges according to the modified Broyden method. + !! + !! Warning: The complex-valued Broyden mixer requires flattened hermitian matrices as input. + !! You are free to permute the individual elements of the flattened arrays as long as the same + !! permutation is applied to qInpResult and qDiff. + !! The restriction arises from the assumption that the dot-products of density matrices are + !! real-valued (imaginary parts add up to zero due to the hermitian property) and the linear + !! system of equations remains real-valued. + subroutine TBroydenMixer_mix(this, qInpResult, qDiff) + + !> The Broyden mixer + type(TBroydenMixer), intent(inout) :: this + + !> Input charges on entry, mixed charges on exit + real(dp), intent(inout) :: qInpResult(:) + + !> Charge difference between output and input charges + real(dp), intent(in) :: qDiff(:) + + this%iIter = this%iIter + 1 + if (this%iIter > this%mIter) then + error stop "Broyden mixer: Maximal nr. of steps exceeded" + end if + + call modifiedBroydenMixing(qInpResult, this%qInpLast, this%qDiffLast, this%aa,& + & this%ww, this%iIter, qDiff, this%alpha, this%omega0, this%minWeight, this%maxWeight,& + & this%weightFac, this%nElem, this%dF, this%uu) + + end subroutine TBroydenMixer_mix + + + !> Does the real work for the Broyden mixer. + subroutine modifiedBroydenMixing(qInpResult, qInpLast, qDiffLast, aa, ww, nn, qDiff,& + & alpha, omega0, minWeight, maxWeight, weightFac, nElem, dF, uu) + + !> Current input charge on entry, mixed charged on exit + real(dp), intent(inout) :: qInpResult(:) + + !> Input charge vector of the previous iterations + real(dp), intent(inout) :: qInpLast(:) + + !> Charge difference of the previous iteration + real(dp), intent(inout) :: qDiffLast(:) + + !> The matrix a (needed for the mixing) + real(dp), intent(inout) :: aa(:,:) + + !> Weighting factors of the iterations + real(dp), intent(inout) :: ww(:) + + !> Current iteration number + integer, intent(in) :: nn + + !> Charge difference of the current iteration + real(dp), intent(in) :: qDiff(:) + + !> Mixing parameter + real(dp), intent(in) :: alpha + + !> Weight for the Jacobi matrix differences + real(dp), intent(in) :: omega0 + + !> Minimal weight allowed + real(dp), intent(in) :: minWeight + + !> Maximal weight allowed + real(dp), intent(in) :: maxWeight + + !> Numerator of the weight + real(dp), intent(in) :: weightFac + + !> Nr. of elements in the vectors + integer, intent(in) :: nElem + + !> Prev. DFs + real(dp), intent(inout) :: dF(:,:) + + !> Prev. U vectors + real(dp), intent(inout) :: uu(:,:) + + real(dp), allocatable :: beta(:,:), cc(:) + + ! Current DF or U-vector + real(dp), allocatable :: dF_uu(:) + + real(dp) :: invNorm + integer :: ii, nn_1 + integer, allocatable :: ipiv(:) + + nn_1 = nn - 1 + + @:ASSERT(nn > 0) + @:ASSERT(size(qInpResult) == nElem) + @:ASSERT(size(qInpLast) == nElem) + @:ASSERT(size(qDiffLast) == nElem) + @:ASSERT(size(qDiff) == nElem) + @:ASSERT(all(shape(aa) >= [nn_1, nn_1])) + @:ASSERT(size(ww) >= nn_1) + + ! First iteration: simple mix and storage of qInp and qDiff + if (nn == 1) then + qInpLast(:) = qInpResult + qDiffLast(:) = qDiff + qInpResult(:) = qInpResult + alpha * qDiff + return + end if + + allocate(beta(nn_1, nn_1)) + allocate(cc(nn_1)) + allocate(dF_uu(nElem)) + allocate(ipiv(nn_1)) + + ! Create weight factor omega for current iteration + ww(nn_1) = sqrt(dot_product(qDiff, qDiff)) + if (ww(nn_1) > weightFac / maxWeight) then + ww(nn_1) = weightFac / ww(nn_1) + else + ww(nn_1) = maxWeight + end if + if (ww(nn_1) < minWeight) then + ww(nn_1) = minWeight + end if + + ! Build |DF(m-1)> and (m is the current iteration number) + dF_uu(:) = qDiff - qDiffLast + invNorm = sqrt(dot_product(dF_uu, dF_uu)) + invNorm = max(invNorm, epsilon(1.0_dp)) + invNorm = 1.0_dp / invNorm + dF_uu(:) = invNorm * dF_uu + + ! Build a, beta, c, and gamma + ! (due to the hermitian property of our density matrices, the dot-products below are real) + do ii = 1, nn - 2 + aa(ii, nn_1) = dot_product(dF(:,ii), dF_uu) + aa(nn_1, ii) = aa(ii, nn_1) + cc(ii) = ww(ii) * dot_product(dF(:,ii), qDiff) + end do + aa(nn_1, nn_1) = 1.0_dp + cc(nn_1) = ww(nn_1) * dot_product(dF_uu, qDiff) + + do ii = 1, nn_1 + beta(:nn_1, ii) = ww(:nn_1) * ww(ii) * aa(:nn_1,ii) + beta(ii, ii) = beta(ii, ii) + omega0**2 + end do + + ! LU decomposition + call getrf(beta, ipiv) + ! Solve system of linear equations by using the LU decomposition + call getrs(beta, ipiv, cc, trans='t') + + ! Store |dF(m-1)> + dF(:, nn_1) = dF_uu + + ! Create |u(m-1)> + dF_uu(:) = alpha * dF_uu + invNorm * (qInpResult - qInpLast) + + ! Save charge vectors before overwriting + qInpLast(:) = qInpResult + qDiffLast(:) = qDiff + + ! Build new vector + qInpResult(:) = qInpResult + alpha * qDiff + do ii = 1, nn-2 + qInpResult(:) = qInpResult - ww(ii) * cc(ii) * uu(:,ii) + end do + qInpResult(:) = qInpResult - ww(nn_1) * cc(nn_1) * dF_uu + + ! Save |u(m-1)> + uu(:, nn_1) = dF_uu + + end subroutine modifiedBroydenMixing + +end module broydenmixer diff --git a/slateratom/lib/globals.f90 b/slateratom/lib/globals.f90 index d7ea616d..6687874d 100644 --- a/slateratom/lib/globals.f90 +++ b/slateratom/lib/globals.f90 @@ -2,6 +2,9 @@ module globals use common_accuracy, only : dp + use mixer, only : TMixer, TMixer_init, TMixer_reset, mixerTypes + use broydenmixer, only : TBroydenMixer, TBroydenMixer_init + use simplemixer, only : TSimpleMixer, TSimpleMixer_init implicit none @@ -167,15 +170,24 @@ module globals !> true, if SCF cycle reached convergency logical :: tConverged - !> true, if Broyden mixing is desired, otherwise simple mixing is applied - logical :: tBroyden + !> identifier of mixer + integer :: mixnr - !> true, if average local, effective potential should be calculated - logical :: isAvgPotNeeded + !> mixer instance + type(TMixer), allocatable :: pMixer + + !> simple mixer (if used) + type(TSimpleMixer), allocatable :: pSimpleMixer + + !> broyden mixer (if used) + type(TBroydenMixer), allocatable :: pBroydenMixer !> mixing factor real(dp) :: mixing_factor + !> true, if average local, effective potential should be calculated + logical :: isAvgPotNeeded + !> zora kinetic energy contribution to total energy real(dp) :: zora_ekin @@ -188,6 +200,9 @@ module globals !> Allocates all the variables in the globals module. subroutine allocate_globals() + !! auxiliary variables to count the elements to mix + integer :: ind, idx1, idx2, idx3, idx4 + allocate(weight(num_mesh_points)) allocate(abcissa(num_mesh_points)) allocate(dzdr(num_mesh_points)) @@ -239,6 +254,37 @@ subroutine allocate_globals() if (isAvgPotNeeded) allocate(avgPot(num_mesh_points, 2), source=0.0_dp) + ! initialize mixer + allocate(pMixer) + select case(mixnr) + case(mixerTypes%simple) + allocate(pSimplemixer) + call TSimpleMixer_init(pSimpleMixer, mixing_factor) + call TMixer_init(pMixer, pSimpleMixer) + case(mixerTypes%broyden) + allocate(pBroydenMixer) + ! defaults taken from DFTB+ + call TBroydenMixer_init(pBroydenMixer, maxiter, mixing_factor, 0.01_dp, 1.0_dp, 1.0e5_dp,& + & 1.0e-2_dp) + call TMixer_init(pMixer, pBroydenMixer) + case default + error stop "Unknown mixer type." + end select + + ! count elements to mix + ind = 0 + do idx1 = 1, 2 + do idx2 = 0, max_l + do idx3 = 1, num_alpha(idx2) * poly_order(idx2) + do idx4 = 1, problemsize + ind = ind + 1 + end do + end do + end do + end do + + call TMixer_reset(pMixer, ind) + end subroutine allocate_globals end module globals diff --git a/slateratom/lib/hamiltonian.f90 b/slateratom/lib/hamiltonian.f90 index fc088c8c..dab8d521 100644 --- a/slateratom/lib/hamiltonian.f90 +++ b/slateratom/lib/hamiltonian.f90 @@ -3,7 +3,7 @@ module hamiltonian use common_accuracy, only : dp use dft, only : dft_exc_matrixelement - use broyden, only : mixing_driver + use mixer, only : TMixer, TMixer_mix use zora_routines, only : zora_t_correction use xcfunctionals, only : xcFunctional @@ -17,9 +17,12 @@ module hamiltonian contains !> Main driver routine for Fock matrix build-up. Also calls mixer with potential matrix. - subroutine build_hamiltonian(iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha,& - & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha, pot_old,& - & pot_new, tZora, tBroyden, mixing_factor, ff, camAlpha, camBeta) + subroutine build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l,& + & num_alpha, poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha,& + & pot_old, pot_new, tZora, ff, camAlpha, camBeta) + + !> mixer instances + type(TMixer), intent(inout) :: pMixer !> current SCF step integer, intent(in) :: iScf @@ -87,12 +90,6 @@ subroutine build_hamiltonian(iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, !> true, if zero-order regular approximation for relativistic effects is desired logical, intent(in) :: tZora - !> true, if Broyden mixing is desired, otherwise simple mixing is applied - logical, intent(in) :: tBroyden - - !> mixing factor - real(dp), intent(in) :: mixing_factor - !> fock matrix supervector real(dp), intent(out) :: ff(:,0:,:,:) @@ -108,8 +105,11 @@ subroutine build_hamiltonian(iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, !> auxiliary matrices for (CAM) hybrids real(dp), allocatable :: k_matrix2(:,:,:,:), k_matrix3(:,:,:,:) + !> potential difference + real(dp), allocatable :: pot_diff(:,:,:,:) + !! auxiliary variables - integer :: ii, jjj, kkk, ll, mm, ss, ttt, iMix + integer :: ii, jjj, kkk, ll, mm, ss, ttt ff(:,:,:,:) = 0.0_dp @@ -180,9 +180,9 @@ subroutine build_hamiltonian(iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, pot_new(2, :,:,:) = - real(nuc, dp) * uu + j_matrix - k_matrix(2, :,:,:) ! mixer - iMix = int(iScf / 40) - call mixing_driver(pot_old, pot_new, max_l, num_alpha, poly_order, problemsize,& - & iScf - iMix * 40, tBroyden, mixing_factor) + allocate(pot_diff, mold=pot_old) + pot_diff(:,0:,:,:) = pot_old - pot_new + call TMixer_mix(pMixer, pot_new, pot_diff) ! Not sure: before or after mixer (potential .ne. Matrix elements)? ! Should be irrelevant once self-consistency is reached. diff --git a/slateratom/lib/input.f90 b/slateratom/lib/input.f90 index a5ddc130..51104403 100644 --- a/slateratom/lib/input.f90 +++ b/slateratom/lib/input.f90 @@ -17,7 +17,7 @@ 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, tBroyden, mixing_factor, xalpha_const, omega,& + & num_alphas, xcnr, tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega,& & camAlpha, camBeta, grid_params) !> nuclear charge, i.e. atomic number @@ -77,8 +77,8 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min !> true, if zero-order regular approximation for relativistic effects is desired logical, intent(out) :: tZora - !> true, if Broyden mixing is desired, otherwise simple mixing is applied - logical, intent(out) :: tBroyden + !> identifier of mixer + integer, intent(out) :: mixnr !> mixing factor real(dp), intent(out) :: mixing_factor @@ -213,8 +213,9 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min write(*, '(A)') 'Print Eigenvectors ? .true./.false.' read(*,*) tPrintEigvecs - write(*, '(A)') ' Use Broyden mixer (.true./.false.)? And mixing parameter <1' - read(*,*) tBroyden, mixing_factor + write(*, '(A)') 'Enter mixer and mixing parameter <1:& + & 1: Simple mixer, 2: Broyden mixer' + read(*,*) mixnr, mixing_factor end subroutine read_input_1 diff --git a/slateratom/lib/lapack.f90 b/slateratom/lib/lapack.f90 new file mode 100644 index 00000000..dd4f9db5 --- /dev/null +++ b/slateratom/lib/lapack.f90 @@ -0,0 +1,38 @@ +!> Interface wrapper for the lapack routines. See the lapack +!! project documentation for more details +module lapack + + use common_accuracy, only : rsp, rdp + implicit none + public + + + !> Computes LU factorization of double precision matrix + interface dgetrf + + !> Computes LU factorization of double precision matrix + subroutine dgetrf(mm, nn, aa, lda, ipiv, info) + import rdp + + !> number of rows of the matrix + integer, intent(in) :: mm + + !> matrix dimension + integer, intent(in) :: nn + + !> Leading dimension of A + integer, intent(in) :: lda + + !> matrix A + real(rdp), intent(inout) :: aa(lda, *) + + !> pivot array + integer, intent(out) :: ipiv(*) + + !> state of routine on return + integer, intent(out) :: info + end subroutine dgetrf + + end interface dgetrf + +end module lapack diff --git a/slateratom/lib/lapackroutines.F90 b/slateratom/lib/lapackroutines.F90 new file mode 100644 index 00000000..390ba7cc --- /dev/null +++ b/slateratom/lib/lapackroutines.F90 @@ -0,0 +1,165 @@ +#:include 'common.fypp' + +!> Contains F90 wrapper functions for some commonly used lapack calls needed in the code. +!! The interface of all LAPACK calls must be defined in the module lapack. +module lapackroutines + + use common_accuracy, only : dp, rdp + implicit none + + private + public :: getrf, getrs + + + !> Computes the LU decomposition of a general rectangular matrix using partial pivoting with row + !! interchanges. + !! The decomposition has the form: A = P*L*U, where P is a permutation matrix, L is a lower + !! triangular matrix with unit diagonal elements and U is an upper triangular matrix. + interface getrf + module procedure getrf_dble + end interface getrf + + + !> Solves a system of linear equations + !! A * X = B or A**T * X = B + !! with a general N-by-N matrix A using the LU factorization computed by getrf. + interface getrs + module procedure :: getrs_dble + module procedure :: getrs1_dble + end interface getrs + + +contains + + !> Double precision version of getrf. + subroutine getrf_dble(aa, ipiv, nRow, nColumn, iError) + + !> Matrix to decompose on entry, L and U on exit. Unit diagonal elements of L are not stored. + real(rdp), intent(inout) :: aa(:,:) + + !> Pivot indices, row i of the matrix was interchanged with row ipiv(i). + integer, intent(out) :: ipiv(:) + + !> Number of rows of the matrix to decomposea. (Necessary if different from the number of rows + !> of the passed matrix) + integer, optional, intent(in) :: nRow + + !> Number of rows of the matrix to decompose. (Necessary if different from the number of columns + !> of the passed matrix) + integer, optional, intent(in) :: nColumn + + !> Error flag. Zero on successful exit. If not present, any lapack error causes program + !> termination. If passed only fatal lapack errors with error flag < 0 cause abort. + integer, optional, intent(out) :: iError + + integer :: mm, nn, lda, info + character(len=100) :: error_string + + lda = size(aa, dim=1) + nn = size(aa, dim=2) + if (present(nRow)) then + @:ASSERT(nRow >= 1 .and. nRow <= lda) + mm = nRow + else + mm = lda + end if + if (present(nColumn)) then + @:ASSERT(nColumn >= 1 .and. nColumn <= nn) + nn = nColumn + end if + @:ASSERT(size(ipiv) == min(mm, nn)) + + call dgetrf(mm, nn, aa, lda, ipiv, info) + + if (info < 0) then +99060 format('Failure in LU factorisation dgetrf,', ' illegal argument at position ', i10) + write(error_string, 99060) info + error stop error_string + else + if (present(iError)) then + iError = info + elseif (info > 0) then +99070 format('Factor U is exactly zero in dgetrf,', ' info flag is ', i10) + write(error_string, 99070) info + error stop error_string + end if + end if + + end subroutine getrf_dble + + + !> Solves a system of linear equations with multiple right hand sides + subroutine getrs_dble(amat, ipiv, bmat, trans, iError) + + !> Matrix of the linear system + real(rdp), intent(in) :: amat(:, :) + + !> Pivot indices, row i of the matrix was interchanged with row ipiv(i). + integer, intent(in) :: ipiv(:) + + !> Matrix of the right hand side vectors + real(rdp), intent(inout) :: bmat(:, :) + + !> Optional transpose (defaults to 'n') + character(len=1), intent(in), optional :: trans + + !> Error flag, zero on successful exit + integer, intent(out), optional :: iError + + character(len=1) :: atr + integer :: info, nn, nrhs, lda, ldb + + @:ASSERT(size(amat, 1) == size(amat, dim=2)) + @:ASSERT(size(amat, 1) == size(bmat, dim=1)) + + if (present(trans)) then + @:ASSERT(any(trans == ['n', 'N', 't', 'T', 'c', 'C'])) + atr = trans + else + atr = 'n' + end if + + lda = max(1, size(amat, 1)) + ldb = max(1, size(bmat, 1)) + nn = size(amat, 2) + nrhs = size(bmat, 2) + + call dgetrs(atr, nn, nrhs, amat, lda, ipiv, bmat, ldb, info) + + if (present(iError)) then + iError = info + else + if (info /= 0) then + error stop "Failed to solve linear system by diagonal pivoting" + end if + end if + + end subroutine getrs_dble + + + !> Solves a system of linear equations with one right hand side + subroutine getrs1_dble(amat, ipiv, bvec, trans, iError) + + !> Matrix of the linear system + real(rdp), intent(in) :: amat(:,:) + + !> Pivot indices, row i of the matrix was interchanged with row ipiv(i). + integer, intent(in) :: ipiv(:) + + !> Right hand side vector + real(rdp), intent(inout), target :: bvec(:) + + !> optional transpose (defaults to 'n') + character(len=1), intent(in), optional :: trans + + !> Error flag, zero on successful exit + integer, intent(out), optional :: iError + + real(rdp), pointer :: bptr(:,:) + + bptr(1:size(bvec, 1), 1:1) => bvec(1:size(bvec, 1)) + call getrs(amat, ipiv, bptr, trans, iError) + + end subroutine getrs1_dble + +end module lapackroutines diff --git a/slateratom/lib/mixer.f90 b/slateratom/lib/mixer.f90 new file mode 100644 index 00000000..1ac8230f --- /dev/null +++ b/slateratom/lib/mixer.f90 @@ -0,0 +1,150 @@ +!> Provides a general mixer which contains the desired actual mixer. +module mixer + + use common_accuracy, only : dp + use broydenmixer, only : TBroydenMixer, TBroydenMixer_mix, TBroydenMixer_reset + use simplemixer, only : TSimpleMixer, TSimpleMixer_mix, TSimpleMixer_reset + implicit none + + private + public :: TMixer, TMixer_init, TMixer_reset, TMixer_mix, mixerTypes + + + !> Interface type for mixers + type TMixer + private + + !> Numerical type of mixer 1:2 + integer :: mixerType + + !> Simple mixer instance + type(TSimpleMixer), allocatable :: pSimpleMixer + + !> Broyden mixer instance + type(TBroydenMixer), allocatable :: pBroydenMixer + + end type TMixer + + + !> Initialises a specific mixer + interface TMixer_init + module procedure TMixer_initSimple + module procedure TMixer_initBroyden + end interface TMixer_init + + + !> Mixes the given quantity + interface TMixer_mix + module procedure TMixer_mix1D + module procedure TMixer_mix4D + end interface TMixer_mix + + + type :: TMixerTypesEnum + integer :: simple = 1 + integer :: broyden = 2 + end type TMixerTypesEnum + + !> Contains mixer types + type(TMixerTypesEnum), parameter :: mixerTypes = TMixerTypesEnum() + + +contains + + !> Initializes a simple mixer. + subroutine TMixer_initSimple(this, pSimple) + + !> Mixer instance + type(TMixer), intent(out) :: this + + !> A valid Simple mixer instance on exit + type(TSimpleMixer), allocatable, intent(inout) :: pSimple + + this%mixerType = mixerTypes%simple + call move_alloc(pSimple, this%pSimpleMixer) + + end subroutine TMixer_initSimple + + + !> Initializes a Broyden mixer. + subroutine TMixer_initBroyden(this, pBroyden) + + !> Mixer instance + type(TMixer), intent(out) :: this + + !> A valid Broyden mixer instance on exit + type(TBroydenMixer), allocatable, intent(inout) :: pBroyden + + this%mixerType = mixerTypes%broyden + call move_alloc(pBroyden, this%pBroydenMixer) + + end subroutine TMixer_initBroyden + + + !> Resets the mixer. + subroutine TMixer_reset(this, nElem) + + !> Mixer instance + type(TMixer), intent(inout) :: this + + !> Size of the vectors to mix + integer, intent(in) :: nElem + + select case (this%mixerType) + case(mixerTypes%simple) + call TSimpleMixer_reset(this%pSimpleMixer, nElem) + case(mixerTypes%broyden) + call TBroydenMixer_reset(this%pBroydenMixer, nElem) + end select + + end subroutine TMixer_reset + + + !> Mixes two vectors. + subroutine TMixer_mix1D(this, inp, diff) + + !> Mixer instance + type(TMixer), intent(inout) :: this + + !> Input vector on entry, result vector on exit + real(dp), intent(inout) :: inp(:) + + !> Difference between input and output vectors (measure of lack of convergence) + real(dp), intent(in) :: diff(:) + + select case (this%mixerType) + case(mixerTypes%simple) + call TSimpleMixer_mix(this%pSimpleMixer, inp, diff) + case(mixerTypes%broyden) + call TBroydenMixer_mix(this%pBroydenMixer, inp, diff) + end select + + end subroutine TMixer_mix1D + + + !> Mixes two 4D matrices. + subroutine TMixer_mix4D(this, inp, diff) + + !> Mixer instance + type(TMixer), intent(inout) :: this + + !> Input vector on entry, result vector on exit + real(dp), intent(inout), contiguous, target :: inp(:,0:,:,:) + + !> Difference between input and output vectors (measure of lack of convergence) + real(dp), intent(in), contiguous, target :: diff(:,0:,:,:) + + !! Difference between input and output vectors (1D pointer) + real(dp), pointer :: pDiff(:) + + !! Input vector on entry, result vector on exit (1D pointer) + real(dp), pointer :: pInp(:) + + pInp(1:size(inp)) => inp + pDiff(1:size(diff)) => diff + + call TMixer_mix1D(this, pInp, pDiff) + + end subroutine TMixer_mix4D + +end module mixer diff --git a/slateratom/lib/simplemixer.F90 b/slateratom/lib/simplemixer.F90 new file mode 100644 index 00000000..fa354afd --- /dev/null +++ b/slateratom/lib/simplemixer.F90 @@ -0,0 +1,74 @@ +#:include 'common.fypp' + +!> Simple mixer for mixing charges. +module simplemixer + + use common_accuracy, only : dp + implicit none + + private + public :: TSimpleMixer + public :: TSimpleMixer_init, TSimpleMixer_reset, TSimpleMixer_mix + + + !> Contains data for a simple mixer + type TSimpleMixer + private + + !> Mixing parameter + real(dp) :: mixParam + + end type TSimpleMixer + + +contains + + !> Creates a simple mixer. + subroutine TSimpleMixer_init(this, mixParam) + + !> Simple mixer instance on exit + type(TSimpleMixer), intent(out) :: this + + !> Mixing parameter + real(dp), intent(in) :: mixParam + + this%mixParam = mixParam + + end subroutine TSimpleMixer_init + + + !> Resets the mixer. + subroutine TSimpleMixer_reset(this, nElem) + + !> Simple mixer instance + type(TSimpleMixer), intent(inout) :: this + + !> Length of the vectors to mix + integer, intent(in) :: nElem + + @:ASSERT(nElem > 0) + + continue + + end subroutine TSimpleMixer_reset + + + !> Does the actual mixing. + subroutine TSimpleMixer_mix(this, qInpResult, qDiff) + + !> SimpleMixer instance + type(TSimpleMixer), intent(inout) :: this + + !> Input charge on entry, mixed charge on exit + real(dp), intent(inout) :: qInpResult(:) + + !> Charge difference + real(dp), intent(in) :: qDiff(:) + + @:ASSERT(size(qInpResult) == size(qDiff)) + + qInpResult(:) = qInpResult + this%mixParam * qDiff + + end subroutine TSimpleMixer_mix + +end module simplemixer diff --git a/slateratom/prog/main.f90 b/slateratom/prog/main.f90 index 04b6fde9..a4143149 100644 --- a/slateratom/prog/main.f90 +++ b/slateratom/prog/main.f90 @@ -55,7 +55,7 @@ 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, tBroyden, mixing_factor, xalpha_const, omega, camAlpha, camBeta,& + & tPrintEigvecs, tZora, mixnr, mixing_factor, xalpha_const, omega, camAlpha, camBeta,& & grid_params) problemsize = num_power * num_alphas @@ -125,9 +125,9 @@ program HFAtom pot_old(:,:,:,:) = 0.0_dp ! kinetic energy, nuclear-electron, and confinement matrix elements which are constant during SCF - call build_hamiltonian(0, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha, poly_order,& - & problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha, pot_old, pot_new, tZora,& - & tBroyden, mixing_factor, ff, camAlpha, camBeta) + call build_hamiltonian(pMixer, 0, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha,& + & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha, pot_old,& + & pot_new, tZora, ff, camAlpha, camBeta) ! self-consistency cycles write(*,*) 'Energies in Hartree' @@ -149,9 +149,9 @@ program HFAtom & dz, xcnr, omega, camAlpha, camBeta, rho, drho, ddrho, vxc, exc, xalpha_const) ! build Fock matrix and get total energy during SCF - call build_hamiltonian(iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha,& + call build_hamiltonian(pMixer, iScf, tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha,& & poly_order, problemsize, xcnr, num_mesh_points, weight, abcissa, vxc, alpha, pot_old,& - & pot_new, tZora, tBroyden, mixing_factor, ff, camAlpha, camBeta) + & pot_new, tZora, ff, camAlpha, camBeta) if (tZora) then call getTotalEnergyZora(tt, uu, nuc, vconf, jj, kk, kk_lr, pp, max_l, num_alpha, poly_order,& diff --git a/test/prog/sktable/bin/testwithworkdir.py b/test/prog/sktable/bin/testwithworkdir.py index 083a99b9..582887f7 100644 --- a/test/prog/sktable/bin/testwithworkdir.py +++ b/test/prog/sktable/bin/testwithworkdir.py @@ -16,8 +16,8 @@ SKDEF = 'skdef.hsd' -ATOL = 1e-10 -RTOL = 1e-09 +ATOL = 1e-08 +RTOL = 1e-07 class TestWithWorkDir: