comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Refactoring of linres solver (#2665)

Browse source
This commit is contained in:
Juerg Hutter 2023-03-09 13:57:14 +01:00 • committed by GitHub
parent 07ce17837e
commit e19d499b64
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
10 changed files with 187 additions and 83 deletions
Download patch file Download diff file Expand all files Collapse all files

8
src/mp2_cphf.F
View file

@ -112,7 +112,8 @@ MODULE mp2_cphf
USE qs_ks_types, ONLY: qs_ks_env_type,&
set_ks_env
USE qs_linres_types, ONLY: linres_control_type
USE qs_mo_types, ONLY: mo_set_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
USE qs_overlap, ONLY: build_overlap_matrix
USE qs_p_env_methods, ONLY: p_env_check_i_alloc,&
@ -1316,7 +1317,7 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'update_mp2_forces'
INTEGER :: alpha, beta, handle, idir, iounit, &
ispin, nspins
ispin, nocc, nspins
INTEGER, DIMENSION(3) :: comp
LOGICAL :: do_exx, do_hfx, use_virial
REAL(KIND=dp) :: e_dummy, e_hartree, e_xc, ehartree, exc
@ -1948,8 +1949,9 @@ CONTAINS
! Finish matrix_w_mp2 with occ-occ block
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), homo=nocc, nmo=alpha)
CALL calculate_whz_matrix(mos(ispin)%mo_coeff, p_env%kpp1(ispin)%matrix, &
p_env%w1(1)%matrix, 1.0_dp)
p_env%w1(1)%matrix, 1.0_dp, nocc)
END DO
IF (debug_forces .AND. use_virial) e_dummy = third_tr(virial%pv_virial)

84
src/preconditioner.F
View file

