2
subroutine pseudv(pstype,icorr,ispp,nr,a,b,r,rab,nameat,norb,
3
& ncore,no,lo,so,zo,znuc,zel,cdd,cdu,cdc,viod,
4
& viou,vid,viu,vod,vou,etot,ev,ek,ep)
6
c *************************************************************
8
c * This routine was written by Norman J. Troullier Jr. *
9
c * Nov. 1989, while at the U. of Minnesota, all *
10
c * comments concerning this routine should be directed *
13
c * troullie@128.101.224.101 *
14
c * troullie@csfsa.cs.umn.edu *
17
c * pseudv generates a pseudopotential using the *
18
c * scheme of D. Vanderbilt, ref. Physical Review B, *
19
c * vol. 32, num 12, page 8412. *
20
c * The general format of this routine is the same as the *
21
c * pseudo, pseudk and pseudt routines. Output/input is *
24
c *************************************************************
27
double precision zero, deltas, tpfive, one, two
28
parameter (zero=0.D0,deltas=5.D-3,tpfive=2.5D0,one=1.D0,two=2.D0)
29
double precision small, small2, small3, pzfive
30
parameter (small=1.D-32,small2=1.D-8,small3=1.D-16,pzfive=0.05D0)
31
double precision pfive, small4, ai
32
parameter (pfive=0.5D0,small4=1.D-6,ai=2*137.0360411D0)
33
double precision onepf, oneh
34
parameter (onepf=1.5D0,oneh=100.D0)
35
integer lmax, norbmx, nrmax
36
parameter (lmax=4,norbmx=40,nrmax=1000)
38
C .. Scalar Arguments ..
39
double precision a, b, zel, znuc
40
integer ncore, norb, nr, pstype
41
character ispp*1, icorr*2, nameat*2
43
C .. Array Arguments ..
44
double precision cdc(*), cdd(*), cdu(*), ek(*), ep(*), etot(10),
45
& ev(*), r(*), rab(*), so(*), vid(*), viod(lmax,*),
46
& viou(lmax,*), viu(*), vod(*), vou(*), zo(*)
50
double precision ac, ag, arp, arpm, b0, b2, b4, bc, c1, c2, c3,
51
& cdcp, cfac, cl, coshxr, dcl, delta, dev, devold,
52
& ecut, eviae, fcut, gamma, gg, pi, rbnew, rbold,
53
& rcfac, rcopf, rextr, rmind, rminu, rrc, rzero,
54
& sinhb2, sinhxr, tanb, vap, vapp, viodj, viouj,
55
& vp2z, vps, vpsdm, vpsum, vrc, xr, zeld, zelt,
56
& zelu, zion, zot, zratio, zval
57
integer i, icore, ifcore, iflag, ifull, ist, j, j3rc, jcut, jrc,
58
& k, ka, ll, llp, lp, ncp, nextr, noi, npotd, npotu
59
character blank*1, irel*3, nicore*4
62
double precision ab(5), aj(3,3), bj(3), coe(5), rc(5), rcut(10),
64
integer indd(5), indu(5)
65
character il(5)*1, iray(6)*10, ititle(7)*10
67
C .. External Subroutines ..
68
external difnrl, difrel, dsolv2, etotal, ext, gaussj, polcoe,
69
& potran, potrv, potrw, velect, wtrans, zedate
71
C .. Intrinsic Functions ..
72
intrinsic abs, atan, cosh, exp, log, sin, sinh, sqrt
74
double precision ar(nrmax), arps(nrmax), br(nrmax), f(nrmax),
75
& g(nrmax), v(nrmax), wk1(nrmax), wk2(nrmax),
76
& wk3(nrmax), wk4(nrmax), wk5(nrmax), wk6(nrmax),
81
C .. Data statements ..
83
data il/'s', 'p', 'd', 'f', 'g'/
89
if (ncore .eq. norb) return
93
c Spin-polarized potentails should be unscreened with
94
c a spin-polarized valence charge. This was not
95
c done in pseudo and pseudk in earlier versions
98
if (ispp .eq. 's') then
104
c read rc(s),rc(p),rc(d),rc(f),cfac,rcfac
106
c cfac is used for the pseudocore - the pseudocore stops where
107
c the core charge density equals cfac times the renormalized
108
c valence charge density (renormalized to make the atom neutral).
109
c If cfac is input as negative, the full core charge is used,
110
c if cfac is input as zero, it is set equal to one.
111
c rcfac is used for the pseudocore cut off radius. If set
112
c to less then or equal to zero cfac is used. cfac must be
113
c set to greater then zero.
115
read(5,9000) (rc(i),i=1,4), cfac, rcfac
117
if (cfac .eq. zero) cfac = one
119
c Reset vod and vou to zero. They are here used to store
120
c the pseudo valence charge density.
130
9010 format(//a2,' Pseudopotential Vanderbilt generation',/1x,50('-'),
131
& //' nl s eigenvalue',6x,'rc',4x,6x,'cl',9x,'gamma',7x,
134
c start loop over valence orbitals
140
if (so(i) .lt. 0.1D0) then
141
if (indd(lp) .ne. 0) then
148
if (indu(lp) .ne. 0) then
156
& 'error in pseudv - two down spin orbitals of the same ',
157
& /'angular momentum (',i1,') exist')
158
9030 format(//'error in pseudv - two up spin orbitals of the same ',
159
& /'angular momentum (',i1,') exist')
161
c find all electron wave function
166
if (so(i) .lt. 0.1D0) then
168
v(j) = viod(lp,j)/r(j) + vid(j)
172
v(j) = viou(lp,j)/r(j) + viu(j)
175
if (ispp .ne. 'r') then
177
v(j) = v(j) + llp/r(j)**2
180
if (ispp .ne. 'r') then
181
call difnrl(0,i,v,ar,br,nr,a,b,r,rab,norb,no,lo,so,znuc,
182
& viod,viou,vid,viu,ev,iflag)
184
call difrel(0,i,v,ar,br,nr,a,b,r,rab,norb,no,lo,so,znuc,
185
& viod,viou,vid,viu,ev)
188
c njtj *** plotting routines ***
189
c potrw is called to save an usefull number of points
190
c of the wave function to make a plot. The info is
191
c written to the current plot.dat file.
194
if (ar(nr-85) .lt. zero) ist = -1
195
call potrw(ar,r,nr-85,lo(i),1,ist,rc(lp))
197
c njtj *** user should adjust for their needs ***
199
c Find the last zero and extremum.
202
if (so(i) .lt. 0.1D0 .and. lo(i) .ne. 0) ka = -lo(i)
203
nextr = no(i) - lo(i)
207
if (ispp .eq. 'r') then
208
if (so(i) .lt. 0.1D0) then
209
arp = ka*ar(2)/r(2) + (ev(i)-viod(lp,2)/r(2)-vid(2)+
212
arp = ka*ar(2)/r(2) + (ev(i)-viou(lp,2)/r(2)-viu(2)+
218
if (nextr .eq. 0) go to 80
219
if (ar(j-1)*ar(j) .le. zero) rzero = (ar(j)*r(j-1)-
220
& ar(j-1)*r(j))/(ar(j)-ar(j-1))
224
if (ispp .eq. 'r') then
225
if (so(i) .lt. 0.1D0) then
226
arp = ka*ar(j)/r(j) + (ev(i)-viod(lp,j)/r(j)-vid(j)+
229
arp = ka*ar(j)/r(j) + (ev(i)-viou(lp,j)/r(j)-viu(j)+
234
if (arp*arpm .gt. zero) go to 70
235
rextr = (arp*r(j-1)-arpm*r(j))/(arp-arpm)
239
c Check rc, if outside bounds reset.
242
if (rzero .lt. r(2)) rzero = r(2)
243
if (rc(lp) .gt. rzero .and. rc(lp) .lt. rextr) go to 90
244
if (rc(lp) .ge. zero) write(6,9040) rc(lp), rextr
245
9040 format(/'Warning, the Core radius =',f5.2,
246
& /' is larger then wave function',' extrema position =',
248
if (rc(lp) .lt. zero) rc(lp) = rzero - rc(lp)*(rextr-rzero)
250
c Find the index for grid point closest to 1.5*rc.
251
c Find the index for 3*rc which is used for matching norms.
256
if (r(j) .le. rcopf) then
259
if (r(j) .lt. 3*rc(lp)) then
264
c Reset the n quantum numbers.
271
c Set up potential vl1, first find true potential,
272
c its first and second derivative at rc. Store new
273
c potential(unscreen it first, screening added back
276
if (so(i) .lt. 0.1D0) then
277
vrc = viod(lp,jrc)/r(jrc) + vid(jrc)
278
ab(1) = viod(lp,jrc-2)/r(jrc-2) + vid(jrc-2)
279
ab(2) = viod(lp,jrc-1)/r(jrc-1) + vid(jrc-1)
281
ab(4) = viod(lp,jrc+1)/r(jrc+1) + vid(jrc+1)
282
ab(5) = viod(lp,jrc+2)/r(jrc+2) + vid(jrc+2)
284
vrc = viou(lp,jrc)/r(jrc) + viu(jrc)
285
ab(1) = viou(lp,jrc-2)/r(jrc-2) + viu(jrc-2)
286
ab(2) = viou(lp,jrc-1)/r(jrc-1) + viu(jrc-1)
288
ab(4) = viou(lp,jrc+1)/r(jrc+1) + viu(jrc+1)
289
ab(5) = viou(lp,jrc+2)/r(jrc+2) + viu(jrc+2)
291
rr(1) = r(jrc-2) - r(jrc)
292
rr(2) = r(jrc-1) - r(jrc)
294
rr(4) = r(jrc+1) - r(jrc)
295
rr(5) = r(jrc+2) - r(jrc)
296
call polcoe(rr,ab,5,coe)
309
aj(2,3) = 4*r(jrc)**3
310
aj(3,3) = 12*r(jrc)**2
311
call gaussj(aj,3,3,bj,1,1)
315
if (so(i) .lt. 0.1D0) then
317
viod(lp,j) = ((b0+b2*r(j)**2+b4*r(j)**4)-vid(j))*r(j)
321
viou(lp,j) = ((b0+b2*r(j)**2+b4*r(j)**4)-viu(j))*r(j)
325
c Set up the functions f(r/rc) and g(r/rc) and modify the ionic potential.
333
sinhb2 = (sinh(one))**2
336
rrc = r(j)/rc(lp)/onepf
337
f(j) = oneh**(-((sinh(rrc))**2)/sinhb2)
338
if (f(j) .lt. small2) f(j) = zero
341
if (so(i) .lt. 0.1D0) then
343
viod(lp,j) = viod(lp,j) + dcl*f(j)*r(j)
347
viou(lp,j) = viou(lp,j) + dcl*f(j)*r(j)
352
c Start the iteration loop to find cl.
357
call dsolv2(j,2,blank,ifcore,nr,a,b,r,rab,norb,ncore,nops,
358
& lo,so,zo,znuc,cdd,cdu,cdc,viod,viou,vid,viu,ev,
362
c The abs(dev-devold) condition was added to eliminate
363
c division by zero errors in the calculation of
364
c dcl = -dev*dcl / (dev-devold).
366
if ((abs(dev).lt.small2.or.abs(dev-devold).lt.small3) .and.
372
if (j .gt. 15 .or. abs(dev) .lt. 0.001D0) then
374
c Use newton raphson iteration to change cl.
376
dcl = -dev*dcl/(dev-devold)
378
if (dev*dcl .le. zero) then
384
c Find the new potential.
386
if (so(i) .lt. 0.1D0) then
388
viod(lp,k) = viod(lp,k) + dcl*r(k)*f(k)
392
viou(lp,k) = viou(lp,k) + dcl*r(k)*f(k)
399
c End the iteration loop for cl.
403
c Find the new pseudo-wavefunction.
406
if (so(i) .lt. 0.1D0) then
408
v(j) = (viod(lp,j)+llp/r(j))/r(j) + vid(j)
412
v(j) = (viou(lp,j)+llp/r(j))/r(j) + viu(j)
418
call difnrl(0,i,v,arps,br,nr,a,b,r,rab,norb,nops,lo,so,znuc,
419
& viod,viou,vid,viu,ev,iflag)
421
c Compute yl store in g, store ln(arps) in br.
427
br(j) = log(arps(j)+small)
430
c Compute delta and gamma.
432
gamma = abs(ar(j3rc)/arps(j3rc)+ar(j3rc+1)/arps(j3rc+1))/2
437
ag = ag + ll*arps(j)*g(j)*rab(j)
438
gg = gg + ll*g(j)*g(j)*rab(j)
443
delta = sqrt((ag/gg)**2+(one/gamma**2-one)/gg) - ag/gg
445
c Modify the pseudo-wavefunction.
448
arps(j) = gamma*arps(j)*(one+delta*f(j))
451
c Find d(ln(wl)/dr and store in g(). Note the use of additional
452
c given information of the Vanderbilt method, i.e. the use of
453
c d(ln(wl)/dr to improve stability.
461
rr(1) = r(j-2) - r(j)
462
rr(2) = r(j-1) - r(j)
464
rr(4) = r(j+1) - r(j)
465
rr(5) = r(j+2) - r(j)
466
call polcoe(rr,ab,5,coe)
481
call polcoe(rr,ab,5,coe)
493
call polcoe(rr,ab,5,coe)
496
c Find constants for inversion.
498
c3 = log(oneh)/onepf/rc(lp)/sinhb2
499
c2 = 2/onepf/rc(lp)*c3
502
c Modify potential and find total charge density.
504
if (so(i) .lt. 0.1D0) then
506
vod(j) = vod(j) + zo(i)*arps(j)*arps(j)
510
vou(j) = vou(j) + zo(i)*arps(j)*arps(j)
513
if (so(i) .lt. 0.1D0) then
515
xr = two*r(j)/rc(lp)/onepf
518
viod(lp,j) = viod(lp,j) + delta*f(j)/(one+delta*f(j))*
519
& (c1*(sinhxr)**2-c2*coshxr-2*c3*sinhxr*g(j))*
524
xr = two*r(j)/rc(lp)/onepf
527
viou(lp,j) = viou(lp,j) + delta*f(j)/(one+delta*f(j))*
528
& (c1*(sinhxr)**2-c2*coshxr-2*c3*sinhxr*g(j))*
533
c njtj *** plotting routines ***
534
c potrw is called to save a usefull number of points
535
c of the pseudowave function to make a plot. The
536
c info is written to the current plot.dat file.
537
c wtrans is called to fourier transform the the pseudo
538
c wave function and save it to the current plot.dat file.
541
if (arps(nr-85) .lt. zero) ist = -1
542
call potrw(arps,r,nr-85,lo(i),0,ist,rc(lp))
543
call wtrans(arps,r,nr,lo(i),ist,wk1,wk2,wk3)
545
c njtj *** user should adjust for their needs ***
547
write(6,9050) nops(i), il(lp), so(i), ev(i), rc(lp), cl,
549
9050 format(1x,i1,a1,f6.1,2f12.6,f12.3,2f12.4)
552
c End the loop over the valence orbitals.
554
c Reset the n quantum numbers to include all valence orbitals.
555
c Compute the ratio between the valence charge present and the
556
c valence charge of a neutral atom.
557
c Transfer pseudo valence charge to charge array
565
zion = zval + znuc - zel
566
if (zval .ne. zero) zratio = zion/zval
576
c If a core correction is indicated construct pseudo core charge
577
c cdc(r) = ac*r * sin(bc*r) inside r(icore)
578
c if cfac < 0 or the valence charge is zero the full core is used
580
if (ifcore .ne. 0) then
584
if (cfac .le. zero .or. zratio .eq. zero) then
585
write(6,9060) r(icore), ac, bc
587
if (rcfac .le. zero) then
589
if (cdc(i) .gt. cfac*zratio*
590
& (cdd(i)+cdu(i))) go to 390
594
if (r(i) .le. rcfac) go to 390
599
cdcp = (cdc(icore+1)-cdc(icore))/(r(icore+1)-r(icore))
600
tanb = cdc(icore)/(r(icore)*cdcp-cdc(icore))
603
rbnew = pi + atan(tanb*rbold)
604
if (abs(rbnew-rbold) .lt. .00001D0) then
606
ac = cdc(icore)/(r(icore)*sin(rbnew))
608
cdc(j) = ac*r(j)*sin(bc*r(j))
610
write(6,9060) r(icore), ac, bc
622
9060 format(//' core correction used',/' pseudo core inside r =',f6.3,
623
& /' ac =',f6.3,' bc =',f6.3,/)
624
9070 format(//' error in pseudv - noncovergence in finding ',
625
& /'pseudo-core values')
627
c End the pseudo core charge.
628
c Compute the potential due to pseudo valence charge.
631
c Spin-polarized potentails should be unscreend with
632
c spin-polarized valence charge. This was not
633
c done in pseudo and pseudok in earlier versions
638
if (ispp .eq. 's') then
643
call velect(0,1,icorr,blank,ifcore,nr,r,rab,zval,cdd,cdu,cdc,vod,
646
c Construct the ionic pseudopotential and find the cutoff,
647
c ecut should be adjusted to give a reassonable ionic cutoff
648
c radius, but should not alter the pseudopotential, ie.,
649
c the ionic cutoff radius should not be inside the pseudopotential
657
if (so(i) .lt. 0.1D0) then
659
viod(lp,j) = viod(lp,j) + (vid(j)-vod(j))*r(j)
660
vp2z = viod(lp,j) + 2*zion
661
if (abs(vp2z) .gt. ecut) jcut = j
663
rcut(i-ncore) = r(jcut)
665
fcut = exp(-5*(r(j)-r(jcut)))
666
viod(lp,j) = -2*zion + fcut*(viod(lp,j)+2*zion)
669
v(j) = viod(lp,j)/r(j)
673
viou(lp,j) = viou(lp,j) + (viu(j)-vou(j))*r(j)
674
vp2z = viou(lp,j) + 2*zion
675
if (abs(vp2z) .gt. ecut) jcut = j
677
rcut(i-ncore) = r(jcut)
679
fcut = exp(-5*(r(j)-r(jcut)))
680
viou(lp,j) = -2*zion + fcut*(viou(lp,j)+2*zion)
683
v(j) = viou(lp,j)/r(j)
687
c njtj *** plotting routines ***
689
call potran(lo(i)+1,v,r,nr,zion,wk1,wk2,wk3)
690
call potrv(v,r,nr-120,lo(i),zion)
692
c njtj *** user should adjust for their needs ***
696
c njtj *** plotting routines ***
697
c The calls to 1)potran take the fourier transform of
698
c the potential and saves it in the current plot.dat file,
699
c 2)potrv saves the potential in the current plot.dat file
700
c 3)zion is saved to the current plot.dat file wtih a
701
c marker 'zio' for latter plotting
705
c 9090 format(1x,'marker zio')
706
c 9100 format(2x,f5.2)
708
c njtj *** user should adjust for their needs ***
710
c Convert spin-polarized potentials back to nonspin-polarized
711
c by occupation weight(zo). Assumes core polarization is
712
c zero, ie. polarization is only a valence effect.
714
if (ispp .eq. 's') then
715
do 520 i = ncp, norb, 2
717
zot = zo(i) + zo(i+1)
718
if (zot .ne. zero) then
720
viod(lp,j) = (viod(lp,j)*zo(i)+viou(lp,j)*zo(i+1))/zot
721
viou(lp,j) = viod(lp,j)
725
viod(lp,j) = viod(lp,j)/2 + viou(lp,j)/2
726
viou(lp,j) = viod(lp,j)
737
c Test the pseudopotential self consistency. Spin-polarized
738
c is tested as spin-polarized(since up/down potentials are
741
call dsolv2(0,1,blank,ifcore,nr,a,b,r,rab,norb,ncore,nops,lo,so,
742
& zo,znuc,cdd,cdu,cdc,viod,viou,vid,viu,ev,ek,ep)
744
c Printout the pseudo eigenvalues after cutoff.
746
write(6,9110) (il(lo(i)+1),rcut(i-ncore),i=ncp,norb)
747
write(6,9120) (ev(i),i=ncp,norb)
748
9110 format(//' test of eigenvalues',//' rcut =',8(2x,a1,f7.2))
749
9120 format(' eval =',8(2x,f8.5))
751
c Printout the data for potentials.
754
9130 format(///' l vps(0) vpsmin at r',/)
756
if (indd(i)+indu(i) .eq. 0) go to 560
757
if (indd(i) .ne. 0) then
760
if (r(j) .lt. .00001D0) go to 540
762
if (vps .lt. vpsdm) then
767
write(6,9140) il(i), viod(i,2)/r(2), vpsdm, rmind
769
if (indu(i) .ne. 0) then
772
if (r(j) .lt. .00001D0) go to 550
774
if (vps .lt. vpsum) then
779
write(6,9140) il(i), viou(i,2)/r(2), vpsum, rminu
781
9140 format(1x,a1,3f10.3)
784
c Print out the energies from etotal.
786
call etotal(pstype,one,nameat,norb-ncore,nops(ncp),lo(ncp),
787
& so(ncp),zo(ncp),etot,ev(ncp),ek(ncp),ep(ncp))
789
c Find the jobname and date, date is a machine
790
c dependent routine and must be chosen/written/
791
c comment in/out in the zedate section.
793
iray(1) = 'atom-lda '
795
iray(3) = 'Vanderbilt'
796
iray(4) = ' Pseudo - '
797
iray(5) = 'potential '
798
iray(6) = 'generation'
800
c Encode the title array.
806
if (indd(i) .eq. 0 .and. indu(i) .eq. 0) go to 580
809
if (indd(i) .ne. 0) then
813
if (indu(i) .ne. 0) then
818
if (ispp .ne. 's') then
819
write(ititle(2*i-1),9150) noi, il(i), zelt
820
write(ititle(2*i),9160) ispp, rc(i)
821
9150 format(' ',i1,a1,'(',f5.2,')')
822
9160 format(a1,' rc=',f5.2)
824
write(ititle(2*i-1),9170) noi, il(i), zeld
825
write(ititle(2*i),9180) zelu, ispp, rc(i)
826
9170 format(i1,a1,' (',f4.2,',')
827
9180 format(f4.2,')',a1,f4.2)
831
c Construct relativistic sum and difference potentials.
833
if (ispp .eq. 'r') then
834
if (indu(1) .eq. 0) go to 600
838
viod(1,j) = viou(1,j)
843
if (indd(i) .eq. 0 .or. indu(i) .eq. 0) go to 620
847
viod(i,j) = ((i-1)*viodj+i*viouj)/(2*i-1)
848
viou(i,j) = 2*(viouj-viodj)/(2*i-1)
853
c Determine the number of potentials. Coded them as
854
c two digits, where the first digit is the number
855
c of down or sum potentials and the second the number of
856
c up or difference potentials.
861
if (indd(i) .ne. 0) npotd = npotd + 1
862
if (indu(i) .ne. 0) npotu = npotu + 1
865
c Write the heading to the current pseudo.dat
869
if (cfac .le. zero .or. zratio .eq. zero) ifull = 1
870
if (ifcore .eq. 1) then
871
if (ifull .eq. 0) then
876
else if (ifcore .eq. 2) then
877
if (ifull .eq. 0) then
885
if (ispp .eq. 's') then
887
else if (ispp .eq. 'r') then
893
write(1) nameat, icorr, irel, nicore, (iray(i),i=1,6),
894
& (ititle(i),i=1,7), npotd, npotu, nr - 1, a, b, zion
895
write(1) (r(i),i=2,nr)
897
c Write the potentials to the current pseudo.dat
901
if (indd(i) .eq. 0) go to 640
902
write(1) i - 1, (viod(i,j),j=2,nr)
905
if (indu(i) .eq. 0) go to 650
906
write(1) i - 1, (viou(i,j),j=2,nr)
909
c Write the charge densities to the current pseudo.dat
912
if (ifcore .ne. 1) then
913
write(1) (zero,i=2,nr)
915
write(1) (cdc(i),i=2,nr)
917
write(1) (zratio*(cdd(i)+cdu(i)),i=2,nr)
added
removed
Pseudo/atom/Contrib/pseudv.for
