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