Merge pull request #403 from OrderN/f-pDOS_MSSFs

PDOS with MSSFs

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,10 +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
-          mx_matrices_tmp = mx_matrices ! = 30
+          mx_matrices_tmp = 32
        else
           mx_matrices_tmp = 30
        end if
+       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)
     endif
 
     !n_my_grid_points = n_pts_in_block * n_blocks    
@@ -371,6 +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)<very_small) rcut(LD_range) = 0.001_double
+       rcut(aSs_in_sSs_range) = rcut(Srange)
     endif
     if(flag_neutral_atom_projector) then
        rcut(aNArange)   = r_s_atomf + r_h_atomf
@@ -418,6 +424,7 @@ subroutine set_dimensions(inode, ionode,HNL_fac,non_local, n_species, non_local_
        mat_name(SFcoeff_range)   = "MS"
        mat_name(SFcoeffTr_range) = "MSt"
        mat_name(LD_range)        = "LD"
+       mat_name(aSs_in_sSs_range) = "aSs_in_sSs"
     endif
     if(inode==ionode.AND.iprint_init>1) then
        do n=1,mx_matrices_tmp

src/initial_read_module.f90

@@ -807,6 +807,8 @@ end subroutine read_and_write
   !!     Update test for solution method (diagon vs ordern) following issue #47
   !!   2024/12/03 lionel
   !!     Added grid specification of EXX coarse/standard/fine
+  !!   2025/02/03 nakata
+  !!     Set flag_out_wf = .true. expricitly when flag_write_projected_DOS is .true.
   !!  TODO
   !!  SOURCE
   !!
@@ -1765,6 +1767,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.
              E_wf_min = fdf_double('IO.min_wf_E',-BIG)
              E_wf_max = fdf_double('IO.max_wf_E',BIG)
           end if

src/matrix_data_module.f90

@@ -62,6 +62,9 @@
 !!    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 and aSs_in_sSs_matind were added
+!!    and mx_matrices was changed from 32 to 33 for pDOS with MSSFs
 !!  SOURCE
 !!
 module matrix_data
@@ -73,7 +76,7 @@ module matrix_data
   save
 
   ! This will need to change if the above parameters are changed
-  integer, parameter :: mx_matrices = 32
+  integer, parameter :: mx_matrices = 33
 
   ! Store ALL indices in a large array
   type(matrix),      allocatable, dimension(:,:), target :: mat
@@ -86,7 +89,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
   integer, dimension(:), pointer :: aNAmatind, NAamatind
 
   ! Parameters for the different matrix ranges
@@ -126,6 +130,7 @@ module matrix_data
   ! Ranges for NA projectors set later also (dimens.module.f90)
   integer :: aNArange        ! 31
   integer :: NAarange        ! 32
+  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 :: max_range ! Indexes matrix with largest range
 

src/mult_module.f90

@@ -266,6 +266,8 @@ module mult_module
   !!    Removed dSrange, dHrange, PAOPrange and PAOP_PS_H which are no longer used
   !!   2017/12/05 10:24 dave with TM and NW (Mizuho)
   !!    Adding initialisation for NA projector matrices
+  !!   2025/02/06 15:00 nakata
+  !!    Added aSs_in_sSs_range
   !!  SOURCE
   !!
   subroutine immi(parts, prim, gcs, myid, partial)
@@ -478,6 +480,11 @@ 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))
+       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))
        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
@@ -1397,6 +1404,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)
     endif
     if( flag_neutral_atom_projector ) then
        call end_ops(prim,aNArange, aNAmatind,aNA_trans)

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,9 @@ subroutine setup_pseudo_info
           call read_ion_ascii_tmp(pseudo(ispecies),pao(ispecies))
 
           npao_species(ispecies) = pao(ispecies)%count
-          ! 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)
+          ! because eigenvectors (ProcessWF and ProcessSijWF) output by CQ are always in PAO basis
+          nsf_species(ispecies) = pao(ispecies)%count
           maxnsf = max(maxnsf,nsf_species(ispecies))
           ! Find radius for atom functions
           do l=0,pao(ispecies)%greatest_angmom

tools/PostProcessing/read_module.f90

@@ -118,15 +118,12 @@ 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")
     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")
        ! Allow user to specify output filename
        root_file = fdf_string(50,'Process.RootFile','STM')
     else if(leqi(job,'stm')) then
        i_job = 5
-       if(flag_Multisite) call cq_abort("Not yet compatible with multi-site support functions")
        ! Allow user to specify output filename
        root_file = fdf_string(50,'Process.RootFile','STM')
     else if(leqi(job,'dos')) then