Merge lp:~nickpapior/siesta/4.1-fdf-blocks into lp:siesta/4.1

[Nick Papior]

Implemented closing of the majority of commonly used blocks

Closing of blocks in the log when using fdf_b* calls was problematic
if the code did not read the last line of the block (i.e. this was
the only way to ensure the %endblock <label> in the fdf-*.log.

There are some problems with this and transiesta because the transiesta
codes are using fdf_bbackspace (which re-writes to the log file).
However, for Siesta the above should work.%

Since this change touches some core places I think it is better if you had a glance at it.

I have tested a few tests and it now writes out the %endblock in the log file.

It does, however, have some problems when using the fdf_bbackspace methods.

[Alberto Garcia (albertog)]

Approved.

review: Approve

=== modified file 'Src/Makefile'
--- Src/Makefile	2018-09-04 11:56:15 +0000
+++ Src/Makefile	2018-09-09 19:50:11 +0000
@@ -865,7 +865,7 @@
 m_ncdf_io.o: class_Sparsity.o m_io_s.o parallel.o precision.o
 m_ncdf_siesta.o: atm_types.o atmparams.o atomlist.o class_Sparsity.o files.o
 m_ncdf_siesta.o: kpoint_grid.o m_energies.o m_forces.o m_kinetic.o m_ncdf_io.o
-m_ncdf_siesta.o: m_ntm.o m_spin.o m_stress.o m_ts_electype.o m_ts_kpoints.o
+m_ncdf_siesta.o: m_spin.o m_stress.o m_ts_electype.o m_ts_kpoints.o
 m_ncdf_siesta.o: m_ts_options.o parallel.o precision.o radial.o siesta_geom.o
 m_ncdf_siesta.o: siesta_options.o sparse_matrices.o timestamp.o
 m_new_dm.o: alloc.o atomlist.o class_Data2D.o class_Fstack_Data1D.o
@@ -1238,9 +1238,9 @@
 siesta_analysis.o: write_subs.o writewave.o zmatrix.o
 siesta_cmlsubs.o: files.o m_uuid.o parallel.o siesta_cml.o timestamp.o
 siesta_cmlsubs.o: 
-siesta_dicts.o: atomlist.o files.o kpoint_grid.o m_energies.o m_forces.o
-siesta_dicts.o: m_mixing_scf.o m_steps.o m_stress.o precision.o siesta_geom.o
-siesta_dicts.o: siesta_options.o
+siesta_dicts.o: atomlist.o class_SpData1D.o files.o kpoint_grid.o m_energies.o
+siesta_dicts.o: m_forces.o m_mixing_scf.o m_steps.o m_stress.o precision.o
+siesta_dicts.o: siesta_geom.o siesta_options.o sparse_matrices.o
 siesta_end.o: alloc.o bands.o densematrix.o diag.o extrae_eventllist.o
 siesta_end.o: flook_siesta.o m_dscfcomm.o m_fixed.o m_io.o m_mixing_scf.o
 siesta_end.o: m_rhog.o m_wallclock.o meshdscf.o meshphi.o moremeshsubs.o
 
=== modified file 'Src/broyden_optim.F'
--- Src/broyden_optim.F	2016-01-25 16:00:16 +0000
+++ Src/broyden_optim.F	2018-09-09 19:50:11 +0000
@@ -187,6 +187,8 @@
             sxy = fdf_bvalues(pline,4)
             sxz = fdf_bvalues(pline,5)
             syz = fdf_bvalues(pline,6)
+            call fdf_bclose(bfdf)
+
             tstres(1,1) = - sxx * tp
             tstres(2,2) = - syy * tp
             tstres(3,3) = - szz * tp
 
=== modified file 'Src/cell_broyden_optim.F'
--- Src/cell_broyden_optim.F	2016-01-25 16:00:16 +0000
+++ Src/cell_broyden_optim.F	2018-09-09 19:50:11 +0000
@@ -180,6 +180,8 @@
             sxy = fdf_bvalues(pline,4)
             sxz = fdf_bvalues(pline,5)
             syz = fdf_bvalues(pline,6)
+            call fdf_bclose(bfdf)
+            
             target_stress(1,1) = - sxx * tp
             target_stress(2,2) = - syy * tp
             target_stress(3,3) = - szz * tp
 
=== modified file 'Src/cgvc.F'
--- Src/cgvc.F	2016-01-25 16:00:16 +0000
+++ Src/cgvc.F	2018-09-09 19:50:11 +0000
@@ -145,6 +145,8 @@
             sxy = fdf_bvalues(pline,4)
             sxz = fdf_bvalues(pline,5)
             syz = fdf_bvalues(pline,6)
+            call fdf_bclose(bfdf)
+            
             tstres(1,1) = - sxx * tp
             tstres(2,2) = - syy * tp
             tstres(3,3) = - szz * tp
 
=== modified file 'Src/cgvc_zmatrix.F'
--- Src/cgvc_zmatrix.F	2016-01-25 16:00:16 +0000
+++ Src/cgvc_zmatrix.F	2018-09-09 19:50:11 +0000
@@ -141,6 +141,8 @@
             sxy = fdf_bvalues(pline,4)
             sxz = fdf_bvalues(pline,5)
             syz = fdf_bvalues(pline,6)
+            call fdf_bclose(bfdf)
+            
             tstres(1,1) = - sxx * tp
             tstres(2,2) = - syy * tp
             tstres(3,3) = - szz * tp
 
=== modified file 'Src/chemical.f'
--- Src/chemical.f	2017-01-25 10:42:39 +0000
+++ Src/chemical.f	2018-09-09 19:50:11 +0000
@@ -115,22 +115,17 @@
       if ( present(silent) ) lsilent = silent
       if ( Node /= 0 ) lsilent = .true.
 
+      ! Default to 0
+      nsp = fdf_integer('Number_of_species',0)
+
       ! The most important thing to find is
       ! the block containing the species
-      found = fdf_block('Chemical_species_label',bfdf)
+      found = fdf_block('Chemical_species_label', bfdf)
       if (.not. found )
      $     call die("Block Chemical_species_label does not exist.")
 
-      ! Default to 0
-      nsp = fdf_integer('Number_of_species',0)
       if ( nsp == 0 ) then
-         ! try and guess the number of species
-         ns_read = 0
-         do while( fdf_bline(bfdf,pline) )
-            if ( fdf_bmatch(pline,'iin') ) then
-               ns_read = ns_read + 1
-            end if
-         end do
+         ns_read = fdf_block_linecount('Chemical_species_label', 'iin')
       else
          ns_read = nsp
       end if
@@ -144,9 +139,6 @@
       allocate(chemical_list%z(nsp))
       chemical_list%no_of_species = nsp
 
-      ! Rewind the block
-      call fdf_brewind(bfdf)
-
       ns_read = 0
       do while( fdf_bline(bfdf,pline) )
          if ( .not. fdf_bmatch(pline,'iin') ) cycle
 
=== modified file 'Src/coor.F'
--- Src/coor.F	2016-09-14 10:57:36 +0000
+++ Src/coor.F	2018-09-09 19:50:11 +0000
@@ -138,13 +138,8 @@
       if ( na == 0 ) then
        ! estimate the number of atoms by reading the block
        ! Currently we do not try to estimate from a Z-matrix
-       if (fdf_block('AtomicCoordinatesAndAtomicSpecies',bfdf)) then
-        na = 0
-        do while ( fdf_bline(bfdf,pline) )
-           ! each line should contain one atom
-           na = na + 1
-        end do
-       end if
+        na = fdf_block_linecount('AtomicCoordinatesAndAtomicSpecies',
+     &      'vvvi')
       end if
       if ( na == 0 ) then
          call die("Must specify number of atoms AND coordinates!")
@@ -170,6 +165,7 @@
          origin(1) = fdf_bvalues(pline,1)
          origin(2) = fdf_bvalues(pline,2)
          origin(3) = fdf_bvalues(pline,3)
+         call fdf_bclose(bfdf)
       else
          origin = 0._dp
       endif
@@ -195,7 +191,7 @@
       else
 
 C If Z matrix hasn't been found, read regular atomic coordinates
-        if (fdf_block('AtomicCoordinatesAndAtomicSpecies',bfdf)) then
+        if ( fdf_block('AtomicCoordinatesAndAtomicSpecies',bfdf) ) then
           do ia= 1, nua
             if (.not. fdf_bline(bfdf,pline))
      .        call die('coor: ERROR in ' //
@@ -209,7 +205,7 @@
             endif
 
           enddo
-
+          call fdf_bclose(bfdf)
         else
           call die("coor: You must specify the atomic coordinates")
         endif
 
=== modified file 'Src/fdf_extra.F90'
--- Src/fdf_extra.F90	2018-07-30 10:26:25 +0000
+++ Src/fdf_extra.F90	2018-09-09 19:50:11 +0000
@@ -73,13 +73,13 @@
           call fdf_blists(pline,il,i1,list(n+1:))
           if ( i1 + n > size(list) ) then
             print *,'Parsed line: ',trim(pline%line)
-            call die('fdf_extra: number of elements in block list &
+            call die('fdf_brange: number of elements in block list &
                 &is too many to fit the maximal range of the &
                 &list. Please correct.')
           end if
           if ( i1 == 0 ) then
             print *,'Parsed line: ',trim(pline%line)
-            call die('fdf_extra: a block list with zero elements is not &
+            call die('fdf_brange: a block list with zero elements is not &
                 &allowed, please correct input.')
           end if
 
@@ -102,14 +102,14 @@
       g = fdf_bnames(pline,2)
       if ( .not. leqi(g,'from') ) then
         print *,'Parsed line: ',trim(pline%line)
-        call die('fdf_extra: error in range block: &
+        call die('fdf_brange: error in range block: &
             &from <int> to/plus/minus <int> is ill formatted')
       end if
 
       g = fdf_bnames(pline,3)
       if ( fdf_bnintegers(pline) < 2 ) then
         print *,'Parsed line: ',trim(pline%line)
-        call die('fdf_extra: error in range block &
+        call die('fdf_brange: error in range block &
             &from <int> to/plus/minus <int> is ill formatted')
       end if
 
@@ -135,7 +135,7 @@
         i2 = i1 - i2 + 1
       else
         print *,'Parsed line: ',trim(pline%line)
-        call die('fdf_extra: unrecognized designator of ending range, &
+        call die('fdf_brange: unrecognized designator of ending range, &
             &[to, plus, minus] accepted.')
       end if
 
@@ -153,7 +153,7 @@
           (i1 > i2 .and. step > 0) ) then
         print *,'Parsed line: ',trim(pline%line)
         print *,i1,i2,step
-        call die('fdf_extra: block range is not consecutive')
+        call die('fdf_brange: block range is not consecutive')
       end if
 
       ! Create list
@@ -166,14 +166,14 @@
     else
 
       print *,'Parsed line: ',trim(pline%line)
-      call die('fdf_extra: error in range block, input not recognized')
+      call die('fdf_brange: error in range block, input not recognized')
 
     end if
 
     do j = 1 , n
       if ( list(j) < low .or. high < list(j) ) then
         print *,'Parsed line: ',trim(pline%line)
-        call die('fdf_extra: error in range block. Input is beyond range')
+        call die('fdf_brange: error in range block. Input is beyond range')
       end if
     end do
 
@@ -240,16 +240,16 @@
     n_r = 0
     if ( allocated(rgns) ) deallocate(rgns)
 
-    ! If the block does not exist, simply return
-    if ( .not. fdf_block(bName,bfdf) ) return
-
-    ! the initial number of regions
-    il = 0
-    do while ( fdf_bnext(bfdf,pline) ) 
-      il = il + 1
-    end do
+    ! Get number of regions
+    il = fdf_block_linecount(bName)
+    if ( il == 0 ) return
+
+    if ( .not. fdf_block(bName,bfdf) ) then
+      call die('fdf_bregions: failed implementation.')
+    end if
+
+    ! allocate
     allocate(rlist(il))
-    call fdf_brewind(bfdf)
 
     ! first count number of differently named regions
     do while ( fdf_bnext(bfdf,pline) ) 
@@ -302,7 +302,7 @@
         call fdf_brange(pline,r1,1,n)
         if ( r1%n == 0 ) then
           print *,'Region: ',trim(g)
-          call die('fdf_extra: Could not read in anything in region!')
+          call die('fdf_bregions: Could not read in anything in region!')
         end if
         call rgn_union(rgns(il),r1,rgns(il))
         rgns(il)%name = trim(g)
@@ -312,7 +312,7 @@
         call fdf_brange(pline,r1,1,n)
         if ( r1%n == 0 ) then
           print *,'Region: ',trim(g)
-          call die('fdf_extra: Could not read in anything in region!')
+          call die('fdf_bregions: Could not read in anything in region!')
         end if
         call rgn_union(rgns(ic),r1,rgns(ic))
         rgns(ic)%name = trim(g)
 
=== modified file 'Src/get_target_stress.f'
--- Src/get_target_stress.f	2016-01-25 16:00:16 +0000
+++ Src/get_target_stress.f	2018-09-09 19:50:11 +0000
@@ -42,6 +42,8 @@
             sxy = fdf_bvalues(pline,4)
             sxz = fdf_bvalues(pline,5)
             syz = fdf_bvalues(pline,6)
+            call fdf_bclose(bfdf)
+            
             tstres(1,1) = - sxx * tp
             tstres(2,2) = - syy * tp
             tstres(3,3) = - szz * tp
 
=== modified file 'Src/kgrid.F'
--- Src/kgrid.F	2016-01-25 16:00:16 +0000
+++ Src/kgrid.F	2018-09-09 19:50:11 +0000
@@ -63,7 +63,7 @@
       parameter (tiny   = 1.d-12)
 
 C Find out if there is spiral arrangement of spins
-      spiral = fdf_block('SpinSpiral',bfdf)
+      spiral = fdf_isblock('SpinSpiral')
 
 C Find total number of points (determinant of kscell)
       nktot = abs( kscell(1,1) * kscell(2,2) * kscell(3,3) +
 
=== modified file 'Src/kgridinit.F'
--- Src/kgridinit.F	2017-10-10 19:47:15 +0000
+++ Src/kgridinit.F	2018-09-09 19:50:11 +0000
@@ -92,7 +92,7 @@
 
       genlogic = .false.
 C Find out if there is spiral arrangement of spins 
-      spiral = fdf_block('SpinSpiral',bfdf)
+      spiral = fdf_isblock('SpinSpiral')
 
 C Find total number of points (determinant of kscell)
       nktot = abs( kscell(1,1) * kscell(2,2) * kscell(3,3) +
@@ -118,6 +118,7 @@
                displ(i) = 0._dp
             end if
           enddo
+          call fdf_bclose(bfdf)
           nktot = abs( kscell(1,1) * kscell(2,2) * kscell(3,3) +
      .                 kscell(2,1) * kscell(3,2) * kscell(1,3) +
      .                 kscell(3,1) * kscell(1,2) * kscell(2,3) -
 
=== modified file 'Src/kpoint_grid.F90'
--- Src/kpoint_grid.F90	2018-04-14 23:14:53 +0000
+++ Src/kpoint_grid.F90	2018-09-09 19:50:11 +0000
@@ -145,6 +145,7 @@
               kdispl(i) = 0._dp
             end if
          enddo
+         call fdf_bclose(bfdf)
          firm_displ = .true.
 
       else
 
=== modified file 'Src/kpoint_pdos.F90'
--- Src/kpoint_pdos.F90	2018-04-14 23:14:53 +0000
+++ Src/kpoint_pdos.F90	2018-09-09 19:50:11 +0000
@@ -159,7 +159,8 @@
             else
                kdispl(i) = 0._dp
             end if
-         enddo
+          enddo
+          call fdf_bclose(bfdf)
          firm_displ = .true.
 
       else
 
=== modified file 'Src/ksvinit.F'
--- Src/ksvinit.F	2016-01-25 16:00:16 +0000
+++ Src/ksvinit.F	2018-09-09 19:50:11 +0000
@@ -229,7 +229,7 @@
           endif
 
         enddo
-
+        call fdf_bclose(bfdf)
       endif
 
       end subroutine repol
 
=== modified file 'Src/ldau_specs.f'
--- Src/ldau_specs.f	2016-11-25 18:04:25 +0000
+++ Src/ldau_specs.f	2018-09-09 19:50:11 +0000
@@ -428,7 +428,7 @@
       nprojsldau(:) = 0
 
 !     Read the LDAU.proj block
-      if (.not. fdf_block('LDAU.proj',bfdf)) RETURN
+      if (.not. fdf_block('LDAU.proj', bfdf)) RETURN
 
       ! Add citation
       if ( IONode ) then
 
=== modified file 'Src/local_DOS.F'
--- Src/local_DOS.F	2018-06-25 11:33:42 +0000
+++ Src/local_DOS.F	2018-09-09 19:50:11 +0000
@@ -70,6 +70,7 @@
         factor = fdf_convfac( fdf_bnames(pline,1), 'Ry' )
         e1 = fdf_bvalues(pline,1)*factor
         e2 = fdf_bvalues(pline,2)*factor
+        call fdf_bclose(bfdf)
 
         !Find the density matrix for states between e1 and e2
         if ((isolve .eq. SOLVE_DIAGON) .or.
 
=== modified file 'Src/m_efield.F'
--- Src/m_efield.F	2018-01-26 13:17:30 +0000
+++ Src/m_efield.F	2018-09-09 19:50:11 +0000
@@ -52,7 +52,6 @@
 
       real(dp),   intent(out)     :: input_field(3)
 
-      logical                    :: found
       type(block_fdf)            :: bfdf
       type(parsed_line), pointer :: pline => null()
 
@@ -63,10 +62,8 @@
 
       input_field(1:3) = 0.0_dp
 
-      found = fdf_block('ExternalElectricField',bfdf)
-      if (.not. found ) RETURN
-      loop: DO
-        if (.not. fdf_bline(bfdf,pline)) exit loop
+      if (.not. fdf_block('ExternalElectricField',bfdf) ) RETURN
+      do while ( fdf_bline(bfdf,pline) )
         if (.not. fdf_bmatch(pline,"vvvn"))
      $     call die("Wrong format in ElectricField block")
         eunits = fdf_bnames(pline,1)
@@ -74,7 +71,7 @@
         do ix = 1,3
            input_field(ix) = fdf_bvalues(pline,ix) * cfactor
         enddo
-      enddo loop
+      enddo
 
       end subroutine get_user_specified_field
 !----------------------------------------------------------------
 
=== modified file 'Src/m_new_dm.F90'
--- Src/m_new_dm.F90	2018-06-28 13:53:17 +0000
+++ Src/m_new_dm.F90	2018-09-09 19:50:11 +0000
@@ -981,7 +981,7 @@
       
       ! Read the data from the block and then we populate DM
       na = 0
-      do while( fdf_bline(bfdf,pline) .and. na < na_u )
+      do while( fdf_bline(bfdf,pline) )
          
          ! Read number of names, integers and reals on this line
          ni = fdf_bnintegers(pline)
 
=== modified file 'Src/m_spin.F90'
--- Src/m_spin.F90	2018-02-27 14:03:49 +0000
+++ Src/m_spin.F90	2018-09-09 19:50:11 +0000
@@ -338,7 +338,7 @@
   subroutine init_spiral( ucell )
     use fdf, only : fdf_get, leqi
     use fdf, only: block_fdf, parsed_line
-    use fdf, only: fdf_block, fdf_bline
+    use fdf, only: fdf_block, fdf_bline, fdf_bclose
     use fdf, only: fdf_bnames, fdf_bvalues
     use units, only: Pi
 
@@ -381,6 +381,9 @@
        call die('init_spiral: ERROR: ReciprocalCoordinates must be' // &
             ' ''Cubic'' or ''ReciprocalLatticeVectors''')
     end if
+
+    call fdf_bclose(bfdf)
+
   end subroutine init_spiral
 
   function fname_spin(nspin,ispin,slabel,suffix,basename) result(fname)
 
=== modified file 'Src/m_target_stress.F'
--- Src/m_target_stress.F	2016-10-31 21:00:50 +0000
+++ Src/m_target_stress.F	2018-09-09 19:50:11 +0000
@@ -73,6 +73,8 @@
             sxy = fdf_bvalues(pline,4)
             sxz = fdf_bvalues(pline,5)
             syz = fdf_bvalues(pline,6)
+            call fdf_bclose(bfdf)
+            
             target_stress(1,1) = - sxx * tp
             target_stress(2,2) = - syy * tp
             target_stress(3,3) = - szz * tp
 
=== modified file 'Src/m_ts_io_ctype.f90'
--- Src/m_ts_io_ctype.f90	2018-05-04 10:51:06 +0000
+++ Src/m_ts_io_ctype.f90	2018-09-09 19:50:11 +0000
@@ -98,26 +98,11 @@
     character(len=*), intent(in) :: suffix
     integer :: n
 
-    ! prepare to read in the data...
-    type(block_fdf) :: bfdf
-    type(parsed_line), pointer :: pline => null()
-    
-    logical :: found
-
-    n = 0
     if ( len_trim(suffix) == 0 ) then
-       found = fdf_block(trim(prefix)//'.Contours',bfdf)
+      n = fdf_block_linecount(trim(prefix)//'.Contours', 'l')
     else
-       found = fdf_block(trim(prefix)//'.Contours.'//trim(suffix),bfdf)
+      n = fdf_block_linecount(trim(prefix)//'.Contours'//trim(suffix), 'l')
     end if
-    if ( .not. found ) return
-
-    ! first count the number of electrodes
-    n = 0
-    do while ( fdf_bline(bfdf,pline) )
-       if ( fdf_bnnames(pline) == 0 ) cycle
-       n = n + 1 
-    end do
 
   end function fdf_nc_iotype
 
@@ -143,15 +128,16 @@
     end if
     if ( .not. found ) return
 
-    ! first count the number of electrodes
+    ! Find the name of the contour
     n = 0
     do while ( fdf_bline(bfdf,pline) )
-       if ( fdf_bnnames(pline) == 0 ) cycle
-       n = n + 1 
-       if ( n == i ) then
-          name = fdf_bnames(pline,1)
-          return
-       end if
+      if ( fdf_bnnames(pline) == 0 ) cycle
+      n = n + 1 
+      if ( n == i ) then
+        name = fdf_bnames(pline,1)
+        call fdf_bclose(bfdf)
+        return
+      end if
     end do
 
   end function fdf_name_c_iotype
@@ -226,7 +212,6 @@
     type(block_fdf), optional :: bfdf
     logical :: exist
 
-    type(block_fdf) :: bfdf_tmp
     character(len=c_N) :: g
 
     ! if the block does not exist, return
@@ -242,7 +227,7 @@
           bfdf%label = trim(g)
        end if
     else
-       exist = fdf_block(trim(g),bfdf_tmp)
+       exist = fdf_isblock(trim(g))
     end if
 
   end function ts_exists_contour_block
 
=== modified file 'Src/m_ts_kpoints.F90'
--- Src/m_ts_kpoints.F90	2017-11-23 14:17:27 +0000
+++ Src/m_ts_kpoints.F90	2018-09-09 19:50:11 +0000
@@ -131,6 +131,7 @@
                   'kgrid_Monkhorst_Pack block' )
           endif
        enddo
+       call fdf_bclose(bfdf)
        ts_firm_displ = .true.
        
     else
 
=== modified file 'Src/meshsubs.F'
--- Src/meshsubs.F	2017-07-01 00:59:38 +0000
+++ Src/meshsubs.F	2018-09-09 19:50:11 +0000
@@ -334,6 +334,7 @@
         else
            call die('initmesh: ERROR in MeshSizes block')
         endif
+        call fdf_bclose(bfdf)
         RealCutoff = huge(1.0_dp)
         call chkgmx( k0, rcell, ntm, RealCutoff )
         if (RealCutoff < G2max) then
 
=== modified file 'Src/metaforce.F'
--- Src/metaforce.F	2016-01-25 16:00:16 +0000
+++ Src/metaforce.F	2018-09-09 19:50:11 +0000
@@ -11,13 +11,12 @@
 
       implicit none
 
-      integer,           save :: maxGauss = 100
       integer,           save :: nGauss = 0
-      logical,           save :: lMetaForce
-      integer,  pointer, save :: nGaussPtr(:,:)
-      real(dp), pointer, save :: GaussK(:)
-      real(dp), pointer, save :: GaussR0(:)
-      real(dp), pointer, save :: GaussZeta(:)
+      logical,           save :: lMetaForce = .false.
+      integer,  pointer, save :: nGaussPtr(:,:) => null()
+      real(dp), pointer, save :: GaussK(:) => null()
+      real(dp), pointer, save :: GaussR0(:) => null()
+      real(dp), pointer, save :: GaussZeta(:) => null()
 
       CONTAINS
 
@@ -41,34 +40,42 @@
       type(block_fdf)            :: bfdf
       type(parsed_line), pointer :: pline
 
-C Find if there are any Gaussians specified
-      lMetaForce = fdf_defined('MetaForce')
-      if (lMetaForce) then
+C     Find if there are any Gaussians specified
+      nGauss = fdf_block_linecount('MetaForce', 'iirrr')
+      lMetaForce = nGauss > 0
+      if ( lMetaForce ) then
 
 C Allocate arrays for Gaussian data
-        call re_alloc( nGaussPtr, 1, 2, 1, maxGauss, 'nGaussPtr',
+        call re_alloc( nGaussPtr, 1, 2, 1, nGauss, 'nGaussPtr',
      &                 'initmeta' )
-        call re_alloc( GaussK, 1, maxGauss, 'GaussK', 'initmeta' )
-        call re_alloc( GaussR0, 1, maxGauss, 'GaussR0', 'initmeta' )
-        call re_alloc( GaussZeta, 1, maxGauss, 'GaussZeta', 'initmeta' )
+        call re_alloc( GaussK, 1, nGauss, 'GaussK', 'initmeta' )
+        call re_alloc( GaussR0, 1, nGauss, 'GaussR0', 'initmeta' )
+        call re_alloc( GaussZeta, 1, nGauss, 'GaussZeta', 'initmeta' )
 
 C Read Gaussians from block
-        lMetaForce = fdf_block('MetaForce',bfdf)
+        lMetaForce = fdf_block('MetaForce', bfdf)
 
-        do k= 1, maxGauss
+        do while ( fdf_bline(bfdf, pline) )
 C Read and parse data line
-          if (.not. fdf_bline(bfdf,pline))
-     .      call die('initmeta: ERROR in MetaForce block')
           ni = fdf_bnintegers(pline)
           nr = fdf_bnreals(pline)
 
 C Check that correct info is given
           if ((ni .ge. 2) .and. (nr .ge. 3)) then
             nGauss = nGauss + 1
+            ! Atom 1
             nGaussPtr(1,nGauss) = abs(fdf_bintegers(pline,1))
+            ! Atom 2
             nGaussPtr(2,nGauss) = abs(fdf_bintegers(pline,2))
+            ! Gaussian pre-factor for the force
+            ! Unit: Ry
             GaussK(nGauss) = fdf_breals(pline,1)
+            ! Distance in r - r_0, where r is the distance between
+            ! atom i and j
+            ! Unit: Ang
             GaussR0(nGauss) = fdf_breals(pline,2)/0.529177_dp
+            ! Gaussian prefactor
+            ! Unit: 1/Bohr**2
             GaussZeta(nGauss) = fdf_breals(pline,3)
           endif
         enddo
 
=== modified file 'Src/projected_DOS.F'
--- Src/projected_DOS.F	2018-06-25 11:33:42 +0000
+++ Src/projected_DOS.F	2018-09-09 19:50:11 +0000
@@ -22,17 +22,15 @@
       subroutine init_projected_DOS()
 
       USE siesta_options
-      use fdf,         only: fdf_block, block_fdf
+      use fdf,         only: fdf_block, fdf_isblock
       use Kpoint_pdos
       use parallel,    only: IOnode, Nodes
       use siesta_geom, only: ucell
       use m_spin,      only: spin
 
-      type(block_fdf) :: bfdf
-      
 !-------------------------------------------------------------------------BEGIN
 !     Compute the projected density of states
-      do_pdos = fdf_block('ProjectedDensityOfStates', bfdf)
+      do_pdos = fdf_isblock('ProjectedDensityOfStates')
       if ( .not. do_pdos ) return
          
       if (isolve.ne.SOLVE_DIAGON) then
@@ -108,6 +106,7 @@
      $        'siesta: e1, e2, sigma, nhist: ',
      $        e1/eV,' eV',e2/eV,' eV',sigma/eV,' eV', nhist
       end if
+      call fdf_bclose(bfdf) 
 
 ! If the k points have been set specifically for the PDOS then use this set
       if (different_pdos_grid) then
 
=== modified file 'Src/read_xc_info.F'
--- Src/read_xc_info.F	2016-01-25 16:00:16 +0000
+++ Src/read_xc_info.F	2018-09-09 19:50:11 +0000
@@ -78,6 +78,7 @@
               XCweightC(n) = 1.0_dp
             endif
           enddo
+          call fdf_bclose(bfdf)
         else
           nXCfunc = 1 
           XCfunc(1) = fdf_string('xc.functional','LDA')
 
=== modified file 'Src/readsp.F'
--- Src/readsp.F	2016-01-25 16:00:16 +0000
+++ Src/readsp.F	2018-09-09 19:50:11 +0000
@@ -84,6 +84,7 @@
           call die('redata: ERROR: ReciprocalCoordinates must be' //
      .             ' ''Cubic'' or ''ReciprocalLatticeVectors''')
         endif
+        call fdf_bclose(bfdf)
       else
         spiral = .false.
         q(1) = 0.d0
 
=== modified file 'Src/redcel.F'
--- Src/redcel.F	2016-01-25 16:00:16 +0000
+++ Src/redcel.F	2018-09-09 19:50:11 +0000
@@ -43,8 +43,8 @@
 
 C Lattice vectors
 
-      if ( fdf_block('LatticeParameters',bfdf) .and.
-     .     fdf_block('LatticeVectors',bfdf) ) then
+      if ( fdf_isblock('LatticeParameters') .and.
+     .    fdf_isblock('LatticeVectors') ) then
          call die('redcel: ERROR: Lattice vectors doubly ' //
      $       'specified: by LatticeVectors and by LatticeParameters.')
       endif
@@ -97,6 +97,7 @@
           cell(i,i) = 1.0_dp
         enddo
       endif
+      call fdf_bclose(bfdf)
 
       ! Note that cell will be zero if alat is not specified
       ! This is not very intuitive
@@ -117,6 +118,7 @@
           mscell(2,i) = fdf_bvalues(pline,2)
           mscell(3,i) = fdf_bvalues(pline,3)
         enddo
+        call fdf_bclose(bfdf)
       else
         do i = 1,3
           do j = 1,3
 
=== modified file 'Src/reoptical.F'
--- Src/reoptical.F	2016-03-21 22:47:29 +0000
+++ Src/reoptical.F	2018-09-09 19:50:11 +0000
@@ -138,6 +138,7 @@
           nmeshk(1) = fdf_bintegers(pline,1)
           nmeshk(2) = fdf_bintegers(pline,2)
           nmeshk(3) = fdf_bintegers(pline,3)
+          call fdf_bclose(bfdf)
         endif
 
 C Read whether this is a polarised, unpolarised or polycrystal calc
 
=== modified file 'Src/zm_broyden_optim.F'
--- Src/zm_broyden_optim.F	2016-03-01 14:24:48 +0000
+++ Src/zm_broyden_optim.F	2018-09-09 19:50:11 +0000
@@ -175,7 +175,8 @@
             sxy = fdf_bvalues(pline,4)
             sxz = fdf_bvalues(pline,5)
             syz = fdf_bvalues(pline,6)
-
+            call fdf_bclose(bfdf)
+            
             tstres(1,1) = - sxx * tp
             tstres(2,2) = - syy * tp
             tstres(3,3) = - szz * tp
 
=== modified file 'version.info'
--- version.info	2018-09-08 04:18:59 +0000
+++ version.info	2018-09-09 19:50:11 +0000
@@ -1,1 +1,1 @@
-siesta-4.1--983
+siesta-4.1--983--fdf-close-1