Merge lp:~nickpapior/siesta/4.1-ncdf-md into lp:siesta/4.1

[Nick Papior]

Fixed NCDF write-out for MD and FC calculations

Previously the NCDF file-write-out would crash when performing
MD/FC calculations because the sparse arrays does not really make sense
to store.

Now the file content is based on the type of run, single/MD/FC.
A relaxation is equivalent to an MD (currently).

One non-CDF related change is the change of MD.NumCGSteps to
MD.Steps keyword. The NumCGSteps is deprecated since many more
options have entered and CG was not necessarily the correct instance.
Same applies to MD.MaxCGDispl -> MD.MaxDispl.

Secondly, the MD.Steps keyword is now also the default for
MD.FinalTimeStep which means that a single keyword determines the number
of steps runned.

The old flags are still used as the default values.

Basically this is one big fix for the NCDF routines in Siesta.

However, there is one change that "affects" siesta.

MD.NumCGSteps is renamed to MD.Steps (please suggest if another keyword is better).
The procedure is this:

i = fdf_get('MD.NumCGSteps', 0)
i = fdf_get('MD.Steps', i)

I think the NumCGsteps is far too confusing in 4.1 (even in 4.0) and it should have been updated previously.
The same changes have been applied to MD.MaxCGDispl (now MD.MaxDispl).

Secondly, since MD.Steps is more generic it is now also the default in MD.FinalTimeStep

j = fdf_get('MD.FinalTimeStep', i) ! i from above.

Hence one single keyword may be used generically throughout.

[Alberto Garcia (albertog)]

OK to the introduction of MD.Steps et al.

This is a minor point:

If someone uses MD.TypeOfRun <some MD option> but forgets to set MD.FinalTimeStep, a fossil MD.NumCGSteps might set it silently. I would put in code to complain if MD.Steps is not set either.

Also:

* In the ncdf module, the KB projectors might (in the trunk) now have an extra "j" index (when offsite SO is used). I put them in the .ion file, but I forgot to put them in the vanilla cdf files. Please take this into account when you port this patch into the trunk.

* Correct the indentation of

! 5 = -z

in final_H_f_stress.

review: Approve

[lp:~nickpapior/siesta/4.1-ncdf-md]

939. By Nick Papior on 2018-07-09

Fixed FinalTimeStep option reads and indentation

Thanks Alberto

940. By Nick Papior on 2018-07-09

Added pjnl_j to the ncdf basis

[Nick Papior (nickpapior)]

Ok, I will insert an fdf_defined before FinalTimestep.

Add 0's to the output of the KB projectors.

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2018-06-26 08:27:41 +0000
+++ Docs/siesta.tex	2018-07-09 07:59:46 +0000
@@ -5796,7 +5796,7 @@
   if the lattice changes are moderate).
 
   By default this will be \fdffalse\ when performing single-point calculations
-  (\fdf{MD.TypeOfRun:CG} with \fdf{MD.NumCGsteps:1} or below). Otherwise
+  (\fdf{MD.TypeOfRun:CG} with \fdf{MD.Steps:1} or below). Otherwise
   this will default to \fdftrue.
 
   \note one may try to set this to \fdffalse\ for better performance,
@@ -9364,15 +9364,20 @@
 
 \end{fdfentry}
 
-\begin{fdfentry}{MD.NumCGsteps}[integer]<0>
+\begin{fdfentry}{MD.Steps}[integer]<0>
+  \fdfdepend{MD.NumCGsteps}
   
-  Maximum number of conjugate gradient (or Broyden) minimization moves
+  Maximum number of steps in a minimization routine
   (the minimization will stop if tolerance is reached before; see
   \fdf{MD.MaxForceTol} below).
 
+  \note The old flag \fdf*{MD.NumCGsteps} will remain for historical
+  reasons.
+
 \end{fdfentry}
 
-\begin{fdfentry}{MD.MaxCGDispl}[length]<$0.2\,\mathrm{Bohr}$>
+\begin{fdfentry}{MD.MaxDispl}[length]<$0.2\,\mathrm{Bohr}$>
+  \fdfdepend{MD.MaxCGDispl}
   
   Maximum atomic displacements in an optimization move.
 
@@ -9383,6 +9388,9 @@
 
   Note that there are Zmatrix-specific options that override this option.
 
+  \note The old flag \fdf*{MD.MaxCGDispl} will remain for historical
+  reasons.
+
 \end{fdfentry}
 
 \begin{fdfentry}{MD.PreconditionVariableCell}[length]<$5\,\mathrm{Ang}$>
@@ -9659,7 +9667,7 @@
 
 \end{fdfentry}
 
-\begin{fdfentry}{MD.FinalTimeStep}[integer]<$1$>
+\begin{fdfentry}{MD.FinalTimeStep}[integer]<\fdfvalue{MD.Steps}>
 
   Final time step of the MD simulation.
 

=== modified file 'Src/Makefile'
--- Src/Makefile	2018-06-28 13:53:17 +0000
+++ Src/Makefile	2018-07-09 07:59:46 +0000
@@ -1293,14 +1293,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o m_energies.o m_eo.o m_mixing.o
-state_init.o: m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o m_pivot_methods.o
-state_init.o: m_rmaxh.o m_sparse.o m_sparsity_handling.o m_spin.o m_steps.o
-state_init.o: m_supercell.o m_test_io.o m_ts_charge.o m_ts_electype.o
-state_init.o: m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o m_ts_options.o
-state_init.o: m_ts_sparse.o m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
-state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
-state_init.o: siesta_options.o sparse_matrices.o sys.o units.o write_subs.o
-state_init.o: zmatrix.o
+state_init.o: m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o m_new_dm.o m_os.o
+state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
+state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
+state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
+state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
+state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
+state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
+state_init.o: sys.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Src/final_H_f_stress.F'
--- Src/final_H_f_stress.F	2018-06-20 10:34:33 +0000
+++ Src/final_H_f_stress.F	2018-07-09 07:59:46 +0000
@@ -17,7 +17,7 @@
 #ifdef NCDF_4
       use siesta_options, only: write_cdf
 #endif
