Dear Tristana and Siesta Users, I want to ask if my problem will be fixed with the recompilation of denchar program using files you attached in the previous mail.
I have successfully computed wave functions in 3D mode, but I'm interested in some specific planes of my system. With the same files, and, obviously changing correctly the input_denchar.fdf file, I tried to compute wave functions in 2D mode. First, I found a segmentation fault error, but now (with the same files as before) I'm finding the following error: "/gpfs/apps/SIESTA/SRC/siesta-3.1/Util/Denchar/Src/../../../Src/fdf/fdf.f", line 792: 1525-038 The REWIND statement on unit 5 cannot be processed because an err or with errno 29 (Illegal seek) was encountered. The program will stop. srun: error: s34c2b09: task 0: Exited with exit code 1 Here are more details: http://www.mail-archive.com/[email protected]/msg04776.html I would be really thankful if I can fix this issue (I can share my files to see if they are correct). Thanks again, Nuria Núria García Castelló Departament d'Electrònica. Laboratori de Síntesi (Planta 1) Facultat de Física. Universitat de Barcelona. Martí i Franqués 1. 08028 Barcelona. Telf: (0034)934034804 Fax: (0034)934021148 Dear Ejaz: It is difficult to reproduce your bug without having all the files you used to run Denchar, but there was a bug corrected in version trunk-385 that *might* be related to your issue. It is not yet incorporated to the stable version of SIESTA (I assume you are using version 3.1). You might try to recompile denchar utility using the two source files that I sent you atached to this email. One should solve your problem, the other one fixes another (irrelevant) issue. If this still doesn't solve your issue, It would be nice to have all your imput files to find and fix the bug properly. Let me know if this solved your problem. Cheers Tris Dr. Tristana Sondon Institut de Ciencia de Materials de Barcelona (ICMAB-CSIC) Campus de la UAB 08193 Bellaterra (Barcelona), Spain ----- Mensaje original ----- De: "Khan Muhammad" <[email protected]> Para: [email protected] Enviados: Martes, 15 de Mayo 2012 11:32:04 Asunto: [SIESTA-L] Problem regarding post processing for wavefunction at specific energy using DENCHAR Dear J. Junquera and P. Ordejón, I have generated the band projected wavafunction using Siesta and then used DENCHAR for post processing and its worked perfectly. However when I tried to plot the wavefunction for small energy interval: WFS.EnergyMin -7.72 eV WFS.EnergyMax -7.62 eV the wavefunction file is generated, but Denchar produced the following error: Species number: 1 Label: C Atomic number: 6 Species number: 2 Label: H Atomic number: 1 forrtl: severe (153): allocatable array or pointer is not allocated Image PC Routine Line Source denchar 000000000053762D Unknown Unknown Unknown denchar 0000000000536135 Unknown Unknown Unknown denchar 00000000004E3789 Unknown Unknown Unknown denchar 0000000000497CBD Unknown Unknown Unknown denchar 00000000004C8161 Unknown Unknown Unknown denchar 000000000047ABD6 Unknown Unknown Unknown denchar 000000000047E207 Unknown Unknown Unknown denchar 000000000040361C Unknown Unknown Unknown libc.so.6 0000003743A1D8B4 Unknown Unknown Unknown denchar 0000000000403529 Unknown Unknown Unknown The input for denchar is as follow: SystemLabel cnt Denchar.TypeOfRun 3D Denchar.PlotCharge .FALSE. Denchar.PlotWaveFunctions .TRUE. Denchar.CoorUnits Ang Denchar.DensityUnits Ele/Ang**3 Denchar.NumberPointsX 50 Denchar.NumberPointsY 50 Denchar.NumberPointsZ 50 Denchar.MinX -10.000 Ang Denchar.MinY -10.000 Ang Denchar.MaxY 10.000 Ang Denchar.MinZ -20.000 Ang Denchar.MaxZ 20.000 Ang Denchar.PlaneGeneration NormalVector %block Denchar.CompNormalVector 0.0 0.0 1.0 %endblock Denchar.CompNormalVector %block Denchar.PlaneOrigin 0.0 0.0 0.0 %endblock Denchar.PlaneOrigin Please, help me to fix this problem. Thanking in anticipation With Best Regards Ejaz Khan PhD Candidate Graduate School of Nanoscience and Technology, Korea Advanced Institute of Science and Technology, (KAIST) SUBROUTINE ATOMPLA( NA, ORIGIN, XA, MROT, IDIMEN, ISCALE, . NATINPLA, INDICES, XAPLANE ) C ********************************************************************** C Find the coordinates of some selected atoms on the C plane-reference-frame. Atoms in plane will have the third coordinate C equal to zero. C For 3D-grid, all the atoms are selected C Coded by J. Junquera May '99 C Modified by P. Ordejon to include 3D capability; June 2003 C ********************************************************************** use precision, only: dp USE FDF, only: block_fdf, fdf_block, fdf_bline use fdf, only: fdf_bmatch, fdf_bintegers USE PARSE, only: parsed_line, digest, match USE SYS, only: die IMPLICIT NONE INTEGER . NA, NATINPLA, INDICES(NA), ISCALE real(dp) . ORIGIN(3), XA(3,NA), MROT(3,3), XAPLANE(3,NA) C ******* INPUT ******************************************************** C INTEGER NA : Number of atoms C REAL*8 ORIGIN(3) : Origin of the plane reference frame C REAL*8 XA(3,NA) : Atomic coordinates in lattice reference frame C REAL*8 MROT(3,3) : Rotation matrix from the lattice reference frame C to the in-plane reference frame C INTEGER IDIMEN : Determines if run is plane or 3D-grid C INTEGER ISCALE : Units of the atomic coordinates C ISCALE = 1 => bohrs, ISCALE = 2 => Ang C ******* OUTPUT ******************************************************* C INTEGER NATINPLA : Number of atoms whose coordinates will be rotated C INTEGER INDICES(NA) : Atomic indices of the atoms whose coordinates C will be rotated C REAL*8 XAPLANE(3,NA) : Atomic coordinates in plane reference frame C ********************************************************************** INTEGER . NATINPL_DEFECT, IDIMEN, IUNIT, IAT, IX real(dp) . VAUX1(3), VAUX2(3) LOGICAL . ATINPLA TYPE(PARSED_LINE), pointer :: line TYPE(BLOCK_fdf) :: BP EXTERNAL . MATVECT C ********************************************************************** C INTEGER NATINPLA : Number of atoms whose coordinates will be C rotated from the lattice reference frame to C the in-plane reference frame C REAL*8 VAUX1(3) : Auxiliar vector C REAL*8 VAUX2(3) : Auxiliar vector C LOGICAL ATINPLA : Is block Denchar.AtomsInPlane present in fdf? C ********************************************************************** IF (IDIMEN .EQ. 2) THEN C Read fdf data block 'Denchar.AtomsInPlane' --------------------------- NATINPL_DEFECT = 0 NATINPLA = 0 IF ( .NOT. FDF_BLOCK('Denchar.AtomsInPlane',BP) ) GOTO 2000 LOOP: DO IF (.NOT. FDF_BLINE(BP,LINE)) EXIT LOOP IF (.NOT. fdf_bmatch(line,"I") ) . CALL DIE("Wrong format in Denchar.AtomsInPlane") NATINPLA = NATINPLA + 1 INDICES(NATINPLA) = fdf_bintegers(line,1) ENDDO LOOP 2000 CONTINUE ELSE IF (IDIMEN .EQ. 3) THEN NATINPLA = NA DO IAT = 1, NA INDICES(IAT) = IAT ENDDO ELSE CALL DIE("Wrong IDIMEN in ATOMPLA") ENDIF C Rotate the coordinates ----------------------------------------------- DO IAT = 1, NATINPLA DO IX = 1,3 VAUX1(IX) = XA(IX,INDICES(IAT)) - ORIGIN(IX) ENDDO CALL MATVECT(MROT, VAUX1, VAUX2) DO IX = 1,3 XAPLANE(IX,INDICES(IAT)) = VAUX2(IX) ENDDO C IF (ISCALE .EQ. 2) THEN C DO IX = 1,3 C XAPLANE(IX,INDICES(IAT))=XAPLANE(IX,INDICES(IAT))*0.529177D0 C ENDDO C ENDIF ENDDO RETURN END 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 IF ( MOD(NPO,NPX) .EQ. 0 ) THEN CC-AG Use * for single line separation, CC-AG since (/) gives two... 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) !==TS== bug fixed for non spin-polarized cases IF (NSPIN .EQ. 2) THEN CALL IO_CLOSE(UNIT3) ENDIF !==TS== IF (IDIMEN .EQ. 2 .AND. NSPIN .EQ. 2) CALL IO_CLOSE(UNIT4) RETURN END
