!!---------- cp2kmate is a collection of utilities of CP2K
subroutine cp2kmate
use defvar
use util
character c200tmp*200

write(*,*)
call menutitle("Auxiliary tools for CP2K (CP2Kmate)",10,1)
write(*,*) "0 Return"
write(*,*) "1 Create CP2K input file"
write(*,*) "2 Convert band structure file (.bs) to plain text file and then plot"
write(*,*) "3 Plot band structure map based on the files generated by option 2"
write(*,"(a)") " 4 Obtain LUCO and HOCO energies and show band gap based on the output file containing energy levels of every k-point"
write(*,*) "5 Obtain exact DOS based on energy levels at all k-points"
write(*,*) "6 Calculate and export overlap matrix to CP2K_overlap.txt in current folder"
write(*,*) "7 Load orbital energies from CP2K output file"
write(*,*) "8 Load orbitals of specific k-point from CP2K output file"
write(*,*) "9 Load overlap matrix from .csr file exported by CP2K"
read(*,*) isel

if (isel==1) then
    call outCP2Kinp_wrapper
else if (isel==2) then
    call CP2K_BS
else if (isel==3) then
    call plotBS
else if (isel==4) then
    call CP2K_bandgap_DOS(1)
else if (isel==5) then
    call CP2K_bandgap_DOS(2)
else if (isel==6) then
    call ask_Sbas_PBC
    write(*,*) "Exporting..."
    open(10,file="CP2K_overlap.txt",status="replace")
	call showmatgau(Sbas,"Overlap matrix",1,fileid=10)
    close(10)
    write(*,*) "Done! Overlap matrix has been exported to CP2K_overlap.txt in current folder"
else if (isel==7) then
    call CP2K_MOene_load
else if (isel==8) then
    call CP2K_loadkpwfn
else if (isel==9) then
    c200tmp=" "
    call CP2K_loadSbas(c200tmp)
end if
end subroutine




!!---------- Interface of outputting CP2K input file
subroutine outCP2Kinp_wrapper
use util
use defvar
character(len=200) outname,c200tmp
call path2filename(filename,c200tmp)
write(*,"(/,a)") " Note: Please cite original papers of Multiwfn if you benefits from this function in your study!"
write(*,*) 
write(*,*) "Input path for generating CP2K input file, e.g. C:\ltwd.inp"
write(*,"(a)") " If press ENTER button directly, will export to "//trim(c200tmp)//".inp"
read(*,"(a)") outname
if (outname==" ") outname=trim(c200tmp)//".inp"
call outCP2Kinp(outname,10)
end subroutine
!!---------- Output current coordinate to CP2K input file
!NOTE: To determine isolated, use (ifPBC==0.or.PBCdir=="NONE"). ifPBC==0 doesn't necessarily mean the system is isolated, because .restart of CP2K carries cell information
!Use PBCdir=="NONE" to determine only when the options are relevant to calculation of electrostatic interaction
subroutine outCP2Kinp(outname,ifileid)
use defvar
use util
implicit real*8 (a-h,o-z)
character(len=*) outname
integer :: ifileid,ibas=2,tmparr(ncenter)
character selectyn,c80tmp*80,c80tmp2*80,c200tmp*200,c2000tmp*2000
character :: method*22="PBE",PBCdir*4="XYZ ",cellfix*4="NONE"
character(len=30) :: basname(-10:30)=" "
integer :: itask=1,idispcorr=0,imolden=0,ioutvibmol=1,ithermostat=0,ibarostat=0,inoSCFinfo=0,iSCCS=0,idipcorr=0,iwfc=0,iHFX=0,iRIHFX=0,imoment=0,ihyperfine=0,ioptmethod=1,iprintlevel=1
integer :: iTDDFT=0,nstates_TD=3,iTDtriplet=0,isTDA=0,iNTO=0,nADDED_MOS=0,icentering=0,itightopt=0,iraman=0,iSOCTDDFT=0
integer :: iMDformat=1,nMDsavefreq=1,ioutcube=0,idiagOT=1,imixing=2,ismear=0,iatomcharge=0,ifineXCgrid=0,iouterSCF=1,iDFTplusU=0,NHOMO=0,NLUMO=0
integer :: natmcons=0,nthermoatm=0,ikpoint1=1,ikpoint2=1,ikpoint3=1,nrep1=1,nrep2=1,nrep3=1,ikeepcell=0
integer,allocatable :: atmcons(:),thermoatm(:)
real*8 :: efieldvec(3)=0,vacsizex=5/b2a,vacsizey=5/b2a,vacsizez=5/b2a
real*8 :: frag1chg,frag2chg
integer :: frag1multi,frag2multi,totalmulti
integer :: iprestype=1,ioutSbas=0,ioutKSbas=0,ioutorbene=0,istate_force=1,idiaglib=1,iGAPW=0,iLSSCF=0,iLRIGPW=0,iPSOLVER=1,niter_evGW=1,niter_scGW0=1,istructfile=0
real*8 :: Piso=1.01325D0,Ptens(3,3)=reshape( [1.01325D0,0D0,0D0, 0D0,1.01325D0,0D0, 0D0,0D0,1.01325D0], shape=shape(Ptens))
real*8 :: PBEh_HFX=45
integer :: CUTOFF=350,REL_CUTOFF=50
real*8 :: cellv1_pseudo(3),cellv2_pseudo(3),cellv3_pseudo(3) !Pseudo cell for low-dimensional system
real*8 :: cellv1_tmp(3),cellv2_tmp(3),cellv3_tmp(3)
!CDFT related arrays
integer :: nCDFTgroup=0 !Number of CDFT groups
integer,allocatable,save :: CDFTatm(:,:),CDFTtype(:),CDFTnatm(:) !Atom indices, constraint type and number of involved atoms of CDFT groups
real*8,allocatable,save :: CDFTtarget(:) !Target value of CDFT groups
!XAS related
integer :: iXAS_SOC=0,iGW2X=0
integer,allocatable,save :: XASatm(:)
character,save :: kindnameXAS*8

!Status information of current system. ",save" is used so that for the same system we can enter this interface multiple times to generate various input files
integer,save :: netchg,multispin
integer,allocatable,save :: atmkind(:) !The kind that atoms belonging to
integer,parameter :: nkindmax=200
integer,save :: nkind=0 !Current number of kinds
character(len=8),save :: kindname(nkindmax) !Name of each kind
integer,save :: kindeleidx(nkindmax) !Element idx of each kind
real*8,save :: kindmag(nkindmax) !Magnetization of each kind
character,save :: lastinpname*200 !Input file of last time in this interface

!Defining array recording number of valence electrons for various elements of MOLOPT-GTH basis set, so that -q? can be automatically added to basis set name (except for 0)
!For some elements, there are small and large core versions. Here records the default one, which corresponds to large core for transition metals, and small core for others
!For CP2K 9.1, the information comes from GTH-PADE of GTH_POTENTIALS, except that Ln and Ac series come from LnPP1_POTENTIALS and AcPP1_POTENTIALS (medium-core version), respectively
integer :: Nval(0:nelesupp)=(/ 0, & !X
1, 2,                                             & !H ~He
3, 4,                                3,4,5,6,7,8, & !Li~Ne
9,10,                                3,4,5,6,7,8, & !Na~Ar
9,10,11,12,13,14,15,16,17,18,11,12, 13,4,5,6,7,8, & !K ~Kr
9,10,11,12,13,14,15,16,17,18,11,12, 13,4,5,6,7,8, & !Rb~Xe
9,10,11,12,13,14,15,16,17,18,29,30,31,32,33,34,35,& !Cs,Ba,La~Lu
        12,13,14,15,16,17,18,11,12, 13,4,5,6,7,8, & !Hf~Rn
0, 0,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,& !Fr,Ra,Ac~Lr
(0,iele=104,nelesupp) /)

iconvtest=0
ikpconvtest=0
natmcons=0 !Disable atom constraint setting if set last time
if (index(outname,"cutconv.inp")/=0) then
    iconvtest=1
    iprintlevel=2 !Medium printing level
    !ifineXCgrid=1 !I found this effectively often improves smoothness of convergence with respect to cutoff
    write(*,*) "Note: The generated file is dedicated to cutoff convergence testing purpose"
else if (index(outname,"kpconv.inp")/=0) then
    ikpconvtest=1
    write(*,*) "Note: The generated file is dedicated to k-point convergence testing purpose"
end if

!Conversion of basis set index to basis set name
basname(-1)="SZV-GTH"
basname(-2)="DZVP-GTH"
basname(-3)="TZVP-GTH"
basname(-4)="TZV2P-GTH"
basname(-5)="QZV2P-GTH"
basname(-6)="QZV3P-GTH"
basname(1)="SZV-MOLOPT-SR-GTH"
basname(2)="DZVP-MOLOPT-SR-GTH"
basname(3)="TZVP-MOLOPT-GTH"
basname(4)="TZV2P-MOLOPT-GTH"
basname(5)="TZV2PX-MOLOPT-GTH"
basname(7)="ccGRB-D"
basname(8)="ccGRB-T"
basname(9)="ccGRB-Q"
basname(10)="6-31G*"
basname(11)="6-311G**"
basname(12)="Ahlrichs-def2-TZVP"
basname(13)="pob-TZVP"
basname(14)="pob-DZVP-rev2"
basname(15)="pob-TZVP-rev2"
basname(16)="Ahlrichs-def2-QZVP"
basname(19)="def2-QZVP with RI_5Z"
basname(20)="cc-DZ with RI_DZ"
basname(21)="cc-TZ with RI_TZ"
basname(22)="cc-QZ with RI_QZ"

!Construct kind information, charge and multiplicity
10 if (allocated(atmkind).and.filename/=lastinpname) deallocate(atmkind) !Has allocated last time, need to regenerate if the last system is different to the present one
if (.not.allocated(atmkind)) then !Haven't been constructed, namely first time use this function for present system, or this system it not identical to the last one
    write(*,*) "Generating KIND information..."
    allocate(atmkind(ncenter))
    !Build initial kind list
    nkind=1
    kindeleidx(1)=a(1)%index
    kindname(1)=a(1)%name
    do iatm=2,ncenter
        if (all(kindeleidx(:nkind)/=a(iatm)%index)) then
            nkind=nkind+1
            kindeleidx(nkind)=a(iatm)%index
            kindname(nkind)=a(iatm)%name
        end if
    end do
    !Assign each atom by a kind
    do iatm=1,ncenter
        do ikind=1,nkind
            if (a(iatm)%index==kindeleidx(ikind)) then
                atmkind(iatm)=ikind
                exit
            end if
        end do
    end do
    kindmag=0
    netchg=sum(a%charge)-nint(nelec)
    multispin=nint(naelec-nbelec)+1
    if (ncenter>500) ioptmethod=2 !LBFGS is more suitable for large system than BFGS
    lastinpname=filename
end if

