      PROGRAM MoS2_continuation
c
c     length unit is Bohr
c     energy unit is either eV or Hartee (so, be careful! )
c
c
c     P as an external, extrapolated input 
c
c     Included the semimetal regime (modified BSE for plasmons), 
c     Thomas-Fermi screening,
c     automated linear regressione extrapolation, evaluation density of
c     states at the Fermi energy.
c     The energy zero is placed at midgap; the indirect bandgap is E_x;
c     in the semimetal E_x < 0 . 
c
c     WARNING: valence band effective masses are extrapolated as
c     positive quantities but used as negative quantities
c
c     Massimo Rontani 23 May 2019
c
c     Added energy cutoff to truncate the Hilbert space 
c     for the semimetal calculation
c
c     Massimo Rontani 10 June 2019
c
c     Using the renormalized vertex evaluated by A. N. Kozlov, L. A. Maksimov 
c     [The metal-dielectric divalent crystal phase transition,
c     JETP 21, 790-795 (1965)] for the two-band, isotropic-mass model
c     with Delta << epsilon_F on the semimetal side of the phase diagram
c
c     Massimo Rontani 26 June 2019
c
c     Freely available version
c
c     Massimo Rontani 28 October 2020
c
      IMPLICIT NONE
      REAL ee
      PARAMETER( ee = 27.210 )  ! ee = 27.21 eV; Hartree energy
      REAL b1(3), b2(3), b3(3)
      INTEGER N_kpts, i,  ii, j
      REAL k_Lambda_mod
      REAL E_qp  
      REAL i_point, i_band, nothing 
      REAL Delta_qx, Delta_qy, Delta_qz, k_x, k_y, k_z, Re_W, Im_W
      REAL Delta_q, Wq0, pi, kappa_r
      REAL eps, kappino, delta_kappino
      INTEGER N_kappino
      PARAMETER( delta_kappino = 0.001, N_kappino = 150 )
      REAL y_lm, y_l, y_lp, x_lm, x_l, x_lp, x_min, y_min
      REAL m_v_K, m_v_perp, m_v_z
      REAL m_c_K, m_c_perp, m_c_z
      REAL mu , mu_K , mu_perp , mu_z, E_x, E_b
      CHARACTER*40 word1, word2, word3
c
c     Effective-mass computation 
c
      REAL Wq0_EM, Delta_qx_EM, Delta_qy_EM, Delta_qz_EM
      REAL Dqx, Dqy, Dqz
      REAL theta_1, theta_2, phi_1, sum_a, sum_b, delta_a, delta_b
      REAL phi
      INTEGER i_phi, nint_phi
      PARAMETER( nint_phi = 1000 )
      REAL I_0, I_1, I_2, I_3, I_4
c     cutoff_i (i = x, y, z) is the cube side sampled 
c     in the BSE equation (in units of Bohr^1) and N_i is
c     the number of points (actually it is 2 * N_i/2 + 1 )
      REAL cutoff_x, cutoff_y, cutoff_z
      INTEGER N_x, N_y, N_z
CCC      PARAMETER( Dqx = 0.12, Dqy = 0.121 , Dqz = 2. * 0.04136 )
c      PARAMETER( N_x = 26, N_y = 26, N_z = 30 )
      PARAMETER( N_x = 13, N_y = 13, N_z = 13 )
      PARAMETER(cutoff_x= .05, cutoff_y= .05, cutoff_z= 0.06 )
CCC      PARAMETER(cutoff_x= .015, cutoff_y= .015, cutoff_z= 0.015 )
      PARAMETER( Dqx = cutoff_x / ( ( N_x /2 )*2 ) )
      PARAMETER( Dqy = cutoff_y / ( ( N_y/2)*2 ) )
      PARAMETER( Dqz = cutoff_z / ( ( N_z/2)*2 ) )
      INTEGER k, l, lpr, level, lstore, ipr, jpr, kpr, i_nn_sqrt
      INTEGER level0, lmax, lmaxx, lmaxy, lmaxz
      PARAMETER( level0 = ( (N_x/2)*2 + 1) * ( (N_y/2)*2 + 1 ) * 
     &                     ( (N_z/2)*2 + 1 ) )
      INTEGER, ALLOCATABLE :: iindex(:), jindex(:), kindex(:)
      REAL E_eh_cutoff, sqrwf_cutoff
      PARAMETER( E_eh_cutoff = 0.065 / ee )
c      PARAMETER( E_eh_cutoff = 0.015 / ee )
      PARAMETER( sqrwf_cutoff = 0.99 )
CCC      INTEGER :: iindex( 1:level ), jindex( 1:level ), kindex( 1:level)
      REAL, ALLOCATABLE :: E_cv( : , : , : )
      REAL, ALLOCATABLE :: EEpsi_c( : , : , : ), EEpsi_v( : , : , : )
      REAL qsq
c     routine DSYEV (Lapack) variables
      INTEGER levelmax, ifail, lwork, itype
      PARAMETER( levelmax = 16500, lwork = 64 * levelmax )
      CHARACTER*1 job , uplo
      REAL*8 work(1:lwork)
      REAL*8 H(1:levelmax,1:levelmax)
      real*8 eigen(1:levelmax)
      real*8 maxpsisqr
c     mean-field variables
      REAL, ALLOCATABLE :: Delta_old(:,:,:), Delta_new(:,:,:)
      REAL, ALLOCATABLE :: phi_sum(:,:,:), Warray(:,:,:)
      REAL epsilon_u, tolerance, maxdiff, psi_u_max, summ
      INTEGER elle, n_iteration
      PARAMETER( n_iteration = 500 )
      PARAMETER( tolerance = 1.e-6 )
      REAL T_min, T_max, deltaT, beta, kBT, Tfact, T
      REAL E_vtop, E_cbottom, epsi_c, epsi_v
      INTEGER iT, nintT
c     extrapolation part
      REAL Eg(2), alpha, x1, x2, P, P0
      PARAMETER( x1 = 25. , x2 = 34. )
      REAL kappa_rx(2), m_cKx(2), m_cperpx(2), m_czx(2)
      REAL m_vKx(2), m_vperpx(2), m_vzx(2)
      INTEGER unit_kpts 
      CHARACTER*12 file_kpts 
c     end extrapolation part
c     semimetal part
      REAL epsilon_F, E_xprime
      REAL delta_e, delta_target, energy, sum_c, sum_v, sigma
      REAL window_e, DOS_F, fFermi_v, fFermi_c
      PARAMETER( delta_target = .0001 / ee ) ! energy resolution, in Hartree
      INTEGER nint_e, ie, ie_min, signum, n_a, n_b
      REAL, ALLOCATABLE :: density_c(:), density_v(:) 
      REAL, ALLOCATABLE :: I_density_c(:), I_density_v(:)
      REAL, ALLOCATABLE :: epsi_c0( : , : , : ), epsi_v0( : , : , : ) 
      REAL factor_Maksimov, G_0_abs, Delta_G_0_abs, G, alpha_Maksimov
      PARAMETER( G_0_abs = 0.00938/ee , Delta_G_0_abs = 0.008/ee )
      REAL Delta_trial, factor_Maksimov_old, factor_Maksimov_new
      INTEGER i_F, j_F, k_F
      REAL, ALLOCATABLE :: temparray(:)
      REAL Delta_k_F_av, EEpsi_cmin, EEpsi_vmax
      INTEGER i_min(3), i_max(3)
      REAL dipole_z, P_z, a_Bohr, eA0	
      PARAMETER( a_Bohr = 0.5292 ) !! Bohr in Angstrom units
      PARAMETER( dipole_z = 3.78 ) !! units of e Bohr 
      PARAMETER( eA0 = 4.802 ) !! e A0 in Debye units
