1! 2! Copyright (C) 2004i-2010 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8!--------------------------------------------------------------- 9SUBROUTINE scf(ic) 10 !--------------------------------------------------------------- 11 ! 12 ! this routine performs the atomic self-consistent procedure 13 ! self-interaction-correction allowed 14 ! 15 USE kinds, ONLY : dp 16 USE funct, ONLY : dft_is_meta 17 USE radial_grids, ONLY : ndmx 18 USE constants, ONLY: e2 19 USE ld1inc, ONLY : grid, zed, psi, isic, vpot, vh, vxt, rho, iter, & 20 lsd, rel, latt, enne, beta, nspin, tr2, eps0, & 21 nwf, nn, ll, jj, enl, oc, isw, core_state, frozen_core, & 22 tau, vtau, vsic, vsicnew, vhn1, egc, relpert, noscf 23 IMPLICIT NONE 24 25 INTEGER, INTENT(in) :: ic 26 27 LOGICAL:: meta, conv 28 INTEGER:: nerr, nstop, n, i, is, id, nin 29 real(DP) :: vnew(ndmx,2), vtaunew(ndmx), rhoc1(ndmx), ze2 30 INTEGER, PARAMETER :: maxter=200 31 real(DP), PARAMETER :: thresh=1.0e-10_dp 32 ! 33 ! 34 meta = dft_is_meta() 35 ze2 = - zed * e2 36 rhoc1=0.0_dp 37 IF (.not.frozen_core.or.ic==1) psi=0.0_dp 38 DO iter = 1,maxter 39 nerr=0 40 vnew=vpot 41 vtaunew=vtau 42 DO n=1,nwf 43 IF (oc(n) >= 0.0_dp) THEN 44 IF (ic==1.or..not.frozen_core.or..not.core_state(n)) THEN 45 is=isw(n) 46 IF (isic /= 0 .and. iter > 1) vnew(:,is)=vpot(:,is)-vsic(:,n) 47 IF (rel == 0) THEN 48 ! nonrelativistic calculation 49 IF ( meta ) THEN 50 ! Meta-GGA version of lschps 51 CALL lschps_meta (2, zed, thresh, grid, nin, nn(n), ll(n),& 52 enl(n), vnew(1,is), vtaunew, psi(1,1,n), nstop) 53 ELSE 54 ! print *, "Solving nonrelativistic nonmeta equation" 55 CALL ascheq (nn(n),ll(n),enl(n),grid%mesh,grid,& 56 vnew(1,is), & ! potential 57 ze2,thresh,psi(1,1,n),nstop) 58 END IF 59 ELSEIF (rel == 1) THEN 60 ! relativistic scalar calculation 61 IF ( meta ) THEN 62 CALL lschps_meta (1, zed, thresh, grid, nin, nn(n), ll(n),& 63 enl(n), vnew(1,is), vtaunew, psi(1,1,n), nstop) 64 ELSE 65 CALL lschps (1, zed, thresh, grid, nin, nn(n), ll(n),& 66 enl(n), vnew(1,is), psi(1,1,n), nstop) 67 END IF 68 IF (nstop>0.and.oc(n)<1.e-10_DP) nstop=0 69 ELSEIF (rel == 2) THEN 70 CALL dirsol (ndmx,grid%mesh,nn(n),ll(n),jj(n),iter,enl(n), & 71 thresh,grid,psi(1,1,n),vnew(1,is),nstop) 72 ELSE 73 CALL errore('scf','relativistic not programmed',1) 74 ENDIF 75 ! write(6,*) nn(n),ll(n),enl(n) 76 ! if (nstop /= 0) write(6,'(4i6)') iter,nn(n),ll(n),nstop 77 nerr=nerr+nstop 78 ENDIF 79 ELSE 80 enl(n)=0.0_dp 81 psi(:,:,n)=0.0_dp 82 ENDIF 83 ENDDO 84 85 ! 86 ! calculate charge density (spherical approximation) 87 ! 88 rho=0.0_dp 89 90 IF (noscf) GOTO 500 91 DO n=1,nwf 92 rho(1:grid%mesh,isw(n))=rho(1:grid%mesh,isw(n)) + & 93 oc(n)*(psi(1:grid%mesh,1,n)**2+psi(1:grid%mesh,2,n)**2) 94 ENDDO 95 96 97 ! 98 ! calculate kinetc energy density (spherical approximation) 99 ! 100 IF ( meta ) CALL kin_e_density (ndmx, grid%mesh, nwf, & 101 ll, oc, psi, grid%r, grid%r2, grid%dx, tau) 102 ! 103 ! calculate new potential 104 ! 105 106 CALL new_potential ( ndmx, grid%mesh, grid, zed, vxt, & 107 lsd, .false., latt, enne, rhoc1, rho, vh, vnew, 1 ) 108 ! 109 ! calculate SIC correction potential (if present) 110 ! 111 IF (isic /= 0) THEN 112 DO n=1,nwf 113 IF (oc(n) >= 0.0_dp) THEN 114 is=isw(n) 115 CALL sic_correction(n,vhn1,vsicnew,egc) 116 ! 117 ! use simple mixing for SIC correction 118 ! 119 vsic(:,n) = (1.0_dp-beta)*vsic(:,n)+beta*vsicnew(:) 120 ENDIF 121 ENDDO 122 ENDIF 123 ! 124 ! mix old and new potential 125 ! 126 id = 3 127 IF (isic /= 0 .and. relpert) id=1 128 ! 129 CALL vpack(grid%mesh,ndmx,nspin,vnew,vpot,1) 130 CALL dmixp(grid%mesh*nspin,vnew,vpot,beta,tr2,iter,id,eps0,conv,maxter) 131 CALL vpack(grid%mesh,ndmx,nspin,vnew,vpot,-1) 132! write(6,*) iter, eps0 133 ! 134 ! mix old and new metaGGA potential - use simple mixing 135 ! 136 IF ( meta ) vtau(:) = (1.0_dp-beta)*vtaunew(:)+beta*vtau(:) 137 ! 138500 IF (noscf) THEN 139 conv=.true. 140 eps0 = 0.0_DP 141 ENDIF 142 IF (conv) THEN 143 IF (nerr /= 0) CALL infomsg ('scf','warning: at least one error in KS equations') 144 EXIT ! exit cycle 145 ENDIF 146 ENDDO 147 IF ( .not. conv ) CALL infomsg('scf','warning: convergence not achieved') 148 149END SUBROUTINE scf 150