implemented: read FPLO grids

src/fieldmod@proc.f90

@@ -292,8 +292,8 @@ module subroutine field_new(f,seed,c,id,sptr,errmsg,ti)
     use param, only: ifformat_unknown, ifformat_wien, ifformat_elk, ifformat_pi,&
        ifformat_cube, ifformat_bincube, ifformat_abinit, ifformat_vasp,&
        ifformat_vaspnov, ifformat_qub,&
-       ifformat_xsf, ifformat_elkgrid, ifformat_siestagrid, ifformat_dftb,&
-       ifformat_pwc, ifformat_fmt,&
+       ifformat_xsf, ifformat_elkgrid, ifformat_siestagrid, ifformat_fplogrid,&
+       ifformat_dftb, ifformat_pwc, ifformat_fmt,&
        ifformat_wfn, ifformat_wfx, ifformat_fchk, ifformat_molden, ifformat_as,&
        ifformat_as_promolecular, ifformat_as_core, ifformat_as_lap, ifformat_as_grad,&
        ifformat_as_pot, ifformat_as_resample,&
@@ -435,6 +435,13 @@ module subroutine field_new(f,seed,c,id,sptr,errmsg,ti)
        f%type = type_grid
        f%file = seed%file(1)
 
+    elseif (seed%iff == ifformat_fplogrid) then
+       if (.not.allocated(f%grid)) allocate(f%grid)
+       call f%grid%end()
+       call f%grid%read_fplo(c_loc(c),seed%file(1),c%m_x2c,errmsg,ti=ti)
+       f%type = type_grid
+       f%file = seed%file(1)
+
     elseif (seed%iff == ifformat_dftb) then
        if (.not.allocated(f%dftb)) allocate(f%dftb)
        call f%dftb%end()

src/fieldseedmod@proc.f90

@@ -64,7 +64,7 @@ module subroutine fieldseed_parse(f,line,withoptions,lp0)
        ifformat_unknown, ifformat_wien, ifformat_elk, ifformat_pi, ifformat_cube,&
        ifformat_bincube, ifformat_abinit, ifformat_fmt,&
        ifformat_vasp, ifformat_vaspnov, ifformat_qub, ifformat_xsf, ifformat_elkgrid,&
-       ifformat_siestagrid, ifformat_dftb, ifformat_pwc,&
+       ifformat_siestagrid, ifformat_fplogrid, ifformat_dftb, ifformat_pwc,&
        ifformat_wfn, ifformat_wfx, ifformat_fchk,&
        ifformat_molden, ifformat_as, ifformat_as_promolecular, ifformat_as_core, ifformat_as_lap,&
        ifformat_as_resample,&
@@ -126,6 +126,9 @@ module subroutine fieldseed_parse(f,line,withoptions,lp0)
     elseif (equal(lfile,"siesta")) then
        f%iff = ifformat_siestagrid
        call read_next_as_file()
+    elseif (equal(lfile,"fplo")) then
+       f%iff = ifformat_fplogrid
+       call read_next_as_file()
     elseif (equal(lfile,"dftb")) then
        f%iff = ifformat_dftb
        call read_next_as_file()
@@ -193,7 +196,7 @@ module subroutine fieldseed_parse(f,line,withoptions,lp0)
        ! no files needed
        nfile = 0
     elseif (f%iff == ifformat_cube .or. f%iff == ifformat_bincube .or.&
-       f%iff == ifformat_abinit .or. f%iff == ifformat_siestagrid .or.&
+       f%iff == ifformat_abinit .or. f%iff == ifformat_siestagrid .or. f%iff == ifformat_fplogrid .or.&
        f%iff == ifformat_vasp .or. f%iff == ifformat_vaspnov .or. f%iff == ifformat_qub .or.&
        f%iff == ifformat_xsf .or. f%iff == ifformat_wfn .or. f%iff == ifformat_wfx .or.&
        f%iff == ifformat_fmt .or.&
@@ -647,7 +650,7 @@ module function field_detect_format(file,molden_type)
     use wfn_private, only: molden_type_unknown, molden_type_orca
     use tools_io, only: ferror, equal
     use param, only: dirsep,&
-       ifformat_cube,ifformat_bincube,ifformat_abinit,ifformat_siestagrid,&
+       ifformat_cube,ifformat_bincube,ifformat_abinit,ifformat_siestagrid,ifformat_fplogrid,&
        ifformat_dftb,ifformat_vasp,ifformat_vaspnov,ifformat_qub,ifformat_fmt,&
        ifformat_xsf,ifformat_wfn,ifformat_wfx,ifformat_fchk,ifformat_molden,&
        ifformat_molden,ifformat_wien,ifformat_elkgrid,ifformat_elk,&
@@ -687,6 +690,8 @@ module function field_detect_format(file,molden_type)
        equal(extdot,'DRHO') .or. equal(extdot,'LDOS') .or.&
        equal(extdot,'VT') .or. equal(extdot,'VH')) then
        field_detect_format = ifformat_siestagrid