do while(.true.)
    !do iatm=1,ncenter
    !    ikind=atmkind(iatm)
    !    write(*,"(' Atom',i6,1x,a,'  Kind index:',i3,'  Kind name: ',a,'  Kind elem idx:',i3)") iatm,a(iatm)%name,ikind,kindname(ikind),kindeleidx(ikind)
    !end do
    write(*,*)
    write(*,*) "-11 Enter the interface for geometry operations"
    write(*,*) "-10 Return"
    write(*,*) "-9 Other settings"
    write(*,*) "-7 Set direction(s) of applying periodic boundary condition, current: "//PBCdir
    if (itask==3.or.itask==4.or.itask==6.or.itask==7) then !Optimizing minimum, MD, TS search
        if (itask==6) write(*,"(a,i6)") " -6 Set frequency of writing molecular dynamics trajectory, current:",nMDsavefreq
        if (iMDformat==1) write(*,*) "-5 Choose format of outputted trajectory, current: xyz"
        if (iMDformat==2) write(*,*) "-5 Choose format of outputted trajectory, current: dcd"
        if (iMDformat==-2) write(*,*) "-5 Choose format of outputted trajectory, current: dcd_aligned_cell"
        if (iMDformat==3) write(*,*) "-5 Choose format of outputted trajectory, current: pdb"
    end if
    if (method/="FIST") then
        if (iatomcharge==0) write(*,*) "-4 Calculate atomic charges, current: None"
        if (iatomcharge==1) write(*,*) "-4 Calculate atomic charges, current: Mulliken"
        if (iatomcharge==2) write(*,*) "-4 Calculate atomic charges, current: Lowdin"
        if (iatomcharge==3) write(*,*) "-4 Calculate atomic charges, current: Hirshfeld"
        if (iatomcharge==4) write(*,*) "-4 Calculate atomic charges, current: Hirshfeld-I"
        if (iatomcharge==5) write(*,*) "-4 Calculate atomic charges, current: Voronoi"
        if (iatomcharge==6) write(*,*) "-4 Calculate atomic charges, current: RESP"
        if (iatomcharge==7) write(*,*) "-4 Calculate atomic charges, current: REPEAT"
    end if
    if (ioutcube==0) write(*,*) "-3 Set exporting cube file, current: None"
    if (ioutcube==1) write(*,*) "-3 Set exporting cube file, current: Electron density"
    if (ioutcube==2) write(*,*) "-3 Set exporting cube file, current: ELF"
    if (ioutcube==3) write(*,*) "-3 Set exporting cube file, current: XC potential"
    if (ioutcube==4) write(*,*) "-3 Set exporting cube file, current: Hartree potential (negative of ESP)"
    if (ioutcube==5) write(*,*) "-3 Set exporting cube file, current: Electric field"
    if (ioutcube==6) write(*,"(a,i6,a,i6)") " -3 Set exporting cube file, current: MOs, with NHOMO=",NHOMO,", NLUMO=",NLUMO
    if (ioutcube==7) write(*,*) "-3 Set exporting cube file, current: Electron density + Hartree pot."
    if (imolden==0) write(*,*) "-2 Toggle exporting .molden file for Multiwfn, current: No"
    if (imolden==1) write(*,*) "-2 Toggle exporting .molden file for Multiwfn, current: Yes"
    if (itask==1) write(*,*) "-1 Choose task, current: Energy"
    if (itask==2) write(*,*) "-1 Choose task, current: Energy + force"
    if (itask==3) write(*,*) "-1 Choose task, current: Optimizing structure"
    if (itask==4) write(*,*) "-1 Choose task, current: Optimizing structure and cell"
    if (itask==5) write(*,*) "-1 Choose task, current: Vibrational analysis"
    if (itask==6) write(*,*) "-1 Choose task, current: Molecular dynamics"
    if (itask==7) write(*,*) "-1 Choose task, current: Searching transition state by dimer algorithm"
    if (itask==8) write(*,*) "-1 Choose task, current: BAND"
    if (itask==9) write(*,*) "-1 Choose task, current: NMR"
    if (itask==10) write(*,*) "-1 Choose task, current: Polarizability"
    if (itask==11) write(*,*) "-1 Choose task, current: BSSE"
    if (itask==13) write(*,*) "-1 Choose task, current: Real-time propagation for electron dynamics"
    if (itask==14) write(*,*) "-1 Choose task, current: Path-integral molecular dynamics (PIMD)"
    if (itask==15) write(*,*) "-1 Choose task, current: X-ray absorption spectroscopy"
    write(*,*) " 0 Generate input file now!"
    c80tmp=" "
    if (iHFX==1) then
        if (iRIHFX==1) then
            c80tmp=" with RI-HFX"
        else
            c80tmp=" without RI-HFX"
        end if
    end if
    if (index(method,"PBEh")==0) then
        write(*,*) " 1 Choose theoretical method, current: "//trim(method)//trim(c80tmp)
    else
        write(*,"(a,f7.3,' %',a)") "  1 Choose theoretical method, current: "//trim(method)//" with HFX of",PBEh_HFX,trim(c80tmp)
    end if
    if (method(1:3)/="GFN".and.method/="PM6".and.method/="SCC-DFTB".and.method/="FIST") write(*,*) " 2 Choose basis set and pseudopotential, current: "//trim(basname(ibas))
    if (index(method,"MP2")==0.and.index(method,"RPA")==0.and.index(method,"GW")==0.and.method(1:3)/="GFN".and.method/="PM6".and.method/="SCC-DFTB".and.method/="BEEFVDW".and.method/="FIST") then
        if (idispcorr==0) write(*,*) " 3 Set dispersion correction, current: None"
        if (idispcorr==1) write(*,*) " 3 Set dispersion correction, current: DFT-D3"
        if (idispcorr==2) write(*,*) " 3 Set dispersion correction, current: DFT-D3(BJ)"
        if (idispcorr==3) write(*,*) " 3 Set dispersion correction, current: DFT-D4"
        if (idispcorr==5) write(*,*) " 3 Set dispersion correction, current: rVV10"
    end if
    if (method/="FIST") then
        if (iLSSCF==0) then
            if (idiagOT==1) then
                write(*,*) " 4 Switching between diagonalization and OT, current: Diagonalization"
                if (method/="PM6") then !I found PM6 can only use Direct mixing
                    if (imixing==1) write(*,*) " 5 Set density matrix mixing, current: Direct mixing + DIIS"
                    if (imixing==2) write(*,*) " 5 Set density matrix mixing, current: Broyden mixing"
                    if (imixing==3) write(*,*) " 5 Set density matrix mixing, current: Pulay mixing"
                end if
                if (ismear==0) write(*,*) " 6 Toggle smearing electron occupation, current: No"
                if (ismear==1) write(*,*) " 6 Toggle smearing electron occupation, current: Yes"
            else if (idiagOT==2) then
                write(*,*) " 4 Switching between diagonalization and OT, current: OT"
                if (iouterSCF==0) write(*,*) " 5 Toggle using outer SCF process, current: No"
                if (iouterSCF==1) write(*,*) " 5 Toggle using outer SCF process, current: Yes"
            end if
        end if
        if (iSCCS==0) write(*,*) " 7 Toggle using self-consistent continuum solvation (SCCS), current: No"
        if (iSCCS==1) write(*,*) " 7 Toggle using self-consistent continuum solvation (SCCS), current: Yes"
        if (PBCdir/="NONE") then
            if (ikpoint1==1.and.ikpoint2==1.and.ikpoint3==1) then
                write(*,*) " 8 Set k-points, current: GAMMA only"
            else
                write(*,"(a,3i3)") "  8 Set k-points, current: MONKHORST-PACK",ikpoint1,ikpoint2,ikpoint3
            end if
        end if
    end if
    !!!! Below are task specific options
    if (itask>=3.and.itask<=8) then
        if (natmcons==0) then
            write(*,*) " 9 Set atom position freeze, current: None"
        else
            write(*,"(a,i6)") "  9 Set atom position freeze, current:",natmcons
        end if
    end if
    if (itask==6) then
        if (ithermostat==0) write(*,*) "10 Set thermostat, current: None"
        if (ithermostat==1) write(*,*) "10 Set thermostat, current: Adaptive-Langevin"
        if (ithermostat==2) write(*,*) "10 Set thermostat, current: Canonical sampling through velocity rescaling"
        if (ithermostat==3) write(*,*) "10 Set thermostat, current: Generalized Langevin Equation (GLE)"
        if (ithermostat==4) write(*,*) "10 Set thermostat, current: Nose-Hoover"
        if (ithermostat>0) then
            !if (nthermoatm==ncenter) then !This option is fully misleading!
            !    write(*,*) "11 Set region for the thermostat, current: All atoms"
            !else
            !    write(*,"(a,i6,' atoms')") " 11 Set region for the thermostat, number of current atoms:",nthermoatm
            !end if
        end if
        if (ibarostat==0) write(*,*) "12 Set barostat, current: None"
        if (ibarostat==1) write(*,*) "12 Set barostat, current: Yes, flexible cell"
        if (ibarostat==2) write(*,*) "12 Set barostat, current: Yes, isotropic cell"
        if (method/="FIST") then
            if (inoSCFinfo==0) write(*,*) "13 Toggle suppressing printing SCF information during MD, current: No"
            if (inoSCFinfo==1) write(*,*) "13 Toggle suppressing printing SCF information during MD, current: Yes"
        end if
    else if (itask==3.or.itask==4) then
        if (ioptmethod==1) write(*,*) "10 Set optimization method, current: BFGS"
        if (ioptmethod==2) write(*,*) "10 Set optimization method, current: LBFGS"
        if (ioptmethod==3) write(*,*) "10 Set optimization method, current: CG"
    else if (itask==5) then
        if (ioutvibmol==0) write(*,*) "10 Toggle exporting Molden file recording vibrational modes, current: No"
        if (ioutvibmol==1) write(*,*) "10 Toggle exporting Molden file recording vibrational modes, current: Yes"
        if (PBCdir/="XYZ") then
            if (ikeepcell==0) write(*,*) "11 Toggle keeping cell and coordinates in input file unchanged, current: No"
            if (ikeepcell==1) write(*,*) "11 Toggle keeping cell and coordinates in input file unchanged, current: Yes"
        end if
    else if (itask==15) then
        if (iXAS_SOC==0) write(*,*) "10 Toggle considering spin-orbit coupling, current: No"
        if (iXAS_SOC==1) write(*,*) "10 Toggle considering spin-orbit coupling, current: Yes"
        if (iGW2X==0) write(*,*) "11 Toggle using GW2X correction, current: No"
        if (iGW2X==1) write(*,*) "11 Toggle using GW2X correction, current: Yes"
    end if
    if (itask==4) then
        write(*,*) "11 Set constraint of cell length(s) during optimization, current: "//trim(cellfix)
        if (iprestype==1) write(*,"(a,1PE13.5,' bar')") " 12 Set external pressure, current:",Piso
        if (iprestype==2) write(*,"(a)") " 12 Set external pressure, current: Anisotropic"
    end if
    if (itask==2.or.itask==3.or.itask==4.or.itask==7) then
        if (itightopt==0) write(*,"(a)") " 13 Toggle using tighter threshold (for subsequent freq), current: No"
        if (itightopt==1) write(*,"(a)") " 13 Toggle using tighter threshold (for subsequent freq), current: Yes"
    else if (itask==5) then
        if (iraman==0) write(*,"(a)") " 13 Toggle calculating Raman activities, current: No"
        if (iraman==1) write(*,"(a)") " 13 Toggle calculating Raman activities, current: Yes"
    end if
    if (method/="FIST".and.itask/=15) then
        if (iTDDFT==0) then
            write(*,*) "15 Toggle calculating excited states via TDDFT, current: No"
        else if (iTDDFT==1) then
            write(*,*) "15 Toggle calculating excited states via TDDFT, current: Yes"
            write(*,"(' 16 Set number of excited states to solve by TDDFT, current:',i5)") nstates_TD
            if (.not.(multispin>1.or.any(kindmag(1:nkind)/=0))) then !Current is closed-shell
                if (iTDtriplet==0) write(*,*) "17 Toggle spin of the excited states to be calculated, current: Singlet"
                if (iTDtriplet==1) write(*,*) "17 Toggle spin of the excited states to be calculated, current: Triplet"
            end if
            if (isTDA==0) write(*,*) "18 Toggle using sTDA approximation in TDDFT, current: No"
            if (isTDA==1) write(*,*) "18 Toggle using sTDA approximation in TDDFT, current: Yes"
            !At least for CP2K 8.1, 9.1, the occupied NTOs are wrong, so do not let user know this feature
            !if (iNTO==0) write(*,*) "19 Toggle if performing NTO analysis, current: No"
            !if (iNTO==1) write(*,*) "19 Toggle if performing NTO analysis, current: Yes"
            if (itask==2.or.itask==3.or.itask==4.or.itask==5.or.itask==6.or.itask==7.or.itask==8) then !Need force
                write(*,"(a,i5)") " 20 Choose the state to evaluate force, current:",istate_force
            end if
            if (iSOCTDDFT==0) write(*,*) "21 Toggle considering spin-orbit coupling effect in TDDFT, current: No"
            if (iSOCTDDFT==1) write(*,*) "21 Toggle considering spin-orbit coupling effect in TDDFT, current: Yes"
        end if
    end if
    read(*,*) isel
    
    if (isel==-11) then
        natmold=ncenter
        call geom_operation
        !User may use geometry operation interface to reorder atoms or construct supercell, so we need to rebuild atom kind information
        if (any(a(1:size(a_org))%index/=a_org%index).or.ncenter/=natmold) then
            deallocate(atmkind)
            if (allocated(atmcons)) deallocate(atmcons)
            natmcons=0
            !if (allocated(thermoatm)) deallocate(thermoatm)
            !nthermoatm=0
            ithermostat=0
            goto 10 
        end if
    else if (isel==-10) then
        return
    else if (isel==-9) then
        do while(.true.)
            write(*,*)
            write(*,*) "                     ---------- Other settings ----------"
            write(*,*) "0 Return"
            write(*,"(a,i5)") " 1 Set net charge, current:",netchg
            write(*,"(a,i5)") " 2 Set spin multiplicity, current:",multispin
            write(*,"(a,3i3)") " 3 Set number of repetitions of the cell in X, Y, Z, current:",nrep1,nrep2,nrep3
            if (ifineXCgrid==0) write(*,"(a)") " 4 Toggle using finer grid for exchange-correlation part, current: No"
            if (ifineXCgrid==1) write(*,"(a)") " 4 Toggle using finer grid for exchange-correlation part, current: Yes"
            write(*,"(a,i5,' and',i4,' Ry')") " 5 Set CUTOFF and REL_CUTOFF, current:",CUTOFF,REL_CUTOFF
            if (imoment==0) write(*,*) "6 Print electric/magnetic moments, current: No"
            if (imoment==1) write(*,*) "6 Print electric/magnetic moments, current: Yes"
            if (PBCdir=="NONE") write(*,"(a,3f8.3,' A')") " 7 Set vacuum size in both sides of X,Y,Z, current:",vacsizex*b2a,vacsizey*b2a,vacsizez*b2a
            if (PBCdir=="X") write(*,"(a,2f8.3,' A')") " 7 Set vacuum size in both sides of Y and Z, current:",vacsizey*b2a,vacsizez*b2a
            if (PBCdir=="Y") write(*,"(a,2f8.3,' A')") " 7 Set vacuum size in both sides of X and Z, current:",vacsizex*b2a,vacsizez*b2a
            if (PBCdir=="Z") write(*,"(a,2f8.3,' A')") " 7 Set vacuum size in both sides of X and Y, current:",vacsizex*b2a,vacsizey*b2a
            if (PBCdir=="XY") write(*,"(a,f8.3,' A')") " 7 Set vacuum size in both sides of Z, current:",vacsizez*b2a
            if (PBCdir=="XZ") write(*,"(a,f8.3,' A')") " 7 Set vacuum size in both sides of Y, current:",vacsizey*b2a
            if (PBCdir=="YZ") write(*,"(a,f8.3,' A')") " 7 Set vacuum size in both sides of X, current:",vacsizex*b2a
            if (iDFTplusU==0) write(*,*) "8 Toggle using DFT+U, current: No"
            if (iDFTplusU==1) write(*,*) "8 Toggle using DFT+U, current: Yes"
            if (all(kindmag(1:nkind)==0)) then
                write(*,*) "9 Define atomic magnetization in initial guess"
            else
                nmagset=count(kindmag(1:nkind)/=0)
                write(*,"(a,i3,a)") " 9 Redefine atomic magnetization in initial guess, current: Manually defined",nmagset," kinds"
            end if
            if (iprintlevel==0) write(*,*) "10 Choose printing level of output information, current: Silent"
            if (iprintlevel==1) write(*,*) "10 Choose printing level of output information, current: Low"
            if (iprintlevel==2) write(*,*) "10 Choose printing level of output information, current: Medium"
            if (iprintlevel==3) write(*,*) "10 Choose printing level of output information, current: High"
            if (all(efieldvec==0)) then
                write(*,"(a)") " 11 Set external electric field vector"
            else
                write(*,"(a,3f8.5,' a.u.')") " 11 Set external electric field vector, current:",efieldvec
            end if
            if (nADDED_MOS==-1) then
                write(*,*) "12 Set number of virtual orbitals to solve, current: All"
            else
                write(*,"(a,i6)") " 12 Set number of virtual orbitals to solve, current:",nADDED_MOS
            end if
            if (icentering==0) write(*,*) "13 Toggle centering the coordinates of the system in the box, current: No"
            if (icentering==1) write(*,*) "13 Toggle centering the coordinates of the system in the box, current: Yes"
            if (ioutorbene==0) write(*,*) "14 Toggle printing orbital energies and occupancies after SCF, current: No"
            if (ioutorbene==1) write(*,*) "14 Toggle printing orbital energies and occupancies after SCF, current: Yes"
            if (ioutSbas==0) write(*,*) "15 Toggle outputting overlap matrix to .csr file, current: No"
            if (ioutSbas==1) write(*,*) "15 Toggle outputting overlap matrix to .csr file, current: Yes"
            if (ioutKSbas==0) write(*,*) "16 Toggle outputting Kohn-Sham matrix to .csr file, current: No"
            if (ioutKSbas==1) write(*,*) "16 Toggle outputting Kohn-Sham matrix to .csr file, current: Yes"
            !if (idiaglib==1) write(*,*) "20 Choose diagonalization library, current: Default"
            !if (idiaglib==2) write(*,*) "20 Choose diagonalization library, current: ELPA"
            !if (idiaglib==3) write(*,*) "20 Choose diagonalization library, current: Scalapack"
            if (ihyperfine==0) write(*,*) "18 Toggle printing EPR hyperfine coupling tensor, current: No"
            if (ihyperfine==1) write(*,*) "18 Toggle printing EPR hyperfine coupling tensor, current: Yes"
            if (nCDFTgroup>0) then
                write(*,"(a,i3,a)") " 19 Redefine constrained DFT (CDFT) groups, current:",nCDFTgroup," groups"
            else
                write(*,"(a)") " 19 Enable constrained DFT (CDFT) and define groups"
            end if
            if (iLRIGPW==0) write(*,*) "20 Toggle using LRIGPW instead of GPW to accelerate calculation, current: No"
            if (iLRIGPW==1) write(*,*) "20 Toggle using LRIGPW instead of GPW to accelerate calculation, current: Yes"
            if (iLSSCF==0) write(*,*) "21 Toggle using Linear Scaling Self Consistent Field Method, current: No"
            if (iLSSCF==1) write(*,*) "21 Toggle using Linear Scaling Self Consistent Field Method, current: Yes"
            if (PBCdir=="XYZ") then
                if (iPSOLVER==1) write(*,*) "22 Set Poisson solver, current: PERIODIC"
                if (iPSOLVER==2) write(*,*) "22 Set Poisson solver, current: ANALYTIC"
                if (iPSOLVER==3) write(*,*) "22 Set Poisson solver, current: MT"
                if (iPSOLVER==4) write(*,*) "22 Set Poisson solver, current: WAVELET"
            else
                if (iPSOLVER==1) write(*,*) "22 Set Poisson solver and automatically set vacuum size, current: PERIODIC"
                if (iPSOLVER==2) write(*,*) "22 Set Poisson solver and automatically set vacuum size, current: ANALYTIC"
                if (iPSOLVER==3) write(*,*) "22 Set Poisson solver and automatically set vacuum size, current: MT"
                if (iPSOLVER==4) write(*,*) "22 Set Poisson solver and automatically set vacuum size, current: WAVELET"
            end if
            if (PBCdir=="XYZ".and.iPSOLVER==1) then
                if (idipcorr==0) write(*,*) "23 Set surface dipole correction, current: None"
                if (idipcorr==1) write(*,*) "23 Set surface dipole correction, current: X direction"
                if (idipcorr==2) write(*,*) "23 Set surface dipole correction, current: Y direction"
                if (idipcorr==3) write(*,*) "23 Set surface dipole correction, current: Z direction"
            end if
            if (istructfile==1) write(*,*) "30 Toggle giving path of geometry file instead of &COORD, current: Yes"
            if (istructfile==0) write(*,*) "30 Toggle giving path of geometry file instead of &COORD, current: No"
            read(*,*) isel2
            if (isel2==0) then
                exit
            else if (isel2==1) then
                write(*,*) "Input net charge of the system, e.g. 1"
                read(*,*) netchg
            else if (isel2==2) then
                write(*,*) "Input spin multiplicity (N_alpha - N_beta +1), e.g. 3"
                read(*,*) multispin
            else if (isel2==3) then
                write(*,*) "Input number of repetitions of the cell in three directions, e.g. 2,1,2"
                read(*,*) nrep1,nrep2,nrep3
            else if (isel2==4) then
                if (ifineXCgrid==0) then
                    ifineXCgrid=1
                else
                    ifineXCgrid=0
                end if
            else if (isel2==5) then
                write(*,*) "Input CUTOFF and REL_CUTOFF in Ry, e.g. 350,50"
                read(*,*) CUTOFF,REL_CUTOFF
            else if (isel2==6) then
                if (imoment==1) then
                    imoment=0
                else
                    imoment=1
                end if
            else if (isel2==7) then
                write(*,*) "Note: The inputted size will be applied to both sides of non-periodic direction(s)"
                if (PBCdir=="NONE") then
                    write(*,*) "Input vacuum size in Angstrom, e.g. 4.2"
                    if (iPSOLVER==4) write(*,"(a)") " Note: WAVELET Poisson solver requires cubic cell, &
                    Multiwfn will determine the longest cell size and set it to all directions"
                    read(*,*) vacsizex
                    vacsizex=vacsizex/b2a
                    vacsizey=vacsizex
                    vacsizez=vacsizex
                else if (PBCdir=="X") then
                    write(*,*) "Input vacuum size in Y and Z in Angstrom, e.g. 5.0,6.5"
                    read(*,*) vacsizey,vacsizez
                    vacsizey=vacsizey/b2a;vacsizez=vacsizez/b2a
                else if (PBCdir=="Y") then
                    write(*,*) "Input vacuum size in X and Z in Angstrom, e.g. 5.0,6.5"
                    read(*,*) vacsizex,vacsizez
                    vacsizex=vacsizex/b2a;vacsizez=vacsizez/b2a
                else if (PBCdir=="Z") then
                    write(*,*) "Input vacuum size in X and Y in Angstrom, e.g. 5.0,6.5"
                    read(*,*) vacsizex,vacsizey
                    vacsizex=vacsizex/b2a;vacsizey=vacsizey/b2a
                else if (PBCdir=="XY") then
                    write(*,*) "Input vacuum size in Z in Angstrom, e.g. 5.0"
                    read(*,*) vacsizez
                    vacsizez=vacsizez/b2a
                else if (PBCdir=="XZ") then
                    write(*,*) "Input vacuum size in Y in Angstrom, e.g. 5.0"
                    read(*,*) vacsizey
                    vacsizey=vacsizey/b2a
                else if (PBCdir=="YZ") then
                    write(*,*) "Input vacuum size in X in Angstrom, e.g. 5.0"
                    read(*,*) vacsizex
                    vacsizex=vacsizex/b2a
                end if
            else if (isel2==8) then
                if (iDFTplusU==1) then
                    iDFTplusU=0
                else
                    iDFTplusU=1
                    write(*,"(a)") " IMPORTANT NOTE: DO NOT forget to manually replace the default DFT+U parameters in the generated input file with proper value!"
                end if
            else if (isel2==9) then
                !do ikind=1,nkind
                !    write(*,"(' #',i3,':  Kind name: ',a5,' Element: ',a,'  Magnetization:',i3,'   Natoms:',i5)") &
                !    ikind,kindname(ikind),ind2name(kindeleidx(ikind)),kindmag(ikind),count(atmkind(:)==ikind)
                !end do
                do while(.true.)
                    write(*,*)
                    write(*,*) "Current magnetization status:"
                    idx=0
                    do ikind=1,nkind
                        ncount=count(atmkind(:)==ikind)
                        if (ncount>0) then
                            idx=idx+1
                            write(*,"(' #',i3,':  Kind name: ',a,' Element: ',a,'  Magnetization:',f5.2,'   Natoms:',i5)") &
                            idx,kindname(ikind),ind2name(kindeleidx(ikind)),kindmag(ikind),ncount
                        end if
                    end do
                    write(*,*)
                    write(*,*) "Input indices of the atoms to define magnetization, e.g. 1,5-10,13,19"
                    write(*,*) "To select all atoms of an element, input element name with *, e.g. *Fe"
                    write(*,*) "To change kind name, input old and new name, e.g. Fe_1 Fe_B"
                    write(*,*) "To exit, inputting ""q"""
                    read(*,"(a)") c2000tmp
                    istar=index(c2000tmp,'*')
                    if (c2000tmp=="q") then
                        multispin=nint(sum(kindmag(atmkind(:))))+1
                        exit
                    else if (istar/=0) then !Select element
                        c80tmp=(c2000tmp(istar+1:))
                        call elename2idx(c80tmp,iele)
                        ntmp=0
                        do iatm=1,ncenter
                            if (a(iatm)%index==iele) then
                                ntmp=ntmp+1
                                tmparr(ntmp)=iatm
                            end if
                        end do
                        write(*,"(i6,' atoms have been selected')") ntmp
                        if (ntmp==0) then
                            write(*,*) "Error: No atom was selected"
                            cycle
                        end if
                    else if (iachar(c2000tmp(1:1))<48.or.iachar(c2000tmp(1:1))>57) then !Change kind name
                        read(c2000tmp,*) c80tmp,c80tmp2
                        do ikind=1,nkind
                            if (trim(kindname(ikind))==trim(c80tmp)) kindname(ikind)=trim(c80tmp2)
                        end do
                        cycle
                    else !Input atom indices
                        call str2arr(c2000tmp,ntmp,tmparr)
                    end if
                    if (all(a(tmparr(1:ntmp))%index==a(tmparr(1))%index)) then
                        nkind=nkind+1
                        write(*,*) "Input magnetization (difference between alpha and beta electrons), e.g. 3.2"
                        read(*,*) kindmag(nkind)
                        atmkind(tmparr(1:ntmp))=nkind
                        iele=a(tmparr(1))%index
                        kindeleidx(nkind)=iele
                        isuffix=0
                        do while(.true.) !Test which suffix can be used
                            if (isuffix>0) then
                                write(c80tmp,"(i5)") isuffix
                                c80tmp=trim(ind2name(iele))//'_'//trim(adjustl(c80tmp))
                            else
                                c80tmp=ind2name(iele)
                            end if
                            if (all(kindname(1:nkind-1)/=trim(c80tmp))) then !New kind name
                                exit
                            else
                                ncount=count(kindname(atmkind(:))==trim(c80tmp))
                                if (ncount==0) exit
                            end if
                            isuffix=isuffix+1
                        end do
                        kindname(nkind)=trim(c80tmp)
                        write(*,*) "Done!"
                    else
                        write(*,*) "Error: Not all atoms you selected belong to the same element!"
                    end if
                end do
            else if (isel2==10) then
                write(*,*) "Choose printing level"
                write(*,*) "0 Silent"
                write(*,*) "1 Low"
                write(*,*) "2 Medium"
                write(*,*) "3 High"
                read(*,*) iprintlevel
            else if (isel2==11) then
                write(*,*) "Input external electric field vector in a.u., e.g. 0.0,0.0,0.025"
                read(*,*) efieldvec
                if (ifPBC>0.and.idiagOT==1) then
                    write(*,"(a)") " NOTE: Because this is a periodic system, &PERIODIC_EFIELD will be used in the generated input file to specify the field, &
                    in this case you must use OT rather than diagonalizaton!"
                end if
            else if (isel2==12) then
                write(*,*) "Input number of virtual orbitals to solve, e.g. 30"
                write(*,*) "If inputting -1, then all virtual orbitals will be solved" !Work since CP2K 9.1
                read(*,*) nADDED_MOS
                if (nADDED_MOS>0.and.idiagOT==2) then
                    write(*,"(a)") " NOTE: OT is changed to diagonalization since virtual orbitals cannot be generated in the case of OT"
                    idiagOT=1
                end if
            else if (isel2==13) then
                if (icentering==0) then
                    icentering=1
                else if (icentering==1) then
                    icentering=0
                end if
            else if (isel2==14) then
                if (ioutorbene==0) then
                    ioutorbene=1
                else
                    ioutorbene=0
                end if
            else if (isel2==15) then
                if (ioutSbas==0) then
                    ioutSbas=1
                else
                    ioutSbas=0
                end if
            else if (isel2==16) then
                if (ioutKSbas==0) then
                    ioutKSbas=1
                else
                    ioutKSbas=0
                end if
            else if (isel2==18) then
                if (ihyperfine==0) then
                    ihyperfine=1
                else
                    ihyperfine=0
                end if
            else if (isel2==19) then
                if (allocated(CDFTatm)) deallocate(CDFTatm,CDFTtype,CDFTnatm,CDFTtarget)
                write(*,*) "Define how many CDFT constraint groups? e.g. 3"
                write(*,*) "If input 0, CDFT will be disabled"
                read(*,*) nCDFTgroup
                if (nCDFTgroup>0) then
                    allocate(CDFTatm(ncenter,nCDFTgroup),CDFTtype(nCDFTgroup),CDFTnatm(nCDFTgroup),CDFTtarget(nCDFTgroup))
                    do igroup=1,nCDFTgroup
                        write(*,"(a,i3,a)") " Input index of the atoms involved in CDFT group",igroup,", e.g. 2,3,7-10"
                        read(*,"(a)") c2000tmp
                        call str2arr(c2000tmp,CDFTnatm(igroup),CDFTatm(:,igroup))
                        write(*,"(a,i3)") " Choose type of constraint of CDFT group",igroup
                        write(*,*) "1 Number of alpha electrons"
                        write(*,*) "2 Number of beta electrons"
                        write(*,*) "3 Number of all electrons"
                        write(*,*) "4 Spin population"
                        read(*,*) CDFTtype(igroup)
                        write(*,"(a,i3,a)") " Choose target constraint value of CDFT group",igroup,", e.g. 0.5"
                        read(*,*) CDFTtarget(igroup)
                    end do
                    write(*,"(a)") " Done! Note that atom coefficients of all CDFT groups in the generated input file will be set to 1.0, please properly modify if need"
                    if (idiagOT==1) then
                        idiagOT=2
                        write(*,"(a)") " Note: OT has been activated because it is needed by CDFT calculation"
                    end if
                end if
            else if (isel2==20) then
                !write(*,*) "Choose diagonalization library"
                !write(*,*) "1 Default"
                !write(*,*) "2 ELPA"
                !write(*,*) "3 Scalapack"
                !read(*,*) idiaglib
                if (iLRIGPW==0) then
                    iLRIGPW=1
                else
                    iLRIGPW=0
                end if
            else if (isel2==21) then
                if (iLSSCF==0) then
                    iLSSCF=1
                    if (idiagOT==2) then
                        idiagOT=1
                        write(*,"(a)") " Note: OT is disabled because it is incompatible with linear scaling self consistent field method"
                    end if
                else
                    iLSSCF=0
                end if
            else if (isel2==22) then
                write(*,*) "Choose Poisson solver:"
                write(*,*) "1 PERIODIC"
                write(*,*) "2 ANALYTIC"
                if (PBCdir=="NONE".or.PBCdir=="XY".or.PBCdir=="XZ".or.PBCdir=="YZ") write(*,*) "3 MT"
                if (PBCdir=="NONE".or.PBCdir=="XZ".or.PBCdir=="XYZ") write(*,*) "4 WAVELET"
                read(*,*) iPSOLVER
                call determine_vacuumsize(itask,iPSOLVER,vacsizex,vacsizey,vacsizez,icentering) !Automatically set proper vacuum size
            else if (isel2==23) then
                write(*,*) "0 Do not use surface dipole correction"
                write(*,*) "1 Use surface dipole correction in X direction"
                write(*,*) "2 Use surface dipole correction in Y direction"
                write(*,*) "3 Use surface dipole correction in Z direction"
                read(*,*) idipcorr
            else if (isel2==30) then
                if (istructfile==1) then
                    istructfile=0
                else
                    istructfile=1
                end if
            end if
        end do
    else if (isel==-7) then
        write(*,*) "Input one of following strings to specify periodic boundary condition (PBC)"
        write(*,*) "NONE, X, XY, XYZ, XZ, Y, YZ, Z"
        write(*,*) "If press ENTER button, NONE will be used"
        read(*,"(a)") c80tmp
        if (c80tmp==" ") then
            PBCdir="NONE"
        else
            read(c80tmp,*) PBCdir
        end if
        call strlc2uc(PBCdir) !Foolish user may input in lower case
        !Automatically set proper Poisson solver and vacuum sizes
        if (PBCdir=="NONE") then !Use WAVELET for 0D, usually best choice
            iPSOLVER=4
            write(*,*) "Note: Poisson solver has been automatically changed to WAVELET"
        else if (PBCdir=="X".or.PBCdir=="Y".or.PBCdir=="Z".or.PBCdir=="XYZ") then !Use PERIODIC for 1D and 3D
            iPSOLVER=1
            write(*,*) "Note: Poisson solver has been automatically changed to PERIODIC"
        else if (PBCdir=="XY".or.PBCdir=="XZ".or.PBCdir=="YZ") then !Use MT for 2D. Needs vaccum size in each side >= half of system
            iPSOLVER=3
            write(*,*) "Note: Poisson solver has been automatically changed to MT"
        else
            write(*,*) "Error: Your input cannot be recognized! Now set to XYZ periodicity"
            PBCdir="XYZ"
            cycle
        end if
        call determine_vacuumsize(itask,iPSOLVER,vacsizex,vacsizey,vacsizez,icentering) !Automatically set proper vacuum size
    else if (isel==-6) then
        write(*,*) "Input frequency of writing molecular dynamics trajectory, 1 means every step"
        read(*,*) nMDsavefreq
    else if (isel==-5) then
        if (itask==4) write(*,*) "Choose the format for recording trajectory of optimization"
        if (itask==6.or.itask==14) write(*,*) "Choose the format for recording trajectory of molecular dynamics"
        write(*,*) " 1 xyz (Simplest. Does not contain cell information)"
        write(*,*) " 2 dcd (Binary file, smallest size. Containing cell information)"
        write(*,"(a)") " -2 dcd_aligned_cell (same as 2, but transform so that cell vector 1 is along X-axis, cell vector 2 is in XY plane)"
        write(*,*) " 3 pdb (Containing cell information, but accuracy of coordinates is limited)"
        read(*,*) iMDformat
    else if (isel==-4) then
        write(*,*) "Printing which kind of atomic charge?"
        write(*,*) "0 None"
        write(*,*) "1 Mulliken"
        write(*,*) "2 Lowdin"
        write(*,*) "3 Hirshfeld"
        write(*,*) "4 Hirshfeld-I (its implementation in CP2K is not rigorous)"
        write(*,*) "5 Voronoi"
        write(*,*) "6 RESP"
        write(*,*) "7 REPEAT"
        read(*,*) iatomcharge
        if ((ifPBC==0.or.PBCdir=="NONE").and.iatomcharge==7) then
            write(*,*) "Error: REPEAT can only be used for periodic system"
            write(*,*) "Press ENTER button to continue"
            read(*,*)
            iatomcharge=0
        end if
        if (iatomcharge==4) then
            write(*,"(a)") " Warning: The Hirshfeld-I charge calculated by CP2K is not in line with its standard definition. Rigorous Hirshfeld-I charges can be calculated by main function 7 of Multiwfn"
            write(*,*) "Press ENTER button to continue"
            read(*,*)
        end if
        if (iatomcharge==6) then
            if (ifPBC/=0) then
                write(*,*) "Warning: REPEAT charge is a much better choice than RESP for periodic systems!"
            else
                write(*,*) "Warning: The RESP charge calculated by CP2K is not in line with its standard definition in J. Phys. Chem., 97, 10269 (1993). Rigorous RESP charges can be calculated by main function 7 of Multiwfn"
            end if
            write(*,*) "Press ENTER button to continue"
            read(*,*)
        end if
    else if (isel==-3) then
        write(*,*) "Output cube file for which real space function?"
        write(*,"(a)") " -1 Just for printing HOMO and LUMO energies as well as HOMO-LUMO gap (i.e. Outputting HOMO and LUMO cubes only)"
        write(*,*) "0 None"
        write(*,*) "1 Electron density (also with spin density for unrestricted calculation)"
        write(*,*) "2 Electron localization function (ELF)"
        write(*,*) "3 Exchange-correlation potential"
        write(*,*) "4 Hartree potential (negative of ESP)"
        write(*,*) "5 Each component of electric field"
        write(*,*) "6 Molecular orbital(s)"
        write(*,*) "7 Electron density + Hartree potential"
        read(*,*) ioutcube
        if (ioutcube==6) then
            write(*,*) "Output how many highest occupied orbitals? e.g. 5"
            write(*,*) "If inputting -1, all occupied orbitals will be outputted"
            read(*,*) NHOMO
            write(*,*) "Output how many lowest unoccupied orbitals? e.g. 5"
            write(*,*) "If inputting -1, all unoccupied orbitals will be outputted"
            read(*,*) NLUMO
        else if (ioutcube==-1) then
            ioutcube=6
            NHOMO=1
            NLUMO=1
        end if
    else if (isel==-2) then
        if (imolden==0) then
            imolden=1
        else
            imolden=0
        end if
    else if (isel==-1) then
        write(*,*) "Please select a task"
        write(*,*) "1 Energy"
        write(*,*) "2 Energy + force"
        if (PBCdir=="NONE") then
            write(*,*) "3 Optimizing structure"
        else
            write(*,*) "3 Optimizing structure (cell is fixed)"
            write(*,*) "4 Optimizing both structure and cell"
        end if
        write(*,*) "5 Vibrational analysis"
        write(*,*) "6 Molecular dynamics (MD)"
        write(*,*) "7 Searching transition state (dimer algorithm)"
        write(*,*) "8 BAND (e.g. CI-NEB)"
        write(*,*) "9 NMR"
        write(*,*) "10 Polarizability"
        write(*,*) "11 Correct for basis set superposition error (BSSE)"
        write(*,*) "13 Real-time propagation for electron dynamics"
        write(*,*) "14 Path-integral molecular dynamics (PIMD)"
        write(*,*) "15 X-ray absorption spectroscopy (XAS)"
        read(*,*) itask
        if (itask==9.and.ibas<=5) then
            ibas=10 !Use 6-31G* if current basis set is a pseudopotential basis set
            iGAPW=1
        else if (itask==4) then !Use pdb to record variable cell during cell optimization
            iMDformat=3
        else if (itask==5) then !Vibrational analysis
            iprintlevel=2 !Use medium printing level, otherwise calculation progress will not be shown
        else if (itask==11) then !BSSE
            write(*,*) "Input atoms in fragment 1, e.g. 1,3-6,10,14"
            read(*,"(a)") c2000tmp
            call str2arr(c2000tmp,nfrag1)
            if (allocated(frag1)) deallocate(frag1)
            allocate(frag1(nfrag1))
            call str2arr(c2000tmp,nfrag1,frag1)
            write(*,*) "Input charge for fragment 1, e.g. 0"
            read(*,*) frag1chg
            write(*,*) "Input spin multiplicity for fragment 1, e.g. 1"
            read(*,*) frag1multi
            write(*,*) "Input atoms in fragment 2, e.g. 2,7-9,11-13"
            read(*,"(a)") c2000tmp
            call str2arr(c2000tmp,nfrag2)
            if (allocated(frag2)) deallocate(frag2)
            allocate(frag2(nfrag2))
            call str2arr(c2000tmp,nfrag2,frag2)
            write(*,*) "Input charge for fragment 2, e.g. 0"
            read(*,*) frag2chg
            write(*,*) "Input spin multiplicity for fragment 2, e.g. 1"
            read(*,*) frag2multi
            write(*,*) "Input spin multiplicity for whole system, e.g. 1"
            read(*,*) totalmulti
        else if (itask==13) then !Real-time propagation of electron
            imoment=1
            !iatomcharge=1 !If enable this, Mulliken charge will print every step
        end if
        if (itask==1) then !Single point, medicore accuracy
            CUTOFF=350
            REL_CUTOFF=50
        else if (itask==6.or.itask==14) then !MD and PIMD, do not need high accuracy (assume without barostat)
            CUTOFF=300
            REL_CUTOFF=40
        else !If task involves energy derivative, or TDDFT, use higher cutoff
            CUTOFF=400
            REL_CUTOFF=55
        end if
        if (itask==9.or.itask==10) call determine_vacuumsize(itask,iPSOLVER,vacsizex,vacsizey,vacsizez,icentering) !Set vacuum size for NMR and polar tasks larger
        if (itask==15) then !PBEh with HFX=45% is recommended for XAS
            method="PBEh_ADMM"
            PBEh_HFX=45
            iHFX=1
            if (allocated(XASatm)) deallocate(XASatm)
            write(*,*) "Input the indices of the atoms for which XAS will be calculated, e.g. 2,4"
            write(*,*) "Note: Multiwfn assumes the atoms correspond to the same element"
            read(*,"(a)") c2000tmp
            call str2arr(c2000tmp,nXASatm)
            allocate(XASatm(nXASatm))
            call str2arr(c2000tmp,nXASatm,XASatm)
            kindnameXAS=trim(a(XASatm(1))%name)//'x'
            write(*,*) "IMPORTANT NOTE:"
            write(*,"(a)") " After choosing option 0, Multiwfn will generate XAS_TDP task. The input file corresponding to PBE45(with ADMM)/DZVP-MOLOPT-SR-GTH level with GAPW, &
            and there is a special &KIND correpsonding to pcseg-2 all-electron basis set accompanied with admm-2 auxiliary ADMM basis set. You should do following things before running the task:"
            write(*,"(a)") " 1 Copy basis set definition of pcseg-2 from BSE basis set library to a file named ""pcseg"" and put it to the ""data"" folder of CP2K"
            write(*,"(a)") " 2 Copy basis set definition of admm-2 from BSE basis set library to a file named ""pcseg-admm"" and put it to the ""data"" folder of CP2K"
            write(*,"(a)") " 3 Properly modify &XAS_TDP / STATE_TYPES according to practical requirement"
        end if
        if (itask==5) then !Vibrational analysis should use same cell as the inputted one, however for low-dimension case it will be changed, so fix using this option
            ikeepcell=1
        else
            ikeepcell=0
        end if
    else if (isel==1) then !Functionals description: https://manual.cp2k.org/trunk/CP2K_INPUT/ATOM/METHOD/XC/XC_FUNCTIONAL.html
        do while(.true.)
            !write(*,*) "-1 Molecular mechanism (MM)"
            if (iRIHFX==0) write(*,*) "0 Toggle using RI-HFX for hybrid functionals, current: No"
            if (iRIHFX==1) write(*,*) "0 Toggle using RI-HFX for hybrid functionals, current: Yes"
            write(*,*) "1 Pade (LDA)"
            write(*,*) "2 PBE               -2 revPBE            -3 PBEsol"
            write(*,*) "3 TPSS (via LibXC)   4 BP86               5 BLYP"
            write(*,*) "6 PBE0        -6 PBE0 with ADMM"
            write(*,*) "7 B3LYP       -7 B3LYP with ADMM"
            write(*,*) "8 HSE06       -8 HSE06 with ADMM"
            write(*,*) "9 BHandHLYP   -9 BHandHLYP with ADMM"
            write(*,*) "10 M06-2X    -10 M06-2X with ADMM"
            write(*,*) "11 B97M-rV (via LibXC)       12 MN15L (via LibXC)"
            write(*,*) "13 SCAN (via LibXC)          14 r2SCAN (via LibXC)"
            write(*,*) "15 RPBE (via LibXC)          16 revTPSS (via LibXC)"
            write(*,*) "17 BEEF-vdW                  18 HLE17 (via LibXC)"
            write(*,*) "20 RI-MP2     21 RI-SCS-MP2    22 RI-(EXX+RPA)@PBE"
            write(*,*) "25 RI-B2PLYP  26 RI-B2GP-PLYP  27 RI-DSD-BLYP  28 RI-revDSD-PBEP86 with ADMM"
            write(*,*) "30 GFN1-xTB   31 GFN0-xTB      32 GFN2-xTB (via tblite)"
            write(*,*) "40 PM6        50 SCC-DFTB + disp. corr."
            write(*,*) "60 GW@BHandHLYP with ADMM      61 GW@MN15L"
            write(*,*) "80 PBEh      -80 PBEh with ADMM  (customize HFX composition)"
            write(*,*) "100 FIST module (molecular mechanics)"
            read(*,*) isel2
            if (isel2==0) then
                if (iRIHFX==0) then
                    iRIHFX=1
                else
                    iRIHFX=0
                end if
                cycle
            end if
            iwfc=0 !Do not involve wavefunction-based correlation
            iHFX=0 !Do not involve HF exchange
            if (isel2==1) method="Pade"
            if (isel2==2) method="PBE"
            if (isel2==-2) method="revPBE"
            if (isel2==-3) method="PBEsol"
            if (isel2==3) method="TPSS_LIBXC"
            if (isel2==4) method="BP"
            if (isel2==5) method="BLYP"
            if (isel2==6) method="PBE0"
            if (isel2==-6) method="PBE0_ADMM"
            if (isel2==7) method="B3LYP"
            if (isel2==-7) method="B3LYP_ADMM"
            if (isel2==8) method="HSE06"
            if (isel2==-8) method="HSE06_ADMM"
            if (isel2==9) method="BHandHLYP"
            if (isel2==-9) method="BHandHLYP_ADMM"
            if (isel2==10) method="M06-2X"
            if (isel2==-10) method="M06-2X_ADMM"
            if (abs(isel2)==80) then
                if (isel2==80) method="PBEh"
                if (isel2==-80) method="PBEh_ADMM"
                write(*,*) "Input HF exchange hybrid composition (%), e.g. 45"
                read(*,*) PBEh_HFX
            end if
            if (isel2==100) method="FIST"
            if (isel2==11) then
                method="B97M-rV_LIBXC"
                idispcorr=5
            end if
            if (isel2==12) method="MN15L_LIBXC"
            if (isel2==13) method="SCAN_LIBXC"
            if (isel2==14) method="r2SCAN_LIBXC"
            if (isel2==15) method="RPBE_LIBXC"
            if (isel2==16) method="revTPSS_LIBXC"
            if (isel2==17) then
                method="BEEFVDW"
                idispcorr=0
            end if
            if (isel2==18) method="HLE17_LIBXC"
            if ((isel2>=20.and.isel2<=29).or.(isel2>=60.and.isel2<=61)) then !Involve wavefunction-based correlation
                if (isel2==20) method="RI-MP2"
                if (isel2==21) method="RI-SCS-MP2"
                if (isel2==22) method="RI-(EXX+RPA)@PBE"
                if (isel2==25) method="RI-B2PLYP"
                if (isel2==26) method="RI-B2GP-PLYP"
                if (isel2==27) method="RI-DSD-BLYP"
                if (isel2==28) method="RI-revDSD-PBEP86_ADMM"
                if (isel2>=60.and.isel2<=61) then
                    if (isel2==60) method="GW@BHandHLYP_ADMM"
                    if (isel2==61) method="GW@MN15L"
                    write(*,*) "Use which type of GW?"
                    write(*,*) "1 G0W0"
                    write(*,*) "2 evGW"
                    write(*,*) "3 scGW0"
                    read(*,*) itmp
                    niter_evGW=1
                    niter_scGW0=1
                    if (itmp==2) then
                        niter_evGW=10
                    else if (itmp==3) then
                        niter_scGW0=10
                    end if
                end if
                iwfc=1
                if (ibas/=20) ibas=21 !Default to cc-TZ with RI_TZ
                if (index(method,"GW")/=0) ibas=22 !Change to QZ level because slow basis set convergence feature of GW
            end if
            if (isel2==30) method="GFN1-xTB"
            if (isel2==31) method="GFN0-xTB"
            if (isel2==32) method="GFN2-xTB"
            if (isel2==40) method="PM6"
            if (isel2==50) method="SCC-DFTB"
            if (index(method,"ADMM")/=0.or.isel2==30.or.isel2==40.or.isel2==50) then !When ADMM is used, OT is suggested to be used. OT is suggested for GFN1-xTB, PM6, SCC-DFTB dealing with large systems
                if (ikpoint1==1.and.ikpoint2==1.and.ikpoint3==1) idiagOT=2
            end if
            if (isel2==32.or.isel==40) imixing=1
            if (index(method,"SCAN")/=0) then
                write(*,"(a)") " NOTE: If you are using CP2K >=9.1, in the generated CP2K input file, it is suggested to replace &
                ""POTENTIAL_FILE_NAME  POTENTIAL"" with ""POTENTIAL_FILE_NAME  POTENTIAL_UZH"", and replace ""BASIS_SET_FILE_NAME  BASIS_MOLOPT"" with ""BASIS_SET_FILE_NAME  BASIS_MOLOPT_UZH"", &
                and manually specify proper GTH potential and corresponding valence basis set optimized for SCAN calculation"
            else if (index(method,"PBE0")/=0) then
                write(*,"(a)") " NOTE: If you are using CP2K >=9.1, in the generated CP2K input file, it is suggested to replace &
                ""POTENTIAL_FILE_NAME  POTENTIAL"" with ""POTENTIAL_FILE_NAME  POTENTIAL_UZH"", and replace ""BASIS_SET_FILE_NAME  BASIS_MOLOPT"" with ""BASIS_SET_FILE_NAME  BASIS_MOLOPT_UZH"", &
                and manually specify proper GTH potential and corresponding valence basis set optimized for PBE0 calculation"
            end if
            if (abs(isel2)==6.or.abs(isel2)==7.or.abs(isel2)==8.or.abs(isel2)==9.or.abs(isel2)==10.or.abs(isel2)==80.or.iwfc==1) then !Hybrid functionals and wavefunction-based correlation methods need HF exchange
                iHFX=1
                if (method=="RI-(EXX+RPA)@PBE".or.method=="GW@MN15L") iHFX=0 !Based on PBE orbitals, HFX is not involved in SCF process
            end if
            exit
        end do
        
    else if (isel==2) then !Select basis set
       write(*,*) "  GTH pseudopotential basis sets:"
        do ibas=-10,9
            if (basname(ibas)/=" ") write(*,"(1x,i2,1x,a)") ibas,trim(basname(ibas))
        end do
        do ibas=20,30
            if (basname(ibas)/=" ") write(*,"(1x,i2,1x,a)") ibas,trim(basname(ibas))
        end do
        write(*,*) "  All-electron basis sets:"
        do ibas=10,19
            if (basname(ibas)/=" ") write(*,"(1x,i2,1x,a)") ibas,trim(basname(ibas))
        end do
        read(*,*) ibassel
        if (iwfc==1.and.ibassel/=20.and.ibassel/=21.and.ibassel/=22) then
            write(*,"(a)") " Error: To perform RI calculation, you must choose 20 or 21, because they are accompanied by auxiliary basis set" 
            write(*,*) "Press ENTER button to continue"
            read(*,*)
        else
            ibas=ibassel
            if (ibas>=10.and.ibas<=16) then
                iGAPW=1
            else
                iGAPW=0
            end if
        end if
    else if (isel==3) then
        write(*,*) "Choose dispersion correction method"
        write(*,*) "0 None"
        write(*,*) "1 DFT-D3"
        write(*,*) "2 DFT-D3(BJ)"
        write(*,*) "3 DFT-D4"
        write(*,*) "5 rVV10"
        read(*,*) idispcorr
    else if (isel==4) then
        if (idiagOT==1) then
            if (ikpoint1/=1.or.ikpoint2/=1.or.ikpoint3/=1) then
                write(*,*) "Error: OT can only be used for Gamma point!"
                write(*,*) "Press ENTER button to continue"
                read(*,*)
                cycle
            end if
            if (ismear==1) then
                write(*,*) "Error: OT cannot be used in combination with smearing!"
                write(*,*) "Press ENTER button to continue"
                read(*,*)
                cycle
            end if
            idiagOT=2
        else if (idiagOT==2) then
            idiagOT=1
        end if
    else if (isel==5) then
        if (idiagOT==1) then
            write(*,*) "Choose how to mixing old and new density matrices"
            write(*,*) "1 Direct mixing with DIIS (default, usually poor)"
            write(*,*) "2 Broyden mixing"
            write(*,*) "3 Pulay mixing"
            read(*,*) imixing
        else if (idiagOT==2) then
            if (iouterSCF==0) then
                iouterSCF=1
            else
                iouterSCF=0
            end if
        end if
    else if (isel==6) then
        if (ismear==0) then
            ismear=1
            nADDED_MOS=nint(ncenter/2D0) !This is usually sufficient for 300 K smearing
            if (nADDED_MOS<30) nADDED_MOS=30
            write(*,"(a,i5,a)") " Note: The number of virtual orbitals to solve has been changed to",nADDED_MOS,", please properly adjust if need"
        else
            ismear=0
            nADDED_MOS=0
            write(*,"(a)") " Note: The number of virtual orbitals to solve has been changed to 0"
        end if
    else if (isel==7) then
        if (iSCCS==0) then
            iSCCS=1
            write(*,"(a)") " Note: In the generated input file, the SCCS parameters correspond to water case given in &
            J. Chem. Phys., 150, 041710 (2019), please manually modify according to actual situation"
        else
            iSCCS=0
        end if
    else if (isel==8) then
        write(*,*) "Input number of k-points of MONKHORST-PACK in three directions, e.g. 8,6,2"
        read(*,*) ikpoint1,ikpoint2,ikpoint3
        if (idiagOT==2.and.(ikpoint1/=1.or.ikpoint2/=1.or.ikpoint3/=1)) then
            write(*,"(a)") " Warning: OT can be used for Gamma point only! Now diagonalization is used instead"
            idiagOT=1
        end if
    else if (isel==9) then
        do while(.true.)
            write(*,*) "Input indices of the atoms to be freezed, e.g. 1,5,9-12,14-18"
            write(*,"(a)") " If inputting ""optH"", then only hydrogens will be optimized while others will be fixed"
            read(*,"(a)") c2000tmp
            if (.not.allocated(atmcons)) allocate(atmcons(ncenter))
            if (index(c2000tmp,"optH")/=0) then
                natmcons=0
                do iatm=1,ncenter
                    if (a(iatm)%index==1) cycle
                    natmcons=natmcons+1
                    atmcons(natmcons)=iatm
                end do
            else
                call str2arr(c2000tmp,natmcons,atmcons)
                if (natmcons>ncenter) then
                    write(*,*) "Error: The indices you inputted is invalid!"
                    cycle
                end if
            end if
            write(*,"(i8,' atoms will be fixed')") natmcons
            write(*,"(a)") " Note that the direction(s) of fixing can be manually set by changing ""COMPONENTS_TO_FIX"" &
            in the generated input file, by default they are fixed in all directions"
            exit
        end do
    else if (isel==10) then
        if (itask==6) then
            write(*,*) "0 Do not use thermostat"
            write(*,*) "1 Adaptive-Langevin thermostat"
            write(*,"(a)") " 2 Canonical sampling through velocity rescaling (CSVR, also known as V-rescale, recommended!)"
            write(*,*) "3 Generalized Langevin Equation (GLE) thermostat"
            write(*,*) "4 Nose-Hoover thermostat"
            read(*,*) ithermostat
            !if (ithermostat>0) then
            !    if (.not.allocated(thermoatm)) allocate(thermoatm(ncenter))
            !    nthermoatm=ncenter
            !    forall(iatm=1:ncenter) thermoatm(iatm)=iatm
            !end if
        else if (itask==3.or.itask==4) then
            write(*,*) "Choose optimization method"
            write(*,*) "1 BFGS (Best choice for most situations)"
            write(*,*) "2 LBFGS (Suitable for very large systems)"
            write(*,"(a)") " 3 Conjugate gradient (More robust than BFGS and LBFGS especially when initial geometry &
            is far from minimum, unfortunately more expensive. Try it for difficult cases)"
            read(*,*) ioptmethod
        else if (itask==5) then
            if (ioutvibmol==0) then
                ioutvibmol=1
            else
                ioutvibmol=0
            end if
        else if (itask==15) then
            if (iXAS_SOC==0) then
                iXAS_SOC=1
            else
                iXAS_SOC=0
            end if
        end if
    else if (isel==11) then
        if (itask==6) then
            !write(*,*) "Input indices of the atoms to whom the thermostat will be applied"
            !write(*,*) "For example: 1,5,9-12,14-18"
            !read(*,"(a)") c2000tmp
            !if (.not.allocated(thermoatm)) allocate(thermoatm(ncenter))
            !call str2arr(c2000tmp,nthermoatm,thermoatm)
        else if (itask==4) then
            write(*,*) "Input one of following constraints on cell optimization"
            write(*,*) "NONE, X, Y, Z, XY, XZ, YZ"
            read(*,"(a)") cellfix
        else if (itask==5) then
            if (ikeepcell==0) then
                ikeepcell=1
            else
                ikeepcell=0
            end if
        else if (itask==15) then
            if (iGW2X==0) then
                iGW2X=1
            else
                iGW2X=0
            end if
        end if
    else if (isel==12) then
        if (itask==6) then
            write(*,*) "0 Do not use barostat"
            write(*,*) "1 Use barostat, flexible cell"
            write(*,*) "2 Use barostat, isotropic cell"
            read(*,*) ibarostat
            if (ibarostat>0) iMDformat=2 !Use dcd format to record variable cell size during MD
        else if (itask==4) then
            write(*,*) "0 Return"
            write(*,*) "1 Set isotropic external pressure"
            write(*,*) "2 Set anisotropic external pressure by specifying 9 pressure tensor components"
            read(*,*) isel2
            if (isel2==0) then
                cycle
            else if (isel2==1) then
                write(*,*) "Input external pressure in bar, e.g. 150"
                read(*,*) Piso
                iprestype=1
            else if (isel2==2) then
                do while(.true.)
                    write(*,*)
                    write(*,*) "Current pressure tensor:"
                    call showmatgau(Ptens,fileid=6,form="1PE14.5")
                    write(*,*) "Input for example 1,3,250.4 to define component XZ as 250.4"
                    write(*,*) "Input ""q"" can return"
                    read(*,"(a)") c80tmp
                    if (index(c80tmp,'q')/=0) then
                        iprestype=2
                        exit
                    else
                        read(c80tmp,*,iostat=ierror) idx,jdx,presval
                        if (ierror/=0) then
                            write(*,*) "Error: Unable to recognize your input!"
                            cycle
                        end if
                        Ptens(idx,jdx)=presval
                    end if
                end do
            end if
        end if
    else if (isel==13) then
        if (itask==2.or.itask==3.or.itask==4.or.itask==7) then
            if (itightopt==0) then
                itightopt=1
            else if (itightopt==1) then
                itightopt=0
            end if
        else if (itask==5) then
            if (iraman==0) then
                iraman=1
            else if (iraman==1) then
                iraman=0
            end if
        else if (itask==6) then
            if (inoSCFinfo==0) then
                inoSCFinfo=1
            else
                inoSCFinfo=0
            end if
        end if
    else if (isel==15) then
        if (iTDDFT==0) then
            iTDDFT=1
            write(*,*) "Note: The TDDFT realized by CP2K employs Tamm-Dancoff approximation"
            write(*,"(/,a)") " If outputting .molden file containing all occupied and a batch of virtual orbitals for post-processing analysis? (y/n)"
            read(*,*) selectyn
            if (selectyn=='y') then
                write(*,"(a)") " How many virtual orbitals to solve and record in the .molden file? e.g. 40"
                write(*,*) "You can input -1 or a very large number to solve all virtual orbitals"
                read(*,*) nADDED_MOS
                imolden=1
                idiagOT=1
            end if
        else
            iTDDFT=0
        end if
    !Below are specific for TDDFT
    else if (isel==16) then
        write(*,*) "Input number of excited states to solve, e.g. 5"
        read(*,*) nstates_TD
    else if (isel==17) then
        if (iTDtriplet==0) then
            iTDtriplet=1
        else
            iTDtriplet=0
        end if
    else if (isel==18) then
        if (isTDA==0) then
            isTDA=1
        else
            isTDA=0
        end if
    else if (isel==19) then
        if (iNTO==0) then
            iNTO=1
        else
            iNTO=0
        end if
    else if (isel==20) then
        if (itask==2.or.itask==3.or.itask==4.or.itask==5.or.itask==6.or.itask==7.or.itask==8) then !Need force
            write(*,*) "Input the index of the excited state for which force will be evaluated, e.g. 2"
            read(*,*) istate_force
        end if
    else if (isel==21) then
        if (iSOCTDDFT==0) then
            iSOCTDDFT=1
        else
            iSOCTDDFT=0
        end if
    
    else if (isel==0) then
        exit
    end if
end do

