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