small restructuring of the routines for the AHC calculation with multiple Fermi energies and adaptive refinement

src/postw90/berry.F90

@@ -152,9 +152,10 @@ subroutine berry_main
     integer           :: n, i, j, k, jk, ikpt, if, ispn, ierr, loop_x, loop_y, loop_z, &
                          loop_xyz, loop_adpt, adpt_counter_list(nfermi), ifreq, &
                          file_unit
-    character(len=24) :: file_name
+    character(len=120) :: file_name
     logical           :: eval_ahc, eval_morb, eval_kubo, not_scannable, eval_sc, eval_shc
     logical           :: ladpt_kmesh
+    logical           :: ladpt(nfermi)
 
     if (nfermi == 0) call io_error( &
       'Must specify one or more Fermi levels when berry=true')
@@ -371,6 +372,7 @@ subroutine berry_main
         !
         if (eval_ahc) then
           call berry_get_imf_klist(kpt, imf_k_list)
+          ladpt = .false.
           do if = 1, nfermi
             vdum(1) = sum(imf_k_list(:, 1, if))
             vdum(2) = sum(imf_k_list(:, 2, if))
@@ -379,18 +381,25 @@ subroutine berry_main
             rdum = sqrt(dot_product(vdum, vdum))
             if (rdum > berry_curv_adpt_kmesh_thresh) then
               adpt_counter_list(if) = adpt_counter_list(if) + 1
-              do loop_adpt = 1, berry_curv_adpt_kmesh**3
-                ! Using imf_k_list here would corrupt values for other
-                ! frequencies, hence dummy. Only if-th element is used
-                call berry_get_imf_klist(kpt(:) + adkpt(:, loop_adpt), &
-                                         imf_k_list_dummy)
-                imf_list(:, :, if) = imf_list(:, :, if) &
-                                     + imf_k_list_dummy(:, :, if)*kweight_adpt
-              end do
+              ladpt(if) = .true.
             else
               imf_list(:, :, if) = imf_list(:, :, if) + imf_k_list(:, :, if)*kweight
             endif
           enddo
+          if (any(ladpt)) then
+            do loop_adpt = 1, berry_curv_adpt_kmesh**3
+              ! Using imf_k_list here would corrupt values for other
+              ! frequencies, hence dummy. Only if-th element is used
+              call berry_get_imf_klist(kpt(:) + adkpt(:, loop_adpt), &
+                                       imf_k_list_dummy, ladpt=ladpt)
+              do if = 1, nfermi
+                if (ladpt(if)) then
+                  imf_list(:, :, if) = imf_list(:, :, if) &
+                                       + imf_k_list_dummy(:, :, if)*kweight_adpt
+                endif
+              enddo
+            end do
+          endif
         end if
 
         if (eval_morb) then
@@ -495,6 +504,7 @@ subroutine berry_main
         !
         if (eval_ahc) then
           call berry_get_imf_klist(kpt, imf_k_list)
+          ladpt = .false.
           do if = 1, nfermi
             vdum(1) = sum(imf_k_list(:, 1, if))
             vdum(2) = sum(imf_k_list(:, 2, if))
@@ -503,18 +513,25 @@ subroutine berry_main
             rdum = sqrt(dot_product(vdum, vdum))
             if (rdum > berry_curv_adpt_kmesh_thresh) then
               adpt_counter_list(if) = adpt_counter_list(if) + 1
-              do loop_adpt = 1, berry_curv_adpt_kmesh**3
-                ! Using imf_k_list here would corrupt values for other
-                ! frequencies, hence dummy. Only if-th element is used
-                call berry_get_imf_klist(kpt(:) + adkpt(:, loop_adpt), &
-                                         imf_k_list_dummy)
-                imf_list(:, :, if) = imf_list(:, :, if) &
-                                     + imf_k_list_dummy(:, :, if)*kweight_adpt
-              end do
+              ladpt(if) = .true.
             else
               imf_list(:, :, if) = imf_list(:, :, if) + imf_k_list(:, :, if)*kweight
             endif
           enddo
+          if (any(ladpt)) then
+            do loop_adpt = 1, berry_curv_adpt_kmesh**3
+              ! Using imf_k_list here would corrupt values for other
+              ! frequencies, hence dummy. Only if-th element is used
+              call berry_get_imf_klist(kpt(:) + adkpt(:, loop_adpt), &
+                                       imf_k_list_dummy, ladpt=ladpt)
+              do if = 1, nfermi
+                if (ladpt(if)) then
+                  imf_list(:, :, if) = imf_list(:, :, if) &
+                                       + imf_k_list_dummy(:, :, if)*kweight_adpt
+                endif
+              enddo
+            end do
+          endif
         end if
 
         if (eval_morb) then
