TABLE OF CONTENTS
ABINIT/mklocl_recipspace [ Functions ]
NAME
mklocl_recipspace
FUNCTION
Optionally compute : option=1 : local ionic potential throughout unit cell option=2 : contribution of local ionic potential to E gradient wrt xred option=3 : contribution of local ionic potential to stress tensor option=4 : contribution of local ionic potential to second derivative of E wrt xred
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
if(option==3) eei=local pseudopotential part of total energy (hartree) gmet(3,3)=reciprocal space metric ($\textrm{Bohr}^{-2}$). gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=cutoff on $|G|^2$: see setup1 for definition (doubled sphere). mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization mqgrid=number of grid pts in q array for f(q) spline. natom=number of atoms in unit cell. nattyp(ntypat)=number of 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 ntypat=number of types of atoms. option= (see above) ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information. qgrid(mqgrid)=q grid for spline from 0 to qmax. qprtrb(3)= integer wavevector of possible perturbing potential in basis of reciprocal lattice translations rhog(2,nfft)=electron density rho(G) (electrons/$\textrm{Bohr}^3$) (needed if option==2 or if option==4) ucvol=unit cell volume ($\textrm{Bohr}^3$). vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom. vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero, perturbing potential is added of the form $V(G)=(vprtrb(1)+I*vprtrb(2))/2$ at the values G=qprtrb and $(vprtrb(1)-I*vprtrb(2))/2$ at $G=-qprtrb$ (integers)
OUTPUT
(if option==1) vpsp(nfft)=local crystal pseudopotential in real space. (if option==2) grtn(3,natom)=grads of Etot wrt tn. (if option==3) lpsstr(6)=components of local psp part of stress tensor (Cartesian coordinates, symmetric tensor) in hartree/$\textrm{bohr}^3$ Store 6 unique components in order 11, 22, 33, 32, 31, 21 (if option==4) dyfrlo(3,3,natom)=d2 Eei/d tn(i)/d tn(j). (Hartrees)
SIDE EFFECTS
NOTES
Note that the present routine is tightly connected to the dfpt_vlocal.f routine, that compute the derivative of the local ionic potential with respect to one atomic displacement. The argument list and the internal loops to be considered were sufficiently different as to make the two routine different.
PARENTS
dfpt_dyfro,mklocl,stress
CHILDREN
fourdp,ptabs_fourdp,timab,wrtout,xmpi_sum
SOURCE
72 #if defined HAVE_CONFIG_H 73 #include "config.h" 74 #endif 75 76 #include "abi_common.h" 77 78 subroutine mklocl_recipspace(dyfrlo,eei,gmet,gprimd,grtn,gsqcut,lpsstr,mgfft,& 79 & mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,ntypat,option,paral_kgb,ph1d,qgrid,qprtrb,& 80 & rhog,ucvol,vlspl,vprtrb,vpsp) 81 82 use defs_basis 83 use defs_abitypes 84 use m_profiling_abi 85 use m_errors 86 use m_xmpi 87 88 use m_geometry, only : xred2xcart 89 use m_mpinfo, only : ptabs_fourdp 90 91 !This section has been created automatically by the script Abilint (TD). 92 !Do not modify the following lines by hand. 93 #undef ABI_FUNC 94 #define ABI_FUNC 'mklocl_recipspace' 95 use interfaces_14_hidewrite 96 use interfaces_18_timing 97 use interfaces_53_ffts 98 !End of the abilint section 99 100 implicit none 101 102 !Arguments ------------------------------------ 103 !scalars 104 integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat,option,paral_kgb 105 real(dp),intent(in) :: eei,gsqcut,ucvol 106 type(MPI_type),intent(in) :: mpi_enreg 107 !arrays 108 integer,intent(in) :: nattyp(ntypat),ngfft(18),qprtrb(3) 109 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom) 110 real(dp),intent(in) :: qgrid(mqgrid),rhog(2,nfft),vlspl(mqgrid,2,ntypat) 111 real(dp),intent(in) :: vprtrb(2) 112 real(dp),intent(out) :: dyfrlo(3,3,natom),grtn(3,natom),lpsstr(6) !vz_i 113 real(dp),intent(inout) :: vpsp(nfft) !vz_i 114 115 !Local variables------------------------------- 116 !scalars 117 integer,parameter :: im=2,re=1 118 integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ierr,ig1,ig2,ig3,ii,itypat 119 integer :: jj,me_fft,me_g0,n1,n2,n3,nproc_fft,shift1 120 integer :: shift2,shift3 121 real(dp),parameter :: tolfix=1.0000001_dp 122 real(dp) :: aa,bb,cc,cutoff,dbl_ig1,dbl_ig2,dbl_ig3,dd,diff,dq,dq2div6,dqdiv6 123 real(dp) :: dqm1,ee,ff,gmag,gsquar,ph12i,ph12r,ph1i,ph1r,ph2i,ph2r 124 real(dp) :: ph3i,ph3r,phimag_igia,phre_igia,sfi,sfr 125 real(dp) :: svion,svioni,svionr,term,vion1,vion2,xnorm 126 character(len=500) :: message 127 !arrays 128 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 129 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 130 real(dp) :: gcart(3),tsec(2) 131 real(dp),allocatable :: work1(:,:) 132 133 ! ************************************************************************* 134 135 !Define G^2 based on G space metric gmet. 136 ! gsq(i1,i2,i3)=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 137 !& dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 138 !& dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 139 140 !Real and imaginary parts of phase--statment functions: 141 ! phr(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 142 ! phi(x1,y1,x2,y2,x3,y3)=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 143 ! ph1(nri,i1,ia)=ph1d(nri,i1+1+n1+(ia-1)*(2*n1+1)) 144 ! ph2(nri,i2,ia)=ph1d(nri,i2+1+n2+(ia-1)*(2*n2+1)+& 145 !& natom*(2*n1+1)) 146 ! ph3(nri,i3,ia)=ph1d(nri,i3+1+n3+(ia-1)*(2*n3+1)+& 147 !& natom*(2*n1+1+2*n2+1)) 148 ! phre(i1,i2,i3,ia)=phr(ph1(re,i1,ia),ph1(im,i1,ia),ph2(re,i2,ia),& 149 !& ph2(im,i2,ia),ph3(re,i3,ia),ph3(im,i3,ia)) 150 ! phimag(i1,i2,i3,ia)=phi(ph1(re,i1,ia),ph1(im,i1,ia),ph2(re,i2,ia),& 151 !& ph2(im,i2,ia),ph3(re,i3,ia),ph3(im,i3,ia)) 152 153 !----- 154 155 !Keep track of total time spent in mklocl 156 if(option==2)then 157 call timab(72,1,tsec) 158 end if 159 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 160 me_fft=ngfft(11) 161 nproc_fft=ngfft(10) 162 163 !Get the distrib associated with this fft_grid 164 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 165 166 !Zero out array to permit accumulation over atom types below: 167 if(option==1)then 168 ABI_ALLOCATE(work1,(2,nfft)) 169 work1(:,:)=zero 170 end if 171 ! 172 dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1) 173 dqm1=1.0_dp/dq 174 dqdiv6=dq/6.0_dp 175 dq2div6=dq**2/6.0_dp 176 cutoff=gsqcut*tolfix 177 id1=n1/2+2 178 id2=n2/2+2 179 id3=n3/2+2 180 grtn(:,:)=zero 181 lpsstr(:)=zero 182 dyfrlo(:,:,:)=zero 183 me_g0=0 184 ia1=1 185 186 do itypat=1,ntypat 187 ! ia1,ia2 sets range of loop over atoms: 188 ia2=ia1+nattyp(itypat)-1 189 190 ii=0 191 do i3=1,n3 192 ig3=i3-(i3/id3)*n3-1 193 do i2=1,n2 194 ig2=i2-(i2/id2)*n2-1 195 if(fftn2_distrib(i2) == me_fft ) then 196 do i1=1,n1 197 ig1=i1-(i1/id1)*n1-1 198 199 ii=ii+1 200 ! *** GET RID OF THIS THESE IF STATEMENTS (if they slow code) 201 ! Skip G=0: 202 ! if (ii==1) cycle 203 if (ig1==0 .and. ig2==0 .and. ig3==0) me_g0=1 204 if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 205 206 gsquar=gsq_mk(ig1,ig2,ig3) 207 ! Skip G**2 outside cutoff: 208 if (gsquar<=cutoff) then 209 gmag=sqrt(gsquar) 210 211 ! Compute vion(G) for given type of atom 212 jj=1+int(gmag*dqm1) 213 diff=gmag-qgrid(jj) 214 215 ! Evaluate spline fit from q^2 V(q) to get V(q): 216 ! (p. 86 Numerical Recipes, Press et al; 217 ! NOTE error in book for sign 218 ! of "aa" term in derivative; also see splfit routine). 219 220 bb = diff*dqm1 221 aa = 1.0_dp-bb 222 cc = aa*(aa**2-1.0_dp)*dq2div6 223 dd = bb*(bb**2-1.0_dp)*dq2div6 224 225 vion1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +& 226 & cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat) ) / gsquar 227 228 if(option==1)then 229 230 ! Assemble structure factor over all atoms of given type: 231 sfr=zero 232 sfi=zero 233 do ia=ia1,ia2 234 sfr=sfr+phre_mk(ig1,ig2,ig3,ia) 235 sfi=sfi-phimag_mk(ig1,ig2,ig3,ia) 236 end do 237 ! Multiply structure factor times vion: 238 work1(re,ii)=work1(re,ii)+sfr*vion1 239 work1(im,ii)=work1(im,ii)+sfi*vion1 240 241 else if(option==2 .or. option==4)then 242 243 ! Compute Re and Im part of (2Pi)*Vion(G)*rho(G): 244 svionr=(two_pi*vion1)*rhog(re,ii) 245 svioni=(two_pi*vion1)*rhog(im,ii) 246 247 ! Loop over all atoms of this type: 248 do ia=ia1,ia2 249 shift1=1+n1+(ia-1)*(2*n1+1) 250 shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1) 251 shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1) 252 ph1r=ph1d(1,ig1+shift1) 253 ph1i=ph1d(2,ig1+shift1) 254 ph2r=ph1d(1,ig2+shift2) 255 ph2i=ph1d(2,ig2+shift2) 256 ph3r=ph1d(1,ig3+shift3) 257 ph3i=ph1d(2,ig3+shift3) 258 ph12r=ph1r*ph2r-ph1i*ph2i 259 ph12i=ph1r*ph2i+ph1i*ph2r 260 phre_igia=ph12r*ph3r-ph12i*ph3i 261 phimag_igia=ph12r*ph3i+ph12i*ph3r 262 263 if(option==2)then 264 265 ! Compute "Vion" part of gradient 266 ! svion=svioni*phre(ig1,ig2,ig3,ia)+svionr*phimag(ig1,ig2,ig3,ia) 267 svion=svioni*phre_igia+svionr*phimag_igia 268 269 ! Open loop over 3-index for speed: 270 grtn(1,ia)=grtn(1,ia)-dble(ig1)*svion 271 grtn(2,ia)=grtn(2,ia)-dble(ig2)*svion 272 grtn(3,ia)=grtn(3,ia)-dble(ig3)*svion 273 274 else 275 276 ! Compute "Vion" part of the second derivative 277 ! svion=two_pi* 278 ! (svionr*phre(ig1,ig2,ig3,ia)-svioni*phimag(ig1,ig2,ig3,ia)) 279 svion=two_pi*(svionr*phre_igia-svioni*phimag_igia) 280 281 ! Open loop over 3-index for speed 282 dbl_ig1=dble(ig1) ; dbl_ig2=dble(ig2) ; dbl_ig3=dble(ig3) 283 dyfrlo(1,1,ia)=dyfrlo(1,1,ia)-dbl_ig1*dbl_ig1*svion 284 dyfrlo(1,2,ia)=dyfrlo(1,2,ia)-dbl_ig1*dbl_ig2*svion 285 dyfrlo(1,3,ia)=dyfrlo(1,3,ia)-dbl_ig1*dbl_ig3*svion 286 dyfrlo(2,2,ia)=dyfrlo(2,2,ia)-dbl_ig2*dbl_ig2*svion 287 dyfrlo(2,3,ia)=dyfrlo(2,3,ia)-dbl_ig2*dbl_ig3*svion 288 dyfrlo(3,3,ia)=dyfrlo(3,3,ia)-dbl_ig3*dbl_ig3*svion 289 290 end if 291 292 end do 293 294 else if(option==3)then 295 296 ! Also get (dV(q)/dq)/q: 297 ! (note correction of Numerical Recipes sign error 298 ! before (3._dp*aa**2-1._dp) 299 ! ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines 300 ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat) 301 ff= (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) & 302 & - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat) 303 vion2 = ( ( ee*dqm1 + ff*dqdiv6 )/gmag& 304 & - 2.0_dp*vion1 ) / gsquar 305 306 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+& 307 & gprimd(1,3)*dble(ig3) 308 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+& 309 & gprimd(2,3)*dble(ig3) 310 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+& 311 & gprimd(3,3)*dble(ig3) 312 ! Assemble structure over all atoms of given type 313 sfr=zero 314 sfi=zero 315 do ia=ia1,ia2 316 sfr=sfr+phre_mk(ig1,ig2,ig3,ia) 317 sfi=sfi-phimag_mk(ig1,ig2,ig3,ia) 318 end do 319 320 ! Compute Re( rho^*(G)* sf ) * [(dV(G)/dG)/|G|] 321 term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)*vion2 322 323 ! Compute contribution to stress tensor 324 lpsstr(1)=lpsstr(1)-term*gcart(1)*gcart(1) 325 lpsstr(2)=lpsstr(2)-term*gcart(2)*gcart(2) 326 lpsstr(3)=lpsstr(3)-term*gcart(3)*gcart(3) 327 lpsstr(4)=lpsstr(4)-term*gcart(3)*gcart(2) 328 lpsstr(5)=lpsstr(5)-term*gcart(3)*gcart(1) 329 lpsstr(6)=lpsstr(6)-term*gcart(2)*gcart(1) 330 331 else 332 write(message, '(a,i0,a)' )' mklocl: Option=',option,' not allowed.' 333 MSG_BUG(message) 334 end if ! End option choice 335 336 ! End skip G**2 outside cutoff: 337 end if 338 339 ! End loop on n1, n2, n3. There is a "cycle" inside the loop 340 end do 341 end if ! this plane is for me_fft 342 end do 343 end do 344 345 ! Symmetrize the dynamical matrix with respect to indices 346 do ia=ia1,ia2 347 dyfrlo(2,1,ia)=dyfrlo(1,2,ia) 348 dyfrlo(3,1,ia)=dyfrlo(1,3,ia) 349 dyfrlo(3,2,ia)=dyfrlo(2,3,ia) 350 end do 351 352 ia1=ia2+1 353 354 ! End loop on type of atoms 355 end do 356 357 if(option==1)then 358 ! Dont't change work1 on g=0 if Poisson solver is used since work1 359 ! hold not the potential but the density generated by the pseudo. 360 if(me_g0 == 1) then 361 ! Set Vloc(G=0)=0: 362 work1(re,1)=zero 363 work1(im,1)=zero 364 end if 365 366 ! DEBUG 367 ! write(std_out,*) ' mklocl_recipspace : will add potential with strength vprtrb(:)=',vprtrb(:) 368 ! ENDDEBUG 369 370 ! Allow for the addition of a perturbing potential 371 if ((vprtrb(1)**2+vprtrb(2)**2) > 1.d-30) then 372 ! Find the linear indices which correspond with the input 373 ! wavevector qprtrb 374 ! The double modulus handles both i>=n and i<0, mapping into [0,n-1]; 375 ! then add 1 to get range [1,n] for each 376 i3=1+mod(n3+mod(qprtrb(3),n3),n3) 377 i2=1+mod(n2+mod(qprtrb(2),n2),n2) 378 i1=1+mod(n1+mod(qprtrb(1),n1),n1) 379 ! Compute the linear index in the 3 dimensional array 380 ii=i1+n1*((ffti2_local(i2)-1)+(n2/nproc_fft)*(i3-1)) 381 ! Add in the perturbation at G=qprtrb 382 work1(re,ii)=work1(re,ii)+0.5_dp*vprtrb(1) 383 work1(im,ii)=work1(im,ii)+0.5_dp*vprtrb(2) 384 ! Same thing for G=-qprtrb 385 i3=1+mod(n3+mod(-qprtrb(3),n3),n3) 386 i2=1+mod(n2+mod(-qprtrb(2),n2),n2) 387 i1=1+mod(n1+mod(-qprtrb(1),n1),n1) 388 ! ii=i1+n1*((i2-1)+n2*(i3-1)) 389 work1(re,ii)=work1(re,ii)+0.5_dp*vprtrb(1) 390 work1(im,ii)=work1(im,ii)-0.5_dp*vprtrb(2) 391 write(message, '(a,1p,2e12.4,a,0p,3i4,a)' )& 392 & ' mklocl: perturbation of vprtrb=', vprtrb,& 393 & ' and q=',qprtrb,' has been added' 394 call wrtout(std_out,message,'COLL') 395 end if 396 397 ! Transform back to real space 398 call fourdp(1,work1,vpsp,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 399 400 ! Divide by unit cell volume 401 xnorm=1.0_dp/ucvol 402 vpsp(:)=vpsp(:)*xnorm 403 404 ABI_DEALLOCATE(work1) 405 406 end if 407 408 if(option==2)then 409 ! Init mpi_comm 410 if(mpi_enreg%nproc_fft>1)then 411 call timab(48,1,tsec) 412 call xmpi_sum(grtn,mpi_enreg%comm_fft ,ierr) 413 call timab(48,2,tsec) 414 end if 415 call timab(72,2,tsec) 416 end if 417 418 if(option==3)then 419 ! Init mpi_comm 420 if(mpi_enreg%nproc_fft>1)then 421 call timab(48,1,tsec) 422 call xmpi_sum(lpsstr,mpi_enreg%comm_fft ,ierr) 423 call timab(48,2,tsec) 424 end if 425 426 ! Normalize and add term -eei/ucvol on diagonal 427 ! (see page 802 of notes) 428 lpsstr(1)=(lpsstr(1)-eei)/ucvol 429 lpsstr(2)=(lpsstr(2)-eei)/ucvol 430 lpsstr(3)=(lpsstr(3)-eei)/ucvol 431 lpsstr(4)=lpsstr(4)/ucvol 432 lpsstr(5)=lpsstr(5)/ucvol 433 lpsstr(6)=lpsstr(6)/ucvol 434 435 end if 436 437 if(option==4)then 438 ! Init mpi_comm 439 if(mpi_enreg%nproc_fft>1)then 440 call timab(48,1,tsec) 441 call xmpi_sum(dyfrlo,mpi_enreg%comm_fft ,ierr) 442 call timab(48,2,tsec) 443 end if 444 end if 445 446 contains 447 448 !Real and imaginary parts of phase--statment functions: 449 function phr_mk(x1,y1,x2,y2,x3,y3) 450 451 452 !This section has been created automatically by the script Abilint (TD). 453 !Do not modify the following lines by hand. 454 #undef ABI_FUNC 455 #define ABI_FUNC 'phr_mk' 456 !End of the abilint section 457 458 real(dp) :: phr_mk,x1,x2,x3,y1,y2,y3 459 phr_mk=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3 460 end function phr_mk 461 462 function phi_mk(x1,y1,x2,y2,x3,y3) 463 464 465 !This section has been created automatically by the script Abilint (TD). 466 !Do not modify the following lines by hand. 467 #undef ABI_FUNC 468 #define ABI_FUNC 'phi_mk' 469 !End of the abilint section 470 471 real(dp):: phi_mk,x1,x2,x3,y1,y2,y3 472 phi_mk=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3 473 end function phi_mk 474 475 function ph1_mk(nri,ig1,ia) 476 477 478 !This section has been created automatically by the script Abilint (TD). 479 !Do not modify the following lines by hand. 480 #undef ABI_FUNC 481 #define ABI_FUNC 'ph1_mk' 482 !End of the abilint section 483 484 real(dp):: ph1_mk 485 integer :: nri,ig1,ia 486 ph1_mk=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1)) 487 end function ph1_mk 488 489 function ph2_mk(nri,ig2,ia) 490 491 492 !This section has been created automatically by the script Abilint (TD). 493 !Do not modify the following lines by hand. 494 #undef ABI_FUNC 495 #define ABI_FUNC 'ph2_mk' 496 !End of the abilint section 497 498 real(dp):: ph2_mk 499 integer :: nri,ig2,ia 500 ph2_mk=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)) 501 end function ph2_mk 502 503 function ph3_mk(nri,ig3,ia) 504 505 506 !This section has been created automatically by the script Abilint (TD). 507 !Do not modify the following lines by hand. 508 #undef ABI_FUNC 509 #define ABI_FUNC 'ph3_mk' 510 !End of the abilint section 511 512 real(dp):: ph3_mk 513 integer :: nri,ig3,ia 514 ph3_mk=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)) 515 end function ph3_mk 516 517 function phre_mk(ig1,ig2,ig3,ia) 518 519 520 !This section has been created automatically by the script Abilint (TD). 521 !Do not modify the following lines by hand. 522 #undef ABI_FUNC 523 #define ABI_FUNC 'phre_mk' 524 !End of the abilint section 525 526 real(dp):: phre_mk 527 integer :: ig1,ig2,ig3,ia 528 phre_mk=phr_mk(ph1_mk(re,ig1,ia),ph1_mk(im,ig1,ia),& 529 & ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia)) 530 end function phre_mk 531 532 function phimag_mk(ig1,ig2,ig3,ia) 533 534 535 !This section has been created automatically by the script Abilint (TD). 536 !Do not modify the following lines by hand. 537 #undef ABI_FUNC 538 #define ABI_FUNC 'phimag_mk' 539 !End of the abilint section 540 541 real(dp) :: phimag_mk 542 integer :: ig1,ig2,ig3,ia 543 phimag_mk=phi_mk(ph1_mk(re,ig1,ia),ph1_mk(im,ig1,ia),& 544 & ph2_mk(re,ig2,ia),ph2_mk(im,ig2,ia),ph3_mk(re,ig3,ia),ph3_mk(im,ig3,ia)) 545 end function phimag_mk 546 547 function gsq_mk(i1,i2,i3) 548 549 550 !This section has been created automatically by the script Abilint (TD). 551 !Do not modify the following lines by hand. 552 #undef ABI_FUNC 553 #define ABI_FUNC 'gsq_mk' 554 !End of the abilint section 555 556 real(dp) :: gsq_mk 557 integer :: i1,i2,i3 558 gsq_mk=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 559 & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 560 & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 561 end function gsq_mk 562 563 end subroutine mklocl_recipspace