Merge lp:~nickpapior/siesta/trunk-k-file into lp:siesta

[Nick Papior]

Restructured k-point sampling

  Now k-points may be user-supplied for further finetuning.
  This may be extremely useful in various situations:
  1) Detailed analysis on PDOS for spin-orbit calculations
  2) Faster convergence with fewer k-points
  3) Improved convergence for bias calculations by adding zoom-in
     k-points in the bias-window.
  4) And others.

  Now the kpoints are associated with a type:
    kpoint_t
  this will allow one to read k-points using kgrid.File flag.
  The precedence is much like the older implementation:

  1) kgrid.MonkhorstPack
  2) kgrid.Cutoff
  3) kgrid.File

  in that order. If none are specified, then the Gamma-point will be used.

  This meant that several places the code needed updating.
  - Projected DOS k-points
  - TranSiesta k-points
  - SCF k-points

  Now all references to the kpoint_grid module is using kpoint_scf_m module
  which uses the type.
  This not only reduces lots of duplicated code, but it also allows
  one implementation to benefit across all code segments.

  Also added a test that reads in the k-points from a file (si_coop_kp_file).

This commit cleans up the k-point read utilities and enables user-defined k-grids.

I.e. SCF, PDOS and TS k-point reads are done using the same code.
Secondly, the k-points are now stored in a type which may be more easily passed around.

[lp:~nickpapior/siesta/trunk-k-file]

700. By Nick Papior on 2018-05-25

Made some memory usages smaller.

701. By Nick Papior on 2018-05-25

Updated Makefiles in Util directories

[Alberto Garcia (albertog)]

Very good overall. Just a couple of points:

1. The documentation for the Kgrid.file format in the manual is missing the first (index) column.

2. Apart from the sum_i{weight_i} = 1 check for the contents of the user-provided file, it might be good to add a sanity check on the sampling quality. For example, simulate the integration over the whole BZ of a simple 3D periodic model function epsilon(k) of the type that comes out of a one-parameter tight-binding. If the function is f, Sum_i{f(k_i)} should be close to 0 if the average band position is at zero.

I have a patch for 1. above and two other minor comments in the code:

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex 2018-05-24 19:58:35 +0000
+++ Docs/siesta.tex 2018-05-26 06:45:39 +0000
@@ -3474,14 +3474,14 @@
   An example input may be (not physically justified in any sense):
   \begin{shellexample}
     4
- 0.0 0.0 0.0 0.25
- 0.5 0.5 0,5 0.25
- 0.2 0.2 0,2 0.25
- 0.3 0.3 0,3 0.25
+ 1 0.0 0.0 0.0 0.25
+ 2 0.5 0.5 0,5 0.25
+ 3 0.2 0.2 0,2 0.25
+ 4 0.3 0.3 0,3 0.25
   \end{shellexample}
- The first 3 columns are the k-point specification for each of the
- reciprocal lattice vectors while the fourth number is the weight for
- the k-point.
+ The first column is an index; the next 3 columns are the k-point
+ specification for each of the reciprocal lattice vectors, and the
+ last number is the weight for the k-point.

   \siesta\ checks whether the sum of weights equals 1. If not,
   \siesta\ will die.

=== modified file 'Src/kpoint_pdos.F90'
--- Src/kpoint_pdos.F90 2018-05-24 19:58:35 +0000
+++ Src/kpoint_pdos.F90 2018-05-26 06:44:10 +0000
@@ -42,7 +42,7 @@
     if ( kpoints_pdos%method == K_METHOD_NONE ) then

       ! The user hasn't specified anything.
- ! This means that we will use the
+ ! This means that we will use the scf sampling grid
       call kpoint_delete(kpoints_pdos)
       call kpoint_read(kpoints_pdos, '', ucell, TrSym)

=== modified file 'Src/kpoint_t.F90'
--- Src/kpoint_t.F90 2018-05-24 19:58:35 +0000
+++ Src/kpoint_t.F90 2018-05-26 06:44:03 +0000
@@ -152,7 +152,7 @@
     call kpoint_fdf_name(prefix, 'kgrid.MonkhorstPack', name_fdf)
     if ( fdf_block(name_fdf, bfdf) ) then

- ! The method is a k-point grid
+ ! The method is a k-point MP grid
       this%method = K_METHOD_MONKHORST_PACK

       ! Now read data

review: Needs Fixing

[Nick Papior (nickpapior)]

Thanks for the comments!

I have just amended a fix for your comments.

As for your suggestion 2 it requires the epsilon(k) to be rotational invariant such that one can be sure that all irreducible wedges (for any crystal) will yield the correct result. Or, did you have something else in mind?

[lp:~nickpapior/siesta/trunk-k-file]

702. By Nick Papior on 2018-05-27

Fixed Alberto's suggestions (except point #2).

Amended documentation for the kgrid.File format.

[Alberto Garcia (albertog)]

OK. The sanity check with the model function could wait, and it should be optional in any case, since you might need a user-defined list of k-points for many things besides full sampling.

Approved.

review: Approve

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2018-05-15 13:01:22 +0000
+++ Docs/siesta.tex	2018-05-27 11:17:14 +0000
@@ -3397,37 +3397,9 @@
 
 \subsection{\texorpdfstring{$k$}{k}-point sampling}
 
-These are options for the k-point grid used in the SCF cycle. For
-other specialized grids, see the Macroscopic Polarization and Density
-of States sections.
-
-\begin{fdfentry}{kgrid!Cutoff}[length]<$0.\,\mathrm{Bohr}$>
-
-  Parameter which determines the fineness of the k-grid used for
-  Brillouin zone sampling.  It is half the length of the smallest
-  lattice vector of the supercell required to obtain the same sampling
-  precision with a single k point.  Ref: Moreno and Soler, PRB 45,
-  13891 (1992).
-
-  \textit{Use:} If it is zero, only the gamma point is used.  The resulting
-  k-grid is chosen in an optimal way, according to the method of Moreno
-  and Soler (using an effective supercell which is as spherical as
-  possible, thus minimizing the number of k-points for a given
-  precision). The grid is displaced for even numbers of effective mesh
-  divisions.  This parameter is not used if \fdf{kgrid!MonkhorstPack}
-  is specified. If the unit cell changes during the calculation (for
-  example, in a cell-optimization run, the k-point
-  grid will change accordingly (see \fdf{ChangeKgridInMD} for the case
-  of variable-cell molecular-dynamics runs, such as Parrinello-Rahman).
-  This is analogous to the changes in the
-  real-space grid, whose fineness is specified by an energy cutoff. If
-  sudden changes in the number of k-points are not desired, then the
-  Monkhorst-Pack data block should be used instead. In this case there
-  will be an implicit change in the quality of the sampling as the cell
-  changes. Both methods should be equivalent for a well-converged
-  sampling.
-  
-\end{fdfentry}
+These are options for the $k$-point grid used in the SCF cycle. For
+other specialized grids, see Secs.~\ref{sec:macroscopic-polarization}
+and \ref{sec:dos}.
 
 \begin{fdfentry}{kgrid!MonkhorstPack}[block]<$\Gamma$-point>
   
@@ -3464,6 +3436,59 @@
 
 \end{fdfentry}
 
+\begin{fdfentry}{kgrid!Cutoff}[length]<$0.\,\mathrm{Bohr}$>
+
+  Parameter which determines the fineness of the $k$-grid used for
+  Brillouin zone sampling. It is half the length of the smallest
+  lattice vector of the supercell required to obtain the same sampling
+  precision with a single k point.  Ref: Moreno and Soler, PRB 45,
+  13891 (1992).
+
+  \textit{Use:} If it is zero, only the gamma point is used.  The resulting
+  k-grid is chosen in an optimal way, according to the method of Moreno
+  and Soler (using an effective supercell which is as spherical as
+  possible, thus minimizing the number of k-points for a given
+  precision). The grid is displaced for even numbers of effective mesh
+  divisions.  This parameter is not used if \fdf{kgrid!MonkhorstPack}
+  is specified. If the unit cell changes during the calculation (for
+  example, in a cell-optimization run, the k-point
+  grid will change accordingly (see \fdf{ChangeKgridInMD} for the case
+  of variable-cell molecular-dynamics runs, such as Parrinello-Rahman).
+  This is analogous to the changes in the
+  real-space grid, whose fineness is specified by an energy cutoff. If
+  sudden changes in the number of k-points are not desired, then the
+  Monkhorst-Pack data block should be used instead. In this case there
+  will be an implicit change in the quality of the sampling as the cell
+  changes. Both methods should be equivalent for a well-converged
+  sampling.
+  
+\end{fdfentry}
+
+\begin{fdfentry}{kgrid!File}[string]<none>
+
+  Specify a file from where the $k$-points are read in.  The format of
+  the file is identical to the \sysfile{KP} file with the exception
+  that the $k$-points are given in units of the reciprocal lattice
+  vectors. I.e. the range of the $k$-points are $]-1/2 ; 1/2]$.
+
+  An example input may be (not physically justified in any sense):
+  \begin{shellexample}
+    4
+    1 0.0 0.0 0.0 0.25
+    2 0.5 0.5 0.5 0.25
+    3 0.2 0.2 0.2 0.25
+    4 0.3 0.3 0.3 0.25
+  \end{shellexample}
+  The first integer specifies the total number of $k$-points in the
+  file.  The first column is an index; the next 3 columns are the
+  $k$-point specification for each of the reciprocal lattice vectors
+  while the fifth column is the weight for the $k$-point.
+
+  \siesta\ checks whether the sum of weights equals 1. If not,
+  \siesta\ will die.
+
+\end{fdfentry}
+
 \begin{fdflogicalF}{ChangeKgridInMD}
 
   If \fdftrue, the $k$-point grid is recomputed at every
@@ -3870,8 +3895,8 @@
 \fdf{SCF.H!Tolerance} and high values of \fdf{MeshCutoff}.  We
 encourage the user to test carefully these options for each system. An
 additional point to take into account is the mixing scheme
-employed. You are encouraged to use \fdf{SCF.Mix:hamiltonian}
-(currently is set up by default) instead of density matrix mixing,
+employed. You are encouraged to use \fdf{SCF.Mix:Hamiltonian}
+(currently this is the default) instead of density matrix mixing,
 since it speeds up the convergence.  The pseudopotentials have to be
 properly generated and tested for each specific system and they have
 to be in their fully relativistic form, together with the non-linear
@@ -4040,14 +4065,17 @@
 
   \begin{fdfoptions}
     \option[Hamiltonian]%
+    \fdfindex*{SCF.Mix:Hamiltonian}%
     \index{SCF!mixing!Hamiltonian}
     Mix the Hamiltonian matrix (default).
 
     \option[density]%
+    \fdfindex*{SCF.Mix:density}%
     \index{SCF!mixing!Density}
     Mix the density matrix.
 
     \option[charge]%
+    \fdfindex*{SCF.Mix:charge}%
     \index{SCF!mixing!Charge}
     Mix the real-space charge density. Note this is an experimental
     feature.
@@ -7558,7 +7586,8 @@
 
 
 
-\subsection{Densities of states}
+\subsection{Density of states}
+\label{sec:dos}
 
 \subsubsection{Total density of states}
 There are several options to obtain the
@@ -8001,7 +8030,7 @@
 
 
 \subsection{Macroscopic polarization}
-
+\label{sec:macroscopic-polarization}
 
 \begin{fdfentry}{PolarizationGrids}[block]
   \index{bulk polarization}%
@@ -10533,18 +10562,20 @@
 
 \begin{fdfentry}{MD.TypeOfRun}[string]<TDED>
 
-A new Molecular Dynamics option, \fdf{TDED}, has been
-included. The second run of \siesta\ uses this option with the files
-\sysfile{TDWF} and \sysfile{TDXV} present in the working directory.  In this option ions
-and electrons are assumed to move simultaneously.  The occupied
-electronic states are time-evolved instead of the usual SCF
-calculations in each step.  Choose this option even if you intend to
-do only-electron dynamics. If you want to do an electron dynamics-only
-calculation set \fdf{MD.FinalTimeStep} equal to $1$. For optical
-response calculations switch off the external field during the second
-run. The \fdf{MD.LengthTimeStep}, unlike in the standard MD simulation,
-is defined by mulitpilication of \fdf{TDED.TimeStep} and \fdf{TDED.Nsteps}.
-In TDDFT calculations, the user defined \fdf{MD.LengthTimeStep} is ignored.
+  A new Molecular Dynamics option, \fdf{TDED}, has been included. The
+  second run of \siesta\ uses this option with the files
+  \sysfile{TDWF} and \sysfile{TDXV} present in the working directory.
+  In this option ions and electrons are assumed to move
+  simultaneously. The occupied electronic states are time-evolved
+  instead of the usual SCF calculations in each step.  Choose this
+  option even if you intend to do only-electron dynamics. If you want
+  to do an electron dynamics-only calculation set
+  \fdf{MD.FinalTimeStep} equal to $1$. For optical response
+  calculations switch off the external field during the second
+  run. The \fdf{MD.LengthTimeStep}, unlike in the standard MD
+  simulation, is defined by mulitpilication of \fdf{TDED.TimeStep} and
+  \fdf{TDED.Nsteps}. In TDDFT calculations, the user defined
+  \fdf{MD.LengthTimeStep} is ignored.
 
 \end{fdfentry}
 
@@ -11535,6 +11566,34 @@
   
 \end{fdflogicalF}
 
+\subsection{\texorpdfstring{$k$}{k}-point sampling}
+
+The options for $k$-point sampling are identical to the \siesta\
+options, \fdf{kgrid!MonkhorstPack}, \fdf{kgrid!Cutoff} or
+\fdf{kgrid!File}.
+
+One may however use specific \tsiesta\ $k$-points by using these
+options:
+
+\begin{fdfentry}{TS.kgrid!MonkhorstPack}[block]<\fdfvalue{kgrid!MonkhorstPack}>%
+
+  See \fdf{kgrid!MonkhorstPack} for details.
+  
+\end{fdfentry}
+
+\begin{fdfentry}{TS.kgrid!Cutoff}[length]<$0.\,\mathrm{Bohr}$>
+
+  See \fdf{kgrid!Cutoff} for details.
+  
+\end{fdfentry}
+
+\begin{fdfentry}{TS.kgrid!File}[string]<none>
+
+  See \fdf{kgrid!File} for details.
+  
+\end{fdfentry}
+
+
 \subsubsection{Algorithm specific options}
 
 These options adhere to the specific solution methods available for
@@ -12576,7 +12635,7 @@
   \sysfile{TSDE}. This file also contains geometry information
   etc. needed by \tsiesta\ and \tbtrans.
 
-  \item[\sysfile{TSKP}]: The $k$-points used in the \tsiesta\ calculation. See
+  \item[\sysfile{TS.KP}]: The $k$-points used in the \tsiesta\ calculation. See
   \siesta\ \sysfile{KP} file for formatting information.
 
   \item[\sysfile{TSCCEQ*}]: The equilibrium complex contour integration paths.
 
=== modified file 'RELEASE_NOTES'
--- RELEASE_NOTES	2018-04-15 15:29:11 +0000
+++ RELEASE_NOTES	2018-05-27 11:17:14 +0000
@@ -18,6 +18,8 @@
 
 *** New features
 
+ - Enabled custom k-points in PDOS and SCF calculations.
+
  - Added the CheSS solver as a new solution method. This is
    a linear scaling solver which calculates the density matrix and the
    energy density matrix.
 
=== modified file 'Src/Makefile'
--- Src/Makefile	2018-05-15 13:01:22 +0000
+++ Src/Makefile	2018-05-27 11:17:14 +0000
@@ -103,9 +103,10 @@
 	savepsi.o shaper.o timer_tree.o timer.o \
 	vmb.o vmat.o vmatsp.o volcel.o \
         cgvc.o cgvc_zmatrix.o m_convergence.o \
-        iocg.o ioeig.o iofa.o iokp.o iomd.o kpoint_pdos.o typecell.o \
+	kpoint_t.o kpoint_scf.o kpoint_pdos.o \
+        iocg.o ioeig.o iofa.o iokp.o iomd.o typecell.o \
         ofc.o poison.o readsp.o radfft.o \
-	write_md_record.o kpoint_grid.o find_kgrid.o proximity_check.o\
+	write_md_record.o find_kgrid.o proximity_check.o\
 	state_init.o siesta_move.o setup_hamiltonian.o compute_dm.o mixer.o\
 	scfconvergence_test.o post_scf_work.o state_analysis.o write_subs.o \
 	siesta_init.o struct_init.o siesta_options.o read_options.o siesta_geom.o \
@@ -443,7 +444,7 @@
 OBJS += m_ts_io.o
 
 # TranSIESTA...
-TS_OBJS= m_ts_global_vars.o m_ts_options.o m_ts_aux.o m_ts_kpoints.o \
+TS_OBJS= m_ts_global_vars.o m_ts_options.o m_ts_aux.o ts_kpoint_scf.o \
          m_ts_iodm.o \
          m_ts_electrode.o m_ts_gf.o \
 	 m_ts_cctype.o ts_init.o ts_show_regions.o m_ts_electype.o \
@@ -629,7 +630,7 @@
 class_TriMat.o: alloc.o intrinsic_missing.o
 coceri.o: files.o periodic_table.o precision.o units.o
 compute_dm.o: atomlist.o class_SpData1D.o compute_ebs_shift.o diagon.o
-compute_dm.o: iodmhs_netcdf.o kpoint_grid.o m_chess.o m_dminim.o m_energies.o
+compute_dm.o: iodmhs_netcdf.o kpoint_scf.o m_chess.o m_dminim.o m_energies.o
 compute_dm.o: m_eo.o m_gamma.o m_hsx.o m_pexsi_driver.o m_rmaxh.o m_spin.o
 compute_dm.o: m_steps.o m_transiesta.o m_ts_global_vars.o m_zminim.o
 compute_dm.o: normalize_dm.o ordern.o parallel.o precision.o siesta_geom.o
@@ -656,7 +657,7 @@
 cranknic_evolg.o: m_matswinvers.o m_spin.o m_steps.o parallel.o parallelsubs.o
 cranknic_evolg.o: precision.o siesta_options.o sparse_matrices.o sys.o units.o
 cranknic_evolg.o: wavefunctions.o
-cranknic_evolk.o: atomlist.o cranknic_evolg.o kpoint_grid.o m_energies.o m_eo.o
+cranknic_evolk.o: atomlist.o cranknic_evolg.o kpoint_scf.o m_energies.o m_eo.o
 cranknic_evolk.o: m_spin.o parallel.o parallelsubs.o precision.o
 cranknic_evolk.o: siesta_options.o sparse_matrices.o units.o wavefunctions.o
 create_Sparsity_SC.o: class_Sparsity.o geom_helper.o intrinsic_missing.o
@@ -731,10 +732,10 @@
 final_H_f_stress.o: ldau.o ldau_specs.o m_dipol.o m_energies.o m_forces.o
 final_H_f_stress.o: m_gamma.o m_hsx.o m_mpi_utils.o m_ncdf_siesta.o m_ntm.o
 final_H_f_stress.o: m_spin.o m_steps.o m_stress.o m_ts_global_vars.o m_ts_io.o
-final_H_f_stress.o: m_ts_kpoints.o m_ts_options.o metaforce.o
-final_H_f_stress.o: molecularmechanics.o naefs.o nlefsm.o overfsm.o parallel.o
-final_H_f_stress.o: siesta_geom.o siesta_options.o sparse_matrices.o
-final_H_f_stress.o: spinorbit.o sys.o units.o
+final_H_f_stress.o: m_ts_options.o metaforce.o molecularmechanics.o naefs.o
+final_H_f_stress.o: nlefsm.o overfsm.o parallel.o siesta_geom.o
+final_H_f_stress.o: siesta_options.o sparse_matrices.o spinorbit.o sys.o
+final_H_f_stress.o: ts_kpoint_scf.o units.o
 find_kgrid.o: alloc.o minvec.o parallel.o precision.o units.o
 fire_optim.o: alloc.o m_fire.o m_mpi_utils.o parallel.o precision.o
 fire_optim.o: siesta_options.o units.o
@@ -785,10 +786,10 @@
 kinefsm.o: alloc.o atmfuncs.o mneighb.o new_matel.o parallel.o parallelsubs.o
 kinefsm.o: precision.o
 kpoint_convert.o: precision.o sys.o units.o
-kpoint_grid.o: find_kgrid.o m_spin.o minvec.o parallel.o precision.o
-kpoint_grid.o: siesta_cml.o siesta_options.o sys.o units.o
-kpoint_pdos.o: find_kgrid.o m_spin.o minvec.o parallel.o precision.o
-kpoint_pdos.o: siesta_options.o sys.o units.o
+kpoint_pdos.o: kpoint_t.o m_spin.o parallel.o precision.o siesta_options.o
+kpoint_scf.o: kpoint_t.o m_spin.o parallel.o precision.o siesta_options.o
+kpoint_t.o: alloc.o files.o find_kgrid.o m_char.o m_io.o m_os.o minvec.o
+kpoint_t.o: parallel.o precision.o siesta_cml.o units.o
 ksv.o: alloc.o atmfuncs.o densematrix.o ksvinit.o parallel.o precision.o sys.o
 ksvinit.o: alloc.o parallel.o precision.o sys.o
 ldau.o: alloc.o atm_types.o atmfuncs.o compute_max_diff.o ldau_specs.o
@@ -798,7 +799,7 @@
 ldau_specs.o: basis_specs.o basis_types.o interpolation.o m_cite.o parallel.o
 ldau_specs.o: precision.o pseudopotential.o radial.o sys.o units.o
 listsc.o: alloc.o
-local_DOS.o: atomlist.o dhscf.o diagon.o files.o kpoint_grid.o m_energies.o
+local_DOS.o: atomlist.o dhscf.o diagon.o files.o kpoint_scf.o m_energies.o
 local_DOS.o: m_eo.o m_forces.o m_gamma.o m_ntm.o m_spin.o parallel.o
 local_DOS.o: siesta_geom.o siesta_options.o sparse_matrices.o sys.o
 m_broyddj.o: alloc.o m_mpi_utils.o parallel.o precision.o sys.o
@@ -847,7 +848,7 @@
 m_hsx.o: atm_types.o atmfuncs.o atomlist.o files.o parallel.o parallelsubs.o
 m_hsx.o: precision.o siesta_geom.o sys.o
 m_initwf.o: alloc.o atomlist.o densematrix.o diag.o diag_option.o fermid.o
-m_initwf.o: kpoint_grid.o m_eo.o m_gamma.o m_memory.o m_spin.o parallel.o
+m_initwf.o: kpoint_scf.o m_eo.o m_gamma.o m_memory.o m_spin.o parallel.o
 m_initwf.o: parallelsubs.o precision.o sparse_matrices.o sys.o wavefunctions.o
 m_integrate.o: precision.o
 m_inversemm.o: precision.o
@@ -879,10 +880,10 @@
 m_ncdf_io.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 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_gamma.o m_kinetic.o
+m_ncdf_siesta.o: kpoint_scf.o m_energies.o m_forces.o m_gamma.o m_kinetic.o
 m_ncdf_siesta.o: m_ncdf_io.o m_ntm.o m_spin.o m_stress.o m_ts_electype.o
-m_ncdf_siesta.o: m_ts_kpoints.o m_ts_options.o parallel.o precision.o radial.o
-m_ncdf_siesta.o: siesta_geom.o siesta_options.o sparse_matrices.o timestamp.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 ts_kpoint_scf.o
 m_new_dm.o: alloc.o atomlist.o class_Data2D.o class_Fstack_Data1D.o
 m_new_dm.o: class_Fstack_Pair_Geometry_SpData2D.o class_Geometry.o
 m_new_dm.o: class_OrbitalDistribution.o class_Pair_Geometry_SpData2D.o
@@ -946,9 +947,9 @@
 m_transiesta.o: class_SpData2D.o class_Sparsity.o files.o m_interpolate.o
 m_transiesta.o: m_ts_charge.o m_ts_contour_eq.o m_ts_contour_neq.o
 m_transiesta.o: m_ts_electype.o m_ts_fullg.o m_ts_fullk.o m_ts_gf.o
-m_transiesta.o: m_ts_kpoints.o m_ts_method.o m_ts_mumpsg.o m_ts_mumpsk.o
-m_transiesta.o: m_ts_options.o m_ts_sparse.o m_ts_tri_common.o m_ts_tri_init.o
-m_transiesta.o: m_ts_trig.o m_ts_trik.o parallel.o precision.o units.o
+m_transiesta.o: m_ts_method.o m_ts_mumpsg.o m_ts_mumpsk.o m_ts_options.o
+m_transiesta.o: m_ts_sparse.o m_ts_tri_common.o m_ts_tri_init.o m_ts_trig.o
+m_transiesta.o: m_ts_trik.o parallel.o precision.o ts_kpoint_scf.o units.o
 m_trialorbitalclass.o: precision.o units.o
 m_trimat_invert.o: class_TriMat.o intrinsic_missing.o m_pivot_array.o
 m_trimat_invert.o: precision.o
@@ -1001,9 +1002,9 @@
 m_ts_fullk.o: class_SpData2D.o class_SpData2D.o class_Sparsity.o m_iterator.o
 m_ts_fullk.o: m_ts_cctype.o m_ts_contour_eq.o m_ts_contour_neq.o m_ts_debug.o
 m_ts_fullk.o: m_ts_dm_update.o m_ts_elec_se.o m_ts_electype.o m_ts_full_scat.o
-m_ts_fullk.o: m_ts_gf.o m_ts_kpoints.o m_ts_method.o m_ts_options.o
-m_ts_fullk.o: m_ts_sparse.o m_ts_sparse_helper.o m_ts_weight.o parallel.o
-m_ts_fullk.o: precision.o units.o
+m_ts_fullk.o: m_ts_gf.o m_ts_method.o m_ts_options.o m_ts_sparse.o
+m_ts_fullk.o: m_ts_sparse_helper.o m_ts_weight.o parallel.o precision.o
+m_ts_fullk.o: ts_kpoint_scf.o units.o
 m_ts_gf.o: m_os.o m_ts_cctype.o m_ts_contour_eq.o m_ts_contour_neq.o
 m_ts_gf.o: m_ts_electrode.o m_ts_electype.o parallel.o precision.o sys.o
 m_ts_gf.o: units.o
@@ -1017,9 +1018,6 @@
 m_ts_io_ctype.o: m_io.o parallel.o precision.o units.o
 m_ts_iodm.o: class_OrbitalDistribution.o class_SpData2D.o class_Sparsity.o
 m_ts_iodm.o: m_io_s.o m_os.o parallel.o precision.o
-m_ts_kpoints.o: files.o find_kgrid.o kpoint_grid.o m_ts_global_vars.o
-m_ts_kpoints.o: m_ts_tdir.o minvec.o parallel.o precision.o siesta_cml.o
-m_ts_kpoints.o: siesta_options.o sys.o
 m_ts_method.o: alloc.o fdf_extra.o geom_helper.o m_region.o m_ts_electype.o
 m_ts_mumps_init.o: class_OrbitalDistribution.o class_Sparsity.o
 m_ts_mumps_init.o: create_Sparsity_Union.o m_ts_electype.o m_ts_method.o
@@ -1036,10 +1034,10 @@
 m_ts_mumpsk.o: class_SpData2D.o class_SpData2D.o class_Sparsity.o
 m_ts_mumpsk.o: intrinsic_missing.o m_iterator.o m_ts_cctype.o m_ts_contour_eq.o
 m_ts_mumpsk.o: m_ts_contour_neq.o m_ts_debug.o m_ts_dm_update.o m_ts_elec_se.o
-m_ts_mumpsk.o: m_ts_electype.o m_ts_gf.o m_ts_kpoints.o m_ts_method.o
-m_ts_mumpsk.o: m_ts_mumps_init.o m_ts_mumps_scat.o m_ts_options.o m_ts_sparse.o
+m_ts_mumpsk.o: m_ts_electype.o m_ts_gf.o m_ts_method.o m_ts_mumps_init.o
+m_ts_mumpsk.o: m_ts_mumps_scat.o m_ts_options.o m_ts_sparse.o
 m_ts_mumpsk.o: m_ts_sparse_helper.o m_ts_weight.o parallel.o precision.o
-m_ts_mumpsk.o: units.o
+m_ts_mumpsk.o: ts_kpoint_scf.o units.o
 m_ts_options.o: files.o intrinsic_missing.o m_mixing.o m_mixing_scf.o m_os.o
 m_ts_options.o: m_ts_charge.o m_ts_chem_pot.o m_ts_contour.o m_ts_contour_eq.o
 m_ts_options.o: m_ts_contour_neq.o m_ts_electype.o m_ts_global_vars.o
@@ -1084,10 +1082,10 @@
 m_ts_trik.o: m_iterator.o m_mat_invert.o m_region.o m_trimat_invert.o
 m_ts_trik.o: m_ts_cctype.o m_ts_charge.o m_ts_contour_eq.o m_ts_contour_neq.o
 m_ts_trik.o: m_ts_debug.o m_ts_dm_update.o m_ts_elec_se.o m_ts_electype.o