@ -16,7 +16,13 @@ MODULE preconditioner
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
USE cp_fm_types, ONLY: cp_fm_get_info,&
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_get_info,&
cp_fm_release,&
cp_fm_to_fm,&
cp_fm_type
USE dbcsr_api, ONLY: dbcsr_p_type,&
dbcsr_type
@ -45,7 +51,8 @@ MODULE preconditioner
qs_environment_type
USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_type
mo_set_type,&
set_mo_set
#include "./base/base_uses.f90"
IMPLICIT NONE
@ -106,20 +113,22 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner'
INTEGER :: handle, k, my_solver_type
INTEGER :: handle, k, my_solver_type, nao, nhomo
LOGICAL :: my_convert_precond_to_dbcsr, &
needs_full_spectrum, needs_homo, &
use_mo_coeff_b
REAL(KIND=dp) :: energy_homo
REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_ot
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_ot
TYPE(cp_fm_struct_type), POINTER :: fm_struct
TYPE(cp_fm_type) :: mo_occ
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(dbcsr_type), POINTER :: mo_coeff_b
CALL timeset(routineN, handle)
CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b)
CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b, homo=nhomo)
use_mo_coeff_b = mo_set%use_mo_coeff_b
CALL cp_fm_get_info(mo_coeff, ncol_global=k)
CALL cp_fm_get_info(mo_coeff, ncol_global=k, nrow_global=nao)
! Starting some matrix mess, check where to store the result in preconditioner_env, fm or dbcsr_matrix
my_convert_precond_to_dbcsr = .FALSE.
@ -152,9 +161,9 @@ CONTAINS
CPABORT("The preconditioner is unknown ...")
END SELECT
ALLOCATE (eigenvalues_ot(k))
energy_homo = 0.0_dp
IF (needs_full_spectrum) THEN
ALLOCATE (eigenvalues_ot(k))
! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either
! the case of non-equivalent MOs but also for the derivate
! we could already have all eigenvalues e.g. full_all and we could skip this
@ -170,7 +179,10 @@ CONTAINS
CALL calculate_subspace_eigenvalues(mo_coeff, matrix_h, &
eigenvalues_ot, do_rotation=.FALSE.)
END IF
IF (k > 0) energy_homo = eigenvalues_ot(k)
IF (k > 0) THEN
CPASSERT(nhomo > 0 .AND. nhomo <= k)
energy_homo = eigenvalues_ot(nhomo)
END IF
ELSE
IF (needs_homo) THEN
CPABORT("Not yet implemented")
@ -183,16 +195,24 @@ CONTAINS
my_solver_type = solver_type
preconditioner_env%in_use = precon_type
preconditioner_env%cholesky_use = cholesky_reduce
!dbg
! write(*,*) PRESENT(chol_type)
! IF(PRESENT(chol_type)) write(*,*) chol_type, cholesky_reduce
! stop
!dbg
IF (PRESENT(chol_type)) preconditioner_env%cholesky_use = chol_type
preconditioner_env%in_use = precon_type
CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
energy_homo, eigenvalues_ot, energy_gap, &
my_solver_type)
IF (nhomo == k) THEN
CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_coeff, &
energy_homo, eigenvalues_ot, energy_gap, my_solver_type)
ELSE
CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nhomo, &
context=preconditioner_env%ctxt, &
para_env=preconditioner_env%para_env)
CALL cp_fm_create(mo_occ, fm_struct)
CALL cp_fm_to_fm(mo_coeff, mo_occ, nhomo)
CALL cp_fm_struct_release(fm_struct)
!
CALL make_preconditioner_matrix(preconditioner_env, matrix_h, matrix_s, matrix_t, mo_occ, &
energy_homo, eigenvalues_ot(1:nhomo), energy_gap, my_solver_type)
!
CALL cp_fm_release(mo_occ)
END IF
CALL solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, matrix_h)
@ -206,9 +226,7 @@ CONTAINS
preconditioner_env%para_env, preconditioner_env%ctxt)
END IF
IF (needs_full_spectrum) THEN
DEALLOCATE (eigenvalues_ot)
END IF
DEALLOCATE (eigenvalues_ot)
CALL timestop(handle)
@ -282,14 +300,15 @@ CONTAINS
!> \param has_unit_metric ...
!> \param convert_to_dbcsr ...
!> \param chol_type ...
!> \param full_mo_set ...
! **************************************************************************************************
SUBROUTINE prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, &
ot_preconditioner, prec_type, solver_type, &
energy_gap, nspins, has_unit_metric, &
convert_to_dbcsr, chol_type)
convert_to_dbcsr, chol_type, full_mo_set)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
TYPE(preconditioner_p_type), DIMENSION(:), POINTER :: ot_preconditioner
INTEGER, INTENT(IN) :: prec_type, solver_type
@ -297,13 +316,16 @@ CONTAINS
INTEGER, INTENT(IN) :: nspins
LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, convert_to_dbcsr
INTEGER, INTENT(IN), OPTIONAL :: chol_type
LOGICAL, INTENT(IN), OPTIONAL :: full_mo_set
CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_preconditioner'
CHARACTER(LEN=default_string_length) :: msg
INTEGER :: handle, icall, ispin, n_loops
INTEGER, DIMENSION(5) :: nocc, norb
LOGICAL :: do_co_rotate, my_convert_to_dbcsr, &
my_has_unit_metric, use_mo_coeff_b
my_full_mo_set, my_has_unit_metric, &
use_mo_coeff_b
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: kinetic
@ -317,6 +339,8 @@ CONTAINS
IF (PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric
my_convert_to_dbcsr = .TRUE.
IF (PRESENT(convert_to_dbcsr)) my_convert_to_dbcsr = convert_to_dbcsr
my_full_mo_set = .FALSE.
IF (PRESENT(full_mo_set)) my_full_mo_set = full_mo_set
CALL get_qs_env(qs_env, &
dft_control=dft_control, &
@ -336,6 +360,15 @@ CONTAINS
matrix_t => kinetic(1)%matrix
END IF
! use full set of MOs or just occupied MOs
nocc = 0
norb = 0
IF (my_full_mo_set) THEN
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), homo=nocc(ispin), nmo=norb(ispin))
CALL set_mo_set(mo_set=mos(ispin), homo=norb(ispin))
END DO
END IF
!determines how often make preconditioner is called, spin dependent methods have to be called twice
n_loops = 1
IF (prec_type == ot_precond_full_single_inverse) n_loops = nspins
@ -399,6 +432,13 @@ CONTAINS
END DO
END SELECT
! reset homo values
IF (my_full_mo_set) THEN
DO ispin = 1, nspins
CALL set_mo_set(mo_set=mos(ispin), homo=nocc(ispin))
END DO
END IF
CALL timestop(handle)
END SUBROUTINE prepare_preconditioner