-      use siesta_options, only: g2cut, savehs, temp
+      use siesta_options, only: g2cut, savehs, temp, nmove
       use siesta_options, only: recompute_H_after_scf
       use sparse_matrices, only: numh, listh, listhptr
       use sparse_matrices, only: H, S, Dscf, Escf, maxnh, xijo
@@ -297,10 +297,23 @@
 #ifdef CDF
 #ifdef NCDF_4
       if ( write_cdf ) then
-         dic_save = ('fa'.kv.1)//('stress'.kv.1)
-         if ( savehs .or. TS_HS_save ) dic_save = dic_save //('H'.kv.1)
-         call cdf_save_state(trim(slabel)//'.nc',dic_save)
-         call delete(dic_save)
+        if ( idyn == 6 ) then
+          call cdf_save_fc(trim(slabel)//'.nc', istep)
+          if ( istep == 0 ) then
+            dic_save = ('fa'.kv.1)//('stress'.kv.1)
+            if ( savehs .or. TS_HS_save )
+     &          dic_save = dic_save//('H'.kv.1)
+            call cdf_save_state(trim(slabel)//'.nc',dic_save)
+            call delete(dic_save)
+          end if
+        else if ( idyn /=0 .or. nmove /= 0 ) then
+          call cdf_save_md(trim(slabel)//'.nc')
+        else
+          dic_save = ('fa'.kv.1)//('stress'.kv.1)
+          if ( savehs .or. TS_HS_save ) dic_save = dic_save //('H'.kv.1)
+          call cdf_save_state(trim(slabel)//'.nc',dic_save)
+          call delete(dic_save)
+        end if
       end if
 #endif
 #endif

=== modified file 'Src/m_ncdf_siesta.F90'
--- Src/m_ncdf_siesta.F90	2018-06-12 22:59:08 +0000
+++ Src/m_ncdf_siesta.F90	2018-07-09 07:59:46 +0000
@@ -29,6 +29,7 @@
   public :: cdf_save_settings
   public :: cdf_save_basis
   public :: cdf_save_md
+  public :: cdf_save_fc
   public :: cdf_save_state
   public :: cdf_save_grid
 
@@ -38,7 +39,7 @@
 
 #ifdef NCDF_4
 
-  subroutine cdf_init_file(fname,is_md)
+  subroutine cdf_init_file(fname, is_md, is_fc)
     use fdf, only : fdf_get, leqi
     use class_Sparsity
     use files, only : slabel
@@ -57,13 +58,15 @@
     use siesta_options, only: savebader
     use siesta_options, only: save_initial_charge_density
     use siesta_options, only: fixspin, total_spin
+    use siesta_options, only: ia1, ia2, dx ! FC information
     use m_timestamp, only: datestring
     use m_ts_electype, elec_name => name
     use m_ts_options, only : Volt, N_Elec, Elecs
     use m_ts_options, only : TS_HS_Save
 
     character(len=*), intent(in) :: fname
-    logical, intent(in) :: is_md
+    logical, intent(in), optional :: is_md
+    logical, intent(in), optional :: is_fc
 
     ! Local variables
     type(hNCDF) :: ncdf, grp, grp2
@@ -72,10 +75,16 @@
     integer :: n_nzs, tmp, i, chks(3), iEl
     integer, allocatable :: ibuf(:)
     logical :: exists
+    logical :: lis_md, lis_fc
 #ifdef MPI
     integer :: MPIerror
 #endif
 
+    lis_md = .false.
+    if ( present(is_md) ) lis_md = is_md
+    lis_fc = .false.
+    if ( present(is_fc) ) lis_fc = is_fc
+
     ! We always re-write the file...
     call ncdf_create(ncdf,fname, &
          mode=ior(NF90_WRITE,NF90_NETCDF4), overwrite=.true.)
@@ -96,112 +105,120 @@
     d = d// ('DIMone'.kv. 1)
     call ncdf_crt(ncdf,d)
     call delete(d)
-    
-    ! Create all necessary containers...
-    dic = ('info'.kv.'Number of supercells in each unit-cell direction') 
-    call ncdf_def_var(ncdf,'nsc',NF90_INT,(/'xyz'/), &
-         compress_lvl=0, atts=dic)
 
     dic = dic//('info'.kv.'Last orbital of equivalent atom')
     call ncdf_def_var(ncdf,'lasto',NF90_INT,(/'na_u'/), &
-         compress_lvl=0, atts=dic)
+        compress_lvl=0, atts=dic)
+    call ncdf_put_var(ncdf,'lasto',lasto(1:na_u))
     
     dic = dic//('info'.kv.'Total charge')
     call ncdf_def_var(ncdf,'Qtot',NF90_DOUBLE,(/'one'/), &
          compress_lvl=0, atts=dic)
+    call ncdf_put_var(ncdf,'Qtot',Qtot)
 
-    dic = dic//('info'.kv.'Fermi level')//('unit'.kv.'Ry')
     if ( fixspin ) then
-      call ncdf_def_var(ncdf,'Ef',NF90_DOUBLE,(/'spin'/), &
-          compress_lvl=0, atts=dic)
-
-      call delete(dic)
       dic = dic//('info'.kv.'Total spin')
       call ncdf_def_var(ncdf,'Qspin',NF90_DOUBLE,(/'one'/), &
           compress_lvl=0, atts=dic)
-    else
-      call ncdf_def_var(ncdf,'Ef',NF90_DOUBLE,(/'one'/), &
-          compress_lvl=0, atts=dic)
-    end if
-    
-    dic = dic//('info'.kv.'Atomic coordinates')//('unit'.kv.'Bohr')
-    call ncdf_def_var(ncdf,'xa',NF90_DOUBLE,(/'xyz ','na_u'/), &
-         compress_lvl=0, atts=dic)
-    
-    dic = dic//('info'.kv.'Unit cell')
-    call ncdf_def_var(ncdf,'cell',NF90_DOUBLE,(/'xyz','xyz'/), &
-         compress_lvl=0, atts=dic)
-    
-    dic = dic//('info'.kv.'Atomic forces')//('unit'.kv.'Ry/Bohr')
-    call ncdf_def_var(ncdf,'fa',NF90_DOUBLE,(/'xyz ','na_u'/), &
-         compress_lvl=0, atts=dic)
-    
-    dic = dic//('info'.kv.'Cell stress')//('unit'.kv.'Ry/Bohr**3')
-    call ncdf_def_var(ncdf,'stress',NF90_DOUBLE,(/'xyz','xyz'/), &
-         compress_lvl=0, atts=dic)
+      call ncdf_put_var(ncdf,'Qspin',total_spin)
+    end if
+    
+    if ( .not. lis_MD ) then
+    
+      ! Create all necessary containers...
+      dic = dic//('info'.kv.'Number of supercells in each unit-cell direction') 
+      call ncdf_def_var(ncdf,'nsc',NF90_INT,(/'xyz'/), &
+          compress_lvl=0, atts=dic)
+      
+      dic = dic//('info'.kv.'Fermi level')//('unit'.kv.'Ry')
+      if ( fixspin ) then
+        call ncdf_def_var(ncdf,'Ef',NF90_DOUBLE,(/'spin'/), &
+            compress_lvl=0, atts=dic)
+      else
+        call ncdf_def_var(ncdf,'Ef',NF90_DOUBLE,(/'one'/), &
+            compress_lvl=0, atts=dic)
+      end if
+
+      dic = dic//('info'.kv.'Atomic coordinates')//('unit'.kv.'Bohr')
+      call ncdf_def_var(ncdf,'xa',NF90_DOUBLE,(/'xyz ','na_u'/), &
+          compress_lvl=0, atts=dic)
+
+      dic = dic//('info'.kv.'Unit cell')
+      call ncdf_def_var(ncdf,'cell',NF90_DOUBLE,(/'xyz','xyz'/), &
+          compress_lvl=0, atts=dic)
+
+      dic = dic//('info'.kv.'Atomic forces')//('unit'.kv.'Ry/Bohr')
+      call ncdf_def_var(ncdf,'fa',NF90_DOUBLE,(/'xyz ','na_u'/), &
+          compress_lvl=0, atts=dic)
+
+      dic = dic//('info'.kv.'Cell stress')//('unit'.kv.'Ry/Bohr**3')
+      call ncdf_def_var(ncdf,'stress',NF90_DOUBLE,(/'xyz','xyz'/), &
+          compress_lvl=0, atts=dic)
+    end if
 
     call delete(dic)
 
-    
-    ! Create matrix group
-    call ncdf_def_grp(ncdf,'SPARSE',grp)
-
-    call ncdf_def_dim(grp,'nnzs',n_nzs)
-
-    ! EDM is required to have its own spin (because of spin-orbit coupling
-    ! where the spin == 4, and not 8)
-    call ncdf_def_dim(grp,'spin_EDM',spin%EDM)
-
-
-    dic = dic//('info'.kv.'Index of supercell coordinates')
-    call ncdf_def_var(grp,'isc_off',NF90_INT,(/'xyz','n_s'/), &
-         compress_lvl=0, atts=dic)
-
-    dic = dic//('info'.kv.'Number of non-zero elements per row')
-    call ncdf_def_var(grp,'n_col',NF90_INT,(/'no_u'/), &
-         compress_lvl=0, atts=dic)
-
-    chks = (/n_nzs,1,1/)
-    
-    dic = dic//('info'.kv. &
-         'Supercell column indices in the sparse format ')
-    call ncdf_def_var(grp,'list_col',NF90_INT,(/'nnzs'/), &
-         compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
-
-    ! Create the overlap matrix (we know it will not change)
-    dic = ('info'.kv.'Overlap matrix')
-    call ncdf_def_var(grp,'S',NF90_DOUBLE,(/'nnzs'/), &
-         compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
-    
-    dic = dic//('info'.kv.'Density matrix')
-    call ncdf_def_var(grp,'DM',NF90_DOUBLE,(/'nnzs','spin'/), &
-         compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
-       
-    dic = dic//('info'.kv.'Energy density matrix')
-    dic = dic//('unit'.kv.'Ry')
-    call ncdf_def_var(grp,'EDM',NF90_DOUBLE,(/'nnzs    ','spin_EDM'/), &
-         compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
-
-    if ( savehs .or. TS_HS_save ) then
-       ! Unit is already present in dictionary
-       dic = dic//('info'.kv.'Hamiltonian')
-       call ncdf_def_var(grp,'H',NF90_DOUBLE,(/'nnzs','spin'/), &
+    if ( .not. lis_MD ) then
+      
+      ! Create matrix group
+      call ncdf_def_grp(ncdf,'SPARSE',grp)
+      
+      call ncdf_def_dim(grp,'nnzs',n_nzs)
+      
+      ! EDM is required to have its own spin (because of spin-orbit coupling
+      ! where the spin == 4, and not 8)
+      call ncdf_def_dim(grp,'spin_EDM',spin%EDM)
+      
+
+      dic = dic//('info'.kv.'Index of supercell coordinates')
+      call ncdf_def_var(grp,'isc_off',NF90_INT,(/'xyz','n_s'/), &
+          compress_lvl=0, atts=dic)
+      
+      dic = dic//('info'.kv.'Number of non-zero elements per row')
+      call ncdf_def_var(grp,'n_col',NF90_INT,(/'no_u'/), &
+          compress_lvl=0, atts=dic)
+      
+      chks = (/n_nzs,1,1/)
+      
+      dic = dic//('info'.kv. &
+          'Supercell column indices in the sparse format ')
+      call ncdf_def_var(grp,'list_col',NF90_INT,(/'nnzs'/), &
+          compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
+      
+      ! Create the overlap matrix (we know it will not change)
+      dic = ('info'.kv.'Overlap matrix')
+      call ncdf_def_var(grp,'S',NF90_DOUBLE,(/'nnzs'/), &
+          compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
+
+      dic = dic//('info'.kv.'Density matrix')
+      call ncdf_def_var(grp,'DM',NF90_DOUBLE,(/'nnzs','spin'/), &
+          compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
+
+      dic = dic//('info'.kv.'Energy density matrix')
+      dic = dic//('unit'.kv.'Ry')
+      call ncdf_def_var(grp,'EDM',NF90_DOUBLE,(/'nnzs    ','spin_EDM'/), &
+          compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
+
+      if ( savehs .or. TS_HS_save ) then
+        ! Unit is already present in dictionary
+        dic = dic//('info'.kv.'Hamiltonian')
+        call ncdf_def_var(grp,'H',NF90_DOUBLE,(/'nnzs','spin'/), &
             compress_lvl=cdf_comp_lvl,atts=dic,chunks=chks)
+      end if
+
+      ! Even though I think we could do without, I add the
+      ! xij array to the file
+      ! Note that xij can be re-created using
+      !    nsc, isc_off and xa
+      !    if ( .not. Gamma ) then
+      !       dic = dic//('info'.kv. &
+      !            'Distance between orbital i and j')
+      !       dic = dic//('unit'.kv.'Bohr')
+      !       call ncdf_def_var(grp,'xij',NF90_DOUBLE,(/'xyz ','nnzs'/), &
+      !            compress_lvl=cdf_comp_lvl,atts=dic)
+      !    end if
     end if
 
-    ! Even though I think we could do without, I add the
-    ! xij array to the file
-    ! Note that xij can be re-created using
-    !    nsc, isc_off and xa
-!    if ( .not. Gamma ) then
-!       dic = dic//('info'.kv. &
-!            'Distance between orbital i and j')
-!       dic = dic//('unit'.kv.'Bohr')
-!       call ncdf_def_var(grp,'xij',NF90_DOUBLE,(/'xyz ','nnzs'/), &
-!            compress_lvl=cdf_comp_lvl,atts=dic)
-!    end if
-
-
     ! Delete the dictionary
     call delete(dic)
 
@@ -336,21 +353,21 @@
     call ncdf_def_var(grp,'MeshCutoff',NF90_DOUBLE,(/'one'/), &
          compress_lvl=0, atts=dic)
 
-    ! Create matrix group
-    if ( is_MD ) then
+    ! Create MD group
+    if ( lis_MD ) then
 
        call ncdf_def_grp(ncdf,'MD',grp)
 
        ! Create arrays for containers
        call ncdf_def_dim(grp,'MD',NF90_UNLIMITED)
 
-       dic = ('info'.kv.'Temperature')//('unit'.kv.'K')
+       dic = ('info'.kv.'Ionic temperature')//('unit'.kv.'K')
        call ncdf_def_var(grp,'T',NF90_DOUBLE,(/'MD'/), &
             compress_lvl=0, atts=dic) ! do not compress unlimited D
 
-       dic = dic//('info'.kv.'Kohn-Sham energy')//('unit'.kv.'Ry')
-       call ncdf_def_var(grp,'E_KS',NF90_DOUBLE,(/'MD'/), &
-            compress_lvl=0, atts=dic) ! do not compress unlimited D
+       !dic = dic//('info'.kv.'Kohn-Sham energy')//('unit'.kv.'Ry')
+       !call ncdf_def_var(grp,'E_KS',NF90_DOUBLE,(/'MD'/), &
+       !     compress_lvl=0, atts=dic) ! do not compress unlimited D
 
        dic = dic//('info'.kv.'Total energy')
        call ncdf_def_var(grp,'E_tot',NF90_DOUBLE,(/'MD'/), &
@@ -360,6 +377,10 @@
        call ncdf_def_var(grp,'E_kin',NF90_DOUBLE,(/'MD'/), &
             compress_lvl=0, atts=dic) ! do not compress unlimited D
 
+       dic = dic//('info'.kv.'XC energy')
+       call ncdf_def_var(grp,'E_xc',NF90_DOUBLE,(/'MD'/), &
+            compress_lvl=0, atts=dic) ! do not compress unlimited D
+
        dic = dic//('info'.kv.'Fermi level')
        call ncdf_def_var(grp,'Ef',NF90_DOUBLE,(/'MD'/), &
             compress_lvl=0, atts=dic) ! do not compress unlimited D
@@ -367,6 +388,10 @@
        dic = dic//('info'.kv.'Atomic coordinates')//('unit'.kv.'Bohr')
        call ncdf_def_var(grp,'xa',NF90_DOUBLE,(/'xyz ','na_u','MD  '/), &
             compress_lvl=0, atts=dic) ! do not compress unlimited D
+       
+       dic = dic//('info'.kv.'Cell stress')//('unit'.kv.'Ry/Bohr**3')
+       call ncdf_def_var(grp,'stress',NF90_DOUBLE,(/'xyz','xyz','MD '/), &
+           compress_lvl=0, atts=dic) ! do not compress unlimited D
 
        dic = dic//('info'.kv.'Atomic forces')//('unit'.kv.'Ry/Bohr')
        call ncdf_def_var(grp,'fa',NF90_DOUBLE,(/'xyz ','na_u','MD  '/), &
@@ -378,6 +403,32 @@
 
     end if
 
+    ! Create FC group
+    if ( lis_FC ) then
+      
+      call ncdf_def_grp(ncdf,'FC',grp)
+
+      ! Create dimension arrays
+      ! Create arrays for containers
+      call ncdf_def_dim(grp,'ia',ia2-ia1+1)
+      call ncdf_def_dim(grp,'m_p',2)
+
+      dic = dic//('info'.kv.'Displacement length')//('unit'.kv.'Bohr')
+      call ncdf_def_var(grp,'disp',NF90_DOUBLE,(/'one'/), &
+          compress_lvl=0, atts=dic) ! do not compress unlimited D
+      call ncdf_put_var(grp,'disp',dx)
+
+      dic = dic//('info'.kv.'Undisplaced atomic forces ')//('unit'.kv.'Ry/Bohr')
+      call ncdf_def_var(grp,'fa0',NF90_DOUBLE,(/'xyz ','na_u'/), &
+          compress_lvl=0, atts=dic) ! do not compress unlimited D
+
+      dic = dic//('info'.kv.'Displaced atomic forces')//('unit'.kv.'Ry/Bohr')
+      call ncdf_def_var(grp,'fa',NF90_DOUBLE,(/'xyz ','na_u','m_p ','xyz ', 'ia  '/), &
+          compress_lvl=0, atts=dic) ! do not compress unlimited D
+      
+    end if
+
+    
     call delete(dic)
 
     if ( isolve == SOLVE_TRANSI ) then
@@ -447,12 +498,6 @@
 
     call delete(dic)
 
-    ! Save the total charge and lasto
-    call ncdf_put_var(ncdf,'Qtot',Qtot)
-    if ( fixspin ) then
-      call ncdf_put_var(ncdf,'Qspin',total_spin)
-    end if
-
     call ncdf_put_var(ncdf,'lasto',lasto(1:na_u))
 
     if ( isolve == SOLVE_TRANSI ) then
@@ -541,7 +586,8 @@
     use siesta_geom, only: na_u, ucell, xa, va
     use siesta_geom, only: nsc, isc_off
     use sparse_matrices, only: sparse_pattern, block_dist
-    use sparse_matrices, only: S_1D, DM_2D, EDM_2D, H_2D
+    use sparse_matrices, only: S_1D
+    use sparse_matrices, only: DM_2D, EDM_2D, H_2D
 !    use sparse_matrices, only: xij_2D
     use m_stress, only : stress
     use m_forces, only: fa
@@ -603,7 +649,7 @@
 !    if ( .not. Gamma .and. ('xij' .in. dic_save) ) then
        ! Write the xij array, it will not change during SCF
 !       call cdf_w_d2D(grp,'xij',xij_2D)
-!    end if
+    !    end if
     if ( 'H' .in. dic_save ) &
          call cdf_w_d2D(grp,'H',H_2D)
     if ( 'DM' .in. dic_save ) &
@@ -679,7 +725,9 @@
 
     ! Used for saving variables
     integer :: no, nk
-    integer, allocatable :: aux(:)
+    integer, allocatable :: iaux(:)
+    real(dp), allocatable :: raux(:)
+
 
     ! Unluckily the new basis saves only
     ! saved on the IO node. 
@@ -751,8 +799,8 @@
 
        ! Create all projector variables...
        if ( nk > 0 ) then
-          dic = ('pjnl_l'.kv.NF90_INT)//('pjnl_n'.kv.NF90_INT)
-          dic = dic//('pjnl_ekb'.kv.NF90_DOUBLE)
+          dic = ('pjnl_l'.kv.NF90_INT)//('pjnl_j'.kv.NF90_DOUBLE)
+          dic = dic//('pjnl_n'.kv.NF90_INT)//('pjnl_ekb'.kv.NF90_DOUBLE)
           dic = dic//('kbcutoff'.kv.NF90_DOUBLE)//('kbdelta'.kv.NF90_DOUBLE)
           d = .first. dic
           do while ( .not. (.empty. d) )
@@ -817,21 +865,27 @@
        call ncdf_put_var(grp,'orbnl_n',spp%orbnl_n(1:no))
        call ncdf_put_var(grp,'orbnl_z',spp%orbnl_z(1:no))
 
-       allocate(aux(no))
+       allocate(iaux(no))
        do i = 1, no
           if ( spp%orbnl_ispol(i) ) then
-             aux(i) = 1
+             iaux(i) = 1
           else
-             aux(i) = 0 
+             iaux(i) = 0 
           end if
        end do
-       call ncdf_put_var(grp,'orbnl_ispol',aux(1:no))
-       deallocate(aux)
+       call ncdf_put_var(grp,'orbnl_ispol',iaux(1:no))
+       deallocate(iaux)
        call ncdf_put_var(grp,'orbnl_pop',spp%orbnl_pop(1:no))
 
        ! Save projector
        if ( nk > 0 ) then
           call ncdf_put_var(grp,'pjnl_l',spp%pjnl_l(1:nk))
+          allocate(raux(nk))
+          do i = 1, nk
+            raux(i) = 0._dp
+          end do
+          call ncdf_put_var(grp,'pjnl_j',raux)
+          deallocate(raux)
           call ncdf_put_var(grp,'pjnl_n',spp%pjnl_n(1:nk))
           call ncdf_put_var(grp,'pjnl_ekb',spp%pjnl_ekb(1:nk))
        end if
@@ -874,9 +928,11 @@
 
   subroutine cdf_save_md(fname)
     use siesta_geom, only: na_u, xa, va
+    use m_stress, only: stress
     use m_forces, only: fa
-    use m_energies, only: Etot, Ef, Ekinion
+    use m_energies, only: Etot, Ekinion, Exc
     use m_kinetic, only : Tempion
+    use m_energies, only: Ef
 
     character(len=*), intent(in) :: fname
 
@@ -894,19 +950,54 @@
     if ( Node == 0 ) then
        ! Save current state
        call ncdf_put_var(ncdf,'T',Tempion,start=(/MD/))
-       call ncdf_put_var(ncdf,'Ef',Ef,start=(/MD/))
        !call ncdf_put_var(ncdf,'E_KS',EKS,start=(/MD/))
        call ncdf_put_var(ncdf,'E_tot',Etot,start=(/MD/))
        call ncdf_put_var(ncdf,'E_kin',Ekinion,start=(/MD/))
+       call ncdf_put_var(ncdf,'E_xc',Exc,start=(/MD/))
+       call ncdf_put_var(ncdf,'Ef',Ef,start=(/MD/))
        
        call ncdf_put_var(ncdf,'xa',xa,count=(/3,na_u/),start=(/1,1,MD/))
        call ncdf_put_var(ncdf,'fa',fa,count=(/3,na_u/),start=(/1,1,MD/))
+       call ncdf_put_var(ncdf,'stress',stress,count=(/3,3/),start=(/1,1,MD/))
        call ncdf_put_var(ncdf,'va',va,count=(/3,na_u/),start=(/1,1,MD/))
     end if
 
     call ncdf_close(ncdf)
 
   end subroutine cdf_save_md
+
+  subroutine cdf_save_fc(fname,istep)
+    use siesta_geom, only: na_u
+    use siesta_options, only: ia1
+    use m_forces, only: fa
+
+    character(len=*), intent(in) :: fname
+    integer, intent(in) :: istep
+
+    type(hNCDF) :: ncdf
+    integer :: ia, dir, pm
+
+    if ( Node /= 0 ) return
+
+    ! open the file...
+    call ncdf_open(ncdf,fname,mode=ior(NF90_WRITE,NF90_NETCDF4),group='FC')
+    
+    if ( istep == 0 ) then
+      call ncdf_put_var(ncdf,'fa0',fa)
+      
+    else
+      ! Atomic index (with respect to ia1)
+      ia = mod(istep - 1, 6) + 1
+      dir = (ia - 1) / 2 + 1
+      pm = mod(ia-1, 2) + 1
+      ia = (istep - 1) / 6 + 1
+      
+      call ncdf_put_var(ncdf,'fa',fa,count=(/3,na_u/),start=(/1,1,pm,dir,ia/))
+    end if
+    
+    call ncdf_close(ncdf)
+
+  end subroutine cdf_save_fc
   
 #else
   subroutine dummy()

=== modified file 'Src/read_options.F90'
--- Src/read_options.F90	2018-06-26 08:27:41 +0000
+++ Src/read_options.F90	2018-07-09 07:59:46 +0000
@@ -994,9 +994,11 @@
 
   ! Maximum number of steps in CG/Broyden coordinate optimization
   nmove = fdf_get('MD.NumCGsteps',0)
+  nmove = fdf_get('MD.Steps',nmove)
 
   ! Maximum atomic displacement in one CG step
   dxmax = fdf_get('MD.MaxCGDispl',0.2_dp,'Bohr')
+  dxmax = fdf_get('MD.MaxDispl',dxmax,'Bohr')
 
   ! Tolerance in the maximum atomic force [0.04 eV/Ang]
   ftol = fdf_get('MD.MaxForceTol', 0.00155574_dp, 'Ry/Bohr')
@@ -1052,11 +1054,19 @@
                    name  = 'MD.NumCGSteps', &
                    value = nmove,           &
                    units = "cmlUnits:countable" )
+              call cmlAddParameter( xf    = mainXML,         &
+                   name  = 'MD.Steps', &
+                   value = nmove,           &
+                   units = "cmlUnits:countable" )
 
               call cmlAddParameter( xf    = mainXML,           &
                    name  = 'MD.MaxCGDispl',   &
                    value = dxmax,             &
                    units = 'siestaUnits:Bohr' )
+              call cmlAddParameter( xf    = mainXML,           &
+                   name  = 'MD.MaxDispl',   &
+                   value = dxmax,             &
+                   units = 'siestaUnits:Bohr' )
 
               call cmlAddParameter( xf=mainXML,                 &
                    name='MD.MaxForceTol',      &
@@ -1156,7 +1166,11 @@
 
   ! Initial and final time steps for MD
   istart = fdf_get('MD.InitialTimeStep',1)
-  ifinal = fdf_get('MD.FinalTimeStep',1)
+  if ( fdf_defined('MD.Steps') ) then
+    ifinal = fdf_get('MD.FinalTimeStep',max(1,nmove - istart + 1))
+  else
+    ifinal = fdf_get('MD.FinalTimeStep',1)
+  end if
 
   ! Length of time step for MD
   dt = fdf_get('MD.LengthTimeStep',1._dp,'fs')
@@ -1583,7 +1597,7 @@
 
 
   ! Default variable depending on the action required
-  if ( idyn == 0 .and. nmove <= 1 ) then
+  if ( idyn == 0 .and. nmove <= 0 ) then
     ! Single-point calculation: use the tight auxiliary supercell by default
     naive_aux_cell = fdf_get( 'NaiveAuxiliaryCell', .false. )
   else

=== modified file 'Src/save_density_matrix.F'
--- Src/save_density_matrix.F	2017-11-23 14:17:27 +0000
+++ Src/save_density_matrix.F	2018-07-09 07:59:46 +0000
@@ -28,7 +28,7 @@
       use atomlist,             only:  no_l
       use siesta_options,       only:  writedm, writedm_cdf
       use siesta_options,       only:  writedm_cdf_history
-      use siesta_options,       only:  idyn
+      use siesta_options,       only:  idyn, nmove
       use m_steps,              only:  istp
       use sparse_matrices,      only:  Dscf, DM_2D
       use fdf,                  only:  fdf_get
@@ -51,7 +51,7 @@
       logical, intent(in), optional  :: when
       character(len=*), intent(in), optional   :: file
 
-      logical :: do_it
+      logical :: do_it, is_MD
 
 #ifdef CDF
 #ifdef NCDF_4
@@ -83,13 +83,17 @@
      &                         overwrite=.true. )
       endif
 #ifdef NCDF_4
-      if ( write_cdf ) then
-
-         dic_save = ('Ef'.kv.1)
-         dic_save = dic_save//('DM'.kv.1)//('EDM'.kv.1)
-         call cdf_save_state(trim(slabel)//'.nc',dic_save)
-
-         call delete(dic_save)
+      if ( write_cdf .and. do_it ) then
+
+        if ( (idyn == 0 .and. nmove == 0) .or.
+     &      (idyn == 6 .and. istp == 1 ) ) then
+          
+          dic_save = ('Ef'.kv.1)
+          dic_save = dic_save//('DM'.kv.1)//('EDM'.kv.1)
+          call cdf_save_state(trim(slabel)//'.nc',dic_save)
+          
+          call delete(dic_save)
+        end if
 
       end if
 #endif

=== modified file 'Src/siesta_forces.F'
--- Src/siesta_forces.F	2018-06-19 11:00:00 +0000
+++ Src/siesta_forces.F	2018-07-09 07:59:46 +0000
@@ -72,9 +72,7 @@
       use write_subs,            only: siesta_write_forces
       use write_subs,            only: siesta_write_stress_pressure
 #ifdef NCDF_4
-      use dictionary
-      use m_ncdf_siesta, only : cdf_init_file, cdf_save_settings
-      use m_ncdf_siesta, only : cdf_save_state, cdf_save_basis
+      use m_ncdf_siesta, only : cdf_save_settings
 #endif
       use m_compute_energies, only: compute_energies
       use m_mpi_utils, only: broadcast
@@ -121,9 +119,6 @@
       external :: die, message
 
       real(dp) :: Qcur
-#ifdef NCDF_4
-      type(dict) :: d_sav
-#endif
 #ifdef SIESTA__FLOOK
       ! len=24 from m_mixing.F90
       character(len=1), target :: next_mixer(24)
@@ -230,27 +225,12 @@
       ! Initialize the NC file
       if ( write_cdf ) then
 
-         ! Initialize the file...
-         call cdf_init_file(trim(slabel)//'.nc',is_MD = .false.)
-#ifdef MPI
-         call MPI_Barrier(MPI_Comm_World,MPIerror)
-#endif
-
-         ! Save the settings
-         call cdf_save_settings(trim(slabel)//'.nc')
-#ifdef MPI
-         call MPI_Barrier(MPI_Comm_World,MPIerror)
-#endif
-
-         d_sav = ('sp'.kv.1)//('S'.kv.1)
-         d_sav = d_sav//('nsc'.kv.1)//('xij'.kv.1)
-         d_sav = d_sav//('xa'.kv.1)//('cell'.kv.1)
-         d_sav = d_sav//('isc_off'.kv.1)
-         call cdf_save_state(trim(slabel)//'.nc',d_sav)
-         call delete(d_sav)
-
-         ! Save the basis set
-         call cdf_save_basis(trim(slabel)//'.nc')
+!       Save the settings (important to do here since mesh-cutoff may
+!       have changed).
+        call cdf_save_settings(trim(slabel)//'.nc')
+#ifdef MPI
+        call MPI_Barrier(MPI_Comm_World,MPIerror)
+#endif
 
       end if
 #endif

=== modified file 'Src/state_init.F'
--- Src/state_init.F	2018-06-26 08:27:41 +0000
+++ Src/state_init.F	2018-07-09 07:59:46 +0000
@@ -22,9 +22,9 @@
       use sparse_matrices, only: maxnh, numh, listh, listhptr
       use sparse_matrices, only: Dold, Dscf, DM_2D
       use sparse_matrices, only: Eold, Escf, EDM_2D
-      use sparse_matrices, only: Hold, H   , H_2D
+      use sparse_matrices, only: Hold, H, H_2D
       use sparse_matrices, only: xijo, xij_2D
-      use sparse_matrices, only: S   , S_1D
+      use sparse_matrices, only: S, S_1D
 
       use sparse_matrices, only: H_kin_1D, H_vkb_1D
       use sparse_matrices, only: H_ldau_2D, H_so_2D
@@ -91,6 +91,15 @@
       use iodm_netcdf, only: setup_dm_netcdf_file
       use iodmhs_netcdf, only: setup_dmhs_netcdf_file
 #endif
+#ifdef NCDF_4
+      use dictionary
+      use m_ncdf_siesta, only : cdf_init_file, cdf_save_settings
+      use m_ncdf_siesta, only : cdf_save_state, cdf_save_basis
+#ifdef MPI
+      use mpi_siesta
+#endif
+#endif
+
       use class_Sparsity
       use class_dSpData1D
       use class_dSpData2D
@@ -123,6 +132,13 @@
       character(len=label_length+6) :: fname
       real(dp)                      :: dummyef=0.0, dummyqtot=0.0
       type(Sparsity) :: g_Sp
+#ifdef NCDF_4
+      type(dict) :: d_sav
+#ifdef MPI
+      integer :: MPIerror
+#endif
+      logical :: is_MD, is_FC
+#endif
       
       character(len=256)            :: oname
 
@@ -555,8 +571,42 @@
       call overlap( na_u, na_s, no_s, scell, xa, indxua, rmaxo, maxnh,
      &              lasto, iphorb, isa, numh, listhptr, listh, S )
 
-
+#ifdef NCDF_4
+!     At this point the sparsity pattern, overlap matrix and some
+!     other details.
+!     Before continuing we will create the CDF file output
 !
+!     Initialize the NC file
+      if ( write_cdf ) then
+        if ( inicoor == istep ) then
+          is_MD = .true.
+          if ( idyn == 0 .and. nmove == 0 ) is_MD = .false.
+          is_FC = (idyn == 6)
+          if ( is_FC ) is_MD = .false.
+          call cdf_init_file(trim(slabel)//'.nc',
+     &        is_MD=is_MD, is_FC=is_FC)
+#ifdef MPI
+          call MPI_Barrier(MPI_Comm_World,MPIerror)
+#endif
+
+!         Save the basis set (only once)
+          call cdf_save_basis(trim(slabel)//'.nc')
+#ifdef MPI
+          call MPI_Barrier(MPI_Comm_World,MPIerror)
+#endif
+
+          if ( .not. is_MD ) then
+            d_sav = d_sav//('xa'.kv.1)//('cell'.kv.1)
+            d_sav = d_sav//('sp'.kv.1)//('S'.kv.1)//('xij'.kv.1)
+            d_sav = d_sav//('isc_off'.kv.1)//('nsc'.kv.1)
+            call cdf_save_state(trim(slabel)//'.nc',d_sav)
+            call delete(d_sav)
+          end if
+        end if
+      end if
+#endif
+
+!     
 !     Here we could also read a Hamiltonian, either to proceed to
 !     the analysis section (with nscf=0) or to start a mix-H scf cycle.
 !

=== modified file 'Util/COOP/Makefile'
--- Util/COOP/Makefile	2018-06-28 13:53:17 +0000
+++ Util/COOP/Makefile	2018-07-09 07:59:46 +0000
@@ -814,14 +814,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/Denchar/Src/Makefile'
--- Util/Denchar/Src/Makefile	2018-06-28 13:53:17 +0000
+++ Util/Denchar/Src/Makefile	2018-07-09 07:59:46 +0000
@@ -888,14 +888,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/Gen-basis/Makefile'
--- Util/Gen-basis/Makefile	2018-06-28 13:53:17 +0000
+++ Util/Gen-basis/Makefile	2018-07-09 07:59:46 +0000
@@ -878,14 +878,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/Grimme/Makefile'
--- Util/Grimme/Makefile	2018-06-28 13:53:17 +0000
+++ Util/Grimme/Makefile	2018-07-09 07:59:46 +0000
@@ -816,14 +816,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/Helpers/Makefile'
--- Util/Helpers/Makefile	2018-06-28 13:53:17 +0000
+++ Util/Helpers/Makefile	2018-07-09 07:59:46 +0000
@@ -819,14 +819,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/STM/ol-stm/Src/Makefile'
--- Util/STM/ol-stm/Src/Makefile	2018-06-28 13:53:17 +0000
+++ Util/STM/ol-stm/Src/Makefile	2018-07-09 07:59:46 +0000
@@ -901,14 +901,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/SpPivot/Makefile'
--- Util/SpPivot/Makefile	2018-06-28 13:53:17 +0000
+++ Util/SpPivot/Makefile	2018-07-09 07:59:46 +0000
@@ -830,14 +830,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/TS/TBtrans/Makefile'
--- Util/TS/TBtrans/Makefile	2018-06-28 13:53:17 +0000
+++ Util/TS/TBtrans/Makefile	2018-07-09 07:59:46 +0000
@@ -1086,14 +1086,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/TS/ts2ts/Makefile'
--- Util/TS/ts2ts/Makefile	2018-06-28 13:53:17 +0000
+++ Util/TS/ts2ts/Makefile	2018-07-09 07:59:46 +0000
@@ -823,14 +823,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/TS/tshs2tshs/Makefile'
--- Util/TS/tshs2tshs/Makefile	2018-06-28 13:53:17 +0000
+++ Util/TS/tshs2tshs/Makefile	2018-07-09 07:59:46 +0000
@@ -889,14 +889,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'Util/VCA/Makefile'
--- Util/VCA/Makefile	2018-06-28 13:53:17 +0000
+++ Util/VCA/Makefile	2018-07-09 07:59:46 +0000
@@ -849,14 +849,14 @@
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
 state_init.o: ioxv.o kpoint_grid.o ldau_specs.o local_sys.o m_energies.o m_eo.o
-state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o m_os.o
-state_init.o: m_pivot_methods.o m_rmaxh.o m_sparse.o m_sparsity_handling.o
-state_init.o: m_spin.o m_steps.o m_supercell.o m_test_io.o m_ts_charge.o
-state_init.o: m_ts_electype.o m_ts_global_vars.o m_ts_io.o m_ts_kpoints.o
-state_init.o: m_ts_options.o m_ts_sparse.o m_ts_tri_init.o normalize_dm.o
-state_init.o: overlap.o parallel.o proximity_check.o siesta_cml.o
-state_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: units.o write_subs.o zmatrix.o
+state_init.o: m_mixing.o m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
+state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
+state_init.o: m_ts_io.o m_ts_kpoints.o m_ts_options.o m_ts_sparse.o
+state_init.o: m_ts_tri_init.o normalize_dm.o overlap.o parallel.o
+state_init.o: proximity_check.o siesta_cml.o siesta_dicts.o siesta_geom.o
+state_init.o: siesta_options.o sparse_matrices.o units.o write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
 struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
 struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o

=== modified file 'version.info'
--- version.info	2018-07-06 20:00:09 +0000
+++ version.info	2018-07-09 07:59:46 +0000
@@ -1,1 +1,5 @@
+<<<<<<< TREE
 siesta-4.1--938
+=======
+siesta-4.1--937--cdf-md-3
+>>>>>>> MERGE-SOURCE