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