TABLE OF CONTENTS


ABINIT/m_paw_optics [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_optics

FUNCTION

  This module contains several routines related to conductivity:
    optical conductivity, X spectroscopy, linear susceptibility, ...

COPYRIGHT

 Copyright (C) 2018-2022 ABINIT group (SM,VR,FJ,MT,NB,PGhosh)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_paw_optics
24 
25  use defs_basis
26  use m_xmpi
27  use m_errors
28  use m_wffile
29  use m_abicore
30  use m_hdr
31  use m_dtset
32  use m_dtfil
33  use m_nctk
34 #ifdef HAVE_NETCDF
35  use netcdf
36 #endif
37 
38  use defs_datatypes, only : pseudopotential_type
39  use defs_abitypes,  only : MPI_type
40  use m_time,         only : timab
41  use m_io_tools,     only : open_file,get_unit,close_unit
42  use m_pawpsp,       only : pawpsp_read_corewf
43  use m_pawrad,       only : pawrad_type,pawrad_deducer0,simp_gen,nderiv_gen,poisson
44  use m_pawtab,       only : pawtab_type
45  use m_pawcprj,      only : pawcprj_type,pawcprj_alloc,pawcprj_get, &
46 &                           pawcprj_free,pawcprj_mpi_allgather
47  use m_pawang,       only : pawang_type
48  use m_paw_denpot,   only : pawdensities,pawkindensities,pawdenpot
49  use m_paw_an,       only : paw_an_type,paw_an_init,paw_an_free,paw_an_copy
50  use m_pawrhoij,     only : pawrhoij_type
51  use m_paw_ij,       only : paw_ij_type
52  use m_paw_onsite,   only : pawnabla_init,pawnabla_core_init
53  use m_paw_sphharm,  only : setnabla_ylm
54  use m_pawxc,        only : pawxc,pawxcm,pawxc_get_xclevel,pawxc_get_usekden
55  use m_mpinfo,       only : destroy_mpi_enreg,nullify_mpi_enreg,initmpi_seq,proc_distrb_cycle
56  use m_numeric_tools,only : kramerskronig
57  use m_geometry,     only : metric
58  use m_hide_lapack,  only : matrginv
59  use m_paral_atom,   only : get_my_atmtab,free_my_atmtab
60 
61  implicit none
62 
63  private
64 
65 !public procedures.
66  public :: optics_paw
67  public :: optics_paw_core
68  public :: linear_optics_paw
69 
70 !I/O parameters
71 !Set to true to force the use of netCDF when available
72 ! overriding the value of dtset%iomode
73  logical,parameter :: use_netcdf_forced=.true.
74 !Set to true to use netcdf-MPIIO when available
75  logical,parameter :: use_netcdf_mpiio=.true.
76 !Set to true to compute/write only half of the (n,m) dipoles
77  logical,parameter :: compute_half_dipoles=.true.
78 !Set to true to use unlimited dimensions in netCDF file (experimental)
79 !  This is not mandatory because we know exactly the amount of data to write
80 !    and seems to impact performances negatively...
81 !    Not compatible with compute_half_dipoles=.true.
82  logical,parameter :: use_netcdf_unlimited=.false.
83 
84 CONTAINS  !========================================================================================

m_paw_optics/linear_optics_paw [ Functions ]

[ Top ] [ m_paw_optics ] [ Functions ]

NAME

 linear_optics_paw

FUNCTION

 This program computes the elements of the optical frequency dependent
 linear susceptiblity using matrix elements <-i Nabla> obtained from a
 PAW ground state calculation. It uses formula 17 from Gadoc et al,
 Phys. Rev. B 73, 045112 (2006) [[cite:Gajdo2006]] together with a scissors correction. It uses
 a Kramers-Kronig transform to compute the real part from the imaginary part, and
 it will work on all types of unit cells. It outputs all tensor elements of
 both the real and imaginary parts.

INPUTS

  filnam: base of file names to read data from
  mpi_enreg: mpi set up variable, not used in this code

OUTPUT

  _real and _imag output files

NOTES

  This routine is not tested

SOURCE

1507  subroutine linear_optics_paw(filnam,filnam_out)
1508 
1509 !Arguments -----------------------------------
1510 !scalars
1511  character(len=fnlen),intent(in) :: filnam,filnam_out
1512 
1513 !Local variables-------------------------------
1514  integer,parameter :: master=0
1515  integer :: iomode,bantot,bdtot_index,fform1,headform
1516  integer :: iband,ierr,ii,ikpt,iom,iout,isppol,isym,jband,jj,me,mband
1517  integer :: method,mom,nband_k,nkpt,nspinor,nsppol,nsym,occopt,only_check
1518  integer :: rdwr,spaceComm,inpunt,reunt,imunt,wfunt
1519  integer,allocatable :: nband(:),symrel(:,:,:)
1520  real(dp) :: del,dom,fij,gdelta,omin,omax,paijpbij(2),mbpt_sciss,wij,ucvol
1521  real(dp) :: diffwp, diffwm
1522  real(dp) :: e2rot(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),rprimdinv(3,3),symd(3,3),symdinv(3,3)
1523  real(dp),allocatable :: e1(:,:,:),e2(:,:,:,:),epsilon_tot(:,:,:,:),eigen0(:),eig0_k(:)
1524  real(dp),allocatable :: kpts(:,:),occ(:),occ_k(:),oml1(:),wtk(:)
1525  complex(dpc),allocatable :: eps_work(:)
1526  character(len=fnlen) :: filnam1,filnam_gen
1527  character(len=500) :: msg
1528  type(hdr_type) :: hdr
1529  type(wffile_type) :: wff1
1530 !arrays
1531  real(dp),allocatable :: psinablapsi(:,:,:,:)
1532 
1533 ! *********************************************************************************
1534 
1535  DBG_ENTER("COLL")
1536 
1537 !write(std_out,'(a)')' Give the name of the output file ...'
1538 !read(std_in, '(a)') filnam_out
1539 !write(std_out,'(a)')' The name of the output file is :',filnam_out
1540 
1541 !Read data file
1542  if (open_file(filnam,msg,newunit=inpunt,form='formatted') /= 0 ) then
1543    ABI_ERROR(msg)
1544  end if
1545 
1546  rewind(inpunt)
1547  read(inpunt,*)
1548  read(inpunt,'(a)')filnam_gen       ! generic name for the files
1549  filnam1=trim(filnam_gen)//'_OPT' ! nabla matrix elements file
1550 
1551 !Open the Wavefunction and optic files
1552 !These default values are typical of sequential use
1553  iomode=IO_MODE_FORTRAN ; spaceComm=xmpi_comm_self; me=0
1554 
1555 ! Read the header of the optic files
1556  call hdr_read_from_fname(hdr, filnam1, fform1, spaceComm)
1557  call hdr%free()
1558  if (fform1 /= 610) then
1559    ABI_ERROR("Abinit8 requires an OPT file with fform = 610")
1560  end if
1561 
1562 !Open the conducti optic files
1563  wfunt = get_unit()
1564  call WffOpen(iomode,spaceComm,filnam1,ierr,wff1,master,me,wfunt)
1565 
1566 !Read the header from Ground state file
1567  rdwr=1
1568  call hdr_io(fform1,hdr,rdwr,wff1)
1569 
1570 !Extract info from the header
1571  headform=hdr%headform
1572  bantot=hdr%bantot
1573  nkpt=hdr%nkpt
1574  ABI_MALLOC(kpts,(3,nkpt))
1575  ABI_MALLOC(wtk,(nkpt))
1576  kpts(:,:)=hdr%kptns(:,:)
1577  wtk(:)=hdr%wtk(:)
1578  nspinor=hdr%nspinor
1579  nsppol=hdr%nsppol
1580  occopt=hdr%occopt
1581  rprimd(:,:)=hdr%rprimd(:,:)
1582  rprimdinv(:,:) = rprimd(:,:)
1583  call matrginv(rprimdinv,3,3) ! need the inverse of rprimd to symmetrize the tensors
1584  ABI_MALLOC(nband,(nkpt*nsppol))
1585  ABI_MALLOC(occ,(bantot))
1586  occ(1:bantot)=hdr%occ(1:bantot)
1587  nband(1:nkpt*nsppol)=hdr%nband(1:nkpt*nsppol)
1588  nsym=hdr%nsym
1589  ABI_MALLOC(symrel,(3,3,nsym))
1590  symrel(:,:,:)=hdr%symrel(:,:,:)
1591 
1592 !Get mband, as the maximum value of nband(nkpt)
1593  mband=maxval(nband(:))
1594 
1595 !get ucvol etc.
1596  iout = -1
1597  call metric(gmet,gprimd,iout,rmet,rprimd,ucvol)
1598 
1599  write(std_out,*)
1600  write(std_out,'(a,3f10.5,a)' )' rprimd(bohr)      =',rprimd(1:3,1)
1601  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,2)
1602  write(std_out,'(a,3f10.5,a)' )'                    ',rprimd(1:3,3)
1603  write(std_out,*)
1604  write(std_out,'(a,3f10.5,a)' )' rprimdinv         =',rprimdinv(1:3,1)
1605  write(std_out,'(a,3f10.5,a)' )'                    ',rprimdinv(1:3,2)
1606  write(std_out,'(a,3f10.5,a)' )'                    ',rprimdinv(1:3,3)
1607  write(std_out,'(a,2i8)')      ' nkpt,mband        =',nkpt,mband
1608 
1609 !get eigen0
1610  ABI_MALLOC(eigen0,(mband*nkpt*nsppol))
1611  read(wfunt)(eigen0(iband),iband=1,mband*nkpt*nsppol)
1612 
1613  read(inpunt,*)mbpt_sciss
1614  read(inpunt,*)dom,omin,omax,mom
1615  close(inpunt)
1616 
1617  ABI_MALLOC(oml1,(mom))
1618  ABI_MALLOC(e1,(3,3,mom))
1619  ABI_MALLOC(e2,(2,3,3,mom))
1620  ABI_MALLOC(epsilon_tot,(2,3,3,mom))
1621  ABI_MALLOC(eps_work,(mom))
1622  del=(omax-omin)/(mom-1)
1623  do iom=1,mom
1624    oml1(iom)=omin+dble(iom-1)*del
1625  end do
1626  write(std_out,'(a,i8,4f10.5,a)')' npts,omin,omax,width,mbpt_sciss      =',mom,omin,omax,dom,mbpt_sciss,' Ha'
1627 
1628  ABI_MALLOC(psinablapsi,(2,3,mband,mband))
1629 
1630 !loop over spin components
1631  do isppol=1,nsppol
1632    bdtot_index = 0
1633 !  loop over k points
1634    do ikpt=1,nkpt
1635 !
1636 !    number of bands for this k point
1637      nband_k=nband(ikpt+(isppol-1)*nkpt)
1638      ABI_MALLOC(eig0_k,(nband_k))
1639      ABI_MALLOC(occ_k,(nband_k))
1640 !    eigenvalues for this k-point
1641      eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
1642 !    occupation numbers for this k-point
1643      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1644 !    values of -i*nabla matrix elements for this k point
1645      psinablapsi=zero
1646      read(wfunt)((psinablapsi(1:2,1,iband,jband),iband=1,nband_k),jband=1,nband_k)
1647      read(wfunt)((psinablapsi(1:2,2,iband,jband),iband=1,nband_k),jband=1,nband_k)
1648      read(wfunt)((psinablapsi(1:2,3,iband,jband),iband=1,nband_k),jband=1,nband_k)
1649 
1650 !    occupation numbers for k-point
1651      occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
1652 !    accumulate e2 for this k point, Eq. 17 from PRB 73, 045112 (2006) [[cite:Gajdo2006]]
1653      do iband = 1, nband_k
1654        do jband = 1, nband_k
1655          fij = occ_k(iband) - occ_k(jband) !occ number difference
1656          wij = eig0_k(iband) - eig0_k(jband) !energy difference
1657          if (abs(fij) > zero) then ! only consider states of differing occupation numbers
1658            do ii = 1, 3
1659              do jj = 1, 3
1660                paijpbij(1) = psinablapsi(1,ii,iband,jband)*psinablapsi(1,jj,iband,jband) + &
1661 &               psinablapsi(2,ii,iband,jband)*psinablapsi(2,jj,iband,jband)
1662                paijpbij(2) = psinablapsi(2,ii,iband,jband)*psinablapsi(1,jj,iband,jband) - &
1663 &               psinablapsi(1,ii,iband,jband)*psinablapsi(2,jj,iband,jband)
1664                do iom = 1, mom
1665 !                original version
1666 !                diffw = wij + mbpt_sciss - oml1(iom) ! apply scissors term here
1667 !                gdelta = exp(-diffw*diffw/(4.0*dom*dom))/(2.0*dom*sqrt(pi)) ! delta fnc resolved as Gaussian
1668 !                e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(oml1(iom)*oml1(iom))
1669 !                e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(oml1(iom)*oml1(iom))
1670                  diffwm = wij - mbpt_sciss + oml1(iom) ! apply scissors term here
1671                  diffwp = wij + mbpt_sciss - oml1(iom) ! apply scissors term here
1672                  gdelta = exp(-diffwp*diffwp/(4.0*dom*dom))/(2.0*dom*sqrt(pi))
1673                  e2(1,ii,jj,iom) = e2(1,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(1)*gdelta/(wij*wij)
1674                  e2(2,ii,jj,iom) = e2(2,ii,jj,iom) - (4.0*pi*pi/ucvol)*wtk(ikpt)*fij*paijpbij(2)*gdelta/(wij*wij)
1675                end do ! end loop over spectral points
1676              end do ! end loop over jj = 1, 3
1677            end do ! end loop over ii = 1, 3
1678          end if ! end selection on fij /= 0
1679        end do ! end loop over jband
1680      end do ! end loop over iband
1681 
1682      ABI_FREE(eig0_k)
1683      ABI_FREE(occ_k)
1684      bdtot_index=bdtot_index+nband_k
1685    end do ! end loop over k points
1686  end do ! end loop over spin polarizations
1687 
1688 !here apply nsym symrel transformations to reconstruct full tensor from IBZ part
1689  epsilon_tot(:,:,:,:) = zero
1690  do isym = 1, nsym
1691    symd(:,:)=matmul(rprimd(:,:),matmul(symrel(:,:,isym),rprimdinv(:,:)))
1692    symdinv(:,:)=symd(:,:)
1693    call matrginv(symdinv,3,3)
1694    do iom = 1, mom
1695      e2rot(:,:)=matmul(symdinv(:,:),matmul(e2(1,:,:,iom),symd(:,:)))
1696      epsilon_tot(2,:,:,iom) = epsilon_tot(2,:,:,iom)+e2rot(:,:)/nsym
1697    end do
1698  end do
1699 
1700 !generate e1 from e2 via KK transforma
1701  method=0 ! use naive integration ( = 1 for simpson)
1702  only_check=0 ! compute real part of eps in kk routine
1703  do ii = 1, 3
1704    do jj = 1, 3
1705      eps_work(:) = cmplx(0.0,epsilon_tot(2,ii,jj,:), kind=dpc)
1706      call kramerskronig(mom,oml1,eps_work,method,only_check)
1707      epsilon_tot(1,ii,jj,:) = real(eps_work(:))
1708      if (ii /= jj) epsilon_tot(1,ii,jj,:) = epsilon_tot(1,ii,jj,:)- 1.0
1709    end do ! end loop over jj
1710  end do ! end loop over ii
1711 
1712  if (open_file(trim(filnam_out)//'_imag',msg,newunit=reunt,form='formatted') /= 0) then
1713    ABI_ERROR(msg)
1714  end if
1715 
1716  if (open_file(trim(filnam_out)//'_real',msg,unit=imunt,form='formatted') /= 0) then
1717    ABI_ERROR(msg)
1718  end if
1719 
1720  write(reunt,'(a12,6a13)')' # Energy/Ha ','eps_2_xx','eps_2_yy','eps_2_zz',&
1721 & 'eps_2_yz','eps_2_xz','eps_2_xy'
1722  write(imunt,'(a12,6a13)')' # Energy/Ha ','eps_1_xx','eps_1_yy','eps_1_zz',&
1723 & 'eps_1_yz','eps_1_xz','eps_1_xy'
1724 
1725  do iom = 1, mom
1726    write(reunt,'(ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4)') oml1(iom),' ',&
1727 &   epsilon_tot(2,1,1,iom),' ',epsilon_tot(2,2,2,iom),' ',epsilon_tot(2,3,3,iom),' ',&
1728 &   epsilon_tot(2,2,3,iom),' ',epsilon_tot(2,1,3,iom),' ',epsilon_tot(2,1,2,iom)
1729    write(imunt,'(ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4,a,ES12.4)') oml1(iom),' ',&
1730 &   epsilon_tot(1,1,1,iom),' ',epsilon_tot(1,2,2,iom),' ',epsilon_tot(1,3,3,iom),' ',&
1731 &   epsilon_tot(1,2,3,iom),' ',epsilon_tot(1,1,3,iom),' ',epsilon_tot(1,1,2,iom)
1732  end do
1733 
1734  close(reunt)
1735  close(imunt)
1736 
1737  ABI_FREE(nband)
1738  ABI_FREE(oml1)
1739  ABI_FREE(e2)
1740  ABI_FREE(e1)
1741  ABI_FREE(occ)
1742  ABI_FREE(psinablapsi)
1743  ABI_FREE(eigen0)
1744  ABI_FREE(wtk)
1745  ABI_FREE(kpts)
1746 
1747  call hdr%free()
1748 
1749  DBG_EXIT("COLL")
1750 
1751  end subroutine linear_optics_paw

m_paw_optics/optics_paw [ Functions ]

[ Top ] [ m_paw_optics ] [ Functions ]

NAME

 optics_paw

FUNCTION

 Compute matrix elements need for optical conductivity (in the PAW context) and store them in a file
  Matrix elements = <Phi_i|Nabla|Phi_j>

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cg(2,mcg)=planewave coefficients of wavefunctions.
  cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk>
                                          and each |p_lmn> non-local projector
  dimcprj(natom)=array of dimensions of array cprj (not ordered)
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  gprimd(3,3)=dimensional reciprocal space primitive translations
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mkmem =number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpsang =1+maximum angular momentum for nonlocal pseudopotentials
  mpw=maximum dimensioned size of npw.
  natom=number of atoms in cell.
  nkpt=number of k points.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nsppol=1 for unpolarized, 2 for spin-polarized
  pawang <type(pawang_type)>= PAW ANGular mesh discretization and related data
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  znucl(ntypat)=atomic number of atom type

OUTPUT

  (only writing in a file)

SIDE EFFECTS

NOTES

SOURCE

133  subroutine optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen0,gprimd,hdr,kg,&
134 &               mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,&
135 &               pawang,pawrad,pawrhoij,pawtab,znucl)
136 
137 !Arguments ------------------------------------
138 !scalars
139  integer,intent(in) :: mband,mcg,mcprj,mkmem,mpsang,mpw,natom,nkpt,nsppol
140  type(MPI_type),intent(in) :: mpi_enreg
141  type(datafiles_type),intent(in) :: dtfil
142  type(dataset_type),intent(in) :: dtset
143  type(hdr_type),intent(inout) :: hdr
144  type(pawang_type),intent(in) :: pawang
145 !arrays
146  integer,intent(in) :: atindx1(natom),dimcprj(natom),npwarr(nkpt)
147  integer,intent(in),target :: kg(3,mpw*mkmem)
148  real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
149  real(dp),intent(in) :: gprimd(3,3),znucl(dtset%ntypat)
150  real(dp),intent(inout) :: cg(2,mcg)
151  type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
152  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
153  type(pawrhoij_type),intent(inout) :: pawrhoij(mpi_enreg%my_natom)
154  type(pawtab_type),target,intent(inout) :: pawtab(dtset%ntypat)
155 
156 !Local variables-------------------------------
157 !scalars
158  integer,parameter :: master=0
159  integer :: bsize,iomode,bdtot_index,cplex,etiq,fformopt,iatom,ib,ibmax,ibg,ibsp
160  integer :: ibshift,icg,ierr,ikg,ikpt,ilmn,ount,ncid,varid,idir
161  integer :: iorder_cprj,ipw,ispinor,isppol,istwf_k,itypat,iwavef
162  integer :: jb,jbshift,jbsp,my_jb,jlmn,jwavef,lmn_size,mband_cprj,option_core
163  integer :: my_nspinor,nband_k,nband_cprj_k,npw_k,sender,me,master_spfftband,pnp_size
164  integer :: spaceComm_band,spaceComm_bandspinorfft,spaceComm_fft,spaceComm_kpt
165  integer :: spaceComm_spinor,spaceComm_bandspinor,spaceComm_spinorfft,spaceComm_w
166  logical :: already_has_nabla,cprj_paral_band,myband,mykpt,iomode_etsf_mpiio
167  logical :: i_am_master,i_am_master_kpt,i_am_master_band,i_am_master_spfft,nc_unlimited,store_half_dipoles
168  real(dp) :: cgnm1,cgnm2,cpnm1,cpnm2,cpnm11,cpnm22,cpnm12,cpnm21,cpnm_11m22,cpnm_21p12,cpnm_21m12
169  character(len=500) :: msg
170  type(nctkdim_t) :: nctkdim
171 !arrays
172  integer :: nc_count_5(5),nc_count_6(6),nc_start_5(5),nc_start_6(6),nc_stride_5(5),nc_stride_6(6),tmp_shape(3)
173  integer, ABI_CONTIGUOUS pointer :: kg_k(:,:)
174  real(dp) :: kpoint(3),tsec(2),nabla_ij(3)
175  real(dp),allocatable :: kpg_k(:,:)
176  real(dp),allocatable :: psinablapsi(:,:,:),psinablapsi_paw(:,:,:),psinablapsi_soc(:,:,:)
177  real(dp),pointer :: soc_ij(:,:,:)
178  type(coeff5_type),allocatable,target :: phisocphj(:)
179  type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
180  type(nctkarr_t) :: nctk_arrays(1)
181 
182 ! ************************************************************************
183 
184  DBG_ENTER("COLL")
185 
186 !Compatibility tests
187  ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
188  ABI_CHECK(mpi_enreg%paral_spinor==0.or.dtset%pawspnorb==0,"spinor parallelization not supported with SOC!")
189 !  MJV 6/12/2008: looks like mpi_enreg may not be completely initialized here
190  if (xmpi_paral==1) then
191    tmp_shape = shape(mpi_enreg%proc_distrb)
192    if (nkpt > tmp_shape(1)) then
193      ABI_BUG('problem with proc_distrb!')
194    end if
195  end if
196 
197 !Init parallelism
198  spaceComm_w=mpi_enreg%comm_cell
199  if (mpi_enreg%paral_kgb==1) then
200    spaceComm_kpt=mpi_enreg%comm_kpt
201    spaceComm_fft=mpi_enreg%comm_fft
202    spaceComm_band=mpi_enreg%comm_band
203    spaceComm_spinor=mpi_enreg%comm_spinor
204    spaceComm_bandspinor=mpi_enreg%comm_bandspinor
205    spaceComm_spinorfft=mpi_enreg%comm_spinorfft
206    spaceComm_bandspinorfft=mpi_enreg%comm_bandspinorfft
207  else
208    spaceComm_kpt=mpi_enreg%comm_kpt
209    spaceComm_fft=xmpi_comm_self
210    spaceComm_band=mpi_enreg%comm_band
211    spaceComm_spinor=xmpi_comm_self
212    spaceComm_bandspinor=spaceComm_band
213    spaceComm_spinorfft=xmpi_comm_self
214    spaceComm_bandspinorfft=xmpi_comm_self
215  end if
216  me=xmpi_comm_rank(spaceComm_w)
217  i_am_master=(me==master)
218  i_am_master_kpt=(xmpi_comm_rank(spaceComm_kpt)==master)
219  i_am_master_band=(xmpi_comm_rank(spaceComm_band)==master)
220  i_am_master_spfft=(xmpi_comm_rank(spaceComm_spinorfft)==master)
221  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
222 
223 !----------------------------------------------------------------------------------
224 !1- Opening of OPT file and header writing
225 !----------------------------------------------------------------------------------
226 
227 !I/O mode is netCDF or Fortran
228  iomode=merge(IO_MODE_ETSF,IO_MODE_FORTRAN_MASTER,dtset%iomode==IO_MODE_ETSF)
229 #ifdef HAVE_NETCDF
230  if (use_netcdf_forced) iomode=IO_MODE_ETSF
231 #endif
232 
233  !(master proc only)
234  if (i_am_master) then
235    fformopt=610 ; if (compute_half_dipoles) fformopt=620
236 !  ====> NETCDF format
237    if (iomode==IO_MODE_ETSF) then
238 #ifdef HAVE_NETCDF
239 !    Open/create nc file
240      NCF_CHECK(nctk_open_create(ncid,nctk_ncify(dtfil%fnameabo_app_opt),xmpi_comm_self))
241 !    Write header data
242      NCF_CHECK(hdr%ncwrite(ncid,fformopt,nc_define=.true.))
243 !    Define dims and array for dipole variables
244      nctk_arrays(1)%name="dipole_valence_valence"
245      nctk_arrays(1)%dtype="dp"
246      nc_unlimited=(use_netcdf_unlimited.and.(.not.(nctk_has_mpiio.and.use_netcdf_mpiio)))
247      if (nc_unlimited) then
248        nctkdim%name="unlimited_bands"
249        nctkdim%value=NF90_UNLIMITED
250        NCF_CHECK(nctk_def_dims(ncid,nctkdim))
251        nctk_arrays(1)%shape_str=&
252 &      "complex,number_of_cartesian_directions,max_number_of_states,number_of_kpoints,number_of_spins,unlimited_bands"
253      else if (compute_half_dipoles) then
254        nctkdim%name="max_number_of_state_pairs"
255        nctkdim%value=(mband*(mband+1))/2
256        NCF_CHECK(nctk_def_dims(ncid,nctkdim))
257        nctk_arrays(1)%shape_str=&
258 &      "complex,number_of_cartesian_directions,max_number_of_state_pairs,number_of_kpoints,number_of_spins"
259      else
260        nctk_arrays(1)%shape_str=&
261 &      "complex,number_of_cartesian_directions,max_number_of_states,max_number_of_states,number_of_kpoints,number_of_spins"
262      end if
263      NCF_CHECK(nctk_def_arrays(ncid, nctk_arrays))
264      NCF_CHECK(nctk_set_atomic_units(ncid, "dipole_valence_valence"))
265 !    Write eigenvalues
266      NCF_CHECK(nctk_set_datamode(ncid))
267      varid=nctk_idname(ncid,"eigenvalues")
268      NCF_CHECK(nf90_put_var(ncid,varid,reshape(eigen0,[mband,nkpt,nsppol])))
269      !Close file here because the rest has possibly to be written with collective I/O
270      NCF_CHECK(nf90_close(ncid))
271 #else
272      msg = "In order to use iomode=3 and prtnabla>0 together, NetCDF support must be enabled!"
273      ABI_ERROR(msg)
274 #endif
275 !  ====> Standard FORTRAN binary format
276    else if (iomode==IO_MODE_FORTRAN_MASTER) then
277      if (open_file(dtfil%fnameabo_app_opt,msg,newunit=ount,form="unformatted",status="unknown")/= 0) then
278        ABI_ERROR(msg)
279      end if
280      call hdr%fort_write(ount,fformopt,ierr,rewind=.true.)
281      write(ount)(eigen0(ib),ib=1,mband*nkpt*nsppol)
282    else
283      msg = "Wrong OPT file format!"
284      ABI_BUG(msg)
285    end if ! File format
286  end if ! master node
287  call xmpi_bcast(iomode,master,spaceComm_w,ierr)  ! Seems mandatory; why ?
288  iomode_etsf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
289  nc_unlimited=(iomode==IO_MODE_ETSF.and.use_netcdf_unlimited.and.(.not.iomode_etsf_mpiio)) ! UNLIMITED not compatible with mpi-io
290  store_half_dipoles=(compute_half_dipoles.and.(.not.nc_unlimited))
291 
292 !----------------------------------------------------------------------------------
293 !2- Computation of on-site contribution: <phi_i|nabla|phi_j>-<tphi_i|nabla|tphi_j>
294 !----------------------------------------------------------------------------------
295 
296  already_has_nabla=all(pawtab(:)%has_nabla==2)
297  call pawnabla_init(mpsang,dtset%ntypat,pawrad,pawtab)
298  
299 !Compute spin-orbit contributions if necessary
300  if (dtset%pawspnorb==1) then
301    option_core=0
302    call pawnabla_soc_init(phisocphj,option_core,dtset%ixc,mpi_enreg%my_natom,natom,&
303 &       dtset%nspden,dtset%ntypat,pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,&
304 &       dtset%spnorbscl,dtset%typat,dtset%xc_denpos,znucl,&
305 &       comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
306  end if
307 
308 !----------------------------------------------------------------------------------
309 !3- Computation of <psi_n|-i.nabla|psi_m> for each k
310 !----------------------------------------------------------------------------------
311 
312 !Prepare valence-valence dipoles writing
313 !In case of netCDF access to OPT file, prepare collective I/O
314  if (iomode == IO_MODE_ETSF) then
315    if (iomode_etsf_mpiio) then
316      if (i_am_master_spfft) then
317        NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt),spaceComm_band))
318        varid=nctk_idname(ncid,"dipole_valence_valence")
319        if (xmpi_comm_size(spaceComm_w)>1) then
320          NCF_CHECK(nctk_set_collective(ncid,varid))
321        end if
322        NCF_CHECK(nctk_set_datamode(ncid))
323      end if
324    else if (i_am_master) then
325      NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt),xmpi_comm_self))
326      varid=nctk_idname(ncid,"dipole_valence_valence")
327      if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
328        NCF_CHECK(nctk_set_collective(ncid,varid))
329      end if
330      NCF_CHECK(nctk_set_datamode(ncid))
331    end if     
332  end if
333  if (iomode_etsf_mpiio) then
334    !If MPI-IO, store only ib elements for each jb
335    ABI_MALLOC(psinablapsi,(2,3,mband))
336    ABI_MALLOC(psinablapsi_paw,(2,3,mband))
337    if (dtset%pawspnorb==1) then
338      ABI_MALLOC(psinablapsi_soc,(2,3,mband))
339    end if
340  else
341    !If not, store all (ib,jb) pairs (or half)
342    bsize=mband**2 ; if (store_half_dipoles) bsize=(mband*(mband+1))/2
343    ABI_MALLOC(psinablapsi,(2,3,bsize))
344    ABI_MALLOC(psinablapsi_paw,(2,3,bsize))
345    if (dtset%pawspnorb==1) then
346      ABI_MALLOC(psinablapsi_soc,(2,3,bsize))
347    end if
348    psinablapsi=zero
349  end if
350  pnp_size=size(psinablapsi)
351  
352 !Determine if cprj datastructure is distributed over bands
353  mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
354  cprj_paral_band=(mband_cprj<mband)
355 
356 !LOOP OVER SPINS
357  ibg=0;icg=0
358  bdtot_index=0
359  do isppol=1,nsppol
360 
361 !  LOOP OVER k POINTS
362    ikg=0
363    do ikpt=1,nkpt
364 
365      etiq=ikpt+(isppol-1)*nkpt
366      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
367      master_spfftband=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
368 
369 !    Select k-points for current proc
370      mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me_kpt))
371      if (mykpt) then
372 
373 !      Data depending on k-point
374        npw_k=npwarr(ikpt)
375        istwf_k=dtset%istwfk(ikpt)
376        cplex=2;if (istwf_k>1) cplex=1
377        kpoint(:)=dtset%kptns(:,ikpt)
378 
379 !      Extract cprj for this k-point
380        nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
381        if (mkmem*nsppol/=1) then
382          iorder_cprj=0
383          ABI_MALLOC(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
384          call pawcprj_alloc(cprj_k_loc,0,dimcprj)
385          call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
386 &         mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
387 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
388        else
389          cprj_k_loc => cprj
390        end if
391 
392 !      if cprj are distributed over bands, gather them (because we need to mix bands)
393        if (cprj_paral_band) then
394          ABI_MALLOC(cprj_k,(natom,my_nspinor*nband_k))
395          call pawcprj_alloc(cprj_k,0,dimcprj)
396          call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom, &
397 &             my_nspinor*nband_cprj_k,my_nspinor*mpi_enreg%bandpp,&
398 &         dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
399        else
400          cprj_k => cprj_k_loc
401        end if
402 
403 !      Compute k+G in cartesian coordinates
404        ABI_MALLOC(kpg_k,(3,npw_k*dtset%nspinor))
405        kg_k => kg(:,1+ikg:npw_k+ikg)
406        do ipw=1,npw_k
407          kpg_k(1,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(1,1) &
408 &                    +(kpoint(2)+kg_k(2,ipw))*gprimd(1,2) &
409 &                    +(kpoint(3)+kg_k(3,ipw))*gprimd(1,3)
410          kpg_k(2,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(2,1) &
411 &                    +(kpoint(2)+kg_k(2,ipw))*gprimd(2,2) &
412 &                    +(kpoint(3)+kg_k(3,ipw))*gprimd(2,3)
413          kpg_k(3,ipw)=(kpoint(1)+kg_k(1,ipw))*gprimd(3,1) &
414 &                    +(kpoint(2)+kg_k(2,ipw))*gprimd(3,2) &
415 &                    +(kpoint(3)+kg_k(3,ipw))*gprimd(3,3)
416        end do
417        kpg_k=two_pi*kpg_k
418        if (dtset%nspinor==2) kpg_k(1:3,npw_k+1:2*npw_k)=kpg_k(1:3,1:npw_k)
419 
420 !      Loops over bands
421        do jb=1,nband_k
422          jwavef=(jb-1)*npw_k*my_nspinor+icg
423          !If MPI-IO, compute all (ib,jb) pairs ; if not, compute only ib<=jb
424          ibmax=merge(nband_k,jb,iomode_etsf_mpiio)
425          !If MPI-IO, store only ib elements for each jb ; if not, store all (ib,jb) pairs
426          my_jb=merge(1,jb,iomode_etsf_mpiio)
427 
428 !        Fill output arrays with zeros
429          if (store_half_dipoles) then
430            jbshift=(my_jb*(my_jb-1))/2 ; bsize=nband_k
431          else
432            jbshift=(my_jb-1)*mband ; bsize=mband
433          end if
434          psinablapsi(:,:,jbshift+1:jbshift+bsize)=zero
435          psinablapsi_paw(:,:,jbshift+1:jbshift+bsize)=zero
436          if (dtset%pawspnorb==1) psinablapsi_soc(:,:,jbshift+1:jbshift+bsize)=zero
437            
438 !        2-A Computation of <psi_tild_n|-i.nabla|psi_tild_m>
439 !        ----------------------------------------------------------------------------------
440 !        Note: <psi_nk|-i.nabla|psi_mk> => Sum_g[ <G|Psi_nk>^* <G|Psi_mk> G ]
441 
442 !        Select bands for current proc
443          myband=.true.
444          if (xmpi_paral==1.and.mpi_enreg%paral_kgb/=1.and.(.not.iomode_etsf_mpiio)) then
445            myband=(abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-mpi_enreg%me_kpt)==0)
446          end if
447          if (myband) then
448 
449            do ib=1,ibmax
450              iwavef=(ib-1)*npw_k*my_nspinor+icg
451 
452 !            (C_nk^*)*C_mk*(k+g) is expressed in cartesian coordinates
453              if (istwf_k>1) then
454                !In this case (istwfk>1): Sum_g>=g0[ 2i.Im(<G|Psi_nk>^* <G|Psi_mk> G) ]
455                !G=k+g=0 term is included but do not contribute
456                do ipw=1,npw_k*my_nspinor
457                  cgnm2=two*(cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef))
458                  psinablapsi(2,1:3,jbshift+ib)=psinablapsi(2,1:3,jbshift+ib)+cgnm2*kpg_k(1:3,ipw)
459                end do
460              else
461                do ipw=1,npw_k*my_nspinor
462                  cgnm1=cg(1,ipw+iwavef)*cg(1,ipw+jwavef)+cg(2,ipw+iwavef)*cg(2,ipw+jwavef)
463                  cgnm2=cg(1,ipw+iwavef)*cg(2,ipw+jwavef)-cg(2,ipw+iwavef)*cg(1,ipw+jwavef)
464                  psinablapsi(1,1:3,jbshift+ib)=psinablapsi(1,1:3,jbshift+ib)+cgnm1*kpg_k(1:3,ipw)
465                  psinablapsi(2,1:3,jbshift+ib)=psinablapsi(2,1:3,jbshift+ib)+cgnm2*kpg_k(1:3,ipw)
466                end do
467              end if
468 
469 !            Second half of the (n,m) matrix
470              if ((ib/=jb).and.(.not.store_half_dipoles).and.(.not.iomode_etsf_mpiio)) then
471                ibshift=(ib-1)*mband
472                psinablapsi(1,1:3,ibshift+jb)= psinablapsi(1,1:3,jbshift+ib)
473                psinablapsi(2,1:3,ibshift+jb)=-psinablapsi(2,1:3,jbshift+ib)
474              end if
475 
476            end do ! ib
477 
478 !          Reduction in case of parallelism
479            if (iomode_etsf_mpiio) then
480              call timab(48,1,tsec)
481              if (mpi_enreg%paral_kgb==1) then
482                call xmpi_sum_master(psinablapsi,master,spaceComm_spinorfft,ierr)
483              end if
484              if (i_am_master_spfft) then
485                call xmpi_sum(psinablapsi,spaceComm_band,ierr)
486              end if
487              call timab(48,2,tsec)
488            end if
489 
490          end if ! myband
491 
492        
493 !        2-B Computation of <psi_n|p_i><p_j|psi_m>(<phi_i|-i.nabla|phi_j>-<tphi_i|-i.nabla|tphi_j>)
494 !        ----------------------------------------------------------------------------------
495 !        Non relativistic contribution
496 !        Note: <psi|-i.nabla|psi_mk>
497 !              => -i Sum_ij[ <p_i|Psi_nk>^* <p_j|Psi_mk> (<Phi_i|Nabla|Phi_j>-<tPhi_i|Nabla|tPhi_j>) ]
498 
499 !        Select bands for current proc
500          myband=.true.
501          if (mpi_enreg%paral_kgb==1) then
502            myband=(mod(jb-1,mpi_enreg%nproc_band)==mpi_enreg%me_band)
503          else if (xmpi_paral==1) then
504            myband=(abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-mpi_enreg%me_kpt)==0)
505          end if
506          if (myband) then
507 
508            do ib=1,ibmax          
509 
510              ibsp=(ib-1)*my_nspinor ; jbsp=(jb-1)*my_nspinor
511              do ispinor=1,my_nspinor
512                ibsp=ibsp+1;jbsp=jbsp+1
513                if (cplex==1) then
514                  do iatom=1,natom
515                    itypat=dtset%typat(iatom)
516                    lmn_size=pawtab(itypat)%lmn_size
517                    do jlmn=1,lmn_size
518                      do ilmn=1,lmn_size                  
519                        nabla_ij(:)=pawtab(itypat)%nabla_ij(:,ilmn,jlmn)
520                        if (ib>jb) nabla_ij(:)=-pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
521                        cpnm1=cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn)
522                        if (dtset%nspinor==2) cpnm1=cpnm1+cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn)
523                        psinablapsi_paw(2,:,jbshift+ib)=psinablapsi_paw(2,:,jbshift+ib)-cpnm1*nabla_ij(:)
524                      end do !ilmn
525                    end do !jlmn
526                  end do !iatom
527                else
528                  do iatom=1,natom
529                    itypat=dtset%typat(iatom)
530                    lmn_size=pawtab(itypat)%lmn_size
531                    do jlmn=1,lmn_size
532                      do ilmn=1,lmn_size
533                        nabla_ij(:)=pawtab(itypat)%nabla_ij(:,ilmn,jlmn)
534                        if (ib>jb) nabla_ij(:)=-pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
535                        cpnm1=(cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn) &
536 &                            +cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn))
537                        cpnm2=(cprj_k(iatom,ibsp)%cp(1,ilmn)*cprj_k(iatom,jbsp)%cp(2,jlmn) &
538 &                            -cprj_k(iatom,ibsp)%cp(2,ilmn)*cprj_k(iatom,jbsp)%cp(1,jlmn))
539                        psinablapsi_paw(1,:,jbshift+ib)=psinablapsi_paw(1,:,jbshift+ib)+cpnm2*nabla_ij(:)
540                        psinablapsi_paw(2,:,jbshift+ib)=psinablapsi_paw(2,:,jbshift+ib)-cpnm1*nabla_ij(:)
541                      end do !ilmn
542                    end do !jlmn
543                  end do !iatom
544                end if
545              end do !ispinor
546 
547 !        2-C Computation of Spin-orbit coupling contribution:
548 !             Sum_ij,ss'[<psi_n,s|p_i><p_j|psi_m,s'>(<phi_i|1/4 Alpha^2 (Sigma^ss' X dV/dr)|phi_j>]
549 !        ----------------------------------------------------------------------------------
550              if (dtset%pawspnorb==1) then
551 !              Add: Sum_ij,ss'[<Psi^s_n|p_i><p_j|Psi^s'_m> (Sigma^ss' X g_ij)]
552 !                  =Sum_ij[ (<Psi^up_n|p_i><p_j|Psi^up_m>-<Psi^dn_n|p_i><p_j|Psi^dn_m>) (Sigma^up-up     X g_ij)
553 !                          +(<Psi^dn_n|p_i><p_j|Psi^up_m>+<Psi^up_n|p_i><p_j|Psi^dn_m>) (Re{Sigma^dn-up} X g_ij)
554 !                          +(<Psi^dn_n|p_i><p_j|Psi^up_m>-<Psi^up_n|p_i><p_j|Psi^dn_m>) (Im{Sigma^dn-up} X g_ij) ]
555 !               where: g_ij = <Phi_i| 1/4 Alpha^2 dV/dr vec(r)/r) |Phi_j>
556 !              Note that:
557 !                phisocphj(:)%value(re:im,1,idir,ilmn,jlmn) is (Sigma^up-up X g_ij)
558 !                phisocphj(:)%value(re:im,2,idir,ilmn,jlmn) is (Sigma^dn-up X g_ij)
559 !              Not compatible with parallelization over spinors
560                ibsp=1+(ib-1)*dtset%nspinor ; jbsp=1+(jb-1)*dtset%nspinor
561                do iatom=1,natom
562                  itypat=dtset%typat(iatom)
563                  lmn_size=pawtab(itypat)%lmn_size
564                  do jlmn=1,lmn_size
565                    do ilmn=1,lmn_size
566                      soc_ij => phisocphj(iatom)%value(:,:,:,ilmn,jlmn)
567                      !Contribution from real part of <Psi^s_n|p_i><p_j|Psi^s'_m>
568                      cpnm11=cprj_k(iatom,ibsp  )%cp(1,ilmn)*cprj_k(iatom,jbsp  )%cp(1,jlmn) &
569 &                          +cprj_k(iatom,ibsp  )%cp(2,ilmn)*cprj_k(iatom,jbsp  )%cp(2,jlmn)
570                      cpnm22=cprj_k(iatom,ibsp+1)%cp(1,ilmn)*cprj_k(iatom,jbsp+1)%cp(1,jlmn) &
571 &                          +cprj_k(iatom,ibsp+1)%cp(2,ilmn)*cprj_k(iatom,jbsp+1)%cp(2,jlmn)
572                      cpnm12=cprj_k(iatom,ibsp  )%cp(1,ilmn)*cprj_k(iatom,jbsp+1)%cp(1,jlmn) &
573 &                          +cprj_k(iatom,ibsp  )%cp(2,ilmn)*cprj_k(iatom,jbsp+1)%cp(2,jlmn)
574                      cpnm21=cprj_k(iatom,ibsp+1)%cp(1,ilmn)*cprj_k(iatom,jbsp  )%cp(1,jlmn) &
575 &                          +cprj_k(iatom,ibsp+1)%cp(2,ilmn)*cprj_k(iatom,jbsp  )%cp(2,jlmn)
576                      cpnm_11m22=cpnm11-cpnm22
577                      cpnm_21p12=cpnm21+cpnm12
578                      cpnm_21m12=cpnm21-cpnm12
579                      do idir=1,3
580                        psinablapsi_soc(1,idir,jbshift+ib)=psinablapsi_soc(1,idir,jbshift+ib) &
581 &                             +soc_ij(1,1,idir)*cpnm_11m22+soc_ij(1,2,idir)*cpnm_21p12
582                        psinablapsi_soc(2,idir,jbshift+ib)=psinablapsi_soc(2,idir,jbshift+ib) &
583 &                             +soc_ij(2,2,idir)*cpnm_21m12
584                      end do
585                      !Contribution from imaginary part of <Psi^s_n|p_i><p_j|Psi^s'_m>
586                      cpnm11=cprj_k(iatom,ibsp  )%cp(1,ilmn)*cprj_k(iatom,jbsp  )%cp(2,jlmn) &
587 &                          -cprj_k(iatom,ibsp  )%cp(2,ilmn)*cprj_k(iatom,jbsp  )%cp(1,jlmn)
588                      cpnm22=cprj_k(iatom,ibsp+1)%cp(1,ilmn)*cprj_k(iatom,jbsp+1)%cp(2,jlmn) &
589 &                          -cprj_k(iatom,ibsp+1)%cp(2,ilmn)*cprj_k(iatom,jbsp+1)%cp(1,jlmn)
590                      cpnm12=cprj_k(iatom,ibsp  )%cp(1,ilmn)*cprj_k(iatom,jbsp+1)%cp(2,jlmn) &
591 &                          -cprj_k(iatom,ibsp  )%cp(2,ilmn)*cprj_k(iatom,jbsp+1)%cp(1,jlmn)
592                      cpnm21=cprj_k(iatom,ibsp+1)%cp(1,ilmn)*cprj_k(iatom,jbsp  )%cp(2,jlmn) &
593 &                          -cprj_k(iatom,ibsp+1)%cp(2,ilmn)*cprj_k(iatom,jbsp  )%cp(1,jlmn)
594                      cpnm_11m22=cpnm11-cpnm22
595                      cpnm_21p12=cpnm21+cpnm12
596                      cpnm_21m12=cpnm21-cpnm12
597                      do idir=1,3
598                        psinablapsi_soc(1,idir,jbshift+ib)=psinablapsi_soc(1,idir,jbshift+ib) &
599 &                          -soc_ij(2,2,idir)*cpnm_21m12
600                        psinablapsi_soc(2,idir,jbshift+ib)=psinablapsi_soc(2,idir,jbshift+ib) &
601 &                          +soc_ij(1,1,idir)*cpnm_11m22+soc_ij(1,2,idir)*cpnm_21p12
602                      end do
603                    end do ! ilmn
604                  end do ! jlmn
605                end do ! iatom
606              end if ! pawspnorb
607 
608 !            Second half of the (n,m) matrix
609              if ((ib/=jb).and.(.not.store_half_dipoles).and.(.not.iomode_etsf_mpiio)) then
610                ibshift=(ib-1)*mband
611                psinablapsi_paw(1,1:3,ibshift+jb)= psinablapsi_paw(1,1:3,jbshift+ib)
612                psinablapsi_paw(2,1:3,ibshift+jb)=-psinablapsi_paw(2,1:3,jbshift+ib)
613                if (dtset%pawspnorb==1) then
614                  psinablapsi_soc(1,1:3,ibshift+jb)= psinablapsi_soc(1,1:3,jbshift+ib)
615                  psinablapsi_soc(2,1:3,ibshift+jb)=-psinablapsi_soc(2,1:3,jbshift+ib)
616                end if
617              end if
618 
619            end do ! ib loop
620 
621            if (iomode_etsf_mpiio.and.mpi_enreg%paral_spinor==1) then
622              call timab(48,1,tsec)
623              call xmpi_sum_master(psinablapsi_paw,master,spaceComm_spinor,ierr)
624              call timab(48,2,tsec)
625            end if
626 
627          end if ! myband
628 
629 !        Write to OPT file in case of MPI-IO
630          if (iomode_etsf_mpiio.and.i_am_master_spfft) then
631            if (myband) then
632              psinablapsi=psinablapsi+psinablapsi_paw
633              if (dtset%pawspnorb==1) psinablapsi=psinablapsi+psinablapsi_soc
634            end if       
635            if (store_half_dipoles) then
636              nc_start_5=[1,1,(jb*(jb-1))/2+1,ikpt,isppol];nc_stride_5=[1,1,1,1,1]
637              nc_count_5=[0,0,0,0,0];if (myband) nc_count_5=[2,3,jb,1,1] 
638 #ifdef HAVE_NETCDF
639              NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
640 #endif
641            else
642              nc_start_6=[1,1,1,jb,ikpt,isppol];nc_stride_6=[1,1,1,1,1,1]           
643              nc_count_6=[0,0,0,0,0,0];if (myband) nc_count_6=[2,3,mband,1,1,1]
644 #ifdef HAVE_NETCDF
645              NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
646 #endif
647            end if
648          end if
649 
650        end do ! jb
651 
652        if (mkmem/=0) then
653          ibg = ibg +       my_nspinor*nband_cprj_k
654          icg = icg + npw_k*my_nspinor*nband_k
655          ikg = ikg + npw_k
656        end if
657 
658        if (cprj_paral_band) then
659          call pawcprj_free(cprj_k)
660          ABI_FREE(cprj_k)
661        end if
662        if (mkmem*nsppol/=1) then
663          call pawcprj_free(cprj_k_loc)
664          ABI_FREE(cprj_k_loc)
665        end if
666        ABI_FREE(kpg_k)
667 
668 !      Write to OPT file if not MPI-IO
669 
670 !      >>> Reduction in case of parallelism
671        if (.not.iomode_etsf_mpiio) then
672          call timab(48,1,tsec)
673          call xmpi_sum_master(psinablapsi,master,spaceComm_bandspinorfft,ierr)
674          call xmpi_sum_master(psinablapsi_paw,master,spaceComm_bandspinor,ierr)
675          call timab(48,2,tsec)
676          psinablapsi=psinablapsi+psinablapsi_paw
677          if (dtset%pawspnorb==1) then
678            call xmpi_sum_master(psinablapsi_soc,master,spaceComm_band,ierr)
679            psinablapsi=psinablapsi+psinablapsi_soc
680          end if
681        end if
682 
683 !      >>> This my kpt and I am the master node: I write the data    
684        if (.not.iomode_etsf_mpiio) then
685          if (i_am_master) then
686            if (iomode==IO_MODE_ETSF) then
687 #ifdef HAVE_NETCDF
688              if (nc_unlimited) then
689                nc_start_6=[1,1,1,ikpt,isppol,1] ; nc_count_6=[2,3,mband,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1] 
690                NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
691              else if (.not.store_half_dipoles) then
692                nc_start_6=[1,1,1,1,ikpt,isppol] ; nc_count_6=[2,3,mband,mband,1,1] ; nc_stride_6=[1,1,1,1,1,1] 
693                NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
694              else
695                nc_start_5=[1,1,1,ikpt,isppol] ; nc_count_5=[2,3,(mband*(mband+1))/2,1,1] ; nc_stride_5=[1,1,1,1,1] 
696                NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
697              end if
698 #endif
699            else
700              bsize=nband_k**2;if (store_half_dipoles) bsize=(nband_k*(nband_k+1))/2
701              write(ount)(psinablapsi(1:2,1,ib),ib=1,bsize)
702              write(ount)(psinablapsi(1:2,2,ib),ib=1,bsize)
703              write(ount)(psinablapsi(1:2,3,ib),ib=1,bsize)
704            end if
705 
706 !        >>> This my kpt and I am not the master node: I send the data    
707          else if (i_am_master_band.and.i_am_master_spfft) then
708            if (mpi_enreg%me_kpt/=master_spfftband) then
709              ABI_BUG('Problem with band communicator!')
710            end if
711            call xmpi_exch(psinablapsi,pnp_size,mpi_enreg%me_kpt,psinablapsi,master,spaceComm_kpt,etiq,ierr)
712          end if
713        end if  
714 
715 !    >>> This is not my kpt and I am the master node: I receive the data and I write    
716      elseif ((.not.iomode_etsf_mpiio).and.i_am_master) then ! mykpt
717        sender=master_spfftband
718        call xmpi_exch(psinablapsi,pnp_size,sender,psinablapsi,master,spaceComm_kpt,etiq,ierr)
719        if (iomode==IO_MODE_ETSF) then
720 #ifdef HAVE_NETCDF
721          if (nc_unlimited) then
722            nc_start_6=[1,1,1,ikpt,isppol,1] ; nc_count_6=[2,3,mband,1,1,mband] ; nc_stride_6=[1,1,1,1,1,1] 
723            NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
724          else if (.not.store_half_dipoles) then
725            nc_start_6=[1,1,1,1,ikpt,isppol] ; nc_count_6=[2,3,mband,mband,1,1] ; nc_stride_6=[1,1,1,1,1,1] 
726            NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_6,stride=nc_stride_6,count=nc_count_6))
727          else
728            nc_start_5=[1,1,1,ikpt,isppol] ; nc_count_5=[2,3,(mband*(mband+1))/2,1,1] ; nc_stride_5=[1,1,1,1,1] 
729            NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start_5,stride=nc_stride_5,count=nc_count_5))
730          end if
731 #endif
732        else
733          bsize=nband_k**2;if (store_half_dipoles) bsize=(nband_k*(nband_k+1))/2
734          write(ount)(psinablapsi(1:2,1,ib),ib=1,bsize)
735          write(ount)(psinablapsi(1:2,2,ib),ib=1,bsize)
736          write(ount)(psinablapsi(1:2,3,ib),ib=1,bsize)
737        end if
738      end if ! mykpt
739 
740      bdtot_index=bdtot_index+nband_k
741 
742 !    End loop on spin,kpt
743    end do ! ikpt
744  end do !isppol
745 
746 !Close file
747  if (i_am_master.or.(iomode_etsf_mpiio.and.i_am_master_spfft)) then
748    if (iomode==IO_MODE_ETSF) then
749 #ifdef HAVE_NETCDF
750      NCF_CHECK(nf90_close(ncid))
751 #endif
752    else
753      ierr=close_unit(ount,msg)
754      ABI_CHECK(ierr==0,"Error while closing OPT file")
755    end if
756  end if
757 
758 !Datastructures deallocations
759  ABI_FREE(psinablapsi)
760  ABI_FREE(psinablapsi_paw)
761  if (dtset%pawspnorb==1) then
762    ABI_FREE(psinablapsi_soc)
763    do iatom=1,natom
764      ABI_FREE(phisocphj(iatom)%value)
765    end do
766    ABI_FREE(phisocphj)
767  end if
768  if (.not.already_has_nabla) then
769    do itypat=1,dtset%ntypat
770      if (allocated(pawtab(itypat)%nabla_ij)) then
771        ABI_FREE(pawtab(itypat)%nabla_ij)
772        pawtab(itypat)%has_nabla=0
773      end if
774    end do
775  end if
776 
777  DBG_EXIT("COLL")
778 
779  end subroutine optics_paw

