[Wien] Bug in latest version of Supercell

Peter Blaha pblaha at zeus.theochem.tuwien.ac.at
Tue Dec 7 09:54:17 CET 2004


> I am using the latest version (04.11) of wien2k compiled on Fedora Core 1
> system using the Lahey compiler. I wanted to generate supercell of the GaN
> structure which is
> 
> Spgr 186, a=b=3.189 angs, c=5.185 angs, 
> Ga at 1/3,2/3,0.3125, 
> N  at 1/3,2/3,0.
> 
> When I generate the (2x2x2 or 3x3x3) supercell by issuing the command 'x
> supercell', the program makes the new struct file by using the H symmetry,
> correct lattice parameters (2a,2c or 3a,3c), and generates the number of
> atoms correctly (16: 8Ga & 8N, or 54: 27Ga & 27N), but puts the fractional
> coordinates of all the 16 or 54 atoms as 0,0,0. It also asks for the
> length of vacuum region in a.u. which I give as 0.

I could verify your problem. 

The H-lattice support worked only with vacuum, but not without additional
vacuum. I've fixed it and supercell.f is included.

It will also be on the web with the next update.



                                      P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-15671             FAX: +43-1-58801-15698
Email: blaha at theochem.tuwien.ac.at    WWW: http://info.tuwien.ac.at/theochem/
--------------------------------------------------------------------------
-------------- next part --------------
      program supercell

!     ########################################      
!
!     Program generates supercell from a WIEN struct file.
!     (c) 2002 Alexander Riss (riss at m3plus.com)
!
!     ########################################


     
      IMPLICIT NONE

      character*60    fname
      character*70    fname2
      character*100   prob
      integer         numx,numy,numz,iatom,i,j,k,NATtar,nattar0
      real*8          vac(3),vacfac(3)

      integer         maxnum
      parameter       (maxnum=500)
      
      character*80    TITLE
      character*4     LATTIC, lattar, IREL,CFORM, UNIT
      character*1     toplayer(3)
      real*8          AA,BB,CC,ALPHA(3),ROTLOC(3,3)
      real*8          POS(3,maxnum),oldpos(3,maxnum),bakpos(3,maxnum)
      real*8          RMT,R0,ZZ
      integer         NAT,IATNR(maxnum),MULT,ISPLIT,JRI
      character*10    ANAME
      integer         ATNR

      
      toplayer='N'      
!     ###
!     ### program input

      write(*,*) 'Program generates supercell from a WIEN struct file.'
      write(*,*)
      write (*,*) 'Filename of struct file: '
      read (*,17) fname
   17 format(A60)

      open (20,FILE=fname,STATUS='old',FORM='formatted',ERR=500)
      read(20,1000,ERR=505) TITLE
      read(20,1010,ERR=505) LATTIC,NAT,CFORM,IREL,UNIT
      cform='    '
      read(20,1020,ERR=505) AA,BB,CC,ALPHA(1),ALPHA(2),ALPHA(3)

      if (LATTIC.ne.'P' .and. LATTIC.ne.'B' .and.LATTIC.ne.'F' .and. & 
          LATTIC.ne.'H'.and. LATTIC.ne.'CXY'.and. LATTIC.ne.'CXZ') goto 510

      write (*,*)
      write (*,*) 'Number of cells in x direction: '
      read (*,*) numx
      write (*,*) 'Number of cells in y direction: '
      read (*,*) numy
      write (*,*) 'Number of cells in z direction: '
      read (*,*) numz

