From 80a1506e1b7667fd311368c61a629e5c0e27cc9c Mon Sep 17 00:00:00 2001 From: ayakon Date: Thu, 6 Feb 2025 20:30:56 +0900 Subject: [PATCH 1/7] pDOS with MSSFs (still there would be some bug) modified: DiagModule.f90 modified: initial_read_module.f90 --- src/DiagModule.f90 | 207 ++++++++++++++++++++++++++++-------- src/initial_read_module.f90 | 1 + 2 files changed, 164 insertions(+), 44 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index 9f272c3a9..ebfa45dce 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -276,6 +276,7 @@ module DiagModule ! These two flags are used to perform specific calculations in buildK ! (projected DOS and Resta polarisation respectively) logical :: flag_pDOS_buildK = .false. + logical :: flag_WFatomf_buildK = .false. logical :: flag_pol_buildS = .false. !logical :: diagon ! Do we diagonalise or use O(N) ? @@ -502,10 +503,11 @@ subroutine FindEvals(electrons) nkpoints_max, pgid, N_procs_in_pg, & N_kpoints_in_pg use mult_module, only: matH, matS, matK, matM12, SF_to_AtomF_transform, & - matrix_scale, matrix_product_trace, allocate_temp_matrix, free_temp_matrix - use matrix_data, only: Hrange, Srange, aHa_range + matrix_scale, matrix_product, matrix_product_trace, allocate_temp_matrix, free_temp_matrix, & !!! 2025.02.03 nakata + matSFcoeff_tran, mult, sCaTr_sSs_aSs !!! 2025.02.03 nakata + use matrix_data, only: Hrange, Srange, aHa_range, aSs_range !!! 2025.02.03 nakata use primary_module, only: bundle - use species_module, only: species, nsf_species, species_label, n_species + use species_module, only: species, nsf_species, natomf_species, species_label, n_species !!! 2025.02.03 nakata use memory_module, only: type_dbl, type_int, type_cplx, & reg_alloc_mem, reg_dealloc_mem use energy, only: entropy @@ -529,13 +531,14 @@ subroutine FindEvals(electrons) bandE_total, coeff, setA, setB, Eband real(double), dimension(nspin) :: locc, bandE, entropy_local real(double), external :: dlamch - complex(double_cplx), dimension(:,:), allocatable :: scaledEig + complex(double_cplx), dimension(:,:), allocatable :: scaledEig, expH_atomf !!! 2025.02.03 nakata complex(double_cplx), dimension(:,:,:), allocatable :: expH complex(double_cplx) :: c_n_alpha2, c_n_setA2, c_n_setB2 - integer :: info, stat, il, iu, i, j, m, mz, prim_size, ng, wf_no, & + integer :: info, stat, il, iu, i, j, m, mz, prim_size, prim_size_atomf, ng, wf_no, & !!! 2025.02.03 nakata kp, spin, spin_SF, iacc, iprim, l, band, cdft_group, atom_fns_K, & n_band_min, n_band_max integer :: iatom_spec, Nangmom ! number of orbital angular momentum to be dumped (ex. (s,p,d)=3) + integer :: matStmp !!! 2025.02.03 nakata logical :: flag_keepexcite @@ -548,6 +551,12 @@ subroutine FindEvals(electrons) do i = 1, bundle%n_prim prim_size = prim_size + nsf_species(bundle%species(i)) end do +!!! 2025.02.03 nakata + prim_size_atomf = 0 + do i = 1, bundle%n_prim + prim_size_atomf = prim_size_atomf + natomf_species(bundle%species(i)) + end do +!!! nakata end ! Initialise - start BLACS, sort out matrices, allocate memory call initDiag @@ -597,13 +606,25 @@ subroutine FindEvals(electrons) if (stat /= 0) & call cq_abort('FindEvals: failed to alloc expH', stat) call reg_alloc_mem(area_DM, matrix_size * prim_size * nspin, type_cplx) - if(wf_self_con .and. flag_write_projected_DOS) then - allocate(scaledEig(matrix_size,prim_size), STAT=stat) - if (stat /= 0) & - call cq_abort('FindEvals: failed to alloc scaledEig', stat) - call reg_alloc_mem(area_DM, matrix_size * prim_size * nspin, type_cplx) - scaledEig = zero +!!! 2025.02.03 nakata + if(wf_self_con) then + if (flag_write_projected_DOS) then + allocate(scaledEig(matrix_size,prim_size_atomf), STAT=stat) + if (stat /= 0) & + call cq_abort('FindEvals: failed to alloc scaledEig', stat) + call reg_alloc_mem(area_DM, matrix_size * prim_size_atomf * nspin, type_cplx) + scaledEig = zero + end if + if(atomf.ne.sf) then + ! expH_atomf: eigenvectors in PAO basis + allocate(expH_atomf(matrix_size,prim_size_atomf), STAT=stat) + if (stat /= 0) & + call cq_abort('FindEvals: failed to alloc expH_atomf', stat) + call reg_alloc_mem(area_DM, matrix_size * prim_size_atomf, type_cplx) + expH_atomf = zero + end if end if +!!! 2025.02.03 nakata end !call gcopy(Efermi) !call gcopy(occ,matrix_size,nkp) ! else @@ -797,6 +818,7 @@ subroutine FindEvals(electrons) ! Second diagonalisation - get eigenvectors and build K entropy = zero flag_pDOS_buildK = .false. + flag_WFatomf_buildK = .false. spin_SF = 1 pol_S_size = maxval(ne_spin_in_cell) ! Size for polarisation matrix do spin = 1, nspin @@ -823,26 +845,57 @@ subroutine FindEvals(electrons) Hrange, matK(spin) ! Output wavefunction coefficients if(wf_self_con .and. (flag_out_wf .or. flag_write_projected_DOS)) then - if(i==1) then - call write_wavefn_coeffs(evals(:,kp,spin),expH(:,:,spin),spin,firstcall=1) - else - call write_wavefn_coeffs(evals(:,kp,spin),expH(:,:,spin),spin) - end if + if (atomf.ne.sf) expH_atomf = zero if(flag_write_projected_DOS) then scaledEig = zero flag_pDOS_buildK = .true. - call buildK(Hrange, matK(spin), occ(:,kp,spin), & - kk(:,kp), wtk(kp), expH(:,:,spin),scaledEig,matS(spin_SF)) + if (atomf.ne.sf) then + flag_WFatomf_buildK = .true. + matStmp = allocate_temp_matrix(aSs_range,0,atomf,sf) + ! call matrix_product(matSatomf, matSFcoeff_tran(spin_SF), matStmp, mult(aSa_sCaTr_aSs)) + call matrix_product(matSFcoeff_tran(spin_SF), matS(spin_SF), matStmp, mult(sCaTr_sSs_aSs)) + call buildK(Hrange, matK(spin), occ(:,kp,spin), & + kk(:,kp), wtk(kp), expH(:,:,spin), & + scaledEig, matStmp, & + Eig_atomf=expH_atomf, matSFcoeffTran=matSFcoeff_tran(spin_SF)) + flag_WFatomf_buildK = .false. + else + call buildK(Hrange, matK(spin), occ(:,kp,spin), & + kk(:,kp), wtk(kp), expH(:,:,spin), & + scaledEig, matS(spin_SF)) + endif flag_pDOS_buildK = .false. if(i==1) then call write_wavefn_coeffs(evals(:,kp,spin),scaledEig,spin,tag="Sij",firstcall=1) else call write_wavefn_coeffs(evals(:,kp,spin),scaledEig,spin,tag="Sij") end if + else if (atomf.ne.sf) then + write(io_lun,*) "test0 before calling buildK for WFs" ! 2025.02.03 nakata + flag_WFatomf_buildK = .true. + call buildK(Hrange, matK(spin), occ(:,kp,spin), & + kk(:,kp), wtk(kp), expH(:,:,spin), & + Eig_atomf=expH_atomf, matSFcoeffTran=matSFcoeff_tran(spin_SF)) + flag_WFatomf_buildK = .false. + write(io_lun,*) "test0 after calling buildK for WFs" ! 2025.02.03 nakata else call buildK(Hrange, matK(spin), occ(:,kp,spin), & kk(:,kp), wtk(kp), expH(:,:,spin)) end if + ! write WFs + if (atomf.ne.sf) then ! MSSFs + if(i==1) then + call write_wavefn_coeffs(evals(:,kp,spin),expH_atomf,spin,firstcall=1) + else + call write_wavefn_coeffs(evals(:,kp,spin),expH_atomf,spin) + end if + else ! primitive PAOs + if(i==1) then + call write_wavefn_coeffs(evals(:,kp,spin),expH(:,:,spin),spin,firstcall=1) + else + call write_wavefn_coeffs(evals(:,kp,spin),expH(:,:,spin),spin) + end if + endif else if(flag_do_pol_calc) then ! Set up polarisation calculation for this spin channel @@ -895,8 +948,10 @@ subroutine FindEvals(electrons) end do ! j = 1, matrix_size ! Now build data_M12_ij (=-\sum_n eps^n c^n_i c^n_j - ! hence scaling occs by eps allows reuse of buildK) + write(io_lun,*) "test 100.0" !!! 2025.02.03 nakata call buildK(Srange, matM12(spin), occ(:,kp,spin), & kk(:,kp), wtk(kp), expH(:,:,spin)) + write(io_lun,*) "test 100.1" !!! 2025.02.03 nakata end if ! End if (i <= N_kpoints_in_pg(ng)) then end do ! End do ng = 1, proc_groups end do ! End do i = 1, nkpoints_max @@ -969,6 +1024,13 @@ subroutine FindEvals(electrons) if (stat /= 0) call cq_abort('FindEvals: failed to deallocacte scaledEig', stat) call reg_dealloc_mem(area_DM, matrix_size * prim_size, type_cplx) end if +!!! 2025.02.03 nakata + if(wf_self_con .and. (atomf.ne.sf)) then + deallocate(expH_atomf, STAT=stat) + if (stat /= 0) call cq_abort('FindEvals: failed to deallocacte expH_atomf', stat) + call reg_dealloc_mem(area_DM, matrix_size * prim_size_atomf, type_cplx) + end if +!!! 2025.02.03 nakata end ! global call endDiag min_layer = min_layer + 1 @@ -3257,7 +3319,7 @@ end function MP_entropy !! Introduce flag_pol_build_S to select polarisation calculation !! SOURCE !! - subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) + subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, Eig_atomf, matSFcoeffTran) !use maxima_module, only: mx_nponn, mx_at_prim use numbers @@ -3272,13 +3334,14 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) use global_module, only: numprocs, iprint_DM, id_glob, & ni_in_cell, x_atom_cell, y_atom_cell, & z_atom_cell, max_wf, min_layer, flag_do_pol_calc, mat_polX_re, mat_polX_im, & - i_pol_dir_st, i_pol_dir_end, wf_self_con, flag_write_projected_DOS, ne_spin_in_cell + i_pol_dir_st, i_pol_dir_end, wf_self_con, flag_write_projected_DOS, ne_spin_in_cell, & + sf, atomf !!! 2025.02.03 nakata use mpi use GenBlas, only: dot use GenComms, only: myid use mult_module, only: store_matrix_value_pos, matrix_pos, matK, return_matrix_value_pos use matrix_data, only: mat, halo - use species_module, only: nsf_species + use species_module, only: nsf_species, natomf_species !!! 2025.02.03 nakata implicit none @@ -3291,6 +3354,11 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) ! For pDOS complex(double_cplx), optional, dimension(:,:) :: overlapEig integer, optional :: matSij +!!! 2025.02.03 nakata + ! For WFs with MSSFs + complex(double_cplx), optional, dimension(:,:) :: Eig_atomf + integer, optional :: matSFcoeffTran +!!! 2025.02.03 nakata end ! Local variables type(Krecv_data), dimension(:), allocatable :: recv_info @@ -3303,14 +3371,14 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) integer :: len, send_size, recv_size, send_proc, recv_proc, nsf1, & sendtag, recvtag integer :: req1, req2, ierr, atom, inter, prim, wheremat, row_sup,& - col_sup + col_sup, col_atomf !!! 2025.02.03 nakata integer, dimension(:,:), allocatable :: ints, atom_list, & send_prim, send_info, send_orbs, send_off integer, dimension(:), allocatable :: current_loc_atoms, & LocalAtom, num_send, norb_send, send_FSC, recv_to_FSC, & - mapchunk, prim_orbs + mapchunk, prim_orbs, prim_orbs_atomf !!! 2025.02.03 nakata integer, dimension(MPI_STATUS_SIZE) :: mpi_stat - real(double) :: phase, rfac, ifac, rcc, icc, rsum, exp_X_value_real, exp_X_value_imag, Siajb + real(double) :: phase, rfac, ifac, rcc, icc, rsum, exp_X_value_real, exp_X_value_imag, Siajb, SFcoeffTran_iajb !!! 2025.02.03 nakata complex(double_cplx) :: zsum, exp_X_value complex(double_cplx), dimension(:,:), allocatable :: RecvBuffer, & SendBuffer @@ -3411,7 +3479,7 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) write(io_lun,fmt='(10x,i6,a,3i6)') myid,' Maxima: ',maxloc, maxint, maxsend ! Allocate recv_info allocate(send_info(numprocs,maxsend),send_orbs(numprocs,maxsend),send_off(numprocs,maxsend), & - prim_orbs(bundle%mx_iprim),STAT=stat) + prim_orbs(bundle%mx_iprim),prim_orbs_atomf(bundle%mx_iprim),STAT=stat) !!! 2025.02.03 nakata if(stat/=0) call cq_abort('buildK: Error allocating send_info !',stat) send_info = 0 send_orbs = 0 @@ -3422,6 +3490,14 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) prim_orbs(j) = orb_count orb_count = orb_count + nsf_species(bundle%species(j)) end do +!!! 2025.02.03 nakata + prim_orbs_atomf = 0 + orb_count = 0 + do j=1,bundle%n_prim + prim_orbs_atomf(j) = orb_count + orb_count = orb_count + natomf_species(bundle%species(j)) + end do +!!! 2025.02.03 nakata end allocate(recv_info(numprocs),STAT=stat) if(stat/=0) call cq_abort('buildK: Error allocating recv_info !',stat) do i=1,numprocs @@ -3504,8 +3580,12 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) if(myid==0.AND.iprint_DM + min_layer>=4) write(io_lun,fmt='(10x,a,f8.4)') 'Occ is ',occs(i) end if end do +!!! 2025.02.03 nakata if(wf_self_con .and. flag_write_projected_DOS) len = matrix_size + if(wf_self_con .and. flag_WFatomf_buildK) len = matrix_size +!!! 2025.02.03 nakata end len_occ = len + write(io_lun,*) "matrix_size, len, len_occ=", matrix_size, len, len_occ !!! 2025.02.03 nakata if(iprint_DM+min_layer>3.AND.myid==0) & write(io_lun,fmt='(10x,a,2i6)') 'buildK: Stage three len:',len, matA if(flag_do_pol_calc.AND.flag_pol_buildS) then @@ -3588,7 +3668,10 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) end do if(iprint_DM + min_layer>=4.AND.myid==0) & write(io_lun,fmt='(10x,a)') 'filling buffer' - if(.not.(wf_self_con .and. flag_write_projected_DOS)) then +!!! 2025.02.03 nakata +! if(.not.(wf_self_con .and. flag_write_projected_DOS)) then + if(.not.(flag_pDOS_buildK .or. flag_WFatomf_buildK)) then +!!! 2025.02.03 nakata end do j=1,len_occ ! This is a loop over eigenstates RecvBuffer(j,1:recv_info(recv_proc+1)%orbs) = RecvBuffer(j,1:recv_info(recv_proc+1)%orbs)*occ_correction*occs(j) end do @@ -3639,20 +3722,52 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) polSloc(1:pol_S_size,1:pol_S_size,dir),pol_S_size) end do end if - ! projected DOS - if(flag_pDOS_buildK .and. wf_self_con .and. flag_write_projected_DOS) then - whereMat = matrix_pos(matSij, prim, jatom, col_sup, row_sup) - Siajb = return_matrix_value_pos(matSij,whereMat) - ! 1:len_occ gives bands; we want c_jb^n * S_iajb - overlapEig(1:len_occ,prim_orbs(prim)+col_sup) = & - overlapEig(1:len_occ,prim_orbs(prim)+col_sup) + & - Siajb*RecvBuffer(1:len_occ,orb_count+row_sup)*cmplx(rfac,ifac,double_cplx) - end if +!!! 2025.02.03 nakata pDOS-MSSF +! ! projected DOS +! if(flag_pDOS_buildK .and. wf_self_con .and. flag_write_projected_DOS) then +! whereMat = matrix_pos(matSij, prim, jatom, col_sup, row_sup) +! Siajb = return_matrix_value_pos(matSij,whereMat) +! ! 1:len_occ gives bands; we want c_jb^n * S_iajb +! overlapEig(1:len_occ,prim_orbs(prim)+col_sup) = & +! overlapEig(1:len_occ,prim_orbs(prim)+col_sup) + & +! Siajb*RecvBuffer(1:len_occ,orb_count+row_sup)*cmplx(rfac,ifac,double_cplx) +! end if +!!! 2025.02.03 nakata end ! Build K or M3 whereMat = matrix_pos(matA, prim, jatom, col_sup, row_sup) zsum = dot(len_occ,localEig(1:len_occ,prim_orbs(prim)+col_sup),1,RecvBuffer(1:len_occ,orb_count+row_sup),1) call store_matrix_value_pos(matA,whereMat,real(zsum*cmplx(rfac,ifac,double_cplx),double)) end do ! col_sup=nsf +!!! 2025.02.03 nakata pDOS-MSSF + ! projected DOS in atomf basis + if(wf_self_con) then + do col_atomf = 1,natomf_species(bundle%species(prim)) + ! c (WF coefficients) in PAO basis + if (flag_WFatomf_buildK .and. (atomf.ne.sf)) then + write(io_lun,*) 'sub:buildK test1' + write(io_lun,'(A,4I3,F20.15)') 'prim, jatom, col_atomf, row_sup', prim, jatom, col_atomf, row_sup + whereMat = matrix_pos(matSFcoeffTran, prim, jatom, col_atomf, row_sup) + write(io_lun,*) 'sub:buildK test1.1' + SFcoeffTran_iajb = return_matrix_value_pos(matSFcoeffTran,whereMat) + write(io_lun,'(A,4I3,F20.15)') 'prim, jatom, col_atomf, row_sup, SFcoeffTran_iajb', prim, jatom, col_atomf, row_sup, SFcoeffTran_iajb + ! 1:len_occ gives bands; we want c_ia(pao)^n = c_jb^n(sf) * SFcoeffTran_iajb(pao,sf) + Eig_atomf(1:len_occ,prim_orbs_atomf(prim)+col_atomf) = & + Eig_atomf(1:len_occ,prim_orbs_atomf(prim)+col_atomf) + & + SFcoeffTran_iajb*RecvBuffer(1:len_occ,orb_count+row_sup)*cmplx(rfac,ifac,double_cplx) + write(io_lun,*) 'sub:buildK test1.2' + endif + ! c*S for pDOS + if (flag_pDOS_buildK .and. flag_write_projected_DOS) then + whereMat = matrix_pos(matSij, prim, jatom, col_atomf, row_sup) + Siajb = return_matrix_value_pos(matSij,whereMat) + ! 1:len_occ gives bands; we want d_ia^n = c_jb^n(sf) * S_iajb(pao,sf) + overlapEig(1:len_occ,prim_orbs_atomf(prim)+col_atomf) = & + overlapEig(1:len_occ,prim_orbs_atomf(prim)+col_atomf) + & + Siajb*RecvBuffer(1:len_occ,orb_count+row_sup)*cmplx(rfac,ifac,double_cplx) + endif + end do ! col_atomf=natomf + end if +!!! 2025.02.03 nakata end end do ! row_sup=nsf end do ! inter=recv_info%ints ! Careful - we only want to increment after ALL interactions done @@ -3685,7 +3800,7 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij) recv_info(i)%dx,recv_info(i)%dy,recv_info(i)%dz,STAT=stat) if(stat/=0) call cq_abort('buildK: Error deallocating recvinfo !',i,stat) end do - deallocate(prim_orbs,send_off,send_orbs,send_info,STAT=stat) + deallocate(prim_orbs_atomf,prim_orbs,send_off,send_orbs,send_info,STAT=stat) !!! 2025.02.03 nakata if(stat/=0) call cq_abort('buildK: Error allocating send_info !',stat) deallocate(recv_info,STAT=stat) if(stat/=0) call cq_abort('buildK: Error allocating recv_info !',stat) @@ -4003,6 +4118,10 @@ end subroutine buildK !! Added support for writing out specific bands !! 2023/05/10 08:24 dave !! Reworked for post-processing requirements and added binary output option + !! 2025/02/03 17:30 nakata + !! Changed nsf_species to natomf_species + !! to write out WFs in PAO-basis rather than SF-basis. + !! WFs with PAO-basis is treatable in PostProcess for both PAOs and MSSFs. !! SOURCE !! subroutine write_wavefn_coeffs(eval, evec, spin, tag, firstcall) @@ -4065,10 +4184,10 @@ subroutine write_wavefn_coeffs(eval, evec, spin, tag, firstcall) acc = 0 do atom=1,bundle%n_prim write(lun) bundle%ig_prim(atom) - do isf1 = 1,nsf_species(bundle%species(atom)) + do isf1 = 1,natomf_species(bundle%species(atom)) write(lun) evec(wf_no,acc+isf1) end do - acc = acc + nsf_species(bundle%species(atom)) + acc = acc + natomf_species(bundle%species(atom)) end do end do ! iwf else @@ -4078,10 +4197,10 @@ subroutine write_wavefn_coeffs(eval, evec, spin, tag, firstcall) acc = 0 do atom=1,bundle%n_prim write(lun) bundle%ig_prim(atom) - do isf1 = 1,nsf_species(bundle%species(atom)) + do isf1 = 1,natomf_species(bundle%species(atom)) write(lun) evec(iwf,acc+isf1) end do - acc = acc + nsf_species(bundle%species(atom)) + acc = acc + natomf_species(bundle%species(atom)) end do end if end do ! iwf @@ -4100,10 +4219,10 @@ subroutine write_wavefn_coeffs(eval, evec, spin, tag, firstcall) acc = 0 do atom=1,bundle%n_prim write(lun,*) bundle%ig_prim(atom) - do isf1 = 1,nsf_species(bundle%species(atom)) + do isf1 = 1,natomf_species(bundle%species(atom)) write(lun,*) evec(wf_no,acc+isf1) end do - acc = acc + nsf_species(bundle%species(atom)) + acc = acc + natomf_species(bundle%species(atom)) end do end do ! iwf else @@ -4113,10 +4232,10 @@ subroutine write_wavefn_coeffs(eval, evec, spin, tag, firstcall) acc = 0 do atom=1,bundle%n_prim write(lun,*) bundle%ig_prim(atom) - do isf1 = 1,nsf_species(bundle%species(atom)) + do isf1 = 1,natomf_species(bundle%species(atom)) write(lun,*) evec(iwf,acc+isf1) end do - acc = acc + nsf_species(bundle%species(atom)) + acc = acc + natomf_species(bundle%species(atom)) end do end if end do ! iwf diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index dc3db8d17..a339a1d28 100644 --- a/src/initial_read_module.f90 +++ b/src/initial_read_module.f90 @@ -1765,6 +1765,7 @@ subroutine read_input(start, start_L, titles, vary_mu,& if(flag_diagonalisation) then flag_write_projected_DOS = fdf_boolean('IO.write_proj_DOS',.false.) if(flag_write_projected_DOS) then + flag_out_wf = .true. !!! 2025.02.03 nakata E_wf_min = fdf_double('IO.min_wf_E',-BIG) E_wf_max = fdf_double('IO.max_wf_E',BIG) end if From 1b991b7116fc0d9a28e323419b6e3cdd0450f203 Mon Sep 17 00:00:00 2001 From: ayakon Date: Thu, 6 Feb 2025 20:34:44 +0900 Subject: [PATCH 2/7] pDOS with MSSFs (in PostProcessing) --- tools/PostProcessing/pseudo_tm_info.f90 | 7 +++++++ tools/PostProcessing/read_module.f90 | 3 ++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/tools/PostProcessing/pseudo_tm_info.f90 b/tools/PostProcessing/pseudo_tm_info.f90 index 5ff6bea62..1dee95a7a 100644 --- a/tools/PostProcessing/pseudo_tm_info.f90 +++ b/tools/PostProcessing/pseudo_tm_info.f90 @@ -151,6 +151,9 @@ module pseudo_tm_info !! Add IF with pseudo_type for calling init_rad when ghost atoms are used !! 2019/12/09 20:16 dave !! Changes to read valence charge from pseudopotential file and test spin-polarised charge + !! 2025/02/04 17:00 nakata + !! nsf_species is set to npao_species anyway, even when MSSFs are used, + !! because evec_coeff was changed to be in the (pao, wf) format from the (sf, wf) format. !! SOURCE !! subroutine setup_pseudo_info @@ -206,8 +209,12 @@ subroutine setup_pseudo_info call read_ion_ascii_tmp(pseudo(ispecies),pao(ispecies)) npao_species(ispecies) = pao(ispecies)%count +!!! 2025.02.03 nakata ! Set NSF if not set by user if(nsf_species(ispecies)==0) nsf_species(ispecies) = pao(ispecies)%count + ! Set NSF to NPAO anyway (even if NSF is set by user for MSSFs) + nsf_species(ispecies) = pao(ispecies)%count +!!! 2025.02.03 nakata end maxnsf = max(maxnsf,nsf_species(ispecies)) ! Find radius for atom functions do l=0,pao(ispecies)%greatest_angmom diff --git a/tools/PostProcessing/read_module.f90 b/tools/PostProcessing/read_module.f90 index 2e381d2e0..c0e03c5b7 100644 --- a/tools/PostProcessing/read_module.f90 +++ b/tools/PostProcessing/read_module.f90 @@ -118,7 +118,8 @@ subroutine read_input i_job = 2 else if(leqi(job,'ban')) then i_job = 3 - if(flag_Multisite) call cq_abort("Not yet compatible with multi-site support functions") +!!! nakata 2025.02.03 comment out +! if(flag_Multisite) call cq_abort("Not yet compatible with multi-site support functions") else if(leqi(job,'ter').or.leqi(job,'th')) then i_job = 4 if(flag_Multisite) call cq_abort("Not yet compatible with multi-site support functions") From 075a65c0904e7bc8784a1a5ca8cccb91195c2955 Mon Sep 17 00:00:00 2001 From: ayakon Date: Fri, 7 Feb 2025 00:26:52 +0900 Subject: [PATCH 3/7] pDOS with MSSFs: introducing aSs_in_sSs_range etc. --- src/DiagModule.f90 | 26 +++++++++++++++++++------- src/dimens_module.f90 | 15 ++++++++++++++- src/matrix_data_module.f90 | 14 ++++++++++++-- src/mult_module.f90 | 38 +++++++++++++++++++++++++++++++++++--- 4 files changed, 80 insertions(+), 13 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index ebfa45dce..94b1cae57 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -503,9 +503,9 @@ subroutine FindEvals(electrons) nkpoints_max, pgid, N_procs_in_pg, & N_kpoints_in_pg use mult_module, only: matH, matS, matK, matM12, SF_to_AtomF_transform, & - matrix_scale, matrix_product, matrix_product_trace, allocate_temp_matrix, free_temp_matrix, & !!! 2025.02.03 nakata + matrix_scale, matrix_product, matrix_product_trace, matrix_sum, allocate_temp_matrix, free_temp_matrix, & !!! 2025.02.03 nakata matSFcoeff_tran, mult, sCaTr_sSs_aSs !!! 2025.02.03 nakata - use matrix_data, only: Hrange, Srange, aHa_range, aSs_range !!! 2025.02.03 nakata + use matrix_data, only: Hrange, Srange, aHa_range, aSs_range, aSs_in_sSs_range !!! 2025.02.03 nakata use primary_module, only: bundle use species_module, only: species, nsf_species, natomf_species, species_label, n_species !!! 2025.02.03 nakata use memory_module, only: type_dbl, type_int, type_cplx, & @@ -538,7 +538,7 @@ subroutine FindEvals(electrons) kp, spin, spin_SF, iacc, iprim, l, band, cdft_group, atom_fns_K, & n_band_min, n_band_max integer :: iatom_spec, Nangmom ! number of orbital angular momentum to be dumped (ex. (s,p,d)=3) - integer :: matStmp !!! 2025.02.03 nakata + integer :: matStmp, matStmp0, matSFcoeffTran_tmp !!! 2025.02.03 nakata logical :: flag_keepexcite @@ -851,13 +851,21 @@ subroutine FindEvals(electrons) flag_pDOS_buildK = .true. if (atomf.ne.sf) then flag_WFatomf_buildK = .true. - matStmp = allocate_temp_matrix(aSs_range,0,atomf,sf) + matStmp0 = allocate_temp_matrix(aSs_range,0,atomf,sf) ! call matrix_product(matSatomf, matSFcoeff_tran(spin_SF), matStmp, mult(aSa_sCaTr_aSs)) - call matrix_product(matSFcoeff_tran(spin_SF), matS(spin_SF), matStmp, mult(sCaTr_sSs_aSs)) + call matrix_product(matSFcoeff_tran(spin_SF), matS(spin_SF), matStmp0, mult(sCaTr_sSs_aSs)) + ! change matSFcoeff_tran and matStmp from aSs_range to sSs_range + matStmp = allocate_temp_matrix(aSs_in_sSs_range,0,atomf,sf) + call matrix_sum(zero,matStmp,one,matStmp0) + matSFcoeffTran_tmp = allocate_temp_matrix(aSs_in_sSs_range,0,atomf,sf) + call matrix_sum(zero,matSFcoeffTran_tmp,one,matSFcoeff_tran(spin_SF)) call buildK(Hrange, matK(spin), occ(:,kp,spin), & kk(:,kp), wtk(kp), expH(:,:,spin), & scaledEig, matStmp, & - Eig_atomf=expH_atomf, matSFcoeffTran=matSFcoeff_tran(spin_SF)) + Eig_atomf=expH_atomf, matSFcoeffTran=matSFcoeffTran_tmp) + call free_temp_matrix(matSFcoeffTran_tmp) + call free_temp_matrix(matStmp) + call free_temp_matrix(matStmp0) flag_WFatomf_buildK = .false. else call buildK(Hrange, matK(spin), occ(:,kp,spin), & @@ -873,9 +881,13 @@ subroutine FindEvals(electrons) else if (atomf.ne.sf) then write(io_lun,*) "test0 before calling buildK for WFs" ! 2025.02.03 nakata flag_WFatomf_buildK = .true. + ! change matSFcoeff_tran from aSs_range to sSs_range + matSFcoeffTran_tmp = allocate_temp_matrix(aSs_in_sSs_range,0,atomf,sf) + call matrix_sum(zero,matSFcoeffTran_tmp,one,matSFcoeff_tran(spin_SF)) call buildK(Hrange, matK(spin), occ(:,kp,spin), & kk(:,kp), wtk(kp), expH(:,:,spin), & - Eig_atomf=expH_atomf, matSFcoeffTran=matSFcoeff_tran(spin_SF)) + Eig_atomf=expH_atomf, matSFcoeffTran=matSFcoeffTran_tmp) + call free_temp_matrix(matSFcoeffTran_tmp) flag_WFatomf_buildK = .false. write(io_lun,*) "test0 after calling buildK for WFs" ! 2025.02.03 nakata else diff --git a/src/dimens_module.f90 b/src/dimens_module.f90 index 467488d91..ad89da916 100644 --- a/src/dimens_module.f90 +++ b/src/dimens_module.f90 @@ -235,10 +235,19 @@ subroutine set_dimensions(inode, ionode,HNL_fac,non_local, n_species, non_local_ if(flag_neutral_atom_projector) then aNArange = 31 NAarange = 32 - mx_matrices_tmp = mx_matrices ! = 30 +!!! 2025.02.03 nakata +! mx_matrices_tmp = mx_matrices ! = 30 + mx_matrices_tmp = 32 +!!! 2025.02.03 nakata end else mx_matrices_tmp = 30 end if +!!! 2025.02.03 nakata + aSs_in_sSs_range = mx_matrices_tmp + 1 + sSa_in_sSs_range = mx_matrices_tmp + 2 + mx_matrices_tmp = mx_matrices_tmp + 2 + if (mx_matrices_tmp > mx_matrices) call cq_abort('ERROR : mx_matrices_tmp is larger than mx_matrices',mx_matrices_tmp) +!!! 2025.02.03 nakata end endif !n_my_grid_points = n_pts_in_block * n_blocks @@ -371,6 +380,8 @@ subroutine set_dimensions(inode, ionode,HNL_fac,non_local, n_species, non_local_ rcut(SFcoeffTr_range) = 0.001_double endif if (abs(r_LD)1) then do n=1,mx_matrices_tmp diff --git a/src/matrix_data_module.f90 b/src/matrix_data_module.f90 index aaa3f3507..a9a29e029 100644 --- a/src/matrix_data_module.f90 +++ b/src/matrix_data_module.f90 @@ -62,6 +62,8 @@ !! which are no longer used !! 2017/12/05 10:20 dave (with TM and NW (Mizuho)) !! Adding new matrix indices (aNA and NAa) for atom function - NA projectors +!! 2025/02/06 14:30 nakata +!! aSs_in_sSs_range, aSs_in_sSs_matind were added for pDOS with MSSFs !! SOURCE !! module matrix_data @@ -73,7 +75,10 @@ module matrix_data save ! This will need to change if the above parameters are changed - integer, parameter :: mx_matrices = 32 +!!! 2025.02.03 nakata +! integer, parameter :: mx_matrices = 32 + integer, parameter :: mx_matrices = 34 +!!! 2025.02.03 nakata end ! Store ALL indices in a large array type(matrix), allocatable, dimension(:,:), target :: mat @@ -86,7 +91,8 @@ module matrix_data SLSmatind, Tmatind, TTrmatind, TSmatind, THmatind, TLmatind, Xmatind, SXmatind integer, dimension(:), pointer :: aSa_matind, aHa_matind, STr_matind, HTr_matind, & aSs_matind, aHs_matind, sSa_matind, sHa_matind, & - SFcoeff_matind, SFcoeffTr_matind, LD_matind + SFcoeff_matind, SFcoeffTr_matind, LD_matind, & + aSs_in_sSs_matind, sSa_in_sSs_matind !!! 2025.02.03 nakata integer, dimension(:), pointer :: aNAmatind, NAamatind ! Parameters for the different matrix ranges @@ -126,6 +132,10 @@ module matrix_data ! Ranges for NA projectors set later also (dimens.module.f90) integer :: aNArange ! 31 integer :: NAarange ! 32 +!!! 2025.02.03 nakata + integer :: aSs_in_sSs_range ! 33 for S(atomf,sf) but with the range of Srange (= r_sf + r_sf, not r_atomf + r_sf) + integer :: sSa_in_sSs_range ! 34 for S(sf,atomf) but with the range of Srange (= r_sf + r_sf, not r_atomf + r_sf) +!!! 2025.02.03 nakata end integer :: max_range ! Indexes matrix with largest range diff --git a/src/mult_module.f90 b/src/mult_module.f90 index 784960228..e515b158a 100644 --- a/src/mult_module.f90 +++ b/src/mult_module.f90 @@ -100,6 +100,8 @@ !! Adding multiplications for NA projectors !! 2018/11/13 17:30 nakata !! Changed matS, matT, matTtran, matKE, matNL and matNA to be spin_SF dependent +!! 2025/02/06 14:00 nakata +!! Added aSs_in_sSs_trans, aSs_in_sSs_pairind, sSa_in_sSs_trans, sSa_in_sSs_pairind for pDOS with MSSFs !! SOURCE !! module mult_module @@ -176,15 +178,22 @@ module mult_module integer :: SFcoeff_trans ! 10 integer :: aNA_trans ! 11 integer :: aNAa_trans ! 12 + integer :: aSs_in_sSs_trans ! 13 ! 2025.02.04 nakata + integer :: sSa_in_sSs_trans ! 14 ! 2025.02.04 nakata - integer(integ), parameter :: mx_trans = 12 +!!! 2025.02.03 nakata +! integer(integ), parameter :: mx_trans = 12 + integer(integ), parameter :: mx_trans = 14 +!!! 2025.02.03 nakata end type(pair_data), allocatable, dimension(:,:) :: pairs integer, dimension(:), pointer :: Spairind, Lpairind, Tpairind, & APpairind, LSpairind, LHpairind, & LSLpairind, & aSs_pairind, aHs_pairind, SFcoeff_pairind, & - aNApairind, aNAapairind + aNApairind, aNAapairind, & + aSs_in_sSs_pairind , & !!! 2025.02.03 nakata + sSa_in_sSs_pairind !!! 2025.02.03 nakata type(matrix_trans), dimension(mx_matrices), target :: ltrans type(trans_remote) :: gtrans(mx_trans) @@ -354,6 +363,8 @@ subroutine immi(parts, prim, gcs, myid, partial) SFcoeff_trans = 10 aNA_trans = 11 aNAa_trans = 12 + aSs_in_sSs_trans = 13 !!! 2025.02.03 nakata + sSa_in_sSs_trans = 14 !!! 2025.02.03 nakata else aNA_NAa_aHa = 23 aHa_aNA_aNA = 24 @@ -478,6 +489,18 @@ subroutine immi(parts, prim, gcs, myid, partial) call matrix_ini(parts, prim, gcs, mat(1:prim%groups_on_node,SFcoeffTr_range), & SFcoeffTr_matind, rcut(SFcoeffTr_range), myid-1, & halo(SFcoeffTr_range), ltrans(SFcoeffTr_range)) +!!! 2025.02.03 nakata + mat(1:prim%groups_on_node,aSs_in_sSs_range)%sf1_type = atomf + mat(1:prim%groups_on_node,aSs_in_sSs_range)%sf2_type = sf + call matrix_ini(parts, prim, gcs, mat(1:prim%groups_on_node,aSs_in_sSs_range), & + aSs_in_sSs_matind, rcut(aSs_in_sSs_range), myid-1, & + halo(aSs_in_sSs_range), ltrans(aSs_in_sSs_range)) + mat(1:prim%groups_on_node,sSa_in_sSs_range)%sf1_type = sf + mat(1:prim%groups_on_node,sSa_in_sSs_range)%sf2_type = atomf + call matrix_ini(parts, prim, gcs, mat(1:prim%groups_on_node,sSa_in_sSs_range), & + sSa_in_sSs_matind, rcut(sSa_in_sSs_range), myid-1, & + halo(sSa_in_sSs_range), ltrans(sSa_in_sSs_range)) +!!! if (flag_LFD) then mat(1:prim%groups_on_node,LD_range)%sf1_type = atomf mat(1:prim%groups_on_node,LD_range)%sf2_type = atomf @@ -551,6 +574,14 @@ subroutine immi(parts, prim, gcs, myid, partial) myid-1, halo(aHa_range), halo(aHa_range), ltrans(aHa_range), & gtrans(aNAa_trans), pairs(:,aNAa_trans), aNAapairind) end if +!!! 2025.02.03 nakata + call trans_ini(parts, prim, gcs, mat(1:prim%groups_on_node,aSs_in_sSs_range), & + myid-1, halo(aSs_in_sSs_range), halo(sSa_in_sSs_range), ltrans(aSs_in_sSs_range), & + gtrans(aSs_in_sSs_trans), pairs(:, aSs_in_sSs_trans), aSs_in_sSs_pairind) + call trans_ini(parts, prim, gcs, mat(1:prim%groups_on_node,sSa_in_sSs_range), & + myid-1, halo(sSa_in_sSs_range), halo(aSs_in_sSs_range), ltrans(sSa_in_sSs_range), & + gtrans(sSa_in_sSs_trans), pairs(:, sSa_in_sSs_trans), sSa_in_sSs_pairind) +!!! 2025.02.03 nakata end ! NA projectors else if(flag_neutral_atom_projector) then call trans_ini(parts, prim, gcs, mat(1:prim%groups_on_node,aNArange), & @@ -1359,7 +1390,7 @@ subroutine fmmi(prim) deallocate(pairs) deallocate(Spairind, Lpairind, APpairind, LSpairind, LHpairind, & LSLpairind, Tpairind) - if (atomf.ne.sf) deallocate(aSs_pairind, aHs_pairind, SFcoeff_pairind) + if (atomf.ne.sf) deallocate(aSs_pairind, aHs_pairind, SFcoeff_pairind, aSs_in_sSs_pairind) if(flag_neutral_atom_projector) then deallocate(aNApairind) if(atomf.ne.sf) deallocate(aNAapairind) @@ -1397,6 +1428,7 @@ subroutine fmmi(prim) call end_ops(prim,SFcoeff_range,SFcoeff_matind,SFcoeff_trans) call end_ops(prim,SFcoeffTr_range,SFcoeffTr_matind) if (flag_LFD) call end_ops(prim,LD_range,LD_matind) + call end_ops(prim,aSs_in_sSs_range,aSs_in_sSs_matind,aSs_in_sSs_trans) endif if( flag_neutral_atom_projector ) then call end_ops(prim,aNArange, aNAmatind,aNA_trans) From d961c169c2c0b6ecd95304a01f56c77a381cb1fe Mon Sep 17 00:00:00 2001 From: David Bowler Date: Thu, 6 Feb 2025 15:57:20 +0000 Subject: [PATCH 4/7] Tidy up unnecessary variables --- src/dimens_module.f90 | 7 ++----- src/matrix_data_module.f90 | 3 +-- src/mult_module.f90 | 25 +++---------------------- 3 files changed, 6 insertions(+), 29 deletions(-) diff --git a/src/dimens_module.f90 b/src/dimens_module.f90 index ad89da916..a296ea04a 100644 --- a/src/dimens_module.f90 +++ b/src/dimens_module.f90 @@ -244,8 +244,7 @@ subroutine set_dimensions(inode, ionode,HNL_fac,non_local, n_species, non_local_ end if !!! 2025.02.03 nakata aSs_in_sSs_range = mx_matrices_tmp + 1 - sSa_in_sSs_range = mx_matrices_tmp + 2 - mx_matrices_tmp = mx_matrices_tmp + 2 + mx_matrices_tmp = mx_matrices_tmp + 1 if (mx_matrices_tmp > mx_matrices) call cq_abort('ERROR : mx_matrices_tmp is larger than mx_matrices',mx_matrices_tmp) !!! 2025.02.03 nakata end endif @@ -380,8 +379,7 @@ subroutine set_dimensions(inode, ionode,HNL_fac,non_local, n_species, non_local_ rcut(SFcoeffTr_range) = 0.001_double endif if (abs(r_LD)1) then do n=1,mx_matrices_tmp diff --git a/src/matrix_data_module.f90 b/src/matrix_data_module.f90 index a9a29e029..ec9b67665 100644 --- a/src/matrix_data_module.f90 +++ b/src/matrix_data_module.f90 @@ -92,7 +92,7 @@ module matrix_data integer, dimension(:), pointer :: aSa_matind, aHa_matind, STr_matind, HTr_matind, & aSs_matind, aHs_matind, sSa_matind, sHa_matind, & SFcoeff_matind, SFcoeffTr_matind, LD_matind, & - aSs_in_sSs_matind, sSa_in_sSs_matind !!! 2025.02.03 nakata + aSs_in_sSs_matind !!! 2025.02.03 nakata integer, dimension(:), pointer :: aNAmatind, NAamatind ! Parameters for the different matrix ranges @@ -134,7 +134,6 @@ module matrix_data integer :: NAarange ! 32 !!! 2025.02.03 nakata integer :: aSs_in_sSs_range ! 33 for S(atomf,sf) but with the range of Srange (= r_sf + r_sf, not r_atomf + r_sf) - integer :: sSa_in_sSs_range ! 34 for S(sf,atomf) but with the range of Srange (= r_sf + r_sf, not r_atomf + r_sf) !!! 2025.02.03 nakata end integer :: max_range ! Indexes matrix with largest range diff --git a/src/mult_module.f90 b/src/mult_module.f90 index e515b158a..335174561 100644 --- a/src/mult_module.f90 +++ b/src/mult_module.f90 @@ -178,8 +178,6 @@ module mult_module integer :: SFcoeff_trans ! 10 integer :: aNA_trans ! 11 integer :: aNAa_trans ! 12 - integer :: aSs_in_sSs_trans ! 13 ! 2025.02.04 nakata - integer :: sSa_in_sSs_trans ! 14 ! 2025.02.04 nakata !!! 2025.02.03 nakata ! integer(integ), parameter :: mx_trans = 12 @@ -191,9 +189,7 @@ module mult_module APpairind, LSpairind, LHpairind, & LSLpairind, & aSs_pairind, aHs_pairind, SFcoeff_pairind, & - aNApairind, aNAapairind, & - aSs_in_sSs_pairind , & !!! 2025.02.03 nakata - sSa_in_sSs_pairind !!! 2025.02.03 nakata + aNApairind, aNAapairind type(matrix_trans), dimension(mx_matrices), target :: ltrans type(trans_remote) :: gtrans(mx_trans) @@ -363,8 +359,6 @@ subroutine immi(parts, prim, gcs, myid, partial) SFcoeff_trans = 10 aNA_trans = 11 aNAa_trans = 12 - aSs_in_sSs_trans = 13 !!! 2025.02.03 nakata - sSa_in_sSs_trans = 14 !!! 2025.02.03 nakata else aNA_NAa_aHa = 23 aHa_aNA_aNA = 24 @@ -495,11 +489,6 @@ subroutine immi(parts, prim, gcs, myid, partial) call matrix_ini(parts, prim, gcs, mat(1:prim%groups_on_node,aSs_in_sSs_range), & aSs_in_sSs_matind, rcut(aSs_in_sSs_range), myid-1, & halo(aSs_in_sSs_range), ltrans(aSs_in_sSs_range)) - mat(1:prim%groups_on_node,sSa_in_sSs_range)%sf1_type = sf - mat(1:prim%groups_on_node,sSa_in_sSs_range)%sf2_type = atomf - call matrix_ini(parts, prim, gcs, mat(1:prim%groups_on_node,sSa_in_sSs_range), & - sSa_in_sSs_matind, rcut(sSa_in_sSs_range), myid-1, & - halo(sSa_in_sSs_range), ltrans(sSa_in_sSs_range)) !!! if (flag_LFD) then mat(1:prim%groups_on_node,LD_range)%sf1_type = atomf @@ -574,14 +563,6 @@ subroutine immi(parts, prim, gcs, myid, partial) myid-1, halo(aHa_range), halo(aHa_range), ltrans(aHa_range), & gtrans(aNAa_trans), pairs(:,aNAa_trans), aNAapairind) end if -!!! 2025.02.03 nakata - call trans_ini(parts, prim, gcs, mat(1:prim%groups_on_node,aSs_in_sSs_range), & - myid-1, halo(aSs_in_sSs_range), halo(sSa_in_sSs_range), ltrans(aSs_in_sSs_range), & - gtrans(aSs_in_sSs_trans), pairs(:, aSs_in_sSs_trans), aSs_in_sSs_pairind) - call trans_ini(parts, prim, gcs, mat(1:prim%groups_on_node,sSa_in_sSs_range), & - myid-1, halo(sSa_in_sSs_range), halo(aSs_in_sSs_range), ltrans(sSa_in_sSs_range), & - gtrans(sSa_in_sSs_trans), pairs(:, sSa_in_sSs_trans), sSa_in_sSs_pairind) -!!! 2025.02.03 nakata end ! NA projectors else if(flag_neutral_atom_projector) then call trans_ini(parts, prim, gcs, mat(1:prim%groups_on_node,aNArange), & @@ -1390,7 +1371,7 @@ subroutine fmmi(prim) deallocate(pairs) deallocate(Spairind, Lpairind, APpairind, LSpairind, LHpairind, & LSLpairind, Tpairind) - if (atomf.ne.sf) deallocate(aSs_pairind, aHs_pairind, SFcoeff_pairind, aSs_in_sSs_pairind) + if (atomf.ne.sf) deallocate(aSs_pairind, aHs_pairind, SFcoeff_pairind) if(flag_neutral_atom_projector) then deallocate(aNApairind) if(atomf.ne.sf) deallocate(aNAapairind) @@ -1428,7 +1409,7 @@ subroutine fmmi(prim) call end_ops(prim,SFcoeff_range,SFcoeff_matind,SFcoeff_trans) call end_ops(prim,SFcoeffTr_range,SFcoeffTr_matind) if (flag_LFD) call end_ops(prim,LD_range,LD_matind) - call end_ops(prim,aSs_in_sSs_range,aSs_in_sSs_matind,aSs_in_sSs_trans) + call end_ops(prim,aSs_in_sSs_range,aSs_in_sSs_matind) endif if( flag_neutral_atom_projector ) then call end_ops(prim,aNArange, aNAmatind,aNA_trans) From 078f82260d4afd522ae1b94130c765a47841be9b Mon Sep 17 00:00:00 2001 From: David Bowler Date: Wed, 12 Feb 2025 09:46:31 +0000 Subject: [PATCH 5/7] Change to overlap matrix used in pDOS for MSSF The PAO-SF overlap matrix is now generated from the PAO-PAO matrix --- src/DiagModule.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index 94b1cae57..056534d4c 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -504,7 +504,7 @@ subroutine FindEvals(electrons) N_kpoints_in_pg use mult_module, only: matH, matS, matK, matM12, SF_to_AtomF_transform, & matrix_scale, matrix_product, matrix_product_trace, matrix_sum, allocate_temp_matrix, free_temp_matrix, & !!! 2025.02.03 nakata - matSFcoeff_tran, mult, sCaTr_sSs_aSs !!! 2025.02.03 nakata + matSFcoeff_tran, mult, sCaTr_sSs_aSs, aSa_sCaTr_aSs, matSatomF !!! 2025.02.03 nakata use matrix_data, only: Hrange, Srange, aHa_range, aSs_range, aSs_in_sSs_range !!! 2025.02.03 nakata use primary_module, only: bundle use species_module, only: species, nsf_species, natomf_species, species_label, n_species !!! 2025.02.03 nakata @@ -852,8 +852,8 @@ subroutine FindEvals(electrons) if (atomf.ne.sf) then flag_WFatomf_buildK = .true. matStmp0 = allocate_temp_matrix(aSs_range,0,atomf,sf) - ! call matrix_product(matSatomf, matSFcoeff_tran(spin_SF), matStmp, mult(aSa_sCaTr_aSs)) - call matrix_product(matSFcoeff_tran(spin_SF), matS(spin_SF), matStmp0, mult(sCaTr_sSs_aSs)) + call matrix_product(matSatomf, matSFcoeff_tran(spin_SF), matStmp0, mult(aSa_sCaTr_aSs)) + !call matrix_product(matSFcoeff_tran(spin_SF), matS(spin_SF), matStmp0, mult(sCaTr_sSs_aSs)) ! change matSFcoeff_tran and matStmp from aSs_range to sSs_range matStmp = allocate_temp_matrix(aSs_in_sSs_range,0,atomf,sf) call matrix_sum(zero,matStmp,one,matStmp0) From 097837f7779a986c756a3ca6fbf03965567ac7b6 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Wed, 12 Feb 2025 13:47:53 +0000 Subject: [PATCH 6/7] Tidying up code Removing debug write statements --- src/DiagModule.f90 | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index 056534d4c..7f7f86b3b 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -502,9 +502,9 @@ subroutine FindEvals(electrons) block_size_c, pg_kpoints, proc_groups, & nkpoints_max, pgid, N_procs_in_pg, & N_kpoints_in_pg - use mult_module, only: matH, matS, matK, matM12, SF_to_AtomF_transform, & + use mult_module, only: matH, matS, matK, matM12, & matrix_scale, matrix_product, matrix_product_trace, matrix_sum, allocate_temp_matrix, free_temp_matrix, & !!! 2025.02.03 nakata - matSFcoeff_tran, mult, sCaTr_sSs_aSs, aSa_sCaTr_aSs, matSatomF !!! 2025.02.03 nakata + matSFcoeff_tran, mult, aSa_sCaTr_aSs, matSatomF !!! 2025.02.03 nakata use matrix_data, only: Hrange, Srange, aHa_range, aSs_range, aSs_in_sSs_range !!! 2025.02.03 nakata use primary_module, only: bundle use species_module, only: species, nsf_species, natomf_species, species_label, n_species !!! 2025.02.03 nakata @@ -853,7 +853,6 @@ subroutine FindEvals(electrons) flag_WFatomf_buildK = .true. matStmp0 = allocate_temp_matrix(aSs_range,0,atomf,sf) call matrix_product(matSatomf, matSFcoeff_tran(spin_SF), matStmp0, mult(aSa_sCaTr_aSs)) - !call matrix_product(matSFcoeff_tran(spin_SF), matS(spin_SF), matStmp0, mult(sCaTr_sSs_aSs)) ! change matSFcoeff_tran and matStmp from aSs_range to sSs_range matStmp = allocate_temp_matrix(aSs_in_sSs_range,0,atomf,sf) call matrix_sum(zero,matStmp,one,matStmp0) @@ -879,7 +878,6 @@ subroutine FindEvals(electrons) call write_wavefn_coeffs(evals(:,kp,spin),scaledEig,spin,tag="Sij") end if else if (atomf.ne.sf) then - write(io_lun,*) "test0 before calling buildK for WFs" ! 2025.02.03 nakata flag_WFatomf_buildK = .true. ! change matSFcoeff_tran from aSs_range to sSs_range matSFcoeffTran_tmp = allocate_temp_matrix(aSs_in_sSs_range,0,atomf,sf) @@ -889,7 +887,6 @@ subroutine FindEvals(electrons) Eig_atomf=expH_atomf, matSFcoeffTran=matSFcoeffTran_tmp) call free_temp_matrix(matSFcoeffTran_tmp) flag_WFatomf_buildK = .false. - write(io_lun,*) "test0 after calling buildK for WFs" ! 2025.02.03 nakata else call buildK(Hrange, matK(spin), occ(:,kp,spin), & kk(:,kp), wtk(kp), expH(:,:,spin)) @@ -960,10 +957,8 @@ subroutine FindEvals(electrons) end do ! j = 1, matrix_size ! Now build data_M12_ij (=-\sum_n eps^n c^n_i c^n_j - ! hence scaling occs by eps allows reuse of buildK) - write(io_lun,*) "test 100.0" !!! 2025.02.03 nakata call buildK(Srange, matM12(spin), occ(:,kp,spin), & kk(:,kp), wtk(kp), expH(:,:,spin)) - write(io_lun,*) "test 100.1" !!! 2025.02.03 nakata end if ! End if (i <= N_kpoints_in_pg(ng)) then end do ! End do ng = 1, proc_groups end do ! End do i = 1, nkpoints_max @@ -3597,7 +3592,6 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, if(wf_self_con .and. flag_WFatomf_buildK) len = matrix_size !!! 2025.02.03 nakata end len_occ = len - write(io_lun,*) "matrix_size, len, len_occ=", matrix_size, len, len_occ !!! 2025.02.03 nakata if(iprint_DM+min_layer>3.AND.myid==0) & write(io_lun,fmt='(10x,a,2i6)') 'buildK: Stage three len:',len, matA if(flag_do_pol_calc.AND.flag_pol_buildS) then @@ -3756,17 +3750,12 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, do col_atomf = 1,natomf_species(bundle%species(prim)) ! c (WF coefficients) in PAO basis if (flag_WFatomf_buildK .and. (atomf.ne.sf)) then - write(io_lun,*) 'sub:buildK test1' - write(io_lun,'(A,4I3,F20.15)') 'prim, jatom, col_atomf, row_sup', prim, jatom, col_atomf, row_sup whereMat = matrix_pos(matSFcoeffTran, prim, jatom, col_atomf, row_sup) - write(io_lun,*) 'sub:buildK test1.1' SFcoeffTran_iajb = return_matrix_value_pos(matSFcoeffTran,whereMat) - write(io_lun,'(A,4I3,F20.15)') 'prim, jatom, col_atomf, row_sup, SFcoeffTran_iajb', prim, jatom, col_atomf, row_sup, SFcoeffTran_iajb ! 1:len_occ gives bands; we want c_ia(pao)^n = c_jb^n(sf) * SFcoeffTran_iajb(pao,sf) Eig_atomf(1:len_occ,prim_orbs_atomf(prim)+col_atomf) = & Eig_atomf(1:len_occ,prim_orbs_atomf(prim)+col_atomf) + & SFcoeffTran_iajb*RecvBuffer(1:len_occ,orb_count+row_sup)*cmplx(rfac,ifac,double_cplx) - write(io_lun,*) 'sub:buildK test1.2' endif ! c*S for pDOS if (flag_pDOS_buildK .and. flag_write_projected_DOS) then From 5a7231a2920f23e2fa81336ea60be41e93cc2eb1 Mon Sep 17 00:00:00 2001 From: ayakon Date: Thu, 13 Feb 2025 20:22:40 +0900 Subject: [PATCH 7/7] Tidying up code Removing comments for debug --- src/DiagModule.f90 | 62 ++++++++----------------- src/dimens_module.f90 | 9 ++-- src/initial_read_module.f90 | 4 +- src/matrix_data_module.f90 | 12 ++--- src/mult_module.f90 | 11 ++--- tools/PostProcessing/pseudo_tm_info.f90 | 5 +- tools/PostProcessing/read_module.f90 | 4 -- 7 files changed, 33 insertions(+), 74 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index 7f7f86b3b..9055a95ba 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -478,6 +478,8 @@ module DiagModule !! 2023/03/15 08:37 dave !! Removed pDOS output (now done in post-processing) and added calls !! for wavefunction coefficients scaled by overlap matrix (for pDOS) + !! 2025/02/13 17:30 dave and nakata + !! Introduced expH_atomf and flag_WFatomf_buildK for pDOS with MSSFs !! SOURCE !! subroutine FindEvals(electrons) @@ -503,11 +505,11 @@ subroutine FindEvals(electrons) nkpoints_max, pgid, N_procs_in_pg, & N_kpoints_in_pg use mult_module, only: matH, matS, matK, matM12, & - matrix_scale, matrix_product, matrix_product_trace, matrix_sum, allocate_temp_matrix, free_temp_matrix, & !!! 2025.02.03 nakata - matSFcoeff_tran, mult, aSa_sCaTr_aSs, matSatomF !!! 2025.02.03 nakata - use matrix_data, only: Hrange, Srange, aHa_range, aSs_range, aSs_in_sSs_range !!! 2025.02.03 nakata + matrix_scale, matrix_product, matrix_product_trace, matrix_sum, allocate_temp_matrix, free_temp_matrix, & + matSFcoeff_tran, mult, aSa_sCaTr_aSs, matSatomf + use matrix_data, only: Hrange, Srange, aHa_range, aSs_range, aSs_in_sSs_range use primary_module, only: bundle - use species_module, only: species, nsf_species, natomf_species, species_label, n_species !!! 2025.02.03 nakata + use species_module, only: species, nsf_species, natomf_species, species_label, n_species use memory_module, only: type_dbl, type_int, type_cplx, & reg_alloc_mem, reg_dealloc_mem use energy, only: entropy @@ -531,14 +533,14 @@ subroutine FindEvals(electrons) bandE_total, coeff, setA, setB, Eband real(double), dimension(nspin) :: locc, bandE, entropy_local real(double), external :: dlamch - complex(double_cplx), dimension(:,:), allocatable :: scaledEig, expH_atomf !!! 2025.02.03 nakata + complex(double_cplx), dimension(:,:), allocatable :: scaledEig, expH_atomf complex(double_cplx), dimension(:,:,:), allocatable :: expH complex(double_cplx) :: c_n_alpha2, c_n_setA2, c_n_setB2 - integer :: info, stat, il, iu, i, j, m, mz, prim_size, prim_size_atomf, ng, wf_no, & !!! 2025.02.03 nakata + integer :: info, stat, il, iu, i, j, m, mz, prim_size, prim_size_atomf, ng, wf_no, & kp, spin, spin_SF, iacc, iprim, l, band, cdft_group, atom_fns_K, & n_band_min, n_band_max integer :: iatom_spec, Nangmom ! number of orbital angular momentum to be dumped (ex. (s,p,d)=3) - integer :: matStmp, matStmp0, matSFcoeffTran_tmp !!! 2025.02.03 nakata + integer :: matStmp, matStmp0, matSFcoeffTran_tmp logical :: flag_keepexcite @@ -551,12 +553,10 @@ subroutine FindEvals(electrons) do i = 1, bundle%n_prim prim_size = prim_size + nsf_species(bundle%species(i)) end do -!!! 2025.02.03 nakata prim_size_atomf = 0 do i = 1, bundle%n_prim prim_size_atomf = prim_size_atomf + natomf_species(bundle%species(i)) end do -!!! nakata end ! Initialise - start BLACS, sort out matrices, allocate memory call initDiag @@ -606,7 +606,6 @@ subroutine FindEvals(electrons) if (stat /= 0) & call cq_abort('FindEvals: failed to alloc expH', stat) call reg_alloc_mem(area_DM, matrix_size * prim_size * nspin, type_cplx) -!!! 2025.02.03 nakata if(wf_self_con) then if (flag_write_projected_DOS) then allocate(scaledEig(matrix_size,prim_size_atomf), STAT=stat) @@ -624,7 +623,6 @@ subroutine FindEvals(electrons) expH_atomf = zero end if end if -!!! 2025.02.03 nakata end !call gcopy(Efermi) !call gcopy(occ,matrix_size,nkp) ! else @@ -1031,13 +1029,11 @@ subroutine FindEvals(electrons) if (stat /= 0) call cq_abort('FindEvals: failed to deallocacte scaledEig', stat) call reg_dealloc_mem(area_DM, matrix_size * prim_size, type_cplx) end if -!!! 2025.02.03 nakata if(wf_self_con .and. (atomf.ne.sf)) then deallocate(expH_atomf, STAT=stat) if (stat /= 0) call cq_abort('FindEvals: failed to deallocacte expH_atomf', stat) call reg_dealloc_mem(area_DM, matrix_size * prim_size_atomf, type_cplx) end if -!!! 2025.02.03 nakata end ! global call endDiag min_layer = min_layer + 1 @@ -3324,6 +3320,8 @@ end function MP_entropy !! Add possibility to calculate S_{ij}.c^n_j for pDOS calculation post-processing !! 2023/06/08 15:38 dave !! Introduce flag_pol_build_S to select polarisation calculation + !! 2025/02/13 17:30 dave and nakata + !! Introduce flag_WFatomf_buildK and Eig_atomf for pDOS with MSSFs !! SOURCE !! subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, Eig_atomf, matSFcoeffTran) @@ -3342,13 +3340,13 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, ni_in_cell, x_atom_cell, y_atom_cell, & z_atom_cell, max_wf, min_layer, flag_do_pol_calc, mat_polX_re, mat_polX_im, & i_pol_dir_st, i_pol_dir_end, wf_self_con, flag_write_projected_DOS, ne_spin_in_cell, & - sf, atomf !!! 2025.02.03 nakata + sf, atomf use mpi use GenBlas, only: dot use GenComms, only: myid use mult_module, only: store_matrix_value_pos, matrix_pos, matK, return_matrix_value_pos use matrix_data, only: mat, halo - use species_module, only: nsf_species, natomf_species !!! 2025.02.03 nakata + use species_module, only: nsf_species, natomf_species implicit none @@ -3361,11 +3359,9 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, ! For pDOS complex(double_cplx), optional, dimension(:,:) :: overlapEig integer, optional :: matSij -!!! 2025.02.03 nakata ! For WFs with MSSFs complex(double_cplx), optional, dimension(:,:) :: Eig_atomf integer, optional :: matSFcoeffTran -!!! 2025.02.03 nakata end ! Local variables type(Krecv_data), dimension(:), allocatable :: recv_info @@ -3378,14 +3374,14 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, integer :: len, send_size, recv_size, send_proc, recv_proc, nsf1, & sendtag, recvtag integer :: req1, req2, ierr, atom, inter, prim, wheremat, row_sup,& - col_sup, col_atomf !!! 2025.02.03 nakata + col_sup, col_atomf integer, dimension(:,:), allocatable :: ints, atom_list, & send_prim, send_info, send_orbs, send_off integer, dimension(:), allocatable :: current_loc_atoms, & LocalAtom, num_send, norb_send, send_FSC, recv_to_FSC, & - mapchunk, prim_orbs, prim_orbs_atomf !!! 2025.02.03 nakata + mapchunk, prim_orbs, prim_orbs_atomf integer, dimension(MPI_STATUS_SIZE) :: mpi_stat - real(double) :: phase, rfac, ifac, rcc, icc, rsum, exp_X_value_real, exp_X_value_imag, Siajb, SFcoeffTran_iajb !!! 2025.02.03 nakata + real(double) :: phase, rfac, ifac, rcc, icc, rsum, exp_X_value_real, exp_X_value_imag, Siajb, SFcoeffTran_iajb complex(double_cplx) :: zsum, exp_X_value complex(double_cplx), dimension(:,:), allocatable :: RecvBuffer, & SendBuffer @@ -3486,7 +3482,7 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, write(io_lun,fmt='(10x,i6,a,3i6)') myid,' Maxima: ',maxloc, maxint, maxsend ! Allocate recv_info allocate(send_info(numprocs,maxsend),send_orbs(numprocs,maxsend),send_off(numprocs,maxsend), & - prim_orbs(bundle%mx_iprim),prim_orbs_atomf(bundle%mx_iprim),STAT=stat) !!! 2025.02.03 nakata + prim_orbs(bundle%mx_iprim),prim_orbs_atomf(bundle%mx_iprim),STAT=stat) if(stat/=0) call cq_abort('buildK: Error allocating send_info !',stat) send_info = 0 send_orbs = 0 @@ -3497,14 +3493,12 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, prim_orbs(j) = orb_count orb_count = orb_count + nsf_species(bundle%species(j)) end do -!!! 2025.02.03 nakata prim_orbs_atomf = 0 orb_count = 0 do j=1,bundle%n_prim prim_orbs_atomf(j) = orb_count orb_count = orb_count + natomf_species(bundle%species(j)) end do -!!! 2025.02.03 nakata end allocate(recv_info(numprocs),STAT=stat) if(stat/=0) call cq_abort('buildK: Error allocating recv_info !',stat) do i=1,numprocs @@ -3587,10 +3581,8 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, if(myid==0.AND.iprint_DM + min_layer>=4) write(io_lun,fmt='(10x,a,f8.4)') 'Occ is ',occs(i) end if end do -!!! 2025.02.03 nakata if(wf_self_con .and. flag_write_projected_DOS) len = matrix_size if(wf_self_con .and. flag_WFatomf_buildK) len = matrix_size -!!! 2025.02.03 nakata end len_occ = len if(iprint_DM+min_layer>3.AND.myid==0) & write(io_lun,fmt='(10x,a,2i6)') 'buildK: Stage three len:',len, matA @@ -3674,10 +3666,7 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, end do if(iprint_DM + min_layer>=4.AND.myid==0) & write(io_lun,fmt='(10x,a)') 'filling buffer' -!!! 2025.02.03 nakata -! if(.not.(wf_self_con .and. flag_write_projected_DOS)) then if(.not.(flag_pDOS_buildK .or. flag_WFatomf_buildK)) then -!!! 2025.02.03 nakata end do j=1,len_occ ! This is a loop over eigenstates RecvBuffer(j,1:recv_info(recv_proc+1)%orbs) = RecvBuffer(j,1:recv_info(recv_proc+1)%orbs)*occ_correction*occs(j) end do @@ -3728,23 +3717,11 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, polSloc(1:pol_S_size,1:pol_S_size,dir),pol_S_size) end do end if -!!! 2025.02.03 nakata pDOS-MSSF -! ! projected DOS -! if(flag_pDOS_buildK .and. wf_self_con .and. flag_write_projected_DOS) then -! whereMat = matrix_pos(matSij, prim, jatom, col_sup, row_sup) -! Siajb = return_matrix_value_pos(matSij,whereMat) -! ! 1:len_occ gives bands; we want c_jb^n * S_iajb -! overlapEig(1:len_occ,prim_orbs(prim)+col_sup) = & -! overlapEig(1:len_occ,prim_orbs(prim)+col_sup) + & -! Siajb*RecvBuffer(1:len_occ,orb_count+row_sup)*cmplx(rfac,ifac,double_cplx) -! end if -!!! 2025.02.03 nakata end ! Build K or M3 whereMat = matrix_pos(matA, prim, jatom, col_sup, row_sup) zsum = dot(len_occ,localEig(1:len_occ,prim_orbs(prim)+col_sup),1,RecvBuffer(1:len_occ,orb_count+row_sup),1) call store_matrix_value_pos(matA,whereMat,real(zsum*cmplx(rfac,ifac,double_cplx),double)) end do ! col_sup=nsf -!!! 2025.02.03 nakata pDOS-MSSF ! projected DOS in atomf basis if(wf_self_con) then do col_atomf = 1,natomf_species(bundle%species(prim)) @@ -3767,8 +3744,7 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, Siajb*RecvBuffer(1:len_occ,orb_count+row_sup)*cmplx(rfac,ifac,double_cplx) endif end do ! col_atomf=natomf - end if -!!! 2025.02.03 nakata end + end if ! wf_self_con (projected DOS) end do ! row_sup=nsf end do ! inter=recv_info%ints ! Careful - we only want to increment after ALL interactions done @@ -3801,7 +3777,7 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, recv_info(i)%dx,recv_info(i)%dy,recv_info(i)%dz,STAT=stat) if(stat/=0) call cq_abort('buildK: Error deallocating recvinfo !',i,stat) end do - deallocate(prim_orbs_atomf,prim_orbs,send_off,send_orbs,send_info,STAT=stat) !!! 2025.02.03 nakata + deallocate(prim_orbs_atomf,prim_orbs,send_off,send_orbs,send_info,STAT=stat) if(stat/=0) call cq_abort('buildK: Error allocating send_info !',stat) deallocate(recv_info,STAT=stat) if(stat/=0) call cq_abort('buildK: Error allocating recv_info !',stat) diff --git a/src/dimens_module.f90 b/src/dimens_module.f90 index a296ea04a..638cfbcba 100644 --- a/src/dimens_module.f90 +++ b/src/dimens_module.f90 @@ -163,6 +163,8 @@ module dimens !! Added checks to round RadiusAtomf and RadiusSupport to safe value (including grid points) !! 2024/07/18 14:18 lionel !! Check consistency of Xrange wrt r_exx read from input +!! 2025/02/03 14:00 nakata +!! Added aSs_in_sSs_range for pDOS with MSSFs !! SOURCE !! subroutine set_dimensions(inode, ionode,HNL_fac,non_local, n_species, non_local_species, core_radius) @@ -235,18 +237,13 @@ subroutine set_dimensions(inode, ionode,HNL_fac,non_local, n_species, non_local_ if(flag_neutral_atom_projector) then aNArange = 31 NAarange = 32 -!!! 2025.02.03 nakata -! mx_matrices_tmp = mx_matrices ! = 30 mx_matrices_tmp = 32 -!!! 2025.02.03 nakata end else mx_matrices_tmp = 30 end if -!!! 2025.02.03 nakata aSs_in_sSs_range = mx_matrices_tmp + 1 mx_matrices_tmp = mx_matrices_tmp + 1 if (mx_matrices_tmp > mx_matrices) call cq_abort('ERROR : mx_matrices_tmp is larger than mx_matrices',mx_matrices_tmp) -!!! 2025.02.03 nakata end endif !n_my_grid_points = n_pts_in_block * n_blocks @@ -379,7 +376,7 @@ subroutine set_dimensions(inode, ionode,HNL_fac,non_local, n_species, non_local_ rcut(SFcoeffTr_range) = 0.001_double endif if (abs(r_LD)