m_paw_optics/optics_paw_core [ Functions ]

[ Top ] [ m_paw_optics ] [ Functions ]

NAME

 optics_paw_core

FUNCTION

 Compute matrix elements need for X spectr. (in the PAW context) and store them in a file
  Matrix elements = <Phi_core|Nabla|Phi_j>

COPYRIGHT

 Copyright (C) 2005-2022 ABINIT group (SM,MT,NB)
 This file is distributed under the terms of the
 GNU General Public License, see ~ABINIT/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
  dimcprj(natom)=array of dimensions of array cprj (not ordered)
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  filpsp(ntypat)=name(s) of the pseudopotential file(s)
  mband=maximum number of bands
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mkmem =number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpsang =1+maximum angular momentum for nonlocal pseudopotentials
  natom=number of atoms in cell.
  nkpt=number of k points.
  nsppol=1 for unpolarized, 2 for spin-polarized
  pawang <type(pawang_type)>= PAW ANGular mesh discretization and related data
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(my_natom) <type(pawrhoij_type)>= PAW rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  znucl(ntypat)=atomic number of atom type

OUTPUT

  (only writing in a file)

SOURCE

 824  subroutine optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen0,filpsp,hdr,&
 825 &               mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,&
 826 &               pawang,pawrad,pawrhoij,pawtab,znucl)
 827 
 828 !Arguments ------------------------------------
 829 !scalars
 830  integer,intent(in) :: mband,mcprj,mkmem,mpsang,natom,nkpt,nsppol
 831  type(MPI_type),intent(in) :: mpi_enreg
 832  type(datafiles_type),intent(in) :: dtfil
 833  type(dataset_type),intent(in) :: dtset
 834  type(hdr_type),intent(inout) :: hdr
 835  type(pawang_type),intent(in) :: pawang
 836 !arrays
 837  integer,intent(in) :: atindx1(natom),dimcprj(natom)
 838  character(len=fnlen),intent(in) :: filpsp(dtset%ntypat)
 839  real(dp),intent(in) :: eigen0(mband*nkpt*nsppol),znucl(dtset%ntypat)
 840  type(pawcprj_type),target,intent(inout) :: cprj(natom,mcprj)
 841  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
 842  type(pawrhoij_type),intent(inout) :: pawrhoij(mpi_enreg%my_natom)
 843  type(pawtab_type),target,intent(inout) :: pawtab(dtset%ntypat)
 844 
 845 !Local variables-------------------------------
 846 !scalars
 847  integer,parameter :: master=0
 848  integer :: bdtot_index,cplex,etiq,iatom,ic,ibg,idir
 849  integer :: ierr,ikpt,ilmn,iln,ount,is,my_jb
 850  integer :: iorder_cprj,ispinor,isppol,istwf_k,itypat
 851  integer :: jb,jbsp,jlmn,lmn_size,lmncmax,mband_cprj,ncid,varid
 852  integer :: me,my_nspinor,nband_cprj_k,option_core,pnp_size
 853  integer :: nband_k,nphicor,ncorespinor,sender,iomode,fformopt,master_spfftband
 854  integer :: spaceComm_band,spaceComm_bandspinorfft,spaceComm_fft,spaceComm_kpt
 855  integer :: spaceComm_spinor,spaceComm_bandspinor,spaceComm_spinorfft,spaceComm_w
 856  logical :: already_has_nabla,cprj_paral_band,ex,mykpt,myband
 857  logical :: iomode_etsf_mpiio,abinitcorewf,use_spinorbit,xmlcorewf
 858  logical :: i_am_master,i_am_master_band,i_am_master_spfft
 859  character(len=fnlen) :: filecore
 860  real(dp) :: cpnm1,cpnm2
 861  character(len=500) :: msg
 862 !arrays
 863  integer :: nc_count(7),nc_start(7),nc_stride(7),tmp_shape(3)
 864  integer,allocatable :: indlmn_cor(:,:),lcor(:),ncor(:),kappacor(:)
 865  real(dp) :: tsec(2)
 866  real(dp),allocatable :: energy_cor(:),phi_cor(:,:)
 867  real(dp),allocatable :: psinablapsi(:,:,:,:,:),psinablapsi_soc(:,:,:,:,:)
 868  real(dp),pointer :: soc_ij(:,:,:)
 869  type(coeff5_type),allocatable,target :: phisocphj(:)
 870  type(pawcprj_type),pointer :: cprj_k(:,:),cprj_k_loc(:,:)
 871  type(nctkdim_t) :: ncdims(3)
 872  type(nctkarr_t) :: nctk_arrays(5)
 873 
 874 ! ************************************************************************
 875 
 876  DBG_ENTER("COLL")
 877 
 878 !Compatibility tests
 879  msg="mkmem==0 not supported anymore!"
 880  ABI_CHECK(mkmem/=0,msg)
 881 !Probably should check for spinor parallelism, because thats likely to not work correctly
 882  msg="Spinor parallelism not implemented for optics_paw_core!"
 883  ABI_CHECK(dtset%npspinor==1,msg)
 884 !Is mpi_enreg initialized?
 885  if (xmpi_paral==1) then
 886    tmp_shape = shape(mpi_enreg%proc_distrb)
 887    if (nkpt > tmp_shape(1)) then
 888      ABI_BUG('problem with proc_distrb!')
 889    end if
 890  end if
 891 
 892 !Init parallelism
 893  spaceComm_w=mpi_enreg%comm_cell
 894  if (mpi_enreg%paral_kgb==1) then
 895    spaceComm_kpt=mpi_enreg%comm_kpt
 896    spaceComm_fft=mpi_enreg%comm_fft
 897    spaceComm_band=mpi_enreg%comm_band
 898    spaceComm_spinor=mpi_enreg%comm_spinor
 899    spaceComm_bandspinor=mpi_enreg%comm_bandspinor
 900    spaceComm_spinorfft=mpi_enreg%comm_spinorfft
 901    spaceComm_bandspinorfft=mpi_enreg%comm_bandspinorfft
 902  else
 903    spaceComm_kpt=mpi_enreg%comm_kpt
 904    spaceComm_fft=xmpi_comm_self
 905    spaceComm_band=mpi_enreg%comm_band
 906    spaceComm_spinor=xmpi_comm_self
 907    spaceComm_bandspinor=spaceComm_band
 908    spaceComm_spinorfft=xmpi_comm_self
 909    spaceComm_bandspinorfft=xmpi_comm_self
 910  end if
 911  me=xmpi_comm_rank(spaceComm_w)
 912  i_am_master=(me==master)
 913  i_am_master_band=(xmpi_comm_rank(spaceComm_band)==master)
 914  i_am_master_spfft=(xmpi_comm_rank(spaceComm_spinorfft)==master)
 915  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 916 
 917 !------------------------------------------------------------------------------------------------
 918 !1- Reading of core wavefunctions
 919 !------------------------------------------------------------------------------------------------
 920 
 921 !Note: core WF is read for itypat=1
 922 !At present, impose 2-spinor simulataneously for core anf valence WF
 923  ncorespinor=merge(2,1,dtset%pawspnorb==1.or.dtset%nspinor==2)
 924  filecore=trim(filpsp(1)) ; iln=len(trim(filecore))
 925  abinitcorewf=.false. ; if (iln>3) abinitcorewf=(filecore(iln-6:iln)=='.abinit')
 926  xmlcorewf=.false. ; if (iln>3) xmlcorewf=(filecore(iln-3:iln)=='.xml')
 927  if ((.not.xmlcorewf).and.(.not.abinitcorewf)) filecore=filecore(1:iln)//'.corewf'
 928  if (abinitcorewf) filecore=filecore(1:iln-6)//'corewf.abinit'
 929  if (xmlcorewf) filecore=filecore(1:iln-3)//'corewf.xml'
 930  inquire(file=filecore,exist=ex)
 931 
 932 !Relativistic case
 933  if (ncorespinor==2) then
 934    if (ex) then
 935      !Use <filepsp>.corewf.xml or <filepsp>.corewf.abinit
 936      call pawpsp_read_corewf(energy_cor,indlmn_cor,lcor,lmncmax,ncor,nphicor,pawrad(1),phi_cor,&
 937 &                            filename=filecore,kappacor=kappacor)
 938    else
 939      !Use default name
 940      call pawpsp_read_corewf(energy_cor,indlmn_cor,lcor,lmncmax,ncor,nphicor,pawrad(1),phi_cor,&
 941 &                            kappacor=kappacor)
 942    end if
 943 
 944 !Non-relativistic case
 945  else
 946    if (ex) then
 947      !Use <filepsp>.corewf.xml or <filepsp>.corewf.abinit
 948      call pawpsp_read_corewf(energy_cor,indlmn_cor,lcor,lmncmax,ncor,nphicor,pawrad(1),phi_cor,&
 949 &                            filename=filecore)
 950    else
 951      !Use default name
 952      call pawpsp_read_corewf(energy_cor,indlmn_cor,lcor,lmncmax,ncor,nphicor,pawrad(1),phi_cor)
 953    end if
 954    ABI_MALLOC(kappacor,(nphicor))
 955    kappacor(:)=zero
 956  endif
 957 
 958 !----------------------------------------------------------------------------------
 959 !2- Computation of phipphj=<phi_i|nabla|phi_core>
 960 !----------------------------------------------------------------------------------
 961 
 962  already_has_nabla=all(pawtab(:)%has_nabla==3)
 963  if (ncorespinor==2) then
 964    already_has_nabla=all(pawtab(:)%has_nabla==4)
 965 !  Should check whether this would work with spinor parallelism
 966    call pawnabla_core_init(mpsang,dtset%ntypat,pawrad,pawtab,phi_cor,indlmn_cor,diracflag=1)
 967  else
 968    call pawnabla_core_init(mpsang,dtset%ntypat,pawrad,pawtab,phi_cor,indlmn_cor)
 969  endif
 970 
 971 !Compute spin-orbit contributions if necessary
 972  use_spinorbit=(dtset%pawspnorb==1.and.dtset%userie/=111) ! For testing purpose
 973  if (use_spinorbit) then
 974    option_core=1
 975    call pawnabla_soc_init(phisocphj,option_core,dtset%ixc,mpi_enreg%my_natom,natom,&
 976 &       dtset%nspden,dtset%ntypat,pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,&
 977 &       dtset%spnorbscl,dtset%typat,dtset%xc_denpos,znucl,&
 978 &       phi_cor=phi_cor,indlmn_cor=indlmn_cor,&
 979 &       comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 980  end if
 981 
 982 !----------------------------------------------------------------------------------
 983 !3- Opening of OPT2 file and header writing
 984 !----------------------------------------------------------------------------------
 985 
 986 !I/O mode is netCDF or Fortran
 987  iomode=merge(IO_MODE_ETSF,IO_MODE_FORTRAN_MASTER,dtset%iomode==IO_MODE_ETSF)
 988 #ifdef HAVE_NETCDF
 989  if (use_netcdf_forced) iomode=IO_MODE_ETSF
 990 #endif
 991 
 992  !(master proc only)
 993  if (i_am_master) then
 994 !  ====> NETCDF format
 995    if (iomode==IO_MODE_ETSF) then
 996      fformopt=611 
 997 #ifdef HAVE_NETCDF
 998 !    Open/create nc file
 999      NCF_CHECK(nctk_open_create(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),xmpi_comm_self))
