TABLE OF CONTENTS
- ABINIT/crystal_from_file
- ABINIT/ebands_from_file
- ABINIT/get_dtsets_pspheads
- ABINIT/m_common
- ABINIT/prteigrs
- ABINIT/prtene
- ABINIT/scprqt
- ABINIT/setup1
ABINIT/crystal_from_file [ Functions ]
NAME
crystal_from_file
FUNCTION
Build crystal_t object from netcdf file
INPUTS
OUTPUT
SOURCE
2017 type(crystal_t) function crystal_from_file(path, comm) result(new) 2018 2019 !Arguments ------------------------------------ 2020 !scalars 2021 character(len=*),intent(in) :: path 2022 integer,intent(in) :: comm 2023 2024 !Local variables------------------------------- 2025 !scalars 2026 integer :: fform 2027 type(hdr_type) :: hdr 2028 2029 ! ************************************************************************* 2030 2031 ! Assume file with Abinit header 2032 ! TODO: Should add routine to read crystal from structure without hdr 2033 call hdr_read_from_fname(hdr, path, fform, comm) 2034 ABI_CHECK(fform /= 0, "fform == 0") 2035 new = hdr%get_crystal() 2036 call hdr%free() 2037 2038 end function crystal_from_file
ABINIT/ebands_from_file [ Functions ]
NAME
ebands_from_file
FUNCTION
Build and ebands_t object from file. Supports Fortran and netcdf files provided they have a Abinit header and obviously GS eigenvalues
INPUTS
path: File name. comm: MPI communicator.
OUTPUT
SOURCE
1963 type(ebands_t) function ebands_from_file(path, comm) result(new) 1964 1965 !Arguments ------------------------------------ 1966 !scalars 1967 character(len=*),intent(in) :: path 1968 integer,intent(in) :: comm 1969 1970 !Local variables------------------------------- 1971 !scalars 1972 integer :: ncid, fform 1973 type(hdr_type) :: hdr 1974 !arrays 1975 real(dp),pointer :: gs_eigen(:,:,:) 1976 1977 ! ************************************************************************* 1978 1979 ! NOTE: Assume file with header. Must use wfk_read_eigenvalues to handle Fortran WFK 1980 if (endswith(path, "_WFK") .or. endswith(path, "_WFK.nc")) then 1981 call wfk_read_eigenvalues(path, gs_eigen, hdr, comm) 1982 new = ebands_from_hdr(hdr, maxval(hdr%nband), gs_eigen) 1983 1984 else if (endswith(path, ".nc")) then 1985 #ifdef HAVE_NETCDF 1986 NCF_CHECK(nctk_open_read(ncid, path, comm)) 1987 call hdr_ncread(hdr, ncid, fform) 1988 ABI_CHECK(fform /= 0, "fform == 0") 1989 ABI_MALLOC(gs_eigen, (hdr%mband, hdr%nkpt, hdr%nsppol)) 1990 NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "eigenvalues"), gs_eigen)) 1991 new = ebands_from_hdr(hdr, maxval(hdr%nband), gs_eigen) 1992 NCF_CHECK(nf90_close(ncid)) 1993 #endif 1994 else 1995 ABI_ERROR(sjoin("Don't know how to construct crystal structure from: ", path, ch10, "Supported extensions: _WFK or .nc")) 1996 end if 1997 1998 ABI_FREE(gs_eigen) 1999 call hdr%free() 2000 2001 end function ebands_from_file
ABINIT/get_dtsets_pspheads [ Functions ]
NAME
get_dtsets_pspheads
FUNCTION
Parse input file, get list of pseudos for files file and build list of datasets pseudopotential headers, maxval of dimensions needed in outvars
INPUTS
input_path: Input filename specifed on the command line. zero lenght if files file syntax is used.
Mainly used to check whether pseudos are defined in the input to avoid entering the files file
branch that prompts for pseudos.
path: Input Filename
comm: MPI communicator
OUTPUT
lenstr= the length of the resulting string. ndtset= the number of declared datasets. string= contains on output the content of the file, ready for parsing. dtsets(0:ndtset): List of datasets dmatpuflag=flag controlling the use of an initial density matrix in PAW+U (max. value over datasets) mx<ab_dimensions>=datatype storing the maximal dimensions. pspheads(npsp)=<type pspheader_type>=all the important information from the pseudopotential file headers, as well as the psp file names
SOURCE
1769 subroutine get_dtsets_pspheads(input_path, path, ndtset, lenstr, string, timopt, dtsets, pspheads, mx, dmatpuflag, comm) 1770 1771 !Arguments ------------------------------------ 1772 !scalars 1773 integer,intent(out) :: lenstr, ndtset 1774 type(ab_dimensions),intent(out) :: mx 1775 character(len=strlen), intent(out) :: string 1776 character(len=*),intent(in) :: input_path, path 1777 integer,intent(in) :: comm 1778 integer,intent(out) :: timopt, dmatpuflag 1779 !arrays 1780 type(dataset_type),allocatable,intent(out) :: dtsets(:) 1781 type(pspheader_type),allocatable,intent(out):: pspheads(:) 1782 1783 !Local variables------------------------------- 1784 !scalars 1785 integer :: ipsp,ios, me, ndtset_alloc, nprocs 1786 integer :: istatr,istatshft, papiopt, npsp, ii, idtset, msym, usepaw 1787 character(len=fnlen) :: filpsp 1788 character(len=500) :: msg 1789 !arrays 1790 integer,allocatable :: mband_upper_(:) 1791 real(dp) :: ecut_tmp(3,2,10),tsec(2) 1792 real(dp),allocatable :: zionpsp(:) 1793 character(len=fnlen), allocatable :: pspfilnam_(:), pseudo_paths(:) 1794 1795 !************************************************************************ 1796 1797 me = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1798 1799 ! Read the file, stringify it and return the number of datasets. 1800 call parsefile(path, lenstr, ndtset, string, comm) 1801 1802 ndtset_alloc = ndtset; if (ndtset == 0) ndtset_alloc=1 1803 ABI_MALLOC(dtsets, (0:ndtset_alloc)) 1804 1805 timopt = 1; if (xmpi_paral==1) timopt = 0 1806 1807 ! Continue to analyze the input string, get upper dimensions, and allocate the remaining arrays. 1808 call invars0(dtsets, istatr, istatshft, lenstr, msym, mx%natom, mx%nimage, mx%ntypat, & 1809 ndtset, ndtset_alloc, npsp, pseudo_paths, papiopt, timopt, string, comm) 1810 1811 ! Enable PAPI timers 1812 call time_set_papiopt(papiopt) 1813 1814 dtsets(:)%timopt = timopt 1815 dtsets(0)%timopt = 1 1816 if (xmpi_paral == 1) dtsets(0)%timopt = 0 1817 1818 call timab(timopt,5,tsec) 1819 1820 ! Initialize pspheads, that contains the important information 1821 ! from the pseudopotential headers, as well as the psp filename 1822 call timab(102,1,tsec) 1823 call timab(1021,3,tsec) 1824 1825 usepaw = 0 1826 ABI_MALLOC(pspheads, (npsp)) 1827 if (npsp > 10) then 1828 ABI_BUG('ecut_tmp is not well defined.') 1829 end if 1830 ecut_tmp = -one 1831 1832 pspheads(:)%usewvl = dtsets(1)%usewvl 1833 1834 if (me == 0) then 1835 ABI_MALLOC(pspfilnam_, (npsp)) 1836 1837 if (len_trim(pseudo_paths(1)) == 0) then 1838 ! Enter Legacy `files file` mode --> Read the name of the psp file from files file. 1839 1840 ! Catch possible mistake done by user (input without pseudos and `abinit t01.in` syntax) 1841 ! else the code starts to prompt for pseudos and execution gets stuck 1842 if (len_trim(input_path) /= 0) then 1843 ABI_ERROR("`pseudos` variable must be specified in input when the code is invoked with the `abinit t01.in` syntax") 1844 end if 1845 1846 ! Finish to read the "file" file completely, as npsp is known, 1847 do ipsp=1,npsp 1848 write(std_out,'(/,a)' )' Please give name of formatted atomic psp file (and finish with a newline character)' 1849 read (std_in, '(a)' , iostat=ios ) filpsp 1850 ! It might be that a file name is missing 1851 if (ios /= 0) then 1852 write(msg, '(5a)' )& 1853 'There are not enough names of pseudopotentials provided in the files file.',ch10,& 1854 'Action: check first the variable ntypat (and/or npsp) in the input file;',ch10,& 1855 'if they are correct, complete your files file.' 1856 ABI_ERROR(msg) 1857 end if 1858 pspfilnam_(ipsp) = trim(filpsp) 1859 write(std_out,'(a,i0,2a)' )' For atom type ',ipsp,', psp file is ',trim(filpsp) 1860 end do ! ipsp 1861 1862 else 1863 ! Get pseudopotential paths from input file. 1864 pspfilnam_ = pseudo_paths 1865 do ipsp=1,npsp 1866 write(std_out,'(a,i0,2a)' )' For atom type ',ipsp,', psp file is ',trim(pspfilnam_(ipsp)) 1867 end do 1868 end if 1869 1870 ! Now read the psp headers 1871 call inpspheads(pspfilnam_, npsp, pspheads, ecut_tmp) 1872 ABI_FREE(pspfilnam_) 1873 1874 if (minval(abs(pspheads(1:npsp)%pspcod - 7)) == 0) usepaw=1 1875 if (minval(abs(pspheads(1:npsp)%pspcod - 17)) == 0) usepaw=1 1876 end if ! me == 0 1877 1878 ABI_FREE(pseudo_paths) 1879 1880 ! Communicate pspheads to all processors 1881 call pspheads_comm(npsp, pspheads, usepaw) 1882 1883 ! If (all) pspcod are 7 then this is a PAW calculation. Initialize (default) the value of ratsph 1884 do idtset=0,ndtset_alloc 1885 dtsets(idtset)%usepaw = usepaw 1886 if (usepaw == 0) then 1887 dtsets(idtset)%ratsph(:)=two 1888 else 1889 ! Note that the following coding assumes that npsp=ntypat for PAW, which is true as of now (XG20101024). 1890 ! dtsets(idtset)%ratsph(1:npsp)=token%pspheads(1:npsp)%pawheader%rpaw 1891 do ipsp=1,npsp 1892 dtsets(idtset)%ratsph(ipsp) = pspheads(ipsp)%pawheader%rpaw 1893 end do 1894 endif 1895 end do 1896 1897 ! Take care of other dimensions, and part of the content of dtsets that is or might be needed early. 1898 ABI_MALLOC(zionpsp, (npsp)) 1899 do ii=1,npsp 1900 zionpsp(ii) = pspheads(ii)%zionpsp 1901 end do 1902 1903 ABI_MALLOC(mband_upper_, (0:ndtset_alloc)) 1904 1905 ! Get MAX dimension over datasets 1906 call invars1m(dmatpuflag, dtsets, ab_out, lenstr, mband_upper_, mx,& 1907 msym, ndtset, ndtset_alloc, string, npsp, zionpsp, comm) 1908 1909 ABI_FREE(zionpsp) 1910 call timab(1021,2,tsec) 1911 call timab(1022,3,tsec) 1912 1913 ! Provide defaults for the variables that have not yet been initialized. 1914 call indefo(dtsets, ndtset_alloc, nprocs) 1915 1916 call timab(1022,2,tsec) 1917 call timab(1023,3,tsec) 1918 1919 ! Perform some global initialization, depending on the value of 1920 ! pseudopotentials, parallelism variables, or macro input variables 1921 call macroin(dtsets, ecut_tmp, lenstr, ndtset_alloc, string) 1922 1923 ! If all the pseudopotentials have the same pspxc, override the default value for dtsets 1 to ndtset 1924 if (minval(abs((pspheads(1:npsp)%pspxc - pspheads(1)%pspxc)))==0) then 1925 dtsets(1:ndtset_alloc)%ixc = pspheads(1)%pspxc 1926 end if 1927 1928 ! Call the main input routine. 1929 call invars2m(dtsets,ab_out,lenstr,mband_upper_,msym,ndtset,ndtset_alloc,npsp,pspheads,string, comm) 1930 1931 call macroin2(dtsets, ndtset_alloc) 1932 1933 mx%mband = dtsets(1)%mband 1934 do ii=1,ndtset_alloc 1935 mx%mband = max(dtsets(ii)%mband, mx%mband) 1936 end do 1937 1938 call timab(1023,2,tsec) 1939 call timab(102,2,tsec) 1940 1941 ABI_FREE(mband_upper_) 1942 1943 end subroutine get_dtsets_pspheads
ABINIT/m_common [ Modules ]
NAME
m_common
FUNCTION
This module gathers routines used by higher-level procedures. Mainly printing routines.
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, AF, GMR, LBoeri, MT) 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
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_common 24 25 use defs_basis 26 use m_errors 27 use m_abicore 28 use m_exit 29 use m_fftcore 30 use m_fock 31 use m_io_tools 32 #if defined DEV_YP_VDWXC 33 use m_xc_vdw 34 #endif 35 #ifdef HAVE_NETCDF 36 use netcdf 37 #endif 38 use m_nctk 39 use m_crystal 40 use m_wfk 41 use m_ebands 42 use m_hdr 43 use m_xmpi 44 use m_dtset 45 use m_xpapi 46 use m_yaml 47 use m_invars2 48 use m_dtset 49 50 use m_fstrings, only : indent, endswith, sjoin, itoa 51 use m_electronpositron, only : electronpositron_type 52 use m_energies, only : energies_type, energies_eval_eint 53 use m_pair_list, only : pair_list 54 use m_geometry, only : mkrdim, metric 55 use m_kg, only : getcut 56 use m_parser, only : parsefile, ab_dimensions 57 use m_invars1, only : invars0, invars1m, indefo 58 use m_time, only : timab, time_set_papiopt 59 use defs_abitypes, only : MPI_type 60 use defs_datatypes, only : pspheader_type, ebands_t 61 use m_pspheads, only : inpspheads, pspheads_comm 62 use m_kpts, only : kpts_timrev_from_kptopt 63 64 implicit none 65 66 private
ABINIT/prteigrs [ Functions ]
NAME
prteigrs
FUNCTION
Print out eigenvalues band by band and k point by k point. If option=1, do it in a standard way, for self-consistent calculations. If option=2, print out residuals and eigenvalues, in a format adapted for nonself-consistent calculations, within the loops. If option=3, print out eigenvalues, in a format adapted for nonself-consistent calculations, at the end of the job. If option=4, print out derivatives of eigenvalues (same format as option==3, except header that is printed) If option=5, print out Fan contribution to zero-point motion correction to eigenvalues (averaged) (same format as option==3, except header that is printed) If option=6, print out DDW contribution to zero-point motion correction to eigenvalues (averaged) (same format as option==3, except header that is printed) If option=7, print out Fan+DDW contribution to zero-point motion correction to eigenvalues (averaged) (same format as option==3, except header that is printed)
INPUTS
eigen(mband*nkpt*nsppol)=eigenvalues (hartree) or, if option==4, diagonal of derivative of eigenvalues or, if option==5...7, zero-point motion correction to eigenvalues (averaged) enunit=choice parameter: 0=>output in hartree; 1=>output in eV; 2=> output in both hartree and eV fermie=fermi energy (Hartree) / for electrons thermalized in the conduction bands when occopt==9 fermih=fermi energy (Hartree) for holes thermalized in the VB when occopt==9 fname_eig=filename of printing the eigenenergies iout=unit number for formatted output file iscf=option for self-consistency kptns(3,nkpt)=k points in reduced coordinates kptopt=option for the generation of k points mband=maximum number of bands nband(nkpt)=number of bands at each k point nbdbuf= number of buffer bands nkpt=number of k points nnsclo_now=number of non-self-consistent loops for the current vtrial (often 1 for SCF calculation, =nstep for non-SCF calculations) nsppol=1 for unpolarized, 2 for spin-polarized occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point occopt=option for occupancies option= (see above) prteig=control print eigenenergies prtvol=control print volume and debugging resid(mband*nkpt*nsppol)=residuals (hartree**2) tolwfr=tolerance on band residual of wf, hartrees**2 (needed when option=2) vxcavg=average of vxc potential wtk(nkpt)=k-point weights
OUTPUT
(only writing)
SOURCE
1183 subroutine prteigrs(eigen,enunit,fermie,fermih,fname_eig,iout,iscf,kptns,kptopt,mband,nband,& 1184 & nbdbuf,nkpt,nnsclo_now,nsppol,occ,occopt,option,prteig,prtvol,resid,tolwfr,vxcavg,wtk) 1185 1186 !Arguments ------------------------------------ 1187 !scalars 1188 integer,intent(in) :: enunit,iout,iscf,kptopt,mband,nbdbuf,nkpt,nnsclo_now,nsppol 1189 integer,intent(in) :: occopt,option,prteig,prtvol 1190 real(dp),intent(in) :: fermie,fermih,tolwfr,vxcavg 1191 character(len=*),intent(in) :: fname_eig 1192 !arrays 1193 integer,intent(in) :: nband(nkpt*nsppol) 1194 real(dp),intent(in) :: eigen(mband*nkpt*nsppol),kptns(3,nkpt) 1195 real(dp),intent(in) :: occ(mband*nkpt*nsppol),resid(mband*nkpt*nsppol) 1196 real(dp),intent(in) :: wtk(nkpt) 1197 1198 !Local variables------------------------------- 1199 !scalars 1200 integer,parameter :: nkpt_max=50 1201 integer :: band_index,iband,ienunit,ii,ikpt,isppol,nband_index,nband_k,nkpt_eff,tmagnet,tmetal,temp_unit 1202 real(dp) :: convrt,magnet,residk,rhodn,rhoup 1203 character(len=2) :: ibnd_fmt,ikpt_fmt 1204 character(len=7) :: strunit1,strunit2 1205 character(len=39) :: kind_of_output 1206 character(len=500) :: msg 1207 1208 ! ************************************************************************* 1209 1210 if (enunit<0.or.enunit>2) then 1211 ABI_BUG(sjoin('enunit must be 0, 1 or 2. Argument was:', itoa(enunit))) 1212 end if 1213 1214 if (prteig > 0) then 1215 call wrtout(iout, sjoin(' prteigrs : about to open file ', fname_eig)) 1216 if (open_file(fname_eig, msg, newunit=temp_unit, status='unknown', form='formatted') /= 0) then 1217 ABI_ERROR(msg) 1218 end if 1219 rewind(temp_unit) ! always rewind disk file and print latest eigenvalues 1220 end if 1221 1222 kind_of_output= ' Eigenvalues ' 1223 if(option==4) kind_of_output=' Expectation of eigenvalue derivatives' 1224 if(option==5) kind_of_output=' Fan corrections to eigenvalues at T=0' 1225 if(option==6) kind_of_output=' DDW corrections to eigenvalues at T=0' 1226 if(option==7) kind_of_output=' Fan+DDW corrs to eigenvalues at T=0' 1227 1228 nkpt_eff=nkpt 1229 1230 !write(msg,'(a,5i5)')' prtvol,iscf,kptopt,nkpt_eff,nkpt_max ',prtvol,iscf,kptopt,nkpt_eff,nkpt_max 1231 !call wrtout(iout,msg) 1232 1233 if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max 1234 if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>1 .and. iout==ab_out)nkpt_eff=1 1235 1236 if(option==1 .or. (option>=3 .and. option<=7))then 1237 1238 do ienunit=0,1 1239 1240 if (enunit==1 .and. ienunit==0)cycle 1241 if (enunit==0 .and. ienunit==1)cycle 1242 ! Print eigenvalues in hartree for enunit=0 or 2 1243 ! The definition of two different strings is quite ridiculous. Historical reasons ... 1244 1245 if (ienunit==0)then 1246 convrt=one 1247 strunit1='hartree' 1248 strunit2='hartree' 1249 end if 1250 if (ienunit==1)then 1251 convrt=Ha_eV 1252 strunit1=' eV ' 1253 strunit2='eV ' 1254 end if 1255 1256 band_index=0 1257 1258 if(ienunit==0)then ! XG20140730 I do not know why this is only done when ienunit==0 1259 tmetal=0 1260 if(option==1 .and. occopt>=3 .and. occopt<=8)tmetal=1 1261 tmagnet=0 1262 if(tmetal==1 .and. nsppol==2)then 1263 tmagnet=1 1264 rhoup = 0._dp 1265 rhodn = 0._dp 1266 nband_index = 1 1267 do isppol=1,nsppol 1268 do ikpt=1,nkpt 1269 nband_k=nband(ikpt+(isppol-1)*nkpt) 1270 do iband=1,nband_k 1271 if(isppol==1) rhoup = rhoup + wtk(ikpt)*occ(nband_index) 1272 if(isppol==2) rhodn = rhodn + wtk(ikpt)*occ(nband_index) 1273 nband_index = nband_index + 1 1274 end do 1275 end do 1276 end do 1277 magnet = abs(rhoup - rhodn) 1278 end if 1279 end if 1280 1281 if(iscf>=0 .and. (ienunit==0 .or. option==1))then 1282 if (occopt == 9) then 1283 write(msg, '(3a,f10.5,a,f10.5,3a,f10.5)' ) & 1284 ' Fermi energy for thermalized electrons and holes (',trim(strunit2),') =',& 1285 convrt*fermie,', ',convrt*fermih,' Average Vxc (',trim(strunit2),')=',convrt*vxcavg 1286 else 1287 write(msg, '(3a,f10.5,3a,f10.5)' ) & 1288 ' Fermi (or HOMO) energy (',trim(strunit2),') =',convrt*fermie,' Average Vxc (',trim(strunit2),')=',convrt*vxcavg 1289 end if 1290 call wrtout(iout,msg) 1291 if (prteig > 0) call wrtout(temp_unit,msg) 1292 end if 1293 1294 ! if( (iscf>=0 .or. iscf==-3) .and. ienunit==0)then ! This is the most correct 1295 if(iscf>=0 .and. ienunit==0)then ! For historical reasons 1296 if(tmagnet==1)then 1297 write(msg, '(a,es16.8,a,a,es16.8,a,es16.8)' )& 1298 & ' Magnetization (Bohr magneton)=',magnet,ch10,& 1299 & ' Total spin up =',rhoup,' Total spin down =',rhodn 1300 call wrtout(iout,msg,'COLL') 1301 if (prteig > 0) call wrtout(temp_unit,msg) 1302 end if 1303 end if 1304 1305 ! Loop over spins (suppress spin data if nsppol not 2) 1306 do isppol=1,nsppol 1307 1308 ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9" 1309 if (nsppol==2.and.isppol==1) then 1310 write(msg, '(4a,'//ikpt_fmt//',2x,a)' )trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points, SPIN UP:' 1311 else if (nsppol==2.and.isppol==2) then 1312 write(msg, '(4a,'//ikpt_fmt//',2x,a)' )trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points, SPIN DOWN:' 1313 else 1314 write(msg, '(4a,'//ikpt_fmt//',2x,a)' )trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points:' 1315 end if 1316 call wrtout(iout,msg) 1317 if (prteig > 0) call wrtout(temp_unit,msg) 1318 1319 if(ienunit==0)then 1320 if(option>=4 .and. option<=7)then 1321 msg = ' (in case of degenerate eigenvalues, averaged derivative)' 1322 call wrtout(iout,msg) 1323 if (prteig > 0) call wrtout(temp_unit,msg) 1324 end if 1325 end if 1326 1327 do ikpt=1,nkpt 1328 nband_k=nband(ikpt+(isppol-1)*nkpt) 1329 ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9" 1330 ibnd_fmt="i3" ; if(nband_k>=1000)ibnd_fmt="i6" ; if(nband_k>=1000000)ibnd_fmt="i9" 1331 if(ikpt<=nkpt_eff)then 1332 write(msg, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) & 1333 & ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',kptns(1:3,ikpt)+tol10,' (reduced coord)' 1334 call wrtout(iout,msg,'COLL') 1335 if (prteig > 0) call wrtout(temp_unit,msg) 1336 do ii=0,(nband_k-1)/8 1337 write(msg, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index), iband=1+ii*8,min(nband_k,8+ii*8)) 1338 call wrtout(iout,msg,'COLL') 1339 if (prteig > 0) call wrtout(temp_unit,msg) 1340 end do 1341 if(ienunit==0 .and. option==1 .and. occopt>=3 .and. occopt<=8)then 1342 write(msg, '(5x,a,'//ikpt_fmt//')' ) ' occupation numbers for kpt#',ikpt 1343 call wrtout(iout,msg) 1344 do ii=0,(nband_k-1)/8 1345 write(msg, '(8(f10.5,1x))' ) (occ(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8)) 1346 call wrtout(iout,msg) 1347 end do 1348 end if 1349 1350 else 1351 if(ikpt==nkpt_eff+1)then 1352 write(msg, '(a,a)' )' prteigrs : prtvol=0 or 1, do not print more k-points.',ch10 1353 call wrtout(iout,msg) 1354 end if 1355 if (prteig > 0) then 1356 write(msg, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) & 1357 & ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',kptns(1:3,ikpt)+tol10,' (reduced coord)' 1358 call wrtout(temp_unit,msg) 1359 do ii=0,(nband_k-1)/8 1360 write(msg, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8)) 1361 call wrtout(temp_unit,msg) 1362 end do 1363 end if 1364 end if 1365 band_index=band_index+nband_k 1366 end do ! do ikpt=1,nkpt 1367 end do ! do isppol=1,nsppol 1368 1369 end do ! End loop over Hartree or eV 1370 1371 else if(option==2)then 1372 1373 band_index=0 1374 do isppol=1,nsppol 1375 1376 if(nsppol==2)then 1377 if(isppol==1)write(msg, '(2a)' ) ch10,' SPIN UP channel ' 1378 if(isppol==2)write(msg, '(2a)' ) ch10,' SPIN DOWN channel ' 1379 call wrtout(iout,msg) 1380 if(prteig>0) call wrtout(temp_unit,msg) 1381 end if 1382 1383 do ikpt=1,nkpt 1384 nband_k=nband(ikpt+(isppol-1)*nkpt) 1385 ikpt_fmt="i5" ; if(nkpt>=10000)ikpt_fmt="i7" ; if(nkpt>=1000000)ikpt_fmt="i9" 1386 1387 if(ikpt<=nkpt_eff)then 1388 write(msg, '(1x,a,'//ikpt_fmt//',a,f9.5,2f9.5,a)' ) & 1389 & 'Non-SCF case, kpt',ikpt,' (',(kptns(ii,ikpt),ii=1,3),'), residuals and eigenvalues=' 1390 call wrtout(iout,msg) 1391 if (prteig > 0) then 1392 write(msg, '(1x,a,'//ikpt_fmt//',a,f9.5,2f9.5,a)' ) & 1393 & 'Non-SCF case, kpt',ikpt,' eig(',(kptns(ii,ikpt),ii=1,3),') ' 1394 call wrtout(temp_unit,msg) 1395 end if 1396 do ii=0,(nband_k-1)/8 1397 write(msg, '(1p,8e10.2)' )(resid(iband+band_index),iband=1+8*ii,min(8+8*ii,nband_k)) 1398 call wrtout(iout,msg) 1399 end do 1400 do ii=0,(nband_k-1)/6 1401 write(msg, '(1p,6e12.4)' )(eigen(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k)) 1402 call wrtout(iout,msg) 1403 if (prteig > 0) call wrtout(temp_unit,msg) 1404 end do 1405 else 1406 if(ikpt==nkpt_eff+1)then 1407 write(msg, '(a,a)' )' prteigrs : prtvol=0 or 1, do not print more k-points.',ch10 1408 call wrtout(iout,msg) 1409 end if 1410 if (prteig > 0) then 1411 write(msg, '(1x,a,i5,a,f9.5,2f9.5,a)' )'Non-SCF kpt',ikpt,' eig(',(kptns(ii,ikpt),ii=1,3),') ' 1412 call wrtout(temp_unit,msg) 1413 do ii=0,(nband_k-1)/6 1414 write(msg, '(1p,6e12.4)' )(eigen(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k)) 1415 call wrtout(temp_unit,msg) 1416 end do 1417 end if 1418 end if 1419 1420 ! Don't include the buffer in the output. 1421 if (nbdbuf>0) then 1422 residk=maxval(resid(band_index+1:band_index+nband_k-nbdbuf)) 1423 else 1424 residk=maxval(resid(band_index+1:band_index+nband_k)) 1425 end if 1426 if (residk>tolwfr) then 1427 write(msg, '(1x,a,2i5,a,1p,e13.5)' ) & 1428 ' prteigrs : nnsclo,ikpt=',nnsclo_now,ikpt,' max resid (excl. the buffer)=',residk 1429 call wrtout(iout,msg) 1430 end if 1431 1432 band_index=band_index+nband_k 1433 end do 1434 end do 1435 call wrtout(iout," ") 1436 1437 else 1438 ABI_BUG(sjoin('option:', itoa(option),', is not allowed.')) 1439 end if 1440 1441 if (prteig > 0) close (temp_unit) 1442 1443 end subroutine prteigrs
ABINIT/prtene [ Functions ]
NAME
prtene
FUNCTION
Print components of total energy in nice format
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset | berryphase | kptopt | occopt | positron=option for electron-positron calculation | tphysel="physical" electronic temperature with FD occupations | tsmear=smearing energy or temperature (if metal) energies <type(energies_type)>=values of parts of total energy iout=unit number to which output is written usepaw= 0 for non paw calculation; =1 for paw calculation
OUTPUT
(only writing)
SOURCE
1471 subroutine prtene(dtset,energies,iout,usepaw) 1472 1473 !Arguments ------------------------------------ 1474 !scalars 1475 integer,intent(in) :: iout,usepaw 1476 type(dataset_type),intent(in) :: dtset 1477 type(energies_type),intent(in) :: energies 1478 1479 !Local variables------------------------------- 1480 !scalars 1481 integer :: ipositron,optdc 1482 logical :: directE_avail,testdmft 1483 real(dp) :: eent,enevalue,etotal,etotaldc,exc_semilocal 1484 ! Do not modify the length of these strings 1485 character(len=14) :: eneName 1486 character(len=500) :: msg 1487 type(yamldoc_t) :: edoc, dc_edoc 1488 !arrays 1489 !character(len=10) :: EPName(1:2)=(/"Positronic","Electronic"/) 1490 1491 ! ************************************************************************* 1492 1493 directE_avail=(usepaw==0.or.dtset%pawspnorb==0.or.dtset%pawcpxocc==2.or.dtset%kptopt==1.or.dtset%kptopt==2) 1494 1495 !============= Evaluate some parts of the energy =========== 1496 1497 optdc=-1;ipositron=merge(0,2,dtset%positron==0) 1498 if (abs(energies%e_ewald)<1.e-15_dp.and.abs(energies%e_hartree)<1.e-15_dp) ipositron=1 1499 call energies_eval_eint(energies,dtset,usepaw,optdc,etotal,etotaldc) 1500 1501 !Here, treat the case of metals 1502 !In re-smeared case the free energy is defined with tphysel 1503 if(dtset%occopt>=3 .and. dtset%occopt<=8)then 1504 if (abs(dtset%tphysel) < tol10) then 1505 eent=-dtset%tsmear * energies%entropy 1506 else 1507 eent=-dtset%tphysel * energies%entropy 1508 end if 1509 else 1510 eent=zero 1511 end if 1512 ! If DMFT is used and DMFT Entropy is not computed, then do not print 1513 ! non interacting entropy 1514 testdmft=(dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(dtset%upawu(:,1))>=tol8.or. & 1515 & sum(dtset%jpawu(:,1))>tol8).and.dtset%dmft_entropy==0) 1516 if(testdmft) eent=zero 1517 1518 etotal = etotal + eent 1519 etotaldc = etotaldc + eent 1520 1521 !============= Printing of Etotal by direct scheme =========== 1522 1523 if (dtset%icoulomb == 1) then 1524 eneName = "Ion-ion energy" 1525 else 1526 eneName = "Ewald energy" 1527 end if 1528 enevalue = energies%e_ewald 1529 1530 1531 if (optdc==0.or.optdc==2) then 1532 1533 if (directE_avail) then 1534 edoc = yamldoc_open('EnergyTerms', info='Components of total free energy in Hartree', & 1535 width=20, real_fmt='(es21.14)') 1536 call edoc%add_real('kinetic', energies%e_kinetic) 1537 if(abs(energies%e_extfpmd)>tiny(0.0_dp)) then 1538 call edoc%add_real('kinetic_extfpmd',energies%e_extfpmd) 1539 call edoc%add_real('total_kinetic',energies%e_extfpmd+energies%e_kinetic) 1540 end if 1541 if (ipositron/=1) then 1542 exc_semilocal=energies%e_xc+energies%e_hybcomp_E0-energies%e_hybcomp_v0+energies%e_hybcomp_v 1543 ! XG20181025 This should NOT be a part of the semilocal XC energy, but treated separately. 1544 ! At present, there is still a problem with the variational formulation for the Fock term with PAW. 1545 ! So, for the time being, keep it inside. 1546 if(usepaw==1)exc_semilocal=exc_semilocal+energies%e_fock 1547 call edoc%add_real('hartree', energies%e_hartree) 1548 call edoc%add_real('xc', exc_semilocal) 1549 call edoc%add_real(eneName, enevalue) 1550 call edoc%add_real('psp_core', energies%e_corepsp) 1551 #if defined DEV_YP_VDWXC 1552 if ( (dtset%vdw_xc > 0) .and. (dtset%vdw_xc < 10) .and. (xc_vdw_status()) ) then 1553 call edoc%add_real('VdWaals_df', energies%e_xc_vdw) 1554 end if 1555 #endif 1556 end if 1557 call edoc%add_real('local_psp', energies%e_localpsp) 1558 if (usepaw==0) then 1559 if(abs(energies%e_fock0)<tol8)then 1560 call edoc%add_real('non_local_psp', energies%e_nlpsp_vfock) 1561 else 1562 call edoc%add_real('non_local_psp+x', energies%e_nlpsp_vfock-energies%e_fock0) 1563 endif 1564 else 1565 call edoc%add_real('spherical_terms', energies%e_paw) 1566 !!!XG20181025 Does not work (yet)... 1567 !!!if(abs(energies%e_nlpsp_vfock)>tol8)then 1568 !!! write(msg, '(a,es21.14)' )' Fock-type term = ',energies%e_nlpsp_vfock 1569 !!! call wrtout(iout,msg,'COLL') 1570 !!! write(msg, '(a,es21.14)' ) ' -frozen Fock en.= ',-energies%e_fock0 1571 !!! call wrtout(iout,msg,'COLL') 1572 !!!endif 1573 end if 1574 if (ANY(ABS(dtset%nucdipmom)>tol8)) then 1575 call edoc%add_real('nucl. magn. dipoles',energies%e_nucdip) 1576 end if 1577 if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).and.ipositron/=1) then 1578 call edoc%add_real('VdWaals_dft_d', energies%e_vdw_dftd) 1579 end if 1580 if (dtset%nzchempot>=1) then 1581 call edoc%add_real('chem_potential', energies%e_chempot) 1582 end if 1583 if(dtset%occopt>=3.and.dtset%occopt<=8.and.ipositron==0) then 1584 call edoc%add_real('internal', etotal-eent) 1585 if(.not.testdmft) then 1586 call edoc%add_real('-kT*entropy', eent) 1587 end if 1588 else if (ipositron/=0) then 1589 if (dtset%occopt>=3.and.dtset%occopt<=8) then 1590 call edoc%add_real('-kT*entropy', eent) 1591 end if 1592 !write(msg, '(3a,es21.14,a)' ) & 1593 ! ' >>> ',EPName(ipositron),' E= ',etotal-energies%e0_electronpositron -energies%e_electronpositron,ch10 1594 !call wrtout(iout,msg,'COLL') 1595 !write(msg, '(3a,es21.14,2a,es21.14)' ) & 1596 ! ' ',EPName(3-ipositron),' ener.= ',energies%e0_electronpositron,ch10,& 1597 ! ' EP interaction E= ' ,energies%e_electronpositron 1598 !call wrtout(iout,msg,'COLL') 1599 if(ipositron == 1) then 1600 call edoc%add_real('positronic', etotal - energies%e0_electronpositron-energies%e_electronpositron) 1601 call edoc%add_real('electronic', energies%e0_electronpositron) 1602 else 1603 call edoc%add_real('electronic', etotal- energies%e0_electronpositron-energies%e_electronpositron) 1604 call edoc%add_real('positronic', energies%e0_electronpositron) 1605 end if 1606 call edoc%add_real('electron_positron_interaction', energies%e_electronpositron) 1607 end if 1608 if ((dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. & 1609 dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17) .and.ipositron/=1) then 1610 call edoc%add_real('electric', energies%e_elecfield) 1611 call edoc%add_real('kohn_sham', etotal-energies%e_elecfield) 1612 end if 1613 call edoc%add_real('total_energy', etotal) 1614 1615 else 1616 write(msg, '(9a)' ) & 1617 ' COMMENT: ',ch10,& 1618 ' "Direct" decomposition of total free energy cannot be printed out !!!',ch10,& 1619 ' PAW contribution due to spin-orbit coupling cannot be evaluated',ch10,& 1620 ' without the knowledge of imaginary part of Rhoij atomic occupancies',ch10,& 1621 ' (computed only when pawcpxocc=2).' 1622 call wrtout(iout,msg,'COLL') 1623 end if 1624 end if 1625 !============= Printing of Etotal by double-counting scheme =========== 1626 1627 if (optdc>=1) then 1628 1629 dc_edoc = yamldoc_open('EnergyTermsDC', info='"Double-counting" decomposition of free energy', & 1630 width=20, real_fmt="(es21.14)") 1631 call dc_edoc%add_real('band_energy', energies%e_eigenvalues) 1632 if(abs(energies%e_extfpmd)>tiny(0.0_dp)) then 1633 call dc_edoc%add_real('kinetic_extfpmd_dc',energies%edc_extfpmd) 1634 end if 1635 if (ipositron/=1) then 1636 !write(msg, '(2(a,es21.14,a),a,es21.14)' ) & 1637 ! ' '//eneName//' =',enevalue,ch10,& 1638 ! ' PspCore energy = ',energies%e_corepsp-energies%e_corepspdc,ch10,& 1639 ! ' Dble-C XC-energy= ',-energies%e_hartree+energies%e_xc-energies%e_xcdc -energies%e_fock0 + & 1640 ! energies%e_hybcomp_E0-energies%e_hybcomp_v0 1641 !call wrtout(iout,msg,'COLL') 1642 call dc_edoc%add_real(eneName, enevalue) 1643 call dc_edoc%add_real('psp_core', energies%e_corepsp-energies%e_corepspdc) 1644 call dc_edoc%add_real('xc_dc', -energies%e_hartree+energies%e_xc-energies%e_xcdc - energies%e_fock0 + & 1645 energies%e_hybcomp_E0-energies%e_hybcomp_v0) 1646 end if 1647 if ((dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. & 1648 dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17).and.ipositron/=1) then 1649 call dc_edoc%add_real('electric_field', energies%e_elecfield) 1650 end if 1651 if (usepaw==1) then 1652 call dc_edoc%add_real('spherical_terms', energies%e_pawdc) 1653 end if 1654 if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).and.ipositron/=1) then 1655 call dc_edoc%add_real('VdWaals_dft_d', energies%e_vdw_dftd) 1656 end if 1657 if (dtset%nzchempot>=1) then 1658 call dc_edoc%add_real('chem_potential', energies%e_chempot) 1659 end if 1660 if(dtset%occopt>=3.and.dtset%occopt<=8.and.ipositron==0) then 1661 if(.not.testdmft) then 1662 !write(msg, '(a,es21.14,a,a,a,es21.14)' ) & 1663 ! ' >>>>> Internal E= ',etotaldc-eent,ch10,ch10,& 1664 ! ' -kT*entropy = ',eent 1665 !call wrtout(iout,msg,'COLL') 1666 call dc_edoc%add_real('internal', etotaldc-eent) 1667 call dc_edoc%add_real('-kT*entropy', eent) 1668 else 1669 call dc_edoc%add_real('internal', etotaldc-eent) 1670 end if 1671 else if (ipositron/=0) then 1672 if (dtset%occopt>=3 .and. dtset%occopt<=8) then 1673 call dc_edoc%add_real('-kT*entropy', eent) 1674 end if 1675 !write(msg, '(a,es21.14,4a,es21.14,a)' ) & 1676 ! ' - EP dble-ct En.= ',-energies%edc_electronpositron,ch10,& 1677 ! ' >>> ',EPName(ipositron),' E= ',etotaldc-energies%e0_electronpositron -energies%e_electronpositron,ch10 1678 !call wrtout(iout,msg,'COLL') 1679 !write(msg, '(3a,es21.14,2a,es21.14)' ) & 1680 ! ' ',EPName(3-ipositron),' ener.= ',energies%e0_electronpositron,ch10,& 1681 ! ' EP interaction E= ' ,energies%e_electronpositron 1682 !call wrtout(iout,msg,'COLL') 1683 call dc_edoc%add_real('electron_positron_dc', -energies%edc_electronpositron) 1684 if(ipositron == 1) then 1685 call dc_edoc%add_real('positronic', etotaldc-energies%e0_electronpositron-energies%e_electronpositron) 1686 call dc_edoc%add_real('electronic', energies%e0_electronpositron) 1687 else 1688 call dc_edoc%add_real('electronic', etotaldc-energies%e0_electronpositron-energies%e_electronpositron) 1689 call dc_edoc%add_real('positronic', energies%e0_electronpositron) 1690 end if 1691 call dc_edoc%add_real('electron_positron_interaction', energies%e_electronpositron) 1692 end if 1693 1694 1695 write(msg, '(a,es21.14)' ) ' >>>> Etotal (DC)= ',etotaldc 1696 !call wrtout(iout,msg,'COLL') 1697 call dc_edoc%add_real('total_energy_dc', etotaldc) 1698 end if 1699 1700 !======= Additional printing for compatibility ========== 1701 1702 if (usepaw==0.and.optdc==0) then 1703 call edoc%add_real('total_energy_eV', etotal*Ha_eV) 1704 call edoc%add_real('band_energy', energies%e_eigenvalues) 1705 end if 1706 1707 if ((optdc==0.or.optdc==2).and.(.not.directE_avail)) then 1708 !write(msg, '(a,a,es18.10)' ) ch10,' Band energy (Ha)= ',energies%e_eigenvalues 1709 !call wrtout(iout,msg,'COLL') 1710 call edoc%add_real('band_energy', energies%e_eigenvalues) 1711 end if 1712 1713 if (usepaw==1) then 1714 if ((optdc==0.or.optdc==2).and.(directE_avail)) then 1715 call edoc%add_real('total_energy_eV', etotal*Ha_eV) 1716 end if 1717 if (optdc>=1) then 1718 !if (optdc==1) write(msg, '(a,a,es21.14)' ) ch10,' >Total DC energy in eV = ',etotaldc*Ha_eV 1719 !if (optdc==2) write(msg, '(a,es21.14)' ) ' >Total DC energy in eV = ',etotaldc*Ha_eV 1720 !call wrtout(iout,msg,'COLL') 1721 call dc_edoc%add_real('total_energy_dc_eV', etotaldc*Ha_eV) 1722 end if 1723 end if 1724 1725 if( dtset%icoulomb/=1.and.abs(dtset%cellcharge(1))>tol8) then 1726 write(msg, '(6a)' ) & 1727 ch10,' Calculation was performed for a charged system with PBC',& 1728 ch10,' You may consider including the monopole correction to the total energy',& 1729 ch10,' The correction is to be divided by the dielectric constant' 1730 call wrtout(iout,msg,'COLL') 1731 call edoc%add_real('monopole_correction', energies%e_monopole) 1732 call edoc%add_real('monopole_correction_eV', energies%e_monopole*Ha_eV) 1733 end if 1734 1735 ! Write components of total energies in Yaml format. 1736 call edoc%write_and_free(iout) 1737 if (optdc >= 1) call dc_edoc%write_and_free(iout) 1738 1739 end subroutine prtene
ABINIT/scprqt [ Functions ]
NAME
scprqt
FUNCTION
Conducts printing inside the scfcv.F90 routine, according to the value of choice. Also checks the convergence with respect to the different criteria. Eventually send a signal to quit the SCF cycle.
INPUTS
choice= if 1 => called at the initialisation of scfcv.f
if 2 => called during the loop in scfcv.f
if 3 => called at the end of scfcv.f
cpus=cpu time limit in seconds
deltae=change in energy between the previous and present SCF cycle
diffor=maximum absolute change in component of fcart between present and previous SCF cycle.
dtset <type(dataset_type)>=all input variables in this dataset
| chkexit= if non-zero, check whether the user wishes to exit
| enunit=parameter determining units of output energies
| ionmov=governs the movement of atoms (see help file)
| kptopt=option for the generation of k points
| mband=maximum number of bands
| natom=number of atoms in cell.
| nnsclo_now=number of non-self-consistent loops for the current vtrial
| (often 1 for SCF calculation, =nstep for non-SCF calculations)
| nsppol=1 for unpolarized, 2 for spin-polarized
| occopt=option for occupancies
| prtxml=1 if values have to be stored in an XML file.
| prteig=
| prtstm=print STM input variable
| prtvol= control print volume
| usedmatpu=DFT+U: number of SCF steps keeping occ. matrix fixed
| usefock=1 if Fock operator is present (hence possibility of a double loop)
| usepawu=0 if no DFT+U; /=0 if DFT+U
eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
etotal=total energy (hartree)
favg(3)=average of forces (ha/bohr)
fcart(3,natom)=cartesian forces (hartree/bohr)
fermie=fermi energy (Hartree) / for electrons thermalized in the conduction bands when occopt==9
fermih=fermi energy (Hartree) for holes thermalized in the VB when occopt==9
fname_eig=filename for printing of the eigenenergies
fock <type(fock_type)>=quantities for the fock operator (optional argument)
character(len=fnlen) :: filnam1=character strings giving input file name
initGS= 1 if one GS SCF cycle has already be done
iscf=( <= 0 =>non-SCF), >0 => SCF)
iscf =1 => determination of the largest eigenvalue of the SCF cycle
iscf =2 => SCF cycle, simple mixing
iscf =3 => SCF cycle, anderson mixing
iscf =5 => SCF cycle, CG based on estimations of gradients of the energy
iscf =6 => SCF cycle, CG based on true minimization of the energy
iscf =-3, although non-SCF, the energy is computed, so print it here.
istep=number of the SCF iteration (needed if choice=2)
istep_fock_outer=number of outer SCF iteration in the double loop approach
istep_mix=number of inner SCF iteration in the double loop approach
kpt(3,nkpt)=reduced coordinates of k points.
maxfor=maximum absolute value of fcart
moved_atm_inside: if==1, the atoms are allowed to move.
mpi_enreg=information about MPI parallelization
nband(nkpt*nsppol)=number of bands at each k point, for each polarization
nkpt=number of k points
nstep=number of steps expected in iterations.
occ(mband*nkpt*nsppol)=occupation number for each band at each k point.
optres=0 if the residual (res2) is a POTENTIAL residual
1 if the residual (res2) is a DENSITY residual
prtfor=1 only if forces have to be printed (0 otherwise)
prtxml=1 if XML file has to be output
res2=square of the density/potential residual
resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins
residm=maximum value from resid array (except for nbdbuf highest bands)
in Wavelets mode, it is used as the maximum value for the gradient norm.
response= if 0, GS case, if 1, RF case.
tollist(12)=tolerance list. Presently, the following are defined :
tollist(1)=tolmxf ; tollist(2)=tolwfr ; tollist(3)=toldff
tollist(4)=toldfe ; tollist(5)=toleig ; tollist(6)=tolvrs
tollist(7)=tolrff
usepaw= 0 for non paw calculation; =1 for paw calculation
vxcavg=mean of the vxc potential
wtk(nkpt)=weight assigned to each k point.
xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
quit= 0 if the SCF cycle is not finished; 1 otherwise. conv_retcode=Only if choice==3, != 0 if convergence is not achieved.
SOURCE
170 subroutine scprqt(choice,cpus,deltae,diffor,dtset,& 171 & eigen,etotal,favg,fcart,fermie,fermih,fname_eig,filnam1,initGS,& 172 & iscf,istep,istep_fock_outer,istep_mix,kpt,maxfor,moved_atm_inside,mpi_enreg,& 173 & nband,nkpt,nstep,occ,optres,& 174 & prtfor,prtxml,quit,res2,resid,residm,response,tollist,usepaw,& 175 & vxcavg,wtk,xred,conv_retcode,& 176 & electronpositron, fock) ! optional arguments) 177 178 !Arguments ------------------------------------ 179 !scalars 180 integer,intent(in) :: choice,initGS,iscf,istep,istep_fock_outer,istep_mix 181 integer,intent(in) :: moved_atm_inside,nkpt,nstep 182 integer,intent(in) :: optres,prtfor,prtxml,response,usepaw 183 integer,intent(out) :: quit,conv_retcode 184 real(dp),intent(in) :: cpus,deltae,diffor,etotal,fermie,fermih,maxfor,res2,residm 185 real(dp),intent(in) :: vxcavg 186 character(len=fnlen),intent(in) :: fname_eig,filnam1 187 type(electronpositron_type),pointer,optional :: electronpositron 188 type(fock_type),pointer,optional :: fock 189 type(MPI_type),intent(in) :: mpi_enreg 190 type(dataset_type),intent(in) :: dtset 191 !arrays 192 integer,intent(in) :: nband(nkpt*dtset%nsppol) 193 real(dp),intent(in) :: eigen(dtset%mband*nkpt*dtset%nsppol),favg(3) 194 real(dp),intent(in) :: fcart(3,dtset%natom),kpt(3,nkpt) 195 real(dp),intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol) 196 real(dp),intent(in) :: resid(dtset%mband*nkpt*dtset%nsppol),tollist(12) 197 real(dp),intent(in) :: wtk(nkpt),xred(3,dtset%natom) 198 199 !Local variables------------------------------- 200 !scalars 201 integer,parameter :: master=0 202 integer,save :: toldfe_ok,toldff_ok,tolrff_ok,ttoldfe,ttoldff,ttolrff,ttolvrs,ttolwfr 203 integer :: iatom,iband,iexit,ikpt,ii,ishift,isppol,my_rank 204 integer :: nband_index,nband_k,nnsclohf 205 integer :: openexit,option,tmagnet,usefock 206 #if defined DEV_YP_VDWXC 207 integer :: ivdw 208 #endif 209 real(dp),save :: toldfe,toldff,tolrff,tolvrs,tolwfr,vdw_df_threshold 210 real(dp) :: diff_e,diff_f,magnet,rhodn,rhoup 211 logical :: noquit,use_dpfft 212 character(len=500) :: message, message2, message3 213 character(len=2) :: format_istep 214 character(len=5) :: format_magnet 215 character(len=8) :: colname 216 character(len=1) :: firstchar 217 type(yamldoc_t) :: ydoc 218 !arrays 219 real(dp) :: residm_band(dtset%mband,dtset%nsppol), f_tmp(3) 220 221 ! ********************************************************************* 222 223 DBG_ENTER("COLL") 224 225 my_rank = mpi_enreg%me_cell 226 227 quit=0; conv_retcode=0 228 usefock=dtset%usefock 229 nnsclohf=dtset%nnsclohf 230 use_dpfft = .False. 231 232 tmagnet=0 233 if(response==0.and.(iscf>0.or.iscf==-3).and.dtset%nsppol==2.and.dtset%occopt>2)tmagnet=1 234 235 ishift=0 236 residm_band = zero 237 do isppol=1, dtset%nsppol 238 do ikpt=1, nkpt 239 do iband=1, nband(ikpt+(isppol-1)*nkpt) 240 ishift = ishift+1 241 residm_band(iband, isppol) = max (resid(ishift), residm_band(iband, isppol)) 242 end do 243 end do 244 end do 245 246 select case (choice) 247 case (1) 248 ! choice= if 1 => called at the initialisation of scfcv.f 249 ! Examine tolerance criteria 250 ! NB: The tests on tolwfr and the presence of tolerances in the SCF case are 251 ! also done at the level of the parser in chkinp. 252 tolwfr=tollist(2) 253 toldff=tollist(3) 254 toldfe=tollist(4) 255 tolvrs=tollist(6) 256 tolrff=tollist(7) 257 vdw_df_threshold=tollist(8) 258 ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0; 259 if(abs(tolwfr)>tiny(zero))ttolwfr=1 260 if(abs(toldff)>tiny(zero))ttoldff=1 261 if(abs(tolrff)>tiny(zero))ttolrff=1 262 if(abs(toldfe)>tiny(zero))ttoldfe=1 263 if(abs(tolvrs)>tiny(zero))ttolvrs=1 264 ! If non-scf calculations, tolwfr must be defined 265 if(ttolwfr /= 1 .and. (iscf<0 .and. iscf/=-3) )then 266 write(message,'(a,a,a,es14.6,a,a)')& 267 'when iscf <0 and /= -3, tolwfr must be strictly',ch10,& 268 'positive, while it is ',tolwfr,ch10,& 269 'Action: change tolwfr in your input file and resubmit the job.' 270 ABI_ERROR(message) 271 end if 272 ! toldff only allowed when prtfor==1 273 ! FIXME: this test should be done on input, not during calculation 274 if((ttoldff == 1 .or. ttolrff == 1) .and. prtfor==0 )then 275 ABI_ERROR('toldff only allowed when prtfor=1!') 276 end if 277 ! If SCF calculations, one and only one of these can differ from zero 278 if( (iscf>0 .or. iscf==-3) .and.(ttolwfr==1.and.ttoldff+ttoldfe+ttolvrs+ttolrff>1) & 279 .and. (ttolwfr==0.and.ttoldff+ttoldfe+ttolvrs+ttolrff/=1) ) then 280 write(message,'(6a,es14.6,a,es14.6,a,es14.6,a,a,es14.6,a,a,a)' )& 281 & 'For the SCF case, one and only one of the input tolerance criteria ',ch10,& 282 & 'toldff, tolrff, toldfe or tolvrs ','must differ from zero, while they are',ch10,& 283 & 'toldff=',toldff,', tolrff=',tolrff,', toldfe=',toldfe,ch10,& 284 & 'and tolvrs=',tolvrs,' .',ch10,& 285 & 'Action: change your input file and resubmit the job.' 286 ABI_ERROR(message) 287 end if 288 289 if (dtset%usewvl == 1) then 290 write(colname, "(A)") "grdnorm " 291 else 292 write(colname, "(A)") "residm " 293 end if 294 if (nstep>0 .and. (iscf>=0 .or.iscf==-3) .and. dtset%prtstm==0) then 295 if(tmagnet==1)then 296 if (prtfor==0) then 297 if (optres==0) then 298 write(message, '(4a)' ) ch10,& 299 & ' iter Etot(hartree) deltaE(h) ',colname,' vres2 magn' 300 else 301 write(message, '(4a)' ) ch10,& 302 & ' iter Etot(hartree) deltaE(h) ',colname,' nres2 magn' 303 end if 304 else 305 if (optres==0) then 306 write(message, '(4a)' ) ch10,& 307 & ' iter Etot(hartree) deltaE(h) ',colname,' vres2 diffor maxfor magn' 308 else 309 write(message, '(4a)' ) ch10,& 310 & ' iter Etot(hartree) deltaE(h) ',colname,' nres2 diffor maxfor magn' 311 end if 312 end if 313 else 314 if(response==0)then 315 if (prtfor==0) then 316 if (optres==0) then 317 write(message, '(4a)' ) ch10,& 318 ' iter Etot(hartree) deltaE(h) ', colname, ' vres2' 319 else 320 write(message, '(4a)' ) ch10,& 321 ' iter Etot(hartree) deltaE(h) ', colname, ' nres2' 322 end if 323 else 324 if (optres==0) then 325 write(message, '(4a)' ) ch10,& 326 ' iter Etot(hartree) deltaE(h) ',colname,' vres2 diffor maxfor ' 327 else 328 write(message, '(4a)' ) ch10,& 329 ' iter Etot(hartree) deltaE(h) ',colname,' nres2 diffor maxfor ' 330 end if 331 end if 332 else 333 if (optres==0) then 334 write(message, '(4a)' ) ch10,& 335 ' iter 2DEtotal(Ha) deltaE(Ha) ', colname, ' vres2' 336 else 337 write(message, '(4a)' ) ch10,& 338 ' iter 2DEtotal(Ha) deltaE(Ha) ', colname, ' nres2' 339 end if 340 end if 341 end if 342 343 ydoc = yamldoc_open('BeginCycle') 344 call ydoc%add_ints("iscf, nstep, nline, wfoptalg", & 345 [dtset%iscf, dtset%nstep, dtset%nline, dtset%wfoptalg], dict_key="solver") 346 call ydoc%add_reals("tolwfr, toldff, toldfe, tolvrs, tolrff", & ! , vdw_df_threshold", & 347 [tolwfr, toldff, toldfe, tolvrs, tolrff], & !, vdw_df_threshold], & 348 real_fmt="(es8.2)", dict_key="tolerances", ignore=zero) 349 350 !if (dtset%use_yaml == 1) then 351 call ydoc%write_and_free(ab_out, newline=.False.) 352 !else 353 !call ydoc%write_and_free(std_out, newline=.False.) 354 !end if 355 356 call wrtout(ab_out,message,'COLL') 357 end if 358 359 case (2) 360 361 ! Examine tolerance criteria 362 tolwfr=tollist(2) 363 toldff=tollist(3) 364 toldfe=tollist(4) 365 tolvrs=tollist(6) 366 tolrff=tollist(7) 367 vdw_df_threshold=tollist(8) 368 ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0; 369 if(abs(tolwfr)>tiny(0.0_dp))ttolwfr=1 370 if(abs(toldff)>tiny(0.0_dp))ttoldff=1 371 if(abs(tolrff)>tiny(0.0_dp))ttolrff=1 372 if(abs(toldfe)>tiny(0.0_dp))ttoldfe=1 373 if(abs(tolvrs)>tiny(0.0_dp))ttolvrs=1 374 375 ! Conduct printing. If extra output follows, then put a blank line into the output here 376 if (dtset%prtvol>=10) call wrtout([std_out, ab_out], ' ') 377 378 ! Calculate up and down charge and magnetization 379 if(tmagnet==1) then 380 rhoup = zero 381 rhodn = zero 382 nband_index = 1 383 do isppol=1,dtset%nsppol 384 do ikpt=1,nkpt 385 nband_k=nband(ikpt+(isppol-1)*nkpt) 386 do iband=1,nband_k 387 if(isppol==1) rhoup = rhoup + wtk(ikpt)*occ(nband_index) 388 if(isppol==2) rhodn = rhodn + wtk(ikpt)*occ(nband_index) 389 nband_index = nband_index + 1 390 end do 391 end do 392 end do 393 magnet = abs(rhoup - rhodn) 394 end if 395 396 if (prtxml == 1) then 397 write(ab_xml_out, "(A)", advance = "NO") ' <scfcvStep' 398 write(message, "(es22.10)") etotal 399 write(ab_xml_out, "(A,A,A)", advance = "NO") ' eTotal="', trim(message) ,'"' 400 write(message, "(es20.8)") deltae 401 write(ab_xml_out, "(A,A,A)", advance = "NO") ' deltaETotal="', trim(message) ,'"' 402 write(message, "(es20.8)") residm 403 write(ab_xml_out, "(A,A,A)", advance = "NO") ' maxResid="', trim(message) ,'"' 404 write(message, "(es20.8)") res2 405 if (optres == 0) then 406 write(ab_xml_out, "(A,A,A)", advance = "NO") ' potResid="', trim(message) ,'"' 407 else 408 write(ab_xml_out, "(A,A,A)", advance = "NO") ' denResid="', trim(message) ,'"' 409 end if 410 if (tmagnet== 1) then 411 write(message, "(es20.8)") magnet 412 write(ab_xml_out, "(A,A,A)", advance = "NO") ' magn="', trim(message) ,'"' 413 end if 414 if (prtfor == 1) then 415 write(message, "(es20.8)") diffor 416 write(ab_xml_out, "(A,A,A)", advance = "NO") ' deltaForces="', trim(message) ,'"' 417 write(message, "(es20.8)") maxfor 418 write(ab_xml_out, "(A,A,A)", advance = "NO") ' maxForces="', trim(message) ,'"' 419 end if 420 write(ab_xml_out, "(A)") " />" 421 end if 422 423 ! Print total (free) energy (hartree) and other convergence measures 424 if(dtset%prtstm==0)then 425 format_istep='i3' 426 if(istep>99)format_istep='i5' 427 if(istep>9999)format_istep='i7' 428 if(tmagnet==1)then 429 if(magnet<10)then 430 format_magnet='f6.3)' 431 else if(magnet<100)then 432 format_magnet='f6.2)' 433 else 434 format_magnet='f6.1)' 435 end if 436 if (prtfor==0) then 437 write(message, '(a,'//format_istep//',1p,g22.14,3es9.2,0p,'//format_magnet ) & 438 & ' ETOT',istep,etotal,deltae,residm,res2,magnet 439 else 440 write(message, '(a,'//format_istep//',1p,g22.14,3es9.2,es8.1,es9.2,0p,'//format_magnet ) & 441 & ' ETOT',istep,etotal,deltae,residm,res2,diffor,maxfor,magnet 442 end if 443 else 444 firstchar=' ' 445 if (response/=0.and.istep==1) firstchar="-" 446 if (response==0) then 447 if (prtfor==0) then 448 write(message, '(2a,'//format_istep//',1p,g22.14,3es10.3)' ) & 449 & firstchar,'ETOT',istep,etotal,deltae,residm,res2 450 else 451 write(message, '(2a,'//format_istep//',1p,g22.14,5es10.3)' ) & 452 & firstchar,'ETOT',istep,etotal,deltae,residm,res2,diffor,maxfor 453 end if 454 else 455 write(message, '(2a,'//format_istep//',1p,g22.14,1x,3es10.3)' ) & 456 & firstchar,'ETOT',istep,etotal,deltae,residm,res2 457 end if 458 end if 459 !if (etot_yaml_doc%stream%length /= 0) call etot_yaml_doc%add_tabular_line(' '//message(6:)) 460 call wrtout(ab_out,message,'COLL') 461 462 if(mpi_enreg%paral_pert==1) then 463 call wrtout(std_out, message,'PERS') 464 elseif(mpi_enreg%paral_pert==0) then 465 call wrtout(std_out, message,'COLL') 466 end if 467 468 end if ! dtset%prtstm==0 469 470 ! Print positions/forces every step if dtset%prtvol>=10 and iscf>0 or -3 and GS case 471 if (dtset%prtvol>=10.and.(iscf>=0.or.iscf==-3).and.response==0.and.dtset%prtstm==0) then 472 call wrtout(ab_out," ",'COLL') 473 474 ! Print up and down charge and magnetization 475 if(tmagnet==1) then 476 write(message,'(a,f11.6,a,f11.6,a,f10.6)')& 477 & ' #electrons spin up=',rhoup,', spin down=',rhodn,', magnetization=',magnet 478 call wrtout([std_out, ab_out], message) 479 end if 480 481 ! Moreover, print atomic positions if dtset%ionmov==4, and moved_atm_inside==1 482 if (dtset%ionmov==4 .and. moved_atm_inside==1)then 483 call wrtout([std_out, ab_out], ' reduced coordinates :') 484 do iatom=1,dtset%natom 485 write(message, '(i5,1x,3es21.11)' ) iatom,xred(:,iatom) 486 call wrtout([std_out, ab_out], message) 487 end do 488 end if 489 490 ! Slightly change favg for printing reasons 491 if (prtfor>0) then 492 f_tmp(:)=favg(:) 493 if(abs(favg(1))<1.0d-13)f_tmp(1)=zero 494 if(abs(favg(2))<1.0d-13)f_tmp(2)=zero 495 if(abs(favg(3))<1.0d-13)f_tmp(3)=zero 496 write(message, '(a,3es10.2)' )' cartesian forces (ha/bohr); non-corrected avg=',f_tmp(:) 497 call wrtout([std_out, ab_out], message) 498 do iatom=1,dtset%natom 499 f_tmp(:)=fcart(:,iatom) 500 if(abs(fcart(1,iatom))<1.0d-13)f_tmp(1)=zero 501 if(abs(fcart(2,iatom))<1.0d-13)f_tmp(2)=zero 502 if(abs(fcart(3,iatom))<1.0d-13)f_tmp(3)=zero 503 write(message, '(i5,1x,3es21.11)' ) iatom,f_tmp(:) 504 call wrtout([std_out, ab_out], message) 505 end do 506 end if 507 508 end if 509 510 ! Print eigenvalues every step if dtset%prtvol>=10 and GS case 511 if (my_rank == master .and. (dtset%prtvol>=10 .and. response==0 .and. dtset%tfkinfunc==0 .and. dtset%usewvl==0)) then 512 option=1 513 call prteigrs(eigen,dtset%enunit,fermie,fermih,fname_eig,ab_out,iscf,kpt,dtset%kptopt,dtset%mband,& 514 & nband,dtset%nbdbuf,nkpt,dtset%nnsclo,dtset%nsppol,occ,dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk) 515 call prteigrs(eigen,dtset%enunit,fermie,fermih,fname_eig,std_out,iscf,kpt,dtset%kptopt,dtset%mband,& 516 & nband,dtset%nbdbuf,nkpt,dtset%nnsclo,dtset%nsppol,occ,dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk) 517 end if 518 519 if(response==0)then 520 write(message, '(a,1p,e15.7,a)' ) ' scprqt: <Vxc>=',vxcavg,' hartree' 521 call wrtout(std_out,message,'COLL') 522 end if 523 524 ! Check whether exiting was required by the user. 525 openexit=1 ; if(dtset%chkexit==0) openexit=0 526 call exit_check(cpus,filnam1,iexit,ab_out,mpi_enreg%comm_cell,openexit) 527 if (iexit/=0) quit=1 528 529 ! In special cases, do not quit even if convergence is reached 530 noquit=((istep<nstep).and.(usepaw==1).and.(dtset%usepawu/=0).and.& 531 & (dtset%usedmatpu/=0).and.(istep<=abs(dtset%usedmatpu)).and.& 532 & (dtset%usedmatpu<0.or.initGS==0)) 533 534 ! Additional stuff for electron/positron 535 if (present(electronpositron)) then 536 if (associated(electronpositron)) then 537 if (electronpositron%istep_scf==1) then 538 toldff_ok=0;tolrff_ok=0;toldfe_ok=0 539 end if 540 end if 541 end if 542 543 ! Stopping criteria in the SCF case 544 if(iscf>1 .or. iscf==-3 .or. iscf == 0) then 545 ! Here treat the vdw_df_threshold criterion : if the change of energy is less than 546 ! input vdw_df_threshold, trigger the calculation of vdW interactions 547 ! write(message,'(1x,a,e10.3,1x,a,e10.3,1x,l1,a)') & 548 ! & '[vdW-DF][DEBUG] deltae=',deltae,'vdw_df_threshold=',vdw_df_threshold, & 549 ! & (abs(deltae)<vdw_df_threshold),ch10 550 ! call wrtout(std_out,message,'COLL') 551 #if defined DEV_YP_VDWXC 552 call xc_vdw_trigger( (abs(deltae)<vdw_df_threshold) ) 553 #endif 554 ! Here treat the tolwfr criterion: if maximum residual is less than 555 ! input tolwfr, stop steps (exit loop here) 556 if (ttolwfr == 1 .and. (ttolvrs+ttoldfe+ttoldff+ttolrff==0) .and. .not. noquit) then 557 if (residm < tolwfr) then 558 if (dtset%usewvl == 0) then 559 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a)' )ch10, & 560 ' At SCF step',istep,' max residual=',residm,' < tolwfr=',tolwfr,' =>converged.' 561 else 562 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a)' )ch10, & 563 ' At SCF step',istep,' max grdnorm=',residm,' < tolwfr=',tolwfr,' =>converged.' 564 end if 565 call wrtout([std_out, ab_out], message) 566 quit=1 567 else 568 use_dpfft = residm < tol7 569 end if 570 end if 571 572 ! Here treat the toldff criterion: if maximum change of fcart is less than 573 ! input toldff twice consecutively, stop steps (exit loop here) 574 if (ttoldff==1) then 575 if (istep==1) then 576 toldff_ok=0 577 else if (diffor < toldff) then 578 toldff_ok=toldff_ok+1 579 ! add warning for forces which are 0 by symmetry. Also added Matteo check below that the wave 580 ! functions are relatively converged as well 581 if (diffor < tol12) then 582 write (message,'(3a)') ' toldff criterion is satisfied, but your forces are suspiciously low.', ch10,& 583 & ' Check if the forces are 0 by symmetry: in that case you can not use the toldff convergence criterion!' 584 ABI_WARNING(message) 585 if (maxfor < tol16 .and. res2 > tol9) tolrff_ok=0 586 end if 587 else 588 toldff_ok=0 589 use_dpfft = diffor < tol6 590 end if 591 592 if(toldff_ok>=2 .and..not.noquit)then 593 if (ttolwfr==0) then 594 write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3)' ) ch10, & 595 & ' At SCF step',istep,', forces are converged : ',ch10,& 596 & ' for the second time, max diff in force=',diffor,' < toldff=',toldff 597 call wrtout([std_out, ab_out], message) 598 quit=1 599 else if (ttolwfr==1 .and. residm < tolwfr )then 600 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a,a,es11.3,a,es11.3)' ) ch10, & 601 & ' At SCF step',istep,', max residual=',residm,' < tolwfr=',tolwfr,' AND forces are converged : ',ch10,& 602 & ' for the second time, max diff in force=',diffor,' < toldff=',toldff 603 call wrtout([std_out, ab_out], message) 604 quit=1 605 end if 606 end if 607 end if 608 609 ! Here treat the tolrff criterion: if maximum change of fcart is less than 610 ! input tolrff times fcart itself twice consecutively, stop steps (exit loop here) 611 if (ttolrff==1) then 612 if (istep==1) then 613 tolrff_ok=0 614 ! 27/7/2009: added test for absolute value of maxfor, otherwise if it is 0 this never exits the scf loop. 615 else if (diffor < tolrff*maxfor .or. (maxfor < tol16 .and. diffor < tol16)) then 616 tolrff_ok=tolrff_ok+1 617 ! Thu Mar 12 19:01:40 MG: added additional check on res2 to make sure the SCF cycle is close to convergence. 618 ! Needed for structural relaxations otherwise the stress tensor is wrong and the relax algo makes wrong moves. 619 if (maxfor < tol16 .and. res2 > tol9) tolrff_ok=0 620 else 621 tolrff_ok=0 622 use_dpfft = diffor < tolrff * maxfor * five 623 end if 624 if(tolrff_ok>=2 .and. (.not.noquit))then 625 if (ttolwfr==0) then 626 write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3,a)' ) ch10, & 627 ' At SCF step',istep,', forces are sufficiently converged : ',ch10,& 628 ' for the second time, max diff in force=',diffor,& 629 ' is less than < tolrff=',tolrff, ' times max force' 630 call wrtout([std_out, ab_out], message) 631 quit=1 632 else if (ttolwfr==1 .and. residm < tolwfr) then 633 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a,a,es11.3,a,es11.3,a)' ) ch10, & 634 ' At SCF step',istep,', max residual=',residm,' < tolwfr=',tolwfr,' AND forces are sufficiently converged : ',ch10,& 635 ' for the second time, max diff in force=',diffor,& 636 ' is less than < tolrff=',tolrff, ' times max force' 637 call wrtout([std_out, ab_out], message) 638 quit=1 639 end if 640 end if 641 end if 642 643 ! Here treat the toldfe criterion: if the change of energy is less than 644 ! input toldfe twice consecutively, stop steps (exit loop here) 645 if (ttoldfe==1) then 646 if (istep==1) then 647 toldfe_ok=0 648 else if (abs(deltae)<toldfe) then 649 toldfe_ok=toldfe_ok+1 650 else 651 toldfe_ok=0 652 use_dpfft = abs(deltae) < tol8 653 end if 654 ! Fock : tolwfr not taken into account 655 if(usefock/=0.and.nnsclohf>=2) then 656 if (toldfe_ok==2 .and. (.not.noquit))then 657 write(message, '(a,i3,a,i3,a,a,a,es11.3,a,es11.3)' ) & 658 ' Outer loop step',istep_fock_outer,' - inner step',istep_mix,' - frozen Fock etot converged : ',ch10,& 659 ' for the second time, diff in etot=',abs(deltae),' < toldfe=',toldfe 660 call wrtout([std_out, ab_out], message) 661 quit=1 662 end if 663 ! No Fock : take into account tolwfr if ttolwfr/=0 664 else if(toldfe_ok>=2 .and. (.not.noquit))then 665 if (ttolwfr==0) then 666 write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3)' ) ch10, & 667 ' At SCF step',istep,', etot is converged : ',ch10,& 668 ' for the second time, diff in etot=',abs(deltae),' < toldfe=',toldfe 669 call wrtout([std_out, ab_out], message) 670 quit=1 671 else if (ttolwfr==1 .and. residm < tolwfr) then 672 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a,a,es11.3,a,es11.3)' ) ch10, & 673 ' At SCF step',istep,', max residual=',residm,' < tolwfr=',tolwfr,' AND etot is converged : ',ch10,& 674 ' for the second time, diff in etot=',abs(deltae),' < toldfe=',toldfe 675 call wrtout([std_out, ab_out], message) 676 quit=1 677 end if 678 end if 679 if(usefock==1 .and. nnsclohf>1)then 680 if(istep_mix==1 .and. (.not.noquit))then 681 ! The change due to the update of the Fock operator is sufficiently small. No need to meet it a second times. 682 if (abs(deltae)<toldfe) then 683 write(message, '(a,i3,a,i3,a,a,a,es11.3,a,es11.3)' ) & 684 ' Outer loop step',istep_fock_outer,' - inner step',istep_mix,' - etot converged : ',ch10,& 685 ' update of Fock operator yields diff in etot=',abs(deltae),' < toldfe=',toldfe 686 call wrtout([std_out, ab_out], message) 687 fock%fock_common%fock_converged=.true. 688 quit=1 689 endif 690 endif 691 !TODO: separate messages: if HF is imposing a continuation of the loop, then abs(deltae) is actually not > toldfe 692 if(istep_mix==nnsclohf .and. quit==0)then 693 write(message, '(a,i3,a,i3,a,a,a,es11.3,a,es11.3)' ) & 694 ' Outer loop step',istep_fock_outer,' - inner step',istep_mix,' - frozen Fock etot NOT converged : ',ch10,& 695 ' diff in etot=',abs(deltae),' > toldfe=',toldfe 696 call wrtout([std_out, ab_out], message) 697 endif 698 endif 699 700 ! Here treat the vdw_df_threshold criterion for non-SCF vdW-DF 701 ! calculations: If input vdw_df_threshold is lesss than toldfe 702 ! then the vdW-DF is triggered once selfconsistency criteria is 703 ! reached for the first time. 704 ! write(message,'(1x,a,e10.3,1x,a,e10.3,1x,l1,a)') & 705 ! & '[vdW-DF][DEBUG] deltae=',deltae,'vdw_df_threshold=',vdw_df_threshold, & 706 ! & (abs(deltae)<toldfe),ch10 707 ! call wrtout(std_out,message,'COLL') 708 #if defined DEV_YP_VDWXC 709 ivdw = 0 710 if ( toldfe > vdw_df_threshold ) then 711 ivdw = ivdw + 1 712 end if 713 call xc_vdw_trigger((toldfe_ok==1 .and. toldfe>vdw_df_threshold)) 714 if ( ivdw == 2) then 715 quit=1 716 end if 717 #endif 718 end if 719 720 ! Here treat the tolvrs criterion: if density/potential residual (squared) 721 ! is less than input tolvrs, stop steps (exit loop here) 722 if (ttolvrs==1 .and. .not. noquit) then 723 if (ttolwfr==0) then 724 if (res2 < tolvrs) then 725 if (optres==0) then 726 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a)' ) ch10,& 727 ' At SCF step',istep,' vres2 =',res2,' < tolvrs=',tolvrs,' =>converged.' 728 else 729 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a)' ) ch10,& 730 ' At SCF step',istep,' nres2 =',res2,' < tolvrs=',tolvrs,' =>converged.' 731 end if 732 call wrtout([std_out, ab_out], message) 733 quit=1 734 else 735 use_dpfft = res2 < tol5 736 end if 737 else if (ttolwfr==1 .and. residm < tolwfr) then 738 if (res2 < tolvrs) then 739 if (optres==0) then 740 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,e10.2,a,e10.2,a)' ) ch10,& 741 ' At SCF step',istep,' max residual=',residm,' < tolwfr=',tolwfr,' AND vres2 =',res2,& 742 & ' < tolvrs=',tolvrs,' =>converged.' 743 else 744 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,e10.2,a,e10.2,a)' ) ch10,& 745 ' At SCF step',istep,' max residual=',residm,' < tolwfr=',tolwfr,' AND nres2 =',res2,& 746 & ' < tolvrs=',tolvrs,' =>converged.' 747 end if 748 call wrtout([std_out, ab_out], message) 749 quit=1 750 else 751 use_dpfft = res2 < tol5 752 end if 753 end if 754 end if 755 756 if (quit==1.and.noquit) then 757 write(message, '(a,a,a)' ) ch10, & 758 ' SCF cycle will continue as it is in an initialization stage',' (occ. matrix was kept constant)...' 759 call wrtout([std_out, ab_out], message) 760 end if 761 762 end if 763 764 ! Activate FFT in double-precision. 765 if (use_dpfft) then 766 if (fftcore_mixprec == 1) call wrtout(std_out, " Approaching convergence. Activating FFT in double-precision") 767 ii = fftcore_set_mixprec(0) 768 end if 769 770 case (3) 771 ! If wavefunction convergence was not reached (for nstep>0) print a warning and return conv_retcode 772 conv_retcode = 0 773 if(nstep>0) then 774 if (.not. converged()) then 775 conv_retcode = 1 776 777 if(iscf>=1 .or. iscf==-3 .or. iscf == 0)then 778 write(message, '(a,a,a,a,i5,a)' ) ch10,& 779 ' scprqt: WARNING -',ch10,& 780 ' nstep=',nstep,' was not enough SCF cycles to converge;' 781 782 write(std_out,'(6a,i0,3a)')ch10,& 783 "--- !ScfConvergenceWarning",ch10,& 784 "message: |",ch10,& 785 ' nstep ',nstep,' was not enough SCF cycles to converge.',ch10,& 786 "..." 787 !ABI_WARNING_CLASS(message, "ScfConvergenceWarning") 788 else 789 write(message, '(a,a,a,a,i5,a)' ) ch10,& 790 ' scprqt: WARNING -',ch10,& 791 ' nstep=',nstep,' was not enough non-SCF iterations to converge;' 792 793 write(std_out,'(8a)')ch10,& 794 "--- !NscfConvergenceWarning",ch10,& 795 "message: |",ch10,TRIM(indent(message)),ch10,& 796 "..." 797 !ABI_WARNING_CLASS(message, "NScfConvergenceWarning") 798 end if 799 call wrtout([std_out, ab_out], message) 800 801 if (ttolwfr==1 .and. residm > tolwfr) then 802 if (dtset%usewvl == 0) then 803 write(message, '(a,es11.3,a,es11.3,a)' ) & 804 ' maximum residual=',residm,' exceeds tolwfr=',tolwfr,ch10 805 806 write(message2, '(a,es11.3,2a)' ) & 807 ' maximum residual each band. tolwfr= ',tolwfr,ch10,& 808 ' iband, isppol, individual band residuals (max over all k-points):' 809 call wrtout(std_out, message2,'COLL') 810 do isppol = 1, dtset%nsppol 811 do iband = 1, dtset%mband 812 write(message3, '(2i6, es11.3)') iband, isppol, residm_band(iband,isppol) 813 call wrtout(std_out,message3,'COLL') 814 end do 815 end do 816 817 else 818 write(message, '(a,es11.3,a,es11.3,a)' ) & 819 ' maximum grdnorm=',residm,' exceeds tolwfr=',tolwfr,ch10 820 end if 821 822 else if (ttoldff==1) then 823 write(message, '(a,es11.3,a,es11.3,a)' ) & 824 ' maximum force difference=',diffor,' exceeds toldff=',toldff,ch10 825 826 else if (ttolrff==1) then 827 write(message, '(a,es11.3,a,es11.3,a)' ) & 828 ' maximum force difference=',diffor,' exceeds tolrff*maxfor=',tolrff*maxfor,ch10 829 830 else if (ttoldfe==1) then 831 write(message, '(a,es11.3,a,es11.3,a)' ) & 832 ' maximum energy difference=',abs(deltae),' exceeds toldfe=',toldfe,ch10 833 834 else if(ttolvrs==1)then 835 if (optres==0) then 836 write(message, '(a,es11.3,a,es11.3,a)' ) & 837 ' potential residual=',res2,' exceeds tolvrs=',tolvrs,ch10 838 else 839 write(message, '(a,es11.3,a,es11.3,a)' ) & 840 ' density residual=',res2,' exceeds tolvrs=',tolvrs,ch10 841 end if 842 end if 843 call wrtout([std_out, ab_out], message) 844 845 if (prtxml == 1) then 846 write(ab_xml_out, "(A)", advance = "NO") ' <status cvState="Failed"' 847 end if 848 849 else 850 ! Convergence is OK 851 if (prtxml == 1) then 852 write(ab_xml_out, "(A)", advance = "NO") ' <status cvState="Ok"' 853 end if 854 end if ! test for convergence reached or not 855 856 if (prtxml == 1) then 857 if (ttoldfe == 1) then 858 if (ttolwfr==0) then 859 write(ab_xml_out, "(A)") ' stop-criterion="toldfe" />' 860 else 861 write(ab_xml_out, "(A)") ' stop-criterion="toldfe+tolwfr" />' 862 end if 863 else if (ttoldff == 1) then 864 if (ttolwfr==0) then 865 write(ab_xml_out, "(A)") ' stop-criterion="toldff" />' 866 else 867 write(ab_xml_out, "(A)") ' stop-criterion="toldff+tolwfr" />' 868 end if 869 else if (ttolrff == 1) then 870 if (ttolwfr==0) then 871 write(ab_xml_out, "(A)") ' stop-criterion="tolrff" />' 872 else 873 write(ab_xml_out, "(A)") ' stop-criterion="tolrff+tolwfr" />' 874 end if 875 else if (ttolvrs == 1) then 876 if (ttolwfr==0) then 877 write(ab_xml_out, "(A)") ' stop-criterion="tolvrs" />' 878 else 879 write(ab_xml_out, "(A)") ' stop-criterion="tolvrs+tolwfr" />' 880 end if 881 else if (ttolwfr == 1) then 882 write(ab_xml_out, "(A)") ' stop-criterion="tolwfr" />' 883 else 884 write(ab_xml_out, "(A)") ' />' 885 end if 886 end if 887 888 ! If enabled, output a YAML document with the ETOT iterations 889 !if (etot_yaml_doc%stream%length > 0) call etot_yaml_doc%write_and_free(ab_out) 890 end if ! nstep == 0 : no output 891 892 case default 893 write(message, '(a,i0,a)' )' choice = ',choice,' is not an allowed value.' 894 ABI_BUG(message) 895 end select 896 897 ! Additional stuff for the Fock+SCF cycle 898 if (present(fock)) then 899 if (associated(fock)) then 900 fock%fock_common%scf_converged=(quit==1) 901 ! At present, the decision that the Fock loop is converged is not taken here 902 if (.not.fock%fock_common%fock_converged)quit=0 903 end if 904 end if 905 906 ! Additional stuff for the two-component DFT SCF cycle (electrons+positron) 907 if (present(electronpositron)) then 908 if (associated(electronpositron)) then 909 electronpositron%scf_converged=(quit==1) 910 if (dtset%positron<0) then 911 diff_e=abs(etotal-electronpositron%etotal_prev) 912 diff_f=abs(maxfor-electronpositron%maxfor_prev) 913 end if 914 if (choice==1) then 915 ttoldff=0;ttoldfe=0 916 if(abs(dtset%postoldff)>tiny(0.0_dp))ttoldff=1 917 if(abs(dtset%postoldfe)>tiny(0.0_dp))ttoldfe=1 918 if (dtset%positron<0.and.ttoldff+ttoldfe/=1.and.iscf>0) then 919 ABI_ERROR('one and only one of toldff or toldfe must differ from zero !') 920 end if 921 end if 922 if (choice==2) then 923 if (dtset%positron<0.and.istep<=nstep) then 924 if (electronpositron%scf_converged) then 925 if (electronpositron%istep/=electronpositron%nstep) then 926 if ((.not.noquit).and.& 927 & (diff_e<electronpositron%postoldfe.or.diff_f<electronpositron%postoldff).and.& 928 & (mod(electronpositron%calctype,2)==0.or.(dtset%positron>-20.and.dtset%positron/=-2))) then 929 if (diff_e<electronpositron%postoldfe) then 930 write(message, '(2a,i5,5a,es11.3,a,es11.3)' ) ch10, & 931 & ' At SCF step',istep,', the difference between',ch10,& 932 & ' etotal from electronic calculation and etotal from positronic calculation',ch10,& 933 & ' is converged : diff(etot_el-etot_pos)=',diff_e,' < postoldfe=',electronpositron%postoldfe 934 else 935 write(message, '(2a,i5,5a,es11.3,a,es11.3)' ) ch10, & 936 & ' At SCF step',istep,', the difference between',ch10,& 937 & ' max. force from electronic calculation and max. force from positronic calculation',ch10,& 938 & ' is converged : diff(maxfor_el-maxfor_pos)=',diff_f,' < postoldff=',electronpositron%postoldff 939 end if 940 call wrtout([std_out, ab_out], message) 941 else 942 quit=0 943 end if 944 end if 945 end if 946 end if 947 end if 948 if (choice==3) then 949 if (dtset%positron<0.and.nstep>0)then 950 if (diff_e>=electronpositron%postoldfe.and.abs(dtset%postoldfe)>tiny(0.0_dp)) then 951 write(message, '(4a,i5,5a,es11.3,a,es11.3)' ) ch10,& 952 & ' scprqt: WARNING -',ch10,& 953 & ' posnstep=',dtset%posnstep,' was not enough SCF cycles to converge difference between',ch10,& 954 & ' etotal from electronic calculation and etotal from positronic calculation;',ch10,& 955 & ' diff=',diff_e,' exceeds postoldfe=',electronpositron%postoldfe 956 call wrtout([std_out, ab_out], message) 957 end if 958 if (diff_f>=electronpositron%postoldff.and.abs(dtset%postoldff)>tiny(0.0_dp)) then 959 write(message, '(4a,i5,5a,es11.3,a,es11.3)' ) ch10,& 960 & ' scprqt: WARNING -',ch10,& 961 & ' posnstep=',dtset%posnstep,' was not enough SCF cycles to converge difference between',ch10,& 962 & ' max. force from electronic calculation and max. force from positronic calculation;',ch10,& 963 & ' diff=',diff_e,' exceeds postoldff=',electronpositron%postoldff 964 call wrtout([std_out, ab_out], message) 965 end if 966 end if 967 end if 968 end if 969 end if 970 971 call flush_unit(ab_out) 972 973 DBG_EXIT("COLL") 974 975 contains 976 977 logical function converged() 978 979 ! LB-02/01/2017: 980 ! This code avoids evaluation of undefined variables (which could happen in respfn, apparently) 981 logical :: loc_conv 982 loc_conv = .true. 983 if (ttolwfr==1) then 984 if (residm > tolwfr) loc_conv=.false. 985 end if 986 if (ttoldff==1) then 987 if (diffor > toldff) loc_conv=.false. 988 end if 989 if (ttolrff==1) then 990 if (diffor > tolrff*maxfor .and. maxfor > tol16) loc_conv=.false. 991 end if 992 if (ttoldfe==1) then 993 if (abs(deltae) > toldfe) loc_conv=.false. 994 end if 995 if (ttolvrs==1) then 996 if (res2 > tolvrs) loc_conv=.false. 997 end if 998 converged = loc_conv 999 1000 end function converged 1001 1002 end subroutine scprqt
ABINIT/setup1 [ Functions ]
NAME
setup1
FUNCTION
Call near top of main routine to handle setup of various arrays, filenames, checking of input data, etc.
INPUTS
acell(3)=length scales (bohr) ecut_eff=effective energy cutoff (hartree) for planewave basis sphere ecutc_eff=- PAW only - effective energy cutoff (hartree) for the coarse grid natom=number of atoms ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ngfftc(18)=contain all needed information about 3D FFT for the coarse grid nkpt=number of k points nsppol=1 for unpolarized, 2 for spin-polarized ntypat=number of types of atoms response=0 if called by gstate, =1 if called by respfn rprim(3,3)=dimensionless real space primitive translations usepaw= 0 for non paw calculation; =1 for paw calculation
OUTPUT
bantot=total number of bands at all k points gmet(3,3)=metric for reciprocal space inner products (bohr^-2) gprimd(3,3)=dimens. primitive translations for reciprocal space (bohr**-1) gsqcut_eff=Fourier cutoff on G^2 for "large sphere" of radius double gsqcutc_eff=(PAW) Fourier cutoff on G^2 for "large sphere" of radius double for the coarse FFT grid that of the basis sphere--appropriate for charge density rho(G), Hartree potential, and pseudopotentials, corresponding to ecut_eff rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive translations in real space (bohr) ucvol=unit cell volume (bohr^3)
NOTES
SHOULD BE CLEANED !
SOURCE
1044 subroutine setup1(acell,bantot,dtset,ecut_eff,ecutc_eff,gmet,& 1045 & gprimd,gsqcut_eff,gsqcutc_eff,ngfft,ngfftc,nkpt,nsppol,& 1046 & response,rmet,rprim,rprimd,ucvol,usepaw) 1047 1048 !Arguments ------------------------------------ 1049 !scalars 1050 type(dataset_type),intent(in) :: dtset 1051 integer,intent(in) :: nkpt,nsppol 1052 integer,intent(in) :: response,usepaw 1053 integer,intent(out) :: bantot 1054 real(dp),intent(in) :: ecut_eff,ecutc_eff 1055 real(dp),intent(out) :: gsqcut_eff,gsqcutc_eff,ucvol 1056 !arrays 1057 integer,intent(in) :: ngfft(18),ngfftc(18) 1058 real(dp),intent(in) :: acell(3),rprim(3,3) 1059 real(dp),intent(out) :: gmet(3,3),gprimd(3,3),rmet(3,3) 1060 real(dp),intent(out) :: rprimd(3,3) 1061 1062 !Local variables------------------------------- 1063 !scalars 1064 integer :: ikpt,isppol 1065 real(dp) :: boxcut,boxcutc 1066 character(len=500) :: message 1067 !arrays 1068 real(dp) :: k0(3) 1069 1070 ! ************************************************************************ 1071 1072 !Compute bantot 1073 bantot=0 1074 do isppol=1,nsppol 1075 do ikpt=1,nkpt 1076 bantot=bantot+dtset%nband(ikpt+(isppol-1)*nkpt) 1077 end do 1078 end do 1079 1080 if(dtset%nqpt>1.or.dtset%nqpt<0) then 1081 write(message,'(a,i0,5a)')& 1082 'nqpt =',dtset%nqpt,' is not allowed',ch10,& 1083 '(only 0 or 1 are allowed).',ch10,& 1084 'Action: correct your input file.' 1085 ABI_ERROR(message) 1086 end if 1087 1088 ! Compute dimensional primitive translations rprimd 1089 call mkrdim(acell,rprim,rprimd) 1090 1091 ! Obtain dimensional translations in reciprocal space gprimd, 1092 ! metrics and unit cell volume, from rprimd. 1093 ! Also output rprimd, gprimd and ucvol 1094 call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol) 1095 1096 ! Get boxcut for given acell, gmet, ngfft, and ecut_eff 1097 ! (center at 000 for groundstate, center at q for respfn): 1098 ! boxcut=ratio of basis sphere diameter to fft box side 1099 k0(:)=0.0_dp 1100 if(response==1 .and. dtset%nqpt==1)then 1101 k0(:)=dtset%qptn(:) 1102 write(message, '(a)' )' setup1 : take into account q-point for computing boxcut.' 1103 call wrtout([std_out, ab_out], message) 1104 end if 1105 if (usepaw==1) then 1106 write(message,'(2a)') ch10,' Coarse grid specifications (used for wave-functions):' 1107 call wrtout([std_out, ab_out], message) 1108 call getcut(boxcutc,ecutc_eff,gmet,gsqcutc_eff,dtset%iboxcut,ab_out,k0,ngfftc) 1109 write(message,'(2a)') ch10,' Fine grid specifications (used for densities):' 1110 call wrtout([std_out, ab_out], message) 1111 call getcut(boxcut,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,ab_out,k0,ngfft) 1112 else 1113 call getcut(boxcut,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,ab_out,k0,ngfft) 1114 gsqcutc_eff=gsqcut_eff 1115 end if 1116 1117 ! Check that boxcut>=2 if dtset%intxc=1; otherwise dtset%intxc must be set=0 1118 if (boxcut<2.0_dp.and.dtset%intxc==1) then 1119 write(message, '(a,es12.4,a,a,a,a,a)' )& 1120 'boxcut= ',boxcut,' is < 2.0 => intxc must be 0;',ch10,& 1121 'Need larger ngfft to use intxc=1.',ch10,& 1122 'Action: you could increase ngfft, or decrease ecut, or put intxcn=0.' 1123 ABI_ERROR(message) 1124 end if 1125 1126 end subroutine setup1