6
src/preconditioner_makes.F
View file

@ -83,7 +83,7 @@ CONTAINS
TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t
TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
REAL(KIND=dp) :: energy_homo
REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_ot
REAL(KIND=dp), DIMENSION(:) :: eigenvalues_ot
REAL(KIND=dp) :: energy_gap
INTEGER :: my_solver_type
@ -377,7 +377,7 @@ CONTAINS
TYPE(preconditioner_type) :: preconditioner_env
TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
TYPE(dbcsr_type), POINTER :: matrix_h, matrix_s
REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals
REAL(KIND=dp), DIMENSION(:) :: c0_evals
REAL(KIND=dp) :: energy_gap
CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all'
@ -552,7 +552,7 @@ CONTAINS
TYPE(preconditioner_type) :: preconditioner_env
TYPE(cp_fm_type), INTENT(IN) :: matrix_c0
TYPE(dbcsr_type), POINTER :: matrix_h
REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals
REAL(KIND=dp), DIMENSION(:) :: c0_evals
REAL(KIND=dp) :: energy_gap
CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all_ortho'

35
src/qs_density_matrices.F
View file

@ -368,25 +368,40 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_wz_matrix'
INTEGER :: handle, ncol, nrow
INTEGER :: handle, ncol, nocc, nrow
TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
TYPE(cp_fm_type) :: ksmat, scrv
CALL timeset(routineN, handle)
! CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
! CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors")
! CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
! para_env=mo_set%mo_coeff%matrix_struct%para_env, &
! context=mo_set%mo_coeff%matrix_struct%context)
! CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
! CALL cp_fm_struct_release(fm_struct_tmp)
! CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, ncol)
! CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
! CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
! CALL dbcsr_set(w_matrix, 0.0_dp)
! CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, &
! ncol=mo_set%homo, symmetry_mode=1)
! CALL cp_fm_release(scrv)
! CALL cp_fm_release(ksmat)
CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
nocc = mo_set%homo
CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors")
CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nocc, ncol_global=nocc, &
para_env=mo_set%mo_coeff%matrix_struct%para_env, &
context=mo_set%mo_coeff%matrix_struct%context)
CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
CALL cp_fm_struct_release(fm_struct_tmp)
CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, ncol)
CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, nocc)
CALL parallel_gemm("T", "N", nocc, nocc, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
CALL parallel_gemm("N", "N", nrow, nocc, nocc, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
CALL dbcsr_set(w_matrix, 0.0_dp)
CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, &
ncol=mo_set%homo, symmetry_mode=1)
CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, ncol=nocc, symmetry_mode=1)
CALL cp_fm_release(scrv)
CALL cp_fm_release(ksmat)
@ -401,15 +416,17 @@ CONTAINS
!> \param hzm ...
!> \param w_matrix sparse matrix
!> \param focc ...
!> \param nocc ...
!> \par History
!> adapted from calculate_w_matrix_1
!> \author JGH
! **************************************************************************************************
SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc)
SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
TYPE(cp_fm_type), INTENT(IN) :: c0vec
TYPE(dbcsr_type), POINTER :: hzm, w_matrix
REAL(KIND=dp), INTENT(IN) :: focc
INTEGER, INTENT(IN) :: nocc
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_whz_matrix'
@ -424,6 +441,8 @@ CONTAINS
CALL cp_fm_create(hcvec, c0vec%matrix_struct, "hcvec")
CALL cp_fm_get_info(hcvec, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb)
CPASSERT(nocc <= norb .AND. nocc > 0)
norb = nocc
CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
ncol_global=norb, para_env=fm_struct%para_env)
CALL cp_fm_create(chcmat, fm_struct_mat)

68
src/qs_linres_methods.F
View file

