From 1fe25884b58a274f3d44fa979318cf2b5b923fc9 Mon Sep 17 00:00:00 2001
From: Philipp Ruessmann <p.ruessmann@fz-juelich.de>
Date: Wed, 1 Dec 2021 17:42:02 +0100
Subject: [PATCH] Add write_double_precision run option for KKRhost and
 read_doubleprecision testflag for kkrimp to use doble precision kkrflex_green
 files instead of single precision files

---
 source/KKRhost/main1b.F90         | 11 +++-
 source/KKRhost/rotgll.F90         |  9 ++-
 source/KKRimp/preconditioning.F90 | 99 +++++++++++++++++--------------
 source/common/runoptions.F90      |  4 ++
 4 files changed, 75 insertions(+), 48 deletions(-)

diff --git a/source/KKRhost/main1b.F90 b/source/KKRhost/main1b.F90
index 657266646..f7dda9d80 100644
--- a/source/KKRhost/main1b.F90
+++ b/source/KKRhost/main1b.F90
@@ -37,7 +37,7 @@ contains
     use :: mod_datatypes, only: dp, sp
     use :: mod_runoptions, only: calc_exchange_couplings, formatted_file, set_gmat_to_zero, use_Chebychev_solver, &
       use_qdos, use_readcpa, write_deci_tmat, write_gmat_plain, write_green_host, write_green_imp, write_kkrimp_input, &
-      write_pkkr_input, write_pkkr_operators, write_rhoq_input, write_gmat_ascii, decouple_spin_cheby
+      write_pkkr_input, write_pkkr_operators, write_rhoq_input, write_gmat_ascii, decouple_spin_cheby, write_double_precision
     use mod_constants, only: czero, cone, pi
     use mod_operators_for_fscode, only: operators_for_fscode
     use mod_getscratch, only: opendafile
@@ -357,8 +357,13 @@ contains
       if ((write_kkrimp_input)) then
         ! ! Green functions has (lmmaxd*natomimp)**2 double complex (i.e. factor '4') values
         ! RECLENGTH = WLENGTH*4*NATOMIMP*LMMAXD*NATOMIMP*LMMAXD
-        ! at the moment kkrflex_green file is only written with single precision (factor'2')
-        reclength = wlength*2*natomimp*lmmaxd*natomimp*lmmaxd
+        if ( write_double_precision ) then
+          ! this is used to write double precision files (factor'4')
+          reclength = wlength*4*natomimp*lmmaxd*natomimp*lmmaxd
+        else
+          ! at the moment kkrflex_green file is only written with single precision (factor'2')
+          reclength = wlength*2*natomimp*lmmaxd*natomimp*lmmaxd
+        end if
         ! sometimes (lmax=2) the record length might be too small to store the parameters, then reclength needs to be bigger
         if (reclength<8*ielast+6) then
           stop '[main1b] record length for kkrflex_green is too small to store parameters, use either more atoms in the cluster, a higher lmax or less energy points'
diff --git a/source/KKRhost/rotgll.F90 b/source/KKRhost/rotgll.F90
index fbb0c3f8f..c396d0629 100644
--- a/source/KKRhost/rotgll.F90
+++ b/source/KKRhost/rotgll.F90
@@ -53,7 +53,7 @@ contains
 
     use :: mod_mympi, only: myrank, master
     use :: mod_datatypes, only: dp, sp
-    use :: mod_runoptions, only: print_Gij, write_gmat_plain, write_green_host, write_kkrimp_input
+    use :: mod_runoptions, only: print_Gij, write_gmat_plain, write_green_host, write_kkrimp_input, write_double_precision
     use :: mod_changerep
     use :: mod_cmatstr
     use :: mod_constants, only: czero, cone
@@ -193,7 +193,12 @@ contains
 #endif
         ! force single precision complex writeout to minimize file size etc.
         ! maybe this can be removed in the future
-        write (888, rec=irec) cmplx(gclust, kind=sp)
+        if ( write_double_precision ) then
+          write (888, rec=irec) gclust
+        else
+          write (888, rec=irec) cmplx(gclust, kind=sp)
+        end if
+
         if ((write_gmat_plain)) then
           write (8888, '(50000E25.16)') gclust
         end if
diff --git a/source/KKRimp/preconditioning.F90 b/source/KKRimp/preconditioning.F90
index 7e1215a12..966d02dfc 100644
--- a/source/KKRimp/preconditioning.F90
+++ b/source/KKRimp/preconditioning.F90
@@ -2,51 +2,51 @@ module mod_preconditioning
   contains
    subroutine preconditioning_start(my_rank,mpi_size,ez, wez, ielast, intercell_ach, alat,vmtzero,lattice_relax,gmatbulk) 
 #ifdef CPP_MPI
-  use mpi
+    use mpi
 #endif