-m_ts_trik.o: m_ts_gf.o m_ts_kpoints.o m_ts_method.o m_ts_options.o
-m_ts_trik.o: m_ts_sparse.o m_ts_sparse_helper.o m_ts_tri_common.o
-m_ts_trik.o: m_ts_tri_init.o m_ts_tri_scat.o m_ts_trimat_invert.o m_ts_weight.o
-m_ts_trik.o: parallel.o precision.o units.o
+m_ts_trik.o: m_ts_gf.o m_ts_method.o m_ts_options.o m_ts_sparse.o
+m_ts_trik.o: m_ts_sparse_helper.o m_ts_tri_common.o m_ts_tri_init.o
+m_ts_trik.o: m_ts_tri_scat.o m_ts_trimat_invert.o m_ts_weight.o parallel.o
+m_ts_trik.o: precision.o ts_kpoint_scf.o units.o
 m_ts_trimat_invert.o: class_TriMat.o m_pivot_array.o m_region.o
 m_ts_trimat_invert.o: m_trimat_invert.o m_ts_electype.o m_ts_method.o
 m_ts_trimat_invert.o: precision.o
@@ -1175,16 +1173,16 @@
 poison.o: alloc.o cellsubs.o chkgmx.o fft.o parallel.o precision.o sys.o
 post_scf_work.o: atomlist.o class_Fstack_Pair_Geometry_SpData2D.o
 post_scf_work.o: class_Geometry.o class_Pair_Geometry_SpData2D.o compute_dm.o
-post_scf_work.o: diagon.o final_H_f_stress.o kpoint_grid.o m_dminim.o
+post_scf_work.o: diagon.o final_H_f_stress.o kpoint_scf.o m_dminim.o
 post_scf_work.o: m_energies.o m_eo.o m_gamma.o m_spin.o m_steps.o m_zminim.o
 post_scf_work.o: mneighb.o parallel.o siesta_geom.o siesta_options.o
 post_scf_work.o: sparse_matrices.o
 print_spin.o: atomlist.o m_mpi_utils.o m_spin.o parallel.o precision.o
 print_spin.o: siesta_cml.o sparse_matrices.o
 printmatrix.o: alloc.o
-projected_DOS.o: alloc.o atomlist.o kpoint_grid.o kpoint_pdos.o m_eo.o
-projected_DOS.o: m_gamma.o m_spin.o parallel.o precision.o siesta_geom.o
-projected_DOS.o: siesta_options.o sparse_matrices.o sys.o units.o
+projected_DOS.o: alloc.o atomlist.o kpoint_pdos.o kpoint_scf.o m_eo.o m_gamma.o
+projected_DOS.o: m_spin.o parallel.o precision.o siesta_geom.o siesta_options.o
+projected_DOS.o: sparse_matrices.o sys.o units.o
 propor.o: precision.o sys.o
 proximity_check.o: chemical.o m_ts_global_vars.o mneighb.o parallel.o
 proximity_check.o: precision.o siesta_geom.o siesta_options.o units.o
@@ -1213,7 +1211,7 @@
 rhoofd.o: meshdscf.o meshphi.o parallel.o parallelsubs.o precision.o sys.o
 rhoofdsp.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o mesh.o
 rhoofdsp.o: meshdscf.o meshphi.o precision.o sys.o
-sankey_change_basis.o: alloc.o atomlist.o diag_option.o kpoint_grid.o m_gamma.o
+sankey_change_basis.o: alloc.o atomlist.o diag_option.o kpoint_scf.o m_gamma.o
 sankey_change_basis.o: m_matdiag.o m_spin.o parallel.o parallelsubs.o
 sankey_change_basis.o: precision.o sparse_matrices.o sys.o wavefunctions.o
 save_density_matrix.o: atomlist.o files.o iodm_netcdf.o m_energies.o m_iodm.o
@@ -1249,7 +1247,7 @@
 siesta2wannier90.o: files.o m_digest_nnkp.o m_spin.o m_trialorbitalclass.o
 siesta2wannier90.o: parallel.o precision.o siesta_options.o sys.o
 siesta_analysis.o: alloc.o atomlist.o bands.o basis_enthalpy.o dhscf.o files.o
-siesta_analysis.o: flook_siesta.o kpoint_grid.o ksv.o ksvinit.o local_DOS.o
+siesta_analysis.o: flook_siesta.o kpoint_scf.o ksv.o ksvinit.o local_DOS.o
 siesta_analysis.o: m_dipol.o m_energies.o m_eo.o m_forces.o m_gamma.o
 siesta_analysis.o: m_iodm_old.o m_mpi_utils.o m_ntm.o m_partial_charges.o
 siesta_analysis.o: m_pexsi_dos.o m_pexsi_local_dos.o m_spin.o m_steps.o
@@ -1258,7 +1256,7 @@
 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: atomlist.o files.o kpoint_scf.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_end.o: alloc.o bands.o densematrix.o diag.o extrae_eventllist.o
@@ -1269,7 +1267,7 @@
 siesta_forces.o: atomlist.o class_Fstack_Data1D.o class_SpData2D.o compute_dm.o
 siesta_forces.o: compute_energies.o compute_max_diff.o densematrix.o
 siesta_forces.o: dm_charge.o files.o final_H_f_stress.o flook_siesta.o
-siesta_forces.o: kpoint_grid.o m_check_walltime.o m_convergence.o m_energies.o
+siesta_forces.o: kpoint_scf.o m_check_walltime.o m_convergence.o m_energies.o
 siesta_forces.o: m_forces.o m_initwf.o m_iodm_old.o m_mixing.o m_mixing_scf.o
 siesta_forces.o: m_mpi_utils.o m_ncdf_siesta.o m_pexsi.o m_pexsi_driver.o
 siesta_forces.o: m_rhog.o m_spin.o m_steps.o m_stress.o m_transiesta.o
@@ -1283,8 +1281,8 @@
 siesta_geom.o: precision.o
 siesta_init.o: alloc.o atomlist.o bands.o bsc_xcmod.o
 siesta_init.o: class_Fstack_Pair_Geometry_SpData2D.o densematrix.o
-siesta_init.o: diag_option.o files.o flook_siesta.o ioxv.o kpoint_grid.o
-siesta_init.o: kpoint_pdos.o ksvinit.o m_check_walltime.o m_cite.o m_energies.o
+siesta_init.o: diag_option.o files.o flook_siesta.o ioxv.o kpoint_pdos.o
+siesta_init.o: kpoint_scf.o ksvinit.o m_check_walltime.o m_cite.o m_energies.o
 siesta_init.o: m_eo.o m_fixed.o m_forces.o m_gamma.o m_iostruct.o m_mpi_utils.o
 siesta_init.o: m_new_dm.o m_rmaxh.o m_spin.o m_steps.o m_supercell.o m_timer.o
 siesta_init.o: m_wallclock.o metaforce.o molecularmechanics.o object_debug.o
@@ -1301,7 +1299,7 @@
 siesta_move.o: siesta_master.o siesta_options.o sys.o units.o write_subs.o
 siesta_move.o: zm_broyden_optim.o zm_fire_optim.o zmatrix.o
 siesta_tddft.o: alloc.o atomlist.o compute_energies.o final_H_f_stress.o
-siesta_tddft.o: kpoint_grid.o m_energies.o m_eo.o m_evolve.o m_gamma.o
+siesta_tddft.o: kpoint_scf.o m_energies.o m_eo.o m_evolve.o m_gamma.o
 siesta_tddft.o: m_initwf.o m_iotddft.o m_mpi_utils.o m_spin.o m_steps.o
 siesta_tddft.o: overfsm.o parallel.o precision.o sankey_change_basis.o
 siesta_tddft.o: setup_H0.o setup_hamiltonian.o siesta_cml.o siesta_options.o
@@ -1323,16 +1321,16 @@
 state_init.o: alloc.o atomlist.o class_Data2D.o class_SpData1D.o
 state_init.o: class_SpData2D.o class_SpData2D.o class_Sparsity.o
 state_init.o: create_Sparsity_SC.o domain_decom.o files.o hsparse.o
-state_init.o: iodm_netcdf.o iodmhs_netcdf.o iotdxv.o ioxv.o kpoint_grid.o
-state_init.o: ldau_specs.o m_chess.o m_energies.o m_eo.o m_gamma.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: iodm_netcdf.o iodmhs_netcdf.o iotdxv.o ioxv.o kpoint_scf.o
+state_init.o: kpoint_t.o ldau_specs.o m_chess.o m_energies.o m_eo.o m_gamma.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_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: siesta_options.o sparse_matrices.o sys.o ts_kpoint_scf.o units.o
+state_init.o: write_subs.o zmatrix.o
 struct_init.o: alloc.o atmfuncs.o atomlist.o files.o iotdxv.o ioxv.o
 struct_init.o: m_exp_coord.o m_iostruct.o m_mpi_utils.o m_steps.o parallel.o
 struct_init.o: periodic_table.o siesta_cml.o siesta_geom.o siesta_master.o
@@ -1342,9 +1340,12 @@
 timer.o: timer_tree.o
 timer_tree.o: m_walltime.o
 transition_rate.o: alloc.o fermid.o parallel.o parallelsubs.o precision.o sys.o
-ts_init.o: m_fixed.o m_os.o m_ts_cctype.o m_ts_charge.o m_ts_electrode.o
-ts_init.o: m_ts_electype.o m_ts_gf.o m_ts_global_vars.o m_ts_kpoints.o
+ts_init.o: kpoint_scf.o m_fixed.o m_os.o m_ts_cctype.o m_ts_charge.o
+ts_init.o: m_ts_electrode.o m_ts_electype.o m_ts_gf.o m_ts_global_vars.o
 ts_init.o: m_ts_method.o m_ts_options.o parallel.o siesta_options.o
+ts_init.o: ts_kpoint_scf.o
+ts_kpoint_scf.o: kpoint_t.o m_spin.o m_ts_global_vars.o m_ts_tdir.o parallel.o
+ts_kpoint_scf.o: precision.o siesta_options.o
 ts_show_regions.o: m_region.o m_ts_electype.o m_ts_method.o parallel.o
 ts_show_regions.o: precision.o units.o
 typecell.o: precision.o
@@ -1356,7 +1357,7 @@
 vmatsp.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o mesh.o meshdscf.o
 vmatsp.o: meshphi.o precision.o
 vmb.o: m_fixed.o parallel.o precision.o sys.o
-wavefunctions.o: atomlist.o kpoint_grid.o m_gamma.o m_matswinvers.o m_spin.o
+wavefunctions.o: atomlist.o kpoint_scf.o m_gamma.o m_matswinvers.o m_spin.o
 wavefunctions.o: parallel.o parallelsubs.o precision.o sparse_matrices.o
 write_inp_wannier.o: alloc.o atmfuncs.o atomlist.o m_ntm.o m_orderbands.o
 write_inp_wannier.o: mneighb.o parallel.o parallelsubs.o precision.o
@@ -1370,7 +1371,7 @@
 write_subs.o: m_steps.o m_stress.o m_ts_global_vars.o parallel.o precision.o
 write_subs.o: siesta_cml.o siesta_geom.o siesta_options.o units.o zmatrix.o
 writewave.o: alloc.o atmfuncs.o atomlist.o densematrix.o diag.o diag_option.o
-writewave.o: get_kpoints_scale.o kpoint_grid.o m_spin.o parallel.o
+writewave.o: get_kpoints_scale.o kpoint_scf.o m_spin.o parallel.o
 writewave.o: parallelsubs.o precision.o siesta_geom.o sys.o units.o
 xc.o: alloc.o bsc_xcmod.o precision.o sys.o
 xml.o: precision.o
@@ -1432,6 +1433,9 @@
 diagmemory.o: memoryinfo.o
 f90sockets.o: fsockets.o
 fsiesta.o: fsiesta_sockets.o
+kpoint_pdos_m.o: kpoint_pdos.o
+kpoint_scf_m.o: kpoint_scf.o
+kpoint_t_m.o: kpoint_t.o
 listsc_module.o: listsc.o
 m_bessph.o: bessph.o
 m_born_charge.o: born_charge.o
@@ -1535,6 +1539,7 @@
 t_spin.o: m_spin.o
 timer_options.o: timer.o
 trialorbitalclass.o: m_trialorbitalclass.o
+ts_kpoint_scf_m.o: ts_kpoint_scf.o
 version_info.o: version.o
 write_subs_energies.o: write_subs.o
 write_subs_positions.o: write_subs.o
 
=== modified file 'Src/compute_dm.F'
--- Src/compute_dm.F	2018-04-19 10:08:47 +0000
+++ Src/compute_dm.F	2018-05-27 11:17:14 +0000
@@ -23,7 +23,7 @@
       use atomlist, only: qa, lasto, iphorb, iaorb, no_u, no_s, indxuo, 
      &                    qtot, Qtots, no_l
       use sys, only: die, bye
-      use Kpoint_grid
+      use kpoint_scf_m, only: kpoints_scf
       use m_energies,   only: Ebs, Ecorrec, Entropy, DE_NEGF
       use m_energies,   only: Ef, Efs
       use m_rmaxh
@@ -177,7 +177,8 @@
      &              no_l, maxnh, maxnh, no_u,
      &              numh, listhptr, listh, numh, listhptr, listh, 
      &              H, S, qtot, fixspin, qtots, temp, e1, e2,
-     &              gamma, xijo, indxuo, nkpnt, kpoint, kweight,
+     &              gamma, xijo, indxuo,
+     &              kpoints_scf%N, kpoints_scf%k, kpoints_scf%w,
      &              eo, qo, Dscf, Escf, ef, efs, Entropy, no_u,
      &              occtol, iscf, neigwanted)
         Ecorrec = 0.0_dp
@@ -204,8 +205,9 @@
         else
           call zminim(.false., PreviousCallDiagon, iscf, istp, no_l,
      &                 spin%H, no_u, maxnh, numh, listhptr, listh, Dscf,
-     &                 eta, qtots, no_s, xijo, indxuo, nkpnt, kpoint,
-     &                 kweight, H, S, H_kin)
+     &                 eta, qtots, no_s, xijo, indxuo,
+     &                 kpoints_scf%N, kpoints_scf%k, kpoints_scf%w,
+     &                 H, S, H_kin)
         end if
         Ecorrec = 0.0_dp
         Entropy = 0.0_dp
@@ -229,7 +231,8 @@
      &              no_l, maxnh, maxnh, no_u,
      &              numh, listhptr, listh, numh, listhptr, listh,
      &              H, S, qtot, fixspin, qtots, temp, e1, e2,
-     &              gamma, xijo, indxuo, nkpnt, kpoint, kweight,
+     &              gamma, xijo, indxuo,
+     &              kpoints_scf%N, kpoints_scf%k, kpoints_scf%w,
      &              eo, qo, Dscf, Escf, ef, efs, Entropy, no_u,
      &              occtol, iscf, neigwanted)
 
 
=== modified file 'Src/cranknic_evolk.F90'
--- Src/cranknic_evolk.F90	2017-10-05 16:44:37 +0000
+++ Src/cranknic_evolk.F90	2018-05-27 11:17:14 +0000
@@ -32,7 +32,7 @@
   USE sparse_matrices,       ONLY: H, S, numh, listh, listhptr, xijo, Dscf
   USE atomlist,              ONLY: no_u, no_l, indxuo
   USE m_spin,                ONLY: nspin
-  USE Kpoint_grid,           ONLY: kpoint, nkpnt
+  USE kpoint_scf_m,          ONLY: kpoints_scf
   USE wavefunctions,         ONLY: compute_tddm, wavef_ms, complx_0, complx_1
   USE m_energies,            ONLY: etot 
   USE m_eo,                  ONLY: qo, eo
@@ -64,7 +64,7 @@
   call timer( 'cn_evolk', 1)
   ! 
   !
-  DO ik = 1,nkpnt
+  DO ik = 1,kpoints_scf%N
     DO ispin =1,nspin
       call m_allocate ( Hauxms, no_u, no_u, m_storage)
       call m_allocate ( Sauxms, no_u, no_u, m_storage)
@@ -76,8 +76,8 @@
           ind  = listhptr(i) + j
           juo  = listh(ind)
           jo   = indxuo (juo)
-          kxij = kpoint(1,ik)*xijo(1,ind) + kpoint(2,ik)*xijo(2,ind) +        &     
-                 kpoint(3,ik)*xijo(3,ind)
+          kxij = kpoints_scf%k(1,ik)*xijo(1,ind) + kpoints_scf%k(2,ik)*xijo(2,ind) + &
+                 kpoints_scf%k(3,ik)*xijo(3,ind)
          ckxij =  cos(kxij)
          skxij = -sin(kxij)
          cvar1 = cmplx(H(ind,ispin)*ckxij,H(ind,ispin)*skxij,dp)
@@ -132,7 +132,7 @@
 
  USE wavefunctions
  USE m_spin,                ONLY: nspin
- USE Kpoint_grid,           ONLY: nkpnt
+ USE kpoint_scf_m,          ONLY: kpoints_scf
  USE cranknic_evolg,        ONLY: Uphi
  USE atomlist,              ONLY: no_u
  USE MatrixSwitch
@@ -165,14 +165,14 @@
    ENDIF
  
    IF(extrapol_H_tdks) THEN
-     ALLOCATE(firstimeK(nkpnt, nspin))
-     ALLOCATE(Hsave(nkpnt, nspin))
-     DO i=1,nkpnt
+     ALLOCATE(firstimeK(kpoints_scf%N, nspin))
+     ALLOCATE(Hsave(kpoints_scf%N, nspin))
+     DO i=1,kpoints_scf%N
        DO j=1,nspin
          call m_allocate (Hsave(i,j),no_u, no_u, m_storage) 
        END DO
      END DO
-   firstimeK(1:nkpnt,1:nspin) = .true.
+   firstimeK(1:kpoints_scf%N,1:nspin) = .true.
    END IF
  
    firsttime = .false.
 
=== modified file 'Src/final_H_f_stress.F'
--- Src/final_H_f_stress.F	2018-04-23 11:09:01 +0000
+++ Src/final_H_f_stress.F	2018-05-27 11:17:14 +0000
@@ -71,7 +71,7 @@
       use files,        only : label_length
       use m_ts_global_vars, only : TSrun
       use m_ts_options, only : TS_HS_save
-      use m_ts_kpoints, only : ts_Gamma, ts_kscell, ts_kdispl
+      use ts_kpoint_scf_m, only: ts_kpoints_scf, ts_gamma_scf
       use m_ts_io,      only : ts_write_TSHS,fname_TSHS, FC_index
 
 #ifdef FINAL_CHECK_HS
@@ -340,9 +340,9 @@
          ! This is "pure" MD and we only write consecutive numbers
          ! Together with this you cannot also save FC
          fname = fname_TSHS(slabel, istep = istep )