open(ifileid,file=outname,status="replace")
write(ifileid,"(a)") "#Generated by Multiwfn (http://sobereva.com/multiwfn)"
call path2filename(filename,c200tmp)
write(ifileid,"(a)") "&GLOBAL"
call path2filename(outname,c200tmp)
write(ifileid,"(a)") "  PROJECT "//trim(c200tmp)
if (idiaglib==2) write(ifileid,"(a)") "  PREFERRED_DIAG_LIBRARY ELPA #Library for diagonalization"
if (idiaglib==3) write(ifileid,"(a)") "  PREFERRED_DIAG_LIBRARY SL #Library for diagonalization"
if (iprintlevel==0) write(ifileid,"(a)") "  PRINT_LEVEL SILENT"
if (iprintlevel==1) write(ifileid,"(a)") "  PRINT_LEVEL LOW"
if (iprintlevel==2) write(ifileid,"(a)") "  PRINT_LEVEL MEDIUM"
if (iprintlevel==3) write(ifileid,"(a)") "  PRINT_LEVEL HIGH"
if (itask==1.or.itask==15) write(ifileid,"(a)") "  RUN_TYPE ENERGY"
if (itask==2) write(ifileid,"(a)") "  RUN_TYPE ENERGY_FORCE"
if (itask==3) write(ifileid,"(a)") "  RUN_TYPE GEO_OPT"
if (itask==4) write(ifileid,"(a)") "  RUN_TYPE CELL_OPT"
if (itask==5) write(ifileid,"(a)") "  RUN_TYPE VIBRATIONAL_ANALYSIS"
if (itask==6) write(ifileid,"(a)") "  RUN_TYPE MD"
if (itask==7) write(ifileid,"(a)") "  RUN_TYPE GEO_OPT"
if (itask==8) write(ifileid,"(a)") "  RUN_TYPE BAND"
if (itask==9.or.itask==10) write(ifileid,"(a)") "  RUN_TYPE LR"
if (itask==11) write(ifileid,"(a)") "  RUN_TYPE BSSE"
if (itask==13) write(ifileid,"(a)") "  RUN_TYPE RT_PROPAGATION"
if (itask==14) write(ifileid,"(a)") "  RUN_TYPE PINT"
write(ifileid,"(a)") "&END GLOBAL"
write(ifileid,"(/,a)") "&FORCE_EVAL"
if (method=="FIST") then
    write(ifileid,"(a)") "  METHOD FIST"
else
    write(ifileid,"(a)") "  METHOD Quickstep"
end if

write(ifileid,"(a)") "  &SUBSYS"
if (nrep1/=1.or.nrep2/=1.or.nrep3/=1.or.icentering==1.or.method=="FIST".or.istructfile==1) then
    write(ifileid,"(a)") "    &TOPOLOGY"
    if (nrep1/=1.or.nrep2/=1.or.nrep3/=1) write(ifileid,"(a,3i3)") "      MULTIPLE_UNIT_CELL",nrep1,nrep2,nrep3
    if (icentering==1.and.ikeepcell==0) then
        write(ifileid,"(a)") "      &CENTER_COORDINATES #Centering atoms in the box"
        write(ifileid,"(a)") "      &END CENTER_COORDINATES"
    end if
    if (istructfile==1) then
        inamelen=len_trim(filename)
        if (filename(inamelen-2:inamelen)=="xyz".or.filename(inamelen-2:inamelen)=="XYZ") then
            write(ifileid,"(a)") "      COORD_FILE_FORMAT XYZ #The format of the file providing coordinates"
        else if (filename(inamelen-2:inamelen)=="pdb".or.filename(inamelen-2:inamelen)=="PDB") then
            write(ifileid,"(a)") "      COORD_FILE_FORMAT PDB #The format of the file providing coordinates"
        else if (filename(inamelen-2:inamelen)=="cif".or.filename(inamelen-2:inamelen)=="CIF") then
            write(ifileid,"(a)") "      COORD_FILE_FORMAT CIF #The format of the file providing coordinates"
        else if (filename(inamelen-2:inamelen)=="gro".or.filename(inamelen-2:inamelen)=="GRO") then
            write(ifileid,"(a)") "      COORD_FILE_FORMAT G96 #The format of the file providing coordinates"
        else if (filename(inamelen-6:inamelen)=="restart") then
            write(ifileid,"(a)") "      COORD_FILE_FORMAT CP2K #The format of the file providing coordinates"
        else
            write(ifileid,"(a)") "      COORD_FILE_FORMAT ? #The format of the file providing coordinates"
        end if
        call path2filename(filename,c200tmp) !Remove folder part
        ipos=index(filename,'.',back=.true.)
        write(ifileid,"(a)") "      COORD_FILE_NAME "//trim(c200tmp)//trim(filename(ipos:))
    end if
    if (method=="FIST") then
        write(ifileid,"(a)") "      #CONN_FILE_FORMAT PSF #Determining connectivity from topology file rather than automatic generating"
        write(ifileid,"(a)") "      #CONN_FILE_NAME foo.psf"
        if (istructfile==0) then
            write(ifileid,"(a)") "      #COORD_FILE_FORMAT PDB #The format of the file providing coordinates"
            write(ifileid,"(a)") "      #COORD_FILE_NAME foo.pdb"
        end if
        write(ifileid,"(a)") "      #&DUMP_PSF #Dumping the .psf file containing connectivity"
        write(ifileid,"(a)") "      #&END DUMP_PSF"
        write(ifileid,"(a)") "      #&DUMP_PDB #Dumping .pdb file at starting geometry"
        write(ifileid,"(a)") "      #  CHARGE_BETA #Write atomic charges to beta field"
        write(ifileid,"(a)") "      #&END DUMP_PDB"
    end if
    write(ifileid,"(a)") "    &END TOPOLOGY"
end if

!---- &CELL
!A pseudo-cell cellv1/2/3_pseudo is the cell actually presented in the input file
!According to PBCdir, the periodic direction will employ the original cell information loaded from input file,&
!while other direction(s) employ system size plus vacuum size, and the vectors are parallel to the corresponding Cartesian axes
cellv1_pseudo=cellv1
cellv2_pseudo=cellv2
cellv3_pseudo=cellv3
xdist=(maxval(a%x)-minval(a%x)+2*vacsizex) !Extended size in X
ydist=(maxval(a%y)-minval(a%y)+2*vacsizey)
zdist=(maxval(a%z)-minval(a%z)+2*vacsizez)
write(ifileid,"(a)") "    &CELL"
if (ikeepcell==0) then
    if (PBCdir=="NONE") then
        if (iPSOLVER==4) then !WAVELET, needs cubic box
            tmp=max(max(xdist,ydist),zdist)
            cellv1_pseudo(:)=(/ tmp,0D0,0D0 /)
            cellv2_pseudo(:)=(/ 0D0,tmp,0D0 /)
            cellv3_pseudo(:)=(/ 0D0,0D0,tmp /)
        else
            cellv1_pseudo(:)=(/ xdist,0D0,0D0 /)
            cellv2_pseudo(:)=(/ 0D0,ydist,0D0 /)
            cellv3_pseudo(:)=(/ 0D0,0D0,zdist /)
        end if
    else if (PBCdir=="X") then
        cellv2_pseudo(:)=(/ 0D0,ydist,0D0 /)
        cellv3_pseudo(:)=(/ 0D0,0D0,zdist /)
    else if (PBCdir=="Y") then
        cellv1_pseudo(:)=(/ xdist,0D0,0D0 /)
        cellv3_pseudo(:)=(/ 0D0,0D0,zdist /)
    else if (PBCdir=="Z") then
        cellv1_pseudo(:)=(/ xdist,0D0,0D0 /)
        cellv2_pseudo(:)=(/ 0D0,ydist,0D0 /)
    else if (PBCdir=="XY") then
        cellv3_pseudo(:)=(/ 0D0,0D0,zdist /)
    else if (PBCdir=="XZ") then
        cellv2_pseudo(:)=(/ 0D0,ydist,0D0 /)
    else if (PBCdir=="YZ") then
        cellv1_pseudo(:)=(/ xdist,0D0,0D0 /)
    else if (PBCdir=="XYZ") then
        !When ifPBC==3, namely the inputted file provides cellv1/2/3, it will be directly used
        !while if the loaded system is nonperiodic while we request to use XYZ periodicity to calculate, then employ extended sizes
        if (ifPBC==0) then
            cellv1_pseudo(:)=(/ xdist,0D0,0D0 /)
            cellv2_pseudo(:)=(/ 0D0,ydist,0D0 /)
            cellv3_pseudo(:)=(/ 0D0,0D0,zdist /)
        else if (ifPBC==1) then
            cellv2_pseudo(:)=(/ 0D0,ydist,0D0 /)
            cellv3_pseudo(:)=(/ 0D0,0D0,zdist /)
            write(*,"(a)") " Note: The translation vector in your input file is assumed to correspond to X direction, the cell sizes of Y and Z are automatically determined"
        else if (ifPBC==2) then
            cellv3_pseudo(:)=(/ 0D0,0D0,zdist /)
            write(*,"(a)") " Note: The translation vectors in your input file are assumed to correspond to the directions in XY plane, the cell size of Z is automatically determined"
        end if
    end if
end if
write(ifileid,"(a,3f15.8)") "      A",cellv1_pseudo(:)*b2a
write(ifileid,"(a,3f15.8)") "      B",cellv2_pseudo(:)*b2a
write(ifileid,"(a,3f15.8)") "      C",cellv3_pseudo(:)*b2a
if (iPSOLVER==1) then !PSOLVER PERIODIC forces to use XYZ periodicity
    write(ifileid,"(a)") "      PERIODIC XYZ #Direction(s) of applied PBC (geometry aspect)"
    if (PBCdir/="XYZ") write(*,"(a)") " Note: PERIODIC in the generated input file is changed to XYZ since current PSOLVER is PERIODIC"
else
    write(ifileid,"(a)") "      PERIODIC "//trim(PBCdir)//" #Direction(s) of applied PBC (geometry aspect)"
end if
if (nrep1/=1.or.nrep2/=1.or.nrep3/=1) write(ifileid,"(a,3i3)") "      MULTIPLE_UNIT_CELL",nrep1,nrep2,nrep3
write(ifileid,"(a)") "    &END CELL"

!---- &COORD
if (istructfile==0) then
    write(ifileid,"(a)") "    &COORD"
    do iatm=1,ncenter
        if (itask==15.or.allocated(XASatm)) then !XAS, use special KIND for excited atom
            if (any(XASatm(:)==iatm)) then
                write(ifileid,"(6x,a,3f14.8)") kindnameXAS,a(iatm)%x*b2a,a(iatm)%y*b2a,a(iatm)%z*b2a
            else
                write(ifileid,"(6x,a,3f14.8)") kindname(atmkind(iatm)),a(iatm)%x*b2a,a(iatm)%y*b2a,a(iatm)%z*b2a
            end if
        else
            write(ifileid,"(6x,a,3f14.8)") kindname(atmkind(iatm)),a(iatm)%x*b2a,a(iatm)%y*b2a,a(iatm)%z*b2a
        end if
    end do
    write(ifileid,"(a)") "    &END COORD"
    if (itask==6) then
        write(ifileid,"(a)") "#   &VELOCITY #You can set initial atomic velocities in this section"
        write(ifileid,"(a)") "#   &END VELOCITY"
    end if
end if

!---- &KIND
if (method(1:3)/="GFN".and.method/="PM6".and.method/="SCC-DFTB".and.method/="FIST") then !Semi-empirical methods and forcefield do not need to define these
    ntime=1
    if (itask==11) ntime=2 !For BSSE task, the second time write same kinds but with _ghost suffix
    do itime=1,ntime
        do ikind=1,nkind
            if (count(atmkind(:)==ikind)==0) cycle
            if (itime==1) then
                write(ifileid,"(a)") "    &KIND "//trim(kindname(ikind))
            else if (itime==2) then
                write(ifileid,"(a)") "    &KIND "//trim(kindname(ikind))//"_ghost"
            end if
            write(ifileid,"(a)") "      ELEMENT "//trim(ind2name(kindeleidx(ikind)))
            if (ibas==19) then
                write(ifileid,"(a)") "      BASIS_SET def2-QZVP"
                write(ifileid,"(a)") "      BASIS_SET RI_AUX RI-5Z"
            else if (ibas==20) then
                write(ifileid,"(a)") "      BASIS_SET cc-DZ"
                write(ifileid,"(a)") "      BASIS_SET RI_AUX RI_DZ"
            else if (ibas==21) then
                write(ifileid,"(a)") "      BASIS_SET cc-TZ"
                write(ifileid,"(a)") "      BASIS_SET RI_AUX RI_TZ"
            else if (ibas==22) then
                write(ifileid,"(a)") "      BASIS_SET cc-QZ"
                write(ifileid,"(a)") "      BASIS_SET RI_AUX RI_QZ"
            else
                if ((ibas>=-6.and.ibas<=-1).or.(ibas>=1.and.ibas<=5).or.ibas==7.or.ibas==8.or.ibas==9) then
                    write(c80tmp,"(i3)") Nval(kindeleidx(ikind))
                    write(ifileid,"(a)") "      BASIS_SET "//trim(basname(ibas))//'-q'//trim(adjustl(c80tmp))
                else
                    write(ifileid,"(a)") "      BASIS_SET "//trim(basname(ibas))
                end if
            end if
            if (index(method,"ADMM")/=0) then !Set auxiliary basis set file for ADMM
                if (iGAPW==1) then !All-electron basis set
                    write(ifileid,"(a)") "      BASIS_SET AUX_FIT pob-DZVP-rev2" !Small but good enough for ADMM
                else !Pseudopotential basis set
                    write(c80tmp,"(i3)") Nval(kindeleidx(ikind))
                    write(ifileid,"(a)") "      BASIS_SET AUX_FIT admm-dzp"//'-q'//trim(adjustl(c80tmp)) !Small but good enough for ADMM
                end if
            end if
            if (iLRIGPW==1) then !Set auxiliary basis set file for LRIGPW
                write(ifileid,"(a)") "      BASIS_SET LRI_AUX LRI-DZVP-MOLOPT-GTH-MEDIUM"
            end if
            if (iGAPW==0) then !GPW
                if (method=="Pade") then
                    write(ifileid,"(a)") "      POTENTIAL GTH-PADE"
                else if (method=="BP") then          
                    write(ifileid,"(a)") "      POTENTIAL GTH-BP"
                else if (method=="BLYP".or.method=="B3LYP") then
                    write(ifileid,"(a)") "      POTENTIAL GTH-BLYP"
                !else if (index(method,"SCAN")/=0) then
                !    write(ifileid,"(a)") "      POTENTIAL GTH-SCAN"
                !else if (index(method,"PBE0")/=0) then
                !    write(ifileid,"(a)") "      POTENTIAL GTH-PBE0"
                else if (index(method,"MP2")/=0) then
                    write(ifileid,"(a)") "      POTENTIAL GTH-HF"
                else if (index(method,"revDSD-PBEP86")/=0) then
                    write(c80tmp,"(i3)") Nval(kindeleidx(ikind))
                    write(ifileid,"(a)") "      POTENTIAL GTH-PBE0"//'-q'//trim(adjustl(c80tmp))
                else                           
                    write(ifileid,"(a)") "      POTENTIAL GTH-PBE"
                end if
            else !GAPW
                write(ifileid,"(a)") "      POTENTIAL ALL"
            end if
            if (itime==2.and.itask==11) write(ifileid,"(a)") "      GHOST"
            if (iDFTplusU==1) then !DFT+U
                ie=kindeleidx(ikind)
                if ((ie>=57.and.ie<=70).or.(ie>=89.and.ie<=102)) then !f elements
                    write(ifileid,"(a)") "      &DFT_PLUS_U"
                    write(ifileid,"(a)") "        L 3 #Quantum number of angular momentum the atomic orbitals to +U. 0=s, 1=p, 2=d, 3=f"
                    write(ifileid,"(a)") "        U_MINUS_J [eV] 2.0 #Effective on-site Coulomb interaction parameter U(eff) = U - J"
                    write(ifileid,"(a)") "      &END DFT_PLUS_U"
                else if ((ie>=21.and.ie<=28).or.(ie>=39.and.ie<=46).or.(ie>=71.and.ie<=78)) then !d elements
                    write(ifileid,"(a)") "      &DFT_PLUS_U"
                    write(ifileid,"(a)") "        L 2 #Quantum number of angular momentum the atomic orbitals to +U. 0=s, 1=p, 2=d, 3=f"
                    write(ifileid,"(a)") "        U_MINUS_J [eV] 2.0 #Effective on-site Coulomb interaction parameter U(eff) = U - J"
                    write(ifileid,"(a)") "      &END DFT_PLUS_U"
                else if ((ie>=5.and.ie<=9).or.(ie>=13.and.ie<=17).or.(ie>=31.and.ie<=35).or.(ie>=49.and.ie<=53).or.(ie>=81.and.ie<=85)) then !p elements
                    write(ifileid,"(a)") "      #&DFT_PLUS_U"
                    write(ifileid,"(a)") "      #  L 1 #Quantum number of angular momentum the atomic orbitals to +U. 0=s, 1=p, 2=d, 3=f"
                    write(ifileid,"(a)") "      #  U_MINUS_J [eV] 2.0 #Effective on-site Coulomb interaction parameter U(eff) = U - J"
                    write(ifileid,"(a)") "      #&END DFT_PLUS_U"
                else if ((ie>=3.and.ie<=4).or.(ie>=11.and.ie<=12).or.(ie>=19.and.ie<=20).or.(ie>=37.and.ie<=38).or.(ie>=55.and.ie<=56).or.(ie>=87.and.ie<=88)) then !s elements
                    write(ifileid,"(a)") "      #&DFT_PLUS_U"
                    write(ifileid,"(a)") "      #  L 0 #Quantum number of angular momentum the atomic orbitals to +U. 0=s, 1=p, 2=d, 3=f"
                    write(ifileid,"(a)") "      #  U_MINUS_J [eV] 2.0 #Effective on-site Coulomb interaction parameter U(eff) = U - J"
                    write(ifileid,"(a)") "      #&END DFT_PLUS_U"
                end if
            end if
            if (kindmag(ikind)/=0) write(ifileid,"(a,f6.2)") "      MAGNETIZATION",kindmag(ikind)
            write(ifileid,"(a)") "    &END KIND"
        end do
        if (itask==15.or.allocated(XASatm)) then !For XAS_TDP, provide a dummy KIND so that user can easily specify the atom to be excited
            write(ifileid,"(a)") "    &KIND "//kindnameXAS//" #For atom to be excited in XAS task"
            write(ifileid,"(a)") "      ELEMENT "//trim(a(XASatm(1))%name)
            if (iGW2X==1) then
                write(ifileid,"(a)") "      BASIS_SET pcseg-3"
                if (itask==15) write(ifileid,"(a)") "      BASIS_SET AUX_FIT admm-3"
            else
                write(ifileid,"(a)") "      BASIS_SET pcseg-2"
                if (itask==15) write(ifileid,"(a)") "      BASIS_SET AUX_FIT admm-2"
            end if
            write(ifileid,"(a)") "      POTENTIAL ALL"
            write(ifileid,"(a)") "    &END KIND"
        end if
    end do
end if

write(ifileid,"(a)") "  &END SUBSYS"


!---- &FIST
if (method=="FIST") then
    write(ifileid,"(/,a)") "  #Ref: https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/MM/FORCEFIELD.html"
    write(ifileid,"(a)") "  &MM"
    write(ifileid,"(a)") "    &FORCEFIELD"
    write(ifileid,"(a)") "      &SPLINE"
    write(ifileid,"(a)") "        #EMAX_SPLINE 5.0 #Maximum value of the potential up to which splines will be constructed, in Hartree"
    write(ifileid,"(a)") "        #RCUT_NB 12.0 #Cutoff radius for nonbonded interactions, overrides that specified in potential definitions and is global for all potentials"
    write(ifileid,"(a)") "      &END SPLINE"
    write(ifileid,"(a)") "      IGNORE_MISSING_CRITICAL_PARAMS F"
    write(ifileid,"(a)") "      #PARM_FILE_NAME sobereva.prm"
    write(ifileid,"(a)") "      #PARMTYPE CHM"
    write(ifileid,"(a)") "      #&CHARGE"
    write(ifileid,"(a)") "      #  ATOM X"
    write(ifileid,"(a)") "      #  CHARGE 0.0"
    write(ifileid,"(a)") "      #&END CHARGE"
    write(ifileid,"(a)") "      #&BOND"
    write(ifileid,"(a)") "      #  ATOMS X X"
    write(ifileid,"(a)") "      #  KIND HARMONIC"
    write(ifileid,"(a)") "      #  R0 [angstrom] 0.0"
    write(ifileid,"(a)") "      #  K [angstrom^-2kcalmol] 0.0"
    write(ifileid,"(a)") "      #&END BOND"
    write(ifileid,"(a)") "      #&BEND"
    write(ifileid,"(a)") "      #  ATOMS X X X"
    write(ifileid,"(a)") "      #  KIND HARMONIC"
    write(ifileid,"(a)") "      #  THETA0 [deg] 0.0"
    write(ifileid,"(a)") "      #  K [rad^-2kcalmol] 0.0"
    write(ifileid,"(a)") "      #&END BEND"
    write(ifileid,"(a)") "      #&NONBONDED"
    write(ifileid,"(a)") "      #  &LENNARD-JONES"
    write(ifileid,"(a)") "      #    ATOMS X X"
    write(ifileid,"(a)") "      #    EPSILON [kcalmol]  0.0"
    write(ifileid,"(a)") "      #    SIGMA   [angstrom] 0.0"
    write(ifileid,"(a)") "      #    RCUT    [angstrom] 12.0"
    write(ifileid,"(a)") "      #  &END LENNARD-JONES"
    write(ifileid,"(a)") "      #&END NONBONDED"
    write(ifileid,"(a)") "    &END FORCEFIELD"
    write(ifileid,"(a)") "    #&PRINT"
    write(ifileid,"(a)") "    #  &FF_INFO #Print employed forcefield parameters, missing non-critical parameters, atomic charges, etc."
    write(ifileid,"(a)") "    #    SPLINE_INFO F"
    write(ifileid,"(a)") "    #  &END FF_INFO"
    write(ifileid,"(a)") "    #&END PRINT"
    write(ifileid,"(a)") "    &POISSON"
    write(ifileid,"(a)") "      &EWALD"
    if (PBCdir=="NONE") then
        write(ifileid,"(a)") "        EWALD_TYPE NONE"
    else
        write(ifileid,"(a)") "        EWALD_TYPE SPME"
        ngmax1=ceiling(cellv1_pseudo(1)*nrep1*b2a) !Make Gmax approximately 1 pt/Angstrom in each direction
        ngmax2=ceiling(cellv2_pseudo(2)*nrep2*b2a)
        ngmax3=ceiling(cellv3_pseudo(3)*nrep3*b2a)
        !Only for Ewald Gmax should be odd, while we default to use SPME
        !if (mod(ngmax1,2)==0) ngmax1=ngmax1+1
        !if (mod(ngmax2,2)==0) ngmax2=ngmax2+1
        !if (mod(ngmax3,2)==0) ngmax3=ngmax3+1
        write(ifileid,"(a,3i4,a)") "        GMAX",ngmax1,ngmax2,ngmax3," #Number of grid points for Ewald/SPME in each direction"
    end if
    write(ifileid,"(a)") "      &END EWALD"
    write(ifileid,"(a)") "    &END POISSON"
    write(ifileid,"(a)") "  &END MM"

else

!---- &DFT
write(ifileid,"(/,a)") "  &DFT"
if (method(1:3)/="GFN".and.method/="PM6".and.method/="SCC-DFTB") then
    if (ibas<0) then
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  GTH_BASIS_SETS"
    else if (ibas<=5) then
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  BASIS_MOLOPT"
    else if (ibas==7.or.ibas==8.or.ibas==9) then
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  BASIS_ccGRB_UZH"
    else if (ibas==10.or.ibas==11.or.ibas==12.or.ibas==16) then
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  EMSL_BASIS_SETS"
    else if (ibas==13.or.ibas==14.or.ibas==15.or. (index(method,"ADMM")/=0.and.iGAPW==1) ) then !This file is also needed when pob-DZVP-rev2 is used as ADMM auxiliary basis set
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  BASIS_pob"
    else if (ibas==19) then
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  BASIS_def2_QZVP_RI_ALL"
    else if (ibas==20.or.ibas==21.or.ibas==22) then
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  BASIS_RI_cc-TZ"
    end if
    if (index(method,"ADMM")/=0.and.iGAPW==0) then !Set basis set file for ADMM
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  BASIS_ADMM_UZH"
    end if
    if (iLRIGPW==1) then !Set basis set file for LRIGPW
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  BASIS_LRIGPW_AUXMOLOPT"
    end if
    if (itask==15.or.allocated(XASatm)) then
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  pcseg"
        write(ifileid,"(a)") "    BASIS_SET_FILE_NAME  pcseg-admm"
        if (itask==15) write(ifileid,"(a)") "    AUTO_BASIS RI_XAS MEDIUM #Generate auxiliary basis set for XAS"
    end if
    if (index(method,"MP2")/=0) then
        write(ifileid,"(a)") "    POTENTIAL_FILE_NAME  HF_POTENTIALS"
    else if (index(method,"revDSD-PBEP86")/=0) then
        write(ifileid,"(a)") "    POTENTIAL_FILE_NAME  POTENTIAL_UZH"
    !else if (index(method,"SCAN")/=0.or.index(method,"r2SCAN")/=0.or.index(method,"PBE0")/=0) then
    !    write(ifileid,"(a)") "    POTENTIAL_FILE_NAME  POTENTIAL_UZH"
    else
        write(ifileid,"(a)") "    POTENTIAL_FILE_NAME  POTENTIAL"
    end if
    if (iRIHFX==1) write(ifileid,"(a)") "    AUTO_BASIS RI_HFX MEDIUM"
end if
if (iRIHFX==1) write(ifileid,"(a)") "    SORT_BASIS EXP"
if (ikpoint1==1.and.ikpoint2==1.and.ikpoint3==1) then
    c80tmp=".wfn"
else
    c80tmp=".kp"
end if
if (iHFX==1) then
    write(ifileid,"(a)") "    WFN_RESTART_FILE_NAME "//trim(c200tmp)//"-RESTART"//trim(c80tmp)
else
    write(ifileid,"(a)") "#   WFN_RESTART_FILE_NAME "//trim(c200tmp)//"-RESTART"//trim(c80tmp)
end if
write(ifileid,"(a,i5,a)") "    CHARGE",netchg," #Net charge"
write(ifileid,"(a,i5,a)") "    MULTIPLICITY",multispin," #Spin multiplicity"
if (multispin>1.or.any(kindmag(1:nkind)/=0)) write(ifileid,"(a)") "    UKS"
if (any(efieldvec/=0)) then
    efieldmag=dsqrt(sum(efieldvec**2))
    if (ifPBC==0.or.PBCdir=="NONE") then
        write(ifileid,"(a)") "    &EFIELD"
        write(ifileid,"(a,f12.6)") "      INTENSITY",efieldmag
        write(ifileid,"(a,3f12.6)") "      POLARISATION",efieldvec/efieldmag
        write(ifileid,"(a)") "    &END EFIELD"
    else
        write(ifileid,"(a)") "    &PERIODIC_EFIELD"
        write(ifileid,"(a,f12.6)") "      INTENSITY",efieldmag
        write(ifileid,"(a,3f12.6)") "      POLARISATION",efieldvec/efieldmag
        write(ifileid,"(a)") "    &END PERIODIC_EFIELD"
    end if
end if
if (ikpconvtest==1) then
    write(ifileid,"(a)") "    &KPOINTS"
    write(ifileid,"(a)") "      SCHEME MONKHORST-PACK kp_test"
    write(ifileid,"(a)") "    &END KPOINTS"
else if (ikpoint1/=1.or.ikpoint2/=1.or.ikpoint3/=1) then
    write(ifileid,"(a)") "    &KPOINTS"
    write(ifileid,"(a,3i3)") "      SCHEME MONKHORST-PACK",ikpoint1,ikpoint2,ikpoint3
    write(ifileid,"(a)") "    &END KPOINTS"
end if
if (iDFTplusU==1) write(ifileid,"(a)") "    PLUS_U_METHOD MULLIKEN #The method used in DFT+U. Can also be Lowdin"
if (iTDDFT==1.and.(itask==2.or.itask==3.or.itask==4.or.itask==5.or.itask==6.or.itask==7.or.itask==8)) then !Evaluate force for excited state
    write(ifileid,"(a)") "    &EXCITED_STATES"
    if (itask==2) then !Only calculate force
        write(ifileid,"(a,i5,a)") "      STATE",istate_force," #For which excited state the force will be evaluated. Negative value indicates state following"
    else !State following
        write(ifileid,"(a,i5,a)") "      STATE",-istate_force," #For which excited state the force will be evaluated. Negative value indicates state following"
    end if
    write(ifileid,"(a)") "    &END EXCITED_STATES"
end if

!---- &QS
write(ifileid,"(a)") "    &QS"

if (iLSSCF==1) write(ifileid,"(a)") "      LS_SCF #Use linear scaling self consistent field method"
!Set proper EPS_SCF (default is 1E-5) and EPS_DEFAULT (default is 1E-10)
if (itask==1.and.iTDDFT==0) then !Single point of ground state, does not need high accuracy
    eps_scf=5D-6
    eps_def=1D-11
    if (iwfc==1) then
        eps_scf=1D-6
        eps_def=1D-12
    end if
else if (itask==5.or.itask==9.or.itask==10.or.((itask==2.or.itask==3.or.itask==4.or.itask==7).and.itightopt==1)) then !freq, NMR, polar, force, opt(structure, cell, TS) for subsequent freq purpose
    eps_scf=1D-7
    eps_def=1D-14
else if (itask==6.or.itask==14) then !For faster MD/PIMD, use even looser threshold
    eps_scf=1D-5
    eps_def=1D-10
else !Other tasks involving energy derivative, use marginally tighter convergence
    eps_scf=1D-6
    eps_def=1D-12
end if
write(ifileid,"(a,1PE8.1,a)") "      EPS_DEFAULT",eps_def," #Set all EPS_xxx to values such that the energy will be correct up to this value"

if (iRIHFX==1) then
    write(ifileid,"(a)") "      EPS_PGF_ORB 1E-5"
else
    if (iHFX==1.or.method=="RI-(EXX+RPA)@PBE".or.index(method,"GW")/=0) then
        write(ifileid,"(a)") "      EPS_PGF_ORB 1E-12 #If warning ""Kohn Sham matrix not 100% occupied"" occurs and meantime calculation is unstable, decrease it"
    end if
end if
if (ikpoint1==1.and.ikpoint2==1.and.ikpoint3==1) then !Gamma
    if (itask==6) then !Explicitly add these, allowing users to adapt
        write(ifileid,"(a)") "      EXTRAPOLATION ASPC #Extrapolation for wavefunction during e.g. MD. ASPC is default, PS can also be used"
        write(ifileid,"(a)") "      EXTRAPOLATION_ORDER 3 #Order for PS or ASPC extrapolation. 3 is default"
    end if
else !k-point case, wavefunction extrapolation cannot be uesd, so add a line so that users can uncomment it to use previous density matrix
    if (itask==3.or.itask==4.or.itask==5.or.itask==6.or.itask==7.or.itask==8) then !Tasks involving geometry change
        write(ifileid,"(a)") "#     EXTRAPOLATION USE_PREV_P #Use converged density matrix of last geometry as initial guess"
    end if
end if
if (method=="GFN1-xTB") then
    write(ifileid,"(a)") "      METHOD xTB"
    write(ifileid,"(a)") "      &xTB"
    if (PBCdir/="NONE") write(ifileid,"(a)") "        DO_EWALD T" !Default is Coulomb way to calculate electrostatic interaction
    write(ifileid,"(a)") "        CHECK_ATOMIC_CHARGES F #xTB calculation often crashes without setting this to false"
    write(ifileid,"(a)") "      &END xTB"
else if (method=="GFN2-xTB") then
    write(ifileid,"(a)") "      METHOD xTB"
    write(ifileid,"(a)") "      &xTB"
    write(ifileid,"(a)") "        &TBLITE"
    write(ifileid,"(a)") "          METHOD GFN2"
    write(ifileid,"(a)") "        &END TBLITE"
    write(ifileid,"(a)") "      &END xTB"
else if (method=="GFN0-xTB") then
    write(ifileid,"(a)") "      METHOD xTB"
    write(ifileid,"(a)") "      &xTB"
    write(ifileid,"(a)") "        GFN_TYPE 0"
    write(ifileid,"(a)") "      &END xTB"
else if (method=="PM6") then
    write(ifileid,"(a)") "      METHOD PM6"
    if (ifPBC>0) then
        write(ifileid,"(a)") "      &SE"
        write(ifileid,"(a)") "        PERIODIC EWALD"
        write(ifileid,"(a)") "      &END SE"
    end if
else if (method=="SCC-DFTB") then
    write(ifileid,"(a)") "      METHOD DFTB"
    write(ifileid,"(a)") "      &DFTB"
    write(ifileid,"(a)") "        SELF_CONSISTENT  T"
    write(ifileid,"(a)") "        DISPERSION       T"
    if (PBCdir/="NONE") write(ifileid,"(a)") "        DO_EWALD         T"
    write(ifileid,"(a)") "        &PARAMETER"
    write(ifileid,"(a)") "          PARAM_FILE_PATH  DFTB/scc"
    write(ifileid,"(a)") "          PARAM_FILE_NAME  scc_parameter"
    write(ifileid,"(a)") "          UFF_FORCE_FIELD  uff_table"
    write(ifileid,"(a)") "        &END PARAMETER"
    write(ifileid,"(a)") "      &END DFTB"
else if (iGAPW==1.or.itask==15.or.allocated(XASatm)) then !XAS always needs GAPW
    write(ifileid,"(a)") "      METHOD GAPW"
else !GPW with GTH pseudopotential
    if (iLRIGPW==1) then
        write(ifileid,"(a)") "      METHOD LRIGPW"
        write(ifileid,"(a)") "      &LRIGPW"
        write(ifileid,"(a)") "        LRI_OVERLAP_MATRIX AUTOSELECT #Choose automatically for each pair whether to calculate the inverse or pseudoinverse"
        write(ifileid,"(a)") "      &END LRIGPW"
    end if
end if
if (nCDFTgroup>0) then !Constrained DFT
    write(ifileid,"(a)") "      &CDFT #Set parameters of constrained DFT"
    write(ifileid,"(a)") "        TYPE_OF_CONSTRAINT Hirshfeld #Type of weighting function for constraint"
    write(ifileid,"(a)") "        ATOMIC_CHARGES F #Print atomic CDFT charges during iterations"
    write(ifileid,"(a)",advance="no") "        STRENGTH"
    do igroup=1,nCDFTgroup
        write(ifileid,"(a)",advance="no") " 0"
    end do
    write(ifileid,"(a)") " #Initial Lagrangian multipliers"
    write(ifileid,"(a)",advance="no") "        TARGET"
    do igroup=1,nCDFTgroup
        write(ifileid,"(f10.3)",advance="no") CDFTtarget(igroup)
    end do
    write(ifileid,"(a)") " #Constraint target values"
    do igroup=1,nCDFTgroup
        write(ifileid,"(a,i3)") "        &ATOM_GROUP #CDFT group",igroup
        write(ifileid,"(a)") "          #Indices of the atoms involved in the constraint"
        call outCP2K_LIST(ifileid,CDFTatm(1:CDFTnatm(igroup),igroup),CDFTnatm(igroup),"          ","ATOMS")
        write(ifileid,"(a)") "          #Coefficients of the atoms involved in the constraint"
        write(ifileid,"(a)",advance="no") "          COEFF"
        do idx=1,CDFTnatm(igroup)
            write(ifileid,"(a)",advance="no") " 1"
        end do
        write(ifileid,*)
        if (CDFTtype(igroup)==1) write(ifileid,"(a)") "          CONSTRAINT_TYPE ALPHA #Type of constraint"
        if (CDFTtype(igroup)==2) write(ifileid,"(a)") "          CONSTRAINT_TYPE BETA #Type of constraint"
        if (CDFTtype(igroup)==3) write(ifileid,"(a)") "          CONSTRAINT_TYPE CHARGE #Type of constraint"
        if (CDFTtype(igroup)==4) write(ifileid,"(a)") "          CONSTRAINT_TYPE MAGNETIZATION #Type of constraint"
        write(ifileid,"(a)") "        &END ATOM_GROUP"
    end do
    write(ifileid,"(a)") "        #&DUMMY_ATOMS #No constraint applied but calculate charges"
    write(ifileid,"(a)") "        #  ATOMS ..."
    write(ifileid,"(a)") "        #&END DUMMY_ATOMS"
    write(ifileid,"(a)") "        &OUTER_SCF #Method of optimizing Lagrangian multiplier"
    write(ifileid,"(a)") "          TYPE CDFT_CONSTRAINT #Kind of outer SCF is CDFT"
    write(ifileid,"(a)") "          EPS_SCF 1E-3 #Convergence threshold Lagrangian multiplier"
    if (nCDFTgroup==1) then
        write(ifileid,"(a)") "          OPTIMIZER SECANT"
    else
        write(ifileid,"(a)") "          OPTIMIZER NEWTON_LS #Newton's method with backtracking line search"
    end if
    write(ifileid,"(a)") "          #STEP_SIZE 0.1 #Initial step size used in optimizer"
    write(ifileid,"(a)") "        &END OUTER_SCF"
    write(ifileid,"(a)") "        &HIRSHFELD_CONSTRAINT"
    write(ifileid,"(a)") "          SHAPE_FUNCTION DENSITY #Use atomic density expanded by multiple Gaussians for Hirshfeld partitioning"
    write(ifileid,"(a)") "        &END HIRSHFELD_CONSTRAINT"
    write(ifileid,"(a)") "      &END CDFT"