+    else if (equal(extdot,'001')) then
+       field_detect_format = ifformat_fplogrid
     else if (equal(extdot,'xml')) then
        field_detect_format = ifformat_dftb
     else if (equal(extdot,'CHG').or.equal(extdot,'CHGCAR').or.&

src/grid3mod.f90

@@ -95,6 +95,7 @@ module grid3mod
      procedure :: read_cube !< grid3 from a Gaussian cube file
      procedure :: read_bincube !< grid3 from a binary cube file
      procedure :: read_siesta !< grid3 from siesta RHO file
+     procedure :: read_fplo !< grid3 from FPLO 001 file
      procedure :: read_abinit !< grid3 from abinit binary file
      procedure :: read_vasp !< grid3 from VASP file (CHG, CHGCAR, etc.)
      procedure :: read_qub !< grid3 from aimpac qub format
@@ -168,6 +169,14 @@ module subroutine read_siesta(f,cptr,file,x2c,errmsg,ti)
        character(len=:), allocatable, intent(out) :: errmsg
        type(thread_info), intent(in), optional :: ti
      end subroutine read_siesta
+     module subroutine read_fplo(f,cptr,file,x2c,errmsg,ti)
+       class(grid3), intent(inout) :: f
+       type(c_ptr), intent(in) :: cptr
+       character*(*), intent(in) :: file
+       real*8, intent(in) :: x2c(3,3)
+       character(len=:), allocatable, intent(out) :: errmsg
+       type(thread_info), intent(in), optional :: ti
+     end subroutine read_fplo
      module subroutine read_abinit(f,cptr,file,x2c,errmsg,ti)
        class(grid3), intent(inout) :: f
        type(c_ptr), intent(in) :: cptr

src/grid3mod@proc.f90

@@ -659,6 +659,86 @@ module subroutine read_siesta(f,cptr,file,x2c,errmsg,ti)
 
   end subroutine read_siesta
 
+  !> Read a grid in FPLO grid format
+  module subroutine read_fplo(f,cptr,file,x2c,errmsg,ti)
+    use tools_io, only: fopen_read, fclose, getline_raw, isinteger
+    use tools_math, only: matinv
+    class(grid3), intent(inout) :: f
+    type(c_ptr), intent(in) :: cptr
+    character*(*), intent(in) :: file !< Input file
+    real*8, intent(in) :: x2c(3,3)
+    character(len=:), allocatable, intent(out) :: errmsg
+    type(thread_info), intent(in), optional :: ti
+
+    logical :: ok
+    integer :: luc, nspin, istat
+    integer :: i, ix, iy, iz, n(3), ll, lp
+    real*8 :: r(3,3)
+    real*4, allocatable :: g(:)
+    character(len=:), allocatable :: line
+
+    errmsg = "Error reading file"
+
+    ! initialize
+    call f%end()
+    f%isinit = .true.
+    f%isqe = .false.
+    f%iswan = .false.
+    f%mode = mode_default
+
+    ! open file
+    luc = fopen_read(file,ti=ti)
+    if (luc < 0) goto 999
+
+    ! read the header and initialize geometry
+    n = 0
+    do while (getline_raw(luc,line))
+       ll = len(line)
+       if (ll >= 17) then
+          if (line(1:17) == "#position-subdiv:") then
+             line = line(18:)
+             lp = 1
+             ok = isinteger(n(1),line,lp)
+             ok = ok .and. isinteger(n(2),line,lp)
+             ok = ok .and. isinteger(n(3),line,lp)
+             if (.not.ok) then
+                errmsg = "Error reading number of grid points."
+                goto 999
+             end if
+          elseif (ll >= 9) then
+             if (line(1:9) == "#columns:") exit
+          end if
+       end if
+    end do
+    if (any(n <= 0)) then
+       errmsg = "Number of grid points not found."
+       goto 999
+    end if
+    call init_geometry(f,x2c,n,cptr)
+
+    ! read the field and close
+    allocate(f%f(n(1),n(2),n(3)),stat=istat)
+    if (istat /= 0) then
+       errmsg = "Error allocating grid"
+       goto 999
+    end if
+    f%f = 0d0
+    do iz = 1, n(3)
+       do iy = 1, n(2)
+          do ix = 1, n(1)
+             read (luc,*,err=999,end=999) f%f(ix,iy,iz)
+          end do
+       end do
+    end do
+    call fclose(luc)
+
+    errmsg = ""
+    return
+999 continue
+    if (luc > 0) call fclose(luc)
+
+  end subroutine read_fplo
+
   !> Read a grid in abinit format
   module subroutine read_abinit(f,cptr,file,x2c,errmsg,ti)
     use tools_math, only: matinv

src/param.F90

@@ -144,29 +144,30 @@ module param
   integer, parameter, public :: ifformat_xsf = 10
   integer, parameter, public :: ifformat_elkgrid = 11
   integer, parameter, public :: ifformat_siestagrid = 12
-  integer, parameter, public :: ifformat_dftb = 13
-  integer, parameter, public :: ifformat_pwc = 14
-  integer, parameter, public :: ifformat_wfn = 15
-  integer, parameter, public :: ifformat_wfx = 16
-  integer, parameter, public :: ifformat_fchk = 17
-  integer, parameter, public :: ifformat_molden = 18
-  integer, parameter, public :: ifformat_fmt = 19
-  integer, parameter, public :: ifformat_as = 20
-  integer, parameter, public :: ifformat_as_promolecular = 21
-  integer, parameter, public :: ifformat_as_core = 22
-  integer, parameter, public :: ifformat_as_lap = 23
-  integer, parameter, public :: ifformat_as_grad = 24
-  integer, parameter, public :: ifformat_as_pot = 25
-  integer, parameter, public :: ifformat_as_resample = 26
-  integer, parameter, public :: ifformat_as_clm = 27
-  integer, parameter, public :: ifformat_as_clm_sub = 28
-  integer, parameter, public :: ifformat_as_ghost = 29
-  integer, parameter, public :: ifformat_copy = 30
-  integer, parameter, public :: ifformat_promolecular = 31
-  integer, parameter, public :: ifformat_promolecular_fragment = 32
-  integer, parameter, public :: ifformat_as_hxx1 = 33
-  integer, parameter, public :: ifformat_as_hxx2 = 34
-  integer, parameter, public :: ifformat_as_hxx3 = 35
+  integer, parameter, public :: ifformat_fplogrid = 13
+  integer, parameter, public :: ifformat_dftb = 14
+  integer, parameter, public :: ifformat_pwc = 15
+  integer, parameter, public :: ifformat_wfn = 16
+  integer, parameter, public :: ifformat_wfx = 17
+  integer, parameter, public :: ifformat_fchk = 18
+  integer, parameter, public :: ifformat_molden = 19
+  integer, parameter, public :: ifformat_fmt = 20
+  integer, parameter, public :: ifformat_as = 21
+  integer, parameter, public :: ifformat_as_promolecular = 22
+  integer, parameter, public :: ifformat_as_core = 23
+  integer, parameter, public :: ifformat_as_lap = 24
+  integer, parameter, public :: ifformat_as_grad = 25
+  integer, parameter, public :: ifformat_as_pot = 26
+  integer, parameter, public :: ifformat_as_resample = 27
+  integer, parameter, public :: ifformat_as_clm = 28
+  integer, parameter, public :: ifformat_as_clm_sub = 29
+  integer, parameter, public :: ifformat_as_ghost = 30
+  integer, parameter, public :: ifformat_copy = 31
+  integer, parameter, public :: ifformat_promolecular = 32
+  integer, parameter, public :: ifformat_promolecular_fragment = 33
+  integer, parameter, public :: ifformat_as_hxx1 = 34
+  integer, parameter, public :: ifformat_as_hxx2 = 35
+  integer, parameter, public :: ifformat_as_hxx3 = 36
 
   ! free atomic polarizabilities from CRC handbook, 88th ed.
   real*8, parameter :: alpha_free(1:maxzat0) = (/  0.6668D0,  0.2051D0, 24.3300D0,  5.6000D0,& ! 1-4