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