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 (KPBASIS .EQ. 0) THEN ! Winkler basis
  88             IF (KCOOR.EQ.0) THEN
  89                 CALL POT_CALC_CAR_WZ(CHI,ENERGY)
  90             ELSE
  91                 CALL POT_CALC_CYL_WZ(CHI,ENERGY)
  92             END IF
  93           ELSE
  94             IF (KCOOR.EQ.0) THEN
  95                 CALL POT_CALC_CAR_CH_WZ(CHI,ENERGY)
  96             ELSE
  97                 CALL POT_CALC_CYL_CH_WZ(CHI,ENERGY)
  98             END IF
  99           END IF
 100
 101            EEL(I_X,I_Y,I_Z)   = ENERGY(1)
 102            EHHUP(I_X,I_Y,I_Z) = ENERGY(2)
 103            ELHUP(I_X,I_Y,I_Z) = ENERGY(3)
 104            ELHDW(I_X,I_Y,I_Z) = ENERGY(4)
 105            EHHDW(I_X,I_Y,I_Z) = ENERGY(5)
 106            ESOUP(I_X,I_Y,I_Z) = ENERGY(6)
 107            ESODW(I_X,I_Y,I_Z) = ENERGY(7)
 108            ELAST(I_X,I_Y,I_Z) = ENERGY(8)
 109
 110 !           WRITE(26,'(10(E15.8,1X))')Z,ENERGY(1:8)
 111
 112       END DO XD
 113       END DO YD
 114       END DO ZD
 115
 116 !     STOP
 117
 118       RETURN
 119
 120       CONTAINS
 121
 122       SUBROUTINE ENER_CONSTANTS_WZ( )
 123       IMPLICIT NONE
 124
 125       REAL, DIMENSION(3) :: EW,EB
 126       REAL :: EVW,EVB,ECW,ECB
 127       INTEGER, DIMENSION(1) :: AUX
 128       INTEGER :: NW,NB
 129
 130          EVW=0.E0; EVB=0.E0
 131
 132          EW(1)=EVW+DW1+DW2
 133          EW(2)=EVW+(DW1-DW2)/2.E0+SQRT(((DW1-DW2)/2.E0)**2+2.E0*DW3**2)
 134          EW(3)=EVW+(DW1-DW2)/2.E0-SQRT(((DW1-DW2)/2.E0)**2+2.E0*DW3**2)
 135
 136          EB(1)=EVB+DB1+DB2
 137          EB(2)=EVB+(DB1-DB2)/2.E0+SQRT(((DB1-DB2)/2.E0)**2+2.E0*DB3**2)
 138          EB(3)=EVB+(DB1-DB2)/2.E0-SQRT(((DB1-DB2)/2.E0)**2+2.E0*DB3**2)
 139
 140          AUX=MAXLOC(EW); NW=AUX(1)
 141          AUX=MAXLOC(EB); NB=AUX(1)
 142
 143          EVB=0.E0
 144          EVW=VBO+EB(NB)-EW(NW)
 145
 146          EW(1)=EVW+DW1+DW2
 147          EW(2)=EVW+(DW1-DW2)/2.E0+SQRT(((DW1-DW2)/2.E0)**2+2.E0*DW3**2)
 148          EW(3)=EVW+(DW1-DW2)/2.E0-SQRT(((DW1-DW2)/2.E0)**2+2.E0*DW3**2)
 149
 150          EB(1)=EVB+DB1+DB2
 151          EB(2)=EVB+(DB1-DB2)/2.E0+SQRT(((DB1-DB2)/2.E0)**2+2.E0*DB3**2)
 152          EB(3)=EVB+(DB1-DB2)/2.E0-SQRT(((DB1-DB2)/2.E0)**2+2.E0*DB3**2)
 153
 154          ECW=EGW+MAXVAL(EW)
 155          ECB=EGB+MAXVAL(EB)
 156
 157          VWE=ECW; VBE=ECB
 158          VWH=EVW; VBH=EVB
 159
 160 ! Definition of parameters for Barrier and Well
 161
 162
 163          POTWE  =  VWE
 164          POTWHH = (VWH+DW1+DW2)
 165          POTBE  =  VBE
 166          POTBHH = (VBH+DB1+DB2)
 167
 168          IF (KPBASIS == 0) THEN
 169            POTWLH = (VWH+(DW1-DW2+4.E0*DW3)/3.E0)
 170            POTWSO = (VWH+2.E0*(DW1-DW2-2.E0*DW3)/3.E0)
 171            POTWLS = (DW1-DW2+DW3)
 172            POTBLH = (VBH+(DB1-DB2+4.E0*DB3)/3.E0)
 173            POTBSO = (VBH+2.E0*(DB1-DB2-2.E0*DB3)/3.E0)
 174            POTBLS = (DB1-DB2+DB3)
 175          ELSE
 176            POTWLH = (VWH+DW1-DW2)
 177            POTWLS = DW3
 178            POTWSO = VWH
 179            POTBLH = (VBH+DB1-DB2)
 180            POTBLS = DB3
 181            POTBSO = VBH
 182          END IF
 183
 184 !!!!! Output Potential Profiles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 185
 186       WRITE(6,*)"PERFILES DE POTENCIAL SIN DEFORMACION:"
 187       WRITE(6,*)"POTWE= ",POTWE," POTBE= ",POTBE
 188       WRITE(6,*)"POTWHH= ",POTWHH," POTBHH= ",POTBHH
 189       WRITE(6,*)"POTWLH= ",POTWLH," POTBLH= ",POTBLH
 190       WRITE(6,*)"POTWSO= ",POTWSO," POTBSO= ",POTBSO
 191       WRITE(6,*)"POTWLS= ",POTWLS," POTBLS= ",POTBLS
 192 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 193
 194
 195 !!!!! Output Potential Profiles including strain effect !!!!!!!!!!!!!!!!!!!!!!!!
 196
 197          IF(STR_Action.EQ.0) THEN
 198             IF (KPBASIS == 0) THEN
 199               VBIEL=           ( C2*BISUM + C1*BIZZ )
 200               VBIHH=      ( (D2+D4)*BISUM + (D1+D3)*BIZZ )
 201               VBILH=   ( (D2+D4/3.)*BISUM + (D1+D3/3.)*BIZZ )
 202               VBISO=( (D2+2.*D4/3.)*BISUM + (D1+2.*D3/3.)*BIZZ )
 203               VBILS=           ( D4*BISUM + D3*BIZZ )
 204             ELSE
 205               VBIEL=        ( C2*BISUM + C1*BIZZ )
 206               VBIHH=   ( (D2+D4)*BISUM + (D1+D3)*BIZZ )
 207               VBILH=   ( (D2+D4)*BISUM + (D1+D3)*BIZZ )
 208               VBISO=         (D2*BISUM + D1*BIZZ )
 209               VBILS= 0.E0
 210             END IF
 211          END IF
 212
 213
 214
 215 !             Potential edges including strain effect
 216          IF(STR_Action.EQ.0) THEN
 217           POTWE= (POTWE+VBIEL)
 218           POTWHH= (POTWHH+VBIHH)
 219           POTWLH= (POTWLH+VBILH)
 220           POTWSO= (POTWSO+VBISO)
 221           POTWLS= (POTWLS+VBILS)
 222           WRITE(6,*)"PERFILES DE POTENCIAL CON DEFORMACION:"
 223           WRITE(6,*)"POTWE= ",POTWE," POTBE= ",POTBE
 224           WRITE(6,*)"POTWHH= ",POTWHH," POTBHH= ",POTBHH
 225           WRITE(6,*)"POTWLH= ",POTWLH," POTBLH= ",POTBLH
 226           WRITE(6,*)"POTWSO= ",POTWSO," POTBSO= ",POTBSO
 227           WRITE(6,*)"POTWLS= ",POTWLS," POTBLS= ",POTBLS
 228          END IF
 229
 230              POTE(0)  = POTBE  ; POTE(1)  = POTWE
 231              POTHH(0) = POTBHH ; POTHH(1) = POTWHH
 232              POTLH(0) = POTBLH ; POTLH(1) = POTWLH
 233              POTLS(0) = POTBLS ; POTLS(1) = POTWLS
 234              POTSO(0) = POTBSO ; POTSO(1) = POTWSO
 235
 236 !!!! Deformation potentials for barrier and well are equal
 237
 238          PC1(0) = C1 ; PC1(1) = C1
 239          PC2(0) = C2 ; PC2(1) = C2
 240          PD1(0) = D1 ; PD1(1) = D1
 241          PD2(0) = D2 ; PD2(1) = D2
 242          PD3(0) = D3 ; PD3(1) = D3
 243          PD4(0) = D4 ; PD4(1) = D4
 244          PD5(0) = D5 ; PD5(1) = D5
 245          PD6(0) = D6 ; PD6(1) = D6
 246
 247          RETURN
 248
 249       END SUBROUTINE ENER_CONSTANTS_WZ
 250
 251       SUBROUTINE POT_CALC_CAR_CH_WZ(CHI,ENERGY)
 252       IMPLICIT NONE
 253
 254         REAL, DIMENSION(:) :: ENERGY
 255         INTEGER, PARAMETER :: DIMM=6
 256         COMPLEX, DIMENSION(1:DIMM,1:DIMM) :: HKANE
 257         REAL :: SXX,SYY,SZZ,SXY,SXZ,SYZ,PZO
 258         INTEGER :: I,J,CHI
 259 !! DIAGONALIZATION VARIABLES
 260
 261         COMPLEX, DIMENSION(1:(2*DIMM-1)) :: WORK
 262         REAL, DIMENSION(1:(3*DIMM-2)) :: RWORK
 263         REAL         W(1:DIMM), ROOT(DIMM), CARAC(DIMM)
 264         INTEGER      INFO, LWORK
 265
 266         LWORK=2*DIMM-1
 267
 268         IF(STR_Action.EQ.0) THEN
 269           SXX=0.E0; SYY=0.E0; SZZ=0.E0
 270           SXY=0.E0; SXZ=0.E0; SYZ=0.E0
 271         ELSE
 272           SXX=EXX(I_X,I_Y,I_Z)
 273           SYY=EYY(I_X,I_Y,I_Z)
 274           SZZ=EZZ(I_X,I_Y,I_Z)
 275           SXY=EXY(I_X,I_Y,I_Z)
 276           SXZ=EXZ(I_X,I_Y,I_Z)
 277           SYZ=EYZ(I_X,I_Y,I_Z)
 278         END IF
 279         IF(PZO_Action.EQ.0.E0) THEN
 280           PZO=0.E0
 281         ELSE
 282           PZO=POT(I_X,I_Y,I_Z)
 283         END IF
 284
 285          HKANE(1,1)=POTHH(CHI)-PZO+ &
 286                      (PD1(CHI)+PD3(CHI))*SZZ+(PD2(CHI)+PD4(CHI))*(SXX+SYY)
 287
 288          HKANE(2,2)=POTLH(CHI)-PZO+ &
 289                      (PD1(CHI)+PD3(CHI))*SZZ+(PD2(CHI)+PD4(CHI))*(SXX+SYY)
 290
 291          HKANE(3,3)=POTSO(CHI)-PZO+ &
 292                      (PD1(CHI))*SZZ+(PD2(CHI))*(SXX+SYY)
 293
 294          HKANE(4,4)=HKANE(1,1)
 295          HKANE(5,5)=HKANE(2,2)
 296          HKANE(6,6)=HKANE(3,3)
 297
 298          IF(DIAG_Action.NE.0) THEN
 299
 300           HKANE(1,2)=-PD5(CHI)*CMPLX(SXX-SYY,-2.E0*SXY)
 301           HKANE(1,3)=-PD6(CHI)*CMPLX(SXZ,-SYZ)
 302           HKANE(1,4)=0.E0
 303           HKANE(1,5)=0.E0
 304           HKANE(1,6)=0.E0
 305
 306           HKANE(2,3)=PD6(CHI)*CMPLX(SXZ,+SYZ)
 307           HKANE(2,4)=0.E0
 308           HKANE(2,5)=0.E0
 309           HKANE(2,6)= SQRT(2.E0)*POTLS(CHI)
 310
 311           HKANE(3,4)= 0.E0
 312           HKANE(3,5)= SQRT(2.E0)*POTLS(CHI)
 313           HKANE(3,6)= 0.E0
 314
 315           HKANE(4,5)=-PD5(CHI)*CMPLX(SXX-SYY,2.E0*SXY)
 316           HKANE(4,6)=HKANE(2,3)
 317           HKANE(5,6)=HKANE(1,3)
 318
 319
 320
 321           DO I=1,DIMM
 322            DO J=I+1,DIMM
 323             HKANE(J,I)=CONJG(HKANE(I,J))
 324            END DO
 325           END DO
 326
 327           CALL CHEEV('V','U',DIMM,HKANE,DIMM,W,WORK,LWORK,RWORK,INFO)
 328
 329           IF (INFO.ne.0) then
 330                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 331                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
 332 !              STOP
 333           END IF
 334
 335 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 336 !!   Ordering Eigenvalues according to the Bloch func. caracter
 337 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 338
 339           ROOT = 0.E0
 340
 341           DO I=1,3
 342 ! Calculo de la componentes dependientes del spin
 343            DO J=1,6
 344             CARAC(J)=HKANE(J,I)*HKANE(J,I)
 345            END DO
 346 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 347            CARAC(1)=CARAC(1)+CARAC(4)
 348            CARAC(4)=CARAC(1)
 349            CARAC(2)=CARAC(2)+CARAC(5)
 350            CARAC(5)=CARAC(2)
 351            CARAC(3)=CARAC(3)+CARAC(6)
 352            CARAC(6)=CARAC(3)
 353 !
 354 ! Peso de las componentes de las que se extraeran los autovalores.
 355 !
 356            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 357                    ROOT(4)=W(I)
 358            END IF
 359            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 360                    ROOT(3)=W(I)
 361            END IF
 362            IF (CARAC(3).GE.CARAC_MAX) THEN ! SO
 363                    ROOT(1)=W(I)
 364            END IF
 365
 366           END DO
 367
 368          ELSE ! Only the diagonal elements are calculated
 369
 370           ROOT(4)=HKANE(1,1)
 371           ROOT(3)=HKANE(2,2)
 372           ROOT(1)=HKANE(3,3)
 373           ROOT(2)=HKANE(6,6) !Redundant
 374
 375
 376          END IF
 377
 378 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 379 !! RESULTS
 380 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 381
 382          ENERGY(1) = POTE(CHI)-PZO+(PC1(CHI)*SZZ+PC2(CHI)*(SXX+SYY))
 383
 384          ENERGY(2) = ROOT(4)
 385          ENERGY(3) = ROOT(3)
 386          ENERGY(4) = ROOT(3)
 387          ENERGY(5) = ROOT(4)
 388
 389          ENERGY(6) = ROOT(1)
 390          ENERGY(7) = ROOT(1)
 391 !        ENERGY(7) = CHI ! To save the Structure-profile
 392
 393 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 394 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 395 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
 396 !!!!! we will use this array to pack the Elastic Energy.
 397 !!!!!
 398 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
 399 !!!!!
 400 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 401 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 402
 403           ENERGY(8) = C11*(SXX**2 + SYY**2)+C33*SZZ**2+ &
 404                       2.E0*(XLAMB*(SXX*SYY)+C13*(SXX+SYY)*SZZ+ &
 405                       (C11-XLAMB)*SXY**2+2.E0*XMU*(SXZ**2+SYZ**2))
 406
 407           ENERGY(8) = ENERGY(8)/2.E0
 408
 409           RETURN
 410
 411       END SUBROUTINE POT_CALC_CAR_CH_WZ
 412
 413       SUBROUTINE POT_CALC_CYL_CH_WZ(CHI,ENERGY)
 414       IMPLICIT NONE
 415
 416         REAL, DIMENSION(:) :: ENERGY
 417         INTEGER, PARAMETER :: DIMM=6
 418         REAL, DIMENSION(1:DIMM,1:DIMM) :: HKANE
 419         REAL :: PZO
 420         REAL :: SSUM,SDIF,SZZ,SRZ,SRR,S00
 421         INTEGER :: I,J,CHI
 422 !! DIAGONALIZATION VARIABLES
 423
 424         REAL ::    W(1:DIMM), WORK(1:10*DIMM), AW(1:DIMM,1:DIMM), &
 425                    ROOT(DIMM), CARAC(DIMM), SLAMCH
 426         INTEGER :: IWORK(DIMM*5), IFAIL(DIMM), INFO, NUM
 427
 428         IF(STR_Action.EQ.0) THEN
 429           SSUM=0.E0; SDIF=0.E0; SZZ=0.E0;
 430           SRZ=0.E0; SRR=0.E0; S00=0.E0
 431         ELSE
 432
 433 !!! FOR CYLINDRICAL COORDINATES WE CALCULATE ONLY THE X-Z PLANE,
 434 !!! THAT MEANS: THETA=0
 435
 436           SRR=EXX(I_X,I_Y,I_Z)
 437           S00=EYY(I_X,I_Y,I_Z)
 438           SSUM=EXX(I_X,I_Y,I_Z)+EYY(I_X,I_Y,I_Z)
 439           SDIF=EXX(I_X,I_Y,I_Z)-EYY(I_X,I_Y,I_Z)
 440           SZZ=EZZ(I_X,I_Y,I_Z)
 441           SRZ=EXZ(I_X,I_Y,I_Z)
 442         END IF
 443         IF(PZO_Action.EQ.0.E0) THEN
 444           PZO=0.E0
 445         ELSE
 446           PZO=POT(I_X,I_Y,I_Z)
 447         END IF
 448
 449          HKANE(1,1)=POTHH(CHI)-PZO+ &
 450                      (PD1(CHI)+PD3(CHI))*SZZ+(PD2(CHI)+PD4(CHI))*SSUM
 451
 452          HKANE(2,2)=POTLH(CHI)-PZO+ &
 453                      (PD1(CHI)+PD3(CHI))*SZZ+(PD2(CHI)+PD4(CHI))*SSUM
 454
 455          HKANE(3,3)=POTSO(CHI)-PZO+ &
 456                      PD1(CHI)*SZZ+PD2(CHI)*SSUM
 457
 458          HKANE(4,4)=HKANE(1,1)
 459          HKANE(5,5)=HKANE(2,2)
 460          HKANE(6,6)=HKANE(3,3)
 461
 462          IF(DIAG_Action.NE.0) THEN
 463
 464
 465           HKANE(1,2)=-PD5(CHI)*SDIF
 466           HKANE(1,3)=-PD6(CHI)*SRZ
 467           HKANE(1,4)=0.E0
 468           HKANE(1,5)=0.E0
 469           HKANE(1,6)=0.E0
 470
 471           HKANE(2,3)=PD6(CHI)*SRZ
 472           HKANE(2,4)=0.E0
 473           HKANE(2,5)=0.E0
 474           HKANE(2,6)= SQRT(2.E0)*POTLS(CHI)
 475
 476           HKANE(3,4)= 0.E0
 477           HKANE(3,5)= SQRT(2.E0)*POTLS(CHI)
 478           HKANE(3,6)= 0.E0
 479
 480           HKANE(4,5)=-PD5(CHI)*SDIF
 481           HKANE(4,6)=HKANE(2,3)
 482           HKANE(5,6)=HKANE(1,3)
 483
 484           DO I=1,6
 485            DO J=I+1,6
 486             HKANE(J,I)=HKANE(I,J)
 487            END DO
 488           END DO
 489
 490 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 491 !!   Analytic solutions for Rho = 0.
 492 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 493 !
 494 !         ZETA = ZMIN+ZINC*IZ ! XZ normalized to ZC
 495 !
 496 !        W(1)=(HKANE(2,2)+HKANE(5,5))+
 497 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 498 !        W(2)=(HKANE(2,2)+HKANE(5,5))-
 499 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 500 !        write(24,'(4(f18.8,1x))')ZETA,HKANE(1,1),W(1)/2.,W(2)/2.
 501 !
 502 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 503
 504
 505           CALL DSYEVX('V','A','U',DIMM,HKANE,DIMM,0.,0.,2,2, &
 506            2*SLAMCH('S'),NUM,W,AW,DIMM,WORK,10*DIMM,IWORK,IFAIL,INFO)
 507
 508           IF (INFO.ne.0) then
 509                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 510                write(6,'(A,10(I2,1X))')"IFAIL=",IFAIL(1:DIMM)
 511                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
 512 !              STOP
 513           END IF
 514
 515 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 516 !!   Ordering Eigenvalues according to the Bloch func. caracter
 517 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 518
 519           ROOT = 0.E0
 520           WRITE(6,*)W(1:6)
 521
 522           DO I=1,6
 523 ! Calculo de la componentes dependientes del spin
 524            DO J=1,6
 525             CARAC(J)=AW(J,I)*AW(J,I)
 526            END DO
 527 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 528            CARAC(1)=CARAC(1)+CARAC(4)
 529            CARAC(4)=CARAC(1)
 530            CARAC(2)=CARAC(2)+CARAC(5)
 531            CARAC(5)=CARAC(2)
 532            CARAC(3)=CARAC(3)+CARAC(6)
 533            CARAC(6)=CARAC(3)
 534 !
 535 ! Peso de las componentes de las que se extraeran los autovalores.
 536 !
 537            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 538                    ROOT(4)=W(I)
 539            END IF
 540            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 541                    ROOT(3)=W(I)
 542            END IF
 543            IF (CARAC(3).GE.CARAC_MAX) THEN ! SO
 544                    ROOT(1)=W(I)
 545            END IF
 546
 547           END DO
 548
 549          ELSE ! Only the diagonal elements are calculated
 550
 551           ROOT(4)=HKANE(1,1)
 552           ROOT(3)=HKANE(2,2)
 553           ROOT(1)=HKANE(3,3)
 554           ROOT(2)=HKANE(6,6) !Redundant
 555
 556          END IF
 557
 558 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 559 !! RESULTS
 560 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 561
 562          ENERGY(1) = POTE(CHI)-PZO+(PC1(CHI)*SZZ+PC2(CHI)*SSUM)
 563
 564          ENERGY(2) = ROOT(4)
 565          ENERGY(3) = ROOT(3)
 566          ENERGY(4) = ROOT(3)
 567          ENERGY(5) = ROOT(4)
 568
 569          ENERGY(6) = ROOT(1)
 570          ENERGY(7) = ROOT(1)
 571 !        ESODW(IR,IZ,IGEO) = CHI ! To save the Structure-profile
 572
 573 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 574 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 575 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
 576 !!!!! we will use this array to pack the Elastic Energy.
 577 !!!!!
 578 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
 579 !!!!!
 580 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 581 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 582
 583
 584           ENERGY(8) = C11*(SRR**2 + S00**2)+C33*SZZ**2+ &
 585                       2.E0*(XLAMB*(SRR*S00)+C13*(SRR+S00)*SZZ+ &
 586                       +4.E0*XMU*(SRZ**2))
 587
 588           ENERGY(8) = ENERGY(8)/2.E0
 589
 590           RETURN
 591
 592       END SUBROUTINE POT_CALC_CYL_CH_WZ
 593
 594       SUBROUTINE POT_CALC_CAR_WZ(CHI,ENERGY)
 595       IMPLICIT NONE
 596
 597         REAL, DIMENSION(:) :: ENERGY
 598         INTEGER, PARAMETER :: DIMM=6
 599         COMPLEX, DIMENSION(1:DIMM,1:DIMM) :: HKANE
 600         REAL :: SXX,SYY,SZZ,SXY,SXZ,SYZ,PZO
 601         INTEGER :: I,J,CHI
 602 !! DIAGONALIZATION VARIABLES
 603
 604         COMPLEX, DIMENSION(1:(2*DIMM-1)) :: WORK
 605         REAL, DIMENSION(1:(3*DIMM-2)) :: RWORK
 606         REAL         W(1:DIMM), ROOT(DIMM), CARAC(DIMM)
 607         INTEGER      INFO, LWORK
 608
 609         LWORK=2*DIMM-1
 610
 611         IF(STR_Action.EQ.0) THEN
 612           SXX=0.E0; SYY=0.E0; SZZ=0.E0
 613           SXY=0.E0; SXZ=0.E0; SYZ=0.E0
 614         ELSE
 615           SXX=EXX(I_X,I_Y,I_Z)
 616           SYY=EYY(I_X,I_Y,I_Z)
 617           SZZ=EZZ(I_X,I_Y,I_Z)
 618           SXY=EXY(I_X,I_Y,I_Z)
 619           SXZ=EXZ(I_X,I_Y,I_Z)
 620           SYZ=EYZ(I_X,I_Y,I_Z)
 621         END IF
 622         IF(PZO_Action.EQ.0.E0) THEN
 623           PZO=0.E0
 624         ELSE
 625           PZO=POT(I_X,I_Y,I_Z)
 626         END IF
 627
 628          HKANE(1,1)=POTHH(CHI)-PZO+ &
 629                      (PD1(CHI)+PD3(CHI))*SZZ+(PD2(CHI)+PD4(CHI))*(SXX+SYY)
 630
 631          HKANE(2,2)=POTLH(CHI)-PZO+(PD1(CHI)+PD3(CHI)/3.E0)*SZZ+ &
 632                      (PD2(CHI)+PD4(CHI)/3.E0)*(SXX+SYY)
 633          HKANE(3,3)=HKANE(2,2)
 634          HKANE(4,4)=HKANE(1,1)
 635          HKANE(5,5)=POTSO(CHI)-PZO+(PD1(CHI)+2.E0*PD3(CHI)/3.E0)*SZZ+ &
 636                      (PD2(CHI)+2.E0*PD4(CHI)/3.E0)*(SXX+SYY)
 637          HKANE(6,6)=HKANE(5,5)
 638
 639          IF(DIAG_Action.NE.0) THEN
 640
 641           HKANE(1,2)=-SQRT(2.E0/3.E0)*PD6(CHI)*CMPLX(SXZ,-SYZ)
 642           HKANE(1,3)=-1.E0/SQRT(3.E0)*PD5(CHI)*CMPLX(SXX-SYY,-2.E0*SXY)
 643           HKANE(1,4)=0.E0
 644           HKANE(1,5)=1.E0/SQRT(3.E0)*PD6(CHI)*CMPLX(SXZ,-SYZ)
 645           HKANE(1,6)=SQRT(2./3.)*PD5(CHI)*CMPLX(SXX-SYY,-2.E0*SXY)
 646
 647           HKANE(2,3)=0.
 648           HKANE(2,4)=HKANE(1,3)
 649           HKANE(2,5)=SQRT(2.E0)/3.E0*( POTLS(CHI)+&
 650                                        (PD3(CHI)*SZZ+PD4(CHI)*(SXX+SYY)) )
 651           HKANE(2,6)=-PD6(CHI)*CMPLX(SXZ,-SYZ)
 652
 653           HKANE(3,4)=-HKANE(1,2)
 654           HKANE(3,5)=-PD6(CHI)*CMPLX(SXZ,SYZ)
 655           HKANE(3,6)=-HKANE(2,5)
 656
 657           HKANE(4,5)=-SQRT(2./3.)*PD5(CHI)*CMPLX(SXX-SYY,+2.E0*SXY)
 658           HKANE(4,6)=PD6(CHI)*CMPLX(SXZ,SYZ)/SQRT(3.E0)
 659
 660           HKANE(5,6)=0
 661
 662
 663
 664           DO I=1,DIMM
 665            DO J=I+1,DIMM
 666             HKANE(J,I)=CONJG(HKANE(I,J))
 667            END DO
 668           END DO
 669
 670           CALL CHEEV('V','U',DIMM,HKANE,DIMM,W,WORK,LWORK,RWORK,INFO)
 671
 672           IF (INFO.ne.0) then
 673                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 674                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
 675 !              STOP
 676           END IF
 677
 678 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 679 !!   Ordering Eigenvalues according to the Bloch func. caracter
 680 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 681
 682           ROOT = 0.E0
 683
 684           DO I=1,6,2
 685 ! Calculo de la componentes dependientes del spin
 686            DO J=1,6
 687             CARAC(J)=HKANE(J,I)*HKANE(J,I)
 688            END DO
 689 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 690            CARAC(1)=CARAC(1)+CARAC(4)
 691            CARAC(4)=CARAC(1)
 692            CARAC(2)=CARAC(2)+CARAC(3)
 693            CARAC(3)=CARAC(2)
 694            CARAC(5)=CARAC(5)+CARAC(6)
 695            CARAC(6)=CARAC(5)
 696 !
 697 ! Peso de las componentes de las que se extraeran los autovalores.
 698 !
 699            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 700                    ROOT(4)=W(I)
 701            END IF
 702            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 703                    ROOT(3)=W(I)
 704            END IF
 705            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
 706                    ROOT(1)=W(I)
 707            END IF
 708
 709           END DO
 710
 711          ELSE ! Only the diagonal elements were calculated
 712
 713           ROOT(4)=HKANE(1,1)
 714           ROOT(3)=HKANE(2,2)
 715           ROOT(1)=HKANE(5,5)
 716           ROOT(2)=HKANE(6,6) !Redundant
 717
 718
 719          END IF
 720
 721 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 722 !! RESULTS
 723 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 724
 725          ENERGY(1) = POTE(CHI)-PZO+(PC1(CHI)*SZZ+PC2(CHI)*(SXX+SYY))
 726
 727          ENERGY(2) = ROOT(4)
 728          ENERGY(3) = ROOT(3)
 729          ENERGY(4) = ROOT(3)
 730          ENERGY(5) = ROOT(4)
 731
 732          ENERGY(6) = ROOT(1)
 733          ENERGY(7) = ROOT(1)
 734 !        ENERGY(7) = CHI ! To save the Structure-profile
 735
 736 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 737 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 738 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
 739 !!!!! we will use this array to pack the Elastic Energy.
 740 !!!!!
 741 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
 742 !!!!!
 743 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 744 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 745
 746           ENERGY(8) = C11*(SXX**2 + SYY**2)+C33*SZZ**2+ &
 747                       2.E0*(XLAMB*(SXX*SYY)+C13*(SXX+SYY)*SZZ+ &
 748                       (C11-XLAMB)*SXY**2+2.E0*XMU*(SXZ**2+SYZ**2))
 749
 750           ENERGY(8) = ENERGY(8)/2.E0
 751
 752           RETURN
 753
 754       END SUBROUTINE POT_CALC_CAR_WZ
 755
 756       SUBROUTINE POT_CALC_CYL_WZ(CHI,ENERGY)
 757       IMPLICIT NONE
 758
 759         REAL, DIMENSION(:) :: ENERGY
 760         INTEGER, PARAMETER :: DIMM=6
 761         REAL, DIMENSION(1:DIMM,1:DIMM) :: HKANE
 762         REAL :: PZO
 763         REAL :: SSUM,SDIF,SZZ,SRZ,SRR,S00
 764         INTEGER :: I,J,CHI
 765 !! DIAGONALIZATION VARIABLES
 766
 767         REAL ::    W(1:DIMM), WORK(1:10*DIMM), AW(1:DIMM,1:DIMM), &
 768                    ROOT(DIMM), CARAC(DIMM), SLAMCH
 769         INTEGER :: IWORK(DIMM*5), IFAIL(DIMM), INFO, NUM
 770
 771         IF(STR_Action.EQ.0) THEN
 772           SSUM=0.E0; SDIF=0.E0; SZZ=0.E0;
 773           SRZ=0.E0; SRR=0.E0; S00=0.E0
 774         ELSE
 775
 776 !!! FOR CYLINDRICAL COORDINATES WE CALCULATE ONLY THE X-Z PLANE,
 777 !!! THAT MEANS: THETA=0
 778
 779           SRR=EXX(I_X,I_Y,I_Z)
 780           S00=EYY(I_X,I_Y,I_Z)
 781           SSUM=EXX(I_X,I_Y,I_Z)+EYY(I_X,I_Y,I_Z)
 782           SDIF=EXX(I_X,I_Y,I_Z)-EYY(I_X,I_Y,I_Z)
 783           SZZ=EZZ(I_X,I_Y,I_Z)
 784           SRZ=EXZ(I_X,I_Y,I_Z)
 785         END IF
 786         IF(PZO_Action.EQ.0.E0) THEN
 787           PZO=0.E0
 788         ELSE
 789           PZO=POT(I_X,I_Y,I_Z)
 790         END IF
 791
 792          HKANE(1,1)=POTHH(CHI)-PZO+ &
 793                      (PD1(CHI)+PD3(CHI))*SZZ+(PD2(CHI)+PD4(CHI))*SSUM
 794
 795          HKANE(2,2)=POTLH(CHI)-PZO+(PD1(CHI)+PD3(CHI)/3.E0)*SZZ+ &
 796                      (PD2(CHI)+PD4(CHI)/3.E0)*SSUM
 797          HKANE(3,3)=HKANE(2,2)
 798          HKANE(4,4)=HKANE(1,1)
 799          HKANE(5,5)=POTSO(CHI)-PZO+(PD1(CHI)+2.E0*PD3(CHI)/3.E0)*SZZ+ &
 800                      (PD2(CHI)+2.E0*PD4(CHI)/3.E0)*SSUM
 801          HKANE(6,6)=HKANE(5,5)
 802
 803          IF(DIAG_Action.NE.0) THEN
 804
 805           HKANE(1,2)=-SQRT(2.E0/3.E0)*PD6(CHI)*SRZ
 806           HKANE(1,3)=-1.E0/SQRT(3.E0)*PD5(CHI)*SDIF
 807           HKANE(1,4)=0.
 808           HKANE(1,5)=1.E0/SQRT(3.E0)*PD6(CHI)*SRZ
 809           HKANE(1,6)=SQRT(2./3.)*PD5(CHI)*SDIF
 810
 811           HKANE(2,3)=0.
 812           HKANE(2,4)=HKANE(1,3)
 813           HKANE(2,5)=SQRT(2.E0)/3.E0*(POTLS(CHI)+(PD3(CHI)*SZZ+PD4(CHI)*SSUM))
 814           HKANE(2,6)=-PD6(CHI)*SRZ
 815
 816           HKANE(3,4)=-HKANE(1,2)
 817           HKANE(3,5)=HKANE(2,6)
 818           HKANE(3,6)=-HKANE(2,5)
 819
 820           HKANE(4,5)=-HKANE(1,6)
 821           HKANE(4,6)=HKANE(1,5)
 822
 823           HKANE(5,6)=0
 824
 825
 826
 827           DO I=1,6
 828            DO J=I+1,6
 829             HKANE(J,I)=HKANE(I,J)
 830            END DO
 831           END DO
 832
 833 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 834 !!   Analytic solutions for Rho = 0.
 835 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 836 !
 837 !         ZETA = ZMIN+ZINC*IZ ! XZ normalized to ZC
 838 !
 839 !        W(1)=(HKANE(2,2)+HKANE(5,5))+
 840 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 841 !        W(2)=(HKANE(2,2)+HKANE(5,5))-
 842 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
 843 !        write(24,'(4(f18.8,1x))')ZETA,HKANE(1,1),W(1)/2.,W(2)/2.
 844 !
 845 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 846
 847
 848           CALL DSYEVX('V','A','U',DIMM,HKANE,DIMM,0.,0.,2,2, &
 849            2*SLAMCH('S'),NUM,W,AW,DIMM,WORK,10*DIMM,IWORK,IFAIL,INFO)
 850
 851           IF (INFO.ne.0) then
 852                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
 853                write(6,'(A,10(I2,1X))')"IFAIL=",IFAIL(1:DIMM)
 854                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
 855 !              STOP
 856           END IF
 857
 858 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 859 !!   Ordering Eigenvalues according to the Bloch func. caracter
 860 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 861
 862           ROOT = 0.E0
 863
 864           DO I=1,6,2
 865 ! Calculo de la componentes dependientes del spin
 866            DO J=1,6
 867             CARAC(J)=AW(J,I)*AW(J,I)
 868            END DO
 869 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
 870            CARAC(1)=CARAC(1)+CARAC(4)
 871            CARAC(4)=CARAC(1)
 872            CARAC(2)=CARAC(2)+CARAC(3)
 873            CARAC(3)=CARAC(2)
 874            CARAC(5)=CARAC(5)+CARAC(6)
 875            CARAC(6)=CARAC(5)
 876 !
 877 ! Peso de las componentes de las que se extraeran los autovalores.
 878 !
 879            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
 880                    ROOT(4)=W(I)
 881            END IF
 882            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
 883                    ROOT(3)=W(I)
 884            END IF
 885            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
 886                    ROOT(1)=W(I)
 887            END IF
 888
 889           END DO
 890
 891          ELSE ! Only the diagonal elements were calculated
 892
 893           ROOT(4)=HKANE(1,1)
 894           ROOT(3)=HKANE(2,2)
 895           ROOT(1)=HKANE(5,5)
 896           ROOT(2)=HKANE(6,6) !Redundant
 897
 898
 899          END IF
 900
 901 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 902 !! RESULTS
 903 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 904
 905          ENERGY(1) = POTE(CHI)-PZO+(PC1(CHI)*SZZ+PC2(CHI)*SSUM)
 906
 907          ENERGY(2) = ROOT(4)
 908          ENERGY(3) = ROOT(3)
 909          ENERGY(4) = ROOT(3)
 910          ENERGY(5) = ROOT(4)
 911
 912          ENERGY(6) = ROOT(1)
 913          ENERGY(7) = ROOT(1)
 914 !        ESODW(IR,IZ,IGEO) = CHI ! To save the Structure-profile
 915
 916 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 917 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 918 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
 919 !!!!! we will use this array to pack the Elastic Energy.
 920 !!!!!
 921 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
 922 !!!!!
 923 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 924 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 925
 926
 927           ENERGY(8) = C11*(SRR**2 + S00**2)+C33*SZZ**2+ &
 928                       2.E0*(XLAMB*(SRR*S00)+C13*(SRR+S00)*SZZ+ &
 929                       +4.E0*XMU*(SRZ**2))
 930
 931           ENERGY(8) = ENERGY(8)/2.E0
 932
 933           RETURN
 934
 935       END SUBROUTINE POT_CALC_CYL_WZ
 936
 937     END SUBROUTINE POTENTIAL_WZ
 938
 939     SUBROUTINE POTENTIAL_ZB(EEL,EHHUP,EHHDW,ELHUP,ELHDW,&
 940                             ESOUP,ESODW,ELAST,EXX,EYY,&
 941                             EZZ,EXY,EXZ,EYZ)
 942
 943
 944               Use Input_Data
 945               Use Dot_Geometry
 946               Use Auxiliar_Procedures, ONLY : AISO
 947
 948       IMPLICIT NONE
 949
 950 !!!!! 'dummy' and local variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 951
 952       REAL,DIMENSION(:,:,:),OPTIONAL ::  EXX,EYY,EZZ,EXY,EXZ,EYZ
 953
 954       REAL,DIMENSION(:,:,:) ::  EEL,EHHUP,EHHDW,ELHUP,ELHDW,&
 955                                 ESOUP,ESODW,ELAST
 956
 957       REAL         ZM,THETA,CTHETA,STHETA, &
 958                    X,Y,Z,ZMAUX,RHO,ZETA
 959
 960       INTEGER      I_X,I_Y,I_Z,I_N1,I_N2,I_N3,CHI
 961
 962       REAL, DIMENSION(3) :: R_SL,X_VEC,XI_VEC
 963
 964       REAL :: POTWE,POTWHH,POTWLH,POTWSO,&
 965               VBIEL,VBIHH,VBILH,VBISO
 966
 967       REAL,DIMENSION(0:1) :: POTE, POTHH, POTLH, POTSO, &
 968                              DVD,DSD,DVU,DVSU,D2VU,D2VSU,DSO,ZBC1
 969
 970       REAL, DIMENSION(1:8) :: ENERGY
 971
 972 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 973
 974       CALL ENER_CONSTANTS_ZB( )
 975
 976 ZD:   DO I_Z=1,ZDim
 977 !       WRITE(16,'(A,I3,A,I3)')"I_Z ",I_Z," of ",ZDIM
 978         Z=Z_Min+REAL(I_Z-1)*Z_Inc
 979 YD:   DO I_Y=1,YDim
 980         Y=Y_Min+REAL(I_Y-1)*Y_Inc
 981 XD:   DO I_X=1,XDim
 982         X=X_Min+REAL(I_X-1)*X_Inc
 983
 984         X_VEC=(/X,Y,Z/)
 985
 986       I_N1=0; I_N2=0; I_N3=0
 987
 988       R_SL=REAL(I_N1)*A1_S+REAL(I_N2)*A2_S+REAL(I_N3)*A3_S
 989
 990       XI_VEC=X_VEC-R_SL
 991
 992       RHO=SQRT(XI_VEC(1)**2+XI_VEC(2)**2)/RC
 993       IF(XI_VEC(1).EQ.0.E0.AND.XI_VEC(2).EQ.0.E0) THEN
 994               CTHETA=1.E0/SQRT(2.E0); STHETA=1.E0/SQRT(2.E0) ! It is not the Mathematical limit
 995       ELSE
 996               THETA=ATAN(XI_VEC(2)/XI_VEC(1))
 997               CTHETA=Cos(THETA); STHETA=Sin(THETA)
 998       END IF
 999       ZETA=XI_VEC(3)/ZC
