TABLE OF CONTENTS
- ABINIT/m_rttddft_tdks
- m_rttddft_tdks/first_setup
- m_rttddft_tdks/read_wfk
- m_rttddft_tdks/second_setup
- m_rttddft_tdks/tdks_free
- m_rttddft_tdks/tdks_init
ABINIT/m_rttddft_tdks [ Modules ]
NAME
m_rttddft_tdks
FUNCTION
Contains the main object (tdks) to propagate the time-dependent Kohn-Sham equations in RT-TDDFT
COPYRIGHT
Copyright (C) 2021-2024 ABINIT group (FB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_rttddft_tdks 24 25 use defs_basis 26 use defs_abitypes, only: MPI_type 27 use defs_datatypes, only: pseudopotential_type, ebands_t 28 use defs_wvltypes, only: wvl_data, nullify_wvl_data 29 30 use libxc_functionals, only: libxc_functionals_get_hybridparams 31 use m_bandfft_kpt, only: bandfft_kpt, bandfft_kpt_init1, & 32 & bandfft_kpt_destroy_array 33 use m_cgprj, only: ctocprj 34 use m_common, only: setup1 35 use m_dtfil, only: datafiles_type 36 use m_dtset, only: dataset_type 37 use m_ebands, only: ebands_from_dtset, ebands_free, unpack_eneocc 38 use m_energies, only: energies_type, energies_init 39 use m_errors, only: msg_hndl, assert 40 use m_extfpmd, only: extfpmd_type 41 use m_gemm_nonlop, only: init_gemm_nonlop, destroy_gemm_nonlop 42 use m_geometry, only: fixsym 43 use m_hdr, only: hdr_type, hdr_init 44 use m_initylmg, only: initylmg 45 use m_invovl, only: init_invovl, destroy_invovl 46 use m_io_tools, only: open_file 47 use m_inwffil, only: inwffil 48 use m_kg, only: kpgio, getph, getcut 49 use m_mpinfo, only: proc_distrb_cycle 50 use m_occ, only: newocc 51 use m_paw_an, only: paw_an_type, paw_an_init, paw_an_free, & 52 & paw_an_nullify 53 use m_pawang, only: pawang_type 54 use m_pawcprj, only: pawcprj_type,pawcprj_free,pawcprj_alloc, & 55 & pawcprj_getdim 56 use m_paw_dmft, only: init_sc_dmft,destroy_sc_dmft,paw_dmft_type 57 use m_pawfgr, only: pawfgr_type, pawfgr_init, pawfgr_destroy 58 use m_pawfgrtab, only: pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 59 use m_paw_init, only: pawinit,paw_gencond 60 use m_paw_ij, only: paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 61 use m_paw_nhat, only: nhatgrid 62 use m_paw_occupancies, only: initrhoij 63 use m_pawrad, only: pawrad_type 64 use m_pawrhoij, only: pawrhoij_type, pawrhoij_copy, pawrhoij_free 65 use m_paw_sphharm, only: setsym_ylm 66 use m_pawtab, only: pawtab_type, pawtab_get_lsize 67 use m_paw_tools, only: chkpawovlp 68 use m_pspini, only: pspini 69 use m_profiling_abi, only: abimem_record 70 use m_spacepar, only: setsym 71 use m_specialmsg, only: wrtout 72 use m_symtk, only: symmetrize_xred 73 use m_wffile, only: wffile_type, WffClose 74 use m_xmpi, only: xmpi_bcast, xmpi_sum 75 76 implicit none 77 78 private
m_rttddft_tdks/first_setup [ Functions ]
[ Top ] [ m_rttddft_tdks ] [ Functions ]
NAME
first_setup
FUNCTION
Intialize many important quantities before running RT-TDDFT (PW, FFT, PSP, Symmetry etc.)
INPUTS
codvsn = code version dtfil <type datafiles_type> = infos about file names, file unit numbers dtset <type(dataset_type)> = all input variables for this dataset mpi_enreg <MPI_type> = MPI-parallelisation information pawrad(ntypat*usepaw) <type(pawrad_type)> = paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)> = paw tabulated starting data psps <type(pseudopotential_type)> = variables related to pseudopotentials tdks <type(tdks_type)> = the tdks object to initialize
OUTPUT
psp_gencond <integer> = store conditions for generating psp ecut_eff <real(dp)> = effective PW cutoff energy
SOURCE
474 subroutine first_setup(codvsn,dtfil,dtset,ecut_eff,mpi_enreg,pawrad,pawtab,psps,psp_gencond,tdks) 475 476 implicit none 477 478 !Arguments ------------------------------------ 479 !scalars 480 character(len=8), intent(in) :: codvsn 481 integer, intent(out) :: psp_gencond 482 real(dp), intent(out) :: ecut_eff 483 type(datafiles_type), intent(in) :: dtfil 484 type(dataset_type), intent(inout) :: dtset 485 type(pseudopotential_type), intent(inout) :: psps 486 type(MPI_type), intent(inout) :: mpi_enreg 487 type(tdks_type), intent(inout) :: tdks 488 !arrays 489 type(pawrad_type), intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 490 type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 491 492 !Local variables------------------------------- 493 !scalars 494 integer,parameter :: response=0, cplex=1 495 integer :: comm_psp 496 integer :: gscase 497 integer :: iatom, ierr, itypat, indx 498 integer :: mgfftf, my_natom 499 integer :: npwmin, nfftot 500 integer :: ylm_option 501 real(dp) :: gsqcut_eff, gsqcutc_eff 502 real(dp) :: ecutdg_eff 503 type(ebands_t) :: bstruct 504 !arrays 505 character(len=500) :: msg 506 integer, allocatable :: npwarr_(:) 507 integer :: ngfft(18) 508 integer :: ngfftf(18) 509 integer :: npwtot(dtset%nkpt) 510 511 ! *********************************************************************** 512 513 my_natom=mpi_enreg%my_natom 514 515 !** Init FFT grid(s) sizes (be careful !) 516 !See NOTES in the comments at the beginning of this file. 517 tdks%nfft = dtset%nfft 518 call pawfgr_init(tdks%pawfgr,dtset,mgfftf,tdks%nfftf,ecut_eff,ecutdg_eff, & 519 & ngfft,ngfftf) 520 521 !** Init to zero different energies 522 call energies_init(tdks%energies) 523 tdks%ecore = zero 524 tdks%etot = zero 525 526 !** various additional setup mostly related to fft grids and the box (rprimd, metric..) 527 call setup1(dtset%acell_orig,tdks%bantot,dtset,ecutdg_eff,ecut_eff,tdks%gmet, & 528 & tdks%gprimd,gsqcut_eff,gsqcutc_eff,ngfftf,ngfft,dtset%nkpt, & 529 & dtset%nsppol,response,tdks%rmet,dtset%rprim_orig,tdks%rprimd, & 530 & tdks%ucvol,psps%usepaw) 531 532 !FB: @MT Needed? 533 !!In some cases (e.g. getcell/=0), the plane wave vectors have 534 !! to be generated from the original simulation cell 535 !rprimd_for_kg=rprimd 536 !if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=args_gs%rprimd_orig 537 !call matr3inv(rprimd_for_kg,gprimd_for_kg) 538 !gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg) 539 540 !** Set up the basis sphere of planewaves 541 ABI_MALLOC(tdks%npwarr,(dtset%nkpt)) 542 ABI_MALLOC(tdks%kg,(3,dtset%mpw*dtset%mkmem)) 543 call kpgio(ecut_eff,dtset%exchn2n3d,tdks%gmet,dtset%istwfk,tdks%kg,dtset%kptns, & 544 & dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw, & 545 & tdks%npwarr,npwtot,dtset%nsppol) 546 call bandfft_kpt_init1(bandfft_kpt,dtset%istwfk,tdks%kg,dtset%mgfft,dtset%mkmem, & 547 & mpi_enreg,dtset%mpw,dtset%nband,dtset%nkpt,tdks%npwarr, & 548 & dtset%nsppol) 549 550 !** Use efficient BLAS calls for computing the non local potential 551 if(dtset%use_gemm_nonlop == 1 .and. dtset%gpu_option/=ABI_GPU_DISABLED) then 552 ! set global variable 553 tdks%gemm_nonlop_use_gemm = .true. 554 call init_gemm_nonlop(dtset%nkpt) 555 else 556 tdks%gemm_nonlop_use_gemm = .false. 557 end if 558 559 !** TODO: uncomment when gemm_nonlop can be used on GPU 560 ! if(dtset%use_gemm_nonlop == 1 .and. dtset%gpu_option/=ABI_GPU_DISABLED) then 561 ! ! set global variable 562 ! tdks%gemm_nonlop_use_gemm_gpu = .true. 563 ! !call init_gemm_nonlop_gpu(dtset%nkpt) 564 ! else 565 ! tdks%gemm_nonlop_use_gemm_gpu = .false. 566 ! end if 567 568 !** Setup the Ylm for each k point 569 if (psps%useylm==1) then 570 ABI_MALLOC(tdks%ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)) 571 ABI_MALLOC(tdks%ylmgr,(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)) 572 ylm_option=0 573 if ((dtset%prtstm==0.and.dtset%iscf>0.and.dtset%positron/=1) .or. & 574 & (dtset%berryopt==4 .and. dtset%optstress /= 0 .and. psps%usepaw==1) .or. & 575 & (dtset%orbmag<0 .and. psps%usepaw==1)) then 576 ylm_option = 1 ! gradients of YLM 577 end if 578 call initylmg(tdks%gprimd,tdks%kg,dtset%kptns,dtset%mkmem,mpi_enreg,& 579 & psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,& 580 & tdks%npwarr,dtset%nsppol,ylm_option,tdks%rprimd,tdks%ylm,tdks%ylmgr) 581 else 582 ABI_MALLOC(tdks%ylm,(0,0)) 583 ABI_MALLOC(tdks%ylmgr,(0,0,0)) 584 end if 585 586 !** Open and read pseudopotential files 587 comm_psp=mpi_enreg%comm_cell 588 call pspini(dtset,dtfil,tdks%ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad, & 589 & pawtab,psps,tdks%rprimd,comm_mpi=comm_psp) 590 591 !In case of isolated computations, ecore must be set to zero 592 !because its contribution is counted in the ewald energy as the ion-ion interaction. 593 if (dtset%icoulomb == 1) tdks%ecore = zero 594 595 !Include core energy? 596 select case(dtset%usepotzero) 597 case(0,1) 598 tdks%energies%e_corepsp = tdks%ecore / tdks%ucvol 599 tdks%energies%e_corepspdc = zero 600 case(2) 601 ! No need to include the PspCore energy since it is already included in the 602 ! local pseudopotential (vpsp) 603 tdks%energies%e_corepsp = zero 604 tdks%energies%e_corepspdc = zero 605 end select 606 607 !** Initialize band structure datatype 608 ABI_MALLOC(npwarr_,(dtset%nkpt)) 609 npwarr_(:)=tdks%npwarr(:) 610 if (dtset%paral_kgb/=0) then 611 call xmpi_sum(npwarr_,mpi_enreg%comm_bandfft,ierr) 612 end if 613 bstruct = ebands_from_dtset(dtset, npwarr_) 614 ABI_FREE(npwarr_) 615 call unpack_eneocc(dtset%nkpt,dtset%nsppol,bstruct%mband,bstruct%nband,dtset%occ_orig(:,1),bstruct%occ,val=zero) 616 617 !** Initialize PAW atomic occupancies 618 ABI_MALLOC(tdks%pawrhoij,(my_natom*psps%usepaw)) 619 if (psps%usepaw == 1) then 620 call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu,my_natom,dtset%natom, & 621 & dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat, & 622 & tdks%pawrhoij,dtset%pawspnorb,pawtab,cplex,dtset%spinat, & 623 & dtset%typat,comm_atom=mpi_enreg%comm_atom, & 624 & mpi_atmtab=mpi_enreg%my_atmtab) 625 end if 626 627 !Nullify wvl_data. It is important to do so irregardless of the value of usewvl 628 !Only needed here because hdr%init requires a wvl object for the wvl%descr input 629 call nullify_wvl_data(tdks%wvl) 630 631 !** Initialize header 632 gscase=0 633 call hdr_init(bstruct,codvsn,dtset,tdks%hdr,pawtab,gscase,psps,tdks%wvl%descr,& 634 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 635 636 !Clean band structure datatype 637 call ebands_free(bstruct) 638 639 !** PW basis set: test if the problem is ill-defined. 640 npwmin=minval(tdks%hdr%npwarr(:)) 641 if (dtset%mband > npwmin) then 642 ! No way we can solve the problem. Abort now! 643 write(msg,"(2(a,i0),4a)") "Number of bands nband= ",dtset%mband, & 644 & " > number of planewaves npw= ",npwmin,ch10, & 645 & "The number of eigenvectors cannot be greater that the size of the Hamiltonian!",& 646 & ch10, "Action: decrease nband or, alternatively, increase ecut" 647 if (dtset%ionmov/=23) then 648 ABI_ERROR(msg) 649 else 650 ABI_WARNING(msg) 651 end if 652 653 else if (dtset%mband >= 0.9 * npwmin) then 654 ! Warn the user 655 write(msg,"(a,i0,a,f6.1,4a)") "Number of bands nband= ",dtset%mband, & 656 & " >= 0.9 * maximum number of planewaves= ",0.9*npwmin,ch10,& 657 & "This could lead to some instabilities, you might want to decrease nband or increase ecut!", & 658 & ch10,"Assume experienced user. Execution will continue." 659 ABI_WARNING(msg) 660 end if 661 662 !** Initialize symmetry 663 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 664 ABI_MALLOC(tdks%irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 665 ABI_MALLOC(tdks%phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 666 ABI_MALLOC(tdks%indsym,(4,dtset%nsym,dtset%natom)) 667 ABI_MALLOC(tdks%symrec,(3,3,dtset%nsym)) 668 tdks%irrzon(:,:,:)=0 669 tdks%phnons(:,:,:)=zero 670 tdks%indsym(:,:,:)=0 671 tdks%symrec(:,:,:)=0 672 673 !FB TODO: Should symmetry be used when ions are moving? Modify if Ehrenfest dynamics 674 !Do symmetry stuff if nsym>1 675 if (dtset%nsym>1) then 676 call setsym(tdks%indsym,tdks%irrzon,dtset%iscf,dtset%natom, & 677 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym, & 678 & tdks%phnons,dtset%symafm,tdks%symrec,dtset%symrel, & 679 & dtset%tnons,dtset%typat,dtset%xred_orig) 680 681 !Make sure dtset%iatfix does not break symmetry 682 call fixsym(dtset%iatfix,tdks%indsym,dtset%natom,dtset%nsym) 683 else 684 !The symrec array is used by initberry even in case nsym = 1 - FB: @MT Needed ? 685 tdks%symrec(:,:,1) = 0 686 tdks%symrec(1,1,1) = 1 ; tdks%symrec(2,2,1) = 1 ; tdks%symrec(3,3,1) = 1 687 end if 688 689 !** Initialize and eventually symmetrize reduced atomic coordinates 690 ABI_MALLOC(tdks%xred,(3,dtset%natom,dtset%nimage)) 691 tdks%xred = dtset%xred_orig 692 !FB: Should we? 693 !Eventually symmetrize atomic coordinates over space group elements 694 call symmetrize_xred(dtset%natom,dtset%nsym,dtset%symrel,dtset%tnons,tdks%xred, & 695 & indsym=tdks%indsym) 696 697 !** Create the atindx array 698 !** index table of atoms, in order for them to be used type after type. 699 ABI_MALLOC(tdks%atindx,(dtset%natom)) 700 ABI_MALLOC(tdks%atindx1,(dtset%natom)) 701 ABI_MALLOC(tdks%nattyp,(psps%ntypat)) 702 indx=1 703 do itypat=1,psps%ntypat 704 tdks%nattyp(itypat)=0 705 do iatom=1,dtset%natom 706 if(dtset%typat(iatom)==itypat)then 707 tdks%atindx(iatom)=indx 708 tdks%atindx1(indx)=iatom 709 indx=indx+1 710 tdks%nattyp(itypat)=tdks%nattyp(itypat)+1 711 end if 712 end do 713 end do 714 715 !** Calculate zion: the total positive charge acting on the valence electrons 716 tdks%zion=zero 717 do iatom=1,dtset%natom 718 tdks%zion=tdks%zion+psps%ziontypat(dtset%typat(iatom)) 719 end do 720 721 !FB: probably not needed 722 !Further setup 723 !call setup2(dtset,npwtot,start,tdks%wvl%wfs,tdks%xred) 724 725 end subroutine first_setup
m_rttddft_tdks/read_wfk [ Functions ]
[ Top ] [ m_rttddft_tdks ] [ Functions ]
NAME
read_wfk
FUNCTION
Reads initial wavefunctions (KS orbitals) in WFK file (call inwffil)
INPUTS
dtfil <type datafiles_type> = infos about file names, file unit numbers dtset <type(dataset_type)> = all input variables for this dataset ecut_eff <real(dp)> = effective PW cutoff energy mpi_enreg <MPI_type> = MPI-parallelisation information tdks <type(tdks_type)> = the tdks object to initialize
OUTPUT
SOURCE
1024 subroutine read_wfk(dtfil, dtset, ecut_eff, fname_wfk, mpi_enreg, tdks) 1025 1026 implicit none 1027 1028 !Arguments ------------------------------------ 1029 !scalars 1030 character(len=fnlen), intent(in) :: fname_wfk 1031 real(dp), intent(in) :: ecut_eff 1032 type(datafiles_type), intent(in) :: dtfil 1033 type(dataset_type), intent(inout) :: dtset 1034 type(MPI_type), intent(inout) :: mpi_enreg 1035 type(tdks_type), intent(inout) :: tdks 1036 1037 !Local variables------------------------------- 1038 !scalars 1039 integer,parameter :: formeig=0 1040 integer :: ask_accurate 1041 integer :: band 1042 integer :: cnt 1043 integer :: ierr, ikpt 1044 integer :: my_nspinor 1045 integer :: optorth 1046 integer :: spin 1047 type(wffile_type) :: wff1, wffnow 1048 !arrays 1049 character(len=500) :: msg 1050 1051 ! *********************************************************************** 1052 1053 !If paral_kgb == 0, it may happen that some processors are idle (no entry in proc_distrb) 1054 !but mkmem == nkpt and this can cause integer overflow in mcg or allocation error. 1055 !Here we count the number of states treated by the proc. if cnt == 0, mcg is then set to 0. 1056 cnt = 0 1057 do spin=1,dtset%nsppol 1058 do ikpt=1,dtset%nkpt 1059 do band=1,dtset%nband(ikpt + (spin-1) * dtset%nkpt) 1060 if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, band, band, spin, mpi_enreg%me_kpt)) cnt = cnt + 1 1061 end do 1062 end do 1063 end do 1064 1065 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 1066 tdks%mcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 1067 if (cnt == 0) then 1068 tdks%mcg = 0 1069 write(msg,"(2(a,i0))")"rank: ",mpi_enreg%me, "does not have wavefunctions to treat. Setting mcg to: ",tdks%mcg 1070 ABI_WARNING(msg) 1071 end if 1072 1073 if (dtset%usewvl == 0 .and. dtset%mpw > 0 .and. cnt /= 0)then 1074 if (my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol > floor(real(HUGE(0))/real(dtset%mpw) )) then 1075 ierr = 0 1076 write (msg,'(9a)')& 1077 & "Default integer is not wide enough to store the size of the wavefunction array (mcg).",ch10,& 1078 & "This usually happens when paral_kgb == 0 and there are not enough procs to distribute kpts and spins",ch10,& 1079 & "Action: if paral_kgb == 0, use nprocs = nkpt * nsppol to reduce the memory per node.",ch10,& 1080 & "If tdks does not solve the problem, use paral_kgb 1 with nprocs > nkpt * nsppol and use npfft/npband/npspinor",ch10,& 1081 & "to decrease the memory requirements. Consider also OpenMP threads." 1082 ABI_ERROR_NOSTOP(msg,ierr) 1083 write (msg,'(5(a,i0), 2a)')& 1084 & "my_nspinor: ",my_nspinor, ", mpw: ",dtset%mpw, ", mband: ",dtset%mband,& 1085 & ", mkmem: ",dtset%mkmem, ", nsppol: ",dtset%nsppol,ch10,& 1086 & 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS...) compiled in int64 mode' 1087 ABI_ERROR(msg) 1088 end if 1089 end if 1090 1091 ! Alloc size for wfk and bands 1092 ABI_MALLOC_OR_DIE(tdks%cg,(2,tdks%mcg),ierr) 1093 ABI_MALLOC(tdks%eigen,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1094 ABI_MALLOC(tdks%eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1095 1096 tdks%eigen(:) = zero 1097 ask_accurate=0 1098 1099 !Actually read the intial KS orbitals here 1100 if (dtset%td_restart /= 1) then 1101 write(msg,'(3a)') ch10,'------------- Reading initial wavefunctions -------------',ch10 1102 else 1103 write(msg,'(3a)') ch10,'------------- Reading wavefunctions for restart -------------',ch10 1104 end if 1105 call wrtout(ab_out,msg) 1106 if (do_write_log) call wrtout(std_out,msg) 1107 wff1%unwff=dtfil%unwff1 1108 optorth=0 !No need to orthogonalize the wfk 1109 tdks%hdr%rprimd=tdks%rprimd 1110 tdks%cg=zero 1111 call inwffil(ask_accurate,tdks%cg,dtset,dtset%ecut,ecut_eff,tdks%eigen, & 1112 & dtset%exchn2n3d,formeig,tdks%hdr,1,dtset%istwfk,tdks%kg, & 1113 & dtset%kptns,dtset%localrdwf,dtset%mband,tdks%mcg,dtset%mkmem, & 1114 & mpi_enreg,dtset%mpw,dtset%nband,tdks%pawfgr%ngfft,dtset%nkpt, & 1115 & tdks%npwarr,dtset%nsppol,dtset%nsym,dtset%occ_orig,optorth, & 1116 & dtset%symafm,dtset%symrel,dtset%tnons,dtfil%unkg,wff1,wffnow, & 1117 & dtfil%unwff1,fname_wfk,tdks%wvl) 1118 1119 !Close wff1 1120 call WffClose(wff1,ierr) 1121 1122 !Keep initial eigenvalues in memory 1123 tdks%eigen0(:) = tdks%eigen(:) 1124 1125 !In case of restart also read wfk file containing wave functions at t=0 1126 if (dtset%td_restart == 1 .and. tdks%fname_wfk0 /= fname_wfk) then 1127 write(msg,'(3a)') ch10,'------------- Reading initial wavefunctions -------------',ch10 1128 ABI_MALLOC_OR_DIE(tdks%cg0,(2,tdks%mcg),ierr) 1129 call wrtout(ab_out,msg) 1130 if (do_write_log) call wrtout(std_out,msg) 1131 tdks%cg0=zero 1132 call inwffil(ask_accurate,tdks%cg0,dtset,dtset%ecut,ecut_eff,tdks%eigen0, & 1133 & dtset%exchn2n3d,formeig,tdks%hdr,1,dtset%istwfk,tdks%kg, & 1134 & dtset%kptns,dtset%localrdwf,dtset%mband,tdks%mcg,dtset%mkmem, & 1135 & mpi_enreg,dtset%mpw,dtset%nband,tdks%pawfgr%ngfft,dtset%nkpt, & 1136 & tdks%npwarr,dtset%nsppol,dtset%nsym,dtset%occ_orig,optorth, & 1137 & dtset%symafm,dtset%symrel,dtset%tnons,dtfil%unkg,wff1,wffnow, & 1138 & dtfil%unwff1,tdks%fname_wfk0,tdks%wvl) 1139 end if 1140 1141 end subroutine read_wfk
m_rttddft_tdks/second_setup [ Functions ]
[ Top ] [ m_rttddft_tdks ] [ Functions ]
NAME
second_setup
FUNCTION
Further important initialization required after reading WFK and computing occupation numbers in paticular related to PAW
INPUTS
dtset <type(dataset_type)> = all input variables for this dataset mpi_enreg <MPI_type> = MPI-parallelisation information pawang <type(pawang_type)> = paw angular mesh and related data pawrad(ntypat*usepaw) <type(pawrad_type)> = paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)> = paw tabulated starting data psps <type(pseudopotential_type)> = variables related to pseudopotentials psp_gencond <integer> = store conditions for generating psp tdks <type(tdks_type)> = the tdks object to initialize
OUTPUT
SOURCE
750 subroutine second_setup(dtset, mpi_enreg, pawang, pawrad, pawtab, psps, psp_gencond, tdks) 751 752 implicit none 753 754 !Arguments ------------------------------------ 755 !scalars 756 integer, intent(in) :: psp_gencond 757 type(pawang_type), intent(inout) :: pawang 758 type(dataset_type), intent(inout) :: dtset 759 type(pseudopotential_type), intent(inout) :: psps 760 type(MPI_type), intent(inout) :: mpi_enreg 761 !arrays 762 type(pawrad_type), intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 763 type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 764 type(tdks_type), intent(inout) :: tdks 765 766 !Local variables------------------------------- 767 !scalars 768 logical :: call_pawinit 769 integer, parameter :: cplex = 1 770 integer :: forces_needed 771 integer :: gnt_option 772 integer :: has_dijhat, has_vhartree, has_dijfock 773 integer :: has_dijnd, has_dijU, has_vxctau 774 integer :: my_natom, my_nspinor 775 integer :: ncpgr 776 integer :: optcut, optgr0, optgr1, optgr2, optrad 777 integer :: stress_needed 778 integer :: use_hybcomp 779 real(dp) :: gsqcut_shp 780 real(dp) :: hyb_range_fock 781 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 782 !arrays 783 integer,allocatable :: l_size_atm(:) 784 785 ! *********************************************************************** 786 787 my_natom=mpi_enreg%my_natom 788 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 789 790 !FB: @MT needed? 791 if (psps%usepaw==1) then 792 call pawrhoij_copy(tdks%hdr%pawrhoij,tdks%pawrhoij,comm_atom=mpi_enreg%comm_atom, & 793 & mpi_atmtab=mpi_enreg%my_atmtab) 794 end if 795 796 !FB: Needed because paw_dmft is required in mkrho 797 !PAW related operations 798 !Initialize paw_dmft, even if neither dmft not paw are used 799 call init_sc_dmft(dtset%nbandkss,dtset%dmftbandi,dtset%dmftbandf, & 800 & dtset%dmft_read_occnd,dtset%mband,dtset%nband,dtset%nkpt, & 801 & dtset%nspden,dtset%nspinor,dtset%nsppol,tdks%occ0,dtset%usedmft,& 802 & tdks%paw_dmft,dtset%usedmft,dtset%dmft_solv,mpi_enreg) 803 804 !*** Main PAW initialization *** 805 tdks%mcprj=0;tdks%mband_cprj=0 806 if(psps%usepaw==1) then 807 gnt_option=1 808 if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 809 810 !** Test if we have to call pawinit 811 ! Some gen-cond have to be added... 812 call paw_gencond(dtset,gnt_option,"test",call_pawinit) 813 814 if (psp_gencond==1.or.call_pawinit) then 815 gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 816 hyb_range_fock=zero 817 if (dtset%ixc<0) then 818 call libxc_functionals_get_hybridparams(hyb_range=hyb_range_fock) 819 end if 820 call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,hyb_range_fock, & 821 & dtset%pawlcutd,dtset%pawlmix,psps%mpsang,dtset%pawnphi, & 822 & dtset%nsym,dtset%pawntheta,pawang,pawrad,dtset%pawspnorb, & 823 & pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero) 824 825 ! Update internal values 826 call paw_gencond(dtset,gnt_option,"save",call_pawinit) 827 end if 828 psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore) 829 call setsym_ylm(tdks%gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol, & 830 & tdks%rprimd,tdks%symrec,pawang%zarot) 831 832 !** Initialisation of cprj 833 tdks%mband_cprj=dtset%mband 834 if (dtset%paral_kgb/=0) tdks%mband_cprj=tdks%mband_cprj/mpi_enreg%nproc_band 835 tdks%mcprj=my_nspinor*tdks%mband_cprj*dtset%mkmem*dtset%nsppol 836 ABI_MALLOC(tdks%cprj,(dtset%natom,tdks%mcprj)) 837 ncpgr=0 838 !FB: @MT dimcprj_srt needed? 839 ABI_MALLOC(tdks%dimcprj,(dtset%natom)) 840 !ABI_MALLOC(dimcprj_srt,(dtset%natom)) 841 call pawcprj_getdim(tdks%dimcprj,dtset%natom,tdks%nattyp,dtset%ntypat, & 842 & dtset%typat,pawtab,'R') 843 !call pawcprj_getdim(dimcprj_srt,dtset%natom,tdks%nattyp,dtset%ntypat, & 844 ! & dtset%typat,pawtab,'O') 845 !call pawcprj_alloc(tdks%cprj,ncpgr,dimcprj_srt) 846 call pawcprj_alloc(tdks%cprj,ncpgr,tdks%dimcprj) 847 !ABI_FREE(dimcprj_srt) 848 849 !** Variables/arrays related to the fine FFT grid 850 ABI_MALLOC(tdks%pawfgrtab,(my_natom)) 851 if (my_natom>0) then 852 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat, & 853 & mpi_atmtab=mpi_enreg%my_atmtab) 854 call pawfgrtab_init(tdks%pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat, & 855 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 856 ABI_FREE(l_size_atm) 857 end if 858 tdks%usexcnhat=maxval(pawtab(:)%usexcnhat) 859 860 !** Variables/arrays related to the PAW spheres 861 ABI_MALLOC(tdks%paw_ij,(my_natom)) 862 ABI_MALLOC(tdks%paw_an,(my_natom)) 863 call paw_an_nullify(tdks%paw_an) 864 call paw_ij_nullify(tdks%paw_ij) 865 has_dijhat=0; if (dtset%iscf==22) has_dijhat=1 866 has_vhartree=0; if (dtset%prtvha > 0 .or. dtset%prtvclmb > 0) has_vhartree=1 867 has_dijnd=0;if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1 868 has_dijfock=0 869 has_dijU=merge(0,1,dtset%usepawu>0) !Be careful on this! 870 has_vxctau=dtset%usekden 871 call paw_an_init(tdks%paw_an,dtset%natom,dtset%ntypat,0,0,dtset%nspden, & 872 & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1, & 873 & has_vxctau=has_vxctau,has_vxc_ex=1,has_vhartree=has_vhartree, & 874 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 875 call paw_ij_init(tdks%paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden, & 876 & dtset%pawspnorb,dtset%natom,dtset%ntypat,dtset%typat,pawtab, & 877 & has_dij=1,has_dijfock=has_dijfock,has_dijhartree=1, & 878 & has_dijnd=has_dijnd,has_dijso=1,has_dijhat=has_dijhat, & 879 & has_dijU=has_dijU,has_pawu_occ=1,has_exexch_pot=1, & 880 & nucdipmom=dtset%nucdipmom,comm_atom=mpi_enreg%comm_atom, & 881 & mpi_atmtab=mpi_enreg%my_atmtab) 882 883 !** Check for non-overlapping spheres 884 call chkpawovlp(dtset%natom,psps%ntypat,dtset%pawovlp,pawtab,tdks%rmet, & 885 & dtset%typat,tdks%xred) 886 887 !** Identify parts of the rectangular grid where the density has to be calculated 888 optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm 889 forces_needed=0 !FB TODO Maybe needs to be changed if Ehrenfest? 890 stress_needed=0 891 if ((forces_needed==1) .or. & 892 & (dtset%xclevel==2 .and. dtset%pawnhatxc>0 .and. tdks%usexcnhat>0) .or. & 893 & (dtset%positron/=0.and.forces_needed==2)) then 894 optgr1=dtset%pawstgylm 895 if (stress_needed==1) optrad=1; if (dtset%pawprtwf==1) optrad=1 896 end if 897 call nhatgrid(tdks%atindx1,tdks%gmet,my_natom,dtset%natom, & 898 & tdks%nattyp,tdks%pawfgr%ngfft,psps%ntypat,optcut,optgr0,optgr1, & 899 & optgr2,optrad,tdks%pawfgrtab,pawtab,tdks%rprimd,dtset%typat, & 900 & tdks%ucvol,tdks%xred,comm_atom=mpi_enreg%comm_atom, & 901 & mpi_atmtab=mpi_enreg%my_atmtab,comm_fft=mpi_enreg%comm_fft, & 902 & distribfft=mpi_enreg%distribfft) 903 904 tdks%nhatgrdim=0;if (dtset%xclevel==2) tdks%nhatgrdim=tdks%usexcnhat*dtset%pawnhatxc 905 if (tdks%nhatgrdim>0) then 906 ABI_MALLOC(tdks%nhatgr,(cplex*tdks%nfftf,dtset%nspden,3*tdks%nhatgrdim)) 907 else 908 ABI_MALLOC(tdks%nhatgr,(0,0,0)) 909 end if 910 911 ABI_MALLOC(tdks%nhat,(tdks%nfftf,dtset%nspden*psps%usepaw)) 912 913 !Required in the PAW case to compute the inverse of the overlap (invovl) operator 914 call init_invovl(dtset%nkpt) 915 else 916 ABI_MALLOC(tdks%nhat,(0,0)) 917 ABI_MALLOC(tdks%nhatgr,(0,0,0)) 918 tdks%nhatgrdim=0 919 end if 920 921 !Allocate various required arrays for calculation of the Hamiltonian 922 !Potentials 923 ABI_MALLOC(tdks%vhartr,(tdks%nfftf)) 924 tdks%vhartr=zero 925 ABI_MALLOC(tdks%vpsp,(tdks%nfftf)) 926 tdks%vpsp=zero 927 ABI_MALLOC(tdks%vtrial,(tdks%nfftf,dtset%nspden)) 928 tdks%vtrial=zero 929 ABI_MALLOC(tdks%vxc,(tdks%nfftf,dtset%nspden)) 930 tdks%vxc=zero 931 if (psps%n1xccc/=0) then 932 ABI_MALLOC(tdks%xccc3d,(tdks%nfftf)) 933 else 934 ABI_MALLOC(tdks%xccc3d,(0)) 935 end if 936 tdks%xccc3d=zero 937 if (psps%usepaw==1) then 938 ABI_MALLOC(tdks%xcctau3d,(tdks%nfftf*dtset%usekden)) 939 tdks%xcctau3d=zero 940 endif 941 !For mGGA 942 ABI_MALLOC(tdks%vxctau,(tdks%nfftf,dtset%nspden,4*dtset%usekden)) 943 tdks%vxctau=zero 944 !For hybrid functionals 945 use_hybcomp=0 946 if(mod(dtset%fockoptmix,100)==11) use_hybcomp=1 947 ABI_MALLOC(tdks%vxc_hybcomp,(tdks%pawfgr%nfft,dtset%nspden*use_hybcomp)) 948 tdks%vxc_hybcomp=zero 949 !For VDW corrected functionals 950 tdks%ngrvdw=0 951 if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7)) then 952 tdks%ngrvdw=dtset%natom 953 end if 954 ABI_MALLOC(tdks%grvdw,(3,tdks%ngrvdw)) 955 tdks%grvdw=zero 956 957 !Compute large sphere G^2 cut-off (gsqcut) and box / sphere ratio 958 if (psps%usepaw==1) then 959 call getcut(tdks%boxcut,dtset%pawecutdg,tdks%gmet,tdks%gsqcut,dtset%iboxcut, & 960 & std_out,k0,tdks%pawfgr%ngfft) 961 else 962 call getcut(tdks%boxcut,dtset%ecut,tdks%gmet,tdks%gsqcut,dtset%iboxcut, & 963 & std_out,k0,tdks%pawfgr%ngfft) 964 end if 965 966 !Compute structure factor phases (exp(2Pi i G.xred)) on coarse and fine grid 967 ABI_MALLOC(tdks%ph1d,(2,3*(2*tdks%pawfgr%mgfftc+1)*dtset%natom)) 968 ABI_MALLOC(tdks%ph1df,(2,3*(2*tdks%pawfgr%mgfft+1)*dtset%natom)) 969 call getph(tdks%atindx,dtset%natom,tdks%pawfgr%ngfftc(1),tdks%pawfgr%ngfftc(2), & 970 & tdks%pawfgr%ngfftc(3),tdks%ph1d,tdks%xred) 971 if (psps%usepaw==1.and.tdks%pawfgr%usefinegrid==1) then 972 call getph(tdks%atindx,dtset%natom,tdks%pawfgr%ngfft(1),tdks%pawfgr%ngfft(2), & 973 & tdks%pawfgr%ngfft(3),tdks%ph1df,tdks%xred) 974 else 975 tdks%ph1df(:,:)=tdks%ph1d(:,:) 976 end if 977 978 !!FB: @MT Needed? If yes, then don't forget to put it back in the begining of 979 !! propagate_ele as well 980 !!if any nuclear dipoles are nonzero, compute the vector potential in real space (depends on 981 !!atomic position so should be done for nstep = 1 and for updated ion positions 982 !if ( any(abs(dtset%nucdipmom(:,:))>tol8) ) then 983 ! with_vectornd = 1 984 !else 985 ! with_vectornd = 0 986 !end if 987 !if(allocated(vectornd)) then 988 ! ABI_FREE(vectornd) 989 !end if 990 !ABI_MALLOC(vectornd,(with_vectornd*nfftf,3)) 991 !vectornd=zero 992 !if(with_vectornd .EQ. 1) then 993 ! call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,ngfftf,dtset%nucdipmom,& 994 ! & rprimd,vectornd,xred) 995 !endif 996 997 !Allocate memory for density 998 ABI_MALLOC(tdks%rhor,(tdks%nfftf,dtset%nspden)) 999 ABI_MALLOC(tdks%taur,(tdks%nfftf,dtset%nspden*dtset%usekden)) 1000 ABI_MALLOC(tdks%rhog,(2,tdks%nfftf)) 1001 ABI_MALLOC(tdks%taug,(2,tdks%nfftf*dtset%usekden)) 1002 1003 end subroutine second_setup
m_rttddft_tdks/tdks_free [ Functions ]
[ Top ] [ m_rttddft_tdks ] [ Functions ]
NAME
tdks_free
FUNCTION
Free all the memory associated with the tdks object
INPUTS
tdks <class(tdks_type)> = the tdks object to free dtset <type(dataset_type)> = all input variables for this dataset mpi_enreg <MPI_type> = MPI-parallelisation information psps <type(pseudopotential_type)> = variables related to pseudopotentials
OUTPUT
SOURCE
355 subroutine tdks_free(tdks,dtset,mpi_enreg,psps) 356 357 implicit none 358 359 !Arguments ------------------------------------ 360 !scalars 361 class(tdks_type), intent(inout) :: tdks 362 type(dataset_type), intent(inout) :: dtset 363 type(MPI_type), intent(inout) :: mpi_enreg 364 type(pseudopotential_type), intent(inout) :: psps 365 366 ! *********************************************************************** 367 368 !Destroy hidden save variables 369 call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg) 370 if (psps%usepaw ==1) then 371 call destroy_invovl(dtset%nkpt,dtset%gpu_option) 372 end if 373 if(tdks%gemm_nonlop_use_gemm) then 374 call destroy_gemm_nonlop(dtset%nkpt) 375 end if 376 377 !Call type destructors 378 call destroy_sc_dmft(tdks%paw_dmft) 379 call pawfgr_destroy(tdks%pawfgr) 380 call tdks%hdr%free() 381 382 !Nullify pointers 383 if(associated(tdks%pawang)) tdks%pawang => null() 384 if(associated(tdks%pawrad)) tdks%pawrad => null() 385 if(associated(tdks%pawtab)) tdks%pawtab => null() 386 if(associated(tdks%pawrhoij)) call pawrhoij_free(tdks%pawrhoij) 387 388 !Deallocate allocatables 389 ABI_SFREE(tdks%atindx) 390 ABI_SFREE(tdks%atindx1) 391 ABI_SFREE(tdks%cg) 392 ABI_SFREE(tdks%cg0) 393 ABI_SFREE(tdks%dimcprj) 394 ABI_SFREE(tdks%eigen) 395 ABI_SFREE(tdks%eigen0) 396 ABI_SFREE(tdks%grvdw) 397 ABI_SFREE(tdks%indsym) 398 ABI_SFREE(tdks%irrzon) 399 ABI_SFREE(tdks%kg) 400 ABI_SFREE(tdks%nattyp) 401 ABI_SFREE(tdks%nhat) 402 ABI_SFREE(tdks%nhatgr) 403 ABI_SFREE(tdks%npwarr) 404 ABI_SFREE(tdks%occ) 405 ABI_SFREE(tdks%occ0) 406 ABI_SFREE(tdks%ph1d) 407 ABI_SFREE(tdks%ph1df) 408 ABI_SFREE(tdks%phnons) 409 ABI_SFREE(tdks%rhog) 410 ABI_SFREE(tdks%rhor) 411 ABI_SFREE(tdks%symrec) 412 ABI_SFREE(tdks%taug) 413 ABI_SFREE(tdks%taur) 414 ABI_SFREE(tdks%vhartr) 415 ABI_SFREE(tdks%vpsp) 416 ABI_SFREE(tdks%vtrial) 417 ABI_SFREE(tdks%vxc) 418 ABI_SFREE(tdks%vxctau) 419 ABI_SFREE(tdks%vxc_hybcomp) 420 ABI_SFREE(tdks%xred) 421 ABI_SFREE(tdks%xccc3d) 422 ABI_SFREE(tdks%xcctau3d) 423 ABI_SFREE(tdks%ylm) 424 ABI_SFREE(tdks%ylmgr) 425 426 if(allocated(tdks%cprj)) then 427 call pawcprj_free(tdks%cprj) 428 ABI_FREE(tdks%cprj) 429 end if 430 if(allocated(tdks%cprj0)) then 431 call pawcprj_free(tdks%cprj0) 432 ABI_FREE(tdks%cprj0) 433 end if 434 if(allocated(tdks%paw_an)) then 435 call paw_an_free(tdks%paw_an) 436 ABI_FREE(tdks%paw_an) 437 end if 438 if(allocated(tdks%pawfgrtab)) then 439 call pawfgrtab_free(tdks%pawfgrtab) 440 ABI_FREE(tdks%pawfgrtab) 441 end if 442 if(allocated(tdks%paw_ij)) then 443 call paw_ij_free(tdks%paw_ij) 444 ABI_FREE(tdks%paw_ij) 445 end if 446 447 end subroutine tdks_free
m_rttddft_tdks/tdks_init [ Functions ]
[ Top ] [ m_rttddft_tdks ] [ Functions ]
NAME
tdks_init
FUNCTION
Initialize the tdks object
INPUTS
codvsn = code version dtfil <type datafiles_type> = infos about file names, file unit numbers dtset <type(dataset_type)> = all input variables for this dataset mpi_enreg <MPI_type> = MPI-parallelisation information pawang <type(pawang_type)> = paw angular mesh and related data pawrad(ntypat*usepaw) <type(pawrad_type)> = paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)> = paw tabulated starting data psps <type(pseudopotential_type)> = variables related to pseudopotentials
OUTPUT
tdks <class(tdks_type)> = the tdks object to initialize
SOURCE
213 subroutine tdks_init(tdks ,codvsn, dtfil, dtset, mpi_enreg, pawang, pawrad, pawtab, psps) 214 215 implicit none 216 217 !Arguments ------------------------------------ 218 !scalars 219 class(tdks_type), intent(inout) :: tdks 220 character(len=8), intent(in) :: codvsn 221 type(datafiles_type), intent(in) :: dtfil 222 type(dataset_type), intent(inout) :: dtset 223 type(MPI_type), intent(inout) :: mpi_enreg 224 type(pawang_type), intent(inout),target :: pawang 225 type(pseudopotential_type), intent(inout) :: psps 226 !arrays 227 type(pawrad_type), intent(inout),target :: pawrad(psps%ntypat*psps%usepaw) 228 type(pawtab_type), intent(inout),target :: pawtab(psps%ntypat*psps%usepaw) 229 230 !Local variables------------------------------- 231 !scalars 232 integer :: ierr 233 integer :: my_natom 234 integer :: ncpgr 235 integer :: psp_gencond 236 real(dp) :: ecut_eff 237 character(len=500) :: msg 238 character(len=fnlen) :: fname_wfk 239 type(extfpmd_type),pointer :: extfpmd => null() 240 !arrays 241 real(dp),allocatable :: doccde(:) 242 243 ! *********************************************************************** 244 245 my_natom=mpi_enreg%my_natom 246 247 !1) Various initializations & checks (MPI, PW, FFT, PSP, Symmetry ...) 248 call first_setup(codvsn,dtfil,dtset,ecut_eff,mpi_enreg,pawrad,pawtab,psps,psp_gencond,tdks) 249 250 !2) Deals with restart 251 !FB: @MT Is this the proper way to read a file in abinit..? 252 tdks%first_step = 1 253 tdks%fname_tdener = dtfil%fnameabo_td_ener 254 tdks%fname_wfk0 = dtfil%fnamewffk 255 fname_wfk = dtfil%fnamewffk 256 if (dtset%td_restart > 0) then 257 if (mpi_enreg%me == 0) then 258 if (open_file('TD_RESTART', msg, newunit=tdks%tdrestart_unit, status='old', form='formatted') /= 0) then 259 write(msg,'(a,a,a)') 'Error while trying to open file TD_RESTART needed to restart the calculation.' 260 ABI_ERROR(msg) 261 end if 262 read(tdks%tdrestart_unit,*) tdks%first_step 263 tdks%first_step = tdks%first_step + 1 264 read(tdks%tdrestart_unit,*) tdks%fname_tdener 265 read(tdks%tdrestart_unit,*) tdks%fname_wfk0 266 read(tdks%tdrestart_unit,*) fname_wfk 267 end if 268 !Send to all procs 269 call xmpi_bcast(tdks%first_step,0,mpi_enreg%comm_world,ierr) 270 call xmpi_bcast(tdks%fname_tdener,0,mpi_enreg%comm_world,ierr) 271 call xmpi_bcast(tdks%fname_wfk0,0,mpi_enreg%comm_world,ierr) 272 call xmpi_bcast(fname_wfk,0,mpi_enreg%comm_world,ierr) 273 else 274 if (mpi_enreg%me == 0) then 275 if (open_file('TD_RESTART', msg, newunit=tdks%tdrestart_unit, status='unknown', form='formatted') /= 0) then 276 write(msg,'(a,a,a)') 'Error while trying to open file TD_RESTART.' 277 ABI_ERROR(msg) 278 end if 279 end if 280 end if 281 282 !3) Reads initial KS orbitals from file (calls inwffil) 283 call read_wfk(dtfil,dtset,ecut_eff,fname_wfk,mpi_enreg,tdks) 284 285 !4) Init occupation numbers 286 ABI_MALLOC(tdks%occ0,(dtset%mband*dtset%nkpt*dtset%nsppol)) 287 tdks%occ0(:)=dtset%occ_orig(:,1) 288 !calc occupation number with metallic occupation using the previously read WF 289 if (dtset%occopt>=3.and.dtset%occopt<=9) then ! allowing for occopt 9 290 ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 291 call newocc(doccde,tdks%eigen0,tdks%energies%entropy,tdks%energies%e_fermie, & 292 & tdks%energies%e_fermih,dtset%ivalence,dtset%spinmagntarget, & 293 & dtset%mband,dtset%nband,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD, & 294 & dtset%nkpt,dtset%nspinor,dtset%nsppol,tdks%occ0,dtset%occopt, & 295 & dtset%prtvol,dtset%tphysel,dtset%tsmear,dtset%wtk,extfpmd=extfpmd) 296 ABI_FREE(doccde) 297 end if 298 299 !5) Some further initialization (Mainly for PAW) 300 call second_setup(dtset,mpi_enreg,pawang,pawrad,pawtab,psps,psp_gencond,tdks) 301 302 !Keep initial wavefunction in memory 303 if (dtset%td_restart == 0) then 304 ABI_MALLOC(tdks%cg0,(2,tdks%mcg)) 305 tdks%cg0(:,:) = tdks%cg(:,:) 306 end if 307 !and associated cprojs to compute occupations 308 if (psps%usepaw ==1) then 309 ncpgr=0 310 ABI_MALLOC(tdks%cprj0,(dtset%natom,tdks%mcprj)) 311 call pawcprj_alloc(tdks%cprj0,ncpgr,tdks%dimcprj) 312 call ctocprj(tdks%atindx,tdks%cg0,1,tdks%cprj0,tdks%gmet,tdks%gprimd,0,0,0, & 313 & dtset%istwfk,tdks%kg,dtset%kptns,tdks%mcg,tdks%mcprj,dtset%mgfft, & 314 & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,tdks%nattyp, & 315 & dtset%nband,dtset%natom,dtset%ngfft,dtset%nkpt,dtset%nloalg, & 316 & tdks%npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,psps%ntypat, & 317 & dtset%paral_kgb,tdks%ph1d,psps,tdks%rmet,dtset%typat,tdks%ucvol, & 318 & tdks%unpaw,tdks%xred,tdks%ylm,tdks%ylmgr) 319 end if 320 ABI_MALLOC(tdks%occ,(dtset%mband*dtset%nkpt*dtset%nsppol)) 321 322 !FB: That should be all for now but there were few more initializations in 323 !g_state.F90 in particular related to electric field, might want to check it out 324 !once we reach the point of including external electric field 325 326 !Keep some additional stuff in memory within the tdks object 327 tdks%unpaw = dtfil%unpaw 328 tdks%dt = dtset%dtele 329 tdks%ntime = dtset%ntime 330 331 tdks%pawang => pawang 332 tdks%pawrad => pawrad 333 tdks%pawtab => pawtab 334 335 end subroutine tdks_init