From 4f9b29877c5b4fc5ce071173ec3b241cd0e20ccc Mon Sep 17 00:00:00 2001 From: Daniel Wortmann Date: Mon, 4 Feb 2019 21:17:46 +0100 Subject: [PATCH] More coding, just fragments... --- docs/ContributorsGuide.md | 2 + main/mix.F90 | 133 ++----------------------- mix/CMakeLists.txt | 2 + mix/distance.F90 | 66 ++++++++++++ mix/kerker.F90 | 26 ++--- mix/type_mixvector.F90 | 205 ++++++++++++++++++++------------------ 6 files changed, 203 insertions(+), 231 deletions(-) create mode 100644 mix/distance.F90 diff --git a/docs/ContributorsGuide.md b/docs/ContributorsGuide.md index 987cba8e..748cd662 100644 --- a/docs/ContributorsGuide.md +++ b/docs/ContributorsGuide.md @@ -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 diff --git a/main/mix.F90 b/main/mix.F90 index a5cf8c6e..c770256e 100644 --- a/main/mix.F90 +++ b/main/mix.F90 @@ -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 - 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) diff --git a/mix/CMakeLists.txt b/mix/CMakeLists.txt index b886e7ae..2dd6e3cf 100644 --- a/mix/CMakeLists.txt +++ b/mix/CMakeLists.txt @@ -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 diff --git a/mix/distance.F90 b/mix/distance.F90 new file mode 100644 index 00000000..c90e63b8 --- /dev/null +++ b/mix/distance.F90 @@ -0,0 +1,66 @@ +!-------------------------------------------------------------------------------- +! 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 + 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') diff --git a/mix/kerker.F90 b/mix/kerker.F90 index da099bae..e1f16e5a 100644 --- a/mix/kerker.F90 +++ b/mix/kerker.F90 @@ -1,4 +1,4 @@ --------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------- ! 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 diff --git a/mix/type_mixvector.F90 b/mix/type_mixvector.F90 index 99a327df..f4d69a09 100644 --- a/mix/type_mixvector.F90 +++ b/mix/type_mixvector.F90 @@ -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 -- GitLab