-     use nrtype
-     use mod_dysonvirtatom, only: dysonvirtatom
-     use mod_read_atominfo
-     use mod_config, only: config_testflag,config_runflag
-     use mod_log, only: log_write
-     use mod_mpienergy, only: mpienergy_distribute
-     use type_gmatbulk
-     use mod_types, only: t_inc
-
-     implicit none
+    use nrtype
+    use mod_dysonvirtatom, only: dysonvirtatom
+    use mod_read_atominfo
+    use mod_config, only: config_testflag,config_runflag
+    use mod_log, only: log_write
+    use mod_mpienergy, only: mpienergy_distribute
+    use type_gmatbulk
+    use mod_types, only: t_inc
+
+    implicit none
 !interface
-     double complex,allocatable,intent(out)  :: ez(:)
-     double complex,allocatable,intent(out)  :: wez(:)
-     integer,intent(out)                     :: ielast
-     real(kind=DP),allocatable,intent(out)   :: intercell_ach(:,:)   ! intercell potential
-                                               !(lmpotd,ntotatom),achnew(lmpotd,ntotatom)
-     real(kind=DP),intent(out)               :: alat
-     real(kind=dp)                           :: vmtzero(2)  
-     integer                                 :: lattice_relax
-     type(gmatbulk_type)                     :: gmatbulk
+    double complex,allocatable,intent(out)  :: ez(:)
+    double complex,allocatable,intent(out)  :: wez(:)
+    integer,intent(out)                     :: ielast
+    real(kind=DP),allocatable,intent(out)   :: intercell_ach(:,:)   ! intercell potential
+    real(kind=DP),intent(out)               :: alat
+    real(kind=dp)                           :: vmtzero(2)  
+    integer                                 :: lattice_relax
+    type(gmatbulk_type)                     :: gmatbulk
 !local
-
-     integer                     :: natom                ! number of impurity atoms
-     integer                     :: ntotatom             ! number of imp atoms+killatoms
-     integer,allocatable         :: lmaxatom(:)          ! lmax for all atoms
-     integer,allocatable         :: isvatom(:)           ! 1=vatom (delta t=0), 0=no vatom
-     integer,allocatable         :: killatom(:)          ! 1=atom will be removed in a dyson step
-     integer                     :: KGREFSOC             ! =1 if SOC for the GREF 
-     integer                     :: NSOC                 ! =2 if SOC for the GREF else: =1 
-
-     integer                     :: ie
-     integer                     :: natomimpd,lmsizehost
-     double complex,allocatable  :: gmathost(:,:)
-     double complex,allocatable  :: gmathostnew(:,:)
-     integer                     :: recl1,recl2,nspin,ispin,iatom
-     double complex,allocatable  ::  tmat(:,:,:,:,:)
-     integer                     :: nlmhostnew
-     real(kind=dp),allocatable   :: RIMPATOM(:,:),zatom(:)
-     integer                      :: lmaxd
+    integer                     :: natom                ! number of impurity atoms
+    integer                     :: ntotatom             ! number of imp atoms+killatoms
+    integer,allocatable         :: lmaxatom(:)          ! lmax for all atoms
+    integer,allocatable         :: isvatom(:)           ! 1=vatom (delta t=0), 0=no vatom
+    integer,allocatable         :: killatom(:)          ! 1=atom will be removed in a dyson step
+    integer                     :: KGREFSOC             ! =1 if SOC for the GREF 
+    integer                     :: NSOC                 ! =2 if SOC for the GREF else: =1 
+
+    integer                     :: ie
+    integer                     :: natomimpd,lmsizehost
+    double complex,allocatable  :: gmathost(:,:)
+    double complex,allocatable  :: gmathostnew(:,:)
+    integer                     :: recl1,recl2,nspin,ispin,iatom
+    double complex,allocatable  ::  tmat(:,:,:,:,:)
+    integer                     :: nlmhostnew
+    real(kind=dp),allocatable   :: RIMPATOM(:,:),zatom(:)
+    integer                     :: lmaxd
+    character (len=15)          :: precision_mode
 !mpi
