1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8 MODULE m_state_init 9 10 private 11 public :: state_init 12 13 CONTAINS 14 15 subroutine state_init( istep ) 16 use Kpoint_grid, only: setup_Kpoint_grid, nkpnt, gamma_scf 17 use m_os, only: file_exist 18 use m_new_dm, only: new_dm 19 use m_proximity_check, only: proximity_check 20 use siesta_options 21 use units, only: Ang 22 use sparse_matrices, only: maxnh, numh, listh, listhptr 23 use sparse_matrices, only: Dold, Dscf, DM_2D 24 use sparse_matrices, only: Eold, Escf, EDM_2D 25 use sparse_matrices, only: Hold, H, H_2D 26 use sparse_matrices, only: xijo, xij_2D 27 use sparse_matrices, only: S, S_1D 28 29 use sparse_matrices, only: H_kin_1D, H_vkb_1D 30 use sparse_matrices, only: H_dftu_2D, H_so_2D 31 32 use sparse_matrices, only: sparse_pattern 33 use sparse_matrices, only: block_dist, single_dist 34 use sparse_matrices, only: DM_history 35 36 37 use create_Sparsity_SC, only: crtSparsity_SC 38 use m_sparsity_handling, only: SpOrb_to_SpAtom 39 use m_sparsity_handling, only: Sp_to_Spglobal 40 use m_pivot_methods, only: sp2graphviz 41 42 use siesta_geom 43 use atomlist, only: iphorb, iphkb, indxua, iaorb, 44 & rmaxo, rmaxkb, rmaxv, rmaxdftu, 45 & lastkb, lasto, superc, indxuo, 46 & no_u, no_s, no_l, iza 47 use alloc, only: re_alloc, de_alloc, alloc_report 48 use m_hsparse, only: hsparse 49 use m_overlap, only: overlap 50 use m_supercell, only: exact_sc_ag 51 use siesta_cml, only: cml_p, cmlStartStep, mainXML 52 use siesta_cml, only: cmlStartPropertyList 53 use siesta_cml, only: cmlEndPropertyList 54 use siesta_cml, only: cmlAddProperty 55 use zmatrix, only: lUseZmatrix, write_zmatrix 56 use m_energies, only: Emad 57 use write_subs 58 use m_ioxv, only: ioxv 59 use m_steps 60 use parallel, only: IOnode, node, nodes 61 use m_spin, only: spin 62 use m_rmaxh 63 use m_mixing, only: mixers_history_init 64 use m_mixing_scf, only: scf_mixs, scf_mix 65 66 use m_normalize_dm, only: normalize_dm 67 68 use m_eo 69 use files, only: slabel 70 use m_mpi_utils, only: globalize_or 71 use m_mpi_utils, only: globalize_sum 72 use domain_decom, only: domainDecom, use_dd, use_dd_perm 73 use dftu_specs, only: switch_dftu, dftu_init 74 use fdf, only: fdf_get 75 use sys, only: message, die, bye 76 use m_sparse, only : xij_offset 77 78 use m_ts_kpoints, only: setup_ts_kpoint_grid 79 use ts_dq_m, only : TS_DQ_METHOD, TS_DQ_METHOD_FERMI 80 use m_ts_options, only : BTD_method 81 use m_ts_options, only : TS_Analyze 82 use m_ts_options, only : N_Elec, Elecs, IsVolt 83 use m_ts_electype 84 use m_ts_global_vars, only: TSrun, TSmode, onlyS 85 use m_ts_io, only : fname_TSHS, ts_write_tshs 86 use m_ts_sparse, only : ts_sparse_init 87 use m_ts_tri_init, only : ts_tri_init, ts_tri_analyze 88 use files, only: slabel, label_length 89#ifdef CDF 90 use iodm_netcdf, only: setup_dm_netcdf_file 91 use iodmhs_netcdf, only: setup_dmhs_netcdf_file 92#endif 93#ifdef NCDF_4 94 use dictionary 95 use m_ncdf_siesta, only : cdf_init_file, cdf_save_settings 96 use m_ncdf_siesta, only : cdf_save_state, cdf_save_basis 97#ifdef MPI 98 use mpi_siesta 99#endif 100#endif 101 102 use class_Sparsity 103 use class_dSpData1D 104 use class_dSpData2D 105 use class_dData2D 106 use class_Pair_Geometry_dSpData2D 107 use class_Fstack_Pair_Geometry_dSpData2D 108 109#ifdef TEST_IO 110 use m_test_io 111#endif 112#ifdef SIESTA__FLOOK 113 use siesta_dicts, only : dict_repopulate_MD 114 use siesta_dicts, only : dict_repopulate_sparse_matrices 115#endif 116 117 use m_handle_sparse, only: correct_supercell_SpD 118 use m_restruct_SpData2D, only: restruct_dSpData2D 119 120 implicit none 121 122 integer :: istep, nnz 123 real(dp) :: veclen ! Length of a unit-cell vector 124 real(dp) :: rmax 125 logical :: cell_can_change 126 integer :: i, ix, iadispl, ixdispl 127 logical :: auxchanged ! Auxiliary supercell changed? 128 logical :: folding, folding1 129 logical :: diag_folding, diag_folding1 130 logical :: foundxv ! dummy for call to ioxv 131 132 external :: madelung, timer 133 real(dp), external :: volcel 134 135 integer :: ts_kscell_file(3,3) = 0 136 real(dp) :: ts_kdispl_file(3) = 0.0 137 logical :: ts_Gamma_file = .true. 138 character(len=label_length+6) :: fname 139 real(dp) :: dummyef=0.0, dummyqtot=0.0 140 type(Sparsity) :: g_Sp 141#ifdef NCDF_4 142 type(dictionary_t) :: d_sav 143#ifdef MPI 144 integer :: MPIerror 145#endif 146#endif 147 148 character(len=256) :: oname 149 150 type(dData2D) :: tmp_2D 151 152! Variables required to correct the DM in the history 153 type(dSpData2D), pointer :: tmp_Sp2D 154 type(Pair_Geometry_dSpData2D), pointer :: pair 155 156 real(dp) :: dummy_qspin(8) 157 158!------------------------------------------------------------------- BEGIN 159 call timer( 'IterGeom', 1 ) 160 161 call timer( 'state_init', 1 ) 162 163 istp = istp + 1 164 165 if (IOnode) then 166 write(6,'(/,t22,a)') repeat('=',36) 167 select case (idyn) 168 case (0) 169 if (nmove == 0) then 170 write(6,'(t25,a)') 'Single-point calculation' 171 if (cml_p) call cmlStartStep(mainXML, type='Single-Point', 172 $ index=istp) 173 else 174 if (broyden_optim) then 175 write(6,'(t25,a,i6)') 'Begin Broyden opt. move = ', 176 $ istep 177 else if (fire_optim) then 178 write(6,'(t25,a,i6)') 'Begin FIRE opt. move = ', 179 $ istep 180 else 181 write(6,'(t25,a,i6)') 'Begin CG opt. move = ', 182 $ istep 183 end if 184 if (cml_p) call cmlStartStep(mainXML, type='Geom. Optim', 185 $ index=istp) 186 endif 187 188! Print Z-matrix coordinates 189 if (lUseZmatrix) then 190 call write_Zmatrix() 191 endif 192 case (1, 3) 193 if (iquench > 0 ) then 194 write(6,'(t25,a,i6)') 'Begin MD quenched step = ', 195 $ istep 196 if (cml_p) call cmlStartStep(mainXML, type='MD-quenched', 197 $ index=istep) 198 else 199 write(6,'(t25,a,i6)') 'Begin MD step = ', 200 $ istep 201 if (cml_p) call cmlStartStep(mainXML, type='MD', 202 $ index=istep) 203 endif 204 case (2,4,5) 205 write(6,'(t25,a,i6)') 'Begin MD step = ', istep 206 if (cml_p) call cmlStartStep(mainXML, type='MD', index=istep) 207 case (6) 208 write(6,'(t25,a,i6)') 'Begin FC step = ',istep 209 if (cml_p) call cmlStartStep(mainXML, type='FC', index=istep) 210 211 if (istep .eq. 0) then 212 write(6,'(t25,a)') 'Undisplaced coordinates' 213 else 214 iadispl = (istep-mod(istep-1,6))/6+ia1 215 ix = mod(istep-1,6)+1 216 ixdispl = (ix - mod(ix-1,2) +1)/2 217 write(6,'(t26,a,i0,/,t26,a,i1,a,f10.6,a)') 'displace atom ', 218 & iadispl,'in direction ',ixdispl,' by', dx/Ang,' Ang' 219 endif 220 221 case (8) 222 write(6,'(t25,a,i6)') 'Begin Server step = ',istep 223 if (cml_p) call cmlStartStep(mainXML, type='FS', index=istep) 224 225 case (9) 226 if ( istep == 0 ) then 227 write(6,'(t25,a,i7)')'Explicit coord. initialization' 228 else 229 write(6,'(t25,a,i7)')'Explicit coord. step =',istep 230 end if 231 if (cml_p) call cmlStartStep(mainXML, type='ECS', index=istep) 232 233 case (10) 234 write(6,'(t25,a,i7)')'LUA coord. step =',istep 235 if (cml_p) call cmlStartStep(mainXML, type='LUA', index=istep) 236 237 end select 238 239 write(6,'(t22,a)') repeat('=',36) 240 241! Print atomic coordinates 242 call outcoor( ucell, xa, na_u, ' ', writec ) 243 244! Save structural information in crystallographic format 245! (in file SystemLabel.STRUCT_OUT), 246! canonical Zmatrix (if applicable), and CML record 247 248 call siesta_write_positions(moved=.false.) 249 250 endif ! IONode 251 252 ! Write the XV file for single-point calculations, so that 253 ! it is there at the end for those users who rely on it 254 call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, 255 & foundxv) 256 257! Actualize things if variable cell 258! These checks are made to ensure that the k-points 259! and Madelung terms are corrected in case the cell is changed. 260 auxchanged = .false. 261 cell_can_change = ( varcel .or. 262 & (idyn .eq. 8) ! Force/stress evaluation 263 & ) 264 if (change_kgrid_in_md) then 265 cell_can_change = cell_can_change .or. 266 & (idyn .eq. 3) .or. ! Parrinello-Rahman 267 & (idyn .eq. 4) .or. ! Nose-Parrinello-Rahman 268 & (idyn .eq. 5) ! Anneal 269 endif 270 271 if ( cell_can_change .and. istep /= inicoor ) then 272 273! Madelung correction for charged systems 274 if (charnet .ne. 0.0_dp) then 275 call madelung(ucell, shape, charnet, Emad) 276 end if 277 278 if ( .not. gamma_scf ) then 279 280! Will print k-points also 281 call setup_Kpoint_grid( ucell ) 282 283 call setup_ts_kpoint_grid( ucell ) 284 285 call re_alloc( eo, 1, no_u, 1, spin%spinor, 1, nkpnt, 'eo', 286 & 'state_init') 287 call re_alloc( qo, 1, no_u, 1, spin%spinor, 1, nkpnt, 'qo', 288 & 'state_init' ) 289 end if 290 291 end if 292! End variable cell actualization 293 294! Always (in case of auxiliary cell) check the new sparse pattern 295 if ( use_aux_cell ) then 296! Determine the tightest auxiliary supercell using 297! also the atomic positions 298 call exact_sc_ag(negl,ucell,na_u,isa,xa,nsc) 299 300 mscell = 0.0_dp 301 do i = 1, 3 302 mscell(i,i) = nsc(i) 303 if (nsc(i) /= nsc_old(i)) auxchanged = .true. 304 end do 305 306 end if 307 308! Auxiliary supercell 309! Do not move from here, as the coordinates might have changed 310! even if not the unit cell 311 call superc(ucell, scell, nsc) 312#ifdef SIESTA__FLOOK 313 call dict_repopulate_MD() 314#endif 315 316! Print unit cell and compute cell volume 317! Possible BUG: 318! Note that this volume is later used in write_subs and the md output 319! routines, even if the cell later changes. 320 if (IOnode) call outcell(ucell) 321 volume_of_some_cell = volcel(ucell) 322 323! Use largest possible range in program, except hsparse... 324! 2 * rmaxv: Vna overlap 325! rmaxo + rmaxkb: Non-local KB action 326! 2 * (rmaxo + rmaxdftu): Interaction through DFTU projector 327! 2.0_dp * (rmaxo+rmaxkb) : Orbital interaction through KB projectors 328 rmax = max( 2._dp*rmaxv, 2._dp*(rmaxo+rmaxdftu), rmaxo+rmaxkb) 329 330 if ( .not. negl ) then 331 rmax = max(rmax, 2.0_dp * (rmaxo+rmaxkb) ) 332 endif 333 334! Check if any two atoms are unreasonably close 335 call proximity_check(rmax) 336 337 ! Clear history of mixing parameters 338 call mixers_history_init( scf_mixs ) 339 scf_mix => scf_mixs(1) 340 341 ! Ensure sparsity pattern is empty 342 call delete(sparse_pattern) 343 ! sadly deleting the sparse pattern does not necessarily 344 ! mean that the arrays are de-associated. 345 ! Remember that the reference counter could (in MD) 346 ! be higher than 1, hence we need to create "fake" 347 ! containers and let the new<class> delete the old 348 ! sparsity pattern 349 nullify(numh,listhptr,listh) 350 allocate(numh(no_l),listhptr(no_l)) 351 ! We do not need to allocate listh 352 ! that will be allocated in hsparse 353 354! List of nonzero Hamiltonian matrix elements 355! and, if applicable, vectors between orbital centers 356 357! Listh and xijo are allocated inside hsparse 358! Note: We always generate xijo now, for COOP and other 359! analyses. 360 call delete(xij_2D) ! as xijo will be reallocated 361 nullify(xijo) 362 call hsparse( negl, scell, nsc, na_s, isa, xa, lasto, 363 & lastkb, iphorb, iphKB, maxnh, .not. use_aux_cell, 364 $ set_xijo=.true., folding=folding1, 365 $ diagonal_folding=diag_folding1, 366 $ debug_folding=fdf_get('debug-folding',.false.)) 367! 368 call globalize_or(diag_folding1,diag_folding) 369 call globalize_or(folding1,folding) 370 if (diag_folding .and. .not. use_aux_cell ) then 371 call message("WARNING","Gamma-point calculation " // 372 $ "with interaction between periodic images") 373 call message("WARNING", 374 $ "Some features might not work optimally:") 375 call message("WARNING", 376 $ "e.g. DM initialization from atomic data") 377 if (harrisfun) call die("Harris functional run needs " // 378 $ "'force-aux-cell T'") 379 380 else if (folding) then 381 if ( .not. use_aux_cell ) then 382 call message("INFO","Gamma-point calculation " // 383 $ "with multiply-connected orbital pairs") 384 call message("INFO", 385 $ "Folding of H and S implicitly performed") 386 call check_cohp() 387 else 388 write(6,"(a,/,a)") "Non Gamma-point calculation " // 389 $ "with multiply-connected orbital pairs " // 390 $ "in auxiliary supercell.", 391 $ "Possible internal error. " // 392 $ "Use 'debug-folding T' to debug." 393 call die("Inadequate auxiliary supercell") 394 endif 395 endif 396! 397 call globalize_sum(maxnh,nnz) 398 if (cml_p) then 399 call cmlStartPropertyList(mainXML,title='Orbital info') 400 call cmlAddProperty(xf=mainXML, value=no_u, 401 $ title='Number of orbitals in unit cell', 402 $ dictref='siesta:no_u', units="cmlUnits:countable") 403 call cmlAddProperty(xf=mainXML, value=nnz, 404 $ title='Number of non-zeros', 405 $ dictref='siesta:nnz', units="cmlUnits:countable") 406 call cmlEndPropertyList(mainXML) 407 endif 408 ! 409 ! 410 ! If using domain decomposition, redistribute orbitals 411 ! for this geometry, based on the hsparse info. 412 ! The first time round, the initial distribution is a 413 ! simple block one (given by preSetOrbitLimits). 414 ! 415 ! Any DM, etc, read from file will be redistributed according 416 ! to the new pattern. 417 ! Inherited DMs from a previous geometry cannot be used if the 418 ! orbital distribution changes. For now, we avoid changing the 419 ! distribution (the variable use_dd_perm is .true. if domain 420 ! decomposition is in effect). Names should be changed... 421 422 if (use_dd .and. (.not. use_dd_perm)) then 423 call domainDecom( no_u, no_l, maxnh ) ! maxnh intent(in) here 424 maxnh = sum(numh(1:no_l)) 425 ! We still need to re-create Julian Gale's 426 ! indexing for O(N) in parallel. 427 print "(a5,i3,a20,3i8)", 428 $ "Node: ", Node, "no_u, no_l, maxnh: ", no_u, no_l, maxnh 429 call setup_ordern_indexes(no_l, no_u, Nodes) 430 endif 431 432 ! I would like to skip this alloc/move/dealloc/attach 433 ! by allowing sparsity to have pointer targets. 434 ! However, this poses a problem with intel compilers, 435 ! as it apparently errors out when de-allocating a target pointer 436 write(oname,"(a,i0)") "sparsity for geom step ", istep 437 call newSparsity(sparse_pattern,no_l,no_u,maxnh, 438 & numh,listhptr,listh, name = oname) 439 deallocate(numh,listhptr,listh) 440 call attach(sparse_pattern, 441 & n_col = numh, list_ptr = listhptr, list_col = listh ) 442 443 ! In case the user requests to create the connectivity graph 444 if ( write_GRAPHVIZ > 0 ) then 445 ! first create the unit-cell sparsity pattern 446 call crtSparsity_SC(sparse_pattern, g_Sp, UC=.true.) 447 ! next move to global sparsity pattern 448 call Sp_to_Spglobal(block_dist, g_Sp, g_Sp) 449 if ( IONode ) then 450 if ( write_GRAPHVIZ /= 2 ) 451 & call sp2graphviz(trim(slabel)//'.ORB.gv', g_Sp) 452 ! Convert to atomic 453 if ( write_GRAPHVIZ /= 1 ) then 454 call SpOrb_to_SpAtom(single_dist,g_Sp,na_u,lasto,g_Sp) 455 call sp2graphviz(trim(slabel)//'.ATOM.gv', g_Sp) 456 end if 457 end if 458 call delete(g_Sp) 459 end if 460 461 462 ! Copy over xijo array (we can first do it here... :( ) 463 call newdData2D(tmp_2D,xijo,'xijo') 464 deallocate(xijo) 465 write(oname,"(a,i0)") "xijo at geom step ", istep 466 call newdSpData2D(sparse_pattern,tmp_2D,block_dist,xij_2D, 467 & name=oname) 468 call delete(tmp_2D) ! decrement container... 469 xijo => val(xij_2D) 470 471! Convert the xijo array into a super cell offset array 472! isc_off (located in siesta_geom). 473! This array is primarily used in TranSiesta but may easily 474! be used elsewhere to figure out orbital/atomic placements 475! in the sparsity pattern. 476 if ( .not. use_aux_cell ) then 477 ! Here we create the super-cell offsets 478 call re_alloc(isc_off,1,3,1,1) 479 isc_off(:,:) = 0 480 else 481 call xij_offset(ucell,nsc, na_u,xa,lasto, 482 & xij_2D, isc_off, Bcast=.true.) 483 end if 484 485 486 ! When the user requests to only do an analyzation, we can call 487 ! appropriate routines and quit 488 if ( TS_Analyze ) then 489 490 ! Force the creation of the full sparsity pattern 491 call ts_sparse_init(slabel,IsVolt, N_Elec, Elecs, 492 & ucell, nsc, na_u, xa, lasto, block_dist, sparse_pattern, 493 & .not. use_aux_cell, isc_off) 494 495 ! create the tri-diagonal matrix 496 call ts_tri_analyze( block_dist, sparse_pattern , N_Elec, 497 & Elecs, ucell, na_u, lasto, nsc, isc_off, 498 & BTD_method ) 499 500 ! Print-out timers 501 call timer('TS-analyze',3) 502 503 ! Bye also waits for all processors 504 call bye('transiesta analyzation performed') 505 506 end if 507 508 509 write(oname,"(a,i0)") "EDM at geom step ", istep 510 call newdSpData2D(sparse_pattern,spin%EDM,block_dist,EDM_2D, 511 & name=oname) 512 !if (ionode) call print_type(EDM_2D) 513 Escf => val(EDM_2D) 514 515 call re_alloc(Dold,1,maxnh,1,spin%DM,name='Dold', 516 . routine='sparseMat',copy=.false.,shrink=.true.) 517 call re_alloc(Hold,1,maxnh,1,spin%H,name='Hold', 518 . routine='sparseMat',copy=.false.,shrink=.true.) 519 if ( converge_EDM ) then 520 call re_alloc(Eold,1,maxnh,1,spin%EDM,name='Eold', 521 . routine='sparseMat',copy=.false.,shrink=.true.) 522 end if 523 524! Allocate/reallocate storage associated with Hamiltonian/Overlap matrix 525 write(oname,"(a,i0)") "H at geom step ", istep 526 call newdSpData2D(sparse_pattern,spin%H,block_dist,H_2D, 527 & name=oname) 528 !if (ionode) call print_type(H_2D) 529 H => val(H_2D) 530 531 write(oname,"(a,i0)") "H_vkb at geom step ", istep 532 call newdSpData1D(sparse_pattern,block_dist,H_vkb_1D,name=oname) 533 !if (ionode) call print_type(H_vkb_1D) 534 535 write(oname,"(a,i0)") "H_kin at geom step ", istep 536 call newdSpData1D(sparse_pattern,block_dist,H_kin_1D,name=oname) 537 !if (ionode) call print_type(H_kin_1D) 538 539 if ( switch_dftu ) then 540 write(oname,"(a,i0)") "H_dftu at geom step ", istep 541 call newdSpData2D(sparse_pattern,spin%spinor, 542 & block_dist,H_dftu_2D,name=oname) 543 544 ! Initialize to 0, LDA+U may re-calculate 545 ! this matrix sporadically doing the SCF. 546 ! Hence initialization MUST be performed upon 547 ! re-allocation. 548 call init_val(H_dftu_2D) 549 if ( inicoor /= istep ) then 550 ! Force initialization of the LDA+U 551 ! when changing geometry 552 ! For the first geometry this is controlled 553 ! by the user via an fdf-key 554 dftu_init = .true. 555 end if 556 end if 557 558 if ( spin%SO ) then 559 write(oname,"(a,i0)") "H_so at geom step ", istep 560 call newdSpData2D(sparse_pattern,spin%H - 2, 561 & block_dist,H_so_2D,name=oname) 562 end if 563 564 write(oname,"(a,i0)") "S at geom step ", istep 565 call newdSpData1D(sparse_pattern,block_dist,S_1D,name=oname) 566 if (ionode) call print_type(S_1D) 567 S => val(S_1D) 568 569!> Before proceeding we need to "fix" a few things before we can 570!> successfully read the new DM. 571!> 1. Cases where the atomic displacements yields a new sparse 572!> pattern it is vital to remove the elements that does not exist. 573!> Such cases are often encountered because atoms move in/out of 574!> neighbouring atoms orbital ranges. Everytime two orbitals meet, 575!> new sparse elements are added, and everytime two orbitals flee 576!> sparse elements are removed. 577!> 2. If the supercell has changed size it is necessary to 578!> remove/translate the old/new supercells such that they are 579!> conforming with the new sparse pattern. 580!> For details regarding the sparsity pattern, see sparse_matrices.F90 581 do i = 1, n_items(DM_history) 582 pair => get_pointer(DM_history,i) 583 call secondp(pair,tmp_Sp2D) 584 if ( auxchanged ) then 585! Transfer the old mixing matrix to the new supercell pattern 586 call correct_supercell_SpD(nsc_old, tmp_Sp2D, nsc) 587 end if 588! Only retain the new orbital interactions 589 call restruct_dSpData2D(tmp_Sp2D, sparse_pattern, tmp_Sp2D, 590 & show_warning = .false.) 591 end do 592 593 594! Find overlap matrix 595 call overlap( na_u, na_s, no_s, scell, xa, indxua, rmaxo, maxnh, 596 & lasto, iphorb, isa, numh, listhptr, listh, S ) 597 598#ifdef NCDF_4 599! At this point the sparsity pattern, overlap matrix and some 600! other details. 601! Before continuing we will create the CDF file output 602! 603! Initialize the NC file 604 if ( write_cdf ) then 605 if ( inicoor == istep ) then 606 call cdf_init_file(trim(slabel)//'.nc', is_FC=(idyn==6)) 607#ifdef MPI 608 call MPI_Barrier(MPI_Comm_World,MPIerror) 609#endif 610 611! Save the basis set (only once) 612 call cdf_save_basis(trim(slabel)//'.nc') 613#ifdef MPI 614 call MPI_Barrier(MPI_Comm_World,MPIerror) 615#endif 616 617 d_sav = d_sav//('xa'.kv.1)//('cell'.kv.1) 618 d_sav = d_sav//('sp'.kv.1)//('S'.kv.1)//('xij'.kv.1) 619 d_sav = d_sav//('isc_off'.kv.1)//('nsc'.kv.1) 620 call cdf_save_state(trim(slabel)//'.nc',d_sav) 621 call delete(d_sav) 622 623 end if 624 end if 625#endif 626 627! 628! Here we could also read a Hamiltonian, either to proceed to 629! the analysis section (with nscf=0) or to start a mix-H scf cycle. 630! 631! Initialize density matrix 632! The resizing of Dscf is done inside new_dm 633 call new_DM(auxchanged, DM_history, DM_2D, EDM_2D) 634 Dscf => val(DM_2D) 635 Escf => val(EDM_2D) 636 if (spin%H > 1) call print_spin(dummy_qspin) 637 638! The number of old supercells is used in new_DM. Hence we may first 639! update the oldnsc here. 640 if ( auxchanged ) nsc_old(:) = nsc(:) 641 642 ! Initialize energy-density matrix to zero for first call to overfsm 643 ! Only part of Escf is updated in TS, so if it is put as zero here 644 ! a continuation run gives bad forces. 645 if ( .not. TSrun ) then 646 call normalize_DM( first= .true. ) 647 Escf(:,:) = 0.0_dp 648 end if 649 650#ifdef TEST_IO 651 ! We test the io-performance here 652 call time_io(spin%H,nsc,H_2D) 653#endif 654 655#ifdef SIESTA__FLOOK 656 call dict_repopulate_sparse_matrices() 657#endif 658 659! Write out the ORB_INDX file (if requested) 660 if ( save_ORB_INDX ) then 661 if ( IOnode .and. istep == inicoor ) then 662 call write_orb_indx( na_u, na_s, no_u, no_s, isa, xa, 663 . iaorb, iphorb, indxuo, nsc, ucell ) 664 end if 665 end if 666 667! If onlyS, Save overlap matrix and exit 668 if ( onlyS ) then 669 fname = fname_TSHS(slabel, onlyS = .true. ) 670 ! We include H as S, well-knowing that we only write one of 671 ! them, there is no need to allocate space for no reason! 672 call ts_write_tshs(fname, 673 & .true., .not. use_aux_cell, ts_Gamma_file, 674 & ucell, nsc, isc_off, na_u, no_s, spin%H, 675 & ts_kscell_file, ts_kdispl_file, 676 & xa, lasto, 677 & H_2D, S_1D, indxuo, 678 & dummyEf, dummyQtot, Temp,0,0) 679 call bye( 'Save overlap matrix and exit' ) ! Exit siesta 680 681 endif 682 683#ifdef CDF 684 if (writedm_cdf) then 685 call setup_dm_netcdf_file( maxnh, no_l, spin%H, 686 & no_s, indxuo, 687 & numh, listhptr, listh) 688 endif 689 if (writedm_cdf_history) then 690 call setup_dm_netcdf_file( maxnh, no_l, spin%H, 691 & no_s, indxuo, 692 & numh, listhptr, listh, 693 & istep) 694 endif 695 if (writedmhs_cdf) then 696 call setup_dmhs_netcdf_file( maxnh, no_l, spin%H, 697 & no_s, indxuo, 698 & numh, listhptr, listh, 699 & s) 700 endif 701 if(writedmhs_cdf_history) then 702 call setup_dmhs_netcdf_file( maxnh, no_l, spin%H, 703 & no_s, indxuo, 704 & numh, listhptr, listh, 705 & s, 706 & istep) 707 endif 708#endif 709 call timer( 'state_init', 2 ) 710 711 END subroutine state_init 712 713 subroutine check_cohp() 714 use siesta_options, only: write_coop 715 use sys, only: message 716 717 if (write_coop) then 718 call message("WARNING","There are multiply-connected "// 719 $ "orbitals.") 720 call message("WARNING","Your COOP/COHP analysis might " // 721 $ "be affected by folding.") 722 call message("WARNING",'Use "force-aux-cell T "' // 723 $ 'or k-point sampling') 724 endif 725 end subroutine check_cohp 726 727 END module m_state_init 728