c
      EXTERNAL DYSEV
      
      IF( Dqx .le. 0. ) STOP 'parameter error'
      IF( Dqy .le. 0. ) STOP 'parameter error'
      IF( Dqz .le. 0. ) STOP 'parameter error'
      IF( cutoff_x .le. 0. ) STOP 'parameter error'
      IF( cutoff_y .le. 0. ) STOP 'parameter error'
      IF( cutoff_z .le. 0. ) STOP 'parameter error'

      pi = 4. * atan( 1.)
c
      nintT = 1
      T_min = 0.0d0  ! Kelvin
      T_max = 80.0d0 ! Kelvin
      deltaT = ( T_max - T_min ) / max( 1, nintT - 1)
c     P = 34.0457611
c     P = 34.22
      P = 34.0
c
      write(*,*)
      write(*,*) '~~~~~~~~~~~~~~MoS2 vs P~~~~~~~~~~~~~~~'
      write(*,*)
      write(*,*)'min value of T (K) =',T_min
      write(*,*)'max value of T (K) =',T_max
      write(*,*)'number of T-points =',nintT
      write(*,*)
      write(*,*) 'Input pressure =',P,'(GPa)'
      write(*,*)
      write(*,*) 'Gamma-Lambda Dipole computed by Daniele  ='
      write(*,*) dipole_z,'(units of electric charge times Bohr)'
      write(*,*) dipole_z * a_Bohr,'(units of electric charge times A0)'
      write(*,*) dipole_z * a_Bohr * eA0 ,'(units of Debye)'
      write(*,*)

            
c
c     reciprocal lattice vectors at P = 25 GPa
c

      b1(1) = 1.10085297
      b1(2) = 0.635577977
      b1(3) = 0.000000

      b2(1) = 1.51099997E-07
      b2(2) = 1.27115595
      b2(3) = 0.000000

      b3(1) = 0.000000
      b3(2) = 0.000000
      b3(3) = 0.302446008 

c
c     reciprocal lattice vectors at P = 34 GPa
c
      b1(1) = 1.116741 
      b1(2) = 0.644751
      b1(3) = 0.000000

      b2(1) = 0.3547E-7
      b2(2) = 1.289502
      b2(3) = 0.000000

      b3(1) = 0.000000
      b3(2) = 0.000000
      b3(3) = 0.307613

      write(*,*)
      write(*,*) 'Length unit = Bohr; Coordinates are cartesian'
      write(*,*)
      write(*,*) 'Reciprocal lattice vectors at 34 GPa'
      write(*,*) 'b1 =', b1
      write(*,*) 'b2 =', b2 
      write(*,*) 'b3 =', b3
      write(*,*)
      write(*,*)
      write(*,*)

c
c     extrapolation
c
      Eg(1) = 0.379127324
      Eg(2) = 0.00191795
      kappa_rx(1) = 16.398
      kappa_rx(2) = 16.782
      m_cKx(1) = 0.467051864
      m_cKx(2) = 0.396498442
      m_cperpx(1) = 0.383014709
      m_cperpx(2) = 0.318608373
      m_czx(1) = 0.395981252
      m_czx(2) = 0.340174794
      m_vKx(1) = 0.297190428
      m_vKx(2) = 0.235990331
      m_vperpx(1) = 0.283419490
      m_vperpx(2) = 0.237984970
      m_vzx(1) = 1.36233139
      m_vzx(2) = 1.54141688
c
      alpha = (Eg(2) - Eg(1) )/(x2 - x1)
      P0 = -Eg(1)/alpha + x1
      E_x = alpha*(P-x1) + Eg(1)
      E_xprime = E_x / ee
      E_cbottom = E_xprime * .5
      E_vtop = - E_xprime * .5

      write(*,*)
      write(*,*) 'Zero-gap pressure = ', P0,'(GPa)'
      write(*,*) 'Indirect gap =',E_x,'(eV)'
      write(*,*)

      alpha = (kappa_rx(2) - kappa_rx(1))/(x2 - x1)
      kappa_r = alpha*(P-x1) + kappa_rx(1)

      write(*,*) 'Zero-gap kappa_r = ',
     &           alpha*(P0-x1) + kappa_rx(1) , '(dimensionless)'
      write(*,*) 'Kappa_r =', kappa_r, '(dimensionless)'
      write(*,*)

      alpha = (m_cKx(2) - m_cKx(1))/(x2 - x1)
      m_c_K = alpha*(P-x1) + m_cKx(1)

      write(*,*) 'Zero-gap m_cK = ',
     &           alpha*(P0-x1) + m_cKx(1) , '(units of m_e)'
      write(*,*) 'm_cK =', m_c_K, '(units of m_e)'
      write(*,*)

      alpha = (m_cperpx(2) - m_cperpx(1))/(x2 - x1)
      m_c_perp = alpha*(P-x1) + m_cperpx(1)

      write(*,*) 'Zero-gap m_cperp = ',
     &           alpha*(P0-x1) + m_cperpx(1) , '(units of m_e)'
      write(*,*) 'm_cperp =', m_c_perp, '(units of m_e)'
      write(*,*)

      alpha = (m_czx(2) - m_czx(1))/(x2 - x1)
      m_c_z = alpha*(P-x1) + m_czx(1)

      write(*,*) 'Zero-gap m_cz = ',
     &           alpha*(P0-x1) + m_czx(1) , '(units of m_e)'
      write(*,*) 'm_cz =', m_c_z, '(units of m_e)'
      write(*,*)

      alpha = (m_vKx(2) - m_vKx(1))/(x2 - x1)
      m_v_K = alpha*(P-x1) + m_vKx(1)

      write(*,*) 'Zero-gap m_vK = ',
     &           alpha*(P0-x1) + m_vKx(1) , '(units of m_e)'
      write(*,*) 'm_vK =', m_v_K, '(units of m_e)'
      write(*,*)

      alpha = (m_vperpx(2) - m_vperpx(1))/(x2 - x1)
      m_v_perp = alpha*(P-x1) + m_vperpx(1)

      write(*,*) 'Zero-gap m_vperp = ',
     &           alpha*(P0-x1) + m_vperpx(1) , '(units of m_e)'
      write(*,*) 'm_vperp =', m_v_perp, '(units of m_e)'
      write(*,*)

      alpha = (m_vzx(2) - m_vzx(1))/(x2 - x1)
      m_v_z = alpha*(P-x1) + m_vzx(1)

      write(*,*) 'Zero-gap m_vz = ',
     &           alpha*(P0-x1) + m_vzx(1) , '(units of m_e)'
      write(*,*) 'm_vz =', m_v_z, '(units of m_e)'
      write(*,*)

      m_v_K = - m_v_K
      m_v_perp = - m_v_perp 
      m_v_z = - m_v_z

      if( m_c_K .lt. 0. ) stop 'wrong sign of m_c_K'
      if( m_c_perp .lt. 0. ) stop 'wrong sign of m_c_perp'
      if( m_c_z .lt. 0. ) stop 'wrong sign of m_c_z'
      if( m_v_K .gt. 0. ) stop 'wrong sign of m_v_K'
      if( m_v_perp .gt. 0. ) stop 'wrong sign of m_v_perp'
      if( m_v_z .gt. 0. ) stop 'wrong sign of m_v_z'

      write(*,*)
      write(*,*) 'Pressure =',P,' GPa'
      write(*,*) 
      write(*,*) 'Effective masses at Lambda (cb) and Gamma (vb) pts'
      write(*,*) '(in units of m_e )'
      write(*,*)
      write(*,*) 'along line Gamma - Kappa'
      write(*,*) ' cb =', m_c_K,'vb =', abs(m_v_K)
      write(*,*)
      write(*,*) 'along line perp to Gamma - Kappa'
      write(*,*) ' cb =', m_c_perp ,'vb =', abs(m_v_perp)
      write(*,*)
      write(*,*) 'along z'
      write(*,*) ' cb =', m_c_z ,'vb =', abs(m_v_z)
      write(*,*)

      mu_K = 1./(1./m_c_K + 1./abs(m_v_K) )
      mu_perp = 1./(1./m_c_perp + 1./abs(m_v_perp))
      mu_z = 1./(1./m_c_z + 1./abs(m_v_z) )
      mu = ( mu_K * mu_perp * mu_z )**(0.333333)

      write(*,*) 'Exciton masses'
      write(*,*) ' mu_K =',mu_K
      write(*,*) ' mu_perp =',mu_perp
      write(*,*) ' mu_z =',mu_z
      write(*,*) 'algebraic mean =',mu

