TABLE OF CONTENTS
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.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG,AF) 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
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=LDA+U: number of SCF steps keeping occ. matrix fixed | usepawu=0 if no LDA+U; 1 if LDA+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) 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) 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.
PARENTS
afterscfloop,dfpt_scfcv,scfcv
CHILDREN
exit_check,flush_unit,prteigrs,wrtout,xc_vdw_trigger
SOURCE
99 #if defined HAVE_CONFIG_H 100 #include "config.h" 101 #endif 102 103 #include "abi_common.h" 104 105 106 subroutine scprqt(choice,cpus,deltae,diffor,dtset,& 107 & eigen,etotal,favg,fcart,fermie,fname_eig,filnam1,initGS,& 108 & iscf,istep,kpt,maxfor,moved_atm_inside,mpi_enreg,& 109 & nband,nkpt,nstep,occ,optres,& 110 & prtfor,prtxml,quit,res2,resid,residm,response,tollist,usepaw,& 111 & vxcavg,wtk,xred,conv_retcode,& 112 & electronpositron, fock) ! optional arguments) 113 114 use defs_basis 115 use defs_abitypes 116 use m_errors 117 use m_profiling_abi 118 use m_exit 119 use m_fock 120 use m_io_tools 121 #if defined DEV_YP_VDWXC 122 use m_xc_vdw 123 #endif 124 125 use m_fstrings, only : indent 126 use m_electronpositron, only : electronpositron_type 127 128 !This section has been created automatically by the script Abilint (TD). 129 !Do not modify the following lines by hand. 130 #undef ABI_FUNC 131 #define ABI_FUNC 'scprqt' 132 use interfaces_14_hidewrite 133 use interfaces_67_common, except_this_one => scprqt 134 !End of the abilint section 135 136 implicit none 137 138 !Arguments ------------------------------------ 139 !scalars 140 integer,intent(in) :: choice,initGS,iscf,istep,moved_atm_inside,nkpt,nstep 141 integer,intent(in) :: optres,prtfor,prtxml,response,usepaw 142 integer,intent(out) :: quit,conv_retcode 143 real(dp),intent(in) :: cpus,deltae,diffor,etotal,fermie,maxfor,res2,residm 144 real(dp),intent(in) :: vxcavg 145 character(len=fnlen),intent(in) :: fname_eig,filnam1 146 type(electronpositron_type),pointer,optional :: electronpositron 147 type(fock_type),pointer,optional :: fock 148 type(MPI_type),intent(in) :: mpi_enreg 149 type(dataset_type),intent(in) :: dtset 150 !arrays 151 integer,intent(in) :: nband(nkpt*dtset%nsppol) 152 real(dp),intent(in) :: eigen(dtset%mband*nkpt*dtset%nsppol),favg(3) 153 real(dp),intent(in) :: fcart(3,dtset%natom),kpt(3,nkpt) 154 real(dp),intent(in) :: occ(dtset%mband*nkpt*dtset%nsppol) 155 real(dp),intent(in) :: resid(dtset%mband*nkpt*dtset%nsppol),tollist(12) 156 real(dp),intent(in) :: wtk(nkpt),xred(3,dtset%natom) 157 158 !Local variables------------------------------- 159 !scalars 160 integer,save :: toldfe_ok,toldff_ok,tolrff_ok,ttoldfe,ttoldff,ttolrff,ttolvrs 161 integer,save :: ttolwfr 162 integer :: iatom,iband,iexit,ikpt,isppol,nband_index,nband_k,openexit,option, ishift 163 integer :: tmagnet 164 #if defined DEV_YP_VDWXC 165 integer :: ivdw 166 #endif 167 real(dp),save :: toldfe,toldff,tolrff,tolvrs,tolwfr,vdw_df_threshold 168 real(dp) :: diff_e,diff_f,magnet,rhodn,rhoup 169 real(dp) :: residm_band(dtset%mband,dtset%nsppol) 170 logical :: noquit 171 character(len=500) :: message, message2, message3 172 character(len=2) :: format_istep 173 character(len=5) :: format_magnet 174 character(len=8) :: colname 175 character(len=1) :: firstchar 176 !arrays 177 real(dp) :: f_tmp(3) 178 179 ! ********************************************************************* 180 181 DBG_ENTER("COLL") 182 183 quit=0; conv_retcode=0 184 185 tmagnet=0 186 if(response==0.and.(iscf>0.or.iscf==-3).and.dtset%nsppol==2.and.dtset%occopt>2)tmagnet=1 187 188 ishift=0 189 residm_band = zero 190 do isppol=1, dtset%nsppol 191 do ikpt=1, nkpt 192 do iband=1, nband(ikpt+(isppol-1)*nkpt) 193 ishift = ishift+1 194 residm_band(iband, isppol) = max (resid(ishift), residm_band(iband, isppol)) 195 end do 196 end do 197 end do 198 199 select case (choice) 200 201 case (1) 202 ! Examine tolerance criteria 203 ! NB: The tests on tolwfr and the presence of tolerances in the SCF case are 204 ! also done at the level of the parser in chkinp. 205 tolwfr=tollist(2) 206 toldff=tollist(3) 207 toldfe=tollist(4) 208 tolvrs=tollist(6) 209 tolrff=tollist(7) 210 vdw_df_threshold=tollist(8) 211 ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0; 212 if(abs(tolwfr)>tiny(zero))ttolwfr=1 213 if(abs(toldff)>tiny(zero))ttoldff=1 214 if(abs(tolrff)>tiny(zero))ttolrff=1 215 if(abs(toldfe)>tiny(zero))ttoldfe=1 216 if(abs(tolvrs)>tiny(zero))ttolvrs=1 217 ! If non-scf calculations, tolwfr must be defined 218 if(ttolwfr /= 1 .and. (iscf<0 .and. iscf/=-3) )then 219 write(message,'(a,a,a,es14.6,a,a)')& 220 & 'when iscf <0 and /= -3, tolwfr must be strictly',ch10,& 221 & 'positive, while it is ',tolwfr,ch10,& 222 & 'Action: change tolwfr in your input file and resubmit the job.' 223 MSG_ERROR(message) 224 end if 225 ! toldff only allowed when prtfor==1 226 ! FIXME: this test should be done on input, not during calculation 227 if((ttoldff == 1 .or. ttolrff == 1) .and. prtfor==0 )then 228 MSG_ERROR('toldff only allowed when prtfor=1!') 229 end if 230 ! If SCF calculations, one and only one of these can differ from zero 231 if(ttolwfr+ttoldff+ttoldfe+ttolvrs+ttolrff /= 1 .and. (iscf>0 .or. iscf==-3))then 232 write(message,'(6a,es14.6,a,es14.6,a,es14.6,a,es14.6,a,a,es14.6,a,a,a)' )& 233 & 'For the SCF case, one and only one of the input tolerance criteria ',ch10,& 234 & 'tolwfr, toldff, tolrff, toldfe or tolvrs ','must differ from zero, while they are',ch10,& 235 & 'tolwfr=',tolwfr,', toldff=',toldff,', tolrff=',tolrff,', toldfe=',toldfe,ch10,& 236 & 'and tolvrs=',tolvrs,' .',ch10,& 237 & 'Action: change your input file and resubmit the job.' 238 MSG_ERROR(message) 239 end if 240 241 if (dtset%usewvl == 1) then 242 write(colname, "(A)") "grdnorm " 243 else 244 write(colname, "(A)") "residm " 245 end if 246 if (nstep>0 .and. (iscf>=0 .or.iscf==-3) .and. dtset%prtstm==0) then 247 if(tmagnet==1)then 248 if (prtfor==0) then 249 if (optres==0) then 250 write(message, '(4a)' ) ch10,& 251 & ' iter Etot(hartree) deltaE(h) ',colname,' vres2 magn' 252 else 253 write(message, '(4a)' ) ch10,& 254 & ' iter Etot(hartree) deltaE(h) ',colname,' nres2 magn' 255 end if 256 else 257 if (optres==0) then 258 write(message, '(4a)' ) ch10,& 259 & ' iter Etot(hartree) deltaE(h) ',colname,' vres2 diffor maxfor magn' 260 else 261 write(message, '(4a)' ) ch10,& 262 & ' iter Etot(hartree) deltaE(h) ',colname,' nres2 diffor maxfor magn' 263 end if 264 end if 265 else 266 if(response==0)then 267 if (prtfor==0) then 268 if (optres==0) then 269 write(message, '(4a)' ) ch10,& 270 & ' iter Etot(hartree) deltaE(h) ', colname, ' vres2' 271 else 272 write(message, '(4a)' ) ch10,& 273 & ' iter Etot(hartree) deltaE(h) ', colname, ' nres2' 274 end if 275 else 276 if (optres==0) then 277 write(message, '(4a)' ) ch10,& 278 & ' iter Etot(hartree) deltaE(h) ',colname,' vres2 diffor maxfor ' 279 else 280 write(message, '(4a)' ) ch10,& 281 & ' iter Etot(hartree) deltaE(h) ',colname,' nres2 diffor maxfor ' 282 end if 283 end if 284 else 285 if (optres==0) then 286 write(message, '(4a)' ) ch10,& 287 & ' iter 2DEtotal(Ha) deltaE(Ha) ', colname, ' vres2' 288 else 289 write(message, '(4a)' ) ch10,& 290 & ' iter 2DEtotal(Ha) deltaE(Ha) ', colname, ' nres2' 291 end if 292 end if 293 end if 294 call wrtout(ab_out,message,'COLL') 295 end if 296 297 case (2) 298 299 300 ! Examine tolerance criteria 301 tolwfr=tollist(2) 302 toldff=tollist(3) 303 toldfe=tollist(4) 304 tolvrs=tollist(6) 305 tolrff=tollist(7) 306 vdw_df_threshold=tollist(8) 307 ttolwfr=0 ; ttoldff=0 ; ttoldfe=0 ; ttolvrs=0; ttolrff=0; 308 if(abs(tolwfr)>tiny(0.0_dp))ttolwfr=1 309 if(abs(toldff)>tiny(0.0_dp))ttoldff=1 310 if(abs(tolrff)>tiny(0.0_dp))ttolrff=1 311 if(abs(toldfe)>tiny(0.0_dp))ttoldfe=1 312 if(abs(tolvrs)>tiny(0.0_dp))ttolvrs=1 313 ! Conduct printing. If extra output follows, then put a blank line into the output here 314 if (dtset%prtvol>=10) then 315 message = ' ' 316 call wrtout(ab_out,message,'COLL') 317 call wrtout(std_out, message,'COLL') 318 end if 319 320 ! Calculate up and down charge and magnetization 321 if(tmagnet==1) then 322 rhoup = zero 323 rhodn = zero 324 nband_index = 1 325 do isppol=1,dtset%nsppol 326 do ikpt=1,nkpt 327 nband_k=nband(ikpt+(isppol-1)*nkpt) 328 do iband=1,nband_k 329 if(isppol==1) rhoup = rhoup + wtk(ikpt)*occ(nband_index) 330 if(isppol==2) rhodn = rhodn + wtk(ikpt)*occ(nband_index) 331 nband_index = nband_index + 1 332 end do 333 end do 334 end do 335 magnet = abs(rhoup - rhodn) 336 end if 337 338 if (prtxml == 1) then 339 write(ab_xml_out, "(A)", advance = "NO") ' <scfcvStep' 340 write(message, "(es22.10)") etotal 341 write(ab_xml_out, "(A,A,A)", advance = "NO") ' eTotal="', trim(message) ,'"' 342 write(message, "(es20.8)") deltae 343 write(ab_xml_out, "(A,A,A)", advance = "NO") ' deltaETotal="', trim(message) ,'"' 344 write(message, "(es20.8)") residm 345 write(ab_xml_out, "(A,A,A)", advance = "NO") ' maxResid="', trim(message) ,'"' 346 write(message, "(es20.8)") res2 347 if (optres == 0) then 348 write(ab_xml_out, "(A,A,A)", advance = "NO") ' potResid="', trim(message) ,'"' 349 else 350 write(ab_xml_out, "(A,A,A)", advance = "NO") ' denResid="', trim(message) ,'"' 351 end if 352 if (tmagnet== 1) then 353 write(message, "(es20.8)") magnet 354 write(ab_xml_out, "(A,A,A)", advance = "NO") ' magn="', trim(message) ,'"' 355 end if 356 if (prtfor == 1) then 357 write(message, "(es20.8)") diffor 358 write(ab_xml_out, "(A,A,A)", advance = "NO") ' deltaForces="', trim(message) ,'"' 359 write(message, "(es20.8)") maxfor 360 write(ab_xml_out, "(A,A,A)", advance = "NO") ' maxForces="', trim(message) ,'"' 361 end if 362 write(ab_xml_out, "(A)") " />" 363 end if 364 365 ! Print total (free) energy (hartree) and other convergence measures 366 if(dtset%prtstm==0)then 367 format_istep='i3' 368 if(istep>99)format_istep='i5' 369 if(istep>9999)format_istep='i7' 370 if(tmagnet==1)then 371 if(magnet<10)then 372 format_magnet='f6.3)' 373 else if(magnet<100)then 374 format_magnet='f6.2)' 375 else 376 format_magnet='f6.1)' 377 end if 378 if (prtfor==0) then 379 write(message, '(a,'//format_istep//',1p,g22.14,3es9.2,0p,'//format_magnet ) & 380 & ' ETOT',istep,etotal,deltae,residm,res2,magnet 381 else 382 write(message, '(a,'//format_istep//',1p,g22.14,3es9.2,es8.1,es9.2,0p,'//format_magnet ) & 383 & ' ETOT',istep,etotal,deltae,residm,res2,diffor,maxfor,magnet 384 end if 385 else 386 firstchar=' ' 387 if (response/=0.and.istep==1) firstchar="-" 388 if (response==0) then 389 if (prtfor==0) then 390 write(message, '(2a,'//format_istep//',1p,g22.14,3es10.3)' ) & 391 & firstchar,'ETOT',istep,etotal,deltae,residm,res2 392 else 393 write(message, '(2a,'//format_istep//',1p,g22.14,5es10.3)' ) & 394 & firstchar,'ETOT',istep,etotal,deltae,residm,res2,diffor,maxfor 395 end if 396 else 397 write(message, '(2a,'//format_istep//',1p,g22.14,1x,3es10.3)' ) & 398 & firstchar,'ETOT',istep,etotal,deltae,residm,res2 399 end if 400 end if 401 call wrtout(ab_out,message,'COLL') 402 403 if(mpi_enreg%paral_pert==1) then 404 call wrtout(std_out, message,'PERS') 405 elseif(mpi_enreg%paral_pert==0) then 406 call wrtout(std_out, message,'COLL') 407 end if 408 409 end if ! dtset%prtstm==0 410 411 ! Print positions/forces every step if dtset%prtvol>=10 and iscf>0 or -3 and GS case 412 if (dtset%prtvol>=10.and.(iscf>=0.or.iscf==-3).and.response==0.and.dtset%prtstm==0) then 413 call wrtout(ab_out," ",'COLL') 414 415 ! Print up and down charge and magnetization 416 if(tmagnet==1) then 417 write(message,'(a,f11.6,a,f11.6,a,f10.6)')& 418 & ' #electrons spin up=',rhoup,', spin down=',rhodn,', magnetization=',magnet 419 call wrtout(ab_out,message,'COLL') 420 call wrtout(std_out, message,'COLL') 421 end if 422 423 ! Moreover, print atomic positions if dtset%ionmov==4, and moved_atm_inside==1 424 if (dtset%ionmov==4 .and. moved_atm_inside==1)then 425 message = ' reduced coordinates :' 426 call wrtout(ab_out,message,'COLL') 427 call wrtout(std_out,message,'COLL') 428 do iatom=1,dtset%natom 429 write(message, '(i5,1x,3es21.11)' ) iatom,xred(:,iatom) 430 call wrtout(ab_out,message,'COLL') 431 call wrtout(std_out,message,'COLL') 432 end do 433 end if 434 435 ! Slightly change favg for printing reasons 436 if (prtfor>0) then 437 f_tmp(:)=favg(:) 438 if(abs(favg(1))<1.0d-13)f_tmp(1)=zero 439 if(abs(favg(2))<1.0d-13)f_tmp(2)=zero 440 if(abs(favg(3))<1.0d-13)f_tmp(3)=zero 441 write(message, '(a,3es10.2)' )' cartesian forces (ha/bohr); non-corrected avg=',f_tmp(:) 442 call wrtout(ab_out,message,'COLL') 443 call wrtout(std_out,message,'COLL') 444 do iatom=1,dtset%natom 445 f_tmp(:)=fcart(:,iatom) 446 if(abs(fcart(1,iatom))<1.0d-13)f_tmp(1)=zero 447 if(abs(fcart(2,iatom))<1.0d-13)f_tmp(2)=zero 448 if(abs(fcart(3,iatom))<1.0d-13)f_tmp(3)=zero 449 write(message, '(i5,1x,3es21.11)' ) iatom,f_tmp(:) 450 call wrtout(ab_out,message,'COLL') 451 call wrtout(std_out,message,'COLL') 452 end do 453 end if 454 455 end if 456 457 ! Print eigenvalues every step if dtset%prtvol>=10 and GS case 458 if (dtset%prtvol>=10 .and. response==0 .and. dtset%tfkinfunc==0 .and. dtset%usewvl==0) then 459 option=1 460 call prteigrs(eigen,dtset%enunit,fermie,fname_eig,ab_out,iscf,kpt,dtset%kptopt,dtset%mband,& 461 & nband,nkpt,dtset%nnsclo,dtset%nsppol,occ,dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk) 462 463 call prteigrs(eigen,dtset%enunit,fermie,fname_eig,std_out,iscf,kpt,dtset%kptopt,dtset%mband,& 464 & nband,nkpt,dtset%nnsclo,dtset%nsppol,occ,dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwfr,vxcavg,wtk) 465 end if 466 467 if(response==0)then 468 write(message, '(a,1p,e15.7,a)' ) ' scprqt: <Vxc>=',vxcavg,' hartree' 469 call wrtout(std_out,message,'COLL') 470 end if 471 472 ! Check whether exiting was required by the user. 473 openexit=1 ; if(dtset%chkexit==0) openexit=0 474 call exit_check(cpus,filnam1,iexit,ab_out,mpi_enreg%comm_cell,openexit) 475 if (iexit/=0) quit=1 476 477 ! In special cases, do not quit even if convergence is reached 478 noquit=((istep<nstep).and.(usepaw==1).and.(dtset%usepawu>0).and.& 479 & (dtset%usedmatpu/=0).and.(istep<=abs(dtset%usedmatpu)).and.& 480 & (dtset%usedmatpu<0.or.initGS==0)) 481 482 ! Additional stuff for electron/positron 483 if (present(electronpositron)) then 484 if (associated(electronpositron)) then 485 if (electronpositron%istep_scf==1) then 486 toldff_ok=0;tolrff_ok=0;toldfe_ok=0 487 end if 488 end if 489 end if 490 491 ! Stopping criteria in the SCF case 492 if(iscf>1 .or. iscf==-3 .or. iscf == 0) then 493 ! Here treat the vdw_df_threshold criterion : if the change of energy is less than 494 ! input vdw_df_threshold, trigger the calculation of vdW interactions 495 ! write(message,'(1x,a,e10.3,1x,a,e10.3,1x,l1,a)') & 496 ! & '[vdW-DF][DEBUG] deltae=',deltae,'vdw_df_threshold=',vdw_df_threshold, & 497 ! & (abs(deltae)<vdw_df_threshold),ch10 498 ! call wrtout(std_out,message,'COLL') 499 #if defined DEV_YP_VDWXC 500 call xc_vdw_trigger( (abs(deltae)<vdw_df_threshold) ) 501 #endif 502 ! Here treat the tolwfr criterion : if maximum residual is less than 503 ! input tolwfr, stop steps (exit loop here) 504 if( ttolwfr==1 .and. residm<tolwfr .and. (.not.noquit)) then 505 if (dtset%usewvl == 0) then 506 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a)' )ch10, & 507 & ' At SCF step',istep,' max residual=',residm,' < tolwfr=',tolwfr,' =>converged.' 508 else 509 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a,a)' )ch10, & 510 & ' At SCF step',istep,' max grdnorm=',residm,' < tolwfr=',tolwfr,' =>converged.' 511 end if 512 call wrtout(ab_out,message,'COLL') 513 call wrtout(std_out,message,'COLL') 514 quit=1 515 end if 516 ! Here treat the toldff criterion : if maximum change of fcart is less than 517 ! input toldff twice consecutively, stop steps (exit loop here) 518 if( ttoldff==1 ) then 519 if( istep==1 )then 520 toldff_ok=0 521 else if (diffor<toldff) then 522 toldff_ok=toldff_ok+1 523 ! add warning for forces which are 0 by symmetry. Also added Matteo check below that the wave 524 ! functions are relatively converged as well 525 if (diffor < tol12) then 526 write (message,'(3a)') ' toldff criterion is satisfied, but your forces are suspiciously low.', ch10,& 527 & ' Check if the forces are 0 by symmetry: in that case you can not use the toldff convergence criterion!' 528 MSG_WARNING(message) 529 if (maxfor < tol16 .and. res2 > tol9) tolrff_ok=0 530 end if 531 else 532 toldff_ok=0 533 end if 534 if(toldff_ok==2 .and. (.not.noquit))then 535 write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3)' ) ch10, & 536 & ' At SCF step',istep,', forces are converged : ',ch10,& 537 & ' for the second time, max diff in force=',diffor,' < toldff=',toldff 538 call wrtout(ab_out,message,'COLL') 539 call wrtout(std_out,message,'COLL') 540 quit=1 541 end if 542 end if 543 ! Here treat the tolrff criterion : if maximum change of fcart is less than 544 ! input tolrff times fcart itself twice consecutively, stop steps (exit loop here) 545 if( ttolrff==1 ) then 546 if( istep==1 )then 547 tolrff_ok=0 548 ! 27/7/2009: added test for absolute value of maxfor, otherwise if it is 0 this never exits the scf loop. 549 else if (diffor<tolrff*maxfor .or. (maxfor < tol16 .and. diffor < tol16)) then 550 tolrff_ok=tolrff_ok+1 551 ! Thu Mar 12 19:01:40 MG: added additional check on res2 to make sure the SCF cycle is close to convergence. 552 ! Needed for structural relaxations otherwise the stress tensor is wrong and the relax algo makes wrong moves. 553 if (maxfor < tol16 .and. res2 > tol9) tolrff_ok=0 554 else 555 tolrff_ok=0 556 end if 557 if(tolrff_ok==2 .and. (.not.noquit))then 558 write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3,a)' ) ch10, & 559 & ' At SCF step',istep,', forces are sufficiently converged : ',ch10,& 560 & ' for the second time, max diff in force=',diffor,& 561 & ' is less than < tolrff=',tolrff, ' times max force' 562 call wrtout(ab_out,message,'COLL') 563 call wrtout(std_out,message,'COLL') 564 quit=1 565 end if 566 end if 567 ! Here treat the toldfe criterion : if the change of energy is less than 568 ! input toldfe twice consecutively, stop steps (exit loop here) 569 if( ttoldfe==1 ) then 570 if( istep==1 )then 571 toldfe_ok=0 572 else if (abs(deltae)<toldfe) then 573 toldfe_ok=toldfe_ok+1 574 else 575 toldfe_ok=0 576 end if 577 if(toldfe_ok==2 .and. (.not.noquit))then 578 write(message, '(a,a,i5,a,a,a,es11.3,a,es11.3)' ) ch10, & 579 & ' At SCF step',istep,', etot is converged : ',ch10,& 580 & ' for the second time, diff in etot=',abs(deltae),' < toldfe=',toldfe 581 call wrtout(ab_out,message,'COLL') 582 call wrtout(std_out,message,'COLL') 583 quit=1 584 end if 585 ! Here treat the vdw_df_threshold criterion for non-SCF vdW-DF 586 ! calculations: If input vdw_df_threshold is lesss than toldfe 587 ! then the vdW-DF is triggered once selfconsistency criteria is 588 ! reached for the first time. 589 ! write(message,'(1x,a,e10.3,1x,a,e10.3,1x,l1,a)') & 590 ! & '[vdW-DF][DEBUG] deltae=',deltae,'vdw_df_threshold=',vdw_df_threshold, & 591 ! & (abs(deltae)<toldfe),ch10 592 ! call wrtout(std_out,message,'COLL') 593 #if defined DEV_YP_VDWXC 594 ivdw = 0 595 if ( toldfe > vdw_df_threshold ) then 596 ivdw = ivdw + 1 597 end if 598 call xc_vdw_trigger((toldfe_ok==1 .and. toldfe>vdw_df_threshold)) 599 if ( ivdw == 2) then 600 quit=1 601 end if 602 #endif 603 end if 604 ! Here treat the tolvrs criterion : if density/potential residual (squared) 605 ! is less than input tolvrs, stop steps (exit loop here) 606 if( ttolvrs==1 .and. res2<tolvrs .and. (.not.noquit)) then 607 if (optres==0) then 608 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a)' ) ch10,& 609 & ' At SCF step',istep,' vres2 =',res2,' < tolvrs=',tolvrs,' =>converged.' 610 else 611 write(message, '(a,a,i5,a,1p,e10.2,a,e10.2,a)' ) ch10,& 612 & ' At SCF step',istep,' nres2 =',res2,' < tolvrs=',tolvrs,' =>converged.' 613 end if 614 call wrtout(ab_out,message,'COLL') 615 call wrtout(std_out,message,'COLL') 616 quit=1 617 end if 618 619 if (quit==1.and.noquit) then 620 write(message, '(a,a,a)' ) ch10, & 621 & ' SCF cycle will continue as it is in an initialization stage',' (occ. matrix was kept constant)...' 622 call wrtout(ab_out,message,'COLL') 623 call wrtout(std_out,message,'COLL') 624 end if 625 626 end if 627 628 case (3) 629 630 ! If wavefunction convergence was not reached (for nstep>0) print a warning and return conv_retcode 631 conv_retcode = 0 632 if(nstep>0) then 633 if (.not. converged()) then 634 conv_retcode = 1 635 636 if(iscf>=1 .or. iscf==-3 .or. iscf == 0)then 637 write(message, '(a,a,a,a,i5,a)' ) ch10,& 638 & ' scprqt: WARNING -',ch10,& 639 & ' nstep=',nstep,' was not enough SCF cycles to converge;' 640 641 write(std_out,'(6a,i0,3a)')ch10,& 642 & "--- !ScfConvergenceWarning",ch10,& 643 & "message: |",ch10,& 644 & ' nstep ',nstep,' was not enough SCF cycles to converge.',ch10,& 645 & "..." 646 !MSG_WARNING_CLASS(message, "ScfConvergenceWarning") 647 else 648 write(message, '(a,a,a,a,i5,a)' ) ch10,& 649 & ' scprqt: WARNING -',ch10,& 650 & ' nstep=',nstep,' was not enough non-SCF iterations to converge;' 651 652 write(std_out,'(8a)')ch10,& 653 & "--- !NscfConvergenceWarning",ch10,& 654 & "message: |",ch10,TRIM(indent(message)),ch10,& 655 & "..." 656 !MSG_WARNING_CLASS(message, "NScfConvergenceWarning") 657 end if 658 659 call wrtout(ab_out,message,'COLL') 660 call wrtout(std_out,message,'COLL') 661 662 if (ttolwfr==1) then 663 if (dtset%usewvl == 0) then 664 write(message, '(a,es11.3,a,es11.3,a)' ) & 665 & ' maximum residual=',residm,' exceeds tolwfr=',tolwfr,ch10 666 667 write(message2, '(a,es11.3,2a)' ) & 668 & ' maximum residual each band. tolwfr= ',tolwfr,ch10,& 669 & ' iband, isppol, individual band residuals (max over all k-points):' 670 call wrtout(std_out, message2,'COLL') 671 do isppol = 1, dtset%nsppol 672 do iband = 1, dtset%mband 673 write(message3, '(2i6, es11.3)') iband, isppol, residm_band(iband,isppol) 674 call wrtout(std_out,message3,'COLL') 675 end do 676 end do 677 678 else 679 write(message, '(a,es11.3,a,es11.3,a)' ) & 680 & ' maximum grdnorm=',residm,' exceeds tolwfr=',tolwfr,ch10 681 end if 682 683 else if (ttoldff==1) then 684 write(message, '(a,es11.3,a,es11.3,a)' ) & 685 & ' maximum force difference=',diffor,' exceeds toldff=',toldff,ch10 686 687 else if (ttolrff==1) then 688 write(message, '(a,es11.3,a,es11.3,a)' ) & 689 & ' maximum force difference=',diffor,' exceeds tolrff*maxfor=',tolrff*maxfor,ch10 690 691 else if (ttoldfe==1) then 692 write(message, '(a,es11.3,a,es11.3,a)' ) & 693 & ' maximum energy difference=',abs(deltae),' exceeds toldfe=',toldfe,ch10 694 695 else if(ttolvrs==1)then 696 if (optres==0) then 697 write(message, '(a,es11.3,a,es11.3,a)' ) & 698 & ' potential residual=',res2,' exceeds tolvrs=',tolvrs,ch10 699 else 700 write(message, '(a,es11.3,a,es11.3,a)' ) & 701 & ' density residual=',res2,' exceeds tolvrs=',tolvrs,ch10 702 end if 703 end if 704 705 call wrtout(ab_out,message,'COLL') 706 call wrtout(std_out,message, 'COLL') 707 708 if (prtxml == 1) then 709 write(ab_xml_out, "(A)", advance = "NO") ' <status cvState="Failed"' 710 end if 711 712 else ! Convergence is OK 713 if (prtxml == 1) then 714 write(ab_xml_out, "(A)", advance = "NO") ' <status cvState="Ok"' 715 end if 716 end if ! test for convergence reached or not 717 718 if (prtxml == 1) then 719 if (ttoldfe == 1) then 720 write(ab_xml_out, "(A)") ' stop-criterion="toldfe" />' 721 else if (ttoldff == 1) then 722 write(ab_xml_out, "(A)") ' stop-criterion="toldff" />' 723 else if (ttolrff == 1) then 724 write(ab_xml_out, "(A)") ' stop-criterion="tolrff" />' 725 else if (ttolvrs == 1) then 726 write(ab_xml_out, "(A)") ' stop-criterion="tolvrs" />' 727 else if (ttolwfr == 1) then 728 write(ab_xml_out, "(A)") ' stop-criterion="tolwfr" />' 729 else 730 write(ab_xml_out, "(A)") ' />' 731 end if 732 end if 733 734 end if ! nstep == 0 : no output 735 736 case default 737 write(message, '(a,i0,a)' )' choice = ',choice,' is not an allowed value.' 738 MSG_BUG(message) 739 end select 740 741 !Additional stuff for the Fock+SCF cycle 742 if (present(fock)) then 743 if (associated(fock)) then 744 fock%fock_common%scf_converged=(quit==1) 745 ! At present, the decision that the Fock loop is converged is not taken here 746 if (.not.fock%fock_common%fock_converged)quit=0 747 end if 748 end if 749 750 !Additional stuff for the two-component DFT SCF cycle (electrons+positron) 751 if (present(electronpositron)) then 752 if (associated(electronpositron)) then 753 electronpositron%scf_converged=(quit==1) 754 if (dtset%positron<0) then 755 diff_e=abs(etotal-electronpositron%etotal_prev) 756 diff_f=abs(maxfor-electronpositron%maxfor_prev) 757 end if 758 if (choice==1) then 759 ttoldff=0;ttoldfe=0 760 if(abs(dtset%postoldff)>tiny(0.0_dp))ttoldff=1 761 if(abs(dtset%postoldfe)>tiny(0.0_dp))ttoldfe=1 762 if (dtset%positron<0.and.ttoldff+ttoldfe/=1.and.iscf>0) then 763 message = 'one and only one of toldff or toldfe must differ from zero !' 764 MSG_ERROR(message) 765 end if 766 end if 767 if (choice==2) then 768 if (dtset%positron<0.and.istep<=nstep) then 769 if (electronpositron%scf_converged) then 770 if (electronpositron%istep/=electronpositron%nstep) then 771 if ((.not.noquit).and.& 772 & (diff_e<electronpositron%postoldfe.or.diff_f<electronpositron%postoldff).and.& 773 & (mod(electronpositron%calctype,2)==0.or.(dtset%positron>-20.and.dtset%positron/=-2))) then 774 if (diff_e<electronpositron%postoldfe) then 775 write(message, '(2a,i5,5a,es11.3,a,es11.3)' ) ch10, & 776 & ' At SCF step',istep,', the difference between',ch10,& 777 & ' etotal from electronic calculation and etotal from positronic calculation',ch10,& 778 & ' is converged : diff(etot_el-etot_pos)=',diff_e,' < postoldfe=',electronpositron%postoldfe 779 else 780 write(message, '(2a,i5,5a,es11.3,a,es11.3)' ) ch10, & 781 & ' At SCF step',istep,', the difference between',ch10,& 782 & ' max. force from electronic calculation and max. force from positronic calculation',ch10,& 783 & ' is converged : diff(maxfor_el-maxfor_pos)=',diff_f,' < postoldff=',electronpositron%postoldff 784 end if 785 call wrtout(ab_out,message,'COLL') 786 call wrtout(std_out,message,'COLL') 787 else 788 quit=0 789 end if 790 end if 791 end if 792 end if 793 end if 794 if (choice==3) then 795 if (dtset%positron<0.and.nstep>0)then 796 if (diff_e>=electronpositron%postoldfe.and.abs(dtset%postoldfe)>tiny(0.0_dp)) then 797 write(message, '(4a,i5,5a,es11.3,a,es11.3)' ) ch10,& 798 & ' scprqt: WARNING -',ch10,& 799 & ' posnstep=',dtset%posnstep,' was not enough SCF cycles to converge difference between',ch10,& 800 & ' etotal from electronic calculation and etotal from positronic calculation;',ch10,& 801 & ' diff=',diff_e,' exceeds postoldfe=',electronpositron%postoldfe 802 call wrtout(ab_out,message,'COLL') 803 call wrtout(std_out,message,'COLL') 804 end if 805 if (diff_f>=electronpositron%postoldff.and.abs(dtset%postoldff)>tiny(0.0_dp)) then 806 write(message, '(4a,i5,5a,es11.3,a,es11.3)' ) ch10,& 807 & ' scprqt: WARNING -',ch10,& 808 & ' posnstep=',dtset%posnstep,' was not enough SCF cycles to converge difference between',ch10,& 809 & ' max. force from electronic calculation and max. force from positronic calculation;',ch10,& 810 & ' diff=',diff_e,' exceeds postoldff=',electronpositron%postoldff 811 call wrtout(ab_out,message,'COLL') 812 call wrtout(std_out,message,'COLL') 813 end if 814 end if 815 end if 816 end if 817 end if 818 819 call flush_unit(ab_out) 820 821 DBG_EXIT("COLL") 822 823 contains 824 825 logical function converged() 826 827 828 !This section has been created automatically by the script Abilint (TD). 829 !Do not modify the following lines by hand. 830 #undef ABI_FUNC 831 #define ABI_FUNC 'converged' 832 !End of the abilint section 833 834 converged = .not.( & 835 & (ttolwfr==1 .and. residm > tolwfr) .or. & 836 & (ttoldff==1 .and. diffor > toldff) .or. & 837 & (ttolrff==1 .and. diffor > tolrff*maxfor .and. maxfor > tol16) .or.& 838 & (ttoldfe==1 .and. abs(deltae) > toldfe) .or. & 839 & (ttolvrs==1 .and. res2 > tolvrs) ) 840 841 end function converged 842 843 end subroutine scprqt