TABLE OF CONTENTS
ABINIT/etotfor [ Functions ]
NAME
etotfor
FUNCTION
This routine is called to compute the total energy and various parts of it. The routine computes -if requested- the forces.
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx dtefield <type(efield_type)> = variables related to Berry phase dtset <type(dataset_type)>=all input variables in this dataset | berryopt = 4: electric field is on -> add the contribution of the | - \Omega E.P term to the total energy | /= 4: electric field is off | bfield = cartesian coordinates of magnetic field in atomic units | efield = cartesian coordinates of the electric field in atomic units | iatfix(3,natom)=1 for frozen atom along some direction, 0 for unfrozen | ionmov=governs the movement of atoms (see help file) | densfor_pred=governs the mixed electronic-atomic part of the preconditioner | natom=number of atoms in cell. | nconeq=number of atomic constraint equations | nspden=number of spin-density components | nsym=number of symmetry elements in space group | occopt=option for occupancies | prtvol=integer controlling volume of printed output | tsmear=smearing energy or temperature (if metal) | typat(natom)=type integer for each atom in cell | wtatcon(3,natom,nconeq)=weights for atomic constraints gmet(3,3)=metric tensor for G vecs (in bohr**-2) fock <type(fock_type)>= quantities to calculate Fock exact exchange grchempottn(3,natom)=grads of spatially-varying chemical potential energy (hartree) grcondft(3,natom)=grads of constrained DFT energy (hartree) grewtn(3,natom)=grads of Ewald energy (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree) gsqcut=cutoff on (k+G)^2 (bohr^-2) indsym(4,nsym,natom)=indirect indexing array for atom labels kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor nattyp(ntypat)=number of atoms of each type nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc nhat(nfft,nspden*usepaw)= -PAW only- compensation density nkxc=second dimension of the array kxc, see rhotoxc.f for a description ntypat=number of types of atoms in unit cell. nvresid(nfft,nspden)=potential or density residual n1xccc=dimension of xccc1d ; 0 if no XC core correction is used n3xccc=dimension of the xccc3d array (0 or nfft). optene=option for the computation of total energy (-1=no computation; 0=direct scheme; 1=double-counting scheme) optforces=option for the computation of forces optres=0 if residual array (nvresid) contains the potential residual =1 if residual array (nvresid) contains the density residual pawang <type(pawang_type)>=paw angular mesh and related data pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information. psps <type(pseudopotential_type)>=variables related to pseudopotentials rhog(2,nfft)=array for Fourier transform of electron density rhor(nfft,nspden)=array for electron density in electrons/bohr**3 rprimd(3,3)=dimensional primitive translations in real space (bohr) symrec(3,3,nsym)=symmetry operations in reciprocal space ucvol=unit cell volume usepaw= 0 for non paw calculation; =1 for paw calculation vhartr(nfft)=array for holding Hartree potential vpsp(nfft)=array for holding local psp vxc(nfft,nspden)=array for holding XC potential vxctau(nfftf,dtset%nspden,4*usekden)]=derivative of XC energy density wrt kinetic energy density (metaGGA cases) xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
deltae=change in total energy between the previous and present SCF cycle etotal=total energy (hartree) ===== if optforces==1 diffor=maximum absolute change in component of forces between present and previous SCF cycle. favg(3)=mean of fcart before correction for translational symmetry fcart(3,natom)=cartesian forces from gred (hartree/bohr) gred(3,natom)=symmetrized form of grtn (grads of Etot) (hartree) gresid(3,natom)=forces due to the residual of the density/potential grhf(3,natom)=Hellman-Feynman derivatives of the total energy grxc(3,natom)=d(Exc)/d(xred) derivatives (0 without core charges) maxfor=maximum absolute value of force synlgr(3,natom)=symmetrized form of grads of Enl (hartree)
SIDE EFFECTS
Input/Output: elast=previous value of the energy, needed to compute deltae, then updated. electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation energies <type(energies_type)>=all part of total energy. | entropy(IN)=entropy due to the occupation number smearing (if metal) | e_localpsp(IN)=local psp energy (hartree) | e_eigenvalues(IN)=Sum of the eigenvalues - Band energy (Hartree) | e_chempot(IN)=energy from spatially varying chemical potential (hartree) | e_ewald(IN)=Ewald energy (hartree) | e_vdw_dftd(IN)=VdW DFT-D energy | e_hartree(IN)=Hartree part of total energy (hartree units) | e_corepsp(IN)=psp core-core energy | e_hybcomp_E0(IN)=energy compensation energy for the hybrid functionals at frozen density | e_hybcomp_v0(IN)=potential compensation energy for the hybrid functionals at frozen density | e_hybcomp_v (IN)=potential compensation energy for the hybrid functionals at self-consistent density | e_kinetic(IN)=kinetic energy part of total energy. | e_nlpsp_vfock(IN)=nonlocal psp + potential Fock ACE part of total energy. | e_nucdip(IN)=energy due to array of nuclear magnetic dipoles | e_xc(IN)=exchange-correlation energy (hartree) | e_xcdc(IN)=exchange-correlation double-counting energy (hartree) | e_paw(IN)=PAW spherical part energy | e_pawdc(IN)=PAW spherical part double-counting energy | e_elecfield(OUT)=the term of the energy functional that depends explicitely | on the electric field: enefield = -ucvol*E*P | e_magfield(OUT)=the term of the energy functional that depends explicitely | on the magnetic field: e_magfield = -ucvol*E*P | e_entropy(OUT)=entropy energy due to the occupation number smearing (if metal) | this value is %entropy * dtset%tsmear (hartree). ===== if optforces==1 forold(3,natom)=cartesian forces of previous SCF cycle (hartree/bohr) grnl(3*natom)=gradients of Etot due to nonlocal contributions Input for norm-conserving psps, output for PAW ===== if psps%usepaw==1 pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data (gradients of rhoij for each atom with respect to atomic positions are computed here)
NOTES
In case of PAW calculations:
All computations are done on the fine FFT grid.
All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
All arrays (densities/potentials...) are computed on this fine FFT grid.
! Developpers have to be careful when introducing others arrays:
they have to be stored on the fine FFT grid.
In case of norm-conserving calculations the FFT grid is the usual FFT grid.
SOURCE
2524 subroutine etotfor(atindx1,deltae,diffor,dtefield,dtset,& 2525 & elast,electronpositron,energies,& 2526 & etotal,favg,fcart,fock,forold,gred,gmet,grchempottn,grcondft,gresid,grewtn,grhf,grnl,grvdw,& 2527 & grxc,gsqcut,extfpmd,indsym,kxc,maxfor,mgfft,mpi_enreg,my_natom,nattyp,& 2528 & nfft,ngfft,ngrvdw,nhat,nkxc,ntypat,nvresid,n1xccc,n3xccc,optene,optforces,optres,& 2529 & pawang,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,red_ptot,psps,rhog,rhor,rmet,rprimd,& 2530 & symrec,synlgr,ucvol,usepaw,vhartr,vpsp,vxc,vxctau,wvl,wvl_den,xccc3d,xred) 2531 2532 !Arguments ------------------------------------ 2533 !scalars 2534 integer,intent(in) :: my_natom,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nkxc,ntypat,optene,optforces 2535 integer,intent(in) :: optres,usepaw 2536 real(dp),intent(in) :: gsqcut 2537 real(dp),intent(inout) :: elast,ucvol 2538 real(dp),intent(out) :: deltae,diffor,etotal,maxfor 2539 type(MPI_type),intent(in) :: mpi_enreg 2540 type(efield_type),intent(in) :: dtefield 2541 type(dataset_type),intent(in) :: dtset 2542 type(electronpositron_type),pointer :: electronpositron 2543 type(energies_type),intent(inout) :: energies 2544 type(extfpmd_type),pointer,intent(inout) :: extfpmd 2545 type(pawang_type),intent(in) :: pawang 2546 type(pseudopotential_type),intent(in) :: psps 2547 type(wvl_internal_type), intent(in) :: wvl 2548 type(wvl_denspot_type), intent(inout) :: wvl_den 2549 type(fock_type),pointer, intent(inout) :: fock 2550 !arrays 2551 integer,intent(in) :: atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom) 2552 integer,intent(in) :: nattyp(ntypat),ngfft(18),symrec(3,3,dtset%nsym) 2553 real(dp),intent(in) :: gmet(3,3),grchempottn(3,dtset%natom),grcondft(3,dtset%natom) 2554 real(dp),intent(in) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc) 2555 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),red_ptot(3) 2556 real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,dtset%nspden),rmet(3,3) 2557 real(dp),intent(in) :: vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden) 2558 real(dp),intent(in) :: vxctau(nfft,dtset%nspden,4*dtset%usekden) 2559 real(dp),intent(in) :: xccc3d(n3xccc) 2560 real(dp),intent(inout) :: forold(3,dtset%natom),grnl(3*dtset%natom) 2561 real(dp),intent(inout) :: nhat(nfft,dtset%nspden*psps%usepaw) 2562 real(dp),intent(inout),target :: nvresid(nfft,dtset%nspden) 2563 real(dp),intent(inout) :: xred(3,dtset%natom) 2564 real(dp),intent(out) :: favg(3),gred(3,dtset%natom) 2565 real(dp),intent(inout) :: fcart(3,dtset%natom) 2566 real(dp),intent(inout) :: rprimd(3,3) 2567 real(dp),intent(out) :: gresid(3,dtset%natom),grhf(3,dtset%natom) 2568 real(dp),intent(inout) :: grxc(3,dtset%natom) 2569 real(dp),intent(out) :: synlgr(3,dtset%natom) 2570 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 2571 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw) 2572 type(pawrad_type),intent(in) :: pawrad(ntypat*psps%usepaw) 2573 type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw) 2574 2575 !Local variables------------------------------- 2576 !scalars 2577 integer :: comm_grid,dimnhat,ifft,ipositron,ispden,optgr,optgr2,option,optnc,optstr,optstr2,iir,jjr,kkr 2578 logical :: apply_residual 2579 real(dp) :: eenth,ucvol_ 2580 !arrays 2581 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 2582 real(dp) :: tsec(2),A(3,3),A1(3,3),A_new(3,3),efield_new(3) 2583 real(dp) :: dummy(0),nhat_dum(0,0) 2584 real(dp),allocatable :: vlocal(:,:) 2585 real(dp), ABI_CONTIGUOUS pointer :: resid(:,:) 2586 2587 ! ********************************************************************* 2588 2589 call timab(80,1,tsec) 2590 2591 ipositron=electronpositron_calctype(electronpositron) 2592 2593 if (optene>-1) then 2594 2595 ! Add the contribution of extfpmd to the entropy 2596 if(associated(extfpmd)) then 2597 energies%entropy=energies%entropy+extfpmd%entropy 2598 end if 2599 2600 ! When the finite-temperature VG broadening scheme is used, 2601 ! the total entropy contribution "tsmear*entropy" has a meaning, 2602 ! and gather the two last terms of Eq.8 of VG paper 2603 ! Warning : might have to be changed for fixed moment calculations 2604 if(dtset%occopt>=3 .and. dtset%occopt<=8) then 2605 if (abs(dtset%tphysel) < tol10) then 2606 energies%e_entropy = - dtset%tsmear * energies%entropy 2607 else 2608 energies%e_entropy = - dtset%tphysel * energies%entropy 2609 end if 2610 else 2611 energies%e_entropy = zero 2612 end if 2613 2614 ! Turn it into an electric enthalpy, refer to Eq.(33) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]], 2615 ! the missing volume is added here 2616 energies%e_elecfield=zero 2617 if ((dtset%berryopt==4.or.dtset%berryopt==14).and.ipositron/=1) then 2618 energies%e_elecfield=-dot_product(dtset%red_efieldbar,red_ptot) !!ebar_i p_i 2619 eenth=zero 2620 do iir=1,3 2621 do jjr=1,3 2622 eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !!g^{-1})_ij ebar_i ebar_j 2623 end do 2624 end do 2625 energies%e_elecfield=energies%e_elecfield-eenth*ucvol/(8._dp*pi) 2626 end if 2627 2628 ! Turn it into an internal energy, refer to Eq.(36) of Suppl. of Nat. Phys. paper (5,304,2009) [[cite:Stengel2009]], 2629 ! but a little different: U=E_ks + (vol/8*pi) * g^{-1})_ij ebar_i ebar_j 2630 if ((dtset%berryopt==6.or.dtset%berryopt==16).and.ipositron/=1) then 2631 energies%e_elecfield=zero 2632 eenth=zero 2633 do iir=1,3 2634 do jjr=1,3 2635 eenth=eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! g^{-1})_ij ebar_i ebar_j 2636 end do 2637 end do 2638 energies%e_elecfield=energies%e_elecfield+eenth*ucvol/(8._dp*pi) 2639 end if 2640 2641 ! Calculate internal energy and electric enthalpy for mixed BC case. 2642 if (dtset%berryopt==17.and.ipositron/=1) then 2643 energies%e_elecfield=zero 2644 A(:,:)=(four_pi/ucvol)*rmet(:,:) 2645 A1(:,:)=A(:,:) ; A_new(:,:)=A(:,:) 2646 efield_new(:)=dtset%red_efield(:) 2647 eenth=zero 2648 do kkr=1,3 2649 if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction 2650 ! step 1 add -ebar*p 2651 eenth=eenth-dtset%red_efieldbar(kkr)*red_ptot(kkr) 2652 ! step 2 chang to e_new (change e to ebar) 2653 efield_new(kkr)=dtset%red_efieldbar(kkr) 2654 ! step 3 chang matrix A to A1 2655 do iir=1,3 2656 do jjr=1,3 2657 if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr) 2658 if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) & 2659 & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr) 2660 if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr) 2661 end do 2662 end do 2663 A(:,:)=A1(:,:) ; A_new(:,:)=A1(:,:) 2664 end if 2665 end do ! end fo kkr 2666 do iir=1,3 2667 do jjr=1,3 2668 eenth= eenth+half*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr) 2669 end do 2670 end do 2671 energies%e_elecfield=energies%e_elecfield+eenth 2672 end if ! berryopt==17 2673 2674 ! Turn it into a magnetic enthalpy, by adding orbital electronic contribution 2675 energies%e_magfield = zero 2676 ! if (dtset%berryopt == 5 .and. ipositron/=1) then 2677 ! emag = dot_product(mag_cart,dtset%bfield) 2678 ! energies%e_magfield = emag 2679 ! end if 2680 2681 ! Compute total (free)- energy by direct scheme 2682 if (optene==0) then 2683 etotal = energies%e_kinetic + energies%e_hartree + energies%e_xc & 2684 & + energies%e_localpsp + energies%e_corepsp & 2685 & + energies%e_entropy + energies%e_elecfield & 2686 & + energies%e_magfield + energies%e_nucdip & 2687 & + energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_hybcomp_v & 2688 & + energies%e_constrained_dft + energies%e_ewald & 2689 & + energies%e_chempot + energies%e_vdw_dftd 2690 2691 ! +two*energies%e_fock-energies%e_fock0 ! The Fock energy is already included in the non-local one 2692 ! +energies%e_nlpsp_vfock - energies%e_fock0 2693 2694 ! See similar section in m_energies.F90 2695 ! XG 20181025 This gives a variational energy in case of NCPP with all bands occupied - not yet for metals. 2696 if (usepaw==0) etotal = etotal + energies%e_nlpsp_vfock - energies%e_fock0 2697 ! XG 20181025 I was expecting the following to give also a variational energy in case of PAW, but this is not true. 2698 ! if (usepaw==1) etotal = etotal + energies%e_paw + energies%e_nlpsp_vfock - energies%e_fock0 2699 ! XG 20181025 So, the following is giving a non-variational expression ... 2700 if (usepaw==1) etotal = etotal + energies%e_paw + energies%e_fock 2701 2702 end if 2703 2704 ! Compute total (free) energy by double-counting scheme 2705 if (optene==1) then 2706 etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc & 2707 & - energies%e_xcdc + energies%e_corepsp - energies%e_corepspdc- energies%e_fock0 & 2708 & + energies%e_entropy + energies%e_elecfield + energies%e_magfield & 2709 & + energies%e_hybcomp_E0 - energies%e_hybcomp_v0 + energies%e_constrained_dft 2710 etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd 2711 if (usepaw/=0) etotal = etotal + energies%e_pawdc 2712 end if 2713 2714 ! Additional stuff for electron-positron 2715 if (dtset%positron/=0) then 2716 if (ipositron==0) then 2717 energies%e_electronpositron =zero 2718 energies%edc_electronpositron=zero 2719 else 2720 energies%e_electronpositron =electronpositron%e_hartree+electronpositron%e_xc 2721 energies%edc_electronpositron=electronpositron%e_hartree+electronpositron%e_xcdc 2722 if (usepaw==1) then 2723 energies%e_electronpositron =energies%e_electronpositron +electronpositron%e_paw 2724 energies%edc_electronpositron=energies%edc_electronpositron+electronpositron%e_pawdc 2725 end if 2726 end if 2727 if (optene==0) electronpositron%e0=etotal 2728 if (optene==1) electronpositron%e0=etotal-energies%edc_electronpositron 2729 etotal=electronpositron%e0+energies%e0_electronpositron+energies%e_electronpositron 2730 end if 2731 2732 ! Add the extfpmd energy contribution to the internal energy 2733 if(associated(extfpmd)) then 2734 energies%e_extfpmd=extfpmd%e_kinetic 2735 energies%edc_extfpmd=extfpmd%edc_kinetic 2736 if(optene==0) etotal=etotal+energies%e_extfpmd 2737 if(optene==1) etotal=etotal+energies%edc_extfpmd 2738 end if 2739 2740 ! Compute energy residual 2741 deltae=etotal-elast 2742 elast=etotal 2743 end if !optene/=-1 2744 2745 call timab(80,2,tsec) 2746 2747 !------Compute forces----------------------------------------------------- 2748 2749 if (optforces==1) then 2750 2751 ! PAW: add gradients due to Dij derivatives to non-local term 2752 if (usepaw==1) then 2753 ABI_MALLOC(vlocal,(nfft,dtset%nspden)) 2754 do ispden=1,min(dtset%nspden,2) 2755 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vhartr,vlocal,vpsp,vxc) 2756 do ifft=1,nfft 2757 vlocal(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden) 2758 end do 2759 end do 2760 2761 if(dtset%nspden==4)then 2762 do ispden=3,4 2763 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vlocal,vxc) 2764 do ifft=1,nfft 2765 vlocal(ifft,ispden)=vxc(ifft,ispden) 2766 end do 2767 end do 2768 end if 2769 ucvol_=ucvol 2770 #if defined HAVE_BIGDFT 2771 if (dtset%usewvl==1) ucvol_=product(wvl_den%denspot%dpbox%hgrids)*real(product(wvl_den%denspot%dpbox%ndims),dp) 2772 #endif 2773 dimnhat=0;optgr=1;optgr2=0;optstr=0;optstr2=0 2774 comm_grid=mpi_enreg%comm_fft;if(dtset%usewvl==1) comm_grid=mpi_enreg%comm_wvl 2775 call pawgrnl(atindx1,dimnhat,dummy,1,dummy,grnl,gsqcut,mgfft,my_natom, & 2776 & dtset%natom, nattyp,nfft,ngfft,nhat_dum,dummy,dtset%nspden,dtset%nsym,ntypat,optgr,optgr2,optstr,optstr2,& 2777 & pawang,pawfgrtab,pawrhoij,pawtab,ph1d,psps,k0,rprimd,symrec,dtset%typat,ucvol_,vlocal,vxc,xred, & 2778 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,mpi_comm_grid=mpi_enreg%comm_fft,& 2779 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,paral_kgb=mpi_enreg%paral_kgb) 2780 ABI_FREE(vlocal) 2781 end if 2782 2783 apply_residual=(optres==1 .and. dtset%usewvl==0.and.abs(dtset%densfor_pred)>=1 .and. & 2784 & abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5) 2785 2786 ! If residual is a density residual (and forces from residual asked), 2787 ! has to convert it into a potential residual before calling forces routine 2788 if (apply_residual) then 2789 ABI_MALLOC(resid,(nfft,dtset%nspden)) 2790 option=0; if (dtset%densfor_pred<0) option=1 2791 optnc=1;if (dtset%nspden==4.and.(abs(dtset%densfor_pred)==4.or.abs(dtset%densfor_pred)==6)) optnc=2 2792 call nres2vres(dtset,gsqcut,usepaw,kxc,mpi_enreg,my_natom,nfft,ngfft,nhat,& 2793 & nkxc,nvresid,n3xccc,optnc,option,pawang,pawfgrtab,pawrhoij,pawtab,& 2794 & rhor,rprimd,usepaw,resid,xccc3d,xred,vxc) 2795 else 2796 resid => nvresid 2797 end if 2798 call forces(atindx1,diffor,dtefield,dtset,favg,fcart,fock,forold,gred,grchempottn,grcondft,gresid,grewtn,& 2799 & grhf,grnl,grvdw,grxc,gsqcut,indsym,maxfor,mgfft,mpi_enreg,& 2800 & n1xccc,n3xccc,nattyp,nfft,ngfft,ngrvdw,ntypat,pawrad,pawtab,& 2801 & ph1d,psps,rhog,rhor,rprimd,symrec,synlgr,dtset%usefock,resid,vxc,vxctau,wvl,wvl_den,xred,& 2802 & electronpositron=electronpositron) 2803 if (apply_residual) then 2804 ABI_FREE(resid) 2805 end if 2806 2807 ! Returned gred are full symmetrized gradients of Etotal 2808 ! wrt reduced coordinates xred, d(Etotal)/d(xred) 2809 ! Forces are contained in array fcart 2810 2811 else ! if optforces==0 2812 fcart=zero 2813 gred=zero 2814 favg=zero 2815 diffor=zero 2816 gresid=zero 2817 grhf=zero 2818 maxfor=zero 2819 synlgr=zero 2820 end if 2821 2822 call timab(80,2,tsec) 2823 2824 end subroutine etotfor
ABINIT/m_scfcv_core [ Modules ]
NAME
m_scfcv_core
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (XG, GMR, AR, MKV, MT, FJ, MB) 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 ! nvtx related macro definition 23 #include "nvtx_macros.h" 24 25 module m_scfcv_core 26 27 use defs_basis 28 use defs_wvltypes 29 use defs_rectypes 30 use m_xmpi 31 use m_abicore 32 use m_wffile 33 use m_rec 34 use m_ab7_mixing 35 use m_errors 36 use m_efield 37 use mod_prc_memory 38 use m_nctk 39 use m_hdr 40 use m_xcdata 41 use m_cgtools 42 use m_dtfil 43 use m_distribfft 44 use m_extfpmd 45 46 use m_nonlop, only : nonlop_counter 47 use defs_datatypes, only : pseudopotential_type 48 use defs_abitypes, only : MPI_type 49 use m_berryphase_new, only : update_e_field_vars 50 use m_dens, only : constrained_dft_t, constrained_dft_ini, constrained_dft_free 51 use m_time, only : timab 52 use m_fstrings, only : int2char4, sjoin, itoa 53 use m_symtk, only : symmetrize_xred 54 use m_geometry, only : metric 55 use m_fftcore, only : getng, sphereboundary 56 use m_time, only : abi_wtime, sec2str 57 use m_exit, only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in 58 use m_mpinfo, only : destroy_mpi_enreg, iwrite_fftdatar, initmpi_seq, proc_distrb_cycle 59 use m_ioarr, only : fftdatar_write_from_hdr 60 use m_results_gs , only : results_gs_type 61 use m_scf_history, only : scf_history_type, scf_history_init, scf_history_free 62 use m_energies, only : energies_type, energies_init, energies_copy 63 use m_electronpositron, only : electronpositron_type, electronpositron_calctype 64 use m_pawang, only : pawang_type 65 use m_pawrad, only : pawrad_type 66 use m_pawtab, only : pawtab_type,pawtab_get_lsize 67 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags 68 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 69 use m_pawrhoij, only : pawrhoij_type 70 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, & 71 & pawcprj_free, pawcprj_axpby, pawcprj_put, pawcprj_getdim, pawcprj_reorder 72 use m_pawdij, only : pawdij, symdij 73 use m_pawfgr, only : pawfgr_type 74 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags 75 use m_paw_dmft, only : paw_dmft_type 76 use m_paw_nhat, only : nhatgrid,wvl_nhatgrid,pawmknhat 77 use m_paw_tools, only : chkpawovlp 78 use m_paw_denpot, only : pawdenpot 79 use m_paw_occupancies, only : pawmkrhoij 80 use m_paw_correlations, only : setnoccmmp,setrhoijpbe0 81 use m_paw_mkrho, only : pawmkrho 82 use m_paw_uj, only : pawuj_red, macro_uj_type 83 use m_paw_dfpt, only : pawgrnl 84 use m_fock, only : fock_type, fock_init, fock_destroy, fock_ACE_destroy, fock_common_destroy, & 85 fock_BZ_destroy, fock_update_exc, fock_updatecwaveocc 86 use m_gwls_hamiltonian, only : build_vxc 87 #if defined HAVE_BIGDFT 88 use BigDFT_API, only : cprj_clean,cprj_paw_alloc 89 #endif 90 use m_outxml, only : out_resultsgs_XML, out_geometry_XML 91 use m_kg, only : getcut, getmpw, kpgio, getph 92 use m_fft, only : fourdp 93 use m_vtorhorec, only : first_rec, vtorhorec 94 use m_vtorhotf, only : vtorhotf 95 use m_outscfcv, only : outscfcv 96 use m_afterscfloop, only : afterscfloop 97 use m_extraprho, only : extraprho 98 use m_spacepar, only : make_vectornd,setsym 99 use m_newrho, only : newrho 100 use m_newvtr, only : newvtr 101 use m_vtorho, only : vtorho 102 use m_setvtr, only : setvtr 103 use m_mkrho, only : mkrho 104 use m_rhotov, only : rhotov 105 use m_forces, only : fresid, forces 106 use m_dft_energy, only : energy 107 use m_initylmg, only : initylmg 108 use m_rhotoxc, only : rhotoxc 109 use m_drivexc, only : check_kxc 110 use m_odamix, only : odamix 111 use m_common, only : scprqt, prtene 112 use m_fourier_interpol, only : transgrid 113 use m_fock_getghc, only : fock2ACE 114 use m_forstr, only : nres2vres 115 use m_positron, only : setup_positron 116 use m_cgprj, only : ctocprj 117 use m_psolver, only : psolver_rhohxc 118 use m_paw2wvl, only : paw2wvl_ij, wvl_cprjreorder 119 120 #if defined(HAVE_GPU) && defined(HAVE_GPU_MARKERS) 121 use m_nvtx_data 122 #endif 123 124 implicit none 125 126 private
ABINIT/scfcv_core [ Functions ]
NAME
scfcv_core
FUNCTION
Self-consistent-field convergence. Conducts set of passes or overall iterations of preconditioned conjugate gradient algorithm to converge wavefunctions to ground state and optionally to compute forces and energy. This routine is called to compute forces for given atomic positions or else to do non-SCF band structures.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f) cpus= cpu time limit in seconds dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only) dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset | mband=maximum number of bands | mgfft=maximum size of 1D FFTs for the "coarse" grid (see NOTES below) | mkmem =number of k points treated by this node. | mpw=maximum dimensioned size of npw. | natom=number of atoms in cell. | nfft=(effective) number of FFT grid points (for this processor) | for the "coarse" grid (see NOTES below) | nkpt=number of k points | nspden=number of spin-density components | nsppol=1 for unpolarized, 2 for spin-polarized | nsym=number of symmetry elements in space group ecore=core psp energy (part of total energy) (hartree) fatvshift=factor to multiply dtset%atvshift itimes(2)=itime array, contain itime=itimes(1) and itimimage_gstate=itimes(2) from outer loops kg(3,mpw*mkmem)=reduced planewave coordinates. mcg=size of wave-functions array (cg) =mpw*my_nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor nattyp(ntypat)= # atoms of each type. ndtpawuj=size of dtpawuj npwarr(nkpt)=number of planewaves in basis at this k point paw_dmft <type(paw_dmft_type)>= paw+dmft related data pawang <type(pawang_type)>=paw angular mesh and related data pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat between k-points (see initberry.f) pwind_alloc = first dimension of pwind pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations (see initberry.f) ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
OUTPUT
resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins
SIDE EFFECTS
cg(2,mcg)=updated wavefunctions; if mkmem>=nkpt, these are kept in a disk file. cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn> dtefield <type(efield_type)> = variables related to Berry phase dtpawuj(ndtpawuj)= data used for the automatic determination of U (relevant only for PAW+U) calculations (see initberry.f) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation hdr <type(hdr_type)>=the header of wf, den and pot files indsym(4,nsym,natom)=indirect indexing array for atom labels initialized= if 0 the initialization of the gstate run is not yet finished irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data nfftf=(effective) number of FFT grid points (for this processor) for the "fine" grid (see NOTES below) occ(mband*nkpt*nsppol)=occupation number for each band (often 2) at each k point pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases results_gs <type(results_gs_type)>=results (energy and its components, forces and its components, the stress tensor) of a ground-state computation (should be made a pure output quantity) rhog(2,nfftf)=array for Fourier transform of electron density rhor(nfftf,nspden)=array for electron density in el./bohr**3 rprimd(3,3)=dimensional primitive translations in real space (bohr) scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles symrec(3,3,nsym)=symmetry operations in reciprocal space taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density wffnew=struct info for wf disk files wvl <type(wvl_data)>=all wavelets data xred(3,natom)=reduced dimensionless atomic coordinates xred_old(3,natom)= at input, previous reduced dimensionless atomic coordinates at output, current xred is transferred to xred_old conv_retcode=return code, 0 if convergence was achieved
NOTES
It is worth to explain THE USE OF FFT GRIDS:
============================================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut)
for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ...
are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg)
for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ...
Total density, potentials, ...
are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used.
It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf)
are set equal to (nfft,ngfft,mgfft) in that case.
SOURCE
259 subroutine scfcv_core(atindx,atindx1,cg,cprj,cpus,dmatpawu,dtefield,dtfil,dtpawuj,& 260 & dtset,ecore,eigen,electronpositron,fatvshift,hdr,extfpmd,indsym,& 261 & initialized,irrzon,itimes,kg,mcg,mcprj,mpi_enreg,my_natom,nattyp,ndtpawuj,nfftf,npwarr,occ,& 262 & paw_dmft,pawang,pawfgr,pawrad,pawrhoij,pawtab,phnons,psps,pwind,& 263 & pwind_alloc,pwnsfac,rec_set,resid,results_gs,rhog,rhor,rprimd,& 264 & scf_history,symrec,taug,taur,wffnew,wvl,xred,xred_old,ylm,ylmgr,conv_retcode) 265 266 !Arguments ------------------------------------ 267 !scalars 268 integer,intent(in) :: mcg,my_natom,ndtpawuj,pwind_alloc 269 integer,intent(inout) :: initialized,nfftf,mcprj 270 integer,intent(out) :: conv_retcode 271 real(dp),intent(in) :: cpus,ecore,fatvshift 272 type(MPI_type),intent(inout) :: mpi_enreg 273 type(datafiles_type),intent(in) :: dtfil 274 type(dataset_type),intent(inout) :: dtset 275 type(efield_type),intent(inout) :: dtefield 276 type(electronpositron_type),pointer:: electronpositron 277 type(hdr_type),intent(inout) :: hdr 278 type(extfpmd_type),pointer,intent(inout) :: extfpmd 279 type(pawang_type),intent(in) :: pawang 280 type(pawfgr_type),intent(inout) :: pawfgr 281 type(pseudopotential_type),intent(in) :: psps 282 type(recursion_type),intent(inout) :: rec_set 283 type(results_gs_type),intent(inout) :: results_gs 284 type(scf_history_type),intent(inout) :: scf_history 285 type(wffile_type),intent(inout) :: wffnew 286 type(wvl_data),intent(inout) :: wvl 287 !arrays 288 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom) 289 integer,intent(in) :: indsym(4,dtset%nsym,dtset%natom),itimes(2) 290 !no_abirules 291 integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 292 !(nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise) 293 integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem) 294 integer, intent(in) :: nattyp(psps%ntypat),npwarr(dtset%nkpt),pwind(pwind_alloc,2,3) 295 integer, intent(in) :: symrec(3,3,dtset%nsym) 296 real(dp), intent(inout) :: cg(2,mcg),dmatpawu(:,:,:,:) 297 real(dp), intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 298 real(dp), intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 299 real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 300 !(nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise) 301 real(dp), intent(in) :: pwnsfac(2,pwind_alloc) 302 real(dp), intent(inout) :: rprimd(3,3) 303 real(dp), pointer :: rhog(:,:),rhor(:,:) 304 real(dp), pointer :: taug(:,:),taur(:,:) 305 real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol) 306 real(dp), intent(inout) :: xred(3,dtset%natom) 307 real(dp), intent(inout) :: xred_old(3,dtset%natom) 308 real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 309 real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 310 type(macro_uj_type),intent(inout) :: dtpawuj(0:ndtpawuj) 311 type(pawrhoij_type), intent(inout) :: pawrhoij(my_natom*psps%usepaw) 312 type(pawrad_type), intent(in) :: pawrad(psps%ntypat*psps%usepaw) 313 type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 314 type(paw_dmft_type), intent(inout) :: paw_dmft 315 type(pawcprj_type),pointer, intent(inout) :: cprj(:,:) 316 !Local variables ------------------------- 317 !scalars 318 integer,parameter :: level=110,response=0,cplex1=1 319 integer :: afford,bantot,choice 320 integer :: computed_forces,cplex,cplex_hf,ctocprj_choice,dbl_nnsclo,dielop,dielstrt,dimdmat 321 integer :: forces_needed,errid,has_vxctau,has_dijhat,has_dijnd,has_dijU,has_vhartree,has_dijfock 322 integer :: history_size 323 integer :: iatom,ider,idir,ierr,ii,ikpt,impose_dmat,denpot 324 integer :: initialized0,iorder_cprj,ipert,ipositron,isave_den,isave_kden,iscf10,ispden 325 integer :: ispmix,istep,istep_fock_outer,istep_mix,istep_updatedfock,itypat,izero,lmax_diel,lpawumax,mband_cprj 326 #if defined HAVE_BIGDFT 327 integer :: mcprj_wvl 328 #endif 329 integer :: me,me_wvl,mgfftdiel,mgfftf,moved_atm_inside,moved_rhor,my_nspinor,n1xccc 330 integer :: n3xccc,ncpgr,nfftdiel,nfftmix,nfftmix_per_nfft,nfftotf,ngrcondft,ngrvdw,nhatgrdim,nk3xc,nkxc 331 integer :: npawmix,npwdiel,nremit,nstep,nzlmopt,optcut,optcut_hf,optene,optgr0,optgr0_hf 332 integer :: optgr1,optgr2,optgr1_hf,optgr2_hf,option,optrad,optrad_hf,optres,optxc,prtfor,prtxml,quit 333 integer :: quit_sum,rdwrpaw,shft,spaceComm,spaceComm_fft,spaceComm_wvl,spaceComm_grid 334 integer :: spare_mem 335 integer :: stress_needed,sz1,sz2,tim_mkrho,unit_out 336 integer :: usecprj,usexcnhat,use_hybcomp 337 integer :: my_quit,quitsum_request,timelimit_exit,usecg,wfmixalg,with_vectornd 338 integer ABI_ASYNC :: quitsum_async 339 real(dp) :: boxcut,compch_fft,compch_sph,deltae,diecut,diffor,ecut 340 real(dp) :: ecutf,ecutsus,edum,elast,etotal,evxc,fermie,fermih,gsqcut,hyb_mixing,hyb_mixing_sr 341 real(dp) :: maxfor,res2,residm,ucvol,ucvol_local,val_max 342 real(dp) :: val_min,vxcavg,vxcavg_dum 343 real(dp) :: zion,wtime_step,now,prev,esum,enonlocalpsp !MRM 344 character(len=10) :: tag 345 character(len=500) :: MY_NAME = "scfcv_core" 346 character(len=1500) :: msg 347 !character(len=500) :: dilatmx_errmsg 348 character(len=fnlen) :: fildata 349 type(MPI_type) :: mpi_enreg_diel 350 type(xcdata_type) :: xcdata 351 type(energies_type) :: energies 352 type(ab7_mixing_object) :: mix,mix_mgga 353 logical,parameter :: VERBOSE=.FALSE. 354 logical :: dummy_nhatgr 355 logical :: finite_efield_flag=.false. 356 logical :: non_magnetic_xc=.false. 357 logical :: recompute_cprj=.false.,reset_mixing=.false. 358 logical,save :: tfw_activated=.false. 359 logical :: wvlbigdft=.false. 360 !type(energies_type),pointer :: energies_wvl ! TO BE ACTIVATED LATER 361 !arrays 362 integer :: ngfft(18),ngfftdiel(18),ngfftf(18),ngfftmix(18),npwarr_diel(1) 363 integer :: npwtot_diel(1) 364 integer, save :: scfcv_jdtset = 0 ! To simulate iapp behavior 365 integer, save :: scfcv_itime = 1 ! To simulate iapp behavior 366 integer,allocatable :: dimcprj(:),dimcprj_srt(:) 367 integer,allocatable :: gbound_diel(:,:),irrzondiel(:,:,:),kg_diel(:,:) 368 integer,allocatable :: l_size_atm(:) 369 integer,allocatable :: indsym_dum(:,:,:),symrec_dum(:,:,:), rmm_diis_status(:,:,:) 370 logical,pointer :: lmselect_ep(:,:) 371 real(dp) :: dielar(7),dphase(3),dummy2(6),favg(3),gmet(3,3),gprimd(3,3) 372 real(dp) :: kpt_diel(3),pel(3),pel_cg(3),pelev(3),pion(3),ptot(3),qpt(3),red_ptot(3) !!REC 373 real(dp) :: rhodum(1),rmet(3,3),strscondft(6),strsxc(6),strten(6),tollist(12) 374 real(dp) :: tsec(2),vnew_mean(dtset%nspden),vres_mean(dtset%nspden) 375 real(dp) :: efield_old_cart(3), ptot_cart(3) 376 real(dp) :: red_efield2(3),red_efield2_old(3) 377 real(dp) :: vpotzero(2) 378 ! red_efield1(3),red_efield2(3) is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) [[cite:Stengel2009]] 379 ! red_efield1(3) for fixed ebar calculation, red_efield2(3) for fixed reduced d calculation, in mixed BC 380 ! red_efieldbar_lc(3) is local reduced electric field, defined by Eq.(28) of Nat. Phys. suppl. (2009) [[cite:Stengel2009]] 381 ! pbar(3) and dbar(3) are reduced polarization and displacement field, 382 ! defined by Eq.(27) and (29) Nat. Phys. suppl. (2009) [[cite:Stengel2009]] 383 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 384 real(dp),allocatable :: dielinv(:,:,:,:,:),dtn_pc(:,:) 385 real(dp),allocatable :: fcart(:,:),forold(:,:),gred(:,:),gresid(:,:) 386 real(dp),allocatable :: grchempottn(:,:),grcondft(:,:),grewtn(:,:) 387 real(dp),allocatable :: grhf(:,:),grnl(:),grvdw(:,:),grxc(:,:) 388 real(dp),allocatable :: intgres(:,:),kxc(:,:),nhat(:,:),nhatgr(:,:,:),nvresid(:,:),nvtauresid(:,:) 389 real(dp),allocatable :: ph1d(:,:),ph1ddiel(:,:),ph1df(:,:) 390 real(dp),allocatable :: phnonsdiel(:,:,:),rhowfg(:,:),rhowfr(:,:),shiftvector(:) 391 real(dp),allocatable :: susmat(:,:,:,:,:),synlgr(:,:) 392 real(dp),allocatable :: vectornd(:,:,:),vhartr(:),vpsp(:),vtrial(:,:) 393 real(dp),allocatable :: vxc(:,:),vxc_hybcomp(:,:),vxctau(:,:,:),workr(:,:),xccc3d(:),xcctau3d(:),ylmdiel(:,:) 394 real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:) 395 type(scf_history_type) :: scf_history_wf 396 type(constrained_dft_t) :: constrained_dft 397 type(paw_an_type),allocatable :: paw_an(:) 398 type(paw_ij_type),allocatable :: paw_ij(:) 399 type(pawfgrtab_type),allocatable,save :: pawfgrtab(:) 400 type(pawrhoij_type),pointer :: pawrhoij_ep(:) 401 type(fock_type),pointer :: fock 402 type(pawcprj_type),allocatable, target :: cprj_local(:,:) 403 404 ! ********************************************************************* 405 406 !write(std_out,'(a,5i4)')' scfcv_core, enter : itimes(1:2)=',itimes(1:2) 407 DBG_ENTER("COLL") 408 409 call timab(1440,1,tsec) 410 call timab(1441,3,tsec) 411 412 ! enable time limit handler if not done in callers. 413 if (enable_timelimit_in(MY_NAME) == MY_NAME) then 414 write(std_out,*)"Enabling timelimit check in function: ",trim(MY_NAME)," with timelimit: ",trim(sec2str(get_timelimit())) 415 end if 416 417 ! Initialise non_magnetic_xc for rhohxc 418 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 419 420 !###################################################################### 421 !Initializations - Memory allocations 422 !---------------------------------------------------------------------- 423 lmax_diel = 0 424 425 !MPI communicators 426 if (xmpi_paral==1.and.mpi_enreg%paral_hf==1) then 427 spaceComm=mpi_enreg%comm_kpt 428 else 429 spaceComm=mpi_enreg%comm_cell 430 end if 431 me=xmpi_comm_rank(spaceComm) 432 spaceComm_fft=mpi_enreg%comm_fft 433 spaceComm_wvl=mpi_enreg%comm_wvl 434 me_wvl=mpi_enreg%me_wvl 435 spaceComm_grid=mpi_enreg%comm_fft 436 if(dtset%usewvl==1) spaceComm_grid=mpi_enreg%comm_wvl 437 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 438 439 !Save some variables from dataset definition 440 nstep=dtset%nstep 441 ecut=dtset%ecut 442 ecutf=ecut; if (psps%usepaw==1) ecutf=dtset%pawecutdg 443 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) ecutf=dtset%pawecutdg 444 iscf10=mod(dtset%iscf,10) 445 tollist(1)=dtset%tolmxf;tollist(2)=dtset%tolwfr 446 tollist(3)=dtset%toldff;tollist(4)=dtset%toldfe 447 tollist(6)=dtset%tolvrs;tollist(7)=dtset%tolrff 448 tollist(8)=dtset%vdw_df_threshold 449 dielstrt=0 450 finite_efield_flag=(dtset%berryopt == 4 .or. & 451 & dtset%berryopt == 6 .or. & 452 & dtset%berryopt == 7 .or. & 453 & dtset%berryopt == 14 .or. & 454 & dtset%berryopt == 16 .or. & 455 & dtset%berryopt == 17) 456 457 !Get FFT grid(s) sizes (be careful !) 458 !See NOTES in the comments at the beginning of this file. 459 ngfft(:)=dtset%ngfft(:) 460 if (psps%usepaw==1) then 461 mgfftf=pawfgr%mgfft;ngfftf(:)=pawfgr%ngfft(:) 462 else 463 mgfftf=dtset%mgfft;ngfftf(:)=ngfft(:) 464 end if 465 466 !Calculate zion: the total positive charge acting on the valence electrons 467 zion=zero 468 do iatom=1,dtset%natom 469 zion=zion+psps%ziontypat(dtset%typat(iatom)) 470 end do 471 472 !Compute different geometric tensor, as well as ucvol, from rprimd 473 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 474 475 !Fock: be sure that the pointer is initialized to Null. 476 nullify(fock) 477 478 !Special care in case of WVL 479 !wvlbigdft indicates that the BigDFT workflow will be followed 480 wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1) 481 !if (wvlbigdft) then ! TO BE ACTIVATED LATER 482 ! ABI_MALLOC(energies_wvl,) 483 !end if 484 ucvol_local = ucvol 485 #if defined HAVE_BIGDFT 486 if (dtset%usewvl == 1) then 487 ! We need to tune the volume when wavelets are used because, not 488 ! all FFT points are used. 489 ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3) 490 ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp) 491 end if 492 #endif 493 494 !Some variables need to be initialized/nullified at start 495 nullify(grhor,lrhor,elfr) 496 quit=0 ; dbl_nnsclo=0 ; conv_retcode=0 497 dielop=0 ; strsxc=zero 498 deltae=zero ; elast=zero ; 499 vpotzero(:)=zero 500 ! JWZ April 12 2018: Intel 18 compiler seems to require maxfor initialized, 501 ! else it dies in scprqt in some scenarios 502 maxfor=zero 503 ! 504 results_gs%residm=zero;results_gs%res2=zero 505 results_gs%deltae=zero;results_gs%diffor=zero 506 call energies_init(energies) 507 if (dtset%positron/=0.and.initialized/=0) then 508 energies%e0_electronpositron =results_gs%energies%e0_electronpositron 509 energies%e_electronpositron =results_gs%energies%e_electronpositron 510 energies%edc_electronpositron=results_gs%energies%edc_electronpositron 511 maxfor=zero 512 end if 513 514 ! Initialize fermi level. 515 if (dtset%nstep==0 .or. dtset%iscf < 0) then 516 !if ((dtset%nstep==0 .or. dtset%iscf < 0) .and. dtset%plowan_compute==0) then 517 energies%e_fermie = results_gs%energies%e_fermie 518 results_gs%fermie = results_gs%energies%e_fermie 519 energies%e_fermih = results_gs%energies%e_fermih 520 results_gs%fermih = results_gs%energies%e_fermih 521 ! End CP addition 522 end if 523 524 select case(dtset%usepotzero) 525 case(0,1) 526 energies%e_corepsp = ecore / ucvol 527 energies%e_corepspdc = zero 528 case(2) 529 ! No need to include the PspCore energy since it is already included in the 530 ! local pseudopotential (vpsp) 531 energies%e_corepsp = zero 532 energies%e_corepspdc = zero 533 end select 534 if(wvlbigdft) energies%e_corepsp = zero 535 if(dtset%icutcoul.ne.3) energies%e_corepsp = zero 536 537 538 fermie=energies%e_fermie 539 fermih=energies%e_fermih 540 isave_den=0; isave_kden=0 !initial index of density protection file 541 optres=merge(0,1,dtset%iscf<10) 542 usexcnhat=0!;mcprj=0 543 initialized0=initialized 544 if (dtset%tfkinfunc==12) tfw_activated=.true. 545 ipert=0;idir=0;cplex=1 546 istep_mix=1 547 istep_fock_outer=1 548 ipositron=electronpositron_calctype(electronpositron) 549 unit_out=0;if (dtset%prtvol >= 10) unit_out=ab_out 550 nfftotf=product(ngfftf(1:3)) 551 552 usecprj=0 553 if (mcprj>0) then 554 usecprj=1 555 end if 556 557 !Stresses and forces flags 558 forces_needed=0;prtfor=0 559 if ((dtset%optforces==1.or.dtset%ionmov==4.or.dtset%ionmov==5.or.& 560 & abs(tollist(3))>tiny(0._dp)).or.abs(tollist(7))>tiny(0._dp)) then 561 if (dtset%iscf>0.and.nstep>0) forces_needed=1 562 if (nstep==0) forces_needed=2 563 prtfor=1 564 else if (dtset%iscf>0.and.dtset%optforces==2) then 565 forces_needed=2 566 end if 567 568 stress_needed=0 569 if (dtset%optstress>0.and.dtset%iscf>0.and.dtset%prtstm==0.and. (nstep>0.or.dtfil%ireadwf==1)) stress_needed=1 570 if (dtset%optstress>0.and.dtset%iscf>0.and.psps%usepaw==1 & 571 & .and.finite_efield_flag.and.(nstep>0.or.dtfil%ireadwf==1)) stress_needed=1 572 573 !This is only needed for the tddft routine, and does not 574 !correspond to the intented use of results_gs (should be only 575 !for output of scfcv_core 576 etotal = results_gs%etotal 577 578 !Entering a scfcv_core loop, printing data to XML file if required. 579 prtxml=0;if (me==0.and.dtset%prtxml==1) prtxml=1 580 if (prtxml == 1) then 581 ! scfcv_core() will handle a scf loop, so we output the scfcv markup. 582 write(ab_xml_out, "(A)") ' <scfcvLoop>' 583 write(ab_xml_out, "(A)") ' <initialConditions>' 584 ! We output the geometry of the dataset given in argument. 585 ! xred and rprimd are given independently since dtset only 586 ! stores original and final values. 587 call out_geometry_XML(dtset, 4, dtset%natom, rprimd, xred) 588 write(ab_xml_out, "(A)") ' </initialConditions>' 589 end if 590 591 !Examine tolerance criteria, and eventually print a line to the output 592 !file (with choice=1, the only non-dummy arguments of scprqt are 593 !nstep, tollist and iscf - still, diffor and res2 are here initialized to 0) 594 choice=1 ; diffor=zero ; res2=zero 595 ABI_MALLOC(fcart,(3,dtset%natom)) 596 ABI_MALLOC(gred,(3,dtset%natom)) 597 gred(:,:)=zero 598 fcart(:,:)=results_gs%fcart(:,:) ! This is a side effect ... 599 !results_gs should not be used as input of scfcv_core 600 !HERE IS PRINTED THE FIRST LINE OF SCFCV 601 602 call scprqt(choice,cpus,deltae,diffor,dtset,& 603 & eigen,etotal,favg,fcart,energies%e_fermie,energies%e_fermih,dtfil%fnameabo_app_eig,& 604 & dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,istep_fock_outer,istep_mix,dtset%kptns,& 605 & maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,& 606 & occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,& 607 & psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode) 608 609 !Various allocations (potentials, gradients, ...) 610 ABI_MALLOC(forold,(3,dtset%natom)) 611 ABI_MALLOC(grchempottn,(3,dtset%natom)) 612 ABI_MALLOC(grcondft,(3,dtset%natom)) 613 ABI_MALLOC(gresid,(3,dtset%natom)) 614 ABI_MALLOC(grewtn,(3,dtset%natom)) 615 ABI_MALLOC(grnl,(3*dtset%natom)) 616 ABI_MALLOC(grxc,(3,dtset%natom)) 617 ABI_MALLOC(synlgr,(3,dtset%natom)) 618 ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom)) 619 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom)) 620 ABI_MALLOC(vhartr,(nfftf)) 621 ABI_MALLOC(vtrial,(nfftf,dtset%nspden)) 622 ABI_MALLOC(vpsp,(nfftf)) 623 ABI_MALLOC(vxc,(nfftf,dtset%nspden)) 624 ABI_MALLOC(vxctau,(nfftf,dtset%nspden,4*dtset%usekden)) 625 626 wfmixalg=dtset%fockoptmix/100 627 use_hybcomp=0 628 if(mod(dtset%fockoptmix,100)==11)use_hybcomp=1 629 ABI_MALLOC(vxc_hybcomp,(nfftf,dtset%nspden*use_hybcomp)) 630 631 ngrvdw=0;if (dtset%vdw_xc>=5.and.dtset%vdw_xc<=7) ngrvdw=dtset%natom 632 ABI_MALLOC(grvdw,(3,ngrvdw)) 633 634 ngrcondft=0 635 if(any(dtset%constraint_kind(:)/=0)) ngrcondft=dtset%natom 636 ABI_MALLOC(intgres,(dtset%nspden,ngrcondft)) 637 if(ngrcondft/=0)then 638 intgres(:,:)=zero 639 endif 640 641 grchempottn(:,:)=zero 642 grcondft(:,:)=zero 643 forold(:,:)=zero ; gresid(:,:)=zero ; pel(:)=zero 644 strscondft(:)=zero 645 vtrial(:,:)=zero; vxc(:,:)=zero 646 n1xccc=0;if (psps%n1xccc/=0) n1xccc=psps%n1xccc 647 n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf 648 ABI_MALLOC(xccc3d,(n3xccc)) 649 650 !Allocations/initializations for PAW only 651 lpawumax=-1 652 if(psps%usepaw==1) then 653 ! Variables/arrays related to the fine FFT grid 654 ABI_MALLOC(xcctau3d,(nfftf*dtset%usekden)) 655 ABI_MALLOC(nhat,(nfftf,dtset%nspden*psps%usepaw)) 656 if (nstep==0) nhat=zero 657 ABI_MALLOC(pawfgrtab,(my_natom)) 658 if (my_natom>0) then 659 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,& 660 & mpi_atmtab=mpi_enreg%my_atmtab) 661 call pawfgrtab_init(pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat,& 662 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 663 ABI_FREE(l_size_atm) 664 end if 665 compch_fft=-1.d5 666 usexcnhat=maxval(pawtab(:)%usexcnhat) 667 if (usexcnhat==0.and.dtset%ionmov==4.and.dtset%iscf<10) then 668 ABI_ERROR('You cannot simultaneously use ionmov=4 and such a PAW psp file !') 669 end if 670 671 ! Variables/arrays related to the PAW spheres 672 ABI_MALLOC(paw_ij,(my_natom)) 673 ABI_MALLOC(paw_an,(my_natom)) 674 call paw_an_nullify(paw_an) 675 call paw_ij_nullify(paw_ij) 676 has_dijhat=0;if (dtset%iscf==22) has_dijhat=1 677 has_vhartree=0; if (dtset%prtvha > 0 .or. dtset%prtvclmb > 0) has_vhartree=1 678 has_dijfock=0; if (dtset%usefock == 1) has_dijfock=1 679 has_dijnd=0;if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1 680 has_dijU=merge(0,1,dtset%usepawu>0) !Be careful on this! 681 has_vxctau=dtset%usekden 682 call paw_an_init(paw_an,dtset%natom,dtset%ntypat,0,0,dtset%nspden,& 683 & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,& 684 & has_vxctau=has_vxctau,has_vxc_ex=1,has_vhartree=has_vhartree,& 685 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 686 call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,& 687 & dtset%pawspnorb,dtset%natom,dtset%ntypat,dtset%typat,pawtab,& 688 & has_dij=1,has_dijfock=has_dijfock,has_dijhartree=1,has_dijnd=has_dijnd,has_dijso=1,& 689 & has_dijhat=has_dijhat,& 690 & has_dijU=has_dijU,has_pawu_occ=1,has_exexch_pot=1,nucdipmom=dtset%nucdipmom,& 691 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 692 if(dtset%usewvl==1) then 693 call paw2wvl_ij(1,paw_ij,wvl%descr) 694 end if 695 compch_sph=-1.d5 696 ABI_MALLOC(dimcprj,(dtset%natom)) 697 ABI_MALLOC(dimcprj_srt,(dtset%natom)) 698 call pawcprj_getdim(dimcprj ,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'R') 699 call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 700 do itypat=1,dtset%ntypat 701 if (pawtab(itypat)%usepawu/=0) lpawumax=max(pawtab(itypat)%lpawu,lpawumax) 702 end do 703 if (dtset%usedmatpu/=0.and.lpawumax>0) then 704 if (2*lpawumax+1/=size(dmatpawu,1).or.2*lpawumax+1/=size(dmatpawu,2)) then 705 ABI_BUG('Incorrect size for dmatpawu!') 706 end if 707 end if 708 709 ! Allocation of projected WF (optional) 710 if (usecprj==1) then 711 iorder_cprj=0 712 if (dtset%usefock==1) then 713 ctocprj_choice = 1 714 if (dtset%optforces == 1) then 715 ctocprj_choice = 2; ! ncpgr = 3 716 end if 717 ! if (dtset%optstress /= 0) then 718 ! ncpgr = 6 ; ctocprj_choice = 3 719 ! end if 720 end if 721 722 #if defined HAVE_BIGDFT 723 if (dtset%usewvl==1) then 724 mband_cprj=dtset%mband;if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band 725 mcprj_wvl=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol 726 ABI_MALLOC(wvl%descr%paw%cprj,(dtset%natom,mcprj_wvl)) 727 call cprj_paw_alloc(wvl%descr%paw%cprj,0,dimcprj_srt) 728 end if 729 #endif 730 end if 731 732 ! Other variables for PAW 733 nullify(pawrhoij_ep);if(associated(electronpositron))pawrhoij_ep=>electronpositron%pawrhoij_ep 734 nullify(lmselect_ep);if(associated(electronpositron))lmselect_ep=>electronpositron%lmselect_ep 735 else 736 ABI_MALLOC(dimcprj,(0)) 737 ABI_MALLOC(dimcprj_srt,(0)) 738 ABI_MALLOC(nhat,(0,0)) 739 ABI_MALLOC(xcctau3d,(0)) 740 ABI_MALLOC(paw_ij,(0)) 741 ABI_MALLOC(paw_an,(0)) 742 ABI_MALLOC(pawfgrtab,(0)) 743 end if ! PAW 744 745 !Several parameters and arrays for the SCF mixing: 746 !These arrays are needed only in the self-consistent case 747 if (dtset%iscf>=0) then 748 dielar(1)=dtset%diecut;dielar(2)=dtset%dielng 749 dielar(3)=dtset%diemac;dielar(4)=dtset%diemix 750 dielar(5)=dtset%diegap;dielar(6)=dtset%dielam 751 dielar(7)=dtset%diemix;if (dtset%iscf>=10) dielar(7)=dtset%diemixmag 752 ABI_MALLOC(nvresid,(nfftf,dtset%nspden)) 753 ABI_MALLOC(nvtauresid,(nfftf,dtset%nspden*dtset%usekden)) 754 if (nstep==0) then 755 nvresid=zero 756 nvtauresid=zero 757 end if 758 ABI_MALLOC(dtn_pc,(3,dtset%natom)) 759 ! The next arrays are needed if iscf==5 and ionmov==4, 760 ! but for the time being, they are always allocated 761 ABI_MALLOC(grhf,(3,dtset%natom)) 762 ! Additional allocation for mixing within PAW 763 npawmix=0 764 if(psps%usepaw==1) then 765 do iatom=1,my_natom 766 itypat=pawrhoij(iatom)%itypat 767 pawrhoij(iatom)%use_rhoijres=1 768 sz1=pawrhoij(iatom)%cplex_rhoij*pawtab(itypat)%lmn2_size 769 sz2=pawrhoij(iatom)%nspden 770 ABI_MALLOC(pawrhoij(iatom)%rhoijres,(sz1,sz2)) 771 do ispden=1,pawrhoij(iatom)%nspden 772 pawrhoij(iatom)%rhoijres(:,ispden)=zero 773 end do 774 ABI_MALLOC(pawrhoij(iatom)%kpawmix,(pawtab(itypat)%lmnmix_sz)) 775 pawrhoij(iatom)%lmnmix_sz=pawtab(itypat)%lmnmix_sz 776 pawrhoij(iatom)%kpawmix=pawtab(itypat)%kmix 777 npawmix=npawmix+pawrhoij(iatom)%nspden*pawtab(itypat)%lmnmix_sz & 778 & *pawrhoij(iatom)%cplex_rhoij*pawrhoij(iatom)%qphase 779 end do 780 end if 781 if (dtset%iscf > 0) then 782 denpot = AB7_MIXING_POTENTIAL 783 if (dtset%iscf > 10) denpot = AB7_MIXING_DENSITY 784 if (psps%usepaw==1.and.dtset%pawmixdg==0 .and. dtset%usewvl==0) then 785 ispmix=AB7_MIXING_FOURRIER_SPACE;nfftmix=dtset%nfft;ngfftmix(:)=ngfft(:) 786 else 787 ispmix=AB7_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:) 788 end if 789 !TRangel: added to avoid segfaults with Wavelets 790 nfftmix_per_nfft=0;if(nfftf>0) nfftmix_per_nfft=(1-nfftmix/nfftf) 791 call ab7_mixing_new(mix, iscf10, denpot, ispmix, nfftmix, dtset%nspden, npawmix, errid, msg, dtset%npulayit) 792 if (errid /= AB7_NO_ERROR) then 793 ABI_ERROR(msg) 794 end if 795 if (dtset%usekden/=0) then 796 if (dtset%useria==12345) then ! This is temporary 797 call ab7_mixing_new(mix_mgga, iscf10, denpot, ispmix, nfftmix, dtset%nspden, 0, errid, msg, dtset%npulayit) 798 else 799 call ab7_mixing_new(mix_mgga, 0, denpot, ispmix, nfftmix, dtset%nspden, 0, errid, msg, dtset%npulayit) 800 end if 801 if (errid /= AB7_NO_ERROR) then 802 ABI_ERROR(msg) 803 end if 804 end if 805 if (dtset%mffmem == 0) then 806 call ab7_mixing_use_disk_cache(mix, dtfil%fnametmp_fft) 807 if (dtset%usekden/=0.and.denpot==AB7_MIXING_DENSITY) & 808 & call ab7_mixing_use_disk_cache(mix, dtfil%fnametmp_fft_mgga) 809 end if 810 ! else if (dtset%iscf==0.and.dtset%usewvl==1) then 811 ! ispmix=AB7_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:) 812 end if 813 else 814 ABI_MALLOC(nvresid,(0,0)) 815 ABI_MALLOC(nvtauresid,(0,0)) 816 ABI_MALLOC(dtn_pc,(0,0)) 817 ABI_MALLOC(grhf,(0,0)) 818 end if ! iscf>0 819 820 ! Here initialize the datastructure constrained_dft, for constrained DFT calculations 821 ! as well as penalty function constrained magnetization 822 if(any(dtset%constraint_kind(:)/=0).or.dtset%magconon/=0)then 823 call constrained_dft_ini(dtset%chrgat,constrained_dft,dtset%constraint_kind,dtset%magconon,dtset%magcon_lambda,& 824 & mpi_enreg,dtset%natom,nfftf,ngfftf,dtset%nspden,dtset%ntypat,& 825 & dtset%ratsm,dtset%ratsph,rprimd,dtset%spinat,dtset%typat,xred,dtset%ziontypat) 826 endif 827 828 !Here, allocate arrays for computation of susceptibility and dielectric matrix or for TDDFT 829 830 if( (nstep>0 .and. dtset%iscf>=0) .or. dtset%iscf==-1 ) then !MF 831 832 ! Here, for TDDFT, artificially set iprcel . Also set a variable to reduce the memory needs. 833 afford=1 834 if(dtset%iscf==-1) then 835 ! dtset%iprcel=21 836 afford=0 837 end if 838 839 ! First compute dimensions 840 if(dtset%iprcel>=21 .or. dtset%iscf==-1)then 841 ! With dielop=1, the matrices will be computed when istep=dielstrt 842 ! With dielop=2, the matrices will be computed when istep=dielstrt and 1 843 dielop=1 844 if(dtset%iprcel>=41)dielop=2 845 if((dtset%iprcel >= 71).and.(dtset%iprcel<=79)) dielop=0 !RSkerker preconditioner do not need the susceptibility matrix 846 ! Immediate computation of dielectric matrix 847 dielstrt=1 848 ! Or delayed computation 849 if(modulo(dtset%iprcel,100)>21 .and. modulo(dtset%iprcel,100)<=29)dielstrt=modulo(dtset%iprcel,100)-20 850 if(modulo(dtset%iprcel,100)>31 .and. modulo(dtset%iprcel,100)<=39)dielstrt=modulo(dtset%iprcel,100)-30 851 if(modulo(dtset%iprcel,100)>41 .and. modulo(dtset%iprcel,100)<=49)dielstrt=modulo(dtset%iprcel,100)-40 852 if(modulo(dtset%iprcel,100)>51 .and. modulo(dtset%iprcel,100)<=59)dielstrt=modulo(dtset%iprcel,100)-50 853 if(modulo(dtset%iprcel,100)>61 .and. modulo(dtset%iprcel,100)<=69)dielstrt=modulo(dtset%iprcel,100)-60 854 ! Get diecut, and the fft grid to be used for the susceptibility computation 855 diecut=abs(dtset%diecut) 856 if( dtset%diecut<0.0_dp )then 857 ecutsus=ecut 858 else 859 ecutsus= ( sqrt(ecut) *0.5_dp + sqrt(diecut) *0.25_dp )**2 860 end if 861 ! Impose sequential calculation 862 ngfftdiel(1:3)=0 ; ngfftdiel(7)=100 ; ngfftdiel(9)=0; ngfftdiel(8)=dtset%ngfft(8);ngfftdiel(10:18)=0 863 if(dtset%iscf==-1)ngfftdiel(7)=102 864 865 ! The dielectric stuff is performed in sequential mode; set mpi_enreg_diel accordingly 866 call initmpi_seq(mpi_enreg_diel) 867 call getng(dtset%boxcutmin,dtset%chksymtnons,ecutsus,gmet,k0,mpi_enreg_diel%me_fft,mgfftdiel,nfftdiel,ngfftdiel,& 868 & mpi_enreg_diel%nproc_fft,dtset%nsym,mpi_enreg_diel%paral_kgb,dtset%symrel,dtset%tnons,& 869 & gpu_option=dtset%gpu_option) 870 ! Update the fft distribution 871 call init_distribfft_seq(mpi_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all') 872 873 ! Compute the size of the dielectric matrix 874 kpt_diel(1:3)=(/ 0.0_dp, 0.0_dp, 0.0_dp /) 875 call getmpw(diecut,dtset%exchn2n3d,gmet,(/1/),kpt_diel,mpi_enreg_diel,npwdiel,1) 876 lmax_diel=0 877 if (psps%usepaw==1) then 878 do ii=1,dtset%ntypat 879 lmax_diel=max(lmax_diel,pawtab(ii)%lcut_size) 880 end do 881 end if 882 else 883 npwdiel=1 884 mgfftdiel=1 885 nfftdiel=1 886 lmax_diel=0 887 afford=0 888 end if 889 890 ! Now, performs allocation 891 ABI_MALLOC(dielinv,(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden)) 892 ABI_MALLOC(susmat,(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden)) 893 ABI_MALLOC(kg_diel,(3,npwdiel)) 894 ABI_MALLOC(gbound_diel,(2*mgfftdiel+8,2)) 895 ABI_MALLOC(irrzondiel,(nfftdiel**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 896 ABI_MALLOC(phnonsdiel,(2,nfftdiel**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 897 ABI_MALLOC(ph1ddiel,(2,3*(2*mgfftdiel+1)*dtset%natom*psps%usepaw)) 898 ABI_MALLOC(ylmdiel,(npwdiel,lmax_diel**2)) 899 ! Then, compute the values of different arrays 900 if(dielop>=1)then 901 ! Note : npwarr_diel is dummy, npwtot_diel is dummy 902 ! This kpgio call for going from the suscep FFT grid to the diel sphere 903 npwarr_diel(1)=npwdiel 904 905 call kpgio(diecut,dtset%exchn2n3d,gmet,(/1/),kg_diel,& 906 & kpt_diel,1,(/1/),1,'COLL',mpi_enreg_diel,npwdiel,& 907 & npwarr_diel,npwtot_diel,dtset%nsppol) 908 call sphereboundary(gbound_diel,1,kg_diel,mgfftdiel,npwdiel) 909 910 if (dtset%nsym>1 .and. dtset%iscf>=0 ) then 911 ! Should replace this initialization of irrzondiel and phnonsdiel through setsym by a direct call to irrzg 912 ABI_MALLOC(indsym_dum,(4,dtset%nsym,dtset%natom)) 913 ABI_MALLOC(symrec_dum,(3,3,dtset%nsym)) 914 call setsym(indsym_dum,irrzondiel,dtset%iscf,dtset%natom,& 915 & nfftdiel,ngfftdiel,dtset%nspden,dtset%nsppol,dtset%nsym,phnonsdiel,& 916 & dtset%symafm,symrec_dum,dtset%symrel,dtset%tnons,dtset%typat,xred) 917 ABI_FREE(indsym_dum) 918 ABI_FREE(symrec_dum) 919 end if 920 if (psps%usepaw==1) then 921 call getph(atindx,dtset%natom,ngfftdiel(1),ngfftdiel(2),& 922 & ngfftdiel(3),ph1ddiel,xred) 923 call initylmg(gprimd,kg_diel,kpt_diel,1,mpi_enreg_diel,& 924 & lmax_diel,npwdiel,dtset%nband,1,npwarr_diel,dtset%nsppol,0,& 925 & rprimd,ylmdiel,rhodum) 926 end if 927 end if 928 929 if(dtset%iprcel>=21 .or. dtset%iscf==-1)then 930 call destroy_mpi_enreg(mpi_enreg_diel) 931 end if 932 933 else 934 npwdiel=1 935 mgfftdiel=1 936 nfftdiel=1 937 afford = 0 938 ABI_MALLOC(susmat,(0,0,0,0,0)) 939 ABI_MALLOC(kg_diel,(0,0)) 940 ABI_MALLOC(gbound_diel,(0,0)) 941 ABI_MALLOC(irrzondiel,(0,0,0)) 942 ABI_MALLOC(phnonsdiel,(0,0,0)) 943 ABI_MALLOC(ph1ddiel,(0,0)) 944 ABI_MALLOC(ylmdiel,(0,0)) 945 end if 946 947 nkxc=0 948 !TDDFT - For a first coding 949 if (dtset%iscf==-1 .and. dtset%nspden==1) nkxc=2 950 if (dtset%iscf==-1 .and. dtset%nspden==2) nkxc=3 951 !Eventually need kxc-LDA when susceptibility matrix has to be computed 952 if (dtset%iscf>0.and.modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79)) nkxc=2*min(dtset%nspden,2)-1 953 !Eventually need kxc-LDA for residual forces (when density mixing is selected) 954 if (dtset%iscf>=10.and.dtset%usewvl==0.and.forces_needed>0 .and. & 955 & abs(dtset%densfor_pred)>=1.and.abs(dtset%densfor_pred)<=6.and.abs(dtset%densfor_pred)/=5) then 956 if (dtset%xclevel==1.or.dtset%densfor_pred>=0) nkxc=2*min(dtset%nspden,2)-1 957 if (dtset%xclevel==2.and.dtset%nspden==1.and.dtset%densfor_pred<0) nkxc=7 ! This is not full kxc for mGGA 958 if (dtset%xclevel==2.and.dtset%nspden==2.and.dtset%densfor_pred<0) nkxc=19 ! This is not full kxc for mGGA 959 end if 960 if (nkxc>0) then 961 call check_kxc(dtset%ixc,dtset%optdriver) 962 end if 963 ABI_MALLOC(kxc,(nfftf,nkxc)) 964 965 !This flag will be set to 1 just before an eventual change of atomic 966 !positions inside the iteration, and set to zero when the consequences 967 !of this change are taken into account. 968 moved_atm_inside=0 969 !This flag will be set to 1 if the forces are computed inside the iteration. 970 computed_forces=0 971 972 if(dtset%wfoptalg==2)then 973 ABI_MALLOC(shiftvector,((dtset%mband+2)*dtset%nkpt)) 974 val_min=-1.0_dp 975 val_max=zero 976 else 977 ABI_MALLOC(shiftvector,(1)) 978 end if 979 980 !!PAW+DMFT: allocate structured datatype paw_dmft if dtset%usedmft=1 981 !call init_sc_dmft(dtset%dmftbandi,dtset%dmftbandf,dtset%mband,dtset%nkpt,& 982 !& dtset%nsppol,dtset%usedmft,paw_dmft,dtset%usedmft) 983 !call print_sc_dmft(paw_dmft) 984 985 !!Electric field initializations: initialize pel_cg(:) and p_ion(:) 986 call update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,& 987 & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,& 988 & dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,ngfft,dtset%nkpt,npwarr,& 989 & dtset%ntypat,pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,& 990 & pwind,pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,& 991 & 0,quit,istep,ucvol,unit_out,psps%usepaw,xred,ylm,ylmgr) 992 993 if (dtset%iscf==22) energies%h0=zero 994 995 call timab(1441,2,tsec) 996 997 !################################################################## 998 !PERFORM ELECTRONIC ITERATIONS 999 !################################################################## 1000 1001 !Offer option of computing total energy with existing 1002 !wavefunctions when nstep<=0, else do nstep iterations 1003 !Note that for non-self-consistent calculations, this loop will be exited 1004 !after the first call to vtorho 1005 !Pass through the first routines even when nstep==0 1006 1007 quitsum_request = xmpi_request_null; timelimit_exit = 0 1008 istep_updatedfock=0 1009 1010 ABI_ICALLOC(rmm_diis_status, (2, dtset%nkpt, dtset%nsppol)) 1011 1012 ! start SCF loop 1013 ABI_NVTX_START_RANGE(NVTX_SCF) 1014 do istep=1,max(1,nstep) 1015 1016 ! Handle time limit condition. 1017 if (istep == 1) prev = abi_wtime() 1018 if (istep > 1) then 1019 now = abi_wtime() 1020 wtime_step = now - prev 1021 prev = now 1022 call wrtout(std_out, sjoin("{SCF_istep:", itoa(istep-1), ", Vnl|psi>:", itoa(nonlop_counter), & 1023 ", wall_time: '", sec2str(wtime_step), "'} <<< TIME")) 1024 nonlop_counter = 0 1025 1026 if (have_timelimit_in(MY_NAME)) then 1027 if (istep > 2) then 1028 call xmpi_wait(quitsum_request,ierr) 1029 if (quitsum_async > 0) then 1030 write(msg,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())),". Will exit istep loop in scfcv_core." 1031 ABI_COMMENT(msg) 1032 call wrtout(ab_out, msg) 1033 timelimit_exit = 1 1034 exit 1035 end if 1036 end if 1037 1038 my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1 1039 call xmpi_isum(my_quit,quitsum_async,spacecomm,quitsum_request,ierr) 1040 end if 1041 end if 1042 1043 call timab(1442,1,tsec) 1044 if (moved_atm_inside==1 .or. istep==1) then 1045 ! ############################################################## 1046 ! The following steps are done once for a given set of atomic 1047 ! coordinates or for the nstep=1 case 1048 ! -------------------------------------------------------------- 1049 1050 ! Eventually symmetrize atomic coordinates over space group elements: 1051 call symmetrize_xred(dtset%natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym) 1052 1053 if (dtset%usewvl == 0) then 1054 ! Get cut-off for g-vectors 1055 if (psps%usepaw==1) call wrtout(std_out,' FFT (fine) grid used in SCF cycle:') 1056 call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf) 1057 1058 ! Compute structure factor phases and large sphere cut-off (gsqcut): 1059 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 1060 1061 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 1062 call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred) 1063 else 1064 ph1df(:,:)=ph1d(:,:) 1065 end if 1066 end if 1067 1068 ! Initialization of atomic data for PAW 1069 if (psps%usepaw==1) then 1070 1071 ! Check for non-overlapping spheres, allow one remittance in case of optimization or MD or image algorithms 1072 nremit=0 1073 if(dtset%ionmov>0 .and. itimes(1)==1)nremit=1 1074 if(dtset%imgmov>0 .and. itimes(2)==1)nremit=mpi_enreg%my_nimage 1075 call chkpawovlp(dtset%natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred,nremit=nremit) 1076 1077 ! Identify parts of the rectangular grid where the density has to be calculated 1078 optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm 1079 if ((forces_needed==1).or. & 1080 & (dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0).or. & 1081 & (dtset%positron/=0.and.forces_needed==2)) then 1082 optgr1=dtset%pawstgylm;if (stress_needed==1) optrad=1; if (dtset%pawprtwf==1) optrad=1 1083 end if 1084 1085 if(dtset%usewvl==0) then 1086 call nhatgrid(atindx1,gmet,my_natom,dtset%natom,& 1087 & nattyp,ngfftf,psps%ntypat,optcut,optgr0,optgr1,optgr2,optrad,& 1088 & pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,& 1089 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 1090 & comm_fft=spaceComm_fft,distribfft=mpi_enreg%distribfft) 1091 else 1092 shft=0 1093 #if defined HAVE_BIGDFT 1094 shft=wvl%descr%Glr%d%n1i*wvl%descr%Glr%d%n2i*wvl%den%denspot%dpbox%nscatterarr(me_wvl,4) 1095 call wvl_nhatgrid(atindx1,wvl%descr%atoms%astruct%geocode,& 1096 & wvl%descr%h,wvl%den%denspot%dpbox%i3s,dtset%natom,dtset%natom,& 1097 & nattyp,psps%ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,& 1098 & wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,& 1099 & wvl%den%denspot%dpbox%n3pi,optcut,optgr0,optgr1,optgr2,optrad,& 1100 & pawfgrtab,pawtab,psps%gth_params%psppar,rprimd,shft,xred) 1101 #endif 1102 end if 1103 end if 1104 1105 ! If we are inside SCF cycle or inside dynamics over ions, 1106 ! we have to translate the density of previous iteration 1107 moved_rhor=0 1108 1109 if (initialized/=0.and.dtset%usewvl == 0.and.ipositron/=1.and. & 1110 & (abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)) then 1111 moved_rhor=1 1112 if (abs(dtset%densfor_pred)==2) then 1113 option=2 1114 ABI_MALLOC(workr,(nfftf,dtset%nspden)) 1115 call fresid(dtset,gresid,mpi_enreg,nfftf,ngfftf,& 1116 & psps%ntypat,option,pawtab,rhor,rprimd,& 1117 & ucvol,workr,xred,xred_old,psps%znuclpsp) 1118 rhor=workr 1119 ABI_FREE(workr) 1120 else if (abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6) then 1121 scf_history%icall=scf_history%icall+1 1122 call extraprho(atindx,atindx1,cg,cprj,dtset,gmet,gprimd,gsqcut,& 1123 & scf_history%icall,kg,mcg,mcprj,mgfftf,mpi_enreg,psps%mqgrid_vl,& 1124 & my_natom,nattyp,nfftf,ngfftf,npwarr,psps%ntypat,pawrhoij,pawtab,& 1125 & ph1df,psps,psps%qgrid_vl,rhor,rprimd,scf_history,ucvol,& 1126 & psps%usepaw,xred,xred_old,ylm,psps%ziontypat,psps%znuclpsp) 1127 end if 1128 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 1129 end if 1130 if (initialized/=0.and.dtset%usewvl == 0.and.ipositron/=1) then 1131 ! In some cases cprj are kept in memory, so we have to update them before the call of vtorho 1132 if (dtset%cprj_in_memory/=0) then 1133 iatom=0 1134 idir=0 1135 iorder_cprj=0 1136 call wrtout(std_out,' Computing cprj from wavefunctions (scfcv_core)') 1137 ABI_NVTX_START_RANGE(NVTX_CTOCPRJ) 1138 call ctocprj(atindx,cg,1,cprj,gmet,gprimd,iatom,idir,& 1139 & iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,& 1140 & dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft, dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,& 1141 & dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,& 1142 & xred,ylm,ylmgr) 1143 ABI_NVTX_END_RANGE() 1144 call wrtout(std_out,' cprj is computed') 1145 end if 1146 end if 1147 1148 ! if any nuclear dipoles are nonzero, compute the vector potential in real space (depends on 1149 ! atomic position so should be done for nstep = 1 and for updated ion positions 1150 if ( any(abs(dtset%nucdipmom(:,:))>tol8) ) then 1151 with_vectornd = 1 1152 else 1153 with_vectornd = 0 1154 end if 1155 if(allocated(vectornd)) then 1156 ABI_FREE(vectornd) 1157 end if 1158 ABI_MALLOC(vectornd,(with_vectornd*nfftf,dtset%nspden,3)) 1159 vectornd=zero 1160 if(with_vectornd .EQ. 1) then 1161 call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,ngfftf,& 1162 & dtset%nspden,dtset%nucdipmom,rprimd,vectornd,xred) 1163 endif 1164 1165 end if ! moved_atm_inside==1 .or. istep==1 1166 1167 call timab(1442,2,tsec) 1168 1169 !Initialize/Update data in the case of an Exact-exchange (Hartree-Fock) or hybrid XC calculation 1170 hyb_mixing=zero;hyb_mixing_sr=zero 1171 if (dtset%usefock==1) then 1172 call timab(1443,1,tsec) 1173 if (istep==1) then 1174 ! Initialize data_type fock for the calculation 1175 cplex_hf=cplex; if (psps%usepaw==1) cplex_hf=dtset%pawcpxocc 1176 call fock_init(atindx,cplex_hf,dtset,fock,gsqcut,kg,mpi_enreg,nattyp,npwarr,pawang,pawfgr,pawtab,rprimd) 1177 1178 if (fock%fock_common%usepaw==1) then 1179 optcut_hf = 0 ! use rpaw to construct local_pawfgrtab 1180 optgr0_hf = 0; optgr1_hf = 0; optgr2_hf = 0 ! dont need gY terms locally 1181 optrad_hf = 1 ! do store r-R 1182 call nhatgrid(atindx1,gmet,dtset%natom,dtset%natom,nattyp,ngfftf,psps%ntypat,& 1183 & optcut_hf,optgr0_hf,optgr1_hf,optgr2_hf,optrad_hf,fock%fock_common%pawfgrtab,pawtab,& 1184 & rprimd,dtset%typat,ucvol,xred,typord=1) 1185 iatom=-1;idir=0 1186 call ctocprj(atindx,cg,ctocprj_choice,cprj,gmet,gprimd,iatom,idir,& 1187 & iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,& 1188 & dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft, dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,& 1189 & dtset%nsppol,dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,& 1190 & xred,ylm,ylmgr) 1191 end if 1192 if(wfmixalg/=0)then 1193 spare_mem=0 1194 if(spare_mem==1)history_size=wfmixalg ! Not yet coded 1195 if(spare_mem==0)history_size=2*(wfmixalg-1)+1 1196 ! Specific case of simple mixing : always history_size=1 1197 if(wfmixalg==2)history_size=1 1198 scf_history_wf%history_size=history_size 1199 usecg=2 1200 call scf_history_init(dtset,mpi_enreg,usecg,scf_history_wf) 1201 end if 1202 end if 1203 1204 !Fock energy 1205 energies%e_exactX=zero 1206 if (fock%fock_common%optfor) fock%fock_common%forces=zero 1207 1208 call timab(1443,2,tsec) 1209 1210 if (istep==1 .or. istep_updatedfock==fock%fock_common%nnsclo_hf .or. & 1211 & (fock%fock_common%nnsclo_hf>1 .and. fock%fock_common%scf_converged) ) then 1212 1213 istep_updatedfock=1 1214 fock%fock_common%scf_converged=.false. 1215 1216 call timab(1444,1,tsec) 1217 1218 !Possibly mix the wavefunctions from different steps before computing the Fock operator 1219 if(wfmixalg/=0 .and. .not. (wfmixalg==2 .and. abs(scf_history_wf%alpha-one)<tol8) )then 1220 call wf_mixing(atindx1,cg,cprj,dtset,istep_fock_outer,mcg,mcprj,mpi_enreg,& 1221 & nattyp,npwarr,pawtab,scf_history_wf) 1222 istep_fock_outer=istep_fock_outer+1 1223 1224 call timab(1444,2,tsec) 1225 1226 !DEBUG 1227 if(.false.)then 1228 !Update the density, from the newly mixed cg and cprj. 1229 !Be careful: in PAW, rho does not include the compensation density (added later) ! 1230 tim_mkrho=6 1231 if (psps%usepaw==1) then 1232 ABI_MALLOC(rhowfg,(2,dtset%nfft)) 1233 ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden)) 1234 ! 1-Compute density from WFs 1235 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,& 1236 & rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,extfpmd=extfpmd) 1237 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor) 1238 ! 2-Compute rhoij 1239 call pawmkrhoij(atindx,atindx1,cprj,dimcprj,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,& 1240 & mcprj,dtset%mkmem,mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspden,dtset%nspinor,& 1241 & dtset%nsppol,occ,dtset%paral_kgb,paw_dmft,pawrhoij,dtfil%unpaw,dtset%usewvl,dtset%wtk) 1242 ! 3-Symetrize rhoij, compute nhat and add it to rhor 1243 ! Note pawrhoij_unsym and pawrhoij are the same, which means that pawrhoij 1244 ! cannot be distributed over different atomic sites. 1245 cplex=1;ipert=0;idir=0;qpt(:)=zero 1246 call pawmkrho(1,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,& 1247 & my_natom,dtset%natom,dtset%nspden,dtset%nsym,dtset%ntypat,& 1248 & dtset%paral_kgb,pawang,pawfgr,pawfgrtab,dtset%pawprtvol,pawrhoij,pawrhoij,& 1249 & pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,symrec,dtset%typat,ucvol,& 1250 & dtset%usewvl,xred,pawnhat=nhat,rhog=rhog) 1251 ! 2-Take care of kinetic energy density 1252 if(dtset%usekden==1)then 1253 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,& 1254 & rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 1255 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,taug,rhowfr,taur) 1256 end if 1257 ABI_FREE(rhowfg) 1258 ABI_FREE(rhowfr) 1259 else 1260 write(std_out,*)' scfcv_core : recompute the density after the wf mixing ' 1261 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 1262 & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,extfpmd=extfpmd) 1263 if(dtset%usekden==1)then 1264 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,taug,taur,& 1265 & rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 1266 end if 1267 end if 1268 end if 1269 end if 1270 !ENDDEBUG 1271 1272 call timab(1445,1,tsec) 1273 1274 ! Update data relative to the occupied states in fock 1275 call fock_updatecwaveocc(cg,cprj,dtset,fock,indsym,mcg,mcprj,mpi_enreg,nattyp,npwarr,occ,ucvol) 1276 1277 call timab(1445,2,tsec) 1278 call timab(1446,3,tsec) 1279 1280 ! Possibly (re)compute the ACE operator 1281 if(fock%fock_common%use_ACE/=0) then 1282 call fock2ACE(cg,cprj,fock,dtset%istwfk,kg,dtset%kptns,dtset%mband,mcg,mcprj,dtset%mgfft,& 1283 & dtset%mkmem,mpi_enreg,psps%mpsang,& 1284 & dtset%mpw,my_natom,dtset%natom,dtset%nband,dtset%nfft,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspden,& 1285 & dtset%nspinor,dtset%nsppol,dtset%ntypat,occ,dtset%optforces,paw_ij,pawtab,ph1d,psps,rprimd,& 1286 & dtset%typat,usecprj,dtset%gpu_option,dtset%wtk,xred,ylm) 1287 energies%e_fock0=fock%fock_common%e_fock0 1288 end if 1289 1290 call timab(1446,2,tsec) 1291 1292 !Should place a test on whether there should be the final exit of the istep loop. 1293 !This test should use focktoldfe. 1294 !This should update the info in fock%fock_common%fock_converged. 1295 !For the time being, fock%fock_common%fock_converged=.false., so the loop end with the maximal value of nstep always, 1296 !except when nnsclo_hf==1 (so the Fock operator is always updated), in which case, the usual exit tests (toldfe, tolvrs, etc) 1297 !work fine. 1298 !if(fock%fock_common%nnsclo_hf==1 .and. fock%fock_common%use_ACE==0)then 1299 if(fock%fock_common%nnsclo_hf==1) fock%fock_common%fock_converged=.TRUE. 1300 1301 !Depending on fockoptmix, possibly restart the mixing procedure for the potential 1302 if(mod(dtset%fockoptmix,10)==1) istep_mix=1 1303 else 1304 istep_updatedfock=istep_updatedfock+1 1305 end if 1306 1307 !Used locally 1308 hyb_mixing=fock%fock_common%hyb_mixing ; hyb_mixing_sr=fock%fock_common%hyb_mixing_sr 1309 end if ! usefock 1310 1311 call timab(1447,1,tsec) 1312 1313 ! Initialize/update data in the electron-positron case 1314 if (dtset%positron<0.or.(dtset%positron>0.and.istep==1)) then 1315 call setup_positron(atindx,atindx1,cg,cprj,dtefield,dtfil,dtset,ecore,eigen,& 1316 & etotal,electronpositron,energies,fock,forces_needed,gred,gmet,gprimd,& 1317 & grchempottn,grcondft,grewtn,grvdw,gsqcut,hdr,extfpmd,initialized0,indsym,istep,istep_mix,kg,& 1318 & kxc,maxfor,mcg,mcprj,mgfftf,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfftf,ngrvdw,nhat,& 1319 & nkxc,npwarr,nvresid,occ,optres,paw_ij,pawang,pawfgr,pawfgrtab,& 1320 & pawrad,pawrhoij,pawtab,ph1df,ph1d,psps,rhog,rhor,rmet,rprimd,& 1321 & stress_needed,strscondft,strsxc,symrec,ucvol,usecprj,vhartr,vpsp,vxc,vxctau,& 1322 & xccc3d,xcctau3d,xred,ylm,ylmgr) 1323 ipositron=electronpositron_calctype(electronpositron) 1324 end if 1325 1326 call timab(1447,2,tsec) 1327 call timab(1448,3,tsec) 1328 1329 if ((moved_atm_inside==1 .or. istep==1).or.& 1330 & (dtset%positron<0.and.istep_mix==1).or.& 1331 & (mod(dtset%fockoptmix,100)==11 .and. istep_updatedfock==1)) then 1332 ! PAW only: we sometimes have to compute compensation density 1333 ! and eventually add it to density from WFs 1334 nhatgrdim=0 1335 dummy_nhatgr = .False. 1336 if (psps%usepaw==1.and.(dtset%positron>=0.or.ipositron/=1) & 1337 & .and.((usexcnhat==0) & 1338 & .or.(dtset%xclevel==2.and.(dtfil%ireadwf/=0.or.dtfil%ireadden/=0.or.initialized/=0)) & 1339 & .or.(dtfil%ireadwf/=0.and.dtfil%ireadden==0.and.initialized==0))) then 1340 call timab(558,1,tsec) 1341 nhatgrdim=0;if (dtset%xclevel==2) nhatgrdim=usexcnhat*dtset%pawnhatxc 1342 ider=2*nhatgrdim;izero=0 1343 if (nhatgrdim>0) then 1344 ABI_MALLOC(nhatgr,(cplex*nfftf,dtset%nspden,3*nhatgrdim)) 1345 else 1346 ABI_MALLOC(nhatgr,(0,0,0)) 1347 dummy_nhatgr = .True. 1348 end if 1349 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,& 1350 & nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,& 1351 & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,k0,rprimd,ucvol_local,dtset%usewvl,xred,& 1352 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 1353 & comm_fft=spaceComm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 1354 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 1355 if (dtfil%ireadwf/=0.and.dtfil%ireadden==0.and.initialized==0) then 1356 rhor(:,:)=rhor(:,:)+nhat(:,:) 1357 if(dtset%usewvl==0) then 1358 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 1359 end if 1360 end if 1361 call timab(558,2,tsec) 1362 end if 1363 1364 ! The following steps have been gathered in the setvtr routine: 1365 ! - get Ewald energy and Ewald forces 1366 ! - compute local ionic pseudopotential vpsp 1367 ! - possibly compute 3D core electron density xccc3d 1368 ! - possibly compute 3D core kinetic energy density 1369 ! - possibly compute vxc and vhartr 1370 ! - set up vtrial 1371 1372 optene = 4 * optres 1373 if(dtset%iscf==-3) optene=4 1374 if (wvlbigdft) optene = 1 ! VH needed for the WF mixing 1375 1376 if (.not.allocated(nhatgr)) then 1377 ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3*nhatgrdim)) 1378 dummy_nhatgr = .True. 1379 end if 1380 1381 ! Compute trial potential 1382 ABI_NVTX_START_RANGE(NVTX_SCFCV_SETVTR) 1383 call setvtr(atindx1,dtset,energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,& 1384 & istep,kxc,mgfftf,moved_atm_inside,moved_rhor,mpi_enreg,& 1385 & nattyp,nfftf,ngfftf,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,psps%ntypat,& 1386 & n1xccc,n3xccc,optene,pawang,pawrad,pawrhoij,pawtab,ph1df,psps,rhog,rhor,& 1387 & rmet,rprimd,strsxc,ucvol,usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,& 1388 & xccc3d,xred,electronpositron=electronpositron,& 1389 & taur=taur,vxc_hybcomp=vxc_hybcomp,vxctau=vxctau,add_tfw=tfw_activated,xcctau3d=xcctau3d) 1390 ABI_NVTX_END_RANGE() 1391 1392 ! set the zero of the potentials here 1393 if(dtset%usepotzero==2) vpsp(:) = vpsp(:) + ecore / ( zion * ucvol ) 1394 1395 if(dtset%optdriver==RUNL_GWLS) call build_vxc(vxc,nfftf,dtset%nspden) 1396 1397 if ((nhatgrdim>0.and.nstep>0).or.dummy_nhatgr) then 1398 ABI_FREE(nhatgr) 1399 end if 1400 1401 ! Recursion Initialisation 1402 if(dtset%userec==1 .and. istep==1) then 1403 rec_set%quitrec = 0 1404 ! --At any step calculate the metric 1405 call Init_MetricRec(rec_set%inf,rec_set%nl%nlpsp,rmet,ucvol,rprimd,xred,dtset%ngfft(1:3),dtset%natom,rec_set%debug) 1406 call destroy_distribfft(rec_set%mpi%distribfft) 1407 call init_distribfft(rec_set%mpi%distribfft,'c',rec_set%mpi%nproc_fft,rec_set%ngfftrec(2),rec_set%ngfftrec(3)) 1408 call init_distribfft(rec_set%mpi%distribfft,'f',rec_set%mpi%nproc_fft,dtset%ngfft(2),dtset%ngfft(3)) 1409 if(initialized==0) call first_rec(dtset,psps,rec_set) 1410 end if 1411 1412 ! End the condition of atomic position change or istep==1 1413 end if 1414 1415 call timab(1448,2,tsec) 1416 call timab(1449,1,tsec) 1417 1418 ! ###################################################################### 1419 ! The following steps are done at every iteration 1420 ! ---------------------------------------------------------------------- 1421 ! PAW: Compute energies and potentials in the augmentation regions (spheres) 1422 ! Compute pseudopotential strengths (Dij quantities) 1423 if (psps%usepaw==1)then 1424 1425 ! Local exact exch.: impose occ. matrix if required 1426 if (dtset%useexexch/=0) then 1427 call setrhoijpbe0(dtset,initialized0,istep,istep_mix,& 1428 & spaceComm,my_natom,dtset%natom,dtset%ntypat,pawrhoij,pawtab,dtset%typat,& 1429 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1430 end if 1431 1432 ! Computation of on-site densities/potentials/energies 1433 nzlmopt=0;if (istep_mix==2.and.dtset%pawnzlm>0) nzlmopt=-1 1434 if (istep_mix>2) nzlmopt=dtset%pawnzlm 1435 call paw_an_reset_flags(paw_an) ! Force the recomputation of on-site potentials 1436 call paw_ij_reset_flags(paw_ij,self_consistent=.true.) ! Force the recomputation of Dij 1437 option=0;if (dtset%iscf>0.and.dtset%iscf<10.and.nstep>0) option=1 1438 ABI_NVTX_START_RANGE(NVTX_SCFCV_PAWDENPOT) 1439 call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,ipert,dtset%ixc,my_natom,dtset%natom,& 1440 & dtset%nspden,psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,pawrad,& 1441 & pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,& 1442 & dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp,& 1443 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 1444 & hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,& 1445 & electronpositron=electronpositron,vpotzero=vpotzero) 1446 ABI_NVTX_END_RANGE() 1447 1448 ! Correct the average potential with the calculated constant vpotzero 1449 ! Correct the total energies accordingly 1450 ! vpotzero(1) = -beta/ucvol 1451 ! vpotzero(2) = -1/ucvol sum_ij rho_ij gamma_ij 1452 write(msg,'(a,f14.6,2x,f14.6)') & 1453 & ' average electrostatic smooth potential [Ha] , [eV]',SUM(vpotzero(:)),SUM(vpotzero(:))*Ha_eV 1454 call wrtout(std_out, msg) 1455 vtrial(:,:)=vtrial(:,:)+SUM(vpotzero(:)) 1456 if(option/=1)then 1457 ! Fix the direct total energy (non-zero only for charged systems) 1458 energies%e_paw=energies%e_paw-SUM(vpotzero(:))*dtset%cellcharge(1) 1459 ! Fix the double counting total energy accordingly (for both charged AND 1460 ! neutral systems) 1461 energies%e_pawdc=energies%e_pawdc-SUM(vpotzero(:))*zion+vpotzero(2)*dtset%cellcharge(1) 1462 end if 1463 1464 ! PAW+U: impose density matrix if required 1465 ! not available if usepawu<0 (PAW+U without occupation matrix) 1466 if (dtset%usepawu>0.and.(ipositron/=1)) then 1467 impose_dmat=0 1468 if ((istep<=abs(dtset%usedmatpu)).and.(dtset%usedmatpu<0.or.initialized0==0)) impose_dmat=1 1469 if (impose_dmat==1.or.dtset%dmatudiag/=0) then 1470 dimdmat=0;if (impose_dmat==1) dimdmat=2*lpawumax+1 1471 call setnoccmmp(0,dimdmat,& 1472 & dmatpawu(1:dimdmat,1:dimdmat,1:dtset%nsppol*dtset%nspinor,1:dtset%natpawu*impose_dmat),& 1473 & dtset%dmatudiag,impose_dmat,indsym,my_natom,dtset%natom,dtset%natpawu,& 1474 & dtset%nspinor,dtset%nsppol,dtset%nsym,dtset%ntypat,paw_ij,pawang,dtset%pawprtvol,& 1475 & pawrhoij,pawtab,dtset%spinat,dtset%symafm,dtset%typat,0,dtset%usepawu,& 1476 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1477 ! Reinitalize mixing if PAW+U and occupation matrix now allowed to change 1478 ! For experimental purpose... 1479 if ((dtset%userib==1234).and.(istep==abs(dtset%usedmatpu)).and. & 1480 & (dtset%usedmatpu<0.or.initialized0==0)) reset_mixing=.true. 1481 end if 1482 end if 1483 1484 ! Write out unperturbed occupancies to dtpawuj-dataset LMac 1485 if (dtset%usepawu/=0.and.dtset%macro_uj>0.and.istep==1.and.ipositron/=1) then 1486 call pawuj_red(istep, 0, dtfil, dtset,dtpawuj,fatvshift,my_natom,dtset%natom,dtset%ntypat,& 1487 paw_ij,pawrad,pawtab,ndtpawuj,spaceComm, comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1488 end if 1489 ABI_NVTX_START_RANGE(NVTX_SCFCV_DIJ) 1490 1491 ! Dij computation 1492 call timab(561,1,tsec) 1493 1494 call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,dtset%natom,nfftf,nfftotf,& 1495 & dtset%nspden,psps%ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,& 1496 & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,dtset%spnorbscl,& 1497 & ucvol_local,dtset%cellcharge(1),vtrial,vxc,xred,& 1498 & natvshift=dtset%natvshift,atvshift=dtset%atvshift,fatvshift=fatvshift,& 1499 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 1500 & mpi_comm_grid=spaceComm_grid,& 1501 & hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,& 1502 & electronpositron_calctype=ipositron,& 1503 & electronpositron_pawrhoij=pawrhoij_ep,& 1504 & electronpositron_lmselect=lmselect_ep,& 1505 & nucdipmom=dtset%nucdipmom) 1506 1507 ! Symetrize Dij 1508 call symdij(gprimd,indsym,ipert,my_natom,dtset%natom,dtset%nsym,& 1509 & psps%ntypat,0,paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,& 1510 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1511 if (has_dijhat==1) then 1512 call symdij(gprimd,indsym,ipert,my_natom,dtset%natom,dtset%nsym,& 1513 & psps%ntypat,1,paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,& 1514 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1515 end if 1516 if(dtset%usewvl==1) then 1517 call paw2wvl_ij(3,paw_ij,wvl%descr) 1518 end if 1519 1520 call timab(561,2,tsec) 1521 ABI_NVTX_END_RANGE() 1522 end if 1523 1524 ! Now that the perturbation has been applied, we harvest occupancies for the perturbed case: LMac 1525 if (dtset%usepawu/=0.and.dtset%macro_uj>0.and.istep>1.and.ipositron/=1) then 1526 call pawuj_red(istep, 1, dtfil, dtset,dtpawuj,fatvshift,my_natom,dtset%natom,dtset%ntypat,& 1527 paw_ij,pawrad,pawtab,ndtpawuj,spaceComm,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1528 end if 1529 1530 call timab(1449,2,tsec) 1531 1532 ! No need to continue and call vtorho, when nstep==0 1533 if(nstep==0)exit 1534 1535 ! ###################################################################### 1536 ! The following steps are done only when nstep>0 1537 ! ---------------------------------------------------------------------- 1538 call timab(1450,1,tsec) 1539 1540 if(dtset%iscf>=0)then 1541 write(msg, '(a,a,i4)' )ch10,' ITER STEP NUMBER ',istep 1542 call wrtout(std_out,msg) 1543 end if 1544 1545 ! The next flag says whether the xred have to be changed in the current iteration 1546 moved_atm_inside=0 1547 ! /< Hack to remove iapp from scfcv_core 1548 ! for ionmov 4|5 ncycle=1 1549 ! Hence iapp = itime 1550 if ( dtset%jdtset /= scfcv_jdtset ) then 1551 ! new dtset -> reinitialize 1552 scfcv_jdtset = dtset%jdtset 1553 scfcv_itime = 0 1554 end if 1555 if ( istep==1 ) scfcv_itime = scfcv_itime + 1 1556 if(dtset%ionmov==4 .and. mod(scfcv_itime,2)/=1 .and. dtset%iscf>=0 ) moved_atm_inside=1 1557 if(dtset%ionmov==5 .and. scfcv_itime/=1 .and. istep==1 .and. dtset%iscf>=0) moved_atm_inside=1 1558 ! /< Hack to remove iapp from scfcv_core 1559 1560 ! Thomas-Fermi scheme might use a different toldfe criterion 1561 if (dtset%tfkinfunc>0.and.dtset%tfkinfunc/=2) then 1562 tollist(4)=dtset%toldfe;if (.not.tfw_activated) tollist(4)=dtset%tfw_toldfe 1563 end if 1564 1565 ! The next flag says whether the forces have to be computed in the current iteration 1566 computed_forces=0 1567 if ((dtset%optforces==1 .and. dtset%usewvl == 0).or.(moved_atm_inside==1)) computed_forces=1 1568 if (abs(tollist(3))>tiny(0._dp)) computed_forces=1 1569 if (dtset%iscf<0) computed_forces=0 1570 if ((istep==1).and.(dtset%optforces/=1)) then 1571 if (moved_atm_inside==1) then 1572 write(msg,'(5a)')& 1573 & 'Although the computation of forces during electronic iterations',ch10,& 1574 & 'was not required by user, it is done (required by the',ch10,& 1575 & 'choice of ionmov input parameter).' 1576 ABI_WARNING(msg) 1577 end if 1578 if (abs(tollist(3))+abs(tollist(7))>tiny(0._dp)) then 1579 write(msg,'(5a)')& 1580 & 'Although the computation of forces during electronic iterations',ch10,& 1581 & 'was not required by user, it is done (required by the',ch10,& 1582 & '"toldff" or "tolrff" tolerance criteria).' 1583 ABI_WARNING(msg) 1584 end if 1585 end if 1586 if ((istep==1).and.(dtset%optforces==1).and. dtset%usewvl == 1) then 1587 write(msg,'(5a)')& 1588 & 'Although the computation of forces during electronic iterations',ch10,& 1589 & 'was required by user, it has been disable since the tolerence',ch10,& 1590 & 'is not on forces (force computation is expensive in wavelets).' 1591 ABI_WARNING(msg) 1592 end if 1593 1594 call timab(1450,2,tsec) 1595 1596 ! ###################################################################### 1597 ! Compute the density rho from the trial potential 1598 ! ---------------------------------------------------------------------- 1599 call timab(1451,1,tsec) 1600 ! Compute the density from the trial potential 1601 if (dtset%tfkinfunc==0) then 1602 if(VERBOSE) call wrtout(std_out,'*. Compute the density from the trial potential (vtorho)') 1603 1604 ABI_NVTX_START_RANGE(NVTX_VTORHO) 1605 1606 call vtorho(afford,atindx,atindx1,cg,compch_fft,cprj,cpus,dbl_nnsclo,& 1607 & dielop,dielstrt,dmatpawu,dphase,dtefield,dtfil,dtset,& 1608 & eigen,electronpositron,energies,etotal,gbound_diel,& 1609 & gmet,gprimd,grnl,gsqcut,hdr,extfpmd,indsym,irrzon,irrzondiel,& 1610 & istep,istep_mix,itimes,kg,kg_diel,kxc,lmax_diel,mcg,mcprj,mgfftdiel,mpi_enreg,& 1611 & my_natom,dtset%natom,nattyp,nfftf,nfftdiel,ngfftdiel,nhat,nkxc,& 1612 & npwarr,npwdiel,res2,psps%ntypat,nvresid,occ,& 1613 & computed_forces,optres,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,& 1614 & pawrhoij,pawtab,phnons,phnonsdiel,ph1d,ph1ddiel,psps,fock,& 1615 & pwind,pwind_alloc,pwnsfac,resid,residm,rhog,rhor,rmet,rprimd,& 1616 & susmat,symrec,taug,taur,nvtauresid,ucvol_local,usecprj,wffnew,with_vectornd,& 1617 & vectornd,vtrial,vxctau,wvl,xred,ylm,ylmgr,ylmdiel, rmm_diis_status) 1618 1619 ABI_NVTX_END_RANGE() 1620 1621 else if (dtset%tfkinfunc==1.or.dtset%tfkinfunc==11.or.dtset%tfkinfunc==12) then 1622 ! CP: occopt 9 not available with Thomas Fermi functionals 1623 ABI_WARNING('THOMAS FERMI') 1624 call vtorhotf(dtset,energies%e_kinetic,energies%e_nlpsp_vfock,& 1625 & energies%entropy,energies%e_fermie,gprimd,grnl,irrzon,mpi_enreg,& 1626 & dtset%natom,nfftf,dtset%nspden,dtset%nsppol,dtset%nsym,phnons,& 1627 & rhog,rhor,rprimd,ucvol,vtrial) 1628 1629 residm=zero 1630 energies%e_eigenvalues=zero 1631 end if 1632 1633 ! Recursion method 1634 if(dtset%userec==1)then 1635 call vtorhorec(dtset,& 1636 & energies%e_kinetic,energies%e_nlpsp_vfock,energies%entropy,energies%e_eigenvalues,& 1637 & energies%e_fermie,grnl,initialized,irrzon,nfftf,phnons,& 1638 & rhog,rhor,vtrial,rec_set,istep-nstep,rprimd,gprimd) 1639 residm=zero 1640 end if 1641 1642 ! Update Fermi level in energies 1643 results_gs%fermie = energies%e_fermie 1644 1645 if(dtset%wfoptalg==2)then 1646 do ikpt=1,dtset%nkpt 1647 shiftvector(1+(ikpt-1)*(dtset%mband+2))=val_min 1648 shiftvector(2+(ikpt-1)*(dtset%mband+2):ikpt*(dtset%mband+2)-1)=& 1649 & eigen((ikpt-1)*dtset%mband+1:ikpt*dtset%mband) 1650 shiftvector(ikpt*(dtset%mband+2))=val_max 1651 end do 1652 end if 1653 1654 call timab(1451,2,tsec) 1655 1656 ! ###################################################################### 1657 ! Skip out of step loop if non-SCF (completed) 1658 ! ---------------------------------------------------------------------- 1659 1660 ! Indeed, nstep loops have been done inside vtorho 1661 if (dtset%iscf<0) exit 1662 1663 ! ###################################################################### 1664 ! In case of density mixing or wavelet handling, compute the total energy 1665 ! ---------------------------------------------------------------------- 1666 1667 call timab(1452,1,tsec) 1668 if (dtset%iscf>=10 .or. wvlbigdft) then 1669 optene = 1 ! use double counting scheme (default) 1670 if (wvlbigdft.and.dtset%iscf==0) optene = 0 ! use direct scheme 1671 if (dtset%iscf==22) optene = -1 1672 1673 ! Add the Fock contribution to E_xc and E_xcdc if required 1674 if (dtset%usefock==1) then 1675 energies%e_fockdc=two*energies%e_fock 1676 end if 1677 1678 ! if the mixing is the ODA mixing, compute energy and new density here 1679 if (dtset%iscf==22) then 1680 call odamix(deltae,dtset,& 1681 & elast,energies,etotal,gprimd,gsqcut,kxc,mpi_enreg,& 1682 & my_natom,nfftf,ngfftf,nhat,nkxc,psps%ntypat,nvresid,n3xccc,optres,& 1683 & paw_ij,paw_an,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,& 1684 & red_ptot,psps,rhog,rhor,rprimd,strsxc,ucvol,psps%usepaw,& 1685 & usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,xccc3d,xred,& 1686 & taur=taur,vxctau=vxctau,add_tfw=tfw_activated) 1687 end if 1688 ! If the density mixing is required, compute the total energy here 1689 ! TODO: add nvtauresid if needed (for forces?) 1690 call etotfor(atindx1,deltae,diffor,dtefield,dtset,& 1691 & elast,electronpositron,energies,& 1692 & etotal,favg,fcart,fock,forold,gred,gmet,grchempottn,grcondft,gresid,grewtn,grhf,grnl,grvdw,& 1693 & grxc,gsqcut,extfpmd,indsym,kxc,maxfor,mgfftf,mpi_enreg,my_natom,& 1694 & nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,psps%ntypat,nvresid,n1xccc,n3xccc,& 1695 & optene,computed_forces,optres,pawang,pawfgrtab,pawrad,pawrhoij,pawtab,& 1696 & ph1df,red_ptot,psps,rhog,rhor,rmet,rprimd,symrec,synlgr,ucvol,& 1697 & psps%usepaw,vhartr,vpsp,vxc,vxctau,wvl%descr,wvl%den,xccc3d,xred) 1698 !if (wvlbigdft) energies_copy(energies,energies_wvl) ! TO BE ACTIVATED LATER 1699 end if 1700 call timab(1452,2,tsec) 1701 1702 ! ###################################################################### 1703 ! In case of density mixing, check the exit criterion 1704 ! ---------------------------------------------------------------------- 1705 if (dtset%iscf>=10.or.(wvlbigdft.and.dtset%iscf>0)) then 1706 ! Check exit criteria 1707 call timab(1453,1,tsec) 1708 choice=2 1709 if(paw_dmft%use_dmft==1) then 1710 call prtene(dtset,energies,std_out,psps%usepaw) 1711 end if 1712 call scprqt(choice,cpus,deltae,diffor,dtset,& 1713 & eigen,etotal,favg,fcart,energies%e_fermie,energies%e_fermih,dtfil%fnameabo_app_eig,& 1714 & dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,istep_fock_outer,istep_mix,dtset%kptns,& 1715 & maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,& 1716 & occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,& 1717 & psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,& 1718 & electronpositron=electronpositron,fock=fock) 1719 call timab(1453,2,tsec) 1720 1721 ! Check if we need to exit the loop 1722 call timab(1454,1,tsec) 1723 if (dtset%tfkinfunc>10.and.(.not.tfw_activated).and.quit==1) then 1724 quit=0;tfw_activated=.true.;reset_mixing=.true. 1725 end if 1726 if(dtset%userec==1.and.rec_set%quitrec==2)quit=1 1727 if (istep==nstep) quit=1 1728 quit_sum=quit 1729 call xmpi_sum(quit_sum,spaceComm,ierr) 1730 if (quit_sum>0) quit=1 1731 call timab(1454,2,tsec) 1732 1733 ! If criteria in scprqt say to quit, then exit the loop over istep. 1734 if (quit==1) exit 1735 end if 1736 1737 1738 ! ###################################################################### 1739 ! Mix the total density (if required) 1740 ! ---------------------------------------------------------------------- 1741 call timab(1455,1,tsec) 1742 1743 if (dtset%iscf>=10 .and.dtset%iscf/=22.and. .not. wvlbigdft ) then 1744 1745 ! If LDA dielectric matrix is used for preconditionning, has to update here Kxc 1746 if (nkxc>0.and.modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79) & 1747 & .and.((istep==1.or.istep==dielstrt).or.(dtset%iprcel>=100))) then 1748 optxc=10 1749 call xcdata_init(xcdata,dtset=dtset) 1750 ! to be adjusted for the call to rhotoxc 1751 nk3xc=1 1752 if(dtset%icoulomb==0 .and. dtset%usewvl==0) then 1753 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 1754 call rhotoxc(edum,kxc,mpi_enreg,nfftf,& 1755 & ngfftf,nhat,psps%usepaw,nhatgr,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 1756 & optxc,rhor,rprimd,dummy2,0,vxc,vxcavg_dum,xccc3d,xcdata,& 1757 & add_tfw=tfw_activated,taur=taur,vhartr=vhartr,vxctau=vxctau,xcctau3d=xcctau3d) 1758 else if(.not. wvlbigdft) then 1759 ! WVL case: 1760 call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, & 1761 & dtset%icoulomb, dtset%ixc, & 1762 & mpi_enreg, nfftf, ngfftf,& 1763 & nhat,psps%usepaw,& 1764 & dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, & 1765 & usexcnhat,psps%usepaw,dtset%usewvl,vhartr, vxc, vxcavg,& 1766 & wvl%descr,wvl%den,& 1767 & wvl%e,xccc3d,dtset%xclevel,dtset%xc_denpos) 1768 end if 1769 end if 1770 1771 ABI_NVTX_START_RANGE(NVTX_SCFCV_NEWRHO) 1772 call newrho(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,dtn_pc,& 1773 & dtset,etotal,fcart,pawfgr%fintocoa,& 1774 & gmet,grhf,gsqcut,initialized,ispmix,istep_mix,kg_diel,kxc,& 1775 & mgfftf,mix,pawfgr%coatofin,moved_atm_inside,mpi_enreg,my_natom,nattyp,nfftf,& 1776 & nfftmix,nfftmix_per_nfft,ngfftf,ngfftmix,nkxc,npawmix,npwdiel,nvresid,psps%ntypat,& 1777 & n1xccc,pawrhoij,pawtab,ph1df,psps,rhog,rhor,& 1778 & rprimd,susmat,psps%usepaw,vtrial,wvl%descr,wvl%den,xred,& 1779 & mix_mgga=mix_mgga,taug=taug,taur=taur,tauresid=nvtauresid) 1780 ABI_NVTX_END_RANGE() 1781 1782 end if ! iscf>=10 1783 1784 call timab(1455,2,tsec) 1785 1786 ! ###################################################################### 1787 ! Additional computation in case of an electric field or electric displacement field 1788 ! ---------------------------------------------------------------------- 1789 1790 call timab(1456,1,tsec) 1791 1792 call update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,& 1793 & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,& 1794 & dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,dtset%natom,nattyp,ngfft,dtset%nkpt,npwarr,& 1795 & dtset%ntypat,pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,& 1796 & pwind,pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,& 1797 & 1,quit,istep,ucvol,unit_out,psps%usepaw,xred,ylm,ylmgr) 1798 1799 call timab(1456,2,tsec) 1800 1801 ! ###################################################################### 1802 ! Compute the new potential from the trial density 1803 ! ---------------------------------------------------------------------- 1804 1805 call timab(1457,1,tsec) 1806 if(VERBOSE) call wrtout(std_out,'*. Compute the new potential from the trial density') 1807 1808 ! Set XC computation flag 1809 optxc=1 1810 if (nkxc>0) then 1811 ! MJV 2017 May 25: you should not be able to get here with iscf < 0 1812 if (dtset%iscf<0) optxc=2 1813 if (modulo(dtset%iprcel,100)>=61.and.(dtset%iprcel<71.or.dtset%iprcel>79).and. & 1814 & dtset%iscf<10.and. & 1815 & (dtset%iprcel>=100.or.istep==1.or.istep==dielstrt)) optxc=2 1816 if (dtset%iscf>=10.and.dtset%densfor_pred/=0.and.abs(dtset%densfor_pred)/=5) optxc=2 1817 if (optxc==2.and.dtset%xclevel==2.and.nkxc==2*min(dtset%nspden,2)-1) optxc=12 1818 end if 1819 1820 if (dtset%iscf/=22) then 1821 ! PAW: eventually recompute compensation density (and gradients) 1822 nhatgrdim=0 1823 if ( allocated(nhatgr) ) then 1824 ABI_FREE(nhatgr) 1825 end if 1826 if (psps%usepaw==1) then 1827 ider=-1;if (dtset%iscf>=10.and.((dtset%xclevel==2.and.dtset%pawnhatxc>0).or.usexcnhat==0)) ider=0 1828 if (dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0) ider=ider+2 1829 if (ipositron==1) ider=-1 1830 if (ider>0) then 1831 nhatgrdim=1 1832 ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3)) 1833 else 1834 ABI_MALLOC(nhatgr,(0,0,0)) 1835 end if 1836 if (ider>=0) then 1837 ABI_NVTX_START_RANGE(NVTX_SCFCV_PAWKNHAT) 1838 call timab(558,1,tsec) 1839 izero=0 1840 1841 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,dtset%natom,nfftf,ngfftf,& 1842 & nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,nhatgr,nhat,& 1843 & pawrhoij,pawrhoij,pawtab,k0,rprimd,ucvol_local,dtset%usewvl,xred,& 1844 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 1845 & comm_fft=spaceComm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,& 1846 & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl) 1847 1848 call timab(558,2,tsec) 1849 ABI_NVTX_END_RANGE() 1850 end if 1851 else 1852 ABI_MALLOC(nhatgr,(0,0,0)) 1853 end if 1854 1855 ! Compute new potential from the trial density 1856 optene=2*optres;if(psps%usepaw==1) optene=2 1857 ABI_NVTX_START_RANGE(NVTX_SCFCV_RHOTOV) 1858 call rhotov(constrained_dft,dtset,energies,gprimd,grcondft,gsqcut,intgres,istep,& 1859 & kxc,mpi_enreg,nfftf,ngfftf,nhat,nhatgr,nhatgrdim,nkxc,nvresid,n3xccc,& 1860 & optene,optres,optxc,pawang,pawrad,pawrhoij,pawtab,& 1861 & rhog,rhor,rprimd,strscondft,strsxc,ucvol_local,psps%usepaw,usexcnhat,& 1862 & vhartr,vnew_mean,vpsp,vres_mean,res2,vtrial,vxcavg,vxc,wvl,xccc3d,xred,& 1863 & electronpositron=electronpositron,taur=taur,vxctau=vxctau,vtauresid=nvtauresid,& 1864 & vxc_hybcomp=vxc_hybcomp,add_tfw=tfw_activated,xcctau3d=xcctau3d) 1865 ABI_NVTX_END_RANGE() 1866 1867 end if 1868 1869 call timab(1457,2,tsec) 1870 call timab(1452,1,tsec) 1871 1872 ! This is inside the loop, its not equivalent to the line 1821 1873 if(moved_atm_inside==1) xred_old(:,:)=xred(:,:) 1874 1875 if (dtset%iscf<10) then 1876 1877 if(VERBOSE) call wrtout(std_out,'Check exit criteria in case of potential mixing') 1878 1879 ! If the potential mixing is required, compute the total energy here 1880 ! PAW: has to compute here spherical terms 1881 if (psps%usepaw==1) then 1882 nzlmopt=0;if (istep_mix==1.and.dtset%pawnzlm>0) nzlmopt=-1 1883 if (istep_mix>1) nzlmopt=dtset%pawnzlm 1884 call paw_an_reset_flags(paw_an) ! Force the recomputation of on-site potentials 1885 option=2 1886 ABI_NVTX_START_RANGE(NVTX_SCFCV_PAWDENPOT) 1887 call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,ipert,& 1888 & dtset%ixc,my_natom,dtset%natom,dtset%nspden,& 1889 & psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,& 1890 & paw_ij,pawang,dtset%pawprtvol,pawrad,pawrhoij,dtset%pawspnorb,& 1891 & pawtab,dtset%pawxcdev,dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp,& 1892 & hyb_mixing=hyb_mixing,hyb_mixing_sr=hyb_mixing_sr,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 1893 & electronpositron=electronpositron) 1894 ABI_NVTX_END_RANGE() 1895 end if 1896 1897 ! Add the Fock contribution to E_xc and E_xcdc if required 1898 if (dtset%usefock==1) energies%e_fockdc=two*energies%e_fock 1899 1900 ABI_NVTX_START_RANGE(NVTX_SCFCV_ETOTFOR) 1901 if (.not.wvlbigdft) then 1902 ! TODO: add nvtauresid if needed (for forces?) 1903 call etotfor(atindx1,deltae,diffor,dtefield,dtset,& 1904 & elast,electronpositron,energies,& 1905 & etotal,favg,fcart,fock,forold,gred,gmet,grchempottn,grcondft,gresid,grewtn,grhf,grnl,grvdw,& 1906 & grxc,gsqcut,extfpmd,indsym,kxc,maxfor,mgfftf,mpi_enreg,my_natom,& 1907 & nattyp,nfftf,ngfftf,ngrvdw,nhat,nkxc,dtset%ntypat,nvresid,n1xccc, & 1908 & n3xccc,0,computed_forces,optres,pawang,pawfgrtab,pawrad,pawrhoij,& 1909 & pawtab,ph1df,red_ptot,psps,rhog,rhor,rmet,rprimd,symrec,synlgr,ucvol,& 1910 & psps%usepaw,vhartr,vpsp,vxc,vxctau,wvl%descr,wvl%den,xccc3d,xred) 1911 end if 1912 ABI_NVTX_END_RANGE() 1913 1914 end if 1915 call timab(1452,2,tsec) 1916 1917 ! ###################################################################### 1918 ! Check exit criteria in case of potential mixing or direct minimization 1919 ! ---------------------------------------------------------------------- 1920 if ((dtset%iscf<10.and.(.not.wvlbigdft)) .or. dtset%iscf == 0) then 1921 ! Check exit criteria 1922 call timab(1453,1,tsec) 1923 choice=2 1924 ABI_NVTX_START_RANGE(NVTX_SCFCV_SCPRQT) 1925 call scprqt(choice,cpus,deltae,diffor,dtset,& 1926 & eigen,etotal,favg,fcart,energies%e_fermie,energies%e_fermih,dtfil%fnameabo_app_eig,& 1927 & dtfil%filnam_ds(1),initialized0,dtset%iscf,istep,istep_fock_outer,istep_mix,dtset%kptns,& 1928 & maxfor,moved_atm_inside,mpi_enreg,dtset%nband,dtset%nkpt,nstep,& 1929 & occ,optres,prtfor,prtxml,quit,res2,resid,residm,response,tollist,& 1930 & psps%usepaw,vxcavg,dtset%wtk,xred,conv_retcode,& 1931 & electronpositron=electronpositron,fock=fock) 1932 ABI_NVTX_END_RANGE() 1933 call timab(1453,2,tsec) 1934 1935 ! Check if we need to exit the loop 1936 call timab(1454,1,tsec) 1937 if (dtset%tfkinfunc>10.and.(.not.tfw_activated).and.quit==1) then 1938 quit=0;tfw_activated=.true.;reset_mixing=.true. 1939 end if 1940 if (istep==nstep.and.psps%usepaw==1) quit=1 1941 if(dtset%userec==1 .and. rec_set%quitrec==2) quit=1 1942 quit_sum=quit 1943 call xmpi_sum(quit_sum,spaceComm,ierr) 1944 if (quit_sum > 0) quit=1 1945 1946 ! If criteria in scprqt say to quit, then exit the loop over istep. 1947 if (quit==1) then 1948 do ispden=1,dtset%nspden 1949 vtrial(:,ispden)=vtrial(:,ispden)+nvresid(:,ispden)+vres_mean(ispden) 1950 end do 1951 !if (dtset%usekden==1) then 1952 ! do ispden=1,dtset%nspden 1953 ! vxctau(:,ispden,1)=vxctau(:,ispden,1)+nvtauresid(:,ispden) 1954 ! end do 1955 !end if 1956 call timab(1454,2,tsec) ! Due to the exit instruction, two timab calls are needed 1957 exit ! exit the loop over istep 1958 end if 1959 call timab(1454,2,tsec) ! Due to the exit instruction, two timab calls are needed 1960 end if 1961 1962 ! ###################################################################### 1963 ! Mix the potential (if required) - Check exit criteria 1964 ! ---------------------------------------------------------------------- 1965 1966 call timab(1458,1,tsec) 1967 if (dtset%iscf<10 .and. dtset%iscf>0 .and. .not. wvlbigdft) then 1968 1969 if(VERBOSE) call wrtout(std_out,'*. Mix the potential (if required) - Check exit criteria') 1970 1971 ! Precondition the residual and forces, then determine the new vtrial 1972 ! (Warning: the (H)xc potential may have been subtracted from vtrial) 1973 1974 call newvtr(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,& 1975 & dtn_pc,dtset,etotal,fcart,pawfgr%fintocoa,& 1976 & gmet,grhf,gsqcut,initialized,ispmix,& 1977 & istep_mix,kg_diel,kxc,mgfftf,mix,pawfgr%coatofin,& 1978 & moved_atm_inside,mpi_enreg,my_natom,nattyp,nfftf,nfftmix,& 1979 & ngfftf,ngfftmix,nkxc,npawmix,npwdiel,& 1980 & nstep,psps%ntypat,n1xccc,& 1981 & pawrhoij,ph1df,psps,rhor,rprimd,susmat,psps%usepaw,& 1982 & vhartr,vnew_mean,vpsp,nvresid,vres_mean,vtrial,vxc,xred,& 1983 & nfftf,pawtab,rhog,wvl,& 1984 & mix_mgga=mix_mgga,vtau=vxctau,vtauresid=nvtauresid) 1985 1986 end if ! iscf<10 1987 1988 ! ###################################################################### 1989 ! END MINIMIZATION ITERATIONS 1990 ! ###################################################################### 1991 1992 if(VERBOSE) call wrtout(std_out,'*. END MINIMIZATION ITERATIONS') 1993 1994 ! The initialisation of the gstate run should be done when this point is reached 1995 initialized=1 1996 1997 ! This is to save the density for restart. 1998 if (iwrite_fftdatar(mpi_enreg)) then 1999 2000 if(dtset%prtden<0.or.dtset%prtkden<0) then 2001 ! Update the content of the header (evolving variables) 2002 ! Don't use parallelism over atoms because only me=0 accesses here 2003 bantot=hdr%bantot 2004 if (dtset%positron==0) then 2005 call hdr%update(bantot,etotal,energies%e_fermie,energies%e_fermih,residm,& 2006 & rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1)) 2007 else 2008 call hdr%update(bantot,electronpositron%e0,energies%e_fermie,energies%e_fermih,residm,& 2009 & rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1)) 2010 end if 2011 end if 2012 2013 if (dtset%prtden<0) then 2014 if (mod(istep-1,abs(dtset%prtden))==0) then 2015 isave_den=isave_den+1 2016 rdwrpaw=0 2017 call int2char4(mod(isave_den,2),tag) 2018 ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!') 2019 fildata=trim(dtfil%fnametmp_app_den)//'_'//trim(tag) 2020 if (dtset%iomode == IO_MODE_ETSF) fildata = nctk_ncify(fildata) 2021 call fftdatar_write_from_hdr("density",fildata,dtset%iomode,hdr,ngfftf,cplex1,nfftf,& 2022 & dtset%nspden,rhor,mpi_enreg,eigen=eigen) 2023 end if 2024 end if 2025 2026 if (dtset%prtkden<0) then 2027 if (mod(istep-1,abs(dtset%prtkden))==0) then 2028 isave_kden=isave_kden+1 2029 rdwrpaw=0 2030 call int2char4(mod(isave_kden,2),tag) 2031 ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!') 2032 fildata=trim(dtfil%fnametmp_app_kden)//'_'//trim(tag) 2033 if (dtset%iomode == IO_MODE_ETSF) fildata = nctk_ncify(fildata) 2034 ! output the Laplacian of density 2035 call fftdatar_write_from_hdr("kinedr",fildata,dtset%iomode,hdr,ngfftf,cplex1,nfftf,& 2036 & dtset%nspden,taur,mpi_enreg,eigen=eigen) 2037 end if 2038 end if 2039 2040 end if 2041 2042 ABI_FREE(nhatgr) 2043 2044 istep_mix=istep_mix+1 2045 if (reset_mixing) then 2046 istep_mix=1;reset_mixing=.false. 2047 end if 2048 if (ipositron/=0) electronpositron%istep_scf=electronpositron%istep_scf+1 2049 2050 call timab(1458,2,tsec) 2051 end do ! istep 2052 ABI_NVTX_END_RANGE() 2053 2054 ABI_FREE(rmm_diis_status) 2055 ABI_SFREE(nhatgr) 2056 2057 ! Avoid pending requests if itime == ntime. 2058 call xmpi_wait(quitsum_request,ierr) 2059 if (timelimit_exit == 1) istep = istep - 1 2060 2061 call timab(1459,1,tsec) 2062 2063 if (dtset%iscf > 0) then 2064 call ab7_mixing_deallocate(mix) 2065 if (dtset%usekden/=0) call ab7_mixing_deallocate(mix_mgga) 2066 end if 2067 2068 if (dtset%usefock==1)then 2069 if(wfmixalg/=0) call scf_history_free(scf_history_wf) 2070 end if 2071 2072 if (quit==1.and.nstep==1) initialized=1 2073 2074 !###################################################################### 2075 !Case nstep==0: compute energy based on incoming wf 2076 !---------------------------------------------------------------------- 2077 2078 if(nstep==0) then 2079 optene=2*psps%usepaw+optres 2080 energies%entropy=results_gs%energies%entropy !MT20070219: entropy is not recomputed in routine energy 2081 if (.not.allocated(nhatgr) ) then 2082 ABI_MALLOC(nhatgr,(0,0,0)) 2083 end if 2084 2085 call energy(cg,compch_fft,constrained_dft,dtset,electronpositron,& 2086 & energies,eigen,etotal,gsqcut,extfpmd,indsym,irrzon,kg,mcg,mpi_enreg,my_natom,& 2087 & nfftf,ngfftf,nhat,nhatgr,nhatgrdim,npwarr,n3xccc,& 2088 & occ,optene,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,pawrhoij,pawtab,& 2089 & phnons,ph1d,psps,resid,rhog,rhor,rprimd,strsxc,symrec,taug,taur,usexcnhat,& 2090 & vhartr,vtrial,vpsp,vxc,wvl%wfs,wvl%descr,wvl%den,wvl%e,xccc3d,xred,ylm,& 2091 & add_tfw=tfw_activated,vxctau=vxctau) 2092 2093 if (nhatgrdim>0) then 2094 ABI_FREE(nhatgr) 2095 end if 2096 2097 end if ! nstep==0 2098 2099 !###################################################################### 2100 !Additional steps after SC iterations, including force, stress, polarization calculation 2101 !---------------------------------------------------------------------- 2102 2103 if (dtset%userec==1) then 2104 call prtene(dtset,energies,ab_out,psps%usepaw) 2105 call prtene(dtset,energies,std_out,psps%usepaw) 2106 end if 2107 2108 !if (wvlbigdft) call energies_copy(energies_wvl,energies) ! TO BE ACTIVATED LATER 2109 2110 !PAW: if cprj=<p_lmn|Cnk> are in memory, 2111 !need to reorder them (from atom-sorted to unsorted) 2112 if (psps%usepaw==1.and.usecprj==1) then 2113 iorder_cprj=1 2114 call pawcprj_reorder(cprj,atindx1) 2115 if (dtset%positron/=0) then 2116 if (electronpositron%dimcprj>0) then 2117 call pawcprj_reorder(electronpositron%cprj_ep,atindx1) 2118 end if 2119 end if 2120 if (dtset%usewvl==1) then 2121 call wvl_cprjreorder(wvl%descr,atindx1) 2122 end if 2123 end if 2124 2125 !PAW: if cprj=<p_lmn|Cnk> are not in memory,need to compute them in some cases 2126 recompute_cprj = psps%usepaw ==1 .and. usecprj==0 .and. & 2127 & (dtset%prtwant ==2 .or. & 2128 & dtset%prtwant ==3 .or. & 2129 & dtset%prtnabla > 0 .or. & 2130 & dtset%prtdos ==3 .or. & 2131 & dtset%berryopt /=0 .or. & 2132 & dtset%kssform ==3 .or. & 2133 & dtset%pawfatbnd> 0 .or. & 2134 & dtset%pawprtwf > 0 .or. & 2135 & dtset%plowan_compute > 0 .or. & 2136 & dtset%userid .EQ. 1 ) 2137 2138 if(ANY(ABS(dtset%nucdipmom)>tol8)) recompute_cprj=.TRUE. 2139 if(dtset%berryopt == -2 .AND. dtset%orbmag /= 0) recompute_cprj=.TRUE. 2140 2141 if (recompute_cprj) then 2142 usecprj=1 2143 mband_cprj=dtset%mband/mpi_enreg%nproc_band 2144 mcprj=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol 2145 ABI_MALLOC(cprj_local,(dtset%natom,mcprj)) 2146 ncpgr = 0 ; ctocprj_choice = 1 2147 if (finite_efield_flag) then 2148 if (forces_needed /= 0 .and. stress_needed == 0) then 2149 ncpgr = 3 ; ctocprj_choice = 2 2150 else if (forces_needed /= 0 .and. stress_needed /= 0) then 2151 ncpgr = 9 ; ctocprj_choice = 23 2152 else if (forces_needed == 0 .and. stress_needed /= 0) then 2153 ncpgr = 6 ; ctocprj_choice = 3 2154 end if 2155 end if 2156 if (dtset%berryopt == -2 .AND. dtset%orbmag /= 0) then 2157 ncpgr=3; ctocprj_choice=5 2158 end if 2159 call pawcprj_alloc(cprj_local,ncpgr,dimcprj) 2160 cprj=> cprj_local 2161 iatom=0 ; iorder_cprj=1 ! cprj are not ordered 2162 call ctocprj(atindx,cg,ctocprj_choice,cprj_local,gmet,gprimd,& 2163 & iatom,idir,iorder_cprj,dtset%istwfk,kg,dtset%kptns,& 2164 & mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,& 2165 & dtset%mpw,dtset%natom,nattyp,dtset%nband,dtset%natom,ngfft,& 2166 & dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,dtset%nsppol,& 2167 & dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,& 2168 & dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr) 2169 end if 2170 2171 ! MRM print final Hartree energy components 2172 enonlocalpsp=energies%e_nlpsp_vfock-2.0d0*energies%e_fock0 2173 esum=energies%e_kinetic+energies%e_ewald+energies%e_corepsp+energies%e_hartree+energies%e_xc& 2174 &+energies%e_localpsp+enonlocalpsp+energies%e_fock0& 2175 &+energies%e_hybcomp_E0+energies%e_hybcomp_v0+energies%e_hybcomp_v+energies%e_vdw_dftd& 2176 &+energies%e_elecfield+energies%e_magfield 2177 2178 if (me == 0) then 2179 write(std_out,'(a1)')' ' 2180 write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------' 2181 write(std_out,'(a,2(es16.6,a))')' Ekinetic = : ',energies%e_kinetic ,' Ha ,',energies%e_kinetic*Ha_eV ,' eV' 2182 write(std_out,'(a,2(es16.6,a))')' Evext_l = : ',energies%e_localpsp ,' Ha ,',energies%e_localpsp*Ha_eV ,' eV' 2183 write(std_out,'(a,2(es16.6,a))')' Evext_nl = : ',enonlocalpsp ,' Ha ,',enonlocalpsp*Ha_eV ,' eV' 2184 write(std_out,'(a,2(es16.6,a))')' Epsp_core = : ',energies%e_corepsp ,' Ha ,',energies%e_corepsp*Ha_eV ,' eV' 2185 write(std_out,'(a,2(es16.6,a))')' Ehartree = : ',energies%e_hartree ,' Ha ,',energies%e_hartree*Ha_eV ,' eV' 2186 if(dtset%usefock==1) then 2187 write(std_out,'(a,2(es16.6,a))')' Efock = : ',energies%e_fock0 ,' Ha ,',energies%e_fock0*Ha_eV ,' eV' 2188 endif 2189 write(std_out,'(a,2(es16.6,a))')' Exc_ks = : ',energies%e_xc ,' Ha ,',energies%e_xc*Ha_eV ,' eV' 2190 if(abs(energies%e_vdw_dftd)>1.0d-6) then 2191 write(std_out,'(a,2(es16.6,a))')' EvdW-D = : ',energies%e_vdw_dftd ,' Ha ,',energies%e_vdw_dftd*Ha_eV ,' eV' 2192 endif 2193 if(abs(energies%e_elecfield)>1.0d-6) then 2194 write(std_out,'(a,2(es16.6,a))')' Eefield = : ',energies%e_elecfield ,' Ha ,',energies%e_elecfield*Ha_eV ,' eV' 2195 endif 2196 if(abs(energies%e_magfield)>1.0d-6) then 2197 write(std_out,'(a,2(es16.6,a))')' Emfield = : ',energies%e_magfield ,' Ha ,',energies%e_magfield*Ha_eV ,' eV' 2198 endif 2199 write(std_out,'(a,2(es16.6,a))')' Enn = : ',energies%e_ewald ,' Ha ,',energies%e_ewald*Ha_eV ,' eV' 2200 write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------' 2201 write(std_out,'(a,2(es16.6,a))')' Etot = : ',esum ,' Ha ,',esum*Ha_eV ,' eV' 2202 write(std_out,'(a98)')'-------------------------------------------------------------------------------------------------' 2203 end if ! end MRM printing energy components 2204 2205 call timab(1459,2,tsec) 2206 call timab(1460,1,tsec) 2207 2208 2209 !SHOULD CLEAN THE ARGS OF THIS ROUTINE 2210 call afterscfloop(atindx,atindx1,cg,computed_forces,cprj,cpus,& 2211 & deltae,diffor,dtefield,dtfil,dtset,eigen,electronpositron,elfr,& 2212 & energies,etotal,favg,fcart,fock,forold,grchempottn,grcondft,& 2213 & gred,gresid,grewtn,grhf,grhor,grvdw,& 2214 & grxc,gsqcut,hdr,extfpmd,indsym,intgres,irrzon,istep,istep_fock_outer,istep_mix,& 2215 & kg,kxc,lrhor,maxfor,mcg,mcprj,mgfftf,& 2216 & moved_atm_inside,mpi_enreg,my_natom,n3xccc,nattyp,nfftf,ngfft,ngfftf,ngrvdw,nhat,& 2217 & nkxc,npwarr,nvresid,occ,optres,paw_an,paw_ij,pawang,pawfgr,& 2218 & pawfgrtab,pawrad,pawrhoij,pawtab,pel,pel_cg,ph1d,ph1df,phnons,pion,prtfor,& 2219 & prtxml,psps,pwind,pwind_alloc,pwnsfac,res2,resid,residm,results_gs,& 2220 & rhog,rhor,rprimd,stress_needed,strscondft,strsxc,strten,symrec,synlgr,taug,& 2221 & taur,tollist,usecprj,vhartr,vpsp,vtrial,vxc,vxctau,vxcavg,wvl,& 2222 & xccc3d,xcctau3d,xred,ylm,ylmgr,dtset%cellcharge(1)*SUM(vpotzero(:)),conv_retcode) 2223 2224 !Before leaving the present routine, save the current value of xred. 2225 xred_old(:,:)=xred(:,:) 2226 2227 call timab(1460,2,tsec) 2228 2229 !###################################################################### 2230 !All calculations in scfcv_core are finished. Printing section 2231 !---------------------------------------------------------------------- 2232 2233 call timab(1461,1,tsec) 2234 2235 call outscfcv(atindx1,cg,compch_fft,compch_sph,cprj,dimcprj,dmatpawu,dtfil,& 2236 & dtset,ecut,eigen,electronpositron,elfr,etotal,& 2237 & gmet,gprimd,grhor,hdr,intgres,kg,lrhor,dtset%mband,mcg,mcprj,dtset%mgfft,& 2238 & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,my_natom,dtset%natom,nattyp,& 2239 & nfftf,ngfftf,nhat,dtset%nkpt,npwarr,dtset%nspden,& 2240 & dtset%nsppol,dtset%nsym,psps%ntypat,n3xccc,occ,paw_dmft,pawang,pawfgr,pawfgrtab,& 2241 & pawrad,pawrhoij,pawtab,paw_an,paw_ij,dtset%prtvol,psps,results_gs,& 2242 & rhor,rprimd,taur,ucvol,usecprj,vhartr,vpsp,vtrial,vxc,wvl%den,xccc3d,xred) 2243 2244 call timab(1461,2,tsec) 2245 call timab(1462,1,tsec) 2246 2247 !Transfer eigenvalues and occupation computed by BigDFT in afterscfloop to eigen. 2248 #if defined HAVE_BIGDFT 2249 if (dtset%usewvl == 1) then 2250 if (dtset%nsppol == 1) then 2251 eigen = wvl%wfs%ks%orbs%eval 2252 occ = wvl%wfs%ks%orbs%occup 2253 else 2254 eigen(1:wvl%wfs%ks%orbs%norbu) = wvl%wfs%ks%orbs%eval(1:wvl%wfs%ks%orbs%norbu) 2255 eigen(dtset%mband + 1:dtset%mband + wvl%wfs%ks%orbs%norbd) = & 2256 & wvl%wfs%ks%orbs%eval(wvl%wfs%ks%orbs%norbu + 1:wvl%wfs%ks%orbs%norb) 2257 occ(1:wvl%wfs%ks%orbs%norbu) = wvl%wfs%ks%orbs%occup(1:wvl%wfs%ks%orbs%norbu) 2258 occ(dtset%mband + 1:dtset%mband + wvl%wfs%ks%orbs%norbd) = & 2259 & wvl%wfs%ks%orbs%occup(wvl%wfs%ks%orbs%norbu + 1:wvl%wfs%ks%orbs%norb) 2260 end if 2261 end if 2262 #endif 2263 !need to reorder cprj (from unsorted to atom-sorted) 2264 if (psps%usepaw==1.and.usecprj==1) then 2265 iorder_cprj=0 2266 call pawcprj_reorder(cprj,atindx) 2267 if (dtset%positron/=0) then 2268 if (electronpositron%dimcprj>0) then 2269 call pawcprj_reorder(electronpositron%cprj_ep,atindx) 2270 end if 2271 end if 2272 end if 2273 !###################################################################### 2274 !Deallocate memory and save results 2275 !---------------------------------------------------------------------- 2276 2277 call prc_mem_free() 2278 2279 ABI_FREE(fcart) 2280 ABI_FREE(gred) 2281 ABI_FREE(forold) 2282 ABI_FREE(grchempottn) 2283 ABI_FREE(grcondft) 2284 ABI_FREE(gresid) 2285 ABI_FREE(grewtn) 2286 ABI_FREE(grnl) 2287 ABI_FREE(grvdw) 2288 ABI_FREE(grxc) 2289 ABI_FREE(intgres) 2290 ABI_FREE(synlgr) 2291 ABI_FREE(ph1d) 2292 ABI_FREE(ph1df) 2293 ABI_FREE(vhartr) 2294 ABI_FREE(vtrial) 2295 ABI_FREE(vpsp) 2296 ABI_FREE(vxc) 2297 ABI_FREE(vxc_hybcomp) 2298 ABI_FREE(vxctau) 2299 ABI_FREE(xccc3d) 2300 ABI_FREE(kxc) 2301 ABI_FREE(shiftvector) 2302 ABI_FREE(dtn_pc) 2303 ABI_FREE(grhf) 2304 ABI_FREE(nvresid) 2305 ABI_FREE(nvtauresid) 2306 2307 if(allocated(vectornd)) then 2308 ABI_FREE(vectornd) 2309 end if 2310 2311 if((nstep>0.and.dtset%iscf>0).or.dtset%iscf==-1) then 2312 ABI_FREE(dielinv) 2313 end if 2314 ABI_FREE(gbound_diel) 2315 ABI_FREE(irrzondiel) 2316 ABI_FREE(kg_diel) 2317 ABI_FREE(phnonsdiel) 2318 ABI_FREE(susmat) 2319 ABI_FREE(ph1ddiel) 2320 ABI_FREE(ylmdiel) 2321 2322 if (psps%usepaw==1) then 2323 if (dtset%iscf>0) then 2324 do iatom=1,my_natom 2325 pawrhoij(iatom)%lmnmix_sz=0 2326 pawrhoij(iatom)%use_rhoijres=0 2327 ABI_FREE(pawrhoij(iatom)%kpawmix) 2328 ABI_FREE(pawrhoij(iatom)%rhoijres) 2329 end do 2330 end if 2331 ! if (recompute_cprj.or.usecprj==1) then 2332 if (recompute_cprj) then 2333 usecprj=0;mcprj=0 2334 call pawcprj_free(cprj) 2335 ABI_FREE(cprj_local) 2336 end if 2337 call paw_an_free(paw_an) 2338 call paw_ij_free(paw_ij) 2339 call pawfgrtab_free(pawfgrtab) 2340 if(dtset%usewvl==1) then 2341 #if defined HAVE_BIGDFT 2342 call cprj_clean(wvl%descr%paw%cprj) 2343 ABI_FREE(wvl%descr%paw%cprj) 2344 #endif 2345 call paw2wvl_ij(2,paw_ij,wvl%descr) 2346 end if 2347 end if 2348 ABI_FREE(pawfgrtab) 2349 ABI_FREE(paw_an) 2350 ABI_FREE(paw_ij) 2351 ABI_FREE(nhat) 2352 ABI_FREE(xcctau3d) 2353 ABI_FREE(dimcprj_srt) 2354 ABI_FREE(dimcprj) 2355 2356 2357 ! Deallocate exact exchange data at the end of the calculation 2358 if (dtset%usefock==1) then 2359 if (fock%fock_common%use_ACE/=0) call fock_ACE_destroy(fock%fockACE) 2360 call fock_common_destroy(fock%fock_common) 2361 call fock_BZ_destroy(fock%fock_BZ) 2362 call fock_destroy(fock) 2363 nullify(fock) 2364 end if 2365 2366 if (prtxml == 1) then 2367 ! We output the final result given in results_gs 2368 write(ab_xml_out, "(A)") ' <finalConditions>' 2369 call out_resultsgs_XML(dtset, 4, results_gs, psps%usepaw) 2370 write(ab_xml_out, "(A)") ' </finalConditions>' 2371 write(ab_xml_out, "(A)") ' </scfcvLoop>' 2372 end if 2373 2374 !Free the datastructure constrained_dft 2375 call constrained_dft_free(constrained_dft) 2376 2377 call timab(1462,2,tsec) 2378 call timab(1440,2,tsec) 2379 2380 DBG_EXIT("COLL") 2381 2382 end subroutine scfcv_core
ABINIT/wf_mixing [ Functions ]
NAME
wf_mixing
FUNCTION
Mixing of wavefunctions in the outer loop of a double loop SCF approach. Different algorithms are implemented, depending on the value of wfmixalg.
INPUTS
atindx1(dtset%natom)=index table for atoms, inverse of atindx dtset <type(dataset_type)>=all input variables in this dataset istep=number of call the routine (usually the outer loop in the SCF double loop) mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of cprj array mpi_enreg=information about MPI parallelization nattyp(dtset%ntypat)=number of atoms of each type in cell. npwarr(nkpt)=number of planewaves in basis at this k point pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
SIDE EFFECTS
cg(2,mcg)= plane wave wavefunction coefficient
Value from previous SCF cycle is input and stored in some form
Extrapolated value is output
cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
Value from previous SCF cycle is input and stored in some form
Extrapolated value is output
scf_history_wf <type(scf_history_type)>=arrays obtained from previous SCF cycles
SOURCE
2857 subroutine wf_mixing(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,& 2858 & nattyp,npwarr,pawtab,scf_history_wf) 2859 2860 use m_cgcprj, only : dotprod_set_cgcprj, dotprodm_sumdiag_cgcprj, lincom_cgcprj, cgcprj_cholesky 2861 2862 !Arguments ------------------------------------ 2863 !scalars 2864 integer,intent(in) :: istep,mcg,mcprj 2865 type(MPI_type),intent(in) :: mpi_enreg 2866 type(dataset_type),intent(in) :: dtset 2867 type(scf_history_type),intent(inout) :: scf_history_wf 2868 !arrays 2869 integer,intent(in) :: atindx1(dtset%natom),nattyp(dtset%ntypat) 2870 integer,intent(in) :: npwarr(dtset%nkpt) 2871 real(dp), intent(inout) :: cg(2,mcg) 2872 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj) 2873 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 2874 2875 !Local variables------------------------------- 2876 !scalars 2877 integer :: hermitian 2878 integer :: ibdmix,ibdsp,ibg,ibg_hist,icg,icg_hist 2879 integer :: ierr,ikpt,indh,ind_biorthog,ind_biorthog_eff,ind_newwf,ind_residual,inplace 2880 integer :: iorder,iset2,isppol,istep_cycle,istep_new,istwf_k,kk,me_distrb,my_nspinor 2881 integer :: nband_k,nbdmix,npw_k,nset1,nset2,ntypat 2882 integer :: shift_set1,shift_set2,spaceComm_band,spare_mem,usepaw,wfmixalg 2883 real(dp) :: alpha,beta 2884 complex(dpc) :: sum_coeffs 2885 !arrays 2886 integer,allocatable :: ipiv(:),dimcprj(:) 2887 real(dp) :: tsec(2) 2888 real(dp),allocatable :: al(:,:),mmn(:,:,:) 2889 real(dp),allocatable :: dotprod_res(:,:,:),dotprod_res_k(:,:,:),res_mn(:,:,:),smn(:,:,:) 2890 complex(dpc),allocatable :: coeffs(:) 2891 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kh(:,:) 2892 2893 ! ************************************************************************* 2894 2895 !DEBUG 2896 !write(std_out,*) 2897 !write(std_out,*)' wf_mixing : enter, istep= ',istep 2898 !call flush(std_out) 2899 !write(std_out,*)' istep,scf_history_wf%alpha=',istep,scf_history_wf%alpha 2900 !write(std_out,*)' cg(1,1)=',cg(1,1) 2901 !write(std_out,*)' scf_history_wf%cg(1,1,1:5)=',scf_history_wf%cg(1,1,1:5) 2902 !ABI_MALLOC(cg_ref,(2,mcg)) 2903 !cg_ref(:,:)=cg(:,:) 2904 !ABI_MALLOC(cprj_ref,(dtset%natom,mcprj)) 2905 !cprj_ref(:,:)=cprj(:,:) 2906 ! write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',& 2907 !& scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2) 2908 ! call flush(std_out) 2909 !ENDDEBUG 2910 2911 if (istep==0) return 2912 2913 ntypat=dtset%ntypat 2914 usepaw=dtset%usepaw 2915 wfmixalg=scf_history_wf%wfmixalg 2916 nbdmix=dtset%nbandhf 2917 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 2918 me_distrb=mpi_enreg%me_kpt 2919 spaceComm_band=xmpi_comm_self 2920 2921 spare_mem=0 2922 if(scf_history_wf%history_size==wfmixalg-1)spare_mem=1 2923 2924 !scf_history_wf%alpha contains dtset%wfmix 2925 alpha=scf_history_wf%alpha 2926 beta=one-scf_history_wf%alpha 2927 icg=0 2928 icg_hist=0 2929 ibg=0 2930 ibg_hist=0 2931 2932 !Useful array 2933 ABI_MALLOC(dimcprj,(dtset%natom)) 2934 if (usepaw==1) then 2935 iorder=0 ! There is no change of ordering in the mixing of wavefunctions 2936 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O') 2937 end if 2938 2939 if(istep==1)then 2940 do indh=1,scf_history_wf%history_size 2941 call pawcprj_alloc(scf_history_wf%cprj(:,:,indh),0,dimcprj) 2942 end do 2943 end if 2944 2945 ABI_MALLOC(cprj_k,(dtset%natom,my_nspinor*nbdmix)) 2946 ABI_MALLOC(cprj_kh,(dtset%natom,my_nspinor*nbdmix)) 2947 if(usepaw==1) then 2948 call pawcprj_alloc(cprj_k,0,dimcprj) 2949 call pawcprj_alloc(cprj_kh,0,dimcprj) 2950 end if 2951 ABI_MALLOC(smn,(2,nbdmix,nbdmix)) 2952 ABI_MALLOC(mmn,(2,nbdmix,nbdmix)) 2953 2954 if(wfmixalg>2)then 2955 nset1=1 2956 nset2=min(istep-1,wfmixalg-1) 2957 ABI_MALLOC(dotprod_res_k,(2,1,nset2)) 2958 ABI_MALLOC(dotprod_res,(2,1,nset2)) 2959 ABI_MALLOC(res_mn,(2,wfmixalg-1,wfmixalg-1)) 2960 dotprod_res=zero 2961 if(istep==1)then 2962 scf_history_wf%dotprod_sumdiag_cgcprj_ij=zero 2963 end if 2964 end if 2965 2966 !Explanation for the index for the wavefunction stored in scf_history_wf 2967 !The reference is the cg+cprj output after the wf optimization at istep 1. 2968 !It comes as input to the present routine as cgcprj input at step 2, and is usually found at indh=1. 2969 2970 !In the simple mixing case (wfmixalg==2), the reference is never stored, because it is used "on-the-fly" to biothogonalize the 2971 !previous input (that was stored in indh=1), then generate the next input, which is stored again in indh=1 2972 2973 !When the storage is not spared: 2974 !- the values of indh from 2 to wfmixalg store the (computed here) biorthogonalized input cgcprj, then the residual 2975 !- the values of indh from wfmixalg+1 to 2*wfmixalg-1 store the biorthogonalized output cgcprj (coming as argument) 2976 2977 !First step 2978 if (istep==1 .or. (wfmixalg==2 .and. abs(scf_history_wf%alpha-one)<tol8) ) then 2979 2980 indh=2 ! This input wavefunction is NOT the reference 2981 if(wfmixalg==2)indh=1 ! But this does not matter in the simple mixing case that has history_size=1 2982 2983 ! Simply store the wavefunctions and cprj. However, nband_k might be different from nbandhf... 2984 ! LOOP OVER SPINS 2985 do isppol=1,dtset%nsppol 2986 2987 ! BIG FAT k POINT LOOP 2988 do ikpt=1,dtset%nkpt 2989 2990 ! Select k point to be treated by this proc 2991 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 2992 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 2993 2994 npw_k=npwarr(ikpt) 2995 2996 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 2997 if(usepaw==1) then 2998 ! scf_history_wf%cprj(:,ibg_hist+1:ibg_hist+my_nspinor*nbdmix,1)=cprj(:,ibg+1:ibg+my_nspinor*nbdmix) 2999 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,dtset%mband,& 3000 & dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,& 3001 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 3002 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,& 3003 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 3004 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 3005 end if 3006 3007 ! Update the counters 3008 ibg=ibg+my_nspinor*nband_k 3009 ibg_hist=ibg_hist+my_nspinor*nbdmix 3010 icg=icg+my_nspinor*nband_k*npw_k 3011 icg_hist=icg_hist+my_nspinor*nbdmix*npw_k 3012 3013 end do 3014 end do 3015 3016 else 3017 ! From istep==2 3018 3019 ! First part of the computation : biorthogonalization, and computation of the residual (possibly, prediction of the next input in the case of simple mixing) 3020 ! Index for the wavefunctions stored in scf_history_wf whose scalar products with the argument cgcprj will have to be computed. 3021 indh=1 ! This input wavefunction is the reference 3022 if(wfmixalg/=2 .and. istep==2)indh=2 ! except for istep=2 in the rmm-diis 3023 3024 if(wfmixalg>2)then 3025 ! istep inside the cycle defined by wfmixalg, and next index. Then, indices of the wavefunction sets. 3026 istep_cycle=mod((istep-2),wfmixalg-1) 3027 istep_new=mod((istep-1),wfmixalg-1) 3028 ind_biorthog=1+wfmixalg+istep_cycle 3029 ind_residual=2+istep_cycle 3030 ind_newwf=2+istep_new 3031 shift_set1=ind_residual-1 3032 shift_set2=1 3033 end if 3034 3035 ! LOOP OVER SPINS 3036 do isppol=1,dtset%nsppol 3037 3038 ! BIG FAT k POINT LOOP 3039 do ikpt=1,dtset%nkpt 3040 3041 ! Select k point to be treated by this proc 3042 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 3043 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 3044 3045 istwf_k=dtset%istwfk(ikpt) 3046 npw_k=npwarr(ikpt) 3047 3048 ! Biorthogonalization 3049 3050 if(usepaw==1) then 3051 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,dtset%mband,& 3052 & dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,& 3053 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 3054 call pawcprj_get(atindx1,cprj_kh,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,& 3055 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,my_nspinor,dtset%nsppol,0,& 3056 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 3057 end if !end usepaw=1 3058 3059 hermitian=0 3060 if(wfmixalg==2 .or. istep==2)then 3061 call dotprod_set_cgcprj(atindx1,cg,scf_history_wf%cg(:,:,indh),cprj_k,cprj_kh,dimcprj,hermitian,& 3062 & 0,0,icg,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,& 3063 & mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw) 3064 else 3065 call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),cg,cprj_kh,cprj_k,dimcprj,hermitian,& 3066 & 0,0,icg_hist,icg,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,& 3067 & mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw) 3068 end if 3069 3070 ! Invert S matrix, that is NOT hermitian. 3071 ! Calculate M=S^-1 3072 mmn=zero 3073 do kk=1,nbdmix 3074 mmn(1,kk,kk)=one 3075 end do 3076 3077 ABI_MALLOC(ipiv,(nbdmix)) 3078 ! The smn is destroyed by the following inverse call 3079 call zgesv(nbdmix,nbdmix,smn,nbdmix,ipiv,mmn,nbdmix,ierr) 3080 ABI_CHECK(ierr == 0, sjoin('zgesv general inversion routine returned ierr:', itoa(ierr))) 3081 ABI_FREE(ipiv) 3082 3083 ! The M matrix is used to compute the biorthogonalized set of wavefunctions, and to store it at the proper place 3084 if(wfmixalg==2 .or. istep==2)then 3085 inplace=1 3086 call lincom_cgcprj(mmn,scf_history_wf%cg(:,:,indh),cprj_kh,dimcprj,& 3087 & icg_hist,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw) 3088 else 3089 inplace=0 3090 call lincom_cgcprj(mmn,cg,cprj_k,dimcprj,& 3091 & icg,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw,& 3092 & cgout=scf_history_wf%cg(:,:,ind_biorthog),cprjout=scf_history_wf%cprj(:,:,ind_biorthog),icgout=icg_hist) 3093 end if 3094 3095 ! The biorthogonalised set of wavefunctions is now stored at the proper place 3096 3097 ! Finalize this first part of the computation, depending on the algorithm and the step. 3098 3099 if(wfmixalg==2)then 3100 3101 ! Wavefunction extrapolation, simple mixing case 3102 ! alpha contains dtset%wfmix, beta contains one-alpha 3103 cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=& 3104 & alpha*cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)& 3105 & +beta*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh) 3106 if(usepaw==1) then 3107 do ibdmix=1,nbdmix 3108 call pawcprj_axpby(beta,alpha,cprj_kh(:,ibdmix:ibdmix),cprj_k(:,ibdmix:ibdmix)) 3109 end do ! end loop on ibdmix 3110 end if 3111 3112 ! Back to usual orthonormalization 3113 call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,& 3114 & mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw) 3115 3116 ! Store the newly extrapolated wavefunctions, orthonormalized, in scf_history_wf 3117 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 3118 if(usepaw==1) then 3119 do ibdmix=1,nbdmix 3120 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,& 3121 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 3122 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 3123 end do ! end loop on ibdmix 3124 end if 3125 3126 else ! wfmixalg/=2 3127 ! RMM-DIIS 3128 3129 if (istep==2)then 3130 ! Store the argument wf as the reference for all future steps, in scf_history_wf with index 1. 3131 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 3132 if(usepaw==1) then 3133 do ibdmix=1,nbdmix 3134 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,1),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,& 3135 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 3136 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 3137 end do ! end loop on ibdmix 3138 end if 3139 end if 3140 3141 ind_biorthog_eff=ind_biorthog 3142 if(istep==2)ind_biorthog_eff=1 ! The argument wf has not been stored in ind_biorthog 3143 ! Compute the residual of the wavefunctions for this istep, 3144 ! that replaces the previously stored set of (biorthogonalized) input wavefunctions 3145 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)=& 3146 & scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)& 3147 & -scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual) 3148 if(usepaw==1) then 3149 do ibdmix=1,nbdmix 3150 call pawcprj_axpby(one,-one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff),& 3151 & scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual)) 3152 end do ! end loop on ibdmix 3153 end if 3154 3155 ! Compute the new scalar products to fill the res_mn matrix 3156 call dotprodm_sumdiag_cgcprj(atindx1,scf_history_wf%cg,scf_history_wf%cprj,dimcprj,& 3157 & ibg_hist,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcprj,dtset%mkmem,& 3158 & mpi_enreg,scf_history_wf%history_size,dtset%natom,nattyp,nbdmix,npw_k,nset1,nset2,my_nspinor,dtset%nsppol,ntypat,& 3159 & shift_set1,shift_set2,pawtab,dotprod_res_k,usepaw) 3160 3161 dotprod_res=dotprod_res+dotprod_res_k 3162 3163 ! scf_history_wf for index ind_biorthog will contain the extrapolated wavefunctions (and no more the output of the SCF loop). 3164 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog)=& 3165 & scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)+& 3166 & (alpha-one)*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual) 3167 if(usepaw==1) then 3168 do ibdmix=1,nbdmix 3169 if(ind_biorthog/=ind_biorthog_eff)then 3170 scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog)=scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff) 3171 end if 3172 call pawcprj_axpby((alpha-one),one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual),& 3173 & scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog)) 3174 end do ! end loop on ibdmix 3175 end if 3176 3177 end if 3178 3179 ibg=ibg+my_nspinor*nband_k 3180 ibg_hist=ibg_hist+my_nspinor*nbdmix 3181 icg=icg+my_nspinor*nband_k*npw_k 3182 icg_hist=icg_hist+my_nspinor*nbdmix*npw_k 3183 3184 end do ! End big k point loop 3185 end do ! End loop over spins 3186 3187 end if ! istep>=2 3188 3189 if(wfmixalg>2 .and. istep>1)then 3190 3191 !DEBUG 3192 ! write(std_out,*)' ' 3193 ! write(std_out,*)' Entering the residual minimisation part ' 3194 ! write(std_out,*)' ' 3195 ! call flush(std_out) 3196 !ENDDEBUG 3197 3198 call timab(48,1,tsec) 3199 call xmpi_sum(dotprod_res,mpi_enreg%comm_kpt,ierr) 3200 call timab(48,2,tsec) 3201 3202 scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set1,1+shift_set2:nset2+shift_set2)=dotprod_res(:,1,1:nset2) 3203 scf_history_wf%dotprod_sumdiag_cgcprj_ij(1,1+shift_set2:nset2+shift_set2,1+shift_set1)=dotprod_res(1,1,1:nset2) 3204 scf_history_wf%dotprod_sumdiag_cgcprj_ij(2,1+shift_set2:nset2+shift_set2,1+shift_set1)=-dotprod_res(2,1,1:nset2) 3205 3206 end if ! wfmixalg>2 and istep>1 3207 3208 if(wfmixalg>2 .and. istep>2)then 3209 3210 ! Extract the relevant matrix R_mn 3211 res_mn(:,1:nset2,1:nset2)=& 3212 & scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set2:nset2+shift_set2,1+shift_set2:nset2+shift_set2) 3213 3214 !DEBUG 3215 ! write(std_out,*)' The matrix res_mn(:,1:nset2,1:nset2) is :' 3216 ! write(std_out,*)res_mn(:,1:nset2,1:nset2) 3217 ! call flush(std_out) 3218 !ENDDEBUG 3219 3220 ! Solve R_mn \alpha_n = 1_m 3221 ABI_MALLOC(ipiv,(nset2)) 3222 ABI_MALLOC(coeffs,(nset2)) 3223 coeffs(:)=cone 3224 ! The res_mn is destroyed by the following inverse call 3225 call zgesv(nset2,1,res_mn,wfmixalg-1,ipiv,coeffs,nset2,ierr) 3226 ABI_CHECK(ierr == 0, sjoin('zgesv general inversion routine returned ierr:', itoa(ierr))) 3227 ABI_FREE(ipiv) 3228 ! The coefficients must sum to one 3229 sum_coeffs=sum(coeffs) 3230 coeffs=coeffs/sum_coeffs 3231 3232 !DEBUG 3233 ! write(std_out,*)' The coefficients that minimize the residual have been found' 3234 ! write(std_out,*)' coeffs =',coeffs 3235 ! call flush(std_out) 3236 !ENDDEBUG 3237 end if ! wfmixalg>2 and istep>2 3238 3239 if(wfmixalg>2 .and. istep>1)then 3240 3241 ! Find the new "input" wavefunction, bi-orthogonalized, and store it replacing the adequate "old" input wavefunction. 3242 3243 icg=0 3244 icg_hist=0 3245 ibg=0 3246 ibg_hist=0 3247 ABI_MALLOC(al,(2,nset2)) 3248 if(istep>2)then 3249 do iset2=1,nset2 3250 al(1,iset2)=real(coeffs(iset2)) ; al(2,iset2)=aimag(coeffs(iset2)) 3251 end do 3252 else 3253 al(1,1)=one ; al(2,1)=zero 3254 end if 3255 3256 !DEBUG 3257 ! write(std_out,*)' Overload the coefficients, in order to simulate a simple mixing with wfmix ' 3258 ! write(std_out,*)' Set al(1,ind_biorthog-3)=one, for ind_biorthog=',ind_biorthog 3259 ! write(std_out,*)' This will feed scf_history for set ind_biorthog-3+wfmixalg=',ind_biorthog-3+wfmixalg 3260 ! al(:,:)=zero 3261 ! al(1,ind_biorthog-3)=one 3262 ! call flush(std_out) 3263 !ENDDEBUG 3264 3265 ! LOOP OVER SPINS 3266 do isppol=1,dtset%nsppol 3267 3268 ! BIG FAT k POINT LOOP 3269 do ikpt=1,dtset%nkpt 3270 3271 ! Select k point to be treated by this proc 3272 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 3273 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 3274 3275 istwf_k=dtset%istwfk(ikpt) 3276 npw_k=npwarr(ikpt) 3277 3278 if(istep>2)then 3279 ! Make the appropriate linear combination (from the extrapolated wfs) 3280 cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=zero 3281 do iset2=1,nset2 3282 cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)& 3283 & +al(1,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)& 3284 & -al(2,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg) 3285 cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)& 3286 & +al(1,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)& 3287 & +al(2,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg) 3288 end do 3289 else ! One needs a simple copy from the extrapolated wavefunctions 3290 cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1+wfmixalg) 3291 end if 3292 ! Note the storage in cprj_k. By the way, a simple copy might also be used in case istep=2. 3293 if(usepaw==1) then 3294 do ibdsp=1,my_nspinor*nbdmix 3295 call pawcprj_lincom(al,scf_history_wf%cprj(:,ibdsp,1+wfmixalg:nset2+wfmixalg),cprj_k(:,ibdsp:ibdsp),nset2) 3296 end do 3297 end if 3298 3299 ! Store the newly extrapolated wavefunctions for this k point, still bi-orthonormalized, in scf_history_wf 3300 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_newwf)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 3301 if(usepaw==1) then 3302 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,ind_newwf),dtset%natom,1,ibg_hist,ikpt,iorder,isppol,& 3303 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 3304 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 3305 end if 3306 3307 ! Back to usual orthonormalization for the cg and cprj_k 3308 call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,& 3309 & mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw) 3310 3311 ! Need to transfer cprj_k to cprj 3312 if(usepaw==1) then 3313 call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,iorder,isppol,& 3314 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 3315 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 3316 end if 3317 3318 ibg=ibg+my_nspinor*nband_k 3319 ibg_hist=ibg_hist+my_nspinor*nbdmix 3320 icg=icg+my_nspinor*nband_k*npw_k 3321 icg_hist=icg_hist+my_nspinor*nbdmix*npw_k 3322 3323 end do ! End big k point loop 3324 end do ! End loop over spins 3325 3326 if(istep>2)then 3327 ABI_FREE(coeffs) 3328 end if 3329 ABI_FREE(al) 3330 3331 end if ! wfmixalg>2 and istep>1 3332 3333 !DEBUG 3334 ! write(std_out,*)' wf_mixing : exit ' 3335 ! write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',& 3336 !& scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2) 3337 ! write(std_out,*)' cg(1:2,1:2)=',cg(1:2,1:2) 3338 ! write(std_out,*)' scf_history_wf%cg(1:2,1:2,1)=',scf_history_wf%cg(1:2,1:2,1) 3339 ! ABI_FREE(cg_ref) 3340 ! ABI_FREE(cprj_ref) 3341 !ENDDEBUG 3342 3343 if(usepaw==1) then 3344 call pawcprj_free(cprj_k) 3345 call pawcprj_free(cprj_kh) 3346 end if 3347 ABI_FREE(cprj_k) 3348 ABI_FREE(cprj_kh) 3349 ABI_FREE(dimcprj) 3350 ABI_FREE(mmn) 3351 ABI_FREE(smn) 3352 if(wfmixalg>2)then 3353 ABI_FREE(dotprod_res_k) 3354 ABI_FREE(dotprod_res) 3355 ABI_FREE(res_mn) 3356 end if 3357 3358 end subroutine wf_mixing