@ -199,6 +199,9 @@ CONTAINS
!> \brief scf loop to optimize the first order wavefunctions (psi1)
!> given a perturbation as an operator applied to the ground
!> state orbitals (h1_psi0)
!> psi1 is defined wrt psi0_order (can be a subset of the occupied space)
!> The reference ground state is defined through qs_env (density and ground state MOs)
!> psi1 is orthogonal to the occupied orbitals in the ground state MO set (qs_env%mos)
!> \param p_env ...
!> \param qs_env ...
!> \param psi1 ...
@ -219,8 +222,8 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'linres_solver'
INTEGER :: handle, ispin, iter, maxnmo, maxnmo_o, &
nao, ncol, nmo, nspins
INTEGER :: handle, ispin, iter, maxnmo, nao, ncol, &
nmo, nocc, nspins
LOGICAL :: restart
REAL(dp) :: alpha, beta, norm_res, t1, t2
REAL(dp), DIMENSION(:), POINTER :: tr_pAp, tr_rz0, tr_rz00, tr_rz1
@ -253,25 +256,18 @@ CONTAINS
mos=mos)
nspins = dft_control%nspins
CALL get_mo_set(mos(1), nao=nao)
CALL cp_fm_get_info(psi1(1), nrow_global=nao)
maxnmo = 0
maxnmo_o = 0
DO ispin = 1, nspins
CALL get_mo_set(mos(ispin), nmo=ncol)
CALL cp_fm_get_info(psi1(ispin), ncol_global=ncol)
maxnmo = MAX(maxnmo, ncol)
CALL cp_fm_get_info(psi0_order(ispin), ncol_global=ncol)
maxnmo_o = MAX(maxnmo_o, ncol)
END DO
!
CALL check_p_env_init(p_env, linres_control, nspins)
!
! allocate the vectors
ALLOCATE (tr_pAp(nspins), tr_rz0(nspins), tr_rz00(nspins), tr_rz1(nspins), &
r(nspins), p(nspins), z(nspins), Ap(nspins), mo_coeff_array(nspins))
DO ispin = 1, nspins
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
mo_coeff_array(ispin) = mo_coeff
END DO
r(nspins), p(nspins), z(nspins), Ap(nspins))
!
DO ispin = 1, nspins
CALL cp_fm_create(r(ispin), psi1(ispin)%matrix_struct)
@ -280,11 +276,24 @@ CONTAINS
CALL cp_fm_create(Ap(ispin), psi1(ispin)%matrix_struct)
END DO
!
! compute S*C0, C0_order'*H*C0_order (this should be done once for all)
ALLOCATE (chc(nspins), Sc(nspins))
! build C0 occupied vectors and S*C0 matrix
ALLOCATE (Sc(nspins), mo_coeff_array(nspins))
DO ispin = 1, nspins
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
CALL cp_fm_create(Sc(ispin), mo_coeff%matrix_struct)
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, homo=nocc)
NULLIFY (tmp_fm_struct)
CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
ncol_global=nocc, para_env=para_env, &
context=mo_coeff%matrix_struct%context)
CALL cp_fm_create(mo_coeff_array(ispin), tmp_fm_struct)
CALL cp_fm_to_fm(mo_coeff, mo_coeff_array(ispin), nocc)
CALL cp_fm_create(Sc(ispin), tmp_fm_struct)
CALL cp_fm_struct_release(tmp_fm_struct)
END DO
!
! Allocate C0_order'*H*C0_order
ALLOCATE (chc(nspins))
DO ispin = 1, nspins
CALL cp_fm_get_info(psi1(ispin), ncol_global=nmo)
NULLIFY (tmp_fm_struct)
CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
ncol_global=nmo, para_env=para_env, &
@ -305,9 +314,8 @@ CONTAINS
END ASSOCIATE
!
! S * C0
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
CALL cp_fm_get_info(mo_coeff, ncol_global=ncol)
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, Sc(ispin), ncol)
CALL cp_fm_get_info(mo_coeff_array(ispin), ncol_global=ncol)
CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff_array(ispin), Sc(ispin), ncol)
END DO
!
! header
@ -384,7 +392,7 @@ CONTAINS
CALL cp_fm_trace(r(ispin), z(ispin), tr_rz0(ispin))
END DO
IF (SUM(tr_rz0) < 0.0_dp) CPABORT("tr(r_j*z_j) < 0")
norm_res = ABS(SUM(tr_rz0))/SQRT(REAL(nspins*nao*maxnmo_o, dp))
norm_res = ABS(SUM(tr_rz0))/SQRT(REAL(nspins*nao*maxnmo, dp))
!
alpha = 0.0_dp
restart = .FALSE.
@ -463,7 +471,6 @@ CONTAINS
CALL apply_op(qs_env, p_env, psi0_order, psi1, Ap, chc)
!
DO ispin = 1, nspins
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
CALL cp_fm_to_fm(h1_psi0(ispin), r(ispin))
CALL cp_fm_scale_and_add(-1.0_dp, r(ispin), -1.0_dp, Ap(ispin))
END DO
@ -506,7 +513,7 @@ CONTAINS
CALL cp_fm_trace(r(ispin), z(ispin), tr_rz1(ispin))
END DO
IF (SUM(tr_rz1) < 0.0_dp) CPABORT("tr(r_j+1*z_j+1) < 0")
norm_res = SUM(tr_rz1)/SQRT(REAL(nspins*nao*maxnmo_o, dp))
norm_res = SUM(tr_rz1)/SQRT(REAL(nspins*nao*maxnmo, dp))
!
! beta = tr(r_j+1*z_j+1) / tr(r_j*z_j)
IF (SUM(tr_rz0) < 1.0e-10_dp) THEN
@ -539,6 +546,7 @@ CONTAINS
CALL cp_fm_release(z(ispin))
CALL cp_fm_release(Ap(ispin))
!
CALL cp_fm_release(mo_coeff_array(ispin))
CALL cp_fm_release(Sc(ispin))
CALL cp_fm_release(chc(ispin))
END DO
@ -565,7 +573,8 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op'
INTEGER :: handle, ispin, nspins
INTEGER :: handle, ispin, nc1, nc2, nc3, nc4, nr1, &
nr2, nr3, nr4, nspins
REAL(dp) :: chksum
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
TYPE(dft_control_type), POINTER :: dft_control
@ -583,6 +592,19 @@ CONTAINS
nspins = dft_control%nspins
DO ispin = 1, nspins
!c0, v, Av, chc
CALL cp_fm_get_info(c0(ispin), ncol_global=nc1, nrow_global=nr1)
CALL cp_fm_get_info(v(ispin), ncol_global=nc2, nrow_global=nr2)
CALL cp_fm_get_info(Av(ispin), ncol_global=nc3, nrow_global=nr3)
CALL cp_fm_get_info(chc(ispin), ncol_global=nc4, nrow_global=nr4)
IF (.NOT. (nc1 == nc2 .AND. nr1 == nr2 .AND. nc1 == nc3 .AND. nr1 == nr3 &
.AND. nc4 == nr4 .AND. nc1 <= nc4)) THEN
CALL cp_abort(__LOCATION__, &
"Number of vectors inconsistent or CHC matrix too small")
END IF
END DO
! apply the uncoupled operator
DO ispin = 1, nspins
CALL apply_op_1(v(ispin), Av(ispin), matrix_ks(ispin)%matrix, &

3
src/qs_scf_diagonalization.F
View file

@ -1641,7 +1641,8 @@ CONTAINS
scf_env%block_davidson_env(1)%prec_type, &
scf_env%block_davidson_env(1)%solver_type, &
scf_env%block_davidson_env(1)%energy_gap, nspins, &
convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos)
convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos, &
full_mo_set=.TRUE.)
END IF
DO ispin = 1, nspins