@@ -1139,7 +1156,7 @@ subroutine berry_main
 
   end subroutine berry_main
 
-  subroutine berry_get_imf_klist(kpt, imf_k_list, occ)
+  subroutine berry_get_imf_klist(kpt, imf_k_list, occ, ladpt)
     !============================================================!
     !                                                            !
     !! Calculates the Berry curvature traced over the occupied
@@ -1152,16 +1169,21 @@ subroutine berry_get_imf_klist(kpt, imf_k_list, occ)
     real(kind=dp), intent(in)                    :: kpt(3)
     real(kind=dp), intent(out), dimension(:, :, :) :: imf_k_list
     real(kind=dp), intent(in), optional, dimension(:) :: occ
+    logical, intent(in), optional, dimension(:) :: ladpt
 
     if (present(occ)) then
       call berry_get_imfgh_klist(kpt, imf_k_list, occ=occ)
     else
-      call berry_get_imfgh_klist(kpt, imf_k_list)
+      if (present(ladpt)) then
+        call berry_get_imfgh_klist(kpt, imf_k_list, ladpt=ladpt)
+      else
+        call berry_get_imfgh_klist(kpt, imf_k_list)
+      endif
     endif
 
   end subroutine berry_get_imf_klist
 
-  subroutine berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list, occ)
+  subroutine berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list, occ, ladpt)
     !=========================================================!
     !
     !! Calculates the three quantities needed for the orbital
@@ -1194,6 +1216,7 @@ subroutine berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list, occ)
     real(kind=dp), intent(out), dimension(:, :, :), optional &
       :: imf_k_list, img_k_list, imh_k_list
     real(kind=dp), intent(in), optional, dimension(:) :: occ
+    logical, intent(in), optional, dimension(:) :: ladpt
 
     complex(kind=dp), allocatable :: HH(:, :)
     complex(kind=dp), allocatable :: UU(:, :)
@@ -1208,6 +1231,7 @@ subroutine berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list, occ)
     real(kind=dp)                 :: eig(num_wann)
     integer                       :: i, j, ife, nfermi_loc
     real(kind=dp)                 :: s
+    logical                       :: todo(nfermi)
 
     ! Temporary space for matrix products
     complex(kind=dp), allocatable, dimension(:, :, :) :: tmp
@@ -1218,6 +1242,12 @@ subroutine berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list, occ)
       nfermi_loc = nfermi
     endif
 
+    if (present(ladpt)) then
+      todo = ladpt
+    else
+      todo = .true.
+    endif
+
     allocate (HH(num_wann, num_wann))
     allocate (UU(num_wann, num_wann))
     allocate (f_list(num_wann, num_wann, nfermi_loc))
@@ -1244,23 +1274,25 @@ subroutine berry_get_imfgh_klist(kpt, imf_k_list, img_k_list, imh_k_list, occ)
       ! Trace formula for -2Im[f], Eq.(51) LVTS12
       !
       do ife = 1, nfermi_loc
-        do i = 1, 3
-          !
-          ! J0 term (Omega_bar term of WYSV06)
-          imf_k_list(1, i, ife) = &
-            utility_re_tr_prod(f_list(:, :, ife), OOmega(:, :, i))
-          !
-          ! J1 term (DA term of WYSV06)
-          imf_k_list(2, i, ife) = -2.0_dp* &
-                                  ( &
-                                  utility_im_tr_prod(AA(:, :, alpha_A(i)), JJp_list(:, :, ife, beta_A(i))) &
-                                  + utility_im_tr_prod(JJm_list(:, :, ife, alpha_A(i)), AA(:, :, beta_A(i))) &
-                                  )
-          !
-          ! J2 term (DD of WYSV06)
-          imf_k_list(3, i, ife) = -2.0_dp* &
-                                  utility_im_tr_prod(JJm_list(:, :, ife, alpha_A(i)), JJp_list(:, :, ife, beta_A(i)))
-        end do
+        if (todo(ife)) then
+          do i = 1, 3
+            !
+            ! J0 term (Omega_bar term of WYSV06)
+            imf_k_list(1, i, ife) = &
+              utility_re_tr_prod(f_list(:, :, ife), OOmega(:, :, i))
+            !
+            ! J1 term (DA term of WYSV06)
+            imf_k_list(2, i, ife) = -2.0_dp* &
+                                    ( &
+                                    utility_im_tr_prod(AA(:, :, alpha_A(i)), JJp_list(:, :, ife, beta_A(i))) &
+                                    + utility_im_tr_prod(JJm_list(:, :, ife, alpha_A(i)), AA(:, :, beta_A(i))) &
+                                    )
+            !
+            ! J2 term (DD of WYSV06)
+            imf_k_list(3, i, ife) = -2.0_dp* &
+                                    utility_im_tr_prod(JJm_list(:, :, ife, alpha_A(i)), JJp_list(:, :, ife, beta_A(i)))
+          end do
+        endif
       end do
     end if
 

