Merge lp:~nickpapior/siesta/tddft-work into lp:~rraffiu/siesta/tddft02-format0.92

[Nick Papior]

Fixes for running TDDFT branch on my box, reduced memory requirement

Made TDDFT branch work, and reduced memory.

I think this branch is a requirement to get it to work.

[Rafi Ullah (rraffiu)]

review: Approve

=== modified file 'Src/Makefile'
--- Src/Makefile	2017-10-03 13:27:08 +0000
+++ Src/Makefile	2017-10-03 17:43:24 +0000
@@ -583,9 +583,6 @@
 cgvc_zmatrix.o: alloc.o conjgr.o m_mpi_utils.o parallel.o precision.o sys.o
 cgvc_zmatrix.o: units.o zmatrix.o
 cgwf.o: egandd.o onmod.o parallel.o precision.o sys.o
-sankey_change_basis.o: alloc.o m_matdiag.o parallel.o parallelsubs.o precision.o
-sankey_change_basis.o: sparse_matrices.o sys.o wavefunctions.o kpoint_grid.o
-sankey_change_basis.o: m_gamma.o m_spin.o atomlist.o diag_option.o
 chemical.o: parallel.o precision.o sys.o
 chempot.o: alloc.o mneighb.o parallel.o parallelsubs.o precision.o sys.o
 chkdim.o: sys.o
@@ -648,6 +645,13 @@
 conjgr_old.o: precision.o
 coor.o: alloc.o parallel.o precision.o siesta_geom.o sys.o units.o zmatrix.o
 coxmol.o: files.o periodic_table.o precision.o
+cranknic_evolg.o: alloc.o atomlist.o m_energies.o m_eo.o m_inversemm.o
+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: 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
 create_Sparsity_Union.o: class_OrbitalDistribution.o class_Sparsity.o
 create_Sparsity_Union.o: m_region.o parallel.o precision.o
@@ -698,7 +702,7 @@
 domain_decom.o: sparse_matrices.o sys.o
 doping_uniform.o: alloc.o m_ntm.o mesh.o parallel.o precision.o sys.o
 dynamics.o: alloc.o atomlist.o files.o ioxv.o m_mpi_utils.o parallel.o
-dynamics.o: precision.o siesta_options.o sys.o units.o
+dynamics.o: precision.o sys.o units.o
 egandd.o: alloc.o denmat.o ener3.o globalise.o gradient.o on_subs.o onmod.o
 egandd.o: onmod.o precision.o sys.o
 eggbox.o: parallel.o precision.o
@@ -706,12 +710,6 @@
 electrostatic.o: precision.o radfft.o radial.o sys.o
 ener3.o: alloc.o globalise.o m_mpi_utils.o onmod.o precision.o
 ener3lomem.o: alloc.o globalise.o m_mpi_utils.o onmod.o precision.o
-cranknic_evolg.o: alloc.o m_eo.o m_matswinvers.o m_steps.o parallel.o parallelsubs.o
-cranknic_evolg.o: precision.o siesta_options.o sparse_matrices.o sys.o wavefunctions.o
-cranknic_evolg.o: m_spin.o atomlist.o m_energies.o units.o m_inversemm.o
-cranknic_evolk.o: parallel.o parallelsubs.o precision.o siesta_options.o m_spin.o
-cranknic_evolk.o: wavefunctions.o sparse_matrices.o atomlist.o kpoint_grid.o
-cranknic_evolk.o: units.o m_energies.o m_eo.o cranknic_evolg.o
 extrapolateSpData2D.o: class_Data2D.o class_OrbitalDistribution.o
 extrapolateSpData2D.o: class_SpData2D.o class_Sparsity.o restructSpData2D.o
 extrapolon.o: parallel.o precision.o sys.o
@@ -768,9 +766,9 @@
 iomd.o: files.o precision.o
 iopipes.o: parallel.o precision.o sys.o
 iosockets.o: cellsubs.o fsockets.o m_mpi_utils.o parallel.o precision.o sys.o
+iotdxv.o: files.o parallel.o precision.o
 iowfs_netcdf.o: alloc.o parallel.o parallelsubs.o precision.o sys.o
 ioxv.o: files.o parallel.o precision.o
