TABLE OF CONTENTS
ABINIT/nonlop_ylm [ Functions ]
NAME
nonlop_ylm
FUNCTION
* Compute application of a nonlocal operator Vnl in order to get:
- contracted elements (energy, forces, stresses, ...), if signs=1
- a function in reciprocal space (|out> = Vnl|in>), if signs=2
Operator Vnl, as the following general form:
$Vnl=sum_{R,lmn,l''m''n''} {|P_{Rlmn}> Enl^{R}_{lmn,l''m''n''} <P_{Rl''m''n''}|}$
Operator Vnl is -- in the typical case -- the nonlocal potential.
- With norm-conserving pseudopots, $Enl^{R}_{lmn,l''m''n''}$ is the
Kleinmann-Bylander energy $Ekb^{R}_{ln}$.
- In a PAW calculation, $Enl^{R}_{lmn,l''m''n''}$ are the nonlocal
coefficients to connect projectors $D_{ij}$.
- The |P_{Rlmn}> are the projector functions.
* Optionnaly, in case of PAW calculation, compute:
- Application of the overlap matrix in reciprocal space
(<in|S|in> or (I+S)|in>).
- Application of (Vnl-lambda.S) in reciprocal space
(<in|Vnl-lambda.S|in> and derivatives or (Vnl-lambda.S)|in>).
* This routine uses spherical harmonics Ylm to express Vnl.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx
choice: chooses possible output:
choice=0 => do nothing (only compute WF projected with NL projectors)
=1 => non-local energy contribution
=2 => 1st derivative(s) with respect to atomic position(s)
=3 => 1st derivative(s) with respect to strain(s)
=23=> 1st derivative(s) with respect to atomic pos. and
1st derivative(s) with respect to atomic pos. and strains
=4 => 2nd derivative(s) with respect to 2 atomic pos.
=24=> 1st derivative(s) with respect to atm. pos. and
2nd derivative(s) with respect to 2 atomic pos.
=5 => 1st derivative(s) with respect to k wavevector, typically
sum_ij [ |p_i> D_ij <dp_j/dk| + |dp_i/dk> D_ij < p_j| ]
=6 => 2nd derivative(s) with respect to 2 strains and
mixed 2nd derivative(s) with respect to strains & atomic pos.
=51 =>right 1st derivative(s) with respect to k wavevector, typically
sum_ij [ |p_i> D_ij <dp_j/dk| ]
=52 =>left 1st derivative(s) with respect to k wavevector, typically
sum_ij [ |dp_i/dk> D_ij < p_j| ]
=53 =>twist 1st derivative(s) with respect to k, typically
sum_ij [ |dp_i/dk_(idir+1)> D_ij <dp_j//dk_(idir-1)|
-|dp_i/dk_(idir-1)> D_ij <dp_j//dk_(idir+1)|]
=54=> mixed 2nd derivative(s) with respect to atomic pos. and left k wavevector
=55=> mixed 2nd derivative(s) with respect to strain and right k wavevector
=7 => apply operator $\sum_i [ |p_i> <p_i| ],
same as overlap operator with s_ij=identity (paw_opt==3 only)
=8 => 2nd derivatives with respect to 2 k wavevectors
=81=> partial 2nd derivatives with respect to 2 k wavevectors,
full derivative with respect to k1, right derivative with respect to k2,
(derivative with respect to k of choice 51), typically
sum_ij [ |dp_i/dk1> D_ij <dp_j/dk2| + |p_i> D_ij < d2p_j/dk1dk2| ]
Only choices 1,2,3,23,4,5,6 are compatible with useylm=0.
Only choices 1,2,3,5,51,52,53,7,8,81 are compatible with signs=2
cpopt=flag defining the status of cprjin%cp(:)=<Proj_i|Cnk> scalars (see below, side effects)
dimenl1,dimenl2=dimensions of enl (see enl)
dimffnlin=second dimension of ffnlin (1+number of derivatives)
dimffnlout=second dimension of ffnlout (1+number of derivatives)
enl(dimenl1,dimenl2,nspinortot**2)=
->Norm conserving : ==== when paw_opt=0 ====
(Real) Kleinman-Bylander energies (hartree)
dimenl1=lmnmax - dimenl2=ntypat
->PAW : ==== when paw_opt=1, 2 or 4 ====
(Real or complex, hermitian) Dij coefs to connect projectors
dimenl1=cplex_enl*lmnmax*(lmnmax+1)/2 - dimenl2=natom
These are complex numbers if cplex_enl=2
enl(:,:,1) contains Dij^up-up
enl(:,:,2) contains Dij^dn-dn
enl(:,:,3) contains Dij^up-dn (only if nspinor=2)
enl(:,:,4) contains Dij^dn-up (only if nspinor=2)
ffnlin(npwin,dimffnlin,lmnmax,ntypat)=nonlocal form factors to be used
for the application of the nonlocal operator to the |in> vector
ffnlout(npwout,dimffnlout,lmnmax,ntypat)=nonlocal form factors to be used
for the application of the nonlocal operator to the |out> vector
-----------------------------------------------------------
gprimd(3,3)=dimensional reciprocal space primitive translations
hermdij,optional,logical :: if true forces Hermitial Dij behavior in opernlc_ylm
idir=direction of the - atom to be moved in the case (choice=2,signs=2),
- k point direction in the case (choice=5,signs=2S)
for choice 53, twisted derivative involves idir+1 and idir-1
- strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=lmn
istwf_k=option parameter that describes the storage of wfs
kgin(3,npwin)=integer coords of planewaves in basis sphere, for the |in> vector
kgout(3,npwout)=integer coords of planewaves in basis sphere, for the |out> vector
kpgin(npw,npkgin)= (k+G) components and related data, for the |in> vector
kpgout(npw,nkpgout)=(k+G) components and related data, for the |out> vector
kptin(3)=k point in terms of recip. translations, for the |in> vector
kptout(3)=k point in terms of recip. translations, for the |out> vector
lambda=factor to be used when computing (Vln-lambda.S) - only for paw_opt=2
Typically lambda is the eigenvalue (or its guess)
lmnmax=max. number of (l,m,n) components over all types of atoms
matblk=dimension of the arrays ph3din and ph3dout
mgfft=maximum size of 1D FFTs
mpi_enreg=information about MPI parallelization
natom=number of atoms in cell
nattyp(ntypat)=number of atoms of each type
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
nkpgin,nkpgout=second sizes of arrays kpgin/kpgout
nloalg(3)=governs the choice of the algorithm for nonlocal operator
nnlout=dimension of enlout (when signs=1 and choice>0):
==== if paw_opt=0, 1 or 2 ====
choice nnlout | choice nnlout
1 1 | 51 6 (complex)
2 3*natom | 52 6 (complex)
3 6 | 53 6
4 6*natom | 54 9*natom
23 6+3*natom | 55 36 (complex)
24 9*natom | 6 36+18*natom
5 3 | 8 6
| 81 18 (complex)
==== if paw_opt=3 ====
choice nnlout
1 1
2 3*natom
5 3
51 3
52 3
54 9*natom
55 36
7 1
8 6
81 9
==== if paw_opt=4 ====
not available
npwin=number of planewaves for given k point, for the |in> vector
npwout=number of planewaves for given k point, for the |out> vector
nspinor=number of spinorial components of the wavefunctions (on current proc)
nspinortot=total number of spinorial components of the wavefunctions
ntypat=number of types of atoms in cell
paw_opt= define the nonlocal operator concerned with:
paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
paw_opt=3 : PAW overlap matrix (Sij)
paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
phkxredin(2,natom)=phase factors exp(2 pi kptin.xred)
phkxredout(2,natom)=phase factors exp(2 pi kptout.xred)
ph1d(2,3*(2*mgfft+1)*natom)=1D structure factors phase information
ph3din(2,npwin,matblk)=3D structure factors, for each atom and plane wave (in)
ph3dout(2,npwout,matblk)=3-dim structure factors, for each atom and plane wave (out)
signs= if 1, get contracted elements (energy, forces, stress, ...)
if 2, applies the non-local operator to a function in reciprocal space
sij(dimenl1,ntypat*(paw_opt/3))=overlap matrix components (only if paw_opt=2, 3 or 4)
ucvol=unit cell volume (bohr^3)
vectin(2,npwin*nspinor)=input cmplx wavefunction coefficients <G|in>
[cprjin_left(natom,nspinor)]=The projected input wave function <p_nlm|in_left>
for the left wavefunction. Data are assumed to be in memory, they are NOT recalculated here.
OUTPUT
==== if (signs==1) ====
--If (paw_opt==0, 1 or 2)
enlout(nnlout)= contribution to the non-local part of the following properties:
if choice=1 : enlout(1) -> the energy
if choice=2 : enlout(3*natom) -> 1st deriv. of energy wrt atm. pos (forces)
if choice=3 : enlout(6) -> 1st deriv. of energy wrt strain (stresses)
if choice=4 : enlout(6*natom) -> 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
if choice=23: enlout(6+3*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
1st deriv. of energy wrt strain (stresses)
if choice=24: enlout(9*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
if choice=5 : enlout(3) -> 1st deriv. of energy wrt k
if choice=51: enlout(3) -> 1st deriv. (right) of energy wrt k
if choice=52: enlout(3) -> 1st deriv. (left) of energy wrt k
if choice=53: enlout(3) -> 1st deriv. (twist) of energy wrt k
if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
if choice=6 : enlout(36+18*natom) -> 2nd deriv. of energy wrt 2 strains (elast. tensor) and
2nd deriv. of energy wrt to atm. pos and strain (internal strain)
if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
if choice=81: enlout(9) -> 2nd deriv.of E: full derivative w.r.t. k1, right derivative w.r.t k2
--If (paw_opt==3)
if choice=1 : enlout(1) -> contribution to <c|S|c> (note: not including <c|c>)
if choice=2 : enlout(3*natom) -> contribution to <c|dS/d_atm.pos|c>
if choice=51: enlout(3) -> contribution to <c|d(right)S/d_k|c>
if choice=52: enlout(3) -> contribution to <c|d(left)S/d_k|c>
if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
if choice=7 : enlout(1) -> contribution to <c|sum_i[p_i><p_i]|c>
if choice=8 : enlout(6) -> contribution to <c|d2S/d_k1d_k2|c>
if choice=81: enlout(9) -> contribution to <c|dS/d_k1[d(right)d_k2]|c>
--If (paw_opt==4)
not available
==== if (signs==2) ====
--if (paw_opt=0, 1 or 4)
vectout(2,npwout*my_nspinor*ndat)=result of the aplication of the concerned operator
or one of its derivatives to the input vect.:
if (choice=1) <G|V_nonlocal|vect_in>
if (choice=2) <G|dV_nonlocal/d(atm. pos)|vect_in>
if (choice=3) <G|dV_nonlocal/d(strain)|vect_in>
if (choice=5) <G|dV_nonlocal/d(k)|vect_in>
if (choice=51) <G|d(right)V_nonlocal/d(k)|vect_in>
if (choice=52) <G|d(left)V_nonlocal/d(k)|vect_in>
if (choice=53) <G|d(twist)V_nonlocal/d(k)|vect_in>
if (choice=8) <G|d2V_nonlocal/d(k)d(k)|vect_in>
if (choice=81) <G|d[d(right)V_nonlocal/d(k)]/d(k)|vect_in>
--if (paw_opt=2)
vectout(2,npwout*my_nspinor*ndat)=final vector in reciprocal space:
if (choice=1) <G|V_nonlocal-lamdba.(I+S)|vect_in>
if (choice=2) <G|d[V_nonlocal-lamdba.(I+S)]/d(atm. pos)|vect_in>
if (choice=3) <G|d[V_nonlocal-lamdba.(I+S)]/d(strain)|vect_in>
if (choice=5) <G|d[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
if (choice=51) <G|d(right)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
if (choice=52) <G|d(left)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
if (choice=53) <G|d(twist)[V_nonlocal-lamdba.(I+S)]/d(k)|vect_in>
if (choice=8) <G|d2[V_nonlocal-lamdba.(I+S)]/d(k)d(k)|vect_in>
if (choice=81) <G|d[d(right[V_nonlocal-lamdba.(I+S)]/d(k)]/d(k)|vect_in>
--if (paw_opt=3 or 4)
svectout(2,npwout*my_nspinor*ndat)=result of the aplication of Sij (overlap matrix)
or one of its derivatives to the input vect.:
if (choice=1) <G|I+S|vect_in>
if (choice=2) <G|dS/d(atm. pos)|vect_in>
if (choice=3) <G|dS/d(strain)|vect_in>
if (choice=5) <G|dS/d(k)|vect_in>
if (choice=51) <G|d(right)S/d(k)|vect_in>
if (choice=52) <G|d(left)S/d(k)|vect_in>
if (choice=53) <G|d(twist)S/d(k)|vect_in>
if (choice=3) <G|d[V_nonlocal-lamdba.(I+S)]/d(strain)|vect_in>
if (choice=7) <G|sum_i[p_i><p_i]|vect_in>
if (choice=8) <G|d2S/d(k)d(k)|vect_in>
if (choice=81) <G|d[d(right)S/d(k)]/d(k)|vect_in>
SIDE EFFECTS
cprjin(natom,nspinor) <type(pawcprj_type)>=projected input wave function |in> on non-local projectors =<p_lmn|in> and derivatives Treatment depends on cpopt parameter: if cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved) if cpopt= 0, <p_lmn|in> are computed here and saved derivatives are eventually computed but not saved if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved other derivatives are eventually computed but not saved if cpopt= 2 <p_lmn|in> are already in memory; first (and 2nd) derivatives are computed here and not saved if cpopt= 3 <p_lmn|in> are already in memory; first derivatives are computed here and saved other derivatives are eventually computed but not saved if cpopt= 4 <p_lmn|in> and first derivatives are already in memory; other derivatives are not computed (except when choice=8 or 81) This option is not compatible with choice=4,24 or 6 Warning: for cpopt= 1 or 3, derivatives wrt strains do not contain the contribution due to the volume change; i.e. <dp_lmn/dEps|in> are incomplete.
NOTES
This application of the nonlocal operator is programmed using a direct implementation of spherical harmonics (Ylm). Abinit used historically Legendre polynomials for the application of nonlocal operator; but the implementation of PAW algorithm enforced the use of Ylm. In the case signs=1, the array vectout is not used, nor modified so that the same array as vectin can be used as a dummy argument; the same is true for the pairs npwin-npwout, ffnlin-ffnlout, kgin-kgout, ph3din-ph3dout, phkredin-phkxredout).
TODO
* Complete implementation of spin-orbit
PARENTS
nonlop
CHILDREN
mkkpg,opernla_ylm,opernlb_ylm,opernlc_ylm,opernld_ylm,ph1d3d,strconv xmpi_sum
SOURCE
279 #if defined HAVE_CONFIG_H 280 #include "config.h" 281 #endif 282 283 #include "abi_common.h" 284 285 subroutine nonlop_ylm(atindx1,choice,cpopt,cprjin,dimenl1,dimenl2,dimffnlin,dimffnlout,& 286 & enl,enlout,ffnlin,ffnlout,gprimd,idir,indlmn,istwf_k,& 287 & kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,lmnmax,matblk,mgfft,& 288 & mpi_enreg,natom,nattyp,ngfft,nkpgin,nkpgout,nloalg,nnlout,& 289 & npwin,npwout,nspinor,nspinortot,ntypat,paw_opt,phkxredin,phkxredout,ph1d,& 290 & ph3din,ph3dout,signs,sij,svectout,ucvol,vectin,vectout,cprjin_left,hermdij) 291 292 use defs_basis 293 use defs_abitypes 294 use m_xmpi 295 use m_profiling_abi 296 use m_errors 297 298 use m_kg, only : ph1d3d 299 use m_pawcprj, only : pawcprj_type 300 301 !This section has been created automatically by the script Abilint (TD). 302 !Do not modify the following lines by hand. 303 #undef ABI_FUNC 304 #define ABI_FUNC 'nonlop_ylm' 305 use interfaces_41_geometry 306 use interfaces_66_nonlocal, except_this_one => nonlop_ylm 307 !End of the abilint section 308 309 implicit none 310 311 !Arguments ------------------------------------ 312 !scalars 313 integer,intent(in) :: choice,cpopt,dimenl1,dimenl2,dimffnlin,dimffnlout,idir 314 integer,intent(in) :: istwf_k,lmnmax,matblk,mgfft,natom,nkpgin,nkpgout,nnlout 315 integer,intent(in) :: npwin,npwout,nspinor,nspinortot,ntypat,paw_opt,signs 316 real(dp),intent(in) :: lambda,ucvol 317 logical,optional,intent(in) :: hermdij 318 type(MPI_type),intent(in) :: mpi_enreg 319 !arrays 320 integer,intent(in) :: atindx1(natom),kgin(3,npwin) 321 integer,intent(in),target :: indlmn(6,lmnmax,ntypat) 322 integer,intent(in) :: kgout(3,npwout),nattyp(ntypat),ngfft(18),nloalg(3) 323 real(dp),intent(in) :: enl(dimenl1,dimenl2,nspinortot**2) 324 real(dp),intent(in),target :: ffnlin(npwin,dimffnlin,lmnmax,ntypat) 325 real(dp),intent(in),target :: ffnlout(npwout,dimffnlout,lmnmax,ntypat) 326 real(dp),intent(in) :: gprimd(3,3) 327 real(dp),intent(in),target :: kpgin(npwin,nkpgin),kpgout(npwout,nkpgout) 328 real(dp),intent(in) :: kptin(3),kptout(3),ph1d(2,3*(2*mgfft+1)*natom) 329 real(dp),intent(in) :: phkxredin(2,natom),phkxredout(2,natom) 330 real(dp),intent(in) :: sij(dimenl1,ntypat*((paw_opt+1)/3)) 331 real(dp),intent(inout) :: ph3din(2,npwin,matblk),ph3dout(2,npwout,matblk) 332 real(dp),intent(inout) :: vectin(:,:) 333 real(dp),intent(out) :: enlout(:) 334 real(dp),intent(out) :: svectout(:,:) 335 real(dp),intent(inout) :: vectout (:,:) 336 type(pawcprj_type),intent(inout) :: cprjin(:,:) 337 type(pawcprj_type),optional,intent(in) :: cprjin_left(:,:) 338 339 !Local variables------------------------------- 340 !scalars 341 integer :: choice_a,choice_b,cplex,cplex_enl,cplex_fac,ia,ia1,ia2,ia3,ia4,ia5 342 integer :: iatm,ic,idir1,idir2,ii,ierr,ilmn,iln,ishift,ispinor,itypat,jc,mincat,mu,mua,mub,mu0 343 integer :: n1,n2,n3,nd2gxdt,ndgxdt,ndgxdt_stored,nd2gxdtfac,ndgxdtfac 344 integer :: nincat,nkpgin_,nkpgout_,nlmn,nu,nua1,nua2,nub1,nub2,optder 345 real(dp) :: enlk 346 logical :: check,hermdij_,testnl 347 character(len=500) :: message 348 !arrays 349 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 350 integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 351 integer,allocatable :: cplex_dgxdt(:),cplex_d2gxdt(:) 352 integer,ABI_CONTIGUOUS pointer :: indlmn_typ(:,:) 353 real(dp),allocatable :: d2gxdt(:,:,:,:,:),d2gxdtfac(:,:,:,:,:),d2gxdtfac_sij(:,:,:,:,:) 354 real(dp),allocatable :: dgxdt(:,:,:,:,:),dgxdtfac(:,:,:,:,:),dgxdtfac_sij(:,:,:,:,:) 355 real(dp),allocatable :: ddkk(:),fnlk(:),gmet(:,:) 356 real(dp),allocatable :: gx(:,:,:,:),gxfac(:,:,:,:),gxfac_sij(:,:,:,:),gx_left(:,:,:,:) 357 real(dp),allocatable :: sij_typ(:),strnlk(:) 358 real(dp),allocatable :: work1(:),work2(:),work3(:,:),work4(:,:),work5(:,:,:),work6(:,:,:),work7(:,:,:) 359 real(dp),ABI_CONTIGUOUS pointer :: ffnlin_typ(:,:,:),ffnlout_typ(:,:,:),kpgin_(:,:),kpgout_(:,:) 360 361 ! ********************************************************************** 362 363 DBG_ENTER("COLL") 364 365 !Check consistency of arguments 366 !============================================================== 367 368 hermdij_=.FALSE.; if(present(hermdij)) hermdij_=hermdij 369 370 !signs=1, almost all choices 371 if (signs==1) then 372 if(paw_opt<3) then 373 check=(choice==0 .or.choice==1 .or.choice==2 .or.choice==3 .or.choice==4 .or.& 374 & choice==23.or.choice==24.or.choice==5 .or.choice==51.or.choice==52.or.& 375 & choice==53.or.choice==54.or.choice==55.or.& 376 & choice==6 .or.choice==8 .or.choice==81) 377 else if (paw_opt==3) then 378 check=(choice== 0.or.choice== 1.or.choice== 2.or.choice==3.or.choice==5.or.& 379 & choice==23.or.choice==51.or.choice==52.or.choice==54.or.choice==55.or.& 380 & choice== 8.or.choice==81) 381 else 382 check = .false. 383 end if 384 ABI_CHECK(check,'BUG: choice not compatible (for signs=1)') 385 end if 386 387 !signs=2, less choices 388 if (signs==2) then 389 check=(choice==0.or.choice==1.or.choice==2.or.choice==3 .or.& 390 & choice==5.or.choice==51.or.choice==52.or.choice==53.or.choice==54.or.& 391 & choice==7.or.choice==8.or.choice==81) 392 ABI_CHECK(check,'BUG: choice not compatible (for signs=2)') 393 end if 394 !1<=idir<=6 is required when choice=3 and signs=2 395 if (choice==3.and.signs==2) then 396 check=(idir>=1.and.idir<=6) 397 ABI_CHECK(check,'BUG: choice=3 and signs=2 requires 1<=idir<=6') 398 !1<=idir<=9 is required when choice==8/81 and signs=2 399 else if ((choice==8.or.choice==81.or.choice==54).and.signs==2) then 400 check=(idir>=1.and.idir<=9) 401 ABI_CHECK(check,'BUG: choice=8/81 and signs=2 requires 1<=idir<=9') 402 else 403 ! signs=2 requires 1<=idir<=3 when choice>1 404 check=(signs/=2.or.choice<=1.or.choice==7.or.(idir>=1.and.idir<=3)) 405 ABI_CHECK(check,'BUG: signs=2 requires 1<=idir<=3') 406 end if 407 !check allowed values for cpopt 408 check=(cpopt>=-1.and.cpopt<=4) 409 ABI_CHECK(check,'bad value for cpopt') 410 check=(cpopt/=4.or.(choice/=4.and.choice/=24.and.choice/=6)) 411 ABI_CHECK(check,'BUG: cpopt=4 not allowed for 2nd derivatives') 412 check=(cpopt/=2.or.(choice/=8.and.choice/=81)) 413 ABI_CHECK(check,'BUG: cpopt=2 not allowed for choice=8,81, use cpopt=4 instead') 414 !check conditions for optional arguments 415 check=((.not.present(cprjin_left)).or.(signs==1.and.choice==1)) 416 ABI_CHECK(check,'BUG: when cprjin_left is present, must have choice=1,signs=1') 417 !protect special case choice==7 418 check=(choice/=7.or.paw_opt==3) 419 ABI_CHECK(check,'BUG: when choice=7, paw_opt must be 3') 420 !spin-orbit not yet allowed 421 check=(maxval(indlmn(6,:,:))<=1) 422 ABI_CHECK(check,'BUG: spin-orbit not yet allowed') 423 424 !Test: size of blocks of atoms 425 mincat=min(NLO_MINCAT,maxval(nattyp)) 426 if (nloalg(2)<=0.and.mincat>matblk) then 427 write(message, '(a,a,a,i4,a,i4,a)' ) & 428 & 'With nloc_mem<=0, mincat must be less than matblk.',ch10,& 429 & 'Their value is ',mincat,' and ',matblk,'.' 430 MSG_BUG(message) 431 end if 432 433 !Define dimensions of projected scalars 434 !============================================================== 435 436 !Define some useful variables 437 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3) 438 choice_a=merge(choice,1,choice/=7);choice_b=choice_a 439 if (cpopt>=2) choice_a=-choice_a 440 cplex=2;if (istwf_k>1) cplex=1 !Take into account TR-symmetry 441 cplex_enl=1;if (paw_opt>0) cplex_enl=2*dimenl1/(lmnmax*(lmnmax+1)) 442 cplex_fac=cplex;if ((nspinortot==2.or.cplex_enl==2).and.paw_opt>0.and.choice/=7) cplex_fac=2 443 444 !Define dimensions of projected scalars 445 ndgxdt=0;ndgxdtfac=0;nd2gxdt=0;nd2gxdtfac=0 446 if (choice==2) then 447 if (signs==1) ndgxdt=3 448 if (signs==2) ndgxdt=1 449 if (signs==2) ndgxdtfac=1 450 end if 451 if (choice==3) then 452 if (signs==1) ndgxdt=6 453 if (signs==2) ndgxdt=1 454 if (signs==2) ndgxdtfac=1 455 end if 456 if (choice==23) then 457 if (signs==1) ndgxdt=9 458 end if 459 if (choice==4) then 460 if(signs==1) ndgxdt=3 461 if(signs==1) ndgxdtfac=3 462 if(signs==1) nd2gxdt=6 463 end if 464 if (choice==24) then 465 if(signs==1) ndgxdt=3 466 if(signs==1) ndgxdtfac=3 467 if(signs==1) nd2gxdt=6 468 end if 469 if (choice==5) then 470 if(signs==1) ndgxdt=3 471 if(signs==2) ndgxdt=1 472 if(signs==2) ndgxdtfac=1 473 end if 474 if (choice==51) then 475 if(signs==1) ndgxdt=3 476 if(signs==2) ndgxdt=1 477 if(signs==2) ndgxdtfac=1 478 end if 479 if (choice==52) then 480 if(signs==1) ndgxdt=3 481 if(signs==2) ndgxdt=1 482 if(signs==2) ndgxdtfac=1 483 end if 484 if (choice==53) then 485 if(signs==1) ndgxdt=3 486 if(signs==1) ndgxdtfac=3 487 if(signs==2) ndgxdt=2 488 if(signs==2) ndgxdtfac=2 489 end if 490 if (choice==54) then 491 if(signs==1) ndgxdt=6 492 if(signs==1) ndgxdtfac=6 493 if(signs==1) nd2gxdt=9 494 if(signs==2) ndgxdt=1 495 if(signs==2) nd2gxdt=1 496 if(signs==2) ndgxdtfac=1 497 if(signs==2) nd2gxdtfac=1 498 end if 499 if (choice==55) then 500 if(signs==1) ndgxdt=9 501 if(signs==1) ndgxdtfac=9 502 if(signs==1) nd2gxdt=18 503 end if 504 if (choice==6) then 505 if(signs==1) ndgxdt=9 506 if(signs==1) ndgxdtfac=9 507 if(signs==1) nd2gxdt=54 508 end if 509 if (choice==8) then 510 if(signs==1) ndgxdt=3 511 if(signs==1) ndgxdtfac=3 512 if(signs==1) nd2gxdt=6 513 if(signs==2) ndgxdt=2 514 if(signs==2) ndgxdtfac=2 515 if(signs==2) nd2gxdt=1 516 if(signs==2) nd2gxdtfac=1 517 end if 518 if (choice==81) then 519 if(signs==1) ndgxdt=3 520 if(signs==1) ndgxdtfac=3 521 if(signs==1) nd2gxdt=6 522 if(signs==2) ndgxdt=1 523 if(signs==2) ndgxdtfac=1 524 if(signs==2) nd2gxdt=1 525 if(signs==2) nd2gxdtfac=1 526 end if 527 ABI_CHECK(ndgxdtfac<=ndgxdt,"BUG: ndgxdtfac>ndgxdt!") 528 optder=0;if (ndgxdtfac>0) optder=1 529 if (nd2gxdtfac>0) optder=2 530 531 !Consistency tests 532 if (cpopt==4) then 533 if (ndgxdt>0.and.cprjin(1,1)%ncpgr<=0) then 534 message='cprjin%ncpgr=0 not allowed with cpopt=4 and these (choice,signs) !' 535 MSG_BUG(message) 536 end if 537 end if 538 if (cpopt==1.or.cpopt==3) then 539 if (cprjin(1,1)%ncpgr<ndgxdt) then 540 message='should have cprjin%ncpgr>=ndgxdt with cpopt=1 or 3 !' 541 MSG_BUG(message) 542 end if 543 end if 544 545 !Additional steps before calculation 546 !============================================================== 547 548 !Initialize output arrays 549 if (signs==1) then 550 ABI_ALLOCATE(fnlk,(3*natom)) 551 ABI_ALLOCATE(ddkk,(6)) 552 ABI_ALLOCATE(strnlk,(6)) 553 enlk=zero;fnlk=zero;ddkk=zero;strnlk=zero 554 enlout(:)=zero 555 end if 556 if (signs==2) then 557 if (paw_opt==0.or.paw_opt==1.or.paw_opt==4) vectout(:,:)=zero 558 if (paw_opt==2.and.choice==1) vectout(:,:)=-lambda*vectin(:,:) 559 if (paw_opt==2.and.choice> 1) vectout(:,:)=zero 560 if (paw_opt==3.or.paw_opt==4) then 561 if (choice==1) svectout(:,:)=vectin(:,:) 562 if (choice> 1) svectout(:,:)=zero 563 end if 564 end if 565 566 !Eventually re-compute (k+G) vectors (and related data) 567 nkpgin_=0 568 if (choice==2.or.choice==54) nkpgin_=3 569 if (signs==1) then 570 if (choice==4.or.choice==24) nkpgin_=9 571 if (choice==3.or.choice==23.or.choice==6) nkpgin_=3 572 if (choice==55) nkpgin_=3 573 end if 574 if (nkpgin<nkpgin_) then 575 ABI_ALLOCATE(kpgin_,(npwin,nkpgin_)) 576 call mkkpg(kgin,kpgin_,kptin,nkpgin_,npwin) 577 else 578 nkpgin_ = nkpgin 579 kpgin_ => kpgin 580 end if 581 nkpgout_=0 582 if ((choice==2.or.choice==54).and.signs==2) nkpgout_=3 583 if (nkpgout<nkpgout_) then 584 ABI_ALLOCATE(kpgout_,(npwout,nkpgout_)) 585 call mkkpg(kgout,kpgout_,kptout,nkpgout_,npwout) 586 else 587 nkpgout_ = nkpgout 588 kpgout_ => kpgout 589 end if 590 591 !Big loop on atom types 592 !============================================================== 593 594 ia1=1;iatm=0 595 do itypat=1,ntypat 596 597 ! Get atom loop indices for different types: 598 ia2=ia1+nattyp(itypat)-1;ia5=1 599 600 ! Select quantities specific to the current type of atom 601 nlmn=count(indlmn(3,:,itypat)>0) 602 603 ! Temporary test on local part 604 testnl=(paw_opt/=0) 605 if (paw_opt==0) then 606 do ispinor=1,nspinortot 607 do ilmn=1,nlmn 608 iln=indlmn(5,ilmn,itypat) 609 if (abs(enl(iln,itypat,ispinor))>tol10) testnl=.true. 610 end do 611 end do 612 end if 613 614 ! Some non-local part is to be applied for that type of atom 615 if (testnl) then 616 617 ! Store some quantities depending only of the atom type 618 ffnlin_typ => ffnlin(:,:,:,itypat) 619 indlmn_typ => indlmn(:,:,itypat) 620 if (signs==2) then 621 ffnlout_typ => ffnlout(:,:,:,itypat) 622 end if 623 if (paw_opt>=2) then 624 ABI_ALLOCATE(sij_typ,(nlmn*(nlmn+1)/2)) 625 if (cplex_enl==1) then 626 do ilmn=1,nlmn*(nlmn+1)/2 627 sij_typ(ilmn)=sij(ilmn,itypat) 628 end do 629 else 630 do ilmn=1,nlmn*(nlmn+1)/2 631 sij_typ(ilmn)=sij(2*ilmn-1,itypat) 632 end do 633 end if 634 else 635 ABI_ALLOCATE(sij_typ,(0)) 636 end if 637 638 ! Loop over atoms of the same type 639 ! ============================================================== 640 641 ! Cut the sum on different atoms in blocks, to allow memory saving. 642 ! Inner summations on atoms will be done from ia3 to ia4. 643 ! Note: the maximum range from ia3 to ia4 is mincat (max. increment of atoms). 644 645 do ia3=ia1,ia2,mincat 646 ia4=min(ia2,ia3+mincat-1) 647 ! Give the increment of number of atoms in this subset. 648 nincat=ia4-ia3+1 649 650 ! Prepare the phase factors if they were not already computed 651 if (nloalg(2)<=0) then 652 call ph1d3d(ia3,ia4,kgin,matblk,natom,npwin,n1,n2,n3,phkxredin,ph1d,ph3din) 653 end if 654 655 ! Allocate memory for projected scalars 656 ABI_ALLOCATE(gx,(cplex,nlmn,nincat,nspinor)) 657 ABI_ALLOCATE(dgxdt,(cplex,ndgxdt,nlmn,nincat,nspinor)) 658 ABI_ALLOCATE(d2gxdt,(cplex,nd2gxdt,nlmn,nincat,nspinor)) 659 ABI_ALLOCATE(d2gxdtfac,(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor)) 660 ABI_ALLOCATE(dgxdtfac,(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)) 661 ABI_ALLOCATE(gxfac,(cplex_fac,nlmn,nincat,nspinor)) 662 gx(:,:,:,:)=zero;gxfac(:,:,:,:)=zero 663 if (ndgxdt>0) dgxdt(:,:,:,:,:)=zero 664 if (ndgxdtfac>0) dgxdtfac(:,:,:,:,:)=zero 665 if (nd2gxdt>0) d2gxdt(:,:,:,:,:)=zero 666 if (nd2gxdtfac>0) d2gxdtfac(:,:,:,:,:)=zero 667 if (paw_opt>=3) then 668 ABI_ALLOCATE(gxfac_sij,(cplex,nlmn,nincat,nspinor)) 669 ABI_ALLOCATE(dgxdtfac_sij,(cplex,ndgxdtfac,nlmn,nincat,nspinor)) 670 ABI_ALLOCATE(d2gxdtfac_sij,(cplex,nd2gxdtfac,nlmn,nincat,nspinor)) 671 gxfac_sij(:,:,:,:)=zero 672 if (ndgxdtfac>0) dgxdtfac_sij(:,:,:,:,:)=zero 673 if (nd2gxdtfac>0) d2gxdtfac_sij(:,:,:,:,:) = zero 674 else 675 ABI_ALLOCATE(gxfac_sij,(0,0,0,0)) 676 ABI_ALLOCATE(dgxdtfac_sij,(0,0,0,0,0)) 677 ABI_ALLOCATE(d2gxdtfac_sij,(0,0,0,0,0)) 678 end if 679 680 ! When istwf_k > 1, gx derivatives can be real or pure imaginary 681 ! cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary 682 ! cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary 683 ABI_ALLOCATE(cplex_dgxdt,(ndgxdt)) 684 ABI_ALLOCATE(cplex_d2gxdt,(nd2gxdt)) 685 cplex_dgxdt(:) = 1 ; cplex_d2gxdt(:) = 1 686 if(ndgxdt > 0) then 687 if (choice==5.or.choice==51.or.choice==52.or.choice==53.or. & 688 & choice==8.or.choice==81) cplex_dgxdt(:) = 2 689 if (choice==54.and.signs==1) cplex_dgxdt(4:6) = 2 690 if (choice==54.and.signs==2) cplex_dgxdt(:) = 2 691 if (choice==55.and.signs==1) cplex_dgxdt(7:9) = 2 692 end if 693 if(nd2gxdt > 0) then 694 if (choice==54) cplex_d2gxdt(:) = 2 695 if (choice==55.and.signs==1) cplex_d2gxdt(1:18)= 2 696 end if 697 698 ! Compute projection of current wave function |c> on each 699 ! non-local projector: <p_lmn|c> 700 ! ============================================================== 701 702 ! Retrieve eventually <p_lmn|c> coeffs (and derivatives) 703 if (cpopt>=2) then 704 do ispinor=1,nspinor 705 do ia=1,nincat 706 gx(1:cplex,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%cp(1:cplex,1:nlmn) 707 end do 708 end do 709 end if 710 if (cpopt==4.and.ndgxdt>0) then 711 ndgxdt_stored = cprjin(1,1)%ncpgr 712 ishift=0 713 if (((choice==2).or.(choice==3)).and.(ndgxdt_stored>ndgxdt).and.(signs==2)) ishift=idir-ndgxdt 714 if ((choice==2).and.(ndgxdt_stored==9).and.(signs==2)) ishift=ishift+6 715 if (choice==2.and.(ndgxdt_stored>ndgxdt).and.(signs==1)) ishift=ndgxdt_stored-ndgxdt 716 if(cplex == 2) then 717 do ispinor=1,nspinor 718 do ia=1,nincat 719 if (ndgxdt_stored==ndgxdt.or.(ndgxdt_stored>ndgxdt.and.((choice==2).or.(choice==3)))) then 720 dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn) 721 else if (signs==2.and.ndgxdt_stored==3) then 722 if (choice==5.or.choice==51.or.choice==52) then ! ndgxdt=1 723 dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir,1:nlmn) 724 else if (choice==8) then ! ndgxdt=2 725 idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1 726 dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir1,1:nlmn) 727 dgxdt(1:2,2,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir2,1:nlmn) 728 else if (choice==81) then ! ndgxdt=1 729 idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1 730 dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir2,1:nlmn) 731 end if 732 end if 733 end do 734 end do 735 else 736 do ispinor=1,nspinor 737 do ia=1,nincat 738 do ilmn=1,nlmn 739 if (ndgxdt_stored==ndgxdt.or.(ndgxdt_stored>ndgxdt.and.((choice==2).or.(choice==3)))) then 740 do ii=1,ndgxdt 741 ic = cplex_dgxdt(ii) 742 dgxdt(1,ii,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn) 743 end do 744 else if (signs==2.and.ndgxdt_stored==3) then 745 if (choice==5.or.choice==51.or.choice==52) then ! ndgxdt=1 746 dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir,ilmn) 747 else if (choice==8) then ! ndgxdt=2 748 idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1 749 dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir1,ilmn) 750 dgxdt(1,2,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(2),idir2,ilmn) 751 else if (choice==81) then ! ndgxdt=1 752 idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1 753 dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir2,ilmn) 754 end if 755 end if 756 end do 757 end do 758 end do 759 end if 760 end if 761 762 ! Computation or <p_lmn|c> (and derivatives) for this block of atoms 763 if ((cpopt<4.and.choice_a/=-1).or.choice==8.or.choice==81) then 764 call opernla_ylm(choice_a,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnlin,d2gxdt,dgxdt,ffnlin_typ,gx,& 765 & ia3,idir,indlmn_typ,istwf_k,kpgin_,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpgin_,nlmn,& 766 & nloalg,npwin,nspinor,ph3din,signs,ucvol,vectin) 767 end if 768 769 ! Transfer result to output variable cprj (if requested) 770 ! cprj(:)%cp receive the <p_i|Psi> factors (p_i: non-local projector) 771 ! Be careful: cprj(:)%dcp does not exactly contain the derivative of cprj(:)%cp. 772 ! - Volume contributions (in case of strain derivative) are not included, 773 ! - Global coordinate transformation are not applied, 774 ! cprj(:)%dcp is meant to be a restart argument of the present nonlop routine. 775 if (cpopt==0.or.cpopt==1) then 776 do ispinor=1,nspinor 777 do ia=1,nincat 778 cprjin(iatm+ia,ispinor)%nlmn=nlmn 779 cprjin(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)=gx(1:cplex,1:nlmn,ia,ispinor) 780 if (cplex==1) cprjin(iatm+ia,ispinor)%cp(2,1:nlmn)=zero 781 end do 782 end do 783 end if 784 if ((cpopt==1.or.cpopt==3).and.ndgxdt>0) then 785 ishift=0 786 if ((choice==2).and.(cprjin(1,1)%ncpgr>ndgxdt)) ishift=cprjin(1,1)%ncpgr-ndgxdt 787 if(cplex==2)then 788 do ispinor=1,nspinor 789 do ia=1,nincat 790 cprjin(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn)=dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor) 791 ! cprjin(iatm+ia,ispinor)%dcp(1:2,1:ndgxdt,1:nlmn)=dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor) 792 end do 793 end do 794 else 795 do ispinor=1,nspinor 796 do ia=1,nincat 797 do ilmn=1,nlmn 798 do ii=1,ndgxdt 799 ic = cplex_dgxdt(ii) ; jc = 3 - ic 800 cprjin(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor) 801 cprjin(iatm+ia,ispinor)%dcp(jc,ii+ishift,ilmn)=zero 802 ! cprjin(iatm+ia,ispinor)%dcp(ic,ii,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor) 803 ! cprjin(iatm+ia,ispinor)%dcp(jc,ii,ilmn)=zero 804 end do 805 end do 806 end do 807 end do 808 end if 809 end if 810 811 ! If choice==0, that's all for these atoms ! 812 if (choice>0) then 813 if(choice/=7) then 814 ! Contraction from <p_i|c> to Sum_j[Dij.<p_j|c>] (and derivatives) 815 call opernlc_ylm(atindx1,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_enl,cplex_fac,dgxdt,dgxdtfac,dgxdtfac_sij,& 816 & d2gxdt,d2gxdtfac,d2gxdtfac_sij,dimenl1,dimenl2,enl,gx,gxfac,gxfac_sij,& 817 & iatm,indlmn_typ,itypat,lambda,mpi_enreg,natom,ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,& 818 & nincat,nlmn,nspinor,nspinortot,optder,paw_opt,sij_typ,hermdij=hermdij_) 819 else 820 gxfac_sij=gx 821 end if 822 823 ! Operate with the non-local potential on the projected scalars, 824 ! in order to get contributions to energy/forces/stress/dyn.mat 825 ! ============================================================== 826 if (signs==1) then 827 if (.not.present(cprjin_left)) then 828 call opernld_ylm(choice_b,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,& 829 & enlk,enlout,fnlk,gx,gxfac,gxfac_sij,ia3,natom,nd2gxdt,ndgxdt,ndgxdtfac,& 830 & nincat,nlmn,nnlout,nspinor,paw_opt,strnlk) 831 else 832 ABI_ALLOCATE(gx_left,(cplex,nlmn,nincat,nspinor)) 833 ! Retrieve <p_lmn|c> coeffs 834 do ispinor=1,nspinor 835 do ia=1,nincat 836 gx_left(1:cplex,1:nlmn,ia,ispinor)=cprjin_left(iatm+ia,ispinor)%cp(1:cplex,1:nlmn) 837 end do 838 end do 839 ! TODO 840 ! if (cpopt==4.and.ndgxdt>0) then 841 ! do ispinor=1,nspinor 842 ! do ia=1,nincat 843 ! dgxdt_left(1:cplex,1:ndgxdt,1:nlmn,ia,ispinor)=cprjin_left(iatm+ia,ispinor)%dcp(1:cplex,1:ndgxdt,1:nlmn) 844 ! end do 845 ! end do 846 ! end if 847 call opernld_ylm(choice_b,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,& 848 & enlk,enlout,fnlk,gx_left,gxfac,gxfac_sij,ia3,natom,nd2gxdt,ndgxdt,ndgxdtfac,& 849 & nincat,nlmn,nnlout,nspinor,paw_opt,strnlk) 850 ABI_DEALLOCATE(gx_left) 851 end if 852 end if 853 854 ! Operate with the non-local potential on the projected scalars, 855 ! in order to get matrix element 856 ! ============================================================== 857 if (signs==2) then 858 ! Prepare the phase factors if they were not already computed 859 if(nloalg(2)<=0) then 860 call ph1d3d(ia3,ia4,kgout,matblk,natom,npwout,n1,n2,n3,phkxredout,ph1d,ph3dout) 861 end if 862 call opernlb_ylm(choice_b,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,& 863 & d2gxdtfac,d2gxdtfac_sij,dgxdtfac,dgxdtfac_sij,dimffnlout,ffnlout_typ,gxfac,gxfac_sij,ia3,& 864 & idir,indlmn_typ,kpgout_,matblk,ndgxdtfac,nd2gxdtfac,nincat,nkpgout_,nlmn,& 865 & nloalg,npwout,nspinor,paw_opt,ph3dout,svectout,ucvol,vectout) 866 end if 867 868 end if ! choice==0 869 870 ! Deallocate temporary projected scalars 871 ABI_DEALLOCATE(gx) 872 ABI_DEALLOCATE(gxfac) 873 ABI_DEALLOCATE(dgxdt) 874 ABI_DEALLOCATE(dgxdtfac) 875 ABI_DEALLOCATE(d2gxdt) 876 ABI_DEALLOCATE(d2gxdtfac) 877 ABI_DEALLOCATE(dgxdtfac_sij) 878 ABI_DEALLOCATE(d2gxdtfac_sij) 879 ABI_DEALLOCATE(gxfac_sij) 880 ABI_DEALLOCATE(cplex_dgxdt) 881 ABI_DEALLOCATE(cplex_d2gxdt) 882 883 ! End sum on atom subset loop 884 iatm=iatm+nincat;ia5=ia5+nincat 885 end do 886 !if (paw_opt>=2) then 887 ABI_DEALLOCATE(sij_typ) 888 !end if 889 890 ! End condition of existence of a non-local part 891 else 892 if (cpopt==0.or.cpopt==1) then 893 do ispinor=1,nspinor 894 do ia=1,nattyp(itypat) 895 cprjin(iatm+ia,ispinor)%cp(:,1:nlmn)=zero 896 end do 897 end do 898 end if 899 if ((cpopt==1.or.cpopt==3).and.ndgxdt>0) then 900 ishift=0 901 if ((choice==2).and.(cprjin(1,1)%ncpgr>ndgxdt)) ishift=cprjin(1,1)%ncpgr-ndgxdt 902 do ispinor=1,nspinor 903 do ia=1,nattyp(itypat) 904 cprjin(iatm+ia,ispinor)%dcp(:,1+ishift:ndgxdt+ishift,1:nlmn)=zero 905 ! cprjin(iatm+ia,ispinor)%dcp(:,1:ndgxdt,1:nlmn)=zero 906 end do 907 end do 908 end if 909 iatm=iatm+nattyp(itypat) 910 end if 911 912 ! End atom type loop 913 ia1=ia2+1 914 end do 915 916 917 !Reduction in case of parallelism 918 !============================================================== 919 920 if (signs==1.and.mpi_enreg%paral_spinor==1) then 921 if (nnlout/=0) then 922 call xmpi_sum(enlout,mpi_enreg%comm_spinor,ierr) 923 end if 924 if (choice==3.or.choice==6.or.choice==23) then 925 call xmpi_sum(enlk,mpi_enreg%comm_spinor,ierr) 926 end if 927 if (choice==6) then 928 call xmpi_sum(fnlk,mpi_enreg%comm_spinor,ierr) 929 call xmpi_sum(strnlk,mpi_enreg%comm_spinor,ierr) 930 end if 931 if (choice==55) then 932 call xmpi_sum(ddkk,mpi_enreg%comm_spinor,ierr) 933 end if 934 end if 935 936 !Coordinate transformations 937 !============================================================== 938 939 !Need sometimes gmet 940 if ((signs==1.and.paw_opt<=3).and. & 941 & (choice==5 .or.choice==51.or.choice==52.or.choice==53.or.& 942 & choice==54.or.choice==55)) then 943 ABI_ALLOCATE(gmet,(3,3)) 944 gmet = MATMUL(TRANSPOSE(gprimd),gprimd) 945 end if 946 947 !1st derivative wrt to strain (stress tensor): 948 ! - convert from reduced to cartesian coordinates 949 ! - substract volume contribution 950 if ((choice==3.or.choice==23).and.signs==1.and.paw_opt<=3) then 951 mu0=0 ! Shift to be applied in enlout array 952 ABI_ALLOCATE(work1,(6)) 953 work1(1:6)=enlout(mu0+1:mu0+6) 954 call strconv(work1,gprimd,work1) 955 enlout(mu0+1:mu0+3)=(work1(1:3)-enlk) 956 enlout(mu0+4:mu0+6)= work1(4:6) 957 ABI_DEALLOCATE(work1) 958 end if 959 960 !1st derivative wrt to k wave vector (ddk): 961 ! - convert from cartesian to reduced coordinates 962 if ((choice==5.or.choice==53).and.signs==1.and.paw_opt<=3) then 963 mu0=0 ! Shift to be applied in enlout array 964 ABI_ALLOCATE(work1,(3)) 965 work1(:)=enlout(mu0+1:mu0+3) 966 enlout(mu0+1:mu0+3)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3) 967 ABI_DEALLOCATE(work1) 968 end if 969 if ((choice==51.or.choice==52).and.signs==1.and.paw_opt<=3) then 970 mu0=0 ! Shift to be applied in enlout array 971 ABI_ALLOCATE(work1,(3)) 972 do mu=1,2 ! Loop for Re,Im 973 work1(1:3)=(/enlout(mu0+1),enlout(mu0+3),enlout(mu0+5)/) 974 enlout(mu0+1)=gmet(1,1)*work1(1)+gmet(1,2)*work1(2)+gmet(1,3)*work1(3) 975 enlout(mu0+3)=gmet(2,1)*work1(1)+gmet(2,2)*work1(2)+gmet(2,3)*work1(3) 976 enlout(mu0+5)=gmet(3,1)*work1(1)+gmet(3,2)*work1(2)+gmet(3,3)*work1(3) 977 mu0=mu0+1 978 end do 979 ABI_DEALLOCATE(work1) 980 end if 981 982 !2nd derivative wrt to k wave vector and atomic position (effective charges): 983 ! - convert from cartesian to reduced coordinates 984 if (choice==54.and.signs==1.and.paw_opt<=3) then 985 mu0=0 ! Shift to be applied in enlout array 986 ABI_ALLOCATE(work1,(3)) 987 ABI_ALLOCATE(work2,(3)) 988 do mu=1,3*natom 989 ! First, real part 990 work1(1)=enlout(mu0+1);work1(2)=enlout(mu0+3);work1(3)=enlout(mu0+5) 991 work2(:)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3) 992 enlout(mu0+1)=work2(1);enlout(mu0+3)=work2(2);enlout(mu0+5)=work2(3) 993 ! Then imaginary part 994 work1(1)=enlout(mu0+2);work1(2)=enlout(mu0+4);work1(3)=enlout(mu0+6) 995 work2(:)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3) 996 enlout(mu0+2)=work2(1);enlout(mu0+4)=work2(2);enlout(mu0+6)=work2(3) 997 mu0=mu0+6 998 end do 999 ABI_DEALLOCATE(work1) 1000 ABI_DEALLOCATE(work2) 1001 end if 1002 1003 !2nd derivative wrt to k wave vector and strain (piezoelectric tensor): 1004 ! - convert from cartesian to reduced coordinates (k point) 1005 ! - convert from reduced to cartesian coordinates (strain) 1006 ! - substract volume contribution 1007 ! - symetrize strain components 1008 if (choice==55.and.signs==1.and.paw_opt<=3) then 1009 ABI_ALLOCATE(work3,(2,3)) 1010 ABI_ALLOCATE(work4,(2,3)) 1011 ABI_ALLOCATE(work5,(2,3,6)) 1012 ABI_ALLOCATE(work7,(2,3,6)) 1013 ABI_ALLOCATE(work6,(2,3,3)) 1014 do ic=1,3 ! gamma 1015 work5=zero 1016 do jc=1,3 ! nu 1017 do ii=1,3 ! lambda 1018 mu=(gamma(jc,ii)-1)*3+1 1019 work5(1,jc,ii)=gmet(ic,1)*enlout(2*mu-1)+gmet(ic,2)*enlout(2*mu+1) & 1020 & +gmet(ic,3)*enlout(2*mu+3) 1021 work5(2,jc,ii)=gmet(ic,1)*enlout(2*mu )+gmet(ic,2)*enlout(2*mu+2) & 1022 & +gmet(ic,3)*enlout(2*mu+4) 1023 end do 1024 end do 1025 work6=zero 1026 do jc=1,3 ! nu 1027 do ii=1,3 ! beta 1028 work6(1:cplex,ii,jc)=gprimd(ii,1)*work5(1:cplex,jc,1)+gprimd(ii,2)*work5(1:cplex,jc,2) & 1029 & +gprimd(ii,3)*work5(1:cplex,jc,3) 1030 end do 1031 end do 1032 do jc=1,3 ! alpha 1033 do ii=1,3 ! beta 1034 mu=gamma(jc,ii) 1035 work7(1:cplex,ic,mu)=gprimd(jc,1)*work6(1:cplex,ii,1)+gprimd(jc,2)*work6(1:cplex,ii,2) & 1036 & +gprimd(jc,3)*work6(1:cplex,ii,3) 1037 end do 1038 end do 1039 end do ! gamma 1040 1041 do ii=1,3 ! alpha 1042 work3(1,ii)=gprimd(ii,1)*ddkk(2*1-1)+gprimd(ii,2)*ddkk(2*2-1) & 1043 & +gprimd(ii,3)*ddkk(2*3-1) 1044 work3(2,ii)=gprimd(ii,1)*ddkk(2*1 )+gprimd(ii,2)*ddkk(2*2 ) & 1045 & +gprimd(ii,3)*ddkk(2*3 ) 1046 end do 1047 do ii=1,3 ! gamma 1048 work4(1,ii)=gmet(ii,1)*ddkk(2*1-1)+gmet(ii,2)*ddkk(2*2-1) & 1049 & +gmet(ii,3)*ddkk(2*3-1) 1050 work4(2,ii)=gmet(ii,1)*ddkk(2*1 )+gmet(ii,2)*ddkk(2*2 ) & 1051 & +gmet(ii,3)*ddkk(2*3 ) 1052 end do 1053 1054 do mu=1,6 1055 ii=alpha(mu) ! alpha 1056 ic=beta(mu) ! beta 1057 do jc=1,3 ! gamma 1058 work7(1:cplex,jc,mu)=work7(1:cplex,jc,mu)-half & 1059 & *(gprimd(ic,jc)*work3(1:cplex,ii)+gprimd(ii,jc)*work3(1:cplex,ic)) 1060 if (ii==ic) work7(1:cplex,jc,mu)=work7(1:cplex,jc,mu)-work4(1:cplex,jc) 1061 end do 1062 end do 1063 do mu=1,6 ! alpha,beta 1064 do nu=1,3 ! gamma 1065 mu0=3*(mu-1)+nu 1066 enlout(2*mu0-1)=work7(1,nu,mu) 1067 enlout(2*mu0 )=work7(2,nu,mu) 1068 end do 1069 end do 1070 ABI_DEALLOCATE(gmet) 1071 ABI_DEALLOCATE(work3) 1072 ABI_DEALLOCATE(work4) 1073 ABI_DEALLOCATE(work5) 1074 ABI_DEALLOCATE(work6) 1075 ABI_DEALLOCATE(work7) 1076 end if 1077 1078 !2nd derivative wrt to 2 k wave vectors (effective mass): 1079 ! - convert from cartesian to reduced coordinates 1080 if ((choice==8.or.choice==81).and.signs==1.and.paw_opt<=3) then 1081 mu0=0 ! Shift to be applied in enlout array 1082 ABI_ALLOCATE(work3,(3,3)) 1083 ABI_ALLOCATE(work4,(3,3)) 1084 mua=1;if (choice==81) mua=2 1085 do ii=1,mua ! Loop Re,Im 1086 if (choice==8) then ! enlout is real in Voigt notation 1087 work3(1,1)=enlout(mu0+1) ; work3(1,2)=enlout(mu0+6) ; work3(1,3)=enlout(mu0+5) 1088 work3(2,1)=enlout(mu0+6) ; work3(2,2)=enlout(mu0+2) ; work3(2,3)=enlout(mu0+4) 1089 work3(3,1)=enlout(mu0+5) ; work3(3,2)=enlout(mu0+4) ; work3(3,3)=enlout(mu0+3) 1090 else ! enlout is complex in matrix notation 1091 work3(1,1)=enlout(mu0+1 ) ; work3(1,2)=enlout(mu0+3 ) ; work3(1,3)=enlout(mu0+5 ) 1092 work3(2,1)=enlout(mu0+7 ) ; work3(2,2)=enlout(mu0+9 ) ; work3(2,3)=enlout(mu0+11) 1093 work3(3,1)=enlout(mu0+13) ; work3(3,2)=enlout(mu0+15) ; work3(3,3)=enlout(mu0+17) 1094 end if 1095 do mu=1,3 1096 work4(:,mu)=gprimd(:,1)*work3(mu,1)+gprimd(:,2)*work3(mu,2)+gprimd(:,3)*work3(mu,3) 1097 end do 1098 do mu=1,3 1099 work3(:,mu)=gprimd(:,1)*work4(mu,1)+gprimd(:,2)*work4(mu,2)+gprimd(:,3)*work4(mu,3) 1100 end do 1101 do mu=1,3 1102 work4(:,mu)=gprimd(1,:)*work3(mu,1)+gprimd(2,:)*work3(mu,2)+gprimd(3,:)*work3(mu,3) 1103 end do 1104 do mu=1,3 1105 work3(:,mu)=gprimd(1,:)*work4(mu,1)+gprimd(2,:)*work4(mu,2)+gprimd(3,:)*work4(mu,3) 1106 end do 1107 if (choice==8) then ! enlout is real in Voigt notation 1108 enlout(mu0+1) = work3(1,1) ; enlout(mu0+2) = work3(2,2) ; enlout(mu0+3) = work3(3,3) 1109 enlout(mu0+4) = work3(3,2) ; enlout(mu0+5) = work3(1,3) ; enlout(mu0+6) = work3(2,1) 1110 else ! enlout is complex in matrix notation 1111 enlout(mu0+1 )=work3(1,1) ; enlout(mu0+3 )=work3(1,2) ; enlout(mu0+5 )=work3(1,3) 1112 enlout(mu0+7 )=work3(2,1) ; enlout(mu0+9 )=work3(2,2) ; enlout(mu0+11)=work3(2,3) 1113 enlout(mu0+13)=work3(3,1) ; enlout(mu0+15)=work3(3,2) ; enlout(mu0+17)=work3(3,3) 1114 end if 1115 mu0=mu0+1 1116 end do 1117 ABI_DEALLOCATE(work3) 1118 ABI_DEALLOCATE(work4) 1119 end if 1120 1121 !2nd derivative wrt to 2 strains (elastic tensor): 1122 ! - convert from reduced to cartesian coordinates 1123 ! - substract volume contribution 1124 if (choice==6.and.signs==1.and.paw_opt<=3) then 1125 mu0=0 ! Shift to be applied in enlout array 1126 ABI_ALLOCATE(work1,(6)) 1127 ABI_ALLOCATE(work2,(6)) 1128 ABI_ALLOCATE(work3,(6+3*natom,6)) 1129 work3(:,:)=reshape(enlout(mu0+1:mu0+6*(6+3*natom)),(/6+3*natom,6/)) 1130 do mu=1,6 1131 call strconv(work3(1:6,mu),gprimd,work3(1:6,mu)) 1132 end do 1133 do mu=1,6+3*natom 1134 work1(1:6)=work3(mu,1:6) 1135 call strconv(work1,gprimd,work2) 1136 work3(mu,1:6)=work2(1:6) 1137 end do 1138 enlout(mu0+1:mu0+6*(6+3*natom))=reshape(work3(:,:),(/6*(6+3*natom)/)) 1139 ABI_DEALLOCATE(work1) 1140 ABI_DEALLOCATE(work2) 1141 ABI_DEALLOCATE(work3) 1142 call strconv(strnlk,gprimd,strnlk) 1143 do mub=1,6 1144 nub1=alpha(mub);nub2=beta(mub) 1145 do mua=1,6 1146 mu=mu0+mua+(3*natom+6)*(mub-1) 1147 nua1=alpha(mua);nua2=beta(mua) 1148 if (mua<=3.and.mub<=3) enlout(mu)=enlout(mu)+enlk 1149 if (mua<=3) enlout(mu)=enlout(mu)-strnlk(mub) 1150 if (mub<=3) enlout(mu)=enlout(mu)-strnlk(mua) 1151 if (nub1==nua2) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua1,nub2)) 1152 if (nub2==nua2) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua1,nub1)) 1153 if (nub1==nua1) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua2,nub2)) 1154 if (nub2==nua1) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua2,nub1)) 1155 end do 1156 if (mub<=3) then 1157 do nua1=1,natom 1158 nua2=3*(nua1-1);mu=mu0+nua2+6+(3*natom+6)*(mub-1) 1159 enlout(mu+1:mu+3)=enlout(mu+1:mu+3)-fnlk(nua2+1:nua2+3) 1160 end do 1161 end if 1162 end do 1163 end if 1164 1165 if (allocated(gmet)) then 1166 ABI_DEALLOCATE(gmet) 1167 end if 1168 1169 !Final deallocations 1170 !============================================================== 1171 1172 if (signs==1) then 1173 ABI_DEALLOCATE(fnlk) 1174 ABI_DEALLOCATE(ddkk) 1175 ABI_DEALLOCATE(strnlk) 1176 end if 1177 if (nkpgin<nkpgin_) then 1178 ABI_DEALLOCATE(kpgin_) 1179 end if 1180 if (nkpgout<nkpgout_) then 1181 ABI_DEALLOCATE(kpgout_) 1182 end if 1183 1184 DBG_EXIT("COLL") 1185 1186 end subroutine nonlop_ylm