Commit 5b72e8c9 authored by Daniel Wortmann's avatar Daniel Wortmann

Changes to vgen now compile, tests do not work yet...

parent 0725c6e9
......@@ -58,9 +58,10 @@ CONTAINS
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(INOUT) :: atoms!in u_setup n_u might be modified
TYPE(t_potden),INTENT(IN) :: inden
TYPE(t_potden),INTENT(INOUT) :: v,vx
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: inden,vx
TYPE(t_potden),INTENT(INOUT) :: v !u_setup will modify the potential matrix
#ifdef CPP_MPI
INCLUDE 'mpif.h'
#endif
......
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_constants
IMPLICIT NONE
INTEGER, PARAMETER :: noState_const = 0
INTEGER, PARAMETER :: coreState_const = 1
INTEGER, PARAMETER :: valenceState_const = 2
REAL, PARAMETER :: pi_const=3.1415926535897932
REAL, PARAMETER :: tpi_const=2.*3.1415926535897932
REAL, PARAMETER :: fpi_const=4.*3.1415926535897932
REAL, PARAMETER :: sfp_const=sqrt(4.*3.1415926535897932)
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(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',
& ' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se',
& 'Br','Kr','Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd',
& 'Ag','Cd','In','Sn','Sb','Te',' I','Xe','Cs','Ba','La','Ce',
& 'Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb',
& '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)',
+ '(4p1/2)','(4p3/2)','(5s1/2)','(4d3/2)','(4d5/2)',
+ '(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()
!
! 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)
!
! This subprogram supplies the value of c according to
! 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
RETURN
END FUNCTION c_light
END MODULE m_constants
......@@ -55,7 +55,7 @@ CONTAINS
!The total nucleii charge
zc=SUM(atoms%neq(:)*atoms%zatom(:))
zc = zc + 2*input%efield%sigma
zc = zc + 2*input%sigma
IF (fixtotal) THEN
!-roa
......
......@@ -3,9 +3,9 @@
USE m_constants
IMPLICIT NONE
PRIVATE
PUBLIC :: efield
PUBLIC :: e_field
CONTAINS
SUBROUTINE efield(atoms, dimension, stars, sym, vacuum, cell, input)
SUBROUTINE e_field(atoms, DIMENSION, stars, sym, vacuum, cell, input,efield)
!
!*********************************************************************
! sets the values of the sheets of charge for external electric
......@@ -28,7 +28,8 @@
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_input),INTENT(INOUT) :: input
TYPE(t_input),INTENT(IN) :: input
TYPE(t_efield),INTENT(INOUT) :: efield
! ..
! ..
! .. Local Scalars ..
......@@ -74,7 +75,7 @@
WRITE (6, '(3x,a,f12.5)') 'total electronic charge =', qe
WRITE (6, '(3x,a,f12.5)') 'total nuclear charge =', qn
CALL read_efield (input%efield, stars%mx1, stars%mx2, vacuum%nvac, cell%area)
CALL read_efield (efield, stars%mx1, stars%mx2, vacuum%nvac, cell%area)
! Sign convention of electric field: E>0 repels electrons,
! consistent with conventional choice that current is
......@@ -83,13 +84,13 @@
! In case of Dirichlet boundary conditions, we ignore the
! value of "sigma" and take the surplus charge into account
! in vvac.
if (input%efield%autocomp .or. input%efield%dirichlet) input%efield%sigma = 0.5*(qe-qn)
if (efield%autocomp .or. efield%dirichlet) efield%sigma = 0.5*(qe-qn)
CALL print_efield (6, input%efield, cell%area, vacuum%nvac, cell%amat,vacuum%dvac)
CALL print_efield (6, efield, cell%area, vacuum%nvac, cell%amat,vacuum%dvac)
IF (.NOT. input%efield%dirichlet&
& .AND. ABS (SUM (input%efield%sig_b) + 2*input%efield%sigma - (qe-qn)) > eps) THEN
IF (ABS (SUM (input%efield%sig_b) - (qe-qn)) < eps) THEN
IF (.NOT. efield%dirichlet&
& .AND. ABS (SUM (efield%sig_b) + 2*efield%sigma - (qe-qn)) > eps) THEN
IF (ABS (SUM (efield%sig_b) - (qe-qn)) < eps) THEN
CALL juDFT_error&
& ("E field: top+bottom+film charge does not add up to "&
& //"zero.",&
......@@ -104,29 +105,29 @@
& ,calledby ="efield")
ENDIF
IF (ABS (input%efield%sigma) > 0.49 .OR. ANY (ABS (input%efield%sig_b) > 0.49)) THEN
IF (ABS (efield%sigma) > 0.49 .OR. ANY (ABS (efield%sig_b) > 0.49)) THEN
WRITE ( 6,*) 'If you really want to calculate an e-field this'
WRITE ( 6,*) 'big, uncomment STOP in efield.f !'
CALL juDFT_error("E-field too big or No. of e- not correct"&
& ,calledby ="efield")
ENDIF
IF (input%efield%l_segmented) THEN
IF (efield%l_segmented) THEN
CALL V_seg_EF(&
& input,&
& efield,&
& vacuum, stars)
IF (input%efield%plot_rho)&
IF (efield%plot_rho)&
& CALL print_rhoEF(&
& input%efield, stars%mx1, stars%mx2, vacuum%nvac, stars%ng2, sym%nop, sym%nop2,&
& efield, stars%mx1, stars%mx2, vacuum%nvac, stars%ng2, sym%nop, sym%nop2,&
& stars%ng2, stars%kv2, sym%mrot, sym%symor, sym%tau, sym%invtab,&
& cell%area, stars%nstr2, cell%amat)
END IF
IF (ALLOCATED (input%efield%sigEF)) DEALLOCATE (input%efield%sigEF)
IF (ALLOCATED (efield%sigEF)) DEALLOCATE (efield%sigEF)
IF (input%efield%dirichlet .AND. ALLOCATED (input%efield%rhoEF))&
& call set_dirchlet_coeff (input%efield, vacuum%dvac, stars%ng2, stars%sk2, vacuum%nvac)
IF (efield%dirichlet .AND. ALLOCATED (efield%rhoEF))&
& call set_dirchlet_coeff (efield, vacuum%dvac, stars%ng2, stars%sk2, vacuum%nvac)
CONTAINS
......@@ -797,13 +798,12 @@
END SUBROUTINE print_efield
SUBROUTINE V_seg_EF(&
& input,&
& vacuum, stars&
& )
& efield,&
& vacuum, stars)
USE m_fft2d
use m_types
! Dummy variables:
TYPE(t_input), INTENT(INOUT) :: input
TYPE(t_efield), INTENT(INOUT) :: efield
TYPE(t_vacuum), INTENT(IN) :: vacuum
TYPE(t_stars), INTENT(IN) :: stars
......@@ -812,17 +812,17 @@
REAL :: fg, fgi
REAL :: rhoRS(3*stars%mx1,3*stars%mx2), rhoRSimag(3*stars%mx1,3*stars%mx2) ! Real space density
ALLOCATE (input%efield%rhoEF (stars%ng2-1, vacuum%nvac))
ALLOCATE (efield%rhoEF (stars%ng2-1, vacuum%nvac))
DO ivac = 1, vacuum%nvac
rhoRSimag = 0.0 ! Required in the loop. fft2d overrides this
! The fft2d algorithm puts the normalization to the isn=-1
! (r->k) transformation, thus we need to multiply by
! ifft2 = 3*k1d*3*k2d to compensate.
IF (input%efield%dirichlet) THEN
rhoRS(:,:) = input%efield%sigEF(:,:,ivac)
IF (efield%dirichlet) THEN
rhoRS(:,:) = efield%sigEF(:,:,ivac)
ELSE
rhoRS(:,:) = input%efield%sigEF(:,:,ivac)*(3*stars%mx1*3*stars%mx2)
rhoRS(:,:) = efield%sigEF(:,:,ivac)*(3*stars%mx1*3*stars%mx2)
END IF
! FFT rhoRS(r_2d) -> rhoEF(g_2d)
......@@ -830,16 +830,16 @@
& stars,&
& rhoRS, rhoRSimag,&
& fg, fgi,&
& input%efield%rhoEF(:,ivac), 1, -1)
& efield%rhoEF(:,ivac), 1, -1)
! FFT gives the the average charge per grid point
! while sig_b stores the (total) charge per sheet
IF (input%efield%dirichlet .and. ABS (fg) > 1.0e-15) THEN
IF (efield%dirichlet .and. ABS (fg) > 1.0e-15) THEN
PRINT *, 'INFO: Difference of average potential: fg=',&
& fg,', sig_b=', input%efield%sig_b(ivac),&
& fg,', sig_b=', efield%sig_b(ivac),&
& ", ivac=", ivac
ELSE IF (ABS (fg/(3*stars%mx1*3*stars%mx2)) > 1.0e-15) THEN
PRINT *, 'INFO: Difference of average potential: fg=',&
& fg/(3*stars%mx1*3*stars%mx2),', sig_b=', input%efield%sig_b(ivac),&
& fg/(3*stars%mx1*3*stars%mx2),', sig_b=', efield%sig_b(ivac),&
& ", ivac=", ivac
END IF
END DO ! ivac
......@@ -928,7 +928,7 @@
CLOSE (754)
END DO ! ivac
END SUBROUTINE print_rhoEF
END SUBROUTINE efield
END SUBROUTINE e_field
SUBROUTINE read_namelist (iou, E, eV)
USE m_types, only: t_efield
......@@ -978,4 +978,4 @@
END IF
END do
END FUNCTION lower_case
END MODULE m_efield
END MODULE m_efield
......@@ -185,7 +185,7 @@
!
!---> set up electric field parameters (if needed)
!
CALL efield(atoms, DIMENSION, stars, sym, vacuum, cell, input)
! CALL e_field(atoms, DIMENSION, stars, sym, vacuum, cell, input,field)
ENDIF
ENDIF
......
......@@ -8,7 +8,7 @@ MODULE m_postprocessInput
CONTAINS
SUBROUTINE postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts,&
oneD,hybrid,cell,banddos,sliceplot,xcpot,forcetheo,&
noco,DIMENSION,enpara,sphhar,l_opti,noel,l_kpts)
......@@ -62,6 +62,7 @@ SUBROUTINE postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
TYPE(t_dimension),INTENT(INOUT) :: dimension
TYPE(t_enpara) ,INTENT(INOUT) :: enpara
TYPE(t_sphhar) ,INTENT (OUT) :: sphhar
TYPE(t_field), INTENT(INOUT) :: field
LOGICAL, INTENT (OUT) :: l_opti
LOGICAL, INTENT (IN) :: l_kpts
CHARACTER(len=3), ALLOCATABLE, INTENT(IN) :: noel(:)
......@@ -555,7 +556,7 @@ SUBROUTINE postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
CALL stepf(sym,stars,atoms,oneD,input,cell,vacuum,mpi)
IF (mpi%irank.EQ.0) THEN
CALL convn(DIMENSION,atoms,stars)
CALL efield(atoms,DIMENSION,stars,sym,vacuum,cell,input)
CALL e_field(atoms,DIMENSION,stars,sym,vacuum,cell,input,field%efield)
END IF !(mpi%irank.EQ.0)
END IF
......
......@@ -639,22 +639,22 @@ SUBROUTINE r_inpXML(&
numberNodes = xmlGetNumberOfNodes(xPathA)
IF (numberNodes.EQ.1) THEN
input%efield%zsigma = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zsigma'))
input%efield%sig_b(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_1'))
input%efield%sig_b(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_2'))
input%efield%plot_charge = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_charge'))
input%efield%plot_rho = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_rho'))
input%efield%autocomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@autocomp'))
input%efield%dirichlet = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dirichlet'))
l_eV = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eV'))
!input%efield%zsigma = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zsigma'))
!input%efield%sig_b(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_1'))
!input%efield%sig_b(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_2'))
!input%efield%plot_charge = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_charge'))
!input%efield%plot_rho = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_rho'))
!input%efield%autocomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@autocomp'))
!input%efield%dirichlet = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dirichlet'))
!l_eV = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eV'))
STOP 'Error: Reading input for E-Fields not yet implemented completely!'
! ALLOCATE(input%efield%sigEF(3*k1d, 3*k2d, nvac))
! input%efield%sigEF = 0.0
IF (l_eV) THEN
input%efield%sig_b(:) = input%efield%sig_b/hartree_to_ev_const
!IF (l_eV) THEN
! input%efield%sig_b(:) = input%efield%sig_b/hartree_to_ev_const
! input%efield%sigEF(:,:,:) = input%efield%sigEF/hartree_to_ev_const
END IF
!END IF
END IF
! Read in optional energy parameter limits
......
......@@ -76,6 +76,7 @@ CONTAINS
! Types, these variables contain a lot of data!
TYPE(t_input) :: input
TYPE(t_field) :: field
TYPE(t_dimension):: DIMENSION
TYPE(t_atoms) :: atoms
TYPE(t_sphhar) :: sphhar
......@@ -113,7 +114,7 @@ CONTAINS
mpi%mpi_comm = mpi_comm
CALL timestart("Initialization")
CALL fleur_init(mpi,input,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,&
CALL fleur_init(mpi,input,field,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,&
sliceplot,banddos,obsolete,enpara,xcpot,results,kpts,hybrid,&
oneD,coreSpecInput,wann,l_opti)
CALL timestop("Initialization")
......@@ -145,8 +146,8 @@ CONTAINS
IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('scfLoop')
! Initialize and load inDen density (start)
CALL inDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
CALL inDenRot%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
CALL inDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
CALL inDenRot%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
archiveType = CDN_ARCHIVE_TYPE_CDN1_const
IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
......@@ -163,10 +164,10 @@ CONTAINS
! Initialize and load inDen density (end)
! Initialize potentials (start)
CALL vTot%init(stars,atoms,sphhar,vacuum,noco,oneD,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTTOT)
CALL vCoul%init(stars,atoms,sphhar,vacuum,noco,oneD,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTCOUL)
CALL vx%init(stars,atoms,sphhar,vacuum,noco,oneD,DIMENSION%jspd,.FALSE.,POTDEN_TYPE_POTCOUL)
CALL vTemp%init(stars,atoms,sphhar,vacuum,noco,oneD,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTTOT)
CALL vTot%init(stars,atoms,sphhar,vacuum,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTTOT)
CALL vCoul%init(stars,atoms,sphhar,vacuum,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTCOUL)
CALL vx%init(stars,atoms,sphhar,vacuum,DIMENSION%jspd,.FALSE.,POTDEN_TYPE_POTCOUL)
CALL vTemp%init(stars,atoms,sphhar,vacuum,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTTOT)
! Initialize potentials (end)
scfloop:DO WHILE (l_cont)
......@@ -244,7 +245,7 @@ CONTAINS
!---< gwf
CALL timestart("generation of potential")
CALL vgen(hybrid,reap,input,xcpot,DIMENSION, atoms,sphhar,stars,vacuum,&
CALL vgen(hybrid,field,input,xcpot,DIMENSION, atoms,sphhar,stars,vacuum,&
sym,obsolete,cell, oneD,sliceplot,mpi ,results,noco,inDen,inDenRot,vTot,vx,vCoul)
CALL timestop("generation of potential")
......@@ -380,7 +381,7 @@ CONTAINS
!!$ !-Wannier
CALL timestart("generation of new charge density (total)")
CALL outDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
CALL outDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
outDen%iter = inDen%iter
CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
DIMENSION,kpts,atoms,sphhar,stars,sym,obsolete,&
......
......@@ -7,7 +7,7 @@
IMPLICIT NONE
CONTAINS
SUBROUTINE fleur_init(mpi,&
input,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,&
input,field,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,&
sliceplot,banddos,obsolete,enpara,xcpot,results,kpts,hybrid,&
oneD,coreSpecInput,wann,l_opti)
USE m_judft
......@@ -43,6 +43,7 @@
! Types, these variables contain a lot of data!
TYPE(t_mpi) ,INTENT(INOUT):: mpi
TYPE(t_input) ,INTENT(OUT):: input
TYPE(t_field), INTENT(OUT) :: field
TYPE(t_dimension),INTENT(OUT):: DIMENSION
TYPE(t_atoms) ,INTENT(OUT):: atoms
TYPE(t_sphhar) ,INTENT(OUT):: sphhar
......@@ -185,7 +186,7 @@
results%force(:,:,:) = 0.0
END IF
CALL postprocessInput(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,&
CALL postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts,&
oneD,hybrid,cell,banddos,sliceplot,xcpot,forcetheo,&
noco,dimension,enpara,sphhar,l_opti,noel,l_kpts)
......
......@@ -6,8 +6,8 @@
MODULE m_vgen
USE m_juDFT
CONTAINS
SUBROUTINE vgen(hybrid,reap,input,xcpot,DIMENSION, atoms,sphhar,stars,&
vacuum,sym, obsolete,cell,oneD,sliceplot,mpi, results,noco,den,denRot,vTot,vx,vCoul)
SUBROUTINE vgen(hybrid,field,input,xcpot,DIMENSION, atoms,sphhar,stars,&
vacuum,sym,obsolete,cell,oneD,sliceplot,mpi, results,noco,den,denRot,vTot,vx,vCoul)
! ***********************************************************
! FLAPW potential generator *
! ***********************************************************
......@@ -17,34 +17,11 @@ CONTAINS
! TE_VEFF: charge density-effective potential integral
! TE_EXC : charge density-ex-corr.energy density integral
! ***********************************************************
#include"cpp_double.h"
USE m_constants
USE m_vmts
USE m_intnv
USE m_vmtxcg
USE m_vmtxc
USE m_vvacxc
USE m_vvacxcg
USE m_visxc
USE m_visxcg
USE m_vvac
USE m_vvacis
USE m_vvacxy
USE m_vintcz
USE m_checkdopall
USE m_wrtdop
USE m_cdn_io
USE m_qfix
USE m_bfield
USE m_vgen_coulomb
USE m_vgen_xcpot
USE m_vgen_finalize
USE m_types
USE m_od_vvac
USE m_od_vvacis
USE m_convol
USE m_xyavden
USE m_psqpw
USE m_potmod
USE m_intgr, ONLY : intgr3
USE m_cfft
USE m_fleur_vdw
#ifdef CPP_MPI
USE m_mpi_bc_potden
#endif
......@@ -57,68 +34,20 @@ CONTAINS
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_obsolete),INTENT(IN) :: obsolete
TYPE(t_sliceplot),INTENT(IN) :: sliceplot
TYPE(t_input),INTENT(INOUT) :: input !efield can be modified
TYPE(t_input),INTENT(IN) :: input
TYPE(t_field),INTENT(INOUT) :: field !efield can be modified
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(INOUT) :: atoms !vr0 is updated
TYPE(t_potden), INTENT(IN) :: den, denRot
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden), INTENT(INOUT) :: den, denRot
TYPE(t_potden),INTENT(INOUT) :: vTot,vx,vCoul
! ..
! .. Scalar Arguments ..
LOGICAL, INTENT (IN) :: reap
! .. Local type instances ..
TYPE(t_potden) :: workDen
! .. Local Scalars ..
COMPLEX vintcza,xint,rhobar
INTEGER i,i3,irec2,irec3,ivac,j,js,k,k3,lh,n,nzst1
INTEGER imz,imzxy,ichsmrg,ivfft
INTEGER ifftd,ifftd2, ifftxc3d,datend
INTEGER itypsym,itype,jsp,l,nat
! INTEGER i_sm,n_sm,i_sta,i_end
REAL ani,g3,z,rhmn,mfie,fermiEnergyTemp
REAL sig1dh,vz1dh,zat_l(atoms%ntype),rdum,dpdot ! ,delta,deltb,corr
LOGICAL l_pottot,l_vdw,l_qfix
LOGICAL exi
LOGICAL, PARAMETER :: l_xyav=.FALSE.
! ..
! .. Local Arrays ..
COMPLEX, ALLOCATABLE :: alphm(:,:)
COMPLEX, ALLOCATABLE :: excpw(:),excxy(:,:,:),vpw_w(:,:),psq(:)
REAL, ALLOCATABLE :: vbar(:),af1(:),bf1(:)
REAL, ALLOCATABLE :: rhoc(:,:,:),rhoc_vx(:)
REAL, ALLOCATABLE :: tec(:,:), qintc(:,:)
!.....potential
REAL, ALLOCATABLE :: vr_exx(:,:,:,:)
REAL, ALLOCATABLE :: veffr(:,:,:,:)
COMPLEX, ALLOCATABLE :: vpw_exx(:,:),vpw_wexx(:,:)
COMPLEX, ALLOCATABLE :: vxpw_w(:,:),veffpw_w(:,:)
!.....energy density
REAL, ALLOCATABLE :: excz(:,:),excr(:,:,:)
#ifdef CPP_MPI
include 'mpif.h'
integer:: ierr
#endif
!
! if you want to calculate potential gradients
!
!pg COMPLEX vlm(-lmaxd:lmaxd,0:lmaxd,ntypd)
!pg REAL fint,f(jmtd)
!pg INTEGER ns, nl, l, jm, m, mb
! ..
! ..
! ..
! ----> note the following conventions:
! ivac=2: lower (negative z) vacuum
! ivac=1: upper (positive z) vacuum
! units: hartrees
!
TYPE(t_potden) :: workden
WRITE (6,FMT=8000)
8000 FORMAT (/,/,t10,' p o t e n t i a l g e n e r a t o r',/)
......@@ -128,438 +57,29 @@ CONTAINS
CALL vCoul%resetPotDen()
CALL vx%resetPotDen()
ALLOCATE ( alphm(stars%ng2,2),excpw(stars%ng3),excxy(vacuum%nmzxyd,oneD%odi%n2d-1,2),&
vbar(dimension%jspd),af1(3*stars%mx3),bf1(3*stars%mx3),&
vpw_exx(stars%ng3,dimension%jspd),vpw_wexx(stars%ng3,dimension%jspd),&
excz(vacuum%nmzd,2),excr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype),&
vpw_w(stars%ng3,dimension%jspd),vxpw_w(stars%ng3,dimension%jspd),psq(stars%ng3) )
vTot%iter = den%iter
CALL workDen%init(stars,atoms,sphhar,vacuum,noco,oneD,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
workDen = den
IF (mpi%irank == 0) THEN
vxpw_w = CMPLX(0.,0.)
! ---> perform spin summation of charge densities
! ---> for the calculation of the coulomb potentials
IF ( (input%jspins.EQ.2).AND.(.NOT.l_xyav) ) THEN
workDen%mt(:,0:,:,1)=workDen%mt(:,0:,:,1)+workDen%mt(:,0:,:,2)
workDen%pw(:,1)=workDen%pw(:,1)+workDen%pw(:,2)
IF (input%film) THEN
workDen%vacxy(:,:,:,1)=workDen%vacxy(:,:,:,1)+workDen%vacxy(:,:,:,2)
workDen%vacz(:,:,1) = workDen%vacz(:,:,1) + workDen%vacz(:,:,2)
END IF
END IF
!
! ************** coulomb potential ***********************
CALL workDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,0)
! ----> create pesudo-charge density coefficients
ENDIF ! (mpi%irank == 0)
CALL vgen_coulomb(1,mpi,DIMENSION,oneD,input,vacuum,sym,stars,cell,sphhar,atoms,den,vCoul,results)
!sum up both spins in den into workden
CALL den%sum_both_spin(workden)
CALL vgen_coulomb(1,mpi,DIMENSION,oneD,input,field,vacuum,sym,stars,cell,sphhar,atoms,den,vCoul,results)
CALL vCoul%copy_both_spin(vTot)
call vgen_xcpot(hybrid,input,xcpot,DIMENSION, atoms,sphhar,stars,&
vacuum,sym, obsolete,cell,oneD,sliceplot,mpi,noco,den,