TABLE OF CONTENTS


ABINIT/initmv [ Functions ]

[ Top ] [ Functions ]

NAME

 initmv

FUNCTION

 Initialize finite difference calculation of the ddk im dfptnl_mv.f

INPUTS

  dtset <type(dataset_type)> = all input variables in this dataset
  gmet(3,3) = reciprocal space metric tensor in bohr**-2
  kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
  kneigh(30,nkpt2) = index of the neighbours of each k-point
  kg_neigh(30,nkpt2,3) = necessary to construct the vector joining a k-point
                         to its nearest neighbour in case of a single k-point,
                         a line of k-points or a plane of k-points.
                         See getshell.F90 for details
  kptindex(2,nkpt3)= index of the k-points in the reduced BZ
                     related to a k-point in the full BZ
  kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ
  mband = maximum number of bands
  mkmem = number of k points which can fit in memory
  mpi_enreg = informations about MPI parallelization
  mpw = maximum number of plane waves
  nband(nkpt*nsppol)=number of bands at each k point, for each polarization
  nkpt2 = number of k-points in the reduced BZ
  nkpt3 = number of k-points in the full BZ
  nneigh = total number of neighbours required to evaluate the finite
          difference formula
  npwarr(nkpt2)=number of planewaves at each k point
  nsppol = number of spin polarizations
  occ(mband*nkpt*nsppol) = occupation number for each band for each k

OUTPUT

 cgindex(nkpt2,nsppol) = for each k-point, cgindex tores the location of the WF in the cg array
        me = index of the current processor
        ineigh = index of a neighbour
        ikpt_loc = index of the iteration on ikpt on the current processor
        ikpt_rbz = index of a k-point in the reduced BZ
 pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat between k-points
                           (see initberry.f for more explanations)

PARENTS

      nonlinear

CHILDREN

      kpgio,xmpi_sum

SOURCE

1542 subroutine initmv(cgindex,dtset,gmet,kg,kneigh,kg_neigh,kptindex,&
1543 &  kpt3,mband,mkmem,mpi_enreg,mpw,nband,nkpt2,&
1544 &  nkpt3,nneigh,npwarr,nsppol,occ,pwind)
1545 
1546 
1547 !This section has been created automatically by the script Abilint (TD).
1548 !Do not modify the following lines by hand.
1549 #undef ABI_FUNC
1550 #define ABI_FUNC 'initmv'
1551 !End of the abilint section
1552 
1553  implicit none
1554 
1555 !Arguments ------------------------------------
1556 !scalars
1557  integer,intent(in) :: mband,mkmem,mpw,nkpt2,nkpt3,nneigh,nsppol
1558  type(MPI_type),intent(inout) :: mpi_enreg
1559  type(dataset_type),intent(in) :: dtset
1560 !arrays
1561  integer,intent(in) :: kg(3,mpw*mkmem),kneigh(30,nkpt2),kg_neigh(30,nkpt2,3)
1562  integer,intent(in) :: nband(nkpt2*nsppol),npwarr(nkpt2),kptindex(2,nkpt3)
1563  integer,intent(out) :: cgindex(nkpt2,nsppol),pwind(mpw,nneigh,mkmem)
1564  real(dp),intent(in) :: gmet(3,3),kpt3(3,nkpt3),occ(mband*nkpt2*nsppol)
1565 
1566 !Local variables-------------------------------
1567 !scalars
1568  integer :: flag,iband,icg,ierr,ikg,ikg1,ikpt,ikpt2,ikpt_loc,ikpt_rbz
1569  integer :: index,ineigh,ipw,isppol,jpw,nband_k,mband_occ,mband_occ_k,npw_k
1570  integer :: npw_k1,orig,spaceComm
1571  real(dp) :: ecut_eff,sdeg
1572  character(len=500) :: message
1573 !arrays
1574  integer :: dg(3)
1575  integer,allocatable :: kg1(:,:),kg1_k(:,:),npwar1(:),npwtot(:)
1576  real(dp) :: dk(3),dk_(3)
1577  real(dp),allocatable :: kpt1(:,:)
1578 
1579 !************************************************************************
1580 
1581 !DEBUG
1582 !write(std_out,*)' initmv : enter '
1583 !ENDDEBUG
1584 
1585  if (xmpi_paral== 1) then
1586 !  BEGIN TF_CHANGES
1587    spaceComm=mpi_enreg%comm_cell
1588 !  END TF_CHANGES
1589    mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0
1590    mpi_enreg%mkmem(:) = 0
1591  end if
1592 
1593  ecut_eff = dtset%ecut*(dtset%dilatmx)**2
1594  ABI_ALLOCATE(kg1_k,(3,mpw))
1595  ABI_ALLOCATE(kg1,(3,mkmem*mpw))
1596  ABI_ALLOCATE(kpt1,(3,nkpt2))
1597  ABI_ALLOCATE(npwar1,(nkpt2))
1598  ABI_ALLOCATE(npwtot,(nkpt2))
1599  kg1_k(:,:) = 0
1600  pwind(:,:,:) = 0
1601  cgindex(:,:) = 0
1602 
1603 !Compute the number of occupied bands.
1604 !Check that it is the same for every k-point and that
1605 !nband(ikpt) is equal to this value
1606 
1607  if (nsppol == 1) then
1608    sdeg = two
1609  else if (nsppol == 2) then
1610    sdeg = one
1611  end if
1612 
1613 !DEBUG
1614 !write(std_out,*)' list of nband '
1615 !do isppol = 1, nsppol
1616 !do ikpt = 1, nkpt2
1617 !
1618 !nband_k = nband(ikpt + (isppol - 1)*nkpt2)
1619 !write(std_out,*)' isppol, ikpt, nband_k=',isppol, ikpt, nband_k
1620 !end do
1621 !end do
1622 !ENDDEBUG
1623 
1624  index = 0
1625  do isppol = 1, nsppol
1626    do ikpt = 1, nkpt2
1627 
1628      mband_occ_k = 0
1629      nband_k = nband(ikpt + (isppol - 1)*nkpt2)
1630 
1631      do iband = 1, nband_k
1632        index = index + 1
1633        if (abs(occ(index) - sdeg) < tol8) mband_occ_k = mband_occ_k + 1
1634      end do
1635 
1636      if (nband_k /= mband_occ_k) then
1637        write(message,'(a,a,a)')&
1638 &       '  In a non-linear response calculation, nband must be equal ',ch10,&
1639 &       '  to the number of valence bands.'
1640        MSG_ERROR(message)
1641      end if
1642 
1643 !    Note that the number of bands can be different for spin up and spin down
1644      if (ikpt > 1) then
1645        if (mband_occ /= mband_occ_k) then
1646          message = 'The number of valence bands is not the same for every k-point'
1647          MSG_ERROR(message)
1648        end if
1649      else
1650        mband_occ = mband_occ_k
1651      end if
1652 
1653    end do                ! close loop over ikpt
1654  end do                ! close loop over isppol
1655 
1656 !Find the location of each wavefunction
1657 
1658  icg = 0
1659  do isppol = 1, nsppol
1660    do ikpt = 1, nkpt2
1661 
1662 
1663 !    fab: inserted the shift due to the spin...
1664 
1665      nband_k = dtset%nband(ikpt+(isppol - 1)*nkpt2)
1666      npw_k = npwarr(ikpt)
1667 
1668      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me)) cycle
1669 
1670      cgindex(ikpt,isppol) = icg
1671      icg = icg + dtset%nspinor*npw_k*nband_k
1672 
1673    end do
1674  end do
1675 
1676 
1677 !Build pwind
1678 
1679  do ineigh = 1, nneigh
1680 
1681    do ikpt = 1, nkpt2
1682      ikpt2  = kneigh(ineigh,ikpt)
1683      ikpt_rbz = kptindex(1,ikpt2)   ! index of the k-point in the reduced BZ
1684      kpt1(:,ikpt) = dtset%kptns(:,ikpt_rbz)
1685    end do
1686 
1687 !  Set up the basis sphere of plane waves at kpt1
1688    kg1(:,:) = 0
1689    call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg1,&
1690 &   kpt1,mkmem,dtset%nband,nkpt2,'PERS',mpi_enreg,mpw,&
1691 &   npwar1,npwtot,dtset%nsppol)
1692 
1693    ikg = 0 ; ikg1 = 0 ; ikpt_loc = 0
1694 
1695    if(dtset%nsppol/=1)then
1696      if(mpi_enreg%nproc/=1)then
1697        message = ' At present, non-linear response calculations for spin-polarized system cannot be done in parallel.'
1698        MSG_ERROR(message)
1699      else
1700        isppol=1
1701      end if
1702    else
1703      isppol=1
1704    end if
1705 
1706    do ikpt = 1, nkpt2
1707 
1708      nband_k = dtset%nband(ikpt+(isppol - 1)*nkpt2)
1709      ikpt2  = kneigh(ineigh,ikpt)
1710      ikpt_rbz = kptindex(1,ikpt2)   ! index of the k-point in the reduced BZ
1711 
1712      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,-1,mpi_enreg%me))  cycle
1713 
1714      ikpt_loc = ikpt_loc + 1
1715 
1716      mpi_enreg%kpt_loc2ibz_sp(mpi_enreg%me, ikpt_loc, 1) = ikpt
1717 
1718      flag = 0
1719      npw_k = npwarr(ikpt)
1720      npw_k1 = npwarr(ikpt_rbz)
1721      dk_(:) = kpt3(:,ikpt2) - dtset%kptns(:,ikpt)
1722      dk(:)  = dk_(:) - nint(dk_(:)) + real(kg_neigh(ineigh,ikpt,:),dp)
1723      dg(:)  = nint(dk(:) - dk_(:))
1724 
1725 
1726      if (kptindex(2,ikpt2) == 0) then
1727        kg1_k(:,1:npw_k1) = kg1(:,ikg1+1:ikg1+npw_k1)
1728        if (dg(1)==0.and.dg(2)==0.and.dg(3)==0) flag = 1
1729      else
1730        kg1_k(:,1:npw_k1) = -1*kg1(:,ikg1+1:ikg1+npw_k1)
1731      end if
1732 
1733      orig = 1
1734      do ipw = 1, npw_k
1735        do jpw = orig, npw_k1
1736 
1737          if ((kg(1,ikg + ipw) == kg1_k(1,jpw) - dg(1)).and. &
1738 &         (kg(2,ikg + ipw) == kg1_k(2,jpw) - dg(2)).and. &
1739 &         (kg(3,ikg + ipw) == kg1_k(3,jpw) - dg(3)))  then
1740 
1741            pwind(ipw,ineigh,ikpt_loc) = jpw
1742            if (flag == 1)  orig = jpw + 1
1743            exit
1744 
1745          end if
1746 
1747        end do
1748      end do
1749 
1750      ikg = ikg + npw_k
1751      ikg1 = ikg1 + npw_k1
1752 
1753    end do     ! close loop over k-points
1754 
1755  end do    ! close loop over ineigh
1756  mpi_enreg%mkmem(mpi_enreg%me) = mkmem
1757 
1758  call xmpi_sum(mpi_enreg%kpt_loc2ibz_sp,spaceComm,ierr)
1759  call xmpi_sum(mpi_enreg%mkmem,spaceComm,ierr)
1760 
1761 !----------------------------------------------------------------------------
1762 
1763  ABI_DEALLOCATE(kg1)
1764  ABI_DEALLOCATE(kg1_k)
1765  ABI_DEALLOCATE(kpt1)
1766  ABI_DEALLOCATE(npwar1)
1767  ABI_DEALLOCATE(npwtot)
1768 
1769 end subroutine initmv

