1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8 9 SUBROUTINE READSTS(VOLUME,EMIN, NE, DELTAE, ETA, XSTS, BFUNC, 10 . IUNITCD, ARMUNI) 11 12 13C ********************************************************************** 14C Reads from the input file the data to calculate STS spectra 15C 16C Coded by P. Ordejon, July 2004 17C ********************************************************************** 18 19 USE FDF 20 21 IMPLICIT NONE 22 23 DOUBLE PRECISION, INTENT(IN) :: 24 . VOLUME 25 26 INTEGER, INTENT(OUT) :: 27 . NE, IUNITCD, BFUNC 28 29 DOUBLE PRECISION, INTENT(OUT) :: 30 . EMIN, DELTAE, ETA, XSTS(3), ARMUNI 31 32 33C ****** INPUT ******************************************************** 34C REAL*8 VOLUME : Volume of the unit cell 35C ****** OUTPUT ******************************************************** 36C REAL*8 EMIN : Minimum energy in STS spectrum 37C INTEGER NE : Number of energy points 38C REAL*8 DELTAE : Step in energy plot 39C REAL*8 ETA : Lorenzian broadening of levels 40C REAL*8 XSTS(3) : Position where LDOS is computed 41C INTEGER BFUNC : Broadening function for STS spectra 42C 1 = Lorentzian, 0 = Gaussian 43C INTEGER IUNITCD : Units for the electron density 44C IUNITCD = 1 => Ele/(bohr)**3 45C IUNITCD = 2 => Ele/(Ang)**3 46C IUNITCD = 3 => Ele/(unitcell) 47C REAL*8 ARMUNI : Conversion factors for the charge density 48C ********************************************************************** 49 50C Internal variables --------------------------------------------------- 51C INTEGER ISCALE : Units for the atomic positions 52C (ISCALE = 1 => Bohrs, ISCALE = 2 => Ang) 53 CHARACTER 54 . CPF*22, CPF_DEFECT*22, 55 . UCD*22, UCD_DEFECT*22, 56 . BROAD*22, BROAD_DEFECT*22 57 58 INTEGER 59 . IUNIT, IX, JX, ISCALE 60 61 DOUBLE PRECISION EMAX 62 63 LOGICAL 64 . LEQI 65 66 EXTERNAL 67 . LEQI 68C ---- 69 70 CPF_DEFECT = 'Bohr' 71 72 CPF = FDF_STRING('Denchar.CoorUnits',CPF_DEFECT) 73 74 IF (LEQI(CPF,'bohr')) then 75 ISCALE = 1 76 ELSEIF (LEQI(CPF,'ang')) then 77 ISCALE = 2 78 ELSE 79 WRITE(6,*)'readsts: ERROR Denchar.CoorUnits must be Ang or Bohr' 80 STOP 81 ENDIF 82 83 UCD_DEFECT = 'Ele/bohr**3' 84 UCD = FDF_STRING('Denchar.DensityUnits',UCD_DEFECT) 85 IF (LEQI(UCD,'ele/bohr**3')) then 86 IUNITCD = 1 87 ELSEIF (LEQI(UCD,'ele/ang**3')) then 88 IUNITCD = 2 89 ELSEIF (LEQI(UCD,'ele/unitcell')) then 90 IUNITCD = 3 91 ELSE 92 WRITE(6,'(A)')' readsts: ERROR Wrong Option in Units of ' 93 WRITE(6,'(A)')' readsts: Charge Density ' 94 WRITE(6,'(A)')' readsts: You must choose one of the following:' 95 WRITE(6,'(A)')' readsts: ' 96 WRITE(6,'(A)')' readsts: - Ele/bohr**3 ' 97 WRITE(6,'(A)')' readsts: - Ele/ang**3 ' 98 WRITE(6,'(A)')' readsts: - Ele/unitcell ' 99 STOP 100 ENDIF 101 102 103 BROAD_DEFECT = 'Lorentzian' 104 105 BROAD = FDF_STRING('Denchar.STSBroadeningFunc',BROAD_DEFECT) 106 107 IF (LEQI(BROAD,'lorentzian')) then 108 BFUNC = 1 109 ELSEIF (LEQI(BROAD,'gaussian')) then 110 BFUNC = 2 111 ELSE 112 WRITE(6,*)'readsts: ERROR Denchar.STSBroadeningFunc must be' 113 WRITE(6,*)'readsts: Lorentzian or Gaussian' 114 STOP 115 ENDIF 116 117 NE = FDF_INTEGER('Denchar.STSEnergyPoints',100) 118 119 IF ( FDF_BLOCK('Denchar.STSPosition',IUNIT) ) THEN 120 READ(IUNIT,*)(XSTS(IX),IX=1,3) 121 ELSE 122 WRITE(6,*)'To calculate STS spectra, you must' 123 WRITE(6,*)'specify a position (Denchar.STSPosition)' 124 WRITE(6,*)'in the input file' 125 STOP 126 ENDIF 127 128 EMIN = FDF_PHYSICAL('Denchar.STSEmin',-1.0D0,'eV') 129 EMAX = FDF_PHYSICAL('Denchar.STSEmax',1.0D0,'eV') 130 ETA = FDF_PHYSICAL('Denchar.STSEta',0.1D0,'eV') 131 132 DELTAE = (EMAX-EMIN)/DFLOAT(NE-1) 133 134C Scale points coordinates 135C Iscale = 1 => Do nothing 136C Iscale = 2 => Multiply by 1./0.529177 (Ang --> Bohr) 137 138 IF (ISCALE .EQ. 2) THEN 139 DO IX = 1,3 140 XSTS(IX) = 1.D0 / 0.529177D0 * XSTS(IX) 141 ENDDO 142 ENDIF 143 144C Units of Charge Density 145C Iunitcd = 1 => Do nothing 146C Iunitcd = 2 => Multiply by (1.d0 / 0.529177d0) **3 (bohr**3 --> Ang**3) 147C Iunitcd = 3 => Multiply by volume unit cell (in bohrs**3) 148 149 IF (IUNITCD .EQ. 1) THEN 150 ARMUNI = 1.D0 151 ELSEIF( IUNITCD .EQ. 2 ) THEN 152 ARMUNI = (1.D0 / 0.529177D0)**3 153 ELSEIF( IUNITCD .EQ. 3 ) THEN 154 ARMUNI = VOLUME 155 ENDIF 156 157 END 158 159