-iotdxv.o: files.o parallel.o precision.o
 iozm.o: files.o parallel.o precision.o siesta_geom.o zmatrix.o
 ipack.o: sys.o
 kgrid.o: parallel.o precision.o units.o
@@ -815,9 +813,8 @@
 m_efield.o: siesta_geom.o sys.o units.o
 m_energies.o: precision.o
 m_eo.o: precision.o
-m_evolve.o: alloc.o m_memory.o parallel.o parallelsubs.o precision.o sys.o
-m_evolve.o: atomlist.o m_spin.o m_gamma.o sparse_matrices.o kpoint_grid.o
-m_evolve.o: cranknic_evolg.o cranknic_evolk.o
+m_evolve.o: cranknic_evolg.o cranknic_evolk.o m_gamma.o m_spin.o precision.o
+m_evolve.o: sys.o
 m_exp_coord.o: files.o m_os.o parallel.o precision.o units.o
 m_filter.o: bessph.o precision.o radfft.o sorting.o sys.o
 m_fire.o: parallel.o precision.o
@@ -840,10 +837,9 @@
 m_hs_matrix.o: alloc.o cellsubs.o geom_helper.o precision.o sys.o
 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 densematrix.o fermid.o diag.o m_eo.o m_memory.o
-m_initwf.o: parallel.o parallelsubs.o precision.o sparse_matrices.o sys.o
-m_initwf.o: wavefunctions.o kpoint_grid.o atomlist.o m_spin.o m_gamma.o
-m_initwf.o: diag_option.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: parallelsubs.o precision.o sparse_matrices.o sys.o wavefunctions.o
 m_integrate.o: precision.o
 m_inversemm.o: precision.o
 m_io.o: sys.o
@@ -857,8 +853,8 @@
 m_iorho.o: alloc.o parallel.o parallelsubs.o precision.o sys.o
 m_iostruct.o: alloc.o files.o m_mpi_utils.o parallel.o precision.o
 m_iostruct.o: siesta_geom.o sys.o units.o
-m_iotddft.o: files.o m_dipol.o m_steps.o parallel.o siesta_options.o
-m_iotddft.o: wavefunctions.o units.o m_io.o
+m_iotddft.o: files.o m_dipol.o m_io.o m_steps.o parallel.o siesta_options.o
+m_iotddft.o: units.o wavefunctions.o
 m_kinetic.o: precision.o
 m_mat_invert.o: intrinsic_missing.o m_pivot_array.o precision.o
 m_matdiag.o: precision.o
@@ -1207,6 +1203,9 @@
 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: 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
 save_density_matrix.o: m_matio.o m_ncdf_siesta.o m_spin.o m_steps.o
 save_density_matrix.o: m_ts_global_vars.o m_ts_iodm.o m_ts_options.o
@@ -1258,17 +1257,17 @@
 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 files.o
 siesta_forces.o: final_H_f_stress.o flook_siesta.o kpoint_grid.o
-siesta_forces.o: m_check_walltime.o m_convergence.o m_energies.o m_eo.o
-siesta_forces.o: m_forces.o m_gamma.o m_initwf.o m_iodm_old.o m_mixing.o
-siesta_forces.o: m_mixing_scf.o m_mpi_utils.o m_ncdf_siesta.o m_pexsi.o
-siesta_forces.o: m_pexsi_driver.o m_rhog.o m_spin.o m_steps.o m_stress.o
-siesta_forces.o: m_transiesta.o m_ts_charge.o m_ts_electype.o
-siesta_forces.o: m_ts_global_vars.o m_ts_method.o m_ts_options.o mixer.o
-siesta_forces.o: parallel.o post_scf_work.o precision.o save_density_matrix.o
-siesta_forces.o: scfconvergence_test.o setup_H0.o setup_hamiltonian.o
-siesta_forces.o: siesta_cml.o siesta_dicts.o siesta_geom.o siesta_master.o
-siesta_forces.o: siesta_options.o sparse_matrices.o state_analysis.o
-siesta_forces.o: state_init.o sys.o timer.o units.o write_subs.o
+siesta_forces.o: m_check_walltime.o m_convergence.o m_energies.o m_forces.o
+siesta_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
+siesta_forces.o: m_ts_charge.o m_ts_electype.o m_ts_global_vars.o m_ts_method.o
+siesta_forces.o: m_ts_options.o mixer.o parallel.o post_scf_work.o precision.o
+siesta_forces.o: save_density_matrix.o scfconvergence_test.o setup_H0.o
+siesta_forces.o: setup_hamiltonian.o siesta_cml.o siesta_dicts.o siesta_geom.o
+siesta_forces.o: siesta_master.o siesta_options.o sparse_matrices.o
+siesta_forces.o: state_analysis.o state_init.o sys.o timer.o units.o
+siesta_forces.o: write_subs.o
 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
