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