Commit a7ce652d authored by Gustav Bihlmayer's avatar Gustav Bihlmayer

Corrected the treatment of centered monoclinic lattices in lattice2.f.

To transform the atom positions in atom_sym.f into the centered basis,
a variable (i_c) is passed through to indicate the type of centering.
parent 2b763377
......@@ -5,7 +5,7 @@
> dispfh,outfh,errfh,dispfn,natmax,
X natin,atomid,atompos,
X ngen,mmrot,ttr,
> cartesian,symor,as,bs,nop48,
> cartesian,i_c,symor,as,bs,nop48,
< ntype,nat,nops,mrot,tau,
< neq,ntyrep,zatom,natype,natrep,natmap,pos)
......@@ -53,6 +53,7 @@
LOGICAL, INTENT (IN) :: symor ! whether to reduce to symmorphic subgroup
INTEGER, INTENT (INOUT) :: natin ! formerly 'ntype0'
INTEGER, INTENT (IN) :: ngen ! Number of generators
INTEGER, INTENT (IN) :: i_c ! centering of lattice
INTEGER, INTENT (IN) :: nop48, natmax ! dimensioning
INTEGER, INTENT (INOUT) :: mmrot(3,3,nop48)
REAL, INTENT (INOUT) :: ttr(3,nop48)
......@@ -76,7 +77,7 @@
REAL tr(3),tt(3),disp(3,natmax)
INTEGER mp(3,3),mtmp(3,3)
REAL ttau(3),orth(3,3)
REAL ttau(3),orth(3,3),tc(3,3),td(3,3)
INTEGER mmrot2(3,3,ngen)
REAL ttr2(3,ngen)
......@@ -84,20 +85,49 @@
LOGICAL l_exist,lclose,l_inipos
INTEGER n,na,ng,ncyl,nc,no,nop0,nn,nt,i,j,mops
INTEGER ios,istep0
REAL eps7
CHARACTER(len=30) :: filen
REAL, ALLOCATABLE :: inipos(:,:)
eps7= 1.0e-7 ; istep0 = 0
REAL, PARAMETER :: eps = 1.0e-7, isqrt3 = 1.0/sqrt(3.0),
& thrd = 1.0/3.0, mtthrd = -2.0/3.0
REAL :: lmat(3,3,8)
DATA lmat / 1.0, 0.0, 0.0, ! 1: primitive : P
& 0.0, 1.0, 0.0,
& 0.0, 0.0, 1.0,
+ -1.0, 1.0, 1.0, ! 2: Inverse (F)
& 1.0, -1.0, 1.0,
& 1.0, 1.0, -1.0,
+ 0.0, 1.0, 1.0, ! 3: Inverse (I)
& 1.0, 0.0, 1.0,
& 1.0, 1.0, 0.0,
+ 1.0, 1.0, 0.0, ! 4: Inverse (hP)
& -isqrt3, isqrt3, 0.0,
& 0.0, 0.0, 1.0,
+ 0.0, isqrt3, -isqrt3, ! 5: Inverse (hR)
& mtthrd, thrd, thrd,
& thrd, thrd, thrd,
+ 1.0, 1.0, 0.0, ! 6: Inverse ( S (C) )
& -1.0, 1.0, 0.0,
& 0.0, 0.0, 1.0,
+ 1.0, 0.0, 1.0, ! 7: Inverse (B)
& 0.0, 1.0, 0.0,
& -1.0, 0.0, 1.0,
+ 1.0, 0.0, 0.0, ! 8: Inverse (A)
& 0.0, 1.0, -1.0,
& 0.0, 1.0, 1.0/
istep0 = 0
!
!---> take atomic positions and shift to (-1/2,1/2] in lattice coords.
!
natin = abs(natin)
DO n=1,natin
IF (cartesian) THEN ! convert to lattice coords. if necessary
atompos(:,n) = matmul( bs, atompos(:,n) )
! atompos(:,n) = matmul( bs, atompos(:,n) )
atompos(:,n) = matmul( lmat(:,:,i_c), atompos(:,n) )
ENDIF
atompos(:,n) = atompos(:,n) - anint( atompos(:,n) - eps7 )
atompos(:,n) = atompos(:,n) - anint( atompos(:,n) - eps )
ENDDO
!---> store the positions (in lattice coord.s) given in the input file
......@@ -124,7 +154,7 @@
tr = matmul( bs, tr )
ENDIF
atompos(:,n) = atompos(:,n) + tr(:)
atompos(:,n) = atompos(:,n) - anint( atompos(:,n)- eps7 )
atompos(:,n) = atompos(:,n) - anint( atompos(:,n)- eps )
ENDDO
CLOSE (dispfh)
IF ( ios==0 ) THEN
......@@ -142,10 +172,16 @@
!---> save generators
IF (cartesian) THEN ! convert to lattice coords. if necessary
DO ng = 2, ngen+1
mmrot2(:,:,1) = matmul( bs, mmrot(:,:,ng) )
mmrot(:,:,ng) = matmul( mmrot2(:,:,1), as )
ttr2(:,1) = matmul( bs, ttr(:,ng) )
! mmrot2(:,:,1) = matmul( bs, mmrot(:,:,ng) )
! mmrot(:,:,ng) = matmul( mmrot2(:,:,1), as )
tc = mmrot(:,:,ng)
td = matmul( bs, tc )
tc = matmul( td, as )
mmrot(:,:,ng) = NINT(tc)
write(*,*) i_c, ttr(:,ng)
ttr2(:,1) = matmul( lmat(:,:,i_c), ttr(:,ng) )
ttr(:,ng) = ttr2(:,1)
write(*,*) mmrot(:,:,ng),ttr(:,ng)
ENDDO
ENDIF
mmrot2(:,:,1:ngen) = mmrot(:,:,2:ngen+1)
......@@ -214,7 +250,7 @@
ENDIF
!---> rewrite all the non-primitive translations so in (-1/2,1/2]
ttr(:,1:nops) = ttr(:,1:nops) - anint( ttr(:,1:nops)- eps7 )
ttr(:,1:nops) = ttr(:,1:nops) - anint( ttr(:,1:nops)- eps )
!---> allocate arrays for space group information (mod_spgsym)
! if( nopd < nops )then
......@@ -226,6 +262,7 @@
mrot(:,:,n) = mmrot(:,:,n)
tau(:,n) = ttr(:,n)
index_op(n) = n
write(*,*) n,mrot(:,:,n),tau(:,n)
ENDDO
!---> check that the group is closed, etc.
......@@ -245,7 +282,7 @@
!---> reduce symmetry to the largest symmorphic subgroup
j = 1
DO i = 1, nops
IF ( all ( abs( tau(:,i) ) < eps7 ) ) THEN
IF ( all ( abs( tau(:,i) ) < eps ) ) THEN
IF ( j<i ) then
mrot(:,:,j) = mrot(:,:,i)
ENDIF
......@@ -276,10 +313,10 @@
DO n = 1, nat
tt = ( atompos(:,nt) - tpos(:,n) )
& - anint( atompos(:,nt) - tpos(:,n) )
IF ( all( abs(tt) < eps7 ) ) THEN
IF ( all( abs(tt) < eps ) ) THEN
icount(n) = icount(n) + 1
imap(nt) = n
IF ( abs( atomid(nt)-atomid(ity(n)) ) < eps7 )
IF ( abs( atomid(nt)-atomid(ity(n)) ) < eps )
& CYCLE repres_atoms
CALL juDFT_error("ERROR! mismatch between atoms."
+ ,calledby ="atom_sym")
......@@ -297,14 +334,14 @@
!---> loop over operations
opts: DO no = 2, nops
tr = matmul( mrot(:,:,no) , atompos(:,nt) ) + tau(:,no)
tr = tr - anint( tr - eps7 )
tr = tr - anint( tr - eps )
!---> check whether this is a new atom
DO n = 1, nneq(ntype)
tt = ( tr-tpos(:,nat+n) ) - anint( tr-tpos(:,nat+n) )
IF ( all( abs(tt) < eps7 ) ) THEN
IF ( all( abs(tt) < eps ) ) THEN
nn = ity(nat+n)
IF ( abs( atomid(nt)-atomid(nn) ) < eps7 ) CYCLE opts
IF ( abs( atomid(nt)-atomid(nn) ) < eps ) CYCLE opts
WRITE (6,'(" Mismatch between atoms and",
& " symmetry input")')
CALL juDFT_error("atom_sym: mismatch rotated",calledby
......
......@@ -8,7 +8,7 @@
> dbgfh,errfh,outfh,dispfh,dispfn,
> cal_symm,cartesian,symor,oldfleur,
> natin,natmax,nop48,
> atomid,atompos,a1,a2,a3,aa,scale,noangles,
> atomid,atompos,a1,a2,a3,aa,scale,noangles,i_c,
< invs,zrfs,invs2,nop,nop2,
< ngen,mmrot,ttr,ntype,nat,nops,
< neq,ntyrep,zatom,natype,natrep,natmap,
......@@ -26,7 +26,7 @@
!===> Arguments
LOGICAL, INTENT(IN) :: cal_symm,cartesian,oldfleur,noangles
INTEGER, INTENT(IN) :: ngen,natmax,nop48
INTEGER, INTENT(IN) :: ngen,natmax,nop48,i_c
INTEGER, INTENT(IN) :: dbgfh,errfh,outfh,dispfh ! file handles, mainly 6
REAL, INTENT(IN) :: aa
LOGICAL, INTENT(INOUT) :: symor ! on input: if true, reduce symmetry if oldfleur
......@@ -154,7 +154,7 @@
> dispfh,outfh,errfh,dispfn,natmax,
X natin,atomid,atompos,
X ngen,mmrot,ttr,
> cartesian,symor,as,bs,nop48,
> cartesian,i_c,symor,as,bs,nop48,
< ntype,nat,nops,mrot,tau,
< neq,ntyrep,zatom,natype,natrep,natmap,pos)
......
......@@ -20,7 +20,7 @@ PROGRAM inpgen
IMPLICIT NONE
INTEGER natmax,nop48,nline,natin,ngen,i,j
INTEGER nops,no3,no2,na,numSpecies
INTEGER nops,no3,no2,na,numSpecies,i_c
INTEGER infh,errfh,bfh,warnfh,symfh,dbgfh,outfh,dispfh
LOGICAL cal_symm,checkinp,newSpecies,noangles
LOGICAL cartesian,oldfleur,l_hyb ,inistop
......@@ -73,7 +73,7 @@ PROGRAM inpgen
& natmax,nop48,&
& nline,xl_buffer,buffer,&
& title,input%film,cal_symm,checkinp,sym%symor,&
& cartesian,oldfleur,a1,a2,a3,vacuum%dvac,aa,scale,noangles,&
& cartesian,oldfleur,a1,a2,a3,vacuum%dvac,aa,scale,noangles,i_c,&
& factor,natin,atomid,atompos,ngen,mmrot,ttr,&
& l_hyb,noco%l_soc,noco%l_ss,noco%theta,noco%phi,noco%qss,inistop)!keep
......@@ -113,7 +113,7 @@ PROGRAM inpgen
& dbgfh,errfh,outfh,dispfh,dispfn,&
& cal_symm,cartesian,sym%symor,input%film,&
& natin,natmax,nop48,&
& atomid,atompos,a1,a2,a3,aa,scale,noangles,&
& atomid,atompos,a1,a2,a3,aa,scale,noangles,i_c,&
& sym%invs,sym%zrfs,sym%invs2,sym%nop,sym%nop2,&
& ngen,mmrot,ttr,atoms%ntype,atoms%nat,nops,&
& atoms%neq,ntyrep,atoms%zatom,natype,natrep,natmap,&
......
......@@ -8,7 +8,7 @@
CONTAINS
SUBROUTINE lattice2(
> buffer,xl_buffer,errfh,bfh,nline,
< a1,a2,a3,aa,scale,noangles,ios )
< a1,a2,a3,aa,scale,noangles,i_c,ios )
USE m_constants
IMPLICIT NONE
......@@ -20,7 +20,7 @@
REAL, INTENT (OUT) :: a1(3),a2(3),a3(3)
REAL, INTENT (OUT) :: aa
REAL, INTENT (OUT) :: scale(3)
INTEGER, INTENT (OUT) :: ios
INTEGER, INTENT (OUT) :: i_c,ios
LOGICAL, INTENT (OUT) :: noangles
!==> Local Variables
......@@ -63,8 +63,8 @@
& 0.0, 1.0, 0.0,
& 0.5, 0.0, 0.5,
+ 1.0, 0.0, 0.0, ! 8: base-centered: A
& 0.0, 0.5, -0.5,
& 0.0, 0.5, 0.5/
& 0.0, 0.5, 0.5,
& 0.0, -0.5, 0.5/
!===> 12: monoclinic-P (mP)
!===> 13: monoclinic-P (mS) (mA) (mB) (mC)
......@@ -76,6 +76,7 @@
latsys = ' ' ; a0 = 0.0
a = 0.0 ; b = 0.0 ; c = 0.0
alpha = 0.0 ; beta = 0.0 ; gamma = 0.0
scale = 0.0
READ (bfh,lattice,err=911,end=911,iostat=ios)
......@@ -105,10 +106,10 @@
& latsys =='simple-cubic' ) THEN
noangles=.true.
i = 1
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 1
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
IF ( a.NE.b .OR. a.NE.c ) err = 11
IF ( ar.NE.br .OR. ar.NE.cr .OR. ar.NE.(pi_const/2.0) ) err = 12
......@@ -119,10 +120,10 @@
& latsys =='face-centered-cubic' ) THEN
noangles=.true.
i = 2
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 2
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
IF ( a.NE.b .OR. a.NE.c ) err = 21
......@@ -132,10 +133,10 @@
& latsys =='body-centered-cubic' ) THEN
noangles=.true.
i = 3
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 3
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
IF ( a.NE.b .OR. a.NE.c ) err = 31
......@@ -145,10 +146,10 @@
& .OR.latsys =='hexagonal' ) THEN
noangles=.true.
i = 4
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 4
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
IF ( a.NE.b ) err = 41
......@@ -157,10 +158,11 @@
ELSEIF ( latsys =='hdp' ) THEN
noangles=.true.
i = 4
a1 = lmat((/2,1,3/),1,i)
a2 = -lmat((/2,1,3/),2,i)
a3 = lmat(:,3,i)
i_c = 4
a1 = lmat((/2,1,3/),1,i_c)
a2 = -lmat((/2,1,3/),2,i_c)
a3 = lmat(:,3,i_c)
i_c = 9
IF ( a.NE.b ) err = 41
......@@ -171,10 +173,10 @@
& latsys =='rho'.OR.latsys =='trigonal' ) THEN
noangles=.false.
i = 5
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 5
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
IF ( a.NE.b ) err = 51
IF ( alpha.EQ.0.0 .OR.
......@@ -196,6 +198,7 @@
& latsys =='trigonal2' ) THEN
noangles=.false.
i_c = 5
IF ( a.NE.b .OR. a.NE.c ) err = 53
IF ( alpha.NE.beta .OR. alpha.NE.gamma ) err = 54
......@@ -223,10 +226,10 @@
& .OR.latsys =='simple-tetragonal' ) THEN
noangles=.true.
i = 1
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 1
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
IF ( a.NE.b ) err = 61
IF ( ar.NE.br .OR. ar.NE.cr .OR. ar.NE.(pi_const/2.0) ) err= 62
......@@ -237,10 +240,10 @@
& .OR.latsys =='body-centered-tetragonal' ) THEN
noangles=.true.
i = 3
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 3
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
IF ( a.NE.b ) err = 61
......@@ -250,10 +253,10 @@
& latsys =='simple-orthorhombic' ) THEN
noangles=.true.
i = 1
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 1
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
!===> 9: orthorhombic-F (oF)
......@@ -262,10 +265,10 @@
& latsys =='face-centered-orthorhombic' ) THEN
noangles=.true.
i = 2
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 2
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
!===> 10: orthorhombic-I (oI)
......@@ -274,10 +277,10 @@
& latsys =='body-centered-orthorhombic' ) THEN
noangles=.true.
i = 3
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 3
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
!===> 11: orthorhombic-S (oS) (oC)
......@@ -286,10 +289,10 @@
& latsys =='base-centered-orthorhombic' ) THEN
noangles=.true.
i = 6
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 6
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
!===> 11a: orthorhombic-A (oA)
......@@ -298,10 +301,10 @@
& latsys =='base-centered-orthorhombic2' ) THEN
noangles=.true.
i = 8
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 8
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
!===> 11b: orthorhombic-B (oB)
......@@ -310,10 +313,10 @@
& latsys =='base-centered-orthorhombic3' ) THEN
noangles=.true.
i = 7
a1 = lmat(:,1,i)
a2 = lmat(:,2,i)
a3 = lmat(:,3,i)
i_c = 7
a1 = lmat(:,1,i_c)
a2 = lmat(:,2,i_c)
a3 = lmat(:,3,i_c)
!===> 12: monoclinic-P (mP)
ELSEIF ( latsys =='monoclinic-P'.OR.latsys =='mP'.OR
......@@ -321,6 +324,8 @@
& latsys =='simple-monoclinic' ) THEN
noangles=.false.
i_c = 1
IF ( (abs(alpha-90.0)<eps).AND.(abs(beta-90.0)<eps) ) THEN
IF ( ABS(gamma - 90.0) <eps ) CALL juDFT_error
+ ("no monoclinic angle!",calledby ="lattice2")
......@@ -359,12 +364,16 @@
+ ,calledby ="lattice2")
ENDIF
CALL brvmat ( alpha, beta, gamma, am )
i = 8
am = matmul ( am, lmat(:,:,i) )
am(:,1) = a * am(:,1)
am(:,2) = b * am(:,2)
am(:,3) = c * am(:,3)
i_c = 8
am = matmul ( am, lmat(:,:,i_c) )
a1 = am(:,1)
a2 = am(:,2)
a3 = am(:,3)
CALL angles( am )
scale = 1.0
!===> 13b monoclinic-B (mB)
ELSEIF ( latsys =='monoclinic-B'.OR.latsys =='mB'.OR
......@@ -380,12 +389,40 @@
+ ,calledby ="lattice2")
ENDIF
CALL brvmat ( alpha, beta, gamma, am )
i = 7
am = matmul ( am, lmat(:,:,i) )
am(:,1) = a * am(:,1)
am(:,2) = b * am(:,2)
am(:,3) = c * am(:,3)
i_c = 7
am = matmul ( am, lmat(:,:,i_c) )
a1 = am(:,1)
a2 = am(:,2)
a3 = am(:,3)
CALL angles( am )
scale = 1.0
!===> 13c monoclinic-I (mI)
ELSEIF ( latsys =='monoclinic-I'.OR.latsys =='mI'.OR
& .latsys =='moI' ) THEN
noangles=.false.
IF ( (abs(alpha-90.0)<eps).AND.(abs(beta-90.0)<eps) ) THEN
IF ( ABS(gamma - 90.0) <eps ) CALL juDFT_error
+ ("no monoclinic angle!",calledby ="lattice2")
ELSE
CALL juDFT_error("Please take gamma as monoclinic angle!"
+ ,calledby ="lattice2")
ENDIF
CALL brvmat ( alpha, beta, gamma, am )
am(:,1) = a * am(:,1)
am(:,2) = b * am(:,2)
am(:,3) = c * am(:,3)
i_c = 3
am = matmul ( am, lmat(:,:,i_c) )
a1 = am(:,1)
a2 = am(:,2)
a3 = am(:,3)
CALL angles( am )
scale = 1.0
!===> 14: triclinic (aP)
......@@ -393,6 +430,7 @@
& latsys =='tcl' ) THEN
noangles=.false.
i_c = 1
CALL brvmat ( alpha, beta, gamma, am )
a1 = am(:,1)
a2 = am(:,2)
......@@ -421,9 +459,11 @@
WRITE (errfh,*)
ENDIF
IF (abs(scale(1)) < eps) THEN
scale(1) = a
scale(2) = b
scale(3) = c
ENDIF
aa = a0
911 CONTINUE
......
......@@ -9,7 +9,7 @@
> natmax,nop48,
X nline,xl_buffer,buffer,
< title,film,cal_symm,checkinp,symor,
< cartesian,oldfleur,a1,a2,a3,dvac,aa,scale,noangles,
< cartesian,oldfleur,a1,a2,a3,dvac,aa,scale,noangles,i_c,
< factor,natin,atomid,atompos,ngen,mmrot,ttr,
< l_hyb,l_soc,l_ss,theta,phi,qss,inistop)
......@@ -27,7 +27,7 @@
LOGICAL :: cal_symm, checkinp, symor, film
LOGICAL :: cartesian,oldfleur,inistop
LOGICAL, INTENT (OUT) :: l_hyb,l_soc,l_ss,noangles
INTEGER, INTENT (OUT) :: natin
INTEGER, INTENT (OUT) :: natin,i_c
INTEGER, INTENT (OUT) :: ngen
REAL, INTENT (OUT) :: aa,theta,phi
REAL, INTENT (OUT) :: dvac
......@@ -140,7 +140,7 @@
IF ( buffer(1:8)=='&lattice' ) THEN
CALL lattice2(
> buffer,xl_buffer,errfh,bfh,nline,
< a1,a2,a3,aa,scale,noangles,ios )
< a1,a2,a3,aa,scale,noangles,i_c,ios )
dvac = 0.00
IF ( ios.NE.0 ) THEN
WRITE (errfh,*)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment