TABLE OF CONTENTS
- ABINIT/dfpt_vlocal
- ABINIT/dfpt_vlocaldq
- ABINIT/dfpt_vlocaldqdq
- ABINIT/dfpt_vmetdqdq
- ABINIT/m_mklocl
- ABINIT/mklocl
- ABINIT/mklocl_recipspace
- ABINIT/vlocalstr
ABINIT/dfpt_vlocal [ Functions ]
NAME
dfpt_vlocal
FUNCTION
Compute local part of 1st-order potential from the appropriate atomic pseudopotential with structure and derivative factor. In case of derivative with respect to k or electric (magnetic Zeeman) field perturbation, the 1st-order local potential vanishes.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX gmet(3,3)=reciprocal space metric (Bohr**-2) gsqcut=cutoff G**2 for included G s in fft box. idir=direction of atomic displacement (=1,2 or 3 : displacement of atom ipert along the 1st, 2nd or 3rd axis). ipert=number of the atom being displaced in the frozen-phonon mpi_enreg=information about MPI parallelization mqgrid=dimension of q grid for pseudopotentials natom=number of atoms in cell. nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ntypat=number of types of atoms in cell. n1,n2,n3=fft grid. ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. qgrid(mqgrid)=grid of q points from 0 to qmax. qphon(3)=wavevector of the phonon ucvol=unit cell volume (Bohr**3). vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom. xred(3,natom)=reduced atomic coordinates
OUTPUT
vpsp1(cplex*nfft)=first-order local crystal pseudopotential in real space (including the minus sign, forgotten in the paper non-linear..
SOURCE
793 subroutine dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir,ipert,& 794 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,& 795 & ntypat,n1,n2,n3,ph1d,qgrid,qphon,ucvol,vlspl,vpsp1,xred) 796 797 !Arguments ------------------------------- 798 !scalars 799 integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,ntypat 800 real(dp),intent(in) :: gsqcut,ucvol 801 type(MPI_type),intent(in) :: mpi_enreg 802 !arrays 803 integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18) 804 real(dp),intent(in) :: gmet(3,3),ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom) 805 real(dp),intent(in) :: qgrid(mqgrid),qphon(3),vlspl(mqgrid,2,ntypat) 806 real(dp),intent(in) :: xred(3,natom) 807 real(dp),intent(out) :: vpsp1(cplex*nfft) 808 809 !Local variables ------------------------- 810 !scalars 811 integer :: i1,i2,i3,ia1,iatom,id1,id2,id3,ig1,ig2,ig3,ii,ii1,im=2 812 integer :: itypat,jj,re=1 813 real(dp),parameter :: tolfix=1.000000001_dp 814 real(dp) :: aa,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,gmag,gq1 815 real(dp) :: gq2,gq3,gsquar,phqim,phqre 816 real(dp) :: qxred2pi,sfi,sfr,vion1,xnorm 817 logical :: qeq0 818 !arrays 819 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 820 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 821 real(dp) :: gq(3) 822 real(dp),allocatable :: work1(:,:) 823 824 ! ********************************************************************* 825 826 iatom=ipert 827 828 if(iatom==natom+1 .or. iatom==natom+2 .or. iatom==natom+10 .or. iatom==natom+11 .or. iatom==natom+5)then 829 830 ! (In case of d/dk or an electric field, or magnetic (Zeeman) field->[natom+5] SPr deb ) 831 vpsp1(1:cplex*nfft)=zero 832 833 else 834 835 ! (In case of a phonon perturbation) 836 ABI_MALLOC(work1,(2,nfft)) 837 work1(1:2,1:nfft)=0.0_dp 838 839 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 840 dqm1=1.0_dp/dq 841 dqdiv6=dq/6.0_dp 842 dq2div6=dq**2/6.0_dp 843 cutoff=gsqcut*tolfix 844 id1=n1/2+2 845 id2=n2/2+2 846 id3=n3/2+2 847 848 ! Get the distrib associated with this fft_grid 849 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 850 851 ! This is to allow q=0 852 qeq0=.false. 853 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)qeq0=.true. 854 855 ! Determination of the atom type 856 ia1=0 857 itypat=0 858 do ii=1,ntypat 859 ia1=ia1+nattyp(ii) 860 if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii 861 end do 862 863 ! Determination of phase qxred* 864 qxred2pi=2.0_dp*pi*(qphon(1)*xred(1,iatom)+ & 865 & qphon(2)*xred(2,iatom)+ & 866 & qphon(3)*xred(3,iatom) ) 867 phqre=cos(qxred2pi) 868 phqim=sin(qxred2pi) 869 ii=0 870 871 do i3=1,n3 872 ig3=i3-(i3/id3)*n3-1 873 gq3=dble(ig3)+qphon(3) 874 gq(3)=gq3 875 do i2=1,n2 876 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 877 ig2=i2-(i2/id2)*n2-1 878 gq2=dble(ig2)+qphon(2) 879 gq(2)=gq2 880 881 ! Note the lower limit of the next loop 882 ii1=1 883 if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then 884 ii1=2 885 ii=ii+1 886 end if 887 do i1=ii1,n1 888 ig1=i1-(i1/id1)*n1-1 889 gq1=dble(ig1)+qphon(1) 890 gq(1)=gq1 891 ii=ii+1 892 gsquar=gsq_vl3(gq1,gq2,gq3) 893 ! Skip G**2 outside cutoff: 894 if (gsquar<=cutoff) then 895 gmag=sqrt(gsquar) 896 897 ! Compute vion(G) for given type of atom 898 jj=1+int(gmag*dqm1) 899 diff=gmag-qgrid(jj) 900 901 ! Evaluate spline fit from q^2 V(q) to get V(q): 902 ! (p. 86 Numerical Recipes, Press et al; NOTE error in book for sign 903 ! of "aa" term in derivative; also see splfit routine. 904 ! This bug fixed here 27 Jan 1992.) 905 906 bb = diff*dqm1 907 aa = 1.0_dp-bb 908 cc = aa*(aa**2-1.0_dp)*dq2div6 909 dd = bb*(bb**2-1.0_dp)*dq2div6 910 vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) + & 911 & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) & 912 & / gsquar 913 914 ! Phase G*xred (complex conjugate) * -i *2pi*(g+q)*vion 915 sfr=-phimag_vl3(ig1,ig2,ig3,iatom)*2.0_dp*pi*gq(idir)*vion1 916 sfi=-phre_vl3(ig1,ig2,ig3,iatom)*2.0_dp*pi*gq(idir)*vion1 917 918 ! Phase q*xred (complex conjugate) 919 work1(re,ii)=sfr*phqre+sfi*phqim 920 work1(im,ii)=-sfr*phqim+sfi*phqre 921 end if 922 923 end do 924 end if 925 end do 926 end do 927 928 ! Transform back to real space 929 call fourdp(cplex,work1,vpsp1,1,mpi_enreg,nfft,1,ngfft,0) 930 931 xnorm=1.0_dp/ucvol 932 vpsp1(1:cplex*nfft)=vpsp1(1:cplex*nfft)*xnorm 933 934 ABI_FREE(work1) 935 936 ! End the condition of non-electric-field 937 end if 938 939 contains 940 941 !Real and imaginary parts of phase. 942 function phr_vl3(x1,y1,x2,y2,x3,y3) 943 944 real(dp) :: phr_vl3 945 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 946 phr_vl3=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 947 end function phr_vl3 948 949 function phi_vl3(x1,y1,x2,y2,x3,y3) 950 951 real(dp) :: phi_vl3 952 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 953 phi_vl3=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 954 end function phi_vl3 955 956 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 957 function ph1_vl3(nri,ig1,ia) 958 959 real(dp) :: ph1_vl3 960 integer,intent(in) :: nri,ig1,ia 961 ph1_vl3=ph1d(nri,ig1+1+n1+(atindx(ia)-1)*(2*n1+1)) 962 end function ph1_vl3 963 964 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 965 function ph2_vl3(nri,ig2,ia) 966 967 real(dp) :: ph2_vl3 968 integer,intent(in) :: nri,ig2,ia 969 ph2_vl3=ph1d(nri,ig2+1+n2+(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1)) 970 end function ph2_vl3 971 972 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 973 function ph3_vl3(nri,ig3,ia) 974 975 real(dp) :: ph3_vl3 976 integer,intent(in) :: nri,ig3,ia 977 ph3_vl3=ph1d(nri,ig3+1+n3+(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 978 end function ph3_vl3 979 980 function phre_vl3(ig1,ig2,ig3,ia) 981 982 real(dp) :: phre_vl3 983 integer,intent(in) :: ig1,ig2,ig3,ia 984 phre_vl3=phr_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),& 985 & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia)) 986 end function phre_vl3 987 988 function phimag_vl3(ig1,ig2,ig3,ia) 989 990 real(dp) :: phimag_vl3 991 integer,intent(in) :: ig1,ig2,ig3,ia 992 phimag_vl3=phi_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),& 993 & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia)) 994 end function phimag_vl3 995 996 function gsq_vl3(g1,g2,g3) 997 998 real(dp) :: gsq_vl3 999 real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions 1000 !Define G^2 based on G space metric gmet. 1001 gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+& 1002 & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+& 1003 & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1) 1004 end function gsq_vl3 1005 1006 end subroutine dfpt_vlocal
ABINIT/dfpt_vlocaldq [ Functions ]
NAME
dfpt_vlocaldq
FUNCTION
Compute q-gradient (at q=0) of the local part of 1st-order atomic displacement potential from the appropriate atomic pseudopotential with structure and derivative factor.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX gmet(3,3)=reciprocal space metric (Bohr**-2) gsqcut=cutoff G**2 for included G s in fft box. idir=direction of atomic displacement (=1,2 or 3 : displacement of atom ipert along the 1st, 2nd or 3rd axis). ipert=number of the atom being displaced in the frozen-phonon mpi_enreg=information about MPI parallelization mqgrid=dimension of q grid for pseudopotentials natom=number of atoms in cell. nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft ntypat=number of types of atoms in cell. n1,n2,n3=fft grid. ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. qdir=direction of the q-gradient qgrid(mqgrid)=grid of q points from 0 to qmax. qphon(3)=wavevector of the phonon ucvol=unit cell volume (Bohr**3). vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom.
OUTPUT
vpsp1dq(cplex*nfft)=q-gradient (at q=0) of the first-order local crystal pseudopotential in real space (including the minus sign, forgotten in the paper non-linear..
NOTES
* IMPORTANT: the formalism followed in this routine assumes a phase factor for the perturbation that is different to the one used elsewhere in the code (See M.Stengel paper): here: e^{i q (R_l + \tau_{\kappa})} rest of ABINIT: e^{i q R_l} **A -i factor has been factorized out in all the contributions of the first q-gradient of the atomic displacement Hamiltonian. This is lately included in the matrix element calculation.
SOURCE
1391 subroutine dfpt_vlocaldq(atindx,cplex,gmet,gsqcut,idir,ipert,& 1392 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,& 1393 & ntypat,n1,n2,n3,ph1d,qdir,qgrid,qphon,ucvol,vlspl,vpsp1dq) 1394 1395 implicit none 1396 1397 !Arguments ------------------------------- 1398 !scalars 1399 integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,ntypat 1400 integer,intent(in) :: qdir 1401 real(dp),intent(in) :: gsqcut,ucvol 1402 type(MPI_type),intent(in) :: mpi_enreg 1403 !arrays 1404 integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18) 1405 real(dp),intent(in) :: gmet(3,3),ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom) 1406 real(dp),intent(in) :: qgrid(mqgrid),qphon(3),vlspl(mqgrid,2,ntypat) 1407 real(dp),intent(out) :: vpsp1dq(cplex*nfft) 1408 1409 !Local variables ------------------------- 1410 !scalars 1411 integer :: i1,i2,i3,ia1,iatom,id1,id2,id3,ig1,ig2,ig3,ii,ii1,im=2 1412 integer :: itypat,re=1 1413 real(dp),parameter :: tolfix=1.000000001_dp 1414 real(dp) :: cutoff,gfact,gmag,gq1 1415 real(dp) :: gq2,gq3,gsquar 1416 real(dp) :: sfi,sfr,xnorm 1417 logical :: qeq0 1418 character(len=500) :: msg 1419 !arrays 1420 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1421 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1422 real(dp) :: gq(3),gvec(3),vion1(1),vion1dq(1) 1423 real(dp),allocatable :: work1(:,:) 1424 1425 ! ********************************************************************* 1426 1427 iatom=ipert 1428 1429 if(iatom==natom+1 .or. iatom==natom+2 .or. iatom==natom+10 .or. iatom==natom+11 .or. iatom==natom+5)then 1430 1431 ! (In case of d/dk or an electric field, or magnetic (Zeeman) field->[natom+5] SPr deb ) 1432 vpsp1dq(1:cplex*nfft)=zero 1433 1434 else 1435 1436 ! (In case of a phonon perturbation) 1437 ABI_MALLOC(work1,(2,nfft)) 1438 work1(1:2,1:nfft)=0.0_dp 1439 1440 cutoff=gsqcut*tolfix 1441 id1=n1/2+2 1442 id2=n2/2+2 1443 id3=n3/2+2 1444 1445 ! Get the distrib associated with this fft_grid 1446 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1447 1448 ! This is to allow q=0 1449 qeq0=.false. 1450 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) then 1451 qeq0=.true. 1452 else 1453 msg='This routine cannot be used for q/=0' 1454 ABI_BUG(msg) 1455 end if 1456 1457 ! Determination of the atom type 1458 ia1=0 1459 itypat=0 1460 do ii=1,ntypat 1461 ia1=ia1+nattyp(ii) 1462 if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii 1463 end do 1464 1465 ii=0 1466 1467 do i3=1,n3 1468 ig3=i3-(i3/id3)*n3-1 1469 gq3=dble(ig3)+qphon(3) 1470 gq(3)=gq3 1471 do i2=1,n2 1472 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 1473 ig2=i2-(i2/id2)*n2-1 1474 gq2=dble(ig2)+qphon(2) 1475 gq(2)=gq2 1476 1477 ! Note the lower limit of the next loop 1478 ii1=1 1479 if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then 1480 ii1=2 1481 ii=ii+1 1482 end if 1483 do i1=ii1,n1 1484 ig1=i1-(i1/id1)*n1-1 1485 gq1=dble(ig1)+qphon(1) 1486 gq(1)=gq1 1487 ii=ii+1 1488 gsquar=gsq_vl3(gq1,gq2,gq3) 1489 ! Skip G**2 outside cutoff: 1490 if (gsquar<=cutoff) then 1491 gmag=sqrt(gsquar) 1492 1493 ! Evaluate spline fit to get V(q) and V(q)': 1494 call splfit(qgrid,vion1dq,vlspl(:,:,itypat),1,(/gmag/),vion1,mqgrid,1) 1495 1496 vion1=vion1/gsquar 1497 vion1dq=(vion1dq-2.0_dp*gmag*vion1)/gsquar 1498 1499 gvec=(/ig1,ig2,ig3/) 1500 gfact=dot_product(gmet(qdir,:),gvec(:))/gmag 1501 1502 ! Phase G*xred (complex conjugate) *2*pi*(g_idir)*vion1dq*gfact 1503 sfr=phre_vl3(ig1,ig2,ig3,iatom)*two_pi*gq(idir)*vion1dq(1)*gfact 1504 sfi=-phimag_vl3(ig1,ig2,ig3,iatom)*two_pi*gq(idir)*vion1dq(1)*gfact 1505 ! Phase G*xred (complex conjugate) *2*pi*(\delta_{idir,qdir})*vion1 1506 if (idir==qdir) then 1507 sfr=sfr + phre_vl3(ig1,ig2,ig3,iatom)*two_pi*vion1(1) 1508 sfi=sfi - phimag_vl3(ig1,ig2,ig3,iatom)*two_pi*vion1(1) 1509 end if 1510 1511 work1(re,ii)=sfr 1512 work1(im,ii)=sfi 1513 end if 1514 1515 end do 1516 end if 1517 end do 1518 end do 1519 1520 ! Transform back to real space 1521 call fourdp(cplex,work1,vpsp1dq,1,mpi_enreg,nfft,1,ngfft,0) 1522 1523 xnorm=1.0_dp/ucvol 1524 vpsp1dq(1:cplex*nfft)=vpsp1dq(1:cplex*nfft)*xnorm 1525 1526 ABI_FREE(work1) 1527 1528 ! End the condition of non-electric-field 1529 end if 1530 1531 contains 1532 1533 !Real and imaginary parts of phase. 1534 function phr_vl3(x1,y1,x2,y2,x3,y3) 1535 real(dp) :: phr_vl3 1536 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 1537 phr_vl3=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 1538 end function phr_vl3 1539 1540 function phi_vl3(x1,y1,x2,y2,x3,y3) 1541 real(dp) :: phi_vl3 1542 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 1543 phi_vl3=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 1544 end function phi_vl3 1545 1546 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 1547 function ph1_vl3(nri,ig1,ia) 1548 real(dp) :: ph1_vl3 1549 integer,intent(in) :: nri,ig1,ia 1550 ph1_vl3=ph1d(nri,ig1+1+n1+(atindx(ia)-1)*(2*n1+1)) 1551 end function ph1_vl3 1552 1553 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 1554 function ph2_vl3(nri,ig2,ia) 1555 real(dp) :: ph2_vl3 1556 integer,intent(in) :: nri,ig2,ia 1557 ph2_vl3=ph1d(nri,ig2+1+n2+(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1)) 1558 end function ph2_vl3 1559 1560 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 1561 function ph3_vl3(nri,ig3,ia) 1562 real(dp) :: ph3_vl3 1563 integer,intent(in) :: nri,ig3,ia 1564 ph3_vl3=ph1d(nri,ig3+1+n3+(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 1565 end function ph3_vl3 1566 1567 function phre_vl3(ig1,ig2,ig3,ia) 1568 real(dp) :: phre_vl3 1569 integer,intent(in) :: ig1,ig2,ig3,ia 1570 phre_vl3=phr_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),& 1571 & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia)) 1572 end function phre_vl3 1573 1574 function phimag_vl3(ig1,ig2,ig3,ia) 1575 real(dp) :: phimag_vl3 1576 integer,intent(in) :: ig1,ig2,ig3,ia 1577 phimag_vl3=phi_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),& 1578 & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia)) 1579 end function phimag_vl3 1580 1581 function gsq_vl3(g1,g2,g3) 1582 real(dp) :: gsq_vl3 1583 real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions 1584 !Define G^2 based on G space metric gmet. 1585 gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+& 1586 & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+& 1587 & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1) 1588 end function gsq_vl3 1589 1590 end subroutine dfpt_vlocaldq
ABINIT/dfpt_vlocaldqdq [ Functions ]
NAME
dfpt_vlocaldqdq
FUNCTION
Compute 2nd q-gradient (at q=0) of the local part of 1st-order atomic displacement potential from the appropriate atomic pseudopotential with structure and derivative factor.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX gmet(3,3)=reciprocal space metric (Bohr**-2) gsqcut=cutoff G**2 for included G s in fft box. idir=direction of atomic displacement (=1,2 or 3 : displacement of atom ipert along the 1st, 2nd or 3rd axis). ipert=number of the atom being displaced in the frozen-phonon mpi_enreg=information about MPI parallelization mqgrid=dimension of q grid for pseudopotentials natom=number of atoms in cell. nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft ntypat=number of types of atoms in cell. n1,n2,n3=fft grid. ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. qdir1=direction of the first q-gradient qdir2=direction of the second q-gradient qgrid(mqgrid)=grid of q points from 0 to qmax. qphon(3)=wavevector of the phonon ucvol=unit cell volume (Bohr**3). vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom.
OUTPUT
vpsp1dqdq(cplex*nfft)=2nd q-gradient (at q=0) of the first-order local crystal pseudopotential in real space
NOTES
* IMPORTANT: the formalism followed in this routine assumes a phase factor for the perturbation that is different to the one used elsewhere in the code (See M.Stengel paper): here: e^{i q (R_l + \tau_{\kappa})} rest of ABINIT: e^{i q R_l} **A -i factor has been factorized out in all the contributions of the second q-gradient of the atomic displacement Hamiltonian. This is lately included in the whole frozen contribution to the q-gradient of the 2nd order energy wrt an atomic displacement and a strain: \Delta E^{\tau_{\kappa\alpha}^* (\beta)}_{m\kvec,\gamma\delta}
SOURCE
1648 subroutine dfpt_vlocaldqdq(atindx,cplex,gmet,gsqcut,idir,ipert,& 1649 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,& 1650 & ntypat,n1,n2,n3,ph1d,qdir1,qdir2,qgrid,qphon,ucvol,vlspl,vpsp1dqdq) 1651 1652 implicit none 1653 1654 !Arguments ------------------------------- 1655 !scalars 1656 integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,ntypat 1657 integer,intent(in) :: qdir1,qdir2 1658 real(dp),intent(in) :: gsqcut,ucvol 1659 type(MPI_type),intent(in) :: mpi_enreg 1660 !arrays 1661 integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18) 1662 real(dp),intent(in) :: gmet(3,3),ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom) 1663 real(dp),intent(in) :: qgrid(mqgrid),qphon(3),vlspl(mqgrid,2,ntypat) 1664 real(dp),intent(out) :: vpsp1dqdq(cplex*nfft) 1665 1666 !Local variables ------------------------- 1667 !scalars 1668 integer :: alpha, delta, gamma 1669 integer :: i1,i2,i3,ia1,iatom,id1,id2,id3,ig1,ig2,ig3,ii,ii1,im=2 1670 integer :: itypat,re=1 1671 real(dp),parameter :: tolfix=1.000000001_dp 1672 real(dp) :: cutoff,delad,delag,gfact,gfact1,gfact2,gmag,gq1 1673 real(dp) :: gq2,gq3,gsquar 1674 real(dp) :: sfi,sfr,term1,term2,xnorm 1675 logical :: qeq0 1676 character(len=500) :: msg 1677 !arrays 1678 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1679 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1680 real(dp) :: gq(3),gvec(3),vion1(1),vion1dq(1),vion1dqdq(1) 1681 real(dp),allocatable :: work1(:,:) 1682 1683 ! ********************************************************************* 1684 1685 iatom=ipert 1686 1687 if(iatom==natom+1 .or. iatom==natom+2 .or. iatom==natom+10 .or. iatom==natom+11 .or. iatom==natom+5)then 1688 1689 ! (In case of d/dk or an electric field, or magnetic (Zeeman) field->[natom+5] SPr deb ) 1690 vpsp1dqdq(1:cplex*nfft)=zero 1691 1692 else 1693 1694 alpha=idir; delta=qdir2; gamma=qdir1 1695 1696 !Kronecker deltas 1697 delad=0.0_dp; delag=0.0_dp 1698 if (alpha==delta) delad=1.0_dp 1699 if (alpha==gamma) delag=1.0_dp 1700 1701 ! (In case of a phonon perturbation) 1702 ABI_MALLOC(work1,(2,nfft)) 1703 work1(1:2,1:nfft)=0.0_dp 1704 1705 cutoff=gsqcut*tolfix 1706 id1=n1/2+2 1707 id2=n2/2+2 1708 id3=n3/2+2 1709 1710 ! Get the distrib associated with this fft_grid 1711 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1712 1713 ! This is to allow q=0 1714 qeq0=.false. 1715 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) then 1716 qeq0=.true. 1717 else 1718 msg='This routine cannot be used for q/=0' 1719 ABI_BUG(msg) 1720 end if 1721 1722 ! Determination of the atom type 1723 ia1=0 1724 itypat=0 1725 do ii=1,ntypat 1726 ia1=ia1+nattyp(ii) 1727 if(atindx(iatom)<=ia1.and.itypat==0)itypat=ii 1728 end do 1729 1730 ii=0 1731 1732 do i3=1,n3 1733 ig3=i3-(i3/id3)*n3-1 1734 gq3=dble(ig3)+qphon(3) 1735 gq(3)=gq3 1736 do i2=1,n2 1737 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 1738 ig2=i2-(i2/id2)*n2-1 1739 gq2=dble(ig2)+qphon(2) 1740 gq(2)=gq2 1741 1742 ! Note the lower limit of the next loop 1743 ii1=1 1744 if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then 1745 ii1=2 1746 ii=ii+1 1747 end if 1748 do i1=ii1,n1 1749 ig1=i1-(i1/id1)*n1-1 1750 gq1=dble(ig1)+qphon(1) 1751 gq(1)=gq1 1752 ii=ii+1 1753 gsquar=gsq_vl3(gq1,gq2,gq3) 1754 ! Skip G**2 outside cutoff: 1755 if (gsquar<=cutoff) then 1756 gmag=sqrt(gsquar) 1757 1758 ! Evaluate spline fit to get first V(q) and V(q)' and later V(q)'': 1759 call splfit(qgrid,vion1dq,vlspl(:,:,itypat),1,(/gmag/),vion1,mqgrid,1) 1760 vion1=vion1/gsquar 1761 vion1dq=(vion1dq-2.0_dp*gmag*vion1)/gsquar 1762 1763 call splfit(qgrid,vion1dqdq,vlspl(:,:,itypat),2,(/gmag/),vion1,mqgrid,1) 1764 vion1dqdq=(vion1dqdq-4.0_dp*gmag*vion1dq-2.0_dp*vion1)/gsquar 1765 1766 gvec=(/ig1,ig2,ig3/) 1767 gfact1=dot_product(gmet(gamma,:),gvec(:)) 1768 gfact2=dot_product(gmet(delta,:),gvec(:)) 1769 gfact=gvec(alpha)*gfact1*gfact2/gsquar 1770 1771 term1=delag*gfact2+delad*gfact1+gvec(alpha)*gmet(gamma,delta) 1772 term1=term1-gfact 1773 term1=term1*vion1dq(1)/gmag 1774 1775 term2=vion1dqdq(1)*gfact 1776 1777 ! structure factors 1778 sfr=phre_vl3(ig1,ig2,ig3,iatom) 1779 sfi=-phimag_vl3(ig1,ig2,ig3,iatom) 1780 1781 ! Multiply structure factor times vion derivatives: 1782 work1(re,ii)=sfr*(term1+term2)*two_pi 1783 work1(im,ii)=sfi*(term1+term2)*two_pi 1784 1785 end if 1786 1787 end do 1788 end if 1789 end do 1790 end do 1791 1792 ! Transform back to real space 1793 call fourdp(cplex,work1,vpsp1dqdq,1,mpi_enreg,nfft,1,ngfft,0) 1794 1795 xnorm=1.0_dp/ucvol 1796 vpsp1dqdq(1:cplex*nfft)=vpsp1dqdq(1:cplex*nfft)*xnorm 1797 1798 ABI_FREE(work1) 1799 1800 ! End the condition of non-electric-field 1801 end if 1802 1803 contains 1804 1805 !Real and imaginary parts of phase. 1806 function phr_vl3(x1,y1,x2,y2,x3,y3) 1807 real(dp) :: phr_vl3 1808 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 1809 phr_vl3=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 1810 end function phr_vl3 1811 1812 function phi_vl3(x1,y1,x2,y2,x3,y3) 1813 real(dp) :: phi_vl3 1814 real(dp),intent(in) :: x1,x2,x3,y1,y2,y3 1815 phi_vl3=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 1816 end function phi_vl3 1817 1818 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 1819 function ph1_vl3(nri,ig1,ia) 1820 real(dp) :: ph1_vl3 1821 integer,intent(in) :: nri,ig1,ia 1822 ph1_vl3=ph1d(nri,ig1+1+n1+(atindx(ia)-1)*(2*n1+1)) 1823 end function ph1_vl3 1824 1825 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 1826 function ph2_vl3(nri,ig2,ia) 1827 real(dp) :: ph2_vl3 1828 integer,intent(in) :: nri,ig2,ia 1829 ph2_vl3=ph1d(nri,ig2+1+n2+(atindx(ia)-1)*(2*n2+1)+natom*(2*n1+1)) 1830 end function ph2_vl3 1831 1832 ! Warning : this function differ from similar ones for ground-state calculations : note the atindx !! 1833 function ph3_vl3(nri,ig3,ia) 1834 real(dp) :: ph3_vl3 1835 integer,intent(in) :: nri,ig3,ia 1836 ph3_vl3=ph1d(nri,ig3+1+n3+(atindx(ia)-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 1837 end function ph3_vl3 1838 1839 function phre_vl3(ig1,ig2,ig3,ia) 1840 real(dp) :: phre_vl3 1841 integer,intent(in) :: ig1,ig2,ig3,ia 1842 phre_vl3=phr_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),& 1843 & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia)) 1844 end function phre_vl3 1845 1846 function phimag_vl3(ig1,ig2,ig3,ia) 1847 real(dp) :: phimag_vl3 1848 integer,intent(in) :: ig1,ig2,ig3,ia 1849 phimag_vl3=phi_vl3(ph1_vl3(re,ig1,ia),ph1_vl3(im,ig1,ia),& 1850 & ph2_vl3(re,ig2,ia),ph2_vl3(im,ig2,ia),ph3_vl3(re,ig3,ia),ph3_vl3(im,ig3,ia)) 1851 end function phimag_vl3 1852 1853 function gsq_vl3(g1,g2,g3) 1854 real(dp) :: gsq_vl3 1855 real(dp),intent(in) :: g1,g2,g3 ! Note that they are real, unlike in other similar function definitions 1856 !Define G^2 based on G space metric gmet. 1857 gsq_vl3=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+& 1858 & g3*g3*gmet(3,3)+2.0_dp*g1*g2*gmet(1,2)+& 1859 & 2.0_dp*g2*g3*gmet(2,3)+2.0_dp*g3*g1*gmet(3,1) 1860 end function gsq_vl3 1861 1862 end subroutine dfpt_vlocaldqdq
ABINIT/dfpt_vmetdqdq [ Functions ]
NAME
dfpt_vmetdqdq
FUNCTION
Compute second q-gradient (at q=0) of the local part of 1st-order metric potential from the appropriate atomic pseudopotential with structure and derivative factor. Additionaly, compute the second q-gradient (at q=0) of the Hartree and XC (if GGA) potentials of the metric perturbation. Cartesian coordinates are employed to define the direction of the metric perturbation and the two q-gradients.
INPUTS
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX gmet(3,3)=reciprocal space metric (Bohr**-2) gsqcut=cutoff G**2 for included G s in fft box. gprimd(3,3)=dimensional reciprocal space primitive translations idir= strain perturbation direction ipert=number of the atom being displaced in the frozen-phonon kxc(nfft,nkxc)=exchange and correlation kernel mpi_enreg=information about MPI parallelization mqgrid=dimension of q grid for pseudopotentials natom=number of atoms in cell. nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed. nspden=number of spin-density components ntypat=number of types of atoms in cell. n1,n2,n3=fft grid. opthartdqdq= if 1 activates the calculation 2nd q-gradient of the Hartree potential ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. qdir=direction of the q-gradient qgrid(mqgrid)=grid of q points from 0 to qmax. qphon(3)=wavevector of the phonon rhog(2,nfft)=array for Fourier transform of GS electron density rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3. ucvol=unit cell volume (Bohr**3). vlspl(mqgrid,2,ntypat)=spline fit of q^2 V(q) for each type of atom.
OUTPUT
vhart1dqdq(cplex*nfft)=2nd q-gradient (at q=0) of the GS density Hartree potential from the metric perturbation vpsp1dqdq(cplex*nfft)=2nd q-gradient (at q=0) of the first-order metric local crystal pseudopotential in real space vxc1dqdq(cplex*nfft)=2nd q-gradient (at q=0) of the GS density XC potential from the metric perturbation (only finite if GGA)
NOTES
** IMPORTANT: the formalism followed in this routine
assumes a phase factor for the perturbation that
is different to the one used elsewhere in the code (See M.Stengel paper):
here: e^{i q (R_l + \tau_{\kappa})}
rest of ABINIT: e^{i q R_l}
**Since the 2nd derivative w.r.t q-vector is calculated along cartesian
directions, the 1/twopi**2 factor (that in the rest of the code is applied
in the reduced to cartesian derivative conversion process) is here
explicictly included in the formulas.
**Notice that idir=1-9, in contrast to the strain perturbation (idir=1-6),
because this term is not symmetric w.r.t permutations of the two strain
indices.
**A -i factor has been factorized out in all the contributions of the second
q-gradient of the metric Hamiltonian. This is lately included in the contribution
of the corresponing term (T4) to the flexoelectric tensor in dfpt_flexoout.F90
SOURCE
1937 subroutine dfpt_vmetdqdq(cplex,gmet,gprimd,gsqcut,idir,ipert,& 1938 & kxc,mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,& 1939 & ntypat,n1,n2,n3,nkxc,nspden,opthartdqdq,ph1d,qdir,qgrid,qphon,rhog,rhor,& 1940 & ucvol,vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq) 1941 1942 implicit none 1943 1944 !Arguments ------------------------------- 1945 !scalars 1946 integer,intent(in) :: cplex,idir,ipert,mqgrid,n1,n2,n3,natom,nfft,nkxc,ntypat 1947 integer,intent(in) :: nspden,opthartdqdq,qdir 1948 real(dp),intent(in) :: gsqcut,ucvol 1949 type(MPI_type),intent(in) :: mpi_enreg 1950 !arrays 1951 integer,intent(in) :: nattyp(ntypat),ngfft(18) 1952 real(dp),intent(in) :: gmet(3,3),gprimd(3,3), kxc(nfft,nkxc) 1953 real(dp),intent(in) :: ph1d(2,(2*n1+1+2*n2+1+2*n3+1)*natom) 1954 real(dp),intent(in) :: qgrid(mqgrid),qphon(3),rhog(2,nfft),rhor(cplex*nfft,nspden) 1955 real(dp),intent(in) :: vlspl(mqgrid,2,ntypat) 1956 real(dp),intent(out) :: vhart1dqdq(cplex*nfft),vpsp1dqdq(cplex*nfft) 1957 real(dp),intent(out) :: vxc1dqdq(cplex*nfft) 1958 1959 !Local variables ------------------------- 1960 !scalars 1961 integer :: beta, delta, gamma 1962 integer :: ia,i1,i2,i3,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,ii1 1963 integer :: itypat,jj 1964 integer, parameter :: im=2, re=1 1965 real(dp),parameter :: tolfix=1.000000001_dp 1966 real(dp) :: cutoff,delbd,delbg,deldg,gfact,gmag,gq1 1967 real(dp) :: gq2,gq3,gsquar,pisqrinv 1968 real(dp) :: sfi,sfr,term1,term2,tmpre,tmpim,uogsquar,work1re,xnorm 1969 logical :: qeq0 1970 character(len=500) :: msg 1971 !arrays 1972 integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 1973 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1974 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1975 real(dp) :: gq(3),gqc(3),vion1(1),vion1dq(1),vion1dqdq(1) 1976 real(dp),allocatable :: work1(:,:) 1977 1978 ! ********************************************************************* 1979 1980 1981 if(ipert/=natom+3 .and. ipert/=natom+4)then 1982 1983 vpsp1dqdq(1:cplex*nfft)=zero 1984 1985 else 1986 1987 beta=idx(2*idir-1); delta=idx(2*idir); gamma=qdir 1988 1989 !Kronecker deltas 1990 delbd=0.0_dp; delbg=0.0_dp; deldg=0.0_dp 1991 if (beta==delta) delbd=1.0_dp 1992 if (beta==gamma) delbg=1.0_dp 1993 if (delta==gamma) deldg=1.0_dp 1994 1995 ABI_MALLOC(work1,(2,nfft)) 1996 work1(1:2,1:nfft)=0.0_dp 1997 1998 cutoff=gsqcut*tolfix 1999 id1=n1/2+2 2000 id2=n2/2+2 2001 id3=n3/2+2 2002 2003 !Get the distrib associated with this fft_grid 2004 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 2005 2006 ! This is to allow q=0 2007 qeq0=.false. 2008 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) then 2009 qeq0=.true. 2010 else 2011 msg='This routine cannot be used for q/=0' 2012 ABI_BUG(msg) 2013 end if 2014 2015 ia1=1 2016 do itypat=1,ntypat 2017 ! ia1,ia2 sets range of loop over atoms: 2018 ia2=ia1+nattyp(itypat)-1 2019 2020 ii=0 2021 2022 do i3=1,n3 2023 ig3=i3-(i3/id3)*n3-1 2024 gq3=dble(ig3)+qphon(3) 2025 gq(3)=gq3 2026 do i2=1,n2 2027 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 2028 ig2=i2-(i2/id2)*n2-1 2029 gq2=dble(ig2)+qphon(2) 2030 gq(2)=gq2 2031 2032 ! Note the lower limit of the next loop 2033 ii1=1 2034 if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then 2035 ii1=2 2036 ii=ii+1 2037 end if 2038 do i1=ii1,n1 2039 ig1=i1-(i1/id1)*n1-1 2040 gq1=dble(ig1)+qphon(1) 2041 gq(1)=gq1 2042 ii=ii+1 2043 gsquar=gsq_vl(ig1,ig2,ig3) 2044 ! Skip G**2 outside cutoff: 2045 if (gsquar<=cutoff) then 2046 gmag=sqrt(gsquar) 2047 2048 ! Obtain G in cartesian coordinates 2049 gqc(1)=gprimd(1,1)*gq(1)+gprimd(1,2)*gq(2)+gprimd(1,3)*gq(3) 2050 gqc(2)=gprimd(2,1)*gq(1)+gprimd(2,2)*gq(2)+gprimd(2,3)*gq(3) 2051 gqc(3)=gprimd(3,1)*gq(1)+gprimd(3,2)*gq(2)+gprimd(3,3)*gq(3) 2052 2053 ! Evaluate spline fit to get first V(q) and V(q)' and later V(q)'': 2054 call splfit(qgrid,vion1dq,vlspl(:,:,itypat),1,(/gmag/),vion1,mqgrid,1) 2055 vion1=vion1/gsquar 2056 vion1dq=(vion1dq-2.0_dp*gmag*vion1)/gsquar 2057 2058 call splfit(qgrid,vion1dqdq,vlspl(:,:,itypat),2,(/gmag/),vion1,mqgrid,1) 2059 vion1dqdq=(vion1dqdq-4.0_dp*gmag*vion1dq-2.0_dp*vion1)/gsquar 2060 2061 ! Assemble structure factor over all atoms of given type: 2062 sfr=0.0_dp 2063 sfi=0.0_dp 2064 do ia=ia1,ia2 2065 sfr=sfr+phre_vl(ig1,ig2,ig3,ia) 2066 sfi=sfi-phimag_vl(ig1,ig2,ig3,ia) 2067 end do 2068 2069 gfact=gqc(beta)*gqc(delta)*gqc(gamma)/gsquar 2070 2071 term1=delbd*gqc(gamma)+delbg*gqc(delta)+deldg*gqc(beta) 2072 term1=term1-gfact 2073 term1=term1*vion1dq(1)/gmag 2074 2075 term2=vion1dqdq(1)*gfact 2076 2077 ! Multiply structure factor times vion derivatives: 2078 work1(re,ii)=work1(re,ii)+sfr*(term1+term2) 2079 work1(im,ii)=work1(im,ii)+sfi*(term1+term2) 2080 2081 ! End skip G**2 outside cutoff: 2082 end if 2083 2084 end do 2085 end if 2086 end do 2087 end do 2088 2089 ia1=ia2+1 2090 2091 ! End loop on type of atoms 2092 end do 2093 2094 ! Set Vloc(G=0)=0: 2095 work1(re,1)=0.0_dp 2096 work1(im,1)=0.0_dp 2097 2098 ! Transform back to real space 2099 call fourdp(cplex,work1,vpsp1dqdq,1,mpi_enreg,nfft,1,ngfft,0) 2100 2101 xnorm=1.0_dp/ucvol/two_pi 2102 vpsp1dqdq(1:cplex*nfft)=vpsp1dqdq(1:cplex*nfft)*xnorm 2103 2104 work1=0.0_dp 2105 2106 ! Calculate the GS density Hartree contribution 2107 if (opthartdqdq==1) then 2108 2109 pisqrinv=1.0_dp/pi**2 2110 2111 ii=0 2112 do i3=1,n3 2113 ig3=i3-(i3/id3)*n3-1 2114 gq3=dble(ig3)+qphon(3) 2115 gq(3)=gq3 2116 do i2=1,n2 2117 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 2118 ig2=i2-(i2/id2)*n2-1 2119 gq2=dble(ig2)+qphon(2) 2120 gq(2)=gq2 2121 2122 ! Note the lower limit of the next loop 2123 ii1=1 2124 if(i3==1 .and. i2==1 .and. qeq0 .and. ig2==0 .and. ig3==0)then 2125 ii1=2 2126 ii=ii+1 2127 end if 2128 do i1=ii1,n1 2129 ig1=i1-(i1/id1)*n1-1 2130 gq1=dble(ig1)+qphon(1) 2131 gq(1)=gq1 2132 ii=ii+1 2133 gsquar=gsq_vl(ig1,ig2,ig3) 2134 ! Skip G**2 outside cutoff: 2135 if (gsquar<=cutoff) then 2136 2137 ! Precalculate quotient of G powers 2138 uogsquar= 1.0_dp/gsquar 2139 2140 ! Obtain G in cartesian coordinates 2141 gqc(1)=gprimd(1,1)*gq(1)+gprimd(1,2)*gq(2)+gprimd(1,3)*gq(3) 2142 gqc(2)=gprimd(2,1)*gq(1)+gprimd(2,2)*gq(2)+gprimd(2,3)*gq(3) 2143 gqc(3)=gprimd(3,1)*gq(1)+gprimd(3,2)*gq(2)+gprimd(3,3)*gq(3) 2144 2145 term1=4.0_dp*gqc(beta)*gqc(gamma)*gqc(delta)*uogsquar*uogsquar 2146 term2=delbd*gqc(gamma)+delbg*gqc(delta)+deldg*gqc(beta) 2147 term2=-term2*uogsquar 2148 2149 work1re=pisqrinv*uogsquar*(term1+term2) 2150 work1(re,ii)=rhog(re,ii)*work1re 2151 work1(im,ii)=rhog(im,ii)*work1re 2152 2153 ! End skip G**2 outside cutoff: 2154 end if 2155 2156 end do 2157 end if 2158 end do 2159 end do 2160 2161 ! Set V(G=0)=0: 2162 work1(re,1)=0.0_dp 2163 work1(im,1)=0.0_dp 2164 2165 ! Transform back to real space 2166 call fourdp(cplex,work1,vhart1dqdq,1,mpi_enreg,nfft,1,ngfft,0) 2167 2168 ! End the calculation of the Hartree contribution 2169 end if 2170 2171 ABI_FREE(work1) 2172 2173 ! Calculate the GS density XC contribution (if GGA) 2174 vxc1dqdq(:)=zero 2175 if (nkxc == 7) then 2176 call dfpt_mkvxcgga_n0met(beta,1,delta,gamma,gprimd,kxc,mpi_enreg, & 2177 & nfft,ngfft,nkxc,nspden,rhor,vxc1dqdq) 2178 2179 !Fictitious i factor temporarily applied. 2180 !It is later canceled by the (-i) factor of the total matrix element 2181 do ii=1,nfft 2182 jj=ii*2 2183 tmpre=vxc1dqdq(jj-1); tmpim=vxc1dqdq(jj) 2184 vxc1dqdq(jj-1)=-tmpim; vxc1dqdq(jj)=tmpre 2185 end do 2186 end if 2187 2188 !End the condition of non-electric-field 2189 end if 2190 2191 contains 2192 2193 !Real and imaginary parts of phase. 2194 function phr_vl(x1,y1,x2,y2,x3,y3) 2195 real(dp) :: phr_vl,x1,x2,x3,y1,y2,y3 2196 phr_vl=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 2197 end function phr_vl 2198 2199 function phi_vl(x1,y1,x2,y2,x3,y3) 2200 real(dp):: phi_vl,x1,x2,x3,y1,y2,y3 2201 phi_vl=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 2202 end function phi_vl 2203 2204 function ph1_vl(nri,ig1,ia) 2205 real(dp):: ph1_vl 2206 integer :: nri,ig1,ia 2207 ph1_vl=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1)) 2208 end function ph1_vl 2209 2210 function ph2_vl(nri,ig2,ia) 2211 real(dp):: ph2_vl 2212 integer :: nri,ig2,ia 2213 ph2_vl=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)) 2214 end function ph2_vl 2215 2216 function ph3_vl(nri,ig3,ia) 2217 real(dp):: ph3_vl 2218 integer :: nri,ig3,ia 2219 ph3_vl=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 2220 end function ph3_vl 2221 2222 function phre_vl(ig1,ig2,ig3,ia) 2223 real(dp):: phre_vl 2224 integer :: ig1,ig2,ig3,ia 2225 phre_vl=phr_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),& 2226 & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia)) 2227 end function phre_vl 2228 2229 function phimag_vl(ig1,ig2,ig3,ia) 2230 real(dp) :: phimag_vl 2231 integer :: ig1,ig2,ig3,ia 2232 phimag_vl=phi_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),& 2233 & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia)) 2234 end function phimag_vl 2235 2236 function gsq_vl(i1,i2,i3) 2237 real(dp) :: gsq_vl 2238 integer :: i1,i2,i3 2239 !Define G^2 based on G space metric gmet. 2240 gsq_vl=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 2241 & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 2242 & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 2243 end function gsq_vl 2244 2245 end subroutine dfpt_vmetdqdq 2246 2247 end module m_mklocl
ABINIT/m_mklocl [ Modules ]
NAME
m_mklocl
FUNCTION
Routines related to the local part of the pseudopotentials.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, MM, DRH) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_mklocl 23 24 use defs_basis 25 use defs_wvltypes 26 use m_abicore 27 use m_errors 28 use m_xmpi 29 use m_dtset 30 31 use defs_datatypes, only : pseudopotential_type 32 use defs_abitypes, only : MPI_type 33 use m_time, only : timab 34 use m_geometry, only : xred2xcart 35 use m_mpinfo, only : ptabs_fourdp 36 use m_pawtab, only : pawtab_type 37 use m_mklocl_realspace, only : mklocl_realspace, mklocl_wavelets 38 use m_fft, only : fourdp 39 use m_gtermcutoff,only : termcutoff 40 41 use m_splines, only : splfit 42 use m_dfpt_mkvxc, only : dfpt_mkvxcgga_n0met 43 44 #if defined HAVE_BIGDFT 45 use BigDFT_API, only : ELECTRONIC_DENSITY 46 use m_abi2big, only : wvl_rho_abi2big 47 #endif 48 49 implicit none 50 51 private
ABINIT/mklocl [ Functions ]
NAME
mklocl
FUNCTION
This method is a wrapper for mklocl_recipspace and mklocl_realspace. It does some consistency checks before calling one of the two methods. Optionally compute : option=1 : local ionic potential throughout unit cell option=2 : contribution of local ionic potential to E gradient wrt xred option=3 : contribution of local ionic potential to stress tensor (only with reciprocal space computations) option=4 : contribution of local ionic potential to second derivative of E wrt xred (only with reciprocal space computations)
INPUTS
if(option==3) eei=local pseudopotential part of total energy (hartree) gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$). gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere). mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization natom=number of atoms in unit cell. nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nspden=number of spin-density components ntypat=number of types of atoms. option= (see above) pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. psps <type(pseudopotential_type)>=variables related to pseudopotentials qprtrb(3)= integer wavevector of possible perturbing potential in basis of reciprocal lattice translations rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$) (needed if option==2 or if option==4) rhor(nfft,nspden)=electron density in electrons/bohr**3. (needed if option==2 or if option==4) rprimd(3,3)=dimensional primitive translations in real space (bohr) ucvol=unit cell volume ($\textrm{Bohr}^3$). vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero, perturbing potential is added of the form $V(G)=(vprtrb(1)+I*vprtrb(2))/2$ at the values G=qprtrb and $(vprtrb(1)-I*vprtrb(2))/2$ at $G=-qprtrb$ (integers) xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
(if option==1) vpsp(nfft)=local crystal pseudopotential in real space. (if option==2) grtn(3,natom)=grads of Etot wrt tn. (if option==3) lpsstr(6)=components of local psp part of stress tensor (Cartesian coordinates, symmetric tensor) in hartree/$\textrm{bohr}^3$ Store 6 unique components in order 11, 22, 33, 32, 31, 21 (if option==4) dyfrlo(3,3,natom)=d2 Eei/d tn(i)/d tn(j). (Hartrees)
NOTES
Note that the present routine is tightly connected to the dfpt_vlocal.f routine, that compute the derivative of the local ionic potential with respect to one atomic displacement. The argument list and the internal loops to be considered were sufficiently different as to make the two routine different.
SOURCE
130 subroutine mklocl(dtset, dyfrlo,eei,gmet,gprimd,grtn,gsqcut,lpsstr,mgfft,& 131 & mpi_enreg,natom,nattyp,nfft,ngfft,nspden,ntypat,option,pawtab,ph1d,psps,qprtrb,& 132 & rhog,rhor,rprimd,ucvol,vprtrb,vpsp,wvl,wvl_den,xred) 133 134 !Arguments ------------------------------------ 135 !scalars 136 integer,intent(in) :: mgfft,natom,nfft,nspden,ntypat,option 137 real(dp),intent(in) :: eei,gsqcut,ucvol 138 type(MPI_type),intent(in) :: mpi_enreg 139 type(dataset_type),intent(in) :: dtset 140 type(pseudopotential_type),intent(in) :: psps 141 type(wvl_internal_type), intent(in) :: wvl 142 type(wvl_denspot_type), intent(inout) :: wvl_den 143 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 144 !arrays 145 integer,intent(in) :: nattyp(ntypat),ngfft(18),qprtrb(3) 146 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom) 147 real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3) 148 real(dp),intent(in) :: vprtrb(2),xred(3,natom) 149 real(dp),intent(in),target :: rhor(nfft,nspden) 150 real(dp),intent(out) :: dyfrlo(3,3,natom),grtn(3,natom),lpsstr(6) 151 real(dp),intent(inout) :: vpsp(nfft) 152 153 !Local variables------------------------------- 154 !scalars 155 character(len=500) :: message 156 !arrays 157 real(dp),allocatable :: xcart(:,:) 158 #if defined HAVE_BIGDFT 159 real(dp),pointer :: rhor_ptr(:,:) 160 #endif 161 162 ! ************************************************************************* 163 164 if (option < 1 .or. option > 4) then 165 write(message,'(a,i0,a,a)')& 166 & 'From the calling routine, option=',option,ch10,& 167 & 'The only allowed values are between 1 and 4.' 168 ABI_ERROR(message) 169 end if 170 if (option > 2 .and. .not.psps%vlspl_recipSpace) then 171 write(message,'(a,i0,a,a,a,a)')& 172 & 'From the calling routine, option=',option,ch10,& 173 & 'but the local part of the pseudo-potential is in real space.',ch10,& 174 & 'Action: set icoulomb = 0 to turn-off real space computations.' 175 ABI_ERROR(message) 176 end if 177 if (option > 2 .and. dtset%usewvl == 1) then 178 write(message,'(a,i0,a,a)')& 179 & 'From the calling routine, option=',option,ch10,& 180 & 'but this is not implemented yet from wavelets.' 181 ABI_ERROR(message) 182 end if 183 184 if (dtset%usewvl == 0) then 185 ! Plane wave case 186 if (psps%vlspl_recipSpace) then 187 call mklocl_recipspace(dyfrlo,eei,gmet,gprimd,grtn,gsqcut,& 188 & dtset%icutcoul,lpsstr,mgfft,mpi_enreg,psps%mqgrid_vl,natom,nattyp, & 189 & nfft,ngfft,dtset%nkpt,ntypat,option,ph1d,psps%qgrid_vl,qprtrb,dtset%rcut,& 190 & rhog,rprimd,ucvol,dtset%vcutgeo,psps%vlspl,vprtrb,vpsp) 191 else 192 call mklocl_realspace(grtn,dtset%icoulomb,mpi_enreg,natom,nattyp,nfft, & 193 & ngfft,dtset%nscforder,nspden,ntypat,option,pawtab,psps,rhog,rhor, & 194 & rprimd,dtset%typat,ucvol,dtset%usewvl,vpsp,xred) 195 end if 196 else 197 ! Store xcart for each atom 198 ABI_MALLOC(xcart,(3, dtset%natom)) 199 call xred2xcart(dtset%natom, rprimd, xcart, xred) 200 ! Eventually retrieve density 201 #if defined HAVE_BIGDFT 202 if (option>1.and.wvl_den%denspot%rhov_is/=ELECTRONIC_DENSITY) then 203 rhor_ptr => rhor ! Just to bypass intent(inout) 204 call wvl_rho_abi2big(1,rhor_ptr,wvl_den) 205 end if 206 #endif 207 ! Wavelets case 208 call mklocl_wavelets(dtset%efield, grtn, mpi_enreg, dtset%natom, & 209 & nfft, nspden, option, rprimd, vpsp, & 210 & wvl_den, wvl, xcart) 211 ABI_FREE(xcart) 212 end if 213 214 end subroutine mklocl
ABINIT/mklocl_recipspace [ Functions ]
NAME
mklocl_recipspace
FUNCTION
Optionally compute : option=1 : local ionic potential throughout unit cell option=2 : contribution of local ionic potential to E gradient wrt xred option=3 : contribution of local ionic potential to stress tensor option=4 : contribution of local ionic potential to second derivative of E wrt xred
INPUTS
if(option==3) eei=local pseudopotential part of total energy (hartree) gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$). gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere). mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization mqgrid=number of grid pts in q array for f(q) spline. natom=number of atoms in unit cell. nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ntypat=number of types of atoms. option= (see above) ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. qgrid(mqgrid)=q grid for spline from 0 to qmax. qprtrb(3)= integer wavevector of possible perturbing potential in basis of reciprocal lattice translations rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$) (needed if option==2 or if option==4) ucvol=unit cell volume ($\textrm{Bohr}^3$). vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom. vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero, perturbing potential is added of the form $V(G)=(vprtrb(1)+I*vprtrb(2))/2$ at the values G=qprtrb and $(vprtrb(1)-I*vprtrb(2))/2$ at $G=-qprtrb$ (integers)
OUTPUT
(if option==1) vpsp(nfft)=local crystal pseudopotential in real space. (if option==2) grtn(3,natom)=grads of Etot wrt tn. (if option==3) lpsstr(6)=components of local psp part of stress tensor (Cartesian coordinates, symmetric tensor) in hartree/$\textrm{bohr}^3$ Store 6 unique components in order 11, 22, 33, 32, 31, 21 (if option==4) dyfrlo(3,3,natom)=d2 Eei/d tn(i)/d tn(j). (Hartrees)
NOTES
Note that the present routine is tightly connected to the dfpt_vlocal.f routine, that compute the derivative of the local ionic potential with respect to one atomic displacement. The argument list and the internal loops to be considered were sufficiently different as to make the two routine different.
SOURCE
273 subroutine mklocl_recipspace(dyfrlo,eei,gmet,gprimd,grtn,gsqcut,icutcoul,lpsstr,mgfft,& 274 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,nkpt,ntypat,option,ph1d,qgrid,qprtrb,& 275 & rcut,rhog,rprimd,ucvol,vcutgeo,vlspl,vprtrb,vpsp) 276 277 !Arguments ------------------------------------ 278 !scalars 279 integer,intent(in) :: mgfft,mqgrid,natom,nfft,nkpt,ntypat,option,icutcoul 280 real(dp),intent(in) :: eei,gsqcut,rcut,rprimd(3,3),ucvol,vcutgeo(3) 281 type(MPI_type),intent(in) :: mpi_enreg 282 !arrays 283 integer,intent(in) :: nattyp(ntypat),ngfft(18),qprtrb(3) 284 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom) 285 real(dp),intent(in) :: qgrid(mqgrid),rhog(2,nfft),vlspl(mqgrid,2,ntypat) 286 real(dp),intent(in) :: vprtrb(2) 287 real(dp),intent(out) :: dyfrlo(3,3,natom),grtn(3,natom),lpsstr(6) !vz_i 288 real(dp),intent(inout) :: vpsp(nfft) !vz_i 289 290 !Local variables------------------------------- 291 !scalars 292 integer,parameter :: im=2,re=1 293 integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ierr,ig1,ig2,ig3,ii,itypat 294 integer :: jj,me_fft,me_g0,n1,n2,n3,nproc_fft,shift1 295 integer :: shift2,shift3 296 #ifdef FC_NVHPC 297 !Silly trick to prevent NVHPC optimization issue 298 logical :: nothing=.false. 299 #endif 300 real(dp),parameter :: tolfix=1.0000001_dp 301 real(dp) :: aa,bb,cc,cutoff,dbl_ig1,dbl_ig2,dbl_ig3,dd,diff,dq,dq2div6,dqdiv6 302 real(dp) :: dqm1,ee,ff,gmag,gsquar!beta,gcart_para,gcart_perp 303 real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r 304 real(dp) :: ph3i,ph3r,phimag_igia,phre_igia,rcut_loc,sfi,sfr 305 real(dp) :: svion,svioni,svionr,term,vion1,vion2,xnorm 306 character(len=500) :: message 307 !arrays 308 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 309 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 310 real(dp) :: gcart(3),tsec(2) 311 real(dp),allocatable :: gcutoff(:) 312 real(dp),allocatable :: work1(:,:) 313 314 ! ************************************************************************* 315 316 !Define G^2 based on G space metric gmet. 317 ! gsq(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 318 !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 319 !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 320 321 !Real and imaginary parts of phase--statment functions: 322 ! phr(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 323 ! phi(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 324 ! ph1(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1)) 325 ! ph2(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+& 326 !& natom*(2*n1+1)) 327 ! ph3(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+& 328 !& natom*(2*n1+1+2*n2+1)) 329 ! phre(i1,i2,i3,ia)=phr(ph1(re,i1,ia),ph1(im,i1,ia),ph2(re,i2,ia),& 330 !& ph2(im,i2,ia),ph3(re,i3,ia),ph3(im,i3,ia)) 331 ! phimag(i1,i2,i3,ia)=phi(ph1(re,i1,ia),ph1(im,i1,ia),ph2(re,i2,ia),& 332 !& ph2(im,i2,ia),ph3(re,i3,ia),ph3(im,i3,ia)) 333 334 !----- 335 336 !Keep track of total time spent in mklocl 337 if(option==2)then 338 call timab(72,1,tsec) 339 end if 340 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 341 me_fft=ngfft(11) 342 nproc_fft=ngfft(10) 343 344 !Get the distrib associated with this fft_grid 345 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 346 347 !Zero out array to permit accumulation over atom types below: 348 if(option==1)then 349 ABI_MALLOC(work1,(2,nfft)) 350 work1(:,:)=zero 351 end if 352 353 ! 354 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 355 dqm1=1.0_dp/dq 356 dqdiv6=dq/6.0_dp 357 dq2div6=dq**2/6.0_dp 358 cutoff=gsqcut*tolfix 359 id1=n1/2+2 360 id2=n2/2+2 361 id3=n3/2+2 362 grtn(:,:)=zero 363 lpsstr(:)=zero 364 dyfrlo(:,:,:)=zero 365 me_g0=0 366 ia1=1 367 368 !Initialize Gcut-off array from m_gtermcutoff 369 call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo) 370 rcut_loc = half*SQRT(DOT_PRODUCT(rprimd(:,3),rprimd(:,3))) 371 372 do itypat=1,ntypat 373 ! ia1,ia2 sets range of loop over atoms: 374 ia2=ia1+nattyp(itypat)-1 375 376 ii=0 377 do i3=1,n3 378 ig3=i3-(i3/id3)*n3-1 379 do i2=1,n2 380 ig2=i2-(i2/id2)*n2-1 381 if(fftn2_distrib(i2) == me_fft ) then 382 do i1=1,n1 383 ig1=i1-(i1/id1)*n1-1 384 385 ii=ii+1 386 ! *** GET RID OF THIS THESE IF STATEMENTS (if they slow code) 387 ! Skip G=0: 388 ! if (ii==1) cycle 389 if (ig1==0 .and. ig2==0 .and. ig3==0) me_g0=1 390 if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 391 392 gsquar=gsq_mk(ig1,ig2,ig3) 393 ! Skip G**2 outside cutoff: 394 if (gsquar<=cutoff) then 395 gmag=sqrt(gsquar) 396 397 ! Compute vion(G) for given type of atom 398 jj=1+int(gmag*dqm1) 399 diff=gmag-qgrid(jj) 400 401 ! Evaluate spline fit from q^2 V(q) to get V(q): 402 ! (p. 86 Numerical Recipes, Press et al; 403 ! NOTE error in book for sign 404 ! of "aa" term in derivative; also see splfit routine). 405 406 bb = diff*dqm1 407 aa = 1.0_dp-bb 408 cc = aa*(aa**2-1.0_dp)*dq2div6 409 dd = bb*(bb**2-1.0_dp)*dq2div6 410 411 vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +& 412 & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) / gsquar * gcutoff(ii) 413 414 if(option==1)then 415 416 ! Assemble structure factor over all atoms of given type: 417 sfr=zero 418 sfi=zero 419 do ia=ia1,ia2 420 sfr=sfr+phre_mk(ig1,ig2,ig3,ia) 421 sfi=sfi-phimag_mk(ig1,ig2,ig3,ia) 422 end do 423 ! Multiply structure factor times vion: 424 work1(re,ii)=work1(re,ii)+sfr*vion1 425 work1(im,ii)=work1(im,ii)+sfi*vion1 426 427 else if(option==2 .or. option==4)then 428 429 ! Compute Re and Im part of (2Pi)*Vion(G)*rho(G): 430 svionr=(two_pi*vion1)*rhog(re,ii) 431 svioni=(two_pi*vion1)*rhog(im,ii) 432 433 ! Loop over all atoms of this type: 434 do ia=ia1,ia2 435 #ifdef FC_NVHPC 436 !Silly trick to prevent NVHPC optimization issue 437 if(nothing) write(100,*) shift1,shift2,shift3 438 #endif 439 440 shift1=1+n1+(ia-1)*(2*n1+1) 441 shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1) 442 shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1) 443 ph1r=ph1d(1,ig1+shift1) 444 ph1i=ph1d(2,ig1+shift1) 445 ph2r=ph1d(1,ig2+shift2) 446 ph2i=ph1d(2,ig2+shift2) 447 ph3r=ph1d(1,ig3+shift3) 448 ph3i=ph1d(2,ig3+shift3) 449 ph12r=ph1r*ph2r-ph1i*ph2i 450 ph12i=ph1r*ph2i+ph1i*ph2r 451 phre_igia=ph12r*ph3r-ph12i*ph3i 452 phimag_igia=ph12r*ph3i+ph12i*ph3r 453 454 if(option==2)then 455 456 ! Compute "Vion" part of gradient 457 ! svion=svioni*phre(ig1,ig2,ig3,ia)+svionr*phimag(ig1,ig2,ig3,ia) 458 svion=svioni*phre_igia+svionr*phimag_igia 459 460 ! Open loop over 3-index for speed: 461 grtn(1,ia)=grtn(1,ia)-dble(ig1)*svion 462 grtn(2,ia)=grtn(2,ia)-dble(ig2)*svion 463 grtn(3,ia)=grtn(3,ia)-dble(ig3)*svion 464 465 else 466 467 ! Compute "Vion" part of the second derivative 468 ! svion=two_pi* 469 ! (svionr*phre(ig1,ig2,ig3,ia)-svioni*phimag(ig1,ig2,ig3,ia)) 470 svion=two_pi*(svionr*phre_igia-svioni*phimag_igia) 471 472 ! Open loop over 3-index for speed 473 dbl_ig1=dble(ig1) ; dbl_ig2=dble(ig2) ; dbl_ig3=dble(ig3) 474 dyfrlo(1,1,ia)=dyfrlo(1,1,ia)-dbl_ig1*dbl_ig1*svion 475 dyfrlo(1,2,ia)=dyfrlo(1,2,ia)-dbl_ig1*dbl_ig2*svion 476 dyfrlo(1,3,ia)=dyfrlo(1,3,ia)-dbl_ig1*dbl_ig3*svion 477 dyfrlo(2,2,ia)=dyfrlo(2,2,ia)-dbl_ig2*dbl_ig2*svion 478 dyfrlo(2,3,ia)=dyfrlo(2,3,ia)-dbl_ig2*dbl_ig3*svion 479 dyfrlo(3,3,ia)=dyfrlo(3,3,ia)-dbl_ig3*dbl_ig3*svion 480 481 end if 482 483 end do 484 485 else if(option==3)then 486 ! if(icutcoul .ne. 2) then 487 ! Also get (dV(q)/dq)/q: 488 ! (note correction of Numerical Recipes sign error 489 ! before (3._dp*aa**2-1._dp) 490 ! ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines 491 ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat) 492 ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) & 493 & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat) 494 vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag& 495 & - 2.0_dp*vion1 ) / gsquar 496 497 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+& 498 & gprimd(1,3)*dble(ig3) 499 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+& 500 & gprimd(2,3)*dble(ig3) 501 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+& 502 & gprimd(3,3)*dble(ig3) 503 ! Assemble structure over all atoms of given type 504 sfr=zero 505 sfi=zero 506 do ia=ia1,ia2 507 sfr=sfr+phre_mk(ig1,ig2,ig3,ia) 508 sfi=sfi-phimag_mk(ig1,ig2,ig3,ia) 509 end do 510 ! Compute Re( rho^*(G)* sf ) * [(dV(G)/dG)/|G|] 511 term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)*vion2 512 ! Compute contribution to stress tensor 513 lpsstr(1)=lpsstr(1)-term*gcart(1)*gcart(1) 514 lpsstr(2)=lpsstr(2)-term*gcart(2)*gcart(2) 515 lpsstr(3)=lpsstr(3)-term*gcart(3)*gcart(3) 516 lpsstr(4)=lpsstr(4)-term*gcart(3)*gcart(2) 517 lpsstr(5)=lpsstr(5)-term*gcart(3)*gcart(1) 518 lpsstr(6)=lpsstr(6)-term*gcart(2)*gcart(1) 519 ! else if (icutcoul .eq. 2) then 520 !! Also get (dV(q)/dq)/q: 521 !! (note correction of Numerical Recipes sign error 522 !! before (3._dp*aa**2-1._dp) 523 !! ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines 524 ! ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat) 525 ! ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) & 526 !& - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat) 527 ! vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag& 528 !& - 2.0_dp*vion1 ) / gsquar 529 ! 530 ! gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+& 531 !& gprimd(1,3)*dble(ig3) 532 ! gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+& 533 !& gprimd(2,3)*dble(ig3) 534 ! gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+& 535 !& gprimd(3,3)*dble(ig3) 536 !! Assemble structure over all atoms of given type 537 ! sfr=zero 538 ! sfi=zero 539 ! do ia=ia1,ia2 540 ! sfr=sfr+phre_mk(ig1,ig2,ig3,ia) 541 ! sfi=sfi-phimag_mk(ig1,ig2,ig3,ia) 542 ! end do 543 ! !Implement beta correction as in eq. 62 (PRB 96 075448 2017) 544 ! gcart_para = sqrt(gcart(1)**2+gcart(2)**2) 545 ! gcart_perp = gcart(3) 546 ! gsquar = gcart(1)**2+gcart(2)**2+gcart(3)**2 547 ! if(gcart_para .gt. tol12) then 548 ! beta = gsquar*rcut_loc/(two*gcart_para)* & 549 ! & exp(-gcart_para*rcut_loc)* & 550 ! &cos(gcart_perp*rcut_loc)/(one-exp(-gcart_para*rcut_loc)*cos(gcart_perp*rcut_loc)) 551 ! else 552 ! beta = zero 553 ! end if 554 !! Compute Re( rho^*(G)* sf ) * [(dV(G)/dG)/|G|] 555 ! term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)*vion2 556 !! Compute contribution to stress tensor 557 ! lpsstr(1)=lpsstr(1)-term*(gcart(1)*gcart(1))*(1+beta) 558 ! lpsstr(2)=lpsstr(2)-term*(gcart(2)*gcart(2))*(1+beta) 559 ! lpsstr(3)=lpsstr(3)-term*(gcart(3)*gcart(3)-gsquar) 560 ! lpsstr(4)=lpsstr(4)-term*gcart(3)*gcart(2) 561 ! lpsstr(5)=lpsstr(5)-term*gcart(3)*gcart(1) 562 ! lpsstr(6)=lpsstr(6)-term*gcart(2)*gcart(1) 563 ! endif 564 565 else 566 write(message, '(a,i0,a)' )' mklocl: Option=',option,' not allowed.' 567 ABI_BUG(message) 568 end if ! End option choice 569 570 ! End skip G**2 outside cutoff: 571 end if 572 573 ! End loop on n1, n2, n3. There is a "cycle" inside the loop 574 end do 575 end if ! this plane is for me_fft 576 end do 577 end do 578 579 ! Symmetrize the dynamical matrix with respect to indices 580 do ia=ia1,ia2 581 dyfrlo(2,1,ia)=dyfrlo(1,2,ia) 582 dyfrlo(3,1,ia)=dyfrlo(1,3,ia) 583 dyfrlo(3,2,ia)=dyfrlo(2,3,ia) 584 end do 585 586 ia1=ia2+1 587 588 ! End loop on type of atoms 589 end do 590 591 if(option==1)then 592 ! Dont't change work1 on g=0 if Poisson solver is used since work1 593 ! hold not the potential but the density generated by the pseudo. 594 if(me_g0 == 1) then 595 ! Set Vloc(G=0)=0: 596 work1(re,1)=zero 597 work1(im,1)=zero 598 end if 599 ! write(std_out,*) ' mklocl_recipspace : will add potential with strength vprtrb(:)=',vprtrb(:) 600 601 ! Allow for the addition of a perturbing potential 602 if ((vprtrb(1)**2+vprtrb(2)**2) > 1.d-30) then 603 ! Find the linear indices which correspond with the input 604 ! wavevector qprtrb 605 ! The double modulus handles both i>=n and i<0, mapping into [0,n-1]; 606 ! then add 1 to get range [1,n] for each 607 i3=1+mod(n3+mod(qprtrb(3),n3),n3) 608 i2=1+mod(n2+mod(qprtrb(2),n2),n2) 609 i1=1+mod(n1+mod(qprtrb(1),n1),n1) 610 ! Compute the linear index in the 3 dimensional array 611 ii=i1+n1*((ffti2_local(i2)-1)+(n2/nproc_fft)*(i3-1)) 612 ! Add in the perturbation at G=qprtrb 613 work1(re,ii)=work1(re,ii)+0.5_dp*vprtrb(1) 614 work1(im,ii)=work1(im,ii)+0.5_dp*vprtrb(2) 615 ! Same thing for G=-qprtrb 616 i3=1+mod(n3+mod(-qprtrb(3),n3),n3) 617 i2=1+mod(n2+mod(-qprtrb(2),n2),n2) 618 i1=1+mod(n1+mod(-qprtrb(1),n1),n1) 619 ! ii=i1+n1*((i2-1)+n2*(i3-1)) 620 work1(re,ii)=work1(re,ii)+0.5_dp*vprtrb(1) 621 work1(im,ii)=work1(im,ii)-0.5_dp*vprtrb(2) 622 write(message, '(a,1p,2e12.4,a,0p,3i4,a)' )& 623 & ' mklocl: perturbation of vprtrb=', vprtrb,& 624 & ' and q=',qprtrb,' has been added' 625 call wrtout(std_out,message,'COLL') 626 end if 627 628 ! Transform back to real space 629 call fourdp(1,work1,vpsp,1,mpi_enreg,nfft,1,ngfft,0) 630 631 ! Divide by unit cell volume 632 xnorm=1.0_dp/ucvol 633 vpsp(:)=vpsp(:)*xnorm 634 635 ABI_FREE(work1) 636 637 end if 638 639 ABI_FREE(gcutoff) 640 641 if(option==2)then 642 ! Init mpi_comm 643 if(mpi_enreg%nproc_fft>1)then 644 call timab(48,1,tsec) 645 call xmpi_sum(grtn,mpi_enreg%comm_fft ,ierr) 646 call timab(48,2,tsec) 647 end if 648 call timab(72,2,tsec) 649 end if 650 651 if(option==3)then 652 ! Init mpi_comm 653 if(mpi_enreg%nproc_fft>1)then 654 call timab(48,1,tsec) 655 call xmpi_sum(lpsstr,mpi_enreg%comm_fft ,ierr) 656 call timab(48,2,tsec) 657 end if 658 659 ! Normalize and add term -eei/ucvol on diagonal 660 ! (see page 802 of notes) 661 ! if(icutcoul .ne. 2) then 662 lpsstr(1)=(lpsstr(1)-eei)/ucvol 663 lpsstr(2)=(lpsstr(2)-eei)/ucvol 664 lpsstr(3)=(lpsstr(3)-eei)/ucvol 665 lpsstr(4)=lpsstr(4)/ucvol 666 lpsstr(5)=lpsstr(5)/ucvol 667 lpsstr(6)=lpsstr(6)/ucvol 668 ! elseif (icutcoul .eq. 2) then 669 ! lpsstr(1)=(lpsstr(1)-eei)/ucvol 670 ! lpsstr(2)=(lpsstr(2)-eei)/ucvol 671 ! lpsstr(3)=(lpsstr(3)-eei)/ucvol 672 ! lpsstr(4)=lpsstr(4)/ucvol 673 ! lpsstr(5)=lpsstr(5)/ucvol 674 ! lpsstr(6)=lpsstr(6)/ucvol 675 !lpsstr=lpsstr/ucvol 676 ! endif 677 678 end if 679 680 if(option==4)then 681 ! Init mpi_comm 682 if(mpi_enreg%nproc_fft>1)then 683 call timab(48,1,tsec) 684 call xmpi_sum(dyfrlo,mpi_enreg%comm_fft ,ierr) 685 call timab(48,2,tsec) 686 end if 687 end if 688 689 contains 690 691 !Real and imaginary parts of phase--statment functions: 692 function phr_mk(x1,y1,x2,y2,x3,y3) 693 694 real(dp) :: phr_mk,x1,x2,x3,y1,y2,y3 695 phr_mk=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 696 end function phr_mk 697 698 function phi_mk(x1,y1,x2,y2,x3,y3) 699 700 real(dp):: phi_mk,x1,x2,x3,y1,y2,y3 701 phi_mk=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 702 end function phi_mk 703 704 function ph1_mk(nri,ig1,ia) 705 706 real(dp):: ph1_mk 707 integer :: nri,ig1,ia 708 ph1_mk=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1)) 709 end function ph1_mk 710 711 function ph2_mk(nri,ig2,ia) 712 713 real(dp):: ph2_mk 714 integer :: nri,ig2,ia 715 ph2_mk=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)) 716 end function ph2_mk 717 718 function ph3_mk(nri,ig3,ia) 719 720 real(dp):: ph3_mk 721 integer :: nri,ig3,ia 722 ph3_mk=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 723 end function ph3_mk 724 725 function phre_mk(ig1,ig2,ig3,ia) 726 727 real(dp):: phre_mk 728 integer :: ig1,ig2,ig3,ia 729 phre_mk=phr_mk(ph1_mk(re,ig1,ia),ph1_mk(im,ig1,ia),& 730 & ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia)) 731 end function phre_mk 732 733 function phimag_mk(ig1,ig2,ig3,ia) 734 735 real(dp) :: phimag_mk 736 integer :: ig1,ig2,ig3,ia 737 phimag_mk=phi_mk(ph1_mk(re,ig1,ia),ph1_mk(im,ig1,ia),& 738 & ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia)) 739 end function phimag_mk 740 741 function gsq_mk(i1,i2,i3) 742 743 real(dp) :: gsq_mk 744 integer :: i1,i2,i3 745 gsq_mk=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 746 & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 747 & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 748 end function gsq_mk 749 750 end subroutine mklocl_recipspace
ABINIT/vlocalstr [ Functions ]
NAME
vlocalstr
FUNCTION
Compute strain derivatives of local ionic potential
second derivative of E wrt xred
INPUTS
gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$). gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere). istr=1,...6 specifies cartesian strain component 11,22,33,32,31,21 mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization mqgrid=number of grid pts in q array for f(q) spline. natom=number of atoms in unit cell. nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ntypat=number of types of atoms. ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. qgrid(mqgrid)=q grid for spline from 0 to qmax. ucvol=unit cell volume ($\textrm{Bohr}^3$). vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom. [g0term]= optional, if present an alternative treatment of the G=0 term, adopoted for the flexoelectric tensor calculation, is performed.
OUTPUT
vpsp1(nfft)=first-order local crystal pseudopotential in real space.
NOTES
* Note that the present routine is tightly connected to the dfpt_vlocal.f routine, that compute the derivative of the local ionic potential with respect to one atomic displacement. The argument list and the internal loops to be considered were sufficiently different as to make the two routines different. * The routine was adapted from mklocl.F90
SOURCE
1051 subroutine vlocalstr(gmet,gprimd,gsqcut,istr,mgfft,mpi_enreg,& 1052 & mqgrid,natom,nattyp,nfft,ngfft,ntypat,ph1d,qgrid,& 1053 & ucvol,vlspl,vpsp1,g0term) 1054 1055 !Arguments ------------------------------------ 1056 !scalars 1057 integer,intent(in) :: istr,mgfft,mqgrid,natom,nfft,ntypat 1058 integer,optional,intent(in) :: g0term 1059 real(dp),intent(in) :: gsqcut,ucvol 1060 type(MPI_type),intent(in) :: mpi_enreg 1061 !arrays 1062 integer,intent(in) :: nattyp(ntypat),ngfft(18) 1063 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom) 1064 real(dp),intent(in) :: qgrid(mqgrid),vlspl(mqgrid,2,ntypat) 1065 real(dp),intent(out) :: vpsp1(nfft) 1066 1067 !Local variables------------------------------- 1068 !scalars 1069 integer,parameter :: im=2,re=1 1070 integer :: g0term_ 1071 integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ig1,ig2,ig3,ii,itypat,jj 1072 integer :: ka,kb,n1,n2,n3 1073 real(dp),parameter :: tolfix=1.0000001_dp 1074 real(dp) :: aa,bb,cc,cutoff,dd,dgsquards,diff 1075 real(dp) :: dq,dq2div6,dqdiv6,dqm1,ee,ff,gmag,gsquar 1076 real(dp) :: sfi,sfr,term,vion1,vion2,vlocg0 1077 real(dp) :: xnorm 1078 character(len=500) :: message 1079 !arrays 1080 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 1081 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1082 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1083 real(dp) :: dgmetds(3,3) 1084 real(dp),allocatable :: work1(:,:) 1085 1086 ! ************************************************************************* 1087 1088 !Define G^2 based on G space metric gmet. 1089 ! gsq_vl(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 1090 !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 1091 !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 1092 1093 !Define dG^2/ds based on G space metric derivative dgmetds. 1094 ! dgsqds_vl(i1,i2,i3)=dble(i1*i1)*dgmetds(1,1)+dble(i2*i2)*dgmetds(2,2)+& 1095 !& dble(i3*i3)*dgmetds(3,3)+& 1096 !& dble(i1*i2)*(dgmetds(1,2)+dgmetds(2,1))+& 1097 !& dble(i1*i3)*(dgmetds(1,3)+dgmetds(3,1))+& 1098 !& dble(i2*i3)*(dgmetds(2,3)+dgmetds(3,2)) 1099 1100 !Real and imaginary parts of phase--statment functions: 1101 ! phr_vl(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 1102 ! phi_vl(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 1103 ! ph1_vl(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1)) 1104 ! ph2_vl(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+& 1105 !& natom*(2*n1+1)) 1106 ! ph3_vl(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+& 1107 !& natom*(2*n1+1+2*n2+1)) 1108 ! phre_vl(i1,i2,i3,ia)=phr_vl(ph1_vl(re,i1,ia),ph1_vl(im,i1,ia),ph2_vl(re,i2,ia),& 1109 !& ph2_vl(im,i2,ia),ph3_vl(re,i3,ia),ph3_vl(im,i3,ia)) 1110 ! phimag_vl(i1,i2,i3,ia)=phi_vl(ph1_vl(re,i1,ia),ph1_vl(im,i1,ia),ph2_vl(re,i2,ia),& 1111 !& ph2_vl(im,i2,ia),ph3_vl(re,i3,ia),ph3_vl(im,i3,ia)) 1112 1113 !----- 1114 !Compute derivative of metric tensor wrt strain component istr 1115 if(istr<1 .or. istr>6)then 1116 write(message, '(a,i10,a,a,a)' )& 1117 & ' Input istr=',istr,' not allowed.',ch10,& 1118 & ' Possible values are 1,2,3,4,5,6 only.' 1119 ABI_BUG(message) 1120 end if 1121 1122 ka=idx(2*istr-1);kb=idx(2*istr) 1123 do ii = 1,3 1124 dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii)) 1125 end do 1126 !For historical reasons: 1127 dgmetds(:,:)=0.5_dp*dgmetds(:,:) 1128 1129 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1130 1131 !Get the distrib associated with this fft_grid 1132 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1133 1134 !Zero out array to permit accumulation over atom types below: 1135 ABI_MALLOC(work1,(2,nfft)) 1136 work1(:,:)=0.0_dp 1137 ! 1138 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 1139 dqm1=1.0_dp/dq 1140 dqdiv6=dq/6.0_dp 1141 dq2div6=dq**2/6.0_dp 1142 cutoff=gsqcut*tolfix 1143 id1=n1/2+2 1144 id2=n2/2+2 1145 id3=n3/2+2 1146 1147 ia1=1 1148 do itypat=1,ntypat 1149 ! ia1,ia2 sets range of loop over atoms: 1150 ia2=ia1+nattyp(itypat)-1 1151 1152 ii=0 1153 do i3=1,n3 1154 ig3=i3-(i3/id3)*n3-1 1155 do i2=1,n2 1156 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 1157 ig2=i2-(i2/id2)*n2-1 1158 do i1=1,n1 1159 ig1=i1-(i1/id1)*n1-1 1160 ii=ii+1 1161 ! *** GET RID OF THIS THESE IF STATEMENTS (if they slow code) 1162 ! Skip G=0: 1163 ! if (ii==1) cycle 1164 if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 1165 gsquar=gsq_vl(ig1,ig2,ig3) 1166 1167 ! Skip G**2 outside cutoff: 1168 if (gsquar<=cutoff) then 1169 gmag=sqrt(gsquar) 1170 dgsquards=dgsqds_vl(ig1,ig2,ig3) 1171 ! Compute vion(G) for given type of atom 1172 jj=1+int(gmag*dqm1) 1173 diff=gmag-qgrid(jj) 1174 1175 ! Evaluate spline fit from q^2 V(q) to get V(q): 1176 ! (p. 86 Numerical Recipes, Press et al; 1177 ! NOTE error in book for sign 1178 ! of "aa" term in derivative; also see splfit routine). 1179 1180 bb = diff*dqm1 1181 aa = 1.0_dp-bb 1182 cc = aa*(aa**2-1.0_dp)*dq2div6 1183 dd = bb*(bb**2-1.0_dp)*dq2div6 1184 1185 vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +& 1186 & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) & 1187 & / gsquar 1188 1189 ! Also get (dV(q)/dq)/q: 1190 ! (note correction of Numerical Recipes sign error 1191 ! before (3._dp*aa**2-1._dp) 1192 ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat) 1193 ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) & 1194 & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat) 1195 vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag& 1196 & - 2.0_dp*vion1 ) / gsquar 1197 1198 1199 ! Assemble structure factor over all atoms of given type: 1200 sfr=0.0_dp 1201 sfi=0.0_dp 1202 do ia=ia1,ia2 1203 sfr=sfr+phre_vl(ig1,ig2,ig3,ia) 1204 sfi=sfi-phimag_vl(ig1,ig2,ig3,ia) 1205 end do 1206 1207 term=dgsquards*vion2 1208 ! Add potential for diagonal strain components 1209 if(istr <=3) then 1210 term=term-vion1 1211 end if 1212 1213 ! Multiply structure factor times vion derivatives: 1214 work1(re,ii)=work1(re,ii)+sfr*term 1215 work1(im,ii)=work1(im,ii)+sfi*term 1216 1217 ! End skip G**2 outside cutoff: 1218 end if 1219 ! End loop on n1, n2, n3. There is a "cycle" inside the loop 1220 end do 1221 end if 1222 end do 1223 end do 1224 1225 ia1=ia2+1 1226 1227 ! End loop on type of atoms 1228 end do 1229 1230 1231 !Set Vloc(G=0)=0: 1232 work1(re,1)=0.0_dp 1233 work1(im,1)=0.0_dp 1234 1235 !Alternative treatment of Vloc(G=0) for the flexoelectric tensor calculation 1236 g0term_=0; if (present(g0term)) g0term_=g0term 1237 if (g0term_==1) then 1238 vlocg0=zero 1239 if (istr<=3) then 1240 ia1=1 1241 do itypat=1,ntypat 1242 ! ia1,ia2 sets range of loop over atoms: 1243 1244 ia2=ia1+nattyp(itypat)-1 1245 do ia=ia1,ia2 1246 vlocg0=vlocg0+vlspl(1,2,itypat) 1247 end do 1248 end do 1249 work1(re,1)=-half*vlocg0 1250 end if 1251 end if 1252 1253 !Transform back to real space 1254 call fourdp(1,work1,vpsp1,1,mpi_enreg,nfft,1,ngfft,0) 1255 1256 !Divide by unit cell volume 1257 xnorm=1.0_dp/ucvol 1258 vpsp1(:)=vpsp1(:)*xnorm 1259 1260 ABI_FREE(work1) 1261 1262 contains 1263 1264 !Real and imaginary parts of phase. 1265 function phr_vl(x1,y1,x2,y2,x3,y3) 1266 1267 real(dp) :: phr_vl,x1,x2,x3,y1,y2,y3 1268 phr_vl=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 1269 end function phr_vl 1270 1271 function phi_vl(x1,y1,x2,y2,x3,y3) 1272 1273 real(dp):: phi_vl,x1,x2,x3,y1,y2,y3 1274 phi_vl=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 1275 end function phi_vl 1276 1277 function ph1_vl(nri,ig1,ia) 1278 1279 real(dp):: ph1_vl 1280 integer :: nri,ig1,ia 1281 ph1_vl=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1)) 1282 end function ph1_vl 1283 1284 function ph2_vl(nri,ig2,ia) 1285 1286 real(dp):: ph2_vl 1287 integer :: nri,ig2,ia 1288 ph2_vl=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)) 1289 end function ph2_vl 1290 1291 function ph3_vl(nri,ig3,ia) 1292 1293 real(dp):: ph3_vl 1294 integer :: nri,ig3,ia 1295 ph3_vl=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 1296 end function ph3_vl 1297 1298 function phre_vl(ig1,ig2,ig3,ia) 1299 1300 real(dp):: phre_vl 1301 integer :: ig1,ig2,ig3,ia 1302 phre_vl=phr_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),& 1303 & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia)) 1304 end function phre_vl 1305 1306 function phimag_vl(ig1,ig2,ig3,ia) 1307 1308 real(dp) :: phimag_vl 1309 integer :: ig1,ig2,ig3,ia 1310 phimag_vl=phi_vl(ph1_vl(re,ig1,ia),ph1_vl(im,ig1,ia),& 1311 & ph2_vl(re,ig2,ia),ph2_vl(im,ig2,ia),ph3_vl(re,ig3,ia),ph3_vl(im,ig3,ia)) 1312 end function phimag_vl 1313 1314 function gsq_vl(i1,i2,i3) 1315 1316 real(dp) :: gsq_vl 1317 integer :: i1,i2,i3 1318 !Define G^2 based on G space metric gmet. 1319 gsq_vl=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 1320 & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 1321 & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 1322 end function gsq_vl 1323 1324 function dgsqds_vl(i1,i2,i3) 1325 1326 real(dp) :: dgsqds_vl 1327 integer :: i1,i2,i3 1328 !Define dG^2/ds based on G space metric derivative dgmetds. 1329 dgsqds_vl=dble(i1*i1)*dgmetds(1,1)+dble(i2*i2)*dgmetds(2,2)+& 1330 & dble(i3*i3)*dgmetds(3,3)+& 1331 & dble(i1*i2)*(dgmetds(1,2)+dgmetds(2,1))+& 1332 & dble(i1*i3)*(dgmetds(1,3)+dgmetds(3,1))+& 1333 & dble(i2*i3)*(dgmetds(2,3)+dgmetds(3,2)) 1334 end function dgsqds_vl 1335 1336 end subroutine vlocalstr