end if
write(ifileid,"(a)") "    &END QS"

!---- &POISSON
write(ifileid,"(a)") "    &POISSON" !How to deal with electrostatic part
if (method=="PM6") then !Special for semi-empirical
    if (PBCdir/="NONE") then
        write(ifileid,"(a)") "      &EWALD"
        write(ifileid,"(a)") "        &MULTIPOLES"
        write(ifileid,"(a)") "          MAX_MULTIPOLE_EXPANSION QUADRUPOLE"
        write(ifileid,"(a)") "        &END MULTIPOLES"
        write(ifileid,"(a)") "        EWALD_TYPE EWALD"
        !write(ifileid,"(a)") "        ALPHA  0.5" !See e.g. https://github.com/misteliy/cp2k/blob/master/tests/SE/regtest-3/Al2O3.inp
        ngmax1=nint(cellv1_pseudo(1)*nrep1*b2a) !Make gmax approximately 1 pt/Angstrom in each direction
        if (mod(ngmax1,2)==0) ngmax1=ngmax1+1 !Must be odd number
        ngmax2=nint(cellv2_pseudo(2)*nrep2*b2a)
        if (mod(ngmax2,2)==0) ngmax2=ngmax2+1 !Must be odd number
        ngmax3=nint(cellv3_pseudo(3)*nrep3*b2a)
        if (mod(ngmax3,2)==0) ngmax3=ngmax3+1
        write(ifileid,"(a,3i4)") "        GMAX",ngmax1,ngmax2,ngmax3
        write(ifileid,"(a)") "      &END EWALD"
    end if
else if (method=="SCC-DFTB") then !Special for DFTB
    if (PBCdir/="NONE") then
        write(ifileid,"(a)") "      &EWALD"
        write(ifileid,"(a)") "        EWALD_TYPE SPME"
        ngmax1=2*nint(cellv1_pseudo(1)*nrep1*b2a) !Make gmax approximately 1 pt/Angstrom in each direction
        ngmax2=2*nint(cellv2_pseudo(2)*nrep2*b2a)
        ngmax3=2*nint(cellv3_pseudo(3)*nrep3*b2a)
        write(ifileid,"(a,3i4)") "        GMAX",ngmax1,ngmax2,ngmax3
        write(ifileid,"(a)") "      &END EWALD"
    end if
else !Common case
    if (iPSOLVER==1) then
        write(ifileid,"(a)") "      PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics"
    else
        write(ifileid,"(a)") "      PERIODIC "//trim(PBCdir)//" #Direction(s) of PBC for calculating electrostatics"
    end if
    if (iPSOLVER==1) write(ifileid,"(a)") "      PSOLVER PERIODIC #The way to solve Poisson equation"
    if (iPSOLVER==2) write(ifileid,"(a)") "      PSOLVER ANALYTIC #The way to solve Poisson equation"
    if (iPSOLVER==3) write(ifileid,"(a)") "      PSOLVER MT #The way to solve Poisson equation"
    if (iPSOLVER==4) write(ifileid,"(a)") "      PSOLVER WAVELET #The way to solve Poisson equation"
end if
write(ifileid,"(a)") "    &END POISSON"

if (index(method,"ADMM")/=0) then !Use ADMM
    write(ifileid,"(a)") "    &AUXILIARY_DENSITY_MATRIX_METHOD"
    if (iRIHFX==1) then
        write(ifileid,"(a)") "      ADMM_TYPE ADMMS #Best choice for RI-HFX with ADMM"
    else
        write(ifileid,"(a)") "      METHOD BASIS_PROJECTION #Method used for wavefunction fitting"
        if (iTDDFT==1) then
            write(ifileid,"(a)") "      ADMM_PURIFICATION_METHOD NONE #NONE is the only choice for TDDFT with ADMM"
        else if (itask==15) then !It is recommended by XAS_TDP tutorial to use NONE, although purification version seems also to be supported by XAS_TDP
            write(ifileid,"(a)") "      ADMM_PURIFICATION_METHOD NONE"
        else
            if (idiagOT==1) then !Diagonalization cannot be used with purification
                write(ifileid,"(a)") "      ADMM_PURIFICATION_METHOD NONE"
            else !When use OT, the default MO_DIAG is the best
                write(ifileid,"(a)") "      ADMM_PURIFICATION_METHOD MO_DIAG"
            end if
        end if
    end if
    if (iTDDFT==1.and.index(method,"_ADMM")/=0) then !Needed otherwise cannot run. Suggested by https://www.cp2k.org/howto:tddft
        write(ifileid,"(a)") "      EXCH_SCALING_MODEL NONE"
        if (method=="B3LYP_ADMM".or.method=="BHandHLYP_ADMM") then
            write(ifileid,"(a)") "      EXCH_CORRECTION_FUNC BECKE88X"
        else
            write(ifileid,"(a)") "      EXCH_CORRECTION_FUNC PBEX"
        end if
    end if
    write(ifileid,"(a)") "    &END AUXILIARY_DENSITY_MATRIX_METHOD"
end if

!---- &XC
if (method=="PM6".or.method(1:3)=="GFN".or.method=="SCC-DFTB") goto 100
write(ifileid,"(a)") "    &XC"
if (index(method,"LIBXC")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    if (method=="B97M-rV_LIBXC") then !Non-separable XC
        write(ifileid,"(a)") "        &MGGA_XC_B97M_V"
        write(ifileid,"(a)") "        &END MGGA_XC_B97M_V"
    else !X-C separable
        if (method=="MN15L_LIBXC") then
            write(ifileid,"(a)") "        &MGGA_X_MN15_L"
            write(ifileid,"(a)") "        &END MGGA_X_MN15_L"
            write(ifileid,"(a)") "        &MGGA_C_MN15_L"
            write(ifileid,"(a)") "        &END MGGA_C_MN15_L"
        else if (method=="SCAN_LIBXC") then
            write(ifileid,"(a)") "        &MGGA_X_SCAN"
            write(ifileid,"(a)") "        &END MGGA_X_SCAN"
            write(ifileid,"(a)") "        &MGGA_C_SCAN"
            write(ifileid,"(a)") "        &END MGGA_C_SCAN"
        else if (method=="r2SCAN_LIBXC") then
            write(ifileid,"(a)") "        &MGGA_X_R2SCAN"
            write(ifileid,"(a)") "        &END MGGA_X_R2SCAN"
            write(ifileid,"(a)") "        &MGGA_C_R2SCAN"
            write(ifileid,"(a)") "        &END MGGA_C_R2SCAN"
        else if (method=="RPBE_LIBXC") then
            write(ifileid,"(a)") "        &GGA_X_RPBE"
            write(ifileid,"(a)") "        &END GGA_X_RPBE"
            write(ifileid,"(a)") "        &GGA_C_PBE"
            write(ifileid,"(a)") "        &END GGA_C_PBE"
        else if (method=="TPSS_LIBXC") then
            write(ifileid,"(a)") "        &MGGA_X_TPSS"
            write(ifileid,"(a)") "        &END MGGA_X_TPSS"
            write(ifileid,"(a)") "        &MGGA_C_TPSS"
            write(ifileid,"(a)") "        &END MGGA_C_TPSS"
        else if (method=="revTPSS_LIBXC") then
            write(ifileid,"(a)") "        &MGGA_X_REVTPSS"
            write(ifileid,"(a)") "        &END MGGA_X_REVTPSS"
            write(ifileid,"(a)") "        &MGGA_C_REVTPSS"
            write(ifileid,"(a)") "        &END MGGA_C_REVTPSS"
        else if (method=="HLE17_LIBXC") then
            write(ifileid,"(a)") "        &MGGA_XC_HLE17"
            write(ifileid,"(a)") "        &END MGGA_XC_HLE17"
        end if
    end if
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (method(1:3)=="GFN".or.method=="PM6".or.method=="SCC-DFTB") then
    continue
else if (index(method,"PBE0")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &PBE"
    write(ifileid,"(a)") "          SCALE_X 0.75"
    write(ifileid,"(a)") "          SCALE_C 1.0"
    write(ifileid,"(a)") "        &END PBE"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"PBEh")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &PBE"
    write(ifileid,"(a,f10.5)") "          SCALE_X",(100-PBEh_HFX)/100
    write(ifileid,"(a)") "        &END PBE"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"B3LYP")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &LYP"
    write(ifileid,"(a)") "          SCALE_C 0.81"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "        &BECKE88"
    write(ifileid,"(a)") "          SCALE_X 0.72"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "        &VWN"
    write(ifileid,"(a)") "          FUNCTIONAL_TYPE VWN3 #Gaussian's B3LYP definition"
    write(ifileid,"(a)") "          SCALE_C 0.19"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "        &XALPHA"
    write(ifileid,"(a)") "          SCALE_X 0.08"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"BHandHLYP")/=0) then !Including case of GW@BHandHLYP_ADMM
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &HYB_GGA_XC_BHANDHLYP"
    write(ifileid,"(a)") "        &END HYB_GGA_XC_BHANDHLYP"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"GW@MN15L")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &MGGA_X_MN15_L"
    write(ifileid,"(a)") "        &END MGGA_X_MN15_L"
    write(ifileid,"(a)") "        &MGGA_C_MN15_L"
    write(ifileid,"(a)") "        &END MGGA_C_MN15_L"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"M06-2X")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &HYB_MGGA_X_M06_2X"
    write(ifileid,"(a)") "        &END HYB_MGGA_X_M06_2X"
    write(ifileid,"(a)") "        &MGGA_C_M06_2X"
    write(ifileid,"(a)") "        &END MGGA_C_M06_2X"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (method=="revPBE".or.method=="PBEsol") then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &PBE"
    if (method=="revPBE") write(ifileid,"(a)") "          PARAMETRIZATION REVPBE"
    if (method=="PBEsol") write(ifileid,"(a)") "          PARAMETRIZATION PBESOL"
    write(ifileid,"(a)") "        &END PBE"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"HSE")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &XWPBE"
    write(ifileid,"(a)") "          SCALE_X -0.25"
    write(ifileid,"(a)") "          SCALE_X0 1.0"
    write(ifileid,"(a)") "          OMEGA 0.11"
    write(ifileid,"(a)") "        &END XWPBE"
    write(ifileid,"(a)") "        &PBE"
    write(ifileid,"(a)") "          SCALE_X 0.0"
    write(ifileid,"(a)") "          SCALE_C 1.0"
    write(ifileid,"(a)") "        &END PBE"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"MP2")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL NONE"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (method=="RI-(EXX+RPA)@PBE") then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL PBE"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"B2PLYP")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &LYP"
    write(ifileid,"(a)") "          SCALE_C  0.73"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "        &BECKE88"
    write(ifileid,"(a)") "          SCALE_X  0.47"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"B2GP-PLYP")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &LYP"
    write(ifileid,"(a)") "          SCALE_C  0.64"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "        &BECKE88"
    write(ifileid,"(a)") "          SCALE_X  0.35"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"DSD-BLYP")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &LYP"
    write(ifileid,"(a)") "          SCALE_C  0.69"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "        &BECKE88"
    write(ifileid,"(a)") "          SCALE_X  0.31"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else if (index(method,"revDSD-PBEP86")/=0) then
    write(ifileid,"(a)") "      &XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &GGA_X_PBE"
    write(ifileid,"(a)") "          SCALE 0.31"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "        &P86C"
    write(ifileid,"(a)") "          SCALE_C 0.4296"
    write(ifileid,"(a)") "        &END"
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
else !Common native GGA functionals
    write(ifileid,"(a)") "      &XC_FUNCTIONAL "//trim(method)
    write(ifileid,"(a)") "      &END XC_FUNCTIONAL"
end if

!Coulomb truncation radius (CUTOFF_RADIUS) may be used later, now calculate (recorded in Angstrom!)
!Set it to 1/2.01 of shortest Cartesian length for small cell
cellv1_tmp=cellv1;cellv2_tmp=cellv2;cellv3_tmp=cellv3
i1=ifdoPBCx; i2=ifdoPBCy; i3=ifdoPBCz
cellv1=cellv1_pseudo*nrep1
cellv2=cellv2_pseudo*nrep2
cellv3=cellv3_pseudo*nrep3
ifdoPBCx=1;ifdoPBCy=1;ifdoPBCz=1
call cellplane_spacing(0,0,1,d1)
call cellplane_spacing(0,1,0,d2)
call cellplane_spacing(1,0,0,d3)
trunc_rad=min(min(d1,d2),d3)/2.01D0 !Maximum truncation radius in Angstrom
trunc_rad_kp=min(min(d1*ikpoint1,d2*ikpoint2),d3*ikpoint3)/2.01D0 !Maximum truncation radius in Angstrom with considering k-points
ifdoPBCx=i1; ifdoPBCy=i2; ifdoPBCz=i3
cellv1=cellv1_tmp;cellv2=cellv2_tmp;cellv3=cellv3_tmp

!HF part for hybrid functionals
if (iHFX==1) then !HFX potential
    write(ifileid,"(a)") "      &HF"
    if (index(method,"PBE0")/=0.or.index(method,"HSE")/=0) write(ifileid,"(a)") "        FRACTION 0.25 #HF composition"
    if (index(method,"B3LYP")/=0) write(ifileid,"(a)") "        FRACTION 0.2 #HF composition"
    if (index(method,"BHandHLYP")/=0) write(ifileid,"(a)") "        FRACTION 0.5 #HF composition"
    if (index(method,"PBEh")/=0) write(ifileid,"(a,f10.5,' #HF composition')") "        FRACTION",PBEh_HFx/100
    if (index(method,"M06-2X")/=0) write(ifileid,"(a)") "        FRACTION 0.54 #HF composition"
    if (index(method,"MP2")/=0) write(ifileid,"(a)") "        FRACTION 1.0 #HF composition"
    if (index(method,"B2PLYP")/=0) write(ifileid,"(a)") "        FRACTION 0.53 #HF composition"
    if (index(method,"B2GP-PLYP")/=0) write(ifileid,"(a)") "        FRACTION 0.65 #HF composition"
    if (index(method,"DSD-BLYP")/=0) write(ifileid,"(a)") "        FRACTION 0.69 #HF composition"
    if (index(method,"revDSD-PBEP86")/=0) write(ifileid,"(a)") "        FRACTION 0.69 #HF composition"
    write(ifileid,"(a)") "        &SCREENING"
    if (iwfc==1) then !For wavefunction-based correlation calculation, use tighter threshold for screening
        write(ifileid,"(a)") "          EPS_SCHWARZ 1E-7 #The larger the value, the lower the cost and lower the accuracy"
    else !Hybrid functionals
        write(ifileid,"(a)") "          EPS_SCHWARZ 1E-6 #The larger the value, the lower the cost and lower the accuracy"
    end if
    if (iRIHFX==0) write(ifileid,"(a)") "          SCREEN_ON_INITIAL_P T #Screening ERI based on initial density matrix, need to provide wavefunction restart file"
    write(ifileid,"(a)") "        &END SCREENING"
    if (index(method,"HSE")/=0) then
        write(ifileid,"(a)") "        &INTERACTION_POTENTIAL"
        write(ifileid,"(a)") "          POTENTIAL_TYPE SHORTRANGE"
        write(ifileid,"(a)") "          OMEGA 0.11"
        write(ifileid,"(a)") "        &END INTERACTION_POTENTIAL"
    else
        if (PBCdir/="NONE") then !PBC system needs Coulomb truncation for common hybrid functionals. By default the Coulomb interaction is not truncated
            write(ifileid,"(a)") "        &INTERACTION_POTENTIAL"
            write(ifileid,"(a)") "          POTENTIAL_TYPE TRUNCATED"
            if (ikpoint1==1.and.ikpoint2==1.and.ikpoint3==1) then !Gamma point
                if (trunc_rad>6) then !If half of shortest box length is larger than 6 Angstrom, simply use 6, this is usually adequate
                    write(ifileid,"(a,f8.4)") "          CUTOFF_RADIUS 6.0 #Cutoff radius (Angstrom) for truncated 1/r Coulomb operator"
                else
                    write(ifileid,"(a,f8.4,a)") "          CUTOFF_RADIUS",trunc_rad," #Cutoff radius (Angstrom) for truncated 1/r Coulomb operator"
                end if
            else !RI-HFXk case
                if (trunc_rad_kp>6) then !If half of shortest box length is larger than 6 Angstrom, simply use 6, this is usually adequate
                    write(ifileid,"(a,f8.4)") "          CUTOFF_RADIUS 6.0 #Cutoff radius (Angstrom) for truncated 1/r Coulomb operator"
                else
                    write(ifileid,"(a,f8.4,a)") "          CUTOFF_RADIUS",trunc_rad_kp," #Cutoff radius (Angstrom) for truncated 1/r Coulomb operator"
                end if
            end if
            write(ifileid,"(a)") "        &END INTERACTION_POTENTIAL"
        end if
    end if
    if (iRIHFX==0) then !Calculate 4-center ERI, takes lots of memory
        write(ifileid,"(a)") "        &MEMORY"
        write(ifileid,"(a)") "          MAX_MEMORY 3000 #Memory(MB) per MPI process for calculating HF exchange"
        !Scaling factor to scale EPS_SCHWARZ. Storage threshold for compression will be EPS_SCHWARZ*EPS_STORAGE_SCALING
        write(ifileid,"(a)") "          EPS_STORAGE_SCALING 0.1"
        write(ifileid,"(a)") "        &END MEMORY"
        write(ifileid,"(a)") "      &END HF"
    else
        write(ifileid,"(a)") "        &RI #Activate and set RI-HFX"
        write(ifileid,"(a)") "          KP_NGROUPS 1 #Number of MPI subgroup. Larger leads to evidently faster calculation but takes more memory. Total MPI ranks must be divisible by this"
        if (ikpoint1==1.and.ikpoint2==1.and.ikpoint3==1) then !Gamma point RI-HFX
            if (ifPBC==0.or.PBCdir=="NONE") then !Isolated. Full Coulomb
                write(ifileid,"(a)") "          RI_METRIC HFX #Default, using same Coulomb operator as &INTERACTION_POTENTIAL for RI metric"
            else !Periodic, this setting is robust and efficient. When calculating gamma point, it is unlikely cell is very small, so 1.5 angstrom is always applicable
                write(ifileid,"(a)") "          RI_METRIC TRUNCATED #Type of RI metric operator"
                write(ifileid,"(a)") "          CUTOFF_RADIUS 1.5 #Cutoff radius (in Angstrom) for truncated RI Coulomb operator"
            end if
        else !RI-HFXk, this setting is saftest and recommended
            write(ifileid,"(a)") "          RI_METRIC HFX #Default, using same Coulomb operator as &INTERACTION_POTENTIAL for RI metric"
        end if
        write(ifileid,"(a)") "          KP_USE_DELTA_P T #Using incremental Fock, set to F if SCF is difficult to converge"
        write(ifileid,"(a)") "          EPS_PGF_ORB 1E-5 #Set precision of integral tensors, the default 1E-5 is usually fine enough"
        !write(ifileid,"(a)") "          MEMORY_CUT 3 #Memory reduction factor. Larger saves more memory but more time-consuming"
        write(ifileid,"(a)") "        &END RI"
        write(ifileid,"(a)") "      &END HF"
    end if
end if

!Wavefunction-based correlation part
if (iwfc==1) then
    write(ifileid,"(a)") "      &WF_CORRELATION"
    if (method=="RI-(EXX+RPA)@PBE".or.index(method,"GW")/=0) then
        write(ifileid,"(a)") "        &RI_RPA"
        if (index(method,"GW")/=0) then
            if (ifPBC==0.or.PBCdir=="NONE") then
                write(ifileid,"(a)") "          QUADRATURE_POINTS  60  #Number of quadrature points for the numerical integration in the GW"
            else !Use larger value for PBC system, because QUADRATURE_POINTS doesn't affect cost significantly, so using larger value to ensure numerical accuracy
                write(ifileid,"(a)") "          QUADRATURE_POINTS  100  #Number of quadrature points for the numerical integration in the GW"
            end if
        else
            write(ifileid,"(a)") "          QUADRATURE_POINTS  10  #Number of quadrature points for the numerical integration in the RI-RPA method"
            write(ifileid,"(a)") "          MINIMAX  #Use Minimax quadrature scheme for performing the numerical integration"
        end if
        write(ifileid,"(a)") "          &HF"
        write(ifileid,"(a)") "            FRACTION 1.0"
        write(ifileid,"(a)") "            &SCREENING"
        write(ifileid,"(a)") "              EPS_SCHWARZ 1E-7"
        write(ifileid,"(a)") "              SCREEN_ON_INITIAL_P F"
        write(ifileid,"(a)") "            &END SCREENING"
        if (ifPBC/=0) then !PBC system needs Coulomb truncation for evaluating HF exchange of RPA energy
            write(ifileid,"(a)") "            &INTERACTION_POTENTIAL"
            write(ifileid,"(a)") "              POTENTIAL_TYPE TRUNCATED"
            if (trunc_rad>6) then !If half of shortest box length is larger than 6 Angstrom, simply use 6, this is usually adequate
                write(ifileid,"(a,f8.4)") "              CUTOFF_RADIUS 6.0 #Cutoff radius (Angstrom) for truncated 1/r Coulomb operator"
            else
                write(ifileid,"(a,f8.4,a)") "              CUTOFF_RADIUS",trunc_rad," #Cutoff radius (Angstrom) for truncated 1/r Coulomb operator"
            end if
            write(ifileid,"(a)") "            &END INTERACTION_POTENTIAL"
        end if
        write(ifileid,"(a)") "          &END HF"
        if (index(method,"GW")/=0) then !Set GW
            write(ifileid,"(a)") "          &GW"
            write(ifileid,"(a)") "            CORR_MOS_OCC   5  #Number of occupied MOs whose energies are to be corrected"
            write(ifileid,"(a)") "            CORR_MOS_VIRT  5  #Number of virtual MOs whose energies are to be corrected"
            !write(ifileid,"(a)") "            ANALYTIC_CONTINUATION PADE  #Type of analytic continuation for the self energy to be used"
            !write(ifileid,"(a)") "            NUMB_POLES            2"
            !write(ifileid,"(a)") "            CROSSING_SEARCH       Z_SHOT"
            !write(ifileid,"(a)") "            FERMI_LEVEL_OFFSET    2.0E-2"
            write(ifileid,"(a,i3,a)") "            EV_GW_ITER",niter_evGW,"  #Maximum number of iterations for eigenvalue self-consistency of evGW"
            write(ifileid,"(a,i3,a)") "            SC_GW0_ITER",niter_scGW0,"   #Maximum number of iterations for GW0 self-consistency of scGW0"
            write(ifileid,"(a)") "            UPDATE_XC_ENERGY F  #If total energy will be corrected using exact exchange and the RPA correlation energy"
            if (ifPBC==0.or.PBCdir=="NONE") then
                write(ifileid,"(a)") "            RI_SIGMA_X T  #If exchange self-energy will be calculated approximatively with RI"
            else
                write(ifileid,"(a)") "            RI_SIGMA_X F  #If exchange self-energy will be calculated approximatively with RI"
            end if
            if (PBCdir=="XYZ") write(ifileid,"(a)") "            PERIODIC_CORRECTION F  #If periodic correction scheme will be used"
            if (ikpoint1>1.or.ikpoint2>1.or.ikpoint3>1) write(ifileid,"(a,3i3,a)") "            KPOINTS_SELF_ENERGY",ikpoint1,ikpoint2,ikpoint3,"  #Number of k-points for self-energy"
            write(ifileid,"(a)") "          &END GW"
            !write(ifileid,"(a)") "          SIZE_FREQ_INTEG_GROUP  1  #Number of processes for computing each integration point, must be multiple of GROUP_SIZE in &WF_CORRELATION"
        end if
        write(ifileid,"(a)") "        &END RI_RPA"
    else !Other case, all involve MP2
        write(ifileid,"(a)") "        &RI_MP2"
        write(ifileid,"(a)") "        &END RI_MP2"
        if (index(method,"SCS-MP2")/=0) then
            write(ifileid,"(a)") "        SCALE_S 1.2"
            write(ifileid,"(a)") "        SCALE_T 0.3333333"
        else if (index(method,"B2PLYP")/=0) then
            write(ifileid,"(a)") "        SCALE_S 0.27"
            write(ifileid,"(a)") "        SCALE_T 0.27"
        else if (index(method,"B2GP-PLYP")/=0) then
            write(ifileid,"(a)") "        SCALE_S 0.36"
            write(ifileid,"(a)") "        SCALE_T 0.36"
        else if (index(method,"DSD-BLYP")/=0) then
            write(ifileid,"(a)") "        SCALE_S 0.46"
            write(ifileid,"(a)") "        SCALE_T 0.37"
        else if (index(method,"revDSD-PBEP86")/=0) then
            write(ifileid,"(a)") "        SCALE_S 0.5785"
            write(ifileid,"(a)") "        SCALE_T 0.0799"
        end if
    end if
    write(ifileid,"(a)") "        &INTEGRALS"
    write(ifileid,"(a)") "          &WFC_GPW"
    write(ifileid,"(a)") "            CUTOFF      300"
    write(ifileid,"(a)") "            REL_CUTOFF  50"
    write(ifileid,"(a)") "            EPS_FILTER  1E-12"
    write(ifileid,"(a)") "            EPS_GRID    1E-8"
    write(ifileid,"(a)") "          &END WFC_GPW"
    write(ifileid,"(a)") "        &END INTEGRALS"
    write(ifileid,"(a)") "        MEMORY    3000 #Maximum allowed total memory usage (MB) during wavefunction-based correlation"
    write(ifileid,"(a)") "        GROUP_SIZE  1 #Default. Also known as NUMBER_PROC"
    write(ifileid,"(a)") "      &END WF_CORRELATION"
end if

!--- Dispersion correction
if (idispcorr>0.or.method=="BEEFVDW") then
    write(ifileid,"(a)") "      &VDW_POTENTIAL"
    if (idispcorr==1.or.idispcorr==2.or.idispcorr==3) then
        write(ifileid,"(a)") "        POTENTIAL_TYPE PAIR_POTENTIAL"
        write(ifileid,"(a)") "        &PAIR_POTENTIAL"
        if (idispcorr/=3) write(ifileid,"(a)") "          PARAMETER_FILE_NAME dftd3.dat"
        if (idispcorr==1) write(ifileid,"(a)") "          TYPE DFTD3"
        if (idispcorr==2) write(ifileid,"(a)") "          TYPE DFTD3(BJ)"
        if (idispcorr==3) write(ifileid,"(a)") "          TYPE DFTD4"
        !See qs_dispersion_pairpot.F on how to write functional name
        !Special cases are explicitly list here, used to change name, remove RI-, remove _LIBXC, etc.
        if (method=="BP") then
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL BP86"
        else if (method=="MN15L_LIBXC") then !i.e. Remove _LIBXC suffix
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL MN15L"
        else if (method=="SCAN_LIBXC") then !i.e. Remove _LIBXC suffix
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL SCAN"
        else if (method=="RPBE_LIBXC") then !i.e. Remove _LIBXC suffix
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL RPBE"
        else if (method=="TPSS_LIBXC") then !i.e. Remove _LIBXC suffix
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL TPSS"
        else if (method=="revTPSS_LIBXC") then !i.e. Remove _LIBXC suffix
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL revTPSS"
        else if (method=="HLE17_LIBXC") then !i.e. Remove _LIBXC suffix
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL HLE17"
        else if (index(method,"B2PLYP")/=0) then
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL B2PLYP"
        else if (index(method,"B2GP-PLYP")/=0) then
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL B2GPPLYP"
        else if (index(method,"DSD-BLYP")/=0) then
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL DSD-BLYP"
        else if (index(method,"revDSD-PBEP86")/=0) then
            write(ifileid,"(a)") "          D3BJ_SCALING 0.4377,0,0,5.5 #s6,a1,s8,a2"
        else
            c80tmp=trim(method)
            ipos=index(c80tmp,"_ADMM")
            if (ipos/=0) c80tmp(ipos:ipos+4)="" !Remove _ADMM suffix
            write(ifileid,"(a)") "          REFERENCE_FUNCTIONAL "//trim(c80tmp)
        end if
        !write(ifileid,"(a)") "          R_CUTOFF 10.5835442" !Default DFT-D potential range, cutoff will be 2 times this value
        if (idispcorr/=3) write(ifileid,"(a)") "          #CALCULATE_C9_TERM T #Calculate C9-related three-body term, more accurate for large system"
        write(ifileid,"(a)") "        &END PAIR_POTENTIAL"
    else if (idispcorr==5.or.method=="BEEFVDW") then         
        write(ifileid,"(a)") "        POTENTIAL_TYPE NON_LOCAL"
        write(ifileid,"(a)") "        &NON_LOCAL"
        if (idispcorr==5) then
            write(ifileid,"(a)") "          TYPE RVV10"
            if (method=="B97M-rV_LIBXC") then !See: Ab initio molecular dynamics simulations of liquid water using high quality meta-GGA functionals
                write(ifileid,"(a)") "          PARAMETERS 6.0 0.01"
            else
                write(ifileid,"(a)") "    #The default rVV10 b and C parameters are given below. They should be replaced by proper values for current XC functional"
                write(ifileid,"(a)") "          PARAMETERS 6.3 9.3E-3"
            end if
            write(ifileid,"(a)") "          KERNEL_FILE_NAME rVV10_kernel_table.dat"
        else if (method=="BEEFVDW") then
            write(ifileid,"(a)") "          TYPE LMKLL"
            write(ifileid,"(a)") "          KERNEL_FILE_NAME vdW_kernel_table.dat"
        end if
        write(ifileid,"(a)") "        &END NON_LOCAL"
    end if
    write(ifileid,"(a)") "      &END VDW_POTENTIAL"
end if
if (ifineXCgrid==1) then
    write(ifileid,"(a)") "      &XC_GRID"
    write(ifileid,"(a)") "        USE_FINER_GRID T #Use finer grid for calculating XC"
    write(ifileid,"(a)") "      &END XC_GRID"
end if
write(ifileid,"(a)") "    &END XC"
100 continue

!--- &MGRID
if (method(1:3)=="GFN".or.method=="PM6".or.method=="SCC-DFTB") then
    continue !Semi-empirical methods do not need to set CUTOFF
else
    write(ifileid,"(a)") "    &MGRID"
    if (iconvtest==1) then
        write(ifileid,"(a)") "      CUTOFF LT_cutoff"
        write(ifileid,"(a)") "      REL_CUTOFF LT_rel_cutoff"
    else
        write(ifileid,"(a,i5)") "      CUTOFF",CUTOFF
        write(ifileid,"(a,i4)") "      REL_CUTOFF",REL_CUTOFF
    end if
    if (ibas==3.or.ibas==4.or.ibas==5) write(ifileid,"(a)") "      NGRIDS 5 #The number of multigrids to use. 5 is optimal for MOLOPT-GTH basis sets"
    write(ifileid,"(a)") "    &END MGRID"
end if

!--- &SCF or &LC_SCF
if (iLSSCF==1) then !&LS_SCF
    write(ifileid,"(a)") "    &LS_SCF"
    write(ifileid,"(a)") "      PURIFICATION_METHOD TRS4  #Scheme used to purify Kohn-Sham matrix into density matrix"
    write(ifileid,"(a)") "      #DYNAMIC_THRESHOLD T  #Should the threshold for the purification be chosen dynamically"
    write(ifileid,"(a)") "      EPS_FILTER 1E-7  #Threshold used to determine sparsity and thus speed and accuracy"
    write(ifileid,"(a)") "      EPS_SCF    5E-6  #Target accuracy for SCF convergence in terms of change of total energy per electron"
    write(ifileid,"(a)") "      MAX_SCF 40  #Maximum number of SCF iteration to be performed"
    write(ifileid,"(a)") "      S_PRECONDITIONER ATOMIC  #Preconditions S with some appropriate form. The default ATOMIC is suitable for most cases"
    write(ifileid,"(a)") "      #MIXING_FRACTION 0.45  #Fraction of mixing new density matrix. A value smaller than the default 0.45 may stablize SCF convergence"
    write(ifileid,"(a)") "      #MU -0.15  #Chemical potential in a.u., does not need to set if using TRS4"
    write(ifileid,"(a)") "    &END LS_SCF"
