TABLE OF CONTENTS
- ABINIT/m_ifc
- m_ifc/corsifc9
- m_ifc/ifc_autocutoff
- m_ifc/ifc_calcnwrite_nana_terms
- m_ifc/ifc_fourq
- m_ifc/ifc_free
- m_ifc/ifc_from_file
- m_ifc/ifc_get_dwdq
- m_ifc/ifc_get_phmesh
- m_ifc/ifc_getiaf
- m_ifc/ifc_init
- m_ifc/ifc_outphbtrap
- m_ifc/ifc_print
- m_ifc/ifc_printbxsf
- m_ifc/ifc_speedofsound
- m_ifc/ifc_type
- m_ifc/ifc_write
- m_ifc/omega_decomp
ABINIT/m_ifc [ Modules ]
NAME
m_ifc
FUNCTION
This module contains the declaration of data types and methods used to handle interatomic force constant sets
COPYRIGHT
Copyright (C) 2011-2024 ABINIT group (XG,MJV,EB,MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_ifc 24 25 use defs_basis 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_sort 30 use m_cgtools 31 use m_nctk 32 use m_ddb 33 use m_ddb_hdr 34 use m_symkpt 35 #ifdef HAVE_NETCDF 36 use netcdf 37 #endif 38 39 use m_io_tools, only : open_file 40 use m_numeric_tools, only : arth 41 use m_fstrings, only : ktoa, int2char4, sjoin, itoa, ltoa, ftoa 42 use m_symtk, only : matr3inv 43 use m_special_funcs, only : abi_derfc 44 use m_time, only : cwtime, cwtime_report, timab 45 use m_copy, only : alloc_copy 46 use m_pptools, only : printbxsf 47 use m_lebedev, only : lebedev_t, lebedev_ngrids 48 use m_ewald, only : ewald9 49 use m_crystal, only : crystal_t 50 use m_geometry, only : phdispl_cart2red, normv, mkrdim 51 use m_kpts, only : kpts_ibz_from_kptrlatt, smpbz 52 use m_dynmat, only : canct9, dist9 , ifclo9, axial9, q0dy3_apply, q0dy3_calc, asrif9, dynmat_dq, & 53 make_bigbox, canat9, chkrp9, ftifc_q2r, wght9, nanal9, gtdyn9, dymfz9, & 54 massmult_and_breaksym, dfpt_phfrq, dfpt_prtph 55 56 implicit none 57 58 private
m_ifc/corsifc9 [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
corsifc9
FUNCTION
Applies a cutoff on the ifc in real space
INPUTS
acell(3)=length scales by which rprim is to be multiplied gprim(3,3)=dimensionless primitive translations in reciprocal space natom=number of atoms in unit cell nrpt= Number of R points in the Big Box rcan(3,natom)=canonical coordinates of atoms rprim(3,3)=dimensionless primitive translations in real space rpt(3,nrpt)=canonical coordinates of the points in the BigBox. nsphere=number of atoms to be included in the cut-off sphere for interatomic force constant; if = 0 : maximum extent allowed by the grid. rifcsph=radius for cutoff of IFC wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector.
OUTPUT
wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector with the required cutoff applied. rcut_min=Effective cutoff. Defined by the minimum cutoff radius over the natom sites.
SOURCE
1548 subroutine corsifc9(acell,gprim,natom,nrpt,nsphere,rifcsph,rcan,rprim,rpt,rcut_min,wghatm) 1549 1550 !Arguments ------------------------------- 1551 !scalars 1552 integer,intent(in) :: natom,nrpt,nsphere 1553 real(dp),intent(in) :: rifcsph 1554 real(dp),intent(out) :: rcut_min 1555 !arrays 1556 real(dp),intent(in) :: acell(3) 1557 real(dp),intent(in) :: gprim(3,3),rcan(3,natom) 1558 real(dp),intent(in) :: rprim(3,3),rpt(3,nrpt) 1559 real(dp),intent(inout) :: wghatm(natom,natom,nrpt) 1560 1561 !Local variables ------------------------- 1562 !scalars 1563 integer :: ia,ib,ii,index,irpt 1564 real(dp) :: rmax,rsigma,r0 1565 !arrays 1566 integer,allocatable :: list(:) 1567 real(dp),allocatable :: dist(:,:,:),wkdist(:) 1568 1569 ! ********************************************************************* 1570 1571 ! Compute the distances between atoms 1572 ! dist(ia,ib,irpt) contains the distance from atom ia to atom ib in unit cell irpt. 1573 ABI_MALLOC(dist,(natom,natom,nrpt)) 1574 call dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt) 1575 1576 ABI_MALLOC(list,(natom*nrpt)) 1577 ABI_MALLOC(wkdist,(natom*nrpt)) 1578 1579 ! loop on all generic atoms. 1580 rcut_min = huge(one) 1581 do ia=1,natom 1582 1583 wkdist = reshape(dist(ia,:,:), [natom*nrpt]) 1584 do ii=1,natom*nrpt 1585 list(ii)=ii 1586 end do 1587 ! This sorting algorithm is slow ... 1588 call sort_dp(natom*nrpt,wkdist,list,tol14) 1589 rmax = wkdist(natom*nrpt) 1590 1591 ! zero the outside IFCs: act on wghatm 1592 1593 ! fix number of spheres 1594 if(nsphere/=0.and.nsphere<natom*nrpt)then 1595 rcut_min = min(rcut_min, wkdist(nsphere+1)) 1596 do ii=nsphere+1,natom*nrpt 1597 index=list(ii) 1598 irpt=(index-1)/natom+1 1599 ib=index-natom*(irpt-1) 1600 wghatm(ia,ib,irpt)=zero 1601 end do 1602 end if 1603 1604 ! or fix radius of maximum ifc 1605 if(rifcsph>tol10)then 1606 do ii=nsphere+1,natom*nrpt 1607 index=list(ii) 1608 ! preserve weights for atoms inside sphere of radius rifcsph 1609 if (wkdist(ii) < rifcsph) cycle 1610 rcut_min = min(rcut_min, wkdist(ii)) 1611 irpt=(index-1)/natom+1 1612 ib=index-natom*(irpt-1) 1613 wghatm(ia,ib,irpt)=zero 1614 end do 1615 end if 1616 1617 ! filter smoothly to 0 at edge of WS cell 1618 if (rifcsph < -tol10) then 1619 ! Use different filter 1620 r0 = abs(rifcsph) * rmax; rsigma = half*(rmax-r0) !one 1621 rcut_min = r0 ! Set it to r0 1622 do ii=nsphere+1,natom*nrpt 1623 index=list(ii) 1624 irpt=(index-1)/natom+1 1625 ib=index-natom*(irpt-1) 1626 wghatm(ia,ib,irpt) = wghatm(ia,ib,irpt) * half * abi_derfc((wkdist(ii) - r0) / rsigma) 1627 end do 1628 end if 1629 1630 end do 1631 1632 ABI_FREE(dist) 1633 ABI_FREE(list) 1634 ABI_FREE(wkdist) 1635 1636 end subroutine corsifc9
m_ifc/ifc_autocutoff [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_autocutoff
FUNCTION
Find the value of nsphere that gives non-negative frequencies around Gamma in a small sphere of radius qrad. Use bisection to reduce the number of attemps although there's no guarantee that the number of negative frequencies is monotonic.
INPUTS
crystal<crystal_t> = Information on the crystalline structure. comm=MPI communicator
SIDE EFFECTS
ifc%wghatm(natom,natom,nrpt) = Weights associated to a pair of atoms and to a R vector
with the last cutoff found by the bisection algorithm applied.
ifc%atmfrc(2,3,natom,3,natom,nrpt)= ASR-imposed Interatomic Forces
SOURCE
1401 subroutine ifc_autocutoff(ifc, crystal, comm) 1402 1403 !Arguments ------------------------------------ 1404 !scalars 1405 class(ifc_type),intent(inout) :: ifc 1406 type(crystal_t),intent(in) :: crystal 1407 integer,intent(in) :: comm 1408 1409 !Local variables------------------------------- 1410 !scalars 1411 integer,parameter :: master=0 1412 integer :: iq_ibz,ierr,my_rank,nprocs,ii,nsphere,num_negw,jl,ju,jm,natom,nrpt 1413 real(dp),parameter :: rifcsph0=zero 1414 real(dp) :: adiff,qrad,min_negw,xval,rcut_min 1415 type(lebedev_t) :: lgrid 1416 !arrays 1417 real(dp) :: displ_cart(2*3*ifc%natom*3*ifc%natom) 1418 real(dp) :: qred(3),qred_vers(3),phfrqs(3*ifc%natom) !,dwdq(3,3*ifc%natom) 1419 real(dp),allocatable :: ref_phfrq(:,:),cut_phfrq(:,:) 1420 real(dp),allocatable :: save_wghatm(:,:,:),save_atmfrc(:,:,:,:,:) 1421 1422 ! ********************************************************************* 1423 1424 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1425 natom = ifc%natom; nrpt = ifc%nrpt 1426 1427 ! Compute frequencies on the ab-initio q-mesh without cutoff. 1428 ABI_CALLOC(ref_phfrq, (3*natom, ifc%nqibz)) 1429 do iq_ibz=1,ifc%nqibz 1430 if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism 1431 call ifc%fourq(crystal, ifc%qibz(:,iq_ibz), ref_phfrq(:,iq_ibz), displ_cart) 1432 end do 1433 call xmpi_sum(ref_phfrq, comm, ierr) 1434 1435 ABI_MALLOC(save_wghatm, (natom,natom,nrpt)) 1436 ABI_MALLOC(save_atmfrc, (3,natom,3,natom,ifc%nrpt)) 1437 save_wghatm = ifc%wghatm; save_atmfrc = ifc%atmfrc 1438 1439 ABI_MALLOC(cut_phfrq, (3*natom, ifc%nqibz)) 1440 qrad = 0.01 1441 call lgrid%from_index(16) 1442 1443 if (my_rank == master) then 1444 write(ab_out, "(a)")" Apply cutoff on IFCs. Using bisection algorithm to find initial guess for nsphere." 1445 write(ab_out, "(a,i0)")" Maximum nuber of atom-centered spheres: ",natom * nrpt 1446 write(ab_out, "(a,i0,a,f5.3)")" Using Lebedev-Laikov grid with npts: ",lgrid%npts, ", qrad: ",qrad 1447 write(ab_out, "(/,a)")" <adiff>: Average difference between ab-initio frequencies and frequencies with cutoff." 1448 write(ab_out, "(a)")" num_negw: Number of negative freqs detected in small sphere around Gamma." 1449 write(ab_out, "(a)")" min_negw: Min negative frequency on the small sphere." 1450 write(ab_out, "(a,/,/)")" rifcsph: Effective cutoff radius corresponding to nsphere." 1451 write(ab_out, "(a)")" nsphere <adiff>[meV] num_negw min_negw[meV] rifcsph" 1452 end if 1453 1454 jl = 0; ju = natom * nrpt + 1 ! Initialize lower and upper limits. 1455 do 1456 if (ju - jl <= 1) then 1457 exit 1458 end if 1459 jm = (ju + jl) / 2 ! Compute a midpoint 1460 nsphere = jm 1461 1462 ifc%wghatm = save_wghatm; ifc%atmfrc = save_atmfrc 1463 call corsifc9(ifc%acell,ifc%gprim, natom, nrpt,nsphere,rifcsph0,ifc%rcan,ifc%rprim,ifc%rpt,rcut_min,ifc%wghatm) 1464 if (ifc%asr > 0) call asrif9(ifc%asr,ifc%atmfrc,ifc%natom,ifc%nrpt,ifc%rpt,ifc%wghatm) 1465 1466 cut_phfrq = zero 1467 do iq_ibz=1,ifc%nqibz 1468 if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism 1469 call ifc%fourq(crystal, ifc%qibz(:,iq_ibz), cut_phfrq(:,iq_ibz), displ_cart) 1470 !write(std_out,*)cut_phfrq(1,iq_ibz),ref_phfrq(1,iq_ibz) 1471 end do 1472 call xmpi_sum(cut_phfrq, comm, ierr) 1473 1474 ! Test wether there are negative frequencies around gamma, including reciprocal lattice vectors. 1475 num_negw = 0; min_negw = zero 1476 do ii=1,lgrid%npts+3 1477 if (mod(ii, nprocs) /= my_rank) cycle ! mpi-parallelism 1478 if (ii <= 3) then 1479 qred = zero; qred(ii) = one 1480 else 1481 qred = matmul(crystal%rprimd, lgrid%versors(:, ii-3)) 1482 end if 1483 qred_vers = (qred / normv(qred, crystal%gmet, "G")) 1484 qred = qrad * qred_vers 1485 call ifc%fourq(crystal, qred, phfrqs, displ_cart) !, dwdq=dwdq) 1486 if (any(phfrqs < +tol8)) then 1487 num_negw = num_negw + 1; min_negw = min(min_negw, minval(phfrqs)) 1488 end if 1489 !do jj=1,3 1490 ! xval = dot_product(dwdq(:,jj), matmul(crystal%gprimd, qred_vers)) 1491 ! if (xval < zero) num_negw = num_negw + 1 1492 !end do 1493 end do 1494 call xmpi_sum(num_negw, comm, ierr) 1495 xval = min_negw; call xmpi_min(xval, min_negw, comm, ierr) 1496 1497 adiff = sum(abs(cut_phfrq - ref_phfrq)) / (ifc%nqibz * 3 * natom) 1498 if (my_rank == master) then 1499 write(ab_out,"(a,i7,1x,es13.4,4x,i8,1x,es13.4,2x,es13.4)") & 1500 "-",nsphere, adiff * Ha_meV, num_negw, min_negw * Ha_meV, rcut_min 1501 end if 1502 1503 if (num_negw == 0) then 1504 jl = jm ! Replace lower limit 1505 else 1506 ju = jm ! Replace upper limit 1507 end if 1508 end do 1509 1510 ABI_FREE(ref_phfrq) 1511 ABI_FREE(cut_phfrq) 1512 ABI_FREE(save_wghatm) 1513 ABI_FREE(save_atmfrc) 1514 call lgrid%free() 1515 1516 end subroutine ifc_autocutoff
m_ifc/ifc_calcnwrite_nana_terms [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_calcnwrite_nana_terms
FUNCTION
Compute frequencies and phonon displacement for q-->0 in the presence of non-analytical behaviour.
INPUTS
nph2l=Number of qpoints.
qph2l(3,nph2l)=List of phonon wavevector directions along which the non-analytical correction
to the Gamma-point phonon frequencies will be calculated
The direction is in CARTESIAN COORDINATES
qnrml2(nph2l)=Normalization factor.
OUTPUT
(Optional)
phfrq2l(3*crystal%natom,nph2l)=List of phonon frequencies
polarity2l(3,3*crystal%natom,nph2l)=List of mode-polarities
(see Eq.(41) of Veithen et al, PRB71, 125107 (2005) [[cite:Veithen2005]])
NOTES:
This routine should be called by master node and when ifcflag == 1.
SOURCE
2725 subroutine ifc_calcnwrite_nana_terms(ifc, crystal, nph2l, qph2l, & 2726 qnrml2, ncid, phfrq2l, polarity2l) ! optional arguments 2727 2728 !Arguments ------------------------------------ 2729 class(ifc_type),intent(in) :: ifc 2730 integer,intent(in) :: nph2l 2731 integer,optional,intent(in) :: ncid 2732 type(crystal_t),intent(in) :: crystal 2733 !arrays 2734 real(dp),intent(in) :: qph2l(3, nph2l) 2735 real(dp),optional,intent(in) :: qnrml2(nph2l) 2736 real(dp),optional,intent(out) :: phfrq2l(3*crystal%natom,nph2l), polarity2l(3,3*crystal%natom,nph2l) 2737 2738 !Local variables------------------------------- 2739 !scalars 2740 integer :: iatom,idir,imode,iphl2 2741 #ifdef HAVE_NETCDF 2742 integer :: ncerr 2743 #endif 2744 !arrays 2745 real(dp) :: qphnrm(3),qphon(3,3) 2746 real(dp),allocatable :: displ_cart(:,:,:),phfrq(:),d2cart(:,:,:),eigvec(:,:,:),eigval(:) 2747 2748 ! ************************************************************************ 2749 2750 if (nph2l == 0) return 2751 2752 !Now treat the second list of vectors (only at the Gamma point, but can include non-analyticities) 2753 ABI_MALLOC(phfrq, (3*crystal%natom)) 2754 ABI_MALLOC(displ_cart, (2, 3*crystal%natom, 3*crystal%natom)) 2755 ABI_MALLOC(d2cart, (2, 3*ifc%mpert, 3*ifc%mpert)) 2756 ABI_MALLOC(eigvec, (2, 3*crystal%natom, 3*crystal%natom)) 2757 ABI_MALLOC(eigval, (3*crystal%natom)) 2758 2759 ! Before examining every direction or the dielectric tensor, generates the dynamical matrix at gamma 2760 qphon(:,1)=zero; qphnrm(1)=zero 2761 2762 ! Generation of the dynamical matrix in cartesian coordinates 2763 ! Get d2cart using the interatomic forces and the 2764 ! long-range coulomb interaction through Ewald summation 2765 call gtdyn9(ifc%acell,ifc%atmfrc,ifc%dielt,ifc%dipdip, & 2766 ifc%dyewq0,d2cart,crystal%gmet,ifc%gprim,ifc%mpert,crystal%natom, & 2767 ifc%nrpt,qphnrm(1),qphon,crystal%rmet,ifc%rprim,ifc%rpt, & 2768 ifc%trans,crystal%ucvol,ifc%wghatm,crystal%xred,ifc%zeff,ifc%qdrp_cart,ifc%ewald_option,xmpi_comm_self) 2769 2770 #ifdef HAVE_NETCDF 2771 if (present(ncid)) then 2772 iphl2 = 0 2773 call nctk_defwrite_nonana_terms(ncid, iphl2, nph2l, qph2l, crystal%natom, phfrq, displ_cart, "define") 2774 ! Add epsinf, Born effective charges and some useful metadata. 2775 ncerr = nctk_def_arrays(ncid, [ & 2776 nctkarr_t('emacro_cart', "dp", 'number_of_cartesian_directions, number_of_cartesian_directions'), & 2777 nctkarr_t('becs_cart', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, number_of_atoms")], & 2778 defmode=.True.) 2779 NCF_CHECK(ncerr) 2780 ! TODO chneut is missing 2781 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: & 2782 "asr", "dipdip", "symdynmat"]) 2783 NCF_CHECK(ncerr) 2784 NCF_CHECK(nctk_set_datamode(ncid)) 2785 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'emacro_cart'), ifc%dielt)) 2786 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'becs_cart'), ifc%zeff)) 2787 ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: & 2788 "asr", "dipdip", "symdynmat"], & 2789 [ifc%asr, ifc%dipdip, ifc%symdynmat]) 2790 NCF_CHECK(ncerr) 2791 end if 2792 #endif 2793 2794 ! Examine every wavevector of this list 2795 do iphl2=1,nph2l 2796 ! Initialisation of the phonon wavevector 2797 qphon(:,1) = qph2l(:,iphl2) 2798 qphnrm(1)=zero 2799 if(present(qnrml2)) qphnrm(1) = qnrml2(iphl2) 2800 2801 ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix 2802 call dfpt_phfrq(ifc%amu,displ_cart,d2cart,eigval,eigvec,crystal%indsym, & 2803 ifc%mpert,crystal%nsym,crystal%natom,crystal%nsym,crystal%ntypat,phfrq,qphnrm(1),qphon, & 2804 crystal%rprimd,ifc%symdynmat,crystal%symrel,crystal%symafm,crystal%typat,crystal%ucvol) 2805 2806 ! Write the phonon frequencies 2807 !call dfpt_prtph(displ_cart,inp%eivec,inp%enunit,ab_out,natom,phfrq,qphnrm(1),qphon) 2808 2809 if(present(phfrq2l)) phfrq2l(:,iphl2)=phfrq(:) 2810 2811 if(present(polarity2l))then 2812 polarity2l(:,:,iphl2)=zero 2813 do imode=1,3*crystal%natom 2814 do iatom=1,crystal%natom 2815 do idir=1,3 2816 polarity2l(:,imode,iphl2)=polarity2l(:,imode,iphl2)+ifc%zeff(:,idir,iatom)*displ_cart(1,idir+(iatom-1)*3,imode) 2817 enddo 2818 enddo 2819 enddo 2820 endif 2821 2822 #ifdef HAVE_NETCDF 2823 if(present(ncid))then 2824 ! Loop is not MPI-parallelized --> no need for MPI-IO API. 2825 call nctk_defwrite_nonana_terms(ncid, iphl2, nph2l, qph2l, crystal%natom, phfrq, displ_cart, "write") 2826 endif 2827 #endif 2828 end do ! iphl2 2829 2830 ABI_FREE(phfrq) 2831 ABI_FREE(displ_cart) 2832 ABI_FREE(d2cart) 2833 ABI_FREE(eigvec) 2834 ABI_FREE(eigval) 2835 2836 end subroutine ifc_calcnwrite_nana_terms
m_ifc/ifc_fourq [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_fourq
FUNCTION
Compute the phonon frequencies and the group velocities at the specified q-point by performing a Fourier transform on the IFCs matrix in real space.
INPUTS
Ifc<type(ifc_type)>=Object containing the dynamical matrix and the IFCs. Crystal<type(crystal_t)> = Information on the crystalline structure. qpt(3)=q-point in reduced coordinates (unless nanaqdir is specified) [nanaqdir]=If present, the qpt will be treated as a vector specifying the direction in q-space along which the non-analytic behaviour of the dynamical matrix will be treated. Possible values: "cart" if qpt defines a direction in Cartesian coordinates "reduced" if qpt defines a direction in reduced coordinates [comm]: MPI communicator
OUTPUT
phfrq(3*natom) = Phonon frequencies in Hartree displ_cart(2,3,natom,3*natom) = Phonon displacement in Cartesian coordinates [out_d2cart(2,3,3*natom,3,3*natom)] = The (interpolated) dynamical matrix for this q-point [out_eigvec(2*3*natom*3*natom) = The (interpolated) eigenvectors of the dynamical matrix in Cartesian coords. [out_displ_red(2*3*natom*3*natom) = The (interpolated) displacement in reduced coordinates. [dwdq(3,3*natom)] = Group velocities i.e. d(omega(q))/dq in Cartesian coordinates.
SOURCE
923 subroutine ifc_fourq(ifc, crystal, qpt, phfrq, displ_cart, & 924 nanaqdir, comm, & ! Optional [in] 925 out_d2cart, out_eigvec, out_displ_red, dwdq) ! Optional [out] 926 927 !Arguments ------------------------------------ 928 !scalars 929 class(ifc_type),intent(in) :: Ifc 930 character(len=*),optional,intent(in) :: nanaqdir 931 type(crystal_t),intent(in) :: Crystal 932 integer,optional,intent(in) :: comm 933 !arrays 934 real(dp),intent(in) :: qpt(3) 935 real(dp),intent(out) :: displ_cart(2,3,Crystal%natom,3*Crystal%natom) 936 real(dp),intent(out) :: phfrq(3*Crystal%natom) 937 real(dp),optional,intent(out) :: out_d2cart(2,3,Crystal%natom,3,Crystal%natom) 938 real(dp),optional,intent(out) :: out_eigvec(2,3,Crystal%natom,3*Crystal%natom) 939 real(dp),optional,intent(out) :: out_displ_red(2,3,Crystal%natom,3*Crystal%natom) 940 real(dp),optional,intent(out) :: dwdq(3,3*crystal%natom) 941 942 !Local variables------------------------------- 943 !scalars 944 integer :: natom, comm_ 945 real(dp) :: qphnrm 946 !arrays 947 real(dp) :: my_qpt(3),eigvec(2,3,Crystal%natom,3*Crystal%natom),eigval(3*Crystal%natom) 948 real(dp) :: d2cart(2,3,Ifc%mpert,3,Ifc%mpert),tsec(2) 949 950 ! ************************************************************************ 951 952 ! Keep track of total time spent. 953 call timab(1748, 1, tsec) 954 955 natom = Crystal%natom 956 ! TODO: Rewrite and Parallelize ewald9 in gtdyn9 957 comm_ = xmpi_comm_self; if (present(comm)) comm_ = comm 958 959 ! Use my_qpt because dfpt_phfrq can change the q-point (very bad design) 960 qphnrm = one; my_qpt = qpt 961 962 if (present(nanaqdir)) then 963 ! This will break backward compatibility because qpt is **always** in reduced coordinates. 964 ! while dfpt_phfrq assume cartesian coordinates !!!!!!!!!!! 965 ! It does not make sense to change API just to treat this particular case 966 ! We should **alwayse use q-points in reduced coordinates. 967 qphnrm = zero 968 select case (nanaqdir) 969 case ("reduced") 970 ! Convert to Cartesian. 971 my_qpt = matmul(Crystal%gprimd, qpt) 972 case ("cart") 973 continue 974 case default 975 ABI_ERROR(sjoin("Wrong value for nanaqdir:", nanaqdir)) 976 end select 977 end if 978 979 ! The dynamical matrix d2cart is calculated here: 980 call gtdyn9(Ifc%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart,Crystal%gmet,Ifc%gprim,Ifc%mpert,natom,& 981 Ifc%nrpt,qphnrm,my_qpt,Crystal%rmet,Ifc%rprim,Ifc%rpt,Ifc%trans,Crystal%ucvol,Ifc%wghatm,Crystal%xred,Ifc%zeff,& 982 Ifc%qdrp_cart,Ifc%ewald_option,comm_,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad) 983 984 ! Calculate the eigenvectors and eigenvalues of the dynamical matrix 985 call dfpt_phfrq(Ifc%amu,displ_cart,d2cart,eigval,eigvec,Crystal%indsym,& 986 Ifc%mpert,Crystal%nsym,natom,Crystal%nsym,Crystal%ntypat,phfrq,qphnrm,my_qpt,& 987 Crystal%rprimd,Ifc%symdynmat,Crystal%symrel,Crystal%symafm,Crystal%typat,Crystal%ucvol) 988 989 ! OmegaSRLR: Perform decomposition of dynamical matrix 990 !if (srlr==1) call omega_decomp(amu,natom,ntypat,typat,dynmatfull,dynmatsr,dynmatlr,iqpt,nqpt,eigvec) 991 992 ! Return the interpolated dynamical matrix and the eigenvector for this q-point 993 if (present(out_d2cart)) out_d2cart = d2cart(:,:,:natom,:,:natom) 994 if (present(out_eigvec)) out_eigvec = eigvec 995 996 ! Return phonon displacement in reduced coordinates. 997 if (present(out_displ_red)) call phdispl_cart2red(natom, crystal%gprimd, displ_cart, out_displ_red) 998 999 ! Option to get vectors in reduced coordinates? 1000 !call phdispl_cart2red(natom, crystal%gprimd, out_eigvec, out_eigvec_red) 1001 1002 ! Compute group velocities. 1003 if (present(dwdq)) call ifc%get_dwdq(crystal, my_qpt, phfrq, eigvec, dwdq, comm_) 1004 1005 call timab(1748, 2, tsec) 1006 1007 end subroutine ifc_fourq
m_ifc/ifc_free [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_free
FUNCTION
Deallocate memory for the ifc_type structure
SOURCE
261 subroutine ifc_free(ifc) 262 263 !Arguments ------------------------------------ 264 class(ifc_type),intent(inout) :: ifc 265 266 ! ************************************************************************ 267 268 ABI_SFREE(ifc%amu) 269 ABI_SFREE(ifc%atmfrc) 270 ABI_SFREE(ifc%cell) 271 ABI_SFREE(ifc%ewald_atmfrc) 272 ABI_SFREE(ifc%short_atmfrc) 273 ABI_SFREE(ifc%qshft) 274 ABI_SFREE(ifc%rpt) 275 ABI_SFREE(ifc%wghatm) 276 ABI_SFREE(ifc%rcan) 277 ABI_SFREE(ifc%trans) 278 ABI_SFREE(ifc%dyewq0) 279 ABI_SFREE(ifc%qibz) 280 ABI_SFREE(ifc%wtq) 281 ABI_SFREE(ifc%qbz) 282 ABI_SFREE(ifc%zeff) 283 ABI_SFREE(ifc%qdrp_cart) 284 ABI_SFREE(ifc%dynmat) 285 !ABI_SFREE(ifc%dynmat_lr) 286 287 end subroutine ifc_free
m_ifc/ifc_from_file [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_from_file
FUNCTION
Need to be updated
INPUTS
OUTPUT
SOURCE
742 subroutine ifc_from_file(ifc, dielt,filename,natom,ngqpt,nqshift,qshift,ucell_ddb,zeff,qdrp_cart,comm) 743 744 !Arguments ------------------------------------ 745 !scalars 746 class(ifc_type),intent(out) :: Ifc 747 integer,intent(in) :: nqshift,comm 748 integer,intent(inout) :: natom 749 !arrays 750 integer,intent(in) :: ngqpt(3) 751 real(dp),intent(in) :: qshift(3,nqshift) 752 character(len=*),intent(in) :: filename 753 real(dp),intent(inout) :: dielt(3,3) 754 real(dp),allocatable,intent(inout) :: zeff(:,:,:) 755 real(dp),allocatable,intent(inout) :: qdrp_cart(:,:,:,:) 756 type(crystal_t),intent(out) :: ucell_ddb 757 758 !Local variables ------------------------- 759 !scalars 760 integer :: dipdip,i,iblok,iblok_tmp 761 logical :: file_exists 762 character(len=500) :: msg 763 type(ddb_type) :: ddb 764 type(ddb_hdr_type) :: ddb_hdr 765 766 !****************************************************************** 767 768 !check if ddb file exists 769 inquire(file=filename, exist=file_exists) 770 771 if (file_exists .eqv. .true.)then 772 773 !Reading the ddb 774 call ddb%from_file(filename, ddb_hdr, ucell_ddb, comm) 775 call ddb_hdr%free() 776 777 natom = ddb%natom 778 779 else 780 ABI_ERROR(sjoin("File:", filename, "is not present in the directory")) 781 end if 782 783 ! Get Dielectric Tensor and Effective Charges 784 ! (initialized to one_3D and zero if the derivatives are not available in the DDB file) 785 ABI_MALLOC(zeff,(3,3,natom)) 786 iblok = ddb%get_dielt_zeff(ucell_ddb,1,1,0,dielt,zeff) 787 788 ! Try to get dielt, in case just the DDE are present 789 if (iblok == 0) then 790 iblok_tmp = ddb%get_dielt(1,dielt) 791 end if 792 793 ABI_MALLOC(qdrp_cart,(3,3,3,natom)) 794 iblok = ddb%get_quadrupoles(ddb_hdr%ddb_version,1,3,qdrp_cart) 795 796 ! ifc to be calculated for interpolation 797 write(msg, '(a,a,(80a),a,a,a,a)' ) ch10,('=',i=1,80),ch10,ch10,' Calculation of the interatomic forces ',ch10 798 call wrtout(std_out,msg,'COLL') 799 call wrtout(ab_out,msg,'COLL') 800 if ((maxval(abs(zeff)) .lt. tol10) .OR. (maxval(dielt) .gt. 100000.0)) then 801 dipdip=0 802 else 803 dipdip=1 804 end if 805 call ifc%init(ucell_ddb,ddb,1,1,1,dipdip,1,ngqpt,nqshift,qshift,dielt,zeff,qdrp_cart,0,0.0_dp,0,1,comm) 806 807 ! Free them all 808 call ddb%free() 809 810 end subroutine ifc_from_file
m_ifc/ifc_get_dwdq [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_get_dwdq
FUNCTION
Compute phonon group velocities at an arbitrary q-point.
INPUTS
ifc<ifc_type>=Object containing the dynamical matrix and the IFCs. crystal<crystal_t> = Information on the crystalline structure. qpt(3)=q-point in reduced coordinates. eigvec(2*3*natom*3*natom) = The eigenvectors of the dynamical matrix. comm: MPI communicator
OUTPUT
dwdq(3,3*natom) = Group velocities e.g. d(omega(q))/dq in Cartesian coordinates.
NOTES
Using:
D(q) u(q,nu) = w(q, nu)**2 and <u(q,nu) | u(q,nu')> = \delta_{nu, nu'}
one can show, using the Hellman-Feynman theorem, that:
\nabla_q w(q, nu) = 1/(2 w(q, nu)) <u(q, nu)| \nabla_q D(q) | u(q, nu)>
SOURCE
1038 subroutine ifc_get_dwdq(ifc, cryst, qpt, phfrq, eigvec, dwdq, comm) 1039 1040 !Arguments ------------------------------------ 1041 !scalars 1042 class(ifc_type),intent(in) :: ifc 1043 type(crystal_t),intent(in) :: cryst 1044 integer,intent(in) :: comm 1045 !arrays 1046 real(dp),intent(in) :: qpt(3) 1047 real(dp),intent(in) :: phfrq(3*cryst%natom) 1048 real(dp),intent(in) :: eigvec(2,3*cryst%natom,3*cryst%natom) 1049 real(dp),intent(out) :: dwdq(3,3*cryst%natom) 1050 1051 !Local variables------------------------------- 1052 !scalars 1053 !integer,save :: enough=0 1054 integer,parameter :: nqpt1=1,option2=2,sumg0=0 1055 integer :: ii,nu,natom3,jj 1056 real(dp) :: hh 1057 !arrays 1058 real(dp) :: dddq(2,3*cryst%natom,3*cryst%natom,3),dot(2),qfd(3) 1059 real(dp) :: omat(2,3*cryst%natom,3*cryst%natom) 1060 real(dp) :: dyew(2,3*cryst%natom,3*cryst%natom) 1061 1062 ! ************************************************************************ 1063 1064 ABI_UNUSED((/comm/)) 1065 1066 natom3 = cryst%natom * 3 1067 1068 ! Generate the analytical part from the interatomic forces 1069 call dynmat_dq(qpt, cryst%natom, ifc%gprim, ifc%nrpt, ifc%rpt, ifc%atmfrc, ifc%wghatm, dddq) 1070 1071 ! The analytical dynamical matrix dq has been generated 1072 ! in the normalized canonical coordinate system. Now, the 1073 ! phase is modified, in order to recover the usual (xred) coordinate of atoms. 1074 do ii=1,3 1075 call dymfz9(dddq(:,:,:,ii), cryst%natom, nqpt1, ifc%gprim, option2, qpt, ifc%trans) 1076 dddq(:,:,:,ii) = dddq(:,:,:,ii) * ifc%acell(ii) 1077 end do 1078 1079 if (ifc%dipdip == 1.or.ifc%dipquad == 1.or.ifc%quadquad == 1) then 1080 ! Add the gradient of the non-analytical part. 1081 ! Note that dddq is in cartesian cordinates. 1082 ! For the time being, the gradient is computed with finite difference and step hh. 1083 ! TODO: should generalize ewald9 to compute dq. 1084 !enough = enough + 1 1085 !if (enough <= 5) ABI_WARNING("phonon velocities with dipdip==1 not yet tested.") 1086 hh = 0.01_dp 1087 do ii=1,3 1088 do jj=-1,1,2 1089 ! qcart --> qred 1090 qfd = zero; qfd(ii) = jj 1091 qfd = matmul(cryst%rprimd, qfd); qfd = qfd / normv(qfd, cryst%gmet, "G") 1092 !write(std_out,*)"normv:",normv(qfd, cryst%gmet, "G") 1093 qfd = qpt + hh * qfd 1094 1095 call ewald9(ifc%acell,ifc%dielt,dyew,cryst%gmet,ifc%gprim,cryst%natom,qfd,& 1096 cryst%rmet,ifc%rprim,sumg0,cryst%ucvol,cryst%xred,ifc%zeff,ifc%qdrp_cart,& 1097 ifc%ewald_option,dipquad=ifc%dipquad,quadquad=ifc%quadquad) 1098 call q0dy3_apply(cryst%natom,ifc%dyewq0,dyew) 1099 dddq(:,:,:,ii) = dddq(:,:,:,ii) + (jj * half / hh) * dyew 1100 end do 1101 end do 1102 end if 1103 1104 do ii=1,3 1105 call massmult_and_breaksym(cryst%natom, cryst%ntypat, cryst%typat, ifc%amu, dddq(:,:,:,ii)) 1106 end do 1107 1108 ! Compute 1/(2w(q)) <u(q)|dD(q)/dq|u(q)> 1109 do ii=1,3 1110 call zgemm('N','N',natom3,natom3,natom3,cone,dddq(:,:,:,ii),natom3,eigvec,natom3,czero,omat,natom3) 1111 do nu=1,natom3 1112 if (abs(phfrq(nu)) > tol12) then 1113 dot = cg_zdotc(natom3, eigvec(1,1,nu), omat(1,1,nu)) 1114 ! abs(w) is needed to get the correct derivative if we have a purely imaginary solution. 1115 dwdq(ii, nu) = dot(1) / (two * abs(phfrq(nu))) 1116 else 1117 dwdq(ii, nu) = zero 1118 end if 1119 end do 1120 end do 1121 1122 end subroutine ifc_get_dwdq
m_ifc/ifc_get_phmesh [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_get_phmesh
FUNCTION
Build linear mesh for phonons.
INPUTS
OUTPUT
SOURCE
1140 subroutine ifc_get_phmesh(ifc, ph_wstep, phmesh_size, phmesh) 1141 1142 !Arguments ------------------------------------ 1143 !scalars 1144 class(ifc_type),intent(in) :: ifc 1145 real(dp),intent(in) :: ph_wstep 1146 integer,intent(out) :: phmesh_size 1147 !arrays 1148 real(dp),allocatable,intent(out) :: phmesh(:) 1149 1150 !****************************************************************** 1151 1152 phmesh_size = nint((ifc%omega_minmax(2) - ifc%omega_minmax(1) ) / ph_wstep) + 1 1153 ABI_MALLOC(phmesh, (phmesh_size)) 1154 phmesh = arth(ifc%omega_minmax(1), ph_wstep, phmesh_size) 1155 1156 end subroutine ifc_get_phmesh
m_ifc/ifc_getiaf [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_getiaf
FUNCTION
Extracts the IFCs needed for the output for one atom in the unit cell. Accumulates the results for writing in the NetCDF file. Prints to the output file
INPUTS
Ifc<type(ifc_type)>=Object containing the dynamical matrix and the IFCs. ifcana= 0 => no analysis of ifc ; 1 => full analysis ifcout= Number of interatomic force constants written in the output file iout=unit number for nice output zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement ia=index of the atom in the unit cell for which the IFCs are being analyzed ra(3)=position of atom ia in cartesian coordinates list(ifcout)=index of permutation for distances from atom ia in ascending order dist(natom,natom,nrpt)=distance from atom ia to atom ib in unit cell irpt. invdlt(3,3)=inverse (transpose) of the dielectric tensor detdlt=determinant of the dielectric tensor
OUTPUT
rsiaf(3,3,ifcout)=list of real space IFCs sriaf(3,3,ifcout)=list of the short range part of the real space IFCs vect(3,3,ifcout)=base vectors for local coordinates (longitudinal/transverse), if ifc_getiaf is able to find a third atom not aligned with the two atoms characterizing the IFC. If no, the second and third vectors are set to zero. indngb(ifcout)=indices in the unit cell of the neighbouring atoms posngb(3,ifcout)=position of the neighbouring atoms in cartesian coordinates output file
NOTES
This routine should be executed by one processor only
SOURCE
2031 subroutine ifc_getiaf(Ifc,ifcana,ifcout,iout,zeff,ia,ra,list,& 2032 & dist,invdlt,detdlt,rsiaf,sriaf,vect,indngb,posngb) 2033 2034 !Arguments ------------------------------- 2035 !scalars 2036 class(ifc_type),intent(inout) :: Ifc 2037 integer,intent(in) :: ia,ifcana,ifcout,iout 2038 real(dp), intent(in) :: detdlt 2039 !arrays 2040 real(dp),intent(in) :: invdlt(3,3),ra(3) 2041 real(dp),intent(in) :: dist(Ifc%natom,Ifc%natom,Ifc%nrpt) 2042 real(dp),intent(in) :: zeff(3,3,Ifc%natom) 2043 integer,intent(in) :: list(Ifc%natom*Ifc%nrpt) 2044 integer,intent(out) :: indngb(ifcout) 2045 real(dp),intent(out) :: rsiaf(3,3,ifcout),sriaf(3,3,ifcout),vect(3,3,ifcout),posngb(3,ifcout) 2046 2047 !Local variables ------------------------- 2048 !scalars 2049 integer :: flag,ib,ii,index,jj,kk,mu,nu,irpt 2050 real(dp) :: ew1,rsq,scprod,trace1,trace2,trace3 2051 real(dp) :: yy,dist1 2052 character(len=500) :: message 2053 !arrays 2054 real(dp) :: ewiaf0(3,3),ewiaf1(3,3),ewloc(3,3),ifcloc(3,3) 2055 real(dp) :: rcart(3),rdiff(3),rsloc(3,3) 2056 real(dp) :: srloc(3,3),vect1(3),vect2(3),vect3(3),work(3),xx(3) 2057 2058 ! ********************************************************************* 2059 2060 if(ifcana==1)then 2061 ! Generate the local coordinate system for the atom ia 2062 index=list(2) 2063 write(std_out,*)index 2064 call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,rcart,Ifc%rprim,Ifc%rpt) 2065 vect2(1)=rcart(1)-ra(1) 2066 vect2(2)=rcart(2)-ra(2) 2067 vect2(3)=rcart(3)-ra(3) 2068 flag=0 2069 do ii=3,Ifc%natom*Ifc%nrpt 2070 index=list(ii) 2071 call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,rcart,Ifc%rprim,Ifc%rpt) 2072 vect1(1)=(rcart(1)-ra(1))-vect2(1) 2073 vect1(2)=(rcart(2)-ra(2))-vect2(2) 2074 vect1(3)=(rcart(3)-ra(3))-vect2(3) 2075 scprod=0.0_dp 2076 do jj=1,3 2077 scprod=scprod+vect1(jj)**2 2078 end do 2079 do jj=1,3 2080 vect1(jj)=vect1(jj)/scprod**0.5 2081 end do 2082 scprod=0.0_dp 2083 do jj=1,3 2084 scprod=scprod+vect2(jj)*vect1(jj) 2085 end do 2086 do jj=1,3 2087 work(jj)=vect2(jj)-vect1(jj)*scprod 2088 end do 2089 scprod=0.0_dp 2090 do jj=1,3 2091 scprod=scprod+work(jj)**2 2092 end do 2093 if(scprod>1.0d-10)then 2094 flag=1 2095 end if 2096 if(flag==1)exit 2097 end do 2098 if(flag==0)then 2099 write(message, '(3a)' )& 2100 & 'Unable to find a third atom not aligned with the two selected ones.',ch10,& 2101 & 'The local analysis (longitudinal/transverse) will not be done. The two transverse vectors are set to zero.' 2102 ABI_WARNING(message) 2103 vect2(:)=zero ; vect3(:)=zero 2104 else 2105 vect2(1)=work(1)/scprod**0.5 2106 vect2(2)=work(2)/scprod**0.5 2107 vect2(3)=work(3)/scprod**0.5 2108 vect3(1)=vect1(2)*vect2(3)-vect1(3)*vect2(2) 2109 vect3(2)=vect1(3)*vect2(1)-vect1(1)*vect2(3) 2110 vect3(3)=vect1(1)*vect2(2)-vect1(2)*vect2(1) 2111 if (iout > 0) then 2112 write(iout, '(a)' )' Third atom defining local coordinates : ' 2113 write(iout, '(a,i4,a,i4)' )' ib = ',ib,' irpt = ',irpt 2114 end if 2115 endif 2116 end if 2117 2118 ! Analysis and output of force constants, ordered with respect to the distance from atom ia 2119 do ii=1,ifcout 2120 index=list(ii) 2121 call canct9(Ifc%acell,Ifc%gprim,ib,index,irpt,Ifc%natom,Ifc%nrpt,Ifc%rcan,posngb(:,ii),Ifc%rprim,Ifc%rpt) 2122 indngb(ii)=ib 2123 dist1=dist(ia,ib,irpt) 2124 if (iout > 0) then 2125 write(iout, '(a)' ) 2126 write(iout, '(i4,a,i6,a,i8)' )ii,' interaction with atom',ib,' cell',irpt 2127 write(iout, '(a,3es16.6)' )' with coordinates ',posngb(1:3,ii)*(one+tol8) 2128 write(iout, '(a,es16.6)' )' and distance ',dist1 2129 end if 2130 2131 if(ifcana==1.and.ii/=1)then 2132 vect1(1)=(posngb(1,ii)-ra(1))/dist1 2133 vect1(2)=(posngb(2,ii)-ra(2))/dist1 2134 vect1(3)=(posngb(3,ii)-ra(3))/dist1 2135 end if 2136 2137 if(Ifc%dipdip==0.or.Ifc%dipquad==1.or.Ifc%quadquad==1)then 2138 ! Get the "total" force constants (=real space FC) 2139 ! without taking into account the dipole-dipole interaction 2140 do mu=1,3 2141 do nu=1,3 2142 rsiaf(mu,nu,ii)=Ifc%atmfrc(mu,ia,nu,ib,irpt) * Ifc%wghatm(ia,ib,irpt) 2143 end do 2144 end do 2145 ! Output of the ifcs in cartesian coordinates 2146 if (iout > 0) then 2147 do nu=1,3 2148 write(iout, '(1x,3f9.5)' )(rsiaf(mu,nu,ii)+tol10,mu=1,3) 2149 ! transfer short range and long range 2150 do mu=1,3 2151 Ifc%short_atmfrc(mu,ia,nu,ib,irpt) = rsiaf(mu,nu,ii) + tol10 2152 end do 2153 2154 end do 2155 end if 2156 2157 if(ifcana==1)then 2158 ! Further analysis 2159 trace1=rsiaf(1,1,ii)+rsiaf(2,2,ii)+rsiaf(3,3,ii) 2160 if (iout > 0) then 2161 write(iout, '(a,f9.5)' ) ' Trace ',trace1+tol10 2162 end if 2163 if(flag==1)then 2164 if(ii/=1)then 2165 call axial9(rsiaf(:,:,ii),vect1,vect2,vect3) 2166 end if 2167 if (iout > 0) then 2168 write(iout, '(a)' )' Transformation to local coordinates ' 2169 write(iout, '(a,3f16.6)' ) ' First local vector :',vect1 2170 write(iout, '(a,3f16.6)' ) ' Second local vector :',vect2 2171 write(iout, '(a,3f16.6)' ) ' Third local vector :',vect3 2172 end if 2173 call ifclo9(rsiaf(:,:,ii),ifcloc,vect1,vect2,vect3) 2174 if (iout > 0) then 2175 do nu=1,3 2176 write(iout, '(1x,3f9.5)' )(ifcloc(mu,nu)+tol10,mu=1,3) 2177 end do 2178 end if 2179 endif ! flag==1 2180 2181 vect(:,1,ii) = vect1 2182 vect(:,2,ii) = vect2 2183 vect(:,3,ii) = vect3 2184 2185 end if ! Further analysis finished 2186 2187 else if(Ifc%dipdip==1)then 2188 2189 !write(iout,'(a)') 2190 !write(iout,'(a)')' Enter dipdip section, for debugging' 2191 !write(iout,'(a)') 2192 2193 ! Get the Coulomb part 2194 do jj=1,3 2195 rdiff(jj)=ra(jj)-posngb(jj,ii) 2196 end do 2197 rsq=0.0_dp 2198 xx(1:3)=0.0_dp 2199 do jj=1,3 2200 do kk=1,3 2201 ewiaf0(jj,kk)=0.0_dp 2202 rsq=rsq+rdiff(jj)*invdlt(kk,jj)*rdiff(kk) 2203 xx(kk)=xx(kk)+invdlt(kk,jj)*rdiff(jj) 2204 end do 2205 end do 2206 yy=sqrt(rsq) 2207 ! Avoid zero denominators in term: 2208 if (sqrt(rsq)>=tol12) then 2209 do mu=1,3 2210 do nu=1,3 2211 ewiaf0(mu,nu)=(-3*xx(nu)*xx(mu)+invdlt(nu,mu)*yy**2)/yy**5/sqrt(detdlt) 2212 end do 2213 end do 2214 else 2215 if (ia/=ib)then 2216 write(message, '(a,a,a,a,a,i5,a,i5,a)' )& 2217 & 'The distance between two atoms vanishes.',ch10,& 2218 & 'This is not allowed.',ch10,& 2219 & 'Action: check the input for the atoms number',ia,' and',ib,'.' 2220 ABI_ERROR(message) 2221 end if 2222 end if 2223 2224 ! Take into account the effective charge tensor 2225 do mu=1,3 2226 do nu=1,3 2227 ew1=zero 2228 if(ii==1)then 2229 ew1=-Ifc%dyewq0(mu,nu,ia) 2230 end if 2231 do jj=1,3 2232 do kk=1,3 2233 ew1=ew1+zeff(jj,mu,ia)*(zeff(kk,nu,ib)* ewiaf0(jj,kk)) 2234 end do 2235 end do 2236 ewiaf1(mu,nu)=ew1 2237 end do 2238 end do 2239 ! Get the short-range force constants and the 2240 ! "total" force constants (=real space FC) 2241 do mu=1,3 2242 do nu=1,3 2243 sriaf(mu,nu,ii)=Ifc%atmfrc(mu,ia,nu,ib,irpt)* Ifc%wghatm(ia,ib,irpt) 2244 rsiaf(mu,nu,ii)=ewiaf1(mu,nu)+sriaf(mu,nu,ii) 2245 end do 2246 end do 2247 2248 ! Output of the results 2249 if (iout > 0) then 2250 do nu=1,3 2251 write(iout, '(1x,3(3f9.5,1x))' )& 2252 & (rsiaf(mu,nu,ii) +tol10,mu=1,3),& 2253 & (ewiaf1(mu,nu)+tol10,mu=1,3),& 2254 & (sriaf(mu,nu,ii) +tol10,mu=1,3) 2255 2256 ! transfer short range and long range 2257 do mu=1,3 2258 Ifc%short_atmfrc(mu,ia,nu,ib,irpt) = sriaf(mu,nu,ii) + tol10 2259 Ifc%ewald_atmfrc(mu,ia,nu,ib,irpt) = ewiaf1(mu,nu) + tol10 2260 end do 2261 end do 2262 end if 2263 2264 if(ifcana==1)then 2265 ! Further analysis 2266 if (iout > 0) then 2267 write(iout, '(a)' )' Traces (and ratios) :' 2268 end if 2269 trace1=rsiaf(1,1,ii)+rsiaf(2,2,ii)+rsiaf(3,3,ii) 2270 trace2=ewiaf1(1,1)+ewiaf1(2,2)+ewiaf1(3,3) 2271 trace3=sriaf(1,1,ii)+sriaf(2,2,ii)+sriaf(3,3,ii) 2272 if (iout > 0) then 2273 write(iout,'(3(f9.5,17x))')trace1+tol10,trace2+tol10,trace3+tol10 2274 write(iout,'(3(f9.5,17x))')1.0,(trace2+tol10)/(trace1+tol10),(trace3+tol10)/(trace1+tol10) ! 2275 end if 2276 2277 if(flag==1)then 2278 if(ii/=1)then 2279 call axial9(rsiaf(:,:,ii),vect1,vect2,vect3) 2280 end if 2281 if (iout > 0) then 2282 write(iout, '(a)' )' Transformation to local coordinates ' 2283 write(iout, '(a,3f16.6)' )' First local vector :',vect1 2284 write(iout, '(a,3f16.6)' )' Second local vector :',vect2 2285 write(iout, '(a,3f16.6)' )' Third local vector :',vect3 2286 end if 2287 call ifclo9(rsiaf(:,:,ii),rsloc,vect1,vect2,vect3) 2288 call ifclo9(ewiaf1,ewloc,vect1,vect2,vect3) 2289 call ifclo9(sriaf(:,:,ii),srloc,vect1,vect2,vect3) 2290 if (iout > 0) then 2291 do nu=1,3 2292 write(iout, '(1x,3(3f9.5,1x))' )& 2293 & (rsloc(mu,nu)+tol10,mu=1,3),& 2294 & (ewloc(mu,nu)+tol10,mu=1,3),& 2295 & (srloc(mu,nu)+tol10,mu=1,3) 2296 end do 2297 if(ii/=1)then 2298 write(iout, '(a)' )' Ratio with respect to the longitudinal ifc' 2299 else 2300 write(iout, '(a)' )' Ratio with respect to the (1,1) element' 2301 end if 2302 do nu=1,3 2303 write(iout, '(1x,3(3f9.5,1x))' )& 2304 & (rsloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3),& 2305 & (ewloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3),& 2306 & (srloc(mu,nu)/rsloc(1,1)+tol10,mu=1,3) 2307 end do 2308 end if 2309 endif ! flag==1 2310 2311 vect(:,1,ii) = vect1 2312 vect(:,2,ii) = vect2 2313 vect(:,3,ii) = vect3 2314 2315 end if ! Further analysis finished 2316 end if ! End the condition on dipdip 2317 end do ! End loop over all atoms in BigBox: 2318 2319 end subroutine ifc_getiaf
m_ifc/ifc_init [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_init
FUNCTION
Initialize the dynamical matrix as well as the IFCs. taking into account the dipole-dipole, dipole-quadrupole and quadrupole-quadrupole interaction.
INPUTS
crystal<type(crystal_t)> = Information on the crystalline structure. ddb<type(ddb_type)> = Database with derivatives. brav=bravais lattice (1 or -1=simple lattice, 2=face centered lattice, 3=centered lattice, 4=hexagonal lattice) asr= Option for the imposition of the ASR 0 => no ASR, 1 => modify "asymmetrically" the diagonal element 2 => modify "symmetrically" the diagonal element symdynmat=if 1, (re)symmetrize the dynamical matrix, except if Gamma wavevector with electric field added. dipdip= if 0, no dipole-dipole interaction was subtracted in atmfrc if 1, atmfrc has been build without dipole-dipole part rfmeth = 1 if non-stationary block 2 if stationary block 3 if third order derivatives dielt(3,3)=dielectric tensor. zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement prtsrlr: TODO: TO BE REMOVED enunit: TODO: TO BE REMOVED dielt(3,3)=dielectric tensor ngqpt_in = input values of ngqpt nqshft=Number of shifths in q-grid. q1shft(3,nqshft)=Shifts for q-grid nsphere=number of atoms to be included in the cut-off sphere for interatomic force constant. 0: maximum extent allowed by the grid. > 0: Apply cutoff -1: Analyze the effect of different nsphere values on the phonon spectrum, in particular the frequencies around gamma. rifcsph=radius for cutoff of IFC. comm=MPI communicator. [Ifc_coarse]=Optional. [dipquad] = if 1, atmfrc has been build without dipole-quadrupole part [quadquad] = if 1, atmfrc has been build without quadrupole-quadrupole part
OUTPUT
Ifc<ifc_type>=Object containing the dynamical matrix and the IFCs.
SOURCE
341 subroutine ifc_init(ifc,crystal,ddb,brav,asr,symdynmat,dipdip,& 342 rfmeth,ngqpt_in,nqshft,q1shft,dielt,zeff,qdrp_cart,nsphere,rifcsph,& 343 prtsrlr,enunit, & ! TODO: TO BE REMOVED 344 comm, & 345 Ifc_coarse,dipquad,quadquad) ! Optional 346 347 !Arguments ------------------------------------ 348 class(ifc_type),intent(inout) :: Ifc 349 integer,intent(in) :: asr,brav,dipdip,symdynmat,nqshft,rfmeth,nsphere,comm 350 real(dp),intent(in) :: rifcsph 351 type(crystal_t),intent(in) :: Crystal 352 type(ddb_type),intent(in) :: ddb 353 type(ifc_type),optional,intent(in) :: Ifc_coarse 354 integer,optional,intent(in) :: dipquad, quadquad 355 356 !arrays 357 integer,intent(in) :: ngqpt_in(3) 358 real(dp),intent(in) :: q1shft(3,nqshft) 359 real(dp),intent(in) :: dielt(3,3),zeff(3,3,Crystal%natom) 360 real(dp),intent(in) :: qdrp_cart(3,3,3,Crystal%natom) 361 !anaddb variables (TO BE REMOVED) 362 integer,intent(in) :: prtsrlr,enunit 363 !end anaddb variables 364 365 !Local variables ------------------------- 366 !scalars 367 integer,parameter :: timrev1=1,iout0=0,chksymbreak0=0 368 integer :: mpert,iout,iqpt,mqpt,nsym,ntypat,iq_ibz,iq_bz,ii,natom 369 integer :: nqbz,option,plus,sumg0,irpt,irpt_new 370 integer :: nprocs,my_rank,my_ierr,ierr 371 real(dp),parameter :: qphnrm=one 372 real(dp) :: xval,cpu,wall,gflops,rcut_min 373 real(dp) :: r_inscribed_sphere,toldist 374 character(len=500*4) :: msg 375 type(ifc_type) :: ifc_tmp 376 !arrays 377 integer :: ngqpt(9),qptrlatt(3,3) 378 integer,allocatable :: qmissing(:),ibz2bz(:),bz2ibz_smap(:,:) 379 real(dp) :: gprim(3,3),rprim(3,3),qpt(3),rprimd(3,3), gprim_tmp(3,3), rprim_tmp(3,3) 380 real(dp):: rcan(3,Crystal%natom),trans(3,Crystal%natom),dyewq0(3,3,Crystal%natom) 381 real(dp) :: displ_cart(2*3*Crystal%natom*3*Crystal%natom) 382 real(dp) :: phfrq(3*Crystal%natom) 383 real(dp) :: eigvec(2,3,Crystal%natom,3,Crystal%natom) 384 real(dp),allocatable :: dyew(:,:,:,:,:),out_d2cart(:,:,:,:,:) 385 real(dp),allocatable :: dynmatfull(:,:,:,:,:,:),dynmat_sr(:,:,:,:,:,:),dynmat_lr(:,:,:,:,:,:) ! for OmegaSRLR 386 real(dp),allocatable :: wtq(:),wtq_folded(:),qbz(:,:) 387 388 !****************************************************************** 389 390 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 391 call cwtime(cpu, wall, gflops, "start") 392 393 mpert = ddb%mpert 394 iout = ab_out 395 396 rprim = ddb%rprim; gprim = ddb%gprim 397 398 nsym = Crystal%nsym 399 natom = Crystal%natom 400 ntypat = Crystal%ntypat 401 rprimd = Crystal%rprimd 402 403 ngqpt=0; ngqpt(1:3)=ngqpt_in(1:3) 404 405 ! Copy important parameters in Ifc 406 Ifc%natom = natom 407 Ifc%mpert = mpert 408 Ifc%asr = asr 409 Ifc%brav = brav 410 Ifc%dipdip = abs(dipdip) 411 Ifc%dipquad=0; if (present(dipquad)) Ifc%dipquad = dipquad 412 Ifc%quadquad=0; if (present(quadquad)) Ifc%quadquad = quadquad 413 Ifc%symdynmat = symdynmat 414 Ifc%nqshft = nqshft 415 call alloc_copy(q1shft(:,1:Ifc%nqshft),Ifc%qshft) 416 Ifc%ngqpt = ngqpt_in(1:3) 417 Ifc%rprim = ddb%rprim 418 Ifc%gprim = ddb%gprim 419 Ifc%acell = ddb%acell 420 Ifc%ewald_option = 0; if (dipdip < 0) Ifc%ewald_option = 1 !HM TODO: expose this in the init? 421 422 ! Check if the rprim are coherent with the choice used in the interatomic forces generation 423 call chkrp9(Ifc%brav,rprim) 424 425 ! Compute dyewq0, the correction to be applied to the Ewald, see Eq.(71) of PRB55, 10355 (1997). 426 dyewq0 = zero 427 if ((Ifc%dipdip==1.or.Ifc%dipquad==1.or.Ifc%quadquad==1).and. (Ifc%asr==1.or.Ifc%asr==2)) then 428 ! Calculation of the non-analytical part for q=0 429 sumg0=0 430 qpt(:)=zero 431 ABI_MALLOC(dyew,(2,3,natom,3,natom)) 432 if (Ifc%dipquad==1.or.Ifc%quadquad==1) then 433 call ewald9(ddb%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,& 434 Crystal%xred,zeff,qdrp_cart,option=ifc%ewald_option,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad) 435 else 436 call ewald9(ddb%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,& 437 Crystal%xred,zeff,qdrp_cart,option=ifc%ewald_option) 438 end if 439 440 call q0dy3_calc(natom,dyewq0,dyew,Ifc%asr) 441 ABI_FREE(dyew) 442 end if 443 444 ! Sample the Brillouin zone 445 option=1 446 qptrlatt = 0; qptrlatt(1,1)=ngqpt(1); qptrlatt(2,2)=ngqpt(2); qptrlatt(3,3)=ngqpt(3) 447 mqpt=ngqpt(1)*ngqpt(2)*ngqpt(3)*nqshft 448 if (Ifc%brav==2) mqpt=mqpt/2 449 if (Ifc%brav==3) mqpt=mqpt/4 450 451 ABI_MALLOC(qbz,(3,mqpt)) 452 call smpbz(Ifc%brav,ab_out,qptrlatt,mqpt,nqbz,nqshft,option,q1shft,qbz) 453 454 ! Find the irreducible zone (qibz) 455 ABI_MALLOC(ibz2bz, (nqbz)) 456 ABI_MALLOC(wtq_folded, (nqbz)) 457 ABI_MALLOC(wtq, (nqbz)) 458 wtq = one / nqbz ! Weights sum up to one 459 ABI_MALLOC(bz2ibz_smap, (6, nqbz)) 460 461 ! FIXME: timrev1 should be set to 0 if TR cannot be used 462 call symkpt(chksymbreak0,crystal%gmet,ibz2bz,iout0,qbz,nqbz,ifc%nqibz,crystal%nsym,& 463 crystal%symrec,timrev1,wtq,wtq_folded, bz2ibz_smap, xmpi_comm_self) 464 465 ABI_FREE(bz2ibz_smap) 466 467 ABI_MALLOC(ifc%qibz, (3,ifc%nqibz)) 468 ABI_MALLOC(ifc%wtq, (ifc%nqibz)) 469 do iq_ibz=1,ifc%nqibz 470 ifc%qibz(:,iq_ibz) = qbz(:, ibz2bz(iq_ibz)) 471 ifc%wtq(iq_ibz) = wtq_folded(ibz2bz(iq_ibz)) 472 end do 473 ABI_FREE(ibz2bz) 474 ABI_FREE(wtq_folded) 475 ABI_FREE(wtq) 476 477 ABI_MALLOC(Ifc%dynmat,(2,3,natom,3,natom,nqbz)) 478 479 ! This is needed to preserve the behavior of the old implementation with canonical coordinate. 480 if (Ifc%brav == 1) then 481 gprim_tmp = Crystal%gprimd 482 rprim_tmp = rprimd 483 else 484 gprim_tmp = gprim 485 rprim_tmp = rprim 486 endif 487 488 ! Find symmetrical dynamical matrices 489 if (.not.present(Ifc_coarse)) then 490 491 ! Each q-point in the BZ mush be the symmetrical of one of the qpts in the ddb file. 492 ! SP - gprimd and rprimd is required instead of gprim and rprim for non-diagonal supercells. 493 call symdm9(ddb, & 494 Ifc%dynmat,gprim_tmp,Crystal%indsym,mpert,natom,nqbz,nsym,rfmeth,rprim_tmp,qbz,& 495 Crystal%symrec, Crystal%symrel, comm) 496 497 else 498 ! Symmetrize the qpts in the BZ using the q-points in the ddb. 499 ! Then use Ifc_coarse to fill the missing entries with Fourier interpolated matrices. 500 ! 501 ! TODO: The previous version of refineblk was hacking the DDB database to add the q-points in the **IBZ** 502 ! Then D(q) was symmetrized in symdm9. This version avoids the symmetrization: the q-points 503 ! in the BZ that are not in the coarse q-mesh are obtained by an explicit FT. 504 ! This means that the final D(q) may break some symmetry in q-space if the FT does not preserve it. 505 ! The most elegant approach would be to get D(q_ibz) via FT if q_ibz is not in the coarse mesh and then 506 ! call symdm9 to get D(q) for each q point in the star of q_ibz. 507 call wrtout(std_out,"Will fill missing qpoints in the full BZ using the coarse q-mesh","COLL") 508 509 call symdm9(ddb, & 510 Ifc%dynmat,gprim_tmp,Crystal%indsym,mpert,natom,nqbz,nsym,rfmeth,rprim_tmp,qbz,& 511 Crystal%symrec,Crystal%symrel,comm, qmissing=qmissing) 512 513 ! Compute dynamical matrix with Fourier interpolation on the coarse q-mesh. 514 write(msg,"(a,i0,a)")"Will use Fourier interpolation to construct D(q) for ",size(qmissing)," q-points" 515 call wrtout(std_out,msg) 516 517 ABI_MALLOC(out_d2cart, (2,3,natom,3,natom)) 518 do ii=1,size(qmissing) 519 iq_bz = qmissing(ii) 520 qpt = qbz(:,iq_bz) 521 ! TODO: check dipdip option and phase, but I think this is correct! 522 call Ifc_coarse%fourq(Crystal,qpt,phfrq,displ_cart,out_d2cart=out_d2cart) 523 Ifc%dynmat(:,:,:,:,:,iq_bz) = out_d2cart 524 end do 525 526 ABI_FREE(qmissing) 527 ABI_FREE(out_d2cart) 528 end if 529 530 ! OmegaSRLR: Store full dynamical matrix for decomposition into short- and long-range parts 531 ABI_MALLOC(dynmatfull,(2,3,natom,3,natom,nqbz)) 532 dynmatfull=Ifc%dynmat 533 534 if (Ifc%dipdip==1.or.Ifc%dipquad==1.or.Ifc%quadquad==1) then 535 ! Take off the dipole-dipole part of the dynamical matrix 536 call wrtout(std_out, " Will extract the dipole-dipole part for every wavevector") 537 ABI_MALLOC(dyew,(2,3,natom,3,natom)) 538 539 do iqpt=1,nqbz 540 if (mod(iqpt, nprocs) /= my_rank) then ! mpi-parallelism 541 ifc%dynmat(:,:,:,:,:,iqpt) = zero; cycle 542 end if 543 qpt(:)=qbz(:,iqpt) 544 sumg0=0 545 if (Ifc%dipquad==1.or.Ifc%quadquad==1) then 546 call ewald9(ddb%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,& 547 Crystal%xred,zeff,qdrp_cart,option=ifc%ewald_option,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad) 548 else 549 call ewald9(ddb%acell,dielt,dyew,Crystal%gmet,gprim,natom,qpt,Crystal%rmet,rprim,sumg0,Crystal%ucvol,& 550 Crystal%xred,zeff,qdrp_cart,option=ifc%ewald_option) 551 end if 552 call q0dy3_apply(natom,dyewq0,dyew) 553 plus=0 554 ! Implement Eq.(76) of Gonze&Lee PRB 55, 10355 (1997) [[cite:Gonze1997a]], possibly generalized for quadrupoles 555 call nanal9(dyew,Ifc%dynmat,iqpt,natom,nqbz,plus) 556 end do 557 558 call xmpi_sum(ifc%dynmat, comm, ierr) 559 ABI_FREE(dyew) 560 end if 561 562 ! OmegaSRLR: Store the short-range dynmat and compute long-range as difference 563 ABI_MALLOC(dynmat_sr,(2,3,natom,3,natom,nqbz)) 564 ABI_MALLOC(dynmat_lr,(2,3,natom,3,natom,nqbz)) 565 dynmat_sr=Ifc%dynmat 566 dynmat_lr=dynmatfull-dynmat_sr 567 568 ! Now, take care of the remaining part of the dynamical matrix 569 ! Move to canonical normalized coordinates 570 call canat9(Ifc%brav,natom,rcan,rprim,trans,Crystal%xred) 571 572 ! Multiply the dynamical matrix by a phase shift 573 option=1 574 call dymfz9(Ifc%dynmat,natom,nqbz,gprim,option,qbz,trans) 575 576 ! Create the Big Box of R vectors in real space and compute the number of points (cells) in real space 577 call make_bigbox(Ifc%brav,ifc_tmp%cell,ngqpt,nqshft,rprim,ifc_tmp%nrpt,ifc_tmp%rpt) 578 579 ! Weights associated to these R points and to atomic pairs 580 ABI_MALLOC(ifc_tmp%wghatm, (natom, natom, ifc_tmp%nrpt)) 581 582 ! HM: this tolerance is highly dependent on the compilation/architecture 583 ! numeric errors in the DDB text file. Try a few tolerances and check whether all the weights are found. 584 toldist = tol8 585 do while (toldist <= tol6) 586 ! Note ngqpt(9) with intent(inout)! 587 call wght9(Ifc%brav,gprim,natom,ngqpt,nqbz,nqshft,ifc_tmp%nrpt,q1shft,rcan,& 588 ifc_tmp%rpt,rprimd,toldist,r_inscribed_sphere,ifc_tmp%wghatm,my_ierr) 589 call xmpi_max(my_ierr, ierr, comm, ii) 590 if (ierr > 0) toldist = toldist * 10 591 if (ierr == 0) exit 592 end do 593 594 if (ierr > 0) then 595 write(msg, '(3a,es14.4,2a,i0, 14a)' ) & 596 'The sum of the weight is not equal to nqpt.',ch10,& 597 'The sum of the weights is: ',sum(ifc_tmp%wghatm),ch10,& 598 'The number of q points is: ',nqbz, ch10, & 599 'This might have several sources.',ch10,& 600 'If toldist is larger than 1.0e-8, the atom positions might be loose.',ch10,& 601 'and the q point weights not computed properly.',ch10,& 602 'Action: make input atomic positions more symmetric.',ch10,& 603 'Otherwise, you might increase "buffer" in m_dynmat.F90 see bigbx9 subroutine and recompile.',ch10,& 604 'Actually, this can also happen when ngqpt is 0 0 0,',ch10,& 605 'if abs(brav) /= 1, in this case you should change brav to 1. If brav is already set to 1 (default) try -1.' 606 ABI_ERROR(msg) 607 end if 608 609 ! Fourier transform of the dynamical matrices (q-->R) 610 ABI_MALLOC(ifc_tmp%atmfrc, (3,natom,3,natom,ifc_tmp%nrpt)) 611 call ftifc_q2r(ifc_tmp%atmfrc,Ifc%dynmat,gprim,natom,nqbz,ifc_tmp%nrpt,ifc_tmp%rpt,qbz, comm) 612 613 ! Eventually impose Acoustic Sum Rule on the interatomic forces 614 if (Ifc%asr > 0) call asrif9(Ifc%asr,ifc_tmp%atmfrc,natom,ifc_tmp%nrpt,ifc_tmp%rpt,ifc_tmp%wghatm) 615 616 ! The interatomic forces have been calculated 617 write(msg, '(2a)')ch10,' The interatomic forces have been obtained ' 618 call wrtout([std_out, ab_out], msg,'COLL') 619 call cwtime_report(" ifc_init1", cpu, wall, gflops) 620 621 ! Apply cutoff on ifc if needed 622 if (nsphere > 0 .or. abs(rifcsph) > tol10) then 623 call wrtout(std_out, ' Apply cutoff on IFCs.') 624 call wrtout(std_out, sjoin(" nsphere:", itoa(nsphere), ", rifcsph:", ftoa(rifcsph))) 625 call wrtout(std_out, sjoin(" Radius of biggest sphere inscribed in the WS supercell: ", ftoa(r_inscribed_sphere))) 626 call corsifc9(ddb%acell,gprim,natom,ifc_tmp%nrpt,nsphere,rifcsph,rcan,rprim,ifc_tmp%rpt,rcut_min,ifc_tmp%wghatm) 627 if (Ifc%asr > 0) then 628 call wrtout(std_out, ' Enforcing ASR on cutoffed IFCs.') 629 call asrif9(Ifc%asr,ifc_tmp%atmfrc,natom,ifc_tmp%nrpt,ifc_tmp%rpt,ifc_tmp%wghatm) 630 end if 631 end if 632 633 ! Only conserve the necessary points in rpt: in the FT algorithm the order of the points is unimportant 634 ! In the case of effective potential, we need to keep all the points 635 Ifc%nrpt = 0 636 do irpt=1,ifc_tmp%nrpt 637 if (sum(ifc_tmp%wghatm(:,:,irpt)) /= 0) Ifc%nrpt = Ifc%nrpt + 1 638 end do 639 640 ABI_CALLOC(Ifc%atmfrc,(3,natom,3,natom,Ifc%nrpt)) 641 ABI_CALLOC(Ifc%rpt,(3,Ifc%nrpt)) 642 ABI_CALLOC(Ifc%cell,(3,Ifc%nrpt)) 643 ABI_CALLOC(Ifc%wghatm,(natom,natom,Ifc%nrpt)) 644 ABI_CALLOC(Ifc%short_atmfrc,(3,natom,3,natom,Ifc%nrpt)) 645 ABI_CALLOC(Ifc%ewald_atmfrc,(3,natom,3,natom,Ifc%nrpt)) 646 647 irpt_new = 1 648 do irpt = 1, ifc_tmp%nrpt 649 if (sum(ifc_tmp%wghatm(:,:,irpt)) /= 0) then 650 Ifc%atmfrc(:,:,:,:,irpt_new) = ifc_tmp%atmfrc(:,:,:,:,irpt) 651 Ifc%rpt(:,irpt_new) = ifc_tmp%rpt(:,irpt) 652 Ifc%wghatm(:,:,irpt_new) = ifc_tmp%wghatm(:,:,irpt) 653 Ifc%cell(:,irpt_new) = ifc_tmp%cell(:,irpt) 654 Ifc%r_inscribed_sphere = r_inscribed_sphere 655 irpt_new = irpt_new + 1 656 end if 657 end do 658 659 !write(std_out,*)"nrpt before filter:", ifc_tmp%nrpt, ", after: ", ifc%nrpt 660 !do irpt=1,ifc%nrpt 661 ! write(std_out,*)ifc%rpt(:,irpt), (ifc%wghatm(ii,ii,irpt), ii=1,natom) 662 !end do 663 664 ! Copy other useful arrays. 665 Ifc%dielt = dielt 666 Ifc%nqbz = nqbz 667 668 call alloc_copy(rcan, Ifc%rcan) 669 call alloc_copy(trans, Ifc%trans) 670 call alloc_copy(dyewq0, Ifc%dyewq0) 671 call alloc_copy(qbz(:,1:nqbz), Ifc%qbz) 672 call alloc_copy(zeff, Ifc%zeff) 673 call alloc_copy(qdrp_cart, Ifc%qdrp_cart) 674 call alloc_copy(ddb%amu, Ifc%amu) 675 676 call ifc_tmp%free() 677 678 ! Compute min/max ph frequency with ab-initio q-mesh. 679 ifc%omega_minmax(1) = huge(one); ifc%omega_minmax(2) = -huge(one) 680 do iq_ibz=1,ifc%nqibz 681 if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi-parallelism 682 call ifc%fourq(crystal, ifc%qibz(:,iq_ibz), phfrq, displ_cart) 683 ifc%omega_minmax(1) = min(ifc%omega_minmax(1), minval(phfrq)) 684 ifc%omega_minmax(2) = max(ifc%omega_minmax(2), maxval(phfrq)) 685 end do 686 xval = ifc%omega_minmax(1); call xmpi_min(xval, ifc%omega_minmax(1), comm, ierr) 687 xval = ifc%omega_minmax(2); call xmpi_max(xval, ifc%omega_minmax(2), comm, ierr) 688 ! Enlarge boundaries by 30 cm-1 689 ifc%omega_minmax(1) = ifc%omega_minmax(1) - 30.0_dp/Ha_cmm1 690 ifc%omega_minmax(2) = ifc%omega_minmax(2) + 30.0_dp/Ha_cmm1 691 692 ! TODO (This is to be suppressed in a future version) 693 if (prtsrlr == 1) then 694 ! Check that the starting values are well reproduced. 695 write(msg, '(2a)' )' mkifc9 : now check that the starting values ',& 696 ' are reproduced after the use of interatomic forces ' 697 call wrtout(std_out, msg) 698 do iqpt=1,nqbz 699 qpt(:)=Ifc%qbz(:,iqpt) 700 call ifc%fourq(Crystal,qpt,phfrq,displ_cart,out_eigvec=eigvec) 701 702 ! OmegaSRLR: Perform decomposition of dynamical matrix 703 ! MG: FIXME I don't think the implementation is correct when q !=0 704 if (prtsrlr==1) then 705 call omega_decomp(ddb%amu,natom,ntypat,Crystal%typat,dynmatfull,dynmat_sr,dynmat_lr,iqpt,nqbz,eigvec) 706 end if 707 708 ! Write the phonon frequencies (this is for checking purposes). 709 ! Note: these phonon frequencies are not written on unit iout, only on unit std_out. 710 call dfpt_prtph(displ_cart,0,enunit,-1,natom,phfrq,qphnrm,qpt) 711 end do 712 end if 713 714 ! OmegaSRLR: deallocate memory used by dynmat decomposition 715 ABI_FREE(dynmatfull) 716 ABI_FREE(dynmat_sr) 717 ABI_FREE(dynmat_lr) 718 ABI_FREE(qbz) 719 720 if (nsphere == -1) call ifc_autocutoff(ifc, crystal, comm) 721 722 call cwtime_report(" ifc_init2", cpu, wall, gflops) 723 724 end subroutine ifc_init
m_ifc/ifc_outphbtrap [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_outphbtrap
FUNCTION
Print out phonon frequencies on regular grid for BoltzTrap Flag in input file is outboltztrap=1
INPUTS
ifc<ifc_type>=Stores data related to interatomic force constants. Crystal<crystal_t>=Info on the crystal structure basename = file name for output to disk ngqpt(3)=Divisions of the q-mesh nqshft=Number of shifts qshft(3,nqshft)=Shifts of the q-mesh.
OUTPUT
only write to file. This routine should be called by a single processor.
SOURCE
2532 subroutine ifc_outphbtrap(ifc, cryst, ngqpt, nqshft, qshft, basename) 2533 2534 !Arguments ------------------------------- 2535 !scalars 2536 class(ifc_type),intent(in) :: ifc 2537 integer,intent(in) :: nqshft 2538 character(len=*),intent(in) :: basename 2539 type(crystal_t),intent(in) :: cryst 2540 !arrays 2541 integer,intent(in) :: ngqpt(3) 2542 real(dp),intent(in) :: qshft(3,nqshft) 2543 2544 !Local variables ------------------------- 2545 !scalars 2546 integer,parameter :: qptopt1=1 2547 integer :: natom,imode,iq_ibz,nqbz,nqibz, nreals,unit_btrap,iatom,idir 2548 character(len=500) :: msg,format_nreals,format_line_btrap 2549 character(len=fnlen) :: outfile 2550 !arrays 2551 integer :: qptrlatt(3,3) 2552 real(dp) :: d2cart(2,3,cryst%natom,3,cryst%natom),displ(2*3*cryst%natom*3*cryst%natom) 2553 real(dp) :: phfrq(3*cryst%natom),qphon(3) 2554 real(dp),allocatable :: qbz(:,:),qibz(:,:),wtq(:) 2555 2556 ! ********************************************************************* 2557 2558 DBG_ENTER("COLL") 2559 2560 natom = cryst%natom 2561 2562 ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion 2563 qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3) 2564 call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, nqshft, qshft, nqibz, qibz, wtq, nqbz, qbz) 2565 2566 outfile = trim(basename) // '_BTRAP' 2567 write(msg, '(3a)')ch10,' Will write phonon FREQS in BoltzTrap format to file ',trim(outfile) 2568 call wrtout([std_out, ab_out], msg) 2569 2570 if (open_file(outfile,msg,newunit=unit_btrap,status="replace") /= 0) then 2571 ABI_ERROR(msg) 2572 end if 2573 2574 write (unit_btrap,'(a)') '#' 2575 write (unit_btrap,'(a)') '# ABINIT package : Boltztrap phonon file. With old BT versions remove this header before feeding to BT' 2576 write (unit_btrap,'(a)') '# for compatibility with PHON output the freq are in Ry (before the square)' 2577 write (unit_btrap,'(a)') '#' 2578 write (unit_btrap,'(a)') '# nq, nband ' 2579 write (unit_btrap,'(a)') '# qx, qy, qz ' 2580 write (unit_btrap,'(a)') '# qpt weight ' 2581 write (unit_btrap,'(a)') '# freq_1^2, dynmat column for mode 1 ' 2582 write (unit_btrap,'(a)') '# etc for mode 2,3,4... qpt 2,3,4... ' 2583 write (unit_btrap,'(2I6)') nqibz, 3*natom 2584 2585 ! Loop over irreducible q-points 2586 do iq_ibz=1,nqibz 2587 qphon(:)=qibz(:,iq_ibz) 2588 2589 call ifc%fourq(cryst, qphon, phfrq, displ, out_d2cart=d2cart) 2590 2591 write (unit_btrap,'(3E20.10)') qphon 2592 write (unit_btrap,'(E20.10)') wtq(iq_ibz) 2593 nreals=1+2*3*natom 2594 call appdig(nreals,'(',format_nreals) 2595 format_line_btrap=trim(format_nreals)//'E20.10)' 2596 do iatom = 1, natom 2597 do idir = 1, 3 2598 imode = idir + 3*(iatom-1) 2599 ! factor two for Ry output - this may change in definitive BT and abinit formats 2600 write (unit_btrap,trim(format_line_btrap))phfrq(imode)*two,d2cart(1:2,1:3,1:natom,idir,iatom) 2601 end do 2602 end do 2603 2604 end do !irred q-points 2605 close (unit_btrap) 2606 2607 ABI_FREE(qibz) 2608 ABI_FREE(qbz) 2609 ABI_FREE(wtq) 2610 2611 DBG_EXIT("COLL") 2612 2613 end subroutine ifc_outphbtrap
m_ifc/ifc_print [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_print
FUNCTION
Print info on the object
INPUTS
[unit]=Unit number for output. Defaults to std_out [prtvol]=Verbosity level. [header]=String to be printed as header for additional info.
OUTPUT
Only printing
SOURCE
832 subroutine ifc_print(ifc, header, unit, prtvol) 833 834 !Arguments ------------------------------------ 835 !scalars 836 class(ifc_type),intent(in) :: ifc 837 integer,optional,intent(in) :: unit,prtvol 838 character(len=*),optional,intent(in) :: header 839 840 !Local variables------------------------------- 841 integer :: unt,my_prtvol,iatom,ii,idir 842 character(len=500) :: msg 843 ! ********************************************************************* 844 845 unt = std_out; if (present(unit)) unt = unit 846 my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol 847 848 msg = ' ==== Info on the interatomic force constants ==== ' 849 if (present(header)) msg = ' ==== '//trim(adjustl(header))//' ==== ' 850 call wrtout(unt, msg) 851 852 call wrtout(unt,' Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):') 853 do ii=1,3 854 write(msg,'(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)')& 855 'R(',ii,')=',ifc%rprim(:,ii),'G(',ii,')=',ifc%gprim(:,ii) 856 call wrtout(unt,msg) 857 end do 858 call wrtout(unt, sjoin(" acell:", ltoa(ifc%acell))) 859 call wrtout(unt, sjoin(" Acoustic Sum Rule option (asr):", itoa(ifc%asr))) 860 call wrtout(unt, sjoin(" Option for the sampling of the BZ (brav):", itoa(ifc%brav))) 861 call wrtout(unt, sjoin(" Symmetrization flag (symdynmat):", itoa(ifc%symdynmat))) 862 call wrtout(unt, sjoin(" Dipole-dipole interaction flag (dipdip):", itoa(ifc%dipdip))) 863 call wrtout(unt, sjoin(" Dipole-quadrupole interaction flag (dipquad):", itoa(ifc%dipquad))) 864 call wrtout(unt, sjoin(" quadrupole-quadrupole interaction flag (quadquad):", itoa(ifc%quadquad))) 865 call wrtout(unt, sjoin(" Ewald option:", itoa(ifc%ewald_option))) 866 call wrtout(unt, sjoin(" Dielectric tensor: ", ch10, ltoa(reshape(ifc%dielt, [9]), fmt="f10.2"))) 867 call wrtout(unt, " Effective charges:") 868 do iatom=1,ifc%natom 869 call wrtout(unt, ltoa(reshape(ifc%zeff(:,:,iatom), [3*3]), fmt="f10.2")) 870 end do 871 call wrtout(unt, " Quadrupolar terms:") 872 do iatom=1,ifc%natom 873 do idir=1,3 874 call wrtout(unt, ltoa(reshape(ifc%qdrp_cart(:,:,idir,iatom), [3*3]), fmt="f10.2")) 875 end do 876 end do 877 878 call wrtout(unt, sjoin(" Mass of the atoms (atomic mass unit): ", ltoa(ifc%amu))) 879 call wrtout(unt, sjoin(" Number of real-space points for IFC(R): ", itoa(ifc%nrpt))) 880 call wrtout(std_out, sjoin(" Radius of biggest sphere inscribed in the WS supercell: ", ftoa(ifc%r_inscribed_sphere))) 881 call wrtout(unt, " ") 882 883 call wrtout(unt, " Q-mesh:") 884 call wrtout(unt, sjoin(" ngqpt:", ltoa(ifc%ngqpt),", nqshft:", itoa(ifc%nqshft))) 885 do ii=1,ifc%nqshft 886 call wrtout(unt, sjoin(" ", ktoa(ifc%qshft(:,ii)))) 887 end do 888 889 end subroutine ifc_print
m_ifc/ifc_printbxsf [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_printbxsf
FUNCTION
Output phonon isosurface in Xcrysden format.
INPUTS
ifc<ifc_type>=Stores data related to interatomic force constants. crystal<crystal_t>=Info on the crystal structure ngqpt(3)=Divisions of the q-mesh nqshft=Number of shifts qshft(3,nqshft)=Shifts of the q-mesh. path=File name for output to disk comm=MPI communicator.
OUTPUT
Only write to file
SOURCE
2639 subroutine ifc_printbxsf(ifc, cryst, ngqpt, nqshft, qshft, path, comm) 2640 2641 !Arguments ------------------------------- 2642 !scalars 2643 class(ifc_type),intent(in) :: ifc 2644 integer,intent(in) :: nqshft,comm 2645 character(len=*),intent(in) :: path 2646 type(crystal_t),intent(in) :: cryst 2647 !arrays 2648 integer,intent(in) :: ngqpt(3) 2649 real(dp),intent(in) :: qshft(3,nqshft) 2650 2651 !Local variables ------------------------- 2652 !scalars 2653 integer,parameter :: nsppol1=1,master=0,qptopt1=1 2654 integer :: my_rank,nprocs,iq_ibz,nqibz,nqbz,ierr 2655 character(len=500) :: msg 2656 !arrays 2657 integer :: qptrlatt(3,3),dummy_symafm(cryst%nsym) 2658 real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom) 2659 real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:),freqs_qibz(:,:) 2660 2661 ! ********************************************************************* 2662 2663 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 2664 2665 ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion 2666 qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3) 2667 call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt1, nqshft, qshft, nqibz, qibz, wtq, nqbz, qbz) 2668 ABI_FREE(qbz) 2669 ABI_FREE(wtq) 2670 2671 ! Compute phonon frequencies in the irreducible wedge. 2672 ABI_CALLOC(freqs_qibz, (3*cryst%natom, nqibz)) 2673 2674 do iq_ibz=1,nqibz 2675 if (mod(iq_ibz, nprocs) /= my_rank) cycle ! mpi parallelism. 2676 call ifc%fourq(cryst, qibz(:,iq_ibz), freqs_qibz(:,iq_ibz), displ_cart) 2677 end do 2678 call xmpi_sum(freqs_qibz, comm, ierr) 2679 2680 ! Output phonon isosurface. 2681 if (my_rank == master) then 2682 dummy_symafm = 1 2683 call printbxsf(freqs_qibz, zero, zero, cryst%gprimd, qptrlatt, 3*cryst%natom,& 2684 nqibz, qibz, cryst%nsym, .False., cryst%symrec, dummy_symafm, .True., nsppol1, qshft, nqshft, path, ierr) 2685 if (ierr /=0) then 2686 msg = "Cannot produce BXSF file with phonon isosurface, see log file for more info" 2687 ABI_WARNING(msg) 2688 call wrtout(ab_out, msg, 'COLL') 2689 end if 2690 end if 2691 2692 ABI_FREE(freqs_qibz) 2693 ABI_FREE(qibz) 2694 2695 end subroutine ifc_printbxsf
m_ifc/ifc_speedofsound [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_speedofsound
FUNCTION
Calculate the speed of sound by averaging the phonon group velocities of the three acoustic modes on a small sphere of radius qrad centered around Gamma. Perform spherical integration with Lebedev-Laikov grids
INPUTS
ifc<ifc_type>=Object containing the dynamical matrix and the IFCs. crystal<crystal_t> = Information on the crystalline structure. qrad_tolkms(2): qrad=Radius of the sphere in reciprocal space atols_kms=Absolute tolerance in kilometer/second. The code generates spherical meshes until the results are converged twice within atols_kms. ncid=the id of the open NetCDF file. Use nctk_noid to disable netcdf output. comm=MPI communicator.
OUTPUT
SOURCE
1185 subroutine ifc_speedofsound(ifc, crystal, qrad_tolkms, ncid, comm) 1186 1187 !Arguments ------------------------------- 1188 !scalars 1189 class(ifc_type),intent(in) :: ifc 1190 integer,intent(in) :: comm,ncid 1191 type(crystal_t),intent(in) :: crystal 1192 !arrays 1193 real(dp),intent(in) :: qrad_tolkms(2) 1194 1195 !Local variables ------------------------- 1196 !scalars 1197 integer,parameter :: master=0 1198 integer :: ii,nu,igrid,my_rank,nprocs,ierr,converged,npts,num_negw,vs_ierr,ncerr 1199 integer :: iatom,iatref,num_acoustic,isacoustic 1200 real(dp) :: min_negw,cpu,wall,gflops 1201 real(dp) :: qrad,tolkms,diff 1202 character(len=500) :: msg 1203 type(lebedev_t) :: lgrid 1204 !arrays 1205 integer :: asnu(3) 1206 real(dp) :: qred(3),qvers_cart(3),qvers_red(3),quad(3),prev_quad(3),vs(7,3) 1207 real(dp) :: phfrqs(3*crystal%natom),dwdq(3,3*crystal%natom) 1208 real(dp) :: displ_cart(2,3*crystal%natom,3*crystal%natom),eigvec(2,3*crystal%natom,3*crystal%natom) 1209 1210 ! ********************************************************************* 1211 1212 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1213 1214 if (ifc%asr == 0) ABI_WARNING("Computing speed of sound with asr == 0! Use asr > 0") 1215 qrad = qrad_tolkms(1); tolkms = qrad_tolkms(2) 1216 ABI_CHECK(qrad > zero, "vs_qrad <= 0") 1217 ABI_CHECK(tolkms > zero, "vs_tolkms <= 0") 1218 1219 call cwtime(cpu, wall, gflops, "start") 1220 1221 ! Find the index of the first acoustic modes (needed to handle systems with unstable branches at Gamma 1222 ! In this case, indeed, we end up computing derivatives for the wrong branches, the integration becomes 1223 ! unstable and difficult to converge. 1224 qred = zero 1225 call ifc%fourq(crystal, qred, phfrqs, displ_cart,out_eigvec=eigvec) 1226 !write(std_out,*)"omega(q==Gamma): ",phfrqs 1227 1228 num_acoustic = 0 1229 do nu = 1, 3*crystal%natom 1230 ! Check if this mode is acoustic like: scalar product of all displacement vectors are collinear 1231 isacoustic = 1 1232 ! Find reference atom with non-zero displacement 1233 do iatom=1,crystal%natom 1234 if(sum(displ_cart(:,(iatom-1)*3+1:(iatom-1)*3+3,nu)**2) >tol16)iatref=iatom 1235 end do 1236 ! Now compute scalar product, and check they are all positive 1237 do iatom = 1, crystal%natom 1238 if (sum(eigvec(:,(iatom-1)*3+1:(iatom-1)*3+3, nu)*eigvec(:,(iatref-1)*3+1:(iatref-1)*3+3, nu)) < tol16 ) isacoustic = 0 1239 end do 1240 if (isacoustic == 1) then 1241 num_acoustic=num_acoustic+1 1242 asnu(num_acoustic)=nu 1243 if (num_acoustic==3) exit 1244 end if 1245 end do 1246 1247 ABI_CHECK(num_acoustic == 3, sjoin("Wrong number of acoustic modes:", itoa(num_acoustic))) 1248 1249 write(std_out,"(a,3i2,a)") "The bands with indices ",asnu(:)," will be used to calculate the sound velocities" 1250 1251 ! Speed of sound along reduced directions. 1252 do ii=1,6 1253 qred = zero; qred(MOD(ii-1,3)+1) = one 1254 if (ii >= 4 .and. ii <= 6) qred = matmul(crystal%rprimd, qred) ! Cartesian directions. 1255 !if (ii >= 7 .and. ii <= 9) qred = matmul(crystal%rprimd, qred) ! Cartesian directions. 1256 qvers_red = (qred / normv(qred, crystal%gmet, "G")) 1257 qred = qrad * qvers_red 1258 !write(std_out,*)"dir",normv(qred, crystal%gmet, "G"), qrad 1259 call ifc%fourq(crystal, qred, phfrqs, displ_cart, dwdq=dwdq) 1260 1261 do nu=1,3 1262 vs(ii, nu) = sqrt(sum(dwdq(1:3,asnu(nu)) ** 2)) * Bohr_meter * 0.001_dp / Time_Sec 1263 end do 1264 write(std_out,"(a,3es12.4,a)")" ||vs(nu)||:",vs(ii,:), " [km/s]" 1265 1266 qvers_cart = matmul(crystal%gprimd, qvers_red) * two_pi 1267 do nu=1,3 1268 vs(ii, nu) = dot_product(dwdq(1:3,asnu(nu)), qvers_cart) * Bohr_meter * 0.001_dp / Time_Sec 1269 end do 1270 write(std_out,"(a,3es12.4,a)")" <q|vs(nu)>:",vs(ii,:), " [km/s]" 1271 1272 !do nu=1,3 1273 ! write(std_out,"(a,3es12.4,a)")" vs(nu)_vect_red:",& 1274 ! matmul(crystal%gprimd, dwdq(1:3,asnu(nu))) * Bohr_meter * 0.001_dp / Time_Sec, " [km/s]" 1275 !end do 1276 end do 1277 1278 ! Spherical average with Lebedev-Laikov grids. 1279 converged = 0 1280 do igrid=1,lebedev_ngrids 1281 call lgrid%from_index(igrid) 1282 npts = lgrid%npts; quad = zero; num_negw = 0; min_negw = zero 1283 do ii=1,npts 1284 if (mod(ii, nprocs) /= my_rank) cycle ! mpi-parallelism 1285 1286 ! Build q-point on sphere of radius qrad. qcart --> qred 1287 qred = matmul(crystal%rprimd, lgrid%versors(:, ii)) 1288 qred = qrad * (qred / normv(qred, crystal%gmet, "G")) 1289 !write(std_out,*)"lebe",normv(qred, crystal%gmet, "G"), qrad 1290 call ifc%fourq(crystal, qred, phfrqs, displ_cart, dwdq=dwdq) 1291 if (any(phfrqs(asnu) < -tol8)) then 1292 num_negw = num_negw + 1; min_negw = min(min_negw, minval(phfrqs(asnu))) 1293 end if 1294 1295 do nu=1,3 1296 quad(nu) = quad(nu) + lgrid%weights(ii) * sqrt(sum(dwdq(1:3,asnu(nu)) ** 2)) 1297 !quad(nu) = quad(nu) + lgrid%weights(ii) * abs(dot_product(lgrid%versors(:,ii), dwdq(:,asnu(nu)))) 1298 end do 1299 end do 1300 1301 ! Will use km/sec unit for echo purposes 1302 quad = quad * Bohr_meter * 0.001_dp / Time_Sec 1303 call xmpi_sum(quad, comm, ierr) 1304 call xmpi_sum(num_negw, comm, ierr) 1305 call lgrid%free() 1306 1307 write(std_out,'(2(a,i6),a,3es12.4,a,es12.4)') & 1308 " Lebedev-Laikov grid: ",igrid,", npts: ", npts, " vs_sphavg(ac_modes): ",quad, " <vs>: ",sum(quad)/3 1309 1310 if (igrid > 1) then 1311 diff = zero 1312 do nu=1,3 1313 diff = diff + abs(quad(nu) - prev_quad(nu)) / 3 1314 end do 1315 !if (abs(sum(quad - prev_quad)/3) < tolkms) then 1316 if (diff < tolkms) then 1317 converged = converged + 1 1318 else 1319 converged = 0 1320 end if 1321 end if 1322 prev_quad = quad 1323 vs(7, :) = quad 1324 if (converged == 2) exit 1325 end do ! igrid 1326 1327 if (my_rank == master) then 1328 ! vs_err: 1 if not converged, < 0 if negative freqs, == 0 if success. 1329 vs_ierr = 0 1330 do ii=1,3 1331 write(ab_out,"(a,3es12.4,a,i1)")" Speed of sound:",vs(ii,:)," [km/s] along reduced direction: ",ii 1332 end do 1333 write(ab_out,'(2(a,es12.4),a,i0)') & 1334 " Lebedev-Laikov integration with qradius: ", qrad, " tolkms: ",tolkms, " [km/s], npts: ", npts 1335 write(ab_out,"(a,3es12.4,a,es12.4)")" Spherical average:",vs(7,:)," [km/s], ",sum(vs(7,:))/3 1336 if (converged /= 2) then 1337 vs_ierr = 1 1338 write(msg,'(a,es12.4,a)')" WARNING: Results are not converged within: ",tolkms, " [km/s]" 1339 call wrtout(ab_out, msg) 1340 ABI_WARNING(msg) 1341 end if 1342 if (num_negw > 0) then 1343 vs_ierr = -num_negw 1344 write(msg,'(a,i0,a,es12.4,3a)') & 1345 " WARNING: Detected ",num_negw, " negative frequencies. Minimum was: ",min_negw * Ha_meV, "[meV]",ch10,& 1346 " Speed of sound could be wrong" 1347 call wrtout(ab_out, msg) 1348 ABI_WARNING(msg) 1349 end if 1350 1351 ! Dump results to netcdf file. 1352 if (ncid /= nctk_noid) then 1353 #ifdef HAVE_NETCDF 1354 ncerr = nctk_def_arrays(ncid, [nctkarr_t("vsound", "dp", "seven, three")], defmode=.True.) 1355 NCF_CHECK(ncerr) 1356 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "vsound_ierr"]) 1357 NCF_CHECK(ncerr) 1358 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "vsound_qrad", "vsound_tolkms"]) 1359 NCF_CHECK(ncerr) 1360 ncerr = nctk_def_arrays(ncid, [nctkarr_t("asnu", "i", "three")], defmode=.True.) 1361 NCF_CHECK(ncerr) 1362 ! Write data. 1363 NCF_CHECK(nctk_set_datamode(ncid)) 1364 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_ierr"), vs_ierr)) 1365 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_qrad"), qrad)) 1366 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound_tolkms"), tolkms)) 1367 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vsound"), vs)) 1368 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "asnu"), asnu)) 1369 #endif 1370 end if 1371 end if 1372 1373 call cwtime_report(" ifc_speedofsound", cpu, wall, gflops) 1374 1375 end subroutine ifc_speedofsound
m_ifc/ifc_type [ Types ]
NAME
ifc_type
FUNCTION
Contains the necessary data to interpolate the phonon bandstructure and eigenvectors in reciprocal space (ie. interatomic force constants and corresponding real space grid info).
SOURCE
72 type,public :: ifc_type 73 74 integer :: natom 75 ! Number of atoms in the unit cell. 76 77 integer :: mpert 78 ! Maximum number of ipert. 79 80 integer :: asr 81 ! Option for the treatment of the Acoustic Sum Rule. 82 83 integer :: brav 84 ! Option for the sampling of the BZ (anaddb input variable) 85 86 integer :: dipdip 87 ! dipole dipole interaction flag. 88 89 integer :: dipquad 90 ! dipole quadrupole interaction flag. 91 92 integer :: quadquad 93 ! dipole quadrupole interaction flag. 94 95 integer :: symdynmat 96 ! If equal to 1, the dynamical matrix is symmetrized in dfpt_phfrq before the diagonalization. 97 98 integer :: nqshft 99 ! Number of shifts in the q-mesh (usually 1 since the mesh is gamma-centered!) 100 101 integer :: nqibz 102 ! Number of points in the IBZ 103 104 integer :: nqbz 105 ! Number of points in the full BZ 106 107 integer :: nrpt 108 ! Number of real space points used to integrate IFC (for interpolation of dynamical matrices) 109 110 integer :: ngqpt(3) 111 ! Number of division in the Q mesh. 112 113 integer :: ewald_option 114 ! Option for the ewald sum 115 116 real(dp) :: rprim(3,3),gprim(3,3),acell(3) 117 ! These values are used to call anaddb routines that don't use rprimd, gprimd 118 119 real(dp) :: dielt(3,3) 120 ! Dielectric tensor (Cartesian coordinates) 121 122 real(dp) :: omega_minmax(2) 123 ! Min and max frequency obtained on the initial ab-initio q-mesh (-+ 30 cmm1) 124 ! Used to generate frequency meshes for DOSes. 125 126 real(dp) :: r_inscribed_sphere 127 ! radius of biggest sphere inscribed in the WS supercell 128 129 real(dp),allocatable :: amu(:) 130 ! amu(ntypat) 131 ! mass of the atoms (atomic mass unit) 132 133 real(dp),allocatable :: atmfrc(:,:,:,:,:) 134 ! atmfrc(3,natom,3,natom,nrpt) 135 ! Inter atomic forces in real space 136 137 integer,allocatable :: cell(:,:) 138 ! cell(nrpt,3) 139 ! Give the index of the cell and irpt 140 141 real(dp),allocatable :: ewald_atmfrc(:,:,:,:,:) 142 ! Ewald_atmfrc(3,natom,3,natom,nrpt) 143 ! Ewald Inter atomic forces in real space 144 145 real(dp),allocatable :: short_atmfrc(:,:,:,:,:) 146 ! short_atmfrc(3,natom,3,natom,nrpt) 147 ! Short range part of Inter atomic forces in real space 148 149 real(dp),allocatable :: qshft(:,:) 150 ! qshft(3,nqshft) 151 ! The shifts of the q-mesh 152 153 real(dp), allocatable :: rpt(:,:) 154 ! rpt(3,nrpt) 155 ! Real space points in canonical type coordinates. 156 157 real(dp),allocatable :: wghatm(:,:,:) 158 ! wghatm(natom,natom,nrpt) 159 ! Weights for each point and atom in the Wigner Seitz supercell in real space. 160 161 real(dp),allocatable :: rcan(:,:) 162 ! rcan(3,natom) 163 ! Atomic position in canonical coordinates. 164 165 real(dp),allocatable :: trans(:,:) 166 ! trans(3,natom) 167 ! Atomic translations: xred = rcan + trans 168 169 real(dp),allocatable :: dyewq0(:,:,:) 170 ! dyewq0(3,3,natom) 171 ! Atomic self-interaction correction to the dynamical matrix (only when dipdip = 1). 172 173 real(dp),allocatable :: zeff(:,:,:) 174 ! zeff(3,3,natom) 175 ! Effective charge on each atom, versus electric field and atomic displacement. 176 ! Cartesian coordinates 177 178 real(dp),allocatable :: qdrp_cart(:,:,:,:) 179 ! qdrp_cart(3,3,3,natom) 180 ! Quadrupole tensor on each atom 181 ! Cartesian coordinates 182 183 real(dp),allocatable :: qibz(:,:) 184 ! qibz(3,nqibz)) 185 ! List of q-points in the IBZ 186 187 real(dp),allocatable :: wtq(:) 188 ! wtq(nqibz)) 189 ! q-point Weights. 190 191 real(dp),allocatable :: qbz(:,:) 192 ! qbz(3,nqbz)) 193 ! List of q-points in the full BZ 194 195 real(dp),allocatable :: dynmat(:,:,:,:,:,:) 196 ! dynmat(2,3,natom,3,natom,nqbz)) 197 ! dynamical matrices relative to the q points of the B.Z. sampling 198 ! Note that the long-range dip-dip part has been removed if dipdip = 1 199 ! Moreover the array is multiplied by a phase shift in mkifc9. 200 201 !real(dp),allocatable :: dynmat_lr(:,:,:,:,:,:) 202 ! dynmat_lr(2,3,natom,3,natom,nqbz)) 203 ! Long-range part of dynmat in q-space 204 205 contains 206 207 procedure :: free => ifc_free 208 ! Release memory 209 210 procedure :: print => ifc_print 211 ! Print info on the object. 212 213 procedure :: fourq => ifc_fourq 214 ! Use Fourier interpolation to compute interpolated frequencies w(q) and eigenvectors e(q) 215 216 procedure :: get_dwdq => ifc_get_dwdq 217 ! Compute phonon group velocities at an arbitrary q-point. 218 219 procedure :: get_phmesh => ifc_get_phmesh 220 ! Build linear mesh for phonons. 221 222 procedure :: speedofsound => ifc_speedofsound 223 ! Compute the speed of sound by averaging phonon group velocities. 224 225 procedure :: write => ifc_write 226 ! Print the ifc (output, netcdf and text file) 227 228 procedure :: outphbtrap => ifc_outphbtrap 229 ! Print out phonon frequencies on regular grid for BoltzTrap code. 230 231 procedure :: printbxsf => ifc_printbxsf 232 ! Output phonon isosurface in Xcrysden format. 233 234 procedure :: calcnwrite_nana_terms => ifc_calcnwrite_nana_terms 235 ! Compute phonons for q--> 0 with LO-TO 236 237 procedure :: init => ifc_init 238 ! Constructor from DDB datatype 239 240 procedure :: from_file => ifc_from_file 241 ! Constructor from filename 242 end type ifc_type
m_ifc/ifc_write [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
ifc_write
FUNCTION
Adds the real-space interatomic force constants to: the output file, a NetCDF file which is already open on ncid if prt_ifc==1, to the ifcinfo.out file to a TDEP file named outfile.forceconstants_ABINIT
INPUTS
Ifc<type(ifc_type)>=Object containing the dynamical matrix and the IFCs. ifcana= 0 => no analysis of ifc ; 1 => full analysis atifc(natom) = atifc(ia) equals 1 if the analysis of ifc has to be done for atom ia; otherwise 0. ifcout= Number of interatomic force constants written in the output file prt_ifc = flag to print out ifc information for dynamical matrix (AI2PS) ncid=the id of the open NetCDF file. Set to nctk_noid if netcdf output is not wanted.
OUTPUT
written in the output file and in the NetCDF file
NOTES
This routine should be executed by one processor only TODO: 1) ifc_write should not have side effects 2) the code is unreadable and horrible - 3/4 different file formats for the same stuff. We should make different subroutines, even if it duplicates some code
SOURCE
1675 subroutine ifc_write(Ifc,ifcana,atifc,ifcout,prt_ifc,ncid, & 1676 unit_out) ! optional arguments 1677 1678 !Arguments ------------------------------- 1679 !scalars 1680 class(ifc_type),intent(inout) :: Ifc 1681 integer,intent(in) :: ifcout,ifcana,prt_ifc,ncid 1682 integer,optional,intent(in) :: unit_out 1683 !arrays 1684 integer,intent(in) :: atifc(Ifc%natom) 1685 1686 !Local variables ------------------------- 1687 !scalars 1688 integer :: ia,ib,ii,ncerr,iatifc,ifcout1,mu,nu,iout, irpt 1689 ! unit number to print out ifc information for dynamical matrix (AI2PS) 1690 integer :: unit_ifc, unit_tdep 1691 real(dp) :: detdlt 1692 real(dp) :: maxdist_tdep 1693 character(len=500) :: message 1694 character(len=4) :: str1, str2 1695 !arrays 1696 integer,allocatable :: list(:),indngb(:) 1697 real(dp) :: invdlt(3,3),ra(3),xred(3),dielt(3,3) 1698 real(dp),allocatable :: dist(:,:,:),wkdist(:),rsiaf(:,:,:),sriaf(:,:,:),vect(:,:,:) 1699 real(dp),allocatable :: posngb(:,:),wghia(:) 1700 real(dp) :: gprimd(3,3),rprimd(3,3) 1701 1702 ! ********************************************************************* 1703 1704 iout = ab_out 1705 if (present(unit_out)) then 1706 iout = unit_out 1707 end if 1708 dielt = ifc%dielt 1709 1710 ! Compute the distances between atoms 1711 ABI_MALLOC(dist,(Ifc%natom,Ifc%natom,Ifc%nrpt)) 1712 call dist9(Ifc%acell,dist,Ifc%gprim,Ifc%natom,Ifc%nrpt,Ifc%rcan,Ifc%rprim,Ifc%rpt) 1713 ! Now dist(ia,ib,irpt) contains the distance from atom ia to atom ib in unit cell irpt. 1714 1715 ABI_MALLOC(list,(Ifc%natom*Ifc%nrpt)) 1716 ABI_MALLOC(wkdist,(Ifc%natom*Ifc%nrpt)) 1717 1718 ! Calculating the inverse (transpose) of the dielectric tensor 1719 call matr3inv(dielt,invdlt) 1720 1721 ! Calculating the determinant of the dielectric tensor 1722 detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*& 1723 & dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*& 1724 & dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-& 1725 & dielt(1,2)*dielt(2,1)*dielt(3,3) 1726 1727 ! echo to log file 1728 write(std_out,'(a)' )' ifc_write: analysis of interatomic force constants ' 1729 call mkrdim(Ifc%acell,Ifc%rprim,rprimd) 1730 call matr3inv(rprimd,gprimd) 1731 1732 if (iout > 0) then 1733 write(iout, '(/,a,/)' )' Analysis of interatomic force constants ' 1734 if(Ifc%dipdip==1.and.Ifc%dipquad==0.and.Ifc%quadquad==0)then 1735 write(iout, '(a)' )' Are given : column(1-3), the total force constant' 1736 write(iout, '(a)' )' then column(4-6), the Ewald part' 1737 write(iout, '(a)' )' then column(7-9), the short-range part' 1738 write(iout, '(a)' )' Column 1, 4 and 7 are related to the displacement' 1739 write(iout, '(a)' )' of the generic atom along x, ' 1740 write(iout, '(a)' )' column 2, 5 and 8 are related to the displacement' 1741 write(iout, '(a)' )' of the generic atom along y, ' 1742 write(iout, '(a)' )' column 3, 6 and 9 are related to the displacement' 1743 write(iout, '(a)')' of the generic atom along z. ' 1744 else if(Ifc%dipquad==1.or.Ifc%quadquad==1)then 1745 write(iout, '(a)' )' Are given : column(1-3), ONLY the short-range part!!!!' 1746 write(iout, '(a)' )' column 1 is related to the displacement' 1747 write(iout, '(a)' )' of the generic atom along x, ' 1748 write(iout, '(a)' )' column 2 is related to the displacement' 1749 write(iout, '(a)' )' of the generic atom along y, ' 1750 write(iout, '(a)' )' column 3 is related to the displacement' 1751 write(iout, '(a)' )' of the generic atom along z, ' 1752 else if(Ifc%dipdip==0)then 1753 write(iout, '(a)' )' column 1 is related to the displacement' 1754 write(iout, '(a)' )' of the generic atom along x, ' 1755 write(iout, '(a)' )' column 2 is related to the displacement' 1756 write(iout, '(a)' )' of the generic atom along y, ' 1757 write(iout, '(a)' )' column 3 is related to the displacement' 1758 write(iout, '(a)' )' of the generic atom along z, ' 1759 end if 1760 end if 1761 1762 if (ifcout>Ifc%natom*Ifc%nrpt .or. ifcout == -1) then 1763 ifcout1=Ifc%natom*Ifc%nrpt 1764 write(message, '(3a,i0,a)' )& 1765 & 'The value of ifcout exceeds the number of atoms in the big box.', ch10, & 1766 & 'Output limited to ',Ifc%natom*Ifc%nrpt,' atoms.' 1767 ABI_WARNING(message) 1768 else 1769 ifcout1=ifcout 1770 end if 1771 1772 ! set up file for real space ifc output, if required 1773 if (prt_ifc == 1) then 1774 if (open_file('ifcinfo.out', message, newunit=unit_ifc, status="replace") /= 0) then 1775 ABI_ERROR(message) 1776 end if 1777 write(iout, '(a,a)' )ch10,& 1778 & ' NOTE: Open file ifcinfo.out, for the output of interatomic force constants. This is because prt_ifc==1. ' 1779 1780 if (open_file('outfile.forceconstants_ABINIT', message, newunit=unit_tdep, status="replace") /= 0) then 1781 ABI_ERROR(message) 1782 end if 1783 write(iout, '(a,a,a)' )ch10,& 1784 & ' NOTE: Open file outfile.forceconstants_ABINIT, for the output of interatomic force',& 1785 & ' constants in TDEP format. This is because prt_ifc==1. ' 1786 ! Print necessary stuff for TDEP 1787 write(unit_tdep,"(1X,I10,15X,'How many atoms per unit cell')") Ifc%natom 1788 1789 ! look at all pairs, find furthest one with weight 1 1790 ! do ia 1791 !Ifc%wghatm(ia,ib,irpt) 1792 maxdist_tdep = Ifc%r_inscribed_sphere !maxval(dist)*0.8_dp 1793 write(unit_tdep,"(1X,F20.15,5X,'Realspace cutoff (A)')") maxdist_tdep*Bohr_Ang 1794 end if 1795 1796 #ifdef HAVE_NETCDF 1797 if (ncid /= nctk_noid) then 1798 ! initialize netcdf variables 1799 ncerr = nctk_def_dims(ncid, [nctkdim_t("natifc", SUM(atifc)), nctkdim_t("number_of_r_points_big_box", Ifc%nrpt), & 1800 nctkdim_t("number_of_atoms_big_box", Ifc%natom*Ifc%nrpt), nctkdim_t("ifcout", ifcout1)], defmode=.True.) 1801 NCF_CHECK(ncerr) 1802 1803 ncerr = nctk_def_arrays(ncid, [& 1804 nctkarr_t('ifc_atoms_indices', "i", "natifc"),& 1805 nctkarr_t('ifc_neighbours_indices', "i", "ifcout, natifc"),& 1806 nctkarr_t('ifc_distances', "dp", "ifcout, natifc "),& 1807 nctkarr_t('ifc_matrix_cart_coord', "dp", "number_of_cartesian_directions,number_of_cartesian_directions, ifcout, natifc"),& 1808 nctkarr_t('ifc_atoms_cart_coord', "dp", "number_of_cartesian_directions,ifcout, natifc"),& 1809 nctkarr_t('ifc_weights', "dp", "ifcout, natifc")]) 1810 NCF_CHECK(ncerr) 1811 1812 if (Ifc%dipdip==1) then 1813 ncerr = nctk_def_arrays(ncid, [& 1814 nctkarr_t('ifc_matrix_cart_coord_short_range', "dp", & 1815 "number_of_cartesian_directions, number_of_cartesian_directions, ifcout, natifc")]) 1816 NCF_CHECK(ncerr) 1817 end if 1818 1819 if (ifcana==1) then 1820 ncerr = nctk_def_arrays(ncid, [& 1821 nctkarr_t('ifc_local_vectors', "dp", "number_of_cartesian_directions, number_of_cartesian_directions, ifcout, natifc")]) 1822 NCF_CHECK(ncerr) 1823 end if 1824 1825 NCF_CHECK(nctk_set_datamode(ncid)) 1826 end if 1827 #endif 1828 1829 ABI_MALLOC(rsiaf,(3,3,ifcout1)) 1830 ABI_MALLOC(sriaf,(3,3,ifcout1)) 1831 ABI_MALLOC(vect,(3,3,ifcout1)) 1832 ABI_MALLOC(indngb,(ifcout1)) 1833 ABI_MALLOC(posngb,(3,ifcout1)) 1834 ABI_MALLOC(wghia,(ifcout1)) 1835 1836 iatifc=0 1837 1838 ! BIG loop on all generic atoms 1839 do ia=1,Ifc%natom 1840 if(atifc(ia)==1)then 1841 1842 iatifc=iatifc+1 1843 1844 ! First transform canonical coordinates to reduced coordinates 1845 do ii=1,3 1846 xred(ii)=Ifc%gprim(1,ii)*Ifc%rcan(1,ia)+Ifc%gprim(2,ii)*Ifc%rcan(2,ia)+Ifc%gprim(3,ii)*Ifc%rcan(3,ia) 1847 end do 1848 1849 ! Then to cartesian coordinates 1850 ra(:)=xred(1)*Ifc%acell(1)*Ifc%rprim(:,1)+ xred(2)*Ifc%acell(2)*Ifc%rprim(:,2)+ xred(3)*Ifc%acell(3)*Ifc%rprim(:,3) 1851 1852 ! This sorting algorithm is slow ... 1853 wkdist(:)=reshape(dist(ia,:,:),(/Ifc%natom*Ifc%nrpt/)) 1854 do ii=1,Ifc%natom*Ifc%nrpt 1855 list(ii)=ii 1856 end do 1857 call sort_dp(Ifc%natom*Ifc%nrpt,wkdist,list,tol14) 1858 1859 if (iout > 0) then 1860 write(iout, '(a)' ) 1861 write(std_out,'(a,i4)' )' generic atom number',ia 1862 write(iout, '(a,i4)' )' generic atom number',ia 1863 write(std_out,'(a,3es16.8)' ) ' with cartesian coordinates',ra(1:3) 1864 write(iout,'(a,3es16.8)' ) ' with cartesian coordinates',ra(1:3) 1865 write(iout, '(a)' ) 1866 end if 1867 1868 call ifc_getiaf(Ifc,ifcana,ifcout1,iout,ifc%zeff,ia,ra,list,dist,invdlt,& 1869 detdlt,rsiaf,sriaf,vect,indngb,posngb) 1870 1871 if (prt_ifc == 1) then 1872 do ii=1,ifcout1 1873 if (wkdist(ii) > maxdist_tdep) exit 1874 end do 1875 ii = ii - 1 1876 write(unit_tdep,"(1X,I10,15X,'How many neighbours does atom ',I3,' have')") ii, ia 1877 1878 do ii=1,ifcout1 1879 ib = indngb(ii) 1880 irpt = (list(ii)-1)/Ifc%natom+1 1881 wghia(ii) = Ifc%wghatm(ia,ib,irpt) 1882 ! limit printing to maximum distance for tdep 1883 if (wkdist(ii) > maxdist_tdep) cycle 1884 1885 !TDEP 1886 call int2char4(ii, str1) 1887 call int2char4(ia, str2) 1888 write(unit_tdep,"(1X,I10,15X,a,a,a,a)") ib, & 1889 & 'In the unit cell, what is the index of neighbour ', & 1890 & trim(str1), " of atom ", trim(str2) 1891 ! The lattice vector needs to be in reduced coordinates? 1892 ! TODO: check whether this is correct for TDEP: might need just lattice 1893 ! vector part and not full vector, and could be in integers instead of 1894 ! cartesian vector... 1895 write (unit_tdep,'(3es28.16)') matmul(Ifc%rpt(1:3,irpt),Ifc%gprim) 1896 1897 !AI2PS 1898 write(unit_ifc,'(i6,i6)') ia,ii 1899 write(unit_ifc,'(3es28.16)') posngb(1:3,ii) 1900 do nu=1,3 1901 !TDEp 1902 ! And the actual short ranged forceconstant: TODO: check if 1903 ! a transpose is needed or a swap between the nu and the mu 1904 !write(unit_tdep,'(3f28.16)') (sriaf(nu,mu,ii)*Ha_eV/amu_emass, mu=1, 3) 1905 write(unit_tdep,'(3f28.16)') (Ifc%short_atmfrc(mu,ia,nu,ib,irpt)*Ha_eV/Bohr_Ang**2, mu=1, 3) 1906 1907 !AI2PS 1908 write(unit_ifc,'(3f28.16)')(rsiaf(nu,mu,ii),mu=1,3) 1909 end do 1910 end do 1911 1912 #ifdef HAVE_NETCDF 1913 if (ncid /= nctk_noid) then 1914 NCF_CHECK(nf90_put_var(ncid, vid("ifc_atoms_indices"), ia, start=[iatifc])) 1915 NCF_CHECK(nf90_put_var(ncid, vid("ifc_neighbours_indices"), indngb, start=[1,iatifc], count=[ifcout1,1])) 1916 NCF_CHECK(nf90_put_var(ncid, vid("ifc_distances"), wkdist(:ifcout1), start=[1,iatifc],count=[ifcout1,1])) 1917 ncerr = nf90_put_var(ncid, vid("ifc_matrix_cart_coord"), rsiaf, start=[1,1,1,iatifc], count=[3,3,ifcout1,1]) 1918 NCF_CHECK(ncerr) 1919 NCF_CHECK(nf90_put_var(ncid, vid("ifc_atoms_cart_coord"), posngb, start=[1,1,iatifc], count=[3,ifcout1,1])) 1920 NCF_CHECK(nf90_put_var(ncid, vid("ifc_weights"), wghia, start=[1,iatifc], count=[ifcout1,1])) 1921 1922 if (Ifc%dipdip==1) then 1923 ncerr = nf90_put_var(ncid, vid("ifc_matrix_cart_coord_short_range"), sriaf, & 1924 start=[1,1,1,iatifc], count=[3,3,ifcout1,1]) 1925 NCF_CHECK(ncerr) 1926 end if 1927 if (ifcana==1) then 1928 ncerr = nf90_put_var(ncid, vid("ifc_local_vectors"), vect, start=[1,1,1,iatifc], count=[3,3,ifcout1,1]) 1929 NCF_CHECK(ncerr) 1930 end if 1931 end if 1932 #endif 1933 end if 1934 end if ! End the condition on atifc 1935 end do ! End Big loop on atoms in the unit cell, and corresponding test 1936 1937 1938 ! NB for future use: in TDEP the following can also be provided. 1939 ! ! Print some auxiliary information, if it is there. Such as norm of 1940 ! ! forceconstant per shell, which shells there are and so on. 1941 ! if ( fc%npairshells .gt. 0 .and. allocated(fc%pairshell) ) then 1942 ! write(u,"(1X,I10,15X,'Number of irreducible coordination shells')") fc%npairshells 1943 ! do i=1,fc%npairshells 1944 ! write(u,"(1X,I10,1X,F16.10,1X,F16.10,15X,'number atoms in shell, radius, norm of forceconstant',I0)") fc%pairshell(i)%n,fc%pairshell(i)%rad,fc%pairshell(i)%norm,i 1945 ! enddo 1946 ! do i=1,fc%npairshells 1947 ! do j=1,fc%pairshell(i)%n 1948 ! write(u,"(1X,3(1X,F18.12),2(1X,I0))") lo_chop(matmul(p%inv_latticevectors,fc%pairshell(i)%vec(:,j)),lo_sqtol),fc%pairshell(i)%atind(j),fc%pairshell(i)%pairind(j) 1949 ! enddo 1950 ! enddo 1951 ! endif 1952 1953 if (prt_ifc == 1) then 1954 close(unit_ifc) 1955 close(unit_tdep) 1956 1957 if (open_file('infile.lotosplitting_ABINIT', message, newunit=unit_tdep, status="replace") /= 0) then 1958 ABI_ERROR(message) 1959 end if 1960 write(unit_tdep,'(3es28.16)') dielt(:,1) 1961 write(unit_tdep,'(3es28.16)') dielt(:,2) 1962 write(unit_tdep,'(3es28.16)') dielt(:,3) 1963 do ia = 1, Ifc%natom 1964 do ii = 1, 3 1965 write(unit_tdep,'(3es28.16)') ifc%zeff(:,ii,ia) 1966 end do 1967 end do 1968 close(unit_tdep) 1969 end if 1970 1971 ABI_FREE(rsiaf) 1972 ABI_FREE(sriaf) 1973 ABI_FREE(vect) 1974 ABI_FREE(indngb) 1975 ABI_FREE(posngb) 1976 ABI_FREE(dist) 1977 ABI_FREE(list) 1978 ABI_FREE(wkdist) 1979 ABI_FREE(wghia) 1980 1981 #ifdef HAVE_NETCDF 1982 contains 1983 integer function vid(vname) 1984 character(len=*),intent(in) :: vname 1985 vid = nctk_idname(ncid, vname) 1986 end function vid 1987 #endif 1988 1989 end subroutine ifc_write
m_ifc/omega_decomp [ Functions ]
[ Top ] [ m_ifc ] [ Functions ]
NAME
omega_decomp
FUNCTION
Compute and return the eigenvalues (frequencies) of the short-range and long-range part of the dynamical matrix See Europhys. Lett. 33 p.713 (1996) for details. (included by U. Aschauer and EB)
INPUTS
argin(sizein)=description
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
SOURCE
2344 subroutine omega_decomp(amu,natom,ntypat,typat,dynmatfl,dynmatsr,dynmatlr,iqpt,nqpt,eigenvec) 2345 2346 !Arguments ------------------------------- 2347 !scalars 2348 integer,intent(in) :: natom,ntypat 2349 integer,intent(in) :: iqpt,nqpt 2350 !arrays 2351 integer,intent(in) :: typat(natom) 2352 real(dp),intent(in) :: amu(ntypat) 2353 real(dp),intent(inout) :: dynmatfl(2,3,natom,3,natom,nqpt) 2354 real(dp),intent(inout) :: dynmatsr(2,3,natom,3,natom,nqpt) 2355 real(dp),intent(inout) :: dynmatlr(2,3,natom,3,natom,nqpt) 2356 real(dp),intent(in) :: eigenvec(2*3*natom*3*natom) 2357 2358 !Local variables ------------------------- 2359 !scalars 2360 integer :: i1,i2,idir1,idir2,imode,ipert1,ipert2,index1,index2 2361 real(dp),parameter :: break_symm=1.0d-12 2362 real(dp) :: fac 2363 !arrays 2364 real(dp) :: nearidentity(3,3) 2365 real(dp) :: omegafl, omegasr, omegalr 2366 real(dp) :: sumfl,sumlr,sumsr,asr 2367 ! ********************************************************************* 2368 2369 !write(ab_out,*)'' 2370 !write(std_out,*) 'SR/LR decomposition: enter for wavevector number :',iqpt 2371 2372 !apply asr (note the lr part is already asred by construction in mkifc9) 2373 do ipert1=1,natom 2374 do idir1=1,3 2375 do idir2=1,3 2376 asr=0.0d0 2377 do ipert2=1,natom 2378 asr=asr+dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt) 2379 end do 2380 dynmatfl(1,idir1,ipert1,idir2,ipert1,iqpt)=dynmatfl(1,idir1,ipert1,idir2,ipert1,iqpt)-asr 2381 dynmatsr(1,idir1,ipert1,idir2,ipert1,iqpt)=dynmatsr(1,idir1,ipert1,idir2,ipert1,iqpt)-asr 2382 end do 2383 end do 2384 end do 2385 2386 2387 !This slight breaking of the symmetry allows the 2388 !results to be more portable between machines 2389 nearidentity(:,:)=1.0 2390 nearidentity(1,1)=1.0+break_symm 2391 nearidentity(3,3)=1.0-break_symm 2392 2393 2394 !Include Mass 2395 do ipert1=1,natom 2396 do ipert2=1,natom 2397 2398 fac=1.0d0/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass 2399 2400 do idir1=1,3 2401 do idir2=1,3 2402 2403 dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)=& 2404 & dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)*& 2405 & fac*nearidentity(idir1,idir2) 2406 2407 dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)=& 2408 & dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)*& 2409 & fac*nearidentity(idir1,idir2) 2410 2411 dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)=& 2412 & dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)*& 2413 & fac*nearidentity(idir1,idir2) 2414 2415 ! This is to break slightly the translation invariance, and make 2416 ! the automatic tests more portable 2417 if(ipert1==ipert2 .and. idir1==idir2)then 2418 dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)=& 2419 & dynmatfl(1,idir1,ipert1,idir2,ipert2,iqpt)+& 2420 & break_symm*natom/amu_emass/idir1*0.01d0 2421 2422 dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)=& 2423 & dynmatsr(1,idir1,ipert1,idir2,ipert2,iqpt)+& 2424 & break_symm*natom/amu_emass/idir1*0.01d0 2425 2426 dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)=& 2427 & dynmatlr(1,idir1,ipert1,idir2,ipert2,iqpt)+& 2428 & break_symm*natom/amu_emass/idir1*0.01d0 2429 end if 2430 2431 end do 2432 end do 2433 end do 2434 end do 2435 2436 !Calculation of <eigvec|Dyn_tot,Dyn_SR,Dyn_LR|eigenvec>=omega**2 2437 2438 !write(ab_out,*)'' 2439 !write(ab_out,*)'===============================================================================' 2440 write(ab_out,*)'' 2441 write(ab_out,*) 'Long-Range/Short-Range decomposed phonon freq. (cm-1)**2' 2442 write(ab_out,*) 'at wavevector number:',iqpt 2443 write(ab_out,*)'' 2444 write(ab_out,'(a13,1x,a16,2x,a16,2x,a16)') ' Mode number.','tot**2','SR**2','LR**2' 2445 write(std_out,'(a13,1x,a16,2x,a16,2x,a16)') ' Mode number.','tot**2','SR**2','LR**2' 2446 !write(ab_out,'(a12,2x,a10,2x,a10,2x,a10,2x,a16,2x,a16,2x,a16)') 'Mode number.','tot','SR','LR','tot**2','SR**2','LR**2' 2447 !write(std_out,'(a12,2x,a10,2x,a10,2x,a10,2x,a16,2x,a16,2x,a16)') 'Mode number.','tot','SR','LR','tot**2','SR**2','LR**2' 2448 2449 do imode=1,3*natom 2450 sumfl=zero; sumlr=zero; sumsr=zero 2451 2452 do ipert1=1,natom 2453 do ipert2=1,natom 2454 do i1=1,3 2455 do i2=1,3 2456 2457 index1=i1+(ipert1-1)*3+3*natom*(imode-1) 2458 index2=i2+(ipert2-1)*3+3*natom*(imode-1) 2459 ! MG FIXME: I don't think these expressions are correct when q != 0 2460 ! We should also include the imaginary part 2461 2462 sumfl = sumfl + eigenvec(2*index1-1) * dynmatfl(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1) 2463 sumlr = sumlr + eigenvec(2*index1-1) * dynmatlr(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1) 2464 sumsr = sumsr + eigenvec(2*index1-1) * dynmatsr(1,i1,ipert1,i2,ipert2,iqpt) * eigenvec(2*index2-1) 2465 end do 2466 end do 2467 end do 2468 end do 2469 2470 sumfl = sumfl * Ha_cmm1 * Ha_cmm1 2471 sumsr = sumsr * Ha_cmm1 * Ha_cmm1 2472 sumlr = sumlr * Ha_cmm1 * Ha_cmm1 2473 2474 ! Compute omega=sqrt(omega**2) 2475 if(sumfl>=1.0d-16)then 2476 omegafl=sqrt(sumfl) 2477 else if(sumfl>=-1.0d-16)then 2478 omegafl=zero 2479 else 2480 omegafl=-sqrt(-sumfl) 2481 end if 2482 2483 if(sumsr>=1.0d-16)then 2484 omegasr=sqrt(sumsr) 2485 else if(sumsr>=-1.0d-16)then 2486 omegasr=zero 2487 else 2488 omegasr=-sqrt(-sumsr) 2489 end if 2490 2491 if(sumlr>=1.0d-16)then 2492 omegalr=sqrt(sumlr) 2493 else if(sumlr>=-1.0d-16)then 2494 omegalr=zero 2495 else 2496 omegalr=-sqrt(-sumlr) 2497 end if 2498 2499 ! Output 2500 write(ab_out,'(i4,10x,s,f16.4,2x,f16.4,2x,f16.4)') imode,sumfl,sumsr,sumlr !vz_d 2501 write(std_out,'(i4,10x,s,f16.4,2x,f16.4,2x,f16.4)') imode,sumfl,sumsr,sumlr !vz_d 2502 ! write(ab_out,'(i4,8x,f10.4,2x,f10.4,2x,f10.4,2x,s,f16.6,2x,f16.6,2x,f16.6)') imode,omegafl,omegasr,omegalr,sumfl,sumsr,sumlr 2503 ! write(std_out,'(i4,8x,f10.4,2x,f10.4,2x,f10.4,2x,s,f16.6,2x,f16.6,2x,f16.6)') imode,omegafl,omegasr,omegalr,sumfl,sumsr,sumlr 2504 end do 2505 2506 end subroutine omega_decomp