potdrv.f

   1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   2 ! Program to draw the potential energy from a file containig the strain
   3 ! profiles
   4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   5
   6
   7
   8
   9       PROGRAM POT_PROFILE
  10
  11            Use Dot_Geometry
  12
  13       IMPLICIT NONE
  14
  15       INTEGER,PARAMETER :: DIMM=6
  16
  17 ***** 'dummy' and local variables **************************************
  18
  19       REAL            VWE,VWH,DW1,DW2,DW3,VBE,VBH,DB1,DB2,DB3
  20
  21       INTEGER         Number_Param,IGEO,NWF
  22       REAL            VBIEL,VBIHH,VBILH,VBILS,VBISO
  23       REAL            POTWE,POTWHH,POTWLH,POTWLS,POTWSO,
  24      &                POTBE,POTBHH,POTBLH,POTBLS,POTBSO
  25       REAL            C1,C2,D1,D2,D3,D4,D5,D6,BISUM,BIZZ,BIAUX
  26       REAL            XMU,XLAMB,EPS0
  27       REAL            RC,ZC,DWL,RHO,ZETA,ZM
  28
  29       INTEGER I,IIR,J,K,IZ,IR
  30
  31 !! GRID RESULTS OF STRAIN
  32
  33       REAL,ALLOCATABLE,DIMENSION(:,:,:) :: MSUM, MDIF, MRR,
  34      &                                     M00,  MZZ,  MRZ,
  35      &                                     MHID, MTIL
  36
  37 !! GRID RESULTS OF PIEZOELECTRIC POTENTIAL
  38
  39       REAL,ALLOCATABLE,DIMENSION(:,:,:) :: PTOT
  40
  41 !! GRID RESULTS OF POTENTIALS
  42
  43       REAL,ALLOCATABLE,DIMENSION(:,:,:) :: EEL,   EHHUP, EHHDW,
  44      &                                     ELHUP, ELHDW, ESOUP,
  45      &                                     ESODW, ELAST
  46
  47 !! FRACTIONAL COORDINATES AND GEOM. PARAMETERS
  48
  49       REAL,ALLOCATABLE :: XVALS(:), YVALS(:), RADI(:,:), ZETI(:,:), PAR(:)
  50
  51 !! HAMILTONIAN
  52
  53       REAL         HKANE(1:DIMM,1:DIMM)
  54       INTEGER      CHI
  55
  56 !! DIAGONALIZATION VARIABLES
  57
  58       REAL         W(1:DIMM), WORK(1:10*DIMM), AW(1:DIMM,1:DIMM),
  59      &             ROOT(DIMM), CARAC(DIMM), SLAMCH
  60       INTEGER      IWORK(DIMM*5), IFAIL(DIMM), INFO, NUM
  61
  62
  63 !! DEFROMATION POTENTIALS AND ENERGY PROFILES
  64       REAL,DIMENSION(0:1) :: POTE, POTHH, POTLH, POTSO, POTLS,
  65      &                       PC1, PC2, PD1, PD2, PD3, PD4, PD5, PD6
  66
  67       REAL         ERR,E00,EZZ,ERZ,ESUM,EDIF,PZO
  68
  69 !!!!!!! COMMON BLOCKS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  70
  71       CHARACTER(LEN=120) :: input_str,input_pzo,output_eng
  72       COMMON /FILES/        input_str,input_pzo,output_eng
  73
  74       REAL                  CARAC_MAX
  75       COMMON /CARAC/        CARAC_MAX
  76
  77 !! MORE INOUT.INC VARIABLES
  78
  79       LOGICAL            :: DIAGLOG,PZOLOG,STRLOG,DWLLOG
  80       COMMON /LOGICS/       DIAGLOG,PZOLOG,STRLOG
  81
  82 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  83
  84       INCLUDE    'INC/inout.inc'
  85
  86 ! Allocating Arrays
  87
  88       ALLOCATE(MSUM(0:NWF,0:NWF,1:Number_Param),
  89      &         MDIF(0:NWF,0:NWF,1:Number_Param),
  90      &         MRR(0:NWF,0:NWF,1:Number_Param),
  91      &         M00(0:NWF,0:NWF,1:Number_Param),
  92      &         MZZ(0:NWF,0:NWF,1:Number_Param),
  93      &         MRZ(0:NWF,0:NWF,1:Number_Param),
  94      &         MHID(0:NWF,0:NWF,1:Number_Param),
  95      &         MTIL(0:NWF,0:NWF,1:Number_Param),
  96      &         PTOT(0:NWF,0:NWF,1:Number_Param))
  97
  98       ALLOCATE(EEL(0:NWF,0:NWF,1:Number_Param),
  99      &         EHHUP(0:NWF,0:NWF,1:Number_Param),
 100      &         EHHDW(0:NWF,0:NWF,1:Number_Param),
 101      &         ELHUP(0:NWF,0:NWF,1:Number_Param),
 102      &         ELHDW(0:NWF,0:NWF,1:Number_Param),
 103      &         ESOUP(0:NWF,0:NWF,1:Number_Param),
 104      &         ELAST(0:NWF,0:NWF,1:Number_Param),
 105      &         ESODW(0:NWF,0:NWF,1:Number_Param))
 106
 107       ALLOCATE(XVALS(0:NWF),YVALS(0:NWF),RADI(1:Number_Param,1:3),
 108      &         ZETI(1:Number_Param,1:2),PAR(1:Number_Param))
 109
 110 ! Definition of parameters for Barrier and Well
 111
 112       VWE=VWE+DW1+DW2
 113       VBE=VBE+DB1+DB2
 114
 115       POTWE  =  VWE
 116       POTWHH = (VWH+DW1+DW2)
 117       POTWLH = (VWH+(DW1-DW2+4.E0*DW3)/3.E0)
 118       POTWSO = (VWH+2.E0*(DW1-DW2-2.E0*DW3)/3.E0)
 119       POTWLS = (DW1-DW2+DW3)
 120       POTBE  =  VBE
 121       POTBHH = (VBH+DB1+DB2)
 122       POTBLH = (VBH+(DB1-DB2+4.E0*DB3)/3.E0)
 123       POTBSO = (VBH+2.E0*(DB1-DB2-2.E0*DB3)/3.E0)
 124       POTBLS = (DB1-DB2+DB3)
 125
 126       IF(.NOT.STRLOG) THEN
 127        VBIEL=           ( C2*BISUM + C1*BIZZ )
 128        VBIHH=      ( (D2+D4)*BISUM + (D1+D3)*BIZZ )
 129        VBILH=   ( (D2+D4/3.)*BISUM + (D1+D3/3.)*BIZZ )
 130        VBISO=( (D2+2.*D4/3.)*BISUM + (D1+2.*D3/3.)*BIZZ )
 131        VBILS=           ( D4*BISUM + D3*BIZZ )
 132
 133 !      Potential edges including strain effect
 134        POTWE= (POTWE+VBIEL)
 135        POTWHH= (POTWHH+VBIHH)
 136        POTWLH= (POTWLH+VBILH)
 137        POTWSO= (POTWSO+VBISO)
 138        POTWLS= (POTWLS+VBILS)
 139       END IF
 140
 141       POTE(0)  = POTBE  ; POTE(1)  = POTWE
 142       POTHH(0) = POTBHH ; POTHH(1) = POTWHH
 143       POTLH(0) = POTBLH ; POTLH(1) = POTWLH
 144       POTLS(0) = POTBLS ; POTLS(1) = POTWLS
 145       POTSO(0) = POTBSO ; POTSO(1) = POTWSO
 146
 147 ! Deformation potentials for barrier and well are equal
 148
 149       PC1(0) = C1 ; PC1(1) = C1
 150       PC2(0) = C2 ; PC2(1) = C2
 151       PD1(0) = D1 ; PD1(1) = D1
 152       PD2(0) = D2 ; PD2(1) = D2
 153       PD3(0) = D3 ; PD3(1) = D3
 154       PD4(0) = D4 ; PD4(1) = D4
 155       PD5(0) = D5 ; PD5(1) = D5
 156       PD6(0) = D6 ; PD6(1) = D6
 157
 158       IF(STRLOG) THEN
 159
 160        CALL STR_UNPACK(MSUM,MDIF,MRR,M00,MZZ,MRZ,MHID,MTIL,Number_Param,NWF,
 161      &               RADI,ZETI,PAR,XVALS,YVALS)
 162       ELSE
 163
 164         MSUM=0.E0; MDIF=0.E0; MRR=0.E0; M00=0.E0; MZZ=0.E0; MRZ=0.E0
 165         MHID=0.E0; MTIL=0.E0
 166
 167       END IF
 168
 169       IF(PZOLOG) THEN
 170
 171        CALL PZO_UNPACK(Number_Param,NWF,PTOT,RADI,ZETI,PAR,XVALS,YVALS)
 172
 173       ELSE
 174
 175        PTOT=0.E0
 176
 177       END IF
 178
 179       DO IGEO=1,Number_Param
 180
 181         IIR=0
 182
 183         Rqd_Base = RADI(IGEO,1)
 184         Rqd_Top  = RADI(IGEO,2)
 185         Hqd      = ZETI(IGEO,1)
 186         ZC       = ZETI(IGEO,2)
 187         RC       = RADI(IGEO,3)
 188         write(12,*)Rqd_Base,Rqd_Top,Hqd,ZC,RC
 189
 190 !! In case that the parameter would be DWL:
 191         IF(DWLLOG) DWL = PAR(IGEO)
 192
 193         DO IR=NWF/2,NWF
 194
 195          RHO     = XVALS(IR)*RC
 196
 197          IF (RHO.LE.Rqd_Base) THEN
 198             CALL SHAPERTOZ(RHO,ZM)
 199             write(12,*)RHO,ZM
 200          ELSE
 201             ZM = 0.E0
 202             write(12,*)RHO,ZM
 203          END IF
 204
 205
 206         DO IZ=0,NWF
 207
 208          ZETA   = YVALS(IZ)*ZC
 209
 210          ERR    = MRR(IR,IZ,Number_Param)
 211          E00    = M00(IR,IZ,Number_Param)
 212          EZZ    = MZZ(IR,IZ,Number_Param)
 213          ERZ    = MRZ(IR,IZ,Number_Param)
 214          ESUM   = MSUM(IR,IZ,Number_Param)
 215          EDIF   = MDIF(IR,IZ,Number_Param)
 216          PZO    = PTOT(IR,IZ,Number_Param)
 217
 218          CHI=0
 219          IF (ZETA.GE.-DWL .AND. ZETA.LE.ZM) CHI=1
 220
 221          HKANE(1,1)=POTHH(CHI)-PZO+
 222      &               (PD1(CHI)+PD3(CHI))*EZZ+(PD2(CHI)+PD4(CHI))*ESUM
 223
 224          HKANE(2,2)=POTLH(CHI)-PZO+(PD1(CHI)+PD3(CHI)/3.E0)*EZZ+
 225      &               (PD2(CHI)+PD4(CHI)/3.E0)*ESUM
 226          HKANE(3,3)=HKANE(2,2)
 227          HKANE(4,4)=HKANE(1,1)
 228          HKANE(5,5)=POTSO(CHI)-PZO+(PD1(CHI)+2.E0*PD3(CHI)/3.E0)*EZZ+
 229      &               (PD2(CHI)+2.E0*PD4(CHI)/3.E0)*ESUM
 230          HKANE(6,6)=HKANE(5,5)
 231
 232          IF(.NOT.DIAGLOG.AND.STRLOG) THEN
 233
 234           HKANE(1,2)=-SQRT(2.E0/3.E0)*PD6(CHI)*ERZ
 235           HKANE(1,3)=-1.E0/SQRT(3.E0)*PD5(CHI)*EDIF
 236           HKANE(1,4)=0.
 237           HKANE(1,5)=1.E0/SQRT(3.E0)*PD6(CHI)*ERZ
 238           HKANE(1,6)=SQRT(2./3.)*PD5(CHI)*EDIF
 239
 240           HKANE(2,3)=0.
 241           HKANE(2,4)=HKANE(1,3)
 242           HKANE(2,5)=SQRT(2.E0)/3.E0*(POTLS(CHI)+(PD3(CHI)*EZZ+PD4(CHI)*ESUM))
 243           HKANE(2,6)=-PD6(CHI)*ERZ
 244
 245           HKANE(3,4)=-HKANE(1,2)
 246           HKANE(3,5)=HKANE(2,6)
 247           HKANE(3,6)=-HKANE(2,5)
 248
 249           HKANE(4,5)=-HKANE(1,6)
 250           HKANE(4,6)=HKANE(1,5)
 251
 252           HKANE(5,6)=0
 253
 254
 255
 256           DO I=1,6
 257            DO J=I+1,6
 258             HKANE(J,I)=HKANE(I,J)
 259            END DO
 260           END DO
 261
 262 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 263 !!   Analytic solutions for Rho = 0.
 264 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 265 !
 266 !         ZETA = ZMIN+ZINC*IZ ! XZ normalized to ZC
 267 !
 268 !        W(1)=(HKANE(2,2)+HKANE(5,5))+
 269 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 270 !        W(2)=(HKANE(2,2)+HKANE(5,5))-
 271 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 272 !        write(24,'(4(f18.8,1x))')ZETA,HKANE(1,1),W(1)/2.,W(2)/2.
 273 !
 274 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 275
 276
 277           CALL SSYEVX('V','A','U',DIMM,HKANE,DIMM,0.,0.,2,2,
 278      &     2*SLAMCH('S'),NUM,W,AW,DIMM,WORK,10*DIMM,IWORK,IFAIL,INFO)
 279
 280           IF (INFO.ne.0) then
 281                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 282                write(6,'(A,10(I2,1X))')"IFAIL=",IFAIL(1:DIMM)
 283                write(6,*)"IR=",IR," IZ=",IZ
 284 !              STOP
 285           END IF
 286
 287 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 288 !!   Ordering Eigenvalues according to the Bloch func. caracter
 289 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 290
 291           ROOT = 0.E0
 292
 293           DO I=1,6,2
 294 ! Calculo de la componentes dependientes del spin
 295            DO J=1,6
 296             CARAC(J)=AW(J,I)*AW(J,I)
 297            END DO
 298 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 299            CARAC(1)=CARAC(1)+CARAC(4)
 300            CARAC(4)=CARAC(1)
 301            CARAC(2)=CARAC(2)+CARAC(3)
 302            CARAC(3)=CARAC(2)
 303            CARAC(5)=CARAC(5)+CARAC(6)
 304            CARAC(6)=CARAC(5)
 305 !
 306 ! Peso de las componentes de las que se extraeran los autovalores.
 307 !
 308            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 309                    ROOT(4)=W(I)
 310            END IF
 311            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 312                    ROOT(3)=W(I)
 313            END IF
 314            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
 315                    ROOT(1)=W(I)
 316            END IF
 317
 318           END DO
 319
 320          ELSE ! Only the diagonal elements were calculated
 321
 322           ROOT(4)=HKANE(1,1)
 323           ROOT(3)=HKANE(2,2)
 324           ROOT(1)=HKANE(5,5)
 325           ROOT(2)=HKANE(6,6) !Redundant
 326
 327
 328          END IF
 329
 330 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 331 !! RESULTS
 332 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 333
 334          EEL(IR,IZ,IGEO)   = POTE(CHI)-PZO+(PC1(CHI)*EZZ+PC2(CHI)*ESUM)
 335
 336          EHHUP(IR,IZ,IGEO) = ROOT(4)
 337          ELHUP(IR,IZ,IGEO) = ROOT(3)
 338          ELHDW(IR,IZ,IGEO) = ROOT(3)
 339          EHHDW(IR,IZ,IGEO) = ROOT(4)
 340
 341          ESOUP(IR,IZ,IGEO) = ROOT(1)
 342          ESODW(IR,IZ,IGEO) = ROOT(1)
 343 !        ESODW(IR,IZ,IGEO) = CHI ! To save the Structure-profile
 344
 345 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 346 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 347 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
 348 !!!!! we will use this array to pack the Elastic Energy.
 349 !!!!!
 350 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
 351 !!!!!
 352 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 353 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 354
 355           ELAST(IR,IZ,IGEO)=XMU*(ERR**2 + E00**2 + EZZ**2 + ERZ**2) +
 356      &                      0.5*XLAMB*(ERR + E00 + EZZ)**2
 357
 358          IIR=1
 359         END DO
 360
 361        END DO
 362
 363       END DO
 364
 365
 366        DO IGEO=1,Number_Param
 367         DO IZ=0,NWF
 368          DO IR=0,NWF/2-1
 369           EEL(ir,iz,IGEO)=EEL(NWF-IR,iz,IGEO)
 370           EHHUP(ir,iz,IGEO)=EHHUP(NWF-IR,iz,IGEO)
 371           EHHDW(ir,iz,IGEO)=EHHDW(NWF-IR,iz,IGEO)
 372           ELHUP(ir,iz,IGEO)=ELHUP(NWF-IR,iz,IGEO)
 373           ELHDW(ir,iz,IGEO)=ELHDW(NWF-IR,iz,IGEO)
 374           ESOUP(ir,iz,IGEO)=ESOUP(NWF-IR,iz,IGEO)
 375           ESODW(ir,iz,IGEO)=ESODW(NWF-IR,iz,IGEO)
 376           ELAST(ir,iz,IGEO)=ELAST(NWF-IR,iz,IGEO)
 377          END DO
 378         END DO
 379        END DO
 380
 381       CALL POT_PACK(EEL,EHHUP,ELHUP,ELHDW,EHHDW,
 382      &           ESOUP,ESODW,ELAST,Number_Param,NWF,RADI,ZETI,PAR,XVals,YVals)
 383
 384       END PROGRAM POT_PROFILE