else !&SCF
    write(ifileid,"(a)") "    &SCF"
    if (method/="GFN0-xTB") then !GFN0-xTB doesn't involve SCF iterations
        if (iconvtest==1) then
            write(ifileid,"(a)") "      MAX_SCF 1"
        else
            if (idiagOT==1) then !Diagonalization
                write(ifileid,"(a)") "      MAX_SCF 128"
            else !OT, usually use more cycles
                if (iouterSCF==0) write(ifileid,"(a)") "      MAX_SCF 200 #Should be set to a small value (e.g. 20) if enabling outer SCF"
                if (iouterSCF==1) write(ifileid,"(a)") "      MAX_SCF 25 #Maximum number of steps of inner SCF"
            end if
        end if
        if (imixing==1) then
            write(ifileid,"(a)") "      MAX_DIIS 7 #Maximum number of DIIS vectors to be used" !The default 4 is too small
            write(ifileid,"(a)") "      EPS_DIIS 0.3 #Threshold on the convergence to start using DIAG/DIIS" !The default 0.1 is too small
        end if
        if (iouterSCF==0) then
            write(ifileid,"(a,1PE8.1,a)") "      EPS_SCF",eps_scf," #Convergence threshold of density matrix during SCF"
        else if (iouterSCF==1) then
            write(ifileid,"(a,1PE8.1,a)") "      EPS_SCF",eps_scf," #Convergence threshold of density matrix of inner SCF"
        end if
        !if (method=="GFN1-xTB".or.method=="PM6") write(ifileid,"(a)") "      SCF_GUESS MOPAC" !Seems they benefit from this
        if (iHFX==1) then
            write(ifileid,"(a)") "      SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess"
        else
            write(ifileid,"(a)") "#     SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess"
        end if
        write(ifileid,"(a)") "#     IGNORE_CONVERGENCE_FAILURE #Continue calculation even if SCF not converged, works for version >= 2024.1"
        if (idiagOT==1) then
            write(ifileid,"(a)") "      &DIAGONALIZATION"
            write(ifileid,"(a)") "        ALGORITHM STANDARD #Algorithm for diagonalization"
            write(ifileid,"(a)") "      &END DIAGONALIZATION"
        else if (idiagOT==2) then
            write(ifileid,"(a)") "      &OT"
            if (method(1:3)=="GFN".or.method=="PM6".or.method=="SCC-DFTB") then !Semi-empirical cannot use FULL_KINETIC. For large system FULL_SINGLE_INVERSE is the only good choice
                write(ifileid,"(a)") "        PRECONDITIONER FULL_SINGLE_INVERSE"
            else
                if (ncenter<300.or.iGAPW==1.or.itask==15) then !GAPW (XAS task implies it) should use FULL_ALL even if the system is large, because it converges much better than FULL_KINETIC according to my test and suggestion by Hutter 
                    write(ifileid,"(a)") "        PRECONDITIONER FULL_ALL #Usually best but expensive for large system. Cheaper: FULL_SINGLE_INVERSE and FULL_KINETIC"
                else !GPW for large systems, using FULL_ALL will cause too high cost at the first step
                    write(ifileid,"(a)") "        PRECONDITIONER FULL_KINETIC #FULL_SINGLE_INVERSE is also worth to try. FULL_ALL is better but quite expensive for large system"
                end if
            end if
            write(ifileid,"(a)") "        MINIMIZER DIIS #CG is worth to consider in difficult cases" !BROYDEN in fact can also be used, but quite poor!
            write(ifileid,"(a)") "        LINESEARCH 2PNT #1D line search algorithm for CG. 2PNT is default. 3PNT is more expensive but may be better. GOLD is best but very expensive"
            if (nCDFTgroup>0.or.ibas==13.or.ibas==14.or.ibas==15) then
                write(ifileid,"(a)") "        ALGORITHM IRAC #Algorithm of OT. Can be STRICT (default) or IRAC" !For CDFT, this is much easier to converge than the default STRICT when lambda>0; For pob, often IRAC converges more smoothly
            else
                write(ifileid,"(a)") "        ALGORITHM STRICT #Algorithm of OT. Can be STRICT (default) or IRAC"
            end if
            write(ifileid,"(a)") "      &END OT"
            if (iouterSCF==0) then
                write(ifileid,"(a)") "      #Uncomment following lines can enable outer SCF, important for difficult convergence case"
                write(ifileid,"(a)") "      #&OUTER_SCF"
                write(ifileid,"(a)") "      #  MAX_SCF 20 #Maximum number of steps of outer SCF"
                write(ifileid,"(a,1PE8.1,a)") "      #  EPS_SCF",eps_scf," #Convergence threshold of outer SCF"
                write(ifileid,"(a)") "      #&END OUTER_SCF"
            else if (iouterSCF==1) then
                write(ifileid,"(a)") "      &OUTER_SCF"
                write(ifileid,"(a)") "        MAX_SCF 20 #Maximum number of steps of outer SCF"
                write(ifileid,"(a,1PE8.1,a)") "        EPS_SCF",eps_scf," #Convergence threshold of outer SCF"
                write(ifileid,"(a)") "      &END OUTER_SCF"
            end if
        end if
    end if
    !Case of using diagonalization
    if (idiagOT==1) then
        if (method/="GFN0-xTB") then !GFN0-xTB doesn't involve SCF iterations
            !--- &SCF \ &MIXING
            write(ifileid,"(a)") "      &MIXING #How to mix old and new density matrices"
            if (imixing==1) then !PM6 and only use this
                write(ifileid,"(a)") "        METHOD DIRECT_P_MIXING"
                if (method=="GFN2-xTB") then !The best choice for this method
                    write(ifileid,"(a)") "        ALPHA 0.2 #Default. Mixing 20% of new density matrix with the old one"
                else
                    write(ifileid,"(a)") "        ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one"
                end if
            else if (imixing==2) then
                write(ifileid,"(a)") "        METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative"
                write(ifileid,"(a)") "        ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one"
                write(ifileid,"(a)") "        NBROYDEN 8 #Default is 4. Number of previous steps stored for the actual mixing scheme" !Equivalent to NBUFFER
            else if (imixing==3) then
                write(ifileid,"(a)") "        METHOD PULAY_MIXING #BROYDEN_MIXING is also a good alternative"
                write(ifileid,"(a)") "        NPULAY 8 #Default is 4. Number of previous steps stored for the actual mixing scheme" !Equivalent to NBUFFER
            end if
            write(ifileid,"(a)") "      &END MIXING"
        end if
        !--- &SCF \ &SMEAR
        if (ismear==1) then
            write(ifileid,"(a)") "      &SMEAR"
            write(ifileid,"(a)") "        METHOD FERMI_DIRAC" !Can also be ENERGY_WINDOW, LIST
            write(ifileid,"(a)") "        ELECTRONIC_TEMPERATURE 300 #Electronic temperature of Fermi-Dirac smearing in K"
            write(ifileid,"(a)") "      &END SMEAR"
        end if
        if (nADDED_MOS/=0) then
            if (multispin>1.or.any(kindmag(1:nkind)/=0)) then
                write(ifileid,"(a,i6,i6,a)") "      ADDED_MOS",nADDED_MOS,nADDED_MOS," #Number of virtual MOs to solve for alpha and beta spins"
            else
                write(ifileid,"(a,i6,a)") "      ADDED_MOS",nADDED_MOS," #Number of virtual MOs to solve"
            end if
        end if
    end if
    !--- &SCF \ &PRINT
    write(ifileid,"(a)") "      &PRINT"
    if (itask==5.or.itask==6.or.itask==14) then !Freq, MD, PIMD
        write(ifileid,"(a)") "        &RESTART OFF #Do not generate wfn file to suppress meaningless I/O cost"
        write(ifileid,"(a)") "        &END RESTART"
    else
        write(ifileid,"(a)") "        &RESTART #Note: Use ""&RESTART OFF"" can prevent generating .wfn file"
        write(ifileid,"(a)") "          BACKUP_COPIES 0 #Maximum number of backup copies of wfn file. 0 means never"
        write(ifileid,"(a)") "        &END RESTART"
    end if
    if (itask==6.and.inoSCFinfo==1) then !MD
        write(ifileid,"(a)") "        &PROGRAM_RUN_INFO"
        write(ifileid,"(a)") "          &EACH"
        write(ifileid,"(a)") "            MD 0 #Frequency of printing SCF process during MD. 0 means never"
        write(ifileid,"(a)") "          &END EACH"
        write(ifileid,"(a)") "        &END PROGRAM_RUN_INFO"
    end if
    write(ifileid,"(a)") "      &END PRINT"
    write(ifileid,"(a)") "    &END SCF"
end if

!--- &PRINT of DFT level, FORCE_EVAL/DFT/PRINT
if (imolden==1.or.ioutSbas==1.or.ioutKSbas==1.or.ioutcube>0.or.iatomcharge>0.or.itask==5.or.imoment==1.or.ihyperfine==1.or.ioutorbene==1.or.iSCCS==1.or.iDFTplusU==1) then
    write(ifileid,"(a)") "    &PRINT"
    if (ioutSbas==1) then
        write(ifileid,"(a)") "      &S_CSR_WRITE #Exporting .csr file containing overlap matrix"
        write(ifileid,"(a)") "        REAL_SPACE T #Print the overlap matrix in real-space instead of k-space"
        write(ifileid,"(a)") "        UPPER_TRIANGULAR T #Print the matrix in upper triangular form"
        write(ifileid,"(a)") "        THRESHOLD 0 #Threshold on the absolute value of the elements to be printed out"
        write(ifileid,"(a)") "      &END S_CSR_WRITE"
    end if
    if (ioutKSbas==1) then
        write(ifileid,"(a)") "      &KS_CSR_WRITE #Exporting .csr file containing Kohn-Sham matrix"
        write(ifileid,"(a)") "        REAL_SPACE T #Print the Kohn-Sham matrix in real-space instead of k-space"
        write(ifileid,"(a)") "        UPPER_TRIANGULAR T #Print the matrix in upper triangular form"
        write(ifileid,"(a)") "        THRESHOLD 0 #Threshold on the absolute value of the elements to be printed out"
        write(ifileid,"(a)") "      &END KS_CSR_WRITE"
    end if
    if (ioutorbene==1) then
        write(ifileid,"(a)") "      &MO"
        write(ifileid,"(a)") "        ENERGIES T #Print orbital energies"
        write(ifileid,"(a)") "        OCCUPATION_NUMBERS T #Print orbital occupation numbers"
        write(ifileid,"(a)") "        COEFFICIENTS F #Print orbital coefficients"
        write(ifileid,"(a,2i6)") "        #MO_INDEX_RANGE",1,nint(max(naelec,nbelec))
        write(ifileid,"(a)") "        &EACH"
        write(ifileid,"(a)") "          QS_SCF 0"
        write(ifileid,"(a)") "        &END EACH"
        write(ifileid,"(a)") "      &END MO"
    end if
    if (imolden==1) then
        write(ifileid,"(a)") "      &MO_MOLDEN #Exporting .molden file containing wavefunction information"
        write(ifileid,"(a)") "        NDIGITS 9 #Output orbital coefficients if absolute value is larger than 1E-9"
        write(ifileid,"(a)") "      &END MO_MOLDEN"
    end if
    if (iDFTplusU==1) then
        write(ifileid,"(a)") "      #Uncomment following lines can print occupation for which +U is applied when PRINT_LEVEL is medium"
        write(ifileid,"(a)") "      #&PLUS_U"
        write(ifileid,"(a)") "      #  &EACH"
        write(ifileid,"(a)") "      #    QS_SCF 1"
        write(ifileid,"(a)") "      #  &END EACH"
        write(ifileid,"(a)") "      #&END PLUS_U"
    end if
    if (ioutcube>0) then
        if (ioutcube==1) then
            write(ifileid,"(a)") "      &E_DENSITY_CUBE"
            write(ifileid,"(a)") "        STRIDE 1 #Stride of exported cube file"
            write(ifileid,"(a)") "      &END E_DENSITY_CUBE"
        else if (ioutcube==2) then
            write(ifileid,"(a)") "      &ELF_CUBE"
            write(ifileid,"(a)") "        STRIDE 1 #Stride of exported cube file"
            write(ifileid,"(a)") "      &END ELF_CUBE"
        else if (ioutcube==3) then
            write(ifileid,"(a)") "      &V_XC_CUBE"
            write(ifileid,"(a)") "        STRIDE 1 #Stride of exported cube file"
            write(ifileid,"(a)") "      &END V_XC_CUBE"
        else if (ioutcube==4) then
            write(ifileid,"(a)") "      &V_HARTREE_CUBE"
            write(ifileid,"(a)") "        STRIDE 1 #Stride of exported cube file"
            write(ifileid,"(a)") "      &END V_HARTREE_CUBE"
        else if (ioutcube==5) then
            write(ifileid,"(a)") "      &EFIELD_CUBE"
            write(ifileid,"(a)") "        STRIDE 1 #Stride of exported cube file"
            write(ifileid,"(a)") "      &END EFIELD_CUBE"
        else if (ioutcube==6) then
            write(ifileid,"(a)") "      &MO_CUBES"
            write(ifileid,"(a)") "        STRIDE 1 #Stride of exported cube file"
            write(ifileid,"(a,i6)") "        NHOMO",NHOMO
            write(ifileid,"(a,i6)") "        NLUMO",NLUMO
            write(ifileid,"(a)") "      &END MO_CUBES"
        else if (ioutcube==7) then
            write(ifileid,"(a)") "      &E_DENSITY_CUBE"
            write(ifileid,"(a)") "        STRIDE 1 #Stride of exported cube file"
            write(ifileid,"(a)") "      &END E_DENSITY_CUBE"
            write(ifileid,"(a)") "      &V_HARTREE_CUBE"
            write(ifileid,"(a)") "        STRIDE 1 #Stride of exported cube file"
            write(ifileid,"(a)") "      &END V_HARTREE_CUBE"
        end if
    end if
    if (iatomcharge>0) then
        if (iatomcharge==1) then
            write(ifileid,"(a)") "      &MULLIKEN"
            write(ifileid,"(a)") "        PRINT_ALL F #If T, then printing full net AO and overlap population matrix"
            write(ifileid,"(a)") "      &END MULLIKEN"
        else if (iatomcharge==2) then
            write(ifileid,"(a)") "      &LOWDIN"
            write(ifileid,"(a)") "        PRINT_ALL F #If T, then printing full net AO and overlap population matrix"
            write(ifileid,"(a)") "      &END LOWDIN"
        else if (iatomcharge==3.or.iatomcharge==4) then
            write(ifileid,"(a)") "      &HIRSHFELD"
            write(ifileid,"(a)") "        SHAPE_FUNCTION DENSITY"
            if (iatomcharge==4) write(ifileid,"(a)") "        SELF_CONSISTENT T"
            write(ifileid,"(a)") "      &END HIRSHFELD"
        else if (iatomcharge==5) then
            write(ifileid,"(a)") "      &VORONOI"
            write(ifileid,"(a)") "        VORONOI_RADII Covalent" !Better than default of using vdW radii
            write(ifileid,"(a)") "      &END VORONOI"
        end if
    end if
    if (ihyperfine==1) then
        write(ifileid,"(a)") "      &HYPERFINE_COUPLING_TENSOR"
        write(ifileid,"(a)") "      &END HYPERFINE_COUPLING_TENSOR"
    end if
    if ((itask==5.and.ikpoint1==1.and.ikpoint2==1.and.ikpoint3==1).or.imoment==1) then !Note that vibration analysis cannot produce dipole moment when considering k-points
        write(ifileid,"(a)") "      &MOMENTS"
        if (ifPBC==0.or.PBCdir=="NONE") then
            write(ifileid,"(a)") "        PERIODIC F #Use Berry phase formula (T) or simple operator (F), the latter normally applies to isolated systems"
            write(ifileid,"(a)") "        MAGNETIC F #If also printing magnetic moments. Can only be used for isolated systems"
            write(ifileid,"(a)") "        REFERENCE COM #Reference point for calculating electric moment. COM=center of mass"
            write(ifileid,"(a)") "        MAX_MOMENT 1 #Maximum order of moments to be calculated. 1/2/3: Dipole/Quadrupole/Octapole moments. Up to 5"
        else
            write(ifileid,"(a)") "        PERIODIC T #Use Berry phase formula (T) or simple operator (F), the latter normally applies to isolated systems"
        end if
        if (itask==13) then
            write(ifileid,"(a)") "        FILENAME electric"
            write(ifileid,"(a)") "        COMMON_ITERATION_LEVELS 5"
        end if
        write(ifileid,"(a)") "      &END MOMENTS"
    end if
    if (iSCCS==1) then
        write(ifileid,"(a)") "      @IF 1 #Printing SCCS information in each SCF iteration" !When print level is medium, will also print SCCS iteration information
        write(ifileid,"(a)") "      &SCCS"
        write(ifileid,"(a)") "        &EACH"
        write(ifileid,"(a)") "          QS_SCF 1"
        write(ifileid,"(a)") "        &END EACH"
        write(ifileid,"(a)") "        &POLARISATION_CHARGE_DENSITY"
        write(ifileid,"(a)") "          &EACH"
        write(ifileid,"(a)") "            QS_SCF 0"
        write(ifileid,"(a)") "          &END EACH"
        write(ifileid,"(a)") "          STRIDE 2"
        write(ifileid,"(a)") "        &END POLARISATION_CHARGE_DENSITY"
        write(ifileid,"(a)") "        &DIELECTRIC_FUNCTION"
        write(ifileid,"(a)") "          &EACH"
        write(ifileid,"(a)") "            QS_SCF 0"
        write(ifileid,"(a)") "          &END EACH"
        write(ifileid,"(a)") "          STRIDE 2"
        write(ifileid,"(a)") "        &END DIELECTRIC_FUNCTION"
        write(ifileid,"(a)") "      &END SCCS"
        write(ifileid,"(a)") "      @ENDIF"
    end if
    write(ifileid,"(a)") "    &END PRINT"
end if
if (iSCCS==1) then
    write(ifileid,"(a)") "    &SCCS"
    write(ifileid,"(a)") "      ALPHA [mN/m] 57.2"
    write(ifileid,"(a)") "      BETA [GPa] -0.5"
    write(ifileid,"(a)") "      GAMMA [mN/m] 0.0"
    write(ifileid,"(a)") "      DIELECTRIC_CONSTANT 78.36"
    write(ifileid,"(a)") "      EPS_SCF 0.2 #SCCS is activated only if SCF iteration is converged to this threshold"
    write(ifileid,"(a)") "      EPS_SCCS 1E-6 #Requested accuracy for convergence of polarization charge density iteration"
    write(ifileid,"(a)") "      MAX_ITER 150 #Maximum number of polarization charge density iterations"
    write(ifileid,"(a)") "      DERIVATIVE_METHOD CD5 #Method for calculation of numerical derivatives. Can be FFT, CD3, CD5, CD7"
    write(ifileid,"(a)") "      &ANDREUSSI"
    write(ifileid,"(a)") "        RHO_MIN 0.0001841"
    write(ifileid,"(a)") "        RHO_MAX 0.0013604"
    write(ifileid,"(a)") "      &END ANDREUSSI"
    write(ifileid,"(a)") "    &END SCCS"
end if
if (idipcorr>0) then
    write(ifileid,"(a)") "    SURFACE_DIPOLE_CORRECTION T"
    if (idipcorr==1) write(ifileid,"(a)") "    SURF_DIP_DIR X"
    if (idipcorr==2) write(ifileid,"(a)") "    SURF_DIP_DIR Y"
    if (idipcorr==3) write(ifileid,"(a)") "    SURF_DIP_DIR Z"
end if
if (itask==13) then !Real-time propagation
    write(ifileid,"(a)") "    &REAL_TIME_PROPAGATION"
    write(ifileid,"(a)") "      INITIAL_WFN SCF_WFN #Initial wavefunction used for propagation is obtained by SCF. Can also be RESTART_WFN and RT_RESTART"
    write(ifileid,"(a)") "      EPS_ITER 1E-7 #Convergence criterion for the self consistent propagator loop. This is default value"
    write(ifileid,"(a)") "      MAX_ITER 50 #Maximal number of iterations for the self consistent propagator loop"
    write(ifileid,"(a)") "      APPLY_DELTA_PULSE #Applying a delta kick to the initial wavefunction"
    write(ifileid,"(a)") "      DELTA_PULSE_DIRECTION 0 0 1 #Direction of the applied electric field"
    write(ifileid,"(a)") "      &PRINT"
    write(ifileid,"(a)") "        &RESTART"
    write(ifileid,"(a)") "          BACKUP_COPIES 0 #Never generate backed up .rtpwfn files"
    write(ifileid,"(a)") "        &END RESTART"
    write(ifileid,"(a)") "      &END PRINT"
    write(ifileid,"(a)") "    &END REAL_TIME_PROPAGATION"
end if
if (itask==15) then !XAS
    write(ifileid,"(a)") "    &XAS_TDP"
    write(ifileid,"(a)") "      GRID "//kindnameXAS//" 150 200"
    write(ifileid,"(a)") "      ENERGY_RANGE 30 #The energy range in eV for which excitations are considered"
    write(ifileid,"(a)") "      &DONOR_STATES"
    write(ifileid,"(a)") "        DEFINE_EXCITED BY_INDEX"
    call outCP2K_LIST(ifileid,XASatm,nXASatm,"        ","ATOM_LIST")
    write(ifileid,"(a)",advance="no") "        STATE_TYPES"
    do iatm=1,nXASatm
        write(ifileid,"(a)",advance="no") " 1S"
    end do
    write(ifileid,*) "#Types of orbitals that are excited"
    if (iGW2X==1) then !GW2X requires LOCALIZE
        write(ifileid,"(a)") "        LOCALIZE #Perform orbital localization before searching potential donor states"
        write(ifileid,"(a,i6)") "        N_SEARCH 1"
    else
        write(ifileid,"(a)") "        N_SEARCH -1 #-1 means searching donor orbitals from all occupied orbitals from low to high energy"
        write(ifileid,"(a)") "        #LOCALIZE #Perform orbital localization before searching potential donor states"
    end if
    write(ifileid,"(a)") "      &END DONOR_STATES"
    write(ifileid,"(a)") "      &KERNEL"
    write(ifileid,"(a)") "        RI_REGION 2.0"
    write(ifileid,"(a)") "        &XC_FUNCTIONAL PBE"
    write(ifileid,"(a)") "          &PBE"
    write(ifileid,"(a,f10.5)") "            SCALE_X",(100-PBEh_HFX)/100
    write(ifileid,"(a)") "          &END PBE"
    write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
    write(ifileid,"(a)") "        &EXACT_EXCHANGE"
    write(ifileid,"(a,f10.5,a)") "          FRACTION ",PBEh_HFX/100," #HF composition"
    if (PBCdir/="NONE") then
        write(ifileid,"(a)") "          POTENTIAL_TYPE TRUNCATED"
        if (trunc_rad>6) then !If half of shortest box length is larger than 6 Angstrom, simply use 6, this is usually adequate
            write(ifileid,"(a,f8.4)") "          CUTOFF_RADIUS 6.0 #Cutoff radius (Angstrom) for truncated 1/r Coulomb operator"
        else
            write(ifileid,"(a,f8.4,a)") "          CUTOFF_RADIUS",trunc_rad," #Cutoff radius (Angstrom) for truncated 1/r Coulomb operator"
        end if
    end if
    write(ifileid,"(a)") "        &END EXACT_EXCHANGE"
    write(ifileid,"(a)") "      &END KERNEL"
    if (iXAS_SOC==1) then
        if (multispin>1.or.any(kindmag(1:nkind)/=0)) then !Didn't test
            write(ifileid,"(a)") "      EXCITATIONS OS_SPIN_CONS #Spin-conserving excitations on top of open-shell ground state"
            write(ifileid,"(a)") "      EXCITATIONS OS_SPIN_FLIP #Spin-flip excitation on top of open-shell ground state"
        else !Closed-shell
            write(ifileid,"(a)") "      EXCITATIONS RCS_SINGLET #Calculate singlet excitation on top of restricted closed-shell ground state"
            write(ifileid,"(a)") "      EXCITATIONS RCS_TRIPLET #Calculate triplet excitation on top of restricted closed-shell ground state"
        end if
        write(ifileid,"(a)") "      SOC #Consider spin-orbit coupling"
    end if
    if (iGW2X==1) then
        write(ifileid,"(a)") "      &GW2X #Enable GW2X correction"
        write(ifileid,"(a)") "      &END GW2X"
    end if
    write(ifileid,"(a)") "    &END XAS_TDP"
end if
write(ifileid,"(a)") "  &END DFT"

end if !END distinguishing FIST and DFT


if (itask==2.or.itask==4.or.inoSCFinfo==1.or.method=="FIST") then !&FORCE_EVAL / &PRINT
    write(ifileid,"(a)") "  &PRINT"
    if (itask==2) then
        write(ifileid,"(a)") "    &FORCES ON #Print atomic forces"
        write(ifileid,"(a)") "    &END FORCES"
    end if
    if (itask==4) then !CELL_OPT
        write(ifileid,"(a)") "    &STRESS_TENSOR ON #Print stress tensor"
        write(ifileid,"(a)") "    &END STRESS_TENSOR"
    end if
    if (inoSCFinfo==1.or.method=="FIST") then
        write(ifileid,"(a)") "    &PROGRAM_RUN_INFO"
        write(ifileid,"(a)") "      &EACH"
        write(ifileid,"(a)") "        MD 0 #Frequency of printing evaluated energies during MD. 0 means never"
        write(ifileid,"(a)") "      &END EACH"
        write(ifileid,"(a)") "    &END PROGRAM_RUN_INFO"
    end if
    write(ifileid,"(a)") "  &END PRINT"
end if

if (itask==4.or.ibarostat>0) write(ifileid,"(a)") "  STRESS_TENSOR ANALYTICAL #Compute full stress tensor analytically" !Compute for CELL_OPT and MD with barostat
if ((itask==5.and.iraman==1).or.itask==9.or.itask==10.or.iTDDFT==1) then !Raman, NMR, polar, TDDFT
    write(ifileid,"(a)") "  &PROPERTIES"
    if ((itask==5.and.iraman==1).or.itask==9.or.itask==10) then !Raman, NMR, polar
        write(ifileid,"(a)") "    &LINRES #Activate linear response calculation"
        if (itask==5) then !Raman
            write(ifileid,"(a)") "      PRECONDITIONER FULL_ALL #Preconditioner for conjugate gradient minimization"
            write(ifileid,"(a)") "      EPS 1E-6 #Target accuracy for the convergence of the conjugate gradient" !Identical to default, already sufficient for Raman spectrum
        else if (itask==9) then !NMR
            write(ifileid,"(a)") "      PRECONDITIONER FULL_KINETIC #Preconditioner for conjugate gradient minimization"
            write(ifileid,"(a)") "      EPS 1E-8 #Target accuracy for the convergence of the conjugate gradient" !Tigher than default 1E-6
            write(ifileid,"(a)") "      MAX_ITER 400 #Maximum number of conjugate gradient iterations to be performed" !Sometimes converges very slowly
        else if (itask==10) then !Polar
            write(ifileid,"(a)") "      PRECONDITIONER FULL_ALL #Preconditioner for conjugate gradient minimization"
            write(ifileid,"(a)") "      EPS 1E-7 #Target accuracy for the convergence of the conjugate gradient" !Tigher than default 1E-6
            write(ifileid,"(a)") "      MAX_ITER 100 #Maximum number of conjugate gradient iterations to be performed" !Larger than the default 50
        end if
        if (itask==9) then !NMR
            write(ifileid,"(a)") "      &CURRENT"
            write(ifileid,"(a)") "        GAUGE R_AND_STEP_FUNCTION #The gauge used to compute induced current: ATOM, R, R_AND_STEP_FUNCTION"
            write(ifileid,"(a)") "        ORBITAL_CENTER ATOM #Orbital center: WANNIER, ATOM, BOX, COMMON"
            write(ifileid,"(a)") "        CHI_PBC T #Calculate the succeptibility correction to the shift with PBC"
            write(ifileid,"(a)") "      &END CURRENT"
            write(ifileid,"(a)") "      &LOCALIZE"
            !I found CRAZY doesn't properly work, so do not explicitly mention, just use default JACOBI
            !write(ifileid,"(a)") "	      METHOD JACOBI #Localization optimization method. JACOBI=2x2 orbital rotations. CRAZY is less robust but usually much faster"
            write(ifileid,"(a)") "	      MAX_ITER 20000 #Maximum number of iterations used for localization methods"
            write(ifileid,"(a)") "	      EPS_LOCALIZATION 1E-5 #Convergence criterion of orbital localization"
            !I found BOYS and PIPEK do not work, so do not explicitly mention them, just use default BERRY
            !write(ifileid,"(a)") "	      OPERATOR BERRY #The quantity to be minimized in localization. Can also be BOYS and PIPEK"
            write(ifileid,"(a)") "      &END LOCALIZE"
            write(ifileid,"(a)") "      &NMR"
            write(ifileid,"(a)") "      #  NICS T #Calculate NICS"
            write(ifileid,"(a)") "      #  NICS_FILE_NAME filepath #Path of the file containing NICS points coordinates"
            write(ifileid,"(a)") "      &END NMR"
        else if (itask==5.or.itask==10) then !Raman or polar
            write(ifileid,"(a)") "      &POLAR"
            if (itask==5) write(ifileid,"(a)") "        DO_RAMAN T #Compute the electric-dipole--electric-dipole polarizability" !This is default
            if (ifPBC==0.or.PBCdir=="NONE") then
                write(ifileid,"(a)") "        PERIODIC_DIPOLE_OPERATOR F #Type of dipole operator: Berry phase(T) or Local(F)"
            else
                write(ifileid,"(a)") "        PERIODIC_DIPOLE_OPERATOR T #Type of dipole operator: Berry phase(T) or Local(F)"
            end if
            write(ifileid,"(a)") "      &END POLAR"
        end if
        write(ifileid,"(a)") "    &END LINRES"
    end if
    if (iTDDFT==1) then !TDDFT
        write(ifileid,"(a)") "    &TDDFPT #TDDFT calculation with Tamm-Dancoff approximation"
        write(ifileid,"(a,i5,a)") "      NSTATES",nstates_TD," #Number of excited states to solve"
        if (isTDA==1) then
            write(ifileid,"(a)") "      KERNEL STDA #Using sTDA approximation"
            write(ifileid,"(a)") "      &STDA"
            if (ifPBC==0.or.PBCdir=="NONE") then
                write(ifileid,"(a)") "        DO_EWALD F #Use Ewald type method for periodic Coulomb interaction"
            else
                write(ifileid,"(a)") "        DO_EWALD T #Use Ewald type method for periodic Coulomb interaction"
            end if
            if (index(method,"PBE0")/=0) then
                write(ifileid,"(a)") "        FRACTION 0.25 #Fraction of TB Hartree-Fock exchange to use in the kernel"
            else if (index(method,"BHandHLYP")/=0) then
                write(ifileid,"(a)") "        FRACTION 0.5 #Fraction of TB Hartree-Fock exchange to use in the kernel"
            else if (index(method,"PBEh")/=0) then
                write(ifileid,"(a,f10.5,a)") "        FRACTION",PBEh_HFX/100," #Fraction of TB Hartree-Fock exchange to use in the kernel"
            else if (index(method,"M06-2X")/=0) then
                write(ifileid,"(a)") "        FRACTION 0.54 #Fraction of TB Hartree-Fock exchange to use in the kernel"
            else
                write(ifileid,"(a)") "        FRACTION 0.2 #Fraction of TB Hartree-Fock exchange to use in the kernel"
            end if
            write(ifileid,"(a)") "      &END STDA"
        end if
        if (iTDtriplet==1.or.iSOCTDDFT==1) then
            write(ifileid,"(a)") "      RKS_TRIPLETS T #If calculating triplet rather than singlet excited states"
        else
            write(ifileid,"(a,i5)") "      RKS_TRIPLETS F #If calculating triplet rather than singlet excited states"
        end if
        if (itask==2.or.itask==3.or.itask==4.or.itask==5.or.itask==6.or.itask==7.or.itask==8) then !Involve force
            write(ifileid,"(a)") "      CONVERGENCE [eV] 1E-6 #Convergence criterion of all excitation energies"
        else
            write(ifileid,"(a)") "      CONVERGENCE [eV] 1E-4 #Convergence criterion of all excitation energies"
        end if
        write(ifileid,"(a)") "      MIN_AMPLITUDE 0.01 #The smallest excitation amplitude to print"
        write(ifileid,"(a)") "#     RESTART T #If restarting TDDFT calculation. If true, WFN_RESTART_FILE_NAME should be set to previous .tdwfn file"
        write(ifileid,"(a)") "#     WFN_RESTART_FILE_NAME "//trim(c200tmp)//"-RESTART.tdwfn"
        !if (ifPBC==3) then !The default VELOCITY is also correct for both periodic and isolated systems
        !    write(ifileid,"(a)") "      &DIPOLE_MOMENTS"
        !    write(ifileid,"(a)") "        DIPOLE_FORM BERRY"
        !    write(ifileid,"(a)") "      &END DIPOLE_MOMENTS"
        !end if
        if (index(method,"_ADMM")/=0.and.(itask==2.or.itask==3.or.itask==4.or.itask==5.or.itask==6.or.itask==7.or.itask==8)) then !Need force
            write(ifileid,"(a)") "      ADMM_KERNEL_CORRECTION_SYMMETRIC T"
        end if
        if (PBCdir=="NONE") then
            write(ifileid,"(a)") "      &DIPOLE_MOMENTS"
            write(ifileid,"(a)") "        DIPOLE_FORM LENGTH"
            write(ifileid,"(a)") "      &END DIPOLE_MOMENTS"
        end if
        if (iSOCTDDFT==1) then
            write(ifileid,"(a)") "      &SOC #Consider SOC effect in TDDFT calculation"
            write(ifileid,"(a)") "      &END SOC"
        end if
        write(ifileid,"(a)") "      &PRINT"
        write(ifileid,"(a)") "        #&DETAILED_ENERGY #Print excitation energies at every Davidson iteration"
        write(ifileid,"(a)") "        #&END DETAILED_ENERGY"
        if (iSOCTDDFT==1) then
            write(ifileid,"(a)") "        &SOC_PRINT"
            write(ifileid,"(a)") "          SPLITTING #Print SOC-splitting"
            write(ifileid,"(a)") "          SOME #Print SOC-matrix elements"
            write(ifileid,"(a)") "        &END SOC_PRINT"
        end if
        if (iNTO==1) then
            write(ifileid,"(a)") "        &NTO_ANALYSIS #Do NTO analysis for all excited states"
            !write(ifileid,"(a)") "          FILENAME NTO" !Seems not to affect name of actually exported .molden file
            write(ifileid,"(a)") "        &END NTO_ANALYSIS"
            write(ifileid,"(a)") "        &MOS_MOLDEN #Output .molden file containing NTO of the ""NSTATES""th state"
            write(ifileid,"(a)") "          NDIGITS 9"
            write(ifileid,"(a)") "          FILENAME NTO #Filename of NTO .molden file"
            write(ifileid,"(a)") "        &END MOS_MOLDEN"
        end if
        write(ifileid,"(a)") "      &END PRINT"
        !Current version, at least 2022.1, does not force to specify &TDDFPT / &XC
        !if (index(method,"_ADMM")==0) then !When ADMM is used for TDDFT, &XC should not appear, otherwise error shows: "ADMM is not implemented for a TDDFT kernel XC-functional which is different from the one used for the ground-state calculation"
        !    write(ifileid,"(a)") "      &XC"
        !    write(ifileid,"(a)") "        &XC_GRID"
        !    write(ifileid,"(a)") "          XC_DERIV SPLINE2_SMOOTH #The method used to compute the derivatives"
        !    write(ifileid,"(a)") "        &END XC_GRID"
        !    !XC functional for TDDFT
        !    if (index(method,"LIBXC")/=0) then
        !        write(ifileid,"(a)") "        &XC_FUNCTIONAL"
        !        if (method=="B97M-rV_LIBXC") then !Non-separable XC
        !            write(ifileid,"(a)") "          &MGGA_XC_B97M_V"
        !            write(ifileid,"(a)") "          &END MGGA_XC_B97M_V"
        !        else !X-C separable
        !            if (method=="MN15L_LIBXC") then
        !                write(ifileid,"(a)") "          &MGGA_X_MN15_L"
        !                write(ifileid,"(a)") "          &END MGGA_X_MN15_L"
        !                write(ifileid,"(a)") "          &MGGA_C_MN15_L"
        !                write(ifileid,"(a)") "          &END MGGA_C_MN15_L"
        !            else if (method=="SCAN_LIBXC") then
        !                write(ifileid,"(a)") "          &MGGA_X_SCAN"
        !                write(ifileid,"(a)") "          &END MGGA_X_SCAN"
        !                write(ifileid,"(a)") "          &MGGA_C_SCAN"
        !                write(ifileid,"(a)") "          &END MGGA_C_SCAN"
        !            else if (method=="r2SCAN_LIBXC") then
        !                write(ifileid,"(a)") "          &MGGA_X_R2SCAN"
        !                write(ifileid,"(a)") "          &END MGGA_X_R2SCAN"
        !                write(ifileid,"(a)") "          &MGGA_C_R2SCAN"
        !                write(ifileid,"(a)") "          &END MGGA_C_R2SCAN"
        !            else if (method=="RPBE_LIBXC") then
        !                write(ifileid,"(a)") "          &GGA_X_RPBE"
        !                write(ifileid,"(a)") "          &END GGA_X_RPBE"
        !                write(ifileid,"(a)") "          &GGA_C_PBE"
        !                write(ifileid,"(a)") "          &END GGA_C_PBE"
        !            else if (method=="TPSS_LIBXC") then
        !                write(ifileid,"(a)") "          &MGGA_X_TPSS"
        !                write(ifileid,"(a)") "          &END MGGA_X_TPSS"
        !                write(ifileid,"(a)") "          &MGGA_C_TPSS"
        !                write(ifileid,"(a)") "          &END MGGA_C_TPSS"
        !            else if (method=="revTPSS_LIBXC") then
        !                write(ifileid,"(a)") "          &MGGA_X_REVTPSS"
        !                write(ifileid,"(a)") "          &END MGGA_X_REVTPSS"
        !                write(ifileid,"(a)") "          &MGGA_C_REVTPSS"
        !                write(ifileid,"(a)") "          &END MGGA_C_REVTPSS"
        !            else if (method=="HLE17_LIBXC") then
        !                write(ifileid,"(a)") "          &MGGA_XC_HLE17"
        !                write(ifileid,"(a)") "          &END MGGA_XC_HLE17"
        !            end if
        !        end if
        !        write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
        !    else if (index(method,"PBE0")/=0) then
        !        write(ifileid,"(a)") "        &XC_FUNCTIONAL PBE0"
        !        write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
        !    else if (index(method,"B3LYP")/=0) then
        !        write(ifileid,"(a)") "        &XC_FUNCTIONAL B3LYP"
        !        write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
        !    else if (index(method,"BHandHLYP")/=0) then
        !        write(ifileid,"(a)") "        &XC_FUNCTIONAL"
        !        write(ifileid,"(a)") "          &HYB_GGA_XC_BHANDHLYP"
        !        write(ifileid,"(a)") "          &END HYB_GGA_XC_BHANDHLYP"
        !        write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
        !    else if (index(method,"M06-2X")/=0) then
        !        write(ifileid,"(a)") "        &XC_FUNCTIONAL"
        !        write(ifileid,"(a)") "          &HYB_MGGA_X_M06_2X"
        !        write(ifileid,"(a)") "          &END HYB_MGGA_X_M06_2X"
        !        write(ifileid,"(a)") "          &MGGA_C_M06_2X"
        !        write(ifileid,"(a)") "          &END MGGA_C_M06_2X"
        !        write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
        !    else if (method=="revPBE".or.method=="PBEsol") then
        !        write(ifileid,"(a)") "        &XC_FUNCTIONAL PBE"
        !        write(ifileid,"(a)") "          &PBE"
        !        if (method=="revPBE") write(ifileid,"(a)") "          PARAMETRIZATION REVPBE"
        !        if (method=="PBEsol") write(ifileid,"(a)") "          PARAMETRIZATION PBESOL"
        !        write(ifileid,"(a)") "          &END PBE"
        !        write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
        !    else if (index(method,"HSE")/=0) then
        !        write(ifileid,"(a)") "        &XC_FUNCTIONAL"
        !        write(ifileid,"(a)") "          &XWPBE"
        !        write(ifileid,"(a)") "            SCALE_X -0.25"
        !        write(ifileid,"(a)") "            SCALE_X0 1.0"
        !        write(ifileid,"(a)") "            OMEGA 0.11"
        !        write(ifileid,"(a)") "          &END XWPBE"
        !        write(ifileid,"(a)") "          &PBE"
        !        write(ifileid,"(a)") "            SCALE_X 0.0"
        !        write(ifileid,"(a)") "            SCALE_C 1.0"
        !        write(ifileid,"(a)") "          &END PBE"
        !        write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
        !    else !Common native GGA functionals
        !        write(ifileid,"(a)") "        &XC_FUNCTIONAL "//trim(method)
        !        write(ifileid,"(a)") "        &END XC_FUNCTIONAL"
        !    end if
        !    write(ifileid,"(a)") "      &END XC"
        !end if
        write(ifileid,"(a)") "    &END TDDFPT"
    end if
    write(ifileid,"(a)") "  &END PROPERTIES"