!     ### read in target lattice type
      write (*,*)
      write (*,*) 'Current structure has lattice type ' // LATTIC
  101 if (LATTIC.eq.'H') then
        lattar=LATTIC
        if(numx.ne.numy) lattar='P'
      else if (LATTIC.eq.'CXY') then
        lattar='P'
      else if (LATTIC.eq.'CXZ') then
        lattar='P'
      else if((numx.eq.numy).and.(numx.eq.numz).and.(numx/2e0.eq.numx/2)) then
        write (*,*) 'Enter your target lattice type: (P,B,F)'
	      read (*,'(A4)') lattar
        if (lattar.eq.'') lattar=LATTIC
        if (lattar.eq.'p') lattar='P'
        if (lattar.eq.'f') lattar='F'
        if (lattar.eq.'b') lattar='B'
        if ((lattar.ne.'P').and.(lattar.ne.'B').and. (lattar.ne.'F')) goto 101
      else if ((numx.eq.numy).and.(numx.eq.numz)) then
        write (*,*) 'Enter your target lattice type: (P,B)'
        read (*,'(A4)') lattar
        if (lattar.eq.'') lattar=LATTIC
        if (lattar.eq.'p') lattar='P'
        if (lattar.eq.'b') lattar='B'
        if ((lattar.ne.'P').and.(lattar.ne.'B')) lattar='P'
!        if ((lattar.ne.'P').and.(lattar.ne.'B')) goto 101
      else
        lattar='P'
      end if
      write (*,*) 'Target lattice type will be ' // lattar


!     ### user may add vacuum  in x-, y- or z-direction
      if (lattar.eq.'P') then
        write(*,*)
 102    write (*,*) 'Add vacuum in x-direction for surface-slab [bohr]:'
        read (*,*,err=102) vac(1)
        if(vac(1).gt.0.00001d0) then
          write(*,*) 'Repeat atoms at x=0 at the top (y/n)'
          read(*,'(a)') toplayer(1)
          if(toplayer(1).eq.'y') toplayer(1)='Y'
        endif
 103    write (*,*) 'Add vacuum in y-direction for surface-slab [bohr]:'
        read (*,*,err=103) vac(2)
        if(vac(2).gt.0.00001d0) then
          write(*,*) 'Repeat atoms at y=0 at the top (y/n)'
          read(*,'(a)') toplayer(2)
          if(toplayer(2).eq.'y') toplayer(2)='Y'
        endif
 104    write (*,*) 'Add vacuum in z-direction for surface slab [bohr]:'
        read (*,*,err=104) vac(3)
        if(vac(3).gt.0.00001d0) then
          write(*,*) 'Repeat atoms at z=0 at the top (y/n)'
          read(*,'(a)') toplayer(3)
          if(toplayer(3).eq.'y') toplayer(3)='Y'
        endif
        vacfac(1)=AA*numx/(vac(1)+AA*numx)
        vacfac(2)=BB*numy/(vac(2)+BB*numy)
        vacfac(3)=CC*numz/(vac(3)+CC*numz)
      else if (lattar.eq.'H') then
        write(*,*)
 204    write (*,*) 'Add vacuum in z-direction for surface slab [bohr]:'
        read (*,*,err=204) vac(3)
        if(vac(3).gt.0.00001d0) then
          write(*,*) 'Repeat atoms at z=0 at the top (y/n)'
          read(*,'(a)') toplayer(3)
          if(toplayer(3).eq.'y') toplayer(3)='Y'
        endif
        vacfac(1)=1.d0
        vacfac(2)=1.d0
        vacfac(3)=CC*numz/(vac(3)+CC*numz)
      else
        vacfac(1)=1.d0
        vacfac(2)=1.d0
        vacfac(3)=1.d0
      end if
      
!     ### calculation of number of atoms in target file
      NATtar=NAT*numx*numy*numz
      if (LATTIC.eq.'P' .and. lattar.eq.'B') NATtar=NATtar/2     
      if (LATTIC.eq.'P' .and. lattar.eq.'F') NATtar=NATtar/4     
      if (LATTIC.eq.'B' .and. lattar.eq.'P') NATtar=NATtar*2     
      if (LATTIC.eq.'B' .and. lattar.eq.'F') NATtar=NATtar/2     
      if (LATTIC.eq.'F' .and. lattar.eq.'P') NATtar=NATtar*4     
      if (LATTIC.eq.'F' .and. lattar.eq.'B') NATtar=NATtar*2     
      if (LATTIC.eq.'CXY' .and. lattar.eq.'P') NATtar=NATtar*2     
      if (LATTIC.eq.'CXZ' .and. lattar.eq.'P') NATtar=NATtar*2     
      
      i=index(fname,'.')
      if (i.eq.0) i=1
      j=len(fname)
      fname2=fname(1:i-1) // '_super' // fname(i:j)
      
      AA=(AA*numx)+vac(1)
      BB=(BB*numy)+vac(2)
      CC=(CC*numz)+vac(3)
      write(21,1000,ERR=525) TITLE
      write(21,2010,ERR=525) lattar,NATtar,CFORM,IREL,UNIT
      write(21,2020,ERR=525) AA,BB,CC,ALPHA(1),ALPHA(2),ALPHA(3)

