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_siesta_init 9 private 10 public :: siesta_init 11 12 CONTAINS 13 14 subroutine siesta_init() 15 use Kpoint_grid, only: setup_Kpoint_grid, gamma_scf, nkpnt 16 USE Kpoint_pdos, only: gamma_pdos 17 use Band, only: gamma_bands, setup_bands 18 use m_ksvinit, only: gamma_polarization, 19 & estimate_pol_kpoints 20 use m_struct_init, only: struct_init 21 use units, only: Kelvin 22 use siesta_options 23 use timer_options, only: use_tree_timer, use_parallel_timer 24 use sparse_matrices 25 use class_Fstack_Pair_Geometry_dSpData2D, only: new, max_size 26 use siesta_geom 27 use atomlist, only: no_u, rmaxkb, amass, lasto, qtot, 28 & qtots, 29 & iza, rmaxo, zvaltot, superc, 30 & initatomlists, no_l, rmaxdftu 31 use fdf 32 use sys, only: die, bye 33 use molecularmechanics, only: inittwobody 34 use metaforce, only: initmeta 35 use m_mpi_utils, only: broadcast 36 use alloc, only: re_alloc, alloc_report 37 use parallelsubs, only: getnodeorbs 38 use m_iostruct, only: write_struct, read_struct 39 use zmatrix, only: lUseZmatrix 40 use zmatrix, only: write_canonical_ucell_and_Zmatrix 41 use m_supercell, only: exact_sc_ag 42 use files, only: slabel 43 use siesta_cmlsubs, only: siesta_cml_init 44 use m_timestamp, only: timestamp 45 use m_wallclock, only: wallclock 46 use parallel, only: Node, Nodes, PEXSINodes, IOnode 47 use parallel, only: SIESTA_worker 48 use parallel, only: SIESTA_Group, SIESTA_Comm 49 use m_energies 50 use m_steps 51 use m_spin, only: init_spin, spin 52 use m_spin, only: init_spiral 53 use m_rmaxh 54 use m_forces 55 use m_eo 56 use m_fixed, only: init_fixed, print_fixed 57 use m_ioxv, only: xv_file_read 58 use m_projected_DOS, only: init_projected_DOS 59 use writewave, only: gamma_wavefunctions, 60 & setup_wf_kpoints 61 use m_new_dm, only: get_allowed_history_depth 62 use m_object_debug, only: set_object_debug_level 63#ifdef DEBUG_XC 64 USE siestaXC, only: setDebugOutputUnit 65#endif /* DEBUG_XC */ 66 USE m_timer, only: timer_report ! Sets options for report of CPU times 67 use m_check_walltime 68#ifdef MPI 69 use mpi_siesta 70#endif 71 use m_cite, only: add_citation, init_citation 72 73 use m_ts_init, only : ts_init 74 75 use m_diag_option, only: read_diag, print_diag 76 use m_diag_option, only: diag_recognizes_neigwanted 77 78#ifdef SIESTA__FLOOK 79 use siesta_dicts, only : dict_populate 80 use flook_siesta, only : slua_init, slua_call, LUA_INITIALIZE 81#endif 82 83#ifdef NCDF_4 84 use netcdf_ncdf 85#endif 86#ifdef _OPENMP 87 use omp_lib, only : omp_get_num_threads 88 use omp_lib, only : omp_get_schedule, omp_set_schedule 89#if _OPENMP >= 201307 90 use omp_lib, only : omp_get_proc_bind 91 use omp_lib, only : OMP_PROC_BIND_FALSE, OMP_PROC_BIND_TRUE 92 use omp_lib, only : OMP_PROC_BIND_MASTER 93 use omp_lib, only : OMP_PROC_BIND_CLOSE, OMP_PROC_BIND_SPREAD 94#endif 95 use omp_lib, only : OMP_SCHED_STATIC, OMP_SCHED_DYNAMIC 96 use omp_lib, only : OMP_SCHED_GUIDED, OMP_SCHED_AUTO 97#else 98!$ use omp_lib, only : omp_get_num_threads 99!$ use omp_lib, only : omp_get_schedule, omp_set_schedule 100!$ use omp_lib, only : OMP_SCHED_STATIC, OMP_SCHED_DYNAMIC 101!$ use omp_lib, only : OMP_SCHED_GUIDED, OMP_SCHED_AUTO 102#endif 103 104 implicit none 105 106#ifdef MPI 107 logical :: initialized 108 integer :: MPIerror ! Return error code in MPI routines 109#endif 110 real(dp):: veclen ! Length of a unit-cell vector 111 integer :: i, is, ispin, n_dm_items 112 integer :: neigmin ! Min. number of eigenstates (per k point) 113 integer :: ns ! Number of species 114 logical :: not_using_auxcell 115 real(dp) :: walltime_m, walltime_s 116 integer :: nstates_spin(2) ! dummy intermediate variable 117 integer :: World_Group 118 integer :: npSIESTA 119 external :: read_options 120 external :: reset_messages_file 121 122!---------------------------------------------------------------------- BEGIN 123! Initialise MPI unless siesta is running as a subroutine 124! of a driver program that may have initialized MPI already 125 126#ifdef MPI 127C Initialise MPI and set processor number 128 129 call MPI_Initialized( initialized, MPIerror ) 130 if (.not.initialized) then 131#ifdef _OPENMP 132 call MPI_Init_Thread(MPI_Thread_Funneled, i, MPIerror) 133 if ( MPI_Thread_Funneled /= i ) then 134 ! the requested threading level cannot be asserted 135 ! Notify the user 136 write(*,'(a)') '!!! Could not assert funneled threads' 137 end if 138#else 139 call MPI_Init( MPIerror ) 140#endif 141#ifdef _TRACE_ 142 call MPItrace_shutdown( ) 143#endif 144 endif ! (.not.initialized) 145 146 call MPI_Comm_Rank( MPI_Comm_World, Node, MPIerror ) 147 call MPI_Comm_Size( MPI_Comm_World, Nodes, MPIerror ) 148 149 ! Keeper for PEXSI-runs.. 150 PEXSINodes = Nodes 151 152 call reset_messages_file() 153 154 155 if (.not. fdf_parallel()) then 156 call die('siesta_init: ERROR: FDF module ' // 157 & 'doesn''t have parallel support') 158 endif 159#else 160 call reset_messages_file() 161#endif 162 if (Node==0) then 163 call system("/bin/rm -f 0_NORMAL_EXIT") 164 endif 165 166 ! Initialize the output file 167 ! This will close and open unit=6 in case a command-line 168 ! option is specified 169 call init_output(Node == 0) 170 171#ifdef DEBUG_XC 172! Set output file unit for debug info 173 call setDebugOutputUnit() 174#endif /* DEBUG_XC */ 175 176 ! We have to immediately write out the options used to compile SIESTA 177 if ( Node == 0 ) then 178 ! Print version information 179 call prversion() 180#ifdef MPI 181 ! Sadly PEXSI workers can only be determined after 182 ! fdf has been opened, so we can't print-out that information 183 ! (for now) 184 if ( Nodes > 1 ) then 185 write(6,'(/,a,i0,a)') 186 & '* Running on ', Nodes, ' nodes in parallel' 187 else 188 write(6,'(/,a)') '* Running in serial mode with MPI' 189 endif 190#else 191 write(6,'(/,a)') '* Running in serial mode' 192#endif 193!$OMP parallel default(shared) 194!$OMP master 195!$ i = omp_get_num_threads() 196!$ write(*,'(a,i0,a)') '* Running ',i,' OpenMP threads.' 197!$ write(*,'(a,i0,a)') '* Running ',Nodes*i,' processes.' 198#ifdef _OPENMP 199!$ write(*,'(a,i0)') '* OpenMP version ', _OPENMP 200#endif 201#if _OPENMP >= 201307 202!$ i = omp_get_proc_bind() 203!$ select case ( i ) 204!$ case ( OMP_PROC_BIND_FALSE ) 205!$ write(*,'(a)') '* OpenMP NOT bound (please bind threads!)' 206!$ case ( OMP_PROC_BIND_TRUE ) 207!$ write(*,'(a)') '* OpenMP bound' 208!$ case ( OMP_PROC_BIND_MASTER ) 209!$ write(*,'(a)') '* OpenMP bound (master)' 210!$ case ( OMP_PROC_BIND_CLOSE ) 211!$ write(*,'(a)') '* OpenMP bound (close)' 212!$ case ( OMP_PROC_BIND_SPREAD ) 213!$ write(*,'(a)') '* OpenMP bound (spread)' 214!$ case default 215!$ write(*,'(a)') '* OpenMP bound (unknown)' 216!$ end select 217#endif 218!$ call omp_get_schedule(i,is) 219!$ select case ( i ) 220!$ case ( OMP_SCHED_STATIC ) 221!$ write(*,'(a,i0)') '* OpenMP runtime schedule STATIC, chunks ',is 222!$ case ( OMP_SCHED_DYNAMIC ) 223!$ write(*,'(a,i0)') '* OpenMP runtime schedule DYNAMIC, chunks ',is 224!$ if ( is == 1 ) then 225!$ ! this is the default scheduling, probably the user 226!$ ! have not set the value, predefine it to 32 227!$ is = 32 228!$ write(*,'(a,i0)')'** OpenMP runtime schedule DYNAMIC, chunks ',is 229!$ end if 230!$ case ( OMP_SCHED_GUIDED ) 231!$ write(*,'(a,i0)') '* OpenMP runtime schedule GUIDED, chunks ',is 232!$ case ( OMP_SCHED_AUTO ) 233!$ write(*,'(a,i0)') '* OpenMP runtime schedule AUTO, chunks ',is 234!$ case default 235!$ write(*,'(a,i0)') '* OpenMP runtime schedule UNKNOWN, chunks ',is 236!$ end select 237!$OMP end master 238!$OMP end parallel 239!$ call omp_set_schedule(i,is) 240 241 ! Write time-stamp at the top of the output 242 call timestamp('Start of run') 243 end if 244 245 ! Initialise input/output........................................... 246 call reinit( sname ) 247 248#ifdef MPI 249 ! Optionally, restrict the number of processors that SIESTA uses, out 250 ! of the total pool 251 npSIESTA = fdf_get("MPI.Nprocs.SIESTA",Nodes) 252 call MPI_Comm_Group(MPI_Comm_World, World_Group, MPIerror) 253 call MPI_Group_incl(World_Group, npSIESTA, 254 $ (/ (i,i=0,npSIESTA-1) /), 255 $ SIESTA_Group, MPIerror) 256 call MPI_Comm_create(MPI_Comm_World, SIESTA_Group, 257 $ SIESTA_Comm, MPIerror) 258 259 SIESTA_worker = (Node < npSIESTA) 260 261 ! Swap communicator 262 ! This is needed since MPI_COMM_WORLD is used implicitly by 263 ! Siesta for all operations. 264 MPI_COMM_WORLD = SIESTA_Comm 265 if (SIESTA_worker) then 266 call MPI_Comm_Rank( MPI_Comm_World, Node, MPIerror ) 267 call MPI_Comm_Size( MPI_Comm_World, Nodes, MPIerror ) 268 endif 269#else 270 SIESTA_worker = .true. 271 Node = 0 272 Nodes = 1 273 PEXSINodes = 1 274#endif 275 IOnode = SIESTA_worker .and. (Node .eq. 0) 276 277 if ( IONode ) then 278 ! Add citations 279 call init_citation(trim(slabel)) 280 call add_citation("10.1088/0953-8984/14/11/302") 281 end if 282 283 ! Be sure to initialize the spin-configuration data 284 ! for all processors. 285 call init_spin( ) 286 287 if (.not. SIESTA_worker) RETURN 288 289#ifdef DEBUG 290! Generates a debug file for every process to track the execution. 291! The file is called debug.$(PID) Where PID is the process number. 292! It also works in secuencial mode. 293 call debugMpiOn( ) 294#endif 295#ifdef NCDF_4 296 call ncdf_IONode( IONode ) 297#endif 298 299! Print version information 300 if ( IOnode ) then 301#ifdef MPI 302 if ( PEXSINodes /= Nodes ) then 303 write(6,'(/,a,2(i0,a))') 304 & '* Running on ', Nodes, ' SIESTA-nodes and ', 305 & PEXSINodes, ' PEXSI-nodes in parallel' 306 end if 307#endif 308 call wallclock('Start of run') 309 end if 310 311! Initialize some variables 312 call init_Energies() 313 314 315! Set object debugging level 316 call set_object_debug_level(0) 317 i = fdf_get('DebugObjects.Node',0) 318 if (Node==i) then 319 if (fdf_get('DebugObjects',.false.)) then 320 call set_object_debug_level(1) 321 endif 322 endif 323 324! Initialize CML (relies on reinit) 325 call siesta_cml_init( ) 326 327! Set timer report file and threshold ................................. 328 threshold = fdf_get('timer_report_threshold', 0._dp) 329 call timer_report( file=trim(slabel)//'.times', 330 . threshold=threshold ) 331 ! Note that the parallel timer might be inefficient for reports 332 ! when large numbers of processors are used 333 if ( PEXSINodes /= Nodes ) then 334 use_parallel_timer = fdf_get('UseParallelTimer', .false.) 335 use_tree_timer = fdf_get('UseTreeTimer', .true.) 336 else 337 use_parallel_timer = fdf_get('UseParallelTimer', .true.) 338 use_tree_timer = fdf_get('UseTreeTimer', .false.) 339 end if 340 341! Start time counter 342! Note new placement of this first use, so that 343! initialization is done after the setup of relevant variables 344 345 call timer( 'siesta', 0 ) 346 call timer( 'siesta', 1 ) 347 call timer( 'Setup', 1 ) 348 349! Set allocation report level ......................................... 350! variables level and threshold imported from module siesta_options 351 level = fdf_get('alloc_report_level', 0) 352 threshold = fdf_get('alloc_report_threshold', 0._dp) 353 call alloc_report( level=level, file=trim(slabel)//'.alloc', 354 . threshold=threshold, 355 & printNow=.false. ) 356 357 358! Initialise exchange-correlation functional information 359 call read_xc_info() 360 361! Initialise force field component 362 call inittwobody( ) 363 364! Initialize pseudopotentials and atomic orbitals 365 if (IOnode) call initatom( ns ) 366 call broadcast( ns ) 367 368 call broadcast_basis( ) 369 370 call register_rfs() 371 372 atmonly = fdf_get( 'Atom-Setup-Only', .false. ) 373 if (atmonly) call bye( 'End of atom setup' ) 374 375! Read geometry 376 call struct_init( ) ! Sets na_u, isa, ucell 377 378 ! Initialize the spin-spiral settings 379 call init_spiral( ucell ) 380 381! Initialize atom lists 382 call initatomlists( ) ! Sets iza 383 384! early exit if only checking the structure 385 struct_only = fdf_get( 'Output-Structure-Only', .false. ) 386 if (IONode) then 387 call write_struct( ucell, na_u, isa, iza, xa ) 388 if (fdf_boolean('WriteCoorXmol',.false.)) then 389 call coxmol(iza, xa, na_u) 390 endif 391 if (lUseZmatrix) then 392 call write_canonical_ucell_and_Zmatrix( 393 & filename="OUT.UCELL.ZMATRIX.INITIAL") 394 endif 395 endif 396 if (struct_only) then 397 call bye('End of structure processing') 398 endif 399 400! Walltime control (numbers in seconds) 401 if ( fdf_isphysical("MaxWalltime") ) then 402 walltime_m = fdf_get("MaxWalltime", huge(1._dp), 's') 403 else 404 walltime_m = fdf_get("MaxWalltime", huge(1._dp)) 405 end if 406 ! Period for clean-up operations 407 if ( fdf_isphysical("MaxWalltime.Slack") ) then 408 walltime_s = fdf_get("MaxWalltime.Slack", 5._dp, 's') 409 else 410 walltime_s = fdf_get("MaxWalltime.Slack", 5._dp) 411 end if 412 ! Note that the slack detracts from net available time 413 walltime_warning = walltime_m - walltime_s 414 415!-------------- Now we have the initial geometry 416 417! End of Initial Structure Processing 418 if (Node.eq.0) then 419 write(6,'(/,a,20("*"),a,28("*"))') 420 & 'siesta: ', ' Simulation parameters ' 421 write(6,'(a)') 'siesta:' 422 write(6,'(a)') 'siesta: The following are some of the '// 423 & 'parameters of the simulation.' 424 write(6,'(a)') 'siesta: A complete list of the parameters '// 425 & 'used, including default values,' 426 write(6,'(a,a)')'siesta: can be found in file out.fdf' 427 write(6,'(a)') 'siesta:' 428 endif 429 430! Allocate other arrays based on read sizes 431 nullify(fa,cfa) 432 call re_alloc( fa, 1, 3, 1, na_u, 'fa', 'siesta_init' ) 433 call re_alloc( cfa, 1, 3, 1, na_u, 'cfa', 'siesta_init' ) 434 435! Read simulation data 436 call read_options( na_u, ns, spin%H ) 437 438 call get_allowed_history_depth(n_dm_items) 439 call new(DM_history,n_dm_items,"(DM history stack)") 440 if (ionode) print "(a,i0)", "Size of DM history Fstack: ", 441 $ max_size(DM_history) 442 443 qtot = qtot - charnet ! qtot set in initatomlists 444 ! charnet set in redata 445 if (IOnode) then 446 write(6,fmt="(a,f12.6)") 'Total number of electrons: ', qtot 447 write(6,fmt="(a,f12.6)") 'Total ionic charge: ', zvaltot 448 endif 449! 450! Warn the user: if not doing a direct optimization, the Zmatrix 451! coordinates are no longer updated. Only coordinates are treated. 452! 453 if (lUseZmatrix) then 454 if (idyn .ne. 0) then 455 write(6,"(a)") 456 & " *** WARNING: Zmatrix form will be used only for input !!" 457 write(0,"(a)") 458 & " *** WARNING: Zmatrix form will be used only for input !!" 459 endif 460 endif 461 462! Calculate spin populations for fixed spin case... 463 if (fixspin) then 464 if (.not.spin%Col) 465 & call die( 'siesta: ERROR: ' // 466 & 'You can only fix the spin of the system' // 467 & ' for collinear spin polarized calculations.' ) 468 qtots(1) = (qtot + total_spin) / 2.0_dp 469 qtots(2) = (qtot - total_spin) / 2.0_dp 470 if (IOnode) then 471 write(6,"(a,f12.6)") 'Total number of UP electrons: ', 472 & qtots(1) 473 write(6,"(a,f12.6)") 'Total number of DOWN electrons: ', 474 & qtots(2) 475 end if 476 else 477 ! Distribute electrons equally per spin 478 ! (For spin%spinor = 1, qtots(2) should not be referenced) 479 qtots(:) = qtot / spin%spinor 480 end if 481 482 ! Get number of eigenstates that need to be calculated This 483 ! generally needs to be set to half the number of charges plus a 484 ! few extra states (such that the temperature tail can be found). 485 ! When states are spinors (for collinear or non-collinear spin, 486 ! i.e., for spin%spinor=2), to be occupied by a single electron 487 ! each, this number refers to half the number of spinors. 488 neigwanted = fdf_get('NumberOfEigenStates',no_u) 489 490 ! This is the rough number of states to be occupied (per spin channel) 491 ! For spin%spinor=1, this is half the number of electrons 492 nstates_spin(1:2) = ceiling ( (qtots(1:2)/2.0_dp) * spin%spinor ) 493 494 ! Calculate the number of eigenstates 495 ! dependent on the # of spin-components. 496 neigmin = 0 497 do is = 1 , spin%spinor 498 neigmin = max(neigmin, nstates_spin(is)) 499 end do 500 501 ! Always at least have one additional state per 200 K. This 502 ! choice is rather arbitrary. If the system is large it is likely 503 ! to have many states "close" to the Fermi level, and then the 504 ! number of extra states should be higher. See below 505 506 neigmin = neigmin + ceiling((temp / Kelvin) / 200._dp) 507 508 ! Correct the number of calculated eigenstates 509 if ( neigwanted < 0 ) then 510 511 ! If this number is negative it corresponds to the number of 512 ! states ABOVE the charge neutrality point to be included. For 513 ! normal temperatures and small/medium systems, a number such 514 ! as 10 might be enough, but For very high temperatures and/or 515 ! large systems, this number might need to be increased significantly. 516 517 ! Add the user requested states. 518 neigwanted = neigmin + abs(neigwanted) 519 520 end if 521 522 ! Check number of eigenstates - cannot be larger than number of 523 ! basis functions or smaller than number of occupied states + 1 524 ! so that the Fermi level can be estimated 525 neigwanted = max(neigwanted, neigmin) 526 ! Truncate value to number of basis orbitals 527 neigwanted = min(neigwanted,no_u) 528 529 530! Find maximum interaction range 531 if ( negl ) then 532 rmaxh = 2.0_dp*(rmaxo + rmaxdftu) 533 else 534 rmaxh = 2.0_dp*(rmaxo + max(rmaxkb,rmaxdftu)) 535 endif 536 537! Madelung correction for charged systems 538 if (charnet .ne. 0.0_dp) then 539 call madelung(ucell, shape, charnet, Emad) 540 endif 541 542! Parallel initialisation 543 call initparallel( no_u, na_u, lasto, xa, ucell, rmaxh ) 544 if (IOnode) call show_distribution( ) 545 546! Find number of locally stored orbitals and allocated related arrays 547 call GetNodeOrbs(no_u,Node,Nodes,no_l) 548 549! Find k-grid for Brillouin zone integration 550! NOTE: We need to know whether gamma is .true. or 551! not early, in order to decide whether to use an 552! auxiliary supercell for the calculation of matrix elements. 553 call setup_Kpoint_grid( ucell ) 554 not_using_auxcell = gamma_scf 555 556 ! Read in diagonalization routines 557 ! Note that only the sampled BZ is responsible for 558 ! ParallelOverK option, the remaning options are 559 ! taking care of this option 560 call read_diag(Gamma_scf, spin%H) 561 call print_diag() 562 563 if ( diag_recognizes_neigwanted() ) then 564 if ( IONode .and. neigwanted /= no_u ) 565 & write(6,"(a,t53,'= ',i0,' / ',i0)") 566 & 'diag: Number of calculated states [no_calc / no_u]', 567 & neigwanted, no_u 568 else 569 neigwanted = no_u 570 end if 571 572! Call initialisation of PDOS here since we need to check if 573! the auxiliary supercell is needed for a non-gamma calculation 574 call init_projected_DOS( ) 575 if (do_pdos) then 576 not_using_auxcell = not_using_auxcell.and. gamma_pdos 577 endif 578 579 nullify(eo,qo) 580 call re_alloc(eo, 1, no_u, 1, spin%spinor, 1, nkpnt, 'eo', 581 & 'siesta_init') 582 call re_alloc(qo, 1, no_u, 1, spin%spinor, 1, nkpnt, 'qo', 583 & 'siesta_init') 584 585 call setup_bands( ) 586 not_using_auxcell = not_using_auxcell .and. gamma_bands 587 588 call setup_wf_kpoints( ) 589 not_using_auxcell = not_using_auxcell .and. gamma_wavefunctions 590 591 call estimate_pol_kpoints( ucell ) 592 not_using_auxcell = not_using_auxcell .and. gamma_polarization 593! 594! User can request that the calculation is done with an explicit 595! auxiliary supercell and hermitian version, even if using only 596! the gamma point 597! 598! In case the user has requested the use of auxiliary 599! cell we simply use that option. 600 if ( .not. not_using_auxcell ) use_aux_cell = .true. 601 602! Find required supercell 603! 2*rmaxh is used to guarantee that two given orbitals in the 604! supercell can only overlap once 605 if ( .not. use_aux_cell ) then 606 nsc(1:3) = 1 607 else 608 ! Determine the tightest auxiliary supercell using 609 ! also the atomic positions 610 call exact_sc_ag(negl,ucell,na_u,isa,xa,nsc) 611 end if 612 613 mscell = 0.0_dp 614 do i = 1, 3 615 mscell(i,i) = nsc(i) 616 nsc_old(i) = nsc(i) 617 enddo 618 619! Find auxiliary supercell (required only for k sampling) 620 call superc( ucell, scell, nsc ) 621 622! Initialise metadynamic forces if required 623 call initmeta( ) 624 625 if (idyn .eq. 0) then 626 inicoor = 0 627 fincoor = nmove 628 else if (idyn .ge. 1 .and. idyn .le. 5) then 629 inicoor = istart 630 fincoor = ifinal 631 else if (idyn .eq. 6) then 632 inicoor = 0 633 fincoor = (ia2-ia1+1)*3*2 634 else if (idyn .eq. 7) then 635 call die( "'PHONON' support is deprecated" ) 636 else if (idyn .eq. 8) then 637 inicoor = 0 638 fincoor = huge(1) 639#ifdef NCDF_4 640 else if (idyn == 9) then 641 if ( IONode ) then 642 write(*,'(/,a)')'expcoord: '//repeat('*',69) 643 write(*,'(a,i0)') 644 & 'expcoord: Currently reached coordinate step = ', 645 & inicoor 646 write(*,'(a,i0)') 647 & 'expcoord: Final coordinate step = ', 648 & fincoor 649 write(*,'(a,/)')'expcoord: '//repeat('*',69) 650 end if 651#endif 652#ifdef SIESTA__FLOOK 653 else if (idyn == 10) then 654 inicoor = 0 655 ! Controlled in external LUA machine 656 fincoor = huge(1) 657#endif 658 else 659 call die( 'siesta: wrong idyn' ) 660 endif 661 662 ! initialize the fixed positions (we need to check whether the 663 ! fixed positions make sense) 664 call init_fixed( ucell, na_u , isa, iza ) 665 call print_fixed( ) 666 667 ! Build initial velocities according to Maxwell-Bolzmann distribution.... 668 if (idyn .ne. 0 .and. idyn .ne. 6 .and. (.not. xv_file_read)) 669 & call vmb( na_u, tempinit, amass, xa, isa, va ) 670 671 istp = 0 672 673 call ts_init(spin%Grid,ucell,na_u,xa,lasto,no_u,inicoor,fincoor) 674 675#ifdef SIESTA__FLOOK 676 677 ! Populate dictionaries 678 call dict_populate() 679 680 ! Initialize LUA handle, also die if the Lua relaxation 681 ! is requested 682 call slua_init(LUA, idyn == 10) 683 684 ! Call a preprocess step right after initialization 685 call slua_call(LUA,LUA_INITIALIZE) 686 687#endif 688 689! Output memory use before main loop 690!! call printmemory( 6, 0 ) 691 692! Initialization now complete. Flush stdout. 693 if (ionode) call pxfflush( 6 ) 694 695 call timer( 'Setup', 2 ) 696 697!--------------------------------------------------------------------------- END 698 699 END subroutine siesta_init 700 701 END MODULE m_siesta_init 702