1! 2! Copyright (C) 2004 PWSCF 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! 9!--------------------------------------------------------------- 10SUBROUTINE run_pseudo 11 !--------------------------------------------------------------- 12 ! 13 ! this routine is a driver to a pseudopotential calculation 14 ! with the parameters given in input 15 ! 16 ! 17 USE kinds, ONLY : dp 18 USE radial_grids, ONLY : ndmx 19 USE ld1_parameters, ONLY : nwfsx 20 USE ld1inc, ONLY : enl, lpaw, nlcc, lsd, latt, pawsetup, & 21 nstoaets, grid, nspin, iter, rhos, rhoc, & 22 nwfts, enlts, llts, jjts, iswts, octs, phits, & 23 vxt, enne, vh, vpsloc, file_potscf, beta, tr2, & 24 eps0, file_recon, deld, vpstot, nbeta, ddd, etots, & 25 paw_energy, iswitch, lgipaw_reconstruction, use_paw_as_gipaw 26 USE atomic_paw, ONLY : new_paw_hamiltonian 27 IMPLICIT NONE 28 29 INTEGER :: & 30 ns, & ! counter on pseudowavefunctions 31 n, & ! counter on mesh 32 is ! counter on spin 33 34 real(DP) :: & 35 vnew(ndmx,2) ! the potential 36 37 INTEGER :: & 38 iunps, ios 39 40 LOGICAL :: & 41 conv ! if true convergence reached 42 43 real(DP) :: & 44 dddnew(nwfsx,nwfsx,2), & ! the new D coefficients 45 vd(2*(ndmx+nwfsx+nwfsx)), & ! Vloc and D in one array for mixing 46 vdnew(2*(ndmx+nwfsx+nwfsx)) ! the new vd array 47 INTEGER :: & 48 nerr ! error message 49 50 real(DP), PARAMETER :: thresh=1.e-10_dp 51 INTEGER, PARAMETER :: itmax=200 52 CHARACTER(len=256) :: nomefile 53 54 ! 55 ! initial estimate of the eigenvalues 56 ! 57 DO ns=1,nwfts 58 enlts(ns)=enl(nstoaets(ns)) 59 ENDDO 60 ! 61 ! compute an initial estimate of the potential 62 ! 63 CALL guess_initial_wfc() 64 IF (.not.lpaw) THEN 65 CALL start_potps ( ) 66 ELSE 67 CALL new_paw_hamiltonian (vpstot, ddd, etots,pawsetup, nwfts, & 68 llts, jjts, nspin, iswts, octs, phits, enlts) 69 DO is=1,nspin 70 vpstot(1:grid%mesh,is)=vpstot(1:grid%mesh,is)-pawsetup%psloc(1:grid%mesh) 71 ENDDO 72 CALL vdpack (grid%mesh, ndmx, nbeta, nwfsx, nspin, vpstot, ddd, vd, "PACK") 73 DO is=1,nspin 74 vpstot(1:grid%mesh,is)=vpstot(1:grid%mesh,is)+pawsetup%psloc(1:grid%mesh) 75 ENDDO 76 ENDIF 77 ! 78 ! iterate to self-consistency 79 ! 80 DO iter=1,itmax 81 CALL ascheqps_drv(vpstot, nspin, thresh, .false., nerr) 82 83 IF (.not.lpaw) THEN 84 ! 85 CALL chargeps(rhos,phits,nwfts,llts,jjts,octs,iswts) 86 CALL new_potential(ndmx,grid%mesh,grid,0.0_dp,vxt,lsd,& 87 nlcc,latt,enne,rhoc,rhos,vh,vnew,1) 88 89 DO is=1,nspin 90 vpstot(:,is)=vpstot(:,is)-vpsloc(:) 91 ENDDO 92 93 IF (file_potscf/=' ') THEN 94 IF (iter<10) THEN 95 WRITE(nomefile,'(a,"_",i1)') trim(file_potscf), iter 96 ELSEIF(iter<100) THEN 97 WRITE(nomefile,'(a,"_",i2)') trim(file_potscf), iter 98 ELSEIF(iter<1000) THEN 99 WRITE(nomefile,'(a,"_",i3)') trim(file_potscf), iter 100 ELSE 101 CALL errore('run_pseudo','problem with iteration',1) 102 ENDIF 103 OPEN(unit=18,file=trim(nomefile), status='unknown', & 104 err=100, iostat=ios) 105100 CALL errore('run_pseudo','opening file' // nomefile,abs(ios)) 106 IF (lsd==1) THEN 107 DO n=1,grid%mesh 108 WRITE(18,'(5e20.12)') grid%r(n),vnew(n,1)-vpstot(n,1), & 109 vnew(n,1), vnew(n,2)-vpstot(n,2), vnew(n,2) 110 ENDDO 111 ELSE 112 DO n=1,grid%mesh 113 WRITE(18,'(3e26.15)') grid%r(n),vnew(n,1)-vpstot(n,1), & 114 vnew(n,1) 115 ENDDO 116 ENDIF 117 CLOSE(18) 118 ENDIF 119 120 CALL vpack(grid%mesh,ndmx,nspin,vnew,vpstot,1) 121 CALL dmixp(grid%mesh*nspin,vnew,vpstot,beta,tr2,iter,3,eps0,conv,itmax) 122 CALL vpack(grid%mesh,ndmx,nspin,vnew,vpstot,-1) 123 124 DO is=1,nspin 125 DO n=1,grid%mesh 126 vpstot(n,is)=vpstot(n,is)+vpsloc(n) 127 ENDDO 128 ENDDO 129 CALL newd_at 130 ! 131 ELSE 132 ! 133 CALL new_paw_hamiltonian (vnew, dddnew, etots, & 134 pawsetup, nwfts, llts, jjts, nspin, iswts, octs, phits, enlts) 135 DO is=1,nspin 136 vnew(1:grid%mesh,is)=vnew(1:grid%mesh,is)-pawsetup%psloc(1:grid%mesh) 137 ENDDO 138 CALL vdpack (grid%mesh, ndmx, nbeta, nwfsx, nspin, vnew, dddnew, vdnew, "PACK") 139 CALL dmixp((grid%mesh+nbeta*nbeta)*nspin,vdnew,vd,beta,tr2,iter,3,eps0,conv,itmax) 140 CALL vdpack (grid%mesh, ndmx, nbeta, nwfsx, nspin, vpstot, ddd, vd, "UNDO") 141 DO is=1,nspin 142 vpstot(1:grid%mesh,is)=vpstot(1:grid%mesh,is)+pawsetup%psloc(1:grid%mesh) 143 ENDDO 144 ! 145 ENDIF 146 147 IF (conv) THEN 148 IF (nerr /= 0) THEN 149 IF (iswitch==2) THEN 150 CALL infomsg ('run_pseudo','BEWARE! Errors in PS-KS equations') 151 ELSE 152 CALL errore ('run_pseudo','Errors in PS-KS equation', 1) 153 ENDIF 154 ENDIF 155 GOTO 900 156 ENDIF 157 ENDDO 158 CALL infomsg('run_pseudo','Warning: convergence not achieved') 159 ! 160 ! final calculation with all states 161 ! 162900 CONTINUE 163 164 CALL ascheqps_drv(vpstot, nspin, thresh, .true., nerr) 165 166 IF (.not.lpaw) THEN 167 CALL elsdps ( ) 168 ELSE 169 CALL new_paw_hamiltonian (vnew, dddnew, etots, pawsetup, nwfts, & 170 llts, jjts, nspin, iswts, octs, phits, enlts, paw_energy) 171 CALL elsdps_paw() 172 ENDIF 173 174 IF ( lgipaw_reconstruction.and.(.not.use_paw_as_gipaw) ) & 175 CALL calculate_gipaw_orbitals() 176 177 IF (file_recon/=' ') CALL write_paw_recon ( ) 178 ! 179 ! compute logarithmic derivatives 180 ! 181 IF ( deld > 0.0_dp) CALL lderivps ( ) 182 ! 183 ! compute expansion in partial waves 184 ! 185 IF ( deld > 0.0_dp) CALL partial_wave_expansion ( ) 186 187 RETURN 188END SUBROUTINE run_pseudo 189