!    ### read atom information, call subroutine to write information every atom
      ATNR=0
      do i=1,NAT
        iatom=1
        read(20,1030,ERR=505) IATNR(iatom),(POS(k,iatom),k=1,3),MULT,ISPLIT
        do j=1,MULT-1
          iatom=iatom+1
          if (iatom.gt.maxnum) goto 530
          read(20,1031,ERR=505) IATNR(iatom),(POS(k,iatom),k=1,3)
        end do
        read(20,1050,ERR=505) ANAME,JRI,R0,RMT,ZZ 	
        aname(3:10)='        '
        read(20,1051,ERR=505) ((ROTLOC(k,j),k=1,3),j=1,3)
        call startcell(fname2,maxnum,numx,numy,numz,iatom,NAT, &
          ROTLOC,ISPLIT,MULT,IATNR,POS,oldpos,bakpos,vacfac, &
          ANAME,JRI,R0,RMT,ZZ,ATNR,LATTIC,lattar,toplayer)
      end do

!     correct for actual number of atoms, reread new struct file
      rewind 21
      read(21,1000,ERR=525) TITLE
      read(21,1010,ERR=525) lattar,NATtar,CFORM,IREL,UNIT
      read(21,2020,ERR=525) AA,BB,CC,ALPHA(1),ALPHA(2),ALPHA(3)
      do i=1,1000000
         iatom=1
         read(21,1030,end=600) IATNR(iatom),(POS(k,iatom),k=1,3),MULT,ISPLIT
        do j=1,MULT-1
          iatom=iatom+1
          if (iatom.gt.maxnum) goto 530
          read(21,1031,ERR=505) IATNR(iatom),(POS(k,iatom),k=1,3)
        end do
        read(21,1050,ERR=505) ANAME,JRI,R0,RMT,ZZ 	
        read(21,1051,ERR=505) ((ROTLOC(k,j),k=1,3),j=1,3)
      enddo
 600  nattar=i-1
      rewind 21
      open (22,FILE=fname2,STATUS='unknown',FORM='formatted',ERR=520)
      read(21,1000,ERR=525) TITLE
      read(21,1010,ERR=525) lattar,NATtar0,CFORM,IREL,UNIT
      read(21,2020,ERR=525) AA,BB,CC,ALPHA(1),ALPHA(2),ALPHA(3)
      write(22,1000) TITLE
      write(22,2010) lattar,NATtar,CFORM,IREL,UNIT
      write(22,2020) AA,BB,CC,ALPHA(1),ALPHA(2),ALPHA(3)
      do i=1,nattar
         iatom=1
         read(21,1030,end=600) IATNR(iatom),(POS(k,iatom),k=1,3),MULT,ISPLIT
         write(22,2030) IATNR(iatom),(POS(k,iatom),k=1,3),MULT,ISPLIT
        do j=1,MULT-1
          iatom=iatom+1
          if (iatom.gt.maxnum) goto 530
          read(21,1031,ERR=505) IATNR(iatom),(POS(k,iatom),k=1,3)
          write(22,2031) IATNR(iatom),(POS(k,iatom),k=1,3)
        end do
        read(21,1050,ERR=505) ANAME,JRI,R0,RMT,ZZ 	
        read(21,1051,ERR=505) ((ROTLOC(k,j),k=1,3),j=1,3)
        write(22,2050) ANAME,JRI,R0,RMT,ZZ 	
        write(22,2051) ((ROTLOC(k,j),k=1,3),j=1,3)
      enddo
      
