Line data Source code
1 : !!****m* ABINIT/m_dfpt_loopert
2 : !! NAME
3 : !! m_dfpt_loopert
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1999-2026 ABINIT group (XG, DRH, MB, XW, MT, SPr, MJV)
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 : ! nvtx related macro definition
22 : #include "nvtx_macros.h"
23 :
24 : module m_dfpt_loopert
25 :
26 : use defs_basis
27 : use m_dtset
28 : use m_dtfil
29 : use defs_wvltypes
30 : use m_efmas_defs
31 : use m_abicore
32 : use m_xmpi
33 : use m_errors
34 : use m_wfk
35 : use m_wffile
36 : use m_io_redirect
37 : use m_paral_pert
38 : use m_nctk
39 : use m_ddb
40 : use m_wfd
41 : use m_ddb_hdr
42 : use netcdf
43 : use m_hdr
44 : use m_ebands
45 :
46 : use defs_datatypes, only : pseudopotential_type
47 : use defs_abitypes, only : MPI_type
48 : use m_occ, only : getnel
49 : use m_io_tools, only : file_exists
50 : use m_time, only : timab
51 : use m_fstrings, only : strcat, sjoin, ftoa
52 : use m_geometry, only : mkrdim, metric, littlegroup_pert
53 : use m_exit, only : exit_check, disable_timelimit
54 : use m_atomdata, only : atom_gauss
55 : use m_eig2d, only : eigr2d_init,eigr2d_t, eigr2d_ncwrite,eigr2d_free, &
56 : gkk_t, gkk_init, gkk_ncwrite,gkk_free, outbsd, eig2stern
57 : use m_crystal, only : crystal_t
58 : use m_efmas, only : efmas_main, efmas_analysis, print_efmas
59 : use m_fft, only : fourdp
60 : use m_fftcore, only : fftcore_set_mixprec
61 : use m_kg, only : getcut, getmpw, kpgio, getph
62 : use m_iowf, only : outwf, outresid
63 : use m_ioarr, only : read_rhor
64 : use m_orbmag, only : orbmag, orbmag_ncpp
65 : use m_pawang, only : pawang_type, pawang_init, pawang_free
66 : use m_pawrad, only : pawrad_type
67 : use m_pawtab, only : pawtab_type
68 : use m_paw_an, only : paw_an_type
69 : use m_paw_ij, only : paw_ij_type
70 : use m_pawfgrtab, only : pawfgrtab_type
71 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_bcast, pawrhoij_copy, &
72 : pawrhoij_nullify, pawrhoij_redistribute, pawrhoij_inquire_dim
73 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, pawcprj_getdim , pawcprj_output
74 : use m_pawfgr, only : pawfgr_type
75 : use m_paw_sphharm,only : setsym_ylm
76 : use m_rf2, only : rf2_getidirs
77 : use m_iogkk, only : outgkk
78 : use m_spacepar, only : rotate_rho, setsym
79 : use m_initylmg, only : initylmg
80 : use m_dfpt_scfcv, only : dfpt_scfcv
81 : use m_dfpt_mkrho, only : dfpt_mkrho
82 : use m_mpinfo, only : initmpi_band, distrb2, proc_distrb_cycle
83 : use m_atm2fft, only : dfpt_atm2fft
84 : use m_berrytk, only : smatrix
85 : use m_common, only : prteigrs
86 : use m_fourier_interpol, only : transgrid
87 : use m_mkcore, only : dfpt_mkcore
88 : use m_mklocl, only : dfpt_vlocal, vlocalstr
89 : use m_cgprj, only : ctocprj
90 : use m_symkpt, only : symkpt
91 : use m_pstat, only : pstat_proc
92 :
93 : #if defined(HAVE_GPU_MARKERS)
94 : use m_nvtx_data
95 : #endif
96 :
97 : implicit none
98 :
99 : private
100 : !!***
101 :
102 : public :: dfpt_looppert
103 : public :: eigen_meandege
104 : !!***
105 :
106 : contains
107 : !!***
108 :
109 : !!****f* ABINIT/dfpt_looppert
110 : !! NAME
111 : !! dfpt_looppert
112 : !!
113 : !! FUNCTION
114 : !! Loop over perturbations
115 : !!
116 : !! INPUTS
117 : !! atindx(natom)=index table for atoms (see gstate.f)
118 : !! codvsn=code version
119 : !! cpus=cpu time limit in seconds
120 : !! dim_eigbrd=1 if eigbrd is to be computed
121 : !! dim_eig2nkq=1 if eig2nkq is to be computed
122 : !! doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy
123 : !! dtfil <type(datafiles_type)>=variables related to files
124 : !! dtset <type(dataset_type)>=all input variables for this dataset
125 : !! dyew(2,3,natom,3,natom)=Ewald part of the dynamical matrix
126 : !! dyfrlo(3,3,natom)=frozen wavefunctions local part of the dynamical matrix
127 : !! dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wavefunctions non-local part of the dynamical matrix
128 : !! dyfrx1(2,3,natom,3,natom)=frozen wf nonlin. xc core corr.(2) part of the dynamical matrix
129 : !! dyfrx2(3,3,natom)=frozen wf nonlin. xc core corr.(2) part of the dynamical matrix
130 : !! dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix
131 : !! dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
132 : !! dyfr_nondiag=1 if dyfrnl is non diagonal with respect to atoms; 0 otherwise
133 : !! eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eigbrd)=broadening factors for the electronic eigenvalues
134 : !! eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq)=second derivatives of the electronic eigenvalues
135 : !! eltcore(6,6)=core contribution to the elastic tensor
136 : !! elteew(6+3*natom,6)=Ewald contribution to the elastic tensor
137 : !! eltfrhar(6,6)=Hartree contribution to the elastic tensor
138 : !! eltfrkin(6,6)=kinetic contribution to the elastic tensor
139 : !! eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor
140 : !! eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor
141 : !! eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor
142 : !! eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor
143 : !! fermie=fermi energy (Hartree)
144 : !! iexit=index of "exit" on first line of file (0 if not found)
145 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
146 : !! kxc(nfftf,nkxc)=exchange and correlation kernel (see rhotoxc.f)
147 : !! mkmem =Number of k points treated by this node (GS data)
148 : !! mkqmem=Number of k+q points treated by this node (GS data)
149 : !! mk1mem=Number of k points treated by this node (RF data)
150 : !! mpert=maximum number of ipert
151 : !! mpi_enreg=information about MPI parallelization
152 : !! my_natom=number of atoms treated by current processor
153 : !! nattyp(ntypat)= # atoms of each type.
154 : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
155 : !! nkpt=number of k points
156 : !! nkxc=second dimension of the kxc array
157 : !! nspden=number of spin-density components
158 : !! nsym=number of symmetry elements in space group
159 : !! occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
160 : !! paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS
161 : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
162 : !! pawang <type(pawang_type)>=paw angular mesh and related data
163 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
164 : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
165 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
166 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
167 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
168 : !! pertsy(3,mpert)=set of perturbations that form a basis for all other perturbations
169 : !! prtbbb=if 1, bbb decomposition, also dimension d2bbb
170 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
171 : !! rfpert(mpert)=array defining the type of perturbations that have to be computed
172 : !! 1 -> element has to be computed explicitly
173 : !! -1 -> use symmetry operations to obtain the corresponding element
174 : !! rhog(2,nfftf)=array for Fourier transform of GS electron density
175 : !! rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3.
176 : !! symq(4,2,nsym)=1 if symmetry preserves present qpoint. From littlegroup_q
177 : !! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
178 : !! timrev=1 if time-reversal preserves the q wavevector; 0 otherwise.
179 : !! usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory
180 : !! usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
181 : !! usevxctau=1 if if XC functional depends on kinetic energy density
182 : !! vtrial(nfftf,nspden)=GS potential (Hartree)
183 : !! vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree)
184 : !! vxcavg=average of vxc potential
185 : !! vxctau(nfftf,nspden,4*usevxctau)=derivative of e_xc with respect to kinetic energy density, for mGGA
186 : !! xred(3,natom)=reduced dimensionless atomic coordinates
187 : !!
188 : !! OUTPUT
189 : !! blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
190 : !! ddkfil(3)=unit numbers for the three possible ddk files
191 : !! d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some second order derivatives
192 : !! d2lo(2,mpert,3,mpert)=local contributions to the 2DTEs
193 : !! d2nl(2,mpert,3,mpert)=non-local contributions to the 2DTEs
194 : !! d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs
195 : !! etotal=total energy (sum of 8 contributions) (hartree)
196 : !!
197 : !! SOURCE
198 :
199 716 : subroutine dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,&
200 716 : & ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,&
201 1432 : & dyfr_cplex,dyfr_nondiag,d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasval,eigbrd,eig2nkq,&
202 716 : & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
203 1432 : & etotal,fermie,iexit,indsym,kxc,&
204 716 : & mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,&
205 1432 : & nfftf,nhat,nkpt,nkxc,nspden,nsym,occ,&
206 1432 : & paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
207 1432 : & pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,&
208 1331 : & usecprj,usevdw,usevxctau,vtrial,vxc,vxcavg,vxctau,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,&
209 : & eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0)
210 :
211 : !Arguments ------------------------------------
212 : integer, intent(in) :: dim_eigbrd,dim_eig2nkq,dyfr_cplex,dyfr_nondiag,mk1mem,mkmem,mkqmem,mpert
213 : integer, intent(in) :: nfftf,nkpt,nkxc,nspden,nsym,prtbbb,timrev,usecprj,usevdw,usevxctau
214 : integer, intent(out) :: iexit
215 : integer, intent(inout) :: my_natom
216 : real(dp), intent(in) :: cpus,vxcavg
217 : real(dp), intent(inout) :: fermie
218 : real(dp), intent(inout) :: etotal
219 : character(len=8), intent(in) :: codvsn
220 : type(MPI_type), intent(inout) :: mpi_enreg
221 : type(datafiles_type), intent(in) :: dtfil
222 : type(dataset_type), intent(in), target :: dtset
223 : type(pawang_type),intent(in) :: pawang
224 : type(pawfgr_type),intent(in) :: pawfgr
225 : type(pseudopotential_type), intent(inout) :: psps
226 : integer, intent(in) :: atindx(dtset%natom),indsym(4,nsym,dtset%natom)
227 : integer, intent(in) :: nattyp(dtset%ntypat),pertsy(3,mpert)
228 : integer, intent(in) :: rfpert(mpert),rf2_dirs_from_rfpert_nl(3,3),symq(4,2,nsym),symrec(3,3,nsym)
229 : integer, intent(out) :: ddkfil(3)
230 : integer, intent(inout) :: blkflg(3,mpert,3,mpert)
231 : integer, intent(out) :: clflg(3,mpert)
232 : real(dp), intent(in) :: doccde(dtset%mband*nkpt*dtset%nsppol)
233 : real(dp), intent(in) :: dyew(2,3,dtset%natom,3,dtset%natom)
234 : real(dp), intent(in) :: dyfrlo(3,3,dtset%natom)
235 : real(dp), intent(in) :: dyfrnl(dyfr_cplex,3,3,dtset%natom,1+(dtset%natom-1)*dyfr_nondiag)
236 : real(dp), intent(in) :: dyfrx1(2,3,dtset%natom,3,dtset%natom)
237 : real(dp), intent(in) :: dyfrx2(3,3,dtset%natom),dyvdw(2,3,dtset%natom,3,dtset%natom*usevdw)
238 : real(dp), intent(in) :: eltcore(6,6),elteew(6+3*dtset%natom,6),eltfrhar(6,6)
239 : real(dp), intent(in) :: eltfrkin(6,6),eltfrloc(6+3*dtset%natom,6)
240 : real(dp), intent(in) :: eltfrnl(6+3*dtset%natom,6)
241 : real(dp), intent(in) :: eltfrxc(6+3*dtset%natom,6),eltvdw(6+3*dtset%natom,6*usevdw)
242 : real(dp), intent(in) :: kxc(nfftf,nkxc),nhat(nfftf,nspden)
243 : real(dp), intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol)
244 : real(dp), intent(in) :: rhog(2,nfftf),rhor(nfftf,nspden),vxc(nfftf,nspden)
245 : real(dp), intent(in) :: vtrial(nfftf,nspden)
246 : real(dp), intent(inout) :: vxctau(nfftf,dtset%nspden,4*usevxctau)
247 : real(dp), intent(inout) :: xred(3,dtset%natom)
248 : real(dp), intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)!vz_i
249 : real(dp), intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i
250 : real(dp), intent(inout) :: d2ovl(2,3,mpert,3,mpert*psps%usepaw) !vz_i
251 : real(dp), intent(out) :: eigbrd(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eigbrd)
252 : real(dp), intent(out) :: eig2nkq(2,dtset%mband*dtset%nsppol,nkpt,3,dtset%natom,3,dtset%natom*dim_eig2nkq)
253 : type(efmasdeg_type),allocatable,intent(inout) :: efmasdeg(:)
254 : type(efmasval_type),allocatable,intent(inout) :: efmasval(:,:)
255 : type(paw_an_type),allocatable,target,intent(inout) :: paw_an(:)
256 : type(paw_ij_type),allocatable,target,intent(inout) :: paw_ij(:)
257 : type(pawfgrtab_type),allocatable,target,intent(inout) :: pawfgrtab(:)
258 : type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
259 : type(pawrhoij_type),allocatable,target,intent(inout) :: pawrhoij(:)
260 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
261 : real(dp),pointer :: eigen1_pert(:,:,:)
262 : real(dp),intent(out) :: occ_rbz_pert(:),eigen0_pert(:),eigenq_pert(:)
263 : real(dp),pointer :: eigenq_fine(:,:,:)
264 : integer, intent(out) :: nkpt_rbz
265 : type(hdr_type),intent(out) :: hdr0,hdr_fine
266 :
267 : !Local variables-------------------------------
268 : !scalars
269 : integer,parameter :: level=11,response=1,formeig1=1,master=0,fake_unit=-666
270 : integer :: ask_accurate,band_index,bantot,bantot_rbz,bdeigrf,bdtot1_index,nsppol,nspinor,band2tot_index
271 : integer :: bdtot_index,choice,cplex,cplex_rhoij,dim_eig2rf,formeig
272 : integer :: gscase,g0term,iband,iblok,icase,icase_eq,idir,idir0,idir1,idir2,idir_eq,idir_dkdk,ierr
273 : integer :: ifft,ii,ikpt,ikpt1,jband,initialized,iorder_cprj,ipert,ipert_cnt,ipert_eq,ipert_me,ireadwf0
274 : integer :: iscf_mod,iscf_mod_save,isppol,istr,isym,mcg,mcgq,mcg1,mcprj,mcprjq,mband
275 : integer :: mband_mem_rbz
276 : integer :: mcgmq,mcg1mq,mpw1_mq !+/-q duplicates
277 : integer :: maxidir,me,mgfftf,mkmem_rbz,mk1mem_rbz,mkqmem_rbz,mpw,mpw1,my_nkpt_rbz
278 : integer :: n3xccc,nband_k,ncpgr,ndir,nkpt_eff,nkpt_max,nline_save,nmatel,npert_io,npert_me,nspden_rhoij
279 : integer :: nstep_save,nsym1,ntypat,nwffile,nylmgr,nylmgr1,old_comm_atom,openexit,option,optorth,optthm,pertcase,pertcase_mq
280 : integer :: qphase_rhoij,rdwr,rdwrpaw,spaceComm,smdelta,timrev_pert,timrev_kpt,to_compute_this_pert
281 : integer :: useylmgr,useylmgr1,dfpt_scfcv_retcode,optn2
282 : integer :: ncerr,ncid
283 : real(dp) :: boxcut,dosdeltae,eberry,ecore,ecut_eff,ecutf,edocc,eei,eeig0,eew,efrhar,efrkin,efrloc
284 : real(dp) :: efrnl,efrx1,efrx2,ehart,ehart01,ehart1,eii,ek,ek0,ek1,ek2,elmag1,eloc0
285 : real(dp) :: elpsp1,emagpen1,enl,enl0,enl1,end0,end1,entropy,enxc,eovl1,epaw1,evxctau0,evxctau1,evdw,exc1
286 : real(dp) :: residm,tolwfr,tolwfr_save,toldfe_save,toldff_save,tolrff_save,tolvrs_save,toldmag_save
287 : real(dp) :: fsum,gsqcut,maxocc,nelectkq
288 : real(dp) :: ucvol, eig1_r, eig1_i
289 : real(dp) :: residm_mq !+/-q duplicates
290 : logical,parameter :: paral_pert_inplace=.true.,remove_inv=.false.
291 : logical :: first_entry,found_eq_gkk,has_nd,t_exist,paral_atom,write_1wfk,init_rhor1,with_vxctau
292 : logical :: kramers_deg
293 : character(len=fnlen) :: fiden1i,fiwf1i,fiwf1i_mq,fiwf1o,fiwf1o_mq,fiwfddk,fnamewff(4),gkkfilnam,fname,filnam, fnamewffmq_
294 : character(len=500) :: msg
295 74464 : type(crystal_t) :: crystal, ddb_crystal
296 : type(dataset_type), pointer :: dtset_tmp
297 716 : type(ebands_t) :: ebands_k,ebands_kq,gkk_ebands, ebands_kmq !+/-q duplicates
298 716 : type(gkk_t) :: gkk2d
299 716 : type(hdr_type) :: hdr,hdr_mq,hdr_den,hdr_tmp
300 37232 : type(ddb_hdr_type) :: ddb_hdr
301 716 : type(pawang_type) :: pawang1
302 92364 : type(wfk_t) :: ddk_f(4)
303 716 : type(wvl_data) :: wvl
304 : !arrays
305 : integer :: eq_symop(3,3),ngfftf(18),file_index(4),rfdir(9),rf2dir(9),rf2_dir1(3),rf2_dir2(3) !, units(2)
306 1432 : integer,allocatable :: blkflg_save(:,:,:,:),dimcprj_srt(:),dyn(:),indkpt1(:),indkpt1_tmp(:)
307 716 : integer,allocatable :: indsy1(:,:,:),irrzon1(:,:,:),istwfk_rbz(:),istwfk_pert(:,:,:)
308 716 : integer,allocatable :: kg(:,:),kg1(:,:),nband_rbz(:),npwar1(:),npwarr(:),npwtot(:)
309 716 : integer,allocatable :: kg1_mq(:,:),npwar1_mq(:),npwtot1_mq(:) !+q/-q duplicates
310 716 : integer,allocatable :: npwtot1(:),npwar1_pert(:,:),npwarr_pert(:,:),npwtot_pert(:,:)
311 1432 : integer,allocatable :: pert_calc(:,:),pert_tmp(:,:),bz2ibz_smap(:,:)
312 716 : integer,allocatable :: symaf1(:),symaf1_tmp(:),symrc1(:,:,:),symrl1(:,:,:),symrl1_tmp(:,:,:)
313 716 : integer, pointer :: old_atmtab(:)
314 716 : logical, allocatable :: distrb_flags(:,:,:)
315 : real(dp) :: dielt(3,3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),tsec(2)
316 716 : real(dp),allocatable :: buffer1(:,:,:,:,:),cg(:,:),cg1(:,:),cg1_active(:,:),cg1_3(:,:,:),cg0_pert(:,:)
317 1432 : real(dp),allocatable :: cg1_pert(:,:,:,:),cgq(:,:),gh0c1_pert(:,:,:,:)
318 716 : real(dp),allocatable :: doccde_rbz(:),docckqde(:)
319 716 : real(dp),allocatable :: gh1c_pert(:,:,:,:),eigen0(:),eigen0_copy(:),eigen1(:),eigen1_mean(:)
320 716 : real(dp),allocatable :: eigenq(:),gh1c_set(:,:),gh0c1_set(:,:),kpq(:,:)
321 716 : real(dp),allocatable :: kpq_rbz(:,:),kpt_rbz(:,:),occ_pert(:),occ_rbz(:),occkq(:),kpt_rbz_pert(:,:)
322 716 : real(dp),allocatable :: occ_disk(:)
323 716 : real(dp),allocatable :: vtrial_local(:,:)
324 716 : real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons1(:,:,:),resid(:),rhog1(:,:)
325 716 : real(dp),allocatable :: rhor1_save(:,:,:)
326 716 : real(dp),allocatable :: rhor1(:,:),rho1wfg(:,:),rho1wfr(:,:),tnons1(:,:),tnons1_tmp(:,:)
327 716 : real(dp),allocatable :: rhor1_pq(:,:),rhor1_mq(:,:),rhog1_pq(:,:),rhog1_mq(:,:) !+q/-q duplicates
328 1432 : real(dp),allocatable :: cg_mq(:,:),cg1_mq(:,:),resid_mq(:) !
329 1432 : real(dp),allocatable :: cg1_active_mq(:,:),occk_mq(:) !
330 716 : real(dp),allocatable :: kmq(:,:),kmq_rbz(:,:),gh0c1_set_mq(:,:) !
331 716 : real(dp),allocatable :: eigen_mq(:),gh1c_set_mq(:,:),docckde_mq(:),eigen1_mq(:) !
332 716 : real(dp),allocatable :: vpsp1(:),work(:),wtk_folded(:),wtk_rbz(:),xccc3d1(:)
333 716 : real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:),zeff(:,:,:)
334 716 : real(dp),allocatable :: ylm1_mq(:,:),ylmgr1_mq(:,:,:)
335 716 : real(dp),allocatable :: phasecg(:,:),gauss(:,:)
336 716 : real(dp),allocatable :: gkk(:,:,:,:,:)
337 : logical :: has_cg1_3(3)
338 716 : type(pawcprj_type),allocatable :: cprj(:,:),cprjq(:,:)
339 716 : type(paw_ij_type),pointer :: paw_ij_pert(:)
340 716 : type(paw_an_type),pointer :: paw_an_pert(:)
341 716 : type(pawfgrtab_type),pointer :: pawfgrtab_pert(:)
342 716 : type(pawrhoij_type),allocatable :: pawrhoij1(:)
343 716 : type(pawrhoij_type),pointer :: pawrhoij_pert(:)
344 716 : type(ddb_type) :: ddb
345 716 : real(dp),allocatable :: doccde_tmp(:)
346 : ! ***********************************************************************
347 :
348 : DBG_ENTER("COLL")
349 :
350 716 : call timab(141,1,tsec)
351 :
352 : !Structured debugging if prtvol==-level
353 716 : if(dtset%prtvol==-level)then
354 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' dfpt_looppert : enter , debug mode '
355 0 : call wrtout(std_out,msg)
356 : end if
357 :
358 716 : dfpt_scfcv_retcode = -1
359 716 : nsppol = dtset%nsppol; nspinor = dtset%nspinor
360 :
361 : !Deactivate kramer's degeneracy for calculations that break TRS
362 716 : kramers_deg=.true.
363 716 : if (dtset%tim1rev==0) kramers_deg=.false.
364 :
365 716 : if (.not.kramers_deg.and.dtset%npfft/=1) then
366 : write(msg,'(a,a)')' RF routines are not adapted to fft parallelization'//char(10)//&
367 0 : & ' when tim1rev=0. However, npfft is different from 1 in the present run.'
368 0 : ABI_BUG(msg)
369 0 : call wrtout(ab_out,msg,'COLL')
370 : end if
371 :
372 716 : with_vxctau = (usevxctau > 0)
373 716 : evxctau0=zero; evxctau1=0
374 6792 : has_nd = ANY(ABS(dtset%nucdipmom(:,:))>tol8)
375 716 : end0=zero; end1=zero
376 :
377 : !Obtain dimensional translations in reciprocal space gprimd,
378 : !metrics and unit cell volume, from rprimd. Also output rprimd, gprimd and ucvol
379 716 : call mkrdim(dtset%acell_orig(1:3,1),dtset%rprim_orig(1:3,1:3,1),rprimd)
380 716 : call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol)
381 :
382 : call crystal%init(dtset%amu_orig(:,1),dtset%spgroup,dtset%natom,dtset%npsp,&
383 : psps%ntypat,dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
384 : dtset%nspden==2.and.dtset%nsppol==1,remove_inv,psps%title,&
385 1432 : symrel=dtset%symrel,tnons=dtset%tnons,symafm=dtset%symafm)
386 :
387 : !Get FFT grid(s) sizes (be careful !) See NOTES in the comments at the beginning of respfn.F90
388 716 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
389 1786 : mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:)
390 94 : ecutf=dtset%pawecutdg
391 : else
392 11818 : mgfftf=dtset%mgfft;ngfftf(:)=dtset%ngfft(:)
393 622 : ecutf=dtset%ecut
394 : end if
395 716 : ecut_eff=dtset%ecut*(dtset%dilatmx)**2
396 :
397 : !Compute large sphere cut-off gsqcut
398 716 : if (psps%usepaw==1) then
399 112 : call wrtout(std_out,ch10//' FFT (fine) grid used for densities/potentials:')
400 : end if
401 716 : call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,ngfftf)
402 :
403 : !Various initializations/allocations
404 716 : iscf_mod=dtset%iscf
405 716 : ntypat=psps%ntypat
406 716 : nkpt_max=50;if (xmpi_paral==1) nkpt_max=-1
407 : !TODO: this flag for paral_atom is ignored below
408 716 : paral_atom=(dtset%natom/=my_natom)
409 716 : cplex=2-timrev !cplex=2 ! DEBUG: impose cplex=2
410 716 : first_entry=.true.
411 716 : initialized=0
412 716 : ecore=zero ; ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero ; eii=zero
413 72428 : d2bbb = zero
414 2000742 : d2nl = zero
415 2000742 : d2lo = zero
416 335358 : d2ovl = zero
417 27552 : clflg(:,:)=0 ! Array on calculated perturbations for eig2rf
418 716 : if (psps%usepaw==1) then
419 336 : ABI_MALLOC(dimcprj_srt,(dtset%natom))
420 112 : call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
421 : end if
422 716 : mband_mem_rbz = dtset%mband_mem
423 :
424 : !Save values of SCF cycle parameters
425 716 : iscf_mod_save = iscf_mod
426 716 : nstep_save = dtset%nstep
427 716 : nline_save = dtset%nline
428 716 : tolwfr_save = dtset%tolwfr
429 716 : toldfe_save = dtset%toldfe
430 716 : toldff_save = dtset%toldff
431 716 : toldmag_save = dtset%toldmag
432 716 : tolrff_save = dtset%tolrff
433 716 : tolvrs_save = dtset%tolvrs
434 :
435 : !This dtset will be used in dfpt_scfcv to force non scf calculations for equivalent perturbations
436 716 : nullify(dtset_tmp)
437 716 : if (dtset%prepgkk/=0) then ! .and. dtset%use_nonscf_gkk==1) then !Later uncomment this - in scf case rhor1_save is used below only for testing
438 44 : ABI_MALLOC(dtset_tmp,)
439 44 : dtset_tmp = dtset%copy()
440 : else
441 : dtset_tmp => dtset
442 : end if
443 :
444 : !Generate the 1-dimensional phases
445 2148 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
446 2148 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
447 716 : call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred)
448 716 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
449 94 : call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
450 : else
451 493732 : ph1df(:,:)=ph1d(:,:)
452 : end if
453 :
454 : !!Determine existence of perturbations and of perturbation symmetries
455 : !!Create array with perturbations which have to be calculated
456 : ! ABI_MALLOC(pert_tmp,(3*mpert))
457 : ! ipert_cnt=0
458 : ! do ipert=1,mpert
459 : ! do idir=1,3
460 : ! if( rfpert(ipert)==1 .and. dtset%rfdir(idir) == 1 )then
461 : ! if ((pertsy(idir,ipert)==1).or.&
462 : !& ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.&
463 : !& ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then
464 : ! ipert_cnt = ipert_cnt+1;
465 : ! pert_tmp(ipert_cnt) = idir+(ipert-1)*3
466 : ! else
467 : ! write(msg, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,&
468 : !& ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,&
469 : !& ' symmetric of a previously calculated perturbation.',ch10,&
470 : !& ' So, its SCF calculation is not needed.',ch10
471 : ! call wrtout(std_out,msg,'COLL')
472 : ! call wrtout(ab_out,msg,'COLL')
473 : ! end if ! Test of existence of symmetry of perturbation
474 : ! end if ! Test of existence of perturbation
475 : ! end do
476 : ! end do
477 : ! ABI_MALLOC(pert_calc,(ipert_cnt))
478 : ! do icase=1,ipert_cnt
479 : ! pert_calc(icase)=pert_tmp(icase)
480 : ! end do
481 : ! ABI_FREE(pert_tmp)
482 :
483 : !Initialize rf2dir :
484 2864 : rf2_dir1(1:3)=dtset%rf2_pert1_dir(1:3)
485 2864 : rf2_dir2(1:3)=dtset%rf2_pert2_dir(1:3)
486 5012 : if (sum(rf2_dir1)==3.and.sum(rf2_dir2)==3.and.dtset%prepanl==1) then
487 : ! Diagonal terms :
488 38 : rf2dir(1) = rf2_dirs_from_rfpert_nl(1,1)
489 38 : rf2dir(2) = rf2_dirs_from_rfpert_nl(2,2)
490 38 : rf2dir(3) = rf2_dirs_from_rfpert_nl(3,3)
491 : ! Upper triangular terms :
492 38 : rf2dir(4) = rf2_dirs_from_rfpert_nl(2,3)
493 38 : rf2dir(5) = rf2_dirs_from_rfpert_nl(1,3)
494 38 : rf2dir(6) = rf2_dirs_from_rfpert_nl(1,2)
495 : ! Lower triangular terms :
496 38 : rf2dir(7) = rf2_dirs_from_rfpert_nl(3,2)
497 38 : rf2dir(8) = rf2_dirs_from_rfpert_nl(3,1)
498 38 : rf2dir(9) = rf2_dirs_from_rfpert_nl(2,1)
499 : else
500 : ! Diagonal terms :
501 678 : rf2dir(1) = rf2_dir1(1)*rf2_dir2(1)
502 678 : rf2dir(2) = rf2_dir1(2)*rf2_dir2(2)
503 678 : rf2dir(3) = rf2_dir1(3)*rf2_dir2(3)
504 : ! Upper triangular terms :
505 678 : rf2dir(4) = rf2_dir1(2)*rf2_dir2(3)
506 678 : rf2dir(5) = rf2_dir1(1)*rf2_dir2(3)
507 678 : rf2dir(6) = rf2_dir1(1)*rf2_dir2(2)
508 : ! Lower triangular terms :
509 678 : rf2dir(7) = rf2_dir1(3)*rf2_dir2(2)
510 678 : rf2dir(8) = rf2_dir1(3)*rf2_dir2(1)
511 678 : rf2dir(9) = rf2_dir1(2)*rf2_dir2(1)
512 : end if
513 :
514 : !Determine existence of pertubations and of pertubation symmetries
515 : !Create array with pertubations which have to be calculated
516 2148 : ABI_MALLOC(pert_tmp,(3,3*(dtset%natom+6)+18+3*dtset%natom))
517 716 : ipert_cnt=0
518 7425 : do ipert=1,mpert
519 6709 : if (ipert<dtset%natom+10.or.ipert>dtset%natom+11) then
520 26580 : maxidir = 3
521 26580 : rfdir(1:3) = dtset%rfdir(:)
522 46515 : rfdir(4:9) = 0
523 : else
524 64 : maxidir = 9
525 64 : rfdir(1:9) = rf2dir(:)
526 : end if
527 : !Scalar potential case
528 6709 : if (ipert==dtset%natom+6) then
529 716 : maxidir = 1
530 716 : rfdir(:) = 0
531 716 : rfdir(1) = 1
532 : end if
533 26504 : do idir=1,maxidir
534 19079 : to_compute_this_pert = 0
535 19079 : if(ipert<dtset%natom+10 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then
536 : if ((pertsy(idir,ipert)==1).or.&
537 3064 : & ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.&
538 : & ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then
539 : to_compute_this_pert = 1
540 : else
541 869 : write(msg, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,&
542 869 : ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,&
543 869 : ' symmetric of a previously calculated perturbation.',ch10,&
544 1738 : ' So, its SCF calculation is not needed.',ch10
545 2607 : call wrtout([std_out, ab_out], msg)
546 : end if ! Test of existence of symmetry of perturbation
547 16015 : else if (ipert==dtset%natom+11 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then
548 : to_compute_this_pert = 1
549 15972 : else if (ipert==dtset%natom+10 .and. rfpert(ipert)==1) then
550 162 : if (dtset%rf2_dkdk==2 .or. dtset%rf2_dkdk==3) then
551 99 : if (idir <= 3 .and. rfdir(idir) == 1) then
552 : to_compute_this_pert = 1
553 66 : else if (idir>=4.and.idir<=6) then
554 33 : if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1
555 33 : else if (idir>=7.and.idir<=9) then
556 33 : if (rfdir(idir) == 1 .or. rfdir(idir-3) == 1) to_compute_this_pert = 1
557 : end if
558 : else
559 63 : if (idir<=6) then
560 42 : if (idir<=3) then
561 21 : if (rfdir(idir) == 1) to_compute_this_pert = 1
562 : else
563 21 : if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1
564 : end if
565 : end if
566 : end if
567 15810 : else if (ipert==dtset%natom+11 .and. rfpert(ipert)==1) then
568 20 : if (idir <= 3 .and. rfdir(idir) == 1) then
569 : to_compute_this_pert = 1
570 20 : else if (idir>=4.and.idir<=6) then
571 10 : if (rfdir(idir) == 1 .or. rfdir(idir+3) == 1) to_compute_this_pert = 1
572 10 : else if (idir>=7.and.idir<=9) then
573 5 : if (rfdir(idir) == 1 .or. rfdir(idir-3) == 1) to_compute_this_pert = 1
574 : end if
575 : else if (ipert>dtset%natom+11 .and. ipert<=2*dtset%natom+11 .and. &
576 15790 : & rfpert(ipert)==1 .and. rfdir(idir) == 1) then
577 : to_compute_this_pert = 1
578 : end if
579 6709 : if (to_compute_this_pert /= 0) then
580 2377 : ipert_cnt = ipert_cnt+1;
581 2377 : pert_tmp(1,ipert_cnt) = ipert
582 2377 : pert_tmp(2,ipert_cnt) = idir
583 : ! Store "pertcase" in pert_tmp(3,ipert_cnt)
584 2377 : if (ipert<dtset%natom+10) then
585 2195 : pert_tmp(3,ipert_cnt) = idir + (ipert-1)*3
586 182 : else if (ipert==dtset%natom+10.or.ipert<=dtset%natom+11) then
587 179 : pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-10)*9 + (dtset%natom+6)*3
588 3 : else if (ipert>dtset%natom+11 .and. ipert<=2*dtset%natom+11) then
589 3 : pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-12)*3 + (dtset%natom+6)*3+18
590 : end if
591 : end if
592 : end do ! idir
593 : end do !ipert
594 : ! do ipert=1,mpert
595 : ! if (ipert<dtset%natom+10) then
596 : ! maxidir = 3
597 : ! rfdir(1:3) = dtset%rfdir(:)
598 : ! rfdir(4:9) = 0
599 : ! else
600 : ! maxidir = 9
601 : ! rfdir(1:9) = rf2dir(:)
602 : ! end if
603 : ! do idir=1,maxidir
604 : ! if( rfpert(ipert)==1 .and. rfdir(idir) == 1 )then
605 : ! to_compute_this_pert = 0
606 : ! if (ipert>=dtset%natom+10) then
607 : ! to_compute_this_pert = 1
608 : ! else if ((pertsy(idir,ipert)==1).or.&
609 : !& ((dtset%prepanl == 1).and.(ipert == dtset%natom+2)).or.&
610 : !& ((dtset%prepgkk == 1).and.(ipert <= dtset%natom)) ) then
611 : ! to_compute_this_pert = 1
612 : ! end if
613 : ! if (to_compute_this_pert /= 0) then
614 : ! ipert_cnt = ipert_cnt+1;
615 : ! pert_tmp(1,ipert_cnt) = ipert
616 : ! pert_tmp(2,ipert_cnt) = idir
617 : !! Store "pertcase" in pert_tmp(3,ipert_cnt)
618 : ! if (ipert<dtset%natom+10) then
619 : ! pert_tmp(3,ipert_cnt) = idir + (ipert-1)*3
620 : ! else
621 : ! pert_tmp(3,ipert_cnt) = idir + (ipert-dtset%natom-10)*9 + (dtset%natom+6)*3
622 : ! end if
623 : ! else
624 : ! write(msg, '(a,a,i4,a,i4,a,a,a,a,a,a)' )ch10,&
625 : !& ' The perturbation idir=',idir,' ipert=',ipert,' is',ch10,&
626 : !& ' symmetric of a previously calculated perturbation.',ch10,&
627 : !& ' So, its SCF calculation is not needed.',ch10
628 : ! call wrtout(std_out,msg,'COLL')
629 : ! call wrtout(ab_out,msg,'COLL')
630 : ! end if ! Test of existence of symmetry of perturbation
631 : ! end if ! Test of existence of perturbation
632 : ! end do
633 : ! end do
634 2148 : ABI_MALLOC(pert_calc,(3,ipert_cnt))
635 3093 : do icase=1,ipert_cnt
636 10224 : pert_calc(:,icase)=pert_tmp(:,icase)
637 : end do
638 716 : ABI_FREE(pert_tmp)
639 :
640 716 : if (dtset%prepgkk/=0) then ! .and. dtset%use_nonscf_gkk==1) then !Later uncomment this - in scf case rhor1_save is used below only for testing
641 220 : ABI_MALLOC(rhor1_save,(cplex*nfftf,nspden,ipert_cnt))
642 321524 : rhor1_save=zero
643 220 : ABI_MALLOC(blkflg_save,(3,mpert,3,mpert))
644 : end if
645 :
646 : ! Initialize quantities for netcdf print
647 2148 : ABI_MALLOC(eigen0_copy,(dtset%mband*nkpt*dtset%nsppol))
648 339337 : eigen0_copy(:)=zero
649 :
650 : ! SP : Retrieval of the DDB information and computing of effective charge and
651 : ! dielectric tensor
652 2148 : ABI_MALLOC(zeff,(3,3,dtset%natom))
653 716 : if (dtset%getddb .ne. 0 .or. dtset%irdddb .ne. 0 ) then
654 2 : filnam = dtfil%filddbsin
655 2 : call ddb%from_file(filnam, ddb_hdr, ddb_crystal, mpi_enreg%comm_world)
656 2 : call ddb_hdr%free()
657 : ! Get Dielectric Tensor and Effective Charges
658 : ! (initialized to one_3D and zero if the derivatives are not available in the DDB file)
659 2 : iblok = ddb%get_dielt_zeff(ddb_crystal,1,0,0,dielt,zeff)
660 2 : call ddb_crystal%free()
661 2 : call ddb%free()
662 : end if
663 :
664 : !%%%% Parallelization over perturbations %%%%%
665 : !*Define file output/log file names
666 716 : npert_io=ipert_cnt;if (dtset%nppert<=1) npert_io=0
667 716 : call localfilnam(mpi_enreg%comm_pert,mpi_enreg%comm_cell_pert,mpi_enreg%comm_world,dtfil%filnam_ds,'_PRT',npert_io)
668 : !Compute the number of perturbation done by the current cpu
669 716 : if(mpi_enreg%paral_pert==1) then
670 26 : npert_me = 0 ; ipert_me = 0
671 100 : do icase=1,ipert_cnt
672 100 : if (mpi_enreg%distrb_pert(icase)==mpi_enreg%me_pert) npert_me=npert_me +1
673 : end do
674 : end if
675 :
676 : !*Redefine communicators
677 716 : call set_pert_comm(mpi_enreg,dtset%nppert)
678 :
679 : !*Redistribute PAW on-site data
680 716 : nullify(old_atmtab,pawfgrtab_pert,pawrhoij_pert,paw_an_pert,paw_ij_pert)
681 : if (paral_pert_inplace) then
682 : call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
683 716 : & paw_an,paw_ij,pawfgrtab,pawrhoij)
684 716 : pawfgrtab_pert=>pawfgrtab ; pawrhoij_pert=>pawrhoij
685 716 : paw_an_pert =>paw_an ; paw_ij_pert =>paw_ij
686 :
687 : else
688 : call set_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
689 : & paw_an,paw_ij,pawfgrtab,pawrhoij,&
690 : & paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,&
691 : & pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert)
692 :
693 : end if
694 :
695 : ! We can handle the time limit in dfpt_scfcv in a robust manner only if we have one perturbation.
696 716 : if (ipert_cnt > 1 .or. mpi_enreg%paral_pert == 1) call disable_timelimit()
697 :
698 : !Loop on perturbations
699 : !==========================================================================
700 : ABI_NVTX_START_RANGE(NVTX_DFPT_LOOP)
701 3093 : do icase=1,ipert_cnt
702 :
703 : ! %%%% Parallelization over perturbations %%%%%
704 : ! Select the perturbations treated by current processor
705 2377 : if(mpi_enreg%paral_pert==1) then
706 74 : if (mpi_enreg%distrb_pert(icase)/=mpi_enreg%me_pert) cycle
707 : end if
708 :
709 : ! Redefine output/log files
710 2342 : call localwrfile(mpi_enreg%comm_cell,icase,npert_io,mpi_enreg%paral_pert,0)
711 :
712 : ! Set precision for FFT libs.
713 2342 : ii = fftcore_set_mixprec(dtset%mixprec)
714 :
715 : !! Retrieve type and direction of the perturbation
716 : ! if (pert_calc(icase) <= dtset%natom*3) then
717 : ! idir = mod(pert_calc(icase),3)
718 : ! if (idir==0) idir=3
719 : ! ipert=( (pert_calc(icase)-idir) / 3 + 1)
720 : ! else if (pert_calc(icase) <= dtset%natom*3+4) then
721 : ! ipert = dtset%natom + ((pert_calc(icase) - 3*dtset%natom - 1) / 3) + 1
722 : ! idir = mod(pert_calc(icase),3)
723 : ! if (idir==0) idir=3
724 : ! else
725 : ! ipert = dtset%natom + ((pert_calc(icase) - 3*(dtset%natom+4) - 1) / 9) + 1
726 : ! end if
727 : ! pertcase=idir+(ipert-1)*3
728 :
729 : ! Retrieve type and direction of the perturbation
730 2342 : ipert = pert_calc(1,icase)
731 2342 : idir = pert_calc(2,icase)
732 2342 : istr=idir
733 2342 : pertcase = pert_calc(3,icase)
734 :
735 : if(ipert <= dtset%natom) then
736 : ABI_NVTX_START_RANGE(NVTX_DFPT_LOOP_PHONON)
737 : end if
738 : if(ipert == dtset%natom+1) then
739 : ABI_NVTX_START_RANGE(NVTX_DFPT_LOOP_DDK)
740 : end if
741 : if(ipert == dtset%natom+2) then
742 : ABI_NVTX_START_RANGE(NVTX_DFPT_LOOP_EFELD)
743 : end if
744 : if(ipert == dtset%natom+3 .or. ipert == dtset%natom+4) then
745 : ABI_NVTX_START_RANGE(NVTX_DFPT_LOOP_STRAIN)
746 : end if
747 :
748 :
749 : ! Init MPI communicator
750 2342 : spaceComm=mpi_enreg%comm_cell
751 2342 : me=mpi_enreg%me_cell
752 :
753 : ! ===== Describe the perturbation in output/log file
754 :
755 192044 : write(msg, '(a,80a,a,a,3f10.6)' ) ch10,('-',ii=1,80),ch10,&
756 194386 : ' Perturbation wavevector (in red.coord.) ',dtset%qptn(:)
757 7026 : call wrtout([std_out, ab_out],msg)
758 2342 : if (abs(dtset%rfomega) > tol10) then
759 : write(msg, '(a,3f10.6)' ) &
760 4 : ' Perturbation frequency (in hartree a.u.) ',dtset%rfomega
761 12 : call wrtout([std_out, ab_out],msg)
762 : end if
763 2342 : if(ipert>=1 .and. ipert<=dtset%natom)then
764 1353 : write(msg, '(a,i4,a,i4)' )' Perturbation : displacement of atom',ipert,' along direction',idir
765 4059 : call wrtout([std_out, ab_out], msg)
766 1353 : if(iscf_mod == -3)then
767 2 : write(msg, '(a,a,a,a,a,a,a,a)' )ch10,&
768 2 : ' dfpt_looppert : COMMENT -',ch10,&
769 2 : ' The first-order density is imposed to be zero (iscf=-3).',ch10,&
770 2 : ' Although this is strange in the case of phonons,',ch10,&
771 4 : ' you are allowed to do so.'
772 6 : call wrtout([std_out, ab_out], msg)
773 : end if
774 989 : else if(ipert==dtset%natom+1)then
775 436 : write(msg,'(a,i4)')' Perturbation : derivative vs k along direction',idir
776 1308 : call wrtout([std_out, ab_out], msg)
777 436 : if( iscf_mod /= -3 )then
778 46 : write(msg, '(4a)' )ch10,&
779 46 : ' dfpt_looppert : COMMENT -',ch10,&
780 92 : ' In a d/dk calculation, iscf is set to -3 automatically.'
781 138 : call wrtout([std_out, ab_out], msg)
782 46 : iscf_mod=-3
783 : end if
784 436 : if( abs(dtset%dfpt_sciss) > 1.0d-8 )then
785 0 : write(msg, '(a,a,a,a,f14.8,a,a)' )ch10,&
786 0 : ' dfpt_looppert : WARNING -',ch10,&
787 0 : ' Value of dfpt_sciss=',dtset%dfpt_sciss,ch10,&
788 0 : ' Scissor with d/dk calculation : you are using a "naive" approach !'
789 0 : call wrtout([std_out, ab_out], msg)
790 : end if
791 553 : else if(ipert==dtset%natom+2)then
792 194 : write(msg, '(a,i4)' )' Perturbation : homogeneous electric field along direction',idir
793 582 : call wrtout([std_out, ab_out], msg)
794 194 : if( iscf_mod == -3 )then
795 1 : write(msg, '(a,a,a,a,a,a)' )ch10,&
796 1 : ' dfpt_looppert : COMMENT -',ch10,&
797 1 : ' The first-order density is imposed to be zero (iscf=-3).',ch10,&
798 2 : ' This corresponds to a calculation without local fields.'
799 3 : call wrtout([std_out, ab_out], msg)
800 : end if
801 359 : else if(ipert==dtset%natom+5)then
802 8 : write(msg, '(a,i4)' )' Perturbation : homogeneous spin magnetic field along direction',idir
803 24 : call wrtout([std_out, ab_out], msg)
804 351 : else if(ipert==dtset%natom+6)then
805 1 : write(msg, '(a)' )' Perturbation : homogeneous scalar potential'
806 3 : call wrtout([std_out, ab_out], msg)
807 350 : else if(ipert==dtset%natom+10.or.ipert==dtset%natom+11)then
808 179 : call rf2_getidirs(idir,idir1,idir2)
809 179 : if(ipert==dtset%natom+10)then
810 131 : write(msg,'(2(a,i1))') ' Perturbation : 2nd derivative wrt k, idir1 = ',idir1,&
811 262 : ' idir2 = ',idir2
812 : else
813 48 : write(msg,'(2(a,i1),a)') ' Perturbation : 2nd derivative wrt k (idir1 =',idir1,&
814 96 : ') and Efield (idir2 =',idir2,')'
815 : end if
816 537 : call wrtout([std_out, ab_out], msg)
817 179 : if( iscf_mod /= -3 )then
818 14 : write(msg, '(4a)' )ch10,&
819 14 : ' dfpt_looppert : COMMENT -',ch10,&
820 28 : ' In this case, iscf is set to -3 automatically.'
821 42 : call wrtout([std_out, ab_out], msg)
822 14 : iscf_mod=-3
823 : end if
824 179 : if( abs(dtset%dfpt_sciss) > 1.0d-8 )then
825 0 : write(msg, '(a,a,a,a,f14.8,a,a)' )ch10,&
826 0 : ' dfpt_looppert : WARNING -',ch10,&
827 0 : ' Value of dfpt_sciss=',dtset%dfpt_sciss,ch10,&
828 0 : ' Scissor with d/dk calculation : you are using a "naive" approach !'
829 0 : call wrtout([std_out, ab_out], msg)
830 : end if
831 537 : ABI_MALLOC(occ_pert,(dtset%mband*nkpt*dtset%nsppol))
832 66339 : occ_pert(:) = occ(:) - occ(1)
833 66518 : maxocc = maxval(abs(occ_pert))
834 179 : if (maxocc>1.0d-6.and.abs(maxocc-occ(1))>1.0d-6) then ! True if non-zero occupation numbers are not equal
835 0 : write(msg, '(3a)' ) ' ipert=natom+10 or 11 does not work for a metallic system.',ch10,&
836 0 : ' This perturbation will not be computed.'
837 0 : ABI_WARNING(msg)
838 0 : ABI_FREE(occ_pert)
839 0 : cycle
840 : end if
841 358 : ABI_FREE(occ_pert)
842 171 : else if(ipert>dtset%natom+11.and.ipert<=2*dtset%natom+11)then
843 3 : write(msg, '(a,i4,a,i4)' )' Perturbation : local spin magnetic field from atom', ipert-dtset%natom-11,&
844 6 : & ' along direction', idir
845 9 : call wrtout([std_out, ab_out], msg)
846 168 : else if(ipert>2*dtset%natom+11 .or. ipert<=0 )then
847 : write(msg, '(a,i0,3a)' ) &
848 0 : 'ipert= ',ipert,' is outside the [1,2*natom+11] interval.',ch10,&
849 0 : 'This perturbation is not (yet) allowed.'
850 0 : ABI_BUG(msg)
851 : end if
852 :
853 2342 : if (abs(dtset%magpen) > tol6) then
854 11 : if (dtset%magpen<zero) then
855 0 : write(msg, '(a,3f10.6)' ) ' Homogeneous spin magnetic penalty term is applied.'
856 0 : call wrtout([std_out, ab_out],msg)
857 : else
858 11 : write(msg, '(a,3f10.6)' ) ' Local spin magnetic penalty term from atomic spheres is applied.'
859 33 : call wrtout([std_out, ab_out],msg)
860 : end if
861 : end if
862 :
863 : ! Initialize the diverse parts of energy :
864 2342 : eew=zero ; evdw=zero ; efrloc=zero ; efrnl=zero ; efrx1=zero ; efrx2=zero
865 2342 : efrhar=zero ; efrkin=zero
866 2342 : if(ipert<=dtset%natom)then
867 1353 : eew=dyew(1,idir,ipert,idir,ipert)
868 1353 : if (usevdw==1) evdw=dyvdw(1,idir,ipert,idir,ipert)
869 1353 : efrloc=dyfrlo(idir,idir,ipert)
870 1353 : if (dyfr_nondiag==0) efrnl=dyfrnl(1,idir,idir,ipert,1)
871 1353 : if (dyfr_nondiag/=0) efrnl=dyfrnl(1,idir,idir,ipert,ipert)
872 1353 : efrx1=dyfrx1(1,idir,ipert,idir,ipert)
873 1353 : efrx2=dyfrx2(idir,idir,ipert)
874 989 : else if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then
875 : ! istr = 1,2,...,6 and indicates the cartesian strain component
876 168 : if(ipert==dtset%natom+4) istr=idir+3
877 168 : eii=eltcore(istr,istr)
878 168 : eew=elteew(istr,istr)
879 168 : if (usevdw==1) evdw=eltvdw(istr,istr)
880 168 : efrhar=eltfrhar(istr,istr)
881 168 : efrkin=eltfrkin(istr,istr)
882 168 : efrloc=eltfrloc(istr,istr)
883 168 : efrnl=eltfrnl(istr,istr)
884 168 : efrx1=eltfrxc(istr,istr)
885 : end if
886 :
887 : ! Determine the subset of symmetry operations (nsym1 operations)
888 : ! that leaves the perturbation invariant, and initialize corresponding arrays
889 : ! symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW)..
890 7026 : ABI_MALLOC(symaf1_tmp,(nsym))
891 7026 : ABI_MALLOC(symrl1_tmp,(3,3,nsym))
892 7026 : ABI_MALLOC(tnons1_tmp,(3,nsym))
893 :
894 : if (dtset%prepanl/=1.and.&
895 : & dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.&
896 2342 : & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17) then
897 : call littlegroup_pert(gprimd,idir,indsym,ab_out,ipert,dtset%natom,nsym,nsym1,2,&
898 2087 : & dtset%symafm,symaf1_tmp,symq,symrec,dtset%symrel,symrl1_tmp,0,dtset%tnons,tnons1_tmp)
899 : else
900 255 : nsym1 = 1
901 255 : symaf1_tmp(1) = 1
902 3315 : symrl1_tmp(:,:,1) = dtset%symrel(:,:,1)
903 1020 : tnons1_tmp(:,1) = 0_dp
904 : end if
905 9368 : ABI_MALLOC(indsy1,(4,nsym1,dtset%natom))
906 7026 : ABI_MALLOC(symrc1,(3,3,nsym1))
907 7026 : ABI_MALLOC(symaf1,(nsym1))
908 4684 : ABI_MALLOC(symrl1,(3,3,nsym1))
909 7026 : ABI_MALLOC(tnons1,(3,nsym1))
910 6389 : symaf1(1:nsym1)=symaf1_tmp(1:nsym1)
911 54953 : symrl1(:,:,1:nsym1)=symrl1_tmp(:,:,1:nsym1)
912 18530 : tnons1(:,1:nsym1)=tnons1_tmp(:,1:nsym1)
913 2342 : ABI_FREE(symaf1_tmp)
914 2342 : ABI_FREE(symrl1_tmp)
915 2342 : ABI_FREE(tnons1_tmp)
916 :
917 : ! Set up corresponding symmetry data
918 9368 : ABI_MALLOC(irrzon1,(dtset%nfft**(1-1/nsym1),2,(nspden/dtset%nsppol)-3*(nspden/4)))
919 9368 : ABI_MALLOC(phnons1,(2,dtset%nfft**(1-1/nsym1),(nspden/dtset%nsppol)-3*(nspden/4)))
920 : call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,nspden,dtset%nsppol,&
921 2342 : & nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred)
922 2342 : if (psps%usepaw==1) then
923 : ! Allocate/initialize only zarot in pawang1 datastructure
924 397 : call pawang_init(pawang1,0,0,pawang%l_max-1,0,0,nsym1,0,0,0,0)
925 397 : call setsym_ylm(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot)
926 : end if
927 :
928 : ! Initialize k+q array
929 7026 : ABI_MALLOC(kpq,(3,nkpt))
930 2342 : if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
931 77784 : kpq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran
932 : else
933 155071 : do ikpt=1,nkpt
934 613762 : kpq(:,ikpt)=dtset%qptn(:)+dtset%kptns(:,ikpt)
935 : end do
936 : end if
937 : ! In case wf1 at +q and -q are not related by time inversion symmetry, compute k-q as well for initializations
938 2342 : if (.not.kramers_deg) then
939 33 : ABI_MALLOC(kmq,(3,nkpt))
940 11 : if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
941 0 : kmq(:,1:nkpt)=dtset%kptns(:,1:nkpt) ! Do not modify, needed for gfortran
942 : else
943 131 : do ikpt=1,nkpt
944 491 : kmq(:,ikpt)=-dtset%qptn(:)+dtset%kptns(:,ikpt) ! kmq <= k-q
945 : end do
946 : end if
947 : end if
948 :
949 : ! Determine the subset of k-points needed in the "reduced Brillouin zone" and initialize other quantities
950 7026 : ABI_MALLOC(indkpt1_tmp,(nkpt))
951 7026 : ABI_MALLOC(wtk_folded,(nkpt))
952 7026 : ABI_MALLOC(bz2ibz_smap, (6, nkpt))
953 174643 : indkpt1_tmp(:)=0 ; optthm=0
954 2342 : timrev_pert=timrev
955 2342 : if(dtset%ieig2rf>0) then
956 494 : timrev_pert=0
957 : call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,&
958 494 : & 1,symrc1,timrev_pert,dtset%wtk,wtk_folded, bz2ibz_smap, xmpi_comm_self)
959 : else
960 : ! For the time being, the time reversal symmetry is not used
961 : ! for ddk, elfd, mgfd, sclrpt perturbations.
962 : timrev_pert=timrev
963 : if(ipert==dtset%natom+1.or.ipert==dtset%natom+2.or.&
964 : & ipert==dtset%natom+10.or.ipert==dtset%natom+11.or. &
965 : & dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or. &
966 : & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17.or. &
967 : & ipert==dtset%natom+5.or.ipert==dtset%natom+6.or. &
968 1848 : & (ipert>dtset%natom+11.and.ipert<=2*dtset%natom+11).or. &
969 828 : & dtset%prtfull1wf==1) timrev_pert=0
970 1848 : timrev_kpt = timrev_pert
971 :
972 : !MR: Modified to agree with longwave driver
973 1848 : if(dtset%prepalw/=0) then
974 98 : if (dtset%kptopt==2) timrev_pert=1
975 98 : if (dtset%kptopt==3) timrev_pert=0
976 98 : timrev_kpt = timrev_pert
977 : !MR tmp: this has to be removed if perturbation-dependent spatial symmetries are
978 : !implemented in the spatial-dispersion routines
979 98 : nsym1=1
980 :
981 98 : if (dtset%rfstrs/=0.and.dtset%rfstrs_ref==0) then
982 : write(msg,'(9a)')&
983 0 : 'If the outputs of this strain response function calculation are to be',ch10,&
984 0 : 'subsequently used as inputs of a longwave calculation the same energy',ch10,&
985 0 : 'reference has to be used in both cases. Otherwise wrong spatial-dispersion',ch10,&
986 0 : 'coefficients will be obtained',ch10,&
987 0 : 'Action: Put rfstrs_ref=1 '
988 0 : ABI_WARNING(msg)
989 : end if
990 : end if
991 :
992 : ! The time reversal symmetry is not used for the BZ sampling when kptopt=3 or 4
993 1848 : if (dtset%kptopt==3.or.dtset%kptopt==4) timrev_kpt = 0
994 : call symkpt(0,gmet,indkpt1_tmp,ab_out,dtset%kptns,nkpt,nkpt_rbz,&
995 1848 : nsym1,symrc1,timrev_kpt,dtset%wtk,wtk_folded, bz2ibz_smap, xmpi_comm_self)
996 : end if
997 :
998 2342 : write (msg, "(a,i10,a,i10,a)") " dfpt_looppert : the # of k-points has been reduced from ", &
999 4684 : nkpt, " to ", nkpt_rbz, " using the little group "
1000 2342 : call wrtout(std_out, msg)
1001 :
1002 7026 : ABI_MALLOC(doccde_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol))
1003 7026 : ABI_MALLOC(indkpt1,(nkpt_rbz))
1004 4684 : ABI_MALLOC(istwfk_rbz,(nkpt_rbz))
1005 7026 : ABI_MALLOC(kpq_rbz,(3,nkpt_rbz))
1006 2342 : if (.not.kramers_deg) then
1007 22 : ABI_MALLOC(kmq_rbz,(3,nkpt_rbz))
1008 : end if
1009 4684 : ABI_MALLOC(kpt_rbz,(3,nkpt_rbz))
1010 7026 : ABI_MALLOC(nband_rbz,(nkpt_rbz*dtset%nsppol))
1011 4684 : ABI_MALLOC(occ_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol))
1012 4684 : ABI_MALLOC(occ_disk,(dtset%mband*nkpt_rbz*dtset%nsppol))
1013 7026 : ABI_MALLOC(wtk_rbz,(nkpt_rbz))
1014 150028 : indkpt1(:)=indkpt1_tmp(1:nkpt_rbz)
1015 150028 : do ikpt=1,nkpt_rbz
1016 147686 : istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt))
1017 590744 : kpq_rbz(:,ikpt)=kpq(:,indkpt1(ikpt))
1018 590744 : kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt))
1019 150028 : wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt))
1020 : end do
1021 2342 : if (.not.kramers_deg) then
1022 131 : do ikpt=1,nkpt_rbz
1023 491 : kmq_rbz(:,ikpt)=kmq(:,indkpt1(ikpt))
1024 : end do
1025 : end if
1026 2342 : ABI_FREE(indkpt1_tmp)
1027 2342 : ABI_FREE(wtk_folded)
1028 :
1029 : ! Transfer occ to occ_rbz and doccde to doccde_rbz :
1030 : ! this is a more delicate issue
1031 : ! NOTE : this takes into account that indkpt1 is ordered
1032 : ! MG: What about using occ(band,kpt,spin) ???
1033 : ! GA: Would be better indeed, but I think indices neededed to be consistent
1034 : ! with second derivative of eigenvalues, which has spin index wrapped up
1035 : ! to avoid arrays with rank larger than 7.
1036 2342 : bdtot_index=0;bdtot1_index=0
1037 4851 : do isppol=1,dtset%nsppol
1038 2509 : ikpt1=1
1039 179066 : do ikpt=1,nkpt
1040 174215 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
1041 : ! Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1)
1042 174215 : if(ikpt1/=nkpt_rbz+1)then
1043 167873 : if(ikpt==indkpt1(ikpt1))then
1044 149450 : nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k
1045 987050 : occ_rbz(1+bdtot1_index:nband_k+bdtot1_index) = occ(1+bdtot_index:nband_k+bdtot_index)
1046 987050 : doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index) = doccde(1+bdtot_index:nband_k+bdtot_index)
1047 149450 : ikpt1=ikpt1+1
1048 149450 : bdtot1_index=bdtot1_index+nband_k
1049 : end if
1050 : end if
1051 176724 : bdtot_index=bdtot_index+nband_k
1052 : end do
1053 : end do
1054 :
1055 : ! Compute maximum number of planewaves at k
1056 2342 : call timab(142,1,tsec)
1057 2342 : call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpt_rbz,mpi_enreg,mpw,nkpt_rbz)
1058 2342 : call timab(142,2,tsec)
1059 :
1060 : ! Allocate some k-dependent arrays at k
1061 7026 : ABI_MALLOC(npwarr,(nkpt_rbz))
1062 4684 : ABI_MALLOC(npwtot,(nkpt_rbz))
1063 :
1064 : ! Determine distribution of k-points/bands over MPI processes
1065 2342 : if (allocated(mpi_enreg%my_kpttab)) then
1066 716 : ABI_FREE(mpi_enreg%my_kpttab)
1067 : end if
1068 7026 : ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz))
1069 : if(xmpi_paral==1) then
1070 11710 : ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,dtset%mband,dtset%nsppol))
1071 2342 : call distrb2(dtset%mband,mband_mem_rbz,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,dtset%nsppol,mpi_enreg)
1072 : else
1073 : mpi_enreg%my_kpttab(:)=(/(ii,ii=1,nkpt_rbz)/)
1074 : end if
1075 150028 : my_nkpt_rbz=maxval(mpi_enreg%my_kpttab)
1076 2342 : mkmem_rbz =my_nkpt_rbz ; mkqmem_rbz=my_nkpt_rbz ; mk1mem_rbz=my_nkpt_rbz
1077 : ABI_UNUSED((/mkmem,mk1mem,mkqmem/))
1078 :
1079 2342 : call initmpi_band(mkmem_rbz,mpi_enreg,nband_rbz,nkpt_rbz,dtset%nsppol)
1080 :
1081 : ! given number of reduced kpt, store distribution of bands across procs
1082 11710 : ABI_MALLOC(distrb_flags,(nkpt_rbz,dtset%mband,dtset%nsppol))
1083 866088 : distrb_flags = (mpi_enreg%proc_distrb == mpi_enreg%me_kpt)
1084 :
1085 : ! Set up the basis sphere of planewaves at k
1086 7026 : ABI_MALLOC(kg,(3,mpw*mkmem_rbz))
1087 2342 : call timab(143,1,tsec)
1088 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg,&
1089 : kpt_rbz,mkmem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,&
1090 2342 : mpw,npwarr,npwtot,dtset%nsppol)
1091 2342 : call timab(143,2,tsec)
1092 :
1093 : ! Set up the spherical harmonics (Ylm) at k
1094 2342 : useylmgr=0; option=0 ; nylmgr=0
1095 2342 : if (psps%useylm==1.and. &
1096 : (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. &
1097 : (psps%usepaw==1.and.ipert==dtset%natom+2))) then
1098 318 : useylmgr=1; option=1 ; nylmgr=3
1099 2024 : else if (psps%useylm==1.and.(ipert==dtset%natom+10.or.ipert==dtset%natom+11)) then
1100 179 : useylmgr=1; option=2 ; nylmgr=9
1101 : end if
1102 9368 : ABI_MALLOC(ylm,(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
1103 11710 : ABI_MALLOC(ylmgr,(mpw*mkmem_rbz,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
1104 2342 : if (psps%useylm==1) then
1105 : call initylmg(gprimd,kg,kpt_rbz,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,nband_rbz,nkpt_rbz,&
1106 708 : npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr)
1107 : end if
1108 :
1109 : ! Set up occupations for this perturbation
1110 2342 : if (dtset%ieig2rf>0) then
1111 494 : if (.not.allocated(istwfk_pert)) then
1112 336 : ABI_MALLOC(istwfk_pert,(nkpt,3,mpert))
1113 252 : ABI_MALLOC(occ_pert,(dtset%mband*nkpt*dtset%nsppol))
1114 58134 : istwfk_pert(:,:,:)=0 ; occ_pert(:)=zero
1115 : end if
1116 9220 : istwfk_pert(:,idir,ipert)=istwfk_rbz(:)
1117 91022 : occ_pert(:)=occ_rbz(:)
1118 : end if
1119 2342 : if (dtset%efmas>0) then
1120 51 : if (.not.allocated(istwfk_pert)) then
1121 68 : ABI_MALLOC(istwfk_pert,(nkpt,3,mpert))
1122 1287 : istwfk_pert(:,:,:)=0
1123 : end if
1124 123 : istwfk_pert(:,idir,ipert)=istwfk_rbz(:)
1125 : end if
1126 :
1127 : ! Print a separator in output file
1128 2342 : write(msg,'(3a)')ch10,'--------------------------------------------------------------------------------',ch10
1129 2342 : call wrtout(ab_out,msg)
1130 :
1131 : ! Initialize band structure datatype at k
1132 : ! TODO: occ_rbz and doccde_rbz are overwritten here, though they could be set from the input file, above
1133 : ! need to check if this is correct (or indifferent) in all cases
1134 151792 : bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol))
1135 7026 : ABI_MALLOC(eigen0,(bantot_rbz))
1136 839942 : eigen0(:)=zero
1137 : call ebands_k%init(bantot_rbz, dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
1138 : doccde_rbz,eigen0,istwfk_rbz,kpt_rbz,&
1139 : nband_rbz,nkpt_rbz,npwarr,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,&
1140 : dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
1141 2342 : dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1142 2342 : ABI_FREE(eigen0)
1143 :
1144 : ! Initialize header, update it with evolving variables
1145 2342 : gscase=0 ! A GS WF file is read
1146 : call hdr0%init(ebands_k,codvsn,dtset,pawtab,gscase,psps,wvl%descr,&
1147 2342 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1148 :
1149 : call hdr0%update(bantot_rbz,etotal,fermie,fermie,& ! CP: duplicated fermie to fit new def of hdr_update
1150 : residm,rprimd,occ_rbz,pawrhoij_pert,xred,dtset%amu_orig(:,1),&
1151 2342 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1152 :
1153 : ! Initialize GS wavefunctions at k
1154 2342 : ireadwf0=1; formeig=0 ; ask_accurate=1 ; optorth=0
1155 2342 : mcg=mpw*dtset%nspinor*mband_mem_rbz*mkmem_rbz*dtset%nsppol
1156 :
1157 : call wrtout(std_out, sjoin(" Memory required for psi0_k, psi0_kq, psi1_kq: ", &
1158 2342 : ftoa(3 * two * mcg * dp * b2Mb, fmt="f8.1"), "[Mb] <<< MEM"))
1159 :
1160 2342 : if (one*mpw*dtset%nspinor*mband_mem_rbz*mkmem_rbz*dtset%nsppol > huge(1)) then
1161 : write (msg,'(4a, 5(a,i0), 2a)')&
1162 0 : "Default integer is not wide enough to store the size of the GS wavefunction array (WF0, mcg).",ch10,&
1163 0 : "Action: increase the number of processors. Consider also OpenMP threads.",ch10,&
1164 0 : "nspinor: ",dtset%nspinor, "mpw: ",mpw, "mband_mem_rbz: ",mband_mem_rbz, "mkmem_rbz: ",&
1165 0 : mkmem_rbz, "nsppol: ",dtset%nsppol,ch10,&
1166 0 : 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode'
1167 0 : ABI_ERROR(msg)
1168 : end if
1169 7026 : ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr)
1170 :
1171 7026 : ABI_MALLOC(eigen0,(dtset%mband*nkpt_rbz*dtset%nsppol))
1172 2342 : call timab(144,1,tsec)
1173 2342 : call pstat_proc%print(_PSTAT_ARGS_)
1174 :
1175 : ! Initialize the wave function type and read GS WFK
1176 : call wfk_read_my_kptbands(dtfil%fnamewffk, distrb_flags, spacecomm, dtset%ecut*(dtset%dilatmx)**2,&
1177 : formeig, istwfk_rbz, kpt_rbz, mcg, dtset%mband, mband_mem_rbz, mkmem_rbz, mpw,&
1178 : dtset%natom, nkpt_rbz, npwarr, dtset%nspinor, dtset%nsppol, dtset%usepaw,&
1179 2342 : cg, eigen=eigen0, occ=occ_disk)
1180 :
1181 2342 : call timab(144,2,tsec)
1182 :
1183 : ! Update energies GS energies at k
1184 2342 : call put_eneocc_vect(ebands_k, "eig", eigen0)
1185 :
1186 : ! PAW: compute on-site projections of GS wavefunctions (cprj) (and derivatives) at k
1187 2342 : ncpgr=0
1188 2342 : ABI_MALLOC(cprj,(0,0))
1189 2342 : if (psps%usepaw==1) then
1190 397 : ncpgr=3 ! Valid for ipert<=natom (phonons), ipert=natom+2 (elec. field)
1191 : ! or for ipert==natom+10,11
1192 397 : if (ipert==dtset%natom+1) then
1193 129 : if (dtset%orbmag.NE.0) then
1194 : ncpgr=3
1195 : else
1196 84 : ncpgr=1
1197 : end if
1198 : end if
1199 397 : if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) ncpgr=1
1200 397 : if (usecprj==1) then
1201 : ! distribute cprj by band as well
1202 : ! NB: currently nsppol=2 is distributed in data (0s saved for other spin)
1203 : ! but not in memory: all procs have nsppol 2 below
1204 397 : mcprj=dtset%nspinor*mband_mem_rbz*mkmem_rbz*dtset%nsppol
1205 : !mcprj=dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol
1206 397 : ABI_FREE(cprj)
1207 243120 : ABI_MALLOC(cprj,(dtset%natom,mcprj))
1208 397 : call pawcprj_alloc(cprj,ncpgr,dimcprj_srt)
1209 397 : if (ipert<=dtset%natom) then
1210 114 : choice=2; iorder_cprj=0; idir0=0
1211 283 : else if (ipert==dtset%natom+1) then
1212 129 : if (dtset%orbmag.NE.0) then
1213 45 : choice=5; iorder_cprj=0; idir0=0
1214 : else
1215 84 : choice=5; iorder_cprj=0; idir0=idir
1216 : end if
1217 154 : else if (ipert==dtset%natom+2) then
1218 32 : choice=5; iorder_cprj=0; idir0=0
1219 122 : else if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) then
1220 72 : choice=3; iorder_cprj=0; idir0=istr
1221 50 : else if (ipert==dtset%natom+10.or.ipert==dtset%natom+11) then
1222 50 : choice=5; iorder_cprj=0; idir0=0 ! Compute all first derivatives
1223 : else
1224 0 : choice=1; iorder_cprj=0; idir0=0
1225 : end if
1226 : call ctocprj(atindx,cg,choice,cprj,gmet,gprimd,-1,idir0,iorder_cprj,istwfk_rbz,&
1227 : & kg,kpt_rbz,mcg,mcprj,dtset%mgfft,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,&
1228 : & dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,&
1229 : & npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,&
1230 397 : & rmet,dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr)
1231 : end if
1232 : end if
1233 :
1234 : ! Compute maximum number of planewaves at k+q
1235 : ! Will be useful for both GS wfs at k+q and RF wavefunctions
1236 2342 : call timab(143,1,tsec)
1237 2342 : call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpq_rbz,mpi_enreg,mpw1,nkpt_rbz)
1238 2342 : if (.not.kramers_deg) then
1239 11 : call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kmq_rbz,mpi_enreg,mpw1_mq,nkpt_rbz)
1240 : !number of plane waves at k+q and k-q should be in principle the same to reconstruct rhor1_pq (?)
1241 11 : mpw1=max(mpw1,mpw1_mq)
1242 11 : mpw1_mq=mpw1
1243 : else
1244 2331 : mpw1_mq=0
1245 : end if
1246 2342 : call timab(143,2,tsec)
1247 :
1248 : ! Allocate some arrays at k+q
1249 7026 : ABI_MALLOC(kg1,(3,mpw1*mk1mem_rbz))
1250 7026 : ABI_MALLOC(npwar1,(nkpt_rbz))
1251 4684 : ABI_MALLOC(npwtot1,(nkpt_rbz))
1252 : ! In case Kramers degeneracy is broken, do the same for k-q
1253 2342 : if (.not.kramers_deg) then
1254 33 : ABI_MALLOC(kg1_mq,(3,mpw1_mq*mk1mem_rbz))
1255 22 : ABI_MALLOC(npwar1_mq,(nkpt_rbz))
1256 22 : ABI_MALLOC(npwtot1_mq,(nkpt_rbz))
1257 : end if
1258 :
1259 : ! Set up the basis sphere of planewaves at k+q
1260 : ! Will be useful for both GS wfs at k+q and RF wavefunctions
1261 2342 : call timab(142,1,tsec)
1262 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1,&
1263 : kpq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1,&
1264 2342 : npwar1,npwtot1,dtset%nsppol)
1265 2342 : if (.not.kramers_deg) then
1266 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kg1_mq,&
1267 : kmq_rbz,mk1mem_rbz,nband_rbz,nkpt_rbz,'PERS',mpi_enreg,mpw1_mq,&
1268 11 : npwar1_mq,npwtot1_mq,dtset%nsppol)
1269 : end if
1270 2342 : call timab(142,2,tsec)
1271 :
1272 : ! Set up the spherical harmonics (Ylm) at k+q
1273 2342 : useylmgr1=0; option=0 ; nylmgr1=0
1274 2342 : if (psps%useylm==1.and. &
1275 : (ipert==dtset%natom+1.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or. &
1276 : (psps%usepaw==1.and.ipert==dtset%natom+2))) then
1277 318 : useylmgr1=1; option=1; nylmgr1=3
1278 2024 : else if (psps%useylm==1.and.(ipert==dtset%natom+10.or.ipert==dtset%natom+11)) then
1279 179 : useylmgr1=1; option=2; nylmgr1=9
1280 : end if
1281 9368 : ABI_MALLOC(ylm1,(mpw1*mk1mem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
1282 11710 : ABI_MALLOC(ylmgr1,(mpw1*mk1mem_rbz,nylmgr1,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
1283 2342 : if (psps%useylm==1) then
1284 : call initylmg(gprimd,kg1,kpq_rbz,mk1mem_rbz,mpi_enreg,psps%mpsang,mpw1,nband_rbz,nkpt_rbz,&
1285 708 : npwar1,dtset%nsppol,option,rprimd,ylm1,ylmgr1)
1286 : end if
1287 2342 : if (.not.kramers_deg) then
1288 44 : ABI_MALLOC(ylm1_mq,(mpw1_mq*mk1mem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
1289 55 : ABI_MALLOC(ylmgr1_mq,(mpw1_mq*mk1mem_rbz,nylmgr1,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
1290 11 : if (psps%useylm==1) then
1291 0 : if (abs(sum(dtset%qptn(1:3)**2)) < 1.d-14) then
1292 0 : ylm1_mq(:,:)=ylm1(:,:)
1293 0 : ylmgr1_mq(:,:,:)=ylmgr1(:,:,:)
1294 : else
1295 : call initylmg(gprimd,kg1_mq,kmq_rbz,mk1mem_rbz,mpi_enreg,psps%mpsang,mpw1_mq,nband_rbz,nkpt_rbz,&
1296 0 : & npwar1_mq,dtset%nsppol,option,rprimd,ylm1_mq,ylmgr1_mq)
1297 : end if
1298 : end if
1299 : end if
1300 :
1301 : ! Print a separator in output file
1302 2342 : write(msg, '(a,a)' )'--------------------------------------------------------------------------------',ch10
1303 2342 : call wrtout(ab_out,msg)
1304 :
1305 : ! Initialize band structure datatype at k+q
1306 4684 : ABI_MALLOC(eigenq,(bantot_rbz))
1307 7026 : ABI_MALLOC(doccde_tmp,(dtset%mband*nkpt_rbz*dtset%nsppol))
1308 839942 : eigenq(:)=zero
1309 : call ebands_kq%init(bantot_rbz, dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
1310 : & doccde_tmp,eigenq,istwfk_rbz,kpq_rbz,&
1311 : & nband_rbz,nkpt_rbz,npwar1,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,&
1312 : & dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
1313 2342 : & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1314 2342 : if (.not.kramers_deg) then
1315 4427 : eigenq(:)=zero
1316 : call ebands_kmq%init(bantot_rbz, dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
1317 : & doccde_tmp,eigenq,istwfk_rbz,kmq_rbz,&
1318 : & nband_rbz,nkpt_rbz,npwar1_mq,dtset%nsppol,dtset%nspinor,dtset%tphysel,dtset%tsmear,dtset%occopt,occ_rbz,wtk_rbz,&
1319 : & dtset%cellcharge(1), dtset%kptopt, dtset%kptrlatt_orig, dtset%nshiftk_orig, dtset%shiftk_orig, &
1320 11 : & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1321 : end if
1322 2342 : ABI_FREE(eigenq)
1323 2342 : ABI_FREE(doccde_tmp)
1324 :
1325 : ! Initialize header
1326 : call hdr%init(ebands_kq,codvsn,dtset,pawtab,pertcase,psps,wvl%descr, &
1327 2342 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
1328 2342 : if (.not.kramers_deg) then
1329 : !MR: this is to not interfere with the second-order derivatives
1330 : !neither with the local spin perturbation
1331 11 : pertcase_mq=pertcase+(2*dtset%natom+11)*3
1332 : call hdr_mq%init(ebands_kmq,codvsn,dtset,pawtab,pertcase_mq,psps,wvl%descr, &
1333 11 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
1334 : end if
1335 :
1336 : ! Initialize wavefunctions at k+q
1337 : ! MG: Here it is possible to avoid the extra reading if the same k mesh can be used.
1338 2342 : ireadwf0=1 ; formeig=0 ; ask_accurate=1 ; optorth=0
1339 2342 : mcgq=mpw1*dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol
1340 : !SPr: verified until here, add mcgq for -q
1341 2342 : if (one*mpw1*dtset%nspinor*dtset%mband*mkqmem_rbz*dtset%nsppol > huge(1)) then
1342 : write (msg,'(4a, 5(a,i0), 2a)')&
1343 0 : "Default integer is not wide enough to store the size of the GS wavefunction array (WFKQ, mcgq).",ch10,&
1344 0 : "Action: increase the number of processors. Consider also OpenMP threads.",ch10,&
1345 0 : "nspinor: ",dtset%nspinor, "mpw1: ",mpw1, "mband: ",dtset%mband, "mkqmem_rbz: ",&
1346 0 : mkqmem_rbz, "nsppol: ",dtset%nsppol,ch10,&
1347 0 : 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode'
1348 0 : ABI_ERROR(msg)
1349 : end if
1350 7026 : ABI_MALLOC_OR_DIE(cgq,(2,mcgq), ierr)
1351 :
1352 7026 : ABI_MALLOC(eigenq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1353 2342 : if (.not.kramers_deg) then
1354 : !ABI_MALLOC_OR_DIE(cg_pq,(2,mcgq), ierr)
1355 : !ABI_MALLOC(eigen_pq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1356 11 : mcgmq=mpw1_mq*dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol
1357 33 : ABI_MALLOC_OR_DIE(cg_mq,(2,mcgmq), ierr)
1358 :
1359 33 : ABI_MALLOC(eigen_mq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1360 : end if
1361 2342 : call pstat_proc%print(_PSTAT_ARGS_)
1362 :
1363 : !if (sum(dtset%qptn(1:3)**2)>=1.d-14) then ! non-zero q
1364 : !TODO: for many other q this should be avoidable, in principle all if qptrlatt is a subgrid of kprtlatt
1365 : ! Or at very least make a pointer instead of a full copy!!!
1366 9368 : if (dtfil%fnamewffq == dtfil%fnamewffk .and. sum(dtset%qptn(1:3)**2) < 1.d-14) then
1367 1533 : call wrtout(std_out, " qpt is Gamma --> psi_k+q initialized from psi_k in memory")
1368 160919646 : cgq = cg
1369 407504 : eigenq = eigen0
1370 : else
1371 809 : call timab(144,1,tsec)
1372 : call wfk_read_my_kptbands(dtfil%fnamewffq, distrb_flags, spacecomm, dtset%ecut*(dtset%dilatmx)**2,&
1373 : formeig, istwfk_rbz, kpq_rbz, mcgq, dtset%mband, mband_mem_rbz, mkqmem_rbz, mpw1,&
1374 : dtset%natom, nkpt_rbz, npwar1, dtset%nspinor, dtset%nsppol, dtset%usepaw,&
1375 809 : cgq, eigen=eigenq, occ=occ_disk)
1376 809 : call timab(144,2,tsec)
1377 : end if
1378 :
1379 2342 : if (.not.kramers_deg) then
1380 44 : if (dtfil%fnamewffq == dtfil%fnamewffk .and. sum(dtset%qptn(1:3)**2) < 1.d-14) then
1381 4 : call wrtout(std_out, " qpt is Gamma, psi_k-q initialized from psi_k in memory")
1382 1608776 : cg_mq = cg
1383 2024 : eigen_mq = eigen0
1384 : else
1385 7 : call timab(144,1,tsec)
1386 7 : fnamewffmq_=dtfil%fnamewffq
1387 7 : if (dtset%getwfmq/=0) fnamewffmq_=dtfil%fnamewffmq
1388 : call wfk_read_my_kptbands(fnamewffmq_, distrb_flags, spacecomm,dtset%ecut*(dtset%dilatmx)**2, &
1389 : formeig, istwfk_rbz, kmq_rbz, mcgmq, dtset%mband, mband_mem_rbz, mkqmem_rbz, mpw1_mq,&
1390 : dtset%natom, nkpt_rbz, npwar1_mq, dtset%nspinor, dtset%nsppol, dtset%usepaw,&
1391 7 : cg_mq, eigen=eigen_mq, occ=occ_disk)
1392 7 : call timab(144,2,tsec)
1393 : end if
1394 : end if
1395 2342 : ABI_FREE(occ_disk)
1396 :
1397 : ! Update energies GS energies at k + q
1398 2342 : call put_eneocc_vect(ebands_kq, "eig", eigenq)
1399 2342 : if (.not.kramers_deg) call put_eneocc_vect(ebands_kmq, "eig", eigen_mq)
1400 :
1401 : ! PAW: compute on-site projections of GS wavefunctions (cprjq) (and derivatives) at k+q
1402 2342 : ABI_MALLOC(cprjq,(0,0))
1403 2342 : if (psps%usepaw==1) then
1404 397 : if (usecprj==1) then
1405 397 : mcprjq=dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol
1406 : !TODO: distribute cprj by band as well?
1407 : !mcprjq=dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol
1408 397 : ABI_FREE(cprjq)
1409 243120 : ABI_MALLOC(cprjq,(dtset%natom,mcprjq))
1410 397 : call pawcprj_alloc(cprjq,0,dimcprj_srt)
1411 1588 : if (ipert<=dtset%natom.and.(sum(dtset%qptn(1:3)**2)>=1.d-14)) then ! phonons at non-zero q
1412 36 : choice=1 ; iorder_cprj=0 ; idir0=0
1413 : call ctocprj(atindx,cgq,choice,cprjq,gmet,gprimd,-1,idir0,0,istwfk_rbz,&
1414 : & kg1,kpq_rbz,mcgq,mcprjq,dtset%mgfft,mkqmem_rbz,mpi_enreg,psps%mpsang,mpw1,&
1415 : & dtset%natom,nattyp,nband_rbz,dtset%natom,dtset%ngfft,nkpt_rbz,dtset%nloalg,&
1416 : & npwar1,dtset%nspinor,dtset%nsppol,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,&
1417 36 : & psps,rmet,dtset%typat,ucvol,dtfil%unpawq,xred,ylm1,ylmgr1)
1418 361 : else if (mcprjq>0) then
1419 361 : call pawcprj_copy(cprj,cprjq)
1420 : end if
1421 : end if
1422 : end if
1423 :
1424 : ! ===== Report on eigenq values
1425 2342 : if (dtset%ieig2rf>0.and.icase==ipert_cnt) then
1426 15224 : eigen0_pert(:) = eigen0(:)
1427 15224 : eigenq_pert(:) = eigenq(:)
1428 15224 : occ_rbz_pert(:) = occ_rbz(:)
1429 : end if
1430 2342 : if (dtset%efmas>0.and.icase==ipert_cnt) then
1431 355 : eigen0_pert(:) = eigen0(:)
1432 : end if
1433 : !call wrtout(std_out,ch10//' dfpt_looppert: eigenq array',"COLL")
1434 2342 : nkpt_eff=nkpt
1435 2342 : if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. nkpt>nkpt_max ) nkpt_eff=nkpt_max
1436 2342 : band_index=0
1437 4851 : do isppol=1,dtset%nsppol
1438 154301 : do ikpt=1,nkpt_rbz
1439 149450 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
1440 149450 : if(ikpt<=nkpt_eff)then
1441 22618 : write(msg, '(a,i2,a,i5)' )' isppol=',isppol,', k point number',ikpt
1442 22618 : call wrtout(std_out,msg)
1443 67758 : do iband=1,nband_k,4
1444 45140 : write(msg, '(a,4es16.6)')' ',eigenq(iband+band_index:min(iband+3,nband_k)+band_index)
1445 67758 : call wrtout(std_out,msg)
1446 : end do
1447 126832 : else if(ikpt==nkpt_eff+1)then
1448 0 : write(msg,'(a,a)' )' respfn : prtvol=0, 1 or 2, stop printing eigenq.',ch10
1449 0 : call wrtout(std_out,msg)
1450 : end if
1451 151959 : band_index=band_index+nband_k
1452 : end do
1453 : end do
1454 :
1455 : ! Generate occupation numbers for the reduced BZ at k+q
1456 7026 : ABI_MALLOC(docckqde,(dtset%mband*nkpt_rbz*dtset%nsppol))
1457 4684 : ABI_MALLOC(occkq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1458 2342 : if (.not.kramers_deg) then
1459 22 : ABI_MALLOC(docckde_mq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1460 22 : ABI_MALLOC(occk_mq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1461 : end if
1462 :
1463 2342 : if(0<=dtset%occopt .and. dtset%occopt<=2)then
1464 : ! Same occupation numbers at k and k+q (usually, insulating)
1465 537641 : occkq(:)=occ_rbz(:)
1466 537641 : docckqde(:)=zero ! docckqde is irrelevant in this case
1467 1855 : if(.not.kramers_deg) then
1468 4040 : occk_mq(:)=occ_rbz(:)
1469 4040 : docckde_mq(:)=zero
1470 : end if
1471 : else
1472 : ! Metallic occupation numbers
1473 487 : option=1
1474 487 : dosdeltae=zero ! the DOS is not computed with option=1
1475 487 : maxocc=two/(dtset%nspinor*dtset%nsppol)
1476 : call getnel(docckqde,dosdeltae,eigenq,entropy,fermie,fermie,maxocc,dtset%mband,&
1477 : nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occkq,dtset%occopt,option,&
1478 487 : dtset%tphysel,dtset%tsmear,fake_unit,wtk_rbz,1,dtset%nband(1)) ! CP: added 1, nband(0) to fit new def of getnel
1479 : ! Compare nelect at k and nelelect at k+q
1480 : write(msg, '(a,a,a,es16.6,a,es16.6,a)')&
1481 487 : ' dfpt_looppert : total number of electrons, from k and k+q',ch10,&
1482 974 : ' fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.'
1483 1461 : call wrtout([std_out, ab_out], msg)
1484 487 : if (.not.kramers_deg) then
1485 : call getnel(docckde_mq,dosdeltae,eigen_mq,entropy,fermie,fermie,maxocc,dtset%mband,&
1486 : nband_rbz,nelectkq,nkpt_rbz,dtset%nsppol,occk_mq,dtset%occopt,option,&
1487 3 : dtset%tphysel,dtset%tsmear,fake_unit,wtk_rbz,1,dtset%nband(1))
1488 : ! Compare nelect at k and nelelect at k-q
1489 : write(msg, '(a,a,a,es16.6,a,es16.6,a)')&
1490 3 : ' dfpt_looppert : total number of electrons, from k and k-q',ch10,&
1491 6 : ' fully or partially occupied states are',dtset%nelect,' and',nelectkq,'.'
1492 9 : call wrtout([std_out, ab_out], msg)
1493 : end if
1494 : end if
1495 :
1496 : ! Debug message
1497 2342 : if(dtset%prtvol==-level) call wrtout(std_out,'dfpt_looppert: initialisation of q part done.')
1498 :
1499 : ! Initialisation of first-order wavefunctions
1500 2342 : write(msg,'(3a,i4)')' Initialisation of the first-order wave-functions :',ch10,' ireadwf=',dtfil%ireadwf
1501 7026 : call wrtout([std_out, ab_out], msg)
1502 2342 : call appdig(pertcase,dtfil%fnamewff1,fiwf1i)
1503 2342 : call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o)
1504 2342 : if (.not.kramers_deg) then
1505 11 : call appdig(pertcase_mq,dtfil%fnamewff1,fiwf1i_mq)
1506 11 : call appdig(pertcase_mq,dtfil%fnameabo_1wf,fiwf1o_mq)
1507 : end if
1508 :
1509 : ! Allocate 1st-order PAW occupancies (rhoij1)
1510 2342 : if (psps%usepaw==1) then
1511 1963 : ABI_MALLOC(pawrhoij1,(my_natom))
1512 397 : call pawrhoij_nullify(pawrhoij1)
1513 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
1514 397 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc)
1515 : call pawrhoij_alloc(pawrhoij1,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
1516 : & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,&
1517 397 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1518 794 : if (cplex_rhoij/=hdr%pawrhoij(1)%cplex_rhoij.or.qphase_rhoij/=hdr%pawrhoij(1)%qphase) then
1519 : ! Eventually reallocate hdr%pawrhoij
1520 0 : call pawrhoij_free(hdr%pawrhoij)
1521 : call pawrhoij_alloc(hdr%pawrhoij,cplex_rhoij,nspden_rhoij,dtset%nspinor,dtset%nsppol,&
1522 : & dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,&
1523 0 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1524 : end if
1525 : else
1526 1945 : ABI_MALLOC(pawrhoij1,(0))
1527 : end if
1528 :
1529 : ! Initialize 1st-order wavefunctions
1530 2342 : formeig=1; ask_accurate=0; optorth=0
1531 : ! NB: 4 Sept 2013: this was introducing a bug - for ieig2rf ==0 the dim was being set to 1 and passed to dfpt_vtowfk
1532 2342 : dim_eig2rf=0
1533 2342 : if ((dtset%ieig2rf > 0 .and. dtset%ieig2rf/=2) .or. dtset%efmas > 0) then
1534 543 : dim_eig2rf=1
1535 : end if
1536 2342 : mcg1=mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol
1537 2342 : if (one*mpw1*dtset%nspinor*dtset%mband*mk1mem_rbz*dtset%nsppol > huge(1)) then
1538 : write (msg,'(4a, 5(a,i0), 2a)')&
1539 0 : "Default integer is not wide enough to store the size of the GS wavefunction array (WFK1, mcg1).",ch10,&
1540 0 : "Action: increase the number of processors. Consider also OpenMP threads.",ch10,&
1541 0 : "nspinor: ",dtset%nspinor, "mpw1: ",mpw1, "mband: ",dtset%mband, "mk1mem_rbz: ",&
1542 0 : mk1mem_rbz, "nsppol: ",dtset%nsppol,ch10,&
1543 0 : 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS/LAPACK...) compiled in int64 mode'
1544 0 : ABI_ERROR(msg)
1545 : end if
1546 7026 : ABI_MALLOC_OR_DIE(cg1,(2,mcg1), ierr)
1547 : ! space for all 3 ddk wavefunctions if call to orbmag will be needed
1548 : if ( (dtset%orbmag .NE. 0) .AND. &
1549 : & ( (dtset%rfddk .EQ. 1) .OR. (dtset%rfelfd .EQ. 2) ) &
1550 2342 : & .AND. (.NOT. ALLOCATED(cg1_3)) ) then
1551 64 : ABI_MALLOC(cg1_3,(2,mcg1,3))
1552 16 : has_cg1_3(:) = .FALSE.
1553 : end if
1554 2342 : if (.not.kramers_deg) then
1555 11 : mcg1mq=mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol
1556 33 : ABI_MALLOC_OR_DIE(cg1_mq,(2,mcg1mq), ierr)
1557 : end if
1558 :
1559 7026 : ABI_MALLOC(cg1_active,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1560 4684 : ABI_MALLOC(gh1c_set,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1561 4684 : ABI_MALLOC(gh0c1_set,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1562 2342 : if (.not.kramers_deg) then
1563 33 : ABI_MALLOC(cg1_active_mq,(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1564 22 : ABI_MALLOC(gh1c_set_mq,(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1565 22 : ABI_MALLOC(gh0c1_set_mq,(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
1566 : end if
1567 : ! XG090606 This is needed in the present 5.8.2 , for portability for the pathscale machine.
1568 : ! However, it is due to a bug to be corrected by Paul Boulanger. When the bug will be corrected,
1569 : ! this line should be removed.
1570 2342 : if(mk1mem_rbz/=0 .and. dtset%ieig2rf/=0)then
1571 155063906 : cg1_active=zero
1572 155063906 : gh1c_set=zero
1573 155063906 : gh0c1_set=zero
1574 494 : if (.not.kramers_deg) then
1575 0 : cg1_active_mq=zero
1576 0 : gh1c_set_mq=zero
1577 0 : gh0c1_set_mq=zero
1578 : end if
1579 : end if
1580 : ! MG TODO:
1581 : ! eigen1 is not MPI-distributed, it scales as mband^2 and requires a lot of memory for dense k-meshes e.g. metals:
1582 : ! this is what you get with nband 26 and ngkpt 42**3:
1583 : !
1584 : ! <var=eigen1, A@m_dfpt_looppert.F90:1499, addr=0x150d756c0010, size_mb=391.204>
1585 : !
1586 : ! we should try to MPI-distribute it over nkpt_rbs/nsppol or, at the very least, tell users to activate OpenMP.
1587 :
1588 7026 : ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol))
1589 7026 : ABI_MALLOC(resid,(dtset%mband*nkpt_rbz*dtset%nsppol))
1590 2342 : call timab(144,1,tsec)
1591 2342 : if ((file_exists(nctk_ncify(fiwf1i)) .or. file_exists(fiwf1i)) .and. (dtset%get1wf /= 0 .or. dtset%ird1wf /= 0)) then
1592 : call wfk_read_my_kptbands(fiwf1i, distrb_flags, spacecomm, dtset%ecut*(dtset%dilatmx)**2,&
1593 : formeig, istwfk_rbz, kpq_rbz, mcg1, dtset%mband, mband_mem_rbz, mk1mem_rbz, mpw1,&
1594 : dtset%natom, nkpt_rbz, npwar1, dtset%nspinor, dtset%nsppol, dtset%usepaw,&
1595 17 : cg1, eigen=eigen1, ask_accurate_=0)
1596 : else
1597 428786313 : cg1 = zero
1598 12951427 : eigen1 = zero
1599 : end if
1600 :
1601 2342 : call timab(144,2,tsec)
1602 2342 : if(.not.kramers_deg) then
1603 33 : ABI_MALLOC(eigen1_mq,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol))
1604 33 : ABI_MALLOC(resid_mq,(dtset%mband*nkpt_rbz*dtset%nsppol))
1605 : !initialize cg1_mq:
1606 11 : call timab(144,1,tsec)
1607 11 : if ((file_exists(nctk_ncify(fiwf1i_mq)) .or. file_exists(fiwf1i_mq)) .and. (dtset%get1wf > 0 .or. dtset%ird1wf > 0)) then
1608 : call wfk_read_my_kptbands(fiwf1i_mq, distrb_flags, spacecomm, dtset%ecut*(dtset%dilatmx)**2, &
1609 : formeig, istwfk_rbz, kmq_rbz, mcg1mq, dtset%mband, mband_mem_rbz, mk1mem_rbz, mpw1_mq,&
1610 : dtset%natom, nkpt_rbz, npwar1_mq, dtset%nspinor, dtset%nsppol, dtset%usepaw,&
1611 0 : cg1_mq, eigen=eigen1_mq, ask_accurate_=0)
1612 : else
1613 4765835 : cg1_mq = zero
1614 350998 : eigen1_mq = zero
1615 : end if
1616 11 : call timab(144,2,tsec)
1617 : end if
1618 :
1619 : ! Eventually reytrieve 1st-order PAW occupancies from file header
1620 2342 : if (psps%usepaw==1.and.dtfil%ireadwf/=0) then
1621 0 : call pawrhoij_copy(hdr%pawrhoij,pawrhoij1,comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab)
1622 : end if
1623 :
1624 : ! In case of electric field, or 2nd order perturbation: open the ddk (or dE) wf file(s)
1625 : if ((ipert==dtset%natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and. &
1626 : (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and. &
1627 : dtset%berryopt/= 7.and.dtset%berryopt/=14.and. &
1628 : dtset%berryopt/=16.and.dtset%berryopt/=17)) .or. &
1629 9368 : ipert==dtset%natom+10.or.ipert==dtset%natom+11) then
1630 :
1631 370 : if (ipert<dtset%natom+10) then
1632 : ! 1st order or berry phase, one direction
1633 191 : nwffile=1 ! one direction needed
1634 191 : file_index(1)=idir+dtset%natom*3
1635 191 : fnamewff(1)=dtfil%fnamewffddk
1636 179 : else if (ipert==dtset%natom+10) then
1637 : ! 2nd order (k,k)
1638 131 : if (idir<=3) then ! one direction needed
1639 49 : nwffile=1
1640 49 : file_index(1)=idir+dtset%natom*3
1641 49 : fnamewff(1)=dtfil%fnamewffddk
1642 : else ! two directions needed
1643 82 : nwffile=2
1644 82 : file_index(1)=idir1+dtset%natom*3
1645 82 : file_index(2)=idir2+dtset%natom*3
1646 82 : fnamewff(1)=dtfil%fnamewffddk
1647 82 : fnamewff(2)=dtfil%fnamewffddk
1648 : end if
1649 48 : else if (ipert==dtset%natom+11) then
1650 : ! 2nd order (k,E)
1651 48 : nwffile=3 ! dk, dE and dkdk
1652 48 : idir_dkdk = idir
1653 48 : if(idir_dkdk>6) idir_dkdk = idir_dkdk - 3
1654 48 : file_index(1)=idir_dkdk +(dtset%natom+6)*3 ! dkdk
1655 48 : file_index(2)=idir2+(dtset%natom+1)*3 ! defld file (dir2)
1656 48 : file_index(3)=idir1+dtset%natom*3 ! ddk file (dir1)
1657 48 : fnamewff(1)=dtfil%fnamewffdkdk
1658 48 : fnamewff(2)=dtfil%fnamewffdelfd
1659 48 : fnamewff(3)=dtfil%fnamewffddk
1660 48 : if (idir>3) then
1661 32 : nwffile=4
1662 32 : file_index(4)=idir2+dtset%natom*3 ! ddk file (dir2)
1663 32 : fnamewff(4)=dtfil%fnamewffddk
1664 : end if
1665 : end if
1666 950 : do ii=1,nwffile
1667 580 : call appdig(file_index(ii),fnamewff(ii),fiwfddk)
1668 : ! Checking the existence of data file
1669 580 : if (.not. file_exists(fiwfddk)) then
1670 : ! Trick needed to run Abinit test suite in netcdf mode.
1671 580 : if (file_exists(nctk_ncify(fiwfddk))) then
1672 580 : write(msg,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name."
1673 580 : call wrtout(std_out,msg)
1674 1160 : fiwfddk = nctk_ncify(fiwfddk)
1675 : end if
1676 580 : if (.not. file_exists(fiwfddk)) then
1677 0 : ABI_ERROR('Missing file: '//TRIM(fiwfddk))
1678 : end if
1679 : end if
1680 580 : write(msg,'(2a)')'- dfpt_looppert: read the DDK wavefunctions from file: ',trim(fiwfddk)
1681 1740 : call wrtout([std_out, ab_out],msg)
1682 : ! Note that the unit number for these files is 50,51,52 or 53 (dtfil%unddk=50)
1683 2922 : call ddk_f(ii)%open_read(fiwfddk,formeig1,dtset%iomode,dtfil%unddk+(ii-1), spacecomm) !xmpi_comm_self)
1684 : end do
1685 : end if
1686 :
1687 : ! Get first-order local potentials and 1st-order core correction density change
1688 : ! (do NOT include xccc3d1 in vpsp1 : this will be done in dfpt_scfcv because vpsp1
1689 : ! might become spin-polarized)
1690 :
1691 2342 : n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf
1692 7026 : ABI_MALLOC(xccc3d1,(cplex*n3xccc))
1693 7026 : ABI_MALLOC(vpsp1,(cplex*nfftf))
1694 :
1695 : ! PAW: compute Vloc(1) and core(1) together in reciprocal space
1696 : ! --------------------------------------------------------------
1697 2342 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
1698 399 : ndir=1
1699 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,istr,ipert,&
1700 : & mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,&
1701 : & ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
1702 399 : & atmrhor1=xccc3d1,atmvlocr1=vpsp1,optn_in=n3xccc/nfftf,optn2_in=1,vspl=psps%vlspl)
1703 : else
1704 :
1705 : ! Norm-conserving psp: compute Vloc(1) in reciprocal sp. and core(1) in real sp.
1706 : ! ------------------------------------------------------------------------------
1707 :
1708 1943 : if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then
1709 : ! Section for strain perturbation
1710 :
1711 : !To compute Absolute Deformation Potentials together with FxE tensor
1712 : !the reference has to be the same as in the FxE routines
1713 96 : g0term=0; if (dtset%rfstrs_ref==1) g0term=1
1714 :
1715 : call vlocalstr(gmet,gprimd,gsqcut,istr,mgfftf,mpi_enreg,&
1716 : psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,ntypat,ph1df,psps%qgrid_vl,&
1717 96 : ucvol,psps%vlspl,vpsp1,g0term=g0term)
1718 : else
1719 : call dfpt_vlocal(atindx,cplex,gmet,gsqcut,dtset%icutcoul,idir,ipert,mpi_enreg,psps%mqgrid_vl,dtset%natom,&
1720 : nattyp,nfftf,ngfftf,nkpt,ntypat,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,psps%qgrid_vl,&
1721 : dtset%qptn,dtset%rcut,rprimd,ucvol,dtset%vcutgeo,psps%vlspl,vpsp1,xred,&
1722 1847 : zion=dtset%ziontypat)
1723 : !SPr: need vpsp1 for -q as well, but for magnetic field it's zero, to be done later
1724 : end if
1725 :
1726 1943 : if(psps%n1xccc/=0)then
1727 : call dfpt_mkcore(cplex,idir,ipert,dtset%natom,ntypat,ngfftf(1),psps%n1xccc,&
1728 1191 : ngfftf(2),ngfftf(3),dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
1729 : !SPr: same here, need xccc3d1 for -q as well for phonon pert.. to be done later
1730 : end if ! psps%n1xccc/=0
1731 : end if ! usepaw
1732 :
1733 13876454 : eigen1(:)=zero; resid(:)=zero
1734 2342 : if(.not.kramers_deg) then
1735 350987 : eigen1_mq(:)=zero
1736 4427 : resid_mq(:)=zero
1737 : end if
1738 : ! Get starting charge density and Hartree + xc potential
1739 9368 : ABI_MALLOC(rhor1,(cplex*nfftf,nspden))
1740 7026 : ABI_MALLOC(rhog1,(2,nfftf))
1741 :
1742 : ! MR:do we really need to define the partial rhos here?
1743 2342 : if(.not.kramers_deg) then
1744 : ! Case when first order spinors at both +q and -q are not related by symmetry (time and/or space inversion)
1745 33 : ABI_MALLOC(rhor1_pq,(cplex*nfftf,nspden))
1746 22 : ABI_MALLOC(rhog1_pq,(2,nfftf))
1747 33 : ABI_MALLOC(rhor1_mq,(cplex*nfftf,nspden))
1748 22 : ABI_MALLOC(rhog1_mq,(2,nfftf))
1749 : end if
1750 :
1751 : ! can we get this set of gkk matrices from previously calculated rhog1 through a non-scf calculation?
1752 2342 : found_eq_gkk=.false.
1753 2342 : if (dtset%prepgkk == 1 .and. ipert <= dtset%natom) then
1754 : ! NOTE: this does not take into account combinations e.g. x+y -> z
1755 : ! if rhor1 add linearly this could be done...
1756 401 : do icase_eq = 1, icase-1
1757 276 : idir_eq = mod(icase_eq,3)
1758 276 : if (idir_eq==0) idir_eq=3
1759 276 : ipert_eq = ( (icase_eq-idir_eq) / 3 + 1)
1760 :
1761 : ! find sym which links old perturbation to present one
1762 6066 : do isym=1, nsym
1763 : ! check that isym preserves qpt to begin with
1764 : if (symq(4,1,isym) /= 1 .or. &
1765 : symq(1,1,isym) /= 0 .or. &
1766 5815 : symq(2,1,isym) /= 0 .or. &
1767 : symq(3,1,isym) /= 0 ) cycle
1768 :
1769 12298 : eq_symop = dtset%symrel(:,:,isym)
1770 : if (indsym(4,isym,ipert) == ipert_eq .and. &
1771 3784 : & abs(eq_symop(idir,idir_eq)) == 1 .and. &
1772 251 : & sum(abs(eq_symop(:,idir_eq))) == 1) then
1773 25 : found_eq_gkk = .true.
1774 25 : exit
1775 : end if
1776 : end do ! isym
1777 401 : if (found_eq_gkk) exit
1778 : end do ! icase_eq
1779 : end if ! check for prepgkk with symmetric pert
1780 :
1781 2342 : if (found_eq_gkk) then
1782 : write (msg, '(a,l6,i6,2a,3i6,2a,3i6,a)') &
1783 25 : ' found_eq_gkk,isym = ', found_eq_gkk, isym, ch10, &
1784 25 : ' idir, ipert, icase = ', idir, ipert, icase, ch10, &
1785 50 : ' idireq iperteq icaseeq = ', idir_eq, ipert_eq, icase_eq, ch10
1786 25 : call wrtout(std_out,msg)
1787 : !
1788 : ! Make density for present perturbation, which is symmetric of 1 or more previous perturbations:
1789 : ! rotate 1DEN arrays with symrel(isym) to produce rhog1_eq rhor1_eq
1790 : !
1791 25 : if (dtset%use_nonscf_gkk == 1) then
1792 : call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, &
1793 0 : & rhor1_save(:,:,icase_eq), rhog1, rhor1, eq_symop, dtset%tnons(:,isym))
1794 :
1795 0 : rhor1 = rhor1 * eq_symop(idir,idir_eq)
1796 :
1797 : ! TODO: rotate rhoij in PAW case
1798 :
1799 0 : blkflg_save = blkflg
1800 0 : dtset_tmp%iscf = -2
1801 0 : iscf_mod = -2
1802 0 : dtset_tmp%nstep = 1
1803 0 : dtset_tmp%nline = 1
1804 0 : if (abs(dtset_tmp%tolwfr) < 1.e-24) dtset_tmp%tolwfr = 1.e-24
1805 0 : dtset_tmp%toldfe = zero
1806 0 : dtset_tmp%toldff = zero
1807 0 : dtset_tmp%toldmag = zero
1808 0 : dtset_tmp%tolrff = zero
1809 0 : dtset_tmp%tolvrs = zero
1810 0 : write (msg, '(a,i6,a)') ' NOTE: doing GKK calculation for icase ', icase, ' with non-SCF calculation'
1811 0 : call wrtout(std_out,msg)
1812 : !call wrtout(ab_out,msg,'COLL') ! decomment and update output files
1813 :
1814 : else ! do not use non-scf shortcut, but save rotated 1DEN for comparison
1815 : ! saves the rotated rho, for later comparison with the full SCF rhor1: comment lines below for iscf = -2
1816 : call rotate_rho(cplex, timrev_pert, mpi_enreg, nfftf, ngfftf, nspden, &
1817 : & rhor1_save(:,:,icase_eq), rhog1, rhor1_save(:,:,icase), eq_symop, &
1818 25 : & dtset%tnons(:,isym))
1819 43858 : rhor1_save(:,:,icase) = rhor1_save(:,:,icase) * eq_symop(idir,idir_eq)
1820 :
1821 : end if ! force non scf calculation of other gkk, or not
1822 : end if ! found an equiv perturbation for the gkk
1823 :
1824 : if ( (dtfil%ireadwf==0 .and. iscf_mod/=-4 .and. dtset%get1den==0 .and. dtset%ird1den==0) &
1825 2342 : .or. (iscf_mod== -3 .and. ipert/=dtset%natom+11) ) then
1826 : ! NOTE : For ipert==natom+11, we want to read the 1st order density from a previous calculation
1827 98824914 : rhor1(:,:)=zero ; rhog1(:,:)=zero
1828 : ! PAW: rhoij have been set to zero in call to pawrhoij_alloc above
1829 :
1830 : init_rhor1 = ((ipert>=1 .and. ipert<=dtset%natom).or.ipert==dtset%natom+5 &
1831 2162 : & .or.ipert==dtset%natom+6.or.ipert==dtset%natom+7)
1832 : ! This section is needed in order to maintain the old behavior and pass the automatic tests
1833 2162 : if (psps%usepaw == 0) then
1834 2438 : init_rhor1 = init_rhor1 .and. all(psps%nctab(:ntypat)%has_tvale)
1835 : else
1836 : init_rhor1 = .False.
1837 : end if
1838 :
1839 1795 : if (init_rhor1) then
1840 :
1841 213 : if(ipert/=dtset%natom+5) then
1842 :
1843 : ! Initialize rhor1 and rhog1 from the derivative of atomic densities/gaussians.
1844 210 : ndir = 1; optn2 = 3
1845 : if (psps%usepaw==1) then
1846 : ! FIXME: Here there's a bug because has_tvale == 0 or 1 instead of 2
1847 : if (all(pawtab(:ntypat)%has_tvale/=0)) optn2=2
1848 : else if (psps%usepaw==0) then
1849 487 : if (all(psps%nctab(:ntypat)%has_tvale)) optn2=2
1850 : end if
1851 :
1852 : if (optn2 == 3) then
1853 0 : call wrtout(std_out," Initializing rhor1 from atom-centered gaussians")
1854 0 : ABI_MALLOC(gauss,(2,ntypat))
1855 0 : call atom_gauss(ntypat, dtset%densty, psps%ziontypat, psps%znucltypat, gauss)
1856 :
1857 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,&
1858 : mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,&
1859 : ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
1860 0 : atmrhor1=rhor1,optn_in=1,optn2_in=3,gauss=gauss)
1861 :
1862 0 : ABI_FREE(gauss)
1863 : else
1864 210 : call wrtout(std_out," Initializing rhor1 from valence densities taken from pseudopotential files")
1865 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,&
1866 : mgfftf,psps%mqgrid_vl,dtset%natom,ndir,nfftf,ngfftf,ntypat,&
1867 : ph1df,psps%qgrid_vl,dtset%qptn,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
1868 210 : atmrhor1=rhor1,optn_in=1,optn2_in=2)
1869 : end if
1870 :
1871 : else
1872 : ! Magnetic field perturbation
1873 3 : call wrtout(std_out," Initializing rhor1 guess based on the ground state XC magnetic field")
1874 :
1875 3 : call dfpt_init_mag1(ipert,idir,rhor1,rhor,cplex,nfftf,nspden,vxc,kxc,nkxc)
1876 :
1877 : end if
1878 :
1879 213 : call fourdp(cplex,rhog1,rhor1,-1,mpi_enreg,nfftf,1,ngfftf,0)
1880 : end if
1881 :
1882 : else
1883 : ! rhor1 not being forced to 0.0
1884 180 : if(iscf_mod>0) then
1885 : ! cplex=2 gets the complex density, =1 only real part
1886 8 : if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
1887 0 : call appdig(pertcase,dtfil%fildens1in,fiden1i)
1888 : call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor1, &
1889 0 : hdr_den, pawrhoij1, spaceComm, check_hdr=hdr)
1890 0 : etotal = hdr_den%etot; call hdr_den%free()
1891 :
1892 : ! Compute up+down rho1(G) by fft
1893 0 : ABI_MALLOC(work,(cplex*nfftf))
1894 0 : work(:)=rhor1(:,1)
1895 0 : call fourdp(cplex,rhog1,work,-1,mpi_enreg,nfftf,1,ngfftf,0)
1896 0 : ABI_FREE(work)
1897 8 : else if (dtset%get1wf /= 0 .or. dtset%ird1wf /= 0 ) then
1898 8 : if (psps%usepaw==1) then
1899 : ! Be careful: in PAW, rho does not include the 1st-order compensation density (to be added in dfpt_scfcv.F90) !
1900 0 : ABI_MALLOC(rho1wfg,(2,dtset%nfft))
1901 0 : ABI_MALLOC(rho1wfr,(dtset%nfft,nspden))
1902 : call dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,&
1903 : kg,kg1,dtset%mband,mband_mem_rbz,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,&
1904 : dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,&
1905 0 : occ_rbz,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1,tnons1,ucvol,wtk_rbz)
1906 0 : call transgrid(cplex,mpi_enreg,nspden,+1,1,1,dtset%paral_kgb,pawfgr,rho1wfg,rhog1,rho1wfr,rhor1)
1907 0 : ABI_FREE(rho1wfg)
1908 0 : ABI_FREE(rho1wfr)
1909 : else
1910 : call dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon1,istwfk_rbz,&
1911 : kg,kg1,dtset%mband,mband_mem_rbz,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1,nband_rbz,&
1912 : dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,nsym1,&
1913 8 : occ_rbz,phnons1,rhog1,rhor1,rprimd,symaf1,symrl1,tnons1,ucvol,wtk_rbz)
1914 :
1915 8 : if (.not.kramers_deg) then
1916 0 : rhor1_pq(:,:)=rhor1(:,:)
1917 0 : rhog1_pq(:,:)=rhog1(:,:)
1918 : call dfpt_mkrho(cg,cg1_mq,cplex,gprimd,irrzon1,istwfk_rbz,&
1919 : kg,kg1_mq,dtset%mband,mband_mem_rbz,dtset%mgfft,mkmem_rbz,mk1mem_rbz,mpi_enreg,mpw,mpw1_mq,nband_rbz,&
1920 : dtset%nfft,dtset%ngfft,nkpt_rbz,npwarr,npwar1_mq,nspden,dtset%nspinor,dtset%nsppol,nsym1,&
1921 0 : occ_rbz,phnons1,rhog1_mq,rhor1_mq,rprimd,symaf1,symrl1,tnons1,ucvol,wtk_rbz)
1922 :
1923 0 : do ifft=1,nfftf
1924 0 : rhor1(2*ifft-1,:) = half*(rhor1_pq(2*ifft-1,:)+rhor1_mq(2*ifft-1,:))
1925 0 : rhor1(2*ifft ,:) = half*(rhor1_pq(2*ifft ,:)-rhor1_mq(2*ifft ,:))
1926 : end do
1927 0 : call fourdp(cplex,rhog1,rhor1(:,1),-1,mpi_enreg,nfftf,1, ngfftf, 0)
1928 : end if
1929 : end if
1930 : end if
1931 :
1932 172 : else if (.not. found_eq_gkk) then
1933 : ! negative iscf_mod and no symmetric rotation of rhor1
1934 : ! Read rho1(r) from a disk file and broadcast data.
1935 172 : rdwr=1;rdwrpaw=psps%usepaw;if(dtfil%ireadwf/=0) rdwrpaw=0
1936 172 : if (ipert/=dtset%natom+11) then
1937 124 : call appdig(pertcase,dtfil%fildens1in,fiden1i)
1938 : else
1939 : ! For ipert==natom+11, we want to read the 1st order density from a previous calculation
1940 48 : call appdig(idir2+(dtset%natom+1)*3,dtfil%fildens1in,fiden1i)
1941 : end if
1942 : ! call appdig(pertcase,dtfil%fildens1in,fiden1i)
1943 : call read_rhor(fiden1i, cplex, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor1, &
1944 172 : hdr_den, pawrhoij1, spaceComm, check_hdr=hdr)
1945 172 : etotal = hdr_den%etot; call hdr_den%free()
1946 :
1947 : ! Compute up+down rho1(G) by fft
1948 344 : ABI_MALLOC(work,(cplex*nfftf))
1949 2128912 : work(:)=rhor1(:,1)
1950 172 : call fourdp(cplex,rhog1,work,-1,mpi_enreg,nfftf,1,ngfftf,0)
1951 172 : ABI_FREE(work)
1952 : end if ! rhor1 generated or read in from file
1953 :
1954 : end if ! rhor1 set to 0 or read in from file
1955 :
1956 : ! Check whether exiting was required by the user.
1957 : ! If found then do not start minimization steps
1958 2342 : openexit=1 ; if(dtset%chkexit==0) openexit=0
1959 2342 : call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit)
1960 : ! If immediate exit, and wavefunctions were not read, must zero eigenvalues
1961 2342 : if (iexit/=0) eigen1(:)=zero
1962 2342 : if (iexit/=0.and.(.not.kramers_deg)) eigen1_mq(:)=zero
1963 2342 : call pstat_proc%print(_PSTAT_ARGS_)
1964 :
1965 2342 : if (iexit==0) then
1966 :
1967 : ! Main calculation: get 1st-order wavefunctions from Sternheimer equation (SCF cycle)
1968 : ! if ipert==natom+10 or natom+11 : get 2nd-order wavefunctions
1969 2342 : if (kramers_deg) then
1970 : call dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,&
1971 : & dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,&
1972 : & d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
1973 : & ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,elmag1,eloc0,elpsp1,emagpen1,&
1974 : & end0,end1,enl0,enl1,eovl1,epaw1,etotal,evxctau0,evxctau1,evdw,exc1,fermie,gh0c1_set,gh1c_set,hdr,idir,&
1975 : & indkpt1,indsy1,initialized,ipert,irrzon1,istwfk_rbz,&
1976 : & kg,kg1,kpt_rbz,kxc,mband_mem_rbz,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,&
1977 : & mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,&
1978 : & nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,&
1979 : & npwarr,npwar1,nspden,&
1980 : & nsym1,n3xccc,occkq,occ_rbz,&
1981 : & paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,&
1982 : & pertcase,phnons1,ph1d,ph1df,prtbbb,psps,&
1983 : & dtset%qptn,resid,residm,rhog,rhog1,&
1984 : & rhor,rhor1,rprimd,symaf1,symrc1,symrl1,tnons1,&
1985 : & usecprj,useylmgr,useylmgr1,usevxctau,ddk_f,vpsp1,vtrial,vxc,vxctau,&
1986 : & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,dfpt_scfcv_retcode,&
1987 2331 : & kramers_deg)
1988 : else
1989 : call dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,&
1990 : & dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset_tmp,&
1991 : & d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
1992 : & ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,elmag1,eloc0,elpsp1,emagpen1,&
1993 : & end0,end1,enl0,enl1,eovl1,epaw1,etotal,evxctau0,evxctau1,evdw,exc1,fermie,gh0c1_set,gh1c_set,hdr_mq,idir,&
1994 : & indkpt1,indsy1,initialized,ipert,irrzon1,istwfk_rbz,&
1995 : & kg,kg1,kpt_rbz,kxc,mband_mem_rbz,mgfftf,mkmem_rbz,mkqmem_rbz,mk1mem_rbz,&
1996 : & mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,&
1997 : & nattyp,nband_rbz,ncpgr,nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,&
1998 : & npwarr,npwar1,nspden,&
1999 : & nsym1,n3xccc,occkq,occ_rbz,&
2000 : & paw_an_pert,paw_ij_pert,pawang,pawang1,pawfgr,pawfgrtab_pert,pawrad,pawrhoij_pert,pawrhoij1,pawtab,&
2001 : & pertcase,phnons1,ph1d,ph1df,prtbbb,psps,&
2002 : & dtset%qptn,resid,residm,rhog,rhog1,&
2003 : & rhor,rhor1,rprimd,symaf1,symrc1,symrl1,tnons1,&
2004 : & usecprj,useylmgr,useylmgr1,usevxctau,ddk_f,vpsp1,vtrial,vxc,vxctau,&
2005 : & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,dfpt_scfcv_retcode,&
2006 : & kramers_deg,&
2007 : & cg_mq=cg_mq,cg1_mq=cg1_mq,cg1_active_mq=cg1_active_mq,docckde_mq=docckde_mq,eigen_mq=eigen_mq,&
2008 : & eigen1_mq=eigen1_mq,gh0c1_set_mq=gh0c1_set_mq,gh1c_set_mq=gh1c_set_mq,&
2009 : & kg1_mq=kg1_mq,npwar1_mq=npwar1_mq,occk_mq=occk_mq,resid_mq=resid_mq,residm_mq=residm_mq,&
2010 : & rhog1_pq=rhog1_pq,rhog1_mq=rhog1_mq,rhor1_pq=rhor1_pq,rhor1_mq=rhor1_mq,&
2011 11 : & ylm1_mq=ylm1_mq,ylmgr1_mq=ylmgr1_mq)
2012 : end if
2013 :
2014 : ! 2nd-order eigenvalues stuff
2015 2342 : if (dtset%ieig2rf>0) then
2016 494 : if (first_entry) then
2017 84 : nullify(eigen1_pert)
2018 84 : first_entry = .false.
2019 : end if
2020 494 : if (.not.associated(eigen1_pert)) then
2021 420 : ABI_MALLOC(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert))
2022 420 : ABI_MALLOC_OR_DIE(cg1_pert,(2,mpw1*nspinor*mband_mem_rbz*mk1mem_rbz*nsppol*dim_eig2rf,3,mpert),ierr)
2023 420 : ABI_MALLOC(gh0c1_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
2024 252 : ABI_MALLOC(gh1c_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
2025 252 : ABI_MALLOC(kpt_rbz_pert,(3,nkpt_rbz))
2026 336 : ABI_MALLOC(npwarr_pert,(nkpt_rbz,mpert))
2027 252 : ABI_MALLOC(npwar1_pert,(nkpt_rbz,mpert))
2028 252 : ABI_MALLOC(npwtot_pert,(nkpt_rbz,mpert))
2029 8913148 : eigen1_pert(:,:,:) = zero
2030 695929702 : cg1_pert(:,:,:,:) = zero
2031 695929702 : gh0c1_pert(:,:,:,:) = zero
2032 695929702 : gh1c_pert(:,:,:,:) = zero
2033 14136 : npwar1_pert (:,:) = 0
2034 14136 : npwarr_pert (:,:) = 0
2035 6104 : kpt_rbz_pert = kpt_rbz
2036 : end if
2037 494 : clflg(idir,ipert)=1
2038 1997198 : eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:)
2039 494 : if(dtset%ieig2rf==1.or.dtset%ieig2rf==3.or.dtset%ieig2rf==4.or.dtset%ieig2rf==5) then
2040 155063904 : cg1_pert(:,:,idir,ipert)=cg1_active(:,:)
2041 155063904 : gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:)
2042 155063904 : gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:)
2043 : end if
2044 9220 : npwarr_pert(:,ipert)=npwarr(:)
2045 9220 : npwar1_pert(:,ipert)=npwar1(:)
2046 9220 : npwtot_pert(:,ipert)=npwtot(:)
2047 : end if ! eig2rf
2048 :
2049 : ! 2nd-order eigenvalues stuff for EFMAS
2050 2342 : if (dtset%efmas>0) then
2051 51 : if (first_entry) then
2052 17 : nullify(eigen1_pert)
2053 17 : first_entry = .false.
2054 : end if
2055 51 : if (.not.associated(eigen1_pert)) then
2056 85 : ABI_MALLOC(eigen1_pert,(2*dtset%mband**2*nkpt*dtset%nsppol,3,mpert))
2057 85 : ABI_MALLOC(cg1_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
2058 51 : ABI_MALLOC(gh0c1_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
2059 51 : ABI_MALLOC(gh1c_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf,3,mpert))
2060 51 : ABI_MALLOC(kpt_rbz_pert,(3,nkpt_rbz))
2061 68 : ABI_MALLOC(npwarr_pert,(nkpt_rbz,mpert))
2062 286725 : eigen1_pert(:,:,:) = zero
2063 13632366 : cg1_pert(:,:,:,:) = zero
2064 13632366 : gh0c1_pert(:,:,:,:) = zero
2065 13632366 : gh1c_pert(:,:,:,:) = zero
2066 389 : npwarr_pert (:,:) = 0
2067 130 : kpt_rbz_pert = kpt_rbz
2068 51 : ABI_MALLOC(cg0_pert,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
2069 501937 : cg0_pert = cg
2070 : end if
2071 31647 : eigen1_pert(1:2*dtset%mband**2*nkpt_rbz*dtset%nsppol,idir,ipert) = eigen1(:)
2072 1505760 : cg1_pert(:,:,idir,ipert)=cg1_active(:,:)
2073 1505760 : gh0c1_pert(:,:,idir,ipert)=gh0c1_set(:,:)
2074 1505760 : gh1c_pert(:,:,idir,ipert)=gh1c_set(:,:)
2075 123 : npwarr_pert(:,ipert)=npwarr(:)
2076 : end if ! efmas
2077 :
2078 2342 : ABI_FREE(gh1c_set)
2079 2342 : ABI_FREE(gh0c1_set)
2080 2342 : ABI_FREE(cg1_active)
2081 :
2082 : !deallocate big arrays for case without Kramers' degeneracy
2083 2342 : if(.not.kramers_deg) then
2084 11 : ABI_FREE(gh1c_set_mq)
2085 11 : ABI_FREE(gh0c1_set_mq)
2086 11 : ABI_FREE(cg1_active_mq)
2087 : end if
2088 :
2089 : end if ! End of the check of hasty exit
2090 :
2091 2342 : call pstat_proc%print(_PSTAT_ARGS_)
2092 2342 : call timab(146,1,tsec)
2093 :
2094 : ! Print out message at the end of the iterations
2095 189702 : write(msg, '(80a,a,a,a,a)' ) ('=',ii=1,80),ch10,ch10,&
2096 192044 : & ' ----iterations are completed or convergence reached----',ch10
2097 7026 : call wrtout([std_out, ab_out], msg)
2098 :
2099 : ! Print _gkk file for this perturbation
2100 2342 : if (dtset%prtgkk == 1) then
2101 90 : call appdig(3*(ipert-1)+idir,dtfil%fnameabo_gkk,gkkfilnam)
2102 90 : nmatel = dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol
2103 270 : ABI_MALLOC(phasecg, (2, nmatel))
2104 : !call getcgqphase(dtset, timrev, cg, mcg, cgq, mcgq, mpi_enreg, nkpt_rbz, npwarr, npwar1, phasecg)
2105 585690 : phasecg(1,:) = one
2106 585690 : phasecg(2,:) = zero
2107 : ! NB: phasecg not actually used in outgkk for the moment (2013/08/15)
2108 90 : call outgkk(bantot_rbz, nmatel,gkkfilnam,eigen0,eigen1,hdr0,hdr,mpi_enreg,phasecg)
2109 90 : ABI_FREE(phasecg)
2110 :
2111 : ! Reshape eigen1 into gkk for netCDF output
2112 450 : ABI_MALLOC_OR_DIE(gkk,(2*dtset%mband*dtset%nsppol,dtset%nkpt,1,1,dtset%mband), ierr)
2113 1234236 : gkk(:,:,:,:,:) = zero
2114 90 : mband = dtset%mband
2115 90 : band_index = 0
2116 90 : band2tot_index = 0
2117 180 : do isppol=1,dtset%nsppol
2118 6378 : do ikpt =1,nkpt_rbz
2119 66768 : do iband=1,dtset%mband
2120 652368 : do jband=1,dtset%mband
2121 585600 : eig1_r = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index)
2122 585600 : eig1_i = eigen1(2*jband+(iband-1)*2*mband+band2tot_index)
2123 : gkk(2*iband-1+2*band_index,ikpt,1,1,jband) = &
2124 585600 : & gkk(2*iband-1+2*band_index,ikpt,1,1,jband) + eig1_r
2125 : gkk(2*iband+2*band_index,ikpt,1,1,jband) = &
2126 646080 : & gkk(2*iband+2*band_index,ikpt,1,1,jband) + eig1_i
2127 : end do !jband
2128 : end do !iband
2129 6378 : band2tot_index = band2tot_index + 2*mband**2
2130 : end do !ikpt
2131 180 : band_index = band_index + mband
2132 : end do !isppol
2133 :
2134 : ! Initialize ggk_ebands to write in the GKK.nc file
2135 : ! MG FIXME: Here there's a bug because eigen0 is dimensioned with nkpt_rbz i.e. IBZ(q)
2136 : ! but the ebands_t object is constructed with dimensions taken from hdr0 i.e. the IBZ(q=0).
2137 90 : bantot= dtset%mband*dtset%nkpt*dtset%nsppol
2138 : call gkk_ebands%init(bantot,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
2139 : & doccde,eigen0,hdr0%istwfk,hdr0%kptns,&
2140 : & hdr0%nband, hdr0%nkpt,hdr0%npwarr,hdr0%nsppol,hdr0%nspinor,&
2141 : & hdr0%tphysel,hdr0%tsmear,hdr0%occopt,hdr0%occ,hdr0%wtk,&
2142 : & hdr0%cellcharge, hdr0%kptopt, hdr0%kptrlatt_orig, hdr0%nshiftk_orig, hdr0%shiftk_orig, &
2143 90 : & hdr0%kptrlatt, hdr0%nshiftk, hdr0%shiftk)
2144 :
2145 : ! Init a gkk_t object
2146 90 : call gkk_init(gkk,gkk2d,dtset%mband,dtset%nsppol,nkpt_rbz,1,1)
2147 :
2148 : ! Write the netCDF file.
2149 90 : if (me == master) then
2150 90 : fname = strcat(gkkfilnam,".nc")
2151 90 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating GKK file")
2152 90 : NCF_CHECK(crystal%ncwrite(ncid))
2153 90 : NCF_CHECK(gkk_ebands%ncwrite(ncid))
2154 90 : call gkk_ncwrite(gkk2d,dtset%qptn(:),dtset%wtq, ncid)
2155 90 : NCF_CHECK(nf90_close(ncid))
2156 : end if
2157 :
2158 : ! Free memory
2159 90 : ABI_FREE(gkk)
2160 90 : call gkk_free(gkk2d)
2161 90 : call gkk_ebands%free()
2162 : end if
2163 :
2164 2342 : if (dtset%prepgkk == 1 .and. found_eq_gkk) then
2165 25 : if (dtset%use_nonscf_gkk == 1) then
2166 : ! Restore old values of SCF cycle parameters
2167 0 : iscf_mod = iscf_mod_save
2168 0 : dtset_tmp%iscf = iscf_mod_save
2169 0 : dtset_tmp%nstep = nstep_save
2170 0 : dtset_tmp%nline = nline_save
2171 0 : dtset_tmp%tolwfr = tolwfr_save
2172 0 : dtset_tmp%toldfe = toldfe_save
2173 0 : dtset_tmp%toldff = toldff_save
2174 0 : dtset_tmp%toldmag = toldmag_save
2175 0 : dtset_tmp%tolrff = tolrff_save
2176 0 : dtset_tmp%tolvrs = tolvrs_save
2177 0 : blkflg = blkflg_save ! this ensures we do not use the (unconverged) 2DTE from this non scf run
2178 : ! Save density for present perturbation, for future use in symmetric perturbations
2179 0 : rhor1_save(:,:,icase) = rhor1
2180 : else
2181 43858 : write (msg, '(a,3E20.10)') 'norm diff = ', sum(abs(rhor1_save(:,:,icase) - rhor1)), &
2182 131574 : sum(abs(rhor1)), sum(abs(rhor1_save(:,:,icase) - rhor1))/sum(abs(rhor1))
2183 25 : call wrtout(std_out, msg)
2184 : end if
2185 : end if
2186 :
2187 : ! Write wavefunctions file only if convergence was not achieved.
2188 2342 : write_1wfk = .True.
2189 2342 : if (dtset%prtwf == 0) write_1wfk = .False.
2190 2342 : if (dtset%prtwf == -1 .and. dfpt_scfcv_retcode == 0) then
2191 6 : write_1wfk = .False.
2192 6 : call wrtout(ab_out," dfpt_looppert: DFPT cycle converged with prtwf=-1. Will skip output of the 1st-order WFK file.")
2193 : end if
2194 :
2195 : ! store DDK wavefunctions in memory for later call to orbmag
2196 : ! only relevant for DDK pert with orbmag calculation
2197 2342 : if( (dtset%orbmag .NE. 0) .AND. (ipert .EQ. dtset%natom+1) ) then
2198 8181876 : cg1_3(:,:,idir) = cg1(:,:)
2199 48 : has_cg1_3(idir) = .TRUE.
2200 : end if
2201 :
2202 2342 : call outresid(dtset,kpt_rbz,dtset%mband, nband_rbz,nkpt_rbz,dtset%nsppol,resid)
2203 :
2204 2342 : if (write_1wfk) then
2205 : ! Output 1st-order wavefunctions in file
2206 : call wfk_write_my_kptbands(fiwf1o, distrb_flags, spacecomm, formeig, hdr, dtset%iomode, &
2207 : & dtset%mband, mband_mem_rbz, mk1mem_rbz, dtset%mpw, nkpt_rbz, dtset%nspinor, dtset%nsppol, &
2208 1870 : & cg1, kg1, eigen1)
2209 :
2210 1870 : if (.not.kramers_deg) then
2211 : call wfk_write_my_kptbands(fiwf1o_mq, distrb_flags, spacecomm, formeig, hdr_mq, dtset%iomode, &
2212 : & dtset%mband, mband_mem_rbz, mk1mem_rbz, dtset%mpw, nkpt_rbz, dtset%nspinor, dtset%nsppol, &
2213 6 : & cg1_mq, kg1_mq, eigen1_mq)
2214 : end if
2215 :
2216 : ! call outwf(cg1,dtset,psps,eigen1,fiwf1o,hdr,kg1,kpt_rbz,&
2217 : !& dtset%mband,mcg1,mk1mem_rbz,mpi_enreg,mpw1,dtset%natom,nband_rbz,&
2218 : !& nkpt_rbz,npwar1,dtset%nsppol,&
2219 : !& occ_rbz,resid,response,dtfil%unwff2,wvl%wfs,wvl%descr)
2220 : end if
2221 :
2222 : ! Output DDK file in netcdf format.
2223 : ! Can be used by optic instead of the 1WF file that is really huge.
2224 2342 : if (me == master .and. ipert == dtset%natom + 1 .and. dtset%prtevk == 1) then
2225 15 : fname = strcat(dtfil%filnam_ds(4), "_EVK.nc" )
2226 15 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EVK.nc file")
2227 : ! Have to build hdr on k-grid with info about perturbation.
2228 15 : call hdr0%copy(hdr_tmp)
2229 15 : hdr_tmp%kptopt = dtset%kptopt
2230 15 : hdr_tmp%pertcase = pertcase
2231 60 : hdr_tmp%qptn = dtset%qptn(1:3)
2232 15 : NCF_CHECK(hdr_tmp%ncwrite(ncid, 43, nc_define=.True.))
2233 15 : call hdr_tmp%free()
2234 15 : NCF_CHECK(crystal%ncwrite(ncid))
2235 15 : NCF_CHECK(ebands_k%ncwrite(ncid))
2236 : ncerr = nctk_def_arrays(ncid, [ &
2237 : nctkarr_t('h1_matrix_elements', "dp", &
2238 30 : "two, max_number_of_states, max_number_of_states, number_of_kpoints, number_of_spins")], defmode=.True.)
2239 15 : NCF_CHECK(ncerr)
2240 15 : NCF_CHECK(nctk_set_datamode(ncid))
2241 : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "h1_matrix_elements"), eigen1, &
2242 90 : count=[2, dtset%mband, dtset%mband, nkpt_rbz, dtset%nsppol])
2243 15 : NCF_CHECK(ncerr)
2244 15 : NCF_CHECK(nf90_close(ncid))
2245 : end if
2246 :
2247 : ! If the perturbation is d/dk, evaluate the f-sum rule.
2248 2342 : if (ipert==dtset%natom+1 )then
2249 : ! Note : the factor of two is related to the difference
2250 : ! between Taylor expansion and perturbation expansion
2251 : ! Note : this expression should be modified for ecutsm.
2252 : ! Indeed, the present one will NOT tend to 1.0_dp.
2253 436 : ek2=gmet(idir,idir)*(two_pi**2)*2.0_dp*dtset%nelect
2254 436 : fsum=-ek1/ek2
2255 436 : if(dtset%ecutsm<tol6)then
2256 : write(msg, '(a,es20.10,a,a,es20.10)' ) &
2257 352 : ' dfpt_looppert : ek2=',ek2,ch10,&
2258 704 : ' f-sum rule ratio=',fsum
2259 : else
2260 : write(msg, '(a,es20.10,a,a,es20.10,a)' ) &
2261 84 : ' dfpt_looppert : ek2=',ek2,ch10,&
2262 168 : ' f-sum rule ratio=',fsum,' (note : ecutsm/=0)'
2263 : end if
2264 1308 : call wrtout([std_out, ab_out] , msg)
2265 : ! Write the diagonal elements of the dH/dk operator, after averaging over degenerate states
2266 1308 : ABI_MALLOC(eigen1_mean,(dtset%mband*nkpt_rbz*dtset%nsppol))
2267 436 : call eigen_meandege(eigen0,eigen1,eigen1_mean,dtset%mband,nband_rbz,nkpt_rbz,dtset%nsppol,1)
2268 436 : option=4
2269 436 : if (me == master) then
2270 : call prteigrs(eigen1_mean,dtset%enunit,fermie,fermie,dtfil%fnametmp_1wf1_eig,ab_out,iscf_mod,kpt_rbz,dtset%kptopt,&
2271 : & dtset%mband,nband_rbz,dtset%nbdbuf,nkpt_rbz,dtset%nnsclo,dtset%nsppol,occ_rbz,dtset%occopt,&
2272 397 : & option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk_rbz)
2273 : end if
2274 436 : ABI_FREE(eigen1_mean)
2275 : end if
2276 :
2277 : ! Print the energies
2278 2342 : if (dtset%nline/=0 .or. dtset%nstep/=0)then
2279 : call dfpt_prtene(dtset%berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
2280 : & ehart01,ehart1,eii,ek0,ek1,elmag1,eloc0,elpsp1,emagpen1,end0,end1,enl0,enl1,eovl1,epaw1,evxctau0,evxctau1,&
2281 2342 : & evdw,exc1,has_nd,with_vxctau,ab_out,ipert,dtset%natom,psps%usepaw,usevdw)
2282 : end if
2283 :
2284 : ! call orbmag if needed
2285 : if ( (dtset%orbmag .NE. 0) .AND. &
2286 9368 : & ( (dtset%rfddk .EQ. 1) .OR. (dtset%rfelfd .EQ. 2) ) .AND. &
2287 : & (COUNT(has_cg1_3) .EQ. 3) ) then
2288 :
2289 : if ( .NOT. ALLOCATED(vtrial_local)) then
2290 64 : ABI_MALLOC(vtrial_local,(nfftf,dtset%nspden))
2291 : end if
2292 1399217 : vtrial_local = vtrial
2293 16 : if (psps%usepaw .EQ. 1) then
2294 : call orbmag(cg,cg1_3,cprj,crystal,dtfil,dtset,ebands_k,gsqcut,hdr0,kg,mcg,mcg1,&
2295 : & mcprj,mkmem_rbz,mpi_enreg,mpw,nfftf,ngfftf,paw_ij,pawfgr,pawrad,pawtab,psps,&
2296 15 : & usevxctau,vtrial_local,vxctau,ylm,ylmgr)
2297 : else
2298 : call orbmag_ncpp(cg,cg1_3,dtfil,dtset,crystal,ebands_k,kg,gsqcut,hdr0,&
2299 : & mcg,mcg1,mkmem_rbz,mpi_enreg,mpw,nfftf,ngfftf,pawfgr,pawtab,psps,&
2300 1 : & usevxctau,vtrial_local,vxctau)
2301 : end if
2302 16 : ABI_SFREE(vtrial_local)
2303 16 : ABI_SFREE(cg1_3)
2304 16 : has_cg1_3(:) = .FALSE.
2305 : end if ! end call orbmag
2306 :
2307 2342 : if(mpi_enreg%paral_pert==1) then
2308 39 : if (ipert_me < npert_me -1) then
2309 13 : call hdr0%free()
2310 : else
2311 13194 : eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0
2312 : end if
2313 39 : ipert_me = ipert_me +1
2314 : else
2315 2303 : if (icase == ipert_cnt) then
2316 258851 : eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0
2317 : else
2318 1613 : call hdr0%free()
2319 : end if
2320 : end if
2321 :
2322 : if(ipert <= dtset%natom+4) then
2323 : ABI_NVTX_END_RANGE()
2324 : end if
2325 :
2326 : ! Release the temporary arrays (for k, k+q and 1st-order)
2327 2342 : ABI_FREE(cg)
2328 2342 : ABI_FREE(cgq)
2329 2342 : ABI_FREE(cg1)
2330 2342 : ABI_FREE(docckqde)
2331 2342 : if(.not.kramers_deg) then
2332 11 : ABI_FREE(cg_mq)
2333 11 : ABI_FREE(cg1_mq)
2334 11 : ABI_FREE(docckde_mq)
2335 : end if
2336 2342 : ABI_FREE(doccde_rbz)
2337 2342 : ABI_FREE(eigen0)
2338 2342 : ABI_FREE(eigenq)
2339 2342 : ABI_FREE(eigen1)
2340 2342 : ABI_FREE(kpq)
2341 2342 : if(.not.kramers_deg) then
2342 11 : ABI_FREE(eigen_mq)
2343 11 : ABI_FREE(eigen1_mq)
2344 11 : ABI_FREE(kmq)
2345 : end if
2346 2342 : ABI_FREE(indkpt1)
2347 2342 : ABI_FREE(indsy1)
2348 2342 : ABI_FREE(istwfk_rbz)
2349 2342 : ABI_FREE(irrzon1)
2350 2342 : ABI_FREE(kg)
2351 2342 : ABI_FREE(kg1)
2352 2342 : ABI_FREE(kpq_rbz)
2353 2342 : if(.not.kramers_deg) then
2354 11 : ABI_FREE(kg1_mq)
2355 11 : ABI_FREE(kmq_rbz)
2356 : end if
2357 2342 : ABI_FREE(kpt_rbz)
2358 2342 : ABI_FREE(nband_rbz)
2359 2342 : ABI_FREE(npwarr)
2360 2342 : ABI_FREE(npwar1)
2361 2342 : ABI_FREE(npwtot)
2362 2342 : ABI_FREE(npwtot1)
2363 2342 : ABI_FREE(occkq)
2364 2342 : ABI_FREE(occ_rbz)
2365 2342 : ABI_FREE(phnons1)
2366 2342 : ABI_FREE(resid)
2367 2342 : ABI_FREE(rhog1)
2368 2342 : ABI_FREE(rhor1)
2369 2342 : ABI_FREE(symaf1)
2370 2342 : ABI_FREE(symrc1)
2371 2342 : ABI_FREE(symrl1)
2372 2342 : ABI_FREE(tnons1)
2373 2342 : ABI_FREE(wtk_rbz)
2374 2342 : ABI_FREE(xccc3d1)
2375 2342 : ABI_FREE(vpsp1)
2376 2342 : ABI_FREE(ylm)
2377 2342 : ABI_FREE(ylm1)
2378 2342 : ABI_FREE(ylmgr)
2379 2342 : ABI_FREE(ylmgr1)
2380 2342 : if(.not.kramers_deg) then
2381 11 : ABI_FREE(npwar1_mq)
2382 11 : ABI_FREE(npwtot1_mq)
2383 11 : ABI_FREE(occk_mq)
2384 11 : ABI_FREE(resid_mq)
2385 11 : ABI_FREE(rhor1_pq)
2386 11 : ABI_FREE(rhor1_mq)
2387 11 : ABI_FREE(rhog1_pq)
2388 11 : ABI_FREE(rhog1_mq)
2389 11 : ABI_FREE(ylm1_mq)
2390 11 : ABI_FREE(ylmgr1_mq)
2391 : end if
2392 2342 : if (psps%usepaw==1) then
2393 397 : call pawang_free(pawang1)
2394 397 : call pawrhoij_free(pawrhoij1)
2395 397 : if (usecprj==1) then
2396 397 : call pawcprj_free(cprj)
2397 397 : call pawcprj_free(cprjq)
2398 : end if
2399 : end if
2400 3114 : ABI_FREE(pawrhoij1)
2401 161206 : ABI_FREE(cprjq)
2402 161206 : ABI_FREE(cprj)
2403 : if(xmpi_paral==1) then
2404 2342 : ABI_FREE(mpi_enreg%proc_distrb)
2405 2342 : ABI_FREE(mpi_enreg%my_kpttab)
2406 : end if
2407 2342 : call hdr%free()
2408 2342 : if (.not.kramers_deg) then
2409 11 : call hdr_mq%free()
2410 : end if
2411 :
2412 : ! Clean band structure datatypes (should use it more in the future !)
2413 2342 : call ebands_k%free()
2414 2342 : call ebands_kq%free()
2415 2342 : if(.not.kramers_deg) call ebands_kmq%free()
2416 :
2417 : ! %%%% Parallelization over perturbations %%%%%
2418 : ! Redefine output/log files
2419 2342 : call localredirect(mpi_enreg%comm_cell,mpi_enreg%comm_world,npert_io,mpi_enreg%paral_pert,0)
2420 :
2421 2342 : ABI_FREE(bz2ibz_smap)
2422 2342 : ABI_FREE(distrb_flags)
2423 :
2424 2342 : call timab(146,2,tsec)
2425 10084 : if(iexit/=0) exit
2426 : end do ! End loop on perturbations
2427 :
2428 : ABI_NVTX_END_RANGE()
2429 716 : ABI_FREE(zeff)
2430 :
2431 : !%%%% Parallelization over perturbations %%%%%
2432 : !*Restore default communicators
2433 716 : call unset_pert_comm(mpi_enreg)
2434 : ! Gather output/log files
2435 2148 : ABI_MALLOC(dyn,(npert_io))
2436 786 : if (npert_io>0) dyn=1
2437 716 : call localrdfile(mpi_enreg%comm_pert,mpi_enreg%comm_world,.true.,npert_io,mpi_enreg%paral_pert,0,dyn)
2438 716 : ABI_FREE(dyn)
2439 :
2440 : ! Restore PAW on-site data
2441 716 : if (paral_pert_inplace) then
2442 : call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,paw_an,paw_ij,pawfgrtab,pawrhoij)
2443 : else
2444 : call unset_pert_paw(dtset,mpi_enreg,my_natom,old_atmtab,old_comm_atom,&
2445 : & paw_an,paw_ij,pawfgrtab,pawrhoij,&
2446 : & paw_an_out=paw_an_pert,paw_ij_out=paw_ij_pert,&
2447 : & pawfgrtab_out=pawfgrtab_pert,pawrhoij_out=pawrhoij_pert)
2448 : end if
2449 :
2450 : !#################################################################################
2451 : !Calculate the second-order eigenvalues for a wavevector Q
2452 :
2453 716 : call timab(147,1,tsec)
2454 716 : smdelta = dtset%smdelta
2455 716 : bdeigrf = dtset%bdeigrf
2456 716 : if(dtset%bdeigrf == -1) bdeigrf = dtset%mband
2457 :
2458 716 : if(dtset%ieig2rf > 0) then
2459 :
2460 84 : if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
2461 : (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then
2462 3 : call wrtout(std_out,'Reading the dense grid WF file')
2463 : ! We get the Abinit header of the file hdr_fine as output
2464 : ! We get eigenq_fine(mband,hdr_fine%nkpt,hdr_fine%nsppol) as output
2465 3 : fname = dtfil%fnameabi_wfkfine
2466 3 : if (dtset%iomode == IO_MODE_ETSF) fname = nctk_ncify(dtfil%fnameabi_wfkfine)
2467 :
2468 3 : call wfk_read_eigenvalues(fname,eigenq_fine,hdr_fine,mpi_enreg%comm_world)
2469 3 : ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband")
2470 : end if
2471 : ! DBSP ==> Has been changed to be able to make Bandstructure calculation
2472 : ! if(dtset%kptopt==3 .or. dtset%kptopt==0)then
2473 84 : if(dtset%kptopt==3 .or. dtset%kptopt==0 .or. dtset%kptopt < -4 .or. dtset%nsym==1) then
2474 : !END
2475 84 : if (dtset%nsym > 1) then ! .and. dtset%efmas==0) then
2476 0 : ABI_ERROR("Symmetries are not implemented for temperature dependence calculations")
2477 : end if
2478 84 : write(std_out,*) 'Entering: eig2stern'
2479 84 : if(smdelta>0)then
2480 82 : if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
2481 : (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then
2482 : call eig2stern(dtfil,occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
2483 : & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
2484 : & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
2485 3 : & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,xred,pawtab,psps,eigbrd,eigenq_fine,hdr_fine,hdr0)
2486 : else
2487 : call eig2stern(dtfil,occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
2488 : & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
2489 : & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
2490 79 : & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,xred,pawtab,psps,eigbrd)
2491 : end if
2492 : else
2493 2 : if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.&
2494 : & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then
2495 : call eig2stern(dtfil,occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
2496 : & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
2497 : & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
2498 0 : & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,xred,pawtab,psps,eigbrd,eigenq_fine,hdr_fine,hdr0)
2499 : else
2500 : call eig2stern(dtfil,occ_pert,bdeigrf,clflg,cg1_pert,dim_eig2nkq,dim_eig2rf,eigen0_pert,eigenq_pert,&
2501 : & eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,gh0c1_pert,gh1c_pert,&
2502 : & dtset%ieig2rf,istwfk_pert,dtset%mband,mk1mem_rbz,mpert,dtset%natom,mpi_enreg,mpw1,nkpt_rbz,&
2503 2 : & npwar1_pert,dtset%nspinor,dtset%nsppol,smdelta,dtset,xred,pawtab,psps)
2504 : end if
2505 : end if
2506 84 : call wrtout(std_out, 'Leaving: eig2stern')
2507 : !
2508 : else
2509 : write(msg,'(3a)')&
2510 0 : 'K point grids must be the same for every perturbation: eig2stern not called',ch10,&
2511 0 : 'Action: Put kptopt=3 '
2512 0 : ABI_WARNING(msg)
2513 : end if !kptopt
2514 84 : ABI_FREE(gh1c_pert)
2515 84 : ABI_FREE(gh0c1_pert)
2516 84 : ABI_FREE(cg1_pert)
2517 84 : ABI_FREE(kpt_rbz_pert)
2518 84 : ABI_FREE(istwfk_pert)
2519 84 : ABI_FREE(npwarr_pert)
2520 84 : ABI_FREE(npwar1_pert)
2521 84 : ABI_FREE(npwtot_pert)
2522 84 : ABI_FREE(occ_pert)
2523 : end if !if dtset%ieig2rf
2524 :
2525 : ! Calculation of effective masses.
2526 716 : if(dtset%efmas == 1) then
2527 : call efmas_main(cg0_pert,cg1_pert,dim_eig2rf,dtset,efmasdeg,efmasval,eigen0_pert,&
2528 17 : & eigen1_pert,gh0c1_pert,gh1c_pert,istwfk_pert,mpert,mpi_enreg,nkpt_rbz,npwarr_pert,rprimd)
2529 :
2530 17 : ABI_FREE(gh1c_pert)
2531 17 : ABI_FREE(gh0c1_pert)
2532 17 : ABI_FREE(cg1_pert)
2533 17 : ABI_FREE(istwfk_pert)
2534 17 : ABI_FREE(npwarr_pert)
2535 17 : ABI_FREE(cg0_pert)
2536 :
2537 17 : if (dtset%prtefmas == 1 .and. me == master) then
2538 17 : fname = strcat(dtfil%filnam_ds(4), "_EFMAS.nc")
2539 17 : NCF_CHECK_MSG(nctk_open_create(ncid, fname, xmpi_comm_self), "Creating EFMAS file")
2540 17 : NCF_CHECK(crystal%ncwrite(ncid))
2541 : !NCF_CHECK(ebands_k%ncwrite(ncid)) ! At this stage, ebands_k is not available
2542 17 : call print_efmas(efmasdeg, efmasval, kpt_rbz_pert, ncid)
2543 17 : NCF_CHECK(nf90_close(ncid))
2544 : endif
2545 :
2546 17 : call efmas_analysis(dtset,efmasdeg,efmasval,kpt_rbz_pert,mpi_enreg,nkpt_rbz,rprimd)
2547 17 : ABI_FREE(kpt_rbz_pert)
2548 : end if
2549 :
2550 : !Free memory.
2551 716 : if(dtset%ieig2rf /= 3 .and. dtset%ieig2rf /= 4 .and. dtset%ieig2rf /= 5) call hdr0%free()
2552 716 : ABI_FREE(eigen0_copy)
2553 716 : call crystal%free()
2554 :
2555 716 : call timab(147,2,tsec)
2556 : !######################################################################################
2557 :
2558 : !Get ddk file information, for later use in dfpt_dyout
2559 716 : ddkfil(:)=0
2560 2864 : do idir=1,3
2561 2148 : file_index(1)=idir+dtset%natom*3
2562 2148 : call appdig(file_index(1),dtfil%fnamewffddk,fiwfddk)
2563 : ! Check that ddk file exists
2564 2148 : t_exist = file_exists(fiwfddk)
2565 2148 : if (.not. t_exist) then
2566 : ! Trick needed to run Abinit test suite in netcdf mode.
2567 2148 : t_exist = file_exists(nctk_ncify(fiwfddk))
2568 2148 : if (t_exist) then
2569 413 : write(msg,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name."
2570 413 : call wrtout(std_out, msg)
2571 413 : fiwfddk = nctk_ncify(fiwfddk)
2572 : end if
2573 : end if
2574 :
2575 : ! If the file exists set ddkfil to a non-zero value
2576 1129 : if (t_exist) ddkfil(idir)=20+idir
2577 : end do
2578 :
2579 716 : ABI_FREE(ph1d)
2580 716 : ABI_FREE(ph1df)
2581 716 : ABI_FREE(pert_calc)
2582 716 : if (psps%usepaw==1) then
2583 112 : ABI_FREE(dimcprj_srt)
2584 : end if
2585 :
2586 : !destroy dtset_tmp
2587 716 : if (dtset%prepgkk /= 0) then ! .and. dtset%use_nonscf_gkk == 1) then !Later uncomment this - in scf case rhor1_save is used below only for testing
2588 44 : ABI_FREE(rhor1_save)
2589 44 : ABI_FREE(blkflg_save)
2590 44 : call dtset_tmp%free()
2591 44 : ABI_FREE(dtset_tmp)
2592 : end if
2593 :
2594 : ! In paral_pert-case some array's have to be reconstructed
2595 716 : if(mpi_enreg%paral_pert==1) then
2596 130 : ABI_MALLOC(buffer1,(2,3,mpert,3,mpert*(2+psps%usepaw)))
2597 64142 : buffer1(:,:,:,:,1:mpert)=d2lo(:,:,:,:,:)
2598 64142 : buffer1(:,:,:,:,1+mpert:2*mpert)=d2nl(:,:,:,:,:)
2599 26 : if (psps%usepaw==1) then
2600 49340 : buffer1(:,:,:,:,1+2*mpert:3*mpert)=d2ovl(:,:,:,:,:)
2601 : end if
2602 26 : call xmpi_sum(buffer1,mpi_enreg%comm_pert,ierr)
2603 26 : call xmpi_sum(blkflg,mpi_enreg%comm_pert,ierr)
2604 26 : if(dtset%prtbbb==1) then
2605 0 : call xmpi_sum(d2bbb,mpi_enreg%comm_pert,ierr)
2606 : end if
2607 64142 : d2lo(:,:,:,:,:)=buffer1(:,:,:,:,1:mpert)
2608 64142 : d2nl(:,:,:,:,:)=buffer1(:,:,:,:,1+mpert:2*mpert)
2609 26 : if (psps%usepaw==1) then
2610 49340 : d2ovl(:,:,:,:,:)=buffer1(:,:,:,:,1+2*mpert:3*mpert)
2611 : end if
2612 52 : ABI_FREE(buffer1)
2613 : end if
2614 :
2615 716 : if (associated(old_atmtab)) then
2616 36 : ABI_FREE(old_atmtab)
2617 : nullify(old_atmtab)
2618 : end if
2619 :
2620 716 : call timab(141,2,tsec)
2621 :
2622 : DBG_EXIT("COLL")
2623 :
2624 2864 : end subroutine dfpt_looppert
2625 : !!***
2626 :
2627 : !!****f* ABINIT/getcgqphase
2628 : !! NAME
2629 : !! getcgqphase
2630 : !!
2631 : !! FUNCTION
2632 : !! extract phases from wave functions, to cancel contributions to gkk matrix elements
2633 : !!
2634 : !! INPUTS
2635 : !! dtset <type(dataset_type)>=all input variables for this dataset
2636 : !! timrev = flag for use of time reversal symmetry
2637 : !! cg = input wavefunctions
2638 : !! mcg = dimension of cg = nspinor*mband*mpw*mkmem
2639 : !! cgq = input wavefunctions at k+q
2640 : !! mcgq = dimension of cgq = nspinor*mband*mpw*mkmem
2641 : !! mpi_enreg = datastructure for mpi communication
2642 : !! nkpt_rbz = number of k-points in reduced zone for present q point
2643 : !! npwarr = array of numbers of plane waves for each k-point
2644 : !! npwar1 = array of numbers of plane waves for each k+q point
2645 : !!
2646 : !! OUTPUT
2647 : !! phasecg = phase of different wavefunction products <k,n | k+q,n'>
2648 : !!
2649 : !! SOURCE
2650 :
2651 : subroutine getcgqphase(dtset, timrev, cg, mcg, cgq, mcgq, mpi_enreg, nkpt_rbz, npwarr, npwar1, phasecg)
2652 :
2653 : !Arguments -------------------------------
2654 : ! scalars
2655 : integer, intent(in) :: mcg, mcgq, timrev
2656 : integer, intent(in) :: nkpt_rbz
2657 : type(dataset_type), intent(in) :: dtset
2658 : ! arrays
2659 : integer, intent(in) :: npwarr(nkpt_rbz)
2660 : integer, intent(in) :: npwar1(nkpt_rbz)
2661 : real(dp), intent(in) :: cg(2,mcg)
2662 : real(dp), intent(in) :: cgq(2,mcgq)
2663 : type(MPI_type), intent(in) :: mpi_enreg
2664 : real(dp),intent(out) :: phasecg(2, dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)
2665 :
2666 : !Local variables -------------------------
2667 : !scalars
2668 : integer :: icg, icgq, isppol, ikpt, ipw
2669 : integer :: istate, iband1, iband2, npw_k, npw_q
2670 : integer :: me, ierr, master, spaceComm, nprocs, usepaw
2671 : integer :: ddkflag, itrs, job, maxbd, mcg1_k, minbd, shiftbd
2672 : real(dp) :: normsmat
2673 : !arrays
2674 : real(dp) :: cg1_dummy(1,1), smat_inv_dummy(1,1,1), smat_k_paw_dummy(1,1,1), dtm_k_dummy(2)
2675 : integer, allocatable :: sflag_k(:), pwind_k(:)
2676 : real(dp), allocatable :: smat_k(:,:,:), pwnsfac_k(:,:)
2677 : logical, allocatable :: my_kpt(:,:)
2678 : !character(len=500) :: msg
2679 : ! *********************************************************************
2680 :
2681 : ABI_MALLOC(smat_k,(2,dtset%mband,dtset%mband))
2682 : ABI_MALLOC(sflag_k,(dtset%mband))
2683 :
2684 : !dummy use of timrev so abirules stops complaining.
2685 : icg = timrev
2686 :
2687 : !!MPI data for future use
2688 : spaceComm=mpi_enreg%comm_cell
2689 : nprocs=xmpi_comm_size(spaceComm)
2690 : master=0
2691 : me=mpi_enreg%me_kpt
2692 :
2693 : ABI_MALLOC(my_kpt, (nkpt_rbz, dtset%nsppol))
2694 : my_kpt = .true.
2695 : if (mpi_enreg%nproc_spkpt > 1) then
2696 : do isppol = 1, dtset%nsppol
2697 : do ikpt = 1, nkpt_rbz
2698 : my_kpt(ikpt, isppol) = .not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,&
2699 : & dtset%nband(ikpt),isppol,me))
2700 : end do
2701 : end do
2702 : end if
2703 :
2704 :
2705 : !make trivial association of G vectors: we just want <psi_k| psi_k+q>
2706 : !TODO: check this is correct wrt arrangement of kg vectors for k+q
2707 : !looks ok : usually made in initberry, from the translations associated
2708 : !to the symops, scalar product with the G vectors. The symop is the one
2709 : !used to go from the irreducible k to the full zone k. In present context
2710 : !we should be using only the reduced zone, and anyhow have the same k-grid
2711 : !for the gkk matrix elements and for the cg here...
2712 : ABI_MALLOC(pwind_k,(dtset%mpw))
2713 : ABI_MALLOC(pwnsfac_k,(4,dtset%mpw))
2714 : do ipw = 1, dtset%mpw
2715 : pwind_k(ipw) = ipw
2716 : pwnsfac_k(1,ipw) = one
2717 : pwnsfac_k(2,ipw) = zero
2718 : pwnsfac_k(3,ipw) = one
2719 : pwnsfac_k(4,ipw) = zero
2720 : end do
2721 :
2722 : !flags for call to smatrix
2723 : usepaw = 0 ! for now
2724 : ddkflag = 0
2725 : itrs = 0
2726 : job = 0
2727 : maxbd = 1
2728 : mcg1_k = 1
2729 : minbd = 1
2730 : shiftbd = 1
2731 :
2732 : !from overlap matrix for each wavefunction, extract phase
2733 : icg = 0
2734 : icgq = 0
2735 : istate = 0
2736 :
2737 : phasecg = zero
2738 : do isppol = 1, dtset%nsppol
2739 : do ikpt = 1, nkpt_rbz
2740 : ! each proc only has certain k
2741 : if (.not. my_kpt(ikpt, isppol)) then
2742 : istate = istate + dtset%nband(ikpt)*dtset%nband(ikpt)
2743 : cycle
2744 : end if
2745 :
2746 : npw_k = npwarr(ikpt)
2747 : npw_q= npwar1(ikpt)
2748 :
2749 : ! TODO: question: are the k-points in the ibz correctly ordered in cg and cgq? if not the icg below have to be adapted.
2750 : sflag_k = 0 ! make sure all elements are calculated
2751 : smat_k = zero
2752 :
2753 : call smatrix(cg, cgq, cg1_dummy, ddkflag, dtm_k_dummy, icg, icgq,&
2754 : & itrs, job, maxbd, mcg, mcgq, mcg1_k, minbd,dtset%mpw, dtset%mband, dtset%mband,&
2755 : & npw_k, npw_q, dtset%nspinor, pwind_k, pwnsfac_k, sflag_k, shiftbd,&
2756 : & smat_inv_dummy, smat_k, smat_k_paw_dummy, usepaw)
2757 :
2758 : icg = icg + npw_k*dtset%nspinor*dtset%nband(ikpt)
2759 : icgq = icgq + npw_q*dtset%nspinor*dtset%nband(ikpt)
2760 :
2761 : do iband1 = 1, dtset%nband(ikpt)
2762 : do iband2 = 1, dtset%nband(ikpt)
2763 : istate = istate + 1
2764 : ! normalise the overlap matrix element to get just the phase difference phi_k - phi_k+q
2765 : normsmat = sqrt(smat_k(1,iband2, iband1)**2 &
2766 : & + smat_k(2,iband2, iband1)**2)
2767 : if (normsmat > tol12) then
2768 : phasecg(:,istate) = smat_k(:,iband2, iband1) / normsmat
2769 : ! NOTE: 21/9/2011 these appear to be always 1, i, or -i, to within 1.e-5 at worst!
2770 : end if
2771 : end do
2772 : end do
2773 : end do
2774 : end do
2775 :
2776 : !eventually do an mpi allreduce over the k-points for phasecg
2777 : if (nprocs>1) then
2778 : call xmpi_barrier(spaceComm)
2779 : call xmpi_sum_master(phasecg,master,spaceComm,ierr)
2780 : call xmpi_barrier(spaceComm)
2781 : if (1==1) then
2782 : call wrtout(std_out, 'In getcgqphase - contributions to phasecg collected')
2783 : end if
2784 : end if
2785 :
2786 : ABI_FREE(sflag_k)
2787 : ABI_FREE(smat_k)
2788 : ABI_FREE(pwind_k)
2789 : ABI_FREE(pwnsfac_k)
2790 : ABI_FREE(my_kpt)
2791 :
2792 : end subroutine getcgqphase
2793 : !!***
2794 :
2795 : !!****f* ABINIT/dfpt_prtene
2796 : !!
2797 : !! NAME
2798 : !! dfpt_prtene
2799 : !!
2800 : !! FUNCTION
2801 : !! Print components of second derivative of total energy in nice format
2802 : !!
2803 : !! INPUTS
2804 : !! eberry=energy associated with Berry phase
2805 : !! edocc=correction to 2nd-order total energy coming from changes of occupation
2806 : !! eeig0=0th-order eigenenergies part of 2nd-order total energy
2807 : !! eew=Ewald part of 2nd-order total energy
2808 : !! efrhar=hartree frozen-wavefunction part of 2nd-order tot. en.
2809 : !! efrkin=kinetic frozen-wavefunction part of 2nd-order tot. en.
2810 : !! efrloc=local psp. frozen-wavefunction part of 2nd-order tot. en.
2811 : !! efrnl=nonlocal psp. frozen-wavefunction part of 2nd-order tot. en
2812 : !! efrx1=xc core corr.(1) frozen-wavefunction part of 2nd-order tot. en
2813 : !! efrx2=xc core corr.(2) frozen-wavefunction part of 2nd-order tot. en
2814 : !! ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy
2815 : !! for strain perturbation only (zero otherwise, and not used)
2816 : !! ehart1=1st-order Hartree part of 2nd-order total energy
2817 : !! eii=pseudopotential core part of 2nd-order total energy
2818 : !! ek0=0th-order kinetic energy part of 2nd-order total energy.
2819 : !! ek1=1st-order kinetic energy part of 2nd-order total energy.
2820 : !! elmag1=1st-order hspinfield part of 2nd-order total energy.
2821 : !! eloc0=0th-order local (psp+vxc+Hartree) part of 2nd-order total energy
2822 : !! elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
2823 : !! emagpen1=1st-order magnetic penalty part of 2nd-order total energy.
2824 : !! end0=0th-order nuclear dipole part of 2nd-order total energy.
2825 : !! end1=1st-order nuclear dipole part of 2nd-order total energy.
2826 : !! enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
2827 : !! enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
2828 : !! eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy
2829 : !! PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]]
2830 : !! epaw1=1st-order PAW on-site part of 2nd-order total energy.
2831 : !! evxctau0=0th-order energy due to vxctau
2832 : !! evxctau1=1st order energy due to vxctau
2833 : !! evdw=DFT-D semi-empirical part of 2nd-order total energy
2834 : !! exc1=1st-order exchange-correlation part of 2nd-order total energy
2835 : !! has_nd=logical, whether nuclear dipole energies are present
2836 : !! with_vxctau=logical, whether mgga vxctau energies are present
2837 : !! iout=unit number to which output is written
2838 : !! ipert=type of the perturbation
2839 : !! natom=number of atoms in unit cell
2840 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
2841 : !! usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use
2842 : !!
2843 : !! OUTPUT
2844 : !! (only writing)
2845 : !!
2846 : !! NOTES
2847 : !! all energies in Hartree
2848 : !!
2849 : !! SOURCE
2850 :
2851 2342 : subroutine dfpt_prtene(berryopt,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
2852 : & ehart01,ehart1,eii,ek0,ek1,elmag1,eloc0,elpsp1,emagpen1,end0,end1,enl0,enl1,eovl1,epaw1,evxctau0,evxctau1,&
2853 : & evdw,exc1,has_nd,with_vxctau,iout,ipert,natom,usepaw,usevdw)
2854 :
2855 : !Arguments -------------------------------
2856 : !scalars
2857 : integer,intent(in) :: berryopt,iout,ipert,natom,usepaw,usevdw
2858 : real(dp),intent(in) :: eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1
2859 : real(dp),intent(in) :: efrx2,ehart01,ehart1,eii,ek0,ek1,elmag1,eloc0,elpsp1,emagpen1,end0,end1,enl0,enl1
2860 : real(dp),intent(in) :: eovl1,epaw1,evxctau0,evxctau1,evdw,exc1
2861 : logical,intent(in) :: has_nd,with_vxctau
2862 :
2863 : !Local variables -------------------------
2864 : !scalars
2865 : integer :: nn
2866 : logical :: berry_activated
2867 : real(dp) :: enl1_effective,erelax,etotal
2868 : character(len=10) :: numb
2869 : character(len=10),parameter :: numbstr(20) = &
2870 : & (/'One ','Two ','Three ','Four ','Five ', &
2871 : & 'Six ','Seven ','Eight ','Nine ','Ten ', &
2872 : & 'Eleven ','Twelve ','Thirteen ','Fourteen ','Fifteen ', &
2873 : & 'Sixteen ','Seventeen ','Eighteen ','Nineteen ','Twenty '/)
2874 : character(len=500) :: msg
2875 : ! *********************************************************************
2876 :
2877 : !Count and print the number of components of 2nd-order energy
2878 : !MT feb 2015: this number is wrong! Should change it but
2879 : ! need to change a lot of ref. files
2880 : berry_activated=(berryopt== 4.or.berryopt== 6.or.berryopt== 7.or. &
2881 2342 : & berryopt==14.or.berryopt==16.or.berryopt==17)
2882 2342 : if (ipert==natom+1) nn=8
2883 2342 : if (ipert==natom+5.or.ipert==natom+6) nn=7
2884 2342 : if (ipert==natom+7) nn=8
2885 2342 : if (ipert==natom+2) nn=7
2886 2342 : if (ipert>=1.and.ipert<=natom) nn=13
2887 2342 : if (ipert==natom+3.or.ipert==natom+4) nn=17
2888 2342 : if (ipert==natom+2.and.berry_activated) nn=nn+1
2889 2342 : if (ipert==natom+10.or.ipert==natom+11) nn=1 ! means nothing,
2890 2342 : if (ipert>natom+11.and.ipert<=2*natom+11) nn=7
2891 2342 : if (emagpen1>tol8) nn=nn+1
2892 : ! because we do not compute derivatives of the energy in this case
2893 2342 : if (usepaw==1) nn=nn+1
2894 2342 : if (usevdw==1) nn=nn+1
2895 2342 : write(msg, '(4a)' ) ch10,&
2896 4684 : & ' ',trim(numbstr(nn)),' components of 2nd-order total energy (hartree) are '
2897 2342 : call wrtout(iout,msg)
2898 2342 : call wrtout(std_out,msg)
2899 :
2900 2342 : numb='1,2,3'
2901 : write(msg, '(3a)' )&
2902 2342 : & ' ',trim(numb),': 0th-order hamiltonian combined with 1st-order wavefunctions'
2903 2342 : call wrtout(iout,msg)
2904 2342 : call wrtout(std_out,msg)
2905 : write(msg, '(a,es17.8,a,es17.8,a,es17.8)' )&
2906 2342 : & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0
2907 2342 : if (has_nd) then
2908 : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' )&
2909 48 : & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0,' nclr dpl0=',end0
2910 : end if
2911 2342 : if (with_vxctau) then
2912 : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' )&
2913 6 : & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0,' evxctau0=',evxctau0
2914 : end if
2915 2342 : call wrtout(iout,msg)
2916 2342 : call wrtout(std_out,msg)
2917 :
2918 2342 : numb='4,5,6';if( ipert==natom+3.or.ipert==natom+4) numb='4,5,6,7'
2919 : write(msg, '(3a)' )&
2920 2342 : & ' ',trim(numb),': 1st-order hamiltonian combined with 1st and 0th-order wfs'
2921 2342 : call wrtout(iout,msg)
2922 2342 : call wrtout(std_out,msg)
2923 2342 : if(ipert/=natom+1.and.ipert/=natom+2.and.ipert/=natom+5.and.ipert/=natom+6.and.ipert<natom+11)then
2924 : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,a)' ) &
2925 1652 : & ' loc psp =',elpsp1,' Hartree=',ehart1,' xc=',exc1,ch10,&
2926 3304 : & ' note that "loc psp" includes a xc core correction that could be resolved'
2927 690 : else if(ipert==natom+1) then
2928 : write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) &
2929 436 : & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1
2930 436 : if (has_nd .and. (.not.with_vxctau) ) then
2931 : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' ) &
2932 42 : & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1,' nclr dpl1=',end1
2933 : end if
2934 436 : if (with_vxctau .and. (.not.has_nd) ) then
2935 : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' ) &
2936 0 : & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1,' evxctau1=',evxctau1
2937 : end if
2938 436 : if (with_vxctau .and. has_nd ) then
2939 : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8,a,es17.8)' ) &
2940 6 : & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1,' evxctau1=',evxctau1,' nclr dpl1=',end1
2941 : end if
2942 254 : else if(ipert==natom+2) then
2943 : write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) &
2944 194 : & ' dotwf=',enl1, ' Hartree=',ehart1,' xc=',exc1
2945 60 : else if(ipert==natom+5.or.ipert==natom+6.or.(ipert>natom+1.and.ipert<=2*natom+11)) then
2946 : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,a)' ) &
2947 60 : & ' Hspinfield=',elmag1,' Hartree=',ehart1,' xc=',exc1,ch10,&
2948 120 : & ' note that "loc psp" includes a xc core correction that could be resolved'
2949 : end if
2950 2342 : if(ipert==natom+3 .or. ipert==natom+4) then
2951 : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,a,es17.8)' ) &
2952 168 : & ' loc psp =',elpsp1,' Hartree=',ehart1,' xc=',exc1,ch10,&
2953 336 : & ' kin1=',ek1
2954 : end if
2955 2342 : call wrtout(iout,msg)
2956 2342 : call wrtout(std_out,msg)
2957 :
2958 2342 : enl1_effective=enl1;if (ipert==natom+2) enl1_effective=zero
2959 2342 : numb='7,8,9';if( ipert==natom+3.or.ipert==natom+4) numb='8,9,10'
2960 : write(msg, '(5a,es17.8,a,es17.8,a,es17.8)' )&
2961 2342 : & ' ',trim(numb),': eventually, occupation + non-local contributions',ch10,&
2962 4684 : & ' edocc=',edocc,' enl0=',enl0,' enl1=',enl1_effective
2963 2342 : call wrtout(iout,msg)
2964 2342 : call wrtout(std_out,msg)
2965 :
2966 2342 : if (abs(emagpen1)>tol8) then
2967 : write(msg,'(a,es17.8)') &
2968 6 : & ' Magnetic penalty contribution=', emagpen1
2969 6 : call wrtout(iout,msg)
2970 6 : call wrtout(std_out,msg)
2971 : end if
2972 :
2973 2342 : if (usepaw==1) then
2974 397 : numb='10';if( ipert==natom+3.or.ipert==natom+4) numb='11'
2975 : write(msg, '(3a,es17.8)' )&
2976 397 : & ' ',trim(numb),': eventually, PAW "on-site" Hxc contribution: epaw1=',epaw1
2977 397 : call wrtout(iout,msg)
2978 397 : call wrtout(std_out,msg)
2979 : end if
2980 :
2981 2342 : if(ipert/=natom+10 .and.ipert/=natom+11) then
2982 2163 : erelax=0.0_dp
2983 2163 : if(ipert>=1.and.ipert<=natom)then
2984 1353 : erelax=ek0+edocc+eeig0+eloc0+elpsp1+ehart1+exc1+enl0+enl1+epaw1
2985 810 : else if(ipert==natom+1.or.ipert==natom+2)then
2986 : ! NOTE: end0 and evxctau0 are included as "local" in getghc, while
2987 : ! end1 and evxctau1 are included in gvnlx1 (non local) in getgh1c
2988 : ! but in dfpt_vtowfk, end1 and evxctau1 have been subtracted from enl1,
2989 : ! at the same point where ek1 is subtracted from enl1, leaving enl1 "pure"
2990 630 : erelax=ek0+edocc+eeig0+eloc0+ek1+ehart1+exc1+enl0+enl1+epaw1+end0+end1+evxctau0+evxctau1
2991 180 : else if(ipert==natom+3.or.ipert==natom+4)then
2992 168 : erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1
2993 12 : else if(ipert==natom+5.or.ipert==natom+6)then
2994 9 : erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1+elmag1
2995 3 : else if(ipert>natom+11.and.ipert<=2*natom+11)then
2996 3 : erelax=ek0+edocc+eeig0+eloc0+ek1+elpsp1+ehart1+exc1+enl0+enl1+epaw1+elmag1
2997 : end if
2998 2163 : if (abs(emagpen1)>tol8) erelax=erelax+emagpen1
2999 2163 : enl1_effective=enl1
3000 2163 : if (ipert==natom+1.or.ipert==natom+2) then
3001 630 : if (1.0_dp+enl1/10.0_dp==1.0_dp) enl1_effective=zero
3002 : end if
3003 :
3004 2163 : numb='1-9';if (usepaw==1) numb='1-10'
3005 2163 : if( ipert==natom+3.or.ipert==natom+4) then
3006 168 : numb='1-10';if (usepaw==1) numb='1-11'
3007 : end if
3008 : write(msg, '(5a,es17.8)' )&
3009 2163 : & ' ',trim(numb),' gives the relaxation energy (to be shifted if some occ is /=2.0)',&
3010 4326 : & ch10,' erelax=',erelax
3011 2163 : call wrtout(iout,msg)
3012 2163 : call wrtout(std_out,msg)
3013 : end if
3014 :
3015 2342 : if(ipert>=1.and.ipert<=natom)then
3016 :
3017 1353 : numb='10,11,12';if (usepaw==1) numb='11,12,13'
3018 : write(msg, '(4a)' )&
3019 1353 : & ' ',trim(numb),' Non-relaxation contributions : ',&
3020 2706 : & 'frozen-wavefunctions and Ewald'
3021 1353 : call wrtout(iout,msg)
3022 1353 : call wrtout(std_out,msg)
3023 : write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) &
3024 1353 : & ' fr.local=',efrloc,' fr.nonlo=',efrnl,' Ewald=',eew
3025 1353 : call wrtout(iout,msg)
3026 1353 : call wrtout(std_out,msg)
3027 :
3028 1353 : write(msg, '(a,es16.6)' )' dfpt_prtene : non-relax=',efrloc+efrnl+eew
3029 1353 : call wrtout(std_out,msg)
3030 :
3031 1353 : numb='13,14';if (usepaw==1) numb='14,15'
3032 : write(msg, '(3a)' )&
3033 1353 : & ' ',trim(numb),' Frozen wf xc core corrections (1) and (2)'
3034 1353 : call wrtout(iout,msg)
3035 1353 : call wrtout(std_out,msg)
3036 : write(msg, '(a,es17.8,a,es17.8)' ) &
3037 1353 : & ' frxc 1 =',efrx1,' frxc 2 =',efrx2
3038 1353 : call wrtout(iout,msg)
3039 1353 : call wrtout(std_out,msg)
3040 1353 : if (usepaw==1) then
3041 114 : numb='16'
3042 : write(msg, '(5a,es17.8)' )&
3043 114 : & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',&
3044 228 : & ch10,' eovl1 =',eovl1
3045 114 : call wrtout(iout,msg)
3046 114 : call wrtout(std_out,msg)
3047 : end if
3048 1353 : if (usevdw==1) then
3049 7 : numb='15';if (usepaw==1) numb='17'
3050 : write(msg, '(3a,es17.8)' )&
3051 7 : & ' ',trim(numb),' Contribution from van der Waals DFT-D: evdw =',evdw
3052 7 : call wrtout(iout,msg)
3053 7 : call wrtout(std_out,msg)
3054 : end if
3055 :
3056 1353 : write(msg, '(a)' )' Resulting in : '
3057 1353 : call wrtout(iout,msg)
3058 1353 : call wrtout(std_out,msg)
3059 1353 : etotal=erelax+eew+efrloc+efrnl+efrx1+efrx2+evdw
3060 : write(msg, '(a,e20.10,a,e22.12,a)' ) &
3061 1353 : & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV'
3062 1353 : call wrtout(iout,msg)
3063 1353 : call wrtout(std_out,msg)
3064 : write(msg, '(a,es20.10,a,es20.10,a)' ) &
3065 1353 : & ' (2DErelax=',erelax,' Ha. 2DEnonrelax=',etotal-erelax,' Ha)'
3066 1353 : call wrtout(iout,msg)
3067 1353 : call wrtout(std_out,msg)
3068 : write(msg, '(a,es20.10,a,a)' ) &
3069 1353 : & ' ( non-var. 2DEtotal :',&
3070 2706 : & 0.5_dp*(elpsp1+enl1)+eovl1+eew+efrloc+efrnl+efrx1+efrx2+evdw,' Ha)',ch10
3071 1353 : call wrtout(iout,msg)
3072 1353 : call wrtout(std_out,msg)
3073 :
3074 989 : else if(ipert==natom+1.or.ipert==natom+2)then
3075 630 : if (ipert==natom+1.and.usepaw==1) then
3076 129 : numb='11'
3077 : write(msg, '(5a,es17.8)' )&
3078 129 : & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',&
3079 258 : & ch10,' eovl1 =',eovl1
3080 129 : call wrtout(iout,msg)
3081 129 : call wrtout(std_out,msg)
3082 : end if
3083 630 : write(msg,*)' No Ewald or frozen-wf contrib.:',' the relaxation energy is the total one'
3084 630 : if(berry_activated) then
3085 3 : numb='10';
3086 3 : write(msg,'(3a,es20.10)')' ',trim(numb),' Berry phase energy :',eberry
3087 : end if
3088 630 : call wrtout(iout,msg)
3089 630 : call wrtout(std_out,msg)
3090 630 : etotal=erelax
3091 : write(msg, '(a,e20.10,a,e22.12,a)' ) &
3092 630 : & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV'
3093 630 : call wrtout(iout,msg)
3094 630 : call wrtout(std_out,msg)
3095 : write(msg, '(a,es20.10,a)' ) &
3096 630 : & ' ( non-var. 2DEtotal :',0.5_dp*(ek1+enl1_effective+end1+evxctau1)+eovl1,' Ha)'
3097 630 : call wrtout(iout,msg)
3098 630 : call wrtout(std_out,msg)
3099 :
3100 359 : else if(ipert==natom+3 .or. ipert==natom+4) then
3101 168 : numb='11,12,13';if (usepaw==1) numb='12,13,14'
3102 : write(msg, '(4a)' )&
3103 168 : & ' ',trim(numb),' Non-relaxation contributions : ','frozen-wavefunctions and Ewald'
3104 168 : call wrtout(iout,msg)
3105 168 : call wrtout(std_out,msg)
3106 : write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) &
3107 168 : & ' fr.hart=',efrhar,' fr.kin=',efrkin,' fr.loc=',efrloc
3108 168 : call wrtout(iout,msg)
3109 168 : call wrtout(std_out,msg)
3110 :
3111 168 : numb='14,15,16';if (usepaw==1) numb='15,16,17'
3112 : write(msg, '(4a)' )&
3113 168 : & ' ',trim(numb),' Non-relaxation contributions : ','frozen-wavefunctions and Ewald'
3114 168 : call wrtout(iout,msg)
3115 168 : call wrtout(std_out,msg)
3116 : write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) &
3117 168 : & ' fr.nonl=',efrnl,' fr.xc=',efrx1,' Ewald=',eew
3118 168 : call wrtout(iout,msg)
3119 168 : call wrtout(std_out,msg)
3120 :
3121 168 : numb='17';if (usepaw==1) numb='18'
3122 : write(msg, '(4a)' )&
3123 168 : & ' ',trim(numb),' Non-relaxation contributions : ','pseudopotential core energy'
3124 168 : call wrtout(iout,msg)
3125 168 : call wrtout(std_out,msg)
3126 168 : write(msg, '(a,es17.8)' ) ' pspcore=',eii
3127 168 : call wrtout(iout,msg)
3128 168 : call wrtout(std_out,msg)
3129 168 : if (usepaw==1) then
3130 72 : numb='19'
3131 : write(msg, '(5a,es17.8)' )&
3132 72 : & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',&
3133 144 : & ch10,' eovl1 =',eovl1
3134 72 : call wrtout(iout,msg)
3135 72 : call wrtout(std_out,msg)
3136 : end if
3137 168 : if (usevdw==1) then
3138 6 : numb='18';if (usepaw==1) numb='20'
3139 : write(msg, '(3a,es17.8)' )&
3140 6 : & ' ',trim(numb),' Contribution from van der Waals DFT-D: evdw =',evdw
3141 6 : call wrtout(iout,msg)
3142 6 : call wrtout(std_out,msg)
3143 : end if
3144 :
3145 168 : write(msg, '(a,es16.6)' )' dfpt_prtene : non-relax=',&
3146 336 : & efrhar+efrkin+efrloc+efrnl+efrx1+eew+evdw
3147 168 : call wrtout(std_out,msg)
3148 168 : write(msg, '(a)' )' Resulting in : '
3149 168 : call wrtout(iout,msg)
3150 168 : call wrtout(std_out,msg)
3151 168 : etotal=erelax+efrhar+efrkin+efrloc+efrnl+efrx1+eew+eii+evdw
3152 : write(msg, '(a,e20.10,a,e22.12,a)' ) &
3153 168 : & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV'
3154 168 : call wrtout(iout,msg)
3155 168 : call wrtout(std_out,msg)
3156 : write(msg, '(a,es20.10,a,es20.10,a)' ) &
3157 168 : & ' (2DErelax=',erelax,' Ha. 2DEnonrelax=',etotal-erelax,' Ha)'
3158 168 : call wrtout(iout,msg)
3159 168 : call wrtout(std_out,msg)
3160 : write(msg, '(a,es20.10,a,a)' ) &
3161 168 : & ' ( non-var. 2DEtotal :',&
3162 : & 0.5_dp*(elpsp1+enl1+ek1+ehart01)+eovl1+&
3163 336 : & efrhar+efrkin+efrloc+efrnl+efrx1+eew+eii+evdw,' Ha)',ch10
3164 168 : call wrtout(iout,msg)
3165 168 : call wrtout(std_out,msg)
3166 : end if
3167 :
3168 2342 : end subroutine dfpt_prtene
3169 : !!***
3170 :
3171 : !!****f* ABINIT/eigen_meandege
3172 : !! NAME
3173 : !! eigen_meandege
3174 : !!
3175 : !! FUNCTION
3176 : !! This routine takes the mean values of the responses
3177 : !! for the eigenstates that are degenerate in energy.
3178 : !!
3179 : !! INPUTS
3180 : !! eigenresp((3-option)*mband**(3-option)*nkpt*nsppol)= input eigenresp
3181 : !! eigenrep(2*mband**2*nkpt*nsppol) for first-order derivatives of the eigenvalues
3182 : !! eigenrep(mband*nkpt*nsppol) for Fan or Debye-Waller second-order derivatives of the eigenvalues
3183 : !! mband= maximum number of bands
3184 : !! natom= number of atoms in the unit cell
3185 : !! nkpt= number of k-points
3186 : !! nsppol= 1 for unpolarized, 2 for spin-polarized
3187 : !! option= 1 for eigen(1), 2 for eigen(2) - Fan or Debye-Waller
3188 : !!
3189 : !! OUTPUT
3190 : !! eigenresp_mean(mband*nkpt*nsppol)= eigenresp, averaged over degenerate states
3191 : !!
3192 : !! SOURCE
3193 :
3194 533 : subroutine eigen_meandege(eigen0,eigenresp,eigenresp_mean,mband,nband,nkpt,nsppol,option)
3195 :
3196 : !Arguments ------------------------------------
3197 : !scalars
3198 : integer,intent(in) :: mband,nkpt,nsppol,option
3199 : integer,intent(in) :: nband(nkpt*nsppol)
3200 : !arrays
3201 : real(dp),intent(in) :: eigen0(mband*nkpt*nsppol)
3202 : real(dp),intent(in) :: eigenresp((3-option)*mband**(3-option)*nkpt*nsppol)
3203 : real(dp),intent(out) :: eigenresp_mean(mband*nkpt*nsppol)
3204 :
3205 : !Local variables-------------------------------
3206 : !scalars
3207 : integer :: bdtot_index,bd2tot_index,iband,ii,ikpt,isppol,nband_k
3208 : real(dp) :: eig0,mean
3209 : character(len=500) :: msg
3210 : ! *********************************************************************
3211 :
3212 533 : if(option/=1 .and. option/=2)then
3213 0 : write(msg, '(a,i0)' )' The argument option should be 1 or 2, while it is found that option=',option
3214 0 : ABI_BUG(msg)
3215 : end if
3216 :
3217 : bdtot_index=0 ; bd2tot_index=0
3218 1098 : do isppol=1,nsppol
3219 25030 : do ikpt=1,nkpt
3220 23932 : nband_k=nband(ikpt+(isppol-1)*nkpt)
3221 23932 : if(option==1)then
3222 170476 : do iband=1,nband_k
3223 : eigenresp_mean(iband+bdtot_index)=&
3224 170476 : & eigenresp(2*iband-1 + (iband-1)*2*nband_k + bd2tot_index)
3225 : end do
3226 1659 : else if(option==2)then
3227 18095 : do iband=1,nband_k
3228 18095 : eigenresp_mean(iband+bdtot_index)=eigenresp(iband+bdtot_index)
3229 : end do
3230 : end if
3231 :
3232 : ! Treat the case of degeneracies : take the mean of degenerate states
3233 23932 : if(nband_k>1)then
3234 23902 : eig0=eigen0(1+bdtot_index)
3235 23902 : ii=1
3236 164609 : do iband=2,nband_k
3237 140707 : if(eigen0(iband+bdtot_index)-eig0<tol8)then
3238 8124 : ii=ii+1
3239 : else
3240 271285 : mean=sum(eigenresp_mean(iband-ii+bdtot_index:iband-1+bdtot_index))/ii
3241 271285 : eigenresp_mean(iband-ii+bdtot_index:iband-1+bdtot_index)=mean
3242 : ii=1
3243 : end if
3244 140707 : eig0=eigen0(iband+bdtot_index)
3245 164609 : if(iband==nband_k)then
3246 49809 : mean=sum(eigenresp_mean(iband-ii+1+bdtot_index:iband+bdtot_index))/ii
3247 49809 : eigenresp_mean(iband-ii+1+bdtot_index:iband+bdtot_index)=mean
3248 : end if
3249 : end do
3250 : end if
3251 :
3252 23932 : bdtot_index=bdtot_index+nband_k
3253 24497 : bd2tot_index=bd2tot_index+2*nband_k**2
3254 : end do
3255 : end do
3256 :
3257 533 : end subroutine eigen_meandege
3258 : !!***
3259 :
3260 : !!****f* ABINIT/dfpt_init_mag1
3261 : !! NAME
3262 : !! dfpt_init_mag1
3263 : !!
3264 : !! FUNCTION
3265 : !! Initial guess of the first order magnetization/density for magnetic field perturbation.
3266 : !! The first order magnetization is set so as to zero out the first order XC magnetic field, which
3267 : !! should minimize the second order XC energy (without taking self-consistency into account).
3268 : !!
3269 : !! INPUTS
3270 : !! ipert = perturbation type (works only for ipert==natom+5)
3271 : !! idir = direction of the applied magnetic field
3272 : !! cplex = complex or real first order density and magnetization
3273 : !! nfft = dimension of the fft grid
3274 : !! nspden= number of density matrix components
3275 : !! nkxc = number of kxc components
3276 : !! vxc0(nfft,nspden) = GS XC potential
3277 : !! kxc0(nfft,nspden) = GS XC derivatives
3278 : !! rhor0(nfft,nspden) = GS density matrix
3279 : !!
3280 : !! OUTPUT
3281 : !! rhor1(cplex*nfft) = first order density magnetization guess
3282 : !!
3283 : !! SOURCE
3284 :
3285 3 : subroutine dfpt_init_mag1(ipert,idir,rhor1,rhor0,cplex,nfft,nspden,vxc0,kxc0,nkxc)
3286 :
3287 : !Arguments ------------------------------------
3288 : integer, intent(in) :: ipert,idir,cplex,nfft,nspden,nkxc
3289 : real(dp), intent(in) :: vxc0(nfft,nspden),rhor0(nfft,nspden)
3290 : real(dp), intent(in) :: kxc0(nfft,nkxc)
3291 : real(dp), intent(out) :: rhor1(cplex*nfft,nspden)
3292 :
3293 : !Local variables-------------------------------
3294 : integer :: ipt
3295 : real(dp) :: bxc0,bxc1, m1_norm,m0_norm, f_dot_m
3296 : real(dp) :: mdir(3),fdir(3)
3297 : ! *************************************************************************
3298 :
3299 : ABI_UNUSED(ipert)
3300 :
3301 3 : if (nspden==2) then
3302 :
3303 0 : if(cplex==1) then
3304 0 : do ipt=1,nfft
3305 0 : bxc1=half*(half*(kxc0(ipt,1)+kxc0(ipt,3))-kxc0(ipt,2)) ! d/dm Bxc
3306 : !this overestimates the first order magnetization because of n1 not taken into account
3307 0 : m1_norm=-half*(1/bxc1)
3308 0 : rhor1(ipt,1)=zero ! rho_up+rho_dwn => charge density
3309 0 : rhor1(ipt,2)=half*m1_norm ! rho_up=1/2(rho+m) => half*m
3310 : end do
3311 : else
3312 0 : do ipt=1,cplex*nfft
3313 0 : rhor1(ipt,:)=zero
3314 : end do
3315 : end if
3316 :
3317 3 : else if(nspden==4) then
3318 :
3319 : fdir=zero
3320 : fdir(idir)= 1.0d0
3321 7683 : do ipt=1,nfft
3322 7680 : m0_norm=sqrt(rhor0(ipt,2)**2+rhor0(ipt,3)**2+rhor0(ipt,4)**2)
3323 7680 : mdir(1)=rhor0(ipt,2)/m0_norm
3324 7680 : mdir(2)=rhor0(ipt,3)/m0_norm
3325 7680 : mdir(3)=rhor0(ipt,4)/m0_norm
3326 7680 : f_dot_m=fdir(1)*mdir(1)+fdir(2)*mdir(2)+fdir(3)*mdir(3) ! projection of the field direction on m0
3327 :
3328 7680 : bxc1=half*(half*(kxc0(ipt,1)+kxc0(ipt,3))-kxc0(ipt,2)) ! d/dm Bxc
3329 7680 : m1_norm=(-half/bxc1)*f_dot_m ! get an estimate of the norm of m1
3330 :
3331 : bxc0=-sqrt((half*(vxc0(ipt,1)-vxc0(ipt,2)))**2+vxc0(ipt,3)**2+vxc0(ipt,4)**2)
3332 7683 : if(cplex==1) then
3333 : ! TODO:
3334 : ! rhor1(ipt,1)=zero ! rho_up+rho_dwn => charge density
3335 : ! rhor1(ipt,2)=m1_norm*mdir(1)-half*m0_norm/bxc0*(fdir(1)-f_dot_m*mdir(1)) ! m1x
3336 : ! rhor1(ipt,3)=m1_norm*mdir(2)-half*m0_norm/bxc0*(fdir(2)-f_dot_m*mdir(2)) ! m1y
3337 : ! rhor1(ipt,4)=m1_norm*mdir(3)-half*m0_norm/bxc0*(fdir(3)-f_dot_m*mdir(3)) ! m1z
3338 12800 : rhor1(ipt,:)=zero
3339 : else
3340 : ! rhor1(2*ipt-1,1)=zero ! Re rho_up+rho_dwn
3341 : ! rhor1(2*ipt-1,2)=m1_norm*mdir(1)-half*m0_norm/bxc0*(fdir(1)-f_dot_m*mdir(1)) ! m1x
3342 : ! rhor1(2*ipt-1,3)=m1_norm*mdir(2)-half*m0_norm/bxc0*(fdir(2)-f_dot_m*mdir(2)) ! m1x
3343 : ! rhor1(2*ipt-1,4)=m1_norm*mdir(3)-half*m0_norm/bxc0*(fdir(3)-f_dot_m*mdir(3)) ! m1x
3344 : ! rhor1(2*ipt ,1)=zero ! Im rho_up+rho_dwn
3345 : ! rhor1(2*ipt ,2)=zero
3346 : ! rhor1(2*ipt ,3)=zero
3347 : ! rhor1(2*ipt ,4)=zero
3348 :
3349 5120 : rhor1(2*ipt-1,1)=zero; rhor1(2*ipt,1)=zero
3350 5120 : rhor1(2*ipt-1,2)=zero; rhor1(2*ipt,2)=zero
3351 5120 : rhor1(2*ipt-1,3)=zero; rhor1(2*ipt,3)=zero
3352 5120 : rhor1(2*ipt-1,4)=zero; rhor1(2*ipt,4)=zero
3353 :
3354 : end if
3355 : end do
3356 : end if
3357 :
3358 3 : end subroutine dfpt_init_mag1
3359 : !!***
3360 :
3361 : end module m_dfpt_loopert
3362 : !!***
|
