c---------------------------------------------------------------------------- c program reform c c converts pseudopotential files between formatted and unformatted c forms c c this version is compatible only with pseudopotential program c version 3.0.0 or higher. it assumes logarithmic mesh. c c======================================================================= c c directions for reforming potentials: c c a typical set of input parameters for reform might be: c c 1 # itype c 022-Ti-ca-sp-vgrp.uspp # file to be read from c 022-Ti-ca-sp-vgrp.txt # file to be written to c c The first number 'itype' is a key: c c itype=1 means convert unformatted to formatted c itype=2 means convert formatted to unformatted c c---------------------------------------------------------------------------- c implicit double precision(a-h,o-z) c parameter(idim1=1000,idim3=10,idim4=5,idim5=4,idim6=2*idim5-1) parameter(idim8=20) c c.....radial mesh information dimension r(idim1),rab(idim1) c.....charge densities dimension rsatom(idim1),rspsco(idim1) c.....pseudo quantum numbers, energies occupancies and cutoff radii dimension nnlzps(idim4),eeps(idim4),wwnlps(idim4),rc(idim5) dimension wf(idim1,idim4) c.....beta and q functions and q pseudization coefficients dimension beta(idim1,idim3),qfunc(idim1,idim3,idim3), +qfcoef(idim8,idim6,idim3,idim3),vloc(idim1),vloc0(idim1), +rinner(idim6) c.....indexing for the beta functions dimension nbl0(idim5),nbl(idim5) c.....angular momenta, reference energy, qij and dij of vanderbilt scheme dimension lll(idim3),eee(idim3),iptype(idim3),qqq(idim3,idim3), +ddd0(idim3,idim3),ddd(idim3,idim3) c.....logic to keep track of logarithmic vs herman skillman meshes logical tlog c.....version numbers and date dimension iver(3),idmy(3) c c.....title character*20 title,xctype c.....file name array character*36 infil,outfil c data input,iout,iops /5,6,14/ c c----------------------------------------------------------------------- c c r e a d i n t h e p a r a m e t e r s f o r r u n c c itype (operation type) c infil (input file name) c outfil (output file name) c c write(iout,10) 10 format(' input parameter itype',/, +' ( itype = 1 : unformatted --- > formatted',/, +' itype = 2 : formatted --- > unformatted',/) read (input,*) itype write(iout,*) 'selected itype =',itype write(iout,*) 'name of input file?' read (input,20) infil write(iout,*) 'reading from ',infil write(iout,*) 'name of output file?' read (input,20) outfil write(iout,*) 'writing to ',outfil 20 format (a36) c write(iout,*) 'working on basis that file has version prefix' c c c----------------------------------------------------------------------- c c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c u n f o r m a t t e d t o f o r m a t t e d c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c if ( itype .eq. + 1 ) then c c c r e a d p s e u d o p o t e n t i a l f r o m f i l e c open( unit = iops , file = infil , status = 'old' , + form = 'unformatted' ) c read (iops) (iver(i),i=1,3),(idmy(i),i=1,3) c c specify logarithmic mesh tlog = .true. c read (iops) title,zps,zvps,exftps,nvalps,mesh,etot c c check that the array dimensions are sufficient if (mesh.gt.idim1) call error(' reform',' mesh.gt.idim1',mesh) if (nvalps.gt.idim4) call error(' reform',' nvalps.gt.idim4', + nvalps) c read (iops) (nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps) read (iops) keyps,ifpcor,rinner(1) c if ( iver(1) .ge. 3 ) then read(iops) nang,lloc,eloc,ifqopt,nqf,qtryc endif if (10*iver(1)+iver(2).ge.51) then read (iops) (rinner(i),i=1,nang*2-1) endif c if ( iver(1) .ge. 4 ) then read(iops) irel endif c c set the number of angular momentum terms in q_ij to read in if ( iver(1) .eq. 1 ) then c no distinction between nang and nvalps nang = nvalps c no optimisation of q_ij so 3 term taylor series nqf = 3 nlc = 5 elseif ( iver(1) .eq. 2 ) then c no distinction between nang and nvalps nang = nvalps c no optimisation of q_ij so 3 term taylor series nqf = 3 nlc = 2 * nvalps - 1 else nlc = 2 * nang - 1 endif c if (nang.gt.idim5) call error(' reform',' nang.gt.idim5',nang) if (nqf.gt.idim8) call error(' reform',' nqf.gt.idim8',nqf) c read (iops) (rc(i),i=1,nang) read (iops) nbeta,kkbeta c if (nbeta.gt.idim3) call error(' reform',' nbeta.gt.idim3', + nbeta) c do 100 j=1,nbeta read (iops) lll(j),eee(j),(beta(i,j),i=1,kkbeta) do 100 k=j,nbeta read (iops) ddd0(j,k),ddd(j,k),qqq(j,k), + (qfunc(i,j,k),i=1,kkbeta), + ((qfcoef(i,lp,j,k),i=1,nqf),lp=1,nlc) c c fill transpose parts of matrices (actually, never used) if (k.ne.j) then ddd0(k,j)=ddd0(j,k) ddd (k,j)=ddd (j,k) qqq (k,j)=qqq (j,k) do i=1,kkbeta qfunc(i,k,j)=qfunc(i,j,k) end do do i=1,nqf do lp = 1,nlc qfcoef(i,lp,k,j)=qfcoef(i,lp,j,k) end do end do endif c 100 continue c if (10*iver(1)+iver(2).ge.72) then read (iops) (iptype(j),j=1,nbeta),npf,ptryc endif c read (iops) rcloc,(vloc0(i),i=1,mesh) c c set index arrays nbl and nbl0 do 150 lp = 1,nang nbl(lp) = 0 150 continue do 160 ib = 1,nbeta lp = lll(ib) + 1 nbl(lp) = nbl(lp) + 1 160 continue lmaxps = lll(nbeta) nbl0(1) = 0 do 170 i = 1,lmaxps nbl0(i+1) = nbl0(i) + nbl(i) 170 continue c c possible readin of the frozen core correction if (ifpcor.gt.0) then read (iops) rpcor read (iops) (rspsco(i),i=1,mesh) endif c read (iops) (vloc(i),i=1,mesh) read (iops) (rsatom(i),i=1,mesh) c c possible further read in of information in log mesh case if ( tlog ) then read (iops) (r(i),i=1,mesh) read (iops) (rab(i),i=1,mesh) endif c if (iver(1) .ge. 6) then c nchi = nvales nchi = nvalps if (iver(1) .ge. 7) read (iops) nchi if (nchi.gt.idim4) + call error(' reform',' nchi.gt.idim4',nchi) read (iops) ((wf(i,j), i=1,mesh),j=1,nchi) endif c close (iops) c write (iout,190) 190 format ('1pseudopotential has been read in') c c read in complete c c c w r i t e p s e u d o p o t e n t i a l t o f i l e c open( unit = iops , file = outfil , status = 'new' , + form = 'formatted') c write (iops,890) (iver(i),i=1,3),(idmy(i),i=1,3) write (iops,900) title,zps,zvps,exftps,nvalps,mesh,etot write (iops,910) (nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps) write (iops,920) keyps,ifpcor,rinner(1) if ( iver(1) .ge. 3 ) then write(iops,925) nang,lloc,eloc,ifqopt,nqf,qtryc endif if (10*iver(1)+iver(2).ge.51) then write (iops,930) (rinner(i),i=1,nang*2-1) endif if ( iver(1) .ge. 4 ) then write(iops,926) irel endif write (iops,930) (rc(i),i=1,nang) write (iops,920) nbeta,kkbeta do 300 j=1,nbeta write (iops,940) lll(j),eee(j),(beta(i,j),i=1,kkbeta) do 300 k=j,nbeta write (iops,930) ddd0(j,k),ddd(j,k),qqq(j,k), + (qfunc(i,j,k),i=1,kkbeta), + ((qfcoef(i,lp,j,k),i=1,nqf),lp=1,nlc) 300 continue c if (10*iver(1)+iver(2).ge.72) then write (iops,890) (iptype(j),j=1,nbeta) write (iops,910) npf,ptryc endif c write (iops,930) rcloc,(vloc0(i),i=1,mesh) if (ifpcor.gt.0) then write (iops,930) rpcor write (iops,930) (rspsco(i),i=1,mesh) endif write (iops,930) (vloc(i),i=1,mesh) write (iops,930) (rsatom(i),i=1,mesh) c if ( tlog ) then write (iops,930) (r(i),i=1,mesh) write (iops,930) (rab(i),i=1,mesh) endif if (iver(1) .ge. 6) then if (iver(1) .ge. 7) write (iops,*) nchi write (iops,930) ((wf(i,j), i=1,mesh),j=1,nchi) endif c close (iops) c write (iout,310) 310 format ('1pseudopotential has been written out') c 890 format (6i5) 900 format (a20,3f15.9/2i5,1pe19.11) 910 format (i5,2f15.9) 920 format (2i5,f15.9) 925 format (2i5,f9.5,2i5,f9.5) 926 format (i5) 930 format (1p4e19.11) 940 format (i5/(1p4e19.11)) c endif c c----------------------------------------------------------------------- c c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c f o r m a t t e d t o u n f o r m a t t e d c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c if ( itype .eq. + 2 ) then c c c r e a d p s e u d o p o t e n t i a l f r o m f i l e c open( unit = iops , file = infil , status = 'old' , + form = 'formatted' ) c read (iops,945) (iver(i),i=1,3),(idmy(i),i=1,3) c c specify logarithmic mesh tlog = .true. c read (iops,950) title,zps,zvps,exftps,nvalps,mesh,etot c c check that the array dimensions are sufficient if (mesh.gt.idim1) call error(' reform',' mesh.gt.idim1',mesh) if (nvalps.gt.idim4) call error(' reform',' nvalps.gt.idim4', + nvalps) c read (iops,960) (nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps) read (iops,970) keyps,ifpcor,rinner(1) c if ( iver(1) .ge. 3 ) then read(iops,975) nang,lloc,eloc,ifqopt,nqf,qtryc endif c if (10*iver(1)+iver(2).ge.51) then read (iops,980) (rinner(i),i=1,nang*2-1) endif if ( iver(1) .ge. 4 ) then read(iops,976) irel endif c c set the number of angular momentum terms in q_ij to read in if ( iver(1) .eq. 1 ) then c no distinction between nang and nvalps nang = nvalps c no optimisation of q_ij so 3 term taylor series nqf = 3 nlc = 5 elseif ( iver(1) .eq. 2 ) then c no distinction between nang and nvalps nang = nvalps c no optimisation of q_ij so 3 term taylor series nqf = 3 nlc = 2 * nvalps - 1 else nlc = 2 * nang - 1 endif c if (nang.gt.idim5) call error(' reform',' nang.gt.idim5',nang) if (nqf.gt.idim8) call error(' reform',' nqf.gt.idim8',nqf) c read (iops,980) (rc(i),i=1,nang) read (iops,970) nbeta,kkbeta c if (nbeta.gt.idim3) call error(' reform',' nbeta.gt.idim3', + nbeta) c do 400 j=1,nbeta read (iops,990) lll(j),eee(j),(beta(i,j),i=1,kkbeta) do 400 k=j,nbeta read (iops,980) ddd0(j,k),ddd(j,k),qqq(j,k), + (qfunc(i,j,k),i=1,kkbeta), + ((qfcoef(i,lp,j,k),i=1,nqf),lp=1,nlc) c c fill transpose parts of matrices (actually, never used) if (k.ne.j) then ddd0(k,j)=ddd0(j,k) ddd (k,j)=ddd (j,k) qqq (k,j)=qqq (j,k) do i=1,kkbeta qfunc(i,k,j)=qfunc(i,j,k) end do do i=1,nqf do lp = 1,nlc qfcoef(i,lp,k,j)=qfcoef(i,lp,j,k) end do end do endif c 400 continue c if (10*iver(1)+iver(2).ge.72) then read (iops,945) (iptype(j),j=1,nbeta) read (iops,960) npf,ptryc endif c read (iops,980) rcloc,(vloc0(i),i=1,mesh) c c set index arrays nbl and nbl0 do 450 lp = 1,nang nbl(lp) = 0 450 continue do 460 ib = 1,nbeta lp = lll(ib) + 1 nbl(lp) = nbl(lp) + 1 460 continue lmaxps = lll(nbeta) nbl0(1) = 0 do 470 i = 1,lmaxps nbl0(i+1) = nbl0(i) + nbl(i) 470 continue c c possible readin of the frozen core correction if (ifpcor.gt.0) then read (iops,980) rpcor read (iops,980) (rspsco(i),i=1,mesh) endif c read (iops,980) (vloc(i),i=1,mesh) read (iops,980) (rsatom(i),i=1,mesh) c c possible further read in of information in log mesh case if ( tlog ) then read (iops,980) (r(i),i=1,mesh) read (iops,980) (rab(i),i=1,mesh) endif if (iver(1) .ge. 6) then c nchi = nvales nchi = nvalps if (iver(1) .ge. 7) read (iops,*) nchi read (iops,980) ((wf(i,j), i=1,mesh),j=1,nchi) endif c c close (iops) c write (iout,490) 490 format ('1pseudopotential has been read in') c 945 format (6i5) 950 format (a20,3f15.9/2i5,1pe19.11) 960 format (i5,2f15.9) 970 format (2i5,f15.9) 975 format (2i5,f9.5,2i5,f9.5) 976 format (i5) 980 format (1p4e19.11) 990 format (i5/(1p4e19.11)) c c read in complete c c c w r i t e p s e u d o p o t e n t i a l t o f i l e c open( unit = iops , file = outfil , status = 'new' , + form = 'unformatted') c write (iops) (iver(i),i=1,3),(idmy(i),i=1,3) write (iops) title,zps,zvps,exftps,nvalps,mesh,etot write (iops) (nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps) write (iops) keyps,ifpcor,rinner(1) if ( iver(1) .ge. 3 ) then write(iops) nang,lloc,eloc,ifqopt,nqf,qtryc endif if (10*iver(1)+iver(2).ge.51) then write (iops) (rinner(i),i=1,nang*2-1) endif if ( iver(1) .ge. 4 ) then write(iops) irel endif write (iops) (rc(i),i=1,nang) write (iops) nbeta,kkbeta do 500 j=1,nbeta write (iops) lll(j),eee(j),(beta(i,j),i=1,kkbeta) do 500 k=j,nbeta write (iops) ddd0(j,k),ddd(j,k),qqq(j,k), + (qfunc(i,j,k),i=1,kkbeta), + ((qfcoef(i,lp,j,k),i=1,nqf),lp=1,nlc) 500 continue c if (10*iver(1)+iver(2).ge.72) then write (iops) (iptype(j),j=1,nbeta),npf,ptryc endif c write (iops) rcloc,(vloc0(i),i=1,mesh) if (ifpcor.gt.0) then write (iops) rpcor write (iops) (rspsco(i),i=1,mesh) endif write (iops) (vloc(i),i=1,mesh) write (iops) (rsatom(i),i=1,mesh) c if ( tlog ) then write (iops) (r(i),i=1,mesh) write (iops) (rab(i),i=1,mesh) endif if (iver(1) .ge. 6) then if (iver(1) .ge. 7) write (iops) nchi write (iops) ((wf(i,j), i=1,mesh),j=1,nchi) endif c c close (iops) c write (iout,510) 510 format ('1pseudopotential has been written out') c endif c c----------------------------------------------------------------------- c c p s e u d o p o t e n t i a l r e p o r t c if (exftps.eq. 0.) xctype = ' ceperley-alder' if (exftps.eq.-1.) xctype = ' wigner' if (exftps.eq.-2.) xctype = ' hedin-lundqvist' if (exftps.eq.-3.) xctype = ' gunnarson-lundqvist' if (exftps.eq. 1.) xctype = ' C-A + B88gx + LYPgc' if (exftps.eq. 2.) xctype = ' C-A + B88gx ' if (exftps.eq. 3.) xctype = ' C-A + B88gx + P86gc' if (exftps.eq. 4.) xctype = ' Perdew Wang 1991 ' if (exftps.eq. 5.) xctype = ' PBE - GGA ' c write (iout,1000) (iver(i),i=1,3),idmy(2),idmy(1),idmy(3) 1000 format (/4x,60(1h=)/4x,'| pseudopotential report: version', + i3,'.',i1,'.',i1,' date',i3,'-',i2,'-',i4,2x,'|', + /4x,60(1h-)) write (iout,1010) title,xctype 1010 format (4x,'| ',2a20,' exchange-corr |') write (iout,1020) zps,zvps,exftps 1020 format (4x,'| z =',f5.0,4x,'zv =',f5.0,4x,'exfact =',f10.5, + 13x,'|') write (iout,1030) etot 1030 format (4x,'| ',9x,' ',9x,' etot =',f10.5,13x,'|') write (iout,1040) 1040 format(4x,'| index orbital occupation energy',14x,'|') write (iout,1050) (i,nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps) 1050 format(4x,'|',i5,i11,5x,f10.2,f12.2,15x,'|') if (10*iver(1)+iver(2).ge.51) then write (iout,1061) keyps,ifpcor,rpcor 1061 format(4x,'| keyps =',i2,5x,'ifpcor =',i2,5x,'rpcor =', + f10.5,10x,'|') do 1064 i = 1,2*nang-1 write (iout,1065) rinner(i),i 1064 continue 1065 format(4x,'| rinner =',f10.2,5x,'for L=',i5,22x,'|') else write (iout,1060) keyps,ifpcor,rinner(1) 1060 format(4x,'| keyps =',i2,5x,'ifpcor =',i2,5x,'rinner =', + f10.4,9x,'|') endif c if (keyps.le.2) then c c standard hsc or with vanderbilt smoothness c write (iout,1070) c1070 format(4x,'| hsc-type pseudo generation scheme:',22x,'|') c write (iout,1080) aa c1080 format(4x,'| aa = ',f8.3,41x,'|') c write (iout,1090) key c1090 format(4x,'| key =',3i8,25x,'|') c write (iout,1100) aaa c1100 format(4x,'| aaa =',3f8.3,25x,'|' c + /4x,'| l rc',41x,'|') c write (iout,1110) (i,rc(i),i=1,nvalps) c1110 format(4x,'|',4x,i6,f8.3,40x,'|') c else if (keyps.eq.3) then c c new scheme write (iout,1120) 1120 format(4x,'| new generation scheme:',32x,'|') write (iout,1130) nbeta,kkbeta,rcloc 1130 format(4x,'| nbeta = ',i2,5x,'kkbeta =',i5,5x, + 'rcloc =',f10.4,4x,'|'/ + 4x,'| ibeta l epsilon rcut iptype',17x,'|') do 1140 ib=1,nbeta lp=lll(ib)+1 if (10*iver(1)+iver(2).lt.72) iptype(ib)=-1 c -1 means iptype was not recorded and is unknown write (iout,1150) ib,lll(ib),eee(ib),rc(lp),iptype(ib) 1140 continue 1150 format(4x,'|',6x,i2,5x,i2,4x,2f7.2,6x,i1,18x,'|') c if (10*iver(1)+iver(2).ge.72) then ifip2=0 do ib=1,nbeta if (iptype(ib).eq.2) ifip2=1 end do if (ifip2.eq.1) then write (iout,1155) npf,ptryc 1155 format(4x,'| npf =',i2,' ptryc =',f8.3,29x,'|') endif endif c if ( iver(1) .gt. 2 ) then write(iout,1160) lloc,eloc 1160 format(4x,'| lloc =',i2,' eloc =',f8.3,28x,'|') write(iout,1170) ifqopt,nqf,qtryc 1170 format(4x,'| ifqopt =',i2,' nqf =',i2,' qtryc =',f8.3, + 17x,'|') endif c if ( iver(1) .gt. 3 .and. irel .eq. 2 ) then write(iout,1180) 'koelling-harmon equation' else write(iout,1180) 'schroedinger equation ' endif 1180 format(4x,'|',2x,'all electron calculation used ',a24,2x,'|') c if ( .not. tlog ) then write(iout,2000) 'herman skillman mesh' else write(iout,2000) '**logarithmic mesh**' endif c 2000 format (4x,'|',9x,10('*'),a20,10('*'),9x,'|') c else if (keyps.eq.4) then c c frozen core c write (iout,2500) c2500 format(4x,'| frozen core all-electron case',25x,'|') c endif c write (iout,3000) 3000 format (4x,60(1h=)) c end c c----------------------------------------------------------------------- subroutine error(a,b,n) c----------------------------------------------------------------------- character*(*) a,b write(6,1) a,b,n 1 format(//' program ',a,':',a,'.',8x,i8,8x,'stop') stop end c