!     ### set symmetry operations to zero
      write(22,2060,ERR=525) 0

      write(*,*)
      write(*,*) 'Supercell generated sucessfully.'
      write(*,*) 'Stored in struct file: ' // fname2
      write(*,*) 'You may need to replace an atom by an impurity or', &
                 ' distort the positions, ....'
      
      stop


!     ### errors leading to termination of program
  500 prob='Error accessing file: ' // fname
      call faterr(prob)
      stop
  505 prob='Error reading file: ' // fname
      call faterr(prob)
      stop
  510 prob='Unknown lattice type: ' // LATTIC
      call faterr(prob)
      stop
  520 prob='Error creating file: ' // fname2
      call faterr(prob)
      stop
  525 prob='Error writing to file: ' // fname2
      call faterr(prob)
      stop
  530 prob='Number of Atoms exceeds array memory.'
      call faterr(prob)
      stop


 1000 FORMAT(A80)
 1010 FORMAT(A4,23X,I3,1X,A4,/,13X,A4,6X,A4)
 2010 FORMAT(A4,'LATTICE,NONEQUIV. ATOMS',I3,1X,A4,/,'MODE OF CALC=',A4,' unit=',A4)
 1020 FORMAT(6F10.6)
 2020 FORMAT(6F10.6)
 1030 FORMAT(4X,I4,4X,F10.8,3X,F10.8,3X,F10.8,/,15X,I2,17X,I2)
 1031 FORMAT(4X,I4,4X,F10.8,3X,F10.8,3X,F10.8)
 1050 FORMAT(A10,5X,I5,5X,F10.9,5X,F10.5,5X,F10.5)
 1051 FORMAT(20X,3F10.8)
 2060 FORMAT(I4,6X,'NUMBER OF SYMMETRY OPERATIONS')
 2030 FORMAT('ATOM',I4,': X=',F10.8,' Y=',F10.8,' Z=',F10.8,/,10X,'MULT=',I2,10X,'ISPLIT=',I2)
 2031 FORMAT('ATOM',I4,': X=',F10.8,' Y=',F10.8,' Z=',F10.8)
 2050 FORMAT(A10,1X,'NPT=',I5,2X,'R0=',F10.8,1X,'RMT=',F10.4,3X,'Z:',F5.1)
 2051 FORMAT('LOCAL ROT MATRIX:',3X,3F10.7,/,20X,3F10.7,/,20X,3F10.7)
 
      end

      subroutine faterr(prob)
        character*80 prob
        write(*,*)
        write(*,*) 'Fatal Error occured:'
        write(*,*) prob
        write(*,*) 'Program terminated.'
        write(*,*)
        return
      end

      
      
!     ###      
!     ### converts integer to a string
      character*4 function intstr(i)
        integer i
        if (i.lt.10) then
          write(intstr,'(I1)') i
        else if (i.lt.100) then
          write(intstr,'(I2)') i
        else if (i.lt.1000) then
          write(intstr,'(i3)') i
        else
          write(intstr,'(i4)') i
        end if
      end
      

      
!     ###
!     ### transform body-centered or face-centered cells to primitive cells,
!     ### then call supercell for multiplication and reducing to target lattice type as necessary
      subroutine startcell(fname2,maxnum,numx,numy,numz,iatom,  &
         NAT,ROTLOC,ISPLIT,MULT,IATNR,POS,oldpos,bakpos,vacfac, &
         ANAME,JRI,R0,RMT,ZZ,ATNR,LATTIC,lattar,toplayer)
        integer       maxnum,IATNR(maxnum),MULT,ISPLIT,JRI,NAT,ATNR
        integer       i,j,iatom,nlatA,nlatB
        integer       numx,numy,numz,iname
        real*8        POS(3,maxnum),oldpos(3,maxnum),bakpos(3,maxnum)
        real*8        RMT,R0,ZZ,ROTLOC(3,3)
        real*8        vacfac(3),latAdd(3,11)
        character*10  ANAME
        character*70  fname2
        character*4   LATTIC,lattar
        character*1   toplayer(3)

!       ### get from body-cenetred or face-centered lattice type to primitive lattice type by generating appropriate atoms 
        
