Skip to content

Commit

Permalink
Fixed orbital angular momentum calculation for RTP: (#139)
Browse files Browse the repository at this point in the history 
* Fixed orbital angular momentum calculation for RTP:
1. use non-symmetric neighbour list in qs_moments::build_local_magmom_matrix
2. track imaginary part of the density if magnetic moments are requested
3. use imaginary part of the density to calculate the expectation value
4. fix regtest-ot-1/H20-magnetic.inp
  • Loading branch information
@mattiatj @dev-zero
mattiatj authored and dev-zero committed Dec 21, 2018
1 parent 4ed24fb commit f20b940
Show file tree
Hide file tree
Showing 5 changed files with 67 additions and 56 deletions.
6 changes: 3 additions & 3 deletions src/emd/rt_propagation_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -396,7 +396,7 @@ SUBROUTINE calc_update_rho(qs_env)
END DO
CALL qs_rho_update_rho(rho, qs_env)

IF (rtp%do_hfx) THEN
IF (rtp%do_hfx .OR. rtp%magnetic) THEN
CALL qs_rho_get(rho_struct=rho, rho_ao_im=rho_ao_im)
CALL calculate_P_imaginary(qs_env, rtp, rho_ao_im, keep_sparsity=.TRUE.)
CALL qs_rho_set(rho, rho_ao_im=rho_ao_im)
Expand Down Expand Up @@ -440,11 +440,11 @@ SUBROUTINE calc_update_rho_sparse(qs_env)
rtp_control => dft_control%rtp_control
CALL get_rtp(rtp=rtp, rho_new=rho_new)
CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao)
IF (rtp%do_hfx) CALL qs_rho_get(rho_struct=rho, rho_ao_im=rho_ao_im)
IF (rtp%do_hfx .OR. rtp%magnetic) CALL qs_rho_get(rho_struct=rho, rho_ao_im=rho_ao_im)
DO ispin = 1, SIZE(rho_ao)
CALL dbcsr_set(rho_ao(ispin)%matrix, zero)
CALL dbcsr_copy_into_existing(rho_ao(ispin)%matrix, rho_new(ispin*2-1)%matrix)
IF (rtp%do_hfx) CALL dbcsr_copy_into_existing(rho_ao_im(ispin)%matrix, rho_new(ispin*2)%matrix)
IF (rtp%do_hfx .OR. rtp%magnetic) CALL dbcsr_copy_into_existing(rho_ao_im(ispin)%matrix, rho_new(ispin*2)%matrix)
END DO

