TABLE OF CONTENTS
- ABINIT/m_inwffil
- m_inwffil/cg_from_atoms
- m_inwffil/initwf
- m_inwffil/inwffil
- m_inwffil/newkpt
- m_inwffil/pareigocc
- m_inwffil/wfconv
- m_inwffil/wfsinp
ABINIT/m_inwffil [ Modules ]
NAME
m_inwffil
FUNCTION
Initialization of wavefunctions.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, AR, MB, MVer, ZL, MB, TD, MG) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_inwffil 23 24 use defs_basis 25 use defs_wvltypes 26 use m_abicore 27 use m_wffile 28 use m_wfk 29 use m_errors 30 use m_xomp 31 use m_xmpi 32 use m_nctk 33 use m_hdr 34 use m_dtset 35 #if defined HAVE_MPI2 36 use mpi 37 #endif 38 39 use defs_abitypes, only : MPI_type 40 use m_fstrings, only : sjoin, itoa 41 use m_time, only : timab, cwtime, cwtime_report 42 use m_io_tools, only : file_exists, get_unit 43 use m_geometry, only : getspinrot 44 use m_pptools, only : prmat 45 use m_symtk, only : matr3inv, mati3inv 46 use m_cgtools, only : cg_envlop, pw_orthon 47 use m_fftcore, only : kpgsph, sphere, sphereboundary 48 use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_io 49 use m_mpinfo, only : destroy_mpi_enreg, copy_mpi_enreg, proc_distrb_cycle 50 use m_kg, only : kpgio, ph1d3d, getph 51 use m_kpts, only : listkk 52 use m_rwwf, only : rwwf, WffReadSkipK 53 use m_wvl_wfsinp, only : wvl_wfsinp_disk, wvl_wfsinp_scratch 54 55 implicit none 56 57 #if defined HAVE_MPI1 58 include 'mpif.h' 59 #endif 60 61 private
m_inwffil/cg_from_atoms [ Functions ]
[ Top ] [ m_inwffil ] [ Functions ]
NAME
cg_from_atoms
FUNCTION
Initialize wave functions at a given (k-point, spin) using Bloch sums of atomic orbitals.
INPUTS
ikpt,isppol=k-point index, spin index rprimd(3,3)=Direct lattice vectors in Bohr. xred(3,natom)=Atomic positions. kg_k(3,npw_k)=reduced planewave coordinates. dtset <type(dataset_type)>=all input variables for this dataset gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k mpi_enreg=information about MPI parallelization nband=number of bands at this k point for that spin polarization npw=number of plane waves at this k point my_nspinor=number of spinors treated by this MPI proc
OUTPUT
eig(nband)=array for holding eigenvalues (hartree)
SIDE EFFECTS
cg(2,*)=updated wavefunctions
SOURCE
3423 subroutine cg_from_atoms(ikpt, isppol, rprimd, xred, kg_k, cg, dtset, psps, eig, gs_hamk, & 3424 mpi_enreg, nband, npw, my_nspinor) 3425 3426 use defs_datatypes, only : pseudopotential_type 3427 use m_geometry, only : metric 3428 use m_splines, only : splfit 3429 use m_initylmg, only : initylmg_k 3430 use m_hamiltonian, only : gs_hamiltonian_type 3431 use m_getghc, only : getghc 3432 use m_nonlop, only : nonlop 3433 use m_pawcprj, only : pawcprj_type 3434 use m_rmm_diis, only : subspace_rotation 3435 use m_cgtools, only : cgpaw_cholesky, cgnc_cholesky 3436 3437 !Arguments ------------------------------------ 3438 integer,intent(in) :: ikpt, isppol, nband, npw, my_nspinor 3439 type(gs_hamiltonian_type),intent(inout) :: gs_hamk 3440 type(dataset_type),intent(in) :: dtset 3441 type(pseudopotential_type),intent(in) :: psps 3442 type(mpi_type),intent(inout) :: mpi_enreg 3443 real(dp),intent(in) :: rprimd(3,3), xred(3,dtset%natom) 3444 integer, intent(in) :: kg_k(3,npw) 3445 real(dp),intent(inout) :: cg(2,npw*my_nspinor,nband) 3446 real(dp),intent(out) :: eig(nband) 3447 3448 !Local variables------------------------------- 3449 integer,parameter :: optder0 = 0, ider0 = 0, icg0 = 0 3450 !integer,parameter :: idir0 = 0, !, type_calc0 = 0, option1 = 1, option2 = 2, tim_getghc = 0 3451 integer :: npwsp !, ortalgo, ierr, 3452 integer :: istwf_k, usepaw, mcg !, mgsc 3453 integer :: me_g0, me_cell !prev_mixprec, 3454 integer :: comm_bsf, savemem, ll 3455 integer :: iatom, itypat, iln, ig, iband, ilmn, ilm, im, lnmax 3456 real(dp) :: kpg1, kpg2, kpg3, kpgc1, kpgc2, kpgc3 3457 complex(dp) :: cfact 3458 logical :: supported !, use_fft_mixprec 3459 real(dp) :: ucvol, arg ! cpu, wall, gflops, 3460 !character(len=500) :: msg 3461 !arrays 3462 real(dp) :: gmet(3,3), gprimd(3,3), rmet(3,3), kpt(3), phase_l(2), ri(2) 3463 real(dp) :: enlx(nband) !, tsec(2) 3464 real(dp),allocatable :: ghc(:,:), gvnlxc(:,:) 3465 real(dp),allocatable :: kpg_k(:,:), tphiq(:,:,:), sf(:,:) !,ph3d(:,:,:) 3466 real(dp),allocatable :: ylm(:,:), ylm_gr(:,:,:), gsc(:,:) 3467 real(dp),allocatable :: kpgnorm(:), wk_ffnl2(:) 3468 3469 ! ************************************************************************* 3470 3471 ! Define useful vars. 3472 usepaw = dtset%usepaw; istwf_k = gs_hamk%istwf_k 3473 me_g0 = mpi_enreg%me_g0; comm_bsf = mpi_enreg%comm_bandspinorfft 3474 npwsp = npw * my_nspinor; mcg = npwsp * nband !; mgsc = npwsp * nband * usepaw 3475 me_cell = mpi_enreg%me_cell 3476 kpt = dtset%kptns(:,ikpt) 3477 3478 supported = .True. 3479 if (dtset%usepaw == 0) then 3480 lnmax = maxval(psps%nctab(:)%num_tphi) 3481 supported = supported .and. minval(psps%nctab(:)%num_tphi) > 0 3482 end if 3483 3484 ! Test whether cg initialization from ps atomic orbitals is coded/supported. 3485 if (dtset%nspinor == 2) supported = .False. 3486 if (dtset%usepaw /= 0) supported = .False. 3487 if (.not. supported) then 3488 if (me_cell == 0 .and. ikpt == 1) then 3489 call wrtout(std_out, " cg initialization from atomic orbitals not available. returning") 3490 end if 3491 return 3492 end if 3493 3494 if (me_cell == 0 .and. ikpt == 1) then 3495 call wrtout(std_out, sjoin(" Initializing cg from atomic orbitals for ikpt:", itoa(ikpt), ", spin:", itoa(isppol))) 3496 end if 3497 !call cwtime(cpu, wall, gflops, "start") 3498 3499 call metric(gmet, gprimd, -1, rmet, rprimd, ucvol) 3500 3501 ABI_MALLOC(ylm, (npw, psps%mpsang**2)) 3502 ABI_MALLOC(ylm_gr, (npw, 3+6*(optder0/2), psps%mpsang**2)) 3503 3504 call initylmg_k(npw, psps%mpsang, optder0, rprimd, gprimd, kpt, kg_k, ylm, ylm_gr) 3505 ABI_SFREE(ylm_gr) 3506 3507 ! Compute nonlocal form factors at (k+G) 3508 ! Note that we need to work with useylm = 1 to keep the m-dependency 3509 ! even when Vnl is applied with Legendre polynomials (useylm = 0) 3510 3511 ! Get |k+G| 3512 ABI_MALLOC(kpgnorm, (npw)) 3513 3514 !$OMP PARALLEL DO PRIVATE(kpg1, kpg2, kpg3, kpgc1, kpgc2, kpgc3) 3515 do ig=1,npw 3516 kpg1=kpt(1)+dble(kg_k(1,ig)) 3517 kpg2=kpt(2)+dble(kg_k(2,ig)) 3518 kpg3=kpt(3)+dble(kg_k(3,ig)) 3519 kpgc1=kpg1*gprimd(1,1)+kpg2*gprimd(1,2)+kpg3*gprimd(1,3) 3520 kpgc2=kpg1*gprimd(2,1)+kpg2*gprimd(2,2)+kpg3*gprimd(2,3) 3521 kpgc3=kpg1*gprimd(3,1)+kpg2*gprimd(3,2)+kpg3*gprimd(3,3) 3522 kpgnorm(ig)=sqrt(kpgc1*kpgc1+kpgc2*kpgc2+kpgc3*kpgc3) 3523 end do 3524 3525 ABI_MALLOC(tphiq, (npw, lnmax, psps%ntypat)) 3526 ABI_MALLOC(wk_ffnl2, (npw)) 3527 3528 do itypat=1,psps%ntypat 3529 do iln=1,psps%nctab(itypat)%num_tphi 3530 call splfit(psps%qgrid_ff, wk_ffnl2, psps%nctab(itypat)%tphi_qspl(:,:,iln), & 3531 ider0, kpgnorm, tphiq(:,iln,itypat), psps%mqgrid_ff, npw) 3532 end do 3533 end do 3534 3535 ABI_FREE(kpgnorm) 3536 ABI_FREE(wk_ffnl2) 3537 3538 !call getph(atindx, natom, n1, n2, n3, ph1d, xred) 3539 !call ph1d3d(iatom, jatom, kg_k, matblk, natom, npw_k, n1, n2, n3, phkxred, ph1d, ph3d) 3540 3541 ! Now init cg. Assuming cg has been already filled with random numbers previously 3542 ! so we only need to init the first states. We don't take into account the occupancies 3543 ! in the isolated atom. We just loop over all nlm states until we have filled max nband states. 3544 ABI_MALLOC(sf, (2, npw)) 3545 iband = 0 3546 iatom_loop: do iatom=1,dtset%natom 3547 itypat = dtset%typat(iatom) 3548 3549 ! Structure factor. 3550 do ig=1,npw 3551 arg = -two_pi * dot_product(xred(:,iatom), kpt + kg_k(:,ig)) 3552 sf(1,ig) = cos(arg) 3553 sf(2,ig) = sin(arg) 3554 end do 3555 3556 ilmn = 0 3557 do iln=1,psps%nctab(itypat)%num_tphi 3558 !if (psps%nctab(itypat)%tphi_occ(iln) < zero) cycle 3559 ll = psps%nctab(itypat)%tphi_l(iln) 3560 !cfact = (j_dpc ** ll) * four_pi / sqrt(ucvol) 3561 cfact = (-j_dpc ** ll) * four_pi / sqrt(ucvol) 3562 phase_l(1) = dble(cfact) 3563 phase_l(2) = aimag(cfact) 3564 do im=1, 2*ll+1 3565 ilmn = ilmn + 1 3566 ilm = im + ll**2 3567 iband = iband + 1 3568 ! Another good reason why nband should be > nbocc. 3569 if (iband > nband) exit iatom_loop 3570 3571 ! NB: Assuming nspinor == 1 3572 do ig=1,npw ! *my_nspinor 3573 ri(1) = phase_l(1) * sf(1, ig) - phase_l(2) * sf(2, ig) 3574 ri(2) = phase_l(1) * sf(2, ig) + phase_l(2) * sf(1, ig) 3575 if (dtset%wfinit == 1) call randomize(ri) 3576 cg(:, ig, iband) = ri(:) * ylm(ig, ilm) * tphiq(ig, iln, itypat) 3577 !wfcatom (ig, 1, n_starting_wfc) = phase_l * sf(1, ig) * ylm(ig, ilm) * chiq(ig, nb, nt) 3578 end do ! ig 3579 3580 ! XG030513: Time-reversal symmetry for k=gamma imposes zero imaginary part at G=0 3581 ! XG: I do not know what happens for spin-orbit here. 3582 if (istwf_k == 2 .and. mpi_enreg%me_g0 == 1) cg(2, 1, iband) = zero 3583 end do ! im 3584 end do ! iln 3585 end do iatom_loop 3586 3587 !call cg_envlop(cg, dtset%ecut, gmet, icg0, kg_k, kpt, mcg, nband, npw, my_nspinor) 3588 3589 ABI_FREE(sf) 3590 ABI_FREE(tphiq) 3591 ABI_SFREE(kpg_k) 3592 ABI_FREE(ylm) 3593 3594 ! Use mixed precisions if requested by the user but only for low accuracy_level 3595 !use_fft_mixprec = dtset%mixprec == 1 .and. accuracy_level < 2 3596 !if (use_fft_mixprec) prev_mixprec = fftcore_set_mixprec(1) 3597 3598 ! ========================= 3599 ! === Subspace rotation === 3600 ! ========================= 3601 savemem = 1 3602 ABI_MALLOC(gsc, (2, npw*my_nspinor*nband*dtset%usepaw)) 3603 call subspace_rotation(gs_hamk, dtset%prtvol, mpi_enreg, nband, npw, my_nspinor, savemem, & 3604 enlx, eig, cg, gsc, ghc, gvnlxc) 3605 3606 ABI_SFREE(ghc) 3607 ABI_SFREE(gvnlxc) 3608 3609 ! Revert mixprec to previous status before returning. 3610 !if (use_fft_mixprec) prev_mixprec = fftcore_set_mixprec(prev_mixprec) 3611 3612 ! Ortoghonalization is in principle not needed but it seems to improve a bit. 3613 if (dtset%wfinit < 0) then 3614 !ortalgo = 3 !; ortalgo = mpi_enreg%paral_kgb 3615 !call pw_orthon(0, 0, istwf_k, mcg, mgsc, npwsp, nband, ortalgo, gsc, usepaw, cg, me_g0, comm_bsf) 3616 3617 ! TODO: Merge the two routines. 3618 if (usepaw == 1) then 3619 call cgpaw_cholesky(npwsp, nband, cg, gsc, istwf_k, me_g0, comm_bsf) 3620 else 3621 call cgnc_cholesky(npwsp, nband, cg, istwf_k, me_g0, comm_bsf, use_gemm=.False.) 3622 end if 3623 end if 3624 3625 ABI_FREE(gsc) 3626 !call cwtime_report(" cg_from_atoms:", cpu, wall, gflops) 3627 3628 contains 3629 subroutine randomize(ri) 3630 real(dp),intent(inout) :: ri(2) 3631 3632 !Local variables------------------------------- 3633 real(dp) :: arg, rr 3634 complex(dp) :: c2, c1 3635 ! ************************************************************************* 3636 3637 call random_number(arg) 3638 arg = two_pi * arg 3639 call random_number(rr) 3640 c1 = cmplx(ri(1), ri(2), kind=dp) 3641 c2 = one + 0.05_dp * cmplx(rr*cos(arg), rr*sin(arg), kind=dp) 3642 c2 = c1 * c2 3643 ri(1) = real(c2) 3644 ri(2) = aimag(c2) 3645 3646 end subroutine randomize 3647 3648 end subroutine cg_from_atoms
m_inwffil/initwf [ Functions ]
[ Top ] [ m_inwffil ] [ Functions ]
NAME
initwf
FUNCTION
Initialization of wavefunctions.
If formeig==1, and partially filled case, I am not sure that the eig_k are initialized properly ...
formeig option (format of the eigenvalues and eigenvector) :
0 => ground-state format (initialisation of
eigenvectors with random numbers, vector of eigenvalues)
1 => respfn format (initialisation of
eigenvectors with 0 s, hermitian matrix of eigenvalues)
INPUTS
formeig=see above headform=header format (might be needed to read the block of wfs) icg=shift to be given to the location of the data in the array cg ikpt= number of the k point of which the wf is initialised spin=spin index mcg=dimension of the cg array mpi_enreg=information about MPI parallelization nband_k=number of bands at this particular k point nkpt=number of k points npw=number of plane waves nspinor=number of spinorial components of the wavefunctions (on current proc) wff1=structure info for file containing wavefunctions (when needed)
OUTPUT
cg(2,mcg)=complex wf array
if ground state format (formeig=0):
eig_k(nband_k)=list of eigenvalues (input or init to large number), hartree
if respfn format (formeig=1):
eig_k(2*nband_k*nband_k)= matrix of eigenvalues (input or init to large number), hartree
SIDE EFFECTS
Input/output: occ_k(nband_k)=list of occupations (input or left to their initial value) ikptsp_old=number of the previous spin-k point, or 0 if first call of present file
SOURCE
1820 subroutine initwf(cg,eig_k,formeig,headform,icg,ikpt,ikptsp_old,& 1821 & spin,mcg,mpi_enreg,nband_k,nkpt,npw,nspinor,occ_k,wff1) 1822 1823 !Arguments ------------------------------------ 1824 !scalars 1825 integer,intent(in) :: formeig,headform,icg,ikpt,spin,mcg,nband_k,nkpt,npw,nspinor 1826 integer,intent(inout) :: ikptsp_old 1827 type(MPI_type),intent(in) :: mpi_enreg 1828 type(wffile_type),intent(inout) :: wff1 1829 !arrays 1830 real(dp),intent(inout) :: occ_k(nband_k) 1831 real(dp),intent(inout) :: cg(2,mcg),eig_k((2*nband_k)**formeig*nband_k) !vz_i 1832 1833 !Local variables------------------------------- 1834 !scalars 1835 integer,parameter :: nkpt_max=50 1836 integer :: ikpt0,nband_disk,tim_rwwf 1837 character(len=500) :: msg 1838 !arrays 1839 integer,allocatable :: kg_dum(:,:) 1840 real(dp) :: tsec(2) 1841 #if 0 1842 integer :: iomode,comm,funt,ierr 1843 type(wfk_t) :: Wfk 1844 #endif 1845 1846 ! ************************************************************************* 1847 1848 !write(std_out,*)' initwf : enter, ikptsp_old,ikpt,spin,nkpt= ',ikptsp_old,ikpt,spin,nkpt 1849 !stop 1850 1851 #if 0 1852 ABI_WARNING("Entering new IO section") 1853 !call WffClose(wff1,ierr) 1854 comm = MPI_enreg%comm_cell 1855 iomode = iomode_from_fname(wff1%fname) 1856 call wfk_open_read(Wfk,wff1%fname,formeig,iomode,get_unit(),comm) 1857 call wfk_read_band_block(Wfk,(/1,nband_k/),ikpt,spin,xmpio_at,cg_k=cg(1:,icg:),eig_k=eig_k,occ_k=occ_k) 1858 call wfk_close(Wfk) 1859 !call clsopn(wff1) 1860 RETURN 1861 #endif 1862 1863 call timab(770,1,tsec) 1864 call timab(771,1,tsec) 1865 1866 ABI_MALLOC(kg_dum,(3,0)) 1867 1868 !Skip wavefunctions for k-points not treated by this proc. 1869 !(from ikptsp_old+1 to ikpt+(spin-1)*nkpt-1) 1870 if (ikptsp_old<ikpt+(spin-1)*nkpt-1) then 1871 do ikpt0=ikptsp_old+1,ikpt+(spin-1)*nkpt-1 1872 call WffReadSkipK(formeig,headform,ikpt0,spin,mpi_enreg,wff1) 1873 end do 1874 end if 1875 1876 !write(std_out,*)' initwf : before rwwf' 1877 !write(std_out,*)' formeig,icg,ikpt,spin=',formeig,icg,ikpt,spin 1878 !write(std_out,*)' nband_k,nband_disk,npw,nspinor=',nband_k,nband_disk,npw,nspinor 1879 !write(std_out,*)' unwff1=',unwff1 1880 !stop 1881 1882 if(mpi_enreg%paralbd==0)tim_rwwf=2 1883 if(mpi_enreg%paralbd==1)tim_rwwf=20 1884 1885 call timab(771,2,tsec) 1886 1887 call rwwf(cg,eig_k,formeig,headform,icg,ikpt,spin,kg_dum,nband_k,mcg,mpi_enreg,nband_k,nband_disk,& 1888 & npw,nspinor,occ_k,1,0,tim_rwwf,wff1) 1889 1890 call timab(772,1,tsec) 1891 1892 if (ikpt<=nkpt_max) then 1893 write(msg,'(3(a,i0))')' initwf: disk file gives npw= ',npw,' nband= ',nband_disk,' for kpt number= ',ikpt 1894 call wrtout(std_out,msg) 1895 else if (ikpt==nkpt_max+1) then 1896 call wrtout(std_out,' initwf: the number of similar message is sufficient... stop printing them') 1897 end if 1898 1899 ! Check the number of bands on disk file against desired number. These are not required to agree) 1900 if (nband_disk /= nband_k .and. ikpt<=nkpt_max) then 1901 write(msg,'(3(a,i0),3a)')& 1902 'For kpt number: ',ikpt,' disk file has: ',nband_disk,' bands but input file gave nband: ',nband_k,'.',ch10,& 1903 'This is not fatal. Bands are skipped or filled with random numbers.' 1904 ABI_COMMENT(msg) 1905 end if 1906 1907 if (ikpt<=nkpt_max) then 1908 write(msg,'(a,i0,a)')' initwf: ',nband_disk,' bands have been initialized from disk' 1909 call wrtout(std_out,msg) 1910 end if 1911 1912 ikptsp_old=ikpt+(spin-1)*nkpt 1913 1914 ABI_FREE(kg_dum) 1915 1916 call timab(772,2,tsec) 1917 call timab(770,2,tsec) 1918 1919 end subroutine initwf
m_inwffil/inwffil [ Functions ]
[ Top ] [ m_inwffil ] [ Functions ]
NAME
inwffil
FUNCTION
Do initialization of wavefunctions. Also call other relevant routines for this initialisation (initialization of wavefunctions from scratch or from file, translations of wavefunctions, ...)
INPUTS
ask_accurate= if 1, the wavefunctions and eigenvalues must be
accurate, that is, they must come from a k point that is
symmetric of the needed k point, with a very small tolerance,
the disk file contained sufficient bands to initialize all of them,
the spinor and spin-polarisation characteristics must be identical
dtset <type(dataset_type)>=all input variables for this dataset
ecut=effective kinetic energy planewave cutoff (hartree), beyond
which the coefficients of plane waves are zero
ecut_eff=effective kinetic energy planewave cutoff (hartree), needed
to generate the sphere of plane wave
exchn2n3d=if 1, n2 and n3 are exchanged
formeig=explained above
hdr <type(hdr_type)>=the header of wf, den and pot files
ireadwf=option parameter described above for wf initialization
istwfk(nkpt)=input option parameter that describes the storage of wfs to be initialized here.
kg(3,mpw*my_nkpt)=dimensionless coords of G vecs in basis sphere at k point
kptns(3,nkpt)=reduced coords of k points
localrdwf=(for parallel case) if 1, the wffnm file is local to each machine
mband=maximum number of bands
mband_mem=maximum number of bands for this cpu
mcg=size of wave-functions array (cg) =mpw*nspinor*mband_mem*mkmem*nsppol
mkmem=number of k-points in core memory
mpi_enreg=information about MPI parallelization
mpw=maximum number of planewaves as dimensioned in calling routine
nband(nkpt*nsppol)=number of bands at each k point
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
nkpt=number of k points
npwarr(nkpt)=array holding npw for each k point.
nsppol=1 for unpolarized, 2 for spin-polarized
nsym=number of symmetry elements in space group
occ(mband*nkpt*nsppol)=occupations (from disk or left at their initial value)
optorth= 1 if the WFS have to be orthogonalized; 0 otherwise
prtvol=control print volume and debugging
symafm(nsym)=(anti)ferromagnetic part of symmetry operations
symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
tnons(3,nsym)=nonsymmorphic translations for symmetry operations
unkg=unit number for storage of basis sphere data: stores indirect
indexing array and integer coordinates for all planewaves in basis
sphere for each k point being considered
unwff1,unwfnow= unit numbers for files wffnm and wft1nm.
wffnm=name (character data) of file for input wavefunctions.
OUTPUT
wff1 = structure information for files wffnm .
wffnow= structure information for wf file wft1nm
if ground state format (formeig=0):
eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
if respfn format (formeig=1):
eigen(2*mband*mband*nkpt*nsppol)=matrix of eigenvalues
(input or init to large number), (Ha)
Conditional output (returned if mkmem/=0):
cg(2,mcg)=complex wf array
be careful : an array of size cg(2,npw*nspinor), as used
in the response function code, is not enough !
wvl <type(wvl_data)>=all wavelets data.
NOTES
Detailed description:
Initialize unit wff1%unwff for input of wf data if ireadwf=1
Opens file on unit wffnow%unwff
if the storage on disk is needed (mkmem==0)
Initializes wf data on wffnow%unwff, by calling the appropriate routine.
formeig option (format of the eigenvalues and occupations) :
0 => ground-state format (initialisation of
eigenvectors with random numbers, vector of eigenvalues,
occupations are present)
1 => respfn format (initialisation of
eigenvectors with 0 s, hermitian matrix of eigenvalues)
ireadwf options:
0 => initialize with random numbers or 0 s
1 => read from disk file wff1, initializing higher bands
with random numbers or 0 s if not provided in disk file
The wavefunctions after this initialisation are stored in unit wffnow%unwff
WARNINGS
* The symmetry operations are used to translate the data from one k point to another, symmetric, k point. They can be completely different from the symmetry operations contained on the disk file. No check is performed between the two sets. * Occupations will not be modified nor output, in the present status of this routine. * If ground state format (formeig=0) occ(mband*nkpt*nsppol) was output. NOT OUTPUT NOW!
SOURCE
172 subroutine inwffil(ask_accurate,cg,dtset,ecut,ecut_eff,eigen,exchn2n3d,& 173 & formeig,hdr,ireadwf,istwfk,kg,kptns,localrdwf,mband,& 174 & mcg,mkmem,mpi_enreg,mpw,nband,ngfft,nkpt,npwarr,& 175 & nsppol,nsym,occ,optorth,symafm,symrel,tnons,unkg,wff1,& 176 & wffnow,unwff1,wffnm,wvl) 177 178 !Arguments ------------------------------------ 179 integer,intent(in) :: ask_accurate,exchn2n3d,formeig,ireadwf,localrdwf,mband,mcg,mkmem,mpw 180 integer,intent(in) :: nkpt,nsppol,nsym,optorth,unkg,unwff1 181 real(dp),intent(in) :: ecut,ecut_eff 182 character(len=*),intent(in) :: wffnm 183 type(MPI_type),intent(inout),target :: mpi_enreg 184 type(dataset_type),intent(in) :: dtset 185 type(hdr_type),intent(inout) :: hdr 186 type(wffile_type),intent(inout) :: wff1 187 type(wffile_type),intent(inout) :: wffnow 188 type(wvl_data),intent(inout) :: wvl 189 integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),ngfft(18) 190 integer,intent(in) :: npwarr(nkpt),symafm(nsym),symrel(3,3,nsym) 191 integer,intent(in),target :: nband(nkpt*nsppol) 192 real(dp),intent(inout),target :: cg(2,mcg),eigen((2*mband)**formeig*mband*nkpt*nsppol) 193 real(dp),intent(in) :: kptns(3,nkpt),tnons(3,nsym) 194 real(dp),intent(inout) :: occ(mband*nkpt*nsppol) 195 196 !Local variables------------------------------- 197 integer,parameter :: master=0 198 integer :: iomode,accurate,ceksp,debug,doorth,fform,fform_dum,fill 199 integer :: headform0,iband,ibg,ibg0,icg,icg0,icgsft,ieigsft,ierr,ii 200 integer :: ikassoc,ikpt,ikpt0,ikptsp,ikptsp0,imax,increase_nkassoc,isppol,isppol0 201 integer :: mband0,mband0_rd,mband_eff,mcg_disk,me,me0,mkmem0,mpw0 202 integer :: my_nkpt,my_nspinor,my_nspinor0,nband_k,nband0_k 203 integer :: nkassoc,nkpt0,npw,npw0,nspinor0,nspinor_eff,nsppol0,nsppol_eff,nsppol2nspinor 204 integer :: rdwr,randalg,restart,restartpaw,spaceComm,spaceComm_io,sppoldbl,sppoldbl_eff,squeeze 205 logical :: out_of_core 206 real(dp) :: dksqmax,ecut0 207 character(len=500) :: message 208 type(hdr_type) :: hdr0 209 integer :: ngfft0(18) 210 integer,allocatable :: indkk0(:,:),indx(:),istwfk0(:),kg0(:,:) 211 integer,allocatable :: nband0_rd(:),npwarr0(:),npwi(:),npwtot0(:) 212 integer,allocatable,target :: indkk(:,:),nband0(:) 213 integer, pointer :: indkk_eff(:,:),nband_eff(:) 214 logical,allocatable :: my_kpt(:) 215 real(dp) :: gmet(3,3),gmet0(3,3),gprim0(3,3),rprim0(3,3),tsec(2) 216 real(dp),allocatable :: cg_disk(:,:),kptns0(:,:) 217 real(dp),pointer :: cg_eff(:,:),eigen_eff(:) 218 type(MPI_type),pointer :: mpi_enreg0 219 220 ! ************************************************************************* 221 222 DBG_ENTER("COLL") 223 224 !Keep track of total time spent in inwffil 225 call timab(710,1,tsec) 226 call timab(711,1,tsec) 227 228 !Check the validity of formeig 229 if (formeig/=0.and.formeig/=1) then 230 write(message,'(a,i0,a)')' formeig = ',formeig,', but the only allowed values are 0 or 1.' 231 ABI_BUG(message) 232 end if 233 234 !Init mpi_comm 235 spaceComm=mpi_enreg%comm_cell 236 spaceComm_io=xmpi_comm_self 237 if (mpi_enreg%paral_kgb==1) spaceComm_io= mpi_enreg%comm_bandspinorfft 238 if (mpi_enreg%paral_hf ==1) spaceComm_io= mpi_enreg%comm_hf 239 me=xmpi_comm_rank(spaceComm) 240 241 !Determine number of k points processed by current node 242 my_nkpt=nkpt;if (size(mpi_enreg%my_kpttab)>0) my_nkpt=maxval(mpi_enreg%my_kpttab) 243 out_of_core=(mkmem==0.and.my_nkpt/=0) 244 245 ngfft0(:)=ngfft(:) 246 headform0=0 !Default value for headform0 (will be needed later, to read wf blocks) 247 248 !Chebyshev is more sensitive to the quality of input random numbers, so use a new algorithm 249 if(dtset%wfoptalg == 1 .or. dtset%wfoptalg == 111) then 250 randalg = 1 251 else 252 ! Otherwise, use compatibility mode 253 randalg = 0 254 end if 255 256 !If the input data are on disk, determine the kind of restart 257 wff1%fname = wffnm 258 259 !Checking the existence of data file 260 if (ireadwf==1 .and. .not.file_exists(wff1%fname)) then 261 ! Trick needed to run Abinit test suite in netcdf mode. 262 if (file_exists(nctk_ncify(wff1%fname))) then 263 write(std_out,"(3a)")"- File: ",trim(wff1%fname)," does not exist but found netcdf file with similar name." 264 wff1%fname = nctk_ncify(wff1%fname) 265 end if 266 if (localrdwf/=0 .and. .not. file_exists(wff1%fname)) then 267 ABI_ERROR('Missing data file: '//TRIM(wff1%fname)) 268 end if 269 end if 270 271 !Compute reciprocal space metric gmet 272 call matr3inv(hdr%rprimd,gprim0) ! gprim0 is used as temporary storage 273 gmet=matmul(transpose(gprim0),gprim0) 274 275 if (ireadwf==1)then 276 277 iomode=dtset%iomode 278 if (localrdwf==0) then 279 ! This is in case the wff file must be read by only the master proc 280 if (iomode /= IO_MODE_ETSF) iomode=IO_MODE_FORTRAN_MASTER 281 !iomode=IO_MODE_FORTRAN_MASTER 282 end if 283 284 call WffOpen(iomode,spaceComm,wff1%fname,ierr,wff1,master,me,unwff1,spaceComm_io) 285 286 ! Initialize hdr0 (sent to all procs), thanks to reading of wff1 287 rdwr=1 288 if ( ANY(wff1%iomode == (/IO_MODE_FORTRAN_MASTER, IO_MODE_FORTRAN, IO_MODE_MPI/) )) then 289 call hdr_io(fform_dum,hdr0,rdwr,wff1) 290 #ifdef HAVE_NETCDF 291 else if (wff1%iomode == IO_MODE_ETSF) then 292 call hdr_ncread(hdr0, wff1%unwff, fform_dum) 293 #endif 294 end if 295 296 ! Handle IO Error. 297 if (fform_dum == 0) then 298 write(message,"(4a)")& 299 & "hdr_io returned fform == 0 while trying to read the wavefunctions from file: ",trim(wff1%fname),ch10,& 300 & "This usually means that the file does not exist or that you don't have enough privileges to read it" 301 ABI_ERROR(message) 302 end if 303 304 call wrtout(std_out,' inwffil: examining the header of disk file: '//trim(wff1%fname),'COLL') 305 306 ! Check hdr0 versus hdr (and from now on ignore header consistency and write new info to header for each file) 307 if (dtset%usewvl == 0) then 308 ! wait for plane waves. 309 fform=2 310 else 311 ! wait for wavelets. 312 fform = 200 313 end if 314 call hdr_check(fform,fform_dum,hdr,hdr0,'PERS',restart,restartpaw) 315 316 nkpt0=hdr0%nkpt 317 nsppol0=hdr0%nsppol 318 headform0=hdr0%headform 319 320 write(message,'(2a)')'-inwffil : will read wavefunctions from disk file ',trim(wff1%fname) 321 call wrtout(std_out,message,'COLL') 322 call wrtout(ab_out,message,'COLL') 323 324 else 325 restart=1; restartpaw=0 326 327 ! Fill some data concerning an hypothetical file to be read 328 ! This is to allow the safe use of same routines than with ireadwf==1. 329 nkpt0=nkpt ; nsppol0=nsppol 330 end if ! end ireadwf 331 332 sppoldbl=1 333 if(minval(symafm(:))==-1)then 334 if(nsppol0==1 .and. nsppol==2)sppoldbl=2 335 end if 336 337 ABI_MALLOC(indkk,(nkpt*sppoldbl,6)) 338 ABI_MALLOC(istwfk0,(nkpt0)) 339 ABI_MALLOC(kptns0,(3,nkpt0)) 340 ABI_MALLOC(nband0,(nkpt0*nsppol0)) 341 ABI_MALLOC(npwarr0,(nkpt0)) 342 343 if(restart==2)then ! restart with translations 344 345 ecut0=hdr0%ecut_eff 346 istwfk0(1:nkpt0)=hdr0%istwfk(1:nkpt0) 347 kptns0(1:3,1:nkpt0)=hdr0%kptns(1:3,1:nkpt0) 348 nband0(1:nkpt0*nsppol0)=hdr0%nband(1:nkpt0*nsppol0) 349 ngfft0(1:3)=hdr0%ngfft(1:3) 350 npwarr0(1:nkpt0)=hdr0%npwarr(1:nkpt0) 351 nspinor0=hdr0%nspinor 352 rprim0(:,:)=hdr0%rprimd(:,:) 353 mpw0=maxval(npwarr0(:)) 354 355 ! Compute reciprocal space metric gmet for unit cell of disk wf 356 call matr3inv(rprim0,gprim0) 357 gmet0=matmul(transpose(gprim0),gprim0) 358 359 if ((mpi_enreg%paral_kgb==1).or.(mpi_enreg%paral_hf==1)) then 360 ABI_MALLOC(mpi_enreg0,) 361 call copy_mpi_enreg(mpi_enreg,mpi_enreg0) 362 ABI_MALLOC(kg0,(3,mpw0*nkpt0)) 363 ABI_MALLOC(npwtot0,(nkpt0)) 364 message="tmpfil" 365 call kpgio(ecut0,dtset%exchn2n3d,gmet0,istwfk0,kg0, & 366 & kptns0,nkpt0,nband0,nkpt0,'PERS',mpi_enreg0,& 367 & mpw0,npwarr0,npwtot0,nsppol0) 368 369 ABI_FREE(kg0) 370 ABI_FREE(npwtot0) 371 else 372 mpi_enreg0 => mpi_enreg 373 end if 374 375 ! At this stage, the header of the file wff1i%unwff is read, and 376 ! the pointer is ready to read the first wavefunction block. 377 378 ! Compute k points from input file closest to the output file 379 call listkk(dksqmax,gmet0,indkk,kptns0,kptns,nkpt0,nkpt,nsym,sppoldbl,symafm,symrel,1,spaceComm) 380 381 else if (restart==1) then ! direct restart 382 383 ! Fill variables that must be the same, as determined by hdr_check.f 384 ! This is to allow the safe use of the same routines than with restart==2. 385 nspinor0=dtset%nspinor 386 ecut0=ecut_eff 387 gmet0(:,:)=gmet(:,:) 388 istwfk0(:)=istwfk(:) 389 kptns0(:,:)=kptns(:,:) 390 npwarr0(:)=npwarr(:) 391 mpw0=mpw 392 393 do isppol=1,sppoldbl 394 do ikpt=1,nkpt 395 indkk(ikpt+(isppol-1)*nkpt,1)=ikpt 396 indkk(ikpt+(isppol-1)*nkpt,2:6)=0 397 end do 398 end do 399 dksqmax=0.0_dp 400 401 ! The treatment of nband0 asks for some care 402 if(ireadwf==0)then 403 nband0(:)=0 404 else 405 nband0(1:nkpt0*nsppol0)=hdr0%nband(1:nkpt0*nsppol0) 406 end if 407 408 mpi_enreg0 => mpi_enreg 409 410 else 411 mpi_enreg0 => mpi_enreg 412 end if 413 414 if(mpi_enreg0%paral_pert == 1.and.mpi_enreg0%me_pert/=-1) then 415 me0 = mpi_enreg0%me_pert 416 else 417 me0 = mpi_enreg0%me_cell 418 end if 419 420 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 421 !Before hdr_free: 422 !If restartpaw==1, store hdr0%pawrhoij in hdr%pawrhoij; else if restartpaw==0, 423 !hdr%pawrhoij(:)has been initialized in hdr_init. 424 if(restartpaw==1) then 425 call pawrhoij_copy(hdr0%pawrhoij,hdr%pawrhoij,keep_itypat=.true.) 426 end if 427 428 call timab(711,2,tsec) 429 call timab(712,1,tsec) 430 431 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 432 !At this stage, all the relevant information from the header of the disk file, 433 !has been exploited, and stored in variables, on all processors. 434 !It is also contained in hdr0 435 !(on all processors, except if restart=1 and localrdwf=0, 436 !in which case it is only on the master) 437 !These information might be changed later, while processing the 438 !wavefunction data, and converting it. The variable hdr0 might be kept 439 !for further checking, or reference, or debugging, but at present, 440 !it is simpler to close it. The other header, hdr, will be used for the new file, if any. 441 442 if(ask_accurate==1)then 443 444 ! Check whether the accuracy requirements might be fulfilled 445 if(ireadwf==0)then 446 write(message,'(9a)')& 447 & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,& 448 & 'present calculation. It was asked that the wavefunctions be accurate,',ch10,& 449 & 'but they were not even read.',ch10,& 450 & 'Action: use a wf file, with ireadwf/=0.' 451 ABI_ERROR(message) 452 end if 453 if(dksqmax>tol12)then 454 write(message, '(9a,es16.6,4a)' )& 455 & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,& 456 & 'present calculation. It was asked that the wavefunctions be accurate, but',ch10,& 457 & 'at least one of the k points could not be generated from a symmetrical one.',ch10,& 458 & 'dksqmax=',dksqmax,ch10,& 459 & 'Action: check your wf file and k point input variables',ch10,& 460 & ' (e.g. kptopt or shiftk might be wrong in the present dataset or the preparatory one.' 461 ABI_ERROR(message) 462 end if 463 if(dtset%nspinor/=nspinor0)then 464 write(message,'(a,a, a,a,a,a,a, a,a,2i5,a,a)')& 465 & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,& 466 & 'present calculation. It was asked that the wavefunctions be accurate, but',ch10,& 467 & 'nspinor differs in the file from the actual nspinor.',ch10,& 468 & 'nspinor,nspinor0=',dtset%nspinor,nspinor0,ch10,& 469 & 'Action: check your wf file, and nspinor input variables.' 470 ABI_ERROR(message) 471 end if 472 if((nsppol>nsppol0 .and. sppoldbl==1) .or. nsppol<nsppol0 ) then 473 write(message,'(a,a, a,a,a,a,a, a,a,3i5,a,a)')& 474 & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,& 475 & 'present calculation. It was asked that the wavefunctions be accurate, but',ch10,& 476 & 'the nsppol variables do not match in the file and in the actual calculation',ch10,& 477 & 'nsppol,nsppol,sppoldbl=',dtset%nspinor,nspinor0,sppoldbl,ch10,& 478 & 'Action: check your wf file, and nsppol input variables.' 479 ABI_ERROR(message) 480 end if 481 482 ! Now, check the number of bands 483 accurate=1 484 do isppol=1,nsppol 485 do ikpt=1,nkpt 486 ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1) 487 ikptsp =ikpt +(isppol-1)*nkpt 488 ikptsp0=ikpt0+(isppol-1)*(2-sppoldbl)*nkpt0 489 if(nband0(ikptsp0)<nband(ikptsp))accurate=0 490 end do 491 end do 492 if(accurate==0)then 493 write(message,'(a,a, a,a,a,a,a, a,a)')& 494 & 'The file ',trim(wff1%fname),' cannot be used to start the ',ch10,& 495 & 'present calculation. It was asked that the wavefunctions be accurate,',ch10,& 496 & 'but the number of bands differ in the file and in the actual calculation.',ch10,& 497 & 'Action: use a wf file with the correct characteristics.' 498 ABI_ERROR(message) 499 end if 500 501 end if 502 503 !Flag: do we need to translate WF to (from) spinors ? 504 nsppol2nspinor=0 505 if (nsppol0==2.and.dtset%nspinor==2) nsppol2nspinor=+1 506 if (nspinor0==2.and.nsppol==2) nsppol2nspinor=-1 507 508 !Take into account parallism over spinors 509 my_nspinor =max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 510 my_nspinor0=max(1,nspinor0/mpi_enreg0%nproc_spinor) 511 512 !Not all bands might be read, if not needed to fill the wavefunctions 513 mband0=maxval(nband0(1:nkpt0*nsppol0)) 514 mband0_rd=min(mband0,(mband/dtset%nspinor)*nspinor0) 515 516 !**************************************************************************** 517 !If needed, transfer the input wf from disk to core memory 518 !(in the parallel case, it allows to change localrdwf=0 in localrdwf=1) 519 520 mkmem0=0 521 522 if(xmpi_paral == 1 .or. mpi_enreg%paral_kgb == 1 .or. mpi_enreg%paral_hf == 1) then 523 if(localrdwf==0 .and. out_of_core)then 524 ABI_BUG('localrdwf==0 and mkmem==0 (out-of-core solution) are not allowed together (yet)') 525 end if 526 end if 527 528 call timab(712,2,tsec) 529 530 !Here, treat reading wavefunctions with mkmem/=0, first step 531 if(ireadwf==1 .and. (.not.out_of_core))then 532 533 call timab(713,1,tsec) 534 535 ! if(restart==1 .and. ireadwf==1 .and. mkmem/=0)then 536 537 ! Compute table of k point associations. Make a trial choice for nkassoc. 538 nkassoc=(nkpt/nkpt0+1)*2 539 ABI_MALLOC(indkk0,(nkpt0,nkassoc)) 540 ! Infinite loops are allowed in F90 541 do 542 indkk0(:,:)=0 543 increase_nkassoc=0 544 do ikpt=1,nkpt*sppoldbl 545 ikpt0=indkk(ikpt,1) 546 do ikassoc=1,nkassoc 547 if(indkk0(ikpt0,ikassoc)==0)then 548 indkk0(ikpt0,ikassoc)=ikpt 549 exit 550 end if 551 if(nkassoc==ikassoc)increase_nkassoc=1 552 end do 553 if(increase_nkassoc==1)then 554 ABI_FREE(indkk0) 555 nkassoc=2*nkassoc 556 ABI_MALLOC(indkk0,(nkpt0,nkassoc)) 557 exit 558 end if 559 end do 560 if(increase_nkassoc==0)exit 561 end do 562 563 ! DEBUG 564 ! write(std_out,*)' inwffil: indkk0, nkassoc=',nkassoc 565 ! do ikpt0=1,nkpt0 566 ! write(std_out,*)' ikpt0,indkk0(ikpt0,1)=',ikpt0,indkk0(ikpt0,1) 567 ! end do 568 ! ENDDEBUG 569 570 ! DEBUG 571 ! write(std_out,*)' inwffil : indkk(:,1)=',indkk(:,1) 572 ! write(std_out,*)' inwffil : sppoldbl=',sppoldbl 573 ! ENDDEBUG 574 575 ! To treat the case (nsppol0=2,nspinor0=1)<->(nsppol=1,nspinor=2), 576 ! apply the following trick: 577 ! 1- We call wfsinp with fake arguments (nsppol_eff and nspinor_eff) 578 ! 2- We transform collinear polarized WF into spinors 579 ! or spinors into collinear polarized WF 580 if (nsppol2nspinor/=0.and.out_of_core.and.dtset%usewvl==0) then 581 write(message, '(7a)')& 582 & 'When mkmem=0 (out-of-core), the wavefunction translator is unable',ch10,& 583 & 'to interchange spin-polarized wfs and spinor wfs.',ch10,& 584 & 'Action: use a non-spin-polarized wf to start a spinor wf,',ch10,& 585 & ' and a non-spinor wf to start a spin-polarized wf.' 586 ABI_ERROR(message) 587 end if 588 589 ! === Fake arguments definition for wfsinp 590 if (nsppol2nspinor==0.or.dtset%usewvl/=0) then 591 indkk_eff => indkk 592 nband_eff => nband 593 eigen_eff => eigen 594 cg_eff => cg 595 nspinor_eff=dtset%nspinor;nsppol_eff=nsppol;sppoldbl_eff=sppoldbl 596 mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff)) 597 else if (nsppol2nspinor==1.and.(.not.out_of_core)) then 598 nsppol_eff=2;nspinor_eff=1;sppoldbl_eff=1 599 ABI_MALLOC(indkk_eff,(nkpt*sppoldbl_eff,6)) 600 ABI_MALLOC(nband_eff,(nkpt*nsppol_eff)) 601 indkk_eff(1:nkpt,1:6) =indkk(1:nkpt,1:6) 602 nband_eff(1:nkpt) =nband(1:nkpt)/2 603 nband_eff(1+nkpt:2*nkpt)=nband(1:nkpt)/2 604 mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff)) 605 eigen_eff => eigen 606 cg_eff => cg 607 else if (nsppol2nspinor==-1.and.(.not.out_of_core)) then 608 ! WARNING: MT 07072011 -> this is memory consuming 609 ! A copy a spinorial WF (and eigenvalues) is temporary kept in memory; 610 ! But the case (nspinor=2 => nsppol=2) might be rare 611 ! and only useful for testing purposes. 612 ! => print a warning for the user 613 ! NOTE: in that case (nsppol=2), parallelization over spinors is not activated 614 615 write(message,'(5a)')& 616 & 'In the case of spinor WF read from disk and converted into',ch10,& 617 & 'spin-polarized non-spinor WF, the WF translator is memory',ch10,& 618 & 'consuming (a copy of the spinor WF is temporarily stored in memory).' 619 ABI_WARNING(message) 620 621 nsppol_eff=1;nspinor_eff=2;sppoldbl_eff=1 622 ABI_MALLOC(indkk_eff,(nkpt*sppoldbl_eff,6)) 623 ABI_MALLOC(nband_eff,(nkpt*nsppol_eff)) 624 indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6) 625 nband_eff(1:nkpt) =2*nband(1:nkpt) 626 mband_eff=maxval(nband_eff(1:nkpt*nsppol_eff)) 627 ABI_MALLOC(eigen_eff,((2*mband_eff)**formeig*mband_eff*nkpt*nsppol_eff)) 628 ABI_MALLOC(cg_eff,(2,mpw0*nspinor_eff*mband_eff*mkmem*nsppol_eff)) 629 end if 630 631 ! === nband0 argument definition for wfsinp 632 squeeze=0 633 ABI_MALLOC(cg_disk,(0,0)) 634 if(.not.out_of_core)then 635 ABI_MALLOC(nband0_rd,(nkpt0*nsppol0)) 636 nband0_rd(:)=0 637 do isppol=1,nsppol_eff 638 do ikpt=1,nkpt 639 ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1) 640 isppol0=min(isppol,nsppol0) 641 ikptsp =ikpt +(isppol -1)*nkpt 642 ikptsp0=ikpt0+(isppol0-1)*(2-sppoldbl)*nkpt0 643 nband0_k=min(nband0(ikptsp0),(nband_eff(ikptsp)/nspinor_eff)*nspinor0) 644 nband0_rd(ikptsp0)=max(nband0_rd(ikptsp0),nband0_k) 645 npw0=npwarr0(ikpt0) 646 npw =npwarr (ikpt) 647 if(npw0*nspinor0*nband0_k > npw*nspinor_eff*nband_eff(ikptsp))squeeze=1 648 end do 649 end do 650 if(squeeze==1)then 651 mcg_disk=mpw0*my_nspinor0*mband0_rd 652 ABI_FREE(cg_disk) 653 ABI_MALLOC(cg_disk,(2,mcg_disk)) 654 else 655 if(xmpi_paral == 1 .or. mpi_enreg0%paral_kgb == 1 .or. mpi_enreg0%paral_hf == 1)then 656 if(localrdwf==0)then 657 mcg_disk=mpw0*my_nspinor0*mband0_rd 658 ABI_FREE(cg_disk) 659 ABI_MALLOC(cg_disk,(2,mcg_disk)) 660 end if 661 end if 662 end if 663 end if 664 665 call timab(713,2,tsec) 666 call timab(714,1,tsec) 667 668 ! === call to wfsinp 669 if (dtset%usewvl == 0) then 670 call wfsinp(cg_eff,cg_disk,ecut,ecut0,ecut_eff,eigen,& 671 & exchn2n3d,formeig,gmet,gmet0,headform0,& 672 & indkk_eff,indkk0,istwfk,istwfk0,kptns,kptns0,localrdwf,& 673 & mband_eff,mcg,mcg_disk,mpi_enreg,mpi_enreg0,mpw,mpw0,& 674 & nband_eff,nband0_rd,ngfft,nkassoc,nkpt,nkpt0,npwarr,npwarr0,nspinor_eff,nspinor0,& 675 & nsppol_eff,nsppol0,nsym,occ,optorth,dtset%prtvol,randalg,restart,hdr%rprimd,sppoldbl_eff,squeeze,& 676 & symrel,tnons,wff1) 677 if (nsppol2nspinor/=0) then 678 ABI_FREE(indkk_eff) 679 ABI_FREE(nband_eff) 680 end if 681 else 682 ! Read wavefunctions from file. 683 call wvl_wfsinp_disk(dtset, hdr0, hdr, mpi_enreg, occ, 1, & 684 & hdr%rprimd, wff1, wvl%wfs, wvl%descr, hdr%xred) 685 end if 686 687 call timab(714,2,tsec) 688 call timab(715,1,tsec) 689 690 ! Now, update xyz0 variables, for use in newkpt 691 nband0(:)=nband0_rd(:) 692 693 ! If squeeze, the conversion was done in wfsinp, so no conversion left. 694 if(squeeze==1)then 695 ecut0=ecut_eff 696 gmet0(:,:)=gmet(:,:) 697 ABI_FREE(kptns0) 698 ABI_FREE(istwfk0) 699 ABI_FREE(nband0) 700 ABI_FREE(npwarr0) 701 ABI_MALLOC(kptns0,(3,nkpt)) 702 ABI_MALLOC(istwfk0,(nkpt)) 703 ABI_MALLOC(nband0,(nkpt*nsppol)) 704 ABI_MALLOC(npwarr0,(nkpt)) 705 kptns0(:,:)=kptns(:,:) 706 istwfk0(:)=istwfk(:) 707 npwarr0(:)=npwarr(:) 708 nband0(:)=0 709 do isppol=1,nsppol 710 do ikpt=1,nkpt 711 ikpt0=indkk(ikpt+(isppol-1)*(sppoldbl-1)*nkpt,1) 712 isppol0=min(isppol,nsppol0) 713 ikptsp =ikpt +(isppol -1)*nkpt 714 ikptsp0=ikpt0+(isppol0-1)*(sppoldbl-1)*nkpt0 715 nband0(ikptsp)=(nband0_rd(ikptsp0)/nspinor0)*dtset%nspinor 716 end do 717 end do 718 do ikpt=1,nkpt 719 indkk(ikpt,1)=ikpt 720 indkk(ikpt,2:6)=0 721 end do 722 ! This transfer must come after the nband0 transfer 723 nspinor0=dtset%nspinor 724 nkpt0=nkpt 725 nsppol0=nsppol 726 end if ! end squeeze == 1 727 728 ! The input wavefunctions have been transferred from disk to core memory 729 mkmem0=mkmem 730 731 ABI_FREE(indkk0) 732 ABI_FREE(nband0_rd) 733 ABI_FREE(cg_disk) 734 735 call timab(715,2,tsec) 736 737 else !ireadwf == 0 738 if (dtset%usewvl == 1) then 739 740 call timab(714,1,tsec) 741 ! Compute wavefunctions from input guess. 742 call wvl_wfsinp_scratch(dtset, mpi_enreg, occ, hdr%rprimd, wvl, hdr%xred) 743 call timab(714,2,tsec) 744 end if 745 end if 746 747 call timab(716,1,tsec) 748 749 !=== Eventual conversion of WF into (from) spinors 750 if (dtset%usewvl==0) then 751 752 ! ***** No conversion (standard case) **** 753 if (nsppol2nspinor==0) then 754 nspinor_eff=nspinor0;nsppol_eff=nsppol0;sppoldbl_eff=sppoldbl 755 indkk_eff => indkk 756 nband_eff => nband0 757 758 ! ***** Conversion from collinear to spinorial WF **** 759 else if (nsppol2nspinor==1.and.(.not.out_of_core)) then 760 ! Translate the WF and eigenvalues from nsppol=2 to nspinor=2 761 ! This is tricky (because we do not want to create a temporary array for cg) 762 nsppol_eff=1;nspinor_eff=2;sppoldbl_eff=1 763 ABI_MALLOC(indkk_eff,(nkpt*sppoldbl_eff,6)) 764 ABI_MALLOC(nband_eff,(nkpt0*nsppol_eff)) 765 indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6) 766 nband_eff(1:nkpt0)=2*nband0(1:nkpt0) 767 ! Compute some shifts from isspol0=1 to isppol0=2 768 imax=0;icgsft=0;ieigsft=0 769 ABI_MALLOC(my_kpt,(nkpt0)) 770 do ikpt0=1,nkpt0 771 nband0_k=nband0(ikpt0);nband_k=nband(ikpt0) 772 my_kpt(ikpt0)=(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me0))) 773 ieigsft=ieigsft+(2*nband0_k)**formeig*nband0_k 774 if(my_kpt(ikpt0)) then 775 imax=imax+nband0_k;icgsft=icgsft+nband0_k*npwarr0(ikpt0) 776 end if 777 end do 778 ! --- First version: no parallelization over spinors 779 if (mpi_enreg0%paral_spinor==0) then 780 ! Compute some useful indexes 781 ABI_MALLOC(indx,(2*imax)) 782 ABI_MALLOC(npwi,(imax)) 783 ii=0;icg=0 784 do ikpt0=1,nkpt0 785 if(my_kpt(ikpt0)) then 786 nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0) 787 do iband=1,nband0_k 788 ii=ii+1;npwi(ii)=npw0 789 indx(2*ii-1)=icg+mpw0;indx(2*ii)=icg+2*mpw0 790 icg=icg+4*mpw0 791 end do 792 end if 793 end do 794 ! Expand WF in cg (try to use the whole array) 795 ii=nsppol0*imax;icg0=nsppol0*icgsft 796 do isppol=nsppol0,1,-1 797 do ikpt0=nkpt0,1,-1 798 if(my_kpt(ikpt0)) then 799 nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0) 800 do iband=nband0_k,1,-1 801 icg0=icg0-npw0 802 if (indx(ii)<icg0) then 803 ABI_BUG("Unable to read WF!") 804 end if 805 cg(:,indx(ii)+1:indx(ii)+npw0)=cg(:,icg0+1:icg0+npw0) 806 ii=ii-1 807 end do 808 end if 809 end do 810 end do 811 ! Convert polarized WF into spinors 812 ii=1 813 do ikpt0=1,nkpt0 814 if(my_kpt(ikpt0)) then 815 nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0) 816 do iband=1,nband0_k 817 npw0=npwi(ii) 818 cg(:,indx(2*ii-1)-mpw0+1:indx(2*ii-1)-mpw0+npw0)=cg(:,indx(ii)+1:indx(ii)+npw0) 819 cg(:,indx(2*ii )+mpw0+1:indx(2*ii )+mpw0+npw0)=cg(:,indx(ii+imax)+1:indx(ii+imax)+npw0) 820 ii=ii+1 821 end do 822 end if 823 end do 824 ! Compress new cg array (from mpw to npw) and cancel zero-components 825 icg0=0;icg=0 826 do ikpt0=1,nkpt0 827 if(my_kpt(ikpt0)) then 828 nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0) 829 do iband=1,nband0_k 830 cg(:,icg0 +1:icg0+ npw0)=cg(:,icg+1:icg+npw0) 831 cg(:,icg0+ npw0+1:icg0+2*npw0)=zero 832 cg(:,icg0+2*npw0+1:icg0+3*npw0)=zero 833 cg(:,icg0+3*npw0+1:icg0+4*npw0)=cg(:,icg+3*mpw0+1:icg+3*mpw0+npw0) 834 icg0=icg0+4*npw0;icg=icg+4*mpw0 835 end do 836 end if 837 end do 838 ! --- Second version: parallelization over spinors 839 else 840 ! Compute some useful indexes 841 ABI_MALLOC(indx,(imax)) 842 ABI_MALLOC(npwi,(imax)) 843 ii=0;icg=0 844 do ikpt0=1,nkpt0 845 if(my_kpt(ikpt0)) then 846 nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0) 847 do iband=1,nband0_k 848 ii=ii+1;npwi(ii)=npw0 849 indx(ii)=icg+mpi_enreg0%me_spinor*mpw0 850 icg=icg+2*mpw0 851 end do 852 end if 853 end do 854 ! Expand WF in cg 855 ii=(mpi_enreg0%me_spinor+1)*imax;icg0=(mpi_enreg0%me_spinor+1)*icgsft 856 do ikpt0=nkpt0,1,-1 857 if(my_kpt(ikpt0)) then 858 nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0) 859 do iband=nband0_k,1,-1 860 icg0=icg0-npw0 861 if (indx(ii)<icg0) then 862 ABI_BUG("Unable to read WF!") 863 end if 864 cg(:,indx(ii)+1:indx(ii)+npw0)=cg(:,icg0+1:icg0+npw0) 865 ii=ii-1 866 end do 867 end if 868 end do 869 ! Compress new cg array (from mpw to npw) and cancel zero-components 870 icg0=0;icg=0 871 do ikpt0=1,nkpt0 872 if(my_kpt(ikpt0)) then 873 nband0_k=nband0(ikpt0);npw0=npwarr0(ikpt0) 874 do iband=1,nband0_k 875 if (mpi_enreg0%me_spinor==0) then 876 cg(:,icg0 +1:icg0+ npw0)=cg(:,icg+1:icg+npw0) 877 cg(:,icg0+npw0+1:icg0+2*npw0)=zero 878 else 879 cg(:,icg0 +1:icg0+ npw0)=zero 880 cg(:,icg0+npw0+1:icg0+2*npw0)=cg(:,icg+mpw0+1:icg+mpw0+npw0) 881 end if 882 icg0=icg0+2*npw0;icg=icg+2*mpw0 883 end do 884 end if 885 end do 886 end if 887 ! Translate eigenvalues 888 ibg0=2*ieigsft;ibg=2*ieigsft 889 do ikpt0=nkpt0,1,-1 890 nband0_k=nband0(ikpt0) 891 ibg0=ibg0- nband0_k*(2*nband0_k)**formeig 892 ibg =ibg -2*nband0_k*(2*nband0_k)**formeig 893 if(my_kpt(ikpt0)) then 894 do iband=nband0_k*(2*nband0_k)**formeig,1,-1 895 eigen(2*iband-1+ibg)=eigen(iband+ibg0-ieigsft) 896 eigen(2*iband +ibg)=eigen(iband+ibg0) 897 end do 898 end if 899 end do 900 ABI_FREE(indx) 901 ABI_FREE(npwi) 902 ABI_FREE(my_kpt) 903 904 ! ***** Conversion from spinorial to collinear WF **** 905 else if (nsppol2nspinor==-1.and.(.not.out_of_core)) then 906 ! In that case parallelization over spinors is never activated 907 nsppol_eff=2;nspinor_eff=1;sppoldbl_eff=1 908 ABI_MALLOC(indkk_eff,(nkpt*sppoldbl_eff,6)) 909 ABI_MALLOC(nband_eff,(nkpt0*nsppol_eff)) 910 indkk_eff(1:nkpt,1:6)=indkk(1:nkpt,1:6) 911 nband_eff(1:nkpt0) =nband0(1:nkpt0)/2 912 nband_eff(1+nkpt0:2*nkpt0)=nband0(1:nkpt0)/2 913 ! Compute shifts from isspol0=1 to isppol0=2 914 icgsft=0;ieigsft=0 915 do ikpt0=1,nkpt0 916 nband0_k=nband0(ikpt0);nband_k=nband(ikpt0) 917 ieigsft=ieigsft+(nband0_k/2)*(nband0_k)**formeig 918 if(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me))) & 919 & icgsft=icgsft+(nband0_k/2)*npwarr0(ikpt0) 920 end do 921 ! Translate the WF and eigenvalues from nspinor=2 to nsppol=2 922 icg0=0;icg=0;ibg=0 923 do ikpt0=1,nkpt0 924 nband0_k=nband0(ikpt0);nband_k=nband(ikpt0);npw0=npwarr0(ikpt0) 925 if(.not.(proc_distrb_cycle(mpi_enreg0%proc_distrb,ikpt0,1,nband_k,1,me))) then 926 do iband=1,nband0_k/2 927 do ii=1,npw0 928 cg(:,ii+icg) =cg_eff(:,ii+icg0) 929 cg(:,ii+icg+icgsft)=cg_eff(:,ii+icg0+3*npw0) 930 end do 931 icg0=icg0+4*npw0;icg=icg+npw0 932 end do 933 do iband=(nband0_k/2)*(nband0_k)**formeig,1,-1 934 eigen(iband+ibg) =eigen_eff(2*iband-1+2*ibg) 935 eigen(iband+ibg+ieigsft)=eigen_eff(2*iband +2*ibg) 936 ! occ(iband+ibg) =occ_eff(2*iband-1+2*ibg) 937 ! occ(iband+ibg+ieigsft)=occ_eff(2*iband +2*ibg) 938 end do 939 end if 940 ibg=ibg+(nband0_k/2)*(nband0_k)**formeig 941 end do 942 ABI_FREE(cg_eff) 943 ABI_FREE(eigen_eff) 944 945 else 946 ABI_BUG('unable to interchange nsppol and nspinor when mkmem=0') 947 end if 948 end if 949 950 !Clean hdr0 951 call hdr0%free() 952 953 call timab(716,2,tsec) 954 call timab(717,1,tsec) 955 956 957 !**************************************************************************** 958 !Now, treat translation of wavefunctions if wavefunctions are planewaves 959 960 ceksp=0; debug=0; doorth=1; fill=1 961 if (dtset%usewvl == 0) then 962 963 call newkpt(ceksp,cg,debug,ecut0,ecut,ecut_eff,eigen,exchn2n3d,& 964 & fill,formeig,gmet0,gmet,headform0,indkk_eff,& 965 & ab_out,ireadwf,istwfk0,istwfk,kg,kptns0,kptns,& 966 & mband,mcg,mkmem0,mkmem,mpi_enreg0,mpi_enreg,& 967 & mpw0,mpw,my_nkpt,nband_eff,nband,ngfft0,ngfft,nkpt0,nkpt,npwarr0,npwarr,& 968 & nspinor_eff,dtset%nspinor,nsppol_eff,nsppol,nsym,occ,optorth,& 969 & dtset%prtvol,randalg,restart,hdr%rprimd,sppoldbl_eff,symrel,tnons,unkg,wff1,wffnow) 970 971 if (nsppol2nspinor/=0) then 972 ABI_FREE(indkk_eff) 973 ABI_FREE(nband_eff) 974 end if 975 976 end if ! dtset%usewvl == 0 977 978 !**************************************************************************** 979 980 ABI_FREE(indkk) 981 ABI_FREE(istwfk0) 982 ABI_FREE(kptns0) 983 ABI_FREE(nband0) 984 ABI_FREE(npwarr0) 985 if (restart==2 .and.(mpi_enreg0%paral_kgb==1 .or. mpi_enreg0%paral_hf == 1)) then 986 call destroy_mpi_enreg(mpi_enreg0) 987 ABI_FREE(mpi_enreg0) 988 else 989 nullify(mpi_enreg0) 990 end if 991 992 call timab(717,2,tsec) 993 call timab(710,2,tsec) 994 995 DBG_EXIT("COLL") 996 997 end subroutine inwffil
m_inwffil/newkpt [ Functions ]
[ Top ] [ m_inwffil ] [ Functions ]
NAME
newkpt
FUNCTION
This subroutine writes a starting guess for wave function (set 2) It performs a "zero order" interpolation, ie simply searches the nearest available k-point. The data (set 1) associated with this point is either read from a disk file (with a random access reading routine), or input as argument.
INPUTS
ceksp2=if 1, center the sphere of pw on Gamma; if 0, on each k-point. doorth=1 to do orthogonalization debug=>0 for debugging output ecut1=kinetic energy cutoffs for basis sphere 1 (hartree) ecut2=kinetic energy cutoffs beyond which the coefficients of wf2 vanish (Ha) ecut2_eff=kinetic energy cut-off for basis sphere 2 (hartree) exchn2n3d=if 1, n2 and n3 are exchanged fill=if 1, fill the supplementary bands ; if 0, reduce the number of bands Note : must have fill/=0 in the parallel execution formeig=if 0, GS format for wfs, eig and occ ; if 1, RF format. gmet1(3,3), gmet2(3,3)=reciprocal space metrics (bohr^-2) headform1=header format (might be needed to read the block of wfs) indkk(nkpt2*sppoldbl,6)=describe k point number of kptns1 that allows to generate wavefunctions closest to given kpt2 (and possibly isppol2=2) indkk(:,1)=k point number of kpt1 indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a (if 0, means no symmetry operation, equivalent to identity ) indkk(:,3:5)=shift in reciprocal space to be given to kpt1a, to give kpt1b, that is the closest to ikpt2. indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise iout=unit number for output file ireadwf=if 0, no reading of disk wavefunction file (random or 0.0 initialisation) istwfk1(nkpt1)=input parameter that describes the storage of wfs in set1 istwfk2(nkpt2)=input parameter that describes the storage of wfs in set2 kg2(3,mpw2*mkmem2)=dimensionless coords of G vecs in basis sphere at k point kptns1(3,nkpt1), kptns2(3,nkpt2)=k point sets (reduced coordinates) mband2= maximum number of bands of the output wavefunctions mcg=dimension of the cg array In case mkmem2/=0, all the output data must find their place in cg, so that mcg must be at least Sum(ikpt,isppol) [npw*nspinor*nband](ikpt,isppol) where these data are related to the output parameters In case mkmem1/=0, the same is true, for the input parameters, however, the maximum number of bands that will be read will be at most (mband2/nspinor2)*nspinor1 In case mkmem1==0 and mkmem2==0, one must have at least mpw*nspinor*mband for BOTH the input and output parameters, taking into account the maximal number of band to be read, described above. In case mkmem1/=0 and mkmem2/=0, it is expected that the input cg array is organised using the output parameters nkpt2, nband2 ... This is needed, in order to use the same pointer. mkmem1= if 0, the input wf, eig, occ are available from disk mkmem2= if 0, the output wf, eig, occ must be written onto disk mpi_enreg1=information about MPI parallelization, for the input wf file mpi_enreg2=information about MPI parallelization, for the output wf file mpw1=maximum allowed number of planewaves at any k, for the input wf file mpw2=maximum allowed number of planewaves at any k, for the output wf file my_nkpt2= number of k points for the output wf file, handled by current processus nband1(nkpt1*nsppol1)=number of bands, at each k point, on disk nband2(nkpt2*nsppol2)=desired number of bands at each k point ngfft1(18)=all needed information about 3D FFT, for the input wf file ngfft2(18)=all needed information about 3D FFT, for the output wf file see ~abinit/doc/variables/vargs.htm#ngfft nkpt1, nkpt2=number of k points in each set npwarr1(nkpt1)=array holding npw for each k point (input wf file). npwarr2(nkpt2)=array holding npw for each k point (output wf file). nspinor1,nspinor2=number of spinorial components of the wavefunctions for each wf file (input or output) nsppol1=1 for unpolarized, 2 for spin-polarized, input wf file nsppol2=1 for unpolarized, 2 for spin-polarized, output wf file nsym=number of symmetry elements in space group optorth= 1 if the WFS have to be orthogonalized; 0 otherwise prtvol=control print volume and debugging randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility restart= if 2, conversion between wavefunctions if 1, direct restart is allowed (see hdr_check.f) rprimd(3,3)=dimensional primitive translations for real space (bohr) sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations tnons(3,nsym)=nonsymmorphic translations for symmetry operations unkg2=unit number for storage of basis sphere data: stores indirect indexing array and integer coordinates for all planewaves in basis sphere for each k point being considered (kptns2 set) wffinp=structure info of input wf file unit number wffout=structure info of output wf file unit number dtset <type(dataset_type)>=all input variables for this dataset
OUTPUT
(see side effects)
SIDE EFFECTS
The following arrays are input if mkmem1/=0, otherwise their input
values are taken from disk, and are output if mkmem2/=0, otherwise
their output values are written on disk.
The location of the block for a given spin-k point at input MUST
be the same as the location of the corresponding spin-k point at output.
cg(2,mcg)=complex wf array
eigen(mband2*(2*mband2)**formeig *nkpt2*nsppol2)=
eigenvalues (input or init to large number for GS or init to 0.0 for RF), (Ha)
occ(mband2*nkpt2*nsppol2)=occupation (input or init to 0.0) NOT USED NOW
NOTES
* When reading from disk, it is expected that the next record of the wffinp%unwff disk unit is the first record of the first wavefunction block. * When the data is input as argument, it is assumed that the data for each spin- k wavefunction block is located at the proper corresponding location of the output array (this is to be described). * The information is pumped onto an fft box for the conversion. This allows for changing the number of plane waves. * In the present status of this routine, occ is not output.
SOURCE
2042 subroutine newkpt(ceksp2,cg,debug,ecut1,ecut2,ecut2_eff,eigen,exchn2n3d,fill,& 2043 & formeig,gmet1,gmet2,headform1,indkk,iout,ireadwf,& 2044 & istwfk1,istwfk2,kg2,kptns1,kptns2,mband2,mcg,mkmem1,mkmem2,& 2045 & mpi_enreg1,mpi_enreg2,mpw1,mpw2,my_nkpt2,nband1,nband2,& 2046 & ngfft1,ngfft2,nkpt1,nkpt2,npwarr1,npwarr2,nspinor1,nspinor2,& 2047 & nsppol1,nsppol2,nsym,occ,optorth,prtvol,randalg,restart,rprimd,& 2048 & sppoldbl,symrel,tnons,unkg2,wffinp,wffout) 2049 2050 !Arguments ------------------------------------ 2051 !scalars 2052 integer,intent(in) :: ceksp2,debug,exchn2n3d,fill,formeig,headform1,iout 2053 integer,intent(in) :: ireadwf,mband2,mcg,mkmem1,mkmem2,mpw1,mpw2,my_nkpt2,nkpt1,nkpt2 2054 integer,intent(in) :: nspinor1,nspinor2,nsppol1,nsppol2,nsym,optorth,prtvol,restart 2055 integer,intent(in) :: randalg,sppoldbl,unkg2 2056 real(dp),intent(in) :: ecut1,ecut2,ecut2_eff 2057 type(MPI_type),intent(inout) :: mpi_enreg1,mpi_enreg2 2058 type(wffile_type),intent(inout) :: wffinp,wffout 2059 !arrays 2060 integer,intent(in) :: indkk(nkpt2*sppoldbl,6),istwfk1(nkpt1),istwfk2(nkpt2) 2061 integer,intent(in) :: kg2(3,mpw2*mkmem2),nband1(nkpt1*nsppol1) 2062 integer,intent(in) :: nband2(nkpt2*nsppol2),ngfft1(18),ngfft2(18) 2063 integer,intent(in) :: npwarr1(nkpt1),npwarr2(nkpt2),symrel(3,3,nsym) 2064 real(dp),intent(in) :: gmet1(3,3),gmet2(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2) 2065 real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym) 2066 real(dp),intent(inout) :: cg(2,mcg) !vz_i pw_orthon vecnm 2067 real(dp),intent(inout) :: eigen(mband2*(2*mband2)**formeig*nkpt2*nsppol2)!vz_i newocc 2068 real(dp),intent(inout) :: occ(mband2*nkpt2*nsppol2) !vz_i 2069 2070 !Local variables------------------------------- 2071 !scalars 2072 integer,parameter :: init_random=-5,nkpt_max=50,tobox=1,tosph=-1,wr=2 2073 integer :: aux_stor,band_index,iband,icg,icg_aux,idum 2074 integer :: ii,ikg2,ikpt1,ikpt10,ikpt2,ikptsp_prev,inplace,iproc 2075 integer :: isppol1,isppol2,istwf10_k,localrdwf 2076 integer :: mband1,mband_rd,mband_rw,mcg_aux,me1,me2,mgfft1,mgfft2 2077 integer :: my_nspinor1,my_nspinor2 2078 integer :: nb_band,nbd1,nbd1_rd,nbd2,nkpt_eff,nproc2,npw1,npw2,nsp 2079 integer :: test_cycle,tim_rwwf 2080 logical :: out_of_core2 2081 character(len=500) :: message 2082 !arrays 2083 integer,allocatable :: kg1(:,:),kg2_k(:,:),kg_dum(:,:) 2084 real(dp) :: kpoint(3),tsec(2) 2085 real(dp),allocatable :: cg_aux(:,:),eig_k(:),occ_k(:) 2086 2087 ! ************************************************************************* 2088 2089 call timab(780,1,tsec) 2090 call timab(781,1,tsec) 2091 2092 icg=0 2093 2094 !Init MPI data 2095 me1=mpi_enreg1%me_kpt 2096 me2=mpi_enreg2%me_kpt 2097 nproc2 = mpi_enreg2%nproc_cell 2098 out_of_core2=(my_nkpt2/=0.and.mkmem2==0) 2099 2100 2101 if((nsppol1==2.and.nspinor2==2).or.(nspinor1==2.and. nsppol2==2))then 2102 ! This is not yet possible. See later for a message about where to make the needed modifs. 2103 ! EDIT MT 20110707: these modifs are no more needed as they are now done in inwffil 2104 write(message, '(5a,i2,a,i2,2a,i2,a,i2,4a)' ) & 2105 & 'The wavefunction translator is (still) unable to interchange',ch10,& 2106 & 'spin-polarized wfs and spinor wfs. However,',ch10,& 2107 & 'the input variables are nsppol1=',nsppol1,', and nspinor1=',nspinor1,ch10,& 2108 & 'the output variables are nsppol2=',nsppol2,', and nspinor2=',nspinor2,ch10,& 2109 & 'Action: use a non-spin-polarized wf to start a spinor wf,',ch10,& 2110 & ' and a non-spinor wf to start a spin-polarized wf.' 2111 ABI_ERROR(message) 2112 end if 2113 2114 my_nspinor1=max(1,nspinor1/mpi_enreg1%nproc_spinor) 2115 my_nspinor2=max(1,nspinor2/mpi_enreg2%nproc_spinor) 2116 mband1=maxval(nband1(1:nkpt1*nsppol1)) 2117 2118 if(mkmem1==0 .and. out_of_core2)then 2119 mband_rd=min(mband1,(mband2/nspinor2)*nspinor1) 2120 if(mcg<mpw1*my_nspinor1*mband_rd)then 2121 write(message,'(2(a,i0))')' The dimension mcg= ',mcg,', should be larger than mband_rd= ',mband_rd 2122 ABI_BUG(message) 2123 end if 2124 if(mcg<mband2*mpw2*my_nspinor2)then 2125 write(message,'(a,i0,a,a,a,i0,a,i0,a,i2)' )& 2126 & 'The dimension mcg= ',mcg,', should be larger than',ch10,& 2127 & 'the product of mband2= ',mband2,', mpw2= ',mpw2,', and nspinor2= ',my_nspinor2 2128 ABI_BUG(message) 2129 end if 2130 end if 2131 2132 idum=init_random 2133 ikpt10 = 0 2134 istwf10_k=0 2135 band_index=0 2136 icg=0 2137 2138 nkpt_eff=nkpt2 2139 if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max ) nkpt_eff=nkpt_max 2140 2141 mgfft1=maxval(ngfft1(1:3)) 2142 mgfft2=maxval(ngfft2(1:3)) 2143 ABI_MALLOC(kg1,(3,mpw1)) 2144 ABI_MALLOC(kg2_k,(3,mpw2)) 2145 ABI_MALLOC(kg_dum,(3,0)) 2146 2147 if (debug>0) then 2148 if (me1==0) then 2149 write(std_out,'(a)' ) ' newkpt: kptns1' 2150 call prmat (kptns1, 3, nkpt1, 3) 2151 end if 2152 if (me2==0) then 2153 write(std_out,'(a)' ) ' newkpt: kptns2' 2154 call prmat (kptns2, 3, nkpt2, 3) 2155 end if 2156 end if 2157 2158 ikptsp_prev=0 2159 2160 call timab(781,2,tsec) 2161 2162 !Do outer loop over spins 2163 do isppol2=1,nsppol2 2164 2165 if (nsppol2==2 .and. me2==0) then 2166 write(std_out,'(a,i5)' ) ' newkpt: spin channel isppol2 = ',isppol2 2167 end if 2168 2169 if (restart==1 .and. out_of_core2) rewind (unkg2) 2170 ikg2=0 2171 2172 ! Do loop over new k point set 2173 do ikpt2=1,nkpt2 2174 2175 call timab(782,1,tsec) 2176 2177 nbd2=nband2(ikpt2+(isppol2-1)*nkpt2) 2178 npw2=npwarr2(ikpt2) 2179 2180 if(restart==1)then 2181 2182 ! Announce the treatment of k point ikpt 2183 if(ikpt2<=nkpt_eff)then 2184 ! This message might be overwritten in parallel 2185 write(message, '(a,i6,a,i8,a,i4)' )'P newkpt: treating ',nbd2,' bands with npw=',npw2,' for ikpt=',ikpt2 2186 ! This message might be overwritten in parallel 2187 if(mpi_enreg2%paralbd==1)then 2188 do iproc=0,nproc2-1 2189 nb_band=0 2190 do iband=1,nbd2 2191 if(mpi_enreg2%proc_distrb(ikpt2,iband,isppol2) == iproc)nb_band=nb_band+1 2192 end do 2193 if(nb_band/=0)then 2194 write(message, '(a,i6,a,i8,a,i4,a,i4)' ) & 2195 & 'P newkpt: treating ',nb_band,' bands with npw=',npw2,' for ikpt=',ikpt2,' by node ',iproc 2196 end if 2197 end do 2198 end if 2199 if(mpi_enreg2%paralbd==0) then 2200 write(message, '(a,i6,a,i8,a,i4,a,i4)' )& 2201 & 'P newkpt: treating ',nbd2,' bands with npw=',npw2,& 2202 & ' for ikpt=',ikpt2,' by node ',mpi_enreg2%proc_distrb(ikpt2,1,isppol2) 2203 end if 2204 if(prtvol>0)then 2205 call wrtout(iout,message,'COLL') 2206 end if 2207 end if 2208 2209 ! Cut the writing if the limit is reached 2210 if(ikpt2==nkpt_eff+1)then 2211 if(prtvol>0)then 2212 call wrtout(iout,' newkpt: prtvol=0 or 1, do not print more k-points.','COLL') 2213 end if 2214 end if 2215 2216 ! End of restart==1 2217 end if 2218 2219 test_cycle=0 2220 if(proc_distrb_cycle(mpi_enreg2%proc_distrb,ikpt2,1,nbd2,isppol2,me2)) test_cycle=1 2221 if(test_cycle==1)then 2222 if(formeig==0)then 2223 eigen(1+band_index : nbd2+band_index) = zero 2224 ! occ(1+band_index : nbd2+band_index) = zero 2225 band_index=band_index+nbd2 2226 else 2227 eigen(1+band_index : 2*nbd2**2+band_index) = 0.0_dp 2228 band_index=band_index+2*nbd2**2 2229 end if 2230 ! In the case this k point does not belong to me, cycle 2231 if (my_nkpt2==0) cycle 2232 if ((mkmem1==0) .and. (ireadwf==1) .and. (mpi_enreg2%paralbd==1))then 2233 call WffReadSkipK(formeig,headform1,ikpt2,isppol2,mpi_enreg2,wffinp) 2234 ikptsp_prev=ikptsp_prev+1 2235 end if 2236 cycle 2237 end if 2238 2239 if(restart==1)then 2240 2241 if(mkmem2/=0)then 2242 kg2_k(:,1:npw2)=kg2(:,1+ikg2:npw2+ikg2) 2243 else if(mkmem2==0)then 2244 ! Read the first line of a block and performs some checks on the unkg file. 2245 ABI_ERROR("mkmem2 == 0 and rdnpw are not supported anymore.") 2246 nsp=nspinor2 2247 !call rdnpw(ikpt2,isppol2,nbd2,npw2,nsp,0,unkg2) 2248 ! Read k+g data 2249 read (unkg2) kg2_k(1:3,1:npw2) 2250 end if 2251 2252 end if 2253 2254 ! Get ikpt1, the closest k from original set, from indkk 2255 ikpt1=indkk(ikpt2,1) 2256 if(sppoldbl==2 .and. isppol2==2)ikpt1=indkk(ikpt2+nkpt2,1) 2257 2258 npw1=npwarr1(ikpt1) 2259 kpoint(:)=kptns1(:,ikpt1) 2260 2261 ! Determine the spin polarization of the input data 2262 isppol1=isppol2 2263 if(nsppol2==2 .and. nsppol1==1)isppol1=1 2264 2265 if(restart==2)then 2266 if(ikpt2<=nkpt_eff)then 2267 write(message,'(a,i4,i8,a,i4,i8)')'- newkpt: read input wf with ikpt,npw=',ikpt1,npw1,', make ikpt,npw=',ikpt2,npw2 2268 call wrtout(std_out,message) 2269 if(iout/=6 .and. me2==0 .and. prtvol>0)then 2270 call wrtout(iout,message) 2271 end if 2272 else if(ikpt2==nkpt_eff+1)then 2273 call wrtout(std_out, '- newkpt: prtvol=0 or 1, do not print more k-points.') 2274 if(iout/=6 .and. me2==0 .and. prtvol>0)then 2275 call wrtout(iout,message) 2276 end if 2277 end if 2278 end if 2279 2280 ! Set up the number of bands to be read 2281 nbd1=nband1(ikpt1+(isppol1-1)*nkpt1) 2282 nbd1_rd=min(nbd1,(nbd2/nspinor2)*nspinor1) 2283 2284 ! Check that number of bands is not being increased if fill==0 --if so 2285 ! print warning and reset new wf file nband2 to only allowed number 2286 if ( nbd2/nspinor2 > nbd1/nspinor1 .and. fill==0) then 2287 if(ikpt2<=nkpt_eff)then 2288 write(message, '(a,i8,a,i8,a,i8)' )' newkpt: nband2=',nbd2,' < nband1=',nbd1,' => reset nband2 to ',nbd1 2289 call wrtout(std_out,message) 2290 end if 2291 nbd2=nbd1 2292 end if 2293 2294 ! Prepare the reading of the wavefunctions: the correct record is selected 2295 ! WARNING : works only for GS - for RF the number of record differs 2296 if(restart==2 .and. mkmem1==0)then 2297 ABI_ERROR("mkmem1 == 0 has been removed.") 2298 2299 if(debug>0)then 2300 write(message, '(a,a,a,a,i5,a,i5,a,a,i5,a,i5)' ) ch10,& 2301 ' newkpt: about to call randac',ch10,& 2302 ' for ikpt1=',ikpt1,', ikpt2=',ikpt2,ch10,& 2303 ' and isppol1=',isppol1,', isppol2=',isppol2 2304 call wrtout(std_out,message) 2305 end if 2306 2307 !call randac(debug,headform1,ikptsp_prev,ikpt1,isppol1,nband1,nkpt1,nsppol1,wffinp) 2308 end if 2309 2310 ! Read the data for nbd2 bands at this k point 2311 ! Must decide whether an auxiliary storage is needed 2312 ! When mkmem1==0 and mkmem2==0 , the cg array should be large enough ... 2313 ! When mkmem1==0 and mkmem2/=0 , each k-point block in cg might not be large enough 2314 ! however, will read at most (nbd2/nspinor2)*nspinor1 bands from disk 2315 ! When mkmem1/=0 , it is supposed that each input k-point block is smaller 2316 ! than the corresponding output k-point block, so that the input data 2317 ! have been placed already in cg, at the k-point location where they are needed 2318 aux_stor=0 2319 if(mkmem2/=0 .and. mkmem1==0)then 2320 mcg_aux=npw1*my_nspinor1*nbd1 2321 if(nbd1_rd<nbd1)mcg_aux=npw1*my_nspinor1*nbd1_rd 2322 if( mcg_aux > npw2*my_nspinor2*nbd2 )then 2323 aux_stor=1 ; icg_aux=0 2324 ABI_MALLOC(cg_aux,(2,mcg_aux)) 2325 end if 2326 end if 2327 2328 mband_rw=max(nbd1_rd,nbd2) 2329 ABI_MALLOC(eig_k,(mband_rw*(2*mband_rw)**formeig)) 2330 if(formeig==0) then 2331 ABI_MALLOC(occ_k,(mband_rw)) 2332 else 2333 ABI_MALLOC(occ_k,(0)) 2334 end if 2335 2336 if(mkmem1/=0 .and. ireadwf==1)then 2337 ! Checks that nbd1 and nbd1_rd are equal if eig and occ are input 2338 if(nbd1/=nbd1_rd)then 2339 write(message,'(a,a,a,i6,a,i6)')& 2340 & 'When mkmem1/=0, one must have nbd1=nbd1_rd, while',ch10,& 2341 & 'nbd1 = ',nbd1,', and nbd1_rd = ',nbd1_rd 2342 ABI_BUG(message) 2343 end if 2344 ! Need to put eigenvalues in eig_k, same for occ 2345 ! Note use of band_index, since it is assumed that eigen and occ 2346 ! already have spin-k point location structure than output. 2347 if(formeig==0)then 2348 eig_k(1:nbd1_rd)=eigen(1+band_index : nbd1_rd+band_index) 2349 ! occ_k(1:nbd1_rd)=occ(1+band_index : nbd1_rd+band_index) 2350 else if(formeig==1)then 2351 ! The matrix of eigenvalues has size nbd1 , that must be equal 2352 ! to nbd1_rd in the case mkmem1/=0) 2353 eig_k(1:2*nbd1_rd**2)=eigen(1+band_index : 2*nbd1_rd**2+band_index) 2354 end if 2355 end if 2356 2357 call timab(782,2,tsec) 2358 2359 ! Must read the wavefunctions if they are not yet in place 2360 if(mkmem1==0 .and. ireadwf==1)then 2361 2362 if (debug>0 .and. restart==2) then 2363 write(message,'(a,i5,a,a,i5,a,i5,a)' ) & 2364 & ' newkpt: about to call rwwf with ikpt1=',ikpt1,ch10,& 2365 & ' and nband(ikpt1)=',nband1(ikpt1),' nbd2=',nbd2,'.' 2366 call wrtout(std_out,message) 2367 end if 2368 2369 if(mpi_enreg1%paralbd==0)tim_rwwf=21 2370 if(mpi_enreg1%paralbd==1)tim_rwwf=22 2371 2372 if(aux_stor==0)then 2373 call rwwf(cg,eig_k,formeig,headform1,icg,ikpt1,isppol1,kg_dum,mband_rw,mcg,mpi_enreg1,& 2374 & nbd1_rd,nbd1,npw1,my_nspinor1,occ_k,1,0,tim_rwwf,wffinp) 2375 else 2376 icg_aux=0 2377 call rwwf(cg_aux,eig_k,formeig,headform1,icg_aux,ikpt1,isppol1,kg_dum,mband_rw,mcg_aux,& 2378 & mpi_enreg1,nbd1_rd,nbd1,npw1,my_nspinor1,occ_k,1,0,tim_rwwf,wffinp) 2379 end if 2380 end if 2381 2382 call timab(783,1,tsec) 2383 2384 if(formeig==1 .and. nbd2/=nbd1_rd .and. ireadwf==1)then 2385 ! Change the storage of eig_k 2386 if(nbd1_rd<nbd2)then 2387 do iband=nbd1_rd,1,-1 2388 ! The factor of two is for complex eigenvalues 2389 do ii=2*nbd2,2*nbd1_rd+1,-1 2390 eig_k(ii+(iband-1)*2*nbd2)=huge(zero)/10.0_dp 2391 end do 2392 do ii=2*nbd1_rd,1,-1 2393 eig_k(ii+(iband-1)*2*nbd2)=eig_k(ii+(iband-1)*2*nbd1_rd) 2394 end do 2395 end do 2396 else if(nbd1_rd>nbd2)then 2397 do iband=1,nbd2 2398 ! The factor of two is for complex eigenvalues 2399 do ii=1,2*nbd2 2400 eig_k(ii+(iband-1)*2*nbd2)=eig_k(ii+(iband-1)*2*nbd1_rd) 2401 end do 2402 end do 2403 end if 2404 end if 2405 2406 ! If change nsppol, must adapt the occupation numbers 2407 ! if(nsppol1/=nsppol2)then 2408 ! occ_k(1:nbd2)=occ_k(1:nbd2)*nsppol1/dbl(nsppol2) 2409 ! then 2410 2411 ! In case nsppol1=2 and nspinor2=2, one should read 2412 ! the other spin-component, and form a spinor wf here, before calling 2413 ! wfconv. One should treat eig_k and occ_k as well. 2414 ! A similar operation is to be performed when nspino1=2 and nsppol2=2 2415 ! EDIT - MT 20110707: the building of the spinor wf is now done in wfffil 2416 ! no need to make it here.... 2417 2418 ! DEBUG 2419 ! write(std_out,*)' newkpt: before wfconv' 2420 ! write(std_out,*)' newkpt: mkmem2=',mkmem2 2421 ! stop 2422 ! ENDDEBUG 2423 2424 call timab(783,2,tsec) 2425 call timab(784,1,tsec) 2426 2427 ! Note the use of mband2, while mband is used inside 2428 ! write(std_out,*) 'in newkpt,before wfconv,npw1,npw2',npw1,npw2 2429 inplace=1 2430 if(aux_stor==0)then 2431 call wfconv(ceksp2,cg,cg,debug,ecut1,ecut2,ecut2_eff,& 2432 & eig_k,eig_k,exchn2n3d,formeig,gmet1,gmet2,icg,icg,& 2433 & ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,& 2434 & kg1,kg2_k,kptns1,kptns2,mband_rw,mband_rw,mcg,mcg,& 2435 & mpi_enreg1,mpi_enreg2,mpw1,mpw2,nbd1_rd,nbd2,& 2436 & ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,nsym,& 2437 & occ_k,occ_k,optorth,randalg,restart,rprimd,sppoldbl,symrel,tnons) 2438 else 2439 call wfconv(ceksp2,cg_aux,cg_aux,debug,ecut1,ecut2,ecut2_eff,& 2440 & eig_k,eig_k,exchn2n3d,formeig,gmet1,gmet2,icg_aux,icg_aux,& 2441 & ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,& 2442 & kg1,kg2_k,kptns1,kptns2,mband_rw,mband_rw,mcg,mcg,& 2443 & mpi_enreg1,mpi_enreg2,mpw1,mpw2,nbd1_rd,nbd2,& 2444 & ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,nsym,& 2445 & occ_k,occ_k,optorth,randalg,restart,rprimd,sppoldbl,symrel,tnons) 2446 end if 2447 2448 call timab(784,2,tsec) 2449 2450 ! Finally write new wf to disk file or save in permanent file 2451 if(mkmem2==0)then 2452 2453 ! Note that in this case, we are sure aux_stor==0 2454 if(mpi_enreg2%paralbd==0)tim_rwwf=21 2455 if(mpi_enreg2%paralbd==1)tim_rwwf=22 2456 call rwwf(cg,eig_k,formeig,0,0,ikpt2,isppol2,kg2_k,nbd2,mcg,mpi_enreg2,& 2457 & nbd2,nbd2,npw2,my_nspinor2,occ_k,wr,1,tim_rwwf,wffout) 2458 2459 end if 2460 2461 call timab(785,1,tsec) 2462 2463 if(mkmem2/=0)then 2464 if(aux_stor==1)then 2465 cg(:,1+icg:npw2*nbd2*my_nspinor2+icg)=cg_aux(:,1:npw2*nbd2*my_nspinor2) 2466 ABI_FREE(cg_aux) 2467 end if 2468 2469 icg=icg+npw2*nbd2*my_nspinor2 2470 ikg2=ikg2+npw2 2471 end if 2472 2473 eigen(1+band_index:nbd2*(2*nbd2)**formeig+band_index) = eig_k(1:nbd2*(2*nbd2)**formeig) 2474 ! occ(1+band_index:nbd2+band_index)=occ_k(1:nbd2) 2475 2476 if(formeig==0)then 2477 band_index=band_index+nbd2 2478 else if(formeig==1)then 2479 band_index=band_index+2*nbd2**2 2480 end if 2481 2482 ABI_FREE(eig_k) 2483 ABI_FREE(occ_k) 2484 2485 call timab(785,2,tsec) 2486 2487 end do ! ikpt2 2488 end do ! isppol2 2489 2490 call timab(786,1,tsec) 2491 2492 if(xmpi_paral==1)then 2493 ! Transmit eigenvalues (not yet occupation numbers) 2494 ! newkpt.F90 is not yet suited for RF format 2495 ! This routine works in both localrdwf=0 or 1 cases. 2496 ! However, in the present routine, localrdwf is to be considered 2497 ! as 1 always, since the transfer has been made in wfsinp . 2498 localrdwf=1 2499 call pareigocc(eigen,formeig,localrdwf,mpi_enreg2,mband2,nband2,nkpt2,nsppol2,occ,1) 2500 end if 2501 2502 ABI_FREE(kg1) 2503 ABI_FREE(kg2_k) 2504 ABI_FREE(kg_dum) 2505 2506 call timab(786,2,tsec) 2507 call timab(780,2,tsec) 2508 2509 end subroutine newkpt
m_inwffil/pareigocc [ Functions ]
[ Top ] [ m_inwffil ] [ Functions ]
NAME
pareigocc
FUNCTION
This subroutine transmit to all processors, using MPI:
- the eigenvalues and,
- if ground-state, the occupation numbers
(In fact, in the present status of the routine,
occupation numbers are NOT transmitted)
transmit_occ = 2 is used in case the occ should be transmitted.
Yet the code is not already written.
INPUTS
formeig=format of eigenvalues (0 for GS, 1 for RF)
localrdwf=(for parallel case) if 1, the eig and occ initial values
are local to each machine, if 0, they are on proc me=0.
mband=maximum number of bands of the output wavefunctions
mpi_enreg=information about MPI parallelization
nband(nkpt*nsppol)=desired number of bands at each k point
nkpt=number of k points
nsppol=1 for unpolarized, 2 for spin-polarized, output wf file processors,
Warning : defined only when paralbd=1
transmit_occ/=2 transmit only eigenvalues, =2 for transmission of occ also
(yet transmit_occ=2 is not safe or finished at all)
OUTPUT
(see side effects)
SIDE EFFECTS
eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha) occ(mband*nkpt*nsppol)=occupation (input or init to 0.0) NOT USED NOW
NOTES
* The case paralbd=1 with formeig=0 is implemented, but not yet used. * The transmission of occ is not activated yet ! * The routine takes the eigenvalues in the eigen array on one of the processors that possess the wavefunctions, and transmit it to all procs. If localrdwf==0, me=0 has the full array at start, If localrdwf==1, the transfer might be more complex. * This routine should not be used for RF wavefunctions, since it does not treat the eigenvalues as a matrix.
SOURCE
3274 subroutine pareigocc(eigen,formeig,localrdwf,mpi_enreg,mband,nband,nkpt,nsppol,occ,transmit_occ) 3275 3276 !Arguments ------------------------------------ 3277 !scalars 3278 integer,intent(in) :: formeig,localrdwf,mband,nkpt,nsppol,transmit_occ 3279 type(MPI_type),intent(in) :: mpi_enreg 3280 !arrays 3281 integer,intent(in) :: nband(nkpt*nsppol) 3282 real(dp),intent(inout) :: eigen(mband*(2*mband)**formeig*nkpt*nsppol) 3283 real(dp),intent(inout) :: occ(mband*nkpt*nsppol) 3284 3285 !Local variables------------------------------- 3286 !scalars 3287 integer :: band_index,iband,ierr,ikpt,isppol,me,nbks,spaceComm 3288 !character(len=500) :: msg 3289 !arrays 3290 real(dp) :: tsec(2) 3291 real(dp),allocatable :: buffer1(:),buffer2(:) 3292 3293 ! ************************************************************************* 3294 3295 if(xmpi_paral==1)then 3296 3297 ! Init mpi_comm 3298 spaceComm=mpi_enreg%comm_cell 3299 if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt 3300 if(mpi_enreg%paral_hf==1) spaceComm=mpi_enreg%comm_kpt 3301 ! Init me 3302 me=mpi_enreg%me_kpt 3303 3304 if(localrdwf==0)then 3305 call xmpi_bcast(eigen,0,spaceComm,ierr) 3306 3307 else if(localrdwf==1)then 3308 3309 ! Prepare transmission of eigen (and occ) 3310 ABI_MALLOC(buffer1,(2*mband**(formeig+1)*nkpt*nsppol)) 3311 ABI_MALLOC(buffer2,(2*mband**(formeig+1)*nkpt*nsppol)) 3312 buffer1(:)=zero 3313 buffer2(:)=zero 3314 3315 band_index=0 3316 do isppol=1,nsppol 3317 do ikpt=1,nkpt 3318 nbks=nband(ikpt+(isppol-1)*nkpt) 3319 3320 if(mpi_enreg%paralbd==0)then 3321 3322 if(formeig==0)then 3323 buffer1(2*band_index+1:2*band_index+nbks) = eigen(band_index+1:band_index+nbks) 3324 if(transmit_occ==2) then 3325 buffer1(2*band_index+nbks+1:2*band_index+2*nbks) = occ(band_index+1:band_index+nbks) 3326 end if 3327 band_index=band_index+nbks 3328 else if(formeig==1)then 3329 buffer1(band_index+1:band_index+2*nbks**2) = eigen(band_index+1:band_index+2*nbks**2) 3330 band_index=band_index+2*nbks**2 3331 end if 3332 3333 else if(mpi_enreg%paralbd==1)then 3334 3335 ! Skip this k-point if not the proper processor 3336 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nbks,isppol,me)) then 3337 if(formeig==0) then 3338 band_index=band_index+nbks 3339 else 3340 band_index=band_index+2*nbks**2 3341 end if 3342 cycle 3343 end if 3344 ! Loop on bands 3345 do iband=1,nbks 3346 if(mpi_enreg%proc_distrb(ikpt, iband,isppol) /= me)cycle 3347 if(formeig==0)then 3348 buffer1(2*band_index+iband)=eigen(band_index+iband) 3349 ! if(transmit_occ==2) buffer1(2*band_index+iband+nbdks)=occ(band_index+iband) 3350 else if (formeig==1)then 3351 buffer1(band_index+(iband-1)*2*nbks+1:band_index+(iband-1)*2*nbks+2*nbks) = & 3352 & eigen(band_index+(iband-1)*2*nbks+1:band_index+(iband-1)*2*nbks+2*nbks) 3353 end if 3354 end do 3355 if(formeig==0)then 3356 band_index=band_index+nbks 3357 else 3358 band_index=band_index+2*nbks**2 3359 end if 3360 end if 3361 3362 end do 3363 end do 3364 3365 ! Build sum of everything 3366 call timab(48,1,tsec) 3367 if(formeig==0)band_index=band_index*2 3368 call xmpi_sum(buffer1,buffer2,band_index,spaceComm,ierr) 3369 call timab(48,2,tsec) 3370 3371 band_index=0 3372 do isppol=1,nsppol 3373 do ikpt=1,nkpt 3374 nbks=nband(ikpt+(isppol-1)*nkpt) 3375 if(formeig==0)then 3376 eigen(band_index+1:band_index+nbks) = buffer2(2*band_index+1:2*band_index+nbks) 3377 if(transmit_occ==2) then 3378 occ(band_index+1:band_index+nbks) = buffer2(2*band_index+nbks+1:2*band_index+2*nbks) 3379 end if 3380 band_index=band_index+nbks 3381 else if(formeig==1)then 3382 eigen(band_index+1:band_index+2*nbks**2) = buffer1(band_index+1:band_index+2*nbks**2) 3383 band_index=band_index+2*nbks**2 3384 end if 3385 end do 3386 end do 3387 3388 ABI_FREE(buffer1) 3389 ABI_FREE(buffer2) 3390 end if 3391 end if 3392 3393 end subroutine pareigocc
m_inwffil/wfconv [ Functions ]
[ Top ] [ m_inwffil ] [ Functions ]
NAME
wfconv
FUNCTION
This subroutine treats the wavefunctions for one k point, and converts them to other parameters.
INPUTS
ceksp2=if 1, center the output sphere of pw on Gamma; if 0, on each k-point (usual). cg1(2,mcg1)=wavefunction array debug= if 1, print some messages ; otherwise, 0. ecut1=kinetic energy cutoffs for basis sphere 1 (hartree) ecut2=kinetic energy cutoff beyond which the coefficients of wf2 vanish (Ha) ecut2_eff=kinetic energy cut-off for basis sphere 2 (hartree) eig_k1(mband1*(2*mband1)**formeig)=eigenvalues exchn2n3d=if 1, n2 and n3 are exchanged formeig option (format of the eigenvalues and eigenvector) : 0 => ground-state format (initialisation of eigenvectors with random numbers, vector of eigenvalues) 1 => respfn format (initialisation of eigenvectors with 0 s, hermitian matrix of eigenvalues) gmet1(3,3)=reciprocal space metric (bohr^-2) for input wf gmet2(3,3)=reciprocal space metric (bohr^-2) for output wf icg1=shift to be given to the location of the data in the array cg1 icg2=shift to be given to the location of the data in the array cg2 ikpt1=number of the k point actually treated (input wf numbering) ikpt10=number of the k point previously treated (input wf numbering) ikpt2=number of the k point actually treated (output numbering) indkk(nkpt2*sppoldbl,6)=describe k point number of kptns1 that allows to generate wavefunctions closest to given kpt2 (and possibly isppol2=2) indkk(:,1)=k point number of kpt1 indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a (if 0, means no symmetry operation, equivalent to identity ) indkk(:,3:5)=shift in reciprocal space to be given to kpt1a, to give kpt1b, that is the closest to kpt2. indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise inplace= if 0, cg1 and cg2 are different in the calling routine, if 1, cg1 and cg2 are identical (they have the same memory location) This is also true for the pairs (eig_k1,eig_k2) and (occ_k1,occ_k2) isppol2=spin variable for output wavefunctions istwfk1(nkpt1)=input parameter that describes the storage of wfs in set1 istwfk2(nkpt2)=input parameter that describes the storage of wfs in set2 kg1(3,mpw1)=dimensionless coords of G vecs in basis sphere at k point (input wf) kg2(3,mpw2)=dimensionless coords of G vecs in basis sphere at k point (output wf) kptns1(3,nkpt1)=k point set for input wavefunctions kptns2(3,nkpt2)=k point set for output wavefunctions mband1=dimension of eig_k1 and occ_k1 arrays mband2=dimension of eig_k2 and occ_k2 arrays mcg1=dimension of cg1 array (at least npw1*nspinor1*nbd1) mcg2=dimension of cg2 array (at least npw2*nspinor2*nbd2) mpi_enreg1=information about MPI parallelization for set 1 mpi_enreg2=information about MPI parallelization for set 2 mpw1=dimension of kg1, can be set to 0 if not needed mpw2=dimension of kg2, can be set to 0 if not needed nbd1=number of bands contained in cg1,eig_k1,occ_k1 at this k-point - spin (at input) nbd2=number of bands contained in cg2,eig_k2,occ_k2 at this k-point - spin (at output) ngfft1(18)=all needed information about 3D FFT, for input wavefunctions ngfft2(18)=all needed information about 3D FFT, for output wavefunctions see ~abinit/doc/variables/vargs.htm#ngfft nkpt1=number of k points for input wavefunctions nkpt2=number of k points for output wavefunctions npw1=number of planewaves for input wavefunctions npw2=number of planewaves for output wavefunctions nspinor1=number of spinors for input wavefunctions nspinor2=number of spinors for output wavefunctions nsym=number of symmetry elements in space group occ_k1(mband1)=occupation numbers optorth=1 if the WFs are orthogonalized before leaving the routine randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility restart=if 2, conversion between wavefunctions if 1, direct restart is allowed (see hdr_check.f) rprimd2(3,3)=dimensional primitive translations for real space (bohr) needed only for the spinor rotation sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations tnons(3,nsym)=nonsymmorphic translations for symmetry operations dtset <type(dataset_type)>=all input variables for this dataset
OUTPUT
cg2(2,mcg2)=wavefunction array eig_k2(mband2*(2*mband2)**formeig)=eigenvalues occ_k2(mband2)=occupation (completed with zeros)
SIDE EFFECTS
Input/Output:
ikpt10=at input, number of the k point previously treated (input wf numbering)
(if this is the first call for the present k point set, ikpt10 should be 0)
at output, number of the k point just treated (input wf numbering)
kg1, kg2, npw1 and npw2 should not be modified by kpgsph (TD).
NOTES
Note that this routine can make an in-place conversion (see the input variable "inplace"), if cg1 and cg2 are equal, as well as the pairs (icg1,icg2), (eig_k1,eig_k2),(occ_k1,occ_k2) and (mband1,mband2) It can also be used to fill or to initialize wavefunctions at one k point (filling with random numbers or 0''s, according to the value of formeig), if the input number of bands (nbd1) is 0. In the latter case, one should use the same values of input wavefunction parameters than for output wavefunction parameters, except nbd1. The input parameters are indexed with 1, the output parameters are indexed with 2. Some of the arguments are arrays dimensioned with nkpt1 or nkpt2. Note that for these, only the elements for ikpt1 or ikpt2 will be used. The number of input bands must already be minimal at the input. This means, when input and output nspinor are equal : nbd1<nbd2 When the two nspinor differ, one must have nbd1/nspinor1<nbd2/nspinor2
SOURCE
2631 subroutine wfconv(ceksp2,cg1,cg2,debug,ecut1,ecut2,ecut2_eff,& 2632 & eig_k1,eig_k2,exchn2n3d,formeig,gmet1,gmet2,icg1,icg2,& 2633 & ikpt1,ikpt10,ikpt2,indkk,inplace,isppol2,istwfk1,istwfk2,& 2634 & kg1,kg2,kptns1,kptns2,mband1,mband2,mcg1,mcg2,mpi_enreg1,mpi_enreg2,& 2635 & mpw1,mpw2,nbd1,nbd2,ngfft1,ngfft2,nkpt1,nkpt2,npw1,npw2,nspinor1,nspinor2,& 2636 & nsym,occ_k1,occ_k2,optorth,randalg,restart,rprimd2,sppoldbl,symrel,tnons) 2637 2638 !Arguments ------------------------------------ 2639 !scalars 2640 integer,intent(in) :: ceksp2,debug,exchn2n3d,formeig,icg1,icg2,ikpt1 2641 integer,intent(in) :: ikpt2,inplace,isppol2,mband1,mband2,mcg1,mcg2,mpw1,mpw2 2642 integer,intent(in) :: nbd1,nbd2,nkpt1,nkpt2,nspinor1,nspinor2,nsym 2643 integer,intent(in) :: optorth,randalg,restart,sppoldbl 2644 integer,intent(inout) :: ikpt10,npw1,npw2 2645 real(dp),intent(in) :: ecut1,ecut2,ecut2_eff 2646 type(MPI_type),intent(inout) :: mpi_enreg1,mpi_enreg2 2647 !arrays 2648 integer,intent(in) :: indkk(nkpt2*sppoldbl,6),istwfk1(nkpt1),istwfk2(nkpt2) 2649 integer,intent(in) :: ngfft1(18),ngfft2(18),symrel(3,3,nsym) 2650 integer,intent(inout) :: kg1(3,mpw1),kg2(3,mpw2) 2651 real(dp),intent(in) :: gmet1(3,3),gmet2(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2) 2652 real(dp),intent(in) :: rprimd2(3,3),tnons(3,nsym) 2653 real(dp),intent(inout) :: cg1(2,mcg1),cg2(2,mcg2) 2654 real(dp),intent(inout) :: eig_k1(mband1*(2*mband1)**formeig) 2655 real(dp),intent(inout) :: eig_k2(mband2*(2*mband2)**formeig),occ_k1(mband1) 2656 real(dp),intent(inout) :: occ_k2(mband2) 2657 2658 !Local variables ------------------------------ 2659 !scalars 2660 integer,parameter :: nkpt_max=50,tobox=1,tosph=-1 2661 integer :: conv_tnons,convert,fftalg,fold1,fold2,foldim,foldre,i1,i2,iband 2662 integer :: iband_first,iband_last,icgmod,ierr,index,ipw 2663 integer :: ispinor,ispinor1,ispinor2,ispinor_first,ispinor_last 2664 integer :: istwf10_k,istwf1_k,istwf2_k,isym,itimrev 2665 integer :: mgfft1,mgfft2,n1,n2,n3,n4,n5,n6 2666 integer :: nbremn,npwtot,nspinor_index,nspinor1_this_proc,nspinor2_this_proc 2667 integer :: order,ortalgo 2668 real(dp) :: ai,ar,arg,bi,br,eig_tmp,spinrots,spinrotx,spinroty,spinrotz 2669 character(len=500) :: message 2670 integer, parameter :: int64 = selected_int_kind(18) 2671 integer(KIND=int64) :: seed 2672 !arrays 2673 integer :: atindx(1),identity(3,3),ngfft_now(18),no_shift(3),shiftg(3) 2674 integer :: symm(3,3),symrel_conv(3,3) 2675 integer,allocatable :: gbound1(:,:),gbound2(:,:) 2676 real(dp) :: kpoint1(3),kpoint2_sph(3),phktnons(2,1),spinrot(4),tnons_conv(3),tsec(2) 2677 real(dp),allocatable :: cfft(:,:,:,:),dum(:,:),phase1d(:,:),phase3d(:,:) 2678 real(dp),allocatable :: wavef1(:,:),wavef2(:,:),wavefspinor(:,:) 2679 2680 ! ************************************************************************* 2681 2682 mgfft1=maxval(ngfft1(1:3)) 2683 mgfft2=maxval(ngfft2(1:3)) 2684 if(.false.)write(std_out,*)occ_k1 ! just to keep occ_k1 as an argument before resolving the issue of its transfer 2685 2686 if(nspinor1/=1 .and. nspinor1/=2)then 2687 write(message,'(a,i0)')'The argument nspinor1 must be 1 or 2, while it is nspinor1 = ',nspinor1 2688 ABI_BUG(message) 2689 end if 2690 2691 if(nspinor2/=1 .and. nspinor2/=2)then 2692 write(message,'(a,i0)')' The argument nspinor2 must be 1 or 2, while it is nspinor2=',nspinor2 2693 ABI_BUG(message) 2694 end if 2695 2696 if(nspinor1==2 .and. mod(nbd1,2)/=0)then 2697 write(message,'(a,i0)')' When nspinor1 is 2, nbd1 must be even, while it is nbd1 = ',nbd1 2698 ABI_BUG(message) 2699 end if 2700 2701 if(nspinor2==2 .and. mod(nbd2,2)/=0)then 2702 write(message,'(a,i0)')' When nspinor2 is 2, nbd2 must be even, while it is nbd2=',nbd2 2703 ABI_BUG(message) 2704 end if 2705 2706 if(nbd1/nspinor1>nbd2/nspinor2)then 2707 write(message, '(3a,2i6,3a,2i6,a)' )& 2708 & 'In wfconv, the nbd/nspinor ratio cannot decrease. However,',ch10,& 2709 & 'the initial quantities are nbd1,nspinor1=',nbd1,nspinor1,', and',ch10,& 2710 & 'the requested final quantities are nbd2,nspinor2=',nbd2,nspinor2,'.' 2711 ABI_BUG(message) 2712 end if 2713 2714 ngfft_now(1:3)=ngfft1(1:3) 2715 ngfft_now(8:18)=ngfft1(8:18) 2716 !This line is the reason why ngfft_now has to be introduced 2717 ngfft_now(7)=101 2718 ngfft_now(4:6)=ngfft_now(1:3) 2719 n1=ngfft_now(1) ; n2=ngfft_now(2) ; n3=ngfft_now(3) 2720 n4=ngfft_now(4) ; n5=ngfft_now(5) ; n6=ngfft_now(6) 2721 fftalg=ngfft_now(7) 2722 2723 !Parallelization over spinors management 2724 nspinor1_this_proc=max(1,nspinor1/mpi_enreg1%nproc_spinor) 2725 nspinor2_this_proc=max(1,nspinor2/mpi_enreg2%nproc_spinor) 2726 2727 !In order to generate IN PLACE new wfs from old wfs, the loop 2728 !over bands and spinors must be done in one direction or the other, 2729 !depending on npw1 and npw2, nspinor1 and nspinor2. 2730 !If nspinor1=1 and nspinor2=2 , note that one will generate 2731 !from nbd1 states of npw1 coefficients, 2732 !2*nbd1 states of 2*npw2 coefficients. nbd1 cancels in comparing 2733 !these expressions, but nspinor2 appears squared. 2734 !The same line of thought works for the case nspinor1=2 and nspinor2=1 2735 order=1 2736 iband_first=1 ; iband_last=nbd1 2737 ispinor_first=1 ; ispinor_last=nspinor1 2738 if(nspinor1==2 .and. nspinor2==1)then 2739 order=2 ; iband_last=nbd1-1 ; ispinor_last=1 2740 end if 2741 !Here, reverse the order if needed 2742 if( npw2*nspinor2**2 > npw1*nspinor1**2 )then 2743 order=-order 2744 iband_first=iband_last ; iband_last=1 2745 ispinor_first=ispinor_last ; ispinor_last=1 2746 end if 2747 2748 kpoint1(:)=kptns1(:,ikpt1) 2749 istwf1_k=istwfk1(ikpt1) 2750 2751 kpoint2_sph(:)=0.0_dp 2752 if(ceksp2==0)kpoint2_sph(:)=kptns2(:,ikpt2) 2753 istwf2_k=istwfk2(ikpt2) 2754 2755 !DEBUG 2756 !write(std_out,*)'ecut1,ecut2_eff=',ecut1,ecut2_eff 2757 !write(std_out,*)'gmet1,gmet2=',gmet1,gmet2 2758 !write(std_out,*)'kpoint1,kpoint2_sph=',kpoint1,kpoint2_sph 2759 !write(std_out,*)'nspinor1,nspinor2',nspinor1,nspinor2 2760 !write(std_out,*)'istwf1_k,istwf2_k=',istwf1_k,istwf2_k 2761 !write(std_out,*)'nbd1,tol8=',nbd1,tol8 2762 !ENDDEBUG 2763 2764 !Determine whether it will be needed to convert the existing 2765 !wavefunctions, or simply to complete them. 2766 2767 convert=0 2768 if(nbd1/=0)then 2769 if(abs(ecut2_eff-ecut1)>tol8)convert=convert+1 2770 if(sum(abs(gmet2(:,:)-gmet1(:,:)))>tol8)convert=convert+2 2771 if(sum(abs(kpoint2_sph(:)-kpoint1(:)))>tol8)convert=convert+4 2772 if(nspinor2/=nspinor1)convert=convert+8 2773 if(istwf2_k/=istwf1_k)convert=convert+16 2774 end if 2775 2776 !This is a supplementary check 2777 if(restart==1 .and. convert/=0)then 2778 ABI_BUG('Restart==1 and convert/=0 are exclusive') 2779 end if 2780 2781 !Determine whether symmetries must be used 2782 conv_tnons=0 2783 no_shift(:)=0 2784 identity(:,:)=0 2785 identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1 2786 isym=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,2) 2787 !write(std_out,*)' wfconv : isym=',isym 2788 itimrev=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,6) 2789 if(isym/=0)then 2790 symrel_conv(:,:)=symrel(:,:,isym) 2791 call mati3inv(symrel_conv,symm) 2792 shiftg(:)=indkk(ikpt2+(sppoldbl-1)*(isppol2-1)*nkpt2,3:5) 2793 tnons_conv(:)=tnons(:,isym) 2794 if(sum(tnons_conv(:)**2)>tol8)then 2795 ! Need to compute phase factors associated with nonsymmorphic translations. 2796 conv_tnons=1 2797 ABI_MALLOC(phase3d,(2,npw1)) 2798 ABI_MALLOC(phase1d,(2,(2*n1+1)+(2*n2+1)+(2*n3+1))) 2799 ! Although the routine getph is originally written for 2800 ! atomic phase factors, it does precisely what we want 2801 atindx(1)=1 2802 call getph(atindx,1,n1,n2,n3,phase1d,tnons_conv) 2803 end if 2804 if(nspinor1==2 .and. nspinor2==2)then 2805 ! Compute rotation in spinor space 2806 call getspinrot(rprimd2,spinrot,symrel_conv) 2807 end if 2808 else 2809 shiftg(:)=0 2810 symm(:,:)=identity(:,:) 2811 spinrot(:)=zero 2812 spinrot(1)=one 2813 end if 2814 if(itimrev/=0)then 2815 symm(:,:)=-symm(:,:) 2816 end if 2817 2818 !DEBUG 2819 !write(std_out,'(a,i3,2x,3i3,2x,9i3)')' wfconv : isym,shiftg,symm=',isym,shiftg,symm 2820 !write(std_out,*)' wfconv : ecut2_eff,ecut1=',ecut2_eff,ecut1 2821 !write(std_out,*)' wfconv : istwf1_k,istwf2_k=',istwf1_k,istwf2_k 2822 !write(std_out,*)' wfconv : kpoint1(:),kpoint2_sph(:)=',& 2823 !& kpoint1(:),kpoint2_sph(:) 2824 !write(std_out,*)' wfconv : nspinor1,nspinor2=',nspinor1,nspinor2 2825 !ENDDEBUG 2826 2827 !if (mpi_enreg1%fft_option_lob==0) mpi_enreg1%fft_option_lob=1 2828 !if (mpi_enreg2%fft_option_lob==0) mpi_enreg2%fft_option_lob=1 2829 2830 if (restart==2.and.(convert/=0.or.(nbd2/nspinor2>nbd1/nspinor1.and.formeig==0))) then 2831 ! kg2 is needed both for FFT grid conversion and for envlop 2832 ! Choose the center of the sphere : either gamma, or each k-point 2833 kpoint2_sph(:)=0.0_dp 2834 if(ceksp2==0)kpoint2_sph(:)=kptns2(:,ikpt2) 2835 istwf2_k=istwfk2(ikpt2) 2836 call kpgsph(ecut2_eff,exchn2n3d,gmet2,0,ikpt2,istwf2_k,kg2,kpoint2_sph,1,mpi_enreg2,mpw2,npw2) 2837 end if 2838 2839 if(convert/=0)then 2840 istwf10_k=0 2841 if(ikpt10/=0)istwf10_k=istwfk1(ikpt10) 2842 2843 ! Only need G sphere if different from last time 2844 if ( ikpt1/=ikpt10 .or. istwf1_k/=istwf10_k ) then 2845 2846 call kpgsph (ecut1,exchn2n3d,gmet1,0,ikpt1,istwf1_k,kg1,kpoint1,1,mpi_enreg1,mpw1,npw1) 2847 if (debug>0) then 2848 write(message, '(a,f8.3,a,a,3f8.5,a,a,i3,a,3(a,3es16.8,a),a,3i4,a,i5,a)' )& 2849 & ' wfconv: called kpgsph with ecut1=',ecut1,ch10,& 2850 & ' kpt1=',kptns1(1:3,ikpt1),ch10,& 2851 & ' istwf1_k=',istwf1_k,ch10,& 2852 & ' gmet1= ',gmet1(1:3,1),ch10,& 2853 & ' ',gmet1(1:3,2),ch10,& 2854 & ' ',gmet1(1:3,3),ch10,& 2855 & ' ngfft=',ngfft_now(1:3),' giving npw1=',npw1,'.' 2856 call wrtout(std_out,message) 2857 end if 2858 ikpt10 = ikpt1 2859 istwf10_k=istwf1_k 2860 end if 2861 2862 if(conv_tnons==1)then 2863 arg=two_pi*(kpoint1(1)*tnons_conv(1)+ kpoint1(2)*tnons_conv(2)+ kpoint1(3)*tnons_conv(3) ) 2864 phktnons(1,1)=cos(arg) 2865 phktnons(2,1)=sin(arg) 2866 ! Convert 1D phase factors to 3D phase factors exp(i 2 pi (k+G).tnons ) 2867 call ph1d3d(1,1,kg1,1,1,npw1,n1,n2,n3,phktnons,phase1d,phase3d) 2868 end if 2869 2870 ABI_MALLOC(cfft,(2,n4,n5,n6)) 2871 ABI_MALLOC(wavef1,(2,npw1)) 2872 ABI_MALLOC(wavef2,(2,npw2)) 2873 if(nspinor1==2 .and. nspinor2==2) then 2874 ABI_MALLOC(wavefspinor,(2,2*npw2)) 2875 end if 2876 ABI_MALLOC(gbound1,(2*mgfft1+8,2)) 2877 ABI_MALLOC(gbound2,(2*mgfft2+8,2)) 2878 call sphereboundary(gbound1,istwf1_k,kg1,mgfft1,npw1) 2879 call sphereboundary(gbound2,istwf2_k,kg2,mgfft2,npw2) 2880 2881 ! Take old wf from sphere->box, the new from box->sphere 2882 ! One pays attention not to have a problem of erasing data when replacing 2883 ! a small set of coefficient by a large set, or the reverse. 2884 ! This is the reason of the use of order, _first and _last variables, 2885 ! defined earlier. 2886 nspinor_index=mpi_enreg1%me_spinor+1 2887 do iband=iband_first,iband_last,order 2888 do ispinor1=ispinor_first,ispinor_last,order 2889 ispinor=ispinor1 2890 if (mpi_enreg1%paral_spinor==1) then 2891 if (ispinor1==nspinor_index) then 2892 ispinor=1 2893 else 2894 if (nspinor1==2.and.nspinor2==2) wavefspinor(:,(ispinor1-1)*npw2+1:ispinor1*npw2)=zero 2895 cycle 2896 end if 2897 end if 2898 2899 ! Copy input wf 2900 i1=(ispinor-1)*npw1+(iband-1)*nspinor1_this_proc*npw1+icg1 2901 wavef1(:,1:npw1)=cg1(:,i1+1:i1+npw1) 2902 2903 ! Make symmetry-induced conversion, if needed (translation part) 2904 if(conv_tnons==1)then 2905 !$OMP PARALLEL DO PRIVATE(ai,ar) 2906 do ipw=1,npw1 2907 ar=phase3d(1,ipw)*wavef1(1,ipw)-phase3d(2,ipw)*wavef1(2,ipw) 2908 ai=phase3d(2,ipw)*wavef1(1,ipw)+phase3d(1,ipw)*wavef1(2,ipw) 2909 wavef1(1,ipw)=ar 2910 wavef1(2,ipw)=ai 2911 end do 2912 end if 2913 2914 ! Take into account time-reversal symmetry, if needed, in the scalar case 2915 if(itimrev==1 .and. (nspinor1==1 .or. nspinor2==1))then 2916 !$OMP PARALLEL DO 2917 do ipw=1,npw1 2918 wavef1(2,ipw)=-wavef1(2,ipw) 2919 end do 2920 end if 2921 2922 ! DEBUG 2923 ! write(std_out,*)' wfconv : before sphere, isym,ispinor=',isym,ispinor 2924 ! write(std_out,*)' no_shift,identity=',no_shift,identity 2925 ! write(std_out,*)' shiftg,symm=',shiftg,symm 2926 ! stop 2927 ! This debugging sequence is an attempt to rotate spinors, 2928 ! and works indeed for test13, when symmetry 9 is used ... 2929 ! if(isym==9 .and. ispinor==1)then 2930 ! write(std_out,*)' wfconv : gives a 120 degree rotation to first component' 2931 ! do ipw=1,npw1 2932 ! ar=- half*wavef1(1,ipw)-sqrt(three)*half*wavef1(2,ipw) 2933 ! ai= sqrt(three)*half*wavef1(1,ipw)- half*wavef1(2,ipw) 2934 ! wavef1(1,ipw)=ar 2935 ! wavef1(2,ipw)=ai 2936 ! end do 2937 ! end if 2938 ! ENDDEBUG 2939 2940 ! Convert wf, and also include the symmetry operation and shiftg. 2941 call sphere(wavef1,1,npw1,cfft,n1,n2,n3,n4,n5,n6,kg1,istwf1_k,tobox,& 2942 & mpi_enreg1%me_g0,no_shift,identity,one) 2943 2944 call sphere(wavef2,1,npw2,cfft,n1,n2,n3,n4,n5,n6,kg2,istwf2_k,tosph,& 2945 & mpi_enreg2%me_g0,shiftg,symm,one) 2946 2947 if(nspinor2==1 )then 2948 i2=(ispinor-1)*npw2+(iband-1)*nspinor2_this_proc*npw2+icg2 2949 cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2) 2950 else if(nspinor1==2.and.nspinor2==2)then 2951 ! Will treat this case outside of the ispinor loop 2952 i2=(ispinor1-1)*npw2 2953 wavefspinor(:,i2+1:i2+npw2)=wavef2(:,1:npw2) 2954 else if(nspinor1==1 .and. nspinor2==2)then 2955 ! The number of bands is doubled, and the number of coefficients 2956 ! is doubled also 2957 if (mpi_enreg2%paral_spinor==0) then 2958 i2=(iband-1)*nspinor2_this_proc*nspinor2_this_proc*npw2+icg2 2959 cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2) 2960 cg2(:,i2+npw2+1:i2+2*npw2)=zero 2961 cg2(:,i2+2*npw2+1:i2+3*npw2)=zero 2962 cg2(:,i2+3*npw2+1:i2+4*npw2)=wavef2(:,1:npw2) 2963 else 2964 i2=(iband-1)*nspinor2_this_proc*npw2+icg2 2965 if (nspinor_index==1) then 2966 cg2(:,i2+1:i2+npw2)=wavef2(:,1:npw2) 2967 cg2(:,i2+npw2+1:i2+2*npw2)=zero 2968 else 2969 cg2(:,i2+1:i2+npw2)=zero 2970 cg2(:,i2+npw2+1:i2+2*npw2)=wavef2(:,1:npw2) 2971 end if 2972 end if 2973 end if 2974 end do ! ispinor=ispinor_first,ispinor_last,order 2975 2976 if(nspinor1==2.and.nspinor2==2)then 2977 ! Take care of possible parallelization over spinors 2978 if (mpi_enreg2%paral_spinor==1) then 2979 call xmpi_sum(wavefspinor,mpi_enreg2%comm_spinor,ierr) 2980 end if 2981 ! Take care of time-reversal symmetry, if needed 2982 if(itimrev==1)then 2983 ! Exchange spin-up and spin-down 2984 ! Make complex conjugate of one component, 2985 ! and change sign of other component 2986 !$OMP PARALLEL DO PRIVATE(ipw,ar,ai) SHARED(wavefspinor,npw2) 2987 do ipw=1,npw2 2988 ! Here, change sign of real part 2989 ar=-wavefspinor(1,ipw) 2990 ai= wavefspinor(2,ipw) 2991 wavefspinor(1,ipw)= wavefspinor(1,npw2+ipw) 2992 ! Here, change sign of imaginary part 2993 wavefspinor(2,ipw)=-wavefspinor(2,npw2+ipw) 2994 wavefspinor(1,npw2+ipw)=ar 2995 wavefspinor(2,npw2+ipw)=ai 2996 end do 2997 end if ! itimrev==1 2998 2999 ! Rotation in spinor space 3000 !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(npw2,spinrot,wavefspinor) 3001 spinrots=spinrot(1) 3002 spinrotx=spinrot(2) 3003 spinroty=spinrot(3) 3004 spinrotz=spinrot(4) 3005 !$OMP DO 3006 do ipw=1,npw2 3007 ar=wavefspinor(1,ipw) 3008 ai=wavefspinor(2,ipw) 3009 br=wavefspinor(1,npw2+ipw) 3010 bi=wavefspinor(2,npw2+ipw) 3011 wavefspinor(1,ipw) = spinrots*ar - spinrotz*ai + spinroty*br - spinrotx*bi 3012 wavefspinor(2,ipw) = spinrots*ai + spinrotz*ar + spinroty*bi + spinrotx*br 3013 wavefspinor(1,npw2+ipw)= -spinroty*ar - spinrotx*ai + spinrots*br + spinrotz*bi 3014 wavefspinor(2,npw2+ipw)= -spinroty*ai + spinrotx*ar + spinrots*bi - spinrotz*br 3015 end do 3016 !$OMP END DO 3017 !$OMP END PARALLEL 3018 3019 ! Save wavefunction 3020 i2=(iband-1)*nspinor2_this_proc*npw2+icg2 3021 if (mpi_enreg2%paral_spinor==0) then 3022 cg2(:,i2 +1:i2+ npw2)=wavefspinor(:,1:npw2) 3023 cg2(:,i2+npw2+1:i2+2*npw2)=wavefspinor(:,npw2+1:2*npw2) 3024 else 3025 if (nspinor_index==1) then 3026 cg2(:,i2+1:i2+npw2)=wavefspinor(:,1:npw2) 3027 else 3028 cg2(:,i2+1:i2+npw2)=wavefspinor(:,npw2+1:2*npw2) 3029 end if 3030 end if 3031 end if ! nspinor1==2 .and. nspinor2==2 3032 3033 end do 3034 3035 ! Take care of copying eig and occ when nspinor increases or decreases 3036 if(nspinor1==1.and.nspinor2==2)then 3037 if(formeig==0)then 3038 ! Note the reverse order, needed in case inplace=1 3039 do iband=nbd1,1,-1 3040 ! use eig_tmp to avoid bug on ifort10.1 x86_64 3041 eig_tmp=eig_k1(iband) 3042 eig_k2(2*iband-1:2*iband)=eig_tmp 3043 ! occ_tmp=occ_k1(iband)*0.5_dp 3044 ! occ_k2(2*iband-1:2*iband )=occ_tmp 3045 end do 3046 else 3047 call wrtout(std_out,' wfconv: not yet coded, formeig=1!',"COLL") 3048 end if 3049 end if 3050 if(nspinor1==2 .and. nspinor2==1)then 3051 if(formeig==0)then 3052 do iband=1,nbd1 3053 ! use eig_tmp to avoid bug on ifort10.1 x86_64 3054 eig_tmp=eig_k1(2*iband-1) 3055 eig_k2(iband)=eig_tmp 3056 ! occ_tmp=occ_k1(2*iband-1)*2.0_dp 3057 ! occ_k2(iband)=occ_tmp 3058 end do 3059 else 3060 call wrtout(std_out,' wfconv: not yet coded, formeig=1!',"COLL") 3061 end if 3062 end if 3063 3064 ABI_FREE(cfft) 3065 ABI_FREE(gbound1) 3066 ABI_FREE(gbound2) 3067 ABI_FREE(wavef1) 3068 ABI_FREE(wavef2) 3069 if(nspinor1==2 .and. nspinor2==2) then 3070 ABI_FREE(wavefspinor) 3071 end if 3072 3073 else if(convert==0)then 3074 3075 if(inplace==0)then 3076 ! Must copy cg, eig and occ if not in-place while convert==0 3077 ! Note that npw1=npw2, nspinor1=nspinor2 3078 cg2(:,1+icg2:npw1*nspinor1_this_proc*nbd1+icg2)=& 3079 & cg1(:,1+icg1:npw1*nspinor1_this_proc*nbd1+icg1) 3080 eig_k2(:)=eig_k1(:) 3081 ! occ_k2(:)=occ_k1(:) 3082 end if 3083 3084 end if ! End of if convert/=0 3085 3086 if(conv_tnons==1) then 3087 ABI_FREE(phase1d) 3088 ABI_FREE(phase3d) 3089 end if 3090 3091 3092 !If not enough bands, complete with random numbers or zeros 3093 if(nbd2/nspinor2>nbd1/nspinor1)then 3094 if(formeig==0)then 3095 3096 ! Ground state wf and eig case 3097 eig_k2((nbd1/nspinor1)*nspinor2+1:nbd2)=huge(zero)/10.0_dp 3098 occ_k2((nbd1/nspinor1)*nspinor2+1:nbd2)=0.0_dp 3099 index=(nbd1/nspinor1)*nspinor2*npw2*nspinor2_this_proc 3100 3101 ! Initialisation of wavefunctions 3102 ! One needs to initialize wfs in such a way to avoid symmetry traps, 3103 ! and to avoid linear dependencies between wavefunctions 3104 ! No need for a difference for different k points and/or spin-polarization 3105 3106 npwtot=npw2 3107 if (mpi_enreg1%paral_kgb == 1) then 3108 call timab(539,1,tsec) 3109 call xmpi_sum(npwtot, mpi_enreg1%comm_bandfft, ierr) 3110 call timab(539,2,tsec) 3111 end if 3112 3113 do iband=(nbd1/nspinor1)*nspinor2+1,nbd2 3114 do ispinor2=1,nspinor2_this_proc 3115 ispinor=ispinor2;if (nspinor2_this_proc/=nspinor2) ispinor=mpi_enreg2%me_spinor+1 3116 3117 do ipw=1,npw2 3118 index=index+1 3119 ! Different seed for different planewave and band 3120 ! DEBUG seq==par 3121 ! if(.false.) then 3122 ! ENDDEBUG seq==par 3123 3124 if ( mpi_enreg2%paral_kgb /= 1.or.mpi_enreg2%nproc_cell == 1) then 3125 seed=(iband-1)*npw2*nspinor2 + (ispinor-1)*npw2 + ipw 3126 else 3127 seed=kg2(1,ipw)*npwtot*npwtot + kg2(2,ipw)*npwtot + kg2(3,ipw) 3128 seed=(iband*nspinor2+ispinor-1)*seed 3129 end if 3130 3131 if(randalg == 0) then 3132 ! For portability, use only integer numbers 3133 ! The series of couples (fold1,fold2) is periodic with a period of 3134 ! 3x5x7x11x13x17x19x23x29x31, that is, larger than 2**32, the largest integer*4 3135 ! fold1 is between 0 and 34, fold2 is between 0 and 114. As sums of five 3136 ! uniform random variables, their distribution is close to a gaussian 3137 fold1=modulo(seed,3)+modulo(seed,5)+modulo(seed,7)+modulo(seed,11)+modulo(seed,13) 3138 fold2=modulo(seed,17)+modulo(seed,19)+modulo(seed,23)+modulo(seed,29)+modulo(seed,31) 3139 3140 ! The gaussian distributions are folded, in order to be back to a uniform distribution 3141 ! foldre is between 0 and 20, foldim is between 0 and 18 3142 foldre=mod(fold1+fold2,21) 3143 foldim=mod(3*fold1+2*fold2,19) 3144 3145 cg2(1,index+icg2)=dble(foldre) 3146 cg2(2,index+icg2)=dble(foldim) 3147 else 3148 ! (Antoine Levitt) Simple linear congruential generator from 3149 ! numerical recipes, modulo'ed and 64bit'ed to avoid 3150 ! overflows (NAG doesn't like overflows, even though 3151 ! they are perfectly legitimate here). Then, we get some 3152 ! lowest order bits and sum them, as the previous 3153 ! generator, to get quasi-normal numbers. 3154 ! This is clearly suboptimal and might cause problems, 3155 ! but at least it doesn't seem to create linear 3156 ! dependencies and local minima like the previous one. 3157 ! it's not trivial to generate good reproductible random 3158 ! numbers in parallel. Patches welcome ! 3159 ! Note a fun fortran fact : MOD simply ignores 64 bits integer 3160 ! and casts them into 32bits, so we use MODULO. 3161 fold1 = modulo(1664525_int64 * seed + 1013904223_int64, 2147483648_int64) 3162 fold2 = modulo(1664525_int64 * fold1 + 1013904223_int64, 2147483648_int64) 3163 fold1=modulo(fold1,3)+modulo(fold1,5)+modulo(fold1,7)+modulo(fold1,11)+modulo(fold1,13) 3164 fold2=modulo(fold2,3)+modulo(fold2,5)+modulo(fold2,7)+modulo(fold2,11)+modulo(fold2,13) 3165 3166 cg2(1,index+icg2)=dble(fold1)/34-0.5 3167 cg2(2,index+icg2)=dble(fold2)/34-0.5 3168 end if 3169 end do 3170 end do 3171 3172 ! XG030513: Time-reversal symmetry for k=gamma imposes zero imaginary part at G=0 3173 ! XG: I do not know what happens for spin-orbit here. 3174 if (istwf2_k == 2 .and. mpi_enreg2%me_g0 == 1) then 3175 cg2(2,1+(iband-1)*npw2*nspinor2_this_proc+icg2)=zero 3176 end if 3177 end do ! iband 3178 3179 ! Multiply with envelope function to reduce kinetic energy 3180 icgmod=icg2+npw2*nspinor2_this_proc*(nbd1/nspinor1) 3181 nbremn=nbd2-nbd1 3182 call cg_envlop(cg2,ecut2,gmet2,icgmod,kg2,kpoint2_sph,mcg2,nbremn,npw2,nspinor2_this_proc) 3183 3184 if(ikpt2<=nkpt_max)then 3185 write(message,'(3(a,i6))')' wfconv:',nbremn,' bands initialized randomly with npw=',npw2,', for ikpt=',ikpt2 3186 call wrtout(std_out,message) 3187 end if 3188 3189 else if(formeig==1)then 3190 3191 ! For response function, put large numbers in the remaining of the 3192 ! eigenvalue array (part of it was already filled in calling routine) 3193 ! WARNING : Change of nspinor not yet coded 3194 eig_k2(1+2*nbd1*nbd2 : 2*nbd2*nbd2)=huge(zero)/10.0_dp 3195 ! Initialisation of wfs with 0 s 3196 index=npw2*nbd1*nspinor2_this_proc 3197 do iband=nbd1+1,nbd2 3198 do ipw=1,npw2*nspinor2_this_proc 3199 index=index+1 3200 cg2(:,index+icg2)=zero 3201 end do 3202 end do 3203 3204 if(ikpt2<=nkpt_max)then 3205 nbremn=nbd2-nbd1 3206 write(message,'(3(a,i0))')' wfconv:',nbremn,' bands set=0 with npw=',npw2,', for ikpt=',ikpt2 3207 call wrtout(std_out,message) 3208 end if 3209 3210 end if ! End of initialisation to 0 3211 end if 3212 3213 !Orthogonalize GS wfs 3214 if (optorth==1.and.formeig==0.and.(mpi_enreg2%paral_kgb/=1.or.mpi_enreg2%nproc_cell==1)) then 3215 !if (.True.) then 3216 ABI_MALLOC(dum,(2,0)) 3217 ortalgo=0 !;ortalgo=3 3218 call pw_orthon(icg2,0,istwf2_k,mcg2,0,npw2*nspinor2_this_proc,nbd2,ortalgo,dum,0,cg2,& 3219 & mpi_enreg2%me_g0,mpi_enreg2%comm_bandspinorfft) 3220 ABI_FREE(dum) 3221 end if 3222 3223 end subroutine wfconv
m_inwffil/wfsinp [ Functions ]
[ Top ] [ m_inwffil ] [ Functions ]
NAME
wfsinp
FUNCTION
Do initialization of wavefunction files.
Also call other relevant routines for this initialisation.
Detailed description :
- Initialize unit wff1 for input of wf data
formeig option (format of the eigenvalues and occupations) :
0 => ground-state format (initialisation of
eigenvectors with random numbers, vector of eigenvalues,
occupations are present)
1 => respfn format (initialisation of
eigenvectors with 0 s, hermitian matrix of eigenvalues)
INPUTS
ecut0=kinetic energy cutoffs for basis sphere 0 (hartree) (if squeeze=1) ecut=kinetic energy cutoffs beyond which the coefficients of cg vanish (Ha) (needed only if squeeze=1) ecut_eff=effective kinetic energy planewave cutoff (hartree), needed to generate the sphere of plane wave exchn2n3d=if 1, n2 and n3 are exchanged formeig=explained above gmet(3,3), gmet0(3,3)=reciprocal space metrics (bohr^-2) headform0=header format (might be needed to read the block of wfs) indkk(nkpt*sppoldbl,6)=describe k point number of kptns0 that allows to generate wavefunctions closest to given kpt indkk(:,1)=k point number of kptns0 indkk(:,2)=symmetry operation to be applied to kpt0, to give kpt0a (if 0, means no symmetry operation, equivalent to identity ) indkk(:,3:5)=shift in reciprocal space to be given to kpt0a, to give kpt0b, that is the closest to kpt. indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise indkk0(nkpt0,nkassoc)=list of k points that will be generated by k point number ikpt0 istwfk(nkpt)=input parameter that describes the storage of wfs istwfk0(nkpt0)=input parameter that describes the storage of wfs in set0 kptns(3,nkpt),kptns0(3,nkpt0)=k point sets (reduced coordinates) localrdwf=(for parallel case) if 1, the wff1%unwff file is local to each machine mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*nsppol mpi_enreg=information about MPI parallelization mpi_enreg0=information about MPI parallelization in set0 mpw=maximum number of planewaves as dimensioned in calling routine mpw0=maximum number of planewaves on disk file nban_dp_rd(nkpt0*nsppol0)=number of bands to be read at each k point nband(nkpt*nsppol)=number of bands at each k point ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkassoc=dimension of indkk0 array nkpt=number of k points expected nkpt0=number of k points on disk npwarr(nkpt)=array holding npw for each k point. npwarr0(nkpt0)=array holding npw for each k point, disk format. nspinor=number of spinorial components of the wavefunctions nspinor0=number of spinorial components of the wavefunctions on disk nsppol=1 for unpolarized, 2 for spin-polarized nsppol0=1 for unpolarized, 2 for spin-polarized, when on disk nsym=number of symmetry elements in space group optorth= 1 if the WFS have to be orthogonalized; 0 otherwise prtvol=control print volume and debugging randalg=1 if "good" (but non-portable) random numbers should be used, 0 for compatibility restart= if 2, want conversion between wavefunctions if 1, direct restart is allowed (see hdr_check.f) rprimd(3,3)=dimensional primitive translations for real space (bohr) sppoldbl= if 1, no doubling of the number if spins thanks to antiferromagn if 2, deduce nsppol=2 from nsppol=1, using Shubnikov symmetries squeeze=1 if cg_disk is to be used, even with mkmem/=0 symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations tnons(3,nsym)=nonsymmorphic translations for symmetry operations wff1, structure information for input and output files dtset <type(dataset_type)>=all input variables for this dataset
OUTPUT
if ground state format (formeig=0):
eigen(mband*nkpt*nsppol)=eigenvalues (input or init to large number), (Ha)
if respfn format (formeig=1):
eigen(2*mband*mband*nkpt*nsppol)=
matrix of eigenvalues (input or init to large number), (Ha)
Conditional output:
cg_disk(2,mpw*nspinor*mband*mkmem*nsppol)=complex wf array
be careful : an array of size cg(2,npw*nspinor), as used
in the response function code, is not enough !
SIDE EFFECTS
if ground state format (formeig=0):
occ(mband*nkpt*nsppol)=occupations (from disk or left at their initial value)
NOT OUTPUT NOW !
NOTES
occ will not be modified nor output, in the present status of this routine.
WARNINGS
For parallelism : no distinction yet between nban_dp_rd and nband
TODO
THE DESCRIPTION IS TO BE COMPLETELY REVISED, AS THIS ONE COMES FROM inwffil.f
SOURCE
1101 subroutine wfsinp(cg,cg_disk,ecut,ecut0,ecut_eff,eigen,exchn2n3d,& 1102 & formeig,gmet,gmet0,headform0,indkk,indkk0,istwfk,& 1103 & istwfk0,kptns,kptns0,localrdwf,mband,& 1104 & mcg,mcg_disk,mpi_enreg,mpi_enreg0,mpw,mpw0,nband,nban_dp_rd,& 1105 & ngfft,nkassoc,nkpt,nkpt0,npwarr,npwarr0,nspinor,& 1106 & nspinor0,nsppol,nsppol0,nsym,occ,optorth,prtvol,randalg,restart,rprimd,& 1107 & sppoldbl,squeeze,symrel,tnons,wff1) 1108 1109 !Arguments ------------------------------------ 1110 integer, intent(in) :: exchn2n3d,formeig,headform0,localrdwf,mband,mcg,mcg_disk 1111 integer, intent(in) :: mpw,mpw0,nkassoc,nkpt,nkpt0,nspinor,nspinor0,nsppol,nsppol0,nsym 1112 integer, intent(in) :: optorth,prtvol,randalg,restart,sppoldbl,squeeze 1113 real(dp), intent(in) :: ecut,ecut0,ecut_eff 1114 type(MPI_type), intent(inout) :: mpi_enreg,mpi_enreg0 1115 type(wffile_type), intent(inout) :: wff1 1116 integer, intent(in) :: indkk(nkpt*sppoldbl,6),indkk0(nkpt0,nkassoc),istwfk(nkpt) 1117 integer, intent(in) :: istwfk0(nkpt0),nband(nkpt*nsppol),nban_dp_rd(nkpt0*nsppol0) 1118 integer, intent(in) :: ngfft(18),npwarr(nkpt),npwarr0(nkpt0),symrel(3,3,nsym) 1119 real(dp), intent(in) :: gmet(3,3),gmet0(3,3),kptns(3,nkpt),kptns0(3,nkpt0),rprimd(3,3) 1120 real(dp), intent(in) :: tnons(3,nsym) 1121 real(dp), intent(out) :: eigen((2*mband)**formeig*mband*nkpt*nsppol) 1122 real(dp), intent(inout) :: cg(2,mcg),cg_disk(2,mcg_disk) !vz_i pw_ortho 1123 real(dp), intent(inout) :: occ(mband*nkpt*nsppol) 1124 1125 !Local variables------------------------------- 1126 integer :: band_index,band_index_trial,ceksp,debug,dim_eig_k,iband,icg 1127 integer :: icg_disk,icg_trial,idum,ierr,ii,ikassoc,ikassoc_trial,ikpt,ikpt0 1128 integer :: ikpt10,ikpt_trial,ikptsp,ikptsp_old,inplace,isp,isp_max,isppol,isppol0 1129 ! integer :: ipw ! commented out below 1130 integer :: isppol_trial,me,mgfft,my_nspinor,my_nspinor0 1131 integer :: nban_dp_k,nban_dp_rdk,nband_k,nband_rdk,nband_trial,nbd,nbd_max 1132 integer :: ncopy,nkpt_eff,nproc_max,npw0_k,npw_k,npw_ktrial 1133 integer :: read_cg,read_cg_disk,sender,spaceComm 1134 character(len=500) :: message 1135 integer,allocatable :: band_index_k(:,:),icg_k(:,:),kg0_k(:,:),kg_k(:,:) 1136 real(dp) :: tsec(2) 1137 real(dp),allocatable :: eig0_k(:),eig_k(:),occ0_k(:),occ_k(:) 1138 #if defined HAVE_MPI 1139 integer :: iproc,my_ikpt 1140 integer :: tag,test_cycle 1141 integer :: statux(MPI_STATUS_SIZE) 1142 integer,allocatable :: ikassoc_me(:),ikpt_me(:),isppol_me(:),nband_k_me(:) 1143 #endif 1144 integer :: nkpt_max=50 1145 1146 ! ************************************************************************* 1147 1148 !DEBUG 1149 !write(std_out,*)' wfsinp : enter' 1150 !write(std_out,*)' wfsinp : nband=',nband(:) 1151 !write(std_out,*)' wfsinp : nban_dp_rd=',nban_dp_rd(:) 1152 !write(std_out,*)' wfsinp : localrdwf=',localrdwf 1153 !write(std_out,*)' wfsinp : paralbd,formeig=',mpi_enreg%paralbd,formeig 1154 !write(std_out,*)' wfsinp : indkk0(:,1)=',indkk0(:,1) 1155 !ENDDEBUG 1156 1157 call timab(720,1,tsec) 1158 call timab(721,3,tsec) 1159 1160 nkpt_max=50; if(xmpi_paral==1)nkpt_max=-1 1161 nbd_max=size(mpi_enreg%proc_distrb,2) 1162 isp_max=size(mpi_enreg%proc_distrb,3) 1163 1164 !Init mpi_comm 1165 spaceComm=mpi_enreg%comm_cell 1166 nproc_max=xmpi_comm_size(spaceComm) 1167 me=mpi_enreg%me_kpt 1168 sender = 0 1169 1170 #if defined HAVE_MPI 1171 if(localrdwf==0)then 1172 ABI_MALLOC(ikpt_me,(nproc_max)) 1173 ABI_MALLOC(nband_k_me,(nproc_max)) 1174 ABI_MALLOC(ikassoc_me,(nproc_max)) 1175 ABI_MALLOC(isppol_me,(nproc_max)) 1176 end if 1177 #endif 1178 1179 !Check the validity of formeig 1180 if(formeig/=0.and.formeig/=1)then 1181 write(message, '(a,i0,a)' )' formeig=',formeig,' , but the only allowed values are 0 or 1.' 1182 ABI_BUG(message) 1183 end if 1184 1185 my_nspinor =max(1,nspinor /mpi_enreg%nproc_spinor) 1186 my_nspinor0=max(1,nspinor0/mpi_enreg%nproc_spinor) 1187 1188 nkpt_eff=max(nkpt0,nkpt) 1189 if( (prtvol==0.or.prtvol==1) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max 1190 1191 ABI_MALLOC(icg_k,(nkpt,nsppol)) 1192 ABI_MALLOC(band_index_k,(nkpt,nsppol)) 1193 1194 !write(std_out,*)' wfsinp : me,isppol,ikpt,icg_k,band_index_k' 1195 1196 !Compute the locations of the blocks in cg, eig and occ 1197 icg=0 1198 band_index=0 1199 ikpt10=0 1200 1201 do isppol=1,nsppol 1202 do ikpt=1,nkpt 1203 1204 nband_k=nband(ikpt+(isppol-1)*nkpt) 1205 band_index_k(ikpt,isppol)=band_index 1206 1207 #if defined HAVE_MPI 1208 test_cycle=0;nbd=min(nband_k,nbd_max);isp=min(isppol,isp_max) 1209 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nbd,isp,me))test_cycle=1 1210 if(test_cycle==1)then 1211 band_index=band_index+nband_k*(2*nband_k)**formeig 1212 ! In the case this k point does not belong to me, cycle 1213 cycle 1214 end if 1215 #endif 1216 1217 npw_k=npwarr(ikpt) 1218 icg_k(ikpt,isppol)=icg 1219 icg=icg+npw_k*my_nspinor*nband_k 1220 1221 band_index=band_index+nband_k*(2*nband_k)**formeig 1222 ! write(std_out,'(5i8)' )me,isppol,ikpt,icg_k(ikpt,isppol),band_index_k(ikpt,isppol) 1223 end do ! End k point loop 1224 end do ! End spin loop 1225 1226 band_index=0 1227 ikptsp_old=0 1228 1229 !DEBUG 1230 !write(std_out,*)' wfsinp: before loop' 1231 !write(std_out,*)' nsppol0,nsppol,nkpt0',nsppol0,nsppol,nkpt0 1232 !write(std_out,*)' mpw,mgfft,mpw,mpw0',mpw,mgfft,mpw,mpw0 1233 !ENDDEBUG 1234 1235 mgfft=maxval(ngfft(1:3)) 1236 if(squeeze==1)then 1237 ABI_MALLOC(kg_k,(3,mpw)) 1238 ABI_MALLOC(kg0_k,(3,mpw0)) 1239 end if 1240 1241 eigen(:)=0.0_dp 1242 !occ(:)=0.0_dp 1243 1244 call timab(721,2,tsec) 1245 1246 !Loop over spins 1247 !For the time being, do not allow nsppol=2 to nspinor=2 conversion 1248 !MT 20110707: this can be done by a fake call to the routine: see inwffil 1249 do isppol0=1,min(nsppol0,nsppol) 1250 1251 ! Loop on k points : get the cg then eventually write on unwfnow 1252 do ikpt0=1,nkpt0 1253 1254 call timab(722,1,tsec) 1255 1256 nban_dp_rdk=nban_dp_rd(ikpt0+(isppol0-1)*nkpt0) 1257 1258 ! DEBUG 1259 ! write(std_out,*)' wfsinp: ikpt0,isppol0,nkpt0=',ikpt0,isppol0,nkpt0 1260 ! write(std_out,*)' nban_dp_rdk=',nban_dp_rdk 1261 ! ENDDEBUG 1262 1263 npw0_k=npwarr0(ikpt0) 1264 if(ikpt0<=nkpt_eff)then 1265 write(message,'(a,a,2i4)')ch10,' wfsinp: inside loop, init ikpt0,isppol0=',ikpt0,isppol0 1266 call wrtout(std_out,message) 1267 end if 1268 1269 ! Must know whether this k point is needed, and in which 1270 ! block (ikpt, isppol), the wavefunction is to be placed. 1271 ! Select the one for which the number of bands is the biggest. 1272 ikpt=0 1273 isppol=0 1274 ikassoc=0 1275 nband_k=0 1276 #if defined HAVE_MPI 1277 if(localrdwf==0)then 1278 nband_k_me(:)=0 1279 ikpt_me(:)=0 1280 isppol_me(:)=0 1281 ikassoc_me(:)=0 1282 nband_k_me(:)=0 1283 end if 1284 #endif 1285 1286 do isppol_trial=1,nsppol 1287 1288 if(nsppol==2 .and. nsppol0==2 .and. isppol0/=isppol_trial)cycle 1289 1290 do ikassoc_trial=1,nkassoc 1291 1292 ikpt_trial=indkk0(ikpt0,ikassoc_trial) 1293 if(sppoldbl==2)then 1294 if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle 1295 if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle 1296 if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt 1297 end if 1298 1299 #if defined HAVE_MPI 1300 if(localrdwf==1)then 1301 if(ikpt_trial/=0)then 1302 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 1303 nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max) 1304 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me))ikpt_trial=0 1305 end if 1306 end if 1307 #endif 1308 1309 if(ikpt_trial/=0)then 1310 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 1311 if(nband_k<nband_trial)then 1312 nband_k=nband_trial ; ikpt=ikpt_trial ; isppol=isppol_trial 1313 ikassoc=ikassoc_trial 1314 end if 1315 1316 #if defined HAVE_MPI 1317 if(localrdwf==0)then 1318 do iproc=1,nproc_max 1319 my_ikpt=1 1320 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 1321 nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max) 1322 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,(iproc-1))) my_ikpt=0 1323 if(my_ikpt/=0)then 1324 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 1325 if(nband_k_me(iproc)<nband_trial)then 1326 nband_k_me(iproc)=nband_trial 1327 ikpt_me(iproc)=ikpt_trial 1328 isppol_me(iproc)=isppol_trial 1329 ikassoc_me(iproc)=ikassoc_trial 1330 end if 1331 end if 1332 end do 1333 end if 1334 #endif 1335 1336 end if 1337 1338 end do ! ikassoc_trial 1339 end do ! isppol_trial 1340 1341 ! DEBUG 1342 ! write(std_out,*)' wfsinp : me,select isppol,ikpt=',me,isppol,ikpt 1343 #if defined HAVE_MPI 1344 ! write(std_out,*)' wfsinp : me,ikpt_me(:)=',me,ikpt_me(:) 1345 #endif 1346 ! write(std_out,*)' wfsinp : me,isppol_me(:)=',me,isppol_me(:) 1347 ! stop 1348 ! ENDDEBUG 1349 1350 call timab(722,2,tsec) 1351 1352 ! If the wavefunction block to be read is interesting ... 1353 if (ikpt/=0)then 1354 1355 call timab(723,3,tsec) 1356 sender = me 1357 nband_k=nband(ikpt+(isppol-1)*nkpt) 1358 npw_k=npwarr(ikpt) 1359 1360 #if defined HAVE_MPI 1361 if (localrdwf==1.or.(localrdwf==0.and.me==0)) then 1362 1363 if(ikpt<=nkpt_eff)then 1364 write(message,'(a,i6,a,i8,a,i4,a,i4)') & 1365 & ' wfsinp: treating ',nband_k,' bands with npw=',npw_k,' for ikpt=',ikpt,' by node ',me 1366 call wrtout(std_out,message) 1367 else if(ikpt==nkpt_eff+1)then 1368 call wrtout(std_out,' wfsinp: prtvol=0 or 1, do not print more k-points.') 1369 end if 1370 1371 end if 1372 #endif 1373 1374 nband_rdk=nban_dp_rdk 1375 if(squeeze==1)nband_rdk=(nban_dp_rdk/nspinor0)*nspinor 1376 if(formeig==0)then 1377 ABI_MALLOC(eig_k,(nband_rdk)) 1378 ABI_MALLOC(occ_k,(nband_rdk)) 1379 ABI_MALLOC(eig0_k,(nban_dp_rdk)) 1380 ABI_MALLOC(occ0_k,(nban_dp_rdk)) 1381 dim_eig_k=nband_rdk 1382 else if(formeig==1)then 1383 ABI_MALLOC(eig_k,(2*nband_rdk*nband_rdk)) 1384 ABI_MALLOC(eig0_k,(2*nban_dp_rdk*nban_dp_rdk)) 1385 dim_eig_k=2*nband_rdk*nband_rdk 1386 ABI_MALLOC(occ0_k,(0)) 1387 ABI_MALLOC(occ_k,(0)) 1388 else 1389 ABI_MALLOC(occ0_k,(0)) 1390 ABI_MALLOC(eig0_k,(0)) 1391 ABI_MALLOC(occ_k,(0)) 1392 ABI_MALLOC(eig_k,(0)) 1393 end if 1394 eig_k(:)=0.0_dp 1395 eig0_k(:)=0.0_dp 1396 1397 ! Generate or read the cg for this k point 1398 ! Either read into cg, or read into cg_disk 1399 read_cg=1 ; read_cg_disk=0 1400 if(squeeze==1)then 1401 read_cg=0 ; read_cg_disk=1 1402 end if 1403 #if defined HAVE_MPI 1404 if(localrdwf==0)then 1405 read_cg=0 1406 read_cg_disk=0 1407 ! XG20040106 The following condition is correct 1408 if(me==0)read_cg_disk=1 1409 end if 1410 #endif 1411 1412 icg=0 1413 if(read_cg==1)icg=icg_k(ikpt,isppol) 1414 1415 ! DEBUG 1416 ! write(std_out,*)' wfsinp: before initwf',wff1%offwff 1417 ! write(std_out,*)' wfsinp: me,read_cg,read_cg_disk=',me,read_cg,read_cg_disk 1418 ! write(std_out,*)' wfsinp: nban_dp_rdk=',nban_dp_rdk 1419 ! ENDDEBUG 1420 call timab(723,2,tsec) 1421 call timab(724,3,tsec) 1422 1423 if(read_cg_disk==1)then 1424 call initwf (cg_disk,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,& 1425 & isppol0,mcg_disk,mpi_enreg0, & 1426 & nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1) 1427 end if 1428 1429 if(read_cg==1)then 1430 call initwf (cg,eig0_k,formeig,headform0,icg,ikpt0,ikptsp_old,& 1431 & isppol0,mcg,mpi_enreg0,& 1432 & nban_dp_rdk,nkpt0,npw0_k,my_nspinor0,occ0_k,wff1) 1433 end if 1434 1435 call timab(724,2,tsec) 1436 call timab(725,3,tsec) 1437 1438 nban_dp_k=min(nban_dp_rdk,(nband_k/nspinor)*nspinor0) 1439 ! This band_index is defined BEFORE the eventual redefinition 1440 ! of ikpt and isppol, needed when localrdwf==0 in parallel 1441 band_index=band_index_k(ikpt,isppol) 1442 ! DEBUG 1443 ! write(std_out,*)' wfsinp: me,cg_disk(:,1)=',me,cg_disk(:,1) 1444 ! ENDDEBUG 1445 1446 ! DEBUG 1447 ! if(me==0 .and. ikpt0==1)then 1448 ! write(std_out,*)' wfsinp : cg array, before trial, ikpt0=',ikpt0 1449 ! do ipw=1,15 1450 ! write(std_out,'(i4,2es20.10)' )ipw,cg(:,ipw) 1451 ! end do 1452 ! end if 1453 ! ENDDEBUG 1454 1455 1456 #if defined HAVE_MPI 1457 if(localrdwf==0)then 1458 ! Warning: In that case , not yet // on nspinors 1459 ! Transmit to each of the other processors, when needed 1460 if(nproc_max>=2)then 1461 do iproc=2,nproc_max 1462 ! Only me=0 and me=iproc-1 are concerned by this 1463 if(me==0 .or. me==iproc-1)then 1464 1465 ikpt=ikpt_me(iproc) 1466 isppol=isppol_me(iproc) 1467 1468 if(ikpt/=0)then 1469 ! In this case, processor iproc-1 needs the data 1470 ! Generate a common tag 1471 tag=256*(ikpt-1)+iproc+1 1472 if(isppol==2)tag=-tag 1473 nband_k=nband(ikpt+(isppol-1)*nkpt) 1474 npw_k=npwarr(ikpt) 1475 ! SEND 1476 if(me==0)then 1477 write(std_out,*)'SENDWFSINP ',me 1478 call MPI_SEND(cg_disk,2*npw_k*my_nspinor*nband_k,& 1479 & MPI_DOUBLE_PRECISION,iproc-1,tag,spaceComm,ierr) 1480 end if 1481 ! RECEIVE 1482 if(me==iproc-1)then 1483 call MPI_RECV(cg_disk,2*npw_k*my_nspinor*nband_k,& 1484 & MPI_DOUBLE_PRECISION,0,tag,spaceComm,statux,ierr) 1485 icg=icg_k(ikpt,isppol) 1486 if(squeeze==0)then 1487 cg(:,icg+1:icg+npw_k*my_nspinor*nband_k)=& 1488 & cg_disk(:,1:npw_k*my_nspinor*nband_k) 1489 end if 1490 ikassoc=ikassoc_me(iproc) 1491 end if 1492 end if 1493 1494 end if 1495 end do ! iproc 1496 end if 1497 1498 ! DEBUG 1499 ! write(std_out,*)' wfsinp: me, iproc loop finished',me 1500 ! ENDDEBUG 1501 1502 ! Take care of me=0 needing the data 1503 if (me==0) then 1504 ikpt=ikpt_me(me+1) 1505 isppol=isppol_me(me+1) 1506 if(ikpt/=0 )then 1507 nband_k=nband(ikpt+(isppol-1)*nkpt) 1508 npw_k=npwarr(ikpt) 1509 ! I am the master node, and I might need my own data 1510 icg=icg_k(ikpt,isppol) 1511 if(squeeze==0)then 1512 ! Copy from cg_disk to cg 1513 cg(:,1+icg:npw_k*my_nspinor*nband_k+icg)= & 1514 & cg_disk(:,1:npw_k*my_nspinor*nband_k) 1515 end if 1516 ikassoc=ikassoc_me(me+1) 1517 end if 1518 end if 1519 ! For the eigenvalues and occ, the transmission is much easier to write ! 1520 call MPI_BCAST(eig0_k,nban_dp_rdk*(2*nban_dp_rdk)**formeig ,& 1521 & MPI_DOUBLE_PRECISION,0,spaceComm,ierr) 1522 end if 1523 #endif 1524 1525 if(formeig==0)then 1526 ! The transfer from eig0_k to eig_k uses nban_dp_rdk, which contains 1527 ! the maximal information. 1528 if(nspinor0==nspinor .or. squeeze==0)then 1529 eig_k(1:nban_dp_rdk)=eig0_k(1:nban_dp_rdk) 1530 occ_k(1:nban_dp_rdk)=occ0_k(1:nban_dp_rdk) 1531 else if(nspinor0==1 .and. nspinor==2)then 1532 do iband=1,nban_dp_rdk 1533 eig_k(2*iband )=eig0_k(iband) 1534 eig_k(2*iband-1)=eig0_k(iband) 1535 occ_k(2*iband )=occ0_k(iband)*0.5_dp 1536 occ_k(2*iband-1)=occ0_k(iband)*0.5_dp 1537 end do 1538 else if(nspinor0==2 .and. nspinor==1)then 1539 do iband=1,nban_dp_rdk 1540 eig_k(iband)=eig0_k(2*iband-1) 1541 occ_k(iband)=occ0_k(2*iband-1)*2.0_dp 1542 end do 1543 end if 1544 1545 ! DEBUG 1546 ! write(std_out,*)' wfsinp: me,band_index,ikpt,isppol',me,band_index,ikpt,isppol 1547 ! ENDDEBUG 1548 1549 ! The transfer to eigen uses nban_dp_k, that is bound by the number 1550 ! of bands for this k point. 1551 ncopy=min(dim_eig_k,(nban_dp_k/nspinor0)*nspinor) 1552 eigen(1+band_index:ncopy+band_index)=eig_k(1:ncopy) 1553 ! The transfer of occ should be done. 1554 1555 #if defined HAVE_MPI 1556 if(localrdwf==0 .and. ikpt/=0)then 1557 ! Do not forget : ikpt,isppol were redefined ... 1558 band_index=band_index_k(ikpt,isppol) 1559 eigen(1+band_index:(nban_dp_k/nspinor0)*nspinor+band_index) = eig_k(1:(nban_dp_k/nspinor0)*nspinor) 1560 ! The transfer of occ should be done. 1561 end if 1562 #endif 1563 1564 else if(formeig==1)then 1565 call wrtout(std_out,'wfsinp: transfer of first-order eigs not yet coded!',"COLL") 1566 end if 1567 1568 ! DEBUG 1569 ! write(std_out,*)' wfsinp : me,transferred eig_k',me 1570 ! write(std_out,*)' me,mkmem,nsppol,nsppol0,isppol',me,mkmem,nsppol,nsppol0,isppol 1571 ! write(std_out,*)' me,nkassoc,ikassoc',me,nkassoc,ikassoc 1572 ! ENDDEBUG 1573 1574 call timab(725,2,tsec) 1575 1576 ! Write to disk if appropriate 1577 ! The coding has to be done ... here, only fragments ... 1578 call timab(727,3,tsec) 1579 1580 do isppol_trial=1,nsppol 1581 if(nsppol==2 .and. nsppol0==2 .and. isppol_trial/=isppol)cycle 1582 do ikassoc_trial=1,nkassoc 1583 1584 ! DEBUG 1585 ! write(std_out,*)' wfsinp: me, for ikassoc,isppol',& 1586 ! & me,ikassoc,isppol 1587 ! write(std_out,*)' wfsinp: me, try ikassoc_trial,isppol_trial,nband_k',& 1588 ! & me,ikassoc_trial,isppol_trial,nband_k 1589 ! ENDDEBUG 1590 1591 ! No conversion is to be done : it will be converted in newkpt 1592 ! If squeeze==0, the block with the ikpt corresponding to ikassoc, 1593 ! and with isppol, contains the wavefunction already 1594 if( ikassoc_trial/=ikassoc .or. & 1595 & (isppol_trial/=isppol .and. sppoldbl==1 ).or. & 1596 & squeeze==1 )then 1597 1598 ikpt_trial=indkk0(ikpt0,ikassoc_trial) 1599 if(sppoldbl==2)then 1600 if(isppol_trial==1 .and. ikpt_trial>nkpt)cycle 1601 if(isppol_trial==2 .and. ikpt_trial<=nkpt)cycle 1602 if(isppol_trial==2)ikpt_trial=ikpt_trial-nkpt 1603 end if 1604 1605 1606 #if defined HAVE_MPI 1607 if(ikpt_trial/=0)then 1608 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 1609 nbd=min(nband_trial,nbd_max);isp=min(isppol_trial,isp_max) 1610 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt_trial,1,nbd,isp,me)) ikpt_trial=0 1611 end if 1612 #endif 1613 1614 if(ikpt_trial/=0 .and. ikpt_trial<=nkpt_eff)then 1615 write(message,'(2a,2i5)')ch10,' wfsinp: transfer to ikpt_trial,isppol_trial=',ikpt_trial,isppol_trial 1616 call wrtout(std_out,message) 1617 end if 1618 1619 if(ikpt_trial/=0)then 1620 icg_trial=icg_k(ikpt_trial,isppol_trial) 1621 band_index_trial=band_index_k(ikpt_trial,isppol_trial) 1622 nband_trial=nband(ikpt_trial+(isppol_trial-1)*nkpt) 1623 nban_dp_k=min(nban_dp_rdk,(nband_trial/nspinor)*nspinor0) 1624 1625 if(squeeze==0)then 1626 ! GMR: modified to avoid compiler bug 1627 ! cg(:,1+icg_trial:npw0_k*my_nspinor0*nband_trial+icg_trial)=& 1628 ! & cg(:,1+icg:npw0_k*my_nspinor0*nband_trial+icg) 1629 do ii=1,npw0_k*my_nspinor0*nband_trial 1630 cg(:,ii+icg_trial)=cg(:,ii+icg) 1631 end do 1632 ! GMR 1633 if(formeig==0)then 1634 eigen(1+band_index_trial:nban_dp_k+band_index_trial)=eig_k(1:nban_dp_k) 1635 ! occ(1+band_index_trial:nban_dp_k+band_index_trial)=& 1636 ! & occ_k(1:nban_dp_k) 1637 end if 1638 ! RF transfer of eigenvalues still to be coded 1639 else if(squeeze==1)then 1640 npw_ktrial=npwarr(ikpt_trial) 1641 nband_k=(nban_dp_k/nspinor0)*nspinor 1642 ! Conversion to be done 1643 ceksp=0 ; debug=0 ; icg_disk=0 ; idum=0 ; inplace=0 1644 ! Note that this routine also convert eig and occ 1645 ! even if the conversion had already been done 1646 1647 1648 call wfconv(ceksp,cg_disk,cg,debug,ecut0,ecut,ecut_eff,& 1649 & eig0_k,eig_k,exchn2n3d,formeig,gmet0,gmet,& 1650 & icg_disk,icg_trial,ikpt0,ikpt10,ikpt_trial,indkk,& 1651 & inplace,isppol_trial,istwfk0,istwfk,& 1652 & kg0_k,kg_k,kptns0,kptns,nban_dp_rdk,nband_rdk,& 1653 & mcg_disk,mcg,mpi_enreg0,mpi_enreg,mpw0,mpw,& 1654 & nban_dp_rdk,nband_trial,ngfft,ngfft,nkpt0,nkpt,& 1655 & npw0_k,npw_ktrial,nspinor0,nspinor,nsym,& 1656 & occ0_k,occ_k,optorth,randalg,restart,rprimd,& 1657 & sppoldbl,symrel,tnons) 1658 1659 ! DEBUG 1660 ! write(std_out,*)' wfsinp: ikpt_trial=',ikpt_trial 1661 ! write(std_out,*)' nband_k,band_index_trial',nband_k,band_index_trial 1662 ! write(std_out,*)eig0_k(1:nban_dp_k) 1663 ! write(std_out,*)eig_k(1:nband_k) 1664 ! ENDDEBUG 1665 eigen(1+band_index_trial:nband_k+band_index_trial)=eig_k(1:nband_k) 1666 ! occ(1+band_index_trial:nband_k+band_index_trial)=& 1667 ! & occ_k(1:nband_k) 1668 ! RF transfer of eigenvalues still to be coded 1669 1670 ! Endif squeeze==1 1671 end if 1672 1673 ! DEBUG 1674 ! if(ikpt_trial==2)then 1675 ! write(std_out,*)' wfsinp: iband,ipw,cg for ikpt_trial=2' 1676 ! write(std_out,*)' nband_trial,npw_ktrial=',nband_trial,npw_ktrial 1677 ! do iband=1,nband_trial 1678 ! do ipw=1,npw_ktrial 1679 ! write(std_out,'(2i5,2es16.6)' )& 1680 ! & iband,ipw,cg(:,ipw+(iband-1)*npw_ktrial+icg_trial) 1681 ! end do 1682 ! end do 1683 ! end if 1684 ! ENDDEBUG 1685 1686 ! End if ikpt_trial/=0 1687 end if 1688 1689 ! End if ikpt_trial already initialized 1690 end if 1691 1692 end do ! ikassoc_trial 1693 end do ! isppol_trial 1694 1695 call timab(727,2,tsec) 1696 1697 ABI_FREE(eig_k) 1698 ABI_FREE(eig0_k) 1699 !if(formeig==0) then 1700 ABI_FREE(occ_k) 1701 ABI_FREE(occ0_k) 1702 !end if 1703 1704 end if ! End the condition of need of this k point 1705 end do ! End of the k loop 1706 end do ! End of spin loop 1707 1708 #if defined HAVE_MPI 1709 call timab(67,1,tsec) 1710 1711 !Still need to skip last k points to read history (MS) 1712 !WARNING : not yet for formeig=1, but is it needed ? 1713 if(formeig==0)then 1714 do ikptsp=ikptsp_old+1,nkpt0*nsppol0 1715 isppol=1 ; if(ikptsp>nkpt0)isppol=2 1716 ikpt=ikptsp-nkpt0*(isppol-1) 1717 call WffReadSkipK(formeig,headform0,ikpt,isppol,mpi_enreg,wff1) 1718 end do 1719 end if 1720 1721 !Transmit eigenvalues. This routine works in both localrdwf=0 or 1 cases. 1722 call pareigocc(eigen,formeig,localrdwf,mpi_enreg,mband,nband,nkpt,nsppol,occ,1) 1723 1724 1725 if(localrdwf==0)then 1726 ABI_FREE(ikpt_me) 1727 ABI_FREE(nband_k_me) 1728 ABI_FREE(ikassoc_me) 1729 ABI_FREE(isppol_me) 1730 end if 1731 1732 call timab(67,2,tsec) 1733 #endif 1734 1735 !**************************************************************************** 1736 1737 if(squeeze==1)then 1738 ABI_FREE(kg_k) 1739 ABI_FREE(kg0_k) 1740 end if 1741 1742 1743 !DEBUG 1744 !if(me==0)then 1745 !write(std_out,*)' wfsinp : cg array=' 1746 !icg=0 1747 !do isppol=1,nsppol 1748 !do ikpt=1,1 1749 !nband_k=nband(ikpt+(isppol-1)*nkpt) 1750 !npw_k=npwarr(ikpt) 1751 !do iband=1,nband_k 1752 !write(std_out,*)' new band, icg=',icg 1753 !do ipw=1,npw_k 1754 !write(std_out,'(4i4,2es20.10)' )isppol,ikpt,iband,ipw,cg(:,icg+ipw) 1755 !end do 1756 !icg=icg+npw_k 1757 !end do 1758 !end do 1759 !end do 1760 !end if 1761 !if(ireadwf==1)stop 1762 !write(std_out,*)' wfsinp : eigen array=' 1763 !do ikpt=1,nkpt 1764 !do iband=1,mband 1765 !write(std_out,*)'ikpt,iband,eigen',ikpt,iband,eigen(iband+(ikpt-1)*mband) 1766 !end do 1767 !end do 1768 !ENDDEBUG 1769 1770 ABI_FREE(icg_k) 1771 ABI_FREE(band_index_k) 1772 1773 call timab(720,2,tsec) 1774 1775 end subroutine wfsinp