!       ## primitive: atom1=atom 
        latAdd(1,1)=0.d0
        latAdd(2,1)=0.d0
        latAdd(3,1)=0.d0

!       ## body-centered: atom1=atom, atom2=atom+(1/2,1/2,1/2)        
        latAdd(1,2)=0.d0
        latAdd(2,2)=0.d0
        latAdd(3,2)=0.d0
        latAdd(1,3)=0.5d0
        latAdd(2,3)=0.5d0
        latAdd(3,3)=0.5d0
        
!       ## face-centered: atom1=atom, atom2=atom+(1/2,1/2,0), atom3=atom+(1/2,0,1/2), atom4=atom+(0,1/2,1/2)
        latAdd(1,4)=0.d0
        latAdd(2,4)=0.d0
        latAdd(3,4)=0.d0
        latAdd(1,5)=0.5d0
        latAdd(2,5)=0.5d0
        latAdd(3,5)=0.d0
        latAdd(1,6)=0.5d0
        latAdd(2,6)=0.d0
        latAdd(3,6)=0.5d0
        latAdd(1,7)=0.d0
        latAdd(2,7)=0.5d0
        latAdd(3,7)=0.5d0
!       ## cxy centered: atom1=atom, atom2=atom+(1/2,1/2,0)
        latAdd(1,8)=0.d0
        latAdd(2,8)=0.d0
        latAdd(3,8)=0.d0
        latAdd(1,9)=0.5d0
        latAdd(2,9)=0.5d0
        latAdd(3,9)=0.d0
!       ## cxz centered: atom1=atom, atom2=atom+(1/2,0,1/2)
        latAdd(1,10)=0.d0
        latAdd(2,10)=0.d0
        latAdd(3,10)=0.d0
        latAdd(1,11)=0.5d0
        latAdd(2,11)=0.d0
        latAdd(3,11)=0.5d0

        if (LATTIC.eq.'B') then
          nlatA=2
          nlatB=3
        else if (LATTIC.eq.'F') then
          nlatA=4
          nlatB=7
        else if (LATTIC.eq.'CXY') then
          nlatA=8
          nlatB=9
        else if (LATTIC.eq.'CXZ') then
          nlatA=10
          nlatB=11
        else
          nlatA=1
          nlatB=1
        end if

        
!       ### save old coordinates
        do i=1,MULT
          do j=1,3
            bakpos(j,i)=POS(j,i)
          end do
        end do
        
        iname=0 

!       ### add positions and call makecell
        do i=nlatA,nlatB
          do iatom1=1,MULT
            do j=1,3
              POS(j,iatom1)=bakpos(j,iatom1)+latAdd(j,i)
              if (POS(j,iatom1).ge.1) POS(j,iatom1)=POS(j,iatom1)-1 
            end do
          end do
          call makecell(fname2,maxnum,numx,numy,numz,iatom,NAT, &
            ROTLOC,ISPLIT,MULT,IATNR,POS,oldpos,vacfac,         &
            ANAME,JRI,R0,RMT,ZZ,ATNR,LATTIC,lattar,iname,toplayer)
        end do
        
        return
      end
      

      
!     ###
!     ### make supercell
      subroutine makecell(fname2,maxnum,numx,numy,numz,iatom, &
         NAT,ROTLOC,ISPLIT,MULT,IATNR,POS,oldpos,vacfac,      &
         ANAME,JRI,R0,RMT,ZZ,ATNR,LATTIC,lattar,iname,toplayer)
        integer       maxnum,IATNR(maxnum),MULT,ISPLIT,JRI,NAT,ATNR
        integer       i,j,ix,iy,iz,iatom
        integer       laname,iname,curmul
        integer       numx,numy,numz
        real*8        POS(3,maxnum),oldpos(3,maxnum)
        real*8        posx,posy,posz,vacfac(3)
        real*8        RMT,R0,ZZ,ROTLOC(3,3)
        character*10  ANAME,aname2
        character*100 prob
        character*70  fname2
        character*4   intstr
        character*4   LATTIC,lattar
        character*1   toplayer(3)
        logical       cut           

!       ### strip whitspaces from ANAME
        laname=1
        do i=1,len(ANAME)
          if (ANAME(i:i).ne.' ') laname=i
        end do
        if (laname.gt.8) laname=8
!       ### names should not be changed to third position (otherwise they would not be interpeted correctly)	
	      if (laname.lt.3) laname=3

!       ### save old coordinates
        do i=1,iatom
          do j=1,3
            oldpos(j,i)=POS(j,i)
          end do
        end do
        
        iz1=numz-1
        iy1=numy-1
        ix1=numx-1
        if(toplayer(1).eq.'Y') ix1=numx
        if(toplayer(2).eq.'Y') iy1=numy
        if(toplayer(3).eq.'Y') iz1=numz
!       ### calculate new coordinates now
        do iz=0,numz
          do iy=0,numy
            do ix=0,numx
            
!             ### check if some atoms should be "reduced", when transforming to a face-centered or body-centered type lattice
              iatom=1
              iatom1=1
              curmul=MULT
              do i=1,MULT
                cut=.false.
        if((iz.eq.numz).and.(toplayer(3).ne.'Y'.or.oldpos(3,iatom).ne.0.d0)) cut=.true.
        if((iy.eq.numy).and.(toplayer(2).ne.'Y'.or.oldpos(2,iatom).ne.0.d0)) cut=.true.
        if((ix.eq.numx).and.(toplayer(1).ne.'Y'.or.oldpos(1,iatom).ne.0.d0)) cut=.true.
                posx=0.d0+oldpos(1,iatom)/numx+1.d0*ix/numx
                posy=0.d0+oldpos(2,iatom)/numy+1.d0*iy/numy
                posz=0.d0+oldpos(3,iatom)/numz+1.d0*iz/numz
                iatom=iatom+1
                if (lattar.eq.'B') then
                  if (posz.ge.0.5d0) cut=.true.
                else if (lattar.eq.'F') then
                  if (posz.ge.0.5d0 .or. posy.ge.0.5d0) cut=.true.
                end if
                if (cut) then
                  curmul=curmul-1
                  goto 201
                end if
                POS(1,iatom1)=posx*vacfac(1)
                POS(2,iatom1)=posy*vacfac(2)
                POS(3,iatom1)=posz*vacfac(3)
               iatom1=iatom1+1
  201         end do
              
              if (curmul.lt.1) goto 202
              
!             ### get name of atom              
              ATNR=ATNR+1
              iname=iname+1
              if (iname.gt.100) then
                if (laname.gt.6) laname=6
                    end if
              if (iname.gt.10) then
                if (laname.gt.7) laname=7
              end if
              aname2= ANAME(1:laname) // '_' // intstr(iname)

!             ## write to file
              do i=1,curmul
                if (i.eq.1) then
                  write(21,2030,ERR=625) ATNR,(POS(k,i),k=1,3),curmul,ISPLIT
                else
                  write(21,2031,ERR=625)ATNR,(POS(k,i),k=1,3)
                end if
              end do
              write(21,2050,ERR=625) aname,JRI,R0,RMT,ZZ 	
              write(21,2051,ERR=625) ((ROTLOC(k,j),k=1,3),j=1,3)

  202       end do
          end do
        end do

        return
        
!       ### has to be defined here again

 2030   FORMAT('ATOM',I4,': X=',F10.8,' Y=',F10.8,' Z=',F10.8,/,10X, &
               'MULT=',I2,10X,'ISPLIT=',I2)
 2031   FORMAT('ATOM',I4,': X=',F10.8,' Y=',F10.8,' Z=',F10.8)
 2050   FORMAT(A10,1X,'NPT=',I5,2X,'R0=',F10.8,1X,'RMT=',F10.4,3X,'Z:',F5.1)
 2051   FORMAT('LOCAL ROT MATRIX:',3X,3F10.7,/,20X,3F10.7,/,20X,3F10.7)
 
  625   prob='Error writing to file: ' // fname2
        call faterr(prob)
        stop

      end




More information about the Wien mailing list