Commit 6626c30d authored by Daniel Wortmann's avatar Daniel Wortmann

Made fleur compatible with new inpgen

parent 9bf7c6c2
......@@ -13,7 +13,7 @@ SUBROUTINE genNewNocoInp(input,atoms,noco,noco_new)
USE m_juDFT
USE m_types
USE m_constants
USE m_rwnoco
!USE m_rwnoco
IMPLICIT NONE
......@@ -47,9 +47,10 @@ SUBROUTINE genNewNocoInp(input,atoms,noco,noco_new)
iatom= iatom + atoms%neq(iType)
END DO
CALL judft_error("BUG:noco-write feature not implemented at present")
OPEN (24,file='nocoinp',form='formatted', status='unknown')
REWIND (24)
CALL rw_noco_write(atoms,noco_new, input)
!CALL rw_noco_write(atoms,noco_new, input)
CLOSE (24)
END SUBROUTINE genNewNocoInp
......
......@@ -86,8 +86,7 @@
bkTemp(:,n) = syp(:,nosyp)
kpts%specialPointIndices(nosyp) = n
kpts%nkpt = kpts%nkpt
kpts%posScale = 1.0
IF(l_fillArrays) THEN
IF(ALLOCATED(kpts%bk)) THEN
DEALLOCATE(kpts%bk)
......@@ -101,12 +100,12 @@
ELSE
OPEN (41,file='kpts',form='formatted',status='new')
IF (.NOT.input%film) THEN
WRITE(41,'(i5,f20.10)') kpts%nkpt,kpts%posScale
WRITE(41,'(i5,f20.10)') kpts%nkpt
DO n = 1, kpts%nkpt
WRITE (41,'(4f10.5)') bkTemp(:,n),1.0
ENDDO
ELSE
WRITE(41,'(i5,f20.10,3x,l1)') kpts%nkpt,kpts%posScale,.false.
WRITE(41,'(i5,3x,l1)') kpts%nkpt,.false.
DO n = 1, kpts%nkpt
WRITE (41,'(3f10.5)') bkTemp(1:2,n),1.0
ENDDO
......
......@@ -34,9 +34,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,kpts,&
USE m_efield
USE m_od_mapatom
USE m_od_kptsgen
USE m_gen_bz
USE m_nocoInputCheck
USE m_kpoints
USE m_types_forcetheo_extended
USE m_relaxio
......@@ -227,8 +225,6 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,kpts,&
END IF
END IF
! Calculate missing kpts parameters
CALL kpoints(oneD,sym,cell,input,noco,banddos,kpts,l_kpts)
! Generate missing general parameters
......@@ -476,13 +472,7 @@ SUBROUTINE postprocessInput(mpi,input,field,sym,stars,atoms,vacuum,kpts,&
!adjust positions by displacements
CALL apply_displacements(cell,input,vacuum,oneD,sym,noco,atoms)
!Calculate kpoint in the full BZ
IF (kpts%l_gamma.and. banddos%ndir .eq. 0.and.kpts%specificationType==2) THEN
CALL gen_bz(kpts,sym)
ELSE
kpts%nkptf=0
ENDIF
! Missing xc functionals initializations
IF (xcpot%needs_grad()) THEN
ALLOCATE (stars%ft2_gfx(0:stars%kimax2),stars%ft2_gfy(0:stars%kimax2))
......
......@@ -24,7 +24,6 @@ atompar.F90
film_sym.f90
make_atomic_defaults.f90
read_inpgen_input.f90
closure.f90
super_check.f90
make_kpoints.f90
read_old_inp.f90
......@@ -58,7 +57,6 @@ ${FLEUR_SRC}/kpoints/ordstar.f
${FLEUR_SRC}/kpoints/fulstar.f
${FLEUR_SRC}/kpoints/tetcon.f
${FLEUR_SRC}/kpoints/kvecon.f
${FLEUR_SRC}/kpoints/gen_bz.F90
${FLEUR_SRC}/kpoints/kprep.f
${FLEUR_SRC}/math/util.F
......
......@@ -2,6 +2,8 @@
MODULE m_make_spacegroup
USE m_juDFT
PRIVATE
public make_spacegroup
!********************************************************************
! Generate the spacegroup given the cell and the atomic positions/ids
! Takes into account the symmetry reductions needed for films,soc&SSDW
......@@ -13,7 +15,6 @@ CONTAINS
SUBROUTINE make_spacegroup(film,noco,cell,pos,atomid,sym)
USE m_bravaissymm
USE m_closure, ONLY : close_pt, closure
USE m_supercheck
USE m_types_noco
USE m_types_cell
......@@ -397,7 +398,7 @@ CONTAINS
!---> check closure of group
CALL closure(mops,mmrot,ttau,nops,index_op, lclose)
lclose=sym%closure()
IF ( ( ns==1 ) .AND. ( .NOT. lclose ) ) THEN
WRITE (6,'(/," Congratulations, you have found a system (not"," a supercell) that breaks the algorithms. Sorry...")')
......@@ -528,5 +529,61 @@ CONTAINS
END FUNCTION l_rotmatch
SUBROUTINE close_pt(nops,mrot,mtable)
IMPLICIT NONE
INTEGER, INTENT (IN) :: nops,mrot(3,3,nops)
INTEGER, INTENT (OUT) :: mtable(nops,nops) ! table(i,j) = {R_i|0}{R_j|0}
INTEGER :: i,j,k,mp(3,3),map(nops)
! loop over all operations
DO j = 1, nops
map(1:nops) = 0
! multiply {R_j|0}{R_i|0}
DO i = 1, nops
mp = matmul( mrot(:,:,j) , mrot(:,:,i) )
! determine which operation this is
DO k = 1, nops
IF ( all( mp(:,:)==mrot(:,:,k) ) ) THEN
IF ( map(i) .eq. 0 ) THEN
map(i) = k
ELSE
WRITE (6,'(" Symmetry error : multiple ops")')
CALL juDFT_error("close_pt: Multiple ops (Bravais)",calledby ="closure")
END IF
END IF
END DO
IF (map(i).eq.0) THEN
WRITE(6,*) 'Symmetry operations:'
DO k = 1, nops
WRITE(6,*) 'Matrix ', k, ':'
WRITE(6,'(3i7)') mrot(:,1,k)
WRITE(6,'(3i7)') mrot(:,2,k)
WRITE(6,'(3i7)') mrot(:,3,k)
WRITE(6,*) ''
END DO
WRITE (6,'(" Group not closed (Bravais lattice)")')
WRITE (6,'(" operation j=",i2," map=",12i4,:/,(21x,12i4))') j, map(1:nops)
WRITE(6,*) ''
WRITE(6,*) 'Expected product of operations ', j, ' and ', i, ':'
WRITE(6,'(3i7)') mp(:,1)
WRITE(6,'(3i7)') mp(:,2)
WRITE(6,'(3i7)') mp(:,3)
WRITE(6,*) ''
CALL juDFT_error("close_pt:Not closed",calledby="closure")
END IF
END DO
mtable(j,1:nops) = map(1:nops)
END DO
END SUBROUTINE close_pt
END SUBROUTINE make_spacegroup
END MODULE m_make_spacegroup
This diff is collapsed.
......@@ -79,10 +79,9 @@ SUBROUTINE writeOutParameters(mpi,input,sym,stars,atoms,vacuum,kpts,&
CALL closeXMLElement('volumes')
sumWeight = SUM(kpts%wtkpt(:kpts%nkpt))
WRITE(attributes(1),'(f0.8)') kpts%posScale
WRITE(attributes(2),'(f0.8)') sumWeight
WRITE(attributes(3),'(i0)') kpts%nkpt
CALL openXMLElementFormPoly('kPointList',(/'posScale ', 'weightScale', 'count '/),&
WRITE(attributes(1),'(f0.8)') sumWeight
WRITE(attributes(2),'(i0)') kpts%nkpt
CALL openXMLElementFormPoly('kPointList',(/ 'weightScale', 'count '/),&
attributes(:3),reshape((/8,11,5,10,10,5/),(/3,2/)))
DO i = 1, kpts%nkpt
WRITE(attributes(1),'(f12.6)') kpts%wtkpt(i)
......
......@@ -15,9 +15,23 @@ MODULE m_xmlIntWrapFort
USE m_juDFT
PRIVATE :: init_from_command_line
CONTAINS
private
TYPE t_xml
integer:: id
contains
PROCEDURE,NOPASS :: InitInterface
PROCEDURE,NOPASS :: ParseSchema
PROCEDURE,NOPASS :: ParseDoc
PROCEDURE,NOPASS :: ValidateDoc
PROCEDURE,NOPASS :: InitXPath
PROCEDURE,NOPASS :: GetNumberOfNodes
PROCEDURE,NOPASS :: SetAttributeValue
PROCEDURE,NOPASS :: FreeResources
end type
public t_xml
CONTAINS
SUBROUTINE init_from_command_line()
IMPLICIT NONE
......@@ -31,7 +45,7 @@ SUBROUTINE init_from_command_line()
ii=INDEX(xpath,":")
IF (ii==0) ii=LEN(TRIM(xpath))+1
IF (i>100.OR.ii-i>100) CALL judft_error("Too long xmlXPath argument",calledby="xmlIntWarpFort.f90")
CALL xmlSetAttributeValue(xpath(:i-1),xpath(i+1:ii-1))
CALL SetAttributeValue(xpath(:i-1),xpath(i+1:ii-1))
WRITE(*,*) "Set from command line:",TRIM(xpath(:i-1)),"=",TRIM(xpath(i+1:ii-1))
IF (ii+1<len(xpath))THEN
xpath=xpath(ii+1:)
......@@ -43,7 +57,7 @@ SUBROUTINE init_from_command_line()
END SUBROUTINE init_from_command_line
SUBROUTINE xmlInitInterface()
SUBROUTINE InitInterface()
USE iso_c_binding
USE m_types
......@@ -65,9 +79,9 @@ SUBROUTINE xmlInitInterface()
CALL juDFT_error("Could not initialize XML interface.",calledby="xmlInitInterface")
END IF
END SUBROUTINE xmlInitInterface
END SUBROUTINE InitInterface
SUBROUTINE xmlParseSchema(schemaFilename)
SUBROUTINE ParseSchema(schemaFilename)
USE iso_c_binding
USE m_types
......@@ -92,9 +106,9 @@ SUBROUTINE xmlParseSchema(schemaFilename)
CALL juDFT_error("XML Schema file not parsable: "//TRIM(ADJUSTL(schemaFilename)),calledby="xmlParseSchema")
END IF
END SUBROUTINE xmlParseSchema
END SUBROUTINE ParseSchema
SUBROUTINE xmlParseDoc(docFilename)
SUBROUTINE ParseDoc(docFilename)
USE iso_c_binding
USE m_types
......@@ -119,9 +133,9 @@ SUBROUTINE xmlParseDoc(docFilename)
CALL juDFT_error("XML document file not parsable: "//TRIM(ADJUSTL(docFilename)),calledby="xmlParseDoc")
END IF
END SUBROUTINE xmlParseDoc
END SUBROUTINE ParseDoc
SUBROUTINE xmlValidateDoc()
SUBROUTINE ValidateDoc()
USE iso_c_binding
USE m_types
......@@ -143,9 +157,9 @@ SUBROUTINE xmlValidateDoc()
CALL juDFT_error("XML document cannot be validated against Schema.",calledby="xmlValidateDoc")
END IF
END SUBROUTINE xmlValidateDoc
END SUBROUTINE ValidateDoc
SUBROUTINE xmlInitXPath()
SUBROUTINE InitXPath()
USE iso_c_binding
USE m_types
......@@ -164,19 +178,19 @@ SUBROUTINE xmlInitXPath()
errorStatus = 0
errorStatus = initializeXPath()
IF(errorStatus.NE.0) THEN
CALL juDFT_error("Could not initialize XPath.",calledby="xmlInitXPath")
CALL juDFT_error("Could not initialize XPath.",calledby="InitXPath")
END IF
CALL init_from_command_line()
END SUBROUTINE xmlInitXPath
END SUBROUTINE InitXPath
FUNCTION xmlGetNumberOfNodes(xPath)
FUNCTION GetNumberOfNodes(xPath)
USE iso_c_binding
USE m_types
IMPLICIT NONE
INTEGER :: xmlGetNumberOfNodes
INTEGER :: GetNumberOfNodes
CHARACTER(LEN=*,KIND=c_char), INTENT(IN) :: xPath
interface
......@@ -187,18 +201,18 @@ FUNCTION xmlGetNumberOfNodes(xPath)
end function getNumberOfXMLNodes
end interface
xmlGetNumberOfNodes = getNumberOfXMLNodes(TRIM(ADJUSTL(xPath))//C_NULL_CHAR)
GetNumberOfNodes = getNumberOfXMLNodes(TRIM(ADJUSTL(xPath))//C_NULL_CHAR)
END FUNCTION xmlGetNumberOfNodes
END FUNCTION GetNumberOfNodes
FUNCTION xmlGetAttributeValue(xPath)
FUNCTION GetAttributeValue(xPath)
USE iso_c_binding
USE m_types
IMPLICIT NONE
CHARACTER(LEN=255) :: xmlGetAttributeValue
CHARACTER(LEN=:),ALLOCATABLE :: GetAttributeValue
CHARACTER(LEN=*, KIND=c_char), INTENT(IN) :: xPath
......@@ -234,12 +248,12 @@ FUNCTION xmlGetAttributeValue(xPath)
END DO
length = i-1
xmlGetAttributeValue = value(1:length)
GetAttributeValue = trim(adjustl(value(1:length)))
END FUNCTION xmlGetAttributeValue
END FUNCTION GetAttributeValue
SUBROUTINE xmlSetAttributeValue(xPath,VALUE)
SUBROUTINE SetAttributeValue(xPath,VALUE)
USE iso_c_binding
USE m_types
......@@ -268,9 +282,9 @@ SUBROUTINE xmlSetAttributeValue(xPath,VALUE)
CALL juDFT_error("Attribute value could not be set.",calledby="xmlSetAttributeValue")
END IF
END SUBROUTINE xmlSetAttributeValue
END SUBROUTINE SetAttributeValue
SUBROUTINE xmlFreeResources()
SUBROUTINE FreeResources()
USE iso_c_binding
USE m_types
......@@ -293,6 +307,6 @@ SUBROUTINE xmlFreeResources()
STOP 'Error!'
END IF
END SUBROUTINE xmlFreeResources
END SUBROUTINE FreeResources
END MODULE m_xmlIntWrapFort
......@@ -15,9 +15,6 @@ kpoints/tetcon.f
set(fleur_F90 ${fleur_F90}
kpoints/gkptwgt.f90
kpoints/julia.f90
kpoints/kpoints.f90
kpoints/unfoldBandKPTS.f90
kpoints/gen_bz.F90
kpoints/brzone2.f90
)
......@@ -28,7 +28,6 @@ CONTAINS
END DO
CALL inv3(p_cell%amat,p_cell%bmat,p_cell%omtil)
p_cell%bmat=p_cell%bmat*tpi_const
p_cell%latnam=cell%latnam
END SUBROUTINE build_primitive_cell
SUBROUTINE unfold_band_kpts(banddos,p_cell,cell,p_kpts,kpts)
......
......@@ -94,7 +94,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
LOGICAL :: l_error
CALL regCharges%init(input,atoms)
CALL dos%init(input,atoms,dimension,kpts,vacuum)
CALL dos%init(dimension%neigd,input,atoms,kpts,vacuum)
CALL moments%init(input,atoms)
CALL mcd%init1(banddos,dimension,input,atoms,kpts)
CALL slab%init(banddos,dimension,atoms,cell,input,kpts)
......
......@@ -18,7 +18,7 @@
USE m_postprocessInput
USE m_gen_map
USE m_dwigner
USE m_gen_bz
!USE m_gen_bz
USE m_ylm
USE m_InitParallelProcesses
USE m_xmlOutput
......@@ -32,7 +32,7 @@
USE m_checks
USE m_prpqfftmap
USE m_writeOutHeader
USE m_fleur_init_old
!USE m_fleur_init_old
USE m_types_xcpot_inbuild
USE m_mpi_bc_xcpot
......@@ -197,10 +197,8 @@
numSpecies = SIZE(speciesRepAtomType)
CALL w_inpXML(&
atoms,vacuum,input,stars,sliceplot,forcetheo,banddos,&
cell,sym,xcpot,noco,oneD,hybrid,kpts,kpts%nkpt3,kpts%l_gamma,&
namex,relcor,a1,a2,a3,dtild,input%comment,&
.TRUE.,filename,&
.TRUE.,enpara)
cell,sym,xcpot,noco,oneD,hybrid,kpts,enpara,&
.true.,[.true.,.true.,.true.,.true.])
DEALLOCATE(atomTypeSpecies,speciesRepAtomType)
DEALLOCATE(xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs)
......@@ -216,10 +214,10 @@
#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,enpara,xcpot,kpts,hybrid,&
oneD,coreSpecInput,l_opti)
!CALL fleur_init_old(mpi,&
! input,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,&
! sliceplot,banddos,enpara,xcpot,kpts,hybrid,&
! oneD,coreSpecInput,l_opti)
END IF ! end of else branch of "IF (input%l_inpXML) THEN"
!
......@@ -458,7 +456,7 @@
ELSE
IF ( banddos%dos .AND. banddos%ndir == -3 ) THEN
WRITE(*,*) 'Recalculating k point grid to cover the full BZ.'
CALL gen_bz(kpts,sym)
!CALL gen_bz(kpts,sym)
kpts%nkpt = kpts%nkptf
DEALLOCATE(kpts%bk,kpts%wtkpt)
ALLOCATE(kpts%bk(3,kpts%nkptf),kpts%wtkpt(kpts%nkptf))
......
......@@ -252,7 +252,7 @@ SUBROUTINE rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars
ALLOCATE(dEdOcc(MAXVAL(results%neig(1:kpts%nkpt,1:input%jspins)),kpts%nkpt,input%jspins))
CALL regCharges%init(input,atoms)
CALL dos%init(input,atoms,dimension,kpts,vacuum)
CALL dos%init(dimension%neigd,input,atoms,kpts,vacuum)
CALL moments%init(input,atoms)
CALL overallDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
CALL overallVCoul%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTCOUL)
......
......@@ -5,6 +5,7 @@ types/types.F90
types/types_mat.F90
types/types_xcpot.F90
types/types_xcpot_inbuild.F90
types/types_xcpot_inbuild_nofunction.F90
types/types_xcpot_data.F90
types/types_xcpot_libxc.F90
types/types_mpi.F90
......@@ -37,6 +38,9 @@ types/types_oneD.f90
types/types_hybrid.f90
types/types_noco.f90
types/types_banddos.f90
types/types_dimension.f90
types/types_vacuum.f90
)
set(inpgen_F90 ${inpgen_F90}
......
......@@ -101,9 +101,8 @@ MODULE m_types_atoms
TYPE(t_utype), ALLOCATABLE::lda_u(:)
INTEGER, ALLOCATABLE :: relax(:, :) !<(3,ntype)
INTEGER, ALLOCATABLE :: nflip(:) !<flip magnetisation of this atom
! CONTAINS
! procedure :: nsp => calc_nsp_atom
contains
CONTAINS
procedure :: nsp => calc_nsp_atom
PROCEDURE :: same_species
END TYPE t_atoms
......@@ -130,5 +129,11 @@ MODULE m_types_atoms
same_species=same_species.and.trim(atoms%econf(n)%valenceconfig)==trim(atoms%econf(nn)%valenceconfig)
same_species=same_species.and.trim(atoms%econf(n)%valenceconfig)==trim(atoms%econf(nn)%valenceconfig)
end function
pure function calc_nsp_atom(self) result(nsp)
implicit none
CLASS(t_atoms),INTENT(IN) :: self
INTEGER :: nsp
nsp = (self%lmaxd+1+MOD(self%lmaxd+1,2))*(2*self%lmaxd+1)
end function
END MODULE m_types_atoms
......@@ -55,7 +55,7 @@ CONTAINS
this%firstloop=.FALSE.
END FUNCTION forcetheo_next_job
FUNCTION forcetheo_eval(this,eig_id,atoms,kpts,sym,&
FUNCTION forcetheo_eval(this,eig_id,DIMENSION,atoms,kpts,sym,&
cell,noco, input,mpi, oneD,enpara,v,results)RESULT(skip)
USE m_types_atoms
USE m_types_oneD
......@@ -68,12 +68,14 @@ CONTAINS
USE m_types_misc
USE m_types_kpts
USE m_types_enpara
use m_types_dimension
IMPLICIT NONE
CLASS(t_forcetheo),INTENT(INOUT):: this
LOGICAL :: skip
!Stuff that might be used...
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_dimension),INTENT(IN) :: dimension
TYPE(t_oneD),INTENT(IN) :: oneD
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
......
......@@ -164,12 +164,4 @@ MODULE m_types_setup
!---> gwf
END TYPE t_wann
CONTAINS
pure function calc_nsp_atom(self) result(nsp)
implicit none
CLASS(t_atoms),INTENT(IN) :: self
INTEGER :: nsp
nsp = (self%lmaxd+1+MOD(self%lmaxd+1,2))*(2*self%lmaxd+1)
end function
END MODULE m_types_setup
......@@ -5,6 +5,7 @@
!--------------------------------------------------------------------------------
MODULE m_types_sym
USE m_juDFT
IMPLICIT NONE
PRIVATE
!symmetry information
......@@ -56,9 +57,14 @@ MODULE m_types_sym
CONTAINS
PROCEDURE :: init
PROCEDURE :: print_xml
PROCEDURE :: closure
PROCEDURE,PRIVATE :: check_close
END TYPE t_sym
PUBLIC t_sym
CONTAINS
SUBROUTINE print_xml(sym,fh,filename)
CLASS(t_sym),INTENT(IN) :: sym
INTEGER,INTENT(in) ::fh
......@@ -86,7 +92,6 @@ CONTAINS
SUBROUTINE init(sym,cell,film)
!Generates missing symmetry info.
!tau,mrot and nop have to be specified already
USE m_closure
USE m_dwigner
USE m_types_cell
CLASS(t_sym),INTENT(INOUT):: sym
......@@ -106,7 +111,7 @@ CONTAINS
IF (ALLOCATED(sym%invtab)) DEALLOCATE(sym%invtab)
IF (ALLOCATED(sym%multab)) DEALLOCATE(sym%multab)
ALLOCATE ( sym%invtab(sym%nop),sym%multab(sym%nop,sym%nop) )
CALL check_close(sym%nop,sym%mrot,sym%tau, sym%multab,sym%invtab,optype)
CALL sym%check_close(optype)
!---> determine properties of symmmetry operations,
! Code previously in symproperties
......@@ -234,6 +239,116 @@ CONTAINS
CALL d_wigner(sym%nop,sym%mrot,cell%bmat,3,sym%d_wgn)
!---> redo to ensure proper mult. table and mapping functions
CALL check_close(sym%nop,sym%mrot,sym%tau, sym%multab,sym%invtab,optype)
CALL sym%check_close(optype)
END SUBROUTINE init
FUNCTION closure(sym)RESULT(lclose)
CLASS(t_sym),INTENT(IN):: sym
LOGICAL :: lclose
!INTEGER, INTENT (IN) :: mops ! number of operations of the bravais lattice
!INTEGER, INTENT (IN) :: sym%nop ! number of operations in space group
!INTEGER, INTENT (IN) :: sym%mrot(3,3,mops) ! refer to the operations of the
!REAL, INTENT (IN) :: sym%tau(3,mops) ! bravais lattice
!LOGICAL, INTENT (OUT) :: lclose
REAL ttau(3),eps7
INTEGER i,j,k,mp(3,3),map(sym%nop)
eps7 = 1.0e-7
! loop over all operations
DO j = 1, sym%nop
map(1:sym%nop) = 0
! multiply {R_j|t_j}{R_i|t_i}
DO i = 1, sym%nop
mp = matmul( sym%mrot(:,:,j) , sym%mrot(:,:,i) )
ttau = sym%tau(:,j) + matmul( sym%mrot(:,:,j) , sym%tau(:,i) )
ttau = ttau - anint( ttau - eps7 )
! determine which operation this is
DO k=1,sym%nop
IF ( all( mp(:,:) == sym%mrot(:,:,k) ) .AND. all( abs( ttau(:)-sym%tau(:,k) ) < eps7 ) ) THEN
IF ( map(i) .eq. 0 ) THEN
map(i) = k
ELSE
write(6,*)'ERROR Closure: Multiplying ', j,' with ',k, ' and with ',map(i)
write(6,*) 'yields the same matrix'
lclose = .false.
RETURN
END IF
END IF
END DO
IF (map(i).eq.0) THEN
write(6,*)'ERROR Closure:',i,' times',j,' leaves group'
lclose = .false.
RETURN
END IF
END DO
END DO
lclose = .true.
END FUNCTION closure
SUBROUTINE check_close(sym,optype)
CLASS(t_sym),INTENT(inout)::sym
INTEGER, INTENT (OUT) :: optype(sym%nop)
REAL ttau(3)
INTEGER i,j,n,k,mp(3,3),mdet,mtr
REAL, PARAMETER :: eps=1.0e-7
INTEGER, PARAMETER :: cops(-1:3)=(/ 2, 3, 4, 6, 1 /)
sym%invtab(1:sym%nop) = 0
sym%multab = 0
! loop over all operations
DO j = 1, sym%nop
! multiply {R_j|t_j}{R_i|t_i}
DO i = 1, sym%nop
mp = matmul( sym%mrot(:,:,j) , sym%mrot(:,:,i) )
ttau = sym%tau(:,j) + matmul( sym%mrot(:,:,j) , sym%tau(:,i) )
ttau = ttau - anint( ttau - eps )
! determine which operation this is
DO k=1,sym%nop
IF ( all( mp(:,:) == sym%mrot(:,:,k) ) .and. all( abs( ttau(:)-sym%tau(:,k) ) < eps ) ) THEN
IF ( sym%multab(j,i) .eq. 0 ) THEN
sym%multab(j,i) = k
IF (k .eq. 1) sym%invtab(j)=i
ELSE
WRITE(6,'(" Symmetry error: multiple ops")')
CALL juDFT_error("check_close: Multiple ops",calledby ="closure")
END IF
END IF
END DO
IF (sym%multab(j,i).eq.0) THEN
WRITE (6,'(" Group not closed")')
WRITE (6,'(" j , i =",2i4)') j,i
CALL juDFT_error("check_close: Not closed",calledby="closure")
END IF
END DO
END DO
! determine the type of each operation
DO n = 1, sym%nop
mtr = sym%mrot(1,1,n) + sym%mrot(2,2,n) + sym%mrot(3,3,n)
mdet = sym%mrot(1,1,n)*(sym%mrot(2,2,n)*sym%mrot(3,3,n)-sym%mrot(3,2,n)*sym%mrot(2,3,n)) +&
sym%mrot(1,2,n)*(sym%mrot(3,1,n)*sym%mrot(2,3,n)-sym%mrot(2,1,n)*sym%mrot(3,3,n)) +&
sym%mrot(1,3,n)*(sym%mrot(2,1,n)*sym%mrot(3,2,n)-sym%mrot(3,1,n)*sym%mrot(2,2,n))
optype(n) = mdet*cops(mdet*mtr)
END DO
END SUBROUTINE check_close
END MODULE m_types_sym
......@@ -95,7 +95,7 @@ CONTAINS
CALL regCharges%init(input, atoms)
CALL dos%init(input, atoms, DIMENSION, kpts, vacuum)
CALL dos%init(DIMENSION%neigd,input,atoms,kpts, vacuum)
CALL moments%init(input, atoms)
tmp_results = results
......