Skip to content

Commit

Permalink
task_list: Simplify matrix symmetry handling in integrate_v_rspace
Browse files Browse the repository at this point in the history
  • Loading branch information
@oschuett
oschuett committed Oct 3, 2020
1 parent a5986e6 commit 2ab8ed7
Show file tree
Hide file tree
Showing 4 changed files with 47 additions and 80 deletions.
117 changes: 42 additions & 75 deletions src/qs_integrate_potential_product.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,7 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env
pab_required, use_subpatch, use_virial
REAL(KIND=dp) :: admm_scal_fac, eps_rho_rspace, f, &
prefactor, radius, scalef, zetp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_block_set, p_block_set
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab, rab_inv, rb, &
rp
REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b, pv_thread
Expand Down Expand Up @@ -469,8 +470,8 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env
!$OMP PRIVATE(img,brow,bcol,orb_basis_set,first_sgfa,la_max,la_min,npgfa,nseta,nsgfa), &
!$OMP PRIVATE(rpgfa,set_radius_a,sphi_a,zeta,first_sgfb,lb_max,lb_min,npgfb), &
!$OMP PRIVATE(nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,found), &
!$OMP PRIVATE(force_a,force_b,my_virial_a,my_virial_b,atom_pair_changed,h_block), &
!$OMP PRIVATE(p_block,ncoa,sgfa,ncob,sgfb,rab,ra,rb,rp,zetp,f,prefactor,radius,igrid_level), &
!$OMP PRIVATE(force_a,force_b,my_virial_a,my_virial_b,atom_pair_changed,h_block,h_block_set), &
!$OMP PRIVATE(p_block,p_block_set,ncoa,sgfa,ncob,sgfb,rab,ra,rb,rp,zetp,f,prefactor,radius,igrid_level), &
!$OMP PRIVATE(na1,na2,nb1,nb2,use_subpatch,rab_inv,new_set_pair_coming,atom_pair_done), &
!$OMP PRIVATE(iset_new,jset_new,ipgf_new,jpgf_new,scalef), &
!$OMP PRIVATE(itask,dist,has_threads) &
Expand Down Expand Up @@ -621,24 +622,21 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env
sgfa = first_sgfa(1, iset)
ncob = npgfb(jset)*ncoset(lb_max(jset))
sgfb = first_sgfb(1, jset)

IF (ALLOCATED(p_block_set)) DEALLOCATE (p_block_set)
ALLOCATE (p_block_set(nsgfa(iset), nsgfb(jset)))
IF (pab_required) THEN
IF (iatom <= jatom) THEN
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), &
p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1))
pab(1:ncoa, 1:ncob) = MATMUL(work(1:ncoa, 1:nsgfb(jset)), TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
p_block_set(:, :) = p_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1)
ELSE
work(1:ncob, 1:nsgfa(iset)) = MATMUL(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1), &
p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1))
pab(1:ncob, 1:ncoa) = MATMUL(work(1:ncob, 1:nsgfa(iset)), TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)))
END IF
p_block_set(:, :) = TRANSPOSE(p_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1))
ENDIF
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1), p_block_set)
pab(1:ncoa, 1:ncob) = MATMUL(work(1:ncoa, 1:nsgfb(jset)), &
TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)))
END IF

IF (iatom <= jatom) THEN
hab(1:ncoa, 1:ncob) = 0._dp
ELSE
hab(1:ncob, 1:ncoa) = 0._dp
ENDIF

hab(1:ncoa, 1:ncob) = 0._dp
iset_old = iset
jset_old = jset