test-suite/tests/jobconfig

@@ -207,6 +207,12 @@ program = POSTW90_WPOUT_OK
 inputs_args = ('Fe.win', '')
 output = Fe.wpout 
 
+# Test of AHC for bcc iron with adaptive mesh refinement and Fermi scan
+[testpostw90_fe_ahc_adaptandfermi/]
+program = POSTW90_FERMISCAN_OK
+inputs_args = ('Fe.win', '')
+output = Fe-ahc-fermiscan.dat 
+
 # Test of Orbital Magnetisation for bcc iron
 [testpostw90_fe_morb/]
 program = POSTW90_WPOUT_OK

test-suite/tests/testpostw90_fe_ahc_adaptandfermi/Fe.amn

@@ -0,0 +1 @@
+../testpostw90_fe_ahc/Fe.amn

test-suite/tests/testpostw90_fe_ahc_adaptandfermi/Fe.chk.fmt.bz2

@@ -0,0 +1 @@
+../testpostw90_fe_ahc/Fe.chk.fmt.bz2

test-suite/tests/testpostw90_fe_ahc_adaptandfermi/Fe.eig

@@ -0,0 +1 @@
+../testpostw90_fe_ahc/Fe.eig

test-suite/tests/testpostw90_fe_ahc_adaptandfermi/Fe.mmn.bz2

@@ -0,0 +1 @@
+../testpostw90_fe_ahc/Fe.mmn.bz2

test-suite/tests/testpostw90_fe_ahc_adaptandfermi/Fe.win

@@ -0,0 +1,70 @@
+num_bands         =   28        
+num_wann          =   18 
+use_ws_distance = .false.
+search_shells=12
+
+dis_win_min       = -8.0d0       
+dis_win_max       = 70.0d0
+dis_froz_min      = -8.0d0
+dis_froz_max      =  30.0d0
+dis_num_iter      =  500 
+
+num_iter          = 200 
+dis_mix_ratio = 1.0
+
+spinors = true
+begin projections
+Fe: sp3d2;dxy;dxz;dyz
+end projections
+
+fermi_energy_min = 11.6279
+fermi_energy_max = 13.6279
+fermi_energy_step = 0.2
+
+uHu_formatted = .true.
+
+berry = true
+berry_task = ahc!+morb
+berry_kmesh = 10 10 10 
+
+berry_curv_adpt_kmesh = 5 5 5
+berry_curv_adpt_kmesh_thresh = 10
+
+!kpath = true
+!kpath_task = bands+morb
+!kpath_num_points=500
+begin kpoint_path
+G 0.0000 0.0000 0.0000   H 0.500 -0.5000 -0.5000
+H 0.500 -0.5000 -0.5000  P 0.7500 0.2500 -0.2500
+end kpoint_path
+
+#kslice = true
+#kslice_task = morb+fermi_lines
+#kslice_2dkmesh = 50 50
+#kslice_corner = 0.0  0.0  0.0
+#kslice_b1 =     0.5 -0.5 -0.5
+#kslice_b2 =     0.5  0.5  0.5
+
+begin unit_cell_cart
+bohr
+ 2.71175  2.71175 2.71175
+-2.71175  2.71175 2.71175
+-2.71175 -2.71175 2.71175
+end unit_cell_cart
+
+begin atoms_frac
+Fe  0.000  0.000  0.000
+end atoms_frac
+
+mp_grid           = 2 2 2 
+
+begin kpoints
+  0.00000000  0.00000000  0.00000000 
+  0.00000000  0.00000000  0.50000000 
+  0.00000000  0.50000000  0.00000000 
+  0.00000000  0.50000000  0.50000000 
+  0.50000000  0.00000000  0.00000000 
+  0.50000000  0.00000000  0.50000000 
+  0.50000000  0.50000000  0.00000000 
+  0.50000000  0.50000000  0.50000000 
+end kpoints