-      integer,allocatable                      :: mpi_iebounds(:,:)
-      integer                                  :: my_rank
-      integer                                  :: mpi_size,ierror
+    integer, allocatable        :: mpi_iebounds(:,:)
+    integer                     :: my_rank
+    integer                     :: mpi_size,ierror
+
 
 call log_write('>>>>>>>>>>>> preconditioning read_atominfo >>>>>>>>>>>>>>>>>>>>>')
 call read_atominfo('total','kkrflex_atominfo',natom,ntotatom,RIMPATOM,&
@@ -170,6 +170,7 @@ allocate(gmathostnew(nlmhostnew,nlmhostnew))
 ! ###########################################
 
 recl1=wlength*2*natomimpd*lmsizehost*natomimpd*lmsizehost
+if (config_testflag('read_doubleprecision')) recl1 = recl1 * 2
 recl2=wlength*4*nlmhostnew*nlmhostnew
 open (88,access='direct',recl=recl1,file='kkrflex_green',form='unformatted')
 if (lattice_relax==0) then
@@ -207,7 +208,14 @@ do ie=mpi_iebounds(1,my_rank),mpi_iebounds(2,my_rank)
   do ispin=1,nspin-KGREFSOC
     if (t_inc%i_write>0) write(1337,*) 'proc = ',my_rank,' IE = ',ie,' ispin= ',ispin
 
-    call preconditioning_readgreenfn(ie,ispin,ielast,lmsizehost,ntotatom,gmathost,'singleprecision')
+    if ( config_testflag('read_doubleprecision') ) then
+      if (t_inc%i_write>0) write(*,*) 'ATTENTION: read kkrflex_green in double precision'
+      precision_mode = 'doubleprecision'
+    else
+      precision_mode = 'singleprecision'
+    end if
+
+    call preconditioning_readgreenfn(ie,ispin,ielast,lmsizehost,ntotatom,gmathost,precision_mode)
 !                                      in    in         in        in     out
     if(config_testflag('gtest')) write(10000+my_rank,'(832E25.14)') gmathost
 
@@ -492,7 +500,11 @@ subroutine preconditioning_readenergy(my_rank,IELAST,NSPIN,EZ,WEZ,NATOMIMPD,NTOT
   ! *******************************************************************
   if (my_rank==0) then 
     write(*,*) 'my_rank=1 reads kkrflex_green and communicates to other processes'
-    OPEN (88,ACCESS='direct',RECL=wlength*2*16,FILE='kkrflex_green',FORM='unformatted')
+    if (config_testflag('read_doubleprecision')) then
+      OPEN (88,ACCESS='direct',RECL=wlength*4*16,FILE='kkrflex_green',FORM='unformatted')
+    else
+      OPEN (88,ACCESS='direct',RECL=wlength*2*16,FILE='kkrflex_green',FORM='unformatted')
+    end if
       read(88,rec=1) ielasttemp,nspintemp,natomimpd,natomimp,lmsizehosttemp,kgrefsoctemp
     
        if (kgrefsoctemp/=kgrefsoc) then
@@ -540,6 +552,7 @@ subroutine preconditioning_readenergy(my_rank,IELAST,NSPIN,EZ,WEZ,NATOMIMPD,NTOT
 !   *******************************************************************
     allocate( ez(ielast), wez(ielast) )
     RECL1=wlength*2*NATOMIMPD*lmsizehost*NATOMIMPD*lmsizehost
+    if (config_testflag('read_doubleprecision')) RECL1 = 2 * RECL1
     OPEN (88,ACCESS='direct',RECL=RECL1,FILE='kkrflex_green',FORM='unformatted')
     if ( .not. config_runflag('oldJMcode') ) then
       read(88,rec=1) ielasttemp,nspintemp,natomimpd,NATOMIMP,lmsizehosttemp,kgrefsoctemp,(ez(ie),ie=1,ielast),(wez(ie),ie=1,ielast)
diff --git a/source/common/runoptions.F90 b/source/common/runoptions.F90
index 47a17bf93..e4e73da64 100644
--- a/source/common/runoptions.F90
+++ b/source/common/runoptions.F90
@@ -122,6 +122,7 @@ module mod_runoptions
   logical :: use_broyden_spinmix = .false.             !! use broyden spin mixing for noncollinear angles
   logical :: write_angles_alliter= .false.             !! write out noncollinear angles for all iterations
   logical :: write_tmat_all= .false.                   !! write out the tmat for all atoms and energies
+  logical :: write_double_precision= .false.           !! write out kkrflex files in double precision
 
   !some old run and test options have been removed:
   !  'atptshft': replaced by presence or absence of IVSHIFT in inputcard
@@ -245,6 +246,7 @@ module mod_runoptions
     call set_runoption(use_broyden_spinmix           , '<use_broyden_spinmix>'            , '<bryspinmix>')
     call set_runoption(write_angles_alliter          , '<write_angles_alliter>')
     call set_runoption(write_tmat_all                , '<write_tmat_all>')
+    call set_runoption(write_double_precision        , '<write_double_precision>')
 
   end subroutine read_runoptions
 
@@ -789,6 +791,7 @@ module mod_runoptions
     call mpi_bcast(use_broyden_spinmix           , 1, mpi_logical, master, mpi_comm_world, ierr)
     call mpi_bcast(write_angles_alliter          , 1, mpi_logical, master, mpi_comm_world, ierr)
     call mpi_bcast(write_tmat_all                , 1, mpi_logical, master, mpi_comm_world, ierr)
+    call mpi_bcast(write_double_precision        , 1, mpi_logical, master, mpi_comm_world, ierr)
 
   end subroutine bcast_runoptions
 #endif
@@ -905,6 +908,7 @@ module mod_runoptions
     write(iounit, '(A35,1x,1L,3x,A)') '<use_broyden_spinmix>=', use_broyden_spinmix, "use broyden spin mixing for noncollinear angles"
     write(iounit, '(A35,1x,1L,3x,A)') '<write_angles_alliter>=', write_angles_alliter, "write out noncollinear angles for all iterations"
     write(iounit, '(A35,1x,1L,3x,A)') '<write_tmat_all>=', write_tmat_all, "write out the tmat for all energies and all atoms"
+    write(iounit, '(A35,1x,1L,3x,A)') '<write_double_precision>=', write_double_precision, "write out kkrflex files in double precision"
 
   end subroutine print_runoptions
 
-- 
GitLab