c c Copyright (C) 2002 Vanderbilt group c This file is distributed under the terms of the GNU General Public c License as described in the file 'License' in the current directory. c c---------------------------------------------------------------------------- c c ****************************************************** subroutine psrv(rho,rcut,ifv,r,mesh,ifprt,idim1) c ****************************************************** c c this subroutine pseudises the function rho to a taylor c series in r**2 between origin and rcut. the function c and its derivative are made continuous at rcut. c c designed for use either on partial core density (ifv=0) c or local potential (ifv=1) c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c c c.....logarithmic radial mesh information dimension r(idim1) c.....function to be interpolated dimension rho(idim1) c.....scratch dimension targ(2),con(2),bb(2,2) c index on targ and 2nd index of bb is (value,deriv) c index on con, and 1st index of bb taylor coeff of ( 1 , r^2 ) c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c c---------------------------------------------------------------------------- c c i n i t i a l i s a t i o n c if (ifprt.gt.-2) write (iout,'(/a)') ' entering subroutine psrv' c c ifv = 1 if this is really a potential, not a charge density c c rho is really 4*pi*r**2*rho c v is really r*v c body of subroutine assumes proportional to v or rho do 10 i=2,mesh rho(i) = rho(i) / r(i)**(2-ifv) 10 continue c c---------------------------------------------------------------------------- c c f i n d i n t e r p o l a t i o n p o i n t s c do 100 i2 = mesh,1,-1 if ( r(i2) .lt. rcut ) goto 110 100 continue c 110 continue c c check that the value of i2 is reasonable if ( i2 .lt. 3 .or. i2 .gt. mesh-2 ) then write(iout,*) '***error in subroutine psrv' write(iout,*) 'rcut =',rcut,' generated i2 =',i2, + ' which is illegal with mesh =',mesh call exit(0) endif c i1 = i2 - 1 i3 = i2 + 1 i4 = i2 + 2 c if ( ifprt .ge. 5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'rh1,rh2,rh3,rh4 =',(rho(i),i=i1,i4) 120 format(a18,4f8.4) write (iout,130) rcut 130 format(' rcut =',f8.4) endif c c---------------------------------------------------------------------------- c c s e t u p t h e p o l y n o m i a l i n t e r p o l a t i o n c call polyn(rcut,val,r(i1),rho(i1),-1) call polyn(rcut,val,r(i1),rho(i1), 0) call polyn(rcut, d1,r(i1),rho(i1),+1) c c---------------------------------------------------------------------------- c c s o l v e s i m u l t a n e o u s e q u a t i o n s c targ(1) = val targ(2) = d1 c bb(1,1) = 1.0d0 bb(2,1) = 0.0d0 bb(1,2) = rcut**2 bb(2,2) = 2.0d0 * rcut c c if (ifprt .ge. 5) then write (iout,150) val,d1 150 format (' val,d1=',2f15.5) write (iout,155) 'val ',(bb(1,j),j=1,2),targ(1) write (iout,155) 'd1 ',(bb(2,j),j=1,2),targ(2) 155 format(1x,a4,3e15.5) endif c call simeq(2,2,bb,targ,con) c if ( ifprt .ge. 5) write (iout,160) (con(i),i=1,2) 160 format (' constants from simeq are ',2e15.5) c c---------------------------------------------------------------------------- c c f i l l r h o w i t h p s e u d i s e d f o r m c do 200 i=2,i2 rr=r(i)**2 rho(i)=con(1)+rr*con(2) 200 continue c if ( ifprt .gt. 5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'rh1,rh2,rh3,rh4 =',(rho(i),i=i1,i4) endif c do 210 i=2,mesh rho(i) = rho(i) * r(i)**(2-ifv) 210 continue c if (ifprt.gt.-2) write (iout,'(a)') ' exiting subroutine psrv' c return end c c----------------------------------------------------------------------- c subroutine psv(ruae,vloc,rcloc,rcmax,r,rab,sqr,a,b,mesh, + irel,yy,bndmat, + kkbeta,klogd,dum1,dum2,dum3,lloc,eloc,ifprt,idim1,iploctype, + ikeyeeloc,refwf, + nang,npf,ptryc,idim8,idim3) c c pseudization of the local potential to obtain correct c logarithmic derivative in channel lloc at energy eloc c implicit double precision(a-h,o-z) c c c.....logarithmic radial mesh information dimension r(idim1),rab(idim1),sqr(idim1) c.....potentials for the all electron calculation dimension ruae(idim1) c.....local potential dimension vloc(idim1) c.....for the phakoe call dimension bndmat(idim1,5),yy(idim1,2) c.....type of wave function dimension refwf(idim1,idim3) c.....scratch arrays dimension dum1(idim1),dum2(idim1),dum3(idim1),con(10) c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c c----------------------------------------------------------------------- c c i n i t i a l i z a t i o n c if ( ifprt .ge. 1 ) write (iout,'(/3x,a)') + 'subroutine psv: pseudize the local potential' write (iout,'(/3x,a,i2,a,f9.5)') 'generating vloc with lloc ', + lloc,' eloc =',eloc c if ( lloc .lt. 0 .or. lloc .gt. 3 ) then write(iout,*) '***error subroutine psv' write(iout,*) 'lloc is not reasonable, lloc =',lloc call exit(1) endif c c----------------------------------------------------------------------- c c g e n e r a t e t h e l o c a l p o t e n t i a l c c obtain all electron wavefunction at energy eloc in dum2 if (ikeyeeloc .lt. 0) then do i=1,mesh dum2(i) = refwf(i,3-ikeyeeloc) enddo do i=1,mesh if (r(i) .gt. rcloc) goto 203 enddo 203 if (dum2(i) .lt. 0.0) then do i=1,mesh dum2(i) = -1.0*dum2(i) enddo endif elseif ( irel .eq. 0 ) then c write(iout,'(3x,a)') 'calling phase' zz = 0.0d0 vzero = ruae(2) / r(2) call phase(zz,vzero,eloc,lloc,b,r,rab,sqr, + ruae,dum1,dum2,mesh,kkbeta,dlwf,idim1) c c renormalise psi and remove square root type factor factor = 1.0d0 / ( sqr(klogd) * dum2(klogd) ) do i = 1,kkbeta dum2(i) = dum2(i) * sqr(i) * factor enddo c elseif ( irel .eq. 2 ) then c write(iout,'(3x,a)') 'calling phakoe' call phakoe(ruae,r,rab,bndmat,eloc,lloc, + yy,dum1,kkbeta,dlwf,idim1) c c copy major component of koelling equation into psi c and renormalise factor = 1.0d0 / yy(klogd,2) do 205 ir = 1,kkbeta dum2(ir) = yy(ir,2) * factor 205 continue c else c write(iout,*) '***error in psw' write(iout,*) 'irel =',irel,' is not allowed' call exit(1) c endif c c nodeless pseudization of the wavefunction if(iploctype .eq. 1) then write(iout,'(3x,a)') 'calling pswf1' call pswf1(dum2,dum1,rcloc,lloc,r,rab,con,kkbeta, + ifprt,idim1) c c load vloc using analytic inversion of schroedinger equation vloc(1) = 0.0d0 do 120 i = 2,mesh if ( r(i) .lt. rcloc ) then t1 = 2.0d0*con(2)*r(i) + 4.0d0*con(3)*r(i)**3 + + 6.0d0*con(4)*r(i)**5 t2 = 2.0d0*con(2) + 12.0d0*con(3)*r(i)**2 + + 30.0d0*con(4)*r(i)**4 vtemp = eloc + dble(2*lloc+2)*t1/r(i) + t1**2 + t2 c remember that vloc is potential * r vloc(i) = r(i) * vtemp else c no pseudization necessary beyond rcloc vloc(i) = ruae(i) endif 120 continue c elseif(iploctype .eq. 2) then write(iout,'(3x,a)') 'Not implemented' call exit(1) c c this is a little mysterious. it seems kurt stokbro put in c the capability for iploctype=2 in version a6.0 (it was not c there in a5.1.3, where pswf1 was always used to generate the c pseudo wf as the intermediate step to pseudizing the local c potential). but for some reason kurt put this 'Not c implemented' abort into the code. maybe he just never tested c it, or there was a bug he never bothered to track down. c - dv 8/17/97 c c! write(iout,*) ' calling pswf2' c! call pswf2(dum2,dum1,rcloc,lloc,r,rab,con,kkbeta, c! + ifprt,nang,npf,ptryc,idim1,10) c! c! load vloc using analytic inversion of schroedinger equation c! vloc(1) = 0.0d0 c! do i = 2,mesh c! if ( r(i) .lt. rcloc ) then c! t1 = 0.d0 c! t2 = 0.d0 c! do j =2,npf c! t1 = t1 + dble(2*(j-1))*con(j)*r(i)**(2*j-3) c! t2 = t2 + dble((2*j-2)*(2*j-3))*con(j)*r(i)**(2*j-4) c! enddo c! vtemp = eloc + dble(2*lloc+2)*t1/r(i) + t1**2 + t2 c! remember that vloc is potential * r c! vloc(i) = r(i) * vtemp c! else c! no pseudization necessary beyond rcloc c! vloc(i) = ruae(i) c! endif c! enddo c elseif (iploctype .eq. 3) then c do i = 1,mesh if (ikeyeeloc .lt. 0) then vloc(i) = refwf(i,1) else vloc(i) = ruae(i) endif enddo write(iout,'(3x,a)') 'calling pswfnormc' call pswfnormc(dum2,dum1,vloc,rcloc,lloc,eloc,r,rab,b, + con,kkbeta,ifprt,idim1,10) nqtr = 7 c c load vloc using analytic inversion of schroedinger equation write(iout,'(3x,a,f8.3)') 'rcloc is now',rcloc vloc(1) = 0.0d0 do i = 2,mesh if ( r(i) .le. rcloc ) then t1 = 0.d0 t2 = 0.d0 do j =2,nqtr t1 = t1 + dble(2*(j-1))*con(j)*r(i)**(2*j-3) t2 = t2 + dble((2*j-2)*(2*j-3))*con(j)*r(i)**(2*j-4) enddo vtemp = eloc + dble(2*lloc+2)*t1/r(i) + t1**2 + t2 c y = r(i)*r(i) c pr = (((((con(7)*y+con(6))*y+con(5))*y+con(4))*y+ c + con(3))*y+con(2))*y + con(1) c t1=exp(pr)*r(i)**(lloc+1) c write(*,*) r(i),vloc(i),r(i) * vtemp,dum2(i),t1 c remember that vloc is potential * r vloc(i) = r(i) * vtemp if (ikeyeeloc .lt. 0) then vloc(i) = vloc(i)-refwf(i,2)+refwf(i,3) endif else c no pseudization necessary beyond rcloc vloc(i) = ruae(i) c write(*,*) r(i),vloc(i),vloc(i),dum2(i),dum2(i) endif enddo c else write(iout,*) '***error in subroutine psv' write(iout,*) 'iploctype=',iploctype,'not implemented' call exit(1) c endif c c check the results by generating wavefunction with vloc write(iout,'(3x,a)') 'calling phase' zz = 0.0d0 vzero = vloc(2) / r(2) call phase(zz,vzero,eloc,lloc,b,r,rab,sqr, + vloc,dum1,dum3,mesh,kkbeta,dlwf,idim1) c c renormalise psi to remove square root type factor factor = 1.0d0 / ( sqr(klogd) * dum3(klogd) ) do 130 i = 1,kkbeta dum3(i) = dum3(i) * sqr(i) * factor 130 continue c c possible print out of results if ( ifprt .ge. 1) then write(iout,9980) 9980 format(/3x,'results for the local potential vloc',/, + 8x,'r',9x,'vloc',9x,'dum2',9x,'dum3') nind = 55 rinc = rcmax / dble(nind) rcur = -0.75d0 * rinc do 190 i = 1,nind rcur = rcur + rinc in = int( dlog( rcur / a + 1.0d0 ) / b ) + 1 write(iout,9979) r(in),vloc(in),dum2(in),dum3(in) 190 continue endif 9979 format(3x,f9.4,3f13.8) c if ( ifprt .ge. 1 ) write (iout,'(3x,a)') 'leaving psv' return end c c---------------------------------------------------------------------------- c c ************************************************************* subroutine pswf(psi,rho,rcut,lcur,r,rab,con,mesh,ifprt,idim1) c ************************************************************* c this subroutine pseudises the radial wave function to c as a taylor series c c snl(r) = r**(l+1) * ( a + b*r**2 + c*r**4 + d*r**6 ) c c matches value and 1st through third derivatives at rcut c c snl(r) = snl(r) + r**(l+1) * tor * ( r**2 - rcut**2 )**4 c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c c.....logarithmic radial mesh information dimension r(idim1),rab(idim1) c.....wavefunction to be pseudised (psi is really r*psi on input) dimension psi(idim1),rho(idim1) c.....scratch dimension targ(4),con(4),bb(4,4) c index on targ and 2nd index of bb is (val,d1,d2,d3) c index on con, and 1st index of bb taylor coeff of 1, r^2, ... c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c c---------------------------------------------------------------------------- c c i n i t i a l i s a t i o n c if ( ifprt .ge. 0) write (iout,'(/3x,a)') 'subroutine pswf' c c compute the weight under psi rho(1) = 0.0d0 do 10 i = 2,mesh rho(i) = psi(i)**2 10 continue call radin(mesh,rho,rab,asum,idim1) c if (ifprt .ge. 5) write (iout,20) asum 20 format (' integral under unpseudised psi is ',f15.10) c lp = lcur + 1 do 30 i=2,mesh psi(i) = psi(i) / r(i)**lp 30 continue c body of subroutine assumes proportional to psi/r**l c c---------------------------------------------------------------------------- c c f i n d i n t e r p o l a t i o n p o i n t s c do 100 i2 = mesh,1,-1 if ( r(i2) .lt. rcut ) goto 110 100 continue c 110 continue c c check that the value of i2 is reasonable if ( i2 .lt. 3 .or. i2 .gt. mesh-2 ) then write(iout,*) '***error in subroutine pswf' write(iout,*) 'rcut =',rcut,' generated i2 =',i2, + ' which is illegal with mesh =',mesh call exit(0) endif c i1 = i2 - 1 i3 = i2 + 1 i4 = i2 + 2 c if ( ifprt .ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'ps1,ps2,ps3,ps4 =',(psi(i),i=i1,i4) 120 format(a18,4f8.4) write (iout,130) rcut 130 format(' rcut =',f8.4) endif c c---------------------------------------------------------------------------- c c s e t u p t h e p o l y n o m i a l i n t e r p o l a t i o n c call polyn(rcut,val,r(i1),psi(i1),-1) call polyn(rcut,val,r(i1),psi(i1), 0) call polyn(rcut, d1,r(i1),psi(i1),+1) call polyn(rcut, d2,r(i1),psi(i1),+2) call polyn(rcut, d3,r(i1),psi(i1),+3) c if ( ifprt .ge.5) write (iout,140) val,d1,d2,d3 140 format (' val,d1,d2,d3=',4f15.5) c c---------------------------------------------------------------------------- c c s o l v e s i m u l t a n e o u s e q u a t i o n s c targ(1) = val targ(2) = d1 targ(3) = d2 targ(4) = d3 c bb(1,1) = 1.0d0 bb(2,1) = 0.0d0 bb(3,1) = 0.0d0 bb(4,1) = 0.0d0 bb(1,2) = rcut**2 bb(2,2) = 2.0d0 * rcut bb(3,2) = 2.0d0 bb(4,2) = 0.0d0 bb(1,3) = rcut**4 bb(2,3) = 4.0d0 * rcut**3 bb(3,3) = 12.0d0 * rcut**2 bb(4,3) = 24.0d0 * rcut bb(1,4) = rcut**6 bb(2,4) = 6.0d0 * rcut**5 bb(3,4) = 30.0d0 * rcut**4 bb(4,4) = 120.0d0 * rcut**3 c if ( ifprt .ge.5) then write (iout,150) 'val ',(bb(1,j),j=1,4),targ(1) write (iout,150) 'd1 ',(bb(2,j),j=1,4),targ(2) write (iout,150) 'd2 ',(bb(3,j),j=1,4),targ(3) write (iout,150) 'd3 ',(bb(4,j),j=1,4),targ(4) endif 150 format (1x,a4,5e15.5) c call simeq(4,4,bb,targ,con) c if ( ifprt .ge.5) write (iout,160) 'con ',con 160 format (1x,a4,4e15.5) c c---------------------------------------------------------------------------- c c f i l l p s i w i t h p s e u d i s e d f o r m c do 200 i=2,i2 rr = r(i)**2 psi(i) = con(1) + con(2)*rr + con(3)*rr**2 + con(4)*rr**3 200 continue c if ( ifprt .ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'ps1,ps2,ps3,ps4 =',(psi(i),i=i1,i4) endif c c run a self-consistency check call polyn(rcut,val,r(i1),psi(i1),-1) call polyn(rcut,val,r(i1),psi(i1), 0) call polyn(rcut, d1,r(i1),psi(i1),+1) call polyn(rcut, d2,r(i1),psi(i1),+2) call polyn(rcut, d3,r(i1),psi(i1),+3) c if ( ifprt .ge.5) write (iout,140) val,d1,d2,d3 c c----------------------------------------------------------------------- c c r e f o r m r a d i a l f u n c t i o n a n d w e i g h t c psi(1) = 0.0d0 do 210 i=2,mesh psi(i) = psi(i) * r(i)**lp 210 continue c rho(1) = 0.0d0 do 220 i = 2,mesh rho(i) = psi(i)**2 220 continue call radin(mesh,rho,rab,asum,idim1) c if ( ifprt .ge.5) write (iout,230) asum 230 format (' integral under pseudised psi is ',f15.10) c if ( ifprt .ge. 0) write (iout,'(3x,a)') 'leaving pswf' c return end c c---------------------------------------------------------------------------- c subroutine psqf(qfunc,qfcoef,rho,lmin, + rinner,r,rab,a,b,kkbeta,ifprt,nang,idim1,idim6,idim8) c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c c index on targ and 2nd index of bb is (value,deriv,norm) c 1st indices of qcoef and bb are taylor coeff of 1, r^2, r^4 c c ------------------------------------------------------ c notes on qfunc and qfcoef: c ------------------------------------------------------ c since q_ij(r) is the product of two orbitals like c psi_{l1,m1}^star * psi_{l2,m2}, it can be decomposed by c total angular momentum l, where l runs over | l1-l2 | , c | l1-l2 | +2 , ... , l1+l2. (l=0 is the only component c needed by the atomic program, which assumes spherical c charge symmetry.) c c recall qfunc(r) = y1(r) * y2(r) where y1 and y2 are the c radial parts of the wave functions defined according to c c psi(r-vec) = (1/r) * y(r) * y_lm(r-hat) . c c for each total angular momentum l, we pseudize qfunc(r) c inside rc as: c c qfunc(r) = r**(l+2) * [ a_1 + a_2*r**2 + a_3*r**4 ] c c in such a way as to match qfunc and its 1'st derivative at c rc, and to preserve c c integral dr r**l * qfunc(r) , c c i.e., to preserve the l'th moment of the charge. the array c qfunc has been set inside rc to correspond to this pseudized c version using the minimal l, namely l = | l1-l2 | (e.g., l=0 c for diagonal elements). the coefficients a_i (i=1,2,3) c are stored in the array qfcoef(i,l+1,j,k) for each l so that c the correctly pseudized versions of qfunc can be reconstructed c for each l. (note that for given l1 and l2, only the values c l = | l1-l2 | , | l1-l2 | +2 , ... , l1+l2 are ever used.) c ------------------------------------------------------ c c.....logarithmic radial mesh information dimension r(idim1),rab(idim1) c.....q functions and q pseudization coefficients dimension qfunc(idim1),qfcoef(idim8,idim6) c.....scratch dimension rho(idim1) dimension targ(3),bb(3,3) c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c c if ( ifprt .ge. 0) write (iout,'(/3x,a)') 'subroutine psqf' c c---------------------------------------------------------------------------- c c f i n d i n t e r p o l a t i o n p o i n t s c do 100 i2 = kkbeta,1,-1 if ( r(i2) .lt. rinner ) goto 110 100 continue c 110 continue c c check that the value of i2 is reasonable if ( i2 .lt. 3 .or. i2 .gt. kkbeta-2 ) then write(iout,*) '***error in subroutine psqf' write(iout,*) 'rinner =',rinner,' generated i2 =',i2, + ' which is illegal with kkbeta =',kkbeta call exit(0) endif c i1 = i2 - 1 i3 = i2 + 1 i4 = i2 + 2 c if ( ifprt .ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'qf1,qf2,qf3,qf4 =',(qfunc(i),i=i1,i4) 120 format(a18,4f8.4) write (iout,130) rinner 130 format(' rinner =',f8.4) endif c c---------------------------------------------------------------------------- c c b e g i n l o o p o v e r l t o t c do 200 ltot = 0,2*(nang-1) c ltp = ltot + 1 ll2 = 2 * ltot + 2 ll3 = 2 * ltot + 3 ll5 = 2 * ltot + 5 ll7 = 2 * ltot + 7 c do 210 i = 2,kkbeta rho(i) = qfunc(i) / r(i)**(ltot+2) 210 continue c call polyn(rinner,val,r(i1),rho(i1),-1) call polyn(rinner,val,r(i1),rho(i1), 0) call polyn(rinner, d1,r(i1),rho(i1),+1) c targ(1) = val targ(2) = d1 c bb(1,1) = 1.0d0 bb(2,1) = 0.0d0 bb(1,2) = rinner**2 bb(2,2) = 2.0d0 * rinner bb(1,3) = rinner**4 bb(2,3) = 4.0d0*rinner**3 c c integral of r**(2*l+2)*rho should be preserved c c bb(3,j) = integral of r**(2*l+2)*rho bb(3,1) = rinner**ll3/ll3 bb(3,2) = rinner**ll5/ll5 bb(3,3) = rinner**ll7/ll7 c c simpsons rule integration up to rtrue c rho(1)=0.d0 do 220 i = 2,kkbeta rho(i) = qfunc(i) * r(i) ** ltot 220 continue c call simp(rho,r,rab,a,b,rinner,kkbeta,asum,idim1) c targ(3)=asum c if ( ifprt .ge.5) then write (iout,230) val,d1,asum 230 format ('val,d1,asum=',3f15.5) write (iout,240) 'val ',(bb(1,j),j=1,3),targ(1) write (iout,240) 'd1 ',(bb(2,j),j=1,3),targ(2) write (iout,240) 'asum',(bb(3,j),j=1,3),targ(3) 240 format (1x,a4,4e15.5) endif c call simeq(3,3,bb,targ,qfcoef(1,ltp)) c if ( ifprt .ge.5) write (iout,250) 'con ',(qfcoef(i,ltp),i=1,3) 250 format (1x,a4,3e15.5) c c end loop over tot ang momentum components 200 continue c c----------------------------------------------------------------------------- c c s e t q f u n c f o r l t o t = l m i n c lp = lmin+1 do 300 i = 1,i2 rr = r(i)**2 rho(i) = qfcoef(1,lp) + rr*(qfcoef(2,lp) + rr*qfcoef(3,lp)) qfunc(i) = rho(i)*r(i)**(lmin+2) 300 continue c if ( ifprt .ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'qf1,qf2,qf3,qf4 =',(qfunc(i),i=i1,i4) endif c if (ifprt.gt.0) write (iout,*) ' exiting subroutine psqf' c return end c c----------------------------------------------------------------------------- c subroutine simp(rho,r,rab,a,b,rinner,kkbeta,asum,idim1) c c routine for simpson's rule integration up to rinner c implicit double precision (a-h,o-z) c dimension rho(idim1),r(idim1),rab(idim1) c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c x = dlog( rinner / a + 1.0d0 ) / b + 1.0d0 c c check that rinner will not require rho beyond kkbeta if ( int(x) + 2 .gt. kkbeta ) then write(iout,*) '***error in subroutine simp' write(iout,*) 'int(x)+2=',int(x)+2,' .gt. kkbeta =',kkbeta call exit(1) endif c f3 = rho(1) * rab(1) asum = 0.0d0 c c commence the integration loop do 100 i = 1,kkbeta,2 f1 = f3 f2 = rho(i+1) * rab(i+1) f3 = rho(i+2) * rab(i+2) if ( i + 2 .gt. int(x) ) goto 110 asum = asum + ( f1 + 4.0d0 * f2 + f3 ) / 3.0d0 100 continue c 110 continue c x = x - dble(i) aa3 = ( f1 - 2.0d0 * f2 + f3 ) / 6.0d0 aa2 = ( -3.0d0 * f1 + 4.0d0 * f2 - f3 ) / 4.0d0 aa1 = f1 asum = asum + x * ( aa1 + x * ( aa2 + x * aa3 ) ) c return end c c----------------------------------------------------------------------- c subroutine psqf1(qfunc,qfcoef,rho,lmin, + rinner,mesh,r,rab,a,b,kkbeta,ifprt,nang, + nqf,qtryc,idim1,idim6,idim8) c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c c modified by X.-P. Li, in Oct. 1991 c c targ are dimensioned nqf+3, the first nqf of them c are for taylor coeff of 1, r**2, r**4 etc, the last are c (value,deriv,norm) corresponding to 3 lagrange multipliers c c ------------------------------------------------------ c notes on qfunc and qfcoef: c ------------------------------------------------------ c since q_ij(r) is the product of two orbitals like c psi_{l1,m1}^star * psi_{l2,m2}, it can be decomposed by c total angular momentum l, where l runs over | l1-l2 | , c | l1-l2 | +2 , ... , l1+l2. (l=0 is the only component c needed by the atomic program, which assumes spherical c charge symmetry.) c c recall qfunc(r) = y1(r) * y2(r) where y1 and y2 are the c radial parts of the wave functions defined according to c c psi(r-vec) = (1/r) * y(r) * y_lm(r-hat) . c c for each total angular momentum l, we pseudize qfunc(r) c inside rc as: c c qfunc(r) = r**(l+2) * [ a_1 + a_2*r**2 + a_3*r**4 + ...] c c in such a way as to match qfunc and its 1'st derivative at c rc, and to preserve c c integral dr r**l * qfunc(r) , c c i.e., to preserve the l'th moment of the charge. the rest of c the coeffs are determined by minimizing c c integral (from qtryc to infinite) dq q**2 * qfunc(q)**2 c where c qfunc(q)=intgral dr qfunc(r) * j_l(qr) c c where qtryc and the number of coeffs nqf are input. c c qfunc has been set inside rc to correspond to this pseudized c version using the minimal l, namely l = | l1-l2 | (e.g., l=0 c for diagonal elements). the coefficients a_i (i=1,2,3) c are stored in the array qfcoef(i,l+1,j,k) for each l so that c the correctly pseudized versions of qfunc can be reconstructed c for each l. (note that for given l1 and l2, only the values c l = | l1-l2 | , | l1-l2 | +2 , ... , l1+l2 are ever used.) c c ------------------------------------------------------ c c.....in order to avoid passing many unnecessary arrays, use following parameter(idm8=13,idm1=101,idm9=3) c.....logarithmic radial mesh information dimension r(idim1),rab(idim1) c.....q functions and q pseudization coefficients dimension qfunc(idim1),qfcoef(idim8,idim6) c.....scratch dimension rho(idim1) dimension targ(idm8),bb(idm8,idm8,idm9),targinv(idm8) dimension qn(idm1,idm8,idm9),qlf(idm1),g(idm1) c integer stderr logical ifinit c common /files/ stderr,input,iout,ioae,iplot,iologd,iops c c to speed up a little c save qn,bb,g data ifinit/.false./ c if ( ifprt .ge. 0) write (iout,'(/3x,a,l1)') + 'subroutine psqf1, ifinit = ',ifinit c if(idm8.lt.nqf+3) then write(iout,*) '***error: in subroutine psqf1, idm8 too small' write(iout,*) 'idm8, nqf+3 = ',idm8,nqf+3 call exit(1) endif c if(idm9.lt.nang) then write(iout,*) '***error: in subroutine psqf1, idm9 too small' write(iout,*) 'idm9, nang = ',idm9,nang call exit(1) endif c c---------------------------------------------------------------------------- c c f i n d i n t e r p o l a t i o n p o i n t s c do 100 i2 = kkbeta,1,-1 if ( r(i2) .lt. rinner ) goto 110 100 continue c 110 continue c c check that the value of i2 is reasonable if ( i2 .lt. 3 .or. i2 .gt. kkbeta-2 ) then write(iout,*) '***error in subroutine psqf1' write(iout,*) 'rinner =',rinner,' generated i2 =',i2, + ' which is illegal with kkbeta =',kkbeta call exit(0) endif c i1 = i2 - 1 i3 = i2 + 1 i4 = i2 + 2 c if ( ifprt .ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'qf1,qf2,qf3,qf4 =',(qfunc(i),i=i1,i4) 120 format(a18,4f8.4) write (iout,130) rinner 130 format(' rinner =',f8.4) endif c pi=3.14159265358979 tpi=2.0*pi tpi3=tpi**3 fpi=4.0*pi c c minimization if nqf gt 3 c if(nqf.gt.3) then c c---------------------------------------------------------------------------- c c prepare things beforehand c if(.not.ifinit) then c c prepare g's c dg=qtryc/dfloat(idm1-1) do 2010 iq=1,idm1 2010 g(iq)=dble(iq-1)*dg c c prepare qn's c facl=1.0 do 2050 il=1,nang ll=il-1 fll=dble(ll) facl=facl*(2.0*fll+1.0) do 2050 in=1,nqf lp2n=ll+2*in+1 tlnp=dble(2*(ll+in)+1) if(ll.eq.0) then qn(1,in,il)=fpi*rinner**lp2n/facl/tlnp else qn(1,in,il)=0.0d0 endif c do 2050 iq=2,idm1 qdr=g(iq)*rinner hqdr2=0.5*qdr*qdr c termth=1.0d0 ith=0 2020 continue ith=ith+1 di=dble(ith) termth=termth*hqdr2/di/(2.0*(fll+di)+1.0) if(abs(termth).gt.1.0d-20) goto 2020 nterm=ith c termth=0.0d0 do 2030 ith=nterm,1,-1 di=dble(ith) termth=hqdr2/di/(2.0*(fll+di)+1.0) * *(1.0/(tlnp+2.0*di)-termth) 2030 continue qn(iq,in,il)=fpi*rinner**lp2n*qdr**ll/facl * *(1.0/tlnp-termth) 2050 continue c if(ifprt.gt.0) then write(iout,'(5x,a,i5)') 'notice!! idm1 = ',idm1 write(iout,'(5x,a,i5)') 'qns for l = 2, n = ',nqf write(iout,250) (qn(iq,nqf,3),iq=1,idm1) endif c c now construct matrices bb c do 2100 ltp=1,nang c do 2060 iqf=1,nqf+3 do 2060 jqf=1,nqf+3 2060 bb(iqf,jqf,ltp)=0.0 c do 2070 iqf=1,nqf iqft=2*(iqf-1) bb(nqf+1,iqf,ltp)=rinner**iqft bb(nqf+2,iqf,ltp)=dfloat(iqft)*rinner**(iqft-1) c integral of r**(2*l+2)*rho should be preserved lmtp=2*(ltp-1+iqf)+1 bb(nqf+3,iqf,ltp)=rinner**lmtp/dfloat(lmtp) c symmetrical part bb(iqf,nqf+1,ltp)=bb(nqf+1,iqf,ltp) bb(iqf,nqf+2,ltp)=bb(nqf+2,iqf,ltp) bb(iqf,nqf+3,ltp)=bb(nqf+3,iqf,ltp) 2070 continue c do 2090 iqf=1,nqf do 2090 jqf=1,nqf do 2080 iq=1,idm1 2080 rho(iq)=g(iq)*g(iq)*qn(iq,iqf,ltp)*qn(iq,jqf,ltp) call simpagn(idm1,g,rho,asum,0) lijtm=2*(ltp-1+iqf+jqf)-1 bb(iqf,jqf,ltp)=tpi3*rinner**lijtm/dfloat(lijtm)-asum 2090 continue c 2100 continue c endif c c c b e g i n l o o p o v e r l t o t c do 200 ltot = 0,2*(nang-1) c ltp = ltot + 1 c c since only ltot=lmin is used, set the rest to 0 c if(ltot.ne.lmin) then c do 2110 iqf=1,nqf 2110 qfcoef(iqf,ltp)=0.0d0 c else c if(ltp.gt.nang) then write(iout,*) '***error! you must have skipped an l.' write(iout,*) 'not programed. it will work with some change' call exit(1) endif c c construct qlf c do 140 iq=1,idm1 do 150 i=1,mesh 150 rho(i)=qfunc(i)*bessel(g(iq)*r(i),ltot) call simp(rho,r,rab,a,b,rinner,kkbeta,asum2,idim1) call radin(mesh,rho,rab,asum1,idim1) qlf(iq)=fpi*(asum1-asum2) 140 continue c do 152 iq=1,idm1 152 rho(iq)=(g(iq)*qlf(iq))**2 call simpagn(idm1,g,rho,tqlf,0) c rho(1)=0.0 do 154 i=2,mesh 154 rho(i)=(qfunc(i)/r(i))**2 call simp(rho,r,rab,a,b,rinner,kkbeta,asum2,idim1) call radin(mesh,rho,rab,asum1,idim1) tqr=tpi3*(asum1-asum2) c c now construct the vector targ which depends on qfunc c do 210 i = 2,kkbeta rho(i) = qfunc(i) / r(i)**(ltot+2) 210 continue c call polyn(rinner,val,r(i1),rho(i1),-1) call polyn(rinner,val,r(i1),rho(i1), 0) call polyn(rinner, d1,r(i1),rho(i1),+1) c targ(nqf+1) = val targ(nqf+2) = d1 c c simpsons rule integration up to rtrue c rho(1)=0.d0 do 220 i = 2,kkbeta rho(i) = qfunc(i) * r(i) ** ltot 220 continue c call simp(rho,r,rab,a,b,rinner,kkbeta,asum,idim1) c targ(nqf+3)=asum c do 240 iqf=1,nqf do 235 iq=1,idm1 235 rho(iq)=g(iq)*g(iq)*qn(iq,iqf,ltp)*qlf(iq) call simpagn(idm1,g,rho,asum,0) targ(iqf)=asum 240 continue c c now invert bb for qfcoef c call simeq(idm8,nqf+3,bb(1,1,ltp),targ,targinv) do 245 iqf=1,nqf 245 qfcoef(iqf,ltp)=targinv(iqf) c c print out information if necessary c if ( ifprt .ge.4) then write(iout,*) ' the bb matrix = ' do 246 i=1,nqf 246 write(iout,250) (bb(i,j,ltp),j=1,nqf) write(iout,*) ' the targ vector is ' write(iout,250) (targ(i),i=1,nqf) write (iout,*) ' qfcoef for ltp = ',ltp write (iout,250) (qfcoef(i,ltp),i=1,nqf) write (iout,*) ' the lagrangian mutipliers = ' write (iout,250) (targinv(i),i=nqf+1,nqf+3) c calculate and print residual chisq chisq=tqr-tqlf do 247 i=1,nqf chisq=chisq-2.0*targ(i)*targinv(i) do 247 j=1,nqf chisq=chisq+targinv(i)*bb(i,j,ltp)*targinv(j) 247 continue write (iout,*) ' tqr, tqlf = ',tqr,tqlf write (iout,*) ' residual chi-square = ',chisq endif 250 format (5x,8(1pd13.5)) c endif c c end loop over tot ang momentum components 200 continue c c for nqf=3, use the original algorithm c else if(nqf.eq.3) then c c do 500 ltot = 0,2*(nang-1) c ltp = ltot + 1 ll2 = 2 * ltot + 2 ll3 = 2 * ltot + 3 ll5 = 2 * ltot + 5 ll7 = 2 * ltot + 7 c do 510 i = 2,kkbeta rho(i) = qfunc(i) / r(i)**(ltot+2) 510 continue c call polyn(rinner,val,r(i1),rho(i1),-1) call polyn(rinner,val,r(i1),rho(i1), 0) call polyn(rinner, d1,r(i1),rho(i1),+1) c targ(1) = val targ(2) = d1 c bb(1,1,1) = 1.0d0 bb(2,1,1) = 0.0d0 bb(1,2,1) = rinner**2 bb(2,2,1) = 2.0d0 * rinner bb(1,3,1) = rinner**4 bb(2,3,1) = 4.0d0*rinner**3 c c integral of r**(2*l+2)*rho should be preserved c c bb(3,j,1) = integral of r**(2*l+2)*rho bb(3,1,1) = rinner**ll3/ll3 bb(3,2,1) = rinner**ll5/ll5 bb(3,3,1) = rinner**ll7/ll7 c c simpsons rule integration up to rtrue c rho(1)=0.d0 do 520 i = 2,kkbeta rho(i) = qfunc(i) * r(i) ** ltot 520 continue c call simp(rho,r,rab,a,b,rinner,kkbeta,asum,idim1) c targ(3)=asum c if ( ifprt .ge.4) then write (iout,530) val,d1,asum 530 format ('val,d1,asum=',3f15.5) write (iout,540) 'val ',(bb(1,j,1),j=1,3),targ(1) write (iout,540) 'd1 ',(bb(2,j,1),j=1,3),targ(2) write (iout,540) 'asum',(bb(3,j,1),j=1,3),targ(3) 540 format (1x,a4,4e15.5) endif c call simeq(idm8,3,bb(1,1,1),targ,qfcoef(1,ltp)) c if ( ifprt .ge.4) write (iout,550) 'con ',(qfcoef(i,ltp),i=1,3) 550 format (1x,a4,3e15.5) c c end loop over tot ang momentum components 500 continue c c----------------------------------------------------------------------------- c c otherwise, error c else c write(iout,*) '***error!!!' write(iout,*) 'nqf not appropriate, nqf = ',nqf call exit(1) c endif c c s e t q f u n c f o r l t o t = l m i n c lp = lmin+1 do 300 i = 1,i2 rr = r(i)**2 rho(i)=qfcoef(1,lp) do 295 iqf=2,nqf 295 rho(i)=rho(i)+qfcoef(iqf,lp)*rr**(iqf-1) qfunc(i) = rho(i)*r(i)**(lmin+2) 300 continue c if ( ifprt .ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',r(i1),r(i2),r(i3),r(i4) write (iout,120) 'qf1,qf2,qf3,qf4 =',qfunc(i1),qfunc(i2), + qfunc(i3),qfunc(i4) endif c if (ifprt.gt.0) write (iout,*) ' exiting subroutine psqf1' ifinit=.true. c return end c c----------------------------------------------------------------------- c subroutine psqf2(qfunc,qfcoef,rho,lll, + rinner,mesh,r,rab,a,b,kkbeta,ifprt,nang, + nqf,qtryc,ibfix,rfix,qfix,nfix,ifqopt, + phi1,phi2,ib,jb,idim1,idim3,idim6,idim8) c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c c written from psqf1 by dominic king-smith in november 1992 c modified by kurt stokbro c modified further by d vanderbilt 8/97 c c arrays like 'targ' that are dimensioned idm8 have the c following structure: c 1 to nqf taylor functions 1, r**2, r**4 etc c nqf+1 to nqf+ncon-1 normal constraints (see below) c nqf+ncon optional extra constraint set by nfix c where ncon=4 and ncon=6 for ifqopt=2 and 3 respectively: c ifqopt=2: 3 normal constraints are (value,d1,norm) c ifqopt=3: 5 normal constraints are (value,d1,d2,d3,norm) c where d1,d2,d3 are 1st, 2nd, 3rd derivative c c ------------------------------------------------------ c notes on qfunc and qfcoef: c ------------------------------------------------------ c since q_ij(r) is the product of two orbitals like c psi_{l1,m1}^star * psi_{l2,m2}, it can be decomposed by c total angular momentum l, where l runs over | l1-l2 | , c | l1-l2 | +2 , ... , l1+l2. (l=0 is the only component c needed by the atomic program, which assumes spherical c charge symmetry.) c c recall qfunc(r) = y1(r) * y2(r) where y1 and y2 are the c radial parts of the wave functions defined according to c c psi(r-vec) = (1/r) * y(r) * y_lm(r-hat) . c c for each total angular momentum l, we pseudize qfunc(r) c inside rc as: c c qfunc(r) = r**(l+2) * [ a_1 + a_2*r**2 + a_3*r**4 + ...] c c in such a way as to match qfunc and its c 1'st derivative at rc (ifqopt=2) c 1'st-3'rd derivs at rc (ifqopt=3) c and to preserve c c int_{0}^{r_c} dr r**l * qfunc(r) , c c i.e., to preserve the l'th moment of the charge. also, c depending on nfix,ibfix,etc, there may be one more constraint c to obtain a particular value at a particular radius (as a trick c for avoiding negative charge densities). the rest of c the coeffs are determined by minimizing c c integral (from qtryc to infinity) dq q**2 * qfunc(q)**2 c where c qfunc(q)=int_{0}^{r_c} dr qfunc(r) * j_l(qr) c c where qtryc and the number of coeffs nqf are input. c c qfunc has been set inside rc to correspond to this pseudized c version using the minimal l, namely l = | l1-l2 | (e.g., l=0 c for diagonal elements). the coefficients a_i (i=1,2,3) c are stored in the array qfcoef(i,l+1,j,k) for each l so that c the correctly pseudized versions of qfunc can be reconstructed c for each l. (note that for given l1 and l2, only the values c l = | l1-l2 | , | l1-l2 | +2 , ... , l1+l2 are ever used.) c c ------------------------------------------------------ c c.....in order to avoid passing many unnecessary arrays, use following parameter(idm6=7,nqfx=10,nconx=6,ngx=101,idm8=nqfx+nconx) c c.....logarithmic radial mesh information dimension r(idim1),rab(idim1) c.....q functions and q pseudization coefficients dimension qfunc(idim1),qfcoef(idim8,idim6) c.....q functions and q pseudization cutoffs dimension rinner(idim6),icut(idm6) c.....q functions and l-dependent pseudization arrays for plot dimension qfplot(1000,idm6) c.....extra fixed points to avoid negative densities when pseudizing q_ij dimension ibfix(idim3),rfix(idim3),qfix(idim3) c.....the phi functions dimension phi1(idim1),phi2(idim1) c.....angular momentum information dimension lll(idim3) c.....scratch dimension rho(idim1) dimension targ(idm8),bb(idm8,idm8,idm6),targinv(idm8) dimension qn(ngx,nqfx,idm6),qlf(ngx),g(ngx) logical tinit,tcon c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c c data to be saved to avoid reinitialization c save qn,bb,g,icut data tinit/.false./ c c----------------------------------------------------------------------- c c i n i t i a l i z a t i o n c if ( ifprt .ge. 2) write (iout,'(/3x,a,l1)') + 'entering subroutine psqf2, tinit = ',tinit c c set up qfplot for original unpseudized qfs if (ib .eq. jb) then nline = 1 do 90 i=1,idim1 qfplot(i,nline)=qfunc(i) 90 continue endif c c compute maximum and minimum angular momenta c lmin = iabs ( lll(ib) - lll(jb) ) lmax = iabs ( lll(ib) + lll(jb) ) c if ( ifprt .ge.1) write(iout,'(3x,100a)') ('=',i=1,100) if ( ifprt .ge.1) write(iout,'(3x,6(a,i2))') + ' ib =',ib,' , lll(ib) =',lll(ib),' , jb = ',jb, + ' , lll(jb) =',lll(jb),' , lmin =',lmin,' , lmax =',lmax if ( ifprt .ge.2) write(iout,'(3x,100a)') ('-',i=1,100) c c ========================================================= c set variables that control array index idm8 c ========================================================= c if (ifqopt.eq.2) ncon=4 if (ifqopt.eq.3) ncon=6 c nnorm=nqf+ncon-1 nxtra=nqf+ncon c c ================ c ckeck dimensions c ================ c if ( nqf .le. ncon-1 ) then write(iout,*) '***error in subroutine psqf2' write(iout,*) 'nqf must be .gt.',ncon-1,' but nqf=',nqf call exit(1) endif c if(nqfx.lt.nqf) then write(iout,*) '***error: in subroutine psqf2, nqfx too small' write(iout,*) 'nqfx, nqf = ',nqfx,nqf call exit(1) endif c if(nconx.lt.ncon) then write(iout,*) '***error: in subroutine psqf2, nconx too small' write(iout,*) 'nconx, ncon = ',nconx,ncon call exit(1) endif c if(idm6.lt.2*nang-1) then write(iout,*) '***error: in subroutine psqf2, idm6 too small' write(iout,*) 'idm6, 2*nang-1 = ',idm6,2*nang-1 call exit(1) endif c c set the value of pi c pi=4.0d0*atan(1.0d0) tpi=2.0*pi tpi3=tpi**3 fpi=4.0*pi c c ========================= c possible extra constraint c ========================= c tcon = .false. if ( ib .eq. jb ) then c do 10 ifix = 1,nfix c if ( ib .eq. ibfix(ifix) ) then c if ( ifprt .ge.2) + write (iout,'(5x,a/5x,a,i2/5x,a,f12.7,a,f12.7)') + 'extra constraint q_ij pseudization', + 'constraint for ibeta =',ib, + 'will ensure q_ij passes through',qfix(ifix), + ' at r =',rfix(ifix) c ifcur = ifix c tcon = .true. if ( nqf .le. ncon ) then write(iout,*) '***error subroutine psqf2' write(iout,*) 'nqf .le.',ncon,' will produce singular bb' call exit(1) endif c c set up the last row and column of the bb matrix c corresponding to the extra constraint c do 20 iqf = 1,nqf c iqft = 2*(iqf-1) c bb(nxtra,iqf,1) = rfix(ifix) ** iqft c 20 continue c do 30 iqf = nqf+1,nxtra bb(nxtra,iqf,1) = 0.0d0 30 continue c do 40 iqf = 1,nxtra bb(iqf,nxtra,1) = bb(nxtra,iqf,1) 40 continue c targ(nxtra) = qfix(ifix) / rfix(ifix) ** 2 c endif c 10 continue c endif c c ============================= c one time initialization tasks c ============================= c if (.not. tinit) then c c ========================================== c prepare the uniformly spaced g-vector list c ========================================== c dg=qtryc/dfloat(ngx-1) do 200 iq=1,ngx g(iq)=dble(iq-1)*dg 200 continue c c ==================== c compute the qn array c ==================== c c write(iout,*)'nqf,nqfx,idm6',nqf,nqfx,idm6 call setqn(qn,g,rinner,nqf,2*nang-1,ngx,nqfx,idm6) c write(iout,*)'exit setqn psqf2' c do 250 ltp=1,2*nang-1 c c ========================= c find interpolation points c ========================= c do 100 i2 = kkbeta,1,-1 if ( r(i2) .lt. rinner(ltp) ) goto 110 100 continue c 110 icut(ltp) = i2 c c ===================== c construct matrices bb c ===================== c do 260 iqf=1,nnorm do 260 jqf=1,nnorm bb(iqf,jqf,ltp)=0.0d0 260 continue c c ======================================================== c set parts of bb arising from imposing normal constraints c ======================================================== c do 270 iqf=1,nqf c iqft=2*(iqf-1) c c (1) value at rinner is preserved c bb(nqf+1,iqf,ltp)=rinner(ltp)**iqft c c (2) derivative at rinner is preserved c bb(nqf+2,iqf,ltp)=iqft*rinner(ltp)**(iqft-1) c if (ifqopt.eq.3) then c c (3) 2nd derivative at rinner is preserved c bb(nqf+3,iqf,ltp)=(iqft*(iqft-1))*rinner(ltp)**(iqft-2) c c (4) 3rd derivative at rinner is preserved c bb(nqf+4,iqf,ltp)=(iqft*(iqft-1)*(iqft-2))* + rinner(ltp)**(iqft-3) c endif c (5) integral of r**(2*l+2)*rho should be preserved c lmtp=2*(ltp-1+iqf)+1 bb(nnorm,iqf,ltp)=rinner(ltp)**lmtp/dfloat(lmtp) c c symmetrical part do jqf=nqf+1,nnorm bb(iqf,jqf,ltp)=bb(jqf,iqf,ltp) end do c 270 continue c c set parts of bb to do with smoothing of qfunc c do 280 iqf=1,nqf do 280 jqf=1,nqf do 290 iq=1,ngx rho(iq)=g(iq)*g(iq)*qn(iq,iqf,ltp)*qn(iq,jqf,ltp) 290 continue call simpagn(ngx,g,rho,asum,0) lijtm=2*(ltp-1+iqf+jqf)-1 bb(iqf,jqf,ltp)=tpi3*rinner(ltp)**lijtm/dfloat(lijtm)-asum 280 continue c 250 continue c c end one time initialization events c tinit = .true. endif c c----------------------------------------------------------------------- c c b e g i n l o o p o v e r l t o t c do 300 ltot = 0,2*(nang-1) c ltp = ltot + 1 c c only need to pseudize for lmin,lmin+2 ...lmax-2,lmax c if ( ltot .lt. lmin .or. ltot .gt. lmax .or. + mod(ltot-lmin,2) .eq. 1 ) then c do 310 iqf=1,nqf qfcoef(iqf,ltp)=0.0d0 310 continue c else c c ============= c construct qlf c ============= c c note that c c qlf(g) = 4 * pi * int_{rinner}^{r_c} qfunc(r) j_l(gr) dr c do 320 iq=1,ngx do 330 i=1,mesh rho(i)=qfunc(i)*bessel(g(iq)*r(i),ltot) 330 continue call simp(rho,r,rab,a,b,rinner(ltp),kkbeta,asum2,idim1) call radin(mesh,rho,rab,asum1,idim1) qlf(iq)=fpi*(asum1-asum2) 320 continue c c ============ c compute tqlf c ============ c c note that c c tqlf = int_{0}^{qtryc} g**2 qlf(g)**2 dg c do 340 iq=1,ngx rho(iq)=(g(iq)*qlf(iq))**2 340 continue call simpagn(ngx,g,rho,tqlf,0) c c =========== c compute tqr c =========== c c note that c c tqr = 8 * pi**3 int_{rinner}^{r_c} qfunc(r)**2 / r**2 dr c rho(1)=0.0 do 350 i=2,mesh rho(i)=(qfunc(i)/r(i))**2 350 continue call simp(rho,r,rab,a,b,rinner(ltp),kkbeta,asum2,idim1) call radin(mesh,rho,rab,asum1,idim1) tqr=tpi3*(asum1-asum2) c c ================================================ c construct the vector targ which depends on qfunc c ================================================ c do 400 i = 2,kkbeta rho(i) = qfunc(i) / r(i)**(ltot+2) 400 continue c if (ifqopt.eq.2) then c i2 = icut(ltp) i1 = i2 - 1 i3 = i2 + 1 i4 = i2 + 2 c call polyn(rinner(ltp),val,r(i1),rho(i1),-1) call polyn(rinner(ltp),val,r(i1),rho(i1), 0) call polyn(rinner(ltp), d1,r(i1),rho(i1),+1) c targ(nqf+1) = val targ(nqf+2) = d1 c c possibly save val and d1 for later print out if ( ltot .eq. lmin ) then valmin = val d1min = d1 endif c elseif (ifqopt.eq.3) then c i2 = icut(ltp) i1 = i2 - 5 i3 = i2 + 1 i4 = i2 + 6 c call polyn2(rinner(ltp),val,r(i1),rho(i1),-1) call polyn2(rinner(ltp),val,r(i1),rho(i1), 0) targ(nqf+1)=val do i = 1,ncon-3 call polyn2(rinner(ltp), d1,r(i1),rho(i1),i) targ(nqf+1+i) = d1 enddo c c possibly save val and d1 for later print out if ( ltot .eq. lmin ) then valmin = val d1min = targ(nqf+2) d2min = targ(nqf+3) d3min = targ(nqf+4) endif c else call exit(1) endif c rho(1) = 0.d0 do 410 i = 2,kkbeta rho(i) = qfunc(i) * r(i) ** ltot 410 continue c call simp(rho,r,rab,a,b,rinner(ltp),kkbeta,conae,idim1) c targ(nnorm)=conae c do 420 iqf=1,nqf do 430 iq=1,ngx rho(iq)=g(iq)*g(iq)*qn(iq,iqf,ltp)*qlf(iq) 430 continue call simpagn(ngx,g,rho,asum,0) targ(iqf)=asum 420 continue c c ==================== c invert bb for qfcoef c ==================== c if ( tcon .and. ltot .eq. 0 ) then ncons = ncon else ncons = ncon-1 endif c call simeq(idm8,nqf+ncons,bb(1,1,ltp),targ,targinv) do 435 iqf=1,nqf qfcoef(iqf,ltp)=targinv(iqf) 435 continue c c ================================= c possiple print out of information c ================================= c if ( ifprt .ge.5) then write(iout,'(//,a,i3,a,2i3,//)') 'psqf2: results ltot =', + ltot,' with ibeta,jbeta =',ib,jb write(iout,*) ' the bb matrix = ' do 440 i=1,nqf write(iout,490) (bb(i,j,ltp),j=1,nqf) 440 continue write (iout,*) ' the targ vector is ' write (iout,490) (targ(i),i=1,nqf) write (iout,*) ' qfcoef for ltp = ',ltp write (iout,490) (qfcoef(i,ltp),i=1,nqf) write (iout,*) ' the lagrange mutipliers = ' write (iout,490) (targinv(i),i=nqf+1,nqf+ncons) 490 format (8(1pd13.5)) c endif c c ================================== c calculate and print residual chisq c ================================== c c note: chisq = fpi**2 int_{qtryc}^{inf} g**2 qfunc(g)**2 dg c chisq=tqr-tqlf do 450 i=1,nqf chisq=chisq-2.0*targ(i)*targinv(i) do 460 j=1,nqf chisq=chisq+targinv(i)*bb(i,j,ltp)*targinv(j) 460 continue 450 continue c c ===================================== c compute percentage weight above qtryc c ===================================== c c note: total = fpi**2 int_{0}^{inf} g**2 qfunc(g)**2 dg c total = tpi**3 int_{0}^{r_c} qfunc(r)**2 / r**2 dr c c load rho for r <= r(i2) with qfunc call setqf(qfcoef,rho,r,nqf,ltot,i2,idim1,idim6,idim8) c do 600 ir = i3,mesh rho(ir) = qfunc(ir) 600 continue c rho(1) = 0.0d0 do 610 ir = 2,mesh rho(ir) = ( rho(ir) / r(ir) ) ** 2 610 continue call radin(mesh,rho,rab,total,idim1) total = total * tpi**3 c c ================================================== c check constraint int_{0}^{r_c} r**ltot qfunc(r) dr c ================================================== c c load rho with pseudized form for r <= r(i2) call setqf(qfcoef,rho,r,nqf,ltot,i2,idim1,idim6,idim8) c rho(1)=0.d0 do 630 ir = 2,i2 rho(ir) = rho(ir) * r(ir) ** ltot 630 continue c do 620 ir = i2,kkbeta rho(ir) = qfunc(ir) * r(ir) ** ltot 620 continue c call simp(rho,r,rab,a,b,rinner(ltp),kkbeta,conps,idim1) c c summary of results if ( ifprt .ge.2) then write(iout,'(5x,a,8x,a,13x,a,12x,a,11x,a,13x,a)') + 'ltot','chisq','total','residual','conae','conps' write(iout,'(5x,i3,3x,5(g15.7,3x))') + ltot,chisq,total,chisq/total,conae,conps endif c c-------------------------------------------------------------- c set up qfplot for l-dependent pseudized q functions c if (ib .eq. jb) then nline = nline + 1 call setqf(qfcoef,rho,r,nqf,ltot,i2, + idim1,idim6,idim8) do 91 i=1,i2 qfplot(i,nline)=rho(i) 91 continue c write(iout,*) 'psqf2',nline,qfplot(1,nline) do 92 i=i3,idim1 qfplot(i,nline)=qfunc(i) 92 continue endif c endif c c end loop over tot ang momentum components 300 continue c c print the real-space qfunctions and Fourier transforms c if (ib .eq. jb) then call prqf1(ifprt,qfplot,kkbeta,nline,r,a,b,idim1,idim3) c c fourier transform and print reciprocal-space q functions call fanalq1(r,rab,qfplot,nline,mesh,ifprt, + idim1,idim3) c endif c c----------------------------------------------------------------------- c c s e t q f u n c f o r l t o t = l m i n c lminp = lmin +1 i2 = icut(lminp) if (ifqopt.eq.2) then c i1 = i2 - 1 i3 = i2 + 1 i4 = i2 + 2 c c copy elements up to i2 into qfunc call setqf(qfcoef,qfunc,r,nqf,lmin,i2,idim1,idim6,idim8) c c copy elements up to i4 into rho call setqf(qfcoef, rho,r,nqf,lmin,i4,idim1,idim6,idim8) c c run some checks do ir = i1,i4 rho(ir) = rho(ir) / r(ir)**(lmin+2) end do c lminp = lmin +1 call polyn(rinner(lminp),val,r(i1),rho(i1),-1) call polyn(rinner(lminp),val,r(i1),rho(i1), 0) call polyn(rinner(lminp), d1,r(i1),rho(i1),+1) c if ( ifprt .ge. 2) write (iout,'(5x,a,2f9.5)') + 'target and achieved values of qfunc ',valmin,val, + 'target and achieved values of dqfunc',d1min,d1 c elseif (ifqopt.eq.3) then c i1 = i2 - 5 i3 = i2 + 1 i4 = i2 + 6 c c copy elements up to i2 into qfunc call setqf(qfcoef,qfunc,r,nqf,lmin,i2,idim1,idim6,idim8) c c copy elements up to i4 into rho call setqf(qfcoef, rho,r,nqf,lmin,i4,idim1,idim6,idim8) c c run some checks do ir = i1,i4 rho(ir) = rho(ir) / r(ir)**(lmin+2) end do c lminp = lmin +1 call polyn2(rinner(lminp),val,r(i1),rho(i1),-1) call polyn2(rinner(lminp),val,r(i1),rho(i1), 0) call polyn2(rinner(lminp), d1,r(i1),rho(i1),+1) call polyn2(rinner(lminp), d2,r(i1),rho(i1),+2) call polyn2(rinner(lminp), d3,r(i1),rho(i1),+3) c if ( ifprt .ge. 5) then write(iout,*) 'for lp=',lminp write(iout,'(a,2f11.5)') + 'target and achieved values of qfunc',valmin,val write(iout,'(a,2f11.5)') + 'target and achieved values of dqfunc',d1min,d1 write(iout,'(a,2f11.5)') + 'target and achieved values of d2qfunc',d2min,d2 write(iout,'(a,2f11.5)') + 'target and achieved values of d3qfunc',d3min,d3 endif c else call exit(1) endif c c possible test of the constraint agreement if ( ifprt .ge. 2 .and. tcon ) then c do 560 j3 = 3,i2 if ( rfix(ifcur) .lt. r(j3) ) goto 570 560 continue c 570 continue c call polyn(rfix(ifcur),val,r(j3-2),qfunc(j3-2),-1) call polyn(rfix(ifcur),val,r(j3-2),qfunc(j3-2), 0) c write(iout,'(a,4f9.5)') 'values of r about rfix ', + (r(j),j=j3-2,j3+1) write(iout,'(a,4f9.5)') 'values of q_ij about qfix ', + (qfunc(j),j=j3-2,j3+1) write(iout,'(2a,2f12.7)') 'extra constraint:', + ' target and achieved',qfix(ifcur),val c endif c if ( ifprt .ge. 5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',r(i1),r(i2),r(i3),r(i4) write (iout,120) 'qf1,qf2,qf3,qf4 =',qfunc(i1),qfunc(i2), + qfunc(i3),qfunc(i4) 120 format(a18,4f8.4) endif c return end c----------------------------------------------------------------------- c subroutine pspcor(rspsco,qfcoef,rho,rpcor,mesh,r,rab,a,b, + ifprt,nqf,qtryc,idim1,idim8) c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c c written from psqf2 by Kurt Stokbro February 1996 c c targ are dimensioned nqf+nconx-1, the first nqf of them c are for taylor coeff of 1, r**2, r**4 etc, the last are c (value,deriv,norm) corresponding to 5 lagrange multipliers c c ------------------------------------------------------ c notes on rspsco and qfcoef: c ------------------------------------------------------ c since q_ij(r) is the product of two orbitals like c psi_{l1,m1}^star * psi_{l2,m2}, it can be decomposed by c total angular momentum l, where l runs over | l1-l2 | , c | l1-l2 | +2 , ... , l1+l2. (l=0 is the only component c needed by the atomic program, which assumes spherical c charge symmetry.) c c recall rspsco(r) = y1(r) * y2(r) where y1 and y2 are the c radial parts of the wave functions defined according to c c psi(r-vec) = (1/r) * y(r) * y_lm(r-hat) . c c for each total angular momentum l, we pseudize rspsco(r) c inside rc as: c c rspsco(r) = r**(l+2) * [ a_1 + a_2*r**2 + a_3*r**4 + ...] c c in such a way as to match rspsco and its 1'st derivative at c rc, and to preserve c c int_{0}^{r_c} dr r**l * rspsco(r) , c c i.e., to preserve the l'th moment of the charge. the rest of c the coeffs are determined by minimizing c c integral (from qtryc to infinity) dq q**2 * rspsco(q)**2 c where c rspsco(q)=int_{0}^{r_c} dr rspsco(r) * j_l(qr) c c where qtryc and the number of coeffs nqf are input. c c rspsco has been set inside rc to correspond to this pseudized c version using the minimal l, namely l = | l1-l2 | (e.g., l=0 c for diagonal elements). the coefficients a_i (i=1,2,3) c are stored in the array qfcoef(i,l+1,j,k) for each l so that c the correctly pseudized versions of rspsco can be reconstructed c for each l. (note that for given l1 and l2, only the values c l = | l1-l2 | , | l1-l2 | +2 , ... , l1+l2 are ever used.) c c ------------------------------------------------------ c c.....in order to avoid passing many unnecessary arrays, use following parameter(idm6=1,nqfx=10,nconx=4,ngx=101,idm8=nqfx+nconx) c c.....logarithmic radial mesh information dimension r(idim1),rab(idim1) c.....q functions and q pseudization coefficients dimension rspsco(idim1),qfcoef(idim8,idm6) c.....scratch dimension rho(idim1),rinner(1),icut(1) dimension targ(idm8),bb(idm8,idm8,idm6),targinv(idm8) dimension qn(ngx,nqfx,idm6),qlf(ngx),g(ngx) logical tinit c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c c data to be saved to avoid reinitialization c save qn,bb,g,icut data tinit/.false./ c c----------------------------------------------------------------------- c c i n i t i a l i z a t i o n c if (ifprt .ge. 2 .or. (.not. tinit)) + write (iout,'(/3x,a,l1)') + 'subroutine pspcor: tinit = ',tinit c ctest ? c write(6,*)'nqf,idim1,idim8',nqf,idm1,idm8 c c ========================================================= c this routine is not appropriate for case where nqf .le. nconx-1 c ========================================================= c if ( nqf .le. nconx ) then write(iout,*) '***error in subroutine pspcor' write(iout,*) 'nqf must be .gt. nconx, but nqf=',nqf call exit(1) endif c c ================ c ckeck dimensions c ================ c if(nqfx.lt.nqf) then write(iout,*) '***error: in subroutine pspcor, nqfx too small' write(iout,*) 'nqfx, nqf = ',nqfx,nqf call exit(1) endif c c set the value of pi c pi=4.0d0*atan(1.0d0) tpi=2.0*pi tpi3=tpi**3 fpi=4.0*pi c c c ============================= c one time initialization tasks c ============================= c if (.not. tinit) then c c ========================================== c prepare the uniformly spaced g-vector list c ========================================== c dg=qtryc/dfloat(ngx-1) do 200 iq=1,ngx g(iq)=dble(iq-1)*dg 200 continue c c ==================== c compute the qn array c ==================== c rinner(1) = rpcor ctest? c write(6,*)'nqf,nqfx,idm6',nqf,nqfx,idm6 call setqn(qn,g,rinner,nqf,1,ngx,nqfx,1) c write(6,*)'exit setqn in pspcor' c ltp=1 c c ========================= c find interpolation points c ========================= c do 100 i2 = mesh,1,-1 if ( r(i2) .lt. rinner(ltp) ) goto 110 100 continue c 110 icut(ltp) = i2 c c ===================== c construct matrices bb c ===================== c do 260 iqf=1,nqf+nconx do 260 jqf=1,nqf+nconx bb(iqf,jqf,ltp)=0.0d0 260 continue c c =================================================== c set parts of bb arising from imposing nconx constraints c =================================================== c do 270 iqf=1,nqf c iqft=2*(iqf-1) c c (1) value at rinner is preserved c ctest? c write(6,*)'rinner(ltp)**dble(iqft)',rinner(ltp),dble(iqft) bb(nqf+1,iqf,ltp)=rinner(ltp)**dble(iqft) c c (2) derivative at rinner is preserved c bb(nqf+2,iqf,ltp)=dble(iqft)*rinner(ltp)**(iqft-1) c c (3) 2. derivative at rinner is preserved c bb(nqf+3,iqf,ltp)=dble(iqft*(iqft-1))*rinner(ltp)**(iqft-2) c c (4) 3. derivative at rinner is preserved c bb(nqf+4,iqf,ltp)=dble(iqft*(iqft-1)*(iqft-2))* + rinner(ltp)**(iqft-3) c c c symmetrical part do jqf=nqf+1,nqf+nconx bb(iqf,jqf,ltp)=bb(jqf,iqf,ltp) enddo c 270 continue c c set parts of bb to do with smoothing of rspsco c do 280 iqf=1,nqf do 280 jqf=1,nqf do 290 iq=1,ngx rho(iq)=g(iq)*g(iq)*qn(iq,iqf,ltp)*qn(iq,jqf,ltp) 290 continue call simpagn(ngx,g,rho,asum,0) lijtm=2*(ltp-1+iqf+jqf)-1 if ( ifprt .ge. 5) write(iout,*) 'lijtm=',lijtm bb(iqf,jqf,ltp)= 1 tpi3*rinner(ltp)**dble(lijtm)/dfloat(lijtm)-asum 280 continue c c c end one time initialization events c tinit = .true. endif c c----------------------------------------------------------------------- c c b e g i n l o o p o v e r l t o t c ltot=0 c ltp = ltot + 1 c i2 = icut(ltp) i1 = i2 - 5 i3 = i2 + 1 i4 = i2 + 6 c c c ============= c construct qlf c ============= c c note that c c qlf(g) = 4 * pi * int_{rinner}^{r_c} rspsco(r) j_l(gr) dr c do 320 iq=1,ngx do 330 i=1,mesh rho(i)=rspsco(i)*bessel(g(iq)*r(i),ltot) 330 continue call simp(rho,r,rab,a,b,rinner(ltp),mesh,asum2,idim1) call radin(mesh,rho,rab,asum1,idim1) qlf(iq)=fpi*(asum1-asum2) 320 continue c c ============ c compute tqlf c ============ c c note that c c tqlf = int_{0}^{qtryc} g**2 qlf(g)**2 dg c do 340 iq=1,ngx rho(iq)=(g(iq)*qlf(iq))**2.0 340 continue call simpagn(ngx,g,rho,tqlf,0) c c =========== c compute tqr c =========== c c note that c c tqr = 8 * pi**3 int_{rinner}^{r_c} rspsco(r)**2 / r**2 dr c rho(1)=0.0 do 350 i=2,mesh rho(i)=(rspsco(i)/r(i))**2.0 350 continue call simp(rho,r,rab,a,b,rinner(ltp),mesh,asum2,idim1) call radin(mesh,rho,rab,asum1,idim1) tqr=tpi3*(asum1-asum2) c c ================================================ c construct the vector targ which depends on rspsco c ================================================ c do 400 i = 2,mesh rho(i) = rspsco(i) / r(i)**dble(ltot+2) 400 continue c call polyn2(rinner(ltp),val,r(i1),rho(i1),-1) call polyn2(rinner(ltp),val,r(i1),rho(i1), 0) targ(nqf+1)=val do i = 1,nconx-1 call polyn2(rinner(ltp), d1,r(i1),rho(i1),i) targ(nqf+1+i) = d1 enddo c c c possibly save val and d1 for later print out valmin = val d1min = targ(nqf+2) d2min = targ(nqf+3) d3min = targ(nqf+4) c rho(1)=0.0 do 410 i = 2,mesh rho(i) = rspsco(i) * r(i) ** dble(ltot) 410 continue c c do 420 iqf=1,nqf do 430 iq=1,ngx rho(iq)=g(iq)*g(iq)*qn(iq,iqf,ltp)*qlf(iq) 430 continue call simpagn(ngx,g,rho,asum,0) targ(iqf)=asum 420 continue c c ==================== c invert bb for qfcoef c ==================== c ncon = nconx c call simeq(idm8,nqf+ncon,bb(1,1,ltp),targ,targinv) do 435 iqf=1,nqf qfcoef(iqf,ltp)=targinv(iqf) 435 continue c c ================================= c possiple print out of information c ================================= c if ( ifprt .ge. 3) then c write(iout,*) ' the bb matrix = ' do 440 i=1,nqf write(iout,490) (bb(i,j,ltp),j=1,nqf) 440 continue write (iout,*) ' the targ vector is ' write (iout,490) (targ(i),i=1,nqf) write (iout,*) ' rspscoef for ltp = ',ltp write (iout,490) (qfcoef(i,ltp),i=1,nqf) write (iout,*) ' the lagrange mutipliers = ' write (iout,490) (targinv(i),i=nqf+1,nqf+ncon) 490 format (8(1pd13.5)) c endif c c ================================== c calculate and print residual chisq c ================================== c c note: chisq = fpi**2 int_{qtryc}^{inf} g**2 rspsco(g)**2 dg c chisq=tqr-tqlf do 450 i=1,nqf chisq=chisq-2.0*targ(i)*targinv(i) do 460 j=1,nqf chisq=chisq+targinv(i)*bb(i,j,ltp)*targinv(j) 460 continue 450 continue c c ===================================== c compute percentage weight above qtryc c ===================================== c c note: total = fpi**2 int_{0}^{inf} g**2 rspsco(g)**2 dg c total = tpi**3 int_{0}^{r_c} rspsco(r)**2 / r**2 dr c c load rho for r <= r(i2) with rspsco call setqf(qfcoef,rho,r,nqf,ltot,i2,idim1,1,idim8) c do 600 ir = i3,mesh rho(ir) = rspsco(ir) 600 continue c rho(1) = 0.0d0 do 610 ir = 2,mesh rho(ir) = ( rho(ir) / r(ir) ) ** 2.0 610 continue call radin(mesh,rho,rab,total,idim1) total = total * tpi**3 c c ================================================== c check constraint int_{0}^{r_c} r**ltot rspsco(r) dr c ================================================== c c load rho with pseudized form for r <= r(i2) call setqf(qfcoef,rho,r,nqf,ltot,i2,idim1,1,idim8) c rho(1)=0.0 do 630 ir = 2,i2 rho(ir) = rho(ir) * r(ir) ** dble(ltot) 630 continue c do 620 ir = i2,mesh rho(ir) = rspsco(ir) * r(ir) ** dble(ltot) 620 continue c end loop over tot ang momentum components c c----------------------------------------------------------------------- c c s e t q f u n c f o r l t o t = l m i n c lmin=0 lminp = lmin +1 i2 = icut(lminp) i1 = i2 - 5 i3 = i2 + 1 i4 = i2 + 6 call setqf(qfcoef,rspsco,r,nqf,lmin,i2,idim1,1,idim8) c c run some checks do 520 ir = i1,i4 rho(ir) = rspsco(ir) / r(ir)**dble(lmin+2) 520 continue c lminp = lmin +1 call polyn2(rinner(lminp),val,r(i1),rho(i1),-1) call polyn2(rinner(lminp),val,r(i1),rho(i1), 0) call polyn2(rinner(lminp), d1,r(i1),rho(i1),+1) call polyn2(rinner(lminp), d2,r(i1),rho(i1),+2) call polyn2(rinner(lminp), d3,r(i1),rho(i1),+3) c if ( ifprt .ge. 4) then c write(iout,*) 'For lp=',lminp write(iout,'(a,2f11.5)') + 'target and achieved values of rspsco',valmin,val write(iout,'(a,2f11.5)') + 'target and achieved values of drspsco',d1min,d1 write(iout,'(a,2f11.5)') + 'target and achieved values of d2rspsco',d2min,d2 write(iout,'(a,2f11.5)') + 'target and achieved values of d3rspsco',d3min,d3 c endif c if ( ifprt .ge. 5) then c write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',r(i1),r(i2),r(i3),r(i4) write (iout,120) 'qf1,qf2,qf3,qf4 =',rspsco(i1),rspsco(i2), + rspsco(i3),rspsco(i4) 120 format(a18,4f8.4) c endif c if ( ifprt .ge. 2) write (iout,'(3x,a)') 'leaving pspcor' c return end c c c----------------------------------------------------------------------- c subroutine setqn(qn,g,rinner,nqf,lmaxp,ngx,nqfx,idm6) c c----------------------------------------------------------------------- c c ==================== c compute the qn array c ==================== c c qn(g,n,l) is a moment of a spherical bessel c function between 0 and rinner c c qn(g,n,l) = 4 * pi * int_{0}^{rinner} r**(2n+l) j_l(gr) dr c c for n = 1 to nqf and l = 0 to lmaxp-1 c c programmed as a recurrence relation using power series for j_l c c----------------------------------------------------------------------- c implicit double precision(a-h,o-z) c dimension qn(ngx,nqfx,idm6),g(ngx),rinner(idm6) c ctest? c write(6,*)'qn,g,rinner,nqf',qn,g(1),rinner(1),nqf c write(6,*)'lmaxp,ngx,nqfx,idm6',lmaxp,ngx,nqfx,idm6 c fpi = 16.0d0 * atan(1.0d0) c facl=1.0d0 c c loop over angular momenta do 200 lp=1,lmaxp c ll=lp-1 fll=dble(ll) c note: facl is double factorial facl=facl*(2.0d0*fll+1.0d0) c c loop over the terms in the taylor expansion do 210 in=1,nqf c lp2n=ll+2*in+1 tlnp=dble(2*(ll+in)+1) if(ll.eq.0) then ctest? c write(6,*)'rinner(lp),lp2n',rinner(lp),lp2n c write(6,*)'facl/tlnp',facl,tlnp c qn(1,in,lp)=fpi*rinner(lp)**lp2n/facl/tlnp else qn(1,in,lp)=0.0d0 endif c c loop over points in wavevector grid do 220 ig=2,ngx c qdr=g(ig)*rinner(lp) hqdr2=0.5d0*qdr*qdr c termth=1.0d0 ith=0 c c establish the number of terms in power series 250 continue c ith=ith+1 di=dble(ith) termth=termth*hqdr2/di/(2.0*(fll+di)+1.0) ctest? if(abs(termth).gt.1.0d-20) goto 250 nterm=ith c termth=0.0d0 do 240 ith=nterm,1,-1 di=dble(ith) c termth=hqdr2/di/(2.0*(fll+di)+1.0) + *(1.0/(tlnp+2.0*di)-termth) 240 continue c c finally set the value of qn qn(ig,in,lp)=fpi*rinner(lp)**lp2n*qdr**ll/facl + *(1.0/tlnp-termth) c c close loop over ig 220 continue c close loop over in ctest? c write(6,*)'end loop in' 210 continue c close loop over lp ctest? c write(6,*)'end loop lp' 200 continue c return end c c----------------------------------------------------------------------- c subroutine setqf(qfcoef,rho,r,nqf,ltot,i2,idim1,idim6,idim8) c c----------------------------------------------------------------------- c implicit double precision(a-h,o-z) c c.....logarithmic radial mesh information dimension r(idim1) c.....q pseudization coefficients dimension qfcoef(idim8,idim6) c.....scratch dimension rho(idim1) c c----------------------------------------------------------------------- c c s e t q f u n c f o r t h i s l t o t = l m i n c ltotp = ltot + 1 c do 500 ir = 1,i2 rr = r(ir)**2.0 rho(ir)=qfcoef(1,ltotp) do 510 iqf=2,nqf rho(ir)=rho(ir)+qfcoef(iqf,ltotp)*rr**dble(iqf-1) 510 continue rho(ir) = rho(ir)*r(ir)**(ltot+2) 500 continue c c return end c c----------------------------------------------------------------------------- c subroutine polyn(x,val,r,f,imode) c c routine fits a cubic polynomial through function f c returns 0 to 3rd derivative of f at point x c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c dimension r(4),f(4),c(4) c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c save c c c polynomial definitions: p0(x,y1,y2,y3) = (x-y1) * (x-y2) * (x-y3) p1(x,y1,y2,y3) = (x-y2)*(x-y3) + (x-y1)*(x-y3) + (x-y1)*(x-y2) p2(x,y1,y2,y3) = 2.0d0 * ( (x-y1) + (x-y2) + (x-y3) ) p3(x,y1,y2,y3) = 6.0d0 c c---------------------------------------------------------------------------- c if ( imode .eq. - 1 ) then c initialisation c(1) = f(1) / p0(r(1),r(2),r(3),r(4)) c(2) = f(2) / p0(r(2),r(3),r(4),r(1)) c(3) = f(3) / p0(r(3),r(4),r(1),r(2)) c(4) = f(4) / p0(r(4),r(1),r(2),r(3)) elseif ( imode .eq. 0 ) then c return the zeroth derivative val = c(1)*p0(x,r(2),r(3),r(4)) + c(2)*p0(x,r(3),r(4),r(1)) + + c(3)*p0(x,r(4),r(1),r(2)) + c(4)*p0(x,r(1),r(2),r(3)) elseif ( imode .eq. 1 ) then c return the first derivative val = c(1)*p1(x,r(2),r(3),r(4)) + c(2)*p1(x,r(3),r(4),r(1)) + + c(3)*p1(x,r(4),r(1),r(2)) + c(4)*p1(x,r(1),r(2),r(3)) elseif ( imode .eq. 2 ) then c return the second derivative val = c(1)*p2(x,r(2),r(3),r(4)) + c(2)*p2(x,r(3),r(4),r(1)) + + c(3)*p2(x,r(4),r(1),r(2)) + c(4)*p2(x,r(1),r(2),r(3)) elseif ( imode .eq. 3 ) then c return the third derivative val = c(1)*p3(x,r(2),r(3),r(4)) + c(2)*p3(x,r(3),r(4),r(1)) + + c(3)*p3(x,r(4),r(1),r(2)) + c(4)*p3(x,r(1),r(2),r(3)) endif c return end c---------------------------------------------------------------------------- c c ************************************************************** subroutine pswf1(psi,rho,rcut,lcur,r,rab,con,mesh,ifprt,idim1) c ************************************************************** c this subroutine pseudises the radial wave function to c as a taylor series c c snl(r) = r**(l+1) * exp ( a + b*r**2 + c*r**4 + d*r**6 ) c c matches value and 1st through third derivatives at rcut c c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c c.....logarithmic radial mesh information dimension r(idim1),rab(idim1) c.....wavefunction to be pseudised (psi is really r*psi on input) dimension psi(idim1),rho(idim1) c.....scratch dimension targ(4),con(4),bb(4,4) c index on targ and 2nd index of bb is (val,d1,d2,d3) c index on con, and 1st index of bb taylor coeff of 1, r^2, ... c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c c---------------------------------------------------------------------------- c c i n i t i a l i s a t i o n c if ( ifprt .ge.0) write (iout,'(/3x,a)') 'subroutine pswf1' c c compute the weight under psi rho(1) = 0.0d0 do 10 i = 2,mesh rho(i) = psi(i)**2.0 10 continue call radin(mesh,rho,rab,asum,idim1) c if ( ifprt .ge.5) write (iout,20) asum 20 format (' integral under unpseudised psi is ',f15.10) c lp = lcur + 1 do 30 i=2,mesh psi(i) = psi(i) / r(i)**dble(lp) 30 continue c body of subroutine assumes proportional to psi/r**l c c---------------------------------------------------------------------------- c c f i n d i n t e r p o l a t i o n p o i n t s c do 100 i2 = mesh,1,-1 if ( r(i2) .lt. rcut ) goto 110 100 continue c 110 continue c c check that the value of i2 is reasonable if ( i2 .lt. 3 .or. i2 .gt. mesh-2 ) then write(iout,*) '***error in subroutine pswf1' write(iout,*) 'rcut =',rcut,' generated i2 =',i2, + ' which is illegal with mesh =',mesh call exit(0) endif c i1 = i2 - 1 i3 = i2 + 1 i4 = i2 + 2 c if ( ifprt .ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'ps1,ps2,ps3,ps4 =',(psi(i),i=i1,i4) 120 format(a18,4f8.4) write (iout,130) rcut 130 format(' rcut =',f8.4) endif c c---------------------------------------------------------------------------- c c s e t u p t h e p o l y n o m i a l i n t e r p o l a t i o n c call polyn(rcut,val,r(i1),psi(i1),-1) call polyn(rcut,val,r(i1),psi(i1), 0) call polyn(rcut, d1,r(i1),psi(i1),+1) call polyn(rcut, d2,r(i1),psi(i1),+2) call polyn(rcut, d3,r(i1),psi(i1),+3) c if ( ifprt. ge.5) write (iout,140) val,d1,d2,d3 140 format (' val,d1,d2,d3=',4f15.5) c c---------------------------------------------------------------------------- c c s o l v e s i m u l t a n e o u s e q u a t i o n s c c the basis for the following is that if f=exp(g) then c ln f = g c f' = g' f c f'' = g'' f + g' f' c f''' = g''' f + 2 g'' f' + g' f'' c c NOTE BUG FIX BELOW ! see README and top of runatom.f c for discussion of effect of bug c c val must be positive to prevent error when taking log sig = val / abs(val) c c possible flip of signs val = sig * val d1 = sig * d1 d2 = sig * d2 d3 = sig * d3 c targ(1) = log(val) targ(2) = d1 targ(3) = d2 targ(4) = d3 c bb(1,1) = 1.0d0 bb(2,1) = 0.0d0 bb(3,1) = 0.0d0 bb(4,1) = 0.0d0 bb(1,2) = rcut**2 bb(2,2) = 2.0d0*val*rcut bb(3,2) = 2.0d0*val + 2.0d0*d1*rcut c bug! bb(4,2) = 4.0d0*val + 2.0d0*d2*rcut bb(4,2) = 4.0d0*d1 + 2.0d0*d2*rcut bb(1,3) = rcut**4 bb(2,3) = 4.0d0*val*rcut**3 bb(3,3) = 12.0d0*val*rcut**2 + 4.0d0*d1*rcut**3 bb(4,3) = 24.0d0*val*rcut + 24.0d0*d1*rcut**2 + 4.0d0*d2*rcut**3 bb(1,4) = rcut**6 bb(2,4) = 6.0d0*val*rcut**5 bb(3,4) = 30.0d0*val*rcut**4 + 6.0d0*d1*rcut**5 bb(4,4) = 120.0d0*val*rcut**3 + 60.0d0*d1*rcut**4 + + 6.0d0*d2*rcut**5 c if ( ifprt. ge.5) then write (iout,150) 'val ',(bb(1,j),j=1,4),targ(1) write (iout,150) 'd1 ',(bb(2,j),j=1,4),targ(2) write (iout,150) 'd2 ',(bb(3,j),j=1,4),targ(3) write (iout,150) 'd3 ',(bb(4,j),j=1,4),targ(4) endif 150 format (1x,a4,5e15.5) c call simeq(4,4,bb,targ,con) c if ( ifprt. ge.5) write (iout,160) 'con ',con 160 format (1x,a4,4e15.5) c c---------------------------------------------------------------------------- c c f i l l p s i w i t h p s e u d i s e d f o r m c c save two values psii3=psi(i3) psii4=psi(i4) c c do 200 i=2,i2 do 200 i=2,i4 rr = r(i)**2 psi(i) = con(1) + con(2)*rr + con(3)*rr**2 + con(4)*rr**3 if (psi(i).gt.1.0d2) then write(iout,*) '***error in subroutine pswf1' write(iout,*) 'i,r(i),psi(i)=',i,r(i),psi(i) call exit(1) endif psi(i) = sig * exp( psi(i) ) 200 continue c if ( ifprt. ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'ps1,ps2,ps3,ps4 =',(psi(i),i=i1,i4) endif c c run a self-consistency check call polyn(rcut,val,r(i1),psi(i1),-1) call polyn(rcut,val,r(i1),psi(i1), 0) call polyn(rcut, d1,r(i1),psi(i1),+1) call polyn(rcut, d2,r(i1),psi(i1),+2) call polyn(rcut, d3,r(i1),psi(i1),+3) c if ( ifprt. ge.5) write (iout,140) val,d1,d2,d3 c c restore two values psi(i3)=psii3 psi(i4)=psii4 c c----------------------------------------------------------------------- c c r e f o r m r a d i a l f u n c t i o n a n d w e i g h t c psi(1) = 0.0d0 do 210 i=2,mesh psi(i) = psi(i) * r(i)**dble(lp) 210 continue c rho(1) = 0.0d0 do 220 i = 2,mesh rho(i) = psi(i)**2 220 continue call radin(mesh,rho,rab,asum,idim1) c if ( ifprt .ge.5) write (iout,230) asum 230 format (' integral under pseudised psi is ',f15.10) c if ( ifprt .ge. 0) write (iout,'(3x,a)') 'leaving pswf1' c return end c c c---------------------------------------------------------------------------- c c ************************************************************* subroutine pswf2(psi,rho,rcut,lcur,r,rab,dummy,mesh,ifprt + ,nang,npf,ptryc,idim1,idim8) c ************************************************************* c c modified from pswf by X.-P. Li c c this subroutine pseudises the radial wave function to c as a taylor series c c snl(r) = r**(l+1) * ( a + b*r**2 + c*r**4 + d*r**6 + ... ) c c matches value and 1st through third derivatives at rcut c and the rest of the coeffients are determined by optimizing c the smoothness of the pseudowavefunction c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) parameter(idm8=24,idm1=201,idm9=5) c c.....logarithmic radial mesh information dimension r(idim1),rab(idim1) c.....wavefunction to be pseudised (psi is really r*psi on input) dimension psi(idim1),rho(idim1) c.....scratch dimension targ(idm8),bb(idm8,idm8,idm9),qn(idm1,idm8,idm9) + ,qlf(idm1),g(idm1),targin(idm8),ifin(idm9),con(idm8) + ,dvs(4),dummy(4) c c bug fix by d.v. 5/96: argument 'dummy' of subroutine pswf2 c used to be 'con' but really only allocated as con(4) c now actally allocate scratch array in this subroutine, c as for other arrays like targin c c index on con, and 1st index of bb taylor coeff of 1, r^2, ... c integer stderr logical ifinit c common /files/ stderr,input,iout,ioae,iplot,iologd,iops c save qn,bb,g data ifinit,ifin/.false.,idm9*0/ c c if ( ifprt .ge. 0) write (iout,'(/3x,a)') 'subroutine pswf2' c if(idm8.lt.npf+4) then write(iout,*) '***error: in subroutine pswf2, idm8 too small' write(iout,*) 'idm8, npf+4 = ',idm8,npf+4 call exit(1) endif c if(idm9.lt.nang) then write(iout,*) '***error: in subroutine psqf1, idm9 too small' write(iout,*) 'idm9, nang = ',idm9,nang call exit(1) endif c c c compute the weight under psi rho(1) = 0.0d0 do 10 i = 2,mesh rho(i) = psi(i)**2 10 continue call radin(mesh,rho,rab,asum,idim1) c if ( ifprt .ge.5) write (iout,20) asum 20 format (' integral under unpseudised psi is ',f15.10) c c---------------------------------------------------------------------------- c c f i n d i n t e r p o l a t i o n p o i n t s c do 100 i2 = mesh,1,-1 if ( r(i2) .lt. rcut ) goto 110 100 continue c 110 continue c c check that the value of i2 is reasonable if ( i2 .lt. 6 .or. i2 .gt. mesh-6 ) then write(iout,*) '***error in subroutine pswf2' write(iout,*) 'rcut =',rcut,' generated i2 =',i2, + ' which is illegal with mesh =',mesh call exit(0) endif c i1 = i2 - 5 i3 = i2 + 1 i4 = i2 + 6 c if ( ifprt. ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',r(i1),r(i2),r(i3),r(i4) write (iout,120) 'ps1,ps2,ps3,ps4 =' + ,psi(i1),psi(i2),psi(i3),psi(i4) 120 format(a18,4f8.4) write (iout,130) rcut 130 format(' rcut =',f8.4) endif c lp = lcur + 1 c do 30 i=1,12 c psitt(i) = psi(i+i1-1) / r(i+i1-1)**lp c 30 continue c pi=3.14159265358979 tpi=2.0*pi fpi=4.0*pi qwc=ptryc/2.0 npfw=npf if(npfw.lt.4) npfw=4 c c optimization if npfw gt 4 c if(npfw.ge.4) then c c prepare things beforehand c if(.not.ifinit) then c c prepare g's c dg=ptryc/dfloat(idm1-1) do 2010 iq=1,idm1 2010 g(iq)=qwc+dfloat(iq-1)*dg c ifinit=.true. endif c if(ifin(lp).eq.0) then c c prepare qn's c ll=lp-1 fll=dble(ll) if(ll.eq.0) then facl=1.0 else if(ll.eq.1) then facl=3.0 else if(ll.eq.2) then facl=3.0*5.0 else if(ll.eq.3) then facl=3.0*5.0*7.0 endif do 2050 in=1,npfw lp2n=ll+2*in+1 tlnp=dble(2*(ll+in)+1) c do 2050 iq=1,idm1 qdr=g(iq)*rcut hqdr2=0.5*qdr*qdr c termth=1.0d0 ith=0 2020 continue ith=ith+1 di=dble(ith) termth=termth*hqdr2/di/(2.0*(fll+di)+1.0) if(abs(termth).gt.1.0d-20) goto 2020 nterm=ith c termth=0.0d0 do 2030 ith=nterm,1,-1 di=dble(ith) termth=hqdr2/di/(2.0*(fll+di)+1.0) * *(1.0/(tlnp+2.0*di)-termth) 2030 continue qn(iq,in,lp)=fpi*rcut**dble(lp2n)*qdr**dble(ll)/facl * *(1.0/tlnp-termth) 2050 continue c if((ifprt .ge.1).and.(lp.eq.3)) then write(iout,'(5x,a,i5)') 'notice!! idm1 = ',idm1 write(iout,'(5x,a,i5)') 'qns for l = 2, n = ',npfw write(iout,'(5x,8(1pd13.5))') (qn(iq,npfw,3),iq=1,idm1) endif c c now construct matrices bb c do 2060 iqf=1,npfw+4 do 2060 jqf=1,npfw+4 2060 bb(iqf,jqf,lp)=0.0 c do 2070 iqf=1,npfw iqft=2*(iqf-1)+lp bb(npfw+1,iqf,lp)=rcut**dble(iqft) bb(npfw+2,iqf,lp)=dble(iqft)*rcut**(iqft-1) bb(npfw+3,iqf,lp)=dble(iqft*(iqft-1))*rcut**(iqft-2) bb(npfw+4,iqf,lp)=dble(iqft*(iqft-1)*(iqft-2))*rcut**(iqft-3) c symmetrical part bb(iqf,npfw+1,lp)=bb(npfw+1,iqf,lp) bb(iqf,npfw+2,lp)=bb(npfw+2,iqf,lp) bb(iqf,npfw+3,lp)=bb(npfw+3,iqf,lp) bb(iqf,npfw+4,lp)=bb(npfw+4,iqf,lp) 2070 continue c do 2090 iqf=1,npfw do 2090 jqf=1,npfw do 2080 iq=1,idm1 2080 rho(iq)=g(iq)**4*qn(iq,iqf,lp)*qn(iq,jqf,lp) c2080 rho(iq)=g(iq)*g(iq)*qn(iq,iqf,lp)*qn(iq,jqf,lp) call simpagn(idm1,g,rho,asum,0) bb(iqf,jqf,lp)=asum 2090 continue c ifin(lp)=1 endif c c construct the vector targ c call polyn2(rcut,val,r(i1),psi(i1),-1) call polyn2(rcut,val,r(i1),psi(i1), 0) call polyn2(rcut, d1,r(i1),psi(i1),+1) call polyn2(rcut, d2,r(i1),psi(i1),+2) call polyn2(rcut, d3,r(i1),psi(i1),+3) c if ( ifprt. ge.5) write (iout,140) val,d1,d2,d3 140 format (' val,d1,d2,d3=',4f15.5) c targ(npfw+1) = val targ(npfw+2) = d1 targ(npfw+3) = d2 targ(npfw+4) = d3 c dvs(1)=val dvs(2)=d1 dvs(3)=d2 dvs(4)=d3 c do 2100 iq=1,idm1 2100 qlf(iq)=fpi*expjl(rcut,g(iq),dvs,lp) do 2120 in=1,npfw do 2110 iq=1,idm1 2110 rho(iq)=g(iq)**4*qlf(iq)*qn(iq,in,lp) c2110 rho(iq)=g(iq)*g(iq)*qlf(iq)*qn(iq,in,lp) call simpagn(idm1,g,rho,asum,0) targ(in)=-asum 2120 continue c c---------------------------------------------------------------------------- c c s o l v e s i m u l t a n e o u s e q u a t i o n s c call simeq(idm8,npfw+4,bb(1,1,lp),targ,targin) c do 2130 in=1,npfw 2130 con(in)=targin(in) if ( ifprt. ge.5) then write (iout,*) 'con: ' write (iout,160) (con(in),in=1,npfw) write (iout,*) 'Lagrange multipliers: ' write (iout,160) (targin(in),in=npfw+1,npfw+4) 160 format (8(1pd13.5)) c c residual c chisq=0.0 do 2170 in1=1,npfw do 2170 in2=1,npfw 2170 chisq=chisq+con(in1)*bb(in1,in2,lp)*con(in2) do 2180 in=1,npfw 2180 chisq=chisq-2.0*con(in)*targ(in) do 2190 iq=1,idm1 2190 rho(iq)=(qlf(iq)*g(iq)**2)**2 c2190 rho(iq)=(qlf(iq)*g(iq))**2 call simpagn(idm1,g,rho,asum,0) chisq=chisq+asum write(iout,*) 'residual chi-sq = ',chisq endif c else c write(iout,*) 'npfw not correct, stop' call exit(1) c endif c c---------------------------------------------------------------------------- c c f i l l p s i w i t h p s e u d i s e d f o r m c do 210 i=2,i2 rr = r(i)**2 psi(i)=con(1) do 200 in=2,npfw psi(i) = psi(i)+con(in)*rr**(in-1) 200 continue psi(i)=psi(i)*r(i)**dble(lp) 210 continue c if ( ifprt. ge.5) then write (iout,*) ' i1, i2, i3, i4 =',i1,i2,i3,i4 write (iout,120) ' r1, r2, r3, r4 =',r(i1),r(i2),r(i3),r(i4) write (iout,120) 'ps1,ps2,ps3,ps4 =' + ,psi(i1),psi(i2),psi(i3),psi(i4) endif c c run a self-consistency check c if ( ifprt. ge.5) then c c do 2140 i=1,12 c psitt(i) = psi(i+i1-1) / r(i+i1-1)**lp c2140 continue call polyn2(rcut,val,r(i1),psi(i1),-1) call polyn2(rcut,val,r(i1),psi(i1), 0) call polyn2(rcut, d1,r(i1),psi(i1),+1) call polyn2(rcut, d2,r(i1),psi(i1),+2) call polyn2(rcut, d3,r(i1),psi(i1),+3) write (iout,140) val,d1,d2,d3 c endif c rho(1) = 0.0d0 do 220 i = 2,mesh rho(i) = psi(i)**2 220 continue call radin(mesh,rho,rab,asum,idim1) c if ( ifprt .ge.5) write (iout,230) asum 230 format (' integral under pseudised psi is ',f15.10) c if ( ifprt .ge. 0) write (iout,'(3x,a)') 'leaving pswf2' c return end c c----------------------------------------------------------------------------- c subroutine polyn2(x,val,r,f,imode) c c routine fits a 5th order polynomial through function f using 12 points c returns 0 to 3rd derivative of f at point x c c by X.-P. Li c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c dimension r(12),f(12),bb(6,6),c(6),cc(6) c integer stderr common /files/ stderr,input,iout,ioae,iplot,iologd,iops c save c c c---------------------------------------------------------------------------- c if ( imode .eq. - 1 ) then c c initialize bb do n1=1,6 do n2=1,6 npow=n1+n2-2 if (npow.eq.0) then bb(n1,n2)=12.0d0 else bb(n1,n2)=0.0 do i=1,12 bb(n1,n2)=bb(n1,n2)+(r(i)-x)**npow end do endif end do end do c c initialize cc do n1=1,6 if (n1.eq.1) then cc(n1)=0.0 do i=1,12 cc(n1)=cc(n1)+f(i) end do else cc(n1)=0.0 do i=1,12 cc(n1)=cc(n1)+f(i)*(r(i)-x)**(n1-1) end do endif end do c c solve set of equations call simeq(6,6,bb,cc,c) c elseif ( imode .eq. 0 ) then c return the zeroth derivative val=c(1) c elseif ( imode .eq. 1 ) then c return the first derivative val=c(2) c elseif ( imode .eq. 2 ) then c return the second derivative val=2.0*c(3) c elseif ( imode .eq. 3 ) then c return the third derivative val=6.0*c(4) c endif c return end c double precision function expjl(rc,q,dvs,lp) c c the integral from rc to inf. r*psi(r)*jl(q*r) with large q c and assuming that psi(inf.) is not important c expand to q**(-5) c c by X.-P. Li c c l=3 added by Chris J. Pickard February 1998 c implicit double precision (a-h,o-z) dimension dvs(4) c qdr=q*rc cqdr=cos(qdr) sqdr=sin(qdr) qinv=1.0/q qinv2=qinv*qinv if(lp.eq.1) then expjl=cqdr*qinv2*(dvs(1)-dvs(3)*qinv2) + -sqdr*qinv2*qinv*(dvs(2)-dvs(4)*qinv2) else if(lp.eq.2) then expjl=sqdr*qinv2*(dvs(1)-qinv2*(dvs(3)+dvs(2)/rc + -dvs(1)/(rc*rc))) + +cqdr*qinv2*qinv*(dvs(2)+dvs(1)/rc-qinv2*(dvs(4) + +dvs(3)/rc-2.0*(dvs(2)-dvs(1)/rc)/(rc*rc))) else if(lp.eq.3) then expjl=cqdr*qinv2*(-dvs(1)+qinv2*(dvs(3)+3.0*dvs(2)/rc)) + +sqdr*qinv2*qinv*(dvs(2)+3.0*dvs(1)/rc + -qinv2*(dvs(4)+3.0*(dvs(3)-dvs(2)/rc)/rc)) else if(lp.eq.4) then expjl=dvs(1)*(cqdr*(-6.d0*qinv**3/rc-3.d0*qinv**5/rc**3)+ + sqdr*(9.d0*qinv**4/rc**2-qinv2)) + + dvs(2)*(cqdr*(3.d0*qinv**5/rc**2-qinv**3)+ + sqdr*6.d0*qinv**4/rc) + + dvs(3)*(6.d0*cqdr*qinv**5/rc+sqdr*qinv**4)+ + dvs(4)*cqdr*qinv**5 else stop 99 endif c return end