Line data Source code
1 : !!****m* m_paw_dfpt/m_paw_dfpt
2 : !! NAME
3 : !! m_paw_dfpt
4 : !!
5 : !! FUNCTION
6 : !! This module contains several routines related to the 1st and 2nd order derivatives
7 : !! (in the DFPT approach) of PAW on-site quantities.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2018-2026 ABINIT group (MT,AM,FJ,JWZ)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_paw_dfpt
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_xmpi
28 : use m_errors
29 : use m_time, only : timab
30 :
31 : use defs_datatypes, only : pseudopotential_type
32 : use m_pawang, only : pawang_type
33 : use m_pawrad, only : pawrad_type
34 : use m_pawtab, only : pawtab_type
35 : use m_paw_an, only : paw_an_type
36 : use m_paw_ij, only : paw_ij_type
37 : use m_pawcprj, only : pawcprj_type
38 : use m_pawdij, only : pawdijhartree,pawdiju_euijkl
39 : use m_pawrhoij, only : pawrhoij_type,pawrhoij_free,pawrhoij_gather,pawrhoij_nullify
40 : use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_nullify, pawfgrtab_gather
41 : use m_paw_finegrid, only : pawgylm, pawrfgd_fft, pawexpiqr
42 : use m_pawxc, only : pawxc_dfpt, pawxcm_dfpt
43 : use m_paw_denpot, only : pawdensities,pawaccenergy,pawaccenergy_nospin
44 : use m_paral_atom, only : get_my_atmtab,free_my_atmtab
45 : use m_atm2fft, only : dfpt_atm2fft
46 : use m_distribfft, only : distribfft_type
47 : use m_geometry, only : metric, stresssym
48 : use m_efield, only : efield_type
49 :
50 : implicit none
51 :
52 : private
53 :
54 : !public procedures.
55 : public :: pawdfptenergy ! Compute Hartree+XC PAW on-site contrib. to a 1st or 2nd-order energy
56 : public :: pawgrnl ! Compute derivatives of total energy due to NL terms (PAW Dij derivatives)
57 : public :: dsdr_k_paw ! Compute PAW on-site terms for forces/stresses for finite electric fields
58 :
59 : CONTAINS !========================================================================================
60 : !!***
61 :
62 : !----------------------------------------------------------------------
63 :
64 : !!****f* m_paw_dfpt/pawdfptenergy
65 : !! NAME
66 : !! pawdfptenergy
67 : !!
68 : !! FUNCTION
69 : !! This routine compute the Hartree+XC+U PAW on-site contributions to a 1st-order or 2nd-order energy.
70 : !! These contributions are equal to:
71 : !! E_onsite=
72 : !! Int{ VHxc[n1_a^(j1);nc^(j1)].n1_b^(j2) }
73 : !! -Int{ VHxc[tild_n1_a^(j1)+hat_n1_a^(j1);tild_n_c^(j1)].(tild_n1_b+n1_b)^(j2) }
74 : !! Some typical uses:
75 : !! A-Contribution to non-stationary expression of the 2nd-order total energy:
76 : !! In that case, n1_a^(1)[r]=n1^(j1)[r] and n1_b[r]=delta_n1^(j2)[r]
77 : !! where j1 and j2 are two given perturbations,
78 : !! and delta_n1^(j)[r] is the 1s-order density only due to change of WF overlap.
79 : !! See PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq.(80).
80 : !! E_onsite=
81 : !! Int{ VHxc[n1^(j1);nc^(j1)].delta_n1^(j2) }
82 : !! -Int{ VHxc[tild_n1^(j1)+hat_n1^(j1);tild_n_c^(j1)].delta_(tild_n1+hat_n1)^(j2) }
83 : !! B-Contribution to first-order Fermi energy:
84 : !! In that case, n1_a^(1)[r]=n1^(j1)[r] and n1_b[r]=n1[r,EFermi]
85 : !! where j1 is the current perturbation, and n1[r,EFermi] is the density at Fermi level.
86 : !! E_onsite=
87 : !! Int{ VHxc[n1^(j1);nc^(j1)].n1[r,EFermi] }
88 : !! -Int{ VHxc[tild_n1^(j1)+hat_n1^(j1);tild_n_c^(j1)].(tild_n1+hat_n1)[r,EFermi] }
89 : !!
90 : !! INPUTS
91 : !! ipert1,ipert2=indexes of perturbations (j1) and (j2)
92 : !! if ipert2<=0, we compute a first-order energy
93 : !! if ipert2> 0, we compute a second-order energy
94 : !! ixc= choice of exchange-correlation scheme
95 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
96 : !! comm_atom=--optional-- MPI communicator over atoms
97 : !! my_natom=number of atoms treated by current processor
98 : !! natom=total number of atoms in cell
99 : !! ntypat=number of types of atoms in unit cell.
100 : !! nzlmopt_a= For the n1_a density:
101 : !! if -1, compute all LM-moments of the density and use non-zero LM-moments
102 : !! if 0, compute all LM-moments of the density and use all LM-moments
103 : !! if +1, compute only non-zero LM-moments of the density (stored before)
104 : !! nzlmopt_b= For the n1_b density:
105 : !! if -1, compute all LM-moments of the density and use non-zero LM-moments
106 : !! if 0, compute all LM-moments of the density and use all LM-moments
107 : !! if +1, compute only non-zero LM-moments of the density (stored before)
108 : !! paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
109 : !! paw_an1(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh
110 : !! This corresponds to (j1) perturbation
111 : !! paw_ij1(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
112 : !! This corresponds to (j1) perturbation
113 : !! pawang <type(pawang_type)>=paw angular mesh and related data
114 : !! pawprtvol=control print volume and debugging output for PAW
115 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
116 : !! pawrhoij_a(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies for the (j1) perturbation
117 : !! pawrhoij_b(natom) <type(pawrhoij_type)>=
118 : !! if ipert2> 0: paw rhoij 1st-order occupancies for the (j2) perturbation corrsponding to n1_b^(j2)[r]
119 : !! if ipert2<=0: paw rhoij occupancies corresponding to n1_b[r]
120 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
121 : !! pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
122 : !! xclevel= XC functional level
123 : !!
124 : !! OUTPUT
125 : !! delta_energy(2)= real and imaginary parts of contributions to non-stationary expression for the
126 : !! second derivative of the total energy
127 : !!
128 : !! SIDE EFFECTS
129 : !! ==== if paw_an1(:)%has_vxc<2, compute 1st-order XC potentials
130 : !! paw_an1(natom)%vxc1(cplex_a*mesh_size,:,nspden) =AE 1st-order XC potential Vxc^(j1)
131 : !! paw_an1(natom)%vxct1(cplex_a*mesh_size,:,nspden)=PS 1st-order XC potential tVxc^(j1)
132 : !! ==== if paw_ij1(:)%has_dijhartree<2, compute 1st-order Dij_hartree
133 : !! paw_ij1(natom)%dijhartree(cplex_a*lmn2_size)=Hartree contribution to Dij^(j1)
134 : !! ==== if paw_ij1(:)%has_dijU<2, compute 1st-order Dij_U
135 : !! paw_ij1(natom)%diju(cplex_a*lmn2_size)=DFT+U contribution to Dij^(j1)
136 : !!
137 : !! SOURCE
138 :
139 4196 : subroutine pawdfptenergy(delta_energy,ipert1,ipert2,ixc,my_natom,natom,ntypat,nzlmopt_a,nzlmopt_b,&
140 8392 : & paw_an0,paw_an1,paw_ij1,pawang,pawprtvol,pawrad,pawrhoij_a,pawrhoij_b,&
141 4196 : & pawtab,pawxcdev,xclevel, &
142 4196 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
143 :
144 : !Arguments ---------------------------------------------
145 : !scalars
146 : integer,intent(in) :: ipert1,ipert2,ixc,my_natom,natom,ntypat,nzlmopt_a,nzlmopt_b
147 : integer,intent(in) :: pawprtvol,pawxcdev,xclevel
148 : integer,optional,intent(in) :: comm_atom
149 : type(pawang_type),intent(in) :: pawang
150 : !arrays
151 : integer,optional,target,intent(in) :: mpi_atmtab(:)
152 : real(dp),intent(out) :: delta_energy(2)
153 : type(paw_an_type),intent(in) :: paw_an0(my_natom)
154 : type(paw_an_type),intent(inout) :: paw_an1(my_natom)
155 : type(paw_ij_type),intent(inout) :: paw_ij1(my_natom)
156 : type(pawrad_type),intent(in) :: pawrad(ntypat)
157 : type(pawrhoij_type),intent(in) :: pawrhoij_a(my_natom),pawrhoij_b(my_natom)
158 : type(pawtab_type),intent(in) :: pawtab(ntypat)
159 :
160 : !Local variables ---------------------------------------
161 : !scalars
162 : integer, parameter :: PAWU_ALGO_1=1,PAWU_ALGO_2=2
163 : integer :: cplex_a,cplex_b,cplex_vxc1,iatom,iatom_tot,ierr,itypat,lm_size_a,lm_size_b,mesh_size
164 : integer :: my_comm_atom,nspden,opt_compch,optexc,optvxc,pawu_algo,qphase_dijh1,qphase_diju1
165 : integer :: usecore,usepawu,usetcore,usexcnhat
166 : logical :: my_atmtab_allocated,non_magnetic_xc,paral_atom
167 : real(dp) :: compch,eexc,eexc_im
168 : character(len=500) :: msg
169 : !arrays
170 4196 : integer,pointer :: my_atmtab(:)
171 4196 : logical,allocatable :: lmselect_a(:),lmselect_b(:),lmselect_tmp(:)
172 : real(dp) :: delta_energy_h(2),delta_energy_u(2),delta_energy_xc(2),tsec(2)
173 4196 : real(dp),allocatable :: kxc_dum(:,:,:),nhat1(:,:,:),rho1(:,:,:),trho1(:,:,:)
174 :
175 : ! *************************************************************************
176 :
177 : DBG_ENTER("COLL")
178 :
179 4196 : call timab(567,1,tsec)
180 :
181 4196 : if (.not.(ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11 &
182 : & .or.ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11)) then
183 4196 : if((abs(nzlmopt_a)/=1.and.nzlmopt_a/=0).or.(abs(nzlmopt_b)/=1.and.nzlmopt_b/=0)) then
184 0 : msg='invalid value for nzlmopt!'
185 0 : ABI_BUG(msg)
186 : end if
187 4196 : if (my_natom>0) then
188 4172 : if(paw_ij1(1)%has_dijhartree==0) then
189 0 : ABI_BUG('dijhartree must be allocated!')
190 : end if
191 12251 : if (any(pawtab(1:ntypat)%usepawu/=0)) then
192 0 : if(paw_ij1(1)%has_dijU==0) then
193 0 : ABI_BUG('dijU must be allocated!')
194 : end if
195 : end if
196 4172 : if(paw_an1(1)%has_vxc==0) then
197 0 : ABI_BUG('vxc1 and vxct1 must be allocated!')
198 : end if
199 4172 : if(paw_an0(1)%has_kxc==0) then
200 0 : ABI_BUG('kxc1 must be allocated!')
201 : end if
202 4172 : if ((ipert1<=natom.or.ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11).and.paw_an0(1)%has_kxc/=2) then
203 0 : ABI_BUG('XC kernels for ground state must be in memory!')
204 : end if
205 4172 : if (paw_ij1(1)%qphase/=paw_an1(1)%cplex) then
206 0 : msg='paw_ij1()%qphase and paw_an1()%cplex must be equal!'
207 0 : ABI_BUG(msg)
208 : end if
209 4172 : if (pawrhoij_a(1)%qphase<paw_an1(1)%cplex.or.pawrhoij_b(1)%qphase<paw_an1(1)%cplex) then
210 0 : msg='pawrhoij()%qphase must be >=paw_an1()%cplex!'
211 0 : ABI_BUG(msg)
212 : end if
213 4172 : if (paw_ij1(1)%nspden/=paw_an1(1)%nspden) then
214 0 : msg='paw_ij1()%nspden and paw_an1()%nspden must be equal!'
215 0 : ABI_BUG(msg)
216 : end if
217 4172 : if (pawrhoij_a(1)%nspden/=pawrhoij_b(1)%nspden) then
218 0 : msg='pawrhoij_a()%nspden must =pawrhoij_b()%nspden !'
219 0 : ABI_BUG(msg)
220 : end if
221 : end if
222 : end if
223 :
224 : !Set up parallelism over atoms
225 4196 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
226 4196 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
227 4196 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
228 4196 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
229 :
230 : !Init contribution to 1st-order (or 2nd-order) energy
231 4196 : delta_energy(1:2)=zero
232 :
233 : !For some perturbations, nothing else to do
234 : if (ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11 .or. &
235 4196 : & ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11) return
236 :
237 : !Various inits
238 4196 : opt_compch=0;optvxc=1;optexc=3
239 4196 : usecore=0;usetcore=0 ! This is true for phonons and Efield pert.
240 12323 : usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
241 4196 : delta_energy_xc(1:2)=zero
242 4196 : delta_energy_h(1:2)=zero
243 4196 : delta_energy_u(1:2)=zero
244 :
245 : !================ Loop on atomic sites =======================
246 13304 : do iatom=1,my_natom
247 9108 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
248 :
249 9108 : itypat=pawrhoij_a(iatom)%itypat
250 9108 : mesh_size=pawtab(itypat)%mesh_size
251 9108 : nspden=paw_an1(iatom)%nspden
252 9108 : cplex_a=pawrhoij_a(iatom)%qphase
253 9108 : cplex_b=pawrhoij_b(iatom)%qphase
254 9108 : cplex_vxc1=paw_an1(iatom)%cplex
255 9108 : qphase_dijh1=paw_ij1(iatom)%qphase
256 9108 : qphase_diju1=paw_ij1(iatom)%qphase
257 9108 : lm_size_a=paw_an1(iatom)%lm_size
258 9108 : if (ipert2<=0) lm_size_b=paw_an0(iatom)%lm_size
259 9108 : if (ipert2> 0) lm_size_b=paw_an1(iatom)%lm_size
260 9108 : usepawu=pawtab(itypat)%usepawu
261 9108 : pawu_algo=merge(PAWU_ALGO_1,PAWU_ALGO_2,ipert1<=0.and.ipert2<=0.and.usepawu>=0)
262 9108 : non_magnetic_xc=(mod(abs(usepawu),10)==4)
263 :
264 : ! If Vxc potentials are not in memory, compute them
265 9108 : if (paw_an1(iatom)%has_vxc/=2) then
266 2175 : ABI_MALLOC(rho1 ,(cplex_a*mesh_size,lm_size_a,nspden))
267 1740 : ABI_MALLOC(trho1,(cplex_a*mesh_size,lm_size_a,nspden))
268 2175 : ABI_MALLOC(nhat1,(cplex_a*mesh_size,lm_size_a,nspden*usexcnhat))
269 1305 : ABI_MALLOC(lmselect_a,(lm_size_a))
270 4526 : lmselect_a(:)=paw_an1(iatom)%lmselect(:)
271 870 : ABI_MALLOC(lmselect_tmp,(lm_size_a))
272 4526 : lmselect_tmp(:)=.true.
273 4526 : if (nzlmopt_a==1) lmselect_tmp(:)=lmselect_a(:)
274 : ! Compute on-site 1st-order densities
275 : call pawdensities(compch,cplex_a,iatom_tot,lmselect_tmp,lmselect_a,&
276 : & lm_size_a,nhat1,nspden,nzlmopt_a,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
277 435 : & pawrad(itypat),pawrhoij_a(iatom),pawtab(itypat),rho1,trho1)
278 435 : ABI_FREE(lmselect_tmp)
279 : ! Compute on-site 1st-order xc potentials
280 435 : if (pawxcdev/=0) then
281 : call pawxcm_dfpt(pawtab(itypat)%coredens,cplex_a,cplex_vxc1,eexc,ixc,paw_an0(iatom)%kxc1,&
282 : & lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,optvxc,&
283 : & pawang,pawrad(itypat),rho1,usecore,0,&
284 371 : & paw_an1(iatom)%vxc1,xclevel)
285 : call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),cplex_a,cplex_vxc1,eexc,ixc,paw_an0(iatom)%kxct1,&
286 : & lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,optvxc,&
287 : & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,&
288 371 : & paw_an1(iatom)%vxct1,xclevel)
289 : else
290 : call pawxc_dfpt(pawtab(itypat)%coredens,cplex_a,cplex_vxc1,eexc,ixc,paw_an0(iatom)%kxc1,&
291 : & lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,optvxc,&
292 : & pawang,pawrad(itypat),rho1,usecore,0,&
293 64 : & paw_an0(iatom)%vxc1,paw_an1(iatom)%vxc1,xclevel)
294 : call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),cplex_a,cplex_vxc1,eexc,ixc,paw_an0(iatom)%kxct1,&
295 : & lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,non_magnetic_xc,mesh_size,nspden,optvxc,&
296 : & pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,&
297 64 : & paw_an0(iatom)%vxct1,paw_an1(iatom)%vxct1,xclevel)
298 : end if
299 :
300 435 : paw_an1(iatom)%has_vxc=2
301 435 : ABI_FREE(lmselect_a)
302 435 : ABI_FREE(rho1)
303 435 : ABI_FREE(trho1)
304 435 : ABI_FREE(nhat1)
305 : end if ! has_vxc
306 :
307 : ! Compute contribution to 1st-order (or 2nd-order) energy from 1st-order XC potential
308 45540 : ABI_MALLOC(rho1 ,(cplex_b*mesh_size,lm_size_b,nspden))
309 36432 : ABI_MALLOC(trho1,(cplex_b*mesh_size,lm_size_b,nspden))
310 45540 : ABI_MALLOC(nhat1,(cplex_b*mesh_size,lm_size_b,nspden*usexcnhat))
311 27324 : ABI_MALLOC(lmselect_b,(lm_size_b))
312 34569 : if (ipert2<=0) lmselect_b(:)=paw_an0(iatom)%lmselect(:)
313 66579 : if (ipert2> 0) lmselect_b(:)=paw_an1(iatom)%lmselect(:)
314 18216 : ABI_MALLOC(lmselect_tmp,(lm_size_b))
315 92040 : lmselect_tmp(:)=.true.
316 9108 : if (nzlmopt_b==1) lmselect_tmp(:)=lmselect_b(:)
317 : ! Compute on-site 1st-order densities
318 : call pawdensities(compch,cplex_b,iatom_tot,lmselect_tmp,lmselect_b,&
319 : & lm_size_b,nhat1,nspden,nzlmopt_b,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
320 9108 : & pawrad(itypat),pawrhoij_b(iatom),pawtab(itypat),rho1,trho1)
321 9108 : ABI_FREE(lmselect_tmp)
322 : ! Compute contributions to 1st-order (or 2nd-order) energy
323 9108 : if (pawxcdev/=0) then
324 25128 : ABI_MALLOC(kxc_dum,(mesh_size,pawang%angl_size,0))
325 : call pawxcm_dfpt(pawtab(itypat)%coredens,cplex_b,cplex_vxc1,eexc,ixc,kxc_dum,&
326 : & lm_size_b,lmselect_b,nhat1,0,non_magnetic_xc,mesh_size,nspden,optexc,pawang,pawrad(itypat),&
327 8376 : & rho1,usecore,0,paw_an1(iatom)%vxc1,xclevel,d2enxc_im=eexc_im)
328 8376 : delta_energy_xc(1)=delta_energy_xc(1)+eexc
329 8376 : delta_energy_xc(2)=delta_energy_xc(2)+eexc_im
330 : call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),&
331 : & cplex_b,cplex_vxc1,eexc,ixc,kxc_dum,&
332 : & lm_size_b,lmselect_b,nhat1,0,non_magnetic_xc,mesh_size,nspden,optexc,pawang,pawrad(itypat),&
333 : & trho1,usetcore,2*usexcnhat,paw_an1(iatom)%vxct1,xclevel,&
334 8376 : & d2enxc_im=eexc_im)
335 8376 : ABI_FREE(kxc_dum)
336 8376 : delta_energy_xc(1)=delta_energy_xc(1)-eexc
337 8376 : delta_energy_xc(2)=delta_energy_xc(2)-eexc_im
338 : else
339 2196 : ABI_MALLOC(kxc_dum,(mesh_size,lm_size_b,0))
340 : call pawxc_dfpt(pawtab(itypat)%coredens,cplex_b,cplex_vxc1,eexc,ixc,kxc_dum,&
341 : & lm_size_b,lmselect_b,nhat1,0,non_magnetic_xc,mesh_size,nspden,optexc,pawang,pawrad(itypat),&
342 732 : & rho1,usecore,0,paw_an0(iatom)%vxc1,paw_an1(iatom)%vxc1,xclevel,d2enxc_im=eexc_im)
343 732 : delta_energy_xc(1)=delta_energy_xc(1)+eexc
344 732 : delta_energy_xc(2)=delta_energy_xc(2)+eexc_im
345 : call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),&
346 : & cplex_b,cplex_vxc1,eexc,ixc,kxc_dum,&
347 : & lm_size_b,lmselect_b,nhat1,0,non_magnetic_xc,mesh_size,nspden,optexc,pawang,pawrad(itypat),&
348 : & trho1,usetcore,2*usexcnhat,paw_an0(iatom)%vxct1,paw_an1(iatom)%vxct1,xclevel,&
349 732 : & d2enxc_im=eexc_im)
350 732 : ABI_FREE(kxc_dum)
351 732 : delta_energy_xc(1)=delta_energy_xc(1)-eexc
352 732 : delta_energy_xc(2)=delta_energy_xc(2)-eexc_im
353 : end if
354 9108 : ABI_FREE(lmselect_b)
355 9108 : ABI_FREE(rho1)
356 9108 : ABI_FREE(trho1)
357 9108 : ABI_FREE(nhat1)
358 :
359 : ! If Dij_hartree are not in memory, compute them
360 9108 : if (paw_ij1(iatom)%has_dijhartree/=2) then
361 : call pawdijhartree(paw_ij1(iatom)%dijhartree,qphase_dijh1,paw_ij1(iatom)%nspden,&
362 435 : & pawrhoij_a(iatom),pawtab(itypat))
363 435 : paw_ij1(iatom)%has_dijhartree=2
364 : end if
365 :
366 : ! Compute contribution to 1st-order(or 2nd-order) energy from 1st-order Hartree potential
367 : call pawaccenergy_nospin(delta_energy_h(1),pawrhoij_b(iatom),paw_ij1(iatom)%dijhartree, &
368 9108 : & 1,qphase_dijh1,pawtab(itypat),epaw_im=delta_energy_h(2))
369 :
370 : ! Compute contribution to 1st-order(or 2nd-order) energy from 1st-order PAW+U potential
371 13304 : if (usepawu/=0.and.pawu_algo==PAWU_ALGO_2) then
372 : ! If DijU are not in memory, compute them
373 0 : if (paw_ij1(iatom)%has_dijU/=2) then ! We force the recomputation of dijU in when cplex=2 to get diju_im
374 : call pawdiju_euijkl(paw_ij1(iatom)%dijU,paw_ij1(iatom)%cplex_dij,qphase_diju1,&
375 0 : & paw_ij1(iatom)%ndij,pawrhoij_a(iatom),pawtab(itypat))
376 0 : paw_ij1(iatom)%has_dijU=2
377 : end if
378 : ! Compute contribution to 1st-order(or 2nd-order) energy
379 : call pawaccenergy(delta_energy_u(1),pawrhoij_b(iatom),paw_ij1(iatom)%dijU,paw_ij1(iatom)%cplex_dij, &
380 0 : & qphase_diju1,paw_ij1(iatom)%ndij,pawtab(itypat),epaw_im=delta_energy_u(2))
381 : ! Add FLL double-counting contribution
382 0 : if (ipert1==0) then ! If j1/=0, Dij^FLL^(j1)=0 because it is constant
383 : call pawaccenergy_nospin(delta_energy_u(1),pawrhoij_b(iatom),pawtab(itypat)%euij_fll,1,1,&
384 0 : & pawtab(itypat),epaw_im=delta_energy_u(2))
385 : end if
386 : end if
387 :
388 : ! ================ End loop on atomic sites =======================
389 : end do
390 :
391 : !Final building of 1st-order (or 2nd-order) energy
392 12588 : delta_energy(1:2)=delta_energy_xc(1:2)+delta_energy_h(1:2)+delta_energy_u(1:2)
393 :
394 : !Reduction in case of parallelism
395 4196 : if (paral_atom) then
396 2340 : call xmpi_sum(delta_energy,my_comm_atom,ierr)
397 : end if
398 :
399 : !Destroy atom table used for parallelism
400 4196 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
401 :
402 4196 : call timab(567,2,tsec)
403 :
404 : DBG_EXIT("COLL")
405 :
406 4196 : end subroutine pawdfptenergy
407 : !!***
408 :
409 : !----------------------------------------------------------------------
410 :
411 : !!****f* m_paw_dfpt/pawgrnl
412 : !!
413 : !! NAME
414 : !! pawgrnl
415 : !!
416 : !! FUNCTION
417 : !! PAW: Add to GRadients of total energy due to non-local term of Hamiltonian
418 : !! the contribution due to Dij derivatives
419 : !! In particular, compute contribution to forces, stresses, dyn. matrix
420 : !! Remember: Vnl=Sum_ij[|p_i>Dij<p_j|]
421 : !!
422 : !! INPUTS
423 : !! atindx1(natom)=index table for atoms, inverse of atindx
424 : !! dimnhat=second dimension of array nhat (0 or # of spin components)
425 : !! distribfft<type(distribfft_type)>=--optional-- contains all the information related
426 : !! to the FFT parallelism and plane sharing
427 : !! dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
428 : !! gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere
429 : !! mgfft=maximum size of 1D FFTs
430 : !! me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
431 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
432 : !! comm_atom=--optional-- MPI communicator over atoms
433 : !! comm_fft=--optional-- MPI communicator over FFT components (=mpi_comm_grid is not present)
434 : !! mpi_comm_grid=--optional-- MPI communicator over real space grid components (=comm_fft is not present)
435 : !! my_natom=number of atoms treated by current processor
436 : !! natom=total number of atoms in cell
437 : !! nattyp(ntypat)=array describing how many atoms of each type in cell
438 : !! nfft=(effective) number of FFT grid points (for this processor)
439 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
440 : !! nhat(nfft,dimnhat)=compensation charge density on rectangular grid in real space
441 : !! nspden=number of spin-density components
442 : !! nsym=number of symmetries in space group
443 : !! ntypat=number of types of atoms
444 : !! optgr= 1 if gradients with respect to atomic position(s) have to be computed
445 : !! optgr2= 1 if 2nd gradients with respect to atomic position(s) have to be computed
446 : !! optstr= 1 if gradients with respect to strain(s) have to be computed
447 : !! optstr2= 1 if 2nd gradients with respect to strain(s) have to be computed
448 : !! paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
449 : !! pawang <type(pawang_type)>=paw angular mesh and related data
450 : !! pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
451 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
452 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
453 : !! ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information
454 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
455 : !! qphon(3)=wavevector of the phonon
456 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
457 : !! symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
458 : !! typat(natom)=types of atoms
459 : !! ucvol=unit cell volume
460 : !! vtrial(nfft,nspden)= total local potential
461 : !! vxc(nfft,nspden)=XC potential
462 : !! xred(3,natom)=reduced dimensionless atomic coordinates
463 : !!
464 : !! SIDE EFFECTS
465 : !! At input, this terms contain contribution from non-local projectors derivatives
466 : !! At output, they are updated with the contribution of Dij derivatives
467 : !! ==== if optgr=1 ====
468 : !! grnl(3*natom) =gradients of NL energy wrt atomic coordinates
469 : !! ==== if optstr=1 ====
470 : !! nlstr(6) =gradients of NL energy wrt strains
471 : !! ==== if optgr2=1 ====
472 : !! dyfrnl(dyfr_cplex,3,3,natom,natom) =2nd gradients of NL energy wrt atomic coordinates
473 : !! ==== if optstr=2 ====
474 : !! eltfrnl(6+3*natom,6)=non-symmetrized non-local contribution to the elastic tensor
475 : !! NOTES
476 : !! In the case of parallelisation over atoms and calculation of dynamical matrix (optgr2=1)
477 : !! several data are gathered and no more distributed inside this routine.
478 : !!
479 : !! SOURCE
480 :
481 14846 : subroutine pawgrnl(atindx1,dimnhat,dyfrnl,dyfr_cplex,eltfrnl,grnl,gsqcut,mgfft,my_natom,natom,&
482 7423 : & nattyp,nfft,ngfft,nhat,nlstr,nspden,nsym,ntypat,optgr,optgr2,optstr,optstr2,&
483 7423 : & pawang,pawfgrtab,pawrhoij,pawtab,ph1d,psps,qphon,rprimd,symrec,typat,ucvol,vtrial,vxc,xred,&
484 7423 : & mpi_atmtab,comm_atom,comm_fft,mpi_comm_grid,me_g0,paral_kgb,distribfft) ! optional arguments (parallelism)
485 :
486 : !Arguments ------------------------------------
487 : !scalars
488 : integer,intent(in) :: dimnhat,dyfr_cplex,mgfft,my_natom,natom,nfft,nspden,nsym,ntypat
489 : integer,intent(in) :: optgr,optgr2,optstr,optstr2
490 : integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,mpi_comm_grid,paral_kgb
491 : real(dp),intent(in) :: gsqcut,ucvol
492 : type(distribfft_type),optional,target,intent(in) :: distribfft
493 : type(pawang_type),intent(in) :: pawang
494 : type(pseudopotential_type),intent(in) :: psps
495 : !arrays
496 : integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18)
497 : integer,intent(in) :: symrec(3,3,nsym),typat(natom)
498 : integer,optional,target,intent(in) :: mpi_atmtab(:)
499 : real(dp),intent(in) :: nhat(nfft,dimnhat),ph1d(2,3*(2*mgfft+1)*natom),qphon(3)
500 : real(dp),intent(in) :: rprimd(3,3),vxc(nfft,nspden),xred(3,natom)
501 : real(dp),intent(in),target :: vtrial(nfft,nspden)
502 : real(dp),intent(inout) :: dyfrnl(dyfr_cplex,3,3,natom,natom*optgr2)
503 : real(dp),intent(inout) :: eltfrnl(6+3*natom,6),grnl(3*natom*optgr)
504 : real(dp),intent(inout) :: nlstr(6*optstr)
505 : type(pawfgrtab_type),target,intent(inout) :: pawfgrtab(:)
506 : type(pawrhoij_type),target,intent(inout) :: pawrhoij(:)
507 : type(pawtab_type),intent(in) :: pawtab(ntypat)
508 :
509 : !Local variables-------------------------------
510 : !scalars
511 : integer :: bufind,bufsiz,cplex,dimvtrial,eps_alpha,eps_beta,eps_gamma,eps_delta,iatm,iatom
512 : integer :: iatom_pawfgrtab,iatom_pawrhoij,iatom_tot,iatshft,ic,idiag,idir,ier,ilm,indx,irhoij
513 : integer :: isel,ishift_grhoij,ishift_gr,ishift2_gr,ishift_gr2,ishift_str,ishift_str2,ishift_str2is,ispden
514 : integer :: ispvtr,itypat,jatom,jatom_tot,jatm,jc,jrhoij,jtypat,klm,klmn,klmn1,ll,lm_size
515 : integer :: lm_sizej,lmax,lmin,lmn2_size,me_fft,mu,mua,mub,mushift,my_me_g0,my_comm_atom,my_comm_fft
516 : integer :: my_comm_grid,my_paral_kgb,n1,n2,n3,nfftot,nfgd,nfgd_jatom
517 : integer :: ngrad,ngrad_nondiag,ngradp,ngradp_nondiag,ngrhat,nsploop
518 : integer :: opt1,opt2,opt3,qne0,usexcnhat
519 : logical,parameter :: save_memory=.true.
520 : logical :: has_phase,my_atmtab_allocated
521 : logical :: paral_atom,paral_atom_pawfgrtab,paral_atom_pawrhoij,paral_grid
522 : real(dp) :: dlt_tmp,fact_ucvol,grhat_x,hatstr_diag,rcut_jatom,ro,ro_d,ucvol_
523 : character(len=500) :: msg
524 : type(distribfft_type),pointer :: my_distribfft
525 : type(pawfgrtab_type),pointer :: pawfgrtab_iatom,pawfgrtab_jatom
526 : type(pawrhoij_type),pointer :: pawrhoij_iatom,pawrhoij_jatom
527 : !arrays
528 : integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
529 : integer,parameter :: eps1(6)=(/1,2,3,2,3,1/),eps2(6)=(/1,2,3,3,1,2/)
530 : integer,parameter :: mu9(9)=(/1,2,3,4,5,6,4,5,6/)
531 7423 : integer,allocatable :: atindx(:),atm_indx(:),mu4(:)
532 7423 : integer,allocatable,target :: ifftsph_tmp(:)
533 7423 : integer,ABI_CONTIGUOUS pointer :: ffti3_local(:),fftn3_distrib(:),ifft_jatom(:)
534 7423 : integer, pointer :: my_atmtab(:)
535 : real(dp) :: gmet(3,3),gprimd(3,3),hatstr(6),rdum(1),rdum2(1),rmet(3,3),tmp(12)
536 14846 : real(dp) :: work1(dyfr_cplex,3,3),work2(dyfr_cplex,3,3)
537 7423 : real(dp),allocatable :: buf(:,:),buf1(:),dyfr(:,:,:,:,:),eltfr(:,:)
538 7423 : real(dp),allocatable :: grhat_tmp(:,:),grhat_tmp2(:,:),hatgr(:)
539 7423 : real(dp),allocatable :: prod(:,:),prodp(:,:),vloc(:),vpsp1_gr(:,:),vpsp1_str(:,:)
540 7423 : real(dp),allocatable,target :: rfgd_tmp(:,:)
541 7423 : real(dp),ABI_CONTIGUOUS pointer :: gylm_jatom(:,:),gylmgr_jatom(:,:,:),gylmgr2_jatom(:,:,:),expiqr_jatom(:,:)
542 7423 : real(dp),ABI_CONTIGUOUS pointer :: rfgd_jatom(:,:),vtrial_(:,:)
543 7423 : type(coeff2_type),allocatable :: prod_nondiag(:),prodp_nondiag(:)
544 7423 : type(pawfgrtab_type),pointer :: pawfgrtab_(:),pawfgrtab_tot(:)
545 7423 : type(pawrhoij_type),pointer :: pawrhoij_(:),pawrhoij_tot(:)
546 :
547 : ! *************************************************************************
548 :
549 : DBG_ENTER("COLL")
550 :
551 : !Compatibility tests
552 14 : qne0=0;if (qphon(1)**2+qphon(2)**2+qphon(3)**2>=1.d-15) qne0=1
553 7423 : if (my_natom>0) then
554 6379 : if ((optgr2==1.or.optstr2==1).and.pawrhoij(1)%ngrhoij==0) then
555 0 : msg='pawgrnl: inconsistency between variables optgr2/optstr2 and ngrhoij!'
556 0 : ABI_BUG(msg)
557 : end if
558 6379 : if (pawfgrtab(1)%rfgd_allocated==0) then
559 321 : if ((optgr2==1.and.qne0==1).or.optstr2==1) then
560 0 : ABI_BUG('pawgrnl: pawfgrtab()%rfgd array must be allocated!')
561 : end if
562 : end if
563 6379 : if (pawrhoij(1)%qphase/=1) then
564 0 : ABI_BUG('pawgrnl: not supposed to be called with pawrhoij(:)%qphase=2!')
565 : end if
566 : end if
567 :
568 : !----------------------------------------------------------------------
569 : !Parallelism setup
570 :
571 : !Set up parallelism over atoms
572 7423 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
573 7423 : paral_atom_pawfgrtab=(size(pawfgrtab)/=natom)
574 7423 : paral_atom_pawrhoij=(size(pawrhoij)/=natom)
575 7423 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
576 7423 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
577 7423 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
578 7423 : if (paral_atom) then
579 9216 : ABI_MALLOC(atm_indx,(natom))
580 9428 : atm_indx=-1
581 5151 : do iatom=1,my_natom
582 5151 : atm_indx(my_atmtab(iatom))=iatom
583 : end do
584 : end if
585 :
586 : !Set up parallelism over real space grid and/or FFT
587 7423 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nfftot=n1*n2*n3
588 7423 : my_comm_grid=xmpi_comm_self;my_comm_fft=xmpi_comm_self;me_fft=0
589 7423 : my_me_g0=1;my_paral_kgb=0;paral_grid=.false.;nullify(my_distribfft)
590 7423 : if (present(mpi_comm_grid).or.present(comm_fft)) then
591 7366 : if (present(mpi_comm_grid)) my_comm_grid=mpi_comm_grid
592 7366 : if (present(comm_fft)) my_comm_fft=comm_fft
593 7366 : if (.not.present(mpi_comm_grid)) my_comm_grid=comm_fft
594 7366 : if (.not.present(comm_fft)) my_comm_fft=mpi_comm_grid
595 7366 : paral_grid=(xmpi_comm_size(my_comm_grid)>1)
596 7366 : me_fft=xmpi_comm_rank(my_comm_fft)
597 : end if
598 7423 : if (optgr2==1.or.optstr2==1) then
599 44 : if (present(comm_fft)) then
600 0 : if ((.not.present(paral_kgb)).or.(.not.present(me_g0)).or.(.not.present(distribfft))) then
601 0 : ABI_BUG(' Need paral_kgb, me_g0 and distribfft with comm_fft !')
602 : end if
603 0 : my_me_g0=me_g0;my_paral_kgb=paral_kgb
604 0 : my_distribfft => distribfft
605 : else
606 44 : ABI_MALLOC(my_distribfft,)
607 44 : call my_distribfft%init_seq('f',n2,n3,'fourdp')
608 : end if
609 44 : if (n2 == my_distribfft%n2_coarse) then
610 0 : fftn3_distrib => my_distribfft%tab_fftdp3_distrib
611 0 : ffti3_local => my_distribfft%tab_fftdp3_local
612 : else
613 44 : fftn3_distrib => my_distribfft%tab_fftdp3dg_distrib
614 44 : ffti3_local => my_distribfft%tab_fftdp3dg_local
615 : end if
616 : else
617 : nullify(my_distribfft,fftn3_distrib,ffti3_local)
618 : end if
619 :
620 : !----------------------------------------------------------------------
621 : !Initializations
622 :
623 : !Compute different geometric tensors
624 : !ucvol is not computed here but provided as input arg
625 7423 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol_)
626 7423 : fact_ucvol=ucvol/dble(nfftot)
627 :
628 : !Retrieve local potential according to the use of nhat in XC
629 18858 : usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
630 7423 : if (usexcnhat==0) then
631 : ! dimvtrial=nspden
632 5085 : dimvtrial=1
633 15255 : ABI_MALLOC(vtrial_,(nfft,dimvtrial))
634 : !!!$OMP PARALLEL DO PRIVATE(ic) SHARED(nfft,vtrial,vtrial_,vxc)
635 131385402 : do ic=1,nfft
636 525526353 : vtrial_(ic,1:dimvtrial)=vtrial(ic,1:dimvtrial)-vxc(ic,1:dimvtrial)
637 : end do
638 : else
639 2338 : dimvtrial=nspden
640 2338 : vtrial_ => vtrial
641 : end if
642 :
643 : !Initializations and allocations
644 7423 : ngrhat=0;ngrad=0;ngradp=0;ngrad_nondiag=0;ngradp_nondiag=0
645 7423 : ishift_grhoij=0;ishift_gr=0;ishift_gr2=0;ishift_str=0;ishift_str2=0;ishift_str2is=0;ishift2_gr=0
646 7423 : cplex=1;if (qne0==1) cplex=2
647 7423 : if (optgr==1) then
648 26776 : ABI_MALLOC(hatgr,(3*natom))
649 71029 : hatgr=zero
650 : ngrad=ngrad+3
651 : ngrhat=ngrhat+3
652 : ishift_gr2=ishift_gr2+3
653 : end if
654 7423 : if (optgr2==1) then
655 44 : mu=min(dyfr_cplex,cplex)
656 44 : ngrad =ngrad +9
657 44 : ngradp=ngradp+3
658 44 : ngrad_nondiag = ngrad_nondiag +9*mu
659 44 : ngradp_nondiag= ngradp_nondiag+3*mu
660 44 : ngrhat= ngrhat+9*mu
661 : end if
662 7423 : if (optstr==1) then
663 1395 : hatstr=zero
664 1395 : ngrad=ngrad+6
665 1395 : ngrhat=ngrhat+6
666 1395 : ishift_gr=ishift_gr+6
667 1395 : ishift_gr2=ishift_gr2+6
668 1395 : ishift_str2=ishift_str2+6
669 1395 : ishift_str2is = ishift_str2is+6
670 : end if
671 7423 : if (optstr2==1) then
672 12 : ngrad =ngrad+6*(6+3)
673 12 : ngradp=ngradp+(6+3)
674 12 : ngrad_nondiag =ngrad_nondiag+6*(6+3)
675 12 : ngradp_nondiag=ngradp_nondiag+3
676 12 : ishift2_gr=ishift2_gr+3
677 12 : ngrhat=ngrhat+6*(6+3)
678 12 : ishift_gr=ishift_gr+(6+3)
679 12 : ishift_gr2=ishift_gr2+6*(6+3)
680 12 : ishift_str2is=ishift_str2is+36
681 12 : ishift_grhoij = 6
682 : end if
683 :
684 : !DEBUG
685 : ! write(6,*)' preparatory computations : usexcnhat, nspden, dimvtrial=',usexcnhat, nspden, dimvtrial
686 : !ENDDEBUG
687 : !nsploop=nspden;if (dimvtrial<nspden) nsploop=2
688 7423 : nsploop=nspden;if (dimvtrial<nspden .and. nspden==4) nsploop=1
689 7423 : if (optgr2/=1.and.optstr2/=1) then
690 22137 : ABI_MALLOC(grhat_tmp,(ngrhat,1))
691 : else
692 176 : ABI_MALLOC(grhat_tmp,(ngrhat,natom))
693 3866 : grhat_tmp=zero
694 264 : ABI_MALLOC(prod_nondiag,(natom))
695 220 : ABI_MALLOC(prodp_nondiag,(natom))
696 132 : ABI_MALLOC(atindx,(natom))
697 44 : if(optgr2==1.or.optstr2==1)then
698 132 : ABI_MALLOC(vpsp1_gr,(cplex*nfft,3))
699 3449624 : vpsp1_gr(:,:)= zero
700 : end if
701 44 : if (optgr2==1) then
702 264 : ABI_MALLOC(dyfr,(dyfr_cplex,3,3,natom,natom))
703 13650 : dyfr=zero
704 : end if
705 44 : if (optstr2==1) then
706 36 : ABI_MALLOC(vpsp1_str,(cplex*nfft,6))
707 36 : ABI_MALLOC(grhat_tmp2,(18,natom))
708 60 : ABI_MALLOC(eltfr,(6+3*natom,6))
709 1272 : eltfr=zero
710 : end if
711 44 : ABI_MALLOC(mu4,(4))
712 176 : atindx(:)=0
713 176 : do iatom=1,natom
714 : iatm=0
715 538 : do while (atindx(iatom)==0.and.iatm<natom)
716 494 : iatm=iatm+1;if (atindx1(iatm)==iatom) atindx(iatom)=iatm
717 : end do
718 : end do
719 : end if
720 :
721 : !The computation of dynamical matrix and elastic tensor requires the knowledge of
722 : !g_l(r-R).Y_lm(r-R) and derivatives for all atoms
723 : !Compute them here, except memory saving is activated
724 : if ((.not.save_memory).and.(optgr2==1.or.optstr2==1)) then
725 : do jatom=1,size(pawfgrtab)
726 : jatom_tot=jatom;if (paral_atom_pawfgrtab) jatom_tot=my_atmtab(jatom)
727 : pawfgrtab_jatom => pawfgrtab(jatom)
728 : lm_sizej=pawfgrtab_jatom%l_size**2
729 : opt1=0;opt2=0;opt3=0
730 : if (pawfgrtab_jatom%gylm_allocated==0) then
731 : if (allocated(pawfgrtab_jatom%gylm)) then
732 : ABI_FREE(pawfgrtab_jatom%gylm)
733 : end if
734 : ABI_MALLOC(pawfgrtab_jatom%gylm,(pawfgrtab_jatom%nfgd,lm_sizej))
735 : pawfgrtab_jatom%gylm_allocated=2;opt1=1
736 : end if
737 : if (pawfgrtab_jatom%gylmgr_allocated==0) then
738 : if (allocated(pawfgrtab_jatom%gylmgr)) then
739 : ABI_FREE(pawfgrtab_jatom%gylmgr)
740 : end if
741 : ABI_MALLOC(pawfgrtab_jatom%gylmgr,(3,pawfgrtab_jatom%nfgd,lm_sizej))
742 : pawfgrtab_jatom%gylmgr_allocated=2;opt2=1
743 : end if
744 : if (opt1+opt2+opt3>0) then
745 : call pawgylm(pawfgrtab_jatom%gylm,pawfgrtab_jatom%gylmgr,&
746 : & pawfgrtab_jatom%gylmgr2,lm_sizej,pawfgrtab_jatom%nfgd,&
747 : & opt1,opt2,opt3,pawtab(typat(jatom_tot)),pawfgrtab_jatom%rfgd)
748 : end if
749 : if (optgr2==1.and.qne0==1) then
750 : if (pawfgrtab_jatom%expiqr_allocated==0) then
751 : if (allocated(pawfgrtab_jatom%expiqr)) then
752 : ABI_FREE(pawfgrtab_jatom%expiqr)
753 : end if
754 : pawfgrtab_jatom%expiqr_allocated=2
755 : ABI_MALLOC(pawfgrtab_jatom%expiqr,(2,nfgd))
756 : call pawexpiqr(pawfgrtab_jatom%expiqr,gprimd,pawfgrtab_jatom%nfgd,&
757 : & qphon,pawfgrtab_jatom%rfgd,xred(:,jatom_tot))
758 : end if
759 : end if
760 : end do
761 : end if
762 :
763 : !The computation of dynamical matrix and elastic tensor might require some communications
764 : if ((optgr2==1.or.optstr2==1).and.paral_atom.and.paral_atom_pawfgrtab.and.(.not.save_memory)) then
765 : ABI_MALLOC(pawfgrtab_tot,(natom))
766 : call pawfgrtab_nullify(pawfgrtab_tot)
767 : call pawfgrtab_gather(pawfgrtab,pawfgrtab_tot,my_comm_atom,ier,mpi_atmtab=my_atmtab)
768 : else
769 7423 : pawfgrtab_tot => pawfgrtab
770 : end if
771 7423 : if ((optgr2==1.or.optstr2==1).and.paral_atom.and.paral_atom_pawrhoij) then
772 0 : ABI_MALLOC(pawrhoij_tot,(natom))
773 0 : call pawrhoij_nullify(pawrhoij_tot)
774 : call pawrhoij_gather(pawrhoij,pawrhoij_tot,-1,my_comm_atom, &
775 0 : & with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.)
776 : else
777 7423 : pawrhoij_tot => pawrhoij
778 : end if
779 :
780 : if (save_memory) then
781 7423 : pawfgrtab_ => pawfgrtab
782 7423 : pawrhoij_ => pawrhoij
783 : else
784 : pawfgrtab_ => pawfgrtab_tot
785 : pawrhoij_ => pawrhoij_tot
786 : end if
787 :
788 : !----------------------------------------------------------------------
789 : !Loops over types and atoms
790 :
791 7423 : iatshft=0
792 18858 : do itypat=1,ntypat
793 :
794 11435 : lmn2_size=pawtab(itypat)%lmn2_size
795 11435 : lm_size=pawtab(itypat)%lcut_size**2
796 :
797 34803 : do iatm=iatshft+1,iatshft+nattyp(itypat)
798 :
799 23368 : iatom_tot=atindx1(iatm)
800 23368 : iatom=iatom_tot
801 23368 : if (paral_atom) then
802 6356 : if (save_memory.or.(optgr2/=1.and.optstr2/=1)) iatom=atm_indx(iatom_tot)
803 : end if
804 :
805 23368 : if (iatom==-1) cycle
806 19091 : iatom_pawfgrtab=iatom_tot;if (paral_atom_pawfgrtab) iatom_pawfgrtab=iatom
807 19091 : iatom_pawrhoij =iatom_tot;if (paral_atom_pawrhoij) iatom_pawrhoij =iatom
808 19091 : pawfgrtab_iatom => pawfgrtab_(iatom_pawfgrtab)
809 19091 : pawrhoij_iatom => pawrhoij_(iatom_pawrhoij)
810 :
811 19091 : idiag=1;if (optgr2==1.or.optstr2==1) idiag=iatm
812 19091 : nfgd=pawfgrtab_iatom%nfgd
813 :
814 57273 : ABI_MALLOC(vloc,(nfgd))
815 19091 : if (ngrad>0) then
816 76284 : ABI_MALLOC(prod,(ngrad,lm_size))
817 : end if
818 19091 : if (ngradp>0) then
819 284 : ABI_MALLOC(prodp,(ngradp,lm_size))
820 : end if
821 19091 : if (ngrad_nondiag>0.and.ngradp_nondiag>0) then
822 320 : do jatm=1,natom
823 249 : jtypat=typat(atindx1(jatm))
824 249 : lm_sizej=pawtab(jtypat)%lcut_size**2
825 996 : ABI_MALLOC(prod_nondiag(jatm)%value,(ngrad_nondiag,lm_sizej))
826 996 : ABI_MALLOC(prodp_nondiag(jatm)%value,(ngradp_nondiag,lm_sizej))
827 74601 : prod_nondiag(jatm)%value=zero
828 13496 : prodp_nondiag(jatm)%value=zero
829 : end do
830 : end if
831 :
832 114776 : grhat_tmp=zero
833 20979 : if(optstr2==1) grhat_tmp2=zero
834 :
835 : ! ------------------------------------------------------------------
836 : ! Compute some useful data
837 :
838 : ! Eventually compute g_l(r).Y_lm(r) derivatives for the current atom (if not already done)
839 19091 : if ((optgr==1.or.optstr==1).and.(optgr2/=1).and.(optstr2/=1)) then
840 19000 : if (pawfgrtab_iatom%gylmgr_allocated==0) then
841 1 : if (allocated(pawfgrtab_iatom%gylmgr)) then
842 1 : ABI_FREE(pawfgrtab_iatom%gylmgr)
843 : end if
844 4 : ABI_MALLOC(pawfgrtab_iatom%gylmgr,(3,pawfgrtab_iatom%nfgd,lm_size))
845 1 : pawfgrtab_iatom%gylmgr_allocated=2
846 : call pawgylm(rdum,pawfgrtab_iatom%gylmgr,rdum2,lm_size,pawfgrtab_iatom%nfgd,&
847 1 : & 0,1,0,pawtab(itypat),pawfgrtab_iatom%rfgd)
848 : end if
849 :
850 : end if
851 19091 : if (optgr2==1.or.optstr2==1) then
852 71 : opt1=0;opt2=0;opt3=0
853 71 : if (pawfgrtab_iatom%gylm_allocated==0) then
854 0 : if (allocated(pawfgrtab_iatom%gylm)) then
855 0 : ABI_FREE(pawfgrtab_iatom%gylm)
856 : end if
857 0 : ABI_MALLOC(pawfgrtab_iatom%gylm,(pawfgrtab_iatom%nfgd,lm_size))
858 0 : pawfgrtab_iatom%gylm_allocated=2;opt1=1
859 : end if
860 71 : if (pawfgrtab_iatom%gylmgr_allocated==0) then
861 0 : if (allocated(pawfgrtab_iatom%gylmgr)) then
862 0 : ABI_FREE(pawfgrtab_iatom%gylmgr)
863 : end if
864 0 : ABI_MALLOC(pawfgrtab_iatom%gylmgr,(3,pawfgrtab_iatom%nfgd,lm_size))
865 0 : pawfgrtab_iatom%gylmgr_allocated=2;opt2=1
866 : end if
867 71 : if (pawfgrtab_iatom%gylmgr2_allocated==0) then
868 0 : if (allocated(pawfgrtab_iatom%gylmgr2)) then
869 0 : ABI_FREE(pawfgrtab_iatom%gylmgr2)
870 : end if
871 0 : ABI_MALLOC(pawfgrtab_iatom%gylmgr2,(6,pawfgrtab_iatom%nfgd,lm_size))
872 0 : pawfgrtab_iatom%gylmgr2_allocated=2;opt3=1
873 : end if
874 71 : if (opt1+opt2+opt3>0) then
875 : call pawgylm(pawfgrtab_iatom%gylm,pawfgrtab_iatom%gylmgr,&
876 : & pawfgrtab_iatom%gylmgr2,lm_size,pawfgrtab_iatom%nfgd,&
877 0 : & opt1,opt2,opt3,pawtab(itypat),pawfgrtab_iatom%rfgd)
878 : end if
879 : end if
880 :
881 : ! Eventually compute exp(-i.q.r) factors for the current atom (if not already done)
882 19091 : if (optgr2==1.and.qne0==1.and.(pawfgrtab_iatom%expiqr_allocated==0)) then
883 0 : if (allocated(pawfgrtab_iatom%expiqr)) then
884 0 : ABI_FREE(pawfgrtab_iatom%expiqr)
885 : end if
886 0 : ABI_MALLOC(pawfgrtab_iatom%expiqr,(2,nfgd))
887 : call pawexpiqr(pawfgrtab_iatom%expiqr,gprimd,nfgd,qphon,&
888 0 : & pawfgrtab_iatom%rfgd,xred(:,iatom))
889 0 : pawfgrtab_iatom%expiqr_allocated=2
890 : end if
891 19091 : has_phase=(optgr2==1.and.pawfgrtab_iatom%expiqr_allocated/=0)
892 :
893 : ! Eventually compute 1st-order potential
894 19091 : if (optgr2==1.or.optstr2==1) then
895 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,iatom_tot,&
896 : & mgfft,psps%mqgrid_vl,natom,3,nfft,ngfft,ntypat,ph1d,&
897 : & psps%qgrid_vl,qphon,typat,ucvol,psps%usepaw,xred,psps,pawtab,atmvlocr1=vpsp1_gr,&
898 : & vspl=psps%vlspl,comm_fft=my_comm_fft,me_g0=my_me_g0,&
899 71 : & paral_kgb=my_paral_kgb,distribfft=my_distribfft)
900 71 : if (cplex==1) then
901 1630953 : do ic=1,nfft
902 6523584 : tmp(1:3)=vpsp1_gr(ic,1:3)
903 6523641 : do mu=1,3
904 6523584 : vpsp1_gr(ic,mu)=-(gprimd(mu,1)*tmp(1)+gprimd(mu,2)*tmp(2)+gprimd(mu,3)*tmp(3))
905 : end do
906 : end do
907 : else ! cplex=2
908 177566 : do ic=1,nfft
909 1242864 : jc=2*ic;tmp(1:3)=vpsp1_gr(jc-1,1:3);tmp(4:6)=vpsp1_gr(jc,1:3)
910 729242 : do mu=1,3
911 532656 : vpsp1_gr(jc-1,mu)=-(gprimd(mu,1)*tmp(1)+gprimd(mu,2)*tmp(2)+gprimd(mu,3)*tmp(3))
912 710208 : vpsp1_gr(jc ,mu)=-(gprimd(mu,1)*tmp(4)+gprimd(mu,2)*tmp(5)+gprimd(mu,3)*tmp(6))
913 : end do
914 : end do
915 : end if
916 : end if
917 19091 : if (optstr2==1) then
918 1469846 : vpsp1_str(:,:) = zero
919 : call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,natom+3,&
920 : & mgfft,psps%mqgrid_vl,natom,6,nfft,ngfft,ntypat,&
921 : & ph1d,psps%qgrid_vl,qphon,typat,ucvol,psps%usepaw,xred,psps,pawtab,atmvlocr1=vpsp1_str,&
922 : & vspl=psps%vlspl,comm_fft=my_comm_fft,me_g0=my_me_g0,&
923 26 : & paral_kgb=my_paral_kgb,distribfft=my_distribfft)
924 : end if
925 :
926 : ! ------------------------------------------------------------------
927 : ! Loop over spin components
928 :
929 40010 : do ispden=1,nsploop
930 :
931 : ! ----- Retrieve potential (subtle if nspden=4 ;-)
932 20919 : if (nspden/=4) then
933 20034 : ispvtr=min(dimvtrial,ispden)
934 19361864 : do ic=1,nfgd
935 19341830 : jc = pawfgrtab_iatom%ifftsph(ic)
936 19361864 : vloc(ic)=vtrial_(jc,ispvtr)
937 : end do
938 : else
939 885 : if (ispden==1) then
940 237 : ispvtr=min(dimvtrial,2)
941 399818 : do ic=1,nfgd
942 399581 : jc=pawfgrtab_iatom%ifftsph(ic)
943 399818 : vloc(ic)=half*(vtrial_(jc,1)+vtrial_(jc,ispvtr))
944 : end do
945 648 : else if (ispden==4) then
946 216 : ispvtr=min(dimvtrial,2)
947 338150 : do ic=1,nfgd
948 337934 : jc=pawfgrtab_iatom%ifftsph(ic)
949 338150 : vloc(ic)=half*(vtrial_(jc,1)-vtrial_(jc,ispvtr))
950 : end do
951 432 : else if (ispden==2) then
952 216 : ispvtr=min(dimvtrial,3)
953 338150 : do ic=1,nfgd
954 337934 : jc=pawfgrtab_iatom%ifftsph(ic)
955 338150 : vloc(ic)=vtrial_(jc,ispvtr)
956 : end do
957 : else ! ispden=3
958 216 : ispvtr=min(dimvtrial,4)
959 338150 : do ic=1,nfgd
960 337934 : jc=pawfgrtab_iatom%ifftsph(ic)
961 338150 : vloc(ic)=-vtrial_(jc,ispvtr)
962 : end do
963 : end if
964 : end if
965 :
966 : ! -----------------------------------------------------------------------
967 : ! ----- Compute projected scalars (integrals of vloc and Q_ij^hat) ------
968 : ! ----- and/or their derivatives ----------------------------------------
969 :
970 867127 : if (ngrad>0) prod=zero
971 26373 : if (ngradp>0) prodp=zero
972 :
973 : ! ==== Contribution to forces ====
974 20919 : if (optgr==1) then
975 161970 : do ilm=1,lm_size
976 230940249 : do ic=1,pawfgrtab_iatom%nfgd
977 923255896 : do mu=1,3
978 : prod(mu+ishift_gr,ilm)=prod(mu+ishift_gr,ilm)-&
979 923113116 : & vloc(ic)*pawfgrtab_iatom%gylmgr(mu,ic,ilm)
980 : end do
981 : end do
982 : end do
983 : end if ! optgr
984 :
985 : ! ==== Contribution to stresses ====
986 20919 : if (optstr==1) then
987 43360 : do ilm=1,lm_size
988 144250062 : do ic=1,pawfgrtab_iatom%nfgd
989 144206702 : jc=pawfgrtab_iatom%ifftsph(ic)
990 1009487029 : do mu=1,6
991 865240212 : mua=alpha(mu);mub=beta(mu)
992 : prod(mu+ishift_str,ilm)=prod(mu+ishift_str,ilm) &
993 : & +half*vloc(ic)&
994 : & *(pawfgrtab_iatom%gylmgr(mua,ic,ilm)*pawfgrtab_iatom%rfgd(mub,ic)&
995 1009446914 : & +pawfgrtab_iatom%gylmgr(mub,ic,ilm)*pawfgrtab_iatom%rfgd(mua,ic))
996 : end do
997 : end do
998 : end do
999 : end if ! optstr
1000 : !DEBUG
1001 : ! write(6,*)' after loops on ilm, ic, mu : ispden,lm_size, pawfgrtab_iatom%nfgd=',ispden,lm_size, pawfgrtab_iatom%nfgd
1002 : ! write(6,*)' after loops on ilm, ic, mu, writes ilm, prod(1+ishift_str,ilm:lm_size) when bigger than tol10 (ilm between 1 and lm_size)'
1003 : ! do ilm=1, lm_size
1004 : ! if( abs(prod(1+ishift_str,ilm))>tol6 )then
1005 : ! write(6,*)ilm,prod(1+ishift_str,ilm)
1006 : ! endif
1007 : ! enddo
1008 : !ENDDEBUG
1009 :
1010 : ! ==== Diagonal contribution to frozen wf part of dyn. matrix ====
1011 20919 : if (optgr2==1) then
1012 : ! Diagonal contribution
1013 914 : do ilm=1,lm_size
1014 1514434 : do ic=1,pawfgrtab_iatom%nfgd
1015 15135200 : do mu=1,9
1016 : prod(ishift_gr2+mu,ilm)=prod(ishift_gr2+mu,ilm) &
1017 15135200 : & +half*vloc(ic)*pawfgrtab_iatom%gylmgr2(mu9(mu),ic,ilm)
1018 : end do
1019 6054917 : do mu=1,3
1020 : prodp(ishift_gr+mu,ilm)=prodp(ishift_gr+mu,ilm) &
1021 6054080 : & -vloc(ic)*pawfgrtab_iatom%gylmgr(mu,ic,ilm)
1022 : end do
1023 : end do
1024 : end do
1025 : end if ! optgr2
1026 :
1027 : ! ==== Diagonal contribution to elastic tensor ====
1028 20919 : if (optstr2==1) then
1029 260 : do ilm=1,lm_size
1030 216278 : do ic=1,pawfgrtab_iatom%nfgd
1031 216018 : mu=1
1032 216018 : jc=pawfgrtab_iatom%ifftsph(ic)
1033 1512126 : do mua=1,6
1034 1296108 : eps_alpha=eps1(mua);eps_beta=eps2(mua);
1035 9288774 : do mub=1,6
1036 7776648 : eps_gamma=eps1(mub);eps_delta=eps2(mub);
1037 38883240 : mu4 = zero
1038 7776648 : call pawgrnl_convert(mu4,eps_alpha,eps_beta,eps_gamma,eps_delta)
1039 : ! v_loc*d2glylm
1040 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) + half*half*vloc(ic)*( &
1041 : & pawfgrtab_iatom%rfgd(eps_beta,ic)*pawfgrtab_iatom%rfgd(eps_gamma,ic)*&
1042 : & pawfgrtab_iatom%gylmgr2(mu4(2),ic,ilm)&
1043 : & +pawfgrtab_iatom%rfgd(eps_alpha,ic)*pawfgrtab_iatom%rfgd(eps_gamma,ic)*&
1044 : pawfgrtab_iatom%gylmgr2(mu4(4),ic,ilm)&
1045 : & +pawfgrtab_iatom%rfgd(eps_beta,ic) *pawfgrtab_iatom%rfgd(eps_delta,ic)*&
1046 : pawfgrtab_iatom%gylmgr2(mu4(1),ic,ilm)&
1047 : & +pawfgrtab_iatom%rfgd(eps_alpha,ic)*pawfgrtab_iatom%rfgd(eps_delta,ic)*&
1048 7776648 : pawfgrtab_iatom%gylmgr2(mu4(3),ic,ilm))
1049 7776648 : if(eps_gamma==eps_beta)then
1050 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1051 2592216 : & +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic))
1052 : end if
1053 7776648 : if(eps_gamma==eps_alpha)then
1054 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1055 2592216 : & +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic))
1056 : end if
1057 7776648 : if(eps_delta==eps_beta)then
1058 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1059 2592216 : & +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic))
1060 : end if
1061 7776648 : if(eps_delta==eps_alpha)then
1062 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1063 2592216 : & +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic))
1064 : end if
1065 : ! d(vloc)/d(eps_gammadelta) * d(gylm)/d(eps_alphabeta)
1066 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm)&
1067 : & +vpsp1_str(jc,mub)*half*(&
1068 : & (pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic)&
1069 7776648 : & +pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm) *pawfgrtab_iatom%rfgd(eps_alpha,ic)))
1070 : ! d(vloc)/d(eps_alphabeta) * d(gylm)/d(eps_gammadelta)
1071 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm)&
1072 : & +vpsp1_str(jc,mua)*half*(&
1073 : & (pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_delta,ic)&
1074 7776648 : & +pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_gamma,ic)))
1075 : ! delta_alphabeta * dv_loc/depsgammadelta * (gylm)
1076 7776648 : if (mua<=3) then
1077 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1078 3888324 : & +vpsp1_str(jc,mub)*pawfgrtab_iatom%gylm(ic,ilm)
1079 : end if
1080 : ! delta_gammadelta * dv_loc/depsalphabeta * (gylm)
1081 7776648 : if (mub<=3) then
1082 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1083 3888324 : & +vpsp1_str(jc,mua) * pawfgrtab_iatom%gylm(ic,ilm)
1084 : end if
1085 : ! delta_gammadelta * v_loc * d(gylm)/d(eps_alphabeta)
1086 : if (mub<=3) then
1087 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1088 : & +half*vloc(ic)&
1089 : & *(pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic)&
1090 3888324 : & + pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic))
1091 : end if
1092 : ! delta_alphabeta * v_loc * d(gylm)/d(eps_gammadelta)
1093 7776648 : if (mua<=3) then
1094 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1095 : & +half*vloc(ic)&
1096 : & *(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_delta,ic)&
1097 3888324 : & + pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_gamma,ic))
1098 : end if
1099 : ! delta_gammadelta delta_alphabeta * v_loc * (gylm)
1100 7776648 : if (mua<=3.and.mub<=3) then
1101 : prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1102 1944162 : & +vloc(ic)*pawfgrtab_iatom%gylm(ic,ilm)
1103 : end if
1104 9072756 : mu=mu+1
1105 : end do !end loop mub
1106 : end do !end loop mua
1107 : ! vloc * d(gylm)/d(eps_alphabeta)
1108 1512126 : do mu=1,6
1109 1296108 : mua=alpha(mu);mub=beta(mu)
1110 : prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)&
1111 : & +half*vloc(ic)*&
1112 : & (pawfgrtab_iatom%gylmgr(mua,ic,ilm)*pawfgrtab_iatom%rfgd(mub,ic)&
1113 1296108 : & +pawfgrtab_iatom%gylmgr(mub,ic,ilm)*pawfgrtab_iatom%rfgd(mua,ic))
1114 : ! d(vloc)/d(eps_alphabeta or gammadelta) * gylm
1115 : prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)&
1116 1296108 : & +vpsp1_str(jc,mu)*pawfgrtab_iatom%gylm(ic,ilm)
1117 : ! delta_alphabeta * vloc * gylm
1118 1296108 : if (mu<=3) then
1119 : prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)&
1120 648054 : & +vloc(ic)*pawfgrtab_iatom%gylm(ic,ilm)
1121 : end if
1122 :
1123 : ! INTERNAL STRAIN CONTRIBUTION:
1124 5400450 : do idir=1,3
1125 : ! v_loc*d2glylm/dR contribution:
1126 3888324 : eps_alpha=alpha(mu);eps_beta=beta(mu);
1127 3888324 : call pawgrnl_convert(mu4,eps_alpha,eps_beta,idir,idir)
1128 : prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)&
1129 : & -half*vloc(ic)&
1130 : & *(pawfgrtab_iatom%gylmgr2(mu4(3),ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic)&
1131 3888324 : & +pawfgrtab_iatom%gylmgr2(mu4(1),ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic))
1132 3888324 : if (idir==eps_beta)then
1133 : prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)&
1134 1296108 : & -half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm))
1135 : end if
1136 3888324 : if (idir==eps_alpha)then
1137 : prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)&
1138 1296108 : & -half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm))
1139 : end if
1140 : ! delta_gammadelta * v_loc * d(gylm)/dR
1141 3888324 : if (mu<=3) then
1142 : prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)-&
1143 1944162 : vloc(ic)*pawfgrtab_iatom%gylmgr(idir,ic,ilm)
1144 : end if
1145 : ! dv_loc/deps_alph_beta * d(gylm)/dR
1146 : prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)-&
1147 5184432 : vpsp1_str(jc,mu)*pawfgrtab_iatom%gylmgr(idir,ic,ilm)
1148 : end do
1149 : end do
1150 864306 : do idir=1,3
1151 : ! v_loc * d(gylm)/dR
1152 864072 : prodp(6+idir,ilm) = prodp(6+idir,ilm)-vloc(ic)*pawfgrtab_iatom%gylmgr(idir,ic,ilm)
1153 : end do !end loop idir
1154 : ! END INTERNAL STRAIN CONTRIBUTION
1155 :
1156 : end do
1157 : end do
1158 : end if !optstr2
1159 :
1160 : ! Off-diagonal contributions
1161 20919 : if (optgr2==1.or.optstr2==1) then
1162 336 : do jatm=1,natom
1163 259 : jatom_tot=atindx1(jatm);jtypat=typat(jatom_tot)
1164 259 : jatom=jatom_tot;if (paral_atom.and.save_memory) jatom=atm_indx(jatom_tot)
1165 259 : lm_sizej=pawtab(jtypat)%lcut_size**2
1166 :
1167 : ! Retrieve data for the atom j
1168 259 : if (save_memory.and.jatom/=iatom) then
1169 182 : rcut_jatom=pawtab(jtypat)%rshp
1170 : call pawrfgd_fft(ifftsph_tmp,gmet,n1,n2,n3,nfgd_jatom,rcut_jatom,rfgd_tmp,rprimd,&
1171 182 : & ucvol,xred(:,jatom_tot),fft_distrib=fftn3_distrib,fft_index=ffti3_local,me_fft=me_fft)
1172 182 : ifft_jatom => ifftsph_tmp ; rfgd_jatom => rfgd_tmp
1173 728 : ABI_MALLOC(gylm_jatom,(nfgd_jatom,lm_sizej))
1174 728 : ABI_MALLOC(gylmgr_jatom,(3,nfgd_jatom,lm_sizej))
1175 182 : opt1=1;opt2=1;opt3=0;gylmgr2_jatom=>gylmgr_jatom
1176 : call pawgylm(gylm_jatom,gylmgr_jatom,gylmgr2_jatom,lm_sizej,nfgd_jatom,&
1177 182 : & opt1,opt2,opt3,pawtab(typat(jatom_tot)),rfgd_jatom)
1178 182 : if (optgr2==1.and.qne0==1) then
1179 36 : ABI_MALLOC(expiqr_jatom,(2,nfgd_jatom))
1180 12 : call pawexpiqr(expiqr_jatom,gprimd,nfgd_jatom,qphon,rfgd_jatom,xred(:,jatom_tot))
1181 : end if
1182 : else
1183 77 : pawfgrtab_jatom => pawfgrtab_tot(jatom)
1184 77 : nfgd_jatom = pawfgrtab_jatom%nfgd
1185 77 : ifft_jatom => pawfgrtab_jatom%ifftsph
1186 77 : rfgd_jatom => pawfgrtab_jatom%rfgd
1187 77 : gylm_jatom => pawfgrtab_jatom%gylm
1188 77 : gylmgr_jatom => pawfgrtab_jatom%gylmgr
1189 77 : gylmgr2_jatom => pawfgrtab_jatom%gylmgr2
1190 77 : expiqr_jatom => pawfgrtab_jatom%expiqr
1191 : end if
1192 :
1193 : ! ==== Off-diagonal contribution to frozen wf part of dyn. matrix ====
1194 259 : if (optgr2==1) then
1195 259 : mu = min(dyfr_cplex,cplex)
1196 28653 : prod_nondiag(jatm)%value(ishift_gr2+1:ishift_gr2+(9*mu),:) = zero
1197 11427 : prodp_nondiag(jatm)%value(ishift2_gr+1:ishift2_gr +(3*mu),:) = zero
1198 259 : if (has_phase.or.cplex==2) then
1199 28 : if (dyfr_cplex==1.or.cplex==1) then
1200 0 : do ilm=1,lm_sizej
1201 0 : do ic=1,nfgd_jatom
1202 0 : jc=2*ifft_jatom(ic)
1203 : tmp(1:3)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(1,ic) &
1204 0 : & -vpsp1_gr(jc ,1:3)*expiqr_jatom(2,ic)
1205 0 : do mu=1,9
1206 0 : mua=alpha(mu);mub=beta(mu)
1207 : prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) &
1208 0 : & +tmp(mua)*gylmgr_jatom(mub,ic,ilm)
1209 : end do
1210 0 : do mu=1,3
1211 : prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) &
1212 0 : & -tmp(mu)*gylm_jatom(ic,ilm)
1213 : end do
1214 : end do
1215 : end do
1216 : else
1217 344 : do ilm=1,lm_sizej
1218 649584 : do ic=1,nfgd_jatom
1219 649240 : jc=2*ifft_jatom(ic)
1220 : tmp(1:3)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(1,ic) &
1221 2596960 : & -vpsp1_gr(jc ,1:3)*expiqr_jatom(2,ic)
1222 : tmp(4:6)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(2,ic) &
1223 2596960 : & +vpsp1_gr(jc ,1:3)*expiqr_jatom(1,ic)
1224 6492400 : do mu=1,9
1225 5843160 : mua=alpha(mu);mub=beta(mu)
1226 : prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) &
1227 5843160 : & +tmp(mua )*gylmgr_jatom(mub,ic,ilm)
1228 : prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm) &
1229 6492400 : & +tmp(3+mua)*gylmgr_jatom(mub,ic,ilm)
1230 : end do
1231 2597276 : do mu=1,3
1232 : prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) &
1233 1947720 : & -tmp( mu)*gylm_jatom(ic,ilm)
1234 : prodp_nondiag(jatm)%value(ishift2_gr+3+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+3+mu,ilm) &
1235 2596960 : & -tmp(3+mu)*gylm_jatom(ic,ilm)
1236 : end do
1237 : end do
1238 : end do
1239 : end if
1240 : else ! no phase
1241 2470 : do ilm=1,lm_sizej
1242 3028767 : do ic=1,nfgd_jatom
1243 3026297 : jc=ifft_jatom(ic)
1244 30262970 : do mu=1,9
1245 27236673 : mua=alpha(mu);mub=beta(mu)
1246 : prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) &
1247 30262970 : & +vpsp1_gr(jc,mua)*gylmgr_jatom(mub,ic,ilm)
1248 : end do
1249 12107427 : do mu=1,3
1250 : prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) &
1251 12105188 : & -vpsp1_gr(jc,mu)*gylm_jatom(ic,ilm)
1252 : end do
1253 : end do
1254 : end do
1255 : end if
1256 : end if ! optgr2
1257 :
1258 : ! ==== Off-diagonal contribution to elastic tensor ====
1259 259 : if (optstr2==1) then
1260 16856 : prod_nondiag(jatm)%value(ishift_str2is+1:ishift_str2is+18,:)=zero
1261 3626 : prodp_nondiag(jatm)%value(1:3,:)=zero
1262 980 : do ilm=1,lm_sizej
1263 631718 : do ic=1,nfgd_jatom
1264 630738 : mu=1;jc=ifft_jatom(ic)
1265 : ! INTERNAL STRAIN CONTRIBUTION:
1266 4415166 : do mua=1,6
1267 3784428 : eps_alpha=eps1(mua);eps_beta=eps2(mua);
1268 : ! d(-vloc)/dR * d(gylm)/d(eps_gamma_delta)
1269 15768450 : do idir=1,3
1270 : prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)=&
1271 : & prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)&
1272 : & -vpsp1_gr(jc,idir)*half&
1273 : & *(gylmgr_jatom(eps_alpha,ic,ilm) * rfgd_jatom(eps_beta ,ic)&
1274 11353284 : & +gylmgr_jatom(eps_beta ,ic,ilm) * rfgd_jatom(eps_alpha,ic))
1275 : ! delta_alphabeta * d(-v_loc/dr) * gylm
1276 15137712 : if (mua<=3) then
1277 : prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)=&
1278 : & prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)&
1279 5676642 : & -vpsp1_gr(jc,idir)*gylm_jatom(ic,ilm)
1280 : end if
1281 : end do ! dir
1282 : end do ! mua
1283 2523834 : do idir=1,3
1284 : ! d(-v_loc/dr) * gylm
1285 : prodp_nondiag(jatm)%value(idir,ilm) = prodp_nondiag(jatm)%value(idir,ilm)&
1286 2522952 : & -vpsp1_gr(jc,idir)*gylm_jatom(ic,ilm)
1287 : end do !end loop idir
1288 : ! END INTERNAL STRAIN CONTRIBUTION
1289 : end do
1290 : end do
1291 : end if ! optstr2
1292 :
1293 : ! Release temp memory allocated for atom j
1294 21178 : if (save_memory.and.jatom/=iatom) then
1295 182 : ABI_FREE(ifftsph_tmp)
1296 182 : ABI_FREE(rfgd_tmp)
1297 182 : ABI_FREE(gylm_jatom)
1298 182 : ABI_FREE(gylmgr_jatom)
1299 182 : if (optgr2==1.and.qne0==1) then
1300 12 : ABI_FREE(expiqr_jatom)
1301 : end if
1302 : end if
1303 :
1304 : end do ! loop on atoms j
1305 : end if ! optgr2 or optstr2
1306 :
1307 : ! --- Apply scaling factor on integrals ---
1308 867127 : if (ngrad >0) prod (:,:)=prod (:,:)*fact_ucvol
1309 26373 : if (ngradp>0) prodp(:,:)=prodp(:,:)*fact_ucvol
1310 20919 : if (ngrad_nondiag>0) then
1311 336 : do jatm=1,natom
1312 76358 : prod_nondiag(jatm)%value(:,:)=prod_nondiag(jatm)%value(:,:)*fact_ucvol
1313 : end do
1314 : end if
1315 20919 : if (ngradp_nondiag>0) then
1316 336 : do jatm=1,natom
1317 14150 : prodp_nondiag(jatm)%value(:,:)=prodp_nondiag(jatm)%value(:,:)*fact_ucvol
1318 : end do
1319 : end if
1320 :
1321 : ! --- Reduction in case of parallelization ---
1322 20919 : if (paral_grid) then
1323 1040 : if (ngrad>0) then
1324 1040 : call xmpi_sum(prod,my_comm_grid,ier)
1325 : end if
1326 1040 : if (ngradp>0) then
1327 0 : call xmpi_sum(prodp,my_comm_grid,ier)
1328 : end if
1329 1040 : if (ngrad_nondiag>0.or.ngradp_nondiag>0) then
1330 0 : bufsiz=0;bufind=0
1331 0 : do jatm=1,natom
1332 0 : jtypat=typat(atindx1(jatm))
1333 0 : bufsiz=bufsiz+pawtab(jtypat)%lcut_size**2
1334 : end do
1335 0 : ABI_MALLOC(buf,(ngrad_nondiag+ngradp_nondiag,bufsiz))
1336 0 : do jatm=1,natom
1337 0 : jtypat=typat(atindx1(jatm))
1338 0 : lm_sizej=pawtab(jtypat)%lcut_size**2
1339 0 : if (ngrad_nondiag> 0) buf(1:ngrad_nondiag,bufind+1:bufind+lm_sizej)= &
1340 0 : & prod_nondiag(jatm)%value(:,:)
1341 0 : if (ngradp_nondiag>0) buf(ngrad_nondiag+1:ngrad_nondiag+ngradp_nondiag, &
1342 0 : & bufind+1:bufind+lm_sizej)=prodp_nondiag(jatm)%value(:,:)
1343 0 : bufind=bufind+lm_sizej*(ngrad_nondiag+ngradp_nondiag)
1344 : end do
1345 0 : call xmpi_sum(buf,my_comm_grid,ier)
1346 0 : bufind=0
1347 0 : do jatm=1,natom
1348 0 : jtypat=typat(atindx1(jatm))
1349 0 : lm_sizej=pawtab(jtypat)%lcut_size**2
1350 0 : if (ngrad> 0) prod_nondiag(jatm)%value(:,:)= &
1351 0 : & buf(1:ngrad_nondiag,bufind+1:bufind+lm_sizej)
1352 0 : if (ngradp>0) prodp_nondiag(jatm)%value(:,:)= &
1353 0 : & buf(ngrad_nondiag+1:ngrad_nondiag+ngradp_nondiag,bufind+1:bufind+lm_sizej)
1354 0 : bufind=bufind+lm_sizej*(ngrad_nondiag+ngradp_nondiag)
1355 : end do
1356 0 : ABI_FREE(buf)
1357 : end if
1358 : end if
1359 :
1360 : ! ----------------------------------------------------------------
1361 : ! Compute final sums (i.e. derivatives of Sum_ij[rho_ij.Intg{Qij.Vloc}]
1362 :
1363 : ! ---- Compute terms common to all gradients
1364 20919 : jrhoij=1
1365 565368 : do irhoij=1,pawrhoij_iatom%nrhoijsel
1366 544449 : klmn=pawrhoij_iatom%rhoijselect(irhoij)
1367 544449 : klm =pawtab(itypat)%indklmn(1,klmn)
1368 544449 : lmin=pawtab(itypat)%indklmn(3,klmn)
1369 544449 : lmax=pawtab(itypat)%indklmn(4,klmn)
1370 544449 : ro =pawrhoij_iatom%rhoijp(jrhoij,ispden)
1371 544449 : ro_d=ro*pawtab(itypat)%dltij(klmn)
1372 1528528 : do ll=lmin,lmax,2
1373 5186573 : do ilm=ll**2+1,(ll+1)**2
1374 3658045 : isel=pawang%gntselect(ilm,klm)
1375 4642124 : if (isel>0) then
1376 963741 : grhat_x=ro_d*pawtab(itypat)%qijl(ilm,klmn)
1377 5356380 : do mu=1,ngrad
1378 5356380 : grhat_tmp(mu,idiag)=grhat_tmp(mu,idiag)+grhat_x*prod(mu,ilm)
1379 : ! DEBUG
1380 : ! if(mu==ishift_str+1 .and. &
1381 : !& (abs(grhat_x*prod(mu,ilm))>tol6 .or. irhoij==1 ) )then
1382 : ! write(6,'(a,5i4,3es16.6)')&
1383 : !& 'mu,idiag,ilm,irhoij,ll, grhat_tmp(mu,idiag),grhat_x,prod(mu,ilm)=',&
1384 : !& mu,idiag,ilm,irhoij,ll, grhat_tmp(mu,idiag),grhat_x,prod(mu,ilm)
1385 : ! endif
1386 : ! ENDDEBUG
1387 : end do
1388 : end if
1389 : end do
1390 : end do
1391 565368 : jrhoij=jrhoij+pawrhoij_iatom%cplex_rhoij
1392 : end do
1393 :
1394 : ! DEBUG
1395 : ! write(6,*)' Accumulation of grhat_tmp : idiag,grhat_tmp(ishift_str+1,idiag),',idiag,grhat_tmp(ishift_str+1,idiag)
1396 : ! ENDDEBUG
1397 :
1398 : ! ---- Add additional (diagonal) terms for dynamical matrix
1399 : ! ---- Terms including rhoij derivatives
1400 20919 : if (optgr2==1) then
1401 : klmn1=1
1402 3915 : do klmn=1,lmn2_size
1403 3838 : klm =pawtab(itypat)%indklmn(1,klmn)
1404 3838 : lmin=pawtab(itypat)%indklmn(3,klmn)
1405 3838 : lmax=pawtab(itypat)%indklmn(4,klmn)
1406 3838 : dlt_tmp=pawtab(itypat)%dltij(klmn)
1407 10567 : do ll=lmin,lmax,2
1408 34810 : do ilm=ll**2+1,(ll+1)**2
1409 24243 : isel=pawang%gntselect(ilm,klm)
1410 30972 : if (isel>0) then
1411 6247 : ro_d= dlt_tmp*pawtab(itypat)%qijl(ilm,klmn)
1412 62470 : do mu=1,9
1413 56223 : mua=alpha(mu);mub=beta(mu)
1414 : grhat_tmp(ishift_gr2+mu,idiag)=grhat_tmp(ishift_gr2+mu,idiag)&
1415 62470 : & +ro_d*pawrhoij_iatom%grhoij(ishift_grhoij+mua,klmn1,ispden)*prodp(mub+ishift_gr,ilm)
1416 : end do
1417 : end if
1418 : end do
1419 : end do
1420 3915 : klmn1=klmn1+pawrhoij_iatom%cplex_rhoij
1421 : end do ! klmn
1422 : end if ! optgr2
1423 :
1424 : ! ---- Add additional (diagonal) terms for elastic tensor
1425 : ! ---- Terms including rhoij derivatives
1426 20919 : if (optstr2==1)then
1427 : klmn1=1
1428 962 : do klmn=1,lmn2_size
1429 936 : klm =pawtab(itypat)%indklmn(1,klmn)
1430 936 : lmin=pawtab(itypat)%indklmn(3,klmn)
1431 936 : lmax=pawtab(itypat)%indklmn(4,klmn)
1432 936 : dlt_tmp=pawtab(itypat)%dltij(klmn)
1433 2418 : do ll=lmin,lmax,2
1434 6708 : do ilm=ll**2+1,(ll+1)**2
1435 4290 : isel=pawang%gntselect(ilm,klm)
1436 5772 : if (isel>0) then
1437 1326 : ro_d=dlt_tmp*pawtab(itypat)%qijl(ilm,klmn)
1438 1326 : mu=1
1439 9282 : do mua=1,6
1440 55692 : do mub=1,6
1441 : grhat_tmp(ishift_str2+mu,iatm)= grhat_tmp(ishift_str2+mu,iatm)&
1442 47736 : & +ro_d*pawrhoij_iatom%grhoij(mub,klmn1,ispden)*prodp(mua,ilm)
1443 : grhat_tmp(ishift_str2+mu,iatm)= grhat_tmp(ishift_str2+mu,iatm)&
1444 47736 : & +ro_d*pawrhoij_iatom%grhoij(mua,klmn1,ispden)*prodp(mub,ilm)
1445 55692 : mu=mu+1
1446 : end do
1447 : ! INTERNAL STRAIN CONTRIBUTION
1448 33150 : do idir=1,3
1449 : grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm) = grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm)&
1450 23868 : & +ro_d*pawrhoij_iatom%grhoij(ishift_grhoij+idir,klmn1,ispden)*prodp(mua,ilm)
1451 : grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm) = grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm)&
1452 31824 : & +ro_d*pawrhoij_iatom%grhoij(mua,klmn1,ispden)*prodp(6+idir,ilm)
1453 : end do
1454 : end do
1455 : end if
1456 : end do
1457 : end do
1458 962 : klmn1=klmn1+pawrhoij_iatom%cplex_rhoij
1459 : end do
1460 : end if ! optstr2
1461 :
1462 : ! ---- Add off-diagonal additional contributions for second gradients
1463 40010 : if (optgr2==1.or.optstr2==1) then
1464 336 : do jatm=1,natom
1465 259 : jatom_tot=atindx1(jatm);jtypat=typat(jatom_tot)
1466 259 : pawrhoij_jatom => pawrhoij_tot(jatom_tot)
1467 :
1468 : ! ---- Dynamical matrix
1469 259 : if (optgr2==1) then
1470 :
1471 : ! Off-diagonal term including rhoij
1472 259 : if (dyfr_cplex==1.or.cplex==1) then
1473 231 : jrhoij=1
1474 5517 : do irhoij=1,pawrhoij_jatom%nrhoijsel
1475 5286 : klmn=pawrhoij_jatom%rhoijselect(irhoij)
1476 5286 : klm =pawtab(jtypat)%indklmn(1,klmn)
1477 5286 : lmin=pawtab(jtypat)%indklmn(3,klmn)
1478 5286 : lmax=pawtab(jtypat)%indklmn(4,klmn)
1479 5286 : ro =pawrhoij_jatom%rhoijp(jrhoij,ispden)
1480 5286 : ro_d=ro*pawtab(jtypat)%dltij(klmn)
1481 13913 : do ll=lmin,lmax,2
1482 39176 : do ilm=ll**2+1,(ll+1)**2
1483 25263 : isel=pawang%gntselect(ilm,klm)
1484 33890 : if (isel>0) then
1485 9101 : grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn)
1486 91010 : do mu=1,9
1487 : grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) &
1488 91010 : & +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)
1489 : end do
1490 : end if
1491 : end do
1492 : end do
1493 5517 : jrhoij=jrhoij+pawrhoij_jatom%cplex_rhoij
1494 : end do
1495 : else
1496 28 : jrhoij=1;mushift=ishift_gr2+9
1497 424 : do irhoij=1,pawrhoij_jatom%nrhoijsel
1498 396 : klmn=pawrhoij_jatom%rhoijselect(irhoij)
1499 396 : klm =pawtab(jtypat)%indklmn(1,klmn)
1500 396 : lmin=pawtab(jtypat)%indklmn(3,klmn)
1501 396 : lmax=pawtab(jtypat)%indklmn(4,klmn)
1502 396 : ro =pawrhoij_jatom%rhoijp(jrhoij,ispden)
1503 396 : ro_d=ro*pawtab(jtypat)%dltij(klmn)
1504 1164 : do ll=lmin,lmax,2
1505 3660 : do ilm=ll**2+1,(ll+1)**2
1506 2496 : isel=pawang%gntselect(ilm,klm)
1507 3264 : if (isel>0) then
1508 1008 : grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn)
1509 10080 : do mu=1,9
1510 : grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm)&
1511 9072 : & +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)
1512 : grhat_tmp(mushift+mu,jatm)=grhat_tmp(mushift+mu,jatm)&
1513 10080 : & +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm)
1514 : end do
1515 : end if
1516 : end do
1517 : end do
1518 424 : jrhoij=jrhoij+pawrhoij_jatom%cplex_rhoij
1519 : end do
1520 : end if
1521 :
1522 : ! Off-diagonal term including rhoij derivative
1523 259 : if (dyfr_cplex==1.or.cplex==1) then
1524 231 : klmn1=1
1525 9599 : do klmn=1,pawrhoij_jatom%lmn2_size
1526 9368 : klm =pawtab(jtypat)%indklmn(1,klmn)
1527 9368 : lmin=pawtab(jtypat)%indklmn(3,klmn)
1528 9368 : lmax=pawtab(jtypat)%indklmn(4,klmn)
1529 9368 : dlt_tmp=pawtab(jtypat)%dltij(klmn)
1530 24775 : do ll=lmin,lmax,2
1531 73766 : do ilm=ll**2+1,(ll+1)**2
1532 48991 : isel=pawang%gntselect(ilm,klm)
1533 64398 : if (isel>0) then
1534 13997 : ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn)
1535 139970 : do mu=1,9
1536 125973 : mua=alpha(mu);mub=beta(mu)
1537 : grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) &
1538 : & +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) &
1539 139970 : & *prodp_nondiag(jatm)%value(ishift2_gr+mub,ilm)
1540 : end do
1541 : end if
1542 : end do
1543 : end do
1544 9599 : klmn1=klmn1+pawrhoij_jatom%cplex_rhoij
1545 : end do ! klmn
1546 : else ! ngradp_nondiag>=6
1547 28 : klmn1=1;mushift=ishift_gr2+9
1548 1576 : do klmn=1,pawrhoij_jatom%lmn2_size
1549 1548 : klm =pawtab(jtypat)%indklmn(1,klmn)
1550 1548 : lmin=pawtab(jtypat)%indklmn(3,klmn)
1551 1548 : lmax=pawtab(jtypat)%indklmn(4,klmn)
1552 1548 : dlt_tmp=pawtab(jtypat)%dltij(klmn)
1553 4364 : do ll=lmin,lmax,2
1554 15084 : do ilm=ll**2+1,(ll+1)**2
1555 10720 : isel=pawang%gntselect(ilm,klm)
1556 13536 : if (isel>0) then
1557 2640 : ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn)
1558 26400 : do mu=1,9
1559 23760 : mua=alpha(mu);mub=beta(mu)
1560 : grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) &
1561 : & +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) &
1562 23760 : & *prodp_nondiag(jatm)%value(ishift2_gr+mub,ilm)
1563 : grhat_tmp(mushift+mu,jatm)=grhat_tmp(mushift+mu,jatm) &
1564 : & +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) &
1565 26400 : & *prodp_nondiag(jatm)%value(ishift2_gr+3+mub,ilm)
1566 : end do
1567 : end if
1568 : end do
1569 : end do
1570 1576 : klmn1=klmn1+pawrhoij_jatom%cplex_rhoij
1571 : end do
1572 : end if
1573 : end if ! optgr2
1574 :
1575 : ! ---- Elastic tensor
1576 21178 : if (optstr2==1)then
1577 :
1578 : ! Off-diagonal term including rhoij
1579 98 : jrhoij=1;
1580 1658 : do irhoij=1,pawrhoij_jatom%nrhoijsel
1581 1560 : klmn=pawrhoij_jatom%rhoijselect(irhoij)
1582 1560 : klm =pawtab(jtypat)%indklmn(1,klmn)
1583 1560 : lmin=pawtab(jtypat)%indklmn(3,klmn)
1584 1560 : lmax=pawtab(jtypat)%indklmn(4,klmn)
1585 1560 : ro =pawrhoij_jatom%rhoijp(jrhoij,ispden)
1586 1560 : ro_d=ro*pawtab(jtypat)%dltij(klmn)
1587 4002 : do ll=lmin,lmax,2
1588 10740 : do ilm=ll**2+1,(ll+1)**2
1589 6738 : isel=pawang%gntselect(ilm,klm)
1590 9180 : if (isel>0) then
1591 3030 : grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn)
1592 57570 : do mu=1,18
1593 : grhat_tmp2(mu,jatm)=grhat_tmp2(mu,jatm) &
1594 57570 : & +grhat_x*prod_nondiag(jatm)%value(ishift_str2is+mu,ilm)
1595 : end do
1596 : end if
1597 : end do
1598 : end do
1599 1658 : jrhoij=jrhoij+pawrhoij_jatom%cplex_rhoij
1600 : end do
1601 : ! Off-diagonal term including rhoij derivative
1602 98 : klmn1=1
1603 3626 : do klmn=1,pawrhoij_jatom%lmn2_size
1604 3528 : klm =pawtab(jtypat)%indklmn(1,klmn)
1605 3528 : lmin=pawtab(jtypat)%indklmn(3,klmn)
1606 3528 : lmax=pawtab(jtypat)%indklmn(4,klmn)
1607 3528 : dlt_tmp=pawtab(jtypat)%dltij(klmn)
1608 9114 : do ll=lmin,lmax,2
1609 25284 : do ilm=ll**2+1,(ll+1)**2
1610 16170 : isel=pawang%gntselect(ilm,klm)
1611 21756 : if (isel>0) then
1612 4998 : ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn)
1613 4998 : mu=1
1614 34986 : do mua=1,6
1615 124950 : do idir=1,3
1616 : grhat_tmp2((mua-1)*3+idir,jatm) = grhat_tmp2((mua-1)*3+idir,jatm) &
1617 : & +ro_d*pawrhoij_jatom%grhoij(mua,klmn1,ispden) &
1618 119952 : & *prodp_nondiag(jatm)%value(idir,ilm)
1619 : end do
1620 : end do
1621 : end if
1622 : end do
1623 : end do
1624 3626 : klmn1=klmn1+pawrhoij_jatom%cplex_rhoij
1625 : end do
1626 : end if ! optstr2
1627 :
1628 : end do ! jatm
1629 : end if ! optgr2 or optstr2
1630 :
1631 : ! ----------------------------------------------------------------
1632 : ! End of loop over spin components
1633 :
1634 : end do ! ispden
1635 :
1636 : ! Eventually free temporary space for g_l(r).Y_lm(r) factors
1637 19091 : if (pawfgrtab_iatom%gylm_allocated==2) then
1638 0 : ABI_FREE(pawfgrtab_iatom%gylm)
1639 0 : ABI_MALLOC(pawfgrtab_iatom%gylm,(0,0))
1640 0 : pawfgrtab_iatom%gylm_allocated=0
1641 : end if
1642 19091 : if (pawfgrtab_iatom%gylmgr_allocated==2) then
1643 1 : ABI_FREE(pawfgrtab_iatom%gylmgr)
1644 1 : ABI_MALLOC(pawfgrtab_iatom%gylmgr,(0,0,0))
1645 1 : pawfgrtab_iatom%gylmgr_allocated=0
1646 : end if
1647 19091 : if (pawfgrtab_iatom%gylmgr2_allocated==2) then
1648 0 : ABI_FREE(pawfgrtab_iatom%gylmgr2)
1649 0 : ABI_MALLOC(pawfgrtab_iatom%gylmgr2,(0,0,0))
1650 0 : pawfgrtab_iatom%gylmgr2_allocated=0
1651 : end if
1652 19091 : if (pawfgrtab_iatom%expiqr_allocated==2) then
1653 0 : ABI_FREE(pawfgrtab_iatom%expiqr)
1654 0 : ABI_MALLOC(pawfgrtab_iatom%expiqr,(0,0))
1655 0 : pawfgrtab_iatom%expiqr_allocated=0
1656 : end if
1657 :
1658 : ! ----------------------------------------------------------------
1659 : ! Copy results in corresponding arrays
1660 :
1661 : ! ==== Forces ====
1662 : ! Convert from cartesian to reduced coordinates
1663 19091 : if (optgr==1) then
1664 17601 : mushift=3*(iatm-1)
1665 70404 : tmp(1:3)=grhat_tmp(ishift_gr+1:ishift_gr+3,idiag)
1666 70404 : do mu=1,3
1667 70404 : hatgr(mu+mushift)=rprimd(1,mu)*tmp(1)+rprimd(2,mu)*tmp(2)+rprimd(3,mu)*tmp(3)
1668 : end do
1669 : end if
1670 :
1671 : ! ==== Stresses ====
1672 19091 : if (optstr==1) then
1673 : ! This is contribution Eq.(41) of Torrent2008.
1674 : !DEBUG
1675 : ! write(6,*)' after loop on ispden,ishift_str,idiag,hatstr(1),grhat_tmp(ishift_str+1)',hatstr(1),grhat_tmp(ishift_str+1,idiag)
1676 : !ENDDEBUG
1677 18487 : hatstr(1:6)=hatstr(1:6)+grhat_tmp(ishift_str+1:ishift_str+6,idiag)
1678 : end if
1679 :
1680 : ! ==== Frozen wf part of dyn. matrix ====
1681 19091 : if (optgr2==1) then
1682 320 : do jatm=1,natom
1683 2490 : do mu=1,9
1684 2241 : mua=alpha(mu);mub=beta(mu)
1685 2490 : dyfr(1,mub,mua,jatm,iatm)=grhat_tmp(ishift_gr2+mu,jatm)
1686 : end do
1687 320 : if (dyfr_cplex==2.and.cplex==2) then
1688 26 : mushift=ishift_gr2+9
1689 260 : do mu=1,9
1690 234 : mua=alpha(mu);mub=beta(mu)
1691 260 : dyfr(2,mub,mua,jatm,iatm)=grhat_tmp(mushift+mu,jatm)
1692 : end do
1693 : end if
1694 : end do
1695 : end if
1696 :
1697 : ! ==== Elastic tensor ====
1698 19091 : if (optstr2==1) then
1699 : eltfr(1:6,1:6)=eltfr(1:6,1:6)+ &
1700 1118 : & reshape(grhat_tmp(ishift_str2+1:ishift_str2+36,iatm),(/6,6/))
1701 : ! Convert internal Strain in reduced coordinates
1702 182 : do mua = 1,6
1703 624 : tmp(1:3)=grhat_tmp(ishift_str2is+(mua-1)*3+1:ishift_str2is+(mua-1)*3+3,iatm)
1704 624 : do idir=1,3
1705 : eltfr(6+(iatm-1)*3+idir,mua)=eltfr(6+(iatm-1)*3+idir,mua)+ &
1706 624 : & (rprimd(1,idir)*tmp(1)+rprimd(2,idir)*tmp(2)+rprimd(3,idir)*tmp(3))
1707 : end do
1708 770 : do jatm=1,natom
1709 2352 : tmp(1:3)=grhat_tmp2((mua-1)*3+1:(mua-1)*3+3,jatm)
1710 2508 : do idir=1,3
1711 : eltfr(6+(iatm-1)*3+idir,mua)=eltfr(6+(iatm-1)*3+idir,mua)+ &
1712 2352 : & (rprimd(1,idir)*tmp(1)+rprimd(2,idir)*tmp(2)+rprimd(3,idir)*tmp(3))
1713 : end do
1714 : end do
1715 : end do
1716 : end if
1717 :
1718 : ! ----------------------------------------------------------------
1719 : ! End loops on types and atoms
1720 :
1721 19091 : ABI_FREE(vloc)
1722 19091 : if (ngrad>0) then
1723 19071 : ABI_FREE(prod)
1724 : end if
1725 19091 : if (ngradp>0) then
1726 71 : ABI_FREE(prodp)
1727 : end if
1728 30526 : if (optgr2==1.or.optstr2==1) then
1729 320 : do jatm=1,natom
1730 249 : ABI_FREE(prod_nondiag(jatm)%value)
1731 23617 : ABI_FREE(prodp_nondiag(jatm)%value)
1732 : end do
1733 : end if
1734 : end do ! iatm
1735 18858 : iatshft=iatshft+nattyp(itypat)
1736 : end do ! itypat
1737 :
1738 : !DEBUG
1739 : ! write(6,*)' before parallelization over atoms hatstr(1)',hatstr(1)
1740 : !ENDDEBUG
1741 :
1742 : !Reduction in case of parallelisation over atoms
1743 7423 : if (paral_atom) then
1744 3072 : bufsiz=3*natom*optgr+6*optstr
1745 3072 : if (save_memory) bufsiz=bufsiz+9*dyfr_cplex*natom**2*optgr2+6*(6+3*natom)*optstr2
1746 3072 : if (bufsiz>0) then
1747 9204 : ABI_MALLOC(buf1,(bufsiz))
1748 20108 : if (optgr==1) buf1(1:3*natom)=hatgr(1:3*natom)
1749 3068 : indx=optgr*3*natom
1750 6272 : if (optstr==1) buf1(indx+1:indx+6)=hatstr(1:6)
1751 3068 : indx=indx+optstr*6
1752 : if (save_memory) then
1753 3068 : if (optgr2==1) then
1754 : buf1(indx+1:indx+9*dyfr_cplex*natom**2)= &
1755 44 : & reshape(dyfr,(/9*dyfr_cplex*natom**2/))
1756 22 : indx=indx+9*dyfr_cplex*natom**2
1757 : end if
1758 3068 : if (optstr2==1) then
1759 : buf1(indx+1:indx+6*(6+3*natom))= &
1760 12 : & reshape(eltfr,(/6*(6+3*natom)/))
1761 6 : indx=indx+6*(6+3*natom)
1762 : end if
1763 : end if
1764 3068 : call xmpi_sum(buf1,my_comm_atom,ier)
1765 20108 : if (optgr==1) hatgr(1:3*natom)=buf1(1:3*natom)
1766 3068 : indx=optgr*3*natom
1767 6272 : if (optstr==1) hatstr(1:6)=buf1(indx+1:indx+6)
1768 3068 : indx=indx+optstr*6
1769 : if (save_memory) then
1770 3068 : if (optgr2==1) then
1771 : dyfr(1:dyfr_cplex,1:3,1:3,1:natom,1:natom)= &
1772 10800 : & reshape(buf1(indx+1:indx+9*dyfr_cplex*natom**2),(/dyfr_cplex,3,3,natom,natom/))
1773 22 : indx=indx+9*dyfr_cplex*natom**2
1774 : end if
1775 3068 : if (optstr2==1) then
1776 : eltfr(1:6+3*natom,1:6)= &
1777 702 : & reshape(buf1(indx+1:indx+6*(6+3*natom)),(/6+3*natom,6/))
1778 6 : indx=indx+6*(6+3*natom)
1779 : end if
1780 : end if
1781 6136 : ABI_FREE(buf1)
1782 : end if
1783 : end if
1784 :
1785 : !Deallocate additional memory
1786 7423 : ABI_FREE(grhat_tmp)
1787 7423 : if (optgr2==1.or.optstr2==1) then
1788 44 : ABI_FREE(mu4)
1789 44 : ABI_FREE(atindx)
1790 44 : if (optgr2==1.or.optstr2==1) then
1791 44 : ABI_FREE(vpsp1_gr)
1792 : end if
1793 44 : if (optstr2==1) then
1794 12 : ABI_FREE(grhat_tmp2)
1795 12 : ABI_FREE(vpsp1_str)
1796 : end if
1797 176 : ABI_FREE(prod_nondiag)
1798 176 : ABI_FREE(prodp_nondiag)
1799 : if (.not.save_memory) then
1800 : do jatom=1,size(pawfgrtab)
1801 : pawfgrtab_jatom => pawfgrtab(jatom)
1802 : if (pawfgrtab(jatom)%gylm_allocated==2) then
1803 : ABI_FREE(pawfgrtab(jatom)%gylm)
1804 : ABI_MALLOC(pawfgrtab(jatom)%gylm,(0,0))
1805 : pawfgrtab(jatom)%gylm_allocated=0
1806 : end if
1807 : if (pawfgrtab(jatom)%gylmgr_allocated==2) then
1808 : ABI_FREE(pawfgrtab(jatom)%gylmgr)
1809 : ABI_MALLOC(pawfgrtab(jatom)%gylmgr,(0,0,0))
1810 : pawfgrtab(jatom)%gylmgr_allocated=0
1811 : end if
1812 : if (pawfgrtab(jatom)%gylmgr2_allocated==2) then
1813 : ABI_FREE(pawfgrtab(jatom)%gylmgr2)
1814 : ABI_MALLOC(pawfgrtab(jatom)%gylmgr2,(0,0,0))
1815 : pawfgrtab(jatom)%gylmgr2_allocated=0
1816 : end if
1817 : if (pawfgrtab(jatom)%expiqr_allocated==2) then
1818 : ABI_FREE(pawfgrtab(jatom)%expiqr)
1819 : ABI_MALLOC(pawfgrtab(jatom)%expiqr,(0,0))
1820 : pawfgrtab(jatom)%expiqr_allocated=0
1821 : end if
1822 : end do
1823 : end if
1824 44 : if (paral_atom) then
1825 : if ((.not.save_memory).and.paral_atom_pawfgrtab) then
1826 : call pawfgrtab_free(pawfgrtab_tot)
1827 : ABI_FREE(pawfgrtab_tot)
1828 : end if
1829 22 : if (paral_atom_pawrhoij) then
1830 0 : call pawrhoij_free(pawrhoij_tot)
1831 0 : ABI_FREE(pawrhoij_tot)
1832 : end if
1833 : end if
1834 : end if
1835 :
1836 : !----------------------------------------------------------------------
1837 : !Update non-local gradients
1838 :
1839 : !===== Update forces =====
1840 7423 : if (optgr==1) then
1841 71029 : grnl(1:3*natom)=grnl(1:3*natom)+hatgr(1:3*natom)
1842 6694 : ABI_FREE(hatgr)
1843 : end if
1844 :
1845 : !===== Convert stresses (add diag and off-diag contributions) =====
1846 7423 : if (optstr==1) then
1847 :
1848 : ! Has to compute int[nhat*vtrial]. See Eq.(40) in Torrent2008 .
1849 1395 : hatstr_diag=zero
1850 1395 : if (nspden==1.or.dimvtrial==1) then
1851 45018755 : do ic=1,nfft
1852 45018755 : hatstr_diag=hatstr_diag+vtrial_(ic,1)*nhat(ic,1)
1853 : end do
1854 158 : else if (nspden==2) then
1855 4271268 : do ic=1,nfft
1856 4271268 : hatstr_diag=hatstr_diag+vtrial_(ic,1)*nhat(ic,2)+vtrial_(ic,2)*(nhat(ic,1)-nhat(ic,2))
1857 : end do
1858 38 : else if (nspden==4) then
1859 484381 : do ic=1,nfft
1860 : hatstr_diag=hatstr_diag+half*(vtrial_(ic,1)*(nhat(ic,1)+nhat(ic,4)) &
1861 : & +vtrial_(ic,2)*(nhat(ic,1)-nhat(ic,4))) &
1862 484381 : & +vtrial_(ic,3)*nhat(ic,2)-vtrial_(ic,4)*nhat(ic,3)
1863 : end do
1864 : end if
1865 1395 : hatstr_diag=hatstr_diag*fact_ucvol
1866 1395 : if (paral_grid) then
1867 236 : call xmpi_sum(hatstr_diag,my_comm_grid,ier)
1868 : end if
1869 :
1870 : ! Convert hat contribution
1871 :
1872 : !DEBUG
1873 : ! write(6,*)' hatstr(1),hatstr_diag,nlstr(1)=',hatstr(1),hatstr_diag,nlstr(1)
1874 : !ENDDEBUG
1875 :
1876 5580 : hatstr(1:3)=(hatstr(1:3)+hatstr_diag)/ucvol
1877 5580 : hatstr(4:6)= hatstr(4:6)/ucvol
1878 :
1879 : ! Add to already computed NL contrib
1880 9765 : nlstr(1:6)=nlstr(1:6)+hatstr(1:6)
1881 :
1882 : ! Apply symmetries
1883 1395 : call stresssym(gprimd,nsym,nlstr,symrec)
1884 : end if
1885 :
1886 : !===== Convert dynamical matrix (from cartesian to reduced coordinates) =====
1887 7423 : if (optgr2==1) then
1888 176 : do iatm=1,natom
1889 768 : do jatm=1,natom
1890 2368 : do mua=1,3
1891 7696 : do mub=1,3
1892 7104 : work1(1,mua,mub)=dyfr(1,mub,mua,jatm,iatm)+dyfr(1,mua,mub,iatm,jatm)
1893 : end do
1894 : end do
1895 592 : if (dyfr_cplex==2) then
1896 200 : do mua=1,3
1897 650 : do mub=1,3
1898 600 : work1(2,mua,mub)=dyfr(2,mub,mua,jatm,iatm)-dyfr(2,mua,mub,iatm,jatm)
1899 : end do
1900 : end do
1901 : end if
1902 2368 : do mu=1,3
1903 13474 : work2(:,:,mu)=rprimd(1,mu)*work1(:,:,1)+rprimd(2,mu)*work1(:,:,2)+rprimd(3,mu)*work1(:,:,3)
1904 : end do
1905 2500 : do mub=1,3
1906 7696 : do mua=1,3
1907 : dyfrnl(:,mua,mub,jatm,iatm)=dyfrnl(:,mua,mub,jatm,iatm) & ! Already contains NL projectors contribution
1908 : & +rprimd(1,mua)*work2(:,1,mub) &
1909 : & +rprimd(2,mua)*work2(:,2,mub) &
1910 12882 : & +rprimd(3,mua)*work2(:,3,mub)
1911 : end do
1912 : end do
1913 : end do
1914 : end do
1915 44 : ABI_FREE(dyfr)
1916 : end if
1917 :
1918 : !===== Update elastic tensor =====
1919 7423 : if (optstr2==1) then
1920 1272 : eltfrnl(1:6+3*natom,1:6)=eltfrnl(1:6+3*natom,1:6)+eltfr(1:6+3*natom,1:6)
1921 12 : ABI_FREE(eltfr)
1922 : end if
1923 :
1924 : !----------------------------------------------------------------------
1925 : !End
1926 :
1927 : !Destroy temporary space
1928 7423 : if (usexcnhat==0) then
1929 5085 : ABI_FREE(vtrial_)
1930 : end if
1931 :
1932 : !Destroy atom tables used for parallelism
1933 7423 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1934 7423 : if (paral_atom) then
1935 3072 : ABI_FREE(atm_indx)
1936 : end if
1937 :
1938 : !Destroy FFT tables used for parallelism
1939 14890 : if ((optgr2==1.or.optstr2==1).and.(.not.present(comm_fft))) then
1940 44 : call my_distribfft%free()
1941 44 : ABI_FREE(my_distribfft)
1942 : end if
1943 :
1944 : DBG_ENTER("COLL")
1945 :
1946 : CONTAINS
1947 : !!***
1948 :
1949 : ! ------------------------------------------------
1950 : !!****f* pawgrnl/pawgrnl_convert
1951 : !! NAME
1952 : !! pawgrnl_convert
1953 : !!
1954 : !! FUNCTION
1955 : !! notation: Convert index of the elastic tensor:
1956 : !! - voigt notation => 32
1957 : !! - normal notation => 3 3 2 2
1958 : !! - notation for gylmgr2 => 32 32 32 32 => 4 4 4
1959 : !!
1960 : !! INPUTS
1961 : !! eps_alpha, eps_beta, eps_delta, eps_gamma
1962 : !!
1963 : !! OUTPUT
1964 : !! mu4(4) = array with index for the second derivative of gylm
1965 : !!
1966 : !! SIDE EFFECTS
1967 : !! mu4(4) = input : array with index for the second derivative of gylm
1968 : !! output: the 4 indexes for the calculation of the second derivative of gylm
1969 : !!
1970 : !! SOURCE
1971 :
1972 11664972 : subroutine pawgrnl_convert(mu4,eps_alpha,eps_beta,eps_gamma,eps_delta)
1973 :
1974 : !Arguments ------------------------------------
1975 : !scalar
1976 : integer,intent(in) :: eps_alpha,eps_beta
1977 : integer,optional,intent(in) :: eps_gamma,eps_delta
1978 : !array
1979 : integer,intent(inout) :: mu4(4)
1980 :
1981 : !Local variables-------------------------------
1982 : integer :: eps1,eps2,i,j,k
1983 11664972 : integer,allocatable :: mu_temp(:)
1984 :
1985 : ! *************************************************************************
1986 :
1987 11664972 : ABI_MALLOC(mu_temp,(4))
1988 11664972 : if (present(eps_gamma).and.present(eps_delta)) then
1989 11664972 : mu_temp(1)=eps_alpha
1990 11664972 : mu_temp(2)=eps_beta
1991 11664972 : mu_temp(3)=eps_gamma
1992 11664972 : mu_temp(4)=eps_delta
1993 : else
1994 0 : mu_temp(1)=eps_alpha
1995 0 : mu_temp(2)=eps_beta
1996 0 : mu_temp(3)= 0
1997 0 : mu_temp(4)= 0
1998 : end if
1999 : k=1
2000 34994916 : do i=1,2
2001 23329944 : eps1=mu_temp(i)
2002 81654804 : do j=1,2
2003 46659888 : eps2=mu_temp(2+j)
2004 46659888 : if(eps1==eps2) then
2005 15553296 : if(eps1==1) mu4(k)=1;
2006 15553296 : if(eps1==2) mu4(k)=2;
2007 15553296 : if(eps1==3) mu4(k)=3;
2008 : else
2009 31106592 : if((eps1==3.and.eps2==2).or.(eps1==2.and.eps2==3)) mu4(k)=4;
2010 31106592 : if((eps1==3.and.eps2==1).or.(eps1==1.and.eps2==3)) mu4(k)=5;
2011 31106592 : if((eps1==1.and.eps2==2).or.(eps1==2.and.eps2==1)) mu4(k)=6;
2012 : end if
2013 69989832 : k=k+1
2014 : end do
2015 : end do
2016 11664972 : ABI_FREE(mu_temp)
2017 :
2018 11664972 : end subroutine pawgrnl_convert
2019 : ! ------------------------------------------------
2020 :
2021 : end subroutine pawgrnl
2022 : !!***
2023 :
2024 : !----------------------------------------------------------------------
2025 :
2026 : !!****f* m_paw_dfpt/dsdr_k_paw
2027 : !! NAME
2028 : !! dsdr_k_paw
2029 : !!
2030 : !! FUNCTION
2031 : !! compute on-site terms for forces and stresses for finite electric fields with PAW
2032 : !!
2033 : !! INPUTS
2034 : !! cprj_k (pawcprj_type) :: cprj for occupied bands at point k
2035 : !! cprj_kb :: cprj for occupied bands at point k+b
2036 : !! dtefield :: structure referring to all efield and berry's phase variables
2037 : !! kdir :: integer giving direction along which overlap is computed for ket
2038 : !! kfor :: integer indicating whether to compute forward (1) or backward (2)
2039 : !! along kpt string
2040 : !! natom :: number of atoms in cell
2041 : !! typat :: typat(natom) type of each atom
2042 : !!
2043 : !! OUTPUT
2044 : !!
2045 : !! SIDE EFFECTS
2046 : !! dsdr :: array of the on-site PAW parts of the derivatives with respect to atm
2047 : !! positions and/or strains of the overlaps between Bloch states at points
2048 : !! k and k+b, for the various pairs of bands
2049 : !!
2050 : !! NOTES
2051 : !! This routine assumes that the cprj are not explicitly ordered by
2052 : !! atom type.
2053 : !!
2054 : !! SOURCE
2055 :
2056 149568 : subroutine dsdr_k_paw(cprj_k,cprj_kb,dsdr,dtefield,kdir,kfor,mband,natom,ncpgr,typat)
2057 :
2058 : !Arguments---------------------------
2059 : !scalars
2060 : integer,intent(in) :: kdir,kfor,mband,natom,ncpgr
2061 : character(len=500) :: message
2062 : type(efield_type),intent(in) :: dtefield
2063 : type(pawcprj_type),intent(in) :: cprj_k(natom,dtefield%nspinor*mband)
2064 : type(pawcprj_type),intent(in) :: cprj_kb(natom,dtefield%nspinor*mband)
2065 :
2066 : !arrays
2067 : integer,intent(in) :: typat(natom)
2068 : real(dp),intent(inout) :: dsdr(2,natom,ncpgr,dtefield%mband_occ,dtefield%mband_occ)
2069 :
2070 : !Local variables---------------------------
2071 : !scalars
2072 : integer :: iatom,iband,ibs,icpgr,ilmn,ispinor,itypat
2073 : integer :: jband,jbs,jlmn,klmn,nspinor
2074 : complex(dp) :: cpk,cpkb,dcpk,dcpkb,cterm,paw_onsite
2075 : ! *************************************************************************
2076 :
2077 : !initialize dsdr
2078 153905472 : dsdr(:,:,:,:,:) = zero
2079 :
2080 : ! if 3 gradients we are in the ctocprj choice 2 case
2081 : ! and the 3 gradients are due to the atomic displacements
2082 : ! if 6 gradients we are in the ctocprj choice 3 case
2083 : ! and the 6 gradients are due to the strains
2084 : ! if 9 gradients we are in the ctocprj choice 23 case
2085 : ! and the first six are due to strain, last three due to displacements
2086 149568 : if (ncpgr /= 3 .and. ncpgr /= 6 .and. ncpgr /= 9) then
2087 0 : message = ' dsdr_k_paw called with ncpgr /= 3, 6, or 9 (no gradients) '
2088 0 : ABI_BUG(message)
2089 : end if
2090 :
2091 : nspinor = dtefield%nspinor
2092 :
2093 448704 : do iatom = 1, natom
2094 299136 : itypat = typat(iatom)
2095 :
2096 2841792 : do ilmn=1,dtefield%lmn_size(itypat)
2097 21836928 : do jlmn=1,dtefield%lmn_size(itypat)
2098 19144704 : klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn)
2099 : paw_onsite = cmplx(dtefield%qijb_kk(1,klmn,iatom,kdir),&
2100 19144704 : & dtefield%qijb_kk(2,klmn,iatom,kdir))
2101 19144704 : if (kfor > 1) paw_onsite = conjg(paw_onsite)
2102 98116608 : do iband = 1, dtefield%mband_occ
2103 402038784 : do jband = 1, dtefield%mband_occ
2104 689209344 : do ispinor = 1, nspinor
2105 3369467904 : do icpgr = 1, ncpgr
2106 2756837376 : ibs = nspinor*(iband-1) + ispinor
2107 2756837376 : jbs = nspinor*(jband-1) + ispinor
2108 2756837376 : cpk=cmplx(cprj_k(iatom,ibs)%cp(1,ilmn),cprj_k(iatom,ibs)%cp(2,ilmn))
2109 2756837376 : dcpk=cmplx(cprj_k(iatom,ibs)%dcp(1,icpgr,ilmn),cprj_k(iatom,ibs)%dcp(2,icpgr,ilmn))
2110 2756837376 : cpkb=cmplx(cprj_kb(iatom,jbs)%cp(1,jlmn),cprj_kb(iatom,jbs)%cp(2,jlmn))
2111 2756837376 : dcpkb=cmplx(cprj_kb(iatom,jbs)%dcp(1,icpgr,jlmn),cprj_kb(iatom,jbs)%dcp(2,icpgr,jlmn))
2112 2756837376 : cterm=paw_onsite*(conjg(dcpk)*cpkb+conjg(cpk)*dcpkb)
2113 2756837376 : dsdr(1,iatom,icpgr,iband,jband) = dsdr(1,iatom,icpgr,iband,jband)+real(cterm)
2114 3063152640 : dsdr(2,iatom,icpgr,iband,jband) = dsdr(2,iatom,icpgr,iband,jband)+aimag(cterm)
2115 : end do ! end loop over icpgr
2116 : end do ! end loop over ispinor
2117 : end do ! end loop over jband
2118 : end do ! end loop over iband
2119 : end do ! end loop over ilmn
2120 : end do ! end loop over jlmn
2121 :
2122 : end do ! end loop over atoms
2123 :
2124 149568 : end subroutine dsdr_k_paw
2125 : !!***
2126 :
2127 : !----------------------------------------------------------------------
2128 :
2129 : END MODULE m_paw_dfpt
2130 : !!***
|
