TABLE OF CONTENTS
ABINIT/clnup1 [ 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) fnameabo_dos=filename of output DOS file fnameabo_eig=filename of output EIG file fred(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)
PARENTS
gstate
CHILDREN
getnel,metric,prteigrs,prtrhomxmn,prtxf,write_eig,wrtout
SOURCE
1844 subroutine clnup1(acell,dtset,eigen,fermie,& 1845 & fnameabo_dos,fnameabo_eig,fred,& 1846 & mpi_enreg,nfft,ngfft,occ,prtfor,& 1847 & resid,rhor,rprimd,vxcavg,xred) 1848 1849 !use defs_wvltypes 1850 1851 !This section has been created automatically by the script Abilint (TD). 1852 !Do not modify the following lines by hand. 1853 #undef ABI_FUNC 1854 #define ABI_FUNC 'clnup1' 1855 use interfaces_14_hidewrite 1856 use interfaces_59_ionetcdf 1857 use interfaces_67_common 1858 !End of the abilint section 1859 1860 implicit none 1861 1862 !Arguments ------------------------------------ 1863 !scalars 1864 integer,intent(in) :: nfft 1865 integer,intent(in) :: prtfor 1866 real(dp),intent(in) :: fermie 1867 real(dp),intent(in) :: vxcavg 1868 character(len=*),intent(in) :: fnameabo_dos,fnameabo_eig 1869 type(dataset_type),intent(in) :: dtset 1870 type(MPI_type),intent(in) :: mpi_enreg 1871 !arrays 1872 integer,intent(in) :: ngfft(18) 1873 real(dp),intent(in) :: acell(3) 1874 real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 1875 real(dp),intent(in) :: fred(3,dtset%natom) 1876 real(dp),intent(in) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol) 1877 real(dp),intent(in) :: rhor(nfft,dtset%nspden) 1878 real(dp),intent(in) :: rprimd(3,3) 1879 real(dp),intent(in) :: xred(3,dtset%natom) 1880 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 1881 1882 !Local variables------------------------------- 1883 !scalars 1884 integer,parameter :: master=0 1885 integer :: comm,iatom,ii,iscf_dum,iwfrc,me,nnonsc,option,unitdos 1886 real(dp) :: entropy,grmax,grsum,maxocc,nelect,tolwf,ucvol 1887 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 1888 character(len=500) :: message 1889 character(len=fnlen) filename 1890 !arrays 1891 real(dp),allocatable :: doccde(:) 1892 1893 ! **************************************************************** 1894 1895 comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm) 1896 1897 if(dtset%prtstm==0)then ! Write reduced coordinates xred 1898 write(message, '(a,i5,a)' )' reduced coordinates (array xred) for',dtset%natom,' atoms' 1899 call wrtout(ab_out,message,'COLL') 1900 do iatom=1,dtset%natom 1901 write(message, '(1x,3f20.12)' ) xred(:,iatom) 1902 call wrtout(ab_out,message,'COLL') 1903 end do 1904 end if 1905 1906 !Write reduced gradients if iscf > 0 and dtset%nstep>0 and 1907 if (dtset%iscf>=0.and.dtset%nstep>0.and.dtset%prtstm==0) then 1908 1909 ! Compute absolute maximum and root mean square value of gradients 1910 grmax=0.0_dp 1911 grsum=0.0_dp 1912 do iatom=1,dtset%natom 1913 do ii=1,3 1914 ! To be activated in v5.5 1915 ! grmax=max(grmax,abs(fred(ii,iatom))) 1916 grmax=max(grmax,fred(ii,iatom)) 1917 grsum=grsum+fred(ii,iatom)**2 1918 end do 1919 end do 1920 grsum=sqrt(grsum/dble(3*dtset%natom)) 1921 1922 write(message, '(1x,a,1p,e12.4,a,e12.4,a)' )'rms dE/dt=',grsum,'; max dE/dt=',grmax,'; dE/dt below (all hartree)' 1923 call wrtout(ab_out,message,'COLL') 1924 do iatom=1,dtset%natom 1925 write(message, '(i5,1x,3f20.12)' ) iatom,fred(1:3,iatom) 1926 call wrtout(ab_out,message,'COLL') 1927 end do 1928 1929 end if 1930 1931 if(dtset%prtstm==0)then 1932 1933 ! Compute and write out dimensional cartesian coords and forces: 1934 call wrtout(ab_out,' ','COLL') 1935 1936 ! (only write forces if iscf > 0 and dtset%nstep>0) 1937 if (dtset%iscf<0.or.dtset%nstep<=0.or.prtfor==0) then 1938 iwfrc=0 1939 else 1940 iwfrc=1 1941 end if 1942 1943 call prtxf(fred,dtset%iatfix,ab_out,iwfrc,dtset%natom,rprimd,xred) 1944 1945 ! Write length scales 1946 write(message, '(1x,a,3f16.12,a)' )'length scales=',acell,' bohr' 1947 call wrtout(ab_out,message,'COLL') 1948 write(message, '(14x,a,3f16.12,a)' )'=',Bohr_Ang*acell(1:3),' angstroms' 1949 call wrtout(ab_out,message,'COLL') 1950 1951 end if 1952 1953 option=1; nnonsc=0; tolwf=0.0_dp 1954 1955 if(dtset%iscf<0 .and. dtset%iscf/=-3)option=3 1956 iscf_dum=dtset%iscf 1957 if(dtset%nstep==0)iscf_dum=-1 1958 1959 if(dtset%tfkinfunc==0)then 1960 call prteigrs(eigen,dtset%enunit,fermie,fnameabo_eig,ab_out,& 1961 & iscf_dum,dtset%kptns,dtset%kptopt,dtset%mband,& 1962 & dtset%nband,dtset%nkpt,nnonsc,dtset%nsppol,occ,& 1963 & dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwf,& 1964 & vxcavg,dtset%wtk) 1965 call prteigrs(eigen,dtset%enunit,fermie,fnameabo_eig,std_out,& 1966 & iscf_dum,dtset%kptns,dtset%kptopt,dtset%mband,& 1967 & dtset%nband,dtset%nkpt,nnonsc,dtset%nsppol,occ,& 1968 & dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwf,& 1969 & vxcavg,dtset%wtk) 1970 1971 #if defined HAVE_NETCDF 1972 if (dtset%prteig==1 .and. me == master) then 1973 filename=trim(fnameabo_eig)//'.nc' 1974 call write_eig(eigen,filename,dtset%kptns,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol) 1975 end if 1976 #endif 1977 1978 end if 1979 1980 !Compute and print location of maximal and minimal density 1981 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1982 call prtrhomxmn(std_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol) 1983 if( dtset%prtvol>1)then 1984 call prtrhomxmn(ab_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol) 1985 end if 1986 1987 !If needed, print DOS (unitdos is closed in getnel, occ is not changed if option == 2 1988 if (dtset%prtdos==1 .and. me == master) then 1989 if (open_file(fnameabo_dos,message, newunit=unitdos, status='unknown', action="write", form='formatted') /= 0) then 1990 MSG_ERROR(message) 1991 end if 1992 rewind(unitdos) 1993 maxocc=two/(dtset%nspinor*dtset%nsppol) ! Will not work in the fixed moment case 1994 option=2 1995 ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1996 call getnel(doccde,dtset%dosdeltae,eigen,entropy,fermie,& 1997 & maxocc,dtset%mband,dtset%nband,nelect,dtset%nkpt,& 1998 & dtset%nsppol,occ,dtset%occopt,option,dtset%tphysel,& 1999 & dtset%tsmear,unitdos,dtset%wtk) 2000 ABI_DEALLOCATE(doccde) 2001 end if 2002 2003 end subroutine clnup1
ABINIT/clnup2 [ 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
fred(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)
PARENTS
gstate
CHILDREN
wrtout
SOURCE
2045 subroutine clnup2(n1xccc,fred,grchempottn,gresid,grewtn,grvdw,grxc,iscf,natom,ngrvdw,& 2046 & prtfor,prtstr,prtvol,start,strten,synlgr,xred) 2047 2048 2049 !This section has been created automatically by the script Abilint (TD). 2050 !Do not modify the following lines by hand. 2051 #undef ABI_FUNC 2052 #define ABI_FUNC 'clnup2' 2053 use interfaces_14_hidewrite 2054 !End of the abilint section 2055 2056 implicit none 2057 2058 !Arguments ------------------------------------ 2059 !scalars 2060 integer,intent(in) :: iscf,n1xccc,natom,ngrvdw,prtfor,prtstr,prtvol 2061 !arrays 2062 real(dp),intent(in) :: fred(3,natom),grchempottn(3,natom),gresid(3,natom) 2063 real(dp),intent(in) :: grewtn(3,natom),grvdw(3,ngrvdw) 2064 real(dp),intent(in) :: grxc(3,natom),start(3,natom),strten(6),synlgr(3,natom) 2065 real(dp),intent(in) :: xred(3,natom) 2066 2067 !Local variables------------------------------- 2068 character(len=*), parameter :: format01020 ="(i5,1x,3f20.12)" 2069 !scalars 2070 integer :: iatom,mu 2071 real(dp) :: devsqr,grchempot2 2072 character(len=500) :: message 2073 2074 ! ************************************************************************* 2075 ! 2076 !DEBUG 2077 !write(std_out,*)' clnup2 : enter ' 2078 !ENDDEBUG 2079 2080 !Only print additional info for scf calculations 2081 if (iscf>=0) then 2082 2083 if((prtvol>=10).and.(prtfor>0))then 2084 2085 write(message, '(a,10x,a)' ) ch10,& 2086 & '===> extra information on forces <===' 2087 call wrtout(ab_out,message,'COLL') 2088 2089 write(message, '(a)' ) ' ewald contribution to reduced grads' 2090 call wrtout(ab_out,message,'COLL') 2091 do iatom=1,natom 2092 write(message,format01020) iatom,(grewtn(mu,iatom),mu=1,3) 2093 call wrtout(ab_out,message,'COLL') 2094 end do 2095 2096 grchempot2=sum(grchempottn(:,:)**2) 2097 if(grchempot2>tol16)then 2098 write(message, '(a)' ) ' chemical potential contribution to reduced grads' 2099 call wrtout(ab_out,message,'COLL') 2100 do iatom=1,natom 2101 write(message,format01020) iatom,(grchempottn(mu,iatom),mu=1,3) 2102 call wrtout(ab_out,message,'COLL') 2103 end do 2104 end if 2105 2106 write(message, '(a)' ) ' nonlocal contribution to red. grads' 2107 call wrtout(ab_out,message,'COLL') 2108 do iatom=1,natom 2109 write(message,format01020) iatom,(synlgr(mu,iatom),mu=1,3) 2110 call wrtout(ab_out,message,'COLL') 2111 end do 2112 2113 write(message, '(a)' ) ' local psp contribution to red. grads' 2114 call wrtout(ab_out,message,'COLL') 2115 if (n1xccc/=0) then 2116 do iatom=1,natom 2117 write(message,format01020) iatom,fred(:,iatom)-& 2118 & (grewtn(:,iatom)+grchempottn(:,iatom)+synlgr(:,iatom)+grxc(:,iatom)+gresid(:,iatom)) 2119 call wrtout(ab_out,message,'COLL') 2120 end do 2121 else 2122 do iatom=1,natom 2123 write(message,format01020) iatom,fred(:,iatom)-& 2124 & (grewtn(:,iatom)+grchempottn(:,iatom)+synlgr(:,iatom)+gresid(:,iatom)) 2125 call wrtout(ab_out,message,'COLL') 2126 end do 2127 end if 2128 2129 if (n1xccc/=0) then 2130 write(message, '(a)' ) ' core charge xc contribution to reduced grads' 2131 call wrtout(ab_out,message,'COLL') 2132 do iatom=1,natom 2133 write(message,format01020) iatom,(grxc(mu,iatom),mu=1,3) 2134 call wrtout(ab_out,message,'COLL') 2135 end do 2136 end if 2137 2138 if (ngrvdw==natom) then 2139 write(message, '(a)' ) ' Van der Waals DFT-D contribution to reduced grads' 2140 call wrtout(ab_out,message,'COLL') 2141 do iatom=1,natom 2142 write(message,format01020) iatom,(grvdw(mu,iatom),mu=1,3) 2143 call wrtout(ab_out,message,'COLL') 2144 end do 2145 end if 2146 2147 write(message, '(a)' ) ' residual contribution to red. grads' 2148 call wrtout(ab_out,message,'COLL') 2149 do iatom=1,natom 2150 write(message,format01020) iatom,(gresid(mu,iatom),mu=1,3) 2151 call wrtout(ab_out,message,'COLL') 2152 end do 2153 2154 end if 2155 2156 ! Compute mean squared deviation from starting coords 2157 devsqr=0.0_dp 2158 do iatom=1,natom 2159 do mu=1,3 2160 devsqr=devsqr+(xred(mu,iatom)-start(mu,iatom))**2 2161 end do 2162 end do 2163 2164 ! When shift is nonnegligible then print values 2165 if (devsqr>1.d-14) then 2166 write(message, '(a,1p,e12.4,3x,a)' ) & 2167 & ' rms coord change=',sqrt(devsqr/dble(3*natom)),& 2168 & 'atom, delta coord (reduced):' 2169 call wrtout(ab_out,message,'COLL') 2170 do iatom=1,natom 2171 write(message, '(1x,i5,2x,3f20.12)' ) iatom,& 2172 & (xred(mu,iatom)-start(mu,iatom),mu=1,3) 2173 call wrtout(ab_out,message,'COLL') 2174 end do 2175 end if 2176 2177 ! Write out stress results 2178 if (prtstr>0) then 2179 write(message, '(a,a)' ) ch10,& 2180 & ' Cartesian components of stress tensor (hartree/bohr^3)' 2181 call wrtout(ab_out,message,'COLL') 2182 call wrtout(std_out, message,'COLL') 2183 2184 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2185 & ' sigma(1 1)=',strten(1),' sigma(3 2)=',strten(4) 2186 call wrtout(ab_out,message,'COLL') 2187 call wrtout(std_out, message,'COLL') 2188 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2189 & ' sigma(2 2)=',strten(2),' sigma(3 1)=',strten(5) 2190 call wrtout(ab_out,message,'COLL') 2191 call wrtout(std_out, message,'COLL') 2192 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2193 & ' sigma(3 3)=',strten(3),' sigma(2 1)=',strten(6) 2194 call wrtout(ab_out,message,'COLL') 2195 call wrtout(std_out, message,'COLL') 2196 2197 ! Also output the pressure (minus one third the trace of the stress 2198 ! tensor. 2199 write(message, '(a,a,es12.4,a)' ) ch10,& 2200 & '-Cartesian components of stress tensor (GPa) [Pressure=',& 2201 & -(strten(1)+strten(2)+strten(3))*HaBohr3_GPa/3.0_dp,' GPa]' 2202 2203 call wrtout(ab_out,message,'COLL') 2204 call wrtout(std_out, message,'COLL') 2205 2206 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2207 & '- sigma(1 1)=',strten(1)*HaBohr3_GPa,& 2208 & ' sigma(3 2)=',strten(4)*HaBohr3_GPa 2209 call wrtout(ab_out,message,'COLL') 2210 call wrtout(std_out, message,'COLL') 2211 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2212 & '- sigma(2 2)=',strten(2)*HaBohr3_GPa,& 2213 & ' sigma(3 1)=',strten(5)*HaBohr3_GPa 2214 call wrtout(ab_out,message,'COLL') 2215 call wrtout(std_out, message,'COLL') 2216 write(message, '(a,1p,e16.8,a,1p,e16.8)' ) & 2217 & '- sigma(3 3)=',strten(3)*HaBohr3_GPa,& 2218 & ' sigma(2 1)=',strten(6)*HaBohr3_GPa 2219 call wrtout(ab_out,message,'COLL') 2220 call wrtout(std_out, message,'COLL') 2221 end if 2222 2223 ! Last end if above refers to iscf > 0 2224 end if 2225 2226 !DEBUG 2227 !write(std_out,*)' clnup2 : exit ' 2228 !ENDDEBUG 2229 2230 end subroutine clnup2
ABINIT/gstate [ Functions ]
NAME
gstate
FUNCTION
Primary routine for conducting DFT calculations by CG minimization.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JYR, MKV, MT, FJ, MB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
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
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
PARENTS
gstateimg
CHILDREN
wrtout
SOURCE
108 #if defined HAVE_CONFIG_H 109 #include "config.h" 110 #endif 111 112 #include "abi_common.h" 113 114 subroutine gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,initialized,& 115 & mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,& 116 & psps,results_gs,rprim,scf_history,vel,vel_cell,wvl,xred) 117 118 use defs_basis 119 use defs_datatypes 120 use defs_abitypes 121 use defs_parameters 122 use defs_rectypes 123 use m_errors 124 use m_xmpi 125 use m_profiling_abi 126 use libxc_functionals 127 use m_exit 128 use m_crystal 129 use m_crystal_io 130 use m_scf_history 131 use m_abimover 132 use m_wffile 133 use m_rec 134 use m_efield 135 use m_orbmag 136 use m_ddb 137 use m_bandfft_kpt 138 use m_invovl 139 use m_gemm_nonlop 140 use m_tetrahedron 141 use m_wfk 142 use m_nctk 143 use m_hdr 144 use m_ebands 145 146 use m_io_tools, only : open_file 147 use m_occ, only : newocc, getnel 148 use m_ddb_hdr, only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write 149 use m_fstrings, only : strcat, sjoin 150 use m_geometry, only : fixsym, mkradim, metric 151 use m_kpts, only : tetra_from_kptrlatt 152 use m_kg, only : kpgio, getph 153 use m_pawang, only : pawang_type 154 use m_pawrad, only : pawrad_type 155 use m_pawtab, only : pawtab_type 156 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 157 use m_abi2big, only : wvl_occ_abi2big 158 use m_energies, only : energies_type, energies_init 159 use m_args_gs, only : args_gs_type 160 use m_results_gs, only : results_gs_type 161 use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_free 162 use m_paw_dmft, only : init_sc_dmft,destroy_sc_dmft,print_sc_dmft,paw_dmft_type,readocc_dmft 163 use m_data4entropyDMFT, only : data4entropyDMFT_t, data4entropyDMFT_init, data4entropyDMFT_destroy 164 use m_electronpositron, only : electronpositron_type,init_electronpositron,destroy_electronpositron, & 165 & electronpositron_calctype 166 use m_scfcv, only : scfcv_t,scfcv_init, scfcv_destroy, scfcv_run 167 use m_iowf, only : outwf 168 use m_outqmc, only : outqmc 169 use m_ioarr, only : ioarr,read_rhor 170 use defs_wvltypes, only : wvl_data,coulomb_operator,wvl_wf_type 171 #if defined HAVE_BIGDFT 172 use BigDFT_API, only : wvl_timing => timing,xc_init,xc_end,XC_MIXED,XC_ABINIT,& 173 & local_potential_dimensions,nullify_gaussian_basis, & 174 & copy_coulomb_operator,deallocate_coulomb_operator 175 #else 176 use defs_wvltypes, only : coulomb_operator 177 #endif 178 179 #if defined HAVE_LOTF 180 use defs_param_lotf, only : lotfparam_init 181 #endif 182 183 !This section has been created automatically by the script Abilint (TD). 184 !Do not modify the following lines by hand. 185 #undef ABI_FUNC 186 #define ABI_FUNC 'gstate' 187 use interfaces_14_hidewrite 188 use interfaces_18_timing 189 use interfaces_32_util 190 use interfaces_43_wvl_wrappers 191 #if defined HAVE_GPU_CUDA 192 use interfaces_52_manage_cuda 193 #endif 194 use interfaces_53_ffts 195 use interfaces_56_io_mpi 196 use interfaces_56_recipspace 197 use interfaces_62_poisson 198 use interfaces_64_psp 199 use interfaces_65_paw 200 use interfaces_67_common 201 use interfaces_79_seqpar_mpi 202 use interfaces_95_drive, except_this_one => gstate 203 !End of the abilint section 204 205 implicit none 206 207 !Arguments ------------------------------------ 208 !scalars 209 integer,intent(inout) :: iexit,initialized 210 real(dp),intent(in) :: cpui 211 character(len=6),intent(in) :: codvsn 212 type(MPI_type),intent(inout) :: mpi_enreg 213 type(args_gs_type),intent(in) :: args_gs 214 type(datafiles_type),intent(inout) :: dtfil 215 type(dataset_type),intent(inout) :: dtset 216 type(pawang_type),intent(inout) :: pawang 217 type(pseudopotential_type),intent(inout) :: psps 218 type(results_gs_type),intent(inout) :: results_gs 219 type(scf_history_type),target,intent(inout) :: scf_history 220 type(wvl_data),intent(inout) :: wvl 221 !arrays 222 integer,intent(out) :: npwtot(dtset%nkpt) 223 real(dp),intent(inout) :: acell(3),occ(dtset%mband*dtset%nkpt*dtset%nsppol) 224 real(dp),intent(inout) :: rprim(3,3),vel(3,dtset%natom),vel_cell(3,3),xred(3,dtset%natom) 225 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 226 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 227 228 !Local variables------------------------------- 229 !Define file format for different type of files. Presently, 230 !only one file format is supported for each type of files, but this might change soon ... 231 !2 for wavefunction file, new format (version 2.0 and after) (fform) NOT USED 232 !52 for density rho(r) (fformr) 233 !102 for potential V(r) file. (fformv) NOT USED 234 !scalars 235 integer,parameter :: formeig=0,level=101,response=0 ,cplex1=1, master=0 236 integer :: ndtpawuj=0 ! Cannot use parameter because scfargs points to this! Have to get rid of pointers to scalars! 237 #if defined HAVE_BIGDFT 238 integer :: icoulomb 239 #endif 240 integer :: accessfil,ask_accurate,bantot,choice,comm_psp,fform 241 integer :: gnt_option,gscase,iatom,idir,ierr,ii,indx,jj,kk,ios,itypat 242 integer :: ixfh,izero,mcg,me,mgfftf,mpert,msize,mu,my_natom,my_nspinor 243 integer :: nblok,nfftf,nfftot,npwmin 244 integer :: openexit,option,optorth,psp_gencond,conv_retcode 245 integer :: pwind_alloc,rdwrpaw,comm,tim_mkrho,use_sc_dmft 246 integer :: cnt,spin,band,ikpt,usecg 247 real(dp) :: cpus,ecore,ecut_eff,ecutdg_eff,etot,fermie 248 real(dp) :: gsqcut_eff,gsqcut_shp,gsqcutc_eff,hyb_range_fock,residm,ucvol 249 logical :: read_wf_or_den,has_to_init,call_pawinit,write_wfk 250 logical :: wvlbigdft=.false.,wvl_debug=.false. 251 character(len=500) :: message 252 character(len=fnlen) :: ddbnm,dscrpt,filnam,wfkfull_path 253 real(dp) :: fatvshift 254 type(crystal_t) :: cryst 255 type(ebands_t) :: bstruct,ebands 256 type(t_tetrahedron) :: tetra 257 type(efield_type) :: dtefield 258 type(electronpositron_type),pointer :: electronpositron 259 type(hdr_type) :: hdr,hdr_den 260 type(macro_uj_type) :: dtpawuj(0) 261 type(orbmag_type) :: dtorbmag 262 type(paw_dmft_type) :: paw_dmft 263 type(pawfgr_type) :: pawfgr 264 type(recursion_type) ::rec_set 265 type(wffile_type) :: wff1,wffnew,wffnow 266 type(ab_xfh_type) :: ab_xfh 267 type(ddb_type) :: ddb 268 type(ddb_hdr_type) :: ddb_hdr 269 type(scfcv_t) :: scfcv_args 270 !arrays 271 integer :: ngfft(18),ngfftf(18) 272 integer,allocatable :: atindx(:),atindx1(:),indsym(:,:,:) 273 integer,allocatable :: irrzon(:,:,:),kg(:,:),nattyp(:),symrec(:,:,:) 274 integer,allocatable,target :: npwarr(:) 275 integer,pointer :: npwarr_(:),pwind(:,:,:) 276 real(dp) :: efield_band(3),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3) 277 real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3),tsec(2) 278 real(dp),allocatable :: amass(:),cg(:,:),doccde(:) 279 real(dp),allocatable :: eigen(:),ph1df(:,:),phnons(:,:,:),resid(:),rhowfg(:,:) 280 real(dp),allocatable :: rhowfr(:,:),spinat_dum(:,:),start(:,:),work(:) 281 real(dp),allocatable :: ylm(:,:),ylmgr(:,:,:) 282 real(dp),pointer :: pwnsfac(:,:),rhog(:,:),rhor(:,:) 283 real(dp),pointer :: taug(:,:),taur(:,:),xred_old(:,:) 284 type(pawrhoij_type),pointer :: pawrhoij(:) 285 type(coulomb_operator) :: kernel_dummy 286 287 ! *********************************************************************** 288 289 DBG_ENTER("COLL") 290 291 call timab(32,1,tsec) 292 call timab(33,3,tsec) 293 294 call status(0,dtfil%filstat,iexit,level,'enter') 295 296 !########################################################### 297 !### 01. Initializations XML, MPI, WVL, etc 298 299 !Init MPI data 300 comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm) 301 302 !Set up MPI information from the dataset 303 my_natom=mpi_enreg%my_natom 304 305 !Set up information when wavelets are in use 306 if (dtset%usewvl == 1) then 307 308 ! If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed 309 wvlbigdft=(dtset%wvl_bigdft_comp==1) 310 311 ! Default value, to be set-up elsewhere. 312 wvl%descr%h(:) = dtset%wvl_hgrid 313 314 #if defined HAVE_BIGDFT 315 wvl%descr%paw%usepaw=psps%usepaw 316 wvl%descr%paw%natom=dtset%natom 317 #endif 318 319 ! We set the atom-related internal wvl variables. 320 call wvl_descr_atoms_set(acell, dtset%icoulomb, dtset%natom, & 321 & dtset%ntypat, dtset%typat, wvl%descr) 322 if(dtset%usepaw==0) then 323 ! nullify PAW proj_G in NC case: 324 #if defined HAVE_BIGDFT 325 ABI_DATATYPE_ALLOCATE(wvl%projectors%G,(dtset%ntypat)) 326 do itypat=1,dtset%ntypat 327 call nullify_gaussian_basis(wvl%projectors%G(itypat)) 328 end do 329 #endif 330 end if 331 332 wvl%descr%exctxpar = "OP2P" 333 end if 334 335 if (me == master .and. dtset%prtxml == 1) then 336 ! gstate() will handle a dataset, so we output the dataSet markup. 337 write(ab_xml_out, "(A)") ' <dataSet>' 338 end if 339 340 !Define FFT grid(s) sizes (be careful !) 341 !See NOTES in the comments at the beginning of this file. 342 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf) 343 344 !Structured debugging if prtvol==-level 345 if(dtset%prtvol==-level)then 346 write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' gstate : enter , debug mode ' 347 call wrtout(std_out,message,'COLL') 348 end if 349 350 !########################################################### 351 !### 02. Calls setup1, kpgio, initylmg 352 353 ecore=zero 354 results_gs%pel(1:3) =zero 355 results_gs%grchempottn(:,:)=zero 356 results_gs%grewtn(:,:)=zero 357 !MT Feb 2012: I dont know why but grvdw has to be allocated 358 !when using BigDFT to ensure success on inca_gcc44_sdebug 359 if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).or.dtset%usewvl==1) then 360 results_gs%ngrvdw=dtset%natom 361 if (allocated(results_gs%grvdw)) then 362 ABI_DEALLOCATE(results_gs%grvdw) 363 end if 364 ABI_ALLOCATE(results_gs%grvdw,(3,dtset%natom)) 365 results_gs%grvdw(:,:)=zero 366 end if 367 call energies_init(results_gs%energies) 368 369 !Set up for iterations 370 ABI_ALLOCATE(amass,(dtset%natom)) 371 call setup1(acell,amass,args_gs%amu,bantot,dtset,& 372 & ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,& 373 & dtset%natom,ngfftf,ngfft,dtset%nkpt,dtset%nsppol,& 374 & response,rmet,rprim,rprimd,ucvol,psps%usepaw) 375 376 !In some cases (e.g. getcell/=0), the plane wave vectors have 377 ! to be generated from the original simulation cell 378 rprimd_for_kg=rprimd 379 if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=args_gs%rprimd_orig 380 call matr3inv(rprimd_for_kg,gprimd_for_kg) 381 gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg) 382 383 !Set up the basis sphere of planewaves 384 ABI_ALLOCATE(npwarr,(dtset%nkpt)) 385 if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then 386 ABI_ALLOCATE(kg,(3,dtset%mpw*dtset%mkmem)) 387 call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg, & 388 & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,& 389 & dtset%mpw,npwarr,npwtot,dtset%nsppol) 390 call bandfft_kpt_init1(bandfft_kpt,dtset%istwfk,kg,dtset%mgfft,dtset%mkmem,mpi_enreg,& 391 & dtset%mpw,dtset%nband,dtset%nkpt,npwarr,dtset%nsppol) 392 else 393 ABI_ALLOCATE(kg,(0,0)) 394 npwarr(:) = 0 395 npwtot(:) = 0 396 end if 397 398 if(dtset%wfoptalg == 1 .and. psps%usepaw == 1) then 399 call init_invovl(dtset%nkpt) 400 end if 401 402 if(dtset%use_gemm_nonlop == 1 .and. dtset%use_gpu_cuda/=1) then 403 ! set global variable 404 gemm_nonlop_use_gemm = .true. 405 call init_gemm_nonlop(dtset%nkpt) 406 else 407 gemm_nonlop_use_gemm = .false. 408 end if 409 410 !Set up the Ylm for each k point 411 if ( dtset%tfkinfunc /= 2) then 412 ABI_ALLOCATE(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)) 413 ABI_ALLOCATE(ylmgr,(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)) 414 if (psps%useylm==1) then 415 option=0 416 if (dtset%prtstm==0.and.dtset%iscf>0.and.dtset%positron/=1) option=1 ! compute gradients of YLM 417 if (dtset%berryopt==4 .and. dtset%optstress /= 0 .and. psps%usepaw==1) option = 1 ! compute gradients of YLM 418 call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,& 419 & psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,& 420 & npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr) 421 end if 422 else 423 ABI_ALLOCATE(ylm,(0,0)) 424 ABI_ALLOCATE(ylmgr,(0,0,0)) 425 end if 426 427 !SCF history management (allocate it at first call) 428 has_to_init=(initialized==0.or.scf_history%history_size<0) 429 if (initialized==0) then 430 ! This call has to be done before any use of SCF history 431 usecg=0 432 if(dtset%extrapwf>0)usecg=1 433 call scf_history_init(dtset,mpi_enreg,usecg,scf_history) 434 end if 435 436 call timab(33,2,tsec) 437 call timab(701,3,tsec) 438 439 !########################################################### 440 !### 03. Calls pspini 441 442 !Open and read pseudopotential files 443 comm_psp=mpi_enreg%comm_cell;if (dtset%usewvl==1) comm_psp=mpi_enreg%comm_wvl 444 if (dtset%nimage>1) psps%mixalch(:,:)=args_gs%mixalch(:,:) ! mixalch can evolve for some image algos 445 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,& 446 & pawrad,pawtab,psps,rprimd,comm_mpi=comm_psp) 447 448 call timab(701,2,tsec) 449 call timab(33,3,tsec) 450 if (psps%usepaw==1.and.dtset%usefock==1)then 451 do itypat=1,dtset%ntypat 452 if(pawtab(itypat)%has_fock==0) then 453 write(message, '(a)' )& 454 & 'The PAW data file does not contain Fock information. Change the PAW data file.' 455 MSG_BUG(message) 456 end if 457 end do 458 end if 459 !In case of isolated computations, ecore must set to zero 460 !because its contribution is counted in the ewald energy as the ion-ion interaction. 461 if (dtset%icoulomb == 1) ecore = zero 462 463 !WVL - Now that psp data are available, we compute rprimd, acell... from the atomic positions. 464 if (dtset%usewvl == 1) then 465 call wvl_descr_psp_set(trim(dtfil%filnam_ds(3))//"_OCCUP",dtset%nsppol,psps,dtset%spinat,wvl%descr) 466 call wvl_setBoxGeometry(dtset%prtvol, psps%gth_params%radii_cf, rprimd, xred, & 467 & wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult) 468 call mkradim(acell,rprim,rprimd) 469 rprimd_for_kg=rprimd 470 call wvl_denspot_set(wvl%den, psps%gth_params, dtset%ixc, dtset%natom, dtset%nsppol, rprimd, & 471 & wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult, mpi_enreg%comm_wvl, xred) 472 ! TODO: to be moved in a routine. 473 #if defined HAVE_BIGDFT 474 if (wvl%descr%atoms%astruct%geocode == "F") then 475 icoulomb = 1 476 else if (wvl%descr%atoms%astruct%geocode == "S") then 477 icoulomb = 2 478 else 479 icoulomb = 0 480 end if 481 ! calculation of the Poisson kernel anticipated to reduce memory peak for small systems 482 call psolver_kernel( wvl%den%denspot%dpbox%hgrids, 1, icoulomb, mpi_enreg%me_wvl, wvl%den%denspot%pkernel , & 483 & mpi_enreg%comm_wvl, wvl%den%denspot%dpbox%ndims, mpi_enreg%nproc_wvl, dtset%nscforder) 484 nullify(wvl%den%denspot%pkernelseq%kernel) 485 !call copy_coulomb_operator(wvl%den%denspot%pkernel,wvl%den%denspot%pkernelseq,ABI_FUNC) 486 ! Associate the denspot distribution into mpi_enreg. 487 mpi_enreg%nscatterarr => wvl%den%denspot%dpbox%nscatterarr 488 mpi_enreg%ngatherarr => wvl%den%denspot%dpbox%ngatherarr 489 mpi_enreg%ngfft3_ionic = wvl%den%denspot%dpbox%n3pi 490 call wvl_setngfft(mpi_enreg%me_wvl, dtset%mgfft, dtset%nfft, & 491 & dtset%ngfft, mpi_enreg%nproc_wvl, wvl%den%denspot%dpbox%ndims(1), & 492 & wvl%den%denspot%dpbox%ndims(2),wvl%den%denspot%dpbox%ndims(3),& 493 & wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl, 1)) 494 #endif 495 nfftf = dtset%nfft 496 mgfftf = dtset%mgfft 497 ngfftf(:) = dtset%ngfft(:) 498 ! Recalculate gprimd 499 call matr3inv(rprimd,gprimd) 500 ! PAW section 501 if(psps%usepaw==1) then 502 ! Reinitialize Pawfgr with new values of ngfft 503 call pawfgr_destroy(pawfgr) 504 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf) 505 ! fill wvl objects from paw objects 506 ! wvl%descr%paw%usepaw=dtset%usepaw 507 call paw2wvl(pawtab,wvl%projectors,wvl%descr) 508 end if 509 else if (dtset%icoulomb /= 0) then 510 #if defined HAVE_BIGDFT 511 if (dtset%ixc < 0) then 512 call xc_init(wvl%den%denspot%xc, dtset%ixc, XC_MIXED, dtset%nsppol) 513 else 514 call xc_init(wvl%den%denspot%xc, dtset%ixc, XC_ABINIT, dtset%nsppol) 515 end if 516 #endif 517 end if 518 519 !Initialize band structure datatype 520 if (dtset%paral_kgb/=0) then ! We decide to store total npw in bstruct, 521 ABI_ALLOCATE(npwarr_,(dtset%nkpt)) 522 npwarr_(:)=npwarr(:) 523 call xmpi_sum(npwarr_,mpi_enreg%comm_bandfft,ierr) 524 else 525 npwarr_ => npwarr 526 end if 527 528 bstruct = ebands_from_dtset(dtset, npwarr_) 529 530 if (dtset%paral_kgb/=0) then 531 ABI_DEALLOCATE(npwarr_) 532 end if 533 nullify(npwarr_) 534 535 !Initialize PAW atomic occupancies 536 if (scf_history%history_size>=0) then 537 pawrhoij => scf_history%pawrhoij_last 538 else 539 ABI_DATATYPE_ALLOCATE(pawrhoij,(my_natom*psps%usepaw)) 540 end if 541 if (psps%usepaw==1.and.has_to_init) then 542 call initrhoij(dtset%pawcpxocc,dtset%lexexch,& 543 & dtset%lpawu,my_natom,dtset%natom,dtset%nspden,dtset%nspinor,dtset%nsppol,& 544 & dtset%ntypat,pawrhoij,dtset%pawspnorb,pawtab,dtset%spinat,dtset%typat,& 545 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 546 end if 547 548 !Initialize header 549 gscase=0 550 call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr,& 551 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 552 553 !Clean band structure datatype (should use it more in the future !) 554 call ebands_free(bstruct) 555 556 !Update header, with evolving variables, when available 557 !Here, rprimd, xred and occ are available 558 etot=hdr%etot ; fermie=hdr%fermie ; residm=hdr%residm 559 call hdr_update(hdr,bantot,etot,fermie,& 560 & residm,rprimd,occ,pawrhoij,xred,args_gs%amu,& 561 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 562 563 ! PW basis set: test if the problem is ill-defined. 564 if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then 565 npwmin=minval(hdr%npwarr(:)) 566 if (dtset%mband > npwmin) then 567 ! No way we can solve the problem. Abort now! 568 write(message,"(2(a,i0),4a)")& 569 & "Number of bands nband= ",dtset%mband," > number of planewaves npw= ",npwmin,ch10,& 570 & "The number of eigenvectors cannot be greater that the size of the Hamiltonian!",ch10,& 571 & "Action: decrease nband or, alternatively, increase ecut" 572 if (dtset%ionmov/=23) then 573 MSG_ERROR(message) 574 else 575 MSG_WARNING(message) 576 end if 577 578 else if (dtset%mband >= 0.9 * npwmin) then 579 ! Warn the user 580 write(message,"(a,i0,a,f6.1,4a)")& 581 & "Number of bands nband= ",dtset%mband," >= 0.9 * maximum number of planewaves= ",0.9*npwmin,ch10,& 582 & "The problem is ill-defined and the GS algorithm will show numerical instabilities!",ch10,& 583 & "Assume experienced user. Execution will continue." 584 MSG_WARNING(message) 585 end if 586 end if 587 588 !########################################################### 589 !### 04. Symmetry operations when nsym>1 590 591 !Do symmetry stuff only for nsym>1 592 if (dtset%usewvl == 0) then 593 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 594 else 595 #if defined HAVE_BIGDFT 596 nfftot=product(wvl%den%denspot%dpbox%ndims) 597 #else 598 BIGDFT_NOTENABLED_ERROR() 599 #endif 600 end if 601 ABI_ALLOCATE(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 602 ABI_ALLOCATE(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 603 ABI_ALLOCATE(indsym,(4,dtset%nsym,dtset%natom)) 604 ABI_ALLOCATE(symrec,(3,3,dtset%nsym)) 605 irrzon(:,:,:)=0 606 phnons(:,:,:)=zero 607 indsym(:,:,:)=0 608 symrec(:,:,:)=0 609 610 if (dtset%nsym>1) then 611 call setsym(indsym,irrzon,dtset%iscf,dtset%natom,& 612 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 613 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred) 614 615 ! Make sure dtset%iatfix does not break symmetry 616 call fixsym(dtset%iatfix,indsym,dtset%natom,dtset%nsym) 617 else 618 ! The symrec array is used by initberry even in case nsym = 1 619 symrec(:,:,1) = 0 620 symrec(1,1,1) = 1 ; symrec(2,2,1) = 1 ; symrec(3,3,1) = 1 621 end if 622 if (dtset%usewvl == 1) then 623 call wvl_descr_atoms_set_sym(wvl%descr, dtset%efield, irrzon, dtset%nsppol, & 624 & dtset%nsym, phnons, dtset%symafm, dtset%symrel, dtset%tnons, dtset%tolsym) 625 #if defined HAVE_BIGDFT 626 wvl%den%symObj = wvl%descr%atoms%astruct%sym%symObj 627 #endif 628 end if 629 630 !########################################################### 631 !### 05. Calls inwffil 632 633 ! if paral_kgb == 0, it may happen that some processors are idle (no entry in proc_distrb) 634 ! but mkmem == nkpt and this can cause integer overflow in mcg or allocation error. 635 ! Here we count the number of states treated by the proc. if cnt == 0, mcg is then set to 0. 636 cnt = 0 637 do spin=1,dtset%nsppol 638 do ikpt=1,dtset%nkpt 639 do band=1,dtset%nband(ikpt + (spin-1) * dtset%nkpt) 640 if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, band, band, spin, mpi_enreg%me_kpt)) cnt = cnt + 1 641 end do 642 end do 643 end do 644 645 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 646 mcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 647 if (cnt == 0) then 648 mcg = 0 649 write(message,"(2(a,i0))")"rank: ",mpi_enreg%me, "does not have wavefunctions to treat. Setting mcg to: ",mcg 650 MSG_WARNING(message) 651 end if 652 653 if (dtset%usewvl == 0 .and. dtset%mpw > 0 .and. cnt /= 0)then 654 if (my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol > floor(real(HUGE(0))/real(dtset%mpw) )) then 655 write (message,'(9a)')& 656 & "Default integer is not wide enough to store the size of the wavefunction array (mcg).",ch10,& 657 & "This usually happens when paral_kgb == 0 and there are not enough procs to distribute kpts and spins",ch10,& 658 & "Action: if paral_kgb == 0, use nprocs = nkpt * nsppol to reduce the memory per node.",ch10,& 659 & "If this does not solve the problem, use paral_kgb 1 with nprocs > nkpt * nsppol and use npfft/npband/npspinor",ch10,& 660 & "to decrease the memory requirements. Consider also OpenMP threads." 661 MSG_ERROR_NOSTOP(message,ii) 662 write (message,'(5(a,i0), 2a)')& 663 & "my_nspinor: ",my_nspinor, ", mpw: ",dtset%mpw, ", mband: ",dtset%mband,& 664 & ", mkmem: ",dtset%mkmem, ", nsppol: ",dtset%nsppol,ch10,& 665 & 'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS...) compiled in int64 mode' 666 MSG_ERROR(message) 667 end if 668 end if 669 670 ABI_STAT_ALLOCATE(cg,(2,mcg), ierr) 671 ABI_CHECK(ierr==0, "out of memory in cg") 672 673 ABI_ALLOCATE(eigen,(dtset%mband*dtset%nkpt*dtset%nsppol)) 674 ABI_ALLOCATE(resid,(dtset%mband*dtset%nkpt*dtset%nsppol)) 675 eigen(:)=zero ; resid(:)=zero 676 !mpi_enreg%paralbd=0 ; ask_accurate=0 677 ask_accurate=0 678 679 !WVL - Branching, allocating wavefunctions as wavelets. 680 if (dtset%usewvl == 1) then 681 call wvl_wfs_lr_copy(wvl%wfs, wvl%descr) 682 ! Create access arrays for wavefunctions and allocate wvl%wfs%psi (other arrays are left unallocated). 683 call wvl_wfs_set(dtset%strprecon,dtset%spinmagntarget, dtset%kpt, mpi_enreg%me_wvl,& 684 & dtset%natom, sum(dtset%nband), & 685 & dtset%nkpt, mpi_enreg%nproc_wvl, dtset%nspinor, dtset%nsppol, dtset%nwfshist, dtset%occ_orig, & 686 & psps, rprimd, wvl%wfs, dtset%wtk, wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult, & 687 & xred) 688 ! We transfer wavelets information to the hdr structure. 689 #if defined HAVE_BIGDFT 690 call local_potential_dimensions(mpi_enreg%me_wvl,wvl%wfs%ks%lzd,wvl%wfs%ks%orbs,wvl%den%denspot%xc,& 691 & wvl%den%denspot%dpbox%ngatherarr(0,1)) 692 hdr%nwvlarr(1) = wvl%wfs%ks%lzd%Glr%wfd%nvctr_c 693 hdr%nwvlarr(2) = 7 * wvl%wfs%ks%lzd%Glr%wfd%nvctr_f 694 #endif 695 ! Create access arrays for projectors and allocate them. 696 ! Compute projectors from each atom. 697 call wvl_projectors_set(mpi_enreg%me_wvl, dtset%natom, wvl%projectors, psps, rprimd, & 698 & wvl%wfs, wvl%descr, dtset%wvl_frmult, xred) 699 end if 700 701 read_wf_or_den=(dtset%iscf<=0.or.dtfil%ireadwf/=0.or.(dtfil%ireadden/=0.and.dtset%positron<=0)) 702 read_wf_or_den=(read_wf_or_den.and.has_to_init) 703 704 !RECURSION - initialization 705 if(has_to_init .and. dtset%userec==1) then 706 call InitRec(dtset,mpi_enreg,rec_set,rmet,maxval(psps%indlmn(3,:,:))) 707 end if 708 709 !LOTF - initialization 710 #if defined HAVE_LOTF 711 if(has_to_init .and. dtset%ionmov==23) then 712 call lotfparam_init(dtset%natom,dtset%lotf_version,1,& 713 & dtset%lotf_nitex,dtset%lotf_nneigx,& 714 & dtset%lotf_classic,1,1) 715 end if 716 #endif 717 718 !Initialize wavefunctions. 719 if(dtset%tfkinfunc /=2) then 720 wff1%unwff=dtfil%unwff1 721 optorth=1 !if (psps%usepaw==1) optorth=0 722 if(psps%usepaw==1 .and. dtfil%ireadwf==1)optorth=0 723 hdr%rprimd=rprimd_for_kg ! We need the rprimd that was used to generate de G vectors 724 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen,& 725 & dtset%exchn2n3d,formeig,hdr,dtfil%ireadwf,dtset%istwfk,kg,& 726 & dtset%kptns,dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,& 727 & dtset%mpw,dtset%nband,ngfft,dtset%nkpt,npwarr,& 728 & dtset%nsppol,dtset%nsym,occ,optorth,dtset%symafm,& 729 & dtset%symrel,dtset%tnons,dtfil%unkg,wff1,wffnow,dtfil%unwff1,& 730 & dtfil%fnamewffk,wvl) 731 hdr%rprimd=rprimd 732 end if 733 734 if (psps%usepaw==1.and.dtfil%ireadwf==1)then 735 call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 736 ! Has to update header again (because pawrhoij has changed) - MT 2007-10-22: Why ? 737 ! call hdr_update(hdr,bantot,etot,fermie,residm,rprimd,occ,pawrhoij,xred,args_gs%amu, & 738 ! & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 739 end if 740 741 !########################################################### 742 !### 06. Operations related to restartxf (Old version) 743 744 !TODO: Remove ab_xfh 745 !Initialize xf history (should be put in inwffil) 746 ab_xfh%nxfh=0 747 if(dtset%restartxf>=1 .and. dtfil%ireadwf==1)then 748 749 ! Should exchange the data about history in parallel localrdwf==0 750 if(xmpi_paral==1 .and. dtset%localrdwf==0)then 751 write(message, '(a,a,a)' )& 752 & 'It is not yet possible to use non-zero restartxf,',ch10,& 753 & 'in parallel, when localrdwf=0. Sorry for this ...' 754 MSG_BUG(message) 755 end if 756 757 ABI_ALLOCATE(ab_xfh%xfhist,(3,dtset%natom+4,2,0)) 758 call outxfhist(ab_xfh,dtset%natom,2,wff1,ios) 759 ABI_DEALLOCATE(ab_xfh%xfhist) 760 761 if(ios>0)then 762 write(message,'(a,a,a)')& 763 & 'An error occurred reading the input wavefunction file,',ch10,& 764 & 'with restartxf=1.' 765 MSG_ERROR(message) 766 else if(ios==0)then 767 write(message, '(a,a,i4,a)' )ch10,& 768 & ' gstate : reading',ab_xfh%nxfh,' (x,f) history pairs from input wf file.' 769 call wrtout(std_out,message,'COLL') 770 call wrtout(ab_out,message,'COLL') 771 end if 772 ! WARNING : should check that restartxf is not negative 773 ! WARNING : should check that restartxf /= only when dtfil%ireadwf is activated 774 end if 775 776 !Allocate the xf history array : takes into account the existing 777 !pairs, minus those that will be discarded, then those that will 778 !be computed, governed by dtset%ntime, and some additional pairs 779 !(needed when it will be possible to use xfhist for move.f) 780 ab_xfh%mxfh=(ab_xfh%nxfh-dtset%restartxf+1)+dtset%ntime+5 781 ABI_ALLOCATE(ab_xfh%xfhist,(3,dtset%natom+4,2,ab_xfh%mxfh)) 782 ab_xfh%xfhist(:,:,:,:) = zero 783 !WARNING : should check that the number of atoms in the wf file and natom are the same 784 785 !Initialize the xf history array 786 if(ab_xfh%nxfh>=dtset%restartxf .and. ab_xfh%nxfh>0)then 787 ! Eventually skip some of the previous history 788 if(dtset%restartxf>=2)then 789 do ixfh=1,dtset%restartxf-1 790 call WffReadSkipRec(ios,1,wff1) 791 end do 792 end if 793 794 ! Read and store the relevant history 795 ab_xfh%nxfhr=ab_xfh%nxfh-dtset%restartxf+1 796 call outxfhist(ab_xfh,dtset%natom,3,wff1,ios) 797 end if 798 799 !Close wff1, if it was ever opened (in inwffil) 800 if (dtfil%ireadwf==1) then 801 call WffClose(wff1,ierr) 802 end if 803 804 !########################################################### 805 !### 07. Calls setup2 806 807 !Further setup 808 ABI_ALLOCATE(start,(3,dtset%natom)) 809 call setup2(dtset,npwtot,start,wvl%wfs,xred) 810 811 !Allocation of previous atomic positions 812 if (scf_history%history_size>=0) then 813 xred_old => scf_history%xred_last 814 else 815 ABI_ALLOCATE(xred_old,(3,dtset%natom)) 816 end if 817 if (has_to_init) xred_old=xred 818 819 !Timing for initialisation period 820 call timab(33,2,tsec) 821 call timab(34,3,tsec) 822 823 !########################################################### 824 !### 08. Compute new occupation numbers 825 826 !Compute new occupation numbers, in case wavefunctions and eigenenergies 827 !were read from disk, occupation scheme is metallic (this excludes iscf=-1), 828 !and occupation numbers are required by iscf 829 if( dtfil%ireadwf==1 .and. & 830 & (dtset%occopt>=3.and.dtset%occopt<=8) .and. & 831 & (dtset%iscf>0 .or. dtset%iscf==-3) .and. dtset%positron/=1 ) then 832 833 ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 834 ! Warning : ideally, results_gs%entropy should not be set up here XG 20011007 835 ! Do not take into account the possible STM bias 836 call newocc(doccde,eigen,results_gs%energies%entropy,& 837 & results_gs%energies%e_fermie,& 838 & dtset%spinmagntarget,dtset%mband,dtset%nband,& 839 & dtset%nelect,dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,& 840 & dtset%occopt,dtset%prtvol,zero,dtset%tphysel,dtset%tsmear,dtset%wtk) 841 if (dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(args_gs%upawu(:))>=tol8.or. & 842 & sum(args_gs%jpawu(:))>tol8).and.dtset%dmft_entropy==0) results_gs%energies%entropy=zero 843 ABI_DEALLOCATE(doccde) 844 845 ! Transfer occupations to bigdft object: 846 if(dtset%usewvl==1 .and. .not. wvlbigdft) then 847 call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,1,wvl%wfs) 848 ! call wvl_energies_abi2big(results_gs%energies,wvl%wfs,2) 849 end if 850 851 else 852 ! Warning : ideally, results_gs%entropy should not be set up here XG 20011007 853 results_gs%energies%entropy=zero 854 end if 855 856 !########################################################### 857 !### 09. Generate an index table of atoms 858 859 !Definition of atindx array 860 !Generate an index table of atoms, in order for them to be used type after type. 861 ABI_ALLOCATE(atindx,(dtset%natom)) 862 ABI_ALLOCATE(atindx1,(dtset%natom)) 863 ABI_ALLOCATE(nattyp,(psps%ntypat)) 864 indx=1 865 do itypat=1,psps%ntypat 866 nattyp(itypat)=0 867 do iatom=1,dtset%natom 868 if(dtset%typat(iatom)==itypat)then 869 atindx(iatom)=indx 870 atindx1(indx)=iatom 871 indx=indx+1 872 nattyp(itypat)=nattyp(itypat)+1 873 end if 874 end do 875 end do 876 877 !Compute structure factor phases for current atomic pos: 878 if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then 879 ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*dtset%natom)) 880 call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred) 881 end if 882 883 !Here allocation of GPU for vtorho calculations 884 #if defined HAVE_GPU_CUDA 885 if (dtset%use_gpu_cuda==1) then 886 call alloc_hamilt_gpu(atindx1,dtset,gprimd,mpi_enreg,nattyp,npwarr,2,psps,dtset%use_gpu_cuda) 887 end if 888 #endif 889 890 !########################################################### 891 !### 10. PAW related operations 892 893 !Initialize paw_dmft, even if neither dmft not paw are used 894 !write(std_out,*) "dtset%usedmft",dtset%usedmft 895 use_sc_dmft=dtset%usedmft 896 if(dtset%paral_kgb>0) use_sc_dmft=0 897 call init_sc_dmft(dtset%nbandkss,dtset%dmftbandi,dtset%dmftbandf,dtset%dmft_read_occnd,dtset%mband,& 898 & dtset%nband,dtset%nkpt,dtset%nspden, & 899 & dtset%nspinor,dtset%nsppol,occ,dtset%usedmft,paw_dmft,use_sc_dmft,dtset%dmft_solv,mpi_enreg) 900 if (paw_dmft%use_dmft==1.and.me==0) then 901 call readocc_dmft(paw_dmft,dtfil%filnam_ds(3),dtfil%filnam_ds(4)) 902 end if 903 !Should be done inside init_sc_dmft 904 if ( dtset%usedmft /= 0 ) then 905 call data4entropyDMFT_init(paw_dmft%forentropyDMFT,& 906 dtset%natom,& 907 dtset%typat,& 908 dtset%lpawu,& 909 dtset%dmft_t2g==1, & 910 args_gs%upawu,& 911 args_gs%jpawu) 912 end if 913 !write(std_out,*) "paw_dmft%use_dmft",paw_dmft%use_dmft 914 915 !PAW: 1- Initialize values for several arrays unchanged during iterations 916 !2- Initialize data for LDA+U 917 !3- Eventually open temporary storage file 918 if(psps%usepaw==1) then 919 ! 1- 920 gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2 921 922 ! Test if we have to call pawinit 923 ! Some gen-cond have to be added... 924 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 925 926 if (psp_gencond==1.or.call_pawinit) then 927 call timab(553,1,tsec) 928 gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2 929 hyb_range_fock=zero;if (dtset%ixc<0) call libxc_functionals_get_hybridparams(hyb_range=hyb_range_fock) 930 call pawinit(gnt_option,gsqcut_shp,hyb_range_fock,dtset%pawlcutd,dtset%pawlmix,& 931 & psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,& 932 & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%xclevel,dtset%usepotzero) 933 934 ! Update internal values 935 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 936 call timab(553,2,tsec) 937 #if defined HAVE_BIGDFT 938 ! In the PAW+WVL case, copy sij: 939 if(dtset%usewvl==1) then 940 do itypat=1,dtset%ntypat 941 wvl%descr%paw%sij(1:pawtab(itypat)%lmn2_size,itypat)=pawtab(itypat)%sij(:) 942 end do 943 end if 944 #endif 945 end if 946 psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore) 947 call setsymrhoij(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot) 948 949 ! 2-Initialize and compute data for LDA+U, EXX, or LDA+DMFT 950 pawtab(:)%usepawu=0 951 pawtab(:)%useexexch=0 952 pawtab(:)%exchmix=zero 953 if(paw_dmft%use_dmft==1) then 954 call print_sc_dmft(paw_dmft,dtset%pawprtvol) 955 end if 956 if (dtset%usepawu>0.or.dtset%useexexch>0.or.paw_dmft%use_dmft>0) then 957 call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,& 958 & args_gs%jpawu,dtset%lexexch,dtset%lpawu,dtset%ntypat,pawang,dtset%pawprtvol,& 959 & pawrad,pawtab,args_gs%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu,ucrpa=dtset%ucrpa) 960 end if 961 end if 962 963 !########################################################### 964 !### 11. Initialize (eventually) electron-positron data and 965 !### electric and magnetic field data 966 967 !Initialize (eventually) electron-positron data 968 nullify (electronpositron) 969 if (dtset%positron/=0) then 970 call init_electronpositron(dtfil%ireadwf,dtset,electronpositron,mpi_enreg,nfftf,pawrhoij,pawtab) 971 end if 972 973 !########################################################### 974 !### 12. Operations dependent of iscf value 975 976 !Get starting charge density : rhor as well as rhog 977 !Also initialize the kinetic energy density 978 if (scf_history%history_size>=0) then 979 rhor => scf_history%rhor_last 980 taur => scf_history%taur_last 981 else 982 ABI_ALLOCATE(rhor,(nfftf,dtset%nspden)) 983 ABI_ALLOCATE(taur,(nfftf,dtset%nspden*dtset%usekden)) 984 end if 985 ABI_ALLOCATE(rhog,(2,nfftf)) 986 ABI_ALLOCATE(taug,(2,nfftf*dtset%usekden)) 987 988 if (has_to_init) then 989 if (dtset%iscf>0 .or. (dtset%iscf==0 .and. dtset%usewvl==1 )) then ! .and. dtset%usepaw==1)) then 990 991 if(dtfil%ireadden/=0.and.dtset%positron<=0)then 992 ! Read density 993 rdwrpaw=psps%usepaw; if(dtfil%ireadwf/=0) rdwrpaw=0 994 if (dtset%usewvl==0) then 995 call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, & 996 mpi_enreg, rhor, hdr_den, pawrhoij, comm, check_hdr=hdr, allow_interp=.True.) 997 results_gs%etotal = hdr_den%etot; call hdr_free(hdr_den) 998 else 999 fform=52 ; accessfil=0 1000 if (dtset%iomode == IO_MODE_MPI ) accessfil=4 1001 if (dtset%iomode == IO_MODE_ETSF) accessfil=3 1002 call ioarr(accessfil,rhor,dtset,results_gs%etotal,fform,dtfil%fildensin,hdr,& 1003 & mpi_enreg,ngfftf,cplex1,nfftf,pawrhoij,1,rdwrpaw,wvl%den) 1004 end if 1005 1006 if (rdwrpaw/=0) then 1007 call hdr_update(hdr,bantot,etot,fermie,residm,& 1008 & rprimd,occ,pawrhoij,xred,args_gs%amu,& 1009 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1010 end if 1011 1012 ! Read kinetic energy density 1013 if(dtfil%ireadkden/=0 .and. dtset%usekden==1 )then 1014 call read_rhor(dtfil%filkdensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, & 1015 mpi_enreg, taur, hdr_den, pawrhoij, comm, check_hdr=hdr, allow_interp=.True.) 1016 call hdr_free(hdr_den) 1017 end if 1018 1019 ! Compute up+down rho(G) by fft 1020 ABI_ALLOCATE(work,(nfftf)) 1021 work(:)=rhor(:,1) 1022 call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 1023 if(dtset%usekden==1)then 1024 work(:)=taur(:,1) 1025 call fourdp(1,taug,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 1026 end if 1027 ABI_DEALLOCATE(work) 1028 1029 else if(dtfil%ireadwf/=0)then 1030 izero=0 1031 ! Obtain the charge density from wfs that were read previously 1032 ! Be careful: in PAW, rho does not include the compensation 1033 ! density (to be added in scfcv.F90) ! 1034 ! tim_mkrho=1 ; mpi_enreg%paralbd=0 1035 tim_mkrho=1 1036 if (psps%usepaw==1) then 1037 ABI_ALLOCATE(rhowfg,(2,dtset%nfft)) 1038 ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden)) 1039 ! write(std_out,*) "mkrhogstate" 1040 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 1041 & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 1042 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor) 1043 ABI_DEALLOCATE(rhowfg) 1044 ABI_DEALLOCATE(rhowfr) 1045 else 1046 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 1047 & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 1048 if(dtset%usekden==1)then 1049 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 1050 & mpi_enreg,npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 1051 end if 1052 1053 end if 1054 1055 else if(dtfil%ireadwf==0.and.dtset%positron/=1)then 1056 1057 ! Crude, but realistic initialisation of the density 1058 ! There is not point to compute it from random wavefunctions. 1059 if (dtset%usewvl == 0) then 1060 call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,& 1061 & mgfftf,mpi_enreg,psps%mqgrid_vl,dtset%natom,nattyp,nfftf,& 1062 & ngfftf,dtset%nspden,psps%ntypat,dtset%paral_kgb,psps,pawtab,ph1df,& 1063 & psps%qgrid_vl,rhog,rhor,dtset%spinat,ucvol,psps%usepaw,& 1064 & dtset%ziontypat,dtset%znucl) 1065 ! Update initialized density taking into account jellium slab 1066 if(dtset%jellslab/=0) then 1067 option=2 1068 ABI_ALLOCATE(work,(nfftf)) 1069 call jellium(gmet,gsqcut_eff,mpi_enreg,nfftf,ngfftf,dtset%nspden,& 1070 & option,dtset%paral_kgb,dtset%slabwsrad,rhog,rhor,rprimd,work,dtset%slabzbeg,dtset%slabzend) 1071 ABI_DEALLOCATE(work) 1072 end if ! of usejell 1073 ! Kinetic energy density initialized to zero (used only in metaGGAs ... ) 1074 if(dtset%usekden==1)then 1075 taur=zero ; taug=zero 1076 end if 1077 else if (dtset%usewvl/=0) then 1078 ! skip for the moment for wvl+paw, not yet coded 1079 if(dtset%usepaw==0) then 1080 !Wavelet density corresponds exactly to the wavefunctions, 1081 !since wavefunctions are taken from diagonalisation of LCAO. 1082 call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl%wfs, wvl%den) 1083 else !usepaw 1084 #if defined HAVE_BIGDFT 1085 call wvl_initro(atindx1,wvl%descr%atoms%astruct%geocode,wvl%descr%h,mpi_enreg%me_wvl,& 1086 & dtset%natom,nattyp,nfftf,dtset%nspden,psps%ntypat,& 1087 & wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,& 1088 & wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,& 1089 & wvl%descr%Glr%d%n3,& 1090 & pawrad,pawtab,psps%gth_params%psppar,rhor,rprimd,& 1091 & dtset%spinat,wvl%den,dtset%xc_denpos,xred,dtset%ziontypat) 1092 call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl%wfs, wvl%den) 1093 #endif 1094 end if 1095 end if 1096 1097 end if 1098 1099 else if ((dtset%iscf==-1.or.dtset%iscf==-2.or.dtset%iscf==-3).and.dtset%positron<=0) then 1100 1101 call status(0,dtfil%filstat,iexit,level,'call ioarr ') 1102 ! Read rho(r) from a disk file 1103 rdwrpaw=psps%usepaw 1104 ! Note : results_gs%etotal is read here, 1105 ! and might serve in the tddft routine, but it is contrary to the intended use of results_gs ... 1106 ! Warning : should check the use of results_gs%e_fermie 1107 ! Warning : should check the use of results_gs%residm 1108 ! One might make them separate variables. 1109 1110 ! Read density and get Fermi level from hdr_den 1111 call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, & 1112 mpi_enreg, rhor, hdr_den, pawrhoij, comm, check_hdr=hdr) 1113 results_gs%etotal = hdr_den%etot; results_gs%energies%e_fermie = hdr_den%fermie 1114 call hdr_free(hdr_den) 1115 1116 if(dtfil%ireadkden/=0 .and. dtset%usekden==1)then 1117 call read_rhor(dtfil%filkdensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, & 1118 mpi_enreg, taur, hdr_den, pawrhoij, comm, check_hdr=hdr) 1119 call hdr_free(hdr_den) 1120 end if 1121 1122 ! Compute up+down rho(G) by fft 1123 ABI_ALLOCATE(work,(nfftf)) 1124 work(:)=rhor(:,1) 1125 call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 1126 if(dtset%usekden==1)then 1127 work(:)=taur(:,1) 1128 call fourdp(1,taug,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0) 1129 end if 1130 ABI_DEALLOCATE(work) 1131 1132 end if 1133 end if ! has_to_init 1134 1135 !########################################################### 1136 !### 13. If needed, initialize SCF history variables 1137 1138 !If needed, initialize atomic density in SCF history 1139 if (scf_history%history_size>0.and.has_to_init) then 1140 ! If rhor is an atomic density, just store it in history 1141 if (.not.read_wf_or_den) then 1142 scf_history%atmrho_last(:)=rhor(:,1) 1143 else 1144 ! If rhor is not an atomic density, has to compute rho_at(r) 1145 ABI_ALLOCATE(rhowfg,(2,nfftf)) 1146 ABI_ALLOCATE(rhowfr,(nfftf,1)) 1147 ABI_ALLOCATE(spinat_dum,(3,dtset%natom)) 1148 spinat_dum=zero 1149 call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,mgfftf,mpi_enreg,& 1150 & psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,1,psps%ntypat,dtset%paral_kgb,psps,pawtab,& 1151 & ph1df,psps%qgrid_vl,rhowfg,rhowfr,spinat_dum,ucvol,& 1152 & psps%usepaw,dtset%ziontypat,dtset%znucl) 1153 scf_history%atmrho_last(:)=rhowfr(:,1) 1154 ABI_DEALLOCATE(rhowfg) 1155 ABI_DEALLOCATE(rhowfr) 1156 ABI_DEALLOCATE(spinat_dum) 1157 end if 1158 end if 1159 1160 if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then 1161 ABI_DEALLOCATE(ph1df) 1162 end if 1163 1164 !!Electric field: initialization stage 1165 !!further initialization and updates happen in scfcv.F90 1166 call init_e_field_vars(dtefield,dtset,gmet,gprimd,kg,& 1167 & mpi_enreg,npwarr,occ,pawang,pawrad,pawtab,psps,& 1168 & pwind,pwind_alloc,pwnsfac,rprimd,symrec,xred) 1169 1170 !! orbital magnetization initialization 1171 dtorbmag%orbmag = dtset%orbmag 1172 if (dtorbmag%orbmag > 0) then 1173 call initorbmag(dtorbmag,dtset,gmet,gprimd,kg,mpi_enreg,npwarr,occ,& 1174 & pawtab,psps,pwind,pwind_alloc,pwnsfac,& 1175 & rprimd,symrec,xred) 1176 end if 1177 1178 1179 fatvshift=one 1180 1181 if (dtset%usewvl == 1 .and. wvl_debug) then 1182 #if defined HAVE_BIGDFT 1183 call wvl_timing(me,'INIT','PR') 1184 #endif 1185 end if 1186 1187 !Check whether exiting was required by the user. If found then do not start minimization steps 1188 !At this first call to chkexi, initialize cpus, if it 1189 !is non-zero (which would mean that no action has to be taken) 1190 !Should do this in driver ... 1191 cpus=dtset%cpus 1192 if(abs(cpus)>1.0d-5)cpus=cpus+cpui 1193 openexit=1 ; if(dtset%chkexit==0) openexit=0 1194 call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit) 1195 1196 !If immediate exit, and wavefunctions were not read, must zero eigenvalues 1197 if (iexit/=0) eigen(:)=zero 1198 1199 call timab(34,2,tsec) 1200 1201 conv_retcode = 0 1202 1203 if (iexit==0) then 1204 1205 ! ########################################################### 1206 ! ### 14. Move atoms and acell acording to ionmov value 1207 1208 1209 ! Eventually symmetrize atomic coordinates over space group elements: 1210 ! call symmetrize_xred(indsym,dtset%natom,dtset%nsym,dtset%symrel,dtset%tnons,xred) 1211 1212 ! If atoms are not being moved and U should not be determined, 1213 ! use scfcv directly; else 1214 ! call move, pawuj_drive or brdmin which in turn calls scfcv. 1215 1216 call timab(35,3,tsec) 1217 1218 call scfcv_init(scfcv_args,atindx,atindx1,cg,cpus,& 1219 & args_gs%dmatpawu,dtefield,dtfil,dtorbmag,dtpawuj,dtset,ecore,eigen,hdr,& 1220 & indsym,initialized,irrzon,kg,mcg,mpi_enreg,my_natom,nattyp,ndtpawuj,& 1221 & nfftf,npwarr,occ,pawang,pawfgr,pawrad,pawrhoij,& 1222 & pawtab,phnons,psps,pwind,pwind_alloc,pwnsfac,rec_set,& 1223 & resid,results_gs,scf_history,fatvshift,& 1224 & symrec,taug,taur,wvl,ylm,ylmgr,paw_dmft,wffnew,wffnow) 1225 1226 call dtfil_init_time(dtfil,0) 1227 1228 write(message,'(a,80a)')ch10,('=',mu=1,80) 1229 call wrtout(ab_out,message,'COLL') 1230 call wrtout(std_out,message,'COLL') 1231 1232 if (dtset%ionmov==0) then 1233 1234 ! Should merge this call with the call for dtset%ionmov==4 and 5 1235 1236 if (dtset%macro_uj==0) then 1237 1238 call status(0,dtfil%filstat,iexit,level,'call scfcv_run') 1239 call scfcv_run(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode) 1240 1241 else 1242 ! Conduct determination of U 1243 1244 call pawuj_drive(scfcv_args,dtset,electronpositron,rhog,rhor,rprimd,xred,xred_old) 1245 end if 1246 1247 ! ======================================== 1248 ! New structure for geometry optimization 1249 ! ======================================== 1250 else if (dtset%ionmov>50.or.dtset%ionmov<=25) then 1251 1252 ! TODO: return conv_retcode 1253 call mover(scfcv_args,ab_xfh,acell,amass,dtfil,& 1254 & electronpositron,rhog,rhor,rprimd,vel,vel_cell,xred,xred_old) 1255 1256 ! Compute rprim from rprimd and acell 1257 do kk=1,3 1258 do jj=1,3 1259 rprim(jj,kk)=rprimd(jj,kk)/acell(kk) 1260 end do 1261 end do 1262 1263 ! ========================================= 1264 ! New structure for geometry optimization 1265 ! ========================================= 1266 1267 else ! Not an allowed option 1268 write(message, '(a,i0,2a)' )& 1269 & 'Disallowed value for ionmov=',dtset%ionmov,ch10,& 1270 & 'Allowed values are: 1,2,3,4,5,6,7,8,9,10,11,12,13,14,20,21,22,23,24 and 30' 1271 MSG_BUG(message) 1272 end if 1273 1274 call scfcv_destroy(scfcv_args) 1275 1276 call timab(35,2,tsec) 1277 1278 ! ########################################################### 1279 ! ### 15. Final operations and output for gstate 1280 1281 end if ! End of the check of hasty exit 1282 1283 call timab(36,3,tsec) 1284 1285 write(message, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,& 1286 & ' ----iterations are completed or convergence reached----',ch10 1287 call wrtout(ab_out,message,'COLL') 1288 call wrtout(std_out, message,'COLL') 1289 1290 !Mark this GS computation as done 1291 initialized=1 1292 if (dtset%usewvl == 1 .and. wvl_debug) then 1293 #if defined HAVE_BIGDFT 1294 call wvl_timing(me,'WFN_OPT','PR') 1295 #endif 1296 end if 1297 1298 !Will be put here later. 1299 !! ! WVL - maybe compute the tail corrections to energy 1300 !! if (dtset%tl_radius > zero) then 1301 !! ! Store xcart for each atom 1302 !! allocate(xcart(3, dtset%natom)) 1303 !! call xred2xcart(dtset%natom, rprimd, xcart, xred) 1304 !! ! Use the tails to improve energy precision. 1305 !! call wvl_tail_corrections(dtset, results_gs%energies, results_gs%etotal, & 1306 !! & mpi_enreg, occ, psps, vtrial, wvl, xcart) 1307 !! deallocate(xcart) 1308 !! end if 1309 1310 !Update the header, before using it 1311 call hdr_update(hdr,bantot,results_gs%etotal,results_gs%energies%e_fermie,& 1312 & results_gs%residm,rprimd,occ,pawrhoij,xred,args_gs%amu,& 1313 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 1314 1315 ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1316 doccde=zero 1317 1318 call ebands_init(bantot,ebands,dtset%nelect,doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,& 1319 & hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,& 1320 & hdr%charge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, & 1321 & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk) 1322 1323 ebands%fermie = results_gs%energies%e_fermie 1324 ABI_DEALLOCATE(doccde) 1325 !write(std_out,*)"efermi after ebands_init",ebands%fermie 1326 1327 ! Compute and print the gaps. 1328 call ebands_report_gap(ebands,header="Gap info",unit=std_out,mode_paral="COLL",gaps=results_gs%gaps) 1329 1330 if(dtset%nqpt==0)filnam=dtfil%fnameabo_wfk 1331 if(dtset%nqpt==1)filnam=dtfil%fnameabo_wfq 1332 1333 ! Write wavefunctions file only if convergence was not achieved. 1334 !write(std_out,*)"conv_retcode", conv_retcode 1335 write_wfk = .True. 1336 if (dtset%prtwf==-1 .and. conv_retcode == 0) then 1337 write_wfk = .False. 1338 message = "GS calculation converged with prtwf=-1 --> Skipping WFK file output" 1339 call wrtout(ab_out, message, "COLL") 1340 MSG_COMMENT(message) 1341 end if 1342 1343 !To print out the WFs, need the rprimd that was used to generate the G vectors 1344 hdr%rprimd=rprimd_for_kg 1345 1346 if (write_wfk) then 1347 call outwf(cg,dtset,psps,eigen,filnam,hdr,kg,dtset%kptns,& 1348 & dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%natom,& 1349 & dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,& 1350 & occ,resid,response,dtfil%unwff2,wvl%wfs,wvl%descr) 1351 !SPr: add input variable managing the .vtk file OUTPUT (Please don't remove the next commented line) 1352 !call printmagvtk(mpi_enreg,cplex1,dtset%nspden,nfftf,ngfftf,rhor,rprimd,'DEN') 1353 end if 1354 1355 if (dtset%prtwf==2) then 1356 call outqmc(cg,dtset,eigen,gprimd,hdr,kg,mcg,mpi_enreg,npwarr,occ,psps,results_gs) 1357 end if 1358 1359 !Restore the original rprimd in hdr 1360 hdr%rprimd=rprimd 1361 1362 ! Generate WFK in full BZ (needed by LOBSTER) 1363 if (me == master .and. dtset%prtwf == 1 .and. dtset%prtwf_full == 1 .and. dtset%nqpt == 0) then 1364 wfkfull_path = strcat(dtfil%filnam_ds(4), "_FULL_WFK") 1365 if (dtset%iomode == IO_MODE_ETSF) wfkfull_path = nctk_ncify(wfkfull_path) 1366 call wfk_tofullbz(filnam, dtset, psps, pawtab, wfkfull_path) 1367 1368 ! Write tetrahedron tables. 1369 call crystal_from_hdr(cryst, hdr, 2) 1370 tetra = tetra_from_kptrlatt(cryst, dtset%kptopt, dtset%kptrlatt, dtset%nshiftk, & 1371 dtset%shiftk, dtset%nkpt, dtset%kptns, message, ierr) 1372 if (ierr == 0) then 1373 call tetra_write(tetra, dtset%nkpt, dtset%kptns, strcat(dtfil%filnam_ds(4), "_TETRA")) 1374 else 1375 MSG_WARNING(sjoin("Cannot produce TETRA file", ch10, message)) 1376 end if 1377 1378 call destroy_tetra(tetra) 1379 call crystal_free(cryst) 1380 end if 1381 1382 call clnup1(acell,dtset,eigen,results_gs%energies%e_fermie,& 1383 & dtfil%fnameabo_dos,dtfil%fnameabo_eig,results_gs%fred,& 1384 & mpi_enreg,nfftf,ngfftf,occ,dtset%optforces,& 1385 & resid,rhor,rprimd,results_gs%vxcavg,xred) 1386 1387 if ( (dtset%iscf>=0 .or. dtset%iscf==-3) .and. dtset%prtstm==0) then 1388 call prtene(dtset,results_gs%energies,ab_out,psps%usepaw) 1389 end if 1390 1391 !write final electric field components HONG 1392 1393 if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt ==7 .or. & 1394 & dtset%berryopt == 14 .or. dtset%berryopt == 16 .or. dtset%berryopt ==17 ) then ! output final elctric field data !!HONG 1395 if (dtset%berryopt == 4) then 1396 write(message,'(a,a)') ch10, 'Constant unreduced E calculation - final values:' 1397 else if (dtset%berryopt == 6 ) then 1398 write(message,'(a,a)') ch10, 'Constant unreduced D calculation - final values:' 1399 else if (dtset%berryopt == 14) then 1400 write(message,'(a,a)') ch10, 'Constant reduced ebar calculation - final values:' 1401 else if (dtset%berryopt == 16 ) then 1402 write(message,'(a,a)') ch10, 'Constant reduced d calculation - final values:' 1403 else if (dtset%berryopt == 17) then 1404 write(message,'(a,a)') ch10, 'Constant reduced ebar and d calculation - final values:' 1405 end if 1406 1407 call wrtout(ab_out,message,'COLL') 1408 call prtefield(dtset,dtefield,ab_out,rprimd) 1409 1410 call wrtout(std_out,message,'COLL') 1411 call prtefield(dtset,dtefield,std_out,rprimd) 1412 1413 ! To check if the final electric field is below the critical field 1414 do kk = 1, 3 1415 efield_band(kk) = abs(dtset%red_efieldbar(kk))*dtefield%nkstr(kk) 1416 end do 1417 ! eg = maxval(eg_dir) 1418 ! eg_ev = eg*Ha_eV 1419 write(message,'(a,a,a,a,a,a,a,a,f7.2,a,a)')ch10,& 1420 & ' Please check: COMMENT - ',ch10,& 1421 & ' As a rough estimate,',ch10,& 1422 & ' to be below the critical field, the bandgap of your system',ch10,& 1423 & ' should be larger than ',maxval(efield_band)*Ha_eV,' eV.',ch10 1424 call wrtout(ab_out,message,'COLL') 1425 call wrtout(std_out,message,'COLL') 1426 1427 write(message,'(a)') '--------------------------------------------------------------------------------' 1428 call wrtout(ab_out,message,'COLL') 1429 call wrtout(std_out,message,'COLL') 1430 1431 end if 1432 1433 !Open the formatted derivative database file, and write the preliminary information 1434 !In the // case, only one processor writes the energy and the gradients to the DDB 1435 1436 if (me==0.and.dtset%nimage==1.and.((dtset%iscf > 0).or.& 1437 & (dtset%berryopt == -1).or.(dtset%berryopt) == -3)) then 1438 1439 if (dtset%iscf > 0) then 1440 nblok = 2 ! 1st blok = energy, 2nd blok = gradients 1441 else 1442 nblok = 1 1443 end if 1444 1445 dscrpt=' Note : temporary (transfer) database ' 1446 ddbnm=trim(dtfil%filnam_ds(4))//'_DDB' 1447 1448 call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,& 1449 & nblok,xred=xred,occ=occ,ngfft=ngfft) 1450 1451 call ddb_hdr_open_write(ddb_hdr,ddbnm,dtfil%unddb,fullinit=0) 1452 1453 call ddb_hdr_free(ddb_hdr) 1454 1455 choice=2 1456 mpert = dtset%natom + 6 ; msize = 3*mpert 1457 1458 ! create a ddb structure with just one blok 1459 call ddb_malloc(ddb,msize,1,dtset%natom,dtset%ntypat) 1460 1461 ddb%flg = 0 1462 ddb%qpt = zero 1463 ddb%nrm = one 1464 ddb%val = zero 1465 1466 ! Write total energy to the DDB 1467 if (dtset%iscf > 0) then 1468 ddb%typ(1) = 0 1469 ddb%val(1,1,1) = results_gs%etotal 1470 ddb%flg(1,1) = 1 1471 call ddb_write_blok(ddb,1,choice,dtset%mband,mpert,msize,dtset%nkpt,dtfil%unddb) 1472 end if 1473 1474 ! Write gradients to the DDB 1475 ddb%typ = 4 1476 ddb%flg = 0 1477 ddb%val = zero 1478 indx = 0 1479 if (dtset%iscf > 0) then 1480 do iatom = 1, dtset%natom 1481 do idir = 1, 3 1482 indx = indx + 1 1483 ddb%flg(indx,1) = 1 1484 ddb%val(1,indx,1) = results_gs%fred(idir,iatom) 1485 end do 1486 end do 1487 end if 1488 1489 indx = 3*dtset%natom + 3 1490 if ((abs(dtset%berryopt) == 1).or.(abs(dtset%berryopt) == 3)) then 1491 do idir = 1, 3 1492 indx = indx + 1 1493 if (dtset%rfdir(idir) == 1) then 1494 ddb%flg(indx,1) = 1 1495 ddb%val(1,indx,1) = results_gs%pel(idir) 1496 end if 1497 end do 1498 end if 1499 1500 indx = 3*dtset%natom + 6 1501 if (dtset%iscf > 0) then 1502 ddb%flg(indx+1:indx+6,1) = 1 1503 ddb%val(1,indx+1:indx+6,1) = results_gs%strten(1:6) 1504 end if 1505 1506 call ddb_write_blok(ddb,1,choice,dtset%mband,mpert,msize,dtset%nkpt,dtfil%unddb) 1507 1508 call ddb_free(ddb) 1509 1510 ! Close DDB 1511 close(dtfil%unddb) 1512 end if 1513 1514 if (dtset%nstep>0 .and. dtset%prtstm==0 .and. dtset%positron/=1) then 1515 call clnup2(psps%n1xccc,results_gs%fred,results_gs%grchempottn,results_gs%gresid,& 1516 & results_gs%grewtn,results_gs%grvdw,results_gs%grxc,dtset%iscf,dtset%natom,& 1517 & results_gs%ngrvdw,dtset%optforces,dtset%optstress,dtset%prtvol,start,& 1518 & results_gs%strten,results_gs%synlgr,xred) 1519 end if 1520 1521 !Deallocate arrays 1522 ABI_DEALLOCATE(amass) 1523 ABI_DEALLOCATE(atindx) 1524 ABI_DEALLOCATE(atindx1) 1525 ABI_DEALLOCATE(cg) 1526 ABI_DEALLOCATE(eigen) 1527 ABI_DEALLOCATE(indsym) 1528 ABI_DEALLOCATE(npwarr) 1529 ABI_DEALLOCATE(nattyp) 1530 ABI_DEALLOCATE(resid) 1531 ABI_DEALLOCATE(rhog) 1532 ABI_DEALLOCATE(start) 1533 ABI_DEALLOCATE(symrec) 1534 ABI_DEALLOCATE(taug) 1535 ABI_DEALLOCATE(ab_xfh%xfhist) 1536 call pawfgr_destroy(pawfgr) 1537 1538 if (dtset%usewvl == 0 .or. dtset%nsym <= 1) then 1539 ! In wavelet case, irrzon and phnons are deallocated by wavelet object. 1540 ABI_DEALLOCATE(irrzon) 1541 ABI_DEALLOCATE(phnons) 1542 end if 1543 1544 ABI_DEALLOCATE(ylm) 1545 ABI_DEALLOCATE(ylmgr) 1546 1547 if (scf_history%history_size<0) then 1548 if (psps%usepaw==1) then 1549 call pawrhoij_free(pawrhoij) 1550 end if 1551 ABI_DEALLOCATE(rhor) 1552 ABI_DEALLOCATE(taur) 1553 ABI_DATATYPE_DEALLOCATE(pawrhoij) 1554 ABI_DEALLOCATE(xred_old) 1555 else 1556 nullify(rhor,taur,pawrhoij,xred_old) 1557 end if 1558 1559 !PAW+DMFT 1560 call destroy_sc_dmft(paw_dmft) 1561 ! This call should be done inside destroy_sc_dmft 1562 if ( dtset%usedmft /= 0 ) then 1563 call data4entropyDMFT_destroy(paw_dmft%forentropyDMFT) 1564 end if 1565 1566 !Destroy electronpositron datastructure 1567 if (dtset%positron/=0) then 1568 call destroy_electronpositron(electronpositron) 1569 end if 1570 1571 !Deallocating the basis set. 1572 if (dtset%usewvl == 1) then 1573 call wvl_projectors_free(wvl%projectors) 1574 call wvl_wfs_free(wvl%wfs) 1575 call wvl_descr_free(wvl%descr) 1576 call wvl_denspot_free(wvl%den) 1577 if(dtset%usepaw == 1) then 1578 call wvl_paw_free(wvl%descr) 1579 end if 1580 end if 1581 1582 ABI_DEALLOCATE(kg) 1583 1584 if (dtset%icoulomb /= 0) then 1585 call psolver_kernel((/ 0._dp, 0._dp, 0._dp /), 0, dtset%icoulomb, 0, kernel_dummy, & 1586 & 0, dtset%ngfft, 1, dtset%nscforder) 1587 end if 1588 1589 if (associated(pwind)) then 1590 ABI_DEALLOCATE(pwind) 1591 end if 1592 if (associated(pwnsfac)) then 1593 ABI_DEALLOCATE(pwnsfac) 1594 end if 1595 if ((dtset%berryopt<0).or.& 1596 & (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1597 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17)) then 1598 if (xmpi_paral == 1) then 1599 ABI_DEALLOCATE(mpi_enreg%kptdstrb) 1600 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1601 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 1602 ABI_DEALLOCATE(mpi_enreg%kpt_loc2ibz_sp) 1603 end if 1604 end if 1605 if (allocated(mpi_enreg%kpt_loc2ibz_sp)) then 1606 ABI_DEALLOCATE(mpi_enreg%kpt_loc2ibz_sp) 1607 end if 1608 if (allocated(mpi_enreg%kpt_loc2fbz_sp)) then 1609 ABI_DEALLOCATE(mpi_enreg%kpt_loc2fbz_sp) 1610 end if 1611 if (allocated(mpi_enreg%mkmem)) then 1612 ABI_DEALLOCATE(mpi_enreg%mkmem) 1613 end if 1614 end if 1615 1616 ! deallocate efield 1617 call destroy_efield(dtefield) 1618 1619 ! deallocate orbmag 1620 call destroy_orbmag(dtorbmag) 1621 1622 !deallocate Recursion 1623 if (dtset%userec == 1) then 1624 call CleanRec(rec_set) 1625 end if 1626 1627 !Clean the header 1628 call hdr_free(hdr) 1629 call ebands_free(ebands) 1630 1631 if (me == master .and. dtset%prtxml == 1) then 1632 ! The dataset given in argument has been treated, then we output its variables. 1633 ! call outvarsXML() 1634 ! gstate() will handle a dataset, so we output the dataSet markup. 1635 write(ab_xml_out, "(A)") ' </dataSet>' 1636 end if 1637 1638 if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2 .and. dtset%optdriver /= RUNL_GWLS) then 1639 ! Plane-wave case 1640 call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg) 1641 end if 1642 1643 if(dtset%wfoptalg == 1 .and. psps%usepaw == 1) then 1644 call destroy_invovl(dtset%nkpt) 1645 end if 1646 1647 if(gemm_nonlop_use_gemm) then 1648 call destroy_gemm_nonlop(dtset%nkpt) 1649 gemm_nonlop_use_gemm = .false. 1650 end if 1651 1652 !Eventually clean cuda runtime 1653 #if defined HAVE_GPU_CUDA 1654 if (dtset%use_gpu_cuda==1) then 1655 call dealloc_hamilt_gpu(2,dtset%use_gpu_cuda) 1656 end if 1657 #endif 1658 1659 call status(0,dtfil%filstat,iexit,level,'exit') 1660 1661 call timab(36,2,tsec) 1662 call timab(32,2,tsec) 1663 1664 DBG_EXIT("COLL")
ABINIT/setup2 [ 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
PARENTS
gstate
CHILDREN
wrtout
SOURCE
1699 subroutine setup2(dtset,npwtot,start,wfs,xred) 1700 1701 1702 !This section has been created automatically by the script Abilint (TD). 1703 !Do not modify the following lines by hand. 1704 #undef ABI_FUNC 1705 #define ABI_FUNC 'setup2' 1706 use interfaces_14_hidewrite 1707 !End of the abilint section 1708 1709 implicit none 1710 1711 !Arguments ------------------------------------ 1712 !scalars 1713 type(dataset_type),intent(in) :: dtset 1714 type(wvl_wf_type),intent(in) :: wfs 1715 !arrays 1716 integer,intent(in) :: npwtot(dtset%nkpt) 1717 real(dp),intent(in) :: xred(3,dtset%natom) 1718 real(dp),intent(out) :: start(3,dtset%natom) 1719 1720 !Local variables------------------------------- 1721 !scalars 1722 integer :: ikpt,npw 1723 real(dp) :: arith,geom,wtknrm 1724 character(len=500) :: message 1725 1726 ! ************************************************************************* 1727 1728 !DEBUG 1729 !write(std_out,*)' setup2 : enter ' 1730 !ENDDEBUG 1731 1732 if (dtset%iscf>=0) then 1733 1734 ! Copy coordinates into array start 1735 start(:,:)=xred(:,:) 1736 1737 if (dtset%usewvl == 0) then 1738 ! Get average number of planewaves per k point: 1739 ! both arithmetic and GEOMETRIC averages are desired-- 1740 ! need geometric average to use method of Francis and Payne, 1741 ! J. Phys.: Condens. Matter 2, 4395-4404 (1990). 1742 ! Also note: force k point wts to sum to 1 for this averaging. 1743 ! (wtk is not forced to add to 1 in a case with occopt=2) 1744 arith=zero 1745 geom=one 1746 wtknrm=zero 1747 do ikpt=1,dtset%nkpt 1748 npw=npwtot(ikpt) 1749 wtknrm=wtknrm+dtset%wtk(ikpt) 1750 arith=arith+npw*dtset%wtk(ikpt) 1751 geom=geom*npw**dtset%wtk(ikpt) 1752 end do 1753 1754 ! Enforce normalization of weights to 1 1755 arith=arith/wtknrm 1756 geom=geom**(1.0_dp/wtknrm) 1757 1758 end if 1759 1760 ! Ensure portability of output thanks to tol8 1761 if (dtset%usewvl == 0) then 1762 write(message, '(a,2f12.3)' ) & 1763 & '_setup2: Arith. and geom. avg. npw (full set) are',arith+tol8,geom 1764 else 1765 #if defined HAVE_BIGDFT 1766 write(message, '(a,2I8)' ) ' setup2: nwvl coarse and fine are', & 1767 & wfs%ks%lzd%Glr%wfd%nvctr_c, wfs%ks%lzd%Glr%wfd%nvctr_f 1768 #endif 1769 end if 1770 call wrtout(ab_out, message,'COLL') 1771 call wrtout(std_out, message,'COLL') 1772 1773 end if 1774 1775 #if !defined HAVE_BIGDFT 1776 if (.false.) write(std_out,*) wfs%ks 1777 #endif 1778 1779 end subroutine setup2