CALL qs_rho_update_rho(rho, qs_env)
Expand Down
7 changes: 5 additions & 2 deletions src/motion/rt_propagation.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,8 @@ SUBROUTINE rt_prop_setup(force_env)
TYPE(rt_prop_type), POINTER :: rtp
TYPE(rtp_control_type), POINTER :: rtp_control
TYPE(section_vals_type), POINTER :: hfx_sections, input, ls_scf_section, &
md_section, motion_section
md_section, motion_section, &
print_moments_section
NULLIFY (qs_env, rtp_control, dft_control)
Expand All @@ -129,6 +130,7 @@ SUBROUTINE rt_prop_setup(force_env)
motion_section => section_vals_get_subs_vals(force_env%root_section, "MOTION")
md_section => section_vals_get_subs_vals(motion_section, "MD")
hfx_sections => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%DFT%XC%HF")
print_moments_section => section_vals_get_subs_vals(force_env%root_section, "FORCE_EVAL%DFT%PRINT%MOMENTS")
CALL section_vals_val_get(md_section, "TIMESTEP", r_val=qs_env%rtp%dt)
CALL section_vals_val_get(md_section, "STEP_START_VAL", i_val=qs_env%rtp%i_start)
CALL section_vals_val_get(md_section, "STEPS", i_val=rtp%nsteps)
Expand All @@ -143,10 +145,11 @@ SUBROUTINE rt_prop_setup(force_env)
CALL section_vals_val_get(ls_scf_section, "MAX_ITER_LANCZOS", i_val=rtp%lanzcos_max_iter)
CALL section_vals_val_get(ls_scf_section, "SIGN_SQRT_ORDER", i_val=rtp%newton_schulz_order)
CALL section_vals_get(hfx_sections, explicit=rtp%do_hfx)
CALL section_vals_val_get(print_moments_section, "MAGNETIC", l_val=rtp%magnetic)
! Hmm, not really like to initialize with the structure of S but I reckon it is
! done everywhere like this
IF (rtp%do_hfx) &
IF (rtp%do_hfx .OR. rtp%magnetic) &
CALL rtp_hfx_rebuild(qs_env)
CALL init_propagation_run(qs_env)
Expand Down
106 changes: 56 additions & 50 deletions src/qs_moments.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,8 +47,9 @@ MODULE qs_moments
put_results
USE cp_result_types, ONLY: cp_result_type
USE dbcsr_api, ONLY: &
dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, &
dbcsr_get_block_p, dbcsr_p_type, dbcsr_set, dbcsr_trace, dbcsr_type, dbcsr_type_symmetric
dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_multiply, dbcsr_p_type, dbcsr_set, &
dbcsr_trace, dbcsr_type, dbcsr_type_symmetric
USE distribution_1d_types, ONLY: distribution_1d_type
USE kinds, ONLY: default_string_length,&
dp
Expand Down Expand Up @@ -364,7 +365,7 @@ SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_po
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
POINTER :: sab_all
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
Expand All @@ -373,7 +374,7 @@ SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_po
CALL timeset(routineN, handle)
NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
NULLIFY (qs_kind_set, cell, particle_set, sab_all)
NULLIFY (matrix_s)
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
Expand All @@ -384,15 +385,15 @@ SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_po
DO i = 1, nm
ALLOCATE (magmom(i)%matrix)
CALL dbcsr_copy(magmom(i)%matrix, matrix_s(1)%matrix, "Moments")
CALL dbcsr_desymmetrize(matrix_s(1)%matrix, magmom(i)%matrix)
CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
END DO
NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
NULLIFY (qs_kind_set, particle_set, sab_all, cell)
CALL get_qs_env(qs_env=qs_env, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, cell=cell, &
sab_orb=sab_orb)
sab_all=sab_all)
nkind = SIZE(qs_kind_set)
natom = SIZE(particle_set)
Expand Down Expand Up @@ -422,7 +423,7 @@ SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_po
NULLIFY (basis_set_list(ikind)%gto_basis_set)
END IF
END DO
CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
CALL neighbor_list_iterator_create(nl_iterator, sab_all)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
iatom=iatom, jatom=jatom, r=rab)
Expand Down Expand Up @@ -462,19 +463,15 @@ SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_po
last_jatom = jatom
IF (iatom <= jatom) THEN
irow = iatom
icol = jatom
ELSE
irow = jatom
icol = iatom
END IF
irow = iatom
icol = jatom
DO i = 1, nm
NULLIFY (mint(i)%block)
CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
row=irow, col=icol, BLOCK=mint(i)%block, found=found)
mint(i)%block = 0._dp
CPASSERT(ASSOCIATED(mint(i)%block))
END DO
! fold atomic position back into unit cell
Expand Down Expand Up @@ -522,23 +519,11 @@ SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_po
sphi_b(1, sgfb), SIZE(sphi_b, 1), &
0.0_dp, work(1, 1), SIZE(work, 1))

IF (iatom <= jatom) THEN

CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
work(1, 1), SIZE(work, 1), &
1.0_dp, mint(i)%block(sgfa, sgfb), &
SIZE(mint(i)%block, 1))

ELSE

CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1.0_dp, work(1, 1), SIZE(work, 1), &
sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1.0_dp, mint(i)%block(sgfb, sgfa), &
SIZE(mint(i)%block, 1))

END IF
CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
work(1, 1), SIZE(work, 1), &
1.0_dp, mint(i)%block(sgfa, sgfb), &
SIZE(mint(i)%block, 1))

END DO

Expand Down Expand Up @@ -1495,6 +1480,7 @@ SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, uni
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(cp_result_type), POINTER :: results
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, rho_ao
TYPE(dbcsr_type), POINTER :: rho_magmom_ao
TYPE(dft_control_type), POINTER :: dft_control
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
Expand Down Expand Up @@ -1559,23 +1545,6 @@ SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, uni

CALL dbcsr_deallocate_matrix_set(moments)

