TABLE OF CONTENTS


ABINIT/dfptnl_loop [ Functions ]

[ Top ] [ Functions ]

NAME

 dfptnl_loop

FUNCTION

 Loop over the perturbations j1, j2 and j3

COPYRIGHT

 Copyright (C) 2018-2018 ABINIT group (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 .
 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) = array for planewave coefficients of wavefunctions
  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)
  gmet(3,3)=reciprocal space metric tensor in bohr**-2
  gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
  gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double
   that of the basis sphere--appropriate for charge density rho(G),
   Hartree potential, and pseudopotentials
  kg(3,mpw*mkmem)=reduced planewave coordinates
  kxc(nfftf,nkxc)=exchange-correlation kernel
  k3xc(nfftf,nk3xc)=third-order exchange-correlation kernel
  mband = maximum number of bands
  mgfft = maximum single fft dimension
  mkmem = Number of k points treated by this node.
  mk1mem = Number of k points for first-order WF treated by this node.
  mpert =maximum number of ipert
  mpi_enreg=MPI-parallelisation information
  mpw   = maximum number of planewaves in basis sphere (large number)
  natom = number of atoms in unit cell
  nattyp(ntypat)= # atoms of each type.
  nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
  ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90)
  nhat=compensation charge density on fine rectangular grid
  nkpt  = number of k points
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  nk3xc=second dimension of the array k3xc
  nspinor = number of spinorial components of the wavefunctions
  nsppol = number of channels for spin-polarization (1 or 2)
  npwarr(nkpt) = array holding npw for each k point
  occ(mband*nkpt*nsppol) = occupation number for each band and k
  paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
  paw_ij0(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation
  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
  pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  ph1df(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information (fine grid)
  psps <type(pseudopotential_type)> = variables related to pseudopotentials
  rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
       that have to be computed
       1   ->   element has to be computed explicitely
      -1   ->   use symmetry operations to obtain the corresponding element
  rhog(2,nfftf)=array for Fourier transform of GS electron density
  rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3.
  rprimd(3,3)=dimensional primitive translations (bohr)
  ucvol = unit cell volume (bohr^3)
  usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory
  vtrial(nfftf,nspden)=GS Vtrial(r).
  vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree)
  xred(3,natom) = reduced atomic coordinates
  nsym1=number of symmetry elements in space group consistent with perturbation
  indsy1(4,nsym1,natom)=indirect indexing array for atom labels
  symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations
  symrc1(3,3,nsym1)=symmetry operations in reciprocal space

OUTPUT

  blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE
                             (=1 if computed)
  d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor
                                    = \sum_{i=1}^9 d3etot_i
  d3etot_1(2,3,mpert,3,mpert,3,mpert) = 1st term of d3etot
  d3etot_2(2,3,mpert,3,mpert,3,mpert) = 2nd term of d3etot
  d3etot_3(2,3,mpert,3,mpert,3,mpert) = 3rd term of d3etot
  d3etot_4(2,3,mpert,3,mpert,3,mpert) = 4th term of d3etot
  d3etot_5(2,3,mpert,3,mpert,3,mpert) = 5th term of d3etot
  d3etot_6(2,3,mpert,3,mpert,3,mpert) = 6th term of d3etot
  d3etot_7(2,3,mpert,3,mpert,3,mpert) = 7th term of d3etot
  d3etot_8(2,3,mpert,3,mpert,3,mpert) = 8th term of d3etot
  d3etot_9(2,3,mpert,3,mpert,3,mpert) = 9th term of d3etot

SIDE EFFECTS

  hdr <type(hdr_type)>=the header of wf, den and pot files

PARENTS

      nonlinear

CHILDREN

      appdig,dfpt_mkcore,dfpt_mkvxc,dfpt_vlocal,dfptnl_mv,dfptnl_resp
      dotprod_vn,fourdp,getph,hartre,hdr_free,initylmg,inwffil,read_rhor
      status,timab,wffclose,wrtout

SOURCE

143 subroutine dfptnl_loop(atindx,blkflg,cg,dtfil,dtset,d3etot,eigen0,gmet,gprimd,gsqcut, &
144 & hdr,kg,kxc,k3xc,mband,mgfft,mgfftf,mkmem,mk1mem,&
145 & mpert,mpi_enreg,mpw,natom,nattyp,ngfftf,nfftf,nhat,nkpt,nkxc,nk3xc,nspinor,nsppol,&
146 & npwarr,occ,paw_an0,paw_ij0,&
147 & pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
148 & ph1d,ph1df,psps,rfpert,rhog,rhor,rprimd,ucvol,usecprj,vtrial,vxc,xred,&
149 & nsym1,indsy1,symaf1,symrc1,&
150 & d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9)
151 
152  use defs_basis
153  use defs_datatypes
154  use defs_abitypes
155  use defs_wvltypes
156 
157  use m_errors
158  use m_abicore
159  use m_hdr
160  use m_nctk
161  use m_wffile
162  use m_wfk
163 
164  use m_dtfil,       only : status
165  use m_time,        only : timab
166  use m_io_tools,    only : file_exists
167  use m_kg,          only : getph
168  use m_inwffil,     only : inwffil
169  use m_fft,         only : fourdp
170  use m_ioarr,       only : read_rhor
171  use m_hamiltonian, only : destroy_hamiltonian,destroy_rf_hamiltonian,gs_hamiltonian_type,&
172                            init_hamiltonian,init_rf_hamiltonian,rf_hamiltonian_type
173  use m_pawdij,      only : pawdij, pawdijfr, symdij
174  use m_pawfgr,      only : pawfgr_type
175  use m_pawfgrtab,   only : pawfgrtab_type
176  use m_paw_an,      only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags
177  use m_paw_ij,      only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags, paw_ij_print
178  use m_pawang,      only : pawang_type
179  use m_pawrad,      only : pawrad_type
180  use m_pawrhoij,    only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_nullify, pawrhoij_io
181  use m_paw_nhat,    only : pawmknhat,pawnhatfr
182  use m_paw_denpot,  only : pawdenpot
183  use m_pawtab,      only : pawtab_type
184  use m_rf2,         only : rf2_getidir
185  use m_initylmg,    only : initylmg
186  use m_atm2fft,     only : dfpt_atm2fft
187  use m_dfpt_mkvxc,  only : dfpt_mkvxc
188  use m_dfpt_rhotov, only : dfpt_rhotov
189  use m_mkcore,      only : dfpt_mkcore
190  use m_mklocl,      only : dfpt_vlocal
191  use m_dfptnl_pert, only : dfptnl_pert
192 
193 !This section has been created automatically by the script Abilint (TD).
194 !Do not modify the following lines by hand.
195 #undef ABI_FUNC
196 #define ABI_FUNC 'dfptnl_loop'
197 !End of the abilint section
198 
199  implicit none
200 
201 !Arguments ------------------------------------
202 !scalars
203  integer,intent(in) :: mband,mgfft,mgfftf,mk1mem,mkmem,mpert,mpw,natom,nfftf
204  integer,intent(in) :: nk3xc,nkpt,nkxc,nspinor,nsppol,nsym1,usecprj
205  real(dp),intent(in) :: gsqcut,ucvol
206  type(MPI_type),intent(inout) :: mpi_enreg
207  type(datafiles_type),intent(in) :: dtfil
208  type(dataset_type),intent(in) :: dtset
209  type(hdr_type),intent(inout) :: hdr
210  type(pawang_type),intent(inout) :: pawang,pawang1
211  type(pawfgr_type),intent(in) :: pawfgr
212  type(pseudopotential_type),intent(in) :: psps
213 
214 !arrays
215  integer,intent(in) :: atindx(natom),kg(3,mk1mem*mpw)
216  integer,intent(in) :: nattyp(psps%ntypat),ngfftf(18),npwarr(nkpt)
217  integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert)
218  integer,intent(in) :: indsy1(4,nsym1,dtset%natom),symaf1(nsym1),symrc1(3,3,nsym1)
219  integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) !vz_i
220  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3)
221  real(dp),intent(in) :: eigen0(dtset%mband*dtset%nkpt*dtset%nsppol)
222  real(dp),intent(in) :: gprimd(3,3),k3xc(nfftf,nk3xc),kxc(nfftf,nkxc)
223  real(dp),intent(in) :: nhat(nfftf,dtset%nspden)
224  real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),rprimd(3,3)
225  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),ph1df(2,3*(2*mgfftf+1)*natom)
226  real(dp),intent(in) :: vtrial(nfftf,dtset%nspden),xred(3,natom)
227  real(dp),intent(in) :: vxc(nfftf,dtset%nspden)
228  real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
229  real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) !vz_i
230  real(dp),intent(inout) :: d3etot_1(2,3,mpert,3,mpert,3,mpert)
231  real(dp),intent(inout) :: d3etot_2(2,3,mpert,3,mpert,3,mpert)
232  real(dp),intent(inout) :: d3etot_3(2,3,mpert,3,mpert,3,mpert)
233  real(dp),intent(inout) :: d3etot_4(2,3,mpert,3,mpert,3,mpert)
234  real(dp),intent(inout) :: d3etot_5(2,3,mpert,3,mpert,3,mpert)
235  real(dp),intent(inout) :: d3etot_6(2,3,mpert,3,mpert,3,mpert)
236  real(dp),intent(inout) :: d3etot_7(2,3,mpert,3,mpert,3,mpert)
237  real(dp),intent(inout) :: d3etot_8(2,3,mpert,3,mpert,3,mpert)
238  real(dp),intent(inout) :: d3etot_9(2,3,mpert,3,mpert,3,mpert)
239  type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom*psps%usepaw)
240  type(pawrhoij_type),intent(in) :: pawrhoij(natom*psps%usepaw)
241  type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
242  type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
243  type(paw_an_type),intent(in) :: paw_an0(natom*psps%usepaw)
244  type(paw_ij_type),intent(in) :: paw_ij0(natom*psps%usepaw)
245 
246 !Local variables-------------------------------
247 !scalars
248  integer,parameter :: level=51
249  integer :: ask_accurate,comm_cell,counter,cplex,cplex_rhoij,formeig,flag1,flag3
250  integer :: has_dijfr,has_diju
251  integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,iatom,idir_dkde,ierr,iexit,ii,ireadwf
252  integer :: mcg,mpsang,n1,n2,n3,n3xccc,ndir,nfftotf,nhat1grdim,npert_phon,nspden,nspden_rhoij,nwffile
253  integer :: option,optene,optfr,optorth,pert1case,pert2case,pert3case
254  integer :: rdwrpaw,second_idir,timrev,usexcnhat
255  real(dp) :: dummy_real,ecut_eff
256  character(len=500) :: message
257  character(len=fnlen) :: fiden1i,fiwf1i,fiwf2i,fiwf3i,fiwfddk,fnamewff(5)
258  type(gs_hamiltonian_type) :: gs_hamkq
259  type(wffile_type) :: wff1,wff2,wff3,wfft1,wfft2,wfft3
260  type(wfk_t) :: ddk_f(5)
261  type(wvl_data) :: wvl
262  type(hdr_type) :: hdr_den
263 !arrays
264  integer :: file_index(5)
265  real(dp) :: qphon(3),tsec(2)
266  real(dp),allocatable :: cg1(:,:),cg2(:,:),cg3(:,:),eigen1(:),eigen2(:),eigen3(:)
267  real(dp),allocatable :: nhat1_i1pert(:,:),nhat1_i2pert(:,:),nhat1_i3pert(:,:)
268  real(dp),allocatable :: nhat1gr(:,:,:),vresid_dum(:,:)
269  real(dp),allocatable :: rho1r1(:,:)
270  real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:),rho3r1(:,:),vhartr1_i2pert(:)
271  real(dp),allocatable :: vpsp1(:),vxc1_i2pert(:,:),work(:)
272  real(dp),allocatable,target :: vtrial1_i2pert(:,:)
273  real(dp),pointer :: vtrial1_tmp(:,:)
274  real(dp),allocatable :: xccc3d1(:),xccc3d2(:),xccc3d3(:)
275  type(pawrhoij_type),allocatable :: pawrhoij1_i1pert(:),pawrhoij1_i2pert(:),pawrhoij1_i3pert(:)
276  type(paw_an_type),allocatable :: paw_an1_i2pert(:)
277  type(paw_ij_type),allocatable :: paw_ij1_i2pert(:)
278 
279 ! ***********************************************************************
280 
281  DBG_ENTER("COLL")
282 
283  call timab(503,1,tsec)
284  call status(0,dtfil%filstat,iexit,level,'enter         ')
285 
286  comm_cell = mpi_enreg%comm_cell
287 
288  timrev = 1 ! as q=0
289  cplex = 2 - timrev
290  nspden = dtset%nspden
291  ecut_eff = (dtset%ecut)*(dtset%dilatmx)**2
292  mpsang = psps%mpsang
293  optorth=1;if (psps%usepaw==1) optorth=0
294 
295  qphon(:)=zero
296 
297  ABI_ALLOCATE(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
298  ABI_ALLOCATE(cg2,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
299  ABI_ALLOCATE(cg3,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
300  ABI_ALLOCATE(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
301  ABI_ALLOCATE(eigen2,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
302  ABI_ALLOCATE(eigen3,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
303  ABI_ALLOCATE(rho1r1,(cplex*nfftf,dtset%nspden))
304  ABI_ALLOCATE(rho2r1,(cplex*nfftf,dtset%nspden))
305  ABI_ALLOCATE(rho2g1,(2,nfftf))
306  ABI_ALLOCATE(rho3r1,(cplex*nfftf,dtset%nspden))
307 
308  ask_accurate=1 ; formeig = 1 ; ireadwf = 1
309  n1=ngfftf(1) ; n2=ngfftf(2) ; n3=ngfftf(3)
310  nfftotf=n1*n2*n3
311 
312 !==== Initialize most of the Hamiltonian (and derivative) ====
313 !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
314 !2) Perform the setup needed for the non-local factors:
315 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
316 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
317  call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,&
318 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
319 & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda,paw_ij=paw_ij0)
320 
321  ABI_ALLOCATE(vpsp1,(cplex*nfftf))
322  ABI_ALLOCATE(xccc3d1,(cplex*nfftf))
323  ABI_ALLOCATE(xccc3d2,(cplex*nfftf))
324  ABI_ALLOCATE(xccc3d3,(cplex*nfftf))
325  ABI_ALLOCATE(vhartr1_i2pert,(cplex*nfftf))
326  ABI_ALLOCATE(vxc1_i2pert,(cplex*nfftf,dtset%nspden))
327  ABI_ALLOCATE(vtrial1_i2pert,(cplex*nfftf,dtset%nspden))
328 
329  ABI_ALLOCATE(vresid_dum,(0,0))
330 ! PAW stuff
331  usexcnhat = 0
332  nhat1grdim=0
333  ABI_ALLOCATE(nhat1gr,(0,0,0))
334  nhat1gr(:,:,:) = zero
335  rdwrpaw=psps%usepaw
336 !Allocate 1st-order PAW occupancies (rhoij1)
337  if (psps%usepaw==1) then
338    cplex_rhoij=max(cplex,dtset%pawcpxocc);nspden_rhoij=dtset%nspden
339    ABI_DATATYPE_ALLOCATE(pawrhoij1_i1pert,(natom))
340    ABI_DATATYPE_ALLOCATE(pawrhoij1_i2pert,(natom))
341    ABI_DATATYPE_ALLOCATE(pawrhoij1_i3pert,(natom))
342    call pawrhoij_nullify(pawrhoij1_i1pert)
343    call pawrhoij_nullify(pawrhoij1_i2pert)
344    call pawrhoij_nullify(pawrhoij1_i3pert)
345    call pawrhoij_alloc(pawrhoij1_i1pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
346 &   dtset%typat,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
347    call pawrhoij_alloc(pawrhoij1_i2pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
348 &   dtset%typat,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
349    call pawrhoij_alloc(pawrhoij1_i3pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
350 &   dtset%typat,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
351  else
352    ABI_DATATYPE_ALLOCATE(pawrhoij1_i1pert,(0))
353    ABI_DATATYPE_ALLOCATE(pawrhoij1_i2pert,(0))
354    ABI_DATATYPE_ALLOCATE(pawrhoij1_i3pert,(0))
355  end if
356 
357  mcg=mpw*nspinor*mband*mkmem*nsppol
358 
359 !Allocations/initializations for PAW only
360  if(psps%usepaw==1) then
361    usexcnhat=maxval(pawtab(:)%usexcnhat)
362 !  1st-order compensation density
363    ABI_ALLOCATE(nhat1_i1pert,(cplex*nfftf,dtset%nspden))
364    nhat1_i1pert=zero
365    ABI_ALLOCATE(nhat1_i2pert,(cplex*nfftf,dtset%nspden))
366    nhat1_i2pert=zero
367    ABI_ALLOCATE(nhat1_i3pert,(cplex*nfftf,dtset%nspden))
368    nhat1_i3pert=zero
369 
370 !  1st-order arrays/variables related to the PAW spheres
371    ABI_DATATYPE_ALLOCATE(paw_an1_i2pert,(natom))
372    ABI_DATATYPE_ALLOCATE(paw_ij1_i2pert,(natom))
373    call paw_an_nullify(paw_an1_i2pert)
374    call paw_ij_nullify(paw_ij1_i2pert)
375 
376    has_dijfr=1
377    has_diju=0; if(dtset%usepawu==5.or.dtset%usepaw==6) has_diju=1
378 
379    call paw_an_init(paw_an1_i2pert,dtset%natom,dtset%ntypat,0,0,dtset%nspden,cplex,dtset%pawxcdev,&
380 &   dtset%typat,pawang,pawtab,has_vxc=1,&
381 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
382 
383    call paw_ij_init(paw_ij1_i2pert,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,dtset%natom,&
384 &   dtset%ntypat,dtset%typat,pawtab,&
385 &   has_dij=1,has_dijhartree=1,has_dijfr=has_dijfr,has_dijU=has_diju,&
386 &   mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
387  else
388    ABI_ALLOCATE(nhat1_i1pert,(0,0))
389    ABI_ALLOCATE(nhat1_i2pert,(0,0))
390    ABI_ALLOCATE(nhat1_i3pert,(0,0))
391    ABI_DATATYPE_ALLOCATE(paw_an1_i2pert,(0))
392    ABI_DATATYPE_ALLOCATE(paw_ij1_i2pert,(0))
393  end if ! PAW
394 
395  n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf
396 
397 !Loop over the perturbations j1, j2, j3
398 
399  pert1case = 0 ; pert2case = 0 ; pert3case = 0
400 
401  do i1pert = 1, mpert
402    do i1dir = 1, 3
403 
404      if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then
405 
406        pert1case = i1dir + (i1pert-1)*3
407        counter = pert1case
408        call appdig(pert1case,dtfil%fnamewff1,fiwf1i)
409 
410        call status(counter,dtfil%filstat,iexit,level,'call inwffil  ')
411 
412        call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
413 &       formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
414 &       dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
415 &       dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
416 &       dtset%nsppol,dtset%nsym,&
417 &       occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
418 &       dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl)
419 
420        if (ireadwf==1) then
421          call WffClose (wff1,ierr)
422        end if
423 
424        flag1 = 0
425        rho1r1(:,:) = zero
426        if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
427          call appdig(pert1case,dtfil%fildens1in,fiden1i)
428          call status(counter,dtfil%filstat,iexit,level,'call ioarr    ')
429 
430          call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho1r1, &
431          hdr_den, pawrhoij1_i1pert, comm_cell, check_hdr=hdr)
432          call hdr_free(hdr_den)
433        end if
434 
435        xccc3d1(:) = zero
436        if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
437          ndir=1
438          call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i1dir,i1pert,&
439 &         mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,&
440 &         ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
441 &         atmrhor1=xccc3d1,optn_in=n3xccc/nfftf,optn2_in=1,optv_in=0,vspl=psps%vlspl)
442        else
443     !    Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
444     !    ------------------------------------------------------------------------------
445          if(psps%n1xccc/=0)then
446            call dfpt_mkcore(cplex,i1dir,i1pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,&
447 &           n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
448          end if ! psps%n1xccc/=0
449        end if ! usepaw
450 
451        do i3pert = 1, mpert
452          do i3dir = 1, 3
453 
454            if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then
455 
456              pert3case = i3dir + (i3pert-1)*3
457              counter = 100*pert3case + pert1case
458              call appdig(pert3case,dtfil%fnamewff1,fiwf3i)
459 
460              call status(counter,dtfil%filstat,iexit,level,'call inwffil  ')
461              call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen3,dtset%exchn2n3d,&
462 &             formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
463 &             dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
464 &             dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
465 &             dtset%nsppol,dtset%nsym,&
466 &             occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
467 &             dtfil%unkg1,wff3,wfft3,dtfil%unwff3,&
468 &             fiwf3i,wvl)
469              if (ireadwf==1) then
470                call WffClose (wff3,ierr)
471              end if
472 
473              flag3 = 0
474              rho3r1(:,:) = zero
475              if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
476 
477                call appdig(pert3case,dtfil%fildens1in,fiden1i)
478                call status(counter,dtfil%filstat,iexit,level,'call ioarr    ')
479 
480                call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho3r1, &
481                hdr_den, pawrhoij1_i3pert, comm_cell, check_hdr=hdr)
482                call hdr_free(hdr_den)
483              end if
484 
485              xccc3d3(:) = zero
486              if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
487                ndir=1
488                call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i3dir,i3pert,&
489 &               mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,&
490 &               ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
491 &               atmrhor1=xccc3d3,optn_in=n3xccc/nfftf,optn2_in=1,optv_in=0,vspl=psps%vlspl)
492              else
493             !    Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
494             !    ------------------------------------------------------------------------------
495                if(psps%n1xccc/=0)then
496                  call dfpt_mkcore(cplex,i3dir,i3pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,&
497 &                 n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d3,xred)
498                end if ! psps%n1xccc/=0
499              end if ! usepaw
500 
501              do i2pert = 1, mpert
502                do i2dir = 1, 3
503 
504                  if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
505 
506                    blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
507 
508                    npert_phon = 0
509                    if(i1pert<=dtset%natom) npert_phon = npert_phon + 1
510                    if(i2pert<=dtset%natom) npert_phon = npert_phon + 1
511                    if(i3pert<=dtset%natom) npert_phon = npert_phon + 1
512                    if (npert_phon>1) then
513                      MSG_ERROR("dfptnl_loop is available with at most one phonon perturbation. Change your input!")
514                    end if
515 
516                    pert2case = i2dir + (i2pert-1)*3
517                    counter = 100*pert2case + pert2case
518                    call appdig(pert2case,dtfil%fnamewff1,fiwf2i)
519 
520                    call status(counter,dtfil%filstat,iexit,level,'call inwffil  ')
521                    call inwffil(ask_accurate,cg2,dtset,dtset%ecut,ecut_eff,eigen2,dtset%exchn2n3d,&
522 &                   formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
523 &                   dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
524 &                   dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
525 &                   dtset%nsppol,dtset%nsym,&
526 &                   occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
527 &                   dtfil%unkg1,wff2,wfft2,dtfil%unwff2,&
528 &                   fiwf2i,wvl)
529                    if (ireadwf==1) then
530                      call WffClose (wff2,ierr)
531                    end if
532 
533 !                  Read the first-order densities from disk-files
534                    rho2r1(:,:) = zero ; rho2g1(:,:) = zero
535 
536                    if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
537 
538                      call appdig(pert2case,dtfil%fildens1in,fiden1i)
539                      call status(counter,dtfil%filstat,iexit,level,'call ioarr    ')
540 
541                      call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho2r1, &
542                      hdr_den, pawrhoij1_i2pert , comm_cell, check_hdr=hdr)
543                      call hdr_free(hdr_den)
544 
545 !                    Compute up+down rho1(G) by fft
546                      ABI_ALLOCATE(work,(cplex*nfftf))
547                      work(:)=rho2r1(:,1)
548                      call status(counter,dtfil%filstat,iexit,level,'call fourdp   ')
549                      call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
550                      ABI_DEALLOCATE(work)
551 
552                    end if
553 
554                    xccc3d2(:)=zero ; vpsp1(:)=zero
555                    !  PAW: compute Vloc(1) and core(1) together in reciprocal space
556                    !  --------------------------------------------------------------
557                    if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
558                      ndir=1
559                      call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i2dir,i2pert,&
560 &                     mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,&
561 &                     ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
562 &                     atmrhor1=xccc3d2,atmvlocr1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,vspl=psps%vlspl)
563                      !    PAW only: we sometimes have to compute 1st-order compensation density
564                      !    and eventually add it to density from 1st-order WFs
565                      !    ----------------------------------------------------------------------
566                      if (psps%usepaw==1) then
567 
568                        !Force the computation of nhatfr
569                        do iatom=1,dtset%natom
570                          pawfgrtab(iatom)%nhatfr_allocated = 0
571                          pawfgrtab(iatom)%nhatfr = zero
572                        end do
573 
574 !                      This portion of code works only when npert_phon<=1
575                        if (i1pert<=natom.and.usexcnhat==0) then
576                          call pawnhatfr(0,i1dir,i1pert,1,dtset%natom,nspden,psps%ntypat,&
577 &                         pawang,pawfgrtab(i1pert),pawrhoij(i1pert),pawtab,rprimd,&
578 &                         mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
579                        end if
580                        if (i2pert<=natom) then
581                          call pawnhatfr(0,i2dir,i2pert,1,dtset%natom,nspden,psps%ntypat,&
582 &                         pawang,pawfgrtab(i2pert),pawrhoij(i2pert),pawtab,rprimd,&
583 &                         mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
584                        end if
585                        if (i3pert<=natom.and.usexcnhat==0) then
586                          call pawnhatfr(0,i3dir,i3pert,1,dtset%natom,nspden,psps%ntypat,&
587 &                         pawang,pawfgrtab(i3pert),pawrhoij(i3pert),pawtab,rprimd,&
588 &                         mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
589                        end if
590 
591                        if (usexcnhat==0) then
592 
593                          call pawmknhat(dummy_real,cplex,0,i1dir,i1pert,0,gprimd,natom,dtset%natom,&
594 &                         nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i1pert,&
595 &                         pawrhoij1_i1pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,&
596 &                         mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
597                          if (flag1==0) then
598                            rho1r1(:,:) = rho1r1(:,:) - nhat1_i1pert(:,:)
599                            flag1 = 1
600                          end if
601 
602                          call pawmknhat(dummy_real,cplex,0,i3dir,i3pert,0,gprimd,natom,dtset%natom,&
603 &                         nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i3pert,&
604 &                         pawrhoij1_i3pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,&
605 &                         mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
606                          if (flag3==0) then
607                            rho3r1(:,:) = rho3r1(:,:) - nhat1_i3pert(:,:)
608                            flag3 = 1
609                          end if
610 
611                        end if
612 
613                        call pawmknhat(dummy_real,cplex,0,i2dir,i2pert,0,gprimd,natom,dtset%natom,&
614 &                       nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i2pert,&
615 &                       pawrhoij1_i2pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,&
616 &                       mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
617 
618                      end if
619 
620                    else
621 
622                   !    Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
623                   !    ------------------------------------------------------------------------------
624                      if(psps%n1xccc/=0)then
625                        call dfpt_mkcore(cplex,i2dir,i2pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,&
626 &                       n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d2,xred)
627                      end if ! psps%n1xccc/=0
628 
629                      call dfpt_vlocal(atindx,cplex,gmet,gsqcut,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,dtset%natom,&
630 &                     nattyp,nfftf,ngfftf,psps%ntypat,n1,n2,n3,dtset%paral_kgb,ph1df,psps%qgrid_vl,&
631 &                     dtset%qptn,ucvol,psps%vlspl,vpsp1,xred)
632 
633                    end if ! usepaw
634 
635                    call status(counter,dtfil%filstat,iexit,level,'get vtrial1   ')
636                    option=1;optene=0
637                    call dfpt_rhotov(cplex,dummy_real,dummy_real,dummy_real,dummy_real,dummy_real,&
638 &                   gsqcut,i2dir,i2pert,dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,&
639 &                   nhat1_i2pert,nhat1gr,nhat1grdim,nkxc,nspden,n3xccc,optene,option,dtset%paral_kgb,&
640 &                   dtset%qptn,rhog,rho2g1,rhor,rho2r1,rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1_i2pert,&
641 &                   vpsp1,vresid_dum,dummy_real,vtrial1_i2pert,vxc,vxc1_i2pert,xccc3d2,dtset%ixcrot)
642 
643                    if (psps%usepaw==1.and.usexcnhat==0) then
644                      rho2r1(:,:) = rho2r1(:,:) - nhat1_i2pert(:,:)
645                    end if
646 
647                    if (psps%usepaw==1)then
648                      call paw_an_reset_flags(paw_an1_i2pert) ! Force the recomputation of on-site potentials
649                      call paw_ij_reset_flags(paw_ij1_i2pert,all=.true.) ! Force the recomputation of Dij
650                      optfr=0
651                      call pawdijfr(cplex,gprimd,i2dir,i2pert,natom,natom,nfftf,ngfftf,nspden,nsppol,&
652 &                     psps%ntypat,optfr,paw_ij1_i2pert,pawang,pawfgrtab,pawrad,pawtab,qphon,&
653 &                     rprimd,ucvol,vpsp1,vtrial,vxc,xred,&
654 &                     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
655 
656 !                    Computation of "on-site" first-order potentials, first-order densities
657                      option=1
658                      call pawdenpot(dummy_real,dummy_real,dummy_real,i2pert,dtset%ixc,natom,dtset%natom,&
659 &                     nspden,psps%ntypat,dtset%nucdipmom,&
660 &                     0,option,paw_an1_i2pert,paw_an0,paw_ij1_i2pert,pawang,&
661 &                     dtset%pawprtvol,pawrad,pawrhoij1_i2pert,dtset%pawspnorb,pawtab,dtset%pawxcdev,&
662 &                     dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,ucvol,psps%znuclpsp, &
663 &                     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
664                 !    First-order Dij computation
665 !                     call timab(561,1,tsec)
666                      if (has_dijfr>0) then
667                        !vpsp1 contribution to Dij already stored in frozen part of Dij
668                        ABI_ALLOCATE(vtrial1_tmp,(cplex*nfftf,nspden))
669                        vtrial1_tmp=vtrial1_i2pert
670                        do ii=1,min(nspden,2)
671                          vtrial1_tmp(:,ii)=vtrial1_tmp(:,ii)-vpsp1(:)
672                        end do
673                      else
674                        vtrial1_tmp => vtrial1_i2pert
675                      end if
676                      call pawdij(cplex,dtset%enunit,gprimd,i2pert,natom,dtset%natom,&
677 &                     nfftf,nfftotf,dtset%nspden,psps%ntypat,paw_an1_i2pert,paw_ij1_i2pert,pawang,&
678 &                     pawfgrtab,dtset%pawprtvol,pawrad,pawrhoij1_i2pert,dtset%pawspnorb,pawtab,&
679 &                     dtset%pawxcdev,qphon,dtset%spnorbscl,ucvol,dtset%charge,vtrial1_tmp,vxc1_i2pert,xred,&
680 &                     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
681                      if (has_dijfr>0) then
682                        ABI_DEALLOCATE(vtrial1_tmp)
683                      end if
684                      call symdij(gprimd,indsy1,i2pert,natom,dtset%natom,nsym1,psps%ntypat,0,&
685 &                     paw_ij1_i2pert,pawang1,dtset%pawprtvol,pawtab,rprimd,symaf1,symrc1, &
686 &                     mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,&
687 &                     qphon=qphon)
688 !                     call timab(561,2,tsec)
689 
690                    end if ! end usepaw section
691 
692                    nwffile = 1
693                    file_index(1) = i2dir + 3*(i2pert-1)
694                    fnamewff(1) = dtfil%fnamewff1
695 
696                    if (i2pert==natom+2) then
697 
698                      nwffile = 3
699                      file_index(2) = i2dir+natom*3
700                      fnamewff(2) = dtfil%fnamewffddk
701 !                    As npert_phon<=1 and i2pert==natom+2, i1pert or i3pert is necessarly equal to natom+2
702                      if (i3pert==natom+2) then
703                        second_idir = i3dir
704                      else if (i1pert==natom+2) then
705                        second_idir = i1dir
706                      else
707                        MSG_BUG(" i1pert or i3pert is supposed to be equal to natom+2, which is not the case here.")
708                      end if
709                      call rf2_getidir(i2dir,second_idir,idir_dkde)
710                      file_index(3) = idir_dkde+9+(dtset%natom+6)*3
711                      fnamewff(3) = dtfil%fnamewffdkde
712 
713                      if (npert_phon==1.and.psps%usepaw==1.and.second_idir/=i2dir) then
714                        nwffile = 5
715                        file_index(4) = second_idir+natom*3
716                        fnamewff(4) = dtfil%fnamewffddk
717                        call rf2_getidir(second_idir,i2dir,idir_dkde) ! i2dir and second_idir are reversed
718                        file_index(5) = idir_dkde+9+(dtset%natom+6)*3
719                        fnamewff(5) = dtfil%fnamewffdkde
720                      end if
721 
722                    end if
723 
724                    do ii=1,nwffile
725                      call appdig(file_index(ii),fnamewff(ii),fiwfddk)
726                      ! Checking the existence of data file
727                      if (.not. file_exists(fiwfddk)) then
728                        ! Trick needed to run Abinit test suite in netcdf mode.
729                        if (file_exists(nctk_ncify(fiwfddk))) then
730                          write(message,"(3a)")"- File: ",trim(fiwfddk),&
731                          " does not exist but found netcdf file with similar name."
732                          call wrtout(std_out,message,'COLL')
733                          fiwfddk = nctk_ncify(fiwfddk)
734                        end if
735                        if (.not. file_exists(fiwfddk)) then
736                          MSG_ERROR('Missing file: '//TRIM(fiwfddk))
737                        end if
738                      end if
739                      write(message,'(2a)')'-dfptnl_loop : read the wavefunctions from file: ',trim(fiwfddk)
740                      call wrtout(std_out,message,'COLL')
741                      call wrtout(ab_out,message,'COLL')
742 !                    Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50)
743                      call wfk_open_read(ddk_f(ii),fiwfddk,1,dtset%iomode,dtfil%unddk+(ii-1),mpi_enreg%comm_cell)
744                    end do
745 
746 !                  Perform DFPT part of the 3dte calculation
747                    call timab(513,1,tsec)
748                    call status(counter,dtfil%filstat,iexit,level,'call dfptnl_resp ')
749 !                  NOTE : eigen2 equals zero here
750 
751                    call dfptnl_pert(atindx,cg,cg1,cg2,cg3,cplex,dtfil,dtset,d3etot,eigen0,gs_hamkq,k3xc,indsy1,i1dir,&
752 &                   i2dir,i3dir,i1pert,i2pert,i3pert,kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,&
753 &                   mpsang,mpw,natom,nattyp,nfftf,nfftotf,ngfftf,nkpt,nk3xc,nspden,nspinor,nsppol,nsym1,npwarr,occ,&
754 &                   pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawtab,pawrhoij,pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,&
755 &                   paw_an0,paw_an1_i2pert,paw_ij1_i2pert,ph1d,psps,rho1r1,rho2r1,rho3r1,&
756 &                   rprimd,symaf1,symrc1,ucvol,vtrial,vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert,&
757 &                   ddk_f,xccc3d1,xccc3d2,xccc3d3,xred,&
758 &                   d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9)
759                    call timab(513,2,tsec)
760 
761                    call status(counter,dtfil%filstat,iexit,level,'after dfptnl_resp')
762 
763 !                  Eventually close the dot file
764                    do ii=1,nwffile
765                      call wfk_close(ddk_f(ii))
766                    end do
767 
768 !                   if (psps%usepaw==1) then
769 !                     do ii=1,natom
770 !                       pawfgrtab(ii)%nhatfr = zero
771 !                     end do
772 !                   end if
773 
774                  end if   ! rfpert
775                end do    ! i2dir
776              end do     ! i2pert
777 
778            end if   ! rfpert
779          end do    ! i3dir
780        end do     ! i3pert
781 
782      end if   ! rfpert
783    end do    ! i1dir
784  end do     ! i1pert
785 
786  call status(0,dtfil%filstat,iexit,level,'exit          ')
787 
788 !More memory cleaning
789  call destroy_hamiltonian(gs_hamkq)
790 
791  ABI_DEALLOCATE(cg1)
792  ABI_DEALLOCATE(cg2)
793  ABI_DEALLOCATE(cg3)
794  ABI_DEALLOCATE(eigen1)
795  ABI_DEALLOCATE(eigen2)
796  ABI_DEALLOCATE(eigen3)
797  ABI_DEALLOCATE(rho1r1)
798  ABI_DEALLOCATE(rho2r1)
799  ABI_DEALLOCATE(rho2g1)
800  ABI_DEALLOCATE(rho3r1)
801  ABI_DEALLOCATE(nhat1gr)
802  ABI_DEALLOCATE(vresid_dum)
803  ABI_DEALLOCATE(vtrial1_i2pert)
804  ABI_DEALLOCATE(vxc1_i2pert)
805  ABI_DEALLOCATE(vhartr1_i2pert)
806  ABI_DEALLOCATE(vpsp1)
807  ABI_DEALLOCATE(xccc3d1)
808  ABI_DEALLOCATE(xccc3d2)
809  ABI_DEALLOCATE(xccc3d3)
810  if (psps%usepaw==1) then
811    call pawrhoij_free(pawrhoij1_i1pert)
812    call pawrhoij_free(pawrhoij1_i2pert)
813    call pawrhoij_free(pawrhoij1_i3pert)
814    ABI_DEALLOCATE(nhat1_i1pert)
815    ABI_DEALLOCATE(nhat1_i2pert)
816    ABI_DEALLOCATE(nhat1_i3pert)
817    call paw_an_free(paw_an1_i2pert)
818    call paw_ij_free(paw_ij1_i2pert)
819    ABI_DATATYPE_DEALLOCATE(paw_an1_i2pert)
820    ABI_DATATYPE_DEALLOCATE(paw_ij1_i2pert)
821  end if
822  ABI_DATATYPE_DEALLOCATE(pawrhoij1_i1pert)
823  ABI_DATATYPE_DEALLOCATE(pawrhoij1_i2pert)
824  ABI_DATATYPE_DEALLOCATE(pawrhoij1_i3pert)
825 
826  call timab(503,2,tsec)
827 
828  DBG_EXIT("COLL")
829 
830 end subroutine dfptnl_loop

ABINIT/m_dfptnl_loop [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfptnl_loop

FUNCTION

COPYRIGHT

  Copyright (C) 2018-2018 ABINIT group (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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_dfptnl_loop
27 
28  implicit none
29 
30  private