49
src/response_solver.F
View file

@ -32,12 +32,15 @@ MODULE response_solver
cp_dbcsr_sm_fm_multiply,&
dbcsr_allocate_matrix_set,&
dbcsr_deallocate_matrix_set
USE cp_fm_struct, ONLY: cp_fm_struct_type
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_get_info,&
cp_fm_init_random,&
cp_fm_release,&
cp_fm_set_all,&
cp_fm_to_fm,&
cp_fm_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_unit_nr,&
@ -546,7 +549,7 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'response_equation_new'
INTEGER :: handle, ispin, nao, nao_aux, nspins
INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
LOGICAL :: should_stop
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_struct_type), POINTER :: fm_struct
@ -555,13 +558,14 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
CALL timeset(routineN, handle)
NULLIFY (dft_control, matrix_ks, mo_coeff, mos)
CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks, &
matrix_s=matrix_s, mos=mos)
para_env=para_env, matrix_s=matrix_s, mos=mos)
nspins = dft_control%nspins
! Initialize vectors:
@ -569,11 +573,16 @@ CONTAINS
! psi1 : The "perturbed" linear response orbitals
ALLOCATE (psi0(nspins), psi1(nspins))
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
psi0(ispin) = mo_coeff
CALL cp_fm_get_info(mo_coeff, matrix_struct=fm_struct)
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, homo=nocc)
NULLIFY (fm_struct)
CALL cp_fm_struct_create(fm_struct, nrow_global=nao, &
ncol_global=nocc, para_env=para_env, &
context=mo_coeff%matrix_struct%context)
CALL cp_fm_create(psi0(ispin), fm_struct)
CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc)
CALL cp_fm_create(psi1(ispin), fm_struct)
CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
CALL cp_fm_struct_release(fm_struct)
END DO
should_stop = .FALSE.
@ -617,6 +626,7 @@ CONTAINS
DO ispin = 1, nspins
CALL cp_fm_release(cpmos(ispin))
CALL cp_fm_release(psi1(ispin))
CALL cp_fm_release(psi0(ispin))
END DO
DEALLOCATE (psi0, psi1)
@ -643,7 +653,7 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'response_equation'
INTEGER :: handle, ispin, nao, nao_aux, nspins
INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
LOGICAL :: should_stop
TYPE(admm_type), POINTER :: admm_env
TYPE(cp_fm_struct_type), POINTER :: fm_struct
@ -653,6 +663,7 @@ CONTAINS
TYPE(dft_control_type), POINTER :: dft_control
TYPE(linres_control_type), POINTER :: linres_control
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
@ -688,7 +699,7 @@ CONTAINS
CALL p_env_psi0_changed(p_env, qs_env)
p_env%new_preconditioner = .TRUE.
CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env, mos=mos)
!
nspins = dft_control%nspins
@ -697,11 +708,16 @@ CONTAINS
! psi1 : The "perturbed" linear response orbitals
ALLOCATE (psi0(nspins), psi1(nspins))
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
psi0(ispin) = mo_coeff
CALL cp_fm_get_info(mo_coeff, matrix_struct=fm_struct)
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nao=nao, homo=nocc)
NULLIFY (fm_struct)
CALL cp_fm_struct_create(fm_struct, nrow_global=nao, &
ncol_global=nocc, para_env=para_env, &
context=mo_coeff%matrix_struct%context)
CALL cp_fm_create(psi0(ispin), fm_struct)
CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc)
CALL cp_fm_create(psi1(ispin), fm_struct)
CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
CALL cp_fm_struct_release(fm_struct)
END DO
should_stop = .FALSE.
@ -765,6 +781,7 @@ CONTAINS
p_env%w1(ispin)%matrix)
END DO
DO ispin = 1, nspins
CALL cp_fm_release(psi0(ispin))
CALL cp_fm_release(psi1(ispin))
END DO
DEALLOCATE (psi0, psi1)
@ -811,7 +828,7 @@ CONTAINS
INTEGER :: handle, iounit, ispin, mspin, n_rep_hf, &
nao, nao_aux, natom, nder, nimages, &
nspins
nocc, nspins
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: debug_forces, debug_stress, distribute_fock_matrix, do_ex, do_hfx, gapw, gapw_xc, &
hfx_treat_lsd_in_core, resp_only, s_mstruct_changed, use_virial
@ -2036,8 +2053,9 @@ CONTAINS
IF (nspins == 1) focc = 2.0_dp
CALL get_qs_env(qs_env, mos=mos)
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
matrix_wz(ispin)%matrix, focc)
matrix_wz(ispin)%matrix, focc, nocc)
END DO
END IF
IF (nspins == 2) THEN
@ -2110,7 +2128,7 @@ CONTAINS
INTEGER :: atom_a, handle, iatom, ikind, iounit, &
is, ispin, na, natom, natorb, nimages, &
nkind, ns, nsgf, nspins
nkind, nocc, ns, nsgf, nspins
INTEGER, DIMENSION(25) :: lao
INTEGER, DIMENSION(5) :: occ
LOGICAL :: debug_forces, do_ex, use_virial
@ -2259,8 +2277,9 @@ CONTAINS
IF (nspins == 1) focc = 1.0_dp
CALL get_qs_env(qs_env, mos=mos)
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
matrix_wz(ispin)%matrix, focc)
matrix_wz(ispin)%matrix, focc, nocc)
END DO
IF (nspins == 2) THEN
CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &

5
src/rpa_im_time_force_methods.F
View file

@ -3003,7 +3003,7 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'update_im_time_forces'
INTEGER :: handle, i, idens, ispin, n_rep_hf, nao, &
nao_aux, nder, nimages, nspins
nao_aux, nder, nimages, nocc, nspins
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: do_exx, do_hfx, do_tau, do_tau_admm, &
use_virial
@ -3639,8 +3639,9 @@ CONTAINS
IF (nspins == 2) focc = 1.0_dp
CALL get_qs_env(qs_env, mos=mos)
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
p_env%w1(ispin)%matrix, focc)
p_env%w1(ispin)%matrix, focc, nocc)
END DO
IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, 1.0_dp)

2
tests/QS/regtest-mp2-grad-2/TEST_FILES
View file

@ -2,7 +2,7 @@ O2_dyn.inp 11 2e-11 -
O2_dyn_mme.inp 72 1e-8 0.50585576
CH_dyn_screen.inp 11 2e-04 -5.700371305153854
CH_dyn_screen_single_group.inp 11 2e-04 -5.700371305153854
MOM_MP2_geoopt.inp 11 4e-06 -13.969926262647746
MOM_MP2_geoopt.inp 11 4e-06 -13.969416309750072
H2O_MD_mme.inp 11 1e-10 -17.056807598425291
O2_dyn_ri-hfx.inp 11 1e-10 -31.518511363027905
#EOF