ccc      Delta_qx = k_IBZ(1,29)
ccc      Delta_qy = k_IBZ(2,3)
ccc      Delta_qz = k_IBZ(3,2)
ccc      Delta_q = min( Delta_qx, Delta_qy, Delta_qz )
ccc      write(*,*) 'Delta_qx =',Delta_qx,'(Bohr^-1)'
ccc      write(*,*) 'Delta_qy =',Delta_qy,'(Bohr^-1)'
ccc      write(*,*) 'Delta_qz =',Delta_qz,'(Bohr^-1)'
ccc      write(*,*) 'Delta_q =',Delta_q,'(Bohr^-1)'

      write(*,*)
      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      write(*,*) 'Effective mass calculation'
      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      write(*,*)

      write(*,*) 'Long-wavelength dielectric constant =',kappa_r
      write(*,*) ' N_x =',N_x,'N_y =',N_y,'N_z =',N_z
      write(*,*) ' cutoff_x =',cutoff_x,'cutoff_y =',cutoff_y,
     &           'cutoff_z =',cutoff_z,'(units of Bohr)'
      write(*,*) ' Total number of q points without cutoff =',level0
      write(*,*) ' Dqx = ',Dqx,'Dqy =',Dqy,'Dqz =',Dqz,'(Bohr^-1)'

      ALLOCATE( epsi_c0( -N_x/2:N_x/2 , -N_y/2:N_y/2 , -N_z/2:N_z/2) )
      ALLOCATE( epsi_v0( -N_x/2:N_x/2 , -N_y/2:N_y/2 , -N_z/2:N_z/2) )
c               epsi_c0, epsi_v0 are in units of Hartree

      do i = -N_x/2, N_x/2   ! k_x loop
       do j = -N_y/2, N_y/2   ! k_y loop
        do k = -N_z/2, N_z/2   ! k_z loop

         epsi_c0( i, j, k) = E_cbottom + (i*Dqx)**2 / 2. / m_c_K +
     &   (j*Dqy)**2 /  2. / m_c_perp + (k*Dqz)**2 /  2. / m_c_z
         epsi_v0( i, j, k) = E_vtop + (i*Dqx)**2 / 2. / m_v_K +
     &   (j*Dqy)**2 /  2. / m_v_perp + (k*Dqz)**2 /  2. / m_v_z

        enddo
       enddo
      enddo

c
c     semiconductor / semimetal discrimination
c

      IF( E_xprime .ge. 0.0d0 )THEN ! exciton

       write(*,*)
       write(*,*) 'PRISTINE SEMICONDUCTOR SYSTEM'
       write(*,*)
       write(*,*) '3d exciton Bohr radius =',kappa_r/mu,'(Bohr)'
       write(*,*)
       write(*,*) '3d exciton binding energy, with average mass =',
     &             ee / 2. / kappa_r**2 * mu,'(eV)'
       E_b = 4. * ee / 2. / kappa_r**2 * sqrt( mu_K * mu_perp)
       write(*,*) '2d exciton binding energy, with average mass =',
     &             E_b,'(eV)'
       write(*,*)
       write(*,*) '2d exciton excitation energy =',E_x-E_b,'(eV)'
       write(*,*)
       write(*,*)

       epsilon_F = 0.
       write(*,*) 'Fermi energy =',epsilon_F * ee,'(eV)'

      ELSE ! plasmon

       write(*,*)
       write(*,*) 'PRISTINE SEMIMETAL SYSTEM'
       write(*,*)

       window_e = 10. * abs( E_xprime )
       nint_e = window_e / delta_target
       delta_e = window_e / nint_e
       nint_e = nint_e + 1
       sigma = 20. * delta_target
       write(*,*) 'Energy window =',window_e * ee,'(eV)'
       write(*,*) 'Energy resolution = ', delta_e * ee,'(eV)'
       write(*,*) 'sigma = ',sigma * ee,'(eV)'
       write(*,*) 'no. of energy gripoints =', nint_e
       ALLOCATE( density_c( 1:nint_e) , density_v( 1:nint_e) )
       ALLOCATE( I_density_c( 1:nint_e) , I_density_v( 1:nint_e) )

       do ie = 1, nint_e

        energy = - window_e * .5 + delta_e * (ie - 1) 

        sum_c = 0.
        sum_v = 0.
        do i = -N_x/2, N_x/2   ! k_x loop
         do j = -N_y/2, N_y/2   ! k_y loop
          do k = -N_z/2, N_z/2   ! k_z loop
           sum_c = sum_c + exp( -.5 * ( energy - epsi_c0(i,j,k) )**2 /
     &             sigma**2 ) / sigma / sqrt( 2. * pi )
     &             * 12.  
c          spin and valley degeneracy, normalized to level
           sum_v = sum_v + exp( -.5 * ( energy - epsi_v0(i,j,k) )**2 /
     &             sigma**2 ) / sigma / sqrt( 2. * pi )
     &             * 2.   
c          spin degeneracy, normalized to level
c          here I use a Gaussian with sigma as a multiple of delta_target
          enddo
         enddo
        enddo

        density_c( ie ) = sum_c * Dqx * Dqy * Dqz / (2. * pi)**3
        density_v( ie ) = sum_v * Dqx * Dqy * Dqz / (2. * pi)**3

        write(9,*) energy * ee, density_c( ie)/ee, density_v( ie)/ee,
     &  m_c_K**1.5 * sqrt( abs( energy - E_xprime * 0.5 ) ) * sqrt(2.)
     &  / pi**2 * 6. / ee, 
     &  m_c_K**1.5 * sqrt( abs( energy + E_xprime * 0.5 ) ) * sqrt(2.)
     &  / pi**2 / ee 
        write(11,*) energy * ee, (density_c( ie) + density_v( ie) )/ee

       enddo

       do ie = 1, nint_e
        I_density_c( ie ) = SUM( density_c( 1 : ie ) ) * delta_e
        I_density_v( ie ) = SUM( density_v( ie : nint_e ) ) * delta_e
        energy = - window_e * .5 + delta_e * (ie - 1) 
        write(10,*) energy*ee, I_density_c( ie), I_density_v( ie)
       enddo

       ie_min = MINLOC( abs( I_density_c - I_density_v ) , 1 )
       epsilon_F = - window_e * .5 + delta_e * ( ie_min - 1 )
       DOS_F = density_c( ie_min ) + density_v( ie_min )

       if( ( epsilon_F .lt. .5 * E_xprime ) .or.
     &     ( epsilon_F .gt. -.5 * E_xprime ) )then
            write(*,*) 'epsilon_F out of range'
            epsilon_F = -0.54 * abs( E_xprime ) * .5
            DOS_F =
     &       ( m_c_K**1.5 * sqrt( abs( epsilon_F - E_xprime * 0.5 ) ) 
     &        * sqrt(2.) / pi**2 * 6. ) + (
     &          m_c_K**1.5 * sqrt( abs( epsilon_F + E_xprime * 0.5 ) ) 
     &           * sqrt(2.) / pi**2 ) 
       endif

       write(*,*) 'Fermi energy =',epsilon_F * ee,'(eV)'
       write(*,*) 'DOS at Fermi energy =',DOS_F,' (Hartree^-1 Bohr^-3)'
       write(*,*) 'Total number of single-species carriers'
       write(*,*) '(either electrons or holes) present in the'
       write(*,*) 'noninteracting semimetal per volume ='
       write(*,*) I_density_c( ie_min ) ,'(Bohr^-3)'
       write(*,*)
       write(*,*) 'Interaction range in q space =',
     &  sqrt( 4. * pi * DOS_F / kappa_r ),'(units of Bohr^-1)'

       DEALLOCATE( density_c , density_v )
       DEALLOCATE( I_density_c , I_density_v )

       ALLOCATE( temparray( 0 : N_x/2 ) )
       do i = 0, N_x/2
        temparray(i) = abs( epsi_c0(i,0,0) - epsilon_F )
       enddo
       i_F = MINLOC( temparray, 1 )
       write(*,*) 'k_x_Fermi =', i_F * Dqx ,'(Bohr^-1)'
       DEALLOCATE( temparray )

       ALLOCATE( temparray( 0 : N_y/2 ) )
       do j = 0, N_y/2
        temparray(j) = abs( epsi_c0(0,j,0) - epsilon_F )
       enddo
       j_F = MINLOC( temparray, 1 )
       write(*,*) 'k_y_Fermi =', j_F * Dqy ,'(Bohr^-1)'
       DEALLOCATE( temparray )

       ALLOCATE( temparray( 0 : N_z/2 ) )
       do k = 0, N_z/2
        temparray(k) = abs( epsi_c0(0,0,k) - epsilon_F )
       enddo
       k_F = MINLOC( temparray, 1 )
       write(*,*) 'k_z_Fermi =', k_F * Dqz ,'(Bohr^-1)'
       DEALLOCATE( temparray )


      ENDIF ! end discrimination semiconductor / semimetal

      ALLOCATE( E_cv( -N_x/2: N_x/2 , -N_y/2: N_y/2 , -N_z/2: N_z/2) ) 