1000 !    Write header data
1001      NCF_CHECK(hdr%ncwrite(ncid,fformopt,nc_define=.true.))
1002 !    Define additional dimensions
1003      ncdims(1)%name="number_of_core_states"
1004      ncdims(1)%value=nphicor
1005      ncdims(2)%name="number_of_core_spinor_components"
1006      ncdims(2)%value=ncorespinor
1007      ncdims(3)%name="number_of_core_spins"
1008      ncdims(3)%value=3-ncorespinor
1009      NCF_CHECK(nctk_def_dims(ncid,ncdims))
1010 !    Define additional variables
1011      nctk_arrays(1)%name="eigenvalues_core"
1012      nctk_arrays(1)%dtype="dp"
1013      nctk_arrays(1)%shape_str="number_of_core_states"
1014      nctk_arrays(2)%name="dipole_core_valence"
1015      nctk_arrays(2)%dtype="dp"
1016      nctk_arrays(2)%shape_str=&
1017 &     "complex, number_of_cartesian_directions,number_of_core_states,"// &
1018 &     "number_of_atoms,max_number_of_states,number_of_kpoints,number_of_spins"
1019      nctk_arrays(3)%name="n_quantum_number_core"
1020      nctk_arrays(3)%dtype="int"
1021      nctk_arrays(3)%shape_str="number_of_core_states"
1022      nctk_arrays(4)%name="l_quantum_number_core"
1023      nctk_arrays(4)%dtype="int"
1024      nctk_arrays(4)%shape_str="number_of_core_states"
1025      nctk_arrays(5)%name="kappa_core"
1026      nctk_arrays(5)%dtype="int"
1027      nctk_arrays(5)%shape_str="number_of_core_states"
1028      NCF_CHECK(nctk_def_arrays(ncid, nctk_arrays))
1029      NCF_CHECK(nctk_set_atomic_units(ncid, "eigenvalues_core"))
1030      NCF_CHECK(nctk_set_atomic_units(ncid, "dipole_core_valence"))
1031 !    Write core states
1032      NCF_CHECK(nctk_set_datamode(ncid))
1033      varid=nctk_idname(ncid,"eigenvalues_core")
1034      NCF_CHECK(nf90_put_var(ncid,varid,energy_cor))
1035      varid=nctk_idname(ncid,"n_quantum_number_core")
1036      NCF_CHECK(nf90_put_var(ncid,varid,ncor))
1037      varid=nctk_idname(ncid,"l_quantum_number_core")
1038      NCF_CHECK(nf90_put_var(ncid,varid,lcor))
1039      varid=nctk_idname(ncid,"kappa_core")
1040      NCF_CHECK(nf90_put_var(ncid,varid,kappacor))
1041 !    Write eigenvalues
1042      varid=nctk_idname(ncid,"eigenvalues")
1043      NCF_CHECK(nf90_put_var(ncid,varid,reshape(eigen0,[mband,nkpt,nsppol])))
1044      !Close file here because the rest has possibly to be written with collective I/O
1045      NCF_CHECK(nf90_close(ncid))
1046 #else
1047      msg = "In order to use iomode=3 and prtnabla>0 together, NetCDF support must be enabled!"
1048      ABI_ERROR(msg)
1049 #endif
1050 !  ====> Standard FORTRAN binary file format
1051    else if (iomode==IO_MODE_FORTRAN_MASTER) then
1052      fformopt=612  ! MT 12sept21: change the OPT2 Fortran file format
1053      if (2*nphicor*natom*mband>2**30) fformopt=613 ! Format for large file records
1054      if (open_file(dtfil%fnameabo_app_opt2,msg,newunit=ount,form="unformatted",status="unknown")/= 0) then
1055        ABI_ERROR(msg)
1056      end if
1057      call hdr%fort_write(ount,fformopt,ierr,rewind=.true.)
1058      write(ount)(eigen0(jb),jb=1,mband*nkpt*nsppol)
1059      write(ount) nphicor
1060      do iln=1,nphicor
1061        write(ount) ncor(iln),lcor(iln),kappacor(iln),energy_cor(iln)
1062      end do
1063    else
1064      msg = "Wrong OPT2 file format!"
1065      ABI_BUG(msg)
1066    end if ! File format
1067  end if ! master node
1068  call xmpi_bcast(iomode,master,spaceComm_w,ierr)  ! Seems mandatory; why ?
1069  call xmpi_bcast(fformopt,master,spaceComm_w,ierr)
1070  iomode_etsf_mpiio=(iomode==IO_MODE_ETSF.and.nctk_has_mpiio.and.use_netcdf_mpiio)
1071 
1072  ABI_FREE(ncor)
1073  ABI_FREE(lcor)
1074  ABI_FREE(kappacor)
1075  ABI_FREE(phi_cor)
1076  ABI_FREE(energy_cor)
1077 
1078 !----------------------------------------------------------------------------------
1079 !4- Computation of <psi_n|p_i>(<phi_i|-i.nabla|phi_core>)
1080 !----------------------------------------------------------------------------------
1081 
1082 !Prepare core-valence dipoles writing
1083 !In case of netCDF access to OPT2 file, prepare collective I/O
1084  if (iomode == IO_MODE_ETSF) then
1085    if (iomode_etsf_mpiio) then
1086      if (i_am_master_spfft) then
1087        NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),spaceComm_band))
1088        varid=nctk_idname(ncid,"dipole_core_valence")
1089        if (xmpi_comm_size(spaceComm_w)>1) then
1090          NCF_CHECK(nctk_set_collective(ncid,varid))
1091        end if
1092        NCF_CHECK(nctk_set_datamode(ncid))
1093      end if
1094    else if (i_am_master) then
1095      NCF_CHECK(nctk_open_modify(ncid,nctk_ncify(dtfil%fnameabo_app_opt2),xmpi_comm_self))
1096      varid=nctk_idname(ncid,"dipole_core_valence")
1097      if (nctk_has_mpiio.and.(.not.use_netcdf_mpiio)) then
1098        NCF_CHECK(nctk_set_collective(ncid,varid))
1099      end if
1100      NCF_CHECK(nctk_set_datamode(ncid))
1101    end if     
1102  end if
1103  if (iomode_etsf_mpiio) then
1104    !If MPI-IO, store only elements for one band
1105    ABI_MALLOC(psinablapsi,(2,3,nphicor,natom,1))
1106    if (use_spinorbit) then
1107      ABI_MALLOC(psinablapsi_soc,(2,3,nphicor,natom,1))
1108    end if
1109  else
1110    !If not, store the elements for all bands
1111    ABI_MALLOC(psinablapsi,(2,3,nphicor,natom,mband))
1112    if (use_spinorbit) then
1113      ABI_MALLOC(psinablapsi_soc,(2,3,nphicor,natom,mband))
1114    end if
1115  end if
1116  pnp_size=size(psinablapsi)
1117  
1118 !Determine if cprj datastructure is distributed over bands
1119  mband_cprj=mcprj/(my_nspinor*mkmem*nsppol)
1120  cprj_paral_band=(mband_cprj<mband)
1121 
1122 !LOOP OVER SPINS
1123  ibg=0
1124  bdtot_index=0
1125  do isppol=1,nsppol
1126 
1127 !  LOOP OVER k POINTS
1128    do ikpt=1,nkpt
1129 
1130      etiq=ikpt+(isppol-1)*nkpt
1131      nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
1132      master_spfftband=minval(mpi_enreg%proc_distrb(ikpt,1:nband_k,isppol))
1133 
1134 !    Select k-points for current proc
1135      mykpt=.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me_kpt))
1136      if (mykpt) then
1137 
1138 !      Data depending on k-point
1139        istwf_k=dtset%istwfk(ikpt)
1140        cplex=2;if (istwf_k>1) cplex=1
1141 
1142 !      Extract cprj for this k-point
1143        nband_cprj_k=nband_k;if (cprj_paral_band) nband_cprj_k=nband_k/mpi_enreg%nproc_band
1144        if (mkmem*nsppol/=1) then
1145          iorder_cprj=0
1146          ABI_MALLOC(cprj_k_loc,(natom,my_nspinor*nband_cprj_k))
1147          call pawcprj_alloc(cprj_k_loc,0,dimcprj)
1148          call pawcprj_get(atindx1,cprj_k_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isppol,&
1149 &         mband_cprj,mkmem,natom,nband_cprj_k,nband_cprj_k,my_nspinor,nsppol,dtfil%unpaw,&
1150 &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1151        else
1152          cprj_k_loc => cprj
1153        end if
1154 
1155 !      if cprj are distributed over bands, gather them (because we need to mix bands)
1156        if (cprj_paral_band) then
1157          ABI_MALLOC(cprj_k,(natom,my_nspinor*nband_k))
1158          call pawcprj_alloc(cprj_k,0,dimcprj)
1159          call pawcprj_mpi_allgather(cprj_k_loc,cprj_k,natom,my_nspinor*nband_cprj_k, &
1160 &                                   my_nspinor*mpi_enreg%bandpp,&
1161 &         dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.false.)
1162        else
1163          cprj_k => cprj_k_loc
1164        end if
1165 
1166 !      Loops over bands
1167        do jb=1,nband_k
1168          !If MPI-IO, store only ib elements for each jb ; if not, store all (ib,jb) pairs
1169          my_jb=merge(1,jb,iomode_etsf_mpiio)
1170 
1171          psinablapsi(:,:,:,:,my_jb)=zero
1172          if (use_spinorbit) psinablapsi_soc(:,:,:,:,my_jb)=zero
1173 
1174 !        Computation of <psi_n|p_i><phi_i|-i.nabla|phi_core>
1175 !        ----------------------------------------------------------------------------------
1176 
1177 !        Select bands for current proc
1178          myband=.true.
1179          if (mpi_enreg%paral_kgb==1) then
1180            myband=(mod(jb-1,mpi_enreg%nproc_band)==mpi_enreg%me_band)
1181          else if (xmpi_paral==1) then
1182            myband=(abs(mpi_enreg%proc_distrb(ikpt,jb,isppol)-mpi_enreg%me_kpt)==0)
1183          end if
1184          if (myband) then
1185 
1186            jbsp=(jb-1)*my_nspinor
1187 
1188 !          1-spinor case
1189            if (dtset%nspinor==1.and.ncorespinor==1) then
1190              jbsp=jbsp+1
1191              if (cplex==1) then ! Real WF case
1192                do iatom=1,natom
1193                  itypat=dtset%typat(iatom)
1194                  lmn_size=pawtab(itypat)%lmn_size
1195                  do jlmn=1,lmn_size
1196                    do ilmn=1,lmncmax
1197                      ic=indlmn_cor(5,ilmn)
1198                      cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
1199                      psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1200 &                        -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1201                    end do !ilmn
1202                  end do !jlmn
1203                end do !iatom
1204              else ! Complex WF case
1205                do iatom=1,natom
1206                  itypat=dtset%typat(iatom)
1207                  lmn_size=pawtab(itypat)%lmn_size
1208                  do jlmn=1,lmn_size
1209                    do ilmn=1,lmncmax
1210                      ic=indlmn_cor(5,ilmn)
1211                      cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
1212                      cpnm2=cprj_k(iatom,jbsp)%cp(2,jlmn)
1213                      psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
1214 &                        -cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1215                      psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1216 &                        -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1217                    end do !ilmn
1218                  end do !jlmn
1219                end do !iatom
1220              end if
1221 
1222 !          2-spinor case
1223            else if (dtset%nspinor==2.and.ncorespinor==2) then
1224              do ispinor=1,my_nspinor
1225                jbsp=jbsp+1
1226                do iatom=1,natom
1227                  itypat=dtset%typat(iatom)
1228                  lmn_size=pawtab(itypat)%lmn_size
1229                  do jlmn=1,lmn_size
1230                    do ilmn=1,lmncmax
1231                      is=indlmn_cor(6,ilmn)
1232                      if (modulo(jbsp,2)==modulo(is,2)) then ! Nabla is a spin-diagonal operator
1233                        ic=indlmn_cor(5,ilmn)
1234                        if (ic>0) then
1235                          cpnm1=cprj_k(iatom,jbsp)%cp(1,jlmn)
1236                          cpnm2=cprj_k(iatom,jbsp)%cp(2,jlmn)
1237                          psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
1238 &                            +cpnm1*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
1239 &                            -cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1240                          psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1241 &                            -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn) &
1242 &                            -cpnm2*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn)
1243                        end if
1244                      end if
1245                    end do ! ilmn
1246                  end do !jlmn
1247                end do !iatom
1248              end do !ispinor
1249            else
1250              msg="Core and valence WF should have the same spinor representation!"
1251              ABI_BUG(msg)
1252              !N. Brouwer initial coding: should be justified!
1253              !if (dtset%nspinor==1.and.ncorespinor==2) then  ! Average core spinors
1254              !  psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
1255              ! &    +half_sqrt2*(cpnm1*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
1256              ! &                +cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn))
1257              !  psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1258              ! &    +half_sqrt2*(cpnm2*pawtab(itypat)%nabla_im_ij(:,jlmn,ilmn) &
1259              ! &                -cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn))
1260              !else if (dtset%nspinor==2.and.ncorespinor==1) then ! Average valence spinors
1261              !  psinablapsi(1,:,ic,iatom,my_jb)=psinablapsi(1,:,ic,iatom,my_jb) &
1262              ! &    +half_sqrt2*cpnm2*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1263              !  psinablapsi(2,:,ic,iatom,my_jb)=psinablapsi(2,:,ic,iatom,my_jb) &
1264              ! &    -half_sqrt2*cpnm1*pawtab(itypat)%nabla_ij(:,jlmn,ilmn)
1265              !endif
1266            endif
1267 
1268 !          Spin-orbit coupling contribution:
1269 !             Sum_i,ss'[<psi^s_n|p_i><phi_i|1/4 Alpha^2 (Sigma^ss' X dV/dr)|phj_core^s'>]
1270            if (use_spinorbit) then
1271 !            Add: Sum_i,ss'[<Psi^s_n|p_i> (Sigma^ss' X g_ij_core^s']
1272 !             where: g_ij_core^s = <Phi_i| 1/4 Alpha^2 dV/dr vec(r)/r) |Phj_core^s>
1273 !            Note that:
1274 !             if phi_cor_jlmn_cor(:) is up:
1275 !              phisocphj(:)%value(re:im,1,idir,ilmn,jlmn_cor) is (Sigma^up-up X g_ij_core^up)
1276 !              phisocphj(:)%value(re:im,2,idir,ilmn,jlmn_cor) is (Sigma^dn-up X g_ij_core^up)
1277 !             if phi_cor_jlmn_cor(:) is down:
1278 !              phisocphj(:)%value(re:im,1,idir,ilmn,jlmn_cor) is (Sigma^up-dn X g_ij_core^dn)
1279 !              phisocphj(:)%value(re:im,2,idir,ilmn,jlmn_cor) is (Sigma^dn-dn X g_ij_core^dn)
1280 !            Not compatible with parallelization over spinors
1281              jbsp=1+(jb-1)*dtset%nspinor
1282              do iatom=1,natom
1283                itypat=dtset%typat(iatom)
1284                lmn_size=pawtab(itypat)%lmn_size
1285                do jlmn=1,lmn_size
1286                  do ilmn=1,lmncmax
1287                    ic=indlmn_cor(5,ilmn)
1288                    if (ic>0) then
1289                      soc_ij => phisocphj(iatom)%value(:,:,:,jlmn,ilmn)
1290                      !Contribution from real part of <Psi^s_n|p_i>
1291                      cpnm1=cprj_k(iatom,jbsp  )%cp(1,jlmn)
1292                      cpnm2=cprj_k(iatom,jbsp+1)%cp(1,jlmn)
1293                      do idir=1,3
1294                        psinablapsi_soc(1,idir,ic,iatom,my_jb)=psinablapsi_soc(1,idir,ic,iatom,my_jb) &
1295 &                             +soc_ij(1,1,idir)*cpnm1+soc_ij(1,2,idir)*cpnm2
1296                        psinablapsi_soc(2,idir,ic,iatom,my_jb)=psinablapsi_soc(2,idir,ic,iatom,my_jb) &
1297 &                             +soc_ij(2,1,idir)*cpnm1+soc_ij(2,2,idir)*cpnm2
1298                      end do
1299                      !Contribution from imaginary part of <Psi^s_n|p_i>
1300                      cpnm1=cprj_k(iatom,jbsp  )%cp(2,jlmn)
1301                      cpnm2=cprj_k(iatom,jbsp+1)%cp(2,jlmn)
1302                      do idir=1,3
1303                        psinablapsi_soc(1,idir,ic,iatom,my_jb)=psinablapsi_soc(1,idir,ic,iatom,my_jb) &
1304 &                             +soc_ij(2,1,idir)*cpnm1+soc_ij(2,2,idir)*cpnm2
1305                        psinablapsi_soc(2,idir,ic,iatom,my_jb)=psinablapsi_soc(2,idir,ic,iatom,my_jb) &
1306 &                             -soc_ij(1,1,idir)*cpnm1-soc_ij(1,2,idir)*cpnm2
1307                      end do
1308                    end if
1309                  end do ! ilmn
1310                end do ! jlmn
1311              end do ! iatom
1312            end if ! use_spinorbit
1313 
1314            if (iomode_etsf_mpiio.and.mpi_enreg%paral_spinor==1) then
1315              call timab(48,1,tsec)
1316              call xmpi_sum_master(psinablapsi,master,spaceComm_spinor,ierr)
1317              call timab(48,2,tsec)
1318            end if
1319 
1320          end if ! myband
1321 
1322 !        Write to OPT2 file in case of MPI-IO
1323          if (iomode_etsf_mpiio.and.i_am_master_spfft) then
1324            nc_start=[1,1,1,1,jb,ikpt,isppol];nc_stride=[1,1,1,1,1,1,1] 
1325            if (myband) then
1326              if (use_spinorbit) psinablapsi=psinablapsi+psinablapsi_soc
1327              nc_count=[2,3,nphicor,natom,1,1,1]
1328            else
1329              nc_count=[0,0,0,0,0,0,0]
1330            end if
1331 #ifdef HAVE_NETCDF
1332            NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start,stride=nc_stride,count=nc_count))
1333 #endif
1334          end if
1335 
1336        end do ! jb
1337        
1338        if (mkmem/=0) then
1339          ibg = ibg +  my_nspinor*nband_cprj_k
1340        end if
1341 
1342        if (cprj_paral_band) then
1343          call pawcprj_free(cprj_k)
1344          ABI_FREE(cprj_k)
1345        end if
1346        if (mkmem*nsppol/=1) then
1347          call pawcprj_free(cprj_k_loc)
1348          ABI_FREE(cprj_k_loc)
1349        end if
1350 
1351 !      Write to OPT2 file if not MPI-IO
1352 
1353 !      >>> Reduction in case of parallelism
1354        if (.not.iomode_etsf_mpiio) then
1355          call timab(48,1,tsec)
1356          call xmpi_sum_master(psinablapsi,master,spaceComm_bandspinor,ierr)
1357          call timab(48,2,tsec)
1358          if (use_spinorbit) then
1359            call xmpi_sum_master(psinablapsi_soc,master,spaceComm_band,ierr)
1360            psinablapsi=psinablapsi+psinablapsi_soc
1361          end if
1362        end if
1363 
1364 !      >>> This my kpt and I am the master node: I write the data    
1365        if (.not.iomode_etsf_mpiio) then
1366          if (i_am_master) then
1367            if (iomode==IO_MODE_ETSF) then
1368              nc_start=[1,1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1,1] 
1369              nc_count=[2,3,nphicor,natom,mband,1,1]
1370 #ifdef HAVE_NETCDF
1371              NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start,stride=nc_stride,count=nc_count))
1372 #endif
1373            else
1374              if (fformopt==612) then ! New OPT2 file format
1375                write(ount) (((psinablapsi(1:2,1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1376                write(ount) (((psinablapsi(1:2,2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1377                write(ount) (((psinablapsi(1:2,3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1378              else if (fformopt==613) then ! Large OPT2 file format
1379                do jb=1,nband_k
1380                  write(ount) ((psinablapsi(1:2,1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1381                  write(ount) ((psinablapsi(1:2,2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1382                  write(ount) ((psinablapsi(1:2,3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1383                end do
1384              else ! Old OPT2 file format
1385                !The old writing was not efficient (indexes order is bad)
1386                do iatom=1,natom
1387                  write(ount) ((psinablapsi(1:2,1,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1388                  write(ount) ((psinablapsi(1:2,2,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1389                  write(ount) ((psinablapsi(1:2,3,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1390                end do
1391              end if
1392            end if
1393 
1394 !        >>> This my kpt and I am not the master node: I send the data    
1395          else if (i_am_master_band.and.i_am_master_spfft) then
1396            if (mpi_enreg%me_kpt/=master_spfftband) then
1397              ABI_BUG('Problem with band communicator!')
1398            end if
1399            call xmpi_exch(psinablapsi,pnp_size,mpi_enreg%me_kpt,psinablapsi,master,spaceComm_kpt,etiq,ierr)
1400          end if
1401        end if  
1402 
1403 !    >>> This is not my kpt and I am the master node: I receive the data and I write    
1404      elseif ((.not.iomode_etsf_mpiio).and.i_am_master) then ! mykpt
1405        sender=master_spfftband
1406        call xmpi_exch(psinablapsi,pnp_size,sender,psinablapsi,master,spaceComm_kpt,etiq,ierr)
1407        if (iomode==IO_MODE_ETSF) then
1408          nc_start=[1,1,1,1,1,ikpt,isppol];nc_stride=[1,1,1,1,1,1,1] 
1409          nc_count=[2,3,nphicor,natom,mband,1,1]
1410 #ifdef HAVE_NETCDF
1411          NCF_CHECK(nf90_put_var(ncid,varid,psinablapsi,start=nc_start,stride=nc_stride,count=nc_count))
1412 #endif
1413        else
1414          if (fformopt==612) then ! New OPT2 file format
1415            write(ount) (((psinablapsi(1:2,1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1416            write(ount) (((psinablapsi(1:2,2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1417            write(ount) (((psinablapsi(1:2,3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom),jb=1,nband_k)
1418          else if (fformopt==613) then ! Large OPT2 file format
1419            do jb=1,nband_k
1420              write(ount) ((psinablapsi(1:2,1,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1421              write(ount) ((psinablapsi(1:2,2,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1422              write(ount) ((psinablapsi(1:2,3,ic,iatom,jb),ic=1,nphicor),iatom=1,natom)
1423            end do
1424          else ! Old OPT2 file format
1425            !The old writing was not efficient (indexes order is bad)
1426            do iatom=1,natom
1427              write(ount) ((psinablapsi(1:2,1,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1428              write(ount) ((psinablapsi(1:2,2,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1429              write(ount) ((psinablapsi(1:2,3,ic,iatom,jb),jb=1,nband_k),ic=1,nphicor)
1430            end do
1431          end if
1432        end if
1433      end if ! mykpt
1434 
1435      bdtot_index=bdtot_index+nband_k
1436 
1437 !    End loop on spin,kpt
1438    end do ! ikpt
1439  end do !isppol
1440 
1441 !Close file
1442  if (i_am_master.or.(iomode_etsf_mpiio.and.i_am_master_spfft)) then
1443    if (iomode==IO_MODE_ETSF) then
1444 #ifdef HAVE_NETCDF
1445      NCF_CHECK(nf90_close(ncid))
1446 #endif
1447    else
1448      ierr=close_unit(ount,msg)
1449      ABI_CHECK(ierr==0,"Error while closing OPT2 file")
1450    end if
1451  end if
1452 
1453 !Datastructures deallocations
1454  ABI_FREE(indlmn_cor)
1455  ABI_FREE(psinablapsi)
1456  if (use_spinorbit) then
1457    ABI_FREE(psinablapsi_soc)
1458    do iatom=1,natom
1459      ABI_FREE(phisocphj(iatom)%value)
1460    end do
1461    ABI_FREE(phisocphj)
1462  end if
1463  if (.not.already_has_nabla) then
1464    do itypat=1,dtset%ntypat
1465      if (allocated(pawtab(itypat)%nabla_ij)) then
1466        ABI_FREE(pawtab(itypat)%nabla_ij)
1467        pawtab(itypat)%has_nabla=0
1468      end if
1469      if (allocated(pawtab(itypat)%nabla_im_ij)) then
1470        ABI_FREE(pawtab(itypat)%nabla_im_ij)
1471      end if
1472    end do
1473  end if
1474 
1475  DBG_EXIT("COLL")
1476 
1477  end subroutine optics_paw_core

m_paw_optics/pawnabla_soc_init [ Functions ]

[ Top ] [ m_paw_optics ] [ Functions ]

NAME

 pawnabla_soc_init

FUNCTION

 Compute the PAW SOC contribution(s) to the momentum PAW matrix elements,
  i.e. <Phi_i|1/4 Alpha^2 dV/dr (Sigma X vec(r)/r) |Phi_j>
   where:
    {Phi_i}= AE partial waves
    Alpha = inverse of fine structure constant
    Sigma^alpha,beta= Pauli matrices
    X = cross product

 There are 2 typical uses:
   - Valence-valence terms: Phi_i and Phi_j are PAW AE partial waves (unpolarized)
   - Core-valence terms: Phi_i is are AE partial waves and Phi_j are core AE wave-functions (spinors)

 In practice we compute:
  (Sigma^up-up X g_ij) and (Sigma^up-dn X g_ij)    (X = vector cross product)
  (Sigma^dn-up X g_ij) and (Sigma^dn-dn X g_ij)
   where:
    g_ij= 1/4 Alpha^2 Int_[Phi_i(r)/r Phi_j(r)/r dV(r)/dr r^2 dr] . Gvec_ij
        and Gvec_ij= Int[S_limi S_ljmj vec(r)/r dOmega] (Gaunt coefficients)

COPYRIGHT

 Copyright (C) 2021-2022 ABINIT group (NBrouwer,MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~ABINIT/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

  ixc= choice of exchange-correlation scheme (see above, and below)
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  nspden=number of spin-density components
  ntypat=number of types of atoms in unit cell.
  option_core=Type of calculation: 0=valence-valence, 1=core-valence
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
  spnorbscl=scaling factor for spin-orbit coupling
  typat(natom) =Type of each atoms
  xc_denpos= lowest allowe density (usually for the computation of the XC functionals)
  znucl(ntypat)=gives the nuclear charge for all types of atoms
  [phi_cor(mesh_size,nphicor)]=--optional-- core wave-functions for the current type of atoms;
                               only needed when option_core=1
  [indlmn_cor(6,nlmn_core)]=--optional-- array giving l,m,n,lm,ln,s for i=lmn, for the core wave functions;
                            only needed when option_core=1
  [mpi_atmtab(:)]=--optional-- indexes of the atoms treated by current proc
  [comm_atom]=--optional-- MPI communicator over atoms

OUTPUT

     phisocphj(iat)%value(1,1,idir,ilmn,jlmn) is real part of (Sigma^up-up X g_ij)
     phisocphj(iat)%value(2,1,idir,ilmn,jlmn) is imaginary part of (Sigma^up-up X g_ij)
     phisocphj(iat)%value(1,2,idir,ilmn,jlmn) is real part of (Sigma^up-dn X g_ij)
     phisocphj(iat)%value(2,2,idir,ilmn,jlmn) is imaginary part of (Sigma^up-dn X g_ij)
   If option_core==1 and nspinor_cor==2 (core-valence with spinorial core WF):
     phisocphj(iat)%value(1,1,idir,ilmn,2*jlmn-1) is real part of (Sigma^up-up X g_ij^up)
     phisocphj(iat)%value(2,1,idir,ilmn,2*jlmn-1) is imaginary part of (Sigma^up-up X g_ij^up)
     phisocphj(iat)%value(1,2,idir,ilmn,2*jlmn-1) is real part of (Sigma^dn-up X g_ij^dn)
     phisocphj(iat)%value(2,2,idir,ilmn,2*jlmn-1) is imaginary part of (Sigma^dn-up X g_ij^dn)
     phisocphj(iat)%value(1,1,idir,ilmn,2*jlmn  ) is real part of (Sigma^up-dn X g_ij^up)
     phisocphj(iat)%value(2,1,idir,ilmn,2*jlmn  ) is imaginary part of (Sigma^up-dn X g_ij^up)
     phisocphj(iat)%value(1,2,idir,ilmn,2*jlmn  ) is real part of (Sigma^dn-dn X g_ij^dn)
     phisocphj(iat)%value(2,2,idir,ilmn,2*jlmn  ) is imaginary part of (Sigma^dn-dn X g_ij^dn)
 (idir=cartesian direction)

SIDE EFFECTS

NOTES

 If Phi_j is not polarized,
    (Sigma^dn-dn X g_ij)=-(Sigma^up-up X g_ij)
    (Sigma^dn-up X g_ij)= (Sigma^up-dn X g_ij)^*
    So, we store only 2 components.
 If Phi_j is polarized, the spin component is included in the last dimension of
   phisocphj(iat)%value, i.e. lmn_size_cor=2*lmn_size

SOURCE

1839  subroutine pawnabla_soc_init(phisocphj,option_core,ixc,my_natom,natom,nspden,ntypat,pawang, &
1840 &           pawrad,pawrhoij,pawtab,pawxcdev,spnorbscl,typat,xc_denpos,znucl, &
1841 &           phi_cor,indlmn_cor,mpi_atmtab,comm_atom) ! Optional arguments
1842 
1843 !Arguments ------------------------------------
1844 !scalars
1845  integer,intent(in) :: ixc,my_natom,natom,nspden,ntypat,option_core,pawxcdev
1846  integer,optional,intent(in) :: comm_atom
1847  real(dp),intent(in) :: spnorbscl,xc_denpos
1848  type(pawang_type),intent(in) :: pawang
1849 !arrays
1850  integer,intent(in),target,optional :: indlmn_cor(:,:)
1851  integer,intent(in) :: typat(natom)
1852  integer,optional,target,intent(in) :: mpi_atmtab(:)
1853  real(dp),intent(in),target,optional :: phi_cor(:,:)
1854  real(dp),intent(in) :: znucl(ntypat)
1855  type(coeff5_type),allocatable,target,intent(inout) :: phisocphj(:)
1856  type(pawrad_type),target,intent(in) :: pawrad(ntypat)
1857  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
1858  type(pawtab_type),target,intent(in) :: pawtab(ntypat)
1859 
1860 !Local variables-------------------------------
1861 !scalars
1862  real(dp),parameter :: one_over_fourpi   = one/sqrt(four_pi)
1863  real(dp),parameter :: sqr_fourpi_over_3 = sqrt(four_pi/3)
1864  real(dp),parameter :: QuarterFineStruct2=(half/InvFineStruct)**2
1865  real(dp),parameter :: hyb_mixing_ = 0.0_dp   ! Fake value to be updated in the future
1866  integer :: iatom,iatom_tot,itypat,ii,jj,ierr,ipts,ignt,sgnkappa
1867  integer :: idum,option,usenhat,usekden,usecore,xclevel,nkxc,my_comm_atom
1868  integer :: mesh_size,mesh_size_cor,lmn_size,lmn2_size,lmn_size_j,lmn_size_cor
1869  integer :: lm_size,ln_size,ln_size_j,ln_size_cor,nspinor_cor
1870  integer :: ilmn,ilm,iln,jl,jm,jm_re,jm_im,jlmn,jlm,jlm_re,jlm_im,jln,js,klm_re,klm_im
1871  logical :: my_atmtab_allocated,paral_atom
1872  real(dp) :: avg,cgc,compch_sph_dum,eexc_dum,eexcdc_dum,jmj
1873  real(dp) :: fact_re,fact_im,gx_re,gx_im,gy_re,gy_im,gz_re,gz_im,if3
1874  character(len=500) :: msg
1875 !arrays
1876  integer,pointer :: my_atmtab(:),indlmn(:,:),indlmn_j(:,:)
1877  logical,allocatable :: lmselect(:)
1878  real(dp) :: nhat_dum(1,1,1),trho_dum(1,1,1),kxc_dum(1,1,1),k3xc_dum(1,1,1)
1879  real(dp),allocatable :: rho1(:,:,:),tau1(:,:,:),rhosph(:),vhartree(:),vxc(:,:,:)
1880  real(dp),allocatable :: intf3(:,:),potsph(:),dVdr(:),func(:)
1881  real(dp),pointer :: phi_j(:,:),soc_ij(:,:,:,:,:)
1882  type(pawrad_type),pointer :: pawrd
1883  type(pawtab_type),pointer :: pawtb
1884 
1885 ! ************************************************************************
1886 
1887 !Some checks in case of core-valence (option_core=1)
1888  if (option_core/=0.and.option_core/=1) then
1889    msg='Wrong option_core value!'
1890    ABI_BUG(msg)
1891  end if
1892  if (option_core==1) then
1893 !  Check if we have the optional arguments
1894    if ((.not.present(phi_cor)).or.(.not.present(indlmn_cor))) then
1895      msg='For core-valence calculation, need phi_cor and indlmn_cor!'
1896      ABI_BUG(msg)
1897    end if
1898 !  Check if we have relativistic core wave functions
1899    if (size(indlmn_cor,1)<8) then
1900      write(msg,'(a)') 'Wrong 1st dim. of indlmn_cor in pawnabla_soc_init (need spinors)!'
1901      ABI_BUG(msg)
1902    end if
1903  endif
1904 
1905 !Some useful variables
1906  usekden=pawxc_get_usekden(ixc)
1907  usecore=1 ; nkxc=0 ; usenhat=0
1908  xclevel=pawxc_get_xclevel(ixc)
1909  if (option_core==1) then
1910    mesh_size_cor=size(phi_cor,1)
1911    lmn_size_cor=size(indlmn_cor,2) !Includes spinors
1912    ln_size_cor=size(phi_cor,2)
1913    nspinor_cor=maxval(indlmn_cor(6,1:lmn_size_cor))
1914  end if
1915 
1916 !Prepare output arrays
1917  if (allocated(phisocphj)) then
1918    do iatom=1,natom
1919       if (allocated(phisocphj(iatom)%value)) then
1920         ABI_FREE(phisocphj(iatom)%value)
1921       end if
1922    end do
1923    ABI_FREE(phisocphj)
1924  end if
1925  ABI_MALLOC(phisocphj,(natom))
1926  do iatom=1,natom
1927    lmn_size=pawtab(typat(iatom))%lmn_size
1928    if (option_core==0) then
1929      ABI_MALLOC(phisocphj(iatom)%value,(2,2,3,lmn_size,lmn_size))
1930    else
1931      !lmn_size_cor is double because it contains the spin component
1932      ABI_MALLOC(phisocphj(iatom)%value,(2,2,3,lmn_size,lmn_size_cor))
1933    end if
1934    phisocphj(iatom)%value=zero
1935  end do
1936 
1937 !Set up parallelism over atoms
1938  paral_atom=(present(comm_atom).and.(my_natom/=natom))
1939  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1940  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1941  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
1942 &                   my_natom_ref=my_natom)
1943 
1944 !-------------------------------------------------------------------
1945 !Loop over atoms
1946 
1947  do iatom=1,my_natom
1948 
1949 !  Atom-dependent data
1950    itypat=typat(iatom)
1951    pawrd => pawrad(itypat)
1952    pawtb  => pawtab(itypat)
1953    lmn_size=pawtb%lmn_size
1954    lmn2_size=pawtb%lmn2_size
1955    lm_size=pawtb%lcut_size**2
1956    ln_size=pawtb%basis_size
1957    mesh_size=pawtb%mesh_size
1958    indlmn => pawtb%indlmn
1959    ABI_MALLOC(lmselect,(lm_size))
1960    lmselect(:)=.true.
1961 
1962 !  Distinguish valence-valence and core-valence cases
1963    if (option_core==0) then
1964      ln_size_j=ln_size
1965      lmn_size_j=lmn_size
1966      indlmn_j => pawtb%indlmn
1967      phi_j => pawtb%phi
1968    else
1969      ln_size_j=ln_size_cor
1970      lmn_size_j=lmn_size_cor
1971      indlmn_j => indlmn_cor
1972      phi_j => phi_cor
1973      if (mesh_size_cor<mesh_size) then
1974        msg='mesh_size and mesh_size_cor not compatible!'
1975        ABI_BUG(msg)
1976      end if
1977    endif
1978 
1979 !  Manage parallelism
1980    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
1981    soc_ij => phisocphj(iatom_tot)%value(:,:,:,:,:)
1982 
1983 !-------------------------------------------------------------------
1984 !Compute all-electron density
1985 
1986    ABI_MALLOC(rho1,(mesh_size,lm_size,nspden))
1987    option=2 ; rho1=zero
1988    call pawdensities(compch_sph_dum,1,iatom_tot,lmselect,lmselect,lm_size,&
1989 &       nhat_dum,nspden,-1,0,option,-1,0,pawang,0,pawrd,&
1990 &       pawrhoij(iatom),pawtb,rho1,trho_dum)
1991    if (usekden==1) then
1992      ABI_MALLOC(tau1,(mesh_size,lm_size,nspden))
1993      tau1=zero
1994      call pawkindensities(1,lmselect,lm_size,nspden,-1,option,-1,&
1995 &         pawang,pawrd,pawrhoij(iatom),pawtb,tau1,trho_dum)
1996    end if
1997 
1998 !-------------------------------------------------------------------
1999 !Compute spherical potential and compute its first derivative dV/dr
2000 
2001    ABI_MALLOC(potsph,(mesh_size))
2002    ABI_MALLOC(dVdr,(mesh_size))
2003    potsph=zero ; dVdr=zero
2004 
2005 !  Compute XC potential
2006    option=1
2007    if (pawxcdev/=0) then
2008      ABI_MALLOC(vxc,(mesh_size,lm_size,nspden))
2009      vxc=zero
2010      call pawxcm(pawtb%coredens,eexc_dum,eexcdc_dum,idum,hyb_mixing_,ixc,kxc_dum,lm_size,&
2011 &         lmselect,nhat_dum,nkxc,.false.,mesh_size,nspden,option,pawang,pawrd,&
2012 &         pawxcdev,rho1,usecore,usenhat,vxc,xclevel,xc_denpos)
2013      potsph(1:mesh_size)=half*(vxc(1:mesh_size,1,1)+vxc(1:mesh_size,1,nspden))
2014    else
2015      ABI_MALLOC(vxc,(mesh_size,pawang%angl_size,nspden))
2016      vxc=zero
2017      call pawxc(pawtb%coredens,eexc_dum,eexcdc_dum,hyb_mixing_,ixc,kxc_dum,k3xc_dum,&
2018 &         lm_size,lmselect,nhat_dum,nkxc,nkxc,.false.,mesh_size,nspden,option,pawang,&
2019 &         pawrd,rho1,usecore,usenhat,vxc,xclevel,xc_denpos,&
2020 &         coretau=pawtb%coretau,taur=tau1)
2021      potsph(1:mesh_size)=zero
2022      do ipts=1,pawang%angl_size
2023        potsph(1:mesh_size)=potsph(1:mesh_size) &
2024 &        +half*(vxc(1:mesh_size,ipts,1)+vxc(1:mesh_size,ipts,nspden))*pawang%angwgth(ipts)
2025      end do
2026      potsph(1:mesh_size)=sqrt(four_pi)*potsph(1:mesh_size)
2027    end if
2028 
2029 !  Compute Hartree potentialHalfFineStruct2
2030    ABI_MALLOC(vhartree,(mesh_size))
2031    ABI_MALLOC(rhosph,(mesh_size))
2032    vhartree=zero ; rhosph=zero
2033    rhosph(1:mesh_size)=rho1(1:mesh_size,1,1)
2034    if (usecore==1) rhosph(1:mesh_size)=rhosph(1:mesh_size)+sqrt(four_pi)*pawtb%coredens(1:mesh_size)
2035    rhosph(1:mesh_size)=rhosph(1:mesh_size)*four_pi*pawrd%rad(1:mesh_size)**2
2036    call poisson(rhosph,0,pawrd,vhartree)
2037    vhartree(2:mesh_size)=(vhartree(2:mesh_size)-sqrt(four_pi)*znucl(itypat))/pawrd%rad(2:mesh_size)
2038    call pawrad_deducer0(vhartree,mesh_size,pawrd)
2039    potsph(1:mesh_size)=potsph(1:mesh_size)+vhartree(1:mesh_size)
2040 
2041 !  Apply angular scaling factor
2042    potsph(1:mesh_size)=one_over_fourpi*potsph(1:mesh_size)
2043 
2044 !  Compute 1st derivative of potential
2045    call nderiv_gen(dVdr,potsph,pawrd)
2046 
2047 !  Multiply by relativistic factor
2048    dVdr(1:mesh_size)=dVdr(1:mesh_size)*(one/(one-potsph(1:mesh_size)*half/InvFineStruct**2)**2)
2049 
2050    ABI_FREE(vxc)
2051    ABI_FREE(vhartree)
2052    ABI_FREE(potsph)
2053    ABI_FREE(rhosph)
2054    ABI_FREE(lmselect)
2055    ABI_FREE(rho1)
2056    if (usekden==1) then
2057      ABI_FREE(tau1)
2058    end if
2059 
2060 !-------------------------------------------------------------------
2061 !Compute radial and angular contributions
2062 
2063    ABI_MALLOC(intf3,(ln_size,ln_size_j))
2064    intf3=zero
2065 
2066 !  >>>> Calculate f_3= alpha^2/4 int[dr ui*(r) uj(r) dV/dr]
2067    ABI_MALLOC(func,(mesh_size))
2068    do jln=1,ln_size_j
2069      do iln=1,ln_size
2070        func(1:mesh_size)=dVdr(1:mesh_size)*pawtb%phi(1:mesh_size,iln)*phi_j(1:mesh_size,jln)
2071        call simp_gen(intf3(iln,jln),func,pawrd)
2072      end do
2073    end do
2074    intf3(:,:)=QuarterFineStruct2*spnorbscl*intf3(:,:)
2075    ABI_FREE(func)
2076 
2077 !  Loop over initial states (valence or core, according to option_core)
2078    do jlmn=1,lmn_size_j
2079      jl=indlmn_j(1,jlmn)
2080      jm=indlmn_j(2,jlmn)
2081      jlm=indlmn_j(4,jlmn)
2082      jln=indlmn_j(5,jlmn)
2083 
2084 !    In case of spinorial core wave function, we have to handle imaginary
2085 !      spherical harmonics as linear combination of real spherical harmonics
2086 !    See Brouwer et al, CPC 266, 108029 (2021), equation (21)
2087      if (option_core==1) then
2088        jm_re= abs(jm)
2089        jm_im=-abs(jm)
2090        jlm_re=jl*(jl+1)+jm_re+1
2091        jlm_im=jl*(jl+1)+jm_im+1
2092 !      Calculate spinor dependent coefficients
2093        sgnkappa=indlmn_j(3,jlmn)   !sign of kappa
2094        jmj=half*indlmn_j(8,jlmn)   !2mj is stored in indlmn_cor
2095        js=indlmn_cor(6,jlmn)       !1 is up, 2 is down
2096        if (sgnkappa==1) then
2097          if(js==1) then
2098            cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
2099          else
2100            cgc=-sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
2101          endif
2102        else
2103          if(js==1) then
2104            cgc= sqrt((dble(jl)+jmj+half)/dble(2*jl+1))
2105          else
2106            cgc= sqrt((dble(jl)-jmj+half)/dble(2*jl+1))
2107          endif
2108        endif
2109 
2110 !      Calculate factors to convert from complex to real sph. harm.
2111        if (jm<0) then
2112          fact_re=sqr_fourpi_over_3*half_sqrt2*cgc
2113          fact_im=-fact_re
2114        else if (jm>0) then
2115          fact_re=sqr_fourpi_over_3*half_sqrt2*cgc*(-1)**jm
2116          fact_im=fact_re
2117        else
2118          fact_re=sqr_fourpi_over_3
2119          fact_im=0
2120        end if
2121      else ! valence-valence case (real)
2122        js=1
2123        jlm_re=jlm ; jlm_im=jlm
2124        fact_re=sqr_fourpi_over_3
2125        fact_im=0
2126      end if
2127 
2128 !    Loop over final states
2129      do ilmn=1,lmn_size
2130        ilm=indlmn(4,ilmn)
2131        iln=indlmn(5,ilmn)
2132 
2133 !      >>>> Calculate g_ij=(g_x,g_y,g_z) = sqrt(4pi/3) int dOmega Ylm Ylm' S1-1,0,1
2134 !              using real Gaunt coefficients
2135        gx_re=zero;gy_re=zero;gz_re=zero         
2136        gx_im=zero;gy_im=zero;gz_im=zero
2137        if3=zero
2138 
2139 !      jl was set as a flag for invalid combinations
2140 !        i.e. m=-(l+1) or m=(l+1)
2141 !      In these cases, cgc=0 ; so gx=gy=gz=0 
2142        if (jl/=-1) then
2143          if3=intf3(iln,jln)
2144          klm_re=merge((jlm_re*(jlm_re-1))/2+ilm,(ilm*(ilm-1))/2+jlm_re,ilm<=jlm_re)
2145 
2146 !        Real parts
2147 !        M=-1
2148          ignt=pawang%gntselect(2,klm_re) !get index for L=1 M =-1 ilm jlm_re
2149          if (ignt/=0) gy_re=fact_re*pawang%realgnt(ignt)
2150 !        M=0
2151          ignt=pawang%gntselect(3,klm_re) !get index for L=1 M = 0 ilm jlm_re
2152          if (ignt/=0) gz_re=fact_re*pawang%realgnt(ignt)
2153 !        M=1
2154          ignt=pawang%gntselect(4,klm_re) !get index for L=1 M = 1 ilm jlm_re
2155          if (ignt/=0) gx_re=fact_re*pawang%realgnt(ignt)
2156 
2157 !        Imaginary parts
2158          if (option_core==1) then
2159            klm_im=merge((jlm_im*(jlm_im-1))/2+ilm,(ilm*(ilm-1))/2+jlm_im,ilm<=jlm_im)
2160 !          M=-1
2161            ignt=pawang%gntselect(2,klm_im) !get index for L=1 M =-1 ilm jlm_im
2162            if (ignt/=0) gy_im=fact_im*pawang%realgnt(ignt)
2163 !          M=0
2164            ignt=pawang%gntselect(3,klm_im) !get index for L=1 M = 0 ilm jlm_im
2165            if (ignt/=0) gz_im=fact_im*pawang%realgnt(ignt)
2166 !          M=1
2167            ignt=pawang%gntselect(4,klm_im) !get index for L=1 M = 1 ilm jlm_im
2168            if (ignt/=0) gx_im=fact_im*pawang%realgnt(ignt)
2169          end if
2170        end if
2171 
2172 !      >>>> Calculate Sigma X g_ij
2173 
2174        if (option_core==0.or.js==1) then
2175          !(Sigma^up-up X gij)_x = -gy*f_3
2176          soc_ij(1,1,1,ilmn,jlmn)=-if3*gy_re           ! real part
2177          soc_ij(2,1,1,ilmn,jlmn)=-if3*gy_im           ! imag part
2178          !(Sigma^up-up X gij)_y = gx*f_3
2179          soc_ij(1,1,2,ilmn,jlmn)= if3*gx_re           ! real part
2180          soc_ij(2,1,2,ilmn,jlmn)= if3*gx_im           ! imag part
2181          !(Sigma^up-up X gij)_z = 0
2182          soc_ij(1,1,3,ilmn,jlmn)= zero                ! real part
2183          soc_ij(2,1,3,ilmn,jlmn)= zero                ! imag part
2184 
2185          !(Sigma^dn-up X gij)_x = i.gz*f_3
2186          soc_ij(1,2,1,ilmn,jlmn)=-if3*gz_im           ! real part
2187          soc_ij(2,2,1,ilmn,jlmn)= if3*gz_re           ! imag part
2188          !(Sigma^dn-up X gij)_y = -gz*f_3
2189          soc_ij(1,2,2,ilmn,jlmn)=-if3*gz_re           ! real part
2190          soc_ij(2,2,2,ilmn,jlmn)=-if3*gz_im           ! imag part
2191          !(Sigma^dn-up X gij)_z = (gy-i.gx)*f_3
2192          soc_ij(1,2,3,ilmn,jlmn)= if3*(gy_re+gx_im)   ! real part
2193          soc_ij(2,2,3,ilmn,jlmn)= if3*(gy_im-gx_re)   ! imag part
2194 
2195        else if (option_core==1.and.js==2) then
2196          !(Sigma^up-dn X gij^dn)_x = -i.gz^dn*f_3
2197          soc_ij(1,1,1,ilmn,jlmn)= if3*gz_im           ! real part
2198          soc_ij(2,1,1,ilmn,jlmn)=-if3*gz_re           ! imag part
2199          !(Sigma^up-dn X gij^dn)_y = -gz^dn*f_3
2200          soc_ij(1,1,2,ilmn,jlmn)=-if3*gz_re           ! real part
2201          soc_ij(2,1,2,ilmn,jlmn)=-if3*gz_im           ! imag part
2202          !(Sigma^up-dn X gij^dn)_z = (gy^dn+i.gx^dn)*f_3
2203          soc_ij(1,1,3,ilmn,jlmn)= if3*(gy_re-gx_im)   ! real part
2204          soc_ij(2,1,3,ilmn,jlmn)= if3*(gy_im+gx_re)   ! imag part
2205 
2206          !(Sigma^dn-dn X gij^dn)_x =  gy^dn*f_3
2207          soc_ij(1,2,1,ilmn,jlmn)= if3*gy_re           ! real part
2208          soc_ij(2,2,1,ilmn,jlmn)= if3*gy_im           ! imag part
2209          !(Sigma^dn-dn X gij^dn)_y = -gx^dn*f_3
2210          soc_ij(1,2,2,ilmn,jlmn)=-if3*gx_re           ! real part
2211          soc_ij(2,2,2,ilmn,jlmn)=-if3*gx_im           ! imag part
2212          !(Sigma^dn-dn X gij^dn)_z = 0
2213          soc_ij(1,2,3,ilmn,jlmn)= zero                ! real part
2214          soc_ij(2,2,3,ilmn,jlmn)= zero                ! imag part
2215        end if
2216 
2217      end do ! ilmn
2218    end do ! jlmn
2219      
2220 !  Symetrization
2221    if (option_core==0.and.lmn_size>1) then
2222      do jlmn=2,lmn_size
2223        do ilmn=1,jlmn-1
2224          do ii=1,3
2225            do jj=1,2
2226              avg=half*(soc_ij(1,jj,ii,ilmn,jlmn)+soc_ij(1,jj,ii,jlmn,ilmn))
2227              soc_ij(1,jj,ii,ilmn,jlmn)=avg ; soc_ij(1,jj,ii,jlmn,ilmn)=avg
2228              avg=half*(soc_ij(2,jj,ii,ilmn,jlmn)+soc_ij(2,jj,ii,ilmn,jlmn))
2229              soc_ij(2,jj,ii,ilmn,jlmn)=avg ; soc_ij(2,jj,ii,jlmn,ilmn)=avg
2230            end do
2231          end do           
2232        end do
2233      end do
2234    end if
2235 
2236    ABI_FREE(dVdr)
2237    ABI_FREE(intf3)
2238 
2239  end do  ! iatom
2240  
2241 !Reduction in case of parallelism
2242  if (paral_atom) then
2243    call xmpi_sum(phisocphj,my_comm_atom,ierr)
2244  end if
2245  
2246 !Destroy atom table(s) used for parallelism
2247  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2248 
2249  end subroutine pawnabla_soc_init