@@ -1280,7 +1279,7 @@
 siesta_init.o: parallel.o parallelsubs.o projected_DOS.o siesta_cmlsubs.o
 siesta_init.o: siesta_dicts.o siesta_geom.o siesta_options.o sparse_matrices.o
 siesta_init.o: struct_init.o sys.o timer.o timestamp.o ts_init.o writewave.o
-siesta_init.o: zmatrix.o 
+siesta_init.o: zmatrix.o
 siesta_master.o: iopipes.o iosockets.o precision.o sys.o
 siesta_move.o: atomlist.o broyden_optim.o cell_broyden_optim.o
 siesta_move.o: cell_fire_optim.o dynamics.o fire_optim.o flook_siesta.o ioxv.o
@@ -1292,10 +1291,10 @@
 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: 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 setup_H0.o setup_hamiltonian.o
-siesta_tddft.o: siesta_cml.o siesta_options.o sparse_matrices.o
-siesta_tddft.o: state_analysis.o state_init.o sys.o wavefunctions.o
-siesta_tddft.o: sankey_change_basis.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
+siesta_tddft.o: sparse_matrices.o state_analysis.o state_init.o sys.o
+siesta_tddft.o: wavefunctions.o
 sparse_matrices.o: alloc.o class_Fstack_Pair_Geometry_SpData2D.o
 sparse_matrices.o: class_OrbitalDistribution.o class_SpData1D.o
 sparse_matrices.o: class_SpData2D.o class_Sparsity.o precision.o
@@ -1311,19 +1310,19 @@
 state_init.o: alloc.o atomlist.o class_Data2D.o class_SpData1D.o
 state_init.o: class_SpData2D.o class_Sparsity.o create_Sparsity_SC.o
 state_init.o: domain_decom.o files.o hsparse.o iodm_netcdf.o iodmhs_netcdf.o
-state_init.o: ioxv.o kpoint_grid.o ldau_specs.o m_chess.o m_energies.o m_eo.o
-state_init.o: m_gamma.o m_mixing.o m_mixing_scf.o m_mpi_utils.o m_new_dm.o
-state_init.o: m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
+state_init.o: iotdxv.o ioxv.o kpoint_grid.o ldau_specs.o m_chess.o m_energies.o
+state_init.o: m_eo.o m_gamma.o m_mixing.o m_mixing_scf.o m_mpi_utils.o
+state_init.o: m_new_dm.o m_os.o m_pivot_methods.o m_rmaxh.o m_sparse.o
 state_init.o: m_sparsity_handling.o m_spin.o m_steps.o m_supercell.o
 state_init.o: m_test_io.o m_ts_charge.o m_ts_electype.o m_ts_global_vars.o
 state_init.o: m_ts_io.o m_ts_options.o m_ts_sparse.o m_ts_tri_init.o
 state_init.o: normalize_dm.o overlap.o parallel.o proximity_check.o
 state_init.o: siesta_cml.o siesta_geom.o siesta_options.o sparse_matrices.o
-state_init.o: sys.o units.o write_subs.o zmatrix.o iotdxv.o
-struct_init.o: alloc.o atmfuncs.o atomlist.o files.o ioxv.o m_exp_coord.o
-struct_init.o: m_iostruct.o m_mpi_utils.o m_steps.o parallel.o periodic_table.o
-struct_init.o: siesta_cml.o siesta_geom.o siesta_master.o siesta_options.o
-struct_init.o: units.o zmatrix.o iotdxv.o
+state_init.o: sys.o units.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
+struct_init.o: siesta_options.o units.o zmatrix.o
 sys.o: parallel.o siesta_cml.o
 timer.o: extrae_eventllist.o extrae_module.o m_timer.o parallel.o sys.o
 timer.o: timer_tree.o
