Branch data 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 [ - + - + : 720 : subroutine dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,&
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - - - -
- - - - -
- - - - -
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - - -
- - - - -
- - - - -
- ]
200 : 720 : & ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,&
201 : 1440 : & dyfr_cplex,dyfr_nondiag,d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasval,eigbrd,eig2nkq,&
202 : 720 : & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,&
203 : 1440 : & etotal,fermie,iexit,indsym,kxc,&
204 : 720 : & mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,&
205 : 1440 : & nfftf,nhat,nkpt,nkxc,nspden,nsym,occ,&
206 : 1440 : & paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
207 : 1440 : & pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,&
208 [ + + - + : 1339 : & 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 [ + + + + : 74880 : type(crystal_t) :: crystal, ddb_crystal
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + ]
296 : : type(dataset_type), pointer :: dtset_tmp
297 : 720 : type(ebands_t) :: ebands_k,ebands_kq,gkk_ebands, ebands_kmq !+/-q duplicates
298 : 720 : type(gkk_t) :: gkk2d
299 : 720 : type(hdr_type) :: hdr,hdr_mq,hdr_den,hdr_tmp
300 [ + + + + : 37440 : type(ddb_hdr_type) :: ddb_hdr
+ + + + +
+ + + + +
+ + + + ]
301 : 720 : type(pawang_type) :: pawang1
302 [ + + + + : 92880 : type(wfk_t) :: ddk_f(4)
+ + + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - - + -
+ - + - ]
303 : 720 : 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 : 1440 : integer,allocatable :: blkflg_save(:,:,:,:),dimcprj_srt(:),dyn(:),indkpt1(:),indkpt1_tmp(:)
307 : 720 : integer,allocatable :: indsy1(:,:,:),irrzon1(:,:,:),istwfk_rbz(:),istwfk_pert(:,:,:)
308 : 720 : integer,allocatable :: kg(:,:),kg1(:,:),nband_rbz(:),npwar1(:),npwarr(:),npwtot(:)
309 : 720 : integer,allocatable :: kg1_mq(:,:),npwar1_mq(:),npwtot1_mq(:) !+q/-q duplicates
310 : 720 : integer,allocatable :: npwtot1(:),npwar1_pert(:,:),npwarr_pert(:,:),npwtot_pert(:,:)
311 : 1440 : integer,allocatable :: pert_calc(:,:),pert_tmp(:,:),bz2ibz_smap(:,:)
312 : 720 : integer,allocatable :: symaf1(:),symaf1_tmp(:),symrc1(:,:,:),symrl1(:,:,:),symrl1_tmp(:,:,:)
313 : 720 : integer, pointer :: old_atmtab(:)
314 : 720 : 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 : 720 : real(dp),allocatable :: buffer1(:,:,:,:,:),cg(:,:),cg1(:,:),cg1_active(:,:),cg1_3(:,:,:),cg0_pert(:,:)
317 : 1440 : real(dp),allocatable :: cg1_pert(:,:,:,:),cgq(:,:),gh0c1_pert(:,:,:,:)
318 : 720 : real(dp),allocatable :: doccde_rbz(:),docckqde(:)
319 : 720 : real(dp),allocatable :: gh1c_pert(:,:,:,:),eigen0(:),eigen0_copy(:),eigen1(:),eigen1_mean(:)
320 : 720 : real(dp),allocatable :: eigenq(:),gh1c_set(:,:),gh0c1_set(:,:),kpq(:,:)
321 : 720 : real(dp),allocatable :: kpq_rbz(:,:),kpt_rbz(:,:),occ_pert(:),occ_rbz(:),occkq(:),kpt_rbz_pert(:,:)
322 : 720 : real(dp),allocatable :: occ_disk(:)
323 : 720 : real(dp),allocatable :: vtrial_local(:,:)
324 : 720 : real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons1(:,:,:),resid(:),rhog1(:,:)
325 : 720 : real(dp),allocatable :: rhor1_save(:,:,:)
326 : 720 : real(dp),allocatable :: rhor1(:,:),rho1wfg(:,:),rho1wfr(:,:),tnons1(:,:),tnons1_tmp(:,:)
327 : 720 : real(dp),allocatable :: rhor1_pq(:,:),rhor1_mq(:,:),rhog1_pq(:,:),rhog1_mq(:,:) !+q/-q duplicates
328 : 1440 : real(dp),allocatable :: cg_mq(:,:),cg1_mq(:,:),resid_mq(:) !
329 : 1440 : real(dp),allocatable :: cg1_active_mq(:,:),occk_mq(:) !
330 : 720 : real(dp),allocatable :: kmq(:,:),kmq_rbz(:,:),gh0c1_set_mq(:,:) !
331 : 720 : real(dp),allocatable :: eigen_mq(:),gh1c_set_mq(:,:),docckde_mq(:),eigen1_mq(:) !
332 : 720 : real(dp),allocatable :: vpsp1(:),work(:),wtk_folded(:),wtk_rbz(:),xccc3d1(:)
333 : 720 : real(dp),allocatable :: ylm(:,:),ylm1(:,:),ylmgr(:,:,:),ylmgr1(:,:,:),zeff(:,:,:)
334 : 720 : real(dp),allocatable :: ylm1_mq(:,:),ylmgr1_mq(:,:,:)
335 : 720 : real(dp),allocatable :: phasecg(:,:),gauss(:,:)
336 : 720 : real(dp),allocatable :: gkk(:,:,:,:,:)
337 : : logical :: has_cg1_3(3)
338 : 720 : type(pawcprj_type),allocatable :: cprj(:,:),cprjq(:,:)
339 : 720 : type(paw_ij_type),pointer :: paw_ij_pert(:)
340 : 720 : type(paw_an_type),pointer :: paw_an_pert(:)
341 : 720 : type(pawfgrtab_type),pointer :: pawfgrtab_pert(:)
342 : 720 : type(pawrhoij_type),allocatable :: pawrhoij1(:)
343 : 720 : type(pawrhoij_type),pointer :: pawrhoij_pert(:)
344 : 720 : type(ddb_type) :: ddb
345 : 720 : real(dp),allocatable :: doccde_tmp(:)
346 : : ! ***********************************************************************
347 : :
348 : : DBG_ENTER("COLL")
349 : :
350 : 720 : call timab(141,1,tsec)
351 : :
352 : : !Structured debugging if prtvol==-level
353 [ - + ]: 720 : 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 : 720 : dfpt_scfcv_retcode = -1
359 : 720 : nsppol = dtset%nsppol; nspinor = dtset%nspinor
360 : :
361 : : !Deactivate kramer's degeneracy for calculations that break TRS
362 : 720 : kramers_deg=.true.
363 [ + + ]: 720 : if (dtset%tim1rev==0) kramers_deg=.false.
364 : :
365 [ + + - + ]: 720 : 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 : 720 : with_vxctau = (usevxctau > 0)
373 : 720 : evxctau0=zero; evxctau1=0
374 [ + + + + : 6796 : has_nd = ANY(ABS(dtset%nucdipmom(:,:))>tol8)
+ + ]
375 : 720 : 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 [ - + - - : 720 : call mkrdim(dtset%acell_orig(1:3,1),dtset%rprim_orig(1:3,1:3,1),rprimd)
- - ]
380 : 720 : 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 [ + + + - ]: 1440 : 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 [ + + + + ]: 720 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
389 [ + + ]: 1862 : mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:)
390 : 98 : ecutf=dtset%pawecutdg
391 : : else
392 [ + + ]: 11818 : mgfftf=dtset%mgfft;ngfftf(:)=dtset%ngfft(:)
393 : 622 : ecutf=dtset%ecut
394 : : end if
395 : 720 : ecut_eff=dtset%ecut*(dtset%dilatmx)**2
396 : :
397 : : !Compute large sphere cut-off gsqcut
398 [ + + ]: 720 : if (psps%usepaw==1) then
399 : 116 : call wrtout(std_out,ch10//' FFT (fine) grid used for densities/potentials:')
400 : : end if
401 : 720 : call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,dtset%qptn,ngfftf)
402 : :
403 : : !Various initializations/allocations
404 : 720 : iscf_mod=dtset%iscf
405 : 720 : ntypat=psps%ntypat
406 : 720 : nkpt_max=50;if (xmpi_paral==1) nkpt_max=-1
407 : : !TODO: this flag for paral_atom is ignored below
408 : 720 : paral_atom=(dtset%natom/=my_natom)
409 : 720 : cplex=2-timrev !cplex=2 ! DEBUG: impose cplex=2
410 : 720 : first_entry=.true.
411 : 720 : initialized=0
412 : 720 : ecore=zero ; ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero ; eii=zero
413 [ + + + + : 72432 : d2bbb = zero
+ + + + +
+ + + ]
414 [ + + + + : 2008554 : d2nl = zero
+ + + + +
+ ]
415 [ + + + + : 2008554 : d2lo = zero
+ + + + +
+ ]
416 [ + + + + : 343170 : d2ovl = zero
+ + + + +
+ ]
417 [ + + + + ]: 27684 : clflg(:,:)=0 ! Array on calculated perturbations for eig2rf
418 [ + + ]: 720 : if (psps%usepaw==1) then
419 [ - + + - : 348 : ABI_MALLOC(dimcprj_srt,(dtset%natom))
+ - + - -
+ - + ]
420 : 116 : call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
421 : : end if
422 : 720 : mband_mem_rbz = dtset%mband_mem
423 : :
424 : : !Save values of SCF cycle parameters
425 : 720 : iscf_mod_save = iscf_mod
426 : 720 : nstep_save = dtset%nstep
427 : 720 : nline_save = dtset%nline
428 : 720 : tolwfr_save = dtset%tolwfr
429 : 720 : toldfe_save = dtset%toldfe
430 : 720 : toldff_save = dtset%toldff
431 : 720 : toldmag_save = dtset%toldmag
432 : 720 : tolrff_save = dtset%tolrff
433 : 720 : tolvrs_save = dtset%tolvrs
434 : :
435 : : !This dtset will be used in dfpt_scfcv to force non scf calculations for equivalent perturbations
436 : 720 : nullify(dtset_tmp)
437 [ + + ]: 720 : 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 [ + - - + : 2160 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
+ - + - +
- - + - +
- + ]
446 [ + - - + : 2160 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
+ - + - +
- - + - +
- + ]
447 : 720 : call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),ph1d,xred)
448 [ + + + + ]: 720 : if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then
449 : 98 : 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 [ + + ]: 2880 : rf2_dir1(1:3)=dtset%rf2_pert1_dir(1:3)
485 [ + + ]: 2880 : rf2_dir2(1:3)=dtset%rf2_pert2_dir(1:3)
486 [ + + + + : 5040 : 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 : 682 : rf2dir(1) = rf2_dir1(1)*rf2_dir2(1)
502 : 682 : rf2dir(2) = rf2_dir1(2)*rf2_dir2(2)
503 : 682 : rf2dir(3) = rf2_dir1(3)*rf2_dir2(3)
504 : : ! Upper triangular terms :
505 : 682 : rf2dir(4) = rf2_dir1(2)*rf2_dir2(3)
506 : 682 : rf2dir(5) = rf2_dir1(1)*rf2_dir2(3)
507 : 682 : rf2dir(6) = rf2_dir1(1)*rf2_dir2(2)
508 : : ! Lower triangular terms :
509 : 682 : rf2dir(7) = rf2_dir1(3)*rf2_dir2(2)
510 : 682 : rf2dir(8) = rf2_dir1(3)*rf2_dir2(1)
511 : 682 : 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 [ - + + - : 2160 : ABI_MALLOC(pert_tmp,(3,3*(dtset%natom+6)+18+3*dtset%natom))
+ - + - -
+ - + ]
517 : 720 : ipert_cnt=0
518 [ + + ]: 7461 : do ipert=1,mpert
519 [ + + + + ]: 6741 : if (ipert<dtset%natom+10.or.ipert>dtset%natom+11) then
520 : 26708 : maxidir = 3
521 [ + + ]: 26708 : rfdir(1:3) = dtset%rfdir(:)
522 [ + + ]: 46739 : rfdir(4:9) = 0
523 : : else
524 : 64 : maxidir = 9
525 : 64 : rfdir(1:9) = rf2dir(:)
526 : : end if
527 : : !Scalar potential case
528 [ + + ]: 6741 : if (ipert==dtset%natom+6) then
529 : 720 : maxidir = 1
530 : 720 : rfdir(:) = 0
531 : 720 : rfdir(1) = 1
532 : : end if
533 [ + + ]: 26628 : do idir=1,maxidir
534 : 19167 : to_compute_this_pert = 0
535 [ + + + + : 19167 : if(ipert<dtset%natom+10 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then
+ + ]
536 : : if ((pertsy(idir,ipert)==1).or.&
537 [ + + + + : 3076 : & ((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 [ + + + + : 16091 : else if (ipert==dtset%natom+11 .and. rfpert(ipert)==1 .and. rfdir(idir) == 1 ) then
+ + ]
548 : : to_compute_this_pert = 1
549 [ + + + + ]: 16048 : 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 [ + + + + ]: 15886 : 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 [ + + - + : 15866 : & rfpert(ipert)==1 .and. rfdir(idir) == 1) then
+ + + + ]
577 : : to_compute_this_pert = 1
578 : : end if
579 : 6741 : if (to_compute_this_pert /= 0) then
580 : 2389 : ipert_cnt = ipert_cnt+1;
581 : 2389 : pert_tmp(1,ipert_cnt) = ipert
582 : 2389 : pert_tmp(2,ipert_cnt) = idir
583 : : ! Store "pertcase" in pert_tmp(3,ipert_cnt)
584 [ + + ]: 2389 : if (ipert<dtset%natom+10) then
585 : 2207 : 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 [ - + + - : 2160 : ABI_MALLOC(pert_calc,(3,ipert_cnt))
+ - + - -
+ - + ]
635 [ + + ]: 3109 : do icase=1,ipert_cnt
636 [ + + ]: 10276 : pert_calc(:,icase)=pert_tmp(:,icase)
637 : : end do
638 [ - + ]: 720 : ABI_FREE(pert_tmp)
639 : :
640 [ + + ]: 720 : 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 [ - + + - : 2160 : ABI_MALLOC(eigen0_copy,(dtset%mband*nkpt*dtset%nsppol))
+ - + - -
+ - + ]
648 [ + + ]: 339357 : eigen0_copy(:)=zero
649 : :
650 : : ! SP : Retrieval of the DDB information and computing of effective charge and
651 : : ! dielectric tensor
652 [ - + + - : 2160 : ABI_MALLOC(zeff,(3,3,dtset%natom))
+ - + - -
+ - + ]
653 [ + + + + ]: 720 : 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 [ + + ]: 720 : npert_io=ipert_cnt;if (dtset%nppert<=1) npert_io=0
667 : 720 : 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 [ + + ]: 720 : 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 : 720 : call set_pert_comm(mpi_enreg,dtset%nppert)
678 : :
679 : : !*Redistribute PAW on-site data
680 : 720 : 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 : 720 : & paw_an,paw_ij,pawfgrtab,pawrhoij)
684 : 720 : pawfgrtab_pert=>pawfgrtab ; pawrhoij_pert=>pawrhoij
685 : 720 : 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 [ + + + + ]: 720 : 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 [ + + ]: 3109 : do icase=1,ipert_cnt
702 : :
703 : : ! %%%% Parallelization over perturbations %%%%%
704 : : ! Select the perturbations treated by current processor
705 [ + + ]: 2389 : 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 : 2354 : call localwrfile(mpi_enreg%comm_cell,icase,npert_io,mpi_enreg%paral_pert,0)
711 : :
712 : : ! Set precision for FFT libs.
713 : 2354 : 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 : 2354 : ipert = pert_calc(1,icase)
731 : 2354 : idir = pert_calc(2,icase)
732 : 2354 : istr=idir
733 : 2354 : 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 : 2354 : spaceComm=mpi_enreg%comm_cell
751 : 2354 : me=mpi_enreg%me_cell
752 : :
753 : : ! ===== Describe the perturbation in output/log file
754 : :
755 [ + + ]: 193028 : write(msg, '(a,80a,a,a,3f10.6)' ) ch10,('-',ii=1,80),ch10,&
756 [ + - ]: 195382 : ' Perturbation wavevector (in red.coord.) ',dtset%qptn(:)
757 [ + + ]: 7062 : call wrtout([std_out, ab_out],msg)
758 [ + + ]: 2354 : 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 [ + - + + ]: 2354 : 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 [ + + ]: 1001 : else if(ipert==dtset%natom+1)then
775 : 448 : write(msg,'(a,i4)')' Perturbation : derivative vs k along direction',idir
776 [ + + ]: 1344 : call wrtout([std_out, ab_out], msg)
777 [ + + ]: 448 : 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 [ + - ]: 448 : 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 [ + + ]: 2354 : 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 : 2354 : eew=zero ; evdw=zero ; efrloc=zero ; efrnl=zero ; efrx1=zero ; efrx2=zero
865 : 2354 : efrhar=zero ; efrkin=zero
866 [ + + ]: 2354 : 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 [ + + + + ]: 1001 : 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 [ - + + - : 7062 : ABI_MALLOC(symaf1_tmp,(nsym))
+ - + - -
+ - + ]
891 [ - + + - : 7062 : ABI_MALLOC(symrl1_tmp,(3,3,nsym))
+ - + - -
+ - + ]
892 [ - + + - : 7062 : 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 [ + + + + ]: 2354 : & 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 : 2099 : & 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 [ - + + - : 9416 : ABI_MALLOC(indsy1,(4,nsym1,dtset%natom))
- + + - +
- + - - +
- + ]
906 [ - + + - : 7062 : ABI_MALLOC(symrc1,(3,3,nsym1))
+ - + - -
+ - + ]
907 [ - + + - : 7062 : ABI_MALLOC(symaf1,(nsym1))
+ - + - -
+ - + ]
908 [ - + + - : 4708 : ABI_MALLOC(symrl1,(3,3,nsym1))
- + - + ]
909 [ - + + - : 7062 : ABI_MALLOC(tnons1,(3,nsym1))
+ - + - -
+ - + ]
910 [ + + ]: 6413 : symaf1(1:nsym1)=symaf1_tmp(1:nsym1)
911 [ + + + + : 55121 : symrl1(:,:,1:nsym1)=symrl1_tmp(:,:,1:nsym1)
+ + ]
912 [ + + + + ]: 18590 : tnons1(:,1:nsym1)=tnons1_tmp(:,1:nsym1)
913 [ - + ]: 2354 : ABI_FREE(symaf1_tmp)
914 [ - + ]: 2354 : ABI_FREE(symrl1_tmp)
915 [ - + ]: 2354 : ABI_FREE(tnons1_tmp)
916 : :
917 : : ! Set up corresponding symmetry data
918 [ - + + - : 9416 : ABI_MALLOC(irrzon1,(dtset%nfft**(1-1/nsym1),2,(nspden/dtset%nsppol)-3*(nspden/4)))
- + + - +
- + - - +
- + ]
919 [ - + + - : 9416 : 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 : 2354 : & nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred)
922 [ + + ]: 2354 : if (psps%usepaw==1) then
923 : : ! Allocate/initialize only zarot in pawang1 datastructure
924 : 409 : call pawang_init(pawang1,0,0,pawang%l_max-1,0,0,nsym1,0,0,0,0)
925 : 409 : 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 [ - + + - : 7062 : ABI_MALLOC(kpq,(3,nkpt))
+ - + - -
+ - + ]
930 [ + + + + ]: 2354 : 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 [ + + ]: 155095 : do ikpt=1,nkpt
934 [ + + ]: 613822 : 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 [ + + ]: 2354 : 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 [ - + + - : 7062 : ABI_MALLOC(indkpt1_tmp,(nkpt))
+ - + - -
+ - + ]
951 [ - + + - : 7062 : ABI_MALLOC(wtk_folded,(nkpt))
+ - + - -
+ - + ]
952 [ - + + - : 7062 : ABI_MALLOC(bz2ibz_smap, (6, nkpt))
+ - + - -
+ - + ]
953 [ + + ]: 174667 : indkpt1_tmp(:)=0 ; optthm=0
954 : 2354 : timrev_pert=timrev
955 [ + + ]: 2354 : 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 [ + + + + : 1860 : & (ipert>dtset%natom+11.and.ipert<=2*dtset%natom+11).or. &
+ + + + +
+ + + + +
+ + - + +
+ ]
969 : 840 : & dtset%prtfull1wf==1) timrev_pert=0
970 : 1860 : timrev_kpt = timrev_pert
971 : :
972 : : !MR: Modified to agree with longwave driver
973 [ + + ]: 1860 : 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 [ + + ]: 1860 : 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 : 1860 : nsym1,symrc1,timrev_kpt,dtset%wtk,wtk_folded, bz2ibz_smap, xmpi_comm_self)
996 : : end if
997 : :
998 : 2354 : write (msg, "(a,i10,a,i10,a)") " dfpt_looppert : the # of k-points has been reduced from ", &
999 : 4708 : nkpt, " to ", nkpt_rbz, " using the little group "
1000 : 2354 : call wrtout(std_out, msg)
1001 : :
1002 [ - + + - : 7062 : ABI_MALLOC(doccde_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol))
+ - + - -
+ - + ]
1003 [ - + + - : 7062 : ABI_MALLOC(indkpt1,(nkpt_rbz))
+ - + - -
+ - + ]
1004 [ - + + - : 4708 : ABI_MALLOC(istwfk_rbz,(nkpt_rbz))
- + - + ]
1005 [ - + + - : 7062 : ABI_MALLOC(kpq_rbz,(3,nkpt_rbz))
+ - + - -
+ - + ]
1006 [ + + ]: 2354 : if (.not.kramers_deg) then
1007 [ - + + - : 22 : ABI_MALLOC(kmq_rbz,(3,nkpt_rbz))
- + - + -
+ ]
1008 : : end if
1009 [ - + + - : 4708 : ABI_MALLOC(kpt_rbz,(3,nkpt_rbz))
- + - + ]
1010 [ - + + - : 7062 : ABI_MALLOC(nband_rbz,(nkpt_rbz*dtset%nsppol))
+ - + - -
+ - + ]
1011 [ - + + - : 4708 : ABI_MALLOC(occ_rbz,(dtset%mband*nkpt_rbz*dtset%nsppol))
- + - + ]
1012 [ - + + - : 4708 : ABI_MALLOC(occ_disk,(dtset%mband*nkpt_rbz*dtset%nsppol))
- + - + ]
1013 [ - + + - : 7062 : ABI_MALLOC(wtk_rbz,(nkpt_rbz))
+ - + - -
+ - + ]
1014 [ + + ]: 150052 : indkpt1(:)=indkpt1_tmp(1:nkpt_rbz)
1015 [ + + ]: 150052 : do ikpt=1,nkpt_rbz
1016 : 147698 : istwfk_rbz(ikpt)=dtset%istwfk(indkpt1(ikpt))
1017 [ + + ]: 590792 : kpq_rbz(:,ikpt)=kpq(:,indkpt1(ikpt))
1018 [ + + ]: 590792 : kpt_rbz(:,ikpt)=dtset%kptns(:,indkpt1(ikpt))
1019 : 150052 : wtk_rbz(ikpt)=wtk_folded(indkpt1(ikpt))
1020 : : end do
1021 [ + + ]: 2354 : 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 [ - + ]: 2354 : ABI_FREE(indkpt1_tmp)
1027 [ - + ]: 2354 : 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 : 2354 : bdtot_index=0;bdtot1_index=0
1037 [ + + ]: 4875 : do isppol=1,dtset%nsppol
1038 : 2521 : ikpt1=1
1039 [ + + ]: 179102 : do ikpt=1,nkpt
1040 : 174227 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
1041 : : ! Must test against ikpt1/=nkpt_rbz+1, before evaluate indkpt1(ikpt1)
1042 [ + + ]: 174227 : if(ikpt1/=nkpt_rbz+1)then
1043 [ + + ]: 167885 : if(ikpt==indkpt1(ikpt1))then
1044 : 149462 : nband_rbz(ikpt1+(isppol-1)*nkpt_rbz)=nband_k
1045 [ + + ]: 987110 : occ_rbz(1+bdtot1_index:nband_k+bdtot1_index) = occ(1+bdtot_index:nband_k+bdtot_index)
1046 [ + + ]: 987110 : doccde_rbz(1+bdtot1_index:nband_k+bdtot1_index) = doccde(1+bdtot_index:nband_k+bdtot_index)
1047 : 149462 : ikpt1=ikpt1+1
1048 : 149462 : bdtot1_index=bdtot1_index+nband_k
1049 : : end if
1050 : : end if
1051 : 176748 : bdtot_index=bdtot_index+nband_k
1052 : : end do
1053 : : end do
1054 : :
1055 : : ! Compute maximum number of planewaves at k
1056 : 2354 : call timab(142,1,tsec)
1057 : 2354 : call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpt_rbz,mpi_enreg,mpw,nkpt_rbz)
1058 : 2354 : call timab(142,2,tsec)
1059 : :
1060 : : ! Allocate some k-dependent arrays at k
1061 [ - + + - : 7062 : ABI_MALLOC(npwarr,(nkpt_rbz))
+ - + - -
+ - + ]
1062 [ - + + - : 4708 : ABI_MALLOC(npwtot,(nkpt_rbz))
- + - + ]
1063 : :
1064 : : ! Determine distribution of k-points/bands over MPI processes
1065 [ + + ]: 2354 : if (allocated(mpi_enreg%my_kpttab)) then
1066 : 720 : ABI_FREE(mpi_enreg%my_kpttab)
1067 : : end if
1068 [ - + + - : 7062 : ABI_MALLOC(mpi_enreg%my_kpttab,(nkpt_rbz))
+ - + - -
+ - + -
+ ]
1069 : : if(xmpi_paral==1) then
1070 [ - + + - : 11770 : ABI_MALLOC(mpi_enreg%proc_distrb,(nkpt_rbz,dtset%mband,dtset%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
1071 : 2354 : 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 [ + + ]: 150052 : my_nkpt_rbz=maxval(mpi_enreg%my_kpttab)
1076 : 2354 : mkmem_rbz =my_nkpt_rbz ; mkqmem_rbz=my_nkpt_rbz ; mk1mem_rbz=my_nkpt_rbz
1077 : : ABI_UNUSED((/mkmem,mk1mem,mkqmem/))
1078 : :
1079 : 2354 : 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 [ - + + - : 11770 : ABI_MALLOC(distrb_flags,(nkpt_rbz,dtset%mband,dtset%nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
1083 [ + - + - : 866220 : distrb_flags = (mpi_enreg%proc_distrb == mpi_enreg%me_kpt)
+ - - + -
- + + + +
+ + ]
1084 : :
1085 : : ! Set up the basis sphere of planewaves at k
1086 [ - + + - : 7062 : ABI_MALLOC(kg,(3,mpw*mkmem_rbz))
+ - + - -
+ - + ]
1087 : 2354 : 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 : 2354 : mpw,npwarr,npwtot,dtset%nsppol)
1091 : 2354 : call timab(143,2,tsec)
1092 : :
1093 : : ! Set up the spherical harmonics (Ylm) at k
1094 : 2354 : useylmgr=0; option=0 ; nylmgr=0
1095 [ + + + + : 2354 : 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 : 330 : 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 [ - + + - : 9416 : ABI_MALLOC(ylm,(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
+ + + - +
- + + - +
- + - + ]
1103 [ - + + - : 11770 : ABI_MALLOC(ylmgr,(mpw*mkmem_rbz,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
+ + + - +
+ + - + -
+ + - + -
+ - + ]
1104 [ + + ]: 2354 : if (psps%useylm==1) then
1105 : : call initylmg(gprimd,kg,kpt_rbz,mkmem_rbz,mpi_enreg,psps%mpsang,mpw,nband_rbz,nkpt_rbz,&
1106 : 720 : npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr)
1107 : : end if
1108 : :
1109 : : ! Set up occupations for this perturbation
1110 [ + + ]: 2354 : 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 [ + + ]: 2354 : 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 : 2354 : write(msg,'(3a)')ch10,'--------------------------------------------------------------------------------',ch10
1129 : 2354 : 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 [ + + ]: 151816 : bantot_rbz=sum(nband_rbz(1:nkpt_rbz*dtset%nsppol))
1135 [ - + + - : 7062 : ABI_MALLOC(eigen0,(bantot_rbz))
+ - + - -
+ - + -
+ ]
1136 [ + + ]: 840002 : 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 : 2354 : dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1142 [ - + ]: 2354 : ABI_FREE(eigen0)
1143 : :
1144 : : ! Initialize header, update it with evolving variables
1145 : 2354 : gscase=0 ! A GS WF file is read
1146 : : call hdr0%init(ebands_k,codvsn,dtset,pawtab,gscase,psps,wvl%descr,&
1147 : 2354 : 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 : 2354 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1152 : :
1153 : : ! Initialize GS wavefunctions at k
1154 : 2354 : ireadwf0=1; formeig=0 ; ask_accurate=1 ; optorth=0
1155 : 2354 : 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 : 2354 : ftoa(3 * two * mcg * dp * b2Mb, fmt="f8.1"), "[Mb] <<< MEM"))
1159 : :
1160 [ - + ]: 2354 : 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 [ - + + - : 7062 : ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr)
+ - + - +
- - + -
+ ]
1170 : :
1171 [ - + + - : 7062 : ABI_MALLOC(eigen0,(dtset%mband*nkpt_rbz*dtset%nsppol))
+ - + - -
+ - + ]
1172 : 2354 : call timab(144,1,tsec)
1173 : 2354 : 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 : 2354 : cg, eigen=eigen0, occ=occ_disk)
1180 : :
1181 : 2354 : call timab(144,2,tsec)
1182 : :
1183 : : ! Update energies GS energies at k
1184 : 2354 : call put_eneocc_vect(ebands_k, "eig", eigen0)
1185 : :
1186 : : ! PAW: compute on-site projections of GS wavefunctions (cprj) (and derivatives) at k
1187 : 2354 : ncpgr=0
1188 [ - + - + ]: 2354 : ABI_MALLOC(cprj,(0,0))
1189 [ + + ]: 2354 : if (psps%usepaw==1) then
1190 : 409 : ncpgr=3 ! Valid for ipert<=natom (phonons), ipert=natom+2 (elec. field)
1191 : : ! or for ipert==natom+10,11
1192 [ + + ]: 409 : if (ipert==dtset%natom+1) then
1193 [ + + ]: 141 : if (dtset%orbmag.NE.0) then
1194 : : ncpgr=3
1195 : : else
1196 : 84 : ncpgr=1
1197 : : end if
1198 : : end if
1199 [ + + + + ]: 409 : if (ipert==dtset%natom+3.or.ipert==dtset%natom+4) ncpgr=1
1200 [ + - ]: 409 : 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 : 409 : mcprj=dtset%nspinor*mband_mem_rbz*mkmem_rbz*dtset%nsppol
1205 : : !mcprj=dtset%nspinor*dtset%mband*mkmem_rbz*dtset%nsppol
1206 [ - + ]: 409 : ABI_FREE(cprj)
1207 [ - + + - : 243192 : ABI_MALLOC(cprj,(dtset%natom,mcprj))
- + + - +
- + - - +
- + + + +
+ ]
1208 : 409 : call pawcprj_alloc(cprj,ncpgr,dimcprj_srt)
1209 [ + + ]: 409 : if (ipert<=dtset%natom) then
1210 : 114 : choice=2; iorder_cprj=0; idir0=0
1211 [ + + ]: 295 : else if (ipert==dtset%natom+1) then
1212 [ + + ]: 141 : if (dtset%orbmag.NE.0) then
1213 : 57 : 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 : 409 : & 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 : 2354 : call timab(143,1,tsec)
1237 : 2354 : call getmpw(ecut_eff,dtset%exchn2n3d,gmet,istwfk_rbz,kpq_rbz,mpi_enreg,mpw1,nkpt_rbz)
1238 [ + + ]: 2354 : 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 : 2343 : mpw1_mq=0
1245 : : end if
1246 : 2354 : call timab(143,2,tsec)
1247 : :
1248 : : ! Allocate some arrays at k+q
1249 [ - + + - : 7062 : ABI_MALLOC(kg1,(3,mpw1*mk1mem_rbz))
+ - + - -
+ - + ]
1250 [ - + + - : 7062 : ABI_MALLOC(npwar1,(nkpt_rbz))
+ - + - -
+ - + ]
1251 [ - + + - : 4708 : ABI_MALLOC(npwtot1,(nkpt_rbz))
- + - + ]
1252 : : ! In case Kramers degeneracy is broken, do the same for k-q
1253 [ + + ]: 2354 : 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 : 2354 : 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 : 2354 : npwar1,npwtot1,dtset%nsppol)
1265 [ + + ]: 2354 : 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 : 2354 : call timab(142,2,tsec)
1271 : :
1272 : : ! Set up the spherical harmonics (Ylm) at k+q
1273 : 2354 : useylmgr1=0; option=0 ; nylmgr1=0
1274 [ + + + + : 2354 : 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 : 330 : 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 [ - + + - : 9416 : ABI_MALLOC(ylm1,(mpw1*mk1mem_rbz,psps%mpsang*psps%mpsang*psps%useylm))
+ + + - +
- + + - +
- + - + ]
1282 [ - + + - : 11770 : ABI_MALLOC(ylmgr1,(mpw1*mk1mem_rbz,nylmgr1,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
+ + + - +
+ + - + -
+ + - + -
+ - + ]
1283 [ + + ]: 2354 : if (psps%useylm==1) then
1284 : : call initylmg(gprimd,kg1,kpq_rbz,mk1mem_rbz,mpi_enreg,psps%mpsang,mpw1,nband_rbz,nkpt_rbz,&
1285 : 720 : npwar1,dtset%nsppol,option,rprimd,ylm1,ylmgr1)
1286 : : end if
1287 [ + + ]: 2354 : 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 : 2354 : write(msg, '(a,a)' )'--------------------------------------------------------------------------------',ch10
1303 : 2354 : call wrtout(ab_out,msg)
1304 : :
1305 : : ! Initialize band structure datatype at k+q
1306 [ - + + - : 4708 : ABI_MALLOC(eigenq,(bantot_rbz))
- + - + -
+ ]
1307 [ - + + - : 7062 : ABI_MALLOC(doccde_tmp,(dtset%mband*nkpt_rbz*dtset%nsppol))
+ - + - -
+ - + ]
1308 [ + + ]: 840002 : 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 : 2354 : & dtset%kptrlatt, dtset%nshiftk, dtset%shiftk)
1314 [ + + ]: 2354 : 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 [ - + ]: 2354 : ABI_FREE(eigenq)
1323 [ - + ]: 2354 : ABI_FREE(doccde_tmp)
1324 : :
1325 : : ! Initialize header
1326 : : call hdr%init(ebands_kq,codvsn,dtset,pawtab,pertcase,psps,wvl%descr, &
1327 : 2354 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab )
1328 [ + + ]: 2354 : 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 : 2354 : ireadwf0=1 ; formeig=0 ; ask_accurate=1 ; optorth=0
1339 : 2354 : mcgq=mpw1*dtset%nspinor*mband_mem_rbz*mkqmem_rbz*dtset%nsppol
1340 : : !SPr: verified until here, add mcgq for -q
1341 [ - + ]: 2354 : 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 [ - + + - : 7062 : ABI_MALLOC_OR_DIE(cgq,(2,mcgq), ierr)
+ - + - +
- - + -
+ ]
1351 : :
1352 [ - + + - : 7062 : ABI_MALLOC(eigenq,(dtset%mband*nkpt_rbz*dtset%nsppol))
+ - + - -
+ - + ]
1353 [ + + ]: 2354 : 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 : 2354 : 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 [ + + + + : 9416 : if (dtfil%fnamewffq == dtfil%fnamewffk .and. sum(dtset%qptn(1:3)**2) < 1.d-14) then
+ + ]
1367 : 1545 : call wrtout(std_out, " qpt is Gamma --> psi_k+q initialized from psi_k in memory")
1368 [ + - - + : 161715954 : cgq = cg
- - + + +
+ ]
1369 [ + - - + : 407576 : 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 [ + + ]: 2354 : 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 [ - + ]: 2354 : ABI_FREE(occ_disk)
1396 : :
1397 : : ! Update energies GS energies at k + q
1398 : 2354 : call put_eneocc_vect(ebands_kq, "eig", eigenq)
1399 [ + + ]: 2354 : 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 [ - + - + ]: 2354 : ABI_MALLOC(cprjq,(0,0))
1403 [ + + ]: 2354 : if (psps%usepaw==1) then
1404 [ + - ]: 409 : if (usecprj==1) then
1405 : 409 : 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 [ - + ]: 409 : ABI_FREE(cprjq)
1409 [ - + + - : 243192 : ABI_MALLOC(cprjq,(dtset%natom,mcprjq))
- + + - +
- + - - +
- + + + +
+ ]
1410 : 409 : call pawcprj_alloc(cprjq,0,dimcprj_srt)
1411 [ + + + + : 1636 : 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 [ + - ]: 373 : else if (mcprjq>0) then
1419 : 373 : call pawcprj_copy(cprj,cprjq)
1420 : : end if
1421 : : end if
1422 : : end if
1423 : :
1424 : : ! ===== Report on eigenq values
1425 [ + + + + ]: 2354 : 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 [ + + + + ]: 2354 : 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 : 2354 : nkpt_eff=nkpt
1435 [ + + ]: 2354 : if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. nkpt>nkpt_max ) nkpt_eff=nkpt_max
1436 : 2354 : band_index=0
1437 [ + + ]: 4875 : do isppol=1,dtset%nsppol
1438 [ + + ]: 154337 : do ikpt=1,nkpt_rbz
1439 : 149462 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
1440 [ + + ]: 149462 : 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 [ - + ]: 126844 : 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 : 151983 : 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 [ - + + - : 7062 : ABI_MALLOC(docckqde,(dtset%mband*nkpt_rbz*dtset%nsppol))
+ - + - -
+ - + ]
1457 [ - + + - : 4708 : ABI_MALLOC(occkq,(dtset%mband*nkpt_rbz*dtset%nsppol))
- + - + ]
1458 [ + + ]: 2354 : 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 [ + + ]: 2354 : if(0<=dtset%occopt .and. dtset%occopt<=2)then
1464 : : ! Same occupation numbers at k and k+q (usually, insulating)
1465 [ + + ]: 537701 : occkq(:)=occ_rbz(:)
1466 [ + + ]: 537701 : docckqde(:)=zero ! docckqde is irrelevant in this case
1467 [ + + ]: 1867 : 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 [ - + ]: 2354 : if(dtset%prtvol==-level) call wrtout(std_out,'dfpt_looppert: initialisation of q part done.')
1498 : :
1499 : : ! Initialisation of first-order wavefunctions
1500 : 2354 : write(msg,'(3a,i4)')' Initialisation of the first-order wave-functions :',ch10,' ireadwf=',dtfil%ireadwf
1501 [ + + ]: 7062 : call wrtout([std_out, ab_out], msg)
1502 : 2354 : call appdig(pertcase,dtfil%fnamewff1,fiwf1i)
1503 : 2354 : call appdig(pertcase,dtfil%fnameabo_1wf,fiwf1o)
1504 [ + + ]: 2354 : 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 [ + + ]: 2354 : if (psps%usepaw==1) then
1511 [ + + + - : 2011 : ABI_MALLOC(pawrhoij1,(my_natom))
+ - + + -
+ - + - +
+ + + - ]
1512 : 409 : call pawrhoij_nullify(pawrhoij1)
1513 : : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,&
1514 : 409 : & 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 : 409 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1518 [ + - - + ]: 818 : 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 : 2354 : 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 : 2354 : dim_eig2rf=0
1533 [ + + + + ]: 2354 : if ((dtset%ieig2rf > 0 .and. dtset%ieig2rf/=2) .or. dtset%efmas > 0) then
1534 : 543 : dim_eig2rf=1
1535 : : end if
1536 : 2354 : mcg1=mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol
1537 [ - + ]: 2354 : 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 [ - + + - : 7062 : 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 [ + + - + : 2354 : & .AND. (.NOT. ALLOCATED(cg1_3)) ) then
- - + + ]
1551 [ - + + - : 80 : ABI_MALLOC(cg1_3,(2,mcg1,3))
+ - + - +
- - + - +
- + ]
1552 : 20 : has_cg1_3(:) = .FALSE.
1553 : : end if
1554 [ + + ]: 2354 : 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 [ + + + - : 7062 : ABI_MALLOC(cg1_active,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
+ - + + -
+ - + -
+ ]
1560 [ + + + - : 4708 : ABI_MALLOC(gh1c_set,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
- + - + -
+ ]
1561 [ + + + - : 4708 : ABI_MALLOC(gh0c1_set,(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem_rbz*dtset%nsppol*dim_eig2rf))
- + - + -
+ ]
1562 [ + + ]: 2354 : 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 [ + - + + ]: 2354 : 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 [ + - - + : 7062 : ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol))
+ - + - +
- - + - +
- + ]
1589 [ - + + - : 7062 : ABI_MALLOC(resid,(dtset%mband*nkpt_rbz*dtset%nsppol))
+ - + - -
+ - + ]
1590 : 2354 : call timab(144,1,tsec)
1591 [ + + - + : 2354 : 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 [ + + + + ]: 429582609 : cg1 = zero
1598 [ + + ]: 12951835 : eigen1 = zero
1599 : : end if
1600 : :
1601 : 2354 : call timab(144,2,tsec)
1602 [ + + ]: 2354 : 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 [ + + - + ]: 2354 : 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 [ + + + + : 9416 : 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 : 2934 : 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 [ + + ]: 2354 : n3xccc=0;if(psps%n1xccc/=0)n3xccc=nfftf
1692 [ + + + - : 7062 : ABI_MALLOC(xccc3d1,(cplex*n3xccc))
+ - + + -
+ - + ]
1693 [ - + + - : 7062 : ABI_MALLOC(vpsp1,(cplex*nfftf))
+ - + - -
+ - + ]
1694 : :
1695 : : ! PAW: compute Vloc(1) and core(1) together in reciprocal space
1696 : : ! --------------------------------------------------------------
1697 [ + + + + ]: 2354 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
1698 : 411 : 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 : 411 : & 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 [ + + + + ]: 13876898 : eigen1(:)=zero; resid(:)=zero
1734 [ + + ]: 2354 : 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 [ - + + - : 9416 : ABI_MALLOC(rhor1,(cplex*nfftf,nspden))
- + + - +
- + - - +
- + ]
1740 [ - + + - : 7062 : ABI_MALLOC(rhog1,(2,nfftf))
+ - + - -
+ - + ]
1741 : :
1742 : : ! MR:do we really need to define the partial rhos here?
1743 [ + + ]: 2354 : 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 : 2354 : found_eq_gkk=.false.
1753 [ + + + + ]: 2354 : 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 [ + + ]: 2354 : 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 [ + + + - : 2354 : .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 [ + + + + : 146824938 : 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 [ + - + + : 2174 : & .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 [ + + ]: 2174 : 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 [ + - ]: 2354 : openexit=1 ; if(dtset%chkexit==0) openexit=0
1959 : 2354 : 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 [ - + - - ]: 2354 : if (iexit/=0) eigen1(:)=zero
1962 [ - + - - : 2354 : if (iexit/=0.and.(.not.kramers_deg)) eigen1_mq(:)=zero
- - ]
1963 : 2354 : call pstat_proc%print(_PSTAT_ARGS_)
1964 : :
1965 [ + - ]: 2354 : 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 [ + + ]: 2354 : 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 [ - + - + : 2343 : & 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 [ + + ]: 2354 : 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 [ + + ]: 2354 : 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 [ - + ]: 2354 : ABI_FREE(gh1c_set)
2079 [ - + ]: 2354 : ABI_FREE(gh0c1_set)
2080 [ - + ]: 2354 : ABI_FREE(cg1_active)
2081 : :
2082 : : !deallocate big arrays for case without Kramers' degeneracy
2083 [ + + ]: 2354 : 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 : 2354 : call pstat_proc%print(_PSTAT_ARGS_)
2092 : 2354 : call timab(146,1,tsec)
2093 : :
2094 : : ! Print out message at the end of the iterations
2095 [ + + ]: 190674 : write(msg, '(80a,a,a,a,a)' ) ('=',ii=1,80),ch10,ch10,&
2096 [ + - ]: 193028 : & ' ----iterations are completed or convergence reached----',ch10
2097 [ + + ]: 7062 : call wrtout([std_out, ab_out], msg)
2098 : :
2099 : : ! Print _gkk file for this perturbation
2100 [ + + ]: 2354 : 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 [ + + + + ]: 2354 : 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 : 2354 : write_1wfk = .True.
2189 [ + + ]: 2354 : if (dtset%prtwf == 0) write_1wfk = .False.
2190 [ + + + - ]: 2354 : 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 [ + + + - ]: 2354 : if( (dtset%orbmag .NE. 0) .AND. (ipert .EQ. dtset%natom+1) ) then
2198 [ + + + + ]: 8978172 : cg1_3(:,:,idir) = cg1(:,:)
2199 : 60 : has_cg1_3(idir) = .TRUE.
2200 : : end if
2201 : :
2202 : 2354 : call outresid(dtset,kpt_rbz,dtset%mband, nband_rbz,nkpt_rbz,dtset%nsppol,resid)
2203 : :
2204 [ + + ]: 2354 : 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 : 1882 : & cg1, kg1, eigen1)
2209 : :
2210 [ + + ]: 1882 : 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 [ + + + + : 2354 : 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 [ + + ]: 2354 : 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 : 448 : ek2=gmet(idir,idir)*(two_pi**2)*2.0_dp*dtset%nelect
2254 : 448 : fsum=-ek1/ek2
2255 [ + + ]: 448 : if(dtset%ecutsm<tol6)then
2256 : : write(msg, '(a,es20.10,a,a,es20.10)' ) &
2257 : 364 : ' dfpt_looppert : ek2=',ek2,ch10,&
2258 : 728 : ' 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 [ + + ]: 1344 : call wrtout([std_out, ab_out] , msg)
2265 : : ! Write the diagonal elements of the dH/dk operator, after averaging over degenerate states
2266 [ - + + - : 1344 : ABI_MALLOC(eigen1_mean,(dtset%mband*nkpt_rbz*dtset%nsppol))
+ - + - -
+ - + ]
2267 : 448 : call eigen_meandege(eigen0,eigen1,eigen1_mean,dtset%mband,nband_rbz,nkpt_rbz,dtset%nsppol,1)
2268 : 448 : option=4
2269 [ + + ]: 448 : 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 : 400 : & option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk_rbz)
2273 : : end if
2274 [ - + ]: 448 : ABI_FREE(eigen1_mean)
2275 : : end if
2276 : :
2277 : : ! Print the energies
2278 [ + + + - ]: 2354 : 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 : 2354 : & 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 [ + + + + : 9416 : & ( (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 [ - + + - : 80 : ABI_MALLOC(vtrial_local,(nfftf,dtset%nspden))
- + + - +
- + - - +
- + - + ]
2291 : : end if
2292 [ + - + - : 5399229 : vtrial_local = vtrial
- + - - +
+ + + ]
2293 [ + + ]: 20 : 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 : 19 : & 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 [ + - ]: 20 : ABI_SFREE(vtrial_local)
2303 [ + - ]: 20 : ABI_SFREE(cg1_3)
2304 : 20 : has_cg1_3(:) = .FALSE.
2305 : : end if ! end call orbmag
2306 : :
2307 [ + + ]: 2354 : 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 [ + + ]: 2315 : if (icase == ipert_cnt) then
2316 [ + + ]: 258871 : eigen0_copy(1:dtset%mband*nkpt_rbz*dtset%nsppol) = eigen0
2317 : : else
2318 : 1621 : 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 [ - + ]: 2354 : ABI_FREE(cg)
2328 [ - + ]: 2354 : ABI_FREE(cgq)
2329 [ - + ]: 2354 : ABI_FREE(cg1)
2330 [ - + ]: 2354 : ABI_FREE(docckqde)
2331 [ + + ]: 2354 : 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 [ - + ]: 2354 : ABI_FREE(doccde_rbz)
2337 [ - + ]: 2354 : ABI_FREE(eigen0)
2338 [ - + ]: 2354 : ABI_FREE(eigenq)
2339 [ - + ]: 2354 : ABI_FREE(eigen1)
2340 [ - + ]: 2354 : ABI_FREE(kpq)
2341 [ + + ]: 2354 : 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 [ - + ]: 2354 : ABI_FREE(indkpt1)
2347 [ - + ]: 2354 : ABI_FREE(indsy1)
2348 [ - + ]: 2354 : ABI_FREE(istwfk_rbz)
2349 [ - + ]: 2354 : ABI_FREE(irrzon1)
2350 [ - + ]: 2354 : ABI_FREE(kg)
2351 [ - + ]: 2354 : ABI_FREE(kg1)
2352 [ - + ]: 2354 : ABI_FREE(kpq_rbz)
2353 [ + + ]: 2354 : if(.not.kramers_deg) then
2354 [ - + ]: 11 : ABI_FREE(kg1_mq)
2355 [ - + ]: 11 : ABI_FREE(kmq_rbz)
2356 : : end if
2357 [ - + ]: 2354 : ABI_FREE(kpt_rbz)
2358 [ - + ]: 2354 : ABI_FREE(nband_rbz)
2359 [ - + ]: 2354 : ABI_FREE(npwarr)
2360 [ - + ]: 2354 : ABI_FREE(npwar1)
2361 [ - + ]: 2354 : ABI_FREE(npwtot)
2362 [ - + ]: 2354 : ABI_FREE(npwtot1)
2363 [ - + ]: 2354 : ABI_FREE(occkq)
2364 [ - + ]: 2354 : ABI_FREE(occ_rbz)
2365 [ - + ]: 2354 : ABI_FREE(phnons1)
2366 [ - + ]: 2354 : ABI_FREE(resid)
2367 [ - + ]: 2354 : ABI_FREE(rhog1)
2368 [ - + ]: 2354 : ABI_FREE(rhor1)
2369 [ - + ]: 2354 : ABI_FREE(symaf1)
2370 [ - + ]: 2354 : ABI_FREE(symrc1)
2371 [ - + ]: 2354 : ABI_FREE(symrl1)
2372 [ - + ]: 2354 : ABI_FREE(tnons1)
2373 [ - + ]: 2354 : ABI_FREE(wtk_rbz)
2374 [ - + ]: 2354 : ABI_FREE(xccc3d1)
2375 [ - + ]: 2354 : ABI_FREE(vpsp1)
2376 [ - + ]: 2354 : ABI_FREE(ylm)
2377 [ - + ]: 2354 : ABI_FREE(ylm1)
2378 [ - + ]: 2354 : ABI_FREE(ylmgr)
2379 [ - + ]: 2354 : ABI_FREE(ylmgr1)
2380 [ + + ]: 2354 : 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 [ + + ]: 2354 : if (psps%usepaw==1) then
2393 : 409 : call pawang_free(pawang1)
2394 : 409 : call pawrhoij_free(pawrhoij1)
2395 [ + - ]: 409 : if (usecprj==1) then
2396 : 409 : call pawcprj_free(cprj)
2397 : 409 : call pawcprj_free(cprjq)
2398 : : end if
2399 : : end if
2400 [ + - + + : 3138 : ABI_FREE(pawrhoij1)
- + - + -
+ - + - +
- + - + ]
2401 [ + - + + : 161230 : ABI_FREE(cprjq)
- + - + -
+ ]
2402 [ + - + + : 161230 : ABI_FREE(cprj)
- + - + -
+ ]
2403 : : if(xmpi_paral==1) then
2404 [ - + ]: 2354 : ABI_FREE(mpi_enreg%proc_distrb)
2405 [ - + ]: 2354 : ABI_FREE(mpi_enreg%my_kpttab)
2406 : : end if
2407 : 2354 : call hdr%free()
2408 [ + + ]: 2354 : 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 : 2354 : call ebands_k%free()
2414 : 2354 : call ebands_kq%free()
2415 [ + + ]: 2354 : if(.not.kramers_deg) call ebands_kmq%free()
2416 : :
2417 : : ! %%%% Parallelization over perturbations %%%%%
2418 : : ! Redefine output/log files
2419 : 2354 : call localredirect(mpi_enreg%comm_cell,mpi_enreg%comm_world,npert_io,mpi_enreg%paral_pert,0)
2420 : :
2421 [ - + ]: 2354 : ABI_FREE(bz2ibz_smap)
2422 [ - + ]: 2354 : ABI_FREE(distrb_flags)
2423 : :
2424 : 2354 : call timab(146,2,tsec)
2425 [ + - ]: 10136 : if(iexit/=0) exit
2426 : : end do ! End loop on perturbations
2427 : :
2428 : : ABI_NVTX_END_RANGE()
2429 [ - + ]: 720 : ABI_FREE(zeff)
2430 : :
2431 : : !%%%% Parallelization over perturbations %%%%%
2432 : : !*Restore default communicators
2433 : 720 : call unset_pert_comm(mpi_enreg)
2434 : : ! Gather output/log files
2435 [ + + + - : 2160 : ABI_MALLOC(dyn,(npert_io))
+ - + + -
+ - + ]
2436 [ + + + + ]: 790 : if (npert_io>0) dyn=1
2437 : 720 : call localrdfile(mpi_enreg%comm_pert,mpi_enreg%comm_world,.true.,npert_io,mpi_enreg%paral_pert,0,dyn)
2438 [ - + ]: 720 : ABI_FREE(dyn)
2439 : :
2440 : : ! Restore PAW on-site data
2441 : 720 : 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 : 720 : call timab(147,1,tsec)
2454 : 720 : smdelta = dtset%smdelta
2455 : 720 : bdeigrf = dtset%bdeigrf
2456 [ + + ]: 720 : if(dtset%bdeigrf == -1) bdeigrf = dtset%mband
2457 : :
2458 [ + + ]: 720 : 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 [ + + ]: 720 : 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 [ + + ]: 720 : if(dtset%ieig2rf /= 3 .and. dtset%ieig2rf /= 4 .and. dtset%ieig2rf /= 5) call hdr0%free()
2552 [ - + ]: 720 : ABI_FREE(eigen0_copy)
2553 : 720 : call crystal%free()
2554 : :
2555 : 720 : call timab(147,2,tsec)
2556 : : !######################################################################################
2557 : :
2558 : : !Get ddk file information, for later use in dfpt_dyout
2559 : 720 : ddkfil(:)=0
2560 [ + + ]: 2880 : do idir=1,3
2561 : 2160 : file_index(1)=idir+dtset%natom*3
2562 : 2160 : call appdig(file_index(1),dtfil%fnamewffddk,fiwfddk)
2563 : : ! Check that ddk file exists
2564 : 2160 : t_exist = file_exists(fiwfddk)
2565 [ + - ]: 2160 : if (.not. t_exist) then
2566 : : ! Trick needed to run Abinit test suite in netcdf mode.
2567 : 2160 : t_exist = file_exists(nctk_ncify(fiwfddk))
2568 [ + + ]: 2160 : 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 : 1133 : if (t_exist) ddkfil(idir)=20+idir
2577 : : end do
2578 : :
2579 [ - + ]: 720 : ABI_FREE(ph1d)
2580 [ - + ]: 720 : ABI_FREE(ph1df)
2581 [ - + ]: 720 : ABI_FREE(pert_calc)
2582 [ + + ]: 720 : if (psps%usepaw==1) then
2583 [ - + ]: 116 : ABI_FREE(dimcprj_srt)
2584 : : end if
2585 : :
2586 : : !destroy dtset_tmp
2587 [ + + ]: 720 : 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 [ + + ]: 720 : 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 [ + + ]: 720 : if (associated(old_atmtab)) then
2616 : 36 : ABI_FREE(old_atmtab)
2617 : : nullify(old_atmtab)
2618 : : end if
2619 : :
2620 : 720 : call timab(141,2,tsec)
2621 : :
2622 : : DBG_EXIT("COLL")
2623 : :
2624 [ - + - + : 2880 : 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 : 2354 : 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 [ + + ]: 2354 : & berryopt==14.or.berryopt==16.or.berryopt==17)
2882 : 2354 : if (ipert==natom+1) nn=8
2883 [ + + + + ]: 2354 : if (ipert==natom+5.or.ipert==natom+6) nn=7
2884 [ - + ]: 2354 : if (ipert==natom+7) nn=8
2885 [ + + ]: 2354 : if (ipert==natom+2) nn=7
2886 [ - + + + ]: 2354 : if (ipert>=1.and.ipert<=natom) nn=13
2887 [ + + + + ]: 2354 : if (ipert==natom+3.or.ipert==natom+4) nn=17
2888 [ + + + + ]: 2354 : if (ipert==natom+2.and.berry_activated) nn=nn+1
2889 [ + + + + ]: 2354 : if (ipert==natom+10.or.ipert==natom+11) nn=1 ! means nothing,
2890 [ + + - + ]: 2354 : if (ipert>natom+11.and.ipert<=2*natom+11) nn=7
2891 [ + + ]: 2354 : if (emagpen1>tol8) nn=nn+1
2892 : : ! because we do not compute derivatives of the energy in this case
2893 [ + + ]: 2354 : if (usepaw==1) nn=nn+1
2894 [ + + ]: 2354 : if (usevdw==1) nn=nn+1
2895 : 2354 : write(msg, '(4a)' ) ch10,&
2896 [ + - ]: 4708 : & ' ',trim(numbstr(nn)),' components of 2nd-order total energy (hartree) are '
2897 : 2354 : call wrtout(iout,msg)
2898 : 2354 : call wrtout(std_out,msg)
2899 : :
2900 : 2354 : numb='1,2,3'
2901 : : write(msg, '(3a)' )&
2902 : 2354 : & ' ',trim(numb),': 0th-order hamiltonian combined with 1st-order wavefunctions'
2903 : 2354 : call wrtout(iout,msg)
2904 : 2354 : call wrtout(std_out,msg)
2905 : : write(msg, '(a,es17.8,a,es17.8,a,es17.8)' )&
2906 : 2354 : & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0
2907 [ + + ]: 2354 : if (has_nd) then
2908 : : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' )&
2909 : 60 : & ' kin0=',ek0, ' eigvalue=',eeig0,' local=',eloc0,' nclr dpl0=',end0
2910 : : end if
2911 [ + + ]: 2354 : 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 : 2354 : call wrtout(iout,msg)
2916 : 2354 : call wrtout(std_out,msg)
2917 : :
2918 [ + + + + ]: 2354 : numb='4,5,6';if( ipert==natom+3.or.ipert==natom+4) numb='4,5,6,7'
2919 : : write(msg, '(3a)' )&
2920 : 2354 : & ' ',trim(numb),': 1st-order hamiltonian combined with 1st and 0th-order wfs'
2921 : 2354 : call wrtout(iout,msg)
2922 : 2354 : call wrtout(std_out,msg)
2923 [ + + + + : 2354 : 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 [ + + ]: 702 : else if(ipert==natom+1) then
2928 : : write(msg, '(a,es17.8,a,es17.8,a,es17.8)' ) &
2929 : 448 : & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1
2930 [ + + + + ]: 448 : if (has_nd .and. (.not.with_vxctau) ) then
2931 : : write(msg, '(a,es17.8,a,es17.8,a,es17.8,a,es17.8)' ) &
2932 : 54 : & ' kin1=',ek1, ' Hartree=',ehart1,' xc=',exc1,' nclr dpl1=',end1
2933 : : end if
2934 [ + + - + ]: 448 : 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 [ + + + - ]: 448 : 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 [ + + + + ]: 2354 : 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 : 2354 : call wrtout(iout,msg)
2956 : 2354 : call wrtout(std_out,msg)
2957 : :
2958 [ + + ]: 2354 : enl1_effective=enl1;if (ipert==natom+2) enl1_effective=zero
2959 [ + + + + ]: 2354 : 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 : 2354 : & ' ',trim(numb),': eventually, occupation + non-local contributions',ch10,&
2962 : 4708 : & ' edocc=',edocc,' enl0=',enl0,' enl1=',enl1_effective
2963 : 2354 : call wrtout(iout,msg)
2964 : 2354 : call wrtout(std_out,msg)
2965 : :
2966 [ + + ]: 2354 : 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 [ + + ]: 2354 : if (usepaw==1) then
2974 [ + + + + ]: 409 : numb='10';if( ipert==natom+3.or.ipert==natom+4) numb='11'
2975 : : write(msg, '(3a,es17.8)' )&
2976 : 409 : & ' ',trim(numb),': eventually, PAW "on-site" Hxc contribution: epaw1=',epaw1
2977 : 409 : call wrtout(iout,msg)
2978 : 409 : call wrtout(std_out,msg)
2979 : : end if
2980 : :
2981 [ + + + + ]: 2354 : if(ipert/=natom+10 .and.ipert/=natom+11) then
2982 : 2175 : erelax=0.0_dp
2983 [ + - + + ]: 2175 : if(ipert>=1.and.ipert<=natom)then
2984 : 1353 : erelax=ek0+edocc+eeig0+eloc0+elpsp1+ehart1+exc1+enl0+enl1+epaw1
2985 [ + + + + ]: 822 : 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 : 642 : 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 [ + + ]: 2175 : if (abs(emagpen1)>tol8) erelax=erelax+emagpen1
2999 : 2175 : enl1_effective=enl1
3000 [ + + + + ]: 2175 : if (ipert==natom+1.or.ipert==natom+2) then
3001 [ + + ]: 642 : if (1.0_dp+enl1/10.0_dp==1.0_dp) enl1_effective=zero
3002 : : end if
3003 : :
3004 [ + + ]: 2175 : numb='1-9';if (usepaw==1) numb='1-10'
3005 [ + + + + ]: 2175 : 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 : 2175 : & ' ',trim(numb),' gives the relaxation energy (to be shifted if some occ is /=2.0)',&
3010 : 4350 : & ch10,' erelax=',erelax
3011 : 2175 : call wrtout(iout,msg)
3012 : 2175 : call wrtout(std_out,msg)
3013 : : end if
3014 : :
3015 [ + - + + ]: 2354 : 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 [ + + + + ]: 1001 : else if(ipert==natom+1.or.ipert==natom+2)then
3075 [ + + + + ]: 642 : if (ipert==natom+1.and.usepaw==1) then
3076 : 141 : numb='11'
3077 : : write(msg, '(5a,es17.8)' )&
3078 : 141 : & ' ',trim(numb),' Contribution from 1st-order change of wavefunctions overlap',&
3079 : 282 : & ch10,' eovl1 =',eovl1
3080 : 141 : call wrtout(iout,msg)
3081 : 141 : call wrtout(std_out,msg)
3082 : : end if
3083 : 642 : write(msg,*)' No Ewald or frozen-wf contrib.:',' the relaxation energy is the total one'
3084 [ + + ]: 642 : 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 : 642 : call wrtout(iout,msg)
3089 : 642 : call wrtout(std_out,msg)
3090 : 642 : etotal=erelax
3091 : : write(msg, '(a,e20.10,a,e22.12,a)' ) &
3092 : 642 : & ' 2DEtotal=',etotal,' Ha. Also 2DEtotal=',etotal*Ha_eV,' eV'
3093 : 642 : call wrtout(iout,msg)
3094 : 642 : call wrtout(std_out,msg)
3095 : : write(msg, '(a,es20.10,a)' ) &
3096 : 642 : & ' ( non-var. 2DEtotal :',0.5_dp*(ek1+enl1_effective+end1+evxctau1)+eovl1,' Ha)'
3097 : 642 : call wrtout(iout,msg)
3098 : 642 : 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 : 2354 : 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 : 545 : 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 [ - + ]: 545 : 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 [ + + ]: 1122 : do isppol=1,nsppol
3219 [ + + ]: 25066 : do ikpt=1,nkpt
3220 : 23944 : nband_k=nband(ikpt+(isppol-1)*nkpt)
3221 [ + + ]: 23944 : if(option==1)then
3222 [ + + ]: 170536 : do iband=1,nband_k
3223 : : eigenresp_mean(iband+bdtot_index)=&
3224 : 170536 : & 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 [ + + ]: 23944 : if(nband_k>1)then
3234 : 23914 : eig0=eigen0(1+bdtot_index)
3235 : 23914 : ii=1
3236 [ + + ]: 164657 : do iband=2,nband_k
3237 [ + + ]: 140743 : if(eigen0(iband+bdtot_index)-eig0<tol8)then
3238 : 8124 : ii=ii+1
3239 : : else
3240 [ + + ]: 271357 : mean=sum(eigenresp_mean(iband-ii+bdtot_index:iband-1+bdtot_index))/ii
3241 [ + + ]: 271357 : eigenresp_mean(iband-ii+bdtot_index:iband-1+bdtot_index)=mean
3242 : : ii=1
3243 : : end if
3244 : 140743 : eig0=eigen0(iband+bdtot_index)
3245 [ + + ]: 164657 : if(iband==nband_k)then
3246 [ + + ]: 49833 : mean=sum(eigenresp_mean(iband-ii+1+bdtot_index:iband+bdtot_index))/ii
3247 [ + + ]: 49833 : eigenresp_mean(iband-ii+1+bdtot_index:iband+bdtot_index)=mean
3248 : : end if
3249 : : end do
3250 : : end if
3251 : :
3252 : 23944 : bdtot_index=bdtot_index+nband_k
3253 : 24521 : bd2tot_index=bd2tot_index+2*nband_k**2
3254 : : end do
3255 : : end do
3256 : :
3257 : 545 : 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 : : !!***
|