-         call ts_write_tshs(fname, .false., Gamma, ts_Gamma,
+         call ts_write_tshs(fname, .false., Gamma, ts_Gamma_scf,
      &        ucell, nsc, isc_off, na_u, no_s, spin%H, 
-     &        ts_kscell, ts_kdispl, 
+     &        ts_kpoints_scf%k_cell, ts_kpoints_scf%k_displ,
      &        xa, lasto, 
      &        H_2D, S_1D, indxuo,
      &        Ef, Qtot, Temp, istep, 0)
@@ -359,17 +359,17 @@
             !  6 = +z
             call FC_index(istep,ia1,io_istep,io_ia1)
             fname = fname_TSHS(slabel,istep = io_istep, ia1 = io_ia1)
-            call ts_write_tshs(fname, .false., Gamma, ts_Gamma,
+            call ts_write_tshs(fname, .false., Gamma, ts_Gamma_scf,
      &           ucell, nsc, isc_off, na_u, no_s, spin%H, 
-     &           ts_kscell, ts_kdispl, 
+     &           ts_kpoints_scf%k_cell, ts_kpoints_scf%k_displ,
      &           xa, lasto, 
      &           H_2D, S_1D, indxuo,
      &           Ef, Qtot, Temp, io_istep, io_ia1)
          else
             fname = fname_TSHS(slabel)
-            call ts_write_tshs(fname, .false., Gamma, ts_Gamma,
+            call ts_write_tshs(fname, .false., Gamma, ts_Gamma_scf,
      &           ucell, nsc, isc_off, na_u, no_s, spin%H,
-     &           ts_kscell, ts_kdispl, 
+     &           ts_kpoints_scf%k_cell, ts_kpoints_scf%k_displ,
      &           xa, lasto,
      &           H_2D, S_1D, indxuo,
      &           Ef, Qtot, Temp, 0, 0)
 
=== removed file 'Src/kpoint_grid.F90'
--- Src/kpoint_grid.F90	2018-04-14 23:14:53 +0000
+++ Src/kpoint_grid.F90	1970-01-01 00:00:00 +0000
@@ -1,228 +0,0 @@
-! ---
-! Copyright (C) 1996-2016	The SIESTA group
-!  This file is distributed under the terms of the
-!  GNU General Public License: see COPYING in the top directory
-!  or http://www.gnu.org/copyleft/gpl.txt .
-! See Docs/Contributors.txt for a list of contributors.
-! ---
-MODULE Kpoint_grid
-!
-! Contains data structures and routines to deal with the kpoint-grid
-! for the self-consistent calculation
-! Other uses (bands, optical, polarization) have their own structures.
-!
-  USE precision, only : dp
-
-  implicit none
-
-  public :: setup_kpoint_grid, scf_kgrid_first_time, gamma_scf, &
-            nkpnt, kweight, kpoint, kscell, kdispl
-  public :: eff_kgrid_cutoff
-
-  private
-  
-  logical                  :: scf_kgrid_first_time = .true.
-  logical                  :: gamma_scf
-  integer                  :: nkpnt             ! Total number of k-points
-  real(dp)                 :: eff_kgrid_cutoff  ! Effective kgrid_cutoff
-
-  real(dp),        pointer :: kweight(:) => null()
-  real(dp),        pointer :: kpoint(:,:) => null()
-
-  integer,  dimension(3,3) :: kscell = 0
-  real(dp), dimension(3)   :: kdispl = 0.0_dp
-  logical                  :: user_requested_mp = .false.
-  logical                  :: user_requested_cutoff = .false.
-
-  logical, save            :: time_reversal_symmetry = .true.
-  logical                  :: firm_displ = .false.
-
-  CONTAINS
-
-  subroutine setup_Kpoint_grid( ucell )
-  USE parallel, only  : Node
-  USE fdf, only       : fdf_defined, fdf_get
-  USE m_find_kgrid, only : find_kgrid
-  use m_spin, only: TrSym
-
-    implicit none
-    real(dp) :: ucell(3,3)
-    logical  :: spiral
-
-    if (scf_kgrid_first_time) then
-       spiral = fdf_defined('SpinSpiral')
-          ! Allow the user to control the use of time-reversal-symmetry
-          ! By default, it is on, except for "spin-spiral" calculations
-          ! and/or non-collinear calculations
-       time_reversal_symmetry = fdf_get(             &
-            "TimeReversalSymmetryForKpoints",   &
-            (.not. spiral) .and. TrSym)
-       call setup_scf_kscell(ucell, firm_displ)
-
-       scf_kgrid_first_time = .false.
-
-    else
-       if ( user_requested_mp    ) then
-          ! no need to set up the kscell again
-       else
-          ! This was wrong in the old code
-          call setup_scf_kscell(ucell, firm_displ)
-       endif
-    endif
-    
-    call find_kgrid(ucell,kscell,kdispl,firm_displ,     &
-                    time_reversal_symmetry,             &
-                    nkpnt,kpoint,kweight, eff_kgrid_cutoff)
-
-    gamma_scf =  (nkpnt == 1 .and.  &
-                  dot_product(kpoint(:,1),kpoint(:,1)) < 1.0e-20_dp)
-
-    if (Node .eq. 0) call siesta_write_k_points()
-
-  end subroutine setup_Kpoint_grid
-
-!--------------------------------------------------------------------
-  subroutine setup_scf_kscell( cell, firm_displ )
-
-! ***************** INPUT **********************************************
-! real*8  cell(3,3)  : Unit cell vectors in real space cell(ixyz,ivec)
-! ***************** OUTPUT *********************************************
-! logical firm_displ   : User-specified displacements (firm)?
-
-!   The relevant fdf labels are kgrid_cutoff and kgrid_Monkhorst_Pack.
-!   If both are present, kgrid_Monkhorst_Pack has priority. If none is
-!   present, the cutoff default is zero, producing only the gamma point.
-!   Examples of fdf data specifications:
-!     kgrid_cutoff  50. Bohr
-!     %block kgrid_Monkhorst_Pack  # Defines kscell and kdispl
-!     4  0  0   0.50               # (kscell(i,1),i=1,3), kdispl(1)
-!     0  4  0   0.50               # (kscell(i,2),i=1,3), kdispl(2)
-!     0  0  4   0.50               # (kscell(i,3),i=1,3), kdispl(3)
-!     %endblock kgrid_Monkhorst_Pack
-! **********************************************************************
-
-!  Modules
-
-      use precision,  only : dp
-      use parallel,   only : Node
-      use m_minvec,   only : minvec
-      use sys,        only : die
-      use fdf
-
-      implicit          none
-
-! Passed variables
-      real(dp), intent(in)   :: cell(3,3)
-      logical, intent(out)   :: firm_displ
-
-! Internal variables
-      integer           i, j, factor(3,3), expansion_factor
-      real(dp)          scmin(3,3),  vmod, cutoff
-      real(dp)          ctransf(3,3)
-      logical           mp_input
-
-      real(dp), parameter :: defcut = 0.0_dp
-      integer, dimension(3,3), parameter :: unit_matrix =  &
-                         reshape ((/1,0,0,0,1,0,0,0,1/), (/3,3/))
-
-      type(block_fdf)            :: bfdf
-      type(parsed_line), pointer :: pline
-
-
-      mp_input = fdf_block('kgrid_Monkhorst_Pack',bfdf)
-      if ( mp_input ) then
-         user_requested_mp = .true.
-         do i= 1, 3
-            if (.not. fdf_bline(bfdf,pline))            &
-              call die('setup_scf_kscell: ERROR in ' // &
-                       'kgrid_Monkhorst_Pack block')
-            kscell(1,i) = fdf_bintegers(pline,1)
-            kscell(2,i) = fdf_bintegers(pline,2)
-            kscell(3,i) = fdf_bintegers(pline,3)
-            if ( fdf_bnvalues(pline) > 3 ) then
-              kdispl(i) = mod(fdf_bvalues(pline,4), 1._dp)
-            else
-              kdispl(i) = 0._dp
-            end if
-         enddo
-         firm_displ = .true.
-
-      else
-
-         cutoff = fdf_physical('kgrid_cutoff',defcut,'Bohr')
-         if (cutoff /= defcut) then
-         !!  write(6,"(a,f10.5)") "Kgrid cutoff input: ", cutoff
-            user_requested_cutoff = .true.
-         endif
-
-         kdispl(1:3) = 0.0_dp  ! Might be changed later
-         firm_displ = .false.  ! In future we might add new options
-                               ! for user-specified displacements
-         
-         ! Find equivalent rounded unit-cell
-         call minvec( cell, scmin, ctransf )
-
-         expansion_factor = 1
-         do j = 1,3
-            factor(j,1:3) = 0
-            vmod = sqrt(dot_product(scmin(1:3,j),scmin(1:3,j)))
-            factor(j,j) = int(2.0_dp*cutoff/vmod) + 1
-            expansion_factor = expansion_factor * factor(j,j)
-         enddo
-         ! Generate actual supercell skeleton
-         kscell = matmul(ctransf, factor)
-         ! Avoid confusing permutations
-         if (expansion_factor == 1) then
-            kscell = unit_matrix
-         endif
-      endif
-
-    end subroutine setup_scf_kscell
-
-    subroutine siesta_write_k_points()
-      USE siesta_options, only: writek
-      USE units, only: Ang
-      USE siesta_cml
-
-      implicit none
-
-      integer  :: ik, ix, i
-      external :: iokp
-
-      if ( writek ) then
-         write(6,'(/,a)') 'siesta: k-point coordinates (Bohr**-1) and weights:'
-         write(6,'(a,i4,3f12.6,3x,f12.6)')                          &
-              ('siesta: ', ik, (kpoint(ix,ik),ix=1,3), kweight(ik), &
-              ik=1,nkpnt)
-      endif
-      ! Always write KP file
-      call iokp( nkpnt, kpoint, kweight )
-      write(6,'(/a,i6)')  'siesta: k-grid: Number of k-points =', nkpnt
-      write(6,'(a,f10.3,a)')  'siesta: k-grid: Cutoff (effective) =',  &
-           eff_kgrid_cutoff/Ang, ' Ang'
-      write(6,'(a)') 'siesta: k-grid: Supercell and displacements'
-      write(6,'(a,3i4,3x,f8.3)') 'siesta: k-grid: ',        &
-           (kscell(i,1),i=1,3), kdispl(1)
-      write(6,'(a,3i4,3x,f8.3)') 'siesta: k-grid: ',        &
-           (kscell(i,2),i=1,3), kdispl(2)
-      write(6,'(a,3i4,3x,f8.3)') 'siesta: k-grid: ',        &
-           (kscell(i,3),i=1,3), kdispl(3)
-      if (cml_p) then
-         call cmlStartPropertyList(xf=mainXML, title="k-points", &
-                                      dictRef="siesta:kpoints")
-         call cmlAddProperty(xf=mainXML, value=nkpnt, dictref='siesta:nkpnt', &
-                             units="cmlUnits:countable")
-         do ik=1, nkpnt
-           call cmlAddKPoint(xf=mainXML, coords=kpoint(:,ik), weight=kweight(ik))
-         enddo
-         call cmlAddProperty(xf=mainXML, value=eff_kgrid_cutoff/Ang,     & 
-              dictref='siesta:kcutof', units='siestaUnits:angstrom')
-         call cmlEndPropertyList(mainXML)
-         call cmlAddProperty(xf=mainXML, value=kscell, dictref='siesta:kscell', &
-                             units="siestaUnits:Ang")
-         call cmlAddProperty(xf=mainXML, value=kdispl, dictref='siesta:kdispl', &
-                             units="siestaUnits:Ang")
-      endif
-    end subroutine siesta_write_k_points
-
-END MODULE Kpoint_grid
 
=== modified file 'Src/kpoint_pdos.F90'
--- Src/kpoint_pdos.F90	2018-04-14 23:14:53 +0000
+++ Src/kpoint_pdos.F90	2018-05-27 11:17:14 +0000
@@ -5,217 +5,59 @@
 !  or http://www.gnu.org/copyleft/gpl.txt .
 ! See Docs/Contributors.txt for a list of contributors.
 ! ---
-MODULE Kpoint_pdos
-!
-! Contains data structures and routines to deal with the kpoint-grid
-! for the projected density of states calculation. 
-! Modelled on the equivalent module for the SCF calculation.
-!
-  USE precision, only : dp
+module kpoint_pdos_m
+  !
+  ! Contains data structures and routines to deal with the kpoint-grid
+  ! for the self-consistent calculation
+  ! Other uses (bands, optical, polarization) have their own structures.
+  !
+  use precision, only : dp
+  
+  ! The k-point-type
+  use kpoint_t_m
 
   implicit none
 
+  public :: setup_kpoint_pdos
+  public :: kpoints_pdos
+  public :: gamma_pdos
+
   private
+
+  logical, save :: gamma_pdos
+  type(kpoint_t), save :: kpoints_pdos
+
+contains
+
+  subroutine setup_kpoint_pdos( ucell )
+    use parallel, only: Node
+    use siesta_options, only: writek
+    use m_spin, only: TrSym
+
+    real(dp), intent(in) :: ucell(3,3)
+
+    ! First try and read the k-points
+    call kpoint_read(kpoints_pdos, 'PDOS', ucell, TrSym)
+
+    if ( kpoints_pdos%method == K_METHOD_NONE ) then
+
+      ! The user hasn't specified a specific PDOS k-point sampling.
+      call kpoint_delete(kpoints_pdos)
+      call kpoint_read(kpoints_pdos, '', ucell, TrSym)
+
+    end if
+
+    gamma_pdos = (kpoints_pdos%N == 1 .and. &
+        dot_product(kpoints_pdos%k(:,1),kpoints_pdos%k(:,1)) < 1.0e-20_dp)
+
+    ! Quick-return if non-IO
+    if ( Node /= 0 ) return
+
+    ! Write to XML file
+    call kpoint_write_stdout(kpoints_pdos, writek, 'PDOS')
+    call kpoint_write_xml(kpoints_pdos, 'PDOS')
+    call kpoint_write_file(kpoints_pdos, 'PDOS.KP')
+
+  end subroutine setup_kpoint_pdos
   
-  logical, public, save   :: pdos_kgrid_first_time = .true.
-  logical, public, save   :: gamma_pdos
-  integer, public, save   :: maxk_pdos         ! 
-  integer, public, save   :: nkpnt_pdos        ! Total number of k-points
-  real(dp)                :: eff_kgrid_cutoff  ! Effective kgrid_cutoff
-
-  real(dp), pointer, public, save :: kweight_pdos(:) 
-  real(dp), pointer, public, save :: kpoints_pdos(:,:)
-
-  integer,  dimension(3,3), save  :: kscell = 0
-  real(dp), dimension(3), save    :: kdispl = 0.0_dp
-
-  logical, save     :: user_requested_mp = .false.
-  logical, save     :: user_requested_cutoff = .false.
-
-  logical, save     :: time_reversal_symmetry = .true.
-  logical, save     :: firm_displ = .false.
-
-  public :: setup_kpoint_pdos
-
-  CONTAINS
-
-  subroutine setup_Kpoint_pdos( ucell, different_pdos_grid )
-  USE parallel, only  : Node
-  USE fdf, only       : fdf_defined, fdf_get
-  USE m_find_kgrid, only : find_kgrid
-  use m_spin, only: TrSym
-
-  implicit none
-  real(dp), intent(in)  :: ucell(3,3)
-  logical,  intent(out) :: different_pdos_grid
-
-  logical :: spiral
-  if (pdos_kgrid_first_time) then
-    nullify(kweight_pdos,kpoints_pdos)
-    spiral = fdf_defined('SpinSpiral')
-      ! Allow the user to control the use of time-reversal-symmetry
-      ! By default, it is on, except for "spin-spiral" calculations
-      ! and/or non-collinear calculations
-    time_reversal_symmetry = fdf_get(             &
-           "TimeReversalSymmetryForKpoints",   &
-           (.not. spiral) .and. TrSym)
-
-    call setup_pdos_kscell(ucell, firm_displ)
-
-    pdos_kgrid_first_time = .false.
-
-  else
-    if ( user_requested_mp    ) then
-      ! no need to set up the kscell again
-    else
-      ! This was wrong in the old code
-      call setup_pdos_kscell(ucell, firm_displ)
-    endif
-  endif
-
-  if (user_requested_mp.or.user_requested_cutoff) then
-    different_pdos_grid = .true.
-  else
-    different_pdos_grid = .false.
-  endif
-
-! If the grid hasn't been explicit specified then just set dummy values 
-  if (different_pdos_grid) then
-    call find_kgrid(ucell,kscell,kdispl,firm_displ, &
-                    time_reversal_symmetry,         &
-                    nkpnt_pdos,kpoints_pdos,kweight_pdos,eff_kgrid_cutoff)
-
-    maxk_pdos = nkpnt_pdos
-    gamma_pdos = (nkpnt_pdos == 1 .and. dot_product(kpoints_pdos(:,1),kpoints_pdos(:,1)) < 1.0e-20_dp)
-
-    if (Node .eq. 0) call siesta_write_k_points_pdos()
-  else
-    nkpnt_pdos = 0
-    maxk_pdos = nkpnt_pdos
-    gamma_pdos = .true.
-  endif
-
-  end subroutine setup_Kpoint_pdos
-
-!--------------------------------------------------------------------
-  subroutine setup_pdos_kscell( cell, firm_displ )
-
-! ***************** INPUT **********************************************
-! real*8  cell(3,3)  : Unit cell vectors in real space cell(ixyz,ivec)
-! ***************** OUTPUT *********************************************
-! logical firm_displ   : User-specified displacements (firm)?
-
-!   The relevant fdf labels are PDOS.kgrid_cutoff and PDOS.kgrid_Monkhorst_Pack.
-!   If both are present, PDOS.kgrid_Monkhorst_Pack has priority. If neither is
-!   present, the cutoff default is zero, producing only the gamma point.
-!   Examples of fdf data specifications:
-!     PDOS.kgrid_cutoff  50. Bohr
-!     %block PDOS.kgrid_Monkhorst_Pack  # Defines kscell and kdispl
-!     4  0  0   0.50                    # (kscell(i,1),i=1,3), kdispl(1)
-!     0  4  0   0.50                    # (kscell(i,2),i=1,3), kdispl(2)
-!     0  0  4   0.50                    # (kscell(i,3),i=1,3), kdispl(3)
-!     %endblock PDOS.kgrid_Monkhorst_Pack
-! **********************************************************************
-
-!  Modules
-
-    use precision,  only : dp
-    use m_minvec,   only : minvec
-    use sys,        only : die
-    use fdf
-
-    implicit          none
-
-! Passed variables
-    real(dp), intent(in)   :: cell(3,3)
-    logical, intent(out)   :: firm_displ
-
-! Internal variables
-    integer           i, j,  factor(3,3), expansion_factor
-
-    real(dp)          scmin(3,3),  vmod, cutoff
-    real(dp)          ctransf(3,3)
-    logical           mp_input
-
-    real(dp), parameter :: defcut = 0.0_dp
-    integer, dimension(3,3), parameter :: unit_matrix =  &
-                         reshape ((/1,0,0,0,1,0,0,0,1/), (/3,3/))
-
-      type(block_fdf)            :: bfdf
-      type(parsed_line), pointer :: pline
-
-      mp_input = fdf_block('PDOS.kgrid_Monkhorst_Pack',bfdf)
-      if ( mp_input ) then
-         user_requested_mp = .true.
-         do i= 1, 3
-            if (.not. fdf_bline(bfdf,pline))            &
-              call die('setup_pdos_kscell: ERROR in ' // &
-                       'PDOS.kgrid_Monkhorst_Pack block')
-            kscell(1,i) = fdf_bintegers(pline,1)
-            kscell(2,i) = fdf_bintegers(pline,2)
-            kscell(3,i) = fdf_bintegers(pline,3)
-            if ( fdf_bnvalues(pline) > 3 ) then
-               kdispl(i) = mod(fdf_bvalues(pline,4), 1._dp)
-            else
-               kdispl(i) = 0._dp
-            end if
-         enddo
-         firm_displ = .true.
-
-      else
-
-         cutoff = fdf_physical('PDOS.kgrid_cutoff',defcut,'Bohr')
-         if (cutoff /= defcut) then
-         !!  write(6,"(a,f10.5)") "PDOS Kgrid cutoff input: ", cutoff
-            user_requested_cutoff = .true.
-         endif
-
-         kdispl(1:3) = 0.0_dp  ! Might be changed later
-         firm_displ = .false.  ! In future we might add new options
-                               ! for user-specified displacements
-         
-         ! Find equivalent rounded unit-cell
-         call minvec( cell, scmin, ctransf )
-
-         expansion_factor = 1
-         do j = 1,3
-            factor(j,1:3) = 0
-            vmod = sqrt(dot_product(scmin(1:3,j),scmin(1:3,j)))
-            factor(j,j) = int(2.0_dp*cutoff/vmod) + 1
-            expansion_factor = expansion_factor * factor(j,j)
-         enddo
-         ! Generate actual supercell skeleton
-         kscell = matmul(ctransf, factor)
-         ! Avoid confusing permutations
-         if (expansion_factor == 1) then
-            kscell = unit_matrix
-         endif
-      endif
-
-  end subroutine setup_pdos_kscell
-
-  subroutine siesta_write_k_points_pdos()
-
-    USE siesta_options, only: writek
-    USE units, only: Ang
-
-    implicit none
-
-    integer  :: ik, ix, i
-
-    if ( writek ) then
-      write(6,'(/,a)') 'siesta: k-point coordinates (Bohr**-1) and weights for PDOS:'
-      write(6,'(a,i4,3f12.6,3x,f12.6)')                          &
-              ('siesta: ', ik, (kpoints_pdos(ix,ik),ix=1,3), kweight_pdos(ik), ik=1,nkpnt_pdos)
-    endif
-
-    write(6,'(/a,i6)')  'siesta: PDOS.k-grid: Number of k-points =', nkpnt_pdos
-    write(6,'(a,f10.3,a)')  'siesta: PDOS.k-grid: Cutoff (effective) =', eff_kgrid_cutoff/Ang, ' Ang'
-    write(6,'(a)') 'siesta: PDOS.k-grid: Supercell and displacements'
-    write(6,'(a,3i4,3x,f8.3)') 'siesta: PDOS.k-grid: ', (kscell(i,1),i=1,3), kdispl(1)
-    write(6,'(a,3i4,3x,f8.3)') 'siesta: PDOS.k-grid: ', (kscell(i,2),i=1,3), kdispl(2)
-    write(6,'(a,3i4,3x,f8.3)') 'siesta: PDOS.k-grid: ', (kscell(i,3),i=1,3), kdispl(3)
-
-  end subroutine siesta_write_k_points_pdos
-
-END MODULE Kpoint_pdos
+end module kpoint_pdos_m
 
=== added file 'Src/kpoint_scf.F90'
--- Src/kpoint_scf.F90	1970-01-01 00:00:00 +0000
+++ Src/kpoint_scf.F90	2018-05-27 11:17:14 +0000
@@ -0,0 +1,53 @@
+! ---
+! Copyright (C) 1996-2016	The SIESTA group
+!  This file is distributed under the terms of the
+!  GNU General Public License: see COPYING in the top directory
+!  or http://www.gnu.org/copyleft/gpl.txt .
+! See Docs/Contributors.txt for a list of contributors.
+! ---
+module kpoint_scf_m
+  !
+  ! Contains data structures and routines to deal with the kpoint-grid
+  ! for the self-consistent calculation
+  ! Other uses (bands, optical, polarization) have their own structures.
+  !
+  use precision, only : dp
+  
+  ! The k-point-type
+  use kpoint_t_m
+
+  implicit none
+
+  public :: setup_kpoint_scf
+  public :: kpoints_scf
+  public :: gamma_scf
+  
+  private
+
+  logical, save :: gamma_scf
+  type(kpoint_t), save :: kpoints_scf
+
+contains
+
+  subroutine setup_kpoint_scf( ucell )
+    use parallel, only: Node
+    use siesta_options, only: writek
+    use m_spin, only: TrSym
+
+    real(dp), intent(in) :: ucell(3,3)
+
+    call kpoint_read(kpoints_scf, '', ucell, TrSym)
+
+    gamma_scf = (kpoints_scf%N == 1 .and. &
+        dot_product(kpoints_scf%k(:,1),kpoints_scf%k(:,1)) < 1.0e-20_dp)
+
+    ! Quick-return if non-IO
+    if ( Node /= 0 ) return
+
+    call kpoint_write_stdout(kpoints_scf, all=writek)
+    call kpoint_write_xml(kpoints_scf)
+    call kpoint_write_file(kpoints_scf, 'KP')
+
+  end subroutine setup_kpoint_scf
+  
+end module kpoint_scf_m
 
=== added file 'Src/kpoint_t.F90'
--- Src/kpoint_t.F90	1970-01-01 00:00:00 +0000
+++ Src/kpoint_t.F90	2018-05-27 11:17:14 +0000
@@ -0,0 +1,598 @@
+! ---
+! Copyright (C) 1996-2016	The SIESTA group
+!  This file is distributed under the terms of the
+!  GNU General Public License: see COPYING in the top directory
+!  or http://www.gnu.org/copyleft/gpl.txt .
+! See Docs/Contributors.txt for a list of contributors.
+! ---
+module kpoint_t_m
+
+  use precision, only: dp
+
+  implicit none
+
+  public :: kpoint_t
+
+  !< Method specification for non-specified k-list
+  integer, parameter, public :: K_METHOD_NONE = 0
+  !< Method specification for Monkhorst-Pack grid
+  integer, parameter, public :: K_METHOD_MONKHORST_PACK = 1
+  !< Method specification for cutoff designation
+  integer, parameter, public :: K_METHOD_CUTOFF = 2
+  !< Method specification for user-defined list
+  integer, parameter, public :: K_METHOD_LIST = 3
+
+  type :: kpoint_t
+    
+    !< Total number of k-points
+    integer :: N = 0
+    !< K-points
+    real(dp), pointer :: k(:,:) => null()
+    !< weights for k-points
+    real(dp), pointer :: w(:) => null()
+    !< whether these k-points are under time-reversal symmetry
+    logical :: trs = .true.
+    
+    !< The method by which these k-points are generated
+    integer :: method = K_METHOD_MONKHORST_PACK
+
+    ! Specific elements for Monkhorst-Pack, cutoff grids
+    integer :: k_cell(3,3) = 0
+    real(dp) :: k_displ(3) = 0._dp
+
+    ! The cutoff requested
+    real(dp) :: cutoff = 0._dp
+
+  end type kpoint_t
+
+  public :: kpoint_associated
+  public :: kpoint_associate
+  public :: kpoint_nullify
+  public :: kpoint_delete
+  public :: kpoint_read
+  
+  public :: kpoint_read_file
+  public :: kpoint_write_xml
+  public :: kpoint_write_stdout
+  public :: kpoint_write_file
+
+contains
+
+  !< Associate one k-point list with another (no copies made)
+  subroutine kpoint_associate(out, in)
+    type(kpoint_t), intent(inout) :: out
+    type(kpoint_t), intent(in) :: in
+
+    out%N = in%N
+    out%k => in%k
+    out%w => in%w
+    out%trs = in%trs
+    out%method = in%method
+    out%k_cell = in%k_cell
+    out%k_displ = in%k_displ
+
+  end subroutine kpoint_associate
+
+  !< Figure out if two types are associated with each other
+  function kpoint_associated(a, b) result(assoc)
+    type(kpoint_t), intent(in) :: a, b
+    logical :: assoc
+
+    assoc = associated(a%k, b%k) .and. &
+        associated(a%w, b%w)
+
+  end function kpoint_associated
+
+  !< Nullify k-point list
+  subroutine kpoint_nullify(this)
+    type(kpoint_t), intent(inout) :: this
+
+    nullify(this%k, this%w)
+    call kpoint_delete(this)
+
+  end subroutine kpoint_nullify
+
+  !< Delete the k-point list
+  subroutine kpoint_delete(this)
+    use alloc, only: de_alloc
+    type(kpoint_t), intent(inout) :: this
+
+    call de_alloc(this%k, name='kpoint_t%k')
+    call de_alloc(this%w, name='kpoint_t%w')
+    this%N = 0
+    this%trs = .true.
+    this%method = K_METHOD_NONE
+    this%k_cell = 0
+    this%k_displ = 0
+    
+  end subroutine kpoint_delete
+
+  !< Read settings for this k-point grid
+  !!
+  !! The order of reading the k-points are as follows:
+  !!
+  !! 1. kgrid.MonkhorstPack has precedence
+  !!    If this block is found the block will be read as this:
+  !!     %block kgrid.MonkhorstPack  # Defines k_cell and k_displ
+  !!       4  0  0   0.50               # (k_cell(i,1),i=1,3), k_displ(1)
+  !!       0  4  0   0.50               # (k_cell(i,2),i=1,3), k_displ(2)
+  !!       0  0  4   0.50               # (k_cell(i,3),i=1,3), k_displ(3)
+  !!     %endblock kgrid.MonkhorstPack
+  !!    Note that the displacement defaults to zero and is an optional value.
+  !!    One may also write the block like this:
+  !!     %block kgrid.MonkhorstPack  # Defines k_cell and k_displ
+  !!       4 0.50               # k_cell(1,1), k_displ(1)
+  !!       4 0.50               # k_cell(2,2), k_displ(2)
+  !!       4 0.50               # k_cell(3,3), k_displ(3)
+  !!     %endblock kgrid.MonkhorstPack
+  !! 2. kgrid.Cutoff
+  !!    The k_cell is specified using a cutoff parameter and automatically
+  !!    calculated.
+  !! 3. kgrid.file
+  !!    The k-points are user-supplied. In this case the k-points *must*
+  !!    be provided in units of the reciprocal lattice vectors (i.e. in ]-0.5 ; 0.5])
+  subroutine kpoint_read(this, prefix, cell, trs, process_k_cell)
+    use fdf
+    
+    type(kpoint_t), intent(inout) :: this
+    character(len=*), intent(in) :: prefix
+    real(dp), intent(in) :: cell(3,3)
+    logical, intent(in) :: trs
+    abstract interface
+      subroutine my_sub(k_cell, k_displ)
+        use precision, only: dp
+        integer, intent(inout) :: k_cell(3,3)
+        real(dp), intent(inout) :: k_displ(3)
+      end subroutine
+    end interface
+    procedure(my_sub), optional :: process_k_cell
+
+    ! For reading the k-cell etc.
+    type(block_fdf) :: bfdf
+    type(parsed_line), pointer :: pline
+
+    character(len=256) :: name_fdf, file
+    real(dp) :: cutoff
+    integer :: i, nvals
+    logical :: spiral
+
+    ! If this already has k-points associated we continue
+    if ( this%N > 0 ) return
+
+    ! Be sure to delete everything
+    call kpoint_delete(this)
+    this%method = K_METHOD_NONE
+
+    spiral = fdf_get('Spin.Spiral', .false.)
+    ! Allow the user to control the use of time-reversal-symmetry
+    ! By default, it is on, except for "spin-spiral" calculations
+    ! and/or non-collinear calculations
+    this%trs = fdf_get("TimeReversalSymmetryForKpoints", (.not. spiral) .and. trs)
+
+
+    call kpoint_fdf_name(prefix, 'kgrid.MonkhorstPack', name_fdf)
+    if ( fdf_block(name_fdf, bfdf) ) then
+
+      ! The method is a k-point MP grid
+      this%method = K_METHOD_MONKHORST_PACK
+
+      ! Now read data
+      do i = 1, 3
+        if ( .not. fdf_bline(bfdf,pline) ) then
+          call die('kpoint_read: ERROR in ' // trim(name_fdf) // ' block')
+        end if
+
+        ! Read displacement
+        nvals = fdf_bnvalues(pline)
+        if ( nvals > 3 ) then
+          this%k_displ(i) = mod(fdf_bvalues(pline,4), 1._dp)
+        else
+          this%k_displ(i) = 0._dp
+        end if
+
+        if ( nvals == 1 ) then
+          this%k_cell(i,i) = fdf_bintegers(pline,1)
+          
+        else if ( nvals == 2 ) then
+          this%k_cell(i,i) = fdf_bintegers(pline,1)
+          this%k_displ(i) = mod(fdf_bvalues(pline,2), 1._dp)
+          
+        else
+          this%k_cell(1,i) = fdf_bintegers(pline,1)
+          this%k_cell(2,i) = fdf_bintegers(pline,2)
+          this%k_cell(3,i) = fdf_bintegers(pline,3)
+          
+        end if
+        
+      end do
+
+      call setup_mp()
+
+    end if
+
+    ! Quick return
+    if ( this%method /= K_METHOD_NONE ) return
+
+    call kpoint_fdf_name(prefix, 'kgrid.Cutoff', name_fdf)
+    if ( fdf_defined(name_fdf) ) then
+
+      ! The method is cutoff based
+      this%method = K_METHOD_CUTOFF
+      cutoff = fdf_get(name_fdf, -1._dp, 'Bohr')
+      if ( cutoff <= 0._dp ) then
+        this%method = K_METHOD_NONE
+      else
+        call setup_cutoff()
+      end if
+
+    end if
+
+    ! Quick return
+    if ( this%method /= K_METHOD_NONE ) return
+
+    call kpoint_fdf_name(prefix, 'kgrid.File', name_fdf)
+    if ( fdf_defined(name_fdf) ) then
+
+      ! A user-defined list of k-points
+      this%method = K_METHOD_LIST
+      file = fdf_get(name_fdf, '01234567890123456789')
+
+      call setup_file()
+
+    end if
+
+    ! Quick return
+    if ( this%method /= K_METHOD_NONE ) return
+
+    call setup_gamma()
+
+  contains
+
+    subroutine setup_mp()
+      use m_find_kgrid, only: find_kgrid
+      
+      if ( present(process_k_cell) ) then
+        call process_k_cell(this%k_cell, this%k_displ)
+      end if
+
+      ! Find the grid
+      call find_kgrid(cell, this%k_cell, this%k_displ, .true., &
+          this%trs, &
+          this%N, this%k, this%w, this%cutoff)
+      
+    end subroutine setup_mp
+    
+    subroutine setup_cutoff()
+      use m_minvec, only : minvec
+      use m_find_kgrid, only: find_kgrid
+
+      real(dp) :: scmin(3,3), vmod
+      real(dp) :: ctransf(3,3)
+      integer :: factor, expansion_factor
+
+      ! Find equivalent rounded unit-cell
+      call minvec( cell, scmin, ctransf )
+      
+      expansion_factor = 1
+      do i = 1 , 3
+        
+        vmod = sqrt( dot_product(scmin(:,i),scmin(:,i)) )
+        factor = int(2.0_dp * cutoff / vmod) + 1
+        expansion_factor = expansion_factor * factor
+        
+        ! Generate actual supercell skeleton
+        this%k_cell(:,i) = ctransf(:,i) * factor
+        
+      end do
+
+      ! Avoid confusing permutations, revert to identity
+      if ( expansion_factor == 1 ) then
+        
+        this%k_cell(:,:) = 0
+        do i = 1, 3
+          this%k_cell(i,i) = 1
+        end do
+        
+      end if
+
+      if ( present(process_k_cell) ) then
+        call process_k_cell(this%k_cell, this%k_displ)
+      end if
+
+      ! Find the grid
+      call find_kgrid(cell,this%k_cell, this%k_displ, .false., &
+          this%trs, &
+          this%N, this%k, this%w, this%cutoff)
+      
+    end subroutine setup_cutoff
+
+    subroutine setup_file()
+
+      call kpoint_read_file(file, cell, this)
+
+    end subroutine setup_file
+
+    subroutine setup_gamma()
+      use alloc, only: re_alloc
+
+      ! We have a Gamma-only calculation
+      this%method = K_METHOD_NONE
+      this%N = 1
+      call re_alloc(this%k, 1, 3, 1, 1, name='kpoint%k')
+      call re_alloc(this%w, 1, 1, name='kpoint%w')
+      this%k = 0._dp
+      this%w = 1._dp
+      this%k_cell = 0
+      do i = 1, 3
+        this%k_cell(i,i) = 1
+      end do
+      
+    end subroutine setup_gamma
+
+  end subroutine kpoint_read
+
+  ! The user can specify their own k-points
+  subroutine kpoint_read_file(fname, cell, this)
+    use alloc, only: re_alloc
+    use parallel, only : Node
+    use m_os, only : file_exist
+#ifdef MPI
+    use mpi_siesta
+#endif
+    character(len=*), intent(in) :: fname
+    real(dp), intent(in) :: cell(3,3)
+    type(kpoint_t), intent(inout) :: this
+
+    real(dp) :: rcell(3,3), k(3), wsum
+
+#ifdef MPI
+    integer :: MPIerror
+#endif
+
+    ! The user has requested to read in 
+    ! k-points from a specific file...
+    ! We will do that for him/her.
+    integer :: iu, ik, ix, stat
+
+    if ( Node == 0 ) then
+
+      if ( .not. file_exist(trim(fname)) ) then
+        call die('Could not locate file '//trim(fname)// &
+            ' please ensure that the file exists.')
+      end if
+
+      call io_assign( iu )
+      open( iu, file=trim(fname), form='formatted', status='old' )
+
+      ! Read number of k-points
+      read(iu,*,iostat=stat) this%N
+      call kill_iokp(stat,0)
+
+    end if
+
+#ifdef MPI
+    call MPI_Bcast(this%N,1,MPI_Integer,0,MPI_Comm_World,MPIerror)
+#endif
+
+    call re_alloc(this%k, 1, 3, 1, this%N, name='kpoint%k')
+    call re_alloc(this%w, 1, this%N, name='kpoint%w')
+
+    if ( Node == 0 ) then
+
+      call reclat(cell, rcell, 1)
+
+      ! Read in the k-points
+      wsum = 0._dp
+      do ik = 1 , this%N
+        
+        ! read current k-point
+        read(iu,*,iostat=stat) ix, k(:), this%w(ik) ! (i6,3f12.6,3x,f12.6)
+        call kpoint_convert(rcell, k, this%k(:,ik), -2)
+        call kill_iokp(stat,ik)
+        wsum = wsum + this%w(ik)
+
+      end do
+
+      if ( abs(wsum - 1._dp) > 1.e-7_dp ) then
+        write(*,'(a)')'WARNING: Weights for user specified k-points does &
+            &not sum to 1.'
+        call die('User specified k-points does not sum to 1.')
+      end if
+
+      call io_close( iu )
+
+    end if
+
+#ifdef MPI
+    call MPI_Bcast(this%k(1,1),3*this%N,MPI_Double_Precision,0,MPI_Comm_World,MPIerror)
+    call MPI_Bcast(this%w(1),this%N,MPI_Double_Precision,0,MPI_Comm_World,MPIerror)
+#endif
+
+  contains
+
+    subroutine kill_iokp(stat, line)
+      integer, intent(in) :: stat, line
+
+      if ( stat == 0 ) return
+
+      write(*,*) 'Siesta kpoint_read could not read your input file'
+      write(*,*) 'The k-points MUST be in units of reciprocal vectors!'
+      write(*,*) 'Siesta will convert the unit to correct units.'
+      write(*,*) 'Also the sum of weights MUST equal 1.'
+      write(*,*) !
+      if ( line == 0 ) then
+        write(*,*) 'Error occured on reading number of k-points (first line)'
+      else
+        write(*,'(a,i0,a)') 'Error occured on reading the ',line,' kpoint.'
+      end if
+      write(*,*) 'Please format your file like this:'
+      write(*,*) ' $> cat '//trim(fname)
+      write(*,*) ' <nkpt>'
+      write(*,*) '     1  <kpt-A1> <kpt-A2> <kpt-A3> <w-kpt>'
+      write(*,*) '     2  <kpt-A1> <kpt-A2> <kpt-A3> <w-kpt>'
+      write(*,*) ' ....'
+      write(*,*) ' <nkpt> <kpt-A1> <kpt-A2> <kpt-A3> <w-kpt>'
+
+      call die('Siesta reading user specified k-points')
+
+    end subroutine kill_iokp
+
+  end subroutine kpoint_read_file
+
+
+  !< Write k-point list to the XML output file
+  !!
+  !! This will populate a property list with all information about the
+  !! k-point generation.
+  subroutine kpoint_write_xml(this, prefix)
+    use siesta_cml
+    use m_char, only: lcase
+    use units, only: Ang
+    
+    type(kpoint_t), intent(inout) :: this
+    character(len=*), intent(in), optional :: prefix
+    character(len=64) :: lcase_prefix
+    integer :: ik
+
+    ! Quick return, if able
+    if ( .not. cml_p ) return
+
+    if ( present(prefix) ) then
+      lcase_prefix = lcase(prefix)
+      
+      call cmlStartPropertyList(xf=mainXML, title=trim(prefix)//".k-points", &
+          dictRef="siesta:"//trim(lcase_prefix)//".kpoints")
+      call cmlAddProperty(xf=mainXML, value=this%N, &
+          dictref='siesta:'//trim(lcase_prefix)//'.nkpnt', &
+          units="cmlUnits:countable")
+      do ik = 1, this%N
+        call cmlAddKPoint(xf=mainXML, coords=this%k(:,ik), weight=this%w(ik))
+      end do
+
+      ! Add supplementary information for MP and Cutoff methods
+      select case ( this%method )
+      case ( K_METHOD_MONKHORST_PACK, K_METHOD_CUTOFF )
+        call cmlAddProperty(xf=mainXML, value=this%cutoff/Ang, & 
+            dictref='siesta:'//trim(lcase_prefix)//'.kcutoff', units='siestaUnits:angstrom')
+        call cmlAddProperty(xf=mainXML, value=this%k_cell, &
+            dictref='siesta:'//trim(lcase_prefix)//'.kscell', units="cmlUnits:countable")
+        call cmlAddProperty(xf=mainXML, value=this%k_displ, &
+            dictref='siesta:'//trim(lcase_prefix)//'.kdispl')
+      end select
+           
+    else
+      call cmlStartPropertyList(xf=mainXML, title="k-points", &
+          dictRef="siesta:kpoints")
+      call cmlAddProperty(xf=mainXML, value=this%N, dictref='siesta:nkpnt', &
+          units="cmlUnits:countable")
+      do ik = 1, this%N
+        call cmlAddKPoint(xf=mainXML, coords=this%k(:,ik), weight=this%w(ik))
+      end do
+
+      ! Add supplementary information for MP and Cutoff methods
+      select case ( this%method )
+      case ( K_METHOD_MONKHORST_PACK, K_METHOD_CUTOFF )
+        call cmlAddProperty(xf=mainXML, value=this%cutoff/Ang, &
+            dictref='siesta:kcutoff', units='siestaUnits:angstrom')
+        call cmlAddProperty(xf=mainXML, value=this%k_cell, &
+            dictref='siesta:kscell', units="cmlUnits:countable")
+        call cmlAddProperty(xf=mainXML, value=this%k_displ, &
+            dictref='siesta:kdispl')
+      end select
+
+    end if
+
+    call cmlEndPropertyList(mainXML)
+
+  end subroutine kpoint_write_xml
+  
+  !< Write to std-out the k-points and some information regarding the generation of the k-list
+  !!
+  !! The k-points are only written if `all` is `.true.`.
+  !! Otherwise only information regarding the generation will be written.
+  subroutine kpoint_write_stdout(this, all, prefix)
+    use units, only: Ang
+    type(kpoint_t), intent(in) :: this
+    logical, intent(in) :: all
+    character(len=*), intent(in), optional :: prefix
+    character(len=64) :: name
+    integer :: ik
+
+    if ( present(prefix) ) then
+      name = trim(prefix) // ' k-'
+    else
+      name = 'k-'
+    end if
+
+    if ( all ) then
+      write(*,'(/,3a)') 'siesta: ',trim(name), 'point coordinates (Bohr**-1) and weights:'
+      do ik = 1, this%N
+        write(*,'(a,i10,3(tr1,e13.6),tr3,e12.6)') 'siesta: ', ik, this%k(:,ik), this%w(ik)
+      end do
+    end if
+    write(*,'(/3a,i10)')  'siesta: ', trim(name), 'grid: Number of k-points =', this%N
+
+    select case ( this%method )
+    case ( K_METHOD_NONE )
+      write(*,'(3a)')  'siesta: ', trim(name), 'point is Gamma-only'
+    case ( K_METHOD_MONKHORST_PACK )
+      write(*,'(3a)')  'siesta: ', trim(name), 'points from Monkhorst-Pack grid'
+      write(*,'(3a,f10.3,a)')  'siesta: ', trim(name), 'cutoff (effective) =', this%cutoff/Ang, ' Ang'
+      write(*,'(3a)') 'siesta: ', trim(name), 'point supercell and displacements'
+      do ik = 1, 3
+        write(*,'(a,3i4,3x,f8.3)') 'siesta: k-grid: ', this%k_cell(:,ik), this%k_displ(ik)
+      end do
+    case ( K_METHOD_CUTOFF )
+      write(*,'(3a)')  'siesta: ', trim(name), 'points from cutoff'
+      write(*,'(3a,f10.3,a)')  'siesta: ', trim(name), 'cutoff (effective) =', this%cutoff/Ang, ' Ang'
+      write(*,'(3a)') 'siesta: ', trim(name), 'point supercell and displacements'
+      do ik = 1, 3
+        write(*,'(a,3i4,3x,f8.3)') 'siesta: k-grid: ', this%k_cell(:,ik), this%k_displ(ik)
+      end do
+    case ( K_METHOD_LIST )
+      write(*,'(3a)')  'siesta: ', trim(name), 'points from user-defined list'
+    end select
+    
+  end subroutine kpoint_write_stdout
+
+  !< Write to a file the k-point information
+  subroutine kpoint_write_file(this, suffix)
+    use files, only: slabel
+    use m_io, only: io_assign, io_close
+    type(kpoint_t), intent(in) :: this
+    character(len=*), intent(in) :: suffix
+    character(len=256) :: fname
+    integer :: iu, ik
+
+    fname = trim(slabel) // '.' // trim(suffix)
+
+    call io_assign( iu )
+    open( iu, file=fname, form='formatted', status='unknown' )      
+    
+    write(iu,'(i10)') this%N
+    do ik = 1, this%N
+      write(iu,'(i10,3(tr1,e13.6),tr3,e12.6)') ik, this%k(:,ik), this%w(ik)
+    end do
+      
+    call io_close( iu )
+    
+  end subroutine kpoint_write_file
+  
+  subroutine kpoint_fdf_name(prefix, suffix, name)
+    character(len=*), intent(in) :: prefix, suffix
+    character(len=256), intent(out) :: name
+
+    if ( len_trim(prefix) > 0 ) then
+      name = trim(prefix) // '.' // trim(suffix)
+    else
+      name = trim(suffix)
+    end if
+
+  end subroutine kpoint_fdf_name
+
+end module kpoint_t_m
+
+
+
+
+
 
=== modified file 'Src/local_DOS.F'
--- Src/local_DOS.F	2018-02-27 14:03:49 +0000
+++ Src/local_DOS.F	2018-05-27 11:17:14 +0000
@@ -26,7 +26,7 @@
       use sys,            only: die                      
       use files,          only: slabel     ! system label
       use files,          only: filesOut_t ! derived type for output file names
-      use Kpoint_grid
+      use kpoint_scf_m, only: kpoints_scf
       use parallel,       only: IOnode                   
       use files,          only : label_length            
       use m_ntm
@@ -79,7 +79,8 @@
           call diagon(no_s, spinor_dim, no_l, maxnh, maxnh, no_u,
      .                numh, listhptr, listh, numh, listhptr, listh, 
      .                H, S, qtot, fixspin, qtots, temp, e1, e2,
-     .                gamma, xijo, indxuo, nkpnt, kpoint, kweight,
+     .                gamma, xijo, indxuo,
+     .                kpoints_scf%N, kpoints_scf%k, kpoints_scf%w,
      .                eo, qo, Dscf, Escf, ef, efs, dummy_Entrop, no_u,
      .                occtol, dummy_iscf, neigwanted)
 
 
=== modified file 'Src/m_initwf.F90'
--- Src/m_initwf.F90	2017-11-10 16:16:08 +0000
+++ Src/m_initwf.F90	2018-05-27 11:17:14 +0000
@@ -77,7 +77,7 @@
       use fdf
       use densematrix,     only : Haux, Saux, psi
       use sparse_matrices, only : maxnh
-      use Kpoint_grid,     only : nkpnt, kpoint, kweight
+      use kpoint_scf_m,    only : kpoints_scf
       use atomlist,        only : no_s, no_l, no_u, qtot, indxuo
       use m_spin,          only : nspin
       use m_gamma,         only : gamma
@@ -163,10 +163,10 @@
       call re_alloc(psi,1,npsi,name='psi',routine='initwf')
       allocate(muo(nuo),stat=mem_stat)
       call memory('A','I',nuo,'initwf',stat=mem_stat)
-      allocate(nocck(nkpnt,nspin),stat=mem_stat)
-      call memory('A','I',nkpnt*nspin,'initwf',stat=mem_stat)
-      allocate(occup(no_u,nspin,nkpnt),stat=mem_stat)
-      call memory('A','L',nuo*nkpnt*nspin,'initwf',stat=mem_stat)
+      allocate(nocck(kpoints_scf%N,nspin),stat=mem_stat)
+      call memory('A','I',kpoints_scf%N*nspin,'initwf',stat=mem_stat)
+      allocate(occup(no_u,nspin,kpoints_scf%N),stat=mem_stat)
+      call memory('A','L',nuo*kpoints_scf%N*nspin,'initwf',stat=mem_stat)
 !     Check indxuo .......................................................
       do iuo = 1,nuo
         muo(iuo) = 0
@@ -199,7 +199,7 @@
 !     evolved by integrating TDKS equations.                                  !
 ! ............................................................................!
       temp=1.0d-6
-      call fermid( nspin, nspin, nkpnt, kweight, no_u, no_u, eo, &
+      call fermid( nspin, nspin, kpoints_scf%N, kpoints_scf%w, no_u, no_u, eo, &
                    temp, qtot, qo, ef, entrp )
       nocc(1) = 0
       nocc(2) = 0
@@ -207,26 +207,26 @@
       degen= .false.
       !
       !
-      do ik=1,nkpnt
+      do ik=1,kpoints_scf%N
         do ispin=1,nspin
           nocck(ik,ispin)=0
           do io=1,no_u
             occup(io,ispin,ik)=.false.
-            if(dabs(qo(io,ispin,ik)-2.0d0*kweight(ik)/nspin).le.    &
-               1.0d-2*dabs(2.0d0*kweight(ik)/nspin))  then
+            if(dabs(qo(io,ispin,ik)-2.0d0*kpoints_scf%w(ik)/nspin).le.    &
+               1.0d-2*dabs(2.0d0**kpoints_scf%w(ik)/nspin))  then
               nocc(ispin)=nocc(ispin)+1
               nocck(ik,ispin)=nocck(ik,ispin)+1
 !             Accounting the number of electrons corresponding the states being marked
 !             as occupied.
-              nelect=nelect+dabs(2.0d0*kweight(ik)/nspin)
+              nelect=nelect+dabs(2.0d0*kpoints_scf%w(ik)/nspin)
               occup(io,ispin,ik)=.true.
             else
-              if ( dabs( qo(io,ispin,ik)) .gt.1.0d-2*dabs(2.0d0*kweight(ik)/nspin)) then
+              if ( dabs( qo(io,ispin,ik)) .gt.1.0d-2*dabs(2.0d0*kpoints_scf%w(ik)/nspin)) then
                 IF (Node .eq. 0) THEN
                   IF(.not. degen) write(6,fmt="(/,a,tr3,a,tr3,a,tr3,a)") "initwf:","ik", &
                          "occupancy","maximum occupancy"
                   write(6,"(tr2,I10,tr3,f8.6,tr4,f8.6)") ik, qo(io,ispin,ik), &
-                         2.0d0*kweight(ik)/nspin
+                         2.0d0*kpoints_scf%w(ik)/nspin
                 END IF
                 degen = .true.  
               end if
@@ -259,8 +259,8 @@
 #else
       m_storage='szden'
 #endif
-      allocate(wavef_ms(1:nkpnt,1:nspin)) ! allocate (nkpnt*npsin) matrices inside wavef_ms
-      do i=1,nkpnt !for every value of nkpnt and nspin, allocate a matrix of size (no_u x nocck(i,j))
+      allocate(wavef_ms(1:kpoints_scf%N,1:nspin)) ! allocate (nkpnt*npsin) matrices inside wavef_ms
+      do i=1,kpoints_scf%N !for every value of nkpnt and nspin, allocate a matrix of size (no_u x nocck(i,j))
         do j=1,nspin
           call m_allocate(wavef_ms(i,j),no_u,nocck(i,j),m_storage)
         end do
@@ -271,7 +271,7 @@
                     Haux, Saux, psi, no_u, occup)
       else if (nspin.le.2 .and. .not.gamma) then
           call diagkiwf( nspin, nuo, no_s, nspin, no_l, maxnh,                 &
-                         no_u, indxuo, nkpnt, kpoint, Haux, Saux,                 &
+                         no_u, indxuo, kpoints_scf%N, kpoints_scf%k, Haux, Saux, &
                          psi, no_u, occup)
       else 
          call die('initwf: ERROR: non-collinear spin options for TDDFT not yet implemented')
@@ -279,7 +279,7 @@
 !     Write/save wavefunction in .TDWF file to use for TDDFT calculation.
       IF (Node .eq. 0) WRITE(6,'(a)') 'initwf: Saving wavefunctions & 
                   &in <systemlabel>.TDWF file.'
-      call  iowavef('write',wavef_ms,no_u,nkpnt,nspin)
+      call  iowavef('write',wavef_ms,no_u,kpoints_scf%N,nspin)
 !     Free local arrays
       call memory('D','I',size(muo),'initwf',stat=mem_stat)
       deallocate(muo,stat=mem_stat)
 
=== modified file 'Src/m_ncdf_siesta.F90'
--- Src/m_ncdf_siesta.F90	2018-02-27 14:03:49 +0000
+++ Src/m_ncdf_siesta.F90	2018-05-27 11:17:14 +0000
@@ -491,12 +491,12 @@
 
   subroutine cdf_save_settings(fname)
 
-    use kpoint_grid,    only: kscell, kdispl
+    use kpoint_scf_m,   only: kpoints_scf
     use siesta_options, only: cdf_w_parallel
     use siesta_options, only: dDtol, dHtol, charnet, wmix, temp, g2cut
     use siesta_options, only: isolve
     use siesta_options, only: SOLVE_DIAGON, SOLVE_ORDERN, SOLVE_TRANSI
-    use m_ts_kpoints,   only: ts_kscell, ts_kdispl
+    use ts_kpoint_scf_m, only: ts_kpoints_scf
 
     character(len=*), intent(in) :: fname
     
@@ -509,8 +509,8 @@
     call ncdf_open_grp(ncdf,'SETTINGS',grp)
 
     ! Save settings
-    call ncdf_put_var(grp,'BZ',kscell)
-    call ncdf_put_var(grp,'BZ_displ',kdispl)
+    call ncdf_put_var(grp,'BZ',kpoints_scf%k_cell)
+    call ncdf_put_var(grp,'BZ_displ',kpoints_scf%k_displ)
     call ncdf_put_var(grp,'DMTolerance',dDtol)
     call ncdf_put_var(grp,'HTolerance',dHtol)
     call ncdf_put_var(grp,'NetCharge',charnet)
@@ -524,8 +524,8 @@
     if ( isolve == SOLVE_TRANSI ) then
        call ncdf_open_grp(ncdf,'TRANSIESTA',grp)
 
-       call ncdf_put_var(grp,'BZ',ts_kscell)
-       call ncdf_put_var(grp,'BZ_displ',ts_kdispl)
+       call ncdf_put_var(grp,'BZ',ts_kpoints_scf%k_cell)
+       call ncdf_put_var(grp,'BZ_displ',ts_kpoints_scf%k_displ)
 
     end if
 
 
=== modified file 'Src/m_transiesta.F90'
--- Src/m_transiesta.F90	2018-04-04 13:22:33 +0000
+++ Src/m_transiesta.F90	2018-05-27 11:17:14 +0000
@@ -58,7 +58,7 @@
     use class_OrbitalDistribution
     use class_Sparsity
 
-    use m_ts_kpoints, only : ts_nkpnt, ts_Gamma
+    use ts_kpoint_scf_m, only : ts_kpoints_scf, ts_gamma_scf
 
     use m_ts_electype
 
@@ -145,7 +145,7 @@
        end if
 
        ! print out estimated memory usage...
-       call ts_print_memory(ts_Gamma)
+       call ts_print_memory(ts_gamma_scf)
 
        call ts_print_charges(N_Elec,Elecs, Qtot, sp_dist, sparse_pattern, &
             nspin, n_nzs, DM, S)
@@ -247,7 +247,7 @@
        call open_GF(N_Elec,Elecs,uGF,NEn,.false.)
 
        if ( ts_method == TS_FULL ) then
-          if ( ts_Gamma ) then
+          if ( ts_gamma_scf ) then
              call ts_fullg(N_Elec,Elecs, &
                   nq, uGF, nspin, na_u, lasto, &
                   sp_dist, sparse_pattern, &
@@ -262,7 +262,7 @@
                   H, S, DM, EDM, Ef, DE_NEGF)
           end if
        else if ( ts_method == TS_BTD ) then
-          if ( ts_Gamma ) then
+          if ( ts_gamma_scf ) then
              call ts_trig(N_Elec,Elecs, &
                   nq, uGF, nspin, na_u, lasto, &
                   sp_dist, sparse_pattern, &
@@ -278,7 +278,7 @@
           end if
 #ifdef SIESTA__MUMPS
        else if ( ts_method == TS_MUMPS ) then
-          if ( ts_Gamma ) then
+          if ( ts_gamma_scf ) then
              call ts_mumpsg(N_Elec,Elecs, &
                   nq, uGF, nspin, na_u, lasto, &
                   sp_dist, sparse_pattern, &
@@ -335,7 +335,7 @@
                 Ef = Ef + 0.01_dp * eV
              end if
           else if ( ts_method == TS_BTD ) then
-             if ( ts_Gamma ) then
+             if ( ts_gamma_scf ) then
                 call ts_trig_Fermi(N_Elec,Elecs, &
                      nq, uGF, nspin, na_u, lasto, &
                      sp_dist, sparse_pattern, &
@@ -588,7 +588,7 @@
          end if
 
          if ( Elecs(iEl)%out_of_core ) then
-            call read_Green(uGF(iEl),Elecs(iEl), ts_nkpnt, NEn )
+            call read_Green(uGF(iEl),Elecs(iEl), ts_kpoints_scf%N, NEn )
          end if
          
       end do
 
=== modified file 'Src/m_ts_fullk.F90'
--- Src/m_ts_fullk.F90	2017-12-12 09:13:49 +0000
+++ Src/m_ts_fullk.F90	2018-05-27 11:17:14 +0000
@@ -96,7 +96,7 @@
     ! Self-energy expansion
     use m_ts_elec_se
 
-    use m_ts_kpoints, only : ts_nkpnt, ts_kpoint, ts_kweight
+    use ts_kpoint_scf_m, only : ts_kpoints_scf
 
     use m_ts_options, only : Calc_Forces
     use m_ts_options, only : N_mu, mus
@@ -282,7 +282,7 @@
     end if
 
     ! start the itterators
-    call itt_init  (SpKp,end1=nspin,end2=ts_nkpnt)
+    call itt_init  (SpKp,end1=nspin,end2=ts_kpoints_scf%N)
     ! point to the index iterators
     call itt_attach(SpKp,cur1=ispin,cur2=ikpt)
 
@@ -300,10 +300,10 @@
        end if
 
        ! Include spin factor and 1/(2\pi)
-       kpt(:) = ts_kpoint(:,ikpt)
+       kpt(:) = ts_kpoints_scf%k(:,ikpt)
        ! create the k-point in reciprocal space
        call kpoint_convert(ucell,kpt,bkpt,1)
-       kw = 0.5_dp / Pi * ts_kweight(ikpt)
+       kw = 0.5_dp / Pi * ts_kpoints_scf%w(ikpt)
        if ( nspin == 1 ) kw = kw * 2._dp
 
 #ifdef TRANSIESTA_TIMING
 
=== removed file 'Src/m_ts_kpoints.F90'
--- Src/m_ts_kpoints.F90	2017-11-23 14:17:27 +0000
+++ Src/m_ts_kpoints.F90	1970-01-01 00:00:00 +0000
@@ -1,323 +0,0 @@
-! ---
-! Copyright (C) 1996-2016	The SIESTA group
-!  This file is distributed under the terms of the
-!  GNU General Public License: see COPYING in the top directory
-!  or http://www.gnu.org/copyleft/gpl.txt .
-! See Docs/Contributors.txt for a list of contributors.
-! ---
-module m_ts_kpoints
-!
-! Routines that are related to TS kpoint sampling
-!
-!==============================================================================
-! CONTAINS:
-!          1) setup_ts_scf_kscell
-!          2) setup_ts_kpoint_grid
-!          3) ts_write_k_points
-!          4) ts_iokp
-  
-  use precision, only : dp
-
-  implicit none
-
-  private
-  save
-
-!===================== K-POINT RELATED VARIABLES ============================== 
-!==============================================================================
-! Contains data structures and routines to deal with the kpoint-grid
-! for the self-consistent calculation with Green Functions.
-! Original verison modified so that the kpoint sampling is done for the
-! direction perpendicular to the transport directins
-! Version created to integrate TRANSIESTA in siesta2.3
-!==============================================================================
-! The kpoint mesh parameters that are public may be accessed in other 
-! parts of the code. With respesct to the original names of the variables used 
-! in the m_kpoint_grid module, a "ts_" was added. Also in the SIESTA module
-! kscell and kdispl are not public, but are made public (ts_kscell and
-! ts_kdispl) here.
-
-!==============================================================================
-! DETAILS: To obtain the kpoints for the GFs calculations it uses the same 
-! scheme as for SIESTA but puts ts_kscell(3,3)=1 and ts_kdispl(3)=0.0
-!
-!==============================================================================
- 
-  logical, public  :: ts_scf_kgrid_first_time = .true.
-  logical, public  :: ts_Gamma
-  integer, public  :: ts_nkpnt             ! Total number of k-points
-  real(dp), public :: ts_eff_kgrid_cutoff  ! Effective kgrid_cutoff
-
-  real(dp), pointer, public :: ts_kweight(:) 
-  real(dp), pointer, public :: ts_kpoint(:,:)
-
-  integer,  public :: ts_kscell(3,3) = 0
-  real(dp), public :: ts_kdispl(3) = 0.0_dp
-
-  logical, public :: ts_firm_displ = .false.
-  logical, public :: ts_user_requested_mp = .false.
-  logical, public :: ts_user_requested_cutoff = .false.
-
-  public :: setup_ts_scf_kscell, setup_ts_kpoint_grid
-  public :: ts_write_k_points
-
-contains
-
-  subroutine setup_ts_scf_kscell( cell )
-
-! ***************** INPUT **********************************************
-! real*8  cell(3,3)  : Unit cell vectors in real space cell(ixyz,ivec)
-! type(Elec) Elecs(:) : Electrodes
-
-!   The relevant fdf labels are kgrid_cutoff and kgrid_Monkhorst_Pack.
-!   If both are present, kgrid_Monkhorst_Pack has priority. If none is
-!   present, the cutoff default is zero, producing only the gamma point.
-!   Examples of fdf data specifications:
-!     kgrid_cutoff  50. Bohr
-!     %block kgrid_Monkhorst_Pack  # Defines kscell and kdispl
-!     4  0  0   0.50               # (kscell(i,1),i=1,3), kdispl(1)
-!     0  4  0   0.50               # (kscell(i,2),i=1,3), kdispl(2)
-!     0  0  4   0.50               # (kscell(i,3),i=1,3), kdispl(3)
-!     %endblock kgrid_Monkhorst_Pack
-! **********************************************************************
-
-!  Modules
-
-    use parallel,   only : IONode
-    use m_minvec,   only : minvec
-    use fdf
-    use sys,        only : die
-#ifdef MPI
-    use mpi_siesta
-#endif
-
-    use m_ts_tdir, only: ts_tidx
-    use m_ts_global_vars, only : TSmode
-    
-    implicit          none
-
-    ! Passed variables
-    real(dp), intent(in)   :: cell(3,3)
-
-    ! Internal variables
-    integer :: i, j, factor(3,3), expansion_factor
-    real(dp) :: scmin(3,3), vmod, cutoff
-    real(dp) :: ctransf(3,3)
-
-    type(block_fdf)            :: bfdf
-    type(parsed_line), pointer :: pline
-
-    logical :: bool
-    real(dp), parameter :: defcut = 0.0_dp
-
-    bool = fdf_block('TS.kgrid_Monkhorst_Pack',bfdf)
-    if ( .not. bool ) &
-         bool = fdf_block('kgrid_Monkhorst_Pack',bfdf)
-       
-    if ( bool ) then
-       ts_user_requested_mp = .true.
-       do i = 1,3
-          if (fdf_bline(bfdf,pline)) then
-            ts_kscell(1,i) = fdf_bintegers(pline,1)
-            ts_kscell(2,i) = fdf_bintegers(pline,2)
-            ts_kscell(3,i) = fdf_bintegers(pline,3)
-            if ( fdf_bnvalues(pline) > 3 ) then
-              ts_kdispl(i) = mod(fdf_bvalues(pline,4), 1._dp)
-            else
-              ts_kdispl(i) = 0._dp
-            end if
-          else
-             call die( 'setup_ts_scf_kscell: ERROR no data in' // &
-                  'kgrid_Monkhorst_Pack block' )
-          endif
-       enddo
-       ts_firm_displ = .true.
-       
-    else
-
-       if ( IONode .and. TSmode ) then
-          write(*,*) 'WARNING !!!'
-          write(*,*) 'TS kgrid determined first with 3D cell !!!'
-          write(*,*) 'Specifying only cutoff in Electrode AND Scattering calculations might lead to problems !!'
-       end if
-
-       cutoff = fdf_get('kgrid_cutoff',defcut,'Bohr')
-       ts_user_requested_cutoff = (cutoff /= defcut)
-
-       ts_kdispl(1:3) = 0.0_dp  ! Might be changed later
-       ts_firm_displ = .false.
-
-       ! Find equivalent rounded unit-cell
-       call minvec( cell, scmin, ctransf )
-       
-       expansion_factor = 1
-       do j = 1,3
-          factor(j,1:3) = 0
-          vmod = sqrt(dot_product(scmin(1:3,j),scmin(1:3,j)))
-          factor(j,j) = int(2.0_dp*cutoff/vmod) + 1
-          expansion_factor = expansion_factor * factor(j,j)
-       enddo
-       ! Generate actual supercell skeleton
-       ts_kscell = matmul(ctransf, factor)
-
-       if ( expansion_factor == 1 ) then
-          ts_kscell = 0
-          ts_kscell(1,1) = 1
-          ts_kscell(2,2) = 1
-          ts_kscell(3,3) = 1
-       end if
-       
-    end if
-
-    ! In case of TSmode we have a transiesta run.
-    ! This means that we truncate the k-points in the transport direction.
-    ! However, if we are dealing with an electrode calculation we simply allow it
-    ! to have the same cell (the check made will disregard the transport directions k-points)
-    if ( TSmode .and. ts_tidx > 0 ) then
-       i = ts_tidx
-       ts_kscell(:,i) = 0
-       ts_kscell(i,:) = 0
-       ts_kscell(i,i) = 1
-       ts_kdispl(i)   = 0._dp
-    end if
-    
-  end subroutine setup_ts_scf_kscell
-  
-  subroutine setup_ts_kpoint_grid( ucell )
-    
-    ! SIESTA Modules
-    USE precision, only : dp       
-    USE fdf, only       : fdf_get
-    USE m_find_kgrid, only : find_kgrid
-    USE parallel, only  : IONode
-    use m_ts_global_vars, only : TSmode
-    use kpoint_grid
-
-    ! Local Variables
-    real(dp), intent(in) :: ucell(3,3)
-
-    if ( .not. TSMode ) then
-       
-       ! Same as SIESTA (mainly to be able to write the TSHS files).
-       ts_Gamma = gamma_SCF
-       ts_nkpnt = nkpnt
-       ts_eff_kgrid_cutoff = eff_kgrid_cutoff
-       ts_kweight => kweight
-       ts_kpoint => kpoint
-       ts_kscell = kscell
-       ts_kdispl = kdispl
-
-       ! Since this is not transiesta we immediately return
-       ! and then we also skip printing out transiesta information.
-       return
-       
-    end if
-
-    if ( ts_scf_kgrid_first_time ) then
-
-       nullify(ts_kweight,ts_kpoint)
-       if ( fdf_get('SpinSpiral', .false.) ) then
-          call die('transiesta: Does not work with spin-spiral')
-       end if
-
-       call setup_ts_scf_kscell(ucell)
-
-       ts_scf_kgrid_first_time = .false.
-
-    else
-       if ( ts_user_requested_mp    ) then
-          ! no need to set up the kscell again
-       else
-          call setup_ts_scf_kscell(ucell)
-       endif
-    endif
-
-    call find_kgrid(ucell,ts_kscell,ts_kdispl,ts_firm_displ, &
-         .true., &
-         ts_nkpnt,ts_kpoint,ts_kweight, ts_eff_kgrid_cutoff)
-
-    ts_Gamma = ts_nkpnt == 1 .and. &
-         dot_product(ts_kpoint(:,1),ts_kpoint(:,1)) < 1.0e-20_dp
-
-    if (IONode) call ts_write_k_points()
-
-  end subroutine setup_ts_kpoint_grid
-
-  subroutine ts_write_k_points()
-    USE siesta_options, only: writek
-    USE siesta_cml
-
-    implicit none
-
-    integer  :: ik, ix, i
-
-    if ( writek ) then
-       write(*,'(/,a)') 'transiesta: ts_k-point coordinates (Bohr**-1) and weights:'
-       write(*,'(a,i4,3f12.6,3x,f12.6)')                          &
-            ('transiesta: ', ik, (ts_kpoint(ix,ik),ix=1,3), ts_kweight(ik), &
-            ik=1,ts_nkpnt)
-    endif
-
-    ! Always write the TranSIESTA k-points
-    call ts_iokp( ts_nkpnt, ts_kpoint, ts_kweight )
-
-    write(*,'(/,a,i6)')  'transiesta: k-grid: Number of Green function k-points =', ts_nkpnt
-    write(*,'(a)') 'transiesta: k-grid: Supercell and displacements'
-    write(*,'(a,3i4,3x,f8.3)') 'transiesta: k-grid: ',        &
-         (ts_kscell(i,1),i=1,3), ts_kdispl(1)
-    write(*,'(a,3i4,3x,f8.3)') 'transiesta: k-grid: ',        &
-         (ts_kscell(i,2),i=1,3), ts_kdispl(2)
-    write(*,'(a,3i4,3x,f8.3)') 'transiesta: k-grid: ',        &
-         (ts_kscell(i,3),i=1,3), ts_kdispl(3)
-    if (cml_p) then
-       call cmlStartPropertyList(xf=mainXML, title="Transiesta k-points", &
-            dictRef="siesta:ts_kpoints")
-       call cmlAddProperty(xf=mainXML, value=ts_nkpnt, dictref='siesta:ts_nkpnt', &
-            units="cmlUnits:countable")
-       do ik=1, ts_nkpnt
-          call cmlAddKPoint(xf=mainXML, coords=ts_kpoint(:,ik), weight=ts_kweight(ik))
-       enddo
-       call cmlEndPropertyList(mainXML)
-       call cmlAddProperty(xf=mainXML, value=ts_kscell, dictref='siesta:ts_kscell', &
-            units="siestaUnits:Ang")
-       call cmlAddProperty(xf=mainXML, value=ts_kdispl, dictref='siesta:ts_kdispl', &
-            units="siestaUnits:Ang")
-    endif
-  end subroutine ts_write_k_points
-
-  subroutine ts_iokp( nk, points, weight )
-! *******************************************************************
-! Saves TranSIESTA k-points (only writing) Bohr^-1
-! Emilio Artacho, Feb. 1999
-! Modified by Nick Papior Andersen to not overwrite the SIESTA k-points
-! ********** INPUT **************************************************
-! integer nk           : Number of TS k-points
-! real*8  points(3,nk) : TS k-point coordinates
-! real*8  weight(3,nk) : TS k-point weight
-! *******************************************************************
-    use fdf
-    use files,     only : slabel, label_length
-
-    integer  :: nk
-    real(dp) :: points(3,nk), weight(nk)
-    external :: io_assign, io_close
-
-! Internal 
-    character(len=label_length+5) :: fname
-    integer                       :: iu, ik, ix
-! -------------------------------------------------------------------
-
-    fname = trim( slabel ) // '.TSKP'
-
-    call io_assign( iu )
-    open( iu, file=fname, form='formatted', status='unknown' )      
-
-    write(iu,'(i6)') nk
-    write(iu,'(i6,3f12.6,3x,f12.6)') &
-         (ik, (points(ix,ik),ix=1,3), weight(ik), ik=1,nk)
-
-    call io_close( iu )
-
-  end subroutine ts_iokp
-
-end module m_ts_kpoints
 
=== modified file 'Src/m_ts_mumpsk.F90'
--- Src/m_ts_mumpsk.F90	2017-12-12 08:11:00 +0000
+++ Src/m_ts_mumpsk.F90	2018-05-27 11:17:14 +0000
@@ -66,7 +66,7 @@
     ! Self-energy expansion
     use m_ts_elec_se
 
-    use m_ts_kpoints, only : ts_nkpnt, ts_kpoint, ts_kweight
+    use ts_kpoint_scf_m, only : ts_kpoints_scf
 
     use m_ts_options, only : Calc_Forces
     use m_ts_options, only : N_mu, mus
@@ -236,7 +236,7 @@
     end if
 
     ! start the itterators
-    call itt_init  (SpKp,end1=nspin,end2=ts_nkpnt)
+    call itt_init  (SpKp,end1=nspin,end2=ts_kpoints_scf%N)
     ! point to the index iterators
     call itt_attach(SpKp,cur1=ispin,cur2=ikpt)
 
@@ -252,10 +252,10 @@
        end if
 
        ! Include spin factor and 1/(2\pi)
-       kpt(:) = ts_kpoint(:,ikpt)
+       kpt(:) = ts_kpoints_scf%k(:,ikpt)
        ! create the k-point in reciprocal space
        call kpoint_convert(ucell,kpt,bkpt,1)
-       kw = 0.5_dp / Pi * ts_kweight(ikpt)
+       kw = 0.5_dp / Pi * ts_kpoints_scf%w(ikpt)
        if ( nspin == 1 ) kw = kw * 2._dp
 
        write(mum%ICNTL(1),'(/,a,i0,a,3(tr1,g10.4),/)') &
 
=== modified file 'Src/m_ts_trik.F90'
--- Src/m_ts_trik.F90	2017-11-17 20:53:23 +0000
+++ Src/m_ts_trik.F90	2018-05-27 11:17:14 +0000
@@ -72,7 +72,7 @@
     ! Self-energy expansion
     use m_ts_elec_se
 
-    use m_ts_kpoints, only : ts_nkpnt, ts_kpoint, ts_kweight
+    use ts_kpoint_scf_m, only : ts_kpoints_scf
 
     use m_ts_options, only : Calc_Forces
     use m_ts_options, only : N_mu, mus
@@ -298,7 +298,7 @@
 
     
     ! start the itterators
-    call itt_init  (SpKp,end1=nspin,end2=ts_nkpnt)
+    call itt_init  (SpKp,end1=nspin,end2=ts_kpoints_scf%N)
     ! point to the index iterators
     call itt_attach(SpKp,cur1=ispin,cur2=ikpt)
     
@@ -311,10 +311,10 @@
        end if
 
        ! Include spin factor and 1/(2\pi)
-       kpt(:) = ts_kpoint(:,ikpt)
+       kpt(:) = ts_kpoints_scf%k(:,ikpt)
        ! create the k-point in reciprocal space
        call kpoint_convert(ucell,kpt,bkpt,1)
-       kw = 0.5_dp / Pi * ts_kweight(ikpt)
+       kw = 0.5_dp / Pi * ts_kpoints_scf%w(ikpt)
        if ( nspin == 1 ) kw = kw * 2._dp
        
 #ifdef TRANSIESTA_TIMING
@@ -685,7 +685,7 @@
     ! Self-energy expansion
     use m_ts_elec_se
 
-    use m_ts_kpoints, only : ts_nkpnt, ts_kpoint, ts_kweight
+    use ts_kpoint_scf_m, only : ts_kpoints_scf
 
     use m_ts_sparse, only : ts_sp_uc, tsup_sp_uc
     use m_ts_sparse, only : ltsup_sp_sc, sc_off
@@ -811,15 +811,15 @@
     zDM => val(spuDM)
     call newdSpData2D(ltsup_sp_sc,1, sp_dist,spDM   ,name='TS spDM')
     
-    call itt_init  (SpKp,end1=nspin,end2=ts_nkpnt)
+    call itt_init  (SpKp,end1=nspin,end2=ts_kpoints_scf%N)
     call itt_attach(SpKp,cur1=ispin,cur2=ikpt)
     
     call init_val(spDM)
     do while ( .not. itt_step(SpKp) )
 
-       kpt(:) = ts_kpoint(:,ikpt)
+       kpt(:) = ts_kpoints_scf%k(:,ikpt)
        call kpoint_convert(ucell,kpt,bkpt,1)
-       kw = 0.5_dp / Pi * ts_kweight(ikpt)
+       kw = 0.5_dp / Pi * ts_kpoints_scf%w(ikpt)
        if ( nspin == 1 ) kw = kw * 2._dp
        
 #ifdef TRANSIESTA_TIMING
 
=== modified file 'Src/post_scf_work.F'
--- Src/post_scf_work.F	2018-02-27 14:03:49 +0000
+++ Src/post_scf_work.F	2018-05-27 11:17:14 +0000
@@ -40,7 +40,7 @@
       use m_steps,        only : istp, inicoor
       use m_compute_dm,   only : PreviousCallDiagon
       use m_eo
-      use Kpoint_grid
+      use kpoint_scf_m, only: kpoints_scf
       use m_gamma
       implicit none
 
@@ -81,7 +81,8 @@
      &                no_l, maxnh, maxnh, no_u,
      &                numh, listhptr, listh, numh, listhptr, listh, 
      &                H, S, qtot, fixspin, qtots, temp, 1.0_dp, -1.0_dp,
-     &                gamma, xijo, indxuo, nkpnt, kpoint, kweight,
+     &                gamma, xijo, indxuo,
+     &                kpoints_scf%N, kpoints_scf%k, kpoints_scf%w,
      &                eo, qo, Dscf, Escf, ef, efs, Entropy, no_u,
      &                occtol, iscf, neigwanted)
           Ecorrec = 0.0_dp
@@ -92,7 +93,8 @@
           else
             call zminim(.true., .false., iscf, istp, no_l, nspin, no_u,
      &                  maxnh, numh, listhptr, listh, Escf, eta, qtots,
-     &                  no_s, xijo, indxuo, nkpnt, kpoint, kweight)
+     &                  no_s, xijo, indxuo,
+     &                  kpoints_scf%N, kpoints_scf%k, kpoints_scf%w)
           end if
         endif
       endif
 
=== modified file 'Src/projected_DOS.F'
--- Src/projected_DOS.F	2018-01-04 18:48:03 +0000
+++ Src/projected_DOS.F	2018-05-27 11:17:14 +0000
@@ -15,15 +15,15 @@
 
       public :: init_projected_DOS, projected_DOS
 
-      logical  :: different_pdos_grid    ! Indicates if the grid is the same as the SCF one or not
-
       CONTAINS
       
       subroutine init_projected_DOS()
 
       USE siesta_options
       use fdf,         only: fdf_block, block_fdf
-      use Kpoint_pdos
+      ! This is to get the reference kpoints in case PDOS.kgrid* has not
+      ! been specified
+      use kpoint_pdos_m, only: setup_kpoint_pdos
       use parallel,    only: IOnode, Nodes
       use siesta_geom, only: ucell
       use m_spin,      only: spin
@@ -50,7 +50,7 @@
 
 !     Determine whether the projected density of states is to be computed
 !     on a different grid to the SCF calculation
-      call setup_Kpoint_pdos( ucell, different_pdos_grid )
+      call setup_kpoint_pdos( ucell )
       
 !---------------------------------------------------------------------------END
       
@@ -64,8 +64,8 @@
       use atomlist,    only : indxuo, no_s, no_u, no_l
       use fdf
       use sys,         only : die
-      use Kpoint_grid
-      use Kpoint_pdos
+      use kpoint_scf_m, only: kpoints_scf
+      use Kpoint_pdos_m, only: kpoints_pdos, gamma_pdos
       use parallel,    only: IOnode
       use m_eo
       use m_spin,      only: h_spin_dim, spinor_dim
@@ -108,32 +108,19 @@
          write(6,'(a,3(f8.2,a),2x,i5)')
      $        'siesta: e1, e2, sigma, nhist: ',
      $        e1/eV,' eV',e2/eV,' eV',sigma/eV,' eV', nhist
+       end if
+
+      ! If the k points have been set specifically for the PDOS then use this set
+      if ( kpoints_pdos%N > kpoints_scf%N ) then
+        call re_alloc(eo,1,no_u,1,spinor_dim,1,kpoints_pdos%N,name="eo",
+     .      routine="projected_dos")
       end if
-
-! If the k points have been set specifically for the PDOS then use this set
-      if (different_pdos_grid) then
-
-! If the number of k points has increased then reallocate eo and qo
-         if (maxk_pdos.gt.nkpnt) then
-            call re_alloc(eo,1,no_u,1,spinor_dim,1,maxk_pdos,name="eo",
-     .           routine="projected_dos")
-         endif
-         
-         call pdos( no_s, h_spin_dim, spinor_dim, no_l, 
-     .        maxnh,
-     .        no_u, numh, listhptr, listh, H, S,
-     .        e1, e2, sigma, nhist, gamma_pdos, xijo, indxuo, 
-     .        nkpnt_pdos, kpoints_pdos, kweight_pdos, eo,
-     .        no_u)
-      else
-! otherwise use the SCF grid 
-         call pdos( no_s, h_spin_dim, spinor_dim, no_l, 
-     .        maxnh,
-     .        no_u, numh, listhptr, listh, H, S,
-     .        e1, e2, sigma, nhist,
-     .        gamma, xijo, indxuo, nkpnt, kpoint, kweight, eo,
-     .        no_u)
-      endif
+      
+      call pdos( no_s, h_spin_dim, spinor_dim, no_l, 
+     .    maxnh, no_u, numh, listhptr, listh, H, S,
+     .    e1, e2, sigma, nhist, gamma_pdos, xijo, indxuo, 
+     .    kpoints_pdos%N, kpoints_pdos%k, kpoints_pdos%w, eo,
+     .    no_u)
 
       end subroutine projected_DOS
 
 
=== modified file 'Src/sankey_change_basis.F90'
--- Src/sankey_change_basis.F90	2017-08-30 14:09:10 +0000
+++ Src/sankey_change_basis.F90	2018-05-27 11:17:14 +0000
@@ -1,236 +1,236 @@
-!
-! Copyright (C) 1996-2016 The SIESTA group
-!  This file is distributed under the terms of the
-!  GNU General Public License: see COPYING in the top directory
-!  or http://www.gnu.org/copyleft/gpl.txt.
-! See Docs/Contributors.txt for a list of contributors.
-!
-MODULE m_sankey_change_basis 
-
-  IMPLICIT NONE
-  PRIVATE
-
-  PUBLIC :: sankey_change_basis
-
-  CONTAINS
-! *******************************************************************
-
-     SUBROUTINE sankey_change_basis ( istpmove )
-
-!******************************************************************************
-! This subroutine transforms the TDKS wavefunctions into new-basis during
-! Ehrenfest dynamics within TDDFT after atomic positions are changed. The 
-! TDKS are transformed according to the transformation proposed by Tomfohr and 
-! Sankey in Ref. phys. stat. sol. (b) 226 No.1, 115-123 (2001). 
-! 
-! After transforming TDKS states it calculates the density matrix in new basis set. 
-! We use MatrixSwitch to manipulate matrices.
-!
-! This subroutine is roughly based on D. Sanchez-Portal's chgbasis subroutine from, 
-! now absolete, serial version of tddft-siesta.
-!
-! Re-written for parallelization by Adiran Garaizar and Rafi Ullah, October 2015
-! Modified by Rafi Ullah, February 2017
-!********************************************************************************
-
-  use precision
-  use parallel,            only : Node, Nodes,BlockSize, IONode 
-  use parallelsubs,        only : GlobalToLocalOrb, GetNodeOrbs,               &
-                                  LocalToGlobalOrb
-  use m_diag_option,       only : ParallelOverK
-  use fdf
-  use alloc
-  use sys, only: die
-#ifdef MPI
-  use mpi_siesta,          only : mpi_bcast, mpi_comm_world, mpi_logical
-#endif
-  use m_spin,              only: nspin
-  use m_gamma,             only: gamma
-  use Kpoint_grid,         only: nkpnt, kpoint, kweight
-  use atomlist,            only: no_u, indxuo
-  use wavefunctions
-  use sparse_matrices,     only : numh, listhptr, listh, S, xijo, Dscf
-  use MatrixSwitch
-  use matdiagon,           only: geteigen 
-  !
-  implicit none
-  !
-  integer, intent(in)     :: istpmove
-  !
-#ifdef MPI
-  integer                 :: MPIerror,desch(9)
-  external                :: diagkp
-#endif
-  !
-  logical, save           :: frstme = .true.
-  integer                 :: io, iuo, juo, nuo, jo, ind, ispin
-  integer                 :: ik, j,ierror
-  real(dp)                :: skxij,ckxij, kxij, qe 
-  complex(dp)             :: cvar1
-
-  type(matrix)                     :: Maux,invsqS,phi
-  type(matrix)                     :: Sauxms
-  type(matrix),allocatable,save    :: sqrtS(:)
-  character(3)                     :: m_operation
-  character(5)                     :: m_storage
-  !
-#ifdef DEBUG
-  call write_debug( '    PRE sankey_change_basis' )
-#endif
-  !
-#ifdef MPI
-  call GetNodeOrbs(no_u,Node,Nodes,nuo)
-  if (frstme) then
-    if(ParallelOverK) then
-      call die ( "chgbasis: tddft-siesta not parallelized over k-points." )
-    endif
-      call ms_scalapack_setup(mpi_comm_world,1,'c',BlockSize)
-  endif
-#else
-  Node = 0
-  Nodes = 1
-  nuo = no_u
-#endif
-  call timer( 'sankey_change_basis', 1 )
-  !
-#ifdef MPI
-  m_storage='pzdbc'
-  m_operation='lap'
-#else
-  m_storage='szden'
-  m_operation='lap'
-#endif
-  !
-  IF (nspin .eq. 4) THEN
-    call die ('chgbasis: ERROR: EID not yet prepared for non-collinear spin')
-  END IF
-    ! Allocate local arrays
-  if(frstme) then
-    allocate(sqrtS(nkpnt))
-    do ik=1,nkpnt
-      call m_allocate(sqrtS(ik),no_u,no_u,m_storage)
-    end do
-    frstme=.false.
-  endif
-  call m_allocate(Maux,no_u,no_u,m_storage)
-  call m_allocate(invsqS,no_u,no_u,m_storage)
-  !
-  do ik = 1,nkpnt
-  call m_allocate(Sauxms,no_u,no_u,m_storage)
-    do iuo = 1,nuo
-      call LocalToGlobalOrb(iuo, Node, Nodes, io)
-      do j = 1,numh(iuo)
-        ind = listhptr(iuo) + j
-        juo = listh(ind)
-        jo  = indxuo (juo)
-        if(.not.gamma) then 
-          kxij = kpoint(1,ik) * xijo(1,ind) +&
-          kpoint(2,ik) * xijo(2,ind) +&
-          kpoint(3,ik) * xijo(3,ind)
-          ckxij = cos(kxij)
-          skxij = -sin(kxij)
-        else 
-          ckxij=1.0_dp
-          skxij=0.0_dp
-        endif
-        cvar1 =  cmplx(S(ind)*ckxij,S(ind)*skxij,dp)
-        call m_set_element(Sauxms, jo, io, cvar1, complx_1, m_operation)
-      enddo
-    enddo
-    !
-    if(istpmove.eq.1) then   ! istpmove 
-      ! If first step calculate S0^1/2 and save for next step. 
-      call calculatesqrtS(Sauxms,invsqS,sqrtS(ik),nuo,m_storage,m_operation)
-    elseif(istpmove.gt.1) then 
-      ! Calculate both Sn^1/2 and Sn^-1/2 where Sn^1/2 is used in n+1 step. 
-      call calculatesqrtS(Sauxms,invsqS,Maux,nuo,m_storage,m_operation)
-      !Saux= Sn-1^1/2*Sn^-1/2
-      call mm_multiply(invsqS,'n',sqrtS(ik),'n',&
-      Sauxms,cmplx(1.0_dp,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
-      ! Passing Sn^1/2 from Maux to sqrtS for next step usage.
-      call m_add ( Maux,'n',sqrtS(ik),cmplx(1.0_dp,0.0_dp,dp),              &
-                  cmplx(0.0_dp,0.0_dp,dp),m_operation )
-      ! C1=S0^1/2*S1^1/2*C0
-      qe=2.0d0*kweight(ik)/dble(nspin)
-      do ispin=1,nspin
-        ! Cn = Saux*Cn-1 where Saux= Sn-1^1/2*Sn^-1/2
-        call m_allocate ( phi,wavef_ms(ik,ispin)%dim1,                      &
-                          wavef_ms(ik,ispin)%dim2,m_storage )
-        call m_add( wavef_ms(ik,ispin),'n',phi,cmplx(1.0,0.0,dp),           &
-                   cmplx(0.0_dp,0.0_dp,dp),m_operation )
-        call mm_multiply(Sauxms,'n',phi,'n',               &
-             wavef_ms(ik,ispin),cmplx(1.0_dp,0.0_dp,dp),     &
-             cmplx(0.0_dp,0.0_dp,dp),m_operation)
-        call m_deallocate(phi)
-      enddo  
-    endif   !istpmove 
-  call m_deallocate(Sauxms)
-  enddo          ! ik 
-  !
-  IF(istpmove.gt.1) THEN   ! istpmove 
-    IF (IONode) THEN
-      WRITE(*,'(a)') 'chgbasis: Computing DM in new basis'
-    END IF
-      call compute_tddm (Dscf)
-  END IF
-  call m_deallocate(Maux)
-  call m_deallocate(invsqS)
-
-  call timer('sankey_change_basis',2)
-#ifdef DEBUG
-  call write_debug( '    POS sankey_change_basis' )
-#endif
-  END SUBROUTINE sankey_change_basis 
-
-  SUBROUTINE calculatesqrtS(S,invsqS,sqrtS,nu,m_storage,m_operation)
-  
-    use precision 
-    use matdiagon
-    use MatrixSwitch
-    use parallelsubs,          only: LocalToGlobalOrb
-    use parallel,              only: Node, Nodes
-    use wavefunctions,         only: complx_0
-    ! 
-    implicit none
-    ! 
-    character(5), intent(in)                      :: m_storage
-    character(3), intent(in)                      :: m_operation
-    type(matrix), intent(inout)                   :: S,invsqS,sqrtS
-    type(matrix)                                  :: SD01, SD02
-    complex(dp)                                   :: varaux
-    real(dp)                                      :: eig01, eig02
-    real(dp), allocatable                         :: eigen(:)
-    integer                                       :: no,nu, info, i, j,jo
-    real(dp)  tiny
-    data tiny  /1.0d-10/
-    ! 
-    no=S%dim1 
-    allocate(eigen(no))
-    call m_allocate(SD01,no,no,m_storage)
-    call m_allocate(SD02,no,no,m_storage)
-    ! Takes overlap matrix S in dense form and returns its eigenvalues
-    ! in eigen(*) and eigenvectors in S.
-    call geteigen(S,eigen,m_operation)
-    !
-    do j=1,nu
-      call LocalToGlobalOrb(j,Node,Nodes,jo)
-      eig01=dsqrt(dabs(eigen(jo)))
-      eig02=1.0d0/(eig01+tiny)
-      do i=1,no
-        varaux = S%zval(i,j)
-        call m_set_element(SD01,i,jo,eig01*varaux,complx_0,m_operation)
-        call m_set_element(SD02,i,jo,eig02*varaux,complx_0,m_operation)
-      enddo
-    enddo 
-    !
-    deallocate(eigen)
-    
-    call mm_multiply(SD01,'n',S,'c',sqrtS,cmplx(1.0_dp,0.0_dp,dp),&
-                     cmplx(0.0_dp,0.0_dp,dp),m_operation)
-    call mm_multiply(SD02,'n',S,'c',invsqS,cmplx(1.0_dp,0.0_dp,dp),&
-                     cmplx(0.0_dp,0.0_dp,dp),m_operation)
-    call m_deallocate(SD01)
-    call m_deallocate(SD02)
-    !
-  END SUBROUTINE calculatesqrtS
-END MODULE m_sankey_change_basis
+!
+! Copyright (C) 1996-2016 The SIESTA group
+!  This file is distributed under the terms of the
+!  GNU General Public License: see COPYING in the top directory
+!  or http://www.gnu.org/copyleft/gpl.txt.
+! See Docs/Contributors.txt for a list of contributors.
+!
+MODULE m_sankey_change_basis 
+
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: sankey_change_basis
+
+  CONTAINS
+! *******************************************************************
+
+     SUBROUTINE sankey_change_basis ( istpmove )
+
+!******************************************************************************
+! This subroutine transforms the TDKS wavefunctions into new-basis during
+! Ehrenfest dynamics within TDDFT after atomic positions are changed. The 
+! TDKS are transformed according to the transformation proposed by Tomfohr and 
+! Sankey in Ref. phys. stat. sol. (b) 226 No.1, 115-123 (2001). 
+! 
+! After transforming TDKS states it calculates the density matrix in new basis set. 
+! We use MatrixSwitch to manipulate matrices.
+!
+! This subroutine is roughly based on D. Sanchez-Portal's chgbasis subroutine from, 
+! now absolete, serial version of tddft-siesta.
+!
+! Re-written for parallelization by Adiran Garaizar and Rafi Ullah, October 2015
+! Modified by Rafi Ullah, February 2017
+!********************************************************************************
+
+  use precision
+  use parallel,            only : Node, Nodes,BlockSize, IONode 
+  use parallelsubs,        only : GlobalToLocalOrb, GetNodeOrbs,               &
+                                  LocalToGlobalOrb
+  use m_diag_option,       only : ParallelOverK
+  use fdf
+  use alloc
+  use sys, only: die
+#ifdef MPI
+  use mpi_siesta,          only : mpi_bcast, mpi_comm_world, mpi_logical
+#endif
+  use m_spin,              only: nspin
+  use m_gamma,             only: gamma
+  use kpoint_scf_m,        only: kpoints_scf
+  use atomlist,            only: no_u, indxuo
+  use wavefunctions
+  use sparse_matrices,     only : numh, listhptr, listh, S, xijo, Dscf
+  use MatrixSwitch
+  use matdiagon,           only: geteigen 
+  !
+  implicit none
+  !
+  integer, intent(in)     :: istpmove
+  !
+#ifdef MPI
+  integer                 :: MPIerror,desch(9)
+  external                :: diagkp
+#endif
+  !
+  logical, save           :: frstme = .true.
+  integer                 :: io, iuo, juo, nuo, jo, ind, ispin
+  integer                 :: ik, j,ierror
+  real(dp)                :: skxij,ckxij, kxij, qe 
+  complex(dp)             :: cvar1
+
+  type(matrix)                     :: Maux,invsqS,phi
+  type(matrix)                     :: Sauxms
+  type(matrix),allocatable,save    :: sqrtS(:)
+  character(3)                     :: m_operation
+  character(5)                     :: m_storage
+  !
+#ifdef DEBUG
+  call write_debug( '    PRE sankey_change_basis' )
+#endif
+  !
+#ifdef MPI
+  call GetNodeOrbs(no_u,Node,Nodes,nuo)
+  if (frstme) then
+    if(ParallelOverK) then
+      call die ( "chgbasis: tddft-siesta not parallelized over k-points." )
+    endif
+      call ms_scalapack_setup(mpi_comm_world,1,'c',BlockSize)
+  endif
+#else
+  Node = 0
+  Nodes = 1
+  nuo = no_u
+#endif
+  call timer( 'sankey_change_basis', 1 )
+  !
+#ifdef MPI
+  m_storage='pzdbc'
+  m_operation='lap'
+#else
+  m_storage='szden'
+  m_operation='lap'
+#endif
+  !
+  IF (nspin .eq. 4) THEN
+    call die ('chgbasis: ERROR: EID not yet prepared for non-collinear spin')
+  END IF
+    ! Allocate local arrays
+  if(frstme) then
+    allocate(sqrtS(kpoints_scf%N))
+    do ik=1,kpoints_scf%N
+      call m_allocate(sqrtS(ik),no_u,no_u,m_storage)
+    end do
+    frstme=.false.
+  endif
+  call m_allocate(Maux,no_u,no_u,m_storage)
+  call m_allocate(invsqS,no_u,no_u,m_storage)
+  !
+  do ik = 1,kpoints_scf%N
+  call m_allocate(Sauxms,no_u,no_u,m_storage)
+    do iuo = 1,nuo
+      call LocalToGlobalOrb(iuo, Node, Nodes, io)
+      do j = 1,numh(iuo)
+        ind = listhptr(iuo) + j
+        juo = listh(ind)
+        jo  = indxuo (juo)
+        if(.not.gamma) then 
+          kxij = kpoints_scf%k(1,ik) * xijo(1,ind) +&
+          kpoints_scf%k(2,ik) * xijo(2,ind) +&
+          kpoints_scf%k(3,ik) * xijo(3,ind)
+          ckxij = cos(kxij)
+          skxij = -sin(kxij)
+        else 
+          ckxij=1.0_dp
+          skxij=0.0_dp
+        endif
+        cvar1 =  cmplx(S(ind)*ckxij,S(ind)*skxij,dp)
+        call m_set_element(Sauxms, jo, io, cvar1, complx_1, m_operation)
+      enddo
+    enddo
+    !
+    if(istpmove.eq.1) then   ! istpmove 
+      ! If first step calculate S0^1/2 and save for next step. 
+      call calculatesqrtS(Sauxms,invsqS,sqrtS(ik),nuo,m_storage,m_operation)
+    elseif(istpmove.gt.1) then 
+      ! Calculate both Sn^1/2 and Sn^-1/2 where Sn^1/2 is used in n+1 step. 
+      call calculatesqrtS(Sauxms,invsqS,Maux,nuo,m_storage,m_operation)
+      !Saux= Sn-1^1/2*Sn^-1/2
+      call mm_multiply(invsqS,'n',sqrtS(ik),'n',&
+      Sauxms,cmplx(1.0_dp,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
+      ! Passing Sn^1/2 from Maux to sqrtS for next step usage.
+      call m_add ( Maux,'n',sqrtS(ik),cmplx(1.0_dp,0.0_dp,dp),              &
+                  cmplx(0.0_dp,0.0_dp,dp),m_operation )
+      ! C1=S0^1/2*S1^1/2*C0
+      qe=2.0d0*kpoints_scf%w(ik)/dble(nspin)
+      do ispin=1,nspin
+        ! Cn = Saux*Cn-1 where Saux= Sn-1^1/2*Sn^-1/2
+        call m_allocate ( phi,wavef_ms(ik,ispin)%dim1,                      &
+                          wavef_ms(ik,ispin)%dim2,m_storage )
+        call m_add( wavef_ms(ik,ispin),'n',phi,cmplx(1.0,0.0,dp),           &
+                   cmplx(0.0_dp,0.0_dp,dp),m_operation )
+        call mm_multiply(Sauxms,'n',phi,'n',               &
+             wavef_ms(ik,ispin),cmplx(1.0_dp,0.0_dp,dp),     &
+             cmplx(0.0_dp,0.0_dp,dp),m_operation)
+        call m_deallocate(phi)
+      enddo  
+    endif   !istpmove 
+  call m_deallocate(Sauxms)
+  enddo          ! ik 
+  !
+  IF(istpmove.gt.1) THEN   ! istpmove 
+    IF (IONode) THEN
+      WRITE(*,'(a)') 'chgbasis: Computing DM in new basis'
+    END IF
+      call compute_tddm (Dscf)
+  END IF
+  call m_deallocate(Maux)
+  call m_deallocate(invsqS)
+
+  call timer('sankey_change_basis',2)
+#ifdef DEBUG
+  call write_debug( '    POS sankey_change_basis' )
+#endif
+  END SUBROUTINE sankey_change_basis 
+
+  SUBROUTINE calculatesqrtS(S,invsqS,sqrtS,nu,m_storage,m_operation)
+  
+    use precision 
+    use matdiagon
+    use MatrixSwitch
+    use parallelsubs,          only: LocalToGlobalOrb
+    use parallel,              only: Node, Nodes
+    use wavefunctions,         only: complx_0
+    ! 
+    implicit none
+    ! 
+    character(5), intent(in)                      :: m_storage
+    character(3), intent(in)                      :: m_operation
+    type(matrix), intent(inout)                   :: S,invsqS,sqrtS
+    type(matrix)                                  :: SD01, SD02
+    complex(dp)                                   :: varaux
+    real(dp)                                      :: eig01, eig02
+    real(dp), allocatable                         :: eigen(:)
+    integer                                       :: no,nu, info, i, j,jo
+    real(dp)  tiny
+    data tiny  /1.0d-10/
+    ! 
+    no=S%dim1 
+    allocate(eigen(no))
+    call m_allocate(SD01,no,no,m_storage)
+    call m_allocate(SD02,no,no,m_storage)
+    ! Takes overlap matrix S in dense form and returns its eigenvalues
+    ! in eigen(*) and eigenvectors in S.
+    call geteigen(S,eigen,m_operation)
+    !
+    do j=1,nu
+      call LocalToGlobalOrb(j,Node,Nodes,jo)
+      eig01=dsqrt(dabs(eigen(jo)))
+      eig02=1.0d0/(eig01+tiny)
+      do i=1,no
+        varaux = S%zval(i,j)
+        call m_set_element(SD01,i,jo,eig01*varaux,complx_0,m_operation)
+        call m_set_element(SD02,i,jo,eig02*varaux,complx_0,m_operation)
+      enddo
+    enddo 
+    !
+    deallocate(eigen)
+    
+    call mm_multiply(SD01,'n',S,'c',sqrtS,cmplx(1.0_dp,0.0_dp,dp),&
+                     cmplx(0.0_dp,0.0_dp,dp),m_operation)
+    call mm_multiply(SD02,'n',S,'c',invsqS,cmplx(1.0_dp,0.0_dp,dp),&
+                     cmplx(0.0_dp,0.0_dp,dp),m_operation)
+    call m_deallocate(SD01)
+    call m_deallocate(SD02)
+    !
+  END SUBROUTINE calculatesqrtS
+END MODULE m_sankey_change_basis
 
=== modified file 'Src/siesta_analysis.F'
--- Src/siesta_analysis.F	2018-04-15 21:36:51 +0000
+++ Src/siesta_analysis.F	2018-05-27 11:17:14 +0000
@@ -40,7 +40,7 @@
       use files,        only : slabel
       use files,        only : filesOut_t   ! derived type for output file names
       use zmatrix,      only: lUseZmatrix, write_zmatrix
-      use Kpoint_grid
+      use kpoint_scf_m, only: kpoints_scf
       use parallel, only: IOnode
       use parallel, only: SIESTA_worker
       use files,       only : label_length
@@ -286,16 +286,17 @@
         ! The user is responsible for using appropriate values.
         wfs_band_min = fdf_get("WFS.BandMin",1)
         wfs_band_max = fdf_get("WFS.BandMax",no_u)
-        call setup_wfs_list(nkpnt,no_u,wfs_band_min,wfs_band_max,
+        call setup_wfs_list(kpoints_scf%N,no_u,
+     &                      wfs_band_min,wfs_band_max,
      $                      use_scf_weights=.true.,
      $                      use_energy_window=.true.)
          wfs_filename = trim(slabel)//".fullBZ.WFSX"
          if (IONode) print "(a)", "Writing WFSX for COOP/COHP in " 
      $                           // trim(wfs_filename)
          call wwave( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, 
-     .        nkpnt,
+     .        kpoints_scf%N,
      .        numh, listhptr, listh, H, S, Ef, xijo, indxuo,
-     .        nkpnt, kpoint, no_u, gamma, occtol)
+     .        kpoints_scf%N, kpoints_scf%k, no_u, gamma, occtol)
       endif
 
 !     Compute bands
@@ -353,7 +354,7 @@
           if ( h_spin_dim <= 2 ) then
             write(6,'(/,a,/,a4,a3,a7)')
      &       'siesta: Eigenvalues (eV):', 'ik', 'is', 'eps'
-            do ik = 1,nkpnt
+            do ik = 1,kpoints_scf%N
               do ispin = 1,spinor_dim
                 write(6,'(i4,i3,10f7.2)')
      &            ik,ispin,(eo(io,ispin,ik)/eV,io=1,min(10,neigwanted))
@@ -363,7 +364,7 @@
             enddo
           else
             write(6,'(/,a)') 'siesta: Eigenvalues (eV):'
-            do ik = 1,nkpnt
+            do ik = 1,kpoints_scf%N
               write(6,'(a,i6)') 'ik =', ik
               write(6,'(10f7.2)')
      &          ((eo(io,ispin,ik)/eV,io=1,neigwanted),ispin=1,2)
@@ -376,8 +377,9 @@
       if (((isolve.eq.SOLVE_DIAGON).or.
      &     ((isolve.eq.SOLVE_MINIM).and.minim_calc_eigenvalues))
      &     .and.IOnode)
-     &     call ioeig(eo,ef,neigwanted,nspin,nkpnt,no_u,spinor_dim,
-     &     nkpnt,kpoint, kweight)
+     &    call ioeig(eo,ef,neigwanted,nspin,kpoints_scf%N,
+     &     no_u,spinor_dim,
+     &     kpoints_scf%N, kpoints_scf%k, kpoints_scf%w)
 
 
 !**   This uses H,S, and xa, as it diagonalizes them again
 
=== modified file 'Src/siesta_dicts.F90'
--- Src/siesta_dicts.F90	2018-01-25 20:01:19 +0000
+++ Src/siesta_dicts.F90	2018-05-27 11:17:14 +0000
@@ -218,7 +218,7 @@
   subroutine dict_populate_variables()
 
     use siesta_geom
-    use kpoint_grid, only: kscell, kdispl
+    use kpoint_scf_m, only: kpoints_scf
     use m_forces
     use m_energies
     use atomlist
@@ -327,9 +327,9 @@
 
     ! Add the k-point sampling
     variables = variables // &
-         ('BZ.k.Matrix'.kvp.kscell)
+         ('BZ.k.Matrix'.kvp.kpoints_scf%k_cell)
     variables = variables // &
-         ('BZ.k.Displacement'.kvp.kdispl)
+         ('BZ.k.Displacement'.kvp.kpoints_scf%k_displ)
 
   end subroutine dict_populate_variables
 
 
=== modified file 'Src/siesta_forces.F90'
--- Src/siesta_forces.F90	2018-04-24 10:26:39 +0000
+++ Src/siesta_forces.F90	2018-05-27 11:17:14 +0000
@@ -99,7 +99,7 @@
     use m_ts_charge,           only: TS_RHOCORR_FERMI
     use m_ts_charge,           only: TS_RHOCORR_FERMI_TOLERANCE
     use m_transiesta,          only: transiesta
-    use kpoint_grid, only : gamma_scf
+    use kpoint_scf_m, only : gamma_scf
     use m_energies, only : Ef
 
     use m_initwf, only: initwf
 
=== modified file 'Src/siesta_init.F'
--- Src/siesta_init.F	2018-05-15 13:01:22 +0000
+++ Src/siesta_init.F	2018-05-27 11:17:14 +0000
@@ -12,8 +12,9 @@
       CONTAINS
 
       subroutine siesta_init()
-      use Kpoint_grid,        only: setup_Kpoint_grid, gamma_scf, nkpnt
-      USE Kpoint_pdos,        only: gamma_pdos
+      use kpoint_scf_m,       only: setup_kpoint_scf, kpoints_scf
+      use kpoint_scf_m,       only: gamma_scf
+      use kpoint_pdos_m,      only: gamma_pdos
       use Band,               only: gamma_bands, setup_bands
       use m_ksvinit,          only: gamma_polarization,
      &                              estimate_pol_kpoints
@@ -568,9 +569,16 @@
 !     NOTE: We need to know whether gamma is .true. or
 !     not early, in order to decide whether to use an 
 !     auxiliary supercell for the calculation of matrix elements.
-      call setup_Kpoint_grid( ucell )
+      call setup_kpoint_scf( ucell )
       gamma = gamma_scf
 
+!     Call initialisation of PDOS here since we need to check if 
+!     the auxiliary supercell is needed for a non-gamma calculation
+      call init_projected_DOS( )
+      if ( do_pdos ) then
+        gamma = gamma .and. gamma_pdos
+      endif
+
       ! Read in diagonalization routines
       ! Note that only the sampled BZ is responsible for
       ! ParallelOverK option, the remaning options are
@@ -578,17 +586,10 @@
       call read_diag(Gamma_scf, spin%H)
       call print_diag()
 
-!     Call initialisation of PDOS here since we need to check if 
-!     the auxiliary supercell is needed for a non-gamma calculation
-      call init_projected_DOS( )
-      if (do_pdos) then
-        gamma = gamma .and. gamma_pdos
-      endif
-
       nullify(eo,qo)
-      call re_alloc(eo, 1, no_u, 1, spin%spinor, 1, nkpnt, 'eo', 
+      call re_alloc(eo, 1, no_u, 1, spin%spinor, 1, kpoints_scf%N, 'eo', 
      &              'siesta_init')
-      call re_alloc(qo, 1, no_u, 1, spin%spinor, 1, nkpnt, 'qo', 
+      call re_alloc(qo, 1, no_u, 1, spin%spinor, 1, kpoints_scf%N, 'qo', 
      &              'siesta_init')
 
       call setup_bands( )
 
=== modified file 'Src/siesta_tddft.F90'
--- Src/siesta_tddft.F90	2017-08-29 14:04:45 +0000
+++ Src/siesta_tddft.F90	2018-05-27 11:17:14 +0000
@@ -35,7 +35,7 @@
     use atomlist,      only: no_s, no_l, no_u, indxuo
     use m_spin,        only: nspin
     use m_gamma
-    use Kpoint_grid,    only: nkpnt, kpoint, kweight
+    use kpoint_scf_m,    only: kpoints_scf
 
     use m_compute_energies, only: compute_energies
     use m_state_analysis, only: state_analysis
@@ -114,8 +114,8 @@
     ! ms_scalapack_setup. Which is now implicit in changebasis. 
 
     if ( istep == 1 ) then
-       allocate(wavef_ms(nkpnt,nspin))
-       call iowavef('read',wavef_ms,no_u,nkpnt,nspin)
+       allocate(wavef_ms(kpoints_scf%N,nspin))
+       call iowavef('read',wavef_ms,no_u,kpoints_scf%N,nspin)
        IF (IONode) THEN
        write(6,'(a)') 'Computing DM from initial KS wavefunctions'
        END IF
@@ -127,7 +127,7 @@
       ! The wavefunctions are saved after transforming into the current basis
       ! but before evolving them to future.This keeps the wavefunctions
       ! concurrent with atomic position.
-      if(fincoor .gt. 1) call iowavef('write',wavef_ms,no_u,nkpnt,nspin)
+      if(fincoor .gt. 1) call iowavef('write',wavef_ms,no_u,kpoints_scf%N,nspin)
     end if
 
     do itded = 1 , ntded ! TDED loop
@@ -142,7 +142,7 @@
 
        if (tdsavewf) then
          if (fincoor .eq. 1 .and. itded .eq. ntded) then
-           call iowavef('write',wavef_ms,no_u,nkpnt,nspin)
+           call iowavef('write',wavef_ms,no_u,kpoints_scf%N,nspin)
          endif
        end if
 
@@ -155,7 +155,7 @@
        
        call compute_energies (itded)
        call write_tddft(totime, istep, itded, ntded, rstart_time, &
-            etot, eo, no_u,nspin,nkpnt)
+            etot, eo, no_u,nspin,kpoints_scf%N)
        
     end do ! TDED loop
     call compute_tdEdm (Escf)
 
=== modified file 'Src/state_init.F'
--- Src/state_init.F	2018-05-15 13:01:22 +0000
+++ Src/state_init.F	2018-05-27 11:17:14 +0000
@@ -13,7 +13,9 @@
       CONTAINS
 
       subroutine state_init( istep )
-      use Kpoint_grid,       only: setup_Kpoint_grid, nkpnt
+      use kpoint_scf_m,      only: setup_kpoint_scf, kpoints_scf
+      use kpoint_t_m, only: kpoint_delete, kpoint_nullify
+
       use m_os,              only: file_exist
       use m_new_dm,          only: new_dm
       use m_proximity_check, only: proximity_check
@@ -78,7 +80,7 @@
       use sys,               only: message, die
       use m_sparse, only : xij_offset
 
-      use m_ts_kpoints, only: setup_ts_kpoint_grid
+      use ts_kpoint_scf_m, only: setup_ts_kpoint_scf, ts_kpoints_scf
       use m_ts_charge, only : TS_RHOCORR_METHOD, TS_RHOCORR_FERMI
       use m_ts_options, only : BTD_method
       use m_ts_options, only : TS_Analyze
@@ -263,14 +265,20 @@
      &     (istep.ne.inicoor) .and. (.not.gamma) ) then
 
 !       Will print k-points also
-        call setup_Kpoint_grid( ucell )
-
-        call setup_ts_kpoint_grid( ucell )
-
-        call re_alloc( eo, 1, no_u, 1, spin%spinor, 1, nkpnt, 'eo',
-     &                 'state_init')
-        call re_alloc( qo, 1, no_u, 1, spin%spinor, 1, nkpnt, 'qo',
-     &                 'state_init' )
+        call kpoint_delete(kpoints_scf)
+        call setup_kpoint_scf( ucell )
+
+        if ( TSmode ) then
+          call kpoint_delete( ts_kpoints_scf )
+        else
+          call kpoint_nullify( ts_kpoints_scf )
+        end if
+        call setup_ts_kpoint_scf( ucell, kpoints_scf )
+
+        call re_alloc( eo, 1, no_u, 1, spin%spinor, 1, kpoints_scf%N,
+     &                 'eo', 'state_init')
+        call re_alloc( qo, 1, no_u, 1, spin%spinor, 1, kpoints_scf%N, 
+     &                 'qo', 'state_init' )
 
 !       Find required supercell
         if (gamma) then
 
=== modified file 'Src/ts_init.F90'
--- Src/ts_init.F90	2017-11-19 20:18:39 +0000
+++ Src/ts_init.F90	2018-05-27 11:17:14 +0000
@@ -36,9 +36,9 @@
     use m_ts_gf,        only : do_Green, do_Green_Fermi
     use m_ts_electrode, only : init_Electrode_HS
     
-    use m_ts_kpoints, only : setup_ts_kpoint_grid
-    use m_ts_kpoints, only : ts_nkpnt, ts_kpoint, ts_kweight
-    use m_ts_kpoints, only : ts_kscell, ts_kdispl, ts_gamma
+    use kpoint_scf_m, only : kpoints_scf
+    use ts_kpoint_scf_m, only : setup_ts_kpoint_scf
+    use ts_kpoint_scf_m, only : ts_kpoints_scf, ts_Gamma_scf
     use m_ts_cctype
     use m_ts_electype
     use m_ts_options ! Just everything (easier)
@@ -91,17 +91,17 @@
     call read_ts_elec( ucell, na_u, xa, lasto )
 
     ! Read in the k-points
-    call setup_ts_kpoint_grid( ucell )
+    call setup_ts_kpoint_scf( ucell, kpoints_scf )
 
     ! Read after electrode stuff
     call read_ts_after_Elec( ucell, nspin, na_u, xa, lasto, &
-         ts_kscell, ts_kdispl)
+        ts_kpoints_scf%k_cell, ts_kpoints_scf%k_displ)
 
     ! Print the options
     call print_ts_options( ucell )
 
     ! Print all warnings
-    call print_ts_warnings( ts_Gamma, ucell, na_u, xa, Nmove )
+    call print_ts_warnings( ts_Gamma_scf, ucell, na_u, xa, Nmove )
 
     ! If we actually have a transiesta run we need to process accordingly!
     if ( .not. TSmode ) return
@@ -188,13 +188,13 @@
           call init_Electrode_HS(Elecs(i))
 
           call do_Green(Elecs(i), &
-               ucell,ts_nkpnt,ts_kpoint,ts_kweight, &
+               ucell,ts_kpoints_scf%N,ts_kpoints_scf%k,ts_kpoints_scf%w, &
                Elecs_xa_Eps, .false. )
 
           if ( TS_RHOCORR_METHOD == TS_RHOCORR_FERMI ) then
              
              call do_Green_Fermi(mean_kT, Elecs(i), &
-                  ucell,ts_nkpnt,ts_kpoint,ts_kweight, &
+                  ucell,ts_kpoints_scf%N,ts_kpoints_scf%k,ts_kpoints_scf%w, &
                   Elecs_xa_Eps, .false. )
 
           end if
 
=== added file 'Src/ts_kpoint_scf.F90'
--- Src/ts_kpoint_scf.F90	1970-01-01 00:00:00 +0000
+++ Src/ts_kpoint_scf.F90	2018-05-27 11:17:14 +0000
@@ -0,0 +1,89 @@
+! ---
+! Copyright (C) 1996-2016	The SIESTA group
+!  This file is distributed under the terms of the
+!  GNU General Public License: see COPYING in the top directory
+!  or http://www.gnu.org/copyleft/gpl.txt .
+! See Docs/Contributors.txt for a list of contributors.
+! ---
+module ts_kpoint_scf_m
+  !
+  ! Contains data structures and routines to deal with the kpoint-grid
+  ! for the self-consistent calculation
+  ! Other uses (bands, optical, polarization) have their own structures.
+  !
+  use precision, only : dp
+  
+  ! The k-point-type
+  use kpoint_t_m
+
+  implicit none
+
+  public :: setup_ts_kpoint_scf
+  public :: ts_kpoints_scf
+  public :: ts_gamma_scf
+  
+  private
+
+  logical, save :: ts_gamma_scf
+  type(kpoint_t), save :: ts_kpoints_scf
+
+contains
+
+  subroutine setup_ts_kpoint_scf( ucell, kpoints_scf )
+    use parallel, only: Node
+    use siesta_options, only: writek
+    use m_spin, only: TrSym
+    use m_ts_global_vars, only : TSmode
+
+    real(dp), intent(in) :: ucell(3,3)
+    type(kpoint_t), intent(in) :: kpoints_scf
+
+    call kpoint_read(ts_kpoints_scf, 'TS', ucell, TrSym, process_k_cell=process_k_cell_displ)
+
+    if ( ts_kpoints_scf%method == K_METHOD_NONE ) then
+
+      call kpoint_delete(ts_kpoints_scf)
+
+      if ( TSmode ) then
+        ! The user hasn't specified anything.
+        ! This means that we will use the default setting from siesta
+        call kpoint_read(ts_kpoints_scf, '', ucell, TrSym, process_k_cell=process_k_cell_displ)
+      else
+        ! To limit memory usage for very high number of k-points
+        call kpoint_associate(ts_kpoints_scf, kpoints_scf)
+      end if
+
+    end if
+
+    ts_gamma_scf = (ts_kpoints_scf%N == 1 .and. &
+        dot_product(ts_kpoints_scf%k(:,1),ts_kpoints_scf%k(:,1)) < 1.0e-20_dp)
+
+    ! Quick-return if non-IO or not a transiesta run
+    if ( .not. TSmode ) return
+    if ( Node /= 0 ) return
+
+    call kpoint_write_stdout(ts_kpoints_scf, writek, 'transiesta')
+    call kpoint_write_xml(ts_kpoints_scf, 'TS')
+    call kpoint_write_file(ts_kpoints_scf, 'TS.KP')
+
+  end subroutine setup_ts_kpoint_scf
+
+  subroutine process_k_cell_displ(k_cell, k_displ)
+    use m_ts_global_vars, only : TSmode
+    use m_ts_tdir, only: ts_tidx
+
+    integer, intent(inout) :: k_cell(3,3)
+    real(dp), intent(inout) :: k_displ(3)
+    integer :: i
+
+    if ( TSmode .and. ts_tidx > 0 ) then
+      i = ts_tidx
+      k_cell(:,i) = 0
+      k_cell(i,:) = 0
+      k_cell(i,i) = 1
+      k_displ(i) = 0._dp
+    end if
+
+  end subroutine process_k_cell_displ
+  
+end module ts_kpoint_scf_m
 
=== modified file 'Src/wavefunctions.F90'
--- Src/wavefunctions.F90	2017-06-29 12:15:14 +0000
+++ Src/wavefunctions.F90	2018-05-27 11:17:14 +0000
@@ -1,333 +1,333 @@
-MODULE  wavefunctions
-  USE  precision
-  USE  MatrixSwitch
-  !
-  implicit none 
-  !    
-  PRIVATE
-  PUBLIC       :: iowavef, compute_tddm, compute_tdEdm
-  type(matrix), allocatable, save, public :: wavef_ms(:,:)
-  TYPE(matrix), ALLOCATABLE, SAVE, public :: Hsave(:,:)
-  ! 
-  complex(kind=dp), parameter, public :: complx_1 = cmplx(1.0,0.0,dp)
-  complex(kind=dp), parameter, public :: complx_0 = cmplx(0.0,0.0,dp)
-
-CONTAINS
-  !
-  subroutine iowavef(task,wavef_rw,nuotot,nk,nspin)
- 
-    ! To read and write the TDKS orbitals
-    ! 
-    ! Written by Daniel Sanchez-Portal
-    ! Re-written by Rafi Ullah, CIC nanoGUNE, October 16, 2015 to work
-    ! with parallel TDDFT flavor of siesta using MatrixSwtich.
-    !
-    ! Although it is prepared to work with parallel TDDFT-siesta the
-    ! reading/writing process itself is not parallel.
-    use precision
-    use parallel
-    use fdf
-    use MatrixSwitch
-#ifdef MPI
-    use mpi_siesta
-#endif
-    !
-    implicit none
-    !
-    integer, intent(in)             :: nuotot, nk, nspin
-    character*(*), intent(in)       :: task
-    type(matrix), intent(inout)     :: wavef_rw(nk,nspin)
-    ! Internal variables and arrays
-    character         :: fname*33, sname*30, m_storage*5
-    logical           ::  exist1, frstme
-    integer           :: unit1, ie,nuototread,nkread,nspinread,dim2read
-    integer           :: nwf, ik, ispin, mxnwf, io, ix,i,j
-    external          :: chkdim, io_assign, io_close, timer,memory
-    complex(kind=dp)  :: varaux
-    save              :: frstme, fname
-    data                 frstme /.true./
-    !
-#ifdef MPI
-    INTEGER :: MPIerror
-#endif
-    !
-#ifdef MPI
-    call MPI_Comm_Rank(MPI_Comm_World,Node,MPIerror)
-    call MPI_Comm_Size(MPI_Comm_World,Nodes,MPIerror)
-#else
-    Node = 0
-    Nodes = 1
-#endif
-    !
-#ifdef MPI
-      m_storage='pzdbc'
-#else
-      m_storage='szden'
-#endif
-    ! Find file name
-    if (frstme) then
-      if (Node.eq.0) then
-        sname = fdf_string('SystemLabel','siesta')
-      endif
-#ifdef MPI
-      call MPI_Bcast(sname,30,MPI_character,0,MPI_Comm_World,MPIerror)
-#endif
-      fname = trim(sname)//'.TDWF'
-      frstme = .false.
-    endif
-    !
-    if (task.eq.'read' .or. task.eq.'READ') then      ! task read
-      if (Node.eq.0) then
-        inquire (file=fname, exist=exist1)
-        if(exist1) then
-          write(6,'(/,a)')                                            &
-          'iowavef: reading initial KS wavefunctions from file'
-         else
-          write(6,'(/,a)')                                            &
-          'iowavef: file containing the KS orbitals not found'
-          write(6,'(a)')                                              &
-          'iowavef: simulation can not be started or restarted'
-          stop 
-        endif
-        call io_assign(unit1)
-        open( unit1, file=fname, form='unformatted',                 &
-        status='old' )
-      endif
-      if (Node.eq.0) then
-        read(unit1) nuototread, nkread, nspinread
-        if(nkread.ne.nk) stop 'iowavef: Nunber of K-points inconsistent '
-        if(nuototread.ne.nuotot) stop 'iowavef: Number of KS orbitals inconsistent ' 
-        if(nspinread.ne.nspin) stop 'iowavef: Spin inconsistent '
-      endif
-      !   
-      do ispin=1,nspin                                 
-        do ik=1,nk
-          if (Node.eq.0) read(unit1) dim2read
-#ifdef MPI
-          call MPI_Bcast(dim2read,1,MPI_integer,0,MPI_Comm_World,MPIerror)
-#endif
-          call m_allocate(wavef_rw(ik,ispin),nuotot,dim2read,m_storage)
-        enddo
-      enddo
-      !
-      do ispin=1,nspin
-        do ik=1,nk
-          do i=1,nuotot
-            do j=1,wavef_rw(ik,ispin)%dim2
-              if (Node==0) read(unit1) varaux
-#ifdef MPI 
-              call MPI_Bcast(varaux,1,MPI_double_complex,0,MPI_Comm_World,MPIerror)
-#endif
-              call m_set_element(wavef_rw(ik,ispin),i,j,varaux,complx_0,'lap')
-            end do
-          end do
-        enddo
-      enddo
-      if (Node.eq.0)   call io_close(unit1)
-      !......................................................................   
-    elseif(task.eq.'write'.or.task.eq.'WRITE') then     ! task write
-      if (Node.eq.0) then
-        call io_assign(unit1)
-        open( unit1, file=fname, form='unformatted',status='replace' )
-        rewind(unit1)
-        write(unit1) nuotot,nk,nspin
-        !
-        do ispin=1,nspin
-          do ik=1,nk
-            write(unit1) (wavef_rw(ik,ispin)%dim2)
-          enddo
-        enddo
-        !
-      end if
-      !
-      do ispin=1,nspin
-        do ik=1,nk
-          do i=1,nuotot
-            do j=1,wavef_rw(ik,ispin)%dim2
-              call m_get_element(wavef_rw(ik,ispin),i,j,varaux,'lap')
-              if (Node==0) write(unit1) varaux
-            enddo
-          enddo
-        enddo 
-      enddo 
-      !
-      if (Node .eq. 0) then
-        call io_close(unit1)
-      endif
-    endif            ! end task 
-    !
-  end subroutine iowavef
-  !
-  subroutine  compute_tddm (Dnew)
-  
-    ! To calculate the time-dependent density matrix from time-dependent
-    ! wavefunctions. 
-    ! Written by Rafi Ullah January 2017
-    !***********************************
-    use sparse_matrices,      only: numh, maxnh, listh, listhptr,     &
-                                    xijo   
-    use Kpoint_grid,          only: nkpnt, kpoint, gamma_scf, kweight
-    use atomlist,             only: no_l, no_u, indxuo
-    use m_spin,               only: nspin
-    integer                      :: ispin, nuo, nuotot
-    integer                      :: io,jo, juo, j, ind, ik
-    real(dp)                     :: Dnew (maxnh, nspin), wk, kxij
-    real(dp)                     :: ckxij, skxij
-    complex(dp)                  :: varaux
-    type(matrix)                 :: Daux    
-    character                    :: m_storage*5, m_operation*3
-#ifdef MPI
-    m_storage   ='pzdbc'
-    m_operation = 'lap'
-#else
-    m_storage='szden'
-    m_operation = 'lap'
-#endif
-    Dnew(1:maxnh,1:nspin) = 0.0_dp
-    call m_allocate(Daux,no_u,no_u,m_storage)
-    Do ik = 1, nkpnt
-      Do ispin =1, nspin
-      wk = 2.00_dp*kweight(ik)/dble(nspin)
-      ! Calculating density matrix using MatrixSwitch.
-      call mm_multiply(wavef_ms(ik,ispin),'n',wavef_ms(ik,ispin),'c',Daux,           & 
-                    cmplx(wk,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
-      ! Passing DM from dense to sparse form.  
-      DO io=1,no_l
-        DO j=1,numh(io)
-          ind=listhptr(io) + j
-          juo = listh(ind)
-          jo  = indxuo(juo)
-          IF (.NOT. gamma_scf) THEN
-            kxij = kpoint(1,ik) * xijo(1,ind) + &
-            kpoint(2,ik) * xijo(2,ind) +        &
-            kpoint(3,ik) * xijo(3,ind) 
-            ckxij =  cos(kxij)
-            skxij = -sin(kxij)
-          ELSE
-            ckxij = 1.0_dp
-            skxij = 0.0_dp
-          END IF
-          varaux = real(Daux%zval(jo,io))*ckxij +     &
-                   aimag(Daux%zval(jo,io))*skxij
-          Dnew(ind,ispin)  = Dnew(ind,ispin) + varaux
-        END DO
-      END DO
-    END DO
-  END DO
-    call m_deallocate (Daux)
-  end subroutine compute_tddm
-!----------------------------------------------------------------------------------!
-
-
-!----------------------------------------------------------------------------------!
-      subroutine compute_tdEdm(Enew)
-!**********************************************************************************!
-! Re-written by Rafi Ullah in June 2017.                                           !
-! In this subroutine we calculate the energy density matrix in parallel using      !
-! MatrixSwitch.                                                                    !
-!**********************************************************************************!
-      use parallel
-      use precision
-      use parallelsubs,     only: LocalToGlobalOrb
-      use sparse_matrices,  only: numh, maxnh, listh, listhptr, xijo
-      use sparse_matrices,  only: H, S
-      use Kpoint_grid,      only: nkpnt, kpoint, kweight
-      use m_gamma,          only: gamma
-      use atomlist,         only: no_l, no_u, indxuo
-      use m_spin,           only: nspin
-      use MatrixSwitch
-      use matswinversion,   only: getinverse
-      !
-      implicit none
-      !
-      real(dp)               :: Enew(maxnh, nspin)
-      real(dp)               :: wk, kxij, ckxij, skxij
-      integer                :: ik, ispin, i, j, io, jo, ind, juo, nocc
-      logical, save          :: frstime = .true.
-      character              :: m_storage*5, m_operation*3
-      complex(dp)            :: cvar1, cvar2
-      !
-      type(matrix)           :: Sauxms, Hauxms,psi,Eaux
-      type(matrix)           :: S_1
-      !
-#ifdef MPI
-      m_storage='pzdbc'
-      m_operation='lap'
-#else
-      m_storage='szden'
-      m_operation='lap'
-#endif
-      Enew(1:maxnh,1:nspin) = 0.0_dp
-      call m_allocate(S_1,no_u,no_u,m_storage)
-      call m_allocate(Eaux,no_u,no_u,m_storage)
-      Do ik=1,nkpnt
-        Do ispin=1,nspin
-          wk = 2.00_dp*kweight(ik)/dble(nspin)
-          nocc = wavef_ms(ik,ispin)%dim2
-          call m_allocate(psi,nocc,no_u,m_storage)
-          call m_allocate (Hauxms,no_u, no_u, m_storage)
-          call m_allocate (Sauxms,no_u, no_u, m_storage)
-          Do i = 1, no_l
-            call LocalToGlobalOrb(i, Node, Nodes, io)
-            Do j = 1, numh(i)
-              ind = listhptr(i) + j
-              juo = listh(ind)
-              jo  = indxuo(juo)
-              if( .not. gamma ) then
-                kxij = kpoint(1,ik)*xijo(1,ind) + kpoint(2,ik)*xijo(2,ind) +         &
-                       kpoint(3,ik)*xijo(3,ind)
-                ckxij =  cos(kxij)
-                skxij = -sin(kxij)
-              else
-                ckxij = 1.0_dp
-                skxij = 0.0_dp
-              end if
-              cvar1 = cmplx(H(ind,ispin)*ckxij,H(ind,ispin)*skxij,dp)
-              cvar2 = cmplx(S(ind)*ckxij,S(ind)*skxij,dp)
-              call m_set_element(Hauxms,jo,io,cvar1,complx_1,m_operation)
-              call m_set_element(Sauxms,jo,io,cvar2,complx_1,m_operation)
-            end do
-          end do
-          !
-          if (ispin .eq. 1) then
-            ! copy Sauxms to S_1
-            call m_add (Sauxms,'n',S_1,complx_1,complx_0,m_operation)
-            ! Invert the overlap matrix.
-            call getinverse(S_1,m_operation)
-          end if
-          call mm_multiply(Hauxms,'n',S_1,'n',Eaux,complx_1,complx_0,m_operation)
-          call mm_multiply(wavef_ms(ik,ispin),'c',Eaux,'n',psi,complx_1,            &
-                           complx_0,m_operation)
-          call mm_multiply(wavef_ms(ik,ispin),'n',psi,'n',Eaux,cmplx(wk,0.0,dp),    &
-                           complx_0,m_operation)
-          ! Passing Energy-DM from dense to sparse form.  
-          DO io=1,no_l
-            DO j=1,numh(io)
-              ind=listhptr(io) + j
-              juo = listh(ind)
-              jo  = indxuo(juo)
-              IF (.NOT. gamma) THEN
-                kxij = kpoint(1,ik) * xijo(1,ind) + &
-                kpoint(2,ik) * xijo(2,ind) +        &
-                kpoint(3,ik) * xijo(3,ind) 
-                ckxij =  cos(kxij)
-                skxij = -sin(kxij)
-              ELSE
-                ckxij = 1.0_dp
-                skxij = 0.0_dp
-              END IF
-              cvar1 =  real(Eaux%zval(jo,io))*ckxij +     &
-                       aimag(Eaux%zval(jo,io))*skxij
-              Enew(ind,ispin)  = Enew(ind,ispin) + cvar1
-            END DO
-          END DO
-          call m_deallocate(Sauxms)
-          call m_deallocate(Hauxms)
-          call m_deallocate(psi)
-        end do
-      end do
-      call m_deallocate(S_1)
-      call m_deallocate(Eaux)
-      end subroutine compute_tdEdm 
-  !
-end module wavefunctions       
+MODULE  wavefunctions
+  USE  precision
+  USE  MatrixSwitch
+  !
+  implicit none 
+  !    
+  PRIVATE
+  PUBLIC       :: iowavef, compute_tddm, compute_tdEdm
+  type(matrix), allocatable, save, public :: wavef_ms(:,:)
+  TYPE(matrix), ALLOCATABLE, SAVE, public :: Hsave(:,:)
+  ! 
+  complex(kind=dp), parameter, public :: complx_1 = cmplx(1.0,0.0,dp)
+  complex(kind=dp), parameter, public :: complx_0 = cmplx(0.0,0.0,dp)
+
+CONTAINS
+  !
+  subroutine iowavef(task,wavef_rw,nuotot,nk,nspin)
+ 
+    ! To read and write the TDKS orbitals
+    ! 
+    ! Written by Daniel Sanchez-Portal
+    ! Re-written by Rafi Ullah, CIC nanoGUNE, October 16, 2015 to work
+    ! with parallel TDDFT flavor of siesta using MatrixSwtich.
+    !
+    ! Although it is prepared to work with parallel TDDFT-siesta the
+    ! reading/writing process itself is not parallel.
+    use precision
+    use parallel
+    use fdf
+    use MatrixSwitch
+#ifdef MPI
+    use mpi_siesta
+#endif
+    !
+    implicit none
+    !
+    integer, intent(in)             :: nuotot, nk, nspin
+    character*(*), intent(in)       :: task
+    type(matrix), intent(inout)     :: wavef_rw(nk,nspin)
+    ! Internal variables and arrays
+    character         :: fname*33, sname*30, m_storage*5
+    logical           ::  exist1, frstme
+    integer           :: unit1, ie,nuototread,nkread,nspinread,dim2read
+    integer           :: nwf, ik, ispin, mxnwf, io, ix,i,j
+    external          :: chkdim, io_assign, io_close, timer,memory
+    complex(kind=dp)  :: varaux
+    save              :: frstme, fname
+    data                 frstme /.true./
+    !
+#ifdef MPI
+    INTEGER :: MPIerror
+#endif
+    !
+#ifdef MPI
+    call MPI_Comm_Rank(MPI_Comm_World,Node,MPIerror)
+    call MPI_Comm_Size(MPI_Comm_World,Nodes,MPIerror)
+#else
+    Node = 0
+    Nodes = 1
+#endif
+    !
+#ifdef MPI
+      m_storage='pzdbc'
+#else
+      m_storage='szden'
+#endif
+    ! Find file name
+    if (frstme) then
+      if (Node.eq.0) then
+        sname = fdf_string('SystemLabel','siesta')
+      endif
+#ifdef MPI
+      call MPI_Bcast(sname,30,MPI_character,0,MPI_Comm_World,MPIerror)
+#endif
+      fname = trim(sname)//'.TDWF'
+      frstme = .false.
+    endif
+    !
+    if (task.eq.'read' .or. task.eq.'READ') then      ! task read
+      if (Node.eq.0) then
+        inquire (file=fname, exist=exist1)
+        if(exist1) then
+          write(6,'(/,a)')                                            &
+          'iowavef: reading initial KS wavefunctions from file'
+         else
+          write(6,'(/,a)')                                            &
+          'iowavef: file containing the KS orbitals not found'
+          write(6,'(a)')                                              &
+          'iowavef: simulation can not be started or restarted'
+          stop 
+        endif
+        call io_assign(unit1)
+        open( unit1, file=fname, form='unformatted',                 &
+        status='old' )
+      endif
+      if (Node.eq.0) then
+        read(unit1) nuototread, nkread, nspinread
+        if(nkread.ne.nk) stop 'iowavef: Nunber of K-points inconsistent '
+        if(nuototread.ne.nuotot) stop 'iowavef: Number of KS orbitals inconsistent ' 
+        if(nspinread.ne.nspin) stop 'iowavef: Spin inconsistent '
+      endif
+      !   
+      do ispin=1,nspin                                 
+        do ik=1,nk
+          if (Node.eq.0) read(unit1) dim2read
+#ifdef MPI
+          call MPI_Bcast(dim2read,1,MPI_integer,0,MPI_Comm_World,MPIerror)
+#endif
+          call m_allocate(wavef_rw(ik,ispin),nuotot,dim2read,m_storage)
+        enddo
+      enddo
+      !
+      do ispin=1,nspin
+        do ik=1,nk
+          do i=1,nuotot
+            do j=1,wavef_rw(ik,ispin)%dim2
+              if (Node==0) read(unit1) varaux
+#ifdef MPI 
+              call MPI_Bcast(varaux,1,MPI_double_complex,0,MPI_Comm_World,MPIerror)
+#endif
+              call m_set_element(wavef_rw(ik,ispin),i,j,varaux,complx_0,'lap')
+            end do
+          end do
+        enddo
+      enddo
+      if (Node.eq.0)   call io_close(unit1)
+      !......................................................................   
+    elseif(task.eq.'write'.or.task.eq.'WRITE') then     ! task write
+      if (Node.eq.0) then
+        call io_assign(unit1)
+        open( unit1, file=fname, form='unformatted',status='replace' )
+        rewind(unit1)
+        write(unit1) nuotot,nk,nspin
+        !
+        do ispin=1,nspin
+          do ik=1,nk
+            write(unit1) (wavef_rw(ik,ispin)%dim2)
+          enddo
+        enddo
+        !
+      end if
+      !
+      do ispin=1,nspin
+        do ik=1,nk
+          do i=1,nuotot
+            do j=1,wavef_rw(ik,ispin)%dim2
+              call m_get_element(wavef_rw(ik,ispin),i,j,varaux,'lap')
+              if (Node==0) write(unit1) varaux
+            enddo
+          enddo
+        enddo 
+      enddo 
+      !
+      if (Node .eq. 0) then
+        call io_close(unit1)
+      endif
+    endif            ! end task 
+    !
+  end subroutine iowavef
+  !
+  subroutine  compute_tddm (Dnew)
+  
+    ! To calculate the time-dependent density matrix from time-dependent
+    ! wavefunctions. 
+    ! Written by Rafi Ullah January 2017
+    !***********************************
+    use sparse_matrices,      only: numh, maxnh, listh, listhptr,     &
+                                    xijo   
+    use kpoint_scf_m,         only: kpoints_scf, gamma_scf
+    use atomlist,             only: no_l, no_u, indxuo
+    use m_spin,               only: nspin
+    integer                      :: ispin, nuo, nuotot
+    integer                      :: io,jo, juo, j, ind, ik
+    real(dp)                     :: Dnew (maxnh, nspin), wk, kxij
+    real(dp)                     :: ckxij, skxij
+    complex(dp)                  :: varaux
+    type(matrix)                 :: Daux    
+    character                    :: m_storage*5, m_operation*3
+#ifdef MPI
+    m_storage   ='pzdbc'
+    m_operation = 'lap'
+#else
+    m_storage='szden'
+    m_operation = 'lap'
+#endif
+    Dnew(1:maxnh,1:nspin) = 0.0_dp
+    call m_allocate(Daux,no_u,no_u,m_storage)
+    Do ik = 1, kpoints_scf%N
+      Do ispin =1, nspin
+      wk = 2.00_dp*kpoints_scf%w(ik)/dble(nspin)
+      ! Calculating density matrix using MatrixSwitch.
+      call mm_multiply(wavef_ms(ik,ispin),'n',wavef_ms(ik,ispin),'c',Daux,           & 
+                    cmplx(wk,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
+      ! Passing DM from dense to sparse form.  
+      DO io=1,no_l
+        DO j=1,numh(io)
+          ind=listhptr(io) + j
+          juo = listh(ind)
+          jo  = indxuo(juo)
+          IF (.NOT. gamma_scf) THEN
+            kxij = kpoints_scf%k(1,ik) * xijo(1,ind) + &
+            kpoints_scf%k(2,ik) * xijo(2,ind) +        &
+            kpoints_scf%k(3,ik) * xijo(3,ind) 
+            ckxij =  cos(kxij)
+            skxij = -sin(kxij)
+          ELSE
+            ckxij = 1.0_dp
+            skxij = 0.0_dp
+          END IF
+          varaux = real(Daux%zval(jo,io))*ckxij +     &
+                   aimag(Daux%zval(jo,io))*skxij
+          Dnew(ind,ispin)  = Dnew(ind,ispin) + varaux
+        END DO
+      END DO
+    END DO
+  END DO
+    call m_deallocate (Daux)
+  end subroutine compute_tddm
+!----------------------------------------------------------------------------------!
+
+
+!----------------------------------------------------------------------------------!
+      subroutine compute_tdEdm(Enew)
+!**********************************************************************************!
+! Re-written by Rafi Ullah in June 2017.                                           !
+! In this subroutine we calculate the energy density matrix in parallel using      !
+! MatrixSwitch.                                                                    !
+!**********************************************************************************!
+      use parallel
+      use precision
+      use parallelsubs,     only: LocalToGlobalOrb
+      use sparse_matrices,  only: numh, maxnh, listh, listhptr, xijo
+      use sparse_matrices,  only: H, S
+      use kpoint_scf_m,     only: kpoints_scf, gamma_scf
+      use m_gamma,          only: gamma
+      use atomlist,         only: no_l, no_u, indxuo
+      use m_spin,           only: nspin
+      use MatrixSwitch
+      use matswinversion,   only: getinverse
+      !
+      implicit none
+      !
+      real(dp)               :: Enew(maxnh, nspin)
+      real(dp)               :: wk, kxij, ckxij, skxij
+      integer                :: ik, ispin, i, j, io, jo, ind, juo, nocc
+      logical, save          :: frstime = .true.
+      character              :: m_storage*5, m_operation*3
+      complex(dp)            :: cvar1, cvar2
+      !
+      type(matrix)           :: Sauxms, Hauxms,psi,Eaux
+      type(matrix)           :: S_1
+      !
+#ifdef MPI
+      m_storage='pzdbc'
+      m_operation='lap'
+#else
+      m_storage='szden'
+      m_operation='lap'
+#endif
+      Enew(1:maxnh,1:nspin) = 0.0_dp
+      call m_allocate(S_1,no_u,no_u,m_storage)
+      call m_allocate(Eaux,no_u,no_u,m_storage)
+      Do ik=1,kpoints_scf%N
+        Do ispin=1,nspin
+          wk = 2.00_dp*kpoints_scf%w(ik)/dble(nspin)
+          nocc = wavef_ms(ik,ispin)%dim2
+          call m_allocate(psi,nocc,no_u,m_storage)
+          call m_allocate (Hauxms,no_u, no_u, m_storage)
+          call m_allocate (Sauxms,no_u, no_u, m_storage)
+          Do i = 1, no_l
+            call LocalToGlobalOrb(i, Node, Nodes, io)
+            Do j = 1, numh(i)
+              ind = listhptr(i) + j
+              juo = listh(ind)
+              jo  = indxuo(juo)
+              if( .not. gamma_scf ) then
+                kxij = kpoints_scf%k(1,ik)*xijo(1,ind) + kpoints_scf%k(2,ik)*xijo(2,ind) +         &
+                       kpoints_scf%k(3,ik)*xijo(3,ind)
+                ckxij =  cos(kxij)
+                skxij = -sin(kxij)
+              else
+                ckxij = 1.0_dp
+                skxij = 0.0_dp
+              end if
+              cvar1 = cmplx(H(ind,ispin)*ckxij,H(ind,ispin)*skxij,dp)
+              cvar2 = cmplx(S(ind)*ckxij,S(ind)*skxij,dp)
+              call m_set_element(Hauxms,jo,io,cvar1,complx_1,m_operation)
+              call m_set_element(Sauxms,jo,io,cvar2,complx_1,m_operation)
+            end do
+          end do
+          !
+          if (ispin .eq. 1) then
+            ! copy Sauxms to S_1
+            call m_add (Sauxms,'n',S_1,complx_1,complx_0,m_operation)
+            ! Invert the overlap matrix.
+            call getinverse(S_1,m_operation)
+          end if
+          call mm_multiply(Hauxms,'n',S_1,'n',Eaux,complx_1,complx_0,m_operation)
+          call mm_multiply(wavef_ms(ik,ispin),'c',Eaux,'n',psi,complx_1,            &
+                           complx_0,m_operation)
+          call mm_multiply(wavef_ms(ik,ispin),'n',psi,'n',Eaux,cmplx(wk,0.0,dp),    &
+                           complx_0,m_operation)
+          ! Passing Energy-DM from dense to sparse form.  
+          DO io=1,no_l
+            DO j=1,numh(io)
+              ind=listhptr(io) + j
+              juo = listh(ind)
+              jo  = indxuo(juo)
+              IF (.NOT. gamma_scf) THEN
+                kxij = kpoints_scf%k(1,ik) * xijo(1,ind) + &
+                kpoints_scf%k(2,ik) * xijo(2,ind) +        &
+                kpoints_scf%k(3,ik) * xijo(3,ind) 
+                ckxij =  cos(kxij)
+                skxij = -sin(kxij)
+              ELSE
+                ckxij = 1.0_dp
+                skxij = 0.0_dp
+              END IF
+              cvar1 =  real(Eaux%zval(jo,io))*ckxij +     &
+                       aimag(Eaux%zval(jo,io))*skxij
+              Enew(ind,ispin)  = Enew(ind,ispin) + cvar1
+            END DO
+          END DO
+          call m_deallocate(Sauxms)
+          call m_deallocate(Hauxms)
+          call m_deallocate(psi)
+        end do
+      end do
+      call m_deallocate(S_1)
+      call m_deallocate(Eaux)
+      end subroutine compute_tdEdm 
+  !
+end module wavefunctions       
 
=== modified file 'Src/writewave.F'
--- Src/writewave.F	2017-10-10 19:27:53 +0000
+++ Src/writewave.F	2018-05-27 11:17:14 +0000
@@ -693,7 +693,7 @@
       use parallel,     only : Node, Nodes
       use parallelsubs, only : GlobalToLocalOrb, WhichNodeOrb
       use units,        only : eV
-      use kpoint_grid,  only : kweight
+      use kpoint_scf_m, only: kpoints_scf
 
 #ifdef MPI
       use mpi_siesta
@@ -788,7 +788,7 @@
      $                   // "wavefunction coefficients from WFSX file"
         endif
         if (scf_set) then
-           kpoint_weight = kweight(ik)
+           kpoint_weight = kpoints_scf%w(ik)
         else
            kpoint_weight = 1.0_dp
         endif
 
=== added file 'Tests/Reference/si_coop_kp_file.out'
--- Tests/Reference/si_coop_kp_file.out	1970-01-01 00:00:00 +0000
+++ Tests/Reference/si_coop_kp_file.out	2018-05-27 11:17:14 +0000
@@ -0,0 +1,483 @@
+Siesta Version  : trunk-697
+Architecture    : x86_64-linux-gcc
+Compiler version: GNU Fortran (GCC) 7.3.0
+Compiler flags  : mpifort -m64 -fPIC -O3 -ftree-vectorize -fexpensive-optimizations -funroll-loops -fprefetch-loop-arrays -fgraphite -fipa-icf -fipa-pure-const -march=native -fschedule-fusion -fselective-scheduling -fstack-arrays -fmax-stack-var-size=10000 -fipa-sra -fipa-cp -fno-second-underscore
+PP flags        : -DSIESTA__FLOOK -DSIESTA__METIS -DSIESTA__FFTW -DSIESTA__MRRR -DMPI -DFC_HAVE_FLUSH -DFC_HAVE_ABORT -DCDF -DGRID_DP -DMPI_TIMING -DTRANSIESTA_TIMING -DTBTRANS_TIMING -DSIESTA__MUMPS -DNCDF_PARALLEL -DNCDF -DNCDF_4 -DSCALAPACK_DEBUG -D_DIAG_WORK -DSIESTA__ELPA
+Libraries       : libncdf.a libfdict.a  -lzmumps -lmumps_common -lpord -lparmetis -lmetis -lesmumps -lscotch -lscotcherr -L/opt/elpa/2017.05.003/gnu-7.3.0/lib -Wl,-rpath=/opt/elpa/2017.05.003/gnu-7.3.0/lib -lelpa -L/opt/scalapack/204/gnu-7.3.0/lib -Wl,-rpath=/opt/scalapack/204/gnu-7.3.0/lib -lscalapack  -L/opt/openblas/0.2.20/gnu-7.3.0/lib -Wl,-rpath=/opt/openblas/0.2.20/gnu-7.3.0/lib -lopenblas -L/home/nicpa/codes/flook/obj -lflookall -ldl -lnetcdff -lnetcdf -lpnetcdf -lhdf5_hl -lhdf5 -lz -m64 -fPIC -O3 -ftree-vectorize -fexpensive-optimizations -funroll-loops -fprefetch-loop-arrays -fgraphite -fipa-icf -fipa-pure-const -march=native -fschedule-fusion -fselective-scheduling -fstack-arrays -fmax-stack-var-size=10000 -fipa-sra -fipa-cp -fno-second-underscore -L/opt/fftw/3.3.7/gnu-7.3.0/lib -Wl,-rpath=/opt/fftw/3.3.7/gnu-7.3.0/lib -lfftw3
+PARALLEL version
+NetCDF support
+NetCDF-4 support
+NetCDF-4 MPI-IO support
+METIS ordering support
+
+* Running on 2 nodes in parallel
+>> Start of run:  24-MAY-2018  21:57:34
+
+                           ***********************       
+                           *  WELCOME TO SIESTA  *       
+                           ***********************       
+
+reinit: Reading from ../si_coop_kp_file.fdf
+
+reinit: -----------------------------------------------------------------------
+reinit: System Name: Si chain for COOP curves calculation
+reinit: -----------------------------------------------------------------------
+reinit: System Label: si_coop_kp_file
+reinit: -----------------------------------------------------------------------
+
+initatom: Reading input for the pseudopotentials and atomic orbitals ----------
+Species number:   1 Atomic number:   14 Label: Si
+
+Ground state valence configuration:   3s02  3p02
+Reading pseudopotential information in formatted form from Si.psf
+
+Valence configuration for pseudopotential generation:
+3s( 2.00) rc: 1.89
+3p( 2.00) rc: 1.89
+3d( 0.00) rc: 1.89
+4f( 0.00) rc: 1.89
+For Si, standard SIESTA heuristics set lmxkb to 3
+ (one more than the basis l, including polarization orbitals).
+Use PS.lmax or PS.KBprojectors blocks to override.
+
+<basis_specs>
+===============================================================================
+Si                   Z=  14    Mass=  28.090        Charge= 0.17977+309
+Lmxo=1 Lmxkb= 3    BasisType=split      Semic=F
+L=0  Nsemic=0  Cnfigmx=3
+          n=1  nzeta=1  polorb=0
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    0.0000    
+            lambdas:    1.0000    
+L=1  Nsemic=0  Cnfigmx=3
+          n=1  nzeta=1  polorb=1
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    0.0000    
+            lambdas:    1.0000    
+-------------------------------------------------------------------------------
+L=0  Nkbl=1  erefs: 0.17977+309
+L=1  Nkbl=1  erefs: 0.17977+309
+L=2  Nkbl=1  erefs: 0.17977+309
+L=3  Nkbl=1  erefs: 0.17977+309
+===============================================================================
+</basis_specs>
+
+atom: Called for Si                    (Z =  14)
+
+read_vps: Pseudopotential generation method:
+read_vps: ATM3      Troullier-Martins                       
+Total valence charge:    4.00000
+
+xc_check: Exchange-correlation functional:
+xc_check: Ceperley-Alder
+V l=0 = -2*Zval/r beyond r=  2.5494
+V l=1 = -2*Zval/r beyond r=  2.5494
+V l=2 = -2*Zval/r beyond r=  2.5494
+V l=3 = -2*Zval/r beyond r=  2.5494
+All V_l potentials equal beyond r=  1.8652
+This should be close to max(r_c) in ps generation
+All pots = -2*Zval/r beyond r=  2.5494
+Using large-core scheme for Vlocal
+
+atom: Estimated core radius    2.54944
+
+atom: Including non-local core corrections could be a good idea
+atom: Maximum radius for 4*pi*r*r*local-pseudopot. charge    2.85303
+atom: Maximum radius for r*vlocal+2*Zval:    2.58151
+
+KBgen: Kleinman-Bylander projectors: 
+GHOST: No ghost state for L =  0
+   l= 0   rc=  1.936440   el= -0.796617   Ekb=  4.661340   kbcos=  0.299756
+GHOST: No ghost state for L =  1
+   l= 1   rc=  1.936440   el= -0.307040   Ekb=  1.494238   kbcos=  0.301471
+GHOST: No ghost state for L =  2
+   l= 2   rc=  1.936440   el=  0.002313   Ekb= -2.808672   kbcos= -0.054903
+GHOST: No ghost state for L =  3
+   l= 3   rc=  1.936440   el=  0.003402   Ekb= -0.959059   kbcos= -0.005513
+
+KBgen: Total number of  Kleinman-Bylander projectors:   16
+atom: -------------------------------------------------------------------------
+
+atom: SANKEY-TYPE ORBITALS:
+
+SPLIT: Orbitals with angular momentum L= 0
+
+SPLIT: Basis orbitals for state 3s
+
+SPLIT: PAO cut-off radius determined from an
+SPLIT: energy shift=  0.007350 Ry
+
+   izeta = 1
+                 lambda =    1.000000
+                     rc =    5.674097
+                 energy =   -0.790139
+                kinetic =    0.533579
+    potential(screened) =   -1.323718
+       potential(ionic) =   -3.776962
+
+SPLIT: Orbitals with angular momentum L= 1
+
+SPLIT: Basis orbitals for state 3p
+
+SPLIT: PAO cut-off radius determined from an
+SPLIT: energy shift=  0.007350 Ry
+
+   izeta = 1
+                 lambda =    1.000000
+                     rc =    7.105845
+                 energy =   -0.299565
+                kinetic =    0.824289
+    potential(screened) =   -1.123854
+       potential(ionic) =   -3.348521
+
+POLgen: Perturbative polarization orbital with L=  2
+
+POLgen: Polarization orbital for state 3p
+
+   izeta = 1
+                     rc =    7.105845
+                 energy =    0.366373
+                kinetic =    1.162439
+    potential(screened) =   -0.796066
+       potential(ionic) =   -2.795335
+atom: Total number of Sankey-type orbitals:  9
+
+atm_pop: Valence configuration (for local Pseudopot. screening):
+ 3s( 2.00)                                                            
+ 3p( 2.00)                                                            
+Vna: chval, zval:    4.00000   4.00000
+
+Vna:  Cut-off radius for the neutral-atom potential:   7.105845
+
+atom: _________________________________________________________________________
+
+prinput: Basis input ----------------------------------------------------------
+
+PAO.BasisType split     
+
+%block ChemicalSpeciesLabel
+    1   14 Si                      # Species index, atomic number, species label
+%endblock ChemicalSpeciesLabel
+
+%block PAO.Basis                 # Define Basis set
+Si                    2                    # Species label, number of l-shells
+ n=3   0   1                         # n, l, Nzeta 
+   5.674   
+   1.000   
+ n=3   1   1 P   1                   # n, l, Nzeta, Polarization, NzetaPol
+   7.106   
+   1.000   
+%endblock PAO.Basis
+
+prinput: ----------------------------------------------------------------------
+
+Dumping basis to NetCDF file Si.ion.nc
+coor:   Atomic-coordinates input format  =     Cartesian coordinates
+coor:                                          (in Angstroms)
+
+siesta: Atomic coordinates (Bohr) and species
+siesta:      0.00000   0.00000   0.00000  1        1
+siesta:      1.88973   0.00000   0.00000  1        2
+siesta:      3.77945   0.00000   0.00000  1        3
+siesta:      5.66918   0.00000   0.00000  1        4
+
+siesta: System type = chain     
+
+initatomlists: Number of atoms, orbitals, and projectors:      4    36    64
+
+siesta: ******************** Simulation parameters ****************************
+siesta:
+siesta: The following are some of the parameters of the simulation.
+siesta: A complete list of the parameters used, including default values,
+siesta: can be found in file out.fdf
+siesta:
+redata: Spin configuration                          = none
+redata: Number of spin components                   = 1
+redata: Time-Reversal Symmetry                      = T
+redata: Spin-spiral                                 = F
+redata: Long output                                 =   F
+redata: Number of Atomic Species                    =        1
+redata: Charge density info will appear in .RHO file
+redata: Write Mulliken Pop.                         = NO
+redata: Matel table size (NRTAB)                    =     1024
+redata: Mesh Cutoff                                 =   150.0000 Ry
+redata: Net charge of the system                    =     0.0000 |e|
+redata: Min. number of SCF Iter                     =        0
+redata: Max. number of SCF Iter                     =      500
+redata: SCF convergence failure will abort job
+redata: SCF mix quantity                            = Hamiltonian
+redata: Mix DM or H after convergence               =   F
+redata: Recompute H after scf cycle                 =   F
+redata: Mix DM in first SCF step                    =   T
+redata: Write Pulay info on disk                    =   F
+redata: New DM Mixing Weight                        =     0.1000
+redata: New DM Occupancy tolerance                  = 0.000000000001
+redata: No kicks to SCF
+redata: DM Mixing Weight for Kicks                  =     0.5000
+redata: Require Harris convergence for SCF          =   F
+redata: Harris energy tolerance for SCF             =     0.000100 eV
+redata: Require DM convergence for SCF              =   T
+redata: DM tolerance for SCF                        =     0.000100
+redata: Require EDM convergence for SCF             =   F
+redata: EDM tolerance for SCF                       =     0.001000 eV
+redata: Require H convergence for SCF               =   T
+redata: Hamiltonian tolerance for SCF               =     0.001000 eV
+redata: Require (free) Energy convergence for SCF   =   F
+redata: (free) Energy tolerance for SCF             =     0.000100 eV
+redata: Using Saved Data (generic)                  =   F
+redata: Use continuation files for DM               =   F
+redata: Neglect nonoverlap interactions             =   F
+redata: Method of Calculation                       = Diagonalization
+redata: Electronic Temperature                      =   290.1109 K
+redata: Fix the spin of the system                  =   F
+redata: Max. number of TDED Iter                    =        1
+redata: Number of TDED substeps                     =        3
+redata: Dynamics option                             = Single-point calculation
+mix.SCF: Pulay mixing                            = Pulay
+mix.SCF:    Variant                              = stable
+mix.SCF:    History steps                        = 3
+mix.SCF:    Linear mixing weight                 =     0.100000
+mix.SCF:    Mixing weight                        =     0.100000
+mix.SCF:    SVD condition                        = 0.1000E-07
+redata: Save all siesta data in one NC              =   F
+redata: ***********************************************************************
+
+%block SCF.Mixers
+  Pulay
+%endblock SCF.Mixers
+
+%block SCF.Mixer.Pulay
+  # Mixing method
+  method pulay
+  variant stable
+
+  # Mixing options
+  weight 0.1000
+  weight.linear 0.1000
+  history 3
+%endblock SCF.Mixer.Pulay
+
+DM_history_depth set to one: no extrapolation allowed by default for geometry relaxation
+Size of DM history Fstack: 1
+Total number of electrons:    16.000000
+Total ionic charge:    16.000000
+
+* ProcessorY, Blocksize:    1  19
+
+
+* Orbital distribution balance (max,min):    19    17
+
+
+siesta: k-grid: Number of k-points =         5
+siesta: k-points from user-defined list
+
+siesta: PDOS k-grid: Number of k-points =         5
+siesta: PDOS k-points from user-defined list
+
+diag: Algorithm                                     = D&C
+diag: Parallel over k                               =   F
+diag: Use parallel 2D distribution                  =   F
+diag: Parallel block-size                           = 19
+diag: Parallel distribution                         =     1 x     2
+diag: Used triangular part                          = Lower
+diag: Absolute tolerance                            =  0.100E-15
+diag: Orthogonalization factor                      =  0.100E-05
+diag: Memory factor                                 =  1.0000
+
+superc: Internal auxiliary supercell:     5 x     1 x     1  =       5
+superc: Number of atoms, orbitals, and projectors:     20    180    320
+
+
+ts: **************************************************************
+ts: Save H and S matrices                           =    F
+ts: Save DM and EDM matrices                        =    F
+ts: Fix Hartree potential                           =    F
+ts: Only save the overlap matrix S                  =    F
+ts: **************************************************************
+
+************************ Begin: TS CHECKS AND WARNINGS ************************
+************************ End: TS CHECKS AND WARNINGS **************************
+
+
+                     ====================================
+                        Single-point calculation
+                     ====================================
+
+superc: Internal auxiliary supercell:     5 x     1 x     1  =       5
+superc: Number of atoms, orbitals, and projectors:     20    180    320
+
+outcell: Unit cell vectors (Ang):
+        4.000000    0.000000    0.000000
+        0.000000   10.000000    0.000000
+        0.000000    0.000000   10.000000
+
+outcell: Cell vector modules (Ang)   :    4.000000   10.000000   10.000000
+outcell: Cell angles (23,13,12) (deg):     90.0000     90.0000     90.0000
+outcell: Cell volume (Ang**3)        :    400.0000
+<dSpData1D:S at geom step 0
+  <sparsity:sparsity for geom step 0
+    nrows_g=36 nrows=19 sparsity=2.2431 nnzs=2907, refcount: 7>
+  <dData1D:(new from dSpData1D) n=2907, refcount: 1>
+refcount: 1>
+new_DM -- step:     1
+Initializing Density Matrix...
+DM filled with atomic data:
+<dSpData2D:DM initialized from atoms
+  <sparsity:sparsity for geom step 0
+    nrows_g=36 nrows=19 sparsity=2.2431 nnzs=2907, refcount: 8>
+  <dData2D:DM n=2907 m=1, refcount: 1>
+refcount: 1>
+No. of atoms with KB's overlaping orbs in proc 0. Max # of overlaps:      11      81
+New grid distribution:   1
+           1       1:   15    1:   40    1:   20
+           2       1:   15    1:   40   21:   40
+
+InitMesh: MESH =    30 x    80 x    80 =      192000
+InitMesh: (bp) =    15 x    40 x    40 =       24000
+InitMesh: Mesh cutoff (required, used) =   150.000   155.462 Ry
+ExtMesh (bp) on 0 =    75 x   104 x    84 =      655200
+New grid distribution:   2
+           1       1:   15    1:   40    1:   13
+           2       1:   15    1:   40   14:   40
+New grid distribution:   3
+           1       1:   15    1:   40    1:   15
+           2       1:   15    1:   40   16:   40
+Setting up quadratic distribution...
+ExtMesh (bp) on 0 =    75 x   104 x    77 =      600600
+PhiOnMesh: Number of (b)points on node 0 =                 7800
+PhiOnMesh: nlist on node 0 =               242164
+cdiag-debug: jobz=V, algo= 1, Node=    1, work=        2090, rwork=        2161, iwork=         270
+cdiag-debug: jobz=V, algo= 1, Node=    0, work=        2090, rwork=        2377, iwork=         270
+
+stepf: Fermi-Dirac step function
+
+siesta: Program's energy decomposition (eV):
+siesta: Ebs     =      -226.878967
+siesta: Eions   =       761.604247
+siesta: Ena     =       338.976744
+siesta: Ekin    =       265.845531
+siesta: Enl     =       155.324530
+siesta: Eso     =         0.000000
+siesta: Eldau   =         0.000000
+siesta: DEna    =       -17.516845
+siesta: DUscf   =         2.038399
+siesta: DUext   =         0.000000
+siesta: Enegf   =         0.000000
+siesta: Exc     =      -157.428404
+siesta: eta*DQ  =         0.000000
+siesta: Emadel  =         0.000000
+siesta: Emeta   =         0.000000
+siesta: Emolmec =         0.000000
+siesta: Ekinion =         0.000000
+siesta: Eharris =      -121.207230
+siesta: Etot    =      -174.364293
+siesta: FreeEng =      -174.372957
+
+        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV) dHmax(eV)
+   scf:    1     -121.207230     -174.364293     -174.372957  6.055278 -7.994472  2.893698
+timer: Routine,Calls,Time,% = IterSCF        1       0.340  30.58
+   scf:    2     -174.475630     -174.420592     -174.429257  0.109988 -7.666230  2.507378
+   scf:    3     -174.625748     -174.566213     -174.581766  0.701077 -5.533050  0.190181
+   scf:    4     -174.159092     -174.391369     -174.405428  0.352909 -4.746563  1.785466
+   scf:    5     -174.678613     -174.581704     -174.592089  0.348919 -5.541164  0.027093
+   scf:    6     -174.582277     -174.581998     -174.592127  0.006400 -5.561870  0.001956
+   scf:    7     -174.581976     -174.581987     -174.592127  0.000454 -5.561077  0.000760
+   scf:    8     -174.582000     -174.581993     -174.592127  0.000420 -5.561103  0.000663
+   scf:    9     -174.582018     -174.582006     -174.592127  0.000331 -5.561216  0.000500
+   scf:   10     -174.582041     -174.582023     -174.592128  0.000337 -5.561401  0.000441
+   scf:   11     -174.582026     -174.582025     -174.592128  0.000019 -5.561412  0.000442
+
+SCF Convergence by DM+H criterion
+max |DM_out - DM_in|         :     0.0000187439
+max |H_out - H_in|      (eV) :     0.0004420538
+SCF cycle converged after 11 iterations
+
+Using DM_out to compute the final energy and forces
+No. of atoms with KB's overlaping orbs in proc 0. Max # of overlaps:      11      81
+
+siesta: E_KS(eV) =             -174.5820
+
+siesta: E_KS - E_eggbox =      -174.5820
+
+siesta: Atomic forces (eV/Ang):
+----------------------------------------
+   Tot    0.000000   -0.000000    0.000000
+----------------------------------------
+   Max    0.000000
+   Res    0.000000    sqrt( Sum f_i^2 / 3N )
+----------------------------------------
+   Max    0.000000    constrained
+
+Stress-tensor-Voigt (kbar):    -3042.14       -0.01       -0.01        0.00        0.00       -0.00
+(Free)E + p*V (eV/cell)       78.5732
+Target enthalpy (eV/cell)     -174.5921
+Writing WFSX for COOP/COHP in si_coop_kp_file.fullBZ.WFSX
+siesta: PDOS info: 
+siesta: e1, e2, sigma, nhist:   -25.00 eV    5.00 eV    0.20 eV    500
+
+siesta: Program's energy decomposition (eV):
+siesta: Ebs     =      -236.953982
+siesta: Eions   =       761.604247
+siesta: Ena     =       338.976744
+siesta: Ekin    =       261.043605
+siesta: Enl     =       149.913747
+siesta: Eso     =         0.000000
+siesta: Eldau   =         0.000000
+siesta: DEna    =        -8.184581
+siesta: DUscf   =         1.510262
+siesta: DUext   =         0.000000
+siesta: Enegf   =         0.000000
+siesta: Exc     =      -156.237555
+siesta: eta*DQ  =         0.000000
+siesta: Emadel  =         0.000000
+siesta: Emeta   =         0.000000
+siesta: Emolmec =         0.000000
+siesta: Ekinion =         0.000000
+siesta: Eharris =      -174.582026
+siesta: Etot    =      -174.582025
+siesta: FreeEng =      -174.592128
+
+siesta: Final energy (eV):
+siesta:  Band Struct. =    -236.953982
+siesta:       Kinetic =     261.043605
+siesta:       Hartree =     829.334283
+siesta:       Eldau   =       0.000000
+siesta:       Eso     =       0.000000
+siesta:    Ext. field =       0.000000
+siesta:       Enegf   =       0.000000
+siesta:   Exch.-corr. =    -156.237555
+siesta:  Ion-electron =   -1952.658160
+siesta:       Ion-ion =     843.935803
+siesta:       Ekinion =       0.000000
+siesta:         Total =    -174.582025
+siesta:         Fermi =      -5.561412
+
+siesta: Stress tensor (static) (eV/Ang**3):
+siesta:    -1.898732    0.000000    0.000000
+siesta:     0.000000   -0.000004   -0.000000
+siesta:    -0.000000    0.000000   -0.000004
+
+siesta: Cell volume =        400.000000 Ang**3
+
+siesta: Pressure (static):
+siesta:                Solid            Molecule  Units
+siesta:           0.00689323          0.00689323  Ry/Bohr**3
+siesta:           0.63291323          0.63291323  eV/Ang**3
+siesta:        1014.04975274       1014.04975274  kBar
+(Free)E+ p_basis*V_orbitals  =        -171.801899
+(Free)Eharris+ p_basis*V_orbitals  =        -171.801901
+
+siesta: Electric dipole (a.u.)  =   -0.000000    0.000000    0.000000
+siesta: Electric dipole (Debye) =   -0.000000    0.000000    0.000000
+>> End of run:  24-MAY-2018  21:57:38
+Job completed
 
=== added directory 'Tests/si_coop_kp_file'
=== added file 'Tests/si_coop_kp_file/makefile'
--- Tests/si_coop_kp_file/makefile	1970-01-01 00:00:00 +0000
+++ Tests/si_coop_kp_file/makefile	2018-05-27 11:17:14 +0000
@@ -0,0 +1,7 @@
+#
+# Single-test makefile
+#
+name=si_coop_kp_file
+EXTRAFILES=si_coop_kp_file.KP.input
+#
+include ../test.mk
 
=== added file 'Tests/si_coop_kp_file/si_coop_kp_file.KP.input'
--- Tests/si_coop_kp_file/si_coop_kp_file.KP.input	1970-01-01 00:00:00 +0000
+++ Tests/si_coop_kp_file/si_coop_kp_file.KP.input	2018-05-27 11:17:14 +0000
@@ -0,0 +1,6 @@
+         5
+         1  0.000000E+00  0.000000E+00  0.000000E+00   0.125000E+00
+         2  0.125000E+00  0.000000E+00  0.000000E+00   0.250000E+00
+         3  0.250000E+00  0.000000E+00  0.000000E+00   0.250000E+00
+         4  0.375000E+00  0.000000E+00  0.000000E+00   0.250000E+00
+         5  0.500000E+00  0.000000E+00  0.000000E+00   0.125000E+00
 
=== added file 'Tests/si_coop_kp_file/si_coop_kp_file.fdf'
--- Tests/si_coop_kp_file/si_coop_kp_file.fdf	1970-01-01 00:00:00 +0000
+++ Tests/si_coop_kp_file/si_coop_kp_file.fdf	2018-05-27 11:17:14 +0000
@@ -0,0 +1,41 @@
+SystemName  Si chain for COOP curves calculation
+SystemLabel si_coop_kp_file
+NumberOfAtoms       4
+NumberOfSpecies     1
+%block ChemicalSpeciesLabel
+ 1  14  Si
+%endblock ChemicalSpeciesLabel
+#------------6.3 BASIS DEFINITION---------------------------------------
+PAO.BasisSize      SZP
+PAO.EnergyShift   100 meV
+#------------6.4 LATTICE, COORDINATES -----------------------------------
+LatticeConstant     1.0000 Ang
+%block LatticeVectors
+ 4.000  0.000  0.000
+ 0.000  10.0  0.000
+ 0.000  0.000  10.0
+%endblock LatticeVectors
+AtomicCoordinatesFormat Ang
+%block AtomicCoordinatesAndAtomicSpecies
+ 0.0000   0.0000   0.0000  1
+ 1.0000   0.0000   0.0000  1
+ 2.0000   0.0000   0.0000  1
+ 3.0000   0.0000   0.0000  1
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+# User-defined k-point
+kgrid.File si_coop_kp_file.KP.input
+
+MeshCutoff         150.0 Ry
+MaxSCFIterations    500
+DM.MixingWeight      0.1
+DM.NumberPulay       3
+DM.Tolerance         1.d-4
+
+SolutionMethod       diagon
+ElectronicTemperature  25  meV
+
+COOP.write T
+%block ProjectedDensityOfStates
+ -25.  5.  0.2  500   eV
+%endblock ProjectedDensityOfStates
 
=== added file 'Tests/si_coop_kp_file/si_coop_kp_file.pseudos'
--- Tests/si_coop_kp_file/si_coop_kp_file.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/si_coop_kp_file/si_coop_kp_file.pseudos	2018-05-27 11:17:14 +0000
@@ -0,0 +1,1 @@
+Si
 
=== modified file 'Util/COOP/Makefile'
--- Util/COOP/Makefile	2018-04-07 19:24:04 +0000
+++ Util/COOP/Makefile	2018-05-27 11:17:14 +0000
@@ -137,16 +137,17 @@
 class_Sparsity.o: alloc.o
 coceri.o: files.o periodic_table.o precision.o units.o
 compute_dm.o: atomlist.o class_SpData1D.o compute_ebs_shift.o diagon.o
-compute_dm.o: iodmhs_netcdf.o kpoint_grid.o local_sys.o m_chess.o m_dminim.o
+compute_dm.o: iodmhs_netcdf.o kpoint_scf.o local_sys.o m_chess.o m_dminim.o
 compute_dm.o: m_energies.o m_eo.o m_gamma.o m_hsx.o m_pexsi_driver.o m_rmaxh.o
 compute_dm.o: m_spin.o m_steps.o m_transiesta.o m_ts_global_vars.o m_zminim.o
 compute_dm.o: normalize_dm.o ordern.o parallel.o precision.o siesta_geom.o
 compute_dm.o: siesta_options.o sparse_matrices.o units.o
 compute_ebs_shift.o: m_mpi_utils.o parallel.o precision.o
-compute_energies.o: atomlist.o class_SpData1D.o class_SpData2D.o dhscf.o
-compute_energies.o: files.o m_dipol.o m_energies.o m_mpi_utils.o m_ntm.o
-compute_energies.o: m_rhog.o m_spin.o precision.o siesta_geom.o
-compute_energies.o: siesta_options.o sparse_matrices.o
+compute_energies.o: atomlist.o class_SpData1D.o class_SpData2D.o
+compute_energies.o: class_SpData2D.o dhscf.o files.o m_dipol.o m_energies.o
+compute_energies.o: m_mpi_utils.o m_ntm.o m_rhog.o m_spin.o parallel.o
+compute_energies.o: precision.o siesta_geom.o siesta_options.o
+compute_energies.o: sparse_matrices.o
 compute_max_diff.o: m_mpi_utils.o precision.o
 compute_norm.o: m_mpi_utils.o m_spin.o precision.o sparse_matrices.o
 compute_pw_matrix.o: alloc.o m_planewavematrix.o m_planewavematrixvar.o
@@ -164,7 +165,7 @@
 cranknic_evolg.o: m_inversemm.o m_matswinvers.o m_spin.o m_steps.o parallel.o
 cranknic_evolg.o: parallelsubs.o precision.o siesta_options.o sparse_matrices.o
 cranknic_evolg.o: units.o wavefunctions.o
-cranknic_evolk.o: atomlist.o cranknic_evolg.o kpoint_grid.o m_energies.o m_eo.o
+cranknic_evolk.o: atomlist.o cranknic_evolg.o kpoint_scf.o m_energies.o m_eo.o
 cranknic_evolk.o: m_spin.o parallel.o parallelsubs.o precision.o
 cranknic_evolk.o: siesta_options.o sparse_matrices.o units.o wavefunctions.o
 create_Sparsity_SC.o: class_Sparsity.o geom_helper.o intrinsic_missing.o
@@ -178,9 +179,8 @@
 denmatlomem.o: alloc.o globalise.o onmod.o precision.o
 densematrix.o: alloc.o precision.o
 detover.o: alloc.o linpack.o parallel.o parallelsubs.o precision.o
-dfscf.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o local_sys.o
-dfscf.o: m_spin.o mesh.o meshdscf.o meshphi.o parallel.o parallelsubs.o
-dfscf.o: precision.o
+dfscf.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o local_sys.o mesh.o
+dfscf.o: meshdscf.o meshphi.o parallel.o parallelsubs.o precision.o
 dhscf.o: alloc.o atmfuncs.o bsc_xcmod.o cellxc_mod.o delk.o dfscf.o
 dhscf.o: doping_uniform.o files.o forhar.o iogrid_netcdf.o local_sys.o
 dhscf.o: m_charge_add.o m_efield.o m_hartree_add.o m_iorho.o m_iotddft.o
@@ -236,15 +236,15 @@
 fft.o: alloc.o fft1d.o local_sys.o m_timer.o mesh.o parallel.o parallelsubs.o
 fft.o: precision.o
 fft1d.o: local_sys.o parallel.o precision.o
-final_H_f_stress.o: alloc.o atomlist.o class_SpData2D.o compute_max_diff.o
-final_H_f_stress.o: dnaefs.o files.o grdsam.o kinefsm.o ldau.o ldau_specs.o
-final_H_f_stress.o: local_sys.o m_dipol.o m_energies.o m_forces.o m_gamma.o
-final_H_f_stress.o: m_hsx.o m_mpi_utils.o m_ncdf_siesta.o m_ntm.o m_spin.o
-final_H_f_stress.o: m_steps.o m_stress.o m_ts_global_vars.o m_ts_io.o
-final_H_f_stress.o: m_ts_kpoints.o m_ts_options.o metaforce.o
-final_H_f_stress.o: molecularmechanics.o naefs.o nlefsm.o overfsm.o parallel.o
-final_H_f_stress.o: siesta_geom.o siesta_options.o sparse_matrices.o
-final_H_f_stress.o: spinorbit.o units.o
+final_H_f_stress.o: alloc.o atomlist.o class_SpData2D.o class_SpData2D.o
+final_H_f_stress.o: compute_max_diff.o dnaefs.o files.o grdsam.o kinefsm.o
+final_H_f_stress.o: ldau.o ldau_specs.o local_sys.o m_dipol.o m_energies.o
+final_H_f_stress.o: m_forces.o m_gamma.o m_hsx.o m_mpi_utils.o m_ncdf_siesta.o
+final_H_f_stress.o: m_ntm.o m_spin.o m_steps.o m_stress.o m_ts_global_vars.o
+final_H_f_stress.o: m_ts_io.o m_ts_options.o metaforce.o molecularmechanics.o
+final_H_f_stress.o: naefs.o nlefsm.o overfsm.o parallel.o siesta_geom.o
+final_H_f_stress.o: siesta_options.o sparse_matrices.o spinorbit.o
+final_H_f_stress.o: ts_kpoint_scf.o units.o
 find_kgrid.o: alloc.o minvec.o parallel.o precision.o units.o
 fire_optim.o: alloc.o m_fire.o m_mpi_utils.o parallel.o precision.o
 fire_optim.o: siesta_options.o units.o
@@ -263,8 +263,8 @@
 grdsam.o: m_spin.o parallel.o precision.o siesta_geom.o siesta_options.o
 grdsam.o: units.o
 hsparse.o: alloc.o atm_types.o atmfuncs.o atomlist.o ldau_specs.o listsc.o
-hsparse.o: local_sys.o mneighb.o parallel.o parallelsubs.o precision.o radial.o
-hsparse.o: siesta_options.o sorting.o sparse_matrices.o
+hsparse.o: mneighb.o parallel.o parallelsubs.o precision.o radial.o sorting.o
+hsparse.o: sparse_matrices.o
 idiag.o: local_sys.o parallel.o
 init_output.o: files.o
 initatom.o: atm_types.o atmparams.o atom.o atom_options.o basis_io.o
@@ -298,10 +298,10 @@
 kinefsm.o: alloc.o atmfuncs.o mneighb.o new_matel.o parallel.o parallelsubs.o
 kinefsm.o: precision.o
 kpoint_convert.o: local_sys.o precision.o units.o
-kpoint_grid.o: find_kgrid.o local_sys.o m_spin.o minvec.o parallel.o
-kpoint_grid.o: precision.o siesta_cml.o siesta_options.o units.o
-kpoint_pdos.o: find_kgrid.o local_sys.o m_spin.o minvec.o parallel.o
-kpoint_pdos.o: precision.o siesta_options.o units.o
+kpoint_pdos.o: kpoint_t.o m_spin.o parallel.o precision.o siesta_options.o
+kpoint_scf.o: kpoint_t.o m_spin.o parallel.o precision.o siesta_options.o
+kpoint_t.o: alloc.o files.o find_kgrid.o m_char.o m_io.o m_os.o minvec.o
+kpoint_t.o: parallel.o precision.o siesta_cml.o units.o
 ksv.o: alloc.o atmfuncs.o densematrix.o ksvinit.o local_sys.o parallel.o
 ksv.o: precision.o
 ksvinit.o: alloc.o local_sys.o parallel.o precision.o
@@ -312,7 +312,7 @@
 ldau_specs.o: basis_specs.o basis_types.o interpolation.o local_sys.o m_cite.o
 ldau_specs.o: parallel.o precision.o pseudopotential.o radial.o units.o
 listsc.o: alloc.o
-local_DOS.o: atomlist.o dhscf.o diagon.o files.o kpoint_grid.o local_sys.o
+local_DOS.o: atomlist.o dhscf.o diagon.o files.o kpoint_scf.o local_sys.o
 local_DOS.o: m_energies.o m_eo.o m_forces.o m_gamma.o m_ntm.o m_spin.o
 local_DOS.o: parallel.o siesta_geom.o siesta_options.o sparse_matrices.o
 m_broyddj.o: alloc.o local_sys.o m_mpi_utils.o parallel.o precision.o
@@ -362,7 +362,7 @@
 m_hsx.o: atm_types.o atmfuncs.o atomlist.o files.o local_sys.o parallel.o
 m_hsx.o: parallelsubs.o precision.o siesta_geom.o
 m_initwf.o: alloc.o atomlist.o densematrix.o diag.o diag_option.o fermid.o
-m_initwf.o: kpoint_grid.o local_sys.o m_eo.o m_gamma.o m_memory.o m_spin.o
+m_initwf.o: kpoint_scf.o local_sys.o m_eo.o m_gamma.o m_memory.o m_spin.o
 m_initwf.o: parallel.o parallelsubs.o precision.o sparse_matrices.o
 m_initwf.o: wavefunctions.o
 m_integrate.o: precision.o
@@ -395,19 +395,19 @@
 m_ncdf_io.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 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_gamma.o m_kinetic.o
+m_ncdf_siesta.o: kpoint_scf.o m_energies.o m_forces.o m_gamma.o m_kinetic.o
 m_ncdf_siesta.o: m_ncdf_io.o m_ntm.o m_spin.o m_stress.o m_ts_electype.o
-m_ncdf_siesta.o: m_ts_kpoints.o m_ts_options.o parallel.o precision.o radial.o
-m_ncdf_siesta.o: siesta_geom.o siesta_options.o sparse_matrices.o timestamp.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 ts_kpoint_scf.o
 m_new_dm.o: alloc.o atomlist.o class_Data2D.o class_Fstack_Data1D.o
 m_new_dm.o: class_Fstack_Pair_Geometry_SpData2D.o class_Geometry.o
 m_new_dm.o: class_OrbitalDistribution.o class_Pair_Geometry_SpData2D.o
 m_new_dm.o: class_SpData2D.o class_Sparsity.o files.o local_sys.o m_energies.o
-m_new_dm.o: m_handle_sparse.o m_iodm.o m_mixing.o m_mixing_scf.o m_mpi_utils.o
-m_new_dm.o: m_spin.o m_spin.o m_steps.o m_svd.o m_ts_electype.o
-m_new_dm.o: m_ts_global_vars.o m_ts_iodm.o m_ts_method.o m_ts_options.o
-m_new_dm.o: parallel.o parsing.o precision.o restructSpData2D.o siesta_geom.o
-m_new_dm.o: siesta_options.o sparse_matrices.o units.o
+m_new_dm.o: m_handle_sparse.o m_iodm.o m_mixing.o m_mixing_scf.o m_spin.o
+m_new_dm.o: m_spin.o m_steps.o m_svd.o m_ts_electype.o m_ts_global_vars.o
+m_new_dm.o: m_ts_iodm.o m_ts_method.o m_ts_options.o parallel.o parsing.o
+m_new_dm.o: precision.o restructSpData2D.o siesta_geom.o siesta_options.o
+m_new_dm.o: sparse_matrices.o units.o
 m_noccbands.o: alloc.o atmfuncs.o local_sys.o m_spin.o parallel.o precision.o
 m_noccbands.o: siesta_geom.o
 m_options.o: precision.o
@@ -449,7 +449,7 @@
 m_sparsity_handling.o: class_SpData2D.o class_Sparsity.o geom_helper.o
 m_sparsity_handling.o: intrinsic_missing.o m_interpolate.o m_region.o
 m_sparsity_handling.o: precision.o
-m_spin.o: alloc.o local_sys.o parallel.o precision.o units.o
+m_spin.o: alloc.o files.o local_sys.o m_cite.o parallel.o precision.o units.o
 m_stress.o: precision.o
 m_supercell.o: atom_graph.o class_OrbitalDistribution.o class_SpData2D.o
 m_supercell.o: intrinsic_missing.o parallel.o parallelsubs.o precision.o
@@ -463,9 +463,9 @@
 m_transiesta.o: class_SpData2D.o class_Sparsity.o files.o m_interpolate.o
 m_transiesta.o: m_ts_charge.o m_ts_contour_eq.o m_ts_contour_neq.o
 m_transiesta.o: m_ts_electype.o m_ts_fullg.o m_ts_fullk.o m_ts_gf.o
-m_transiesta.o: m_ts_kpoints.o m_ts_method.o m_ts_mumpsg.o m_ts_mumpsk.o
-m_transiesta.o: m_ts_options.o m_ts_sparse.o m_ts_tri_common.o m_ts_tri_init.o