Line data Source code
1 : !!****m* ABINIT/m_dfptnl_loop
2 : !! NAME
3 : !! m_dfptnl_loop
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2018-2026 ABINIT group (LB)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
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
26 : !!***
27 :
28 : public :: dfptnl_loop
29 : !!***
30 :
31 : contains
32 : !!***
33 :
34 : !!****f* ABINIT/dfptnl_loop
35 : !! NAME
36 : !! dfptnl_loop
37 : !!
38 : !! FUNCTION
39 : !! Loop over the perturbations j1, j2 and j3
40 : !!
41 : !! COPYRIGHT
42 : !! Copyright (C) 2018-2026 ABINIT group (LB)
43 : !! This file is distributed under the terms of the
44 : !! GNU General Public License, see ~abinit/COPYING
45 : !! or http://www.gnu.org/copyleft/gpl.txt .
46 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
47 : !!
48 : !! INPUTS
49 : !! atindx(natom)=index table for atoms (see gstate.f)
50 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions
51 : !! dtfil <type(datafiles_type)>=variables related to files
52 : !! dtset <type(dataset_type)>=all input variables for this dataset
53 : !! eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
54 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2
55 : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
56 : !! gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double
57 : !! that of the basis sphere--appropriate for charge density rho(G),
58 : !! Hartree potential, and pseudopotentials
59 : !! kg(3,mpw*mkmem)=reduced planewave coordinates
60 : !! kxc(nfftf,nkxc)=exchange-correlation kernel
61 : !! k3xc(nfftf,nk3xc)=third-order exchange-correlation kernel
62 : !! mband = maximum number of bands
63 : !! mgfft = maximum single fft dimension
64 : !! mkmem = Number of k points treated by this node.
65 : !! mk1mem = Number of k points for first-order WF treated by this node.
66 : !! mpert =maximum number of ipert
67 : !! mpi_enreg=MPI-parallelisation information
68 : !! mpw = maximum number of planewaves in basis sphere (large number)
69 : !! natom = number of atoms in unit cell
70 : !! nattyp(ntypat)= # atoms of each type.
71 : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
72 : !! ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90)
73 : !! nhat=compensation charge density on fine rectangular grid
74 : !! nkpt = number of k points
75 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
76 : !! nk3xc=second dimension of the array k3xc
77 : !! nspinor = number of spinorial components of the wavefunctions
78 : !! nsppol = number of channels for spin-polarization (1 or 2)
79 : !! npwarr(nkpt) = array holding npw for each k point
80 : !! occ(mband*nkpt*nsppol) = occupation number for each band and k
81 : !! paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
82 : !! paw_ij0(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
83 : !! pawang <type(pawang_type)>=paw angular mesh and related data
84 : !! pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation
85 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
86 : !! pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
87 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
88 : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
89 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
90 : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
91 : !! ph1df(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information (fine grid)
92 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
93 : !! rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
94 : !! that have to be computed
95 : !! 1 -> element has to be computed explicitely
96 : !! -1 -> use symmetry operations to obtain the corresponding element
97 : !! rhog(2,nfftf)=array for Fourier transform of GS electron density
98 : !! rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3.
99 : !! rprimd(3,3)=dimensional primitive translations (bohr)
100 : !! ucvol = unit cell volume (bohr^3)
101 : !! usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory
102 : !! vtrial(nfftf,nspden)=GS Vtrial(r).
103 : !! vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree)
104 : !! xred(3,natom) = reduced atomic coordinates
105 : !! nsym1=number of symmetry elements in space group consistent with perturbation
106 : !! indsy1(4,nsym1,natom)=indirect indexing array for atom labels
107 : !! symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations
108 : !! symrc1(3,3,nsym1)=symmetry operations in reciprocal space
109 : !!
110 : !! OUTPUT
111 : !! blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE
112 : !! (=1 if computed)
113 : !! d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor
114 : !! = \sum_{i=1}^9 d3etot_i
115 : !! d3etot_1(2,3,mpert,3,mpert,3,mpert) = 1st term of d3etot
116 : !! d3etot_2(2,3,mpert,3,mpert,3,mpert) = 2nd term of d3etot
117 : !! d3etot_3(2,3,mpert,3,mpert,3,mpert) = 3rd term of d3etot
118 : !! d3etot_4(2,3,mpert,3,mpert,3,mpert) = 4th term of d3etot
119 : !! d3etot_5(2,3,mpert,3,mpert,3,mpert) = 5th term of d3etot
120 : !! d3etot_6(2,3,mpert,3,mpert,3,mpert) = 6th term of d3etot
121 : !! d3etot_7(2,3,mpert,3,mpert,3,mpert) = 7th term of d3etot
122 : !! d3etot_8(2,3,mpert,3,mpert,3,mpert) = 8th term of d3etot
123 : !! d3etot_9(2,3,mpert,3,mpert,3,mpert) = 9th term of d3etot
124 : !!
125 : !! SIDE EFFECTS
126 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
127 : !!
128 : !! SOURCE
129 :
130 7 : subroutine dfptnl_loop(atindx,blkflg,cg,dtfil,dtset,d3etot,eigen0,gmet,gprimd,gsqcut, &
131 7 : & hdr,kg,kxc,k3xc,mband,mgfft,mgfftf,mkmem,mk1mem,&
132 7 : & mpert,mpi_enreg,mpw,natom,nattyp,ngfftf,nfftf,nhat,nkpt,nkxc,nk3xc,nspinor,nsppol,&
133 7 : & npwarr,occ,paw_an0,paw_ij0,&
134 7 : & pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
135 7 : & ph1d,ph1df,psps,rfpert,rhog,rhor,rprimd,ucvol,usecprj,vtrial,vxc,xred,&
136 7 : & nsym1,indsy1,symaf1,symrc1,&
137 7 : & 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
159 : use m_pawdij, only : pawdij, pawdijfr, symdij
160 : use m_paw_energies,only : paw_energies_type
161 : use m_pawfgr, only : pawfgr_type
162 : use m_pawfgrtab, only : pawfgrtab_type
163 : use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags
164 : use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags, paw_ij_print
165 : use m_pawang, only : pawang_type
166 : use m_pawrad, only : pawrad_type
167 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_nullify, &
168 : & pawrhoij_io, pawrhoij_inquire_dim
169 : use m_paw_nhat, only : pawmknhat,pawnhatfr
170 : use m_paw_denpot, only : pawdenpot
171 : use m_pawtab, only : pawtab_type
172 : use m_rf2, only : rf2_getidir
173 : use m_initylmg, only : initylmg
174 : use m_atm2fft, only : dfpt_atm2fft
175 : use m_dfpt_mkvxc, only : dfpt_mkvxc
176 : use m_dfpt_rhotov, only : dfpt_rhotov
177 : use m_mkcore, only : dfpt_mkcore
178 : use m_mklocl, only : dfpt_vlocal
179 : use m_dfptnl_pert, only : dfptnl_pert
180 :
181 : !Arguments ------------------------------------
182 : !scalars
183 : integer,intent(in) :: mband,mgfft,mgfftf,mk1mem,mkmem,mpert,mpw,natom,nfftf
184 : integer,intent(in) :: nk3xc,nkpt,nkxc,nspinor,nsppol,nsym1,usecprj
185 : real(dp),intent(in) :: gsqcut,ucvol
186 : type(MPI_type),intent(inout) :: mpi_enreg
187 : type(datafiles_type),intent(in) :: dtfil
188 : type(dataset_type),intent(in) :: dtset
189 : type(hdr_type),intent(inout) :: hdr
190 : type(pawang_type),intent(inout) :: pawang,pawang1
191 : type(pawfgr_type),intent(in) :: pawfgr
192 : type(pseudopotential_type),intent(in) :: psps
193 :
194 : !arrays
195 : integer,intent(in) :: atindx(natom),kg(3,mk1mem*mpw)
196 : integer,intent(in) :: nattyp(psps%ntypat),ngfftf(18),npwarr(nkpt)
197 : integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert)
198 : integer,intent(in) :: indsy1(4,nsym1,dtset%natom),symaf1(nsym1),symrc1(3,3,nsym1)
199 : integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) !vz_i
200 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3)
201 : real(dp),intent(in) :: eigen0(dtset%mband*dtset%nkpt*dtset%nsppol)
202 : real(dp),intent(in) :: gprimd(3,3),k3xc(nfftf,nk3xc),kxc(nfftf,nkxc)
203 : real(dp),intent(in) :: nhat(nfftf,dtset%nspden)
204 : real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden),rprimd(3,3)
205 : real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),ph1df(2,3*(2*mgfftf+1)*natom)
206 : real(dp),intent(in) :: vtrial(nfftf,dtset%nspden),xred(3,natom)
207 : real(dp),intent(in) :: vxc(nfftf,dtset%nspden)
208 : real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
209 : real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert) !vz_i
210 : real(dp),intent(inout) :: d3etot_1(2,3,mpert,3,mpert,3,mpert)
211 : real(dp),intent(inout) :: d3etot_2(2,3,mpert,3,mpert,3,mpert)
212 : real(dp),intent(inout) :: d3etot_3(2,3,mpert,3,mpert,3,mpert)
213 : real(dp),intent(inout) :: d3etot_4(2,3,mpert,3,mpert,3,mpert)
214 : real(dp),intent(inout) :: d3etot_5(2,3,mpert,3,mpert,3,mpert)
215 : real(dp),intent(inout) :: d3etot_6(2,3,mpert,3,mpert,3,mpert)
216 : real(dp),intent(inout) :: d3etot_7(2,3,mpert,3,mpert,3,mpert)
217 : real(dp),intent(inout) :: d3etot_8(2,3,mpert,3,mpert,3,mpert)
218 : real(dp),intent(inout) :: d3etot_9(2,3,mpert,3,mpert,3,mpert)
219 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(natom*psps%usepaw)
220 : type(pawrhoij_type),intent(in) :: pawrhoij(natom*psps%usepaw)
221 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
222 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
223 : type(paw_an_type),intent(in) :: paw_an0(natom*psps%usepaw)
224 : type(paw_ij_type),intent(in) :: paw_ij0(natom*psps%usepaw)
225 :
226 : !Local variables-------------------------------
227 : !scalars
228 : integer,parameter :: level=51
229 : integer :: ask_accurate,comm_cell,counter,cplex,cplex_rhoij,formeig,flag1,flag3
230 : integer :: has_dijfr,has_diju
231 : integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,iatom,idir_dkde,ierr,ii,ireadwf
232 : integer :: mcg,mpsang,n1,n2,n3,n3xccc,ndir,nfftotf,nhat1grdim,npert_phon,nspden,nspden_rhoij,nwffile
233 : integer :: option,optene,optfr,optorth,pert1case,pert2case,pert3case
234 : integer :: qphase_rhoij,rdwrpaw,second_idir,timrev,usexcnhat
235 : logical :: non_magnetic_xc
236 : real(dp) :: dummy_real,ecut_eff,el_temp
237 : character(len=500) :: message
238 : character(len=fnlen) :: fiden1i,fiwf1i,fiwf2i,fiwf3i,fiwfddk,fnamewff(5)
239 7 : type(gs_hamiltonian_type) :: gs_hamkq
240 : type(paw_energies_type) :: paw_energies_dum
241 : type(wffile_type) :: wff1,wff2,wff3,wfft1,wfft2,wfft3
242 1120 : type(wfk_t) :: ddk_f(5)
243 7 : type(wvl_data) :: wvl
244 7 : type(hdr_type) :: hdr_den
245 : !arrays
246 : integer :: file_index(5)
247 : real(dp) :: qphon(3),tsec(2)
248 7 : real(dp),allocatable :: cg1(:,:),cg2(:,:),cg3(:,:),eigen1(:),eigen2(:),eigen3(:)
249 7 : real(dp),allocatable :: nhat1_i1pert(:,:),nhat1_i2pert(:,:),nhat1_i3pert(:,:)
250 7 : real(dp),allocatable :: nhat1gr(:,:,:),vresid_dum(:,:)
251 7 : real(dp),allocatable :: rho1r1(:,:)
252 7 : real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:),rho3r1(:,:),vhartr1_i2pert(:)
253 7 : real(dp),allocatable :: vpsp1(:),vxc1_i2pert(:,:),work(:)
254 7 : real(dp),allocatable,target :: vtrial1_i2pert(:,:)
255 7 : real(dp),pointer :: vtrial1_tmp(:,:)
256 7 : real(dp),allocatable :: xccc3d1(:),xccc3d2(:),xccc3d3(:)
257 7 : type(pawrhoij_type),allocatable :: pawrhoij1_i1pert(:),pawrhoij1_i2pert(:),pawrhoij1_i3pert(:)
258 7 : type(paw_an_type),allocatable :: paw_an1_i2pert(:)
259 7 : type(paw_ij_type),allocatable :: paw_ij1_i2pert(:)
260 :
261 : ! ***********************************************************************
262 :
263 : DBG_ENTER("COLL")
264 :
265 7 : call timab(503,1,tsec)
266 :
267 7 : comm_cell = mpi_enreg%comm_cell
268 :
269 7 : timrev = 1 ! as q=0
270 7 : cplex = 2 - timrev
271 7 : nspden = dtset%nspden
272 7 : ecut_eff = (dtset%ecut)*(dtset%dilatmx)**2
273 7 : mpsang = psps%mpsang
274 7 : optorth=1;if (psps%usepaw==1) optorth=0
275 :
276 7 : qphon(:)=zero
277 :
278 21 : ABI_MALLOC(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
279 14 : ABI_MALLOC(cg2,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
280 14 : ABI_MALLOC(cg3,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
281 21 : ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
282 14 : ABI_MALLOC(eigen2,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
283 14 : ABI_MALLOC(eigen3,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
284 28 : ABI_MALLOC(rho1r1,(cplex*nfftf,dtset%nspden))
285 21 : ABI_MALLOC(rho2r1,(cplex*nfftf,dtset%nspden))
286 21 : ABI_MALLOC(rho2g1,(2,nfftf))
287 21 : ABI_MALLOC(rho3r1,(cplex*nfftf,dtset%nspden))
288 :
289 7 : ask_accurate=1 ; formeig = 1 ; ireadwf = 1
290 7 : n1=ngfftf(1) ; n2=ngfftf(2) ; n3=ngfftf(3)
291 7 : nfftotf=n1*n2*n3
292 :
293 : !==== Initialize most of the Hamiltonian (and derivative) ====
294 : !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
295 : !2) Perform the setup needed for the non-local factors:
296 : !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
297 : !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
298 : call gs_hamkq%init(psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,&
299 : & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
300 : & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,paw_ij=paw_ij0,&
301 7 : & gpu_option=dtset%gpu_option)
302 :
303 21 : ABI_MALLOC(vpsp1,(cplex*nfftf))
304 14 : ABI_MALLOC(xccc3d1,(cplex*nfftf))
305 14 : ABI_MALLOC(xccc3d2,(cplex*nfftf))
306 14 : ABI_MALLOC(xccc3d3,(cplex*nfftf))
307 14 : ABI_MALLOC(vhartr1_i2pert,(cplex*nfftf))
308 21 : ABI_MALLOC(vxc1_i2pert,(cplex*nfftf,dtset%nspden))
309 21 : ABI_MALLOC(vtrial1_i2pert,(cplex*nfftf,dtset%nspden))
310 :
311 7 : ABI_MALLOC(vresid_dum,(0,0))
312 : ! PAW stuff
313 7 : usexcnhat = 0
314 7 : nhat1grdim=0
315 7 : ABI_MALLOC(nhat1gr,(0,0,0))
316 : nhat1gr(:,:,:) = zero
317 7 : rdwrpaw=psps%usepaw
318 : !Allocate 1st-order PAW occupancies (rhoij1)
319 7 : if (psps%usepaw==1) then
320 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
321 4 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc)
322 20 : ABI_MALLOC(pawrhoij1_i1pert,(natom))
323 16 : ABI_MALLOC(pawrhoij1_i2pert,(natom))
324 16 : ABI_MALLOC(pawrhoij1_i3pert,(natom))
325 4 : call pawrhoij_nullify(pawrhoij1_i1pert)
326 4 : call pawrhoij_nullify(pawrhoij1_i2pert)
327 4 : call pawrhoij_nullify(pawrhoij1_i3pert)
328 : call pawrhoij_alloc(pawrhoij1_i1pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
329 4 : & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
330 : call pawrhoij_alloc(pawrhoij1_i2pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
331 4 : & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
332 : call pawrhoij_alloc(pawrhoij1_i3pert,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
333 8 : & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
334 : else
335 3 : ABI_MALLOC(pawrhoij1_i1pert,(0))
336 3 : ABI_MALLOC(pawrhoij1_i2pert,(0))
337 3 : ABI_MALLOC(pawrhoij1_i3pert,(0))
338 : end if
339 :
340 7 : mcg=mpw*nspinor*mband*mkmem*nsppol
341 :
342 : !Get electronic temperature from dtset
343 7 : el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9)
344 :
345 : !Allocations/initializations for PAW only
346 7 : if(psps%usepaw==1) then
347 12 : usexcnhat=maxval(pawtab(:)%usexcnhat)
348 : ! 1st-order compensation density
349 12 : ABI_MALLOC(nhat1_i1pert,(cplex*nfftf,dtset%nspden))
350 40835 : nhat1_i1pert=zero
351 16 : ABI_MALLOC(nhat1_i2pert,(cplex*nfftf,dtset%nspden))
352 40835 : nhat1_i2pert=zero
353 12 : ABI_MALLOC(nhat1_i3pert,(cplex*nfftf,dtset%nspden))
354 40835 : nhat1_i3pert=zero
355 :
356 : ! 1st-order arrays/variables related to the PAW spheres
357 20 : ABI_MALLOC(paw_an1_i2pert,(natom))
358 20 : ABI_MALLOC(paw_ij1_i2pert,(natom))
359 4 : call paw_an_nullify(paw_an1_i2pert)
360 4 : call paw_ij_nullify(paw_ij1_i2pert)
361 :
362 4 : has_dijfr=1
363 4 : has_diju=merge(0,1,dtset%usepawu==0)
364 :
365 : call paw_an_init(paw_an1_i2pert,dtset%natom,dtset%ntypat,0,0,dtset%nspden,cplex,dtset%pawxcdev,&
366 : & dtset%typat,pawang,pawtab,has_vxc=1,&
367 4 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
368 :
369 : call paw_ij_init(paw_ij1_i2pert,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,dtset%natom,&
370 : & dtset%ntypat,dtset%typat,pawtab,&
371 : & has_dij=1,has_dijhartree=1,has_dijfr=has_dijfr,has_dijU=has_diju,&
372 4 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
373 : else
374 3 : ABI_MALLOC(nhat1_i1pert,(0,0))
375 3 : ABI_MALLOC(nhat1_i2pert,(0,0))
376 3 : ABI_MALLOC(nhat1_i3pert,(0,0))
377 3 : ABI_MALLOC(paw_an1_i2pert,(0))
378 3 : ABI_MALLOC(paw_ij1_i2pert,(0))
379 : end if ! PAW
380 :
381 7 : n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf
382 7 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
383 :
384 : !Loop over the perturbations j1, j2, j3
385 :
386 7 : pert1case = 0 ; pert2case = 0 ; pert3case = 0
387 :
388 63 : do i1pert = 1, mpert
389 231 : do i1dir = 1, 3
390 :
391 134624 : if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then
392 :
393 33 : pert1case = i1dir + (i1pert-1)*3
394 33 : counter = pert1case
395 33 : call appdig(pert1case,dtfil%fnamewff1,fiwf1i)
396 :
397 : call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
398 : & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
399 : & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
400 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
401 : & dtset%nsppol,dtset%nsym,&
402 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
403 33 : & dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl)
404 :
405 33 : if (ireadwf==1) then
406 33 : call WffClose (wff1,ierr)
407 : end if
408 :
409 90 : flag1 = 0
410 283266 : rho1r1(:,:) = zero
411 33 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
412 33 : call appdig(pert1case,dtfil%fildens1in,fiden1i)
413 :
414 : call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho1r1, &
415 33 : hdr_den, pawrhoij1_i1pert, comm_cell, check_hdr=hdr)
416 33 : call hdr_den%free()
417 : end if
418 :
419 155553 : xccc3d1(:) = zero
420 33 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
421 24 : ndir=1
422 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i1dir,i1pert,&
423 : & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,&
424 : & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
425 24 : & atmrhor1=xccc3d1,optn_in=n3xccc/nfftf,optn2_in=1,optv_in=0,vspl=psps%vlspl)
426 : else
427 : ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
428 : ! ------------------------------------------------------------------------------
429 9 : if(psps%n1xccc/=0)then
430 : call dfpt_mkcore(cplex,i1dir,i1pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,&
431 9 : & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
432 : end if ! psps%n1xccc/=0
433 : end if ! usepaw
434 :
435 297 : do i3pert = 1, mpert
436 1089 : do i3dir = 1, 3
437 :
438 26400 : if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then
439 :
440 83 : pert3case = i3dir + (i3pert-1)*3
441 83 : counter = 100*pert3case + pert1case
442 83 : call appdig(pert3case,dtfil%fnamewff1,fiwf3i)
443 :
444 : call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen3,dtset%exchn2n3d,&
445 : & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
446 : & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
447 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
448 : & dtset%nsppol,dtset%nsym,&
449 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
450 : & dtfil%unkg1,wff3,wfft3,dtfil%unwff3,&
451 83 : & fiwf3i,wvl)
452 83 : if (ireadwf==1) then
453 83 : call WffClose (wff3,ierr)
454 : end if
455 :
456 228 : flag3 = 0
457 730956 : rho3r1(:,:) = zero
458 83 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
459 :
460 83 : call appdig(pert3case,dtfil%fildens1in,fiden1i)
461 :
462 : call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho3r1, &
463 83 : hdr_den, pawrhoij1_i3pert, comm_cell, check_hdr=hdr)
464 83 : call hdr_den%free()
465 : end if
466 :
467 397955 : xccc3d3(:) = zero
468 83 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
469 62 : ndir=1
470 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i3dir,i3pert,&
471 : & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,&
472 : & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
473 62 : & atmrhor1=xccc3d3,optn_in=n3xccc/nfftf,optn2_in=1,optv_in=0,vspl=psps%vlspl)
474 : else
475 : ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
476 : ! ------------------------------------------------------------------------------
477 21 : if(psps%n1xccc/=0)then
478 : call dfpt_mkcore(cplex,i3dir,i3pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,&
479 21 : & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d3,xred)
480 : end if ! psps%n1xccc/=0
481 : end if ! usepaw
482 :
483 747 : do i2pert = 1, mpert
484 2739 : do i2dir = 1, 3
485 :
486 2656 : if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
487 :
488 291 : blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
489 :
490 291 : npert_phon = 0
491 291 : if(i1pert<=dtset%natom) npert_phon = npert_phon + 1
492 291 : if(i2pert<=dtset%natom) npert_phon = npert_phon + 1
493 291 : if(i3pert<=dtset%natom) npert_phon = npert_phon + 1
494 291 : if (npert_phon>1) then
495 0 : ABI_ERROR("dfptnl_loop is available with at most one phonon perturbation. Change your input!")
496 : end if
497 :
498 291 : pert2case = i2dir + (i2pert-1)*3
499 291 : counter = 100*pert2case + pert2case
500 291 : call appdig(pert2case,dtfil%fnamewff1,fiwf2i)
501 :
502 : call inwffil(ask_accurate,cg2,dtset,dtset%ecut,ecut_eff,eigen2,dtset%exchn2n3d,&
503 : & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
504 : & dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
505 : & dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
506 : & dtset%nsppol,dtset%nsym,&
507 : & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
508 : & dtfil%unkg1,wff2,wfft2,dtfil%unwff2,&
509 291 : & fiwf2i,wvl)
510 291 : if (ireadwf==1) then
511 291 : call WffClose (wff2,ierr)
512 : end if
513 :
514 : ! Read the first-order densities from disk-files
515 7423668 : rho2r1(:,:) = zero ; rho2g1(:,:) = zero
516 :
517 291 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
518 :
519 291 : call appdig(pert2case,dtfil%fildens1in,fiden1i)
520 :
521 : call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rho2r1, &
522 291 : hdr_den, pawrhoij1_i2pert , comm_cell, check_hdr=hdr)
523 291 : call hdr_den%free()
524 :
525 : ! Compute up+down rho1(G) by fft
526 873 : ABI_MALLOC(work,(cplex*nfftf))
527 1512723 : work(:)=rho2r1(:,1)
528 291 : call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfftf,1,ngfftf,0)
529 291 : ABI_FREE(work)
530 :
531 : end if
532 :
533 3025155 : xccc3d2(:)=zero ; vpsp1(:)=zero
534 : ! PAW: compute Vloc(1) and core(1) together in reciprocal space
535 : ! --------------------------------------------------------------
536 291 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
537 246 : ndir=1
538 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,i2dir,i2pert,&
539 : & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,psps%ntypat,&
540 : & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
541 246 : & atmrhor1=xccc3d2,atmvlocr1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,vspl=psps%vlspl)
542 : ! PAW only: we sometimes have to compute 1st-order compensation density
543 : ! and eventually add it to density from 1st-order WFs
544 : ! ----------------------------------------------------------------------
545 246 : if (psps%usepaw==1) then
546 :
547 : !Force the computation of nhatfr
548 738 : do iatom=1,dtset%natom
549 492 : pawfgrtab(iatom)%nhatfr_allocated = 0
550 382848 : pawfgrtab(iatom)%nhatfr = zero
551 : end do
552 :
553 : ! This portion of code works only when npert_phon<=1
554 246 : if (i1pert<=natom.and.usexcnhat==0) then
555 : call pawnhatfr(0,i1dir,i1pert,1,dtset%natom,nspden,psps%ntypat,&
556 : & pawang,pawfgrtab(i1pert),pawrhoij(i1pert),pawtab,rprimd,&
557 120 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
558 : end if
559 246 : if (i2pert<=natom) then
560 : call pawnhatfr(0,i2dir,i2pert,1,dtset%natom,nspden,psps%ntypat,&
561 : & pawang,pawfgrtab(i2pert),pawrhoij(i2pert),pawtab,rprimd,&
562 80 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
563 : end if
564 246 : if (i3pert<=natom.and.usexcnhat==0) then
565 : call pawnhatfr(0,i3dir,i3pert,1,dtset%natom,nspden,psps%ntypat,&
566 : & pawang,pawfgrtab(i3pert),pawrhoij(i3pert),pawtab,rprimd,&
567 0 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
568 : end if
569 :
570 246 : if (usexcnhat==0) then
571 :
572 : call pawmknhat(dummy_real,cplex,0,i1dir,i1pert,0,gprimd,natom,dtset%natom,&
573 : & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i1pert,&
574 : & pawrhoij1_i1pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,&
575 246 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
576 246 : if (flag1==0) then
577 262509 : rho1r1(:,:) = rho1r1(:,:) - nhat1_i1pert(:,:)
578 : flag1 = 1
579 : end if
580 :
581 : call pawmknhat(dummy_real,cplex,0,i3dir,i3pert,0,gprimd,natom,dtset%natom,&
582 : & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i3pert,&
583 : & pawrhoij1_i3pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,&
584 246 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
585 246 : if (flag3==0) then
586 682523 : rho3r1(:,:) = rho3r1(:,:) - nhat1_i3pert(:,:)
587 : flag3 = 1
588 : end if
589 :
590 : end if
591 :
592 : call pawmknhat(dummy_real,cplex,0,i2dir,i2pert,0,gprimd,natom,dtset%natom,&
593 : & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1_i2pert,&
594 : & pawrhoij1_i2pert,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,&
595 246 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
596 :
597 : end if
598 :
599 : else
600 :
601 : ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
602 : ! ------------------------------------------------------------------------------
603 45 : if(psps%n1xccc/=0)then
604 : call dfpt_mkcore(cplex,i2dir,i2pert,dtset%natom,psps%ntypat,n1,psps%n1xccc,&
605 45 : & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d2,xred)
606 : end if ! psps%n1xccc/=0
607 :
608 : call dfpt_vlocal(atindx,cplex,gmet,gsqcut,dtset%icutcoul,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,dtset%natom,&
609 : & nattyp,nfftf,ngfftf,nkpt,psps%ntypat,n1,n2,n3,ph1df,psps%qgrid_vl,&
610 45 : & dtset%qptn,dtset%rcut,rprimd,ucvol,dtset%vcutgeo,psps%vlspl,vpsp1,xred)
611 :
612 : end if ! usepaw
613 :
614 291 : option=1;optene=0
615 : call dfpt_rhotov(cplex,dummy_real,dummy_real,dummy_real,dummy_real,dummy_real,dummy_real,&
616 : & gsqcut,dtset%icutcoul,i2dir,i2pert,dtset%ixc,kxc,dtset%magpen,dtset%mpatpol,dtset%mpdir,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,&
617 : & nhat1_i2pert,nhat1gr,nhat1grdim,nkxc,nspden,dtset%ntypat,n3xccc,non_magnetic_xc,optene,option,&
618 : & dtset%qptn,dtset%ratsm,dtset%ratsph,rhog,rho2g1,rhor,rho2r1,rprimd,dtset%typat,ucvol,psps%usepaw,&
619 : & usexcnhat,dtset%vcutgeo,vhartr1_i2pert,&
620 291 : & vpsp1,vresid_dum,dummy_real,vtrial1_i2pert,vxc,vxc1_i2pert,xccc3d2,dtset%ixcrot,xred,dtset%qgbt,dtset%use_gbt)
621 :
622 291 : if (psps%usepaw==1.and.usexcnhat==0) then
623 2782587 : rho2r1(:,:) = rho2r1(:,:) - nhat1_i2pert(:,:)
624 : end if
625 :
626 291 : if (psps%usepaw==1)then
627 246 : call paw_an_reset_flags(paw_an1_i2pert) ! Force the recomputation of on-site potentials
628 246 : call paw_ij_reset_flags(paw_ij1_i2pert,all=.true.) ! Force the recomputation of Dij
629 246 : optfr=0
630 : call pawdijfr(gprimd,i2dir,i2pert,natom,natom,nfftf,ngfftf,nspden,nsppol,&
631 : & psps%ntypat,optfr,paw_ij1_i2pert,pawang,pawfgrtab,pawrad,pawtab,cplex,qphon,&
632 : & rprimd,ucvol,vpsp1,vtrial,vxc,xred,&
633 246 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
634 :
635 : ! Computation of "on-site" first-order potentials, first-order densities
636 : option=1
637 : call pawdenpot(dummy_real,el_temp,gprimd,i2pert,dtset%ixc,natom,dtset%natom,&
638 : & nspden,psps%ntypat,dtset%nucdipmom,0,option,paw_an1_i2pert,paw_an0,&
639 : & paw_energies_dum,paw_ij1_i2pert,pawang,dtset%pawprtvol,pawrad,&
640 : & pawrhoij1_i2pert,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,&
641 : & dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,xred,ucvol,psps%znuclpsp,dtset%spinaxis, &
642 246 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
643 : ! First-order Dij computation
644 : ! call timab(561,1,tsec)
645 246 : if (has_dijfr>0) then
646 : !vpsp1 contribution to Dij already stored in frozen part of Dij
647 984 : ABI_MALLOC(vtrial1_tmp,(cplex*nfftf,nspden))
648 5564928 : vtrial1_tmp=vtrial1_i2pert
649 723 : do ii=1,min(nspden,2)
650 2782587 : vtrial1_tmp(:,ii)=vtrial1_tmp(:,ii)-vpsp1(:)
651 : end do
652 : else
653 0 : vtrial1_tmp => vtrial1_i2pert
654 : end if
655 : call pawdij(cplex,dtset%enunit,gprimd,i2pert,natom,dtset%natom,&
656 : & nfftf,nfftotf,dtset%nspden,psps%ntypat,paw_an1_i2pert,paw_ij1_i2pert,pawang,&
657 : & pawfgrtab,dtset%pawprtvol,pawrad,pawrhoij1_i2pert,dtset%pawspnorb,pawtab,&
658 : & dtset%pawxcdev,qphon,dtset%spnorbscl,ucvol,dtset%cellcharge(1),&
659 : & vtrial1_tmp,vxc1_i2pert,xred,dtset%znucl,&
660 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,&
661 246 : & spinaxis=dtset%spinaxis)
662 246 : if (has_dijfr>0) then
663 246 : ABI_FREE(vtrial1_tmp)
664 : end if
665 : call symdij(gprimd,indsy1,i2pert,natom,dtset%natom,nsym1,psps%ntypat,0,&
666 : & paw_ij1_i2pert,pawang1,dtset%pawprtvol,pawtab,rprimd,symaf1,symrc1, &
667 : & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,&
668 246 : & qphon=qphon)
669 : ! call timab(561,2,tsec)
670 :
671 : end if ! end usepaw section
672 :
673 291 : nwffile = 1
674 291 : file_index(1) = i2dir + 3*(i2pert-1)
675 291 : fnamewff(1) = dtfil%fnamewff1
676 :
677 291 : if (i2pert==natom+2) then
678 :
679 199 : nwffile = 3
680 199 : file_index(2) = i2dir+natom*3
681 199 : fnamewff(2) = dtfil%fnamewffddk
682 : ! As npert_phon<=1 and i2pert==natom+2, i1pert or i3pert is necessarly equal to natom+2
683 199 : if (i3pert==natom+2) then
684 199 : second_idir = i3dir
685 0 : else if (i1pert==natom+2) then
686 0 : second_idir = i1dir
687 : else
688 0 : ABI_BUG(" i1pert or i3pert is supposed to be equal to natom+2, which is not the case here.")
689 : end if
690 199 : call rf2_getidir(i2dir,second_idir,idir_dkde)
691 199 : file_index(3) = idir_dkde+9+(dtset%natom+6)*3
692 199 : fnamewff(3) = dtfil%fnamewffdkde
693 :
694 199 : if (npert_phon==1.and.psps%usepaw==1.and.second_idir/=i2dir) then
695 80 : nwffile = 5
696 80 : file_index(4) = second_idir+natom*3
697 80 : fnamewff(4) = dtfil%fnamewffddk
698 80 : call rf2_getidir(second_idir,i2dir,idir_dkde) ! i2dir and second_idir are reversed
699 80 : file_index(5) = idir_dkde+9+(dtset%natom+6)*3
700 80 : fnamewff(5) = dtfil%fnamewffdkde
701 : end if
702 :
703 : end if
704 :
705 1140 : do ii=1,nwffile
706 849 : call appdig(file_index(ii),fnamewff(ii),fiwfddk)
707 : ! Checking the existence of data file
708 849 : if (.not. file_exists(fiwfddk)) then
709 : ! Trick needed to run Abinit test suite in netcdf mode.
710 849 : if (file_exists(nctk_ncify(fiwfddk))) then
711 849 : write(message,"(3a)")"- File: ",trim(fiwfddk),&
712 1698 : " does not exist but found netcdf file with similar name."
713 849 : call wrtout(std_out,message,'COLL')
714 1698 : fiwfddk = nctk_ncify(fiwfddk)
715 : end if
716 849 : if (.not. file_exists(fiwfddk)) then
717 0 : ABI_ERROR('Missing file: '//TRIM(fiwfddk))
718 : end if
719 : end if
720 849 : write(message,'(2a)')'-dfptnl_loop : read the wavefunctions from file: ',trim(fiwfddk)
721 849 : call wrtout(std_out,message,'COLL')
722 849 : call wrtout(ab_out,message,'COLL')
723 : ! Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50)
724 1140 : call ddk_f(ii)%open_read(fiwfddk,1,dtset%iomode,dtfil%unddk+(ii-1),mpi_enreg%comm_cell)
725 : end do
726 :
727 : ! Perform DFPT part of the 3dte calculation
728 291 : call timab(513,1,tsec)
729 : ! NOTE : eigen2 equals zero here
730 :
731 : call dfptnl_pert(atindx,cg,cg1,cg2,cg3,cplex,dtfil,dtset,d3etot,eigen0,gs_hamkq,k3xc,indsy1,i1dir,&
732 : & i2dir,i3dir,i1pert,i2pert,i3pert,kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,&
733 : & mpsang,mpw,natom,nattyp,nfftf,nfftotf,ngfftf,nkpt,nk3xc,nspden,nspinor,nsppol,nsym1,npwarr,occ,&
734 : & pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawtab,pawrhoij,pawrhoij1_i1pert,pawrhoij1_i2pert,pawrhoij1_i3pert,&
735 : & paw_an0,paw_an1_i2pert,paw_ij1_i2pert,ph1d,psps,rho1r1,rho2r1,rho3r1,&
736 : & rprimd,symaf1,symrc1,ucvol,vtrial,vhartr1_i2pert,vtrial1_i2pert,vxc1_i2pert,&
737 : & ddk_f,xccc3d1,xccc3d2,xccc3d3,xred,&
738 291 : & d3etot_1,d3etot_2,d3etot_3,d3etot_4,d3etot_5,d3etot_6,d3etot_7,d3etot_8,d3etot_9)
739 291 : call timab(513,2,tsec)
740 :
741 :
742 : ! Eventually close the dot file
743 1140 : do ii=1,nwffile
744 1140 : call ddk_f(ii)%close()
745 : end do
746 :
747 : ! if (psps%usepaw==1) then
748 : ! do ii=1,natom
749 : ! pawfgrtab(ii)%nhatfr = zero
750 : ! end do
751 : ! end if
752 :
753 : end if ! rfpert
754 : end do ! i2dir
755 : end do ! i2pert
756 :
757 : end if ! rfpert
758 : end do ! i3dir
759 : end do ! i3pert
760 :
761 : end if ! rfpert
762 : end do ! i1dir
763 : end do ! i1pert
764 :
765 : !More memory cleaning
766 7 : call gs_hamkq%free()
767 :
768 7 : ABI_FREE(cg1)
769 7 : ABI_FREE(cg2)
770 7 : ABI_FREE(cg3)
771 7 : ABI_FREE(eigen1)
772 7 : ABI_FREE(eigen2)
773 7 : ABI_FREE(eigen3)
774 7 : ABI_FREE(rho1r1)
775 7 : ABI_FREE(rho2r1)
776 7 : ABI_FREE(rho2g1)
777 7 : ABI_FREE(rho3r1)
778 7 : ABI_FREE(nhat1gr)
779 7 : ABI_FREE(vresid_dum)
780 7 : ABI_FREE(vtrial1_i2pert)
781 7 : ABI_FREE(vxc1_i2pert)
782 7 : ABI_FREE(vhartr1_i2pert)
783 7 : ABI_FREE(vpsp1)
784 7 : ABI_FREE(xccc3d1)
785 7 : ABI_FREE(xccc3d2)
786 7 : ABI_FREE(xccc3d3)
787 7 : if (psps%usepaw==1) then
788 4 : call pawrhoij_free(pawrhoij1_i1pert)
789 4 : call pawrhoij_free(pawrhoij1_i2pert)
790 4 : call pawrhoij_free(pawrhoij1_i3pert)
791 4 : ABI_FREE(nhat1_i1pert)
792 4 : ABI_FREE(nhat1_i2pert)
793 4 : ABI_FREE(nhat1_i3pert)
794 4 : call paw_an_free(paw_an1_i2pert)
795 4 : call paw_ij_free(paw_ij1_i2pert)
796 12 : ABI_FREE(paw_an1_i2pert)
797 12 : ABI_FREE(paw_ij1_i2pert)
798 : end if
799 15 : ABI_FREE(pawrhoij1_i1pert)
800 15 : ABI_FREE(pawrhoij1_i2pert)
801 15 : ABI_FREE(pawrhoij1_i3pert)
802 :
803 7 : call timab(503,2,tsec)
804 :
805 : DBG_EXIT("COLL")
806 :
807 21 : end subroutine dfptnl_loop
808 : !!***
809 :
810 : end module m_dfptnl_loop
811 : !!***
|
