diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 57cd6e1a6730f333573af5ab7487f6e69650e623..927bff0c047243c65f4cc7cdce2e941d8e932729 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -29,7 +29,7 @@ test-gfortran: paths: - build script: - - ulimit -s unlimited ;export juDFT_MPI="mpirun -n 2 --allow-run-as-root ";cd /builds/fleur/fleur/build;ctest + - ulimit -s unlimited ;export juDFT_MPI="mpirun -n 2 --allow-run-as-root ";export OMP_NUM_THREADS=4;cd /builds/fleur/fleur/build;ctest artifacts: when: on_failure paths: diff --git a/eigen_soc/abcof_soc.F90 b/eigen_soc/abcof_soc.F90 index 00f94e15e80051d99d0e9f9b663c115bc6771ce0..876f3cafd3fdb0a726ebb5f36b3050cd0c87f92e 100644 --- a/eigen_soc/abcof_soc.F90 +++ b/eigen_soc/abcof_soc.F90 @@ -114,12 +114,12 @@ CONTAINS ENDIF !---> loop over lapws #ifndef CPP_OLDINTEL - !$OMP PARALLEL DO & - !$OMP& DEFAULT(none)& - !$OMP& PRIVATE(k,i,work_r,work_c,ccchi,kspin,fg,fk,s,r1,fj,dfj,l,df,wronk,tmk,phase,lo,nkvec,& - !$OMP& inap,nap,j,fgr,fgp,s2h,s2h_e,fkr,fkp,ylm,ll1,m,c_0,c_1,c_2,lmp,inv_f,lm)& - !$OMP& SHARED(n,nn,natom,natom_l,noco,atoms,sym,cell,oneD,lapw,nvmax,ne,zMat,usdus,ci,iintsp,eig,l_force,& - !$OMP& alo1,blo1,clo1,jatom,jspin,apw,const,nbasf0,acof,bcof,ccof,force,nat_start,nat_stop) + !!$OMP PARALLEL DO & + !!$OMP& DEFAULT(none)& + !!$OMP& PRIVATE(k,i,work_r,work_c,ccchi,kspin,fg,fk,s,r1,fj,dfj,l,df,wronk,tmk,phase,lo,nkvec,& + !!$OMP& inap,nap,j,fgr,fgp,s2h,s2h_e,fkr,fkp,ylm,ll1,m,c_0,c_1,c_2,lmp,inv_f,lm)& + !!$OMP& SHARED(n,nn,natom,natom_l,noco,atoms,sym,cell,oneD,lapw,nvmax,ne,zMat,usdus,ci,iintsp,eig,l_force,& + !!$OMP& alo1,blo1,clo1,jatom,jspin,apw,const,nbasf0,acof,bcof,ccof,force,nat_start,nat_stop) #endif DO k = 1,nvmax IF (zmat%l_real) THEN @@ -217,7 +217,7 @@ CONTAINS END DO ENDDO ! loop over LAPWs (k) #ifndef CPP_OLDINTEL - !$OMP END PARALLEL DO + !!$OMP END PARALLEL DO #endif IF (zmat%l_real) THEN DEALLOCATE(work_r) diff --git a/io/eig66_data.F90 b/io/eig66_data.F90 index 2c4cc11a854b6cc482cbed0c961ce410f23f5d85..935e887a0952b38490608067a9cabab853301864 100644 --- a/io/eig66_data.F90 +++ b/io/eig66_data.F90 @@ -30,9 +30,9 @@ module m_eig66_data INTEGER,ALLOCATABLE :: pe_basis(:,:),slot_basis(:,:) INTEGER,ALLOCATABLE :: pe_ev(:,:,:),slot_ev(:,:,:) INTEGER :: irank - INTEGER,POINTER :: neig_data(:) - REAL,POINTER :: eig_data(:),zr_data(:), w_iks_data(:) - COMPLEX,POINTER :: zc_data(:) + INTEGER,POINTER :: neig_data(:) + REAL,POINTER :: eig_data(:),zr_data(:), w_iks_data(:) + COMPLEX,POINTER :: zc_data(:) END TYPE TYPE,EXTENDS(t_data):: t_data_hdf #ifdef CPP_HDF @@ -82,7 +82,7 @@ module m_eig66_data INTEGER,INTENT(IN),OPTIONAL :: io_mode CLASS(t_data),pointer::d - TYPE(t_list),POINTER:: listpointer,lastinlist + TYPE(t_list),POINTER,ASYNCHRONOUS:: listpointer,lastinlist lastinlist=>null() listpointer=>linked_list diff --git a/io/eig66_mpi.F90 b/io/eig66_mpi.F90 index 55668317f6ac6e3392ea2f552faa5b0d32e972e3..c35298a4225551da2940e901e71fb3503e8b8361 100644 --- a/io/eig66_mpi.F90 +++ b/io/eig66_mpi.F90 @@ -2,6 +2,7 @@ MODULE m_eig66_mpi #include "juDFT_env.h" USE m_eig66_data USE m_types + USE m_judft #ifdef CPP_MPI use mpi #endif @@ -12,7 +13,7 @@ CONTAINS SUBROUTINE priv_find_data(id,d) INTEGER,INTENT(IN)::id - TYPE(t_data_mpi),POINTER:: d + TYPE(t_data_mpi),POINTER,ASYNCHRONOUS:: d CLASS(t_data),POINTER ::dp CALL eig66_find_data(dp,id) @@ -35,24 +36,15 @@ CONTAINS #ifdef CPP_MPI INTEGER:: isize,e,slot_size,local_slots INTEGER,PARAMETER::mcored=27 !there should not be more that 27 core states - TYPE(t_data_MPI),POINTER :: d + TYPE(t_data_MPI),POINTER,ASYNCHRONOUS :: d CALL priv_find_data(id,d) CALL eig66_data_storedefault(d,jspins,nkpts,nmat,neig,l_real.and..not.l_soc,l_soc) IF (PRESENT(n_size_opt)) d%n_size=n_size_opt IF (ALLOCATED(d%pe_ev)) THEN - IF (create) THEN - d%neig_data=0 - d%eig_data=1E99 - d%w_iks_data=1E99 - if (d%l_real.and..not.l_soc) THEN - d%zr_data=0.0 - else - d%zc_data=0.0 - endif - ENDIF - IF (PRESENT(filename)) CALL priv_readfromfile() + IF (create) CALL reset_eig(id,l_soc) + IF (PRESENT(filename)) CALL judft_error("Storing of data not implemented for MPI case",calledby="eig66_mpi.F") RETURN !everything already done! ENDIF @@ -81,21 +73,21 @@ CONTAINS local_slots=COUNT(d%pe_ev==d%irank) slot_size=nmat - if (l_real.and..not.l_soc) THEN + IF (l_real.AND..NOT.l_soc) THEN CALL priv_create_memory(slot_size,local_slots,d%zr_handle,real_data_ptr=d%zr_data) else CALL priv_create_memory(slot_size,local_slots,d%zc_handle,cmplx_data_ptr=d%zc_data) - endif - IF (PRESENT(filename).AND..NOT.create) CALL priv_readfromfile() + ENDIF + IF (PRESENT(filename).AND..NOT.create) CALL judft_error("Storing of data not implemented for MPI case",calledby="eig66_mpi.F") CALL MPI_BARRIER(MPI_COMM,e) CALL timestop("create data spaces in ei66_mpi") CONTAINS SUBROUTINE priv_create_memory(slot_size,local_slots,handle,int_data_ptr,real_data_ptr,cmplx_data_ptr) IMPLICIT NONE INTEGER,INTENT(IN) :: slot_size,local_slots - INTEGER,POINTER,OPTIONAL :: int_data_ptr(:) - REAL ,POINTER,OPTIONAL :: real_data_ptr(:) - COMPLEX,POINTER,OPTIONAL :: cmplx_data_ptr(:) + INTEGER,POINTER,OPTIONAL,ASYNCHRONOUS :: int_data_ptr(:) + REAL ,POINTER,OPTIONAL,ASYNCHRONOUS :: real_data_ptr(:) + COMPLEX,POINTER,OPTIONAL,ASYNCHRONOUS :: cmplx_data_ptr(:) INTEGER,INTENT(OUT) :: handle #ifdef CPP_MPI TYPE(c_ptr)::ptr @@ -119,91 +111,50 @@ CONTAINS if (length.ne.1) call judft_error("Bug in eig66_mpi:create_memory") length=MAX(1,slot_size*local_slots) +#ifdef CPP_MPI_ALLOC length=length*type_size - CALL MPI_ALLOC_MEM(length,MPI_INFO_NULL,ptr,e) IF (e.NE.0) CPP_error("Could not allocated MPI-Data in eig66_mpi") - - IF (present(real_data_ptr)) THEN - CALL C_F_POINTER(ptr,real_data_ptr,(/length/type_size/)) +#endif + IF (PRESENT(real_data_ptr)) THEN +#ifdef CPP_MPI_ALLOC + CALL C_F_POINTER(ptr,real_data_ptr,(/length/type_size/)) +#else + ALLOCATE(real_data_ptr(length)) +#endif CALL MPI_WIN_CREATE(real_data_ptr, length,slot_size*type_size,Mpi_INFO_NULL, MPI_COMM,handle, e) - ELSEIF(present(int_data_ptr)) THEN - CALL C_F_POINTER(ptr,int_data_ptr,(/length/type_size/)) + ELSEIF(PRESENT(int_data_ptr)) THEN +#ifdef CPP_MPI_ALLOC + CALL C_F_POINTER(ptr,int_data_ptr,(/length/type_size/)) +#else + ALLOCATE(int_data_ptr(length)) +#endif CALL MPI_WIN_CREATE(int_data_ptr, length,slot_size*type_size,Mpi_INFO_NULL, MPI_COMM,handle, e) - ELSE - CALL C_F_POINTER(ptr,cmplx_data_ptr,(/length/type_size/)) - CALL MPI_WIN_CREATE(cmplx_data_ptr, length,slot_size*type_size,Mpi_INFO_NULL, MPI_COMM,handle, e) - ENDIF + ELSE +#ifdef CPP_MPI_ALLOC + CALL C_F_POINTER(ptr,cmplx_data_ptr,(/length/type_size/)) +#else + ALLOCATE(cmplx_data_ptr(length)) +#endif + CALL MPI_WIN_CREATE(cmplx_data_ptr, length,slot_size*type_size,Mpi_INFO_NULL, MPI_COMM,handle, e) + ENDIF #endif END SUBROUTINE priv_create_memory - SUBROUTINE priv_readfromfile() - USE m_eig66_DA,ONLY:open_eig_DA=>open_eig,read_eig_DA=>read_eig,close_eig_da=>close_eig - INTEGER:: jspin,nk,i,ii,iii,nv,tmp_id - REAL :: wk,bk3(3),evac(2) - REAL :: eig(neig),w_iks(neig) - TYPE(t_zmat)::zmat - zmat%l_real=d%l_real - zmat%nbasfcn=nmat - zmat%nbands=neig - allocate(zmat%z_r(nmat,neig),zmat%z_c(nmat,neig)) - !only do this with PE=0 - IF (d%irank==0) THEN - tmp_id=eig66_data_newid(DA_mode) - CALL open_eig_DA(tmp_id,nmat,neig,nkpts,jspins,.FALSE.,d%l_real,l_soc,filename) - DO jspin=1,jspins - DO nk=1,nkpts - !CALL read_eig_DA(tmp_id,nk,jspin,nv,i,bk3,wk,ii,eig,w_iks,el,ello,evac,zmat=zmat) - STOP "code no longer works" - ! CALL write_eig(id,nk,jspin,ii,ii,nv,nmat,bk3,wk,eig,w_iks,el,ello,evac,nlotot,zmat=zmat) - ENDDO - ENDDO - CALL close_eig_DA(tmp_id) - ENDIF - END SUBROUTINE priv_readfromfile #endif END SUBROUTINE open_eig SUBROUTINE close_eig(id,delete,filename) INTEGER,INTENT(IN) :: id LOGICAL,INTENT(IN),OPTIONAL:: delete CHARACTER(LEN=*),INTENT(IN),OPTIONAL::filename - TYPE(t_data_MPI),POINTER :: d + TYPE(t_data_MPI),POINTER,ASYNCHRONOUS :: d CALL priv_find_data(id,d) IF (PRESENT(delete)) THEN IF (delete) WRITE(*,*) "No deallocation of memory implemented in eig66_mpi" ENDIF - IF (PRESENT(filename)) CALL priv_writetofile() - CONTAINS - SUBROUTINE priv_writetofile() - USE m_eig66_DA,ONLY:open_eig_DA=>open_eig,write_eig_DA=>write_eig,close_eig_DA=>close_eig - IMPLICIT NONE - - INTEGER:: nk,jspin,nv,i,ii,tmp_id - REAL :: wk,bk3(3),evac(2) - REAL :: eig(d%neig),w_iks(d%neig) - TYPE(t_zmat)::zmat - zmat%l_real=d%l_real - zmat%nbasfcn=d%nmat - zmat%nbands=d%neig - allocate(zmat%z_r(d%nmat,d%neig),zmat%z_c(d%nmat,d%neig)) - - IF (d%irank==0) THEN - tmp_id=eig66_data_newid(DA_mode) - CALL open_eig_DA(tmp_id,d%nmat,d%neig,d%nkpts,d%jspins,.FALSE.,d%l_real,d%l_soc,filename) - DO jspin=1,d%jspins - DO nk=1,d%nkpts - !CALL read_eig(id,nk,jspin,nv,i,bk3,wk,ii,eig,w_iks,el,ello,evac,zmat=zmat) - stop "CODE no longer working" - !CALL write_eig_DA(tmp_id,nk,jspin,ii,ii,nv,i,bk3,wk,eig,w_iks,el,ello,evac,nlotot,zmat=zmat) - ENDDO - ENDDO - CALL close_eig_DA(tmp_id) - ENDIF - CALL eig66_remove_data(id) - END SUBROUTINE priv_writetofile - + IF (PRESENT(filename)) CALL judft_error("Storing of data not implemented for MPI case",calledby="eig66_mpi.F") END SUBROUTINE close_eig SUBROUTINE read_eig(id,nk,jspin,neig,eig,w_iks,n_start,n_end,zmat) @@ -221,7 +172,7 @@ CONTAINS INTEGER,ALLOCATABLE :: tmp_int(:) REAL,ALLOCATABLE :: tmp_real(:) COMPLEX,ALLOCATABLE :: tmp_cmplx(:) - TYPE(t_data_MPI),POINTER :: d + TYPE(t_data_MPI),POINTER,ASYNCHRONOUS :: d CALL priv_find_data(id,d) pe=d%pe_basis(nk,jspin) slot=d%slot_basis(nk,jspin) @@ -308,7 +259,7 @@ CONTAINS REAL,ALLOCATABLE :: tmp_real(:) COMPLEX,ALLOCATABLE :: tmp_cmplx(:) LOGICAL :: acc - TYPE(t_data_MPI),POINTER :: d + TYPE(t_data_MPI),POINTER,ASYNCHRONOUS :: d CALL priv_find_data(id,d) @@ -394,7 +345,7 @@ CONTAINS INTEGER, INTENT(IN) :: id LOGICAL, INTENT(IN) :: l_soc #ifdef CPP_MPI - TYPE(t_data_MPI),POINTER :: d + TYPE(t_data_MPI),POINTER,ASYNCHRONOUS :: d CALL priv_find_data(id,d) d%neig_data=0 @@ -408,84 +359,11 @@ CONTAINS #endif END SUBROUTINE reset_eig -#ifdef CPP_MPI - SUBROUTINE priv_put_data(pe,slot,DATA,handle) - IMPLICIT NONE - INTEGER,INTENT(IN) :: pe,slot - CLASS(*),INTENT(IN) :: DATA(:) - INTEGER,INTENT(IN) :: handle - - INTEGER :: len,e - INTEGER,ALLOCATABLE :: int_tmp(:) - REAL,ALLOCATABLE :: real_tmp(:) - COMPLEX,ALLOCATABLE:: cmplx_tmp(:) - INCLUDE 'mpif.h' - len=SIZE(DATA) - SELECT TYPE(DATA) - TYPE IS (INTEGER) - ALLOCATE(int_tmp(len)) - int_tmp=DATA - CALL MPI_WIN_LOCK(MPI_LOCK_EXCLUSIVE,pe,0,handle,e) - CALL MPI_PUT(int_tmp,len,MPI_INTEGER,pe,int(slot,MPI_ADDRESS_KIND),len,MPI_INTEGER,handle,e) - CALL MPI_WIN_UNLOCK(pe,handle,e) - TYPE is (REAL) - ALLOCATE(real_tmp(len)) - real_tmp=DATA - CALL MPI_WIN_LOCK(MPI_LOCK_EXCLUSIVE,pe,0,handle,e) - CALL MPI_PUT(real_tmp,len,MPI_DOUBLE_PRECISION,pe,int(slot,MPI_ADDRESS_KIND),len,MPI_DOUBLE_PRECISION,handle,e) - CALL MPI_WIN_UNLOCK(pe,handle,e) - TYPE is (COMPLEX) - ALLOCATE(cmplx_tmp(len)) - cmplx_tmp=DATA - CALL MPI_WIN_LOCK(MPI_LOCK_EXCLUSIVE,pe,0,handle,e) - CALL MPI_PUT(cmplx_tmp,len,MPI_DOUBLE_COMPLEX,pe,int(slot,MPI_ADDRESS_KIND),len,MPI_DOUBLE_COMPLEX,handle,e) - CALL MPI_WIN_UNLOCK(pe,handle,e) - END SELECT - END SUBROUTINE priv_put_data - - SUBROUTINE priv_get_data(pe,slot,len,handle,idata,rdata,cdata) - IMPLICIT NONE - INTEGER,INTENT(IN) :: pe,slot,len - INTEGER,INTENT(OUT),optional :: iDATA(len) - REAL,INTENT(OUT),optional :: rDATA(len) - COMPLEX,INTENT(OUT),optional :: cDATA(len) - INTEGER,INTENT(IN) :: handle - - INTEGER :: e - INTEGER,ALLOCATABLE :: int_tmp(:) - REAL,ALLOCATABLE :: real_tmp(:) - COMPLEX,ALLOCATABLE:: cmplx_tmp(:) - INCLUDE 'mpif.h' - - IF (present(idata)) THEN - ALLOCATE(int_tmp(len)) - CALL MPI_WIN_LOCK(MPI_LOCK_SHARED,pe,0,handle,e) - CALL MPI_GET(int_tmp,len,MPI_INTEGER,pe,int(slot,MPI_ADDRESS_KIND),len,MPI_INTEGER,handle,e) - CALL MPI_WIN_UNLOCK(pe,handle,e) - iDATA=int_tmp - ELSE IF (PRESENT(rdata)) THEN - ALLOCATE(real_tmp(len)) - CALL MPI_WIN_LOCK(MPI_LOCK_SHARED,pe,0,handle,e) - CALL MPI_GET(real_tmp,len,MPI_DOUBLE_PRECISION,pe,int(slot,MPI_ADDRESS_KIND),len,MPI_DOUBLE_PRECISION,handle,e) - CALL MPI_WIN_UNLOCK(pe,handle,e) - rDATA=real_tmp - ELSE IF (PRESENT(cdata)) THEN - ALLOCATE(cmplx_tmp(len)) - CALL MPI_WIN_LOCK(MPI_LOCK_SHARED,pe,0,handle,e) - CALL MPI_GET(cmplx_tmp,len,MPI_DOUBLE_COMPLEX,pe,int(slot,MPI_ADDRESS_KIND),len,MPI_DOUBLE_COMPLEX,handle,e) - CALL MPI_WIN_UNLOCK(pe,handle,e) - cDATA=cmplx_tmp - ELSE - call judft_error("BUG in priv_get_data") - ENDIF - - END SUBROUTINE priv_get_data -#endif #ifdef CPP_MPI SUBROUTINE create_maps(d,isize,nkpts,jspins,neig,n_size) IMPLICIT NONE - TYPE(t_data_MPI),INTENT(INOUT):: d + TYPE(t_data_MPI),INTENT(INOUT),ASYNCHRONOUS:: d INTEGER,INTENT(IN):: isize,nkpts,jspins,neig,n_size INTEGER:: nk,j,n1,n2,n,pe,n_members diff --git a/io/r_inpXML.F90 b/io/r_inpXML.F90 index b1f496977f8f19152f9bf7dc01af43c09a96c994..fde4c4c4f347fe070ff2c9deb65a7925fd68a581 100644 --- a/io/r_inpXML.F90 +++ b/io/r_inpXML.F90 @@ -13,1095 +13,1088 @@ MODULE m_rinpXML !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CONTAINS -SUBROUTINE r_inpXML(& - atoms,obsolete,vacuum,input,stars,sliceplot,banddos,DIMENSION,forcetheo,& - cell,sym,xcpot,noco,oneD,hybrid,kpts,enpara,coreSpecInput,wann,& - noel,namex,relcor,a1,a2,a3,dtild,xmlElectronStates,& - xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType,& - l_kpts) - - USE iso_c_binding - USE m_juDFT - USE m_types - USE m_types_forcetheo_extended - USE m_symdata , ONLY : nammap, ord2, l_c2 - USE m_rwsymfile - USE m_xmlIntWrapFort - USE m_inv3 - USE m_spg2set - USE m_closure, ONLY : check_close - USE m_symproperties - USE m_calculator - USE m_constants - USE m_inpeig - USE m_sort - USE m_types_xcpot_inbuild -#ifdef CPP_LIBXC - USE xc_f03_lib_m -#endif - IMPLICIT NONE - - TYPE(t_input),INTENT(INOUT) :: input - TYPE(t_sym),INTENT(INOUT) :: sym - TYPE(t_stars),INTENT(INOUT) :: stars - TYPE(t_atoms),INTENT(INOUT) :: atoms - TYPE(t_vacuum),INTENT(INOUT) :: vacuum - TYPE(t_obsolete),INTENT(INOUT) :: obsolete - TYPE(t_kpts),INTENT(INOUT) :: kpts - TYPE(t_oneD),INTENT(INOUT) :: oneD - TYPE(t_hybrid),INTENT(INOUT) :: hybrid - TYPE(t_cell),INTENT(INOUT) :: cell - TYPE(t_banddos),INTENT(INOUT) :: banddos - TYPE(t_sliceplot),INTENT(INOUT):: sliceplot - CLASS(t_xcpot),INTENT(INOUT),ALLOCATABLE :: xcpot - TYPE(t_noco),INTENT(INOUT) :: noco - TYPE(t_dimension),INTENT(OUT) :: dimension - TYPE(t_enpara) ,INTENT(OUT) :: enpara - CLASS(t_forcetheo),ALLOCATABLE,INTENT(OUT):: forcetheo - TYPE(t_coreSpecInput),INTENT(OUT) :: coreSpecInput - TYPE(t_wann) ,INTENT(INOUT) :: wann - LOGICAL, INTENT(OUT) :: l_kpts - INTEGER, ALLOCATABLE, INTENT(INOUT) :: xmlElectronStates(:,:) - INTEGER, ALLOCATABLE, INTENT(INOUT) :: atomTypeSpecies(:) - INTEGER, ALLOCATABLE, INTENT(INOUT) :: speciesRepAtomType(:) - REAL, ALLOCATABLE, INTENT(INOUT) :: xmlCoreOccs(:,:,:) - LOGICAL, ALLOCATABLE, INTENT(INOUT) :: xmlPrintCoreStates(:,:) - CHARACTER(len=3), ALLOCATABLE, INTENT(INOUT) :: noel(:) - CHARACTER(len=4), INTENT(OUT) :: namex - CHARACTER(len=12), INTENT(OUT) :: relcor - REAL, INTENT(OUT) :: a1(3),a2(3),a3(3) - REAL, INTENT(OUT) :: dtild - - CHARACTER(len=8) :: name(10) - - !+odim - INTEGER MM,vM,m_cyl - LOGICAL invs1,zrfs1 - INTEGER chi,rot - LOGICAL d1,band - NAMELIST /odim/ d1,MM,vM,m_cyl,chi,rot,invs1,zrfs1 - !-odim - ! .. - ! .. Local Variables - REAL :: scpos ,zc - INTEGER ieq,i,k,na,n,ii,id_c,id_x - REAL s3,ah,a,hs2,rest,thetaj - LOGICAL l_hyb,l_sym,ldum - INTEGER :: ierr - ! .. - !... Local Arrays - ! CHARACTER :: helpchar(atoms%ntype) - CHARACTER(len= 4) :: chntype - CHARACTER(len= 41) :: chform - CHARACTER(len=100) :: line - - CHARACTER(len=20) :: tempNumberString - CHARACTER(len=150) :: format - CHARACTER(len=20) :: mixingScheme - CHARACTER(len=10) :: loType - LOGICAL :: kptGamma, l_relcor,ldummy - INTEGER :: iAtomType, startCoreStates, endCoreStates - CHARACTER(len=100) :: xPosString, yPosString, zPosString - CHARACTER(len=200) :: coreStatesString - ! REAL :: tempTaual(3,atoms%nat) - REAL :: coreStateOccs(29,2) - INTEGER :: coreStateNprnc(29), coreStateKappa(29) - INTEGER :: speciesXMLElectronStates(29) - REAL :: speciesXMLCoreOccs(2,29) - LOGICAL :: speciesXMLPrintCoreStates(29) - - INTEGER :: iType, iLO, iSpecies, lNumCount, nNumCount, iLLO, jsp, j, l, absSum, numTokens - INTEGER :: numberNodes, nodeSum, numSpecies, n2spg, n1, n2, ikpt, iqpt - INTEGER :: atomicNumber, coreStates, gridPoints, lmax, lnonsphr, lmaxAPW - INTEGER :: latticeDef, symmetryDef, nop48, firstAtomOfType, errorStatus - INTEGER :: loEDeriv, ntp1, ios, ntst, jrc, minNeigd, providedCoreStates, providedStates - INTEGER :: nv, nv2, kq1, kq2, kq3, nprncTemp, kappaTemp, tempInt - INTEGER :: ldau_l(4), numVac, numU - INTEGER :: speciesEParams(0:3) - INTEGER :: mrotTemp(3,3,48) - REAL :: tauTemp(3,48) - REAL :: bk(3) - LOGICAL :: flipSpin, l_eV, invSym, l_qfix, relaxX, relaxY, relaxZ - LOGICAL :: l_vca, coreConfigPresent, l_enpara, l_orbcomp, tempBool - REAL :: magMom, radius, logIncrement, qsc(3), latticeScale, dr - REAL :: aTemp, zp, rmtmax, sumWeight, ldau_u(4), ldau_j(4), tempReal - REAL :: weightScale, eParamUp, eParamDown - LOGICAL :: l_amf(4) - REAL, PARAMETER :: boltzmannConst = 3.1668114e-6 ! value is given in Hartree/Kelvin - INTEGER :: lcutm,lcutwf,hybSelect(4) - REAL :: evac0Temp(2,2) - - - CHARACTER(LEN=200,KIND=c_char) :: schemaFilename, docFilename - CHARACTER(LEN=255) :: valueString, lString, nString, token - CHARACTER(LEN=255) :: xPathA, xPathB, xPathC, xPathD, xPathE - CHARACTER(LEN=11) :: latticeType - CHARACTER(LEN=50) :: versionString - - LOGICAL :: ldaSpecies - REAL :: socscaleSpecies - - INTEGER, ALLOCATABLE :: lNumbers(:), nNumbers(:), speciesLLO(:) - INTEGER, ALLOCATABLE :: loOrderList(:) - INTEGER, ALLOCATABLE :: speciesNLO(:) - INTEGER, ALLOCATABLE :: multtab(:,:), invOps(:), optype(:) - INTEGER, ALLOCATABLE :: lmx1(:), nq1(:), nlhtp1(:) - INTEGER, ALLOCATABLE :: speciesLOEDeriv(:) - REAL, ALLOCATABLE :: speciesLOeParams(:), speciesLLOReal(:) - LOGICAL, ALLOCATABLE :: wannAtomList(:) + SUBROUTINE r_inpXML(& + atoms,obsolete,vacuum,input,stars,sliceplot,banddos,DIMENSION,forcetheo,& + cell,sym,xcpot,noco,oneD,hybrid,kpts,enpara,coreSpecInput,wann,& + noel,namex,relcor,a1,a2,a3,dtild,xmlElectronStates,& + xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType,& + l_kpts) + + USE iso_c_binding + USE m_juDFT + USE m_types + USE m_types_forcetheo_extended + USE m_symdata , ONLY : nammap, ord2, l_c2 + USE m_rwsymfile + USE m_xmlIntWrapFort + USE m_inv3 + USE m_spg2set + USE m_closure, ONLY : check_close + USE m_symproperties + USE m_calculator + USE m_constants + USE m_inpeig + USE m_sort + USE m_types_xcpot_inbuild +#ifdef CPP_LIBXC + USE xc_f03_lib_m +#endif + IMPLICIT NONE + + TYPE(t_input),INTENT(INOUT) :: input + TYPE(t_sym),INTENT(INOUT) :: sym + TYPE(t_stars),INTENT(INOUT) :: stars + TYPE(t_atoms),INTENT(INOUT) :: atoms + TYPE(t_vacuum),INTENT(INOUT) :: vacuum + TYPE(t_obsolete),INTENT(INOUT) :: obsolete + TYPE(t_kpts),INTENT(INOUT) :: kpts + TYPE(t_oneD),INTENT(INOUT) :: oneD + TYPE(t_hybrid),INTENT(INOUT) :: hybrid + TYPE(t_cell),INTENT(INOUT) :: cell + TYPE(t_banddos),INTENT(INOUT) :: banddos + TYPE(t_sliceplot),INTENT(INOUT):: sliceplot + CLASS(t_xcpot),INTENT(INOUT),ALLOCATABLE :: xcpot + TYPE(t_noco),INTENT(INOUT) :: noco + TYPE(t_dimension),INTENT(OUT) :: dimension + TYPE(t_enpara) ,INTENT(OUT) :: enpara + CLASS(t_forcetheo),ALLOCATABLE,INTENT(OUT):: forcetheo + TYPE(t_coreSpecInput),INTENT(OUT) :: coreSpecInput + TYPE(t_wann) ,INTENT(INOUT) :: wann + LOGICAL, INTENT(OUT) :: l_kpts + INTEGER, ALLOCATABLE, INTENT(INOUT) :: xmlElectronStates(:,:) + INTEGER, ALLOCATABLE, INTENT(INOUT) :: atomTypeSpecies(:) + INTEGER, ALLOCATABLE, INTENT(INOUT) :: speciesRepAtomType(:) + REAL, ALLOCATABLE, INTENT(INOUT) :: xmlCoreOccs(:,:,:) + LOGICAL, ALLOCATABLE, INTENT(INOUT) :: xmlPrintCoreStates(:,:) + CHARACTER(len=3), ALLOCATABLE, INTENT(INOUT) :: noel(:) + CHARACTER(len=4), INTENT(OUT) :: namex + CHARACTER(len=12), INTENT(OUT) :: relcor + REAL, INTENT(OUT) :: a1(3),a2(3),a3(3) + REAL, INTENT(OUT) :: dtild + + CHARACTER(len=8) :: name(10) + + !+odim + INTEGER MM,vM,m_cyl + LOGICAL invs1,zrfs1 + INTEGER chi,rot + LOGICAL d1,band + NAMELIST /odim/ d1,MM,vM,m_cyl,chi,rot,invs1,zrfs1 + !-odim + ! .. + ! .. Local Variables + REAL :: scpos ,zc + INTEGER ieq,i,k,na,n,ii,id_c,id_x + REAL s3,ah,a,hs2,rest,thetaj + LOGICAL l_hyb,l_sym,ldum + INTEGER :: ierr + ! .. + !... Local Arrays + ! CHARACTER :: helpchar(atoms%ntype) + CHARACTER(len= 4) :: chntype + CHARACTER(len= 41) :: chform + CHARACTER(len=100) :: line + + CHARACTER(len=20) :: tempNumberString + CHARACTER(len=150) :: format + CHARACTER(len=20) :: mixingScheme + CHARACTER(len=10) :: loType + LOGICAL :: kptGamma, l_relcor,ldummy + INTEGER :: iAtomType, startCoreStates, endCoreStates + CHARACTER(len=100) :: xPosString, yPosString, zPosString + CHARACTER(len=200) :: coreStatesString + ! REAL :: tempTaual(3,atoms%nat) + REAL :: coreStateOccs(29,2) + INTEGER :: coreStateNprnc(29), coreStateKappa(29) + INTEGER :: speciesXMLElectronStates(29) + REAL :: speciesXMLCoreOccs(2,29) + LOGICAL :: speciesXMLPrintCoreStates(29) + + INTEGER :: iType, iLO, iSpecies, lNumCount, nNumCount, iLLO, jsp, j, l, absSum, numTokens + INTEGER :: numberNodes, nodeSum, numSpecies, n2spg, n1, n2, ikpt, iqpt + INTEGER :: atomicNumber, coreStates, gridPoints, lmax, lnonsphr, lmaxAPW + INTEGER :: latticeDef, symmetryDef, nop48, firstAtomOfType, errorStatus + INTEGER :: loEDeriv, ntp1, ios, ntst, jrc, minNeigd, providedCoreStates, providedStates + INTEGER :: nv, nv2, kq1, kq2, kq3, nprncTemp, kappaTemp, tempInt + INTEGER :: ldau_l(4), numVac, numU + INTEGER :: speciesEParams(0:3) + INTEGER :: mrotTemp(3,3,48) + REAL :: tauTemp(3,48) + REAL :: bk(3) + LOGICAL :: flipSpin, l_eV, invSym, l_qfix, relaxX, relaxY, relaxZ + LOGICAL :: l_vca, coreConfigPresent, l_enpara, l_orbcomp, tempBool + REAL :: magMom, radius, logIncrement, qsc(3), latticeScale, dr + REAL :: aTemp, zp, rmtmax, sumWeight, ldau_u(4), ldau_j(4), tempReal + REAL :: weightScale, eParamUp, eParamDown + LOGICAL :: l_amf(4) + REAL, PARAMETER :: boltzmannConst = 3.1668114e-6 ! value is given in Hartree/Kelvin + INTEGER :: lcutm,lcutwf,hybSelect(4) + REAL :: evac0Temp(2,2) + + CHARACTER(LEN=200,KIND=c_char) :: schemaFilename, docFilename + CHARACTER(LEN=255) :: valueString, lString, nString, token + CHARACTER(LEN=255) :: xPathA, xPathB, xPathC, xPathD, xPathE + CHARACTER(LEN=11) :: latticeType + CHARACTER(LEN=50) :: versionString + + LOGICAL :: ldaSpecies + REAL :: socscaleSpecies + + INTEGER, ALLOCATABLE :: lNumbers(:), nNumbers(:), speciesLLO(:) + INTEGER, ALLOCATABLE :: loOrderList(:) + INTEGER, ALLOCATABLE :: speciesNLO(:) + INTEGER, ALLOCATABLE :: multtab(:,:), invOps(:), optype(:) + INTEGER, ALLOCATABLE :: lmx1(:), nq1(:), nlhtp1(:) + INTEGER, ALLOCATABLE :: speciesLOEDeriv(:) + REAL, ALLOCATABLE :: speciesLOeParams(:), speciesLLOReal(:) + LOGICAL, ALLOCATABLE :: wannAtomList(:) ! Variables for MT radius testing: - REAL :: dtild1,kmax1,dvac1 - LOGICAL :: l_test - INTEGER, ALLOCATABLE :: jri1(:), lmax1(:) - REAL, ALLOCATABLE :: rmt1(:), dx1(:) - - EXTERNAL prp_xcfft_box - - interface - function dropInputSchema() bind(C, name="dropInputSchema") - use iso_c_binding - INTEGER(c_int) dropInputSchema - end function dropInputSchema - end interface - - errorStatus = 0 - errorStatus = dropInputSchema() - IF(errorStatus.NE.0) THEN - STOP 'Error: Cannot print out FleurInputSchema.xsd' - END IF - - schemaFilename = "FleurInputSchema.xsd"//C_NULL_CHAR - docFilename = "inp.xml"//C_NULL_CHAR - - !TODO! these switches should be in the inp-file - input%l_core_confpot=.true. !former CPP_CORE - input%l_useapw=.false. !former CPP_APW - !WRITE(*,*) 'Start reading of inp.xml file' - CALL xmlInitInterface() - CALL xmlParseSchema(schemaFilename) - CALL xmlParseDoc(docFilename) - CALL xmlValidateDoc() - CALL xmlInitXPath() - - ! Check version of inp.xml - versionString = xmlGetAttributeValue('/fleurInput/@fleurInputVersion') - IF((TRIM(ADJUSTL(versionString)).NE.'0.27').AND.(TRIM(ADJUSTL(versionString)).NE.'0.28').AND.& - (TRIM(ADJUSTL(versionString)).NE.'0.29')) THEN - STOP 'version number of inp.xml file is not compatible with this fleur version' - END IF - - ! Get number of atoms, atom types, and atom species - - numberNodes = xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/relPos') - numberNodes = numberNodes + xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/absPos') - numberNodes = numberNodes + xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/filmPos') - - atoms%nat = numberNodes - - numberNodes = xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup') - - atoms%ntype = numberNodes - - numSpecies = xmlGetNumberOfNodes('/fleurInput/atomSpecies/species') - - ALLOCATE(atoms%nz(atoms%ntype)) !nz and zatom have the same content! - ALLOCATE(atoms%zatom(atoms%ntype)) !nz and zatom have the same content! - ALLOCATE(atoms%jri(atoms%ntype)) - ALLOCATE(atoms%dx(atoms%ntype)) - ALLOCATE(atoms%lmax(atoms%ntype)) - ALLOCATE(atoms%nlo(atoms%ntype)) - ALLOCATE(atoms%ncst(atoms%ntype)) - ALLOCATE(atoms%lnonsph(atoms%ntype)) - ALLOCATE(atoms%nflip(atoms%ntype)) - ALLOCATE(atoms%l_geo(atoms%ntype)) - ALLOCATE(atoms%lda_u(4*atoms%ntype)) - ALLOCATE(atoms%bmu(atoms%ntype)) - ALLOCATE(atoms%relax(3,atoms%ntype)) - ALLOCATE(atoms%neq(atoms%ntype)) - ALLOCATE(atoms%taual(3,atoms%nat)) - ALLOCATE(atoms%label(atoms%nat)) - ALLOCATE(atoms%pos(3,atoms%nat)) - ALLOCATE(atoms%rmt(atoms%ntype)) - ALLOCATE(atoms%numStatesProvided(atoms%ntype)) - ALLOCATE(atoms%namex(atoms%ntype)) - ALLOCATE(atoms%icorr(atoms%ntype)) - ALLOCATE(atoms%igrd(atoms%ntype)) - ALLOCATE(atoms%krla(atoms%ntype)) - ALLOCATE(atoms%relcor(atoms%ntype)) - - atoms%namex = '' - atoms%icorr = -99 - - ALLOCATE(atoms%ncv(atoms%ntype)) ! For what is this? - ALLOCATE(atoms%ngopr(atoms%nat)) ! For what is this? - ALLOCATE(atoms%lapw_l(atoms%ntype)) ! Where do I put this? - ALLOCATE(atoms%invsat(atoms%nat)) ! Where do I put this? - - ALLOCATE(noco%l_relax(atoms%ntype),noco%b_con(2,atoms%ntype)) - ALLOCATE(noco%alphInit(atoms%ntype),noco%alph(atoms%ntype),noco%beta(atoms%ntype)) - ALLOCATE(noco%socscale(atoms%ntype)) - - DEALLOCATE(atomTypeSpecies,speciesRepAtomType) - ALLOCATE(atomTypeSpecies(atoms%ntype)) - ALLOCATE(speciesRepAtomType(numSpecies)) - atomTypeSpecies = -1 - speciesRepAtomType = -1 - - DEALLOCATE(xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs) - ALLOCATE(xmlElectronStates(29,atoms%ntype)) - ALLOCATE(xmlPrintCoreStates(29,atoms%ntype)) - ALLOCATE(xmlCoreOccs(2,29,atoms%ntype)) - xmlElectronStates = noState_const - xmlPrintCoreStates = .FALSE. - xmlCoreOccs = 0.0 - - ALLOCATE (kpts%ntetra(4,kpts%ntet),kpts%voltet(kpts%ntet)) - - ALLOCATE (wannAtomList(atoms%nat)) - - ! Read in constants - - xPathA = '/fleurInput/constants/constant' - numberNodes = xmlGetNumberOfNodes(xPathA) - DO i = 1, numberNodes - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)), '[',i,']' - tempReal = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@value')) - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@name') - CALL ASSIGN_var(valueString,tempReal) - END DO + REAL :: dtild1,kmax1,dvac1 + LOGICAL :: l_test + INTEGER, ALLOCATABLE :: jri1(:), lmax1(:) + REAL, ALLOCATABLE :: rmt1(:), dx1(:) + + EXTERNAL prp_xcfft_box + + interface + function dropInputSchema() bind(C, name="dropInputSchema") + use iso_c_binding + INTEGER(c_int) dropInputSchema + end function dropInputSchema + end interface + + errorStatus = 0 + errorStatus = dropInputSchema() + IF(errorStatus.NE.0) THEN + STOP 'Error: Cannot print out FleurInputSchema.xsd' + END IF + + schemaFilename = "FleurInputSchema.xsd"//C_NULL_CHAR + docFilename = "inp.xml"//C_NULL_CHAR + + !TODO! these switches should be in the inp-file + input%l_core_confpot=.true. !former CPP_CORE + input%l_useapw=.false. !former CPP_APW + !WRITE(*,*) 'Start reading of inp.xml file' + CALL xmlInitInterface() + CALL xmlParseSchema(schemaFilename) + CALL xmlParseDoc(docFilename) + CALL xmlValidateDoc() + CALL xmlInitXPath() + + ! Check version of inp.xml + versionString = xmlGetAttributeValue('/fleurInput/@fleurInputVersion') + IF((TRIM(ADJUSTL(versionString)).NE.'0.27').AND.(TRIM(ADJUSTL(versionString)).NE.'0.28').AND.& + (TRIM(ADJUSTL(versionString)).NE.'0.29')) THEN + STOP 'version number of inp.xml file is not compatible with this fleur version' + END IF + + ! Get number of atoms, atom types, and atom species + + numberNodes = xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/relPos') + numberNodes = numberNodes + xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/absPos') + numberNodes = numberNodes + xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/filmPos') + + atoms%nat = numberNodes + + numberNodes = xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup') + + atoms%ntype = numberNodes + + numSpecies = xmlGetNumberOfNodes('/fleurInput/atomSpecies/species') + + ALLOCATE(atoms%nz(atoms%ntype)) !nz and zatom have the same content! + ALLOCATE(atoms%zatom(atoms%ntype)) !nz and zatom have the same content! + ALLOCATE(atoms%jri(atoms%ntype)) + ALLOCATE(atoms%dx(atoms%ntype)) + ALLOCATE(atoms%lmax(atoms%ntype)) + ALLOCATE(atoms%nlo(atoms%ntype)) + ALLOCATE(atoms%ncst(atoms%ntype)) + ALLOCATE(atoms%lnonsph(atoms%ntype)) + ALLOCATE(atoms%nflip(atoms%ntype)) + ALLOCATE(atoms%l_geo(atoms%ntype)) + ALLOCATE(atoms%lda_u(4*atoms%ntype)) + ALLOCATE(atoms%bmu(atoms%ntype)) + ALLOCATE(atoms%relax(3,atoms%ntype)) + ALLOCATE(atoms%neq(atoms%ntype)) + ALLOCATE(atoms%taual(3,atoms%nat)) + ALLOCATE(atoms%label(atoms%nat)) + ALLOCATE(atoms%pos(3,atoms%nat)) + ALLOCATE(atoms%rmt(atoms%ntype)) + ALLOCATE(atoms%numStatesProvided(atoms%ntype)) + ALLOCATE(atoms%namex(atoms%ntype)) + ALLOCATE(atoms%icorr(atoms%ntype)) + ALLOCATE(atoms%igrd(atoms%ntype)) + ALLOCATE(atoms%krla(atoms%ntype)) + ALLOCATE(atoms%relcor(atoms%ntype)) + + atoms%namex = '' + atoms%icorr = -99 + + ALLOCATE(atoms%ncv(atoms%ntype)) ! For what is this? + ALLOCATE(atoms%ngopr(atoms%nat)) ! For what is this? + ALLOCATE(atoms%lapw_l(atoms%ntype)) ! Where do I put this? + ALLOCATE(atoms%invsat(atoms%nat)) ! Where do I put this? + + ALLOCATE(noco%l_relax(atoms%ntype),noco%b_con(2,atoms%ntype)) + ALLOCATE(noco%alphInit(atoms%ntype),noco%alph(atoms%ntype),noco%beta(atoms%ntype)) + ALLOCATE(noco%socscale(atoms%ntype)) + + DEALLOCATE(atomTypeSpecies,speciesRepAtomType) + ALLOCATE(atomTypeSpecies(atoms%ntype)) + ALLOCATE(speciesRepAtomType(numSpecies)) + atomTypeSpecies = -1 + speciesRepAtomType = -1 + + DEALLOCATE(xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs) + ALLOCATE(xmlElectronStates(29,atoms%ntype)) + ALLOCATE(xmlPrintCoreStates(29,atoms%ntype)) + ALLOCATE(xmlCoreOccs(2,29,atoms%ntype)) + xmlElectronStates = noState_const + xmlPrintCoreStates = .FALSE. + xmlCoreOccs = 0.0 + + ALLOCATE (kpts%ntetra(4,kpts%ntet),kpts%voltet(kpts%ntet)) + + ALLOCATE (wannAtomList(atoms%nat)) + + ! Read in constants + + xPathA = '/fleurInput/constants/constant' + numberNodes = xmlGetNumberOfNodes(xPathA) + DO i = 1, numberNodes + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)), '[',i,']' + tempReal = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@value')) + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@name') + CALL ASSIGN_var(valueString,tempReal) + END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Comment section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - input%comment = ' ' - xPathA = '/fleurInput/comment' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) - DO i = 1, LEN(TRIM(ADJUSTL(valueString))) - IF (valueString(i:i).EQ.achar(10)) valueString(i:i) = ' ' !remove line breaks - END DO - valueString = TRIM(ADJUSTL(valueString)) - DO i = 1, 10 - j = (i-1) * 8 + 1 - input%comment(i) = valueString(j:j+7) - END DO + input%comment = ' ' + xPathA = '/fleurInput/comment' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) + DO i = 1, LEN(TRIM(ADJUSTL(valueString))) + IF (valueString(i:i).EQ.achar(10)) valueString(i:i) = ' ' !remove line breaks + END DO + valueString = TRIM(ADJUSTL(valueString)) + DO i = 1, 10 + j = (i-1) * 8 + 1 + input%comment(i) = valueString(j:j+7) + END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Start of calculationSetup section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Read general cutoff parameters - - input%rkmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Kmax')) - stars%gmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Gmax')) - - - stars%gmaxInit = stars%gmax - - xPathA = '/fleurInput/calculationSetup/cutoffs/@numbands' - numberNodes = xmlGetNumberOfNodes(xPathA) - dimension%neigd = 0 - IF(numberNodes.EQ.1) THEN - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) - IF(TRIM(ADJUSTL(valueString)).EQ.'all') THEN - STOP 'Feature to calculate all eigenfunctions not yet implemented.' - ELSE - READ(valueString,*) dimension%neigd - END IF - END IF - - ! Read SCF loop parametrization - - input%itmax = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@itmax')) - input%minDistance = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@minDistance')) - input%maxiter = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@maxIterBroyd')) - - valueString = TRIM(ADJUSTL(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@imix'))) - SELECT CASE (valueString) - CASE ('straight') - input%imix = 0 - CASE ('Broyden1') - input%imix = 3 - CASE ('Broyden2') - input%imix = 5 - CASE ('Anderson') - input%imix = 7 - CASE DEFAULT - STOP 'Error: unknown mixing scheme selected!' - END SELECT - - input%alpha = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@alpha')) - input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@preconditioning_param')) - input%spinf = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@spinf')) - - ! Get parameters for core electrons - - input%ctail = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@ctail')) - IF((TRIM(ADJUSTL(versionString)).EQ.'0.27')) THEN - input%coretail_lmax = 99 - ELSE - input%coretail_lmax = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@coretail_lmax')) - END IF - input%frcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@frcor')) - input%kcrel = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@kcrel')) - - ! Read in magnetism parameters - - input%jspins = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@jspins')) - noco%l_noco = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@l_noco')) - input%swsp = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@swsp')) - input%lflip = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@lflip')) - input%fixed_moment=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@fixed_moment')) + ! Read general cutoff parameters + + input%rkmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Kmax')) + stars%gmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Gmax')) + + stars%gmaxInit = stars%gmax + + xPathA = '/fleurInput/calculationSetup/cutoffs/@numbands' + numberNodes = xmlGetNumberOfNodes(xPathA) + dimension%neigd = 0 + IF(numberNodes.EQ.1) THEN + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) + IF(TRIM(ADJUSTL(valueString)).EQ.'all') THEN + STOP 'Feature to calculate all eigenfunctions not yet implemented.' + ELSE + READ(valueString,*) dimension%neigd + END IF + END IF + + ! Read SCF loop parametrization + + input%itmax = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@itmax')) + input%minDistance = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@minDistance')) + input%maxiter = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@maxIterBroyd')) + + valueString = TRIM(ADJUSTL(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@imix'))) + SELECT CASE (valueString) + CASE ('straight') + input%imix = 0 + CASE ('Broyden1') + input%imix = 3 + CASE ('Broyden2') + input%imix = 5 + CASE ('Anderson') + input%imix = 7 + CASE DEFAULT + STOP 'Error: unknown mixing scheme selected!' + END SELECT + + input%alpha = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@alpha')) +input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@preconditioning_param')) + input%spinf = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@spinf')) + + ! Get parameters for core electrons + + input%ctail = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@ctail')) + IF((TRIM(ADJUSTL(versionString)).EQ.'0.27')) THEN + input%coretail_lmax = 99 + ELSE + input%coretail_lmax = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@coretail_lmax')) + END IF + input%frcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@frcor')) + input%kcrel = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@kcrel')) + + ! Read in magnetism parameters + + input%jspins = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@jspins')) + noco%l_noco = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@l_noco')) + input%swsp = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@swsp')) + input%lflip = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@lflip')) + input%fixed_moment=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@fixed_moment')) IF (ABS(input%fixed_moment)>1E-8.AND.(input%jspins==1.OR.noco%l_noco)) CALL judft_error("Fixed moment only in collinear calculations with two spins") - dimension%jspd = input%jspins - - ! Read in Brillouin zone integration parameters - - kpts%nkpt3 = 0 - l_kpts = .FALSE. - - valueString = TRIM(ADJUSTL(xmlGetAttributeValue('/fleurInput/calculationSetup/bzIntegration/@mode'))) - SELECT CASE (valueString) - CASE ('hist') - input%gauss = .FALSE. - input%tria = .FALSE. - CASE ('gauss') - input%gauss = .TRUE. - input%tria = .FALSE. - CASE ('tria') - input%gauss = .FALSE. - input%tria = .TRUE. - CASE DEFAULT - STOP 'Invalid bzIntegration mode selected!' - END SELECT - - nodeSum = 0 - xPathA = '/fleurInput/calculationSetup/bzIntegration/@fermiSmearingEnergy' - numberNodes = xmlGetNumberOfNodes(xPathA) - nodeSum = nodeSum + numberNodes - IF (numberNodes.EQ.1) THEN - input%tkb = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) - END IF - xPathA = '/fleurInput/calculationSetup/bzIntegration/@fermiSmearingTemp' - numberNodes = xmlGetNumberOfNodes(xPathA) - nodeSum = nodeSum + numberNodes - IF (numberNodes.EQ.1) THEN - input%tkb = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) - input%tkb = boltzmannConst * input%tkb - END IF - IF(nodeSum.GE.2) THEN - STOP 'Error: Multiple fermi Smearing parameters provided in input file!' - END IF - - xPathA = '/fleurInput/calculationSetup/bzIntegration/@valenceElectrons' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - input%zelec = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) - ELSE - STOP 'Error: Optionality of valence electrons in input file not yet implemented!' - END IF - - ! Option kPointDensity - kpts%kPointDensity(:) = 0.0 - xPathA = '/fleurInput/calculationSetup/bzIntegration/kPointDensity' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - l_kpts = .FALSE. - kpts%kPointDensity(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denX')) - kpts%kPointDensity(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denY')) - kpts%kPointDensity(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denZ')) - kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma')) - kpts%specificationType = 4 - END IF - - ! Option kPointMesh - xPathA = '/fleurInput/calculationSetup/bzIntegration/kPointMesh' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - l_kpts = .FALSE. - kpts%nkpt3(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nx')) - kpts%nkpt3(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ny')) - kpts%nkpt3(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nz')) - kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma')) - kpts%nkpt = kpts%nkpt3(1) * kpts%nkpt3(2) * kpts%nkpt3(3) - kpts%specificationType = 2 - END IF - - ! Option kPointCount - xPathA = '/fleurInput/calculationSetup/bzIntegration/kPointCount' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - l_kpts = .FALSE. - kpts%nkpt = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@count')) - kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma')) - kpts%nkpt = kpts%nkpt - kpts%specificationType = 1 - - ALLOCATE(kpts%bk(3,kpts%nkpt)) - ALLOCATE(kpts%wtkpt(kpts%nkpt)) - kpts%bk = 0.0 - kpts%wtkpt = 0.0 - kpts%posScale = 1.0 - - numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/bzIntegration/kPointCount/specialPoint') - IF(numberNodes.EQ.1) THEN - STOP 'Error: Single special k point provided. This does not make sense!' - END IF - kpts%numSpecialPoints = numberNodes - IF(kpts%numSpecialPoints.GE.2) THEN - DEALLOCATE(kpts%specialPoints) - ALLOCATE(kpts%specialPoints(3,kpts%numSpecialPoints)) - ALLOCATE(kpts%specialPointNames(kpts%numSpecialPoints)) - DO i = 1, kpts%numSpecialPoints - WRITE(xPathA,*) '/fleurInput/calculationSetup/bzIntegration/kPointCount/specialPoint[',i,']' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) - kpts%specialPoints(1,i) = evaluatefirst(valueString) - kpts%specialPoints(2,i) = evaluatefirst(valueString) - kpts%specialPoints(3,i) = evaluatefirst(valueString) - kpts%specialPointNames(i) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@name') - END DO - END IF - ELSE - DEALLOCATE(kpts%specialPoints) - ALLOCATE(kpts%specialPoints(3,kpts%numSpecialPoints)) - ALLOCATE(kpts%specialPointNames(kpts%numSpecialPoints)) - END IF - - ! Option kPointList - numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/bzIntegration/kPointList') - IF (numberNodes.EQ.1) THEN - l_kpts = .TRUE. - numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/bzIntegration/kPointList/kPoint') - kpts%nkpt = numberNodes - kpts%l_gamma = .FALSE. - ALLOCATE(kpts%bk(3,kpts%nkpt)) - ALLOCATE(kpts%wtkpt(kpts%nkpt)) - kpts%bk = 0.0 - kpts%wtkpt = 0.0 - kpts%specificationType = 3 - - kpts%posScale = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/bzIntegration/kPointList/@posScale')) - weightScale = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/bzIntegration/kPointList/@weightScale')) - - DO i = 1, kpts%nkpt - WRITE(xPathA,*) '/fleurInput/calculationSetup/bzIntegration/kPointList/kPoint[',i,']' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) - READ(valueString,*) kpts%bk(1,i), kpts%bk(2,i), kpts%bk(3,i) - kpts%bk(:,i)=kpts%bk(:,i)/kpts%posScale - kpts%wtkpt(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@weight')) - kpts%wtkpt(i) = kpts%wtkpt(i) / weightScale - END DO - kpts%posScale=1.0 - END IF - - ! Read in optional SOC parameters if present - - xPathA = '/fleurInput/calculationSetup/soc' - numberNodes = xmlGetNumberOfNodes(xPathA) - - noco%l_soc = .FALSE. - noco%theta = 0.0 - noco%phi = 0.0 - - IF (numberNodes.EQ.1) THEN - noco%theta=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta')) - noco%phi=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi')) - noco%l_soc = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_soc')) - noco%l_spav = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spav')) - END IF - - ! Read in optional noco parameters if present - - xPathA = '/fleurInput/calculationSetup/nocoParams' - numberNodes = xmlGetNumberOfNodes(xPathA) - - noco%l_ss = .FALSE. - noco%l_mperp = .FALSE. - noco%l_constr = .FALSE. - noco%mix_b = 0.0 - noco%qss = 0.0 - - noco%l_relax(:) = .FALSE. - noco%alphInit(:) = 0.0 - noco%alph(:) = 0.0 - noco%beta(:) = 0.0 - noco%b_con(:,:) = 0.0 - - - IF ((noco%l_noco).AND.(numberNodes.EQ.0)) THEN - STOP 'Error: l_noco is true but no noco parameters set in xml input file!' - END IF - - IF (numberNodes.EQ.1) THEN - noco%l_ss = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_ss')) - noco%l_mperp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_mperp')) - noco%l_constr = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_constr')) - - noco%mix_b = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mix_b')) - - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/qss'))) - READ(valueString,*) noco%qss(1), noco%qss(2), noco%qss(3) - - !WRITE(*,*) 'Note: TODO: Calculation of q points!' - - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/qsc') - IF (numberNodes.EQ.1) THEN - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/qsc'))) - READ(valueString,*) qsc(1), qsc(2), qsc(3) - DO i = 1, 3 - noco%qss(i) = noco%qss(i) / qsc(i) - END DO - !WRITE(*,*) 'Note: TODO: Integrate qsc directly into qss in input file!' - !WRITE(*,*) '(no problem for users)' - END IF - END IF - - ! Read in optional 1D parameters if present - - xPathA = '/fleurInput/calculationSetup/oneDParams' - numberNodes = xmlGetNumberOfNodes(xPathA) - - oneD%odd%d1 = .FALSE. - - IF (numberNodes.EQ.1) THEN - oneD%odd%d1 = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@d1')) - oneD%odd%M = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@MM')) - oneD%odd%mb = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vM')) - oneD%odd%m_cyl = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@m_cyl')) - oneD%odd%chi = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@chi')) - oneD%odd%rot = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@rot')) - oneD%odd%invs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@invs1')) - oneD%odd%zrfs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zrfs1')) - END IF - - ! Read in optional expert modes switches - - xPathA = '/fleurInput/calculationSetup/expertModes' - numberNodes = xmlGetNumberOfNodes(xPathA) - - input%gw = 0 - input%isec1 = 999999 - input%secvar = .FALSE. - - IF (numberNodes.EQ.1) THEN - input%gw = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gw')) - input%isec1 = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@isec1')) - input%secvar = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@secvar')) - END IF - - ! Read in optional geometry optimization parameters - - xPathA = '/fleurInput/calculationSetup/geometryOptimization' - numberNodes = xmlGetNumberOfNodes(xPathA) - - input%l_f = .FALSE. - input%qfix = 0 - - IF (numberNodes.EQ.1) THEN - input%l_f = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_f')) - input%xa = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@xa')) - input%thetad = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@thetad')) - input%epsdisp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsdisp')) - input%epsforce = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsforce')) - - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@qfix') - IF (numberNodes.EQ.1) THEN - input%qfix = 1 - l_qfix = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@qfix')) - IF (l_qfix) THEN - input%qfix = 2 - END IF - END IF - END IF - - ! Read in optional general LDA+U parameters - - xPathA = '/fleurInput/calculationSetup/ldaU' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - input%ldauLinMix = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_linMix')) - input%ldauMixParam = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mixParam')) - input%ldauSpinf = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spinf')) - END IF - - ! Read in optional q point mesh for spin spirals - - xPathA = '/fleurInput/calculationSetup/spinSpiralQPointMesh' - numberNodes = xmlGetNumberOfNodes(xPathA) - - ! IF ((noco%l_ss).AND.(numberNodes.EQ.0)) THEN - ! STOP 'Error: l_ss is true but no q point mesh set in xml input file!' - ! END IF - - - ! Read in optional E-Field parameters - - xPathA = '/fleurInput/calculationSetup/eField' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - !input%efield%zsigma = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zsigma')) - !input%efield%sig_b(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_1')) - !input%efield%sig_b(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_2')) - !input%efield%plot_charge = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_charge')) - !input%efield%plot_rho = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_rho')) - !input%efield%autocomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@autocomp')) - !input%efield%dirichlet = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dirichlet')) - !l_eV = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eV')) - - STOP 'Error: Reading input for E-Fields not yet implemented completely!' - ! ALLOCATE(input%efield%sigEF(3*k1d, 3*k2d, nvac)) - ! input%efield%sigEF = 0.0 - !IF (l_eV) THEN - ! input%efield%sig_b(:) = input%efield%sig_b/hartree_to_ev_const - ! input%efield%sigEF(:,:,:) = input%efield%sigEF/hartree_to_ev_const - !END IF - END IF - - ! Read in optional energy parameter limits - - xPathA = '/fleurInput/calculationSetup/energyParameterLimits' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - input%ellow = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ellow')) - input%elup = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@elup')) - END IF + dimension%jspd = input%jspins + + ! Read in Brillouin zone integration parameters + + kpts%nkpt3 = 0 + l_kpts = .FALSE. + + valueString = TRIM(ADJUSTL(xmlGetAttributeValue('/fleurInput/calculationSetup/bzIntegration/@mode'))) + SELECT CASE (valueString) + CASE ('hist') + input%gauss = .FALSE. + input%tria = .FALSE. + CASE ('gauss') + input%gauss = .TRUE. + input%tria = .FALSE. + CASE ('tria') + input%gauss = .FALSE. + input%tria = .TRUE. + CASE DEFAULT + STOP 'Invalid bzIntegration mode selected!' + END SELECT + + nodeSum = 0 + xPathA = '/fleurInput/calculationSetup/bzIntegration/@fermiSmearingEnergy' + numberNodes = xmlGetNumberOfNodes(xPathA) + nodeSum = nodeSum + numberNodes + IF (numberNodes.EQ.1) THEN + input%tkb = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) + END IF + xPathA = '/fleurInput/calculationSetup/bzIntegration/@fermiSmearingTemp' + numberNodes = xmlGetNumberOfNodes(xPathA) + nodeSum = nodeSum + numberNodes + IF (numberNodes.EQ.1) THEN + input%tkb = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) + input%tkb = boltzmannConst * input%tkb + END IF + IF(nodeSum.GE.2) THEN + STOP 'Error: Multiple fermi Smearing parameters provided in input file!' + END IF + + xPathA = '/fleurInput/calculationSetup/bzIntegration/@valenceElectrons' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + input%zelec = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) + ELSE + STOP 'Error: Optionality of valence electrons in input file not yet implemented!' + END IF + + ! Option kPointDensity + kpts%kPointDensity(:) = 0.0 + xPathA = '/fleurInput/calculationSetup/bzIntegration/kPointDensity' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + l_kpts = .FALSE. + kpts%kPointDensity(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denX')) + kpts%kPointDensity(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denY')) + kpts%kPointDensity(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denZ')) + kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma')) + kpts%specificationType = 4 + END IF + + ! Option kPointMesh + xPathA = '/fleurInput/calculationSetup/bzIntegration/kPointMesh' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + l_kpts = .FALSE. + kpts%nkpt3(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nx')) + kpts%nkpt3(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ny')) + kpts%nkpt3(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nz')) + kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma')) + kpts%nkpt = kpts%nkpt3(1) * kpts%nkpt3(2) * kpts%nkpt3(3) + kpts%specificationType = 2 + END IF + + ! Option kPointCount + xPathA = '/fleurInput/calculationSetup/bzIntegration/kPointCount' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + l_kpts = .FALSE. + kpts%nkpt = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@count')) + kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma')) + kpts%nkpt = kpts%nkpt + kpts%specificationType = 1 + + ALLOCATE(kpts%bk(3,kpts%nkpt)) + ALLOCATE(kpts%wtkpt(kpts%nkpt)) + kpts%bk = 0.0 + kpts%wtkpt = 0.0 + kpts%posScale = 1.0 + + numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/bzIntegration/kPointCount/specialPoint') + IF(numberNodes.EQ.1) THEN + STOP 'Error: Single special k point provided. This does not make sense!' + END IF + kpts%numSpecialPoints = numberNodes + IF(kpts%numSpecialPoints.GE.2) THEN + DEALLOCATE(kpts%specialPoints) + ALLOCATE(kpts%specialPoints(3,kpts%numSpecialPoints)) + ALLOCATE(kpts%specialPointNames(kpts%numSpecialPoints)) + DO i = 1, kpts%numSpecialPoints + WRITE(xPathA,*) '/fleurInput/calculationSetup/bzIntegration/kPointCount/specialPoint[',i,']' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) + kpts%specialPoints(1,i) = evaluatefirst(valueString) + kpts%specialPoints(2,i) = evaluatefirst(valueString) + kpts%specialPoints(3,i) = evaluatefirst(valueString) + kpts%specialPointNames(i) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@name') + END DO + END IF + ELSE + DEALLOCATE(kpts%specialPoints) + ALLOCATE(kpts%specialPoints(3,kpts%numSpecialPoints)) + ALLOCATE(kpts%specialPointNames(kpts%numSpecialPoints)) + END IF + + ! Option kPointList + numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/bzIntegration/kPointList') + IF (numberNodes.EQ.1) THEN + l_kpts = .TRUE. + numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/bzIntegration/kPointList/kPoint') + kpts%nkpt = numberNodes + kpts%l_gamma = .FALSE. + ALLOCATE(kpts%bk(3,kpts%nkpt)) + ALLOCATE(kpts%wtkpt(kpts%nkpt)) + kpts%bk = 0.0 + kpts%wtkpt = 0.0 + kpts%specificationType = 3 + + kpts%posScale = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/bzIntegration/kPointList/@posScale')) + weightScale = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/bzIntegration/kPointList/@weightScale')) + + DO i = 1, kpts%nkpt + WRITE(xPathA,*) '/fleurInput/calculationSetup/bzIntegration/kPointList/kPoint[',i,']' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) + READ(valueString,*) kpts%bk(1,i), kpts%bk(2,i), kpts%bk(3,i) + kpts%bk(:,i)=kpts%bk(:,i)/kpts%posScale + kpts%wtkpt(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@weight')) + kpts%wtkpt(i) = kpts%wtkpt(i) / weightScale + END DO + kpts%posScale=1.0 + END IF + + ! Read in optional SOC parameters if present + + xPathA = '/fleurInput/calculationSetup/soc' + numberNodes = xmlGetNumberOfNodes(xPathA) + + noco%l_soc = .FALSE. + noco%theta = 0.0 + noco%phi = 0.0 + + IF (numberNodes.EQ.1) THEN + noco%theta=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta')) + noco%phi=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi')) + noco%l_soc = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_soc')) + noco%l_spav = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spav')) + END IF + + ! Read in optional noco parameters if present + + xPathA = '/fleurInput/calculationSetup/nocoParams' + numberNodes = xmlGetNumberOfNodes(xPathA) + + noco%l_ss = .FALSE. + noco%l_mperp = .FALSE. + noco%l_constr = .FALSE. + noco%mix_b = 0.0 + noco%qss = 0.0 + + noco%l_relax(:) = .FALSE. + noco%alphInit(:) = 0.0 + noco%alph(:) = 0.0 + noco%beta(:) = 0.0 + noco%b_con(:,:) = 0.0 + + IF ((noco%l_noco).AND.(numberNodes.EQ.0)) THEN + STOP 'Error: l_noco is true but no noco parameters set in xml input file!' + END IF + + IF (numberNodes.EQ.1) THEN + noco%l_ss = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_ss')) + noco%l_mperp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_mperp')) + noco%l_constr = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_constr')) + + noco%mix_b = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mix_b')) + + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/qss'))) + READ(valueString,*) noco%qss(1), noco%qss(2), noco%qss(3) + + !WRITE(*,*) 'Note: TODO: Calculation of q points!' + + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/qsc') + IF (numberNodes.EQ.1) THEN + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/qsc'))) + READ(valueString,*) qsc(1), qsc(2), qsc(3) + DO i = 1, 3 + noco%qss(i) = noco%qss(i) / qsc(i) + END DO + !WRITE(*,*) 'Note: TODO: Integrate qsc directly into qss in input file!' + !WRITE(*,*) '(no problem for users)' + END IF + END IF + + ! Read in optional 1D parameters if present + + xPathA = '/fleurInput/calculationSetup/oneDParams' + numberNodes = xmlGetNumberOfNodes(xPathA) + + oneD%odd%d1 = .FALSE. + + IF (numberNodes.EQ.1) THEN + oneD%odd%d1 = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@d1')) + oneD%odd%M = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@MM')) + oneD%odd%mb = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vM')) + oneD%odd%m_cyl = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@m_cyl')) + oneD%odd%chi = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@chi')) + oneD%odd%rot = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@rot')) + oneD%odd%invs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@invs1')) + oneD%odd%zrfs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zrfs1')) + END IF + + ! Read in optional expert modes switches + + xPathA = '/fleurInput/calculationSetup/expertModes' + numberNodes = xmlGetNumberOfNodes(xPathA) + + input%gw = 0 + input%isec1 = 999999 + input%secvar = .FALSE. + + IF (numberNodes.EQ.1) THEN + input%gw = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gw')) + input%isec1 = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@isec1')) + input%secvar = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@secvar')) + END IF + + ! Read in optional geometry optimization parameters + + xPathA = '/fleurInput/calculationSetup/geometryOptimization' + numberNodes = xmlGetNumberOfNodes(xPathA) + + input%l_f = .FALSE. + input%qfix = 0 + + IF (numberNodes.EQ.1) THEN + input%l_f = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_f')) + input%xa = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@xa')) + input%thetad = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@thetad')) + input%epsdisp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsdisp')) + input%epsforce = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsforce')) + + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@qfix') + IF (numberNodes.EQ.1) THEN + input%qfix = 1 + l_qfix = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@qfix')) + IF (l_qfix) THEN + input%qfix = 2 + END IF + END IF + END IF + + ! Read in optional general LDA+U parameters + + xPathA = '/fleurInput/calculationSetup/ldaU' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + input%ldauLinMix = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_linMix')) + input%ldauMixParam = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mixParam')) + input%ldauSpinf = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spinf')) + END IF + + ! Read in optional q point mesh for spin spirals + + xPathA = '/fleurInput/calculationSetup/spinSpiralQPointMesh' + numberNodes = xmlGetNumberOfNodes(xPathA) + + ! IF ((noco%l_ss).AND.(numberNodes.EQ.0)) THEN + ! STOP 'Error: l_ss is true but no q point mesh set in xml input file!' + ! END IF + + ! Read in optional E-Field parameters + + xPathA = '/fleurInput/calculationSetup/eField' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + !input%efield%zsigma = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zsigma')) + !input%efield%sig_b(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_1')) + !input%efield%sig_b(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_2')) + !input%efield%plot_charge = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_charge')) + !input%efield%plot_rho = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_rho')) + !input%efield%autocomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@autocomp')) + !input%efield%dirichlet = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dirichlet')) + !l_eV = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eV')) + + STOP 'Error: Reading input for E-Fields not yet implemented completely!' + ! ALLOCATE(input%efield%sigEF(3*k1d, 3*k2d, nvac)) + ! input%efield%sigEF = 0.0 + !IF (l_eV) THEN + ! input%efield%sig_b(:) = input%efield%sig_b/hartree_to_ev_const + ! input%efield%sigEF(:,:,:) = input%efield%sigEF/hartree_to_ev_const + !END IF + END IF + + ! Read in optional energy parameter limits + + xPathA = '/fleurInput/calculationSetup/energyParameterLimits' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + input%ellow = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ellow')) + input%elup = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@elup')) + END IF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! End of calculationSetup section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Start of cell section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Read in lattice parameters - - a1 = 0.0 - a2 = 0.0 - a3 = 0.0 - cell%z1 = 0.0 - dtild = 0.0 - input%film = .FALSE. - latticeType = 'bulkLattice' - latticeDef = 0 - symmetryDef = 0 - cell%latnam = 'any' - - numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/filmLattice') - - IF (numberNodes.EQ.1) THEN - input%film = .TRUE. - latticeType = 'filmLattice' - END IF - - xPathA = '/fleurInput/cell/'//latticeType - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - latticeScale = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@scale')) - input%scaleCell = latticeScale - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@latnam'))) - READ(valueString,*) cell%latnam - - IF(input%film) THEN - cell%z1 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dVac')) - dtild = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dTilda')) - vacuum%dvac = cell%z1 - a3(3) = dtild - evac0Temp = eVac0Default_const - xPathB = TRIM(ADJUSTL(xPathA))//'/vacuumEnergyParameters' - numberNodes = xmlGetNumberOfNodes(xPathB) - IF(numberNodes.GE.1) THEN - DO i = 1, numberNodes - xPathC = '' - WRITE(xPathC,'(a,i0,a)') TRIM(ADJUSTL(xPathB))//'[',i,']' - numVac = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@vacuum')) - eParamUp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@spinUp')) - eParamDown = eParamUp - IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathC))//'/@spinDown').GE.1) THEN - eParamDown = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@spinDown')) - END IF - evac0Temp(numVac,1) = eParamUp - IF(input%jspins.GT.1) evac0Temp(numVac,2) = eParamDown - IF(i.EQ.1) THEN - evac0Temp(3-numVac,1) = eParamUp - IF(input%jspins.GT.1) evac0Temp(3-numVac,2) = eParamDown - END IF - END DO - END IF - END IF - - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/a1') - IF (numberNodes.EQ.1) THEN - latticeDef = 1 - input%scaleA1 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a1/@scale')) - a1(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a1')) - a1(1) = a1(1) * input%scaleA1 - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/a2') - IF (numberNodes.EQ.1) THEN - latticeDef = 2 - input%scaleA2 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a2/@scale')) - a2(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a2')) - a2(2) = a2(2) * input%scaleA2 - END IF - IF(.NOT.input%film) THEN - input%scaleC = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c/@scale')) - a3(3) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c')) - a3(3) = a3(3) * input%scaleC - END IF - END IF - - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/row-1') - IF (numberNodes.EQ.1) THEN - latticeDef = 3 - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/row-1'))) - a1(1) = evaluateFirst(valueString) - a1(2) = evaluateFirst(valueString) - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/row-2'))) - a2(1) = evaluateFirst(valueString) - a2(2) = evaluateFirst(valueString) - IF(.NOT.input%film) THEN - input%scaleC = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c/@scale')) - a3(3) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c')) - a3(3) = a3(3) * input%scaleC - END IF - END IF - - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix') - IF (numberNodes.EQ.1) THEN - latticeDef = 4 - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-1'))) - a1(1) = evaluateFirst(valueString) - a1(2) = evaluateFirst(valueString) - IF(.NOT.input%film) THEN - a1(3) = evaluateFirst(valueString) - END IF - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-2'))) - a2(1) = evaluateFirst(valueString) - a2(2) = evaluateFirst(valueString) - IF(.NOT.input%film) THEN - a2(3) = evaluateFirst(valueString) - END IF - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-3'))) - IF(.NOT.input%film) THEN - a3(1) = evaluateFirst(valueString) - a3(2) = evaluateFirst(valueString) - a3(3) = evaluateFirst(valueString) - ELSE - !WRITE(*,*) 'Note: For film calculations only the upper left 2x2 part of the Bravais matrix is considered.' - END IF - END IF - END IF ! Note: Further ways to define lattices might be added later. (1D lattice,...) - - ! Construction of amat requires additional information about the lattice - ! and is done later (scroll down)! - - ! Read in symmetry parameters - - sym%namgrp = 'any' - - xPathA = '/fleurInput/cell/symmetry' - numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/symmetry') - - IF (numberNodes.EQ.1) THEN - sym%symSpecType = 2 - symmetryDef = 1 - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spgrp'))) - READ(valueString,*) sym%namgrp - sym%invs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@invs')) - sym%zrfs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zrfs')) - sym%invs2 = sym%invs.AND.sym%zrfs - - IF (sym%namgrp.EQ.'any ') THEN - sym%nop = 48 - ! Read in sym.out file if sym%namgrp='any' set. - CALL rw_symfile('r',94,'sym.out',48,cell%bmat,& - & mrotTemp,tauTemp,sym%nop,sym%nop2,sym%symor) - IF (ALLOCATED(sym%mrot)) THEN - DEALLOCATE(sym%mrot) - END IF - ALLOCATE(sym%mrot(3,3,sym%nop)) - IF (ALLOCATED(sym%tau)) THEN - DEALLOCATE(sym%tau) - END IF - ALLOCATE(sym%tau(3,sym%nop)) - - DO k = 1, sym%nop - DO i = 1, 3 - DO j = 1, 3 - sym%mrot(j,i,k) = mrotTemp(j,i,k) - END DO - sym%tau(i,k) = tauTemp(i,k) - END DO - END DO - ELSE - n2spg = 0 - DO i = 1, 20 - IF (sym%namgrp.EQ.nammap(i)) n2spg = i - END DO - IF (n2spg == 0 ) THEN - WRITE (*,*) 'Spacegroup ',sym%namgrp,' not known! Choose one of:' - WRITE (*,'(20(a4,1x))') (nammap(i),i=1,20) - CALL juDFT_error("Could not determine spacegroup!", calledby = "r_inpXML") - END IF - IF ((n2spg.GE.13).AND.(n2spg.LE.17)) THEN - IF (.not.((cell%latnam.EQ.'hx3').OR.(cell%latnam.EQ.'hex'))) THEN - CALL juDFT_error("Use only hex or hx3 with p3, p3m1, p31m, p6 or p6m!", calledby ="r_inpXML") - END IF - END IF - sym%nop = ord2(n2spg) - IF (sym%invs) THEN - sym%nop = 2*sym%nop - IF (sym%zrfs.and.(.not.l_c2(n2spg))) sym%nop = 2*sym%nop - ELSE - IF (sym%zrfs) sym%nop = 2*sym%nop - END IF - IF (ALLOCATED(sym%mrot)) THEN - DEALLOCATE(sym%mrot) - END IF - ALLOCATE(sym%mrot(3,3,sym%nop)) - IF (ALLOCATED(sym%tau)) THEN - DEALLOCATE(sym%tau) - END IF - ALLOCATE(sym%tau(3,sym%nop)) - CALL spg2set(sym%nop,sym%zrfs,sym%invs,sym%namgrp,cell%latnam,& - & sym%mrot,sym%tau,sym%nop2,sym%symor) - END IF - END IF - - xPathA = '/fleurInput/cell/symmetryFile' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - symmetryDef = 2 - sym%symSpecType = 1 - sym%nop = 48 - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@filename'))) - - CALL rw_symfile('r',94,TRIM(ADJUSTL(valueString)),48,cell%bmat,& - & mrotTemp,tauTemp,sym%nop,sym%nop2,sym%symor) - - IF (ALLOCATED(sym%mrot)) THEN - DEALLOCATE(sym%mrot) - END IF - ALLOCATE(sym%mrot(3,3,sym%nop)) - IF (ALLOCATED(sym%tau)) THEN - DEALLOCATE(sym%tau) - END IF - ALLOCATE(sym%tau(3,sym%nop)) - - sym%invs = .FALSE. - sym%zrfs = .FALSE. - - DO k = 1, sym%nop - absSum = 0 - DO i = 1, 3 - DO j = 1, 3 - sym%mrot(j,i,k) = mrotTemp(j,i,k) - absSum = absSum + ABS(sym%mrot(j,i,k)) - END DO - sym%tau(i,k) = tauTemp(i,k) - END DO - IF (absSum.EQ.3) THEN - IF (ALL(sym%tau(:,k).EQ.0.0)) THEN - IF ((sym%mrot(1,1,k).EQ.-1).AND.(sym%mrot(2,2,k).EQ.-1).AND.(sym%mrot(3,3,k).EQ.-1)) sym%invs = .TRUE. - IF ((sym%mrot(1,1,k).EQ.1).AND.(sym%mrot(2,2,k).EQ.1).AND.(sym%mrot(3,3,k).EQ.-1)) sym%zrfs = .TRUE. - END IF - END IF - END DO - - sym%invs2 = sym%invs.AND.sym%zrfs - END IF - - xPathA = '/fleurInput/cell/symmetryOperations' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - sym%symSpecType = 3 - symmetryDef = 3 - - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/symOp') - sym%nop = numberNodes - - IF (ALLOCATED(sym%mrot)) DEALLOCATE(sym%mrot) - ALLOCATE(sym%mrot(3,3,sym%nop)) - IF (ALLOCATED(sym%tau)) DEALLOCATE(sym%tau) - ALLOCATE(sym%tau(3,sym%nop)) - sym%symor = .TRUE. - DO i = 1, sym%nop - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-1' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))))) - READ(valueString,*) sym%mrot(1,1,i), sym%mrot(1,2,i), sym%mrot(1,3,i), sym%tau(1,i) - - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-2' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))))) - READ(valueString,*) sym%mrot(2,1,i), sym%mrot(2,2,i), sym%mrot(2,3,i), sym%tau(2,i) - - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-3' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))))) - READ(valueString,*) sym%mrot(3,1,i), sym%mrot(3,2,i), sym%mrot(3,3,i), sym%tau(3,i) - - IF ((sym%tau(1,i)**2 + sym%tau(2,i)**2 + sym%tau(3,i)**2).GT.1.e-8) THEN - sym%symor = .FALSE. - END IF - DO j = 1,3 - IF (ABS(sym%tau(j,i)-0.33333) < 0.00001) THEN - sym%tau(j,i) = 1./3. - ENDIF - IF (ABS(sym%tau(j,i)+0.33333) < 0.00001) THEN - sym%tau(j,i) = -1./3. - ENDIF - IF (ABS(sym%tau(j,i)-0.66667) < 0.00001) THEN - sym%tau(j,i) = 2./3. - ENDIF - IF (ABS(sym%tau(j,i)+0.66667) < 0.00001) THEN - sym%tau(j,i) = -2./3. - ENDIF - ENDDO - END DO - END IF - - ! Calculate missing symmetry and cell properties and check consistency of parameters. - - ! Construction of amat - SELECT CASE (latticeDef) - CASE (1) - IF (cell%latnam.EQ.'squ') THEN - a2(2) = a1(1) - ELSE IF (cell%latnam.EQ.'c-b') THEN - aTemp = a1(1) - a1(1) = aTemp*0.5*sqrt(2.0) - a1(2) = -aTemp*0.5 - a2(1) = aTemp*0.5*sqrt(2.0) - a2(2) = aTemp*0.5 - ELSE IF (cell%latnam.EQ.'hex') THEN - aTemp = 0.5*a1(1) - a1(1) = aTemp*sqrt(3.0) - a1(2) = -aTemp - a2(1) = a1(1) - a2(2) = aTemp - ELSE IF (cell%latnam.EQ.'hx3') THEN - aTemp = 0.5*a1(1) - a1(1) = aTemp - a1(2) = -aTemp*sqrt(3.0) - a2(1) = a1(1) - a2(2) = -a1(2) - ELSE IF (cell%latnam.EQ.'fcc') THEN - aTemp = a1(1) - a1(1) = 0.0 ; a1(2) = 0.5*aTemp ; a1(3) = 0.5*aTemp - a2(1) = 0.5*aTemp ; a2(2) = 0.0 ; a2(3) = 0.5*aTemp - a3(1) = 0.5*aTemp ; a3(2) = 0.5*aTemp ; a3(3) = 0.0 - ELSE IF (cell%latnam.EQ.'bcc') THEN - aTemp = a1(1) - a1(1) =-0.5*aTemp ; a1(2) = 0.5*aTemp ; a1(3) = 0.5*aTemp - a2(1) = 0.5*aTemp ; a2(2) =-0.5*aTemp ; a2(3) = 0.5*aTemp - a3(1) = 0.5*aTemp ; a3(2) = 0.5*aTemp ; a3(3) =-0.5*aTemp - ELSE - CALL juDFT_error("latnam is incompatible to parametrization of lattice (1)", calledby ="r_inpXML") - END IF - CASE (2) - IF ((cell%latnam.EQ.'c-r').OR.(cell%latnam.EQ.'p-r')) THEN - IF (cell%latnam.EQ.'c-r') THEN - a1(2) = -a2(2) - a2(1) = a1(1) - END IF - ELSE - CALL juDFT_error("latnam is incompatible to parametrization of lattice (2)", calledby ="r_inpXML") - END IF - CASE (3) - IF (.NOT.(cell%latnam.EQ.'obl')) THEN - CALL juDFT_error("latnam is incompatible to parametrization of lattice (3)", calledby ="r_inpXML") - END IF - CASE (4) - IF (.NOT.(cell%latnam.EQ.'any')) THEN - CALL juDFT_error("latnam is incompatible to parametrization of lattice (4)", calledby ="r_inpXML") - END IF - CASE DEFAULT - CALL juDFT_error("Illegal lattice definition", calledby ="r_inpXML") - END SELECT - - IF (latticeScale.EQ.0.0) latticeScale = 1.0 - IF (.NOT.input%film) vacuum%dvac = a3(3) - vacuum%dvac = latticeScale*vacuum%dvac - dtild = latticeScale*dtild - - cell%amat(:,1) = a1(:) * latticeScale - cell%amat(:,2) = a2(:) * latticeScale - cell%amat(:,3) = a3(:) * latticeScale - - CALL inv3(cell%amat,cell%bmat,cell%omtil) - cell%bmat(:,:) = tpi_const*cell%bmat(:,:) - cell%omtil = abs(cell%omtil) - - IF (input%film.AND..NOT.oneD%odd%d1) THEN - cell%vol = (cell%omtil/cell%amat(3,3))*vacuum%dvac - cell%area = cell%omtil/cell%amat(3,3) - !-odim - ELSE IF (oneD%odd%d1) THEN - cell%area = tpi_const*cell%amat(3,3) - cell%vol = pi_const*(vacuum%dvac**2)*cell%amat(3,3)/4.0 - !+odim - ELSE - cell%vol = cell%omtil - cell%area = cell%amat(1,1)*cell%amat(2,2)-cell%amat(1,2)*cell%amat(2,1) - IF (cell%area.lt.1.0e-7) THEN - IF (cell%latnam.EQ.'any') THEN - cell%area = 1. - ELSE - CALL juDFT_error("area = 0",calledby ="r_inpXML") - END IF - END IF - END IF - - ! Construction of missing symmetry information - IF ((symmetryDef.EQ.2).OR.(symmetryDef.EQ.3)) THEN - nop48 = 48 - ALLOCATE (invOps(sym%nop),multtab(sym%nop,sym%nop),optype(nop48)) - CALL check_close(sym%nop,sym%mrot,sym%tau,& - & multtab,invOps,optype) - - CALL symproperties(nop48,optype,input%film,sym%nop,multtab,cell%amat,& - & sym%symor,sym%mrot,sym%tau,& - & invSym,sym%invs,sym%zrfs,sym%invs2,sym%nop,sym%nop2) - DEALLOCATE(invOps,multtab,optype) - IF (.not.input%film) sym%nop2=sym%nop - IF (input%film) THEN - DO n = 1, sym%nop - DO i = 1, 3 - IF (ABS(sym%tau(i,n)) > 0.00001) THEN - CALL juDFT_error("nonsymmorphic symmetries not yet implemented for films!",calledby ="r_inpXML") - ENDIF - END DO - END DO - END IF - END IF - sym%invs2 = sym%invs.AND.sym%zrfs - - ALLOCATE (sym%invarop(atoms%nat,sym%nop),sym%invarind(atoms%nat)) - ALLOCATE (sym%multab(sym%nop,sym%nop),sym%invtab(sym%nop)) - ALLOCATE (sym%invsatnr(atoms%nat),sym%d_wgn(-3:3,-3:3,3,sym%nop)) - - !some settings for film calculations - vacuum%nmzd = 250 - vacuum%nmzxyd = 100 - vacuum%nvac = 2 - IF (sym%zrfs.OR.sym%invs) vacuum%nvac = 1 - IF (oneD%odd%d1) vacuum%nvac = 1 - cell%z1 = vacuum%dvac/2 - vacuum%nmz = vacuum%nmzd - vacuum%delz = 25.0/vacuum%nmz - IF (oneD%odd%d1) vacuum%delz = 20.0 / vacuum%nmz - IF (vacuum%nmz.GT.vacuum%nmzd) CALL juDFT_error("nmzd",calledby ="inped") - vacuum%nmzxy = vacuum%nmzxyd - IF (vacuum%nmzxy.GT.vacuum%nmzxyd) CALL juDFT_error("nmzxyd",calledby ="inped") + ! Read in lattice parameters + + a1 = 0.0 + a2 = 0.0 + a3 = 0.0 + cell%z1 = 0.0 + dtild = 0.0 + input%film = .FALSE. + latticeType = 'bulkLattice' + latticeDef = 0 + symmetryDef = 0 + cell%latnam = 'any' + + numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/filmLattice') + + IF (numberNodes.EQ.1) THEN + input%film = .TRUE. + latticeType = 'filmLattice' + END IF + + xPathA = '/fleurInput/cell/'//latticeType + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + latticeScale = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@scale')) + input%scaleCell = latticeScale + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@latnam'))) + READ(valueString,*) cell%latnam + + IF(input%film) THEN + cell%z1 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dVac')) + dtild = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dTilda')) + vacuum%dvac = cell%z1 + a3(3) = dtild + evac0Temp = eVac0Default_const + xPathB = TRIM(ADJUSTL(xPathA))//'/vacuumEnergyParameters' + numberNodes = xmlGetNumberOfNodes(xPathB) + IF(numberNodes.GE.1) THEN + DO i = 1, numberNodes + xPathC = '' + WRITE(xPathC,'(a,i0,a)') TRIM(ADJUSTL(xPathB))//'[',i,']' + numVac = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@vacuum')) + eParamUp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@spinUp')) + eParamDown = eParamUp + IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathC))//'/@spinDown').GE.1) THEN + eParamDown = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@spinDown')) + END IF + evac0Temp(numVac,1) = eParamUp + IF(input%jspins.GT.1) evac0Temp(numVac,2) = eParamDown + IF(i.EQ.1) THEN + evac0Temp(3-numVac,1) = eParamUp + IF(input%jspins.GT.1) evac0Temp(3-numVac,2) = eParamDown + END IF + END DO + END IF + END IF + + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/a1') + IF (numberNodes.EQ.1) THEN + latticeDef = 1 + input%scaleA1 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a1/@scale')) + a1(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a1')) + a1(1) = a1(1) * input%scaleA1 + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/a2') + IF (numberNodes.EQ.1) THEN + latticeDef = 2 + input%scaleA2 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a2/@scale')) + a2(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a2')) + a2(2) = a2(2) * input%scaleA2 + END IF + IF(.NOT.input%film) THEN + input%scaleC = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c/@scale')) + a3(3) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c')) + a3(3) = a3(3) * input%scaleC + END IF + END IF + + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/row-1') + IF (numberNodes.EQ.1) THEN + latticeDef = 3 + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/row-1'))) + a1(1) = evaluateFirst(valueString) + a1(2) = evaluateFirst(valueString) + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/row-2'))) + a2(1) = evaluateFirst(valueString) + a2(2) = evaluateFirst(valueString) + IF(.NOT.input%film) THEN + input%scaleC = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c/@scale')) + a3(3) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c')) + a3(3) = a3(3) * input%scaleC + END IF + END IF + + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix') + IF (numberNodes.EQ.1) THEN + latticeDef = 4 + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-1'))) + a1(1) = evaluateFirst(valueString) + a1(2) = evaluateFirst(valueString) + IF(.NOT.input%film) THEN + a1(3) = evaluateFirst(valueString) + END IF + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-2'))) + a2(1) = evaluateFirst(valueString) + a2(2) = evaluateFirst(valueString) + IF(.NOT.input%film) THEN + a2(3) = evaluateFirst(valueString) + END IF + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-3'))) + IF(.NOT.input%film) THEN + a3(1) = evaluateFirst(valueString) + a3(2) = evaluateFirst(valueString) + a3(3) = evaluateFirst(valueString) + ELSE + !WRITE(*,*) 'Note: For film calculations only the upper left 2x2 part of the Bravais matrix is considered.' + END IF + END IF + END IF ! Note: Further ways to define lattices might be added later. (1D lattice,...) + + ! Construction of amat requires additional information about the lattice + ! and is done later (scroll down)! + + ! Read in symmetry parameters + + sym%namgrp = 'any' + + xPathA = '/fleurInput/cell/symmetry' + numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/symmetry') + + IF (numberNodes.EQ.1) THEN + sym%symSpecType = 2 + symmetryDef = 1 + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spgrp'))) + READ(valueString,*) sym%namgrp + sym%invs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@invs')) + sym%zrfs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zrfs')) + sym%invs2 = sym%invs.AND.sym%zrfs + + IF (sym%namgrp.EQ.'any ') THEN + sym%nop = 48 + ! Read in sym.out file if sym%namgrp='any' set. + CALL rw_symfile('r',94,'sym.out',48,cell%bmat,& + & mrotTemp,tauTemp,sym%nop,sym%nop2,sym%symor) + IF (ALLOCATED(sym%mrot)) THEN + DEALLOCATE(sym%mrot) + END IF + ALLOCATE(sym%mrot(3,3,sym%nop)) + IF (ALLOCATED(sym%tau)) THEN + DEALLOCATE(sym%tau) + END IF + ALLOCATE(sym%tau(3,sym%nop)) + + DO k = 1, sym%nop + DO i = 1, 3 + DO j = 1, 3 + sym%mrot(j,i,k) = mrotTemp(j,i,k) + END DO + sym%tau(i,k) = tauTemp(i,k) + END DO + END DO + ELSE + n2spg = 0 + DO i = 1, 20 + IF (sym%namgrp.EQ.nammap(i)) n2spg = i + END DO + IF (n2spg == 0 ) THEN + WRITE (*,*) 'Spacegroup ',sym%namgrp,' not known! Choose one of:' + WRITE (*,'(20(a4,1x))') (nammap(i),i=1,20) + CALL juDFT_error("Could not determine spacegroup!", calledby = "r_inpXML") + END IF + IF ((n2spg.GE.13).AND.(n2spg.LE.17)) THEN + IF (.not.((cell%latnam.EQ.'hx3').OR.(cell%latnam.EQ.'hex'))) THEN + CALL juDFT_error("Use only hex or hx3 with p3, p3m1, p31m, p6 or p6m!", calledby ="r_inpXML") + END IF + END IF + sym%nop = ord2(n2spg) + IF (sym%invs) THEN + sym%nop = 2*sym%nop + IF (sym%zrfs.and.(.not.l_c2(n2spg))) sym%nop = 2*sym%nop + ELSE + IF (sym%zrfs) sym%nop = 2*sym%nop + END IF + IF (ALLOCATED(sym%mrot)) THEN + DEALLOCATE(sym%mrot) + END IF + ALLOCATE(sym%mrot(3,3,sym%nop)) + IF (ALLOCATED(sym%tau)) THEN + DEALLOCATE(sym%tau) + END IF + ALLOCATE(sym%tau(3,sym%nop)) + CALL spg2set(sym%nop,sym%zrfs,sym%invs,sym%namgrp,cell%latnam,& + & sym%mrot,sym%tau,sym%nop2,sym%symor) + END IF + END IF + + xPathA = '/fleurInput/cell/symmetryFile' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + symmetryDef = 2 + sym%symSpecType = 1 + sym%nop = 48 + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@filename'))) + + CALL rw_symfile('r',94,TRIM(ADJUSTL(valueString)),48,cell%bmat,& + & mrotTemp,tauTemp,sym%nop,sym%nop2,sym%symor) + + IF (ALLOCATED(sym%mrot)) THEN + DEALLOCATE(sym%mrot) + END IF + ALLOCATE(sym%mrot(3,3,sym%nop)) + IF (ALLOCATED(sym%tau)) THEN + DEALLOCATE(sym%tau) + END IF + ALLOCATE(sym%tau(3,sym%nop)) + + sym%invs = .FALSE. + sym%zrfs = .FALSE. + + DO k = 1, sym%nop + absSum = 0 + DO i = 1, 3 + DO j = 1, 3 + sym%mrot(j,i,k) = mrotTemp(j,i,k) + absSum = absSum + ABS(sym%mrot(j,i,k)) + END DO + sym%tau(i,k) = tauTemp(i,k) + END DO + IF (absSum.EQ.3) THEN + IF (ALL(sym%tau(:,k).EQ.0.0)) THEN + IF ((sym%mrot(1,1,k).EQ.-1).AND.(sym%mrot(2,2,k).EQ.-1).AND.(sym%mrot(3,3,k).EQ.-1)) sym%invs = .TRUE. + IF ((sym%mrot(1,1,k).EQ.1).AND.(sym%mrot(2,2,k).EQ.1).AND.(sym%mrot(3,3,k).EQ.-1)) sym%zrfs = .TRUE. + END IF + END IF + END DO + + sym%invs2 = sym%invs.AND.sym%zrfs + END IF + + xPathA = '/fleurInput/cell/symmetryOperations' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + sym%symSpecType = 3 + symmetryDef = 3 + + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/symOp') + sym%nop = numberNodes + + IF (ALLOCATED(sym%mrot)) DEALLOCATE(sym%mrot) + ALLOCATE(sym%mrot(3,3,sym%nop)) + IF (ALLOCATED(sym%tau)) DEALLOCATE(sym%tau) + ALLOCATE(sym%tau(3,sym%nop)) + sym%symor = .TRUE. + DO i = 1, sym%nop + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-1' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))))) + READ(valueString,*) sym%mrot(1,1,i), sym%mrot(1,2,i), sym%mrot(1,3,i), sym%tau(1,i) + + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-2' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))))) + READ(valueString,*) sym%mrot(2,1,i), sym%mrot(2,2,i), sym%mrot(2,3,i), sym%tau(2,i) + + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-3' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))))) + READ(valueString,*) sym%mrot(3,1,i), sym%mrot(3,2,i), sym%mrot(3,3,i), sym%tau(3,i) + + IF ((sym%tau(1,i)**2 + sym%tau(2,i)**2 + sym%tau(3,i)**2).GT.1.e-8) THEN + sym%symor = .FALSE. + END IF + DO j = 1,3 + IF (ABS(sym%tau(j,i)-0.33333) < 0.00001) THEN + sym%tau(j,i) = 1./3. + ENDIF + IF (ABS(sym%tau(j,i)+0.33333) < 0.00001) THEN + sym%tau(j,i) = -1./3. + ENDIF + IF (ABS(sym%tau(j,i)-0.66667) < 0.00001) THEN + sym%tau(j,i) = 2./3. + ENDIF + IF (ABS(sym%tau(j,i)+0.66667) < 0.00001) THEN + sym%tau(j,i) = -2./3. + ENDIF + ENDDO + END DO + END IF + + ! Calculate missing symmetry and cell properties and check consistency of parameters. + + ! Construction of amat + SELECT CASE (latticeDef) + CASE (1) + IF (cell%latnam.EQ.'squ') THEN + a2(2) = a1(1) + ELSE IF (cell%latnam.EQ.'c-b') THEN + aTemp = a1(1) + a1(1) = aTemp*0.5*sqrt(2.0) + a1(2) = -aTemp*0.5 + a2(1) = aTemp*0.5*sqrt(2.0) + a2(2) = aTemp*0.5 + ELSE IF (cell%latnam.EQ.'hex') THEN + aTemp = 0.5*a1(1) + a1(1) = aTemp*sqrt(3.0) + a1(2) = -aTemp + a2(1) = a1(1) + a2(2) = aTemp + ELSE IF (cell%latnam.EQ.'hx3') THEN + aTemp = 0.5*a1(1) + a1(1) = aTemp + a1(2) = -aTemp*sqrt(3.0) + a2(1) = a1(1) + a2(2) = -a1(2) + ELSE IF (cell%latnam.EQ.'fcc') THEN + aTemp = a1(1) + a1(1) = 0.0 ; a1(2) = 0.5*aTemp ; a1(3) = 0.5*aTemp + a2(1) = 0.5*aTemp ; a2(2) = 0.0 ; a2(3) = 0.5*aTemp + a3(1) = 0.5*aTemp ; a3(2) = 0.5*aTemp ; a3(3) = 0.0 + ELSE IF (cell%latnam.EQ.'bcc') THEN + aTemp = a1(1) + a1(1) =-0.5*aTemp ; a1(2) = 0.5*aTemp ; a1(3) = 0.5*aTemp + a2(1) = 0.5*aTemp ; a2(2) =-0.5*aTemp ; a2(3) = 0.5*aTemp + a3(1) = 0.5*aTemp ; a3(2) = 0.5*aTemp ; a3(3) =-0.5*aTemp + ELSE + CALL juDFT_error("latnam is incompatible to parametrization of lattice (1)", calledby ="r_inpXML") + END IF + CASE (2) + IF ((cell%latnam.EQ.'c-r').OR.(cell%latnam.EQ.'p-r')) THEN + IF (cell%latnam.EQ.'c-r') THEN + a1(2) = -a2(2) + a2(1) = a1(1) + END IF + ELSE + CALL juDFT_error("latnam is incompatible to parametrization of lattice (2)", calledby ="r_inpXML") + END IF + CASE (3) + IF (.NOT.(cell%latnam.EQ.'obl')) THEN + CALL juDFT_error("latnam is incompatible to parametrization of lattice (3)", calledby ="r_inpXML") + END IF + CASE (4) + IF (.NOT.(cell%latnam.EQ.'any')) THEN + CALL juDFT_error("latnam is incompatible to parametrization of lattice (4)", calledby ="r_inpXML") + END IF + CASE DEFAULT + CALL juDFT_error("Illegal lattice definition", calledby ="r_inpXML") + END SELECT + + IF (latticeScale.EQ.0.0) latticeScale = 1.0 + IF (.NOT.input%film) vacuum%dvac = a3(3) + vacuum%dvac = latticeScale*vacuum%dvac + dtild = latticeScale*dtild + + cell%amat(:,1) = a1(:) * latticeScale + cell%amat(:,2) = a2(:) * latticeScale + cell%amat(:,3) = a3(:) * latticeScale + + CALL inv3(cell%amat,cell%bmat,cell%omtil) + cell%bmat(:,:) = tpi_const*cell%bmat(:,:) + cell%omtil = abs(cell%omtil) + + IF (input%film.AND..NOT.oneD%odd%d1) THEN + cell%vol = (cell%omtil/cell%amat(3,3))*vacuum%dvac + cell%area = cell%omtil/cell%amat(3,3) + !-odim + ELSE IF (oneD%odd%d1) THEN + cell%area = tpi_const*cell%amat(3,3) + cell%vol = pi_const*(vacuum%dvac**2)*cell%amat(3,3)/4.0 + !+odim + ELSE + cell%vol = cell%omtil + cell%area = cell%amat(1,1)*cell%amat(2,2)-cell%amat(1,2)*cell%amat(2,1) + IF (cell%area.lt.1.0e-7) THEN + IF (cell%latnam.EQ.'any') THEN + cell%area = 1. + ELSE + CALL juDFT_error("area = 0",calledby ="r_inpXML") + END IF + END IF + END IF + + ! Construction of missing symmetry information + IF ((symmetryDef.EQ.2).OR.(symmetryDef.EQ.3)) THEN + nop48 = 48 + ALLOCATE (invOps(sym%nop),multtab(sym%nop,sym%nop),optype(nop48)) + CALL check_close(sym%nop,sym%mrot,sym%tau,& + & multtab,invOps,optype) + + CALL symproperties(nop48,optype,input%film,sym%nop,multtab,cell%amat,& + & sym%symor,sym%mrot,sym%tau,& + & invSym,sym%invs,sym%zrfs,sym%invs2,sym%nop,sym%nop2) + DEALLOCATE(invOps,multtab,optype) + IF (.not.input%film) sym%nop2=sym%nop + IF (input%film) THEN + DO n = 1, sym%nop + DO i = 1, 3 + IF (ABS(sym%tau(i,n)) > 0.00001) THEN + CALL juDFT_error("nonsymmorphic symmetries not yet implemented for films!",calledby ="r_inpXML") + ENDIF + END DO + END DO + END IF + END IF + sym%invs2 = sym%invs.AND.sym%zrfs + + ALLOCATE (sym%invarop(atoms%nat,sym%nop),sym%invarind(atoms%nat)) + ALLOCATE (sym%multab(sym%nop,sym%nop),sym%invtab(sym%nop)) + ALLOCATE (sym%invsatnr(atoms%nat),sym%d_wgn(-3:3,-3:3,3,sym%nop)) + + !some settings for film calculations + vacuum%nmzd = 250 + vacuum%nmzxyd = 100 + vacuum%nvac = 2 + IF (sym%zrfs.OR.sym%invs) vacuum%nvac = 1 + IF (oneD%odd%d1) vacuum%nvac = 1 + cell%z1 = vacuum%dvac/2 + vacuum%nmz = vacuum%nmzd + vacuum%delz = 25.0/vacuum%nmz + IF (oneD%odd%d1) vacuum%delz = 20.0 / vacuum%nmz + IF (vacuum%nmz.GT.vacuum%nmzd) CALL juDFT_error("nmzd",calledby ="inped") + vacuum%nmzxy = vacuum%nmzxyd + IF (vacuum%nmzxy.GT.vacuum%nmzxyd) CALL juDFT_error("nmzxyd",calledby ="inped") !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! End of cell section @@ -1111,463 +1104,462 @@ SUBROUTINE r_inpXML(& !!! Start of XC functional section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Read in libxc parameters if present - if( xmlGetNumberOfNodes('/fleurInput/xcFunctional/LibXCID') == 1 & - .and. xmlGetNumberOfNodes('/fleurInput/xcFunctional/LibXCName') == 1) then - call judft_error("LibXC is given both by Name and ID and is therefore overdetermined", calledby="r_inpXML") - endif - - IF (xmlGetNumberOfNodes('/fleurInput/xcFunctional/LibXCID') == 1) THEN - id_x=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/LibXCID/@exchange')) - id_c=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/LibXCID/@correlation')) - ELSEIF (xmlGetNumberOfNodes('/fleurInput/xcFunctional/LibXCName') == 1) THEN - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/LibXCName/@exchange'))))) + !Read in libxc parameters if present + if( xmlGetNumberOfNodes('/fleurInput/xcFunctional/LibXCID') == 1 & + .and. xmlGetNumberOfNodes('/fleurInput/xcFunctional/LibXCName') == 1) then + call judft_error("LibXC is given both by Name and ID and is therefore overdetermined", calledby="r_inpXML") + endif + + IF (xmlGetNumberOfNodes('/fleurInput/xcFunctional/LibXCID') == 1) THEN + id_x=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/LibXCID/@exchange')) + id_c=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/LibXCID/@correlation')) + ELSEIF (xmlGetNumberOfNodes('/fleurInput/xcFunctional/LibXCName') == 1) THEN + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/LibXCName/@exchange'))))) #ifdef CPP_LIBXC - id_x = xc_f03_functional_get_number(TRIM(valueString)) - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/LibXCName/@correlation'))))) - id_c = xc_f03_functional_get_number(TRIM(valueString)) + id_x = xc_f03_functional_get_number(TRIM(valueString)) + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/LibXCName/@correlation'))))) + id_c = xc_f03_functional_get_number(TRIM(valueString)) #else - CALL judft_error("To use libxc functionals you have to compile with libXC support") -#endif - ELSE - id_x=0;id_c=0 - ENDIF - - write (*,*) "id_x = ", id_x, "id_c = ", id_c - - ! Read in xc functional parameters - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/@name'))))) - namex(1:4) = valueString(1:4) - l_relcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/@relativisticCorrections')) - - relcor = 'non-relativi' - IF (l_relcor) THEN - relcor = 'relativistic' - END IF - - !now initialize the xcpot variable - CALL setXCParameters(atoms,valueString,l_relcor,input%jspins,id_x,id_c,xcpot) - - xPathA = '/fleurInput/calculationSetup/cutoffs/@GmaxXC' - numberNodes = xmlGetNumberOfNodes(xPathA) - xcpot%gmaxxc = stars%gmax - IF(numberNodes.EQ.1) THEN - xcpot%gmaxxc = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) - END IF - hybrid%l_hybrid=xcpot%is_hybrid() - - ALLOCATE(hybrid%lcutm1(atoms%ntype),hybrid%lcutwf(atoms%ntype),hybrid%select1(4,atoms%ntype)) - - obsolete%lwb=.FALSE. - IF (xcpot%is_gga()) THEN - obsolete%ndvgrd=6 - obsolete%chng=-0.1e-11 - END IF - - IF (xcpot%is_gga()) THEN - obsolete%ndvgrd = max(obsolete%ndvgrd,3) - END IF - - - hybrid%gcutm1 = input%rkmax - 0.5 - hybrid%tolerance1 = 1.0e-4 - hybrid%ewaldlambda = 3 - hybrid%lexp = 16 - hybrid%bands1 = dimension%neigd - - numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/prodBasis') - IF (numberNodes==0) THEN - IF (hybrid%l_hybrid) CALL judft_error("Mixed product basis input missing in inp.xml") - ELSE - hybrid%gcutm1=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@gcutm')) - hybrid%tolerance1=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@tolerance')) - hybrid%ewaldlambda=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@ewaldlambda')) - hybrid%lexp=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@lexp')) - hybrid%bands1=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@bands')) - ENDIF + CALL judft_error("To use libxc functionals you have to compile with libXC support") +#endif + ELSE + id_x=0; id_c=0 + ENDIF + + write (*,*) "id_x = ", id_x, "id_c = ", id_c + + ! Read in xc functional parameters + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/@name'))))) + namex(1:4) = valueString(1:4) + l_relcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/@relativisticCorrections')) + + relcor = 'non-relativi' + IF (l_relcor) THEN + relcor = 'relativistic' + END IF + + !now initialize the xcpot variable + CALL setXCParameters(atoms,valueString,l_relcor,input%jspins,id_x,id_c,xcpot) + + xPathA = '/fleurInput/calculationSetup/cutoffs/@GmaxXC' + numberNodes = xmlGetNumberOfNodes(xPathA) + xcpot%gmaxxc = stars%gmax + IF(numberNodes.EQ.1) THEN + xcpot%gmaxxc = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) + END IF + hybrid%l_hybrid=xcpot%is_hybrid() + + ALLOCATE(hybrid%lcutm1(atoms%ntype),hybrid%lcutwf(atoms%ntype),hybrid%select1(4,atoms%ntype)) + + obsolete%lwb=.FALSE. + IF (xcpot%is_gga()) THEN + obsolete%ndvgrd=6 + obsolete%chng=-0.1e-11 + END IF + + IF (xcpot%is_gga()) THEN + obsolete%ndvgrd = max(obsolete%ndvgrd,3) + END IF + + hybrid%gcutm1 = input%rkmax - 0.5 + hybrid%tolerance1 = 1.0e-4 + hybrid%ewaldlambda = 3 + hybrid%lexp = 16 + hybrid%bands1 = dimension%neigd + + numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/prodBasis') + IF (numberNodes==0) THEN + IF (hybrid%l_hybrid) CALL judft_error("Mixed product basis input missing in inp.xml") + ELSE + hybrid%gcutm1=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@gcutm')) + hybrid%tolerance1=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@tolerance')) + hybrid%ewaldlambda=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@ewaldlambda')) + hybrid%lexp=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@lexp')) + hybrid%bands1=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@bands')) + ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! End of XC functional section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Start of species section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ALLOCATE (speciesNLO(numSpecies)) - ALLOCATE(atoms%speciesName(numSpecies)) - - atoms%numStatesProvided = 0 - atoms%lapw_l(:) = -1 - atoms%n_u = 0 - - DEALLOCATE(noel) - ALLOCATE(noel(atoms%ntype)) - - DO iSpecies = 1, numSpecies - ! Attributes of species - WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']' - atoms%speciesName(iSpecies) = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@name'))) - atomicNumber = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomicNumber')) - coreStates = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@coreStates')) - magMom = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@magMom')) - flipSpin = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@flipSpin')) - - ! Attributes of mtSphere element of species - radius = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@radius')) - gridPoints = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@gridPoints')) - logIncrement = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@logIncrement')) - - ! Attributes of atomicCutoffs element of species - lmax = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmax')) - lnonsphr = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lnonsphr')) - lmaxAPW = -1 - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmaxAPW') - IF (numberNodes.EQ.1) THEN - lmaxAPW = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmaxAPW')) - END IF - - numU = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/ldaU') - IF (numU.GT.4) CALL juDFT_error("Too many U parameters provided for a certain species (maximum is 4).",calledby ="r_inpXML") - ldau_l = -1 - ldau_u = 0.0 - ldau_j = 0.0 - l_amf = .FALSE. - DO i = 1, numU - WRITE(xPathB,*) i - ldau_l(i) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@l')) - ldau_u(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@U')) - ldau_j(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@J')) - l_amf(i) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@l_amf')) - END DO - - speciesNLO(iSpecies) = 0 - WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/lo' - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) - DO iLO = 1, numberNodes - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@l' - WRITE(xPathC,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@n' - lString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) - nString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))) - CALL getIntegerSequenceFromString(TRIM(ADJUSTL(lString)), lNumbers, lNumCount) - CALL getIntegerSequenceFromString(TRIM(ADJUSTL(nString)), nNumbers, nNumCount) - IF(lNumCount.NE.nNumCount) THEN - STOP 'Error in LO input: l quantum number count does not equal n quantum number count' - END IF - speciesNLO(iSpecies) = speciesNLO(iSpecies) + lNumCount - DEALLOCATE (lNumbers, nNumbers) - END DO - - DO iType = 1, atoms%ntype - WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']/@species' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) - IF(TRIM(ADJUSTL(atoms%speciesName(iSpecies))).EQ.TRIM(ADJUSTL(valueString))) THEN - atoms%nz(iType) = atomicNumber - IF (atoms%nz(iType).EQ.0) THEN - WRITE(*,*) 'Note: Replacing atomic number 0 by 1.0e-10 on atom type ', iType - atoms%zatom(iType) = 1.0e-10 - END IF - noel(iType) = namat_const(atoms%nz(iType)) - atoms%zatom(iType) = atoms%nz(iType) - atoms%rmt(iType) = radius - atoms%jri(iType) = gridPoints - atoms%dx(iType) = logIncrement - atoms%lmax(iType) = lmax - atoms%nlo(iType) = speciesNLO(iSpecies) - atoms%ncst(iType) = coreStates - atoms%lnonsph(iType) = lnonsphr - atoms%lapw_l(iType) = lmaxAPW - IF (flipSpin) THEN - atoms%nflip(iType) = 1 - ELSE - atoms%nflip(iType) = 0 - ENDIF - atoms%bmu(iType) = magMom - DO i = 1, numU - atoms%n_u = atoms%n_u + 1 - atoms%lda_u(atoms%n_u)%l = ldau_l(i) - atoms%lda_u(atoms%n_u)%u = ldau_u(i) - atoms%lda_u(atoms%n_u)%j = ldau_j(i) - atoms%lda_u(atoms%n_u)%l_amf = l_amf(i) - atoms%lda_u(atoms%n_u)%atomType = iType - END DO - atomTypeSpecies(iType) = iSpecies - IF(speciesRepAtomType(iSpecies).EQ.-1) speciesRepAtomType(iSpecies) = iType - END IF - END DO - END DO - - atoms%lmaxd = maxval(atoms%lmax(:)) - atoms%llod = 0 - atoms%nlod = 0 - DO iType = 1, atoms%ntype - atoms%nlod = max(atoms%nlod,atoms%nlo(iType)) - END DO - atoms%nlod = max(atoms%nlod,2) ! for chkmt - ALLOCATE(atoms%llo(atoms%nlod,atoms%ntype));atoms%llo=-1 - ALLOCATE(atoms%ulo_der(atoms%nlod,atoms%ntype)) - ALLOCATE(atoms%l_dulo(atoms%nlod,atoms%ntype)) ! For what is this? - - dimension%nstd = 29 - - ALLOCATE(atoms%coreStateOccs(dimension%nstd,2,atoms%ntype));atoms%coreStateOccs=0.0 - ALLOCATE(atoms%coreStateNprnc(dimension%nstd,atoms%ntype)) - ALLOCATE(atoms%coreStateKappa(dimension%nstd,atoms%ntype)) - - CALL enpara%init(atoms,input%jspins) - enpara%evac0(:,:) = evac0Temp(:,:) - - DO iSpecies = 1, numSpecies - ALLOCATE(speciesLLO(speciesNLO(iSpecies))) - ALLOCATE(speciesLOeParams(speciesNLO(iSpecies))) - ALLOCATE(speciesLOEDeriv(speciesNLO(iSpecies))) - - ! Attributes of energyParameters element of species - WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']' - speciesEParams(0) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@s')) - speciesEParams(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@p')) - speciesEParams(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@d')) - speciesEParams(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@f')) - - ! Parameters for hybrid functionals - IF (hybrid%l_hybrid) THEN - WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/prodBasis' - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) - IF (numberNodes.NE.1) CALL judft_error("Parameters for mixed basis are missing for some specified") - lcutm =evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lcutm')) - lcutwf=evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lcutwf')) - xPathA=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@select') - hybSelect(1) = NINT(evaluateFirst(xPathA)) - hybSelect(2) = NINT(evaluateFirst(xPathA)) - hybSelect(3) = NINT(evaluateFirst(xPathA)) - hybSelect(4) = NINT(evaluateFirst(xPathA)) - ENDIF - - ! Special switches for species - ldaspecies=.FALSE. - socscalespecies=1.0 - WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/special' - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) - IF (numberNodes==1) THEN - ldaSpecies = evaluateFirstBoolOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lda')))) - socscaleSpecies = evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@socscale')))) - ENDIF - ! Explicitely provided core configurations - - coreConfigPresent = .FALSE. - providedCoreStates = 0 - providedStates = 0 - coreStateOccs = 0.0 - speciesXMLElectronStates = noState_const - speciesXMLCoreOccs = -1.0 - speciesXMLPrintCoreStates = .FALSE. - WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/electronConfig' - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) - IF (numberNodes.EQ.1) THEN - coreConfigPresent = .TRUE. - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/coreConfig') - token = popFirstStringToken(valueString) - DO WHILE (token.NE.' ') - IF (token(1:1).EQ.'[') THEN - DO i = 1, 6 - IF (TRIM(ADJUSTL(token)).EQ.nobleGasConfigList_const(i)) THEN - IF (providedCoreStates+nobleGasNumStatesList_const(i).GT.29) THEN - STOP 'Error: Too many core states provided in xml input file!' - END IF - DO j = providedCoreStates+1, providedCoreStates+nobleGasNumStatesList_const(i) - coreStateOccs(j-providedCoreStates,:) = coreStateNumElecsList_const(j) - coreStateNprnc(j-providedCoreStates) = coreStateNprncList_const(j) - coreStateKappa(j-providedCoreStates) = coreStateKappaList_const(j) - speciesXMLElectronStates(j) = coreState_const - END DO - providedCoreStates = providedCoreStates + nobleGasNumStatesList_const(i) - END IF - END DO - ELSE - DO i = 1, 29 - IF (TRIM(ADJUSTL(token)).EQ.coreStateList_const(i)) THEN - providedCoreStates = providedCoreStates + 1 - IF (providedCoreStates.GT.29) THEN - STOP 'Error: Too many core states provided in xml input file!' - END IF - coreStateOccs(providedCoreStates,:) = coreStateNumElecsList_const(i) - coreStateNprnc(providedCoreStates) = coreStateNprncList_const(i) - coreStateKappa(providedCoreStates) = coreStateKappaList_const(i) - speciesXMLElectronStates(i) = coreState_const - END IF - END DO - END IF - token = popFirstStringToken(valueString) - END DO - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/valenceConfig') - providedStates = providedCoreStates - IF(numberNodes.EQ.1) THEN - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/valenceConfig') - token = popFirstStringToken(valueString) - DO WHILE (token.NE.' ') - DO i = 1, 29 - IF (TRIM(ADJUSTL(token)).EQ.coreStateList_const(i)) THEN - providedStates = providedStates + 1 - IF (providedStates.GT.29) THEN - STOP 'Error: Too many valence states provided in xml input file!' - END IF - coreStateOccs(providedStates,:) = coreStateNumElecsList_const(i) - coreStateNprnc(providedStates) = coreStateNprncList_const(i) - coreStateKappa(providedStates) = coreStateKappaList_const(i) - speciesXMLElectronStates(i) = valenceState_const - END IF - END DO - token = popFirstStringToken(valueString) - END DO - END IF - END IF - - ! Explicitely provided core occupations - - WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/electronConfig/stateOccupation' - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) - IF (numberNodes.GE.1) THEN - IF (.NOT.coreConfigPresent) THEN - WRITE(*,*) 'Note: This just has to be implemented:' - STOP 'Error: Core occupation given while core config not set!' - END IF - DO i = 1, numberNodes - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',i,']' - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@state') - nprncTemp = 0 - kappaTemp = 0 - DO j = 1, 29 - IF (TRIM(ADJUSTL(valueString)).EQ.coreStateList_const(j)) THEN - nprncTemp = coreStateNprncList_const(j) - kappaTemp = coreStateKappaList_const(j) - speciesXMLPrintCoreStates(j) = .TRUE. - speciesXMLCoreOccs(1,j) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinUp')) - speciesXMLCoreOccs(2,j) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinDown')) - END IF - END DO - DO j = 1, providedStates - IF ((nprncTemp.EQ.coreStateNprnc(j)).AND.(kappaTemp.EQ.coreStateKappa(j))) THEN - coreStateOccs(j,1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinUp')) - coreStateOccs(j,2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinDown')) - END IF - END DO - END DO - END IF - - ! local orbitals - - WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/lo' - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) - iLLO = 1 - DO iLO = 1, numberNodes - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@l' - WRITE(xPathC,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@n' - WRITE(xPathD,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@type' - WRITE(xPathE,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@eDeriv' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathD))))) - lString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) - nString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))) - CALL getIntegerSequenceFromString(TRIM(ADJUSTL(lString)), lNumbers, lNumCount) - CALL getIntegerSequenceFromString(TRIM(ADJUSTL(nString)), nNumbers, nNumCount) - IF(lNumCount.NE.nNumCount) THEN - STOP 'Error in LO input: l quantum number count does not equal n quantum number count' - END IF - loEDeriv = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathE)))) - DO i = 1, lNumCount - speciesLLO(iLLO) = lNumbers(i) - speciesLOeParams(iLLO) = nNumbers(i) - IF(TRIM(ADJUSTL(valueString)).EQ.'HELO') THEN - speciesLOeParams(iLLO) = -speciesLOeParams(iLLO) - END IF - speciesLOEDeriv(iLLO) = loEDeriv - iLLO = iLLO + 1 - END DO - DEALLOCATE (lNumbers, nNumbers) - END DO - - ! sort LOs according to l quantum number - - ALLOCATE (loOrderList(speciesNLO(iSpecies)),speciesLLOReal(speciesNLO(iSpecies))) - DO iLLO = 1, speciesNLO(iSpecies) - speciesLLOReal(iLLO) = speciesLLO(iLLO) - END DO - CALL sort(speciesNLO(iSpecies),speciesLLOReal,loOrderList) - DEALLOCATE(speciesLLOReal) - - ! apply species parameters to atom groups - - DO iType = 1, atoms%ntype - WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']/@species' - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) - IF(TRIM(ADJUSTL(atoms%speciesName(iSpecies))).EQ.TRIM(ADJUSTL(valueString))) THEN - atoms%numStatesProvided(iType) = providedStates - IF (coreConfigPresent) THEN - IF (providedCoreStates.NE.atoms%ncst(iType)) THEN - WRITE(6,*) " providedCoreStates:",providedCoreStates - WRITE(6,*) "atoms%ncst(iType) :",atoms%ncst(iType) - STOP 'Wrong number of core states provided!' - END IF - DO k = 1, providedStates !atoms%ncst(iType) - atoms%coreStateOccs(k,1,iType) = coreStateOccs(k,1) - atoms%coreStateOccs(k,2,iType) = coreStateOccs(k,2) - atoms%coreStateNprnc(k,iType) = coreStateNprnc(k) - atoms%coreStateKappa(k,iType) = coreStateKappa(k) - xmlElectronStates(k,iType) = speciesXMLElectronStates(k) - xmlPrintCoreStates(k,iType) = speciesXMLPrintCoreStates(k) - xmlCoreOccs (1,k,iType) = speciesXMLCoreOccs(1,k) - xmlCoreOccs (2,k,iType) = speciesXMLCoreOccs(2,k) - END DO - END IF - DO iLLO = 1, speciesNLO(iSpecies) - atoms%llo(iLLO,iType) = speciesLLO(loOrderList(iLLO)) - atoms%ulo_der(iLLO,iType) = speciesLOEDeriv(loOrderList(iLLO)) - atoms%llod = max(abs(atoms%llo(iLLO,iType)),atoms%llod) - DO jsp = 1, input%jspins - enpara%ello0(iLLO,iType,jsp) = speciesLOeParams(loOrderList(iLLO)) - IF (enpara%ello0(iLLO,iType,jsp)==NINT(enpara%ello0(iLLO,iType,jsp))) THEN - enpara%qn_ello(iLLO,iType,jsp)=NINT(enpara%ello0(iLLO,iType,jsp)) - enpara%ello0(iLLO,iType,jsp)=0 - ELSE - enpara%qn_ello(iLLO,iType,jsp)=0 - ENDIF - enpara%skiplo(iType,jsp)=enpara%skiplo(iType,jsp)+(2*atoms%llo(iLLO,itype)+1) - END DO - END DO - ! Energy parameters - DO jsp = 1, input%jspins - DO l = 0, 3 - enpara%el0(l,iType,jsp) = speciesEParams(l) - IF (enpara%el0(l,iType,jsp)==NINT(enpara%el0(l,iType,jsp))) THEN - enpara%qn_el(l,iType,jsp)=nint(enpara%el0(l,iType,jsp)) - enpara%el0(l,iType,jsp)=0 - ELSE - enpara%qn_el(l,iType,jsp)=0 - ENDIF - END DO - DO l = 4,atoms%lmax(iType) - enpara%el0(l,iType,jsp) = enpara%el0(3,iType,jsp) - END DO - END DO - !Hybrid functional stuff - hybrid%lcutm1(iType) = 4 - hybrid%lcutwf(iType) = atoms%lmax(iType) - atoms%lmax(iType) / 10 - hybrid%select1(:,iType) = (/4, 0, 4, 2 /) - IF (hybrid%l_hybrid) THEN - hybrid%lcutm1(iType)=lcutm - hybrid%lcutwf(iType)=lcutwf - hybrid%select1(:,iType)=hybSelect - ENDIF - ! Explicit xc functional - SELECT TYPE(xcpot) - TYPE IS(t_xcpot_inbuild) - xcpot%lda_atom(iType)=ldaSpecies - END SELECT - noco%socscale(iType)=socscaleSpecies - END IF - END DO - DEALLOCATE(loOrderList) - DEALLOCATE(speciesLLO,speciesLOeParams,speciesLOEDeriv) - END DO + ALLOCATE (speciesNLO(numSpecies)) + ALLOCATE(atoms%speciesName(numSpecies)) + + atoms%numStatesProvided = 0 + atoms%lapw_l(:) = -1 + atoms%n_u = 0 + + DEALLOCATE(noel) + ALLOCATE(noel(atoms%ntype)) + + DO iSpecies = 1, numSpecies + ! Attributes of species + WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']' + atoms%speciesName(iSpecies) = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@name'))) + atomicNumber = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomicNumber')) + coreStates = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@coreStates')) + magMom = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@magMom')) + flipSpin = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@flipSpin')) + + ! Attributes of mtSphere element of species + radius = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@radius')) + gridPoints = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@gridPoints')) + logIncrement = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@logIncrement')) + + ! Attributes of atomicCutoffs element of species + lmax = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmax')) + lnonsphr = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lnonsphr')) + lmaxAPW = -1 + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmaxAPW') + IF (numberNodes.EQ.1) THEN + lmaxAPW = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmaxAPW')) + END IF + + numU = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/ldaU') + IF (numU.GT.4) CALL juDFT_error("Too many U parameters provided for a certain species (maximum is 4).",calledby ="r_inpXML") + ldau_l = -1 + ldau_u = 0.0 + ldau_j = 0.0 + l_amf = .FALSE. + DO i = 1, numU + WRITE(xPathB,*) i + ldau_l(i) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@l')) + ldau_u(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@U')) + ldau_j(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@J')) + l_amf(i) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@l_amf')) + END DO + + speciesNLO(iSpecies) = 0 + WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/lo' + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) + DO iLO = 1, numberNodes + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@l' + WRITE(xPathC,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@n' + lString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) + nString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))) + CALL getIntegerSequenceFromString(TRIM(ADJUSTL(lString)), lNumbers, lNumCount) + CALL getIntegerSequenceFromString(TRIM(ADJUSTL(nString)), nNumbers, nNumCount) + IF(lNumCount.NE.nNumCount) THEN + STOP 'Error in LO input: l quantum number count does not equal n quantum number count' + END IF + speciesNLO(iSpecies) = speciesNLO(iSpecies) + lNumCount + DEALLOCATE (lNumbers, nNumbers) + END DO + + DO iType = 1, atoms%ntype + WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']/@species' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) + IF(TRIM(ADJUSTL(atoms%speciesName(iSpecies))).EQ.TRIM(ADJUSTL(valueString))) THEN + atoms%nz(iType) = atomicNumber + IF (atoms%nz(iType).EQ.0) THEN + WRITE(*,*) 'Note: Replacing atomic number 0 by 1.0e-10 on atom type ', iType + atoms%zatom(iType) = 1.0e-10 + END IF + noel(iType) = namat_const(atoms%nz(iType)) + atoms%zatom(iType) = atoms%nz(iType) + atoms%rmt(iType) = radius + atoms%jri(iType) = gridPoints + atoms%dx(iType) = logIncrement + atoms%lmax(iType) = lmax + atoms%nlo(iType) = speciesNLO(iSpecies) + atoms%ncst(iType) = coreStates + atoms%lnonsph(iType) = lnonsphr + atoms%lapw_l(iType) = lmaxAPW + IF (flipSpin) THEN + atoms%nflip(iType) = 1 + ELSE + atoms%nflip(iType) = 0 + ENDIF + atoms%bmu(iType) = magMom + DO i = 1, numU + atoms%n_u = atoms%n_u + 1 + atoms%lda_u(atoms%n_u)%l = ldau_l(i) + atoms%lda_u(atoms%n_u)%u = ldau_u(i) + atoms%lda_u(atoms%n_u)%j = ldau_j(i) + atoms%lda_u(atoms%n_u)%l_amf = l_amf(i) + atoms%lda_u(atoms%n_u)%atomType = iType + END DO + atomTypeSpecies(iType) = iSpecies + IF(speciesRepAtomType(iSpecies).EQ.-1) speciesRepAtomType(iSpecies) = iType + END IF + END DO + END DO + + atoms%lmaxd = maxval(atoms%lmax(:)) + atoms%llod = 0 + atoms%nlod = 0 + DO iType = 1, atoms%ntype + atoms%nlod = max(atoms%nlod,atoms%nlo(iType)) + END DO + atoms%nlod = max(atoms%nlod,2) ! for chkmt + ALLOCATE(atoms%llo(atoms%nlod,atoms%ntype)); atoms%llo=-1 + ALLOCATE(atoms%ulo_der(atoms%nlod,atoms%ntype)) + ALLOCATE(atoms%l_dulo(atoms%nlod,atoms%ntype)) ! For what is this? + + dimension%nstd = 29 + + ALLOCATE(atoms%coreStateOccs(dimension%nstd,2,atoms%ntype)); atoms%coreStateOccs=0.0 + ALLOCATE(atoms%coreStateNprnc(dimension%nstd,atoms%ntype)) + ALLOCATE(atoms%coreStateKappa(dimension%nstd,atoms%ntype)) + + CALL enpara%init(atoms,input%jspins) + enpara%evac0(:,:) = evac0Temp(:,:) + + DO iSpecies = 1, numSpecies + ALLOCATE(speciesLLO(speciesNLO(iSpecies))) + ALLOCATE(speciesLOeParams(speciesNLO(iSpecies))) + ALLOCATE(speciesLOEDeriv(speciesNLO(iSpecies))) + + ! Attributes of energyParameters element of species + WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']' + speciesEParams(0) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@s')) + speciesEParams(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@p')) + speciesEParams(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@d')) + speciesEParams(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@f')) + + ! Parameters for hybrid functionals + IF (hybrid%l_hybrid) THEN + WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/prodBasis' + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) + IF (numberNodes.NE.1) CALL judft_error("Parameters for mixed basis are missing for some specified") + lcutm =evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lcutm')) + lcutwf=evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lcutwf')) + xPathA=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@select') + hybSelect(1) = NINT(evaluateFirst(xPathA)) + hybSelect(2) = NINT(evaluateFirst(xPathA)) + hybSelect(3) = NINT(evaluateFirst(xPathA)) + hybSelect(4) = NINT(evaluateFirst(xPathA)) + ENDIF + + ! Special switches for species + ldaspecies=.FALSE. + socscalespecies=1.0 + WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/special' + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) + IF (numberNodes==1) THEN + ldaSpecies = evaluateFirstBoolOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lda')))) + socscaleSpecies = evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@socscale')))) + ENDIF + ! Explicitely provided core configurations + + coreConfigPresent = .FALSE. + providedCoreStates = 0 + providedStates = 0 + coreStateOccs = 0.0 + speciesXMLElectronStates = noState_const + speciesXMLCoreOccs = -1.0 + speciesXMLPrintCoreStates = .FALSE. + WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/electronConfig' + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) + IF (numberNodes.EQ.1) THEN + coreConfigPresent = .TRUE. + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/coreConfig') + token = popFirstStringToken(valueString) + DO WHILE (token.NE.' ') + IF (token(1:1).EQ.'[') THEN + DO i = 1, 6 + IF (TRIM(ADJUSTL(token)).EQ.nobleGasConfigList_const(i)) THEN + IF (providedCoreStates+nobleGasNumStatesList_const(i).GT.29) THEN + STOP 'Error: Too many core states provided in xml input file!' + END IF + DO j = providedCoreStates+1, providedCoreStates+nobleGasNumStatesList_const(i) + coreStateOccs(j-providedCoreStates,:) = coreStateNumElecsList_const(j) + coreStateNprnc(j-providedCoreStates) = coreStateNprncList_const(j) + coreStateKappa(j-providedCoreStates) = coreStateKappaList_const(j) + speciesXMLElectronStates(j) = coreState_const + END DO + providedCoreStates = providedCoreStates + nobleGasNumStatesList_const(i) + END IF + END DO + ELSE + DO i = 1, 29 + IF (TRIM(ADJUSTL(token)).EQ.coreStateList_const(i)) THEN + providedCoreStates = providedCoreStates + 1 + IF (providedCoreStates.GT.29) THEN + STOP 'Error: Too many core states provided in xml input file!' + END IF + coreStateOccs(providedCoreStates,:) = coreStateNumElecsList_const(i) + coreStateNprnc(providedCoreStates) = coreStateNprncList_const(i) + coreStateKappa(providedCoreStates) = coreStateKappaList_const(i) + speciesXMLElectronStates(i) = coreState_const + END IF + END DO + END IF + token = popFirstStringToken(valueString) + END DO + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/valenceConfig') + providedStates = providedCoreStates + IF(numberNodes.EQ.1) THEN + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/valenceConfig') + token = popFirstStringToken(valueString) + DO WHILE (token.NE.' ') + DO i = 1, 29 + IF (TRIM(ADJUSTL(token)).EQ.coreStateList_const(i)) THEN + providedStates = providedStates + 1 + IF (providedStates.GT.29) THEN + STOP 'Error: Too many valence states provided in xml input file!' + END IF + coreStateOccs(providedStates,:) = coreStateNumElecsList_const(i) + coreStateNprnc(providedStates) = coreStateNprncList_const(i) + coreStateKappa(providedStates) = coreStateKappaList_const(i) + speciesXMLElectronStates(i) = valenceState_const + END IF + END DO + token = popFirstStringToken(valueString) + END DO + END IF + END IF + + ! Explicitely provided core occupations + + WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/electronConfig/stateOccupation' + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) + IF (numberNodes.GE.1) THEN + IF (.NOT.coreConfigPresent) THEN + WRITE(*,*) 'Note: This just has to be implemented:' + STOP 'Error: Core occupation given while core config not set!' + END IF + DO i = 1, numberNodes + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',i,']' + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@state') + nprncTemp = 0 + kappaTemp = 0 + DO j = 1, 29 + IF (TRIM(ADJUSTL(valueString)).EQ.coreStateList_const(j)) THEN + nprncTemp = coreStateNprncList_const(j) + kappaTemp = coreStateKappaList_const(j) + speciesXMLPrintCoreStates(j) = .TRUE. + speciesXMLCoreOccs(1,j) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinUp')) + speciesXMLCoreOccs(2,j) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinDown')) + END IF + END DO + DO j = 1, providedStates + IF ((nprncTemp.EQ.coreStateNprnc(j)).AND.(kappaTemp.EQ.coreStateKappa(j))) THEN + coreStateOccs(j,1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinUp')) + coreStateOccs(j,2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinDown')) + END IF + END DO + END DO + END IF + + ! local orbitals + + WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/lo' + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))) + iLLO = 1 + DO iLO = 1, numberNodes + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@l' + WRITE(xPathC,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@n' + WRITE(xPathD,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@type' + WRITE(xPathE,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@eDeriv' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathD))))) + lString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) + nString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))) + CALL getIntegerSequenceFromString(TRIM(ADJUSTL(lString)), lNumbers, lNumCount) + CALL getIntegerSequenceFromString(TRIM(ADJUSTL(nString)), nNumbers, nNumCount) + IF(lNumCount.NE.nNumCount) THEN + STOP 'Error in LO input: l quantum number count does not equal n quantum number count' + END IF + loEDeriv = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathE)))) + DO i = 1, lNumCount + speciesLLO(iLLO) = lNumbers(i) + speciesLOeParams(iLLO) = nNumbers(i) + IF(TRIM(ADJUSTL(valueString)).EQ.'HELO') THEN + speciesLOeParams(iLLO) = -speciesLOeParams(iLLO) + END IF + speciesLOEDeriv(iLLO) = loEDeriv + iLLO = iLLO + 1 + END DO + DEALLOCATE (lNumbers, nNumbers) + END DO + + ! sort LOs according to l quantum number + + ALLOCATE (loOrderList(speciesNLO(iSpecies)),speciesLLOReal(speciesNLO(iSpecies))) + DO iLLO = 1, speciesNLO(iSpecies) + speciesLLOReal(iLLO) = speciesLLO(iLLO) + END DO + CALL sort(speciesNLO(iSpecies),speciesLLOReal,loOrderList) + DEALLOCATE(speciesLLOReal) + + ! apply species parameters to atom groups + + DO iType = 1, atoms%ntype + WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']/@species' + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) + IF(TRIM(ADJUSTL(atoms%speciesName(iSpecies))).EQ.TRIM(ADJUSTL(valueString))) THEN + atoms%numStatesProvided(iType) = providedStates + IF (coreConfigPresent) THEN + IF (providedCoreStates.NE.atoms%ncst(iType)) THEN + WRITE(6,*) " providedCoreStates:",providedCoreStates + WRITE(6,*) "atoms%ncst(iType) :",atoms%ncst(iType) + STOP 'Wrong number of core states provided!' + END IF + DO k = 1, providedStates !atoms%ncst(iType) + atoms%coreStateOccs(k,1,iType) = coreStateOccs(k,1) + atoms%coreStateOccs(k,2,iType) = coreStateOccs(k,2) + atoms%coreStateNprnc(k,iType) = coreStateNprnc(k) + atoms%coreStateKappa(k,iType) = coreStateKappa(k) + xmlElectronStates(k,iType) = speciesXMLElectronStates(k) + xmlPrintCoreStates(k,iType) = speciesXMLPrintCoreStates(k) + xmlCoreOccs (1,k,iType) = speciesXMLCoreOccs(1,k) + xmlCoreOccs (2,k,iType) = speciesXMLCoreOccs(2,k) + END DO + END IF + DO iLLO = 1, speciesNLO(iSpecies) + atoms%llo(iLLO,iType) = speciesLLO(loOrderList(iLLO)) + atoms%ulo_der(iLLO,iType) = speciesLOEDeriv(loOrderList(iLLO)) + atoms%llod = max(abs(atoms%llo(iLLO,iType)),atoms%llod) + DO jsp = 1, input%jspins + enpara%ello0(iLLO,iType,jsp) = speciesLOeParams(loOrderList(iLLO)) + IF (enpara%ello0(iLLO,iType,jsp)==NINT(enpara%ello0(iLLO,iType,jsp))) THEN + enpara%qn_ello(iLLO,iType,jsp)=NINT(enpara%ello0(iLLO,iType,jsp)) + enpara%ello0(iLLO,iType,jsp)=0 + ELSE + enpara%qn_ello(iLLO,iType,jsp)=0 + ENDIF + enpara%skiplo(iType,jsp)=enpara%skiplo(iType,jsp)+(2*atoms%llo(iLLO,itype)+1) + END DO + END DO + ! Energy parameters + DO jsp = 1, input%jspins + DO l = 0, 3 + enpara%el0(l,iType,jsp) = speciesEParams(l) + IF (enpara%el0(l,iType,jsp)==NINT(enpara%el0(l,iType,jsp))) THEN + enpara%qn_el(l,iType,jsp)=nint(enpara%el0(l,iType,jsp)) + enpara%el0(l,iType,jsp)=0 + ELSE + enpara%qn_el(l,iType,jsp)=0 + ENDIF + END DO + DO l = 4,atoms%lmax(iType) + enpara%el0(l,iType,jsp) = enpara%el0(3,iType,jsp) + END DO + END DO + !Hybrid functional stuff + hybrid%lcutm1(iType) = 4 + hybrid%lcutwf(iType) = atoms%lmax(iType) - atoms%lmax(iType) / 10 + hybrid%select1(:,iType) = (/4, 0, 4, 2 /) + IF (hybrid%l_hybrid) THEN + hybrid%lcutm1(iType)=lcutm + hybrid%lcutwf(iType)=lcutwf + hybrid%select1(:,iType)=hybSelect + ENDIF + ! Explicit xc functional + SELECT TYPE(xcpot) + TYPE IS(t_xcpot_inbuild) + xcpot%lda_atom(iType)=ldaSpecies + END SELECT + noco%socscale(iType)=socscaleSpecies + END IF + END DO + DEALLOCATE(loOrderList) + DEALLOCATE(speciesLLO,speciesLOeParams,speciesLOEDeriv) + END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! End of species section @@ -1577,106 +1569,106 @@ SUBROUTINE r_inpXML(& !!! Start of atomGroup section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - banddos%l_orb = .FALSE. - banddos%orbCompAtom = 0 - atoms%l_geo = .FALSE. - atoms%relax = 0 - na = 0 - firstAtomOfType = 1 - DO iType = 1, atoms%ntype - WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']' - - ! Read in force parameters - xPathB = TRIM(ADJUSTL(xPathA))//'/force' - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))) - IF (numberNodes.GE.1) THEN - atoms%l_geo(iType) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@calculate')) - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@relaxXYZ') - READ(valueString,'(3l1)') relaxX, relaxY, relaxZ - IF (relaxX) atoms%relax(1,iType) = 1 - IF (relaxY) atoms%relax(2,iType) = 1 - IF (relaxZ) atoms%relax(3,iType) = 1 - END IF - - ! Obtain number of equivalent atoms - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos') - numberNodes = numberNodes + xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/absPos') - numberNodes = numberNodes + xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/filmPos') - atoms%neq(iType) = numberNodes - - IF (iType.GE.2) THEN - firstAtomOfType = firstAtomOfType + atoms%neq(iType-1) - END IF - - ! Read in atom positions - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos') - DO i = 1, numberNodes - na = na + 1 - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'/relPos[',i,']' - IF(xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))//'/@label').NE.0) THEN - atoms%label(na) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@label') - ELSE - WRITE(atoms%label(na),'(i0)') na - END IF - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) - atoms%taual(1,na) = evaluatefirst(valueString) - atoms%taual(2,na) = evaluatefirst(valueString) - atoms%taual(3,na) = evaluatefirst(valueString) - atoms%pos(:,na) = matmul(cell%amat,atoms%taual(:,na)) - l_orbcomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@orbcomp')) - IF(l_orbcomp) THEN - IF(banddos%l_orb) THEN - CALL juDFT_error("Multiple orbcomp flags set.", calledby = "r_inpXML") - END IF - banddos%l_orb = .TRUE. - banddos%orbCompAtom = na - END IF - wannAtomList(na) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@wannier')) - END DO - - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/absPos') - DO i = 1, numberNodes - na = na + 1 - STOP 'absPos not yet implemented!' - END DO - - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/filmPos') - DO i = 1, numberNodes - na = na + 1 - WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'/filmPos[',i,']' - IF(xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))//'/@label').NE.0) THEN - atoms%label(na) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@label') - ELSE - WRITE(atoms%label(na),'(i0)') na - END IF - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) - atoms%taual(1,na) = evaluatefirst(valueString) - atoms%taual(2,na) = evaluatefirst(valueString) - atoms%taual(3,na) = evaluatefirst(valueString) / cell%amat(3,3) - atoms%pos(:,na) = matmul(cell%amat,atoms%taual(:,na)) - l_orbcomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@orbcomp')) - IF(l_orbcomp) THEN - IF(banddos%l_orb) THEN - CALL juDFT_error("Multiple orbcomp flags set.", calledby = "r_inpXML") - END IF - banddos%l_orb = .TRUE. - banddos%orbCompAtom = na - END IF - wannAtomList(na) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@wannier')) - END DO - - !Read in atom group specific noco parameters - xPathB = TRIM(ADJUSTL(xPathA))//'/nocoParams' - numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))) - IF (numberNodes.GE.1) THEN - noco%l_relax(iType) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@l_relax')) - noco%alphInit(iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@alpha')) - noco%alph(iType) = noco%alphInit(iType) - noco%beta(iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@beta')) - noco%b_con(1,iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@b_cons_x')) - noco%b_con(2,iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@b_cons_y')) - END IF - END DO + banddos%l_orb = .FALSE. + banddos%orbCompAtom = 0 + atoms%l_geo = .FALSE. + atoms%relax = 0 + na = 0 + firstAtomOfType = 1 + DO iType = 1, atoms%ntype + WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']' + + ! Read in force parameters + xPathB = TRIM(ADJUSTL(xPathA))//'/force' + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))) + IF (numberNodes.GE.1) THEN + atoms%l_geo(iType) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@calculate')) + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@relaxXYZ') + READ(valueString,'(3l1)') relaxX, relaxY, relaxZ + IF (relaxX) atoms%relax(1,iType) = 1 + IF (relaxY) atoms%relax(2,iType) = 1 + IF (relaxZ) atoms%relax(3,iType) = 1 + END IF + + ! Obtain number of equivalent atoms + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos') + numberNodes = numberNodes + xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/absPos') + numberNodes = numberNodes + xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/filmPos') + atoms%neq(iType) = numberNodes + + IF (iType.GE.2) THEN + firstAtomOfType = firstAtomOfType + atoms%neq(iType-1) + END IF + + ! Read in atom positions + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos') + DO i = 1, numberNodes + na = na + 1 + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'/relPos[',i,']' + IF(xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))//'/@label').NE.0) THEN + atoms%label(na) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@label') + ELSE + WRITE(atoms%label(na),'(i0)') na + END IF + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) + atoms%taual(1,na) = evaluatefirst(valueString) + atoms%taual(2,na) = evaluatefirst(valueString) + atoms%taual(3,na) = evaluatefirst(valueString) + atoms%pos(:,na) = matmul(cell%amat,atoms%taual(:,na)) + l_orbcomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@orbcomp')) + IF(l_orbcomp) THEN + IF(banddos%l_orb) THEN + CALL juDFT_error("Multiple orbcomp flags set.", calledby = "r_inpXML") + END IF + banddos%l_orb = .TRUE. + banddos%orbCompAtom = na + END IF + wannAtomList(na) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@wannier')) + END DO + + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/absPos') + DO i = 1, numberNodes + na = na + 1 + STOP 'absPos not yet implemented!' + END DO + + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/filmPos') + DO i = 1, numberNodes + na = na + 1 + WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'/filmPos[',i,']' + IF(xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))//'/@label').NE.0) THEN + atoms%label(na) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@label') + ELSE + WRITE(atoms%label(na),'(i0)') na + END IF + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) + atoms%taual(1,na) = evaluatefirst(valueString) + atoms%taual(2,na) = evaluatefirst(valueString) + atoms%taual(3,na) = evaluatefirst(valueString) / cell%amat(3,3) + atoms%pos(:,na) = matmul(cell%amat,atoms%taual(:,na)) + l_orbcomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@orbcomp')) + IF(l_orbcomp) THEN + IF(banddos%l_orb) THEN + CALL juDFT_error("Multiple orbcomp flags set.", calledby = "r_inpXML") + END IF + banddos%l_orb = .TRUE. + banddos%orbCompAtom = na + END IF + wannAtomList(na) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@wannier')) + END DO + + !Read in atom group specific noco parameters + xPathB = TRIM(ADJUSTL(xPathA))//'/nocoParams' + numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))) + IF (numberNodes.GE.1) THEN + noco%l_relax(iType) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@l_relax')) + noco%alphInit(iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@alpha')) + noco%alph(iType) = noco%alphInit(iType) + noco%beta(iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@beta')) + noco%b_con(1,iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@b_cons_x')) + noco%b_con(2,iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@b_cons_y')) + END IF + END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! End of atomGroup section @@ -1686,433 +1678,431 @@ SUBROUTINE r_inpXML(& !!! Start of force-theorem section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - xPathA = '/fleurInput/forceTheorem' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - !Magnetic anisotropy... - xPathA = '/fleurInput/forceTheorem/MAE' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - lString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta') - nString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi') - ALLOCATE(t_forcetheo_mae::forcetheo) - SELECT TYPE(forcetheo) - TYPE IS(t_forcetheo_mae) !this is ok, we just allocated the type... - CALL forcetheo%init(lString,nString) - END SELECT - ENDIF - !spin-spiral dispersion - xPathA = '/fleurInput/forceTheorem/spinSpiralDispersion' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - ALLOCATE(t_forcetheo_ssdisp::forcetheo) - SELECT TYPE(forcetheo) - TYPE IS(t_forcetheo_ssdisp) !this is ok, we just allocated the type... - CALL forcetheo%init(priv_read_q_list(xPathA)) - END SELECT - ENDIF - !dmi - xPathA = '/fleurInput/forceTheorem/DMI' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - lString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta') - nString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi') - ALLOCATE(t_forcetheo_dmi::forcetheo) - SELECT TYPE(forcetheo) - TYPE IS(t_forcetheo_dmi) !this is ok, we just allocated the type... - CALL forcetheo%init(priv_read_q_list(TRIM(ADJUSTL(xPathA))//'/qVectors'),lstring,nstring) - END SELECT - ENDIF - !jij - xPathA = '/fleurInput/forceTheorem/Jij' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - thetaj=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@thetaj')) - ALLOCATE(t_forcetheo_jij::forcetheo) - SELECT TYPE(forcetheo) - TYPE IS(t_forcetheo_jij) !this is ok, we just allocated the type... - CALL forcetheo%init(priv_read_q_list(TRIM(ADJUSTL(xPathA))//'/qVectors'),thetaj,atoms) - END SELECT - ENDIF - - ELSE - ALLOCATE(t_forcetheo::forcetheo) !default no forcetheorem type - ENDIF - + xPathA = '/fleurInput/forceTheorem' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + !Magnetic anisotropy... + xPathA = '/fleurInput/forceTheorem/MAE' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + lString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta') + nString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi') + ALLOCATE(t_forcetheo_mae::forcetheo) + SELECT TYPE(forcetheo) + TYPE IS(t_forcetheo_mae) !this is ok, we just allocated the type... + CALL forcetheo%init(lString,nString) + END SELECT + ENDIF + !spin-spiral dispersion + xPathA = '/fleurInput/forceTheorem/spinSpiralDispersion' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + ALLOCATE(t_forcetheo_ssdisp::forcetheo) + SELECT TYPE(forcetheo) + TYPE IS(t_forcetheo_ssdisp) !this is ok, we just allocated the type... + CALL forcetheo%init(priv_read_q_list(xPathA)) + END SELECT + ENDIF + !dmi + xPathA = '/fleurInput/forceTheorem/DMI' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + lString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta') + nString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi') + ALLOCATE(t_forcetheo_dmi::forcetheo) + SELECT TYPE(forcetheo) + TYPE IS(t_forcetheo_dmi) !this is ok, we just allocated the type... + CALL forcetheo%init(priv_read_q_list(TRIM(ADJUSTL(xPathA))//'/qVectors'),lstring,nstring) + END SELECT + ENDIF + !jij + xPathA = '/fleurInput/forceTheorem/Jij' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + thetaj=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@thetaj')) + ALLOCATE(t_forcetheo_jij::forcetheo) + SELECT TYPE(forcetheo) + TYPE IS(t_forcetheo_jij) !this is ok, we just allocated the type... + CALL forcetheo%init(priv_read_q_list(TRIM(ADJUSTL(xPathA))//'/qVectors'),thetaj,atoms) + END SELECT + ENDIF + + ELSE + ALLOCATE(t_forcetheo::forcetheo) !default no forcetheorem type + ENDIF + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! End of force-theorem section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Start of output section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - banddos%dos = .FALSE. - banddos%band = .FALSE. - banddos%vacdos = .FALSE. - sliceplot%slice = .FALSE. - input%l_coreSpec = .FALSE. - input%l_wann = .FALSE. - - input%vchk = .FALSE. - input%cdinf = .FALSE. - - sliceplot%iplot = .FALSE. - input%score = .FALSE. - sliceplot%plpot = .FALSE. - - input%eonly = .FALSE. - input%l_bmt = .FALSE. - - xPathA = '/fleurInput/output' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - - ! Read in general output switches - banddos%dos = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dos')) - banddos%band = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@band')) - banddos%vacdos = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vacdos')) - sliceplot%slice = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@slice')) - input%l_coreSpec = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@coreSpec')) - input%l_wann = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@wannier')) - banddos%l_mcd = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mcd')) - - ! Read in optional switches for checks - - xPathA = '/fleurInput/output/checks' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - input%vchk = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vchk')) - input%cdinf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@cdinf')) - - END IF - - ! Read in optional plotting parameters - - xPathA = '/fleurInput/output/plotting' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - sliceplot%iplot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@iplot')) - input%score = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@score')) - sliceplot%plpot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plplot')) - END IF - - ! Read in optional specialOutput switches - - xPathA = '/fleurInput/output/specialOutput' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - input%eonly = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eonly')) - input%l_bmt = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@bmt')) - END IF - - ! Read in optional densityOfStates output parameters - - xPathA = '/fleurInput/output/densityOfStates' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF ((banddos%dos).AND.(numberNodes.EQ.0)) THEN - CALL juDFT_error("dos is true but densityOfStates parameters are not set!", calledby = "r_inpXML") - END IF - - IF (numberNodes.EQ.1) THEN - banddos%ndir = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ndir')) - banddos%e2_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minEnergy')) - banddos%e1_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxEnergy')) - banddos%sig_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sigma')) - END IF - - ! Read in optional vacuumDOS parameters - - xPathA = '/fleurInput/output/vacuumDOS' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF ((banddos%vacdos).AND.(numberNodes.EQ.0)) THEN - CALL juDFT_error("vacdos is true but vacDOS parameters are not set!", calledby = "r_inpXML") - END IF - - vacuum%layers = 1 - input%integ = .FALSE. - vacuum%starcoeff = .FALSE. - vacuum%nstars = 0 - vacuum%locx = 0.0 - vacuum%locy = 0.0 - vacuum%nstm = 0 - vacuum%tworkf = 0.0 - IF ((banddos%vacdos).AND.(numberNodes.EQ.1)) THEN - vacuum%layers = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@layers')) - input%integ = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@integ')) - vacuum%starcoeff = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@star')) - vacuum%nstars = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nstars')) - vacuum%locx(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locx1')) - vacuum%locx(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locx2')) - vacuum%locy(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locy1')) - vacuum%locy(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locy2')) - vacuum%nstm = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nstm')) - vacuum%tworkf = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@tworkf')) - END IF - vacuum%layerd = vacuum%layers - ALLOCATE(vacuum%izlay(vacuum%layerd,2)) - - ! Read in optional chargeDensitySlicing parameters - - xPathA = '/fleurInput/output/chargeDensitySlicing' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF ((sliceplot%slice).AND.(numberNodes.EQ.0)) THEN - CALL juDFT_error("slice is true but chargeDensitySlicing parameters are not set!", calledby = "r_inpXML") - END IF - - IF (numberNodes.EQ.1) THEN - sliceplot%kk = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@numkpt')) - sliceplot%e1s = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minEigenval')) - sliceplot%e2s = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxEigenval')) - sliceplot%nnne = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nnne')) - input%pallst = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@pallst')) - END IF - - IF (banddos%band) THEN - banddos%dos=.TRUE. - banddos%ndir = -4 - WRITE(*,*) 'band="T" --> Overriding "dos" and "ndir"!' - ENDIF - - ! Read in optional core spectrum (EELS) input parameters - - xPathA = '/fleurInput/output/coreSpectrum' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF ((input%l_coreSpec).AND.(numberNodes.EQ.0)) THEN - CALL juDFT_error("coreSpec is true but coreSpectrum parameters are not set!", calledby = "r_inpXML") - END IF - - IF (numberNodes.EQ.1) THEN - coreSpecInput%verb = 0 - tempBool = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@verbose')) - IF(tempBool) coreSpecInput%verb = 1 - coreSpecInput%ek0 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eKin')) - coreSpecInput%atomType = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomType')) - coreSpecInput%lx = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lmax')) - coreSpecInput%edge = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@edgeType'))) - coreSpecInput%emn = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eMin')) - coreSpecInput%emx = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eMax')) - tempInt = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@numPoints')) - coreSpecInput%ein = (coreSpecInput%emx - coreSpecInput%emn) / (tempInt - 1.0) - xPathB = TRIM(ADJUSTL(xPathA))//'/edgeIndices' - xPathB = TRIM(ADJUSTL(xPathB))//'/text()' - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) - numTokens = countStringTokens(valueString) - coreSpecInput%edgeidx(:) = 0 - IF(numTokens.GT.SIZE(coreSpecInput%edgeidx)) THEN - CALL juDFT_error('More EELS edge indices provided than allowed.',calledby='r_inpXML') - END IF - DO i = 1, MAX(numTokens,SIZE(coreSpecInput%edgeidx)) - coreSpecInput%edgeidx(i) = evaluateFirstIntOnly(popFirstStringToken(valueString)) - END DO - END IF - - ! Read in optional Wannier functions parameters - - xPathA = '/fleurInput/output/wannier' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF ((input%l_wann).AND.(numberNodes.EQ.0)) THEN - CALL juDFT_error("wannier is true but Wannier parameters are not set!", calledby = "r_inpXML") - END IF - - IF (numberNodes.EQ.1) THEN - wann%l_ms = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ms')) - wann%l_sgwf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sgwf')) - wann%l_socgwf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@socgwf')) - wann%l_bs_comf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@bsComf')) - wann%l_atomlist = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomList')) - END IF - - xPathA = '/fleurInput/output/wannier/bandSelection' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - wann%l_byindex=.TRUE. - wann%band_min(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minSpinUp')) - wann%band_max(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxSpinUp')) - xPathA = '/fleurInput/output/wannier/bandSelection/@minSpinDown' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - wann%band_min(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))) - ELSE - wann%band_min(2) = wann%band_min(1) - END IF - xPathA = '/fleurInput/output/wannier/bandSelection/@maxSpinDown' - numberNodes = xmlGetNumberOfNodes(xPathA) - IF (numberNodes.EQ.1) THEN - wann%band_max(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))) - ELSE - wann%band_max(2) = wann%band_max(1) - END IF - wann%l_byindex = .TRUE. - IF(input%l_wann) THEN - IF (dimension%neigd.LT.MAX(wann%band_max(1),wann%band_max(2))) THEN - dimension%neigd = MAX(wann%band_max(1),wann%band_max(2)) - END IF - END IF - END IF - - xPathA = '/fleurInput/output/wannier/jobList' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - xPathA = '/fleurInput/output/wannier/jobList/text()' - - ! Note: At the moment only 255 characters for the text in this node. Maybe this is not enough. - valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))) - numTokens = countStringTokens(valueString) - ALLOCATE(wann%jobList(numTokens)) - DO i = 1, numTokens - wann%jobList(i) = popFirstStringToken(valueString) - END DO - END IF - - ! Read in optional magnetic circular dichroism parameters - xPathA = '/fleurInput/output/magneticCircularDichroism' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF ((banddos%l_mcd).AND.(numberNodes.EQ.0)) THEN - CALL juDFT_error("mcd is true but magneticCircularDichroism parameters are not set!", calledby = "r_inpXML") - END IF - - IF (numberNodes.EQ.1) THEN - banddos%e_mcd_lo = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@energyLo')) - banddos%e_mcd_up = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@energyUp')) - END IF - - ! Read in optional parameter for unfolding bandstructure of supercell - xPathA = '/fleurInput/output/unfoldingBand' - numberNodes = xmlGetNumberOfNodes(xPathA) - - IF (numberNodes.EQ.1) THEN - banddos%unfoldband = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@unfoldband')) - banddos%s_cell_x = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellX')) - banddos%s_cell_y = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellY')) - banddos%s_cell_z = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellZ')) - END IF - - END IF + banddos%dos = .FALSE. + banddos%band = .FALSE. + banddos%vacdos = .FALSE. + sliceplot%slice = .FALSE. + input%l_coreSpec = .FALSE. + input%l_wann = .FALSE. + + input%vchk = .FALSE. + input%cdinf = .FALSE. + + sliceplot%iplot = .FALSE. + input%score = .FALSE. + sliceplot%plpot = .FALSE. + + input%eonly = .FALSE. + input%l_bmt = .FALSE. + + xPathA = '/fleurInput/output' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + + ! Read in general output switches + banddos%dos = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dos')) + banddos%band = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@band')) + banddos%vacdos = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vacdos')) + sliceplot%slice = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@slice')) + input%l_coreSpec = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@coreSpec')) + input%l_wann = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@wannier')) + banddos%l_mcd = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mcd')) + + ! Read in optional switches for checks + + xPathA = '/fleurInput/output/checks' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + input%vchk = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vchk')) + input%cdinf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@cdinf')) + + END IF + + ! Read in optional plotting parameters + + xPathA = '/fleurInput/output/plotting' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + sliceplot%iplot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@iplot')) + input%score = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@score')) + sliceplot%plpot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plplot')) + END IF + + ! Read in optional specialOutput switches + + xPathA = '/fleurInput/output/specialOutput' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + input%eonly = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eonly')) + input%l_bmt = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@bmt')) + END IF + + ! Read in optional densityOfStates output parameters + + xPathA = '/fleurInput/output/densityOfStates' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF ((banddos%dos).AND.(numberNodes.EQ.0)) THEN + CALL juDFT_error("dos is true but densityOfStates parameters are not set!", calledby = "r_inpXML") + END IF + + IF (numberNodes.EQ.1) THEN + banddos%ndir = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ndir')) + banddos%e2_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minEnergy')) + banddos%e1_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxEnergy')) + banddos%sig_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sigma')) + END IF + + ! Read in optional vacuumDOS parameters + + xPathA = '/fleurInput/output/vacuumDOS' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF ((banddos%vacdos).AND.(numberNodes.EQ.0)) THEN + CALL juDFT_error("vacdos is true but vacDOS parameters are not set!", calledby = "r_inpXML") + END IF + + vacuum%layers = 1 + input%integ = .FALSE. + vacuum%starcoeff = .FALSE. + vacuum%nstars = 0 + vacuum%locx = 0.0 + vacuum%locy = 0.0 + vacuum%nstm = 0 + vacuum%tworkf = 0.0 + IF ((banddos%vacdos).AND.(numberNodes.EQ.1)) THEN + vacuum%layers = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@layers')) + input%integ = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@integ')) + vacuum%starcoeff = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@star')) + vacuum%nstars = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nstars')) + vacuum%locx(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locx1')) + vacuum%locx(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locx2')) + vacuum%locy(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locy1')) + vacuum%locy(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locy2')) + vacuum%nstm = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nstm')) + vacuum%tworkf = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@tworkf')) + END IF + vacuum%layerd = vacuum%layers + ALLOCATE(vacuum%izlay(vacuum%layerd,2)) + + ! Read in optional chargeDensitySlicing parameters + + xPathA = '/fleurInput/output/chargeDensitySlicing' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF ((sliceplot%slice).AND.(numberNodes.EQ.0)) THEN + CALL juDFT_error("slice is true but chargeDensitySlicing parameters are not set!", calledby = "r_inpXML") + END IF + + IF (numberNodes.EQ.1) THEN + sliceplot%kk = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@numkpt')) + sliceplot%e1s = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minEigenval')) + sliceplot%e2s = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxEigenval')) + sliceplot%nnne = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nnne')) + input%pallst = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@pallst')) + END IF + + IF (banddos%band) THEN + banddos%dos=.TRUE. + banddos%ndir = -4 + WRITE(*,*) 'band="T" --> Overriding "dos" and "ndir"!' + ENDIF + + ! Read in optional core spectrum (EELS) input parameters + + xPathA = '/fleurInput/output/coreSpectrum' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF ((input%l_coreSpec).AND.(numberNodes.EQ.0)) THEN + CALL juDFT_error("coreSpec is true but coreSpectrum parameters are not set!", calledby = "r_inpXML") + END IF + + IF (numberNodes.EQ.1) THEN + coreSpecInput%verb = 0 + tempBool = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@verbose')) + IF(tempBool) coreSpecInput%verb = 1 + coreSpecInput%ek0 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eKin')) + coreSpecInput%atomType = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomType')) + coreSpecInput%lx = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lmax')) + coreSpecInput%edge = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@edgeType'))) + coreSpecInput%emn = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eMin')) + coreSpecInput%emx = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eMax')) + tempInt = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@numPoints')) + coreSpecInput%ein = (coreSpecInput%emx - coreSpecInput%emn) / (tempInt - 1.0) + xPathB = TRIM(ADJUSTL(xPathA))//'/edgeIndices' + xPathB = TRIM(ADJUSTL(xPathB))//'/text()' + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))) + numTokens = countStringTokens(valueString) + coreSpecInput%edgeidx(:) = 0 + IF(numTokens.GT.SIZE(coreSpecInput%edgeidx)) THEN + CALL juDFT_error('More EELS edge indices provided than allowed.',calledby='r_inpXML') + END IF + DO i = 1, MAX(numTokens,SIZE(coreSpecInput%edgeidx)) + coreSpecInput%edgeidx(i) = evaluateFirstIntOnly(popFirstStringToken(valueString)) + END DO + END IF + + ! Read in optional Wannier functions parameters + + xPathA = '/fleurInput/output/wannier' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF ((input%l_wann).AND.(numberNodes.EQ.0)) THEN + CALL juDFT_error("wannier is true but Wannier parameters are not set!", calledby = "r_inpXML") + END IF + + IF (numberNodes.EQ.1) THEN + wann%l_ms = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ms')) + wann%l_sgwf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sgwf')) + wann%l_socgwf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@socgwf')) + wann%l_bs_comf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@bsComf')) + wann%l_atomlist = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomList')) + END IF + + xPathA = '/fleurInput/output/wannier/bandSelection' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + wann%l_byindex=.TRUE. + wann%band_min(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minSpinUp')) + wann%band_max(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxSpinUp')) + xPathA = '/fleurInput/output/wannier/bandSelection/@minSpinDown' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + wann%band_min(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))) + ELSE + wann%band_min(2) = wann%band_min(1) + END IF + xPathA = '/fleurInput/output/wannier/bandSelection/@maxSpinDown' + numberNodes = xmlGetNumberOfNodes(xPathA) + IF (numberNodes.EQ.1) THEN + wann%band_max(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))) + ELSE + wann%band_max(2) = wann%band_max(1) + END IF + wann%l_byindex = .TRUE. + IF(input%l_wann) THEN + IF (dimension%neigd.LT.MAX(wann%band_max(1),wann%band_max(2))) THEN + dimension%neigd = MAX(wann%band_max(1),wann%band_max(2)) + END IF + END IF + END IF + + xPathA = '/fleurInput/output/wannier/jobList' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + xPathA = '/fleurInput/output/wannier/jobList/text()' + + ! Note: At the moment only 255 characters for the text in this node. Maybe this is not enough. + valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))) + numTokens = countStringTokens(valueString) + ALLOCATE(wann%jobList(numTokens)) + DO i = 1, numTokens + wann%jobList(i) = popFirstStringToken(valueString) + END DO + END IF + + ! Read in optional magnetic circular dichroism parameters + xPathA = '/fleurInput/output/magneticCircularDichroism' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF ((banddos%l_mcd).AND.(numberNodes.EQ.0)) THEN + CALL juDFT_error("mcd is true but magneticCircularDichroism parameters are not set!", calledby = "r_inpXML") + END IF + + IF (numberNodes.EQ.1) THEN + banddos%e_mcd_lo = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@energyLo')) + banddos%e_mcd_up = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@energyUp')) + END IF + + ! Read in optional parameter for unfolding bandstructure of supercell + xPathA = '/fleurInput/output/unfoldingBand' + numberNodes = xmlGetNumberOfNodes(xPathA) + + IF (numberNodes.EQ.1) THEN + banddos%unfoldband = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@unfoldband')) + banddos%s_cell_x = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellX')) + banddos%s_cell_y = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellY')) + banddos%s_cell_z = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellZ')) + END IF + + END IF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! End of output section !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Generate / fill wann%atomlist(:) array - IF (wann%l_atomlist) THEN - absSum = 0 - DO i = 1, atoms%nat - IF (wannAtomList(i)) absSum = absSum + 1 - END DO - wann%atomlist_num = absSum - ALLOCATE(wann%atomlist(wann%atomlist_num)) - j = 1 - DO i = 1, atoms%nat - IF (wannAtomList(i)) THEN - wann%atomlist(j) = i - j = j + 1 - END IF - END DO - ELSE - wann%atomlist_num = atoms%nat - ALLOCATE(wann%atomlist(wann%atomlist_num)) - DO i = 1, atoms%nat - wann%atomlist(i) = i - END DO - END IF - - DEALLOCATE(wannAtomList) + ! Generate / fill wann%atomlist(:) array + IF (wann%l_atomlist) THEN + absSum = 0 + DO i = 1, atoms%nat + IF (wannAtomList(i)) absSum = absSum + 1 + END DO + wann%atomlist_num = absSum + ALLOCATE(wann%atomlist(wann%atomlist_num)) + j = 1 + DO i = 1, atoms%nat + IF (wannAtomList(i)) THEN + wann%atomlist(j) = i + j = j + 1 + END IF + END DO + ELSE + wann%atomlist_num = atoms%nat + ALLOCATE(wann%atomlist(wann%atomlist_num)) + DO i = 1, atoms%nat + wann%atomlist(i) = i + END DO + END IF + + DEALLOCATE(wannAtomList) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Start of non-XML input !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Read in enpara file iff available + ! Read in enpara file iff available - l_enpara = .FALSE. - INQUIRE (file ='enpara',exist= l_enpara) - IF (l_enpara) THEN - CALL enpara%READ(atoms,input%jspins,input%film,.FALSE.) - END IF + l_enpara = .FALSE. + INQUIRE (file ='enpara',exist= l_enpara) + IF (l_enpara) THEN + CALL enpara%READ(atoms,input%jspins,input%film,.FALSE.) + END IF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! End of non-XML input !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - CALL xmlFreeResources() - - !WRITE(*,*) 'Reading of inp.xml file finished' - - DEALLOCATE(speciesNLO) - -END SUBROUTINE r_inpXML - -SUBROUTINE setXCParameters(atoms,namex,relcor,jspins,id_x,id_c,xcpot) - USE m_juDFT - USE m_types - USE m_types_xcpot_inbuild - USE m_types_xcpot_libxc - - IMPLICIT NONE - TYPE(t_atoms),INTENT(IN) :: atoms - CHARACTER(LEN=*), INTENT(IN) :: namex - LOGICAL, INTENT(IN) :: relcor - INTEGER, INTENT(IN) :: jspins,id_c,id_x - CLASS(t_xcpot),INTENT(OUT),ALLOCATABLE :: xcpot - - IF (namex(1:5)=='LibXC') THEN - ALLOCATE(t_xcpot_libxc::xcpot) - ELSE - ALLOCATE(t_xcpot_inbuild::xcpot) - ENDIF - - SELECT TYPE(xcpot) - TYPE IS(t_xcpot_inbuild) - CALL xcpot%init(namex(1:4),relcor,atoms%ntype) - TYPE IS(t_xcpot_libxc) - CALL xcpot%init(jspins,id_x,id_c) - END SELECT - - - END SUBROUTINE setXCParameters - -SUBROUTINE getIntegerSequenceFromString(string, sequence, count) - - IMPLICIT NONE - - CHARACTER(*), INTENT(IN) :: string - INTEGER, ALLOCATABLE, INTENT(OUT) :: sequence(:) - INTEGER, INTENT(OUT) :: count - - INTEGER :: i, length, start, lastNumber, currentNumber, index - LOGICAL singleNumber, comma, dash - - ! 3 cases: 1. a single number - ! 2. number - number - ! 3. comma separated numbers - - length = LEN(string) - count = 0 - start = 1 - singleNumber = .TRUE. - comma = .FALSE. - dash = .FALSE. - lastNumber = 0 - count = 0 - - ! 1. Determine number count - - DO i = 1, length - SELECT CASE (string(i:i)) + CALL xmlFreeResources() + + !WRITE(*,*) 'Reading of inp.xml file finished' + + DEALLOCATE(speciesNLO) + + END SUBROUTINE r_inpXML + + SUBROUTINE setXCParameters(atoms,namex,relcor,jspins,id_x,id_c,xcpot) + USE m_juDFT + USE m_types + USE m_types_xcpot_inbuild + USE m_types_xcpot_libxc + + IMPLICIT NONE + TYPE(t_atoms),INTENT(IN) :: atoms + CHARACTER(LEN=*), INTENT(IN) :: namex + LOGICAL, INTENT(IN) :: relcor + INTEGER, INTENT(IN) :: jspins,id_c,id_x + CLASS(t_xcpot),INTENT(OUT),ALLOCATABLE :: xcpot + + IF (namex(1:5)=='LibXC') THEN + ALLOCATE(t_xcpot_libxc::xcpot) + ELSE + ALLOCATE(t_xcpot_inbuild::xcpot) + ENDIF + + SELECT TYPE(xcpot) + TYPE IS(t_xcpot_inbuild) + CALL xcpot%init(namex(1:4),relcor,atoms%ntype) + TYPE IS(t_xcpot_libxc) + CALL xcpot%init(jspins,id_x,id_c) + END SELECT + + END SUBROUTINE setXCParameters + + SUBROUTINE getIntegerSequenceFromString(string, sequence, count) + + IMPLICIT NONE + + CHARACTER(*), INTENT(IN) :: string + INTEGER, ALLOCATABLE, INTENT(OUT) :: sequence(:) + INTEGER, INTENT(OUT) :: count + + INTEGER :: i, length, start, lastNumber, currentNumber, index + LOGICAL singleNumber, comma, dash + + ! 3 cases: 1. a single number + ! 2. number - number + ! 3. comma separated numbers + + length = LEN(string) + count = 0 + start = 1 + singleNumber = .TRUE. + comma = .FALSE. + dash = .FALSE. + lastNumber = 0 + count = 0 + + ! 1. Determine number count + + DO i = 1, length + SELECT CASE (string(i:i)) CASE ('0':'9') CASE (',') IF ((start.EQ.i).OR.(dash)) THEN @@ -2133,117 +2123,117 @@ SUBROUTINE getIntegerSequenceFromString(string, sequence, count) start = i+1 CASE DEFAULT STOP 'String has wrong syntax (in getIntegerSequenceFromString)' - END SELECT - END DO - IF(start.GT.length) THEN - STOP 'String has wrong syntax (in getIntegerSequenceFromString)' - END IF - READ(string(start:length),*) currentNumber - IF (dash) THEN - count = currentNumber - lastNumber + 1 - ELSE - count = count + 1 - END IF - - IF (ALLOCATED(sequence)) THEN - DEALLOCATE(sequence) - END IF - ALLOCATE(sequence(count)) - - ! 2. Read in numbers iff comma separation ...and store numbers in any case - - IF (singleNumber) THEN - sequence(1) = currentNumber - ELSE IF (dash) THEN - DO i = 1, count - sequence(i) = lastNumber + i - 1 + END SELECT END DO - ELSE - index = 1 - start = 1 - DO i = 1, length - SELECT CASE (string(i:i)) + IF(start.GT.length) THEN + STOP 'String has wrong syntax (in getIntegerSequenceFromString)' + END IF + READ(string(start:length),*) currentNumber + IF (dash) THEN + count = currentNumber - lastNumber + 1 + ELSE + count = count + 1 + END IF + + IF (ALLOCATED(sequence)) THEN + DEALLOCATE(sequence) + END IF + ALLOCATE(sequence(count)) + + ! 2. Read in numbers iff comma separation ...and store numbers in any case + + IF (singleNumber) THEN + sequence(1) = currentNumber + ELSE IF (dash) THEN + DO i = 1, count + sequence(i) = lastNumber + i - 1 + END DO + ELSE + index = 1 + start = 1 + DO i = 1, length + SELECT CASE (string(i:i)) CASE (',') comma = .TRUE. READ(string(start:i-1),*) lastNumber start = i+1 sequence(index) = lastNumber index = index + 1 - END SELECT - END DO - sequence(index) = currentNumber - END IF + END SELECT + END DO + sequence(index) = currentNumber + END IF -END SUBROUTINE getIntegerSequenceFromString + END SUBROUTINE getIntegerSequenceFromString -FUNCTION popFirstStringToken(line) RESULT(firstToken) + FUNCTION popFirstStringToken(line) RESULT(firstToken) - IMPLICIT NONE + IMPLICIT NONE - CHARACTER(*), INTENT(INOUT) :: line - CHARACTER(LEN = LEN(line)) :: firstToken + CHARACTER(*), INTENT(INOUT) :: line + CHARACTER(LEN = LEN(line)) :: firstToken - INTEGER separatorIndex + INTEGER separatorIndex - separatorIndex = 0 - line = TRIM(ADJUSTL(line)) + separatorIndex = 0 + line = TRIM(ADJUSTL(line)) - separatorIndex = INDEX(line,' ') - IF (separatorIndex.LE.1) THEN - firstToken = ' ' - ELSE - firstToken = line(1:separatorIndex-1) - line = line(separatorIndex+1:) - END IF + separatorIndex = INDEX(line,' ') + IF (separatorIndex.LE.1) THEN + firstToken = ' ' + ELSE + firstToken = line(1:separatorIndex-1) + line = line(separatorIndex+1:) + END IF -END FUNCTION popFirstStringToken + END FUNCTION popFirstStringToken -FUNCTION countStringTokens(line) RESULT(tokenCount) + FUNCTION countStringTokens(line) RESULT(tokenCount) - IMPLICIT NONE + IMPLICIT NONE - CHARACTER(*), INTENT(IN) :: line - INTEGER :: tokenCount + CHARACTER(*), INTENT(IN) :: line + INTEGER :: tokenCount - CHARACTER(LEN=LEN(line)) :: tempLine - INTEGER separatorIndex + CHARACTER(LEN=LEN(line)) :: tempLine + INTEGER separatorIndex - tokenCount = 0 + tokenCount = 0 - tempLine = TRIM(ADJUSTL(line)) + tempLine = TRIM(ADJUSTL(line)) - DO WHILE (tempLine.NE.'') - separatorIndex = 0 - separatorIndex = INDEX(tempLine,' ') - IF (separatorIndex.EQ.0) THEN - tempLine = '' - ELSE - tempLine = TRIM(ADJUSTL(tempLine(separatorIndex+1:))) - END IF - tokenCount = tokenCount + 1 - END DO - - END FUNCTION countStringTokens - - FUNCTION priv_read_q_list(path)RESULT(q) - USE m_xmlIntWrapFort - USE m_calculator - IMPLICIT NONE - CHARACTER(len=*),INTENT(in):: path - REAL,ALLOCATABLE::q(:,:) - - INTEGER:: n,i - CHARACTER(len=256):: xpatha,valueString - - n=xmlGetNumberOfNodes(TRIM(ADJUSTL(path))//'/q') - ALLOCATE(q(3,n)) - DO i = 1, n - PRINT *, path,'/q[',i,']' - WRITE(xPathA,"(a,a,i0,a)") TRIM(ADJUSTL(path)),'/q[',i,']' - PRINT *,xpatha - valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) - PRINT *,"Q:",valueString - READ(valueString,*) q(1,i),q(2,i),q(3,i) - END DO - END FUNCTION priv_read_q_list + DO WHILE (tempLine.NE.'') + separatorIndex = 0 + separatorIndex = INDEX(tempLine,' ') + IF (separatorIndex.EQ.0) THEN + tempLine = '' + ELSE + tempLine = TRIM(ADJUSTL(tempLine(separatorIndex+1:))) + END IF + tokenCount = tokenCount + 1 + END DO + + END FUNCTION countStringTokens + + FUNCTION priv_read_q_list(path)RESULT(q) + USE m_xmlIntWrapFort + USE m_calculator + IMPLICIT NONE + CHARACTER(len=*),INTENT(in):: path + REAL,ALLOCATABLE::q(:,:) + + INTEGER:: n,i + CHARACTER(len=256):: xpatha,valueString + + n=xmlGetNumberOfNodes(TRIM(ADJUSTL(path))//'/q') + ALLOCATE(q(3,n)) + DO i = 1, n + PRINT *, path,'/q[',i,']' + WRITE(xPathA,"(a,a,i0,a)") TRIM(ADJUSTL(path)),'/q[',i,']' + PRINT *,xpatha + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))) + PRINT *,"Q:",valueString + READ(valueString,*) q(1,i),q(2,i),q(3,i) + END DO + END FUNCTION priv_read_q_list END MODULE m_rinpXML