TABLE OF CONTENTS


ABINIT/m_common [ Modules ]

[ Top ] [ Modules ]

NAME

  m_common

FUNCTION

  This module gathers routines used by higher-level procedures.
  Mainly printing routines.

COPYRIGHT

  Copyright (C) 1998-2018 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 .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_common
29 
30  use defs_basis
31  use defs_abitypes
32  use m_errors
33  use m_abicore
34  use m_exit
35  use m_fock
36  use m_io_tools
37 #if defined DEV_YP_VDWXC
38  use m_xc_vdw
39 #endif
40 
41  use m_fstrings,         only : indent
42  use m_electronpositron, only : electronpositron_type
43 use m_energies,          only : energies_type, energies_eval_eint
44 
45  implicit none
46 
47  private

ABINIT/prteigrs [ Functions ]

[ Top ] [ 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)
  fname_eig=filename for printing of 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
  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)

PARENTS

      clnup1,dfpt_looppert,respfn,scprqt,vtorho

CHILDREN

      wrtout

SOURCE

1102 subroutine prteigrs(eigen,enunit,fermie,fname_eig,iout,iscf,kptns,kptopt,mband,nband,&
1103 &  nkpt,nnsclo_now,nsppol,occ,occopt,option,prteig,prtvol,resid,tolwfr,vxcavg,wtk)
1104 
1105  use m_io_tools,  only : open_file
1106 
1107 !This section has been created automatically by the script Abilint (TD).
1108 !Do not modify the following lines by hand.
1109 #undef ABI_FUNC
1110 #define ABI_FUNC 'prteigrs'
1111 !End of the abilint section
1112 
1113  implicit none
1114 
1115 !Arguments ------------------------------------
1116 !scalars
1117  integer,intent(in) :: enunit,iout,iscf,kptopt,mband,nkpt,nnsclo_now,nsppol
1118  integer,intent(in) :: occopt,option,prteig,prtvol
1119  real(dp),intent(in) :: fermie,tolwfr,vxcavg
1120  character(len=*),intent(in) :: fname_eig
1121 !arrays
1122  integer,intent(in) :: nband(nkpt*nsppol)
1123  real(dp),intent(in) :: eigen(mband*nkpt*nsppol),kptns(3,nkpt)
1124  real(dp),intent(in) :: occ(mband*nkpt*nsppol),resid(mband*nkpt*nsppol)
1125  real(dp),intent(in) :: wtk(nkpt)
1126 
1127 !Local variables-------------------------------
1128 !scalars
1129  integer,parameter :: nkpt_max=50
1130  integer :: band_index,iband,ienunit,ii,ikpt,isppol,nband_index,nband_k,nkpt_eff,tmagnet,tmetal,temp_unit
1131  real(dp) :: convrt,magnet,residk,rhodn,rhoup
1132  character(len=2) :: ibnd_fmt,ikpt_fmt
1133  character(len=7) :: strunit1,strunit2
1134  character(len=39) :: kind_of_output
1135  character(len=500) :: msg
1136 
1137 ! *************************************************************************
1138 
1139  if (enunit<0.or.enunit>2) then
1140    write(msg, '(a,i0)' )' enunit must be 0, 1 or 2. Argument was ',enunit
1141    MSG_BUG(msg)
1142  end if
1143 
1144  if (prteig > 0) then
1145    write(msg, '(a,a)' ) ' prteigrs : about to open file ',TRIM(fname_eig)
1146    call wrtout(iout,msg,'COLL')
1147    if (open_file(fname_eig, msg, newunit=temp_unit, status='unknown', form='formatted') /= 0) then
1148      MSG_ERROR(msg)
1149    end if
1150    rewind(temp_unit) ! always rewind disk file and print latest eigenvalues
1151  end if
1152 
1153  kind_of_output=              ' Eigenvalues                          '
1154  if(option==4) kind_of_output=' Expectation of eigenvalue derivatives'
1155  if(option==5) kind_of_output=' Fan corrections to eigenvalues at T=0'
1156  if(option==6) kind_of_output=' DDW corrections to eigenvalues at T=0'
1157  if(option==7) kind_of_output=' Fan+DDW corrs   to eigenvalues at T=0'
1158 
1159  nkpt_eff=nkpt
1160 
1161 !write(msg,'(a,5i5)')' prtvol,iscf,kptopt,nkpt_eff,nkpt_max ',prtvol,iscf,kptopt,nkpt_eff,nkpt_max
1162 !call wrtout(iout,msg,'COLL')
1163 
1164  if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>nkpt_max)nkpt_eff=nkpt_max
1165  if( (prtvol==0.or.prtvol==1) .and. (iscf/=-2 .or. kptopt>0) .and. nkpt_eff>1 .and. iout==ab_out)nkpt_eff=1
1166 
1167  if(option==1 .or. (option>=3 .and. option<=7))then
1168 
1169    do ienunit=0,1
1170 
1171      if (enunit==1 .and. ienunit==0)cycle
1172      if (enunit==0 .and. ienunit==1)cycle
1173 !  Print eigenvalues in hartree for enunit=0 or 2
1174 !  The definition of two different strings is quite ridiculous. Historical reasons ...
1175 
1176      if (ienunit==0)then
1177        convrt=one
1178        strunit1='hartree'
1179        strunit2='hartree'
1180      end if
1181      if (ienunit==1)then
1182        convrt=Ha_eV
1183        strunit1='   eV  '
1184        strunit2='eV     '
1185      end if
1186 
1187      band_index=0
1188 
1189      if(ienunit==0)then  ! XG20140730 I do not know why this is only done when ienunit==0
1190        tmetal=0
1191        if(option==1 .and. occopt>=3 .and. occopt<=8)tmetal=1
1192        tmagnet=0
1193        if(tmetal==1 .and. nsppol==2)then
1194          tmagnet=1
1195          rhoup = 0._dp
1196          rhodn = 0._dp
1197          nband_index = 1
1198          do isppol=1,nsppol
1199            do ikpt=1,nkpt
1200              nband_k=nband(ikpt+(isppol-1)*nkpt)
1201              do iband=1,nband_k
1202                if(isppol==1) rhoup = rhoup + wtk(ikpt)*occ(nband_index)
1203                if(isppol==2) rhodn = rhodn + wtk(ikpt)*occ(nband_index)
1204                nband_index = nband_index + 1
1205              end do
1206            end do
1207          end do
1208          magnet = abs(rhoup - rhodn)
1209        end if
1210      end if
1211 
1212      if(iscf>=0 .and. (ienunit==0 .or. option==1))then
1213        write(msg, '(3a,f10.5,3a,f10.5)' ) &
1214         ' Fermi (or HOMO) energy (',trim(strunit2),') =',convrt*fermie,'   Average Vxc (',trim(strunit2),')=',convrt*vxcavg
1215        call wrtout(iout,msg,'COLL')
1216        if (prteig > 0) call wrtout(temp_unit,msg,'COLL')
1217      end if
1218 
1219 
1220 !    if( (iscf>=0 .or. iscf==-3) .and. ienunit==0)then     ! This is the most correct
1221      if(iscf>=0 .and. ienunit==0)then ! For historical reasons
1222        if(tmagnet==1)then
1223          write(msg, '(a,es16.8,a,a,es16.8,a,es16.8)' )&
1224 &         ' Magnetization (Bohr magneton)=',magnet,ch10,&
1225 &         ' Total spin up =',rhoup,'   Total spin down =',rhodn
1226          call wrtout(iout,msg,'COLL')
1227          if (prteig > 0) call wrtout(temp_unit,msg,'COLL')
1228        end if
1229      end if
1230 
1231 !    Loop over spins (suppress spin data if nsppol not 2)
1232      do isppol=1,nsppol
1233 
1234        ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
1235        if (nsppol==2.and.isppol==1) then
1236          write(msg, '(4a,'//ikpt_fmt//',2x,a)' ) &
1237 &         trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points, SPIN UP:'
1238        else if (nsppol==2.and.isppol==2) then
1239          write(msg, '(4a,'//ikpt_fmt//',2x,a)' ) &
1240 &         trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points, SPIN DOWN:'
1241        else
1242          write(msg, '(4a,'//ikpt_fmt//',2x,a)' ) &
1243 &         trim(kind_of_output),' (',strunit1,') for nkpt=',nkpt,'k points:'
1244        end if
1245        call wrtout(iout,msg,'COLL')
1246        if (prteig > 0) call wrtout(temp_unit,msg,'COLL')
1247 
1248        if(ienunit==0)then
1249          if(option>=4 .and. option<=7)then
1250            msg = '  (in case of degenerate eigenvalues, averaged derivative)'
1251            call wrtout(iout,msg,'COLL')
1252            if (prteig > 0) call wrtout(temp_unit,msg,'COLL')
1253          end if
1254        end if
1255 
1256        do ikpt=1,nkpt
1257          nband_k=nband(ikpt+(isppol-1)*nkpt)
1258          ikpt_fmt="i4" ; if(nkpt>=10000)ikpt_fmt="i6" ; if(nkpt>=1000000)ikpt_fmt="i9"
1259          ibnd_fmt="i3" ; if(nband_k>=1000)ibnd_fmt="i6" ; if(nband_k>=1000000)ibnd_fmt="i9"
1260          if(ikpt<=nkpt_eff)then
1261            write(msg, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) &
1262 &           ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',&
1263 &           kptns(1:3,ikpt)+tol10,' (reduced coord)'
1264            call wrtout(iout,msg,'COLL')
1265            if (prteig > 0) call wrtout(temp_unit,msg,'COLL')
1266            do ii=0,(nband_k-1)/8
1267 !            write(msg, '(8f15.10)' ) (convrt*eigen(iband+band_index),&
1268              write(msg, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index),&
1269 &             iband=1+ii*8,min(nband_k,8+ii*8))
1270              call wrtout(iout,msg,'COLL')
1271              if (prteig > 0) call wrtout(temp_unit,msg,'COLL')
1272            end do
1273            if(ienunit==0 .and. option==1 .and. occopt>=3 .and. occopt<=8)then
1274              write(msg, '(5x,a,'//ikpt_fmt//')' )  ' occupation numbers for kpt#',ikpt
1275              call wrtout(iout,msg,'COLL')
1276              do ii=0,(nband_k-1)/8
1277                write(msg, '(8(f10.5,1x))' ) (occ(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8))
1278                call wrtout(iout,msg,'COLL')
1279              end do
1280            end if
1281 
1282          else
1283            if(ikpt==nkpt_eff+1)then
1284              write(msg, '(a,a)' )' prteigrs : prtvol=0 or 1, do not print more k-points.',ch10
1285              call wrtout(iout,msg,'COLL')
1286            end if
1287            if (prteig > 0) then
1288              write(msg, '(a,'//ikpt_fmt//',a,'//ibnd_fmt//',a,f9.5,a,3f8.4,a)' ) &
1289 &             ' kpt#',ikpt,', nband=',nband_k,', wtk=',wtk(ikpt)+tol10,', kpt=',&
1290 &             kptns(1:3,ikpt)+tol10,' (reduced coord)'
1291              call wrtout(temp_unit,msg,'COLL')
1292              do ii=0,(nband_k-1)/8
1293                write(msg, '(8(f10.5,1x))' ) (convrt*eigen(iband+band_index),iband=1+ii*8,min(nband_k,8+ii*8))
1294                call wrtout(temp_unit,msg,'COLL')
1295              end do
1296            end if
1297          end if
1298          band_index=band_index+nband_k
1299        end do ! do ikpt=1,nkpt
1300      end do ! do isppol=1,nsppol
1301 
1302    end do ! End loop over Hartree or eV
1303 
1304  else if(option==2)then
1305 
1306    band_index=0
1307    do isppol=1,nsppol
1308 
1309      if(nsppol==2)then
1310        if(isppol==1)write(msg, '(2a)' ) ch10,' SPIN UP channel '
1311        if(isppol==2)write(msg, '(2a)' ) ch10,' SPIN DOWN channel '
1312        call wrtout(iout,msg,'COLL')
1313        if(prteig>0) call wrtout(temp_unit,msg,'COLL')
1314      end if
1315 
1316      do ikpt=1,nkpt
1317        nband_k=nband(ikpt+(isppol-1)*nkpt)
1318        ikpt_fmt="i5" ; if(nkpt>=10000)ikpt_fmt="i7" ; if(nkpt>=1000000)ikpt_fmt="i9"
1319 
1320        if(ikpt<=nkpt_eff)then
1321          write(msg, '(1x,a,'//ikpt_fmt//',a,f9.5,2f9.5,a)' ) &
1322 &         'Non-SCF case, kpt',ikpt,' (',(kptns(ii,ikpt),ii=1,3),'), residuals and eigenvalues='
1323          call wrtout(iout,msg,'COLL')
1324          if (prteig > 0) then
1325            write(msg, '(1x,a,'//ikpt_fmt//',a,f9.5,2f9.5,a)' ) &
1326 &           'Non-SCF case, kpt',ikpt,' eig(',(kptns(ii,ikpt),ii=1,3),') '
1327            call wrtout(temp_unit,msg,'COLL')
1328          end if
1329          do ii=0,(nband_k-1)/8
1330            write(msg, '(1p,8e10.2)' )(resid(iband+band_index),iband=1+8*ii,min(8+8*ii,nband_k))
1331            call wrtout(iout,msg,'COLL')
1332          end do
1333          do ii=0,(nband_k-1)/6
1334            write(msg, '(1p,6e12.4)' )(eigen(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k))
1335            call wrtout(iout,msg,'COLL')
1336            if (prteig > 0) call wrtout(temp_unit,msg,'COLL')
1337          end do
1338        else
1339          if(ikpt==nkpt_eff+1)then
1340            write(msg, '(a,a)' )' prteigrs : prtvol=0 or 1, do not print more k-points.',ch10
1341            call wrtout(iout,msg,'COLL')
1342          end if
1343          if (prteig > 0) then
1344            write(msg, '(1x,a,i5,a,f9.5,2f9.5,a)' )'Non-SCF kpt',ikpt,' eig(',(kptns(ii,ikpt),ii=1,3),') '
1345            call wrtout(temp_unit,msg,'COLL')
1346            do ii=0,(nband_k-1)/6
1347              write(msg, '(1p,6e12.4)' )(eigen(iband+band_index),iband=1+6*ii,min(6+6*ii,nband_k))
1348              call wrtout(temp_unit,msg,'COLL')
1349            end do
1350          end if
1351        end if
1352 
1353        ! MG: I don't understand why we should include the buffer in the output.
1354        ! It's already difficult to make the tests pass for the residuals without the buffer if nband >> nbocc
1355        residk=maxval(resid(band_index+1:band_index+nband_k))
1356        if (residk>tolwfr) then
1357          write(msg, '(1x,a,2i5,a,1p,e13.5)' ) &
1358 &         ' prteigrs : nnsclo,ikpt=',nnsclo_now,ikpt,' max resid (incl. the buffer)=',residk
1359          call wrtout(iout,msg,'COLL')
1360        end if
1361 
1362        band_index=band_index+nband_k
1363      end do
1364    end do
1365    call wrtout(iout," ",'COLL')
1366 
1367  else
1368    write(msg, '(a,i0,a)' )' option = ',option,', is not an allowed value.'
1369    MSG_BUG(msg)
1370  end if
1371 
1372  if (prteig > 0) close (temp_unit)
1373 
1374 end subroutine prteigrs

ABINIT/prtene [ Functions ]

[ Top ] [ 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)

PARENTS

      gstate,scfcv

CHILDREN

      energies_eval_eint,wrtout

SOURCE

1408 subroutine prtene(dtset,energies,iout,usepaw)
1409 
1410 
1411 !This section has been created automatically by the script Abilint (TD).
1412 !Do not modify the following lines by hand.
1413 #undef ABI_FUNC
1414 #define ABI_FUNC 'prtene'
1415 !End of the abilint section
1416 
1417  implicit none
1418 
1419 !Arguments ------------------------------------
1420 !scalars
1421  integer,intent(in) :: iout,usepaw
1422  type(dataset_type),intent(in) :: dtset
1423  type(energies_type),intent(in) :: energies
1424 
1425 !Local variables-------------------------------
1426 ! Do not modify the length of this string
1427 !scalars
1428  integer :: ipositron,mu,optdc
1429  logical :: directE_avail,testdmft
1430  real(dp) :: eent,enevalue,etotal,etotaldc
1431  character(len=22) :: eneName
1432  character(len=500) :: msg
1433 !arrays
1434  character(len=10) :: EPName(1:2)=(/"Positronic","Electronic"/)
1435 
1436 ! *************************************************************************
1437 
1438  directE_avail=(usepaw==0.or.dtset%pawspnorb==0.or.dtset%pawcpxocc==2.or.dtset%kptopt==1.or.dtset%kptopt==2)
1439 
1440 !============= Evaluate some parts of the energy ===========
1441 
1442  optdc=-1;ipositron=merge(0,2,dtset%positron==0)
1443  if (abs(energies%e_ewald)<1.e-15_dp.and.abs(energies%e_hartree)<1.e-15_dp) ipositron=1
1444  call energies_eval_eint(energies,dtset,usepaw,optdc,etotal,etotaldc)
1445 
1446 !Here, treat the case of metals
1447 !In re-smeared case the free energy is defined with tphysel
1448  if(dtset%occopt>=3 .and. dtset%occopt<=8)then
1449    if (abs(dtset%tphysel) < tol10) then
1450      eent=-dtset%tsmear * energies%entropy
1451    else
1452      eent=-dtset%tphysel * energies%entropy
1453    end if
1454  else
1455    eent=zero
1456  end if
1457 ! If DMFT is used and DMFT Entropy is not computed, then do not print
1458 ! non interacting entropy
1459  testdmft=(dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(dtset%upawu(:,1))>=tol8.or.  &
1460 & sum(dtset%jpawu(:,1))>tol8).and.dtset%dmft_entropy==0)
1461  if(testdmft) eent=zero
1462 
1463  etotal   = etotal   + eent
1464  etotaldc = etotaldc + eent
1465 
1466  write(msg,'(a,80a)') ch10,('-',mu=1,80)
1467  call wrtout(iout,msg,'COLL')
1468 
1469 !============= Printing of Etotal by direct scheme ===========
1470 
1471  if (dtset%icoulomb == 1) then
1472    write(eneName, "(A)") "    Ion-ion energy  = "
1473  else
1474    write(eneName, "(A)") "    Ewald energy    = "
1475  end if
1476  enevalue = energies%e_ewald
1477 
1478  if (optdc==0.or.optdc==2) then
1479 
1480    if (directE_avail) then
1481      write(msg, '(2a)' ) ' Components of total free energy (in Hartree) :',ch10
1482      call wrtout(iout,msg,'COLL')
1483      write(msg, '(a,es21.14)' ) '    Kinetic energy  = ',energies%e_kinetic
1484      call wrtout(iout,msg,'COLL')
1485      if (ipositron/=1) then
1486        write(msg, '(3(a,es21.14,a),a,es21.14)' ) &
1487 &       '    Hartree energy  = ',energies%e_hartree,ch10,&
1488 &       '    XC energy       = ',energies%e_xc+energies%e_fock+&
1489 &       energies%e_hybcomp_E0-energies%e_hybcomp_v0+energies%e_hybcomp_v,ch10,&
1490 &       eneName            ,enevalue,ch10,&
1491 &       '    PspCore energy  = ',energies%e_corepsp
1492        call wrtout(iout,msg,'COLL')
1493 #if defined DEV_YP_VDWXC
1494        if ( (dtset%vdw_xc > 0) .and. (dtset%vdw_xc < 10) .and. (xc_vdw_status()) ) then
1495          write(msg, '(a,es21.14)' )'    vdW-DF energy   = ',energies%e_xc_vdw
1496          call wrtout(iout,msg,'COLL')
1497        end if
1498 #endif
1499      end if
1500      write(msg, '(a,es21.14)' ) '    Loc. psp. energy= ',energies%e_localpsp
1501      call wrtout(iout,msg,'COLL')
1502      if (usepaw==0) then
1503        write(msg, '(a,es21.14)' ) '    NL   psp  energy= ',energies%e_nonlocalpsp
1504      else
1505        write(msg, '(a,es21.14)' ) '    Spherical terms = ',energies%e_paw
1506      end if
1507      call wrtout(iout,msg,'COLL')
1508      if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).and.ipositron/=1) then
1509        write(msg, '(a,es21.14)' ) '    Vd Waals DFT-D = ',energies%e_vdw_dftd
1510        call wrtout(iout,msg,'COLL')
1511      end if
1512      if (dtset%nzchempot>=1) then
1513        write(msg, '(a,es21.14)' ) '    Chem. potential = ',energies%e_chempot
1514        call wrtout(iout,msg,'COLL')
1515      end if
1516      if(dtset%occopt>=3.and.dtset%occopt<=8.and.ipositron==0) then
1517        if(.not.testdmft) then
1518          write(msg, '(a,es21.14,a,a,a,es21.14)' ) &
1519 &         '    >>>>> Internal E= ',etotal-eent,ch10,ch10,&
1520 &         '    -kT*entropy     = ',eent
1521          call wrtout(iout,msg,'COLL')
1522        else if (testdmft) then
1523          write(msg, '(a,es21.14,a)' ) &
1524 &         '    >>>>> Internal E= ',etotal-eent,ch10
1525          call wrtout(iout,msg,'COLL')
1526        end if
1527      else if (ipositron/=0) then
1528        if (dtset%occopt>=3.and.dtset%occopt<=8) then
1529          write(msg, '(a,es21.14)' ) '    -kT*entropy     = ',eent
1530          call wrtout(iout,msg,'COLL')
1531        end if
1532        write(msg, '(3a,es21.14,a)' ) &
1533 &       '    >>> ',EPName(ipositron),' E= ',etotal-energies%e0_electronpositron &
1534 &       -energies%e_electronpositron,ch10
1535        call wrtout(iout,msg,'COLL')
1536        write(msg, '(3a,es21.14,2a,es21.14)' ) &
1537 &       '    ',EPName(3-ipositron),' ener.= ',energies%e0_electronpositron,ch10,&
1538 &       '    EP interaction E= '             ,energies%e_electronpositron
1539        call wrtout(iout,msg,'COLL')
1540      end if
1541      if ((dtset%berryopt==4 .or.  dtset%berryopt==6 .or. dtset%berryopt==7 .or.  &
1542 &     dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17) .and.ipositron/=1) then
1543        write(msg, '(a,es21.14)' ) '    Electric energy = ',energies%e_elecfield
1544        call wrtout(iout,msg,'COLL')
1545        write(msg, '(a,es21.14)' ) '    Kohn-Sham energy= ',etotal-energies%e_elecfield
1546        call wrtout(iout,msg,'COLL')
1547      end if
1548      write(msg, '(a,es21.14)' ) '    >>>>>>>>> Etotal= ',etotal
1549      call wrtout(iout,msg,'COLL')
1550 
1551    else
1552      write(msg, '(9a)' ) &
1553 &     ' COMMENT: ',ch10,&
1554 &     '  "Direct" decomposition of total free energy cannot be printed out !!!',ch10,&
1555 &     '  PAW contribution due to spin-orbit coupling cannot be evaluated',ch10,&
1556 &     '  without the knowledge of imaginary part of Rhoij atomic occupancies',ch10,&
1557 &     '  (computed only when pawcpxocc=2).'
1558      call wrtout(iout,msg,'COLL')
1559    end if
1560 
1561  end if
1562 !============= Printing of Etotal by double-counting scheme ===========
1563 
1564  if (optdc>=1) then
1565 
1566    write(msg, '(4a,es21.14)' ) ch10,&
1567 &   ' "Double-counting" decomposition of free energy:',ch10,&
1568 &   '    Band energy     = ',energies%e_eigenvalues
1569    call wrtout(iout,msg,'COLL')
1570    if (ipositron/=1) then
1571      write(msg, '(2(a,es21.14,a),a,es21.14)' ) &
1572 &     eneName            ,enevalue,ch10,&
1573 &     '    PspCore energy  = ',energies%e_corepsp-energies%e_corepspdc,ch10,&
1574 &     '    Dble-C XC-energy= ',-energies%e_hartree+energies%e_xc-energies%e_xcdc+&
1575 &     energies%e_fock-energies%e_fockdc+&
1576 &     energies%e_hybcomp_E0-energies%e_hybcomp_v0
1577      call wrtout(iout,msg,'COLL')
1578    end if
1579    if ((dtset%berryopt==4 .or.  dtset%berryopt==6 .or. dtset%berryopt==7 .or.  &
1580 &   dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17).and.ipositron/=1) then
1581      write(msg, '(a,es21.14)' ) '    Electric field  = ',energies%e_elecfield
1582      call wrtout(iout,msg,'COLL')
1583    end if
1584    if (usepaw==1) then
1585      write(msg, '(a,es21.14)' ) '    Spherical terms = ',energies%e_pawdc
1586      call wrtout(iout,msg,'COLL')
1587    end if
1588    if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).and.ipositron/=1) then
1589      write(msg, '(a,es21.14)' ) '    Vd Waals DFT-D = ',energies%e_vdw_dftd
1590      call wrtout(iout,msg,'COLL')
1591    end if
1592    if (dtset%nzchempot>=1) then
1593      write(msg, '(a,es21.14)' ) '    Chem. potential = ',energies%e_chempot
1594      call wrtout(iout,msg,'COLL')
1595    end if
1596    if(dtset%occopt>=3.and.dtset%occopt<=8.and.ipositron==0) then
1597      if(.not.testdmft) then
1598        write(msg, '(a,es21.14,a,a,a,es21.14)' ) &
1599 &       '    >>>>> Internal E= ',etotaldc-eent,ch10,ch10,&
1600 &       '    -kT*entropy     = ',eent
1601        call wrtout(iout,msg,'COLL')
1602      else if (testdmft) then
1603        write(msg, '(a,es21.14,a)' ) '    >>>>> Internal E= ',etotaldc-eent,ch10
1604        call wrtout(iout,msg,'COLL')
1605      end if
1606    else if (ipositron/=0) then
1607      if (dtset%occopt>=3 .and. dtset%occopt<=8) then
1608        write(msg, '(a,es21.14)' ) '    -kT*entropy     = ',eent
1609        call wrtout(iout,msg,'COLL')
1610      end if
1611      write(msg, '(a,es21.14,4a,es21.14,a)' ) &
1612 &     '    - EP dble-ct En.= ',-energies%edc_electronpositron,ch10,&
1613 &     '    >>> ',EPName(ipositron),' E= ',etotaldc-energies%e0_electronpositron &
1614 &     -energies%e_electronpositron,ch10
1615      call wrtout(iout,msg,'COLL')
1616      write(msg, '(3a,es21.14,2a,es21.14)' ) &
1617 &     '    ',EPName(3-ipositron),' ener.= ',energies%e0_electronpositron,ch10,&
1618 &     '    EP interaction E= '            ,energies%e_electronpositron
1619      call wrtout(iout,msg,'COLL')
1620    end if
1621    write(msg, '(a,es21.14)' ) '    >>>> Etotal (DC)= ',etotaldc
1622    call wrtout(iout,msg,'COLL')
1623  end if
1624 
1625 !======= Additional printing for compatibility  ==========
1626 
1627  if (usepaw==0.and.optdc==0) then
1628    write(msg, '(a,a,a,a,es21.14,a,es18.10)' ) ch10,&
1629 &   ' Other information on the energy :',ch10,&
1630 &   '    Total energy(eV)= ',etotal*Ha_eV,' ; Band energy (Ha)= ',energies%e_eigenvalues
1631    call wrtout(iout,msg,'COLL')
1632  end if
1633 
1634  if ((optdc==0.or.optdc==2).and.(.not.directE_avail)) then
1635    write(msg, '(a,a,es18.10)' ) ch10,' Band energy (Ha)= ',energies%e_eigenvalues
1636    call wrtout(iout,msg,'COLL')
1637  end if
1638 
1639  if (usepaw==1) then
1640    if ((optdc==0.or.optdc==2).and.(directE_avail)) then
1641      write(msg, '(a,a,es21.14)' ) ch10,'  >Total energy in eV           = ',etotal*Ha_eV
1642      call wrtout(iout,msg,'COLL')
1643    end if
1644    if (optdc>=1) then
1645      if (optdc==1) write(msg, '(a,a,es21.14)' ) ch10,&
1646 &     '  >Total DC energy in eV        = ',etotaldc*Ha_eV
1647      if (optdc==2) write(msg, '(a,es21.14)' ) &
1648 &     '  >Total DC energy in eV        = ',etotaldc*Ha_eV
1649      call wrtout(iout,msg,'COLL')
1650    end if
1651  end if
1652 
1653  if( dtset%icoulomb/=1.and.abs(dtset%charge)>tol8) then
1654    write(msg, '(6a)' ) &
1655 &   ch10,' Calculation was performed for a charged system with PBC',&
1656 &   ch10,' You may consider including the monopole correction to the total energy',&
1657 &   ch10,' The correction is to be divided by the dielectric constant'
1658    call wrtout(iout,msg,'COLL')
1659    write(msg, '(a,es21.14)' ) '    Monopole correction (Ha)=',energies%e_monopole
1660    call wrtout(iout,msg,'COLL')
1661    write(msg, '(a,es21.14)' ) '    Monopole correction (eV)=',energies%e_monopole*Ha_eV
1662    call wrtout(iout,msg,'COLL')
1663  end if
1664 
1665 !=============
1666  write(msg,'(a,80a)')('-',mu=1,80)
1667  call wrtout(iout,msg,'COLL')
1668 
1669 end subroutine prtene

ABINIT/scprqt [ Functions ]

[ Top ] [ 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=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

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

ABINIT/setup1 [ Functions ]

[ Top ] [ 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 !

PARENTS

      gstate,nonlinear,respfn

CHILDREN

      getcut,metric,mkrdim,wrtout

SOURCE

 943 subroutine setup1(acell,bantot,dtset,ecut_eff,ecutc_eff,gmet,&
 944 &  gprimd,gsqcut_eff,gsqcutc_eff,natom,ngfft,ngfftc,nkpt,nsppol,&
 945 &  response,rmet,rprim,rprimd,ucvol,usepaw)
 946 
 947  use m_geometry,   only : mkrdim, metric
 948  use m_kg,         only : getcut
 949 
 950 !This section has been created automatically by the script Abilint (TD).
 951 !Do not modify the following lines by hand.
 952 #undef ABI_FUNC
 953 #define ABI_FUNC 'setup1'
 954 !End of the abilint section
 955 
 956  implicit none
 957 
 958 !Arguments ------------------------------------
 959 !scalars
 960  type(dataset_type),intent(in) :: dtset
 961  integer,intent(in) :: natom,nkpt,nsppol
 962  integer,intent(in) :: response,usepaw
 963  integer,intent(out) :: bantot
 964  real(dp),intent(in) :: ecut_eff,ecutc_eff
 965  real(dp),intent(out) :: gsqcut_eff,gsqcutc_eff,ucvol
 966 !arrays
 967  integer,intent(in) :: ngfft(18),ngfftc(18)
 968  real(dp),intent(in) :: acell(3),rprim(3,3)
 969  real(dp),intent(out) :: gmet(3,3),gprimd(3,3),rmet(3,3)
 970  real(dp),intent(out) :: rprimd(3,3)
 971 
 972 !Local variables-------------------------------
 973 !scalars
 974  integer :: iatom,ikpt,isppol
 975  real(dp) :: boxcut,boxcutc
 976  character(len=500) :: message
 977 !arrays
 978  real(dp) :: k0(3)
 979 
 980 ! ************************************************************************
 981 
 982 !Compute bantot
 983  bantot=0
 984  do isppol=1,nsppol
 985    do ikpt=1,nkpt
 986      bantot=bantot+dtset%nband(ikpt+(isppol-1)*nkpt)
 987    end do
 988  end do
 989 
 990  if(dtset%nqpt>1.or.dtset%nqpt<0) then
 991    write(message,'(a,i0,5a)')&
 992 &   '  nqpt =',dtset%nqpt,' is not allowed',ch10,&
 993 &   '  (only 0 or 1 are allowed).',ch10,&
 994 &   '  Action : correct your input file.'
 995    call wrtout(ab_out,message,'COLL')
 996    MSG_ERROR(message)
 997  end if
 998 
 999 !Compute dimensional primitive translations rprimd
1000  call mkrdim(acell,rprim,rprimd)
1001 
1002 !Obtain dimensional translations in reciprocal space gprimd,
1003 !metrics and unit cell volume, from rprimd.
1004 !Also output rprimd, gprimd and ucvol
1005  call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
1006 
1007 !Get boxcut for given acell, gmet, ngfft, and ecut_eff
1008 !(center at 000 for groundstate, center at q for respfn):
1009 !boxcut=ratio of basis sphere diameter to fft box side
1010  k0(:)=0.0_dp
1011  if(response==1 .and. dtset%nqpt==1)then
1012    k0(:)=dtset%qptn(:)
1013    write(message, '(a)' )' setup1 : take into account q-point for computing boxcut.'
1014    call wrtout(ab_out,message,'COLL')
1015    call wrtout(std_out,message,'COLL')
1016  end if
1017  if (usepaw==1) then
1018    write(message,'(2a)') ch10,' Coarse grid specifications (used for wave-functions):'
1019    call wrtout(ab_out,message,'COLL')
1020    call wrtout(std_out,message,'COLL')
1021    call getcut(boxcutc,ecutc_eff,gmet,gsqcutc_eff,dtset%iboxcut,ab_out,k0,ngfftc)
1022    write(message,'(2a)') ch10,' Fine grid specifications (used for densities):'
1023    call wrtout(ab_out,message,'COLL')
1024    call wrtout(std_out,message,'COLL')
1025    call getcut(boxcut,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,ab_out,k0,ngfft)
1026  else
1027    call getcut(boxcut,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,ab_out,k0,ngfft)
1028    gsqcutc_eff=gsqcut_eff
1029  end if
1030 
1031 !Check that boxcut>=2 if dtset%intxc=1; otherwise dtset%intxc must be set=0
1032  if (boxcut<2.0_dp.and.dtset%intxc==1) then
1033    write(message, '(a,es12.4,a,a,a,a,a)' )&
1034 &   '  boxcut=',boxcut,' is < 2.0  => intxc must be 0;',ch10,&
1035 &   '  Need larger ngfft to use intxc=1.',ch10,&
1036 &   '  Action : you could increase ngfft, or decrease ecut, or put intxcn=0.'
1037    MSG_ERROR(message)
1038  end if
1039 
1040 end subroutine setup1