c               E_cv is in units of Hartree

c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
c     Hilbert space truncation by cutting off eh pair excitation spectrum
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      l = 0
      do i = -N_x/2, N_x/2   ! k_x loop
       do j = -N_y/2, N_y/2   ! k_y loop
        do k = -N_z/2, N_z/2   ! k_z loop

         E_cv( i, j, k) = E_xprime  + (i*Dqx)**2 / ( 2. * mu_K ) +
     &   (j*Dqy)**2 / ( 2. * mu_perp ) + (k*Dqz)**2 / ( 2. * mu_z )
c        it is in Hartree units !

         if( ( abs( E_cv( i, j, k) ) .le. E_eh_cutoff ) )then
          l = l + 1
         endif

        enddo
       enddo
      enddo

      level = l

      write(*,*)
      write(*,*) 'e-h excitation energy cutoff =',E_eh_cutoff*ee,'eV'
      write(*,*) 'level =', level
      write(*,*)


      ALLOCATE( iindex(1:level), jindex(1:level), kindex(1:level) )

      l = 0
      do i = -N_x/2, N_x/2   ! k_x loop
       do j = -N_y/2, N_y/2   ! k_y loop
        do k = -N_z/2, N_z/2   ! k_z loop

         if( ( abs( E_cv( i, j, k) ) .le. E_eh_cutoff ) )then
          l = l + 1
          iindex(l) = i
          jindex(l) = j
          kindex(l) = k
         endif

        enddo
       enddo
      enddo

c     some useful output

      do l = 1, level

       if( (jindex(l).eq.0).and.(kindex(l).eq.0) )then
        write(43,*) iindex(l)*Dqx, 
     &  E_cv( iindex(l), jindex(l), kindex(l) ) * ee,
     &  abs( E_cv( iindex(l), jindex(l), kindex(l) ) * ee )  
       endif
       if( (iindex(l).eq.0).and.(kindex(l).eq.0) )then
        write(44,*) jindex(l)*Dqy, 
     &  E_cv( iindex(l), jindex(l), kindex(l) ) * ee,
     &  abs( E_cv( iindex(l), jindex(l), kindex(l) ) * ee )  
       endif
       if( (iindex(l).eq.0).and.(jindex(l).eq.0) )then
        write(45,*) kindex(l)*Dqz, 
     &  E_cv( iindex(l), jindex(l), kindex(l) ) * ee,
     &  abs( E_cv( iindex(l), jindex(l), kindex(l) ) * ee )  
       endif

      enddo

      do i = -N_x/2, N_x/2 
       write(40,*) i*Dqx, E_cv( i, 0, 0) * ee, 
     &                     abs( E_cv( i, 0, 0) * ee) 
      enddo
      do j = -N_y/2, N_y/2 
       write(41,*) j*Dqy, E_cv( 0, j, 0) * ee,
     &                     abs( E_cv( 0, j, 0) * ee) 
      enddo
      do k = -N_z/2, N_z/2 
       write(42,*) k*Dqz, E_cv( 0, 0, k) * ee,
     &                     abs( E_cv( 0, 0, k) * ee) 
      enddo

c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      write(*,*)
      write(*,*) ' Calculation of W( q = 0 ) for the gapped system'
      write(*,*)
c
c     The Coulomb divergence is integrated over a parallelepipid
c     centered at the origin, whose size is 
c     ( -Dqx : Dqx, -Dqy : Dqy, -Dqz : Dqz ).
c
      if( (Dqz .le. Dqx).and.(Dqz .le. Dqy) )then
       Delta_qx_EM = Dqx
       Delta_qy_EM = Dqy
       Delta_qz_EM = Dqz
      elseif( (Dqx .le. Dqy).and.(Dqx .le. Dqz) )then
       Delta_qx_EM = Dqy
       Delta_qy_EM = Dqz
       Delta_qz_EM = Dqx
      elseif( (Dqy .le. Dqx).and.(Dqy .le. Dqz) )then
       Delta_qx_EM = Dqz
       Delta_qy_EM = Dqx
       Delta_qz_EM = Dqy
      endif

      write(*,*) ' Delta_qx_EM =',Delta_qx_EM,'(Bohr^-1)'
      write(*,*) ' Delta_qy_EM =',Delta_qy_EM,'(Bohr^-1)'
      write(*,*) ' Delta_qz_EM =',Delta_qz_EM,'(Bohr^-1)'
      write(*,*)

      Delta_qx_EM = 2. * Delta_qx_EM ! We redifine here the mesh size since,
c                                      in the previous implementation,
c                                      the integration range was half
c                                      the value
      Delta_qy_EM = 2. * Delta_qy_EM
      Delta_qz_EM = 2. * Delta_qz_EM

      phi_1 = atan( Delta_qy_EM / Delta_qx_EM ) 

c      write(*,*) ' phi_1 = ',phi_1 / pi , '( rad / pi )'

      I_0 = pi / 4. * Delta_qz_EM

      write( *, * ) 'I_0 =',I_0,'( Bohr^-1 )'

      Wq0_EM = 4. / pi**2 / kappa_r * I_0 

      write( *, * ) ' W(q=0) (only I_0) = ', Wq0_EM,'(Hartree)'

c     I_1 ... I_4


c     I_1

      delta_a = phi_1 / nint_phi
      sum_a = 0.

c     i = 0
      phi = 0.
      theta_1 = atan( Delta_qx_EM / Delta_qz_EM / cos( phi ) )
      sum_a = sum_a + .5 * delta_a * ( 
     &         Delta_qx_EM/2. / cos( phi ) * ( pi/2. - theta_1 )
     &       - Delta_qz_EM/2. * cos( theta_1 ) )

      do i = 1, nint_phi - 1
       phi = delta_a *  i 
       theta_1 = atan( Delta_qx_EM / Delta_qz_EM / cos( phi ) )
       sum_a = sum_a + delta_a * ( 
     &         Delta_qx_EM/2. / cos( phi ) * ( pi/2. - theta_1 )
     &       - Delta_qz_EM/2. * cos( theta_1 ) )
