Branch data 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 : 15370 : subroutine pawgrnl(atindx1,dimnhat,dyfrnl,dyfr_cplex,eltfrnl,grnl,gsqcut,mgfft,my_natom,natom,&
482 : 7685 : & nattyp,nfft,ngfft,nhat,nlstr,nspden,nsym,ntypat,optgr,optgr2,optstr,optstr2,&
483 [ - + - + ]: 7685 : & pawang,pawfgrtab,pawrhoij,pawtab,ph1d,psps,qphon,rprimd,symrec,typat,ucvol,vtrial,vxc,xred,&
484 [ + - + + : 7685 : & 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 : 7685 : integer,allocatable :: atindx(:),atm_indx(:),mu4(:)
532 : 7685 : integer,allocatable,target :: ifftsph_tmp(:)
533 : 7685 : integer,ABI_CONTIGUOUS pointer :: ffti3_local(:),fftn3_distrib(:),ifft_jatom(:)
534 : 7685 : 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 : 15370 : real(dp) :: work1(dyfr_cplex,3,3),work2(dyfr_cplex,3,3)
537 : 7685 : real(dp),allocatable :: buf(:,:),buf1(:),dyfr(:,:,:,:,:),eltfr(:,:)
538 : 7685 : real(dp),allocatable :: grhat_tmp(:,:),grhat_tmp2(:,:),hatgr(:)
539 : 7685 : real(dp),allocatable :: prod(:,:),prodp(:,:),vloc(:),vpsp1_gr(:,:),vpsp1_str(:,:)
540 : 7685 : real(dp),allocatable,target :: rfgd_tmp(:,:)
541 : 7685 : real(dp),ABI_CONTIGUOUS pointer :: gylm_jatom(:,:),gylmgr_jatom(:,:,:),gylmgr2_jatom(:,:,:),expiqr_jatom(:,:)
542 : 7685 : real(dp),ABI_CONTIGUOUS pointer :: rfgd_jatom(:,:),vtrial_(:,:)
543 : 7685 : type(coeff2_type),allocatable :: prod_nondiag(:),prodp_nondiag(:)
544 : 7685 : type(pawfgrtab_type),pointer :: pawfgrtab_(:),pawfgrtab_tot(:)
545 : 7685 : 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 [ + + ]: 7685 : if (my_natom>0) then
554 [ + + - + : 6641 : 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 [ + + ]: 6641 : 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 [ - + ]: 6641 : 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 [ + - + + ]: 7685 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
573 : 7685 : paral_atom_pawfgrtab=(size(pawfgrtab)/=natom)
574 : 7685 : paral_atom_pawrhoij=(size(pawrhoij)/=natom)
575 [ + + ]: 7685 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
576 [ + - ]: 7685 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
577 : 7685 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
578 [ + + ]: 7685 : 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 : 7685 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nfftot=n1*n2*n3
588 : 7685 : my_comm_grid=xmpi_comm_self;my_comm_fft=xmpi_comm_self;me_fft=0
589 : 7685 : my_me_g0=1;my_paral_kgb=0;paral_grid=.false.;nullify(my_distribfft)
590 [ + + ]: 7685 : if (present(mpi_comm_grid).or.present(comm_fft)) then
591 [ + - ]: 7628 : if (present(mpi_comm_grid)) my_comm_grid=mpi_comm_grid
592 [ + + ]: 7628 : if (present(comm_fft)) my_comm_fft=comm_fft
593 [ - + ]: 7628 : if (.not.present(mpi_comm_grid)) my_comm_grid=comm_fft
594 [ + + ]: 7628 : if (.not.present(comm_fft)) my_comm_fft=mpi_comm_grid
595 : 7628 : paral_grid=(xmpi_comm_size(my_comm_grid)>1)
596 : 7628 : me_fft=xmpi_comm_rank(my_comm_fft)
597 : : end if
598 [ + + - + ]: 7685 : 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 : 7685 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol_)
626 : 7685 : fact_ucvol=ucvol/dble(nfftot)
627 : :
628 : : !Retrieve local potential according to the use of nhat in XC
629 [ + + ]: 19644 : usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
630 [ + + ]: 7685 : if (usexcnhat==0) then
631 : : ! dimvtrial=nspden
632 : 5347 : dimvtrial=1
633 [ - + + - : 16041 : ABI_MALLOC(vtrial_,(nfft,dimvtrial))
+ - + - -
+ - + ]
634 : : !!!$OMP PARALLEL DO PRIVATE(ic) SHARED(nfft,vtrial,vtrial_,vxc)
635 [ + + ]: 135007552 : do ic=1,nfft
636 [ + + + + ]: 540014167 : 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 : 7685 : ngrhat=0;ngrad=0;ngradp=0;ngrad_nondiag=0;ngradp_nondiag=0
645 : 7685 : ishift_grhoij=0;ishift_gr=0;ishift_gr2=0;ishift_str=0;ishift_str2=0;ishift_str2is=0;ishift2_gr=0
646 [ + + ]: 7685 : cplex=1;if (qne0==1) cplex=2
647 [ + + ]: 7685 : if (optgr==1) then
648 [ + - - + : 27756 : ABI_MALLOC(hatgr,(3*natom))
+ - + - +
- - + -
+ ]
649 [ + + ]: 72744 : hatgr=zero
650 : : ngrad=ngrad+3
651 : : ngrhat=ngrhat+3
652 : : ishift_gr2=ishift_gr2+3
653 : : end if
654 [ + + ]: 7685 : 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 [ + + ]: 7685 : if (optstr==1) then
663 : 1412 : hatstr=zero
664 : 1412 : ngrad=ngrad+6
665 : 1412 : ngrhat=ngrhat+6
666 : 1412 : ishift_gr=ishift_gr+6
667 : 1412 : ishift_gr2=ishift_gr2+6
668 : 1412 : ishift_str2=ishift_str2+6
669 : 1412 : ishift_str2is = ishift_str2is+6
670 : : end if
671 [ + + ]: 7685 : 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 [ + + + + ]: 7685 : nsploop=nspden;if (dimvtrial<nspden .and. nspden==4) nsploop=1
689 [ + + + - ]: 7685 : if (optgr2/=1.and.optstr2/=1) then
690 [ + + + - : 22923 : 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 : 7685 : pawfgrtab_tot => pawfgrtab
770 : : end if
771 [ + + - + : 7685 : 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 : 7685 : pawrhoij_tot => pawrhoij
778 : : end if
779 : :
780 : : if (save_memory) then
781 : 7685 : pawfgrtab_ => pawfgrtab
782 : 7685 : 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 : 7685 : iatshft=0
792 [ + + ]: 19644 : do itypat=1,ntypat
793 : :
794 : 11959 : lmn2_size=pawtab(itypat)%lmn2_size
795 : 11959 : lm_size=pawtab(itypat)%lcut_size**2
796 : :
797 [ + + ]: 35851 : do iatm=iatshft+1,iatshft+nattyp(itypat)
798 : :
799 : 23892 : iatom_tot=atindx1(iatm)
800 : 23892 : iatom=iatom_tot
801 [ + + ]: 23892 : 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 [ + + ]: 23892 : if (iatom==-1) cycle
806 [ + + ]: 19615 : iatom_pawfgrtab=iatom_tot;if (paral_atom_pawfgrtab) iatom_pawfgrtab=iatom
807 [ + + ]: 19615 : iatom_pawrhoij =iatom_tot;if (paral_atom_pawrhoij) iatom_pawrhoij =iatom
808 : 19615 : pawfgrtab_iatom => pawfgrtab_(iatom_pawfgrtab)
809 : 19615 : pawrhoij_iatom => pawrhoij_(iatom_pawrhoij)
810 : :
811 [ + + - + ]: 19615 : idiag=1;if (optgr2==1.or.optstr2==1) idiag=iatm
812 : 19615 : nfgd=pawfgrtab_iatom%nfgd
813 : :
814 [ + + + - : 58845 : ABI_MALLOC(vloc,(nfgd))
+ - + + -
+ - + ]
815 [ + + ]: 19615 : if (ngrad>0) then
816 [ - + + - : 78380 : ABI_MALLOC(prod,(ngrad,lm_size))
- + + - +
- + - - +
- + - + ]
817 : : end if
818 [ + + ]: 19615 : if (ngradp>0) then
819 [ - + + - : 284 : ABI_MALLOC(prodp,(ngradp,lm_size))
- + + - +
- + - - +
- + - + ]
820 : : end if
821 [ + + ]: 19615 : 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 [ + + + + ]: 117498 : grhat_tmp=zero
833 [ + + + + : 21503 : 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 [ + + + + : 19615 : if ((optgr==1.or.optstr==1).and.(optgr2/=1).and.(optstr2/=1)) then
+ - + - ]
840 [ + + ]: 19524 : 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 [ + + - + ]: 19615 : 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 [ + + + + : 19615 : 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 [ + + + + ]: 19615 : has_phase=(optgr2==1.and.pawfgrtab_iatom%expiqr_allocated/=0)
892 : :
893 : : ! Eventually compute 1st-order potential
894 [ + + - + ]: 19615 : 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 [ + + ]: 729766 : 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 [ + + ]: 19615 : 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 [ + + ]: 41058 : do ispden=1,nsploop
930 : :
931 : : ! ----- Retrieve potential (subtle if nspden=4 ;-)
932 [ + + ]: 21443 : if (nspden/=4) then
933 : 20558 : ispvtr=min(dimvtrial,ispden)
934 [ + + ]: 19867262 : do ic=1,nfgd
935 : 19846704 : jc = pawfgrtab_iatom%ifftsph(ic)
936 : 19867262 : 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 [ + + + + : 887433 : if (ngrad>0) prod=zero
+ + ]
971 [ + + + + : 26897 : if (ngradp>0) prodp=zero
+ + ]
972 : :
973 : : ! ==== Contribution to forces ====
974 [ + + ]: 21443 : if (optgr==1) then
975 [ + + ]: 166870 : do ilm=1,lm_size
976 [ + + ]: 235194184 : do ic=1,pawfgrtab_iatom%nfgd
977 [ + + ]: 940256446 : do mu=1,3
978 : : prod(mu+ishift_gr,ilm)=prod(mu+ishift_gr,ilm)-&
979 : 940109256 : & 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 [ + + ]: 21443 : if (optstr==1) then
987 [ + + ]: 43700 : do ilm=1,lm_size
988 [ + + ]: 144545233 : do ic=1,pawfgrtab_iatom%nfgd
989 : 144501533 : jc=pawfgrtab_iatom%ifftsph(ic)
990 [ + + ]: 1011551152 : do mu=1,6
991 : 867009198 : 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 : 1011510731 : & +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 [ + + ]: 21443 : 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 [ + + ]: 21443 : 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 [ + + - + ]: 21443 : 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 [ + + ]: 21702 : 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 [ + + + + : 887433 : if (ngrad >0) prod (:,:)=prod (:,:)*fact_ucvol
+ + ]
1309 [ + + + + : 26897 : if (ngradp>0) prodp(:,:)=prodp(:,:)*fact_ucvol
+ + ]
1310 [ + + ]: 21443 : 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 [ + + ]: 21443 : 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 [ + + ]: 21443 : 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 : 21443 : jrhoij=1
1365 [ + + ]: 584756 : do irhoij=1,pawrhoij_iatom%nrhoijsel
1366 : 563313 : klmn=pawrhoij_iatom%rhoijselect(irhoij)
1367 : 563313 : klm =pawtab(itypat)%indklmn(1,klmn)
1368 : 563313 : lmin=pawtab(itypat)%indklmn(3,klmn)
1369 : 563313 : lmax=pawtab(itypat)%indklmn(4,klmn)
1370 : 563313 : ro =pawrhoij_iatom%rhoijp(jrhoij,ispden)
1371 : 563313 : ro_d=ro*pawtab(itypat)%dltij(klmn)
1372 [ + - + + ]: 1577260 : do ll=lmin,lmax,2
1373 [ + + ]: 5321765 : do ilm=ll**2+1,(ll+1)**2
1374 : 3744505 : isel=pawang%gntselect(ilm,klm)
1375 [ + + ]: 4758452 : if (isel>0) then
1376 : 990465 : grhat_x=ro_d*pawtab(itypat)%qijl(ilm,klmn)
1377 [ + + ]: 5468478 : do mu=1,ngrad
1378 : 5468478 : 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 : 584756 : 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 [ + + ]: 21443 : 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 [ + + ]: 21443 : 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 [ + + - + ]: 41058 : 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 [ + + ]: 21702 : 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 [ - + ]: 19615 : 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 [ + + ]: 19615 : 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 [ - + ]: 19615 : 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 [ - + ]: 19615 : 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 [ + + ]: 19615 : if (optgr==1) then
1664 : 18091 : mushift=3*(iatm-1)
1665 [ + + ]: 72364 : tmp(1:3)=grhat_tmp(ishift_gr+1:ishift_gr+3,idiag)
1666 [ + + ]: 72364 : do mu=1,3
1667 : 72364 : 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 [ + + ]: 19615 : 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 [ + + ]: 18725 : 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 [ + + ]: 19615 : 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 [ + + ]: 19615 : 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 [ - + ]: 19615 : ABI_FREE(vloc)
1722 [ + + ]: 19615 : if (ngrad>0) then
1723 [ - + ]: 19595 : ABI_FREE(prod)
1724 : : end if
1725 [ + + ]: 19615 : if (ngradp>0) then
1726 [ - + ]: 71 : ABI_FREE(prodp)
1727 : : end if
1728 [ + + - + ]: 31574 : if (optgr2==1.or.optstr2==1) then
1729 [ + + ]: 320 : do jatm=1,natom
1730 [ - + ]: 249 : ABI_FREE(prod_nondiag(jatm)%value)
1731 [ - + ]: 24141 : ABI_FREE(prodp_nondiag(jatm)%value)
1732 : : end do
1733 : : end if
1734 : : end do ! iatm
1735 : 19644 : 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 [ + + ]: 7685 : 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 [ - + ]: 7685 : ABI_FREE(grhat_tmp)
1787 [ + + - + ]: 7685 : 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 [ + + ]: 7685 : if (optgr==1) then
1841 [ + + ]: 72744 : grnl(1:3*natom)=grnl(1:3*natom)+hatgr(1:3*natom)
1842 [ - + ]: 6939 : ABI_FREE(hatgr)
1843 : : end if
1844 : :
1845 : : !===== Convert stresses (add diag and off-diag contributions) =====
1846 [ + + ]: 7685 : if (optstr==1) then
1847 : :
1848 : : ! Has to compute int[nhat*vtrial]. See Eq.(40) in Torrent2008 .
1849 : 1412 : hatstr_diag=zero
1850 [ + + + + ]: 1412 : if (nspden==1.or.dimvtrial==1) then
1851 [ + + ]: 45253780 : do ic=1,nfft
1852 : 45253780 : 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 : 1412 : hatstr_diag=hatstr_diag*fact_ucvol
1866 [ + + ]: 1412 : 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 [ + + ]: 5648 : hatstr(1:3)=(hatstr(1:3)+hatstr_diag)/ucvol
1877 [ + + ]: 5648 : hatstr(4:6)= hatstr(4:6)/ucvol
1878 : :
1879 : : ! Add to already computed NL contrib
1880 [ + + ]: 9884 : nlstr(1:6)=nlstr(1:6)+hatstr(1:6)
1881 : :
1882 : : ! Apply symmetries
1883 : 1412 : call stresssym(gprimd,nsym,nlstr,symrec)
1884 : : end if
1885 : :
1886 : : !===== Convert dynamical matrix (from cartesian to reduced coordinates) =====
1887 [ + + ]: 7685 : 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 [ + + ]: 7685 : 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 [ + + ]: 7685 : if (usexcnhat==0) then
1929 [ - + ]: 5347 : ABI_FREE(vtrial_)
1930 : : end if
1931 : :
1932 : : !Destroy atom tables used for parallelism
1933 : 7685 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1934 [ + + ]: 7685 : if (paral_atom) then
1935 [ - + ]: 3072 : ABI_FREE(atm_indx)
1936 : : end if
1937 : :
1938 : : !Destroy FFT tables used for parallelism
1939 [ + + + + : 15414 : 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 : : !!***
|