Expand Down Expand Up @@ -675,61 +673,30 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env
ENDIF
IF (pab_required) THEN
IF (iatom <= jatom) THEN
CALL integrate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, rs_v(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, pab=pab, o1=na1 - 1, o2=nb1 - 1, &
radius=radius, &
calculate_forces=calculate_forces, &
force_a=force_a, force_b=force_b, &
compute_tau=my_compute_tau, &
use_virial=use_virial, my_virial_a=my_virial_a, &
my_virial_b=my_virial_b, use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
ELSE
rab_inv = -rab
CALL integrate_pgf_product( &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
rb, rab_inv, rs_v(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, pab=pab, o1=nb1 - 1, o2=na1 - 1, &
radius=radius, &
calculate_forces=calculate_forces, &
force_a=force_b, force_b=force_a, &
compute_tau=my_compute_tau, &
use_virial=use_virial, my_virial_a=my_virial_b, &
my_virial_b=my_virial_a, use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
END IF
CALL integrate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, rs_v(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, pab=pab, o1=na1 - 1, o2=nb1 - 1, &
radius=radius, &
calculate_forces=calculate_forces, &
force_a=force_a, force_b=force_b, &
compute_tau=my_compute_tau, &
use_virial=use_virial, my_virial_a=my_virial_a, &
my_virial_b=my_virial_b, use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
ELSE
IF (iatom <= jatom) THEN
CALL integrate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, rs_v(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, o1=na1 - 1, o2=nb1 - 1, &
radius=radius, &
calculate_forces=calculate_forces, &
force_a=force_a, force_b=force_b, &
compute_tau=my_compute_tau, &
use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
ELSE
rab_inv = -rab
CALL integrate_pgf_product( &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
rb, rab_inv, rs_v(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, o1=nb1 - 1, o2=na1 - 1, &
radius=radius, &
calculate_forces=calculate_forces, &
force_a=force_b, force_b=force_a, &
compute_tau=my_compute_tau, &
use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
END IF
CALL integrate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
ra, rab, rs_v(igrid_level)%rs_grid, cell, &
cube_info(igrid_level), &
hab, o1=na1 - 1, o2=nb1 - 1, &
radius=radius, &
calculate_forces=calculate_forces, &
force_a=force_a, force_b=force_b, &
compute_tau=my_compute_tau, &
use_subpatch=use_subpatch, subpatch_pattern=tasks(itask)%subpatch_pattern)
END IF
new_set_pair_coming = .FALSE.
Expand All @@ -754,16 +721,16 @@ SUBROUTINE integrate_v_rspace_low(v_rspace, hmat, hmat_kp, pmat, pmat_kp, qs_env
! contract the block into h if we're done with the current set pair
IF (new_set_pair_coming) THEN
IF (ALLOCATED(h_block_set)) DEALLOCATE (h_block_set)
ALLOCATE (h_block_set(nsgfa(iset), nsgfb(jset)))
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab(1:ncoa, 1:ncob), sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
h_block_set(:, :) = MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
IF (iatom <= jatom) THEN
work(1:ncoa, 1:nsgfb(jset)) = MATMUL(hab(1:ncoa, 1:ncob), sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1))
h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) = &
h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + &
MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), work(1:ncoa, 1:nsgfb(jset)))
h_block(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfb(jset) - 1) + h_block_set
ELSE
work(1:ncob, 1:nsgfa(iset)) = MATMUL(hab(1:ncob, 1:ncoa), sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1))
h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + &
MATMUL(TRANSPOSE(sphi_b(1:ncob, sgfb:sgfb + nsgfb(jset) - 1)), work(1:ncob, 1:nsgfa(iset)))
h_block(sgfb:sgfb + nsgfb(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) + TRANSPOSE(h_block_set)
END IF
END IF

Expand Down
4 changes: 2 additions & 2 deletions tests/QS/regtest-kp-1/TEST_FILES
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,8 @@ c_3.inp 1 2e-14
c_4.inp 1 3e-14 -45.66762356227670
c_5.inp 1 3e-14 -46.02680478900018
c_6.inp 1 3e-14 -45.65758013574263
c_gapw.inp 1 2e-14 -302.06260708586098
c_gapwxc.inp 1 2e-14 -45.66746891467081
c_gapw.inp 1 2e-13 -302.06260708586098
c_gapwxc.inp 1 8e-13 -45.66746891467081
cn_1.inp 1 6e-13 -49.27303752713806
c_dos.inp 1 1.0E-14 -45.65699752126638
#EOF
4 changes: 2 additions & 2 deletions tests/QS/regtest-ot-refine-2/TEST_FILES
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@ h2o_ot_chol_1.inp 1 9e-14
h2o_ot_chol_diis_1.inp 1 2e-14 -17.13847256795214
h2o_ot_lwdn_1.inp 1 1e-13 -17.09607220197871
h2o_ot_lwdn_diis_1.inp 1 3e-14 -17.13832134601361
o2_ot_chol_1.inp 1 6e-14 -31.70785812513002
o2_ot_chol_1.inp 1 2e-13 -31.70785812513002
o2_ot_chol_diis_1.inp 1 3e-14 -31.82229826058755
o2_ot_lwdn_1.inp 1 5e-14 -31.70373819277881
o2_ot_lwdn_1.inp 1 2e-13 -31.70373819277881
o2_ot_lwdn_diis_1.inp 1 3e-14 -31.81621403029244
h2o_ot_precond_2.inp 1 1.0E-14 -15.99989407904090
h2o_ot_precond_3.inp 1 1.0E-14 -15.99989407904090
Expand Down
2 changes: 1 addition & 1 deletion tests/QS/regtest-ot-refine-3/TEST_FILES
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ o2_ot_poly_on_the_fly_l1.inp 1 2e-13
h2o_ot_refine_3.inp 1 3e-14 -17.08884672065213
h2o_ot_refine_4.inp 1 2e-14 -17.11857283464189
o2_ot_refine_3.inp 1 5e-13 -31.60858514412410
o2_ot_refine_4.inp 1 4e-14 -31.58768162895743
o2_ot_refine_4.inp 1 2e-13 -31.58768162895743
h2o_ot_precond_1.inp 1 4e-14 -16.05766867329444
h2o_ot_precond_2.inp 1 4e-14 -16.06790080254975
h2o_ot_precond_3.inp 1 4e-14 -16.02559849783371
Expand Down

0 comments on commit 2ab8ed7

Please sign in to comment.