c       write(17,*) phi, Delta_qx_EM/2./cos( phi )*( pi/2. - theta_1)
c     &             - Delta_qz_EM/2. * cos( theta_1 ) 
      enddo

      phi = phi_1
      theta_1 = atan( Delta_qx_EM / Delta_qz_EM / cos( phi ) )
      sum_a = sum_a + .5 * delta_a * ( 
     &         Delta_qx_EM/2. / cos( phi ) * ( pi/2. - theta_1 )
     &       - Delta_qz_EM/2. * cos( theta_1 ) )

      I_1 = sum_a
      write( *, * ) 'I_1 = ',I_1,'(Hartree)'

c     I_2

      delta_b = ( pi/2. - phi_1 ) / nint_phi
      sum_b = 0.

c     i = 0
      phi = phi_1
      theta_2 = atan( Delta_qy_EM / Delta_qz_EM / sin( phi ) )
      sum_b = sum_b + .5 * delta_b * ( 
     &         Delta_qy_EM/2. / sin( phi ) * ( pi/2. - theta_2 )
     &       - Delta_qz_EM/2. * cos( theta_2 ) )

      do i = 1, nint_phi - 1
       phi = delta_b *  i + phi_1
       theta_2 = atan( Delta_qy_EM / Delta_qz_EM / sin( phi ) )
       sum_b = sum_b + delta_b * ( 
     &         Delta_qy_EM/2. / sin( phi ) * ( pi/2. - theta_2 )
     &       - Delta_qz_EM/2. * cos( theta_2 ) )
c       write(17,*) phi, Delta_qy_EM/2./sin( phi )*( pi/2. - theta_2)
c     &             - Delta_qz_EM/2. * cos( theta_2 ) 
      enddo

      phi = pi/2.
      theta_2 = atan( Delta_qy_EM / Delta_qz_EM / sin( phi ) )
      sum_b = sum_b + .5 * delta_b * ( 
     &         Delta_qy_EM/2. / sin( phi ) * ( pi/2. - theta_2 )
     &       - Delta_qz_EM/2. * cos( theta_2 ) )

      I_2 = sum_b
      write( *, * ) 'I_2 = ',I_2,'(Hartree)'

c     I_3

      delta_a = phi_1 / nint_phi
      sum_a = 0.

c     i = 0
      phi = 0.
      theta_1 = atan( Delta_qx_EM / Delta_qz_EM / cos( phi ) )
      sum_a = sum_a + .5 * delta_a * ( 
     &         Delta_qz_EM/2. / cos( phi ) * ( -log(cos(theta_1)) )
     &       - Delta_qz_EM/2. * (1. - cos(theta_1)) )

      do i = 1, nint_phi - 1
       phi = delta_a *  i 
       theta_1 = atan( Delta_qx_EM / Delta_qz_EM / cos( phi ) )
       sum_a = sum_a + delta_a * ( 
     &         Delta_qz_EM/2. / cos( phi ) * ( -log(cos(theta_1)) )
     &       - Delta_qz_EM/2. * (1.- cos(theta_1)) )
c       write(18,*) phi, Delta_qz_EM/2./cos(phi)*( -log(cos(theta_1)) )
c     &             - Delta_qz_EM/2. * (1. - cos(theta_1) )
      enddo

      phi = phi_1
      theta_1 = atan( Delta_qx_EM / Delta_qz_EM / cos( phi ) )
      sum_a = sum_a + .5 * delta_a * ( 
     &         Delta_qx_EM/2. / cos( phi ) * ( -log(cos(theta_1)) )
     &       - Delta_qz_EM/2. * (1.- cos(theta_1)) )

      I_3 = sum_a
      write( *, * ) 'I_3 = ',I_3,'(Hartree)'

c     I_4

      delta_b = ( pi/2. - phi_1 ) / nint_phi
      sum_b = 0.

c     i = 0
      phi = phi_1
      theta_2 = atan( Delta_qy_EM / Delta_qz_EM / sin( phi ) )
      sum_b = sum_b + .5 * delta_b * ( 
     &         Delta_qz_EM/2. / sin( phi ) * ( - log(cos(theta_2)) )
     &       - Delta_qz_EM/2. * (1. - cos(theta_2)) )

      do i = 1, nint_phi - 1
       phi = delta_b *  i + phi_1
       theta_2 = atan( Delta_qy_EM / Delta_qz_EM / sin( phi ) )
       sum_b = sum_b + delta_b * ( 
     &         Delta_qz_EM/2. / sin( phi ) * (  -log(cos(theta_2)) )
     &       - Delta_qz_EM/2. * (1. - cos(theta_2)) )
c       write(18,*) phi, Delta_qz_EM/2./sin(phi)*( -log(cos(theta_2)) )
c     &             - Delta_qz_EM/2. * (1. - cos(theta_2) ) 
      enddo

      phi = pi/2.
      theta_2 = atan( Delta_qy_EM / Delta_qz_EM / sin( phi ) )
      sum_b = sum_b + .5 * delta_b * ( 
     &         Delta_qz_EM/2. / sin( phi ) * ( -log(cos(theta_2)) )
     &       - Delta_qz_EM/2. * (1. - cos(theta_2)) )

      I_4 = sum_b
      write( *, * ) 'I_4 = ',I_4,'(Hartree)'

      Wq0_EM = 4. / pi**2 / kappa_r * ( I_0 + I_1 + I_2 + I_3 + I_4 )

      write( *, * ) ' W(q=0) = ', Wq0_EM,'(Hartree)'
      write(*,*)

