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