pot_slat.f90

   1
   2     SUBROUTINE POTENTIAL_WZ(EEL,EHHUP,EHHDW,ELHUP,ELHDW,&
   3                             ESOUP,ESODW,ELAST,EXX,EYY,&
   4                             EZZ,EXY,EXZ,EYZ,&
   5                             POT)
   6
   7               Use Input_Data
   8               Use Dot_Geometry
   9               Use Auxiliar_Procedures, ONLY : AISO
  10
  11       IMPLICIT NONE
  12
  13 !!!!! 'dummy' and local variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  14
  15       REAL,DIMENSION(:,:,:),OPTIONAL ::  EXX,EYY,EZZ,EXY,EXZ,EYZ,POT
  16
  17       REAL,DIMENSION(:,:,:) ::  EEL,EHHUP,EHHDW,ELHUP,ELHDW,&
  18                                 ESOUP,ESODW,ELAST
  19
  20       REAL         ZM,THETA,CTHETA,STHETA, &
  21                    X,Y,Z,ZMAUX,RHO,ZETA
  22
  23       INTEGER      I_X,I_Y,I_Z,I_N1,I_N2,I_N3,CHI
  24
  25       REAL, DIMENSION(3) :: R_SL,X_VEC,XI_VEC
  26
  27       REAL :: POTBE,POTBHH,POTBLH,POTBSO,POTBLS,&
  28               POTWE,POTWHH,POTWLH,POTWSO,POTWLS,&
  29               VBIEL,VBIHH,VBILH,VBISO,VBILS
  30
  31       REAL,DIMENSION(0:1) :: POTE, POTHH, POTLH, POTSO, POTLS, &
  32                              PC1, PC2, PD1, PD2, PD3, PD4, PD5, PD6
  33
  34       REAL, DIMENSION(1:8) :: ENERGY
  35
  36 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  37
  38       CALL ENER_CONSTANTS_WZ( )
  39
  40 ZD:   DO I_Z=1,ZDim
  41 !       WRITE(16,'(A,I3,A,I3)')"I_Z ",I_Z," of ",ZDIM
  42         Z=Z_Min+REAL(I_Z-1)*Z_Inc
  43 YD:   DO I_Y=1,YDim
  44         Y=Y_Min+REAL(I_Y-1)*Y_Inc
  45 XD:   DO I_X=1,XDim
  46         X=X_Min+REAL(I_X-1)*X_Inc
  47
  48         X_VEC=(/X,Y,Z/)
  49
  50         ENERGY=0.E0
  51
  52       I_N1=0; I_N2=0; I_N3=0
  53
  54       R_SL=REAL(I_N1)*A1_S+REAL(I_N2)*A2_S+REAL(I_N3)*A3_S
  55
  56       XI_VEC=X_VEC-R_SL
  57
  58       RHO=SQRT(XI_VEC(1)**2+XI_VEC(2)**2)/RC
  59       IF(XI_VEC(1).EQ.0.E0.AND.XI_VEC(2).EQ.0.E0) THEN
  60               CTHETA=1.E0/SQRT(2.E0); STHETA=1.E0/SQRT(2.E0) ! It is not the Mathematical limit
  61       ELSE
  62               THETA=ATAN(XI_VEC(2)/XI_VEC(1))
  63               CTHETA=Cos(THETA); STHETA=Sin(THETA)
  64       END IF
  65       ZETA=XI_VEC(3)/ZC
  66
  67           IF (RHO.LE.RD) THEN
  68              CALL SHAPERTOZ(MIN(RHO*RC,Rqd_Base),ZMAUX)
  69              ZM=ZMAUX/ZC
  70           ELSE
  71              ZM = 0.E0
  72           END IF
  73
  74           IF (abs(zeta) .EQ. 0.E0 .OR. ZETA .EQ. ZM) THEN
  75               ZETA=ZETA-1.E-5
  76           END IF
  77
  78           CHI = 0
  79           IF (ZETA.GE.-D.AND.ZETA.LE.ZM) THEN
  80                   CHI = 1
  81                   IF(I_N1.NE.0.OR.I_N2.NE.0.OR.I_N3.NE.0) THEN
  82                           WRITE(16,*)I_N1,I_N2,I_N3
  83                           WRITE(16,*)X_VEC(3),ZETA*ZC,ZM*ZC
  84                   END IF
  85           END IF
  86
  87           IF (KCOOR.EQ.0) THEN
  88               CALL POT_CALC_CAR_WZ(CHI,ENERGY)
  89           ELSE
  90               CALL POT_CALC_CYL_WZ(CHI,ENERGY)
  91           END IF
  92
  93            EEL(I_X,I_Y,I_Z)   = ENERGY(1)
  94            EHHUP(I_X,I_Y,I_Z) = ENERGY(2)
  95            EHHDW(I_X,I_Y,I_Z) = ENERGY(3)
  96            ELHUP(I_X,I_Y,I_Z) = ENERGY(4)
  97            ELHDW(I_X,I_Y,I_Z) = ENERGY(5)
  98            ESOUP(I_X,I_Y,I_Z) = ENERGY(6)
  99            ESODW(I_X,I_Y,I_Z) = ENERGY(7)
 100            ELAST(I_X,I_Y,I_Z) = ENERGY(8)
 101
 102 !           WRITE(26,'(10(E15.8,1X))')Z,ENERGY(1:8)
 103
 104       END DO XD
 105       END DO YD
 106       END DO ZD
 107
 108 !     STOP
 109
 110       RETURN
 111
 112       CONTAINS
 113
 114       SUBROUTINE ENER_CONSTANTS_WZ( )
 115       IMPLICIT NONE
 116
 117 ! Definition of parameters for Barrier and Well
 118
 119          VWE=VWE+DW1+DW2
 120          VBE=VBE+DB1+DB2
 121
 122          POTWE  =  VWE
 123          POTWHH = (VWH+DW1+DW2)
 124          POTWLH = (VWH+(DW1-DW2+4.E0*DW3)/3.E0)
 125          POTWSO = (VWH+2.E0*(DW1-DW2-2.E0*DW3)/3.E0)
 126          POTWLS = (DW1-DW2+DW3)
 127          POTBE  =  VBE
 128          POTBHH = (VBH+DB1+DB2)
 129          POTBLH = (VBH+(DB1-DB2+4.E0*DB3)/3.E0)
 130          POTBSO = (VBH+2.E0*(DB1-DB2-2.E0*DB3)/3.E0)
 131          POTBLS = (DB1-DB2+DB3)
 132
 133          IF(STR_Action.EQ.0) THEN
 134           VBIEL=           ( C2*BISUM + C1*BIZZ )
 135           VBIHH=      ( (D2+D4)*BISUM + (D1+D3)*BIZZ )
 136           VBILH=   ( (D2+D4/3.)*BISUM + (D1+D3/3.)*BIZZ )
 137           VBISO=( (D2+2.*D4/3.)*BISUM + (D1+2.*D3/3.)*BIZZ )
 138           VBILS=           ( D4*BISUM + D3*BIZZ )
 139
 140 !         Potential edges including strain effect
 141           POTWE= (POTWE+VBIEL)
 142           POTWHH= (POTWHH+VBIHH)
 143           POTWLH= (POTWLH+VBILH)
 144           POTWSO= (POTWSO+VBISO)
 145           POTWLS= (POTWLS+VBILS)
 146          END IF
 147
 148          POTE(0)  = POTBE  ; POTE(1)  = POTWE
 149          POTHH(0) = POTBHH ; POTHH(1) = POTWHH
 150          POTLH(0) = POTBLH ; POTLH(1) = POTWLH
 151          POTLS(0) = POTBLS ; POTLS(1) = POTWLS
 152          POTSO(0) = POTBSO ; POTSO(1) = POTWSO
 153
 154 !!!! Deformation potentials for barrier and well are equal
 155
 156          PC1(0) = C1 ; PC1(1) = C1
 157          PC2(0) = C2 ; PC2(1) = C2
 158          PD1(0) = D1 ; PD1(1) = D1
 159          PD2(0) = D2 ; PD2(1) = D2
 160          PD3(0) = D3 ; PD3(1) = D3
 161          PD4(0) = D4 ; PD4(1) = D4
 162          PD5(0) = D5 ; PD5(1) = D5
 163          PD6(0) = D6 ; PD6(1) = D6
 164
 165          RETURN
 166
 167       END SUBROUTINE ENER_CONSTANTS_WZ
 168
 169       SUBROUTINE POT_CALC_CAR_WZ(CHI,ENERGY)
 170       IMPLICIT NONE
 171
 172         REAL, DIMENSION(:) :: ENERGY
 173         INTEGER, PARAMETER :: DIMM=6
 174         COMPLEX, DIMENSION(1:DIMM,1:DIMM) :: HKANE
 175         REAL :: SXX,SYY,SZZ,SXY,SXZ,SYZ,PZO
 176         INTEGER :: I,J,CHI
 177 !! DIAGONALIZATION VARIABLES
 178
 179         COMPLEX, DIMENSION(1:(2*DIMM-1)) :: WORK
 180         REAL, DIMENSION(1:(3*DIMM-2)) :: RWORK
 181         REAL         W(1:DIMM), ROOT(DIMM), CARAC(DIMM)
 182         INTEGER      INFO, LWORK
 183
 184         LWORK=2*DIMM-1
 185
 186         IF(STR_Action.EQ.0) THEN
 187           SXX=0.E0; SYY=0.E0; SZZ=0.E0
 188           SXY=0.E0; SXZ=0.E0; SYZ=0.E0
 189         ELSE
 190           SXX=EXX(I_X,I_Y,I_Z)
 191           SYY=EYY(I_X,I_Y,I_Z)
 192           SZZ=EZZ(I_X,I_Y,I_Z)
 193           SXY=EXY(I_X,I_Y,I_Z)
 194           SXZ=EXZ(I_X,I_Y,I_Z)
 195           SYZ=EYZ(I_X,I_Y,I_Z)
 196         END IF
 197         IF(PZO_Action.EQ.0.E0) THEN
 198           PZO=0.E0
 199         ELSE
 200           PZO=POT(I_X,I_Y,I_Z)
 201         END IF
 202
 203          HKANE(1,1)=POTHH(CHI)-PZO+ &
 204                      (PD1(CHI)+PD3(CHI))*SZZ+(PD2(CHI)+PD4(CHI))*(SXX+SYY)
 205
 206          HKANE(2,2)=POTLH(CHI)-PZO+(PD1(CHI)+PD3(CHI)/3.E0)*SZZ+ &
 207                      (PD2(CHI)+PD4(CHI)/3.E0)*(SXX+SYY)
 208          HKANE(3,3)=HKANE(2,2)
 209          HKANE(4,4)=HKANE(1,1)
 210          HKANE(5,5)=POTSO(CHI)-PZO+(PD1(CHI)+2.E0*PD3(CHI)/3.E0)*SZZ+ &
 211                      (PD2(CHI)+2.E0*PD4(CHI)/3.E0)*(SXX+SYY)
 212          HKANE(6,6)=HKANE(5,5)
 213
 214          IF(DIAG_Action.NE.0.AND.STR_Action.NE.0) THEN
 215
 216           HKANE(1,2)=-SQRT(2.E0/3.E0)*PD6(CHI)*CMPLX(SXZ,-SYZ)
 217           HKANE(1,3)=-1.E0/SQRT(3.E0)*PD5(CHI)*CMPLX(SXX-SYY,-2.E0*SXY)
 218           HKANE(1,4)=0.E0
 219           HKANE(1,5)=1.E0/SQRT(3.E0)*PD6(CHI)*CMPLX(SXZ,-SYZ)
 220           HKANE(1,6)=SQRT(2./3.)*PD5(CHI)*CMPLX(SXX-SYY,-2.E0*SXY)
 221
 222           HKANE(2,3)=0.
 223           HKANE(2,4)=HKANE(1,3)
 224           HKANE(2,5)=SQRT(2.E0)/3.E0*( POTLS(CHI)+&
 225                                        (PD3(CHI)*SZZ+PD4(CHI)*(SXX+SYY)) )
 226           HKANE(2,6)=-PD6(CHI)*CMPLX(SXZ,-SYZ)
 227
 228           HKANE(3,4)=-HKANE(1,2)
 229           HKANE(3,5)=-PD6(CHI)*CMPLX(SXZ,SYZ)
 230           HKANE(3,6)=-HKANE(2,5)
 231
 232           HKANE(4,5)=-SQRT(2./3.)*PD5(CHI)*CMPLX(SXX-SYY,+2.E0*SXY)
 233           HKANE(4,6)=PD6(CHI)*CMPLX(SXZ,SYZ)/SQRT(3.E0)
 234
 235           HKANE(5,6)=0
 236
 237
 238
 239           DO I=1,DIMM
 240            DO J=I+1,DIMM
 241             HKANE(J,I)=CONJG(HKANE(I,J))
 242            END DO
 243           END DO
 244
 245           CALL CHEEV('V','U',DIMM,HKANE,DIMM,W,WORK,LWORK,RWORK,INFO)
 246
 247           IF (INFO.ne.0) then
 248                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 249                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
 250 !              STOP
 251           END IF
 252
 253 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 254 !!   Ordering Eigenvalues according to the Bloch func. caracter
 255 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 256
 257           ROOT = 0.E0
 258
 259           DO I=1,6,2
 260 ! Calculo de la componentes dependientes del spin
 261            DO J=1,6
 262             CARAC(J)=HKANE(J,I)*HKANE(J,I)
 263            END DO
 264 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 265            CARAC(1)=CARAC(1)+CARAC(4)
 266            CARAC(4)=CARAC(1)
 267            CARAC(2)=CARAC(2)+CARAC(3)
 268            CARAC(3)=CARAC(2)
 269            CARAC(5)=CARAC(5)+CARAC(6)
 270            CARAC(6)=CARAC(5)
 271 !
 272 ! Peso de las componentes de las que se extraeran los autovalores.
 273 !
 274            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 275                    ROOT(4)=W(I)
 276            END IF
 277            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 278                    ROOT(3)=W(I)
 279            END IF
 280            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
 281                    ROOT(1)=W(I)
 282            END IF
 283
 284           END DO
 285
 286          ELSE ! Only the diagonal elements were calculated
 287
 288           ROOT(4)=HKANE(1,1)
 289           ROOT(3)=HKANE(2,2)
 290           ROOT(1)=HKANE(5,5)
 291           ROOT(2)=HKANE(6,6) !Redundant
 292
 293
 294          END IF
 295
 296 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 297 !! RESULTS
 298 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 299
 300          ENERGY(1) = POTE(CHI)-PZO+(PC1(CHI)*SZZ+PC2(CHI)*(SXX+SYY))
 301
 302          ENERGY(2) = ROOT(4)
 303          ENERGY(3) = ROOT(3)
 304          ENERGY(4) = ROOT(3)
 305          ENERGY(5) = ROOT(4)
 306
 307          ENERGY(6) = ROOT(1)
 308          ENERGY(7) = ROOT(1)
 309 !        ENERGY(7) = CHI ! To save the Structure-profile
 310
 311 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 312 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 313 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
 314 !!!!! we will use this array to pack the Elastic Energy.
 315 !!!!!
 316 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
 317 !!!!!
 318 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 319 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 320
 321           ENERGY(8) = C11*(SXX**2 + SYY**2)+C33*SZZ**2+ &
 322                       2.E0*(XLAMB*(SXX*SYY)+C13*(SXX+SYY)*SZZ+ &
 323                       (C11-XLAMB)*SXY**2+2.E0*XMU*(SXZ**2+SYZ**2))
 324
 325           ENERGY(8) = ENERGY(8)/2.E0
 326
 327           RETURN
 328
 329       END SUBROUTINE POT_CALC_CAR_WZ
 330
 331       SUBROUTINE POT_CALC_CYL_WZ(CHI,ENERGY)
 332       IMPLICIT NONE
 333
 334         REAL, DIMENSION(:) :: ENERGY
 335         INTEGER, PARAMETER :: DIMM=6
 336         REAL, DIMENSION(1:DIMM,1:DIMM) :: HKANE
 337         REAL :: PZO
 338         REAL :: SSUM,SDIF,SZZ,SRZ,SRR,S00
 339         INTEGER :: I,J,CHI
 340 !! DIAGONALIZATION VARIABLES
 341
 342         REAL ::    W(1:DIMM), WORK(1:10*DIMM), AW(1:DIMM,1:DIMM), &
 343                    ROOT(DIMM), CARAC(DIMM), SLAMCH
 344         INTEGER :: IWORK(DIMM*5), IFAIL(DIMM), INFO, NUM
 345
 346         IF(STR_Action.EQ.0) THEN
 347           SSUM=0.E0; SDIF=0.E0; SZZ=0.E0;
 348           SRZ=0.E0; SRR=0.E0; S00=0.E0
 349         ELSE
 350
 351 !!! FOR CYLINDRICAL COORDINATES WE CALCULATE ONLY THE X-Z PLANE,
 352 !!! THAT MEANS: THETA=0
 353
 354           SRR=EXX(I_X,I_Y,I_Z)
 355           S00=EYY(I_X,I_Y,I_Z)
 356           SSUM=EXX(I_X,I_Y,I_Z)+EYY(I_X,I_Y,I_Z)
 357           SDIF=EXX(I_X,I_Y,I_Z)-EYY(I_X,I_Y,I_Z)
 358           SZZ=EZZ(I_X,I_Y,I_Z)
 359           SRZ=EXZ(I_X,I_Y,I_Z)
 360         END IF
 361         IF(PZO_Action.EQ.0.E0) THEN
 362           PZO=0.E0
 363         ELSE
 364           PZO=POT(I_X,I_Y,I_Z)
 365         END IF
 366
 367          HKANE(1,1)=POTHH(CHI)-PZO+ &
 368                      (PD1(CHI)+PD3(CHI))*SZZ+(PD2(CHI)+PD4(CHI))*SSUM
 369
 370          HKANE(2,2)=POTLH(CHI)-PZO+(PD1(CHI)+PD3(CHI)/3.E0)*SZZ+ &
 371                      (PD2(CHI)+PD4(CHI)/3.E0)*SSUM
 372          HKANE(3,3)=HKANE(2,2)
 373          HKANE(4,4)=HKANE(1,1)
 374          HKANE(5,5)=POTSO(CHI)-PZO+(PD1(CHI)+2.E0*PD3(CHI)/3.E0)*SZZ+ &
 375                      (PD2(CHI)+2.E0*PD4(CHI)/3.E0)*SSUM
 376          HKANE(6,6)=HKANE(5,5)
 377
 378          IF(DIAG_Action.NE.0.AND.STR_Action.NE.0) THEN
 379
 380           HKANE(1,2)=-SQRT(2.E0/3.E0)*PD6(CHI)*SRZ
 381           HKANE(1,3)=-1.E0/SQRT(3.E0)*PD5(CHI)*SDIF
 382           HKANE(1,4)=0.
 383           HKANE(1,5)=1.E0/SQRT(3.E0)*PD6(CHI)*SRZ
 384           HKANE(1,6)=SQRT(2./3.)*PD5(CHI)*SDIF
 385
 386           HKANE(2,3)=0.
 387           HKANE(2,4)=HKANE(1,3)
 388           HKANE(2,5)=SQRT(2.E0)/3.E0*(POTLS(CHI)+(PD3(CHI)*SZZ+PD4(CHI)*SSUM))
 389           HKANE(2,6)=-PD6(CHI)*SRZ
 390
 391           HKANE(3,4)=-HKANE(1,2)
 392           HKANE(3,5)=HKANE(2,6)
 393           HKANE(3,6)=-HKANE(2,5)
 394
 395           HKANE(4,5)=-HKANE(1,6)
 396           HKANE(4,6)=HKANE(1,5)
 397
 398           HKANE(5,6)=0
 399
 400
 401
 402           DO I=1,6
 403            DO J=I+1,6
 404             HKANE(J,I)=HKANE(I,J)
 405            END DO
 406           END DO
 407
 408 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 409 !!   Analytic solutions for Rho = 0.
 410 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 411 !
 412 !         ZETA = ZMIN+ZINC*IZ ! XZ normalized to ZC
 413 !
 414 !        W(1)=(HKANE(2,2)+HKANE(5,5))+
 415 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 416 !        W(2)=(HKANE(2,2)+HKANE(5,5))-
 417 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 418 !        write(24,'(4(f18.8,1x))')ZETA,HKANE(1,1),W(1)/2.,W(2)/2.
 419 !
 420 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 421
 422
 423           CALL DSYEVX('V','A','U',DIMM,HKANE,DIMM,0.,0.,2,2, &
 424            2*SLAMCH('S'),NUM,W,AW,DIMM,WORK,10*DIMM,IWORK,IFAIL,INFO)
 425
 426           IF (INFO.ne.0) then
 427                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 428                write(6,'(A,10(I2,1X))')"IFAIL=",IFAIL(1:DIMM)
 429                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
 430 !              STOP
 431           END IF
 432
 433 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 434 !!   Ordering Eigenvalues according to the Bloch func. caracter
 435 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 436
 437           ROOT = 0.E0
 438
 439           DO I=1,6,2
 440 ! Calculo de la componentes dependientes del spin
 441            DO J=1,6
 442             CARAC(J)=AW(J,I)*AW(J,I)
 443            END DO
 444 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 445            CARAC(1)=CARAC(1)+CARAC(4)
 446            CARAC(4)=CARAC(1)
 447            CARAC(2)=CARAC(2)+CARAC(3)
 448            CARAC(3)=CARAC(2)
 449            CARAC(5)=CARAC(5)+CARAC(6)
 450            CARAC(6)=CARAC(5)
 451 !
 452 ! Peso de las componentes de las que se extraeran los autovalores.
 453 !
 454            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 455                    ROOT(4)=W(I)
 456            END IF
 457            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 458                    ROOT(3)=W(I)
 459            END IF
 460            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
 461                    ROOT(1)=W(I)
 462            END IF
 463
 464           END DO
 465
 466          ELSE ! Only the diagonal elements were calculated
 467
 468           ROOT(4)=HKANE(1,1)
 469           ROOT(3)=HKANE(2,2)
 470           ROOT(1)=HKANE(5,5)
 471           ROOT(2)=HKANE(6,6) !Redundant
 472
 473
 474          END IF
 475
 476 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 477 !! RESULTS
 478 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 479
 480          ENERGY(1) = POTE(CHI)-PZO+(PC1(CHI)*SZZ+PC2(CHI)*SSUM)
 481
 482          ENERGY(2) = ROOT(4)
 483          ENERGY(3) = ROOT(3)
 484          ENERGY(4) = ROOT(3)
 485          ENERGY(5) = ROOT(4)
 486
 487          ENERGY(6) = ROOT(1)
 488          ENERGY(7) = ROOT(1)
 489 !        ESODW(IR,IZ,IGEO) = CHI ! To save the Structure-profile
 490
 491 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 492 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 493 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
 494 !!!!! we will use this array to pack the Elastic Energy.
 495 !!!!!
 496 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
 497 !!!!!
 498 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 499 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 500
 501
 502           ENERGY(8) = C11*(SRR**2 + S00**2)+C33*SZZ**2+ &
 503                       2.E0*(XLAMB*(SRR*S00)+C13*(SRR+S00)*SZZ+ &
 504                       +4.E0*XMU*(SRZ**2))
 505
 506           ENERGY(8) = ENERGY(8)/2.E0
 507
 508           RETURN
 509
 510       END SUBROUTINE POT_CALC_CYL_WZ
 511
 512     END SUBROUTINE POTENTIAL_WZ
 513
 514     SUBROUTINE POTENTIAL_ZB(EEL,EHHUP,EHHDW,ELHUP,ELHDW,&
 515                             ESOUP,ESODW,ELAST,EXX,EYY,&
 516                             EZZ,EXY,EXZ,EYZ)
 517
 518
 519               Use Input_Data
 520               Use Dot_Geometry
 521               Use Auxiliar_Procedures, ONLY : AISO
 522
 523       IMPLICIT NONE
 524
 525 !!!!! 'dummy' and local variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 526
 527       REAL,DIMENSION(:,:,:),OPTIONAL ::  EXX,EYY,EZZ,EXY,EXZ,EYZ
 528
 529       REAL,DIMENSION(:,:,:) ::  EEL,EHHUP,EHHDW,ELHUP,ELHDW,&
 530                                 ESOUP,ESODW,ELAST
 531
 532       REAL         ZM,THETA,CTHETA,STHETA, &
 533                    X,Y,Z,ZMAUX,RHO,ZETA
 534
 535       INTEGER      I_X,I_Y,I_Z,I_N1,I_N2,I_N3,CHI
 536
 537       REAL, DIMENSION(3) :: R_SL,X_VEC,XI_VEC
 538
 539       REAL :: POTWE,POTWHH,POTWLH,POTWSO,&
 540               VBIEL,VBIHH,VBILH,VBISO
 541
 542       REAL,DIMENSION(0:1) :: POTE, POTHH, POTLH, POTSO, &
 543                              DVD,DSD,DVU,DVSU,D2VU,D2VSU,DSO,ZBC1
 544
 545       REAL, DIMENSION(1:8) :: ENERGY
 546
 547 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 548
 549       CALL ENER_CONSTANTS_ZB( )
 550
 551 ZD:   DO I_Z=1,ZDim
 552 !       WRITE(16,'(A,I3,A,I3)')"I_Z ",I_Z," of ",ZDIM
 553         Z=Z_Min+REAL(I_Z-1)*Z_Inc
 554 YD:   DO I_Y=1,YDim
 555         Y=Y_Min+REAL(I_Y-1)*Y_Inc
 556 XD:   DO I_X=1,XDim
 557         X=X_Min+REAL(I_X-1)*X_Inc
 558
 559         X_VEC=(/X,Y,Z/)
 560
 561       I_N1=0; I_N2=0; I_N3=0
 562
 563       R_SL=REAL(I_N1)*A1_S+REAL(I_N2)*A2_S+REAL(I_N3)*A3_S
 564
 565       XI_VEC=X_VEC-R_SL
 566
 567       RHO=SQRT(XI_VEC(1)**2+XI_VEC(2)**2)/RC
 568       IF(XI_VEC(1).EQ.0.E0.AND.XI_VEC(2).EQ.0.E0) THEN
 569               CTHETA=1.E0/SQRT(2.E0); STHETA=1.E0/SQRT(2.E0) ! It is not the Mathematical limit
 570       ELSE
 571               THETA=ATAN(XI_VEC(2)/XI_VEC(1))
 572               CTHETA=Cos(THETA); STHETA=Sin(THETA)
 573       END IF
 574       ZETA=XI_VEC(3)/ZC
 575
 576           IF (RHO.LE.RD) THEN
 577              CALL SHAPERTOZ(MIN(RHO*RC,Rqd_Base),ZMAUX)
 578              ZM=ZMAUX/ZC
 579           ELSE
 580              ZM = 0.E0
 581           END IF
 582
 583           IF (abs(zeta) .EQ. 0.E0 .OR. ZETA .EQ. ZM) THEN
 584               ZETA=ZETA-1.E-5
 585           END IF
 586
 587           CHI = 0
 588           IF (RHO.LE.RD.AND.ZETA.GE.-D.AND.ZETA.LE.ZM) THEN
 589                   CHI = 1
 590                   IF(I_N1.NE.0.OR.I_N2.NE.0.OR.I_N3.NE.0) THEN
 591                           WRITE(16,*)I_N1,I_N2,I_N3
 592                           WRITE(16,*)X_VEC(3),ZETA*ZC,ZM*ZC
 593                   END IF
 594           END IF
 595
 596           IF (KCOOR.EQ.0) THEN
 597               CALL POT_CALC_CAR_ZB(CHI,ENERGY)
 598           ELSE
 599               CALL POT_CALC_CYL_ZB(CHI,ENERGY)
 600           END IF
 601
 602            EEL(I_X,I_Y,I_Z)   = ENERGY(1)
 603            EHHUP(I_X,I_Y,I_Z) = ENERGY(2)
 604            ELHUP(I_X,I_Y,I_Z) = ENERGY(3)
 605            ELHDW(I_X,I_Y,I_Z) = ENERGY(4)
 606            EHHDW(I_X,I_Y,I_Z) = ENERGY(5)
 607            ESOUP(I_X,I_Y,I_Z) = ENERGY(6)
 608            ESODW(I_X,I_Y,I_Z) = ENERGY(7)
 609            ELAST(I_X,I_Y,I_Z) = ENERGY(8)
 610
 611 !           WRITE(26,'(10(E15.8,1X))')Z,ENERGY(1:8)
 612
 613       END DO XD
 614       END DO YD
 615       END DO ZD
 616
 617 !     STOP
 618
 619       RETURN
 620
 621       CONTAINS
 622
 623       SUBROUTINE ENER_CONSTANTS_ZB( )
 624       IMPLICIT NONE
 625
 626 ! Definition of parameters for Barrier and Well
 627
 628 ! If we set both terms equal to zero the spin-interaction in the
 629 ! deformation potentials is removed. The values of AVB,BB,DB are
 630 ! taken from C. Pryor, PRB, 57, 7190 (1998)
 631
 632 !        DSO(0) = VBH - VBSO ; DSO(1) = VWH - VWSO
 633          DSO(0) = 0.0; DSO(1) = 0.0
 634
 635          ZBC1(0) = -9.3 ; ZBC1(1) = ACW
 636
 637          DVD(0) = 0.7E0 - 2./9. * DSO(0)
 638          DVD(1) = AVW - 2./9. * DSO(1)
 639          DSD(0) = 0.7E0 + 4./9. * DSO(0)
 640          DSD(1) = AVW + 4./9. * DSO(1)
 641          DVU(0) = -3./2. * (-2.0) + 1./3. * DSO(0)
 642          DVU(1) = -3./2. * BW + 1./3. * DSO(1)
 643          DVSU(0) = -3./2. * (-2.0) - 1./6. * DSO(0)
 644          DVSU(1) = -3./2. * BW - 1./6. * DSO(1)
 645          D2VU(0) = -SQRT(3.)/2. * (-5.4) + 1./3. * DSO(0)
 646          D2VU(1) = -SQRT(3.)/2. * DW + 1./3. * DSO(1)
 647          D2VSU(0) = -SQRT(3.)/2. * (-5.4) - 1./6. * DSO(0)
 648          D2VSU(1) = -SQRT(3.)/2. * DW - 1./6. * DSO(1)
 649
 650 !!! In the calculation the deformation potentials are equal across the structure
 651          DVD(0)=DVD(1); DSD(0)=DSD(1); DVU(0)=DVU(1)
 652          DVSU(0)=DVSU(1); D2VU(0)=D2VU(1); D2VSU(0)=D2VSU(1)
 653
 654          DSO(0) = VBH - VBSO ; DSO(1) = VWH - VWSO
 655          POTWE= VWE
 656          POTWHH=VWH
 657          POTWLH=VWH
 658          POTWSO=VWSO
 659
 660           IF(STR_Action.EQ.0) THEN
 661            VBIEL  = ACW*(BISUM+BIZZ)
 662            VBIHH  = AVW*(BISUM+BIZZ)-(-3.E0*BW/2.E0)/3.E0*(BISUM-2.E0*BIZZ)
 663            VBILH  = AVW*(BISUM+BIZZ)+(-3.E0*BW/2.E0)/3.E0*(BISUM-2.E0*BIZZ)
 664            VBISO  = AVW*(BISUM+BIZZ)
 665
 666 !          Potential edges including strain effect
 667            POTWE= (VWE+VBIEL)
 668            POTWHH= (VWH+VBIHH)
 669            POTWLH= (VWH+VBILH)
 670            POTWSO= (VWSO+VBISO)
 671           END IF
 672
 673          POTE(0)  = VBE  ; POTE(1)  = POTWE
 674          POTHH(0) = VBH ;  POTHH(1) = POTWHH
 675          POTLH(0) = VBH ;  POTLH(1) = POTWLH
 676          POTSO(0) = VBSO ; POTSO(1) = POTWSO
 677
 678          RETURN
 679
 680       END SUBROUTINE ENER_CONSTANTS_ZB
 681
 682       SUBROUTINE POT_CALC_CAR_ZB(CHI,ENERGY)
 683       IMPLICIT NONE
 684
 685         REAL, DIMENSION(:) :: ENERGY
 686         INTEGER, PARAMETER :: DIMM=6
 687         COMPLEX, DIMENSION(1:DIMM,1:DIMM) :: HKANE
 688         REAL :: SXX,SYY,SZZ,SXY,SXZ,SYZ,SHID,SDIF,STIL
 689         INTEGER :: I,J,CHI
 690 !! DIAGONALIZATION VARIABLES
 691
 692         COMPLEX, DIMENSION(1:(2*DIMM-1)) :: WORK
 693         REAL, DIMENSION(1:(3*DIMM-2)) :: RWORK
 694         REAL         W(1:DIMM), ROOT(DIMM), CARAC(DIMM)
 695         INTEGER      INFO, LWORK
 696
 697
 698         IF(STR_Action.EQ.0) THEN
 699           SXX=0.E0; SYY=0.E0; SZZ=0.E0
 700           SXY=0.E0; SXZ=0.E0; SYZ=0.E0
 701         ELSE
 702           SXX=EXX(I_X,I_Y,I_Z)
 703           SYY=EYY(I_X,I_Y,I_Z)
 704           SZZ=EZZ(I_X,I_Y,I_Z)
 705           SXY=EXY(I_X,I_Y,I_Z)
 706           SXZ=EXZ(I_X,I_Y,I_Z)
 707           SYZ=EYZ(I_X,I_Y,I_Z)
 708           SHID=SXX+SYY+SZZ
 709           SDIF=SXX-SYY
 710           STIL=SXX+SYY-2.E0*SZZ
 711         END IF
 712
 713          HKANE(1,1)=POTHH(CHI)+DVD(CHI)*SHID &
 714                      -1./3.*DVU(CHI)*STIL
 715          HKANE(2,2)=POTLH(CHI)+DVD(CHI)*SHID &
 716                      +1./3.*DVU(CHI)*STIL
 717          HKANE(3,3)=HKANE(2,2)
 718          HKANE(4,4)=HKANE(1,1)
 719          HKANE(5,5)=POTSO(CHI)+DSD(CHI)*SHID
 720
 721          HKANE(6,6)=HKANE(5,5)
 722
 723          IF(DIAG_Action.NE.0.AND.STR_Action.NE.0) THEN
 724
 725            HKANE(1,2)=2./SQRT(3.)*D2VU(CHI)*CMPLX(SXZ,-SYZ)
 726            HKANE(1,3)=1./SQRT(3.)*CMPLX(DVU(CHI)*SDIF,-2.E0*D2VU(CHI)*SXY)
 727            HKANE(1,4)=0.E0
 728            HKANE(1,5)=-SQRT(2./3.)*D2VSU(CHI)*CMPLX(SXZ,-SYZ)
 729            HKANE(1,6)=-SQRT(2./3.)*CMPLX(DVU(CHI)*SDIF,-2.E0*D2VSU(CHI)*SXY)
 730
 731            HKANE(2,3)=0.
 732            HKANE(2,4)=HKANE(1,3)
 733            HKANE(2,5)=-SQRT(2.)/3.*DVSU(CHI)*STIL
 734            HKANE(2,6)=SQRT(2.)*D2VSU(CHI)*CMPLX(SXZ,-SYZ)
 735
 736            HKANE(3,4)=-HKANE(1,2)
 737            HKANE(3,5)=SQRT(2.)*D2VSU(CHI)*CMPLX(SXZ,SYZ)
 738            HKANE(3,6)=-HKANE(2,5)
 739
 740            HKANE(4,5)=SQRT(2./3.)*CMPLX(DVU(CHI)*SDIF,2.E0*D2VSU(CHI)*SXY)
 741            HKANE(4,6)=-HKANE(3,5)/SQRT(3.E0)
 742
 743            HKANE(5,6)=0
 744
 745          DO I=1,6
 746           DO J=I+1,6
 747            HKANE(J,I)=CONJG(HKANE(I,J))
 748           END DO
 749          END DO
 750
 751 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 752 !!   Analytic solutions for Rho = 0.
 753 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 754 !
 755 !         ZETA = ZMIN+ZINC*IZ ! XZ normalized to ZC
 756 !
 757 !        W(1)=(HKANE(2,2)+HKANE(5,5))+
 758 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 759 !        W(2)=(HKANE(2,2)+HKANE(5,5))-
 760 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 761 !        write(24,'(4(f18.8,1x))')ZETA,HKANE(1,1),W(1)/2.,W(2)/2.
 762 !
 763 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 764
 765           CALL CHEEV('V','U',DIMM,HKANE,DIMM,W,WORK,LWORK,RWORK,INFO)
 766
 767           IF (INFO.ne.0) then
 768                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 769                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
 770 !              STOP
 771           END IF
 772
 773 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 774 !!   Ordering Eigenvalues according to the Bloch func. caracter
 775 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 776
 777           ROOT = 0.E0
 778
 779           DO I=1,6,2
 780 ! Calculo de la componentes dependientes del spin
 781            DO J=1,6
 782             CARAC(J)=HKANE(J,I)*HKANE(J,I)
 783            END DO
 784 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 785            CARAC(1)=CARAC(1)+CARAC(4)
 786            CARAC(4)=CARAC(1)
 787            CARAC(2)=CARAC(2)+CARAC(3)
 788            CARAC(3)=CARAC(2)
 789            CARAC(5)=CARAC(5)+CARAC(6)
 790            CARAC(6)=CARAC(5)
 791 !
 792 ! Peso de las componentes de las que se extraeran los autovalores.
 793 !
 794            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 795                    ROOT(4)=W(I)
 796            END IF
 797            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 798                    ROOT(3)=W(I)
 799            END IF
 800            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
 801                    ROOT(1)=W(I)
 802            END IF
 803
 804           END DO
 805
 806          ELSE ! Only the diagonal elements were calculated
 807
 808           ROOT(4)=HKANE(1,1)
 809           ROOT(3)=HKANE(2,2)
 810           ROOT(1)=HKANE(5,5)
 811           ROOT(2)=HKANE(6,6) !Redundant
 812
 813
 814          END IF
 815
 816 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 817 !! RESULTS
 818 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 819
 820          ENERGY(1) = POTE(CHI)+ZBC1(CHI)*SHID
 821
 822          ENERGY(2) = ROOT(4)
 823          ENERGY(3) = ROOT(3)
 824          ENERGY(4) = ROOT(3)
 825          ENERGY(5) = ROOT(4)
 826
 827          ENERGY(6) = ROOT(1)
 828          ENERGY(7) = ROOT(1)
 829 !        ENERGY(7) = CHI ! To save the Structure-profile
 830
 831 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 832 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 833 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
 834 !!!!! we will use this array to pack the Elastic Energy.
 835 !!!!!
 836 !!!!!        U=1/2*(XLAMB*Tr(e)**2+XMU/2*(err**2+e00**2+ezz**2+erz**2)
 837 !!!!!
 838 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 839 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 840
 841           ENERGY(8) = XMU*(SXX**2 + SYY**2 + SZZ**2 + &
 842                            2.E0*(SXY**2+SXZ**2+SYZ**2) ) + &
 843                       XLAMB/2.E0*(SXX + SYY + SZZ)**2
 844           RETURN
 845
 846       END SUBROUTINE POT_CALC_CAR_ZB
 847
 848       SUBROUTINE POT_CALC_CYL_ZB(CHI,ENERGY)
 849       IMPLICIT NONE
 850
 851         REAL, DIMENSION(:) :: ENERGY
 852         INTEGER, PARAMETER :: DIMM=6
 853         REAL, DIMENSION(1:DIMM,1:DIMM) :: HKANE
 854         REAL :: SSUM,SDIF,SZZ,SRZ,SRR,S00,SHID,STIL
 855         INTEGER :: I,J,CHI
 856 !! DIAGONALIZATION VARIABLES
 857
 858         REAL ::    W(1:DIMM), WORK(1:10*DIMM), AW(1:DIMM,1:DIMM), &
 859                    ROOT(DIMM), CARAC(DIMM), SLAMCH
 860         INTEGER :: IWORK(DIMM*5), IFAIL(DIMM), INFO, NUM
 861
 862
 863         IF(STR_Action.EQ.0) THEN
 864           SSUM=0.E0; SDIF=0.E0; SZZ=0.E0;
 865           SRZ=0.E0; SRR=0.E0; S00=0.E0
 866         ELSE
 867           SRR=EXX(I_X,I_Y,I_Z)
 868           S00=EYY(I_X,I_Y,I_Z)
 869           SSUM=EXX(I_X,I_Y,I_Z)+EYY(I_X,I_Y,I_Z)
 870           SDIF=EXX(I_X,I_Y,I_Z)-EYY(I_X,I_Y,I_Z)
 871           SZZ=EZZ(I_X,I_Y,I_Z)
 872           SRZ=EXZ(I_X,I_Y,I_Z)
 873           SHID=SSUM+SZZ
 874           STIL=SSUM-2.E0*SZZ
 875         END IF
 876
 877          HKANE(1,1)=POTHH(CHI)+DVD(CHI)*SHID &
 878                     -1./3.*DVU(CHI)*STIL
 879          HKANE(2,2)=POTLH(CHI)+DVD(CHI)*SHID &
 880                     +1./3.*DVU(CHI)*STIL
 881
 882          HKANE(3,3)=HKANE(2,2)
 883          HKANE(4,4)=HKANE(1,1)
 884          HKANE(5,5)=POTSO(CHI)+DSD(CHI)*SHID
 885          HKANE(6,6)=HKANE(5,5)
 886
 887          IF(DIAG_Action.NE.0.AND.STR_Action.NE.0) THEN
 888             HKANE(1,2)=2./SQRT(3.)*D2VU(CHI)*SRZ
 889             HKANE(1,3)=1./SQRT(3.)*(DVU(CHI)+D2VU(CHI))/2.*SDIF
 890             HKANE(1,4)=0.E0
 891             HKANE(1,5)=-SQRT(2./3.)*D2VSU(CHI)*SRZ
 892             HKANE(1,6)=-SQRT(2./3.)*(DVU(CHI)+D2VSU(CHI))/2.*SDIF
 893
 894             HKANE(2,3)=0.E0
 895             HKANE(2,4)=HKANE(1,3)
 896             HKANE(2,5)=-SQRT(2.)/3.*DVSU(CHI)*STIL
 897             HKANE(2,6)=SQRT(2.)*D2VSU(CHI)*SRZ
 898
 899             HKANE(3,4)=-HKANE(1,2)
 900             HKANE(3,5)=HKANE(2,6)
 901             HKANE(3,6)=-HKANE(2,5)
 902
 903             HKANE(4,5)=-HKANE(1,6)
 904             HKANE(4,6)=HKANE(1,5)
 905
 906             HKANE(5,6)=0
 907
 908           DO I=1,6
 909            DO J=I+1,6
 910             HKANE(J,I)=HKANE(I,J)
 911            END DO
 912           END DO
 913
 914 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 915 !!   Analytic solutions for Rho = 0.
 916 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 917 !
 918 !         ZETA = ZMIN+ZINC*IZ ! XZ normalized to ZC
 919 !
 920 !        W(1)=(HKANE(2,2)+HKANE(5,5))+
 921 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 922 !        W(2)=(HKANE(2,2)+HKANE(5,5))-
 923 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 924 !        write(24,'(4(f18.8,1x))')ZETA,HKANE(1,1),W(1)/2.,W(2)/2.
 925 !
 926 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 927
 928
 929           CALL DSYEVX('V','A','U',DIMM,HKANE,DIMM,0.,0.,2,2, &
 930            2*SLAMCH('S'),NUM,W,AW,DIMM,WORK,10*DIMM,IWORK,IFAIL,INFO)
 931
 932           IF (INFO.ne.0) then
 933                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 934                write(6,'(A,10(I2,1X))')"IFAIL=",IFAIL(1:DIMM)
 935                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
 936 !              STOP
 937           END IF
 938
 939 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 940 !!   Ordering Eigenvalues according to the Bloch func. caracter
 941 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 942
 943           ROOT = 0.E0
 944
 945           DO I=1,6,2
 946 ! Calculo de la componentes dependientes del spin
 947            DO J=1,6
 948             CARAC(J)=AW(J,I)*AW(J,I)
 949            END DO
 950 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 951            CARAC(1)=CARAC(1)+CARAC(4)
 952            CARAC(4)=CARAC(1)
 953            CARAC(2)=CARAC(2)+CARAC(3)
 954            CARAC(3)=CARAC(2)
 955            CARAC(5)=CARAC(5)+CARAC(6)
 956            CARAC(6)=CARAC(5)
 957 !
 958 ! Peso de las componentes de las que se extraeran los autovalores.
 959 !
 960            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 961                    ROOT(4)=W(I)
 962            END IF
 963            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 964                    ROOT(3)=W(I)
 965            END IF
 966            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
 967                    ROOT(1)=W(I)
 968            END IF
 969
 970           END DO
 971
 972          ELSE ! Only the diagonal elements were calculated
 973
 974           ROOT(4)=HKANE(1,1)
 975           ROOT(3)=HKANE(2,2)
 976           ROOT(1)=HKANE(5,5)
 977           ROOT(2)=HKANE(6,6) !Redundant
 978
 979
 980          END IF
 981
 982 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 983 !! RESULTS
 984 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 985
 986          ENERGY(1) = POTE(CHI)+ZBC1(CHI)*SHID
 987
 988          ENERGY(2) = ROOT(4)
 989          ENERGY(3) = ROOT(3)
 990          ENERGY(4) = ROOT(3)
 991          ENERGY(5) = ROOT(4)
 992
 993          ENERGY(6) = ROOT(1)
 994          ENERGY(7) = ROOT(1)
 995 !        ESODW(IR,IZ,IGEO) = CHI ! To save the Structure-profile
 996
 997 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 998 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 999 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
1000 !!!!! we will use this array to pack the Elastic Energy.
1001 !!!!!
1002 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
1003 !!!!!
1004 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1005 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1006
1007           ENERGY(8) = XMU*(SRR**2 + S00**2 + SZZ**2 + SRZ**2) +    &
1008                       XLAMB/2.E0*(SSUM + SZZ)**2
1009
1010           RETURN
1011
1012       END SUBROUTINE POT_CALC_CYL_ZB
1013
1014       END SUBROUTINE POTENTIAL_ZB
1015