c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c     End calculation of W( q = 0 ) for the gapped system
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


      write(*,*)
      write(*,*) 'Computing the Bethe-Salpeter kernel'
      write(*,*)


      IF( E_xprime .ge. 0. )THEN

       write(*,*)
       write(*,*) 'PRISTINE SEMICONDUCTOR SYSTEM'
       write(*,*)

      H = 0.0d0

      do l = 1, level

       i = iindex(l)
       j = jindex(l)
       k = kindex(l)
       H( l , l ) = E_cv( i, j, k ) - Wq0_EM

       if( (i.eq. 0).and.(j .eq. 0).and.(k .eq. 0) )then
        write(12,*) 0, - Wq0_EM
       endif

       do lpr = 1, l - 1

        ipr = iindex( lpr )
        jpr = jindex( lpr )
        kpr = kindex( lpr )
        qsq = Dqx**2 * (i-ipr)**2 + Dqy**2 * (j-jpr)**2 + Dqz**2 *
     &                 (k-kpr)**2
        H( l, lpr ) = - Dqx * Dqy * Dqz / kappa_r / 2. / pi**2 / qsq
        i_nn_sqrt = (i-ipr)**2 + (j-jpr)**2 + (k-kpr)**2
        if ( i_nn_sqrt .eq. 1 )then
         H( l, lpr ) = H( l, lpr ) * 4. / 8.
        elseif( i_nn_sqrt .eq. 2 )then
         H( l, lpr ) = H( l, lpr ) * 6. / 8.
        elseif( i_nn_sqrt .eq. 3 )then
         H( l, lpr ) = H( l, lpr ) * 7. / 8.
        endif

        if( (i .eq. 0).and.(j .eq. 0).and.(k .eq. 0).and.(jpr .eq. 0)
     &      .and.( kpr .eq. 0 ) )then
          write(12,*) i - ipr,  
     &     - Dqx * Dqy * Dqz / kappa_r / 2. / pi**2 / qsq
        endif

       enddo
      enddo

      ENDIF  ! end semiconductor case

      IF( E_xprime .lt. 0. )THEN

       write(*,*)
       write(*,*) 'PRISTINE SEMIMETAL SYSTEM'
       write(*,*)

      H = 0.0d0

      G = E_xprime
      if( abs(E_xprime) .ge. G_0_abs )then
       factor_Maksimov = 0.
      else
       alpha_Maksimov = sqrt( G_0_abs * (epsilon_F - G * 0.5 )**3 )
       Delta_trial = Delta_G_0_abs * sqrt( 1. - abs( G ) / G_0_abs )
       factor_Maksimov = 1. / ( 1. + alpha_Maksimov / Delta_trial**2 )
      endif
      write(*,*)
      write(*,*) 'sqrt(alpha_Maksimov) =',sqrt(alpha_Maksimov) 
     &            * ee,'(eV)' 
      write(*,*) 'Fermi energy (absolute value) =',
     &            (epsilon_F - G * 0.5)*ee,'(eV)'
      write(*,*) 'Delta tentative =',Delta_trial * ee,'(eV)'
      write(*,*) 'factor_Maksimov =',factor_Maksimov
      write(*,*)

      do l = 1, level

       i = iindex(l)
       j = jindex(l)
       k = kindex(l)
       H( l , l ) = E_cv( i, j, k ) 

       if( epsi_v0( i, j, k ) .le. epsilon_F )then
        n_a = 1
       else
        n_a = 0
       endif

       if( epsi_c0( i, j, k ) .le. epsilon_F )then
        n_b = 1
       else
        n_b = 0
       endif
       signum = n_a - n_b

       H( l , l ) = H( l , l ) - Wq0_EM * signum * factor_Maksimov

       if( (i .eq. 0).and.(j .eq. 0).and.(k .eq. 0) )then
        write(12,*) 0, - Wq0_EM * signum * factor_Maksimov
       endif

       if( factor_Maksimov .eq. 0. )then
        H( l, l ) = H( l, l ) - Dqx * Dqy * Dqz / 2. / pi**2 
     &              / ( 4. * pi * DOS_F ) * signum
        if( (i .eq. 0).and.(j .eq. 0).and.(k .eq. 0) )then
         write(12,*) 0, - Dqx * Dqy * Dqz / 2. / pi**2    
     &              / ( 4. * pi * DOS_F ) * signum
        endif
       endif


       do lpr = 1, l - 1 

        ipr = iindex( lpr )
        jpr = jindex( lpr )
        kpr = kindex( lpr )

        if( epsi_v0( ipr, jpr, kpr ) .le. epsilon_F )then
         n_a = 1
        else
         n_a = 0
        endif

        if( epsi_c0( ipr, jpr, kpr ) .le. epsilon_F )then
         n_b = 1
        else
         n_b = 0
        endif

        signum = n_a - n_b

        qsq = Dqx**2 * (i-ipr)**2 + Dqy**2 * (j-jpr)**2 + Dqz**2 *
     &                 (k-kpr)**2

        H( l, lpr ) = - Dqx * Dqy * Dqz / kappa_r / 2. / pi**2 / qsq
     &                * factor_Maksimov
        i_nn_sqrt = (i-ipr)**2 + (j-jpr)**2 + (k-kpr)**2 * signum
        if ( i_nn_sqrt .eq. 1 )then
         H( l, lpr ) = H( l, lpr ) * 4. / 8.
        elseif( i_nn_sqrt .eq. 2 )then
         H( l, lpr ) = H( l, lpr ) * 6. / 8.
        elseif( i_nn_sqrt .eq. 3 )then
         H( l, lpr ) = H( l, lpr ) * 7. / 8.
        endif

        if( (i .eq. 0).and.(j .eq. 0).and.(k .eq. 0).and.(jpr .eq. 0)
     &      .and.( kpr .eq. 0 ) )then
          write(12,*) i - ipr,  
     &     - Dqx * Dqy * Dqz / kappa_r / 2. / pi**2 / qsq
     &     * factor_Maksimov
        endif

        if( factor_Maksimov .eq. 0. )then
         H( l, lpr ) = H( l, lpr ) - Dqx * Dqy * Dqz / 2. / pi**2 
     &               / ( kappa_r * qsq + 4. * pi * DOS_F ) * signum
         if( (i .eq. 0).and.(j .eq. 0).and.(k .eq. 0).and.(jpr .eq. 0)
     &                 .and.( kpr .eq. 0 ) )then
          write(12,*) i - ipr, - Dqx * Dqy * Dqz / 2. / pi**2 
     &                / ( kappa_r * qsq + 4. * pi * DOS_F ) * signum
         endif
        endif

       enddo
      enddo

      ENDIF  ! end semimetal case

c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c     diagonalizing the BSE kernel
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      ifail = 0
      itype = 1
      job   = 'V'
      uplo  = 'L'

      call DSYEV( job, uplo, level, H,
     +            levelmax, eigen, work, lwork, ifail )

      if(ifail.ne.0)then
       WRITE(*,*) 'Error in routine DSYEV: ifail=',ifail
       stop
      endif

      write(*,*)
      write(*,*) 'The lowest one hundred levels are:'
      do l = 1, min( 100 , level )
       write(*,*) l,sngl(eigen(l)) * ee,'(eV)'
      enddo

      IF( E_xprime .lt. 0. )THEN  ! semimetal case

       do l = 1, level

        maxpsisqr = maxval( H(1:level,l)**2 )
        if( maxpsisqr .lt. sqrwf_cutoff )then
         lstore = l
         lmax = maxloc( H(1:level,l)**2 , 1 )
         lmaxx = iindex( lmax )
         lmaxy = jindex( lmax )
         lmaxz = kindex( lmax )
         write(*,*) 'plasmon wave function found'
         write(*,*) 'maximum sqr wf weight =',maxpsisqr
         write(*,*) 'located at k = (', Dqx*lmaxx ,',', Dqy*lmaxy,
     &              ',', Dqz*lmaxz ,')  (units of Bohr^-1 )'
         write(*,*) '(maximum sqr wf weight <',sqrwf_cutoff,')'
         write(*,*) 'excited state no.',lstore

         do lpr = 1, level
          if( (jindex(lpr).eq.lmaxy).and.(kindex(lpr).eq.lmaxz) )then
           write(13,*) iindex(lpr) * Dqx, H(lpr,l)
           write(14,*) iindex(lpr) * Dqx,
     &                 E_cv( iindex(lpr), lmaxy, lmaxz )
          endif
         enddo

         GO TO 20
        endif

       enddo

       write(*,*) 'plasmon wave function not found'
       write(*,*) 'lowest-energy eingenvalue as output'

       lstore = 1

      ELSE  ! semiconductor case

       lstore = 1

      ENDIF

 20   CONTINUE

      
      epsilon_u = sngl(eigen(lstore)) ! BEWARE of Hartree units
      write(*,*) 'energy =', epsilon_u * ee,'(eV)'

c     ground-state / plasmon wave function

      maxpsisqr = maxval( H(1:level,lstore)**2 )
      write(*,*) ' maximum value of psi square =', maxpsisqr

      do l = 1, level
       
       if( (jindex(l) .eq. 0 ).and. (kindex(l) .eq. 0 ) )then
        write(600,*) iindex(l)*Dqx, H(l,lstore)**2, H(l,lstore)
        write(700,*) iindex(l)*Dqx, 1./( (iindex(l)*Dqx)**2 +
     &        (mu/kappa_r)**2 )**4 * (mu/kappa_r)**8  * maxpsisqr
       endif
       if( (iindex(l) .eq. 0 ).and. (kindex(l) .eq. 0 ) )then
        write(601,*) jindex(l)*Dqy, H(l,lstore)**2, H(l,lstore)
        write(701,*) jindex(l)*Dqy, 1./( (jindex(l)*Dqy)**2 +
     &        (mu/kappa_r)**2 )**4 * (mu/kappa_r)**8  * maxpsisqr
       endif
       if( (iindex(l) .eq. 0 ).and. (jindex(l) .eq. 0 ) )then
        write(602,*) kindex(l)*Dqz, H(l,lstore)**2, H(l,lstore)
        write(702,*) kindex(l)*Dqz, 1./( (kindex(l)*Dqz)**2 +
     &        (mu/kappa_r)**2 )**4 * (mu/kappa_r)**8  * maxpsisqr
       endif

      enddo
