TABLE OF CONTENTS
- ABINIT/hartredq.F90
- ABINIT/m_spacepar
- m_spacepar/hartre
- m_spacepar/hartrestr
- m_spacepar/irrzg
- m_spacepar/laplacian
- m_spacepar/make_vectornd
- m_spacepar/meanvalue_g
- m_spacepar/mkunitpawspherepot
- m_spacepar/redgr
- m_spacepar/rotate_rho
- m_spacepar/setsym
- m_spacepar/symrhg
ABINIT/hartredq.F90 [ Functions ]
NAME
hartredq.F90
FUNCTION
Given rho(G), compute the q-gradient of the Hartree potential at q=0 (=FFT of -rho(G)*G_qdir/pi**2/|G|**4 ) -> Cartesian coordinates The calculation is performed in reduced reciprocal space coordinates.
COPYRIGHT
Copyright (C) 2021-2022 ABINIT group (FIXME: add author) 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
cplex= if 1, vqgradhartr is REAL, if 2, vqgradhartr is COMPLEX gmet(3,3)=metrix tensor in G space in Bohr**-2. gprimd(3,3)=reciprocal space dimensional primitive translations gsqcut=cutoff value on G**2 for sphere inside fft box. (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)) mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/input_variables/vargs.htm#ngfft qdir= indicates the direction of the q-gradient (1,2 or 3) rhog(2,nfft)=electron density in G space
OUTPUT
vqgradhart(cplex*nfft)=q-gradient of the Hartree potential at q=0in real space, either REAL or COMPLEX
SIDE EFFECTS
NOTES
SOURCE
2647 subroutine hartredq(cplex,gmet,gsqcut,mpi_enreg,nfft,ngfft,qdir,rhog,vqgradhart) 2648 2649 implicit none 2650 2651 !Arguments ------------------------------------ 2652 !scalars 2653 integer,intent(in) :: cplex,nfft,qdir 2654 real(dp),intent(in) :: gsqcut 2655 type(MPI_type),intent(in) :: mpi_enreg 2656 !arrays 2657 integer,intent(in) :: ngfft(18) 2658 real(dp),intent(in) :: gmet(3,3),rhog(2,nfft) 2659 real(dp),intent(out) :: vqgradhart(cplex*nfft) 2660 2661 !Local variables------------------------------- 2662 !scalars 2663 integer,parameter :: im=2,re=1 2664 integer :: i1,i2,i23,i2_local,i3 2665 integer :: id1,id2,id3,ig1,ig2,ig3,ii,ii1,me_fft,n1,n2,n3,nproc_fft 2666 real(dp) :: cutoff,gfact,gnorm,num 2667 real(dp), parameter :: piinv2= piinv*two 2668 real(dp),parameter :: tolfix=1.000000001e0_dp 2669 !arrays 2670 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 2671 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 2672 real(dp),allocatable :: work1(:,:) 2673 real(dp) :: gvec(3) 2674 2675 ! ************************************************************************* 2676 2677 DBG_ENTER("COLL") 2678 2679 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 2680 nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft 2681 2682 !Get the distrib associated with this fft_grid 2683 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 2684 2685 !Initialize a few quantities 2686 cutoff=gsqcut*tolfix 2687 ABI_MALLOC(work1,(2,nfft)) 2688 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 2689 2690 !Triple loop on each dimension 2691 do i3=1,n3 2692 ig3=i3-(i3/id3)*n3-1 2693 2694 do i2=1,n2 2695 ig2=i2-(i2/id2)*n2-1 2696 2697 if (fftn2_distrib(i2) == me_fft) then 2698 i2_local = ffti2_local(i2) 2699 i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1)) 2700 !Do the test that eliminates the Gamma point outside of the inner loop 2701 ii1=1 2702 if(i23==0 .and. ig2==0 .and. ig3==0)then 2703 ii1=2 2704 work1(re,1+i23)=zero 2705 work1(im,1+i23)=zero 2706 end if 2707 2708 ! Final inner loop on the first dimension (note the lower limit) 2709 do i1=ii1,n1 2710 ig1=i1-(i1/id1)*n1-1 2711 ii=i1+i23 2712 2713 gvec=(/ig1,ig2,ig3/) 2714 gnorm=normv(gvec,gmet,'r') !'r' is to avoid the 2pi scalling 2715 2716 if (gnorm**2<=cutoff) then 2717 num=dot_product(gmet(qdir,:),gvec(:)) 2718 gfact=piinv2*num/gnorm**4 2719 work1(re,ii)=-rhog(re,ii)*gfact 2720 work1(im,ii)=-rhog(im,ii)*gfact 2721 else 2722 work1(re,ii)=zero 2723 work1(im,ii)=zero 2724 end if 2725 2726 end do ! End loop on i1 2727 end if 2728 2729 end do ! End loop on i2 2730 end do ! End loop on i3 2731 2732 ! Fourier Transform the q-gradient of the hartree potential, in reciprocal space it was stored in work1 2733 call fourdp(cplex,work1,vqgradhart,1,mpi_enreg,nfft,1,ngfft,0) 2734 2735 ABI_FREE(work1) 2736 2737 DBG_EXIT("COLL") 2738 2739 end subroutine hartredq
ABINIT/m_spacepar [ Modules ]
NAME
m_spacepar
FUNCTION
Relatively Low-level procedures operating on arrays defined on the FFT box (G- or R- space) Unlike the procedures in m_cgtools, the routines declared in this module can use mpi_type.
COPYRIGHT
Copyright (C) 2008-2022 ABINIT group (XG, BA, MT, DRH, DCA, GMR, MJV, 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_spacepar 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_xmpi 29 use m_sort 30 31 use m_time, only : timab 32 use defs_abitypes, only : MPI_type 33 use m_symtk, only : mati3inv, chkgrp, symdet, symatm, matr3inv 34 use m_geometry, only : metric, normv, symredcart,wedge_basis,wedge_product 35 use m_gtermcutoff, only : termcutoff 36 use m_mpinfo, only : ptabs_fourdp 37 use m_fft, only : zerosym, fourdp 38 39 implicit none 40 41 private
m_spacepar/hartre [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
hartre
FUNCTION
Given rho(G), compute Hartree potential (=FFT of rho(G)/pi/(G+q)**2) When cplex=1, assume q=(0 0 0), and vhartr will be REAL When cplex=2, q must be taken into account, and vhartr will be COMPLEX
NOTES
*Modified code to avoid if statements inside loops to skip G=0. Replaced if statement on G^2>gsqcut to skip G s outside where rho(G) should be 0. Effect is negligible but gsqcut should be used to be strictly consistent with usage elsewhere in code. *The speed-up is provided by doing a few precomputations outside the inner loop. One variable size array is needed for this (gq).
INPUTS
cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX gsqcut=cutoff value on G**2 for sphere inside fft box. (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)) izero=if 1, unbalanced components of Vhartree(g) are set to zero mpi_enreg=information about MPI parallelization 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 [qpt(3)=reduced coordinates for a wavevector to be combined with the G vectors (needed if cplex==2).] rhog(2,nfft)=electron density in G space rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
vhartr(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX
SOURCE
547 subroutine hartre(cplex,gsqcut,icutcoul,izero,mpi_enreg,nfft,ngfft,nkpt,& 548 &rcut,rhog,rprimd,vcutgeo,vhartr,& 549 &qpt) ! Optional arguments 550 551 !Arguments ------------------------------------ 552 !scalars 553 integer,intent(in) :: cplex,icutcoul,izero,nfft,nkpt 554 real(dp),intent(in) :: gsqcut,rcut 555 type(MPI_type),intent(in) :: mpi_enreg 556 !arrays 557 integer,intent(in) :: ngfft(18) 558 real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft),vcutgeo(3) 559 real(dp),intent(in),optional :: qpt(3) 560 real(dp),intent(out) :: vhartr(cplex*nfft) 561 562 !Local variables------------------------------- 563 !scalars 564 integer,parameter :: im=2,re=1 565 integer :: i1,i2,i23,i2_local,i3,id1,id2,id3 566 integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max 567 integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05,me_fft,nproc_fft 568 real(dp),parameter :: tolfix=1.000000001e0_dp 569 real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,ucvol 570 character(len=500) :: message 571 !arrays 572 integer :: id(3) 573 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 574 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 575 real(dp) :: gmet(3,3),gprimd(3,3),qpt_(3),rmet(3,3),tsec(2) 576 real(dp),allocatable :: gcutoff(:) 577 real(dp),allocatable :: gq(:,:),work1(:,:) 578 579 ! ************************************************************************* 580 581 ! Keep track of total time spent in hartre 582 call timab(10,1,tsec) 583 584 ! Check that cplex has an allowed value 585 if(cplex/=1 .and. cplex/=2)then 586 write(message, '(a,i0,a,a)' )& 587 'From the calling routine, cplex=',cplex,ch10,& 588 'but the only value allowed are 1 and 2.' 589 ABI_BUG(message) 590 end if 591 592 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 593 594 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 595 nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft 596 597 ! Get the distrib associated with this fft_grid 598 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 599 600 ! Initialize a few quantities 601 cutoff=gsqcut*tolfix 602 if(present(qpt))then 603 qpt_=qpt 604 else 605 qpt_=zero 606 end if 607 qeq0=0 608 if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1 609 qeq05=0 610 if (qeq0==0) then 611 if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or.abs(abs(qpt_(3))-half)<tol12) qeq05=1 612 end if 613 614 ! If cplex=1 then qpt_ should be 0 0 0 615 if (cplex==1.and. qeq0/=1) then 616 write(message,'(a,3e12.4,a,a)')& 617 'cplex=1 but qpt=',qpt_,ch10,& 618 'qpt should be 0 0 0.' 619 ABI_BUG(message) 620 end if 621 622 ! If FFT parallelism then qpt should not be 1/2 623 if (nproc_fft>1.and.qeq05==1) then 624 write(message, '(a,3e12.4,a,a)' )& 625 'FFT parallelism selected but qpt',qpt_,ch10,& 626 'qpt(i) should not be 1/2...' 627 ABI_ERROR(message) 628 end if 629 630 !Initialize Gcut-off array from m_gtermcutoff 631 !ABI_MALLOC(gcutoff,(ngfft(1)*ngfft(2)*ngfft(3))) 632 call termcutoff(gcutoff,gsqcut,icutcoul,ngfft,nkpt,rcut,rprimd,vcutgeo) 633 634 ! In order to speed the routine, precompute the components of g+q 635 ! Also check if the booked space was large enough... 636 ABI_MALLOC(gq,(3,max(n1,n2,n3))) 637 do ii=1,3 638 id(ii)=ngfft(ii)/2+2 639 do ing=1,ngfft(ii) 640 ig=ing-(ing/id(ii))*ngfft(ii)-1 641 gq(ii,ing)=ig+qpt_(ii) 642 end do 643 end do 644 ig1max=-1;ig2max=-1;ig3max=-1 645 ig1min=n1;ig2min=n2;ig3min=n3 646 647 ABI_MALLOC(work1,(2,nfft)) 648 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 649 650 ! Triple loop on each dimension 651 do i3=1,n3 652 ig3=i3-(i3/id3)*n3-1 653 ! Precompute some products that do not depend on i2 and i1 654 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 655 gqgm23=gq(3,i3)*gmet(2,3)*2 656 gqgm13=gq(3,i3)*gmet(1,3)*2 657 658 do i2=1,n2 659 ig2=i2-(i2/id2)*n2-1 660 if (fftn2_distrib(i2) == me_fft) then 661 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 662 gqgm12=gq(2,i2)*gmet(1,2)*2 663 gqg2p3=gqgm13+gqgm12 664 665 i2_local = ffti2_local(i2) 666 i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1)) 667 ! Do the test that eliminates the Gamma point outside of the inner loop 668 ii1=1 669 if(i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0)then 670 ii1=2 671 work1(re,1+i23)=zero 672 work1(im,1+i23)=zero 673 end if 674 675 ! Final inner loop on the first dimension (note the lower limit) 676 do i1=ii1,n1 677 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 678 ii=i1+i23 679 680 if(gs<=cutoff)then 681 ! Identify min/max indexes (to cancel unbalanced contributions later) 682 ! Count (q+g)-vectors with similar norm 683 if ((qeq05==1).and.(izero==1)) then 684 ig1=i1-(i1/id1)*n1-1 685 ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1) 686 ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2) 687 ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3) 688 end if 689 690 den=piinv/gs*gcutoff(ii) 691 work1(re,ii)=rhog(re,ii)*den 692 work1(im,ii)=rhog(im,ii)*den 693 else 694 ! gs>cutoff 695 work1(re,ii)=zero 696 work1(im,ii)=zero 697 end if 698 699 end do ! End loop on i1 700 end if 701 end do ! End loop on i2 702 end do ! End loop on i3 703 704 ABI_FREE(gq) 705 706 if (izero==1) then 707 ! Set contribution of unbalanced components to zero 708 709 if (qeq0==1) then !q=0 710 call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 711 712 else if (qeq05==1) then 713 !q=1/2; this doesn't work in parallel 714 ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2 715 ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2 716 ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2 717 if (abs(abs(qpt_(1))-half)<tol12) then 718 if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max) 719 if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min) 720 end if 721 if (abs(abs(qpt_(2))-half)<tol12) then 722 if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max) 723 if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min) 724 end if 725 if (abs(abs(qpt_(3))-half)<tol12) then 726 if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max) 727 if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min) 728 end if 729 call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,& 730 comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 731 end if 732 end if 733 734 ! Fourier Transform Vhartree. Vh in reciprocal space was stored in work1 735 call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,1,ngfft,0) 736 737 ABI_FREE(gcutoff) 738 ABI_FREE(work1) 739 740 call timab(10,2,tsec) 741 742 end subroutine hartre
m_spacepar/hartrestr [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
hartrestr
FUNCTION
To be called for strain perturbation only Compute the inhomogenous terms generated by the strain derivative of Hartree potential due to the ground state charge rho(G) FFT of (rho(G)/pi)*[d(1/G**2)/d(strain) - delta(diagonal strain)*(1/G**2)]
NOTES
*based largely on hartre.f *Modified code to avoid if statements inside loops to skip G=0. Replaced if statement on G^2>gsqcut to skip G s outside where rho(G) should be 0. Effect is negligible but gsqcut should be used to be strictly consistent with usage elsewhere in code. *The speed-up is provided by doing a few precomputations outside the inner loop. One variable size array is needed for this (gq).
INPUTS
gsqcut=cutoff value on G**2 for sphere inside fft box. idir=direction of the current perturbation ipert=type of the perturbation mpi_enreg=information about MPI parallelization natom=number of atoms in cell. nfft=number of fft grid points (gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2)) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft rhog(2,nfft)=array for Fourier transform of GS electron density rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
vhartr1(nfft)=Inhomogeneous term in strain-perturbation-induced Hartree potential in real space,
SOURCE
1298 subroutine hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,rhog,rprimd,vhartr1) 1299 1300 !Arguments ------------------------------------ 1301 !scalars 1302 integer,intent(in) :: idir,ipert,natom,nfft 1303 real(dp),intent(in) :: gsqcut 1304 type(MPI_type),intent(in) :: mpi_enreg 1305 !arrays 1306 integer,intent(in) :: ngfft(18) 1307 real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3) 1308 real(dp),intent(out) :: vhartr1(nfft) 1309 1310 !Local variables------------------------------- 1311 !scalars 1312 integer,parameter :: im=2,re=1 1313 integer :: i1,i2,i23,i3,id2,id3,ig,ig2,ig3,ii,ii1,ing,istr,ka,kb,n1,n2,n3 1314 real(dp),parameter :: tolfix=1.000000001_dp 1315 real(dp) :: cutoff,ddends,den,dgsds,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3 1316 real(dp) :: term,ucvol 1317 character(len=500) :: message 1318 !arrays 1319 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 1320 integer :: id(3) 1321 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1322 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1323 real(dp) :: dgmetds(3,3),gmet(3,3),gprimd(3,3),gqr(3),rmet(3,3) 1324 real(dp),allocatable :: gq(:,:),work1(:,:) 1325 1326 ! ************************************************************************* 1327 1328 if( .not. (ipert==natom+3 .or. ipert==natom+4))then 1329 write(message, '(a,i0,a,a)' )& 1330 & 'From the calling routine, ipert=',ipert,ch10,& 1331 & 'so this routine for the strain perturbation should not be called.' 1332 ABI_BUG(message) 1333 end if 1334 1335 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1336 1337 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1338 1339 !Get the distrib associated with this fft_grid 1340 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1341 1342 !Initialize a few quantities 1343 cutoff=gsqcut*tolfix 1344 1345 istr=idir + 3*(ipert-natom-3) 1346 1347 if(istr<1 .or. istr>6)then 1348 write(message, '(a,i10,a,a,a)' )& 1349 & 'Input dir gives istr=',istr,' not allowed.',ch10,& 1350 & 'Possible values are 1,2,3,4,5,6 only.' 1351 ABI_BUG(message) 1352 end if 1353 1354 ka=idx(2*istr-1);kb=idx(2*istr) 1355 do ii = 1,3 1356 dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii)) 1357 end do 1358 !For historical reasons: 1359 dgmetds(:,:)=0.5_dp*dgmetds(:,:) 1360 1361 !In order to speed the routine, precompute the components of g+q 1362 !Also check if the booked space was large enough... 1363 ABI_MALLOC(gq,(3,max(n1,n2,n3))) 1364 do ii=1,3 1365 id(ii)=ngfft(ii)/2+2 1366 do ing=1,ngfft(ii) 1367 ig=ing-(ing/id(ii))*ngfft(ii)-1 1368 gq(ii,ing)=ig 1369 end do 1370 end do 1371 1372 ABI_MALLOC(work1,(2,nfft)) 1373 id2=n2/2+2 1374 id3=n3/2+2 1375 !Triple loop on each dimension 1376 do i3=1,n3 1377 ig3=i3-(i3/id3)*n3-1 1378 ! Precompute some products that do not depend on i2 and i1 1379 gqr(3)=gq(3,i3) 1380 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 1381 gqgm23=gq(3,i3)*gmet(2,3)*2 1382 gqgm13=gq(3,i3)*gmet(1,3)*2 1383 1384 do i2=1,n2 1385 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 1386 gqr(2)=gq(2,i2) 1387 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 1388 gqgm12=gq(2,i2)*gmet(1,2)*2 1389 gqg2p3=gqgm13+gqgm12 1390 ig2=i2-(i2/id2)*n2-1 1391 ! i23=n1*((i2-1)+n2*(i3-1)) 1392 i23=n1*((ffti2_local(i2)-1)+(n2/mpi_enreg%nproc_fft)*(i3-1)) 1393 ! Do the test that eliminates the Gamma point outside 1394 ! of the inner loop 1395 ii1=1 1396 if(i23==0 .and. ig2==0 .and. ig3==0)then 1397 ii1=2 1398 work1(re,1+i23)=0.0_dp 1399 work1(im,1+i23)=0.0_dp 1400 end if 1401 1402 ! Final inner loop on the first dimension 1403 ! (note the lower limit) 1404 do i1=ii1,n1 1405 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 1406 ii=i1+i23 1407 if(gs<=cutoff)then 1408 den=piinv/gs 1409 gqr(1)=gq(1,i1) 1410 dgsds=& 1411 & (gqr(1)*(dgmetds(1,1)*gqr(1)+dgmetds(1,2)*gqr(2)+dgmetds(1,3)*gqr(3))+ & 1412 & gqr(2)*(dgmetds(2,1)*gqr(1)+dgmetds(2,2)*gqr(2)+dgmetds(2,3)*gqr(3))+ & 1413 & gqr(3)*(dgmetds(3,1)*gqr(1)+dgmetds(3,2)*gqr(2)+dgmetds(3,3)*gqr(3)) ) 1414 ddends=-piinv*dgsds/gs**2 1415 if(istr<=3)then 1416 term=2.0_dp*ddends-den 1417 else 1418 term=2.0_dp*ddends 1419 end if 1420 work1(re,ii)=rhog(re,ii)*term 1421 work1(im,ii)=rhog(im,ii)*term 1422 else 1423 work1(re,ii)=0.0_dp 1424 work1(im,ii)=0.0_dp 1425 end if 1426 1427 end do ! End loop on i1 1428 end if 1429 end do ! End loop on i2 1430 end do ! End loop on i3 1431 1432 ABI_FREE(gq) 1433 1434 !Fourier Transform Vhartree. 1435 !Vh in reciprocal space was stored in work1 1436 call fourdp(1,work1,vhartr1,1,mpi_enreg,nfft,1,ngfft,0) 1437 1438 ABI_FREE(work1) 1439 1440 end subroutine hartrestr
m_spacepar/irrzg [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
irrzg
FUNCTION
Find the irreducible zone in reciprocal space under the symmetry group with real space rotations in symrel(3,3,nsym). The (integer) rotation matrices symrel(3,3,nsym) express the new real space positions (e.g. rotated atom positions) in REDUCED coordinates, i.e. in coordinates expressed as fractions of real space primitive translations (atomic coordinates xred). tnons(3,nsym) express the associated nonsymmorphic translations, again in reduced coordinates. Special data structure created in irrzon. First half holds mapping from irr zone to full zone; part of second half holds repetition number info. work1 is a work array to keep track of grid points found so far. In case nspden=2 and nsppol=1, one has to take care of antiferromagnetic operations. The subgroup of non-magnetic operations is used to generate irrzon(:,:,2) and phnons(:,:,2), while the full group is used to generate irrzon(:,:,1) and phnons(:,:,1)
NOTES
for reference in the near future (2018), some notes: this routine should be duplicated for magnetizations in spinorial formalism. The only main difference will be that the density is not simply transported to the image point under symrel but is a vector which has to be transformed as well $S \vec{m} (\vec{x}) = \sigma \vec{m} (S \vec{x} + \tau)$ $S \vec{m} (\vec{G}) = \sigma exp(+ 2 \pi i \vec{G} \vec{tau}) \vec{m}(S^{-1 T} \vec{G})$ S is a symop, sigma the AFM sign flip if any, tau the partial non symmorphic translation x a position, m the magnetization 3 vector The phase factor is the same as below for the density, but the collection of elements which are equal is more complex: for a 3 or 6 axis the m vector could transform one component to a linear combination of several others (I think). Things are not necessarily aligned with the z axis.
INPUTS
nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetry elements in group n1,n2,n3=box dimensions of real space grid (or fft box) symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry matrices in real space (integers) tnons(3,nsym)=reduced nonsymmorphic translations (symrel and tnons are in terms of real space primitive translations)
OUTPUT
irrzon(n1*n2*n3,2+(nspden/4),(nspden/nsppol)-3*(nspden/4))=integer array which contains the locations of related grid points and the repetition number for each symmetry class. phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4))=phases associated with nonsymmorphic translations
SOURCE
2002 subroutine irrzg(irrzon,nspden,nsppol,nsym,n1,n2,n3,phnons,symafm,symrel,tnons) 2003 2004 !Arguments ------------------------------------ 2005 !scalars 2006 integer,intent(in) :: n1,n2,n3,nspden,nsppol,nsym 2007 !arrays 2008 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 2009 integer,intent(out) :: irrzon(n1*n2*n3,2,(nspden/nsppol)-3*(nspden/4)) 2010 real(dp),intent(in) :: tnons(3,nsym) 2011 real(dp),intent(out) :: phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4)) 2012 2013 !Local variables------------------------------- 2014 !scalars 2015 integer :: i1,i2,i3,id1,id2,id3,ifft,imagn,ind1,ind2,ipt,irep,isym,izone 2016 integer :: izonemax,j1,j2,j3,jj,k1,k2,k3,l1,l2,l3,nfftot,npt,nsym_used 2017 integer :: nzone,setzer,sppoldbl 2018 real(dp) :: arg,ph1i,ph1r,ph2i,ph2r,tau1,tau2,tau3 2019 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism 2020 character(len=500) :: message 2021 !arrays 2022 integer,allocatable :: class(:),iperm(:),symafm_used(:),symrel_used(:,:,:) 2023 integer,allocatable :: work1(:) 2024 real(dp),allocatable :: tnons_used(:,:),work2(:,:) 2025 2026 ! ************************************************************************* 2027 2028 ABI_MALLOC(class,(nsym)) 2029 ABI_MALLOC(iperm,(nsym)) 2030 ABI_MALLOC(work1,(n1*n2*n3)) 2031 ABI_MALLOC(work2,(2,n1*n2*n3)) 2032 2033 nfftot=n1*n2*n3 2034 2035 id1=n1/2+2 2036 id2=n2/2+2 2037 id3=n3/2+2 2038 2039 sppoldbl=nspden/nsppol;if (nspden==4) sppoldbl=1 2040 2041 do imagn=1,sppoldbl 2042 2043 ! Treat in a similar way the case of the full group and the non-magnetic subgroup 2044 nsym_used=0 2045 do isym=1,nsym 2046 if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. & 2047 & ((nspden==4).and.afm_noncoll) )then 2048 nsym_used=nsym_used+1 2049 end if 2050 end do 2051 2052 if(imagn==2 .and. nsym_used/=nsym/2)then 2053 write(message, '(a,a,a,a,a,i4,a,i0)' )& 2054 & ' The number of ferromagnetic symmetry operations must be',ch10,& 2055 & ' half the total number of operations, while it is observed that',ch10,& 2056 & ' nsym=',nsym,' and nsym_magn=',nsym_used 2057 ABI_BUG(message) 2058 end if 2059 2060 ABI_MALLOC(symafm_used,(nsym_used)) 2061 ABI_MALLOC(symrel_used,(3,3,nsym_used)) 2062 ABI_MALLOC(tnons_used,(3,nsym_used)) 2063 2064 nsym_used=0 2065 do isym=1,nsym 2066 if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. & 2067 & ((nspden==4).and.afm_noncoll) ) then 2068 nsym_used=nsym_used+1 2069 symrel_used(:,:,nsym_used)=symrel(:,:,isym) 2070 tnons_used(:,nsym_used)=tnons(:,isym) 2071 symafm_used(nsym_used)=symafm(isym) 2072 end if 2073 end do 2074 if ((nspden/=4).or.(.not.afm_noncoll)) symafm_used=1 2075 2076 2077 ! Zero out work array--later on, a zero entry will mean that 2078 ! a given grid point has not yet been assigned to an ibz point 2079 work1(1:nfftot)=0 2080 irrzon(:,2,imagn)=0 2081 2082 ! Initialize at G=0 (always in irreducible zone) 2083 nzone=1 2084 irrzon(1,1,imagn)=1 2085 irrzon(1,2,imagn)=nsym_used 2086 ! Set phase exp(2*Pi*I*G dot tnons) for G=0 2087 phnons(1,1,imagn)=one 2088 phnons(2,1,imagn)=zero 2089 npt=1 2090 ! setting work1(1)=1 indicates that first grid point (G=0) is 2091 ! in the iz (irreducible zone) 2092 work1(1)=1 2093 2094 ind1=0 2095 2096 ! Loop over reciprocal space grid points: 2097 do i3=1,n3 2098 do i2=1,n2 2099 do i1=1,n1 2100 2101 ind1=ind1+1 2102 2103 ! Check to see whether present grid point is equivalent to 2104 ! any previously identified ibz point--if not, a new ibz point 2105 ! has been found 2106 2107 if (work1(ind1)==0) then 2108 2109 ! A new class has been found. 2110 2111 ! Get location of G vector (grid point) centered at 0 0 0 2112 l3=i3-(i3/id3)*n3-1 2113 l2=i2-(i2/id2)*n2-1 2114 l1=i1-(i1/id1)*n1-1 2115 2116 do isym=1,nsym_used 2117 2118 ! Get rotated G vector Gj for each symmetry element 2119 ! -- here we use the TRANSPOSE of symrel_used; assuming symrel_used expresses 2120 ! the rotation in real space, the transpose is then appropriate 2121 ! for G space symmetrization (p. 1172d,e of notes, 2 June 1995). 2122 j1=symrel_used(1,1,isym)*l1+& 2123 & symrel_used(2,1,isym)*l2+symrel_used(3,1,isym)*l3 2124 j2=symrel_used(1,2,isym)*l1+& 2125 & symrel_used(2,2,isym)*l2+symrel_used(3,2,isym)*l3 2126 j3=symrel_used(1,3,isym)*l1+& 2127 & symrel_used(2,3,isym)*l2+symrel_used(3,3,isym)*l3 2128 2129 ! Map into [0,n-1] 2130 k1=mod(n1+mod(j1,n1),n1) 2131 k2=mod(n2+mod(j2,n2),n2) 2132 k3=mod(n3+mod(j3,n3),n3) 2133 2134 ! Get linear index of rotated point Gj 2135 ind2=1+k1+n1*(k2+n2*k3) 2136 2137 ! Store info for new class: 2138 class(isym)=ind2 2139 iperm(isym)=isym 2140 2141 ! Setting work array element to 1 indicates grid point has been 2142 ! identified with iz point 2143 work1(ind2)=1 2144 2145 ! End of loop on isym 2146 end do 2147 2148 ! Sort integers into ascending order in each class 2149 ! (this lumps together G vectors with the same linear index, i.e. 2150 ! groups together symmetries which land on the same Gj) 2151 call sort_int(nsym_used,class,iperm) 2152 2153 ! Check repetition factor (how many identical copies of Gj occur 2154 ! from all symmetries applied to G) 2155 irep=0 2156 do isym=1,nsym_used 2157 if (class(isym)==class(1)) then 2158 irep=irep+1 2159 end if 2160 end do 2161 ipt=nsym_used/irep 2162 2163 ! Repetition factor must be divisor of nsym_used: 2164 if (nsym_used/=(ipt*irep)) then 2165 write(message, '(a,i5,a,i6,a,a,a,a,a,a)' )& 2166 & ' irep=',irep,' not a divisor of nsym_used=',nsym_used,ch10,& 2167 & ' This usually indicates that',& 2168 & ' the input symmetries do not form a group.',ch10,& 2169 & ' Action : check the input symmetries carefully do they',& 2170 & ' form a group ? If they do, there is a code bug.' 2171 ABI_ERROR(message) 2172 end if 2173 2174 ! Compute phases for any nonsymmorphic symmetries 2175 ! exp(-2*Pi*I*G dot tau(j)) for each symmetry j with 2176 ! (possibly zero) nonsymmorphic translation tau(j) 2177 do jj=1,nsym_used 2178 ! First get nonsymmorphic translation and see if nonzero 2179 ! (iperm grabs the symmetries in the new order after sorting) 2180 isym=iperm(jj) 2181 tau1=tnons_used(1,isym) 2182 tau2=tnons_used(2,isym) 2183 tau3=tnons_used(3,isym) 2184 if (abs(tau1)>tol12.or.abs(tau2)>tol12& 2185 & .or.abs(tau3)>tol12) then 2186 ! compute exp(-2*Pi*I*G dot tau) using original G 2187 arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3) 2188 work2(1,jj)=cos(arg) 2189 work2(2,jj)=-sin(arg) 2190 else 2191 work2(1,jj)=one 2192 work2(2,jj)=zero 2193 end if 2194 end do 2195 2196 ! All phases arising from symmetries which map to the same 2197 ! G vector must actually be the same because 2198 ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G) 2199 ! must be satisfied; if exp(2*Pi*I*(G) dot tau_S) can be different 2200 ! for two different symmetries S which both take G to the same St*G, 2201 ! then the related Fourier components rho(St*G) must VANISH. 2202 ! Thus: set "phase" to ZERO here in that case. 2203 ! The G mappings occur in sets of irep members; if irep=1 then 2204 ! all the G are unique. 2205 ! MT 090212: 2206 ! In the case of antiferromagn. symetries, one can have 2207 ! rho(Strans*G)= -exp(2*Pi*I*(G) dot tau_S) rho(G) 2208 ! (look at the minus !) 2209 ! A special treatment is then operated on phons. 2210 ! The later must be consistent with the use of phnons array 2211 ! in symrhg.F90 routine. 2212 ! XG 001108 : 2213 ! Note that there is a tolerance on the 2214 ! accuracy of tnons, especially when they are found from 2215 ! the symmetry finder (with xred that might be a bit inaccurate) 2216 if (irep > 1) then 2217 do jj=1,nsym_used,irep 2218 setzer=0 2219 ph1r=work2(1,jj);ph1i=work2(2,jj) 2220 do j1=jj,jj+irep-1 2221 ph2r=work2(1,j1);ph2i=work2(2,j1) 2222 if (((ph2r+ph1r)**2+(ph2i+ph1i)**2) <= tol14) then 2223 if (setzer/=1) setzer=-1 2224 else if (((ph2r-ph1r)**2+(ph2i-ph1i)**2) > tol14) then 2225 setzer=1 2226 end if 2227 end do 2228 ! Setzer= 0: phnons are all equal 2229 ! Setzer=-1: phnons are equal in absolute value 2230 ! Setzer= 1: some phnons are different 2231 if (setzer/=0) then 2232 if (setzer==-1) then 2233 if (afm_noncoll.and.nspden==4) then 2234 arg=symafm_used(iperm(jj)) 2235 if (all(symafm_used(iperm(jj:jj+irep-1))==arg)) then 2236 setzer=1 2237 else 2238 do j1=jj,jj+irep-1 2239 work2(:,j1)=work2(:,j1)*dble(symafm_used(iperm(j1))) 2240 end do 2241 end if 2242 else 2243 setzer=1 2244 end if 2245 end if 2246 if (setzer==1) work2(:,jj:jj+irep-1)=zero 2247 end if 2248 end do 2249 ! Compress data if irep>1: 2250 jj=0 2251 do isym=1,nsym_used,irep 2252 jj=jj+1 2253 class(jj)=class(isym) 2254 work2(1,jj)=work2(1,isym) 2255 work2(2,jj)=work2(2,isym) 2256 end do 2257 end if 2258 2259 ! Put new unique points into irrzon array: 2260 irrzon(1+npt:ipt+npt,1,imagn)=class(1:ipt) 2261 2262 ! Put repetition number into irrzon array: 2263 irrzon(1+nzone,2,imagn)=irep 2264 2265 ! DEBUG 2266 ! write(std_out,'(a,6i7)' )' irrzg : izone,i1,i2,i3,imagn,irrzon(859,2,1)=',& 2267 ! & 1+nzone,i1,i2,i3,imagn,irrzon(859,2,1) 2268 ! ENDDEBUG 2269 2270 ! Put phases (or 0) in phnons array: 2271 phnons(:,1+npt:ipt+npt,imagn)=work2(:,1:ipt) 2272 2273 ! Update number of points in irrzon array: 2274 ! (irep must divide evenly into nsym_used !) 2275 npt=npt+ipt 2276 2277 ! Update number of classes: 2278 nzone=nzone+1 2279 2280 end if 2281 ! 2282 ! End of loop on reciprocal space points, with indices i1, i2, i3 2283 end do 2284 end do 2285 end do 2286 2287 ABI_SFREE(symafm_used) 2288 ABI_SFREE(symrel_used) 2289 ABI_SFREE(tnons_used) 2290 2291 end do ! imagn 2292 2293 !Make sure number of real space points accounted for equals actual number of grid points 2294 if (npt/=n1*n2*n3) then 2295 write(message, '(a,a,a,a,i10,a,i10,a,a,a,a,a,a,a,a,a)' ) ch10,& 2296 & ' irrzg : ERROR -',ch10,& 2297 & ' npt=',npt,' and n1*n2*n3=',n1*n2*n3,' are not equal',ch10,& 2298 & ' This says that the total of all points in the irreducible',& 2299 & ' sector in real space',ch10,& 2300 & ' and all symmetrically equivalent',& 2301 & ' points, npt, does not equal the actual number',ch10,& 2302 & ' of real space grid points.' 2303 call wrtout(std_out,message,'COLL') 2304 write(message,'(3a)') & 2305 & ' This may mean that the input symmetries do not form a group',ch10,& 2306 & ' Action : check input symmetries carefully for errors.' 2307 ABI_ERROR(message) 2308 end if 2309 2310 !Perform some checks 2311 do imagn=1,sppoldbl 2312 2313 do ifft=1,nfftot 2314 if (irrzon(ifft,1,imagn)<1.or.irrzon(ifft,1,imagn)>nfftot) then 2315 write(message,'(a,4i0,a,a)')& 2316 & ' ifft,irrzon(ifft,1,imagn),nfftot,imagn=',ifft,irrzon(ifft,1,imagn),nfftot,imagn,ch10,& 2317 & ' =>irrzon goes outside acceptable bounds.' 2318 ABI_BUG(message) 2319 end if 2320 end do 2321 2322 izonemax=0 2323 do izone=1,nfftot 2324 ! Global bounds 2325 if (irrzon(izone,2,imagn)<0.or.irrzon(izone,2,imagn)>(nsym/imagn)) then 2326 write(message, '(a,5i7,a,a)' )& 2327 & ' izone,nzone,irrzon(izone,2,imagn),nsym,imagn =',izone,nzone,irrzon(izone,2,imagn),nsym,imagn,ch10,& 2328 & ' =>irrzon goes outside acceptable bounds.' 2329 ABI_BUG(message) 2330 end if 2331 ! Second index only goes up to nzone 2332 if(izonemax==0)then 2333 if (irrzon(izone,2,imagn)==0)izonemax=izone-1 2334 end if 2335 if(izonemax/=0)then 2336 if (irrzon(izone,2,imagn)/=0) then 2337 message = ' beyond izonemax, irrzon(izone,2,imagn) should be zero' 2338 ABI_BUG(message) 2339 end if 2340 end if 2341 end do 2342 2343 end do ! imagn 2344 2345 ABI_FREE(class) 2346 ABI_FREE(iperm) 2347 ABI_FREE(work1) 2348 ABI_FREE(work2) 2349 2350 end subroutine irrzg
m_spacepar/laplacian [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
laplacian
FUNCTION
compute the laplacian of a function defined in real space the code is written in the way of /3xc/xcden.F90
INPUTS
gprimd(3,3)=dimensional reciprocal space primitive translations mpi_enreg=information about MPI parallelization nfft=number of points of the fft grid nfunc=number of functions on the grid for which the laplacian is to be calculated ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft (optional) rdfuncr(nfft,nfunc)=real(dp) discretized functions in real space rdfuncg_in TO BE DESCRIBED SB 090901 laplacerdfuncg_in TO BE DESCRIBED SB 090901 (optional) g2cart_in(nfft) = G**2 on the grid
OUTPUT
(optional) laplacerdfuncr = laplacian in real space of the functions in rdfuncr (optional) rdfuncg = real(dp) discretized functions in fourier space (optional) laplacerdfuncg = real(dp) discretized laplacian of the functions in fourier space (optional) g2cart_out(nfft) = G**2 on the grid rdfuncg_out TO BE DESCRIBED SB 090901 laplacerdfuncg_out TO BE DESCRIBED SB 090901
SOURCE
982 subroutine laplacian(gprimd,mpi_enreg,nfft,nfunc,ngfft,rdfuncr,& 983 & laplacerdfuncr,rdfuncg_out,laplacerdfuncg_out,g2cart_out,rdfuncg_in,g2cart_in) 984 985 !Arguments ------------------------------------ 986 !scalars 987 integer,intent(in) :: nfft,nfunc 988 type(MPI_type),intent(in) :: mpi_enreg 989 !arrays 990 integer,intent(in) :: ngfft(18) 991 real(dp),intent(in) :: gprimd(3,3) 992 real(dp),intent(inout),optional :: laplacerdfuncr(nfft,nfunc) 993 real(dp),intent(inout),optional,target :: rdfuncr(nfft,nfunc) 994 real(dp),intent(in),optional,target :: g2cart_in(nfft) !vz_i 995 real(dp),intent(out),optional,target :: g2cart_out(nfft) !vz_i 996 real(dp),intent(out),optional,target :: laplacerdfuncg_out(2,nfft,nfunc) 997 real(dp),intent(in),optional,target :: rdfuncg_in(2,nfft,nfunc) !vz_i 998 real(dp),intent(out),optional,target :: rdfuncg_out(2,nfft,nfunc) 999 1000 !Local variables------------------------------- 1001 !scalars 1002 integer :: count,i1,i2,i3,id1,id2,id3,ifft,ifunc,ig1,ig2,ig3,ii1,n1,n2 1003 integer :: n3 1004 real(dp) :: b11,b12,b13,b21,b22,b23,b31,b32,b33 1005 !arrays 1006 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1007 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1008 real(dp),ABI_CONTIGUOUS pointer :: g2cart(:),laplacerdfuncg(:,:,:),rdfuncg(:,:,:) 1009 1010 ! ************************************************************************* 1011 1012 !Keep local copy of fft dimensions 1013 n1=ngfft(1) 1014 n2=ngfft(2) 1015 n3=ngfft(3) 1016 1017 if(present(laplacerdfuncg_out)) then 1018 laplacerdfuncg => laplacerdfuncg_out 1019 else 1020 ABI_MALLOC(laplacerdfuncg,(2,nfft,nfunc)) 1021 end if 1022 1023 ! Get the distrib associated with this fft_grid 1024 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1025 1026 !change the real density rdfuncr on real space on the real density 1027 !rdfuncg in reciprocal space 1028 if(.not.present(rdfuncg_in)) then 1029 if(present(rdfuncg_out)) then 1030 rdfuncg => rdfuncg_out 1031 else 1032 ABI_MALLOC(rdfuncg,(2,nfft,nfunc)) 1033 end if 1034 if(present(rdfuncr)) then 1035 do ifunc=1,nfunc 1036 call fourdp(1,rdfuncg(:,:,ifunc),rdfuncr(:,ifunc),-1,mpi_enreg,nfft,1,ngfft,0) 1037 end do 1038 end if 1039 else 1040 rdfuncg => rdfuncg_in 1041 end if 1042 1043 !apply the laplacian on laplacerdfuncr 1044 !code from /3xc/xcden.F90 1045 !see also src/5common/hatre.F90 and src/5common/moddiel.F90 1046 !Keep local copy of fft dimensions 1047 !Initialize computation of G^2 in cartesian coordinates 1048 if(.not.present(g2cart_in)) then 1049 if(present(g2cart_out)) then 1050 g2cart => g2cart_out 1051 else 1052 ABI_MALLOC(g2cart,(nfft)) 1053 end if 1054 id1=int(n1/2)+2 1055 id2=int(n2/2)+2 1056 id3=int(n3/2)+2 1057 count=0 1058 do i3=1,n3 1059 ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft) 1060 ig3=i3-int(i3/id3)*n3-1 1061 do i2=1,n2 1062 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 1063 ig2=i2-int(i2/id2)*n2-1 1064 1065 ii1=1 1066 do i1=ii1,n1 1067 ig1=i1-int(i1/id1)*n1-1 1068 ifft=ifft+1 1069 1070 b11=gprimd(1,1)*real(ig1,dp) 1071 b21=gprimd(2,1)*real(ig1,dp) 1072 b31=gprimd(3,1)*real(ig1,dp) 1073 b12=gprimd(1,2)*real(ig2,dp) 1074 b22=gprimd(2,2)*real(ig2,dp) 1075 b32=gprimd(3,2)*real(ig2,dp) 1076 b13=gprimd(1,3)*real(ig3,dp) 1077 b23=gprimd(2,3)*real(ig3,dp) 1078 b33=gprimd(3,3)*real(ig3,dp) 1079 1080 g2cart(ifft)=( & 1081 & (b11+b12+b13)**2& 1082 & +(b21+b22+b23)**2& 1083 & +(b31+b32+b33)**2& 1084 & ) 1085 do ifunc=1,nfunc 1086 ! compute the laplacian in Fourier space that is * (i x 2pi x G)**2 1087 laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi 1088 laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi 1089 end do 1090 end do 1091 end if 1092 end do 1093 end do 1094 if(.not.present(g2cart_out)) then 1095 ABI_FREE(g2cart) 1096 end if 1097 else 1098 g2cart => g2cart_in 1099 do ifunc=1,nfunc 1100 do ifft=1,nfft 1101 ! compute the laplacian in Fourier space that is * (i x 2pi x G)**2 1102 laplacerdfuncg(1,ifft,ifunc) = -rdfuncg(1,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi 1103 laplacerdfuncg(2,ifft,ifunc) = -rdfuncg(2,ifft,ifunc)*g2cart(ifft)*two_pi*two_pi 1104 end do 1105 end do 1106 end if 1107 1108 !get the result back into real space 1109 if(present(laplacerdfuncr)) then 1110 do ifunc=1,nfunc 1111 call fourdp(1,laplacerdfuncg(:,:,ifunc),laplacerdfuncr(:,ifunc),1,mpi_enreg,nfft,1,ngfft,0) 1112 end do 1113 end if 1114 1115 !deallocate pointers 1116 if((.not.present(rdfuncg_in)).and.(.not.present(rdfuncg_in))) then 1117 ABI_FREE(rdfuncg) 1118 end if 1119 if(.not.present(laplacerdfuncg_out)) then 1120 ABI_FREE(laplacerdfuncg) 1121 end if 1122 1123 end subroutine laplacian
m_spacepar/make_vectornd [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
make_vectornd
FUNCTION
For nuclear dipole moments m, compute vector potential A(r) = (m x (r-R))/|r-R|^3 in r space. This is done by computing A(G) followed by FFT.
NOTES
This code is copied and modified from m_spacepar/hartre where a very similar loop over G is done followed by FFT to real space
INPUTS
OUTPUT
vectornd(3,nfft)=Vector potential in real space, along Cartesian directions
SOURCE
83 subroutine make_vectornd(cplex,gsqcut,izero,mpi_enreg,natom,nfft,ngfft,nspden,nucdipmom,& 84 & rprimd,vectornd,xred) 85 86 !Arguments ------------------------------------ 87 !scalars 88 integer,intent(in) :: cplex,izero,natom,nfft,nspden 89 real(dp),intent(in) :: gsqcut 90 type(MPI_type),intent(in) :: mpi_enreg 91 !arrays 92 integer,intent(in) :: ngfft(18) 93 real(dp),intent(in) :: nucdipmom(3,natom),rprimd(3,3),xred(3,natom) 94 real(dp),intent(out) :: vectornd(nfft,nspden,3) 95 96 !Local variables------------------------------- 97 !scalars 98 integer,parameter :: im=2,re=1 99 integer :: i1,i2,i2_local,i23,i3,iatom,id1,id2,id3,ig,ig1,ig2,ig3,ig1max,ig2max,ig3max 100 integer :: ig1min,ig2min,ig3min 101 integer :: ii,ii1,ing,me_fft,n1,n2,n3,nd_atom,nd_atom_tot,nproc_fft 102 real(dp),parameter :: tolfix=1.000000001e0_dp 103 real(dp) :: cutoff,gqgm12,gqg2p3,gqgm23,gqgm13,gs2,gs3,gs,phase,ucvol 104 complex(dpc) :: prefac,cgr 105 !arrays 106 integer :: id(3) 107 integer,allocatable :: nd_list(:) 108 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 109 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 110 real(dp) :: gmet(3,3),gprimd(3,3),gqred(3),mcgc(3),rmet(3,3) 111 real(dp) :: rgbasis(3,3,3) 112 real(dp),allocatable :: gq(:,:),nd_m(:,:),ndvecr(:),work1(:,:),work2(:,:),work3(:,:) 113 114 115 ! ************************************************************************* 116 117 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 118 119 ! make list of atoms with nonzero nuclear dipole moments 120 ! in typical applications only 0 or 1 atoms have nonzero dipoles. This 121 ! code shouldn't even be called if all dipoles are zero. 122 nd_atom_tot = 0 123 do iatom = 1, natom 124 if (any(abs(nucdipmom(:,iatom))>tol8)) then 125 nd_atom_tot = nd_atom_tot + 1 126 end if 127 end do 128 129 ! construct the basis vectors of the generalized cross product 130 ! real space a, b, c (contained in rprimd) 131 ! reciprocal space a*, b*, c* (contained in gprimd) 132 ! for m x G will need a x a*, a x b* etc (9 a^b type basis vectors) 133 call wedge_basis(gprimd,rprimd,rgbasis) 134 135 ! note that nucdipmom is input as vectors in atomic units referenced 136 ! to cartesian coordinates 137 ABI_MALLOC(nd_list,(nd_atom_tot)) 138 ABI_MALLOC(nd_m,(3,nd_atom_tot)) 139 nd_atom_tot = 0 140 do iatom = 1, natom 141 if (any(abs(nucdipmom(:,iatom))>tol8)) then 142 nd_atom_tot = nd_atom_tot + 1 143 nd_list(nd_atom_tot) = iatom 144 ! the following expresses the dipole moment components in units of rprimd translations 145 nd_m(:,nd_atom_tot) = MATMUL(TRANSPOSE(gprimd),nucdipmom(:,iatom)) 146 end if 147 end do 148 149 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 150 nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft 151 152 prefac = -four_pi*j_dpc/(ucvol*two_pi) 153 154 ! Get the distrib associated with this fft_grid 155 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 156 157 ! Initialize a few quantities 158 cutoff=gsqcut*tolfix 159 160 ! In order to speed the routine, precompute the components of g+q 161 ! Also check if the booked space was large enough... 162 ABI_MALLOC(gq,(3,max(n1,n2,n3))) 163 do ii=1,3 164 id(ii)=ngfft(ii)/2+2 165 do ing=1,ngfft(ii) 166 ig=ing-(ing/id(ii))*ngfft(ii)-1 167 gq(ii,ing)=ig 168 end do 169 end do 170 ig1max=-1;ig2max=-1;ig3max=-1 171 ig1min=n1;ig2min=n2;ig3min=n3 172 173 ABI_MALLOC(work1,(2,nfft)) 174 ABI_MALLOC(work2,(2,nfft)) 175 ABI_MALLOC(work3,(2,nfft)) 176 work1=zero; work2=zero; work3=zero 177 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 178 179 ! Triple loop on each dimension 180 do i3=1,n3 181 ig3=i3-(i3/id3)*n3-1 182 ! Precompute some products that do not depend on i2 and i1 183 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 184 gqgm23=gq(3,i3)*gmet(2,3)*2 185 gqgm13=gq(3,i3)*gmet(1,3)*2 186 187 do i2=1,n2 188 ig2=i2-(i2/id2)*n2-1 189 if (fftn2_distrib(i2) == me_fft) then 190 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 191 gqgm12=gq(2,i2)*gmet(1,2)*2 192 gqg2p3=gqgm13+gqgm12 193 194 i2_local = ffti2_local(i2) 195 i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1)) 196 ! Do the test that eliminates the Gamma point outside of the inner loop 197 ii1=1 198 !if(i23==0 .and. ig2==0 .and. ig3==0)then 199 ! ii1=2 200 ! work1(re,1+i23)=zero 201 ! work1(im,1+i23)=zero 202 !end if 203 204 ! Final inner loop on the first dimension (note the lower limit) 205 do i1=ii1,n1 206 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 207 ig1 = i1 - (i1/id1)*n1 -1 208 ii=i1+i23 209 210 gqred(1) = gq(1,i1); gqred(2) = gq(2,i2); gqred(3) = gq(3,i3) 211 212 if( (gs .LE. cutoff) .AND. (gs .gt. tol8) )then 213 214 do iatom = 1, nd_atom_tot 215 nd_atom = nd_list(iatom) 216 phase = -two_pi*DOT_PRODUCT(xred(:,nd_atom),gqred(:)) 217 cgr = cmplx(cos(phase),sin(phase)) 218 219 ! cross product m x G 220 call wedge_product(mcgc,nd_m(:,iatom),gqred,rgbasis) 221 222 ! express mcgc relative to rprimd translations. This is done because 223 ! we wish ultimately to apply A.p to |cwavef>; in getghc_nucdip, the 224 ! p|cwavef> is done in reduced coordinates so do that here too, because 225 ! r.G has no need of the metric if both terms are in reduced coords 226 mcgc = MATMUL(TRANSPOSE(gprimd),mcgc) 227 228 work1(re,ii) = work1(re,ii) + real(prefac*cgr*mcgc(1)/gs) 229 work2(re,ii) = work2(re,ii) + real(prefac*cgr*mcgc(2)/gs) 230 work3(re,ii) = work3(re,ii) + real(prefac*cgr*mcgc(3)/gs) 231 232 work1(im,ii) = work1(im,ii) + aimag(prefac*cgr*mcgc(1)/gs) 233 work2(im,ii) = work2(im,ii) + aimag(prefac*cgr*mcgc(2)/gs) 234 work3(im,ii) = work3(im,ii) + aimag(prefac*cgr*mcgc(3)/gs) 235 236 end do 237 else 238 ! gs>cutoff 239 work1(re,ii)=zero 240 work1(im,ii)=zero 241 work2(re,ii)=zero 242 work2(im,ii)=zero 243 work3(re,ii)=zero 244 work3(im,ii)=zero 245 end if 246 247 end do ! End loop on i1 248 end if 249 end do ! End loop on i2 250 end do ! End loop on i3 251 252 ABI_FREE(gq) 253 ABI_FREE(nd_list) 254 ABI_FREE(nd_m) 255 256 if ( izero .EQ. 1 ) then 257 ! Set contribution of unbalanced components to zero 258 259 call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 260 call zerosym(work2,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 261 call zerosym(work3,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 262 263 end if 264 265 !note nspden effect--the nuclear vector potential conains no electron spin flip operator, 266 ! so vectornd(:,2,:) = vectornd(:,1,:) and vectornd(:,3:4,:) = zero 267 vectornd = zero 268 ! Fourier Transform 269 ABI_MALLOC(ndvecr,(cplex*nfft)) 270 ndvecr=zero 271 call fourdp(cplex,work1,ndvecr,1,mpi_enreg,nfft,1,ngfft,0) 272 vectornd(:,1,1)=ndvecr(:) 273 if (nspden .GE. 2) vectornd(:,2,1) = ndvecr(:) 274 ABI_FREE(work1) 275 276 ndvecr=zero 277 call fourdp(cplex,work2,ndvecr,1,mpi_enreg,nfft,1,ngfft,0) 278 vectornd(:,1,2) = ndvecr(:) 279 if (nspden .GE. 2) vectornd(:,2,2) = ndvecr(:) 280 ABI_FREE(work2) 281 282 ndvecr=zero 283 call fourdp(cplex,work3,ndvecr,1,mpi_enreg,nfft,1,ngfft,0) 284 vectornd(:,1,3) = ndvecr(:) 285 if (nspden .GE. 2) vectornd(:,2,3) = ndvecr(:) 286 ABI_FREE(work3) 287 ABI_FREE(ndvecr) 288 289 end subroutine make_vectornd
m_spacepar/meanvalue_g [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
meanvalue_g
FUNCTION
Compute the mean value of one wavefunction, in reciprocal space, for an operator that is real, diagonal in G-space: <wf|op|wf> For the time being, only spin-independent operators are treated.
INPUTS
diag(npw)=diagonal operator (real, spin-independent!) filter= if 1, need to filter on the value of diag, that must be less than huge(0.0d0)*1.d-11 otherwise, should be 0 istwf_k=storage mode of the vectors npw=number of planewaves of the vector nspinor=number of spinor components vect(2,npw*nspinor)=vector vect1(2,npw*nspinor*use_ndo)=vector1 (=vector in most of the cases) use_ndo = says if vect=/vect1
OUTPUT
ar=mean value
SOURCE
770 subroutine meanvalue_g(ar,diag,filter,istwf_k,mpi_enreg,npw,nspinor,vect,vect1,use_ndo,ar_im) 771 772 !Arguments ------------------------------------ 773 !scalars 774 integer,intent(in) :: filter,istwf_k,npw,nspinor,use_ndo 775 real(dp),intent(out) :: ar 776 real(dp),intent(out),optional :: ar_im 777 type(MPI_type),intent(in) :: mpi_enreg 778 !arrays 779 real(dp),intent(in) :: diag(npw),vect(2,npw*nspinor) 780 real(dp),intent(in) :: vect1(2,npw*nspinor) 781 782 !Local variables------------------------------- 783 !scalars 784 integer :: i1,ierr,ipw,jpw,me_g0 785 character(len=500) :: message 786 !arrays 787 788 ! ************************************************************************* 789 me_g0 = mpi_enreg%me_g0 790 791 DBG_CHECK(ANY(filter==(/0,1/)),"Wrong filter") 792 DBG_CHECK(ANY(nspinor==(/1,2/)),"Wrong nspinor") 793 DBG_CHECK(ANY(istwf_k==(/(ipw,ipw=1,9)/)),"Wrong istwf_k") 794 795 if(nspinor==2 .and. istwf_k/=1)then 796 write(message,'(a,a,a,i6,a,i6)')& 797 & 'When istwf_k/=1, nspinor must be 1,',ch10,& 798 & 'however, nspinor=',nspinor,', and istwf_k=',istwf_k 799 ABI_BUG(message) 800 end if 801 802 if(use_ndo==1 .and. (istwf_k==2 .and.me_g0==1)) then 803 ABI_BUG('use_ndo==1, not tested, use istwfk=1') 804 end if 805 806 ar=zero 807 if(present(ar_im)) ar_im=zero 808 809 !Normal storage mode 810 if(istwf_k==1)then 811 812 ! No filter 813 if(filter==0)then 814 !$OMP PARALLEL DO REDUCTION(+:ar) 815 do ipw=1,npw 816 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 817 end do 818 if(nspinor==2)then 819 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw) 820 do ipw=1+npw,2*npw 821 jpw=ipw-npw 822 ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 823 end do 824 end if 825 if(use_ndo==1)then 826 !$OMP PARALLEL DO REDUCTION(+:ar_im) 827 do ipw=1,npw 828 ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 829 end do 830 if(nspinor == 2) then 831 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw) 832 do ipw=1+npw,2*npw 833 jpw=ipw-npw 834 ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 835 end do 836 end if 837 end if 838 839 ! !$OMP PARALLEL DO REDUCTION(+:ar,ar_im) 840 ! do ipw=1,npw 841 ! ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 842 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 843 ! end do 844 ! if(nspinor==2)then 845 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 846 ! do ipw=1+npw,2*npw 847 ! ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 848 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw-npw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 849 ! end do 850 ! end if 851 else ! will filter 852 853 !$OMP PARALLEL DO REDUCTION(+:ar) 854 do ipw=1,npw 855 if(diag(ipw)<huge(0.0d0)*1.d-11)then 856 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 857 end if 858 end do 859 if(nspinor==2)then 860 !$OMP PARALLEL DO REDUCTION(+:ar) PRIVATE(jpw) 861 do ipw=1+npw,2*npw 862 jpw=ipw-npw 863 if(diag(jpw)<huge(0.0d0)*1.d-11)then 864 ar=ar+diag(jpw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 865 end if 866 end do 867 end if 868 if(use_ndo==1)then 869 if(.not.present(ar_im)) then 870 ABI_BUG("use_ndo true and ar_im not present") 871 end if 872 !$OMP PARALLEL DO REDUCTION(+:ar_im) 873 do ipw=1,npw 874 if(diag(ipw)<huge(0.0d0)*1.d-11)then 875 ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 876 end if 877 end do 878 if(nspinor == 2) then 879 !$OMP PARALLEL DO REDUCTION(+:ar_im) PRIVATE(jpw) 880 do ipw=1+npw,2*npw 881 jpw=ipw-npw 882 if(diag(jpw)<huge(0.0d0)*1.d-11)then 883 ar_im=ar_im+diag(jpw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 884 end if 885 end do 886 end if 887 end if 888 889 890 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 891 ! do ipw=1,npw 892 ! if(diag(ipw)<huge(0.0d0)*1.d-11)then 893 ! ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 894 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 895 ! end if 896 ! end do 897 ! if(nspinor==2)then 898 ! !$OMP PARALLEL DO PRIVATE(ipw) REDUCTION(+:ar,ar_im) 899 ! do ipw=1+npw,2*npw 900 ! if(diag(ipw-npw)<huge(0.0d0)*1.d-11)then 901 ! ar=ar+diag(ipw-npw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 902 ! if(use_ndo==1.and.nspinor==2) ar_im=ar_im+diag(ipw-npw)*(vect1(1,ipw)*vect(2,ipw)-vect1(2,ipw)*vect(1,ipw)) 903 ! end if 904 ! end do 905 ! end if ! nspinor==2 906 907 end if ! filter==0 908 909 else if(istwf_k>=2)then 910 911 if(filter==0)then 912 i1=1 913 if(istwf_k==2 .and. me_g0==1)then ! MPIWF need to know which proc has G=0 914 ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2 915 end if 916 917 !$OMP PARALLEL DO REDUCTION(+:ar) 918 do ipw=i1,npw 919 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 920 end do 921 922 else ! filter/=0 923 i1=1 924 if(istwf_k==2 .and. me_g0==1)then 925 if(diag(1)<huge(0.0d0)*1.d-11)then 926 ar=half*diag(1)*vect(1,1)*vect1(1,1) ; i1=2 927 end if 928 end if 929 930 !$OMP PARALLEL DO REDUCTION(+:ar) 931 do ipw=i1,npw 932 if(diag(ipw)<huge(0.0d0)*1.d-11)then 933 ar=ar+diag(ipw)*(vect(1,ipw)*vect1(1,ipw)+vect(2,ipw)*vect1(2,ipw)) 934 end if 935 end do 936 end if ! filter==0 937 938 ar=two*ar 939 940 end if ! istwf_k 941 942 !MPIWF need to make reduction on ar and ai . 943 if(mpi_enreg%paral_kgb==1)then 944 call xmpi_sum(ar,mpi_enreg%comm_bandspinorfft ,ierr) 945 if(present(ar_im))then 946 call xmpi_sum(ar_im,mpi_enreg%comm_bandspinorfft,ierr) 947 end if 948 end if 949 950 end subroutine meanvalue_g
m_spacepar/mkunitpawspherepot [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
mkunitpawspherepot
FUNCTION
Compute "potential" due to a sphere of radius r_paw at one of the ions, of strength 1. This is done for testing the completeness of the PAW projectors.
NOTES
INPUTS
OUTPUT
vunitpawspherepot(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX
SOURCE
309 subroutine mkunitpawspherepot(cplex,gsqcut,izero,mpi_enreg,natom,nfft,ngfft,& 310 & rpaw,rprimd,vunitpawspherepot,xred,& 311 the_atom) ! Optional arguments 312 313 !Arguments ------------------------------------ 314 !scalars 315 integer,intent(in) :: cplex,izero,natom,nfft 316 integer,intent(in),optional :: the_atom 317 real(dp),intent(in) :: gsqcut,rpaw 318 type(MPI_type),intent(in) :: mpi_enreg 319 !arrays 320 integer,intent(in) :: ngfft(18) 321 real(dp),intent(in) :: rprimd(3,3),xred(3,natom) 322 real(dp),intent(out) :: vunitpawspherepot(cplex*nfft) 323 324 !Local variables------------------------------- 325 !scalars 326 integer,parameter :: im=2,re=1 327 integer :: i1,i2,i2_local,i23,i3,iatom,id1,id2,id3,ig,ig1,ig2,ig3,ig1max,ig2max,ig3max 328 integer :: ig1min,ig2min,ig3min 329 integer :: ii,ii1,ing,me_fft,n1,n2,n3,nproc_fft,qeq0,qeq05 330 real(dp),parameter :: tolfix=1.000000001e0_dp 331 real(dp) :: cutoff,gqgm12,gqg2p3,gqgm23,gqgm13,gs2,gs3,gs,ogg0,ogr,ogrpre,phgr,phpaw,ucvol 332 character(len=500) :: message 333 !arrays 334 integer :: id(3) 335 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 336 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 337 real(dp) :: gmet(3,3),gprimd(3,3),gqred(3),qpt_(3),rmet(3,3) 338 real(dp),allocatable :: gq(:,:),work1(:,:) 339 340 341 ! ************************************************************************* 342 343 ! Check that cplex has an allowed value 344 if(cplex/=1 .and. cplex/=2)then 345 write(message, '(a,i0,a,a)' )& 346 'From the calling routine, cplex=',cplex,ch10,& 347 'but the only value allowed are 1 and 2.' 348 ABI_BUG(message) 349 end if 350 351 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 352 353 ogrpre = four_pi/(ucvol*(two_pi**3)) 354 ogg0 = four_pi*(rpaw**3)/(three*ucvol) 355 356 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 357 nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft 358 359 ! Get the distrib associated with this fft_grid 360 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 361 362 ! Initialize a few quantities 363 cutoff=gsqcut*tolfix 364 !carrying code over from hartre with minimal changes, in present case qpt always zero 365 qpt_=zero 366 qeq0=0 367 if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1 368 qeq05=0 369 if (qeq0==0) then 370 if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or.abs(abs(qpt_(3))-half)<tol12) qeq05=1 371 end if 372 373 if (present(the_atom)) then 374 iatom = the_atom 375 else 376 iatom = 1 377 end if 378 379 ! If cplex=1 then qpt_ should be 0 0 0 380 if (cplex==1.and. qeq0/=1) then 381 write(message,'(a,3e12.4,a,a)')& 382 'cplex=1 but qpt=',qpt_,ch10,& 383 'qpt should be 0 0 0.' 384 ABI_BUG(message) 385 end if 386 387 ! If FFT parallelism then qpt should not be 1/2 388 if (nproc_fft>1.and.qeq05==1) then 389 write(message, '(a,3e12.4,a,a)' )& 390 'FFT parallelism selected but qpt',qpt_,ch10,& 391 'qpt(i) should not be 1/2...' 392 ABI_ERROR(message) 393 end if 394 395 ! In order to speed the routine, precompute the components of g+q 396 ! Also check if the booked space was large enough... 397 ABI_MALLOC(gq,(3,max(n1,n2,n3))) 398 do ii=1,3 399 id(ii)=ngfft(ii)/2+2 400 do ing=1,ngfft(ii) 401 ig=ing-(ing/id(ii))*ngfft(ii)-1 402 gq(ii,ing)=ig+qpt_(ii) 403 end do 404 end do 405 ig1max=-1;ig2max=-1;ig3max=-1 406 ig1min=n1;ig2min=n2;ig3min=n3 407 408 ABI_MALLOC(work1,(2,nfft)) 409 work1=zero 410 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 411 412 ! Triple loop on each dimension 413 do i3=1,n3 414 ig3=i3-(i3/id3)*n3-1 415 ! Precompute some products that do not depend on i2 and i1 416 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 417 gqgm23=gq(3,i3)*gmet(2,3)*2 418 gqgm13=gq(3,i3)*gmet(1,3)*2 419 420 do i2=1,n2 421 ig2=i2-(i2/id2)*n2-1 422 if (fftn2_distrib(i2) == me_fft) then 423 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 424 gqgm12=gq(2,i2)*gmet(1,2)*2 425 gqg2p3=gqgm13+gqgm12 426 427 i2_local = ffti2_local(i2) 428 i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1)) 429 ! Do the test that eliminates the Gamma point outside of the inner loop 430 ii1=1 431 if(i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0)then 432 ii1=2 433 work1(re,1+i23)=ogg0 434 work1(im,1+i23)=zero 435 end if 436 437 ! Final inner loop on the first dimension (note the lower limit) 438 do i1=ii1,n1 439 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 440 ii=i1+i23 441 442 gqred(1) = gq(1,i1); gqred(2) = gq(2,i2); gqred(3) = gq(3,i3) 443 444 if(gs<=cutoff)then 445 446 447 ! Identify min/max indexes (to cancel unbalanced contributions later) 448 ! Count (q+g)-vectors with similar norm 449 if ((qeq05==1).and.(izero==1)) then 450 ig1=i1-(i1/id1)*n1-1 451 ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1) 452 ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2) 453 ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3) 454 end if 455 456 phpaw=two_pi*rpaw*gs 457 ogr=ogrpre*(sin(phpaw)-phpaw*cos(phpaw))/(gs**3) 458 459 phgr = -two_pi*DOT_PRODUCT(xred(:,iatom),gqred(:)) 460 461 work1(re,ii)=cos(phgr)*ogr 462 work1(im,ii)=sin(phgr)*ogr 463 else 464 ! gs>cutoff 465 work1(re,ii)=zero 466 work1(im,ii)=zero 467 end if 468 469 end do ! End loop on i1 470 end if 471 end do ! End loop on i2 472 end do ! End loop on i3 473 474 ABI_FREE(gq) 475 476 if (izero==1) then 477 ! Set contribution of unbalanced components to zero 478 479 if (qeq0==1) then !q=0 480 call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 481 482 else if (qeq05==1) then 483 !q=1/2; this doesn't work in parallel 484 ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2 485 ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2 486 ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2 487 if (abs(abs(qpt_(1))-half)<tol12) then 488 if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max) 489 if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min) 490 end if 491 if (abs(abs(qpt_(2))-half)<tol12) then 492 if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max) 493 if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min) 494 end if 495 if (abs(abs(qpt_(3))-half)<tol12) then 496 if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max) 497 if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min) 498 end if 499 call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,& 500 comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 501 end if 502 end if 503 504 ! Fourier Transform 505 call fourdp(cplex,work1,vunitpawspherepot,1,mpi_enreg,nfft,1,ngfft,0) 506 507 ABI_FREE(work1) 508 509 end subroutine mkunitpawspherepot
m_spacepar/redgr [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
redgr
FUNCTION
Compute reduced gradients of a real function on the usual unshifted fft grid. The gradient directions are the along the primitive reciprocal lattice vectors. The input function is intended to be a single spin component of the valence charge density, the valence + core charge densities or the first-order core charge density for use in frozen wf elastic tensor calculations within the GGA.
NOTES
Closely linked to xcden, but limited to Q=0, real charge densities, and unshifted grids.
INPUTS
mpi_enreg=information about MPI parallelization 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 frin(nfft)=real space input function
OUTPUT
frredgr(nfft,3)= reduced gradient of input function (same units as frin)
SOURCE
1154 subroutine redgr(frin,frredgr,mpi_enreg,nfft,ngfft) 1155 1156 !Arguments ------------------------------------ 1157 !scalars 1158 integer,intent(in) :: nfft 1159 type(MPI_type),intent(in) :: mpi_enreg 1160 !arrays 1161 integer,intent(in) :: ngfft(18) 1162 real(dp),intent(in) :: frin(nfft) 1163 real(dp),intent(out) :: frredgr(nfft,3) 1164 1165 !Local variables------------------------------- 1166 !scalars 1167 integer :: cplex_tmp,i1,i2,i3,id,idir,ifft,ig,ii,ing,n1,n2,n3 1168 !arrays 1169 real(dp),allocatable :: gg(:,:),wkcmpx(:,:),work(:),workgr(:,:) 1170 1171 ! ************************************************************************* 1172 1173 !Only real arrays are treated 1174 cplex_tmp=1 1175 1176 !Keep local copy of fft dimensions 1177 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1178 1179 !In order to speed the routine, precompute the components of g, including 2pi factor 1180 ABI_MALLOC(gg,(max(n1,n2,n3),3)) 1181 do ii=1,3 1182 id=ngfft(ii)/2+2 1183 do ing=1,ngfft(ii) 1184 ig=ing-(ing/id)*ngfft(ii)-1 1185 gg(ing,ii)=two_pi*ig 1186 end do 1187 ! Note that the G <-> -G symmetry must be maintained 1188 if(mod(ngfft(ii),2)==0)gg(ngfft(ii)/2+1,ii)=zero 1189 end do 1190 1191 ABI_MALLOC(wkcmpx,(2,nfft)) 1192 ABI_MALLOC(work,(nfft)) 1193 ABI_MALLOC(workgr,(2,nfft)) 1194 1195 !Obtain rho(G) in wkcmpx from input rho(r) 1196 work(:)=frin(:) 1197 1198 call fourdp(cplex_tmp,wkcmpx,work,-1,mpi_enreg,nfft,1,ngfft,0) 1199 1200 !Gradient calculation for three reduced components in turn. 1201 !Code duplicated to remove logic from loops. 1202 do idir=1,3 1203 if(idir==1) then 1204 !$OMP PARALLEL DO PRIVATE(ifft) 1205 do i3=1,n3 1206 ifft=(i3-1)*n1*n2 1207 do i2=1,n2 1208 do i1=1,n1 1209 ifft=ifft+1 1210 ! Multiply by i 2pi G(idir) 1211 workgr(2,ifft)= gg(i1,idir)*wkcmpx(1,ifft) 1212 workgr(1,ifft)=-gg(i1,idir)*wkcmpx(2,ifft) 1213 end do 1214 end do 1215 end do 1216 else if(idir==2) then 1217 !$OMP PARALLEL DO PRIVATE(ifft) 1218 do i3=1,n3 1219 ifft=(i3-1)*n1*n2 1220 do i2=1,n2 1221 do i1=1,n1 1222 ifft=ifft+1 1223 ! Multiply by i 2pi G(idir) 1224 workgr(2,ifft)= gg(i2,idir)*wkcmpx(1,ifft) 1225 workgr(1,ifft)=-gg(i2,idir)*wkcmpx(2,ifft) 1226 end do 1227 end do 1228 end do 1229 else 1230 !$OMP PARALLEL DO PRIVATE(ifft) 1231 do i3=1,n3 1232 ifft=(i3-1)*n1*n2 1233 do i2=1,n2 1234 do i1=1,n1 1235 ifft=ifft+1 1236 ! Multiply by i 2pi G(idir) 1237 workgr(2,ifft)= gg(i3,idir)*wkcmpx(1,ifft) 1238 workgr(1,ifft)=-gg(i3,idir)*wkcmpx(2,ifft) 1239 end do 1240 end do 1241 end do 1242 end if !idir 1243 1244 call fourdp(cplex_tmp,workgr,work,1,mpi_enreg,nfft,1,ngfft,0) 1245 1246 !$OMP PARALLEL DO 1247 do ifft=1,nfft 1248 frredgr(ifft,idir)=work(ifft) 1249 end do 1250 1251 end do !idir 1252 1253 ABI_FREE(gg) 1254 ABI_FREE(wkcmpx) 1255 ABI_FREE(work) 1256 ABI_FREE(workgr) 1257 1258 end subroutine redgr
m_spacepar/rotate_rho [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
rotate_rho
FUNCTION
rotate density in real and reciprocal space
INPUTS
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90) ngfft=array of dimensions for different FFT grids nspden=number of spin-density components rhor1(cplex*nfft,nspden)=array for Fourier transform of RF electron density === if psps%usepaw==1 TODO: extend to PAW pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data symrel1=single symmetry operation in real space to apply to rho tnon = eventual translation associated to symrel1
OUTPUT
rhog1_eq(2,nfft)= symmetric density in reciprocal space for equivalent perturbation rhor1_eq(cplex*nfft,nspden) = symmetric density in real space for equivalent perturbation
SOURCE
2378 subroutine rotate_rho(cplex, itirev, mpi_enreg, nfft, ngfft, nspden, & 2379 & rhor1, rhog1_eq, rhor1_eq, symrel1, tnon) 2380 2381 !args 2382 integer,intent(in) :: cplex, nfft, nspden, itirev 2383 integer,intent(in) :: ngfft(18) 2384 2385 integer, intent(in) :: symrel1(3,3) 2386 real(dp),intent(in) :: tnon(3) 2387 real(dp),intent(inout) :: rhor1(cplex*nfft,nspden) 2388 2389 real(dp),intent(out) :: rhog1_eq(2,nfft) 2390 real(dp),intent(out) :: rhor1_eq(cplex*nfft,nspden) 2391 2392 type(MPI_type),intent(in) :: mpi_enreg 2393 2394 ! local vars 2395 integer :: id1, id2, id3 2396 integer :: n1, n2, n3, nd2 2397 integer :: l1, l2, l3 2398 integer :: i1, i2, i3, ind1, ind2 2399 integer :: j1, j2, j3 2400 integer :: k1, k2, k3 2401 integer :: nproc_fft, ispden, me_fft 2402 real(dp) :: arg 2403 logical :: t_tnon_nonzero 2404 2405 real(dp) :: phnon1(2) 2406 real(dp), allocatable :: workg(:,:), workg_eq(:,:) 2407 character(len=500) :: message 2408 2409 ! ************************************************************************* 2410 2411 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft 2412 2413 id1=n1/2+2 2414 id2=n2/2+2 2415 id3=n3/2+2 2416 2417 rhog1_eq = zero 2418 rhor1_eq = zero 2419 2420 if (itirev == 2) then 2421 write (message,'(3a,9I4,1a)') 'using time reversal. ',ch10,'Symrel1 = ', symrel1, ch10 2422 else 2423 write (message,'(3a,9I4,1a)') 'no time reversal. ',ch10,'Symrel1 = ', symrel1, ch10 2424 end if 2425 !call wrtout(std_out, message, 'COLL') 2426 2427 t_tnon_nonzero = (any(abs(tnon) > tol12)) 2428 2429 ! eventually, for FFT parallelization 2430 ! call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 2431 2432 ABI_MALLOC(workg,(2,nfft)) 2433 ABI_MALLOC(workg_eq,(2,nfft)) 2434 do ispden = 1, nspden 2435 2436 ! fft input rhor1 to reciprocal space: uses work* as a buffer 2437 call fourdp(cplex,workg,rhor1(:,ispden),-1,mpi_enreg,nfft,1,ngfft,0) 2438 2439 ! below taken from irrzg and setsym 2440 ! Loop over reciprocal space grid points: 2441 ! loop over local points in workg, and get back transform from rhog1, 2442 ! which is presumed complete on each proc 2443 ind1=0 2444 do i3=1,n3 2445 do i2=1,n2 2446 ! if(fftn2_distrib(i2)/=me_fft) cycle ! this ind is not to be treated by me_fft 2447 do i1=1,n1 2448 2449 ind1=ind1+1 2450 ! r2=ffti2_local(i2+1) - 1 2451 ! ind=n1*(nd2*i3+r2)+i1+1 !this is ind in the current proc 2452 2453 ! Get location of G vector (grid point) centered at 0 0 0 2454 l1=i1-(i1/id1)*n1-1 2455 l2=i2-(i2/id2)*n2-1 2456 l3=i3-(i3/id3)*n3-1 2457 2458 ! Get rotated G vector Gj for each symmetry element 2459 ! -- here we use the TRANSPOSE of symrel1; assuming symrel1 expresses 2460 ! the rotation in real space, the transpose is then appropriate 2461 ! for G space symmetrization (p. 1172d,e of notes, 2 June 1995). 2462 j1=symrel1(1,1)*l1+symrel1(2,1)*l2+symrel1(3,1)*l3 2463 j2=symrel1(1,2)*l1+symrel1(2,2)*l2+symrel1(3,2)*l3 2464 j3=symrel1(1,3)*l1+symrel1(2,3)*l2+symrel1(3,3)*l3 2465 2466 ! Map into [0,n-1] and then add 1 for array index in [1,n] 2467 k1=1+mod(n1+mod(j1,n1),n1) 2468 k2=1+mod(n2+mod(j2,n2),n2) 2469 k3=1+mod(n3+mod(j3,n3),n3) 2470 2471 ! Get linear index of rotated point Gj 2472 ind2=k1+n1*((k2-1)+n2*(k3-1)) 2473 ! r2=ffti2_local(j2+1) - 1 2474 ! ind=n1*(nd2*j3+r2)+j1+1 !this is ind may be in another proc!! 2475 2476 phnon1(1) = one 2477 phnon1(2) = zero 2478 if (t_tnon_nonzero) then 2479 ! compute exp(-2*Pi*I*G dot tau) using original G 2480 ! NB: this phase is same as that in irrzg and phnons1, and corresponds to complex conjugate of phase from G to Gj; 2481 ! we use it immediately below, to go _to_ workg(ind1) 2482 ! TODO : replace this with complex powers of exp(2pi tnon(1)) etc... 2483 arg=two_pi*(dble(l1)*tnon(1)+dble(l2)*tnon(2)+dble(l3)*tnon(3)) 2484 phnon1(1) = cos(arg) 2485 phnon1(2) =-sin(arg) 2486 end if 2487 2488 ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G) 2489 workg_eq (1, ind1) = phnon1(1) * workg(1, ind2) & 2490 & - phnon1(2) * workg(2, ind2) 2491 workg_eq (2, ind1) = phnon1(1) * workg(2, ind2) & 2492 & + phnon1(2) * workg(1, ind2) 2493 2494 end do 2495 end do 2496 end do 2497 2498 ! accumulate rhog1_eq 2499 if (ispden == 1) rhog1_eq = workg_eq 2500 2501 ! FFT back to real space to get rhor1_eq 2502 ! Pull out full or spin up density, now symmetrized 2503 call fourdp(cplex,workg_eq,rhor1_eq(:,ispden),1,mpi_enreg,nfft,1,ngfft,0) 2504 2505 end do !nspden 2506 2507 ABI_FREE(workg) 2508 ABI_FREE(workg_eq) 2509 2510 end subroutine rotate_rho
m_spacepar/setsym [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
setsym
FUNCTION
Set up irreducible zone in G space by direct calculation. Do not call this routine if nsym=1 (only identity symmetry). Only indsym and symrec get returned if iscf=0. symrec needed to symmetrize coordinate gradients in sygrad. (symrec is redundant and could be removed later in favor of symrel)
INPUTS
iscf=(<= 0 =>non-SCF), >0 => SCF natom=number of atoms in unit 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 nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetries in space group (at least 1) symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry operations in terms of real space primitive translations tnons(3,nsym)=nonsymmorphic translations of space group in terms of real space primitive translations (may be 0) typat(natom)=atom type (integer) for each atom xred(3,natom)=atomic coordinates in terms of real space primitive translations
OUTPUT
indsym(4,nsym,natom)=indirect indexing of atom labels--see subroutine symatm for definition (if nsym>1) irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases symrec(3,3,nsym)=symmetry operations in terms of reciprocal space primitive translations (if nsym>1)
NOTES
nsppol and nspden are needed in case of (anti)ferromagnetic symmetry operations
SOURCE
2552 subroutine setsym(indsym,irrzon,iscf,natom,nfft,ngfft,nspden,nsppol,nsym,phnons,& 2553 & symafm,symrec,symrel,tnons,typat,xred) 2554 2555 !Arguments ------------------------------------ 2556 !scalars 2557 integer,intent(in) :: iscf,natom,nfft,nspden,nsppol,nsym 2558 !arrays 2559 integer,intent(in) :: ngfft(18),symafm(nsym),symrel(3,3,nsym),typat(natom) 2560 integer,intent(out) :: indsym(4,nsym,natom) 2561 integer,intent(inout) :: irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4)) !vz_i 2562 integer,intent(out) :: symrec(3,3,nsym) 2563 real(dp),intent(in) :: tnons(3,nsym),xred(3,natom) 2564 real(dp),intent(out) :: phnons(2,nfft,(nspden/nsppol)-3*(nspden/4)) 2565 2566 !Local variables------------------------------- 2567 !scalars 2568 integer :: isym,ierr 2569 real(dp) :: tolsym8 2570 !arrays 2571 integer,allocatable :: determinant(:) 2572 real(dp) :: tsec(2) 2573 2574 ! ************************************************************************* 2575 2576 call timab(6,1,tsec) 2577 2578 !Check that symmetries have unity determinant 2579 ABI_MALLOC(determinant,(nsym)) 2580 call symdet(determinant,nsym,symrel) 2581 ABI_FREE(determinant) 2582 2583 2584 !Get the symmetry matrices in terms of reciprocal basis 2585 do isym=1,nsym 2586 call mati3inv(symrel(:,:,isym),symrec(:,:,isym)) 2587 end do 2588 2589 !Check for group closure 2590 call chkgrp(nsym,symafm,symrel,ierr) 2591 ABI_CHECK(ierr==0,"Error in group closure") 2592 2593 call chkgrp(nsym,symafm,symrec,ierr) 2594 ABI_CHECK(ierr==0,"Error in group closure") 2595 2596 !Obtain a list of rotated atom labels: 2597 tolsym8=tol8 2598 call symatm(indsym,natom,nsym,symrec,tnons,tolsym8,typat,xred,print_indsym=10) 2599 2600 !If non-SCF calculation, or nsym==1, do not need IBZ data 2601 if ( (iscf>0 .or. iscf==-3) .and. nsym>1 ) then 2602 ! Locate irreducible zone in reciprocal space for symmetrization: 2603 call irrzg(irrzon,nspden,nsppol,nsym,ngfft(1),ngfft(2),ngfft(3),phnons,symafm,symrel,tnons) 2604 end if 2605 2606 call timab(6,2,tsec) 2607 2608 end subroutine setsym
m_spacepar/symrhg [ Functions ]
[ Top ] [ m_spacepar ] [ Functions ]
NAME
symrhg
FUNCTION
From rho(r), generate rho(G), symmetrize it, and come back to the real space for a symmetrized rho(r).
INPUTS
cplex=1 if rhor is real, 2 if rhor is complex gprimd(3,3)=dimensional reciprocal space primitive translations irrzon(nfft,2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data mpi_enreg=information about MPI parallelization 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 nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetry elements. phnons(2,nfft,(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases rprimd(3,3)=dimensional real space primitive translations symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry matrices in real space (integers) tnons(3,nsym)=reduced nonsymmorphic translations
OUTPUT
rhog(2,nfft)=symmetrized rho(G) (total) electron density in G space
SIDE EFFECTS
Input/Output rhor(cplex*nfft,nspden)=array for electron density in electrons/bohr**3. Input, but also output, if symmetrization is applied. Also output if nspden > 1 (change spin components)
NOTES
When using spin-polarization (nspden==2), put total density in first half of rhor array and spin up in second half If (nspden=2 and nsppol=2) the density is transformed as (up,down) => (up+down,up) If (nspden=2 and nsppol=1) anti-ferromagnetic symmetry operations must be used, such as to transform (2*up) => (up+down,up) In spin-polarized, and if there is no symmetry to be applied on the system, only the total density is generated in G space
SOURCE
1487 subroutine symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfftot,ngfft,nspden,nsppol,nsym,& 1488 & phnons,rhog,rhor,rprimd,symafm,symrel,tnons) 1489 1490 !Arguments ------------------------------------ 1491 !scalars 1492 integer,intent(in) :: cplex,nfft,nfftot,nspden,nsppol,nsym 1493 type(MPI_type),intent(in) :: mpi_enreg 1494 !arrays 1495 integer,intent(in) :: irrzon(nfftot**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),ngfft(18) 1496 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 1497 real(dp),intent(in) :: gprimd(3,3),phnons(2,nfftot**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),rprimd(3,3) 1498 real(dp),intent(inout) :: rhor(cplex*nfft,nspden) 1499 real(dp),intent(out) :: rhog(2,nfft) 1500 real(dp),intent(in) :: tnons(3,nsym) 1501 1502 !Local variables------------------------------- 1503 !scalars 1504 integer :: id1,id2,id3,ier,imagn,ind,ind2,indsy,ispden,isym,iup,izone,izone_max,j,j1,j2,j3,jsym 1505 integer :: k1,k2,k3,l1,l2,l3,me_fft 1506 integer :: n1,n2,n3,nd2,nproc_fft,nspden_eff,nsym_used,numpt,nup 1507 integer :: r2,rep,spaceComm 1508 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism 1509 real(dp) :: arg,tau1,tau2,tau3 1510 real(dp) :: magxsu1,magxsu2,magysu1,magysu2,magzsu1,magzsu2,mxi,mxr,myi,myr,mzi,mzr,phi,phr,rhosu1,rhosu2 1511 !character(len=500) :: message 1512 !arrays 1513 integer,allocatable :: isymg(:) 1514 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 1515 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 1516 real(dp) :: tsec(2) 1517 real(dp),allocatable :: magngx(:,:),magngy(:,:),magngz(:,:) 1518 real(dp),allocatable :: rhosu1_arr(:),rhosu2_arr(:),work(:) 1519 real(dp),allocatable :: symafm_used(:),symrec_cart(:,:,:),symrel_cart(:,:,:),tnons_used(:,:) 1520 1521 !************************************************************************* 1522 ! 1523 !Note the timing channel 17 excludes the different Fourier transforms 1524 1525 ABI_MALLOC(work,(cplex*nfft)) 1526 1527 !Special treatment for spin-polarized case 1528 if(nspden==2 .and. nsppol==2) then 1529 ! When nspden=2 and nsppol=2, put total density in first half 1530 ! of rhor array and spin up in second half (up,down) => (up+down,up) 1531 call timab(17,1,tsec) 1532 work(:)=rhor(:,1) ! up => work 1533 rhor(:,1)=rhor(:,1)+rhor(:,2) ! up+down 1534 rhor(:,2)=work(:) ! work => up 1535 call timab(17,2,tsec) 1536 end if 1537 1538 !Special treatment for antiferromagnetism case 1539 if(nspden==2 .and. nsppol==1) then 1540 call timab(17,1,tsec) 1541 ! When nspden=2 and nsppol=1, (2*up) => (2*up,up) 1542 ! Indeed, what was delivered to the present routine is a "total" density, 1543 ! obtained from occupation numbers varying between 0 and 2, 1544 ! but for spin up only potential. 1545 rhor(:,2)=half*rhor(:,1) 1546 call timab(17,2,tsec) 1547 end if 1548 1549 !Special treatment for non-collinear magnetism case 1550 if(nspden==4) then 1551 call timab(17,1,tsec) 1552 !FR the half factors missing are recovered in dfpt_mkvxc_noncoll and dfpt_accrho 1553 rhor(:,1)=rhor(:,1)+rhor(:,4) !nup+ndown 1554 rhor(:,2)=rhor(:,2)-rhor(:,1) !mx (n+mx-n) 1555 rhor(:,3)=rhor(:,3)-rhor(:,1) !my (n+my-n) 1556 rhor(:,4)=rhor(:,1)-two*rhor(:,4) !mz=n-2ndown 1557 call timab(17,2,tsec) 1558 end if 1559 1560 1561 if(nsym==1)then 1562 1563 if(nspden==2 .and. nsppol==1) then ! There must be at least one anti-ferromagnetic operation 1564 ABI_BUG('In the antiferromagnetic case, nsym cannot be 1') 1565 end if 1566 1567 ! If not using symmetry, still want total density in G space rho(G). 1568 ! Fourier transform (incl normalization) to get rho(G) 1569 work(:)=rhor(:,1) 1570 call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,1,ngfft,0) 1571 else 1572 1573 ! Treat either full density, spin-up density or magnetization 1574 ! Note the decrease of ispden to the value 1, in order to finish 1575 ! with rhog of the total density (and not the spin-up density or magnetization) 1576 nspden_eff=nspden;if (nspden==4) nspden_eff=1 1577 do ispden=nspden_eff,1,-1 1578 1579 ! Prepare the density to be symmetrized, in the reciprocal space 1580 if(nspden==1 .or. nsppol==2 .or. (nspden==4.and.(.not.afm_noncoll)))then 1581 imagn=1 1582 nsym_used=0 1583 do isym=1,nsym 1584 if(symafm(isym)==1)nsym_used=nsym_used+1 1585 ! DEBUG 1586 ! write(std_out,*)' symrhg : isym,symafm(isym)',isym,symafm(isym) 1587 ! ENDDEBUG 1588 end do 1589 else if(nspden==2 .and. nsppol==1)then ! antiferromagnetic case 1590 imagn=ispden 1591 nsym_used=nsym/ispden 1592 else if (nspden==4) then 1593 imagn=1 1594 nsym_used=nsym/ispden 1595 end if 1596 1597 ! write(std_out,*)' symrhg : nsym_used=',nsym_used 1598 1599 ! rhor -fft-> rhog (rhog is used as work space) 1600 ! Note : it should be possible to reuse rhog in the antiferromagnetic case this would avoid one FFT 1601 work(:)=rhor(:,ispden) 1602 call fourdp(cplex,rhog,work,-1,mpi_enreg,nfft,1,ngfft,0) 1603 if (nspden==4) then 1604 ABI_MALLOC(magngx,(2,nfft)) 1605 ABI_MALLOC(magngy,(2,nfft)) 1606 ABI_MALLOC(magngz,(2,nfft)) 1607 work(:)=rhor(:,2) 1608 call fourdp(cplex,magngx,work,-1,mpi_enreg,nfft,1,ngfft,0) 1609 work(:)=rhor(:,3) 1610 call fourdp(cplex,magngy,work,-1,mpi_enreg,nfft,1,ngfft,0) 1611 work(:)=rhor(:,4) 1612 call fourdp(cplex,magngz,work,-1,mpi_enreg,nfft,1,ngfft,0) 1613 end if 1614 1615 ! Begins the timing here only , to exclude FFTs 1616 call timab(17,1,tsec) 1617 1618 n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nproc_fft=ngfft(10);me_fft=ngfft(11);nd2=n2/nproc_fft 1619 1620 ! Get the distrib associated with this fft_grid 1621 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 1622 1623 ! The following is only valid for total, up or dn density 1624 ! ------------------------------------------------------- 1625 1626 ! Get maxvalue of izone 1627 izone_max=count(irrzon(:,2,imagn)>0) 1628 ABI_MALLOC(rhosu1_arr,(izone_max)) 1629 ABI_MALLOC(rhosu2_arr,(izone_max)) 1630 1631 numpt=0 1632 do izone=1,nfftot 1633 1634 ! Get repetition number 1635 rep=irrzon(izone,2,imagn) 1636 if(rep==0)exit 1637 1638 ! Compute number of unique points in this symm class: 1639 nup=nsym_used/rep 1640 1641 ! Accumulate charge over equivalent points 1642 rhosu1=zero 1643 rhosu2=zero 1644 do iup=1,nup 1645 ind=irrzon(iup+numpt,1,imagn) 1646 j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2); 1647 if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft 1648 r2=ffti2_local(j2+1) - 1 1649 ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc 1650 rhosu1=rhosu1+rhog(1,ind)*phnons(1,iup+numpt,imagn)& 1651 & -rhog(2,ind)*phnons(2,iup+numpt,imagn) 1652 rhosu2=rhosu2+rhog(2,ind)*phnons(1,iup+numpt,imagn)& 1653 & +rhog(1,ind)*phnons(2,iup+numpt,imagn) 1654 end if 1655 1656 end do 1657 rhosu1=rhosu1/dble(nup) 1658 rhosu2=rhosu2/dble(nup) 1659 rhosu1_arr(izone)=rhosu1 1660 rhosu2_arr(izone)=rhosu2 1661 ! Keep index of how many points have been considered: 1662 numpt=numpt+nup 1663 1664 end do ! End loop over izone 1665 1666 ! Reduction in case of FFT parallelization 1667 if(mpi_enreg%nproc_fft>1)then 1668 spaceComm=mpi_enreg%comm_fft 1669 call xmpi_sum(rhosu1_arr,spaceComm,ier) 1670 call xmpi_sum(rhosu2_arr,spaceComm,ier) 1671 end if 1672 1673 ! Now symmetrize the density 1674 numpt=0 1675 do izone=1,nfftot 1676 1677 ! Get repetition number 1678 rep=irrzon(izone,2,imagn) 1679 if(rep==0)exit 1680 1681 ! Compute number of unique points in this symm class: 1682 nup=nsym_used/rep 1683 1684 ! Define symmetrized rho(G) at equivalent points: 1685 do iup=1,nup 1686 ind=irrzon(iup+numpt,1,imagn) 1687 ! decompose ind-1=n1(n2 j3+ j2)+j1 1688 j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2); 1689 if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft 1690 r2=ffti2_local(j2+1) - 1 1691 ! ind in the proc ind-1=n1(nd2 j3+ r2)+j1 1692 ind=n1*(nd2*j3+r2)+j1+1 !this is ind in the current proc 1693 rhog(1,ind)=rhosu1_arr(izone)*phnons(1,iup+numpt,imagn)& 1694 & +rhosu2_arr(izone)*phnons(2,iup+numpt,imagn) 1695 rhog(2,ind)=rhosu2_arr(izone)*phnons(1,iup+numpt,imagn)& 1696 & -rhosu1_arr(izone)*phnons(2,iup+numpt,imagn) 1697 end if 1698 end do 1699 1700 ! Keep index of how many points have been considered: 1701 numpt=numpt+nup 1702 1703 end do ! End loop over izone 1704 1705 ABI_FREE(rhosu1_arr) 1706 ABI_FREE(rhosu2_arr) 1707 1708 ! The following is only valid for magnetization 1709 ! --------------------------------------------- 1710 if (nspden==4) then 1711 1712 id1=n1/2+2 1713 id2=n2/2+2 1714 id3=n3/2+2 1715 1716 ! Transfer symmetries in cartesian coordinates 1717 ! Compute symmetries in reciprocal space in cartesian coordinates 1718 ABI_MALLOC(symrec_cart,(3,3,nsym_used)) 1719 ABI_MALLOC(symrel_cart,(3,3,nsym_used)) 1720 ABI_MALLOC(symafm_used,(nsym_used)) 1721 ABI_MALLOC(tnons_used,(3,nsym_used)) 1722 jsym=0 1723 do isym=1,nsym 1724 if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle 1725 jsym=jsym+1 1726 tnons_used(:,jsym)=tnons(:,isym) 1727 symafm_used(jsym)=dble(symafm(isym)) 1728 call symredcart(rprimd,gprimd,symrel_cart(:,:,jsym),symrel(:,:,isym)) 1729 call matr3inv(symrel_cart(:,:,jsym),symrec_cart(:,:,jsym)) 1730 end do 1731 1732 numpt=count(irrzon(:,1,imagn)>0) 1733 ABI_MALLOC(isymg,(numpt)) 1734 isymg=0 1735 ABI_MALLOC(rhosu1_arr,(3*izone_max)) 1736 ABI_MALLOC(rhosu2_arr,(3*izone_max)) 1737 1738 ! Accumulate magnetization over equivalent points 1739 ! Use all symmetries (not only those linking different g points) 1740 ! Use Inverse[Transpose[symrel]]=symrec 1741 numpt=0 1742 do izone=1,izone_max 1743 magxsu1=zero;magxsu2=zero 1744 magysu1=zero;magysu2=zero 1745 magzsu1=zero;magzsu2=zero 1746 ind=irrzon(1+numpt,1,1) 1747 rep=irrzon(izone,2,1) 1748 nup=nsym_used/rep 1749 ! Get coordinates in the range [0,n-1] 1750 j=ind-1;l1=modulo(j,n1);l2=modulo(j/n1,n2);l3=j/(n1*n2) 1751 ! Get location of G vector (grid point) centered at 0 0 0 1752 !TO BE UNCOMMENTED 1753 l3=l3-(l3/id3)*n3 1754 l2=l2-(l2/id2)*n2 1755 l1=l1-(l1/id1)*n1 1756 1757 jsym=0 1758 do isym=1,nsym 1759 if (symafm(isym)/=1.and.(.not.afm_noncoll)) cycle 1760 jsym=jsym+1 1761 ! The G vectors should transform as vectors in reciprocal space 1762 ! However, one acts with the INVERSE of the symmetry operation => Inverse[symrec]=Transpose[symrel] 1763 j1=symrel(1,1,isym)*l1+symrel(2,1,isym)*l2+symrel(3,1,isym)*l3 1764 j2=symrel(1,2,isym)*l1+symrel(2,2,isym)*l2+symrel(3,2,isym)*l3 1765 j3=symrel(1,3,isym)*l1+symrel(2,3,isym)*l2+symrel(3,3,isym)*l3 1766 k1=map_symrhg(j1,n1);k2=map_symrhg(j2,n2);k3=map_symrhg(j3,n3) 1767 indsy=1+k1+n1*(k2+n2*k3) 1768 ind2=-1;iup=numpt 1769 do while (ind2/=indsy.and.iup<numpt+nup) 1770 iup=iup+1;ind2=irrzon(iup,1,1) 1771 end do 1772 if (ind2/=indsy) then 1773 ABI_ERROR("ind2/=indsy in symrhg !") 1774 end if 1775 if (isymg(iup)==0) isymg(iup)=jsym 1776 if(fftn2_distrib(modulo((indsy-1)/n1,n2) + 1) == me_fft ) then ! this is indsy is to be treated by me_fft 1777 indsy=n1*(nd2*k3+ ffti2_local(k2+1) -1)+k1+1 ! this is indsy in the current proc 1778 1779 ! Working on this: the present coding will be detrimental for speed ! cos and sin are recomputed many times ! 1780 tau1=tnons_used(1,jsym) 1781 tau2=tnons_used(2,jsym) 1782 tau3=tnons_used(3,jsym) 1783 if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then 1784 ! Compute exp(-2*Pi*I*G dot tau) using original G (equivalent of phnons in the collinear case) 1785 arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3) 1786 phr=cos(arg) 1787 phi=-sin(arg) 1788 else 1789 phr=one 1790 phi=zero 1791 end if 1792 phr=phr*symafm_used(jsym) 1793 phi=phi*symafm_used(jsym) 1794 !TO BE COMMENTED 1795 ! phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons 1796 ! phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90) 1797 1798 ! The magnetization should transform as a vector in real space 1799 ! However, one acts with the INVERSE of the symmetry operation. 1800 ! => Inverse[symrel_cart] = Transpose[symrel_cart] because symrel_cart is unitary ?!?!? 1801 mxr=symrel_cart(1,1,jsym)*magngx(1,indsy)+symrel_cart(1,2,jsym)*magngy(1,indsy)+symrel_cart(1,3,jsym)*magngz(1,indsy) 1802 mxi=symrel_cart(1,1,jsym)*magngx(2,indsy)+symrel_cart(1,2,jsym)*magngy(2,indsy)+symrel_cart(1,3,jsym)*magngz(2,indsy) 1803 myr=symrel_cart(2,1,jsym)*magngx(1,indsy)+symrel_cart(2,2,jsym)*magngy(1,indsy)+symrel_cart(2,3,jsym)*magngz(1,indsy) 1804 myi=symrel_cart(2,1,jsym)*magngx(2,indsy)+symrel_cart(2,2,jsym)*magngy(2,indsy)+symrel_cart(2,3,jsym)*magngz(2,indsy) 1805 mzr=symrel_cart(3,1,jsym)*magngx(1,indsy)+symrel_cart(3,2,jsym)*magngy(1,indsy)+symrel_cart(3,3,jsym)*magngz(1,indsy) 1806 mzi=symrel_cart(3,1,jsym)*magngx(2,indsy)+symrel_cart(3,2,jsym)*magngy(2,indsy)+symrel_cart(3,3,jsym)*magngz(2,indsy) 1807 1808 ! mxr=symrel_cart(1,1,jsym)*magngx(1,indsy)+symrel_cart(2,1,jsym)*magngy(1,indsy)+symrel_cart(3,1,jsym)*magngz(1,indsy) 1809 ! mxi=symrel_cart(1,1,jsym)*magngx(2,indsy)+symrel_cart(2,1,jsym)*magngy(2,indsy)+symrel_cart(3,1,jsym)*magngz(2,indsy) 1810 ! myr=symrel_cart(1,2,jsym)*magngx(1,indsy)+symrel_cart(2,2,jsym)*magngy(1,indsy)+symrel_cart(3,2,jsym)*magngz(1,indsy) 1811 ! myi=symrel_cart(1,2,jsym)*magngx(2,indsy)+symrel_cart(2,2,jsym)*magngy(2,indsy)+symrel_cart(3,2,jsym)*magngz(2,indsy) 1812 ! mzr=symrel_cart(1,3,jsym)*magngx(1,indsy)+symrel_cart(2,3,jsym)*magngy(1,indsy)+symrel_cart(3,3,jsym)*magngz(1,indsy) 1813 ! mzi=symrel_cart(1,3,jsym)*magngx(2,indsy)+symrel_cart(2,3,jsym)*magngy(2,indsy)+symrel_cart(3,3,jsym)*magngz(2,indsy) 1814 1815 magxsu1=magxsu1+mxr*phr-mxi*phi;magxsu2=magxsu2+mxi*phr+mxr*phi 1816 magysu1=magysu1+myr*phr-myi*phi;magysu2=magysu2+myi*phr+myr*phi 1817 magzsu1=magzsu1+mzr*phr-mzi*phi;magzsu2=magzsu2+mzi*phr+mzr*phi 1818 end if 1819 end do 1820 rhosu1_arr(3*izone-2)=magxsu1/dble(nsym_used) 1821 rhosu1_arr(3*izone-1)=magysu1/dble(nsym_used) 1822 rhosu1_arr(3*izone )=magzsu1/dble(nsym_used) 1823 rhosu2_arr(3*izone-2)=magxsu2/dble(nsym_used) 1824 rhosu2_arr(3*izone-1)=magysu2/dble(nsym_used) 1825 rhosu2_arr(3*izone )=magzsu2/dble(nsym_used) 1826 numpt=numpt+nup 1827 end do 1828 1829 ! Reduction in case of FFT parallelization 1830 if(mpi_enreg%nproc_fft>1)then 1831 spaceComm=mpi_enreg%comm_fft 1832 call xmpi_sum(rhosu1_arr,spaceComm,ier) 1833 call xmpi_sum(rhosu2_arr,spaceComm,ier) 1834 end if 1835 1836 ! Now symmetrize the magnetization at equivalent points 1837 ! Use Transpose[symrel] 1838 numpt=0 1839 do izone=1,izone_max 1840 rep=irrzon(izone,2,imagn) 1841 nup=nsym_used/rep 1842 do iup=1,nup 1843 ind=irrzon(iup+numpt,1,imagn) 1844 ! Get coordinates in the range [0,n-1] 1845 j=ind-1;j1=modulo(j,n1);j2=modulo(j/n1,n2);j3=j/(n1*n2) 1846 !TO BE UNCOMMENTED 1847 ! Get location of G vector (grid point) centered at 0 0 0 1848 l3=j3-(j3/id3)*n3 1849 l2=j2-(j2/id2)*n2 1850 l1=j1-(j1/id1)*n1 1851 if(fftn2_distrib(j2+1)==me_fft) then ! this ind is to be treated by me_fft 1852 r2=ffti2_local(j2+1) - 1 1853 ind=n1*(nd2*j3+r2)+j1+1 ! this is ind in the current proc 1854 jsym=isymg(iup+numpt) 1855 if (jsym==0) then 1856 ABI_ERROR("jsym=0 in symrhg !") 1857 end if 1858 magxsu1=rhosu1_arr(3*izone-2);magxsu2=rhosu2_arr(3*izone-2) 1859 magysu1=rhosu1_arr(3*izone-1);magysu2=rhosu2_arr(3*izone-1) 1860 magzsu1=rhosu1_arr(3*izone );magzsu2=rhosu2_arr(3*izone ) 1861 ! Working on this: the present coding will be detrimental for speed ! cos and sin are recomputed many times ! 1862 tau1=tnons_used(1,jsym) 1863 tau2=tnons_used(2,jsym) 1864 tau3=tnons_used(3,jsym) 1865 if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then 1866 ! Compute exp(-2*Pi*I*G dot tau) using original G (equivalent of phnons in the collinear case) 1867 arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3) 1868 phr=cos(arg) 1869 phi=-sin(arg) 1870 else 1871 phr=one 1872 phi=zero 1873 end if 1874 phr=phr*symafm_used(jsym) 1875 phi=phi*symafm_used(jsym) 1876 !TO BE COMMENTED 1877 ! phr=phnons(1,iup,imagn);if (rep==1) phr=phr*symafm_used(jsym) !if rep==2, symafm is already included in phnons 1878 ! phi=phnons(2,iup,imagn);if (rep==1) phi=phi*symafm_used(jsym) !(see irrzg.F90) 1879 ! The magnetization should transform as a vector in real space 1880 ! => symrel_cart ?!? 1881 mxr=symrec_cart(1,1,jsym)*magxsu1+symrec_cart(2,1,jsym)*magysu1+symrec_cart(3,1,jsym)*magzsu1 1882 mxi=symrec_cart(1,1,jsym)*magxsu2+symrec_cart(2,1,jsym)*magysu2+symrec_cart(3,1,jsym)*magzsu2 1883 myr=symrec_cart(1,2,jsym)*magxsu1+symrec_cart(2,2,jsym)*magysu1+symrec_cart(3,2,jsym)*magzsu1 1884 myi=symrec_cart(1,2,jsym)*magxsu2+symrec_cart(2,2,jsym)*magysu2+symrec_cart(3,2,jsym)*magzsu2 1885 mzr=symrec_cart(1,3,jsym)*magxsu1+symrec_cart(2,3,jsym)*magysu1+symrec_cart(3,3,jsym)*magzsu1 1886 mzi=symrec_cart(1,3,jsym)*magxsu2+symrec_cart(2,3,jsym)*magysu2+symrec_cart(3,3,jsym)*magzsu2 1887 ! mxr=symrel_cart(1,1,jsym)*magxsu1+symrel_cart(1,2,jsym)*magysu1+symrel_cart(1,3,jsym)*magzsu1 1888 ! mxi=symrel_cart(1,1,jsym)*magxsu2+symrel_cart(1,2,jsym)*magysu2+symrel_cart(1,3,jsym)*magzsu2 1889 ! myr=symrel_cart(2,1,jsym)*magxsu1+symrel_cart(2,2,jsym)*magysu1+symrel_cart(2,3,jsym)*magzsu1 1890 ! myi=symrel_cart(2,1,jsym)*magxsu2+symrel_cart(2,2,jsym)*magysu2+symrel_cart(2,3,jsym)*magzsu2 1891 ! mzr=symrel_cart(3,1,jsym)*magxsu1+symrel_cart(3,2,jsym)*magysu1+symrel_cart(3,3,jsym)*magzsu1 1892 ! mzi=symrel_cart(3,1,jsym)*magxsu2+symrel_cart(3,2,jsym)*magysu2+symrel_cart(3,3,jsym)*magzsu2 1893 magngx(1,ind)=mxr*phr-mxi*phi 1894 magngx(2,ind)=mxi*phr+mxr*phi 1895 magngy(1,ind)=myr*phr-myi*phi 1896 magngy(2,ind)=myi*phr+myr*phi 1897 magngz(1,ind)=mzr*phr-mzi*phi 1898 magngz(2,ind)=mzi*phr+mzr*phi 1899 end if 1900 end do 1901 numpt=numpt+nup 1902 end do 1903 ABI_FREE(isymg) 1904 ABI_FREE(rhosu1_arr) 1905 ABI_FREE(rhosu2_arr) 1906 ABI_FREE(symrec_cart) 1907 ABI_FREE(symrel_cart) 1908 ABI_FREE(symafm_used) 1909 ABI_FREE(tnons_used) 1910 1911 end if ! nspden==4 1912 1913 call timab(17,2,tsec) 1914 1915 ! Pull out full or spin up density, now symmetrized 1916 call fourdp(cplex,rhog,work,1,mpi_enreg,nfft,1,ngfft,0) 1917 rhor(:,ispden)=work(:) 1918 if (nspden==4) then 1919 call fourdp(cplex,magngx,work,1,mpi_enreg,nfft,1,ngfft,0) 1920 rhor(:,2)=work(:) 1921 call fourdp(cplex,magngy,work,1,mpi_enreg,nfft,1,ngfft,0) 1922 rhor(:,3)=work(:) 1923 call fourdp(cplex,magngz,work,1,mpi_enreg,nfft,1,ngfft,0) 1924 rhor(:,4)=work(:) 1925 ABI_FREE(magngx) 1926 ABI_FREE(magngy) 1927 ABI_FREE(magngz) 1928 end if 1929 1930 end do ! ispden 1931 1932 end if ! End on the condition nsym==1 1933 1934 ABI_FREE(work) 1935 1936 contains 1937 1938 function map_symrhg(j1,n1) 1939 1940 integer :: map_symrhg 1941 integer,intent(in) :: j1,n1 1942 ! Map into [0,n-1] 1943 map_symrhg=mod(n1+mod(j1,n1),n1) 1944 end function map_symrhg 1945 1946 end subroutine symrhg