test-suite/tests/testpostw90_fe_ahc_adaptandfermi/Makefile

@@ -0,0 +1,14 @@
+all: $(patsubst %.chk.fmt.bz2, %.chk, $(wildcard *.chk.fmt.bz2)) $(patsubst %.mmn.bz2, %.mmn, $(wildcard *.mmn.bz2)) 
+
+%.chk: %.chk.fmt.bz2
+	$(eval SEEDNAME:=$(patsubst %.chk.fmt.bz2, %, $<))
+	echo $(SEEDNAME)
+	cat $< | bunzip2 - > $(SEEDNAME).chk.fmt && ../../../w90chk2chk.x -f2u $(SEEDNAME) && rm $(SEEDNAME).chk.fmt
+
+%.mmn: %.mmn.bz2
+	$(eval SEEDNAME:=$(patsubst %.mmn.bz2, %, $<))
+	echo $(SEEDNAME)
+	cat $< | bunzip2 - > $(SEEDNAME).mmn
+
+
+.PHONY: all

test-suite/tests/testpostw90_fe_ahc_adaptandfermi/benchmark.out.default.inp=Fe.win

@@ -0,0 +1,11 @@
+   11.627900    50.845203   -50.619458   157.962756
+   11.827900    20.541495   -20.578261   568.011017
+   12.027900    73.684613   -73.701303   533.909220
+   12.227900    32.332682   -32.350845   210.635328
+   12.427900     1.291984    -1.290424    11.323568
+   12.627900    41.926963   -41.913963   722.267035
+   12.827900    35.739944   -35.798745  2804.956845
+   13.027900    -5.553881     5.561349   124.550233
+   13.227900     4.113229    -4.105007   142.046112
+   13.427900    -1.363268     1.404604   -10.104289
+   13.627900     0.001088     0.006733    90.262635

test-suite/tests/userconfig

@@ -103,6 +103,13 @@ tolerance = ( (1.0e-3, 2.0e-3, 'ahc_x'),
                (1.0e-2, 2.0e-2, 'spin_p'),
                (1.0e-2, 2.0e-2, 'spin_a'))
 
+[POSTW90_FERMISCAN_OK]
+exe = ../../postw90.x
+extract_fn = tools parsers.parse_fermiscan_dat.parse
+tolerance = (  (1.0e-6, 5.0e-6, 'fermienergy'),
+               (1.0e-4, 2.0e-4, 'ahc_x'),
+               (1.0e-4, 2.0e-4, 'ahc_y'),
+               (1.0e-4, 2.0e-4, 'ahc_z'))
 
 [POSTW90_JDOS_OK]
 exe = ../../postw90.x

test-suite/tools/parsers/parse_fermiscan_dat.py

@@ -0,0 +1,52 @@
+"""
+Parser function parse() to parse the <seedname>_fermiscan.dat output file of postw90.x.
+"""
+from __future__ import print_function
+
+import inspect
+import re
+from collections import defaultdict
+
+from . import show_output
+
+def parse(fname):
+    """
+    Open the file, parses it and return the values.
+    """
+    retdict = defaultdict(list)
+
+    if show_output:
+        print("[{}.{}] Parsing file '{}'".format(
+            __name__, inspect.currentframe().f_code.co_name, fname))
+
+    with open(fname) as f:
+        lines = f.readlines()
+
+    for lno, l in enumerate(lines):
+
+        if l.strip().startswith('#'):
+            # Skip headers
+            continue
+
+        pieces = l.split()
+
+        if len(pieces) == 0 :
+            # skip blank line
+            continue
+
+        if len(pieces) == 4 : 
+            retdict['fermienergy'].append(float(pieces[0]))
+            retdict['ahc_x'].append(float(pieces[1]))
+            retdict['ahc_y'].append(float(pieces[2]))
+            retdict['ahc_z'].append(float(pieces[3]))
+        else:
+            raise ValueError("Wrong line length ({}, instead of 4); line content: {}".format(
+                len(pieces)), l)
+
+
+    retdict = dict(retdict)
+    if show_output:
+        for k in sorted(retdict):
+            print("  {}: {}".format(k, retdict[k]))
+        print("-"*72)
+    return retdict