Commit 78af809c authored by Daniel Wortmann's avatar Daniel Wortmann

Updated the xml Schema for input and made severl adjustments to force theorem code.

parent a4ea4aa9
......@@ -120,8 +120,8 @@ CONTAINS
CALL mpi_bc_pot(mpi,stars,sphhar,atoms,input,vacuum,&
v%iter,v%mt,v%pw,v%vacz,v%vacxy)
#endif
IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber '/),(/it,v%iter/),&
RESHAPE((/19,13,5,5/),(/2,2/)))
!IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber '/),(/it,v%iter/),&
! RESHAPE((/19,13,5,5/),(/2,2/)))
eig_id=open_eig(&
mpi%mpi_comm,DIMENSION%nbasfcn,DIMENSION%neigd,kpts%nkpt,DIMENSION%jspd,atoms%lmaxd,&
......
......@@ -56,7 +56,6 @@ CONTAINS
TYPE(t_rsoc):: rsoc
TYPE(t_zmat):: zmat
TYPE(t_noco):: noco_tmp
TYPE(t_usdus):: usdus
TYPE(t_lapw) :: lapw
......@@ -84,7 +83,7 @@ CONTAINS
!Calculate radial and angular matrix elements of SOC
!many directions of SOC at once...
CALL spnorb(atoms,noco_tmp,input,mpi, enpara, v%mt, usdus, rsoc,.FALSE.)
CALL spnorb(atoms,noco,input,mpi, enpara, v%mt, usdus, rsoc,.FALSE.)
ALLOCATE(soangl(atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,2,&
atoms%lmaxd,-atoms%lmaxd:atoms%lmaxd,2,SIZE(theta)))
......
......@@ -88,6 +88,8 @@ CONTAINS
!Locals
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,2)/))
......@@ -148,7 +150,9 @@ CONTAINS
TYPE(t_potden),INTENT(IN) :: v
TYPE(t_results),INTENT(IN) :: results
INTEGER,INTENT(IN) :: eig_id
skip=.FALSE.
IF (this%q_done==0) RETURN
this%evsum(0,this%q_done)=results%seigv
CALL ssomat(this%evsum(1:,this%q_done),this%theta,this%phi,eig_id,DIMENSION,atoms,kpts,sym,&
cell,noco, input,mpi, oneD,enpara,v,results)
......
......@@ -154,6 +154,9 @@ CONTAINS
!Locals
INTEGER:: n
CHARACTER(LEN=18):: attributes(6)
IF (this%loopindex==0) RETURN
!Now output the results
call closeXMLElement('Forcetheorem_Loop_JIJ')
CALL openXMLElementPoly('Forcetheorem_JIJ',(/'Configs'/),(/this%no_loops/))
......@@ -192,7 +195,9 @@ CONTAINS
TYPE(t_potden),INTENT(IN) :: v
TYPE(t_results),INTENT(IN) :: results
INTEGER,INTENT(IN) :: eig_id
skip=.FALSE.
IF (this%loopindex==0) RETURN
this%evsum(this%loopindex)=results%seigv
skip=.TRUE.
END FUNCTION jij_eval
......
......@@ -65,6 +65,7 @@ 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/))
END FUNCTION mae_next_job
......@@ -89,7 +90,10 @@ CONTAINS
TYPE(t_potden),INTENT(IN) :: v
TYPE(t_results),INTENT(IN) :: results
INTEGER,INTENT(IN) :: eig_id
IF (this%directions_done==0) THEN
skip=.FALSE.
RETURN
ENDIF
this%evsum(this%directions_done)=results%seigv
skip=.TRUE.
END FUNCTION mae_eval
......@@ -102,6 +106,10 @@ CONTAINS
!Locals
INTEGER:: n
CHARACTER(LEN=12):: attributes(3)
IF (this%directions_done==0) THEN
RETURN
ENDIF
!Now output the results
call closeXMLElement('Forcetheorem_Loop_MAE')
CALL openXMLElementPoly('Forcetheorem_MAE',(/'Angles'/),(/SIZE(this%evsum)/))
......@@ -113,6 +121,7 @@ CONTAINS
reshape((/5,3,6,12,12,12/),(/3,2/)))
END DO
CALL closeXMLElement('Forcetheorem_MAE')
CALL judft_end("Forcetheorem MAE")
END SUBROUTINE mae_postprocess
SUBROUTINE mae_dist(this,mpi)
......
......@@ -74,6 +74,7 @@ CONTAINS
!Locals
INTEGER:: n,q
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)/))
......@@ -84,6 +85,7 @@ CONTAINS
RESHAPE((/1,6,5,12/),(/2,2/)))
ENDDO
CALL closeXMLElement('Forcetheorem_SSDISP')
CALL judft_end("Forcetheorem:SpinSpiralDispersion")
END SUBROUTINE ssdisp_postprocess
SUBROUTINE ssdisp_dist(this,mpi)
......@@ -123,7 +125,9 @@ CONTAINS
TYPE(t_potden),INTENT(IN) :: v
TYPE(t_results),INTENT(IN) :: results
INTEGER,INTENT(IN) :: eig_id
skip=.FALSE.
IF (this%q_done==0) RETURN
this%evsum(this%q_done)=results%seigv
skip=.TRUE.
END FUNCTION ssdisp_eval
......
......@@ -56,4 +56,5 @@ init/brzone2.f90
init/postprocessInput.F90
init/initParallelProcesses.F90
init/old_inp/fleur_init_old.F90
init/lapw_dim.F90
)
!--------------------------------------------------------------------------------
! 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_lapwdim
CONTAINS
SUBROUTINE lapw_dim(kpts,cell,input,noco,oneD,forcetheo,DIMENSION)
!
!*********************************************************************
! determines dimensions of the lapw basis set with |k+G|<rkmax.
! Generalization of the old apws_dim routine
!*********************************************************************
USE m_boxdim
USE m_types
USE m_types_forcetheo_extended
IMPLICIT NONE
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(INOUT) :: noco
TYPE(t_oneD),INTENT(IN) :: oneD
CLASS(t_forcetheo),INTENT(in):: forcetheo
TYPE(t_dimension),INTENT(INOUT)::DIMENSION
INTEGER j1,j2,j3,mk1,mk2,mk3,iofile,ksfft,q,nk,nv,nv2
INTEGER ispin,nvh(2),nv2h(2)
REAL arltv1,arltv2,arltv3,rkm,rk2,r2,s(3),gmaxp,qss(3)
REAL,ALLOCATABLE:: q_vectors(:,:)
REAL :: bkpt(3)
! ..
!
!-------> ABBREVIATIONS
!
! iofile : device number for in and output
! gmax : cut-off wavevector for charge density
! rkmax : cut-off for |g+k|
! gmaxp : gmaxp = gmax/rkmax, ideal: gmaxp=2
! arltv(i) : length of reciprical lattice vector along
! direction (i)
!
!---> Determine rkmax box of size mk1, mk2, mk3,
! for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax
!
CALL boxdim(cell%bmat,arltv1,arltv2,arltv3)
! (add 1+1 due to integer rounding, strange k_vector in BZ)
mk1 = int(input%rkmax/arltv1) + 2
mk2 = int(input%rkmax/arltv2) + 2
mk3 = int(input%rkmax/arltv3) + 2
rkm = input%rkmax
rk2 = rkm*rkm
!Determine the q-vector(s) to use
SELECT TYPE(forcetheo)
TYPE IS (t_forcetheo_ssdisp)
ALLOCATE(q_vectors(3,SIZE(forcetheo%qvec,2)))
q_vectors=forcetheo%qvec
TYPE IS (t_forcetheo_dmi)
ALLOCATE(q_vectors(3,SIZE(forcetheo%qvec,2)))
q_vectors=forcetheo%qvec
TYPE IS (t_forcetheo_jij)
ALLOCATE(q_vectors(3,SIZE(forcetheo%qvec,2)))
q_vectors=forcetheo%qvec
TYPE IS (t_forcetheo) ! DEFAULT
ALLOCATE(q_vectors(3,1))
q_vectors(:,1)=noco%qss
END SELECT
noco%qss=q_vectors(:,1) ! Usually does not do anything, but ensures that in
!force theorem CASE noco%qss is first q-vector in list
DIMENSION%nvd = 0 ; DIMENSION%nv2d = 0
DO q=1,SIZE(q_vectors,2)
qss=q_vectors(:,q)
DO nk=1,kpts%nkpt
bkpt=kpts%bk(:,nk)
!---> obtain vectors
!---> in a spin-spiral calculation different basis sets are used for
!---> the two spin directions, because the cutoff radius is defined
!---> by |G + k +/- qss/2| < rkmax.
DO ispin = 1,2
nv = 0
nv2 = 0
DO j1 = -mk1,mk1
s(1) = bkpt(1) + j1 + (2*ispin - 3)/2.0*noco%qss(1)
DO j2 = -mk2,mk2
s(2) = bkpt(2) + j2 + (2*ispin - 3)/2.0*noco%qss(2)
!---> nv2 for films
s(3) = 0.0
!r2 = dotirp(s,s,cell%bbmat)
r2 = dot_product(matmul(s,cell%bbmat),s)
IF (r2.LE.rk2) nv2 = nv2 + 1
DO j3 = -mk3,mk3
s(3) = bkpt(3) + j3 + (2*ispin - 3)/2.0*noco%qss(3)
!r2 = dotirp(s,s,cell%bbmat)
r2 = dot_product(matmul(s,cell%bbmat),s)
IF (r2.LE.rk2) THEN
nv = nv + 1
END IF
END DO
END DO
END DO
!-odim
IF (oneD%odd%d1) THEN
nv2 = 0
s(1) = 0.0
s(2) = 0.0
DO j3 = -mk3,mk3
s(3) = bkpt(3) + j3 + (2*ispin - 3)/2.0*noco%qss(3)
!r2 = dotirp(s,s,cell%bbmat)
r2 = dot_product(matmul(s,cell%bbmat),s)
IF (r2.LE.rk2) THEN
nv2 = nv2 + 1
END IF
END DO
END IF
!+odim
nvh(ispin) = nv
nv2h(ispin) = nv2
END DO
DIMENSION%nvd=MAX(DIMENSION%nvd,MAX(nvh(1),nvh(2)))
DIMENSION%nv2d=MAX(DIMENSION%nv2d,MAX(nv2h(1),nv2h(2)))
ENDDO !k-loop
ENDDO !q-loop
END SUBROUTINE lapw_dim
SUBROUTINE lapw_fft_dim(cell,input,noco,stars)
!
!*********************************************************************
! determines dimensions of the lapw basis set with |k+G|<rkmax.
! Generalization of the old apws_dim routine
!*********************************************************************
USE m_boxdim
USE m_ifft, ONLY : ifft235
USE m_types
IMPLICIT NONE
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_stars),INTENT(INOUT) :: stars
INTEGER j1,j2,j3,mk1,mk2,mk3,iofile,ksfft,q,nk,nv,nv2
INTEGER ispin,nvh(2),nv2h(2)
REAL arltv1,arltv2,arltv3,rkm,rk2,r2,s(3),gmaxp
REAL,ALLOCATABLE:: q_vectors(:,:)
REAL :: bkpt(3)
! ..
!
!-------> ABBREVIATIONS
!
! iofile : device number for in and output
! gmax : cut-off wavevector for charge density
! rkmax : cut-off for |g+k|
! gmaxp : gmaxp = gmax/rkmax, ideal: gmaxp=2
! arltv(i) : length of reciprical lattice vector along
! direction (i)
!
!---> Determine rkmax box of size mk1, mk2, mk3,
! for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax
!
CALL boxdim(cell%bmat,arltv1,arltv2,arltv3)
! (add 1+1 due to integer rounding, strange k_vector in BZ)
!---> Determine the dimensions kq1d, kq2d, kq3d
! of the dimension of the charge density fft-box
! needed for the fast calculation of pw density
! (add 1 due to integer rounding,
! factor 2 due to positive domain)
!
gmaxp = 2.0
CALL boxdim(cell%bmat,arltv1,arltv2,arltv3)
!
mk1 = int(gmaxp*input%rkmax/arltv1) + 1
mk2 = int(gmaxp*input%rkmax/arltv2) + 1
mk3 = int(gmaxp*input%rkmax/arltv3) + 1
!---> add + 1 in spin spiral calculation, to make sure that all G's are
!---> still within the FFT-box after being shifted by the spin spiral
!---> q-vector.
IF (noco%l_ss) THEN
mk1 = mk1 + 1
mk2 = mk2 + 1
mk3 = mk3 + 1
ENDIF
!
stars%kq1_fft = 2*mk1
stars%kq2_fft = 2*mk2
stars%kq3_fft = 2*mk3
!
!---> fft's are usually fastest for low primes
! (restrict kqid to: kqid= (2**P) * (3**Q) * (5**R)
!
ksfft = 1
! ksfft=(0,1) : KEY OF SELECTING FFT-PRDOGRAM AND RADIX-TYPE
! 0 PROGRAM, RADIX-2 ONLY
! 1 PROGRAM, RADIX-2, RADIX-3,RADIX-5
stars%kq1_fft = ifft235(6,ksfft,stars%kq1_fft,gmaxp)
stars%kq2_fft = ifft235(6,ksfft,stars%kq2_fft,gmaxp)
stars%kq3_fft = ifft235(6,ksfft,stars%kq3_fft,gmaxp)
END SUBROUTINE lapw_fft_dim
END MODULE m_lapwdim
......@@ -9,14 +9,14 @@ MODULE m_postprocessInput
CONTAINS
SUBROUTINE postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
oneD,hybrid,cell,banddos,sliceplot,xcpot,&
noco,dimension,enpara,sphhar,l_opti,noel,l_kpts)
oneD,hybrid,cell,banddos,sliceplot,xcpot,forcetheo,&
noco,DIMENSION,enpara,sphhar,l_opti,noel,l_kpts)
USE m_juDFT
USE m_types
USE m_constants
USE m_julia
USE m_apwsdim
USE m_lapwdim
USE m_ylm
USE m_convndim
USE m_chkmt
......@@ -44,6 +44,7 @@ SUBROUTINE postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
IMPLICIT NONE
TYPE(t_mpi) ,INTENT (IN) :: mpi
CLASS(t_forcetheo),INTENT(IN) :: forcetheo
TYPE(t_input), INTENT(INOUT) :: input
TYPE(t_sym), INTENT(INOUT) :: sym
TYPE(t_stars), INTENT(INOUT) :: stars
......@@ -249,26 +250,14 @@ SUBROUTINE postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
dimension%neigd = minNeigd
END IF
dimension%nvd = 0 ; dimension%nv2d = 0
stars%kq1_fft = 0 ; stars%kq2_fft = 0 ; stars%kq3_fft = 0
!cell%aamat=matmul(transpose(cell%amat),cell%amat)
cell%bbmat=matmul(cell%bmat,transpose(cell%bmat))
DO ikpt = 1,kpts%nkpt
DO i = 1, 3
bk(i) = kpts%bk(i,ikpt)
END DO
!IF (input%film .OR.oneD%odd%d1) THEN
! WRITE(*,*) 'There might be additional work required for the k points here!'
! WRITE(*,*) '...in postprocessInput. See inpeig_dim for comparison!'
!END IF
CALL apws_dim(bk(:),cell,input,noco,oneD,nv,nv2,kq1,kq2,kq3)
stars%kq1_fft = MAX(kq1,stars%kq1_fft)
stars%kq2_fft = MAX(kq2,stars%kq2_fft)
stars%kq3_fft = MAX(kq3,stars%kq3_fft)
DIMENSION%nvd = MAX(DIMENSION%nvd,nv)
DIMENSION%nv2d = MAX(DIMENSION%nv2d,nv2)
END DO ! k-pts
CALL lapw_dim(kpts,cell,input,noco,oneD,forcetheo,DIMENSION)
CALL lapw_fft_dim(cell,input,noco,stars)
obsolete%lepr = 0
......
......@@ -1639,11 +1639,11 @@ SUBROUTINE r_inpXML(&
!!! Start of force-theorem section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
xPathA = '/fleurInput/Forcetheorem'
xPathA = '/fleurInput/forceTheorem'
numberNodes = xmlGetNumberOfNodes(xPathA)
IF (numberNodes.EQ.1) THEN
!Magnetic anisotropy...
xPathA = '/fleurInput/Forcetheorem/MAE'
xPathA = '/fleurInput/forceTheorem/MAE'
numberNodes = xmlGetNumberOfNodes(xPathA)
IF (numberNodes.EQ.1) THEN
lString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta')
......@@ -1655,7 +1655,7 @@ SUBROUTINE r_inpXML(&
END SELECT
ENDIF
!spin-spiral dispersion
xPathA = '/fleurInput/Forcetheorem/SpinSpiralDispersion'
xPathA = '/fleurInput/forceTheorem/spinSpiralDispersion'
numberNodes = xmlGetNumberOfNodes(xPathA)
IF (numberNodes.EQ.1) THEN
ALLOCATE(t_forcetheo_ssdisp::forcetheo)
......@@ -1665,7 +1665,7 @@ SUBROUTINE r_inpXML(&
END SELECT
ENDIF
!dmi
xPathA = '/fleurInput/forcetheorem/DMI'
xPathA = '/fleurInput/forceTheorem/DMI'
numberNodes = xmlGetNumberOfNodes(xPathA)
IF (numberNodes.EQ.1) THEN
lString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta')
......@@ -1673,18 +1673,18 @@ SUBROUTINE r_inpXML(&
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(xPathA//'/QVectors'),lstring,nstring)
CALL forcetheo%init(priv_read_q_list(TRIM(ADJUSTL(xPathA))//'/qVectors'),lstring,nstring)
END SELECT
ENDIF
!jij
xPathA = '/fleurInput/forcetheorem/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(xPathA//'/QVectors'),thetaj,atoms)
CALL forcetheo%init(priv_read_q_list(TRIM(ADJUSTL(xPathA))//'/qVectors'),thetaj,atoms)
END SELECT
ENDIF
......@@ -2157,8 +2157,11 @@ FUNCTION countStringTokens(line) RESULT(tokenCount)
n=evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(path))//'/@count'))
ALLOCATE(q(3,n))
DO i = 1, n
WRITE(xPathA,*) path//'/q[',i,']'
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
......
!--------------------------------------------------------------------------------
! 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
......@@ -139,6 +140,7 @@ CONTAINS
it = 0
ithf = 0
l_cont = ( it < input%itmax )
results%last_distance = -1.0
IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('scfLoop')
......@@ -170,7 +172,8 @@ CONTAINS
scfloop:DO WHILE (l_cont)
it = it + 1
IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber '/)&
,(/it,inden%iter/), RESHAPE((/19,13,5,5/),(/2,2/)))
inDenRot = inDen
!!$ !+t3e
......@@ -358,8 +361,10 @@ CONTAINS
IF (forcetheo%eval(eig_id,DIMENSION,atoms,kpts,sym,&
cell,noco, input,mpi, oneD,enpara,vToT,results)) CYCLE forcetheoloop
cell,noco, input,mpi, oneD,enpara,vToT,results)) THEN
IF ( noco%l_soc .AND. (.NOT. noco%l_noco) ) DIMENSION%neigd=DIMENSION%neigd/2
CYCLE forcetheoloop
ENDIF
CALL force_0(results)! ----> initialise force_old
......
......@@ -186,7 +186,7 @@
END IF
CALL postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
oneD,hybrid,cell,banddos,sliceplot,xcpot,&
oneD,hybrid,cell,banddos,sliceplot,xcpot,forcetheo,&
noco,dimension,enpara,sphhar,l_opti,noel,l_kpts)
IF (mpi%irank.EQ.0) THEN
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment