C C xv2scl converts a given XV into an expanded XV (of the supercell), C according to the translation vectors C defined by the (integer) supercell matrix. C (Note: the input/output format of the structure data C is assumed to be that of the Siesta XV file. C The conversion of the result into the .fdf convention C can be achieved by a subsequent use of my xv2fdf script C on the resulting XV file). C The program strictly checks that all the atoms, C whatever their original placement with respect C to the prototype unit cell, do fall within the supercell C spanned by the resulting translation vectors. C A certain order of the atoms' appearance in the supercell C is enforced, based on their X,Y,Z coordinates; C this heuristics can be easily changed in the code. C C Usage: start executable without passing any parameters, C and type in the necessary when prompted C (here: name of the XV file; 3*3 integer supercell matrix). C C Written by Andrei Postnikov, Nov 2009 C postnikov@univ-metz.de C program xv2scl implicit none integer ii1,io1,ialloc,mode parameter (ii1=11,io1=14) integer ii,jj,kk,ll,iat,nat,nmult,icount,icount1,ncount integer mvec(3,3),mpt(3,8),min_m(3),max_m(3) character inpfil*60,outfil*60,logfil*60,syslab*30 double precision, allocatable :: coord(:,:), coords(:,:), rang(:) integer, allocatable :: nz(:), ityp(:), isort(:) double precision cc_ini(3,3),cc_scl(3,3),cc_inv(3,3),cmult(3), . cvol,base(3),coor_try(3),proj_try(3),rang0 external opnout C C string manipulation functions in Fortran used below: C len_trim(string): returns the length of string C without trailing blank characters, C trim(string) : returns the string with railing blanks removed write (6,701,advance="no") read (5,*) syslab C inpfil = syslab(1:len_trim(syslab))//'.XV' inpfil = trim(syslab)//'.XV' open (ii1,file=inpfil,form='formatted',status='old',err=801) write (6,*) 'Found and opened: ',inpfil C --- open output files: outfil = trim(syslab)//'_expand.XV' call opnout(io1,outfil) C -- read in translation vectors: do ii=1,3 read (ii1,702,end=802,err=802) (cc_ini(jj,ii),jj=1,3) enddo C -- calculate cell volume: cvol = cc_ini(1,1)*cc_ini(2,2)*cc_ini(3,3) + + cc_ini(1,2)*cc_ini(2,3)*cc_ini(3,1) + + cc_ini(1,3)*cc_ini(2,1)*cc_ini(3,2) - - cc_ini(1,1)*cc_ini(2,3)*cc_ini(3,2) - - cc_ini(1,2)*cc_ini(2,1)*cc_ini(3,3) - - cc_ini(1,3)*cc_ini(2,2)*cc_ini(3,1) read (ii1,703,end=803,err=803) nat C--- allocate arrays for atoms in the original cell: allocate (coord(3,nat),STAT=ialloc) allocate (nz(nat),STAT=ialloc) allocate (ityp(nat),STAT=ialloc) if (ialloc.ne.0) then write (6,*) ' Fails to allocate space for ',nat,' atoms.' stop endif do iat=1,nat read (ii1,704,end=804,err=804) ityp(iat), nz(iat), . (coord(jj,iat),jj=1,3) enddo write (6,*) 'Now specify an (integer) 3*3 supercell matrix.' write (6,*) 'Its each line defines one vector of the supercell' write (6,*) 'in terms of three vectors ', . 'of the original simple cell.' 101 write (6,*) 'Type in the coordinates for the first SC vector:' read (5,*,err=101,end=101) mvec(:,1) 102 write (6,*) 'Type in the coordinates for the second SC vector:' read (5,*,err=102,end=102) mvec(:,2) 103 write (6,*) 'Type in the coordinates for the third SC vector:' read (5,*,err=103,end=103) mvec(:,3) nmult = mvec(1,1)*mvec(2,2)*mvec(3,3) + + mvec(1,2)*mvec(2,3)*mvec(3,1) + + mvec(1,3)*mvec(2,1)*mvec(3,2) - - mvec(1,1)*mvec(2,3)*mvec(3,2) - - mvec(1,2)*mvec(2,1)*mvec(3,3) - - mvec(1,3)*mvec(2,2)*mvec(3,1) if (nmult.eq.0) then write (6,*) ' These vectors are linear dependent; try better!' goto 101 elseif (nmult.lt.0) then write (6,*) ' Left-hand vector triplet; try another sequence:' goto 101 endif write (6,*) ' nmult =',nmult C -- construct supercell vectors do ii=1,3 cc_scl(:,ii) = cc_ini(:,1)*mvec(1,ii) + + cc_ini(:,2)*mvec(2,ii) + + cc_ini(:,3)*mvec(3,ii) write (io1,301) cc_scl(:,ii) enddo 301 format(3x,3f18.9) write (io1,'(i12)') nat*nmult C--- allocate arrays for atoms in the supercell: allocate (isort(nmult*nat),STAT=ialloc) allocate (rang(nmult*nat),STAT=ialloc) allocate (coords(3,nmult*nat),STAT=ialloc) if (ialloc.ne.0) then write (6,*) ' Fails to allocate space for ',nmult*nat, . ' atoms in the supercell.' stop endif C -- construct the inverse matrix to cc_scl. Its elements are C minors of cc_scl, divided by nmult*cvol. cc_inv(1,1) = cc_scl(2,2)*cc_scl(3,3) - cc_scl(3,2)*cc_scl(2,3) cc_inv(1,2) = cc_scl(3,2)*cc_scl(1,3) - cc_scl(1,2)*cc_scl(3,3) cc_inv(1,3) = cc_scl(1,2)*cc_scl(2,3) - cc_scl(2,2)*cc_scl(1,3) cc_inv(2,1) = cc_scl(2,3)*cc_scl(3,1) - cc_scl(3,3)*cc_scl(2,1) cc_inv(2,2) = cc_scl(3,3)*cc_scl(1,1) - cc_scl(1,3)*cc_scl(3,1) cc_inv(2,3) = cc_scl(1,3)*cc_scl(2,1) - cc_scl(2,3)*cc_scl(1,1) cc_inv(3,1) = cc_scl(2,1)*cc_scl(3,2) - cc_scl(3,1)*cc_scl(2,2) cc_inv(3,2) = cc_scl(3,1)*cc_scl(1,2) - cc_scl(1,1)*cc_scl(3,2) cc_inv(3,3) = cc_scl(1,1)*cc_scl(2,2) - cc_scl(2,1)*cc_scl(1,2) cc_inv(:,:) = cc_inv(:,:)/(nmult*cvol) C -- debugging C write (6,*) ' inverted matrix:' C do ii = 1,3 C write (6,'(3f13.6)') (cc_inv(ii,jj),jj=1,3) C enddo C -- check multiplication cc_inv * cc_scl C write (6,*) ' Product:' C do ii = 1,3 C cmult(:) = cc_inv(ii,1)*cc_scl(1,:) + C + cc_inv(ii,2)*cc_scl(2,:) + C + cc_inv(ii,3)*cc_scl(3,:) C write (6,'(3f13.6)') cmult C enddo C C -- Find limits of variation of components of supercell vectors C For this, we have to check the coordinates of 8 corner points C of the parallelepiped built by vectors mvec(:,:). min_m = 0 max_m = 0 mpt(:,1) = 0 mpt(:,2) = mvec(:,1) mpt(:,3) = mvec(:,2) mpt(:,4) = mvec(:,3) mpt(:,5) = mvec(:,2) + mvec(:,3) mpt(:,6) = mvec(:,3) + mvec(:,1) mpt(:,7) = mvec(:,1) + mvec(:,2) mpt(:,8) = mvec(:,1) + mvec(:,2) + mvec(:,3) C C write (6,*) ' Corner points:' C do jj = 1,8 C write (6,'(3i6)') (mpt(ii,jj),ii=1,3) C enddo do ii = 1,3 do jj = 1,8 if (mpt(ii,jj).gt.max_m(ii)) max_m(ii) = mpt(ii,jj) if (mpt(ii,jj).lt.min_m(ii)) min_m(ii) = mpt(ii,jj) enddo enddo C write (6,*) ' Min and Max values:' C do ii = 1,3 C write (6,'(2i8)') min_m(ii), max_m(ii) C enddo C write (6,*) ' ----------------------' C -- Run from (min_m - 1) to (max_m + 1) along all three dimensions; C allow +/-1 for the case if atomic coordinates are out of the C primitive cell. Check which atoms which fall within the supercell ncount = 0 do ii = min_m(1)-1, max_m(1)+1 do jj = min_m(2)-1, max_m(2)+1 do kk = min_m(3)-1, max_m(3)+1 base(:) = cc_ini(:,1)*ii + cc_ini(:,2)*jj + cc_ini(:,3)*kk C -- origin of the primitive cell translated by (ii,jj,kk) do iat = 1,nat coor_try(:) = base(:) + coord(:,iat) C -- coordinates of the running atom. Project them onto C the inverse matrix of the supercell. Check if the result C will be within [0..1] proj_try(:) = cc_inv(:,1)*coor_try(1) + + cc_inv(:,2)*coor_try(2) + + cc_inv(:,3)*coor_try(3) if ( (proj_try(1).ge.0.and.proj_try(1).lt.1).and. . (proj_try(2).ge.0.and.proj_try(2).lt.1).and. . (proj_try(3).ge.0.and.proj_try(3).lt.1) ) then C -- ... count this atom rang0 = coor_try(1) + coor_try(2)*500 + coor_try(3)*250000 C -- this construction helps to sort out the atoms, C according to increasing order of rang0. C Can be arbitrarily changed in this place! C -- skip those previously counted atoms whose rang C does not succeed rang0 : if ( (ncount.eq.0) C -- this is the very first atom counted . .or. (rang0.ge.rang(ncount)) C -- the newly counted atom adds to the end of list . ) then ncount = ncount + 1 isort(ncount) = iat coords(:,ncount) = coor_try(:) rang(ncount) = rang0 else C -- insert the newly counted atom some place within the list do icount = 1,ncount if (rang(icount).gt.rang0) then icount1 = icount C -- mark the first member in the list exceeding rang0 exit endif enddo C -- shift the rest upwards: do icount = ncount,icount1,-1 isort(icount+1) = isort(icount) coords(:,icount+1) = coords(:,icount) rang(icount+1) = rang(icount) enddo C -- insert the newly found atom into the place of icount1: isort(icount1) = iat coords(:,icount1) = coor_try(:) rang(icount1) = rang0 ncount = ncount + 1 endif endif enddo enddo enddo enddo if (ncount.ne.nmult*nat) then write (6,*) ' Something wrong: ', . 'Number of replicated atoms =',ncount, . ' .ne. determinant of supercell matrix =',nmult, . ' times nat=',nat stop endif C -- final printout: do icount = 1,ncount write (io1,704) ityp(isort(icount)),nz(isort(icount)), . (coords(ll,icount),ll=1,3), icount, isort(icount) enddo close (ii1) close (io1) stop 701 format(" Specify SystemLabel (or 'siesta' if none): ") 702 format (3x,3f18.9) 703 format (i12) 704 format (i3,i6,3f18.9,i5,i7) 801 write (6,*) ' Error opening file ', . trim(inpfil),' as old formatted' stop 802 write (6,*) ' End/Error reading XV for cell vector ',ii stop 803 write (6,*) ' End/Error reading XV for number of atoms line' stop 804 write (6,*) ' End/Error reading XV for atom number ',iat stop end C C............................................................................... C subroutine opnout(io1,outfil) C C Opens output file, preventing its accidental overrid C implicit none integer io1 character outfil*60,owrite*1 logical filexist inquire (file=outfil, exist=filexist) if (filexist) then write (6,*) ' File ',trim(outfil),' exists.', . ' Overwrite? (Y/N)' read (5,*) owrite if (owrite.eq.'Y'.or.owrite.eq.'y') then open (io1,file=outfil,form='formatted', . status='REPLACE',err=801) else write (6,*) ' Then rename is first. Bye...' stop endif else open (io1,file=outfil,form='formatted',status='NEW',err=802) endif return 801 write (6,*) ' Error replacing formatted file ',trim(outfil) stop 802 write (6,*) ' Error opening file ',trim(outfil), . ' as new formatted' stop end