Skip to content
Snippets Groups Projects
Commit dfa7b15c authored by Rudolf Zeller's avatar Rudolf Zeller
Browse files

New parameter NVAC determining number of empty positions being optimized by centroidal tesselation

parent 0831c8c6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
source/voronoi/maindriver12.f
+ 37
− 4
View file @ dfa7b15c
......@@ -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-------------------------------------------------------------------------------
......
source/voronoi/readinput12.f90
+ 9
− 2
View file @ dfa7b15c
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
......
source/voronoi/voronoi12.f
+ 13
− 2
View file @ dfa7b15c
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
......
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