ABINIT/m_nonlinear [ Modules ]

[ Top ] [ Modules ]

NAME

  m_nonlinear

FUNCTION

 DFT calculations of non linear response functions.

COPYRIGHT

  Copyright (C) 2002-2018 ABINIT group (MVeithen,MB,LB)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_nonlinear
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_wffile
34  use m_errors
35  use m_abicore
36  use m_xmpi
37  use m_hdr
38  use m_ebands
39  use m_xcdata
40 
41  use m_dtfil,    only : status
42  use m_time,     only : timab
43  use m_symtk,    only : symmetrize_xred, littlegroup_q
44  use m_dynmat,   only : d3sym, sytens
45  use m_ddb,      only : nlopt, DDB_VERSION, dfptnl_doutput
46  use m_ddb_hdr,  only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write
47  use m_ioarr,    only : read_rhor
48  use m_kg,       only : getcut, kpgio, getph
49  use m_fft,      only : fourdp
50  use m_kpts,     only : getkgrid
51  use m_inwffil,  only : inwffil
52  use m_spacepar, only : hartre, setsym
53  use m_pawfgr,      only : pawfgr_type,pawfgr_init, pawfgr_destroy
54  use m_pawang,      only : pawang_type, pawang_init, pawang_free
55  use m_pawrad,      only : pawrad_type
56  use m_pawtab,      only : pawtab_type,pawtab_get_lsize
57  use m_paw_an,      only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_print
58  use m_paw_ij,      only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_print
59  use m_pawfgrtab,   only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free
60  use m_pawrhoij,    only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, &
61 &                          pawrhoij_bcast, pawrhoij_nullify
62  use m_pawdij,      only : pawdij, symdij
63  use m_paw_finegrid,only : pawexpiqr
64  use m_pawxc,       only : pawxc_get_nkxc
65  use m_paw_dmft,    only : paw_dmft_type
66  use m_paw_sphharm, only : setsym_ylm
67  use m_paw_nhat,    only : nhatgrid,pawmknhat
68  use m_paw_denpot,  only : pawdenpot
69  use m_paw_init,    only : pawinit,paw_gencond
70  use m_paw_tools,   only : chkpawovlp
71  use m_mkrho,       only : mkrho
72  use m_getshell,    only : getshell
73  use m_pspini,      only : pspini
74  use m_atm2fft,     only : atm2fft
75  use m_rhotoxc,     only : rhotoxc
76  use m_mpinfo,      only : proc_distrb_cycle
77  use m_mklocl,      only : mklocl
78  use m_common,      only : setup1
79  use m_fourier_interpol, only : transgrid
80  use m_paw_occupancies,  only : initrhoij
81  use m_paw_correlations, only : pawpuxinit
82  use m_mkcore,           only : mkcore
83  use m_pead_nl_loop,     only : pead_nl_loop
84  use m_dfptnl_loop,      only : dfptnl_loop
85 
86  implicit none
87 
88  private

ABINIT/nonlinear [ Functions ]

[ Top ] [ Functions ]

NAME

 nonlinear

FUNCTION

 Primary routine for conducting DFT calculations of non linear response functions.

INPUTS

  codvsn = code version
  dtfil <type(datafiles_type)> = variables related to files
  dtset <type(dataset_type)> = all input variables for this dataset
  etotal = new total energy (no meaning at output)
  iexit= exit flag
  mpi_enreg=informations about MPI pnarallelization
  occ(mband*nkpt*nsppol) = occupation number for each band and k
  xred(3,natom) = reduced atomic coordinates

OUTPUT

  npwtot(nkpt) = total number of plane waves at each k point

SIDE EFFECTS

  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)> = variables related to pseudopotentials

NOTES

 USE OF FFT GRIDS:
 =================
 In case of PAW:
 ---------------
    Two FFT grids are used:
    - A "coarse" FFT grid (defined by ecut)
      for the application of the Hamiltonian on the plane waves basis.
      It is defined by nfft, ngfft, mgfft, ...
      Hamiltonian, wave-functions, density related to WFs (rhor here), ...
      are expressed on this grid.
    - A "fine" FFT grid (defined) by ecutdg)
      for the computation of the density inside PAW spheres.
      It is defined by nfftf, ngfftf, mgfftf, ...
      Total density, potentials, ...
      are expressed on this grid.
 In case of norm-conserving:
 ---------------------------
    - Only the usual FFT grid (defined by ecut) is used.
      It is defined by nfft, ngfft, mgfft, ...
      For compatibility reasons, (nfftf,ngfftf,mgfftf)
      are set equal to (nfft,ngfft,mgfft) in that case.

PARENTS

      driver

CHILDREN

      d3sym,dfptnl_doutput,dfptnl_loop,ebands_free,fourdp,getcut
      getkgrid,getshell,hdr_free,hdr_init,hdr_update,initmv,inwffil,kpgio
      mkcore,nlopt,pspini,read_rhor,rhotoxc,setsym,setup1,status
      ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write
      symmetrize_xred,sytens,timab,wffclose,wrtout

