program show_psml ! ! Example of extraction of information from PSML file ! This is quite a "low-level" example, in which the ! program queries all the contents and prints the information. ! ! Examples closer to the use cases in electronic-structure ! codes are in preparation ! ! No examples of extraction of "j" or "s" values yet. ! use m_psml use m_getopts implicit none integer, parameter :: dp = selected_real_kind(10,100) character(len=1), dimension(0:4) :: sym = (/ "s", "p", "d", "f", "g" /) type(ps_t) :: ps type(ps_annotation_t) :: annot character(len=200) :: filename logical :: debug, plot, plot_raw character(len=200) :: opt_arg, mflnm, ref_line character(len=10) :: opt_name character(len=10) :: interp_type integer :: interp_quality integer :: nargs, iostat, n_opts, nlabels, iorb, ikb integer :: i, j, l, n, num, nfun, set, seq, ncont real(dp) :: ekb, rc, jval real(dp), allocatable, dimension(:) :: r real(dp) :: Rmax, delta, val integer :: npts, ir character(len=50) :: outname, target_set, set_str integer, allocatable :: setidx(:), idxl(:) integer :: npots, nprojs, nwfs real(dp), allocatable :: raw_r(:), raw_data(:) type(ps_radfunc_t) :: rfunc character(len=10) :: relativity logical :: core_corrections logical :: meta_gga character(len=100) :: creator, date integer :: prov_depth, ninput character(len=100) :: uuid, version, namespace character(len=100) :: name, type real(dp) :: weight integer :: code logical :: has_chlocal logical :: use_effective_range = .true. ! There is a copy of dpnint1 inside the library ! This external declaration is to support calls ! to ps_SetInterpolator external :: interpolate_drh ! Uncomment the following when a new interpolator ! is added to the source tree. The one used previously ! was removed due to license issues. !! external :: interpolate_other ! ! Process options ! n_opts = 0 debug = .false. plot = .false. plot_raw = .false. Rmax = 4.0_dp ! default maximum radius for plots npts = 200 ! default number of points interp_type = "drh" interp_quality = 7 ! this is npoint=2 in OTHER interpolator do call getopts('dR:n:pi:q:tr',opt_name,opt_arg,n_opts,iostat) if (iostat /= 0) exit select case(opt_name) case ('d') debug = .true. case ('t') use_effective_range = .false. case ('p') plot = .true. case ('r') plot_raw = .true. case ('R') read(opt_arg,*) Rmax case ('n') read(opt_arg,*) npts case ('i') read(opt_arg,*) interp_type case ('q') read(opt_arg,*) interp_quality case ('?',':') write(0,*) "Invalid option: ", opt_arg(1:1) write(0,*) "Usage: show_psml [ -d -p -r -R Rmax -n npts ] FILE" write(0,*) " -d for debug output" write(0,*) " -p for output of evaluated data" write(0,*) " -r for output of raw data" write(0,*) " -R Rmax : maximum range of output grid (def: 4 bohr)" write(0,*) " -n npts : number of points of output grid (def: 200)" write(0,*) " -i type : interpolator type (drh) ('other' not implemented)" write(0,*) " -q nq : interpolator quality (npoly for drh)" STOP end select enddo nargs = command_argument_count() nlabels = nargs - n_opts + 1 if (nlabels /= 1) then write(0,*) "Usage: show_psml [ -d -p -r -R Rmax -n npts ] FILE" write(0,*) " -d for debug output" write(0,*) " -p for output of evaluated data" write(0,*) " -r for output of raw data" write(0,*) " -R Rmax : maximum range of output grid (def: 4 bohr)" write(0,*) " -n npts : number of points of output grid (def: 200)" write(0,*) " -i type : interpolator type (drh) ('other' not implemented)" write(0,*) " -q nq : interpolator quality (npoly for drh)" STOP endif call get_command_argument(n_opts,value=filename,status=iostat) if (iostat /= 0) then STOP "Cannot get filename" endif ! print "(a)", "Processing: " // trim(filename) call psml_reader(filename,ps,debug=debug,stat=iostat) if (iostat == -1) then write(0,*) "Cannot open PSML file " // trim(filename) STOP endif call ps_SetEvaluatorOptions(use_effective_range=use_effective_range) #ifdef __NO_PROC_POINTERS__ if (trim(interp_type)=="other") then print "(a)", "No support for other interpolators if no proc pointers" STOP endif call ps_SetEvaluatorOptions(quality_index=interp_quality,debug=debug) print "(a,i3)", "Using DRH (default) interpolator with npoly:", interp_quality #else if (trim(interp_type)=="other") then ! call ps_SetInterpolator(interpolate_other,interp_quality) ! print "(a,i3)", "Using OTHER interpolator with quality index:",interp_quality print "(a)", "Please implement 'other' interpolator first" STOP else if (trim(interp_type)=="drh") then call ps_SetEvaluatorOptions(custom_interpolator=interpolate_drh,& quality_level=interp_quality,debug=debug) print "(a,i3)", "Using DRH interpolator with npoly:",interp_quality else STOP "unknown interpolator" endif #endif call ps_RootAttributes_Get(ps,uuid=uuid,version=version,& namespace=namespace) print "(a,1x,a,1x,a)", "PSML version, uuid:", version, uuid print "(a,1x,a,1x,a)", "PSML namespace:", trim(namespace) prov_depth = ps_Provenance_Depth(ps) do i = prov_depth, 1, -1 call ps_Provenance_Get(ps,i,creator=creator,date=date,annotation=annot,& number_of_input_files=ninput) print "(a,1x,i1,1x,a,1x,a,1x,i1)", "Prov: (depth)", i, trim(creator), trim(date), & ninput call print_annotation(annot) end do ! ! Set up our grid allocate(r(npts)) delta = Rmax / (npts - 1) do ir = 1, npts r(ir) = (ir-1)*delta enddo call ps_PseudoAtomSpec_Get(ps,relativity=relativity,& core_corrections=core_corrections,meta_gga=meta_gga,annotation=annot) call print_annotation(annot) print *, "Relativity: ", trim(relativity) call ps_ExchangeCorrelation_Get(ps,annotation=annot,n_libxc_functionals=nfun) print *, "Exchange-correlation info:" print *, "Number of functionals: ", nfun do i = 1, nfun call ps_LibxcFunctional_Get(ps,i, & name=name,code=code,weight=weight,type=type) print "(a,1x,a,1x,a,1x,i4,f10.6)", "name, type, id, weight:", & trim(name), trim(type), code, weight enddo call print_annotation(annot) ! ! ================================================ ! call ps_Potential_Filter(ps,set=SET_ALL,indexes=setidx,number=npots) print *, "There are ", npots, " semi-local pots" do i = 1, npots call ps_Potential_Get(ps,setidx(i), & l=l,n=n,j=jval,rc=rc,set=set,annotation=annot,func=rfunc) set_str = str_of_set(set) if (set == SET_LJ) then print "(a,1x,i1,a1,1x,f3.1,f10.3)", trim(set_str), n, sym(l), jval, rc else print "(a,1x,i1,a1,f10.3)", trim(set_str), n, sym(l), rc endif if (plot) then write(outname,"(a,i0,a1)") "Vsl." //trim(set_str) // ".", n, sym(l) open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, npts val = ps_Potential_Value(ps,setidx(i),r(ir)) write(4,*) r(ir), val enddo close(4) endif if (plot_raw) then call ps_GetRawData(rfunc,raw_r,raw_data) write(outname,"(a)") trim(outname) // ".raw" open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, size(raw_data) write(4,*) raw_r(ir), raw_data(ir) enddo close(4) endif enddo if (.not. ps_HasLocalPotential(ps)) then print *, "This file does not have a local potential" else call ps_LocalPotential_Get(ps,annotation=annot,type=type,& has_local_charge=has_chlocal,func=rfunc) call print_annotation(annot) print *, "Vlocal type: " // trim(type) print *, "Has Local Charge: ", has_chlocal if (plot) then write(outname,"(a)") "Vlocal" open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, npts val = ps_LocalPotential_Value(ps,r(ir)) write(4,*) r(ir), val enddo close(4) endif if (plot_raw) then call ps_GetRawData(rfunc,raw_r,raw_data) write(outname,"(a)") trim(outname) // ".raw" open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, size(raw_data) write(4,*) raw_r(ir), raw_data(ir) enddo close(4) endif endif if (.not. ps_HasProjectors(ps)) then print *, "This file does not have non-local projectors" else call ps_Projector_Filter(ps,set=SET_ALL,indexes=setidx,number=nprojs) print *, "There are ", nprojs, " projectors" print "(a20,a8,a4,a12)", "set", "l (j)", "seq", "Ekb" do l = 0, 4 call ps_Projector_Filter(ps,indexes_in=setidx,& l=l,indexes=idxl,number=num) do j = 1, num call ps_Projector_Get(ps,idxl(j),ekb=ekb,set=set,seq=seq,j=jval,func=rfunc) set_str = str_of_set(set) if (set == SET_LJ) then print "(a20,i4,1x,f3.1,i4,f12.6)", set_str, l, jval, seq, ekb write(outname,"(a,i0,a1,f3.1,a1,i0)") & "Proj." //trim(set_str) // ".", & l, ".", jval,".", seq else print "(a20,i4,i4,f12.6)", set_str, l, seq, ekb write(outname,"(a,i0,a1,i0)") & "Proj." //trim(set_str) // ".", l, ".", seq endif if (plot) then open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, npts val = ps_Projector_Value(ps,idxl(j),r(ir)) write(4,*) r(ir), val enddo close(4) if (plot_raw) then call ps_GetRawData(rfunc,raw_r,raw_data) write(outname,"(a)") trim(outname) // ".raw" open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, size(raw_data) write(4,*) raw_r(ir), raw_data(ir) enddo close(4) endif endif enddo enddo endif ! Valence charge density call ps_ValenceCharge_Get(ps,annotation=annot) call print_annotation(annot) if (plot) then write(outname,"(a)") "Valence.charge" open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, npts val = ps_ValenceCharge_Value(ps,r(ir)) write(4,*) r(ir), val enddo close(4) endif ! Pseudo-Core charge if (core_corrections) then call ps_CoreCharge_Get(ps,annotation=annot,& nderivs=ncont) call print_annotation(annot) print *, "No of continuous derivs: ", ncont if (plot) then write(outname,"(a)") "Core.charge" open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, npts val = ps_CoreCharge_Value(ps,r(ir)) write(4,"(1p,2e21.13)") r(ir), val enddo close(4) endif endif if (meta_gga) then call ps_ValenceKineticDensity_Get(ps,annotation=annot) call print_annotation(annot) if (plot) then write(outname,"(a)") "Valence.KE.density" open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, npts val = ps_ValenceKineticDensity_Value(ps,r(ir)) write(4,"(1p,2e21.13)") r(ir), val enddo close(4) endif if (core_corrections) then call ps_CoreKineticDensity_Get(ps,annotation=annot,& nderivs=ncont) call print_annotation(annot) print *, "No of continuous derivs: ", ncont if (plot) then write(outname,"(a)") "Core.KE.density" open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, npts val = ps_CoreKineticDensity_Value(ps,r(ir)) write(4,"(1p,2e21.13)") r(ir), val enddo close(4) endif endif endif call ps_PseudoWaveFunctions_Filter(ps,set=SET_ALL,indexes=setidx,number=nwfs) print *, "There are ", nwfs, " pseudo wave-functions" !call display_annotation(ps,"pseudo-wavefunctions") ! do i = 1, nwfs call ps_PseudoWf_Get(ps,setidx(i),l=l,n=n,set=set,j=jval) set_str = str_of_set(set) if (set == SET_LJ) then print "(a,1x,i1,a1,1x,f3.1)", trim(set_str), n, sym(l), jval else print "(a,1x,i1,a1)", trim(set_str), n, sym(l) endif if (plot) then write(outname,"(a,i0,a1)") "Pswf." //trim(set_str) // ".", n, sym(l) open(4,file=outname,form="formatted",status="unknown",& action="write",position="rewind") do ir = 1, npts val = ps_PseudoWf_Value(ps,setidx(i),r(ir)) write(4,*) r(ir), val enddo close(4) endif enddo call ps_destroy(ps) end program show_psml