diff --git a/eigen/eigen.F90 b/eigen/eigen.F90 index dc75d0ce1e1b0825da738777019e83410556c450..2fedf5d77c522208aecff0329b2720edcc00f44e 100644 --- a/eigen/eigen.F90 +++ b/eigen/eigen.F90 @@ -42,7 +42,7 @@ CONTAINS IMPLICIT NONE TYPE(t_results),INTENT(INOUT):: results - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_dimension),INTENT(IN) :: DIMENSION TYPE(t_oneD),INTENT(IN) :: oneD diff --git a/force/geo.f90 b/force/geo.f90 index ff6565ca4c696dcc045c3e9a218c054e1e22f41c..cba12f633f3ef42c53fbeb51dca718c341d1a768 100644 --- a/force/geo.f90 +++ b/force/geo.f90 @@ -80,7 +80,7 @@ CONTAINS TYPE(t_banddos) :: banddos_temp TYPE(t_obsolete) :: obsolete_temp TYPE(t_enpara) :: enpara_temp - TYPE(t_xcpot) :: xcpot_temp + CLASS(t_xcpot),ALLOCATABLE :: xcpot_temp TYPE(t_results) :: results_temp TYPE(t_kpts) :: kpts_temp TYPE(t_hybrid) :: hybrid_temp diff --git a/hybrid/add_Vnonlocal.F90 b/hybrid/add_Vnonlocal.F90 index 98dc1b85747e160ca3d5410b807dacf6f6d423c5..929935ba9370b42b09c49bcfcb2b6f411818b137 100644 --- a/hybrid/add_Vnonlocal.F90 +++ b/hybrid/add_Vnonlocal.F90 @@ -49,7 +49,7 @@ MODULE m_add_vnonlocal USE m_io_hybrid IMPLICIT NONE TYPE(t_results),INTENT(INOUT) :: results - TYPE(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot TYPE(t_dimension),INTENT(IN) :: dimension TYPE(t_hybrid),INTENT(INOUT) :: hybrid TYPE(t_kpts),INTENT(IN) :: kpts diff --git a/hybrid/coulombmatrix.F90 b/hybrid/coulombmatrix.F90 index 8426ef8bcca06440f432719729dfa3dab0668d80..89253f06c67980807f902a317c175122320284ab 100644 --- a/hybrid/coulombmatrix.F90 +++ b/hybrid/coulombmatrix.F90 @@ -49,7 +49,7 @@ CONTAINS IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_hybrid),INTENT(INOUT) :: hybrid TYPE(t_sym),INTENT(IN) :: sym diff --git a/hybrid/exchange_val_hf.F90 b/hybrid/exchange_val_hf.F90 index 69042867d7d96e2f7f69b4c6bf77a223a12fd036..aa8668b7df04115442708ebb760ac5ad52997551 100644 --- a/hybrid/exchange_val_hf.F90 +++ b/hybrid/exchange_val_hf.F90 @@ -74,7 +74,7 @@ IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_results),INTENT(IN) :: results - TYPE(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_dimension),INTENT(IN) :: dimension TYPE(t_hybrid),INTENT(INOUT) :: hybrid diff --git a/hybrid/hsfock.F90 b/hybrid/hsfock.F90 index 732de8bda70be4273ef443c96de97e4657ede9eb..8ed7d5b562b4d2329e4b9f9fc75224666b86fda1 100644 --- a/hybrid/hsfock.F90 +++ b/hybrid/hsfock.F90 @@ -59,7 +59,7 @@ MODULE m_hsfock IMPLICIT NONE TYPE(t_hybdat),INTENT(IN) :: hybdat TYPE(t_results),INTENT(INOUT) :: results - TYPE(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_dimension),INTENT(IN) :: dimension TYPE(t_hybrid),INTENT(INOUT) :: hybrid diff --git a/hybrid/hybrid.F90 b/hybrid/hybrid.F90 index 289ab03108788d8c7899ea78f20982346d055845..d0736ea611e4fee9a7932e505d70733ef3205b6e 100644 --- a/hybrid/hybrid.F90 +++ b/hybrid/hybrid.F90 @@ -11,7 +11,7 @@ CONTAINS USE m_eig66_io USE m_io_hybrid IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_dimension),INTENT(IN) :: DIMENSION TYPE(t_oneD),INTENT(IN) :: oneD diff --git a/hybrid/mixedbasis.F90 b/hybrid/mixedbasis.F90 index bb6f9d37dfab1be4d011e1dc4a573b03b6f44f13..f3272319d4038938b9a2e7223528da90e916669d 100644 --- a/hybrid/mixedbasis.F90 +++ b/hybrid/mixedbasis.F90 @@ -44,7 +44,7 @@ CONTAINS USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_dimension),INTENT(IN) :: DIMENSION TYPE(t_hybrid),INTENT(INOUT) :: hybrid diff --git a/hybrid/subvxc.F90 b/hybrid/subvxc.F90 index 5c36f0ceb52433491aee7dde39e6fba46290e95a..24d4c3d1123ddf6a4ec425123c08d65c36e8a4de 100644 --- a/hybrid/subvxc.F90 +++ b/hybrid/subvxc.F90 @@ -14,7 +14,7 @@ CONTAINS USE m_abcof3 USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + TYPE(t_xcpot_inbuild),INTENT(IN) :: xcpot TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_dimension),INTENT(IN) :: DIMENSION TYPE(t_oneD),INTENT(IN) :: oneD diff --git a/init/initParallelProcesses.F90 b/init/initParallelProcesses.F90 index 0b647e85e90d34fe0bb4ee530c1c9a189669b744..a3610b2d152e4e9882e6342724ca15f875792c8a 100644 --- a/init/initParallelProcesses.F90 +++ b/init/initParallelProcesses.F90 @@ -36,7 +36,7 @@ SUBROUTINE initParallelProcesses(atoms,vacuum,input,stars,sliceplot,banddos,& TYPE(t_cell), INTENT(INOUT) :: cell TYPE(t_banddos), INTENT(INOUT) :: banddos TYPE(t_sliceplot),INTENT(INOUT) :: sliceplot - TYPE(t_xcpot), INTENT(INOUT) :: xcpot + CLASS(t_xcpot), INTENT(INOUT) :: xcpot TYPE(t_noco), INTENT(INOUT) :: noco TYPE(t_dimension),INTENT(INOUT) :: dimension TYPE(t_enpara), INTENT(INOUT) :: enpara diff --git a/init/od_strgn1.f90 b/init/od_strgn1.f90 index 4ac031ebb734b8f9e6f1ef24cc5d59665c6a0cde..d3205cd0b22504e540cc46a7098b8774e2197cd2 100644 --- a/init/od_strgn1.f90 +++ b/init/od_strgn1.f90 @@ -16,7 +16,7 @@ USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_cell),INTENT(IN) :: cell TYPE(t_sym),INTENT(IN) :: sym TYPE(t_oneD),INTENT(INOUT):: oneD diff --git a/init/old_inp/dimen7.F90 b/init/old_inp/dimen7.F90 index d8cc8e93d9a2adcfac38f1946f0d99c4a74edcad..3d15e0cf3244eccdf6ff250c50be47a1f3ad8115 100644 --- a/init/old_inp/dimen7.F90 +++ b/init/old_inp/dimen7.F90 @@ -49,7 +49,7 @@ TYPE(t_noco) :: noco TYPE(t_sliceplot) :: sliceplot TYPE(t_banddos) :: banddos - TYPE(t_xcpot) :: xcpot + TYPE(t_xcpot_inbuild) :: xcpot ! ! diff --git a/init/old_inp/fleur_init_old.F90 b/init/old_inp/fleur_init_old.F90 index 64dc11878b036d51b514f7c4455314e5cc0d395f..cf9f201be92dd1923180381d421d65e73fd9219c 100644 --- a/init/old_inp/fleur_init_old.F90 +++ b/init/old_inp/fleur_init_old.F90 @@ -38,7 +38,7 @@ CONTAINS TYPE(t_banddos) ,INTENT(OUT) :: banddos TYPE(t_obsolete) ,INTENT(OUT) :: obsolete TYPE(t_enpara) ,INTENT(OUT) :: enpara - TYPE(t_xcpot) ,INTENT(OUT) :: xcpot + TYPE(t_xcpot_inbuild) ,INTENT(OUT) :: xcpot TYPE(t_kpts) ,INTENT(INOUT):: kpts TYPE(t_hybrid) ,INTENT(OUT) :: hybrid TYPE(t_oneD) ,INTENT(OUT) :: oneD @@ -180,7 +180,7 @@ CONTAINS ! namex=xcpot%get_name() - l_krla = xcpot%krla.EQ.1 + l_krla = xcpot%data%krla.EQ.1 END IF ! mpi%irank.eq.0 #ifdef CPP_MPI diff --git a/init/old_inp/inped.F90 b/init/old_inp/inped.F90 index 032334aa8ac3c0e5b61ef9978086b282bbb104c5..b6af1015d3c94fc5163c151767ef6446788055ae 100644 --- a/init/old_inp/inped.F90 +++ b/init/old_inp/inped.F90 @@ -44,7 +44,7 @@ TYPE(t_vacuum), INTENT(INOUT) :: vacuum TYPE(t_input), INTENT(INOUT) :: input TYPE(t_banddos), INTENT(INOUT) :: banddos - TYPE(t_xcpot), INTENT(INOUT) :: xcpot + TYPE(t_xcpot_inbuild), INTENT(INOUT) :: xcpot TYPE(t_sym), INTENT(INOUT) :: sym TYPE(t_cell), INTENT(INOUT) :: cell TYPE(t_sliceplot), INTENT(INOUT) :: sliceplot diff --git a/init/old_inp/setup.f90 b/init/old_inp/setup.f90 index 9cc43cd0640dc9d35593e8506e69fd9b1f18d77c..60af5424938c898b202af01aa41d08f6cc0e6ce0 100644 --- a/init/old_inp/setup.f90 +++ b/init/old_inp/setup.f90 @@ -68,7 +68,7 @@ TYPE(t_noco),INTENT(INOUT) :: noco TYPE(t_vacuum),INTENT(INOUT) :: vacuum TYPE(t_cell),INTENT(INOUT) :: cell - TYPE(t_xcpot),INTENT(INOUT) :: xcpot + CLASS(t_xcpot),INTENT(INOUT) :: xcpot TYPE(t_sliceplot),INTENT(INOUT):: sliceplot TYPE(t_enpara),INTENT(INOUT) :: enpara LOGICAL, INTENT (IN) :: l_opti diff --git a/init/postprocessInput.F90 b/init/postprocessInput.F90 index 287480badf3a60c4b1b2e7672fbb661d2d93133e..41364c9281413a9b55f3acebc16ee8a730921ffa 100644 --- a/init/postprocessInput.F90 +++ b/init/postprocessInput.F90 @@ -57,7 +57,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts TYPE(t_cell), INTENT(INOUT) :: cell TYPE(t_banddos), INTENT(INOUT) :: banddos TYPE(t_sliceplot),INTENT(INOUT) :: sliceplot - TYPE(t_xcpot), INTENT(INOUT) :: xcpot + CLASS(t_xcpot), INTENT(INOUT) :: xcpot TYPE(t_noco), INTENT(INOUT) :: noco TYPE(t_dimension),INTENT(INOUT) :: dimension TYPE(t_enpara) ,INTENT(INOUT) :: enpara @@ -73,9 +73,8 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts REAL :: sumWeight, rmtmax, zp, radius, dr REAL :: kmax1, dtild1, dvac1 REAL :: bk(3) - LOGICAL :: l_vca, l_test,l_gga, l_krla - CHARACTER(len=4) :: namex - + LOGICAL :: l_vca, l_test,l_gga + INTEGER, ALLOCATABLE :: lmx1(:), nq1(:), nlhtp1(:) INTEGER, ALLOCATABLE :: jri1(:), lmax1(:) REAL, ALLOCATABLE :: rmt1(:), dx1(:) @@ -536,22 +535,13 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,obsolete,kpts END IF CALL prp_xcfft(stars,input,cell,xcpot) - - namex = xcpot%get_name() - l_krla = xcpot%krla.EQ.1 - + END IF !(mpi%irank.EQ.0) - + call xcpot%broadcast(mpi) #ifdef CPP_MPI - CALL MPI_BCAST(namex,4,MPI_CHARACTER,0,mpi%mpi_comm,ierr) - CALL MPI_BCAST(l_krla,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr) CALL MPI_BCAST(sliceplot%iplot,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr) #endif - IF (mpi%irank.NE.0) THEN - CALL xcpot%init(namex,l_krla) - END IF - IF (.NOT.sliceplot%iplot) THEN CALL stepf(sym,stars,atoms,oneD,input,cell,vacuum,mpi) IF (mpi%irank.EQ.0) THEN diff --git a/init/prp_xcfft.f90 b/init/prp_xcfft.f90 index 79d3c82fe8851bddcd913a11ac74b636b9eed2de..6989a3dca211925373897dba8d3c95480f37f6bb 100644 --- a/init/prp_xcfft.f90 +++ b/init/prp_xcfft.f90 @@ -32,7 +32,7 @@ TYPE(t_stars),INTENT(INOUT) :: stars TYPE(t_input),INTENT(IN) :: input TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_xcpot),INTENT(INOUT) :: xcpot + CLASS(t_xcpot),INTENT(INOUT) :: xcpot ! !---> local variables diff --git a/init/strgn.f90 b/init/strgn.f90 index ac62fa9d33d7d89e1f1d93127acc7867bc0ff00e..ed39359ec569fadcbd6605ae3bfb876e02daf26f 100644 --- a/init/strgn.f90 +++ b/init/strgn.f90 @@ -28,7 +28,7 @@ CONTAINS TYPE(t_sphhar),INTENT(IN) :: sphhar TYPE(t_input),INTENT(IN) :: input TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot ! .. ! .. @@ -564,7 +564,7 @@ CONTAINS TYPE(t_sphhar),INTENT(IN) :: sphhar TYPE(t_input),INTENT(IN) :: input TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot ! .. ! .. Local Scalars .. REAL arltv1,arltv2,arltv3,s diff --git a/inpgen/set_inp.f90 b/inpgen/set_inp.f90 index e2476739cc0c553a2957a3d57e5f371c1e950c83..0841a20a676688e84529ca9fc8ac3fd5b4b26a64 100644 --- a/inpgen/set_inp.f90 +++ b/inpgen/set_inp.f90 @@ -30,6 +30,7 @@ USE m_juDFT_init USE m_kpoints USE m_inv3 + USE m_types_xcpot_inbuild IMPLICIT NONE TYPE(t_input),INTENT(INOUT) :: input @@ -74,7 +75,7 @@ TYPE(t_oneD)::oneD TYPE(t_stars)::stars TYPE(t_hybrid)::hybrid - TYPE(t_xcpot)::xcpot + TYPE(t_xcpot_inbuild)::xcpot TYPE(t_kpts)::kpts TYPE(t_enpara)::enpara TYPE(t_forcetheo)::forcetheo diff --git a/io/r_inpXML.F90 b/io/r_inpXML.F90 index 693b1b8b5a6c4330b06e32269876ea39b29bd5dd..afafd99bd393525939bf8f19617d4ca63274455e 100644 --- a/io/r_inpXML.F90 +++ b/io/r_inpXML.F90 @@ -35,6 +35,8 @@ SUBROUTINE r_inpXML(& USE m_constants USE m_inpeig USE m_sort + USE m_types_xcpot_inbuild +! USE m_types_xcpot_libxc IMPLICIT NONE TYPE(t_input),INTENT(INOUT) :: input @@ -49,7 +51,7 @@ SUBROUTINE r_inpXML(& TYPE(t_cell),INTENT(INOUT) :: cell TYPE(t_banddos),INTENT(INOUT) :: banddos TYPE(t_sliceplot),INTENT(INOUT):: sliceplot - TYPE(t_xcpot),INTENT(INOUT) :: xcpot + CLASS(t_xcpot),INTENT(INOUT),ALLOCATABLE :: xcpot TYPE(t_noco),INTENT(INOUT) :: noco TYPE(t_dimension),INTENT(OUT) :: dimension TYPE(t_enpara) ,INTENT(OUT) :: enpara @@ -225,7 +227,6 @@ SUBROUTINE r_inpXML(& ALLOCATE(atoms%igrd(atoms%ntype)) ALLOCATE(atoms%krla(atoms%ntype)) ALLOCATE(atoms%relcor(atoms%ntype)) - ALLOCATE(xcpot%lda_atom(atoms%ntype)) atoms%namex = '' atoms%icorr = -99 @@ -293,12 +294,7 @@ SUBROUTINE r_inpXML(& input%rkmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Kmax')) stars%gmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Gmax')) - xPathA = '/fleurInput/calculationSetup/cutoffs/@GmaxXC' - numberNodes = xmlGetNumberOfNodes(xPathA) - xcpot%gmaxxc = stars%gmax - IF(numberNodes.EQ.1) THEN - xcpot%gmaxxc = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) - END IF + stars%gmaxInit = stars%gmax xPathA = '/fleurInput/calculationSetup/cutoffs/@numbands' @@ -1109,16 +1105,24 @@ SUBROUTINE r_inpXML(& ! Read in xc functional parameters - namex = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/@name'))))) + valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/@name'))))) l_relcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/@relativisticCorrections')) relcor = 'non-relativi' IF (l_relcor) THEN relcor = 'relativistic' END IF - - CALL getXCParameters(namex,l_relcor,xcpot,hybrid%l_hybrid) - + !now initialize the xcpot variable + CALL setXCParameters(atoms,valueString,l_relcor,xcpot) + + xPathA = '/fleurInput/calculationSetup/cutoffs/@GmaxXC' + numberNodes = xmlGetNumberOfNodes(xPathA) + xcpot%gmaxxc = stars%gmax + IF(numberNodes.EQ.1) THEN + xcpot%gmaxxc = evaluateFirstOnly(xmlGetAttributeValue(xPathA)) + END IF + hybrid%l_hybrid=xcpot%is_hybrid() + IF (hybrid%l_hybrid) ALLOCATE(hybrid%lcutm1(atoms%ntype),hybrid%lcutwf(atoms%ntype),hybrid%select1(4,atoms%ntype)) obsolete%lwb=.FALSE. @@ -1514,7 +1518,10 @@ SUBROUTINE r_inpXML(& hybrid%select1(:,iType)=SELECT ENDIF ! Explicit xc functional - xcpot%lda_atom(iType)=ldaSpecies + SELECT TYPE(xcpot) + TYPE IS(t_xcpot_inbuild) + xcpot%lda_atom(iType)=ldaSpecies + END SELECT noco%socscale(iType)=socscaleSpecies END IF END DO @@ -2000,24 +2007,34 @@ SUBROUTINE r_inpXML(& END SUBROUTINE r_inpXML -SUBROUTINE getXCParameters(namex,relcor,xcpot,l_hybrid) - +SUBROUTINE setXCParameters(atoms,namex,relcor,xcpot) USE m_juDFT USE m_types + USE m_types_xcpot_inbuild +! USE m_types_xcpot_libxc IMPLICIT NONE - - CHARACTER(LEN=4), INTENT(IN) :: namex + TYPE(t_atoms),INTENT(IN) :: atoms + CHARACTER(LEN=*), INTENT(IN) :: namex LOGICAL, INTENT(IN) :: relcor - TYPE(t_xcpot),INTENT(INOUT) :: xcpot - LOGICAL, INTENT(OUT) :: l_hybrid - - CALL xcpot%init(namex,relcor) - l_hybrid=xcpot%is_hybrid() - + CLASS(t_xcpot),INTENT(OUT),ALLOCATABLE :: xcpot + + IF (namex(1:6)=='libxc:') THEN +! ALLOCATE(t_xcpot_libxc::xcpot) + ELSE + ALLOCATE(t_xcpot_inbuild::xcpot) + ENDIF + + SELECT TYPE(xcpot) + TYPE IS(t_xcpot_inbuild) + CALL xcpot%init(namex(1:4),relcor) + ALLOCATE(xcpot%lda_atom(atoms%ntype)) +! TYPE IS(t_xcpot_libxc) +! CALL xcpot%init(namex) + END SELECT -END SUBROUTINE getXCParameters + END SUBROUTINE setXCParameters SUBROUTINE getIntegerSequenceFromString(string, sequence, count) diff --git a/io/rw_inp.f90 b/io/rw_inp.f90 index b4d83088e0ab3b8c674f7838bcee2176361ed863..fd74a9af34ab257cd94facf9a65fa351899467b3 100644 --- a/io/rw_inp.f90 +++ b/io/rw_inp.f90 @@ -37,7 +37,7 @@ TYPE(t_cell),INTENT(INOUT) :: cell TYPE(t_banddos),INTENT(INOUT) :: banddos TYPE(t_sliceplot),INTENT(INOUT):: sliceplot - TYPE(t_xcpot),INTENT(INOUT) :: xcpot + TYPE(t_xcpot_inbuild),INTENT(INOUT) :: xcpot TYPE(t_noco),INTENT(INOUT) :: noco REAL,INTENT(INOUT) :: a1(3),a2(3),a3(3) diff --git a/io/w_inpXML.f90 b/io/w_inpXML.f90 index 9911def79a05d6fb68e290e66430e0cbbd090587..2b539ab0c49770c8293c11cf54a8a1c476ea668f 100644 --- a/io/w_inpXML.f90 +++ b/io/w_inpXML.f90 @@ -45,7 +45,7 @@ SUBROUTINE w_inpXML(& TYPE(t_cell),INTENT(IN) :: cell TYPE(t_banddos),INTENT(IN) :: banddos TYPE(t_sliceplot),INTENT(IN):: sliceplot - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_noco),INTENT(IN) :: noco TYPE(t_enpara),INTENT(IN) :: enpara CLASS(t_forcetheo),INTENT(IN):: forcetheo !nothing is done here so far.... diff --git a/io/writeOutParameters.f90 b/io/writeOutParameters.f90 index 32d3269dd92d17599688b39ac4a6b242fec0ae61..1197b81d4b288b7e8d32d7e624da75a19910c0e6 100644 --- a/io/writeOutParameters.f90 +++ b/io/writeOutParameters.f90 @@ -24,7 +24,7 @@ SUBROUTINE writeOutParameters(mpi,input,sym,stars,atoms,vacuum,obsolete,kpts,& TYPE(t_cell), INTENT(IN) :: cell TYPE(t_banddos), INTENT(IN) :: banddos TYPE(t_sliceplot), INTENT(IN) :: sliceplot - TYPE(t_xcpot), INTENT(IN) :: xcpot + CLASS(t_xcpot), INTENT(IN) :: xcpot TYPE(t_noco), INTENT(IN) :: noco TYPE(t_dimension), INTENT(IN) :: dimension TYPE(t_enpara), INTENT(IN) :: enpara diff --git a/main/fleur.F90 b/main/fleur.F90 index 2fceb165ef92f769dfa88d850b30659f3c62da1f..c4505de9f954d978144d7b6ab9883b3d0d8f9e41 100644 --- a/main/fleur.F90 +++ b/main/fleur.F90 @@ -88,7 +88,6 @@ CONTAINS TYPE(t_banddos) :: banddos TYPE(t_obsolete) :: obsolete TYPE(t_enpara) :: enpara - TYPE(t_xcpot) :: xcpot TYPE(t_results) :: results TYPE(t_kpts) :: kpts TYPE(t_hybrid) :: hybrid @@ -98,6 +97,7 @@ CONTAINS TYPE(t_wann) :: wann TYPE(t_potden) :: vTot,vx,vCoul,vTemp TYPE(t_potden) :: inDen, outDen + CLASS(t_xcpot),ALLOCATABLE :: xcpot CLASS(t_forcetheo),ALLOCATABLE:: forcetheo ! .. Local Scalars .. @@ -207,8 +207,13 @@ CONTAINS !HF - IF (hybrid%l_hybrid) CALL calc_hybrid(hybrid,kpts,atoms,input,DIMENSION,mpi,noco,& - cell,vacuum,oneD,banddos,results,sym,xcpot,vTot,it) + IF (hybrid%l_hybrid) THEN + SELECT TYPE(xcpot) + TYPE IS(t_xcpot_inbuild) + CALL calc_hybrid(hybrid,kpts,atoms,input,DIMENSION,mpi,noco,& + cell,vacuum,oneD,banddos,results,sym,xcpot,vTot,it) + END SELECT + ENDIF !#endif !!$ DO pc = 1, wann%nparampts diff --git a/main/fleur_init.F90 b/main/fleur_init.F90 index 4a339acebb3d833e830b4cb4a208d0b13d91fadd..923245b8822f4db99060a65f7d26162567deec16 100644 --- a/main/fleur_init.F90 +++ b/main/fleur_init.F90 @@ -32,6 +32,7 @@ USE m_prpqfftmap USE m_writeOutHeader USE m_fleur_init_old + USE m_types_xcpot_inbuild #ifdef CPP_MPI USE m_mpi_bc_all, ONLY : mpi_bc_all USE m_mpi_dist_forcetheorem @@ -56,7 +57,7 @@ TYPE(t_banddos) ,INTENT(OUT):: banddos TYPE(t_obsolete) ,INTENT(OUT):: obsolete TYPE(t_enpara) ,INTENT(OUT):: enpara - TYPE(t_xcpot) ,INTENT(OUT):: xcpot + CLASS(t_xcpot),ALLOCATABLE,INTENT(OUT):: xcpot TYPE(t_results) ,INTENT(OUT):: results TYPE(t_kpts) ,INTENT(OUT):: kpts TYPE(t_hybrid) ,INTENT(OUT):: hybrid @@ -206,10 +207,14 @@ #endif ELSE ! else branch of "IF (input%l_inpXML) THEN" - CALL fleur_init_old(mpi,& - input,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,& - sliceplot,banddos,obsolete,enpara,xcpot,kpts,hybrid,& - oneD,coreSpecInput,l_opti) + ALLOCATE(t_xcpot_inbuild::xcpot) + SELECT TYPE(xcpot) + TYPE IS(t_xcpot_inbuild) + CALL fleur_init_old(mpi,& + input,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,& + sliceplot,banddos,obsolete,enpara,xcpot,kpts,hybrid,& + oneD,coreSpecInput,l_opti) + END SELECT END IF ! end of else branch of "IF (input%l_inpXML) THEN" ! diff --git a/main/optional.F90 b/main/optional.F90 index e8678edbd7e783204bddb8334a1126c7c5127f1d..80ee077250f29f2aac5e88b36fe7da5c119c033e 100644 --- a/main/optional.F90 +++ b/main/optional.F90 @@ -80,7 +80,7 @@ CONTAINS TYPE(t_noco),INTENT(IN) :: noco TYPE(t_vacuum),INTENT(IN) :: vacuum TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_sliceplot),INTENT(IN):: sliceplot ! .. ! .. Local Scalars .. diff --git a/main/totale.f90 b/main/totale.f90 index 9e08427051d556d504b7a7b99de2a962b13b00d5..2cb56ebc00154058d28a2832a9616f30cbe67d86 100644 --- a/main/totale.f90 +++ b/main/totale.f90 @@ -51,7 +51,7 @@ CONTAINS IMPLICIT NONE TYPE(t_results),INTENT(INOUT) :: results - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_oneD),INTENT(IN) :: oneD TYPE(t_hybrid),INTENT(IN) :: hybrid TYPE(t_input),INTENT(IN) :: input @@ -119,11 +119,11 @@ CONTAINS ! ---> Fock exchange contribution ! IF (xcpot%is_hybrid()) THEN - IF (xcpot%is_name("exx")) THEN - results%tote = results%tote + 0.5e0*results%te_hfex%valence - ELSE + !IF (xcpot%is_name("exx")) THEN + ! results%tote = results%tote + 0.5e0*results%te_hfex%valence + !ELSE results%tote = results%tote - 0.5e0*results%te_hfex%valence + 0.5e0*results%te_hfex%core - END IF + !END IF ENDIF WRITE (6,FMT=8100) 0.5e0*results%te_hfex%valence WRITE (16,FMT=8100) 0.5e0*results%te_hfex%valence diff --git a/main/vgen.F90 b/main/vgen.F90 index e3a36e30e47000c72ff1e0ce2b031f221a25deb7..b444d41652d2520931dfa7d4fbcee0a9d97903c2 100644 --- a/main/vgen.F90 +++ b/main/vgen.F90 @@ -29,7 +29,7 @@ CONTAINS #endif IMPLICIT NONE TYPE(t_results),INTENT(INOUT) :: results - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_hybrid),INTENT(IN) :: hybrid TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_dimension),INTENT(IN) :: dimension diff --git a/optional/atom2.f90 b/optional/atom2.f90 index 16416477f832a2480056d4f4a0f5080a1a38b119..eddd19e3cd37d734b741007adc0bbf07217ab5e2 100644 --- a/optional/atom2.f90 +++ b/optional/atom2.f90 @@ -15,7 +15,6 @@ USE m_intgr, ONLY : intgr1,intgr0 USE m_constants - USE m_xcall, ONLY : vxcall USE m_potl0 USE m_stpot1 USE m_setcor @@ -26,7 +25,7 @@ ! .. Scalar Arguments .. TYPE(t_dimension),INTENT(IN) :: dimension TYPE(t_atoms),INTENT(IN) :: atoms - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_input),INTENT(IN) :: input INTEGER,INTENT (IN) :: jrc ,ntyp REAL, INTENT (IN) :: rnot1 ,qdel @@ -224,10 +223,7 @@ xcpot,DIMENSION%msh,DIMENSION%jspd,input%jspins,n,& atoms%dx(ntyp),rad,rhoss, vxc) ELSE - CALL vxcall(6,xcpot,input%jspins,& - SIZE(vx,1),jrc,rhoss,& - vx,vxc) - + CALL xcpot%get_vxc(input%jspins,rhoss,vxc,vx) ENDIF DO ispin = 1, input%jspins DO i = 1,n diff --git a/optional/stden.f90 b/optional/stden.f90 index a013d4c5d81c6a1393c219ec95f6a488924b89c3..765956c2d26367c99345a13cd7c2db8312654bbc 100644 --- a/optional/stden.f90 +++ b/optional/stden.f90 @@ -16,7 +16,6 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,& input,cell,xcpot,obsolete,noco,oneD) USE m_constants - USE m_xcall, ONLY : vxcall USE m_qsf USE m_checkdopall USE m_cdnovlp @@ -24,6 +23,7 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,& USE m_qfix USE m_atom2 USE m_types + USE m_types_xcpot_inbuild USE m_juDFT_init IMPLICIT NONE @@ -40,16 +40,16 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,& TYPE(t_input),INTENT(IN) :: input TYPE(t_vacuum),INTENT(IN) :: vacuum TYPE(t_cell),INTENT(IN) :: cell - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot ! Local type instances TYPE(t_potden) :: den TYPE(t_enpara) :: enpara - TYPE(t_xcpot) :: xcpot_dummy + TYPE(t_xcpot_inbuild) :: xcpot_dummy ! Local Scalars REAL d,del,fix,h,r,rnot,z,bm,qdel - REAL denz1(1),vacxpot(1),vacpot(1) + REAL denz1(1,1),vacxpot(1,1),vacpot(1,1) INTEGER i,ivac,iza,j,jr,k,n,n1,ispin INTEGER nw,ilo,natot,nat @@ -308,17 +308,17 @@ SUBROUTINE stden(mpi,sphhar,stars,atoms,sym,DIMENSION,vacuum,& sigm(i) = (i-1)*vacuum%delz*den%vacz(i,ivac,ispin) END DO CALL qsf(vacuum%delz,sigm,vacpar(ivac),vacuum%nmz,0) - denz1(1) = den%vacz(1,ivac,ispin) ! get estimate for potential at vacuum boundary - CALL vxcall(6,xcpot_dummy,1,1,1,denz1,vacxpot,vacpot) + denz1 = den%vacz(1,ivac,ispin) ! get estimate for potential at vacuum boundary + CALL xcpot%get_vxc(1,denz1,vacpot,vacxpot) ! seems to be the best choice for 1D not to substract vacpar IF (.NOT.oneD%odi%d1) THEN - vacpot(1) = vacpot(1) - fpi_const*vacpar(ivac) + vacpot = vacpot - fpi_const*vacpar(ivac) END IF IF (obsolete%lepr.EQ.1) THEN - vacpar(ivac) = -0.2 - vacpot(1) + vacpar(ivac) = -0.2 - vacpot(1,1) WRITE (6,'(" vacuum",i2," reference energy =",f12.6)') ivac,vacpot ELSE - vacpar(ivac) = vacpot(1) + vacpar(ivac) = vacpot(1,1) END IF END DO IF (vacuum%nvac.EQ.1) vacpar(2) = vacpar(1) diff --git a/types/CMakeLists.txt b/types/CMakeLists.txt index 1bb4290d63dbbbcded12a4a8d59f543feffd12a9..31816fcbba688bb8b58df58f03a119918c74dec4 100644 --- a/types/CMakeLists.txt +++ b/types/CMakeLists.txt @@ -4,6 +4,9 @@ set(fleur_F90 ${fleur_F90} types/types.F90 types/types_mat.F90 types/types_xcpot.F90 +types/types_xcpot_inbuild.F90 +types/types_xcpot_data.F90 +#types/types_xcpot_libxc.F90 types/types_mpi.F90 types/types_lapw.F90 types/types_tlmplm.F90 @@ -24,6 +27,8 @@ types/types.F90 types/types_mat.F90 types/types_field.F90 types/types_xcpot.F90 +types/types_xcpot_data.F90 +types/types_xcpot_inbuild.F90 types/types_mpi.F90 types/types_lapw.F90 types/types_tlmplm.F90 diff --git a/types/types.F90 b/types/types.F90 index f42cd08151d0299067d56dbe1ff1229f93606ffb..0add429942fbb1e3ed486bd61c600c099182afe8 100644 --- a/types/types.F90 +++ b/types/types.F90 @@ -21,4 +21,5 @@ MODULE m_types USE m_types_forcetheo USE m_types_cdnval USE m_types_field + USE m_types_xcpot_inbuild END MODULE m_types diff --git a/types/types_xcpot.F90 b/types/types_xcpot.F90 index f65a7a669aa9bb66a028c46843ae09ac0772c7e3..f3a3f8846e155eab4f565b8eb4000f181862b1a8 100644 --- a/types/types_xcpot.F90 +++ b/types/types_xcpot.F90 @@ -1,140 +1,99 @@ +!-------------------------------------------------------------------------------- +! 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_types_xcpot IMPLICIT NONE PRIVATE - 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_gga(20)=[& - .TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,& - .TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,& - .TRUE.,.TRUE.,.TRUE.,.TRUE.,.FALSE.,.FALSE.] - - 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 + PUBLIC :: t_xcpot,t_gradients - TYPE t_xcpot -#ifdef CPP_MPI - INTEGER :: icorr=0 !not private to allow bcasting it around -#else - INTEGER,PRIVATE :: icorr=0 -#endif - - !in the pbe case (exchpbe.F) lots of test are made - !in addition some constants are set - !to speed up this code precalculate things in init - LOGICAL :: is_rpbe !Rpbe - LOGICAL :: is_wc - LOGICAL :: is_hse !hse,lhse,vhse - REAL :: uk,um - - LOGICAL,ALLOCATABLE :: lda_atom(:) - REAL :: gmaxxc - INTEGER :: krla !relativistic corrections - + TYPE,ABSTRACT :: t_xcpot + REAL :: gmaxxc CONTAINS PROCEDURE :: is_gga=>xcpot_is_gga - PROCEDURE :: get_name=>xcpot_get_name - PROCEDURE :: init=>xcpot_init - PROCEDURE :: is_hybrid=>xcpot_is_hybrid - PROCEDURE :: is_name=>xcpot_is_name + 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 + PROCEDURE :: broadcast=>xcpot_broadcast + PROCEDURE,NOPASS :: alloc_gradients=>xcpot_alloc_gradients END TYPE t_xcpot - PUBLIC t_xcpot -CONTAINS - CHARACTER(len=4) FUNCTION xcpot_get_name(xcpot) - USE m_judft - IMPLICIT NONE - CLASS(t_xcpot),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) - USE m_judft - IMPLICIT NONE - CLASS(t_xcpot),INTENT(INOUT) :: xcpot - CHARACTER(len=*),INTENT(IN) :: namex - LOGICAL,INTENT(IN) :: relcor - INTEGER:: n - !Determine icorr from name - 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%krla=MERGE(1,0,relcor) + TYPE t_gradients + !Naming convention: + !t,u,d as last letter for total,up,down + !agr for absolute value of gradient + !g2r for laplacien of gradient + !+?? + REAL,ALLOCATABLE :: agrt(:),agru(:),agrd(:) + REAL,ALLOCATABLE :: g2ru(:),g2rd(:),gggrt(:) + REAL,ALLOCATABLE :: gggru(:),gzgr(:),g2rt(:) + REAL,ALLOCATABLE :: gggrd(:),grgru(:),grgrd(:) + !These are the contracted Gradients used in libxc + REAL,ALLOCATABLE :: sigma(:,:) + END TYPE t_gradients - - !Code from exchpbe to speed up determination of constants - IF (xcpot%is_name("rpbe")) THEN - xcpot%uk=1.2450 - ELSE - xcpot%uk=0.8040 - ENDIF - IF (xcpot%is_name("PBEs")) THEN ! pbe_sol - xcpot%um=0.123456790123456d0 - ELSE - xcpot%um=0.2195149727645171e0 - ENDIF - xcpot%is_hse=xcpot%is_name("hse").OR.xcpot%is_name("lhse").OR.xcpot%is_name("vhse") - xcpot%is_rpbe=xcpot%is_name("Rpbe") !Rpbe - xcpot%is_wc=xcpot%is_name("wc") - - - - END SUBROUTINE xcpot_init +CONTAINS - LOGICAL FUNCTION xcpot_is_name(xcpot,name) - CLASS(t_xcpot),INTENT(IN):: xcpot - CHARACTER(len=*),INTENT(IN) :: name - xcpot_is_name=(trim(xc_names(xcpot%icorr))==trim((name))) - END FUNCTION xcpot_is_name - - LOGICAL FUNCTION xcpot_is_gga(xcpot,icorr) + LOGICAL FUNCTION xcpot_is_gga(xcpot) IMPLICIT NONE CLASS(t_xcpot),INTENT(IN):: xcpot - INTEGER,OPTIONAL,INTENT(IN):: icorr - IF (PRESENT(icorr)) THEN - xcpot_is_gga=priv_gga(icorr) - ELSE - xcpot_is_gga=priv_gga(xcpot%icorr) - ENDIF + xcpot_is_gga=.false. END FUNCTION xcpot_is_gga - LOGICAL FUNCTION xcpot_is_hybrid(xcpot,icorr) + LOGICAL FUNCTION xcpot_is_hybrid(xcpot) IMPLICIT NONE CLASS(t_xcpot),INTENT(IN):: xcpot - INTEGER,OPTIONAL,INTENT(IN):: icorr - IF (PRESENT(icorr)) THEN - xcpot_is_hybrid=priv_hybrid(icorr) - ELSE - xcpot_is_hybrid=priv_hybrid(xcpot%icorr) - ENDIF + xcpot_is_hybrid=.FALSE. END FUNCTION xcpot_is_hybrid - + FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex) USE m_judft IMPLICIT NONE CLASS(t_xcpot),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 - - IF (a_ex==-1) CALL judft_error('xc functional can not be identified') END FUNCTION xcpot_get_exchange_weight + + SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh,vxc,vx,grad) + CLASS(t_xcpot),INTENT(IN) :: xcpot + INTEGER, INTENT (IN) :: jspins + !--> charge density + REAL,INTENT (IN) :: rh(:,:) + !---> xc potential + REAL, INTENT (OUT) :: vxc (:,:),vx(:,:) + TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad + END SUBROUTINE xcpot_get_vxc + + + SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad) + CLASS(t_xcpot),INTENT(IN) :: xcpot + INTEGER, INTENT (IN) :: jspins + !--> charge density + REAL,INTENT (IN) :: rh(:,:) + !---> xc energy density + REAL, INTENT (OUT) :: exc (:) + TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad + END SUBROUTINE xcpot_get_exc + + SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad) + INTEGER, INTENT (IN) :: jspins,ngrid + TYPE(t_gradients),INTENT(OUT):: grad + + !For the in-build xc-pots + ALLOCATE(grad%agrt(ngrid),grad%agru(ngrid),grad%agrd(ngrid)) + ALLOCATE(grad%g2ru(ngrid),grad%g2rd(ngrid),grad%gggrt(ngrid)) + ALLOCATE(grad%gggru(ngrid),grad%gzgr(ngrid),grad%g2rt(ngrid)) + ALLOCATE(grad%gggrd(ngrid),grad%grgru(ngrid),grad%grgrd(ngrid)) + END SUBROUTINE xcpot_alloc_gradients + + SUBROUTINE xcpot_broadcast(xcpot,mpi) + USE m_types_mpi + IMPLICIT NONE + CLASS(t_xcpot),INTENT(INOUT) :: xcpot + TYPE(t_mpi),INTENT(IN) :: mpi + END SUBROUTINE xcpot_broadcast + END MODULE m_types_xcpot diff --git a/types/types_xcpot_data.F90 b/types/types_xcpot_data.F90 new file mode 100644 index 0000000000000000000000000000000000000000..a694a5b803ccc0c2e305ba03f8854974a7e0a905 --- /dev/null +++ b/types/types_xcpot_data.F90 @@ -0,0 +1,26 @@ +!-------------------------------------------------------------------------------- +! 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_types_xcpot_data + !This module contains the xcpot-type used for the in-build xc-implementations + IMPLICIT NONE + + TYPE t_xcpot_data + !in the pbe case (exchpbe.F) lots of test are made + !in addition some constants are set + !to speed up this code precalculate things in init + LOGICAL :: is_rpbe !Rpbe + LOGICAL :: is_wc + LOGICAL :: is_hse !hse,lhse,vhse + REAL :: uk,um + !many logicals to determine xcpot + LOGICAL :: is_pbes !is pbe-sol + LOGICAL :: is_pbe0 + LOGICAL :: is_bh + LOGICAL :: is_mjw + REAL :: exchange_weight + INTEGER :: krla !relativistic corrections + END TYPE t_xcpot_data +END MODULE m_types_xcpot_data diff --git a/types/types_xcpot_inbuild.F90 b/types/types_xcpot_inbuild.F90 new file mode 100644 index 0000000000000000000000000000000000000000..c9be8a88aae494888cdd11d22be9732f8e06ba7b --- /dev/null +++ b/types/types_xcpot_inbuild.F90 @@ -0,0 +1,338 @@ +!-------------------------------------------------------------------------------- +! 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_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_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_gga(20)=[& + .TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,& + .TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,& + .TRUE.,.TRUE.,.TRUE.,.TRUE.,.FALSE.,.FALSE.] + + 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(:) + + CONTAINS + !overloading t_xcpot: + PROCEDURE :: is_gga=>xcpot_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 + 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) + USE m_judft + IMPLICIT NONE + CLASS(t_xcpot_inbuild),INTENT(INOUT) :: xcpot + CHARACTER(len=*),INTENT(IN) :: namex + LOGICAL,INTENT(IN) :: relcor + INTEGER:: n + !Determine icorr from name + 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) + + + !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.123456790123456d0 + ELSE + xcpot%data%um=0.2195149727645171e0 + 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 + + + LOGICAL FUNCTION xcpot_is_gga(xcpot) + IMPLICIT NONE + CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot + xcpot_is_gga=priv_gga(xcpot%icorr) + END FUNCTION xcpot_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 + + IF (a_ex==-1) CALL judft_error('xc functional can not be identified') + END FUNCTION xcpot_get_exchange_weight + +!*********************************************************************** + SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad) +!*********************************************************************** +! + USE m_xcxal, ONLY : vxcxal + USE m_xcwgn, ONLY : vxcwgn + USE m_xcbh, ONLY : vxcbh + USE m_xcvwn, ONLY : vxcvwn + USE m_xcpz, ONLY : vxcpz + USE m_vxcl91 + USE m_vxcwb91 + USE m_vxcpw91 + USE m_vxcepbe + IMPLICIT NONE +!c +!c---> running mode parameters +!c + CLASS(t_xcpot_inbuild),INTENT(IN) :: xcpot + INTEGER, INTENT (IN) :: jspins +!c +!c---> charge density +!c + REAL,INTENT (IN) :: rh(:,:) +!c +!c---> xc potential +!c + REAL, INTENT (OUT) :: vx (:,:) + REAL, INTENT (OUT) :: vxc(:,:) + + ! optional arguments for GGA + TYPE(t_gradients),INTENT(IN),OPTIONAL::grad +!c +!c ---> local scalars + INTEGER :: ngrid + REAL, PARAMETER :: hrtr_half = 0.5 + + !used to be dummy arguments for testing + INTEGER,PARAMETER :: idsprs=0,isprsv=0,iofile=6 + REAL,PARAMETER :: sprsv=0.0 + LOGICAL,PARAMETER :: lwbc=.false. ! l-white-bird-current (ta) +!c +!c.....------------------------------------------------------------------ +!c +!c-----> determine exchange correlation potential +!c + vx (:,:) = 0.0 + vxc(:,:) = 0.0 + ngrid=SIZE(rh,1) + + IF (xcpot%is_gga()) THEN + IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives") + IF (xcpot%is_name("l91")) THEN ! local pw91 + CALL vxcl91(jspins,ngrid,ngrid,rh,grad%agrt,grad%agru,grad%agrd, grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, vx,vxc, isprsv,sprsv) + ELSEIF (xcpot%is_name("pw91")) THEN ! pw91 + IF (lwbc) THEN + CALL vxcwb91(jspins,ngrid,ngrid,rh,grad%agrt,grad%agru,grad%agrd, grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, vx,vxc, idsprs,isprsv,sprsv) + ELSE + + CALL vxcpw91(jspins,ngrid,ngrid,rh,grad%agrt,grad%agru,grad%agrd, grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, vx,vxc, idsprs,isprsv,sprsv) + + ENDIF + ELSE ! pbe or similar + CALL vxcepbe(xcpot%data,jspins,ngrid,ngrid,rh, grad%agrt,grad%agru,grad%agrd,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd, vx,vxc) + ENDIF + ELSE !LDA potentials + IF (xcpot%is_name("x-a")) THEN ! X-alpha method + CALL vxcxal(xcpot%data%krla,jspins, ngrid,ngrid,rh, vx,vxc) + ELSEIF (xcpot%is_name("wign")) THEN ! Wigner interpolation formula + CALL vxcwgn(xcpot%data%krla,jspins, ngrid,ngrid,rh, vx,vxc) + ELSEIF (xcpot%is_name("mjw").OR.xcpot%is_name("bh")) THEN ! von Barth,Hedin correlation + CALL vxcbh(iofile,xcpot%data,jspins, ngrid,ngrid,rh, vx,vxc) + + ELSEIF (xcpot%is_name("vwn")) THEN ! Vosko,Wilk,Nusair correlation + CALL vxcvwn(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, vx,vxc) + ELSEIF (xcpot%is_name("pz")) THEN ! Perdew,Zunger correlation + CALL vxcpz(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, vx,vxc) + ELSEIF (xcpot%is_name("hf")) THEN + ! Hartree-Fock calculation: X-alpha potential is added to generate a rational local potential, + ! later it is subtracted again + CALL vxcxal(xcpot%data%krla,jspins, ngrid,ngrid,rh, vx,vxc) + ! vxc=0 + ELSEIF (xcpot%is_name("exx")) THEN + ! if exact exchange calculation do nothing + vxc = 0 + ELSE + CALL juDFT_error("Unkown LDA potential",calledby="type xcpot") + ENDIF + ENDIF +! +!-----> hartree units +! + vx=hrtr_half*vx + vxc=hrtr_half*vxc + + END SUBROUTINE xcpot_get_vxc + +!*********************************************************************** + SUBROUTINE xcpot_get_exc(xcpot,jspins,rh, exc,grad) +!*********************************************************************** + USE m_xcxal, ONLY : excxal + USE m_xcwgn, ONLY : excwgn + USE m_xcbh, ONLY : excbh + USE m_xcvwn, ONLY : excvwn + USE m_xcpz, ONLY : excpz + USE m_excl91 + USE m_excwb91 + USE m_excpw91 + USE m_excepbe + IMPLICIT NONE +!c +!c---> running mode parameters +!c + CLASS(t_xcpot_inbuild),INTENT(IN) :: xcpot + INTEGER, INTENT (IN) :: jspins +!c +!c---> charge density +!c + REAL,INTENT (IN) :: rh(:,:) +!c +!c---> xc energy density +!c + REAL, INTENT (OUT) :: exc(:) + + ! optional arguments for GGA + TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad +!c +!c ---> local scalars + INTEGER :: ngrid + REAL, PARAMETER :: hrtr_half = 0.5 + + !used to be dummy arguments for testing + INTEGER,PARAMETER :: idsprs=0,isprsv=0,iofile=6 + REAL,PARAMETER :: sprsv=0.0 + LOGICAL,PARAMETER :: lwbc=.false. ! l-white-bird-current (ta) +!c +!c-----> determine exchange correlation energy density +!c + exc(:) = 0.0 + ngrid=SIZE(rh,1) + IF (xcpot%is_gga()) THEN + IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a GGA potential without providing derivatives") + IF (xcpot%is_name("l91")) THEN ! local pw91 + CALL excl91(jspins,SIZE(rh,1),ngrid,rh,grad%agrt,grad%agru,grad%agrd,grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, exc, isprsv,sprsv) + ELSEIF (xcpot%is_name("pw91")) THEN ! pw91 + IF (lwbc) THEN + CALL excwb91(SIZE(grad%agrt),ngrid,rh(:,1),rh(:,2),grad%agrt,grad%agru,grad%agrd, grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, exc, idsprs,isprsv,sprsv) + ELSE + CALL excpw91(jspins,SIZE(grad%agrt),ngrid,rh,grad%agrt,grad%agru,grad%agrd, grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, exc, idsprs,isprsv,sprsv) + ENDIF + ELSE + CALL excepbe(xcpot%data,jspins,SIZE(rh,1),ngrid, rh,grad%agrt,grad%agru,grad%agrd,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd, exc) + ENDIF + ELSE !LDA + IF (xcpot%is_name("x-a")) THEN ! X-alpha method + CALL excxal(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, exc) + ELSEIF (xcpot%is_name("wign")) THEN ! Wigner interpolation formula + CALL excwgn(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, exc) + ELSEIF (xcpot%is_name("mjw").OR.xcpot%is_name("bh")) THEN ! von Barth,Hedin correlation + CALL excbh(iofile,xcpot%data,jspins, ngrid,ngrid,rh, exc) + ELSEIF (xcpot%is_name("vwn")) THEN ! Vosko,Wilk,Nusair correlation + CALL excvwn(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, exc) + ELSEIF (xcpot%is_name("pz")) THEN ! Perdew,Zunger correlation + CALL excpz(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, exc) + ELSEIF (xcpot%is_name("hf") .OR. xcpot%is_name("exx")) THEN + exc=0 + ELSE + CALL juDFT_error("Unkown LDA potential",calledby="type xcpot") + ENDIF + ENDIF +!c-----> hartree units + exc= hrtr_half*exc + 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 + + SUBROUTINE xcpot_broadcast(xcpot,mpi) + USE m_types_mpi + IMPLICIT NONE + CLASS(t_xcpot),INTENT(INOUT) :: xcpot + TYPE(t_mpi),INTENT(IN) :: mpi + + LOGICAL :: l_relcor + CHARACTER(len=4):: namex + INTEGER :: ierr +#ifdef CPP_MPI + INCLUDE 'mpif.h' + + IF (mpi%irank==0) THEN + namex=xcpot%get_name() + l_relcor=xcpot%krla==1 + ENDIF + CALL MPI_BCAST(namex,4,MPI_CHARACTER,0,mpi%mpi_comm,ierr) + CALL MPI_BCAST(l_relcor,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr) + IF (mpi%irank.NE.0) CALL xcpot%init(namex,l_relcor) +#endif + END SUBROUTINE xcpot_broadcast + + + END MODULE m_types_xcpot_inbuild diff --git a/vgen/lhglptg.f90 b/vgen/lhglptg.f90 index 22d92dcb4efc4df6b4a241b65abbbe7b7ee1f3b0..9acc3ba6deda2fc88a32097da9861d7e8078aaf2 100644 --- a/vgen/lhglptg.f90 +++ b/vgen/lhglptg.f90 @@ -16,7 +16,7 @@ CONTAINS USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_sym),INTENT(IN) :: sym TYPE(t_sphhar),INTENT(IN) :: sphhar TYPE(t_atoms),INTENT(IN) :: atoms diff --git a/vgen/mkgxyz3.f b/vgen/mkgxyz3.f index 40aaf8a04b588f3ba4f96d31def8833b0b9b38ea..a0087766adc33b7b6898cdbf25ea591266ce48ca 100644 --- a/vgen/mkgxyz3.f +++ b/vgen/mkgxyz3.f @@ -1,3 +1,8 @@ +!-------------------------------------------------------------------------------- +! 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_mkgxyz3 c.....------------------------------------------------------------------ c by use of cartesian x,y,z components of charge density gradients, @@ -10,9 +15,8 @@ c.....------------------------------------------------------------------ SUBROUTINE mkgxyz3( > ndm,jsdm,ng3,jspins,vl, > dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy, - < agrt,agru,agrd, g2rt,g2ru,g2rd, - < gggrt,gggru,gggrd,gzgr) - + < grad) + USE m_types IMPLICIT NONE INTEGER, INTENT (IN) :: ndm,ng3,jsdm,jspins REAL, INTENT (IN) :: vl(ndm,jsdm) @@ -20,9 +24,7 @@ c.....------------------------------------------------------------------ REAL, INTENT (IN) :: dvxx(ndm,jsdm),dvyy(ndm,jsdm),dvzz(ndm,jsdm) REAL, INTENT (IN) :: dvyz(ndm,jsdm),dvzx(ndm,jsdm),dvxy(ndm,jsdm) - REAL, INTENT (OUT) :: agrt(ndm),agru(ndm),agrd(ndm) - REAL, INTENT (OUT) :: g2rt(ndm),g2ru(ndm),g2rd(ndm) - REAL, INTENT (OUT) :: gggrt(ndm),gggru(ndm),gggrd(ndm),gzgr(ndm) + TYPE(t_gradients),INTENT(OUT)::grad REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvzxt,dvxyt, & vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvzxu,dvxyu, @@ -34,17 +36,37 @@ c.....------------------------------------------------------------------ sml = 1.e-14 + IF(ALLOCATED(grad%sigma)) THEN +!Use only contracted gradients for libxc + if (jspins==1) THEN + DO i=1,ng3 + grad%sigma(1,i)= + + dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1) + ENDDO + ELSE + DO i=1,ng3 + grad%sigma(1,i)= + + dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1) + grad%sigma(2,i)= + + dvx(i,1)*dvx(i,2)+dvy(i,1)*dvy(i,2)+dvz(i,1)*dvz(i,2) + grad%sigma(3,i)= + + dvx(i,2)*dvx(i,2)+dvy(i,2)*dvy(i,2)+dvz(i,2)*dvz(i,2) + ENDDO + ENDIF + RETURN + ENDIF + DO i = 1,ndm - agrt(i) = 0.0 - agru(i) = 0.0 - agrd(i) = 0.0 - gggrt(i) = 0.0 - gggru(i) = 0.0 - gggrd(i) = 0.0 - gzgr(i) = 0.0 - g2rt(i) = 0.0 - g2ru(i) = 0.0 - g2rd(i) = 0.0 + grad%agrt(i) = 0.0 + grad%agru(i) = 0.0 + grad%agrd(i) = 0.0 + grad%gggrt(i) = 0.0 + grad%gggru(i) = 0.0 + grad%gggrd(i) = 0.0 + grad%gzgr(i) = 0.0 + grad%g2rt(i) = 0.0 + grad%g2ru(i) = 0.0 + grad%g2rd(i) = 0.0 ENDDO IF (jspins.eq.1) THEN @@ -88,25 +110,25 @@ c.....------------------------------------------------------------------ c agr: abs(grad(ro)), t,u,d for total, up and down. - agrt(i) = max(sqrt(dvxt**2+dvyt**2+dvzt**2),sml) - agru(i) = max(sqrt(dvxu**2+dvyu**2+dvzu**2),sml) - agrd(i) = max(sqrt(dvxd**2+dvyd**2+dvzd**2),sml) + grad%agrt(i) = max(sqrt(dvxt**2+dvyt**2+dvzt**2),sml) + grad%agru(i) = max(sqrt(dvxu**2+dvyu**2+dvzu**2),sml) + grad%agrd(i) = max(sqrt(dvxd**2+dvyd**2+dvzd**2),sml) - dagrxt = (dvxt*dvxxt+dvyt*dvxyt+dvzt*dvzxt)/agrt(i) - dagrxu = (dvxu*dvxxu+dvyu*dvxyu+dvzu*dvzxu)/agru(i) - dagrxd = (dvxd*dvxxd+dvyd*dvxyd+dvzd*dvzxd)/agrd(i) + dagrxt = (dvxt*dvxxt+dvyt*dvxyt+dvzt*dvzxt)/grad%agrt(i) + dagrxu = (dvxu*dvxxu+dvyu*dvxyu+dvzu*dvzxu)/grad%agru(i) + dagrxd = (dvxd*dvxxd+dvyd*dvxyd+dvzd*dvzxd)/grad%agrd(i) - dagryt = (dvxt*dvxyt+dvyt*dvyyt+dvzt*dvyzt)/agrt(i) - dagryu = (dvxu*dvxyu+dvyu*dvyyu+dvzu*dvyzu)/agru(i) - dagryd = (dvxd*dvxyd+dvyd*dvyyd+dvzd*dvyzd)/agrd(i) + dagryt = (dvxt*dvxyt+dvyt*dvyyt+dvzt*dvyzt)/grad%agrt(i) + dagryu = (dvxu*dvxyu+dvyu*dvyyu+dvzu*dvyzu)/grad%agru(i) + dagryd = (dvxd*dvxyd+dvyd*dvyyd+dvzd*dvyzd)/grad%agrd(i) - dagrzt = (dvxt*dvzxt+dvyt*dvyzt+dvzt*dvzzt)/agrt(i) - dagrzu = (dvxu*dvzxu+dvyu*dvyzu+dvzu*dvzzu)/agru(i) - dagrzd = (dvxd*dvzxd+dvyd*dvyzd+dvzd*dvzzd)/agrd(i) + dagrzt = (dvxt*dvzxt+dvyt*dvyzt+dvzt*dvzzt)/grad%agrt(i) + dagrzu = (dvxu*dvzxu+dvyu*dvyzu+dvzu*dvzzu)/grad%agru(i) + dagrzd = (dvxd*dvzxd+dvyd*dvyzd+dvzd*dvzzd)/grad%agrd(i) - gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt - gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu - gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd + grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt + grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu + grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd c dzdx=d(zeta)/dx,.. @@ -116,13 +138,13 @@ c dzdx=d(zeta)/dx,.. c gzgr=grad(zeta)*grad(ro). - gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt + grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt c g2r: grad(grad(ro)) - g2rt(i) = dvxxt + dvyyt + dvzzt - g2ru(i) = dvxxu + dvyyu + dvzzu - g2rd(i) = dvxxd + dvyyd + dvzzd + grad%g2rt(i) = dvxxt + dvyyt + dvzzt + grad%g2ru(i) = dvxxu + dvyyu + dvzzu + grad%g2rd(i) = dvxxd + dvyyd + dvzzd 10 ENDDO @@ -167,25 +189,25 @@ c g2r: grad(grad(ro)) c agr: abs(grad(ro)), t,u,d for total, up and down. - agrt(i) = max(sqrt(dvxt**2+dvyt**2+dvzt**2),sml) - agru(i) = max(sqrt(dvxu**2+dvyu**2+dvzu**2),sml) - agrd(i) = max(sqrt(dvxd**2+dvyd**2+dvzd**2),sml) + grad%agrt(i) = max(sqrt(dvxt**2+dvyt**2+dvzt**2),sml) + grad%agru(i) = max(sqrt(dvxu**2+dvyu**2+dvzu**2),sml) + grad%agrd(i) = max(sqrt(dvxd**2+dvyd**2+dvzd**2),sml) - dagrxt = (dvxt*dvxxt+dvyt*dvxyt+dvzt*dvzxt)/agrt(i) - dagrxu = (dvxu*dvxxu+dvyu*dvxyu+dvzu*dvzxu)/agru(i) - dagrxd = (dvxd*dvxxd+dvyd*dvxyd+dvzd*dvzxd)/agrd(i) + dagrxt = (dvxt*dvxxt+dvyt*dvxyt+dvzt*dvzxt)/grad%agrt(i) + dagrxu = (dvxu*dvxxu+dvyu*dvxyu+dvzu*dvzxu)/grad%agru(i) + dagrxd = (dvxd*dvxxd+dvyd*dvxyd+dvzd*dvzxd)/grad%agrd(i) - dagryt = (dvxt*dvxyt+dvyt*dvyyt+dvzt*dvyzt)/agrt(i) - dagryu = (dvxu*dvxyu+dvyu*dvyyu+dvzu*dvyzu)/agru(i) - dagryd = (dvxd*dvxyd+dvyd*dvyyd+dvzd*dvyzd)/agrd(i) + dagryt = (dvxt*dvxyt+dvyt*dvyyt+dvzt*dvyzt)/grad%agrt(i) + dagryu = (dvxu*dvxyu+dvyu*dvyyu+dvzu*dvyzu)/grad%agru(i) + dagryd = (dvxd*dvxyd+dvyd*dvyyd+dvzd*dvyzd)/grad%agrd(i) - dagrzt = (dvxt*dvzxt+dvyt*dvyzt+dvzt*dvzzt)/agrt(i) - dagrzu = (dvxu*dvzxu+dvyu*dvyzu+dvzu*dvzzu)/agru(i) - dagrzd = (dvxd*dvzxd+dvyd*dvyzd+dvzd*dvzzd)/agrd(i) + dagrzt = (dvxt*dvzxt+dvyt*dvyzt+dvzt*dvzzt)/grad%agrt(i) + dagrzu = (dvxu*dvzxu+dvyu*dvyzu+dvzu*dvzzu)/grad%agru(i) + dagrzd = (dvxd*dvzxd+dvyd*dvyzd+dvzd*dvzzd)/grad%agrd(i) - gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt - gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu - gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd + grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt + grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu + grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd c dzdx=d(zeta)/dx,.. @@ -195,13 +217,13 @@ c dzdx=d(zeta)/dx,.. c gzgr=grad(zeta)*grad(ro). - gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt + grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt c g2r: grad(grad(ro)) - g2rt(i) = dvxxt + dvyyt + dvzzt - g2ru(i) = dvxxu + dvyyu + dvzzu - g2rd(i) = dvxxd + dvyyd + dvzzd + grad%g2rt(i) = dvxxt + dvyyt + dvzzt + grad%g2ru(i) = dvxxu + dvyyu + dvzzu + grad%g2rd(i) = dvxxd + dvyyd + dvzzd 20 ENDDO diff --git a/vgen/mkgylm.f b/vgen/mkgylm.f index 901bdd860e69172c21ea18e106170635497375c5..bed6a2affe295a643d4b3b4c9727513224adee8a 100644 --- a/vgen/mkgylm.f +++ b/vgen/mkgylm.f @@ -1,31 +1,35 @@ +!-------------------------------------------------------------------------------- +! 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_mkgylm CONTAINS SUBROUTINE mkgylm( - > jspins,rv,thet,nsp,nspd,jspd, - > rh,rhdr,rhdt,rhdf,rhdrr,rhdtt,rhdff, - > rhdtf,rhdrt,rhdrf, - < agr,agru,agrd,g2r,g2ru,g2rd,gggr, - < gggru,gggrd,grgru,grgrd,gzgr) + > jspins,rv,thet,nsp,nspd,jspd, + > rh,rhdr,rhdt,rhdf,rhdrr,rhdtt,rhdff, + > rhdtf,rhdrt,rhdrf,grad) c.....------------------------------------------------------------------ c by use of charge density and its polar coord. gradient components -cc calculate agr and others used to evaluate gradient -cc contributions to potential and energy. t.a. 1996. +c c calculate agr and others used to evaluate gradient +c c contributions to potential and energy. t.a. 1996. c.....------------------------------------------------------------------ c ro=sum(ro*ylh), rdr=sum(drr*ylh), drdr=sum(ddrr*ylh), -cc rdt=sum(ro*ylht1), rdtt=sum(ro*ylht2), ... -cc rdf=sum(ro*ylhf1), rdff=sum(ro*ylhf2), ... -cc rdtf=sum(ro*ylhtf), rdff=sum(ro*ylhf2), ... -cc rdrt=sum(drr*ylht1),rdrf=sum(drr*ylhf1), +c c rdt=sum(ro*ylht1), rdtt=sum(ro*ylht2), ... +c c rdf=sum(ro*ylhf1), rdff=sum(ro*ylhf2), ... +c c rdtf=sum(ro*ylhtf), rdff=sum(ro*ylhf2), ... +c c rdrt=sum(drr*ylht1),rdrf=sum(drr*ylhf1), c agr: abs(grad(ro)), g2r: laplacian(ro), -cc gggr: grad(ro)*grad(agr), -cc grgru,d: grad(ro)*grad(rou),for rod., gzgr: grad(zeta)*grad(ro). +c c gggr: grad(ro)*grad(agr), +c c grgru,d: grad(ro)*grad(rou),for rod., gzgr: grad(zeta)*grad(ro). c dagrr,-t,-f: d(agr)/dr, d(agr)/dth/r, d(agr)/dfi/r/sint. c.....------------------------------------------------------------------ + USE m_types IMPLICIT NONE -C .. -C .. Arguments .. +C .. +C .. Arguments .. INTEGER, INTENT (IN) :: nspd,jspd REAL, INTENT (IN) :: rv REAL, INTENT (IN) :: thet(nspd) @@ -34,19 +38,17 @@ C .. Arguments .. REAL, INTENT (IN) :: rhdtt(nspd,jspd),rhdff(nspd,jspd) REAL, INTENT (IN) :: rhdtf(nspd,jspd),rhdrt(nspd,jspd) REAL, INTENT (IN) :: rhdrf(nspd,jspd),rhdt(nspd,jspd) - REAL, INTENT (OUT) :: agr(nspd),agru(nspd),agrd(nspd) - REAL, INTENT (OUT) :: g2ru(nspd),g2rd(nspd),gggr(nspd) - REAL, INTENT (OUT) :: gggru(nspd),gzgr(nspd),g2r(nspd) - REAL, INTENT (OUT) :: gggrd(nspd),grgru(nspd),grgrd(nspd) -C .. -C .. Locals .. + TYPE(t_gradients),INTENT(OUT) ::grad + +C .. +C .. Locals .. INTEGER i,jspins,nsp REAL chsml, dagrf,dagrfd,dagrfu,dagrr,dagrrd,dagrru,dagrt, - + dagrtd,dagrtu,drdr,dzdfs,dzdr,dzdtr,grf,grfd,grfu,grr, - + grrd,grru,grt,grtd,grtu,rdf,rdfd,rdff,rdffd,rdffu,rdfu, - + rdr,rdrd,rdrf,rdrfd,rdrfu,rdrrd,rdrru,rdrt,rdrtd,rdrtu, - + rdru,rdt,rdtd,rdtf,rdtfd,rdtfu,rdtt,rdttd,rdttu,rdtu, - + ro,ro2,rod,rou,rv1,rv2,rv3,rvsin1,sint1,sint2,tant1 + + dagrtd,dagrtu,drdr,dzdfs,dzdr,dzdtr,grf,grfd,grfu,grr, + + grrd,grru,grt,grtd,grtu,rdf,rdfd,rdff,rdffd,rdffu,rdfu, + + rdr,rdrd,rdrf,rdrfd,rdrfu,rdrrd,rdrru,rdrt,rdrtd,rdrtu, + + rdru,rdt,rdtd,rdtf,rdtfd,rdtfu,rdtt,rdttd,rdttu,rdtu, + + ro,ro2,rod,rou,rv1,rv2,rv3,rvsin1,sint1,sint2,tant1 chsml = 1.e-10 @@ -54,286 +56,310 @@ C .. Locals .. rv2 = rv1**2 rv3 = rv1**3 - IF (jspins.EQ.1) THEN - - points_1 : DO i = 1,nsp - - agr(i) = 0.0 - agru(i) = 0.0 - agrd(i) = 0.0 - g2r(i) = 0.0 - g2ru(i) = 0.0 - g2rd(i) = 0.0 - gggr(i) = 0.0 - gggru(i) = 0.0 - gggrd(i) = 0.0 - grgru(i) = 0.0 - grgrd(i) = 0.0 - gzgr(i) = 0.0 - - ro = rh(i,1) - - IF (ro.LT.chsml) CYCLE points_1 - - sint1 = sin(thet(i)) - sint2 = sint1**2 - tant1 = tan(thet(i)) - rvsin1 = rv1*sint1 - - - rou = ro/2 - rdru = rhdr(i,1)/2 - rdtu = rhdt(i,1)/2 - rdfu = rhdf(i,1)/2 - rdrru = rhdrr(i,1)/2 - rdttu = rhdtt(i,1)/2 - rdffu = rhdff(i,1)/2 - rdtfu = rhdtf(i,1)/2 - rdrtu = rhdrt(i,1)/2 - rdrfu = rhdrf(i,1)/2 - - rod = rou - rdrd = rdru - rdtd = rdtu - rdfd = rdfu - rdrrd = rdrru - rdttd = rdttu - rdffd = rdffu - rdtfd = rdtfu - rdrtd = rdrtu - rdrfd = rdrfu - - rdr = rdru + rdrd - rdt = rdtu + rdtd - rdf = rdfu + rdfd - drdr = rdrru + rdrrd - rdtt = rdttu + rdttd - rdff = rdffu + rdffd - rdtf = rdtfu + rdtfd - rdrt = rdrtu + rdrtd - rdrf = rdrfu + rdrfd - - ro2 = ro**2 - - grr = rdr - grt = rdt/rv1 - grf = rdf/rvsin1 - - agr(i) = sqrt(grr**2+grt**2+grf**2) - - IF (agr(i).LT.chsml) CYCLE points_1 - - dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ - + rdf* (rdrf*rv1-rdf)/sint2)/agr(i)/rv3 - - dagrt = (rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/ - + (agr(i)*rv3) - - dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/ - + (agr(i)*rv3*sint1) - - g2r(i) = drdr + 2.e0*rdr/rv1 + (rdtt+rdt/tant1+rdff/sint2)/rv2 - - - dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2 - -!--> dzdtr,dzdfs vanish by definition. - - dzdtr = 0.0 - dzdfs = 0.0 - - gggr(i) = grr*dagrr + grt*dagrt + grf*dagrf - - gzgr(i) = dzdr*grr + dzdtr*grt + dzdfs*grf - - grru = rdru - grtu = rdtu/rv1 - grfu = rdfu/rvsin1 - - agru(i) = sqrt(grru**2+grtu**2+grfu**2) - - dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+ - + rdfu* (rdrfu*rv1-rdfu)/sint2)/agru(i)/rv3 - - dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+ - + rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (agru(i)*rv3) - - dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/ - + (agru(i)*rv3*sint1) - - g2ru(i) = rdrru + 2.e0*rdru/rv1 + - + (rdttu+rdtu/tant1+rdffu/sint2)/rv2 - - gggru(i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu - - grgru(i) = grr*grru + grt*grtu + grf*grfu - - grrd = rdrd - grtd = rdtd/rv1 - grfd = rdfd/rvsin1 - - agrd(i) = sqrt(grrd**2+grtd**2+grfd**2) - - dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+ - + rdfd* (rdrfd*rv1-rdfd)/sint2)/agrd(i)/rv3 - - dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+ - + rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (agrd(i)*rv3) - - dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/ - + (agrd(i)*rv3*sint1) - - g2rd(i) = rdrrd + 2*rdrd/rv1 + - + (rdttd+rdtd/tant1+rdffd/sint2)/rv2 - - gggrd(i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd - - grgrd(i) = grr*grrd + grt*grtd + grf*grfd - - - ENDDO points_1 - + IF (ALLOCATED(grad%sigma)) THEN +! Contracted gradients for libxc are needed only + IF (jspins==1) THEN + DO i=1,nsp + rdru = rhdr(i,1) + rdtu = rhdt(i,1) + rdfu = rhdf(i,1) + grad%sigma(1,i)=rdru*rdru+rdtu*rdtu+rdfu*rdfu + ENDDO + ELSE + DO i=1,nsp + rdru = rhdr(i,1) + rdtu = rhdt(i,1) + rdfu = rhdf(i,1) + rdrd = rhdr(i,2) + rdtd = rhdt(i,2) + rdfd = rhdf(i,2) + grad%sigma(1,i)=rdru*rdru+rdtu*rdtu+rdfu*rdfu + grad%sigma(2,i)=rdru*rdrd+rdtu*rdtd+rdfu*rdfd + grad%sigma(3,i)=rdrd*rdrd+rdtd*rdtd+rdfd*rdfd + ENDDO + ENDIF ELSE - - points_2 : DO i = 1,nsp - - agr(i) = 0.0 - agru(i) = 0.0 - agrd(i) = 0.0 - g2r(i) = 0.0 - g2ru(i) = 0.0 - g2rd(i) = 0.0 - gggr(i) = 0.0 - gggru(i) = 0.0 - gggrd(i) = 0.0 - grgru(i) = 0.0 - grgrd(i) = 0.0 - gzgr(i) = 0.0 - - ro = rh(i,1) + rh(i,jspins) - - IF (ro.LT.chsml) CYCLE points_2 - - sint1 = sin(thet(i)) - sint2 = sint1**2 - tant1 = tan(thet(i)) - rvsin1 = rv1*sint1 - - rou = rh(i,1) - rdru = rhdr(i,1) - rdtu = rhdt(i,1) - rdfu = rhdf(i,1) - rdrru = rhdrr(i,1) - rdttu = rhdtt(i,1) - rdffu = rhdff(i,1) - rdtfu = rhdtf(i,1) - rdrtu = rhdrt(i,1) - rdrfu = rhdrf(i,1) - - rod = rh(i,jspins) - rdrd = rhdr(i,jspins) - rdtd = rhdt(i,jspins) - rdfd = rhdf(i,jspins) - rdrrd = rhdrr(i,jspins) - rdttd = rhdtt(i,jspins) - rdffd = rhdff(i,jspins) - rdtfd = rhdtf(i,jspins) - rdrtd = rhdrt(i,jspins) - rdrfd = rhdrf(i,jspins) - - rdr = rdru + rdrd - rdt = rdtu + rdtd - rdf = rdfu + rdfd - drdr = rdrru + rdrrd - rdtt = rdttu + rdttd - rdff = rdffu + rdffd - rdtf = rdtfu + rdtfd - rdrt = rdrtu + rdrtd - rdrf = rdrfu + rdrfd - - ro2 = ro**2 - - grr = rdr - grt = rdt/rv1 - grf = rdf/rvsin1 - - agr(i) = sqrt(grr**2+grt**2+grf**2) - - IF (agr(i).LT.chsml) CYCLE points_2 - - dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ - + rdf* (rdrf*rv1-rdf)/sint2)/agr(i)/rv3 - - dagrt = (rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/ - + (agr(i)*rv3) - - dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/ - + (agr(i)*rv3*sint1) - - g2r(i) = drdr + 2.e0*rdr/rv1 + (rdtt+rdt/tant1+rdff/sint2)/rv2 - - dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2 - -c dzdtr,dzdfs vanish by definition. - dzdtr = 0.0 - dzdfs = 0.0 - - gggr(i) = grr*dagrr + grt*dagrt + grf*dagrf - - gzgr(i) = dzdr*grr + dzdtr*grt + dzdfs*grf - - grru = rdru - grtu = rdtu/rv1 - grfu = rdfu/rvsin1 - - agru(i) = sqrt(grru**2+grtu**2+grfu**2) - - dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+ - + rdfu* (rdrfu*rv1-rdfu)/sint2)/agru(i)/rv3 - - dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+ - + rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (agru(i)*rv3) - - dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/ - + (agru(i)*rv3*sint1) - - g2ru(i) = rdrru + 2.e0*rdru/rv1 + - + (rdttu+rdtu/tant1+rdffu/sint2)/rv2 - - gggru(i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu - - grgru(i) = grr*grru + grt*grtu + grf*grfu - - - grrd = rdrd - grtd = rdtd/rv1 - grfd = rdfd/rvsin1 - - agrd(i) = sqrt(grrd**2+grtd**2+grfd**2) - - dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+ - + rdfd* (rdrfd*rv1-rdfd)/sint2)/agrd(i)/rv3 - - dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+ - + rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (agrd(i)*rv3) - - dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/ - + (agrd(i)*rv3*sint1) - - - - g2rd(i) = rdrrd + 2*rdrd/rv1 + - + (rdttd+rdtd/tant1+rdffd/sint2)/rv2 - - gggrd(i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd - - grgrd(i) = grr*grrd + grt*grtd + grf*grfd - - - ENDDO points_2 - +! Old code for in-build xcpots + IF (jspins.EQ.1) THEN + + points_1 : DO i = 1,nsp + + grad%agrt(i) = 0.0 + grad%agru(i) = 0.0 + grad%agrd(i) = 0.0 + grad%g2rt(i) = 0.0 + grad%g2ru(i) = 0.0 + grad%g2rd(i) = 0.0 + grad%gggrt(i) = 0.0 + grad%gggru(i) = 0.0 + grad%gggrd(i) = 0.0 + grad%grgru(i) = 0.0 + grad%grgrd(i) = 0.0 + grad%gzgr(i) = 0.0 + + ro = rh(i,1) + + IF (ro.LT.chsml) CYCLE points_1 + sint1 = sin(thet(i)) + sint2 = sint1**2 + tant1 = tan(thet(i)) + rvsin1 = rv1*sint1 + + + rou = ro/2 + rdru = rhdr(i,1)/2 + rdtu = rhdt(i,1)/2 + rdfu = rhdf(i,1)/2 + rdrru = rhdrr(i,1)/2 + rdttu = rhdtt(i,1)/2 + rdffu = rhdff(i,1)/2 + rdtfu = rhdtf(i,1)/2 + rdrtu = rhdrt(i,1)/2 + rdrfu = rhdrf(i,1)/2 + + rod = rou + rdrd = rdru + rdtd = rdtu + rdfd = rdfu + rdrrd = rdrru + rdttd = rdttu + rdffd = rdffu + rdtfd = rdtfu + rdrtd = rdrtu + rdrfd = rdrfu + + rdr = rdru + rdrd + rdt = rdtu + rdtd + rdf = rdfu + rdfd + drdr = rdrru + rdrrd + rdtt = rdttu + rdttd + rdff = rdffu + rdffd + rdtf = rdtfu + rdtfd + rdrt = rdrtu + rdrtd + rdrf = rdrfu + rdrfd + + ro2 = ro**2 + + grr = rdr + grt = rdt/rv1 + grf = rdf/rvsin1 + + grad%agrt(i) = sqrt(grr**2+grt**2+grf**2) + + IF (grad%agrt(i).LT.chsml) CYCLE points_1 + + dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ + + rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(i)/rv3 + + dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/ + + (grad%agrt(i)*rv3) + + dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/ + + (grad%agrt(i)*rv3*sint1) + + grad%g2rt(i)=drdr+2.0*rdr/rv1+ + + (rdtt+rdt/tant1+rdff/sint2)/rv2 + + + dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2 + +!-- > dzdtr,dzdfs vanish by definition. + + dzdtr = 0.0 + dzdfs = 0.0 + + grad%gggrt(i) = grr*dagrr + grt*dagrt + grf*dagrf + + grad%gzgr(i) = dzdr*grr + dzdtr*grt + dzdfs*grf + + grru = rdru + grtu = rdtu/rv1 + grfu = rdfu/rvsin1 + + grad%agru(i) = sqrt(grru**2+grtu**2+grfu**2) + + dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+ + + rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(i)/rv3 + + dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+ + + rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(i)*rv3) + + dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/ + + (grad%agru(i)*rv3*sint1) + + grad%g2ru(i) = rdrru + 2.e0*rdru/rv1 + + + (rdttu+rdtu/tant1+rdffu/sint2)/rv2 + + grad%gggru(i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu + + grad%grgru(i) = grr*grru + grt*grtu + grf*grfu + + grrd = rdrd + grtd = rdtd/rv1 + grfd = rdfd/rvsin1 + + grad%agrd(i) = sqrt(grrd**2+grtd**2+grfd**2) + + dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+ + + rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(i)/rv3 + + dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+ + + rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(i)*rv3) + + dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/ + + (grad%agrd(i)*rv3*sint1) + + grad%g2rd(i) = rdrrd + 2*rdrd/rv1 + + + (rdttd+rdtd/tant1+rdffd/sint2)/rv2 + + grad%gggrd(i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd + + grad%grgrd(i) = grr*grrd + grt*grtd + grf*grfd + + + ENDDO points_1 + + ELSE + + points_2 : DO i = 1,nsp + + grad%agrt(i) = 0.0 + grad%agru(i) = 0.0 + grad%agrd(i) = 0.0 + grad%g2rt(i) = 0.0 + grad%g2ru(i) = 0.0 + grad%g2rd(i) = 0.0 + grad%gggrt(i) = 0.0 + grad%gggru(i) = 0.0 + grad%gggrd(i) = 0.0 + grad%grgru(i) = 0.0 + grad%grgrd(i) = 0.0 + grad%gzgr(i) = 0.0 + + ro = rh(i,1) + rh(i,jspins) + + IF (ro.LT.chsml) CYCLE points_2 + + sint1 = sin(thet(i)) + sint2 = sint1**2 + tant1 = tan(thet(i)) + rvsin1 = rv1*sint1 + + rou = rh(i,1) + rdru = rhdr(i,1) + rdtu = rhdt(i,1) + rdfu = rhdf(i,1) + rdrru = rhdrr(i,1) + rdttu = rhdtt(i,1) + rdffu = rhdff(i,1) + rdtfu = rhdtf(i,1) + rdrtu = rhdrt(i,1) + rdrfu = rhdrf(i,1) + + rod = rh(i,jspins) + rdrd = rhdr(i,jspins) + rdtd = rhdt(i,jspins) + rdfd = rhdf(i,jspins) + rdrrd = rhdrr(i,jspins) + rdttd = rhdtt(i,jspins) + rdffd = rhdff(i,jspins) + rdtfd = rhdtf(i,jspins) + rdrtd = rhdrt(i,jspins) + rdrfd = rhdrf(i,jspins) + + rdr = rdru + rdrd + rdt = rdtu + rdtd + rdf = rdfu + rdfd + drdr = rdrru + rdrrd + rdtt = rdttu + rdttd + rdff = rdffu + rdffd + rdtf = rdtfu + rdtfd + rdrt = rdrtu + rdrtd + rdrf = rdrfu + rdrfd + + ro2 = ro**2 + + grr = rdr + grt = rdt/rv1 + grf = rdf/rvsin1 + + grad%agrt(i) = sqrt(grr**2+grt**2+grf**2) + + IF (grad%agrt(i).LT.chsml) CYCLE points_2 + + dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ + + rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(i)/rv3 + + dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/ + + (grad%agrt(i)*rv3) + + dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/ + + (grad%agrt(i)*rv3*sint1) + + grad%g2rt(i)=drdr+2.0*rdr/rv1+(rdtt+rdt/tant1+rdff/sint2)/rv2 + + dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2 + +c dzdtr,dzdfs vanish by definition. + dzdtr = 0.0 + dzdfs = 0.0 + + grad%gggrt(i) = grr*dagrr + grt*dagrt + grf*dagrf + + grad%gzgr(i) = dzdr*grr + dzdtr*grt + dzdfs*grf + + grru = rdru + grtu = rdtu/rv1 + grfu = rdfu/rvsin1 + + grad%agru(i) = sqrt(grru**2+grtu**2+grfu**2) + + dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+ + + rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(i)/rv3 + + dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+ + + rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(i)*rv3) + + dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/ + + (grad%agru(i)*rv3*sint1) + + grad%g2ru(i) = rdrru + 2.e0*rdru/rv1 + + + (rdttu+rdtu/tant1+rdffu/sint2)/rv2 + + grad%gggru(i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu + + grad%grgru(i) = grr*grru + grt*grtu + grf*grfu + + + grrd = rdrd + grtd = rdtd/rv1 + grfd = rdfd/rvsin1 + + grad%agrd(i) = sqrt(grrd**2+grtd**2+grfd**2) + + dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+ + + rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(i)/rv3 + + dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+ + + rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(i)*rv3) + + dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/ + + (grad%agrd(i)*rv3*sint1) + + + + grad%g2rd(i) = rdrrd + 2*rdrd/rv1 + + + (rdttd+rdtd/tant1+rdffd/sint2)/rv2 + + grad%gggrd(i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd + + grad%grgrd(i) = grr*grrd + grt*grtd + grf*grfd + + + ENDDO points_2 + + ENDIF ENDIF - RETURN END SUBROUTINE mkgylm END MODULE m_mkgylm diff --git a/vgen/mkgz.f b/vgen/mkgz.f index dc3fbf130c52399ba569c506c84247e34fefdbf5..5a94a0011787cd979c27973d0929e6090c4e6438 100644 --- a/vgen/mkgz.f +++ b/vgen/mkgz.f @@ -9,17 +9,14 @@ cc energy. c.....------------------------------------------------------------------ CONTAINS SUBROUTINE mkgz( - > nmzdf,jspins,rh1,rh2,rhdz1,rhdz2,rhdzz1,rhdzz2, - < agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr) - + > nmzdf,jspins,rh1,rh2,rhdz1,rhdz2,rhdzz1,rhdzz2, + < grad) + USE m_types IMPLICIT NONE INTEGER, INTENT (IN) :: nmzdf,jspins REAL, INTENT (IN) :: rh1(nmzdf),rhdz1(nmzdf),rhdzz1(nmzdf) REAL, INTENT (IN) :: rh2(nmzdf),rhdz2(nmzdf),rhdzz2(nmzdf) - REAL, INTENT (OUT) :: agrt(nmzdf),agru(nmzdf),agrd(nmzdf) - REAL, INTENT (OUT) :: g2rt(nmzdf),g2ru(nmzdf),g2rd(nmzdf) - REAL, INTENT (OUT) :: gggrt(nmzdf),gggru(nmzdf),gggrd(nmzdf) - REAL, INTENT (OUT) :: gzgr(nmzdf) + TYPE(t_gradients),INTENT(OUT)::grad INTEGER i REAL vlt,dvzt,dvzzt,vlu,dvzu,dvzzu,vld,dvzd,dvzzd @@ -27,6 +24,21 @@ c.....------------------------------------------------------------------ sml = 1.e-14 + IF (ALLOCATED(grad%sigma)) THEN + IF(jspins==1) THEN + DO i=1,nmzdf + grad%sigma(1,i)=rhdz1(i)*rhdz1(i) + ENDDO + ELSE + DO i=1,nmzdf + grad%sigma(1,i)=rhdz1(i)*rhdz1(i) + grad%sigma(2,i)=rhdz1(i)*rhdz2(i) + grad%sigma(3,i)=rhdz2(i)*rhdz2(i) + ENDDO + ENDIF + RETURN + ENDIF + IF (jspins == 1) THEN DO i = 1,nmzdf @@ -44,17 +56,17 @@ c.....------------------------------------------------------------------ c agrt: abs(grad(ro)), u,d for up and down. - agrt(i) = max(abs(dvzt),sml) - agru(i) = max(abs(dvzu),sml) - agrd(i) = max(abs(dvzd),sml) + grad%agrt(i) = max(abs(dvzt),sml) + grad%agru(i) = max(abs(dvzu),sml) + grad%agrd(i) = max(abs(dvzd),sml) - dagrzt= dvzt*dvzzt/agrt(i) - dagrzu= dvzu*dvzzu/agru(i) - dagrzd= dvzd*dvzzd/agrd(i) + dagrzt= dvzt*dvzzt/grad%agrt(i) + dagrzu= dvzu*dvzzu/grad%agru(i) + dagrzd= dvzd*dvzzd/grad%agrd(i) - gggrt(i) = dvzt*dagrzt - gggru(i) = dvzu*dagrzu - gggrd(i) = dvzd*dagrzd + grad%gggrt(i) = dvzt*dagrzt + grad%gggru(i) = dvzu*dagrzu + grad%gggrd(i) = dvzd*dagrzd c dztadz=d(zeta)/dz,.. @@ -62,13 +74,13 @@ c dztadz=d(zeta)/dz,.. c gzgr=grad(zeta)*grad(ro). - gzgr(i) = dztadz*dvzt + grad%gzgr(i) = dztadz*dvzt c g2rt: grad(grad(ro)) - g2rt(i) = dvzzt - g2ru(i) = dvzzu - g2rd(i) = dvzzd + grad%g2rt(i) = dvzzt + grad%g2ru(i) = dvzzu + grad%g2rd(i) = dvzzd ENDDO @@ -89,17 +101,17 @@ c g2rt: grad(grad(ro)) c agrt: abs(grad(ro)), u,d for up and down. - agrt(i) = max(abs(dvzt),sml) - agru(i) = max(abs(dvzu),sml) - agrd(i) = max(abs(dvzd),sml) + grad%agrt(i) = max(abs(dvzt),sml) + grad%agru(i) = max(abs(dvzu),sml) + grad%agrd(i) = max(abs(dvzd),sml) - dagrzt= dvzt*dvzzt/agrt(i) - dagrzu= dvzu*dvzzu/agru(i) - dagrzd= dvzd*dvzzd/agrd(i) + dagrzt= dvzt*dvzzt/grad%agrt(i) + dagrzu= dvzu*dvzzu/grad%agru(i) + dagrzd= dvzd*dvzzd/grad%agrd(i) - gggrt(i) = dvzt*dagrzt - gggru(i) = dvzu*dagrzu - gggrd(i) = dvzd*dagrzd + grad%gggrt(i) = dvzt*dagrzt + grad%gggru(i) = dvzu*dagrzu + grad%gggrd(i) = dvzd*dagrzd c dztadz=d(zeta)/dz,.. @@ -107,13 +119,13 @@ c dztadz=d(zeta)/dz,.. c gzgr=grad(zeta)*grad(ro). - gzgr(i) = dztadz*dvzt + grad%gzgr(i) = dztadz*dvzt c g2rt: grad(grad(ro)) - g2rt(i) = dvzzt - g2ru(i) = dvzzu - g2rd(i) = dvzzd + grad%g2rt(i) = dvzzt + grad%g2ru(i) = dvzzu + grad%g2rd(i) = dvzzd ENDDO diff --git a/vgen/vgen_xcpot.F90 b/vgen/vgen_xcpot.F90 index 25349036667b7db93bf807f80cf8eccef3b19b58..2df938861a41b4bf9e8db729ff7f456ab18b9f57 100644 --- a/vgen/vgen_xcpot.F90 +++ b/vgen/vgen_xcpot.F90 @@ -30,7 +30,7 @@ CONTAINS USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_hybrid),INTENT(IN) :: hybrid TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_dimension),INTENT(IN) :: dimension diff --git a/vgen/visxc.f90 b/vgen/visxc.f90 index ea75bd8dc285208ee8cfda796019326dae02e36c..fee171f8d5a565cd695f48765adf67bfb1a75e7d 100644 --- a/vgen/visxc.f90 +++ b/vgen/visxc.f90 @@ -15,7 +15,6 @@ ! ** r.pentcheva 08.05.96 ! ******************************************************************** USE m_types - USE m_xcall, ONLY : vxcall,excall USE m_fft3d IMPLICIT NONE ! .. @@ -23,7 +22,7 @@ INTEGER, INTENT (IN) :: ifftd TYPE(t_stars),INTENT(IN) :: stars TYPE(t_noco),INTENT(IN) :: noco - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_input),INTENT(IN) :: input TYPE(t_potden),INTENT(IN) :: den TYPE(t_potden),INTENT(INOUT) :: vxc,exc,vx @@ -43,27 +42,21 @@ ! ---> allocate arrays ! ALLOCATE ( e_xc(0:ifftd-1),vcon(0:ifftd-1),v_xc(0:ifftd-1,input%jspins),& - & af3(0:ifftd-1,input%jspins),bf3(0:ifftd-1) ) + af3(0:ifftd-1,input%jspins),bf3(0:ifftd-1) ) ALLOCATE( v_x(0:ifftd-1,input%jspins) ) ! ! transform charge density to real space ! DO js = 1,input%jspins - CALL fft3d(& - & af3(0,js),bf3,& - & den%pw(1,js),& - & stars,+1) + CALL fft3d(af3(0,js),bf3, den%pw(1,js), stars,+1) ENDDO IF (noco%l_noco) THEN ALLOCATE (mx(0:ifftd-1),my(0:ifftd-1)) - CALL fft3d(& - & mx,my,& - & den%pw(:,3),& - & stars,+1) + CALL fft3d(mx,my, den%pw(:,3), stars,+1) DO i=0,27*stars%mx1*stars%mx2*stars%mx3-1 chdens= (af3(i,1)+af3(i,2))/2. magmom= mx(i)**2 + my(i)**2 + ((af3(i,1)-af3(i,2))/2.)**2 @@ -79,17 +72,14 @@ ! ! calculate the exchange-correlation potential in real space ! - nt=ifftd - CALL vxcall& - & (6,xcpot,input%jspins,& - & ifftd,nt,af3,& - & v_x,v_xc) + nt=ifftd + CALL xcpot%get_vxc(input%jspins,af3,v_xc,v_x) ! ! ---> back fft to g space and add to coulomb potential for file nrp ! - IF (input%total) THEN + IF (ALLOCATED(exc%pw_w)) THEN DO js = 1,input%jspins @@ -97,10 +87,7 @@ vcon(i)=stars%ufft(i)*v_xc(i,js) bf3(i)=0.0 ENDDO - CALL fft3d(& - & vcon,bf3,& - & fg3,& - & stars,-1) + CALL fft3d(vcon,bf3, fg3, stars,-1) fg3=fg3*stars%nstr DO k = 1,stars%ng3 @@ -111,10 +98,7 @@ vcon(i)=stars%ufft(i)*v_x(i,js) bf3(i)=0.0 ENDDO - CALL fft3d(& - & vcon,bf3,& - & fg3,& - & stars,-1) + CALL fft3d(vcon,bf3, fg3, stars,-1) fg3=fg3*stars%nstr DO k = 1,stars%ng3 @@ -124,11 +108,8 @@ ENDDO ! ! calculate the ex.-cor energy density in real space -! - CALL excall& - & (6,xcpot,input%jspins,& - & ifftd,nt,af3,& - & e_xc) + ! + CALL xcpot%get_exc(input%jspins,af3,e_xc) DO i=0,ifftd-1 vcon(i)=stars%ufft(i)*e_xc(i) @@ -137,28 +118,19 @@ ! ! ---> back fft to g space ! - CALL fft3d(& - & vcon,bf3,& - & exc%pw_w(:,1),& - & stars,-1) + CALL fft3d(vcon,bf3, exc%pw_w(:,1), stars,-1) exc%pw_w(:,1)=exc%pw_w(:,1)*stars%nstr ! END IF ! input%total DO js = 1,input%jspins bf3(0:ifftd-1)=0.0 - CALL fft3d(& - & v_xc(0,js),bf3,& - & fg3,& - & stars,-1) + CALL fft3d(v_xc(0,js),bf3, fg3, stars,-1) DO k = 1,stars%ng3 vxc%pw(k,js) = vxc%pw(k,js) + fg3(k) ENDDO - CALL fft3d(& - & v_x(0,js),bf3,& - & fg3,& - & stars,-1) + CALL fft3d(v_x(0,js),bf3, fg3, stars,-1) DO k = 1,stars%ng3 vx%pw(k,js) = vx%pw(k,js) + fg3(k) ENDDO diff --git a/vgen/visxcg.f90 b/vgen/visxcg.f90 index c61aee9975538afb7148c45e1bddce0a83f9808f..dff91b53d9764a19a49642d4ae19e8ab9dc7593a 100644 --- a/vgen/visxcg.f90 +++ b/vgen/visxcg.f90 @@ -1,3 +1,8 @@ +!-------------------------------------------------------------------------------- +! 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_visxcg USE m_juDFT ! ****************************************************** @@ -6,13 +11,8 @@ MODULE m_visxcg ! including gradient corrections. t.a. 1996. ! ****************************************************** CONTAINS - SUBROUTINE visxcg(& - & ifftd,stars,sym,& - & ifftxc3d,& - & cell,den,& - & xcpot,input,& - & obsolete,noco,& - & vxc,vx,exc) + SUBROUTINE visxcg(ifftd,stars,sym, ifftxc3d, cell,den,& + xcpot,input, obsolete,noco, vxc,vx,exc) ! ****************************************************** ! instead of visxcor.f: the different exchange-correlation @@ -29,13 +29,12 @@ CONTAINS USE m_grdrsis USE m_prpxcfftmap USE m_mkgxyz3 - USE m_xcallg, ONLY : vxcallg,excallg USE m_fft3d USE m_fft3dxc USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_obsolete),INTENT(IN) :: obsolete TYPE(t_input),INTENT(IN) :: input TYPE(t_noco),INTENT(IN) :: noco @@ -51,9 +50,6 @@ CONTAINS ! .. Array Arguments .. REAL rhmni ,d_15,sprsv - ! lwb: if true, white-bird trick. - ! lwbc: l-white-bird-current. - LOGICAL lwbc ! !-----> fft information for xc potential + energy ! @@ -65,6 +61,7 @@ CONTAINS COMPLEX :: ci REAL :: rhotot INTEGER :: ifftd,ifftxc3d + TYPE(t_gradients)::grad ! .. ! .. Local Arrays .. @@ -75,11 +72,7 @@ CONTAINS REAL, ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:) ! REAL, ALLOCATABLE :: v_x(:,:),v_xc(:,:),e_xc(:),vcon(:) - REAL, ALLOCATABLE :: agr(:),agru(:),agrd(:) - REAL, ALLOCATABLE :: g2r(:),g2ru(:),g2rd(:) - REAL, ALLOCATABLE :: gggr(:),gggru(:),gggrd(:) - REAL, ALLOCATABLE :: gzgr(:) - + ! .. unused input (needed for other noco GGA-implementations) .. !ta+ @@ -122,14 +115,10 @@ CONTAINS ifftd=27*stars%mx1*stars%mx2*stars%mx3 ifftxc3d = stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft ALLOCATE ( igxc_fft(0:ifftxc3d-1),gxc_fft(0:ifftxc3d-1,3) ) - CALL prp_xcfft_map(& - & stars,sym,& - & cell,& - & igxc_fft,gxc_fft) + CALL prp_xcfft_map(stars,sym, cell, igxc_fft,gxc_fft) ! ifftxc3=stars%kxc1_fft*stars%kxc2_fft*stars%kxc3_fft - lwbc=obsolete%lwb - + IF (stars%ng3.GT.stars%ng3) THEN WRITE(6,'(/'' stars%ng3.gt.stars%ng3. stars%ng3,stars%ng3='',2i6)') stars%ng3,stars%ng3 CALL juDFT_error("ng3.gt.n3d",calledby="visxcg") @@ -141,41 +130,30 @@ CONTAINS ! ! Allocate arrays ! ff - ALLOCATE( bf3(0:ifftd-1),ph_wrk(0:ifftxc3d-1), & - & rho(0:ifftxc3d-1,input%jspins),rhd1(0:ifftxc3d-1,input%jspins,3),& - & rhd2(0:ifftxc3d-1,input%jspins,6) ) + ALLOCATE( bf3(0:ifftd-1),ph_wrk(0:ifftxc3d-1), rho(0:ifftxc3d-1,input%jspins),& + rhd1(0:ifftxc3d-1,input%jspins,3), rhd2(0:ifftxc3d-1,input%jspins,6) ) IF (noco%l_noco) THEN ALLOCATE( mx(0:ifftxc3-1),my(0:ifftxc3-1),& - & magmom(0:ifftxc3-1), & - & dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) ) + magmom(0:ifftxc3-1), dmagmom(0:ifftxc3-1,3),ddmagmom(0:ifftxc3-1,3,3) ) END IF !--> transform charge density to real space DO js=1,input%jspins - CALL fft3dxc(& - & rho(0:,js),bf3,& - & den%pw(:,js),& - & stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& - & stars%nxc3_fft,stars%kmxxc_fft,+1,& - & stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) + CALL fft3dxc(rho(0:,js),bf3, den%pw(:,js), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& + stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) END DO IF (noco%l_noco) THEN ! for off-diagonal parts the same - CALL fft3dxc(& - & mx,my,& - & den%pw(:,3),& - & stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& - & stars%nxc3_fft,stars%kmxxc_fft,+1,& - & stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) + CALL fft3dxc(mx,my, den%pw(:,3), stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& + stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) DO i=0,ifftxc3-1 rhotot= 0.5*( rho(i,1) + rho(i,2) ) - magmom(i)= SQRT( (0.5*(rho(i,1)-rho(i,2)))**2 & - & + mx(i)**2 + my(i)**2 ) + magmom(i)= SQRT( (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 ) rho(i,1)= rhotot+magmom(i) rho(i,2)= rhotot-magmom(i) END DO @@ -206,20 +184,15 @@ CONTAINS END DO DO js=1,input%jspins - CALL fft3dxc(& - & rhd1(0:,js,idm),bf3,& - & cqpw(:,js),& - & stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1,& - & stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr) + CALL fft3dxc(rhd1(0:,js,idm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,& + stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr) END DO END DO IF (noco%l_noco) THEN - CALL grdrsis(& - & magmom,cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,obsolete,& - & dmagmom ) + CALL grdrsis(magmom,cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,obsolete, dmagmom ) DO i=0,ifftxc3-1 DO idm=1,3 @@ -231,8 +204,7 @@ CONTAINS END IF - IF (lwbc) GOTO 100 - + !--> for dd(rho)/d(xx,xy,yy,zx,yz,zz) = rhd2(:,:,idm) (idm=1,2,3,4,5,6) ! ! ph_wrk: exp(i*(g_x,g_y,g_z)*tau) * g_(x,y,z) * g_(x,y,z) @@ -253,11 +225,8 @@ CONTAINS ENDDO DO js=1,input%jspins - CALL fft3dxc(& - & rhd2(0:,js,ndm),bf3,& - & cqpw(:,js),& - & stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1,& - & stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr) + CALL fft3dxc(rhd2(0:,js,ndm),bf3, cqpw(:,js), stars%kxc1_fft,stars%kxc2_fft,& + stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,+1, stars%igfft(0:,1),igxc_fft,ph_wrk,stars%nstr) END DO END DO ! jdm END DO ! idm @@ -267,9 +236,7 @@ CONTAINS IF (noco%l_noco) THEN DO idm = 1,3 - CALL grdrsis(& - & dmagmom(0,idm),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,obsolete,& - & ddmagmom(0,1,idm) ) + CALL grdrsis(dmagmom(0,idm),cell,stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,obsolete, ddmagmom(0,1,idm) ) END DO ndm= 0 @@ -279,10 +246,8 @@ CONTAINS DO i=0,ifftxc3-1 rhotot= rhd2(i,1,ndm)/2.+rhd2(i,2,ndm)/2. - rhd2(i,1,ndm)= rhotot +& - & ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. - rhd2(i,2,ndm)= rhotot -& - & ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. + rhd2(i,1,ndm)= rhotot + ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. + rhd2(i,2,ndm)= rhotot - ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2. END DO ENDDO !jdm @@ -306,22 +271,14 @@ CONTAINS bf3=0.0 ! allocate the other arrays ! - ALLOCATE (agr(0:ifftxc3d-1),agru(0:ifftxc3d-1),agrd(0:ifftxc3d-1),& - & g2r(0:ifftxc3d-1),g2ru(0:ifftxc3d-1),g2rd(0:ifftxc3d-1),& - & gggr(0:ifftxc3d-1),gggru(0:ifftxc3d-1),gggrd(0:ifftxc3d-1),& - & gzgr(0:ifftxc3d-1)) - + CALL xcpot%alloc_gradients(ifftxc3d,input%jspins,grad) + ! ! calculate the quantities such as abs(grad(rho)),.. used in ! evaluating the gradient contributions to potential and energy. ! - CALL mkgxyz3& - & (ifftxc3d,input%jspins,ifftxc3,input%jspins,rho,& - & rhd1(0,1,1),rhd1(0,1,2),rhd1(0,1,3),& - & rhd2(0,1,1),rhd2(0,1,3),rhd2(0,1,6),& - & rhd2(0,1,5),rhd2(0,1,4),rhd2(0,1,2),& - & agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,& - & gzgr) + CALL mkgxyz3 (ifftxc3d,input%jspins,ifftxc3,input%jspins,rho, rhd1(0,1,1),rhd1(0,1,2),rhd1(0,1,3),& + rhd2(0,1,1),rhd2(0,1,3),rhd2(0,1,6), rhd2(0,1,5),rhd2(0,1,4),rhd2(0,1,2), grad) DEALLOCATE ( rhd1,rhd2 ) ALLOCATE ( v_xc(0:ifftxc3d-1,input%jspins) ) @@ -344,15 +301,10 @@ CONTAINS IF (rhmni.LT.obsolete%chng) THEN - WRITE(6,'(/'' rhmn.lt.obsolete%chng in visxc. rhmn,obsolete%chng='',& - & 2d9.2)') rhmni,obsolete%chng + WRITE(6,'(/'' rhmn.lt.obsolete%chng in visxc. rhmn,obsolete%chng='',2d9.2)') rhmni,obsolete%chng ! CALL juDFT_error("visxcg: rhmn.lt.chng",calledby="visxcg") ENDIF - - CALL vxcallg(xcpot,lwbc,input%jspins,SIZE(agr),nt,rho,& - & agr,agru,agrd,g2r,g2ru,g2rd,& - & gggr,gggru,gggrd,gzgr,& - & v_x,v_xc) + CALL xcpot%get_vxc(input%jspins,rho,v_xc,v_x,grad) ! !----> back fft to g space !----> perform back fft transform: v_xc(r) --> vxc(star) @@ -361,12 +313,8 @@ CONTAINS DO js = 1,input%jspins bf3=0.0 - CALL fft3dxc(& - & v_xc(0:,js),bf3,& - & fg3,& - & stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& - & stars%nxc3_fft,stars%kmxxc_fft,-1,& - & stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) + CALL fft3dxc(v_xc(0:,js),bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& + stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) ! DO k = 1,stars%nxc3_fft vxc%pw(k,js) = vxc%pw(k,js) + fg3(k) @@ -382,10 +330,7 @@ CONTAINS ! fg3(stars%nxc3_fft+1:)=0.0 ALLOCATE ( vcon(0:ifftd-1) ) - CALL fft3d(& - & vcon(0),bf3,& - & fg3,& - & stars,+1) + CALL fft3d(vcon(0),bf3, fg3, stars,+1) ! !----> Convolute with step function ! @@ -393,10 +338,7 @@ CONTAINS vcon(i)=stars%ufft(i)*vcon(i) ENDDO bf3=0.0 - CALL fft3d(& - & vcon(0),bf3,& - & fg3,& - & stars,-1,.FALSE.) + CALL fft3d(vcon(0),bf3, fg3, stars,-1,.FALSE.) DEALLOCATE ( vcon ) ! !----> add to warped coulomb potential @@ -408,12 +350,8 @@ CONTAINS ENDIF bf3=0.0 - CALL fft3dxc(& - & v_x(0:,js),bf3,& - & fg3,& - & stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& - & stars%nxc3_fft,stars%kmxxc_fft,-1,& - & stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) + CALL fft3dxc(v_x(0:,js),bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft, & + stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) ! DO k = 1,stars%nxc3_fft vx%pw(k,js) = vx%pw(k,js) + fg3(k) @@ -430,10 +368,7 @@ CONTAINS fg3(stars%nxc3_fft+1:)=0.0 ALLOCATE ( vcon(0:ifftd-1) ) - CALL fft3d(& - & vcon(0),bf3,& - & fg3,& - & stars,+1) + CALL fft3d(vcon(0),bf3, fg3, stars,+1) ! !----> Convolute with step function ! @@ -441,10 +376,7 @@ CONTAINS vcon(i)=stars%ufft(i)*vcon(i) ENDDO bf3=0.0 - CALL fft3d(& - & vcon(0),bf3,& - & fg3,& - & stars,-1,.FALSE.) + CALL fft3d(vcon(0),bf3, fg3, stars,-1,.FALSE.) DEALLOCATE ( vcon ) ! !----> add to warped exchange-potential @@ -460,21 +392,15 @@ CONTAINS ! ! calculate the ex.-cor energy density in real space ! - IF (input%total) THEN + IF (ALLOCATED(exc%pw_w)) THEN ALLOCATE ( e_xc(0:ifftxc3d-1) ) - CALL excallg(xcpot,lwbc,input%jspins,nt,rho,& - & agr,agru,agrd,g2r,g2ru,g2rd,& - & gggr,gggru,gggrd,gzgr,& - & e_xc) + CALL xcpot%get_exc(input%jspins,rho,e_xc,grad) ! !----> perform back fft transform: exc(r) --> exc(star) ! bf3=0.0 - CALL fft3dxc(& - & e_xc,bf3,& - & fg3,& - & stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,stars%nxc3_fft,stars%kmxxc_fft,-1,& - & stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) + CALL fft3dxc(e_xc,bf3, fg3, stars%kxc1_fft,stars%kxc2_fft,stars%kxc3_fft,& + stars%nxc3_fft,stars%kmxxc_fft,-1, stars%igfft(0:,1),igxc_fft,stars%pgfft,stars%nstr) DEALLOCATE ( e_xc ) ! !----> Perform fft transform: exc(star) --> exc(r) @@ -483,8 +409,7 @@ CONTAINS fg3(stars%nxc3_fft+1:)=0.0 bf3=0.0 ALLOCATE ( vcon(0:ifftd-1) ) - CALL fft3d(vcon,bf3,fg3,& - & stars,+1) + CALL fft3d(vcon,bf3,fg3, stars,+1) DO i=0,ifftd-1 vcon(i)=stars%ufft(i)*vcon(i) @@ -493,16 +418,14 @@ CONTAINS ! ---> back fft to g space ! bf3=0.0 - CALL fft3d(vcon,bf3,exc%pw_w(:,1),& - & stars,-1,.FALSE.) + CALL fft3d(vcon,bf3,exc%pw_w(:,1), stars,-1,.FALSE.) DEALLOCATE ( vcon ) ! ENDIF DEALLOCATE(fg3) DEALLOCATE ( bf3,rho,igxc_fft,gxc_fft ) - DEALLOCATE ( agr,agru,agrd,g2r,g2ru,g2rd,& - & gggr,gggru,gggrd,gzgr ) + diff --git a/vgen/vmtxc.f90 b/vgen/vmtxc.f90 index 4ac09b99e7eca60b582b0e77363909d310c10f6c..77e0e081e392916020223fe20fb5c3ab5d497bc4 100644 --- a/vgen/vmtxc.f90 +++ b/vgen/vmtxc.f90 @@ -12,17 +12,16 @@ CONTAINS USE m_lhglpts USE m_gaussp - USE m_xcall, ONLY : vxcall,excall USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_dimension),INTENT(IN) :: DIMENSION TYPE(t_input),INTENT(IN) :: input TYPE(t_sym),INTENT(IN) :: sym TYPE(t_sphhar),INTENT(IN) :: sphhar TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_potden),INTENT(IN) :: den - TYPE(t_potden),INTENT(INOUT) :: vxc,exc + TYPE(t_potden),INTENT(INOUT) :: vxc,exc TYPE(t_potden),INTENT(INOUT),optional :: vx ! .. @@ -76,7 +75,7 @@ CONTAINS ENDDO ! calculate the ex.-cor. potential ! - CALL vxcall (6,xcpot,input%jspins, nsp,nsp,rhoxc, v_x,v_xc) + CALL xcpot%get_vxc(input%jspins,rhoxc(:nsp,:), v_xc,v_x) ! ----> now determine the corresponding potential number DO js = 1,input%jspins ! @@ -109,32 +108,32 @@ CONTAINS ENDIF ENDDO ENDDO - IF (input%total) THEN + IF (ALLOCATED(exc%mt)) THEN ! ! calculate the ex.-cor energy density ! - CALL excall(6,xcpot,input%jspins, nsp,nsp,rhoxc, e_ex) - END IF - ! ----> now determine the corresponding energy density number - ! - ! ---> multiplikate exc with the weights of the k-points - ! - DO k = 1,nsp - e_ex(k) = e_ex(k)*wt(k) - ENDDO - DO lh = 0,sphhar%nlh(nd) - + CALL xcpot%get_exc(input%jspins,rhoxc(:nsp,:), e_ex) + ! ----> now determine the corresponding energy density number ! - ! ---> determine the corresponding potential number through gauss integration + ! ---> multiplikate exc with the weights of the k-points ! - elh=DOT_PRODUCT(e_ex(:nsp),ylh(:nsp,lh,nd)) - - ! - ! ---> add to the given potential - - exc%mt(jr,lh,n,1) = elh - ENDDO - + DO k = 1,nsp + e_ex(k) = e_ex(k)*wt(k) + ENDDO + DO lh = 0,sphhar%nlh(nd) + + ! + ! ---> determine the corresponding potential number through gauss integration + ! + elh=DOT_PRODUCT(e_ex(:nsp),ylh(:nsp,lh,nd)) + + ! + ! ---> add to the given potential + + exc%mt(jr,lh,n,1) = elh + ENDDO + END IF + ENDDO ENDDO diff --git a/vgen/vmtxcg.F90 b/vgen/vmtxcg.F90 index a65327cadc748f09a34e121a569a010fb36327b4..e8cf3e4f792625bb3eee29a81c8a15fcb946131b 100644 --- a/vgen/vmtxcg.F90 +++ b/vgen/vmtxcg.F90 @@ -22,23 +22,19 @@ MODULE m_vmtxcg ! ********************************************************* CONTAINS - SUBROUTINE vmtxcg(& - & dimension,mpi,sphhar,atoms,& - & den,xcpot,input,sym,& - & obsolete,& - & vxc,vx,exc) + SUBROUTINE vmtxcg(dimension,mpi,sphhar,atoms,& + den,xcpot,input,sym, obsolete,vxc,vx,exc) #include"cpp_double.h" USE m_lhglptg USE m_grdchlh USE m_mkgylm USE m_gaussp - USE m_xcallg, ONLY : vxcallg,excallg - + USE m_types_xcpot_inbuild USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_dimension),INTENT(IN) :: dimension TYPE(t_mpi),INTENT(IN) :: mpi TYPE(t_obsolete),INTENT(IN) :: obsolete @@ -53,6 +49,7 @@ CONTAINS #endif ! .. ! .. Local Scalars .. + TYPE(t_gradients)::grad INTEGER jr,js,k,lh,n,nd,ist,nsp,ixpm ,i,nat REAL rhmnm,d_15,elh,vlh LOGICAL lwbc ! if true, white-bird trick @@ -60,11 +57,8 @@ CONTAINS ! .. Local Arrays .. REAL v_x(dimension%nspd,dimension%jspd),v_xc(dimension%nspd,dimension%jspd),e_xc(dimension%nspd),rx(3,dimension%nspd) REAL vxcl(DIMENSION%nspd,DIMENSION%jspd),excl(DIMENSION%nspd),divi - TYPE(t_xcpot)::xcpot_tmp + TYPE(t_xcpot_inbuild)::xcpot_tmp REAL wt(dimension%nspd),rr2(atoms%jmtd),thet(dimension%nspd) - REAL agr(dimension%nspd),agru(dimension%nspd),agrd(dimension%nspd),g2r(dimension%nspd),g2ru(dimension%nspd) - REAL g2rd(dimension%nspd),gggr(dimension%nspd),gggru(dimension%nspd),gggrd(dimension%nspd) - REAL grgru(dimension%nspd),grgrd(dimension%nspd),gzgr(dimension%nspd) REAL, ALLOCATABLE :: ylh(:,:,:),ylht(:,:,:),ylhtt(:,:,:) REAL, ALLOCATABLE :: ylhf(:,:,:),ylhff(:,:,:),ylhtf(:,:,:) REAL, ALLOCATABLE :: chlh(:,:,:),chlhdr(:,:,:),chlhdrr(:,:,:) @@ -83,10 +77,8 @@ CONTAINS ! .. ! .. Intrinsic Functions .. INTRINSIC max,mod,min - LOGICAL l_gga !ta+ !.....------------------------------------------------------------------ - l_gga=xcpot%is_gga() #ifdef CPP_MPI CALL MPI_BCAST(obsolete%ndvgrd,1,MPI_INTEGER,0,mpi%mpi_comm,ierr) @@ -105,9 +97,7 @@ CONTAINS ! angular mesh equidistant in phi, ! theta are zeros of the legendre polynomials ! - CALL gaussp(& - & atoms%lmaxd,& - & rx,wt) + CALL gaussp(atoms%lmaxd, rx,wt) nsp = dimension%nspd ! @@ -139,10 +129,15 @@ CONTAINS & chlhdr(atoms%jmtd,0:sphhar%nlhd,dimension%jspd),chlhdrr(atoms%jmtd,0:sphhar%nlhd,dimension%jspd)) DO 200 n = n_start,atoms%ntype,n_stride - IF (xcpot%lda_atom(n))THEN - IF((.NOT.xcpot%is_name("pw91"))) CALL judft_warn("Using locally LDA only possible with pw91 functional") - CALL xcpot_tmp%init("l91",.FALSE.) - ENDIF + SELECT TYPE(xcpot) + TYPE IS(t_xcpot_inbuild) + IF (xcpot%lda_atom(n))THEN + IF((.NOT.xcpot%is_name("pw91"))) & + CALL judft_warn("Using locally LDA only possible with pw91 functional") + CALL xcpot_tmp%init("l91",.FALSE.) + ENDIF + END SELECT + nat=sum(atoms%neq(:n-1))+1 nd = atoms%ntypsy(nat) @@ -163,38 +158,35 @@ CONTAINS chlh(jr,lh,js) = den%mt(jr,lh,n,js)*rr2(jr) ENDDO - IF (l_gga) THEN - CALL grdchlh(& - & ixpm,ist,atoms%jri(n),atoms%dx(n),atoms%rmsh(1,n),& - & chlh(1,lh,js),obsolete%ndvgrd,& - & chlhdr(1,lh,js),chlhdrr(1,lh,js)) - ENDIF - + CALL grdchlh(ixpm,ist,atoms%jri(n),atoms%dx(n),atoms%rmsh(1,n),& + chlh(1,lh,js),obsolete%ndvgrd, chlhdr(1,lh,js),chlhdrr(1,lh,js)) + ENDDO ! js ENDDO ! lh ! !--> loop over radial mesh ! - !$OMP PARALLEL DEFAULT(none) & + !!$OMP PARALLEL DEFAULT(none) & #ifdef CPP_MPI - !$OMP& SHARED(vr_local,vxr_local,excr_local) & + !!$OMP& SHARED(vr_local,vxr_local,excr_local) & #endif - !$OMP& SHARED(vxc,vx,exc,l_gga) & - !$OMP& SHARED(dimension,mpi,sphhar,atoms,den,xcpot,input,sym,obsolete)& - !$OMP& SHARED(n,nd,ist,ixpm,nsp,nat,d_15,lwbc) & - !$OMP& SHARED(rx,wt,rr2,thet,xcpot_tmp) & - !$OMP& SHARED(ylh,ylht,ylhtt,ylhf,ylhff,ylhtf) & - !$OMP& SHARED(chlh,chlhdr,chlhdrr) & - !$OMP& SHARED(ierr,n_start,n_stride) & - !$OMP& PRIVATE(js,k,lh,i,rhmnm,elh,vlh) & - !$OMP& PRIVATE(v_x,v_xc,e_xc,excl,vxcl,divi) & - !$OMP& PRIVATE(agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,grgru,grgrd,gzgr) & - !$OMP& PRIVATE(ch,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf) - ALLOCATE ( ch(dimension%nspd,dimension%jspd),chdr(dimension%nspd,dimension%jspd),chdt(dimension%nspd,dimension%jspd),& - & chdf(dimension%nspd,dimension%jspd),chdrr(dimension%nspd,dimension%jspd),chdtt(dimension%nspd,dimension%jspd),& - & chdff(dimension%nspd,dimension%jspd),chdtf(dimension%nspd,dimension%jspd),chdrt(dimension%nspd,dimension%jspd),& - & chdrf(dimension%nspd,dimension%jspd)) + !!$OMP& SHARED(vxc,vx,exc) & + !!$OMP& SHARED(dimension,mpi,sphhar,atoms,den,xcpot,input,sym,obsolete)& + !!$OMP& SHARED(n,nd,ist,ixpm,nsp,nat,d_15,lwbc) & + !!$OMP& SHARED(rx,wt,rr2,thet,xcpot_tmp) & + !!$OMP& SHARED(ylh,ylht,ylhtt,ylhf,ylhff,ylhtf) & + !!$OMP& SHARED(chlh,chlhdr,chlhdrr) & + !!$OMP& SHARED(ierr,n_start,n_stride) & + !!$OMP& PRIVATE(js,k,lh,i,rhmnm,elh,vlh) & + !!$OMP& PRIVATE(v_x,v_xc,e_xc,excl,vxcl,divi) & + !!$OMP& PRIVATE(grad) & + !!$OMP& PRIVATE(ch,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf) + ALLOCATE ( ch(DIMENSION%nspd,DIMENSION%jspd),chdr(DIMENSION%nspd,DIMENSION%jspd),& + chdt(DIMENSION%nspd,DIMENSION%jspd), chdf(DIMENSION%nspd,DIMENSION%jspd),& + chdrr(DIMENSION%nspd,DIMENSION%jspd),chdtt(DIMENSION%nspd,DIMENSION%jspd), & + chdff(DIMENSION%nspd,DIMENSION%jspd),chdtf(DIMENSION%nspd,DIMENSION%jspd),& + chdrt(DIMENSION%nspd,DIMENSION%jspd), chdrf(DIMENSION%nspd,DIMENSION%jspd)) !$OMP DO DO 190 jr = 1,atoms%jri(n) ! @@ -221,7 +213,6 @@ CONTAINS ch(k,js) = ch(k,js) + ylh(k,lh,nd)*chlh(jr,lh,js) ENDDO - IF (l_gga) THEN ! DO k = 1,nsp chdr(k,js) =chdr(k,js)+ ylh(k,lh,nd)*chlhdr(jr,lh,js) @@ -239,24 +230,13 @@ CONTAINS ENDDO - ENDIF - ENDDO ! lh ENDDO ! js - IF (l_gga) THEN - CALL mkgylm(& - & input%jspins,atoms%rmsh(jr,n),thet,nsp,dimension%nspd,dimension%jspd,ch,chdr,& - & chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf,& - & agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,& - & grgru,grgrd,gzgr)!keep - ELSE - agr(:) = 0.0 ; agru(:) = 0.0 ; agrd(:) = 0.0 - g2r(:) = 0.0 ; g2ru(:) = 0.0 ; g2rd(:) = 0.0 - gggr(:) = 0.0 ; gggru(:) = 0.0 ; gggrd(:) = 0.0 - grgru(:) = 0.0 ; grgrd(:) = 0.0 ; gzgr(:) = 0.0 - ENDIF - + CALL xcpot%alloc_gradients(input%jspins,nsp,grad) + CALL mkgylm(input%jspins,atoms%rmsh(jr,n),thet,nsp,DIMENSION%nspd,DIMENSION%jspd,ch,chdr,& + chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf,grad) + ! ! rhmnm: rho_minimum_muffin-tin.. @@ -279,29 +259,22 @@ CONTAINS ! ! calculate the ex.-cor. potential - - CALL vxcallg(& - & xcpot,lwbc,input%jspins,nsp,nsp,ch,agr,agru,agrd,& - & g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,& - & v_x,v_xc) - - IF (xcpot%lda_atom(n)) THEN - ! Use local part of pw91 for this atom - CALL vxcallg(& - xcpot_tmp,lwbc,input%jspins,nsp,nsp,ch,agr,agru,agrd,& - g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,& - v_x,vxcl) - !Mix the potentials - divi = 1.0 / (atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(1,n)) - v_xc(:,:) = ( vxcl(:,:) * ( atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(jr,n) ) +& + CALL xcpot%get_vxc(input%jspins,ch(:nsp,:),v_xc,v_x,grad) + SELECT TYPE(xcpot) + TYPE IS (t_xcpot_inbuild) + IF (xcpot%lda_atom(n)) THEN + ! Use local part of pw91 for this atom + CALL xcpot_tmp%get_vxc(input%jspins,ch(:nsp,:),vxcl,v_x,grad) + !Mix the potentials + divi = 1.0 / (atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(1,n)) + v_xc(:,:) = ( vxcl(:,:) * ( atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(jr,n) ) +& v_xc(:,:) * ( atoms%rmsh(jr,n) - atoms%rmsh(1,n) ) ) * divi - ENDIF - + ENDIF + END SELECT IF (mpi%irank == 0) THEN - IF (mod(jr,1000).eq.0)& - & WRITE (6,'(/'' 999vxc''/(10d15.7))')& - & ((v_xc(k,js),k=1,nsp),js=1,input%jspins) + IF (mod(jr,1000).eq.0) WRITE (6,'(/'' 999vxc''/(10d15.7))')& + ((v_xc(k,js),k=1,nsp),js=1,input%jspins) ENDIF !irank==0 @@ -317,10 +290,8 @@ CONTAINS ENDDO IF (mpi%irank == 0) THEN - IF (mod(jr,1500).EQ.0)& - & WRITE (6,'('' 999wt''/(10d15.7))') (wt(k),k=1,nsp) - IF (mod(jr,1500).EQ.0)& - & WRITE (6,'('' 999vxc''/(10d15.7))') (v_xc(k,js),k=1,nsp) + IF (mod(jr,1500).EQ.0) WRITE (6,'('' 999wt''/(10d15.7))') (wt(k),k=1,nsp) + IF (mod(jr,1500).EQ.0) WRITE (6,'('' 999vxc''/(10d15.7))') (v_xc(k,js),k=1,nsp) ENDIF !irank==0 DO lh = 0,sphhar%nlh(nd) @@ -354,29 +325,25 @@ CONTAINS ENDDO ! lh ENDDO ! js - IF (input%total) then + IF (ALLOCATED(exc%mt)) THEN ! ! calculate the ex.-cor energy density ! - CALL excallg(xcpot,lwbc,input%jspins,nsp,& - ch,agr,agru,agrd,g2r,g2ru,g2rd,& - gggr,gggru,gggrd,gzgr, e_xc) - - IF (xcpot%lda_atom(n)) THEN - ! Use local part of pw91 for this atom - CALL excallg(xcpot_tmp,lwbc,input%jspins,nsp,& - ch,agr,agru,agrd,g2r,g2ru,g2rd,& - gggr,gggru,gggrd,gzgr, excl) - !Mix the potentials - e_xc(:) = ( excl(:) * ( atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(jr,n) ) +& - e_xc(:) * ( atoms%rmsh(jr,n) - atoms%rmsh(1,n) ) ) * divi - ENDIF + CALL xcpot%get_exc(input%jspins,ch(:nsp,:),e_xc,grad) + SELECT TYPE(xcpot) + TYPE IS(t_xcpot_inbuild) + IF (xcpot%lda_atom(n)) THEN + ! Use local part of pw91 for this atom + CALL xcpot_tmp%get_exc(input%jspins,ch(:nsp,:),excl,grad) + !Mix the potentials + e_xc(:) = ( excl(:) * ( atoms%rmsh(atoms%jri(n),n) - atoms%rmsh(jr,n) ) +& + e_xc(:) * ( atoms%rmsh(jr,n) - atoms%rmsh(1,n) ) ) * divi + ENDIF + END SELECT IF (mpi%irank == 0) THEN - IF (mod(jr,10000).EQ.0)& - & WRITE (6,'(/'' 999exc''/(10d15.7))') (e_xc(k),k=1,nsp) + IF (mod(jr,10000).EQ.0) WRITE (6,'(/'' 999exc''/(10d15.7))') (e_xc(k),k=1,nsp) ENDIF !irank==0 - ENDIF ! now determine the corresponding energy density number ! @@ -402,11 +369,12 @@ CONTAINS #endif ENDDO + ENDIF 190 ENDDO - !$OMP END DO + !!$OMP END DO DEALLOCATE (ch,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf) - !$OMP END PARALLEL + !!$OMP END PARALLEL ! WRITE(6,'(/'' n='',i3/'' 9999vr''/(10d15.7))') n, ! & (((vr(jr,lh,n,js),jr=1,jri(n),100),lh=0,ntypsy(nat)),js=1,jspins) diff --git a/vgen/vvacxc.f90 b/vgen/vvacxc.f90 index 58cabef21cccd8e0830e881db532850e0bf50099..7826f2d7639ba7f2d90f453771ebfa1e3ff98629 100644 --- a/vgen/vvacxc.f90 +++ b/vgen/vvacxc.f90 @@ -16,11 +16,10 @@ CONTAINS ! ** r.pentcheva 08.05.96 ! ******************************************************************** - USE m_xcall, ONLY : vxcall,excall USE m_fft2d USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_input),INTENT(IN) :: input TYPE(t_vacuum),INTENT(IN) :: vacuum TYPE(t_noco),INTENT(IN) :: noco @@ -67,25 +66,17 @@ CONTAINS ! transform charge density to real space: 2-dim FFT ! DO js = 1,input%jspins - CALL fft2d(& - & stars,& - & af2(0,js),bf2,& - & den%vacz(ip,ivac,js),rhti,den%vacxy(ip,1,ivac,js),& - & vacuum%nmzxyd,+1) + CALL fft2d(stars, af2(0,js),bf2, den%vacz(ip,ivac,js),rhti,& + den%vacxy(ip,1,ivac,js), vacuum%nmzxyd,+1) END DO IF (noco%l_noco) THEN - CALL fft2d(& - & stars,& - & mx,my, & - & den%vacz(ip,ivac,3),den%vacz(ip,ivac,4),& - & den%vacxy(ip,1,ivac,3),& - & vacuum%nmzxyd,1) + CALL fft2d(stars, mx,my, den%vacz(ip,ivac,3),den%vacz(ip,ivac,4),& + den%vacxy(ip,1,ivac,3), vacuum%nmzxyd,1) DO i=0,9*stars%mx1*stars%mx2-1 chdens= (af2(i,1)+af2(i,2))/2. - magmom= mx(i)**2 + my(i)**2 +& - & ((af2(i,1)-af2(i,2))/2.)**2 + magmom= mx(i)**2 + my(i)**2 + ((af2(i,1)-af2(i,2))/2.)**2 magmom= SQRT(magmom) af2(i,1)= chdens + magmom af2(i,2)= chdens - magmom @@ -95,21 +86,14 @@ CONTAINS ! ! calculate the exchange-correlation potential in real space ! - CALL vxcall& - & (6,xcpot,input%jspins,& - & ifftd2,nt,af2,& - & v_x,v_xc) + CALL xcpot%get_vxc(input%jspins,af2,v_xc,v_x) DO js = 1,input%jspins ! ! ----> 2-d back fft to g space ! bf2=0.0 - CALL fft2d(& - & stars,& - & v_xc(0,js),bf2,& - & fgz,rhti,fgxy,& - & 1,-1) + CALL fft2d(stars, v_xc(0,js),bf2, fgz,rhti,fgxy, 1,-1) ! ! ----> and add vxc to coulomb potential ! the G||.eq.zero component is added to vz @@ -119,27 +103,19 @@ CONTAINS ! the G||.ne.zero components are added to vxc%vacxy ! DO irec2 = 1,stars%ng2-1 - vxc%vacxy(ip,irec2,ivac,js) = vxc%vacxy(ip,irec2,ivac,js) +& - & fgxy(irec2) + vxc%vacxy(ip,irec2,ivac,js) = vxc%vacxy(ip,irec2,ivac,js) + fgxy(irec2) ENDDO ENDDO ! !i calculate the exchange-correlation energy density in real space ! - IF (input%total) THEN - CALL excall& - & (6,xcpot,input%jspins,& - & ifftd2,nt,af2,& - & e_xc) + IF (ALLOCATED(exc%vacz)) THEN + call xcpot%get_exc(input%jspins,af2,e_xc) ! ! ----> 2-d back fft to g space ! bf2=0.0 - CALL fft2d(& - & stars,& - & e_xc,bf2,& - & exc%vacz(ip,ivac,1),rhti,exc%vacxy(ip,1,ivac,1),& - & vacuum%nmzxyd,-1) + CALL fft2d(stars, e_xc,bf2, exc%vacz(ip,ivac,1),rhti,exc%vacxy(ip,1,ivac,1), vacuum%nmzxyd,-1) ENDIF ENDDO @@ -159,8 +135,7 @@ CONTAINS mx(0)= den%vacz(vacuum%nmzxy+k,ivac,3) my(0)= den%vacz(vacuum%nmzxy+k,ivac,4) chdens= (af2(k-1,1)+af2(k-1,2))/2. - magmom= mx(0)**2 + my(0)**2 +& - & ((af2(k-1,1)-af2(k-1,2))/2.)**2 + magmom= mx(0)**2 + my(0)**2 + ((af2(k-1,1)-af2(k-1,2))/2.)**2 magmom= SQRT(magmom) af2(k-1,1)= chdens + magmom af2(k-1,2)= chdens - magmom @@ -169,10 +144,8 @@ CONTAINS ENDDO - CALL vxcall& - & (6,xcpot,input%jspins,& - & vacuum%nmzd,nmzdiff,af2,& - & vxz,vxcz) + CALL xcpot%get_vxc(input%jspins,af2,vxz,vxcz) + !+gu DO js=1,input%jspins DO k=vacuum%nmzxy+1,vacuum%nmz @@ -186,11 +159,8 @@ CONTAINS ! ! calculate the ex.-corr. energy density now beyond warping region ! - IF (input%total) THEN - CALL excall& - & (6,xcpot,input%jspins,& - & vacuum%nmzd,nmzdiff,af2,& - & exc%vacz(vacuum%nmzxy+1,ivac,1)) + IF (ALLOCATED(exc%vacz)) THEN + CALL xcpot%get_exc(input%jspins,af2,exc%vacz(vacuum%nmzxy+1:,ivac,1)) END IF ENDDO IF (noco%l_noco) THEN diff --git a/vgen/vvacxcg.f90 b/vgen/vvacxcg.f90 index 7d7aeac0f4a45852f4b20f50524cb9acd40fb3b0..dd60692d567005c69525d54d2ca834ac9fdbcaa9 100644 --- a/vgen/vvacxcg.f90 +++ b/vgen/vvacxcg.f90 @@ -7,10 +7,8 @@ MODULE m_vvacxcg ! for the gradient contribution. t.a. 1996 !----------------------------------------------------------------------- CONTAINS - SUBROUTINE vvacxcg(& - & ifftd2,stars,vacuum,noco,oneD,& - & cell,xcpot,input,obsolete,den,& - & vxc,exc) + SUBROUTINE vvacxcg(ifftd2,stars,vacuum,noco,oneD,& + cell,xcpot,input,obsolete,den, vxc,exc) !----------------------------------------------------------------------- ! instead of vvacxcor.f: the different exchange-correlation @@ -27,10 +25,9 @@ CONTAINS USE m_od_mkgxyz3 USE m_od_mkgz USE m_fft2d - USE m_xcallg, ONLY : vxcallg,excallg USE m_types IMPLICIT NONE - TYPE(t_xcpot),INTENT(IN) :: xcpot + CLASS(t_xcpot),INTENT(IN) :: xcpot TYPE(t_oneD),INTENT(IN) :: oneD TYPE(t_obsolete),INTENT(IN) :: obsolete TYPE(t_input),INTENT(IN) :: input @@ -50,12 +47,10 @@ CONTAINS INTEGER :: js,nt,i,iq,irec2,nmz0,nmzdiff,ivac,ip REAL :: rhti,zro,fgz,rhmnv,d_15,bmat1(3,3),rd COMPLEX :: ci - LOGICAL :: lwbc ! if true, white-bird trick. ! .. ! .. Local Arrays .. - REAL, ALLOCATABLE :: af2(:,:),bf2(:),agr(:),agru(:),agrd(:),g2r(:) - REAL, ALLOCATABLE :: g2ru(:),g2rd(:),gggr(:),gggru(:),gggrd(:) - REAL, ALLOCATABLE :: gzgr(:),rhdx(:,:),rhdy(:,:),rhdz(:,:) + REAL, ALLOCATABLE :: af2(:,:),bf2(:) + REAL, ALLOCATABLE :: rhdx(:,:),rhdy(:,:),rhdz(:,:) REAL, ALLOCATABLE :: rhdxx(:,:),rhdyy(:,:),rhtdz(:,:),rhtdzz(:,:) REAL, ALLOCATABLE :: rhdzz(:,:),rhdyz(:,:),rhdzx(:,:),rhdxy(:,:) REAL, ALLOCATABLE :: v_x(:,:),v_xc(:,:),e_xc(:),vxz(:,:),vxcz(:,:) @@ -63,21 +58,21 @@ CONTAINS REAL, ALLOCATABLE :: rxydzzr(:),rxydzzi(:),rhtxyr(:),rhtxyi(:) REAL, ALLOCATABLE :: rhtxc(:,:),rhtz(:,:),dummy(:) COMPLEX, ALLOCATABLE :: fgxy(:),rxydz(:,:,:),rxydzz(:,:,:),cqpw(:) + + TYPE(t_gradients)::grad ! .. ! for the noco-case only REAL :: chdens - REAL, ALLOCATABLE :: magmom(:,:),& - & dxmagmom(:),ddxmagmom(:,:),& - & dymagmom(:),ddymagmom(:,:), & - & dzmagmom(:,:),ddzmagmom(:,:) + REAL, ALLOCATABLE :: magmom(:,:), dxmagmom(:),ddxmagmom(:,:) + REAL, ALLOCATABLE :: dymagmom(:),ddymagmom(:,:), dzmagmom(:,:),ddzmagmom(:,:) REAL, ALLOCATABLE :: mx(:),my(:) REAL, ALLOCATABLE :: af2noco(:,:,:) ! .. unused input (needed for other noco GGA-implementations) .. + - lwbc = .false. d_15 = 1.e-15 zro = 0.0 ci = cmplx(0.,1.) @@ -121,22 +116,14 @@ CONTAINS DO ip=1,vacuum%nmzxy DO js=1,input%jspins - CALL fft2d(& - & stars,& - & af2noco(0,ip,js),bf2,& - & den%vacz(ip,ivac,js),0.,den%vacxy(ip,1,ivac,js),& - & vacuum%nmzxy,+1) + CALL fft2d(stars, af2noco(0,ip,js),bf2, den%vacz(ip,ivac,js),0.,& + den%vacxy(ip,1,ivac,js), vacuum%nmzxy,+1) END DO - CALL fft2d(& - & stars,& - & mx,my,& - & den%vacz(ip,ivac,3),den%vacz(ip,ivac,4),& - & den%vacxy(ip,1,ivac,3),& - & vacuum%nmzxy,+1) + CALL fft2d(stars, mx,my, den%vacz(ip,ivac,3),den%vacz(ip,ivac,4), & + den%vacxy(ip,1,ivac,3), vacuum%nmzxy,+1) DO i=0,9*stars%mx1*stars%mx2-1 - magmom(i,ip)= mx(i)**2 + my(i)**2 +& - & ((af2noco(i,ip,1)-af2noco(i,ip,2))/2.)**2 + magmom(i,ip)= mx(i)**2 + my(i)**2 + ((af2noco(i,ip,1)-af2noco(i,ip,2))/2.)**2 magmom(i,ip)= SQRT(magmom(i,ip)) chdens= af2noco(i,ip,1)/2.+af2noco(i,ip,2)/2. af2noco(i,ip,1)= chdens + magmom(i,ip) @@ -156,9 +143,8 @@ CONTAINS ! calculate first (rhtdz) & second (rhtdzz) derivative of den%vacz(1:nmz) ! ALLOCATE ( dummy(vacuum%nmz) ) - CALL grdchlh(& - & 0,1,vacuum%nmz,vacuum%delz,dummy,den%vacz(1,ivac,js),obsolete%ndvgrd,& - & rhtdz(1,js),rhtdzz(1,js)) + CALL grdchlh(0,1,vacuum%nmz,vacuum%delz,dummy,den%vacz(1,ivac,js),obsolete%ndvgrd,& + rhtdz(1,js),rhtdzz(1,js)) DEALLOCATE ( dummy ) ALLOCATE ( rhtxyr(vacuum%nmzxy), rhtxyi(vacuum%nmzxy),dummy(vacuum%nmzxy) ) ALLOCATE ( rxydzr(vacuum%nmzxy), rxydzi(vacuum%nmzxy) ) @@ -171,16 +157,12 @@ CONTAINS DO ip=1,vacuum%nmzxy rhtxyr(ip)=den%vacxy(ip,iq,ivac,js) ENDDO - CALL grdchlh(& - & 0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyr,obsolete%ndvgrd,& - & rxydzr,rxydzzr) + CALL grdchlh(0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyr,obsolete%ndvgrd, rxydzr,rxydzzr) DO ip=1,vacuum%nmzxy rhtxyi(ip)=aimag(den%vacxy(ip,iq,ivac,js)) ENDDO - CALL grdchlh(& - & 0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyi,obsolete%ndvgrd,& - & rxydzi,rxydzzi) + CALL grdchlh(0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyi,obsolete%ndvgrd, rxydzi,rxydzzi) DO ip=1,vacuum%nmzxy rxydz(ip,iq,js)=cmplx(rxydzr(ip),rxydzi(ip)) @@ -203,9 +185,7 @@ CONTAINS DO ip=1,vacuum%nmzxy rhtxyr(ip)=magmom(i,ip) ENDDO - CALL grdchlh(& - & 0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyr,obsolete%ndvgrd,& - & rxydzr,rxydzzr) + CALL grdchlh(0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyr,obsolete%ndvgrd, rxydzr,rxydzzr) DO ip=1,vacuum%nmzxy dzmagmom(i,ip)= rxydzr(ip) ddzmagmom(i,ip)= rxydzzr(ip) @@ -221,10 +201,9 @@ CONTAINS ALLOCATE ( rhdz(0:ifftd2-1,input%jspins),rhdxx(0:ifftd2-1,input%jspins) ) ALLOCATE ( rhdyy(0:ifftd2-1,input%jspins),rhdzz(0:ifftd2-1,input%jspins) ) ALLOCATE ( rhdyz(0:ifftd2-1,input%jspins),rhdzx(0:ifftd2-1,input%jspins) ) - ALLOCATE ( rhdxy(0:ifftd2-1,input%jspins),gggrd(0:ifftd2-1) ) - ALLOCATE ( agr(0:ifftd2-1),agru(0:ifftd2-1),agrd(0:ifftd2-1) ) - ALLOCATE ( g2r(0:ifftd2-1),g2ru(0:ifftd2-1),g2rd(0:ifftd2-1) ) - ALLOCATE ( gggr(0:ifftd2-1),gggru(0:ifftd2-1),gzgr(0:ifftd2-1) ) + ALLOCATE ( rhdxy(0:ifftd2-1,input%jspins)) + + CALL xcpot%alloc_gradients(ifftd2,input%jspins,grad) ALLOCATE ( v_xc(0:ifftd2-1,input%jspins),e_xc(0:ifftd2-1) ) ALLOCATE ( cqpw(stars%ng2-1),af2(0:ifftd2-1,input%jspins) ) ALLOCATE ( fgxy(stars%ng2-1),bf2(0:ifftd2-1) ) @@ -241,11 +220,8 @@ CONTAINS ! Transform charge and magnetization to real-space. DO js=1,input%jspins - CALL fft2d(& - & stars,& - & af2(0,js),bf2,& - & den%vacz(ip,ivac,js),0.,den%vacxy(ip,1,ivac,js),& - & vacuum%nmzxyd,+1) + CALL fft2d(stars, af2(0,js),bf2, den%vacz(ip,ivac,js),0.,& + den%vacxy(ip,1,ivac,js), vacuum%nmzxyd,+1) END DO ELSE @@ -273,27 +249,16 @@ CONTAINS - CALL fft2d( & - &stars, & - &rhdx(0,js),bf2,& - &zro,rhti,cqpw, & - &1,+1,stars%ft2_gfx) + CALL fft2d(stars, rhdx(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfx) !TODO & pgft2x) rhti = 0.0 CALL fft2d( & ! dn/dy = FFT(0,i*gy*den%vacxy)& - & stars,& - & rhdy(0,js),bf2,& - & zro,rhti,cqpw,& - & 1,+1,stars%ft2_gfy) + stars, rhdy(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfy) rhti = 0.0 CALL fft2d( & ! dn/dz = FFT(rhtdz,rxydz)& - & stars,& - & rhdz(0,js),bf2,& - & rhtdz(ip,js),rhti,rxydz(ip,1,js),& - & vacuum%nmzxyd,+1) - + stars, rhdz(0,js),bf2, rhtdz(ip,js),rhti,rxydz(ip,1,js), vacuum%nmzxyd,+1) DO iq=1,stars%ng2-1 cqpw(iq)=-den%vacxy(ip,iq,ivac,js) @@ -301,24 +266,15 @@ CONTAINS rhti = 0.0 CALL fft2d( & ! d2n/dx2 = FFT(0,-gx^2*den%vacxy)& - & stars,& - & rhdxx(0,js),bf2,& - & zro,rhti,cqpw,& - & 1,+1,stars%ft2_gfx*stars%ft2_gfx) + stars, rhdxx(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfx*stars%ft2_gfx) rhti = 0.0 CALL fft2d( & ! d2n/dy2 = FFT(0,-gy^2*den%vacxy)& - & stars,& - & rhdyy(0,js),bf2,& - & zro,rhti,cqpw,& - & 1,+1,stars%ft2_gfy*stars%ft2_gfy) + stars, rhdyy(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfy*stars%ft2_gfy) rhti = 0.0 CALL fft2d( & ! d2n/dz2 = FFT(rhtdzz,rxydzz)& - & stars,& - & rhdzz(0,js),bf2,& - & rhtdzz(ip,js),rhti,rxydzz(ip,1,js),& - & vacuum%nmzxyd,+1) + stars, rhdzz(0,js),bf2, rhtdzz(ip,js),rhti,rxydzz(ip,1,js), vacuum%nmzxyd,+1) DO iq=1,stars%ng2-1 @@ -327,17 +283,11 @@ CONTAINS rhti = 0.0 CALL fft2d( & ! d2n/dyz = FFT(0,i*gy*rxydz)& - & stars,& - & rhdyz(0,js),bf2,& - & zro,rhti,cqpw,& - & 1,+1,stars%ft2_gfy) + stars, rhdyz(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfy) rhti = 0.0 CALL fft2d( & ! d2n/dzx = FFT(0,i*gx*rxydz)& - & stars,& - & rhdzx(0,js),bf2,& - & zro,rhti,cqpw,& - & 1,+1,stars%ft2_gfx) + stars, rhdzx(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfx) DO iq=1,stars%ng2-1 cqpw(iq)=-den%vacxy(ip,iq,ivac,js) @@ -345,10 +295,7 @@ CONTAINS rhti = 0.0 CALL fft2d( & ! d2n/dxy = FFT(0,-gx*gy*den%vacxy)& - & stars,& - & rhdxy(0,js),bf2,& - & zro,rhti,cqpw,& - & 1,+1,stars%ft2_gfy*stars%ft2_gfx) + stars, rhdxy(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfy*stars%ft2_gfx) END DO ! js=1,input%jspins @@ -358,9 +305,7 @@ CONTAINS ! ! in real-space. The derivatives of the charge density, that are ! ! already calculated in g-space, will be used. - CALL grdrsvac(& - & magmom(0,ip),bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd,& - & dxmagmom,dymagmom) + CALL grdrsvac(magmom(0,ip),bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, dxmagmom,dymagmom) DO i=0,9*stars%mx1*stars%mx2-1 chdens= rhdx(i,1)/2.+rhdx(i,2)/2. rhdx(i,1)= chdens + dxmagmom(i) @@ -373,12 +318,10 @@ CONTAINS rhdz(i,2)= chdens - dzmagmom(i,ip) END DO + CALL grdrsvac(dxmagmom,bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, & + ddxmagmom(0,1),ddymagmom(0,1)) CALL grdrsvac(& - & dxmagmom,bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, & - & ddxmagmom(0,1),ddymagmom(0,1)) - CALL grdrsvac(& - & dymagmom,bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, & - & ddxmagmom(0,2),ddymagmom(0,2)) + dymagmom,bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd,ddxmagmom(0,2),ddymagmom(0,2)) DO i=0,9*stars%mx1*stars%mx2-1 chdens= rhdxx(i,1)/2.+rhdxx(i,2)/2. rhdxx(i,1)= chdens + ddxmagmom(i,1) @@ -387,14 +330,11 @@ CONTAINS rhdyy(i,1)= chdens + ddymagmom(i,2) rhdyy(i,2)= chdens - ddymagmom(i,2) chdens= rhdxy(i,1)/2.+rhdxy(i,2)/2. - rhdxy(i,1)= chdens + & - & ( ddxmagmom(i,2) + ddymagmom(i,1) )/2. - rhdxy(i,2)= chdens - & - & ( ddxmagmom(i,2) + ddymagmom(i,1) )/2. + rhdxy(i,1)= chdens + ( ddxmagmom(i,2) + ddymagmom(i,1) )/2. + rhdxy(i,2)= chdens - ( ddxmagmom(i,2) + ddymagmom(i,1) )/2. END DO - CALL grdrsvac(& - & dzmagmom(0,ip),bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, & - & ddxmagmom(0,1),ddymagmom(0,1)) + CALL grdrsvac(dzmagmom(0,ip),bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, & + ddxmagmom(0,1),ddymagmom(0,1)) DO i=0,9*stars%mx1*stars%mx2-1 chdens= rhdzx(i,1)/2.+rhdzx(i,2)/2. rhdzx(i,1)= chdens + ddxmagmom(i,1) @@ -423,18 +363,16 @@ CONTAINS rd = cell%z1 + vacuum%delz*(ip-1) if(oneD%odi%d1)then - CALL od_mkgxyz3(& - & ifftd2,input%jspins,ifftd2,input%jspins,& - & af2,rd,rhdx,rhdy,rhdz,rhdxx,rhdyy,& - & rhdzz,rhdyz,rhdzx,rhdxy,& - & agr,agru,agrd,g2r,g2ru,g2rd,& - & gggr,gggru,gggrd,gzgr) - else - CALL mkgxyz3(& - & ifftd2,input%jspins,ifftd2,input%jspins,af2,rhdx,rhdy,& - & rhdz,rhdxx,rhdyy,rhdzz,rhdyz,rhdzx,rhdxy,& - & agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,& - & gzgr) +!!$ CALL od_mkgxyz3(& +!!$ & ifftd2,input%jspins,ifftd2,input%jspins,& +!!$ & af2,rd,rhdx,rhdy,rhdz,rhdxx,rhdyy,& +!!$ & rhdzz,rhdyz,rhdzx,rhdxy,& +!!$ & agr,agru,agrd,g2r,g2ru,g2rd,& +!!$ & gggr,gggru,gggrd,gzgr) + CALL judft_error("OneD not implemented") + ELSE + CALL mkgxyz3(ifftd2,input%jspins,ifftd2,input%jspins,af2,rhdx,rhdy,& + rhdz,rhdxx,rhdyy,rhdzz,rhdyz,rhdzx,rhdxy, grad) endif ! rhmnv: rho_minimum_vacuum. @@ -454,24 +392,15 @@ CONTAINS ! calculate the exchange-correlation potential in real space ! - - CALL vxcallg(& - & xcpot,lwbc,input%jspins,nt,nt,af2,agr,agru,agrd,& - & g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,& - & v_x,v_xc) - - + CALL xcpot%get_vxc(input%jspins,af2,v_xc,v_x,grad) + DO js = 1,input%jspins ! ! ----> 2-d back fft to g space ! bf2=0.0 - CALL fft2d(& - & stars,& - & v_xc(0,js),bf2,& - & fgz,rhti,fgxy,& - & 1,-1) + CALL fft2d(stars, v_xc(0,js),bf2, fgz,rhti,fgxy, 1,-1) ! ----> and add vxc to coulomb potential ! the g||.eq.zero component is added to vxc%vacz @@ -489,21 +418,14 @@ CONTAINS ! calculate the exchange-correlation energy density in real space ! - IF (input%total) THEN - - CALL excallg(& - & xcpot,lwbc,input%jspins,nt,af2,agr,agru,agrd,& - & g2r,g2ru,g2rd, gggr,gggru,gggrd,gzgr,& - & e_xc) + IF (ALLOCATED(exc%vacz)) THEN + CALL xcpot%get_exc(input%jspins,af2,e_xc,grad) + ! ----> 2-d back fft to g space ! bf2=0.0 - CALL fft2d(& - & stars,& - & e_xc,bf2,& - & exc%vacz(ip,ivac,1),rhti,exc%vacxy(ip,1,ivac,1),& - & vacuum%nmzxyd,-1) + CALL fft2d(stars, e_xc,bf2, exc%vacz(ip,ivac,1),rhti,exc%vacxy(ip,1,ivac,1), vacuum%nmzxyd,-1) ENDIF @@ -543,8 +465,7 @@ CONTAINS mx(0) = den%vacz(ip,ivac,3) my(0) = den%vacz(ip,ivac,4) chdens= (af2(0,1)+af2(0,2))/2. - magmom(0,1)= mx(0)**2 + my(0)**2 +& - & ((af2(0,1)-af2(0,2))/2.)**2 + magmom(0,1)= mx(0)**2 + my(0)**2 + ((af2(0,1)-af2(0,2))/2.)**2 magmom(0,1)= SQRT(magmom(0,1)) rhtz(ip,1)= chdens + magmom(0,1) rhtz(ip,2)= chdens - magmom(0,1) @@ -555,9 +476,8 @@ CONTAINS DEALLOCATE ( magmom,mx,my ) ALLOCATE ( dummy(vacuum%nmz) ) DO js=1,input%jspins - CALL grdchlh(& - & 0,1,vacuum%nmz-nmz0+1,vacuum%delz,dummy,rhtz(nmz0,js),obsolete%ndvgrd,& - & rhtdz(nmz0,js),rhtdzz(nmz0,js)) + CALL grdchlh(0,1,vacuum%nmz-nmz0+1,vacuum%delz,dummy,rhtz(nmz0,js),obsolete%ndvgrd,& + rhtdz(nmz0,js),rhtdzz(nmz0,js)) END DO DEALLOCATE ( dummy ) @@ -567,28 +487,22 @@ CONTAINS !c evaluating the gradient contributions to potential and !c energy. - agr(:)=0.0 ; agru(:)=0.0 ; agrd(:)=0.0 ; g2r(:)=0.0 - g2ru(:)=0.0 ; g2rd(:)=0.0 ; gggr(:)=0.0 ; gggru(:)=0.0 - gggrd(:)=0.0 ; gzgr(:)=0.0 - + CALL xcpot%alloc_gradients(nmzdiff,input%jspins,grad) IF (xcpot%is_gga()) THEN if(oneD%odi%d1)then - CALL od_mkgz(& - & cell%z1,vacuum%nmzxy,vacuum%delz,& - & nmzdiff,input%jspins,& - & rhtz(vacuum%nmzxy+1,1),rhtz(vacuum%nmzxy+1,input%jspins),& - & rhtdz(vacuum%nmzxy+1,1), rhtdz(vacuum%nmzxy+1,input%jspins),& - & rhtdzz(vacuum%nmzxy+1,1),rhtdzz(vacuum%nmzxy+1,input%jspins),& - & agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,& - & gzgr) - - else - CALL mkgz(nmzdiff,input%jspins,& - & rhtz(vacuum%nmzxy+1,1),rhtz(vacuum%nmzxy+1,input%jspins),& - & rhtdz(vacuum%nmzxy+1,1),rhtdz(vacuum%nmzxy+1,input%jspins),& - & rhtdzz(vacuum%nmzxy+1,1),rhtdzz(vacuum%nmzxy+1,input%jspins),& - & agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,& - & gzgr) +! CALL od_mkgz(& +! cell%z1,vacuum%nmzxy,vacuum%delz,& +! nmzdiff,input%jspins,& +! rhtz(vacuum%nmzxy+1,1),rhtz(vacuum%nmzxy+1,input%jspins),& +! rhtdz(vacuum%nmzxy+1,1), rhtdz(vacuum%nmzxy+1,input%jspins),& +! rhtdzz(vacuum%nmzxy+1,1),rhtdzz(vacuum%nmzxy+1,input%jspins),& +! agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,& +! gzgr) + CALL judft_error("OneD not implemented") + ELSE + CALL mkgz(nmzdiff,input%jspins, rhtz(vacuum%nmzxy+1,1),rhtz(vacuum%nmzxy+1,input%jspins),& + rhtdz(vacuum%nmzxy+1,1),rhtdz(vacuum%nmzxy+1,input%jspins),rhtdzz(vacuum%nmzxy+1,1),& + rhtdzz(vacuum%nmzxy+1,input%jspins),grad) endif ENDIF @@ -613,11 +527,8 @@ CONTAINS ! CALL juDFT_error("vvacxcg: rhmn.lt.chng",calledby="vvacxcg") ENDIF - CALL vxcallg(& - & xcpot,lwbc,input%jspins,vacuum%nmzd,nmzdiff,rhtxc,agr(:vacuum%nmzd),agru,& - & agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,& - & vxz,vxcz) - + CALL xcpot%get_vxc(input%jspins,rhtxc,vxcz,vxz,grad) + DO js = 1,input%jspins DO ip = vacuum%nmzxy + 1,vacuum%nmz vxc%vacz(ip,ivac,js) = vxc%vacz(ip,ivac,js) + vxcz(ip-vacuum%nmzxy,js) @@ -631,24 +542,16 @@ CONTAINS ! ! calculate the ex-corr. energy density now beyond warping region ! - IF (input%total) THEN - CALL excallg(& - & xcpot,lwbc,input%jspins,nmzdiff,rhtxc,agr,agru,& - & agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,& - & exc%vacz(vacuum%nmzxy+1,ivac,1)) + IF (ALLOCATED(exc%vacz)) THEN + CALL xcpot%get_exc(input%jspins,rhtxc,exc%vacz(vacuum%nmzxy+1:,ivac,1),grad) ENDIF DEALLOCATE ( bf2) - DEALLOCATE ( agr,agru ) - DEALLOCATE ( agrd,g2r ) - DEALLOCATE ( g2ru,g2rd ) - DEALLOCATE ( gggr,gggru ) - DEALLOCATE ( gggrd,gzgr,v_x,v_xc,e_xc,vxz,vxcz,rhtz,rhtxc ) + DEALLOCATE ( v_x,v_xc,e_xc,vxz,vxcz,rhtz,rhtxc ) ENDDO ! loop over vacua (ivac) - DEALLOCATE ( rhtdz,rhtdzz,rxydz,rxydzz ) - + END SUBROUTINE vvacxcg diff --git a/xc-pot/CMakeLists.txt b/xc-pot/CMakeLists.txt index 19fbede3ca602838de1978becd8a1cc1baf2514f..c0978f5affa8c26b3792fe08feadf19887d05213 100644 --- a/xc-pot/CMakeLists.txt +++ b/xc-pot/CMakeLists.txt @@ -18,8 +18,6 @@ xc-pot/vxcepbe.f xc-pot/vxcl91.f xc-pot/vxcpw91.f xc-pot/vxcwb91.f -xc-pot/xcall.f -xc-pot/xcallg.f xc-pot/xcbh.f xc-pot/xch91.F xc-pot/xcpz.f @@ -28,3 +26,31 @@ xc-pot/xcwgn.f xc-pot/xcxal.f ) +set(inpgen_F77 ${inpgen_F77} +xc-pot/corg91.F +xc-pot/corl91.f +xc-pot/corpbe.f +xc-pot/easypbe.f +xc-pot/excepbe.f +xc-pot/exchpbe.F +xc-pot/excl91.f +xc-pot/excpw91.f +xc-pot/excwb91.f +xc-pot/gaunt.f +xc-pot/grdchlh.f +xc-pot/pbecor2.f +xc-pot/potl0.f +xc-pot/relcor.f +xc-pot/vxcepbe.f +xc-pot/vxcl91.f +xc-pot/vxcpw91.f +xc-pot/vxcwb91.f +xc-pot/xcbh.f +xc-pot/xch91.F +xc-pot/xcpz.f +xc-pot/xcvwn.f +xc-pot/xcwgn.f +xc-pot/xcxal.f +) + + diff --git a/xc-pot/corpbe.f b/xc-pot/corpbe.f index a14d7350fb999a8a768f1b527c051b11c9ceb680..a8e8e6419444dfd142e7ffbfd6a5a6c1889bb4fd 100644 --- a/xc-pot/corpbe.f +++ b/xc-pot/corpbe.f @@ -14,7 +14,7 @@ c [c] j. p. perdew and y. wang, phys. rev. b {\bf 45}, 13244 (1992). c---------------------------------------------------------------------- CONTAINS SUBROUTINE corpbe( - > xcpot,rs,zet,t,uu,vv,ww,lgga,lpot, + > l_pbes,rs,zet,t,uu,vv,ww,lgga,lpot, < ec,vcup,vcdn,h,dvcup,dvcdn) c---------------------------------------------------------------------- c input: rs=seitz radius=(3/4pi rho)^(1/3) @@ -35,9 +35,8 @@ c : dvcdn=nonlocal correction to vcdn c---------------------------------------------------------------------- USE m_pbecor2 - USE m_types IMPLICIT NONE - type(t_xcpot),intent(in)::xcpot + LOGICAL,INTENT(IN) :: l_pbes INTEGER, INTENT (IN) :: lgga,lpot REAL, INTENT (IN) :: rs,zet,t,uu,vv,ww REAL, INTENT (OUT) :: dvcdn,dvcup,ec,h,vcdn,vcup @@ -79,7 +78,7 @@ c construct ec, using [c](8) + t2,t4,t6,z4,delt,bet c .. - IF (xcpot%is_name("PBEs")) THEN ! PBE_sol + IF (l_pbes) THEN ! PBE_sol bet=0.046e0 ELSE bet=0.06672455060314922e0 diff --git a/xc-pot/easypbe.f b/xc-pot/easypbe.f index aea7119cfc953468b417e134c927f97a0ced6778..f48c615ab1707f51f693c2bfd80d55f665fe7f07 100644 --- a/xc-pot/easypbe.f +++ b/xc-pot/easypbe.f @@ -19,12 +19,12 @@ c---------------------------------------------------------------------- USE m_exchpbe USE m_corpbe - USE m_types + USE m_types_xcpot_data IMPLICIT NONE ! .. Arguments .. - type(t_xcpot),intent(in)::xcpot - INTEGER, INTENT (IN) :: lcor ! flag to do correlation(=0=>don't) + TYPE(t_xcpot_data),INTENT(IN)::xcpot + INTEGER, INTENT (IN) :: lcor ! flag to do correlation(=0=>don't) INTEGER, INTENT (IN) :: lpot ! flag to do potential (=0=>don't) REAL, INTENT (IN) :: up,dn ! density (spin up & down) REAL, INTENT (IN) :: agrup,agrdn ! |grad up|, |grad down| @@ -158,7 +158,7 @@ c9999- ww = (agrup**2-agrdn**2-zet*agr**2)/ (rho*rho*twoksg**2) CALL corpbe( - > xcpot,rs,zet,t,uu,vv,ww,1,lpot, + > xcpot%is_PBEs,rs,zet,t,uu,vv,ww,1,lpot, < ec,vcup,vcdn,h,dvcup,dvcdn) eclsd = ec diff --git a/xc-pot/excepbe.f b/xc-pot/excepbe.f index 8601d8cd8b86ae9ba3bd384cb2c6d284196e8cf2..84fcf03ab46f57426493b8350f7da2398f4cce81 100644 --- a/xc-pot/excepbe.f +++ b/xc-pot/excepbe.f @@ -11,12 +11,12 @@ c.....------------------------------------------------------------------ < exc) USE m_easypbe - USE m_types + USE m_types_xcpot_data IMPLICIT NONE ! .. Arguments .. - TYPE(t_xcpot),INTENT(IN)::xcpot + TYPE(t_xcpot_data),INTENT(IN)::xcpot INTEGER, INTENT (IN) :: irmx,jspins,mirm REAL, INTENT (IN) :: rh(mirm,jspins) REAL, INTENT (IN) :: agr(mirm),agru(mirm),agrd(mirm) @@ -96,7 +96,7 @@ c..... ENDIF ! ro > smlc - IF( xcpot%is_name("pbe0") ) THEN + IF( xcpot%is_pbe0 ) THEN !pbe0: weight exchange energy with factor 0.75 xced = (0.75*xedl+cedl+0.75*xedg+cedg) ELSE diff --git a/xc-pot/exchpbe.F b/xc-pot/exchpbe.F index 42b9e631cce5f35278f78533626dfcc6568732f4..3d1bf9d48a74fca20fb79d5f1f37246eeeb21ba2 100644 --- a/xc-pot/exchpbe.F +++ b/xc-pot/exchpbe.F @@ -16,14 +16,14 @@ c---------------------------------------------------------------------- CONTAINS SUBROUTINE exchpbe(xcpot,rho,s,u,v,lgga,lpot, < ex,vx,vx_sr) - USE m_hsefunctional, ONLY: calculateEnhancementFactor +! USE m_hsefunctional, ONLY: calculateEnhancementFactor USE m_constants, ONLY: pi_const - USE m_types + USE m_types_xcpot_data USE m_judft IMPLICIT NONE ! .. Arguments - TYPE(t_xcpot), INTENT (IN) :: xcpot + TYPE(t_xcpot_data), INTENT (IN) :: xcpot INTEGER, INTENT (IN) :: lgga ! =0=>don't put in gradient corrections, just lda INTEGER, INTENT (IN) :: lpot ! =0=>don't get potential and don't need u and v REAL, INTENT (IN) :: rho ! density @@ -108,9 +108,11 @@ c -gu ! Calculate the enhancement factor fxhse and its derivatives ! as integral over the exchange hole (cf. [d]) kF = (3.0 * pi_const**2 * rho)**thrd - CALL calculateEnhancementFactor(kF, s, fxhse, dFx_ds, d2Fx_ds2, - & dFx_dkF, d2Fx_dsdkF) - ex = exunif * (fxpbe - xcpot%get_exchange_weight() * fxhse ) + CALL judft_error("HSE not implemented",calledby="exchpbe") + !This creates a depency loop +! CALL calculateEnhancementFactor(kF, s, fxhse, dFx_ds, d2Fx_ds2, +! & dFx_dkF, d2Fx_dsdkF) + ex = exunif * (fxpbe - xcpot%exchange_weight * fxhse ) ELSE ex = exunif*fxpbe END IF diff --git a/xc-pot/mkgl0.f b/xc-pot/mkgl0.f index aeab1d7e92e337158b6db5c059dc3d8b0f9b8366..0d29bac0727509eff7ceb497e6c8f1dfc00d3527 100644 --- a/xc-pot/mkgl0.f +++ b/xc-pot/mkgl0.f @@ -11,25 +11,32 @@ c.....------------------------------------------------------------------ CONTAINS SUBROUTINE mkgl0( > mshd,msh,jspd,jspins,rad, densi,drri,ddrri, - < agr,agru,agrd, g2r,g2ru,g2rd, - < gggr,gggru,gggrd, grgru,grgrd, gzgr) + < grad) + USE m_types IMPLICIT NONE INTEGER, INTENT (IN) :: mshd,msh,jspd,jspins REAL, INTENT (IN) :: rad(mshd),densi(mshd,jspd) REAL, INTENT (IN) :: drri(mshd,jspd),ddrri(mshd,jspd) - REAL, INTENT (OUT) :: agr(mshd),agru(mshd),agrd(mshd) - REAL, INTENT (OUT) :: g2r(mshd),g2ru(mshd),g2rd(mshd) - REAL, INTENT (OUT) :: gggr(mshd),gggru(mshd),gggrd(mshd) - REAL, INTENT (OUT) :: grgru(mshd),grgrd(mshd),gzgr(mshd) -c + TYPE(t_gradients)::grad +c REAL dagrr,dagrru,ddrr,ddrrd,ddrru,drr,drrd,drru,dzdr,ro, + rod,rou,rv,spnf INTEGER i c.....------------------------------------------------------------------ c .. spnf = 1./(3-jspins) + IF (allocated(grad%sigma)) THEN + DO i=1,msh + grad%sigma(1,i)=spnf * drri(i,1)*spnf * drri(i,1) + IF (jspins>1) THEN + grad%sigma(2,i)=spnf * drri(i,1)*spnf * drri(i,2) + grad%sigma(3,i)=spnf * drri(i,2)*spnf * drri(i,2) + ENDIF + ENDDO + RETURN + ENDIF DO i = 1,msh rv = rad(i) @@ -44,27 +51,27 @@ c .. ddrru = spnf * ddrri(i,1) ddrrd = spnf * ddrri(i,jspins) - agr(i) = abs(drr) - agru(i) = abs(drru) - agrd(i) = agru(i) + grad%agrt(i) = abs(drr) + grad%agru(i) = abs(drru) + grad%agrd(i) = grad%agru(i) - dagrr = drr*ddrr/agr(i) - dagrru = drru*ddrru/agru(i) + dagrr = drr*ddrr/grad%agrt(i) + dagrru = drru*ddrru/grad%agru(i) - gggr(i) = drr*dagrr - gggru(i) = drru*dagrru - gggrd(i) = gggru(i) + grad%gggrt(i) = drr*dagrr + grad%gggru(i) = drru*dagrru + grad%gggrd(i) = grad%gggru(i) dzdr = ((drru-drrd)*ro- (rou-rod)*drr)/ro**2 - gzgr(i) = dzdr*drr + grad%gzgr(i) = dzdr*drr - g2r(i) = ddrr + 2*drr/rv - g2ru(i) = ddrru + 2*drru/rv - g2rd(i) = ddrrd + 2*drrd/rv + grad%g2rt(i) = ddrr + 2*drr/rv + grad%g2ru(i) = ddrru + 2*drru/rv + grad%g2rd(i) = ddrrd + 2*drrd/rv - grgru(i) = drr*drru - grgrd(i) = drr*drrd + grad%grgru(i) = drr*drru + grad%grgrd(i) = drr*drrd ENDDO diff --git a/xc-pot/potl0.f b/xc-pot/potl0.f index c3860b357abc744f107634c6215d44ef67026f96..198879fb69ad71d83a02154b19a7e6710904b183 100644 --- a/xc-pot/potl0.f +++ b/xc-pot/potl0.f @@ -10,11 +10,10 @@ c ****************************************************************** USE m_grdchlh USE m_mkgl0 - USE m_xcallg , ONLY : vxcallg USE m_types IMPLICIT NONE c .. - TYPE(t_xcpot),intent(in)::xcpot + CLASS(t_xcpot),intent(in)::xcpot INTEGER, INTENT (IN) :: jspins,jspd,mshd,msh REAL, INTENT (IN) :: dx REAL, INTENT (IN) :: rad(msh),dens(mshd,jspd) @@ -23,23 +22,18 @@ c .. c .. previously untyped names .. INTEGER,PARAMETER :: ndvgrd=6 + TYPE(t_gradients)::grad + INTEGER i,ispin - REAL, ALLOCATABLE :: drr(:,:),ddrr(:,:),agrt(:),agrd(:),agru(:) - REAL, ALLOCATABLE :: g2rt(:),g2rd(:),g2ru(:),gggrt(:),gggrd(:) - REAL, ALLOCATABLE :: gggru(:),grgrd(:),grgru(:),gzgr(:) + REAL, ALLOCATABLE :: drr(:,:),ddrr(:,:) REAL :: vx(mshd,jspd) - ALLOCATE ( drr(mshd,jspd),ddrr(mshd,jspd),grgru(mshd),gzgr(mshd) ) - ALLOCATE ( agrt(mshd),agrd(mshd),agru(mshd),g2rt(mshd),g2rd(mshd), - + g2ru(mshd),gggrt(mshd),gggrd(mshd),gggru(mshd),grgrd(mshd) ) - - agrt(:) = 0.0 ; agru(:) = 0.0 ; agrd(:) = 0.0 ; grgrd(:) = 0.0 - g2rt(:) = 0.0 ; g2ru(:) = 0.0 ; g2rd(:) = 0.0 ; gzgr(:) = 0.0 - gggrt(:) = 0.0 ; gggru(:) = 0.0 ; gggrd(:) = 0.0 ; grgru(:) = 0.0 + ALLOCATE ( drr(mshd,jspd),ddrr(mshd,jspd)) ! !--> evaluate gradients of dens. ! + CALL xcpot%alloc_gradients(msh,1,grad) DO ispin = 1, jspins CALL grdchlh( > 1,1,msh,dx,rad,dens(1,ispin),ndvgrd, @@ -48,17 +42,12 @@ c .. previously untyped names .. CALL mkgl0( > mshd,msh,jspd,jspins,rad,dens,drr,ddrr, - < agrt,agru,agrd,g2rt,g2ru,g2rd, - < gggrt,gggru,gggrd,grgru,grgrd,gzgr) + < grad) ! ! --> calculate the potential. ! - CALL vxcallg( - > xcpot,.false.,jspins,mshd,msh,dens, - + agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd, - + gzgr,vx,vxc) + call xcpot%get_vxc(jspins,dens,vxc,vx,grad) - DEALLOCATE ( drr,ddrr,grgru,gzgr,agrt,agrd,agru,g2rt,g2rd, - + g2ru,gggrt,gggrd,gggru,grgrd ) + DEALLOCATE ( drr,ddrr ) END SUBROUTINE potl0 END MODULE m_potl0 diff --git a/xc-pot/vxcepbe.f b/xc-pot/vxcepbe.f index 608af83c5c89ce61327be05152166f77eeae1247..01b673fd35c55771f5c019cfe22660e510c50eb2 100644 --- a/xc-pot/vxcepbe.f +++ b/xc-pot/vxcepbe.f @@ -11,12 +11,12 @@ c.....------------------------------------------------------------------ < vx,vxc) Use m_easypbe - USE m_types + USE m_types_xcpot_data IMPLICIT NONE ! .. Arguments .. - TYPE(t_xcpot),INTENT(IN)::xcpot + TYPE(t_xcpot_data),INTENT(IN)::xcpot INTEGER, INTENT (IN) :: irmx,jspins,mirm REAL, INTENT (IN) :: rh(mirm,jspins) REAL, INTENT (IN) :: agr(mirm),agru(mirm),agrd(mirm) @@ -37,8 +37,7 @@ c.....------------------------------------------------------------------ REAL, PARAMETER :: sml = 1.e-14 REAL, PARAMETER :: smlc = 2.01e-14 LOGICAL :: l_hse - l_hse=(xcpot%is_name("hse").or.xcpot%is_name("vhse").or. - + xcpot%is_name("lhse")) + l_hse=(xcpot%is_hse) !$OMP PARALLEL DEFAULT(none) !$OMP+ SHARED(irmx,rh,xcpot,jspins,l_hse) diff --git a/xc-pot/xcbh.f b/xc-pot/xcbh.f index 66bcc34d7562c98381a44a9f659dbaf74aeb08b1..e3488c9701bf3ae102f75ce6b804802f2ca7cae2 100644 --- a/xc-pot/xcbh.f +++ b/xc-pot/xcbh.f @@ -49,11 +49,11 @@ > mgrid,ngrid,rh, < vx,vxc) !************************************************************************ - USE m_types + USE m_types_xcpot_data ! ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: jspins - TYPE(t_xcpot), INTENT (IN) :: xcpot ! run mode parameters + TYPE(t_xcpot_data), INTENT (IN) :: xcpot ! run mode parameters INTEGER, INTENT (IN) :: iofile ! file number for read and write INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points ! @@ -84,10 +84,10 @@ ! !-----> select exchange correlation potential ! - IF (xcpot%is_name("mjw")) THEN + IF (xcpot%is_mjw) THEN cp = cpmjw ; cf = cfmjw rp = rpmjw ; rf = rfmjw - ELSEIF (xcpot%is_name("bh")) THEN + ELSEIF (xcpot%is_bh) THEN cp = cpvbh ; cf = cfvbh rp = rpvbh ; rf = rfvbh ELSE @@ -158,11 +158,11 @@ C*********************************************************************** > mgrid,ngrid,rh, < exc) C*********************************************************************** - USE m_types + USE m_types_xcpot_data ! ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: jspins - TYPE(t_xcpot), INTENT (IN) :: xcpot ! run mode parameters + TYPE(t_xcpot_data), INTENT (IN) :: xcpot ! run mode parameters INTEGER, INTENT (IN) :: iofile ! file number for read and write INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points ! @@ -192,10 +192,10 @@ C*********************************************************************** ! !-----> select exchange correlation potential ! - IF (xcpot%is_name("mjw")) THEN + IF (xcpot%is_mjw) THEN cp = cpmjw ; cf = cfmjw rp = rpmjw ; rf = rfmjw - ELSEIF (xcpot%is_name("bh")) THEN + ELSEIF (xcpot%is_bh) THEN cp = cpvbh ; cf = cfvbh rp = rpvbh ; rf = rfvbh ELSE