TABLE OF CONTENTS


ABINIT/dfpt_nstdy [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_nstdy

FUNCTION

 This routine compute the non-stationary expression for the
 second derivative of the total energy, for a whole row of
 mixed derivatives.
 Only for norm-conserving pseudopotentials (no PAW)

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (XG, DCA, GMR, MM, AR, MV, MB, MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k
  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
  eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  idir=direction of the perturbation
  indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points
  indsy1(4,nsym1,natom)=indirect indexing array for atom labels
  ipert=type of the perturbation
  istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ
  kxc(nfft,nkxc)=exchange and correlation kernel
  mkmem =number of k points treated by this node (GS data)
  mk1mem =number of k points treated by this node (RF data)
  mpert =maximum number of ipert
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw or wfs at k
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  nattyp(ntypat)= # atoms of each type.
  nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
  nfft=(effective) number of FFT grid points (for this proc)
  ngfft(1:18)=integer array with FFT box dimensions and other
  nkpt=number of k points in the full BZ
  nkpt_rbz=number of k points in the reduced BZ for this perturbation
  nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
  npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
  npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym1=number of symmetry elements in space group consistent with i perturbation
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band
   and k in the reduced Brillouin zone (usually =2)
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhor1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3.
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  symrc1(3,3,nsym1)=symmetry operations in reciprocal space
  ucvol=unit cell volume in bohr**3.
  wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point

OUTPUT

  blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
  d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
                                        second order derivatives
  d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
  d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs

NOTES

 Note that the ddk perturbation should not be treated here.

PARENTS

      dfpt_scfcv

CHILDREN

      appdig,destroy_hamiltonian,dfpt_mkcore,dfpt_mkvxc,dfpt_mkvxc_noncoll
      dfpt_nstwf,dfpt_sygra,dfpt_vlocal,dotprod_vn,init_hamiltonian
      load_spin_hamiltonian,mati3inv,timab,wfk_close,wfk_open_read,wrtout
      xmpi_sum

SOURCE

 92 #if defined HAVE_CONFIG_H
 93 #include "config.h"
 94 #endif
 95 
 96 #include "abi_common.h"
 97 
 98 subroutine dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,&
 99 &          gmet,gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mkmem,mk1mem,&
100 &          mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,nfft,ngfft,nkpt,nkpt_rbz,nkxc,&
101 &          npwarr,npwar1,nspden,nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,&
102 &          symrc1,ucvol,wtk_rbz,xred,ylm,ylm1,rhor,vxc)
103 
104  use defs_basis
105  use defs_datatypes
106  use defs_abitypes
107  use m_profiling_abi
108  use m_xmpi
109  use m_errors
110  use m_wfk
111  use m_nctk
112  use m_hamiltonian
113 
114  use m_io_tools,  only : file_exists
115  use m_cgtools,   only : dotprod_vn
116  use m_hdr,       only : hdr_skip
117  use m_pawtab,    only : pawtab_type
118 
119 !This section has been created automatically by the script Abilint (TD).
120 !Do not modify the following lines by hand.
121 #undef ABI_FUNC
122 #define ABI_FUNC 'dfpt_nstdy'
123  use interfaces_14_hidewrite
124  use interfaces_18_timing
125  use interfaces_32_util
126  use interfaces_56_xc
127  use interfaces_72_response, except_this_one => dfpt_nstdy
128 !End of the abilint section
129 
130  implicit none
131 
132 !Arguments -------------------------------
133 !scalars
134  integer,intent(in) :: cplex,idir,ipert,mk1mem,mkmem,mpert,mpw,mpw1,nfft,nkpt,nkpt_rbz,nkxc,nspden,nsppol,nsym1
135  real(dp),intent(in) :: gsqcut,ucvol
136  type(MPI_type),intent(in) :: mpi_enreg
137  type(datafiles_type),intent(in) :: dtfil
138  type(dataset_type),intent(in) :: dtset
139  type(pseudopotential_type),intent(in) :: psps
140 !arrays
141  integer,intent(in) :: atindx(dtset%natom),indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom)
142  integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
143  integer,intent(in) :: nattyp(dtset%ntypat),nband_rbz(nkpt_rbz*nsppol),ngfft(18)
144  integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz),symrc1(3,3,nsym1)
145  integer,intent(inout) :: blkflg(3,mpert,3,mpert) !vz_i
146  real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
147  real(dp),intent(in) :: cg1(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*nsppol)
148  real(dp),intent(in) :: eigen0(dtset%mband*nkpt_rbz*nsppol)
149  real(dp),intent(in) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*nsppol)
150  real(dp),intent(in) :: gmet(3,3),kpt_rbz(3,nkpt_rbz)
151  real(dp),intent(in) :: kxc(nfft,nkxc),occ_rbz(dtset%mband*nkpt_rbz*nsppol)
152  real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
153  real(dp),intent(in) :: rhor1(cplex*nfft,nspden),rmet(3,3),rprimd(3,3)
154  real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,dtset%natom)
155  real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
156  real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
157  real(dp),intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*dtset%prtbbb)!vz_i
158  real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i
159 ! optional
160  real(dp),optional,intent(in) :: rhor(nfft,nspden)
161  real(dp),optional,intent(in) :: vxc(cplex*nfft,nspden)
162 
163 !Local variables-------------------------------
164 !scalars
165  integer,parameter :: formeig1=1
166  integer :: ban2tot,bantot,bdtot_index,ddkcase,iband,icg,icg1,idir1
167  integer :: ierr,ifft,ii,ikg,ikg1,ikpt,ilm,ipert1,ispden,isppol
168  integer :: istwf_k,isym,jj,master,me,n1,n2,n3,n3xccc,n4,n5,n6
169  integer :: nband_k,nfftot,npw1_k,npw_k,nspinor_,option,spaceworld,optnc
170  real(dp) :: doti,dotr,wtk_k
171  logical :: t_exist
172  character(len=500) :: msg
173  character(len=fnlen) :: fiwfddk
174  type(gs_hamiltonian_type) :: gs_hamkq
175 !arrays
176  integer :: ddkfil(3)
177  integer,allocatable :: kg1_k(:,:),kg_k(:,:),symrl1(:,:,:)
178  real(dp) :: d2nl_elfd(2,3),d2nl_mgfd(2,3),kpoint(3),kpq(3),sumelfd(2),summgfd(2),tsec(2)
179  real(dp),allocatable :: buffer1(:),buffer2(:),d2bbb_k(:,:,:,:),d2nl_k(:,:,:)
180  real(dp),allocatable :: eig1_k(:),eig_k(:),occ_k(:)
181  real(dp) :: rhodummy(0,0)
182  real(dp),allocatable :: vpsp1(:),vxc1(:,:),work1(:,:,:),xccc3d1(:),ylm1_k(:,:),ylm_k(:,:)
183  type(pawtab_type) :: pawtab(dtset%ntypat*psps%usepaw)
184  type(wfk_t) :: ddks(3)
185 
186 ! *********************************************************************
187 
188  ABI_UNUSED(nkpt)
189 
190  DBG_ENTER("COLL")
191 
192 !Not valid for PAW
193  if (psps%usepaw==1) then
194    msg='This routine cannot be used for PAW (use pawnst3 instead) !'
195    MSG_BUG(msg)
196  end if
197 
198 !Keep track of total time spent in dfpt_nstdy
199  call timab(101,1,tsec)
200 
201 !Init parallelism
202  spaceworld=mpi_enreg%comm_cell
203  me=mpi_enreg%me_kpt
204 
205  master =0
206 
207 !Zero only portion of nonlocal matrix to be computed here
208  d2nl(:,:,1:dtset%natom+2,idir,ipert)=zero
209 
210  ABI_ALLOCATE(d2bbb_k,(2,3,dtset%mband,dtset%mband*dtset%prtbbb))
211  ABI_ALLOCATE(d2nl_k,(2,3,mpert))
212  ABI_ALLOCATE(eig_k,(nsppol*dtset%mband))
213  ABI_ALLOCATE(eig1_k,(2*nsppol*dtset%mband**2))
214  ABI_ALLOCATE(kg_k,(3,mpw))
215  ABI_ALLOCATE(kg1_k,(3,mpw1))
216 
217 !Do not try to open electric field file
218  ddkfil(:)=0
219 !The treatment of homogeneous electric field potential need the existence of d/dk files.
220  do idir1=1,3
221    ddkcase=idir1+dtset%natom*3
222    call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk)
223 
224 !  Check that ddk file exists
225    t_exist = file_exists(fiwfddk)
226    if (.not. t_exist) then
227      ! Try netcdf file.
228      t_exist = file_exists(nctk_ncify(fiwfddk))
229      if (t_exist) then
230        fiwfddk = nctk_ncify(fiwfddk)
231        write(msg,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name."
232        call wrtout(std_out,msg,'COLL')
233      end if
234    end if
235 
236    if (t_exist) then
237 !    Note the use of unit numbers 21, 22 and 23
238      ddkfil(idir1)=20+idir1
239      write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk)
240      call wrtout(std_out,msg,'COLL')
241      call wrtout(ab_out,msg,'COLL')
242      call wfk_open_read(ddks(idir1),fiwfddk,formeig1,dtset%iomode,ddkfil(idir1),spaceworld)
243    end if
244  end do
245 
246 !Update list of computed matrix elements
247  if (ipert /= dtset%natom + 1) then
248    do ipert1=1,mpert
249      do idir1=1,3
250        if(ipert1 <= dtset%natom .or. ipert1==dtset%natom+2 &
251 &       .and. ddkfil(idir1)/=0) then
252          blkflg(idir1,ipert1,idir,ipert)=1
253        end if
254      end do
255    end do
256  else
257    ipert1 = dtset%natom + 1
258    do idir1=1,3
259 !    If was already computed in another run or dataset, or if is to be computed in the present one
260      if ((ddkfil(idir1) /= 0).or. (dtset%rfdir(idir1)/=0.and. idir1<=idir) ) then
261 !      if ((ddkfil(idir1) /= 0).or. (idir1==idir) ) then
262        blkflg(idir1,ipert1,idir,ipert)=1
263      end if
264    end do
265  end if
266 
267  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
268  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
269  nspinor_=dtset%nspinor
270 
271  bantot = 0
272  ban2tot = 0
273 
274 !==== Initialize most of the Hamiltonian ====
275 !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
276 !2) Perform the setup needed for the non-local factors:
277 !3) Constant kleimann-Bylander energies are copied from psps to gs_hamk.
278  call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,dtset%natom,&
279 & dtset%typat,xred,nfft,dtset%mgfft,ngfft,rprimd,dtset%nloalg,ph1d=ph1d,&
280 & use_gpu_cuda=dtset%use_gpu_cuda)
281 
282 !LOOP OVER SPINS
283  bdtot_index=0
284  icg=0;icg1=0
285  do isppol=1,nsppol
286 
287    ikg=0;ikg1=0
288 
289 !  Continue to initialize the Hamiltonian
290    call load_spin_hamiltonian(gs_hamkq,isppol,with_nonlocal=.true.)
291 
292 !  BIG FAT k POINT LOOP
293    do ikpt=1,nkpt_rbz
294 
295      nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
296      istwf_k=istwfk_rbz(ikpt)
297      npw_k=npwarr(ikpt)
298      npw1_k=npwar1(ikpt)
299 
300      eig_k(1:nband_k) = eigen0(1+bantot:nband_k+bantot)
301      eig1_k(1:2*nband_k**2) = eigen1(1+ban2tot:2*nband_k**2+ban2tot)
302      bantot = bantot + nband_k
303      ban2tot = ban2tot + 2*nband_k**2
304 
305      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
306        bdtot_index=bdtot_index+nband_k
307 !      The wavefunction blocks for ddk file is skipped elsewhere in the loop
308 !      Skip the rest of the k-point loop
309        cycle
310      end if
311 
312      ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
313      ABI_ALLOCATE(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
314 
315 !    In case of electric field pert1, read ddk wfs file
316 !    Note that the symmetries are not used for ddk, so read each k point
317 !    Also take into account implicitely the parallelism over k points
318 
319      do idir1=1,3
320        if (ddkfil(idir1)/=0) then
321          ii = wfk_findk(ddks(idir1), kpt_rbz(:, ikpt))
322          ABI_CHECK(ii == indkpt1(ikpt),  "ii !=  indkpt1")
323        end if
324      end do
325 
326      ABI_ALLOCATE(occ_k,(nband_k))
327      occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
328      kpoint(:)=kpt_rbz(:,ikpt)
329      kpq(:)=kpoint(:)+dtset%qptn(:)
330      wtk_k=wtk_rbz(ikpt)
331      d2nl_k(:,:,:)=zero
332      if(dtset%prtbbb==1)d2bbb_k(:,:,:,:)=zero
333 
334 !    Get plane-wave vectors and related data at k
335      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
336      if (psps%useylm==1) then
337        do ilm=1,psps%mpsang*psps%mpsang
338          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
339        end do
340      end if
341 
342 !    Get plane-wave vectors and related data at k+q
343      kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
344      if (psps%useylm==1) then
345        do ilm=1,psps%mpsang*psps%mpsang
346          ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
347        end do
348      end if
349 
350 !    Compute the eigenvalues, wavefunction,
351 !    contributions to kinetic energy, nonlocal energy, forces,
352 !    and update of rhor1 to this k-point and this spin polarization.
353 !    Note that dfpt_nstwf is called with kpoint, while kpt is used inside dfpt_vtowfk
354      call dfpt_nstwf(cg,cg1,ddkfil,dtset,d2bbb_k,d2nl_k,eig_k,eig1_k,gs_hamkq,&
355 &     icg,icg1,idir,ikpt,ipert,isppol,istwf_k,kg_k,kg1_k,kpoint,kpq,mkmem,mk1mem,mpert,&
356 &     mpi_enreg,mpw,mpw1,nband_k,npw_k,npw1_k,nsppol,&
357 &     occ_k,psps,rmet,ddks,wtk_k,ylm_k,ylm1_k)
358 
359      d2nl(:,:,:,idir,ipert)=d2nl(:,:,:,idir,ipert)+d2nl_k(:,:,:)
360      if(dtset%prtbbb==1)d2bbb(:,:,idir,ipert,:,:)=d2bbb(:,:,idir,ipert,:,:)+d2bbb_k(:,:,:,:)
361 
362 !    Keep track of total number of bands
363      bdtot_index=bdtot_index+nband_k
364 
365 !    Shift arrays memory
366      if (mkmem/=0) then
367        icg=icg+npw_k*dtset%nspinor*nband_k
368        ikg=ikg+npw_k
369      end if
370      if (mk1mem/=0) then
371        icg1=icg1+npw1_k*dtset%nspinor*nband_k
372        ikg1=ikg1+npw1_k
373      end if
374 
375      ABI_DEALLOCATE(occ_k)
376      ABI_DEALLOCATE(ylm_k)
377      ABI_DEALLOCATE(ylm1_k)
378 
379 !    End big k point loop
380    end do
381 
382 !  End loop over spins
383  end do
384 
385 !if(xmpi_paral==1)then
386 !call timab(161,1,tsec)
387 !call wrtout(std_out,' dfpt_nstdy: loop on k-points and spins done in parallel','COLL')
388 !call xmpi_barrier(spaceworld)
389 !call timab(161,2,tsec)
390 !end if
391 
392  call destroy_hamiltonian(gs_hamkq)
393 
394 !Treat fixed occupation numbers (as in vtorho)
395  if(xmpi_paral==1)then
396    ABI_ALLOCATE(buffer1,(2*3*mpert))
397    ABI_ALLOCATE(buffer2,(2*3*mpert))
398 !  Pack d2nl
399    buffer1(1:2*3*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/2*3*mpert/))
400 !  Build sum of everything
401    call timab(48,1,tsec)
402    call xmpi_sum(buffer1,buffer2,2*3*mpert,spaceworld,ierr)
403    call timab(48,2,tsec)
404 !  Unpack the final result
405    d2nl(:,:,:,idir,ipert)=reshape(buffer2(:),(/2,3,mpert/))
406    ABI_DEALLOCATE(buffer1)
407    ABI_DEALLOCATE(buffer2)
408    if(dtset%prtbbb==1)then
409      ABI_ALLOCATE(buffer1,(2*3*dtset%mband*dtset%mband))
410      ABI_ALLOCATE(buffer2,(2*3*dtset%mband*dtset%mband))
411 !    Pack d2bbb
412      buffer1(1:2*3*dtset%mband*dtset%mband)=reshape(d2bbb(:,:,idir,ipert,:,:),(/2*3*dtset%mband*dtset%mband/))
413 !    Build sum of everything
414      call timab(48,1,tsec)
415      call xmpi_sum(buffer1,buffer2,2*3*dtset%mband*dtset%mband,spaceworld,ierr)
416      call timab(48,2,tsec)
417 !    Unpack the final result
418      d2bbb(:,:,idir,ipert,:,:)=reshape(buffer2(:),(/2,3,dtset%mband,dtset%mband/))
419      ABI_DEALLOCATE(buffer1)
420      ABI_DEALLOCATE(buffer2)
421    end if
422  end if ! xmpi_paral==1
423 
424 !In the case of the strain perturbation time-reversal symmetry will always
425 !be true so imaginary part of d2nl will be must be set to zero here since
426 !the symmetry-reduced kpt set will leave a non-zero imaginary part.
427  if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) d2nl(2,:,:,idir,ipert)=zero
428 
429 !In case of electric field ipert1, close the ddk wf files
430  do idir1=1,3
431    if (ddkfil(idir1)/=0)then
432      call wfk_close(ddks(idir1))
433    end if
434  end do
435 
436 !Symmetrize the non-local contributions,
437 !as was needed for the forces in a ground-state calculation
438 !However, here the quantity is complex, and there are phases !
439 
440 !Do the transform
441  ABI_ALLOCATE(work1,(2,3,dtset%natom))
442  do ipert1=1,dtset%natom
443    do idir1=1,3
444      work1(1,idir1,ipert1)=d2nl(1,idir1,ipert1,idir,ipert)
445      work1(2,idir1,ipert1)=d2nl(2,idir1,ipert1,idir,ipert)
446    end do
447  end do
448  call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work1,indsy1,ipert,nsym1,dtset%qptn,symrc1)
449  ABI_DEALLOCATE(work1)
450 
451 !Must also symmetrize the electric/magnetic field perturbation response !
452 !(XG 000803 This was not implemented until now)
453  if(sum(ddkfil(:))/=0)then
454 !  Get the symmetry matrices in terms of real space basis
455    ABI_ALLOCATE(symrl1,(3,3,nsym1))
456    do isym=1,nsym1
457      call mati3inv(symrc1(:,:,isym),symrl1(:,:,isym))
458    end do
459 !  There should not be any imaginary part, but stay general (for debugging)
460    d2nl_elfd(:,:)=d2nl(:,:,dtset%natom+2,idir,ipert)
461    do ii=1,3
462      sumelfd(:)=zero
463      summgfd(:)=zero
464      do isym=1,nsym1
465        do jj=1,3
466          if(symrl1(ii,jj,isym)/=0)then
467            if(ddkfil(jj)==0)then
468              blkflg(ii,dtset%natom+2,idir,ipert)=0
469            end if
470          end if
471        end do
472        sumelfd(:)=sumelfd(:)+dble(symrl1(ii,1,isym))*d2nl_elfd(:,1)+&
473 &       dble(symrl1(ii,2,isym))*d2nl_elfd(:,2)+&
474 &       dble(symrl1(ii,3,isym))*d2nl_elfd(:,3)
475        summgfd(:)=summgfd(:)+dble(symrl1(ii,1,isym))*d2nl_mgfd(:,1)+&
476 &       dble(symrl1(ii,2,isym))*d2nl_mgfd(:,2)+&
477 &       dble(symrl1(ii,3,isym))*d2nl_mgfd(:,3)
478      end do
479      d2nl(:,ii,dtset%natom+2,idir,ipert)=sumelfd(:)/dble(nsym1)
480    end do
481 
482    if ((dtset%prtbbb==1).and.(ipert<=dtset%natom)) then
483      do iband = 1,dtset%mband
484        d2nl_elfd(:,:)=d2bbb(:,:,idir,ipert,iband,iband)
485        do ii=1,3
486          sumelfd(:)=zero
487          do isym=1,nsym1
488            sumelfd(:)=sumelfd(:)+dble(symrl1(ii,1,isym))*d2nl_elfd(:,1)+&
489 &           dble(symrl1(ii,2,isym))*d2nl_elfd(:,2)+&
490 &           dble(symrl1(ii,3,isym))*d2nl_elfd(:,3)
491          end do
492          d2bbb(:,ii,idir,ipert,iband,iband)=sumelfd(:)/dble(nsym1)
493        end do
494      end do  !iband
495    end if
496 
497    ABI_DEALLOCATE(symrl1)
498  end if
499 
500 !----------------------------------------------------------------------------
501 !Now, treat the local contribution
502 
503  nfftot=ngfft(1)*ngfft(2)*ngfft(3)
504  ABI_ALLOCATE(vpsp1,(cplex*nfft))
505  if (ipert /= dtset%natom + 1) then
506    n3xccc=0;if(psps%n1xccc/=0) n3xccc=nfft
507    ABI_ALLOCATE(xccc3d1,(cplex*n3xccc))
508    ABI_ALLOCATE(vxc1,(cplex*nfft,nspden))
509 
510    do ipert1=1,mpert
511      do idir1=1,3
512        if(ipert1 <= dtset%natom)then
513 
514 !        Get first-order local potential and first-order pseudo core density
515          call dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir1,ipert1,mpi_enreg,psps%mqgrid_ff,dtset%natom,&
516 &         nattyp,nfft,ngfft,dtset%ntypat,n1,n2,n3,dtset%paral_kgb,ph1d,psps%qgrid_ff,&
517 &         dtset%qptn,ucvol,psps%vlspl,vpsp1,xred)
518          if(psps%n1xccc/=0)then
519            call dfpt_mkcore(cplex,idir1,ipert1,dtset%natom,dtset%ntypat,n1,psps%n1xccc,&
520 &           n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
521          end if
522 
523 !        Get first-order exchange-correlation potential (core-correction contribution only !)
524          if(psps%n1xccc/=0)then
525            option=0
526 !FR SPr EB non-collinear magnetism
527            if (nspden==4.and.present(rhor).and.present(vxc)) then
528              optnc=1
529              call dfpt_mkvxc_noncoll(cplex,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,rhodummy,0,rhodummy,0,rhodummy,0,&
530 &             nkxc,nspden,n3xccc,optnc,option,dtset%paral_kgb,dtset%qptn,rhor,rhor1,&
531 &             rprimd,0,vxc,vxc1,xccc3d1)
532            else
533              call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,rhodummy,0,rhodummy,0,&
534 &             nkxc,nspden,n3xccc,option,dtset%paral_kgb,dtset%qptn,rhodummy,&
535 &             rprimd,0,vxc1,xccc3d1)
536            end if
537          else
538            vxc1(:,:)=zero
539          end if
540 
541 !        Norm-conserving pseudpopotential case:
542 !        Combines density j2 with local potential j1 (vpsp1 and vxc1)
543 !        XG030514 : this is a first possible coding, however, each dotprod contains
544 !        a parallel section (reduction), so it is better to use only one dotprod ...
545 !        call dotprod_vn(cplex,rhor1,dr_psp1,di_psp1,mpi_enreg,nfft,nfftot,1,2,vpsp1,ucvol)
546 !        call dotprod_vn(cplex,rhor1,dr_xc1,di_xc1,mpi_enreg,nfft,nfftot,nspden,2,vxc1,ucvol)
547 !        dotr=dr_psp1+dr_xc1;doti=di_psp1+di_xc1... but then, one needs to overload vxc1
548          do ispden=1,min(nspden,2)
549            do ifft=1,cplex*nfft
550              vxc1(ifft,ispden)=vxc1(ifft,ispden)+vpsp1(ifft)
551            end do
552          end do
553          call dotprod_vn(cplex,rhor1,dotr,doti,nfft,nfftot,nspden,2,vxc1,ucvol)
554 
555 !        MVeithen 021212 : in case ipert = 2, these lines compute the local part
556 !        of the Born effective charges from phonon and electric
557 !        field type perturbations, see eq. 43 of
558 !        X. Gonze and C. Lee, PRB 55, 10355 (1997)
559 !        The minus sign is due to the fact that the effective charges
560 !        are minus the second derivatives of the energy
561          if (ipert == dtset%natom+2) then
562            d2lo(1,idir1,ipert1,idir,ipert)=-dotr
563            d2lo(2,idir1,ipert1,idir,ipert)=-doti
564          else
565            d2lo(1,idir1,ipert1,idir,ipert)=dotr
566            d2lo(2,idir1,ipert1,idir,ipert)=doti
567          end if
568 !        Endif ipert1<=natom
569        end if
570      end do
571    end do
572 
573    ABI_DEALLOCATE(vxc1)
574    ABI_DEALLOCATE(xccc3d1)
575 
576  end if ! ipert /= natom +1
577 
578  ABI_DEALLOCATE(d2bbb_k)
579  ABI_DEALLOCATE(d2nl_k)
580  ABI_DEALLOCATE(kg_k)
581  ABI_DEALLOCATE(kg1_k)
582  ABI_DEALLOCATE(vpsp1)
583  ABI_DEALLOCATE(eig_k)
584  ABI_DEALLOCATE(eig1_k)
585 
586  call timab(101,2,tsec)
587 
588  DBG_EXIT("COLL")
589 
590 end subroutine dfpt_nstdy