c
c
c
c     Self-consistent calculation of EI order parameter
c     Preliminaries
c
c
c
      write(*,*)
      write(*,*) 'Self-consistent calculation of EI order parameter'
      write(*,*)
      write(*,*) 'Tolerance in interation cycle (eV) =',tolerance
      write(*,*) 'Number of iteration allowed =',n_iteration
c
c     IMPORTANT:
c     It works on the (k_x, k_y, k_z ) grid,
c     not on the one-dimensional l index
c

c     compute seed based on \psi_u(k_x, k_y, k_z)

      ALLOCATE(Delta_old( -N_x/2:N_x/2 , -N_y/2:N_y/2 , -N_z/2:N_z/2) )
      ALLOCATE(EEpsi_c( -N_x/2:N_x/2 , -N_y/2:N_y/2 , -N_z/2:N_z/2) )
      ALLOCATE(EEpsi_v( -N_x/2:N_x/2 , -N_y/2:N_y/2 , -N_z/2:N_z/2) )

      Delta_old = 0.
      do l = 1, level
       i = iindex(l)
       j = jindex(l)
       k = kindex(l)
       Delta_old( i, j, k ) = H(l,lstore) 
      enddo

      psi_u_max = MAXVAL( abs( Delta_old( -N_x/2 : N_x/2 , 
     &                    -N_y/2 : N_y/2, -N_z/2 : N_z/2 ) ) )
      Delta_old = ( E_cv - epsilon_u) * Delta_old / psi_u_max

      write(*,*) 'Lowest-exciton wf maximum (abs value) =',psi_u_max

      if( epsilon_u .lt. 0.)then
       Delta_trial = - epsilon_u
       factor_Maksimov_old = 1. / (1.+alpha_Maksimov/Delta_trial**2)
      else
       factor_Maksimov_old = 0.
      endif

      write(*,*) 'Trial factor_Maksimov =', factor_Maksimov_old

      ALLOCATE( Delta_new(-N_x/2:N_x/2, -N_y/2:N_y/2, -N_z/2:N_z/2 ) )
      ALLOCATE( phi_sum(-N_x/2:N_x/2, -N_y/2:N_y/2, -N_z/2:N_z/2 ) )
      ALLOCATE( Warray( -N_x : N_x, -N_y : N_y, -N_z : N_z ) )

      Warray = 0.

      do i = -N_x, N_x ! k_x loop
       do j = -N_y , N_y ! k_y loop
        do k = -N_z , N_z ! k_z loop

         if( E_xprime .ge. 0. )then ! semiconductor case

          if( (i.eq.0) .and. (j.eq.0) .and. (k.eq.0) )then
           Warray( i,j,k ) = Wq0_EM
          else
           qsq = Dqx**2 * i**2 + Dqy**2 * j**2 + Dqz**2 * k**2
           Warray( i,j,k ) = Dqx * Dqy * Dqz 
     &                       / kappa_r / 2. / pi**2 / qsq
           i_nn_sqrt = i**2 + j**2 + k**2
           if ( i_nn_sqrt .eq. 1 )then
            Warray( i,j,k ) = Warray( i,j,k ) * 4. / 8.
           elseif( i_nn_sqrt .eq. 2 )then
            Warray( i,j,k ) = Warray( i,j,k ) * 6. / 8.
           elseif( i_nn_sqrt .eq. 3 )then
            Warray( i,j,k ) = Warray( i,j,k ) * 7. / 8.
           endif
          endif

         else ! semimetal case

          if( (i.eq.0) .and. (j.eq.0) .and. (k.eq.0) )then
           Warray( i,j,k ) = Wq0_EM
          else
           qsq = Dqx**2 * i**2 + Dqy**2 * j**2 + Dqz**2 * k**2
           Warray( i,j,k ) = Dqx * Dqy * Dqz 
     &                       / kappa_r / 2. / pi**2 / qsq
           i_nn_sqrt = i**2 + j**2 + k**2
           if ( i_nn_sqrt .eq. 1 )then
            Warray( i,j,k ) = Warray( i,j,k ) * 4. / 8.
           elseif( i_nn_sqrt .eq. 2 )then
            Warray( i,j,k ) = Warray( i,j,k ) * 6. / 8.
           elseif( i_nn_sqrt .eq. 3 )then
            Warray( i,j,k ) = Warray( i,j,k ) * 7. / 8.
           endif
          endif

CCC          qsq = Dqx**2 * i**2 + Dqy**2 * j**2 + Dqz**2 * k**2
CCC          Warray( i,j,k ) = Dqx * Dqy * Dqz / 2. / pi**2 
CCC     &                    / ( kappa_r * qsq + 4. * pi * DOS_F ) 

         endif 

        enddo
       enddo
      enddo

c     Delta seed along x direction
      do i = -N_x/2, N_x/2   ! k_x loop
       write(228,*) Dqx * i, Delta_old( i, 0, 0 ), Warray( i, 0, 0 )
      enddo
c     Delta seed along y direction
      do j = -N_y/2, N_y/2   ! k_y loop
       write(218,*) Dqy * j, Delta_old( 0, j, 0 ), Warray( 0, j, 0 )
      enddo
c     Delta seed along z direction
      do k = -N_z/2, N_z/2   ! k_y loop
       write(208,*) Dqz * k, Delta_old( 0, 0, k ), Warray( 0, 0, k )
      enddo

      DO iT = 1, nintT   ! temperature cycle
       
       T = T_min + deltaT * (it - 1)
       write(*,*)
       write(*,*) ' T (K) = ',T
       kBT = T * 0.0862d0 * 1.e-3
       write(*,*) ' k_B T (eV) =',kBT
       write(*,*) ' k_B T (Hartree) =',kBT / ee
       if( kBT .gt. tolerance )then
               beta = 1. / kBT * ee  ! BEWARE of Hartree units
       else
               beta = 1.e15
       endif

       Delta_new = 0.

       if( epsilon_u .ge. 0. )then
        GO TO 10
       endif

c
c      Self-consistent cycle
c

       write(*,*)
       write(*,*) '...starting the self-consistent cycle'
       write(*,*)

       DO elle = 1, n_iteration ! begin iteration cycle
c
c       compute phi_sum
c    
        do i = -N_x/2, N_x/2   ! k_x loop
         do j = -N_y/2, N_y/2   ! k_y loop
          do k = -N_z/2, N_z/2   ! k_z loop

           epsi_c = .5 * ( epsi_c0(i,j,k) + epsi_v0(i,j,k) )
     &      +   sqrt( E_cv(i,j,k)**2 * .25 + Delta_old(i,j,k)**2 )
           epsi_v = .5 * ( epsi_c0(i,j,k) + epsi_v0(i,j,k) )
     &      -   sqrt( E_cv(i,j,k)**2 * .25 + Delta_old(i,j,k)**2 )

           IF( kBT .le. tolerance )THEN

            if( ( epsi_v - epsilon_F ) .le. 0. )then
             n_a = 1
            else
             n_a = 0
            endif

            if( ( epsi_c - epsilon_F ) .le. 0. )then
             n_b = 1
            else
             n_b = 0
            endif

            Tfact = ( n_a - n_b ) * Factor_Maksimov_old

ccc            Tfact = 1.

           ELSE

            fFermi_v = 1.d0 / 
     &                 ( 1.d0 + exp( beta * ( epsi_v - epsilon_F ) ) )
            fFermi_c = 1.d0 /
     &                 ( 1.d0 + exp( beta * ( epsi_c - epsilon_F ) ) )
            Tfact = ( fFermi_v - fFermi_c ) * Factor_Maksimov_old

