SUBROUTINE WAVOFR( NA, NO, NUO, MAXNA, NSPIN, . ISA, IPHORB, INDXUO, LASTO, XA, CELL, . RPSI, IPSI, E, INDW, NWF, NUMWF, NK, K, . IDIMEN, IOPTION, XMIN, XMAX, YMIN, YMAX, . ZMIN, ZMAX, NPX, NPY, NPZ, COORPO, NORMAL, . DIRVER1, DIRVER2, . ARMUNI, IUNITCD, ISCALE, RMAXO ) C ********************************************************************** C Compute the wave functions at the points of a plane or a 3D grid C in real space C Coded by P. Ordejon, from Junquera's rhoofr. July 2003 C ********************************************************************** use precision USE FDF USE ATMFUNCS USE CHEMICAL USE LISTSC_MODULE, ONLY: LISTSC use planed, only: plane IMPLICIT NONE INTEGER, INTENT(IN) :: . NA, NO, NUO, IOPTION, NPX, NPY, NPZ, ISCALE, IUNITCD, . IDIMEN, NSPIN, MAXNA, NK, NWF(NK), NUMWF, . ISA(NA), IPHORB(NO), INDXUO(NO), LASTO(0:NA), . INDW(NK,NUMWF) real(dp), INTENT(IN) :: . XA(3,NA), XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, . COORPO(3,3), NORMAL(3), DIRVER1(3), DIRVER2(3), ARMUNI, . RMAXO real(dp), INTENT(IN) :: . CELL(3,3), . RPSI(NUO,NK,NUMWF,NSPIN), IPSI(NUO,NK,NUMWF,NSPIN), . E(NK,NUMWF,NSPIN), K(NK,3) C ****** INPUT ********************************************************* C INTEGER NA : Total number of atoms in Supercell C INTEGER NO : Total number of orbitals in Supercell C INTEGER NUO : Total number of orbitals in Unit Cell C INTEGER MAXNA : Maximum number of neighbours of any atom C INTEGER NSPIN : Number of different spin polarizations C Nspin = 1 => unpolarized, Nspin = 2 => polarized C INTEGER ISA(NA) : Species index of each atom C INTEGER IPHORB(NO) : Orital index of each orbital in its atom C INTEGER INDXUO(NO) : Equivalent orbital in unit cell C INTEGER LASTO(0:NA) : Last orbital of each atom in array iphorb C REAL*8 XA(3,NA) : Atomic positions in cartesian coordinates C (in bohr) C REAL*8 CELL(3,3) : Supercell vectors CELL(IXYZ,IVECT) C (in bohr) C REAL*8 RPSI(NUO,NK,NUMWF,NSPIN): Wave function coefficients (real part) C REAL*8 IPSI(NUO,NK,NUMWF,NSPIN): Wave function coefficients (imag part) C REAL*8 E(NK,NUMWF,NSPIN) : Eigen energies C INTEGER INDW(NUMWF,NK) : Index of the wavefunctions C INTEGER NWF(NK) : Number of wavefncts to print for each k-point C INTEGER NUMWF : Max num of wavefncts to print a given k-point C INTEGER NK : Number of k-points C REAL*8 K(NK,3) : k-points C INTEGER IDIMEN : Specify if the run is to plot quantities C in a plane or in a 3D grid (2 or 3, respect) C INTEGER IOPTION : Option to generate the plane C 1 = Normal vector C 2 = Two vectors contained in the plane C 3 = Three points in the plane C INTEGER NPX,NPY,NPZ : Number of points along x and y and z C REAL*8 XMIN, XMAX : Limits of the plane in the PRF for x-direction C REAL*8 YMIN, YMAX : Limits of the plane in the PRF for y-direction C REAL*8 ZMIN, ZMAX : Limits of the z-direction C REAL*8 COORPO(3,3) : Coordinates of the three points used to define C the plane (Only used if ioption = 3) C REAL*8 NORMAL(3) : Components of the normal vector used to define C the plane C REAL*8 DIRVER1(3) : Components of the first vector contained C in the plane C (Only used if ioption = 2) C REAL*8 DIRVER2(3) : Components of the second vector contained C in the plane C (Only used if ioption = 2) C REAL*8 ARMUNI : Conversion factor for the charge density C INTEGER IUNITCD : Unit of the charge density C INTEGER ISCALE : Unit if the points of the plane C REAL*8 RMAXO : Maximum range of basis orbitals C ********************************************************************** INTEGER . NPLAMAX, NAPLA, NAINCELL INTEGER, DIMENSION(:), ALLOCATABLE :: . INDICES, JNA REAL, DIMENSION(:,:), allocatable :: RWF, IMWF, MWF, PWF real(dp), DIMENSION(:), ALLOCATABLE :: . R2IJ real(dp), DIMENSION(:,:), ALLOCATABLE :: . PLAPO, POINRE, XAPLA, XIJ, XAINCELL INTEGER . NPO, IA, ISEL, NNA, UNITRE1, UNITRE2, UNITIM1, UNITIM2, . UNITPH1, UNITPH2, UNITMO1, UNITMO2 INTEGER . I, J, IN, IAT1, IAT2, JO, IO, IUO, IAVEC, IAVEC1, IAVEC2, . IS1, IS2, IPHI1, IPHI2, IND, IX, IY, IZ, NX, NY, NZ, IWF, . INDWF, IZA(NA), IK real(dp) . RMAX, XPO(3), RMAX2, XVEC1(3), . XVEC2(3), PHIMU, GRPHIMU(3), . OCELL(3,3), PHASE, SI, CO, PI real(dp) . RWAVE, RWAVEUP, RWAVEDN, . IWAVE, IWAVEUP, IWAVEDN, . MWAVE, MWAVEUP, MWAVEDN, . PWAVE, PWAVEUP, PWAVEDN LOGICAL FIRST CHARACTER . SNAME*40, FNAMEWFRE*60, FNAMEWFIM*60, . FNAMEWFURE*60, FNAMEWFUIM*60, FNAMEWFDRE*60, FNAMEWFDIM*60, . FNAMEWFMO*60, FNAMEWFPH*60, . FNAMEWFUMO*60, FNAMEWFUPH*60, FNAMEWFDMO*60, FNAMEWFDPH*60, . PASTE*60, CHAR1*10, CHAR2*10, ITOCHAR*10, . EXT*20, EXT2*25 EXTERNAL . IO_ASSIGN, IO_CLOSE, PASTE, . NEIGHB, WROUT, ITOCHAR C ********************************************************************** C INTEGER NPLAMAX : Maximum number of points in the plane C REAL*8 PLAPO(NPLAMAX,3) : Coordinates of the points of the plane in PRF C REAL*8 POINRE(NPLAMAX,3): Coordinates of the points of the plane in Lattice C Reference Frame C INTEGER NAPLA : Number of atoms whose coordinates will be rotated C INTEGER INDICES(NA) : Indices of the atoms whose coordinates will C be rotated from the lattice reference frame C to the in-plane reference frame C REAL*8 XAPLA(3,NA) : Atomic coordinates in plane reference frame C INTEGER IA : Atom whose neighbours are needed. C A routine initialization must be done by C a first call with IA = 0 C If IA0=0, point X0 is used as origin instead C INTEGER ISEL : Single-counting switch (0=No, 1=Yes). If ISEL=1, C only neighbours with JA.LE.IA are included in JNA C INTEGER NNA : Number of non-zero orbitals at a point in C real space C INTEGER JNA(MAXNA) : Atom index of neighbours. The neighbours C atoms might be in the supercell C REAL*8 XIJ(3,MAXNA) : Vectors from point in real space to orbitals C REAL*8 R2IJ(MAXNA) : Squared distance to atomic orbitals C REAL*8 XPO(3) : Coordinates of the point of the plane respect C we are going to calculate the neighbours orbitals C REAL RWF(NPO,NSPIN) : Wave fnctn at each point of the grid (real part) C REAL IMWF(NPO,NSPIN) : Wave fnctn at each point of the grid (imag part) C INTEGER IZA(NA) : Atomic number of each atom C ********************************************************************** C Allocate some variables --------------------------------------------- PI = 4.0D0 * ATAN(1.0D0) NPLAMAX = NPX * NPY * NPZ ALLOCATE(PLAPO(NPLAMAX,3)) CALL MEMORY('A','D',3*NPLAMAX,'wavofr') ALLOCATE(POINRE(NPLAMAX,3)) CALL MEMORY('A','D',3*NPLAMAX,'wavofr') ALLOCATE(INDICES(NA)) CALL MEMORY('A','I',NA,'wavofr') ALLOCATE(XAPLA(3,NA)) CALL MEMORY('A','D',3*NA,'wavofr') ALLOCATE(XAINCELL(3,NA)) CALL MEMORY('A','D',3*NA,'wavofr') C Build the plane ------------------------------------------------------ CALL PLANE( NA, NPLAMAX, IDIMEN, IOPTION, . XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, . NPX, NPY, NPZ, COORPO, NORMAL, . DIRVER1, DIRVER2, . XA, NAPLA, INDICES, ISCALE, . POINRE, PLAPO, XAPLA ) C End building of the plane -------------------------------------------- C Initialize neighbour subroutine -------------------------------------- IA = 0 ISEL = 0 RMAX = RMAXO NNA = MAXNA IF (ALLOCATED(JNA)) THEN CALL MEMORY('D','I',SIZE(JNA),'wavofr') DEALLOCATE(JNA) ENDIF IF (ALLOCATED(R2IJ)) THEN CALL MEMORY('D','D',SIZE(R2IJ),'wavofr') DEALLOCATE(R2IJ) ENDIF IF (ALLOCATED(XIJ)) THEN CALL MEMORY('D','D',SIZE(XIJ),'wavofr') DEALLOCATE(XIJ) ENDIF ALLOCATE(JNA(MAXNA)) CALL MEMORY('A','I',MAXNA,'wavofr') ALLOCATE(R2IJ(MAXNA)) CALL MEMORY('A','D',MAXNA,'wavofr') ALLOCATE(XIJ(3,MAXNA)) CALL MEMORY('A','D',3*MAXNA,'wavofr') FIRST = .TRUE. DO I = 1,3 XPO(I) = 0.D0 ENDDO CALL NEIGHB( CELL, RMAX, NA, XA, XPO, IA, ISEL, . NNA, JNA, XIJ, R2IJ, FIRST ) FIRST = .FALSE. RMAX2 = RMAXO**2 C Loop over wavefunctions to write out DO IK = 1, NK DO IWF = 1,NWF(IK) INDWF = INDW(IK,IWF) CHAR1 = ITOCHAR(INDWF) CHAR2 = ITOCHAR(IK) C Open files to store wave functions ----------------------------------- SNAME = FDF_STRING('SystemLabel','siesta') IF (NSPIN .EQ. 1) THEN IF (IDIMEN .EQ. 2) THEN FNAMEWFRE = PASTE(SNAME,'.CON.K') FNAMEWFRE = PASTE(FNAMEWFRE,CHAR2) FNAMEWFRE = PASTE(FNAMEWFRE,'.WF') FNAMEWFRE = PASTE(FNAMEWFRE,CHAR1) FNAMEWFPH = PASTE(FNAMEWFRE,'.PHASE') FNAMEWFMO = PASTE(FNAMEWFRE,'.MOD') FNAMEWFIM = PASTE(FNAMEWFRE,'.IMAG') FNAMEWFRE = PASTE(FNAMEWFRE,'.REAL') ELSEIF (IDIMEN .EQ. 3) THEN FNAMEWFRE = PASTE(SNAME,'.K') FNAMEWFRE = PASTE(FNAMEWFRE,CHAR2) FNAMEWFRE = PASTE(FNAMEWFRE,'.WF') FNAMEWFRE = PASTE(FNAMEWFRE,CHAR1) FNAMEWFPH = PASTE(FNAMEWFRE,'.PHASE') FNAMEWFPH = PASTE(FNAMEWFPH,'.cube') FNAMEWFMO = PASTE(FNAMEWFRE,'.MOD') FNAMEWFMO = PASTE(FNAMEWFMO,'.cube') FNAMEWFIM = PASTE(FNAMEWFRE,'.IMAG') FNAMEWFIM = PASTE(FNAMEWFIM,'.cube') FNAMEWFRE = PASTE(FNAMEWFRE,'.REAL') FNAMEWFRE = PASTE(FNAMEWFRE,'.cube') ENDIF CALL IO_ASSIGN(UNITRE1) OPEN(UNIT = UNITRE1, FILE = FNAMEWFRE, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITRE1) CALL IO_ASSIGN(UNITIM1) OPEN(UNIT = UNITIM1, FILE = FNAMEWFIM, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITIM1) CALL IO_ASSIGN(UNITMO1) OPEN(UNIT = UNITMO1, FILE = FNAMEWFMO, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITMO1) CALL IO_ASSIGN(UNITPH1) OPEN(UNIT = UNITPH1, FILE = FNAMEWFPH, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITPH1) ELSEIF (NSPIN .EQ. 2) THEN IF (IDIMEN .EQ. 2) THEN FNAMEWFURE = PASTE(SNAME,'.CON.K') FNAMEWFURE = PASTE(FNAMEWFURE,CHAR2) FNAMEWFURE = PASTE(FNAMEWFURE,'.WF') FNAMEWFURE = PASTE(FNAMEWFURE,CHAR1) FNAMEWFDRE = PASTE(FNAMEWFURE,'.DOWN') FNAMEWFDPH = PASTE(FNAMEWFDRE,'.PHASE') FNAMEWFDMO = PASTE(FNAMEWFDRE,'.MOD') FNAMEWFDIM = PASTE(FNAMEWFDRE,'.IMAG') FNAMEWFDRE = PASTE(FNAMEWFDRE,'.REAL') FNAMEWFURE = PASTE(FNAMEWFURE,'.UP') FNAMEWFUPH = PASTE(FNAMEWFURE,'.PHASE') FNAMEWFUMO = PASTE(FNAMEWFURE,'.MOD') FNAMEWFUIM = PASTE(FNAMEWFURE,'.IMAG') FNAMEWFURE = PASTE(FNAMEWFURE,'.REAL') ELSE IF (IDIMEN .EQ. 3) THEN FNAMEWFURE = PASTE(SNAME,'.K') FNAMEWFURE = PASTE(FNAMEWFURE,CHAR2) FNAMEWFURE = PASTE(FNAMEWFURE,'.WF') FNAMEWFURE = PASTE(FNAMEWFURE,CHAR1) FNAMEWFDPH = PASTE(FNAMEWFURE,'.DOWN.PHASE') FNAMEWFDPH = PASTE(FNAMEWFDPH,'.cube') FNAMEWFDMO = PASTE(FNAMEWFURE,'.DOWN.MOD') FNAMEWFDMO = PASTE(FNAMEWFDMO,'.cube') FNAMEWFDIM = PASTE(FNAMEWFURE,'.DOWN.IMAG') FNAMEWFDIM = PASTE(FNAMEWFDIM,'.cube') FNAMEWFDRE = PASTE(FNAMEWFURE,'.DOWN.REAL') FNAMEWFDRE = PASTE(FNAMEWFDRE,'.cube') FNAMEWFUPH = PASTE(FNAMEWFURE,'.UP.PHASE') FNAMEWFUPH = PASTE(FNAMEWFUPH,'.cube') FNAMEWFUMO = PASTE(FNAMEWFURE,'.UP.MOD') FNAMEWFUMO = PASTE(FNAMEWFUMO,'.cube') FNAMEWFUIM = PASTE(FNAMEWFURE,'.UP.IMAG') FNAMEWFUIM = PASTE(FNAMEWFUIM,'.cube') FNAMEWFURE = PASTE(FNAMEWFURE,'.UP.REAL') FNAMEWFURE = PASTE(FNAMEWFURE,'.cube') ENDIF CALL IO_ASSIGN(UNITRE1) OPEN(UNIT = UNITRE1, FILE = FNAMEWFURE, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITRE1) CALL IO_ASSIGN(UNITRE2) OPEN(UNIT = UNITRE2, FILE = FNAMEWFDRE, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITRE2) CALL IO_ASSIGN(UNITIM1) OPEN(UNIT = UNITIM1, FILE = FNAMEWFUIM, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITIM1) CALL IO_ASSIGN(UNITIM2) OPEN(UNIT = UNITIM2, FILE = FNAMEWFDIM, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITIM2) CALL IO_ASSIGN(UNITMO1) OPEN(UNIT = UNITMO1, FILE = FNAMEWFUMO, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITMO1) CALL IO_ASSIGN(UNITMO2) OPEN(UNIT = UNITMO2, FILE = FNAMEWFDMO, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITMO2) CALL IO_ASSIGN(UNITPH1) OPEN(UNIT = UNITPH1, FILE = FNAMEWFUPH, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITPH1) CALL IO_ASSIGN(UNITPH2) OPEN(UNIT = UNITPH2, FILE = FNAMEWFDPH, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNITPH2) ELSE WRITE(6,*)'BAD NUMBER NSPIN IN WAVOFR.F' WRITE(6,*)'NSPIN = ',NSPIN WRITE(6,*)'IT MUST BE 1 => NON POLARIZED, OR 2 = > POLARIZED' STOP ENDIF IF (IDIMEN .EQ. 2) THEN C Select all atoms in list to be printed out NAINCELL=NAPLA DO IA=1,NAPLA DO IX=1,3 XAINCELL(IX,IA)=XAPLA(IX,IA) ENDDO ENDDO ELSE IF (IDIMEN .EQ.3) THEN DO IX = 1,3 DO IY = 1,3 OCELL(IX,IY)=0.D0 ENDDO ENDDO C Determine cell size OCELL(1,1) = DABS(XMAX-XMIN) OCELL(2,2) = DABS(YMAX-YMIN) OCELL(3,3) = DABS(ZMAX-ZMIN) C Determine atoms which are within the plotting box NAINCELL=0 DO IA=1,NA IF ((XAPLA(1,IA).LT.XMIN*1.1).OR.(XAPLA(1,IA).GT.XMAX*1.1) . .OR. (XAPLA(2,IA).LT.YMIN*1.1).OR.(XAPLA(2,IA).GT.YMAX*1.1) . .OR. (XAPLA(3,IA).LT.ZMIN*1.1).OR.(XAPLA(3,IA).GT.ZMAX*1.1)) . GOTO 90 NAINCELL=NAINCELL+1 IZA(NAINCELL) = ATOMIC_NUMBER(ISA(IA)) DO IX=1,3 XAINCELL(IX,NAINCELL)=XAPLA(IX,IA) ENDDO 90 CONTINUE ENDDO IF (NSPIN .EQ. 1) THEN WRITE(UNITRE1,*) FNAMEWFRE WRITE(UNITRE1,*) FNAMEWFRE WRITE(UNITRE1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITRE1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITRE1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITRE1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITRE1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITIM1,*) FNAMEWFIM WRITE(UNITIM1,*) FNAMEWFIM WRITE(UNITIM1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITIM1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITIM1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITIM1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITIM1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITMO1,*) FNAMEWFMO WRITE(UNITMO1,*) FNAMEWFMO WRITE(UNITMO1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITMO1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITMO1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITMO1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITMO1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITPH1,*) FNAMEWFPH WRITE(UNITPH1,*) FNAMEWFPH WRITE(UNITPH1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITPH1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITPH1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITPH1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITPH1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO ELSE IF (NSPIN .EQ. 2) THEN WRITE(UNITRE1,*) FNAMEWFURE WRITE(UNITRE1,*) FNAMEWFURE WRITE(UNITRE1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITRE1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITRE1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITRE1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITRE1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITRE2,*) FNAMEWFDRE WRITE(UNITRE2,*) FNAMEWFDRE WRITE(UNITRE2,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITRE2,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITRE2,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITRE2,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITRE2,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITIM1,*) FNAMEWFUIM WRITE(UNITIM1,*) FNAMEWFUIM WRITE(UNITIM1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITIM1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITIM1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITIM1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITIM1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITIM2,*) FNAMEWFDIM WRITE(UNITIM2,*) FNAMEWFDIM WRITE(UNITIM2,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITIM2,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITIM2,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITIM2,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITIM2,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITMO1,*) FNAMEWFUMO WRITE(UNITMO1,*) FNAMEWFUMO WRITE(UNITMO1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITMO1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITMO1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITMO1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITMO1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITMO2,*) FNAMEWFDMO WRITE(UNITMO2,*) FNAMEWFDMO WRITE(UNITMO2,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITMO2,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITMO2,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITMO2,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITMO2,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITPH1,*) FNAMEWFUPH WRITE(UNITPH1,*) FNAMEWFUPH WRITE(UNITPH1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITPH1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITPH1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITPH1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITPH1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNITPH2,*) FNAMEWFDPH WRITE(UNITPH2,*) FNAMEWFDPH WRITE(UNITPH2,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNITPH2,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNITPH2,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNITPH2,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNITPH2,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO ENDIF ENDIF CALL WROUT(IDIMEN, .FALSE., .TRUE., IOPTION, NORMAL, COORPO, . DIRVER1, DIRVER2, . NPX, NPY, NPZ, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, . IUNITCD, NA, NAINCELL, INDICES, XAINCELL ) WRITE(6,'(A)') WRITE(6,'(A)') . ' Generating values of wavefunction on the grid' WRITE(6,'(A,I6)') . ' for k-point ',IK WRITE(6,'(A,I6)') . ' for eigenstate number ',INDW(IK,IWF) WRITE(6,'(A,I6)') . ' Please, wait....' WRITE(6,'(A)') C Allocate space for wave functions in 3D-grid IF (.NOT.ALLOCATED(RWF)) THEN ALLOCATE(RWF(NPX*NPY*NPZ,NSPIN)) CALL MEMORY('A','D',NPX*NPY*NPZ*NSPIN,'WAVOFR') ALLOCATE(IMWF(NPX*NPY*NPZ,NSPIN)) CALL MEMORY('A','D',NPX*NPY*NPZ*NSPIN,'WAVOFR') ALLOCATE(MWF(NPX*NPY*NPZ,NSPIN)) CALL MEMORY('A','D',NPX*NPY*NPZ*NSPIN,'WAVOFR') ALLOCATE(PWF(NPX*NPY*NPZ,NSPIN)) CALL MEMORY('A','D',NPX*NPY*NPZ*NSPIN,'WAVOFR') ENDIF C Loop over wavefunctions to write out C Loop over all points in real space ----------------------------------- C DO 100 NPO = 1, NPX*NPY*NPZ NPO = 0 DO 102 NZ = 1,NPZ DO 101 NY = 1,NPY DO 100 NX = 1,NPX NPO = NPO + 1 C Initialize the wave function at each point ----------------------- RWAVE = 0.D0 RWAVEUP = 0.D0 RWAVEDN = 0.D0 IWAVE = 0.D0 IWAVEUP = 0.D0 IWAVEDN = 0.D0 MWAVE = 0.D0 MWAVEUP = 0.D0 MWAVEDN = 0.D0 PWAVE = 0.D0 PWAVEUP = 0.D0 PWAVEDN = 0.D0 C Localize non-zero orbitals at each point in real space --------------- DO IX = 1,3 XPO(IX) = POINRE(NPO,IX) ENDDO IA = 0 ISEL = 0 NNA = MAXNA CALL NEIGHB( CELL, RMAX, NA, XA, XPO, IA, ISEL, . NNA, JNA, XIJ, R2IJ, FIRST ) C Loop over Non-zero orbitals ------------------------------------------ DO 110 IAT1 = 1, NNA IF( R2IJ(IAT1) .GT. RMAX2 ) EXIT IAVEC1 = JNA(IAT1) IS1 = ISA(IAVEC1) XVEC1(1) = -XIJ(1,IAT1) XVEC1(2) = -XIJ(2,IAT1) XVEC1(3) = -XIJ(3,IAT1) C We take atom 1 as center of coordinates C PHASE = K(IK,1)*XA(1,IAVEC1)+ C . K(IK,2)*XA(2,IAVEC1)+ C . K(IK,3)*XA(3,IAVEC1) PHASE = K(IK,1)*(XPO(1)+XIJ(1,IAT1))+ . K(IK,2)*(XPO(2)+XIJ(2,IAT1))+ . K(IK,3)*(XPO(3)+XIJ(3,IAT1)) SI=SIN(PHASE) CO=COS(PHASE) DO 120 IO = LASTO(IAVEC1-1) + 1, LASTO(IAVEC1) IPHI1 = IPHORB(IO) IUO = INDXUO(IO) CALL PHIATM( IS1, IPHI1, XVEC1, PHIMU, GRPHIMU ) IF ( NSPIN .EQ. 1 ) THEN RWAVE = RWAVE + PHIMU * . (RPSI(IUO,IK,IWF,1)*CO - IPSI(IUO,IK,IWF,1)*SI) IWAVE = IWAVE + PHIMU * . (RPSI(IUO,IK,IWF,1)*SI + IPSI(IUO,IK,IWF,1)*CO) ELSEIF (NSPIN .EQ. 2) THEN RWAVEUP = RWAVEUP + PHIMU * . (RPSI(IUO,IK,IWF,1)*CO - IPSI(IUO,IK,IWF,1)*SI) IWAVEUP = IWAVEUP + PHIMU * . (RPSI(IUO,IK,IWF,1)*SI + IPSI(IUO,IK,IWF,1)*CO) RWAVEDN = RWAVEDN + PHIMU * . (RPSI(IUO,IK,IWF,2)*CO - IPSI(IUO,IK,IWF,2)*SI) IWAVEDN = IWAVEDN + PHIMU * . (RPSI(IUO,IK,IWF,2)*SI + IPSI(IUO,IK,IWF,2)*CO) ENDIF C calculate module and phase of wavefunction IF ( NSPIN .EQ. 1 ) THEN MWAVE = DSQRT(RWAVE**2 + IWAVE**2) IF (DABS(IWAVE) .LT. 1.D-6) IWAVE = DABS(IWAVE) IF (DABS(RWAVE) .LT. 1.D-12) THEN IF (IWAVE .GT. 0.0D0) PWAVE = PI/2.0D0 IF (IWAVE .LT. 0.0D0) PWAVE = -PI/2.0D0 ELSE PWAVE = DATAN(IWAVE/RWAVE) ENDIF IF (RWAVE .LT. 0.0D0 .AND. IWAVE .GE. 0.0d0) . PWAVE = PWAVE + PI IF (RWAVE .LT. 0.0D0 .AND. IWAVE .LT. 0.0D0) . PWAVE = PWAVE - PI ELSEIF (NSPIN .EQ. 2) THEN MWAVEUP = DSQRT(RWAVEUP**2 + IWAVEUP**2) IF (DABS(IWAVEUP) .LT. 1.D-6) IWAVEUP = DABS(IWAVEUP) IF (DABS(RWAVEUP) .LT. 1.D-12) THEN IF (IWAVEUP .GT. 0.0D0) PWAVEUP = PI/2.0D0 IF (IWAVEUP .LT. 0.0D0) PWAVEUP = -PI/2.0D0 ELSE PWAVEUP = DATAN(IWAVEUP/RWAVEUP) ENDIF IF (RWAVEUP .LT. 0.0D0 .AND. IWAVEUP .GE. 0.0d0) . PWAVEUP = PWAVEUP + PI IF (RWAVEUP .LT. 0.0D0 .AND. IWAVEUP .LT. 0.0D0) . PWAVEUP = PWAVEUP - PI MWAVEDN = DSQRT(RWAVEDN**2 + IWAVEDN**2) IF (DABS(IWAVEDN) .LT. 1.D-6) IWAVEDN = DABS(IWAVEDN) IF (DABS(RWAVEDN) .LT. 1.D-12) THEN IF (IWAVEDN .GT. 0.0D0) PWAVEDN = PI/2.0D0 IF (IWAVEDN .LT. 0.0D0) PWAVEDN = -PI/2.0D0 ELSE PWAVEDN = DATAN(IWAVEDN/RWAVEDN) ENDIF IF (RWAVEDN .LT. 0.0D0 .AND. IWAVEDN .GE. 0.0d0) . PWAVEDN = PWAVEDN + PI IF (RWAVEDN .LT. 0.0D0 .AND. IWAVEDN .LT. 0.0D0) . PWAVEDN = PWAVEDN - PI ENDIF 120 ENDDO 110 ENDDO IF (IDIMEN .EQ. 2) THEN IF ( NSPIN .EQ. 1 ) THEN WRITE(UNITRE1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),RWAVE*SQRT(ARMUNI) WRITE(UNITIM1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),IWAVE*SQRT(ARMUNI) WRITE(UNITMO1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),MWAVE*SQRT(ARMUNI) WRITE(UNITPH1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),PWAVE ELSEIF ( NSPIN .EQ. 2 ) THEN WRITE(UNITRE1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),RWAVEUP*SQRT(ARMUNI) WRITE(UNITRE2,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),RWAVEDN*SQRT(ARMUNI) WRITE(UNITIM1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),IWAVEUP*SQRT(ARMUNI) WRITE(UNITIM2,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),IWAVEDN*SQRT(ARMUNI) WRITE(UNITMO1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),MWAVEUP*SQRT(ARMUNI) WRITE(UNITMO2,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),MWAVEDN*SQRT(ARMUNI) WRITE(UNITPH1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),PWAVEUP WRITE(UNITPH2,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),PWAVEDN ENDIF ELSE IF (IDIMEN .EQ. 3) THEN IF (NSPIN .EQ. 1) THEN RWF(NPO,1) = RWAVE*SQRT(ARMUNI) IMWF(NPO,1) = IWAVE*SQRT(ARMUNI) MWF(NPO,1) = MWAVE*SQRT(ARMUNI) PWF(NPO,1) = PWAVE ELSE IF (NSPIN .EQ.2) THEN RWF(NPO,1) = RWAVEUP*SQRT(ARMUNI) RWF(NPO,2) = RWAVEDN*SQRT(ARMUNI) IMWF(NPO,1) = IWAVEUP*SQRT(ARMUNI) IMWF(NPO,2) = IWAVEDN*SQRT(ARMUNI) MWF(NPO,1) = MWAVEUP*SQRT(ARMUNI) MWF(NPO,2) = MWAVEDN*SQRT(ARMUNI) PWF(NPO,1) = PWAVEUP PWF(NPO,2) = PWAVEDN ENDIF ENDIF IF (IDIMEN .EQ. 2) THEN C TODO: bug fixed here: NPX -> NPY IF ( MOD(NPO,NPY) .EQ. 0 ) THEN C TODO: bug fixed here WRITE(UNITRE1,*) WRITE(UNITIM1,*) WRITE(UNITPH1,*) WRITE(UNITMO1,*) IF ( NSPIN .EQ. 2 ) THEN C TODO: bug fixed here WRITE(UNITRE2,*) WRITE(UNITIM2,*) WRITE(UNITPH2,*) WRITE(UNITMO2,*) ENDIF ENDIF ENDIF C End x loop 100 ENDDO C End y and z loops 101 ENDDO 102 ENDDO IF (IDIMEN .EQ. 3) THEN IF (NSPIN .EQ. 1) THEN DO NX=1,NPX DO NY=1,NPY WRITE(UNITRE1,'(6e13.5)') . (RWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNITIM1,'(6e13.5)') . (IMWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNITMO1,'(6e13.5)') . (MWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNITPH1,'(6e13.5)') . (PWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) ENDDO ENDDO ELSE IF (NSPIN .EQ. 2) THEN DO NX=1,NPX DO NY=1,NPY WRITE(UNITRE1,'(6e13.5)') . (RWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNITRE2,'(6e13.5)') . (RWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ) WRITE(UNITIM1,'(6e13.5)') . (IMWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNITIM2,'(6e13.5)') . (IMWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ) WRITE(UNITMO1,'(6e13.5)') . (MWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNITMO2,'(6e13.5)') . (MWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ) WRITE(UNITPH1,'(6e13.5)') . (PWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNITPH2,'(6e13.5)') . (PWF(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ) ENDDO ENDDO ENDIF WRITE(6,'(A)') WRITE(6,'(A)') . ' Your output files are:' IF (NSPIN .EQ. 1) THEN WRITE(6,'(A,A)') ' ',FNAMEWFRE WRITE(6,'(A,A)') ' ',FNAMEWFIM WRITE(6,'(A,A)') ' ',FNAMEWFMO WRITE(6,'(A,A)') ' ',FNAMEWFPH ELSE IF (NSPIN .EQ. 2) THEN WRITE(6,'(A,A)') ' ',FNAMEWFURE WRITE(6,'(A,A)') ' ',FNAMEWFUIM WRITE(6,'(A,A)') ' ',FNAMEWFDRE WRITE(6,'(A,A)') ' ',FNAMEWFDIM WRITE(6,'(A,A)') ' ',FNAMEWFUMO WRITE(6,'(A,A)') ' ',FNAMEWFUPH WRITE(6,'(A,A)') ' ',FNAMEWFDMO WRITE(6,'(A,A)') ' ',FNAMEWFDPH ENDIF ENDIF CALL IO_CLOSE(UNITRE1) CALL IO_CLOSE(UNITIM1) IF (NSPIN .EQ. 2) CALL IO_CLOSE(UNITRE2) IF (NSPIN .EQ. 2) CALL IO_CLOSE(UNITIM2) CALL IO_CLOSE(UNITMO1) CALL IO_CLOSE(UNITPH1) IF (NSPIN .EQ. 2) CALL IO_CLOSE(UNITMO2) IF (NSPIN .EQ. 2) CALL IO_CLOSE(UNITPH2) ENDDO ENDDO CALL MEMORY('D','D',SIZE(RWF),'wavofr') DEALLOCATE(RWF) CALL MEMORY('D','D',SIZE(IMWF),'wavofr') DEALLOCATE(IMWF) CALL MEMORY('D','D',SIZE(MWF),'wavofr') DEALLOCATE(MWF) CALL MEMORY('D','D',SIZE(PWF),'wavofr') DEALLOCATE(PWF) RETURN END