TABLE OF CONTENTS
- ABINIT/m_gstate
- ABINIT/outxfhist
- m_gstate/clnup1
- m_gstate/clnup2
- m_gstate/gstate
- m_gstate/pawuj_drive
- m_gstate/prtxf
- m_gstate/setup2
ABINIT/m_gstate [ Modules ]
NAME
m_gstate
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, JYR, MKV, MT, FJ, MB, DJA) 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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 ! nvtx related macro definition 22 #include "nvtx_macros.h" 23 24 module m_gstate 25 26 use defs_basis 27 use defs_rectypes 28 use m_errors 29 use m_xmpi 30 use m_abicore 31 use libxc_functionals 32 use m_exit 33 use m_crystal 34 use m_scf_history 35 use m_abimover 36 use m_wffile 37 use m_rec 38 use m_efield 39 use m_ddb 40 use m_bandfft_kpt 41 use m_invovl 42 use m_gemm_nonlop 43 use m_gemm_nonlop_gpu 44 use m_gemm_nonlop_ompgpu 45 use m_wfk 46 use m_nctk 47 use m_hdr 48 use m_ebands 49 use m_dtfil 50 use m_extfpmd 51 52 use defs_datatypes, only : pseudopotential_type, ebands_t 53 use defs_abitypes, only : MPI_type 54 use m_time, only : timab 55 use m_symtk, only : matr3inv 56 use m_io_tools, only : open_file 57 use m_occ, only : newocc, getnel 58 use m_ddb_hdr, only : ddb_hdr_type 59 use m_fstrings, only : strcat, sjoin 60 use m_geometry, only : fixsym, mkradim, metric 61 use m_kpts, only : tetra_from_kptrlatt 62 use m_kg, only : kpgio, getph 63 use m_fft, only : fourdp 64 use m_pawang, only : pawang_type 65 use m_pawrad, only : pawrad_type 66 use m_pawtab, only : pawtab_type, pawtab_print 67 use m_pawcprj, only : pawcprj_type,pawcprj_free,pawcprj_alloc, pawcprj_getdim 68 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 69 use m_abi2big, only : wvl_occ_abi2big, wvl_setngfft, wvl_setBoxGeometry 70 use m_energies, only : energies_type, energies_init 71 use m_args_gs, only : args_gs_type 72 use m_results_gs, only : results_gs_type 73 use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_free 74 use m_paw_dmft, only : init_sc_dmft,destroy_sc_dmft,print_sc_dmft,paw_dmft_type,readocc_dmft 75 use m_paw_sphharm, only : setsym_ylm 76 use m_paw_occupancies, only : initrhoij 77 use m_paw_init, only : pawinit,paw_gencond 78 use m_paw_correlations, only : pawpuxinit 79 use m_paw_uj, only : pawuj_ini,pawuj_free,pawuj_det, macro_uj_type 80 use m_data4entropyDMFT, only : data4entropyDMFT_t, data4entropyDMFT_init, data4entropyDMFT_destroy 81 use m_electronpositron, only : electronpositron_type,init_electronpositron,destroy_electronpositron, & 82 electronpositron_calctype 83 use m_scfcv, only : scfcv_t, scfcv_init, scfcv_destroy, scfcv_run 84 use m_jellium, only : jellium 85 use m_iowf, only : outwf, outresid 86 use m_outqmc, only : outqmc 87 use m_ioarr, only : ioarr,read_rhor 88 use m_inwffil, only : inwffil 89 use m_spacepar, only : setsym 90 use m_mkrho, only : mkrho, initro, prtrhomxmn 91 use m_initylmg, only : initylmg 92 use m_pspini, only : pspini 93 use m_mover, only : mover 94 use m_mpinfo, only : proc_distrb_cycle 95 use m_common, only : setup1, prteigrs, prtene 96 use m_fourier_interpol, only : transgrid 97 use m_psolver, only : psolver_kernel 98 use m_paw2wvl, only : paw2wvl, wvl_paw_free 99 use m_berryphase_new, only : init_e_field_vars,prtefield 100 use m_wvl_wfs, only : wvl_wfs_set, wvl_wfs_free, wvl_wfs_lr_copy 101 use m_wvl_rho, only : wvl_initro, wvl_mkrho 102 use m_wvl_descr_psp, only : wvl_descr_psp_set, wvl_descr_free, wvl_descr_atoms_set, wvl_descr_atoms_set_sym 103 use m_wvl_denspot, only : wvl_denspot_set, wvl_denspot_free 104 use m_wvl_projectors, only : wvl_projectors_set, wvl_projectors_free 105 use m_cgprj, only : ctocprj 106 use m_nonlop_ylm, only : nonlop_ylm_init_counters,nonlop_ylm_output_counters 107 use m_fft, only : fft_init_counters,fft_output_counters 108 use m_getghc_ompgpu, only : free_getghc_ompgpu_buffers 109 110 #if defined HAVE_GPU 111 use m_alloc_hamilt_gpu 112 #endif 113 114 #if defined(HAVE_GPU) && defined(HAVE_GPU_MARKERS) 115 use m_nvtx_data 116 #endif 117 118 #if defined HAVE_YAKL 119 use gator_mod 120 #endif 121 122 use defs_wvltypes, only : wvl_data,coulomb_operator,wvl_wf_type 123 #if defined HAVE_BIGDFT 124 use BigDFT_API, only : wvl_timing => timing,xc_init,xc_end,XC_MIXED,XC_ABINIT,& 125 local_potential_dimensions,nullify_gaussian_basis, & 126 copy_coulomb_operator,deallocate_coulomb_operator 127 #else 128 use defs_wvltypes, only : coulomb_operator 129 #endif 130 131 #if defined HAVE_LOTF 132 use defs_param_lotf, only : lotfparam_init 133 #endif 134 135 implicit none 136 137 private
ABINIT/outxfhist [ Functions ]
NAME
outxfhist
FUNCTION
read/write xfhist
COPYRIGHT
Copyright (C) 2003-2024 ABINIT group (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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
option = 1: write 2: read only nxfh 3: read xfhist response = 0: GS wavefunctions 1: RF wavefunctions natom = number of atoms in unit cell mxfh = last dimension of the xfhist array
OUTPUT
ios = error code returned by read operations
SIDE EFFECTS
nxfh = actual number of (x,f) history pairs, see xfhist array wff2 = structured info for wavefunctions xfhist(3,natom+4,2,ab_xfh%mxfh) = (x,f) history array, also including rprim and stress
SOURCE
2655 subroutine outxfhist(ab_xfh,natom,option,wff2,ios) 2656 2657 use defs_basis 2658 use m_abicore 2659 use m_abimover 2660 use m_xmpi 2661 use m_wffile 2662 use m_errors 2663 #if defined HAVE_NETCDF 2664 use netcdf 2665 #endif 2666 2667 !Arguments ------------------------------------ 2668 integer ,intent(in) :: natom,option 2669 integer ,intent(out) :: ios 2670 type(wffile_type),intent(inout) :: wff2 2671 type(ab_xfh_type),intent(inout) :: ab_xfh 2672 2673 !Local variables------------------------------- 2674 integer :: ierr,ixfh,ncid_hdr,spaceComm,xfdim2 2675 real(dp),allocatable :: xfhist_tmp(:) 2676 character(len=500) :: msg 2677 !no_abirules 2678 #if defined HAVE_NETCDF 2679 integer :: ncerr 2680 integer :: nxfh_id, mxfh_id, xfdim2_id, dim2inout_id, dimr3_id,xfhist_id 2681 integer :: nxfh_tmp,mxfh_tmp,xfdim2_tmp,dim2inout_tmp 2682 #endif 2683 2684 2685 ! ************************************************************************* 2686 2687 ncid_hdr = wff2%unwff 2688 xfdim2 = natom+4 2689 2690 ios = 0 2691 2692 !### (Option=1) Write out content of all iterations 2693 !##################################################################### 2694 if ( option == 1 ) then 2695 2696 ! Write the (x,f) history 2697 if (wff2%iomode == IO_MODE_FORTRAN) then 2698 write(unit=wff2%unwff)ab_xfh%nxfh 2699 do ixfh=1,ab_xfh%nxfh 2700 write(unit=wff2%unwff)ab_xfh%xfhist(:,:,:,ixfh) 2701 end do 2702 2703 else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then 2704 ! FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case 2705 ! if node is master 2706 write(msg, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, & 2707 & 'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.' 2708 ABI_ERROR(msg) 2709 2710 write(unit=wff2%unwff)ab_xfh%nxfh 2711 do ixfh=1,ab_xfh%nxfh 2712 write(unit=wff2%unwff)ab_xfh%xfhist(:,:,:,ixfh) 2713 end do 2714 2715 ! insert mpi broadcast here 2716 2717 else if(wff2%iomode==IO_MODE_MPI)then 2718 ABI_MALLOC(xfhist_tmp,(3*(natom+4)*2)) 2719 spaceComm=xmpi_comm_self 2720 call xderiveWRecInit(wff2,ierr) 2721 call xderiveWrite(wff2,ab_xfh%nxfh,ierr) 2722 call xderiveWRecEnd(wff2,ierr) 2723 do ixfh=1,ab_xfh%nxfh 2724 xfhist_tmp(:)=reshape(ab_xfh%xfhist(:,:,:,ixfh),(/3*(natom+4)*2/)) 2725 call xderiveWRecInit(wff2,ierr) 2726 call xderiveWrite(wff2,xfhist_tmp,3*(natom+4)*2,spaceComm,ierr) 2727 call xderiveWRecEnd(wff2,ierr) 2728 end do 2729 ABI_FREE(xfhist_tmp) 2730 2731 #if defined HAVE_NETCDF 2732 else if (wff2%iomode == IO_MODE_NETCDF) then 2733 ! check if nxfh and xfhist are defined 2734 ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id) 2735 2736 if (ncerr /= NF90_NOERR) then 2737 ! need to define everything 2738 ncerr = nf90_redef (ncid=ncid_hdr) 2739 NCF_CHECK_MSG(ncerr," outxfhist : going to define mode ") 2740 2741 ncerr = nf90_def_dim(ncid=ncid_hdr,name="dim2inout",len=2,dimid=dim2inout_id) 2742 NCF_CHECK_MSG(ncerr," outxfhist : define dim2inout") 2743 ncerr = nf90_def_dim(ncid=ncid_hdr,name="mxfh",len=ab_xfh%mxfh,dimid=mxfh_id) 2744 NCF_CHECK_MSG(ncerr," outxfhist : define mxfh") 2745 ncerr = nf90_def_dim(ncid=ncid_hdr,name="nxfh",len=ab_xfh%nxfh,dimid=nxfh_id) 2746 NCF_CHECK_MSG(ncerr," outxfhist : define nxfh") 2747 ncerr = nf90_def_dim(ncid=ncid_hdr,name="xfdim2",len=xfdim2,dimid=xfdim2_id) 2748 NCF_CHECK_MSG(ncerr," outxfhist : define xfdim2") 2749 2750 ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="dimr3",dimid=dimr3_id) 2751 NCF_CHECK_MSG(ncerr," outxfhist : inquire dimr3") 2752 2753 ! ab_xfh%xfhist(3,natom+4,2,ab_xfh%mxfh) 2754 ncerr = nf90_def_var(ncid=ncid_hdr,name="xfhist",xtype=NF90_DOUBLE,& 2755 & dimids=(/dimr3_id,xfdim2_id,dim2inout_id,mxfh_id/),varid=xfhist_id) 2756 NCF_CHECK_MSG(ncerr," outxfhist : define xfhist") 2757 2758 ! End define mode and go to data mode 2759 ncerr = nf90_enddef(ncid=ncid_hdr) 2760 NCF_CHECK_MSG(ncerr," outxfhist : enddef call ") 2761 else 2762 ! check that the dimensions are correct 2763 ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id) 2764 NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh") 2765 ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,& 2766 & len=nxfh_tmp) 2767 NCF_CHECK_MSG(ncerr," outxfhist : get nxfh") 2768 ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="xfdim2",dimid=xfdim2_id) 2769 NCF_CHECK_MSG(ncerr," outxfhist : inquire xfdim2") 2770 ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=xfdim2_id,& 2771 & len=xfdim2_tmp) 2772 NCF_CHECK_MSG(ncerr," outxfhist : get xfdim2") 2773 ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="mxfh",dimid=mxfh_id) 2774 NCF_CHECK_MSG(ncerr," outxfhist : inquire mxfh") 2775 ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=mxfh_id,& 2776 & len=mxfh_tmp) 2777 NCF_CHECK_MSG(ncerr," outxfhist : get mxfh") 2778 ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="dim2inout",dimid=dim2inout_id) 2779 NCF_CHECK_MSG(ncerr," outxfhist : inquire dim2inout") 2780 ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=dim2inout_id,& 2781 & len=dim2inout_tmp) 2782 NCF_CHECK_MSG(ncerr," outxfhist : get dim2inout") 2783 2784 ncerr = nf90_inq_varid(ncid=ncid_hdr,name="xfhist",varid=xfhist_id) 2785 NCF_CHECK_MSG(ncerr," outxfhist : inquire xfhist") 2786 2787 if (mxfh_tmp /= ab_xfh%mxfh .or. dim2inout_tmp /= 2 .or. xfdim2_tmp /= xfdim2) then 2788 write (msg,"(A)") 'outxfhist : ERROR xfhist has bad dimensions in NetCDF file. Can not re-write it.' 2789 ABI_ERROR(msg) 2790 end if 2791 2792 end if 2793 2794 ! Now fill the data 2795 ncerr = nf90_put_var(ncid=ncid_hdr,varid=xfhist_id,values=ab_xfh%xfhist) 2796 NCF_CHECK_MSG(ncerr," outxfhist : fill xfhist") 2797 2798 ! end NETCDF definition ifdef 2799 #endif 2800 end if ! end iomode if 2801 2802 ! ### (Option=2) Read in number of iterations 2803 ! ##################################################################### 2804 else if ( option == 2 ) then 2805 2806 if (wff2%iomode == IO_MODE_FORTRAN) then 2807 read(unit=wff2%unwff,iostat=ios)ab_xfh%nxfh 2808 2809 else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then 2810 ! FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case 2811 ! if node is master 2812 write(msg, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, & 2813 & 'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.' 2814 ABI_ERROR(msg) 2815 2816 read(unit=wff2%unwff,iostat=ios)ab_xfh%nxfh 2817 2818 else if (wff2%iomode == IO_MODE_MPI) then 2819 call xderiveRRecInit(wff2,ierr) 2820 call xderiveRead(wff2,ab_xfh%nxfh,ierr) 2821 call xderiveRRecEnd(wff2,ierr) 2822 2823 #if defined HAVE_NETCDF 2824 else if (wff2%iomode == IO_MODE_NETCDF) then 2825 ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id) 2826 NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh") 2827 ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,& 2828 & len=ab_xfh%nxfh) 2829 NCF_CHECK_MSG(ncerr," outxfhist : get nxfh") 2830 #endif 2831 end if 2832 2833 ! ### (Option=3) Read in iteration content 2834 ! ##################################################################### 2835 else if ( option == 3 ) then 2836 if (wff2%iomode == IO_MODE_FORTRAN) then 2837 do ixfh=1,ab_xfh%nxfhr 2838 read(unit=wff2%unwff,iostat=ios)ab_xfh%xfhist(:,:,:,ixfh) 2839 end do 2840 else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then 2841 ! FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case 2842 ! if node is master 2843 write(msg, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, & 2844 & 'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.' 2845 ABI_ERROR(msg) 2846 2847 do ixfh=1,ab_xfh%nxfhr 2848 read(unit=wff2%unwff,iostat=ios)ab_xfh%xfhist(:,:,:,ixfh) 2849 end do 2850 2851 else if (wff2%iomode == IO_MODE_MPI) then 2852 ABI_MALLOC(xfhist_tmp,(3*(natom+4)*2)) 2853 spaceComm=xmpi_comm_self 2854 do ixfh=1,ab_xfh%nxfhr 2855 call xderiveRRecInit(wff2,ierr) 2856 call xderiveRead(wff2,xfhist_tmp,3*(natom+4)*2,spaceComm,ierr) 2857 call xderiveRRecEnd(wff2,ierr) 2858 xfhist_tmp(:)=xfhist_tmp(:) 2859 end do 2860 ABI_FREE(xfhist_tmp) 2861 end if 2862 2863 ! FIXME: should this be inside the if not mpi as above for options 1 and 2? 2864 ! it is placed here because the netcdf read is a single operation 2865 #if defined HAVE_NETCDF 2866 if (wff2%iomode == IO_MODE_NETCDF) then 2867 ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id) 2868 NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh") 2869 ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,& 2870 & len=ab_xfh%nxfhr) 2871 NCF_CHECK_MSG(ncerr," outxfhist : get nxfh") 2872 2873 ncerr = nf90_inq_varid(ncid=ncid_hdr,varid=xfhist_id,name="xfhist") 2874 NCF_CHECK_MSG(ncerr," outxfhist : inquire xfhist") 2875 ncerr = nf90_get_var(ncid=ncid_hdr,varid=xfhist_id,values=ab_xfh%xfhist,& 2876 & start=(/1,1,1,1/),count=(/3,natom+4,2,ab_xfh%nxfhr/)) 2877 NCF_CHECK_MSG(ncerr," outxfhist : read xfhist") 2878 end if 2879 #endif 2880 2881 else 2882 ! write(std_out,*)' outxfhist : option ', option , ' not available ' 2883 write(msg, "(A,A,A,A,I3,A)") ch10, "outxfhist: ERROR -", ch10, & 2884 & "option ", option, " not available." 2885 ABI_ERROR(msg) 2886 end if 2887 2888 end subroutine outxfhist
m_gstate/clnup1 [ Functions ]
[ Top ] [ m_gstate ] [ Functions ]
NAME
clnup1
FUNCTION
Perform "cleanup" at end of execution of gstate routine.
INPUTS
acell(3)=length scales of primitive translations (bohr) dosdeltae=DOS delta of Energy dtset <type(dataset_type)>=all input variables in this dataset eigen(mband*nkpt*nsppol)=eigenvalues (hartree) for all bands at each k point enunit=choice for units of output eigenvalues: 0=>hartree, 1=> eV, 2=> hartree and eV fermie=fermi energy (Hartree) fermih=fermi energy for holes (Hartree) (for occopt 9) fnameabo_dos=filename of output DOS file fnameabo_eig=filename of output EIG file gred(3,natom)=d(E)/d(xred) (hartree) iatfix(3,natom)=0 if not fixed along specified direction, 1 if fixed iscf=parameter controlling scf or non-scf choice kptopt=option for the generation of k points kptns(3,nkpt)=k points in terms of recip primitive translations mband=maximum number of bands mpi_enreg=information about MPI parallelization natom=number of atoms in unit cell nband(nkpt*nsppol)=number of bands 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 nkpt=number of k points nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nstep=desired number of electron iteration steps occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point occopt=option for occupancies prtdos= if == 1, will print the density of states prtfor= if >0, will print the forces prtstm= input variable prtstm prtvol=control print volume and debugging resid(mband*nkpt*nsppol)=squared residuals for each band and k point where resid(n,k)=|<C(n,k)|(H-e(n,k))|C(n,k)>|^2 rhor(nfft,nspden)=electron density (electrons/bohr^3) rprimd(3,3)=dimensional real space primitive translations (bohr) tphysel="physical" electronic temperature with FD occupations tsmear=smearing energy or temperature (if metal) vxcavg=average of vxc potential wtk(nkpt)=real(dp) array of k-point weights xred(3,natom)=reduced atomic coordinates
OUTPUT
(only print and write to disk)
SOURCE
1976 subroutine clnup1(acell,dtset,eigen,fermie,fermih, fnameabo_dos,fnameabo_eig,gred,& 1977 mpi_enreg,nfft,ngfft,occ,prtfor, resid,rhor,rprimd,vxcavg,xred) 1978 1979 !Arguments ------------------------------------ 1980 !scalars 1981 integer,intent(in) :: nfft, prtfor 1982 real(dp),intent(in) :: fermie, fermih, vxcavg 1983 character(len=*),intent(in) :: fnameabo_dos,fnameabo_eig 1984 type(dataset_type),intent(in) :: dtset 1985 type(MPI_type),intent(in) :: mpi_enreg 1986 !arrays 1987 integer,intent(in) :: ngfft(18) 1988 real(dp),intent(in) :: acell(3) 1989 real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 1990 real(dp),intent(in) :: gred(3,dtset%natom) 1991 real(dp),intent(in) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol) 1992 real(dp),intent(in) :: rhor(nfft,dtset%nspden) 1993 real(dp),intent(in) :: rprimd(3,3) 1994 real(dp),intent(in) :: xred(3,dtset%natom) 1995 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 1996 1997 !Local variables------------------------------- 1998 !scalars 1999 integer,parameter :: master=0 2000 integer :: comm,iatom,ii,iscf_dum,iwfrc,me,nnonsc,option,unitdos 2001 real(dp) :: entropy,grmax,grsum,maxocc,nelect,tolwf,ucvol 2002 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 2003 character(len=500) :: msg 2004 character(len=fnlen) filename 2005 !arrays 2006 real(dp),allocatable :: doccde(:) 2007 2008 ! **************************************************************** 2009 2010 comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm) 2011 2012 if(dtset%prtstm==0)then ! Write reduced coordinates xred 2013 write(msg, '(a,i5,a)' )' reduced coordinates (array xred) for',dtset%natom,' atoms' 2014 call wrtout(ab_out, msg) 2015 do iatom=1,dtset%natom 2016 write(msg, '(1x,3f20.12)' ) xred(:,iatom) 2017 call wrtout(ab_out, msg) 2018 end do 2019 end if 2020 2021 !Write reduced gradients if iscf > 0 and dtset%nstep>0 and 2022 if (dtset%iscf>=0.and.dtset%nstep>0.and.dtset%prtstm==0) then 2023 2024 ! Compute absolute maximum and root mean square value of gradients 2025 grmax=0.0_dp 2026 grsum=0.0_dp 2027 do iatom=1,dtset%natom 2028 do ii=1,3 2029 ! To be activated in v5.5 2030 ! grmax=max(grmax,abs(gred(ii,iatom))) 2031 grmax=max(grmax,gred(ii,iatom)) 2032 grsum=grsum+gred(ii,iatom)**2 2033 end do 2034 end do 2035 grsum=sqrt(grsum/dble(3*dtset%natom)) 2036 2037 write(msg, '(1x,a,1p,e12.4,a,e12.4,a)' )'rms dE/dt=',grsum,'; max dE/dt=',grmax,'; dE/dt below (all hartree)' 2038 call wrtout(ab_out, msg) 2039 do iatom=1,dtset%natom 2040 write(msg, '(i5,1x,3f20.12)' ) iatom,gred(1:3,iatom) 2041 call wrtout(ab_out, msg) 2042 end do 2043 2044 end if 2045 2046 if(dtset%prtstm==0)then 2047 2048 ! Compute and write out dimensional cartesian coords and forces: 2049 call wrtout(ab_out,' ') 2050 2051 ! (only write forces if iscf > 0 and dtset%nstep>0) 2052 if (dtset%iscf<0.or.dtset%nstep<=0.or.prtfor==0) then 2053 iwfrc=0 2054 else 2055 iwfrc=1 2056 end if 2057 2058 call prtxf(gred,dtset%iatfix,ab_out,iwfrc,dtset%natom,rprimd,xred) 2059 2060 ! Write length scales 2061 write(msg, '(1x,a,3f16.12,a)' )'length scales=',acell,' bohr' 2062 call wrtout(ab_out, msg) 2063 write(msg, '(14x,a,3f16.12,a)' )'=',Bohr_Ang*acell(1:3),' angstroms' 2064 call wrtout(ab_out, msg) 2065 2066 end if 2067 2068 option=1; nnonsc=0; tolwf=0.0_dp 2069 2070 if(dtset%iscf<0 .and. dtset%iscf/=-3)option=3 2071 iscf_dum=dtset%iscf 2072 if(dtset%nstep==0)iscf_dum=-1 2073 2074 if(dtset%tfkinfunc==0)then 2075 if (me == master) then 2076 call prteigrs(eigen,dtset%enunit,fermie,fermih,fnameabo_eig,ab_out,& 2077 & iscf_dum,dtset%kptns,dtset%kptopt,dtset%mband,& 2078 & dtset%nband,dtset%nbdbuf,dtset%nkpt,nnonsc,dtset%nsppol,occ,& 2079 & dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwf,& 2080 & vxcavg,dtset%wtk) 2081 call prteigrs(eigen,dtset%enunit,fermie,fermih,fnameabo_eig,std_out,& 2082 & iscf_dum,dtset%kptns,dtset%kptopt,dtset%mband,& 2083 & dtset%nband,dtset%nbdbuf,dtset%nkpt,nnonsc,dtset%nsppol,occ,& 2084 & dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwf,& 2085 & vxcavg,dtset%wtk) 2086 end if 2087 2088 #if defined HAVE_NETCDF 2089 if (dtset%prteig==1 .and. me == master) then 2090 filename=trim(fnameabo_eig)//'.nc' 2091 call write_eig(eigen,fermie,filename,dtset%kptns,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol) 2092 end if 2093 #endif 2094 end if 2095 2096 !Compute and print location of maximal and minimal density 2097 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2098 call prtrhomxmn(std_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol) 2099 if( dtset%prtvol>1)then 2100 call prtrhomxmn(ab_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol) 2101 end if 2102 2103 !If needed, print DOS (unitdos is closed in getnel, occ is not changed if option == 2 2104 if (dtset%prtdos==1 .and. me == master) then 2105 if (open_file(fnameabo_dos,msg, newunit=unitdos, status='unknown', action="write", form='formatted') /= 0) then 2106 ABI_ERROR(msg) 2107 end if 2108 rewind(unitdos) 2109 maxocc=two/(dtset%nspinor*dtset%nsppol) ! Will not work in the fixed moment case 2110 option=2 2111 ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 2112 call getnel(doccde,dtset%dosdeltae,eigen,entropy,fermie,fermih,& 2113 & maxocc,dtset%mband,dtset%nband,nelect,dtset%nkpt,& 2114 & dtset%nsppol,occ,dtset%occopt,option,dtset%tphysel,& 2115 & dtset%tsmear,unitdos,dtset%wtk,1,dtset%nband(1))!CP: added 1, nband(1) to fit new definition of getnel; parameters only used if 2116 ABI_FREE(doccde) 2117 end if 2118 2119 end subroutine clnup1
m_gstate/clnup2 [ Functions ]
[ Top ] [ m_gstate ] [ Functions ]
NAME
clnup2
FUNCTION
Perform more "cleanup" after completion of iterations. This subroutine prints out more breakdown of force information, shifts of atomic positions, and stresses.
INPUTS
gred(3,natom)=d(E_total)/d(xred) derivatives (hartree) grchempottn(3,natom)=d(E_chempot)/d(xred) derivatives (hartree) grewtn(3,natom)=d(E_Ewald)/d(xred) derivatives (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree) grxc(3,natom)=d(Exc)/d(xred) derivatives (0 without core charges) iscf=parameter controlling scf or non-scf iterations natom=number of atoms in unit cell ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc n1xccc=dimension of xccc1d ; 0 if no XC core correction is used prtfor= >0 if forces have to be printed (0 otherwise) prtstr= >0 if stresses have to be printed (0 otherwise) prtvol=control print volume and debugging output start(3,natom)=starting coordinates in terms of real space primitive translations strten(6)=components of the stress tensor (hartree/bohr^3) synlgr(3,natom)=d(E_nlpsp)/d(xred) derivatives (hartree) xred(3,natom)=final coordinates in terms of primitive translations
OUTPUT
(only print)
SOURCE
2318 subroutine clnup2(n1xccc,gred,grchempottn,gresid,grewtn,grvdw,grxc,iscf,natom,ngrvdw,& 2319 & prtfor,prtstr,prtvol,start,strten,synlgr,xred) 2320 2321 !Arguments ------------------------------------ 2322 !scalars 2323 integer,intent(in) :: iscf,n1xccc,natom,ngrvdw,prtfor,prtstr,prtvol 2324 !arrays 2325 real(dp),intent(in) :: gred(3,natom),grchempottn(3,natom),gresid(3,natom) 2326 real(dp),intent(in) :: grewtn(3,natom),grvdw(3,ngrvdw) 2327 real(dp),intent(in) :: grxc(3,natom),start(3,natom),strten(6),synlgr(3,natom) 2328 real(dp),intent(in) :: xred(3,natom) 2329 2330 !Local variables------------------------------- 2331 character(len=*), parameter :: format01020 ="(i5,1x,3f20.12)" 2332 !scalars 2333 integer :: iatom,mu 2334 real(dp) :: devsqr,grchempot2 2335 character(len=500) :: msg 2336 integer :: units(2) 2337 2338 ! ************************************************************************* 2339 2340 !write(std_out,*)' clnup2 : enter ' 2341 2342 !Only print additional info for scf calculations 2343 if (iscf>=0) then 2344 2345 if(prtvol >= 10 .and. prtfor > 0) then 2346 write(msg, '(a,10x,a)' ) ch10, '===> extra information on forces <===' 2347 call wrtout(ab_out,msg) 2348 2349 call wrtout(ab_out, ' ewald contribution to reduced grads') 2350 do iatom=1,natom 2351 write(msg,format01020) iatom,(grewtn(mu,iatom),mu=1,3) 2352 call wrtout(ab_out,msg) 2353 end do 2354 2355 grchempot2=sum(grchempottn(:,:)**2) 2356 if(grchempot2>tol16)then 2357 call wrtout(ab_out, ' chemical potential contribution to reduced grads') 2358 do iatom=1,natom 2359 write(msg,format01020) iatom,(grchempottn(mu,iatom),mu=1,3) 2360 call wrtout(ab_out,msg) 2361 end do 2362 end if 2363 2364 call wrtout(ab_out,' nonlocal contribution to red. grads') 2365 do iatom=1,natom 2366 write(msg,format01020) iatom,(synlgr(mu,iatom),mu=1,3) 2367 call wrtout(ab_out,msg) 2368 end do 2369 2370 call wrtout(ab_out, ' local psp contribution to red. grads') 2371 if (n1xccc /= 0) then 2372 do iatom=1,natom 2373 write(msg,format01020) iatom,gred(:,iatom) - & 2374 (grewtn(:,iatom)+grchempottn(:,iatom)+synlgr(:,iatom)+grxc(:,iatom)+gresid(:,iatom)) 2375 call wrtout(ab_out,msg) 2376 end do 2377 else 2378 do iatom=1,natom 2379 write(msg,format01020) iatom,gred(:,iatom) - & 2380 (grewtn(:,iatom)+grchempottn(:,iatom)+synlgr(:,iatom)+gresid(:,iatom)) 2381 call wrtout(ab_out,msg) 2382 end do 2383 end if 2384 2385 if (n1xccc /= 0) then 2386 call wrtout(ab_out,' core charge xc contribution to reduced grads') 2387 do iatom=1,natom 2388 write(msg,format01020) iatom,(grxc(mu,iatom),mu=1,3) 2389 call wrtout(ab_out,msg) 2390 end do 2391 end if 2392 2393 if (ngrvdw == natom) then 2394 call wrtout(ab_out,' Van der Waals DFT-D contribution to reduced grads') 2395 do iatom=1,natom 2396 write(msg,format01020) iatom,(grvdw(mu,iatom),mu=1,3) 2397 call wrtout(ab_out,msg) 2398 end do 2399 end if 2400 2401 call wrtout(ab_out,' residual contribution to red. grads') 2402 do iatom=1,natom 2403 write(msg,format01020) iatom,(gresid(mu,iatom),mu=1,3) 2404 call wrtout(ab_out,msg) 2405 end do 2406 2407 end if 2408 2409 ! Compute mean squared deviation from starting coords 2410 devsqr=zero 2411 do iatom=1,natom 2412 do mu=1,3 2413 devsqr=devsqr+(xred(mu,iatom)-start(mu,iatom))**2 2414 end do 2415 end do 2416 2417 ! When shift is nonnegligible then print values 2418 if (devsqr>1.d-14) then 2419 write(msg, '(a,1p,e12.4,3x,a)' )' rms coord change=',sqrt(devsqr/dble(3*natom)),'atom, delta coord (reduced):' 2420 call wrtout(ab_out,msg) 2421 do iatom=1,natom 2422 write(msg, '(1x,i5,2x,3f20.12)' ) iatom, (xred(mu,iatom)-start(mu,iatom),mu=1,3) 2423 call wrtout(ab_out,msg) 2424 end do 2425 end if 2426 2427 if (prtstr > 0) then 2428 ! Write out stress results 2429 units = [std_out, ab_out] 2430 write(msg, '(a,a)' ) ch10,' Cartesian components of stress tensor (hartree/bohr^3)' 2431 call wrtout(units, msg) 2432 2433 write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4) 2434 call wrtout(units, msg) 2435 write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5) 2436 call wrtout(units, msg) 2437 write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6) 2438 call wrtout(units, msg) 2439 2440 ! Also output the pressure (minus one third the trace of the stress tensor). 2441 write(msg, '(a,a,es12.4,a)' ) ch10,& 2442 '-Cartesian components of stress tensor (GPa) [Pressure=',& 2443 -(strten(1)+strten(2)+strten(3))*HaBohr3_GPa/3.0_dp,' GPa]' 2444 call wrtout(units, msg) 2445 2446 write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) & 2447 '- sigma(1 1)=',strten(1)*HaBohr3_GPa,& 2448 ' sigma(3 2)=',strten(4)*HaBohr3_GPa 2449 call wrtout(units, msg) 2450 write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) & 2451 '- sigma(2 2)=',strten(2)*HaBohr3_GPa,& 2452 ' sigma(3 1)=',strten(5)*HaBohr3_GPa 2453 call wrtout(units, msg) 2454 write(msg, '(a,1p,e16.8,a,1p,e16.8)' ) & 2455 '- sigma(3 3)=',strten(3)*HaBohr3_GPa,& 2456 ' sigma(2 1)=',strten(6)*HaBohr3_GPa 2457 call wrtout(units, msg) 2458 end if 2459 2460 end if ! iscf > 0 2461 2462 !write(std_out,*)' clnup2 : exit ' 2463 2464 end subroutine clnup2
m_gstate/gstate [ Functions ]
[ Top ] [ m_gstate ] [ Functions ]
NAME
gstate
FUNCTION
Primary routine for conducting DFT calculations by CG minimization.
INPUTS
args_gs<type(args_gs_type)>=various input arguments for the GS calculation Possibly different from dtset codvsn=code version cpui=initial CPU time itimimage_gstate=counter for calling do loop
OUTPUT
npwtot(nkpt) = total number of plane waves at each k point results_gs <type(results_gs_type)>=results (energy and its components, forces and its components, the stress tensor) of a ground-state computation
SIDE EFFECTS
acell(3)=unit cell length scales (bohr) dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables in this dataset | mband =maximum number of bands (IN) | mgfft =maximum single fft dimension (IN) | mkmem =number of k points treated by this processor (IN) | mpw =maximum number of planewaves in basis sphere (large number) (IN) | natom =number of atoms in unit cell (IN) | nfft =(effective) number of FFT grid points (for this processor) (IN) | nkpt =number of k points (IN) | nspden=number of spin-density components (IN) | nsppol=number of channels for spin-polarization (1 or 2) (IN) | nsym =number of symmetry elements in space group iexit= exit flag initialized= 0 for the first GS calculation (not initialized), else 1 mpi_enreg=MPI-parallelisation information (some already initialized, some others to be initialized here) occ(mband*nkpt*nsppol) = occupation number for each band and k pawang <type(pawang_type)>=paw angular mesh 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 Before entering the first time in gstate, a significant part of psps has been initialized : the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid, ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp All the remaining components of psps are to be initialized in the call to pspini . The next time the code enters gstate, psps might be identical to the one of the previous dtset, in which case, no reinitialisation is scheduled in pspini.f . rprim(3,3)=dimensionless real space primitive translations scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles vel(3,natom)=value of velocity vel_cell(3,3)=value of cell parameters velocity wvl <type(wvl_data)>=all wavelets data xred(3,natom) = reduced atomic coordinates
NOTES
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.
In case of wavelets:
--------------------
- Only the usual FFT grid (defined by wvl_crmult) is used.
It is defined by nfft, ngfft, mgfft, ... This is strictly not
an FFT grid since its dimensions are not suited for FFTs. They are
defined by wvl_setngfft().
For compatibility reasons, (nfftf,ngfftf,mgfftf)
are set equal to (nfft,ngfft,mgfft) in that case.
TODO
Not yet possible to use restartxf in parallel when localrdwf==0
SOURCE
240 subroutine gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,initialized,& 241 & itimimage_gstate,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,& 242 & psps,results_gs,rprim,scf_history,vel,vel_cell,wvl,xred) 243 244 !Arguments ------------------------------------ 245 !scalars 246 integer,intent(inout) :: iexit,initialized 247 integer,intent(in) :: itimimage_gstate 248 real(dp),intent(in) :: cpui 249 character(len=8),intent(in) :: codvsn 250 type(MPI_type),intent(inout) :: mpi_enreg 251 type(args_gs_type),intent(in) :: args_gs 252 type(datafiles_type),intent(inout) :: dtfil 253 type(dataset_type),intent(inout) :: dtset 254 type(pawang_type),intent(inout) :: pawang 255 type(pseudopotential_type),intent(inout) :: psps 256 type(results_gs_type),intent(inout) :: results_gs 257 type(scf_history_type),target,intent(inout) :: scf_history 258 type(wvl_data),intent(inout) :: wvl 259 !arrays 260 integer,intent(out) :: npwtot(dtset%nkpt) 261 real(dp),intent(inout) :: acell(3),occ(dtset%mband*dtset%nkpt*dtset%nsppol) 262 real(dp),intent(inout) :: rprim(3,3),vel(3,dtset%natom),vel_cell(3,3),xred(3,dtset%natom) 263 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 264 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 265 266 !Local variables------------------------------- 267 !Define file format for different type of files. Presently, 268 !only one file format is supported for each type of files, but this might change soon ... 269 !2 for wavefunction file, new format (version 2.0 and after) (fform) NOT USED 270 !52 for density rho(r) (fformr) 271 !102 for potential V(r) file. (fformv) NOT USED 272 !scalars 273 logical :: compute_cprj 274 integer,parameter :: formeig=0, level=101, response=0 ,cplex1=1, master=0, itime0=0 275 integer :: ndtpawuj=0 ! Cannot use parameter because scfargs points to this! Have to get rid of pointers to scalars! 276 #if defined HAVE_BIGDFT 277 integer :: icoulomb 278 #endif 279 integer :: accessfil,ask_accurate,bantot,choice,comm_psp,fform 280 integer :: gnt_option,gscase,iatom,idir,ierr,ii,indx,jj,kk,ios,iorder_cprj,itypat 281 integer :: ixfh,mband_cprj,mcg,mcprj,me,mgfftf,mpert,mu,my_natom,my_nspinor 282 integer :: nband_k,nbandtot,nblok,ncprj,ncpgr,nfftf,nfftot,npwmin 283 integer :: openexit,option,optorth,psp_gencond,conv_retcode 284 integer :: pwind_alloc,rdwrpaw,comm,tim_mkrho,use_sc_dmft 285 integer :: cnt,spin,band,ikpt,usecg,usecprj,ylm_option 286 real(dp) :: cpus,ecore,ecut_eff,ecutdg_eff,etot,fermie,fermih 287 real(dp) :: gsqcut_eff,gsqcut_shp,gsqcutc_eff,hyb_range_fock,residm,ucvol 288 logical :: read_wf_or_den,has_to_init,call_pawinit,write_wfk 289 logical :: is_dfpt=.false.,wvlbigdft=.false. 290 character(len=500) :: msg 291 character(len=fnlen) :: dscrpt,filnam,wfkfull_path 292 real(dp) :: fatvshift 293 type(crystal_t) :: cryst 294 type(ebands_t) :: bstruct,ebands 295 type(efield_type) :: dtefield 296 type(electronpositron_type),pointer :: electronpositron 297 type(hdr_type) :: hdr, hdr_den, hdr_kfull 298 type(extfpmd_type),pointer :: extfpmd => null() 299 type(macro_uj_type) :: dtpawuj(1) 300 type(paw_dmft_type) :: paw_dmft 301 type(pawfgr_type) :: pawfgr 302 type(recursion_type) ::rec_set 303 type(wffile_type) :: wff1,wffnew,wffnow 304 type(ab_xfh_type) :: ab_xfh 305 type(ddb_type) :: ddb 306 type(ddb_hdr_type) :: ddb_hdr 307 type(scfcv_t) :: scfcv_args 308 !arrays 309 integer :: itimes(2),ngfft(18),ngfftf(18) 310 integer,allocatable :: atindx(:),atindx1(:),indsym(:,:,:),dimcprj_srt(:) 311 integer,allocatable :: irrzon(:,:,:),kg(:,:),nattyp(:),symrec(:,:,:) 312 integer,allocatable,target :: npwarr(:) 313 integer,pointer :: npwarr_(:),pwind(:,:,:) 314 real(dp) :: efield_band(3),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3) 315 real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3),tsec(2) 316 real(dp),allocatable :: doccde(:) 317 real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phnons(:,:,:),resid(:),rhowfg(:,:) 318 real(dp),allocatable :: rhowfr(:,:),spinat_dum(:,:),start(:,:),work(:) 319 real(dp),allocatable :: ylm(:,:),ylmgr(:,:,:) 320 real(dp),ABI_CONTIGUOUS pointer :: cg(:,:) => null() 321 real(dp),pointer :: eigen(:),pwnsfac(:,:),rhog(:,:),rhor(:,:) 322 real(dp),pointer :: taug(:,:),taur(:,:),xred_old(:,:) 323 type(pawrhoij_type),pointer :: pawrhoij(:) 324 type(coulomb_operator) :: kernel_dummy 325 type(pawcprj_type),allocatable :: cprj(:,:) 326 327 ! *********************************************************************** 328 329 DBG_ENTER("COLL") 330 331 call timab(1232,1,tsec) 332 call timab(1211,3,tsec) 333 334 !########################################################### 335 !### 01. Initializations XML, MPI, WVL, etc 336 337 !Init MPI data 338 comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm) 339 340 !Set up MPI information from the dataset 341 my_natom=mpi_enreg%my_natom 342 343 !Set up information when wavelets are in use 344 if (dtset%usewvl == 1) then 345 346 ! If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed 347 wvlbigdft=(dtset%wvl_bigdft_comp==1) 348 349 ! Default value, to be set-up elsewhere. 350 wvl%descr%h(:) = dtset%wvl_hgrid 351 352 #if defined HAVE_BIGDFT 353 wvl%descr%paw%usepaw=psps%usepaw 354 wvl%descr%paw%natom=dtset%natom 355 #endif 356 357 ! We set the atom-related internal wvl variables. 358 call wvl_descr_atoms_set(acell, dtset%icoulomb, dtset%natom, & 359 & dtset%ntypat, dtset%typat, wvl%descr) 360 if(dtset%usepaw==0) then 361 ! nullify PAW proj_G in NC case: 362 #if defined HAVE_BIGDFT 363 ABI_MALLOC(wvl%projectors%G,(dtset%ntypat)) 364 do itypat=1,dtset%ntypat 365 call nullify_gaussian_basis(wvl%projectors%G(itypat)) 366 end do 367 #endif 368 end if 369 370 wvl%descr%exctxpar = "OP2P" 371 end if 372 373 if (me == master .and. dtset%prtxml == 1) then 374 ! gstate() will handle a dataset, so we output the dataSet markup. 375 write(ab_xml_out, "(A)") ' <dataSet>' 376 end if 377 378 !Define FFT grid(s) sizes (be careful !) 379 !See NOTES in the comments at the beginning of this file. 380 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf) 381 382 !Structured debugging if prtvol==-level 383 if(dtset%prtvol==-level)then 384 write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' gstate : enter , debug mode ' 385 call wrtout(std_out, msg) 386 end if 387 388 !########################################################### 389 !### 02. Calls setup1, kpgio, initylmg 390 391 ecore=zero 392 results_gs%pel(1:3) =zero 393 results_gs%grchempottn(:,:)=zero 394 results_gs%grewtn(:,:)=zero 395 !MT Feb 2012: I dont know why but grvdw has to be allocated 396 !when using BigDFT to ensure success on inca_gcc44_sdebug 397 if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).or.dtset%usewvl==1) then 398 results_gs%ngrvdw=dtset%natom 399 if (allocated(results_gs%grvdw)) then 400 ABI_FREE(results_gs%grvdw) 401 end if 402 ABI_MALLOC(results_gs%grvdw,(3,dtset%natom)) 403 results_gs%grvdw(:,:)=zero 404 end if 405 call energies_init(results_gs%energies) 406 407 !Set up for iterations 408 call setup1(acell,bantot,dtset,& 409 ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,& 410 ngfftf,ngfft,dtset%nkpt,dtset%nsppol,& 411 response,rmet,rprim,rprimd,ucvol,psps%usepaw) 412 413 !In some cases (e.g. getcell/=0), the plane wave vectors have 414 ! to be generated from the original simulation cell 415 rprimd_for_kg=rprimd 416 if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=args_gs%rprimd_orig 417 if (dtset%optcell/=0.and.dtset%imgmov/=0) rprimd_for_kg=args_gs%rprimd_orig 418 call matr3inv(rprimd_for_kg,gprimd_for_kg) 419 gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg) 420 421 !Set up the basis sphere of planewaves 422 ABI_MALLOC(npwarr,(dtset%nkpt)) 423 if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then 424 ABI_MALLOC(kg,(3,dtset%mpw*dtset%mkmem)) 425 call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg, & 426 & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,& 427 & dtset%mpw,npwarr,npwtot,dtset%nsppol) 428 call bandfft_kpt_init1(bandfft_kpt,dtset%istwfk,kg,dtset%mgfft,dtset%mkmem,mpi_enreg,& 429 & dtset%mpw,dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,gpu_option=dtset%gpu_option) 430 else 431 ABI_MALLOC(kg,(0,0)) 432 npwarr(:) = 0 433 npwtot(:) = 0 434 end if 435 436 if((dtset%wfoptalg == 1 .or. dtset%wfoptalg == 111) .and. psps%usepaw == 1) then 437 call init_invovl(dtset%nkpt) 438 end if 439 440 ! Handling GEMM nonlop use 441 ! Not enabled by default for CPU and CUDA implementations 442 ! Enabled if using OpenMP GPU offload 443 gemm_nonlop_use_gemm = .false. 444 445 ! OpenMP GPU offload case (GEMM nonlop used by default) 446 if(dtset%gpu_option == ABI_GPU_OPENMP) then 447 gemm_nonlop_use_gemm = .true. 448 call init_gemm_nonlop_ompgpu(dtset%nkpt) 449 else if(dtset%use_gemm_nonlop == 1) then 450 gemm_nonlop_use_gemm = .true. 451 ! CUDA & Kokkos case (same routine used) 452 if(dtset%gpu_option == ABI_GPU_LEGACY .or. dtset%gpu_option == ABI_GPU_KOKKOS) then 453 call init_gemm_nonlop(dtset%nkpt) 454 call init_gemm_nonlop_gpu(dtset%nkpt) 455 ! CPU case 456 else if(dtset%gpu_option == ABI_GPU_DISABLED) then 457 call init_gemm_nonlop(dtset%nkpt) 458 end if 459 end if 460 461 gemm_nonlop_is_distributed = .false. 462 if(dtset%gpu_nl_distrib == 1) gemm_nonlop_is_distributed = .true. 463 464 !Set up the Ylm for each k point 465 if ( dtset%tfkinfunc /= 2) then 466 ABI_MALLOC(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)) 467 ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)) 468 if (psps%useylm==1) then 469 ylm_option=0 470 if (dtset%prtstm==0.and.dtset%iscf>0.and.dtset%positron/=1) ylm_option=1 ! compute gradients of YLM 471 if (dtset%berryopt==4 .and. dtset%optstress /= 0 .and. psps%usepaw==1) ylm_option = 1 ! compute gradients of YLM 472 call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,& 473 & psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,& 474 & npwarr,dtset%nsppol,ylm_option,rprimd,ylm,ylmgr) 475 end if 476 else 477 ABI_MALLOC(ylm,(0,0)) 478 ABI_MALLOC(ylmgr,(0,0,0)) 479 end if 480 481 !SCF history management (allocate it at first call) 482 if (initialized==0) then 483 ! This call has to be done before any use of SCF history 484 usecg=0 485 if(dtset%extrapwf>0 .or. dtset%imgwfstor==1)usecg=1 486 call scf_history_init(dtset,mpi_enreg,usecg,scf_history) 487 end if 488 has_to_init=(initialized==0.or.scf_history%history_size<0) 489 490 call timab(1211,2,tsec) 491 call timab(1212,3,tsec) 492 493 !########################################################### 494 !### 03. Calls pspini 495 496 !Open and read pseudopotential files 497 comm_psp=mpi_enreg%comm_cell;if (dtset%usewvl==1) comm_psp=mpi_enreg%comm_wvl 498 if (dtset%nimage>1) psps%mixalch(:,:)=args_gs%mixalch(:,:) ! mixalch can evolve for some image algos 499 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,psps,rprimd,comm_mpi=comm_psp) 500 call timab(1212,2,tsec) 501 call timab(1211,3,tsec) 502 503 !In case of isolated computations, ecore must set to zero 504 !because its contribution is counted in the ewald energy as the ion-ion interaction. 505 if (dtset%icoulomb == 1) ecore = zero 506 507 !WVL - Now that psp data are available, we compute rprimd, acell... from the atomic positions. 508 if (dtset%usewvl == 1) then 509 call wvl_descr_psp_set(trim(dtfil%filnam_ds(3))//"_OCCUP",dtset%nsppol,psps,dtset%spinat,wvl%descr) 510 call wvl_setBoxGeometry(dtset%prtvol, psps%gth_params%radii_cf, rprimd, xred, & 511 & wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult) 512 call mkradim(acell,rprim,rprimd) 513 rprimd_for_kg=rprimd 514 call wvl_denspot_set(wvl%den, psps%gth_params, dtset%ixc, dtset%natom, dtset%nsppol, rprimd, & 515 & wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult, mpi_enreg%comm_wvl, xred) 516 ! TODO: to be moved in a routine. 517 #if defined HAVE_BIGDFT 518 if (wvl%descr%atoms%astruct%geocode == "F") then 519 icoulomb = 1 520 else if (wvl%descr%atoms%astruct%geocode == "S") then 521 icoulomb = 2 522 else 523 icoulomb = 0 524 end if 525 ! calculation of the Poisson kernel anticipated to reduce memory peak for small systems 526 call psolver_kernel( wvl%den%denspot%dpbox%hgrids, 1, icoulomb, mpi_enreg%me_wvl, wvl%den%denspot%pkernel , & 527 & mpi_enreg%comm_wvl, wvl%den%denspot%dpbox%ndims, mpi_enreg%nproc_wvl, dtset%nscforder) 528 nullify(wvl%den%denspot%pkernelseq%kernel) 529 !call copy_coulomb_operator(wvl%den%denspot%pkernel,wvl%den%denspot%pkernelseq, "gstate") 530 ! Associate the denspot distribution into mpi_enreg. 531 mpi_enreg%nscatterarr => wvl%den%denspot%dpbox%nscatterarr 532 mpi_enreg%ngatherarr => wvl%den%denspot%dpbox%ngatherarr 533 mpi_enreg%ngfft3_ionic = wvl%den%denspot%dpbox%n3pi 534 call wvl_setngfft(mpi_enreg%me_wvl, dtset%mgfft, dtset%nfft, & 535 & dtset%ngfft, mpi_enreg%nproc_wvl, wvl%den%denspot%dpbox%ndims(1), & 536 & wvl%den%denspot%dpbox%ndims(2),wvl%den%denspot%dpbox%ndims(3),& 537 & wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl, 1)) 538 #endif 539 nfftf = dtset%nfft 540 mgfftf = dtset%mgfft 541 ngfftf(:) = dtset%ngfft(:) 542 ! Recalculate gprimd 543 call matr3inv(rprimd,gprimd) 544 ! PAW section 545 if(psps%usepaw==1) then 546 ! Reinitialize Pawfgr with new values of ngfft 547 call pawfgr_destroy(pawfgr) 548 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf) 549 ! fill wvl objects from paw objects 550 ! wvl%descr%paw%usepaw=dtset%usepaw 551 call paw2wvl(pawtab,wvl%projectors,wvl%descr) 552 end if 553 else if (dtset%icoulomb /= 0) then 554 #if defined HAVE_BIGDFT 555 if (dtset%ixc < 0) then 556 call xc_init(wvl%den%denspot%xc, dtset%ixc, XC_MIXED, dtset%nsppol) 557 else 558 call xc_init(wvl%den%denspot%xc, dtset%ixc, XC_ABINIT, dtset%nsppol) 559 end if 560 #endif 561 end if 562 563 !Initialize band structure datatype 564 if (dtset%paral_kgb/=0) then ! We decide to store total npw in bstruct, 565 ABI_MALLOC(npwarr_,(dtset%nkpt)) 566 npwarr_(:)=npwarr(:) 567 call xmpi_sum(npwarr_,mpi_enreg%comm_bandfft,ierr) 568 else 569 npwarr_ => npwarr 570 end if 571 572 bstruct = ebands_from_dtset(dtset, npwarr_) 573 574 if (dtset%paral_kgb/=0) then 575 ABI_FREE(npwarr_) 576 end if 577 nullify(npwarr_) 578 579 !Initialize PAW atomic occupancies 580 if (scf_history%history_size>=0) then 581 pawrhoij => scf_history%pawrhoij_last 582 else 583 ABI_MALLOC(pawrhoij,(my_natom*psps%usepaw)) 584 end if 585 if (psps%usepaw==1.and.has_to_init) then 586 call initrhoij(dtset%pawcpxocc,dtset%lexexch,& 587 & dtset%lpawu,my_natom,dtset%natom,dtset%nspden,dtset%nspinor,dtset%nsppol,& 588 & dtset%ntypat,pawrhoij,dtset%pawspnorb,pawtab,cplex1,dtset%spinat,dtset%typat,& 589 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 590 end if 591 592 !Initialize header 593 gscase=0 594 call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr,& 595 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 596 597 !Clean band structure datatype (should use it more in the future !) 598 call ebands_free(bstruct) 599 600 !Update header, with evolving variables, when available 601 !Here, rprimd, xred and occ are available 602 etot=hdr%etot ; fermie=hdr%fermie ; fermih=hdr%fermih ; residm=hdr%residm 603 call hdr%update(bantot,etot,fermie,fermih,& 604 residm,rprimd,occ,pawrhoij,xred,args_gs%amu,& 605 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 606 607 ! PW basis set: test if the problem is ill-defined. 608 if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then 609 npwmin=minval(hdr%npwarr(:)) 610 if (dtset%mband > npwmin) then 611 ! No way we can solve the problem. Abort now! 612 write(msg,"(2(a,i0),4a)")& 613 "Number of bands nband= ",dtset%mband," > number of planewaves npw= ",npwmin,ch10,& 614 "The number of eigenvectors cannot be greater that the size of the Hamiltonian!",ch10,& 615 "Action: decrease nband or, alternatively, increase ecut" 616 if (dtset%ionmov/=23) then 617 ABI_ERROR(msg) 618 else 619 ABI_WARNING(msg) 620 end if 621 622 else if (dtset%mband >= 0.9 * npwmin) then 623 ! Warn the user 624 write(msg,"(a,i0,a,f6.1,4a)")& 625 & "Number of bands nband= ",dtset%mband," >= 0.9 * maximum number of planewaves= ",0.9*npwmin,ch10,& 626 & "The problem is ill-defined and the GS algorithm will show numerical instabilities!",ch10,& 627 & "Assume experienced user. Execution will continue." 628 ABI_WARNING(msg) 629 end if 630 end if 631 632 !########################################################### 633 !### 04. Symmetry operations when nsym>1 634 635 !Do symmetry stuff only for nsym>1 636 if (dtset%usewvl == 0) then 637 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 638 else 639 #if defined HAVE_BIGDFT 640 nfftot=product(wvl%den%denspot%dpbox%ndims) 641 #else 642 BIGDFT_NOTENABLED_ERROR() 643 #endif 644 end if 645 ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 646 ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 647 ABI_MALLOC(indsym,(4,dtset%nsym,dtset%natom)) 648 ABI_MALLOC(symrec,(3,3,dtset%nsym)) 649 irrzon(:,:,:)=0 650 phnons(:,:,:)=zero 651 indsym(:,:,:)=0 652 symrec(:,:,:)=0 653 654 if (dtset%nsym>1) then 655 call setsym(indsym,irrzon,dtset%iscf,dtset%natom,& 656 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 657 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred) 658 659 ! Make sure dtset%iatfix does not break symmetry 660 call fixsym(dtset%iatfix,indsym,dtset%natom,dtset%nsym) 661 else 662 ! The symrec array is used by initberry even in case nsym = 1 663 symrec(:,:,1) = 0 664 symrec(1,1,1) = 1 ; symrec(2,2,1) = 1 ; symrec(3,3,1) = 1 665 end if 666 if (dtset%usewvl == 1) then 667 call wvl_descr_atoms_set_sym(wvl%descr, dtset%efield, irrzon, dtset%nsppol, & 668 & dtset%nsym, phnons, dtset%symafm, dtset%symrel, dtset%tnons, dtset%tolsym) 669 #if defined HAVE_BIGDFT 670 wvl%den%symObj = wvl%descr%atoms%astruct%sym%symObj 671 #endif 672 end if 673 674 !########################################################### 675 !### 05. Calls inwffil 676 ABI_NVTX_START_RANGE(NVTX_INIT_INWFFIL) 677 678 ! if paral_kgb == 0, it may happen that some processors are idle (no entry in proc_distrb) 679 ! but mkmem == nkpt and this can cause integer overflow in mcg or allocation error. 680 ! Here we count the number of states treated by the proc. if cnt == 0, mcg is then set to 0. 681 cnt = 0 682 nbandtot = 0 683 do spin=1,dtset%nsppol 684 do ikpt=1,dtset%nkpt 685 nband_k = dtset%nband(ikpt + (spin-1) * dtset%nkpt) 686 do band=1,nband_k 687 if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, band, band, spin, mpi_enreg%me_kpt)) cnt = cnt + 1 688 end do 689 nbandtot = nbandtot + nband_k 690 end do 691 end do 692 693 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 694 mcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 695 if (cnt == 0) then 696 mcg = 0 697 write(msg,"(2(a,i0))")"rank: ",mpi_enreg%me, "does not have wavefunctions to treat. Setting mcg to: ",mcg 698 ABI_WARNING(msg) 699 end if 700 701 if (dtset%usewvl == 0 .and. dtset%mpw > 0 .and. cnt /= 0)then 702 if (my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol > floor(real(HUGE(0))/real(dtset%mpw) )) then 703 write (msg,'(9a)')& 704 "Default integer is not wide enough to store the size of the wavefunction array (mcg).",ch10,& 705 "This usually happens when paral_kgb == 0 and there are not enough procs to distribute kpts and spins",ch10,& 706 "Action: if paral_kgb == 0, use nprocs = nkpt * nsppol to reduce the memory per node.",ch10,& 707 "If this does not solve the problem, use paral_kgb 1 with nprocs > nkpt * nsppol and use npfft/npband/npspinor",ch10,& 708 "to decrease the memory requirements. Consider also OpenMP threads." 709 ii = 0 710 ABI_ERROR_NOSTOP(msg,ii) 711 write (msg,'(5(a,i0), 2a)')& 712 "my_nspinor: ",my_nspinor, ", mpw: ",dtset%mpw, ", mband: ",dtset%mband,& 713 ", mkmem: ",dtset%mkmem, ", nsppol: ",dtset%nsppol,ch10,& 714 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS...) compiled in int64 mode' 715 ABI_ERROR(msg) 716 end if 717 end if 718 719 if (dtset%imgwfstor==1) then 720 cg => scf_history%cg(:,:,1) 721 eigen => scf_history%eigen(:,1) 722 else 723 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 724 #if defined HAVE_GPU && defined HAVE_YAKL 725 ABI_MALLOC_MANAGED(cg, (/2,mcg/)) 726 #endif 727 else 728 ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr) 729 end if 730 ABI_MALLOC(eigen,(dtset%mband*dtset%nkpt*dtset%nsppol)) 731 end if 732 733 ABI_MALLOC(resid,(dtset%mband*dtset%nkpt*dtset%nsppol)) 734 eigen(:)=zero ; resid(:)=zero 735 !mpi_enreg%paralbd=0 ; ask_accurate=0 736 ask_accurate=0 737 738 !WVL - Branching, allocating wavefunctions as wavelets. 739 if (dtset%usewvl == 1) then 740 call wvl_wfs_lr_copy(wvl%wfs, wvl%descr) 741 ! Create access arrays for wavefunctions and allocate wvl%wfs%psi (other arrays are left unallocated). 742 call wvl_wfs_set(dtset%strprecon,dtset%spinmagntarget, dtset%kpt, mpi_enreg%me_wvl,& 743 & dtset%natom, sum(dtset%nband), & 744 & dtset%nkpt, mpi_enreg%nproc_wvl, dtset%nspinor, dtset%nsppol, dtset%nwfshist, occ, & 745 & psps, rprimd, wvl%wfs, dtset%wtk, wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult, & 746 & xred) 747 ! We transfer wavelets information to the hdr structure. 748 #if defined HAVE_BIGDFT 749 call local_potential_dimensions(mpi_enreg%me_wvl,wvl%wfs%ks%lzd,wvl%wfs%ks%orbs,wvl%den%denspot%xc,& 750 & wvl%den%denspot%dpbox%ngatherarr(0,1)) 751 hdr%nwvlarr(1) = wvl%wfs%ks%lzd%Glr%wfd%nvctr_c 752 hdr%nwvlarr(2) = 7 * wvl%wfs%ks%lzd%Glr%wfd%nvctr_f 753 #endif 754 ! Create access arrays for projectors and allocate them. 755 ! Compute projectors from each atom. 756 call wvl_projectors_set(mpi_enreg%me_wvl, dtset%natom, wvl%projectors, psps, rprimd, & 757 & wvl%wfs, wvl%descr, dtset%wvl_frmult, xred) 758 end if 759 760 read_wf_or_den=(dtset%iscf<=0.or.dtfil%ireadwf/=0.or.(dtfil%ireadden/=0.and.dtset%positron<=0)) 761 read_wf_or_den=(read_wf_or_den.and.has_to_init) 762 763 !RECURSION - initialization 764 if(has_to_init .and. dtset%userec==1) then 765 call InitRec(dtset,mpi_enreg,rec_set,rmet,maxval(psps%indlmn(3,:,:))) 766 end if 767 768 !LOTF - initialization 769 #if defined HAVE_LOTF 770 if(has_to_init .and. dtset%ionmov==23) then 771 call lotfparam_init(dtset%natom,dtset%lotf_version,1,& 772 & dtset%lotf_nitex,dtset%lotf_nneigx,& 773 & dtset%lotf_classic,1,1) 774 end if 775 #endif 776 777 !Initialize wavefunctions. 778 if(dtset%imgwfstor==1 .and. initialized==1)then 779 cg(:,:)=scf_history%cg(:,:,1) 780 eigen(:)=scf_history%eigen(:,1) 781 else if(dtset%tfkinfunc /=2) then 782 !if(dtset%tfkinfunc /=2) then 783 wff1%unwff=dtfil%unwff1 784 optorth=1 !if (psps%usepaw==1) optorth=0 785 if(psps%usepaw==1 .and. dtfil%ireadwf==1)optorth=0 786 hdr%rprimd=rprimd_for_kg ! We need the rprimd that was used to generate de G vectors 787 ABI_NVTX_START_RANGE(NVTX_INIT_INWFFIL2) 788 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen,& 789 & dtset%exchn2n3d,formeig,hdr,dtfil%ireadwf,dtset%istwfk,kg,& 790 & dtset%kptns,dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,& 791 & dtset%mpw,dtset%nband,ngfft,dtset%nkpt,npwarr,& 792 & dtset%nsppol,dtset%nsym,occ,optorth,dtset%symafm,& 793 & dtset%symrel,dtset%tnons,dtfil%unkg,wff1,wffnow,dtfil%unwff1,& 794 & dtfil%fnamewffk,wvl) 795 ABI_NVTX_END_RANGE() 796 hdr%rprimd=rprimd 797 end if 798 799 if (psps%usepaw==1.and.dtfil%ireadwf==1)then 800 call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 801 end if 802 803 ! Now that wavefunctions are initialized, calls of nonlocal operations are possible, so we start the counting (if enabled) 804 if (dtset%useylm==1.and.dtset%nonlop_ylm_count/=0.and.dtset%paral_kgb==0) then 805 call nonlop_ylm_init_counters() 806 end if 807 ! Same for fft counters 808 if (dtset%fft_count/=0.and.dtset%paral_kgb==0) then 809 call fft_init_counters() 810 end if 811 ABI_NVTX_END_RANGE() 812 813 !########################################################### 814 !### 06. Operations related to restartxf (Old version) 815 816 !TODO: Remove ab_xfh 817 !Initialize xf history (should be put in inwffil) 818 ab_xfh%nxfh=0 819 if(dtset%restartxf>=1 .and. dtfil%ireadwf==1)then 820 821 ! Should exchange the data about history in parallel localrdwf==0 822 if(xmpi_paral==1 .and. dtset%localrdwf==0)then 823 write(msg, '(a,a,a)' )& 824 & 'It is not yet possible to use non-zero restartxf,',ch10,& 825 & 'in parallel, when localrdwf=0. Sorry for this ...' 826 ABI_BUG(msg) 827 end if 828 829 ABI_MALLOC(ab_xfh%xfhist,(3,dtset%natom+4,2,0)) 830 call outxfhist(ab_xfh,dtset%natom,2,wff1,ios) 831 ABI_FREE(ab_xfh%xfhist) 832 833 if(ios>0)then 834 write(msg,'(a,a,a)')& 835 & 'An error occurred reading the input wavefunction file,',ch10,& 836 & 'with restartxf=1.' 837 ABI_ERROR(msg) 838 else if(ios==0)then 839 write(msg, '(a,a,i4,a)' )ch10,& 840 & ' gstate : reading',ab_xfh%nxfh,' (x,f) history pairs from input wf file.' 841 call wrtout([std_out, ab_out], msg) 842 end if 843 ! WARNING : should check that restartxf is not negative 844 ! WARNING : should check that restartxf /= only when dtfil%ireadwf is activated 845 end if 846 847 !Allocate the xf history array : takes into account the existing 848 !pairs, minus those that will be discarded, then those that will 849 !be computed, governed by dtset%ntime, and some additional pairs 850 !(needed when it will be possible to use xfhist for move.f) 851 ab_xfh%mxfh=(ab_xfh%nxfh-dtset%restartxf+1)+dtset%ntime+5 852 ABI_MALLOC(ab_xfh%xfhist,(3,dtset%natom+4,2,ab_xfh%mxfh)) 853 ab_xfh%xfhist(:,:,:,:) = zero 854 !WARNING : should check that the number of atoms in the wf file and natom are the same 855 856 !Initialize the xf history array 857 if(ab_xfh%nxfh>=dtset%restartxf .and. ab_xfh%nxfh>0)then 858 ! Eventually skip some of the previous history 859 if(dtset%restartxf>=2)then 860 do ixfh=1,dtset%restartxf-1 861 call WffReadSkipRec(ios,1,wff1) 862 end do 863 end if 864 865 ! Read and store the relevant history 866 ab_xfh%nxfhr=ab_xfh%nxfh-dtset%restartxf+1 867 call outxfhist(ab_xfh,dtset%natom,3,wff1,ios) 868 end if 869 870 !Close wff1, if it was ever opened (in inwffil) 871 if (dtfil%ireadwf==1) then 872 call WffClose(wff1,ierr) 873 end if 874 875 !########################################################### 876 !### 07. Calls setup2 877 878 !Further setup 879 ABI_MALLOC(start,(3,dtset%natom)) 880 call setup2(dtset,npwtot,start,wvl%wfs,xred) 881 882 !Allocation of previous atomic positions 883 if (scf_history%history_size>=0) then 884 xred_old => scf_history%xred_last 885 else 886 ABI_MALLOC(xred_old,(3,dtset%natom)) 887 end if 888 if (has_to_init) xred_old=xred 889 890 !Initialize (eventually) extfpmd object 891 if(dtset%useextfpmd>=1.and.dtset%occopt==3) then 892 if(dtset%useextfpmd/=1.and.dtset%extfpmd_nbcut>dtset%mband) then 893 write(msg,'(3a,i0,a,i0,3a)') "Not enough bands to activate ExtFPMD routines.",ch10,& 894 & "extfpmd_nbcut = ",dtset%extfpmd_nbcut," should be less than or equal to nband = ",dtset%mband,".",ch10,& 895 & "Action: Increase nband or decrease extfpmd_nbcut." 896 ABI_ERROR(msg) 897 else 898 if(dtset%useextfpmd/=1.and.(dtset%extfpmd_nbdbuf+dtset%extfpmd_nbcut)>dtset%mband) then 899 write(msg,'(a,i0,a,i0,a,i0,2a,i0,3a)') "(extfpmd_nbdbuf = ",dtset%extfpmd_nbdbuf," + extfpmd_nbcut = ",& 900 & dtset%extfpmd_nbcut,") = ",dtset%extfpmd_nbdbuf+dtset%extfpmd_nbcut,ch10,& 901 & "should be less than or equal to nband = ",dtset%mband,".",ch10,& 902 & "Assume experienced user. Execution will continue with extfpmd_nbdbuf = 0." 903 ABI_WARNING(msg) 904 dtset%extfpmd_nbdbuf = 0 905 else if(dtset%extfpmd_nbdbuf>dtset%mband) then 906 write(msg,'(a,i0,a,i0,3a)') "extfpmd_nbdbuf = ",dtset%extfpmd_nbdbuf,& 907 & " should be less than or equal to nband = ",dtset%mband,".",ch10,& 908 & "Assume experienced user. Execution will continue with extfpmd_nbdbuf = 0." 909 ABI_WARNING(msg) 910 dtset%extfpmd_nbdbuf = 0 911 end if 912 ABI_MALLOC(extfpmd,) 913 call extfpmd%init(dtset%mband,dtset%extfpmd_nbcut,dtset%extfpmd_nbdbuf,& 914 & dtset%nfft,dtset%nspden,rprimd,dtset%useextfpmd) 915 end if 916 end if 917 918 !Timing for initialisation period 919 call timab(1211,2,tsec) 920 call timab(1213,3,tsec) 921 922 923 !########################################################### 924 !### 08. Compute new occupation numbers 925 926 !Compute new occupation numbers, in case wavefunctions and eigenenergies 927 !were read from disk, occupation scheme is metallic (this excludes iscf=-1), 928 !and occupation numbers are required by iscf 929 if( dtfil%ireadwf==1 .and. & 930 & (dtset%occopt>=3.and.dtset%occopt<=9) .and. & 931 & (dtset%iscf>0 .or. dtset%iscf==-3) .and. dtset%positron/=1 ) then 932 933 ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 934 ! Warning : ideally, results_gs%entropy should not be set up here XG 20011007 935 ! Do not take into account the possible STM bias 936 call newocc(doccde,eigen,results_gs%energies%entropy,& 937 & results_gs%energies%e_fermie,results_gs%energies%e_fermih,dtset%ivalence,& 938 & dtset%spinmagntarget,dtset%mband,dtset%nband,& 939 & dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,& 940 & dtset%occopt,dtset%prtvol,dtset%tphysel,dtset%tsmear,dtset%wtk,& 941 & extfpmd=extfpmd) 942 if (dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(args_gs%upawu(:))>=tol8.or. & 943 & sum(args_gs%jpawu(:))>tol8).and.dtset%dmft_entropy==0) results_gs%energies%entropy=zero 944 ABI_FREE(doccde) 945 946 if(associated(extfpmd)) then 947 extfpmd%nelect=zero 948 call extfpmd%compute_nelect(results_gs%energies%e_fermie,extfpmd%nelect,& 949 & dtset%tsmear) 950 call extfpmd%compute_e_kinetic(results_gs%energies%e_fermie,dtset%tsmear) 951 end if 952 953 ! Transfer occupations to bigdft object: 954 if(dtset%usewvl==1 .and. .not. wvlbigdft) then 955 call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,1,wvl%wfs) 956 ! call wvl_energies_abi2big(results_gs%energies,wvl%wfs,2) 957 end if 958 959 else 960 ! Warning : ideally, results_gs%entropy should not be set up here XG 20011007 961 results_gs%energies%entropy=zero 962 end if 963 964 !########################################################### 965 !### 09. Generate an index table of atoms 966 967 !Definition of atindx array 968 !Generate an index table of atoms, in order for them to be used type after type. 969 ABI_MALLOC(atindx,(dtset%natom)) 970 ABI_MALLOC(atindx1,(dtset%natom)) 971 ABI_MALLOC(nattyp,(psps%ntypat)) 972 indx=1 973 do itypat=1,psps%ntypat 974 nattyp(itypat)=0 975 do iatom=1,dtset%natom 976 if(dtset%typat(iatom)==itypat)then 977 atindx(iatom)=indx 978 atindx1(indx)=iatom 979 indx=indx+1 980 nattyp(itypat)=nattyp(itypat)+1 981 end if 982 end do 983 end do 984 985 !Compute structure factor phases for current atomic pos: 986 if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then 987 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*dtset%natom)) 988 call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred) 989 end if 990 991 !Here allocation of GPU for vtorho calculations 992 #if defined HAVE_GPU 993 if (dtset%gpu_option/=ABI_GPU_DISABLED) then 994 call alloc_hamilt_gpu(atindx1,dtset,gprimd,mpi_enreg,nattyp,npwarr,2,psps,dtset%gpu_option) 995 end if 996 #endif 997 998 !########################################################### 999 !### 10. PAW related operations 1000 1001 !Initialize paw_dmft, even if neither dmft not paw are used 1002 !write(std_out,*) "dtset%usedmft",dtset%usedmft 1003 use_sc_dmft=dtset%usedmft 1004 ! if(dtset%paral_kgb>0) use_sc_dmft=0 1005 call init_sc_dmft(dtset%nbandkss,dtset%dmftbandi,dtset%dmftbandf,dtset%dmft_read_occnd,dtset%mband,& 1006 & dtset%nband,dtset%nkpt,dtset%nspden, & 1007 & dtset%nspinor,dtset%nsppol,occ,dtset%usedmft,paw_dmft,use_sc_dmft,dtset%dmft_solv,mpi_enreg) 1008 if (paw_dmft%use_dmft==1.and.me==0) then 1009 call readocc_dmft(paw_dmft,dtfil%filnam_ds(3),dtfil%filnam_ds(4)) 1010 end if 1011 !Should be done inside init_sc_dmft 1012 if ( dtset%usedmft /= 0 ) then 1013 call data4entropyDMFT_init(paw_dmft%forentropyDMFT,& 1014 dtset%natom,& 1015 dtset%typat,& 1016 dtset%lpawu,& 1017 dtset%dmft_t2g==1, & 1018 args_gs%upawu,& 1019 args_gs%jpawu) 1020 end if 1021 !write(std_out,*) "paw_dmft%use_dmft",paw_dmft%use_dmft 1022 1023 !PAW: 1- Initialize values for several arrays unchanged during iterations 1024 !2- Initialize data for DFT+U 1025 !3- Eventually open temporary storage file 1026 if(psps%usepaw==1) then 1027 ! 1- 1028 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 1029 1030 ! Test if we have to call pawinit 1031 ! Some gen-cond have to be added... 1032 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 1033 1034 if (psp_gencond==1.or.call_pawinit) then 1035 call timab(553,1,tsec) 1036 gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 1037 hyb_range_fock=zero;if (dtset%ixc<0) call libxc_functionals_get_hybridparams(hyb_range=hyb_range_fock) 1038 call pawinit(dtset%effmass_free,gnt_option,gsqcut_shp,hyb_range_fock,dtset%pawlcutd,dtset%pawlmix,& 1039 & psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,& 1040 & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero) 1041 1042 ! Update internal values 1043 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 1044 call timab(553,2,tsec) 1045 #if defined HAVE_BIGDFT 1046 ! In the PAW+WVL case, copy sij: 1047 if(dtset%usewvl==1) then 1048 do itypat=1,dtset%ntypat 1049 wvl%descr%paw%sij(1:pawtab(itypat)%lmn2_size,itypat)=pawtab(itypat)%sij(:) 1050 end do 1051 end if 1052 #endif 1053 end if 1054 psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore) 1055 call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot) 1056 1057 ! 2-Initialize and compute data for DFT+U, EXX, or DFT+DMFT 1058 if(paw_dmft%use_dmft==1) call print_sc_dmft(paw_dmft,dtset%pawprtvol) 1059 call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,& 1060 & is_dfpt,args_gs%jpawu,dtset%lexexch,dtset%lpawu,dtset%nspinor,dtset%ntypat,dtset%optdcmagpawu,pawang,dtset%pawprtvol,& 1061 & pawrad,pawtab,args_gs%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu,ucrpa=dtset%ucrpa) 1062 1063 ! DEBUG: 1064 !if (me == master) call pawtab_print(Pawtab) 1065 end if 1066 1067 !########################################################### 1068 !### 11. Initialize (eventually) electron-positron data and 1069 !### electric and magnetic field data 1070 1071 !Initialize (eventually) electron-positron data 1072 nullify (electronpositron) 1073 if (dtset%positron/=0) then 1074 call init_electronpositron(dtfil%ireadwf,dtset,electronpositron,mpi_enreg,nfftf,pawrhoij,pawtab) 1075 end if 1076 1077 !########################################################### 1078 ! Initialisation of cprj 1079 usecprj=0; mcprj=0;mband_cprj=0 1080 compute_cprj=.false. 1081 ! PAW keeping cprj in memory : some cases are excluded for now... 1082 if (dtset%cprj_in_memory/=0) then 1083 compute_cprj=.true. 1084 usecprj=1 1085 else 1086 if (dtset%usepaw==1) then 1087 if (associated(electronpositron)) then 1088 if (dtset%positron/=0.and.electronpositron%dimcprj>0) usecprj=1 1089 end if 1090 if (dtset%prtnabla>0) usecprj=1 1091 if (dtset%extrapwf>0) usecprj=1 1092 if (dtset%pawfatbnd>0)usecprj=1 1093 if (dtset%prtdos==3) usecprj=1 1094 if (dtset%usewvl==1) usecprj=1 1095 if (dtset%nstep==0) usecprj=0 1096 if (dtset%usefock==1) usecprj=1 1097 end if 1098 end if 1099 if (usecprj==0) then 1100 ABI_MALLOC(cprj,(0,0)) 1101 end if 1102 if (usecprj==1) then 1103 mband_cprj=dtset%mband;if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band 1104 mcprj=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol 1105 !Was allocated above for valgrind sake so should always be true (safety) 1106 if (allocated(cprj)) then 1107 call pawcprj_free(cprj) 1108 ABI_FREE(cprj) 1109 end if 1110 ABI_MALLOC(cprj,(dtset%natom,mcprj)) 1111 ncpgr=0 1112 if (dtset%usefock==1) then ! Note that compute_cprj = false if usefock/=0 1113 if (dtset%optforces == 1) then 1114 ncpgr = 3 1115 end if 1116 ! if (dtset%optstress /= 0) then 1117 ! ncpgr = 6 ; ctocprj_choice = 3 1118 ! end if 1119 end if 1120 ABI_MALLOC(dimcprj_srt,(dtset%natom)) 1121 call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 1122 call pawcprj_alloc(cprj,ncpgr,dimcprj_srt) 1123 if (compute_cprj) then 1124 choice = 1 ! no derivative... 1125 idir = 0 ! ...so no direction 1126 iatom = 0 ! all atoms 1127 iorder_cprj = 0 ! ordered by atom types 1128 ncprj = dtset%natom 1129 ! Compute structure factor phases and large sphere cut-off (gsqcut): 1130 ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom)) 1131 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 1132 call wrtout(std_out,' Computing cprj from initial wavefunctions (gstate)') 1133 call ctocprj(atindx,cg,choice,cprj,gmet,gprimd,iatom,idir,& 1134 & iorder_cprj,dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,& 1135 & dtset%mpw,dtset%natom,nattyp,dtset%nband,ncprj,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,& 1136 & dtset%nsppol,dtset%nsppol,psps%ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,ucvol,dtfil%unpaw,& 1137 & xred,ylm,ylmgr) 1138 call wrtout(std_out,' cprj is computed') 1139 ABI_FREE(ph1d) 1140 end if 1141 end if 1142 1143 !Timing for initialisation period 1144 call timab(1213,2,tsec) 1145 call timab(1214,3,tsec) 1146 1147 !########################################################### 1148 !### 12. Operations dependent of iscf value 1149 1150 !Get starting charge density : rhor as well as rhog 1151 !Also initialize the kinetic energy density 1152 if (scf_history%history_size>=0) then 1153 rhor => scf_history%rhor_last 1154 taur => scf_history%taur_last 1155 else 1156 ABI_MALLOC(rhor,(nfftf,dtset%nspden)) 1157 ABI_MALLOC(taur,(nfftf,dtset%nspden*dtset%usekden)) 1158 end if 1159 ABI_MALLOC(rhog,(2,nfftf)) 1160 ABI_MALLOC(taug,(2,nfftf*dtset%usekden)) 1161 1162 if (has_to_init) then 1163 1164 ! === Self-consistent case 1165 if (dtset%iscf>0 .or. (dtset%iscf==0 .and. dtset%usewvl==1 )) then 1166 1167 ! >>> Initialize charge density 1168 1169 if (dtfil%ireadden/=0.and.dtset%positron<=0) then 1170 ! Choice 1: read charge density from file 1171 rdwrpaw=psps%usepaw ; if(dtfil%ireadwf/=0) rdwrpaw=0 1172 if (dtset%usewvl==0) then 1173 call read_rhor(dtfil%fildensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw, & 1174 mpi_enreg,rhor,hdr_den,pawrhoij,comm,check_hdr=hdr,allow_interp=.True.) 1175 results_gs%etotal = hdr_den%etot 1176 call hdr_den%free() 1177 else 1178 fform=52 ; accessfil=0 1179 if (dtset%iomode == IO_MODE_MPI ) accessfil=4 1180 if (dtset%iomode == IO_MODE_ETSF) accessfil=3 1181 call ioarr(accessfil,rhor,dtset,results_gs%etotal,fform,dtfil%fildensin,hdr,& 1182 mpi_enreg,ngfftf,cplex1,nfftf,pawrhoij,1,rdwrpaw,wvl%den) 1183 end if 1184 if (rdwrpaw/=0) then 1185 call hdr%update(bantot,etot,fermie,fermih,residm,& 1186 rprimd,occ,pawrhoij,xred,args_gs%amu,& 1187 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1188 end if 1189 ! Compute up+down rho(G) by fft 1190 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 1191 1192 else if (dtfil%ireadwf/=0) then 1193 ! Choice 2: obtain charge density from wfs that were read previously 1194 ! Warning: in PAW, rho does not include the compensation density (added later) 1195 tim_mkrho=1 1196 if (psps%usepaw==1) then 1197 ABI_MALLOC(rhowfg,(2,dtset%nfft)) 1198 ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden)) 1199 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 1200 & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,extfpmd=extfpmd) 1201 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor) 1202 ABI_FREE(rhowfg) 1203 ABI_FREE(rhowfr) 1204 else 1205 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 1206 & rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,extfpmd=extfpmd) 1207 end if 1208 1209 else if (dtfil%ireadwf==0.and.dtset%positron/=1) then 1210 ! Choice 3: crude, but realistic initialisation of the charge density 1211 ! There is not point to compute it from random wavefunctions 1212 if (dtset%usewvl == 0) then 1213 call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,& 1214 & mgfftf,mpi_enreg,psps%mqgrid_vl,dtset%natom,nattyp,nfftf,& 1215 & ngfftf,dtset%nspden,psps%ntypat,psps,pawtab,ph1df,& 1216 & psps%qgrid_vl,rhog,rhor,dtset%spinat,ucvol,psps%usepaw,& 1217 & dtset%ziontypat,dtset%znucl) 1218 else 1219 if (dtset%usepaw==0) then 1220 !Wavelet density corresponds exactly to the wavefunctions, 1221 !since wavefunctions are taken from diagonalisation of LCAO. 1222 call wvl_mkrho(dtset,irrzon,mpi_enreg,phnons,rhor,wvl%wfs,wvl%den) 1223 else 1224 #if defined HAVE_BIGDFT 1225 call wvl_initro(atindx1,wvl%descr%atoms%astruct%geocode,wvl%descr%h,& 1226 & mpi_enreg%me_wvl,dtset%natom,nattyp,nfftf,dtset%nspden,psps%ntypat,& 1227 & wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,& 1228 & wvl%descr%Glr%d%n3,pawrad,pawtab,psps%gth_params%psppar,rhor,rprimd,& 1229 & dtset%spinat,wvl%den,dtset%xc_denpos,xred,dtset%ziontypat) 1230 call wvl_mkrho(dtset,irrzon,mpi_enreg,phnons,rhor,wvl%wfs,wvl%den) 1231 #endif 1232 end if ! usepaw 1233 end if ! usewvl 1234 ! Update initialized density taking into account jellium slab 1235 if (dtset%jellslab/=0) then 1236 option=2 1237 ABI_MALLOC(work,(nfftf)) 1238 call jellium(gmet,gsqcut_eff,mpi_enreg,nfftf,ngfftf,dtset%nspden,option,& 1239 & dtset%slabwsrad,rhog,rhor,rprimd,work,dtset%slabzbeg,dtset%slabzend) 1240 ABI_FREE(work) 1241 end if ! of usejell 1242 1243 end if ! choice for charge density initialization 1244 1245 ! >>> Initialize kinetic energy density 1246 if (dtset%usekden==1) then 1247 1248 if (dtfil%ireadkden/=0.and.dtset%positron<=0) then 1249 ! Choice 1: read kinetic energy density from file 1250 rdwrpaw=0 1251 call read_rhor(dtfil%filkdensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw, & 1252 mpi_enreg,taur,hdr_den,pawrhoij,comm,check_hdr=hdr,allow_interp=.True.) 1253 call hdr_den%free() 1254 ! Compute up+down tau(G) by fft 1255 call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 1256 1257 else if (dtfil%ireadwf/=0) then 1258 ! Choice 2: obtain kinetic energy density from wfs that were read previously 1259 tim_mkrho=1 1260 if (psps%usepaw==1) then 1261 ABI_MALLOC(rhowfg,(2,dtset%nfft)) 1262 ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden)) 1263 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 1264 & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 1265 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,taug,rhowfr,taur) 1266 ABI_FREE(rhowfg) 1267 ABI_FREE(rhowfr) 1268 else 1269 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 1270 & taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 1271 end if 1272 1273 else if(dtfil%ireadwf==0.and.dtset%positron/=1)then 1274 ! Choice 3: kinetic energy density initialized to zero (?) 1275 taur=zero ; taug=zero 1276 1277 end if ! choice for kinetic energy density initialization 1278 end if ! usekden 1279 1280 ! === Non self-consistent case 1281 else if ((dtset%iscf==-1.or.dtset%iscf==-2.or.dtset%iscf==-3).and.dtset%positron<=0) then 1282 1283 ! Read density from a disk file (this is mandatory for non-self-consistent calculations) 1284 ! Note : results_gs%etotal is read here, 1285 ! and might serve in the tddft routine, but it is contrary to the intended use of results_gs ... 1286 ! Warning : should check the use of results_gs%e_fermie 1287 ! Warning : should check the use of results_gs%residm 1288 ! One might make them separate variables. 1289 1290 ! Read charge density and get Fermi level from hdr_den 1291 rdwrpaw=psps%usepaw 1292 call read_rhor(dtfil%fildensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw,& 1293 & mpi_enreg,rhor,hdr_den,pawrhoij,comm,check_hdr=hdr) 1294 results_gs%etotal = hdr_den%etot; 1295 results_gs%energies%e_fermie = hdr_den%fermie 1296 results_gs%energies%e_fermih = hdr_den%fermih 1297 ! Compute up+down rho(G) by fft 1298 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 1299 call hdr_den%free() 1300 1301 ! Read kinetic energy density 1302 if(dtset%usekden==1)then 1303 rdwrpaw=0 1304 call read_rhor(dtfil%filkdensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw, & 1305 mpi_enreg,taur,hdr_den,pawrhoij,comm,check_hdr=hdr) 1306 call hdr_den%free() 1307 ! Compute up+down tau(G) by fft 1308 call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 1309 end if 1310 1311 end if ! self-consistent/non self-consistent 1312 end if ! has_to_init 1313 1314 !Timing for initialisation period 1315 call timab(1214,2,tsec) 1316 call timab(1215,3,tsec) 1317 1318 !########################################################### 1319 !### 13. If needed, initialize SCF history variables 1320 1321 !If needed, initialize atomic density in SCF history 1322 if (scf_history%history_size>0.and.has_to_init) then 1323 ! If rhor is an atomic density, just store it in history 1324 if (.not.read_wf_or_den) then 1325 scf_history%atmrho_last(:)=rhor(:,1) 1326 else 1327 ! If rhor is not an atomic density, has to compute rho_at(r) 1328 ABI_MALLOC(rhowfg,(2,nfftf)) 1329 ABI_MALLOC(rhowfr,(nfftf,1)) 1330 ABI_MALLOC(spinat_dum,(3,dtset%natom)) 1331 spinat_dum=zero 1332 call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,mgfftf,mpi_enreg,& 1333 & psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,1,psps%ntypat,psps,pawtab,& 1334 & ph1df,psps%qgrid_vl,rhowfg,rhowfr,spinat_dum,ucvol,& 1335 & psps%usepaw,dtset%ziontypat,dtset%znucl) 1336 scf_history%atmrho_last(:)=rhowfr(:,1) 1337 ABI_FREE(rhowfg) 1338 ABI_FREE(rhowfr) 1339 ABI_FREE(spinat_dum) 1340 end if 1341 end if 1342 1343 if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then 1344 ABI_FREE(ph1df) 1345 end if 1346 1347 !!Electric field: initialization stage 1348 !!further initialization and updates happen in scfcv.F90 1349 call init_e_field_vars(dtefield,dtset,gmet,gprimd,kg,& 1350 & mpi_enreg,npwarr,occ,pawang,pawrad,pawtab,psps,& 1351 & pwind,pwind_alloc,pwnsfac,rprimd,symrec,xred) 1352 1353 fatvshift=one 1354 1355 !Check whether exiting was required by the user. If found then do not start minimization steps 1356 !At this first call to chkexi, initialize cpus, if it 1357 !is non-zero (which would mean that no action has to be taken) 1358 !Should do this in driver ... 1359 cpus=dtset%cpus 1360 if(abs(cpus)>1.0d-5)cpus=cpus+cpui 1361 openexit=1 ; if(dtset%chkexit==0) openexit=0 1362 call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit) 1363 1364 !If immediate exit, and wavefunctions were not read, must zero eigenvalues 1365 if (iexit/=0) eigen(:)=zero 1366 1367 #if defined HAVE_BIGDFT 1368 if (dtset%usewvl == 1 .and. dtset%timopt==10) then 1369 call wvl_timing(xmpi_world,'== INITS','PR') 1370 end if 1371 #endif 1372 1373 call timab(1215,2,tsec) 1374 1375 conv_retcode = 0 1376 1377 if (iexit==0) then 1378 1379 ! ########################################################### 1380 ! ### 14. Move atoms and acell according to ionmov value 1381 1382 call timab(1225,3,tsec) 1383 1384 call scfcv_init(scfcv_args,atindx,atindx1,cg,cprj,cpus,& 1385 & args_gs%dmatpawu,dtefield,dtfil,dtpawuj,dtset,ecore,eigen,hdr,extfpmd,& 1386 & indsym,initialized,irrzon,kg,mcg,mcprj,mpi_enreg,my_natom,nattyp,ndtpawuj,& 1387 & nfftf,npwarr,occ,pawang,pawfgr,pawrad,pawrhoij,& 1388 & pawtab,phnons,psps,pwind,pwind_alloc,pwnsfac,rec_set,& 1389 & resid,results_gs,scf_history,fatvshift,& 1390 & symrec,taug,taur,wvl,ylm,ylmgr,paw_dmft,wffnew,wffnow) 1391 1392 call dtfil_init_time(dtfil,0) 1393 1394 write(msg,'(a,80a)')ch10,('=',mu=1,80) 1395 call wrtout([std_out, ab_out], msg) 1396 1397 if (dtset%ionmov==0 .or. dtset%imgmov==6) then 1398 1399 ! Should merge this call with the call for dtset%ionmov==4 and 5 1400 if (dtset%macro_uj==0) then 1401 itimes(1)=itime0 ; itimes(2)=itimimage_gstate 1402 call scfcv_run(scfcv_args,electronpositron,itimes,rhog,rhor,rprimd,xred,xred_old,conv_retcode) 1403 else 1404 ! Conduct determination of U 1405 call pawuj_drive(scfcv_args,dtset,electronpositron,rhog,rhor,rprimd,xred,xred_old) 1406 end if 1407 1408 ! ======================================== 1409 ! New structure for geometry optimization 1410 ! ======================================== 1411 else if (dtset%ionmov>50.or.dtset%ionmov<=28) then 1412 1413 ! TODO: return conv_retcode 1414 call mover(scfcv_args,ab_xfh,acell,args_gs%amu,dtfil,& 1415 & electronpositron,rhog,rhor,rprimd,vel,vel_cell,xred,xred_old,itimimage_gstate=itimimage_gstate) 1416 1417 ! Compute rprim from rprimd and acell 1418 do kk=1,3 1419 do jj=1,3 1420 rprim(jj,kk)=rprimd(jj,kk)/acell(kk) 1421 end do 1422 end do 1423 1424 ! ========================================= 1425 ! New structure for geometry optimization 1426 ! ========================================= 1427 1428 else ! Not an allowed option 1429 write(msg, '(a,i0,2a)' )& 1430 'Disallowed value for ionmov=',dtset%ionmov,ch10,& 1431 'Allowed values are: 1,2,3,4,5,6,7,8,9,10,11,12,13,14,20,21,22,23,24,28 and 30' 1432 ABI_BUG(msg) 1433 end if 1434 1435 call scfcv_destroy(scfcv_args) 1436 1437 call timab(1225,2,tsec) 1438 1439 ! ########################################################### 1440 ! ### 15. Final operations and output for gstate 1441 1442 end if ! End of the check of hasty exit 1443 1444 call timab(1226,3,tsec) 1445 1446 write(msg, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,& 1447 & ' ----iterations are completed or convergence reached----',ch10 1448 call wrtout([std_out, ab_out], msg) 1449 1450 !Mark this GS computation as done 1451 initialized=1 1452 1453 !Update the header, before using it 1454 call hdr%update(bantot,results_gs%etotal,results_gs%energies%e_fermie,results_gs%energies%e_fermih,& 1455 results_gs%residm,rprimd,occ,pawrhoij,xred,args_gs%amu,& 1456 comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1457 1458 ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1459 doccde=zero 1460 1461 call ebands_init(bantot,ebands,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,& 1462 & doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,& 1463 & hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,& 1464 & hdr%cellcharge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, & 1465 & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk) 1466 1467 ebands%fermie = results_gs%energies%e_fermie 1468 ebands%fermih = results_gs%energies%e_fermih 1469 ABI_FREE(doccde) 1470 1471 ! Compute and print the gaps. 1472 call ebands_report_gap(ebands,header="Gap info",unit=std_out,mode_paral="COLL",gaps=results_gs%gaps) 1473 1474 call timab(1226,2,tsec) 1475 call timab(1227,3,tsec) 1476 1477 if(dtset%nqpt==0)filnam=dtfil%fnameabo_wfk 1478 if(dtset%nqpt==1)filnam=dtfil%fnameabo_wfq 1479 1480 ! Write wavefunctions file only if convergence was not achieved. 1481 !write(std_out,*)"conv_retcode", conv_retcode 1482 write_wfk = .True. 1483 if (dtset%prtwf==-1 .and. conv_retcode == 0) then 1484 write_wfk = .False. 1485 msg = "GS calculation converged with prtwf=-1 --> Skipping WFK file output" 1486 call wrtout(ab_out, msg) 1487 ABI_COMMENT(msg) 1488 end if 1489 1490 !To print out the WFs, need the rprimd that was used to generate the G vectors 1491 hdr%rprimd=rprimd_for_kg 1492 1493 if (write_wfk) then 1494 call outresid(dtset,dtset%kptns,dtset%mband,dtset%nband,dtset%nkpt, dtset%nsppol,resid) 1495 1496 call outwf(cg,dtset,psps,eigen,filnam,hdr,kg,dtset%kptns,& 1497 dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%natom,& 1498 dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,& 1499 occ,response,dtfil%unwff2,wvl%wfs,wvl%descr) 1500 1501 ! Generate WFK with k-mesh from WFK containing list of k-points inside pockets. 1502 if (dtset%getkerange_filepath /= ABI_NOFILE) then 1503 call wfk_klist2mesh(dtfil%fnameabo_wfk, dtset%getkerange_filepath, dtset, comm) 1504 end if 1505 1506 !SPr: add input variable managing the .vtk file OUTPUT (Please don't remove the next commented line) 1507 !call printmagvtk(mpi_enreg,cplex1,dtset%nspden,nfftf,ngfftf,rhor,rprimd,'DEN') 1508 end if 1509 1510 if (dtset%prtwf==2) call outqmc(cg,dtset,eigen,gprimd,hdr,kg,mcg,mpi_enreg,npwarr,occ,psps,results_gs) 1511 1512 !Restore the original rprimd in hdr 1513 hdr%rprimd=rprimd 1514 1515 ! Generate WFK in full BZ (needed by LOBSTER) 1516 if (me == master .and. dtset%prtwf == 1 .and. dtset%prtwf_full == 1 .and. dtset%nqpt == 0) then 1517 wfkfull_path = strcat(dtfil%filnam_ds(4), "_FULL_WFK") 1518 if (dtset%iomode == IO_MODE_ETSF) wfkfull_path = nctk_ncify(wfkfull_path) 1519 call wfk_tofullbz(filnam, dtset, psps, pawtab, wfkfull_path, hdr_kfull) 1520 call hdr_kfull%free() 1521 call cryst%free() 1522 end if 1523 1524 call timab(1227,2,tsec) 1525 call timab(1228,3,tsec) 1526 1527 call clnup1(acell,dtset,eigen,results_gs%energies%e_fermie,results_gs%energies%e_fermih,& 1528 & dtfil%fnameabo_dos,dtfil%fnameabo_eig,results_gs%gred,& 1529 & mpi_enreg,nfftf,ngfftf,occ,dtset%optforces,& 1530 & resid,rhor,rprimd,results_gs%vxcavg,xred) 1531 1532 if ( (dtset%iscf>=0 .or. dtset%iscf==-3) .and. dtset%prtstm==0) then 1533 call prtene(dtset,results_gs%energies,ab_out,psps%usepaw) 1534 end if 1535 1536 call timab(1228,2,tsec) 1537 call timab(1229,3,tsec) 1538 1539 !write final electric field components HONG 1540 1541 if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt ==7 .or. & 1542 & dtset%berryopt == 14 .or. dtset%berryopt == 16 .or. dtset%berryopt ==17 ) then ! output final electric field data !!HONG 1543 if (dtset%berryopt == 4) then 1544 write(msg,'(a,a)') ch10, 'Constant unreduced E calculation - final values:' 1545 else if (dtset%berryopt == 6 ) then 1546 write(msg,'(a,a)') ch10, 'Constant unreduced D calculation - final values:' 1547 else if (dtset%berryopt == 14) then 1548 write(msg,'(a,a)') ch10, 'Constant reduced ebar calculation - final values:' 1549 else if (dtset%berryopt == 16 ) then 1550 write(msg,'(a,a)') ch10, 'Constant reduced d calculation - final values:' 1551 else if (dtset%berryopt == 17) then 1552 write(msg,'(a,a)') ch10, 'Constant reduced ebar and d calculation - final values:' 1553 end if 1554 1555 call wrtout([std_out, ab_out], msg) 1556 call prtefield(dtset,dtefield,ab_out,rprimd) 1557 call prtefield(dtset,dtefield,std_out,rprimd) 1558 1559 ! To check if the final electric field is below the critical field 1560 do kk = 1, 3 1561 efield_band(kk) = abs(dtset%red_efieldbar(kk))*dtefield%nkstr(kk) 1562 end do 1563 ! eg = maxval(eg_dir) 1564 ! eg_ev = eg*Ha_eV 1565 write(msg,'(a,a,a,a,a,a,a,a,f7.2,a,a)')ch10,& 1566 & ' Please check: COMMENT - ',ch10,& 1567 & ' As a rough estimate,',ch10,& 1568 & ' to be below the critical field, the bandgap of your system',ch10,& 1569 & ' should be larger than ',maxval(efield_band)*Ha_eV,' eV.',ch10 1570 call wrtout([std_out, ab_out], msg) 1571 1572 write(msg,'(a)') '--------------------------------------------------------------------------------' 1573 call wrtout([std_out, ab_out], msg) 1574 end if 1575 1576 call timab(1229,2,tsec) 1577 call timab(1230,3,tsec) 1578 1579 !In the // case, only master writes the energy and the gradients to the DDB 1580 if (me==0.and.dtset%nimage==1.and.((dtset%iscf > 0).or.& 1581 & (dtset%berryopt == -1).or.(dtset%berryopt) == -3)) then 1582 1583 ! DDB dimensions 1584 if (dtset%iscf > 0) then 1585 nblok = 2 ! 1st blok = gradients, 2nd blok = energy 1586 else 1587 nblok = 1 ! 1st blok = gradients 1588 end if 1589 mpert = dtset%natom + 6 1590 1591 ! Create header and ddb objects 1592 dscrpt=' Note : temporary (transfer) database ' 1593 call ddb_hdr%init(dtset,psps,pawtab,dscrpt,nblok,& 1594 & xred=xred,occ=occ,ngfft=ngfft) 1595 1596 call ddb%init(dtset, nblok, mpert, with_d1E=.true.) 1597 1598 ! Set the electronic polarization 1599 if ((abs(dtset%berryopt) == 1).or.(abs(dtset%berryopt) == 3)) then 1600 call ddb%set_pel(results_gs%pel, dtset%rfdir, 1) 1601 end if 1602 1603 if (dtset%iscf > 0) then 1604 1605 ! Set the gradients (forces) in reduced coordinates 1606 call ddb%set_gred(results_gs%gred, 1) 1607 1608 ! Set the stress tensor 1609 call ddb%set_strten(results_gs%strten, 1) 1610 1611 ! Set the total energy 1612 call ddb%set_etotal(results_gs%etotal, 2) 1613 end if 1614 1615 if (dtset%prtddb==1) then 1616 ! Write the DDB 1617 call ddb%write(ddb_hdr, dtfil%fnameabo_ddb, with_psps=0) 1618 end if 1619 1620 ! Free memory 1621 call ddb_hdr%free() 1622 call ddb%free() 1623 1624 end if 1625 1626 call timab(1230,2,tsec) 1627 call timab(1231,3,tsec) 1628 1629 1630 if (dtset%nstep>0 .and. dtset%prtstm==0 .and. dtset%positron/=1) then 1631 call clnup2(psps%n1xccc,results_gs%gred,results_gs%grchempottn,results_gs%gresid,& 1632 & results_gs%grewtn,results_gs%grvdw,results_gs%grxc,dtset%iscf,dtset%natom,& 1633 & results_gs%ngrvdw,dtset%optforces,dtset%optstress,dtset%prtvol,start,& 1634 & results_gs%strten,results_gs%synlgr,xred) 1635 end if 1636 1637 ! Write nonlop_ylm_counters (if enabled) in outputs 1638 if (dtset%useylm==1.and.dtset%nonlop_ylm_count/=0.and.dtset%paral_kgb==0) then 1639 call nonlop_ylm_output_counters(dtset%natom,nbandtot,dtset%ntypat,dtset%typat,mpi_enreg) 1640 end if 1641 ! Write fft_counters (if enabled) in output 1642 if (dtset%fft_count/=0.and.dtset%paral_kgb==0) then 1643 call fft_output_counters(nbandtot,mpi_enreg) 1644 end if 1645 1646 if(dtset%imgwfstor==1)then 1647 scf_history%cg(:,:,1)=cg(:,:) 1648 scf_history%eigen(:,1)=eigen(:) 1649 endif 1650 1651 !Deallocate arrays 1652 ABI_FREE(atindx) 1653 ABI_FREE(atindx1) 1654 ABI_FREE(indsym) 1655 ABI_FREE(npwarr) 1656 ABI_FREE(nattyp) 1657 ABI_FREE(resid) 1658 ABI_FREE(rhog) 1659 ABI_FREE(start) 1660 ABI_FREE(symrec) 1661 ABI_FREE(taug) 1662 ABI_FREE(ab_xfh%xfhist) 1663 call pawfgr_destroy(pawfgr) 1664 1665 if(dtset%imgwfstor==0)then 1666 if(dtset%gpu_option == ABI_GPU_KOKKOS) then 1667 #if defined HAVE_GPU && defined HAVE_YAKL 1668 ABI_FREE_MANAGED(cg) 1669 #endif 1670 else 1671 ABI_FREE(cg) 1672 end if 1673 ABI_FREE(eigen) 1674 else 1675 nullify(cg,eigen) 1676 endif 1677 1678 if (dtset%usewvl == 0) then 1679 ! In wavelet case, irrzon and phnons are deallocated by wavelet object. 1680 ABI_FREE(irrzon) 1681 ABI_FREE(phnons) 1682 end if 1683 1684 ABI_FREE(ylm) 1685 ABI_FREE(ylmgr) 1686 1687 if (scf_history%history_size<0) then 1688 if (psps%usepaw==1) then 1689 call pawrhoij_free(pawrhoij) 1690 end if 1691 ABI_FREE(rhor) 1692 ABI_FREE(taur) 1693 ABI_FREE(pawrhoij) 1694 ABI_FREE(xred_old) 1695 else 1696 nullify(rhor,taur,pawrhoij,xred_old) 1697 end if 1698 1699 !PAW+DMFT 1700 call destroy_sc_dmft(paw_dmft) 1701 ! This call should be done inside destroy_sc_dmft 1702 if ( dtset%usedmft /= 0 ) then 1703 call data4entropyDMFT_destroy(paw_dmft%forentropyDMFT) 1704 end if 1705 1706 !Destroy extfpmd datastructure 1707 if(associated(extfpmd)) then 1708 call extfpmd%destroy() 1709 ABI_FREE(extfpmd) 1710 end if 1711 1712 !Destroy electronpositron datastructure 1713 if (dtset%positron/=0) then 1714 call destroy_electronpositron(electronpositron) 1715 end if 1716 1717 !Deallocating the basis set. 1718 if (dtset%usewvl == 1) then 1719 call wvl_projectors_free(wvl%projectors) 1720 call wvl_wfs_free(wvl%wfs) 1721 call wvl_descr_free(wvl%descr) 1722 call wvl_denspot_free(wvl%den) 1723 if(dtset%usepaw == 1) then 1724 call wvl_paw_free(wvl%descr) 1725 end if 1726 end if 1727 1728 ABI_FREE(kg) 1729 1730 if (dtset%icoulomb /= 0) then 1731 call psolver_kernel((/ 0._dp, 0._dp, 0._dp /), 0, dtset%icoulomb, 0, kernel_dummy, & 1732 & 0, dtset%ngfft, 1, dtset%nscforder) 1733 end if 1734 1735 if (associated(pwind)) then 1736 ABI_FREE(pwind) 1737 end if 1738 if (associated(pwnsfac)) then 1739 ABI_FREE(pwnsfac) 1740 end if 1741 if ((dtset%berryopt<0).or.& 1742 & (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1743 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17)) then 1744 if (xmpi_paral == 1) then 1745 ABI_FREE(mpi_enreg%kptdstrb) 1746 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1747 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 1748 ABI_FREE(mpi_enreg%kpt_loc2ibz_sp) 1749 end if 1750 end if 1751 if (allocated(mpi_enreg%kpt_loc2ibz_sp)) then 1752 ABI_FREE(mpi_enreg%kpt_loc2ibz_sp) 1753 end if 1754 if (allocated(mpi_enreg%kpt_loc2fbz_sp)) then 1755 ABI_FREE(mpi_enreg%kpt_loc2fbz_sp) 1756 end if 1757 if (allocated(mpi_enreg%mkmem)) then 1758 ABI_FREE(mpi_enreg%mkmem) 1759 end if 1760 end if 1761 ! deallocate cprj 1762 if(usecprj==1) then 1763 ABI_FREE(dimcprj_srt) 1764 call pawcprj_free(cprj) 1765 end if 1766 ABI_FREE(cprj) 1767 1768 ! deallocate efield 1769 call destroy_efield(dtefield) 1770 1771 !deallocate Recursion 1772 if (dtset%userec == 1) then 1773 call CleanRec(rec_set) 1774 end if 1775 1776 call hdr%free() 1777 call ebands_free(ebands) 1778 1779 if (me == master .and. dtset%prtxml == 1) then 1780 ! The dataset given in argument has been treated, then we output its variables. 1781 ! call outvarsXML() 1782 ! gstate() will handle a dataset, so we output the dataSet markup. 1783 write(ab_xml_out, "(A)") ' </dataSet>' 1784 end if 1785 1786 if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2 .and. dtset%optdriver /= RUNL_GWLS) then 1787 ! Plane-wave case 1788 call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg) 1789 end if 1790 1791 if((dtset%wfoptalg == 1 .or. dtset%wfoptalg == 111) .and. psps%usepaw == 1) then 1792 call destroy_invovl(dtset%nkpt,dtset%gpu_option) 1793 end if 1794 1795 !Clean gemm_nonlop work spaces 1796 if(gemm_nonlop_use_gemm) then 1797 if(dtset%gpu_option == ABI_GPU_OPENMP) then 1798 call destroy_gemm_nonlop_ompgpu() 1799 else if(dtset%gpu_option==ABI_GPU_LEGACY .or. dtset%gpu_option==ABI_GPU_KOKKOS) then 1800 call destroy_gemm_nonlop_gpu(dtset%nkpt) 1801 call destroy_gemm_nonlop(dtset%nkpt) 1802 else if(dtset%gpu_option==ABI_GPU_DISABLED) then 1803 call destroy_gemm_nonlop(dtset%nkpt) 1804 end if 1805 gemm_nonlop_use_gemm = .false. 1806 end if 1807 1808 !Clean GPU work spaces 1809 if(dtset%gpu_option == ABI_GPU_OPENMP) then 1810 call free_getghc_ompgpu_buffers() 1811 end if 1812 #if defined HAVE_GPU 1813 if (dtset%gpu_option/=ABI_GPU_DISABLED) then 1814 call dealloc_hamilt_gpu(2,dtset%gpu_option) 1815 end if 1816 #endif 1817 1818 #if defined HAVE_BIGDFT 1819 if (dtset%usewvl == 1 .and. dtset%timopt==10) then 1820 call wvl_timing(xmpi_world,'== WFN OPT','PR') 1821 end if 1822 #endif 1823 1824 call timab(1231,2,tsec) 1825 call timab(1232,2,tsec) 1826 1827 DBG_EXIT("COLL") 1828 1829 end subroutine gstate
m_gstate/pawuj_drive [ Functions ]
[ Top ] [ m_gstate ] [ Functions ]
NAME
pawuj_drive
FUNCTION
Drive for automatic determination of U Relevant only in PAW+U context
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 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) kg(3,mpw*mkmem)=reduced planewave coordinates. mpi_enreg=information about MPI parallelization nattyp(ntypat)= # atoms of each type. npwarr(nkpt)=number of planewaves in basis at this k point nspinor=number of spinorial components of the wavefunctions 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,mpw*nspinor*mband*mkmem*nsppol)=updated wavefunctions. dtefield <type(efield_type)> = variables related to Berry phase calculations (see initberry.f) eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation 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(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,wffnow=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
SOURCE
2542 subroutine pawuj_drive(scfcv_args, dtset,electronpositron,rhog,rhor,rprimd, xred,xred_old) 2543 2544 !Arguments ------------------------------------ 2545 !scalars 2546 type(scfcv_t), intent(inout) :: scfcv_args 2547 type(dataset_type),intent(inout) :: dtset 2548 type(electronpositron_type),pointer :: electronpositron 2549 !type(wffile_type),intent(inout) :: wffnew,wffnow 2550 !arrays 2551 real(dp), intent(inout) :: rprimd(3,3) 2552 real(dp), pointer :: rhog(:,:),rhor(:,:) 2553 real(dp), intent(inout) :: xred(3,dtset%natom),xred_old(3,dtset%natom) 2554 2555 !Local variables ------------------------- 2556 !scalars 2557 integer,parameter :: itime0 = 0 2558 integer,target :: ndtpawuj=4 2559 integer :: iuj,conv_retcode 2560 integer :: itimes(2) 2561 real(dp) :: ures 2562 !character(len=500) :: msg 2563 !arrays 2564 !real(dp),allocatable :: cgstart(:,:) 2565 type(macro_uj_type),allocatable,target :: dtpawuj(:) 2566 ! ********************************************************************* 2567 2568 2569 DBG_ENTER("COLL") 2570 2571 if (dtset%macro_uj==0) then 2572 ABI_BUG('Macro_uj must be set !') 2573 end if 2574 2575 ABI_MALLOC(dtpawuj,(0:ndtpawuj)) 2576 !ABI_MALLOC(cgstart,(2,scfcv_args%mcg)) 2577 2578 call pawuj_ini(dtpawuj,ndtpawuj) 2579 2580 !cgstart=scfcv_args%cg 2581 do iuj=1,ndtpawuj 2582 ! allocate(dtpawuj(iuj)%rprimd(3,3)) ! this has already been done in pawuj_ini 2583 dtpawuj(iuj)%macro_uj=dtset%macro_uj 2584 dtpawuj(iuj)%pawprtvol=dtset%pawprtvol 2585 dtpawuj(iuj)%diemix=dtset%diemix 2586 dtpawuj(iuj)%diemixmag=dtset%diemixmag 2587 dtpawuj(iuj)%pawujat=dtset%pawujat 2588 dtpawuj(iuj)%nspden=dtset%nspden 2589 dtpawuj(iuj)%rprimd=dtset%rprimd_orig(1:3,1:3,1) 2590 dtpawuj(iuj)%dmatpuopt=dtset%dmatpuopt 2591 end do 2592 2593 iuj=1 !LMac Flag to collect occupancies for unperturbed calculation 2594 dtpawuj(iuj)%iuj=iuj 2595 2596 scfcv_args%ndtpawuj=>ndtpawuj 2597 scfcv_args%dtpawuj=>dtpawuj 2598 2599 itimes(1)=itime0 ; itimes(2)=1 2600 call scfcv_run(scfcv_args,electronpositron,itimes,rhog,rhor,rprimd,xred,xred_old,conv_retcode) 2601 2602 !Calculate Hubbard U (or J) 2603 call pawuj_det(dtpawuj, ndtpawuj, dtset, scfcv_args%dtfil, ures, scfcv_args%mpi_enreg%comm_cell) 2604 dtset%upawu(dtset%typat(dtset%pawujat),1)=ures/Ha_eV 2605 2606 !Deallocations 2607 do iuj=0,ndtpawuj 2608 call pawuj_free(dtpawuj(iuj)) 2609 end do 2610 2611 ABI_FREE(dtpawuj) 2612 !ABI_FREE(cgstart) 2613 2614 DBG_EXIT("COLL") 2615 2616 end subroutine pawuj_drive
m_gstate/prtxf [ Functions ]
[ Top ] [ m_gstate ] [ Functions ]
NAME
prtxf
FUNCTION
Compute and print out dimensional cartesian coordinates and forces. Note: for x=cartesian coordinates, t=reduced coordinates (xred), => $ x= R t $ => $ x(1)=rprimd(1,1) t(1)+rprimd(2,1) t(2)+rprimd(3,1) t(3)$ etc. Also $ t = (R^{-1}) x$ . To convert gradients, $d(E)/dx(n) = [d(E)/dt(m)] [dt(m)/dx(n)]$ and $ dt(m)/dx(n) = (R^{-1})_{mn} = G_{nm}$ because G is the inverse transpose of R. Finally then $d(E)/dx(n) = G_{nm} [d(E)/dt(m)]$. The vector $d(E)/dt(m)$ for each atom is input in gred (grad. wrt xred).
INPUTS
gred(3,natom)=gradients of Etot (hartree) wrt xred(3,natom) iatfix(3,natom)=1 for each fixed atom along specified direction, else 0 iout=unit number for output file iwfrc=controls force output: 0=> no forces output, 1=>forces out in eV/A and Ha/bohr, 2=>forces out in Ha/bohr natom=number of atoms in unit cell rprimd(3,3)=dimensional real space primitive translations (bohr) xred(3,natom)=relative coordinates of atoms (in terms of prim. transl.)
OUTPUT
(data written to unit iout)
SOURCE
2158 subroutine prtxf(gred,iatfix,iout,iwfrc,natom,rprimd,xred) 2159 2160 !Arguments ------------------------------------ 2161 !scalars 2162 integer,intent(in) :: iout,iwfrc,natom 2163 !arrays 2164 integer,intent(in) :: iatfix(3,natom) 2165 real(dp),intent(in) :: gred(3,natom),rprimd(3,3),xred(3,natom) 2166 2167 !Local variables------------------------------- 2168 !scalars 2169 integer :: iatom,mu,unfixd 2170 real(dp) :: convt,fmax,frms 2171 character(len=15) :: format_line21 2172 character(len=15) :: format_line25 2173 character(len=15) :: format_line 2174 character(len=500) :: msg 2175 !arrays 2176 real(dp) :: favg(3),favg_out(3),ff(3),gprimd(3,3),xx(3) 2177 2178 ! **************************************************************** 2179 2180 format_line21='(i5,1x,3f21.14)' 2181 format_line25='(i5,1x,3f25.14)' 2182 2183 !Write cartesian coordinates in angstroms 2184 call wrtout(iout,' cartesian coordinates (angstrom) at end:') 2185 do iatom=1,natom 2186 format_line=format_line21 2187 do mu=1,3 2188 xx(mu)=(rprimd(mu,1)*xred(1,iatom)+& 2189 & rprimd(mu,2)*xred(2,iatom)+& 2190 & rprimd(mu,3)*xred(3,iatom))*Bohr_Ang 2191 if(xx(mu)>99999 .or. xx(mu)<-9999)format_line=format_line25 2192 end do 2193 write(msg,format_line) iatom,xx 2194 call wrtout(iout, msg) 2195 end do 2196 2197 !Optionally write cartesian forces in eV/Angstrom (also provide same in hartree/bohr) 2198 if (iwfrc/=0) then 2199 ! First, provide results in hartree/bohr 2200 write(msg, '(a,a)' ) ch10,' cartesian forces (hartree/bohr) at end:' 2201 call wrtout(iout, msg) 2202 frms=zero 2203 fmax=zero 2204 favg(1)=zero 2205 favg(2)=zero 2206 favg(3)=zero 2207 ! To get cartesian forces from input gradients with respect to 2208 ! dimensionless coordinates xred, multiply by G and negate 2209 ! (see notes at top of this subroutine) 2210 call matr3inv(rprimd,gprimd) 2211 ! First compute (spurious) average force favg 2212 do iatom=1,natom 2213 do mu=1,3 2214 ff(mu)=-(gprimd(mu,1)*gred(1,iatom)+& 2215 & gprimd(mu,2)*gred(2,iatom)+& 2216 & gprimd(mu,3)*gred(3,iatom)) 2217 favg(mu)=favg(mu)+ff(mu) 2218 end do 2219 end do 2220 favg(1) = favg(1)/dble(natom) 2221 favg(2) = favg(2)/dble(natom) 2222 favg(3) = favg(3)/dble(natom) 2223 2224 ! Subtract off average force in what follows 2225 ! (avg is also subtracted off in carfor, called by loopcv, 2226 ! called by grad) 2227 unfixd=0 2228 do iatom=1,natom 2229 format_line=format_line21 2230 do mu=1,3 2231 ff(mu)=-(gprimd(mu,1)*gred(1,iatom)+& 2232 & gprimd(mu,2)*gred(2,iatom)+& 2233 & gprimd(mu,3)*gred(3,iatom))-favg(mu) 2234 if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25 2235 ! For rms and max force, include only unfixed components 2236 if (iatfix(mu,iatom) /= 1) then 2237 unfixd=unfixd+1 2238 frms=frms+ff(mu)**2 2239 fmax=max(fmax,abs(ff(mu))) 2240 end if 2241 end do 2242 write(msg, format_line) iatom,ff 2243 call wrtout(iout, msg) 2244 end do 2245 if ( unfixd /= 0 ) frms = sqrt(frms/dble(unfixd)) 2246 2247 ! The average force is obtained from the cancellation of numbers 2248 ! of typical size unity, so an absolute value lower 2249 ! than tol14 is meaningless for the output file. 2250 favg_out(:)=favg(:) 2251 if(abs(favg_out(1))<tol14)favg_out(1)=zero 2252 if(abs(favg_out(2))<tol14)favg_out(2)=zero 2253 if(abs(favg_out(3))<tol14)favg_out(3)=zero 2254 2255 write(msg, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',frms,fmax,favg_out(1:3),' h/b' 2256 call wrtout(iout, msg) 2257 2258 if (iwfrc==1) then 2259 2260 write(msg, '(a,a)' )ch10,' cartesian forces (eV/Angstrom) at end:' 2261 call wrtout(iout, msg) 2262 convt=Ha_eV/Bohr_Ang 2263 2264 ! Note: subtract off average force 2265 do iatom=1,natom 2266 format_line=format_line21 2267 do mu=1,3 2268 ff(mu)=(-(gprimd(mu,1)*gred(1,iatom)+& 2269 & gprimd(mu,2)*gred(2,iatom)+& 2270 & gprimd(mu,3)*gred(3,iatom))-favg(mu))*convt 2271 if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25 2272 end do 2273 write(msg, format_line) iatom,ff 2274 call wrtout(iout, msg) 2275 end do 2276 write(msg, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',convt*frms,convt*fmax,convt*favg_out(1:3),' e/A' 2277 call wrtout(iout, msg) 2278 2279 end if 2280 end if 2281 2282 end subroutine prtxf
m_gstate/setup2 [ Functions ]
[ Top ] [ m_gstate ] [ Functions ]
NAME
setup2
FUNCTION
Call within main routine for setup of various arrays.
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset | ecut=kinetic energy cutoff for planewave basis (hartree) | natom=number of atoms in unit cell | nkpt=number of k points | wtk(nkpt)=integration weight associated with each k point | iscf=parameter controlling scf or non-scf choice npwtot(nkpt)=number of planewaves in basis and boundary at each k point xred(3,natom)=starting reduced atomic coordinates
OUTPUT
start(3,natom)=copy of starting xred
SOURCE
1855 subroutine setup2(dtset,npwtot,start,wfs,xred) 1856 1857 !Arguments ------------------------------------ 1858 !scalars 1859 type(dataset_type),intent(in) :: dtset 1860 type(wvl_wf_type),intent(in) :: wfs 1861 !arrays 1862 integer,intent(in) :: npwtot(dtset%nkpt) 1863 real(dp),intent(in) :: xred(3,dtset%natom) 1864 real(dp),intent(out) :: start(3,dtset%natom) 1865 1866 !Local variables------------------------------- 1867 !scalars 1868 integer :: ikpt,npw 1869 real(dp) :: arith,geom,wtknrm 1870 character(len=500) :: msg 1871 1872 ! ************************************************************************* 1873 1874 if (dtset%iscf>=0) then 1875 1876 ! Copy coordinates into array start 1877 start(:,:)=xred(:,:) 1878 1879 if (dtset%usewvl == 0) then 1880 ! Get average number of planewaves per k point: 1881 ! both arithmetic and GEOMETRIC averages are desired-- 1882 ! need geometric average to use method of Francis and Payne, 1883 ! J. Phys.: Condens. Matter 2, 4395-4404 (1990) [[cite:Francis1990]]. 1884 ! Also note: force k point wts to sum to 1 for this averaging. 1885 ! (wtk is not forced to add to 1 in a case with occopt=2) 1886 arith=zero 1887 geom=one 1888 wtknrm=zero 1889 do ikpt=1,dtset%nkpt 1890 npw=npwtot(ikpt) 1891 wtknrm=wtknrm+dtset%wtk(ikpt) 1892 arith=arith+npw*dtset%wtk(ikpt) 1893 geom=geom*npw**dtset%wtk(ikpt) 1894 end do 1895 1896 ! Enforce normalization of weights to 1 1897 arith=arith/wtknrm 1898 geom=geom**(1.0_dp/wtknrm) 1899 1900 end if 1901 1902 ! Ensure portability of output thanks to tol8 1903 if (dtset%usewvl == 0) then 1904 write(msg, '(a,2f12.3)' ) '_setup2: Arith. and geom. avg. npw (full set) are',arith+tol8,geom 1905 else 1906 #if defined HAVE_BIGDFT 1907 write(msg, '(a,2I8)' ) ' setup2: nwvl coarse and fine are', & 1908 & wfs%ks%lzd%Glr%wfd%nvctr_c, wfs%ks%lzd%Glr%wfd%nvctr_f 1909 #endif 1910 end if 1911 call wrtout([std_out, ab_out], msg) 1912 end if 1913 1914 #if !defined HAVE_BIGDFT 1915 if (.false.) write(std_out,*) wfs%ks 1916 #endif 1917 1918 end subroutine setup2