end if
if (itask==11) then !BSSE setting
    write(ifileid,"(/,a)") "  &BSSE"
    write(ifileid,"(a)") "    &FRAGMENT"
    call outCP2K_LIST(ifileid,frag1,nfrag1,"      ","LIST")
    write(ifileid,"(a)") "    &END FRAGMENT"
    write(ifileid,"(a)") "    &FRAGMENT"
    call outCP2K_LIST(ifileid,frag2,nfrag2,"      ","LIST")
    write(ifileid,"(a)") "    &END FRAGMENT"
    write(ifileid,"(a)") "    &CONFIGURATION # real(A)+real(B)"
    write(ifileid,"(a)") "      GLB_CONF 1 1"
    write(ifileid,"(a)") "      SUB_CONF 1 1"
    write(ifileid,"(a,i3)") "      CHARGE",frag1chg+frag2chg
    write(ifileid,"(a,i3)") "      MULTIPLICITY ",totalmulti
    write(ifileid,"(a)") "    &END CONFIGURATION"
    write(ifileid,"(a)") "    &CONFIGURATION # real(A)"
    write(ifileid,"(a)") "      GLB_CONF 1 0"
    write(ifileid,"(a)") "      SUB_CONF 1 0"
    write(ifileid,"(a,i3)") "      CHARGE",frag1chg
    write(ifileid,"(a,i3)") "      MULTIPLICITY ",frag1multi
    write(ifileid,"(a)") "    &END CONFIGURATION"
    write(ifileid,"(a)") "    &CONFIGURATION # real(B)"
    write(ifileid,"(a)") "      GLB_CONF 0 1"
    write(ifileid,"(a)") "      SUB_CONF 0 1"
    write(ifileid,"(a,i3)") "      CHARGE",frag2chg
    write(ifileid,"(a,i3)") "      MULTIPLICITY ",frag2multi
    write(ifileid,"(a)") "    &END CONFIGURATION"
    write(ifileid,"(a)") "    &CONFIGURATION # real(A)+ghost(B)"
    write(ifileid,"(a)") "      GLB_CONF 1 1"
    write(ifileid,"(a)") "      SUB_CONF 1 0"
    write(ifileid,"(a,i3)") "      CHARGE",frag1chg
    write(ifileid,"(a,i3)") "      MULTIPLICITY ",frag1multi
    write(ifileid,"(a)") "    &END CONFIGURATION"
    write(ifileid,"(a)") "    &CONFIGURATION # ghost(A)+real(B)"
    write(ifileid,"(a)") "      GLB_CONF 1 1"
    write(ifileid,"(a)") "      SUB_CONF 0 1"
    write(ifileid,"(a,i3)") "      CHARGE",frag2chg
    write(ifileid,"(a,i3)") "      MULTIPLICITY ",frag2multi
    write(ifileid,"(a)") "    &END CONFIGURATION"
    write(ifileid,"(a)") "  &END BSSE"
end if
if (iatomcharge==6.or.iatomcharge==7) then !Atomic charges
    write(ifileid,"(/,a)") "  &PROPERTIES"
    write(ifileid,"(a)") "    &RESP"
    if (iatomcharge==7) write(ifileid,"(a)") "      USE_REPEAT_METHOD T"
    write(ifileid,"(a)") "      &SPHERE_SAMPLING"
    write(ifileid,"(a)") "        AUTO_VDW_RADII_TABLE CAMBRIDGE #vdW radii type. This is default. Can also be UFF"
    write(ifileid,"(a)") "        AUTO_RMIN_SCALE 1.0 #Scaled factor of vdW radii determining the inner boundary of sampling"
    write(ifileid,"(a)") "        AUTO_RMAX_SCALE 2.0 #Scaled factor of vdW radii determining the outer boundary of sampling"
    write(ifileid,"(a)") "      &END SPHERE_SAMPLING"
    if (ifPBC>0) then
        write(ifileid,"(a)") "      #Uncomment following lines can use slab sampling of fitting points"
        write(ifileid,"(a)") "      #&SLAB_SAMPLING #The fitting points will sampled above a slab"
        write(ifileid,"(a)") "      #  RANGE 2.0 4.0"
        write(ifileid,"(a)") "      #  LENGTH 3.0"
        write(ifileid,"(a)") "      #  ATOM_LIST 1..32 #List of considered atoms"
        write(ifileid,"(a)") "      #  SURF_DIRECTION Z #What above the surface means. Can also be e.g. -Z, X, Z..."
        write(ifileid,"(a)") "      #&END SLAB_SAMPLING"
    end if
    write(ifileid,"(a)") "      &PRINT"
    write(ifileid,"(a)") "        &COORD_FIT_POINTS"
    write(ifileid,"(a)") "        &END COORD_FIT_POINTS"
    write(ifileid,"(a)") "        &RESP_CHARGES_TO_FILE"
    write(ifileid,"(a)") "        &END RESP_CHARGES_TO_FILE"
    write(ifileid,"(a)") "      &END PRINT"
    write(ifileid,"(a)") "    &END RESP"
    write(ifileid,"(a)") "  &END PROPERTIES"
end if
write(ifileid,"(a)") "&END FORCE_EVAL"

!--- &MOTION
if (itask==3.or.itask==4.or.itask==5.or.itask==6.or.itask==7.or.itask==8.or.itask==13.or.itask==14) then
    write(ifileid,"(/,a)") "&MOTION"
    if (itask==3.or.itask==7) then !Optimizing atoms for minimum or TS
        write(ifileid,"(a)") "  &GEO_OPT"
        if (itask==3) then
            write(ifileid,"(a)") "    TYPE MINIMIZATION #Search for minimum"
            write(ifileid,"(a)") "    KEEP_SPACE_GROUP F #If T, then space group will be detected and preserved"
            if (ioptmethod==1) then
                write(ifileid,"(a)") "    OPTIMIZER BFGS #Can also be CG (more robust for difficult cases) or LBFGS"
                write(ifileid,"(a)") "    &BFGS"
                write(ifileid,"(a)") "      TRUST_RADIUS 0.2 #Trust radius (maximum stepsize) in Angstrom"
                write(ifileid,"(a)") "#     RESTART_HESSIAN T #If read initial Hessian, uncomment this line and specify the file in the next line"
                write(ifileid,"(a)") "#     RESTART_FILE_NAME to_be_specified"
                write(ifileid,"(a)") "    &END BFGS"
            else if (ioptmethod==2) then
                write(ifileid,"(a)") "    OPTIMIZER LBFGS #Can also be CG (more robust for difficult cases) or BFGS"
                write(ifileid,"(a)") "    &LBFGS"
                write(ifileid,"(a)") "      TRUST_RADIUS 0.2 #Trust radius (maximum stepsize) in Angstrom"
                write(ifileid,"(a)") "      MAX_H_RANK 5 #Larger values (e.g. 30) will accelerate convergence behaviour at the cost of a larger memory consumption"
                write(ifileid,"(a)") "    &END LBFGS"
            else if (ioptmethod==3) then
                write(ifileid,"(a)") "    OPTIMIZER CG #Can also be BFGS or LBFGS"
                write(ifileid,"(a)") "    &CG"
                write(ifileid,"(a)") "      &LINE_SEARCH"
                write(ifileid,"(a)") "        TYPE 2PNT #Two-point extrapolation, cheap while acceptable. Can also be FIT, GOLD"
                write(ifileid,"(a)") "      &END LINE_SEARCH"
                write(ifileid,"(a)") "    &END CG"
            end if
        else if (itask==7) then
            write(ifileid,"(a)") "    TYPE TRANSITION_STATE #Optimizing TS using dimer algorithm"
            write(ifileid,"(a)") "    OPTIMIZER CG" !CG is the only choice for dimer
            write(ifileid,"(a)") "    &CG"
            write(ifileid,"(a)") "      &LINE_SEARCH"
            write(ifileid,"(a)") "        TYPE 2PNT"
            write(ifileid,"(a)") "      &END LINE_SEARCH"
            write(ifileid,"(a)") "    &END CG"
            write(ifileid,"(a)") "    &TRANSITION_STATE"
            write(ifileid,"(a)") "      &DIMER"
            !write(ifileid,"(a)") "        DR 0.01 #Default. DR parameter"
            write(ifileid,"(a)") "        ANGLE_TOLERANCE [deg] 5.0 #Tolerance angle for line search performed to optimize dimer orientation"
            write(ifileid,"(a)") "        &ROT_OPT #How to optimize dimer rotation"
            write(ifileid,"(a)") "          OPTIMIZER CG"
            write(ifileid,"(a)") "          MAX_ITER 200 #Maximum number of optimization steps"
            write(ifileid,"(a)") "          MAX_DR 3E-3 #Maximum geometry change"
            write(ifileid,"(a)") "          RMS_DR 1.5E-3 #RMS geometry change"
            if (itightopt==1) then
                write(ifileid,"(a)") "          MAX_FORCE 1E-4 #Maximum force"
            else
                write(ifileid,"(a)") "          MAX_FORCE 4.5E-4 #Maximum force"
            end if
            write(ifileid,"(a)") "          RMS_FORCE 3E-4 #RMS force"
            write(ifileid,"(a)") "          &CG"
            write(ifileid,"(a)") "            &LINE_SEARCH"
            write(ifileid,"(a)") "              TYPE 2PNT"
            write(ifileid,"(a)") "            &END LINE_SEARCH"
            write(ifileid,"(a)") "          &END CG"
            write(ifileid,"(a)") "        &END ROT_OPT"
            write(ifileid,"(a)") "      &END DIMER"
            write(ifileid,"(a)") "    &END TRANSITION_STATE"
        end if
        if (method=="FIST") then
            write(ifileid,"(a)") "    MAX_ITER 3000 #Maximum number of geometry optimization"
        else
            write(ifileid,"(a)") "    MAX_ITER 500 #Maximum number of geometry optimization"
        end if
        write(ifileid,"(a)") "    MAX_DR 3E-3 #Maximum geometry change"
        write(ifileid,"(a)") "    RMS_DR 1.5E-3 #RMS geometry change"
        if (itightopt==1) then
            write(ifileid,"(a)") "    MAX_FORCE 1E-4 #Maximum force"
        else
            write(ifileid,"(a)") "    MAX_FORCE 4.5E-4 #Maximum force"
        end if
        write(ifileid,"(a)") "    RMS_FORCE 3E-4 #RMS force"
        write(ifileid,"(a)") "  &END GEO_OPT"
    else if (itask==4) then
        write(ifileid,"(a)") "  &CELL_OPT"
        write(ifileid,"(a)") "    MAX_ITER 400 #Maximum number of geometry optimization"
        if (iprestype==1) then
            write(ifileid,"(a,1PE13.5,a)") "    EXTERNAL_PRESSURE",Piso," #External pressure for cell optimization (bar)"
        else if (iprestype==2) then
            write(ifileid,"(a,9(1PE13.5),a)") "    EXTERNAL_PRESSURE",Ptens(1,1:3),Ptens(2,1:3),Ptens(3,1:3)," #External pressure for cell optimization (bar)"
        end if
        write(ifileid,"(a)") "    CONSTRAINT "//trim(cellfix)//" #Constraint of cell length, can be: NONE, X, Y, Z, XY, XZ, YZ"
        write(ifileid,"(a)") "    KEEP_ANGLES F #If T, then cell angles will be kepted"
        write(ifileid,"(a)") "    KEEP_SYMMETRY F #If T, then cell symmetry specified by &CELL / SYMMETRY will be kepted"
        write(ifileid,"(a)") "    KEEP_SPACE_GROUP F #If T, then space group will be detected and preserved"
        write(ifileid,"(a)") "    TYPE DIRECT_CELL_OPT #Geometry and cell are optimized at the same time. Can also be GEO_OPT, MD"
        write(ifileid,"(a)") "    #The following thresholds of optimization convergence are the default ones"
        write(ifileid,"(a)") "    MAX_DR 3E-3 #Maximum geometry change"
        write(ifileid,"(a)") "    RMS_DR 1.5E-3 #RMS geometry change"
        if (itightopt==1) then
            write(ifileid,"(a)") "    MAX_FORCE 1E-4 #Maximum force"
        else
            write(ifileid,"(a)") "    MAX_FORCE 4.5E-4 #Maximum force"
        end if
        write(ifileid,"(a)") "    RMS_FORCE 3E-4 #RMS force"
        if (itightopt==0) then
            write(ifileid,"(a)") "    PRESSURE_TOLERANCE 100 #Pressure tolerance (w.r.t EXTERNAL_PRESSURE)"
        else if (itightopt==1) then
            write(ifileid,"(a)") "    PRESSURE_TOLERANCE 50 #Pressure tolerance (w.r.t EXTERNAL_PRESSURE)"
        end if
        if (ioptmethod==1) then
            write(ifileid,"(a)") "    OPTIMIZER BFGS #Can also be CG (more robust for difficult cases) or LBFGS"
            write(ifileid,"(a)") "    &BFGS"
            write(ifileid,"(a)") "      TRUST_RADIUS 0.2 #Trust radius (maximum stepsize) in Angstrom"
            write(ifileid,"(a)") "#     RESTART_HESSIAN T #If read initial Hessian, uncomment this line and specify the file in the next line"
            write(ifileid,"(a)") "#     RESTART_FILE_NAME to_be_specified"
            write(ifileid,"(a)") "    &END BFGS"
        else if (ioptmethod==2) then
            write(ifileid,"(a)") "    OPTIMIZER LBFGS #Can also be CG (more robust for difficult cases) or BFGS"
        else if (ioptmethod==3) then
            write(ifileid,"(a)") "    OPTIMIZER CG #Can also be BFGS or LBFGS"
            write(ifileid,"(a)") "    &CG"
            write(ifileid,"(a)") "      &LINE_SEARCH"
            write(ifileid,"(a)") "        TYPE 2PNT #Two-point extrapolation, cheap while acceptable. Can also be FIT, GOLD"
            write(ifileid,"(a)") "      &END LINE_SEARCH"
            write(ifileid,"(a)") "    &END CG"
        end if
        write(ifileid,"(a)") "  &END CELL_OPT"
    else if (itask==8) then !NEB
        write(ifileid,"(a)") "  &BAND"
        write(ifileid,"(a)") "    K_SPRING 0.05 #Spring constant"
        write(ifileid,"(a)") "    BAND_TYPE CI-NEB #The type of BAND calculation"
        write(ifileid,"(a)") "    NUMBER_OF_REPLICA 10 #The number of replica to use in the BAND"
        write(ifileid,"(a)") "    NPROC_REP 1 #The number of processors to be used per replica"
        if (ifPBC==0) then
            write(ifileid,"(a)") "    ALIGN_FRAMES T #Enables the alignment of the frames at the beginning of the BAND calculation"
            write(ifileid,"(a)") "    ROTATE_FRAMES T #Compute at each BAND step the RMSD and rotate the frames in order to minimize it"
        else
            write(ifileid,"(a)") "    ALIGN_FRAMES F #Enables the alignment of the frames at the beginning of the BAND calculation"
            write(ifileid,"(a)") "    ROTATE_FRAMES F #Compute at each BAND step the RMSD and rotate the frames in order to minimize it"
        end if
        write(ifileid,"(a)") "    &CI_NEB"
        write(ifileid,"(a)") "      NSTEPS_IT  5 #The number of steps of IT-NEB to perform before switching on CI-NEB"
        write(ifileid,"(a)") "    &END"
        write(ifileid,"(a)") "    &OPTIMIZE_BAND"
        write(ifileid,"(a)") "      OPTIMIZE_END_POINTS F #If also optimizing the end points of the band"
        write(ifileid,"(a)") "      &DIIS #Parameters of optimizing band via DIIS method"
        write(ifileid,"(a)") "        MAX_STEPS 300 #Maximum number of optimization steps"
        write(ifileid,"(a)") "        MAX_STEPSIZE 2.0 #Maximum stepsize used for the line search, may be reduced to stabilize line search for difficult initial geometries"
        write(ifileid,"(a)") "      &END DIIS"
        write(ifileid,"(a)") "    &END OPTIMIZE_BAND"
        write(ifileid,"(a)") "    &PROGRAM_RUN_INFO #Print basic information during running the BAND task"
        write(ifileid,"(a)") "      INITIAL_CONFIGURATION_INFO F #Print details of constructing initial configurations"
        write(ifileid,"(a)") "    &END PROGRAM_RUN_INFO"
        write(ifileid,"(a)") "    &CONVERGENCE_INFO #Print convergence criteria during the BAND run in [proj]-BAND**.out files"
        write(ifileid,"(a)") "    &END CONVERGENCE_INFO"
        write(ifileid,"(a)") "    &CONVERGENCE_CONTROL #Setup parameters to control the convergence criteria for BAND"
        write(ifileid,"(a)") "      MAX_DR 2E-4"
        write(ifileid,"(a)") "      MAX_FORCE 4.5E-4"
        write(ifileid,"(a)") "      RMS_DR 1E-4"
        write(ifileid,"(a)") "      RMS_FORCE 3E-4"
        write(ifileid,"(a)") "    &END CONVERGENCE_CONTROL"
        write(ifileid,"(a)") "    #Specify coordinates of reactant and product structures, and possibly one or more intermediate structure(s) between them"
        write(ifileid,"(a)") "    &REPLICA"
        write(ifileid,"(a)") "      COORD_FILE_NAME reactant.xyz"
        write(ifileid,"(a)") "    &END"
        write(ifileid,"(a)") "    &REPLICA"
        write(ifileid,"(a)") "      COORD_FILE_NAME intermediate.xyz"
        write(ifileid,"(a)") "    &END"
        write(ifileid,"(a)") "    &REPLICA"
        write(ifileid,"(a)") "      COORD_FILE_NAME product.xyz"
        write(ifileid,"(a)") "    &END"
        write(ifileid,"(a)") "  &END BAND"

    else if (itask==6.or.itask==13) then !MD or real-time TDDFT
        write(ifileid,"(a)") "  &MD"
        if (ithermostat==0.and.ibarostat==0) then
            write(ifileid,"(a)") "    ENSEMBLE NVE"
        else if (ithermostat>1.and.ibarostat==0) then
            write(ifileid,"(a)") "    ENSEMBLE NVT"
        else if (ithermostat==0) then
            if (ibarostat==1) write(ifileid,"(a)") "    ENSEMBLE NPE_F"
            if (ibarostat==2) write(ifileid,"(a)") "    ENSEMBLE NPE_I"
        else if (ithermostat>1) then
            if (ibarostat==1) write(ifileid,"(a)") "    ENSEMBLE NPT_F"
            if (ibarostat==2) write(ifileid,"(a)") "    ENSEMBLE NPT_I"
        end if
        if (itask==6) then !MD
            write(ifileid,"(a)") "    STEPS 200 #Number of steps to run"
            write(ifileid,"(a)") "    TIMESTEP 1.0 #Step size in fs. Decrease it properly for high temperature simulation"
        else if (itask==13) then !RT-TDDFT needs very small step
            write(ifileid,"(a)") "    STEPS 20000 #Number of steps to run"
            write(ifileid,"(a)") "    TIMESTEP 0.025 #Step size in fs"
        end if
        write(ifileid,"(a)") "    TEMPERATURE 298.15 #Initial and maintained temperature (K)"
        write(ifileid,"(a)") "#   COMVEL_TOL 0 #Uncomment this can remove translation motion of center-of-mass every step"
        if (ifPBC==0.or.PBCdir=="NONE") then
            write(ifileid,"(a)") "#   ANGVEL_TOL 0 #Uncomment this can remove overall rotation every step"
            write(ifileid,"(a)") "    ANGVEL_ZERO T #Eliminate overall rotation component from initial velocity"
        end if
        if (ithermostat>0) then
            write(ifileid,"(a)") "    &THERMOSTAT"
            if (ithermostat==1) then
                write(ifileid,"(a)") "      TYPE AD_LANGEVIN"
            else if (ithermostat==2) then
                write(ifileid,"(a)") "      TYPE CSVR"
                write(ifileid,"(a)") "      &CSVR"
                write(ifileid,"(a)") "        TIMECON 200 #Time constant in fs. Smaller/larger results in stronger/weaker temperature coupling"
                write(ifileid,"(a)") "      &END CSVR"
            else if (ithermostat==3) then
                write(ifileid,"(a)") "      TYPE GLE"
            else if (ithermostat==4) then
                write(ifileid,"(a)") "      TYPE NOSE"
            end if
            !if (nthermoatm<ncenter) then !Misleading
            !    write(ifileid,"(a)") "      &DEFINE_REGION"
            !    call outCP2K_LIST(ifileid,thermoatm(1:nthermoatm),nthermoatm,"        ","LIST")
            !    write(ifileid,"(a)") "      &END DEFINE_REGION"
            !end if
            write(ifileid,"(a)") "    &END THERMOSTAT"
        end if
        if (ibarostat/=0) then
            write(ifileid,"(a)") "    &BAROSTAT"
            write(ifileid,"(a)") "      PRESSURE 1.01325 #Initial and maintained pressure (bar)"
            write(ifileid,"(a)") "      TIMECON 1000 #Barostat time constant (fs)"
            if (ibarostat==1) write(ifileid,"(a)") "      VIRIAL XYZ #Relax the cell along which cartesian axes"
            write(ifileid,"(a)") "    &END BAROSTAT"
        end if
        if (itask==6) then
            write(ifileid,"(a)") "    &PRINT"
            write(ifileid,"(a)") "      &PROGRAM_RUN_INFO"
            write(ifileid,"(a)") "        &EACH"
            write(ifileid,"(a,i6,a)") "          MD",1," #Output frequency of MD information, 0 means never"
            write(ifileid,"(a)") "        &END EACH"
            write(ifileid,"(a)") "      &END PROGRAM_RUN_INFO"
            write(ifileid,"(a)") "    &END PRINT"
        end if
        write(ifileid,"(a)") "  &END MD"
    else if (itask==14) then !PINT
        write(ifileid,"(a)") "  &PINT"
        write(ifileid,"(a)") "    PROPAGATOR PIMD #Type of propagator: PIMD, CMD, RPMD"
        write(ifileid,"(a)") "    DT 0.5 #Stepsize in fs"
        write(ifileid,"(a)") "    P 8 #Number of beads"
        write(ifileid,"(a)") "    NUM_STEPS 1000 #Number of dynamics steps"
        write(ifileid,"(a)") "    TEMP 298.15 #Simulation temperature (K)"
        write(ifileid,"(a)") "    T_TOL 50.0 #Threshold for the oscillations of the temperature excedeed which the temperature is rescaled. 0 means no rescaling"
        write(ifileid,"(a)") "    TRANSFORMATION NORMAL #Coordinate transformation method: NORMAL or STAGE"
        write(ifileid,"(a)") "    HARM_INT NUMERIC #Integrator scheme for integrating the harmonic bead springs: EXACT or NUMERIC"
        write(ifileid,"(a)") "    NRESPA 1 #Number of RESPA steps for the bead for each MD step"
        write(ifileid,"(a)") "    &NOSE #Use Nose-Hoover chain thermostat"
        write(ifileid,"(a)") "      NNOS 3 #Nose-Hoover chain length"
        write(ifileid,"(a)") "    &END NOSE"
        write(ifileid,"(a)") "  &END PINT"
    end if
    if (natmcons>0) then
        write(ifileid,"(a)") "  &CONSTRAINT"
        write(ifileid,"(a)") "    &FIXED_ATOMS #Set atoms to be fixed"
        write(ifileid,"(a)") "      COMPONENTS_TO_FIX XYZ #Which fractional components will be fixed, can be X, Y, Z, XY, XZ, YZ, XYZ"
        call outCP2K_LIST(ifileid,atmcons(1:natmcons),natmcons,"      ","LIST")
        write(ifileid,"(a)") "    &END FIXED_ATOMS"
        write(ifileid,"(a)") "  &END CONSTRAINT"
    end if
    
    !https://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/PRINT.html
    
    !Control output frequency of various properties
    write(ifileid,"(a)") "  &PRINT"
    if (itask==3.or.itask==4.or.itask==7) then !Optimizing minimum, TS
        write(ifileid,"(a)") "    &TRAJECTORY"
        if (iMDformat==1) write(ifileid,"(a)") "      FORMAT xyz"
        if (iMDformat==2) write(ifileid,"(a)") "      FORMAT dcd"
        if (iMDformat==-2) write(ifileid,"(a)") "      FORMAT DCD_ALIGNED_CELL"
        if (iMDformat==3) write(ifileid,"(a)") "      FORMAT pdb"
        write(ifileid,"(a)") "    &END TRAJECTORY"
    else if (itask==6.or.itask==14) then !MD or PIMD
        write(ifileid,"(a)") "    &TRAJECTORY"
        write(ifileid,"(a)") "      &EACH"
        if (itask==6) write(ifileid,"(a,i4,a)") "        MD",nMDsavefreq," #Output frequency of coordinates, 0 means never"
        if (itask==14) write(ifileid,"(a,i4,a)") "        PINT",nMDsavefreq," #Output frequency of coordinates, 0 means never"
        write(ifileid,"(a)") "      &END EACH"
        if (iMDformat==1) write(ifileid,"(a)") "      FORMAT xyz"
        if (iMDformat==2) write(ifileid,"(a)") "      FORMAT dcd"
        if (iMDformat==-2) write(ifileid,"(a)") "      FORMAT DCD_ALIGNED_CELL"
        if (iMDformat==3) write(ifileid,"(a)") "      FORMAT pdb"
        write(ifileid,"(a)") "    &END TRAJECTORY"
        write(ifileid,"(a)") "    &VELOCITIES"
        write(ifileid,"(a)") "      &EACH"
        if (itask==6) write(ifileid,"(a,i6,a)") "        MD",0," #Output frequency of velocities, 0 means never"
        if (itask==14) write(ifileid,"(a,i6,a)") "        PINT",0," #Output frequency of velocities, 0 means never"
        write(ifileid,"(a)") "      &END EACH"
        write(ifileid,"(a)") "    &END VELOCITIES"
        write(ifileid,"(a)") "    &FORCES"
        write(ifileid,"(a)") "      &EACH"
        if (itask==6) write(ifileid,"(a,i6,a)") "        MD",0," #Output frequency of forces, 0 means never"
        if (itask==14) write(ifileid,"(a,i6,a)") "        PINT",0," #Output frequency of forces, 0 means never"
        write(ifileid,"(a)") "      &END EACH"
        write(ifileid,"(a)") "    &END FORCES"
    end if
    
    write(ifileid,"(a)") "    &RESTART"
    write(ifileid,"(a)") "      BACKUP_COPIES 0 #Maximum number of backing up restart file, 0 means never" !Do not generate annoying .restart.bak file
    !For other tasks, by default, restart file is updated every step. Only for MD it is default to 20, I explicitly provide option to change it
    if (itask==6.or.itask==14) then
        write(ifileid,"(a)") "      &EACH"
        if (itask==6) write(ifileid,"(a)") "        MD  1 #Frequency of updating last restart file, 0 means never"
        if (itask==14) write(ifileid,"(a)") "        PINT  1 #Frequency of updating last restart file, 0 means never"
        write(ifileid,"(a)") "      &END EACH"
    end if
    write(ifileid,"(a)") "    &END RESTART"
    
    !Control how to generate history .restart files
    !For GEO_OPT and MD, default is 500. For other tasks, default is every step
    !Because it is useless, so I completely suppress it
    if (itask==4.or.itask==6.or.itask==14) then !Cell opt, MD, PINT
        write(ifileid,"(a)") "    &RESTART_HISTORY OFF "
        write(ifileid,"(a)") "    &END RESTART_HISTORY"
    end if
    !write(ifileid,"(a)") "    &RESTART_HISTORY"
    !write(ifileid,"(a)") "      &EACH"
    !write(ifileid,"(a)") "        CELL_OPT 0 #How often a history .restart file is generated, 0 means never"
    !write(ifileid,"(a)") "      &END EACH"
    !write(ifileid,"(a)") "    &END RESTART_HISTORY"
    write(ifileid,"(a)") "  &END PRINT"
    write(ifileid,"(a)") "&END MOTION"
end if
if (itask==5) then
    write(ifileid,"(a)") "&VIBRATIONAL_ANALYSIS"
    write(ifileid,"(a)") "  DX 0.01 #Step size of finite difference. This is default (Bohr)"
    write(ifileid,"(a)") "  NPROC_REP 1 #Number of processors to be used per replica. This is default"
    write(ifileid,"(a)") "  TC_PRESSURE 101325 #1 atm. Pressure for calculate thermodynamic data (Pa)"
    write(ifileid,"(a)") "  TC_TEMPERATURE 298.15 #Temperature for calculate thermodynamic data (K)"
    write(ifileid,"(a)") "  THERMOCHEMISTRY #Print thermochemistry information (only valid for isolated systems)"
    if (ikpoint1==1.and.ikpoint2==1.and.ikpoint3==1.and.iTDDFT==0) then
        write(ifileid,"(a)") "  INTENSITIES T #Calculate IR/Raman intensities"
    else !Cannot produce moment/polarizability and hence intensities when k-points are considered or TDDFT
        write(ifileid,"(a)") "  INTENSITIES F #Calculate IR/Raman intensities"
    end if
    if (ifPBC==0.or.PBCdir=="NONE") then
        write(ifileid,"(a)") "  FULLY_PERIODIC F #If T, avoiding to project out rotation component from Hessian matrix"
    else
        write(ifileid,"(a)") "  FULLY_PERIODIC T #Avoiding to project out rotation component from Hessian matrix"
    end if
    if (ioutvibmol==1) then
        write(ifileid,"(a)") "  &PRINT"
        write(ifileid,"(a)") "    &MOLDEN_VIB #Output .mol (Molden file) for visualizing vibrational modes"
        write(ifileid,"(a)") "    &END MOLDEN_VIB"
        write(ifileid,"(a)") "  &END PRINT"
    end if
    write(ifileid,"(a)") "&END VIBRATIONAL_ANALYSIS"
end if

close(ifileid)

write(*,"(a)") " CP2K input file has been exported to "//trim(outname)
write(*,"(a)") " Note: Please cite original papers of Multiwfn in your work if this function brings you convenience!"
end subroutine

!----- Test if a given integer array contains contiguous numbers (index order is unimportant)
!If yes, ifconti=1, else =0
subroutine testidx_contiguous(array,narray,ifconti)
implicit real*8 (a-h,o-z)
integer narray,array(narray)
ifconti=0
do i=minval(array),maxval(array)
    if (all(array/=i)) return
end do
ifconti=1
end subroutine

!----- Output integer arrays in CP2K "LIST" convention
!list(1:nlist) are indices to be exported to ifileid. spacestr is string containing proper number of spaces in front of data, header is label string
subroutine outCP2K_LIST(ifileid,list,nlist,spacestr,header)
character(len=*) spacestr,header
character c80tmp*80
integer ifileid,nlist,ifconti
integer list(nlist)
call testidx_contiguous(list,nlist,ifconti)
if (ifconti==1) then
    write(c80tmp,*) maxval(list(1:nlist))
    if (nlist==1) then
        write(ifileid,"(a)") spacestr//header//' '//trim(adjustl(c80tmp))
    else
        write(ifileid,"(a,i7,'..',a)") spacestr//header//' ',minval(list(1:nlist)),trim(adjustl(c80tmp))
    end if
else
    iline=1
    i=0
    do while(.true.) !Change line ended with \ when outputting every 12 terms
        if (iline==1) then
            write(ifileid,"(a)",advance='no') spacestr//header//' '
        else
            write(ifileid,"(a)",advance='no') spacestr
        end if
        iline=iline+1
        if (i*12<nlist-12) then !Not the last line, the remaining terms is more than 12
            write(ifileid,"(12i6,' \')") list(12*i+1:12*i+12)
            i=i+1
        else !The last line
            write(ifileid,"(12i6)") list(12*i+1:nlist)
            exit
        end if
    end do
end if
end subroutine



!!---------- Automatically set proper vacuum sizes
!For NMR (9), polar(10), RT-TDDFT(13), use ~20% larger size
subroutine determine_vacuumsize(itask,iPSOLVER,vacsizex,vacsizey,vacsizez,icentering)
use defvar
integer itask,iPSOLVER,icentering
real*8 vacsizex,vacsizey,vacsizez,buff
if (iPSOLVER==1) then !Usually adequate for PERIODIC
    buff=5
    if (itask==9.or.itask==10.or.itask==13) buff=6
    vacsizex=buff/b2a
    vacsizey=buff/b2a
    vacsizez=buff/b2a
else if (iPSOLVER==2) then !ANALYTIC converges quite slow
    buff=10
    if (itask==9.or.itask==10.or.itask==13) buff=12
    vacsizex=buff/b2a
    vacsizey=buff/b2a
    vacsizez=buff/b2a
else if (iPSOLVER==3) then !MT needs vaccum size in each side is >= half of system
    buff=4 !4 A is extension distance for electron tail in each side
    if (itask==9.or.itask==10.or.itask==13) buff=5
    vacsizex=((maxval(a%x)-minval(a%x))*b2a+2*buff)/2/b2a
    vacsizey=((maxval(a%y)-minval(a%y))*b2a+2*buff)/2/b2a
    vacsizez=((maxval(a%z)-minval(a%z))*b2a+2*buff)/2/b2a
else if (iPSOLVER==4) then
    icentering=1
    buff=3.5D0 !WAVELET converges quite fast, 3.5 A is adequate for any case
    if (itask==9.or.itask==10.or.itask==13) buff=4.2D0
    vacsizex=buff/b2a 
    vacsizey=buff/b2a
    vacsizez=buff/b2a
end if
end subroutine



!!--------- Convert band structure file of CP2K (.bs) to multiple column file so that band map can be directly plotted by Origin or gnuplot
subroutine CP2K_BS
use defvar
use util
implicit real*8 (a-h,o-z)
character c200tmp*200,c10tmp*10
integer,parameter :: nSPmax=1000,nkpmax=10000
character SPlabel(nSPmax)*10 !Label of special points
real*8 SPvec(3,nSPmax) !XYZ of special points
integer*2 SPpath(nSPmax) !The path number that this special point belongs to. Each path consists of connected special points
integer*2 SPkp(nSPmax) !Index of k-point that each special point corresponds to
real*8 tmpvec(3),kpvec(3,nkpmax)
integer kpplot(nkpmax) !kpplot(i) is k-point index for plotting band map of actual k-point i (two neighbouring k-points with same coordinate is merged to single index in kpplot. kpplot corresponds to X-axis position of band structure map)
E_VBT=-99999 !Valence band top
E_CBB=99999 !Conduction band bottom

eshift=0 !Don't shift here, because plotting band structure module is directly able to do shift
!write(*,*) "Input value for shifting energy levels, e.g. 3.42"
!write(*,*) "If inputting ENTER button directly, 0 will be used"
!read(*,"(a)") c200tmp
!if (c200tmp==" ") then
!    eshift=0
!else
!    read(c200tmp,*) eshift
!end if

open(10,file=filename,status="old")

iopsh=0 !Closed-shell
call loclabel(10,"Spin 2",iopsh)
if (iopsh==0) then
    open(11,file="BS_occ.txt",status="replace")
    open(12,file="BS_vir.txt",status="replace")
else if (iopsh==1) then
    open(11,file="BS_occ.txt",status="replace")
    open(12,file="BS_vir.txt",status="replace")
    open(13,file="BS_occ_B.txt",status="replace")
    open(14,file="BS_vir_B.txt",status="replace")
end if

