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-2022 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

SOURCE

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

ABINIT/m_dfptnl_loop [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfptnl_loop

FUNCTION

COPYRIGHT

  Copyright (C) 2018-2022 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 .

SOURCE

15 #if defined HAVE_CONFIG_H
16 #include "config.h"
17 #endif
18 
19 #include "abi_common.h"
20 
21 module m_dfptnl_loop
22 
23  implicit none
24 
25  private