! magnetic moments
IF (magnetic) THEN
NULLIFY (magmom)
CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
nm = SIZE(magmom)
ALLOCATE (mmom(nm))
DO i = 1, SIZE(magmom)
strace = 0._dp
DO ispin = 1, dft_control%nspins
CALL dbcsr_trace(rho_ao(ispin)%matrix, magmom(i)%matrix, trace)
strace = strace+trace
END DO
mmom(i) = strace
END DO
CALL dbcsr_deallocate_matrix_set(magmom)
END IF

! nuclear contribution
CALL get_qs_env(qs_env=qs_env, &
local_particles=local_particles)
Expand Down Expand Up @@ -1606,6 +1575,43 @@ SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, uni
rmom(:, :) = -rmom(:, :)
rmom(:, 3) = rmom(:, 1)+rmom(:, 2)

! magnetic moments
IF (magnetic) THEN
NULLIFY (magmom)
CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
nm = SIZE(magmom)
ALLOCATE (mmom(nm))

! Allocate matrices to store the matrix product to be traced (dbcsr_trace only works for products of
! symmetric matrices)
NULLIFY (rho_magmom_ao)
ALLOCATE (rho_magmom_ao)
CALL dbcsr_desymmetrize(matrix_s(1)%matrix, rho_magmom_ao)

IF (qs_env%run_rtp) THEN
! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
! There may be other cases, where the imaginary part of the density is relevant
NULLIFY (rho_ao)
CALL qs_rho_get(rho, rho_ao_im=rho_ao)
ENDIF
! if the density is purely real this is an expensive way to calculate zero
DO i = 1, SIZE(magmom)
strace = 0._dp
DO ispin = 1, dft_control%nspins
CALL dbcsr_set(rho_magmom_ao, 0.0_dp)
CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
0.0_dp, rho_magmom_ao)
CALL dbcsr_trace(rho_magmom_ao, trace)
strace = strace+trace
END DO
mmom(i) = strace
END DO

CALL dbcsr_deallocate_matrix_set(magmom)
CALL dbcsr_deallocate_matrix(rho_magmom_ao)
END IF

description = "[DIPOLE]"
CALL cp_results_erase(results=results, description=description)
CALL put_results(results=results, description=description, &
Expand Down Expand Up @@ -1736,7 +1742,7 @@ SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmo
IF (PRESENT(mmom)) THEN
IF (nmom >= 1) THEN
dd = SQRT(SUM(mmom(1:3)**2))
WRITE (unit_number, "(T3,A)") "Magnetic Dipole Moment [WRONG]"
WRITE (unit_number, "(T3,A)") "Magnetic Dipole Moment (only orbital contrib.) [a.u.]"
WRITE (unit_number, "(T5,3(A,A,F14.8,1X),T60,A,T68,F14.8)") &
(TRIM(rlab(i+1)), "=", mmom(i), i=1, 3), "Total=", dd
END IF
Expand Down
2 changes: 2 additions & 0 deletions src/rt_propagation_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,7 @@ MODULE rt_propagation_types
REAL(KIND=dp) :: mixing_factor
LOGICAL :: mixing
LOGICAL :: do_hfx
LOGICAL :: magnetic
INTEGER, DIMENSION(:, :), ALLOCATABLE :: orders
INTEGER :: nsteps, istep, i_start, max_steps
INTEGER :: iter
Expand Down Expand Up @@ -282,6 +283,7 @@ SUBROUTINE rt_prop_create(rtp, mos, mpools, dft_control, template, linear_scalin
rtp%istep = 0
rtp%iter = 0
rtp%do_hfx = .FALSE.
rtp%magnetic = .FALSE.

END SUBROUTINE rt_prop_create

Expand Down
2 changes: 1 addition & 1 deletion tests/QS/regtest-ot-1/TEST_FILES
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,5 +50,5 @@ H2O-broyden-2.inp 1 2e-05 -
#inverse_update_preconditioner solve
H2O-inverse_up.inp 1 1.0E-14 -17.14017488226970
#magnetic dipole moment
H2O-magnetic.inp 17 1.0E-14 0.65963143
H2O-magnetic.inp 17 1.0E-14 0.00000000
#EOF

0 comments on commit f20b940

Please sign in to comment.