CCC            Tfact = tanh( .5 * beta * 
CCC     &      sqrt( E_cv(i,j,k)**2 * .25 + Delta_old(i,j,k)**2 )  )

           ENDIF

           phi_sum(i,j,k) = Delta_old(i,j,k) * Tfact / 2. / 
     &       sqrt( E_cv(i,j,k)**2 * .25 + Delta_old(i,j,k)**2 )

          enddo
         enddo
        enddo

        do i = -N_x/2, N_x/2   ! k_x loop
         do j = -N_y/2, N_y/2   ! k_y loop
          do k = -N_z/2, N_z/2   ! k_z loop

           summ = 0.

           do ipr = -N_x/2, N_x/2
            do jpr = -N_y/2, N_y/2
             do kpr = -N_z/2, N_z/2
              
              summ = summ + Warray( i-ipr, j-jpr, k-kpr ) * 
     &              phi_sum( ipr, jpr, kpr )

             enddo
            enddo
           enddo

           Delta_new( i, j, k ) = summ

          enddo
         enddo
        enddo

        maxdiff = maxval( abs( abs(Delta_new) - abs(Delta_old) ) )
        write(*,*) 'Max. diff. in iteration (Hartree) =', maxdiff

        if( maxdiff .le. tolerance )then

         GO TO 10

        else

         Delta_old = Delta_new

         Delta_k_F_av = ( abs(Delta_new(i_F,0,0)) +
     &    abs(Delta_new(0,j_F,0)) +  abs(Delta_new(0,0,k_F)) ) *.3333
         Delta_trial = Delta_k_F_av
         factor_Maksimov_new = 1. / (1.+alpha_Maksimov/Delta_trial**2)
         factor_Maksimov_old = factor_Maksimov_new
         write(*,*)
         write(*,*) 'Current Delta(k_F) =',Delta_k_F_av*ee,'(eV)'
         write(*,*) 'Current factor_Maksimov =',factor_Maksimov_old

        endif

       ENDDO  !end iteration cycle

 10    CONTINUE

       IF( iT .eq. 1 )THEN

        do i = -N_x/2, N_x/2   ! k_x loop

         kappino = Dqx * i
         epsi_c = .5 * ( E_cbottom + E_vtop + kappino**2 / 2. * (
     &                    1./m_c_K + 1./m_v_K ) )
     &          +   sqrt( E_cv(i,0,0)**2 * .25 + Delta_new(i,0,0)**2 )
         epsi_v = .5 * ( E_cbottom + E_vtop + kappino**2 / 2. * (
     &                    1./m_c_K + 1./m_v_K ) )
     &          -   sqrt( E_cv(i,0,0)**2 * .25 + Delta_new(i,0,0)**2 )

         write(400,*) kappino, Delta_new( i, 0, 0 ) * ee
         write(500,*) kappino, epsi_c * ee, epsi_v * ee, 
     &                epsi_c0(i,0,0) * ee, epsi_v0(i,0,0) * ee

        enddo

        do j = -N_y/2, N_y/2   ! k_y loop

         kappino = Dqy * j
         epsi_c = .5 * ( E_cbottom + E_vtop + kappino**2 / 2. * (
     &                    1./m_c_perp + 1./m_v_perp ) )
     &          +   sqrt( E_cv(0,j,0)**2 * .25 + Delta_new(0,j,0)**2 )
         epsi_v = .5 * ( E_cbottom + E_vtop + kappino**2 / 2. * (
     &                    1./m_c_perp + 1./m_v_perp ) )
     &          -   sqrt( E_cv(0,j,0)**2 * .25 + Delta_new(0,j,0)**2 )

         write(401,*) kappino, Delta_new( 0, j, 0 ) * ee
         write(501,*) kappino, epsi_c * ee, epsi_v * ee, 
     &                epsi_c0(0,j,0) * ee, epsi_v0(0,j,0) * ee

        enddo

        do k = -N_z/2, N_z/2   ! k_z loop

         kappino = Dqz * k
         epsi_c = .5 * ( E_cbottom + E_vtop + kappino**2 / 2. * (
     &                    1./m_c_z + 1./m_v_z ) )
     &          +   sqrt( E_cv(0,0,k)**2 * .25 + Delta_new(0,0,k)**2 )
         epsi_v = .5 * ( E_cbottom + E_vtop + kappino**2 / 2. * (
     &                    1./m_c_z + 1./m_v_z ) )
     &          -   sqrt( E_cv(0,0,k)**2 * .25 + Delta_new(0,0,k)**2 )

         write(402,*) Dqz * k, Delta_new( 0, 0, k ) * ee
         write(502,*) kappino, epsi_c * ee, epsi_v * ee,
     &                epsi_c0(0,0,k) * ee, epsi_v0(0,0,k) * ee

        enddo

        do k = -N_z/2, N_z/2   ! k_z loop

         write(*,*) Dqz * k
         write(*,*)

         do i = -N_x/2, N_x/2
          do j = -N_y/2, N_y/2   ! k_y loop

           write(403,*) Dqx*i, Dqy*j, Delta_new( i, j, k ) * ee

          enddo
         enddo
        enddo




       ENDIF ! end check iT = 1

       EEpsi_c = .5 * ( epsi_c0 + epsi_v0 ) + sqrt( E_cv**2 * .25  +
     &           Delta_new**2 )
       EEpsi_v = .5 * ( epsi_c0 + epsi_v0 ) - sqrt( E_cv**2 * .25  +
     &           Delta_new**2 )

       write(*,*) 
       EEpsi_cmin = MINVAL( EEpsi_c ) * ee
       write(*,*) 'Conduction band min = ',EEpsi_cmin,'(eV)'
       i_min = MINLOC( EEpsi_c  )
       write(*,*) 'k_min = (', i_min(1) * Dqx , i_min(2) * Dqy ,
     &  i_min(3) * Dqz ,') (units of Bohr^-1)'
       EEpsi_vmax = MAXVAL( EEpsi_v ) * ee
       write(*,*) 'Valence band max = ',EEpsi_vmax,'(eV)'
       i_max = MAXLOC( EEpsi_v )
       write(*,*) 'k_max = (', i_max(1) * Dqx , i_max(2) * Dqy ,
     &  i_max(3) * Dqz ,') (units of Bohr^-1)'
       write(*,*)'Energy gap =',EEpsi_cmin - EEpsi_vmax,'(eV)'


       write(*,*)
       write(*,*) 'Average exciton number density per volume ='
       write(*,*) SUM( phi_sum**2 ) * Dqx * Dqy * Dqz * 12. /
     &            (2. * pi)**3, '(Bohr^-3)'
       write(*,*) 'including spin and valley degeneracy'

c
c      electric dipole calculation
c
       P_z = dipole_z * SUM( phi_sum ) * Dqx * Dqy * Dqz / pi**3 

       write(*,*)
       write(*,*) 'Electric dipole maximum value'
       write(*,*) 'within the pristine unitary cell at the origin = '
       write(*,*)  P_z, '(electric charge / Bohr^2)'
       write(*,*)  P_z / a_Bohr**2 , '(electric charge / A0^2)'
       write(*,*)  P_z / a_Bohr**2 * eA0 , '(Debye / A0^3)'

       write(333,*) T, maxval ( abs( Delta_new ) ) * ee  !! output in eV units
       write(334,*) T, P_z, P_z / a_Bohr**2 , P_z / a_Bohr**2 * eA0 
c      output in e / Bohr^2 ; e / A0^2 ; Debye / A0^3 
       write(335,*) T, EEpsi_cmin - EEpsi_vmax

      ENDDO  ! end temperature cycle

      DEALLOCATE( Delta_old, Delta_new, phi_sum, E_cv, Warray )
      DEALLOCATE( epsi_c0, epsi_v0, EEpsi_c, EEpsi_v )
      DEALLOCATE( iindex, jindex, kindex )

      write(*,*) ' Bye...'




      stop
 27   FORMAT(a12,f9.6,x,f9.6,x,f9.6,x,a30)
      END

