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-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