1! 2! Copyright (C) 2004-2008 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! 9!--------------------------------------------------------------- 10SUBROUTINE run_test 11 ! 12 ! This routine is a driver to the tests of the pseudopotential 13 !--------------------------------------------------------------- 14 ! 15 use kinds, only : dp 16 use io_global, only : ionode, ionode_id 17 use mp, only : mp_bcast 18 use mp_world, only : world_comm 19 use radial_grids, only : ndmx 20 use ld1_parameters, only : nwfx 21 use ld1inc, only : file_tests, prefix, nconf, rel, etot0, & 22 nbeta, grid, psi, pseudotype, els, zed, & 23 rcut, rcutus, rcutts, rcutusts, etot, etots0, etots, & 24 nwf, ll, oc, el, & 25 nwfts, nnts, llts, jjts, iswts, octs, elts, nstoaets, & 26 nwftsc, nntsc, lltsc, jjtsc, iswtsc, octsc, eltsc,nstoaec, & 27 file_wavefunctions, file_logder, & 28 file_wavefunctionsps, file_logderps, & 29 core_state, use_paw_as_gipaw !EMINE 30 implicit none 31 32 integer & 33 n, & ! counter on wavefunctions 34 n1,& ! counter on mesh points 35 ir,& ! counter on mesh points 36 im,& ! position of the maximum 37 nc ! counter on configurations 38 integer :: & 39 nn_old(nwfx), ll_old(nwfx), nwf_old, isw_old(nwfx), lsd_old 40 real(DP) :: & 41 jj_old(nwfx), oc_old(nwfx), enl_old(nwfx), psi_old(ndmx,2,nwfx) 42 logical :: & 43 core_state_old(nwfx) 44 integer :: ios 45 character(len=1) :: nch 46 real(DP) :: dum 47 48 file_tests = trim(prefix)//'.test' 49 if (ionode) & 50 open(unit=13, file=file_tests, iostat=ios, err=1111, status='unknown') 511111 call mp_bcast(ios, ionode_id, world_comm) 52 call errore('run_test','opening file_tests',abs(ios)) 53 54 do nc=1,nconf 55 if (nconf == 1) then 56 file_wavefunctions = trim(prefix)//'.wfc' 57 file_wavefunctionsps= trim(prefix)//'ps.wfc' 58 file_logder = trim(prefix)//'.dlog' 59 file_logderps = trim(prefix)//'ps.dlog' 60 else 61 if (nc < 10) then 62 write (nch, '(i1)') nc 63 else 64 nch='0' 65 call errore ('run_test', & 66 'results for some configs not written to file',-1) 67 endif 68 file_wavefunctions = trim(prefix)//nch//'.wfc' 69 file_wavefunctionsps= trim(prefix)//nch//'ps.wfc' 70 file_logder = trim(prefix)//nch//'.dlog' 71 file_logderps = trim(prefix)//nch//'ps.dlog' 72 endif 73 nwfts=nwftsc(nc) 74 if (nc>1) call save_ae(nwf_old,nn_old,ll_old,jj_old,enl_old, & 75 oc_old,isw_old, core_state_old,psi_old,lsd_old,-1) 76 !EMINE 77 core_state_old = core_state 78 call set_conf(nc) 79 call all_electron(.true.,nc) 80 if (nc.eq.1) then 81 etot0=etot 82 call save_ae(nwf_old,nn_old,ll_old,jj_old,enl_old,oc_old,isw_old, & 83 core_state_old,psi_old,lsd_old,1) 84 endif 85 ! 86 ! choose the cut-off radius for the initial estimate of the wavefunctions 87 ! find the maximum of the all electron wavefunction 88 ! 89 do n=1,nwfts 90 do n1=1,nbeta 91 if (els(n1).eq.elts(n).and.(rcut(n1).gt.1.e-3_dp.or.& 92 rcutus(n1)>1.e-3_DP)) then 93 rcutts(n)=rcut(n1) 94 rcutusts(n)=rcutus(n1) 95 goto 20 96 endif 97 enddo 98 dum=0.0_dp 99 im=2 100 do ir=1,grid%mesh-1 101 dum=abs(psi(ir+1,1,nstoaets(n))) 102 if(dum.gt.abs(psi(ir,1,nstoaets(n)))) im=ir+1 103 enddo 104 if (pseudotype.lt.3) then 105 rcutts(n)=grid%r(im)*1.1_dp 106 rcutusts(n)=grid%r(im)*1.1_dp 107 if (el(nstoaets(n))=='6S'.or.el(nstoaets(n))=='5S') then 108 rcutts(n)=grid%r(im)*1.2_dp 109 rcutusts(n)=grid%r(im)*1.2_dp 110 endif 111 else 112 if (ll(nstoaets(n)).eq.0) then 113 rcutts(n)=grid%r(im)*1.6_dp 114 rcutusts(n)=grid%r(im)*1.7_dp 115 elseif (ll(nstoaets(n)).eq.1) then 116 rcutts(n)=grid%r(im)*1.6_dp 117 rcutusts(n)=grid%r(im)*1.7_dp 118 if (el(nstoaets(n)).eq.'2P') then 119 rcutts(n)=grid%r(im)*1.7_dp 120 rcutusts(n)=grid%r(im)*1.8_dp 121 endif 122 elseif (ll(nstoaets(n)).eq.2) then 123 rcutts(n)=grid%r(im)*2.0_dp 124 rcutusts(n)=grid%r(im)*2.2_dp 125 if (el(nstoaets(n)).eq.'3D') then 126 rcutts(n)=grid%r(im)*2.5_dp 127 if (zed>28) then 128 rcutusts(n)=grid%r(im)*3.4_dp 129 else 130 rcutusts(n)=grid%r(im)*3.0_dp 131 endif 132 endif 133 endif 134 endif 13520 continue 136! write(6,*) n, rcutts(n), rcutusts(n) 137 enddo 138 ! 139 ! and run the pseudopotential test 140 ! 141 call run_pseudo 142 ! 143 if (nc.eq.1) etots0=etots 144 ! 145 ! print results 146 ! 147 call write_resultsps 148 ! 149 call test_bessel ( ) 150 ! 151 enddo 152 if (ionode) close (unit = 13) 153 154 !EMINE 155 if(use_paw_as_gipaw)core_state = core_state_old 156 157END SUBROUTINE run_test 158 159 160SUBROUTINE save_ae(nwf_old,nn_old,ll_old,jj_old,enl_old,oc_old,isw_old, & 161 core_state_old,psi_old,lsd_old,iflag) 162! 163! This routine saves the all-electron configuration, or copy it in the 164! all-electron variables 165 166use kinds, only : dp 167use ld1_parameters, only : nwfx 168use radial_grids, only : ndmx 169use ld1inc, only : nwf, nn, ll, jj, enl, oc, isw, core_state, psi, lsd 170implicit none 171integer, intent(in) :: iflag 172integer :: & 173 nn_old(nwfx), & ! the main quantum number 174 ll_old(nwfx), & ! the orbital angular momentum 175 nwf_old, & ! the number of wavefunctions 176 lsd_old, & ! if 1 the calculation has spin 177 isw_old(nwfx) ! spin of the wfc. if(.not.lsd) all 1 (default) 178 179real(DP) :: & 180 jj_old(nwfx), & ! the total angular momentum 181 oc_old(nwfx), & ! the occupations of the all-electron atom 182 enl_old(nwfx), & ! the energies of the all-electron atom 183 psi_old(ndmx,2,nwfx) ! the all-electron (dirac) wavefunctions 184 185logical :: & 186 core_state_old(nwfx) 187 188if (iflag==1) then 189 nwf_old=nwf 190 lsd_old=lsd 191 nn_old=nn 192 ll_old=ll 193 jj_old=jj 194 oc_old=oc 195 isw_old=isw 196 enl_old=enl 197 psi_old=psi 198else 199 nwf=nwf_old 200 lsd=lsd_old 201 nn=nn_old 202 ll=ll_old 203 jj=jj_old 204 oc=oc_old 205 isw=isw_old 206 enl=enl_old 207 psi=psi_old 208endif 209 210END SUBROUTINE save_ae 211 212 213SUBROUTINE set_conf(nc) 214! 215! This routine copy the variables of the current configuration in the 216! test variables. If we pass from a non-spin polarized to a spin 217! polarized calculation it splits the occupations of the all-electron states. 218! 219use ld1_parameters, only : nwfx 220use ld1inc, only : nwf, nn, ll, oc, isw, el, enl, psi, nstoaets, nwftsc, & 221 core_state, lsdts, eltsc, iswtsc, nnts, llts, jjts, & 222 octs, elts, iswts, nntsc, lltsc, jjtsc, octsc, & 223 jj, frozen_core, lsd, nwfts, nspin 224implicit none 225integer, intent(in) :: nc 226integer :: n, n1 227 228if (lsdts(nc)==1) then 229 if (frozen_core.and.nc>1) then 230 call occ_spin_tot(nwf,nwfx,el,nn,ll,oc,isw,enl,psi) 231 else 232 call occ_spin(nwf,nwfx,el,nn,ll,oc,isw) 233 endif 234 lsd=1 235 nspin=2 236else 237 lsd=0 238 nspin=1 239endif 240 241do n=1,nwf 242 core_state(n)=.true. 243enddo 244do n=1,nwftsc(nc) 245 nstoaets(n)=0 246 do n1=1,nwf 247 if (lsdts(nc).eq.1) then 248 if (eltsc(n,nc).eq.el(n1) & 249 .and.iswtsc(n,nc).eq.isw(n1)) then 250 nstoaets(n)=n1 251 core_state(n1)=.false. 252 endif 253 else 254 if (eltsc(n,nc).eq.el(n1).and.jjtsc(n,nc).eq.jj(n1)) then 255 nstoaets(n)=n1 256 core_state(n1)=.false. 257 endif 258 endif 259 enddo 260 if (nstoaets(n).eq.0) call errore('set_conf', & 261 'all electron wfc corresponding to pseudo-state ' & 262 & //eltsc(n,nc)//' not found',nc) 263enddo 264 265do n=1,nwfts 266 nnts(n)=nntsc(n,nc) 267 llts(n)=lltsc(n,nc) 268 elts(n)=eltsc(n,nc) 269 jjts(n)=jjtsc(n,nc) 270 iswts(n)=iswtsc(n,nc) 271 octs(n)=octsc(n,nc) 272 oc(nstoaets(n))=octs(n) 273enddo 274 275END SUBROUTINE set_conf 276 277