Branch data 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 : : !!***
|
