diff --git a/cdn/cdntot.f90 b/cdn/cdntot.f90 index 993d146edda8d7477e9b9392b07de57111fae23b..2622c6278e2f48e88e4fc5d1bbd693ce666dc2d7 100644 --- a/cdn/cdntot.f90 +++ b/cdn/cdntot.f90 @@ -27,7 +27,7 @@ CONTAINS INTEGER :: jsp, j, ivac, nz, n REAL :: q2(vacuum%nmz), w, rht1(vacuum%nmzd,2,input%jspins) COMPLEX :: x(stars%ng3) - + qtot = 0.0 qistot = 0.0 DO jsp = 1,input%jspins @@ -103,7 +103,7 @@ CONTAINS call tmp_potden%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN) call init_mt_grid(input%jspins, atoms, sphhar, xcpot, sym) do n_atm =1,atoms%ntype - call mt_from_grid(atoms, sphhar, n_atm, input%jspins, mt(:,:,n_atm), & + call mt_from_grid(atoms, sym, sphhar, n_atm, input%jspins, mt(:,:,n_atm), & tmp_potden%mt(:,0:,n_atm,:)) do i=1,atoms%jri(n_atm) @@ -156,11 +156,11 @@ CONTAINS REAL qmt(atoms%ntype,input%jspins),qvac(2,input%jspins) INTEGER, ALLOCATABLE :: lengths(:,:) CHARACTER(LEN=20) :: attributes(6), names(6) - + CALL timestart("cdntot") call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, den, & q, qis, qmt, qvac, qtot, qistot) - + IF (input%film) THEN ALLOCATE(lengths(4+vacuum%nvac,2)) ELSE @@ -209,7 +209,7 @@ CONTAINS REAL, INTENT(in) :: q(:), qis(:), qmt(:,:), qvac(:,:), qtot, qistot character(len=*), intent(in), optional :: hint integer :: n_mt - + if(present(hint)) write (*,*) "DEN of ", hint write (*,*) "q = ", q diff --git a/fleurinput/constants.f90 b/fleurinput/constants.f90 index c7ee1badb93f216c595b9d6bcf95327005308480..1e9e5ad8333708ff9060c3558e515968091718de 100755 --- a/fleurinput/constants.f90 +++ b/fleurinput/constants.f90 @@ -13,7 +13,7 @@ MODULE m_constants PROCEDURE:: read_xml=>read_xml_constants PROCEDURE:: mpi_bc=>mpi_bc_constants END TYPE t_constants - + REAL :: warp_factor=1.0 !should be set from input later INTEGER, PARAMETER :: noState_const = 0 INTEGER, PARAMETER :: coreState_const = 1 @@ -26,9 +26,10 @@ MODULE m_constants COMPLEX, PARAMETER :: ImagUnit=(0.0,1.0) REAL, PARAMETER :: hartree_to_ev_const=27.21138602 ! value from 2014 CODATA recommended values. Uncertainty is 0.00000017 REAL, PARAMETER :: eVac0Default_const = -0.25 - CHARACTER(len=9), PARAMETER :: version_const = 'fleur 27' + CHARACTER(len=9), PARAMETER :: version_const = 'fleur 30' + CHARACTER(len=49), PARAMETER :: version_const_MaX = ' MaX-Release 3.1 (www.max-centre.eu)' REAL, PARAMETER :: boltzmann_const = 3.1668114e-6 ! value is given in Hartree/Kelvin - + INTEGER, PARAMETER :: POTDEN_TYPE_OTHER = 0 ! POTDEN_TYPE <= 0 ==> undefined INTEGER, PARAMETER :: POTDEN_TYPE_POTTOT = 1 ! 0 < POTDEN_TYPE <= 1000 ==> potential INTEGER, PARAMETER :: POTDEN_TYPE_POTCOUL = 2 @@ -36,7 +37,7 @@ MODULE m_constants INTEGER, PARAMETER :: POTDEN_TYPE_POTYUK = 4 INTEGER, PARAMETER :: POTDEN_TYPE_EnergyDen = 5 INTEGER, PARAMETER :: POTDEN_TYPE_DEN = 1001 ! 1000 < POTDEN_TYPE ==> density - + CHARACTER(2),DIMENSION(0:103),PARAMETER :: namat_const=(/& 'va',' H','He','Li','Be',' B',' C',' N',' O',' F','Ne',& 'Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca','Sc','Ti',& @@ -47,7 +48,7 @@ MODULE m_constants 'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb',& 'Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa',' U','Np','Pu',& 'Am','Cm','Bk','Cf','Es','Fm','Md','No','Lw'/) - + CHARACTER(7),DIMENSION(29),PARAMETER :: coreStateList_const=(/& '(1s1/2)','(2s1/2)','(2p1/2)','(2p3/2)','(3s1/2)',& '(3p1/2)','(3p3/2)','(4s1/2)','(3d3/2)','(3d5/2)',& @@ -55,50 +56,50 @@ MODULE m_constants '(5p1/2)','(5p3/2)','(6s1/2)','(4f5/2)','(4f7/2)',& '(5d3/2)','(5d5/2)','(6p1/2)','(6p3/2)','(7s1/2)',& '(5f5/2)','(5f7/2)','(6d3/2)','(6d5/2)' /) - + INTEGER,DIMENSION(29),PARAMETER :: coreStateNumElecsList_const=(/& ! This is the number of electrons per spin 1, 1, 1, 2, 1, 1, 2, 1, 2, 3, 1, 2, 1, 2,& 3, 1, 2, 1, 3, 4, 2, 3, 1, 2, 1, 3, 4, 2, 3/) - + INTEGER,DIMENSION(29),PARAMETER :: coreStateNprncList_const=(/& 1, 2, 2, 2, 3, 3, 3, 4, 3, 3, 4, 4, 5, 4, 4,& 5, 5, 6, 4, 4, 5, 5, 6, 6, 7, 5, 5, 6, 6/) INTEGER,DIMENSION(29),PARAMETER :: coreStateKappaList_const=(/& -1,-1, 1,-2,-1, 1,-2,-1, 2,-3, 1,-2,-1, 2,-3,& 1,-2,-1, 3,-4, 2,-3, 1,-2,-1, 3,-4, 2,-3/) - + CHARACTER(4),DIMENSION(6),PARAMETER :: nobleGasConfigList_const=(/'[He]','[Ne]','[Ar]','[Kr]','[Xe]','[Rn]'/) - + INTEGER,DIMENSION(6),PARAMETER :: nobleGasNumStatesList_const=(/1, 4, 7, 12, 17, 24/) - + CONTAINS - + REAL PURE FUNCTION pimach() IMPLICIT NONE ! This subprogram supplies the value of the constant PI correct to ! machine precision where - + ! PI=3.1415926535897932384626433832795028841971693993751058209749446 - + pimach = 3.1415926535897932 END FUNCTION pimach - + REAL ELEMENTAL FUNCTION c_light(fac) IMPLICIT NONE ! This subprogram supplies the value of c according to - ! NIST standard 13.1.99 + ! NIST standard 13.1.99 ! Hartree and Rydbergs changed by fac = 1.0 or 2.0 - + REAL, INTENT (IN) :: fac c_light = 137.0359895e0 * fac * warp_factor !c_light = 1e6*fac END FUNCTION c_light - + SUBROUTINE read_xml_constants(this,xml) USE m_types_xml CLASS(t_constants),INTENT(INout)::this TYPE(t_xml),INTENT(in) ::xml - + IF (xml%GetNumberOfNodes('/fleurInput/calculationSetup/expertModes/@warp_factor')==1)& warp_factor=evaluateFirstOnly(xml%GetAttributeValue('/fleurInput/calculationSetup/expertModes/@warp_factor')) END SUBROUTINE read_xml_constants diff --git a/fleurinput/types_wannier.f90 b/fleurinput/types_wannier.f90 index d37a405087a6378fc36b8041f63baa050ae90d21..4e0223314a3bdd6463135ef940ddb6786972c76b 100644 --- a/fleurinput/types_wannier.f90 +++ b/fleurinput/types_wannier.f90 @@ -13,6 +13,29 @@ MODULE m_types_wannier ! type for wannier-functions ! TYPE,EXTENDS(t_fleurinput_base):: t_wann + !New parameters not handled correctly yet... + LOGICAL :: l_socmatvec + LOGICAL :: l_socmatvecrs + LOGICAL :: l_mmn0_unf_to_spn_unf + LOGICAL :: l_mmn0_to_spn_unf + LOGICAL :: l_mmn0_to_spn + LOGICAL :: l_mmn0_to_spn2 + LOGICAL :: l_mmn0_unf_to_spn + LOGICAL :: l_perpmag_unf_to_tor_unf + LOGICAL :: l_perpmag_to_tor_unf + LOGICAL :: l_perpmag_to_tor + LOGICAL :: l_perpmag_unf_to_tor + LOGICAL :: l_hsomtxvec_unf_to_lmpzsoc_unf + LOGICAL :: l_hsomtxvec_to_lmpzsoc_unf + LOGICAL :: l_hsomtxvec_to_lmpzsoc + LOGICAL :: l_hsomtxvec_unf_to_lmpzsoc + LOGICAL :: l_hsomtx_unf_to_hsoc_unf + LOGICAL :: l_hsomtx_to_hsoc_unf + LOGICAL :: l_hsomtx_to_hsoc + LOGICAL :: l_hsomtx_unf_to_hsoc + INTEGER :: perpmagl + LOGICAL :: l_perpmagatlres + INTEGER :: wan90version =3 INTEGER :: oc_num_orbs =0 INTEGER, ALLOCATABLE :: oc_orbs(:) @@ -133,7 +156,7 @@ MODULE m_types_wannier PROCEDURE :: read_xml => read_xml_wannier PROCEDURE :: mpi_bc => mpi_bc_wannier END TYPE t_wann - + PUBLIC t_wann CONTAINS diff --git a/global/checkdop.F90 b/global/checkdop.F90 index d7f87b6749934a8f881d43752b37d68ecc94cc0b..2e5e4b7aca22ba773820eb06281b1dc18cbe54ce 100644 --- a/global/checkdop.F90 +++ b/global/checkdop.F90 @@ -24,7 +24,7 @@ ! .. ! .. Scalar Arguments .. TYPE(t_dimension),INTENT(IN) :: dimension - type(t_sphhar),intent(in) :: sphhar + type(t_sphhar),intent(in) :: sphhar TYPE(t_stars),INTENT(IN) :: stars TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_sym),INTENT(IN) :: sym @@ -52,7 +52,7 @@ l_cdn = .FALSE. ! By default we assume that the input is a potential. IF (potden%potdenType.LE.0) CALL juDFT_error('unknown potden type', calledby='checkdop') IF (potden%potdenType.GT.1000) l_cdn = .TRUE. ! potdenTypes > 1000 are reserved for densities - + ! .. ! .. @@ -92,7 +92,7 @@ DO j = 1,np IF (.NOT.oneD%odi%d1) THEN CALL starf2(& - & sym%nop2,stars%ng2,stars%kv2,sym%mrot,sym%symor,sym%tau,p(1,j),sym%invtab,& + & sym%nop2,stars%ng2,stars%kv2,sym%mrot,sym%symor,sym%tau,p(1:3,j),sym%invtab,& & sf2)!keep v2(j) = potden%vacz(1,ivac,jsp) DO k = 2,stars%ng2 @@ -190,7 +190,7 @@ ENDDO help = help + potden%mt(atoms%jri(n),lh,n,jsp) * s ENDDO - v2(j) = help * ir2 + v2(j) = help * ir2 IF (j.LE.8) THEN !CALL cotra1(p(1,j),rcc,cell%bmat) rcc=MATMUL(cell%bmat,p(:,j))/tpi_const diff --git a/global/find_enpara.f90 b/global/find_enpara.f90 index 2612066540328762d5531d1edb8cf4c693c69c17..11cd9f85f12a06531f6ea0f3b72d548b1f7de2f7 100644 --- a/global/find_enpara.f90 +++ b/global/find_enpara.f90 @@ -11,31 +11,29 @@ MODULE m_find_enpara PUBLIC:: find_enpara CONTAINS - + !> Function to determine the energy parameter given the quantum number and the potential !! Different schemes are implemented. Nqn (main quantum number) is used as a switch. !! This code was previously in lodpot.f - REAL FUNCTION find_enpara(lo,l,n,jsp,nqn,atoms,irank,vr)RESULT(e) - USE m_types_atoms + REAL FUNCTION find_enpara(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)RESULT(e) + USE m_types_setup USE m_radsra USE m_differ - Use m_xmlOutput USE m_constants IMPLICIT NONE LOGICAL,INTENT(IN):: lo INTEGER,INTENT(IN):: l,n,nqn,jsp REAL,INTENT(OUT) :: e_lo,e_up TYPE(t_atoms),INTENT(IN)::atoms - INTEGER,INTENT(IN) ::irank REAL,INTENT(IN):: vr(:) - IF (nqn>0) e=priv_method1(lo,l,n,jsp,nqn,atoms,irank,vr) - IF (nqn<0) e=priv_method2(lo,l,n,jsp,nqn,atoms,irank,vr) + IF (nqn>0) e=priv_method1(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up) + IF (nqn<0) e=priv_method2(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up) END FUNCTION find_enpara - REAL FUNCTION priv_method1(lo,l,n,jsp,nqn,atoms,irank,vr)RESULT(e) - USE m_types_atoms + REAL FUNCTION priv_method1(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)RESULT(e) + USE m_types_setup USE m_radsra USE m_differ USE m_constants @@ -44,7 +42,6 @@ CONTAINS INTEGER,INTENT(IN):: l,n,nqn,jsp REAL,INTENT(OUT) :: e_lo,e_up TYPE(t_atoms),INTENT(IN)::atoms - INTEGER,INTENT(IN) ::irank REAL,INTENT(IN):: vr(:) @@ -54,12 +51,12 @@ CONTAINS REAL d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c LOGICAL start ! .. - ! .. Local Arrays .. + ! .. Local Arrays .. REAL, ALLOCATABLE :: f(:,:),vrd(:) CHARACTER(LEN=20) :: attributes(6) c=c_light(1.0) - !Core potential setup done for each n,l now + !Core potential setup done for each n,l now d = EXP(atoms%dx(n)) ! set up core-mesh rn = atoms%rmt(n) @@ -84,10 +81,10 @@ CONTAINS node = nqn - (l+1) IF (node<0) CALL judft_error("Error in setup of energy-parameters",hint="This could e.g. happen if you try to use 1p-states") - e = 0.0 + e = 0.0 ! determine upper edge nodeu = -1 ; start = .TRUE. - DO WHILE ( nodeu <= node ) + DO WHILE ( nodeu <= node ) CALL radsra(e,l,vr(:),atoms%rmsh(1,n),& atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2)) IF ( ( nodeu > node ) .AND. start ) THEN @@ -105,7 +102,7 @@ CONTAINS IF (node /= 0) THEN ! determine lower edge nodeu = node + 1 - DO WHILE ( nodeu >= node ) + DO WHILE ( nodeu >= node ) CALL radsra(e,l,vr(:),atoms%rmsh(1,n),& atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2)) e = e - 0.01 @@ -113,7 +110,7 @@ CONTAINS ENDDO e_lo = e ELSE - e_lo = -9.99 + e_lo = -9.99 ENDIF ! calculate core e = (e_up+e_lo)/2 @@ -125,34 +122,15 @@ CONTAINS fn = REAL(nqn) ; fl = REAL(l) ; fj = fl-0.5 CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),& rn,d,msh,vrd, e1, f(:,1),f(:,2),ierr) - e = (2.0*e + e1 ) / 3.0 + e = (2.0*e + e1 ) / 3.0 ENDIF - - IF (irank == 0) THEN - attributes = '' - WRITE(attributes(1),'(i0)') n - WRITE(attributes(2),'(i0)') jsp - WRITE(attributes(3),'(i0,a1)') nqn, ch(l) - WRITE(attributes(4),'(f8.2)') e_lo - WRITE(attributes(5),'(f8.2)') e_up - WRITE(attributes(6),'(f16.10)') e - IF (lo) THEN - CALL writeXMLElementForm('loAtomicEP',(/'atomType ','spin ','branch ',& - 'branchLowest ','branchHighest','value '/),& - attributes,RESHAPE((/10,4,6,12,13,5,6,1,3,8,8,16/),(/6,2/))) - ELSE - CALL writeXMLElementForm('atomicEP',(/'atomType ','spin ','branch ',& - 'branchLowest ','branchHighest','value '/),& - attributes,RESHAPE((/12,4,6,12,13,5,6,1,3,8,8,16/),(/6,2/))) - ENDIF - WRITE(6,'(a6,i5,i2,a1,a12,f6.2,a3,f6.2,a13,f8.4)') ' Atom',n,nqn,ch(l),' branch from',& - e_lo, ' to',e_up,' htr. ; e_l =',e - ENDIF + + END FUNCTION priv_method1 - REAL FUNCTION priv_method2(lo,l,n,jsp,nqn,atoms,irank,vr)RESULT(e) - USE m_types_atoms + REAL FUNCTION priv_method2(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)RESULT(e) + USE m_types_setup USE m_radsra USE m_differ USE m_constants @@ -161,7 +139,6 @@ CONTAINS INTEGER,INTENT(IN):: l,n,nqn,jsp REAL,INTENT(OUT) :: e_lo,e_up TYPE(t_atoms),INTENT(IN)::atoms - INTEGER,INTENT(IN) ::irank REAL,INTENT(IN):: vr(:) INTEGER j,ilo,i @@ -170,14 +147,14 @@ CONTAINS REAL d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c LOGICAL start ! .. - ! .. Local Arrays .. + ! .. Local Arrays .. REAL, ALLOCATABLE :: f(:,:),vrd(:) CHARACTER(LEN=20) :: attributes(6) - + c=c_light(1.0) - !Core potential setup done for each n,l now + !Core potential setup done for each n,l now d = EXP(atoms%dx(n)) ! set up core-mesh rn = atoms%rmt(n) @@ -266,7 +243,7 @@ CONTAINS ! determince notches by intersection ldmt= -99.0 !ldmt = logarithmic derivative @ MT boundary lnd = -l-1 - DO WHILE ( ABS(ldmt-lnd) .GE. 1E-07) + DO WHILE ( ABS(ldmt-lnd) .GE. 1E-07) e = (e_up+e_lo)/2 CALL radsra(e,l,vr(:),atoms%rmsh(1,n),& atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2)) @@ -279,25 +256,6 @@ CONTAINS END IF END DO - IF (irank == 0) THEN - attributes = '' - WRITE(attributes(1),'(i0)') n - WRITE(attributes(2),'(i0)') jsp - WRITE(attributes(3),'(i0,a1)') ABS(nqn), ch(l) - WRITE(attributes(4),'(f16.10)') ldmt - WRITE(attributes(5),'(f16.10)') e - IF (lo) THEN - CALL writeXMLElementForm('heloAtomicEP',(/'atomType ','spin ','branch ',& - 'logDerivMT ','value '/),& - attributes(1:5),reshape((/8,4,6,12,5+17,6,1,3,16,16/),(/5,2/))) - ELSE - CALL writeXMLElementForm('heAtomicEP',(/'atomType ','spin ','branch ',& - 'logDerivMT ','value '/),& - attributes(1:5),reshape((/10,4,6,12,5+17,6,1,3,16,16/),(/5,2/))) - ENDIF - WRITE (6,'(a7,i3,i2,a1,a12,f7.2,a4,f7.2,a5)') " Atom ",n,nqn,ch(l)," branch, D = ",& - ldmt, " at ",e," htr." - ENDIF - END FUNCTION priv_method2 - + END FUNCTION priv_method2 + END MODULE m_find_enpara diff --git a/hybrid/exchange_val_hf.F90 b/hybrid/exchange_val_hf.F90 index 05c3c17ef4f01d236478a66c7876758c8e2d2566..2ca46bce0f654777cc40bee187bfa5ef50fe8b29 100644 --- a/hybrid/exchange_val_hf.F90 +++ b/hybrid/exchange_val_hf.F90 @@ -50,7 +50,7 @@ ! converges well with q0. (Should be the default.) MODULE m_exchange_valence_hf - + use m_judft LOGICAL, PARAMETER:: zero_order = .false., ibs_corr = .false. INTEGER, PARAMETER:: maxmem = 600 @@ -287,7 +287,7 @@ CONTAINS hybrid%lcutm1, hybrid%maxlcutm1, hybrid%nindxm1, hybrid%maxindxm1, hybrid%gptm, & hybrid%ngptm(ikpt0), hybrid%pgptm(:, ikpt0), hybrid%gptmd, hybrid%basm1, & hybrid%nbasm(ikpt0), iband1, hybrid%nbands(nk), nsest, ibando, psize, indx_sest, & - atoms%invsat, sym%invsatnr, mpi%irank, cprod_vv_r(:hybrid%nbasm(ikpt0), :, :), & + sym%invsat, sym%invsatnr, mpi%irank, cprod_vv_r(:hybrid%nbasm(ikpt0), :, :), & cprod_vv_c(:hybrid%nbasm(ikpt0), :, :), mat_ex%l_real, wl_iks(:iband1, nkqpt), n_q(ikpt)) END IF #endif @@ -506,9 +506,10 @@ CONTAINS cdum = sqrt(expo)*rrad divergence = cell%omtil/(tpi_const**2)*sqrt(pi_const/expo)*cerf(cdum) rrad = rrad**2 - kv1 = cell%bmat(1, :)/kpts%nkpt3(1) - kv2 = cell%bmat(2, :)/kpts%nkpt3(2) - kv3 = cell%bmat(3, :)/kpts%nkpt3(3) + call judft_error("Missing functionality") + !kv1 = cell%bmat(1, :)/kpts%nkpt3(1) + !kv2 = cell%bmat(2, :)/kpts%nkpt3(2) + !kv3 = cell%bmat(3, :)/kpts%nkpt3(3) n = 1 found = .true. diff --git a/hybrid/hf_setup.F90 b/hybrid/hf_setup.F90 index 8cce209877a454a568f7f5a3f977000863b39b93..bb973cbee84802e94a4873003f76ad7f57310fdb 100644 --- a/hybrid/hf_setup.F90 +++ b/hybrid/hf_setup.F90 @@ -51,7 +51,7 @@ CONTAINS ! local arrays REAL, ALLOCATABLE :: basprod(:) - INTEGER :: degenerat(DIMENSION%neigd2 + 1, kpts%nkpt) + INTEGER :: degenerat(merge(dimension%neigd*2,dimension%neigd,noco%l_soc) + 1, kpts%nkpt) LOGICAL :: skip_kpt(kpts%nkpt) INTEGER :: g(3) @@ -63,7 +63,7 @@ CONTAINS ALLOCATE (zmat(kpts%nkptf), stat=ok) IF (ok /= 0) STOP 'eigen_hf: failure allocation z_c' - ALLOCATE (eig_irr(DIMENSION%neigd2, kpts%nkpt), stat=ok) + ALLOCATE (eig_irr(merge(dimension%neigd*2,dimension%neigd,noco%l_soc), kpts%nkpt), stat=ok) IF (ok /= 0) STOP 'eigen_hf: failure allocation eig_irr' ALLOCATE (hybdat%kveclo_eig(atoms%nlotot, kpts%nkpt), stat=ok) IF (ok /= 0) STOP 'eigen_hf: failure allocation hybdat%kveclo_eig' @@ -75,7 +75,7 @@ CONTAINS nrec1 = kpts%nkpt*(jsp - 1) + nk CALL lapw%init(input, noco, kpts, atoms, sym, nk, cell, sym%zrfs) nbasfcn = MERGE(lapw%nv(1) + lapw%nv(2) + 2*atoms%nlotot, lapw%nv(1) + atoms%nlotot, noco%l_noco) - CALL zMat(nk)%init(l_real, nbasfcn, dimension%neigd2) + CALL zMat(nk)%init(l_real, nbasfcn, merge(dimension%neigd*2,dimension%neigd,noco%l_soc)) CALL read_eig(eig_id_hf, nk, jsp, zmat=zMat(nk)) eig_irr(:, nk) = results%eig(:, nk, jsp) hybrid%ne_eig(nk) = results%neig(nk, jsp) @@ -83,7 +83,7 @@ CONTAINS !Allocate further space DO nk = kpts%nkpt + 1, kpts%nkptf nbasfcn = zMat(kpts%bkp(nk))%matsize1 - CALL zMat(nk)%init(l_real, nbasfcn, dimension%neigd2) + CALL zMat(nk)%init(l_real, nbasfcn, merge(dimension%neigd*2,dimension%neigd,noco%l_soc)) END DO !determine degenerate states at each k-point diff --git a/hybrid/hyb_abcrot.F90 b/hybrid/hyb_abcrot.F90 index fa39cc8556dea037e531cf8c9a8d0cbac87ba939..e2c9b594f265188b226fa01e825ebdb1826e19e1 100644 --- a/hybrid/hyb_abcrot.F90 +++ b/hybrid/hyb_abcrot.F90 @@ -52,7 +52,7 @@ CONTAINS DO ineq = 1, atoms%neq(itype) iatom = iatom + 1 IF (.NOT. oneD%odi%d1) THEN - iop = atoms%ngopr(iatom) + iop = sym%ngopr(iatom) ELSE iop = oneD%ods%ngopr(iatom) ENDIF @@ -60,9 +60,9 @@ CONTAINS ! inversion of spherical harmonics: Y (pi-theta,pi+phi) = (-1) * Y (theta,phi) ! m m ifac = 1 - IF (atoms%invsat(iatom) == 2) THEN + IF (sym%invsat(iatom) == 2) THEN IF (.NOT. oneD%odi%d1) THEN - iop = atoms%ngopr(sym%invsatnr(iatom)) + iop = sym%ngopr(sym%invsatnr(iatom)) ELSE iop = oneD%ods%ngopr(sym%invsatnr(iatom)) ENDIF diff --git a/hybrid/hybrid.F90 b/hybrid/hybrid.F90 index c0ef6ef4ab80e0947233ba6ac18c75f28244bee9..d7b119318808ccf68e912cc81ba8a23e60d5a865 100644 --- a/hybrid/hybrid.F90 +++ b/hybrid/hybrid.F90 @@ -88,7 +88,7 @@ CONTAINS END IF hybrid%l_subvxc = (hybrid%l_subvxc .AND. hybrid%l_addhf) - IF (.NOT. ALLOCATED(results%w_iks)) ALLOCATE (results%w_iks(DIMENSION%neigd2, kpts%nkpt, input%jspins)) + IF (.NOT. ALLOCATED(results%w_iks)) ALLOCATE (results%w_iks(merge(dimension%neigd*2,dimension%neigd,noco%l_soc), kpts%nkpt, input%jspins)) IF (hybrid%l_calhf) THEN iterHF = iterHF + 1 diff --git a/hybrid/kp_perturbation.F90 b/hybrid/kp_perturbation.F90 index 2aedc5056f0cdc6f2fcaba572c240e3584f21692..7c6244cb84e0f3bab9a1a7abcdb9de8826b28f8a 100644 --- a/hybrid/kp_perturbation.F90 +++ b/hybrid/kp_perturbation.F90 @@ -157,9 +157,9 @@ MODULE m_kp_perturbation const = fpi_const*(atoms%rmt(itype)**2)/2/sqrt(cell%omtil) DO ieq = 1, atoms%neq(itype) iatom = iatom + 1 - IF ((atoms%invsat(iatom) == 0) .or. (atoms%invsat(iatom) == 1)) THEN - IF (atoms%invsat(iatom) == 0) invsfct = 1 - IF (atoms%invsat(iatom) == 1) THEN + IF ((sym%invsat(iatom) == 0) .or. (sym%invsat(iatom) == 1)) THEN + IF (sym%invsat(iatom) == 0) invsfct = 1 + IF (sym%invsat(iatom) == 1) THEN invsfct = 2 iatom1 = sym%invsatnr(iatom) END IF diff --git a/hybrid/read_core.F90 b/hybrid/read_core.F90 index e8f07b946c57c3fce5de26f37658d5184bb68314..b94a28ec422099a0fe72bf6ff26f9592ed2df711 100644 --- a/hybrid/read_core.F90 +++ b/hybrid/read_core.F90 @@ -1,512 +1,549 @@ -MODULE m_read_core + MODULE m_read_core - ! read core radial wavefunctions from corebas - ! corebas is written in cored.F - ! (core basis functions can be read in once during an iteration) + ! read core radial wavefunctions from corebas + ! corebas is written in cored.F + ! (core basis functions can be read in once during an iteration) -CONTAINS + CONTAINS - SUBROUTINE read_core(atoms, hybdat) + SUBROUTINE read_core( atoms, hybdat ) - USE m_types + + USE m_types IMPLICIT NONE - TYPE(t_hybdat), INTENT(INOUT) :: hybdat - TYPE(t_atoms), INTENT(IN) :: atoms + TYPE(t_hybdat),INTENT(INOUT) :: hybdat + TYPE(t_atoms),INTENT(IN) :: atoms + ! - local scalars - INTEGER :: ncst - INTEGER :: ok, itype, i, idum, l + INTEGER :: ok,itype,i,idum,l REAL :: rdum - REAL :: weight1, weight2 + REAL :: weight1,weight2 ! - local arrays - - INTEGER, ALLOCATABLE :: nindxcr(:, :) - INTEGER, ALLOCATABLE :: l_qn(:, :), n_qn(:, :) + INTEGER,ALLOCATABLE :: nindxcr(:,:) + INTEGER,ALLOCATABLE :: l_qn(:,:),n_qn(:,:) - REAL, ALLOCATABLE :: j_qn(:, :) - REAL, ALLOCATABLE :: core1r(:, :, :, :),& - & core2r(:, :, :, :) + REAL,ALLOCATABLE :: j_qn(:,:) + REAL,ALLOCATABLE :: core1r(:,:,:,:),& + & core2r(:,:,:,:) - OPEN (UNIT=77, FILE='corebas', FORM='unformatted') - READ (77) ncst - ALLOCATE (n_qn(0:ncst, atoms%ntype), l_qn(0:ncst, atoms%ntype), & - & j_qn(0:ncst, atoms%ntype), stat=ok) - IF (ok /= 0) STOP 'mhsfock: failure allocation n_qn,l_qn,j_qn' - ALLOCATE (nindxcr(0:ncst, atoms%ntype), stat=ok) - IF (ok /= 0) STOP 'mhsfock: failure allocation nindxcr' + OPEN(UNIT=77,FILE='corebas',FORM='unformatted') + READ(77) ncst + + + ALLOCATE( n_qn(0:ncst,atoms%ntype) , l_qn(0:ncst,atoms%ntype) , & + & j_qn(0:ncst,atoms%ntype) , stat=ok ) + IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation n_qn,l_qn,j_qn' + ALLOCATE( nindxcr(0:ncst,atoms%ntype),stat = ok ) + IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation nindxcr' nindxcr = 0 - hybdat%lmaxc = 0 - l_qn = 0 - ncst = 0 - DO itype = 1, atoms%ntype - READ (77) ncst - IF (ncst == 0) CYCLE - DO i = 1, ncst - READ (77) n_qn(i, itype), l_qn(i, itype), j_qn(i, itype) - nindxcr(l_qn(i, itype), itype) = nindxcr(l_qn(i, itype), itype) + 1 - END DO - hybdat%lmaxc(itype) = maxval(l_qn(:, itype)) + hybdat%lmaxc = 0 + l_qn = 0 + ncst = 0 + DO itype=1,atoms%ntype + READ(77) ncst + IF( ncst .eq. 0 ) CYCLE + DO i=1,ncst + READ(77) n_qn(i,itype),l_qn(i,itype),j_qn(i,itype) + nindxcr(l_qn(i,itype),itype) = nindxcr(l_qn(i,itype),itype)+1 + END DO + hybdat%lmaxc(itype) = maxval(l_qn(:,itype)) END DO - ALLOCATE (core1r(atoms%jmtd, 0:maxval(hybdat%lmaxc), maxval(nindxcr), atoms%ntype)) - ALLOCATE (core2r(atoms%jmtd, 0:maxval(hybdat%lmaxc), maxval(nindxcr), atoms%ntype)) + + + ALLOCATE( core1r(atoms%jmtd,0:maxval(hybdat%lmaxc),maxval(nindxcr),atoms%ntype) ) + ALLOCATE( core2r(atoms%jmtd,0:maxval(hybdat%lmaxc),maxval(nindxcr),atoms%ntype) ) core1r = 0 core2r = 0 - REWIND (77) + REWIND(77) - DEALLOCATE (nindxcr) - ALLOCATE (nindxcr(0:maxval(l_qn), atoms%ntype)) + DEALLOCATE( nindxcr ) + ALLOCATE( nindxcr(0:maxval(l_qn),atoms%ntype)) nindxcr = 0 - READ (77) idum - DO itype = 1, atoms%ntype - READ (77) ncst - DO i = 1, ncst - nindxcr(l_qn(i, itype), itype) = nindxcr(l_qn(i, itype), itype) + 1 - READ (77) idum, idum, rdum, core1r(:atoms%jri(itype), l_qn(i, itype),& - & nindxcr(l_qn(i, itype), itype), itype),& - & core2r(:atoms%jri(itype), l_qn(i, itype),& - & nindxcr(l_qn(i, itype), itype), itype) - END DO + READ(77) idum + DO itype=1,atoms%ntype + READ(77) ncst + DO i=1,ncst + nindxcr(l_qn(i,itype),itype) = nindxcr(l_qn(i,itype),itype)+1 + READ(77) idum,idum,rdum,core1r(:atoms%jri(itype),l_qn(i,itype),& + & nindxcr(l_qn(i,itype),itype),itype),& + & core2r(:atoms%jri(itype),l_qn(i,itype),& + & nindxcr(l_qn(i,itype),itype),itype) + END DO END DO - ALLOCATE (hybdat%nindxc(0:maxval(hybdat%lmaxc), atoms%ntype), stat=ok) - IF (ok /= 0) STOP 'mhsfock: failure allocation nindxc' - ALLOCATE (hybdat%core1(atoms%jmtd, 0:maxval(hybdat%lmaxc), maxval(nindxcr(0, :),& - & nint((maxval(nindxcr)/2.0))), atoms%ntype), stat=ok) - IF (ok /= 0) STOP 'mhsfock: failure allocation core1' - ALLOCATE (hybdat%core2(atoms%jmtd, 0:maxval(hybdat%lmaxc), maxval(nindxcr(0, :),& - & nint((maxval(nindxcr)/2.0))), atoms%ntype), stat=ok) - IF (ok /= 0) STOP 'mhsfock: failure allocation core2' - hybdat%nindxc = 0; hybdat%core1 = 0; hybdat%core2 = 0 + + + ALLOCATE( hybdat%nindxc(0:maxval(hybdat%lmaxc),atoms%ntype),stat=ok ) + IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation nindxc' + ALLOCATE( hybdat%core1(atoms%jmtd,0:maxval(hybdat%lmaxc),maxval(nindxcr(0,:),& + & nint((maxval(nindxcr)/2.0))),atoms%ntype),stat=ok ) + IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation core1' + ALLOCATE( hybdat%core2(atoms%jmtd,0:maxval(hybdat%lmaxc),maxval(nindxcr(0,:),& + & nint((maxval(nindxcr)/2.0))),atoms%ntype),stat=ok ) + IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation core2' + hybdat%nindxc = 0 ; hybdat%core1 = 0 ; hybdat%core2 = 0 ! average over core states that only differ in j ! core functions with l-qn equal 0 doesnot change during averaging - hybdat%nindxc(0, :) = nindxcr(0, :) - DO itype = 1, atoms%ntype - hybdat%core1(:, 0, :hybdat%nindxc(0, itype), itype)& - & = core1r(:, 0, :hybdat%nindxc(0, itype), itype) - hybdat%core2(:, 0, :hybdat%nindxc(0, itype), itype)& - & = core2r(:, 0, :hybdat%nindxc(0, itype), itype) + hybdat%nindxc(0,:) = nindxcr(0,:) + DO itype=1,atoms%ntype + hybdat%core1(:,0,:hybdat%nindxc(0,itype),itype)& + & = core1r(:,0,:hybdat%nindxc(0,itype),itype) + hybdat%core2(:,0,:hybdat%nindxc(0,itype),itype)& + & = core2r(:,0,:hybdat%nindxc(0,itype),itype) END DO - DO itype = 1, atoms%ntype - DO l = 1, hybdat%lmaxc(itype) - weight1 = 2*(l - 0.5) + 1 - weight2 = 2*(l + 0.5) + 1 - IF (modulo(nindxcr(l, itype), 2) == 0) THEN - DO i = 1, nindxcr(l, itype), 2 - hybdat%nindxc(l, itype) = hybdat%nindxc(l, itype) + 1 - hybdat%core1(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype) =& - & (weight1*core1r(:atoms%jri(itype), l, i, itype) +& - & weight2*core1r(:atoms%jri(itype), l, i + 1, itype))& - & /(weight1 + weight2) - hybdat%core2(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype) =& - & (weight1*core2r(:atoms%jri(itype), l, i, itype) +& - & weight2*core2r(:atoms%jri(itype), l, i + 1, itype))& - & /(weight1 + weight2) - END DO - ELSE - DO i = 1, nindxcr(l, itype) - 1, 2 - hybdat%nindxc(l, itype) = hybdat%nindxc(l, itype) + 1 - hybdat%core1(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype) =& - & (weight1*core1r(:atoms%jri(itype), l, i, itype) +& - & weight2*core1r(:atoms%jri(itype), l, i + 1, itype))& - & /(weight1 + weight2) - hybdat%core2(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype) =& - & (weight1*core2r(:atoms%jri(itype), l, i, itype) +& - & weight2*core2r(:atoms%jri(itype), l, i + 1, itype))& - & /(weight1 + weight2) - END DO - hybdat%nindxc(l, itype) = hybdat%nindxc(l, itype) + 1 - hybdat%core1(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype)& - & = core1r(:atoms%jri(itype), l, nindxcr(l, itype), itype) - hybdat%core2(:atoms%jri(itype), l, hybdat%nindxc(l, itype), itype)& - & = core2r(:atoms%jri(itype), l, nindxcr(l, itype), itype) - END IF - END DO + DO itype=1,atoms%ntype + DO l=1,hybdat%lmaxc(itype) + weight1 = 2*(l-0.5) + 1 + weight2 = 2*(l+0.5) + 1 + IF( modulo(nindxcr(l,itype),2) .eq. 0 ) THEN + DO i=1,nindxcr(l,itype),2 + hybdat%nindxc(l,itype) = hybdat%nindxc(l,itype) + 1 + hybdat%core1(:atoms%jri(itype),l,hybdat%nindxc(l,itype),itype) =& + & (weight1*core1r(:atoms%jri(itype),l,i ,itype) +& + & weight2*core1r(:atoms%jri(itype),l,i+1,itype))& + & / (weight1+weight2) + hybdat%core2(:atoms%jri(itype),l,hybdat%nindxc(l,itype),itype) =& + & (weight1*core2r(:atoms%jri(itype),l,i ,itype) +& + & weight2*core2r(:atoms%jri(itype),l,i+1,itype))& + & / (weight1+weight2) + END DO + ELSE + DO i=1,nindxcr(l,itype)-1,2 + hybdat%nindxc(l,itype) = hybdat%nindxc(l,itype) + 1 + hybdat%core1(:atoms%jri(itype),l,hybdat%nindxc(l,itype),itype) =& + & (weight1*core1r(:atoms%jri(itype),l,i ,itype) +& + & weight2*core1r(:atoms%jri(itype),l,i+1,itype))& + & / (weight1+weight2) + hybdat%core2(:atoms%jri(itype),l,hybdat%nindxc(l,itype),itype) =& + & (weight1*core2r(:atoms%jri(itype),l,i ,itype) +& + & weight2*core2r(:atoms%jri(itype),l,i+1,itype))& + & / (weight1+weight2) + END DO + hybdat%nindxc(l,itype) = hybdat%nindxc(l,itype) + 1 + hybdat%core1(:atoms%jri(itype),l,hybdat%nindxc(l,itype),itype)& + & = core1r(:atoms%jri(itype),l,nindxcr(l,itype),itype) + hybdat%core2(:atoms%jri(itype),l,hybdat%nindxc(l,itype),itype)& + & = core2r(:atoms%jri(itype),l,nindxcr(l,itype),itype) + END IF + + END DO END DO - DEALLOCATE (nindxcr, core1r, core2r) + DEALLOCATE( nindxcr, core1r,core2r ) hybdat%maxindxc = maxval(hybdat%nindxc) - CLOSE (77) + CLOSE(77) + + END SUBROUTINE read_core + - END SUBROUTINE read_core - SUBROUTINE corewf(atoms, jsp, input, dimension,& - & vr, lmaxcd, maxindxc, mpi, lmaxc, nindxc, core1, core2, eig_c) + + + + + SUBROUTINE corewf( atoms,jsp,input,dimension,& + & vr, lmaxcd,maxindxc,mpi, lmaxc,nindxc,core1,core2,eig_c) USE m_types IMPLICIT NONE - TYPE(t_mpi), INTENT(IN) :: mpi - TYPE(t_dimension), INTENT(IN) :: dimension - TYPE(t_input), INTENT(IN) :: input - TYPE(t_atoms), INTENT(IN) :: atoms + TYPE(t_mpi),INTENT(IN) :: mpi + TYPE(t_dimension),INTENT(IN) :: dimension + TYPE(t_input),INTENT(IN) :: input + TYPE(t_atoms),INTENT(IN) :: atoms ! - scalars - - INTEGER, INTENT(IN) :: jsp - INTEGER, INTENT(IN) :: lmaxcd - INTEGER, INTENT(INOUT) :: maxindxc + INTEGER,INTENT(IN) :: jsp + INTEGER,INTENT(IN) :: lmaxcd + INTEGER,INTENT(INOUT) :: maxindxc ! -arrays - - INTEGER, INTENT(INOUT) :: lmaxc(:) - INTEGER, INTENT(INOUT) :: nindxc(0:lmaxcd, atoms%ntype) + INTEGER,INTENT(INOUT) :: lmaxc(:) + INTEGER,INTENT(INOUT) :: nindxc(0:lmaxcd,atoms%ntype) - REAL, INTENT(IN) :: vr(:, :, :)!(atoms%jmtd,atoms%ntypd,input%jspins) - REAL, INTENT(INOUT) :: core1(:, :, 0:, :) !(atoms%jmtd,maxindxc,0:lmaxcd,atoms%ntype) - REAL, INTENT(INOUT) :: core2(:, :, 0:, :) !(jmtd,maxindxc,0:lmaxcd,ntype) + REAL , INTENT(IN) :: vr(:,:,:)!(atoms%jmtd,atoms%ntypd,input%jspins) + REAL , INTENT(INOUT) :: core1(:,:,0:,:) !(atoms%jmtd,maxindxc,0:lmaxcd,atoms%ntype) + REAL , INTENT(INOUT) :: core2(:,:,0:,:) !(jmtd,maxindxc,0:lmaxcd,ntype) - REAL, INTENT(INOUT) :: eig_c(maxindxc, 0:lmaxcd, atoms%ntype) + REAL , INTENT(INOUT) :: eig_c(maxindxc,0:lmaxcd,atoms%ntype) ! - local scalars - INTEGER :: ncstd - INTEGER :: ok, itype, i, j, idum, l + INTEGER :: ok,itype,i,j,idum,l REAL :: rdum - REAL :: weight1, weight2 + REAL :: weight1,weight2 ! - local arrays - - INTEGER, ALLOCATABLE :: nindxcr(:, :) - INTEGER, ALLOCATABLE :: l_qn(:, :), n_qn(:, :) + INTEGER,ALLOCATABLE :: nindxcr(:,:) + INTEGER,ALLOCATABLE :: l_qn(:,:),n_qn(:,:) - REAL, ALLOCATABLE :: j_qn(:, :) - REAL, ALLOCATABLE :: core1r(:, :, :, :), core2r(:, :, :, :) - REAL, ALLOCATABLE :: core11r(:, :, :, :), core22r(:, :, :, :) - REAL, ALLOCATABLE :: eig_cr(:, :, :) + REAL,ALLOCATABLE :: j_qn(:,:) + REAL,ALLOCATABLE :: core1r(:,:,:,:),core2r(:,:,:,:) + REAL,ALLOCATABLE :: core11r(:,:,:,:),core22r(:,:,:,:) + REAL,ALLOCATABLE :: eig_cr(:,:,:) - ncstd = maxval(atoms%ncst) - ALLOCATE (nindxcr(0:ncstd, atoms%ntype), stat=ok) + ncstd = maxval(atoms%econf%num_core_states) + ALLOCATE( nindxcr(0:ncstd,atoms%ntype),stat = ok ) ! generate relativistic core wave functions( ->core1r,core2r ) - CALL calcorewf(dimension, input, jsp, atoms,& - & ncstd, vr,& - & lmaxc, nindxcr, core1r, core2r, eig_cr, mpi) + CALL calcorewf( dimension,input,jsp,atoms,& + & ncstd,vr,& + & lmaxc,nindxcr,core1r,core2r,eig_cr,mpi) + nindxc = 0 ! average over core states that only differ in j ! core functions with l-qn equal 0 doesnot change during the average process - - nindxc(0, :) = nindxcr(0, :) - DO itype = 1, atoms%ntype - core1(:, :nindxc(0, itype), 0, itype)& - & = core1r(:, 0, :nindxc(0, itype), itype) - core2(:, :nindxc(0, itype), 0, itype)& - & = core2r(:, 0, :nindxc(0, itype), itype) - eig_c(:nindxc(0, itype), 0, itype)& - & = eig_cr(0, :nindxc(0, itype), itype) + + nindxc(0,:) = nindxcr(0,:) + DO itype=1,atoms%ntype + core1(:,:nindxc(0,itype),0,itype)& + & = core1r(:,0,:nindxc(0,itype),itype) + core2(:,:nindxc(0,itype),0,itype)& + & = core2r(:,0,:nindxc(0,itype),itype) + eig_c(:nindxc(0,itype),0,itype)& + & = eig_cr(0,:nindxc(0,itype),itype) END DO - DO itype = 1, atoms%ntype - DO l = 1, lmaxc(itype) - weight1 = 2*(l - 0.5) + 1 - weight2 = 2*(l + 0.5) + 1 - IF (modulo(nindxcr(l, itype), 2) == 0) THEN - DO i = 1, nindxcr(l, itype), 2 - nindxc(l, itype) = nindxc(l, itype) + 1 - core1(:atoms%jri(itype), nindxc(l, itype), l, itype) =& - & (weight1*core1r(:atoms%jri(itype), l, i, itype) +& - & weight2*core1r(:atoms%jri(itype), l, i + 1, itype))& - & /(weight1 + weight2) - core2(:atoms%jri(itype), nindxc(l, itype), l, itype) =& - & (weight1*core2r(:atoms%jri(itype), l, i, itype) + & - & weight2*core2r(:atoms%jri(itype), l, i + 1, itype))& - & /(weight1 + weight2) - - eig_c(nindxc(l, itype), l, itype) =& - & (weight1*eig_cr(l, i, itype) +& - & weight2*eig_cr(l, i + 1, itype))& - & /(weight1 + weight2) - END DO - ELSE - DO i = 1, nindxcr(l, itype) - 1, 2 - nindxc(l, itype) = nindxc(l, itype) + 1 - core1(:atoms%jri(itype), nindxc(l, itype), l, itype) =& - & (weight1*core1r(:atoms%jri(itype), l, i, itype) +& - & weight2*core1r(:atoms%jri(itype), l, i + 1, itype))& - & /(weight1 + weight2) - core2(:atoms%jri(itype), nindxc(l, itype), l, itype) =& - & (weight1*core2r(:atoms%jri(itype), l, i, itype) + & - & weight2*core2r(:atoms%jri(itype), l, i + 1, itype))& - & /(weight1 + weight2) - - eig_c(nindxc(l, itype), l, itype) =& - & (weight1*eig_cr(l, i, itype) +& - & weight2*eig_cr(l, i + 1, itype))& - & /(weight1 + weight2) - END DO - nindxc(l, itype) = nindxc(l, itype) + 1 - core1(:atoms%jri(itype), nindxc(l, itype), l, itype)& - & = core1r(:atoms%jri(itype), l, nindxcr(l, itype), itype) - core2(:atoms%jri(itype), nindxc(l, itype), l, itype)& - & = core2r(:atoms%jri(itype), l, nindxcr(l, itype), itype) - eig_c(nindxc(l, itype), l, itype)& - & = eig_cr(l, nindxcr(l, itype), itype) - END IF - END DO + DO itype=1,atoms%ntype + DO l=1,lmaxc(itype) + weight1 = 2*(l-0.5) + 1 + weight2 = 2*(l+0.5) + 1 + IF( modulo(nindxcr(l,itype),2) .eq. 0 ) THEN + DO i=1,nindxcr(l,itype),2 + nindxc(l,itype) = nindxc(l,itype) + 1 + core1(:atoms%jri(itype),nindxc(l,itype),l,itype) =& + & (weight1*core1r(:atoms%jri(itype),l,i ,itype) +& + & weight2*core1r(:atoms%jri(itype),l,i+1,itype))& + & / (weight1+weight2) + core2(:atoms%jri(itype),nindxc(l,itype),l,itype) =& + & (weight1*core2r(:atoms%jri(itype),l,i ,itype) + & + & weight2*core2r(:atoms%jri(itype),l,i+1,itype))& + & / (weight1+weight2) + + eig_c(nindxc(l,itype),l,itype) =& + & (weight1*eig_cr(l,i,itype) +& + & weight2*eig_cr(l,i+1,itype))& + & / (weight1+weight2) + END DO + ELSE + DO i=1,nindxcr(l,itype)-1,2 + nindxc(l,itype) = nindxc(l,itype) + 1 + core1(:atoms%jri(itype),nindxc(l,itype),l,itype) =& + & (weight1*core1r(:atoms%jri(itype),l,i ,itype) +& + & weight2*core1r(:atoms%jri(itype),l,i+1,itype))& + & / (weight1+weight2) + core2(:atoms%jri(itype),nindxc(l,itype),l,itype) =& + & (weight1*core2r(:atoms%jri(itype),l,i ,itype) + & + & weight2*core2r(:atoms%jri(itype),l,i+1,itype))& + & / (weight1+weight2) + + eig_c(nindxc(l,itype),l,itype) =& + & (weight1*eig_cr(l,i,itype) +& + & weight2*eig_cr(l,i+1,itype))& + & / (weight1+weight2) + END DO + nindxc(l,itype) = nindxc(l,itype) + 1 + core1(:atoms%jri(itype),nindxc(l,itype),l,itype)& + & = core1r(:atoms%jri(itype),l,nindxcr(l,itype),itype) + core2(:atoms%jri(itype),nindxc(l,itype),l,itype)& + & = core2r(:atoms%jri(itype),l,nindxcr(l,itype),itype) + eig_c(nindxc(l,itype),l,itype)& + & = eig_cr(l,nindxcr(l,itype),itype) + END IF + + END DO END DO - DEALLOCATE (nindxcr, core1r, core2r, eig_cr) + DEALLOCATE( nindxcr,core1r,core2r,eig_cr ) - IF (maxindxc /= maxval(nindxc))& + IF( maxindxc .ne. maxval(nindxc) )& & STOP 'corewf: counting error nindxc' - END SUBROUTINE corewf - SUBROUTINE calcorewf(dimension, input, jspin, atoms,& - & ncstd, vr,& - & lmaxc, nindxcr, core1, core2, eig_c, mpi) + END SUBROUTINE corewf + - USE m_intgr, ONLY: intgr3, intgr0, intgr1 - USE m_constants, ONLY: c_light - USE m_setcor + SUBROUTINE calcorewf(dimension,input,jspin,atoms,& + & ncstd,vr,& + & lmaxc,nindxcr,core1,core2,eig_c,mpi) + + + USE m_intgr, ONLY : intgr3,intgr0,intgr1 + USE m_constants, ONLY : c_light + !USE m_setcor USE m_differ USE m_types IMPLICIT NONE - TYPE(t_mpi), INTENT(IN) :: mpi - TYPE(t_dimension), INTENT(IN) :: dimension - TYPE(t_input), INTENT(IN) :: input - TYPE(t_atoms), INTENT(IN) :: atoms + TYPE(t_mpi),INTENT(IN) :: mpi + TYPE(t_dimension),INTENT(IN) :: dimension + TYPE(t_input),INTENT(IN) :: input + TYPE(t_atoms),INTENT(IN) :: atoms ! - scalars - - INTEGER, INTENT(IN) :: ncstd - INTEGER, INTENT(IN) :: jspin - INTEGER, INTENT(OUT):: lmaxc(:) + INTEGER, INTENT (IN) :: ncstd + INTEGER, INTENT (IN) :: jspin + INTEGER, INTENT (OUT):: lmaxc(:) - ! - arrays - - INTEGER, INTENT(OUT):: nindxcr(0:ncstd, atoms%ntype) - REAL, INTENT(IN) :: vr(:, :, :)!(atoms%jmtd,atoms%ntypd,input%jspins) - REAL, ALLOCATABLE :: core1(:, :, :, :), core2(:, :, :, :) - REAL, ALLOCATABLE :: eig_c(:, :, :) + ! - arrays - + INTEGER, INTENT (OUT):: nindxcr(0:ncstd,atoms%ntype) + REAL , INTENT (IN) :: vr(:,:,:)!(atoms%jmtd,atoms%ntypd,input%jspins) + REAL , ALLOCATABLE :: core1(:,:,:,:),core2(:,:,:,:) + REAL , ALLOCATABLE :: eig_c(:,:,:) ! - local scalars - - INTEGER :: i, j, itype, korb, ncmsh, nst, ierr - REAL :: e, fj, fl, fn, t2, c, bmu, weight - REAL :: d, dxx, rn, rnot, z, t1, rr + INTEGER :: i,j,itype,korb,ncmsh,nst,ierr + REAL :: e,fj,fl,fn,t2,c,bmu,weight + REAL :: d,dxx,rn,rnot,z,t1,rr LOGICAL, SAVE :: first = .true. ! - local arrays - - INTEGER :: kappa(dimension%nstd), nprnc(dimension%nstd) - REAL :: vrd(dimension%msh) - REAL :: occ(dimension%nstd), occ_h(dimension%nstd, 2), a(dimension%msh), b(dimension%msh) - REAL, ALLOCATABLE, SAVE:: vr0(:, :, :) + INTEGER :: kappa(29),nprnc(29) + REAL :: vrd(atoms%msh) + REAL :: occ(29),occ_h(29,2),a(atoms%msh),b(atoms%msh) + REAL,ALLOCATABLE,SAVE:: vr0(:,:,:) ! - intrinsic functions - - INTRINSIC exp, iabs, isign + INTRINSIC exp,iabs,isign c = c_light(1.0) - IF (first) THEN - ALLOCATE (vr0(atoms%jmtd, atoms%ntype, input%jspins)) + IF( first ) THEN + ALLOCATE( vr0(atoms%jmtd,atoms%ntype,input%jspins) ) END IF - IF (input%frcor) THEN - IF (first) THEN - vr0 = vr - first = .false. - END IF + IF ( input%frcor ) THEN + IF ( first ) THEN + vr0 = vr + first = .false. + END IF ELSE - vr0 = vr - first = .false. + vr0 = vr + first = .false. END IF + ! this loop determines the dimensions lmaxc = 0; nindxcr = 0 - DO itype = 1, atoms%ntype - z = atoms%zatom(itype) - dxx = atoms%dx(itype) - bmu = 0.0 - CALL setcor(itype, input%jspins, atoms, input, bmu,& - & nst, kappa, nprnc, occ_h) - - IF ((bmu > 99.)) THEN - occ(1:nst) = input%jspins*occ_h(1:nst, jspin) - ELSE - occ(1:nst) = occ_h(1:nst, 1) - END IF - rnot = atoms%rmsh(1, itype) - d = exp(atoms%dx(itype)) - ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1) - ncmsh = min(ncmsh, dimension%msh) - rn = rnot*(d**(ncmsh - 1)) - - nst = atoms%ncst(itype) - - DO 80 korb = 1, nst - IF (occ(korb) == 0) GOTO 80 + DO itype = 1,atoms%ntype + z = atoms%zatom(itype) + dxx = atoms%dx(itype) + bmu = 0.0 +! CALL setcor( itype,input%jspins,atoms,input,bmu,& +! & nst,kappa,nprnc,occ_h) + call atoms%econf(itype)%get_core(nst,kappa,nprnc,occ_h) + + + IF ((bmu > 99.)) THEN + occ(1:nst) = input%jspins * occ_h(1:nst,jspin) + ELSE + occ(1:nst) = occ_h(1:nst,1) + END IF + rnot = atoms%rmsh(1,itype) + d = exp(atoms%dx(itype)) + ncmsh = nint( log( (atoms%rmt(itype)+10.0)/rnot ) / dxx + 1 ) + ncmsh = min( ncmsh, atoms%msh ) + rn = rnot* (d** (ncmsh-1)) + + nst = atoms%econf(itype)%num_core_states + + DO 80 korb = 1,nst + IF (occ(korb).EQ.0) GOTO 80 fn = nprnc(korb) fj = iabs(kappa(korb)) - .5e0 weight = 2*fj + 1.e0 IF (bmu > 99.) weight = occ(korb) - fl = fj + (.5e0)*isign(1, kappa(korb)) - e = -2*(z/(fn + fl))**2 + fl = fj + (.5e0)*isign(1,kappa(korb)) + e = -2* (z/ (fn+fl))**2 - nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1 - lmaxc(itype) = max(lmaxc(itype), NINT(fl)) -80 END DO + nindxcr(NINT(fl),itype) = nindxcr(NINT(fl),itype) + 1 + lmaxc(itype) = max(lmaxc(itype),NINT(fl)) + 80 END DO END DO - ALLOCATE (core1(atoms%jmtd, 0:maxval(lmaxc), maxval(nindxcr), atoms%ntype)) - ALLOCATE (core2(atoms%jmtd, 0:maxval(lmaxc), maxval(nindxcr), atoms%ntype)) - ALLOCATE (eig_c(0:maxval(lmaxc), maxval(nindxcr), atoms%ntype)) - core1 = 0; core2 = 0 - nindxcr = 0 - DO itype = 1, atoms%ntype - z = atoms%zatom(itype) - dxx = atoms%dx(itype) - bmu = 0.0 - CALL setcor(itype, input%jspins, atoms, input, bmu, nst, kappa, nprnc, occ_h) - - IF ((bmu > 99.)) THEN - occ(1:nst) = input%jspins*occ_h(1:nst, jspin) - ELSE - occ(1:nst) = occ_h(1:nst, 1) - END IF - rnot = atoms%rmsh(1, itype) - d = exp(atoms%dx(itype)) - ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1) - ncmsh = min(ncmsh, dimension%msh) - rn = rnot*(d**(ncmsh - 1)) - IF (mpi%irank == 0) THEN - WRITE (6, FMT=8000) z, rnot, dxx, atoms%jri(itype) - END IF - DO j = 1, atoms%jri(itype) - vrd(j) = vr0(j, itype, jspin) - END DO - if (input%l_core_confpot) THEN + ALLOCATE( core1(atoms%jmtd,0:maxval(lmaxc),maxval(nindxcr),atoms%ntype) ) + ALLOCATE( core2(atoms%jmtd,0:maxval(lmaxc),maxval(nindxcr),atoms%ntype) ) + ALLOCATE( eig_c(0:maxval(lmaxc),maxval(nindxcr),atoms%ntype) ) - ! linear extension of the potential with slope t1 / a.u. - t1 = 0.125 - t2 = vrd(atoms%jri(itype))/atoms%rmt(itype) - atoms%rmt(itype)*t1 - rr = atoms%rmt(itype) - else - t2 = vrd(atoms%jri(itype))/(atoms%jri(itype) - ncmsh) + core1 = 0; core2 =0 + nindxcr = 0 + DO itype = 1,atoms%ntype + z = atoms%zatom(itype) + dxx = atoms%dx(itype) + bmu = 0.0 + !CALL setcor(itype,input%jspins,atoms,input,bmu,nst,kappa,nprnc,occ_h) + call atoms%econf(itype)%get_core(nst,nprnc,kappa,occ_h) + + IF ((bmu > 99.)) THEN + occ(1:nst) = input%jspins * occ_h(1:nst,jspin) + ELSE + occ(1:nst) = occ_h(1:nst,1) + END IF + rnot = atoms%rmsh(1,itype) + d = exp(atoms%dx(itype)) + ncmsh = nint( log( (atoms%rmt(itype)+10.0)/rnot ) / dxx + 1 ) + ncmsh = min( ncmsh, atoms%msh ) + rn = rnot* (d** (ncmsh-1)) + IF ( mpi%irank == 0 ) THEN + WRITE(6 ,FMT=8000) z,rnot,dxx,atoms%jri(itype) + END IF + DO j = 1,atoms%jri(itype) + vrd(j) = vr0(j,itype,jspin) + END DO + + if (input%l_core_confpot) THEN + + ! linear extension of the potential with slope t1 / a.u. + t1=0.125 + t2=vrd(atoms%jri(itype))/atoms%rmt(itype)-atoms%rmt(itype)*t1 + rr = atoms%rmt(itype) + else + t2 = vrd(atoms%jri(itype)) / ( atoms%jri(itype) - ncmsh ) + endif + IF ( atoms%jri(itype) .LT. ncmsh) THEN + DO i = atoms%jri(itype) + 1,ncmsh + if (input%l_core_confpot) THEN + rr = d*rr + vrd(i) = rr*( t2 + rr*t1 ) + else + vrd(i) = vrd(atoms%jri(itype)) + t2* (i-atoms%jri(itype)) endif - IF (atoms%jri(itype) < ncmsh) THEN - DO i = atoms%jri(itype) + 1, ncmsh - if (input%l_core_confpot) THEN - rr = d*rr - vrd(i) = rr*(t2 + rr*t1) - else - vrd(i) = vrd(atoms%jri(itype)) + t2*(i - atoms%jri(itype)) - endif ! - END DO - END IF + END DO + END IF + + nst = atoms%econf(itype)%num_core_states - nst = atoms%ncst(itype) - DO 90 korb = 1, nst - IF (occ(korb) == 0) GOTO 90 + DO 90 korb = 1,nst + IF (occ(korb).EQ.0) GOTO 90 fn = nprnc(korb) fj = iabs(kappa(korb)) - .5e0 weight = 2*fj + 1.e0 IF (bmu > 99.) weight = occ(korb) - fl = fj + (.5e0)*isign(1, kappa(korb)) - e = -2*(z/(fn + fl))**2 - CALL differ(fn, fl, fj, c, z, dxx, rnot, rn, d, ncmsh, vrd, e, a, b, ierr) + fl = fj + (.5e0)*isign(1,kappa(korb)) + e = -2* (z/ (fn+fl))**2 + CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, e, a,b,ierr) - nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1 + nindxcr(NINT(fl),itype) = nindxcr(NINT(fl),itype) + 1 - core1(:atoms%jri(itype), NINT(fl), nindxcr(NINT(fl), itype), itype)& - & = a(:atoms%jri(itype)) - core2(:atoms%jri(itype), NINT(fl), nindxcr(NINT(fl), itype), itype)& - & = b(:atoms%jri(itype)) - eig_c(NINT(fl), nindxcr(NINT(fl), itype), itype) = e + core1(:atoms%jri(itype),NINT(fl),nindxcr(NINT(fl),itype),itype)& + & = a(:atoms%jri(itype)) + core2(:atoms%jri(itype),NINT(fl),nindxcr(NINT(fl),itype),itype)& + & = b(:atoms%jri(itype)) - IF (mpi%irank == 0) THEN - WRITE (6, FMT=8010) fn, fl, fj, e, weight + eig_c(NINT(fl),nindxcr(NINT(fl),itype),itype) = e + + IF ( mpi%irank == 0 ) THEN + WRITE (6,FMT=8010) fn,fl,fj,e,weight END IF - IF (ierr /= 0) STOP 'error in core-level routine' + IF (ierr.NE.0) STOP 'error in core-level routine' + -90 END DO + 90 END DO END DO -8000 FORMAT(/, /, 10x, 'z=', f4.0, 5x, 'r(1)=', e14.6, 5x, 'dx=', f8.6, 5x,& - & 'm.t.index=', i4, /, 15x, 'n', 4x, 'l', 5x, 'j', 4x, 'energy', 7x,& - & 'weight') -8010 FORMAT(12x, 2f5.0, f6.1, f10.4, f10.0) - END SUBROUTINE calcorewf + 8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f8.6,5x,& + & 'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,& + & 'weight') + 8010 FORMAT (12x,2f5.0,f6.1,f10.4,f10.0) + + END SUBROUTINE calcorewf - SUBROUTINE core_init(dimension, input, atoms, lmaxcd, maxindxc) + SUBROUTINE core_init( dimension,input,atoms, lmaxcd,maxindxc) - USE m_intgr, ONLY: intgr3, intgr0, intgr1 - USE m_constants, ONLY: c_light - USE m_setcor + + USE m_intgr, ONLY : intgr3,intgr0,intgr1 + USE m_constants, ONLY : c_light + !USE m_setcor USE m_differ USE m_types IMPLICIT NONE - TYPE(t_dimension), INTENT(IN) :: dimension - TYPE(t_input), INTENT(IN) :: input - TYPE(t_atoms), INTENT(IN) :: atoms - INTEGER, INTENT(OUT) :: maxindxc, lmaxcd + TYPE(t_dimension),INTENT(IN) :: dimension + TYPE(t_input),INTENT(IN) :: input + TYPE(t_atoms),INTENT(IN) :: atoms + INTEGER,INTENT(OUT) :: maxindxc,lmaxcd + ! - local scalars - - INTEGER :: i, j, itype, korb, ncmsh, nst, ierr - REAL :: e, fj, fl, fn, t, bmu, c - REAL :: d, dxx, rn, rnot, z, t1, rr + INTEGER :: i,j,itype,korb,ncmsh,nst,ierr + REAL :: e,fj,fl,fn,t,bmu,c + REAL :: d,dxx,rn,rnot,z,t1,rr ! - local arrays - - INTEGER :: kappa(dimension%nstd), nprnc(dimension%nstd) - INTEGER :: nindxcr(0:dimension%nstd, atoms%ntype) - REAL :: occ(dimension%nstd), occ_h(dimension%nstd, 2), a(dimension%msh), b(dimension%msh) + INTEGER :: kappa(29),nprnc(29) + INTEGER :: nindxcr(0:29,atoms%ntype) + REAL :: occ(29),occ_h(29,2),a(atoms%msh),b(atoms%msh) INTEGER :: lmaxc(atoms%ntype) + ! - intrinsic functions - - INTRINSIC exp, iabs, isign + INTRINSIC exp,iabs,isign c = c_light(1.0) + ! this loop determines the dimensions lmaxc = 0; nindxcr = 0 - DO itype = 1, atoms%ntype - z = atoms%zatom(itype) - dxx = atoms%dx(itype) - bmu = 0.0 - CALL setcor(itype, input%jspins, atoms, input, bmu, nst, kappa, nprnc, occ_h) + DO itype = 1,atoms%ntype + z = atoms%zatom(itype) + dxx = atoms%dx(itype) + bmu = 0.0 + !CALL setcor(itype,input%jspins,atoms,input,bmu, nst,kappa,nprnc,occ_h) + call atoms%econf(itype)%get_core(nst,nprnc,kappa,occ_h) - occ(1:nst) = occ_h(1:nst, 1) + occ(1:nst) = occ_h(1:nst,1) - rnot = atoms%rmsh(1, itype) - d = exp(atoms%dx(itype)) - ncmsh = nint(log((atoms%rmt(itype) + 10.0)/rnot)/dxx + 1) - ncmsh = min(ncmsh, dimension%msh) - rn = rnot*(d**(ncmsh - 1)) + rnot = atoms%rmsh(1,itype) + d = exp(atoms%dx(itype)) + ncmsh = nint( log( (atoms%rmt(itype)+10.0)/rnot ) / dxx + 1 ) + ncmsh = min( ncmsh, atoms%msh ) + rn = rnot* (d** (ncmsh-1)) - nst = atoms%ncst(itype) + nst = atoms%econf(itype)%num_core_states - DO korb = 1, nst - IF (occ(korb) == 0) CYCLE + DO korb = 1,nst + IF (occ(korb).EQ.0) CYCLE fn = nprnc(korb) fj = iabs(kappa(korb)) - .5e0 - fl = fj + (.5e0)*isign(1, kappa(korb)) - e = -2*(z/(fn + fl))**2 + fl = fj + (.5e0)*isign(1,kappa(korb)) + e = -2* (z/ (fn+fl))**2 - nindxcr(NINT(fl), itype) = nindxcr(NINT(fl), itype) + 1 - lmaxc(itype) = max(lmaxc(itype), NINT(fl)) + nindxcr(NINT(fl),itype) = nindxcr(NINT(fl),itype) + 1 + lmaxc(itype) = max(lmaxc(itype),NINT(fl)) END DO END DO - lmaxcd = maxval(lmaxc) - maxindxc = maxval(nindxcr(0, :), nint((maxval(nindxcr)/2.0))) + lmaxcd = maxval(lmaxc) + maxindxc = maxval(nindxcr(0,:),nint((maxval(nindxcr)/2.0))) + + - END SUBROUTINE core_init + END SUBROUTINE core_init -END MODULE m_read_core + END MODULE m_read_core diff --git a/hybrid/subvxc.F90 b/hybrid/subvxc.F90 index 2b79a0ca11ad3c28b95b9df179f35eaaa7f9f3b8..d49747fddc11bf655ef7eedace9c7f4041e5305f 100644 --- a/hybrid/subvxc.F90 +++ b/hybrid/subvxc.F90 @@ -168,7 +168,7 @@ CONTAINS iatom = 0 DO itype = 1, atoms%ntype - typsym = atoms%ntypsy(SUM(atoms%neq(:itype - 1)) + 1) + typsym = sym%ntypsy(SUM(atoms%neq(:itype - 1)) + 1) nlharm = sphhar%nlh(typsym) ! Calculate vxc = vtot - vcoul @@ -298,7 +298,7 @@ CONTAINS DO itype = 1, atoms%ntype - typsym = atoms%ntypsy(SUM(atoms%neq(:itype - 1)) + 1) + typsym = sym%ntypsy(SUM(atoms%neq(:itype - 1)) + 1) nlharm = sphhar%nlh(typsym) ! Calculate vxc = vtot - vcoul @@ -340,10 +340,10 @@ CONTAINS DO ieq = 1, atoms%neq(itype) iatom = iatom + 1 - IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN + IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN - IF (atoms%invsat(iatom) == 0) invsfct = 1 - IF (atoms%invsat(iatom) == 1) invsfct = 2 + IF (sym%invsat(iatom) == 0) invsfct = 1 + IF (sym%invsat(iatom) == 1) invsfct = 2 DO ilo = 1, atoms%nlo(itype) #ifdef CPP_OLDINTEL @@ -484,7 +484,7 @@ CONTAINS END DO ! ilo ikvecprevat = ikvecprevat + ikvecat ikvecat = 0 - END IF ! atoms%invsat(iatom) + END IF ! sym%invsat(iatom) END DO ! ieq END DO !itype END IF ! if any atoms%llo @@ -513,4 +513,3 @@ CONTAINS END SUBROUTINE subvxc END MODULE m_subvxc - diff --git a/hybrid/symm_hf.F90 b/hybrid/symm_hf.F90 index 8bbe41925f20835f8e1c81af7897c10fc7a1a50c..30bd63cd5bf6add79d73bd2cb0ec25903f06d7c6 100644 --- a/hybrid/symm_hf.F90 +++ b/hybrid/symm_hf.F90 @@ -10,7 +10,7 @@ ! M.Betzinger (09/07) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MODULE m_symm_hf - +use m_judft #define irreps .false. CONTAINS @@ -52,7 +52,8 @@ CONTAINS rotkpt = matmul(rrot(:, :, i), kpts%bkf(:, nk)) !transfer rotkpt into BZ - rotkpt = modulo1(rotkpt, kpts%nkpt3) + call judft_error("Missing functionality") + !rotkpt = modulo1(rotkpt, kpts%nkpt3) !check if rotkpt is identical to bk(:,nk) IF (maxval(abs(rotkpt - kpts%bkf(:, nk))) <= 1E-07) THEN @@ -162,7 +163,8 @@ CONTAINS rotkpt = matmul(rrot(:, :, psym(iop)), kpts%bkf(:, ikpt)) !transfer rotkpt into BZ - rotkpt = modulo1(rotkpt, kpts%nkpt3) + call judft_error("Missing functionality") + !rotkpt = modulo1(rotkpt, kpts%nkpt3) !determine number of rotkpt nrkpt = 0 @@ -221,7 +223,8 @@ CONTAINS rotkpt = matmul(rrot(:, :, isym), kpts%bkf(:, ikpt)) !transfer rotkpt into BZ - rotkpt = modulo1(rotkpt, kpts%nkpt3) + call judft_error("Missing functionality") + !rotkpt = modulo1(rotkpt, kpts%nkpt3) !check if rotkpt is identical to bk(:,ikpt) IF (maxval(abs(rotkpt - kpts%bkf(:, ikpt))) <= 1E-06) THEN @@ -573,7 +576,8 @@ CONTAINS rotkpt = matmul(rrot(:, :, iop), kpts%bkf(:, nk)) !transfer rotkpt into BZ - rotkpt = modulo1(rotkpt, kpts%nkpt3) + call judft_error("Missing functionality") + !rotkpt = modulo1(rotkpt, kpts%nkpt3) !check if rotkpt is identical to bk(:,nk) IF (maxval(abs(rotkpt - kpts%bkf(:, nk))) <= 1E-07) THEN @@ -608,7 +612,8 @@ CONTAINS rotkpt = matmul(rrot(:, :, psym(iop)), kpts%bkf(:, ikpt)) !transfer rotkpt into BZ - rotkpt = modulo1(rotkpt, kpts%nkpt3) + call judft_error("Missing functionality") + !rotkpt = modulo1(rotkpt, kpts%nkpt3) !determine number of rotkpt nrkpt = 0 diff --git a/hybrid/symmetrizeh.F90 b/hybrid/symmetrizeh.F90 index 42c4ef1df57986854a999fba3244c3a2bc881da1..48832e54e0ab506cb832ffaed3e8a2f8117d4014 100644 --- a/hybrid/symmetrizeh.F90 +++ b/hybrid/symmetrizeh.F90 @@ -233,9 +233,9 @@ CONTAINS DO itype = 1, atoms%ntype DO ieq = 1, atoms%neq(itype) iatom = iatom + 1 - IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN - IF (atoms%invsat(iatom) == 0) invsfct = 1 - IF (atoms%invsat(iatom) == 1) invsfct = 2 + IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN + IF (sym%invsat(iatom) == 0) invsfct = 1 + IF (sym%invsat(iatom) == 1) invsfct = 2 DO ilo = 1, atoms%nlo(itype) l = atoms%llo(ilo, itype) DO m = 1, invsfct*(2*l + 1) @@ -264,9 +264,9 @@ CONTAINS DO itype = 1, atoms%ntype DO ieq = 1, atoms%neq(itype) iatom = iatom + 1 - IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN - IF (atoms%invsat(iatom) == 0) invsfct = 1 - IF (atoms%invsat(iatom) == 1) invsfct = 2 + IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN + IF (sym%invsat(iatom) == 0) invsfct = 1 + IF (sym%invsat(iatom) == 1) invsfct = 2 DO ilo = 1, atoms%nlo(itype) l = atoms%llo(ilo, itype) @@ -325,10 +325,10 @@ CONTAINS iatom = iatom + 1 ratom = map(isym, iatom) - IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN - IF (atoms%invsat(iatom) == 0) invsfct = 1 - IF (atoms%invsat(iatom) == 1) THEN - IF (atoms%invsat(ratom) == 2) THEN + IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN + IF (sym%invsat(iatom) == 0) invsfct = 1 + IF (sym%invsat(iatom) == 1) THEN + IF (sym%invsat(ratom) == 2) THEN ratom = sym%invsatnr(ratom) END IF invsfct = 2 @@ -387,9 +387,9 @@ CONTAINS DO itype = 1, atoms%ntype DO ieq = 1, atoms%neq(itype) iatom = iatom + 1 - IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN - IF (atoms%invsat(iatom) == 0) invsfct = 1 - IF (atoms%invsat(iatom) == 1) invsfct = 2 + IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN + IF (sym%invsat(iatom) == 0) invsfct = 1 + IF (sym%invsat(iatom) == 1) invsfct = 2 DO ilo = 1, atoms%nlo(itype) l = atoms%llo(ilo, itype) @@ -403,7 +403,7 @@ CONTAINS iop = psym(isym) ratom = map(isym, iatom) IF (invsfct == 2) THEN - IF (atoms%invsat(ratom) == 2) THEN + IF (sym%invsat(ratom) == 2) THEN ratom = sym%invsatnr(ratom) END IF END IF @@ -448,9 +448,9 @@ CONTAINS DO itype = 1, atoms%ntype DO ieq = 1, atoms%neq(itype) iatom = iatom + 1 - IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN - IF (atoms%invsat(iatom) == 0) invsfct = 1 - IF (atoms%invsat(iatom) == 1) invsfct = 2 + IF ((sym%invsat(iatom) == 0) .OR. (sym%invsat(iatom) == 1)) THEN + IF (sym%invsat(iatom) == 0) invsfct = 1 + IF (sym%invsat(iatom) == 1) invsfct = 2 DO ilo = 1, atoms%nlo(itype) l = atoms%llo(ilo, itype) @@ -461,9 +461,9 @@ CONTAINS DO itype1 = 1, atoms%ntype DO ieq1 = 1, atoms%neq(itype1) iatom1 = iatom1 + 1 - IF ((atoms%invsat(iatom1) == 0) .OR. (atoms%invsat(iatom1) == 1)) THEN - IF (atoms%invsat(iatom1) == 0) invsfct1 = 1 - IF (atoms%invsat(iatom1) == 1) invsfct1 = 2 + IF ((sym%invsat(iatom1) == 0) .OR. (sym%invsat(iatom1) == 1)) THEN + IF (sym%invsat(iatom1) == 0) invsfct1 = 1 + IF (sym%invsat(iatom1) == 1) invsfct1 = 2 DO ilo1 = 1, atoms%nlo(itype1) l1 = atoms%llo(ilo1, itype1) @@ -478,12 +478,12 @@ CONTAINS ratom1 = map(isym, iatom1) IF (invsfct == 2) THEN - IF (atoms%invsat(ratom) == 2) THEN + IF (sym%invsat(ratom) == 2) THEN ratom = sym%invsatnr(ratom) END IF END IF IF (invsfct1 == 2) THEN - IF (atoms%invsat(ratom1) == 2) THEN + IF (sym%invsat(ratom1) == 2) THEN ratom1 = sym%invsatnr(ratom1) END IF END IF diff --git a/hybrid/trafo.F90 b/hybrid/trafo.F90 index 1f03d0fa5c57a4cce4570930b8bba083b5306f2c..ae53a6580f142f46ea7d3ee564b6454e807a1352 100644 --- a/hybrid/trafo.F90 +++ b/hybrid/trafo.F90 @@ -5,7 +5,7 @@ !-------------------------------------------------------------------------------- MODULE m_trafo - + use m_judft CONTAINS SUBROUTINE waveftrafo_symm(cmt_out, z_out, cmt, l_real, z_r, z_c, bandi, ndb, & @@ -76,7 +76,8 @@ CONTAINS rkpt = matmul(rrot, kpts%bk(:, nk)) rkpthlp = rkpt - rkpt = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + !rkpt = modulo1(rkpt, kpts%nkpt3) g1 = nint(rkpt - rkpthlp) ! MT coefficients @@ -223,7 +224,8 @@ CONTAINS rkpt = matmul(rrot, kpts%bk(:, nk)) rkpthlp = rkpt - rkpt = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + !rkpt = modulo1(rkpt, kpts%nkpt3) g1 = nint(rkpt - rkpthlp) ! MT coefficients @@ -369,7 +371,7 @@ CONTAINS nn = sum((/((2*l + 1)*nindxm(l, itype), l=0, lcutm(itype))/)) DO ieq = 1, atoms%neq(itype) ic = ic + 1 - IF (atoms%invsat(ic) == 0) THEN + IF (sym%invsat(ic) == 0) THEN ! if the structure is inversion-symmetric, but the equivalent atom belongs to a different unit cell ! invsat(atom) = 0, invsatnr(atom) = 0 ! but we need invsatnr(atom) = natom @@ -481,7 +483,7 @@ CONTAINS nn = sum((/((2*l + 1)*nindxm(l, itype), l=0, lcutm(itype))/)) DO ieq = 1, atoms%neq(itype) ic = ic + 1 - IF (atoms%invsat(ic) == 0) THEN + IF (sym%invsat(ic) == 0) THEN ! if the structure is inversion-symmetric, but the equivalent atom belongs to a different unit cell ! invsat(atom) = 0, invsatnr(atom) =0 ! but we need invsatnr(atom) = natom @@ -633,7 +635,8 @@ CONTAINS rkpt = matmul(rrot, kpts%bkf(:, ikpt0)) rkpthlp = rkpt - rkpt = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + !rkpt = modulo1(rkpt, kpts%nkpt3) g = nint(rkpthlp - rkpt) #ifdef CPP_DEBUG @@ -838,7 +841,8 @@ CONTAINS rrot = transpose(sym%mrot(:, :, sym%invtab(iisym))) invrrot = transpose(sym%mrot(:, :, iisym)) rkpt = matmul(rrot, kpts%bk(:, ikpt0)) - rkpthlp = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + ! rkpthlp = modulo1(rkpt, kpts%nkpt3) g = nint(rkpt - rkpthlp) CALL d_wigner(invrot, cell%bmat, maxlcutm, dwgn(:, :, 1:maxlcutm)) @@ -854,7 +858,8 @@ CONTAINS rrot = -transpose(sym%mrot(:, :, sym%invtab(iisym))) invrrot = -transpose(sym%mrot(:, :, iisym)) rkpt = matmul(rrot, kpts%bk(:, ikpt0)) - rkpthlp = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + !rkpthlp = modulo1(rkpt, kpts%nkpt3) g = nint(rkpt - rkpthlp) matin1 = conjg(matin1) @@ -1111,7 +1116,8 @@ CONTAINS END DO rkpt = matmul(rrot, kpts%bk(:, ikpt0)) - rkpthlp = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + !rkpthlp = modulo1(rkpt, kpts%nkpt3) g = nint(rkpt - rkpthlp) ! determine number of rotated k-point bk(:,ikpt) -> ikpt1 @@ -1332,7 +1338,8 @@ CONTAINS rrot = transpose(sym%mrot(:, :, sym%invtab(iisym))) invrrot = transpose(sym%mrot(:, :, iisym)) rkpt = matmul(rrot, kpts%bk(:, ikpt0)) - rkpthlp = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + !rkpthlp = modulo1(rkpt, kpts%nkpt3) g = nint(rkpt - rkpthlp) CALL d_wigner(invrot, cell%bmat, maxlcutm, dwgn(:, :, 1:maxlcutm)) @@ -1348,7 +1355,8 @@ CONTAINS invrot = sym%mrot(:, :, sym%invtab(iisym)) invrrot = -transpose(sym%mrot(:, :, iisym)) rkpt = matmul(rrot, kpts%bk(:, ikpt0)) - rkpthlp = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + !rkpthlp = modulo1(rkpt, kpts%nkpt3) g = nint(rkpt - rkpthlp) vecin1 = conjg(vecin1) @@ -1534,7 +1542,8 @@ CONTAINS END IF rkpt = matmul(rrot, kpts%bk(:, ikpt0)) - rkpthlp = modulo1(rkpt, kpts%nkpt3) + call judft_error("Missing functionality here") + !rkpthlp = modulo1(rkpt, kpts%nkpt3) g = nint(rkpt - rkpthlp) DO l = 0, maxlcutm @@ -1716,7 +1725,8 @@ CONTAINS END IF rkpthlp = matmul(rrot, kpts%bk(:, ikpt0)) - rkpt = modulo1(rkpthlp, kpts%nkpt3) + call judft_error("Missing functionality here") + ! rkpt = modulo1(rkpthlp, kpts%nkpt3) g = nint(rkpthlp - rkpt) ! ! determine number of rotated k-point bk(:,ikpt) -> ikpt1 diff --git a/hybrid/wavefproducts.F90 b/hybrid/wavefproducts.F90 index 361c702a406d7797a127cc56ba855ec4df5221b1..9cc8f8a6ddb9fdf4680bd97b89840e8a73fd2361 100644 --- a/hybrid/wavefproducts.F90 +++ b/hybrid/wavefproducts.F90 @@ -94,7 +94,8 @@ CONTAINS kqpthlp = kpts%bkf(:, nk) + kpts%bkf(:, iq) ! k+q can lie outside the first BZ, transfer ! it back into the 1. BZ - kqpt = modulo1(kqpthlp, kpts%nkpt3) + call judft_error("Missing functionality") + !kqpt = modulo1(kqpthlp, kpts%nkpt3) g_t(:) = nint(kqpt - kqpthlp) ! determine number of kqpt nkqpt = 0 @@ -368,7 +369,8 @@ CONTAINS kqpthlp = kpts%bkf(:, nk) + kpts%bkf(:, iq) ! kqpt can lie outside the first BZ, transfer it back - kqpt = modulo1(kqpthlp, kpts%nkpt3) + call judft_error("Missing functionality") + !kqpt = modulo1(kqpthlp, kpts%nkpt3) g_t(:) = nint(kqpt - kqpthlp) ! determine number of kqpt nkqpt = 0 @@ -1293,7 +1295,8 @@ CONTAINS kqpthlp = kpts%bkf(:, nk) + kpts%bkf(:, iq) ! kqpt can lie outside the first BZ, transfer it back - kqpt = modulo1(kqpthlp, kpts%nkpt3) + call judft_error("Missing functionality") + !kqpt = modulo1(kqpthlp, kpts%nkpt3) g_t(:) = nint(kqpt - kqpthlp) ! determine number of kqpt nkqpt = 0 @@ -2305,7 +2308,8 @@ CONTAINS kqpthlp = kpts%bkf(:, nk) + kpts%bkf(:, iq) ! k+q can lie outside the first BZ, transfer ! it back into the 1. BZ - kqpt = modulo1(kqpthlp, kpts%nkpt3) + call judft_error("Missing functionality") + !kqpt = modulo1(kqpthlp, kpts%nkpt3) g_t(:) = nint(kqpt - kqpthlp) ! determine number of kqpt nkqpt = 0 diff --git a/inpgen2/old_inp/dimen7.F90 b/inpgen2/old_inp/dimen7.F90 index e195dbe6e9b49c4677ff0aafb63939f16722b819..7716854beaf770b322f59350a14f7bc063bfb7f7 100644 --- a/inpgen2/old_inp/dimen7.F90 +++ b/inpgen2/old_inp/dimen7.F90 @@ -8,7 +8,7 @@ ! ! This program reads the input files of the flapw-programm (inp & kpts) -! and creates a file 'fl7para' that contains dimensions +! and creates a file 'fl7para' that contains dimensions ! for the main flapw-programm. ! @@ -50,7 +50,7 @@ ! TYPE(t_input),INTENT(INOUT) :: input TYPE(t_sym),INTENT(INOUT) :: sym - TYPE(t_stars),INTENT(INOUT) :: stars + TYPE(t_stars),INTENT(INOUT) :: stars TYPE(t_atoms),INTENT(INOUT) :: atoms TYPE(t_sphhar),INTENT(INOUT) :: sphhar TYPE(t_dimension),INTENT(INOUT) :: dimension @@ -59,7 +59,7 @@ TYPE(t_oneD),INTENT(INOUT) :: oneD TYPE(t_hybrid),INTENT(INOUT) :: hybrid TYPE(t_cell),INTENT(INOUT) :: cell - + TYPE(t_noco) :: noco TYPE(t_sliceplot) :: sliceplot TYPE(t_banddos) :: banddos @@ -69,23 +69,23 @@ ! !------------------------------------------------------------------- ! .. Local Scalars .. - REAL :: thetad,xa,epsdisp,epsforce ,rmtmax,arltv1,arltv2,arltv3 + REAL :: thetad,xa,epsdisp,epsforce ,rmtmax,arltv1,arltv2,arltv3 REAL :: s,r,d ,idsprs INTEGER :: ok,ilo,n,nstate,i,j,na,n1,n2,jrc,nopd,symfh INTEGER :: nmopq(3) - + CHARACTER(len=1) :: rw - CHARACTER(len=4) :: namex + CHARACTER(len=4) :: namex CHARACTER(len=7) :: symfn CHARACTER(len=12):: relcor LOGICAL ::l_kpts,l_qpts,l_inpexist,l_tmp(2) ! .. - REAL :: a1(3),a2(3),a3(3) + REAL :: a1(3),a2(3),a3(3) REAL :: q(3) CHARACTER(len=3), ALLOCATABLE :: noel(:) - LOGICAL, ALLOCATABLE :: error(:) - + LOGICAL, ALLOCATABLE :: error(:) + INTEGER ntp1,ii,grid(3) INTEGER, ALLOCATABLE :: lmx1(:), nq1(:), nlhtp1(:) @@ -95,8 +95,8 @@ real :: scalecell EXTERNAL prp_xcfft_box!,parawrite ! .. - - + + !---> First, check whether an inp-file exists ! INQUIRE (file='inp',exist=l_inpexist) @@ -125,7 +125,7 @@ & hybrid%select1(4,atoms%ntype),hybrid%lcutm1(atoms%ntype),& & hybrid%lcutwf(atoms%ntype), STAT=ok) ! -!---> read complete input and calculate nvacd,llod,lmaxd,jmtd,neigd and +!---> read complete input and calculate nvacd,llod,lmaxd,jmtd,neigd and ! CALL rw_inp('r',& & atoms,vacuum,input,stars,sliceplot,banddos,& @@ -136,7 +136,7 @@ !---> read the angle and spin-spiral information from nocoinp noco%qss = 0.0 noco%l_ss = .false. - IF (noco%l_noco) THEN + IF (noco%l_noco) THEN CALL inpnoco(atoms,input,vacuum,noco) ENDIF @@ -147,7 +147,6 @@ atoms%jmtd = 0 rmtmax = 0.0 dimension%neigd = 0 - 29 = maxval(atoms%econf%num_core_states) atoms%lmaxd = maxval(atoms%lmax) atoms%jmtd = maxval(atoms%jri) rmtmax = maxval(atoms%rmt) @@ -178,8 +177,8 @@ DIMENSION%neigd = MAX(5,NINT(0.75*input%zelec)) WRITE(6,*) 'changed dimension%neigd to ',dimension%neigd ENDIF - IF (noco%l_soc .and. (.not. noco%l_noco)) dimension%neigd=2*dimension%neigd - IF (noco%l_soc .and. noco%l_ss) dimension%neigd=(3*dimension%neigd)/2 + IF (noco%l_soc .and. (.not. noco%l_noco)) dimension%neigd=2*dimension%neigd + IF (noco%l_soc .and. noco%l_ss) dimension%neigd=(3*dimension%neigd)/2 ! not as accurate, but saves much time rmtmax = rmtmax*stars%gmax @@ -195,7 +194,7 @@ jrc = jrc + 1 r = r*d ENDDO - atoms%msh = max( atoms%msh, jrc ) + atoms%msh = max( atoms%msh, jrc ) ENDDO ! ! ---> now, set the lattice harmonics, determine nlhd @@ -209,7 +208,7 @@ IF (oneD%odd%d1) cell%omtil = cell%amat(3,3)*pimach()*(vacuum%dvac**2)/4. !+odim cell%bmat=tpi_const*cell%bmat - + na = 0 DO n = 1,atoms%ntype DO n1 = 1,atoms%neq(n) @@ -250,7 +249,7 @@ atoms%nat,ntp1,nq1,cell%amat,cell%bmat,atoms%taual,& sphhar%nlhd,sphhar%memd,sphhar%ntypsd,.true.,nlhtp1,& sphhar%nlh,sphhar%llh,sphhar%nmem,& - sphhar%mlh,sphhar%clnu) + sphhar%mlh,sphhar%clnu) ii = 1 DO i = 1,atoms%ntype atoms%nlhtyp(i) = nlhtp1(ii) @@ -262,23 +261,23 @@ ! ! Check if symmetry is compatible with SOC or SSDW ! - IF (noco%l_soc .and. (.not.noco%l_noco)) THEN + IF (noco%l_soc .and. (.not.noco%l_noco)) THEN ! test symmetry for spin-orbit coupling ALLOCATE ( error(sym%nop) ) CALL soc_sym(sym%nop,sym%mrot,noco%theta,noco%phi,cell%amat,error) IF ( ANY(error(:)) ) THEN WRITE(*,fmt='(1x)') - WRITE(*,fmt='(A)') 'Symmetry incompatible with SOC spin-quantization axis ,' - WRITE(*,fmt='(A)') 'do not perform self-consistent calculations !' + WRITE(*,fmt='(A)') 'Symmetry incompatible with SOC spin-quantization axis ,' + WRITE(*,fmt='(A)') 'do not perform self-consistent calculations !' WRITE(*,fmt='(1x)') IF ( input%eonly .or. (noco%l_soc.and.noco%l_ss) .or. input%gw.ne.0 ) THEN ! .or. . - CONTINUE - ELSE + CONTINUE + ELSE IF (input%itmax>1) THEN CALL juDFT_error("symmetry & SOC",calledby ="dimen7") - ENDIF - ENDIF - ENDIF + ENDIF + ENDIF + ENDIF DEALLOCATE ( error ) ENDIF IF (noco%l_ss) THEN ! test symmetry for spin-spiral @@ -344,7 +343,7 @@ ENDIF ENDIF ENDIF - + dimension%neigd = max(dimension%neigd,input%gw_neigd) ! @@ -364,7 +363,6 @@ ! CALL inpeig_dim(input,cell,noco,oneD,kpts,dimension,stars,latnam) vacuum%layerd = max(vacuum%layerd,1) - 29 = max(29,30) atoms%ntype = atoms%ntype IF (noco%l_noco) dimension%neigd = 2*dimension%neigd diff --git a/inpgen2/old_inp/dimens.F90 b/inpgen2/old_inp/dimens.F90 index 8ebca6ba5cb5777ec4f0b5a903220bac2ade6f88..6a72e5b61fc47c43d33385b82251cc7a10d26380 100644 --- a/inpgen2/old_inp/dimens.F90 +++ b/inpgen2/old_inp/dimens.F90 @@ -1,4 +1,4 @@ -pd-------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------- ! 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. @@ -30,7 +30,7 @@ CONTAINS IMPLICIT NONE TYPE(t_input),INTENT(INOUT) :: input TYPE(t_sym),INTENT(INOUT) :: sym - TYPE(t_stars),INTENT(INOUT) :: stars + TYPE(t_stars),INTENT(INOUT) :: stars TYPE(t_atoms),INTENT(INOUT) :: atoms TYPE(t_sphhar),INTENT(INOUT) :: sphhar TYPE(t_dimension),INTENT(INOUT) :: dimension @@ -38,7 +38,7 @@ CONTAINS TYPE(t_kpts),INTENT(INOUT) :: kpts TYPE(t_oneD),INTENT(INOUT) :: oneD TYPE(t_hybrid),INTENT(INOUT) :: hybrid - + TYPE(t_cell) :: cell LOGICAL l_kpts,l_qpts,l_inpexist,ldum diff --git a/inpgen2/old_inp/rw_inp.f90 b/inpgen2/old_inp/rw_inp.f90 index 63e0b1ee19b63c9a5edc7c686808960024333907..686b6b562453d68cafc388f873fed425aa80968c 100644 --- a/inpgen2/old_inp/rw_inp.f90 +++ b/inpgen2/old_inp/rw_inp.f90 @@ -32,7 +32,7 @@ USE m_types_xcpot_inbuild_nofunction USE m_types_noco - + IMPLICIT NONE ! .. ! .. Arguments .. @@ -40,7 +40,7 @@ TYPE(t_input),INTENT(INOUT) :: input TYPE(t_sym),INTENT(INOUT) :: sym - TYPE(t_stars),INTENT(INOUT) :: stars + TYPE(t_stars),INTENT(INOUT) :: stars TYPE(t_atoms),INTENT(INOUT) :: atoms TYPE(t_vacuum),INTENT(INOUT) :: vacuum TYPE(t_kpts),INTENT(INOUT) :: kpts @@ -51,10 +51,10 @@ TYPE(t_sliceplot),INTENT(INOUT):: sliceplot TYPE(t_xcpot_inbuild_nf),INTENT(INOUT) :: xcpot TYPE(t_noco),INTENT(INOUT) :: noco - + REAL,INTENT(INOUT) :: a1(3),a2(3),a3(3),scalecell CHARACTER(len=3),INTENT(OUT) :: noel(atoms%ntype) - CHARACTER(len=4),INTENT(OUT) :: namex + CHARACTER(len=4),INTENT(OUT) :: namex CHARACTER(len=12),INTENT(OUT):: relcor CHARACTER(len=*),INTENT(INOUT)::latnam,namgrp INTEGER,INTENT(OUT)::grid(3) @@ -80,7 +80,7 @@ !-odim ! .. ! .. Local Variables - REAL ::scpos ,zc,dtild + REAL ::scpos ,zc,dtild INTEGER ::nw,idsprs,ncst INTEGER ieq,i,k,na,n,ilo REAL s3,ah,a,hs2,rest @@ -107,14 +107,14 @@ IF (ch_rw.eq.'r') THEN !-------------------------------------------------------------------- OPEN (5,file='inp',form='formatted',status='old') - + !default not read in in old inp-file input%qfix=2 ! a1(:) = 0 a2(:) = 0 a3(:) = 0 - + WRITE (6,*) '-------- dump of inp-file ------------' ! @@ -259,7 +259,7 @@ a2(2) = -a1(2) END IF - + IF (namgrp.EQ.'any ') THEN INQUIRE (file='sym.out',exist=l_sym) IF (.not.l_sym)& @@ -301,11 +301,11 @@ & (namex.EQ.'Rpbe').OR.(namex.EQ.'wc') .OR.& & (namex.EQ.'pbe0').OR.(namex.EQ.'hse ').OR.& & (namex.EQ.'lhse').OR.(namex.EQ.'vhse')) THEN ! some defaults - idsprs=0 + idsprs=0 ENDIF ! set mixing and screening for variable HSE functional WRITE (6,9040) namex,relcor - + ! look what comes in the next two lines ! READ (UNIT=5,FMT=7182,END=77,ERR=77) ch_test @@ -420,7 +420,7 @@ 7161 FORMAT (i2,8x,l1,5x,i2,5x,60i3) WRITE (6,9091) atoms%neq(n),atoms%l_geo(n),atoms%nlo(n),& & (atoms%llo(ilo,n),ilo=1,atoms%nlo(n)) - END IF + END IF ! DO ieq=1,atoms%neq(n) na = na + 1 @@ -484,7 +484,7 @@ IF(ldum) CALL juDFT_error& & ("QGpsi exists but gw /= 2 in inp.",calledby ="rw_inp") ENDIF - + BACKSPACE(5) ! Make sure that input%vchk,input%cdinf,obsolete%pot8 are all given. READ (UNIT=5,FMT=7220,END=99,ERR=99) input%vchk,input%cdinf,ldum if (ldum) call judft_error("pot8 not longer supported") @@ -492,12 +492,12 @@ 7220 FORMAT (5x,l1,1x,6x,l1,1x,5x,l1,1x,3x,i1,1x,9x,i4) ! DO i=1,100 ; line(i:i)=' ' ; ENDDO - - input%eig66(2)=.false. - + + !input%eig66(2)=.false. + READ (UNIT=5,FMT=6000,END=99,ERR=99)& - & idum,ldum,input%l_f,input%eonly,input%eig66(1)!,input%eig66(2) - WRITE (6,9130) 0,.false.,input%l_f,input%eonly,input%eig66(1)!,input%eig66(2) + & idum,ldum,input%l_f,input%eonly!,input%eig66(1)!,input%eig66(2) + WRITE (6,9130) 0,.false.,input%l_f,input%eonly!,input%eig66(1)!,input%eig66(2) 6000 FORMAT (4x,i1,8x,l1,5x,l1,7x,l1,7x,l1) ! !+roa @@ -523,7 +523,7 @@ ! READ (UNIT=5,FMT='(f10.5)',END=99,ERR=99) input%rkmax WRITE (6,FMT='(f10.5,1x,A)') input%rkmax, '=kmax' - + READ (UNIT=5,FMT=8010,END=99,ERR=99) input%gauss,input%delgau,input%tria WRITE (6,9160) input%gauss,input%delgau,input%tria 8010 FORMAT (6x,l1,f10.5,5x,l1) @@ -559,7 +559,7 @@ !!$ chform= '(40x,l1,1x,'//chntype//'a1)' !!$ CALL judft_error("soc_opt no longer supported") !!$ ENDIF - + READ (UNIT=5,FMT=8050,END=99,ERR=99)& & input%frcor,sliceplot%slice,input%ctail input%coretail_lmax=99 @@ -569,7 +569,7 @@ input%l_bmt= ( line(52:56)=='bmt=T' ).or.( line(52:56)=='bmt=t' ) WRITE (6,9170) input%frcor,sliceplot%slice,input%ctail 8050 FORMAT (6x,l1,7x,l1,7x,l1,6x,l1,7x,i1,5x,l1,5x,l1) - + ! check if itmax consists of 2 or 3 digits READ(unit=5,FMT='(8x,a)') check BACKSPACE 5 @@ -585,8 +585,8 @@ WRITE (6,9180) input%itmax,input%maxiter,input%imix,input%alpha,input%spinf 8061 FORMAT (6x,i3,9x,i3,6x,i2,7x,f6.2,7x,f6.2) END IF - - input%preconditioning_param = 0.0 + + input%preconditioning_param = 0.0 chform = '(5x,l1,'//chntype//'f6.2)' ! chform = '(5x,l1,23f6.2)' @@ -646,7 +646,7 @@ ! READ (UNIT=5,FMT=8090,END=99,ERR=99) ! !input%xa,input%thetad,input%epsdisp,input%epsforce - WRITE (6,*) "No relaxation with old input anymore" + WRITE (6,*) "No relaxation with old input anymore" !input%xa,input%thetad,input%epsdisp,input%epsforce 8090 FORMAT (3x,f10.5,8x,f10.5,9x,f10.5,10x,f10.5) ! @@ -665,7 +665,7 @@ READ (UNIT=5,FMT='(9x,f10.5,10x,f10.5,9x,f10.5)',& & END=98,ERR=98) banddos%e2_dos,banddos%e1_dos,banddos%sig_dos - + ! added for exact-exchange or hybrid functional calculations: ! read in the number of k-points and nx,ny and nz given in the last line ! of the input file, @@ -682,7 +682,7 @@ WRITE(*,*) '' CALL juDFT_error("Invalid declaration of k-point set (1)",calledby="rw_inp") END IF - + !IF( kpts%nkpt3(1)*kpts%nkpt3(2)*kpts%nkpt3(3) .ne. idum ) THEN ! WRITE(*,*) '' ! WRITE(*,*) 'nx*ny*nz is not equal to nkpt.' @@ -704,7 +704,7 @@ !END IF=98,ERR=98) hybrid%lcutm2(i),hybrid%select2(1,i),hybrid%select2(2,i),& ! hybrid%select2(3,i),hybrid%select2(4,i) !END DO - + !ALLOCATE( hybrid%l_exxc(maxval(atoms%ncst),atoms%ntype) ) !DO i=1,atoms%ntype ! READ(UNIT=5,FMT='(60(2x,l1))',END=98,ERR=98)(hybrid%l_exxc(k,i),k=1,atoms%ncst(i)) @@ -904,7 +904,7 @@ WRITE (5,9140) 1,0 WRITE (5,'(a)') 'ellow, elup, valence electrons:' - + WRITE (5,9150) input%ellow,input%elup,input%zelec 9150 FORMAT (4f10.5) WRITE (5,fmt='(f10.5,1x,A)') input%rkmax, '=kmax' diff --git a/kpoints/kptmop.f b/kpoints/kptmop.f index aa75dd2d24606a5e599544070101bb4c447d2aa0..5dc3c92c382f91311b5841146b73d10820d3823d 100644 --- a/kpoints/kptmop.f +++ b/kpoints/kptmop.f @@ -1,79 +1,70 @@ MODULE m_kptmop use m_juDFT !----------------------------------------------------------------------- -! ---> This program generates k-points -! in irreducible wedge of BZ (for nreg=0) -! in total BZ (for nreg=1) -! (BZ = 1. Brillouin-zone) for all canonical Bravais lattices -! in 2 and 3 dimensions, -! using the basis vectors of the reciprocal lattice -! and the bordering planes of the irreducible wedge. -! -! The k-points are generated by the Monkhorst--Pack method. -! The information on the bordering planes of the irr wedge -! is taken from BRZONE. -! -! The program checks the compatibility of the dimension and -! symmetry and of the provided Monkhorst-Pack-parameters. +! ---> This program generates k-points +! in irreducible wedge of BZ +! (BZ = 1. Brillouin-zone) for all canonical Bravais lattices +! in 2 and 3 dimensions, +! using the basis vectors of the reciprocal lattice +! and the bordering planes of the irreducible wedge. +! +! The k-points are generated by the Monkhorst--Pack method. +! The information on the bordering planes of the irr wedge +! is taken from BRZONE. +! +! The program checks the compatibility of the dimension and +! symmetry and of the provided Monkhorst-Pack-parameters. !----------------------------------------------------------------------- CONTAINS SUBROUTINE kptmop( - > iofile,iokpt,kpri,ktest, - > idsyst,idtype,nmop,ikzero,kzero, - > rltv,bltv,nreg,nfulst,nbound,idimens, - > xvec,fnorm,fdist,ncorn,nface,nedge,cpoint, - > nsym,ccr,rlsymr,talfa,mkpt,mface,mdir, - < nkpt,divis,vkxyz,nkstar,wghtkp) + > idsyst,idtype,nmop, + > rltv,bltv,nbound,idimens, + > xvec,fnorm,fdist,ncorn,nface,nedge,cpoint, + > nsym,ccr,rlsymr,talfa,mkpt,mface,mdir, + < nkpt,vkxyz,wghtkp) c -c Meaning of variables: -c INPUT: +c Meaning of variables: +c INPUT: c -c Symmetry of lattice: -c idsyst : crystal system identification in MDDFT programs -c idtype : lattice type identification in MDDFT programs -c bltv : cartesian coordinates of basis vectors for -c Bravais lattice: bltv(ix,jn), ix=1,3; jn=1,3 -c rltv : cartesian coordinates of basis vectors for -c reciprocal lattice: rltv(ix,jn), ix=1,3; jn=1,3 -c nsym : number of symmetry elements of points group -c ccr : rotation matrix for symmetry element -c in cartesian representation -c rlsymr : rotation matrix for symmetry element -c in reciprocal lattice basis representation -c talfa : translation vector associated with (non-symmorphic) -c symmetry elements in Bravais lattice representation +c Symmetry of lattice: +c idsyst : crystal system identification in MDDFT programs +c idtype : lattice type identification in MDDFT programs +c bltv : cartesian coordinates of basis vectors for +c Bravais lattice: bltv(ix,jn), ix=1,3; jn=1,3 +c rltv : cartesian coordinates of basis vectors for +c reciprocal lattice: rltv(ix,jn), ix=1,3; jn=1,3 +c nsym : number of symmetry elements of points group +c ccr : rotation matrix for symmetry element +c in cartesian representation +c rlsymr : rotation matrix for symmetry element +c in reciprocal lattice basis representation +c talfa : translation vector associated with (non-symmorphic) +c symmetry elements in Bravais lattice representation c -c representation of the irreducible part of the BZ: -c fnorm : normal vector of the planes bordering the irrBZ -c fdist : distance vector of the planes bordering the irrBZ -c ncorn : number of corners of the irrBZ -c nface : number of faces of the irrBZ -c nedge : number of edges of the irrBZ -c xvec : arbitrary vector lying in the irrBZ (FOR SURE!!) -c components are: +c representation of the irreducible part of the BZ: +c fnorm : normal vector of the planes bordering the irrBZ +c fdist : distance vector of the planes bordering the irrBZ +c ncorn : number of corners of the irrBZ +c nface : number of faces of the irrBZ +c nedge : number of edges of the irrBZ +c xvec : arbitrary vector lying in the irrBZ (FOR SURE!!) +c components are: c -c characterization of the Monkhorst-Pack k-point set: -c idimens : number of dimensions for k-point set (2 or 3) -c nreg : 1 kpoints in full BZ; 0 kpoints in irrBZ -c nfulst : 1 kpoints ordered in full stars -c (meaningful only for nreg =1; full BZ) -c nbound : 0 no primary points on BZ boundary; -c 1 with boundary points (not for BZ integration!!!) -c ikzero : 0 no shift of k-points; -c 1 shift of k-points for better use of sym in irrBZ -c kzero : shifting vector to bring one k-point to (0,0,0)(for even nmop) -c away from (0,0,0)(for odd nmop) -c nmop : integer number triple: nmop(i), i=1,3; nmop(i) -c determines number of k-points in direction of rltv(ix,i) +c characterization of the Monkhorst-Pack k-point set: +c idimens : number of dimensions for k-point set (2 or 3) +c nbound : 0 no primary points on BZ boundary; +c 1 with boundary points (not for BZ integration!!!) +c nmop : integer number triple: nmop(i), i=1,3; nmop(i) +c determines number of k-points in direction of rltv(ix,i) c -c OUTPUT: k-point set -c nkpt : number of k-points generated in set -c vkxyz : vector of kpoint generated; in cartesian representation -c wghtkp : weight associated with k-points for BZ integration -c divis : integer triple divis(i); i=1,4. -c Used to find more accurate representation of k-points -c vklmn(i,kpt)/divis(i) and weights as wght(kpt)/divis(4) -c nkstar : number of stars for k-points generated in full stars +c OUTPUT: k-point set +c nkpt : number of k-points generated in set +c vkxyz : vector of kpoint generated; in cartesian representation +c wghtkp : weight associated with k-points for BZ integration +c divis : integer triple divis(i); i=1,4. +c Used to find more accurate representation of k-points +c vklmn(i,kpt)/divis(i) and weights as wght(kpt)/divis(4) +c nkstar : number of stars for k-points generated in full stars c c----------------------------------------------------------------------- USE m_constants, ONLY : pimach @@ -85,14 +76,6 @@ C-----> PARAMETER STATEMENTS C INTEGER, INTENT (IN) :: mkpt,mface,mdir c -c ---> file number for read and write -c - INTEGER, INTENT (IN) :: iofile,iokpt -c -c ---> running mode parameter -c - INTEGER, INTENT (IN) :: kpri,ktest -C C-----> Symmetry information C INTEGER, INTENT (IN) :: nsym,idsyst,idtype @@ -111,15 +94,15 @@ C C C-----> BRILLOUINE ZONE INTEGRATION C - INTEGER, INTENT (IN) :: nreg,ikzero,nfulst,nbound,idimens + INTEGER, INTENT (IN) :: nbound,idimens INTEGER, INTENT (INOUT) :: nmop(3) - REAL, INTENT (INOUT) :: kzero(3) - INTEGER, INTENT (OUT):: nkpt,nkstar - REAL, INTENT (OUT):: vkxyz(3,mkpt),wghtkp(mkpt),divis(4) + INTEGER, INTENT (OUT):: nkpt + REAL, INTENT (OUT):: vkxyz(3,mkpt),wghtkp(mkpt) C -C ---> local variables +C ---> local variables c CHARACTER*80 blank + INTEGER nkstar,divis(4) INTEGER i,idim,i1,i2,i3,ii,ij,ik,is,isym,ifac, iik,iiik INTEGER ikc, i1red,nred,isumkpt,nleft,nirrbz INTEGER dirmin,dirmax,ndir1,ndir2,idir @@ -134,260 +117,165 @@ c INTEGER, ALLOCATABLE :: iside(:),iostar(:) REAL, ALLOCATABLE :: fract(:,:), vkrep(:,:) C -C ---> intrinsic functions +C ---> intrinsic functions c INTRINSIC real,abs C -C ---> save and data statements +C ---> save and data statements c SAVE one,zero,half,eps,eps1,iplus,iminus DATA zero/0.00/,one/1.00/,half/0.50/, - + eps/1.0e-8/,eps1/1.0e-6/,iplus/1/,iminus/-1/ + + eps/1.0e-8/,eps1/1.0e-6/,iplus/1/,iminus/-1/ c c----------------------------------------------------------------------- c ALLOCATE (fract(mkpt,3),vkrep(3,mkpt),ikpn(48,mkpt),irrkpn(mkpt)) ALLOCATE (nkrep(mkpt),iostar(mkpt),iside(mface)) - if (kpri .ge. 1) then -c write(iofile,'(/)') - write(iofile,'(3x,'' *<* kptmop *>* '')') - write(iofile,'(3x,'' generate k-vectors'')') - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~'')') - write(iofile,'(3x,'' by Monkhorst-Pack-method'')') - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~'')') - if (idimens .eq. 2) then - write(iofile,'(3x,'' in 2 dimensions'')') - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~'')') - write(iokpt,'(''#'',i4,20x,'' idimens: '', - > ''k-points for 2-dimensional lattice'')') idimens - else if (idimens .lt. 2 .or. idimens.gt.3) then - write(iofile,'(1x,i4,'' idimens wrong, choose 2 or 3'')') - > idimens - CALL juDFT_error("idimens",calledby="kptmop") - end if - if (nreg .eq. 1) then - if (nfulst .eq. 1) then - write(iofile,'(3x,'' in 1. Brillouin zone'')') - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~'')') - write(iofile,'(3x,'' full stars generated'')') - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~'')') - write(iokpt,'(''#'',2(i4,1x),14x,'' nreg,nfulst: '', - > ''k-points in totBZ, ordered in full stars'')') nreg,nfulst - else if (nfulst .eq. 0) then - write(iofile,'(3x,'' in parallelepiped'')') - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~'')') - write(iofile,'(3x,'' full stars not generated'')') - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~'')') - write(iokpt,'(''#'',2(i4,1x),14x,'' nreg,nfulst: '', - > ''k-points in par-epiped, not in full stars'')') nreg,nfulst - else - write(iofile,'(2(1x,i4),15x,'' nreg,nfulst: wrong choice '', - > ''of parameters; allowed combinations: (1,1); (1,0)'')' ) - > nreg,nfulst - CALL juDFT_error("nfulst",calledby="kptmop") - end if - else if (nreg .eq. 0) then - write(iofile,'(3x,'' in irred wedge of 1. Brillouin zone'')') - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')') - write(iokpt,'(''#'',i4,21x,''nreg: k-points in irrBZ'')') nreg - - else - write(iofile,'(3x,'' wrong choice of nreg: '', i4)') nreg - - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')') - CALL juDFT_error("nreg",calledby="kptmop") - end if - write(iokpt,'(''#'',3(i4,1x),10x,''nmop(1), nmop(2), nmop(3)'', - + '' as read in'')') (nmop(i), i=1,3) - write(iokpt,'(''#'',2(i4,1x),14x,'' ikzero, nbound'')') - + ikzero, nbound - write(iofile,'(/)') - if (nbound .eq. 1) then - write(iofile,'(3x,'' k-points on boundary included'')') - write(iofile,'(3x,'' irregular Monkhorst-Pack-ratios'')') - write(iofile,'(3x,'' cannot be used for BZ-integration'')') - else if (nbound.eq. 0) then - write(iofile,'(3x,'' no k-points on boundary of BZ'')') - write(iofile,'(3x,'' regular Monkhorst-Pack-ratios'')') - write(iofile,'(3x,'' can be used for BZ-integration'')') - else - write(iofile,'(3x,'' wrong choice of nbound: '', i4)') nbound - write(iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')') - CALL juDFT_error("nbound",calledby="kptmop") - end if -c - write(iofile,'(1x,i4,10x,''iofile'')') iofile - write(iofile,'(1x,i4,10x,''iokpt'')') iokpt - write(iofile,'(1x,i4,10x,''kpri'')') kpri - write(iofile,'(1x,i4,10x,''ktest'')') ktest - write(iofile,'(1x,3(f10.7,1x),10x,''xvec'')') (xvec(ii),ii=1,3) - write(iofile,'(1x,i4,10x,''ncorn'')') ncorn - write(iofile,'(1x,i4,10x,''nedge'')') nedge - write(iofile,'(1x,i4,10x,''nface'')') nface - do 5 ifac = 1,nface - write(iofile,'(1x,i4,1x,3(f10.7,1x),10x,''fnorm'')') - + ifac,(fnorm(ii,ifac),ii=1,3) - write(iofile,'(1x,i4,1x,f10.7,1x,10x,''fdist'')') - + ifac,fdist(ifac) - 5 continue -c - end if + ! -! ---> for 2 dimensions only the following Bravais lattices exist: -! TYPE EQUIVALENT 3-DIM idsyst/idtype -! square = p-tetragonal ( 1+2 axis ) 2/1 -! rectangular = p-orthorhomb ( 1+2 axis ) 3/1 -! centered rectangular = c-face-orthorhomb( 1+2 axis) 3/6 -! hexagonal = p-hexagonal ( 1+2 axis ) 4/1 -! oblique = p-monoclinic ( 1+2 axis ) 6/1 +! ---> for 2 dimensions only the following Bravais lattices exist: +! TYPE EQUIVALENT 3-DIM idsyst/idtype +! square = p-tetragonal ( 1+2 axis ) 2/1 +! rectangular = p-orthorhomb ( 1+2 axis ) 3/1 +! centered rectangular = c-face-orthorhomb( 1+2 axis) 3/6 +! hexagonal = p-hexagonal ( 1+2 axis ) 4/1 +! oblique = p-monoclinic ( 1+2 axis ) 6/1 ! IF (idimens .EQ. 2) THEN ! -! ---> identify the allowed symmetries -! and check the consistency of the Monkhorst-Pack-parameters -! - IF (idsyst.EQ.2 .OR. idsyst.EQ.4) THEN - IF (idtype.EQ.1) THEN - IF (nmop(3).EQ.1) nmop(3) = 0 - IF (nmop(1).NE.nmop(2) .OR. nmop(3).NE.0) THEN - nmop(2) = nmop(1) - nmop(3) = 0 - WRITE (iofile,'(1x,''WARNING!!!!!!!'',/, - + ''nmop-Parameters not in accordance with symmetry'',/, - + 2(1x,i4),/, - + '' we have set nmop(2) = nmop(1)'',/, - + '' and/or nmop(3) = 0'')') idsyst, idtype - WRITE (iofile,'(3(1x,i4),'' new val for nmop: '')') - + (nmop(i),i=1,3) - CALL juDFT_warn( - + "k point mesh not compatible with symmetry (1)", - + calledby='kptmop') - ELSE - WRITE (iofile,'('' values accepted unchanged'')') - WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')') - + (nmop(i),i=1,3) - ENDIF - ENDIF - ELSEIF (idsyst.EQ.3) THEN - IF (idtype.EQ.1 .OR. idtype.EQ.6) THEN - IF (nmop(3).EQ.1) nmop(3) = 0 - IF (nmop(3).NE.0) THEN - nmop(3) = 0 - WRITE (iofile,'(1x,''WARNING!!!!!!!'',/, - + ''nmop-Parameters not in accordance with symmetry'',/, - + 2(1x,i4),/, - + '' we have set nmop(3) = 0'')') idsyst, idtype - WRITE (iofile,'(3(1x,i4),'' new val for nmop: '')') - + (nmop(i),i=1,3) - CALL juDFT_warn( - + "k point mesh not compatible with symmetry (2)", - + calledby='kptmop') - ELSE - WRITE (iofile,'('' values accepted unchanged'')') - WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')') - + (nmop(i),i=1,3) - ENDIF - ENDIF - ELSEIF (idsyst.EQ.6) THEN - IF (idtype.EQ.1) THEN - IF (nmop(3).EQ.1) nmop(3) = 0 - IF (nmop(3).NE.0) THEN - nmop(3) = 0 - WRITE (iofile,'(1x,''WARNING!!!!!!!'',/, - + ''nmop-Parameters not in accordance with symmetry'',/, - + 2(1x,i4),/, - + '' we have set nmop(3) = 0'')') idsyst, idtype - WRITE (iofile,'(3(1x,i4),'' new val for nmop: '')') - + (nmop(i),i=1,3) - CALL juDFT_warn( - + "k point mesh not compatible with symmetry (3)", - + calledby='kptmop') - ELSE - WRITE (iofile,'('' values accepted unchanged'')') - WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')') - + (nmop(i),i=1,3) - ENDIF - ENDIF - ELSE +! ---> identify the allowed symmetries +! and check the consistency of the Monkhorst-Pack-parameters +! + IF (idsyst.EQ.2 .OR. idsyst.EQ.4) THEN + IF (idtype.EQ.1) THEN + If (Nmop(1)==1) Nmop(1)=0 + IF (nmop(1).NE.nmop(2) .OR. nmop(3).NE.0) THEN + nmop(2) = nmop(1) + nmop(3) = 0 + WRITE (6,'(1x,''WARNING!!!!!!!'',/, + +''nmop-Parameters not in accordance with symmetry'',/, + +2(1x,i4),/, + +'' we have set nmop(2) = nmop(1)'',/, + +'' and/or nmop(3) = 0'')') idsyst, idtype + WRITE (6,'(3(1x,i4),'' new val for nmop: '')') + + (nmop(i),i=1,3) + ELSE + WRITE (6,'('' values accepted unchanged'')') + WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')') + + (nmop(i),i=1,3) + ENDIF + ENDIF + ELSEIF (idsyst.EQ.3) THEN + IF (idtype.EQ.1 .OR. idtype.EQ.6) THEN + If (Nmop(3)==1) Nmop(3)=0 + IF (nmop(3).NE.0) THEN + nmop(3) = 0 + WRITE (6,'(1x,''WARNING!!!!!!!'',/, + +''nmop-Parameters not in accordance with symmetry'',/, + +2(1x,i4),/, + +'' we have set nmop(3) = 0'')') idsyst, idtype + WRITE (6,'(3(1x,i4),'' new val for nmop: '')') + + (nmop(i),i=1,3) + ELSE + WRITE (6,'('' values accepted unchanged'')') + WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')') + + (nmop(i),i=1,3) + ENDIF + ENDIF + ELSEIF (idsyst.EQ.6) THEN + IF (idtype.EQ.1) THEN + If (Nmop(3)==1) Nmop(3)=0 + IF (nmop(3).NE.0) THEN + nmop(3) = 0 + WRITE (6,'(1x,''WARNING!!!!!!!'',/, + +''nmop-Parameters not in accordance with symmetry'',/, + +2(1x,i4),/, + +'' we have set nmop(3) = 0'')') idsyst, idtype + WRITE (6,'(3(1x,i4),'' new val for nmop: '')') + + (nmop(i),i=1,3) + ELSE + WRITE (6,'('' values accepted unchanged'')') + WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')') + + (nmop(i),i=1,3) + ENDIF + ENDIF + ELSE ! -! ---> in all other cases: +! ---> in all other cases: ! - WRITE (iofile,'(3(1x,i4),20x,'' idimens,idsyst,idtype: '', - > ''wrong choice for 2-dimensional crystal structure'')') - > idimens,idsyst,idtype - CALL juDFT_error("2-dim crystal",calledby="kptmop") - ENDIF + WRITE (6,'(3(1x,i4),20x,'' idimens,idsyst,idtype: '', + >''wrong choice for 2-dimensional crystal structure'')') + > idimens,idsyst,idtype + CALL juDFT_error("2-dim crystal",calledby="kptmop") + ENDIF ! -! ---> check consistency of nmop-parameters with crystal symmetry +! ---> check consistency of nmop-parameters with crystal symmetry ! ELSEIF (idimens.EQ.3) THEN IF (idsyst.EQ.1 .OR. idsyst.EQ.5) THEN - IF (nmop(1).NE.nmop(2) .OR. nmop(1).NE.nmop(3) - + .OR. nmop(2).NE.nmop(3)) THEN - nmop(3) = nmop(1) - nmop(2) = nmop(1) - WRITE (iofile,'(1x,''WARNING!!!!!!!'',/, - + ''nmop-Parameters not in accordance with symmetry'',/, - + 2(1x,i4),/, - + '' we have set all nmop(i) = nmop(1)'')') idsyst, idtype - WRITE (iofile,'(3(1x,i4),'' new val for nmop(i): '')') - + (nmop(i),i=1,3) - CALL juDFT_warn( - + "k point mesh not compatible with symmetry (4)", - + calledby='kptmop') - ELSE - WRITE (iofile,'('' values accepted unchanged'')') - WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')') - + (nmop(i),i=1,3) - ENDIF + IF (nmop(1).NE.nmop(2) .OR. nmop(1).NE.nmop(3) + + .OR. nmop(2).NE.nmop(3)) THEN + nmop(3) = nmop(1) + nmop(2) = nmop(1) + WRITE (6,'(1x,''WARNING!!!!!!!'',/, + +''nmop-Parameters not in accordance with symmetry'',/, + +2(1x,i4),/, + +'' we have set all nmop(i) = nmop(1)'')') idsyst, idtype + WRITE (6,'(3(1x,i4),'' new val for nmop(i): '')') + + (nmop(i),i=1,3) + ELSE + WRITE (6,'('' values accepted unchanged'')') + WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')') + + (nmop(i),i=1,3) + ENDIF ELSEIF (idsyst.EQ.2 .OR. idsyst.eq.4) THEN - if((nmop(3).eq.nmop(2)).and.idsyst.eq.2)then - WRITE (iofile,'('' values accepted unchanged'')') - WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')') - + (nmop(i),i=1,3) - elseif(nmop(1).NE.nmop(2)) THEN - nmop(2) = nmop(1) - WRITE (iofile,'(1x,''WARNING!!!!!!!'',/, - + ''nmop-Parameters not in accordance with symmetry'',/, - + 2(1x,i4),/, - + '' we have set nmop(2) = nmop(1)'')') idsyst, idtype - WRITE (iofile,'(3(1x,i4),'' new val for nmop: '')') - + (nmop(i),i=1,3) - CALL juDFT_warn( - + "k point mesh not compatible with symmetry (5)", + if((nmop(3).eq.nmop(2)).and.idsyst.eq.2)then + WRITE (6,'('' values accepted unchanged'')') + WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')') + + (nmop(i),i=1,3) + elseif(nmop(1).NE.nmop(2)) THEN + nmop(2) = nmop(1) + WRITE (6,'(1x,''WARNING!!!!!!!'',/, + +''nmop-Parameters not in accordance with symmetry'',/, + +2(1x,i4),/, + +'' we have set nmop(2) = nmop(1)'')') idsyst, idtype + WRITE (6,'(3(1x,i4),'' new val for nmop: '')') + + (nmop(i),i=1,3) + CALL juDFT_warn( + + "k point mesh not compatible with symmetry (1)", + calledby='kptmop') - ELSE - WRITE (iofile,'('' values accepted unchanged'')') - WRITE (iofile,'(3(1x,i4),14x,''nmop(i),i=1,3'')') - + (nmop(i),i=1,3) - ENDIF + ELSE + WRITE (6,'('' values accepted unchanged'')') + WRITE (6,'(3(1x,i4),14x,''nmop(i),i=1,3'')') + + (nmop(i),i=1,3) + ENDIF ELSEIF (idsyst.LT.1 .OR. idsyst.GT.7) THEN - WRITE (iofile,'(1x,''wrong choice of symmetry'',/, - + 2(1x,i4))') idsyst, idtype - WRITE (iofile,'(''only values 1.le.idsyst.le.7 allowed'')') - CALL juDFT_error("wrong idsyst",calledby="kptmop") + WRITE (6,'(1x,''wrong choice of symmetry'',/, + +2(1x,i4))') idsyst, idtype + WRITE (6,'(''only values 1.le.idsyst.le.7 allowed'')') + CALL juDFT_error("wrong idsyst",calledby="kptmop") ELSE - WRITE (iofile,'('' values accepted unchanged'')') - WRITE (iofile,'(3(1x,i4),11x,''nmop(i),i=1,3'')') - + (nmop(i),i=1,3) + WRITE (6,'('' values accepted unchanged'')') + WRITE (6,'(3(1x,i4),11x,''nmop(i),i=1,3'')') + + (nmop(i),i=1,3) ENDIF ELSE - CALL juDFT_error("idimens =/= 2,3 ",calledby="kptmop") + CALL juDFT_error("idimens =/= 2,3 ",calledby="kptmop") ENDIF ! -! ---> start calculation -! ===================================================================== +! ---> start calculation +! ===================================================================== ! -! ---> set sign constants - isi(1) = 0 - isi(2) = iminus - isi(3) = iplus +! ---> set sign constants + isi(1) = 0 + isi(2) = iminus + isi(3) = iplus ! -! ---> calc orientation of boundary faces of irr wedge of BZ -! characterized by -! iside(i)= sign( (xvec,fnorm(i))-fdist(i) ) ;(i=1,nface ) +! ---> calc orientation of boundary faces of irr wedge of BZ +! characterized by +! iside(i)= sign( (xvec,fnorm(i))-fdist(i) ) ;(i=1,nface ) ! - WRITE (iofile,'(1x,''orientation of boundary faces'')') + WRITE (6,'(1x,''orientation of boundary faces'')') DO ifac = 1, nface orient = zero iside(ifac) = iplus @@ -396,100 +284,85 @@ c ENDDO orient = orient - fdist(ifac) IF (orient .LT. 0) iside(ifac) = iminus - WRITE (iofile,'(1x,2(i4,2x),f10.7,10x,''ifac,iside,orient'', - + '' for xvec'')') ifac,iside(ifac),orient + WRITE (6,'(1x,2(i4,2x),f10.7,10x,''ifac,iside,orient'', + +'' for xvec'')') ifac,iside(ifac),orient ENDDO invtpi = one / ( 2.0 * pimach() ) - WRITE (iofile,'(''Bravais lattice vectors'')' ) + WRITE (6,'(''Bravais lattice vectors'')' ) DO ii = 1, 3 - WRITE (iofile,'(43x,3(1x,f11.6))') (bltv(ii,ikc), ikc=1,3) + WRITE (6,'(43x,3(1x,f11.6))') (bltv(ii,ikc), ikc=1,3) ENDDO - WRITE (iofile,'(''reciprocal lattice vectors'')' ) + WRITE (6,'(''reciprocal lattice vectors'')' ) DO ii = 1, 3 - WRITE (iofile,'(43x,3(1x,f11.6))' ) (rltv(ii,ikc), ikc=1,3) + WRITE (6,'(43x,3(1x,f11.6))' ) (rltv(ii,ikc), ikc=1,3) ENDDO ! -! ---> nmop(i) are Monkhorst-Pack parameters; they determine the -! number of k-points in i-direction -! if basis vector lengths are not related by symmetry, -! we can use independent fractions for each direction +! ---> nmop(i) are Monkhorst-Pack parameters; they determine the +! number of k-points in i-direction +! if basis vector lengths are not related by symmetry, +! we can use independent fractions for each direction ! - WRITE (iofile,'(3(1x,i4),10x,'' Monkhorst-Pack-parameters'')') - + (nmop(i1),i1=1,3) + WRITE (6,'(3(1x,i4),10x,'' Monkhorst-Pack-parameters'')') + + (nmop(i1),i1=1,3) DO idim = 1, idimens IF (nmop(idim).GT.0) THEN ainvnmop(idim) = one/ real(nmop(idim)) ELSE - WRITE (iofile,'('' nmop('',i4,'') ='',i4, - + '' not allowed'')') idim, nmop(idim) - CALL juDFT_error("nmop wrong",calledby="kptmop") + WRITE (6,'('' nmop('',i4,'') ='',i4, + +'' not allowed'')') idim, nmop(idim) + CALL juDFT_error("nmop wrong",calledby="kptmop") ENDIF ENDDO -! -! --->nreg determines region of BZ in which k-points are generated -! - IF ( nreg .EQ. 1) THEN - IF (nfulst .EQ. 1) THEN - WRITE (iofile,'(2(1x,i4),15x,'' nreg,nfulst: '', - > ''k-points in totBZ, ordered in full stars'')') nreg,nfulst - ELSEIF (nfulst .EQ. 0) then - WRITE (iofile,'(2(1x,i4),15x,'' nreg,nfulst: '', - > ''k-points in par-epiped, not in full stars'')') nreg,nfulst - ENDIF - ELSEIF ( nreg .EQ. 0) THEN - WRITE (iofile,'(1x,i4,10x,''nreg; k-points in'', - + '' irreducible wedge of BZ'')' ) nreg - ENDIF - WRITE (iofile,'(1x,''Monkhorst-Pack-fractions'')' ) + WRITE (6,'(1x,''Monkhorst-Pack-fractions'')' ) ! -! ---> nbound=1: k-points are generated on boundary of BZ -! include fract(1) = -1/2 -! and fract(2*nmop+1) = 1/2 for surface points of BZ +! ---> nbound=1: k-points are generated on boundary of BZ +! include fract(1) = -1/2 +! and fract(2*nmop+1) = 1/2 for surface points of BZ ! IF ( nbound .EQ. 1) THEN - WRITE (iofile,'(1x,i4,10x,''nbound; k-points on boundary'', - + '' of BZ included'')' ) nbound + WRITE (6,'(1x,i4,10x,''nbound; k-points on boundary'', + +'' of BZ included'')' ) nbound ! -! ---> irregular Monkhorst--Pack--fractions -! fract(r) = r / (2*nmop) +! ---> irregular Monkhorst--Pack--fractions +! fract(r) = r / (2*nmop) ! - DO idim = 1,idimens - denom = half*ainvnmop(idim) - divis(idim) = one / denom + DO idim = 1,idimens + denom = half*ainvnmop(idim) + divis(idim) = one / denom - DO kpn = -nmop(idim),nmop(idim) - fract(kpn+nmop(idim)+1,idim) = denom * real (kpn) - WRITE (iofile,'(10x,f10.7)' ) fract(kpn+nmop(idim)+1,idim) - ENDDO - nfract(idim) = 2*nmop(idim) + 1 - ENDDO - IF (idimens .eq. 2) THEN - nfract(3) = 1 - fract(1,3) = 0 - divis(3) = one - END IF -! -! ---> nbound=0: k-points are NOT generated on boundary of BZ -! This is the regular Monkhorst-Pack-method + DO kpn = -nmop(idim),nmop(idim) + fract(kpn+nmop(idim)+1,idim) = denom * real (kpn) + WRITE (6,'(10x,f10.7)' ) fract(kpn+nmop(idim)+1,idim) + ENDDO + nfract(idim) = 2*nmop(idim) + 1 + ENDDO + IF (idimens .eq. 2) THEN + nfract(3) = 1 + fract(1,3) = 0 + divis(3) = one + END IF +! +! ---> nbound=0: k-points are NOT generated on boundary of BZ +! This is the regular Monkhorst-Pack-method ! ELSEIF ( nbound .eq. 0) then - WRITE (iofile,'(1x,i4,10x,''nbound; no k-points '', - + '' on boundary of BZ'')' ) nbound + WRITE (6,'(1x,i4,10x,''nbound; no k-points '', + +'' on boundary of BZ'')' ) nbound ! -! ---> regular Monkhorst--Pack--fractions -! fract(r) =(2*r-nmop-1) / (2*nmop) +! ---> regular Monkhorst--Pack--fractions +! fract(r) =(2*r-nmop-1) / (2*nmop) ! DO idim = 1,idimens denom = half*ainvnmop(idim) divis(idim) = one / denom - WRITE(iofile,'(5x,i4,5x,''idim'')' ) idim + WRITE(6,'(5x,i4,5x,''idim'')' ) idim DO kpn = 1,nmop(idim) - fract(kpn,idim) = denom * real (2*kpn -nmop(idim)-1) - write(iofile,'(10x,f10.7)' ) fract(kpn,idim) + fract(kpn,idim) = denom * real (2*kpn -nmop(idim)-1) + write(6,'(10x,f10.7)' ) fract(kpn,idim) ENDDO nfract(idim) = nmop(idim) ENDDO @@ -501,76 +374,14 @@ c ENDIF ELSE - WRITE (iofile,'(3x,'' wrong choice of nbound:'', i4)') nbound - WRITE (iofile,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')') - CALL juDFT_error("nbound",calledby="kptmop") + WRITE (6,'(3x,'' wrong choice of nbound:'', i4)') nbound + WRITE (6,'(3x,'' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'')') + CALL juDFT_error("nbound",calledby="kptmop") ENDIF -! -! ---> set kzero -! - kzero(1) = zero - kzero(2) = zero - kzero(3) = zero - - IF (ikzero.NE.0) THEN -! -! ---> set kzero .ne. zero -! for nmop even and non-orthogonal rec lattice vectors -! shift all k-points to bring one to (0,0,0) -! for nmop odd and non-orthogonal rec lattice vectors -! shift k-points to avoid (0,0,0) -! implemented for hexagon (idsyst =4, idtype arbitrary) -! and ort(2-dim)(idsyst =3, idtype =1) -! and fcc (idsyst =1, idtype =3) -! and p-monocli (idsyst =6, idtype =1) -! - denom = half*ainvnmop(1) - - IF( idsyst .EQ. 4) THEN ! hexagonal Bravais lattice - - DO ii = 1, 3 - kzero(ii) = denom * (rltv(ii,1)+rltv(ii,2)) - ENDDO - ELSEIF ( idsyst .EQ. 3 .and. idtype .EQ. 1) THEN ! orthorhombic (2D) - - IF (nmop(1).EQ.1) THEN - denom = half*ainvnmop(3) - DO ii = 1, 3 - kzero(ii) = denom * (rltv(ii,2)+rltv(ii,3)) - ENDDO - ELSE - WRITE (iofile,'(''ikzero.ne.0 is NOT applied '', - > ''for present choice of parameters'')') - WRITE (iofile,'(5(1x,i4), - > ''; idsyst,idtype,nmop(1),nmop(2),nmop(3)'')') - > idsyst,idtype,nmop(1),nmop(2),nmop(3) - ENDIF - - ELSEIF ( idsyst .eq. 1 .and. idtype .eq. 3) then ! face centered cubic - - DO ii = 1, 3 - kzero(ii) = denom * (rltv(ii,1)+rltv(ii,2)+rltv(ii,3)) - ENDDO - - ELSEIF ( idsyst .eq. 6 .and. idtype .eq. 1) then ! p-monoclinic - - DO ii = 1, 3 - kzero(ii) = half*ainvnmop(1) * rltv(ii,1) - + + half*ainvnmop(2) * rltv(ii,2) - ENDDO - - ELSE - WRITE (iofile,'(''ikzero.ne.0 is NOT applied '', - > ''for this system'')') - WRITE (iofile,'(5(1x,i4), - > ''; idsyst,idtype,nmop(1),nmop(2),nmop(3)'')') - > idsyst,idtype,nmop(1),nmop(2),nmop(3) - ENDIF - ENDIF ! ! -! ---> initialize k-points = zero and weights = 1.0 +! ---> initialize k-points = zero and weights = 1.0 ! DO kpn = 1,mkpt vkxyz(1,kpn) = zero @@ -579,167 +390,104 @@ c wghtkp(kpn) = one ENDDO ! -! ---> generate equidistant k-vectors in cartesian coordinates -! with off-set kzero +! ---> generate equidistant k-vectors in cartesian coordinates ! nkpt = 0 DO i3 = 1,nfract(3) DO i2 = 1,nfract(2) DO i1 = 1,nfract(1) - nkpt = nkpt + 1 - IF (nkpt>mkpt) CALL juDFT_error("nkpt > mkpt",calledby - + ="kptmop") - vkxyz(1,nkpt) = kzero(1) + rltv(1,1)*fract(i1,1) - + + rltv(1,2)*fract(i2,2) - + + rltv(1,3)*fract(i3,3) - vkxyz(2,nkpt) = kzero(2) + rltv(2,1)*fract(i1,1) - + + rltv(2,2)*fract(i2,2) - + + rltv(2,3)*fract(i3,3) - vkxyz(3,nkpt) = kzero(3) + rltv(3,1)*fract(i1,1) - + + rltv(3,2)*fract(i2,2) - + + rltv(3,3)*fract(i3,3) + nkpt = nkpt + 1 + IF (nkpt>mkpt) CALL juDFT_error("nkpt > mkpt",calledby + + ="kptmop") + vkxyz(1,nkpt) = rltv(1,1)*fract(i1,1) + + + rltv(1,2)*fract(i2,2) + + + rltv(1,3)*fract(i3,3) + vkxyz(2,nkpt) = rltv(2,1)*fract(i1,1) + + + rltv(2,2)*fract(i2,2) + + + rltv(2,3)*fract(i3,3) + vkxyz(3,nkpt) = rltv(3,1)*fract(i1,1) + + + rltv(3,2)*fract(i2,2) + + + rltv(3,3)*fract(i3,3) ENDDO ENDDO ENDDO ! -! ---> calculate weights of k-points and print out k-points -! wghtkp = 1/nkpt -! ( = 1/(nmop(1)*nmop(2)*nmop(3)) for reg Monk-Pack-method) +! ---> calculate weights of k-points and print out k-points +! wghtkp = 1/nkpt +! ( = 1/(nmop(1)*nmop(2)*nmop(3)) for reg Monk-Pack-method) ! - divis(4) = real(nkpt) - aivnkpt = one/real(nkpt) +! divis(4) = real(nkpt) +! aivnkpt = one/real(nkpt) - DO kpn= 1,nkpt - wghtkp(kpn) = wghtkp(kpn)*aivnkpt - ENDDO +! DO kpn= 1,nkpt +! wghtkp(kpn) = wghtkp(kpn)*aivnkpt +! ENDDO - IF (kpri .ge. 1) THEN - WRITE (iofile,'(''generated k-points vkxyz'')') - DO kpn= 1,nkpt - WRITE (iofile,'(1x,i4,1x,4(f13.10,1x),10x,''vkxyz, wghtkp'')') - + kpn,(vkxyz(ii,kpn),ii=1,3),wghtkp(kpn) - ENDDO - ENDIF - -! -! ---> for nreg=1 and nfulst=0 -! (k-points in total zone of arbitrary shape) -! k-point generation finished -! (k-points are actually generated in parallel-epiped -! (r1,r2,r3) which is an equivalent unit cell) -! The generated set does not contain all points of the -! symmetry related stars. ! - IF (nreg.EQ.1 .AND. nfulst.EQ.0) GOTO 1000 +! ==================================================================== ! -! ==================================================================== +! ---> order generated k-points in stars by applying symmetry: +! - determine number of different stars nkstar .le. nkpt +! - determine order of star iostar(kpn) .le. nsym +! - assign pointer ikpn(i,ik); i=1,iostar(ik); ik=1,nkstar +! - determine representative vector in irrBZ for each star: +! vkrep(ix,ik); ix=1,3; ik=1,nkstar ! -! ---> order generated k-points in stars by applying symmetry: -! - determine number of different stars nkstar .le. nkpt -! - determine order of star iostar(kpn) .le. nsym -! - assign pointer ikpn(i,ik); i=1,iostar(ik); ik=1,nkstar -! - determine representative vector in irrBZ for each star: -! vkrep(ix,ik); ix=1,3; ik=1,nkstar + CALL ordstar( + > 6,0,0, + > fnorm,fdist,nface,iside, + > nsym,ccr,rltv,mkpt,mface,mdir, + = nkpt,vkxyz, + < nkstar,iostar,ikpn,vkrep,nkrep) ! - CALL ordstar( - > iokpt,kpri,ktest, - > fnorm,fdist,nface,iside, - > nsym,ccr,rltv,mkpt,mface,mdir, - = nkpt,vkxyz, - < nkstar,iostar,ikpn,vkrep,nkrep) ! -! ---> for nreg=0 (k-points in irr part of BZ) +! (a) calculate weights for k-points in irrBZ +! - wghtkp(ik)=iostar(ik)/nkpt_old ; ik=1,nkstar ! -! (a) calculate weights for k-points in irrBZ -! - wghtkp(ik)=iostar(ik)/nkpt_old ; ik=1,nkstar -! - DO ik = 1, nkstar - wghtkp(ik) = wghtkp(ik)*iostar(ik) - ENDDO + DO ik = 1, nkstar + wghtkp(ik) = wghtkp(ik)*iostar(ik) + ENDDO ! -! (b) final preparation of k-points for transfer to file -! - assign nkpt= nkstar -! - assign vkxyz(ix,ik) = vkrep(ix,ik); ix=1,3; ik=1,nkstar +! (b) final preparation of k-points for transfer to file +! - assign nkpt= nkstar +! - assign vkxyz(ix,ik) = vkrep(ix,ik); ix=1,3; ik=1,nkstar ! - IF (nreg.EQ.0) THEN - DO i1 = 1,3 + DO i1 = 1,3 DO ik = 1,nkstar - vkxyz(i1,ik) = vkrep(i1,ik) + vkxyz(i1,ik) = vkrep(i1,ik) ENDDO - ENDDO - nkpt = nkstar - GOTO 1000 - ENDIF -c -c ---> for nreg.eq.0: k-point generation finished -c -c ================================================================ -c -c - IF (nreg.EQ.1 .AND. nfulst.EQ.1) THEN -c -c ---> generate full stars for all representative k-points -c - for nreg=1 and nfulst=1: -c - determine order of full star ifstar(kpn).le.nsym -c - assign nkpt= sum {ifstar(ik)} (ik=1,nkstar) -c - assign vkxyz(ix,kpn) = vkstar(ix,ikpn(is,ik)); -c ix=1,3; kpn=1,nkpt; ik=1,nstar; is=1,ifstar(ik) -c - calculate wghtkp(kpn)=wghtkp_old(ik)/(ifstar(ik) -c kpn=1,nkpt; ik=1,nstar -c -c - CALL fulstar( - > iofile,iokpt,kpri,ktest, - > ccr,nsym, - > vkrep,nkstar,mkpt,mface,mdir, - = nkpt,vkxyz,wghtkp) - - divis(4) = divis(4) * nsym - - ENDIF - - 1000 CONTINUE -! -! --> check for corner points, include them into k-point set: + ENDDO + nkpt = nkstar +! +! --> check for corner points, include them into k-point set: ! IF (nbound.EQ.1) THEN - n = 1 - nc2d = 1 ! determine 2D corner points - cp2d(:,nc2d) = cpoint(:,n) - corn: DO n = 2, ncorn - DO i = 1, n-1 + n = 1 + nc2d = 1 ! determine 2D corner points + cp2d(:,nc2d) = cpoint(:,n) + corn: DO n = 2, ncorn + DO i = 1, n-1 IF ((abs(cpoint(1,n)-cpoint(1,i)).LT.0.0001).AND. - + (abs(cpoint(2,n)-cpoint(2,i)).LT.0.0001)) CYCLE corn - ENDDO - nc2d = nc2d + 1 - cp2d(:,nc2d) = cpoint(:,n) - ENDDO corn - WRITE (iofile,'(''2D corner points in internal units'')') - corn2d: DO n = 1, nc2d - WRITE (iofile,'(i3,3x,2(f10.7,1x))') n,cp2d(1,n),cp2d(2,n) - DO i = 1, nkpt - IF ((abs(cp2d(1,n)-vkxyz(1,i)).LT.0.0001).AND. - + (abs(cp2d(2,n)-vkxyz(2,i)).LT.0.0001)) CYCLE corn2d - ENDDO - nkpt = nkpt + 1 - vkxyz(:,nkpt) = cp2d(:,n) - ENDDO corn2d - ENDIF -! -! ---> print out k-points and weights -! - IF (ktest .GE. 1) THEN - WRITE (iofile,'(''generated k-points vkxyz'')') - WRITE (iofile,'(1x,i4,20x, - + ''nkpt,number of generated k-points'')') nkpt - - DO kpn = 1, nkpt - WRITE (iofile,'(3(f10.7,1x),f12.10,1x,i4,3x, - + ''vkxyz, wghtkp'')') (vkxyz(ii,kpn),ii=1,3),wghtkp(kpn), kpn - ENDDO - + + (abs(cpoint(2,n)-cpoint(2,i)).LT.0.0001)) CYCLE corn + ENDDO + nc2d = nc2d + 1 + cp2d(:,nc2d) = cpoint(:,n) + ENDDO corn + WRITE (6,'(''2D corner points in internal units'')') + corn2d: DO n = 1, nc2d + WRITE (6,'(i3,3x,2(f10.7,1x))') n,cp2d(1,n),cp2d(2,n) + DO i = 1, nkpt + IF ((abs(cp2d(1,n)-vkxyz(1,i)).LT.0.0001).AND. + + (abs(cp2d(2,n)-vkxyz(2,i)).LT.0.0001)) CYCLE corn2d + ENDDO + nkpt = nkpt + 1 + vkxyz(:,nkpt) = cp2d(:,n) + ENDDO corn2d ENDIF +! +! ---> print out k-points and weights +! DEALLOCATE (fract,vkrep,ikpn,irrkpn,nkrep,iostar,iside) RETURN diff --git a/main/cdngen.F90 b/main/cdngen.F90 index df7d2597809c599d5ec8f9249fd194518063eb58..5d78359250b9b4e328d700cdda32c97e1c04c96d 100644 --- a/main/cdngen.F90 +++ b/main/cdngen.F90 @@ -95,7 +95,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,& LOGICAL :: l_error, perform_MetaGGA CALL regCharges%init(input,atoms) - CALL dos%init(input,atoms,dimension,kpts,vacuum) + CALL dos%init(dimension%neigd,input,atoms,kpts,vacuum) CALL moments%init(input,atoms) CALL mcd%init1(banddos,dimension,input,atoms,kpts) CALL slab%init(banddos,dimension,atoms,cell,input,kpts) @@ -197,10 +197,10 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,& END IF IF (noco%l_soc) CALL orbMagMoms(input,atoms,noco,moments%clmom) - + END IF END IF ! mpi%irank == 0 - + perform_MetaGGA = ALLOCATED(EnergyDen%mt) & .AND. (xcpot%exc_is_MetaGGA() .or. xcpot%vx_is_MetaGGA()) if(perform_MetaGGA) then diff --git a/main/fleur.F90 b/main/fleur.F90 index 09372588a4cc970c884cf200aa856977abaa3223..430a29374d58250166bbf74f2e734d6e7e4fc200 100644 --- a/main/fleur.F90 +++ b/main/fleur.F90 @@ -102,7 +102,7 @@ CONTAINS ! local scalars INTEGER :: eig_id,archiveType - INTEGER :: n,iter,iterHF + INTEGER :: n,iter,iterHF,num_threads LOGICAL :: l_cont,l_qfix,l_real REAL :: fix #ifdef CPP_MPI @@ -128,11 +128,11 @@ CONTAINS IF(mpi%isize.NE.1) CALL juDFT_error('No Wannier+MPI at the moment',calledby = 'fleur') CALL wann_optional(input,kpts,atoms,sym,cell,oneD,noco,wann) END IF - + iter = 0 iterHF = 0 l_cont = (iter < input%itmax) - + IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('scfLoop') ! Initialize and load inDen density (start) @@ -160,7 +160,7 @@ CONTAINS ! Open/allocate eigenvector storage (start) l_real=sym%invs.AND..NOT.noco%l_noco eig_id=open_eig(mpi%mpi_comm,DIMENSION%nbasfcn,DIMENSION%neigd,kpts%nkpt,input%jspins,& - noco%l_noco,.NOT.INPUT%eig66(1),l_real,noco%l_soc,INPUT%eig66(1),mpi%n_size) + noco%l_noco,.true.,l_real,noco%l_soc,.false.,mpi%n_size) #ifdef CPP_CHASE CALL init_chase(mpi,dimension,input,atoms,kpts,noco,sym%invs.AND..NOT.noco%l_noco) @@ -196,7 +196,7 @@ CONTAINS CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,inDen) #endif - + !HF !$ num_threads = omp_get_max_threads() !$ call omp_set_num_threads(1) @@ -215,9 +215,9 @@ CONTAINS IF(input%l_rdmft) THEN CALL open_hybrid_io1(DIMENSION,sym%invs) END IF - IF(.not.input%eig66(1))THEN + !IF(.not.input%eig66(1))THEN CALL reset_eig(eig_id,noco%l_soc) ! This has to be placed after the calc_hybrid call but before eigen - END IF + !END IF !$ call omp_set_num_threads(num_threads) !#endif @@ -251,10 +251,10 @@ CONTAINS CALL timestart("Updating energy parameters") CALL enpara%update(mpi%mpi_comm,atoms,vacuum,input,vToT) CALL timestop("Updating energy parameters") - IF(.not.input%eig66(1))THEN + !IF(.not.input%eig66(1))THEN CALL eigen(mpi,stars,sphhar,atoms,xcpot,sym,kpts,DIMENSION,vacuum,input,& cell,enpara,banddos,noco,oneD,hybrid,iter,eig_id,results,inDen,vTemp,vx) - ENDIF + !ENDIF vTot%mmpMat = vTemp%mmpMat !!$ eig_idList(pc) = eig_id CALL timestop("eigen") @@ -277,7 +277,7 @@ CONTAINS #endif ! WRITE(6,fmt='(A)') 'Starting 2nd variation ...' - IF (noco%l_soc.AND..NOT.noco%l_noco.AND..NOT.INPUT%eig66(1)) & + IF (noco%l_soc.AND..NOT.noco%l_noco) & CALL eigenso(eig_id,mpi,DIMENSION,stars,vacuum,atoms,sphhar,& sym,cell,noco,input,kpts, oneD,vTot,enpara,results) CALL timestop("gen. of hamil. and diag. (total)") @@ -344,7 +344,7 @@ CONTAINS CYCLE forcetheoloop ENDIF - + !+Wannier functions IF ((input%l_wann).AND.(.NOT.wann%l_bs_comf)) THEN CALL wannier(DIMENSION,mpi,input,kpts,sym,atoms,stars,vacuum,sphhar,oneD,& @@ -392,7 +392,7 @@ CONTAINS #endif CALL timestop("generation of new charge density (total)") - + !!$ !----> output potential and potential difference !!$ IF (disp) THEN !!$ reap = .FALSE. @@ -404,7 +404,7 @@ CONTAINS !!$ !!$ CALL potdis(stars,vacuum,atoms,sphhar, input,cell,sym) !!$ END IF - + ! total energy CALL timestart('determination of total energy') CALL totale(mpi,atoms,sphhar,stars,vacuum,DIMENSION,sym,input,noco,cell,oneD,& @@ -416,21 +416,21 @@ CONTAINS CALL forcetheo%postprocess() - CALL enpara%mix(mpi,atoms,vacuum,input,vTot%mt(:,0,:,:),vtot%vacz) + CALL enpara%mix(mpi%mpi_comm,atoms,vacuum,input,vTot%mt(:,0,:,:),vtot%vacz) field2 = field ! mix input and output densities CALL mix_charge(field2,DIMENSION,mpi,(iter==input%itmax.OR.judft_was_argument("-mix_io")),& stars,atoms,sphhar,vacuum,input,& sym,cell,noco,oneD,archiveType,xcpot,iter,inDen,outDen,results) - + IF(mpi%irank == 0) THEN WRITE (6,FMT=8130) iter 8130 FORMAT (/,5x,'******* it=',i3,' is completed********',/,/) WRITE(*,*) "Iteration:",iter," Distance:",results%last_distance CALL timestop("Iteration") END IF ! mpi%irank.EQ.0 - + #ifdef CPP_MPI CALL MPI_BCAST(results%last_distance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr) CALL MPI_BARRIER(mpi%mpi_comm,ierr) @@ -452,7 +452,7 @@ CONTAINS ELSE l_cont = l_cont.AND.(iter < input%itmax) ! MetaGGAs need a at least 2 iterations - l_cont = l_cont.AND.((input%mindistance<=results%last_distance).OR.input%l_f & + l_cont = l_cont.AND.((input%mindistance<=results%last_distance).OR.input%l_f & .OR. (xcpot%exc_is_MetaGGA() .and. iter == 1)) CALL check_time_for_next_iteration(iter,l_cont) END IF @@ -464,7 +464,7 @@ CONTAINS END IF END DO scfloop ! DO WHILE (l_cont) - + CALL add_usage_data("Iterations",iter) IF (mpi%irank.EQ.0) CALL closeXMLElement('scfLoop') @@ -472,7 +472,7 @@ CONTAINS CALL close_eig(eig_id) CALL juDFT_end("all done",mpi%irank) - + CONTAINS SUBROUTINE priv_geo_end(mpi) TYPE(t_mpi),INTENT(IN)::mpi @@ -499,6 +499,6 @@ CONTAINS ENDIF CALL juDFT_end(" GEO new inp.xml created ! ",mpi%irank) END SUBROUTINE priv_geo_end - + END SUBROUTINE fleur_execute END MODULE m_fleur diff --git a/main/fleur_init.F90 b/main/fleur_init.F90 index a991e431ccb3145bdeec6c61d0a2398b2e6179f5..6aa39e6a7e7180813fc3da0e42aa7dd71d5a864f 100644 --- a/main/fleur_init.F90 +++ b/main/fleur_init.F90 @@ -36,7 +36,7 @@ CONTAINS USE m_types_xcpot_inbuild USE m_mpi_bc_xcpot - + #ifdef CPP_MPI USE m_mpi_bc_all, ONLY : mpi_bc_all #ifndef CPP_OLDINTEL @@ -134,17 +134,17 @@ CONTAINS xcpot,forcetheo_data,kpts,enparaXML,mpi%mpi_comm) - CALL timestart("postprocessInput") + CALL timestart("postprocessInput") CALL postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,kpts,& oneD,hybrid,cell,banddos,sliceplot,xcpot,forcetheo,forcetheo_data,& noco,DIMENSION,enpara,enparaxml,sphhar,l_kpts) - CALL timestop("postprocessInput") + CALL timestop("postprocessInput") IF (mpi%irank.EQ.0) THEN CALL w_inpXML(& atoms,vacuum,input,stars,sliceplot,forcetheo,banddos,& cell,sym,xcpot,noco,oneD,hybrid,kpts,enpara,& - .TRUE.,[.TRUE.,.TRUE.,.TRUE.,.TRUE.]) + .TRUE.,[.TRUE.,.TRUE.,.TRUE.,.TRUE.]) END IF @@ -166,7 +166,7 @@ CONTAINS END IF !Finalize the MPI setup - CALL setupMPI(kpts%nkpt,mpi) + CALL setupMPI(kpts%nkpt,dimension%neigd,mpi) !Collect some usage info CALL add_usage_data("A-Types",atoms%ntype) diff --git a/types/types_enpara.F90 b/types/types_enpara.F90 index 20bfa6070d5b0950bfa38b2953ba63153e0455ab..c44fcb12c73c865a3e5b1b25989c8befcd1e89cd 100644 --- a/types/types_enpara.F90 +++ b/types/types_enpara.F90 @@ -46,9 +46,9 @@ CONTAINS INTEGER,INTENT(IN) :: jspins LOGICAL,INTENT(IN) :: film TYPE(t_enparaxml),OPTIONAL :: enparaxml - + LOGICAL :: l_enpara - + ALLOCATE(this%el0(0:atoms%lmaxd,atoms%ntype,jspins),this%el1(0:atoms%lmaxd,atoms%ntype,jspins)) ALLOCATE(this%ello0(atoms%nlod,atoms%ntype,jspins),this%ello1(atoms%nlod,atoms%ntype,jspins)) this%el0=-1E99 @@ -83,8 +83,8 @@ CONTAINS this%evac0=enparaxml%evac0 END IF endif - - + + END SUBROUTINE init_enpara !> This subroutine adjusts the energy parameters to the potential. In particular, it @@ -111,6 +111,11 @@ CONTAINS REAL :: vbar,vz0,rj INTEGER :: n,jsp,l,ilo,j,ivac,irank CHARACTER(LEN=20) :: attributes(5) + REAL :: e_lo(0:3,atoms%ntype)!Store info on branches to do IO after OMP + REAL :: e_up(0:3,atoms%ntype) + REAL :: elo_lo(atoms%nlod,atoms%ntype) + REAL :: elo_up(atoms%nlod,atoms%ntype) + #ifdef CPP_MPI INCLUDE 'mpif.h' INTEGER :: ierr @@ -123,20 +128,20 @@ CONTAINS l_done = .FALSE.;lo_done=.FALSE. DO jsp = 1,input%jspins !$OMP PARALLEL DO DEFAULT(none) & - !$OMP SHARED(atoms,enpara,jsp,l_done,irank,v,lo_done) & + !$OMP SHARED(atoms,enpara,jsp,l_done,v,lo_done,e_lo,e_up,elo_lo,elo_up) & !$OMP PRIVATE(n,l,ilo) !! First calculate energy paramter from quantum numbers if these are given... !! l_done stores the index of those energy parameter updated DO n = 1, atoms%ntype DO l = 0,3 - IF( enpara%qn_el(l,n,jsp).ne.0)THEN + IF( enpara%qn_el(l,n,jsp).ne.0)THEN l_done(l,n,jsp) = .TRUE. - enpara%el0(l,n,jsp)=find_enpara(.FALSE.,l,n,jsp,enpara%qn_el(l,n,jsp),atoms,irank,v%mt(:,0,n,jsp)) + enpara%el0(l,n,jsp)=find_enpara(.FALSE.,l,n,jsp,enpara%qn_el(l,n,jsp),atoms,v%mt(:,0,n,jsp),e_lo(l,n),e_up(l,n)) IF( l .EQ. 3 ) THEN enpara%el0(4:,n,jsp) = enpara%el0(3,n,jsp) l_done(4:,n,jsp) = .TRUE. END IF - ELSE + ELSE l_done(l,n,jsp) = .FALSE. END IF ENDDO ! l @@ -145,32 +150,32 @@ CONTAINS l = atoms%llo(ilo,n) IF( enpara%qn_ello(ilo,n,jsp).NE.0) THEN lo_done(ilo,n,jsp) = .TRUE. - enpara%ello0(ilo,n,jsp)=find_enpara(.TRUE.,l,n,jsp,enpara%qn_ello(ilo,n,jsp),atoms,irank,v%mt(:,0,n,jsp)) + enpara%ello0(ilo,n,jsp)=find_enpara(.TRUE.,l,n,jsp,enpara%qn_ello(ilo,n,jsp),atoms,v%mt(:,0,n,jsp),e_lo(l,n),e_up(l,n)) ELSE lo_done(ilo,n,jsp) = .FALSE. ENDIF ENDDO ENDDO ! n !$OMP END PARALLEL DO - IF (mpi%irank==0) THEN - WRITE(6,*) - WRITE(6,*) "Updated energy parameters for spin:",jsp - !Same loop for IO - DO n = 1, atoms%ntype - DO l = 0,3 - IF( l_done(l,n,jsp)) CALL priv_write(.FALSE.,l,n,jsp,enpara%qn_el(l,n,jsp),e_lo(l,n),e_up(l,n),enpara%el0(l,n,jsp)) - ENDDO ! l - ! Now for the lo's - DO ilo = 1, atoms%nlo(n) - l = atoms%llo(ilo,n) - IF( lo_done(ilo,n,jsp)) CALL priv_write(.TRUE.,l,n,jsp,enpara%qn_ello(ilo,n,jsp),elo_lo(ilo,n),elo_up(ilo,n),enpara%ello0(ilo,n,jsp)) - END DO - END DO + IF (irank==0) THEN + WRITE(6,*) + WRITE(6,*) "Updated energy parameters for spin:",jsp + !Same loop for IO + DO n = 1, atoms%ntype + DO l = 0,3 + IF( l_done(l,n,jsp)) CALL priv_write(.FALSE.,l,n,jsp,enpara%qn_el(l,n,jsp),e_lo(l,n),e_up(l,n),enpara%el0(l,n,jsp)) + ENDDO ! l + ! Now for the lo's + DO ilo = 1, atoms%nlo(n) + l = atoms%llo(ilo,n) + IF( lo_done(ilo,n,jsp)) CALL priv_write(.TRUE.,l,n,jsp,enpara%qn_ello(ilo,n,jsp),elo_lo(ilo,n),elo_up(ilo,n),enpara%ello0(ilo,n,jsp)) + END DO + END DO ENDIF !! Now check for floating energy parameters (not for those with l_done=T) IF (enpara%floating) THEN - types_loop: DO n = 1,atoms%ntype + types_loop: DO n = 1,atoms%ntype ! !---> determine energy parameters if lepr=1. the reference energy !---> is the value of the l=0 potential at approximately rmt/4. @@ -213,7 +218,7 @@ CONTAINS IF (enpara%floating) THEN vz0 = v%vacz(1,ivac,jsp) IF (irank.EQ.0) THEN - WRITE ( 6,'('' spin'',i2,'', vacuum '',i3,'' ='',f12.6)') jsp,ivac,vz0 + WRITE ( 6,'('' spin'',i2,'', vacuum '',i3,'' ='',f12.6)') jsp,ivac,vz0 ENDIF ENDIF enpara%evac(ivac,jsp) = enpara%evac0(ivac,jsp) + vz0 @@ -239,7 +244,7 @@ CONTAINS ! enpara%ready=(ALL(enpara%el0>-1E99).AND.ALL(enpara%ello0>-1E99)) enpara%epara_min=MIN(MINVAL(enpara%el0),MINVAL(enpara%ello0)) - + IF (irank == 0) CALL closeXMLElement('energyParameters') END SUBROUTINE update @@ -272,9 +277,9 @@ CONTAINS skip_t = 0 DO n = 1,atoms%ntype READ (40,FMT=8040,END=200) (enpara%el0(l,n,jsp),l=0,3),& - (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp) + (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp) WRITE (6,FMT=8140) n,(enpara%el0(l,n,jsp),l=0,3),& - (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp) + (enpara%lchange(l,n,jsp),l=0,3),enpara%skiplo(n,jsp) ! !---> energy parameters for the local orbitals ! @@ -286,7 +291,7 @@ CONTAINS WRITE (6,FMT=8138) (enpara%llochg(lo,n,jsp),lo=1,atoms%nlo(n)) ELSEIF (enpara%skiplo(n,jsp).GT.0) THEN WRITE (6,*) "for atom",n," no LO's were specified" - WRITE (6,*) 'but skiplo was set to',enpara%skiplo + WRITE (6,*) 'but skiplo was set to',enpara%skiplo CALL juDFT_error("No LO's but skiplo",calledby ="enpara",& hint="If no LO's are set skiplo must be 0 in enpara") END IF @@ -310,9 +315,9 @@ CONTAINS READ (40,FMT=8050,END=200) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp) WRITE (6,FMT=8150) enpara%evac0(1,jsp),enpara%lchg_v(1,jsp),enpara%evac0(2,jsp) ENDIF - ! IF (atoms%nlod.GE.1) THEN + ! IF (atoms%nlod.GE.1) THEN ! WRITE (6,FMT=8090) jsp,skip_t - ! WRITE (6,FMT=8091) + ! WRITE (6,FMT=8091) ! END IF END DO @@ -501,8 +506,8 @@ CONTAINS enpara%el0(:,n,jsp)=enpara%el0(:,n,jsp)-vbar ENDIF END DO - - + + IF (input%film) THEN WRITE(6,*) 'Vacuum:' DO n=1,vacuum%nvac @@ -531,7 +536,7 @@ CONTAINS CALL MPI_BCAST(enpara%el0,SIZE(enpara%el0),MPI_DOUBLE_PRECISION,0,mpi_comm,ierr) CALL MPI_BCAST(enpara%ello0,SIZE(enpara%ello0),MPI_DOUBLE_PRECISION,0,mpi_comm,ierr) CALL MPI_BCAST(enpara%evac,SIZE(enpara%evac),MPI_DOUBLE_PRECISION,0,mpi_comm,ierr) -#endif +#endif RETURN 777 FORMAT('Old:',f8.5,' new:',f8.5,' diff:',f8.5) @@ -557,8 +562,7 @@ CONTAINS IF (input%film) enpara%evac1(:vacuum%nvac,ispin)=regCharges%pvac(:vacuum%nvac,ispin)/regCharges%svac(:vacuum%nvac,ispin) END DO END SUBROUTINE calcOutParams - - SUBROUTINE priv_write(lo,l,n,jsp,nqn,e_lo,e_up,e) +SUBROUTINE priv_write(lo,l,n,jsp,nqn,e_lo,e_up,e) !subroutine to write energy parameters to output USE m_xmlOutput IMPLICIT NONE @@ -603,5 +607,4 @@ CONTAINS END IF END SUBROUTINE priv_write - END MODULE m_types_enpara diff --git a/types/types_xcpot_inbuild.F90 b/types/types_xcpot_inbuild.F90 index 517bc65b0d3b24f1150c695ad169374813b704a7..24c63495f9c4eb66d4debf78fc6176da4eb6fd55 100644 --- a/types/types_xcpot_inbuild.F90 +++ b/types/types_xcpot_inbuild.F90 @@ -7,177 +7,26 @@ MODULE m_types_xcpot_inbuild !This module contains the xcpot-type used for the in-build xc-implementations USE m_types_xcpot_data USE m_types_xcpot + USE m_types_xcpot_inbuild_nofunction USE m_judft IMPLICIT NONE PRIVATE - REAL, PARAMETER, PRIVATE :: hrtr_half = 0.5 CHARACTER(len=4),PARAMETER:: xc_names(20)=[& 'l91 ','x-a ','wign','mjw ','hl ','bh ','vwn ','pz ', & 'pw91','pbe ','rpbe','Rpbe','wc ','PBEs', & 'pbe0','hse ','vhse','lhse','exx ','hf '] - LOGICAL,PARAMETER:: priv_LDA(20)=[& - .FALSE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,& - .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,& - .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.] - - - LOGICAL,PARAMETER:: priv_gga(20)=[& - .TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,& - .TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,& - .TRUE.,.TRUE.,.TRUE.,.TRUE.,.FALSE.,.TRUE.] - - LOGICAL,PARAMETER:: priv_hybrid(20)=[& - .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,& - .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,& - .TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.] - - REAL, PARAMETER :: amix_pbe0 = 0.25 - REAL, PARAMETER :: amix_hse = 0.25 - REAL, PARAMETER :: amix_hf = 1.00 - - TYPE, EXTENDS(t_xcpot):: t_xcpot_inbuild -#ifdef CPP_MPI - INTEGER :: icorr=0 !not private to allow bcasting it around -#else - INTEGER,PRIVATE :: icorr=0 -#endif - - TYPE(t_xcpot_data) :: DATA - - LOGICAL,ALLOCATABLE :: lda_atom(:) - + TYPE, EXTENDS(t_xcpot_inbuild_nf):: t_xcpot_inbuild CONTAINS !overloading t_xcpot: - PROCEDURE :: vx_is_LDA => xcpot_vx_is_LDA - PROCEDURE :: vx_is_GGA => xcpot_vx_is_GGA - - PROCEDURE :: vc_is_LDA => xcpot_vc_is_LDA - PROCEDURE :: vc_is_GGA => xcpot_vc_is_GGA - - PROCEDURE :: exc_is_LDA => xcpot_exc_is_LDA - PROCEDURE :: exc_is_gga => xcpot_exc_is_gga - PROCEDURE :: is_hybrid => xcpot_is_hybrid - - PROCEDURE :: get_exchange_weight => xcpot_get_exchange_weight PROCEDURE :: get_vxc => xcpot_get_vxc PROCEDURE :: get_exc => xcpot_get_exc - !not overloaded - PROCEDURE :: get_name => xcpot_get_name - PROCEDURE :: is_name => xcpot_is_name - PROCEDURE :: init => xcpot_init - END TYPE t_xcpot_inbuild + END TYPE t_xcpot_inbuild PUBLIC t_xcpot_inbuild -CONTAINS - CHARACTER(len=4) FUNCTION xcpot_get_name(xcpot) - USE m_judft - IMPLICIT NONE - CLASS(t_xcpot_inbuild),INTENT(IN) :: xcpot - IF (xcpot%icorr==0) CALL judft_error("xc-potential not initialized",calledby="types_xcpot.F90") - xcpot_get_name=xc_names(xcpot%icorr) - END FUNCTION xcpot_get_name - - SUBROUTINE xcpot_init(xcpot,namex,relcor,ntype) - USE m_judft - IMPLICIT NONE - CLASS(t_xcpot_inbuild),INTENT(INOUT) :: xcpot - CHARACTER(len=*),INTENT(IN) :: namex - LOGICAL,INTENT(IN) :: relcor - INTEGER,INTENT(IN) :: ntype - INTEGER:: n - !Determine icorr from name - - ALLOCATE(xcpot%lda_atom(ntype)) - xcpot%lda_atom=.FALSE. - xcpot%icorr=0 - DO n=1,SIZE(xc_names) - IF (TRIM(ADJUSTL(namex))==TRIM(xc_names(n))) THEN - xcpot%icorr=n - ENDIF - ENDDO - if (xcpot%icorr==0) CALL judft_error("Unkown xc-potential:"//namex,calledby="types_xcpot.F90") - xcpot%data%krla=MERGE(1,0,relcor) + CONTAINS - !Code from exchpbe to speed up determination of constants - IF (xcpot%is_name("rpbe")) THEN - xcpot%data%uk=1.2450 - ELSE - xcpot%data%uk=0.8040 - ENDIF - IF (xcpot%is_name("PBEs")) THEN ! pbe_sol - xcpot%data%um=0.123456790123456 - ELSE - xcpot%data%um=0.2195149727645171 - ENDIF - xcpot%data%is_hse=xcpot%is_name("hse").OR.xcpot%is_name("lhse").OR.xcpot%is_name("vhse") - xcpot%data%is_rpbe=xcpot%is_name("Rpbe") !Rpbe - xcpot%data%is_wc=xcpot%is_name("wc") - xcpot%data%is_pbes=xcpot%is_name("PBEs") - xcpot%data%is_pbe0=xcpot%is_name("pbe0") - xcpot%data%is_mjw=xcpot%is_name("mjw") - xcpot%data%is_bh=xcpot%is_name("bh") - xcpot%DATA%exchange_weight=xcpot%get_exchange_weight() - - END SUBROUTINE xcpot_init - - !! LDA - logical function xcpot_exc_is_lda(xcpot) - implicit none - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - xcpot_exc_is_lda= xcpot%vxc_is_lda() - end function xcpot_exc_is_lda - - logical function xcpot_vx_is_lda(xcpot) - implicit none - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - xcpot_vx_is_lda=(.not. xcpot%vxc_is_gga()) .and. (.not. xcpot%is_hybrid()) - end function xcpot_vx_is_lda - - logical function xcpot_vc_is_lda(xcpot) - implicit none - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - xcpot_vc_is_lda=(.not. xcpot%vxc_is_gga()) .and. (.not. xcpot%is_hybrid()) - end function xcpot_vc_is_lda - - logical function xcpot_vx_is_gga(xcpot) - implicit none - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - xcpot_vx_is_gga=priv_gga(xcpot%icorr) - end function xcpot_vx_is_gga - - logical function xcpot_vc_is_gga(xcpot) - implicit none - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - xcpot_vc_is_gga=priv_gga(xcpot%icorr) - end function xcpot_vc_is_gga - - LOGICAL FUNCTION xcpot_exc_is_gga(xcpot) - IMPLICIT NONE - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - xcpot_exc_is_gga = xcpot%vxc_is_gga() - END FUNCTION xcpot_exc_is_gga - - LOGICAL FUNCTION xcpot_is_hybrid(xcpot) - IMPLICIT NONE - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - xcpot_is_hybrid=priv_hybrid(xcpot%icorr) - END FUNCTION xcpot_is_hybrid - - FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex) - USE m_judft - IMPLICIT NONE - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - - REAL:: a_ex - - a_ex=-1 - IF (xcpot%is_name("pbe0")) a_ex=amix_pbe0 - IF (xcpot%is_name("hf")) a_ex=amix_hf - IF (xcpot%is_name("hse")) a_ex=amix_hse - IF (xcpot%is_name("vhse")) a_ex=amix_hse - END FUNCTION xcpot_get_exchange_weight - - SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad, kinED_KS) + + SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad) ! USE m_xcxal, ONLY : vxcxal USE m_xcwgn, ONLY : vxcwgn @@ -203,7 +52,6 @@ CONTAINS !c REAL, INTENT (OUT) :: vx (:,:) REAL, INTENT (OUT) :: vxc(:,:) - REAL, INTENT(IN),OPTIONAL:: kinED_KS(:,:) ! optional arguments for GGA TYPE(t_gradients),INTENT(INOUT),OPTIONAL::grad @@ -280,7 +128,7 @@ CONTAINS END SUBROUTINE xcpot_get_vxc !*********************************************************************** - SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinED_KS, mt_call) + SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinEnergyDen_KS, mt_call) !*********************************************************************** USE m_xcxal, ONLY : excxal USE m_xcwgn, ONLY : excwgn @@ -299,7 +147,7 @@ CONTAINS REAL, INTENT (OUT) :: exc(:) TYPE(t_gradients),OPTIONAL,INTENT(IN) ::grad LOGICAL, OPTIONAL, INTENT(IN) :: mt_call - REAL, INTENT(IN), OPTIONAL :: kinED_KS(:,:) + REAL, INTENT(IN), OPTIONAL :: kinEnergyDen_KS(:,:) !c !c ---> local scalars @@ -351,10 +199,4 @@ CONTAINS END SUBROUTINE xcpot_get_exc - LOGICAL FUNCTION xcpot_is_name(xcpot,name) - CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot - CHARACTER(len=*),INTENT(IN) :: name - xcpot_is_name=(TRIM(xc_names(xcpot%icorr))==TRIM((name))) - END FUNCTION xcpot_is_name - -END MODULE m_types_xcpot_inbuild + END MODULE m_types_xcpot_inbuild diff --git a/types/types_xcpot_libxc.F90 b/types/types_xcpot_libxc.F90 index f2bddca0507fe18a03d0918367a8ef9d0141d033..91b3a7dac3d37fc5aefbf194e4927b88bd66cbb9 100644 --- a/types/types_xcpot_libxc.F90 +++ b/types/types_xcpot_libxc.F90 @@ -24,14 +24,16 @@ MODULE m_types_xcpot_libxc TYPE(xc_f03_func_t) :: vxc_func_x, vxc_func_c TYPE(xc_f03_func_t) :: exc_func_x, exc_func_c #endif - INTEGER :: func_vxc_id_c, func_vxc_id_x !> functionals to be used for potential & density convergence - INTEGER :: func_exc_id_c, func_exc_id_x !> functionals to be used in exc- & totale-calculations INTEGER :: jspins + + CONTAINS + !PROCEDURE :: vxc_is_LDA => xcpot_vxc_is_LDA + !PROCEDURE :: vxc_is_gga => xcpot_vxc_is_gga + PROCEDURE :: vx_is_LDA => xcpot_vx_is_LDA PROCEDURE :: vx_is_GGA => xcpot_vx_is_GGA PROCEDURE :: vx_is_MetaGGA => xcpot_vx_is_MetaGGA - PROCEDURE :: vx_is_XC => xcpot_vx_is_XC ! exchange and correlation in one functional PROCEDURE :: vc_is_LDA => xcpot_vc_is_LDA PROCEDURE :: vc_is_GGA => xcpot_vc_is_GGA @@ -39,7 +41,6 @@ MODULE m_types_xcpot_libxc PROCEDURE :: exc_is_LDA => xcpot_exc_is_LDA PROCEDURE :: exc_is_gga => xcpot_exc_is_gga PROCEDURE :: exc_is_MetaGGA => xcpot_exc_is_MetaGGA - PROCEDURE :: ex_is_XC => xcpot_ex_is_XC ! exchange and correlation in one functional PROCEDURE :: is_hybrid => xcpot_is_hybrid PROCEDURE :: get_exchange_weight => xcpot_get_exchange_weight @@ -52,78 +53,81 @@ MODULE m_types_xcpot_libxc PUBLIC t_xcpot_libxc CONTAINS - SUBROUTINE xcpot_init(xcpot,jspins,vxc_id_x,vxc_id_c,exc_id_x,exc_id_c) + + + + SUBROUTINE xcpot_init(xcpot,func_vxc_id_x,func_vxc_id_c,func_exc_id_x,func_exc_id_c,jspins) USE m_judft IMPLICIT NONE - CLASS(t_xcpot_libxc),INTENT(OUT) :: xcpot - INTEGER, INTENT(IN) :: jspins - INTEGER, INTENT(IN) :: vxc_id_x, vxc_id_c ! potential functional - INTEGER, INTENT(IN) :: exc_id_x, exc_id_c ! energy functionals + CLASS(t_xcpot_libxc),INTENT(INOUT) :: xcpot + INTEGER, INTENT(IN) :: jspins,func_vxc_id_x,func_vxc_id_c,func_exc_id_x,func_exc_id_c LOGICAL :: same_functionals ! are vxc and exc equal INTEGER :: errors(4) #ifdef CPP_LIBXC INTEGER :: err xcpot%jspins = jspins + xcpot%l_libxc=.TRUE. + xcpot%func_vxc_id_x=func_vxc_id_x + xcpot%func_exc_id_x=func_exc_id_x + xcpot%func_vxc_id_c=func_vxc_id_c + xcpot%func_exc_id_c=func_exc_id_c - xcpot%func_vxc_id_x = vxc_id_x - xcpot%func_vxc_id_c = vxc_id_c - - xcpot%func_exc_id_x = exc_id_x - xcpot%func_exc_id_c = exc_id_c - - if(xcpot%func_vxc_id_x == 0 .or. xcpot%func_exc_id_x == 0 ) then + IF(xcpot%func_vxc_id_x == 0 .OR. xcpot%func_exc_id_x == 0 ) THEN CALL judft_error("LibXC exchange- and correlation-function indicies need to be set"& - ,hint='Try this: ' // ACHAR(10) //& - '' // ACHAR(10) //& - ' ' // ACHAR(10) //& - ' ') - endif - + ,hint='Try this: ' // ACHAR(10) //& + '' // ACHAR(10) //& + ' ' // ACHAR(10) //& + ' ') + ENDIF + IF (jspins==1) THEN ! potential functionals CALL xc_f03_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_UNPOLARIZED, err=errors(1)) IF (xcpot%func_vxc_id_c>0) CALL xc_f03_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, & - XC_UNPOLARIZED, err=errors(2)) - + XC_UNPOLARIZED, err=errors(2)) ! energy functionals CALL xc_f03_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_UNPOLARIZED, err=errors(3)) IF (xcpot%func_exc_id_c>0) CALL xc_f03_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, & - XC_UNPOLARIZED, err=errors(4)) + XC_UNPOLARIZED, err=errors(4)) ELSE ! potential functionals CALL xc_f03_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_POLARIZED, err=errors(1)) IF (xcpot%func_vxc_id_c>0) CALL xc_f03_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, & - XC_POLARIZED, err=errors(2)) - + XC_POLARIZED, err=errors(2)) + !energy functionals CALL xc_f03_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_POLARIZED, err=errors(3)) IF (xcpot%func_exc_id_c>0) CALL xc_f03_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, & - XC_POLARIZED, err=errors(4)) + XC_POLARIZED, err=errors(4)) END IF - - IF(errors(1) /= 0) call juDFT_error("Exchange potential functional not in LibXC") - IF(errors(2) /= 0) call juDFT_error("Correlation potential functional not in LibXC") - IF(errors(3) /= 0) call juDFT_error("Exchange energy functional not in LibXC") - IF(errors(4) /= 0) call juDFT_error("Correlation energy functional not in LibXC") - + + IF(errors(1) /= 0) CALL juDFT_error("Exchange potential functional not in LibXC") + IF(errors(2) /= 0) CALL juDFT_error("Correlation potential functional not in LibXC") + IF(errors(3) /= 0) CALL juDFT_error("Exchange energy functional not in LibXC") + IF(errors(4) /= 0) CALL juDFT_error("Correlation energy functional not in LibXC") + !check if any potental is a MetaGGA - IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_c))) THEN - call juDFT_error("vxc_x: MetaGGA is not implemented for correclation potentials") + IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_x))) THEN + CALL juDFT_error("vxc_x: MetaGGA is not implemented for potentials") + ELSEIF(xcpot%func_vxc_id_c>0) THEN + IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_c))) THEN + CALL juDFT_error("vxc_x: MetaGGA is not implemented for potentials") + ENDIF ENDIF - + CALL write_xc_info(xcpot%vxc_func_x) - + IF (xcpot%func_vxc_id_c>0) THEN CALL write_xc_info(xcpot%vxc_func_c) ELSE WRITE(*,*) "No Correlation functional" END IF - + same_functionals = (xcpot%func_vxc_id_x == xcpot%func_exc_id_x) & - .and. (xcpot%func_vxc_id_c == xcpot%func_exc_id_c) - IF(.not. same_functionals) THEN + .AND. (xcpot%func_vxc_id_c == xcpot%func_exc_id_c) + IF(.NOT. same_functionals) THEN CALL write_xc_info(xcpot%exc_func_x) IF (xcpot%func_exc_id_c>0) THEN CALL write_xc_info(xcpot%exc_func_c) @@ -131,18 +135,16 @@ CONTAINS WRITE(*,*) "No Correlation functional for TotalE" ENDIF ELSE - write (*,*) "Using same functional for VXC and EXC" + WRITE (*,*) "Using same functional for VXC and EXC" END IF - - xcpot%kinED%set = .False. #else CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", & - hint="Please recompile FLEUR with libxc support") + hint="Please recompile FLEUR with libxc support") #endif - END SUBROUTINE xcpot_init - - ! LDA - LOGICAL FUNCTION xcpot_vx_is_LDA(xcpot) + END SUBROUTINE xcpot_init + + ! LDA + LOGICAL FUNCTION xcpot_vx_is_LDA(xcpot) IMPLICIT NONE CLASS(t_xcpot_libxc),INTENT(IN):: xcpot #ifdef CPP_LIBXC @@ -174,8 +176,8 @@ CONTAINS #ifdef CPP_LIBXC TYPE(xc_f03_func_info_t) :: xc_info - xc_info = xc_f03_func_get_info(xcpot%exc_func_x) - xcpot_exc_is_LDA = (XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)) + xc_info = xc_f03_func_get_info(xcpot%vxc_func_x) + xcpot_exc_is_LDA = XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info) #else xcpot_exc_is_LDA = .false. #endif @@ -201,7 +203,7 @@ CONTAINS #ifdef CPP_LIBXC TYPE(xc_f03_func_info_t) :: xc_info - xc_info = xc_f03_func_get_info(xcpot%vxc_func_x) + xc_info = xc_f03_func_get_info(xcpot%vxc_func_c) xcpot_vx_is_gga = ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info)) #else xcpot_vx_is_gga=.false. @@ -214,38 +216,12 @@ CONTAINS #ifdef CPP_LIBXC TYPE(xc_f03_func_info_t) :: xc_info - xc_info = xc_f03_func_get_info(xcpot%vxc_func_x) + xc_info = xc_f03_func_get_info(xcpot%vxc_func_c) xcpot_vx_is_MetaGGA = ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info)) #else xcpot_vx_is_MetaGGA=.false. #endif END FUNCTION xcpot_vx_is_MetaGGA - - LOGICAL FUNCTION xcpot_vx_is_XC(xcpot) - IMPLICIT NONE - CLASS(t_xcpot_libxc),INTENT(IN):: xcpot -#ifdef CPP_LIBXC - TYPE(xc_f03_func_info_t) :: xc_info - - xc_info = xc_f03_func_get_info(xcpot%vxc_func_x) - xcpot_vx_is_XC = xc_f03_func_info_get_kind(xc_info) == XC_EXCHANGE_CORRELATION -#else - xcpot_vx_is_XC=.false. -#endif - END FUNCTION xcpot_vx_is_XC - - LOGICAL FUNCTION xcpot_ex_is_XC(xcpot) - IMPLICIT NONE - CLASS(t_xcpot_libxc),INTENT(IN):: xcpot -#ifdef CPP_LIBXC - TYPE(xc_f03_func_info_t) :: xc_info - - xc_info = xc_f03_func_get_info(xcpot%exc_func_x) - xcpot_ex_is_XC = xc_f03_func_info_get_kind(xc_info) == XC_EXCHANGE_CORRELATION -#else - xcpot_ex_is_XC=.false. -#endif - END FUNCTION xcpot_ex_is_XC LOGICAL FUNCTION xcpot_exc_is_gga(xcpot) IMPLICIT NONE @@ -298,82 +274,49 @@ CONTAINS END FUNCTION xcpot_get_exchange_weight !*********************************************************************** - SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad, kinED_KS) - USE, INTRINSIC :: IEEE_ARITHMETIC + SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad) + !*********************************************************************** IMPLICIT NONE CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot INTEGER, INTENT (IN) :: jspins REAL,INTENT (IN) :: rh(:,:) !Dimensions here REAL, INTENT (OUT) :: vx (:,:) !points,spin REAL, INTENT (OUT ) :: vxc(:,:) ! - REAL, INTENT(IN),OPTIONAL:: kinED_KS(:,:) ! optional arguments for GGA TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad #ifdef CPP_LIBXC - REAL,ALLOCATABLE :: vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:), & - tmp_vsig(:,:), tmp_vlapl(:,:), tmp_vtau(:,:), & - kinED_libxc(:,:) - integer :: idx + REAL,ALLOCATABLE::vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:) !libxc uses the spin as a first index, hence we have to transpose.... ALLOCATE(vxc_tmp(SIZE(vxc,2),SIZE(vxc,1)));vxc_tmp=0.0 ALLOCATE(vx_tmp(SIZE(vx,2),SIZE(vx,1)));vx_tmp=0.0 - IF (xcpot%vc_is_GGA()) THEN + IF (xcpot%needs_grad()) THEN IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives") ALLOCATE(vsigma,mold=grad%vsigma) !where(abs(grad%sigma)<1E-9) grad%sigma=1E-9 CALL xc_f03_gga_vxc(xcpot%vxc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,vsigma) - - ! check if vx has both exchange and correlation - IF (.not. xcpot%vx_is_XC()) THEN + IF (xcpot%func_vxc_id_c>0) THEN CALL xc_f03_gga_vxc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vxc_tmp,grad%vsigma) grad%vsigma=grad%vsigma+vsigma vxc_tmp=vxc_tmp+vx_tmp + ELSE + vxc_tmp=vx_tmp ENDIF ELSE !LDA potentials - if(xcpot%vx_is_MetaGGA() .and. present(kinED_KS)) then - kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace) - - allocate(tmp_vsig, mold=vx_tmp) - allocate(tmp_vlapl, mold=vx_tmp) - allocate(tmp_vtau, mold=vx_tmp) - - call xc_f03_mgga_vxc(xcpot%vxc_func_x, size(rh,1), transpose(rh), & - grad%sigma, transpose(grad%laplace), kinED_libxc,& - vx_tmp, tmp_vsig, tmp_vlapl, tmp_vtau) - - idx = find_first_normal(vx_tmp)+1 - vx_tmp(:,:idx) = 0.0 - - CALL xc_f03_lda_vxc(initial_lda_func(jspins), idx, TRANSPOSE(rh(:idx,:)), vx_tmp(:,:idx)) - else - CALL xc_f03_lda_vxc(initial_lda_func(jspins), SIZE(rh,1), TRANSPOSE(rh), vx_tmp) - endif - IF (.not. xcpot%vx_is_XC()) THEN + CALL xc_f03_lda_vxc(xcpot%vxc_func_x, SIZE(rh,1), TRANSPOSE(rh), vx_tmp) + IF (xcpot%func_vxc_id_c>0) THEN CALL xc_f03_lda_vxc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh), vxc_tmp) vxc_tmp=vxc_tmp+vx_tmp ENDIF ENDIF vx=TRANSPOSE(vx_tmp) vxc=TRANSPOSE(vxc_tmp) + #endif END SUBROUTINE xcpot_get_vxc - FUNCTION find_first_normal(vec) result(idx) - USE, INTRINSIC :: IEEE_ARITHMETIC - implicit none - real, intent(in) :: vec(:,:) - integer :: idx - - idx = size(vec, dim=2) - do while(all(.not. ieee_is_nan(vec(:,idx)))) - idx = idx - 1 - if(idx == 0) exit - enddo - END FUNCTION find_first_normal - - SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinED_KS, mt_call) + SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinEnergyDen_KS, mt_call) use m_constants IMPLICIT NONE CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot @@ -388,7 +331,7 @@ CONTAINS ! tau = sum[phi_i(r)^dag nabla phi_i(r)] ! see eq (2) in https://doi.org/10.1063/1.1565316 ! (-0.5 is applied below) - REAL, INTENT(IN), OPTIONAL :: kinED_KS(:,:) + REAL, INTENT(IN), OPTIONAL :: kinEnergyDen_KS(:,:) #ifdef CPP_LIBXC TYPE(xc_f03_func_info_t) :: xc_info @@ -399,30 +342,26 @@ CONTAINS ! tau = 0.5 * sum[|grad phi_i(r)|²] ! see eq (3) in https://doi.org/10.1063/1.1565316 - REAL, ALLOCATABLE :: kinED_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:) + REAL, ALLOCATABLE :: kinEnergyDen_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:) is_mt = merge(mt_call, .False., present(mt_call)) IF (xcpot%exc_is_gga()) THEN IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a GGA potential without providing derivatives") CALL xc_f03_gga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc) - - !in case ex is both exchange and correlation - IF (.not. xcpot%ex_is_XC()) THEN + IF (xcpot%func_exc_id_c>0) THEN CALL xc_f03_gga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc) exc=exc+excc END IF ELSEIF(xcpot%exc_is_LDA()) THEN !LDA potentials CALL xc_f03_lda_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc) - - !in case ex is both exchange and correlation - IF (.not. xcpot%ex_is_XC()) THEN + IF (xcpot%func_exc_id_c>0) THEN CALL xc_f03_lda_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc) exc=exc+excc END IF ELSEIF(xcpot%exc_is_MetaGGA()) THEN - IF(PRESENT(kinED_KS)) THEN + IF(PRESENT(kinEnergyDen_KS)) THEN ! apply correction in eq (4) in https://doi.org/10.1063/1.1565316 - kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace) + kinEnergyDen_libXC = transpose(kinEnergyDen_KS + 0.25 * grad%laplace) !only cut core of muffin tin cut_idx = MERGE(NINT(size(rh,1) * cut_ratio), 0, is_mt) @@ -433,7 +372,7 @@ CONTAINS TRANSPOSE(rh(cut_idx+1:,:)), & grad%sigma(:,cut_idx+1:), & transpose(grad%laplace(cut_idx+1:,:)), & - kinED_libXC(:,cut_idx+1:), & + kinEnergyDen_libXC(:,cut_idx+1:), & exc(cut_idx+1:)) call xc_f03_gga_exc(xcpot%vxc_func_x, SIZE(rh(:cut_idx,:),1), & @@ -446,7 +385,7 @@ CONTAINS TRANSPOSE(rh(cut_idx+1:,:)), & grad%sigma(:,cut_idx+1:), & transpose(grad%laplace(cut_idx+1:,:)), & - kinED_libXC(:,cut_idx+1:), & + kinEnergyDen_libXC(:,cut_idx+1:), & excc(cut_idx+1:)) call xc_f03_gga_exc(xcpot%vxc_func_c, SIZE(rh(:cut_idx,:),1), & @@ -482,10 +421,6 @@ CONTAINS ALLOCATE(grad%laplace(ngrid,jspins)) ALLOCATE(grad%vsigma(MERGE(1,3,jspins==1),ngrid)) - grad%gr = 0.0 - grad%laplace = 0.0 - grad%sigma = 0.0 - grad%vsigma = 0.0 END SUBROUTINE xcpot_alloc_gradients #ifdef CPP_LIBXC @@ -550,23 +485,6 @@ CONTAINS family = xc_f03_func_info_get_family(xc_f03_func_get_info(xc_func)) END FUNCTION xc_get_family - function initial_lda_func(jspins) result(lda) - use, intrinsic :: iso_c_binding - implicit none - INTEGER, intent(in) :: jspins - TYPE(xc_f03_func_t) :: lda - integer :: ierr - - if(jspins == 1) then - CALL xc_f03_func_init(lda, 1, XC_UNPOLARIZED, err=ierr) - else - CALL xc_f03_func_init(lda, 1, XC_POLARIZED, err=ierr) - ENDIF - - if(ierr /= 0) call juDFT_error("can't find lda_x for init") - - end function initial_lda_func - #endif END MODULE m_types_xcpot_libxc diff --git a/vgen/vis_xc.F90 b/vgen/vis_xc.F90 index 0218dc4ea32f23730427cc7506c99d16a0e5e867..edf0b179d89cf3946b62b1fbcb24f73f52b590cc 100644 --- a/vgen/vis_xc.F90 +++ b/vgen/vis_xc.F90 @@ -62,13 +62,15 @@ CONTAINS ALLOCATE(v_xc,mold=rho) ALLOCATE(v_x,mold=rho) - +#ifdef CPP_LIBXC if(perform_MetaGGA .and. xcpot%kinED%set) then CALL xcpot%get_vxc(input%jspins,rho,v_xc, v_x,grad, kinED_KS=xcpot%kinED%is) else CALL xcpot%get_vxc(input%jspins,rho,v_xc,v_x,grad) endif - +#else + CALL xcpot%get_vxc(input%jspins,rho,v_xc,v_x,grad) +#endif IF (xcpot%needs_grad()) THEN SELECT TYPE(xcpot) TYPE IS (t_xcpot_libxc) @@ -82,12 +84,15 @@ CONTAINS !calculate the ex.-cor energy density IF (ALLOCATED(exc%pw_w)) THEN ALLOCATE ( e_xc(SIZE(rho,1),1) ); e_xc=0.0 - +#ifdef CPP_LIBXC IF(xcpot%kinED%set) THEN CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, xcpot%kinED%is, mt_call=.False.) ELSE CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, mt_call=.False.) ENDIF +#else + CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, mt_call=.False.) +#endif CALL pw_from_grid(xcpot,stars,.TRUE.,e_xc,exc%pw,exc%pw_w) ENDIF diff --git a/vgen/vmt_xc.F90 b/vgen/vmt_xc.F90 index 34212c9922748273d3fa2fbb431b728db24fc8f3..67a95e900053db5c9653b6df37631aec83754473 100644 --- a/vgen/vmt_xc.F90 +++ b/vgen/vmt_xc.F90 @@ -113,13 +113,18 @@ ! ! calculate the ex.-cor. potential +#ifdef CPP_LIBXC if(perform_MetaGGA .and. xcpot%kinED%set) then - CALL xcpot%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),v_xc(:nsp*atoms%jri(n),:)& + CALL xcpot%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),v_xc(:nsp*atoms%jri(n),:)& , v_x(:nsp*atoms%jri(n),:),grad, kinED_KS=xcpot%kinED%mt(:,:,loc_n)) else CALL xcpot%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),v_xc(:nsp*atoms%jri(n),:)& , v_x(:nsp*atoms%jri(n),:),grad) endif +#else + CALL xcpot%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),v_xc(:nsp*atoms%jri(n),:)& + , v_x(:nsp*atoms%jri(n),:),grad) +#endif IF (lda_atom(n)) THEN ! Use local part of pw91 for this atom CALL xcpot_tmp%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),xcl(:nsp*atoms%jri(n),:),v_x(:nsp*atoms%jri(n),:),grad) @@ -147,7 +152,7 @@ ! ! calculate the ex.-cor energy density ! - +#ifdef CPP_LIBXC IF(perform_MetaGGA .and. xcpot%kinED%set) THEN CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),& e_xc(:nsp*atoms%jri(n),1),grad, & @@ -156,7 +161,10 @@ CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),& e_xc(:nsp*atoms%jri(n),1),grad, mt_call=.True.) ENDIF - +#else + CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),& + e_xc(:nsp*atoms%jri(n),1),grad, mt_call=.True.) +#endif !write (*,*) "cut first ", cut_ratio, " number of points" !where(cut_mask) e_xc(:,1) = 0.0 diff --git a/xc-pot/metagga.F90 b/xc-pot/metagga.F90 index ee4357c158021b0756a1962183e4e6baf0f44623..80f66cbc80b465cd1dbec6f75ec04cf05febfd7d 100644 --- a/xc-pot/metagga.F90 +++ b/xc-pot/metagga.F90 @@ -145,7 +145,7 @@ CONTAINS implicit none INTEGER, intent(in) :: dim_idx TYPE (t_mat), intent(in) :: zMat - REAL, intent(in) :: kpt(3) + REAL, intent(in) :: kpt(3) TYPE(t_lapw), intent(in) :: lapw TYPE(t_cell), intent(in) :: cell TYPE (t_mat) :: zPrime @@ -162,9 +162,9 @@ CONTAINS fac = k_plus_g(dim_idx) if(zPrime%l_real) then - zPrime%data_r(basis_idx,:) = fac * zMat%data_r(basis_idx,:) + zPrime%data_r(basis_idx,:) = fac * zMat%data_r(basis_idx,:) else - zPrime%data_c(basis_idx,:) = ImagUnit * fac * zMat%data_c(basis_idx,:) + zPrime%data_c(basis_idx,:) = ImagUnit * fac * zMat%data_c(basis_idx,:) endif enddo end subroutine set_zPrime @@ -188,7 +188,7 @@ CONTAINS type(t_atoms), intent(in) :: atoms type(t_noco), intent(in) :: noco type(t_stars), intent(in) :: stars - + integer :: js, n, st do js = 1,size(vtot%mt,4) @@ -224,17 +224,18 @@ CONTAINS TYPE(t_stars),INTENT(IN) :: stars TYPE(t_cell),INTENT(IN) :: cell TYPE(t_potden),INTENT(IN) :: den, EnergyDen, vTot - + TYPE(t_potden) :: vTot_corrected - +#ifdef CPP_LIBXC call vTot_corrected%copyPotDen(vTot) call undo_vgen_finalize(vTot_corrected, atoms, noco, stars) call set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot_corrected) call set_kinED_mt(mpi, sphhar, atoms, sym, core_den, val_den, & xcpot, EnergyDen, input, vTot_corrected) +#endif end subroutine set_kinED - +#ifdef CPP_LIBXC subroutine set_kinED_is(xcpot, input, noco, stars, sym, cell, den, EnergyDen, vTot) use m_types use m_pw_tofrom_grid @@ -250,7 +251,7 @@ CONTAINS !local arrays REAL, ALLOCATABLE :: den_rs(:,:), ED_rs(:,:), vTot_rs(:,:) TYPE(t_gradients) :: tmp_grad - + CALL init_pw_grid(xcpot,stars,sym,cell) CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, & @@ -261,7 +262,7 @@ CONTAINS cell, den%pw, tmp_grad, den_rs) CALL finish_pw_grid() - + call calc_kinEnergyDen_pw(ED_rs, vTot_rs, den_rs, xcpot%kinED%is) !xcpot%kinED%is = ED_RS - vTot_RS * den_RS xcpot%kinED%set = .True. @@ -312,7 +313,7 @@ CONTAINS lbound(vTot%mt, dim=2):ubound(vTot%mt, dim=2),& lbound(vTot%mt, dim=4):ubound(vTot%mt, dim=4))) endif - + do jr=1,atoms%jri(n) vTot_mt(jr,0:,:) = vTot%mt(jr,0:,n,:) * atoms%rmsh(jr,n)**2 enddo @@ -320,7 +321,7 @@ CONTAINS n, tmp_grad, ED_rs) CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, vTot_mt(:,0:,:), & n, tmp_grad, vTot_rs) - + tmp_sphhar%nlhd = sphhar%nlhd tmp_sphhar%nlh = [(0, cnt=1,size(sphhar%nlh))] @@ -330,7 +331,7 @@ CONTAINS core_den%mt(:,0:,n,:), n, tmp_grad, core_den_rs) CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, & val_den%mt(:,0:,n,:), n, tmp_grad, val_den_rs) - + call calc_kinEnergyDen_mt(ED_RS, vTot_rs, vTot0_rs, core_den_rs, val_den_rs, & xcpot%kinED%mt(:,:,loc_n)) !xcpot%kinED%mt(:,:,loc_n) = ED_RS - (vTot0_rs * core_den_rs + vTot_rs * val_den_rs) @@ -338,4 +339,5 @@ CONTAINS xcpot%kinED%set = .True. CALL finish_mt_grid() end subroutine set_kinED_mt +#endif END MODULE m_metagga