rewind(10)
nSP=0 !Number of currently loaded special points
nkp=0 !Current k-point index
nkpplot=0
ipath=1 !Current path index
xkp_last=-999
ykp_last=-999
zkp_last=-999
do while(.true.)
    !Loop all sets
    read(10,"(a)",iostat=ierror) c200tmp
    if (ierror/=0.or.c200tmp==" ") exit
    if (index(c200tmp,"# Set")/=0) then
        !read(c200tmp,"(5x,i2,1x,i2,16x,i3,10x,i3)") iset,nSPthis,nkpthis,nlevelthis !This often does not work, because the value columns are not fixed
        ipos=index(c200tmp,'t')+1
        jpos=index(c200tmp,':')-1
        read(c200tmp(ipos:jpos),*) iset
        ipos=index(c200tmp,':')+1
        read(c200tmp(ipos:),*) nSPthis
        ipos=index(c200tmp,',')+1
        read(c200tmp(ipos:),*) nkpthis
        ipos=index(c200tmp,',',back=.true.)+1
        read(c200tmp(ipos:),*) nlevelthis
        write(*,"(a,i3,'...')") " Loading &KPOINT_SET",iset
    end if
    
    !Read position and label of special points in current set
    do iSP=1,nSPthis
        read(10,"(a)") c200tmp
        read(c200tmp(24:),*) tmpvec(:),c10tmp
        if (iSP==1.and.nSP/=0) then !A new set but not the first set
            if (c10tmp/=SPlabel(nSP)) ipath=ipath+1 !The first special point in this set has different label to the last one, so begins a new path
        end if
        iadd=0
        if (nSP==0) then
            iadd=1
        else if (c10tmp/=SPlabel(nSP)) then
            iadd=1
        end if
        if (iadd==1) then !Add this special point to list
            nSP=nSP+1
            SPvec(:,nSP)=tmpvec(:)
            SPlabel(nSP)=c10tmp
            SPpath(nSP)=ipath
        end if
    end do
    
    !Loop k-points between neighbouring special points in current set
    do ikp=1,nkpthis
        !Closed-shell or alpha part of open-shell
        read(10,"(a)") c200tmp
        read(10,*)
        read(c200tmp(25:),*) xkp,ykp,zkp !XYZ of k-point
        !write(*,"(2i5,6f10.6)") iset,ikp,xkp,ykp,zkp,xkp_last,ykp_last,zkp_last
        
        !The first point of this set is different to final point of the last set, make they share the same k-point index
        if (ikp==1.and.iset>1.and.(xkp/=xkp_last.or.ykp/=ykp_last.or.zkp/=zkp_last)) nkpplot=nkpplot-1
        
        if (nkp>1.and.xkp==xkp_last.and.ykp==ykp_last.and.zkp==zkp_last) then !This k-point is identical to the last one, so skip this one
            call skiplines(10,nlevelthis)
            if (iopsh==1) call skiplines(10,nlevelthis+2)
            cycle
        else !This k-point will be actually loaded
            nkp=nkp+1
            kpvec(1,nkp)=xkp
            kpvec(2,nkp)=ykp
            kpvec(3,nkp)=zkp
            nkpplot=nkpplot+1
            kpplot(nkp)=nkpplot
            xkp_last=xkp
            ykp_last=ykp
            zkp_last=zkp
        end if
        !write(11,"(i4,3f8.4)",advance="no") nkpplot,xkp,ykp,zkp !For checking xyz of k-point
        !write(12,"(i4,3f8.4)",advance="no") nkpplot,xkp,ykp,zkp
        write(11,"(i4)",advance="no") nkpplot
        write(12,"(i4)",advance="no") nkpplot
        !Loop energy levels in current k-point
        do ilevel=1,nlevelthis
            read(10,*) idx,ene,occ
            if (occ/=0) then
                write(11,"(f16.8)",advance="no") ene+eshift
                if (ene>E_VBT) E_VBT=ene
            else
                write(12,"(f16.8)",advance="no") ene+eshift
                if (ene<E_CBB) E_CBB=ene
            end if
        end do
        write(11,*)
        write(12,*)
        
        !Beta part
        if (iopsh==1) then
            read(10,*);read(10,*)
            write(13,"(i4)",advance="no") nkpplot
            write(14,"(i4)",advance="no") nkpplot
            do ilevel=1,nlevelthis
                read(10,*) idx,ene,occ
                if (occ/=0) then
                    write(13,"(f16.8)",advance="no") ene+eshift
                    if (ene>E_VBT) E_VBT=ene
                else
                    write(14,"(f16.8)",advance="no") ene+eshift
                    if (ene<E_CBB) E_CBB=ene
                end if
            end do
            write(13,*)
            write(14,*)
        end if
    end do
    nSP_last=nSP
end do

!Determine number of occupied and virtual levels
nocclevel=0 !Number of occupied alpha level
call loclabel(10,"#   Band")
read(10,*)
do ilevel=1,nlevelthis
    read(10,*) itmp,tmp,occ
    if (occ/=0) then
        nocclevel=nocclevel+1
    else
        exit
    end if
end do
nvirlevel=nlevelthis-nocclevel
if (iopsh==1) then
    nocclevelB=0 !Number of occupied beta level
    call loclabel(10,"#   Band")
    read(10,*)
    call loclabel(10,"#   Band",irewind=0) !Locate to spin 2
    read(10,*)
    do ilevel=1,nlevelthis
        read(10,*) itmp,tmp,occ
        if (occ/=0) then
            nocclevelB=nocclevelB+1
        else
            exit
        end if
    end do
    nvirlevelB=nlevelthis-nocclevelB
end if

close(10)
close(11)
close(12)
if (iopsh==1) then
    close(13)
    close(14)
end if

write(*,*)
write(*,*) "List of k-points: (Original index, actual index, X, Y, Z)"
do ikp=1,nkp
    write(*,"(2i4,3x,3f16.12)") ikp,kpplot(ikp),kpvec(:,ikp)
end do

write(*,*)
write(*,*) "Determine correspondence between special points and original k-point index..."
SPkp(:)=0
do iSP=1,nSP
    if (iSP==1) then
        ibeg=1
    else
        ibeg=SPkp(iSP-1)+1
    end if
    !write(*,*) iSP,ibeg,nkp
    do ikp=ibeg,nkp
        if (SPvec(1,iSP)==kpvec(1,ikp).and.SPvec(2,iSP)==kpvec(2,ikp).and.SPvec(3,iSP)==kpvec(3,ikp)) then
            SPkp(iSP)=ikp
            exit
        end if
    end do
    if (ikp==nkp+1) then
        write(*,"(a,i6,1x,'(',a,')','!')") " Warning: Unable to determine k-point attribution of special point",iSP,trim(SPlabel(iSP))
        write(*,"(' X,Y,Z of this special point:',3f10.6)") SPvec(:,iSP)
        write(*,*) "Press ENTER button to continue"
        read(*,*)
    else
        write(*,"(a,i4,1x,'(',a,') at',3f10.6,' attributes to k-point',i6)") " SP",iSP,SPlabel(iSP),SPvec(:,iSP),SPkp(iSP)
    end if
end do

write(*,"(/,' Number of total levels:',i7)") nlevelthis
!write(*,*)
!write(*,"(' Valence band top:      ',f16.8,' eV')") E_VBT
!write(*,"(' Conduction band bottom:',f16.8,' eV')") E_CBB
!write(*,"(' Band gap:              ',f16.8,' eV')") E_CBB-E_VBT
write(*,*)
write(*,*) "Special points for plotting band structure:"
write(*,*) "SP#   kp#   Label                Coordinates          Path"
do iSP=1,nSP
    if (iSP>1) then
        if (SPpath(iSP)/=SPpath(iSP-1)) write(*,*)
    end if
    write(*,"(i4,i5,4x,a,3f10.6,i5)") iSP,kpplot(SPkp(iSP)),SPlabel(iSP),SPvec(:,iSP),SPpath(iSP)
end do

!Output BS_info.txt
open(10,file="BS_info.txt",status="replace")
write(10,"('iopsh',i3)") iopsh
write(10,"('N_level',i8)") nlevelthis
if (iopsh==0) then
    write(10,"('N_occ',i8)") nocclevel
    write(10,"('N_vir',i8)") nvirlevel
else
    write(10,"('N_occ',2i8)") nocclevel,nocclevelB
    write(10,"('N_vir',2i8)") nvirlevel,nvirlevelB
end if
write(10,"('N_SP',i5)") nSP
do iSP=1,nSP
    write(10,"(i4,i5,4x,a,3f10.6,i5)") iSP,kpplot(SPkp(iSP)),SPlabel(iSP),SPvec(:,iSP),SPpath(iSP)
end do
close(10)

write(*,"(/,a)") " Basic information of band structure has been exported to BS_info.txt in current folder"
if (iopsh==0) then
    write(*,"(a)") " Occupied and virtual levels of all k-points involved in band structure map have been written to BS_occ.txt and BS_vir.txt, respectively"
else if (iopsh==1) then
    write(*,"(a)") " Alpha occupied and virtual levels of all k-points involved in band structure map have been written to BS_occ.txt and BS_vir.txt, respectively"
    write(*,"(a)") " Beta occupied and virtual levels of all k-points involved in band structure map have been written to BS_occ_B.txt and BS_vir_B.txt, respectively"
end if

!Automatically enter band plot interface
write(*,*)
write(*,*) "Now enter the interface of plotting band structure map"
call plotBS
end subroutine




!!------- Plot band structure based on BS_info.txt, BS_occ/vir.txt
subroutine plotBS
use dislin
use util
use defvar
implicit real*8 (a-h,o-z)
character c80tmp*80,outname*80
integer nocc,noccB !Number of alpha/beta occupied levels. For closed-shell, nocc is total number of occupied levels
integer nvir,nvirB !Number of virtual levels
integer nSP !Number of special points (SP)
integer,allocatable :: SPkp(:) !The k-point index in the whole path that each SP corresponds to
integer,allocatable :: SPpath(:) !The path index that each SP belongs to
character,allocatable :: SPlabel(:)*15,SPlabel_final(:)*15 !SP labels and its plotting form
real*8,allocatable :: SPxyz(:,:) !SP(1/2/3,i) is X,Y,Z of SP i in reciprocal space
real*8,allocatable :: SPXpos(:) !Position in X-axis of band structure map of every SP
real*8,allocatable :: kpXpos(:) !Position in X-axis of band structure map of every k-point
integer npath !Number of paths
integer,allocatable :: pathnkp(:) !Number of k-points of each path
integer,allocatable :: pathkpbeg(:),pathkpend(:) !Index of beginning and ending k-points of each path
real*8,allocatable :: occene(:,:,:),virene(:,:,:) !occene(ikp,ilevel,ipath) is energy of "ilevel" at "ikp" point of "ipath" path
real*8,allocatable :: occeneB(:,:,:),vireneB(:,:,:) !Beta part
!Plotting parameters
integer :: ivertline=1,iedgeline=1,ieneline=0,thkcurve=4,textsize=50
real*8 :: BSxyratio=0.618D0
real*8 :: linex(2),liney(2)
integer :: icurveclr_occ=5,icurveclr_vir=15,icurveclr_occB=11,icurveclr_virB=1,ishowtype=0
integer :: ivertlineclr=14,ihorilineclr=6,thkvertline=2,thkhoriline=4,thkaxis=4
eneshift=0
eneshiftB=0

write(*,*)

!Load BS_info.txt
inquire(file="BS_info.txt",exist=alive)
if (.not.alive) then
	write(*,*) "Cannot find BS_info.txt in current folder! Press ENTER button to exit"
    read(*,*)
    return
end if
write(*,*) "Loading BS_info.txt..."
open(10,file="BS_info.txt",status="old")
read(10,*) c80tmp,iopsh
read(10,*) c80tmp,nlevel
if (iopsh==0) then
    read(10,*) c80tmp,nocc
else
    read(10,*) c80tmp,nocc,noccB
end if
if (iopsh==0) then
    read(10,*) c80tmp,nvir
else
    read(10,*) c80tmp,nvir,nvirB
end if
read(10,*) c80tmp,nSP
allocate(SPkp(nSP),SPpath(nSP),SPlabel(nSP),SPlabel_final(nSP),SPxyz(3,nSP),SPXpos(nSP))
do iSP=1,nSP
    read(10,*) itmp,SPkp(iSP),SPlabel(iSP),SPxyz(:,iSP),SPpath(iSP)
end do
close(10)

!Generate SPXpos
!Total X-axis range of BS map is normalized to 1.0, and thereby we determine position of SP in X-axis. Then determining position of kp in X-axis
totdist=0
SPXpos(1)=0
do iSP=1,nSP-1
    dist=dsqrt(sum((SPxyz(:,iSP+1)-SPxyz(:,iSP))**2))
    if (SPpath(iSP)==SPpath(iSP+1)) then !They are in the same path
        totdist=totdist+dist
        SPXpos(iSP+1)=SPXpos(iSP)+dist
    else
        SPXpos(iSP+1)=SPXpos(iSP)
    end if
end do
!write(*,*) totdist
!do iSP=1,nSP
!    write(*,"(i5,f12.6)") iSP,SPXpos(iSP)
!end do
SPXpos(:)=SPXpos(:)/totdist

!Generate kpXpos
!k-points are evenly distributed between two special points
nkp=SPkp(nSP)
allocate(kpXpos(nkp))
do iSP=1,nSP-1
    if (SPkp(iSP+1)==SPkp(iSP)) cycle
    nspc=SPkp(iSP+1)-SPkp(iSP)
    dist=SPXpos(iSP+1)-SPXpos(iSP)
    spc=dist/nspc
    do itmp=1,nspc+1
        ikp=SPkp(iSP)+itmp-1
        kpXpos(ikp)=SPXpos(iSP)+spc*(itmp-1)
    end do
end do
!do ikp=1,nkp
!    write(*,"(i5,f12.6)") ikp,kpXpos(ikp)
!end do

!Calculate number of k-points in each path
npath=SPpath(nSP)
allocate(pathnkp(npath),pathkpbeg(npath),pathkpend(npath))
ibeg=1
do iSP=2,nSP
    ipath=SPpath(iSP)
    jpath=SPpath(iSP-1)
    if (ipath/=jpath) then
        pathkpbeg(jpath)=ibeg
        pathkpend(jpath)=SPkp(iSP-1)
        ibeg=SPkp(iSP)
    else if (iSP==nSP) then
        pathkpbeg(ipath)=ibeg
        pathkpend(ipath)=SPkp(nSP)
    end if
end do
do ipath=1,npath
    pathnkp(ipath)=pathkpend(ipath)-pathkpbeg(ipath)+1
    !write(*,"(' #Path',i3,'   From kp',i5,' to kp',i5,'   Number of kps:',i5)") ipath,pathkpbeg(ipath),pathkpend(ipath),pathnkp(ipath)
end do

!Load energy levels from BS_occ.txt and BS_vir.txt
allocate(occene(nkp,nocc,npath))
allocate(virene(nkp,nvir,npath))
inquire(file="BS_occ.txt",exist=alive)
if (.not.alive) then
	write(*,*) "Cannot find BS_occ.txt in current folder! Press ENTER button to exit"
    read(*,*)
    return
end if
inquire(file="BS_vir.txt",exist=alive)
if (.not.alive) then
	write(*,*) "Cannot find BS_vir.txt in current folder! Press ENTER button to exit"
    read(*,*)
    return
end if
write(*,*) "Loading BS_occ.txt and BS_vir.txt..."
open(10,file="BS_occ.txt",status="old")
open(11,file="BS_vir.txt",status="old")
do ipath=1,npath
    do ikp=pathkpbeg(ipath),pathkpend(ipath)
        read(10,*) itmp,occene(ikp,:,ipath)
        read(11,*) itmp,virene(ikp,:,ipath)
    end do
end do
close(10)
close(11)

!Load energy levels from BS_occ_B.txt and BS_vir_B.txt for beta part of open-shell case
if (iopsh==1) then
    write(*,*) "Loading BS_occ_B.txt and BS_vir_B.txt..."
    allocate(occeneB(nkp,noccB,npath))
    allocate(vireneB(nkp,nvirB,npath))
    inquire(file="BS_occ_B.txt",exist=alive)
    if (.not.alive) then
	    write(*,*) "Cannot find BS_occ_B.txt in current folder! Press ENTER button to exit"
        read(*,*)
        return
    end if
    inquire(file="BS_vir_B.txt",exist=alive)
    if (.not.alive) then
	    write(*,*) "Cannot find BS_vir_B.txt in current folder! Press ENTER button to exit"
        read(*,*)
        return
    end if
    open(10,file="BS_occ_B.txt",status="old")
    open(11,file="BS_vir_B.txt",status="old")
    do ipath=1,npath
        do ikp=pathkpbeg(ipath),pathkpend(ipath)
            read(10,*) itmp,occeneB(ikp,:,ipath)
            read(11,*) itmp,vireneB(ikp,:,ipath)
        end do
    end do
    close(10)
    close(11)
end if

ispin=1
if (iopsh==1) then
    write(*,*)
    write(*,*) "Plot which spin?"
    write(*,*) "1 Both spins"
    write(*,*) "2 Alpha spin"
    write(*,*) "3 Beta spin"
    read(*,*) ispin
end if

!Determine band edge and proper Y-axis range
do ipath=1,npath
    ibeg=pathkpbeg(ipath)
    iend=pathkpend(ipath)
    if (ipath==1) then
        E_HOCO=maxval(occene(ibeg:iend,:,ipath))
        E_LUCO=minval(virene(ibeg:iend,:,ipath))
        if (iopsh==1) then
            E_HOCO_B=maxval(occeneB(ibeg:iend,:,ipath))
            E_LUCO_B=minval(vireneB(ibeg:iend,:,ipath))
        end if
    else
        tmpval=maxval(occene(ibeg:iend,:,ipath))
        if (tmpval>E_HOCO) E_HOCO=tmpval
        tmpval=minval(virene(ibeg:iend,:,ipath))
        if (tmpval<E_LUCO) E_LUCO=tmpval
        if (iopsh==1) then
            tmpval=maxval(occeneB(ibeg:iend,:,ipath))
            if (tmpval>E_HOCO_B) E_HOCO_B=tmpval
            tmpval=minval(vireneB(ibeg:iend,:,ipath))
            if (tmpval<E_LUCO_B) E_LUCO_B=tmpval
        end if
    end if
end do
!When both alpha and beta are available, _A and _B distinguish spin. If consider only one spin, no suffix is used
if (iopsh==1) then
    if (ispin==1) then
        E_HOCO_A=E_HOCO
        E_LUCO_A=E_LUCO
        E_HOCO=max(E_HOCO,E_HOCO_B)
        E_LUCO=min(E_LUCO,E_LUCO_B)
    else if (ispin==3) then
        E_HOCO=E_HOCO_B
        E_LUCO=E_LUCO_B
    end if
end if

!Determine how many spins are considered in this time
if (iopsh==0.or.(iopsh==1.and.ispin/=1)) then !One spin
    nspin=1
else !Two spins
    nspin=2
    write(*,"(a)") " Note: By default, occupied and virtual levels of alpha spin will be drawn by red and crimson, respectively; &
    occupied and virtual levels of beta spin will be drawn by blue and dark blue, respectively"
    icurveclr_occ=1
    icurveclr_vir=11
    icurveclr_occB=3
    icurveclr_virB=15
end if

!Determine X-pos of HOCO and LUCO
do ipath=1,npath
    do ikp=pathkpbeg(ipath),pathkpend(ipath)
        if (nspin==1.and.ispin==3) then !Consider only beta spin
            do iocc=1,noccB
                if (occeneB(ikp,iocc,ipath)==E_HOCO) iHOCO=ikp
            end do
            do ivir=1,nvirB
                if (vireneB(ikp,ivir,ipath)==E_LUCO) iLUCO=ikp
            end do
        else !Alpha spin only, or alpha part of both spin
            do iocc=1,nocc
                if (occene(ikp,iocc,ipath)==E_HOCO) iHOCO=ikp
            end do
            do ivir=1,nvir
                if (virene(ikp,ivir,ipath)==E_LUCO) iLUCO=ikp
            end do
        end if
        if (nspin==2) then !Beta part of both spin case
            do iocc=1,noccB
                if (occeneB(ikp,iocc,ipath)==E_HOCO_B) iHOCOB=ikp
            end do
            do ivir=1,nvirB
                if (vireneB(ikp,ivir,ipath)==E_LUCO_B) iLUCOB=ikp
            end do
        end if
    end do
end do
!Show LUCO, HOCO, gap
if (nvir==0) then
    write(*,"(/,a)") " Warning: There is no virtual level! Please check if you have set ADDED_MOS in CP2K input file"
    write(*,*) "Press ENTER button to continue"
    read(*,*)
    write(*,*)
    if (nspin==1) then !One spin
        write(*,"(' K-point index of HOCO:',i5)") iHOCO
        write(*,"(/,a,f12.6,' eV')") " Energy of HOCO:",E_HOCO
    else !Two spins
        write(*,"(' K-point index of alpha HOCO:',i5)") iHOCO
        write(*,"(' K-point index of beta HOCO: ',i5)") iHOCOB
        write(*,"(/,a,f12.6,' eV')") " Energy of alpha HOCO:",E_HOCO_A
        write(*,"(a,f12.6,' eV')") " Energy of beta HOCO: ",E_HOCO_B
    end if
    E_LUCO=E_HOCO
else
    write(*,*)
    if (nspin==1) then !One spin
        write(*,"(' K-point index of HOCO:',i5)") iHOCO
        write(*,"(' K-point index of LUCO:',i5)") iLUCO
        write(*,"(/,a,f12.6,' eV')") " Energy of HOCO:",E_HOCO
        write(*,"(a,f12.6,' eV')") " Energy of LUCO:",E_LUCO
        E_gap=E_LUCO-E_HOCO
        if (iHOCO==iLUCO) then
            write(*,"(a,f12.6,' eV')") " Band gap (direct):",max(E_gap,0D0)
        else
            write(*,"(a,f12.6,' eV')") " Band gap (indirect):",max(E_gap,0D0)
        end if
        if (E_gap<=0) then
            write(*,*) "Note: This is a metal"
            iedgeline=0
        end if
    else !Two spins
        write(*,"(' K-point index of alpha HOCO:',i5)") iHOCO
        write(*,"(' K-point index of alpha LUCO:',i5)") iLUCO
        write(*,"(' K-point index of beta HOCO: ',i5)") iHOCOB
        write(*,"(' K-point index of beta LUCO: ',i5)") iLUCOB
        write(*,"(/,a,f12.6,' eV')") " Energy of alpha HOCO:",E_HOCO_A
        write(*,"(a,f12.6,' eV')") " Energy of alpha LUCO:",E_LUCO_A
        E_gap_A=E_LUCO_A-E_HOCO_A
        if (iHOCO==iLUCO) then
            write(*,"(a,f12.6,' eV')") " Alpha band gap (direct):",max(E_gap_A,0D0)
        else
            write(*,"(a,f12.6,' eV')") " Alpha band gap (indirect):",max(E_gap_A,0D0)
        end if
        write(*,"(a,f12.6,' eV')") " Energy of beta HOCO: ",E_HOCO_B
        write(*,"(a,f12.6,' eV')") " Energy of beta LUCO: ",E_LUCO_B
        E_gap_B=E_LUCO_B-E_HOCO_B
        if (iHOCOB==iLUCOB) then
            write(*,"(a,f12.6,' eV')") " Beta band gap (direct): ",max(E_gap_B,0D0)
        else
            write(*,"(a,f12.6,' eV')") " Beta band gap (indirect): ",max(E_gap_B,0D0)
        end if
        if (E_gap_A<=0.and.E_gap_B<=0) then
            write(*,*) "Note: This is a metal"
        else if (E_gap_A*E_gap_B<0) then
            write(*,*) "Note: This is a half-metal"
        end if
        iedgeline=0
    end if
end if

ymin=floor(E_HOCO-8)
ymax=ceiling(E_LUCO+5)
ystep=ceiling((ymax-ymin)/10D0)

!Convert special point labels to plotting form
do iSP=1,nSP
    SPlabel_final(iSP)=SPlabel(iSP)
    if (SPlabel(iSP)=="GAMMA") SPlabel_final(iSP)="\Gamma"
    iunder=index(SPlabel(iSP),'_') !Convert content after _ to subscript form
    if (iunder/=0) then
        write(SPlabel_final(iSP),*) SPlabel(iSP)(:iunder-1)//"$_{"//trim(SPlabel(iSP)(iunder+1:))//"}$"
    end if
end do

!!!! Enter interface
graphformat="pdf"
do while(.true.)
    write(*,*)
    call menutitle("Plotting band structure",10,1)
    write(*,*) "-3 Export curve data to plain text file"
    write(*,*) "-2 Set format of saving image file, current: "//trim(graphformat)
    write(*,*) "-1 Return"
    write(*,*) "0 Plot band structure map on screen"
    write(*,*) "1 Save band structure map to image file in current folder"
    write(*,"(a,f7.1,' to',f7.1,' with step',f4.1,' eV')") " 2 Set range and step of Y-axis, current:",ymin,ymax,ystep
    if (ivertline==0) write(*,*) "3 Toggle showing vertical lines highlighting special points, current: No"
    if (ivertline==1) write(*,*) "3 Toggle showing vertical lines highlighting special points, current: Yes"
    if (iedgeline==0) write(*,*) "4 Toggle showing horizontal lines highlighting band edges, current: No"
    if (iedgeline==1) write(*,*) "4 Toggle showing horizontal lines highlighting band edges, current: Yes"
    if (ieneline==0) write(*,*) "5 Toggle showing horizontal lines highlighting specific level, current: No"
    if (ieneline==1) write(*,*) "5 Toggle showing horizontal lines highlighting specific level, current: Yes"
    write(*,*) "6 Set colors"
    write(*,*) "7 Set thickness"
    write(*,"(a,f8.3)") " 8 Set ratio of X and Y axes, current:",BSxyratio
    write(*,"(a,i3)") " 9 Set label size, current:",textsize
    if (nspin==1) then !One spin
        write(*,"(a,f8.3,' eV')") " 10 Set shift of energy levels, current:",eneshift
    else !Two spins
        write(*,"(a,2f8.3,' eV')") " 10 Set shift of energy levels, current (Alpha and Beta):",eneshift,eneshiftB
    end if
    if (ishowtype==0) write(*,*) "11 Choose level type to be plotted, current: Occupied and virtual"
    if (ishowtype==1) write(*,*) "11 Choose level type to be plotted, current: Occupied"
    if (ishowtype==2) write(*,*) "11 Choose level type to be plotted, current: Virtual"
    read(*,*) isel
    
    if (isel==-3) then
        write(*,*) "Note: Range of X-axis of structure map is 0 to 1"
        do iSP=1,nSP
            write(*,"(' Special point',i5,4x,a,'X-Position:',f7.4,'    Path index:',i3)") iSP,SPlabel(iSP),SPXpos(iSP),SPpath(iSP)
        end do
        write(*,*)
        do ipath=1,npath
            write(c80tmp,*) ipath
            outname="path"//trim(adjustl(c80tmp))//".txt"
            write(*,*) "Outputting "//trim(outname)
            !non-spin or alpha spin
            open(10,file=outname,status="replace")
            if (nspin==1.and.ispin==3) then !Beta only case
                do ikp=pathkpbeg(ipath),pathkpend(ipath)
                    write(10,"(f7.4)",advance="no") kpXpos(ikp)
                    do iorb=1,noccB
                        write(10,"(f12.6)",advance="no") occeneB(ikp,iorb,ipath)+eneshift
                    end do
                    do iorb=1,nvirB
                        write(10,"(f12.6)",advance="no") vireneB(ikp,iorb,ipath)+eneshift
                    end do
                    write(10,*)
                end do
            else !Alpha only or alpha part of both spin
                do ikp=pathkpbeg(ipath),pathkpend(ipath)
                    write(10,"(f7.4)",advance="no") kpXpos(ikp)
                    do iorb=1,nocc
                        write(10,"(f12.6)",advance="no") occene(ikp,iorb,ipath)+eneshift
                    end do
                    do iorb=1,nvir
                        write(10,"(f12.6)",advance="no") virene(ikp,iorb,ipath)+eneshift
                    end do
                    write(10,*)
                end do
            end if
            close(10)
            if (nspin==2) then !Beta spin
                outname="path"//trim(adjustl(c80tmp))//"_B.txt"
                write(*,*) "Outputting "//trim(outname)
                open(10,file=outname,status="replace")
                do ikp=pathkpbeg(ipath),pathkpend(ipath)
                    write(10,"(f7.4)",advance="no") kpXpos(ikp)
                    do iorb=1,noccB
                        write(10,"(f12.6)",advance="no") occeneB(ikp,iorb,ipath)+eneshiftB
                    end do
                    do iorb=1,nvirB
                        write(10,"(f12.6)",advance="no") vireneB(ikp,iorb,ipath)+eneshiftB
                    end do
                    write(10,*)
                end do
                close(10)
            end if
        end do
        if (iopsh==0) then
            write(*,"(/,a)") " Done! Curve data of each path has been exported to current folder with ""path"" prefix"
        else
            write(*,"(/,a)") " Done! Curve data of each path has been exported to current folder with ""path"" prefix, &
            the files with and without ""_B"" suffix correspond to alpha and beta spins, respectively"
        end if
        write(*,"(a)") " Each row corresponds to a k-point, the first column is X-position of the k-points in band structure map, other columns correspond to orbital energies"
        if (eneshift/=0) write(*,*) "The exported orbital energies are the values after shifting"
    else if (isel==-2) then
        call setgraphformat
    else if (isel==-1) then
        return
    else if (isel==2) then
        write(*,*) "Set lower limit, upper limit, step size of Y-axis, e.g. -8,5,1"
        read(*,*) ymin,ymax,ystep
    else if (isel==3) then
        if (ivertline==0) then
            ivertline=1
        else
            ivertline=0
        end if
    else if (isel==4) then
        if (iedgeline==0) then
            iedgeline=1
        else
            iedgeline=0
        end if
    else if (isel==5) then
        if (ieneline==0) then
            ieneline=1
            write(*,*) "Input energy of the energy level to be highlighted, e.g. 0"
            read(*,*) E_highlight
        else
            ieneline=0
        end if
    else if (isel==6) then
        do while(.true.)
            write(*,*) "0 Return"
            if (iopsh==0) then
                write(*,*) "1 Set curve color of occupied levels, current: "//colorname(icurveclr_occ)
                write(*,*) "2 Set curve color of virtual levels, current: "//colorname(icurveclr_vir)
            else
                write(*,*) "1 Set curve color of alpha occupied levels, current: "//colorname(icurveclr_occ)
                write(*,*) "2 Set curve color of alpha virtual levels, current: "//colorname(icurveclr_vir)
                write(*,*) "3 Set curve color of beta occupied levels, current: "//colorname(icurveclr_occB)
                write(*,*) "4 Set curve color of beta virtual levels, current:  "//colorname(icurveclr_virB)
            end if
            write(*,*) "5 Set color of the vertical lines for highlighting, current: "//colorname(ivertlineclr)
            write(*,*) "6 Set color of the horizontal lines for highlighting, current: "//colorname(ihorilineclr)
            read(*,*) isel2
            if (isel2==0) then
                exit
            else if (isel2>=1.and.isel2<=4) then
                write(*,*) "Use which color for the curve?"
		        if (isel2==1) call selcolor(icurveclr_occ)
		        if (isel2==2) call selcolor(icurveclr_vir)
		        if (isel2==3) call selcolor(icurveclr_occB)
		        if (isel2==4) call selcolor(icurveclr_virB)
            else if (isel2==5) then
                write(*,*) "Use which color for the lines?"
                call selcolor(ivertlineclr)
            else if (isel2==6) then
                write(*,*) "Use which color for the lines?"
                call selcolor(ihorilineclr)
            end if
        end do
    else if (isel==7) then
        do while(.true.)
            write(*,*) "0 Return"
            write(*,"(a,i3)") " 1 Set thickness of curves, current:",thkcurve
            write(*,"(a,i3)") " 2 Set thickness of the vertical lines for highlighting, current:",thkvertline
            write(*,"(a,i3)") " 3 Set thickness of the horizontal lines for highlighting, current:",thkhoriline
            write(*,"(a,i3)") " 4 Set thickness of axis, current:",thkaxis
            read(*,*) isel2
            if (isel2==0) then
                exit
            else if (isel2==1) then
                write(*,*) "Input thickness, e.g. 3"
                read(*,*) thkcurve
            else if (isel2==2) then
                write(*,*) "Input thickness, e.g. 3"
                read(*,*) thkvertline
            else if (isel2==3) then
                write(*,*) "Input thickness, e.g. 3"
                read(*,*) thkhoriline
            else if (isel2==4) then
                write(*,*) "Input thickness, e.g. 3"
                read(*,*) thkaxis
            end if
        end do
    else if (isel==8) then
        write(*,*) "Input ratio of Y-axis : X-axis, e.g. 0.618"
        read(*,*) BSxyratio
    else if (isel==9) then
        write(*,*) "Input text size, e.g. 60"
        read(*,*) textsize
    else if (isel==10) then
        if (nspin==1) then
            write(*,*) "Input energy shift value (can be negative) in eV, e.g. 0.23"
        else
            write(*,"(a)") " Input energy shift value (can be negative) in eV for alpha and beta spins, respectively"
            write(*,*) "For example 0.23,0.45"
        end if
        write(*,*) "If you input H, then the levels will be shifted so that HOCO at 0 eV"
        read(*,"(a)") c80tmp
        if (c80tmp=='H'.or.c80tmp=='h') then
            if (nspin==1) then
                eneshift=-E_HOCO
            else
                write(*,*) "1 Shift levels of both spins by -E_HOCO(alpha)"
                write(*,*) "2 Shift levels of both spins by -E_HOCO(beta)"
                write(*,*) "3 Shift alpha levels by -E_HOCO(alpha) and shift beta levels by -E_HOCO(beta)"
                read(*,*) itype
                if (itype==1) then
                    eneshift=-E_HOCO_A
                    eneshiftB=-E_HOCO_A
                else if (itype==2) then
                    eneshift=-E_HOCO_B
                    eneshiftB=-E_HOCO_B
                else if (itype==3) then
                    eneshift=-E_HOCO_A
                    eneshiftB=-E_HOCO_B
                end if
            end if
        else
            if (nspin==1) then !One spin
                read(c80tmp,*) eneshift
            else !Two spins
                read(c80tmp,*) eneshift,eneshiftB
            end if
        end if
    else if (isel==11) then
        write(*,*) "Choose the level type to be plotted"
        write(*,*) "0 Occupied and virtual levels"
        write(*,*) "1 Occupied levels only"
        write(*,*) "2 Virtual levels only"
        read(*,*) ishowtype
    end if
    if (isel/=0.and.isel/=1) cycle
    
    !! Start plotting
    call SCRMOD('REVERSE')
    CALL PAGE(2970,2100)
    CALL setxid(0,'NONE') !If we don't set this, after we draw a graph embedded in GUI (e.g. relif map), curve map will not be shown 
    if (isel==0) then
	    call METAFL('xwin')
	    call window(200,100,900,600)
    else if (isel==1) then
        CALL IMGFMT("RGB")
	    call METAFL(graphformat)
	    call winsiz(graph1Dwidth,graph1Dheight)
    end if
    CALL DISINI
    call ERRMOD("ALL","OFF")
    if (isel==0) call WINTIT("Band structure, click right mouse button to continue...")
    call center
    nysize=nint(2300*BSxyratio)
    call AXSLEN(2300,nysize)
    if (nysize>1800) call AXSLEN(nint(2300/(nysize/1800D0)),1800)
    
    CALL HNAME(textsize+5)
    CALL height(textsize)
    CALL HWFONT
    call namdis(150,"X")
    CALL TICKS(0,"X")
    CALL NAME('k-points','X')
    call labels("NONE","X")
    CALL LABDIG(1,"Y")
    CALL TICPOS("REVERS","Y")
    CALL NAME('Energy (eV)','Y')
    
    !Draw axis
    CALL LINWID(thkaxis)
    CALL GRAF(0D0,1D0,0D0,1D0, ymin,ymax,ymin,ystep)
    
    !Draw special point labels
    do iSP=1,nSP
        if (iSP>1) then
            if (SPkp(iSP)==SPkp(iSP-1)) cycle !This label has already been plotted together last time
        end if
        c80tmp=SPlabel_final(iSP)
        if (iSP<nSP) then
            if (SPkp(iSP)==SPkp(iSP+1)) c80tmp=trim(SPlabel_final(iSP))//'|'//trim(SPlabel_final(iSP+1))
        end if
        CALL TEXMOD("ON")
        CALL ADDLAB(trim(c80tmp),SPXpos(iSP),0,"X")
    end do
    
    if (ivertline==1) then !Highlight position of special points 
        CALL SOLID
        call setcolor(ivertlineclr)
        CALL LINWID(thkvertline)
        liney(1)=ymin
        liney(2)=ymax
        do iSP=1,nSP
            linex(:)=SPXpos(iSP)
	        CALL CURVE(linex,liney,2)
        end do
    end if
    if (iedgeline==1) then !Highlight LUCO and HOCO positions
        CALL DASH
        call setcolor(ihorilineclr)
        CALL LINWID(thkhoriline)
        linex(1)=0
        linex(2)=1
        liney(:)=E_LUCO+eneshift
	    CALL CURVE(linex,liney,2)
        liney(:)=E_HOCO+eneshift
	    CALL CURVE(linex,liney,2)
    end if
    if (ieneline==1) then !Highlight a level
        CALL DASH
        call setcolor(ihorilineclr)
        CALL LINWID(thkhoriline)
        linex(1)=0
        linex(2)=1
        liney(:)=E_highlight
	    CALL CURVE(linex,liney,2)
    end if
    
    !Draw curves
    CALL SOLID
    CALL LINWID(thkcurve)
    do ipath=1,npath
        nkptmp=pathnkp(ipath)
        ibeg=pathkpbeg(ipath)
        iend=pathkpend(ipath)
        !Only one spin case, plot total for closed-shell, or selected spin of open-shell
        if (ishowtype==0.or.ishowtype==1) then !Plot occupied levels
            call setcolor(icurveclr_occ)
            if (nspin==1.and.ispin==3) then !Only consider Beta
                do iocc=1,noccB
                    CALL CURVE(kpXpos(ibeg:iend),occeneB(ibeg:iend,iocc,ipath)+eneshift,nkptmp)
                end do
            else !Alpha only, or alpha part of both spin case
                do iocc=1,nocc
                    CALL CURVE(kpXpos(ibeg:iend),occene(ibeg:iend,iocc,ipath)+eneshift,nkptmp)
                end do
            end if
        end if
        if (ishowtype==0.or.ishowtype==2) then !Plot virtual levels
            call setcolor(icurveclr_vir)
            if (nspin==1.and.ispin==3) then !Only consider Beta
                do ivir=1,nvirB
                    CALL CURVE(kpXpos(ibeg:iend),vireneB(ibeg:iend,ivir,ipath)+eneshift,nkptmp)
                end do
            else !Alpha only, or alpha part of both spin case
                do ivir=1,nvir
                    CALL CURVE(kpXpos(ibeg:iend),virene(ibeg:iend,ivir,ipath)+eneshift,nkptmp)
                end do
            end if
        end if
        if (nspin==2) then !Two spins case, plot beta part
            if (ishowtype==0.or.ishowtype==1) then !Plot occupied levels
                call setcolor(icurveclr_occB)
                do iocc=1,noccB !Plot occupied levels
                    CALL CURVE(kpXpos(ibeg:iend),occeneB(ibeg:iend,iocc,ipath)+eneshiftB,nkptmp)
                end do
            end if
            if (ishowtype==0.or.ishowtype==2) then !Plot virtual levels
                call setcolor(icurveclr_virB)
                do ivir=1,nvirB !Plot virtual levels
                    CALL CURVE(kpXpos(ibeg:iend),vireneB(ibeg:iend,ivir,ipath)+eneshiftB,nkptmp)
                end do
            end if
        end if
    end do

    CALL DISFIN
    if (isel==1) write(*,"(a)") " Done! Graphical file has been saved to current folder with ""dislin"" prefix"
