Commit d1e6ef06 authored by Robin Hilgers's avatar Robin Hilgers

Added alignSpinAxisMagn.f90 routine.

parent 271fe1de
......@@ -25,5 +25,6 @@ cdn_mt/setabc1locdn.f90
cdn_mt/setabc1locdn1.f90
cdn_mt/calcDenCoeffs.f90
cdn_mt/magnMomFromDen.f90
cdn_mt/alignSpinAxisMagn.f90
)
!--------------------------------------------------------------------------------
! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and avhttps://gcc.gnu.org/onlinedocs/gfortran/SQRT.htmlailable as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!------------------------------------------------------------------------------
! This routine
!
!
!
! Robin Hilgers, Nov '19
MODULE m_alignSpinAxisMagn
USE m_magnMomFromDen
USE m_types
USE m_flipcdn
CONTAINS
SUBROUTINE rotateMagnetToSpinAxis(noco,vacuum,sphhar,stars&
,sym,oneD,cell,moments,input,atoms,den)
TYPE(t_input), INTENT(INOUT) :: input
TYPE(t_atoms), INTENT(INOUT) :: atoms
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_potden), INTENT(INOUT) :: den
REAL :: moments(3)
CALL magnMomFromDen(input,atoms,noco,den,moments)
CALL flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,atoms%flipSpinPhi,atoms%flipSpinTheta,den)
END SUBROUTINE rotateMagnetToSpinAxis
SUBROUTINE rotateMagnetFromSpinAxis(noco,vacuum,sphhar,stars&
,sym,oneD,cell,moments,input,atoms,den)
TYPE(t_input), INTENT(INOUT) :: input
TYPE(t_atoms), INTENT(INOUT) :: atoms
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_potden), INTENT(INOUT) :: den
CALL flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,-atoms%flipSpinPhi,-atoms%flipSpinTheta,den)
atoms%flipSpinPhi=0
atoms%flipSpinTheta=0
END SUBROUTINE rotateMagnetFromSpinAxis
END MODULE m_alignSpinAxisMagn
......@@ -31,13 +31,6 @@ SUBROUTINE magnMomFromDen(input,atoms,noco,den,moments)
REAL :: eps=1E-10
REAL, ALLOCATABLE :: dummyResults(:,:)
IF(.NOT.ALLOCATED(atoms%phi_mt_avg).AND..NOT.ALLOCATED(atoms%theta_mt_avg)) THEN
ALLOCATE(atoms%phi_mt_avg(atoms%ntype),atoms%theta_mt_avg(atoms%ntype))
ELSE IF (ALLOCATED(atoms%phi_mt_avg).AND.ALLOCATED(atoms%theta_mt_avg))THEN
ELSE
CALL judft_error("Allocation error", calledby="magnMomFromDen.f90")
END IF
ALLOCATE(dummyResults(SIZE(den%mt,3),SIZE(den%mt,4)))
......@@ -98,8 +91,5 @@ DEALLOCATE(dummyResults)
atoms%theta_mt_avg(i)=theta
ENDDO
END SUBROUTINE magnMomFromDen
END MODULE m_magnMomFromDen
......@@ -223,6 +223,8 @@ CONTAINS
ALLOCATE(atoms%flipSpinPhi(atoms%ntype))
ALLOCATE(atoms%flipSpinTheta(atoms%ntype))
ALLOCATE(atoms%flipSpinScale(atoms%ntype))
ALLOCATE(atoms%phi_mt_avg(atoms%ntype))
ALLOCATE(atoms%theta_mt_avg(atoms%ntype))
ALLOCATE(atoms%l_geo(atoms%ntype))
ALLOCATE(atoms%lda_u(4*atoms%ntype))
ALLOCATE(atoms%bmu(atoms%ntype))
......
......@@ -113,7 +113,7 @@ CONTAINS
INTEGER :: iter,iterHF,i
LOGICAL :: l_opti,l_cont,l_qfix,l_real
REAL :: fix
REAL :: moments(3)
REAL, ALLOCATABLE :: moments(:,:)
#ifdef CPP_MPI
INCLUDE 'mpif.h'
......@@ -126,7 +126,7 @@ CONTAINS
CALL fleur_init(mpi,input,field,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,sliceplot,&
banddos,obsolete,enpara,xcpot,results,kpts,hybrid,oneD,coreSpecInput,wann,l_opti)
CALL timestop("Initialization")
IF ( ( input%preconditioning_param /= 0 ) .AND. oneD%odi%d1 ) THEN
CALL juDFT_error('Currently no preconditioner for 1D calculations', calledby = 'fleur')
END IF
......@@ -388,8 +388,12 @@ CONTAINS
dimension,kpts,atoms,sphhar,stars,sym,&
enpara,cell,noco,vTot,results,oneD,coreSpecInput,&
archiveType,xcpot,outDen,EnergyDen)
ALLOCATE(moments(atoms%ntype,3))
CALL magnMomFromDen(input,atoms,noco,Outden,moments)
write(*,*)atoms%phi_mt_avg
write(*,*)atoms%theta_mt_avg
write(*,*) moments
DEALLOCATE(moments)
IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN
! CDN including core charge
! CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, &
......
......@@ -26,7 +26,7 @@ MODULE m_flipcdn
! *******************************************************
CONTAINS
SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell)
SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell,phi,theta,optDen)
USE m_rotdenmat
USE m_constants
USE m_cdn_io
......@@ -36,56 +36,82 @@ SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell)
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_input),INTENT(INOUT) :: input
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_cell),INTENT(IN) :: cell
REAL, OPTIONAL, INTENT(IN) :: phi(atoms%ntype)
REAL, OPTIONAL, INTENT(IN) :: theta(atoms%ntype)
TYPE(t_potden), OPTIONAL,INTENT(INOUT) :: optDen
! Local type instance
TYPE(t_potden) :: den
! Local Scalars
COMPLEX :: rhodummy, imPart12, realPart12
REAL :: rhodumms,fermiEnergyTemp, realPart1, realPart2, imPart1,imPart2, rhodummyR
INTEGER :: i,nt,j,lh,na,mp,ispin,urec,itype,m,i_u
REAL :: rhodumms,fermiEnergyTemp, realPart1, realPart2, imPart1,imPart2, rhodummyR, rotAnglePhi(atoms%ntype),rotAngleTheta(atoms%ntype)
INTEGER :: i,nt,j,lh,na,mp,ispin,urec,itype,m,i_u,k
INTEGER :: archiveType
LOGICAL :: n_exist,l_qfix,l_error, l_flip(atoms%ntype)
LOGICAL :: n_exist,l_qfix,l_error, l_flip(atoms%ntype), scaleSpin(atoms%ntype)
! Local Arrays
CHARACTER(len=80), ALLOCATABLE :: clines(:)
!Flipcdn by optional given angle if lflip is false but routine is called.
DO k=1, atoms%ntype
IF(.NOT.input%lflip.AND.PRESENT(phi).AND.present(theta)) THEN
rotAnglePhi(k)=phi(k)
rotAngleTheta(k)=theta(k)
scaleSpin(k)=.FALSE.
ELSE IF (input%lflip) THEN
!Rotation triggerd by lflip.
rotAnglePhi(k)=atoms%flipSpinPhi(k)
rotAngleTheta(k)=atoms%flipSpinTheta(k)
scaleSpin(k)=atoms%flipSpinScale(k)
ELSE
CALL judft_error("You shouldn't be here. There went something wrong.",calledby="flipcdn")
END IF
END DO
DO itype=1, atoms%ntype
l_flip(itype)=MERGE(.TRUE.,.FALSE.,(atoms%flipSpinPhi(itype).NE.0.0) .OR.(atoms%flipSpinTheta(itype).NE.0.0))
l_flip(itype)=MERGE(.TRUE.,.FALSE.,(rotAnglePhi(itype).NE.0.0) .OR.(rotAngleTheta(itype).NE.0.0))
END DO
!rot_den_mat(alph,beta,rho11,rho22,rho21)
CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
IF(noco%l_noco) THEN
IF(.NOT.PRESENT(optDen)) THEN
CALL den%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
IF(noco%l_noco) THEN
archiveType = CDN_ARCHIVE_TYPE_NOCO_const
ELSE
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
END IF
! read the charge density
CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
! read the charge density
CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,&
CDN_INPUT_DEN_const,0,fermiEnergyTemp,l_qfix,den)
ELSE
den=optDen
END IF
! flip cdn for each atom with rotation angles given
na = 1
DO itype = 1, atoms%ntype
IF (l_flip(itype).AND.(.NOT.atoms%flipSpinScale(itype))) THEN
IF (l_flip(itype).AND.(.NOT.scaleSpin(itype))) THEN
! spherical and non-spherical m.t. charge density
DO lh = 0,sphhar%nlh(atoms%ntypsy(na))
DO j = 1,atoms%jri(itype)
IF (noco%l_mtNocoPot) THEN
rhodummy=CMPLX(den%mt(j,lh,itype,3),den%mt(j,lh,itype,4))
CALL rot_den_mat(atoms%flipSpinPhi(itype),atoms%flipSpinTheta(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy)
CALL rot_den_mat(rotAnglePhi(itype),rotAngleTheta(itype),den%mt(j,lh,itype,1),den%mt(j,lh,itype,2),rhodummy)
den%mt(j,lh,itype,3)=REAL(rhodummy)
den%mt(j,lh,itype,4)=AIMAG(rhodummy)
ELSE
IF (atoms%flipSpinTheta(itype).EQ.(pimach()).AND.atoms%flipSpinPhi(itype).EQ.0) THEN
IF (rotAngleTheta(itype).EQ.(pimach()).AND.rotAnglePhi(itype).EQ.0) THEN
rhodummyR = den%mt(j,lh,itype,1)
den%mt(j,lh,itype,1) = den%mt(j,lh,itype,input%jspins)
den%mt(j,lh,itype,input%jspins) = rhodummyR
......@@ -97,8 +123,8 @@ SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell)
END IF
END DO
END DO
ELSE IF (l_flip(itype).AND.atoms%flipSpinScale(itype)) THEN
IF((atoms%flipSpinTheta(itype).NE.(pimach()) .OR.atoms%flipSpinPhi(itype).NE.0.0)) CALL judft_error("Spinscaling in combination with flipSpin is currently only implemented using flipSpinTheta=Pi and flipSpinPhi=0.0.",calledby="flipcdn")
ELSE IF (l_flip(itype).AND.scaleSpin(itype)) THEN
IF((rotAngleTheta(itype).NE.(pimach()) .OR.rotAnglePhi(itype).NE.0.0)) CALL judft_error("Spinscaling in combination with flipSpin is currently only implemented using flipSpinTheta=Pi and flipSpinPhi=0.0.",calledby="flipcdn")
DO lh = 0,sphhar%nlh(atoms%ntypsy(na))
DO j = 1,atoms%jri(itype)
rhodummy = den%mt(j,lh,itype,1) + den%mt(j,lh,itype,input%jspins)
......@@ -111,9 +137,8 @@ SUBROUTINE flipcdn(atoms,input,vacuum,sphhar,stars,sym,noco,oneD,cell)
na = na + atoms%neq(itype)
END DO
IF (input%l_removeMagnetisationFromInterstitial) THEN
!!This Segment takes care that no interstitial magnetization is written in the the density. Meaning: Off diagonal elements of density matrix set to 0 and diagonal elements of density matrix are equal to their mean value.
IF (input%l_removeMagnetisationFromInterstitial) THEN
!!This Segment takes care that no interstitial magnetization is written in the the density. Meaning: Off diagonal elements of density matrix set to 0 and diagonal elements of density matrix are equal to their mean value.
den%pw(:,2)=(den%pw(:,1)+den%pw(:,2))*0.5 !mean value
den%pw(:,1)=den%pw(:,2)
IF (noco%l_noco) THEN
......@@ -126,7 +151,7 @@ IF (input%l_removeMagnetisationFromInterstitial) THEN
IF (ANY(den%mmpMat(:,:,:,:).NE.0.0).AND.atoms%n_u>0) THEN
DO i_u = 1, atoms%n_u
itype = atoms%lda_u(i_u)%atomType
IF (l_flip(itype).AND.(.NOT.atoms%flipSpinScale(itype))) THEN
IF (l_flip(itype).AND.(.NOT.scaleSpin(itype))) THEN
DO m = -3,3
DO mp = -3,3
IF (noco%l_mtNocoPot) THEN
......@@ -137,15 +162,15 @@ IF (input%l_removeMagnetisationFromInterstitial) THEN
imPart1=AIMAG(den%mmpMat(m,mp,i_u,1))
imPart2=AIMAG(den%mmpMat(m,mp,i_u,2))
imPart12=CMPLX(0,AIMAG(den%mmpMat(m,mp,i_u,3)))
CALL rot_den_mat(atoms%flipSpinPhi(itype),atoms%flipSpinTheta(itype),realPart1,realPart2,&
CALL rot_den_mat(rotAnglePhi(itype),rotAngleTheta(itype),realPart1,realPart2,&
realPart12)
CALL rot_den_mat(atoms%flipSpinPhi(itype),atoms%flipSpinTheta(itype),imPart1,imPart2,&
CALL rot_den_mat(rotAnglePhi(itype),rotAngleTheta(itype),imPart1,imPart2,&
imPart12)
den%mmpMat(m,mp,i_u,1)=CMPLX(realPart1,imPart1)
den%mmpMat(m,mp,i_u,2)=CMPLX(realPart2,imPart2)
den%mmpMat(m,mp,i_u,3)=realPart12+imPart12
ELSE
IF (atoms%flipSpinTheta(itype).EQ.(pimach()).AND.atoms%flipSpinPhi(itype).EQ.0) THEN
IF (rotAngleTheta(itype).EQ.(pimach()).AND.rotAnglePhi(itype).EQ.0) THEN
rhodummyR = den%mmpMat(m,mp,i_u,1)
den%mmpMat(m,mp,i_u,1) = den%mmpMat(m,mp,i_u,input%jspins)
den%mmpMat(m,mp,i_u,input%jspins) = rhodummyR
......@@ -158,10 +183,10 @@ IF (input%l_removeMagnetisationFromInterstitial) THEN
END DO
END DO
ELSE IF (l_flip(itype).AND.(atoms%flipSpinScale(itype))) THEN
ELSE IF (l_flip(itype).AND.(scaleSpin(itype))) THEN
DO m = -3,3
DO mp = -3,3
IF((atoms%flipSpinTheta(itype).NE.pimach() .OR.atoms%flipSpinPhi(itype).NE.0.0)) CALL judft_error("Spinscaling in combination with flipSpin is currently only implemented using flipSpinTheta=Pi and flipSpinPhi=0.0",calledby="flipcdn")
IF((rotAngleTheta(itype).NE.pimach() .OR.rotAnglePhi(itype).NE.0.0)) CALL judft_error("Spinscaling in combination with flipSpin is currently only implemented using flipSpinTheta=Pi and flipSpinPhi=0.0",calledby="flipcdn")
rhodummy = den%mmpMat(m,mp,i_u,1) + den%mmpMat(m,mp,i_u,input%jspins)
rhodumms = den%mmpMat(m,mp,i_u,1) - den%mmpMat(m,mp,i_u,input%jspins)
den%mmpMat(m,mp,i_u,1) = 0.5 * (rhodummy + atoms%bmu(itype) * rhodumms)
......@@ -219,6 +244,7 @@ END DO
CLOSE(40)
END IF
END SUBROUTINE flipcdn
END MODULE m_flipcdn
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