@@ -1343,8 +1342,8 @@
 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: parallel.o precision.o sparse_matrices.o kpoint_grid.o
-wavefunctions.o: atomlist.o m_spin.o m_matswinvers.o
+wavefunctions.o: atomlist.o kpoint_grid.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
 write_inp_wannier.o: siesta2wannier90.o siesta_geom.o
@@ -1459,8 +1458,8 @@
 m_grdsam.o: grdsam.o
 m_hsparse.o: hsparse.o
 m_intramol_pressure.o: remove_intramol_pressure.o
+m_iotdxv.o: iotdxv.o
 m_ioxv.o: ioxv.o
-m_iotdxv.o: iotdxv.o
 m_kinefsm.o: kinefsm.o
 m_ksv.o: ksv.o
 m_ksvinit.o: ksvinit.o
@@ -1489,6 +1488,7 @@
 m_rhofft.o: rhofft.o
 m_rhoofd.o: rhoofd.o
 m_rusage.o: rusage.o
+m_sankey_change_basis.o: sankey_change_basis.o
 m_save_density_matrix.o: save_density_matrix.o
 m_scf_options.o: m_options.o
 m_scfconvergence_test.o: scfconvergence_test.o

=== modified file 'Src/cranknic_evolg.F90'
--- Src/cranknic_evolg.F90	2017-10-03 13:27:08 +0000
+++ Src/cranknic_evolg.F90	2017-10-03 17:43:24 +0000
@@ -56,21 +56,22 @@
       real(dp)             :: delt
       !
       type(matrix)         :: Hauxms,Sauxms, wfaux1, wfaux2
-      character(3)         :: m_operation
-      character(5)         :: m_storage
       complex(dp)          :: cvar1, cvar2
+
+#ifdef MPI
+      character(len=5), parameter :: m_storage = 'pzdbc'
+      character(len=3), parameter :: m_operation = 'lap'
+#else
+      character(len=5), parameter :: m_storage = 'szden'
+      character(len=5), parameter :: m_operation = 'lap'
+#endif
+
       ! 
       integer              :: i, j, io, jo, ie, ispin, ind, nocc 
       !
       real(dp)             :: eigv
       logical, save        :: firsttime = .true.
-#ifdef MPI
-      m_storage='pzdbc'
-      m_operation='lap'
-#else
-      m_storage='szden'
-      m_operation='lap'
-#endif
+
 #ifdef DEBUG
       call write_debug( '    PRE cn_evolg' )
 #endif
@@ -165,38 +166,46 @@
 !*************************************************************************   
   
  implicit none 
- !     
+
  integer               :: no
  real(kind=dp)         :: deltat
  type(matrix)          :: H,S,phi
+ 
  ! Internal variables 
- type(matrix)          :: aux1,aux2,aux3,aux4
+ type(matrix)          :: LHS, RHS
  complex(kind=dp)      :: alpha
- character             :: m_storage*5, m_operation*3
- logical, save         :: inversemm_linear = .true.
- !
+
+ complex(kind=dp), parameter :: cZERO = cmplx(0._dp, 0._dp, dp)
+ complex(kind=dp), parameter :: cONE = cmplx(1._dp, 0._dp, dp)
+ 
 #ifdef MPI
- m_storage='pzdbc'
- m_operation='lap'
+ character(len=5), parameter :: m_storage = 'pzdbc'
+ character(len=3), parameter :: m_operation = 'lap'
 #else
- m_storage='szden'
- m_operation='lap'
+ character(len=5), parameter :: m_storage = 'szden'
+ character(len=5), parameter :: m_operation = 'lap'
 #endif
- !  
- call m_allocate(aux1,no,no,m_storage)
- call m_allocate(aux2,no,no,m_storage)
- call m_allocate(aux3,phi%dim1,phi%dim2,m_storage)
+ 
+ logical, parameter :: inversemm_linear = .true.
+ 
  ! First order expansion for the evolution operator
- alpha=-0.5_dp*cmplx(0.0_dp,1.0_dp,dp)*deltat
- ! Copying S to aux1 and aux2
- call m_add(S,'n',aux1,cmplx(1.0_dp,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
- call m_add(S,'n',aux2,cmplx(1.0_dp,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
- ! Calculating S - alpha * H
- call m_add(H,'n',aux1,alpha,cmplx(1.0_dp,0.0_dp,dp),m_operation)
- ! Calculating (S - alpha * H) * phi
- call mm_multiply(aux1,'n',phi,'n',aux3,cmplx(1.0,0.0,dp),cmplx(0.0,0.0,dp),m_operation)
- ! Calculating S + alpha * H  
- call m_add(H,'n',aux2,-1.0_dp*alpha,cmplx(1.0_dp,0.0_dp,dp),m_operation)
+ alpha = -0.5_dp * cmplx(0.0_dp,1.0_dp,dp) * deltat
+
+ ! Allocate work arrays
+ call m_allocate(LHS,no,no,m_storage)
+ call m_allocate(RHS,phi%dim1,phi%dim2,m_storage)
+
+ ! Setup S - alpha * H
+ call m_add(S, 'n', LHS, cone, cZERO, m_operation)
+ call m_add(H, 'n', LHS, alpha, cONE, m_operation)
+ 
+ ! Calculate
+ !   (S - alpha * H) psi
+ call mm_multiply(LHS, 'n', phi, 'n', RHS, cONE, cZERO, m_operation)
+
+ ! Setup S + alpha * H
+ call m_add(S, 'n', LHS, cONE, cZERO, m_operation)
+ call m_add(H, 'n', LHS, -alpha, cONE, m_operation)
 
  !---------------------------------------------------------------------------------------!
  ! There are two ways to compute inverse. One is to first obtain
@@ -209,23 +218,21 @@
  ! hard-wired by the inversemm_linear flag, while keeping first one for testing
  ! etc.
  !----------------------------------------------------------------------------------------!
-
- if (inversemm_linear) then
+ if ( inversemm_linear ) then
    ! Calculating (S + alpha * H)^-1 * (S - alpha * H) * phi 
-   call inversemm(aux2,aux3) 
-   ! Compying phi_evolved back to phi
-   call m_add(aux3,'n',phi,cmplx(1.0_dp,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
+   call inversemm(LHS, RHS)
+   ! Copying phi_evolved back to phi
+   call m_add(RHS, 'n', phi, cONE, cZERO, m_operation)
  else
    ! Calculating inverse of (S + alpha * H)
-   call getinverse(aux2)
-   !(S + alpha * H)^-1 * (S - alpha * H) * phi 
-   call mm_multiply(aux2,'n',aux3,'n',phi,cmplx(1.0,0.0,dp),cmplx(0.0,0.0,dp),m_operation)
+   call getinverse(LHS)
+   ! (S + alpha * H)^-1 * (S - alpha * H) * phi 
+   call mm_multiply(LHS, 'n', RHS, 'n', phi, cONE, cZERO, m_operation)
  endif
- !
- call m_deallocate(aux1)
- call m_deallocate(aux2)
- call m_deallocate(aux3)
- !
+
+ call m_deallocate(LHS)
+ call m_deallocate(RHS)
+
  END SUBROUTINE Uphi
  !------------------------------------------------------------------------------------!
  SUBROUTINE evol1new(Hauxms, Sauxms, no, ispin,         & 
@@ -259,8 +266,13 @@
   !
   type(matrix),intent(in)         :: Hauxms, Sauxms 
   type(matrix),allocatable,save   :: Hsve(:)
-  character(5)                    :: m_storage
-  character(3)                    :: m_operation
+#ifdef MPI
+  character(len=5), parameter :: m_storage = 'pzdbc'
+  character(len=3), parameter :: m_operation = 'lap'
+#else
+  character(len=5), parameter :: m_storage = 'szden'
+  character(len=5), parameter :: m_operation = 'lap'
+#endif
   logical                         :: extrapol
   ! Internal variables ...
   integer                :: i, l
@@ -268,14 +280,7 @@
   logical, save          :: fsttim(2) = (/.true. , .true./)
   logical, save          :: frsttime = .true.
   save                   ::  deltat
-  !
-#ifdef MPI
-  m_storage='pzdbc'
-  m_operation='lap'
-#else
-  m_storage='szden'
-  m_operation='lap'
-#endif
+
   if (frsttime) then
     !nstp is the number of "substeps" in the electronic evolution
     !the evolution operator is applied in each substep although
@@ -283,10 +288,10 @@
     !a SCF Hamiltonian
     deltat=delt/0.04837769d0/dble(nstp)
     if (Node.eq.0) then
-      write(6,*) 'cn_evolg: TDED time step (fs)      = ',delt
-      if (extrapol) then
-        write(6,*) 'cn_evolg: TDED time sub-step (fs)  = ',delt/nstp
-      end if
+       write(6,'(/a,f16.6)') 'cn_evolg: TDED time step (fs)      = ',delt
+       if (extrapol) then
+          write(6,'(a,f16.6)') 'cn_evolg: TDED time sub-step (fs)  = ',delt/nstp
+       end if
     end if
     allocate(Hsve(nspin))
     do i=1, nspin

=== modified file 'Src/cranknic_evolk.F90'
--- Src/cranknic_evolk.F90	2017-08-18 18:36:36 +0000
+++ Src/cranknic_evolk.F90	2017-10-03 17:43:24 +0000
@@ -43,21 +43,20 @@
 
   TYPE(matrix)             :: Hauxms, Sauxms, wfaux1, wfaux2
 
-  CHARACTER(LEN=3)         :: m_operation
-  CHARACTER(LEN=5)         :: m_storage
+#ifdef MPI
+  character(len=5), parameter :: m_storage = 'pzdbc'
+  character(len=3), parameter :: m_operation = 'lap'
+#else
+  character(len=5), parameter :: m_storage = 'szden'
+  character(len=5), parameter :: m_operation = 'lap'
+#endif
 
   COMPLEX(dp)              :: cvar1, cvar2
 
   REAL(dp)                 :: kxij, ckxij, skxij
   INTEGER                  :: ik, ispin, i, j, io, jo, ind, juo, nocc
   LOGICAL, SAVE            :: firstime = .true.
-#ifdef MPI
-  m_storage    = 'pzdbc'
-  m_operation  = 'lap'
-#else
-  m_storage    = 'szden'
-  m_operation  = 'lap'
-#endif
+
 #ifdef DEBUG
   call write_debug( '    PRE cn_evolk' )
 #endif
@@ -143,21 +142,19 @@
  TYPE(matrix)                  :: Hauxms, Sauxms
  INTEGER                       :: i, j, l, ik, ispin
  REAL(dp)                      :: delt, deltat, rvar1
- CHARACTER(3)                  :: m_operation
- CHARACTER(5)                  :: m_storage
+#ifdef MPI
+ character(len=5), parameter :: m_storage = 'pzdbc'
+ character(len=3), parameter :: m_operation = 'lap'
+#else
+ character(len=5), parameter :: m_storage = 'szden'
+ character(len=5), parameter :: m_operation = 'lap'
+#endif
 
  COMPLEX(dp)                   :: alpha
 
  LOGICAL, SAVE                                  :: firsttime  = .true.
  LOGICAL, DIMENSION (:,:), ALLOCATABLE, SAVE    :: firstimeK
-#ifdef MPI
- m_storage   = 'pzdbc'
- m_operation = 'lap'
-#else
- m_storage   = 'szden'
- m_operation = 'lap'
-#endif
- !
+
  IF(firsttime) THEN 
    deltat = delt/0.04837769d0/dble(ntded_sub)
    IF (IOnode) THEN

=== modified file 'Src/diag.F90'
--- Src/diag.F90	2017-08-07 10:57:47 +0000
+++ Src/diag.F90	2017-10-03 17:43:24 +0000
@@ -176,6 +176,9 @@
   end subroutine diag_exit
 
   subroutine diag_descinit(N, NR, BlockSize, desc, iC)
+#ifdef MPI
+    use mpi_siesta, only: MPI_Comm_World
+#endif
     integer, intent(in) :: N, NR, BlockSize
     integer, intent(inout) :: desc(9)
     integer, intent(in), optional :: iC
@@ -184,6 +187,7 @@
     if ( present(iC) ) then
        call descinit(desc,N,N,BlockSize,BlockSize,0,0,iC,NR,ierr)
     else
+       if ( iCTXT < 0 ) iCTXT = MPI_Comm_World
        call descinit(desc,N,N,BlockSize,BlockSize,0,0,iCTXT,NR,ierr)
     end if
     if ( ierr /= 0 ) &

=== modified file 'Src/m_initwf.F90'
--- Src/m_initwf.F90	2017-08-30 14:09:10 +0000
+++ Src/m_initwf.F90	2017-10-03 17:43:24 +0000
@@ -13,8 +13,6 @@
   private
   !
   public :: initwf
-  !
-  integer :: ictxt
 
   CONTAINS
   !
@@ -127,7 +125,6 @@
         if(ParallelOverk) then
           call die ('initwf: TDDFT is not parallelized over k-points.')
         end if
-        ictxt = MPI_Comm_World
 #endif
 !       Read spin-spiral wavevector (if defined)
         call readsp( qspiral, spiral )
@@ -303,8 +300,8 @@
   subroutine diaggiwf(nspin,nuo,maxuo,maxnh, maxo,Haux,Saux,psi,           &
                       nuotot,occup)
 #ifdef MPI
-      use parallel, only : BlockSize,Node
-      use m_diag, only: diag_descinit
+    use parallel, only : BlockSize,Node
+    use m_diag, only: diag_descinit
 #endif
       !
       implicit none
@@ -322,7 +319,7 @@
 #endif
       !
 #ifdef MPI
-      call diag_descinit(nuotot,nuotot,BlockSize,desch,ictxt)
+      call diag_descinit(nuotot,nuotot,BlockSize,desch)
 #endif
       !
       indwf=0
@@ -338,7 +335,7 @@
               Haux(jo,io)=Haux(jo,io)+H(ind,ispin)
             end do
           end do
-          call rdiag(Haux,Saux,nuotot,nuo,nuotot,eo,psi(1,1,ispin),nuotot,1,ierror)
+          call rdiag(Haux,Saux,nuotot,nuo,nuotot,eo,psi(1,1,ispin),nuotot,1,ierror,BlockSize)
           if (ierror .eq. 0) then
             exit
           else if ((ierror .ne. -1) .or. (ie .eq. 10)) then
@@ -388,7 +385,7 @@
 #endif
       !
 #ifdef MPI
-      call diag_descinit(nuotot,nuotot,BlockSize,desch,ictxt)
+      call diag_descinit(nuotot,nuotot,BlockSize,desch)
 #endif
       !
       indwf=0
@@ -426,10 +423,7 @@
               indwf=indwf+1
               do j=1,nuotot
 #ifdef MPI
-                do iuo=1,2
-                call pdelget('a',' ',varaux(iuo),psi(iuo,:,:),j,ie,desch)
-                enddo
-                element=cmplx(varaux(1),varaux(2),dp)
+                call pzelget('a',' ',element,psi(:,:,:),j,ie,desch)
 #else
                 element=cmplx(psi(1,j,ie),psi(2,j,ie),dp)
 #endif

=== modified file 'Src/m_inversemm.F90'
--- Src/m_inversemm.F90	2017-10-03 13:27:08 +0000
+++ Src/m_inversemm.F90	2017-10-03 17:43:24 +0000
@@ -20,7 +20,7 @@
 IMPLICIT NONE
 !**************** INPUT ***********************************!
 !
-TYPE(matrix), INTENT(IN) :: C  
+TYPE(matrix), INTENT(INOUT) :: C  
 ! C: Matrix to be inverted before multiplication with D.
 !
 !**************** INOUT ***********************************!
@@ -42,15 +42,13 @@
  ALLOCATE(ipiv(C%iaux1(3)+C%iaux1(5)))
  CALL pzgesv(C%dim1,D%dim2,C%zval,1,1,C%iaux1,ipiv,D%zval,1,1,D%iaux1,info)
  IF (info .NE. 0) CALL die('ERROR: error in pzgesv')
- DEALLOCATE(ipiv)      
- !
 #else
- !
  ALLOCATE(ipiv(C%dim1))
  CALL zgesv(C%dim1,D%dim2,C%zval,C%dim1,ipiv,D%zval,D%dim1,info)
  IF (info .NE. 0) CALL die('ERROR: error in zgsesv')
+#endif
+
  DEALLOCATE(ipiv)
-#endif
 !
 END SUBROUTINE inversemm 
 !

=== modified file 'version.info'
--- version.info	2017-10-03 13:27:08 +0000
+++ version.info	2017-10-03 17:43:24 +0000
@@ -1,1 +1,1 @@
-trunk-611-tddft02-40
+trunk-611-tddft02-40--np-1