TABLE OF CONTENTS
- ABINIT/m_gsphere
- m_gpshere/setshells
- m_gsphere/get_irredg
- m_gsphere/getfullg
- m_gsphere/getkpgnorm
- m_gsphere/gsph_extend
- m_gsphere/gsph_fft_tabs
- m_gsphere/gsph_free
- m_gsphere/gsph_g_idx
- m_gsphere/gsph_gmg_fftidx
- m_gsphere/gsph_gmg_idx
- m_gsphere/gsph_in_fftbox
- m_gsphere/gsph_init
- m_gsphere/gsph_print
- m_gsphere/gsphere_t
- m_gsphere/kg_map
- m_gsphere/make_istwk_table
- m_gsphere/merge_and_sort_kg
- m_gsphere/merge_kgirr
- m_gsphere/setup_G_rotation
- m_gsphere/symg
- m_gsphere/table_gbig2kg
ABINIT/m_gsphere [ Modules ]
NAME
m_gsphere
FUNCTION
The Gsphere data type defines the set of G-vectors centered on Gamma used to describe (chi0|epsilon|W) in the GW code. Note that, unlike the kg_k arrays used for wavefunctions, here the G-vectors are ordered in shells (increasing length). Moreover the sphere can be enlarged to take into account umklapps for which one need the knowledge of several quantities at G-G0.
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (MG, GMR, VO, LR, RWG, MT, XG) 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
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 module m_gsphere 28 29 use defs_basis 30 use m_abicore 31 use m_errors 32 use m_distribfft 33 use m_sort 34 35 use defs_abitypes, only : MPI_type 36 use m_fstrings, only : sjoin, itoa 37 use m_numeric_tools, only : bisect 38 use m_geometry, only : normv 39 use m_crystal, only : crystal_t 40 use m_fftcore, only : kpgsph, kgindex, sphereboundary 41 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 42 43 implicit none 44 45 private 46 47 ! Low-level procedures. 48 public :: merge_and_sort_kg ! Merges a set of k-centered G-spheres of cutoff ecut. Return a Gamma-centered G-spheres. 49 public :: table_gbig2kg ! Associate the kg_k set of G-vectors with Gamma-centered G-sphere. 50 public :: get_irredg ! Given a set of G vectors, find the set of G"s generating the others by symmetry. 51 public :: merge_kgirr ! Merge a list of irreducible G vectors (see routine for more info) 52 public :: setshells ! Set consistently the number of shells, the number of plane-waves, and the energy cut-off 53 public :: kg_map ! Compute the mapping between two lists of g-vectors. 54 public :: make_istwfk_table 55 public :: getkpgnorm ! Compute the norms of the k+G vectors 56 public :: symg
m_gpshere/setshells [ Functions ]
NAME
setshells
FUNCTION
Set consistently the number of shells, the number of plane-waves, and the energy cut-off
INPUTS
nsym=number of symmetry operations gmet(3,3)=metric tensor in reciprocal space gprimd(3,3)=dimensional primitive vectors in reciprocal space symrel(3,3,nsym)=symmetry operations in real space tag=suffix to account for the different possibilities for these variables (npw, ecut or nsh ..) ucvol=unit cell volume
OUTPUT
(see side effects)
SIDE EFFECTS
ecut,npw,nsh=one of them is an input, the two others are output ecut=cut-off energy for plane wave basis sphere (Ha) npw=number of plane waves nsh=number of shells
SOURCE
1542 subroutine setshells(ecut,npw,nsh,nsym,gmet,gprimd,symrel,tag,ucvol) 1543 1544 !Arguments ------------------------------------ 1545 !scalars 1546 integer,intent(in) :: nsym 1547 integer,intent(inout) :: npw,nsh 1548 real(dp),intent(in) :: ucvol 1549 real(dp),intent(inout) :: ecut 1550 character(len=*),intent(in) :: tag 1551 !arrays 1552 integer,intent(in) :: symrel(3,3,nsym) 1553 real(dp),intent(in) :: gmet(3,3),gprimd(3,3) 1554 1555 !Local variables------------------------------- 1556 !scalars 1557 integer :: exchn2n3d,ifound,ig,ii,ish,isym,npw_found,npwave 1558 integer :: npwwrk,nsh_found,pad=50 1559 real(dp) :: ecut_found,ecut_trial,eps,scale=1.3_dp 1560 logical :: found 1561 character(len=500) :: msg 1562 type(MPI_type) :: MPI_enreg_seq 1563 !arrays 1564 integer :: geq(3) 1565 integer,allocatable :: gvec(:,:),gvec_sh(:,:),insort(:),npw_sh(:) 1566 real(dp) :: gctr(3) 1567 real(dp),allocatable :: gnorm(:),gnorm_sh(:) 1568 1569 !****************************************************************** 1570 1571 ! Check coherence of input variables ecut, npw, and nsh. 1572 ! 1-> one at least should be non-null 1573 if (npw==0.and.nsh==0.and.ecut<=tol6) then 1574 write(msg,'(8a)')& 1575 'One of the three variables ecut',TRIM(tag),', npw',TRIM(tag),', or nsh',TRIM(tag),ch10,& 1576 'must be non-null. Returning.' 1577 ABI_COMMENT(msg) 1578 RETURN 1579 end if 1580 ! 2-> one and only one should be non-null 1581 if (npw/=0.and.nsh/=0) then 1582 write(msg,'(6a)')& 1583 'Only one of the two variables npw',TRIM(tag),' and nsh',TRIM(tag),ch10,& 1584 'can be non-null. Modify the value of one of these in input file.' 1585 ABI_ERROR(msg) 1586 end if 1587 if (ecut>tol6.and.npw/=0) then 1588 write(msg,'(6a)')& 1589 'Only one of the two variables ecut',TRIM(tag),' and npw',TRIM(tag),ch10,& 1590 'can be non-null. Modify the value of one of these in input file.' 1591 ABI_ERROR(msg) 1592 end if 1593 if (ecut>tol6.and.nsh/=0) then 1594 write(msg,'(6a)')& 1595 'Only one of the two variables ecut',TRIM(tag),' and nsh',TRIM(tag),ch10,& 1596 'can be non-null Action : modify the value of one of these in input file.' 1597 ABI_ERROR(msg) 1598 end if 1599 1600 ! Calculate an upper bound for npw. 1601 ! gctr is center of the g-vector sphere 1602 gctr(:)= [zero,zero,zero] 1603 if (ecut>tol6) then 1604 ! The average number of plane-waves in the cutoff sphere is given by: 1605 ! npwave = (2*ecut)**(3/2)*ucvol/(6*pi**2) 1606 ! The upper bound is calculated as npwwrk=int(scale * npwave) + pad 1607 npwave=NINT(ucvol*(two*ecut)**1.5_dp/(six*pi**2)) 1608 npwwrk=NINT(DBLE(npwave)*scale)+pad 1609 ecut_trial=ecut 1610 else if (npw/=0) then 1611 ! npw is given in the input 1612 npwwrk=NINT(DBLE(npw)*scale)+pad 1613 ecut_trial=(six*pi**2*npw/ucvol)**two_thirds/two 1614 else 1615 ! If nsh is given in the input 1616 npwwrk=nsh*18+2*pad 1617 ecut_trial=(six*pi**2*nsh*18/ucvol)**two_thirds/two 1618 end if 1619 1620 call initmpi_seq(MPI_enreg_seq) 1621 1622 ABI_MALLOC(gvec,(3,npwwrk)) 1623 ifound=0 1624 do while (ifound==0) 1625 !write(msg,'(a,f8.2)')' setshells : ecut_trial = ',ecut_trial 1626 !call wrtout(std_out,msg,'COLL') 1627 exchn2n3d=0 ! For the time being, no exchange of n2 and n3 1628 1629 call kpgsph(ecut_trial,exchn2n3d,gmet,0,1,1,gvec,gctr,1,MPI_enreg_seq,npwwrk,npw_found) 1630 1631 ABI_MALLOC(gnorm,(npw_found)) 1632 ABI_MALLOC(insort,(npw_found)) 1633 1634 do ig=1,npw_found 1635 insort(ig)=ig 1636 gnorm(ig)=zero 1637 do ii=1,3 1638 gnorm(ig)=gnorm(ig)+(gvec(1,ig)*gprimd(ii,1)+& 1639 gvec(2,ig)*gprimd(ii,2)+& 1640 gvec(3,ig)*gprimd(ii,3))**2 1641 end do 1642 end do 1643 call sort_dp(npw_found,gnorm,insort,tol14) 1644 1645 ABI_MALLOC(npw_sh,(npw_found)) 1646 ABI_MALLOC(gnorm_sh,(npw_found)) 1647 ABI_MALLOC(gvec_sh,(3,npw_found)) 1648 npw_sh(:)=0 1649 gnorm_sh(:)=zero 1650 gvec_sh(:,:)=0 1651 ! Count the number of shells: 1652 ! (search for the G-vectors generating the others by symmetry) 1653 nsh_found=0 1654 1655 do ig=1,npw_found 1656 eps=1.d-8*gnorm(ig) 1657 found=.FALSE. 1658 ish=1 1659 do while ((.not.found).and.(ish<=nsh_found)) 1660 if (ABS(gnorm(ig)-gnorm_sh(ish))<=eps) then 1661 isym=1 1662 do while ((.not.found).and.(isym<=nsym)) 1663 geq(:)=(symrel(1,:,isym)*gvec(1,insort(ig))+& 1664 symrel(2,:,isym)*gvec(2,insort(ig))+& 1665 symrel(3,:,isym)*gvec(3,insort(ig))) 1666 1667 found=((geq(1)==gvec_sh(1,ish)).and.& 1668 (geq(2)==gvec_sh(2,ish)).and.& 1669 (geq(3)==gvec_sh(3,ish))) 1670 isym=isym+1 1671 end do 1672 end if 1673 ish=ish+1 1674 end do 1675 if (.not.found) then 1676 nsh_found=nsh_found+1 1677 gnorm_sh(nsh_found)=gnorm(ig) 1678 gvec_sh(:,nsh_found)=gvec(:,insort(ig)) 1679 npw_sh(nsh_found)=1 1680 else 1681 ish=ish-1 1682 npw_sh(ish)=npw_sh(ish)+1 1683 end if 1684 end do 1685 1686 ecut_found=two*pi**2*gnorm(npw_found) 1687 1688 if(ecut>tol6) then 1689 ! ecut is given in the input 1690 !if (ecut_found<ecut-0.1) then 1691 ! write(msg,'(3a,e14.6,9a,e14.6,3a)')& 1692 ! 'The value ecut',TRIM(tag),'=',ecut,' given in the input file leads to',ch10,& 1693 ! 'the same values for nsh',TRIM(tag),' and npw',TRIM(tag),' as ecut',TRIM(tag),'=',ecut_found,ch10 1694 ! ABI_COMMENT(msg) 1695 !end if 1696 ifound=1 1697 else if (npw/=0) then 1698 ! If npw is given in the input 1699 if (npw_found==npw) then 1700 ecut_found=two*pi**2*gnorm(npw_found) 1701 ifound=1 1702 else if (npw_found>npw) then 1703 npw_found=0 1704 nsh_found=0 1705 do while (npw_found<npw) 1706 nsh_found=nsh_found+1 1707 npw_found=npw_found+npw_sh(nsh_found) 1708 end do 1709 ! check that the shell is closed 1710 if(npw_found>npw) then 1711 ! shell not closed 1712 npw_found=npw_found-npw_sh(nsh_found) 1713 nsh_found=nsh_found-1 1714 do while (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001) 1715 npw_found=npw_found-npw_sh(nsh_found) 1716 nsh_found=nsh_found-1 1717 end do 1718 write(msg,'(3a,i6,5a,i6,3a)')& 1719 'The value npw',TRIM(tag),'=',npw,' given in the input file does not close the shell',ch10,& 1720 'The lower closed-shell is obtained for a value npw',TRIM(tag),'=',npw_found,ch10,& 1721 'This value will be adopted for the calculation.',ch10 1722 ABI_WARNING(msg) 1723 end if 1724 ecut_found=two*pi**2*gnorm(npw_found) 1725 ifound=1 1726 end if 1727 else if (nsh/=0) then 1728 ! If nsh is given in the input 1729 if (nsh_found==nsh) then 1730 ecut_found=two*pi**2*gnorm(npw_found) 1731 ifound=1 1732 else if (nsh_found>nsh) then 1733 npw_found=0 1734 nsh_found=0 1735 do ish=1,nsh 1736 npw_found=npw_found+npw_sh(ish) 1737 nsh_found=nsh_found+1 1738 end do 1739 if (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001) then 1740 do while (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001) 1741 nsh_found=nsh_found+1 1742 npw_found=npw_found+npw_sh(nsh_found) 1743 end do 1744 write(msg,'(3a,i6,5a,i6,3a)')& 1745 'The value nsh',TRIM(tag),'=',nsh,' given in the input file corresponds to the same',ch10,& 1746 'cut-off energy as for closed-shell upto nsh',TRIM(tag),'=',nsh_found,ch10,& 1747 'This value will be adopted for the calculation.',ch10 1748 ABI_WARNING(msg) 1749 end if 1750 ecut_found=two*pi**2*gnorm(npw_found) 1751 ifound=1 1752 end if 1753 end if 1754 1755 if (ifound==0) then 1756 ecut_trial=1.1*ecut_trial 1757 ABI_FREE(gnorm) 1758 ABI_FREE(gnorm_sh) 1759 ABI_FREE(gvec_sh) 1760 ABI_FREE(insort) 1761 ABI_FREE(npw_sh) 1762 else 1763 ! ecut was not provided as an input, then set it now! 1764 if (ecut<tol6) then 1765 ecut=ecut_found 1766 end if 1767 npw=npw_found 1768 nsh=nsh_found 1769 end if 1770 1771 end do ! while(ifound==0) 1772 1773 call destroy_mpi_enreg(MPI_enreg_seq) 1774 1775 ABI_FREE(gnorm) 1776 ABI_FREE(gnorm_sh) 1777 ABI_FREE(gvec) 1778 ABI_FREE(gvec_sh) 1779 ABI_FREE(insort) 1780 ABI_FREE(npw_sh) 1781 1782 end subroutine setshells
m_gsphere/get_irredg [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
get_irredg
FUNCTION
Given a set of reciprocal lattice vectors, find the set of G"s generating the others by symmetry.
INPUTS
nsym=number of symmetry operations pinv=-1 if time-reversal can be used, 1 otherwise npw_k=number of G vectors (for this k-point, as the set of G is k-centered) gcurr(3,npw_k)=the list of G vectors gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) symrec(3,3,nsym)=symmetry operations in terms of reciprocal space primitive translations.
OUTPUT
nbasek=number of irreducible G vectors found cnormk(npw_k)=first nbasek elements are the norm of each irreducible G-vector gbasek(3,npw_k)=first nbasek elements are the irreducible G vectors
NOTES
The search can be optimized by looping over shells. See m_skw for a faster algo
SOURCE
1353 subroutine get_irredg(npw_k,nsym,pinv,gprimd,symrec,gcurr,nbasek,gbasek,cnormk) 1354 1355 !Arguments ------------------------------------ 1356 !scalars 1357 integer,intent(in) :: npw_k,nsym,pinv 1358 integer,intent(out) :: nbasek 1359 !arrays 1360 integer,intent(in) :: gcurr(3,npw_k),symrec(3,3,nsym) 1361 integer,intent(out) :: gbasek(3,npw_k) 1362 real(dp),intent(in) :: gprimd(3,3) 1363 real(dp),intent(out) :: cnormk(npw_k) 1364 1365 !Local variables------------------------------- 1366 !scalars 1367 integer :: ig,irr,isym,jj 1368 real(dp) :: eps,norm 1369 logical :: found 1370 !arrays 1371 integer :: gbas(3),gcur(3),geq(3) 1372 real(dp) :: gcar(3) 1373 1374 ! ************************************************************************* 1375 1376 DBG_ENTER("COLL") 1377 1378 if (pinv/=1.and.pinv/=-1) then 1379 ABI_BUG(sjoin('pinv should be -1 or 1, however, pinv =', itoa(pinv))) 1380 end if 1381 1382 ! Zero irred G vectors found, zeroing output arrays. 1383 nbasek = 0; cnormk(:) = zero; gbasek(:,:) = 0 1384 1385 do ig=1,npw_k 1386 gcur(:) = gcurr(:,ig); norm = zero 1387 do jj=1,3 1388 gcar(jj)=gcur(1)*gprimd(jj,1)+gcur(2)*gprimd(jj,2)+gcur(3)*gprimd(jj,3) 1389 norm=norm+gcar(jj)**2 1390 end do 1391 eps = tol8 * norm; found = .False.; irr = 1 1392 do while (.not.found .and. irr <= nbasek) ! This loop can be optimized by looping inside the shell. 1393 if (abs(norm - cnormk(irr)) <= eps) then 1394 gbas(:) = gbasek(:,irr); isym = 1 1395 do while (.not.found .and. isym <= nsym) 1396 geq(:) = matmul(symrec(:,:,isym),gcur) 1397 found = all(geq(:) == gbas(:)) 1398 if (pinv == -1) found = (found .or. all(geq == -gbas)) ! For time-reversal 1399 isym = isym + 1 1400 end do 1401 end if 1402 irr = irr + 1 1403 end do 1404 if (.not. found) then 1405 nbasek = nbasek + 1; cnormk(nbasek) = norm; gbasek(:,nbasek) = gcur(:) 1406 end if 1407 end do 1408 1409 DBG_EXIT("COLL") 1410 1411 end subroutine get_irredg
m_gsphere/getfullg [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
getfullg
FUNCTION
Reconstruct a G-sphere starting from a set of irreducible lattice vectors
INPUTS
pinv=-1 if time-reversal can be used, 1 otherwise nsym=number of symmetry operations sizepw=Max expected number of G vectors in the shere symrec(3,3,nsym)=symmetry operation in reciprocal space nbase=number of irreducible G vectors gbase(3,nbase)=irreducible G-vectors cnorm(nbase)=norm of the irreducible G vectors (supposed not yet sorted)
OUTPUT
maxpw=Number of G vectors found gbig(3,sizepw)=G vectors in the sphere packed in the first maxpw columns shlim(nbase)=number of G vectors within each shell ierr= Exit status, if /=0 the number of G vectors found exceeds sizepw
SIDE EFFECTS
NOTES
cnorm is a bit redundant since it can be calculated from gbase. However this procedure is called by outkss in which cnorm is already calculated and we dont want to do it twice
SOURCE
1215 subroutine getfullg(nbase,nsym,pinv,sizepw,gbase,symrec,cnorm,maxpw,gbig,shlim,ierr) 1216 1217 !Arguments ------------------------------------ 1218 !scalars 1219 integer,intent(in) :: nbase,nsym,pinv,sizepw 1220 integer,intent(out) :: ierr,maxpw 1221 !arrays 1222 integer,intent(in) :: gbase(3,nbase),symrec(3,3,nsym) 1223 integer,intent(out) :: gbig(3,sizepw),shlim(nbase) 1224 real(dp),intent(inout) :: cnorm(nbase) !sort_dp can change cnorm 1225 1226 !Local variables------------------------------- 1227 !scalars 1228 integer :: ibase,ig,ilim,ish,isym,itim 1229 logical :: found 1230 character(len=500) :: msg 1231 !arrays 1232 integer :: gcur(3),geq(3) 1233 integer,allocatable :: gshell(:,:),insort(:),nshell(:) 1234 1235 ! ************************************************************************* 1236 1237 if (pinv/=1.and.pinv/=-1) then 1238 write(msg,'(a,i6)')& 1239 & ' The argument pinv should be -1 or 1, however, pinv =',pinv 1240 ABI_BUG(msg) 1241 end if 1242 ! 1243 ! === Reorder base g-vectors in order of increasing module === 1244 ABI_MALLOC(insort,(nbase)) 1245 do ibase=1,nbase 1246 insort(ibase)=ibase 1247 end do 1248 call sort_dp(nbase,cnorm,insort,tol14) 1249 ! 1250 ! === Generate all stars of G-vectors === 1251 ! Star of G is the set of all symetrical images of the vector 1252 ! gshell contains the symmetrical G at fixed gbase. No need to add an additional dimension 1253 ! or initialize to zero the array inside the loop over nbase as we loop over (ish<=nshell(ibase)) 1254 ABI_MALLOC(nshell,(nbase)) 1255 ABI_MALLOC(gshell,(3,2*nsym)) 1256 ! 1257 ! === Start with zero number of G vectors found === 1258 maxpw=0 ; ierr=0 1259 do ibase=1,nbase 1260 ! 1261 ! === Loop over all different modules of G === 1262 ! * Start with zero G vectors found in this star 1263 nshell(ibase)=0 1264 gcur(:)=gbase(:,insort(ibase)) 1265 ! 1266 ! === Loop over symmetries === 1267 do isym=1,nsym 1268 do itim=pinv,1,2 1269 geq(:)=itim*MATMUL(symrec(:,:,isym),gcur) 1270 ! 1271 ! * Search for symetric of g and eventually add it: 1272 found=.FALSE. ; ish=1 1273 do while ((.not.found).and. (ish<=nshell(ibase))) 1274 found=ALL(geq(:)==gshell(:,ish)) 1275 ish=ish+1 1276 end do 1277 if (.not.found) then 1278 nshell(ibase)=nshell(ibase)+1 1279 gshell(:,nshell(ibase))=geq(:) 1280 end if 1281 end do 1282 end do 1283 ! 1284 ! * Was sizepw large enough? 1285 if ((maxpw+nshell(ibase))>sizepw) then 1286 write(msg,'(a,i6,2a)')& 1287 & ' Number of G in sphere exceeds maximum allowed value =',sizepw,ch10,& 1288 & ' check the value of sizepw in calling routine ' 1289 ABI_WARNING(msg) 1290 ierr=1; RETURN 1291 end if 1292 ! 1293 ! === Store this shell of Gs in a big array (gbig) === 1294 do ig=1,nshell(ibase) 1295 gbig(:,ig+maxpw)=gshell(:,ig) 1296 end do 1297 maxpw=maxpw+nshell(ibase) 1298 end do ! ibase 1299 ! 1300 ! === Compute number of G"s within each shell === 1301 ilim=0 1302 do ibase=1,nbase 1303 ilim=ilim+nshell(ibase) 1304 shlim(ibase)=ilim 1305 end do 1306 ! 1307 ! === Print out shell limits === 1308 write(msg,'(3a)')& 1309 & ' Shells found:',ch10,& 1310 & ' number of shell number of G vectors cut-off energy [Ha] ' 1311 call wrtout(std_out,msg,'COLL') 1312 1313 do ibase=1,nbase 1314 write(msg,'(12x,i4,17x,i6,12x,f8.3)')ibase,shlim(ibase),two*pi**2*cnorm(ibase) 1315 call wrtout(std_out,msg,'COLL') 1316 end do 1317 write(msg,'(a)')ch10 1318 call wrtout(std_out,msg,'COLL') 1319 ABI_FREE(gshell) 1320 ABI_FREE(insort) 1321 ABI_FREE(nshell) 1322 1323 end subroutine getfullg
m_gsphere/getkpgnorm [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
getkpgnorm
FUNCTION
compute the norms of the k+G vectors
INPUTS
gprimd(3,3)=metric tensor kg_k(3,npw_k)= G vectors, in reduced coordinates kpt(3)=k vector, in reduced coordinates npw_k=size of the G-vector set
OUTPUT
kpgnorm(npw_k)=norms of the k+G vectors
SOURCE
2119 subroutine getkpgnorm(gprimd,kpt,kg_k,kpgnorm,npw_k) 2120 2121 !Arguments ------------------------------------ 2122 !scalars 2123 integer,intent(in) :: npw_k 2124 !arrays 2125 integer,intent(in) :: kg_k(3,npw_k) 2126 real(dp),intent(in) :: gprimd(3,3),kpt(3) 2127 real(dp),intent(out) :: kpgnorm(npw_k) 2128 2129 !Local variables------------------------------- 2130 !scalars 2131 integer :: ipw 2132 real(dp) :: g11,g12,g13,g21,g22,g23,g31,g32,g33,k1,k2,k3,kpg1,kpg2,kpg3,rr,xx 2133 real(dp) :: yy,zz 2134 2135 ! ************************************************************************* 2136 2137 k1=kpt(1) ; k2=kpt(2) ; k3=kpt(3) 2138 g11=gprimd(1,1) 2139 g12=gprimd(1,2) 2140 g13=gprimd(1,3) 2141 g21=gprimd(2,1) 2142 g22=gprimd(2,2) 2143 g23=gprimd(2,3) 2144 g31=gprimd(3,1) 2145 g32=gprimd(3,2) 2146 g33=gprimd(3,3) 2147 2148 !Loop over all k+G 2149 do ipw=1,npw_k 2150 2151 ! Load k+G 2152 kpg1=k1+dble(kg_k(1,ipw)) 2153 kpg2=k2+dble(kg_k(2,ipw)) 2154 kpg3=k3+dble(kg_k(3,ipw)) 2155 2156 ! Calculate module of k+G 2157 xx=g11*kpg1+g12*kpg2+g13*kpg3 2158 yy=g21*kpg1+g22*kpg2+g23*kpg3 2159 zz=g31*kpg1+g32*kpg2+g33*kpg3 2160 rr=sqrt(xx**2+yy**2+zz**2) 2161 kpgnorm(ipw) = rr 2162 2163 end do ! ipw 2164 2165 end subroutine getkpgnorm
m_gsphere/gsph_extend [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_extend
FUNCTION
Construct a new gsphere_t with a larger cutoff energy while preserving the ordering of the first G-vectors stored in in_Gsph
INPUTS
OUTPUT
SOURCE
2032 subroutine gsph_extend(in_Gsph, Cryst, new_ecut, new_Gsph) 2033 2034 !Arguments ------------------------------------ 2035 !scalars 2036 class(gsphere_t),intent(in) :: in_Gsph 2037 type(crystal_t),intent(in) :: Cryst 2038 real(dp),intent(in) :: new_ecut 2039 class(gsphere_t),intent(out) :: new_Gsph 2040 2041 !Local variables------------------------------- 2042 !scalars 2043 integer :: new_ng,in_ng,ig,ierr,sh 2044 !arrays 2045 integer,allocatable :: new_gvec(:,:) 2046 2047 ! ********************************************************************* 2048 2049 call new_Gsph%init(Cryst, 0, ecut=new_ecut) 2050 2051 if (new_Gsph%ng > in_Gsph%ng) then 2052 2053 new_ng = new_Gsph%ng 2054 in_ng = in_Gsph%ng 2055 2056 ierr = 0 2057 do ig=1,in_ng 2058 if (ANY(new_Gsph%gvec(:,ig) /= in_Gsph%gvec(:,ig)) ) then 2059 ierr = ierr + 1 2060 write(std_out,*)" new_gvec, in_gvec",ig,new_Gsph%gvec(:,ig),in_Gsph%gvec(:,ig) 2061 end if 2062 end do 2063 2064 if (ierr==0) RETURN 2065 2066 ierr = 0 2067 do sh=1,in_Gsph%nsh 2068 if ( new_Gsph%shlim(sh) /= in_Gsph%shlim(sh) .or. & 2069 & ABS(new_Gsph%shlen(sh)-in_Gsph%shlen(sh)) > tol12 ) then 2070 ierr = ierr + 1 2071 write(std_out,*)"new_shlim, in_shlim",sh,new_Gsph%shlim(sh),in_Gsph%shlim(sh) 2072 write(std_out,*)"new_shlen, in_shlen",sh,new_Gsph%shlen(sh),in_Gsph%shlen(sh) 2073 end if 2074 end do 2075 ABI_CHECK(ierr==0,"Wrong shells") 2076 2077 ABI_MALLOC(new_gvec,(3,new_ng)) 2078 new_gvec = new_Gsph%gvec 2079 new_gvec(:,1:in_ng) = in_Gsph%gvec 2080 2081 call new_Gsph%free() 2082 call new_Gsph%init(Cryst, new_ng, gvec=new_gvec) 2083 ABI_FREE(new_gvec) 2084 2085 else 2086 ierr = 0 2087 do ig=1,MIN(new_Gsph%ng,in_Gsph%ng) 2088 if (ANY(new_Gsph%gvec(:,ig) /= in_Gsph%gvec(:,ig)) ) then 2089 ierr = ierr + 1 2090 write(std_out,*)" new_gvec, in_gvec",ig,new_Gsph%gvec(:,ig),in_Gsph%gvec(:,ig) 2091 end if 2092 end do 2093 ABI_CHECK(ierr==0,"Fatal error") 2094 end if 2095 2096 end subroutine gsph_extend
m_gsphere/gsph_fft_tabs [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_fft_tabs
FUNCTION
INPUTS
Gsph<gsphere_t>=Info on the G-sphere g0(3) mgfft=MAXVAL(ngfft(1:3)) ngfftf(18)=Info on the FFT mesh.
OUTPUT
use_padfft=1 if padded FFT can be used, 0 otherwise. gmg0_gbound(2*mgfft+8,2)=Tables for improved zero-padded FFTS. Calculated only if use_padfft==1 gmg0_ifft(Gsph%ng)=Index of G-G0 in the FFT mesh defined by ngfft.
NOTES
The routine will stop if any G-G0 happens to be outside the FFT box.
SOURCE
473 subroutine gsph_fft_tabs(Gsph, g0, mgfft, ngfft, use_padfft, gmg0_gbound, gmg0_ifft) 474 475 !Arguments ------------------------------------ 476 !scalars 477 class(gsphere_t),intent(in) :: Gsph 478 integer,intent(in) :: mgfft 479 integer,intent(out) :: use_padfft 480 !arrays 481 integer,intent(in) :: g0(3),ngfft(18) 482 integer,intent(out) :: gmg0_gbound(2*mgfft+8,2),gmg0_ifft(Gsph%ng) 483 484 !Local variables------------------------------- 485 !scalars 486 integer :: ig,ng,ierr 487 character(len=500) :: msg 488 type(MPI_type) :: MPI_enreg_seq 489 !arrays 490 integer,allocatable :: gmg0(:,:) 491 logical,allocatable :: kg_mask(:) 492 493 ! ************************************************************************* 494 495 if (mgfft/=MAXVAL(ngfft(1:3))) then 496 ABI_ERROR("mgfft/-MAXVAL(ngfft(1:3)") 497 end if 498 499 ng = Gsph%ng 500 501 ierr=0; use_padfft=0 502 ABI_MALLOC(gmg0,(3,ng)) 503 do ig=1,ng 504 gmg0(:,ig) = Gsph%gvec(:,ig)-g0 505 ! Consider possible wrap around errors. 506 if ( ANY(gmg0(:,ig)>ngfft(1:3)/2) .or. ANY(gmg0(:,ig)<-(ngfft(1:3)-1)/2) ) then 507 !gmg0_ifft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 0 508 write(std_out,*)" outside FFT box ",gmg0(:,ig) 509 ierr=ierr+1 510 end if 511 if (ALL(gmg0(:,ig) == 0)) use_padfft=1 512 end do 513 514 if (ierr/=0) then 515 write(msg,'(a,i0,a)')'Found ',ierr,' G-G0 vectors falling outside the FFT box. This is not allowed ' 516 ABI_ERROR(msg) 517 end if 518 ! 519 ! Evaluate the tables needed for the padded FFT performed in rhotwg. Note that we have 520 ! to pass G-G0 to sphereboundary instead of G as we need FFT results on the shifted G-sphere, 521 ! If Gamma is not inside G-G0 one has to disable FFT padding as sphereboundary will give wrong tables. 522 if (use_padfft == 1) call sphereboundary(gmg0_gbound,1,gmg0,mgfft,ng) 523 524 call initmpi_seq(MPI_enreg_seq) ! No FFT parallelism. 525 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all') 526 527 ABI_MALLOC(kg_mask, (ng)) 528 call kgindex(gmg0_ifft, gmg0, kg_mask, MPI_enreg_seq, ngfft, ng) 529 530 ABI_CHECK(ALL(kg_mask),"FFT para not yet implemented") 531 ABI_FREE(kg_mask) 532 533 ABI_FREE(gmg0) 534 call destroy_mpi_enreg(MPI_enreg_seq) 535 536 end subroutine gsph_fft_tabs
m_gsphere/gsph_free [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_free
FUNCTION
Deallocate the memory in a gsphere_t data type.
INPUTS
Gsph = datatype to be freed
SOURCE
719 subroutine gsph_free(Gsph) 720 721 !Arguments ------------------------------------ 722 class(gsphere_t),intent(inout) :: Gsph 723 724 ! ************************************************************************* 725 726 DBG_ENTER("COLL") 727 728 !@gsphere_t 729 730 ! integer arrays. 731 ABI_SFREE(Gsph%g2sh) 732 ABI_SFREE(Gsph%gvec) 733 ABI_SFREE(Gsph%g2mg) 734 ABI_SFREE(Gsph%rottb) 735 ABI_SFREE(Gsph%rottbm1) 736 ABI_SFREE(Gsph%shlim) 737 738 ABI_SFREE(Gsph%shlen) 739 740 ! complex arrays 741 ABI_SFREE(Gsph%phmGt) 742 ABI_SFREE(Gsph%phmSGt) 743 744 DBG_EXIT("COLL") 745 746 end subroutine gsph_free
m_gsphere/gsph_g_idx [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_g_idx
FUNCTION
Return the index of G in the sphere. zero if not in the sphere
INPUTS
Gsph<gsphere_t>=Info on the G-sphere gg(3)=Reduced coordinates of the G-vector.
NOTES
The function assumes that the G-vectors are ordered with increasing length.
SOURCE
767 pure function gsph_g_idx(Gsph, gg) result(g_idx) 768 769 !Arguments ------------------------------------ 770 !scalars 771 class(gsphere_t),intent(in) :: Gsph 772 integer :: g_idx 773 !arrays 774 integer,intent(in) :: gg(3) 775 776 !Local variables------------------------------- 777 !scalars 778 integer :: ishbsc,igs,ige 779 real(dp) :: glen 780 logical :: found 781 782 ! ************************************************************************* 783 784 ! * Use shells and bisect to find the star and stop index thus avoiding the storage of a table (ig1,ig2) 785 glen = two_pi*SQRT(DOT_PRODUCT(gg,MATMUL(Gsph%gmet,gg))) 786 787 ishbsc = bisect(Gsph%shlen,glen) 788 if ( ANY(ishbsc==(/0,Gsph%nsh/)) ) then ! glen out of range. 789 g_idx=0; RETURN 790 end if 791 792 igs = Gsph%shlim(ishbsc) 793 ige = Gsph%shlim(MIN(ishbsc+2,Gsph%nsh+1))-1 794 795 g_idx=igs-1; found=.FALSE. 796 do while (.not.found .and. g_idx<ige) 797 g_idx=g_idx+1 798 found=(ALL(Gsph%gvec(:,g_idx)==gg(:))) 799 end do 800 if (.not.found) g_idx=0 801 802 end function gsph_g_idx
m_gsphere/gsph_gmg_fftidx [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_gmg_fftidx
FUNCTION
Return the index of G1-G2 in the FFT mesh defined by ngfft. zero if not found.
INPUTS
Gsph<gsphere_t>=Info on the G-sphere ig1,ig2 index of g1 and g2 in the G-sphere. ngfft(18)=Info on the FFT mesh.
SOURCE
881 pure function gsph_gmg_fftidx(Gsph, ig1, ig2, ngfft) result(fft_idx) 882 883 !Arguments ------------------------------------ 884 !scalars 885 class(gsphere_t),intent(in) :: Gsph 886 integer,intent(in) :: ig1,ig2 887 integer :: fft_idx 888 !arrays 889 integer,intent(in) :: ngfft(18) 890 891 !Local variables------------------------------- 892 !scalars 893 integer :: id1,id2,id3 894 !arrays 895 integer :: g1mg2(3) 896 897 ! ************************************************************************* 898 899 g1mg2(:)=Gsph%gvec(:,ig1)-Gsph%gvec(:,ig2) 900 901 ! Make sure G1-G2 is still in the FFT mesh. 902 ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic! 903 if (ANY(g1mg2(:)>ngfft(1:3)/2) .or. ANY(g1mg2(:)<-(ngfft(1:3)-1)/2)) then 904 fft_idx=0; RETURN 905 end if 906 907 id1=MODULO(g1mg2(1),ngfft(1)) 908 id2=MODULO(g1mg2(2),ngfft(2)) 909 id3=MODULO(g1mg2(3),ngfft(3)) 910 fft_idx= 1 + id1 + id2*ngfft(1) + id3*ngfft(1)*ngfft(2) 911 912 end function gsph_gmg_fftidx
m_gsphere/gsph_gmg_idx [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_gmg_idx
FUNCTION
Return the index of G1-G2 in the sphere. zero if not in the sphere
INPUTS
Gsph<gsphere_t>=Info on the G-sphere ig1,ig2 index of g1 and g2 in the G-sphere.
NOTES
The function assumes that the G-vectors are ordered with increasing length.
SOURCE
823 pure function gsph_gmg_idx(Gsph, ig1, ig2) result(ig1mg2) 824 825 !Arguments ------------------------------------ 826 !scalars 827 class(gsphere_t),intent(in) :: Gsph 828 integer,intent(in) :: ig1,ig2 829 integer :: ig1mg2 830 831 !Local variables------------------------------- 832 !scalars 833 integer :: ishbsc,igs,ige 834 real(dp) :: difflen 835 logical :: found 836 !arrays 837 integer :: g1mg2(3) 838 839 ! ************************************************************************* 840 841 g1mg2 = Gsph%gvec(:,ig1)-Gsph%gvec(:,ig2) 842 843 ! * Use shells and bisect to find the star and stop index thus avoiding the storage of a table (ig1,ig2) 844 difflen = two_pi*SQRT(DOT_PRODUCT(g1mg2,MATMUL(Gsph%gmet,g1mg2))) 845 846 ! FIXME It seems bisect is not portable, on my laptop test v5/t72 the number of skipped G-vectors is > 0 847 ishbsc = bisect(Gsph%shlen,difflen) 848 if ( ANY(ishbsc==(/0,Gsph%nsh/)) ) then ! difflen out of range. 849 ig1mg2=0; RETURN 850 end if 851 852 igs = Gsph%shlim(ishbsc) 853 ige = Gsph%shlim(MIN(ishbsc+2,Gsph%nsh+1))-1 854 855 ig1mg2=igs-1; found=.FALSE. 856 do while (.not.found .and. ig1mg2<ige) 857 ig1mg2=ig1mg2+1 858 found=(ALL(Gsph%gvec(:,ig1mg2)==g1mg2(:))) 859 end do 860 if (.not.found) ig1mg2=0 861 862 end function gsph_gmg_idx
m_gsphere/gsph_in_fftbox [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_in_fftbox
FUNCTION
Initialize the largest Gsphere contained in the FFT box.
INPUTS
Cryst<crystal_t> = Info on unit cell and its symmetries. ngfft(18)=Info on the FFT box.
OUTPUT
Gsph<gsphere_t>=Data type containing information related to the set of G vectors completetly initialized in output.
SOURCE
558 subroutine gsph_in_fftbox(Gsph, Cryst, ngfft) 559 560 !Arguments ------------------------------------ 561 !scalars 562 class(gsphere_t),intent(out) :: Gsph 563 type(crystal_t),intent(in) :: Cryst 564 !arrays 565 integer,intent(in) :: ngfft(18) 566 567 !Local variables------------------------------- 568 !scalars 569 integer :: dir1,dir2,dir3,npw,ig,ist 570 real(dp) :: ecut,trial_ene 571 !arrays 572 integer :: n1_max(3),n2_max(3),n3_max(3),vec(3) 573 integer,allocatable :: gvec(:,:) 574 575 !************************************************************************ 576 577 ! Find ecut for the largest G-sphere contained in the FFT box. 578 n1_max(1) = -(ngfft(1)-1)/2 579 n2_max(1) = -(ngfft(2)-1)/2 580 n3_max(1) = -(ngfft(3)-1)/2 581 582 n1_max(2) = 0 583 n2_max(2) = 0 584 n3_max(2) = 0 585 586 n1_max(3) = ngfft(1)/2 587 n2_max(3) = ngfft(2)/2 588 n3_max(3) = ngfft(3)/2 589 590 ecut = HUGE(one) 591 do dir3=1,3 592 vec(3) = n1_max(dir3) 593 do dir2=1,3 594 vec(2) = n2_max(dir2) 595 do dir1=1,3 596 vec(1) = n1_max(dir1) 597 if (ANY(vec/=0)) then 598 trial_ene = half * normv(vec,Cryst%gmet,"G")**2 599 ecut = MIN(ecut,trial_ene) 600 !write(std_out,*)vec(:),trial_ene 601 end if 602 end do 603 end do 604 end do 605 ! 606 ! Init sphere from ecut. 607 call Gsph%init(Cryst, 0, ecut=ecut) 608 ! 609 ! Make sure that Gsph does not contain G vectors outside the FFT box. 610 ! kpgsph might return G whose energy is larger than the input ecut. 611 npw = Gsph%ng 612 star_loop: do ist=1,Gsph%nsh-1 613 do ig=Gsph%shlim(ist),Gsph%shlim(ist+1) 614 if ( ANY(Gsph%gvec(:,ig)>ngfft(1:3)/2) .or. ANY(Gsph%gvec(:,ig)<-(ngfft(1:3)-1)/2) ) then 615 npw = Gsph%shlim(ist)-1 ! Gsph exceeds the FFT box. Only the shells up to npw will be used. 616 EXIT star_loop 617 end if 618 end do 619 end do star_loop 620 621 if (npw<Gsph%ng) then 622 ABI_COMMENT("Have to reinit Gpshere") 623 ABI_MALLOC(gvec,(3,npw)) 624 gvec = Gsph%gvec(:,1:npw) 625 call Gsph%free() 626 call Gsph%init(Cryst, npw, gvec=gvec) 627 ABI_FREE(gvec) 628 end if 629 630 end subroutine gsph_in_fftbox
m_gsphere/gsph_init [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_init
FUNCTION
Main creation method for the Gvectors data type
INPUTS
Cryst<crystal_t> = Info on unit cell and its symmetries ng=number of G vectors, needed only if gvec is passed. [gvec(3,ng)]=coordinates of G vectors [ecut]=Cutoff energy for G-sphere. gvec and ecut are mutually exclusive.
OUTPUT
Gsph<gsphere_t>=Data type containing information related to the set of G vectors completetly initialized in output.
NOTES
gvec are supposed to be ordered with increasing norm.
SOURCE
281 subroutine gsph_init(Gsph, Cryst, ng, gvec, ecut) 282 283 !Arguments ------------------------------------ 284 !scalars 285 class(gsphere_t),intent(out) :: Gsph 286 integer,intent(in) :: ng 287 real(dp),optional,intent(in) :: ecut 288 type(crystal_t),target,intent(in) :: Cryst 289 290 !arrays 291 integer,optional,intent(in) :: gvec(3,ng) 292 !Local variables------------------------------- 293 !scalars 294 integer,parameter :: nkpt1=1 295 integer :: ig,isearch,img,ish,isym,nsh,nsym,timrev,pinv,g1,g2,g3,ss,ee 296 real(dp) :: eps,norm,norm_old,max_ecut,gsq 297 !arrays 298 real(dp),parameter :: k_gamma(3)=(/zero,zero,zero/) 299 integer :: sg(3),gsearch(3) 300 integer,allocatable :: shlim(:) 301 integer,pointer :: symrec(:,:,:),gvec_ptr(:,:) 302 real(dp) :: kptns1(3,nkpt1) 303 real(dp),allocatable :: shlen(:) 304 real(dp),pointer :: tnons(:,:) 305 306 !************************************************************************ 307 308 DBG_ENTER("COLL") 309 310 !@gsphere_t 311 ! === Copy info on symmetries === 312 nsym = Cryst%nsym 313 timrev = Cryst%timrev 314 symrec => Cryst%symrec 315 tnons => Cryst%tnons 316 ! 317 ! === Initialize the object === 318 Gsph%istwfk = 1 ! Time reversal is not used here. 319 Gsph%nsym = nsym 320 Gsph%timrev = timrev 321 322 Gsph%gmet = Cryst%gmet 323 Gsph%gprimd = Cryst%gprimd 324 325 if (PRESENT(gvec)) then 326 if (PRESENT(ecut)) then 327 ABI_BUG("ecut cannot be present when gvec is used") 328 end if 329 Gsph%ng= ng 330 ABI_MALLOC(Gsph%gvec,(3,ng)) 331 Gsph%gvec=gvec 332 ! 333 ! Calculate cutoff energy.of the sphere. 334 max_ecut=-one 335 do ig=1,ng 336 g1=gvec(1,ig) 337 g2=gvec(2,ig) 338 g3=gvec(3,ig) 339 gsq= Cryst%gmet(1,1)*g1**2+Cryst%gmet(2,2)*g2**2+Cryst%gmet(3,3)*g3**2+ & 340 & two*(Cryst%gmet(1,2)*g1*g2+Cryst%gmet(1,3)*g1*g3+Cryst%gmet(2,3)*g2*g3) 341 max_ecut=MAX(max_ecut,gsq) 342 end do 343 max_ecut=two*max_ecut*pi**2 344 Gsph%ecut= max_ecut 345 346 else 347 ! To be consistent with the previous implementation. 348 ABI_WARNING("Init from ecut has to be tested") 349 !call setshells(ecut,npw,nsh,nsym,Cryst%gmet,Cryst%gprimd,Cryst%symrel,tag,Cryst%ucvol) 350 Gsph%ecut = ecut 351 pinv=+1; kptns1(:,1)=k_gamma 352 call merge_and_sort_kg(nkpt1,kptns1,ecut,Cryst%nsym,pinv,Cryst%symrel,Cryst%gprimd,gvec_ptr,0) 353 Gsph%ng = SIZE(gvec_ptr,DIM=2) 354 ABI_MALLOC(Gsph%gvec, (3,Gsph%ng)) 355 Gsph%gvec = gvec_ptr 356 ABI_FREE(gvec_ptr) 357 end if 358 ! 359 ! Calculate phase exp{-i2\pi G.\tau} 360 ABI_MALLOC(Gsph%phmGt, (Gsph%ng,nsym)) 361 do isym=1,nsym 362 do ig=1,Gsph%ng 363 Gsph%phmGt(ig, isym) = EXP(-j_dpc*two_pi*DOT_PRODUCT(Gsph%gvec(:,ig), tnons(:,isym))) 364 end do 365 end do 366 ! 367 ! === Calculate phase phsgt= exp{-i2\pi SG\cdot t} === 368 ! TODO Here we can store only one of this arrays but I have to rewrite screeening! 369 ABI_MALLOC(Gsph%phmSGt,(Gsph%ng,nsym)) 370 do ig=1,Gsph%ng 371 do isym=1,nsym 372 sg=MATMUL(symrec(:,:,isym),Gsph%gvec(:,ig)) 373 Gsph%phmSGt(ig,isym)=EXP(-j_dpc*two_pi*DOT_PRODUCT(sg,tnons(:,isym))) 374 end do 375 end do 376 ! 377 ! === Calculate number of shells and corresponding starting index === 378 ! * Shells are useful to speed up search algorithms see e.g setup_G_rotation. 379 ! * The last shell ends at ng+1, thus gvec is supposed to be closed. 380 381 ABI_CHECK(ALL(Gsph%gvec(1:3,1)==0),'First G must be 0') 382 383 ABI_MALLOC(Gsph%g2sh,(Gsph%ng)) 384 Gsph%g2sh(1)=1 ! This table is useful if we dont loop over shell 385 386 ! For each shell, gives the index of the initial G-vector. 387 ABI_MALLOC(shlim,(Gsph%ng+1)) 388 shlim(1)=1 389 390 ! For each shell, gives the radius of the shell. 391 ABI_MALLOC(shlen,(Gsph%ng)) 392 shlen(1)=zero 393 394 nsh=1; norm_old=zero 395 do ig=2,Gsph%ng 396 norm=two_pi*SQRT(DOT_PRODUCT(Gsph%gvec(:,ig),MATMUL(Cryst%gmet,Gsph%gvec(:,ig)))) 397 eps=norm*tol8 398 if (ABS(norm-norm_old)>eps) then 399 norm_old = norm; nsh = nsh + 1 400 shlim(nsh)=ig 401 shlen(nsh)=norm 402 end if 403 Gsph%g2sh(ig)=nsh 404 end do 405 shlim(nsh+1)=Gsph%ng+1 406 407 ! === Save info on the shells === 408 Gsph%nsh=nsh 409 ABI_MALLOC(Gsph%shlim,(nsh+1)) 410 Gsph%shlim=shlim(1:nsh+1) 411 ABI_MALLOC(Gsph%shlen,(nsh )) 412 Gsph%shlen=shlen(1:nsh) 413 ABI_FREE(shlim) 414 ABI_FREE(shlen) 415 416 ! Calculate tables for rotated G"s 417 ABI_MALLOC(Gsph%rottb , (Gsph%ng,timrev,nsym)) 418 ABI_MALLOC(Gsph%rottbm1, (Gsph%ng,timrev,nsym)) 419 420 call setup_G_rotation(nsym, symrec, timrev, Gsph%ng, Gsph%gvec,& 421 Gsph%g2sh, Gsph%nsh, Gsph%shlim, Gsph%rottb, Gsph%rottbm1) 422 423 ! Store Mapping G --> -G 424 ! (we use a specialized table instead of rootb since rottb assumes time-reversal symmetry. 425 ABI_MALLOC(gsph%g2mg, (gsph%ng)) 426 427 do ig=1,gsph%ng 428 ish=gsph%g2sh(ig) 429 ss=gsph%shlim(ish); ee=gsph%shlim(ish+1)-1 430 gsearch = -gsph%gvec(:,ig) 431 img = 0 432 ! Loop on shells to speed up the search. 433 do isearch=ss,ee 434 if (all(gsph%gvec(:,isearch) == gsearch)) then 435 img = isearch; exit 436 end if 437 end do 438 if (img==0) ABI_ERROR("Cannot find -G in G-sphere!") 439 gsph%g2mg(ig) = img 440 end do 441 442 !call Gsph%print(unit=std_out,prtvol=1) 443 444 DBG_EXIT("COLL") 445 446 end subroutine gsph_init
m_gsphere/gsph_print [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
gsph_print
FUNCTION
Print the content of a gvectors data type
INPUTS
Gsph<gsphere_t>=Info on the G-sphere unit=the unit number for output prtvol = verbosity level mode_paral =either "COLL" or "PERS"
OUTPUT
Only writing.
SOURCE
653 subroutine gsph_print(Gsph, unit, prtvol, mode_paral) 654 655 !Arguments ------------------------------------ 656 !scalars 657 class(gsphere_t),intent(in) :: Gsph 658 integer,intent(in),optional :: prtvol,unit 659 character(len=4),intent(in),optional :: mode_paral 660 661 !Local variables------------------------------- 662 !scalars 663 integer :: ish,nsc,my_unt,my_prtvol 664 real(dp) :: fact,kin 665 character(len=4) :: my_mode 666 character(len=500) :: msg 667 668 ! ************************************************************************* 669 670 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 671 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 672 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 673 674 write(msg,'(3a,2(a,i8,a))')ch10,& 675 ' ==== Info on the G-sphere ==== ',ch10,& 676 ' Number of G vectors ... ',Gsph%ng,ch10,& 677 ' Number of shells ...... ',Gsph%nsh,ch10 678 call wrtout(my_unt,msg,my_mode) 679 680 SELECT CASE (Gsph%timrev) 681 CASE (1) 682 call wrtout(my_unt,' Time reversal symmetry cannot be used',my_mode) 683 CASE (2) 684 call wrtout(my_unt,' Time reversal symmetry is used',my_mode) 685 CASE DEFAULT 686 ABI_BUG("Wrong timrev") 687 END SELECT 688 689 if (my_prtvol/=0) then 690 fact=half*two_pi**2 691 write(msg,'(a)') 692 call wrtout(my_unt,' Shell Tot no. of Gs Cutoff [Ha]',my_mode) 693 do ish=1,Gsph%nsh 694 nsc=Gsph%shlim(ish+1)-1 695 kin=half*Gsph%shlen(ish)**2 696 write(msg,'(2x,i4,10x,i6,5x,f8.3)')ish,nsc,kin 697 call wrtout(my_unt,msg,'COLL') 698 end do 699 call wrtout(my_unt,ch10,my_mode) 700 end if 701 702 end subroutine gsph_print
m_gsphere/gsphere_t [ Types ]
[ Top ] [ m_gsphere ] [ Types ]
NAME
gsphere_t
FUNCTION
The gsphere_t data type contains information related to the set of G vectors used during a screening or a GW calculation, as well as symmetry tables relating these vectors. Presently the following quantities are stored 1) The reduced coordinates of the G vectors (arrays gvec) 2) Tables giving the correspondence between a G-vector and its rotated image through a symmetry operation in reciprocal space. 3) List of the irreducible G pairs 4) Tables giving, for each pair in the full reciprocal space, the corresponding irreducible pair as well as the symmetry operation in reciprocal space Note that, unlike the GS part, the basis set does not depend on the k-point.
NOTES
To indicate the indices in the arrays grottb, grottbm1 we use the following notation: g defines the index of the reciprocal lattice vector in the array gvec s indicates the index of the symmetry operation in reciprocal space i can be one or two. 1 is used to indicate the identity operator
SOURCE
89 type,public :: gsphere_t 90 91 integer :: ng 92 ! Total number of G vectors in the sphere taking into account umklapps 93 ! it defines the size of the array gvec and it accounts for possible umklapps for which 94 ! one has to shift the sphere. 95 96 ! TODO: The sphere should be enlarged in gpairs_init using mg0 and (ecut|ng) as input. 97 ! For the time being we keep the old implementation. 98 !%integer :: ng_eff 99 ! Effective number of G vectors, i.e. the number of G in the smaller sphere without umklapps. 100 ! ng_eff<=ng and should be used to loop over the elements of (chi0|epsilon|W). 101 ! 102 ! TODO: Add info on FFT including zero-padding algorithm. 103 ! table must be recalculated for each G0 and rho_tw_g should accept gpshere in input. 104 105 integer :: nsh ! Number of shells 106 integer :: nsym ! The number of symmetry operations 107 integer :: timrev ! 2 if time-reversal is used, 1 otherwise 108 integer :: istwfk=1 ! Storage mode. At present time-reversal is not used. 109 110 !integer :: mg0(3)=0 111 ! For each reduced direction gives the max G0 component to account for umklapp processes. 112 113 real(dp) :: ecut 114 ! Cutoff energy of the sphere. 115 116 real(dp) :: gmet(3,3) 117 ! Reciprocal space metric ($\textrm{bohr}^{-2}$). 118 119 real(dp) :: gprimd(3,3) 120 ! Dimensional reciprocal space primitive translations (bohr^{-1}) 121 122 integer,allocatable :: g2sh(:) 123 ! g2sh(ng) 124 ! For each G, it gives the index of the shell to which it belongs. 125 126 integer,allocatable :: gvec(:,:) 127 ! gvec(3,ng) 128 ! Reduced coordinates of G vectors. 129 130 integer,allocatable :: g2mg(:) 131 ! g2mg(ng) 132 ! Correspondence G --> -G 133 134 integer,allocatable :: rottb(:,:,:) 135 ! rottb(ng,timrev,nsym) 136 ! rottb(G,I,S) is the index of (SI) G in the array gvec 137 ! where I is either the identity or the inversion. 138 139 integer,allocatable :: rottbm1(:,:,:) 140 ! rottb(ng,timrev,nsym) 141 ! rottbm1(G,I,S) is the index of IS{^-1} G in the array gvec 142 143 integer,allocatable :: shlim(:) 144 ! shlim(nsh+1) 145 ! Index of the first G vector in each shell, =ng+1 for nsh+1 146 147 real(dp),allocatable :: shlen(:) 148 ! shlen(nsh) 149 ! Radius of each shell. 150 151 !TODO switch to dpc 152 complex(gwpc),allocatable :: phmGt(:,:) 153 ! phmGt(ng,nsym) 154 ! Phase factor e^{-i2\pi(G.\tau)} where $\tau$ is the fractional translation associated to isym. 155 156 complex(gwpc),allocatable :: phmSGt(:,:) 157 ! phmSGt(ng,nsym) 158 ! Phase factor e^{-i2\pi(SG.\tau)} where S is one of the symmetry properties in reciprocal space. 159 160 contains 161 162 procedure :: init => gsph_init ! Initialize the G-sphere. 163 procedure :: fft_tabs => gsph_fft_tabs ! Returns useful tables for FFT (with or without padding). 164 procedure :: in_fftbox => gsph_in_fftbox ! Initialize the largest Gsphere contained in the FFT box. 165 procedure :: print => gsph_print ! Printout of basic dimensions. 166 procedure :: free => gsph_free ! Free memory allocated in the object. 167 procedure :: g_idx => gsph_g_idx ! Returns the index of G from its reduced coordinates. 168 procedure :: gmg_idx => gsph_gmg_idx ! Returns the index of G1-G2 from their indices 169 procedure :: gmg_fftidx => gsph_gmg_fftidx ! Returns the index of G1-G2 in the FFT mesh defined by ngfft. 170 procedure :: extend => gsph_extend ! Construct a new gsphere_t with a larger cutoff energy 171 172 end type gsphere_t
m_gsphere/kg_map [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
kg_map
FUNCTION
Compute the mapping between two lists of g-vectors.
INPUTS
npw1, kg1(3,npw1)=First list of G-vectors npw2, kg2(3,npw2)=Second list of G-vectors
OUTPUT
g2g1(npw2) = Mapping kg2 index --> kg1 index.
Set to 0 if kg2(:,ig) not in kg1
nmiss = Number of G-vectors in kg2 not found in kg1
SOURCE
1805 subroutine kg_map(npw1, kg1, npw2, kg2, g2g1, nmiss) 1806 1807 !Arguments ------------------------------------ 1808 !scalars 1809 integer,intent(in) :: npw1,npw2 1810 integer,intent(out) :: nmiss 1811 !arrays 1812 integer,intent(in) :: kg1(3,npw1),kg2(3,npw2) 1813 integer,intent(out) :: g2g1(npw2) 1814 1815 !Local variables ------------------------------ 1816 !scalars 1817 integer :: ii,ipw,i1,i2,i3,n1,n2,n3 1818 !arrays 1819 integer :: gmax(3),g1_max(3),g2_max(3) 1820 integer,allocatable :: iwork(:,:,:) 1821 1822 !************************************************************************ 1823 1824 g1_max = maxval(abs(kg1)) 1825 g2_max = maxval(abs(kg2)) 1826 do ii=1,3 1827 gmax(ii) = max(g1_max(ii), g2_max(ii)) 1828 end do 1829 gmax = 2*gmax + 1 1830 n1 = gmax(1); n2 = gmax(2); n3 = gmax(3) 1831 1832 ABI_MALLOC(iwork, (n1, n2, n3)) 1833 1834 ! Insert kg1 into work with extra 0 s around outside: 1835 iwork = 0 1836 do ipw=1,npw1 1837 i1 = kg1(1,ipw); if (i1<0) i1=i1+n1; i1=i1+1 1838 i2 = kg1(2,ipw); if (i2<0) i2=i2+n2; i2=i2+1 1839 i3 = kg1(3,ipw); if (i3<0) i3=i3+n3; i3=i3+1 1840 iwork(i1,i2,i3) = ipw 1841 end do 1842 1843 g2g1 = 0; nmiss = 0 1844 do ipw=1,npw2 1845 i1 = kg2(1,ipw); if (i1<0) i1=i1+n1; i1=i1+1 1846 i2 = kg2(2,ipw); if (i2<0) i2=i2+n2; i2=i2+1 1847 i3 = kg2(3,ipw); if (i3<0) i3=i3+n3; i3=i3+1 1848 g2g1(ipw) = iwork(i1,i2,i3) 1849 if (g2g1(ipw) == 0) nmiss = nmiss + 1 1850 end do 1851 1852 ABI_FREE(iwork) 1853 1854 end subroutine kg_map
m_gsphere/make_istwk_table [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
make_istwfk_table
FUNCTION
INPUTS
ng1,ng2,ng3
OUTPUT
NOTES
Useful relations:
u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^*
and therefore:
u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
SOURCE
1878 subroutine make_istwfk_table(istwf_k,ng1,ng2,ng3,ig1_inver,ig2_inver,ig3_inver) 1879 1880 !Arguments ------------------------------------ 1881 !scalars 1882 integer,intent(in) :: ng1,ng2,ng3,istwf_k 1883 !arrays 1884 integer,intent(out) :: ig1_inver(ng1),ig2_inver(ng2),ig3_inver(ng3) 1885 1886 !Local variables ------------------------------ 1887 !scalars 1888 integer :: i1,i2,i3 1889 character(len=500) :: msg 1890 1891 !************************************************************************ 1892 1893 ! Initialize the inverse coordinates 1894 select case (istwf_k) 1895 1896 case (1) 1897 ig1_inver(1)=1 1898 do i1=2,ng1 1899 ig1_inver(i1)=ng1+2-i1 1900 end do 1901 ig2_inver(1)=1 1902 do i2=2,ng2 1903 ig2_inver(i2)=ng2+2-i2 1904 end do 1905 ig3_inver(1)=1 1906 do i3=2,ng3 1907 ig3_inver(i3)=ng3+2-i3 1908 end do 1909 1910 case (2:8) 1911 if (istwf_k==2 .or. istwf_k==4 .or. istwf_k==6 .or. istwf_k==8) then 1912 ig1_inver(1)=1 1913 do i1=2,ng1 1914 ig1_inver(i1)=ng1+2-i1 1915 end do 1916 else 1917 do i1=1,ng1 1918 ig1_inver(i1)=ng1+1-i1 1919 end do 1920 end if 1921 if (istwf_k>=2 .and. istwf_k<=5) then 1922 ig2_inver(1)=1 1923 do i2=2,ng2 1924 ig2_inver(i2)=ng2+2-i2 1925 end do 1926 else 1927 do i2=1,ng2 1928 ig2_inver(i2)=ng2+1-i2 1929 end do 1930 end if 1931 if (istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7) then 1932 ig3_inver(1)=1 1933 do i3=2,ng3 1934 ig3_inver(i3)=ng3+2-i3 1935 end do 1936 else 1937 do i3=1,ng3 1938 ig3_inver(i3)=ng3+1-i3 1939 end do 1940 end if 1941 1942 case default 1943 write(msg,'(a,i0)')" Wrong value for istwf_k: ",istwf_k 1944 ABI_ERROR(msg) 1945 end select 1946 1947 end subroutine make_istwfk_table
m_gsphere/merge_and_sort_kg [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
merge_and_sort_kg
FUNCTION
This routine merges a set of k-centered G-spheres of cutoff energy ecut and returns a Gamma-centered G-spheres. The elements in the final G-spheres are packed with increasing module.
INPUTS
nkpt=Number of k-points kptns(3,nkpt)=The k-points in reduced coordinates defining the k-centered G-spheres. ecut=Cutoff energy for the k-centered G-spheres. nsym2=Number of symmetry operations. pinv=-1 if time-reversal can be used, 1 otherwise symrel2(3,3,nsym2)=symmetry operations in real space. gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) prtvol=Flag defining the verbosity level.
SIDE EFFECTS
gbig(:,:)
in input : pointer to NULL
in output: gbig(3,1:npw) contains the set of G-vectors ordered by shell obtained by
merging the k-centered sphere.
shlim_p(:)
in input : pointer to NULL
in output: shlim_p(nbase)=Cumulative number of G-vectors for each shell.
where nbase is the number of irreducible G"s found.
SOURCE
947 subroutine merge_and_sort_kg(nkpt,kptns,ecut,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p) 948 949 !Arguments ------------------------------------ 950 !scalars 951 integer,intent(in) :: nkpt,nsym2,pinv,prtvol 952 real(dp),intent(in) :: ecut 953 !arrays 954 integer,intent(in) :: symrel2(3,3,nsym2) 955 real(dp),intent(in) :: kptns(3,nkpt),gprimd(3,3) 956 integer,pointer :: gbig(:,:) 957 integer,optional,pointer :: shlim_p(:) 958 959 !Local variables------------------------------- 960 !scalars 961 integer,parameter :: mkmem_=1 962 integer :: ikg,ig,ikpt,nbase,sizepw,in,maxpw,is,iinv,ish,ilim,mpw 963 integer :: exchn2n3d,istwf_k,onpw_k,ierr,npw_k,ii,isym,sizeold 964 logical :: found 965 character(len=500) :: msg 966 type(MPI_type) :: MPI_enreg_seq 967 !arrays 968 integer :: gcur(3),geq(3) 969 integer :: symrec2t(3,3,nsym2) 970 integer :: dum_kg(3,0) 971 integer,allocatable :: gbase(:,:),gbasek(:,:,:) 972 integer,allocatable :: gcurr(:,:),gshell(:,:),insort(:),gtmp(:,:) 973 integer,allocatable :: nbasek(:),nshell(:),shlim(:) 974 integer,allocatable :: npwarr(:) 975 real(dp) :: kpoint(3),gmet(3,3) 976 real(dp),allocatable :: cnorm(:),cnormk(:,:),ctmp(:) 977 978 ! ********************************************************************* 979 980 ! * Fake MPI_type for the sequential part. 981 ! This routine should not be parallelized as communicating gbig and other 982 ! tables takes more time than recalculating them in sequential. 983 call initmpi_seq(MPI_enreg_seq) 984 985 !Compute reciprocal space metrics 986 do ii=1,3 987 gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+& 988 & gprimd(2,ii)*gprimd(2,:)+& 989 & gprimd(3,ii)*gprimd(3,:) 990 end do 991 992 ! * Here we use TRANSPOSE(symrel2) instead of the more intuitive symrel2^{-1t} for historical reasons 993 ! It does not affect the results since in the code below we only check the module of G 994 do isym=1,nsym2 995 symrec2t(:,:,isym)=TRANSPOSE(symrel2(:,:,isym)) 996 end do 997 ! 998 ! ============================================== 999 ! ==== Find irreducible G-vectors at each k ==== 1000 ! ============================================== 1001 1002 ABI_MALLOC(npwarr,(nkpt)) 1003 exchn2n3d=0; ikg=0 1004 do ikpt=1,nkpt 1005 kpoint=kptns(:,ikpt); istwf_k=1 1006 call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,dum_kg,kpoint,0,MPI_enreg_seq,0,npwarr(ikpt)) 1007 end do 1008 mpw = MAXVAL(npwarr) 1009 1010 ABI_MALLOC(nbasek,(nkpt)) 1011 ABI_MALLOC(gbasek,(3,mpw,nkpt)) 1012 ABI_MALLOC(cnormk,(mpw,nkpt)) 1013 nbasek=0 ! # of irreducible G at each k. 1014 cnormk=zero ! Norm of each irreducible G. 1015 gbasek=0 ! The set of irreducible G"s at each k. 1016 1017 do ikpt=1,nkpt 1018 1019 kpoint = kptns(:,ikpt) 1020 npw_k = npwarr(ikpt) 1021 1022 exchn2n3d=0; ikg=0; istwf_k=1 1023 ABI_MALLOC(gcurr,(3,npw_k)) 1024 call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,gcurr,kpoint,mkmem_,MPI_enreg_seq,npw_k,onpw_k) 1025 1026 if (ANY(gcurr(:,1)/=0)) then 1027 ABI_BUG("gcurr(:,1)/=0") 1028 end if 1029 ! 1030 ! * Search for the G"s generating the others by symmetry. 1031 ! NB: Here we use symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above 1032 call get_irredg(npw_k,nsym2,pinv,gprimd,symrec2t,gcurr,nbasek(ikpt),gbasek(:,:,ikpt),cnormk(:,ikpt)) 1033 1034 ABI_FREE(gcurr) 1035 end do 1036 ! 1037 ! === Reduce info over k-points === 1038 ! * Here symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above 1039 sizepw=2*mpw 1040 ABI_MALLOC(gbase,(3,sizepw)) 1041 ABI_MALLOC(cnorm,(sizepw)) 1042 nbase=0 ! # of irred G found. 1043 1044 call merge_kgirr(nsym2,pinv,nkpt,mpw,sizepw,symrec2t,nbasek,cnormk,gbasek,nbase,gbase,cnorm,ierr) 1045 if (ierr/=0) then 1046 ABI_ERROR('merge_kgirr returned a non-zero status error') 1047 end if 1048 1049 ABI_FREE(nbasek) 1050 ABI_FREE(cnormk) 1051 ABI_FREE(gbasek) 1052 ! 1053 !=== Reorder base G-vectors in order of increasing module === 1054 ! 1055 !Generate all shells of G-vectors: star of a g==set of all symetrics of this g 1056 ABI_MALLOC(gshell,(3,2*nsym2)) 1057 ABI_MALLOC(shlim,(nbase)) 1058 1059 ABI_MALLOC(gbig,(3,sizepw)) 1060 ! 1061 !TODO 1062 #if 0 1063 !* Here symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above 1064 call getfullg(nbase,nsym2,pinv,sizepw,gbase,symrec2t,cnorm,maxpw,gbig,shlim,ierr) 1065 if (ierr/0) RETURN 1066 1067 #else 1068 ABI_MALLOC(insort,(nbase)) 1069 ABI_MALLOC(nshell,(nbase)) 1070 do in=1,nbase 1071 insort(in)=in 1072 end do 1073 call sort_dp(nbase,cnorm,insort,tol14) 1074 ! 1075 !Loop over all different modules of g''s (=shells): 1076 maxpw=0 1077 do in=1,nbase 1078 nshell(in)=0 1079 gcur(:)=gbase(:,insort(in)) 1080 1081 do is=1,nsym2 ! Loop over all symetries: 1082 do iinv=pinv,1,2 1083 geq(:)=iinv*(symrel2(1,:,is)*gcur(1)+symrel2(2,:,is)*gcur(2)+symrel2(3,:,is)*gcur(3)) 1084 1085 found=.FALSE.; ish=1 1086 do while ((.not.found) .and. (ish<=nshell(in))) ! Search for symetric of g and eventually add it: 1087 found=ALL(geq(:)==gshell(:,ish)) 1088 ish=ish+1 1089 end do 1090 if (.not.found) then 1091 nshell(in)=nshell(in)+1 1092 gshell(:,nshell(in))=geq(:) 1093 end if 1094 end do 1095 end do 1096 1097 if ((maxpw+nshell(in)) > sizepw) then 1098 ! We need to increase the size of the gbase, gbig and cnorm arrays while still keeping their content. 1099 ! This is done using two temporary arrays gtmp and ctmp 1100 ABI_WARNING("Had to reallocate gbase, gbig, cnorm. Perhaps geometry too inaccurate. Possible fix: correct your input file.") 1101 ABI_MALLOC(ctmp,(sizepw)) 1102 ABI_MALLOC(gtmp,(3,sizepw)) 1103 sizeold=sizepw 1104 sizepw=maxpw+nshell(in) 1105 1106 ctmp(:)=cnorm(:) 1107 gtmp(:,:)=gbase(:,:) 1108 1109 ABI_FREE(cnorm) 1110 ABI_MALLOC(cnorm,(sizepw)) 1111 cnorm(1:sizeold)=ctmp(1:sizeold) 1112 cnorm(sizeold+1:sizepw)=zero 1113 ABI_FREE(ctmp) 1114 1115 ! MG why this? gbase should not be changed! 1116 ABI_FREE(gbase) 1117 ABI_MALLOC(gbase,(3,sizepw)) 1118 gbase(:,:sizeold)=gtmp(:,:sizeold) 1119 gbase(:,sizeold+1:sizepw)=0 1120 gtmp(:,:)=gbig(:,:) 1121 1122 ABI_FREE(gbig) 1123 ABI_MALLOC(gbig,(3,sizepw)) 1124 gbig(:,:sizeold)=gtmp(:,:sizeold) 1125 gbig(:,sizeold+1:sizepw)=0 1126 ABI_FREE(gtmp) 1127 end if 1128 ! 1129 ! Store this shell of g''s in a big array of g (gbig): 1130 do ig=1,nshell(in) 1131 gbig(:,ig+maxpw)=gshell(:,ig) 1132 end do 1133 maxpw=maxpw+nshell(in) 1134 1135 end do ! End loop over shells 1136 ! 1137 ! * Compute shell limits 1138 ilim=0 1139 do in=1,nbase 1140 ilim=ilim+nshell(in) 1141 shlim(in)=ilim 1142 end do 1143 1144 if (PRESENT(shlim_p)) then ! Return shlim_p 1145 ABI_MALLOC(shlim_p,(nbase)) 1146 shlim_p = shlim 1147 end if 1148 1149 ! Re-allocate gbig with correct sizes so that caller can inquire the size 1150 ABI_MALLOC(gtmp,(3,ilim)) 1151 gtmp = gbig(:,1:ilim) 1152 ABI_FREE(gbig) 1153 ABI_MALLOC(gbig,(3,ilim)) 1154 gbig=gtmp 1155 ABI_FREE(gtmp) 1156 1157 if (prtvol>10) then ! Print out shell limits 1158 write(msg,'(3a)')& 1159 & ' Shells found:',ch10,& 1160 & ' number of shell number of G vectors cut-off energy [Ha} ' 1161 call wrtout(std_out,msg,'COLL') 1162 do in=1,nbase 1163 write(msg,'(12x,i4,17x,i6,12x,f8.3)')in,shlim(in),2*pi**2*cnorm(in) 1164 call wrtout(std_out,msg,'COLL') 1165 end do 1166 call wrtout(std_out,ch10,'COLL') 1167 end if 1168 1169 ABI_FREE(gshell) 1170 ABI_FREE(insort) 1171 ABI_FREE(nshell) 1172 #endif 1173 1174 call destroy_mpi_enreg(MPI_enreg_seq) 1175 ABI_FREE(gbase) 1176 ABI_FREE(shlim) 1177 ABI_FREE(cnorm) 1178 ABI_FREE(npwarr) 1179 1180 end subroutine merge_and_sort_kg
m_gsphere/merge_kgirr [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
merge_kgirr
FUNCTION
Given a list of irreducible reciprocal vectors associated to different k-centered spheres, this subroutine finds the minimal set of G vectors needed to reconstruct the union of the spheres through symmetry operations.
INPUTS
nsym=number of symmetry operations pinv=-1 if time-reversal can be used, 0 otherwise nkpt=number of k-points for k-centered spheres mpw=Max number of G vectors for each k-point sizepw=Max expected number of G vectors found symrec(3,3,nsym)=symmetries in reciprocal space given in reduced coordinates nbasek(nkpt)=number of irred G for each k-point cnormk(mpw,nkpt)=the norm of each k-centered G (only 1:nbase(ik)) is used gbasek(3,mpw,nkpt)
OUTPUT
nbase=number of irreducible G needed to reconstruct the initial set of spheres gbase(3,sizepw)=irreducible G found in reciprocal coordinates cnorm(sizepw)=Norm of each irred G vector ierr= Exit status, if /=0 the number of G vectors found exceeds sizepw
SOURCE
1444 subroutine merge_kgirr(nsym,pinv,nkpt,mpw,sizepw,symrec,nbasek,cnormk,gbasek,nbase,gbase,cnorm,ierr) 1445 1446 !Arguments ------------------------------------ 1447 !scalars 1448 integer,intent(in) :: mpw,nkpt,nsym,pinv,sizepw 1449 integer,intent(out) :: ierr,nbase 1450 !arrays 1451 integer,intent(in) :: gbasek(3,mpw,nkpt),nbasek(nkpt),symrec(3,3,nsym) 1452 integer,intent(inout) :: gbase(3,sizepw) !vz_i 1453 real(dp),intent(in) :: cnormk(mpw,nkpt) 1454 real(dp),intent(inout) :: cnorm(sizepw) !vz_i 1455 1456 !Local variables------------------------------- 1457 !scalars 1458 integer :: ikpt,inb,irgk,isym 1459 real(dp) :: eps,norm 1460 logical :: found 1461 character(len=500) :: msg 1462 !arrays 1463 integer :: gbas(3),gcur(3),geq(3) 1464 1465 ! ************************************************************************* 1466 1467 DBG_ENTER("COLL") 1468 1469 if (pinv/=1.and.pinv/=-1) then 1470 write(msg,'(a,i6)')' The argument pinv should be -1 or 1, however, pinv =',pinv 1471 ABI_BUG(msg) 1472 end if 1473 ! 1474 ! === Start with zero number of G found === 1475 nbase=0 ; ierr=0 1476 do ikpt=1,nkpt 1477 do irgk=1,nbasek(ikpt) 1478 gcur(:)=gbasek(:,irgk,ikpt) 1479 norm=cnormk(irgk,ikpt) ; eps=tol8*norm 1480 found=.FALSE. ; inb=1 1481 do while ((.not.found).and.(inb<=nbase)) 1482 if (ABS(norm-cnorm(inb))<=eps) then 1483 gbas(:)=gbase(:,inb) 1484 isym=1 1485 do while ((.not.found).and.(isym<=nsym)) 1486 geq(:)=MATMUL(symrec(:,:,isym),gcur) 1487 found=ALL(geq(:)==gbas(:)) 1488 if (pinv==-1) found= (found.or.ALL(geq(:)==-gbas(:)) ) ! For time-reversal 1489 isym=isym+1 1490 end do 1491 end if 1492 inb=inb+1 1493 end do 1494 if (.not.found) then 1495 ! === Add to the list === 1496 nbase=nbase+1 1497 if (nbase>sizepw) then 1498 write(msg,'(2(a,i5),a)')& 1499 & ' nbase (',nbase,') became greater than sizepw = ',sizepw,' returning ierr=1 ' 1500 ABI_WARNING(msg) 1501 ierr=1; RETURN 1502 end if 1503 cnorm(nbase)=cnormk(irgk,ikpt) 1504 gbase(:,nbase)=gcur(:) 1505 end if 1506 end do 1507 end do 1508 1509 DBG_EXIT("COLL") 1510 1511 end subroutine merge_kgirr
m_gsphere/setup_G_rotation [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
setup_G_rotation
FUNCTION
Set up tables indicating rotation of G-vectors.
INPUTS
nsym=Number of symmetry operations. symrec(3,3,nsym)=Symmetry operations in reciprocal space. timrev=2 if time reversal can be used, 1 otherwise. npw=Number of planewaves in the sphere. gvec(3,npw)=Coordinates of plane waves, supposed to be ordered in increasing modulus g2sh(npw)=For each G, it gives the index of the shell to which it belongs. nsh=Number of shells shlim(nsh+1)=Index of the first G vector in each shell, =npw+1 for nsh+1
OUTPUT
grottb (npw,2,nsym)= grottb(G,I,S) is the index of (SI) G in the array gvec.
grottbm1(npw,2,nsym)= index of IS^{-1} G.
NOTES:
I is either the identity or the inversion (time reversal in reciprocal space).
S is one of the symmetry operation in reciprocal space belonging to the Space group.
SOURCE
205 subroutine setup_G_rotation(nsym,symrec,timrev,npw,gvec,g2sh,nsh,shlim,grottb,grottbm1) 206 207 !Arguments ------------------------------------ 208 !scalars 209 integer,intent(in) :: npw,nsh,nsym,timrev 210 !arrays 211 integer,intent(in) :: g2sh(npw),gvec(3,npw),shlim(nsh+1),symrec(3,3,nsym) 212 integer,intent(inout) :: grottb (npw,timrev,nsym) 213 integer,intent(inout) :: grottbm1(npw,timrev,nsym) 214 215 !Local variables ------------------------------ 216 !scalars 217 integer :: ee,ig1,ig2,ish1,isym,itim,ss 218 logical :: found 219 character(len=500) :: msg 220 !arrays 221 integer :: gbase(3),grot(3) 222 223 !************************************************************************ 224 ! 225 ! === Set up G-rotation table === 226 do ig1=1,npw 227 ish1=g2sh(ig1) ; ss=shlim(ish1) ; ee=shlim(ish1+1)-1 228 gbase(:)=gvec(:,ig1) 229 230 do itim=1,timrev 231 do isym=1,nsym 232 grot=(3-2*itim)*MATMUL(symrec(:,:,isym),gbase) 233 found=.FALSE. 234 ! * Loop on the shell of ig1 to speed up the search. 235 do ig2=ss,ee 236 if (ALL(ABS(grot(:)-gvec(:,ig2))==0)) then 237 found=.TRUE. 238 grottb (ig1,itim,isym)=ig2 239 grottbm1(ig2,itim,isym)=ig1 240 end if 241 end do 242 if (.not.found) then 243 write(msg,'(3a,i5,a,i5,1x,2(3i5,a),a,i3,a,i3)')& 244 'G-shell not closed',ch10,& 245 ' Initial G vector ',ig1,'/',npw,gbase(:),' Rotated G vector ',grot(:),ch10,& 246 ' Through sym ',isym,' and itim ',itim 247 ABI_ERROR(msg) 248 end if 249 end do 250 end do 251 252 end do !ig1 253 254 end subroutine setup_G_rotation
m_gsphere/symg [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
symg
FUNCTION
Treat symmetries applied to the G vectors, in view of the application to symmetrization of the dielectric matrix. Generate a list of time-reversed G vectors, as well as a list of spatially-symmetric G vectors.
INPUTS
kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. npwdiel=number of planewaves for the dielectric matrix nsym=number of symmetry symrel(3,3,nsym)=symmetry matrices in real space (integers) tnons(3,nsym)=reduced nonsymmorphic translations (symrel and tnons are in terms of real space primitive translations)
OUTPUT
phdiel(2,npwdiel,nsym)=phase associated with point symmetries applied to G sym_g(npwdiel,nsym)=index list of symmetric G vectors (could save a bit of space by suppressing isym=1, since the corresponding symmetry is the identity) tmrev_g(npwdiel)=index list of inverted G vectors (time-reversed)
SOURCE
2195 subroutine symg(kg_diel,npwdiel,nsym,phdiel,sym_g,symrel,tmrev_g,tnons) 2196 2197 !Arguments ------------------------------------ 2198 !scalars 2199 integer,intent(in) :: npwdiel,nsym 2200 !arrays 2201 integer,intent(in) :: kg_diel(3,npwdiel),symrel(3,3,nsym) 2202 integer,intent(out) :: sym_g(npwdiel,nsym),tmrev_g(npwdiel) 2203 real(dp),intent(in) :: tnons(3,nsym) 2204 real(dp),intent(out) :: phdiel(2,npwdiel,nsym) 2205 2206 !Local variables------------------------------- 2207 !scalars 2208 integer :: g1,g2,g3,ipw,isym,j1,j2,j3,m1m,m1p,m2m,m2p,m3m,m3p,symmg,trevg 2209 real(dp) :: arg,tau1,tau2,tau3 2210 !character(len=500) :: msg 2211 !arrays 2212 integer,allocatable :: grid(:,:,:) 2213 2214 ! ************************************************************************* 2215 2216 !Determines maximal bounds of the zone spanned by the planewaves 2217 m1m=0 ; m2m=0 ; m3m=0 ; m1p=0 ; m2p=0 ; m3p=0 2218 do ipw=1,npwdiel 2219 g1=kg_diel(1,ipw) 2220 g2=kg_diel(2,ipw) 2221 g3=kg_diel(3,ipw) 2222 if(g1<m1m)m1m=g1 ; if(g1>m1p)m1p=g1 2223 if(g2<m2m)m2m=g2 ; if(g2>m2p)m2p=g2 2224 if(g3<m3m)m3m=g3 ; if(g3>m3p)m3p=g3 2225 end do 2226 2227 !Set up grid, that associate to each point the index of the 2228 !corresponding planewave, if there is one 2229 ABI_MALLOC(grid, (m1m:m1p,m2m:m2p,m3m:m3p)) 2230 grid(:,:,:)=0 2231 do ipw=1,npwdiel 2232 g1=kg_diel(1,ipw) 2233 g2=kg_diel(2,ipw) 2234 g3=kg_diel(3,ipw) 2235 grid(g1,g2,g3)=ipw 2236 end do 2237 2238 !Set up tmrev_g and sym_g arrays 2239 do ipw=1,npwdiel 2240 g1=kg_diel(1,ipw) 2241 g2=kg_diel(2,ipw) 2242 g3=kg_diel(3,ipw) 2243 2244 ! Treat first time-reversal symmetry 2245 trevg=grid(-g1,-g2,-g3) 2246 if(trevg==0)then 2247 ABI_BUG('Do not find the time-reversed symmetric of a G-vector.') 2248 end if 2249 tmrev_g(ipw)=trevg 2250 2251 ! Treat now spatial symmetries 2252 do isym=1,nsym 2253 2254 ! Get rotated G vector Gj for each symmetry element 2255 ! -- here we use the TRANSPOSE of symrel; assuming symrel expresses 2256 ! the rotation in real space, the transpose is then appropriate 2257 ! for G space symmetrization (according to Doug : see routine irrzg.f) 2258 j1=symrel(1,1,isym)*g1+& 2259 & symrel(2,1,isym)*g2+symrel(3,1,isym)*g3 2260 j2=symrel(1,2,isym)*g1+& 2261 & symrel(2,2,isym)*g2+symrel(3,2,isym)*g3 2262 j3=symrel(1,3,isym)*g1+& 2263 & symrel(2,3,isym)*g2+symrel(3,3,isym)*g3 2264 symmg=grid(j1,j2,j3) 2265 if(symmg==0)then 2266 ABI_BUG('Do not find the spatially symmetric of a G-vector.') 2267 end if 2268 sym_g(ipw,isym)=symmg 2269 2270 ! Get associated phase 2271 tau1=tnons(1,isym) 2272 tau2=tnons(2,isym) 2273 tau3=tnons(3,isym) 2274 if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then 2275 ! compute exp(-2*Pi*I*G dot tau) using original G 2276 arg=two_pi*(dble(g1)*tau1+dble(g2)*tau2+dble(g3)*tau3) 2277 phdiel(1,ipw,isym)=cos(arg) 2278 phdiel(2,ipw,isym)=-sin(arg) 2279 else 2280 phdiel(1,ipw,isym)=1._dp 2281 phdiel(2,ipw,isym)=0._dp 2282 end if 2283 2284 end do 2285 end do 2286 2287 ABI_FREE(grid) 2288 2289 end subroutine symg
m_gsphere/table_gbig2kg [ Functions ]
[ Top ] [ m_gsphere ] [ Functions ]
NAME
table_gbig2kg
FUNCTION
Associate the kg_k set of g-vectors with the big array of gbig The array gbig(3,maxpw) contains all g-vectors used for all k-points, in order of increasing shells. For a each k-point, the wave-functions are defined only on a particular set of g-vectors kg_k (included in gbig). This set is defined by array gamma2k: The array gamma2k(ig=1,maxpw) translates the index of the gbig (from 1 to maxpw) into the corresponding index in array kg_k. If gbig(ig) does not exist in kg_k, gamma2k(ig) contains npw_k+1.
INPUTS
npw_k=Number of planewaves in the k-centered basis set kg_k(3,npw_k)=The k-centered basis set maxpw=Number of G in gbig gbig(3,maxpw)=The union of the G-spheres at different k-points.
OUTPUT
ierr=Status error. It gives the number of G of kg_k not contained in gbig. gamma2k(maxpw)=Mapping gbig -> kg_k
SOURCE
1976 pure subroutine table_gbig2kg(npw_k,kg_k,maxpw,gbig,gamma2k,ierr) 1977 1978 !Arguments ------------------------------------ 1979 !scalars 1980 integer,intent(in) :: npw_k,maxpw 1981 integer,intent(out) :: ierr 1982 !arrays 1983 integer,intent(in) :: kg_k(3,npw_k) 1984 integer,intent(in) :: gbig(3,maxpw) 1985 integer,intent(out) :: gamma2k(maxpw) 1986 1987 !Local variables------------------------------- 1988 !scalars 1989 integer :: ig,igp 1990 logical :: found 1991 !arrays 1992 integer :: gcur(3) 1993 1994 ! ********************************************************************* 1995 1996 ierr=0 1997 gamma2k(:)=npw_k+1 ! Initialize array gamma2k 1998 1999 do ig=1,npw_k ! Loop over g-vectors, for this k point. 2000 gcur(:)=kg_k(:,ig) 2001 igp=0; found=.FALSE. 2002 do while ((.not.found) .and. igp<maxpw) ! Search selected vector in array gbig: TODO this part can be optimized 2003 igp=igp+1 2004 found=ALL(gcur(:)==gbig(:,igp)) 2005 end do 2006 if (found) then ! Store it if found: 2007 gamma2k(igp)=ig 2008 else 2009 ierr=ierr+1 2010 end if 2011 end do 2012 2013 end subroutine table_gbig2kg