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