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