SOURCE

 159 subroutine nonlinear(codvsn,dtfil,dtset,etotal,iexit,mpi_enreg,npwtot,occ,&
 160 &                    pawang,pawrad,pawtab,psps,xred)
 161 
 162 
 163 !This section has been created automatically by the script Abilint (TD).
 164 !Do not modify the following lines by hand.
 165 #undef ABI_FUNC
 166 #define ABI_FUNC 'nonlinear'
 167 !End of the abilint section
 168 
 169  implicit none
 170 
 171 !Arguments ------------------------------------
 172 !scalars
 173  integer,intent(in) :: iexit
 174  logical :: non_magnetic_xc
 175  real(dp),intent(inout) :: etotal
 176  character(len=6),intent(in) :: codvsn
 177  type(MPI_type),intent(inout) :: mpi_enreg
 178  type(datafiles_type),intent(in) :: dtfil
 179  type(dataset_type),intent(inout) :: dtset
 180  type(pawang_type),intent(inout) :: pawang
 181  type(pseudopotential_type),intent(inout) :: psps
 182 !arrays
 183  integer,intent(out) :: npwtot(dtset%nkpt)
 184  real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),xred(3,dtset%natom)
 185  type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
 186  type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
 187 
 188 !Local variables-------------------------------
 189 !scalars
 190  logical :: paral_atom,call_pawinit,qeq0
 191  integer,parameter :: level=50,formeig=0,response=1,cplex1=1
 192  integer :: ask_accurate,band_index,bantot,cplex,dum_nshiftk,flag,gnt_option,gscase
 193  integer :: has_dijnd,has_diju,has_kxc,has_k3xc
 194  integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert
 195  integer :: iatom,indx,iband,ider,idir,ierr,ifft,ikpt,ipert,isppol
 196  integer :: ireadwf0,iscf_eff,ispden,itypat,izero,mcg,me,mgfftf,mkmem_max,mpert,my_natom
 197  integer :: n1,n3xccc,natom,nband_k,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim
 198  integer :: nkpt_eff,nkpt_max,nkpt3,nkxc,nkxc1,nk3xc,nk3xc1,nneigh,ntypat,nsym1,nspden_rhoij,nzlmopt
 199  integer :: optcut,optgr0,optgr1,optgr2,optrad,option,optorth
 200  integer :: optatm,optdyfr,opteltfr,optgr,optstr,optv,optn,optn2
 201  integer :: psp_gencond,pead,rdwr,rdwrpaw,spaceworld,tim_mkrho,timrev
 202  integer :: use_sym,usecprj,usexcnhat
 203  real(dp),parameter :: k0(3)=(/zero,zero,zero/)
 204  real(dp) :: boxcut,compch_fft,compch_sph,ecore,ecut_eff,ecutdg_eff,ecutf
 205  real(dp) :: eei,epaw,epawdc,enxc,etot,fermie
 206  real(dp) :: gsqcut,gsqcut_eff,gsqcutc_eff
 207  real(dp) :: rdum,residm,ucvol,vxcavg
 208  character(len=500) :: message
 209  character(len=30) :: small_msg
 210  character(len=fnlen) :: dscrpt
 211  type(pawang_type) :: pawang1
 212  type(ebands_t) :: bstruct
 213  type(hdr_type) :: hdr,hdr_den
 214  type(ddb_hdr_type) :: ddb_hdr
 215  type(wffile_type) :: wffgs,wfftgs
 216  type(wvl_data) :: wvl
 217  type(xcdata_type) :: xcdata
 218 !arrays
 219  integer :: dum_kptrlatt(3,3),dum_vacuum(3),ngfft(18),ngfftf(18),perm(6),ii,theunit
 220  integer,allocatable :: atindx(:),atindx1(:),blkflg(:,:,:,:,:,:),carflg(:,:,:,:,:,:),cgindex(:,:)
 221  integer,allocatable :: blkflg_tmp(:,:,:,:,:,:),blkflg_sav(:,:,:,:,:,:)
 222  integer,allocatable :: carflg_tmp(:,:,:,:,:,:),carflg_sav(:,:,:,:,:,:)
 223  integer,allocatable :: d3e_pert1(:),d3e_pert2(:),d3e_pert3(:)
 224  integer,allocatable :: indsym(:,:,:),indsy1(:,:,:),irrzon(:,:,:),irrzon1(:,:,:)
 225  integer,allocatable :: kg(:,:),kneigh(:,:),kg_neigh(:,:,:)
 226  integer,allocatable :: kptindex(:,:),l_size_atm(:)
 227  integer,allocatable :: npwarr(:),nattyp(:),pwind(:,:,:),rfpert(:,:,:,:,:,:)
 228  integer,allocatable :: symq(:,:,:),symrec(:,:,:),symaf1(:),symrc1(:,:,:),symrl1(:,:,:)
 229  real(dp) :: dum_gauss(0),dum_dyfrn(0),dum_dyfrv(0),dum_eltfrxc(0)
 230  real(dp) :: dum_grn(0),dum_grv(0),dum_rhog(0),dum_vg(0)
 231  real(dp) :: dum_shiftk(3,210),dummy6(6),gmet(3,3),gprimd(3,3)
 232  real(dp) :: qphon(3),rmet(3,3),rprimd(3,3),strsxc(6),tsec(2)
 233  real(dp),allocatable :: cg(:,:),d3cart(:,:,:,:,:,:,:)
 234  real(dp),allocatable :: d3etot(:,:,:,:,:,:,:),dum_kptns(:,:)
 235 ! We need all these arrays instead of one because in Fortran the maximum number of dimensions is 7...
 236  real(dp),allocatable :: d3e_1(:,:,:,:,:,:,:),d3cart_1(:,:,:,:,:,:,:)
 237  real(dp),allocatable :: d3e_2(:,:,:,:,:,:,:),d3cart_2(:,:,:,:,:,:,:)
 238  real(dp),allocatable :: d3e_3(:,:,:,:,:,:,:),d3cart_3(:,:,:,:,:,:,:)
 239  real(dp),allocatable :: d3e_4(:,:,:,:,:,:,:),d3cart_4(:,:,:,:,:,:,:)
 240  real(dp),allocatable :: d3e_5(:,:,:,:,:,:,:),d3cart_5(:,:,:,:,:,:,:)
 241  real(dp),allocatable :: d3e_6(:,:,:,:,:,:,:),d3cart_6(:,:,:,:,:,:,:)
 242  real(dp),allocatable :: d3e_7(:,:,:,:,:,:,:),d3cart_7(:,:,:,:,:,:,:)
 243  real(dp),allocatable :: d3e_8(:,:,:,:,:,:,:),d3cart_8(:,:,:,:,:,:,:)
 244  real(dp),allocatable :: d3e_9(:,:,:,:,:,:,:),d3cart_9(:,:,:,:,:,:,:)
 245  real(dp),allocatable :: dum_wtk(:),dyfrlo_indx(:,:,:),dyfrx2(:,:,:),eigen0(:)
 246  real(dp),allocatable :: grtn_indx(:,:),grxc(:,:),k3xc(:,:),kpt3(:,:),kxc(:,:)
 247  real(dp),allocatable :: mvwtk(:,:),nhat(:,:),nhatgr(:,:,:),ph1d(:,:),ph1df(:,:),phnons(:,:,:),phnons1(:,:,:)
 248  real(dp),allocatable :: rhog(:,:),rhor(:,:),rhowfg(:,:),rhowfr(:,:),tnons1(:,:)
 249  real(dp),allocatable :: vhartr(:),vpsp(:),vtrial(:,:),vxc(:,:),work(:),xccc3d(:)
 250  type(pawfgr_type) :: pawfgr
 251  type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:)
 252  type(pawfgrtab_type),allocatable,save :: pawfgrtab(:)
 253  type(paw_an_type),allocatable :: paw_an(:)
 254  type(paw_ij_type),allocatable :: paw_ij(:)
 255  type(paw_dmft_type) :: paw_dmft
 256 
 257 ! ***********************************************************************
 258 
 259  DBG_ENTER("COLL")
 260 
 261  call timab(501,1,tsec)
 262  call status(0,dtfil%filstat,iexit,level,'enter         ')
 263 
 264 !Structured debugging if dtset%prtvol==-level
 265  if(dtset%prtvol==-level)then
 266    write(message,'(80a,a,a)')  ('=',ii=1,80),ch10,' nonlinear : enter , debug mode '
 267    call wrtout(std_out,message,'COLL')
 268  end if
 269 
 270 ! Initialise non_magnetic_xc for rhohxc
 271  non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14)
 272 
 273 !Check if the perturbations asked in the input file can be computed
 274 
 275  if (((dtset%d3e_pert1_phon == 1).and.(dtset%d3e_pert2_phon == 1)).or. &
 276 & ((dtset%d3e_pert1_phon == 1).and.(dtset%d3e_pert3_phon == 1)).or. &
 277 & ((dtset%d3e_pert2_phon == 1).and.(dtset%d3e_pert3_phon == 1))) then
 278    write(message,'(7a)')&
 279 &   'You have asked for a third-order derivative with respect to',ch10,&
 280 &   '2 or more atomic displacements.',ch10,&
 281 &   'This is not allowed yet.',ch10,&
 282 &   'Action : change d3e_pert1_phon, d3e_pert2_phon or d3e_pert3_phon in your input file.'
 283    MSG_ERROR(message)
 284  end if
 285 
 286 !Computation of third order derivatives from PEAD (pead=1) or full DPFT formalism (pead=0):
 287  pead = dtset%usepead
 288  if (pead==0) then
 289    write(message, '(2a)' ) ch10,'NONLINEAR : PEAD=0, full DFPT computation of third order derivatives'
 290    call wrtout(ab_out,message,'COLL')
 291    call wrtout(std_out,message,'COLL')
 292  end if
 293 
 294 !Some data for parallelism
 295  nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1
 296  my_natom=mpi_enreg%my_natom
 297  paral_atom=(my_natom/=dtset%natom)
 298  if (paral_atom) then
 299    MSG_BUG(" Nonlinear routine is not available yet with parallelization over atoms...")
 300  end if
 301 
 302 !Init spaceworld
 303  spaceworld=mpi_enreg%comm_cell
 304  me = xmpi_comm_rank(spaceworld)
 305 
 306 !Define FFT grid(s) sizes (be careful !)
 307 !See NOTES in the comments at the beginning of this file.
 308  call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
 309 
 310  ntypat=psps%ntypat
 311  natom=dtset%natom
 312  nfftot=product(ngfft(1:3))
 313  nfftotf=product(ngfftf(1:3))
 314 
 315 !Define the set of admitted perturbations taking into account
 316 !the possible permutations
 317  mpert=natom+6
 318  ABI_ALLOCATE(blkflg,(3,mpert,3,mpert,3,mpert))
 319  ABI_ALLOCATE(carflg,(3,mpert,3,mpert,3,mpert))
 320  ABI_ALLOCATE(rfpert,(3,mpert,3,mpert,3,mpert))
 321  ABI_ALLOCATE(d3e_pert1,(mpert))
 322  ABI_ALLOCATE(d3e_pert2,(mpert))
 323  ABI_ALLOCATE(d3e_pert3,(mpert))
 324  ABI_ALLOCATE(d3etot,(2,3,mpert,3,mpert,3,mpert))
 325  ABI_ALLOCATE(d3cart,(2,3,mpert,3,mpert,3,mpert))
 326  ABI_ALLOCATE(blkflg_tmp,(3,mpert,3,mpert,3,mpert))
 327  ABI_ALLOCATE(blkflg_sav,(3,mpert,3,mpert,3,mpert))
 328  ABI_ALLOCATE(carflg_tmp,(3,mpert,3,mpert,3,mpert))
 329  ABI_ALLOCATE(carflg_sav,(3,mpert,3,mpert,3,mpert))
 330  ABI_ALLOCATE(d3e_1,(2,3,mpert,3,mpert,3,mpert))
 331  ABI_ALLOCATE(d3e_2,(2,3,mpert,3,mpert,3,mpert))
 332  ABI_ALLOCATE(d3e_3,(2,3,mpert,3,mpert,3,mpert))
 333  ABI_ALLOCATE(d3e_4,(2,3,mpert,3,mpert,3,mpert))
 334  ABI_ALLOCATE(d3e_5,(2,3,mpert,3,mpert,3,mpert))
 335  ABI_ALLOCATE(d3e_6,(2,3,mpert,3,mpert,3,mpert))
 336  ABI_ALLOCATE(d3e_7,(2,3,mpert,3,mpert,3,mpert))
 337  ABI_ALLOCATE(d3e_8,(2,3,mpert,3,mpert,3,mpert))
 338  ABI_ALLOCATE(d3e_9,(2,3,mpert,3,mpert,3,mpert))
 339  d3e_1(:,:,:,:,:,:,:) = 0_dp
 340  d3e_2(:,:,:,:,:,:,:) = 0_dp
 341  d3e_3(:,:,:,:,:,:,:) = 0_dp
 342  d3e_4(:,:,:,:,:,:,:) = 0_dp
 343  d3e_5(:,:,:,:,:,:,:) = 0_dp
 344  d3e_6(:,:,:,:,:,:,:) = 0_dp
 345  d3e_7(:,:,:,:,:,:,:) = 0_dp
 346  d3e_8(:,:,:,:,:,:,:) = 0_dp
 347  d3e_9(:,:,:,:,:,:,:) = 0_dp
 348  ABI_ALLOCATE(d3cart_1,(2,3,mpert,3,mpert,3,mpert))
 349  ABI_ALLOCATE(d3cart_2,(2,3,mpert,3,mpert,3,mpert))
 350  ABI_ALLOCATE(d3cart_3,(2,3,mpert,3,mpert,3,mpert))
 351  ABI_ALLOCATE(d3cart_4,(2,3,mpert,3,mpert,3,mpert))
 352  ABI_ALLOCATE(d3cart_5,(2,3,mpert,3,mpert,3,mpert))
 353  ABI_ALLOCATE(d3cart_6,(2,3,mpert,3,mpert,3,mpert))
 354  ABI_ALLOCATE(d3cart_7,(2,3,mpert,3,mpert,3,mpert))
 355  ABI_ALLOCATE(d3cart_8,(2,3,mpert,3,mpert,3,mpert))
 356  ABI_ALLOCATE(d3cart_9,(2,3,mpert,3,mpert,3,mpert))
 357  blkflg(:,:,:,:,:,:) = 0
 358  d3etot(:,:,:,:,:,:,:) = 0_dp
 359  rfpert(:,:,:,:,:,:) = 0
 360  d3e_pert1(:) = 0 ; d3e_pert2(:) = 0 ; d3e_pert3(:) = 0
 361 
 362  if (dtset%d3e_pert1_phon==1) d3e_pert1(dtset%d3e_pert1_atpol(1):dtset%d3e_pert1_atpol(2))=1
 363  if (dtset%d3e_pert2_phon==1) d3e_pert2(dtset%d3e_pert2_atpol(1):dtset%d3e_pert2_atpol(2))=1
 364  if (dtset%d3e_pert3_phon==1) d3e_pert3(dtset%d3e_pert3_atpol(1):dtset%d3e_pert3_atpol(2))=1
 365  if (dtset%d3e_pert1_elfd/=0) d3e_pert1(natom+2)=1
 366  if (dtset%d3e_pert2_elfd/=0) d3e_pert2(natom+2)=1
 367  if (dtset%d3e_pert3_elfd/=0) d3e_pert3(natom+2)=1
 368 
 369  do i1pert = 1, mpert
 370    do i1dir = 1, 3
 371      do i2pert = 1, mpert
 372        do i2dir = 1, 3
 373          do i3pert = 1, mpert
 374            do i3dir = 1, 3
 375              perm(1) = &
 376 &              d3e_pert1(i1pert)*dtset%d3e_pert1_dir(i1dir) &
 377 &             *d3e_pert2(i2pert)*dtset%d3e_pert2_dir(i2dir) &
 378 &             *d3e_pert3(i3pert)*dtset%d3e_pert3_dir(i3dir)
 379              perm(2) = &
 380 &              d3e_pert1(i1pert)*dtset%d3e_pert1_dir(i1dir) &
 381 &             *d3e_pert2(i3pert)*dtset%d3e_pert2_dir(i3dir) &
 382 &             *d3e_pert3(i2pert)*dtset%d3e_pert3_dir(i2dir)
 383              perm(3) = &
 384 &              d3e_pert1(i2pert)*dtset%d3e_pert1_dir(i2dir) &
 385 &             *d3e_pert2(i1pert)*dtset%d3e_pert2_dir(i1dir) &
 386 &             *d3e_pert3(i3pert)*dtset%d3e_pert3_dir(i3dir)
 387              perm(4) = &
 388 &              d3e_pert1(i2pert)*dtset%d3e_pert1_dir(i2dir) &
 389 &             *d3e_pert2(i3pert)*dtset%d3e_pert2_dir(i3dir) &
 390 &             *d3e_pert3(i1pert)*dtset%d3e_pert3_dir(i1dir)
 391              perm(5) = &
 392 &              d3e_pert1(i3pert)*dtset%d3e_pert1_dir(i3dir) &
 393 &             *d3e_pert2(i2pert)*dtset%d3e_pert2_dir(i2dir) &
 394 &             *d3e_pert3(i1pert)*dtset%d3e_pert3_dir(i1dir)
 395              perm(6) = &
 396 &              d3e_pert1(i3pert)*dtset%d3e_pert1_dir(i3dir) &
 397 &             *d3e_pert2(i1pert)*dtset%d3e_pert2_dir(i1dir) &
 398 &             *d3e_pert3(i2pert)*dtset%d3e_pert3_dir(i2dir)
 399              if (sum(perm(:)) > 0) rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
 400            end do
 401          end do
 402        end do
 403      end do
 404    end do
 405  end do
 406 
 407 ! call timab(134,2,tsec)
 408 ! call timab(135,1,tsec)
 409 
 410 !Do symmetry stuff
 411  ABI_ALLOCATE(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 412  ABI_ALLOCATE(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 413  ABI_ALLOCATE(indsym,(4,dtset%nsym,natom))
 414  ABI_ALLOCATE(symrec,(3,3,dtset%nsym))
 415  irrzon=0;indsym=0;symrec=0;phnons=zero
 416 !If the density is to be computed by mkrho, need irrzon and phnons
 417  iscf_eff=0;if(dtset%getden==0)iscf_eff=1
 418  call setsym(indsym,irrzon,iscf_eff,natom,&
 419 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
 420 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)
 421 
 422 !Symmetrize atomic coordinates over space group elements:
 423  call symmetrize_xred(indsym,natom,dtset%nsym,dtset%symrel,dtset%tnons,xred)
 424 
 425  call status(0,dtfil%filstat,iexit,level,'call sytens   ')
 426  call sytens(indsym,mpert,natom,dtset%nsym,rfpert,symrec,dtset%symrel)
 427 
 428  write(message, '(a,a,a,a,a)' ) ch10, &
 429 & ' The list of irreducible elements of the Raman and non-linear',&
 430 & ch10,' optical susceptibility tensors is:',ch10
 431  call wrtout(ab_out,message,'COLL')
 432  call wrtout(std_out,message,'COLL')
 433 
 434  write(message,'(12x,a)')&
 435 & 'i1pert  i1dir   i2pert  i2dir   i3pert  i3dir'
 436  call wrtout(ab_out,message,'COLL')
 437  call wrtout(std_out,message,'COLL')
 438  n1 = 0
 439  do i1pert = 1, natom + 2
 440    do i1dir = 1, 3
 441      do i2pert = 1, natom + 2
 442        do i2dir = 1,3
 443          do i3pert = 1, natom + 2
 444            do i3dir = 1, 3
 445              if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
 446                n1 = n1 + 1
 447                write(message,'(2x,i4,a,6(5x,i3))') n1,')', &
 448 &               i1pert,i1dir,i2pert,i2dir,i3pert,i3dir
 449                call wrtout(ab_out,message,'COLL')
 450                call wrtout(std_out,message,'COLL')
 451              else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-2) then
 452                blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
 453                if (dtset%nonlinear_info>0) then
 454 !                 n1 = n1 + 1
 455                  write(message,'(2x,i4,a,6(5x,i3),a)') n1,')', &
 456   &               i1pert,i1dir,i2pert,i2dir,i3pert,i3dir,' => must be zero, not computed'
 457                  call wrtout(ab_out,message,'COLL')
 458                  call wrtout(std_out,message,'COLL')
 459                end if
 460              else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-1) then
 461                if (dtset%nonlinear_info>0) then
 462 !                 n1 = n1 + 1
 463                  write(message,'(2x,i4,a,6(5x,i3),a)') n1,')', &
 464   &               i1pert,i1dir,i2pert,i2dir,i3pert,i3dir,' => symmetric of an other element, not computed'
 465                  call wrtout(ab_out,message,'COLL')
 466                  call wrtout(std_out,message,'COLL')
 467                end if
 468              end if
 469            end do
 470          end do
 471        end do
 472      end do
 473    end do
 474  end do
 475  write(message,'(a,a)') ch10,ch10
 476  call wrtout(ab_out,message,'COLL')
 477  call wrtout(std_out,message,'COLL')
 478 
 479 ! For abipy :
 480  if (dtset%paral_rf == -1) then
 481    write(std_out,'(a)')"--- !IrredPerts"
 482    write(std_out,'(a)')'# List of irreducible perturbations for nonlinear'
 483    write(std_out,'(a)')'irred_perts:'
 484 
 485    n1 = 0
 486    do i1pert = 1, natom + 2
 487      do i1dir = 1, 3
 488        do i2pert = 1, natom + 2
 489          do i2dir = 1, 3
 490            do i3pert = 1, natom + 2
 491              do i3dir = 1,3
 492                if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
 493                  n1 = n1 + 1
 494                  write(std_out,'(a,i0)')"   - i1pert: ",i1pert
 495                  write(std_out,'(a,i0)')"     i1dir: ",i1dir
 496                  write(std_out,'(a,i0)')"     i2pert: ",i2pert
 497                  write(std_out,'(a,i0)')"     i2dir: ",i2dir
 498                  write(std_out,'(a,i0)')"     i3pert: ",i3pert
 499                  write(std_out,'(a,i0)')"     i3dir: ",i3dir
 500                end if
 501              end do
 502            end do
 503          end do
 504        end do
 505      end do
 506    end do
 507    write(std_out,'(a)')"..."
 508    MSG_ERROR_NODUMP("aborting now")
 509  end if
 510 
 511  call status(0,dtfil%filstat,iexit,level,'call setup1   ')
 512 
 513 !Set up for iterations
 514  call setup1(dtset%acell_orig(1:3,1),bantot,dtset,&
 515 & ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,&
 516 & natom,ngfftf,ngfft,dtset%nkpt,dtset%nsppol,&
 517 & response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw)
 518 
 519 !Set up the basis sphere of planewaves
 520  ABI_ALLOCATE(kg,(3,dtset%mpw*dtset%mkmem))
 521  ABI_ALLOCATE(npwarr,(dtset%nkpt))
 522  call status(0,dtfil%filstat,iexit,level,'call kpgio(1) ')
 523  call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg,&
 524 & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,&
 525 & dtset%nsppol)
 526 
 527 !Recompute first large sphere cut-off gsqcut, without taking into account dilatmx
 528  ecutf=dtset%ecut
 529  if (psps%usepaw==1) then
 530    ecutf=dtset%pawecutdg
 531    call wrtout(std_out,ch10//' FFT (fine) grid used in SCF cycle:','COLL')
 532  end if
 533  call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf)
 534 
 535 !Open and read pseudopotential files
 536  ecore = 0_dp
 537  call status(0,dtfil%filstat,iexit,level,'call pspini(1)')
 538  call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,&
 539 & psps,rprimd,comm_mpi=mpi_enreg%comm_cell)
 540 
 541 !Initialize band structure datatype
 542  bstruct = ebands_from_dtset(dtset, npwarr)
 543 
 544 !Initialize PAW atomic occupancies
 545  if (psps%usepaw==1) then
 546    ABI_DATATYPE_ALLOCATE(pawrhoij,(my_natom))
 547    call pawrhoij_nullify(pawrhoij)
 548    call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu, &
 549 &   my_natom,natom,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,&
 550 &   pawrhoij,dtset%pawspnorb,pawtab,dtset%spinat,dtset%typat,&
 551 &   comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
 552  else
 553    ABI_DATATYPE_ALLOCATE(pawrhoij,(0))
 554  end if
 555 
 556 !Initialize header
 557  gscase=0
 558  call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr, &
 559 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
 560 
 561 !Update header, with evolving variables, when available
 562 !Here, rprimd, xred and occ are available
 563  etot=hdr%etot ; fermie=hdr%fermie ; residm=hdr%residm
 564 !If parallelism over atom, hdr is distributed
 565  call hdr_update(hdr,bantot,etot,fermie,&
 566 & residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), &
 567 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
 568 
 569 !Clean band structure datatype (should use it more in the future !)
 570  call ebands_free(bstruct)
 571 
 572  call status(0,dtfil%filstat,iexit,level,'call inwffil(1')
 573 
 574 !Initialize wavefunction files and wavefunctions.
 575  ireadwf0=1
 576 
 577  mcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
 578  ABI_STAT_ALLOCATE(cg,(2,mcg), ierr)
 579  ABI_CHECK(ierr==0, "out-of-memory in cg")
 580 
 581  ABI_ALLOCATE(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol))
 582  eigen0(:)=zero ; ask_accurate=1
 583  optorth=0
 584 
 585  call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,&
 586 & formeig,hdr,ireadwf0,dtset%istwfk,kg,dtset%kptns,&
 587 & dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,&
 588 & dtset%nband,ngfft,dtset%nkpt,npwarr,dtset%nsppol,dtset%nsym,&
 589 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
 590 & dtfil%unkg,wffgs,wfftgs,dtfil%unwffgs,dtfil%fnamewffk,wvl)
 591 
 592 !Close wffgs, if it was ever opened (in inwffil)
 593  if (ireadwf0==1) then
 594    call WffClose(wffgs,ierr)
 595  end if
 596 
 597  if (psps%usepaw==1.and.ireadwf0==1) then
 598 !  if parallelism, pawrhoij is distributed, hdr%pawrhoij is not
 599    call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,&
 600 &   mpi_atmtab=mpi_enreg%my_atmtab)
 601  end if
 602 
 603 ! call timab(135,2,tsec)
 604 ! call timab(136,1,tsec)
 605 
 606 !Report on eigen0 values   ! Should use prteigrs.F90
 607  write(message, '(a,a)' )
 608  call wrtout(std_out,ch10//' respfn : eigen0 array','COLL')
 609  nkpt_eff=dtset%nkpt
 610  if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. dtset%nkpt>nkpt_max ) nkpt_eff=nkpt_max
 611  band_index=0
 612  do isppol=1,dtset%nsppol
 613    do ikpt=1,dtset%nkpt
 614      nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
 615      if(ikpt<=nkpt_eff)then
 616        write(message, '(a,i2,a,i5)' )'  isppol=',isppol,', k point number',ikpt
 617        call wrtout(std_out,message,'COLL')
 618        do iband=1,nband_k,4
 619          write(message, '(a,4es16.6)')'  ',eigen0(iband+band_index:min(iband+3,nband_k)+band_index)
 620          call wrtout(std_out,message,'COLL')
 621        end do
 622      else if(ikpt==nkpt_eff+1)then
 623        write(message,'(a,a)' )'  respfn : prtvol=0, 1 or 2, stop printing eigen0.',ch10
 624        call wrtout(std_out,message,'COLL')
 625      end if
 626      band_index=band_index+nband_k
 627    end do
 628  end do
 629 
 630 !Allocation for forces and atomic positions (should be taken away, also argument ... )
 631  ABI_ALLOCATE(grxc,(3,natom))
 632 
 633 !Examine the symmetries of the q wavevector
 634  ABI_ALLOCATE(symq,(4,2,dtset%nsym))
 635  timrev=1
 636 
 637 ! By default use symmetries.
 638  use_sym = 1
 639  if (dtset%prtgkk == 1)then
 640    use_sym = 0
 641    call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol,use_sym=use_sym)
 642  else
 643    call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol)
 644  end if
 645 
 646 
 647 
 648 !Generate an index table of atoms, in order for them to be used
 649 !type after type.
 650  ABI_ALLOCATE(atindx,(natom))
 651  ABI_ALLOCATE(atindx1,(natom))
 652  ABI_ALLOCATE(nattyp,(ntypat))
 653  indx=1
 654  do itypat=1,ntypat
 655    nattyp(itypat)=0
 656    do iatom=1,natom
 657      if(dtset%typat(iatom)==itypat)then
 658        atindx(iatom)=indx
 659        atindx1(indx)=iatom
 660        indx=indx+1
 661        nattyp(itypat)=nattyp(itypat)+1
 662      end if
 663    end do
 664  end do
 665 
 666 !Compute structure factor phases for current atomic pos:
 667  ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*natom))
 668  ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*natom))
 669  call getph(atindx,natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred)
 670 
 671  if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
 672    call getph(atindx,natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
 673  else
 674    ph1df(:,:)=ph1d(:,:)
 675  end if
 676 
 677 qeq0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2<1.d-14)
 678 if (.not.qeq0) then
 679   MSG_BUG('NONLINEAR with dtset%qptn!=0 is not implemented yet')
 680 end if
 681 
 682 !PAW: 1- Initialize values for several arrays depending only on atomic data
 683 !2- Check overlap
 684 !3- Identify FFT points in spheres and compute g_l(r).Y_lm(r) (and exp(-i.q.r) if needed)
 685 !4- Allocate PAW specific arrays
 686 !5- Compute perturbed local potential inside spheres
 687 !6- Eventually open temporary storage files
 688  if(psps%usepaw==1) then
 689 !  1-Initialize values for several arrays depending only on atomic data
 690 
 691    gnt_option=2
 692 
 693    ! Test if we have to call pawinit
 694    call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
 695 
 696    if (psp_gencond==1.or.call_pawinit) then
 697 !    Some gen-cond have to be added...
 698      call timab(553,1,tsec)
 699      call pawinit(gnt_option,zero,zero,dtset%pawlcutd,dtset%pawlmix,&
 700 &     psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,&
 701 &     pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%xclevel,dtset%usepotzero)
 702      call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,&
 703 &     rprimd,symrec,pawang%zarot)
 704 
 705      ! Update internal values
 706      call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
 707 
 708      call timab(553,2,tsec)
 709    else
 710      if (pawtab(1)%has_kij  ==1) pawtab(1:psps%ntypat)%has_kij  =2
 711      if (pawtab(1)%has_nabla==1) pawtab(1:psps%ntypat)%has_nabla=2
 712    end if
 713    psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore)
 714    call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot)
 715    pawtab(:)%usepawu=0
 716    pawtab(:)%useexexch=0
 717    pawtab(:)%exchmix=zero
 718    call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
 719 &   dtset%jpawu,dtset%lexexch,dtset%lpawu,ntypat,pawang,dtset%pawprtvol,pawrad,&
 720 &   pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu)
 721    compch_fft=-1.d5;compch_sph=-1.d5
 722    usexcnhat=maxval(pawtab(:)%usexcnhat)
 723 
 724 !  Note: many derivatives of cprj are needed and used a few times only, so for simplicity the
 725 !  computation of all needed derivatives will be done on-the-fly.
 726    usecprj=0
 727 
 728 !  2-Check overlap
 729    call chkpawovlp(natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred)
 730 !  3-Identify FFT points in spheres and compute g_l(r).Y_lm(r) and exp(-i.q.r)
 731    ABI_DATATYPE_ALLOCATE(pawfgrtab,(my_natom))
 732    if (my_natom>0) then
 733      call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,&
 734 &     mpi_atmtab=mpi_enreg%my_atmtab)
 735      call pawfgrtab_init(pawfgrtab,1,l_size_atm,pawrhoij(1)%nspden,dtset%typat,&
 736 &     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 737      ABI_DEALLOCATE(l_size_atm)
 738    end if
 739    optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm
 740    optgr1=dtset%pawstgylm
 741    optgr2=dtset%pawstgylm
 742    call status(0,dtfil%filstat,iexit,level,'call nhatgrid ')
 743    call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfftf,psps%ntypat,&
 744 &   optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,&
 745 &   comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab )
 746    ABI_DATATYPE_ALLOCATE(paw_an,(my_natom))
 747    ABI_DATATYPE_ALLOCATE(paw_ij,(my_natom))
 748    call paw_an_nullify(paw_an)
 749    call paw_ij_nullify(paw_ij)
 750    has_kxc=0;nkxc1=0;cplex=1
 751    has_dijnd=0;if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1
 752    has_diju=0;if(dtset%usepawu==5.or.dtset%usepawu==6) has_diju=1
 753    has_kxc=1;nkxc1=2*dtset%nspden-1 ! LDA only
 754    call pawxc_get_nkxc(nkxc1,dtset%nspden,dtset%xclevel)
 755    has_k3xc=1; nk3xc1=3*min(dtset%nspden,2)-2 ! LDA only
 756    call paw_an_init(paw_an,dtset%natom,dtset%ntypat,nkxc1,nk3xc1,dtset%nspden,cplex,dtset%pawxcdev,&
 757 &   dtset%typat,pawang,pawtab,has_vxc=1,has_vxc_ex=1,has_kxc=has_kxc,has_k3xc=has_k3xc,&
 758 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 759    call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%pawspnorb,&
 760 &   natom,dtset%ntypat,dtset%typat,pawtab,has_dij=1,has_dijhartree=1,has_dijnd=has_dijnd,&
 761 &   has_dijso=1,has_dijU=has_diju,has_pawu_occ=1,has_exexch_pot=1,nucdipmom=dtset%nucdipmom,&
 762 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 763  else ! PAW vs NCPP
 764    usexcnhat=0;usecprj=0
 765    ABI_DATATYPE_ALLOCATE(paw_an,(0))
 766    ABI_DATATYPE_ALLOCATE(paw_ij,(0))
 767    ABI_DATATYPE_ALLOCATE(pawfgrtab,(0))
 768  end if
 769 
 770  ABI_ALLOCATE(rhog,(2,nfftf))
 771  ABI_ALLOCATE(rhor,(nfftf,dtset%nspden))
 772 
 773 !Read ground-state charge density from diskfile in case getden /= 0
 774 !or compute it from wfs that were read previously : rhor as well as rhog
 775 
 776  if (dtset%getden /= 0 .or. dtset%irdden /= 0) then
 777    ! Read rho1(r) from a disk file and broadcast data.
 778    ! This part is not compatible with MPI-FFT (note single_proc=.True. below)
 779 
 780    rdwr=1;rdwrpaw=psps%usepaw;if(ireadwf0/=0) rdwrpaw=0
 781    if (rdwrpaw/=0) then
 782      ABI_DATATYPE_ALLOCATE(pawrhoij_read,(natom))
 783      call pawrhoij_nullify(pawrhoij_read)
 784      nspden_rhoij=dtset%nspden;if (dtset%pawspnorb>0.and.dtset%nspinor==2) nspden_rhoij=4
 785      call pawrhoij_alloc(pawrhoij_read,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,&
 786 &     dtset%nsppol,dtset%typat,pawtab=pawtab)
 787    else
 788      ABI_DATATYPE_ALLOCATE(pawrhoij_read,(0))
 789    end if
 790 
 791 !    MT july 2013: Should we read rhoij from the density file ?
 792    call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor, &
 793    hdr_den, pawrhoij_read, spaceworld, check_hdr=hdr)
 794    call hdr_free(hdr_den)
 795 
 796    if (rdwrpaw/=0) then
 797      call pawrhoij_bcast(pawrhoij_read,hdr%pawrhoij,0,spaceworld)
 798      call pawrhoij_free(pawrhoij_read)
 799    end if
 800    ABI_DATATYPE_DEALLOCATE(pawrhoij_read)
 801 
 802 !  Compute up+down rho(G) by fft
 803    ABI_ALLOCATE(work,(nfftf))
 804    work(:)=rhor(:,1)
 805    call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
 806    ABI_DEALLOCATE(work)
 807 
 808  else
 809    izero=0
 810 !  Obtain the charge density from read wfs
 811 !  Be careful: in PAW, compensation density has to be added !
 812    tim_mkrho=4
 813    paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented
 814    paw_dmft%use_dmft=0 ! respfn with dmft not implemented
 815    if (psps%usepaw==1) then
 816      ABI_ALLOCATE(rhowfg,(2,dtset%nfft))
 817      ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden))
 818      call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
 819 &     mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
 820      call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor)
 821      ABI_DEALLOCATE(rhowfg)
 822      ABI_DEALLOCATE(rhowfr)
 823    else
 824      call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
 825 &     mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
 826    end if
 827  end if ! getden
 828 
 829 !In PAW, compensation density has eventually to be added
 830  nhatgrdim=0;nhatdim=0
 831  ABI_ALLOCATE(nhatgr,(0,0,0))
 832  if (psps%usepaw==1.and. ((usexcnhat==0).or.(dtset%getden==0).or.dtset%xclevel==2)) then
 833    nhatdim=1
 834    ABI_ALLOCATE(nhat,(nfftf,dtset%nspden))
 835    call timab(558,1,tsec)
 836    nhatgrdim=0;if (dtset%xclevel==2.and.dtset%pawnhatxc>0) nhatgrdim=usexcnhat
 837    ider=2*nhatgrdim
 838    if (nhatgrdim>0)  then
 839      ABI_DEALLOCATE(nhatgr)
 840      ABI_ALLOCATE(nhatgr,(nfftf,dtset%nspden,3))
 841    end if
 842    izero=0;cplex=1;ipert=0;idir=0;qphon(:)=zero
 843    call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,&
 844 &   nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,&
 845 &   nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, &
 846 &   mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom)
 847    if (dtset%getden==0) then
 848      rhor(:,:)=rhor(:,:)+nhat(:,:)
 849      call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
 850    end if
 851    call timab(558,2,tsec)
 852  else
 853    ABI_ALLOCATE(nhat,(0,0))
 854  end if
 855 
 856 !The GS irrzon and phnons were only needed to symmetrize the GS density
 857  ABI_DEALLOCATE(irrzon)
 858  ABI_DEALLOCATE(phnons)
 859 
 860 !!jmb 2012 write(std_out,'(a)')' ' ! needed to make ibm6_xlf12 pass tests. No idea why this works. JWZ 5 Sept 2011
 861 !!Will compute now the total potential
 862 
 863 !Compute local ionic pseudopotential vpsp and core electron density xccc3d:
 864  n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf
 865  ABI_ALLOCATE(xccc3d,(n3xccc))
 866  ABI_ALLOCATE(vpsp,(nfftf))
 867 
 868  eei = zero
 869  if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
 870 !  PAW or NC with nc_xccc_gspace: compute Vloc and core charge together in reciprocal space
 871    call timab(562,1,tsec)
 872    optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optv=1;optn=n3xccc/nfftf;optn2=1
 873    call atm2fft(atindx1,xccc3d,vpsp,dum_dyfrn,dum_dyfrv,dum_eltfrxc,dum_gauss,gmet,gprimd,&
 874 &   dum_grn,dum_grv,gsqcut,mgfftf,psps%mqgrid_vl,natom,nattyp,nfftf,ngfftf,&
 875 &   ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,&
 876 &   dtset%qprtrb,dum_rhog,dummy6,dummy6,ucvol,psps%usepaw,dum_vg,dum_vg,dum_vg,dtset%vprtrb,psps%vlspl)
 877    call timab(562,2,tsec)
 878  else
 879 !  Norm-cons.: compute Vloc in reciprocal space and core charge in real space
 880    option=1
 881    ABI_ALLOCATE(dyfrlo_indx,(3,3,natom))
 882    ABI_ALLOCATE(grtn_indx,(3,natom))
 883    call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,&
 884 &   grtn_indx,gsqcut,dummy6,mgfftf,mpi_enreg,natom,nattyp,&
 885 &   nfftf,ngfftf,dtset%nspden,ntypat,option,pawtab,ph1df,psps,&
 886 &   dtset%qprtrb,rhog,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
 887    ABI_DEALLOCATE(dyfrlo_indx)
 888    ABI_DEALLOCATE(grtn_indx)
 889    if (psps%n1xccc/=0) then
 890      ABI_ALLOCATE(dyfrx2,(3,3,natom))
 891      ABI_ALLOCATE(vxc,(0,0)) ! dummy
 892      call mkcore(dummy6,dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,&
 893 &     ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,&
 894 &     psps%xcccrc,psps%xccc1d,xccc3d,xred)
 895      ABI_DEALLOCATE(dyfrx2)
 896      ABI_DEALLOCATE(vxc) ! dummy
 897    end if
 898  end if
 899 
 900 !Set up hartree and xc potential. Compute kxc here.
 901  ABI_ALLOCATE(vhartr,(nfftf))
 902  call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,rhog,rprimd,vhartr)
 903 
 904  option=3
 905  nkxc=2*dtset%nspden-1 ! LDA
 906  if(dtset%xclevel==2.and.dtset%nspden==1) nkxc=7  ! non-polarized GGA
 907  if(dtset%xclevel==2.and.dtset%nspden==2) nkxc=19 ! polarized GGA
 908  nk3xc=3*dtset%nspden-2
 909  ABI_ALLOCATE(kxc,(nfftf,nkxc))
 910  ABI_ALLOCATE(k3xc,(nfftf,nk3xc))
 911  ABI_ALLOCATE(vxc,(nfftf,dtset%nspden))
 912 
 913  _IBM6("Before rhotoxc")
 914 
 915  call status(0,dtfil%filstat,iexit,level,'call rhotoxc   ')
 916  call xcdata_init(xcdata,dtset=dtset)
 917  call rhotoxc(enxc,kxc,mpi_enreg,nfftf,ngfftf,&
 918 & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor,&
 919 & rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,k3xc=k3xc,vhartr=vhartr)
 920 
 921 !Compute local + Hxc potential, and subtract mean potential.
 922  ABI_ALLOCATE(vtrial,(nfftf,dtset%nspden))
 923  do ispden=1,min(dtset%nspden,2)
 924    do ifft=1,nfftf
 925      vtrial(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft)
 926    end do
 927  end do
 928  if (dtset%nspden==4) then
 929    do ispden=3,4
 930      do ifft=1,nfftf
 931        vtrial(ifft,ispden)=vxc(ifft,ispden)
 932      end do
 933    end do
 934  end if
 935  ABI_DEALLOCATE(vpsp)
 936  ABI_DEALLOCATE(vhartr)
 937 
 938  if(dtset%prtvol==-level)then
 939    call wrtout(std_out,' nonlinear : ground-state density and potential set up.','COLL')
 940  end if
 941 
 942  epaw = zero ; epawdc = zero
 943 !PAW: compute Dij quantities (psp strengths)
 944  if (psps%usepaw==1)then
 945    cplex=1;ipert=0;option=1
 946    nzlmopt=0;if (dtset%pawnzlm>0) nzlmopt=-1
 947    call pawdenpot(compch_sph,epaw,epawdc,ipert,dtset%ixc,my_natom,natom,dtset%nspden,&
 948 &   ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,&
 949 &   pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,&
 950 &   dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp, &
 951 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 952 
 953    call timab(561,1,tsec)
 954    call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,natom,nfftf,nfftotf,&
 955 &   dtset%nspden,ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,&
 956 &   pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,&
 957 &   dtset%spnorbscl,ucvol,dtset%charge,vtrial,vxc,xred,nucdipmom=dtset%nucdipmom,&
 958 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 959    call symdij(gprimd,indsym,ipert,my_natom,natom,dtset%nsym,ntypat,0,&
 960 &   paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,&
 961 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
 962    call timab(561,2,tsec)
 963  end if
 964 
 965  ABI_DEALLOCATE(xccc3d)
 966 
 967 !  Determine the subset of symmetry operations (nsym1 operations)
 968 !  that leaves the perturbation invariant, and initialize corresponding arrays
 969 !  symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW)..
 970  nsym1 = 1
 971 ! symaf1_tmp(1) = 1
 972 ! symrl1_tmp(:,:,1) = dtset%symrel(:,:,1)
 973 ! tnons1_tmp(:,1) = 0_dp
 974    ABI_ALLOCATE(indsy1,(4,nsym1,dtset%natom))
 975    ABI_ALLOCATE(symrc1,(3,3,nsym1))
 976    ABI_ALLOCATE(symaf1,(nsym1))
 977    ABI_ALLOCATE(symrl1,(3,3,nsym1))
 978    ABI_ALLOCATE(tnons1,(3,nsym1))
 979    symaf1(1)= 1 !symaf1_tmp(1:nsym1)
 980    symrl1(:,:,1)= dtset%symrel(:,:,1) !symrl1_tmp(:,:,1:nsym1)
 981    tnons1(:,1)= 0_dp !tnons1_tmp(:,1:nsym1)
 982 !   ABI_DEALLOCATE(symaf1_tmp)
 983 !   ABI_DEALLOCATE(symrl1_tmp)
 984 !   ABI_DEALLOCATE(tnons1_tmp)
 985 
 986 !  Set up corresponding symmetry data
 987  ABI_ALLOCATE(irrzon1,(dtset%nfft**(1-1/nsym1),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 988  ABI_ALLOCATE(phnons1,(2,dtset%nfft**(1-1/nsym1),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 989  call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,dtset%nspden,dtset%nsppol,&
 990 & nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred)
 991  if (psps%usepaw==1) then
 992 !  Allocate/initialize only zarot in pawang1 datastructure
 993    call pawang_init(pawang1,0,pawang%l_max-1,0,nsym1,0,1,0,0,0)
 994    call setsym_ylm(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot)
 995  end if
 996 
 997  if (pead/=0) then
 998 !  Initialize finite difference calculation of the ddk
 999 
1000    call status(0,dtfil%filstat,iexit,level,'call getshell ')
1001    nkpt3 = 0
1002 
1003 !  Prepare first call to getkgrid (obtain number of k points in FBZ)
1004    dum_kptrlatt(:,:) = dtset%kptrlatt(:,:)
1005    dum_nshiftk = dtset%nshiftk
1006    ABI_CHECK(dum_nshiftk<=210,"dum_nshiftk must be <= 210!")
1007    dum_shiftk(:,:) = zero
1008    dum_shiftk(:,1:dtset%nshiftk) = dtset%shiftk(:,1:dtset%nshiftk)
1009    dum_vacuum(:) = 0
1010 
1011    ABI_ALLOCATE(dum_kptns,(3,0))
1012    ABI_ALLOCATE(dum_wtk,(0))
1013    call getkgrid(0,0,dtset%iscf,dum_kptns,3,dum_kptrlatt,&
1014 &   rdum,dtset%nsym,0,nkpt3,dum_nshiftk,dtset%nsym,&
1015 &   rprimd,dum_shiftk,dtset%symafm,dtset%symrel,&
1016 &   dum_vacuum,dum_wtk)
1017    ABI_DEALLOCATE(dum_kptns)
1018    ABI_DEALLOCATE(dum_wtk)
1019 
1020 !   write(std_out,*) 'nonlinear : nkpt, nkpt3 = ',dtset%nkpt,nkpt3
1021 !call flush(6)
1022 !jmb : malloc() problem with gcc461_openmpi under max2 : change order of allocations works ?!?
1023 !allocate(kneigh(30,nkpt),kg_neigh(30,nkpt,3),mvwtk(30,nkpt))
1024    ABI_ALLOCATE(kg_neigh,(30,dtset%nkpt,3))
1025    ABI_ALLOCATE(mvwtk,(30,dtset%nkpt))
1026    ABI_ALLOCATE(kneigh,(30,dtset%nkpt))
1027 
1028    ABI_ALLOCATE(kptindex,(2,nkpt3))
1029    ABI_ALLOCATE(kpt3,(3,nkpt3))
1030 
1031    call getshell(gmet,kneigh,kg_neigh,kptindex,dtset%kptopt,&
1032 &   dtset%kptrlatt,dtset%kptns,kpt3,dtset%mkmem,mkmem_max,mvwtk,&
1033 &   dtset%nkpt,nkpt3,nneigh,dtset%nshiftk,rmet,rprimd,dtset%shiftk,dtset%wtk, mpi_enreg%comm_cell)
1034 
1035    ABI_ALLOCATE(pwind,(dtset%mpw,nneigh,dtset%mkmem))
1036    ABI_ALLOCATE(cgindex,(dtset%nkpt,dtset%nsppol))
1037    ABI_ALLOCATE(mpi_enreg%kpt_loc2ibz_sp,(0:mpi_enreg%nproc-1,1:mkmem_max, 1:2))
1038    ABI_ALLOCATE(mpi_enreg%mkmem,(0:mpi_enreg%nproc-1))
1039 
1040    call status(0,dtfil%filstat,iexit,level,'call initmv   ')
1041    call initmv(cgindex,dtset,gmet,kg,kneigh,kg_neigh,kptindex,&
1042 &   kpt3,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nband,dtset%nkpt,&
1043 &   nkpt3,nneigh,npwarr,dtset%nsppol,occ,pwind)
1044 
1045   call status(0,dtfil%filstat,iexit,level,'call pead_nl_loop ')
1046   call pead_nl_loop(blkflg,cg,cgindex,dtfil,dtset,d3etot,gmet,gprimd,gsqcut,&
1047 &  hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,dtset%mband,dtset%mgfft,&
1048 &  dtset%mkmem,mkmem_max,dtset%mk1mem,mpert,mpi_enreg,dtset%mpw,mvwtk,natom,nfftf,&
1049 &  dtset%nkpt,nkpt3,nkxc,nk3xc,nneigh,dtset%nspinor,dtset%nsppol,npwarr,occ,psps,pwind,&
1050 &  rfpert,rprimd,ucvol,xred)
1051 
1052  else ! pead=0 in this case
1053 
1054    call status(0,dtfil%filstat,iexit,level,'call dfptnl_loop ')
1055    call dfptnl_loop(atindx,blkflg,cg,dtfil,dtset,d3etot,eigen0,gmet,gprimd,gsqcut,&
1056 &   hdr,kg,kxc,k3xc,dtset%mband,dtset%mgfft,mgfftf,&
1057 &   dtset%mkmem,dtset%mk1mem,mpert,mpi_enreg,dtset%mpw,natom,nattyp,ngfftf,nfftf,nhat,&
1058 &   dtset%nkpt,nkxc,nk3xc,dtset%nspinor,dtset%nsppol,npwarr,occ,&
1059 &   paw_an,paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
1060 &   ph1d,ph1df,psps,rfpert,rhog,rhor,rprimd,ucvol,usecprj,vtrial,vxc,xred,&
1061 &   nsym1,indsy1,symaf1,symrc1,&
1062 &   d3e_1,d3e_2,d3e_3,d3e_4,d3e_5,d3e_6,d3e_7,d3e_8,d3e_9)
1063 
1064  end if ! end pead/=0
1065 
1066  write(message,'(a,a,a)')ch10,&
1067 & ' --- Third order energy calculation completed --- ',ch10
1068  call wrtout(ab_out,message,'COLL')
1069 
1070 !Complete missing elements using symmetry operations
1071  blkflg_sav = blkflg
1072  call status(0,dtfil%filstat,iexit,level,'call d3sym    ')
1073  call d3sym(blkflg,d3etot,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1074 
1075  blkflg_tmp = blkflg_sav
1076  call d3sym(blkflg_tmp,d3e_1,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1077  blkflg_tmp = blkflg_sav
1078  call d3sym(blkflg_tmp,d3e_2,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1079  blkflg_tmp = blkflg_sav
1080  call d3sym(blkflg_tmp,d3e_3,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1081  blkflg_tmp = blkflg_sav
1082  call d3sym(blkflg_tmp,d3e_4,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1083  blkflg_tmp = blkflg_sav
1084  call d3sym(blkflg_tmp,d3e_5,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1085  blkflg_tmp = blkflg_sav
1086  call d3sym(blkflg_tmp,d3e_6,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1087  blkflg_tmp = blkflg_sav
1088  call d3sym(blkflg_tmp,d3e_7,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1089  blkflg_tmp = blkflg_sav
1090  call d3sym(blkflg_tmp,d3e_8,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1091  blkflg_tmp = blkflg_sav
1092  call d3sym(blkflg_tmp,d3e_9,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel)
1093 
1094 !Open the formatted derivative database file, and write the
1095 !preliminary information
1096  if (mpi_enreg%me == 0) then
1097    call status(0,dtfil%filstat,iexit,level,'call ioddb8_ou')
1098 
1099    dscrpt=' Note : temporary (transfer) database '
1100 
1101    call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,&
1102 &                    1,xred=xred,occ=occ)
1103 
1104    call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_ddb, dtfil%unddb)
1105 
1106    call ddb_hdr_free(ddb_hdr)
1107 
1108 !  Call main output routine
1109    call dfptnl_doutput(blkflg,d3etot,dtset%mband,mpert,dtset%nkpt,dtset%natom,dtset%ntypat,dtfil%unddb)
1110 
1111 !  Close DDB
1112    close(dtfil%unddb)
1113 
1114 !  Compute tensors related to third-order derivatives
1115    call nlopt(blkflg,carflg,d3etot,d3cart,gprimd,mpert,natom,rprimd,ucvol)
1116 !  Note that the imaginary part is not transformed into cartesian coordinates
1117    call nlopt(blkflg,carflg_tmp,d3e_1,d3cart_1,gprimd,mpert,natom,rprimd,ucvol)
1118    call nlopt(blkflg,carflg_tmp,d3e_2,d3cart_2,gprimd,mpert,natom,rprimd,ucvol)
1119    call nlopt(blkflg,carflg_tmp,d3e_3,d3cart_3,gprimd,mpert,natom,rprimd,ucvol)
1120    call nlopt(blkflg,carflg_tmp,d3e_4,d3cart_4,gprimd,mpert,natom,rprimd,ucvol)
1121    call nlopt(blkflg,carflg_tmp,d3e_5,d3cart_5,gprimd,mpert,natom,rprimd,ucvol)
1122    call nlopt(blkflg,carflg_tmp,d3e_6,d3cart_6,gprimd,mpert,natom,rprimd,ucvol)
1123    call nlopt(blkflg,carflg_tmp,d3e_7,d3cart_7,gprimd,mpert,natom,rprimd,ucvol)
1124    call nlopt(blkflg,carflg_tmp,d3e_8,d3cart_8,gprimd,mpert,natom,rprimd,ucvol)
1125    call nlopt(blkflg,carflg_tmp,d3e_9,d3cart_9,gprimd,mpert,natom,rprimd,ucvol)
1126 
1127    if ((d3e_pert1(natom+2)==1).and.(d3e_pert2(natom+2)==1).and. &
1128 &   (d3e_pert3(natom+2)==1)) then
1129 
1130      flag = 1
1131      i1pert = natom+2
1132 
1133      d3cart(:,:,i1pert,:,i1pert,:,i1pert) = &
1134 &     d3cart(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1135 
1136      write(ab_out,*)ch10
1137      write(ab_out,*)' Non-linear optical susceptibility tensor d (pm/V)'
1138      write(ab_out,*)' in cartesian coordinates'
1139      write(ab_out,*)'  i1dir  i2dir  i3dir             d'
1140 
1141      do i1dir = 1, 3
1142        do i2dir = 1, 3
1143          do i3dir = 1, 3
1144            write(ab_out,'(3(5x,i2),5x,f16.9)') i1dir,i2dir,i3dir,&
1145 &           d3cart(1,i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)
1146            if ((blkflg(i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)/=1).or.&
1147 &           (carflg(i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)/=1)) flag = 0
1148          end do
1149        end do
1150      end do
1151 
1152      if (pead==0.and.(dtset%nonlinear_info>0)) then
1153 
1154        d3cart_1(:,:,i1pert,:,i1pert,:,i1pert) = &
1155 &       d3cart_1(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1156        d3cart_2(:,:,i1pert,:,i1pert,:,i1pert) = &
1157 &       d3cart_2(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1158        d3cart_8(:,:,i1pert,:,i1pert,:,i1pert) = &
1159 &       d3cart_8(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1160        d3cart_9(:,:,i1pert,:,i1pert,:,i1pert) = &
1161 &       d3cart_9(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb
1162 
1163        theunit = ab_out
1164 
1165        write(small_msg,'(a)') ' ** Total :'
1166        call print_chi2(d3cart,small_msg,theunit)
1167 
1168        write(small_msg,'(a)') ' ** sum_psi1H1psi1 :'
1169        call print_chi2(d3cart_1,small_msg,theunit)
1170 
1171        write(small_msg,'(a)') ' ** sum_lambda1psi1psi1 :'
1172        call print_chi2(d3cart_2,small_msg,theunit)
1173 
1174        write(small_msg,'(a)') ' ** exc3 :'
1175        call print_chi2(d3cart_8,small_msg,theunit)
1176 
1177        write(small_msg,'(a)') ' ** exc3_paw :'
1178        call print_chi2(d3cart_9,small_msg,theunit)
1179 
1180      end if ! nonlinear_info > 0
1181 
1182      if (flag == 0) then
1183        write(message,'(a,a,a,a,a,a)')ch10,&
1184 &       ' dfptnl_doutput: WARNING -',ch10,&
1185 &       '  matrix of third-order energies incomplete,',ch10,&
1186 &       '  non-linear optical coefficients may be wrong.'
1187        call wrtout(ab_out,message,'COLL')
1188        call wrtout(std_out,message,'COLL')
1189      end if
1190 
1191    end if  ! d3e_pert1,d3e_pert2,d3e_pert3
1192 
1193    if (((maxval(d3e_pert1(1:natom))/=0).and.(d3e_pert2(natom+2)/=0).and. &
1194 &   (d3e_pert3(natom+2)/=0)).or.&
1195    ((maxval(d3e_pert2(1:natom))/=0).and.(d3e_pert1(natom+2)/=0).and. &
1196 &   (d3e_pert3(natom+2)/=0)).or.&
1197    ((maxval(d3e_pert3(1:natom))/=0).and.(d3e_pert2(natom+2)/=0).and. &
1198 &   (d3e_pert1(natom+2)/=0))) then
1199 !    Perform a check if all relevant elements are available
1200 
1201      flag = 1
1202      do i1pert = 1, natom
1203        do i1dir = 1, 3
1204          do i2dir = 1, 3
1205            do i3dir = 1, 3
1206              if ((blkflg(i1dir,i1pert,i2dir,natom+2,i3dir,natom+2) /= 1).or.&
1207              (blkflg(i1dir,natom+2,i2dir,i1pert,i3dir,natom+2) /= 1).or.&
1208              (blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,i1pert) /= 1)) flag = 0
1209              if ((carflg(i1dir,i1pert,i2dir,natom+2,i3dir,natom+2) /= 1).or.&
1210              (carflg(i1dir,natom+2,i2dir,i1pert,i3dir,natom+2) /= 1).or.&
1211              (carflg(i1dir,natom+2,i2dir,natom+2,i3dir,i1pert) /= 1)) flag = 0
1212            end do
1213          end do
1214        end do
1215      end do
1216 
1217      write(ab_out,*)ch10
1218      write(ab_out,*)' First-order change in the electronic dielectric '
1219      write(ab_out,*)' susceptibility tensor (Bohr^-1)'
1220      write(ab_out,*)' induced by an atomic displacement'
1221      write(ab_out,*)'  atom  displacement'
1222 
1223      do i1pert = 1,natom
1224        do i1dir = 1,3
1225          write(ab_out,'(1x,i4,9x,i2,3(3x,f16.9))')i1pert,i1dir,&
1226 &         d3cart(1,i1dir,i1pert,1,natom+2,:,natom+2)
1227          write(ab_out,'(16x,3(3x,f16.9))')&
1228 &         d3cart(1,i1dir,i1pert,2,natom+2,:,natom+2)
1229          write(ab_out,'(16x,3(3x,f16.9))')&
1230 &         d3cart(1,i1dir,i1pert,3,natom+2,:,natom+2)
1231        end do
1232        write(ab_out,*)
1233      end do
1234 
1235      if (flag == 0) then
1236        write(message,'(a,a,a,a,a,a)')ch10,&
1237 &       ' dfptnl_doutput: WARNING -',ch10,&
1238 &       '  matrix of third-order energies incomplete,',ch10,&
1239 &       '  changes in the dielectric susceptibility may be wrong.'
1240        call wrtout(ab_out,message,'COLL')
1241        call wrtout(std_out,message,'COLL')
1242      end if
1243 
1244      if (pead==0.and.(dtset%nonlinear_info>0)) then
1245        theunit = ab_out
1246 
1247        write(small_msg,'(a)') ' ** Total :'
1248        call print_dchidtau(d3cart,small_msg,theunit)
1249 
1250        write(small_msg,'(a)') ' ** sum_psi1H1psi1 :'
1251        call print_dchidtau(d3cart_1,small_msg,theunit)
1252 
1253        write(small_msg,'(a)') ' ** sum_lambda1psi1psi1 :'
1254        call print_dchidtau(d3cart_2,small_msg,theunit)
1255 
1256        write(small_msg,'(a)') ' ** sum_lambda1psi0S1psi1 :'
1257        call print_dchidtau(d3cart_3,small_msg,theunit)
1258 
1259        write(small_msg,'(a)') ' ** sum_psi0H2psi1a :'
1260        call print_dchidtau(d3cart_4,small_msg,theunit)
1261 
1262        write(small_msg,'(a)') ' ** sum_psi0H2psi1b :'
1263        call print_dchidtau(d3cart_5,small_msg,theunit)
1264 
1265        write(small_msg,'(a)') ' ** eHxc21_paw :'
1266        call print_dchidtau(d3cart_6,small_msg,theunit)
1267 
1268        write(small_msg,'(a)') ' ** eHxc21_nhat :'
1269        call print_dchidtau(d3cart_7,small_msg,theunit)
1270 
1271        write(small_msg,'(a)') ' ** exc3 :'
1272        call print_dchidtau(d3cart_8,small_msg,theunit)
1273 
1274        write(small_msg,'(a)') ' ** exc3_paw :'
1275        call print_dchidtau(d3cart_9,small_msg,theunit)
1276 
1277      end if ! nonlinear_info > 0
1278 
1279    end if  ! d3e_pert1,d3e_pert2,d3e_pert3
1280  end if   ! mpi_enreg%me
1281 
1282 ! TO OPTIMIZE DEALLOCATION !
1283  if (pead/=0) then
1284    ABI_DEALLOCATE(cgindex)
1285    ABI_DEALLOCATE(kg_neigh)
1286    ABI_DEALLOCATE(kneigh)
1287    ABI_DEALLOCATE(kptindex)
1288    ABI_DEALLOCATE(kpt3)
1289    ABI_DEALLOCATE(mpi_enreg%kpt_loc2ibz_sp)
1290    ABI_DEALLOCATE(mpi_enreg%mkmem)
1291    ABI_DEALLOCATE(mvwtk)
1292    ABI_DEALLOCATE(pwind)
1293  end if
1294  ABI_DEALLOCATE(atindx)
1295  ABI_DEALLOCATE(atindx1)
1296  ABI_DEALLOCATE(blkflg)
1297  ABI_DEALLOCATE(blkflg_tmp)
1298  ABI_DEALLOCATE(blkflg_sav)
1299  ABI_DEALLOCATE(carflg)
1300  ABI_DEALLOCATE(carflg_tmp)
1301  ABI_DEALLOCATE(carflg_sav)
1302  ABI_DEALLOCATE(cg)
1303  ABI_DEALLOCATE(d3cart)
1304  ABI_DEALLOCATE(d3etot)
1305  ABI_DEALLOCATE(d3cart_1)
1306  ABI_DEALLOCATE(d3cart_2)
1307  ABI_DEALLOCATE(d3cart_3)
1308  ABI_DEALLOCATE(d3cart_4)
1309  ABI_DEALLOCATE(d3cart_5)
1310  ABI_DEALLOCATE(d3cart_6)
1311  ABI_DEALLOCATE(d3cart_7)
1312  ABI_DEALLOCATE(d3cart_8)
1313  ABI_DEALLOCATE(d3cart_9)
1314  ABI_DEALLOCATE(d3e_1)
1315  ABI_DEALLOCATE(d3e_2)
1316  ABI_DEALLOCATE(d3e_3)
1317  ABI_DEALLOCATE(d3e_4)
1318  ABI_DEALLOCATE(d3e_5)
1319  ABI_DEALLOCATE(d3e_6)
1320  ABI_DEALLOCATE(d3e_7)
1321  ABI_DEALLOCATE(d3e_8)
1322  ABI_DEALLOCATE(d3e_9)
1323  ABI_DEALLOCATE(d3e_pert1)
1324  ABI_DEALLOCATE(d3e_pert2)
1325  ABI_DEALLOCATE(d3e_pert3)
1326  ABI_DEALLOCATE(eigen0)
1327  ABI_DEALLOCATE(rhog)
1328  ABI_DEALLOCATE(rhor)
1329  ABI_DEALLOCATE(nhat)
1330  ABI_DEALLOCATE(nhatgr)
1331  ABI_DEALLOCATE(rfpert)
1332  ABI_DEALLOCATE(grxc)
1333  ABI_DEALLOCATE(kg)
1334  ABI_DEALLOCATE(kxc)
1335  ABI_DEALLOCATE(k3xc)
1336  ABI_DEALLOCATE(indsym)
1337  ABI_DEALLOCATE(indsy1)
1338  ABI_DEALLOCATE(nattyp)
1339  ABI_DEALLOCATE(npwarr)
1340  ABI_DEALLOCATE(symrec)
1341  ABI_DEALLOCATE(symrc1)
1342  ABI_DEALLOCATE(symaf1)
1343  ABI_DEALLOCATE(symrl1)
1344  ABI_DEALLOCATE(tnons1)
1345  ABI_DEALLOCATE(irrzon1)
1346  ABI_DEALLOCATE(phnons1)
1347  ABI_DEALLOCATE(symq)
1348  ABI_DEALLOCATE(ph1d)
1349  ABI_DEALLOCATE(ph1df)
1350  ABI_DEALLOCATE(vtrial)
1351  ABI_DEALLOCATE(vxc)
1352  call pawfgr_destroy(pawfgr)
1353  if (psps%usepaw==1) then
1354    call pawang_free(pawang1)
1355    call pawrhoij_free(pawrhoij)
1356    call paw_an_free(paw_an)
1357    call paw_ij_free(paw_ij)
1358    call pawfgrtab_free(pawfgrtab)
1359  end if
1360  ABI_DATATYPE_DEALLOCATE(pawrhoij)
1361  ABI_DATATYPE_DEALLOCATE(paw_an)
1362  ABI_DATATYPE_DEALLOCATE(paw_ij)
1363  ABI_DATATYPE_DEALLOCATE(pawfgrtab)
1364 
1365 !Clean the header
1366  call hdr_free(hdr)
1367 
1368 !As the etotal energy has no meaning here, we set it to zero
1369 !(to avoid meaningless side-effects when comparing ouputs...)
1370  etotal = zero
1371 
1372  call status(0,dtfil%filstat,iexit,level,' exit         ')
1373  call timab(501,2,tsec)
1374 
1375  DBG_EXIT("COLL")
1376 
1377  contains

nonlinear/print_chi2 [ Functions ]

[ Top ] [ nonlinear ] [ Functions ]

NAME

 print_chi2

FUNCTION

 Print a third derivative tensor. Used only in nonlinear

INPUTS

  d3cart0 = the tensor to print
  msg = a short message printed before the tensor
  theunit = unit where the tensor is written

OUTPUT

SIDE EFFECTS

PARENTS

      nonlinear

CHILDREN

SOURCE

1403    subroutine print_chi2(d3cart0,msg,theunit)
1404 
1405 
1406 !This section has been created automatically by the script Abilint (TD).
1407 !Do not modify the following lines by hand.
1408 #undef ABI_FUNC
1409 #define ABI_FUNC 'print_chi2'
1410 !End of the abilint section
1411 
1412      implicit none
1413 
1414      integer,intent(in) :: theunit
1415      character(len=30) :: msg
1416      real(dp) :: elem1,elem2
1417      real(dp),intent(in) :: d3cart0(2,3,mpert,3,mpert,3,mpert)
1418 
1419 ! *************************************************************************
1420 
1421      write(theunit,'(2a)') ch10,msg
1422      do i1dir = 1, 3
1423        do i2dir = 1, 3
1424          do i3dir = 1, 3
1425            elem1 = d3cart0(1,i1dir,natom+2,i2dir,natom+2,i3dir,natom+2)
1426            elem2 = d3cart0(2,i1dir,natom+2,i2dir,natom+2,i3dir,natom+2)
1427            write(theunit,'(3(5x,i2),5x,f16.9,2x,f16.9)') i1dir,i2dir,i3dir,elem1,elem2
1428          end do
1429        end do
1430      end do
1431 
1432    end subroutine print_chi2

nonlinear/print_dchidtau [ Functions ]

[ Top ] [ nonlinear ] [ Functions ]

NAME

 print_dchidtau

FUNCTION

 Print a third derivative tensor. Used only in nonlinear

INPUTS

  d3cart0 = the tensor to print
  msg = a short message printed before the tensor
  theunit = unit where the tensor is written

OUTPUT

SIDE EFFECTS

PARENTS

      nonlinear

CHILDREN

SOURCE

1458    subroutine print_dchidtau(d3cart0,msg,theunit)
1459 
1460 
1461 !This section has been created automatically by the script Abilint (TD).
1462 !Do not modify the following lines by hand.
1463 #undef ABI_FUNC
1464 #define ABI_FUNC 'print_dchidtau'
1465 !End of the abilint section
1466 
1467      implicit none
1468 
1469      integer,intent(in) :: theunit
1470      character(len=30) :: msg
1471      real(dp),intent(in) :: d3cart0(2,3,mpert,3,mpert,3,mpert)
1472 
1473 ! *************************************************************************
1474 
1475      write(theunit,'(a)') msg
1476      do i1pert = 1,natom
1477        do i1dir = 1,3
1478          write(theunit,'(1x,i4,9x,i2,3(3x,f16.9),3(3x,f16.9))')i1pert,i1dir,&
1479   &      d3cart0(1,i1dir,i1pert,1,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,1,natom+2,:,natom+2)
1480          write(theunit,'(16x,3(3x,f16.9),3(3x,f16.9))')&
1481   &      d3cart0(1,i1dir,i1pert,2,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,2,natom+2,:,natom+2)
1482          write(theunit,'(16x,3(3x,f16.9),3(3x,f16.9))')&
1483   &      d3cart0(1,i1dir,i1pert,3,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,3,natom+2,:,natom+2)
1484        end do
1485      end do
1486 
1487   end subroutine print_dchidtau