end do

end subroutine





!!---------- Obtain E(LUCO), E(HOCO), band gap and exact DOS curve, based on the output file containing energy levels of every k-point
!itype=1: Only show E(LUCO), E(HOCO), band gap
!itype=2: Evaluate DOS, need PRINT_LEVEL >= medium
subroutine CP2K_bandgap_DOS(itype)
use defvar
use util
implicit real*8 (a-h,o-z)
character c200tmp*200
real*8,allocatable :: eneall(:,:),eneallB(:,:),kpweight(:) !eneall(imo,ikp) is energy of imo at ikp k-point. kpweight is weight of k-point
real*8,allocatable :: xpos(:),DOS(:)
integer :: ispin=0 !0=Both, 1=Alpha, 2=Beta

open(10,file=filename,status="old")
call loclabel(10,"Spin 1",iopsh)

!Total or alpha spin
call loclabel(10,"Number of occupied orbitals:")
read(10,"(a)") c200tmp
call readaftersign_int(c200tmp,':',noccorb)
read(10,"(a)") c200tmp
call readaftersign_int(c200tmp,':',nallorb)
nvirorb=nallorb-noccorb
if (iopsh==1) then !Beta spin
    call loclabel(10,"Number of occupied orbitals:",irewind=0)
    read(10,"(a)") c200tmp
    call readaftersign_int(c200tmp,':',noccorbB)
    read(10,"(a)") c200tmp
    call readaftersign_int(c200tmp,':',nallorbB)
    nvirorbB=nallorbB-noccorbB
end if

if (itype==2) then
    call loclabel(10,"List of Kpoints",ifound)
    if (ifound==0) then
        write(*,*) "Error: Unable to find k-points information! PRINT_LEVEL should be >= medium"
        close(10)
        return
    end if
    read(10,"(a)") c200tmp
    read(c200tmp(70:),*) nkp
    read(10,*)
    allocate(kpweight(nkp))
    do ikp=1,nkp
        read(10,*) c200tmp,itmp,kpweight(ikp)
    end do
    write(*,"(a,f8.5)") " Sum of k-point weights:",sum(kpweight(:))
    allocate(eneall(nallorb,nkp))
    if (iopsh==1) allocate(eneallB(nallorbB,nkp))
end if

if (iopsh==0) then
    write(*,"(' Number of occupied levels:',i8)") noccorb
    write(*,"(' Number of virtual levels: ',i8)") nvirorb
else
    write(*,"(' Number of alpha occupied levels:',i8)") noccorb
    write(*,"(' Number of alpha virtual levels: ',i8)") nvirorb
    write(*,"(' Number of beta occupied levels: ',i8)") noccorbB
    write(*,"(' Number of beta virtual levels:  ',i8)") nvirorbB
end if

nkp=0
do while(.true.)
    call loclabel(10,"EIGENVALUES AND OCCUPATION NUMBERS FOR K POINT",ifound,0)
    if (ifound==0) then
        if (nkp==0) then
            write(*,"(a)") " Error: Unable to find ""EIGENVALUES AND OCCUPATION NUMBERS FOR K POINT..."" information, &
            to use this function you should ask CP2K to print orbital information at every k-point!"
            close(10)
            return
        else
            exit
        end if
    end if
    nkp=nkp+1
    call skiplines(10,3)
    !Load occupied orbitals
    do iorb=1,noccorb
        read(10,*) c200tmp,itmp,tmpval,ene
        if (nkp==1) then
            E_HOCO=ene
            iHOCOkp=nkp
        else
            if (ene>E_HOCO) then
                E_HOCO=ene
                iHOCOkp=nkp
            end if
        end if
        if (itype==2) eneall(iorb,nkp)=ene
    end do
    !Load virtual orbitals
    do iorb=1,nvirorb
        read(10,*) c200tmp,itmp,tmpval,ene
        if (nkp==1) then
            E_LUCO=ene
            iLUCOkp=nkp
        else
            if (ene<E_LUCO) then
                E_LUCO=ene
                iLUCOkp=nkp
            end if
        end if
        if (itype==2) eneall(noccorb+iorb,nkp)=ene
    end do
    if (iopsh==1) then !Load beta part
        call loclabel(10,"EIGENVALUES AND OCCUPATION NUMBERS FOR K POINT",ifound,0)
        call skiplines(10,3)
        do iorb=1,noccorbB
            read(10,*) c200tmp,itmp,tmpval,ene
            if (nkp==1) then
                E_HOCO_B=ene
                iHOCOkp_B=nkp
            else
                if (ene>E_HOCO_B) then
                    E_HOCO_B=ene
                    iHOCOkp_B=nkp
                end if
            end if
            if (itype==2) eneallB(iorb,nkp)=ene
        end do
        do iorb=1,nvirorbB
            read(10,*) c200tmp,itmp,tmpval,ene
            if (nkp==1) then
                E_LUCO_B=ene
                iLUCOkp_B=nkp
            else
                if (ene<E_LUCO_B) then
                    E_LUCO_B=ene
                    iLUCOkp_B=nkp
                end if
            end if
            if (itype==2) eneallB(noccorbB+iorb,nkp)=ene
        end do
    end if
end do

close(10)

write(*,*)
write(*,"(' Totally processed',i8,' k-points')") nkp
if (iopsh==0) then
    write(*,"(' HOCO is located at k-point',i8)") iHOCOkp
    if (nvirorb>0) write(*,"(' LUCO is located at k-point',i8)") iLUCOkp
    write(*,"(/,' HOCO energy:',f12.6' eV')") E_HOCO
    if (nvirorb==0) then
        write(*,*) "There is no virtual orbitals"
    else
        write(*,"(' LUCO energy:',f12.6' eV')") E_LUCO
        write(*,"(' Band gap:   ',f12.6' eV')") E_LUCO-E_HOCO
    end if
else
    write(*,"(' Alpha HOCO is located at k-point',i8)") iHOCOkp
    if (nvirorb>0) write(*,"(' Alpha LUCO is located at k-point',i8)") iLUCOkp
    write(*,"(' Beta HOCO is located at k-point',i8)") iHOCOkp_B
    if (nvirorbB>0) write(*,"(' Beta LUCO is located at k-point',i8)") iLUCOkp_B
    write(*,"(/,' Alpha HOCO energy:',f12.6' eV')") E_HOCO
    if (nvirorb==0) then
        write(*,*) "There is no alpha virtual orbitals"
    else
        write(*,"(' Alpha LUCO energy:',f12.6' eV')") E_LUCO
        write(*,"(' Alpha Band gap:   ',f12.6' eV')") E_LUCO-E_HOCO
    end if
    write(*,"(/,' Beta HOCO energy: ',f12.6' eV')") E_HOCO_B
    if (nvirorbB==0) then
        write(*,*) "There is no beta virtual orbitals"
    else
        write(*,"(' Beta LUCO energy: ',f12.6' eV')") E_LUCO_B
        write(*,"(' Beta Band gap:    ',f12.6' eV')") E_LUCO_B-E_HOCO_B
    end if
end if

if (itype==1) return

!Obtain exact DOS
DOS_FWHM=0.35D0
npt=1000
xlow=-7
xhigh=5
eneshift=-E_HOCO
write(*,*)
if (iopsh==0) then
    write(*,*) "Note: Shift of energy levels has been set to -E(HOCO)"
else
    write(*,*) "Note: Shift of energy levels has been set to -E(HOCO) of alpha spin"
end if

do while(.true.)
    write(*,*)
    call menutitle("Generate exact DOS",10,1)
    write(*,*) "0 Return"
    write(*,*) "1 Generate DOS curve and export to DOS.txt in current folder!"
    write(*,"(a,f8.4,' eV')") " 2 Set FWHM, current:",DOS_FWHM
    write(*,"(a,f10.4,' to',f10.4,' eV')") " 3 Set energy range, current: ",xlow,xhigh
    write(*,"(a,i6)") " 4 Set number of points, current:",npt
    write(*,"(a,f9.4,' eV')") " 5 Set shift of energy levels, current:",eneshift
    if (iopsh==1) then
        if (ispin==0) write(*,*) "6 Choose spin, current: Alpha+Beta"
        if (ispin==1) write(*,*) "6 Choose spin, current: Alpha"
        if (ispin==2) write(*,*) "6 Choose spin, current: Beta"
    end if
    read(*,*) isel
    
    if (isel==0) then
        return
    else if (isel==2) then
        write(*,*) "Input FWHM in eV, e.g. 0.3"
        read(*,*) DOS_FWHM
    else if (isel==3) then
        write(*,*) "Input lower and upper limit of DOS, e.g. -8,5"
        read(*,*) xlow,xhigh
    else if (isel==4) then
        write(*,*) "Input number of points, e.g. 2000"
        read(*,*) npt
    else if (isel==5) then
        write(*,*) "Input shift of energy levels, e.g. 0.34"
        if (iopsh==0) then
            write(*,*) "If input ""H"", then shift value will be -E(HOCO)"
        else
            if (ispin==0.or.ispin==1) write(*,*) "If input ""H"", then shift value will be -E(HOCO) of Alpha spin"
            if (ispin==2) write(*,*) "If input ""H"", then shift value will be -E(HOCO) of Beta spin"
        end if
        read(*,*) c200tmp
        if (c200tmp=='h'.or.c200tmp=='H') then
            if (iopsh==0) then
                eneshift=-E_HOCO
            else
                if (ispin==0.or.ispin==1) then
                    eneshift=-E_HOCO
                else
                    eneshift=-E_HOCO_B
                end if
            end if
        else
            read(c200tmp,*) eneshift
        end if
    else if (isel==6) then
        write(*,*) "Choose spin:"
        write(*,*) "0 Alpha+Beta"
        write(*,*) "1 Alpha"
        write(*,*) "2 Beta"
        read(*,*) ispin
    end if
    if (isel/=1) cycle
    
    write(*,*)
    write(*,*) "Generating DOS data..."
    allocate(xpos(npt),DOS(npt))
    step=(xhigh-xlow)/(npt-1)
    do ipt=1,npt
        xpos(ipt)=xlow+step*(ipt-1)
    end do
    DOS=0
    gauss_c=DOS_FWHM/2D0/sqrt(2*dlog(2D0))
    fac=2*gauss_c**2
    do ikp=1,nkp !Cycle k-point
        gauss_a=kpweight(ikp)/(gauss_c*sqrt(2D0*pi))
        if (iopsh==0.or.ispin==0.or.ispin==1) then
            do imo=1,nallorb !Cycle each orbital
                ene=eneall(imo,ikp)+eneshift
	            do ipt=1,npt !Cycle each point of DOS curve
		            tmp=gauss_a*dexp( -(xpos(ipt)-ene)**2/fac )
		            DOS(ipt)=DOS(ipt)+tmp
	            end do
            end do
        end if
        if (iopsh==1.and.(ispin==0.or.ispin==2)) then
            do imo=1,nallorbB !Cycle each orbital
                ene=eneallb(imo,ikp)+eneshift
	            do ipt=1,npt !Cycle each point of DOS curve
		            tmp=gauss_a*dexp( -(xpos(ipt)-ene)**2/fac )
		            DOS(ipt)=DOS(ipt)+tmp
	            end do
            end do
        end if
    end do
    open(11,file="DOS.txt",status="replace")
    do ipt=1,npt
        write(11,"(2f12.6)") xpos(ipt),DOS(ipt)
    end do
    close(11)
    write(*,*) "Done! DOS curve has been exported to DOS.txt in current folder"
    write(*,*) "Column 1: Energy (eV)"
    write(*,*) "Column 2: DOS"
    deallocate(xpos,DOS)
end do
end subroutine



!!---------- Load orbital energies printed by &PRINT/&MO/ENERGIES T
!Mostly used for OT calculation, for which the printed .molden doesn't contain orbital energies
subroutine CP2K_MOene_load
use defvar
use util
implicit real*8 (a-h,o-z)
character c2000tmp*2000,c200tmp*200,c80tmp*80

if (.not.allocated(CObasa)) then
    write(*,*) "Error: Orbital information must be available to use this function!"
    return
end if

write(*,*)
write(*,*) "1 Load orbital energies from output file with &PRINT/&MO/ENERGIES T"
write(*,*) "2 Load orbital energies from output file with &PRINT/&MO_CUBES"
read(*,*) itype

write(*,"(a)") " Input the path of CP2K output file containing orbital energies, e.g. D:\waifu.out"
do while(.true.)
	read(*,"(a)") c2000tmp
	inquire(file=c2000tmp,exist=alive)
	if (alive) exit
	write(*,*) "Cannot find the file, input again!"
end do
open(10,file=trim(c2000tmp),status="old")

if (itype==1) then !&PRINT/&MO/ENERGIES T
    if (wfntype==0.or.wfntype==3) then !Closed-shell
        call loclabelfinal(10," MO|  Index",ifound)
        read(10,*)
        imo=0
        do while(.true.)
            read(10,"(a)") c200tmp
            if (index(c200tmp,"Sum")/=0) then
                exit
            else
                imo=imo+1
                read(c200tmp,*) c80tmp,inouse,MOene(imo)
            end if
        end do
    else !Open-shell
        call loclabelfinal(10," MO| ALPHA EIGENVALUES",ifound)
        call skiplines(10,3)
        imo=0
        do while(.true.)
            read(10,"(a)") c200tmp
            if (index(c200tmp,"Sum")/=0) then
                exit
            else
                imo=imo+1
                read(c200tmp,*) c80tmp,inouse,MOene(imo)
            end if
        end do
        call loclabel(10," MO|  Index",ifound,0)
        read(10,*)
        imo=0
        do while(.true.)
            read(10,"(a)") c200tmp
            if (index(c200tmp,"Sum")/=0) then
                exit
            else
                imo=imo+1
                read(c200tmp,*) c80tmp,inouse,MOene(nbasis+imo)
            end if
        end do
    end if
    
else !&PRINT/&MO_CUBES
    write(*,*)
    if (wfntype==0.or.wfntype==3) then !Closed-shell
        write(*,*) "Load how many unoccupied orbital energies? e.g. 3"
        write(*,*) "If press ENTER button directly, no unoccupied orbital energies will be loaded"
        read(*,"(a)") c80tmp
        nvirload=0
        if (c80tmp/=" ") read(c80tmp,*) nvirload
        call loclabelfinal(10,"Eigenvalues of the occupied subspace spin",ifound)
        call skiplines(10,2)
        read(10,*) MOene(1:nint(nelec/2))
        if (nvirload/=0) then
            call loclabel(10,"------",ifound,0)
            read(10,*)
            read(10,*) c80tmp
            if (index(c80tmp,"OT|")==0) backspace(10)
            read(10,*) MOene(nint(nelec/2)+1:nint(nelec/2)+nvirload)
        end if
    else !Open-shell
        write(*,*) "Load how many unoccupied orbital energies for alpha and beta spins? e.g. 3,4"
        write(*,*) "If press ENTER button directly, no unoccupied orbital energies will be loaded"
        read(*,"(a)") c80tmp
        navirload=0
        nbvirload=0
        if (c80tmp/=" ") read(c80tmp,*) navirload,nbvirload
        call loclabelfinal(10,"Eigenvalues of the occupied subspace spin            1",ifound)
        call skiplines(10,2)
        read(10,*) MOene(1:nint(naelec)) !Alpha occupied
        call loclabel(10,"------",ifound,0)
        read(10,*)
        read(10,*) MOene(nbasis+1:nbasis+nint(nbelec)) !Beta occupied
        call loclabel(10,"------",ifound,0)
        read(10,*)
        if (navirload/=0) then
            read(10,*) c80tmp
            if (index(c80tmp,"OT|")==0) backspace(10)
            if (index(c80tmp,"WARNING")/=0) call skiplines(10,3) !May encounter "WARNING : did not converge in ot_eigensolver", and the subsequent two lines are detailed information    
            read(10,*) MOene(nint(naelec)+1:nint(naelec)+navirload)
        end if
        if (nbvirload/=0) then
            call loclabel(10,"------",ifound,0)
            read(10,*)
            read(10,*) c80tmp
            if (index(c80tmp,"OT|")==0) backspace(10)
            if (index(c80tmp,"WARNING")/=0) call skiplines(10,3)
            read(10,*) MOene(nbasis+nint(nbelec)+1:nbasis+nint(nbelec)+nbvirload)
        end if
    end if
end if
    
close(10)
write(*,*) "Loading finished!"
end subroutine





!!---------- Load orbital coefficients of specific k-point from CP2K output file with &DFT/&PRINT/&MO/COEFFICIENTS T    
subroutine CP2K_loadkpwfn
use defvar
use util
implicit real*8 (a-h,o-z)
character c2000tmp*2000,c80tmp*80,c200tmp*200
real*8,allocatable :: CObasa_tmp(:,:),CObasb_tmp(:,:)

write(*,"(/,a)") " Input the path of CP2K output file containing orbital coefficients, e.g. D:\sobereva.out"
do while(.true.)
	read(*,"(a)") c2000tmp
	inquire(file=c2000tmp,exist=alive)
	if (alive) exit
	write(*,*) "Cannot find the file, input again!"
end do

ikp=10
write(*,*)
write(*,*) "Load which k-point? e.g. 3"
read(*,*) ikp

open(10,file=trim(c2000tmp),status="old")
write(c80tmp,"(i6)") ikp
call loclabelfinal(10,"FOR K POINT "//trim(adjustl(c80tmp)),ifound)
if (ifound==0) then
    write(*,*) "Error: Unable to find orbital information for this k-point!"
    write(*,*) "Press ENTER button to return"
    read(*,*)
    close(10)
    return
end if

!Try to find k-point coordinate, available when PRINT_LEVEL >= MEDIUM
call loclabel(10," BRILLOUIN| Number ",ifound)
if (ifound==1) then
    call skiplines(10,ikp)
    read(10,*) c200tmp,c200tmp,c200tmp,kp1crd,kp2crd,kp3crd
    write(*,"(' Coordinate of this k-point:',3f12.6)") kp1crd,kp2crd,kp3crd
else
    write(*,*) "Input fractional coordinate of this k-point in reciprocal space"
    write(*,*) "e.g. 0,0.5,0"
    write(*,*) "If press ENTER button directly, 0,0,0 will be used"
    read(*,"(a)") c80tmp
    if (c80tmp==" ") then
        kp1crd=0
        kp2crd=0
        kp3crd=0
    else
        read(c80tmp,*) kp1crd,kp2crd,kp3crd
    end if
end if

!Try to find number of orbitals, available when PRINT_LEVEL >= MEDIUM
call loclabel(10,"Number of molecular orbitals:",ifound)
if (ifound==1) then
    read(10,"(a)") c200tmp
    call readaftersign_int(c200tmp,':',nprintorb)
    write(*,"(' Number of orbitals to read:',i10)") nprintorb
else
    write(*,*)
    write(*,*) "How many orbitals were printed in the CP2K output file? e.g. 18"
    read(*,*) nprintorb
end if

call loclabelfinal(10,"FOR K POINT "//trim(adjustl(c80tmp)),ifound)
call skiplines(10,2)

!Load CP2K coefficient matrix to CObasa_tmp
write(*,*) "Loading..."
allocate(CObasa_tmp(nbasis,nbasis))
CObasa_tmp=0
neach=3 !Number of columns each frame. This is dynamic, when &MO / NDIGITS 8, it should be 3; while NDIGITS 6, it should be 4
nframe=ceiling(nprintorb/dfloat(neach))
do iframe=1,nframe
    ibeg=(iframe-1)*neach+1
    if (iframe==nframe) then
        iend=(iframe-1)*neach+mod(nprintorb,neach)
    else
        iend=(iframe-1)*neach+neach
    end if
    read(10,*)
    read(10,*) c80tmp,MOene(ibeg:iend)
    read(10,*)
    read(10,*) c80tmp,MOocc(ibeg:iend)
    read(10,*)
    do iatm=1,ncenter
        do ibas=basstart(iatm),basend(iatm)
            read(10,"(a)") c200tmp
            !write(*,*) trim(c200tmp)
            read(c200tmp(27:),*) CObasa_tmp(ibas,ibeg:iend)
        end do
        read(10,*)
    end do
end do
close(10)

!Reorder coefficients, map CObasa_tmp to CObasa
ibas=1
CObasa=0
do while(ibas<=nbasis)
    if (bastype(ibas)==1) then !S
        CObasa(ibas,:)=CObasa_tmp(ibas,:)
        ibas=ibas+1
    else if (bastype(ibas)==2) then !PX
        !CP2K: Y,Z,X
        CObasa(ibas,:)=CObasa_tmp(ibas+2,:)   !X
        CObasa(ibas+1,:)=CObasa_tmp(ibas,:)   !Y
        CObasa(ibas+2,:)=CObasa_tmp(ibas+1,:) !Z
        ibas=ibas+3
    else if (bastype(ibas)==-5) then !D0
        !CP2K: D-2,D-1,D0,D+1,D+2
        CObasa(ibas,:)=CObasa_tmp(ibas+2,:)
        CObasa(ibas+1,:)=CObasa_tmp(ibas+3,:)
        CObasa(ibas+2,:)=CObasa_tmp(ibas+1,:)
        CObasa(ibas+3,:)=CObasa_tmp(ibas+4,:)
        CObasa(ibas+4,:)=CObasa_tmp(ibas,:)
        ibas=ibas+5
    else if (bastype(ibas)==-12) then !F0
        !CP2K: F-3,F-2,F-1,F0,F+1,F+2,F+3
        CObasa(ibas,:)=CObasa_tmp(ibas+3,:)
        CObasa(ibas+1,:)=CObasa_tmp(ibas+4,:)
        CObasa(ibas+2,:)=CObasa_tmp(ibas+2,:)
        CObasa(ibas+3,:)=CObasa_tmp(ibas+5,:)
        CObasa(ibas+4,:)=CObasa_tmp(ibas+1,:)
        CObasa(ibas+5,:)=CObasa_tmp(ibas+6,:)
        CObasa(ibas+6,:)=CObasa_tmp(ibas,:)
        ibas=ibas+7
    else if (bastype(ibas)==-21) then !G0
        !CP2K: G-4,G-3,G-2,G-1,G0,G+1,G+2,G+3,G+4
        CObasa(ibas,:)=CObasa_tmp(ibas+4,:)   !G0
        CObasa(ibas+1,:)=CObasa_tmp(ibas+5,:) !G+1
        CObasa(ibas+2,:)=CObasa_tmp(ibas+3,:) !G-1
        CObasa(ibas+3,:)=CObasa_tmp(ibas+6,:) !G+2
        CObasa(ibas+4,:)=CObasa_tmp(ibas+2,:) !G-2
        CObasa(ibas+5,:)=CObasa_tmp(ibas+7,:) !G+3
        CObasa(ibas+6,:)=CObasa_tmp(ibas+1,:) !G-3
        CObasa(ibas+7,:)=CObasa_tmp(ibas+8,:) !G+4
        CObasa(ibas+8,:)=CObasa_tmp(ibas,:)   !G-4
        ibas=ibas+9
    else if (bastype(ibas)==-32) then !H0
        !CP2K: H-5,H-4,H-3,H-2,H-1,H0,H+1,H+2,H+3,H+4,H+5
        CObasa(ibas,:)=CObasa_tmp(ibas+5,:)    !H0
        CObasa(ibas+1,:)=CObasa_tmp(ibas+6,:)  !H+1
        CObasa(ibas+2,:)=CObasa_tmp(ibas+4,:)  !H-1
        CObasa(ibas+3,:)=CObasa_tmp(ibas+7,:)  !H+2
        CObasa(ibas+4,:)=CObasa_tmp(ibas+3,:)  !H-2
        CObasa(ibas+5,:)=CObasa_tmp(ibas+8,:)  !H+3
        CObasa(ibas+6,:)=CObasa_tmp(ibas+2,:)  !H-3
        CObasa(ibas+7,:)=CObasa_tmp(ibas+9,:)  !H+4
        CObasa(ibas+8,:)=CObasa_tmp(ibas+1,:)  !H-4
        CObasa(ibas+9,:)=CObasa_tmp(ibas+10,:) !H+5
        CObasa(ibas+10,:)=CObasa_tmp(ibas,:)   !H-5
        ibas=ibas+11
    end if
end do

write(*,*) "Loading finished!"

if (wfntype==0.or.wfntype==3) then
    call CObas2CO(1)
else
    call CObas2CO(3)
end if
call genP
end subroutine




!!--------- Load overlap matrix from .csr file exported by CP2K
!If infilepath is not " ", then directly load overlap matrix from it
subroutine CP2K_loadSbas(infilepath)
use defvar
implicit real*8 (a-h,o-z)
character c200tmp*200,infilepath*200

if (infilepath==" ") then
    write(*,"(a)") " Input the .csr file exported by CP2K containing overlap matrix of present system, e.g. D:\Palaio\Faliro.csr"
    write(*,*) "Note: The matrix should be in real space and only upper triangular part is recorded" 
    do while(.true.)
	    read(*,"(a)") c200tmp
	    inquire(file=c200tmp,exist=alive)
	    if (alive) exit
	    write(*,*) "Cannot find the file, input again!"
    end do
else
    c200tmp=infilepath
end if

if (allocated(Sbas)) deallocate(Sbas)
allocate(Sbas(nbasis,nbasis))
Sbas=0

write(*,*) "Loading..."
open(10,file=c200tmp,status="old")
do while(.true.) !Note that when CP2K outputting upper triangular part, very few elements (I think should be very small) are not printed, very strange
    read(10,*,iostat=ierror) ibas,jbas,Sbas(ibas,jbas)
    if (ierror/=0) exit
end do
do ibas=1,nbasis !Fill lower triangular part
    do jbas=ibas,nbasis
        Sbas(jbas,ibas)=Sbas(ibas,jbas)
    end do
end do
close(10)
write(*,*) "Overlap matrix has been successfully loaded!"

write(*,*) "Reordering matrix..."
call CP2K_mat_reorder(Sbas)
write(*,*) "Done!"
end subroutine




!!--------- Reordering basis function of loaded CP2K matrix to Multiwfn convention
subroutine CP2K_mat_reorder(mat)
use defvar
integer ibas
real*8 mat(nbasis,nbasis),tmpmat(nbasis,nbasis)

tmpmat=0
ibas=1
!Reorder rows
do while(ibas<=nbasis)
    if (bastype(ibas)==1) then !S
        tmpmat(ibas,:)=mat(ibas,:)
        ibas=ibas+1
    else if (bastype(ibas)==2) then !PX
        !CP2K: Y,Z,X
        tmpmat(ibas,:)=mat(ibas+2,:)   !X
        tmpmat(ibas+1,:)=mat(ibas,:)   !Y
        tmpmat(ibas+2,:)=mat(ibas+1,:) !Z
        ibas=ibas+3
    else if (bastype(ibas)==-5) then !D0
        !CP2K: D-2,D-1,D0,D+1,D+2
        tmpmat(ibas,:)=mat(ibas+2,:)
        tmpmat(ibas+1,:)=mat(ibas+3,:)
        tmpmat(ibas+2,:)=mat(ibas+1,:)
        tmpmat(ibas+3,:)=mat(ibas+4,:)
        tmpmat(ibas+4,:)=mat(ibas,:)
        ibas=ibas+5
    else if (bastype(ibas)==-12) then !F0
        !CP2K: F-3,F-2,F-1,F0,F+1,F+2,F+3
        tmpmat(ibas,:)=mat(ibas+3,:)
        tmpmat(ibas+1,:)=mat(ibas+4,:)
        tmpmat(ibas+2,:)=mat(ibas+2,:)
        tmpmat(ibas+3,:)=mat(ibas+5,:)
        tmpmat(ibas+4,:)=mat(ibas+1,:)
        tmpmat(ibas+5,:)=mat(ibas+6,:)
        tmpmat(ibas+6,:)=mat(ibas,:)
        ibas=ibas+7
    else if (bastype(ibas)==-21) then !G0
        !CP2K: G-4,G-3,G-2,G-1,G0,G+1,G+2,G+3,G+4
        tmpmat(ibas,:)=mat(ibas+4,:)   !G0
        tmpmat(ibas+1,:)=mat(ibas+5,:) !G+1
        tmpmat(ibas+2,:)=mat(ibas+3,:) !G-1
        tmpmat(ibas+3,:)=mat(ibas+6,:) !G+2
        tmpmat(ibas+4,:)=mat(ibas+2,:) !G-2
        tmpmat(ibas+5,:)=mat(ibas+7,:) !G+3
        tmpmat(ibas+6,:)=mat(ibas+1,:) !G-3
        tmpmat(ibas+7,:)=mat(ibas+8,:) !G+4
        tmpmat(ibas+8,:)=mat(ibas,:)   !G-4
        ibas=ibas+9
    else if (bastype(ibas)==-32) then !H0
        !CP2K: H-5,H-4,H-3,H-2,H-1,H0,H+1,H+2,H+3,H+4,H+5
        tmpmat(ibas,:)=mat(ibas+5,:)    !H0
        tmpmat(ibas+1,:)=mat(ibas+6,:)  !H+1
        tmpmat(ibas+2,:)=mat(ibas+4,:)  !H-1
        tmpmat(ibas+3,:)=mat(ibas+7,:)  !H+2
        tmpmat(ibas+4,:)=mat(ibas+3,:)  !H-2
        tmpmat(ibas+5,:)=mat(ibas+8,:)  !H+3
        tmpmat(ibas+6,:)=mat(ibas+2,:)  !H-3
        tmpmat(ibas+7,:)=mat(ibas+9,:)  !H+4
        tmpmat(ibas+8,:)=mat(ibas+1,:)  !H-4
        tmpmat(ibas+9,:)=mat(ibas+10,:) !H+5
        tmpmat(ibas+10,:)=mat(ibas,:)   !H-5
        ibas=ibas+11
    end if
end do
!Reorder columns
mat=0
ibas=1
do while(ibas<=nbasis)
    if (bastype(ibas)==1) then !S
        mat(:,ibas)=tmpmat(:,ibas)
        ibas=ibas+1
    else if (bastype(ibas)==2) then !PX
        !CP2K: Y,Z,X
        mat(:,ibas)=tmpmat(:,ibas+2)   !X
        mat(:,ibas+1)=tmpmat(:,ibas)   !Y
        mat(:,ibas+2)=tmpmat(:,ibas+1) !Z
        ibas=ibas+3
    else if (bastype(ibas)==-5) then !D0
        !CP2K: D-2,D-1,D0,D+1,D+2
        mat(:,ibas)=tmpmat(:,ibas+2)
        mat(:,ibas+1)=tmpmat(:,ibas+3)
        mat(:,ibas+2)=tmpmat(:,ibas+1)
        mat(:,ibas+3)=tmpmat(:,ibas+4)
        mat(:,ibas+4)=tmpmat(:,ibas)
        ibas=ibas+5
    else if (bastype(ibas)==-12) then !F0
        !CP2K: F-3,F-2,F-1,F0,F+1,F+2,F+3
        mat(:,ibas)=tmpmat(:,ibas+3)
        mat(:,ibas+1)=tmpmat(:,ibas+4)
        mat(:,ibas+2)=tmpmat(:,ibas+2)
        mat(:,ibas+3)=tmpmat(:,ibas+5)
        mat(:,ibas+4)=tmpmat(:,ibas+1)
        mat(:,ibas+5)=tmpmat(:,ibas+6)
        mat(:,ibas+6)=tmpmat(:,ibas)
        ibas=ibas+7
    else if (bastype(ibas)==-21) then !G0
        !CP2K: G-4,G-3,G-2,G-1,G0,G+1,G+2,G+3,G+4
        mat(:,ibas)=tmpmat(:,ibas+4)   !G0
        mat(:,ibas+1)=tmpmat(:,ibas+5) !G+1
        mat(:,ibas+2)=tmpmat(:,ibas+3) !G-1
        mat(:,ibas+3)=tmpmat(:,ibas+6) !G+2
        mat(:,ibas+4)=tmpmat(:,ibas+2) !G-2
        mat(:,ibas+5)=tmpmat(:,ibas+7) !G+3
        mat(:,ibas+6)=tmpmat(:,ibas+1) !G-3
        mat(:,ibas+7)=tmpmat(:,ibas+8) !G+4
        mat(:,ibas+8)=tmpmat(:,ibas)   !G-4
        ibas=ibas+9
    else if (bastype(ibas)==-32) then !H0
        !CP2K: H-5,H-4,H-3,H-2,H-1,H0,H+1,H+2,H+3,H+4,H+5
        mat(:,ibas)=tmpmat(:,ibas+5)    !H0
        mat(:,ibas+1)=tmpmat(:,ibas+6)  !H+1
        mat(:,ibas+2)=tmpmat(:,ibas+4)  !H-1
        mat(:,ibas+3)=tmpmat(:,ibas+7)  !H+2
        mat(:,ibas+4)=tmpmat(:,ibas+3)  !H-2
        mat(:,ibas+5)=tmpmat(:,ibas+8)  !H+3
        mat(:,ibas+6)=tmpmat(:,ibas+2)  !H-3
        mat(:,ibas+7)=tmpmat(:,ibas+9)  !H+4
        mat(:,ibas+8)=tmpmat(:,ibas+1)  !H-4
        mat(:,ibas+9)=tmpmat(:,ibas+10) !H+5
        mat(:,ibas+10)=tmpmat(:,ibas)   !H-5
        ibas=ibas+11
    end if
end do
end subroutine