10
tests/TEST_DIRS
View file

@ -120,6 +120,8 @@ QS/regtest-gpw-3
QS/regtest-ps-implicit-1-3 fftw3
QS/regtest-ps-implicit-2-2
QS/regtest-rtp-3
QS/regtest-rtp-4
QS/regtest-rtp-5
QS/regtest-pao-1
QMMM/QS/regtest-1
QS/regtest-ps-implicit-2-1
@ -240,6 +242,8 @@ QS/regtest-hybrid-2 libint
QS/regtest-chi-1
QS/regtest-cdft-4-2
DFTB/regtest-nonscc
QS/regtest-as libint
QS/regtest-as-dft libint
QS/regtest-double-hybrid-3 libint
QS/regtest-dft-filtermat
SE/regtest-3-3
@ -285,8 +289,6 @@ QS/regtest-hfx-2 libint
QS/regtest-nmr-3
QS/regtest-hybrid-1 libint
QS/regtest-rtp-1
QS/regtest-as libint
QS/regtest-as-dft libint
QS/regtest-tddfpt-prop
QS/regtest-sym-3
QS/regtest-gapw-ext
@ -297,6 +299,7 @@ QS/regtest-sym-1
QS/regtest-sym-2
QS/regtest-ot-refine-2
NEB/regtest-2
QS/regtest-lsroks libint
QS/regtest-nmr-2
QS/regtest-cdft-3-1
QS/regtest-optbas
@ -320,7 +323,6 @@ QS/regtest-pao-3
DFTB/regtest-debug
optimize_input/regtest-1
FE/regtest-2
QS/regtest-lsroks libint
NEB/regtest-3
QS/regtest-linearscaling
SE/regtest-3-1
@ -335,7 +337,5 @@ Fist/regtest-plumed2 plumed2
Fist/regtest-quip quip
Fist/regtest-nequip libtorch
Fist/regtest-16
QS/regtest-rtp-4
QS/regtest-rtp-5
DFTB/regtest-vdw
QS/regtest-tddfpt-force-gapw