1000
1001           IF (RHO.LE.RD) THEN
1002              CALL SHAPERTOZ(MIN(RHO*RC,Rqd_Base),ZMAUX)
1003              ZM=ZMAUX/ZC
1004           ELSE
1005              ZM = 0.E0
1006           END IF
1007
1008           IF (abs(zeta) .EQ. 0.E0 .OR. ZETA .EQ. ZM) THEN
1009               ZETA=ZETA-1.E-5
1010           END IF
1011
1012           CHI = 0
1013           IF (ZETA.GE.-D.AND.ZETA.LE.ZM) THEN
1014                   CHI = 1
1015                   IF(I_N1.NE.0.OR.I_N2.NE.0.OR.I_N3.NE.0) THEN
1016                           WRITE(16,*)I_N1,I_N2,I_N3
1017                           WRITE(16,*)X_VEC(3),ZETA*ZC,ZM*ZC
1018                   END IF
1019           END IF
1020
1021           IF (KCOOR.EQ.0) THEN
1022               CALL POT_CALC_CAR_ZB(CHI,ENERGY)
1023           ELSE
1024               CALL POT_CALC_CYL_ZB(CHI,ENERGY)
1025           END IF
1026
1027            EEL(I_X,I_Y,I_Z)   = ENERGY(1)
1028            EHHUP(I_X,I_Y,I_Z) = ENERGY(2)
1029            ELHUP(I_X,I_Y,I_Z) = ENERGY(3)
1030            ELHDW(I_X,I_Y,I_Z) = ENERGY(4)
1031            EHHDW(I_X,I_Y,I_Z) = ENERGY(5)
1032            ESOUP(I_X,I_Y,I_Z) = ENERGY(6)
1033            ESODW(I_X,I_Y,I_Z) = ENERGY(7)
1034            ELAST(I_X,I_Y,I_Z) = ENERGY(8)
1035
1036 !           WRITE(26,'(10(E15.8,1X))')Z,ENERGY(1:8)
1037
1038       END DO XD
1039       END DO YD
1040       END DO ZD
1041
1042 !     STOP
1043
1044       RETURN
1045
1046       CONTAINS
1047
1048       SUBROUTINE ENER_CONSTANTS_ZB( )
1049       IMPLICIT NONE
1050
1051 ! Definition of parameters for Barrier and Well
1052
1053 ! If we set both terms equal to zero the spin-interaction in the
1054 ! deformation potentials is removed. The values of AVB,BB,DB are
1055 ! taken from C. Pryor, PRB, 57, 7190 (1998)
1056
1057 !        DSO(0) = VBH - VBSO ; DSO(1) = VWH - VWSO
1058          DSO(0) = 0.0; DSO(1) = 0.0
1059
1060          ZBC1(0) = -9.3 ; ZBC1(1) = ACW
1061
1062          DVD(0) = 0.7E0 - 2./9. * DSO(0)
1063          DVD(1) = AVW - 2./9. * DSO(1)
1064          DSD(0) = 0.7E0 + 4./9. * DSO(0)
1065          DSD(1) = AVW + 4./9. * DSO(1)
1066          DVU(0) = -3./2. * (-2.0) + 1./3. * DSO(0)
1067          DVU(1) = -3./2. * BW + 1./3. * DSO(1)
1068          DVSU(0) = -3./2. * (-2.0) - 1./6. * DSO(0)
1069          DVSU(1) = -3./2. * BW - 1./6. * DSO(1)
1070          D2VU(0) = -SQRT(3.)/2. * (-5.4) + 1./3. * DSO(0)
1071          D2VU(1) = -SQRT(3.)/2. * DW + 1./3. * DSO(1)
1072          D2VSU(0) = -SQRT(3.)/2. * (-5.4) - 1./6. * DSO(0)
1073          D2VSU(1) = -SQRT(3.)/2. * DW - 1./6. * DSO(1)
1074
1075 !!! In the calculation the deformation potentials are equal across the structure
1076          DVD(0)=DVD(1); DSD(0)=DSD(1); DVU(0)=DVU(1)
1077          DVSU(0)=DVSU(1); D2VU(0)=D2VU(1); D2VSU(0)=D2VSU(1)
1078
1079          DSO(0) = VBH - VBSO ; DSO(1) = VWH - VWSO
1080          POTWE= VWE
1081          POTWHH=VWH
1082          POTWLH=VWH
1083          POTWSO=VWSO
1084
1085           IF(STR_Action.EQ.0) THEN
1086            VBIEL  = ACW*(BISUM+BIZZ)
1087            VBIHH  = AVW*(BISUM+BIZZ)-(-3.E0*BW/2.E0)/3.E0*(BISUM-2.E0*BIZZ)
1088            VBILH  = AVW*(BISUM+BIZZ)+(-3.E0*BW/2.E0)/3.E0*(BISUM-2.E0*BIZZ)
1089            VBISO  = AVW*(BISUM+BIZZ)
1090
1091 !          Potential edges including strain effect
1092            POTWE= (VWE+VBIEL)
1093            POTWHH= (VWH+VBIHH)
1094            POTWLH= (VWH+VBILH)
1095            POTWSO= (VWSO+VBISO)
1096           END IF
1097
1098          POTE(0)  = VBE  ; POTE(1)  = POTWE
1099          POTHH(0) = VBH ;  POTHH(1) = POTWHH
1100          POTLH(0) = VBH ;  POTLH(1) = POTWLH
1101          POTSO(0) = VBSO ; POTSO(1) = POTWSO
1102
1103          RETURN
1104
1105       END SUBROUTINE ENER_CONSTANTS_ZB
1106
1107       SUBROUTINE POT_CALC_CAR_ZB(CHI,ENERGY)
1108       IMPLICIT NONE
1109
1110         REAL, DIMENSION(:) :: ENERGY
1111         INTEGER, PARAMETER :: DIMM=6
1112         COMPLEX, DIMENSION(1:DIMM,1:DIMM) :: HKANE
1113         REAL :: SXX,SYY,SZZ,SXY,SXZ,SYZ,SHID,SDIF,STIL
1114         INTEGER :: I,J,CHI
1115 !! DIAGONALIZATION VARIABLES
1116
1117         COMPLEX, DIMENSION(1:(2*DIMM-1)) :: WORK
1118         REAL, DIMENSION(1:(3*DIMM-2)) :: RWORK
1119         REAL         W(1:DIMM), ROOT(DIMM), CARAC(DIMM)
1120         INTEGER      INFO, LWORK
1121
1122
1123         IF(STR_Action.EQ.0) THEN
1124           SXX=0.E0; SYY=0.E0; SZZ=0.E0
1125           SXY=0.E0; SXZ=0.E0; SYZ=0.E0
1126         ELSE
1127           SXX=EXX(I_X,I_Y,I_Z)
1128           SYY=EYY(I_X,I_Y,I_Z)
1129           SZZ=EZZ(I_X,I_Y,I_Z)
1130           SXY=EXY(I_X,I_Y,I_Z)
1131           SXZ=EXZ(I_X,I_Y,I_Z)
1132           SYZ=EYZ(I_X,I_Y,I_Z)
1133           SHID=SXX+SYY+SZZ
1134           SDIF=SXX-SYY
1135           STIL=SXX+SYY-2.E0*SZZ
1136         END IF
1137
1138          HKANE(1,1)=POTHH(CHI)+DVD(CHI)*SHID &
1139                      -1./3.*DVU(CHI)*STIL
1140          HKANE(2,2)=POTLH(CHI)+DVD(CHI)*SHID &
1141                      +1./3.*DVU(CHI)*STIL
1142          HKANE(3,3)=HKANE(2,2)
1143          HKANE(4,4)=HKANE(1,1)
1144          HKANE(5,5)=POTSO(CHI)+DSD(CHI)*SHID
1145
1146          HKANE(6,6)=HKANE(5,5)
1147
1148          IF(DIAG_Action.NE.0.AND.STR_Action.NE.0) THEN
1149
1150            HKANE(1,2)=2./SQRT(3.)*D2VU(CHI)*CMPLX(SXZ,-SYZ)
1151            HKANE(1,3)=1./SQRT(3.)*CMPLX(DVU(CHI)*SDIF,-2.E0*D2VU(CHI)*SXY)
1152            HKANE(1,4)=0.E0
1153            HKANE(1,5)=-SQRT(2./3.)*D2VSU(CHI)*CMPLX(SXZ,-SYZ)
1154            HKANE(1,6)=-SQRT(2./3.)*CMPLX(DVU(CHI)*SDIF,-2.E0*D2VSU(CHI)*SXY)
1155
1156            HKANE(2,3)=0.
1157            HKANE(2,4)=HKANE(1,3)
1158            HKANE(2,5)=-SQRT(2.)/3.*DVSU(CHI)*STIL
1159            HKANE(2,6)=SQRT(2.)*D2VSU(CHI)*CMPLX(SXZ,-SYZ)
1160
1161            HKANE(3,4)=-HKANE(1,2)
1162            HKANE(3,5)=SQRT(2.)*D2VSU(CHI)*CMPLX(SXZ,SYZ)
1163            HKANE(3,6)=-HKANE(2,5)
1164
1165            HKANE(4,5)=SQRT(2./3.)*CMPLX(DVU(CHI)*SDIF,2.E0*D2VSU(CHI)*SXY)
1166            HKANE(4,6)=-HKANE(3,5)/SQRT(3.E0)
1167
1168            HKANE(5,6)=0
1169
1170          DO I=1,6
1171           DO J=I+1,6
1172            HKANE(J,I)=CONJG(HKANE(I,J))
1173           END DO
1174          END DO
1175
1176 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1177 !!   Analytic solutions for Rho = 0.
1178 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1179 !
1180 !         ZETA = ZMIN+ZINC*IZ ! XZ normalized to ZC
1181 !
1182 !        W(1)=(HKANE(2,2)+HKANE(5,5))+
1183 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
1184 !        W(2)=(HKANE(2,2)+HKANE(5,5))-
1185 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
1186 !        write(24,'(4(f18.8,1x))')ZETA,HKANE(1,1),W(1)/2.,W(2)/2.
1187 !
1188 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1189
1190           CALL CHEEV('V','U',DIMM,HKANE,DIMM,W,WORK,LWORK,RWORK,INFO)
1191
1192           IF (INFO.ne.0) then
1193                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
1194                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
1195 !              STOP
1196           END IF
1197
1198 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1199 !!   Ordering Eigenvalues according to the Bloch func. caracter
1200 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1201
1202           ROOT = 0.E0
1203
1204           DO I=1,6,2
1205 ! Calculo de la componentes dependientes del spin
1206            DO J=1,6
1207             CARAC(J)=HKANE(J,I)*HKANE(J,I)
1208            END DO
1209 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
1210            CARAC(1)=CARAC(1)+CARAC(4)
1211            CARAC(4)=CARAC(1)
1212            CARAC(2)=CARAC(2)+CARAC(3)
1213            CARAC(3)=CARAC(2)
1214            CARAC(5)=CARAC(5)+CARAC(6)
1215            CARAC(6)=CARAC(5)
1216 !
1217 ! Peso de las componentes de las que se extraeran los autovalores.
1218 !
1219            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
1220                    ROOT(4)=W(I)
1221            END IF
1222            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
1223                    ROOT(3)=W(I)
1224            END IF
1225            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
1226                    ROOT(1)=W(I)
1227            END IF
1228
1229           END DO
1230
1231          ELSE ! Only the diagonal elements were calculated
1232
1233           ROOT(4)=HKANE(1,1)
1234           ROOT(3)=HKANE(2,2)
1235           ROOT(1)=HKANE(5,5)
1236           ROOT(2)=HKANE(6,6) !Redundant
1237
1238
1239          END IF
1240
1241 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1242 !! RESULTS
1243 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1244
1245          ENERGY(1) = POTE(CHI)+ZBC1(CHI)*SHID
1246
1247          ENERGY(2) = ROOT(4)
1248          ENERGY(3) = ROOT(3)
1249          ENERGY(4) = ROOT(3)
1250          ENERGY(5) = ROOT(4)
1251
1252          ENERGY(6) = ROOT(1)
1253          ENERGY(7) = ROOT(1)
1254 !        ENERGY(7) = CHI ! To save the Structure-profile
1255
1256 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1257 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1258 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
1259 !!!!! we will use this array to pack the Elastic Energy.
1260 !!!!!
1261 !!!!!        U=1/2*(XLAMB*Tr(e)**2+XMU/2*(err**2+e00**2+ezz**2+erz**2)
1262 !!!!!
1263 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1264 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1265
1266           ENERGY(8) = XMU*(SXX**2 + SYY**2 + SZZ**2 + &
1267                            2.E0*(SXY**2+SXZ**2+SYZ**2) ) + &
1268                       XLAMB/2.E0*(SXX + SYY + SZZ)**2
1269           RETURN
1270
1271       END SUBROUTINE POT_CALC_CAR_ZB
1272
1273       SUBROUTINE POT_CALC_CYL_ZB(CHI,ENERGY)
1274       IMPLICIT NONE
1275
1276         REAL, DIMENSION(:) :: ENERGY
1277         INTEGER, PARAMETER :: DIMM=6
1278         REAL, DIMENSION(1:DIMM,1:DIMM) :: HKANE
1279         REAL :: SSUM,SDIF,SZZ,SRZ,SRR,S00,SHID,STIL
1280         INTEGER :: I,J,CHI
1281 !! DIAGONALIZATION VARIABLES
1282
1283         REAL ::    W(1:DIMM), WORK(1:10*DIMM), AW(1:DIMM,1:DIMM), &
1284                    ROOT(DIMM), CARAC(DIMM), SLAMCH
1285         INTEGER :: IWORK(DIMM*5), IFAIL(DIMM), INFO, NUM
1286
1287
1288         IF(STR_Action.EQ.0) THEN
1289           SSUM=0.E0; SDIF=0.E0; SZZ=0.E0;
1290           SRZ=0.E0; SRR=0.E0; S00=0.E0
1291         ELSE
1292           SRR=EXX(I_X,I_Y,I_Z)
1293           S00=EYY(I_X,I_Y,I_Z)
1294           SSUM=EXX(I_X,I_Y,I_Z)+EYY(I_X,I_Y,I_Z)
1295           SDIF=EXX(I_X,I_Y,I_Z)-EYY(I_X,I_Y,I_Z)
1296           SZZ=EZZ(I_X,I_Y,I_Z)
1297           SRZ=EXZ(I_X,I_Y,I_Z)
1298           SHID=SSUM+SZZ
1299           STIL=SSUM-2.E0*SZZ
1300         END IF
1301
1302          HKANE(1,1)=POTHH(CHI)+DVD(CHI)*SHID &
1303                     -1./3.*DVU(CHI)*STIL
1304          HKANE(2,2)=POTLH(CHI)+DVD(CHI)*SHID &
1305                     +1./3.*DVU(CHI)*STIL
1306
1307          HKANE(3,3)=HKANE(2,2)
1308          HKANE(4,4)=HKANE(1,1)
1309          HKANE(5,5)=POTSO(CHI)+DSD(CHI)*SHID
1310          HKANE(6,6)=HKANE(5,5)
1311
1312          IF(DIAG_Action.NE.0.AND.STR_Action.NE.0) THEN
1313             HKANE(1,2)=2./SQRT(3.)*D2VU(CHI)*SRZ
1314             HKANE(1,3)=1./SQRT(3.)*(DVU(CHI)+D2VU(CHI))/2.*SDIF
1315             HKANE(1,4)=0.E0
1316             HKANE(1,5)=-SQRT(2./3.)*D2VSU(CHI)*SRZ
1317             HKANE(1,6)=-SQRT(2./3.)*(DVU(CHI)+D2VSU(CHI))/2.*SDIF
1318
1319             HKANE(2,3)=0.E0
1320             HKANE(2,4)=HKANE(1,3)
1321             HKANE(2,5)=-SQRT(2.)/3.*DVSU(CHI)*STIL
1322             HKANE(2,6)=SQRT(2.)*D2VSU(CHI)*SRZ
1323
1324             HKANE(3,4)=-HKANE(1,2)
1325             HKANE(3,5)=HKANE(2,6)
1326             HKANE(3,6)=-HKANE(2,5)
1327
1328             HKANE(4,5)=-HKANE(1,6)
1329             HKANE(4,6)=HKANE(1,5)
1330
1331             HKANE(5,6)=0
1332
1333           DO I=1,6
1334            DO J=I+1,6
1335             HKANE(J,I)=HKANE(I,J)
1336            END DO
1337           END DO
1338
1339 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1340 !!   Analytic solutions for Rho = 0.
1341 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1342 !
1343 !         ZETA = ZMIN+ZINC*IZ ! XZ normalized to ZC
1344 !
1345 !        W(1)=(HKANE(2,2)+HKANE(5,5))+
1346 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
1347 !        W(2)=(HKANE(2,2)+HKANE(5,5))-
1348 !    & SQRT(  (HKANE(2,2)-HKANE(5,5))**2+4.*HKANE(3,6)**2  )
1349 !        write(24,'(4(f18.8,1x))')ZETA,HKANE(1,1),W(1)/2.,W(2)/2.
1350 !
1351 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1352
1353
1354           CALL DSYEVX('V','A','U',DIMM,HKANE,DIMM,0.,0.,2,2, &
1355            2*SLAMCH('S'),NUM,W,AW,DIMM,WORK,10*DIMM,IWORK,IFAIL,INFO)
1356
1357           IF (INFO.ne.0) then
1358                write(6,*)"Not Succesful Exit. Stopping. INFO=",INFO
1359                write(6,'(A,10(I2,1X))')"IFAIL=",IFAIL(1:DIMM)
1360                write(6,'(2(A,F6.3,1X))')"RHO=",RC*RHO,"ZETA=",ZETA*ZC
1361 !              STOP
1362           END IF
1363
1364 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1365 !!   Ordering Eigenvalues according to the Bloch func. caracter
1366 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1367
1368           ROOT = 0.E0
1369
1370           DO I=1,6,2
1371 ! Calculo de la componentes dependientes del spin
1372            DO J=1,6
1373             CARAC(J)=AW(J,I)*AW(J,I)
1374            END DO
1375 ! Suma de componentes independientemente del spin (hh_up + hh_dw, lh_up + lh_dw)
1376            CARAC(1)=CARAC(1)+CARAC(4)
1377            CARAC(4)=CARAC(1)
1378            CARAC(2)=CARAC(2)+CARAC(3)
1379            CARAC(3)=CARAC(2)
1380            CARAC(5)=CARAC(5)+CARAC(6)
1381            CARAC(6)=CARAC(5)
1382 !
1383 ! Peso de las componentes de las que se extraeran los autovalores.
1384 !
1385            IF (CARAC(1).GE.CARAC_MAX) THEN ! HH
1386                    ROOT(4)=W(I)
1387            END IF
1388            IF (CARAC(2).GE.CARAC_MAX) THEN ! LH
1389                    ROOT(3)=W(I)
1390            END IF
1391            IF (CARAC(5).GE.CARAC_MAX) THEN ! SO
1392                    ROOT(1)=W(I)
1393            END IF
1394
1395           END DO
1396
1397          ELSE ! Only the diagonal elements were calculated
1398
1399           ROOT(4)=HKANE(1,1)
1400           ROOT(3)=HKANE(2,2)
1401           ROOT(1)=HKANE(5,5)
1402           ROOT(2)=HKANE(6,6) !Redundant
1403
1404
1405          END IF
1406
1407 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1408 !! RESULTS
1409 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1410
1411          ENERGY(1) = POTE(CHI)+ZBC1(CHI)*SHID
1412
1413          ENERGY(2) = ROOT(4)
1414          ENERGY(3) = ROOT(3)
1415          ENERGY(4) = ROOT(3)
1416          ENERGY(5) = ROOT(4)
1417
1418          ENERGY(6) = ROOT(1)
1419          ENERGY(7) = ROOT(1)
1420 !        ESODW(IR,IZ,IGEO) = CHI ! To save the Structure-profile
1421
1422 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1423 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1424 !!!!! Calculation of the Elastic Energy. Since EEL contains no information
1425 !!!!! we will use this array to pack the Elastic Energy.
1426 !!!!!
1427 !!!!!        U=1/2*XLAMB*Tr(e)**2+XMU*(err**2+e00**2+ezz**2+erz**2)
1428 !!!!!
1429 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1430 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1431
1432           ENERGY(8) = XMU*(SRR**2 + S00**2 + SZZ**2 + SRZ**2) +    &
1433                       XLAMB/2.E0*(SSUM + SZZ)**2
1434
1435           RETURN
1436
1437       END SUBROUTINE POT_CALC_CYL_ZB
1438
1439       END SUBROUTINE POTENTIAL_ZB
1440