pot_slat.f90

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