SUBROUTINE RHOOFR( NA, NO, NUO, MAXND, MAXNA, NSPIN, . ISA, IPHORB, INDXUO, LASTO, . XA, CELL, NUMD, LISTD, LISTDPTR, DSCF, DATM, . IDIMEN, IOPTION, XMIN, XMAX, YMIN, YMAX, . ZMIN, ZMAX, NPX, NPY, NPZ, COORPO, NORMAL, . DIRVER1, DIRVER2, . ARMUNI, IUNITCD, ISCALE, RMAXO ) C ********************************************************************** C Compute the density of charge at the points of a plane or a 3D grid C in real space C Coded by J. Junquera November'98 C Modified by P. Ordejon to include 3D capabilities, June 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, MAXND, NSPIN, MAXNA, . ISA(NA), IPHORB(NO), INDXUO(NO), LASTO(0:NA), . NUMD(NUO), LISTDPTR(NUO), LISTD(MAXND) 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), DSCF(MAXND,NSPIN), . DATM(NO) 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 MAXND : Maximum number C of basis orbitals interacting, either directly C or through a KB projector, with any orbital 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 INTEGER NUMD(NUO) : Control vector of Density Matrix C (number of non-zero elements of eaach row) C INTEGER LISTDPTR(NUO) : Pointer to where each row of listh starts - 1 C The reason for pointing to the element before C the first one is so that when looping over the C elements of a row there is no need to shift by C minus one. C INTEGER LISTD(MAXND) : Control vector of Density Matrix C (list of non-zero elements of each row) C REAL*8 DSCF(MAXND,NSPIN): Density Matrix C REAL*8 DATM(NO) : Neutral atom charge density of each orbital 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 : Number of points generated along x and y C direction in a system of reference in which C the third components od the points of the plane is C zero (Plane Reference Frame; PRF) 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 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 COORPO(3,3) : Coordinates of the three points used to define C the plane (Only used if ioption = 3) C INTEGER IUNITCD : Unit of the charge density C INTEGER ISCALE : Unit if the points of the plane C REAL*8 ARMUNI : Conversion factor for the charge density C REAL*8 RMAXO : Maximum range of basis orbitals C ********************************************************************** INTEGER . NPLAMAX, NAPLA, NAINCELL INTEGER, DIMENSION(:), ALLOCATABLE :: . INDICES, JNA REAL, DIMENSION(:,:), allocatable :: RHO REAL, DIMENSION(:), allocatable :: DRHO real(dp), DIMENSION(:), ALLOCATABLE :: . R2IJ, DISCF, DIATM, DIUP, DIDOWN real(dp), DIMENSION(:,:), ALLOCATABLE :: . PLAPO, POINRE, XAPLA, XIJ, XAINCELL INTEGER . NPO, IA, ISEL, NNA, UNIT1, UNIT2, UNIT3, UNIT4 INTEGER . I, J, IN, IAT1, IAT2, JO, IO, IUO, IAVEC, IAVEC1, IAVEC2, . IS1, IS2, IPHI1, IPHI2, IND, IX, IY, IZ, NX, NY, NZ, IZA(NA) real(dp) . RMAX, XPO(3), RMAX2, XVEC1(3), . XVEC2(3), PHIMU, PHINU, GRPHIMU(3), GRPHINU(3), . OCELL(3,3) real(dp) . DENCHAR, DENCHAR0, DENUP, DENDOWN LOGICAL FIRST CHARACTER . SNAME*30, FNAMESCF*38, FNAMEDEL*38, FNAMEUP*38, FNAMEDOWN*38, . PASTE*38 EXTERNAL . IO_ASSIGN, IO_CLOSE, PASTE . NEIGHB, WROUT 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*8 DENCHAR : Self-Consistent Density of Charge at a given C point in real space C REAL*8 DENCHAR0 : Harris Density of Charge at a given point C in real space C INTEGER IZA(NA) : Atomic number of each atom C ********************************************************************** C Allocate some variables --------------------------------------------- NPLAMAX = NPX * NPY * NPZ ALLOCATE(PLAPO(NPLAMAX,3)) CALL MEMORY('A','D',3*NPLAMAX,'rhoofr') ALLOCATE(POINRE(NPLAMAX,3)) CALL MEMORY('A','D',3*NPLAMAX,'rhoofr') ALLOCATE(INDICES(NA)) CALL MEMORY('A','I',NA,'rhoofr') ALLOCATE(XAPLA(3,NA)) CALL MEMORY('A','D',3*NA,'rhoofr') ALLOCATE(XAINCELL(3,NA)) CALL MEMORY('A','D',3*NA,'rhoofr') ALLOCATE(DISCF(NO)) CALL MEMORY('A','D',NO,'rhoofr') ALLOCATE(DIATM(NO)) CALL MEMORY('A','D',NO,'rhoofr') ALLOCATE(DIUP(NO)) CALL MEMORY('A','D',NO,'rhoofr') ALLOCATE(DIDOWN(NO)) CALL MEMORY('A','D',NO,'rhoofr') 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),'rhoofr') DEALLOCATE(JNA) ENDIF IF (ALLOCATED(R2IJ)) THEN CALL MEMORY('D','D',SIZE(R2IJ),'rhoofr') DEALLOCATE(R2IJ) ENDIF IF (ALLOCATED(XIJ)) THEN CALL MEMORY('D','D',SIZE(XIJ),'rhoofr') DEALLOCATE(XIJ) ENDIF ALLOCATE(JNA(MAXNA)) CALL MEMORY('A','I',MAXNA,'rhoofr') ALLOCATE(R2IJ(MAXNA)) CALL MEMORY('A','D',MAXNA,'rhoofr') ALLOCATE(XIJ(3,MAXNA)) CALL MEMORY('A','D',3*MAXNA,'rhoofr') 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 Open files to store charge density ----------------------------------- SNAME = FDF_STRING('SystemLabel','siesta') IF (NSPIN .EQ. 1) THEN IF (IDIMEN .EQ. 2) THEN FNAMESCF = PASTE(SNAME,'.CON.SCF') FNAMEDEL = PASTE(SNAME,'.CON.DEL') CALL IO_ASSIGN(UNIT1) OPEN(UNIT = UNIT1, FILE = FNAMESCF, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT1) CALL IO_ASSIGN(UNIT2) OPEN(UNIT = UNIT2, FILE = FNAMEDEL, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT2) ELSEIF (IDIMEN .EQ. 3) THEN FNAMESCF = PASTE(SNAME,'.RHO.cube') FNAMEDEL = PASTE(SNAME,'.DRHO.cube') CALL IO_ASSIGN(UNIT1) OPEN(UNIT = UNIT1, FILE = FNAMESCF, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT1) CALL IO_ASSIGN(UNIT2) OPEN(UNIT = UNIT2, FILE = FNAMEDEL, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT2) ENDIF ELSEIF (NSPIN .EQ. 2) THEN IF (IDIMEN .EQ. 2) THEN FNAMESCF = PASTE(SNAME,'.CON.MAG' ) FNAMEDEL = PASTE(SNAME,'.CON.DEL' ) FNAMEUP = PASTE(SNAME,'.CON.UP' ) FNAMEDOWN= PASTE(SNAME,'.CON.DOWN') CALL IO_ASSIGN(UNIT1) OPEN(UNIT = UNIT1, FILE = FNAMESCF, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT1) CALL IO_ASSIGN(UNIT2) OPEN(UNIT = UNIT2, FILE = FNAMEDEL, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT2) CALL IO_ASSIGN(UNIT3) OPEN(UNIT = UNIT3, FILE = FNAMEUP, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT3) CALL IO_ASSIGN(UNIT4) OPEN(UNIT = UNIT4, FILE = FNAMEDOWN, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT4) ELSE IF (IDIMEN .EQ. 3) THEN FNAMEUP = PASTE(SNAME,'.RHO.UP.cube' ) FNAMEDOWN = PASTE(SNAME,'.RHO.DOWN.cube' ) FNAMEDEL = PASTE(SNAME,'.DRHO.cube' ) CALL IO_ASSIGN(UNIT1) OPEN(UNIT = UNIT1, FILE = FNAMEUP, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT1) CALL IO_ASSIGN(UNIT2) OPEN(UNIT = UNIT2, FILE = FNAMEDOWN, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT2) CALL IO_ASSIGN(UNIT3) OPEN(UNIT = UNIT3, FILE = FNAMEDEL, STATUS = 'UNKNOWN', . FORM = 'FORMATTED') REWIND(UNIT3) ENDIF ELSE WRITE(6,*)'BAD NUMBER NSPIN IN RHOOFR.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 C so that only they will be printed 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(UNIT1,*) FNAMESCF WRITE(UNIT1,*) FNAMESCF WRITE(UNIT1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNIT1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNIT1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNIT1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNIT1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNIT2,*) FNAMEDEL WRITE(UNIT2,*) FNAMEDEL WRITE(UNIT2,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNIT2,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNIT2,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNIT2,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNIT2,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO ELSE IF (NSPIN .EQ. 2) THEN WRITE(UNIT1,*) FNAMEUP WRITE(UNIT1,*) FNAMEUP WRITE(UNIT1,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNIT1,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNIT1,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNIT1,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNIT1,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNIT2,*) FNAMEDOWN WRITE(UNIT2,*) FNAMEDOWN WRITE(UNIT2,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNIT2,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNIT2,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNIT2,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNIT2,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO WRITE(UNIT3,*) FNAMEDEL WRITE(UNIT3,*) FNAMEDEL WRITE(UNIT3,'(i5,4f12.6)') NAINCELL, XMIN, YMIN, ZMIN WRITE(UNIT3,'(i5,4f12.6)') NPX,(OCELL(1,J)/(NPX-1),J=1,3) WRITE(UNIT3,'(i5,4f12.6)') NPY,(OCELL(2,J)/(NPY-1),J=1,3) WRITE(UNIT3,'(i5,4f12.6)') NPZ,(OCELL(3,J)/(NPZ-1),J=1,3) DO IA = 1,NAINCELL WRITE(UNIT3,'(i5,4f12.6)') IZA(IA),0.0, . (XAINCELL(IX,IA),IX=1,3) ENDDO ENDIF ENDIF CALL WROUT(IDIMEN, .TRUE., .FALSE., 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 Charge Density values on the grid now....' WRITE(6,'(A)') . ' Please, wait....' WRITE(6,'(A)') C Allocate space for density in 3D-grid ALLOCATE(RHO(NPX*NPY*NPZ,NSPIN)) CALL MEMORY('A','S',NPX*NPY*NPZ*NSPIN,'RHOOFR') ALLOCATE(DRHO(NPX*NPY*NPZ)) CALL MEMORY('A','S',NPX*NPY*NPZ,'RHOOFR') 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 density of charge at each point ----------------------- DENCHAR = 0.D0 DENCHAR0 = 0.D0 DENUP = 0.D0 DENDOWN = 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) DO 120 IO = LASTO(IAVEC1-1) + 1, LASTO(IAVEC1) IPHI1 = IPHORB(IO) IUO = INDXUO(IO) CALL PHIATM( IS1, IPHI1, XVEC1, PHIMU, GRPHIMU ) C Copy full row IUO of Density Matrix to DI(j) -------------------------- IF (IO . EQ. IUO) THEN DO 130 IN = 1, NUMD(IUO) IND = IN + LISTDPTR(IUO) J = LISTD(IND) IF (NSPIN .EQ. 1) THEN DISCF(J) = DSCF(IND,1) ELSEIF (NSPIN .EQ. 2) THEN DIUP(J) = DSCF(IND,1) DIDOWN(J) = DSCF(IND,2) ENDIF 130 ENDDO ELSE DO 135 IN = 1, NUMD(IUO) IND = IN + LISTDPTR(IUO) J = LISTSC( IO, IUO, LISTD(IND) ) IF (NSPIN .EQ. 1) THEN DISCF(J) = DSCF(IND,1) ELSEIF (NSPIN .EQ. 2) THEN DIUP(J) = DSCF(IND,1) DIDOWN(J) = DSCF(IND,2) ENDIF 135 ENDDO ENDIF DENCHAR0 = DENCHAR0 + PHIMU*PHIMU*DATM(IUO) C WRITE(6,*) IO,DATM(IO) C Loop again over neighbours ------------------------------------------- DO 140 IAT2 = 1, NNA IAVEC2 = JNA(IAT2) IS2 = ISA(IAVEC2) XVEC2(1) = -XIJ(1,IAT2) XVEC2(2) = -XIJ(2,IAT2) XVEC2(3) = -XIJ(3,IAT2) DO 150 JO = LASTO(IAVEC2-1) + 1, LASTO(IAVEC2) IPHI2 = IPHORB(JO) CALL PHIATM( IS2, IPHI2, XVEC2, PHINU, GRPHINU ) IF ( NSPIN .EQ. 1 ) THEN DENCHAR = DENCHAR + PHINU*PHIMU*DISCF(JO) ELSEIF (NSPIN .EQ. 2) THEN DENUP = DENUP + PHINU*PHIMU*DIUP(JO) DENDOWN = DENDOWN + PHINU*PHIMU*DIDOWN(JO) ENDIF 150 ENDDO 140 ENDDO 120 ENDDO 110 ENDDO IF (IDIMEN .EQ. 2) THEN IF ( NSPIN .EQ. 1 ) THEN WRITE(UNIT1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),DENCHAR*ARMUNI WRITE(UNIT2,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),(DENCHAR-DENCHAR0)*ARMUNI ELSEIF ( NSPIN .EQ. 2 ) THEN WRITE(UNIT1,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),(DENUP-DENDOWN)*ARMUNI WRITE(UNIT2,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2), . (DENUP+DENDOWN-DENCHAR0)*ARMUNI WRITE(UNIT3,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),DENUP*ARMUNI WRITE(UNIT4,'(3F12.5)') . PLAPO(NPO,1),PLAPO(NPO,2),DENDOWN*ARMUNI ENDIF ELSE IF (IDIMEN .EQ. 3) THEN IF (NSPIN .EQ. 1) THEN RHO(NPO,1) = DENCHAR*ARMUNI DRHO(NPO) = (DENCHAR-DENCHAR0)*ARMUNI ELSE IF (NSPIN .EQ.2) THEN RHO(NPO,1) = DENUP*ARMUNI RHO(NPO,2) = DENDOWN*ARMUNI DRHO(NPO) = (DENUP+DENDOWN-DENCHAR0)*ARMUNI ENDIF ENDIF IF (IDIMEN .EQ. 2) THEN C TODO: bug fixed here: NPX -> NPY IF ( MOD(NPO,NPY) .EQ. 0 ) THEN CC-AG WRITE(UNIT1,'(/)') CC-AG WRITE(UNIT2,'(/)') WRITE(UNIT1,*) WRITE(UNIT2,*) IF ( NSPIN .EQ. 2 ) THEN WRITE(UNIT3,'(/)') WRITE(UNIT4,'(/)') 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(UNIT1,'(6e13.5)') . (RHO(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNIT2,'(6e13.5)') . (DRHO(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY),NZ=1,NPZ) ENDDO ENDDO ELSE IF (NSPIN .EQ. 2) THEN DO NX=1,NPX DO NY=1,NPY WRITE(UNIT1,'(6e13.5)') . (RHO(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,1),NZ=1,NPZ) WRITE(UNIT2,'(6e13.5)') . (RHO(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY,2),NZ=1,NPZ) WRITE(UNIT3,'(6e13.5)') . (DRHO(NX+(NY-1)*NPX+(NZ-1)*NPX*NPY),NZ=1,NPZ) ENDDO ENDDO ENDIF ENDIF WRITE(6,'(A)') WRITE(6,'(A)') . ' Your output files are:' IF (NSPIN .EQ. 1) THEN WRITE(6,'(A,A)') ' ',FNAMESCF WRITE(6,'(A,A)') ' ',FNAMEDEL ELSE IF (NSPIN .EQ. 2) THEN WRITE(6,'(A,A)') ' ',FNAMESCF WRITE(6,'(A,A)') ' ',FNAMEDEL WRITE(6,'(A,A)') ' ',FNAMEUP WRITE(6,'(A,A)') ' ',FNAMEDOWN ENDIF CALL IO_CLOSE(UNIT1) CALL IO_CLOSE(UNIT2) CALL IO_CLOSE(UNIT3) IF (IDIMEN .EQ. 2 .AND. NSPIN .EQ. 2) CALL IO_CLOSE(UNIT4) RETURN END