Commit 4f9b2987 authored by Daniel Wortmann's avatar Daniel Wortmann

More coding, just fragments...

parent 561144c4
......@@ -4,6 +4,8 @@ Contributors guide
Everyone is very welcome to contribute to the enhancement of FLEUR.
Please use the [gitlab service] (https://iffgit.fz-juelich.de/fleur/fleur) to obtain the
latest development version of FLEUR.
- The
##Coding rules for FLEUR:
In no particular order we try to collect items to consider when writing code for FLEUR
......
......@@ -67,7 +67,6 @@ contains
integer :: mmap, mmaph, nmaph, nmap, mapmt, mapvac, mapvac2
integer :: iofl, n_u_keep
logical :: l_exist, l_ldaU, l_densityMatrixPresent, l_pot
real :: dist(6)
real, allocatable :: sm(:), fsm(:), fmMet(:), smMet(:)
character(len=20) :: attributes(2)
complex :: n_mmpTemp(-3:3,-3:3,max(1,atoms%n_u),input%jspins)
......@@ -75,64 +74,9 @@ contains
integer :: ierr(2)
!External functions
real :: CPP_BLAS_sdot
external :: CPP_BLAS_sdot
! YM: I have exported 'vol' from outside, be aware
! IF (film) THEN
! vol = 2.0 * z1 * area
! ELSE
! vol = omtil
! ENDIF
MPI0_a: if( mpi%irank == 0 ) then
l_densityMatrixPresent = any( inDen%mmpMat(:,:,:,:) /= 0.0 )
!In systems without inversions symmetry the interstitial star-
!coefficients are complex. Thus twice as many numbers have to be
!stored.
intfac = 2.0
if ( sym%invs ) intfac = 1.0
!The corresponding is true for the coeff. of the warping vacuum
!density depending on the two dimensional inversion.
vacfac = 2.0
if ( sym%invs2 ) vacfac = 1.0
mmaph = intfac * stars%ng3 + atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd + &
vacfac * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + vacuum%nmzd * vacuum%nvac
mmap = mmaph * input%jspins
!in a non-collinear calculations extra space is needed for the
!off-diag. part of the density matrix. these coeff. are generally
!complex independ of invs and invs2.
if ( noco%l_noco ) then
mmap = mmap + 2 * stars%ng3 + 2 * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + &
2 * vacuum%nmzd * vacuum%nvac
IF (noco%l_mtnocopot) mmap= mmap+ 2*atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd
end if
! LDA+U (start)
n_mmpTemp = inDen%mmpMat
n_u_keep = atoms%n_u
if ( atoms%n_u > 0 ) call u_mix( input, atoms, inDen%mmpMat, outDen%mmpMat )
if ( l_densityMatrixPresent ) then
!In an LDA+U caclulation, also the density matrix is included in the
!supervectors (sm,fsm) if no linear mixing is performed on it.
if ( input%ldauLinMix ) then
atoms%n_u = 0
else
mmap = mmap + 7 * 7 * 2 * atoms%n_u * input%jspins ! add 7*7 complex numbers per atoms%n_u and spin
end if
else
atoms%n_u = 0
end if
! LDA+U (end)
allocate( sm(mmap), fsm(mmap) )
allocate( smMet(mmap), fmMet(mmap) )
dist(:) = 0.0
!determine type of mixing:
!imix=0:straight, imix=o broyden first, imix=5:broyden second
......@@ -154,72 +98,19 @@ contains
write( 6, '(''WARNING : for QUASI-NEWTON METHODS SPINF=1'')' )
end if
call sm%init(oneD,input,vacuum,noco,sym,stars,cell,sphhar,atoms)
call fsm%alloc()
call fmMet%alloc()
!put input charge density into array sm
!(in the spin polarized case the arrays sm and fsm consist of spin up and spin down densities)
call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
intfac, vacfac, inDen, nmap, nmaph, mapmt, mapvac, mapvac2, sm )
!put output charge density into array fsm
call brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
intfac, vacfac, outDen, nmap, nmaph, mapmt, mapvac, mapvac2, fsm )
call sm%from_density(inDen)
call fsm%from_density(outDen)
!store the difference fsm - sm in fsm
fsm(:nmap) = fsm(:nmap) - sm(:nmap)
l_pot = .false.
fsm = fsm - sm
! Apply metric w to fsm and store in fmMet: w |fsm>
call metric( cell, atoms, vacuum, sphhar, input, noco, stars, sym, oneD, &
mmap, nmaph, mapmt, mapvac2, fsm, fmMet, l_pot )
!calculate the distance of charge densities for each spin
IF(hybrid%l_calhf) THEN
CALL openXMLElement('densityConvergence',(/'units ','comment'/),(/'me/bohr^3','HF '/))
ELSE
CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
END IF
DO js = 1,input%jspins
dist(js) = CPP_BLAS_sdot(nmaph,fsm(nmaph*(js-1)+1),1, fmMet(nmaph*(js-1)+1),1)
attributes = ''
WRITE(attributes(1),'(i0)') js
WRITE(attributes(2),'(f20.10)') 1000*SQRT(ABS(dist(js)/cell%vol))
CALL writeXMLElementForm('chargeDensity',(/'spin ','distance'/),attributes,reshape((/4,8,1,20/),(/2,2/)))
IF( hybrid%l_calhf ) THEN
WRITE ( 6,FMT=7901) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
ELSE
WRITE ( 6,FMT=7900) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
END IF
END DO
IF (noco%l_noco) dist(6) = CPP_BLAS_sdot((nmap-2*nmaph), fsm(nmaph*2+1),1,fmMet(nmaph*2+1),1)
IF (noco%l_noco) WRITE (6,FMT=7900) 3,inDen%iter,1000*SQRT(ABS(dist(6)/cell%vol))
fmMet=fsm%apply_metric()
!calculate the distance of total charge and spin density
!|rho/m(o) - rho/m(i)| = |rh1(o) -rh1(i)|+ |rh2(o) -rh2(i)| +/_
! +/_2<rh2(o) -rh2(i)|rh1(o) -rh1(i)>
IF (input%jspins.EQ.2) THEN
dist(3) = CPP_BLAS_sdot(nmaph,fsm,1,fmMet(nmaph+1),1)
dist(4) = dist(1) + dist(2) + 2.0e0*dist(3)
dist(5) = dist(1) + dist(2) - 2.0e0*dist(3)
CALL writeXMLElementFormPoly('overallChargeDensity',(/'distance'/),&
(/1000*SQRT(ABS(dist(4)/cell%vol))/),reshape((/10,20/),(/1,2/)))
CALL writeXMLElementFormPoly('spinDensity',(/'distance'/),&
(/1000*SQRT(ABS(dist(5)/cell%vol))/),reshape((/19,20/),(/1,2/)))
IF( hybrid%l_calhf ) THEN
WRITE ( 6,FMT=8001) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
WRITE ( 6,FMT=8011) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
ELSE
WRITE ( 6,FMT=8000) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
WRITE ( 6,FMT=8010) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
END IF
!dist/vol should always be >= 0 ,
!but for dist=0 numerically you might obtain dist/vol < 0
!(e.g. when calculating non-magnetic systems with jspins=2).
END IF
results%last_distance=maxval(1000*SQRT(ABS(dist/cell%vol)))
DEALLOCATE (smMet,fmMet)
CALL closeXMLElement('densityConvergence')
call distance()
end if MPI0_a
......@@ -227,7 +118,7 @@ contains
IF( input%preconditioning_param /= 0 ) THEN
call kerker(field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, inDen, outDen, fsm ,mapmt,mapvac,mapvac2,nmap,nmaph )
oneD, inDen, outDen, fsm )
ENDIF
MPI0_c: if( mpi%irank == 0 ) then
......@@ -250,10 +141,8 @@ contains
!initiatlize mixed density and extract it with brysh2 call
inDen%mmpMat = CMPLX(0.0,0.0)
CALL brysh2(input,stars,atoms,sphhar,noco,vacuum,sym,sm,oneD,inDen)
DEALLOCATE (sm,fsm)
call sm%to_density(inDen)
!fix charge of the new density
CALL qfix(mpi,stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false., fix)
......
......@@ -2,6 +2,8 @@ set(fleur_F77 ${fleur_F77}
mix/metr_z0.f
)
set(fleur_F90 ${fleur_F90}
mix/type_mixvector.F90
mix/kerker.F90
mix/broyden.F90
mix/broyden2.F90
mix/brysh1.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_distance
contains
subroutine distance(fsm)
real :: dist(6)
dist(:) = 0.0
! Apply metric w to fsm and store in fmMet: w |fsm>
!calculate the distance of charge densities for each spin
IF(hybrid%l_calhf) THEN
CALL openXMLElement('densityConvergence',(/'units ','comment'/),(/'me/bohr^3','HF '/))
ELSE
CALL openXMLElement('densityConvergence',(/'units'/),(/'me/bohr^3'/))
END IF
DO js = 1,input%jspins
dist(js) = dot_product(fsm%vec(nmaph*(js-1)+1:nmaph*(js-1)+nmaph),fmMet%vec(nmaph*(js-1)+1:nmaph*(js-1)+nmaph))
attributes = ''
WRITE(attributes(1),'(i0)') js
WRITE(attributes(2),'(f20.10)') 1000*SQRT(ABS(dist(js)/cell%vol))
CALL writeXMLElementForm('chargeDensity',(/'spin ','distance'/),attributes,reshape((/4,8,1,20/),(/2,2/)))
IF( hybrid%l_calhf ) THEN
WRITE ( 6,FMT=7901) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
ELSE
WRITE ( 6,FMT=7900) js,inDen%iter,1000*SQRT(ABS(dist(js)/cell%vol))
END IF
END DO
IF (noco%l_noco) dist(6) = dot_product(fsm%vec(nmaph*2+1:nmaph*2+nmap-2*nmaph),fmMet%vec(nmaph*2+1:nmaph*2+nmap-2*nmaph))
IF (noco%l_noco) WRITE (6,FMT=7900) 3,inDen%iter,1000*SQRT(ABS(dist(6)/cell%vol))
!calculate the distance of total charge and spin density
!|rho/m(o) - rho/m(i)| = |rh1(o) -rh1(i)|+ |rh2(o) -rh2(i)| +/_
! +/_2<rh2(o) -rh2(i)|rh1(o) -rh1(i)>
IF (input%jspins.EQ.2) THEN
dist(3) = fsm.dot.fmMet
dist(4) = dist(1) + dist(2) + 2.0e0*dist(3)
dist(5) = dist(1) + dist(2) - 2.0e0*dist(3)
CALL writeXMLElementFormPoly('overallChargeDensity',(/'distance'/),&
(/1000*SQRT(ABS(dist(4)/cell%vol))/),reshape((/10,20/),(/1,2/)))
CALL writeXMLElementFormPoly('spinDensity',(/'distance'/),&
(/1000*SQRT(ABS(dist(5)/cell%vol))/),reshape((/19,20/),(/1,2/)))
IF( hybrid%l_calhf ) THEN
WRITE ( 6,FMT=8001) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
WRITE ( 6,FMT=8011) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
ELSE
WRITE ( 6,FMT=8000) inDen%iter,1000*SQRT(ABS(dist(4)/cell%vol))
WRITE ( 6,FMT=8010) inDen%iter,1000*SQRT(ABS(dist(5)/cell%vol))
END IF
!dist/vol should always be >= 0 ,
!but for dist=0 numerically you might obtain dist/vol < 0
!(e.g. when calculating non-magnetic systems with jspins=2).
END IF
results%last_distance=maxval(1000*SQRT(ABS(dist/cell%vol)))
DEALLOCATE (smMet,fmMet)
CALL closeXMLElement('densityConvergence')
--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
! 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.
......@@ -8,12 +8,14 @@ MODULE m_kerker
SUBROUTINE kerker(field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, inDen, outDen, fsm ,mapmt,mapvac,mapvac2,nmap,nmaph )
oneD, inDen, outDen, precon_v )
!Implementation of the Kerker preconditioner by M.Hinzen
USE m_vgen_coulomb
USE m_VYukawaFilm
USE m_juDFT
USE m_types
USE m_types_mixvector
USE m_constants
IMPLICIT NONE
TYPE(t_oneD), INTENT(in) :: oneD
TYPE(t_input), INTENT(in) :: input
......@@ -23,20 +25,18 @@ MODULE m_kerker
TYPE(t_stars), INTENT(in) :: stars
TYPE(t_cell), INTENT(in) :: cell
TYPE(t_sphhar), INTENT(in) :: sphhar
TYPE(t_field), INTENT(in) :: field
TYPE(t_field), INTENT(inout) :: field
TYPE(t_dimension), INTENT(in) :: DIMENSION
TYPE(t_mpi), INTENT(in) :: mpi
TYPE(t_atoms), INTENT(in) :: atoms
TYPE(t_potden), INTENT(in) :: outDen
TYPE(t_potden), INTENT(inout) :: inDen
TYPE(t_potden), INTENT(in) :: outDen,inDen
REAL, ALLOCATABLE,INTENT(OUT) :: fsm(:)
INTEGER, INTENT (OUT) :: mapmt,mapvac,mapvac2,nmap,nmaph
TYPE(t_potden), INTENT(inout) :: outDen
TYPE(t_potden), INTENT(in) :: inDen
TYPE(t_mixvector), INTENT(INOUT) :: precon_v
!Locals
type(t_potden) :: resDen, vYukawa
real :: fix
integer :: lh,n
CALL resDen%init( stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN )
CALL vYukawa%init( stars, atoms, sphhar, vacuum, noco, input%jspins, 4 )
......@@ -73,9 +73,9 @@ MODULE m_kerker
CALL outDen%addPotDen( resDen, inDen )
CALL qfix(mpi,stars, atoms, sym, vacuum, sphhar, input, cell, oneD, outDen, noco%l_noco, .FALSE., .TRUE., fix )
CALL resDen%subPotDen( outDen, inDen )
CALL brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
intfac, vacfac, resDen, nmap, nmaph, mapmt, mapvac, mapvac2, fsm )
call precon_v%from_density(resden)
END IF
! end of preconditioner
END IF MPI0_c
END SUBROUTINE kerker
end MODULE m_kerker
......@@ -2,32 +2,38 @@
! 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_type_mixvector
TYPE t_mixvector
REAL,ALLOCATABLE(:) :: vec(:)
INTEGER :: mmap, mmaph, nmaph, nmap, mapmt, mapvac, mapvac2,intfac
LOGICAL :: l_pot !Is this a potential?
!Here we store the pointers used for metric
TYPE(t_oneD),POINTER :: oneD
TYPE(t_input),POINTER :: input
TYPE(t_vacuum),POINTER :: vacuum
TYPE(t_noco),POINTER :: noco
TYPE(t_sym),POINTER :: sym
TYPE(t_stars),POINTER :: stars
TYPE(t_cell),POINTER :: cell
TYPE(t_sphhar),POINTER :: sphhar
TYPE(t_atoms),POINTER :: atoms
!--------------------------------------------------------------------------------
MODULE m_types_mixvector
use m_types
implicit none
PRIVATE
!Here we store the pointers used for metric
TYPE(t_oneD),POINTER :: oneD
TYPE(t_input),POINTER :: input
TYPE(t_vacuum),POINTER :: vacuum
TYPE(t_noco),POINTER :: noco
TYPE(t_sym),POINTER :: sym
TYPE(t_stars),POINTER :: stars
TYPE(t_cell),POINTER :: cell
TYPE(t_sphhar),POINTER :: sphhar
TYPE(t_atoms),POINTER :: atoms =>null()
INTEGER :: mmap, mmaph, nmaph, nmap, mapmt, mapvac, mapvac2
real :: intfac,vacfac
TYPE,PUBLIC:: t_mixvector
REAL,ALLOCATABLE :: vec(:)
LOGICAL :: l_pot=.false. !Is this a potential?
CONTAINS
PROCEDURE init=>mixvector_init
PROCEDURE from_density=>mixvector_from_density
PROCEDURE to_density=>mixvector_to_density
PROCEDURE apply_metric=>mixvector_metric
PROCEDURE :: init=>mixvector_init
procedure :: alloc=>mixvector_alloc
PROCEDURE :: from_density=>mixvector_from_density
PROCEDURE :: to_density=>mixvector_to_density
PROCEDURE :: apply_metric=>mixvector_metric
END TYPE t_mixvector
INTERFACE assignement(=)
INTERFACE assignment(=)
MODULE PROCEDURE assign_vectors
END INTERFACE assignement
END INTERFACE assignment(=)
INTERFACE OPERATOR (*)
MODULE PROCEDURE multiply_scalar
......@@ -41,8 +47,11 @@ MODULE m_type_mixvector
INTERFACE OPERATOR (.dot.)
MODULE PROCEDURE multiply_dot
END INTERFACE OPERATOR (.dot.)
PRIVATE
public :: Operator(+),operator(-),operator(*),operator(.dot.)
public :: assignment(=)
CONTAINS
SUBROUTINE mixvector_from_density(vec,den)
......@@ -52,8 +61,8 @@ CONTAINS
CLASS(t_mixvector),INTENT(INOUT) :: vec
TYPE(t_potden), INTENT(in) :: Den
CALL brysh1( vec%input, vec%stars, vec%atoms, vec%sphhar, vec%noco, vec%vacuum, vec%sym, vec%oneD, &
vec%intfac, vec%vacfac, Den, vec%nmap, vec%nmaph, vec%mapmt, vec%mapvac, vec%mapvac2, vec%vec)
CALL brysh1( input, stars, atoms, sphhar, noco, vacuum, sym, oneD, &
intfac, vacfac, Den, nmap, nmaph, mapmt, mapvac, mapvac2, vec%vec)
END SUBROUTINE mixvector_from_density
SUBROUTINE mixvector_to_density(vec,den)
......@@ -63,7 +72,7 @@ CONTAINS
CLASS(t_mixvector),INTENT(IN) :: vec
TYPE(t_potden), INTENT(OUT) :: Den
CALL brysh2( vec%input, vec%stars, vec%atoms, vec%sphhar, vec%noco, vec%vacuum, vec%sym, vec%vec,vec%oneD,den)
CALL brysh2( input, stars, atoms, sphhar, noco, vacuum, sym, vec%vec,oneD,den)
END SUBROUTINE mixvector_to_density
......@@ -75,50 +84,50 @@ CONTAINS
TYPE(t_mixvector) :: mvec
mvec=vec
CALL metric( vec%cell, vec%atoms, vec%vacuum, vec%sphhar, vec%input, vec%noco, vec%stars, vec%sym, vec%oneD, &
vec%mmap, vec%nmaph, vec%mapmt, vec%mapvac2, vec%vec, mvec%vec, vec%l_pot )
CALL metric( cell, atoms, vacuum, sphhar, input, noco, stars, sym, oneD, &
mmap, nmaph, mapmt, mapvac2, vec%vec, mvec%vec, vec%l_pot )
END FUNCTION mixvector_metric
SUBROUTINE mixvector_init(vec,oneD,input,vacuum,noco,sym,stars,cell,sphhar,atoms)
SUBROUTINE mixvector_init(vec,oneD_i,input_i,vacuum_i,noco_i,sym_i,stars_i,cell_i,sphhar_i,atoms_i)
USE m_types
IMPLICIT NONE
CLASS(t_mixvector),INTENT(OUT) :: vec
TYPE(t_oneD),INTENT(IN),TARGET :: oneD
TYPE(t_input),INTENT(IN),TARGET :: input
TYPE(t_vacuum),INTENT(IN),TARGET :: vacuum
TYPE(t_noco),INTENT(IN),TARGET :: noco
TYPE(t_sym),INTENT(IN),TARGET :: sym
TYPE(t_stars),INTENT(IN),TARGET :: stars
TYPE(t_cell),INTENT(IN),TARGET :: cell
TYPE(t_sphhar),INTENT(IN),TARGET :: sphhar
TYPE(t_atoms),INTENT(IN),TARGET :: atoms
TYPE(t_oneD),INTENT(IN),TARGET :: oneD_i
TYPE(t_input),INTENT(IN),TARGET :: input_i
TYPE(t_vacuum),INTENT(IN),TARGET :: vacuum_i
TYPE(t_noco),INTENT(IN),TARGET :: noco_i
TYPE(t_sym),INTENT(IN),TARGET :: sym_i
TYPE(t_stars),INTENT(IN),TARGET :: stars_i
TYPE(t_cell),INTENT(IN),TARGET :: cell_i
TYPE(t_sphhar),INTENT(IN),TARGET :: sphhar_i
TYPE(t_atoms),INTENT(IN),TARGET :: atoms_i
if(.not.associated(atoms)) then
!Store pointers to data-types
vec%oneD=>oneD;vec%input=>input;vec%vacuum=>vacuum;vec%noco=>noco
vec%sym=>sym;vec%stars=>stars;vec%cell=>cell;vec%sphhar=>sphhar;vec%atoms=>atoms
oneD=>oneD_i;input=>input_i;vacuum=>vacuum_i;noco=>noco_i
sym=>sym_i;stars=>stars_i;cell=>cell_i;sphhar=>sphhar_i;atoms=>atoms_i
!In systems without inversions symmetry the interstitial star-
!coefficients are complex. Thus twice as many numbers have to be
!stored.
vec%intfac=MERGE(1.0,2.0,sym%invs)
intfac=MERGE(1,2,sym%invs)
!The corresponding is true for the coeff. of the warping vacuum
!density depending on the two dimensional inversion.
vec%vacfac=MERGE(1.0,2.0,sym%invs2)
vec%mmaph = vec%intfac * stars%ng3 + atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd + &
vacfac=MERGE(1,2,sym%invs2)
mmaph = intfac * stars%ng3 + atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd + &
vacfac * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + vacuum%nmzd * vacuum%nvac
vec%mmap =vec%mmaph * input%jspins
mmap =mmaph * input%jspins
!in a non-collinear calculations extra space is needed for the
!off-diag. part of the density matrix. these coeff. are generally
!complex independ of invs and invs2.
IF ( noco%l_noco ) THEN
vec%mmap = vec%mmap + 2 * stars%ng3 + 2 * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + &
mmap = mmap + 2 * stars%ng3 + 2 * vacuum%nmzxyd * ( oneD%odi%n2d - 1 ) * vacuum%nvac + &
2 * vacuum%nmzd * vacuum%nvac
IF (noco%l_mtnocopot) vec%mmap= vec%mmap+ 2*atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd
IF (noco%l_mtnocopot) mmap= mmap+ 2*atoms%ntype * ( sphhar%nlhd + 1 ) * atoms%jmtd
END IF
! LDA+U (start)
PRINT *,"MIXING of LDA+U missing....."
!n_mmpTemp = inDen%mmpMat
......@@ -135,46 +144,50 @@ CONTAINS
! ELSE
! atoms%n_u = 0
! END IF
endif
call vec%alloc()
SUBROUTINE mixvector_alloc(vec)
IMPLICIT NONE
CLASS(t_mixvector),INTENT(OUT) :: vec
ALLOCATE( vec%vec(mmap) )
END SUBROUTINE mixvector_init
!The operators
SUBROUTINE assign_vectors(vec,vecin)
TYPE(t_mixvector),INTENT(OUT)::vec
TYPE(t_mixvector),INTENT(IN) ::vecin
vec=vecin
END SUBROUTINE assign_vectors
FUNCTION multiply_scalar(scalar,vec)RESULT(vecout)
TYPE(t_mixvector),INTENT(IN)::vec
REAL,INTENT(IN) ::scalar
TYPE(t_mixvector) ::vecout
vecout=vecin
vecout%vec=vecout%vec*scalar
END FUNCTION multiply_scalar
FUNCTION add_vectors(vec1,vec2)RESULT(vecout)
TYPE(t_mixvector),INTENT(IN)::vec1,vec2
TYPE(t_mixvector) ::vecout
vecout=vec1
vecout%vec=vec1%vec+vec2%vec
END FUNCTION add_vectors
FUNCTION multiply_dot(vec1,vec2)RESULT(dprod)
TYPE(t_mixvector),INTENT(IN)::vec1,vec2
REAL ::dprod
dprod=dot_PRODUCT(vec1%vec,vec2%vec)
END FUNCTION multiply_dot
FUNCTION add_vectors(vec1,vec2)RESULT(vecout)
TYPE(t_mixvector),INTENT(IN)::vec1,vec2
TYPE(t_mixvector) ::vecout
vecout=vec1
vecout%vec=vec1%vec-vec2%vec
END FUNCTION add_vectors
END SUBROUTINE mixvector_alloc
!The operators
SUBROUTINE assign_vectors(vec,vecin)
TYPE(t_mixvector),INTENT(OUT)::vec
TYPE(t_mixvector),INTENT(IN) ::vecin
vec=vecin
END SUBROUTINE assign_vectors
FUNCTION multiply_scalar(scalar,vec)RESULT(vecout)
TYPE(t_mixvector),INTENT(IN)::vec
REAL,INTENT(IN) ::scalar
TYPE(t_mixvector) ::vecout
vecout=vec
vecout%vec=vecout%vec*scalar
END FUNCTION multiply_scalar
FUNCTION add_vectors(vec1,vec2)RESULT(vecout)
TYPE(t_mixvector),INTENT(IN)::vec1,vec2
TYPE(t_mixvector) ::vecout
vecout=vec1
vecout%vec=vec1%vec+vec2%vec
END FUNCTION add_vectors
FUNCTION multiply_dot(vec1,vec2)RESULT(dprod)
TYPE(t_mixvector),INTENT(IN)::vec1,vec2
REAL ::dprod
dprod=dot_PRODUCT(vec1%vec,vec2%vec)
END FUNCTION multiply_dot
FUNCTION subtract_vectors(vec1,vec2)RESULT(vecout)
TYPE(t_mixvector),INTENT(IN)::vec1,vec2
TYPE(t_mixvector) ::vecout
vecout=vec1
vecout%vec=vec1%vec-vec2%vec
END FUNCTION subtract_vectors
end MODULE m_types_mixvector