1!
2! Copyright (C) 2020 Quantum ESPRESSO Foundation
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
9MODULE cpmd_module
10  !
11  ! Read and convert a pseudopotential written in the CPMD format, TYPE:
12  ! NORMCONSERVING [ NUMERIC, LOGARITHMIC, CAR ], single radial grid
13  ! to unified pseudopotential format (v.2)
14  !
15  USE upf_kinds, ONLY: dp
16  !
17  IMPLICIT NONE
18  PRIVATE
19  PUBLIC :: read_cpmd
20  !
21  CHARACTER (len=80) title
22  !
23  ! All variables read from CPMD file format
24  !
25  INTEGER :: ixc, pstype = 0
26  REAL(dp) :: alphaxc
27  REAL(dp) :: z, zv
28  ! Grid variables
29  INTEGER :: mesh
30  REAL(dp) :: amesh, rmax, xmin
31  REAL(dp), ALLOCATABLE :: r(:)
32  ! PP variables
33  INTEGER, PARAMETER :: lmaxx=3
34  INTEGER ::lmax, nwfc=0
35  ! Car PP variables
36  REAL(dp) :: alphaloc, alpha(0:lmaxx), a(0:lmaxx), b(0:lmaxx)
37  ! Goedecker PP variables
38  INTEGER, PARAMETER :: ncmax=4, nlmax=3
39  INTEGER :: nc, nl(0:lmaxx)
40  REAL(dp) :: rc, rl(0:lmaxx), c(ncmax), h(0:lmaxx, nlmax*(nlmax+1)/2 )
41  ! Numeric PP variables
42  REAL(dp), ALLOCATABLE :: vnl(:,:)
43  REAL(dp), ALLOCATABLE :: chi(:,:)
44  ! Core correction variables
45  LOGICAL :: nlcc=.false.
46  REAL(dp), ALLOCATABLE :: rho_atc(:)
47  ! Variables used for reading and for checks
48  INTEGER :: maxinfo, info_lines
49  PARAMETER (maxinfo = 100)
50  CHARACTER (len=80), ALLOCATABLE :: info_sect(:)
51  !------------------------------
52  !
53CONTAINS
54  !
55  !     ----------------------------------------------------------
56SUBROUTINE read_cpmd(iunps, upf)
57  !     ----------------------------------------------------------
58  !
59  USE pseudo_types, ONLY : pseudo_upf
60  USE upf_utils, ONLY : matches
61  !
62  IMPLICIT NONE
63  INTEGER, INTENT(in) :: iunps
64  TYPE (pseudo_upf) :: upf
65  !
66  INTEGER :: found = 0, closed = 0, unknown = 0
67  INTEGER :: i, l, dum, ios
68  CHARACTER (len=256) line
69  CHARACTER (len=4) token
70  REAL (dp) :: amesh_, vnl0(0:3)
71  LOGICAL :: grid_read = .false., wfc_read=.false.
72  REAL(dp), EXTERNAL :: upf_erf
73  !
74  info_lines = 0
7510 READ (iunps,'(A)',end=20,err=20) line
76  IF (matches ("&ATOM", trim(line)) ) THEN
77     found = found + 1
78     ! Z
79     READ (iunps,'(a)',end=200,err=200) line
80     l = len_trim(line)
81     i = INDEX(line,'=')
82     READ (line(i+1:l),*) z
83     ! Zv
84     READ (iunps,'(a)',end=200,err=200) line
85     l = len_trim(line)
86     i = INDEX(line,'=')
87     READ (line(i+1:l),*) zv
88     ! XC
89     READ (iunps,'(a)',end=200,err=200) line
90     l = len_trim(line)
91     i = INDEX(line,'=')
92     READ (line(i+1:l),*) ixc, alphaxc
93     ! TYPE
94     READ (iunps,'(a)',end=200,err=200) line
95     IF ( matches("NORMCONSERVING",line) ) THEN
96        IF ( matches("NUMERIC",line) .or. matches("LOGARITHMIC",line)  ) THEN
97           pstype = 1
98        ELSEIF ( matches("CAR",line) ) THEN
99           pstype = 2
100        ELSEIF ( matches("GOEDECKER",line) ) THEN
101           pstype = 3
102        ENDIF
103     ENDIF
104     IF (pstype == 0 ) CALL upf_error('read_cpmd','unknown type: '//line,1)
105  ELSEIF (matches ("&INFO", trim(line)) ) THEN
106     found = found + 1
107     ! read (iunps,'(a)') title
108     ! store info section for later perusal
109     ALLOCATE (info_sect(maxinfo))
110     DO i=1,maxinfo
111        READ (iunps,'(a)',end=20,err=20) title
112        IF (matches ("&END", trim(title)) ) THEN
113           closed = closed + 1
114           GOTO 10
115        ELSE
116           info_sect(i) = trim(title)
117           info_lines = i
118        ENDIF
119     ENDDO
120  ELSEIF (matches ("&POTENTIAL", trim(line)) ) THEN
121     found = found + 1
122     READ (iunps,'(a)') line
123     IF ( pstype == 1 ) THEN
124        !
125        ! NORMCONSERVING NUMERIC
126        !
127        READ (line,*,iostat=ios) mesh, amesh_
128        IF ( ios /= 0) THEN
129           READ (line,*,iostat=ios) mesh
130           amesh_ = -1.0d0
131        ENDIF
132        IF ( .not. grid_read ) ALLOCATE (r(mesh))
133        !
134        ! determine the number of angular momenta
135        !
136        READ (iunps, '(a)') line
137        ios = 1
138        lmax= 4
139        DO WHILE (ios /= 0)
140           lmax = lmax - 1
141           READ(line,*,iostat=ios) r(1),(vnl0(l),l=0,lmax)
142        ENDDO
143        ALLOCATE (vnl(mesh,0:lmax))
144        vnl(1,0:lmax) = vnl0(0:lmax)
145        DO i=2,mesh
146           READ(iunps, *) r(i),(vnl(i,l),l=0,lmax)
147        ENDDO
148        IF ( .not.grid_read ) THEN
149           CALL check_radial_grid ( amesh_, mesh, r, amesh )
150           grid_read = .true.
151        ENDIF
152     ELSEIF ( pstype == 2 ) THEN
153     !
154     ! NORMCONSERVING CAR
155     !
156        READ(iunps, *) alphaloc
157        ! convert r_c's written in file to alpha's: alpha = 1/r_c^2
158        alphaloc = 1.d0/alphaloc**2
159        DO lmax=-1,2
160           READ(iunps, '(A)') line
161           IF (matches ("&END", trim(line)) ) THEN
162              closed = closed + 1
163              exit
164           ENDIF
165           READ(line, *) alpha(lmax+1), a(lmax+1), b(lmax+1)
166           alpha(lmax+1) = 1.d0/alpha(lmax+1)**2
167        ENDDO
168     ELSEIF ( pstype == 3 ) THEN
169     !
170     ! NORMCONSERVING GOEDECKER
171     !
172        c(:) = 0.d0
173        rl(:) = 0.d0
174        nl(:) = 0
175        h(:,:) = 0.d0
176        READ(iunps, *) lmax
177        lmax = lmax - 1
178        IF ( lmax > lmaxx ) &
179          CALL upf_error('read_cpmd','incorrect parameter read',1)
180        READ(iunps, *) rc
181        READ(iunps, '(A)') line
182        READ(line, *) nc
183        IF ( nc > ncmax ) &
184          CALL upf_error('read_cpmd','incorrect parameter read',2)
185        ! I am not sure if it is possible to use nc in the same line
186        ! where it is read. Just in case, better to read twice
187        READ(line, *) dum, (c(i), i=1,nc)
188        DO l=0,lmax+1
189           READ(iunps, '(A)') line
190           IF ( matches ("&END", trim(line)) ) THEN
191              closed = closed + 1
192              exit
193           ENDIF
194           READ(line, *) rl(l), nl(l)
195           IF ( nl(l) > nlmax ) &
196             CALL upf_error('read_cpmd','incorrect parameter read',3)
197           IF ( nl(l) > 0 ) &
198              READ(line, *) rl(l), dum, ( h(l,i), i=1,nl(l)*(nl(l)+1)/2)
199        ENDDO
200     ENDIF
201  ELSEIF (matches ("&WAVEFUNCTION", trim(line)) ) THEN
202     wfc_read=.true.
203     found = found + 1
204     READ (iunps,'(A)') line
205     READ (line,*,iostat=ios) mesh, amesh_
206     IF ( ios /= 0) THEN
207        READ (line,*,iostat=ios) mesh
208        amesh_ = -1.0d0
209     ENDIF
210     IF ( .not. grid_read )  ALLOCATE(r(mesh))
211     ! find number of atomic wavefunctions
212     READ (iunps,'(A)') line
213     DO nwfc = lmax+1,1,-1
214        READ(line,*,iostat=ios) r(1),(vnl0(l),l=0,nwfc-1)
215        IF ( ios == 0 ) exit
216     ENDDO
217     IF ( ios /= 0 ) &
218        CALL upf_error('read_cpmd','at least one atomic wvfct should be present',1)
219     ALLOCATE(chi(mesh,nwfc))
220     chi(1,1:nwfc) = vnl0(0:nwfc-1)
221     DO i=2,mesh
222        READ(iunps, *) r(i),(chi(i,l),l=1,nwfc)
223     ENDDO
224     IF ( .not.grid_read ) THEN
225        CALL check_radial_grid ( amesh_, mesh, r, amesh )
226        grid_read = .true.
227     ENDIF
228  ELSEIF (matches ("&NLCC", trim(line)) ) THEN
229     found = found + 1
230     READ (iunps, '(a)') line
231     READ(iunps, *) mesh
232     nlcc = ( mesh > 0 )
233     IF (nlcc) THEN
234        IF ( .not. matches ("NUMERIC", trim(line)) ) &
235          CALL upf_error('read_cpmd',' only NUMERIC core-correction supported',1)
236        ALLOCATE (rho_atc(mesh))
237        READ(iunps, * ) (r(i), rho_atc(i), i=1,mesh)
238     ENDIF
239  ELSEIF (matches ("&ATDENS", trim(line)) ) THEN
240     ! skip over &ATDENS section, add others here, if there are more.
241     DO WHILE(.not. matches("&END", trim(line)))
242        READ (iunps,'(a)') line
243     ENDDO
244  ELSEIF (matches ("&END", trim(line)) ) THEN
245     closed = closed + 1
246  ELSE
247     PRINT*, 'line ignored: ', line
248     unknown = unknown + 1
249  ENDIF
250  GOTO 10
251
25220 CONTINUE
253
254  IF ( pstype /= 3 ) THEN
255     IF (nlcc .and. found /= 5 .or. .not.nlcc .and. found /= 4) &
256         CALL upf_error('read_cpmd','some &FIELD card missing',found)
257  ELSE
258     IF (found /= 3) &
259         CALL upf_error('read_cpmd','some &FIELD card missing',found)
260  ENDIF
261  IF (closed /= found) &
262       CALL upf_error('read_cpmd','some &END card missing',closed)
263  IF (unknown /= 0 ) PRINT '("WARNING: ",i3," cards not read")', unknown
264  !
265  IF ( .not. grid_read ) THEN
266     xmin = -7.0d0
267     amesh=0.0125d0
268     rmax =100.0d0
269     PRINT '("A radial grid must be provided. We use the following one:")'
270     PRINT '("r_i = e^{xmin+(i-1)*dx}/Z, i=1,mesh, with parameters:")'
271     PRINT '("Z=",f6.2,", xmin=",f6.2," dx=",f8.4," rmax=",f6.1)', &
272           z, xmin, amesh, rmax
273     mesh = 1 + (log(z*rmax)-xmin)/amesh
274     mesh = (mesh/2)*2+1 ! mesh is odd (for historical reasons?)
275     ALLOCATE (r(mesh))
276     DO i=1, mesh
277        r(i) = exp (xmin+(i-1)*amesh)/z
278     ENDDO
279     PRINT '(I4," grid points, rmax=",f8.4)', mesh, r(mesh)
280     grid_read = .true.
281  ENDIF
282  rmax = r(mesh)
283  xmin = log(z*r(1))
284  !
285  IF ( .not. wfc_read ) PRINT '("Notice: atomic wfcs not found")'
286  !
287  IF ( pstype == 2 ) THEN
288     ALLOCATE (vnl(mesh,0:lmax))
289     DO l=0, lmax
290        DO i=1, mesh
291           vnl(i,l) = ( a(l) + b(l)*r(i)**2 ) * exp (-alpha(l)*r(i)**2) - &
292                       zv * upf_erf (sqrt(alphaloc)*r(i))/r(i)
293        ENDDO
294     ENDDO
295  ENDIF
296
297  CALL convert_cpmd(upf)
298  RETURN
299200 CALL upf_error('read_cpmd','error in reading file',1)
300
301END SUBROUTINE read_cpmd
302
303!     ----------------------------------------------------------
304SUBROUTINE convert_cpmd(upf)
305  !     ----------------------------------------------------------
306  !
307  USE pseudo_types, ONLY : pseudo_upf
308  USE upf_const, ONLY : e2
309  !
310  IMPLICIT NONE
311  !
312  TYPE(pseudo_upf) :: upf
313  !
314  REAL(dp), ALLOCATABLE :: aux(:)
315  REAL(dp) :: x, x2, vll, rcloc, fac
316  REAL(dp), EXTERNAL :: upf_erf
317  CHARACTER (len=20):: dft
318  CHARACTER (len=2):: label
319  CHARACTER (len=1):: spdf(0:3) = ['S','P','D','F']
320  CHARACTER (len=2), EXTERNAL :: atom_name
321  INTEGER :: lloc, my_lmax, l, i, j, ij, ir, iv, jv
322  !
323  ! NOTE: many CPMD pseudopotentials created with the 'Hamann' code
324  ! from Juerg Hutter's homepage have additional (bogus) entries for
325  ! pseudo-potential and wavefunction. In the 'report' they have
326  ! the same rc and energy eigenvalue than the previous angular momentum.
327  ! we need to be able to ignore that part or the resulting UPF file
328  ! will be useless. so we first print the info section and ask
329  ! for the LMAX to really use. AK 2005/03/30.
330  !
331  DO i=1,info_lines
332     PRINT '(A)', info_sect(i)
333  ENDDO
334  IF ( pstype == 3 ) THEN
335     ! not actually used, except by write_upf to write a meaningful message
336     lloc = -3
337     rcloc=0.0
338  ELSE
339     WRITE(*,'("max L to use ( <= ",I1," ) > ")', advance="NO") lmax
340     READ (5,*) my_lmax
341     IF ((my_lmax <= lmax) .and. (my_lmax >= 0)) lmax = my_lmax
342     WRITE(*,'("local L ( <= ",I1," ), Rc for local pot (au) > ")', advance="NO") lmax
343     READ (5,*) lloc, rcloc
344  ENDIF
345  !
346  IF ( pstype == 3 ) THEN
347     upf%generated= "Generated in analytical, separable form"
348     upf%author   = "Goedecker/Hartwigsen/Hutter/Teter"
349     upf%date     = "Phys.Rev.B58, 3641 (1998); B54, 1703 (1996)"
350  ELSE
351     upf%generated= "Generated using unknown code"
352     upf%author   = "unknown"
353     upf%date     = "unknown"
354  ENDIF
355  upf%nv       = "2.0.1"
356  upf%comment  = "Info: automatically converted from CPMD format"
357  upf%psd      = atom_name ( nint(z) )
358  ! reasonable assumption
359  IF (z > 18) THEN
360     upf%rel = 'no'
361  ELSE
362     upf%rel = 'scalar'
363  ENDIF
364  IF ( pstype == 3 ) THEN
365     upf%typ = 'NC'
366  ELSE
367     upf%typ = 'SL'
368  ENDIF
369  upf%tvanp = .false.
370  upf%tpawp = .false.
371  upf%tcoulombp=.false.
372  upf%nlcc = nlcc
373  !
374  IF (ixc==900) THEN
375     PRINT '("Pade approx. not implemented! assuming Perdew-Zunger LDA")'
376     upf%dft='SLA-PZ-NOGX-NOGC'
377  ELSEIF (ixc==1100) THEN
378     upf%dft='SLA-PZ-NOGX-NOGC'
379  ELSEIF (ixc==1111) THEN
380     upf%dft='SLA-PZ-B86-P88'
381  ELSEIF (ixc==1134 .or. ixc==1434) THEN
382     upf%dft='SLA-PW-PBX-PBC'
383  ELSEIF (ixc==1134) THEN
384     upf%dft='revPBE'
385  ELSEIF (ixc==1197) THEN
386     upf%dft='PBESOL'
387  ELSEIF (ixc==1312) THEN
388     upf%dft='BLYP'
389  ELSEIF (ixc==362) THEN
390     upf%dft='OLYP'
391  ELSEIF (ixc==1372) THEN
392     upf%dft='XLYP'
393  ELSEIF (ixc==55) THEN
394     upf%dft='HCTH'
395  ELSE
396     WRITE(*,'("Unknown DFT ixc=",i4,". Please provide a DFT name > ")', advance="NO") ixc
397     READ *, upf%dft
398  ENDIF
399  PRINT '("Assuming DFT: ",A," . Please check this is what you want")', &
400          trim(upf%dft)
401  !
402  upf%zp = zv
403  upf%etotps =0.0d0
404  upf%ecutrho=0.0d0
405  upf%ecutwfc=0.0d0
406  IF ( lmax == lloc) THEN
407     upf%lmax = lmax-1
408  ELSE
409     upf%lmax = lmax
410  ENDIF
411  upf%lloc = lloc
412  upf%lmax_rho = 0
413  upf%nwfc = nwfc
414  !
415  ALLOCATE( upf%els(upf%nwfc) )
416  ALLOCATE( upf%oc(upf%nwfc) )
417  ALLOCATE( upf%epseu(upf%nwfc) )
418  ALLOCATE( upf%lchi(upf%nwfc) )
419  ALLOCATE( upf%nchi(upf%nwfc) )
420  ALLOCATE( upf%rcut_chi (upf%nwfc) )
421  ALLOCATE( upf%rcutus_chi(upf%nwfc) )
422
423  DO i=1, upf%nwfc
42410   WRITE(*,'("Wavefunction # ",i1,": label (e.g. 4s), occupancy > ")', advance="NO") i
425     READ (5,*) label, upf%oc(i)
426     READ (label(1:1),*, err=10) l
427     upf%els(i)  = label
428     upf%nchi(i)  = l
429     IF ( label(2:2) == 's' .or. label(2:2) == 'S') THEN
430        l=0
431     ELSEIF ( label(2:2) == 'p' .or. label(2:2) == 'P') THEN
432        l=1
433     ELSEIF ( label(2:2) == 'd' .or. label(2:2) == 'D') THEN
434        l=2
435     ELSEIF ( label(2:2) == 'f' .or. label(2:2) == 'F') THEN
436        l=3
437     ELSE
438        l=i-1
439     ENDIF
440     upf%lchi(i)  = l
441     upf%rcut_chi(i)  = 0.0d0
442     upf%rcutus_chi(i)= 0.0d0
443     upf%epseu(i) = 0.0d0
444  ENDDO
445
446  upf%mesh = mesh
447  upf%dx   = amesh
448  upf%rmax = rmax
449  upf%xmin = xmin
450  upf%zmesh= z
451  ALLOCATE(upf%rab(upf%mesh))
452  ALLOCATE(upf%r(upf%mesh))
453  upf%r(:) = r(1:upf%mesh)
454  upf%rab(:)=upf%r(:)*amesh
455
456  ALLOCATE (upf%rho_atc(upf%mesh))
457  IF (upf%nlcc) upf%rho_atc(:) = rho_atc(1:upf%mesh)
458
459  upf%rcloc = rcloc
460  ALLOCATE (upf%vloc(upf%mesh))
461  !
462  ! the factor e2=2 converts from Hartree to Rydberg
463  !
464  IF ( upf%typ == "SL" ) THEN
465     upf%vloc(:) = vnl(1:upf%mesh,upf%lloc)*e2
466     ALLOCATE(upf%vnl(upf%mesh,0:upf%lmax,1))
467     upf%vnl(:,:,1) = vnl(1:upf%mesh,0:upf%lmax)*e2
468     upf%nbeta= lmax
469  ELSE
470     DO i=1,upf%mesh
471        x = upf%r(i)/rc
472        x2=x**2
473        upf%vloc(i) = e2 * ( -upf%zp*upf_erf(x/sqrt(2.d0))/upf%r(i) + &
474             exp ( -0.5d0*x2 ) * (c(1) + x2*( c(2) + x2*( c(3) + x2*c(4) ) ) ) )
475     ENDDO
476     upf%nbeta=0
477     DO l=0,upf%lmax
478        upf%nbeta = upf%nbeta + nl(l)
479     ENDDO
480  ENDIF
481
482  IF (upf%nbeta > 0) THEN
483     ALLOCATE(upf%els_beta(upf%nbeta) )
484     ALLOCATE(upf%lll(upf%nbeta))
485     ALLOCATE(upf%kbeta(upf%nbeta))
486     IF ( pstype == 3 ) THEN
487        upf%kbeta(:) = upf%mesh
488     ELSE
489        iv=0  ! counter on beta functions
490        DO i=1,upf%nwfc
491           l=upf%lchi(i)
492           IF (l/=lloc) THEN
493              iv=iv+1
494              upf%kbeta(iv)=upf%mesh
495              DO ir = upf%mesh,1,-1
496                 IF ( abs ( vnl(ir,l) - vnl(ir,lloc) ) > 1.0E-6 ) THEN
497                    ! include points up to the last with nonzero value
498                    upf%kbeta(iv)=ir+1
499                    exit
500                 ENDIF
501              ENDDO
502           ENDIF
503        ENDDO
504        ! the number of points used in the evaluation of integrals
505        ! should be even (for simpson integration)
506        DO i=1,upf%nbeta
507           IF ( mod (upf%kbeta(i),2) == 0 ) upf%kbeta(i)=upf%kbeta(i)+1
508           upf%kbeta(i)=MIN(upf%mesh,upf%kbeta(i))
509        ENDDO
510        upf%kkbeta = maxval(upf%kbeta(:))
511     ENDIF
512     ALLOCATE(upf%beta(upf%mesh,upf%nbeta))
513     ALLOCATE(upf%dion(upf%nbeta,upf%nbeta))
514     upf%beta(:,:) =0.d0
515     upf%dion(:,:) =0.d0
516     ALLOCATE(upf%rcut  (upf%nbeta))
517     ALLOCATE(upf%rcutus(upf%nbeta))
518
519     IF ( pstype == 3 ) THEN
520        iv=0  ! counter on beta functions
521        DO l=0,upf%lmax
522           ij = 0
523           DO i=1, nl(l)
524              iv = iv+1
525              upf%lll(iv)=l
526              WRITE (upf%els_beta(iv), '(I1,A1)' ) i, spdf(l)
527              DO j=i, nl(l)
528                 jv = iv+j-i
529                 ij=ij+1
530                 upf%dion(iv,jv) = h(l,ij)/e2
531                 IF ( j > i ) upf%dion(jv,iv) = upf%dion(iv,jv)
532              ENDDO
533              fac= sqrt(2d0*rl(l)) / ( rl(l)**(l+2*i) * sqrt(mygamma(l+2*i)) )
534              DO ir=1,upf%mesh
535                 x2 = (upf%r(ir)/rl(l))**2
536                 upf%beta(ir,iv) = upf%r(ir)**(l+2*(i-1)) * &
537                                     exp ( -0.5d0*x2 ) * fac * e2
538                 ! ...remember: the beta functions in the UPF format
539                 ! ...have to be multiplied by a factor r !!!
540                 upf%beta(ir,iv) = upf%beta(ir,iv)*upf%r(ir)
541                 !
542              ENDDO
543              ! look for index kbeta such that v(i)=0 if i>kbeta
544              DO ir=upf%mesh,1,-1
545                 IF ( abs(upf%beta(ir,iv)) > 1.D-12 ) exit
546              ENDDO
547              IF ( ir < 2 ) THEN
548                 CALL upf_error('cpmd2upf','zero beta function?!?',iv)
549              ELSEIF ( mod(ir,2) /= 0 ) THEN
550                 ! even index
551                 upf%kbeta(iv) = ir
552              ELSEIF ( ir < upf%mesh .and. mod(ir,2) == 0 ) THEN
553                 ! odd index
554                 upf%kbeta(iv) = ir+1
555              ELSE
556                 upf%kbeta(iv) = upf%mesh
557              ENDIF
558              ! not really the same thing as rc in PP generation
559              upf%rcut  (iv) = upf%r(upf%kbeta(iv))
560              upf%rcutus(iv) = 0.0
561           ENDDO
562        ENDDO
563        upf%kkbeta = maxval(upf%kbeta(:))
564     ELSE
565        ALLOCATE(aux(upf%kkbeta))
566        iv=0  ! counter on beta functions
567        DO i=1,upf%nwfc
568           l=upf%lchi(i)
569           IF (l/=lloc) THEN
570              iv=iv+1
571              upf%lll(iv)=l
572              upf%els_beta(iv)=upf%els(i)
573              DO ir=1,upf%kbeta(iv)
574                 ! the factor e2 converts from Hartree to Rydberg
575                 upf%beta(ir,iv) = e2 * chi(ir,l+1) * &
576                      ( vnl(ir,l) - vnl(ir,lloc) )
577                 aux(ir) = chi(ir,l+1) * upf%beta(ir,iv)
578              ENDDO
579              upf%rcut  (iv) = upf%r(upf%kbeta(iv))
580              upf%rcutus(iv) = 0.0
581              CALL simpson(upf%kbeta(iv),aux,upf%rab,vll)
582              upf%dion(iv,iv) = 1.0d0/vll
583           ENDIF
584        ENDDO
585        DEALLOCATE(aux)
586     ENDIF
587  ELSE
588     ! prevents funny errors when writing file
589     ALLOCATE(upf%dion(upf%nbeta,upf%nbeta))
590  ENDIF
591
592  ALLOCATE (upf%chi(upf%mesh,upf%nwfc))
593  upf%chi(:,:) = chi(1:upf%mesh,1:upf%nwfc)
594
595  ALLOCATE (upf%rho_at(upf%mesh))
596  upf%rho_at(:) = 0.d0
597  DO i=1,upf%nwfc
598     upf%rho_at(:) = upf%rho_at(:) + upf%oc(i) * upf%chi(:,i) ** 2
599  ENDDO
600
601  !     ----------------------------------------------------------
602  WRITE (6,'(a)') 'Pseudopotential successfully converted'
603  !     ----------------------------------------------------------
604
605  RETURN
606END SUBROUTINE convert_cpmd
607!
608! ------------------------------------------------------------------
609SUBROUTINE check_radial_grid ( amesh_, mesh, r, amesh )
610! ------------------------------------------------------------------
611!
612   IMPLICIT NONE
613   INTEGER, INTENT (in) :: mesh
614   REAL(dp), INTENT (in) :: amesh_, r(mesh)
615   REAL(dp), INTENT (out) :: amesh
616   INTEGER :: i
617   !
618   ! get amesh if not available directly, check its value otherwise
619   PRINT  "('Radial grid r(i) has ',i4,' points')", mesh
620   PRINT  "('Assuming log radial grid: r(i)=exp[(i-1)*amesh]*r(1), with:')"
621   IF (amesh_ < 0.0d0) THEN
622      amesh = log (r(mesh)/r(1))/(mesh-1)
623      PRINT  "('amesh = log (r(mesh)/r(1))/(mesh-1) = ',f10.6)",amesh
624   ELSE
625   ! not clear whether the value of amesh read from file
626   ! matches the above definition, or if it is exp(amesh) ...
627      amesh = log (r(mesh)/r(1))/(mesh-1)
628      IF ( abs ( amesh - amesh_ ) > 1.0d-5 ) THEN
629         IF ( abs ( amesh - exp(amesh_) ) < 1.0d-5 ) THEN
630            amesh = log(amesh_)
631            PRINT  "('amesh = log (value read from file) = ',f10.6)",amesh
632         ELSE
633            CALL upf_error ('cpmd2upf', 'unknown real-space grid',2)
634         ENDIF
635      ELSE
636         amesh = amesh_
637         PRINT  "('amesh = value read from file = ',f10.6)",amesh
638      ENDIF
639   ENDIF
640   ! check if the grid is what we expect
641   DO i=2,mesh
642      IF ( abs(r(i) - exp((i-1)*amesh)*r(1)) > 1.0d-5) THEN
643         PRINT  "('grid point ',i4,': found ',f10.6,', expected ',f10.6)",&
644              i, r(i),  exp((i-1)*amesh)*r(1)
645         CALL upf_error ('cpmd2upf', 'unknown real-space grid',1)
646      ENDIF
647   ENDDO
648   RETURN
649 END SUBROUTINE check_radial_grid
650! ------------------------------------------------------------------
651FUNCTION mygamma ( n )
652  !------------------------------------------------------------------
653  !
654  ! mygamma(n)  = \Gamma(n-1/2) = sqrt(pi)*(2n-3)!!/2**(n-1)
655  !
656  USE upf_const, ONLY : pi
657  IMPLICIT NONE
658  REAL(dp) :: mygamma
659  INTEGER, INTENT(in) :: n
660  INTEGER :: i
661  !
662  IF ( n < 2 ) CALL upf_error('mygamma','unexpected input argument',1)
663  mygamma = sqrt(pi) / 2.d0**(n-1)
664  ! next factor is (2n-3)!!
665  do i = 2*n-3, 1, -2
666     mygamma = i*mygamma
667  end do
668  !
669END FUNCTION mygamma
670
671END MODULE cpmd_module
672