diff --git a/source/voronoi/maindriver12.f b/source/voronoi/maindriver12.f
index 05ea5c1291c097be7c2a34ee9d1ad0970d7a56c8..f1e651c53dc3f5520f9247c5bcc75e0c57d6d379 100644
--- a/source/voronoi/maindriver12.f
+++ b/source/voronoi/maindriver12.f
@@ -119,6 +119,8 @@ c
+ ZATOM(NTOTD) ! Nuclear charge
INTEGER ICC, ! center of cluster for output of GF
+ ICLS, NAEZ, ! number of atoms in unit cell
+ + NVAC, ! number of empty cells
+ + NVAC_IT,
+ NATYP, ! number of kinds of atoms in unit cell
+ NCLS, ! number of reference clusters
+ NEMB, ! number of 'embedding' positions
@@ -145,6 +147,8 @@ c
& ROUT_ALL(NTOTD), ! Outer cell-radius per atom
& DISTNN(NAEZD+NIMPD), ! Distance from cell center to nearest-neighbor cell center (2*RMTHLF)
& VOLUME_ALL(NTOTD), ! Volume per atom
+ & VCENTER_ALL(3,NTOTD), ! Center of the voronoi cells
+ & VCENTER_SQSUM,
& A3_ALL(NFACED,NTOTD), ! A3,B3,C3,D3: Defining the faces per atom
& B3_ALL(NFACED,NTOTD),
& C3_ALL(NFACED,NTOTD),
@@ -226,7 +230,7 @@ c
& VOLUMECL(NSHAPED),RWSCL(NSHAPED),RMTCL(NSHAPED)
REAL*8 DX(NTOTD),DY(NTOTD),DZ(NTOTD)
REAL*8 ROUT,RTEST,DLT,CRAD,RX,RY,RZ,MTRADIUS,VTOT,
- & SHAPESHIFT(3,NTOTD)
+ & SHAPESHIFT(3,NTOTD),VCENTER(3)
CHARACTER*256 UIO
INTEGER NATOMS,NSITES,NSHAPE ! Number of atoms, sites, shapes
INTEGER LMAX,KEYPAN,NPOI,NA,IAT,JAT,ICL,N1A,I2,II,ISITE
@@ -269,7 +273,7 @@ c
c -----------------------------------------------------------------------
DATA BBOX/2.0d0,2.0d0,3.0d0/
DATA DLT/0.05d0/ ! Parameter for theta-integration (Gauss-Legendre rule). Usually 0.05
- DATA NPOI/125/ ! Total number of shapefunction points
+ DATA NPOI/555/ ! Total number of shapefunction points
DATA NRAD/10/ ! Muffintinization points
DATA NMIN/7/ ! Minimum number of points in panel
DATA NSMALL/10000/ ! A large number to start (See subr. divpanels)
@@ -319,7 +323,7 @@ c
CALL READINPUT(BRAVAIS,LCARTESIAN,RBASIS,ABASIS,BBASIS,CBASIS,
& DX,DY,DZ,
& ALATC,BLATC,CLATC,
- & IRNS,NAEZ,NEMB,KAOEZ,IRM,ZATOM,SITEAT,
+ & IRNS,NAEZ,NVAC,NEMB,KAOEZ,IRM,ZATOM,SITEAT,
& INS,KSHAPE,
& LMAX,LMMAX,LPOT,
& NATYP,NSPIN,
@@ -356,6 +360,9 @@ c Rationalise basis vectors
X TRIGHT)
ENDIF
c
+ DO NVAC_IT = 1,20
+c the number 20 is an empirical value for the number of iterations used
+c to update the empty-cell positions
CALL CLSGEN_VORONOI(NATYP,NAEZ,NEMB,RR,NR,RBASIS,
& KAOEZ,ZATOM,CLS,NCLS,
& NACLS,ATOM,EZOA,
@@ -363,6 +370,7 @@ c
& ZPERIGHT,TLEFT,TRIGHT,
& RCLS,RMTHLF,RCUTZ,RCUTXY,LINTERFACE,
& ALATC)
+ CLOSE (8)
DISTNN(1:NAEZ) = 2.D0*RMTHLF(1:NAEZ)
@@ -525,12 +533,14 @@ c Therefore, sizefac(0) = 1.0 is defined earlier.
WRITE(6,*) 'Entering VORONOI12 for atom=',IAT
CALL VORONOI12(
> NVEC,RVEC,NVERTD,NFACED,WEIGHT0,WEIGHT,TOLVDIST,TOLAREA,TOLHS,
- < RMT0,ROUT,VOLUME,NFACE,A3,B3,C3,D3,NVERT,XVERT,YVERT,ZVERT)
+ < RMT0,ROUT,VOLUME,NFACE,A3,B3,C3,D3,NVERT,XVERT,YVERT,ZVERT,
+ < VCENTER)
c Now store results in atom-dependent array.
RMT0_ALL(IAT) = RMT0
ROUT_ALL(IAT) = ROUT
VOLUME_ALL(IAT) = VOLUME
+ VCENTER_ALL(:,IAT) = VCENTER(:)
NFACE_ALL(IAT) = NFACE
A3_ALL(:,IAT) = A3(:)
B3_ALL(:,IAT) = B3(:)
@@ -543,6 +553,29 @@ c Now store results in atom-dependent array.
20 ENDDO ! DO 20 IAT = 1,NSITES
+
+ IF(NVAC.GT.0) THEN
+ OPEN(333,file='empty_cell.dat',form='formatted')
+ WRITE(6,FMT='(I6,A)') NVAC,
+ + ' empty cell positions will be updated'
+ VCENTER_SQSUM = 0.0D0
+ DO IAT = 1,NVAC
+ VCENTER_SQSUM = VCENTER_SQSUM + SQRT(VCENTER_ALL(1,IAT)**2+
+ + VCENTER_ALL(2,IAT)**2+VCENTER_ALL(3,IAT)**2)
+c an empirical factor 0.2D0 is used to avoid overshooting of the iterations
+ DO J = 1,3
+ VCENTER_ALL(J,IAT) = VCENTER_ALL(J,IAT)*0.2D0
+ END DO
+ RBASIS(:,IAT) = RBASIS(:,IAT) + VCENTER_ALL(:,IAT)
+ WRITE(6,FMT='(3F16.9)') RBASIS(:,IAT)
+ WRITE(333,FMT='(3F16.9)') RBASIS(:,IAT)
+ END DO
+ WRITE(6,FMT='(F16.9,A)') VCENTER_SQSUM,
+ + ' is a quality measure for the empty cell positions'
+ CLOSE(333)
+ END IF
+ IF(NVAC.EQ.0.OR.ABS(VCENTER_SQSUM).LT.1.D-6) EXIT
+ END DO
c-------------------------------------------------------------------------------
diff --git a/source/voronoi/readinput12.f90 b/source/voronoi/readinput12.f90
index ad79b19414150f440bec8bfcc18d418252194eb3..88f000580e51756b07d60a05c2beea211f8013db 100644
--- a/source/voronoi/readinput12.f90
+++ b/source/voronoi/readinput12.f90
@@ -1,7 +1,7 @@
SUBROUTINE READINPUT(BRAVAIS,LCARTESIAN,RBASIS,ABASIS,BBASIS,CBASIS, &
& DX,DY,DZ, &
& ALATC,BLATC,CLATC, &
- & IRNS,NAEZ,NEMB,KAOEZ,IRM,ZAT,SITEAT, &
+ & IRNS,NAEZ,NVAC,NEMB,KAOEZ,IRM,ZAT,SITEAT, &
& INS,KSHAPE, &
& LMAX,LMMAX,LPOT, &
& NATYP,NSPIN, &
@@ -55,7 +55,7 @@
& KVREL,KWS,KXC,LMAX,LMMAX,LMPOT,LPOT,MD, &
& NATYP,NPNT1,NPNT2,NPNT3,NPOL,NSPIN,INDX,IAT
INTEGER NMIN,NSMALL,NRAD,NFACELIM,NBR
- INTEGER NSTEPS,KMT,NAEZ,NEMB
+ INTEGER NSTEPS,KMT,NAEZ,NVAC,NEMB
INTEGER NINEQ,NEMBZ,NZ,CENTEROFINV(3)
REAL*8 ALATC,BLATC,CLATC
INTEGER MMIN,MMAX,SINN,SOUT,RIN,ROUT
@@ -149,6 +149,13 @@
WRITE(*,*) 'readinput: CARTESIAN=',LCARTESIAN
WRITE(111,*) 'CARTESIAN= ',LCARTESIAN
+ CALL IoInput('NVAC ',UIO,1,7,IER)
+ IF (IER.EQ.0) THEN
+ READ (UNIT=UIO,FMT=*) NVAC
+ WRITE(*,*) 'NVAC=', NVAC
+ ELSE
+ NVAC = 0
+ ENDIF
CALL IoInput('NAEZ ',UIO,1,7,IER)
IF (IER.EQ.0) THEN
diff --git a/source/voronoi/voronoi12.f b/source/voronoi/voronoi12.f
index 4c26d721479e7ae83251b451b7bff0844e9e4bec..35a6b69dd3e611eb9dd7a20d07ebb09be1f255f2 100644
--- a/source/voronoi/voronoi12.f
+++ b/source/voronoi/voronoi12.f
@@ -1,7 +1,8 @@
c***********************************************************************
SUBROUTINE VORONOI12(
> NVEC,RVEC,NVERTMAX,NFACED,WEIGHT0,WEIGHT,TOLVDIST,TOLAREA,TOLHS,
- < RMT,ROUT,VOLUME,NFACE,A3,B3,C3,D3,NVERT,XVERT,YVERT,ZVERT)
+ < RMT,ROUT,VOLUME,NFACE,A3,B3,C3,D3,NVERT,XVERT,YVERT,ZVERT,
+ < VCENTER)
c Given a cluster of atomic positions at RVEC(3,NVEC), this subroutine
c returns information about the Voronoi cell around the origin. It is
c supposed, of course, that the origin corresponds to an atomic position
@@ -73,7 +74,7 @@ c ! and of all others (dimensioned as RVEC).
c Output:
INTEGER NFACE
INTEGER NVERT(*)
- REAL*8 VOLUME
+ REAL*8 VOLUME,VCENTER(3)
REAL*8 A3(NFACED),B3(NFACED),C3(NFACED),D3(NFACED)
REAL*8 XVERT(NVERTMAX,NFACED),YVERT(NVERTMAX,NFACED),
& ZVERT(NVERTMAX,NFACED)
@@ -133,6 +134,9 @@ c origin and r1,r2,r3 the vectors of the 3 other vertices.
c Algorithm requires that the face is a convex polygon with the
c vertices ordered (doesn't matter if they are clock- or anticlockwise).
VOLUME = 0.d0
+ VCENTER(1) = 0.d0
+ VCENTER(2) = 0.d0
+ VCENTER(3) = 0.d0
DO IFACE = 1,NFACE
X1 = XVERT(1,IFACE)
Y1 = YVERT(1,IFACE)
@@ -146,6 +150,9 @@ c vertices ordered (doesn't matter if they are clock- or anticlockwise).
Y3 = YVERT(IVERT+1,IFACE)
Z3 = ZVERT(IVERT+1,IFACE)
TETRVOL = X1*(Y2*Z3-Y3*Z2)+X2*(Y3*Z1-Y1*Z3)+X3*(Y1*Z2-Y2*Z1)
+ VCENTER(1) = VCENTER(1) + 0.25d0*DABS(TETRVOL)*(X1+X2+X3)
+ VCENTER(2) = VCENTER(2) + 0.25d0*DABS(TETRVOL)*(Y1+Y2+Y3)
+ VCENTER(3) = VCENTER(3) + 0.25d0*DABS(TETRVOL)*(Z1+Z2+Z3)
VOLUME = VOLUME + DABS(TETRVOL)
TRIANGLEAREA = 0.5d0 * DSQRT(
& ( X1*Y2 + X2*Y3 + X3*Y1 - Y2*X3 - Y3*X1 - Y1*X2)**2
@@ -157,8 +164,12 @@ c vertices ordered (doesn't matter if they are clock- or anticlockwise).
ENDDO
ENDDO
VOLUME = VOLUME/6.D0
+ VCENTER(1) = VCENTER(1)/VOLUME
+ VCENTER(2) = VCENTER(2)/VOLUME
+ VCENTER(3) = VCENTER(3)/VOLUME
WRITE(6,*) ' Polyhedron properties '
+
WRITE(6,*) ' Number of faces : ',nface
IF (TEST('verb0 ')) THEN