diis Subroutine

private subroutine diis(n, nmat, ndiis, x, e, b, xnew)

Uses

  • proc~~diis~~UsesGraph proc~diis diis module~mod_memory mod_memory proc~diis->module~mod_memory module~mod_io mod_io module~mod_memory->module~mod_io module~mod_constants mod_constants module~mod_memory->module~mod_constants iso_c_binding iso_c_binding module~mod_memory->iso_c_binding module~mod_io->module~mod_constants module~mod_constants->iso_c_binding

perform Pulay's direct inversion in the iterative subspace extrapolation:

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: n
integer(kind=ip), intent(inout) :: nmat
integer(kind=ip), intent(in) :: ndiis
real(kind=rp), intent(inout), dimension(n, ndiis) :: x
real(kind=rp), intent(inout), dimension(n, ndiis) :: e
real(kind=rp), intent(inout), dimension(ndiis+1, ndiis+1) :: b
real(kind=rp), intent(inout), dimension(n) :: xnew

Calls

proc~~diis~~CallsGraph proc~diis diis interface~mallocate mallocate proc~diis->interface~mallocate dgesv dgesv proc~diis->dgesv proc~makeb makeb proc~diis->proc~makeb interface~mfree mfree proc~diis->interface~mfree proc~r_alloc1 r_alloc1 interface~mallocate->proc~r_alloc1 proc~r_alloc3 r_alloc3 interface~mallocate->proc~r_alloc3 proc~i_alloc2 i_alloc2 interface~mallocate->proc~i_alloc2 proc~i_alloc1 i_alloc1 interface~mallocate->proc~i_alloc1 proc~r_alloc2 r_alloc2 interface~mallocate->proc~r_alloc2 proc~i_alloc3 i_alloc3 interface~mallocate->proc~i_alloc3 proc~l_alloc1 l_alloc1 interface~mallocate->proc~l_alloc1 proc~l_alloc2 l_alloc2 interface~mallocate->proc~l_alloc2 proc~i_free1 i_free1 interface~mfree->proc~i_free1 proc~l_free1 l_free1 interface~mfree->proc~l_free1 proc~r_free3 r_free3 interface~mfree->proc~r_free3 proc~r_free1 r_free1 interface~mfree->proc~r_free1 proc~i_free2 i_free2 interface~mfree->proc~i_free2 proc~r_free2 r_free2 interface~mfree->proc~r_free2 proc~l_free2 l_free2 interface~mfree->proc~l_free2 proc~i_free3 i_free3 interface~mfree->proc~i_free3 proc~chk_free chk_free proc~i_free1->proc~chk_free proc~l_free1->proc~chk_free proc~r_free3->proc~chk_free proc~chk_alloc chk_alloc proc~r_alloc1->proc~chk_alloc proc~memory_init memory_init proc~r_alloc1->proc~memory_init proc~r_alloc3->proc~chk_alloc proc~r_alloc3->proc~memory_init proc~i_alloc2->proc~chk_alloc proc~i_alloc2->proc~memory_init proc~i_alloc1->proc~chk_alloc proc~i_alloc1->proc~memory_init proc~r_alloc2->proc~chk_alloc proc~r_alloc2->proc~memory_init proc~i_alloc3->proc~chk_alloc proc~i_alloc3->proc~memory_init proc~l_alloc1->proc~chk_alloc proc~l_alloc1->proc~memory_init proc~l_alloc2->proc~chk_alloc proc~l_alloc2->proc~memory_init proc~r_free1->proc~chk_free proc~i_free2->proc~chk_free proc~r_free2->proc~chk_free proc~l_free2->proc~chk_free proc~i_free3->proc~chk_free proc~fatal_error fatal_error proc~chk_free->proc~fatal_error proc~chk_alloc->proc~fatal_error proc~ommp_message ommp_message proc~fatal_error->proc~ommp_message proc~close_output close_output proc~fatal_error->proc~close_output proc~close_output->proc~ommp_message

Called by

proc~~diis~~CalledByGraph proc~diis diis proc~jacobi_diis_solver jacobi_diis_solver proc~jacobi_diis_solver->proc~diis proc~polarization polarization proc~polarization->proc~jacobi_diis_solver proc~polelec_geomgrad polelec_geomgrad proc~polelec_geomgrad->proc~polarization proc~ommp_get_polelec_energy ommp_get_polelec_energy proc~ommp_get_polelec_energy->proc~polarization proc~ommp_set_external_field ommp_set_external_field proc~ommp_set_external_field->proc~polarization proc~ommp_polelec_geomgrad ommp_polelec_geomgrad proc~ommp_polelec_geomgrad->proc~polelec_geomgrad proc~c_ommp_get_polelec_energy C_ommp_get_polelec_energy proc~c_ommp_get_polelec_energy->proc~ommp_get_polelec_energy proc~ommp_get_full_ele_energy ommp_get_full_ele_energy proc~ommp_get_full_ele_energy->proc~ommp_get_polelec_energy proc~ommp_set_external_field_nomm ommp_set_external_field_nomm proc~ommp_set_external_field_nomm->proc~ommp_set_external_field proc~ommp_full_geomgrad ommp_full_geomgrad proc~ommp_full_geomgrad->proc~polelec_geomgrad proc~c_ommp_set_external_field_nomm C_ommp_set_external_field_nomm proc~c_ommp_set_external_field_nomm->proc~ommp_set_external_field proc~c_ommp_set_external_field C_ommp_set_external_field proc~c_ommp_set_external_field->proc~ommp_set_external_field proc~c_ommp_polelec_geomgrad C_ommp_polelec_geomgrad proc~c_ommp_polelec_geomgrad->proc~ommp_polelec_geomgrad proc~c_ommp_full_geomgrad C_ommp_full_geomgrad proc~c_ommp_full_geomgrad->proc~ommp_full_geomgrad proc~c_ommp_get_full_ele_energy C_ommp_get_full_ele_energy proc~c_ommp_get_full_ele_energy->proc~ommp_get_full_ele_energy proc~ommp_get_full_energy ommp_get_full_energy proc~ommp_get_full_energy->proc~ommp_get_full_ele_energy proc~c_ommp_get_full_energy C_ommp_get_full_energy proc~c_ommp_get_full_energy->proc~ommp_get_full_energy

Contents

Source Code

diis

Source Code

    subroutine diis(n,nmat,ndiis,x,e,b,xnew)
        !! perform Pulay's direct inversion in the iterative subspace extrapolation:
        use mod_memory, only: mallocate, mfree

        implicit none
        ! TODO doc
        integer(ip), intent(in) :: n, ndiis
        integer(ip), intent(inout) :: nmat
        real(rp), dimension(n, ndiis), intent(inout) :: x, e
        real(rp), dimension(ndiis+1, ndiis+1), intent(inout) :: b
        real(rp), dimension(n), intent(inout) :: xnew
        
        integer(ip) :: nmat1, i, info
        integer(ip) :: j, k

        real(rp),    allocatable :: bloc(:,:), cex(:)
        integer(ip), allocatable :: ipiv(:)

        if (nmat.ge.ndiis) then
            do j = 2, nmat - 10
                do k = 2, nmat - 10
                    b(j,k) = b(j+10,k+10)
                end do
            end do
      
            do j = 1, nmat - 10
                x(:,j) = x(:,j+10)
                e(:,j) = e(:,j+10)
            end do
            
            nmat = nmat - 10
        end if

        nmat1 = nmat + 1
    
        call mallocate('diis [bloc]', nmat1, nmat1, bloc)
        call mallocate('diis [cex]', nmat1, cex)
        call mallocate('diis [ipiv]', nmat1, ipiv)
    
        call makeb(n, nmat, ndiis, e, b)
        bloc = b(1:nmat1,1:nmat1)
        cex = 0.0_rp
        cex(1) = 1.0_rp

        call dgesv(nmat1, 1, bloc, nmat1, ipiv, cex, nmat1, info)
        
        if(info /= 0) then
            ! inversion failed. discard the previous points and restart.
            nmat = 1
            call mfree('diis [bloc]', bloc)
            call mfree('diis [cex]', cex)
            call mfree('diis [ipiv]', ipiv)
            return
        end if

        xnew = 0.0_rp
        do i = 1, nmat
            xnew = xnew + cex(i+1)*x(:,i)
        end do
        nmat = nmat + 1
        
        call mfree('diis [bloc]', bloc)
        call mfree('diis [cex]', cex)
        call mfree('diis [ipiv]', ipiv)
  
    end subroutine diis