Skip to content

GitLab

fleur
fleur
Commit d970cdc8 authored by Gregor Michalicek's avatar Gregor Michalicek
Browse files
Options

More subroutines for the t_force type

parent 74cf9fcc
Hide whitespace changes
Inline Side-by-side
Showing with 307 additions and 236 deletions
cdn/cdnval.F90
View file @ d970cdc8
......@@ -56,8 +56,6 @@ CONTAINS
USE m_vacden
USE m_pwden
USE m_forcea8
USE m_forcea12
USE m_forcea21
USE m_checkdopall
USE m_cdnmt ! calculate the density and orbital moments etc.
USE m_orbmom ! coeffd for orbital moments
......@@ -100,21 +98,21 @@ CONTAINS
TYPE(t_mcd), INTENT(INOUT) :: mcd
TYPE(t_slab), INTENT(INOUT) :: slab
! .. Scalar Arguments ..
! Scalar Arguments
INTEGER, INTENT(IN) :: eig_id,jspin
#ifdef CPP_MPI
INCLUDE 'mpif.h'
#endif
! .. Local Scalars ..
! Local Scalars
INTEGER :: ikpt,jsp_start,jsp_end,ispin,jsp
INTEGER :: iErr,ivac,nbands,noccbd,iType
INTEGER :: skip_t,skip_tt
INTEGER :: nStart,nEnd,nbasfcn
LOGICAL :: l_orbcomprot,l_real, l_write
! ...Local Arrays ..
! Local Arrays
INTEGER, ALLOCATABLE :: jsym(:),ksym(:)
REAL, ALLOCATABLE :: we(:)
REAL, ALLOCATABLE :: eig(:)
......@@ -136,28 +134,28 @@ CONTAINS
l_real = sym%invs.AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco)
IF (noco%l_mperp) THEN
!---> when the off-diag. part of the desinsity matrix, i.e. m_x and
!---> m_y, is calculated inside the muffin-tins (l_mperp = T), cdnval
!---> is called only once. therefore, several spin loops have been
!---> added. if l_mperp = F, these loops run only from jspin - jspin.
! when the off-diag. part of the desinsity matrix, i.e. m_x and
! m_y, is calculated inside the muffin-tins (l_mperp = T), cdnval
! is called only once. therefore, several spin loops have been
! added. if l_mperp = F, these loops run only from jspin - jspin.
jsp_start = 1
jsp_end = 2
ELSE
jsp_start = jspin
jsp_end = jspin
ENDIF
!---> if l_mperp = F, these variables are only needed for one spin
!---> at a time, otherwise for both spins:
ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) ) ! Deallocation before mpi_col_den
ALLOCATE ( g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
ALLOCATE ( flo(atoms%jmtd,2,atoms%nlod,dimension%jspd) )
ALLOCATE ( jsym(dimension%neigd),ksym(dimension%neigd) )
! --> Initializations
! if l_mperp = F, these variables are only needed for one spin
! at a time, otherwise for both spins:
ALLOCATE (f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end)) ! Deallocation before mpi_col_den
ALLOCATE (g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end))
ALLOCATE (flo(atoms%jmtd,2,atoms%nlod,dimension%jspd))
ALLOCATE (jsym(dimension%neigd),ksym(dimension%neigd))
! Initializations
CALL usdus%init(atoms,input%jspins)
CALL denCoeffs%init(atoms,sphhar,jsp_start,jsp_end)
!---> The last entry in denCoeffsOffdiag%init is l_fmpl. It is meant as a switch to a plot of the full magnet.
!---> density without the atomic sphere approximation for the magnet. density. It is not completely implemented (lo's missing).
! The last entry in denCoeffsOffdiag%init is l_fmpl. It is meant as a switch to a plot of the full magnet.
! density without the atomic sphere approximation for the magnet. density. It is not completely implemented (lo's missing).
CALL denCoeffsOffdiag%init(atoms,noco,sphhar,.FALSE.)
CALL force%init1(input,atoms)
CALL orb%init(atoms,noco,jsp_start,jsp_end)
......@@ -168,11 +166,11 @@ CONTAINS
IF ((denCoeffsOffdiag%l_fmpl).AND.(.not.noco%l_mperp)) CALL juDFT_error("for fmpl set noco%l_mperp = T!" ,calledby ="cdnval")
IF ((banddos%ndir.EQ.-3).AND.banddos%dos.AND.oneD%odi%d1) CALL juDFT_error("layer-resolved feature does not work with 1D",calledby ="cdnval")
! calculation of core spectra (EELS) initializations -start-
! calculation of core spectra (EELS) initializations -start-
CALL corespec_init(input,atoms,coreSpecInput)
IF(l_cs.AND.(mpi%isize.NE.1)) CALL juDFT_error('EELS + MPI not implemented', calledby = 'cdnval')
IF(l_cs.AND.jspin.EQ.1) CALL corespec_gaunt()
! calculation of core spectra (EELS) initializations -end-
! calculation of core spectra (EELS) initializations -end-
IF (mpi%irank==0) THEN
WRITE (6,FMT=8000) jspin
......@@ -215,9 +213,9 @@ CONTAINS
nStart = cdnvalKLoop%nStart(ikpt)
nEnd = cdnvalKLoop%nEnd(ikpt)
we=0.0
IF(noccbd.GT.0) we(1:noccbd) = results%w_iks(nStart:nEnd,ikpt,jsp)
IF ((sliceplot%slice).AND.(input%pallst)) we(:) = kpts%wtkpt(ikpt)
we = 0.0
IF (noccbd.GT.0) we(1:noccbd) = results%w_iks(nStart:nEnd,ikpt,jsp)
IF (sliceplot%slice.AND.input%pallst) we(:) = kpts%wtkpt(ikpt)
IF (cdnvalKLoop%l_evp) THEN
IF (nStart > skip_tt) skip_t = 0
......@@ -238,36 +236,32 @@ CONTAINS
IF (noccbd.EQ.0) GO TO 199
! ----> add in spin-doubling factor
we(:noccbd) = 2.0 * we(:noccbd) / input%jspins
we(:noccbd) = 2.0 * we(:noccbd) / input%jspins ! add in spin-doubling factor
!---> valence density in the interstitial and vacuum region
!---> has to be called only once (if jspin=1) in the non-collinear case
! valence density in the interstitial and vacuum region
! has to be called only once (if jspin=1) in the non-collinear case
! ----> valence density in the interstitial region
IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
! valence density in the interstitial region
IF (.NOT.((jspin.EQ.2).AND.noco%l_noco)) THEN
CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,&
jspin,lapw,noccbd,we,eig,den,regCharges%qis,results,force%f_b8,zMat)
END IF
!---> charge of each valence state in this k-point of the SBZ
!---> in the layer interstitial region of the film
IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
! charge of each valence state in this k-point of the SBZ
! in the layer interstitial region of the film
IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
CALL q_int_sl(jspin,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat)
END IF
END IF
!---> valence density in the vacuum region
! valence density in the vacuum region
IF (input%film) THEN
IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
IF (.NOT.((jspin.EQ.2).AND.noco%l_noco)) THEN
CALL vacden(vacuum,dimension,stars,oneD, kpts,input,sym,cell,atoms,noco,banddos,&
gVacMap,we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac0,eig,&
den,regCharges%qvac,regCharges%qvlay,regCharges%qstars,zMat)
END IF
!---> perform Brillouin zone integration and summation over the
!---> bands in order to determine the vacuum energy parameters.
DO ispin = jsp_start,jsp_end
! perform Brillouin zone integration and summation over the
! bands in order to determine the vacuum energy parameters.
DO ispin = jsp_start, jsp_end
DO ivac = 1,vacuum%nvac
regCharges%pvac(ivac,ispin)=regCharges%pvac(ivac,ispin)+dot_product(eig(:noccbd)*regCharges%qvac(:noccbd,ivac,ikpt,ispin),we(:noccbd))
regCharges%svac(ivac,ispin)=regCharges%svac(ivac,ispin)+dot_product(regCharges%qvac(:noccbd,ivac,ikpt,ispin),we(:noccbd))
......@@ -275,9 +269,8 @@ CONTAINS
END DO
END IF
!---> valence density in the atomic spheres
! valence density in the atomic spheres
CALL eigVecCoeffs%init(dimension,atoms,noco,jspin,noccbd)
DO ispin = jsp_start,jsp_end
IF (input%l_f) CALL force%init2(noccbd,input,atoms)
CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,noco,ispin,oneD,&
......@@ -285,21 +278,19 @@ CONTAINS
eigVecCoeffs%ccof(-atoms%llod:,:,:,:,ispin),zMat,eig,force)
IF (atoms%n_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,den%mmpMat(:,:,:,jspin))
!---> perform Brillouin zone integration and summation over the
!---> bands in order to determine the energy parameters for each
!---> atom and angular momentum
IF (.not.sliceplot%slice) THEN
CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,&
skip_t,cdnvalKLoop%l_evp,eigVecCoeffs,usdus,regCharges,mcd,banddos%l_mcd)
! perform Brillouin zone integration and summation over the
! bands in order to determine the energy parameters for each
! atom and angular momentum
CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,&
skip_t,cdnvalKLoop%l_evp,eigVecCoeffs,usdus,regCharges,mcd,banddos%l_mcd)
IF (noco%l_mperp.AND.(ispin == jsp_end)) THEN
CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,regCharges)
END IF
IF (noco%l_mperp.AND.(ispin == jsp_end)) THEN
CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,regCharges)
END IF
!---> layer charge of each valence state in this k-point of the SBZ
!---> from the mt-sphere region of the film
IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
! layer charge of each valence state in this k-point of the SBZ
! from the mt-sphere region of the film
IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
INQUIRE (file='orbcomprot',exist=l_orbcomprot)
......@@ -312,29 +303,22 @@ CONTAINS
IF (noco%l_soc) CALL orbmom(atoms,noccbd,we,ispin,eigVecCoeffs,orb)
IF (input%l_f) THEN
IF (.not.input%l_useapw) THEN
CALL force_a12(atoms,noccbd,sym,dimension,cell,oneD,&
we,ispin,noccbd,usdus,eigVecCoeffs,force,results)
ENDIF
CALL force_a21(input,atoms,dimension,noccbd,sym,oneD,cell,we,ispin,&
enpara%el0(0:,:,ispin),noccbd,eig,usdus,eigVecCoeffs,force,results)
END IF
IF (input%l_f) CALL force%addContribsA21A12(input,atoms,dimension,sym,cell,oneD,enpara,&
usdus,eigVecCoeffs,noccbd,ispin,eig,we,results)
IF(l_cs) CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,dimension%neigd,&
noccbd,results%ef,banddos%sig_dos,eig,we,eigVecCoeffs)
END DO !---> end loop over ispin
END DO ! end loop over ispin
IF (noco%l_mperp) CALL denCoeffsOffdiag%calcCoefficients(atoms,sphhar,sym,eigVecCoeffs,we,noccbd)
199 CONTINUE
IF ((banddos%dos.OR.banddos%vacdos.OR.input%cdinf)) THEN
!--dw since z is no longer an argument of cdninf sympsi has to be called here!
! since z is no longer an argument of cdninf sympsi has to be called here!
IF (banddos%ndir.GT.0) CALL sympsi(lapw,jspin,sym,dimension,nbands,cell,eig,noco,ksym,jsym,zMat)
CALL write_dos(eig_id,ikpt,jspin,regCharges,slab,orbcomp,ksym,jsym,mcd%mcd)
END IF
END DO !---> end of k-point loop
END DO ! end of k-point loop
#ifdef CPP_MPI
DO ispin = jsp_start,jsp_end
......@@ -346,14 +330,14 @@ CONTAINS
IF (mpi%irank==0) THEN
CALL cdnmt(dimension%jspd,atoms,sphhar,noco,jsp_start,jsp_end,&
enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt)
IF(l_cs) CALL corespec_ddscs(jspin,input%jspins)
IF (l_cs) CALL corespec_ddscs(jspin,input%jspins)
DO ispin = jsp_start,jsp_end
IF (input%cdinf) THEN
WRITE (6,FMT=8210) ispin
8210 FORMAT (/,5x,'check continuity of cdn for spin=',i2)
CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,den,ispin)
END IF
IF ((input%l_f)) CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,force,results)
IF (input%l_f) CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,force,results)
END DO
CALL closeXMLElement('mtCharges')
END IF
......
cmake/Files_and_Targets.txt
View file @ d970cdc8
......@@ -55,7 +55,8 @@ eigen/vec_for_lo.f90 eigen/orthoglo.F90 juDFT/usage_data.F90
global/enpara.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
juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 io/w_inpXML.f90 init/julia.f90 global/utility.F90
init/compile_descr.F90 init/kpoints.f90 io/xmlOutput.F90 init/brzone2.f90 cdn/slab_dim.f90 cdn/slabgeom.f90 dos/nstm3.f90 cdn/int_21.f90
cdn/int_21lo.f90 cdn_mt/rhomt21.f90 cdn_mt/rhonmt21.f90)
cdn/int_21lo.f90 cdn_mt/rhomt21.f90 cdn_mt/rhonmt21.f90 force/force_a21.F90 force/force_a21_lo.f90 force/force_a21_U.f90 force/force_a12.f90
eigen/tlmplm_store.F90)
set(fleur_SRC ${fleur_F90} ${fleur_F77})
......
eigen/tlmplm_store.F90
View file @ d970cdc8
......@@ -9,7 +9,7 @@ MODULE m_tlmplm_store
! used to transfer the results from tlmplm&density matrix in case of lda+u from eigen
! into force_a21
! D.W 2014
USE m_types
USE m_types_tlmplm
IMPLICIT NONE
PRIVATE
TYPE(t_tlmplm) :: td_stored
......
force/force_a12.f90
View file @ d970cdc8
......@@ -6,13 +6,15 @@ MODULE m_forcea12
!
CONTAINS
SUBROUTINE force_a12(atoms,nobd,sym, DIMENSION, cell,oneD,&
we,jsp,ne,usdus,eigVecCoeffs,force,results)
USE m_types
we,jsp,ne,usdus,eigVecCoeffs,acoflo,bcoflo,e1cof,e2cof,f_a12,results)
USE m_types_setup
USE m_types_misc
USE m_types_usdus
USE m_types_cdnval
USE m_constants
USE m_juDFT
IMPLICIT NONE
TYPE(t_force),INTENT(INOUT) :: force
TYPE(t_results),INTENT(INOUT) :: results
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
......@@ -27,7 +29,12 @@ CONTAINS
INTEGER, INTENT (IN) :: ne ,jsp
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: we(nobd)
REAL, INTENT(IN) :: we(nobd)
COMPLEX, INTENT(IN) :: acoflo(-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
COMPLEX, INTENT(IN) :: bcoflo(-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
COMPLEX, INTENT(IN) :: e1cof(ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
COMPLEX, INTENT(IN) :: e2cof(ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
COMPLEX, INTENT(INOUT) :: f_a12(3,atoms%ntype)
! ..
! .. Local Scalars ..
COMPLEX a12,cil1,cil2
......@@ -87,8 +94,8 @@ CONTAINS
DO m1 = -l1,l1
lm1 = l1* (l1+1) + m1
DO ie = 1,ne
acof_flapw(ie,lm1) = acof_flapw(ie,lm1) - force%acoflo(m1,ie,ilo,natrun)
bcof_flapw(ie,lm1) = bcof_flapw(ie,lm1) - force%bcoflo(m1,ie,ilo,natrun)
acof_flapw(ie,lm1) = acof_flapw(ie,lm1) - acoflo(m1,ie,ilo,natrun)
bcof_flapw(ie,lm1) = bcof_flapw(ie,lm1) - bcoflo(m1,ie,ilo,natrun)
ENDDO
ENDDO
ENDDO
......@@ -107,7 +114,7 @@ CONTAINS
!
a12 = a12 + CONJG(cil1*&
(acof_flapw(ie,lm1)*usdus%us(l1,n,jsp) + bcof_flapw(ie,lm1)*usdus%uds(l1,n,jsp) ))*cil2*&
(force%e1cof(ie,lm2,natrun)*usdus%us(l2,n,jsp)+ force%e2cof(ie,lm2,natrun)*usdus%uds(l2,n,jsp))*we(ie)
(e1cof(ie,lm2,natrun)*usdus%us(l2,n,jsp)+ e2cof(ie,lm2,natrun)*usdus%uds(l2,n,jsp))*we(ie)
END DO
aaa(1) = alpha(l1,m1)*krondel(l2,l1-1)* krondel(m2,m1+1)
......@@ -224,7 +231,7 @@ CONTAINS
! is also a solution of Schr. equ. if psi is one.
DO i = 1,3
results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a12(i))
force%f_a12(i,n) = force%f_a12(i,n) + forc_a12(i)
f_a12(i,n) = f_a12(i,n) + forc_a12(i)
END DO
!
! write result moved to force_a8
......
force/force_a21.F90
View file @ d970cdc8
MODULE m_forcea21
CONTAINS
SUBROUTINE force_a21(input,atoms,DIMENSION,nobd,sym,oneD,cell,&
we,jsp,epar,ne,eig,usdus,eigVecCoeffs,force,results)
SUBROUTINE force_a21(input,atoms,DIMENSION,sym,oneD,cell,&
we,jsp,epar,ne,eig,usdus,eigVecCoeffs,aveccof,bveccof,cveccof,f_a21,f_b4,results)
! ************************************************************
! Pulay 2nd and 3rd (A17+A20) term force contribution a la Rici
......@@ -24,12 +24,15 @@ CONTAINS
USE m_forcea21lo
USE m_forcea21U
USE m_tlmplm_store
USE m_types
USE m_types_setup
USE m_types_misc
USE m_types_usdus
USE m_types_tlmplm
USE m_types_cdnval
USE m_constants
USE m_juDFT
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_force),INTENT(INOUT) :: force
TYPE(t_results),INTENT(INOUT) :: results
TYPE(t_dimension),INTENT(IN) :: DIMENSION
TYPE(t_oneD),INTENT(IN) :: oneD
......@@ -40,12 +43,16 @@ CONTAINS
TYPE(t_eigVecCoeffs),INTENT(IN) :: eigVecCoeffs
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: nobd
INTEGER, INTENT (IN) :: ne,jsp
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: we(nobd),epar(0:atoms%lmaxd,atoms%ntype)
REAL, INTENT (IN) :: eig(DIMENSION%neigd)
REAL, INTENT(IN) :: we(ne),epar(0:atoms%lmaxd,atoms%ntype)
REAL, INTENT(IN) :: eig(DIMENSION%neigd)
COMPLEX, INTENT(IN) :: aveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
COMPLEX, INTENT(IN) :: bveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
COMPLEX, INTENT(IN) :: cveccof(3,-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
COMPLEX, INTENT(INOUT) :: f_a21(3,atoms%ntype)
COMPLEX, INTENT(INOUT) :: f_b4(3,atoms%ntype)
! ..
! .. Local Scalars ..
COMPLEX dtd,dtu,utd,utu
......@@ -145,10 +152,10 @@ CONTAINS
END IF
DO i = 1,3
a21(i,natrun) = a21(i,natrun) + 2.0*&
AIMAG( CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utu*force%aveccof(i,ie,lm2,natrun)&
+CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utd*force%bveccof(i,ie,lm2,natrun)&
+CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtu*force%aveccof(i,ie,lm2,natrun)&
+CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtd*force%bveccof(i,ie,lm2,natrun))*we(ie)/atoms%neq(n)
AIMAG( CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utu*aveccof(i,ie,lm2,natrun)&
+CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utd*bveccof(i,ie,lm2,natrun)&
+CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtu*aveccof(i,ie,lm2,natrun)&
+CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtd*bveccof(i,ie,lm2,natrun))*we(ie)/atoms%neq(n)
! END i loop
END DO
END IF
......@@ -175,10 +182,10 @@ CONTAINS
DO i = 1,3
DO natrun = natom,natom + atoms%neq(n) - 1
a21(i,natrun) = a21(i,natrun) + 2.0*&
AIMAG(CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utu*force%aveccof(i,ie,lm1,natrun)&
+CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utd*force%bveccof(i,ie,lm1,natrun)&
+CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtu*force%aveccof(i,ie,lm1,natrun)&
+CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtd*force%bveccof(i,ie,lm1,natrun)&
AIMAG(CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utu*aveccof(i,ie,lm1,natrun)&
+CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)) *utd*bveccof(i,ie,lm1,natrun)&
+CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtu*aveccof(i,ie,lm1,natrun)&
+CONJG(eigVecCoeffs%bcof(ie,lm1,natrun,jsp)) *dtd*bveccof(i,ie,lm1,natrun)&
)*we(ie) /atoms%neq(n)
END DO
!
......@@ -194,10 +201,10 @@ CONTAINS
!
!---> add the local orbital and U contribution to a21
!
CALL force_a21_lo(nobd,atoms,jsp,n,we,eig,ne,eigVecCoeffs,force,tlmplm,usdus,a21)
CALL force_a21_lo(atoms,jsp,n,we,eig,ne,eigVecCoeffs,aveccof,bveccof,cveccof,tlmplm,usdus,a21)
IF ((atoms%n_u.GT.0).AND.(i_u.LE.atoms%n_u)) THEN
CALL force_a21_U(nobd,atoms,i_u,n,jsp,we,ne,usdus,v_mmp,eigVecCoeffs,force,a21)
CALL force_a21_U(atoms,i_u,n,jsp,we,ne,usdus,v_mmp,eigVecCoeffs,aveccof,bveccof,cveccof,a21)
END IF
IF (input%l_useapw) THEN
! -> B4 force
......@@ -212,10 +219,10 @@ CONTAINS
we(ie)/atoms%neq(n)*atoms%rmt(n)**2*AIMAG(&
CONJG(eigVecCoeffs%acof(ie,lm1,natrun,jsp)*usdus%us(l1,n,jsp)&
+eigVecCoeffs%bcof(ie,lm1,natrun,jsp)*usdus%uds(l1,n,jsp))*&
(force%aveccof(i,ie,lm1,natrun)*usdus%dus(l1,n,jsp)&
+force%bveccof(i,ie,lm1,natrun)*usdus%duds(l1,n,jsp) )&
-CONJG(force%aveccof(i,ie,lm1,natrun)*usdus%us(l1,n,jsp)&
+force%bveccof(i,ie,lm1,natrun)*usdus%uds(l1,n,jsp) )*&
(aveccof(i,ie,lm1,natrun)*usdus%dus(l1,n,jsp)&
+bveccof(i,ie,lm1,natrun)*usdus%duds(l1,n,jsp) )&
-CONJG(aveccof(i,ie,lm1,natrun)*usdus%us(l1,n,jsp)&
+bveccof(i,ie,lm1,natrun)*usdus%uds(l1,n,jsp) )*&
(eigVecCoeffs%acof(ie,lm1,natrun,jsp)*usdus%dus(l1,n,jsp)&
+eigVecCoeffs%bcof(ie,lm1,natrun,jsp)*usdus%duds(l1,n,jsp)) )
END DO
......@@ -232,15 +239,15 @@ CONTAINS
we(ie)/atoms%neq(n)*atoms%rmt(n)**2*AIMAG(&
CONJG( eigVecCoeffs%acof(ie,lm1,natrun,jsp)* usdus%us(l1,n,jsp)&
+ eigVecCoeffs%bcof(ie,lm1,natrun,jsp)* usdus%uds(l1,n,jsp) ) *&
force%cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp)&
cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp)&
+ CONJG(eigVecCoeffs%ccof(m,ie,lo,natrun,jsp)*usdus%ulos(lo,n,jsp)) *&
( force%aveccof(i,ie,lm1,natrun)* usdus%dus(l1,n,jsp)&
+ force%bveccof(i,ie,lm1,natrun)* usdus%duds(l1,n,jsp)&
+ force%cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp) ) &
- (CONJG( force%aveccof(i,ie,lm1,natrun) *usdus%us(l1,n,jsp)&
+ force%bveccof(i,ie,lm1,natrun) *usdus%uds(l1,n,jsp) ) *&
( aveccof(i,ie,lm1,natrun)* usdus%dus(l1,n,jsp)&
+ bveccof(i,ie,lm1,natrun)* usdus%duds(l1,n,jsp)&
+ cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp) ) &
- (CONJG( aveccof(i,ie,lm1,natrun) *usdus%us(l1,n,jsp)&
+ bveccof(i,ie,lm1,natrun) *usdus%uds(l1,n,jsp) ) *&
eigVecCoeffs%ccof(m,ie,lo,natrun,jsp) *usdus%dulos(lo,n,jsp)&
+ CONJG(force%cveccof(i,m,ie,lo,natrun)*usdus%ulos(lo,n,jsp)) *&
+ CONJG(cveccof(i,m,ie,lo,natrun)*usdus%ulos(lo,n,jsp)) *&
( eigVecCoeffs%acof(ie,lm1,natrun,jsp)*usdus%dus(l1,n,jsp)&
+ eigVecCoeffs%bcof(ie,lm1,natrun,jsp)*usdus%duds(l1,n,jsp)&
+ eigVecCoeffs%ccof(m,ie,lo,natrun,jsp)*usdus%dulos(lo,n,jsp) ) ) )
......@@ -347,8 +354,8 @@ CONTAINS
! IS ALSO A SOLUTION OF SCHR. EQU. IF PSI IS ONE.
DO i = 1,3
results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a21(i) + forc_b4(i))
force%f_a21(i,n) = force%f_a21(i,n) + forc_a21(i)
force%f_b4(i,n) = force%f_b4(i,n) + forc_b4(i)
f_a21(i,n) = f_a21(i,n) + forc_a21(i)
f_b4(i,n) = f_b4(i,n) + forc_b4(i)
END DO
!
! write result moved to force_a8
......
force/force_a21_U.f90
View file @ d970cdc8
MODULE m_forcea21U
CONTAINS
SUBROUTINE force_a21_U(nobd,atoms,i_u,itype,isp,we,ne,usdus,v_mmp,eigVecCoeffs,force,a21)
SUBROUTINE force_a21_U(atoms,i_u,itype,isp,we,ne,usdus,v_mmp,eigVecCoeffs,aveccof,bveccof,cveccof,a21)
!
!***********************************************************************
! This subroutine calculates the lda+U contribution to the HF forces,
......@@ -9,24 +9,27 @@ CONTAINS
!***********************************************************************
!
USE m_constants
USE m_types
USE m_types_setup
USE m_types_usdus
USE m_types_cdnval
IMPLICIT NONE
TYPE(t_usdus),INTENT(IN) :: usdus
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_eigVecCoeffs),INTENT(IN) :: eigVecCoeffs
TYPE(t_force),INTENT(IN) :: force
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: nobd
INTEGER, INTENT (IN) :: itype,isp,ne
INTEGER, INTENT (INOUT) :: i_u ! on input: index for the first U for atom type "itype or higher"
! on exit: index for the first U for atom type "itype+1 or higher"
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: we(nobd)
COMPLEX, INTENT (IN) :: v_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)
REAL, INTENT (INOUT) :: a21(3,atoms%nat)
REAL, INTENT(IN) :: we(ne)
COMPLEX, INTENT(IN) :: v_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)
COMPLEX, INTENT(IN) :: aveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
COMPLEX, INTENT(IN) :: bveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
COMPLEX, INTENT(IN) :: cveccof(3,-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
REAL, INTENT(INOUT) :: a21(3,atoms%nat)
! ..
! .. Local Scalars ..
COMPLEX v_a,v_b,v_c,p1,p2,p3
......@@ -59,8 +62,8 @@ CONTAINS
DO iatom = sum(atoms%neq(:itype-1))+1,sum(atoms%neq(:itype))
DO ie = 1,ne
DO i = 1,3
p1 = (CONJG(eigVecCoeffs%acof(ie,lm,iatom,isp)) * v_a) * force%aveccof(i,ie,lmp,iatom)
p2 = (CONJG(eigVecCoeffs%bcof(ie,lm,iatom,isp)) * v_b) * force%bveccof(i,ie,lmp,iatom)
p1 = (CONJG(eigVecCoeffs%acof(ie,lm,iatom,isp)) * v_a) * aveccof(i,ie,lmp,iatom)
p2 = (CONJG(eigVecCoeffs%bcof(ie,lm,iatom,isp)) * v_b) * bveccof(i,ie,lmp,iatom)
a21(i,iatom) = a21(i,iatom) + 2.0*AIMAG(p1 + p2) * we(ie)/atoms%neq(itype)
END DO
END DO
......@@ -84,11 +87,11 @@ CONTAINS
DO iatom = sum(atoms%neq(:itype-1))+1,sum(atoms%neq(:itype))
DO ie = 1,ne
DO i = 1,3
p1 = v_a * (CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * force%cveccof(i,mp,ie,lo,iatom))
p2 = v_b * (CONJG(eigVecCoeffs%acof(ie,lm,iatom,isp)) * force%cveccof(i,mp,ie,lo,iatom) + &
CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * force%aveccof(i,ie,lmp,iatom))
p3 = v_c * (CONJG(eigVecCoeffs%bcof(ie,lm,iatom,isp)) * force%cveccof(i,mp,ie,lo,iatom) + &
CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * force%bveccof(i,ie,lmp,iatom))
p1 = v_a * (CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * cveccof(i,mp,ie,lo,iatom))
p2 = v_b * (CONJG(eigVecCoeffs%acof(ie,lm,iatom,isp)) * cveccof(i,mp,ie,lo,iatom) + &
CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * aveccof(i,ie,lmp,iatom))
p3 = v_c * (CONJG(eigVecCoeffs%bcof(ie,lm,iatom,isp)) * cveccof(i,mp,ie,lo,iatom) + &
CONJG(eigVecCoeffs%ccof(m,ie,lo,iatom,isp)) * bveccof(i,ie,lmp,iatom))
a21(i,iatom) = a21(i,iatom) + 2.0*AIMAG(p1 + p2 + p3)*we(ie)/atoms%neq(itype)
END DO
END DO
......
force/force_a21_lo.f90
View file @ d970cdc8
......@@ -6,7 +6,8 @@
MODULE m_forcea21lo
CONTAINS
SUBROUTINE force_a21_lo(nobd,atoms,isp,itype,we,eig,ne,eigVecCoeffs,force,tlmplm,usdus,a21)
SUBROUTINE force_a21_lo(atoms,isp,itype,we,eig,ne,eigVecCoeffs,&
aveccof,bveccof,cveccof,tlmplm,usdus,a21)
!
!***********************************************************************
! This subroutine calculates the local orbital contribution to A21,
......@@ -15,21 +16,25 @@ CONTAINS
! p.kurz nov. 1997
!***********************************************************************
!
USE m_types
USE m_types_setup
USE m_types_usdus
USE m_types_tlmplm
USE m_types_cdnval
IMPLICIT NONE
TYPE(t_usdus),INTENT(IN) :: usdus
TYPE(t_tlmplm),INTENT(IN) :: tlmplm
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_eigVecCoeffs),INTENT(IN) :: eigVecCoeffs
TYPE(t_force),INTENT(IN) :: force
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: nobd
INTEGER, INTENT (IN) :: itype,ne,isp
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: we(nobd),eig(:)!(dimension%neigd)
REAL, INTENT (INOUT) :: a21(3,atoms%nat)
REAL, INTENT(IN) :: we(ne),eig(:)!(dimension%neigd)
REAL, INTENT(INOUT) :: a21(3,atoms%nat)
COMPLEX, INTENT(IN) :: aveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
COMPLEX, INTENT(IN) :: bveccof(3,ne,0:atoms%lmaxd