Commit 5d8d29d0 authored by Daniel Wortmann's avatar Daniel Wortmann

Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

parents 4ab04e9c 88cf4538
......@@ -324,6 +324,7 @@
IF (rkappa.GT.zero) rkappa=MAX(rkappa,tol_14)
! gb works also around
zvac = - LOG( tol_14/cabs(VALUE) ) / rkappa
zvac = MIN (2.*vacuum%nmz,abs(zvac)) ! avoid int-overflow in next line
nzvac = INT( zvac/vacuum%delz ) + 1
! IF ( rkappa.GT.zero .AND. real(value).GT.zero ) THEN
IF ( rkappa.GT.zero ) THEN
......
......@@ -52,7 +52,7 @@ inpgen/atom_sym.f inpgen/generator.f inpgen/read_record.f inpgen/soc_or_ssdw.f i
inpgen/bravais_symm.f inpgen/set_atom_core.f inpgen/spg_gen.f global/triang.f
inpgen/lapw_input.f inpgen/struct_input.f inpgen/write_struct.f
io/calculator.f global/ss_sym.f global/soc_sym.f math/inv3.f io/rw_symfile.f
global/sort.f kpoints/kptgen_hybrid.f kpoints/od_kptsgen.f kpoints/bravais.f kpoints/divi.f kpoints/brzone.f
kpoints/kptgen_hybrid.f kpoints/od_kptsgen.f kpoints/bravais.f kpoints/divi.f kpoints/brzone.f
kpoints/kptmop.f kpoints/kpttet.f init/bandstr1.F kpoints/ordstar.f kpoints/fulstar.f kpoints/kprep.f
kpoints/tetcon.f kpoints/kvecon.f init/boxdim.f global/radsra.f math/intgr.F global/differ.f math/inwint.f
math/outint.f xc-pot/gaunt.f math/grule.f
......@@ -60,7 +60,7 @@ kpoints/tetcon.f kpoints/kvecon.f init/boxdim.f global/radsra.f math/intgr.F glo
set(inpgen_F90 ${inpgen_F90} global/constants.f90 io/xsf_io.f90
eigen/orthoglo.F90 juDFT/usage_data.F90 math/ylm4.F90
global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/find_enpara.f90
global/sort.f90 global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/find_enpara.f90
inpgen/closure.f90 inpgen/inpgen_arguments.F90
juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 io/w_inpXML.f90 kpoints/julia.f90 global/utility.F90
init/compile_descr.F90 kpoints/kpoints.f90 io/xmlOutput.F90 kpoints/brzone2.f90 cdn/slab_dim.f90 cdn/slabgeom.f90 dos/nstm3.f90 cdn/int_21.f90
......
......@@ -12,14 +12,17 @@
export FC=$MPIFC
export CC=$MPICC
fi
#ELPA
if [ $ELPA_MODULES ]
then
CLI_ELPA_OPENMP=1
FLEUR_LIBDIR="$FLEUR_LIBDIR $ELPA_LIB"
FLEUR_INCLUDEDIR="$FLEUR_INCLUDEDIR $ELPA_MODULES"
fi
elif module list 2>&1 | grep -q pgi
then
echo "found PGI compiler"
fi
fi
#ELPA
if [ $ELPA_MODULES ]
then
CLI_ELPA_OPENMP=1
FLEUR_LIBDIR="$FLEUR_LIBDIR $ELPA_LIB"
FLEUR_INCLUDEDIR="$FLEUR_INCLUDEDIR $ELPA_MODULES"
fi
......@@ -4,9 +4,11 @@ LINK_LIBRARIES ${FLEUR_LIBRARIES})
if (NOT FLEUR_USE_ELPA_ONENODE)
if (DEFINED CLI_ELPA_OPENMP)
if (FLEUR_USES_GPU)
set(TEST_LIBRARIES "-lelpa_onenode;${FLEUR_LIBRARIES}")
else()
set(TEST_LIBRARIES "-lelpa_onenode_openmp;${FLEUR_LIBRARIES}")
#else()
# set(TEST_LIBRARIES "-lelpa;${FLEUR_LIBRARIES}")
endif()
endif()
try_compile(FLEUR_USE_ELPA_ONENODE ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_ELPA.f90
LINK_LIBRARIES ${TEST_LIBRARIES})
......
......@@ -11,9 +11,11 @@ if (CLI_FLEUR_USE_GPU)
elseif(${CLI_FLEUR_USE_GPU} MATCHES "cuda9.1")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Mcuda=cuda9.1,cc60 -Mcuda=rdc -Mcudalib=cublas")
elseif(${CLI_FLEUR_USE_GPU} MATCHES "nvtx")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Mcuda=cuda9.1,cc60 -Mcuda=rdc -Mcudalib=cublas -lnvToolsExt ")
#set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Mcuda=cuda9.1,cc60,ptxinfo,lineinfo -Mcuda=rdc -Mcudalib=cublas -lnvToolsExt ")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Mcuda=cuda9.1,cc60 -Mcuda=rdc -Mcudalib=cublas -lnvToolsExt ")
elseif(${CLI_FLEUR_USE_GPU} MATCHES "emu")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Mcuda=emu -Mcudalib=cublas -Minfo=accel ")
#set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Mcuda=emu -Mcudalib=cublas -Minfo=accel ")
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Mcuda=emu -Mcudalib=cublas ")
endif()
set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_GPU" "CPP_MANAGED=,MANAGED")
set(FLEUR_DEFINITIONS ${FLEUR_DEFINITIONS} "CPP_GPU" "CPP_MANAGED=,MANAGED")
......
......@@ -45,6 +45,7 @@ CONTAINS
INTEGER :: kernel
CLASS(elpa_t),pointer :: elpa_obj
print*, "ELPA 20180525 started"
err = elpa_init(20180525)
elpa_obj => elpa_allocate()
......@@ -60,6 +61,9 @@ CONTAINS
CALL elpa_obj%set("local_ncols", hmat%matsize2, err)
CALL elpa_obj%set("nblk",hmat%matsize1, err)
CALL elpa_obj%set("blacs_context", -1, err)
#ifdef CPP_GPU
CALL elpa_obj%set("gpu",1,err)
#endif
err = elpa_obj%setup()
CALL hmat%add_transpose(hmat)
......
......@@ -61,7 +61,7 @@ CONTAINS
END DO
gvacl(n2) = SQRT(REAL(gvac(1)**2+gvac(2)**2))
ENDDO k_loop
CALL sort(n2,gvacl,gindex)
CALL sort(gindex(:n2),gvacl)
DO j = 1,n2
! gvac1d, gvac2d are now ordered by increasing length
gvac1d(j)=gvac1(gindex(j))
......
......@@ -70,15 +70,15 @@ CONTAINS
!$OMP MASTER
IF ((atoms%invsat(na).EQ.0) .OR. (atoms%invsat(na).EQ.1)) THEN
IF ((atoms%invsat(na) == 0) .OR. (atoms%invsat(na) == 1)) THEN
!---> if this atom is the first of two atoms related by inversion,
!---> the contributions to the overlap matrix of both atoms are added
!---> at once. where it is made use of the fact, that the sum of
!---> these contributions is twice the real part of the contribution
!---> of each atom. note, that in this case there are twice as many
!---> (2*(2*l+1)) k-vectors (compare abccoflo and comments there).
IF (atoms%invsat(na).EQ.0) invsfct = 1
IF (atoms%invsat(na).EQ.1) invsfct = 2
IF (atoms%invsat(na) == 0) invsfct = 1
IF (atoms%invsat(na) == 1) invsfct = 2
!
DO lo = 1,atoms%nlo(ntyp)
......@@ -127,7 +127,7 @@ CONTAINS
!+t3e
DO nkvec = 1,invsfct* (2*l+1)
locol= lapw%nv(jintsp)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix
IF (MOD(locol-1,mpi%n_size).EQ.mpi%n_rank) THEN
IF (MOD(locol-1,mpi%n_size) == mpi%n_rank) THEN
locol=(locol-1)/mpi%n_size+1 !this is the column in local storage
!-t3e
IF (hmat%l_real) THEN
......@@ -174,7 +174,7 @@ CONTAINS
!---> local orbitals at the same atom and with itself
DO nkvec = 1,invsfct* (2*l+1)
locol = lapw%nv(jintsp)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix
IF (MOD(locol-1,mpi%n_size).EQ.mpi%n_rank) THEN
IF (MOD(locol-1,mpi%n_size) == mpi%n_rank) THEN
locol=(locol-1)/mpi%n_size+1 !this is the column in local storage
!-t3e
!---> calculate the hamiltonian matrix elements with other
......
This diff is collapsed.
......@@ -121,14 +121,14 @@ CONTAINS
ENDIF
!---> update hamiltonian and overlap matrices
IF (jspin1==jspin2) THEN
DO i = mpi%n_rank+1,lapw%nv(jspin1),mpi%n_size
DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
i0=(i-1)/mpi%n_size+1 !local column index
ik = map2(i,jspin1)
ik = map2(i,jspin2)
DO j = 1,i - 1 !TODO check noco case
!---> overlap: only (g-g') parallel=0 '
IF (map2(j,jspin1).EQ.ik) THEN
sij = CONJG(a(i,jspin1))*a(j,jspin1) + &
CONJG(b(i,jspin1))*b(j,jspin1)*ddnv(ik,jspin1)
sij = CONJG(a(i,jspin2))*a(j,jspin2) + &
CONJG(b(i,jspin2))*b(j,jspin2)*ddnv(ik,jspin2)
!+APW_LO
IF (input%l_useapw) THEN
apw_lo = CONJG(a(i,jspin1)* uz(ik,jspin1) + b(i,jspin1)* udz(ik,jspin1) ) &
......@@ -151,7 +151,7 @@ CONTAINS
END IF
ENDDO
!Diagonal term of Overlapp matrix, Hamiltonian later
sij = CONJG(a(i,jspin1))*a(i,jspin1) + CONJG(b(i,jspin1))*b(i,jspin1)*ddnv(ik,jspin1)
sij = CONJG(a(i,jspin2))*a(i,jspin2) + CONJG(b(i,jspin2))*b(i,jspin2)*ddnv(ik,jspin2)
IF (hmat(1,1)%l_real) THEN
smat(s1,s2)%data_r(j,i0) = smat(s1,s2)%data_r(j,i0) + REAL(sij)
ELSE
......@@ -163,15 +163,15 @@ CONTAINS
!---> hamiltonian update
DO i = mpi%n_rank+1,lapw%nv(jspin2),mpi%n_size
i0=(i-1)/mpi%n_size+1 !local column index
ik = map2(i,jspin1)
DO j = 1,i
jk = map2(j,jspin2)
ik = map2(i,jspin2)
DO j = 1,MIN(i,lapw%nv(jspin1))
jk = map2(j,jspin1)
IF (jspin2>jspin1) THEN
hij = CONJG(CONJG(a(j,jspin2))* (tuuv(jk,ik)*a(i,jspin1) +tudv(jk,ik)*b(i,jspin1))&
+ CONJG(b(j,jspin2))* (tddv(jk,ik)*b(i,jspin1) +tduv(jk,ik)*a(i,jspin1)))
hij = CONJG(CONJG(a(j,jspin1))* (tuuv(jk,ik)*a(i,jspin2) +tudv(jk,ik)*b(i,jspin2))&
+ CONJG(b(j,jspin1))* (tddv(jk,ik)*b(i,jspin2) +tduv(jk,ik)*a(i,jspin2)))
ELSE
hij = CONJG(a(i,jspin1))* (tuuv(ik,jk)*a(j,jspin2) +tudv(ik,jk)*b(j,jspin2))&
+ CONJG(b(i,jspin1))* (tddv(ik,jk)*b(j,jspin2) +tduv(ik,jk)*a(j,jspin2))
hij = CONJG(a(i,jspin2))* (tuuv(ik,jk)*a(j,jspin1) +tudv(ik,jk)*b(j,jspin1))&
+ CONJG(b(i,jspin2))* (tddv(ik,jk)*b(j,jspin1) +tduv(ik,jk)*a(j,jspin1))
ENDIF
IF (hmat(1,1)%l_real) THEN
hmat(s1,s2)%data_r(j,i0) = hmat(s1,s2)%data_r(j,i0) + REAL(hij)
......
......@@ -58,15 +58,15 @@ CONTAINS
CALL lapw%phase_factors(i,atoms%taual(:,na),noco%qss,cph(:,i))
ENDDO
IF ((atoms%invsat(na).EQ.0) .OR. (atoms%invsat(na).EQ.1)) THEN
IF ((atoms%invsat(na) == 0) .OR. (atoms%invsat(na) == 1)) THEN
!---> if this atom is the first of two atoms related by inversion,
!---> the contributions to the overlap matrix of both atoms are added
!---> at once. where it is made use of the fact, that the sum of
!---> these contributions is twice the real part of the contribution
!---> of each atom. note, that in this case there are twice as many
!---> (2*(2*l+1)) k-vectors (compare abccoflo and comments there).
IF (atoms%invsat(na).EQ.0) invsfct = 1
IF (atoms%invsat(na).EQ.1) invsfct = 2
IF (atoms%invsat(na) == 0) invsfct = 1
IF (atoms%invsat(na) == 1) invsfct = 2
con = fpi_const/SQRT(cell%omtil)* ((atoms%rmt(ntyp))**2)/2.0
......@@ -82,7 +82,7 @@ CONTAINS
DO nkvec = 1,invsfct* (2*l+1) !Each LO can have several functions
!+t3e
locol = lapw%nv(jintsp)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix
IF (MOD(locol-1,mpi%n_size).EQ.mpi%n_rank) THEN
IF (MOD(locol-1,mpi%n_size) == mpi%n_rank) THEN
locol=(locol-1)/mpi%n_size+1 !this is the column in local storage
!-t3e
k = lapw%kvec(nkvec,lo,na)
......@@ -107,7 +107,7 @@ CONTAINS
!---> orbitals at the same atom, if they have the same l
DO lop = 1, (lo-1)
lp = atoms%llo(lop,ntyp)
IF (l.EQ.lp) THEN
IF (l == lp) THEN
fact3 = con**2 * fl2p1 * (&
alo1(lop)*(alo1(lo) + &
clo1(lo)*ud%uulon(lo,ntyp,isp))+&
......
......@@ -53,6 +53,8 @@ CONTAINS
REAL u(vacuum%nmzd,dimension%nv2d,input%jspins),ud(vacuum%nmzd,dimension%nv2d,input%jspins)
REAL v(3),x(vacuum%nmzd), qssbti(2,2)
! ..
tuuv=0.0;tudv=0.0;tddv=0.0;tduv=0.0
udz=0.0;duz=0.0;ddnv=0.0;udz=0.;uz=0.
tail = .true.
np1 = vacuum%nmzxy + 1
!---> wronksian for the schrodinger equation given by an identity
......
......@@ -179,7 +179,7 @@ CONTAINS
END DO
END DO
CALL sort(n,e,index)
CALL sort(index(:n),e)
! Check if no deep eigenvalue is found
IF (e_min-MINVAL(e(1:n))>1.0) THEN
......
......@@ -52,13 +52,14 @@ CONTAINS
this%evsum=0
END SUBROUTINE dmi_init
SUBROUTINE dmi_start(this,potden)
SUBROUTINE dmi_start(this,potden,l_io)
USE m_types_potden
IMPLICIT NONE
CLASS(t_forcetheo_dmi),INTENT(INOUT):: this
TYPE(t_potden) ,INTENT(INOUT) :: potden
LOGICAL,INTENT(IN) :: l_io
this%q_done=0
CALL this%t_forcetheo%start(potden) !call routine of basis type
CALL this%t_forcetheo%start(potden,l_io) !call routine of basis type
END SUBROUTINE dmi_start
LOGICAL FUNCTION dmi_next_job(this,lastiter,noco)
......@@ -80,6 +81,8 @@ CONTAINS
!Now modify the noco-file
noco%qss=this%qvec(:,this%q_done)
IF (.NOT.this%l_io) RETURN
IF (this%q_done.NE.1) CALL closeXMLElement('Forcetheorem_Loop_DMI')
CALL openXMLElementPoly('Forcetheorem_Loop_DMI',(/'Q-vec:'/),(/this%q_done/))
END FUNCTION dmi_next_job
......@@ -93,24 +96,26 @@ CONTAINS
INTEGER:: n,q
CHARACTER(LEN=12):: attributes(4)
IF (this%q_done==0) RETURN
!Now output the results
call closeXMLElement('Forcetheorem_Loop_DMI')
CALL openXMLElementPoly('Forcetheorem_DMI',(/'qPoints','Angles '/),(/SIZE(this%evsum,2),SIZE(this%evsum,1)/))
DO q=1,SIZE(this%evsum,2)
WRITE(attributes(1),'(i5)') q
WRITE(attributes(2),'(f12.7)') this%evsum(0,q)
CALL writeXMLElementForm('Entry',(/'q ','ev-sum'/),attributes(1:2),&
RESHAPE((/1,6,5,12/),(/2,2/)))
DO n=1,SIZE(this%evsum,1)-1
WRITE(attributes(2),'(f12.7)') this%theta(n)
WRITE(attributes(3),'(f12.7)') this%phi(n)
WRITE(attributes(4),'(f12.7)') this%evsum(n,q)
CALL writeXMLElementForm('Entry',(/'q ','theta ','phi ','ev-sum'/),attributes,&
RESHAPE((/1,5,3,6,5,12,12,12/),(/4,2/)))
END DO
ENDDO
CALL closeXMLElement('Forcetheorem_DMI')
IF (this%l_io) THEN
!Now output the results
CALL closeXMLElement('Forcetheorem_Loop_DMI')
CALL openXMLElementPoly('Forcetheorem_DMI',(/'qPoints','Angles '/),(/SIZE(this%evsum,2),SIZE(this%evsum,1)/))
DO q=1,SIZE(this%evsum,2)
WRITE(attributes(1),'(i5)') q
WRITE(attributes(2),'(f12.7)') this%evsum(0,q)
CALL writeXMLElementForm('Entry',(/'q ','ev-sum'/),attributes(1:2),&
RESHAPE((/1,6,5,12/),(/2,2/)))
DO n=1,SIZE(this%evsum,1)-1
WRITE(attributes(2),'(f12.7)') this%theta(n)
WRITE(attributes(3),'(f12.7)') this%phi(n)
WRITE(attributes(4),'(f12.7)') this%evsum(n,q)
CALL writeXMLElementForm('Entry',(/'q ','theta ','phi ','ev-sum'/),attributes,&
RESHAPE((/1,5,3,6,5,12,12,12/),(/4,2/)))
END DO
ENDDO
CALL closeXMLElement('Forcetheorem_DMI')
ENDIF
CALL judft_end("Forcetheorem DMI")
END SUBROUTINE dmi_postprocess
......
......@@ -93,13 +93,14 @@ CONTAINS
SUBROUTINE jij_start(this,potden)
SUBROUTINE jij_start(this,potden,l_io)
USE m_types_potden
IMPLICIT NONE
CLASS(t_forcetheo_jij),INTENT(INOUT):: this
TYPE(t_potden) ,INTENT(INOUT) :: potden
LOGICAL,INTENT(IN) :: l_io
this%loopindex=0
CALL this%t_forcetheo%start(potden) !call routine of basis type
CALL this%t_forcetheo%start(potden,l_io) !call routine of basis type
END SUBROUTINE jij_start
LOGICAL FUNCTION jij_next_job(this,lastiter,noco)
......@@ -145,7 +146,7 @@ CONTAINS
noco%alph(n) = noco%alph(n) + tpi_const*dot_PRODUCT(noco%qss,this%taual_types(:,n))
ENDDO
IF (.NOT.this%l_io) RETURN
IF (this%loopindex.NE.1) CALL closeXMLElement('Forcetheorem_Loop_JIJ')
CALL openXMLElementPoly('Forcetheorem_Loop_JIJ',(/'Loop index:'/),(/this%loopindex/))
END FUNCTION jij_next_job
......@@ -160,6 +161,8 @@ CONTAINS
CHARACTER(LEN=18):: attributes(6)
IF (this%loopindex==0) RETURN
IF (.NOT.this%l_io) RETURN
!Now output the results
call closeXMLElement('Forcetheorem_Loop_JIJ')
......
......@@ -39,13 +39,14 @@ CONTAINS
END SUBROUTINE mae_init
SUBROUTINE mae_start(this,potden)
SUBROUTINE mae_start(this,potden,l_io)
USE m_types_potden
IMPLICIT NONE
CLASS(t_forcetheo_mae),INTENT(INOUT):: this
TYPE(t_potden) ,INTENT(INOUT) :: potden
LOGICAL,INTENT(IN) :: l_io
this%directions_done=0
CALL this%t_forcetheo%start(potden) !call routine of basis type
CALL this%t_forcetheo%start(potden,l_io) !call routine of basis type
END SUBROUTINE mae_start
......@@ -69,8 +70,8 @@ CONTAINS
noco%theta=this%theta(this%directions_done)
noco%phi=this%phi(this%directions_done)
noco%l_soc=.true.
IF (this%directions_done.NE.1) CALL closeXMLElement('Forcetheorem_Loop_MAE')
CALL openXMLElementPoly('Forcetheorem_Loop_MAE',(/'No'/),(/this%directions_done/))
IF (this%directions_done.NE.1.AND.this%l_io) CALL closeXMLElement('Forcetheorem_Loop_MAE')
IF (this%l_io) CALL openXMLElementPoly('Forcetheorem_Loop_MAE',(/'No'/),(/this%directions_done/))
END FUNCTION mae_next_job
FUNCTION mae_eval(this,eig_id,DIMENSION,atoms,kpts,sym,&
......@@ -113,17 +114,19 @@ CONTAINS
RETURN
ENDIF
!Now output the results
call closeXMLElement('Forcetheorem_Loop_MAE')
CALL openXMLElementPoly('Forcetheorem_MAE',(/'Angles'/),(/SIZE(this%evsum)/))
DO n=1,SIZE(this%evsum)
WRITE(attributes(1),'(f12.7)') this%theta(n)
WRITE(attributes(2),'(f12.7)') this%phi(n)
WRITE(attributes(3),'(f12.7)') this%evsum(n)
CALL writeXMLElementForm('Angle',(/'theta ','phi ','ev-sum'/),attributes,&
reshape((/5,3,6,12,12,12/),(/3,2/)))
END DO
CALL closeXMLElement('Forcetheorem_MAE')
IF (this%l_io) THEN
!Now output the results
CALL closeXMLElement('Forcetheorem_Loop_MAE')
CALL openXMLElementPoly('Forcetheorem_MAE',(/'Angles'/),(/SIZE(this%evsum)/))
DO n=1,SIZE(this%evsum)
WRITE(attributes(1),'(f12.7)') this%theta(n)
WRITE(attributes(2),'(f12.7)') this%phi(n)
WRITE(attributes(3),'(f12.7)') this%evsum(n)
CALL writeXMLElementForm('Angle',(/'theta ','phi ','ev-sum'/),attributes,&
RESHAPE((/5,3,6,12,12,12/),(/3,2/)))
END DO
CALL closeXMLElement('Forcetheorem_MAE')
ENDIF
CALL judft_end("Forcetheorem MAE")
END SUBROUTINE mae_postprocess
......
......@@ -38,24 +38,30 @@ CONTAINS
this%evsum=0
END SUBROUTINE ssdisp_init
SUBROUTINE ssdisp_start(this,potden)
SUBROUTINE ssdisp_start(this,potden,l_io)
USE m_types_potden
IMPLICIT NONE
CLASS(t_forcetheo_ssdisp),INTENT(INOUT):: this
TYPE(t_potden) ,INTENT(INOUT) :: potden
LOGICAL,INTENT(IN) :: l_io
this%q_done=0
CALL this%t_forcetheo%start(potden) !call routine of basis type
CALL this%t_forcetheo%start(potden,l_io) !call routine of basis type
IF (SIZE(potden%pw,2)<2) RETURN
!Average out magnetic part of potential/charge in INT+Vacuum
potden%pw(:,1)=(potden%pw(:,1)+potden%pw(:,2))/2.0
potden%pw(:,2)=potden%pw(:,1)
IF (SIZE(potden%pw,2)==3) potden%pw(:,3)=0.0
potden%vacz(:,:,1)=(potden%vacz(:,:,1)+potden%vacz(:,:,2))/2.0
potden%vacxy(:,:,:,1)=(potden%vacxy(:,:,:,1)+potden%vacxy(:,:,:,2))/2.0
potden%vacz(:,:,2)=potden%vacz(:,:,1)
potden%vacxy(:,:,:,2)=potden%vacxy(:,:,:,1)
!Off diagonal part
IF (SIZE(potden%pw,2)==3) THEN
potden%pw(:,3)=0.0
potden%vacz(:,:,3:)=0.0
potden%vacxy(:,:,:,3)=0.0
END IF
END SUBROUTINE ssdisp_start
......@@ -78,6 +84,7 @@ CONTAINS
!Now modify the noco-file
noco%qss=this%qvec(:,this%q_done)
IF (.NOT.this%l_io) RETURN
IF (this%q_done.NE.1) CALL closeXMLElement('Forcetheorem_Loop_SSDISP')
CALL openXMLElementPoly('Forcetheorem_Loop_SSDISP',(/'Q-vec:'/),(/this%q_done/))
END FUNCTION ssdisp_next_job
......@@ -92,15 +99,17 @@ CONTAINS
CHARACTER(LEN=12):: attributes(4)
IF (this%q_done==0) RETURN
!Now output the results
CALL closeXMLElement('Forcetheorem_Loop_SSDISP')
CALL openXMLElementPoly('Forcetheorem_SSDISP',(/'qvectors'/),(/SIZE(this%evsum)/))
DO q=1,SIZE(this%evsum)
WRITE(attributes(1),'(i5)') q
WRITE(attributes(2),'(f12.7)') this%evsum(q)
CALL writeXMLElementForm('Entry',(/'q ','ev-sum'/),attributes(1:2),&
RESHAPE((/1,6,5,12/),(/2,2/)))
ENDDO
CALL closeXMLElement('Forcetheorem_SSDISP')
IF (this%l_io) THEN
CALL closeXMLElement('Forcetheorem_Loop_SSDISP')
CALL openXMLElementPoly('Forcetheorem_SSDISP',(/'qvectors'/),(/SIZE(this%evsum)/))
DO q=1,SIZE(this%evsum)
WRITE(attributes(1),'(i5)') q
WRITE(attributes(2),'(f12.7)') this%evsum(q)
CALL writeXMLElementForm('Entry',(/'q ','ev-sum'/),attributes(1:2),&
RESHAPE((/1,6,5,12/),(/2,2/)))
ENDDO
CALL closeXMLElement('Forcetheorem_SSDISP')
ENDIF
CALL judft_end("Forcetheorem:SpinSpiralDispersion")
END SUBROUTINE ssdisp_postprocess
......
......@@ -7,7 +7,6 @@ global/radsra.f
global/radsrd.F
global/radsrdn.f
global/soc_sym.f
global/sort.f
global/ss_sym.f
global/starf.f
global/triang.f
......@@ -23,6 +22,7 @@ global/constants.f90
global/matrix_copy.F90
global/checkdop.F90
global/checkdopall.f90
global/sort.f90
global/genMTBasis.f90
global/chkmt.f90
global/convn.f90
......
......@@ -274,7 +274,7 @@
CALL juDFT_error("Too many atoms at same position.",calledby ="chkmt")
END IF
numNearestNeighbors(n) = MIN(maxCubeAtoms,iNeighborAtom)
CALL sort(iNeighborAtom,sqrDistances,distIndexList)
CALL sort(distIndexList(:iNeighborAtom),sqrDistances(:iNeighborAtom))
DO i = 1, numNearestNeighbors(n)
nearestNeighbors(i,n) = neighborAtoms(distIndexList(i))
nearestNeighborDists(i,n) = SQRT(sqrDistances(distIndexList(i)))
......@@ -293,7 +293,7 @@
! Sort distances and set MT radii for the atoms
CALL sort(atoms%ntype,nearestAtomDists,sortedDistList)
CALL sort(sortedDistList,nearestAtomDists)
rmt1 = -1.0
minRmts = -1.0
DO i = 1, atoms%ntype
......
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_sort
CONTAINS
SUBROUTINE sort(
> n,sk,
< js)
c********************************************************************
c heapsort routine
c input: n = number of objects
c sk = array of objects to be sorted
c output: js(i) = index of i'th smallest object m.w.
c modified to avoid numerical uncertainties for equal-lentgh
c vectors and to make sorting machine independent aug. 90, r.p.
c********************************************************************
IMPLICIT NONE
c
C .. Scalar Arguments ..
INTEGER, INTENT (IN) :: n
C ..
C .. Array Arguments ..
REAL, INTENT (IN) :: sk(n)
INTEGER, INTENT (OUT) :: js(n)
C ..
C .. Local Scalars ..
REAL eps,q
INTEGER i,ind,ir,j,l
C ..
C .. Data statements ..
DATA eps/1.e-10/
C ..
c
IF (n == 0) RETURN ! Nothing to do
IF (n == 1) THEN ! Not much to do
js(1) = 1
RETURN
END IF
DO i = 1,n
js(i) = i
ENDDO
c
l = n/2 + 1
ir = n
20 CONTINUE
IF (l.GT.1) THEN
l = l - 1
ind = js(l)
q = sk(ind)
ELSE
ind = js(ir)
q = sk(ind)
js(ir) = js(1)
ir = ir - 1
IF (ir.EQ.1) THEN
js(1) = ind
RETURN
END IF
END IF
i = l
j = l + l
30 IF (j.LE.ir) THEN
IF (j.LT.ir) THEN
c if(sk(js(j)).lt.sk(js(j+1))) j=j+1
IF ((sk(js(j+1))-sk(js(j))).GT.eps) j = j + 1
END IF
c if(q.lt.sk(js(j))) then
IF ((sk(js(j))-q).GT.eps) THEN
js(i) = js(j)
i = j
j = j + j
ELSE
j = ir + 1
END IF
GO TO 30
END IF
js(i) = ind
GO TO 20
END SUBROUTINE
END
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_sort
USE m_judft
CONTAINS
SUBROUTINE sort(ind,lv,lv1)
!********************************************************************
! heapsort routine
! input: lv = array of objects to be sorted
! lv1 = second array to use as secondary sort key
! output: ind(i) = index of i'th smallest object
!********************************************************************
IMPLICIT NONE
!
REAL, INTENT (IN) :: lv(:)
INTEGER, INTENT (OUT) :: ind(:)
REAL,INTENT(IN),OPTIONAL :: lv1(:)
! ..
! .. Local Scalars ..
REAL eps,q,q1
INTEGER i,idx,ir,j,l,n
REAL,ALLOCATABLE :: llv(:)
! ..
! .. Data statements ..
DATA eps/1.e-10/
! ..
!
n=SIZE(ind)
IF (n>SIZE(lv)) CALL judft_error("BUG: incosistent dimensions")
ALLOCATE(llv(n))
IF (PRESENT(lv1)) THEN
IF (n>SIZE(lv1)) CALL judft_error("BUG: incosistent dimensions")
llv=lv1
ELSE
llv=(/(1.*i,i=1,n)/)
END IF
IF (n == 0) RETURN ! Nothing to do
IF (n == 1) THEN ! Not much to do
ind(1) = 1
RETURN
END IF
DO i = 1,n
ind(i) = i
ENDDO
!
l = n/2 + 1
ir = n
DO
IF (l.GT.1) THEN
l = l - 1
idx = ind(l)
q = lv(idx)
q1= llv(idx)
ELSE
idx = ind(ir)
q = lv(idx)
q1= llv(idx)
ind(ir) = ind(1)
ir = ir - 1
IF (ir.EQ.1) THEN
ind(1) = idx
RETURN
END IF