text:
*----------------------------------------------------------------------*
      subroutine optimize_kernel(pass,lti_cnt,success,fl_opt,
     &     op_info,str_info,orb_info)
*----------------------------------------------------------------------*
*     automatic optimization of the sequence of binary contractions
*     on fl_opt
*
*     andreas, oct. 2012
*
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'ioparam.h'
      include 'stdunit.h'
      include 'mdef_operator_info.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_orbinf.h'
      include 'def_contraction.h'
      include 'def_formula_item.h'
      include 'def_formula_item_list.h'
      include 'ifc_fl_aux.h'

      integer, parameter ::
     &     ntest = 000

      integer, parameter ::
     &    min_replace = 10,
     &    min_scal = 6

      character(len=15), parameter ::
     &     i_am = 'optimize_kernel'

      type(operator_info), intent(inout) ::
     &     op_info
      type(strinf), intent(inout) ::
     &     str_info
      type(orbinf), intent(in) ::
     &     orb_info
      type(formula_item), intent(inout), target ::
     &     fl_opt
      integer, intent(in) ::
     &     pass
      integer, intent(inout) ::
     &     lti_cnt
      logical, intent(out) ::
     &     success

      character(len=len_opname) ::
     &     label_old, label_new, label_op
      logical ::
     &     ok, tra_op, ldummy, skip
      integer ::
     &     idxfl, ngas, nspin, nlist, ilist, nlist2, rank, nj,
     &     cost_intm, size_intm, cost_ref, idx_rpl, iblk_op,
     &     idummy, icmp
      real(8) ::
     &     fact, fact_itf,
     &     cpu, sys, wall, cpu0, sys0, wall0
      type(formula_item), pointer ::
     &     fl_pnt_o, fl_pnt_i,
     &     fl_pnt_mark1, fl_pnt_mark2, fl_pnt_mark3, fl_pnt_mark4
c dbg
     &     ,fl_pnt_watch
c dbg
      type(formula_item_list), target ::
     &     fpl_marks, fpl_marks2
      type(formula_item_list), pointer ::
     &     fpl_marks_pnt, fpl_marks2_pnt
      type(operator), pointer ::
     &     intm_new, op

      integer, pointer ::
     &     occ(:,:,:), rst(:,:,:,:,:,:), svdummy(:)

      logical, external ::
     &     list_cmp
      integer, external ::
     &     cmp_bcontr_spc, ielsum, len_blk_est
c dbg
      logical watch
c
c
      watch = .false.
c dbg

      if (ntest.gt.0) call write_title(lulog,wst_dbg_subr,i_am)
      if (ntest.gt.0) write(lulog,*) ' pass ',pass
c dbg
c      print *,'check at entry'
c      call check_formula_list(fl_opt,op_info)
c dbg

      ngas = orb_info%ngas
      nspin = orb_info%nspin

      success = .false.
      fl_pnt_o => fl_opt
      idxfl = 0

      ! reference for cost estimate
      ! n^6 with at least two virtual or ext, or n^8 in active
      cost_ref = 3000
      cost_ref = max(cost_ref,
     &           len_blk_est((/2,1,0,0,2,1,0,0/),
     &                       1,orb_info))
      cost_ref = max(cost_ref,
     &           len_blk_est((/1,1,1,0,1,1,1,0/),
     &                       1,orb_info))
      cost_ref = max(cost_ref,
     &           len_blk_est((/0,1,2,0,0,1,2,0/),
     &                       1,orb_info))
      cost_ref = max(cost_ref,
     &           len_blk_est((/0,1,3,0,0,1,3,0/),
     &                       1,orb_info))

      outer_loop: do

C        if (ntest.ge.1000)
C     &      call print_form_item(lulog,idxfl,fl_pnt_o,op_info)

        if (fl_pnt_o%command.eq.command_end_of_formula.or.
     &      fl_pnt_o%command.eq.command_set_target_init)
     &    exit outer_loop

        ! for pass=0: look for [NEW INTERMEDIATE]
        if (pass.eq.0) then
          if (fl_pnt_o%command.ne.command_new_intermediate) then
            fl_pnt_o => fl_pnt_o%next
            cycle outer_loop
          end if

          ! none of the constituents must be a _STIN
          if (fl_pnt_o%parent1(1:5).eq.'_STIN'.or.
     &        fl_pnt_o%parent2(1:5).eq.'_STIN') then
            fl_pnt_o => fl_pnt_o%next
            cycle outer_loop
          end if

          ! if the next entry is not a [CONTRACT] we skip here
          ! should be adapted for [REO]
          ! we also skip for traces
          if (fl_pnt_o%next%command.ne.command_bc.and.
     &        fl_pnt_o%next%command.ne.command_bc_reo) then
            fl_pnt_o => fl_pnt_o%next
            cycle outer_loop
          end if
          ! we also skip for traces
          if (fl_pnt_o%next%bcontr%n_operands.ne.2) then
            fl_pnt_o => fl_pnt_o%next
            cycle outer_loop
          end if

          if (ntest.ge.100)  then
            write(lulog,*) 'now examining: '
            call print_form_item(lulog,idxfl,fl_pnt_o,op_info)
            call print_form_item(lulog,idxfl,fl_pnt_o%next,op_info)
          end if

          ! mark this entry and proceed to next entry (should be
          ! the contraction that has the new intm. as result)
          fl_pnt_mark1 => fl_pnt_o
          label_old = fl_pnt_mark1%interm%name
          nj   = fl_pnt_mark1%interm%njoined
          rank = ielsum(fl_pnt_mark1%interm%ica_occ,2*nj)
          fl_pnt_o => fl_pnt_o%next
          fl_pnt_mark2 => fl_pnt_o

          ! check this assumption
          if (trim(fl_pnt_mark2%bcontr%label_res).ne.label_old)
     &      goto 1001
c dbg
c          print *,'OK'
c dbg

        else

          ! other passes: look for contractions
          ! we exclude reorderings for the moment
          if (fl_pnt_o%command.ne.command_bc.and.
     &        fl_pnt_o%command.ne.command_add_bc)!.and.
!     &        fl_pnt_o%command.ne.command_bc_reo.and.
!     &        fl_pnt_o%command.ne.command_add_bc_reo) then
     &     then
            fl_pnt_o => fl_pnt_o%next
            cycle outer_loop
          end if

          ! should be a binary contr.
          ! and result must not be an _STIN
          if (fl_pnt_o%bcontr%n_operands.ne.2.or.
     &        fl_pnt_o%bcontr%label_res(1:5).eq.'_STIN') then
            fl_pnt_o => fl_pnt_o%next
            cycle outer_loop
          end if

          if (ntest.ge.100)  then
c            write(lulog,*) 'now examining: '
c            call print_form_item(lulog,idxfl,fl_pnt_o,op_info)
          end if

          ! mark this entry
          fl_pnt_mark1 => fl_pnt_o

        end if

        fl_pnt_i => fl_pnt_o

        call init_formula_plist(fpl_marks)
        fpl_marks_pnt => fpl_marks
        nlist = 0
        if (pass.gt.0) then
          call init_formula_plist(fpl_marks2)
          fpl_marks2_pnt => fpl_marks2
          nlist2 = 0
        end if

c dbg
c        print *,'starting inner loop '
c        call print_form_item(lulog,idxfl,fl_pnt_i%next,op_info)
c        print *,'form until end:'
c        call print_form_list(lulog,fl_pnt_i%next,op_info)
c dbg

        inner_loop: do

          if (.not.associated(fl_pnt_i%next)) goto 1003
          fl_pnt_i => fl_pnt_i%next
          if (fl_pnt_i%command.eq.command_end_of_formula.or.
     &        fl_pnt_i%command.eq.command_set_target_init)
     &      exit inner_loop

          if (pass.eq.0) then
            ! pass=0: look for [NEW INTERMEDIATE]
            if (fl_pnt_i%command.ne.command_new_intermediate)
     &         cycle inner_loop
            ! interm. must have the same parents:
            if (trim(fl_pnt_i%parent1).ne.trim(fl_pnt_mark1%parent1).or.
     &          trim(fl_pnt_i%parent2).ne.trim(fl_pnt_mark1%parent2))
     &         cycle inner_loop
            ! transposition of parents must be the same
            if (fl_pnt_i%tra1.neqv.fl_pnt_mark1%tra1.or.
     &          fl_pnt_i%tra2.neqv.fl_pnt_mark1%tra2)
     &         cycle inner_loop
            ! the occupation of the intermediates must match
            if (fl_pnt_i%interm%njoined.ne.fl_pnt_mark1%interm%njoined)
     &         cycle inner_loop
            if (.not.list_cmp(fl_pnt_i%interm%ihpvca_occ,
     &                    fl_pnt_mark1%interm%ihpvca_occ,
     &                        fl_pnt_i%interm%njoined*ngastp*2))
     &         cycle inner_loop
            if (.not.list_cmp(fl_pnt_i%interm%igasca_restr,
     &                    fl_pnt_mark1%interm%igasca_restr,
     &                    fl_pnt_mark1%interm%njoined*8*ngas*nspin))
     &         cycle inner_loop
            ! mark entry
            fl_pnt_mark3 => fl_pnt_i
            label_old = fl_pnt_mark3%interm%name
            fl_pnt_i => fl_pnt_i%next
            ! the next entry should be the corresponding contraction
            if (fl_pnt_i%command.ne.command_bc.and.
     &          fl_pnt_i%command.ne.command_bc_reo) goto 1001
            if (trim(fl_pnt_mark3%interm%name).ne.
     &                                 trim(fl_pnt_i%bcontr%label_res))
     &          goto 1001
            ! last check: the contractions must match
            if (fl_pnt_i%command.ne.fl_pnt_mark2%command)
     &         cycle inner_loop
            ! must check that reorderings are equal ... not yet done so:
            if (fl_pnt_i%command.eq.command_bc_reo)
     &         cycle inner_loop
            if (fl_pnt_i%bcontr%n_cnt.ne.fl_pnt_mark2%bcontr%n_cnt)
     &         cycle inner_loop
            if (abs(fl_pnt_i%bcontr%fact-fl_pnt_mark2%bcontr%fact)
     &                                                       .gt.1d-12)
     &         cycle inner_loop
            if (fl_pnt_i%bcontr%iblk_op1.ne.
     &                                    fl_pnt_mark2%bcontr%iblk_op1)
     &         cycle inner_loop
            if (fl_pnt_i%bcontr%iblk_op2.ne.
     &                                    fl_pnt_mark2%bcontr%iblk_op2)
     &         cycle inner_loop
            if (fl_pnt_i%bcontr%tra_op1.neqv.
     &                                     fl_pnt_mark2%bcontr%tra_op1)
     &         cycle inner_loop
            if (fl_pnt_i%bcontr%tra_op2.neqv.
     &                                     fl_pnt_mark2%bcontr%tra_op2)
     &         cycle inner_loop
            if (fl_pnt_i%bcontr%tra_res.neqv.
     &                                     fl_pnt_mark2%bcontr%tra_res)
     &         cycle inner_loop

            if (.not.list_cmp(fl_pnt_i%bcontr%occ_cnt,
     &                    fl_pnt_mark2%bcontr%occ_cnt,
     &                    fl_pnt_mark2%bcontr%n_cnt*ngastp*2 ))
     &         cycle inner_loop
            if (.not.list_cmp(fl_pnt_i%bcontr%rst_cnt,
     &                    fl_pnt_mark2%bcontr%rst_cnt,
     &                    fl_pnt_mark2%bcontr%n_cnt*8*ngas*nspin))
     &         cycle inner_loop
            ! match: put the two entries to list
            if (ntest.ge.100) then
              write(lulog,*) 'found matching intermediate'
              call print_form_item(lulog,idxfl,fl_pnt_mark3,op_info)
              call print_form_item(lulog,idxfl,fl_pnt_i,op_info)
            end if
            nlist = nlist+1
            fpl_marks_pnt%item => fl_pnt_mark3
            call new_formula_plist_entry(fpl_marks_pnt)
            fpl_marks_pnt => fpl_marks_pnt%next
            fpl_marks_pnt%item => fl_pnt_i
            call new_formula_plist_entry(fpl_marks_pnt)
            fpl_marks_pnt => fpl_marks_pnt%next

          else
            ! other passes:
            ! compare binary contraction, must have common result op.
c            call quit(1,i_am,'baustelle')

            if (fl_pnt_i%command.ne.command_bc.and.
     &        fl_pnt_i%command.ne.command_add_bc) !.and.
     &       ! fl_pnt_i%command.ne.command_bc_reo.and.
     &       ! fl_pnt_i%command.ne.command_add_bc_reo)
     &      then
              cycle inner_loop
            end if

            ! TO DO: if we allow REO here, check that these are the same!
            icmp = cmp_bcontr_spc(fl_pnt_mark1%bcontr,fl_pnt_i%bcontr)

            ! the above will only be 1 or 2 if 1 or 2 was not
            ! a short-time intermediate
            ! we must make sure, however, that the other operator
            ! (if it is a short-time intermediate) has the same
            ! occupation
            if (icmp.eq.1) then
              if (.not.list_cmp(fl_pnt_mark1%bcontr%occ_op2,
     &                          fl_pnt_i%bcontr%occ_op2,
     &                          2*ngastp*fl_pnt_i%bcontr%nj_op2).or.
     &            .not.list_cmp(fl_pnt_mark1%bcontr%rst_op2,
     &                          fl_pnt_i%bcontr%rst_op2,
     &                          8*fl_pnt_i%bcontr%ngas*
     &                            fl_pnt_i%bcontr%nspin*
     &                            fl_pnt_i%bcontr%nj_op2)) icmp = 0
            else if (icmp.eq.2) then
              if (.not.list_cmp(fl_pnt_mark1%bcontr%occ_op1,
     &                          fl_pnt_i%bcontr%occ_op1,
     &                          2*ngastp*fl_pnt_i%bcontr%nj_op1).or.
     &            .not.list_cmp(fl_pnt_mark1%bcontr%rst_op1,
     &                          fl_pnt_i%bcontr%rst_op1,
     &                          8*fl_pnt_i%bcontr%ngas*
     &                            fl_pnt_i%bcontr%nspin*
     &                            fl_pnt_i%bcontr%nj_op1)) icmp = 0
            end if
c dbg
c            if (fl_pnt_mark1%bcontr%nj_op1.eq.2.and.
c     &          trim(fl_pnt_mark1%bcontr%label_op1).eq.'_STIN0001') then
c              print *,'comparing'
c              call print_form_item(lulog,idxfl,fl_pnt_mark1,op_info)
c              call print_form_item(lulog,idxfl,fl_pnt_i,op_info)
c              print *,'icmp = ',icmp
c            end if
c dbg
            ! same second op: put to list_m1
            if (icmp.eq.2) then
              nlist = nlist+1
              fpl_marks_pnt%item => fl_pnt_i
              call new_formula_plist_entry(fpl_marks_pnt)
              fpl_marks_pnt => fpl_marks_pnt%next
            end if
            ! same first op: put to list_m2
            if (icmp.eq.1) then
              nlist2 = nlist2+1
              fpl_marks2_pnt%item => fl_pnt_i
              call new_formula_plist_entry(fpl_marks2_pnt)
              fpl_marks2_pnt => fpl_marks2_pnt%next
            end if
          end if

        end do inner_loop

        if (ntest.ge.100) write(lulog,*)'found matching entries: ',nlist
        if (ntest.ge.100.and.pass.gt.0)
     &                   write(lulog,*)'found also:             ',nlist2

        ! pass == 0
        if (pass.eq.0.and.nlist.gt.0) then
          size_intm = len_blk_est(fl_pnt_mark1%interm%ihpvca_occ,
     &                            fl_pnt_mark1%interm%njoined,
     &                            orb_info)
          ! TO BE FIXED FOR N_OPERAND==1
          cost_intm = len_blk_est(fl_pnt_mark2%bcontr%occ_ex1,
     &                            fl_pnt_mark2%bcontr%nj_op1,
     &                            orb_info)
     &              * len_blk_est(fl_pnt_mark2%bcontr%occ_ex2,
     &                            fl_pnt_mark2%bcontr%nj_op2,
     &                            orb_info)
     &              * len_blk_est(fl_pnt_mark2%bcontr%occ_cnt,
     &                            fl_pnt_mark2%bcontr%n_cnt,
     &                            orb_info)

          if (nlist.ge.100.or.cost_intm.gt.cost_ref) then
          ! if the list is non-empty
          ! rename first entry to _LTIN
          lti_cnt = lti_cnt+1
          if (lti_cnt.ge.10000)
     &          call quit(1,i_am,'LTIN counter exceeded')
          write(label_new,'("_LTIN",i4.4)') lti_cnt

          if (ntest.ge.100) write(lulog,*)'new interm: ',trim(label_new)
c dbg
c          print *,trim(label_new),' -> ',nlist,'  size ',size_intm,
c     &          cost_intm
c dbg
          label_old = fl_pnt_mark1%interm%name
          ! allocate new
          call fl_item_rename_intm(fl_pnt_mark1,label_new,
     &                             size_intm.le.lblk_da,orb_info)
C          allocate(intm_new)
C          intm_new%name = label_new
C          call clone_operator(intm_new,fl_pnt_mark1%interm,.false.,
C     &                                                         orb_info)
C          call dealloc_operator(fl_pnt_mark1%interm)
C          deallocate(fl_pnt_mark1%interm)
C          fl_pnt_mark1%interm => intm_new
C          if (size_intm.le.lblk_da) fl_pnt_mark1%incore = 1 ! keep these incore
          ! look for next appearance of _STIN and rename
          fl_pnt_mark2%bcontr%label_res = label_new

          fpl_marks_pnt => fpl_marks
          do ilist = 0, nlist
c dbg
c            print *,'ilist = ',ilist
c dbg
            if (ilist.gt.0) then
              ! the first entry is the definition
              label_old = fpl_marks_pnt%item%interm%name  ! the _STIN name may differ
              fl_pnt_i => fpl_marks_pnt%item%prev ! start search from here
c dbg
c              print *,'original start: '
c              call print_form_list(lulog,fl_pnt_i,op_info)
c dbg
              ! delete this entry from the formula list
c dbg
c              print *,'deleting 1'
c dbg
              call delete_fl_node(fpl_marks_pnt%item)
              deallocate(fpl_marks_pnt%item)
              fpl_marks_pnt => fpl_marks_pnt%next
              ! the second entry is the contraction
              ! remember the pointer to the formula item after that contr.
              !fl_pnt_i => fpl_marks_pnt%item%next
              ! and delete
c dbg
c              print *,'deleting 2'
c dbg
              call delete_fl_node(fpl_marks_pnt%item)
              deallocate(fpl_marks_pnt%item)
              fpl_marks_pnt => fpl_marks_pnt%next
            else
              ! first time: just ...
              !fl_pnt_i => fl_pnt_mark2
              ! better use:
              fl_pnt_i => fl_pnt_mark2
            end if
            ! now look for the first (and only) occurrence of the _STIN
            ! in a contraction
c dbg
c              print *,'replacing ',trim(label_old),' by ',
c     &                             trim(label_new)
c              print *,'start from '
c              !call print_form_item(lulog,idxfl,fl_pnt_i,op_info)
c              call print_form_list(lulog,fl_pnt_i,op_info)
c dbg
            replace_loop: do
              fl_pnt_i => fl_pnt_i%next
c dbg
c              print *,'now at'
c              call print_form_item(lulog,idxfl,fl_pnt_i,op_info)
c dbg
              if (fl_pnt_i%command.eq.command_end_of_formula.or.
     &            fl_pnt_i%command.eq.command_set_target_init)
     &              goto 1002
              ! possibly, the intermediate occurred in def of further
              ! intermediate, so check for that:
              if (fl_pnt_i%command.eq.command_new_intermediate) then
                if (trim(fl_pnt_i%parent1).eq.trim(label_old)) then
                  fl_pnt_i%parent1 = label_new
                end if
                if (trim(fl_pnt_i%parent2).eq.trim(label_old)) then
                  fl_pnt_i%parent2 = label_new
                end if
              end if
              if (fl_pnt_i%command.eq.command_add_bc.or.
     &            fl_pnt_i%command.eq.command_add_bc_reo.or.
     &            fl_pnt_i%command.eq.command_bc.or.
     &            fl_pnt_i%command.eq.command_bc_reo) then
                if (trim(fl_pnt_i%bcontr%label_op1)
     &                          .eq.trim(label_old)) then
                  fl_pnt_i%bcontr%label_op1 = label_new
                  fl_pnt_mark3 => fl_pnt_i
                  exit replace_loop
                end if
                if (trim(fl_pnt_i%bcontr%label_op2)
     &                          .eq.trim(label_old)) then
                  fl_pnt_i%bcontr%label_op2 = label_new
                  fl_pnt_mark3 => fl_pnt_i
                  exit replace_loop
                end if
              end if
            end do replace_loop
          end do ! nlist
          ! -----------------------------------------------
          ! now, the fl_pnt_mark3 should point to the last
          ! formula item in which the intermediate was used
          ! -----------------------------------------------

          ! insert a [DELETE INTERMEDIATE]
c dbg
c          print *,'now inserting a DEL for '//trim(label_new)
c dbg
          call insert_fl_item(fl_pnt_mark3,command_del_intermediate,-1)
          fl_pnt_i%next%label = label_new

          end if
        ! pass > 0
        else if (pass.gt.0.and.(nlist.gt.0.or.nlist2.gt.0)) then
          ! at this line, fl_pnt_mark1 points to the item
          ! where the first contraction is
          ! we will now define the new intermediate, with
          ! fl_pnt_mark2 pointing at the most recent entry
          ! fl_pnt_mark2 should by construction preceed fl_pnt_mark1
          !    (in the list)

          ! if the list is non-empty
          ! process the longer list of list_m1, list_m2
c dbg
c          if (nlist.ge.nlist2) print *,'chose list for op1'
c          if (nlist.lt.nlist2) print *,'chose list for op2'
c dbg
          if (nlist.ge.nlist2) then
            idx_rpl = 1
            fpl_marks_pnt => fpl_marks
          else
            nlist = nlist2
            idx_rpl = 2
            fpl_marks_pnt => fpl_marks2
          end if

c dbg
c          print *,'idx_rpl = ',idx_rpl
c dbg
          ! define new intermediate
          lti_cnt = lti_cnt+1
          if (lti_cnt.gt.10000)
     &          call quit(1,i_am,'LTIN counter exceeded')

          write(label_new,'("_LTIN",i4.4)') lti_cnt
          if (ntest.ge.100) write(lulog,*)'new interm: ',trim(label_new)

c dbg
          !watch = lti_cnt.eq.10
          if (watch) print *,'watching!'
c          print *,'intermediate def'
c dbg
          if (idx_rpl.eq.1) label_op = fl_pnt_mark1%bcontr%label_op1
          if (idx_rpl.eq.2) label_op = fl_pnt_mark1%bcontr%label_op2
          if (idx_rpl.eq.1) iblk_op = fl_pnt_mark1%bcontr%iblk_op1
          if (idx_rpl.eq.2) iblk_op = fl_pnt_mark1%bcontr%iblk_op2
          if (idx_rpl.eq.1) tra_op = fl_pnt_mark1%bcontr%tra_op1
          if (idx_rpl.eq.2) tra_op = fl_pnt_mark1%bcontr%tra_op2
          if (idx_rpl.eq.1) nj = fl_pnt_mark1%bcontr%nj_op1
          if (idx_rpl.eq.2) nj = fl_pnt_mark1%bcontr%nj_op2
          if (idx_rpl.eq.1) occ => fl_pnt_mark1%bcontr%occ_op1
          if (idx_rpl.eq.2) occ => fl_pnt_mark1%bcontr%occ_op2
          if (idx_rpl.eq.1) rst => fl_pnt_mark1%bcontr%rst_op1
          if (idx_rpl.eq.2) rst => fl_pnt_mark1%bcontr%rst_op2
          skip = .false.
          if (label_op(1:5).eq.'_STIN') then
            ! just look for definition and replace
c dbg
            if (watch) print *,'  case STIN first --- '//trim(label_op)
            if (watch) print *,'  look backward from: '
            if (watch)
     &         call print_form_item(lulog,idxfl,fl_pnt_mark1,op_info)
            if (watch)
     &      call print_form_item(lulog,idxfl,fl_pnt_mark1%prev,op_info)
c dbg
            fl_pnt_mark2 => find_fl_item(fl_pnt_mark1,
     &                                command=command_new_intermediate,
     &                                label=label_op,
     &                                backward=.true.)
            if (.not.associated(fl_pnt_mark2))
     &            call quit(1,i_am,'problem')
c dbg
            fl_pnt_watch => fl_pnt_mark2
c dbg
            ! rename here
            label_old = fl_pnt_mark2%interm%name
            ! if one of the parents is a _STIN, we should use a different
            ! operator as parent
            if (fl_pnt_mark2%parent1(1:5).eq.'_STIN'.or.
     &          fl_pnt_mark2%parent2(1:5).eq.'_STIN') then
c dbg
c              print *,'need better parents: '
c dbg
              ! loop over the other terms in the list and take a
              ! better choice
              fpl_marks2_pnt => fpl_marks_pnt
              do
                if (idx_rpl.eq.1) then
                  label_op = fpl_marks2_pnt%item%bcontr%label_op1
                  tra_op = fpl_marks2_pnt%item%bcontr%tra_op1
                else if (idx_rpl.eq.2) then
                  label_op = fpl_marks2_pnt%item%bcontr%label_op2
                  tra_op = fpl_marks2_pnt%item%bcontr%tra_op2
                end if
                ! a real operator? great! take it!
                if (label_op(1:5).ne.'_STIN') then
                  fl_pnt_mark2%parent1 = label_op
                  fl_pnt_mark2%parent2 = '---'
                  fl_pnt_mark2%tra1 = tra_op
                  fl_pnt_mark2%tra2 = .false.
                  exit
                else
c dbg
c                  print *,'resolving another STIN ...'
c dbg
                  ! another STIN ??? find its definition
                  fl_pnt_mark3 => find_fl_item(fpl_marks2_pnt%item,
     &                                command=command_new_intermediate,
     &                                label=label_op,
     &                                backward=.true.)
                  if (.not.associated(fl_pnt_mark3))
     &                    call quit(1,i_am,'problem (2)')
                  ! hopefully this one is defined w/o _STINs
                  if (fl_pnt_mark3%parent1(1:5).ne.'_STIN'.and.
     &                fl_pnt_mark3%parent2(1:5).ne.'_STIN') then
                    fl_pnt_mark2%parent1 = fl_pnt_mark3%parent1
                    fl_pnt_mark2%parent2 = fl_pnt_mark3%parent2
                    fl_pnt_mark2%tra1 = fl_pnt_mark3%tra1
                    fl_pnt_mark2%tra2 = fl_pnt_mark3%tra2
                    exit
                  end if
                end if
                ! the list has always one more entry ...
                if (.not.associated(fpl_marks2_pnt%next%next)) then
c                  call quit(1,i_am,'damn! trapped in _STINs ...')
                   skip = .true.
                   lti_cnt = lti_cnt-1 ! decrease counter again
                   exit
                end if
                fpl_marks2_pnt => fpl_marks2_pnt%next
              end do
            end if
            if (.not.skip) then
c dbg
c              print *,'   renaming'
c dbg
              call fl_item_rename_intm(fl_pnt_mark2,label_new,.false.,
     &                               orb_info)
            else
              nlist = -1
c dbg
c              print *,'giving up ...'
c dbg
            end if

          else
c dbg
            if (watch) print *,'  case OP first --- '//trim(label_op)
c dbg
            ! insert before first occurrence:
            fl_pnt_mark2 => fl_pnt_mark1%prev
            call insert_fl_item(fl_pnt_mark2,
     &                          command_new_intermediate,-1)
            ! store info
            fl_pnt_mark2 => fl_pnt_mark2%next
            call store_def_intm(fl_pnt_mark2,
     &            label_new,occ,rst,nj,1,
     &            label_op,'---',.false.,tra_op,.false.,
     &            orb_info)
c dbg
            fl_pnt_watch => fl_pnt_mark2
c dbg

          end if
          ! define as sum of second/first operators
c dbg
c          if (nlist.gt.-1) print *,'now processing list of terms'
c dbg
          list_loop: do ilist = 0, nlist
c dbg
            if (watch) print *,'  ilist = ',ilist,'/',nlist
c dbg
            if (ilist.eq.0) then
              fl_pnt_mark3 => fl_pnt_mark1
            else
              fl_pnt_mark3 => fpl_marks_pnt%item
              fpl_marks_pnt => fpl_marks_pnt%next
            end if
            if (idx_rpl.eq.1) label_op = fl_pnt_mark3%bcontr%label_op1
            if (idx_rpl.eq.2) label_op = fl_pnt_mark3%bcontr%label_op2
            if (idx_rpl.eq.1) iblk_op = fl_pnt_mark3%bcontr%iblk_op1
            if (idx_rpl.eq.2) iblk_op = fl_pnt_mark3%bcontr%iblk_op2
            if (idx_rpl.eq.1) tra_op = fl_pnt_mark3%bcontr%tra_op1
            if (idx_rpl.eq.2) tra_op = fl_pnt_mark3%bcontr%tra_op2
            if (idx_rpl.eq.1) nj = fl_pnt_mark3%bcontr%nj_op1
            if (idx_rpl.eq.2) nj = fl_pnt_mark3%bcontr%nj_op2
            if (idx_rpl.eq.1) occ => fl_pnt_mark3%bcontr%occ_op1
            if (idx_rpl.eq.2) occ => fl_pnt_mark3%bcontr%occ_op2
            if (idx_rpl.eq.1) rst => fl_pnt_mark3%bcontr%rst_op1
            if (idx_rpl.eq.2) rst => fl_pnt_mark3%bcontr%rst_op2

            fact = fl_pnt_mark3%bcontr%fact
            fact_itf = fl_pnt_mark3%bcontr%fact_itf
c dbg
c            print *,'  fact = ',fact
c dbg
            ! operator: add to intermediate
            if (label_op(1:1).ne.'_') then
c dbg
              if (watch) print *,'  case: OP --- '//trim(label_op)
c dbg
              if (ilist.eq.0) then
                call insert_fl_item(fl_pnt_mark2,
     &                          command_cp_intm,-1)
              else
                call insert_fl_item(fl_pnt_mark2,
     &                          command_add_intm,-1)
              end if
c dbg
              if (watch) print *,'item inserted'
c dbg
              ! dummy map for itf
              allocate(svdummy(3))
              svdummy(1:3) = 0
              fl_pnt_mark2 => fl_pnt_mark2%next
              call store_bc(fl_pnt_mark2,
     &                  fact,fact_itf,
     &                  label_new,label_op,'---',
     &                  1,iblk_op,idummy,
     &                  .false.,tra_op,ldummy,
     &                  nj,nj,0,
     &                  occ,occ,idummy,
     &                  rst,rst,idummy,
     &                  idummy,idummy,idummy,
     &                  idummy,idummy,idummy,0,
     &                  idummy,idummy,
     &                  idummy,idummy,
     &                  svdummy,
     &                  orb_info)

              deallocate(svdummy)
c dbg
              if (watch) print *,' ... and written'
c dbg
            else if (label_op(1:2).eq.'_S') then
c dbg
              if (watch) print *,'  case: _STIN'
c
c dbg
              if (ilist.gt.0) then
                ! _STI: find definition and remove
                fl_pnt_mark4 => find_fl_item(fl_pnt_mark3,
     &                                command=command_new_intermediate,
     &                                label=label_op,
     &                                backward=.true. )
                if (.not.associated(fl_pnt_mark4))
     &             call quit(1,i_am,'mistmist')
                ! delete the definition
                call delete_fl_node(fl_pnt_mark4)
                deallocate(fl_pnt_mark4)
              end if
              ! find the contraction and relink
              fl_pnt_mark4 => find_fl_item(fl_pnt_mark3,
     &                     command_list=(/command_bc,command_bc_reo,
     &                                    command_reorder/),
     &                                nlist = 3,
     &                                label_res=label_op,
     &                                backward=.true. )
              if (.not.associated(fl_pnt_mark4)) then
                 write(lulog,*) 'Did not find defining [CONTR] for '//
     &                       trim(label_op)
                 write(lulog,*) 'I started from ... (search backwards)'
                 call print_form_item(lulog,idxfl,fl_pnt_mark3,op_info)
                 write(lulog,*) 'Preceeding item ...'
                 call print_form_item(lulog,idxfl,
     &                        fl_pnt_mark3%prev,op_info)
                 call quit(1,i_am,'mistmistmist')
              end if

              ! insert contraction at mark 4 behind mark 2
              ! take care to  move all required terms
              call optk_move_bcontr(fl_pnt_mark4,fl_pnt_mark2)
              ! the result must be renamed and the fact must be adapted
              fl_pnt_mark4%bcontr%fact = fl_pnt_mark4%bcontr%fact*fact
              fl_pnt_mark4%bcontr%fact_itf
     %                           = fl_pnt_mark4%bcontr%fact_itf*fact_itf
              fl_pnt_mark4%bcontr%label_res = label_new
              fl_pnt_mark4%bcontr%iblk_res = 1
              ! if not the first term: contract and add:
              if (ilist.gt.0.and.fl_pnt_mark4%command.eq.command_bc)
     &            fl_pnt_mark4%command = command_add_bc
              if (ilist.gt.0.and.fl_pnt_mark4%command.eq.command_bc_reo)
     &            fl_pnt_mark4%command = command_add_bc_reo
              ! if this is a REO rename here, too
              if (fl_pnt_mark4%command.eq.command_bc_reo.or.
     &            fl_pnt_mark4%command.eq.command_add_bc_reo) then
                fl_pnt_mark4%reo%label_in  = label_new
                fl_pnt_mark4%reo%label_out = label_new
              end if

              ! shift mark2 pointer such that further items are appended
              ! after this contraction
c dbg
c             print *,'after inserting: formula starting from old mark 2'
c             call print_form_list(lulog,fl_pnt_mark2,op_info)
c             print *,'--- finito ---'
c dbg
              fl_pnt_mark2 => fl_pnt_mark4
            else if (label_op(1:2).eq.'_L') then
              ! _LTI: find definition and add, make sure that defined before (relink)
               call warn(i_am,'avoiding LTINs presently ...')
               call quit(1,i_am,'There is a bug here ... try '//
     &                          '"routes auto_opt=F"')
               cycle list_loop
            else
               call quit(1,i_am,'illegal operator name: '//
     &                   trim(label_op))
            end if
            ! delete the node
c dbg
            if (watch) then
              print *,'Formula starting from beg. of new inter:'
              call print_form_list(lulog,fl_pnt_watch,op_info)
              print *,'--- finito ---'
            end if
c dbg
            ! the item for ilist==0 is needed below for defining
            ! the new contraction with the present intermediate ...!!
            if (ilist.gt.0) then
c dbg
            if (watch) print *,'  deleting '//
     &                 trim(fl_pnt_mark3%bcontr%label_op1)//'*'//
     &                 trim(fl_pnt_mark3%bcontr%label_op2)
c dbg
              call delete_fl_node(fl_pnt_mark3)
              deallocate(fl_pnt_mark3)
            end if
          end do list_loop

          if (.not.skip) then
c dbg
c          print *,'defining new bin. c. '
c dbg
            ! rename in binary contr.
            if (idx_rpl.eq.1)
     &          fl_pnt_mark1%bcontr%label_op1 = label_new
            if (idx_rpl.eq.2)
     &          fl_pnt_mark1%bcontr%label_op2 = label_new
            if (idx_rpl.eq.1)
     &          fl_pnt_mark1%bcontr%iblk_op1 = 1
            if (idx_rpl.eq.2)
     &          fl_pnt_mark1%bcontr%iblk_op2 = 1
            ! the intermediate is never transposed;
            ! future improvement: it might sometimes be better to
            ! define the transpose of the intermediate, instead.
            if (idx_rpl.eq.1)
     &          fl_pnt_mark1%bcontr%tra_op1 = .false.
            if (idx_rpl.eq.2)
     &          fl_pnt_mark1%bcontr%tra_op2 = .false.
            ! the factors were absorbed into the intermediate
            fl_pnt_mark1%bcontr%fact = 1d0
            fl_pnt_mark1%bcontr%fact_itf = 1d0
            ! delete intm. thereafter
            call insert_fl_item(fl_pnt_mark1,
     &                          command_del_intermediate,-1)
            fl_pnt_mark1%next%label=label_new
            success = .true.
         end if
c dbg
c         print *,'processing done for this term ...'
c dbg

        end if

        call dealloc_formula_plist(fpl_marks)
        if (pass.gt.0) call dealloc_formula_plist(fpl_marks2)

        fl_pnt_o => fl_pnt_o%next

      end do outer_loop

      if (ntest.gt.0) write(lulog,*) 'end of '//trim(i_am)

      return

 1001 write(lulog,*) 'I am still assuming that '//
     &          '[NEW INTERMEDIATE] is followed by its definition'
      write(lulog,*) 'Expected result: '//trim(label_old)
      write(lulog,*) 'But I found: '
      call print_form_item(lulog,idxfl,fl_pnt_o,op_info)
      write(lulog,*) 'the index is possibly wrong, ignore that'
      call quit(1,i_am,'assumption not true? buggy formula?')
 1002 write(lulog,*) '[NEW TARGET] or [END] occured while '//
     &     'trying to replace '//trim(label_old)
      call quit(1,i_am,'buggy formula?')
 1003 write(lulog,*) 'formula ended unexpectedly (without [END])'
      write(lulog,*) 'last item was:'
      call print_form_item(lulog,idxfl,fl_pnt_i,op_info)
      call quit(1,i_am,'buggy formula?')
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
! -*- mode: F90; mode: font-lock -*-
! ------------------------------------------------------------------------------
! $Id$
! ------------------------------------------------------------------------------
! Module pseudo_tm_info
! ------------------------------------------------------------------------------
! Code area 10: Pseudopotentials
! ------------------------------------------------------------------------------

!!****h* Conquest/pseudo_tm_info
!!
!!  NAME
!!   pseudo_tm_info
!!  PURPOSE
!!    Sets up and stores the information of pseudopotentials for
!!   each species of atom.
!!    For using TM's pseudopotetials,
!!
!!  USES
!!   datatypes, GenComms, numbers, splines, global_module
!!  AUTHOR
!!   T. Miyazaki
!!  CREATION DATE
!!   18/06/2002
!!  MODIFICATION HISTORY
!!   2008/02/06 08:32 dave
!!    Changed for output to file not stdout
!!   2008/05/25
!!    Added timers
!!   2011/03/30 19:07 M.Arita
!!    Added statements for P.C.C.
!!   2014/09/15 18:30 lat
!!    fixed call start/stop_timer to timer_module (not timer_stdlocks_module !)
!!   2015/11/09 17:21 dave with TM, NW (Mizuho)
!!    Adding new members to pseudo_info derived type for neutral atom
!!   2017/11/27 15:52 dave with TM, NW
!!    NA projectors added to pseudo_info type
!!   2019/12/03 08:21 dave
!!    Added functional to ps_info structure
!!  SOURCE
!!
module pseudo_tm_info

  use datatypes,              only: double
  use GenComms,               only: cq_abort
  use global_module,          only: area_pseudo, io_lun, flag_pcc_global
  use timer_module,           only: start_timer,stop_timer
  use timer_stdclocks_module, only: tmr_std_allocation

  implicit none

  save
  private
  public :: rad_func, pseudo_info, pseudo, setup_pseudo_info, &
       rad_alloc, rad_dealloc, loc_pot, loc_chg

  ! Choose between representing local part of pseudopotential as a local potential or charge density
  integer :: loc_pot = 1
  integer :: loc_chg = 2

  !type rad_func is same as in TM
  !  see radial.f of SIESTA
  type rad_func
     integer      :: n               ! number of grids along radial direction
     real(double) :: cutoff          ! cutoff radius
     real(double) :: delta           ! step
     real(double), pointer :: f(:)   ! Actual data
     real(double), pointer :: d2(:)  ! Second derivative
     ! Reciprocal space for integrals
     integer :: k_length
     real(double) :: k_delta
     real(double), pointer, dimension(:) :: k_table
  end type rad_func

  type pseudo_info
     character(len=80) :: filename   ! name of *.ion file
     integer :: lmax                 ! maximum of the angular momentum
     integer :: n_pjnl               ! number of projector functions
     logical :: flag_pcc             ! flag for partial core correction
     integer :: z                    ! atomic weight
     integer :: functional           ! Code for functional if set in ion file
     integer :: ps_type              ! Type of pseudopotential (hamann/siesta)
     real(double) :: zval            ! number of valence electrons
     real(double) :: alpha           ! exponent of gaussian for long range term
     ! of local pseudopotential (might not be used)
     real(double) :: prefac          ! (alpha/pi)**1.5 stored for efficiency
     integer :: tm_loc_pot
     type(rad_func) :: vlocal            ! local potential
     type(rad_func) :: chlocal            ! local charge
     ! for Neutral atom potential
     type(rad_func) :: chna               ! neutral atom charge
     type(rad_func) :: vna                ! neutral atom potential
     real(double)   :: eshift             ! shift of vchlocal-zval/r

     type(rad_func) :: chpcc              ! partial core charge
     type(rad_func), pointer :: pjnl(:)   ! size = n_pjnl, projector functions
     integer, pointer :: pjnl_l(:)        ! size = n_pjnl, angular momentum
     integer, pointer :: pjnl_n(:)        ! size = n_pjnl, index for projectors
     real(double), pointer :: pjnl_ekb(:) ! size = n_pjnl, <phi_ln|dV_l|phi_ln>

     ! Neutral atom projectors
     integer :: n_pjna               ! number of NA projector functions
     type(rad_func), pointer :: pjna(:)   ! size = n_pjna, projector functions
     integer, pointer :: pjna_l(:)        ! size = n_pjna, angular momentum
     integer, pointer :: pjna_n(:)        ! size = n_pjna, index for projectors
     real(double), pointer :: pjna_ekb(:) ! size = n_pjna, <phi_ln|V_na|phi_ln>

  end type pseudo_info

  type(pseudo_info), allocatable :: pseudo(:)

  !!***

contains

  !---------------------------------------------------------------
  ! sbrt setup_pseudo_info
  !---------------------------------------------------------------

  !!****f* pseudo_tm_info/setup_pseudo_info *
  !!
  !!  NAME
  !!   setup_pseudo_info
  !!  USAGE
  !!
  !!  PURPOSE
  !!   allocates pseudo(:), whose size is nspecies.
  !!   sets up the name of files *.ion and reads the files.
  !!  INPUTS
  !!   nspecies, species_label
  !!  USES
  !!
  !!  AUTHOR
  !!   T. Miyazaki
  !!  CREATION DATE
  !!   18/06/2002
  !!  MODIFICATION HISTORY
  !!   2008/05/25
  !!    Added timers
  !!   2011/03/30 M.Arita
  !!    Added statements for P.C.C.
  !!   2011/09/16 11:00 dave
  !!    Changed to use n_species and species_label from species_module, and changed atomicrad to atomicnum
  !!   2014/10/12 12:21 lat
  !!    Added possibility of user defined filename for PAO/pseudo *ion
  !!   2017/02/17 15:59 dave
  !!    Changes to allow reading of Hamann's optimised norm-conserving pseudopotentials (remove ABINIT)
  !!   2018/06/19 20:30 nakata
  !!    Add IF with pseudo_type for calling init_rad when ghost atoms are used
  !!   2019/12/09 20:16 dave
  !!    Changes to read valence charge from pseudopotential file and test spin-polarised charge
  !!  SOURCE
  !!
  subroutine setup_pseudo_info

    use datatypes
    use numbers,        ONLY: zero, RD_ERR, two
    use pao_format,     ONLY: pao
    use species_module, ONLY: n_species, species_label, species_file, species_from_files
    use species_module, ONLY: npao_species, nsf_species, type_species, charge, charge_up, charge_dn, mass
    use global_module,  ONLY: iprint_pseudo, flag_Multisite
    use dimens,         ONLY: RadiusSupport, RadiusAtomf, RadiusMS, InvSRange, atomicnum
    use GenComms,       ONLY: inode, ionode, cq_abort, gcopy
    use pseudopotential_common, ONLY: pseudo_type, SIESTA, ABINIT
    use sfc_partitions_module, only: gap_threshold
    use maxima_module,  ONLY: maxnsf

    implicit none

    !local
    integer :: stat, ispecies, l, zeta
    real(double) :: max_rc
    character(len=80) :: filename
    integer :: ii

    if(allocated(pseudo)) then
       if(iprint_pseudo>2.AND.inode==ionode) write(io_lun,fmt='(10x," setup_pseudo_info is skipped because it is already called")')
    else
       call start_timer(tmr_std_allocation)
       allocate(pseudo(n_species),STAT=stat)
       if(stat /= 0) call cq_abort ('allocating pseudo in setup_pseudo_info',stat)
       allocate(pao(n_species),STAT=stat)
       if(stat /= 0) call cq_abort ('allocating pao in setup_pseudo_info',stat)
       call stop_timer(tmr_std_allocation)
       flag_pcc_global = .false.
       max_rc = zero
       do ispecies=1,n_species
          if(pseudo_type==SIESTA) then
             pseudo(ispecies)%tm_loc_pot = loc_chg
          else if(pseudo_type==ABINIT) then
             pseudo(ispecies)%tm_loc_pot = loc_pot
          else
             call cq_abort("Error in pseudopotential type: ",pseudo_type)
          end if
          !**<lat>** 2014/10/12
          if ( species_from_files ) then
             write(filename,'(a,a)') trim(species_file(ispecies))
          else
             write(filename,'(a,a)') trim(species_label(ispecies)),".ion"
          end if
          if(iprint_pseudo>3.AND.inode==ionode) &
               write(io_lun,fmt='(10x,"ispecies = ",i5," file = ",a)') ispecies, filename

          pseudo(ispecies)%filename = filename
          call read_ion_ascii_tmp(pseudo(ispecies),pao(ispecies))

          npao_species(ispecies) = pao(ispecies)%count
          ! Set NSF if not set by user
          if(nsf_species(ispecies)==0) nsf_species(ispecies) = pao(ispecies)%count
          maxnsf = max(maxnsf,nsf_species(ispecies))
          ! Find radius for atom functions
          do l=0,pao(ispecies)%greatest_angmom
             do zeta=1,pao(ispecies)%angmom(l)%n_zeta_in_angmom
                RadiusAtomf(ispecies) = max(RadiusAtomf(ispecies), &
                     pao(ispecies)%angmom(l)%zeta(zeta)%cutoff)
             end do
          end do
          ! Set SF radius if not set by user
          if(abs(RadiusSupport(ispecies))<RD_ERR) RadiusSupport(ispecies) = RadiusAtomf(ispecies)
          if(abs(InvSRange(ispecies))<RD_ERR) InvSRange(ispecies) = RadiusSupport(ispecies)
          if (flag_Multisite) RadiusSupport(ispecies) = RadiusSupport(ispecies) + RadiusMS(ispecies)
          max_rc = max(RadiusSupport(ispecies),max_rc)
          atomicnum(ispecies)    = pseudo(ispecies)%z
          ! Valence charge
          charge(ispecies)       = pseudo(ispecies)%zval
          if(mass(ispecies) < zero) charge(ispecies) = zero
          ! Test spin polarised initialisation
          if (abs(charge_up(ispecies)+charge_dn(ispecies))>RD_ERR) then
             if (abs(charge_up(ispecies)+charge_dn(ispecies)-charge(ispecies))>RD_ERR) &
                  call cq_abort('read_input: sum of number of electrons &
                  &in spin channels is different from total &
                  &number of electrons for this species ', &
                  charge_up(ispecies)+charge_dn(ispecies),charge(ispecies))
          end if
          ! For P.C.C.
          if (pseudo(ispecies)%flag_pcc) flag_pcc_global = .true.
          !For Ghost atoms
          if(type_species(ispecies) < 0) then
             pseudo(ispecies)%zval = zero
             if(pseudo(ispecies)%n_pjnl > 0) then
                do ii=1, pseudo(ispecies)%n_pjnl
                   call init_rad(pseudo(ispecies)%pjnl(ii))
                enddo
             endif
             if(pseudo_type==ABINIT) call init_rad(pseudo(ispecies)%vlocal)
             if(pseudo_type==SIESTA) call init_rad(pseudo(ispecies)%chlocal)
             if(pseudo(ispecies)%flag_pcc) call init_rad(pseudo(ispecies)%chpcc)
             pseudo(ispecies)%n_pjnl      = 0
             pseudo(ispecies)%flag_pcc    = .false.
             pseudo(ispecies)%pjnl_l(:)   = 0
             pseudo(ispecies)%pjnl_n(:)   = 0
             pseudo(ispecies)%pjnl_ekb(:) = zero
             ! pseudo(ispecies)%alpha
             ! pseudo(ispecies)%prefac
          endif

          if(npao_species(ispecies)<nsf_species(ispecies)) &
               call cq_abort("Error ! Fewer PAOs than SFs.  Decrease NumberOfSupports: ", &
               npao_species(ispecies),nsf_species(ispecies))
          !%%! cutoff = zero
          !%%! do l=0,pao(ispecies)%greatest_angmom
          !%%!    if(pao(ispecies)%angmom(l)%n_zeta_in_angmom>0) then
          !%%!       do zeta = 1,pao(ispecies)%angmom(l)%n_zeta_in_angmom
          !%%!          cutoff = max(cutoff,pao(ispecies)%angmom(l)%zeta(zeta)%cutoff)
          !%%!       end do
          !%%!    end if
          !%%! end do
          !%%! if(cutoff>RadiusSupport(ispecies)) then
          !%%!    if(inode==ionode.AND.iprint_pseudo>0) &
          !%%!         write(io_lun,fmt='(10x,"Warning ! Species ",i3," support radius less than PAO radius ",f8.3)') ispecies,cutoff
          !%%!    !RadiusSupport(ispecies) = cutoff
          !%%! endif
       enddo
       ! This isn't ideal, but replaces a read for gap_threshold with default of two*max_rc
       if(abs(gap_threshold)<RD_ERR) gap_threshold = two*max_rc
    endif
    !call gcopy(flag_pcc_global)
    !call gcopy(gap_threshold)
    !call gcopy(maxnsf)
    if (iprint_pseudo>0.AND.inode==ionode .AND.flag_pcc_global) &
         write (io_lun,fmt='(10x,a)') "P.C.C. is taken into account."
    return
  end subroutine setup_pseudo_info
  !!***

  !---------------------------------------------------------------
  ! sbrt alloc_pseudo_info
  !---------------------------------------------------------------

  !!****f* pseudo_tm_info/alloc_pseudo_info *
  !!
  !!  NAME
  !!   alloc_pseudo_info
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Allocates memory for ps_info type
  !!  INPUTS
  !!
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   T.Miyazaki
  !!  CREATION DATE
  !!   May 2003
  !!  MODIFICATION HISTORY
  !!   2008/05/25
  !!    Added timers
  !!  SOURCE
  !!
  subroutine alloc_pseudo_info(ps_info,n)

    use memory_module, ONLY: reg_alloc_mem, type_int, type_dbl

    implicit none

    integer, intent(in) :: n
    type(pseudo_info),intent(out):: ps_info

    integer :: stat

    ps_info%n_pjnl = n
    call start_timer(tmr_std_allocation)
    allocate(ps_info%pjnl(n), STAT=stat)
    if(stat /= 0) call cq_abort('alloc pjnl in alloc_pseudo_info',stat)
    allocate(ps_info%pjnl_l(n), ps_info%pjnl_n(n), STAT=stat)
    if(stat /= 0) call cq_abort('alloc pjnl_l&n in alloc_pseudo_info',stat)
    call reg_alloc_mem(area_pseudo,n,type_int)
    allocate(ps_info%pjnl_ekb(n), STAT=stat)
    if(stat /= 0) call cq_abort('alloc pjnl_ekb in alloc_pseudo_info',stat)
    call reg_alloc_mem(area_pseudo,n,type_dbl)
    call stop_timer(tmr_std_allocation)

    return
  end subroutine alloc_pseudo_info
  !!***

  !---------------------------------------------------------------
  ! sbrt rad_alloc : from TM
  !---------------------------------------------------------------

  !!****f* pseudo_tm_info/rad_alloc *
  !!
  !!  NAME
  !!   rad_alloc
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Sets the 'size' n of the arrays and allocates f and d2.
  !!
  !!  INPUTS
  !!
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   T.Miyazaki
  !!  CREATION DATE
  !!
  !!  MODIFICATION HISTORY
  !!   2008/05/25
  !!    Added timers
  !!  SOURCE
  !!
  subroutine rad_alloc(func,n)

    use memory_module, ONLY: reg_alloc_mem, type_dbl

    implicit none

    type(rad_func), intent(out)    :: func
    integer, intent(in)        :: n
    integer :: stat

    func%n = n
    call start_timer(tmr_std_allocation)
    allocate(func%f(n),func%d2(n),STAT=stat)
    if(stat /= 0) call cq_abort('allocating radial func in rad_alloc',stat)
    call reg_alloc_mem(area_pseudo,2*n,type_dbl)
    call stop_timer(tmr_std_allocation)
    return
  end subroutine rad_alloc
  !!***

  !---------------------------------------------------------------
  ! sbrt rad_dealloc
  !---------------------------------------------------------------

  !!****f* pseudo_tm_info/rad_dealloc *
  !!
  !!  NAME
  !!   rad_dealloc
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Deallocates f and d2
  !!  INPUTS
  !!
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   T.Miyazaki
  !!  CREATION DATE
  !!
  !!  MODIFICATION HISTORY
  !!   2008/05/25
  !!    Added timers
  !!  SOURCE
  !!
  subroutine rad_dealloc(func)

    use memory_module, ONLY: reg_dealloc_mem, type_dbl

    implicit none

    type(rad_func), intent(inout)    :: func

    integer :: stat

    call start_timer(tmr_std_allocation)
    deallocate(func%f,func%d2,STAT=stat)
    if(stat /= 0) call cq_abort('deallocating radial func in rad_alloc',stat)
    call reg_dealloc_mem(area_pseudo,2*func%n,type_dbl)
    nullify(func%f, func%d2)
    call stop_timer(tmr_std_allocation)
    return
  end subroutine rad_dealloc
  !!***

  !---------------------------------------------------------------
  ! sbrt init_rad: from TM
  !---------------------------------------------------------------

  !!****f* pseudo_tm_info/init_rad *
  !!
  !!  NAME
  !!   init_rad
  !!  USAGE
  !!
  !!  PURPOSE
  !!   initilise all members of rad_func
  !!
  !!  INPUTS
  !!
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   T.Miyazaki
  !!  CREATION DATE
  !!
  !!  MODIFICATION HISTORY
  !!    2012/06/17 L.Tong
  !!    - Changed func to intent(inout) from intent(out). If set to out,
  !!      then all members of func upon entering the subroutine
  !!      becomes undefined! This means f and d2 becomes disassociated
  !!      and non-allocated.
  !!  SOURCE
  !!
  subroutine init_rad(func)

    use numbers, only: zero
    implicit none

    type(rad_func), intent(inout)    :: func

    func%n = 0
    func%cutoff = zero
    func%delta = zero
    func%f(:) = zero
    func%d2(:) = zero
    return
  end subroutine init_rad
  !!***

  !---------------------------------------------------------------
  ! sbrt radial_read_ascii : from SIESTA, with small changes
  !---------------------------------------------------------------

  !!****f* pseudo_tm_info/radial_read_ascii *
  !!
  !!  NAME
  !!   radial_read_ascii
  !!  USAGE
  !!
  !!  PURPOSE
  !!   Reads radial table
  !!  INPUTS
  !!
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   T.Miyazaki
  !!  CREATION DATE
  !!
  !!  MODIFICATION HISTORY
  !!   10:17, 2004/01/12 dave
  !!    Changed formatting for reads to read new Siesta (1.3) pseudos
  !!   2017/11/10 14:24 dave
  !!    Changed yp1 and ypn to use first derivative (simple FD) instead of zero
  !!   2018/03/08 09:54 dave
  !!    Added consistency check between cutoff and step size/number of points
  !!  SOURCE
  !!
  subroutine radial_read_ascii(op,lun)

    use numbers, ONLY: BIG, zero, RD_ERR
    use splines, ONLY: spline
    use global_module, ONLY: iprint_pseudo
    use GenComms, ONLY: myid, cq_warn

    implicit none

    ! Passed variables
    type(rad_func),intent(out) :: op
    integer,intent(in)         :: lun

    ! Local variables
    character(len=80) :: sub_name = "radial_read_ascii"
    integer :: j, npts
    real(double) :: dummy, r0, r1
    real(double) :: yp1, ypn
    real(double) :: delta, cutoff

    !ori read(lun,'(i4,2g25.15)') npts, op%delta, op%cutoff
    !    read(lun,'(i4,2g25.15)') npts, delta, cutoff
    read(lun,*) npts, delta, cutoff
    ! DRB 2018/03/08 Adding consistency check between cutoff and delta
    if(abs(cutoff - delta*real(npts-1,double))>RD_ERR) then
       call cq_warn(sub_name, "In ion file, cutoff and step inconsistent, cutoff will be adjusted: ",&
            cutoff,delta*real(npts-1,double))
       cutoff=delta*real(npts-1,double)
    end if
    if(myid==0.AND.iprint_pseudo>3) write(io_lun,fmt='(10x,"Radius: ",f15.10)') cutoff
    op%delta = delta
    op%cutoff= cutoff
    call rad_alloc(op,npts)
    ! Test step size in radial grid is consistent
    read(lun,*) r0, op%f(1)
    read(lun,*) r1, op%f(2)
    if(abs(r1-r0-delta)>RD_ERR) then
       call cq_warn(sub_name, "Radial grid and step inconsistent ! ",r1-r0,delta)
       ! DRB 2018/03/08 10:55
       ! I can see an argument to abort here - if the step size and grid are inconsistent we may have problems
    end if
    do j=3,npts
       read(lun,*) dummy, op%f(j)
    enddo
    !ori call rad_setup_d2(op)
    ! conquest version
    !  yp1= BIG * 1.1
    !  ypn= BIG * 1.1
    yp1= (op%f(2)-op%f(1))/delta!zero
    ypn= (op%f(npts)-op%f(npts-1))/delta!zero
    call spline(op%n, op%delta, op%f, yp1, ypn, op%d2)
    return
  end subroutine radial_read_ascii
  !!***

  !---------------------------------------------------------------
  ! sbrt read_ion_ascii_tmp
  !---------------------------------------------------------------

  !!****f* pseudo_tm_info/read_ion_ascii_tmp *
  !!
  !!  NAME
  !!   read_ion_ascii_tmp
  !!  USAGE
  !!
  !!  PURPOSE
  !!   this subroutine only reads KB projectors,
  !!   chlocal and partial core charge
  !!  INPUTS
  !!
  !!
  !!  USES
  !!
  !!  AUTHOR
  !!   T.Miyazaki
  !!  CREATION DATE
  !!
  !!  MODIFICATION HISTORY
  !!   13:04, 31/10/2003 drb
  !!    Added routine to create lmax for pseudos
  !!   11:00, 13/02/2006 drb
  !!    Changed line length
  !!   2006/10/17 08:18 dave
  !!    Rewrote to read once on inode and broadcast
  !!   2007/02/12 08:21 dave
  !!    Added changes for TM potentials from ABINIT
  !!   2007/08/16 15:36 dave
  !!    More tweaks for local potential (particularly short/long range division)
  !!   09:30, 27/11/2007 drb
  !!    Added calculation of prefac on all processors for local potential
  !!   2008/05/25
  !!    Added timers
  !!   2008/09/01 08:21 dave
  !!    Added io_ routines from input_module
  !!   2013/07/05 dave
  !!    Added reading and copying of z, charge state of ion
  !!    Changed to read ALL radial functions for a given l (now picks up "semi-core" PAOs)
  !!   2016/01/28 16:46 dave
  !!    Updated module name to ion_electrostatic
  !!   2016/02/09 08:17 dave
  !!    Changed to use erfc from functions module
  !!   2016/01/09 20:00 nakata
  !!    Added pao_info%angmom%prncpl
  !!   2017/02/24 10:06 dave & tsuyoshi
  !!    Bug fix for read_header_tmp: removed formatting for zval
  !!   2017/03/13 dave
  !!    Repeated bug fix for other if branch in read_header_tmp (credit Jac van Driel for finding issue)
  !!   2017/03/14 dave
  !!    Changed string comparisons to use leqi from input_module (direct comparison breaks on Cray)
  !!   2017/03/23 dave
  !!    Added storage of spacing of PAO table
  !!   2017/11/10 14:22 dave
  !!    Bug fix: added Ry->Ha conversion for Siesta VNA d2 table (as well as f table)
  !!   2018/01/22 14:39 JST dave
  !!    Added test for lmax_pao/lmax_ps to find maximum PAO/PP angular momentum
  !!   2018/10/30 11:02 dave
  !!    Added semicore flag for each zeta
  !!   2018/11/02 16:30 nakata
  !!    Bug fix: set semicore when numprocs>1
  !!   2020/01/22 16:59 dave
  !!    Bug fix: change header to read Hamann code version line if present
  !!  SOURCE
  !!
  subroutine read_ion_ascii_tmp(ps_info,pao_info)

    use numbers, ONLY: six, half, zero, pi, one
    use global_module, ONLY: numprocs, iprint_pseudo
    use GenComms, ONLY: inode, ionode, gcopy, cq_warn, cq_abort
    use pao_format, ONLY: species_pao
    use splines, ONLY : spline
    use memory_module, ONLY: reg_alloc_mem, type_dbl
    use functions, ONLY: erfc_cq
    use input_module, ONLY: io_assign, io_close
    use pseudopotential_common, ONLY: pseudo_type, SIESTA, ABINIT
    use maxima_module, only: lmax_pao, lmax_ps

    implicit none

    ! Passed variables
    type(pseudo_info),intent(inout) :: ps_info
    type(species_pao),intent(inout) :: pao_info

    ! Local variables
    character(len=80) :: sub_name = "read_ion_ascii_tmp"
    type(rad_func), dimension(:), allocatable :: dummy_rada
    type(rad_func) :: dummy_rad

    character(len=80) :: filename
    integer :: i, lun , i1, i2, i3, i4, z
    real(double) :: dummy, a, r
    integer :: n_orbnl, n_pjnl
    real(double) :: zval, yp1, ypn, erfarg, tmpv
    real(double), parameter :: ln10 = 2.302585092994_double

    integer :: iproc, lmax, maxz, alls, nzeta, l, count,tzl, xc_func, ps_type
    real(double), allocatable :: thispop(:)
    integer, allocatable :: thisl(:), thisn(:), thisz(:), zl(:), indexlz(:,:)

    if(inode==ionode) then
       filename=ps_info%filename
       call io_assign(lun)
       open(lun,file=filename,status='old',form='formatted')
       rewind(lun)

       xc_func = 0
       ps_type = 0
       call read_header_tmp(n_orbnl,lmax,n_pjnl, zval, z, lun, xc_func, ps_type)
       if(xc_func/=0) then
          ps_info%functional = xc_func
       end if
       if(ps_type/=0) then
          ps_info%ps_type = ps_type
          if(ps_info%ps_type/=pseudo_type) &
               call cq_abort("Error: ion file pseudopotential type incompatible with input file ", &
               ps_info%ps_type, pseudo_type)
       else
          call cq_warn(sub_name, "Pseudopotential type not detected from ion file")
       end if
       call alloc_pseudo_info(ps_info, n_pjnl)
       call start_timer(tmr_std_allocation)
       allocate(dummy_rada(n_orbnl),thisl(n_orbnl),thisn(n_orbnl),thisz(n_orbnl),thispop(n_orbnl),zl(0:lmax))
       call stop_timer(tmr_std_allocation)

       ps_info%z = z
       ps_info%zval = zval

       read(lun,*)
       zl = 0
       maxz = 0
       if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"Reading PAOs")')
       do i=1,n_orbnl
          read(lun,*) i1,i2,i3,i4, dummy
          thisl(i)=i1
          thisn(i)=i2
          !thisz(i)=i3
          thispop(i)=dummy
          zl(thisl(i)) = zl(thisl(i))+1
          thisz(i)=zl(thisl(i))
          !if(thisz(i)>zl(thisl(i))) zl(thisl(i))=thisz(i)
          maxz = max(maxz,thisz(i))
          if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"l: ",i3," z: ",i3)') i1,i3
          call radial_read_ascii(dummy_rada(i),lun)
       enddo
       allocate(indexlz(maxz,0:lmax))
       indexlz = 0
       do i=1,n_orbnl
          indexlz(thisz(i),thisl(i)) = i
          if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"indexlz: ",3i5)') thisz(i),thisl(i),i
       end do
       ! Now store data
       pao_info%greatest_angmom = lmax
       call start_timer(tmr_std_allocation)
       allocate(pao_info%angmom(0:lmax),STAT=alls)
       call stop_timer(tmr_std_allocation)
       count = 0
       if(alls/=0) call cq_abort('Failed to allocate PAOs')
       if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"Storing PAOs lmax: ",i5)') lmax
       do l=0,lmax
          pao_info%angmom(l)%n_zeta_in_angmom = zl(l)
          if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"l, zl: ",2i5)') l,zl(l)
          if(zl(l)>0) then
             call start_timer(tmr_std_allocation)
             allocate(pao_info%angmom(l)%zeta(zl(l)),pao_info%angmom(l)%prncpl(zl(l)), &
                  pao_info%angmom(l)%occ(zl(l)),pao_info%angmom(l)%semicore(zl(l)),STAT=alls)
             pao_info%angmom(l)%semicore(:) = 0
             if(alls/=0) call cq_abort('Failed to allocate PAOs zeta')
             call stop_timer(tmr_std_allocation)
             count = count + zl(l)*(2*l+1)
             do nzeta = 1,zl(l)
                i = indexlz(nzeta,l)
                if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"i,z,l: ",3i5)') i,nzeta,l
                pao_info%angmom(l)%zeta(nzeta)%length = dummy_rada(i)%n
                pao_info%angmom(l)%zeta(nzeta)%cutoff = dummy_rada(i)%cutoff
                pao_info%angmom(l)%zeta(nzeta)%delta = dummy_rada(i)%delta
                pao_info%angmom(l)%prncpl(nzeta) = thisn(i)
                pao_info%angmom(l)%occ(nzeta) = thispop(i)
                if(pao_info%angmom(l)%zeta(nzeta)%length>=1) then
                   call start_timer(tmr_std_allocation)
                   allocate(pao_info%angmom(l)%zeta(nzeta)%table(pao_info%angmom(l)%zeta(nzeta)%length),STAT=alls)
                   if(alls/=0) call cq_abort('Failed to allocate PAOs zeta')
                   call reg_alloc_mem(area_pseudo,pao_info%angmom(l)%zeta(nzeta)%length,type_dbl)
                   call stop_timer(tmr_std_allocation)
                   pao_info%angmom(l)%zeta(nzeta)%table(1:dummy_rada(i)%n) = dummy_rada(i)%f(1:dummy_rada(i)%n)
                   call start_timer(tmr_std_allocation)
                   allocate(pao_info%angmom(l)%zeta(nzeta)%table2(pao_info%angmom(l)%zeta(nzeta)%length),STAT=alls)
                   if(alls/=0) call cq_abort('Failed to allocate PAOs zeta')
                   call reg_alloc_mem(area_pseudo,pao_info%angmom(l)%zeta(nzeta)%length,type_dbl)
                   call stop_timer(tmr_std_allocation)
                else
                   call cq_abort('PAO with zero length: ',l,nzeta)
                end if
                yp1 = (pao_info%angmom(l)%zeta(nzeta)%table(2)-pao_info%angmom(l)%zeta(nzeta)%table(1))/ &
                     dummy_rada(i)%delta
                ypn = (pao_info%angmom(l)%zeta(nzeta)%table(dummy_rada(i)%n)-&
                     pao_info%angmom(l)%zeta(nzeta)%table(dummy_rada(i)%n-1))/ &
                     dummy_rada(i)%delta
                call spline(dummy_rada(i)%n,dummy_rada(i)%delta,pao_info%angmom(l)%zeta(nzeta)%table, &
                     yp1,ypn,pao_info%angmom(l)%zeta(nzeta)%table2)
             end do
          end if
       end do
       pao_info%count = count
       deallocate(indexlz)
       do i=1,n_orbnl
          call rad_dealloc(dummy_rada(i))
       end do
       deallocate(dummy_rada,thisl,thisn,thisz,thispop,zl)

       ! KBs
       if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"Reading KB projectors ")')
       read(lun,*)
       ps_info%lmax = 0
       do i=1,ps_info%n_pjnl
          !          read(lun,'(2i3,f22.16)') &
          read(lun,*) &
               ps_info%pjnl_l(i), ps_info%pjnl_n(i), ps_info%pjnl_ekb(i)
          !!   25/Jul/2002 TM : Rydberg units -> Hartree units
          ps_info%pjnl_ekb(i) = half* ps_info%pjnl_ekb(i)
          if(iprint_pseudo>1.AND.inode==ionode) write(io_lun,fmt='(10x,"PSEUDO_TM: i, l,n, ekb = ",2i3,f22.16)') &
               ps_info%pjnl_l(i), ps_info%pjnl_n(i),ps_info%pjnl_ekb(i)
          if(ps_info%pjnl_l(i)>ps_info%lmax) ps_info%lmax = ps_info%pjnl_l(i)
          call radial_read_ascii(ps_info%pjnl(i),lun)
       enddo
       !Vlocal
       if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"Reading Vna ")')
       read(lun,*)
       ! Read and store local potential; it may be local part of pseudo, or neutral atom
       ! DRB 2017/02/21 switch here: after vna, get vlocal for Hamann, chlocal for siesta
       ! DRB 2017/02/21 This is the NA potential
       if(pseudo_type==SIESTA) then
          call radial_read_ascii(ps_info%vna,lun)
          ps_info%vna%f  = half*ps_info%vna%f  ! Siesta works in Ry
          ps_info%vna%d2 = half*ps_info%vna%d2 ! Siesta works in Ry
          if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"Reading Chlocal ")')
          read(lun,*)
          call radial_read_ascii(ps_info%chlocal,lun)
       else if(pseudo_type==ABINIT) then
          call radial_read_ascii(ps_info%vna,lun)
          if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"Reading Vlocal ")')
          read(lun,*)
          call radial_read_ascii(ps_info%vlocal,lun)
       end if
       !sets up gaussian which is used to calculate G=0 term
       !  exp( - alpha * cutoff **2 ) = 10 **(-six)
       if(ps_info%tm_loc_pot==loc_pot) then
          ps_info%alpha = six * ln10 / (ps_info%vlocal%cutoff ** 2)
          ps_info%prefac = (ps_info%alpha/pi)**1.5_double
          ! Remove long-range part from short-range part
          a = sqrt(ps_info%alpha)
          !r = ps_info%vlocal%delta
          r = 1.0e-7_double
          erfarg = one - erfc_cq(a*r)
          !erfarg = derf(a*r)
          tmpv = ps_info%vlocal%f(1)
          ps_info%vlocal%f(1) = ps_info%vlocal%f(1) + ps_info%zval * erfarg/r
          !write(52,*) r,ps_info%vlocal%f(1),tmpv,ps_info%zval * erfarg/r!ps_info%zval * derf(a*r)/r
          do i=2,ps_info%vlocal%n
             r = ps_info%vlocal%delta*real(i-1,double)
             erfarg = one - erfc_cq(a*r)
             tmpv = ps_info%vlocal%f(i)
             ps_info%vlocal%f(i) = ps_info%vlocal%f(i) + ps_info%zval * erfarg/r
             !write(52,*) r,ps_info%vlocal%f(i), tmpv,ps_info%zval * erfarg/r!,ps_info%zval * derf(a*r)/r
          end do
          call spline(ps_info%vlocal%n, ps_info%vlocal%delta, ps_info%vlocal%f, zero, zero, ps_info%vlocal%d2)
       else
          ps_info%alpha = six * ln10 / (ps_info%chlocal%cutoff ** 2)
       endif
       !Core
       ! 2017/03/17 dave
       ! I'm removing this because we read the PCC flag from the pseudopotential header
       !ps_info%flag_pcc = .false.
       !! I *really* don't like this, but we'll stay with it for now
       !read(lun,*,end=9999)
       if(ps_info%flag_pcc) then
          read(lun,*)
          if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x,"Reading pcc ")')
          call radial_read_ascii(ps_info%chpcc,lun)
       end if
       !ps_info%flag_pcc = .true.

       !9999   continue
       call io_close(lun)
    endif !  (inode == ionode) then
    ! Now broadcast the information
    call gcopy(lmax   )
    if(lmax>lmax_pao) lmax_pao = lmax
    call gcopy(n_pjnl )
    call gcopy(zval)
    call gcopy(xc_func)
    call gcopy(z)
    if(inode/=ionode) then
       call alloc_pseudo_info(ps_info, n_pjnl)
       ps_info%zval = zval
       ps_info%z = z
       if(xc_func/=0) ps_info%functional = xc_func
       pao_info%greatest_angmom = lmax
       allocate(pao_info%angmom(0:lmax),STAT=alls)
    end if
    if(numprocs>1) then
       count = 0
       if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,*) 'PAOs'
       do l=0,lmax
          if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(10x," l is ",i5)') l
          call gcopy(pao_info%angmom(l)%n_zeta_in_angmom)
          tzl = pao_info%angmom(l)%n_zeta_in_angmom
          if(tzl>0) then
             if(inode/=ionode) then
                call start_timer(tmr_std_allocation)
                allocate(pao_info%angmom(l)%zeta(tzl),pao_info%angmom(l)%prncpl(tzl),&
                     pao_info%angmom(l)%occ(tzl),pao_info%angmom(l)%semicore(tzl),STAT=alls)
                if(alls/=0) call cq_abort('Failed to allocate PAOs zeta')
                call stop_timer(tmr_std_allocation)
             end if
             count = count + tzl*(2*l+1)
             do nzeta = 1,tzl
                call gcopy(pao_info%angmom(l)%zeta(nzeta)%length)
                call gcopy(pao_info%angmom(l)%zeta(nzeta)%cutoff)
                call gcopy(pao_info%angmom(l)%zeta(nzeta)%delta)
                call gcopy(pao_info%angmom(l)%prncpl(nzeta))
                call gcopy(pao_info%angmom(l)%occ(nzeta))
                call gcopy(pao_info%angmom(l)%semicore(nzeta))
                if(inode/=ionode) then
                   if(pao_info%angmom(l)%zeta(nzeta)%length>=1) then
                      call start_timer(tmr_std_allocation)
                      allocate(pao_info%angmom(l)%zeta(nzeta)%table(pao_info%angmom(l)%zeta(nzeta)%length),STAT=alls)
                      if(alls/=0) call cq_abort('Failed to allocate PAOs zeta',pao_info%angmom(l)%zeta(nzeta)%length)
                      call reg_alloc_mem(area_pseudo,pao_info%angmom(l)%zeta(nzeta)%length,type_dbl)

                      allocate(pao_info%angmom(l)%zeta(nzeta)%table2(pao_info%angmom(l)%zeta(nzeta)%length),STAT=alls)
                      if(alls/=0) call cq_abort('Failed to allocate PAOs zeta ',pao_info%angmom(l)%zeta(nzeta)%length)
                      call reg_alloc_mem(area_pseudo,pao_info%angmom(l)%zeta(nzeta)%length,type_dbl)
                      call stop_timer(tmr_std_allocation)
                   end if
                end if
                call gcopy(pao_info%angmom(l)%zeta(nzeta)%table,pao_info%angmom(l)%zeta(nzeta)%length)
                call gcopy(pao_info%angmom(l)%zeta(nzeta)%table2,pao_info%angmom(l)%zeta(nzeta)%length)
             end do
          end if
       end do
       pao_info%count = count
       if(inode/=ionode) ps_info%lmax = 0
       if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,*) 'PPs'
       do i=1,ps_info%n_pjnl
          if(iprint_pseudo>3.AND.inode==ionode) write(io_lun,fmt='(2x,"Projector number: ",i3)') i
          call gcopy(ps_info%pjnl_l(i))
          call gcopy(ps_info%pjnl_n(i))
          call gcopy(ps_info%pjnl_ekb(i))
          if(inode/=ionode.AND.ps_info%pjnl_l(i)>ps_info%lmax) ps_info%lmax = ps_info%pjnl_l(i)
          call gcopy(ps_info%pjnl(i)%delta)
          call gcopy(ps_info%pjnl(i)%cutoff)
          call gcopy(ps_info%pjnl(i)%n)
          if(inode/=ionode) call rad_alloc(ps_info%pjnl(i),ps_info%pjnl(i)%n)
          call gcopy(ps_info%pjnl(i)%f,ps_info%pjnl(i)%n)
          call gcopy(ps_info%pjnl(i)%d2,ps_info%pjnl(i)%n)
       end do
       if(ps_info%lmax>lmax_ps) lmax_ps = ps_info%lmax
       if(pseudo_type==SIESTA) then
          !Chlocal
          call gcopy(ps_info%tm_loc_pot)
          call gcopy(ps_info%chlocal%n)
          call gcopy(ps_info%chlocal%cutoff)
          call gcopy(ps_info%chlocal%delta)
          if(inode/=ionode) call rad_alloc(ps_info%chlocal,ps_info%chlocal%n)
          call gcopy(ps_info%chlocal%f,ps_info%chlocal%n)
          call gcopy(ps_info%chlocal%d2,ps_info%chlocal%n)
          ! Vna
          call gcopy(ps_info%vna%n)
          call gcopy(ps_info%vna%cutoff)
          call gcopy(ps_info%vna%delta)
          if(inode/=ionode) call rad_alloc(ps_info%vna,ps_info%vna%n)
          call gcopy(ps_info%vna%f,ps_info%vna%n)
          call gcopy(ps_info%vna%d2,ps_info%vna%n)
       else if(pseudo_type==ABINIT) then
          !Vlocal
          call gcopy(ps_info%tm_loc_pot)
          call gcopy(ps_info%vlocal%n)
          call gcopy(ps_info%vlocal%cutoff)
          call gcopy(ps_info%vlocal%delta)
          if(inode/=ionode) call rad_alloc(ps_info%vlocal,ps_info%vlocal%n)
          call gcopy(ps_info%vlocal%f,ps_info%vlocal%n)
          call gcopy(ps_info%vlocal%d2,ps_info%vlocal%n)
          ! Vna
          call gcopy(ps_info%vna%n)
          call gcopy(ps_info%vna%cutoff)
          call gcopy(ps_info%vna%delta)
          if(inode/=ionode) call rad_alloc(ps_info%vna,ps_info%vna%n)
          call gcopy(ps_info%vna%f,ps_info%vna%n)
          call gcopy(ps_info%vna%d2,ps_info%vna%n)
       end if
       if(inode/=ionode) then
          if(ps_info%tm_loc_pot==loc_pot) then
             ps_info%alpha = six * ln10 / (ps_info%vlocal%cutoff ** 2)
             ps_info%prefac = (ps_info%alpha/pi)**1.5_double
          else
             ps_info%alpha = six * ln10 / (ps_info%chlocal%cutoff ** 2)
          end if
       end if
       call gcopy(ps_info%flag_pcc)
       if(ps_info%flag_pcc) then
          call gcopy(ps_info%chpcc%n)
          call gcopy(ps_info%chpcc%cutoff)
          call gcopy(ps_info%chpcc%delta)
          if(inode/=ionode) call rad_alloc(ps_info%chpcc,ps_info%chpcc%n)
          call gcopy(ps_info%chpcc%f,ps_info%chpcc%n)
          call gcopy(ps_info%chpcc%d2,ps_info%chpcc%n)
       end if
    end if ! numprocs>1
  contains

    subroutine read_header_tmp(n_orbnl, lmax_basis, n_pjnl, zval, z, unit, xc_func, pseudo_type)

      use global_module, ONLY: iprint_pseudo
      use input_module, ONLY: leqi
      use XC, ONLY: functional_lda_pw92, functional_gga_pbe96, functional_gga_pbe96_rev98, &
           functional_gga_pbe96_r99, functional_gga_pbe96_wc

      implicit none

      integer, intent(in)         :: unit
      integer, intent(out) :: n_orbnl, n_pjnl, z, xc_func, pseudo_type
      real(double), intent(out) :: zval

      character(len=78) :: line, trim_line

      character(len=2)  :: symbol
      character(len=20) :: label
      integer :: lmax_basis, lmax_projs, zval_int
      real(double) :: mass, self_energy

      ! For judging P.C.C.
      character(len=2) :: symbol2, xc_code
      character(len=3) :: rel
      character(len=4) :: pcc

      xc_func = 0
      pseudo_type = 0
      ! Judge if P.C.C. is considered
      read(unit,'(a)') line
      trim_line = trim(line)
      if (leqi(trim_line(1:10),'<preamble>')) then
         read(unit,'(a)') line
         trim_line = trim(line)
         do while(.NOT.leqi(trim_line(1:11),'</preamble>'))
            read(unit,'(a)') line
            trim_line = trim(line)
            if (leqi(trim_line(1:24),'<pseudopotential_header>')) then
               read (unit, '(1x,a2,1x,a2,1x,a3,1x,a4)') symbol2, xc_code, rel, pcc
               if (leqi(pcc,'pcec')) then
                  ps_info%flag_pcc = .true.
               else !if (pcc .NE. 'pcec') then
                  ps_info%flag_pcc = .false.
               endif
               ! Set XC functional
               if(xc_func==0) then
                  if(leqi(xc_code,'ca')) then      ! LDA; ATOM code uses PZ81 but use PW92
                     xc_func = functional_lda_pw92
                  else if(leqi(xc_code,'pb')) then ! PBE
                     xc_func = functional_gga_pbe96
                  else if(leqi(xc_code,'rv')) then ! RevPBE
                     xc_func = functional_gga_pbe96_rev98
                  else if(leqi(xc_code,'rp')) then ! RPBE
                     xc_func = functional_gga_pbe96_r99
                  else if(leqi(xc_code,'wc')) then ! Wu-Cohen
                     xc_func = functional_gga_pbe96_wc
                  end if
               end if
               read(unit,'(a)') line
               trim_line = trim(line)
               if(.NOT.leqi(trim_line(1:25),'</pseudopotential_header>')) then
                  if(leqi(trim_line(2:4),'ATM')) pseudo_type = 1 ! Siesta-type
               end if
            else if (leqi(trim_line(1:14),'<Conquest_pseu')) then
               read(unit,'(a)') line ! Check this for Hamann
               if(leqi(line(3:8),'Hamann')) pseudo_type = 3
               if(leqi(line(3:10),'Hamann c')) read(unit,'(a)') line
               read(unit,'(a)') line ! Core radii
               read(unit,'(a)') line ! Valence shells
               read(unit,'(a26,i7)') line, xc_func
            endif
         end do
      endif

      if(iprint_pseudo>1) then
         read(unit,'(a2)') symbol
         write(io_lun,fmt='(10x,"symbol = ",a2)') symbol
         read(unit,'(a20)') label
         write(io_lun,fmt='(10x,"label = ",a20)') label
         read(unit,*) z
         write(io_lun,fmt='(10x,"z = ",i5)') z
         read(unit,*) zval
         write(io_lun,fmt='(10x,"zval_int, zval = ",f12.5)') zval
         read(unit,*) mass
         write(io_lun,fmt='(10x,"mass = ",g25.15)') mass
         read(unit,*) self_energy
         write(io_lun,fmt='(10x,"self_energy = ",g25.15)') self_energy
         read(unit,'(2i4)') lmax_basis, n_orbnl
         write(io_lun,fmt='(10x,"lmax_basis, n_orbnl = ",2i4)') lmax_basis, n_orbnl
         read(unit,'(2i4)') lmax_projs, n_pjnl
         write(io_lun,fmt='(10x,"lmax_projs, n_pjnl = ",2i4)') lmax_projs, n_pjnl
         if(xc_func/=0) write(io_lun,fmt='(10x,"XC functional code = ",i8)') xc_func
      else
         read(unit,'(a2)') symbol
         read(unit,'(a20)') label
         read(unit,*) z
         read(unit,*) zval!_int
         !zval=zval_int
         read(unit,'(g25.15)') mass
         read(unit,'(g25.15)') self_energy
         read(unit,'(2i4)') lmax_basis, n_orbnl
         read(unit,'(2i4)') lmax_projs, n_pjnl
      end if
    end subroutine read_header_tmp
  end subroutine read_ion_ascii_tmp

  !!***
end module pseudo_tm_info

max_stars_repo_name: ferroteam/CONQUEST-release

max_issues_repo_name: ferroteam/CONQUEST-release

max_forks_repo_name: ferroteam/CONQUEST-release


text:
C$Procedure SPAT2B ( SPACIT, convert transfer to binary )

      SUBROUTINE SPAT2B ( LOGFIL, LOGLUN )

C$ Abstract
C
C     SPACIT utility subroutine used to manage the conversion of SPICE
C     transfer format files into their equivalent binary formats. This
C     subroutine is for use only be the SPACIT program. Use it at your
C     own risk.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     None.
C
C$ Declarations

      IMPLICIT NONE

      LOGICAL               LOGFIL
      INTEGER               LOGLUN

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C      LOGFIL    I    Logical flag indicating a log file is being kept.
C      LOGLUN    I    The logical unit of the log file.
C
C$ Detailed_Input
C
C      LOGFIL   Logical flag indicating a log file is being kept. This
C               Variable has the value of .TRUE. if a log file is being
C               written, and a value of .FALSE. otherwise.
C
C      LOGLUN   The logical unit of the log file. If LOGFIL has the
C               value .TRUE. then LOGLUN will be the Fortran logical
C               unit of the log file.
C
C$ Detailed_Output
C
C     None.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     None.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     xxx
C
C$ Examples
C
C     xxx
C
C$ Restrictions
C
C     xxx
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     K.R. Gehringer     (JPL)
C
C$ Version
C
C
C-    SPICE Version 2.0.0, 13-MAR-1999 (EDW)
C
C        Replaced WRITLN's to standard output with TOSTDO,
C        calls to GETFAT with ZZGETFAT, calls to CONVTB to
C        ZZCONVTB.
C
C        Added a condition to the ARCH-TYPE test block following
C        the call to ZZGETFAT.  The block test for both ARCH and
C        TYPE as unknows, then for ARCH as XFR, then either ARCH
C        or TYPE as unknown.
C
C-    Beta Version 1.0.0, 13-NOV-1994 (KRG)
C
C-&

C$ Index_Entries
C
C     spacit convert transfer to binary
C
C-&

C
C     SPICELIB functions
C
      INTEGER               CARDI

      LOGICAL               EXISTS
      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local parameters
C
C     Lower bound for a SPICELIB CELL data structure.
C
      INTEGER               LBCELL
      PARAMETER           ( LBCELL = -5 )
C
C     Maximum number of open binary files allowed.
C
      INTEGER               MAXOPN
      PARAMETER           ( MAXOPN = 1 )
C
C     Set a value for the logical unit which represents the standard
C     output device, commonly a terminal. A value of 6 is widely used,
C     but the Fortran standard does not specify a value, so it may be
C     different for different Fortran implementations.
C
      INTEGER               STDOUT
      PARAMETER           ( STDOUT = 6 )
C
C     Set a value for the length of an input text line.
C
      INTEGER               LINLEN
      PARAMETER           ( LINLEN = 255 )
C
C     Set a value for the length of a filename.
C
      INTEGER               FNMLEN
      PARAMETER           ( FNMLEN = 128 )
C
C     Set a length for the prompt.
C
      INTEGER               PRMLEN
      PARAMETER           ( PRMLEN = 80 )
C
C     Set a length for the status of a file: 'OLD' or 'NEW'.
C
      INTEGER               STALEN
      PARAMETER           ( STALEN = 3 )
C
C     Set the length for the type of a file.
C
      INTEGER               TYPLEN
      PARAMETER           ( TYPLEN = 4 )
C
C     Set the length for the architecture of a file.
C
      INTEGER               ARCLEN
      PARAMETER           ( ARCLEN = 3 )
C
C     Local variables
C
      CHARACTER*(ARCLEN)    ARCH
      CHARACTER*(FNMLEN)    BINFNM
      CHARACTER*(STALEN)    BFSTAT
      CHARACTER*(LINLEN)    LINE
      CHARACTER*(LINLEN)    MESSAG
      CHARACTER*(PRMLEN)    PRMPT
      CHARACTER*(FNMLEN)    XFRFNM
      CHARACTER*(STALEN)    XFSTAT
      CHARACTER*(TYPLEN)    TYPE

      INTEGER               I
      INTEGER               NUMOPN
      INTEGER               OPNSET(LBCELL:MAXOPN)
      INTEGER               NUMBER

      LOGICAL               FILEOK

C
C     Save everything to keep configuration control happy.
C
      SAVE

C
C     Initial values
C
      DATA MESSAG / 'Converting from transfer file to binary file.' /

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'SPAT2B' )
      END IF
C
C
C
      BFSTAT = 'NEW'
      XFSTAT = 'OLD'

      FILEOK = .FALSE.
      PRMPT  = '   Transfer file: '

      CALL GETFNM_1 ( PRMPT, XFSTAT, XFRFNM, FILEOK )

      IF ( FAILED() ) THEN
         CALL CHKOUT ( 'SPAT2B' )
         RETURN
      END IF

C
C     Get the architecture, file type and logical unit number for
C     the transfer file.
C
      CALL ZZGETFAT ( XFRFNM, ARCH, TYPE, NUMBER )

C
C     If something went wrong, close the transfer file and checkout.
C
      IF ( FAILED() ) THEN
         CLOSE ( NUMBER )
         CALL CHKOUT ( 'SPAT2B' )
         RETURN
      END IF

C
C     The file should be a transfer file.  If not, flag an error
C     then exit.
C
      IF ( ( ARCH .EQ. '?' ) .AND. ( TYPE .EQ. '?' ) ) THEN

         CALL SETMSG ( 'The architecture and type of the'
     .   //            ' file ''#'' could not be'
     .   //            ' determined.'                      )
         CALL ERRCH  ( '#', XFRFNM                         )
         CALL SIGERR ( 'SPICE(BADFILEFORMAT)'              )
         CALL CHKOUT ( 'SPAT2B' )
         RETURN

      ELSE IF ( (ARCH .NE. 'XFR') .AND. (ARCH .NE. 'DEC') ) THEN

         CALL SETMSG ( 'The file ''#'' was not a transfer'
     .   //            ' format SPICE file. In order to'
     .   //            ' convert a transfer file to binary'
     .   //            ' format, it must be a CK, EK, PCK,'
     .   //            ' or SPK file in one of the'
     .   //            ' supported transfer file formats.' )
         CALL ERRCH  ( '#', XFRFNM                         )
         CALL SIGERR ( 'SPICE(IMPROPERFILE)'               )
         CALL CHKOUT ( 'SPAT2B' )
         RETURN

      ELSE IF ( ( ARCH .EQ. '?' ) .OR. ( TYPE .EQ. '?' ) ) THEN

         CALL SETMSG ( 'The architecture or type of the'
     .   //            ' file ''#'' could not be'
     .   //            ' determined.  The file did not'
     .   //            ' read as a transfer file'          )
         CALL ERRCH  ( '#', XFRFNM                         )
         CALL SIGERR ( 'SPICE(BADFILEFORMAT)'              )
         CALL CHKOUT ( 'SPAT2B' )
         RETURN

      END IF

C
C     No get the binary file name.
C
      FILEOK = .FALSE.
      PRMPT  = '   Binary file  : '

      CALL GETFNM_1 ( PRMPT, BFSTAT, BINFNM, FILEOK )

      IF ( .NOT. FILEOK ) THEN
         CALL CHKOUT ( 'SPAT2B' )
         RETURN
      END IF

      LINE = '   #'
      CALL REPMC ( LINE, '#', MESSAG, LINE )
      CALL TOSTDO ( ' ' )
      CALL TOSTDO (  LINE )
      IF ( LOGFIL ) THEN
         CALL WRITLN ( ' ',  LOGLUN )
         CALL WRITLN ( LINE, LOGLUN )
      END IF

      LINE = '   Converting Transfer File: #'
      CALL REPMC  ( LINE, '#', XFRFNM, LINE )
      CALL TOSTDO ( ' ')
      CALL TOSTDO ( LINE )
      IF ( LOGFIL ) THEN
         CALL WRITLN ( ' ',  LOGLUN )
         CALL WRITLN ( LINE, LOGLUN )
      END IF

      LINE = '   To Binary File          : #'
      CALL REPMC  ( LINE, '#', BINFNM, LINE )
      CALL TOSTDO ( LINE )
      CALL TOSTDO ( ' ')
      IF ( LOGFIL ) THEN
         CALL WRITLN ( LINE, LOGLUN )
         CALL WRITLN ( ' ',  LOGLUN )
      END IF

      CALL TOSTDO ( '   Please wait...' )

C
C     Now call the converted.  The real work starts here.
C
      CALL ZZCONVTB ( XFRFNM, ARCH, TYPE, BINFNM, NUMBER )

C
C      Important question, did the conversion work?
C
      IF ( FAILED() ) THEN

C
C        Something failed.  Reset the error system.
C
         CALL RESET

         CALL SCARDI ( 0, OPNSET         )
         CALL CLEARI ( MAXOPN, OPNSET(1) )
         CALL DAFHOF ( OPNSET )
         NUMOPN = CARDI ( OPNSET )

C
C        Close the DAF files.
C
         IF ( NUMOPN .GT. 0 ) THEN
            DO I = 1, NUMOPN
               CALL DAFCLS ( OPNSET(I) )
            END DO
         END IF

         CALL SCARDI ( 0, OPNSET )
         CALL CLEARI ( MAXOPN, OPNSET(1) )
         CALL DASHOF ( OPNSET )

C
C        Close the DAS files.
C
         NUMOPN = CARDI ( OPNSET )
         IF ( NUMOPN .GT. 0 ) THEN
            DO I = 1, NUMOPN
               CALL DASCLS ( OPNSET(I) )
            END DO
         END IF

C
C        Erase any partially converted binary file, then checkout.
C
         IF ( EXISTS( BINFNM ) ) THEN
            CALL DELFIL( BINFNM )
         END IF

         CALL CHKOUT ( 'SPAT2B' )
         RETURN

      END IF

C
C     If we got this far, the conversion worked (yippie!).
C
      LINE = '   Binary file ''#'' created.'
      CALL REPMC ( LINE, '#', BINFNM, LINE )
      CALL TOSTDO ( ' ' )
      CALL TOSTDO ( LINE)
      CALL TOSTDO ( ' ' )
      IF ( LOGFIL ) THEN
         CALL WRITLN ( ' ',  LOGLUN )
         CALL WRITLN ( LINE, LOGLUN )
         CALL WRITLN ( ' ',  LOGLUN )
      END IF

      CALL CHKOUT ( 'SPAT2B' )
      RETURN
      END

ext: f

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
module burner_module

  use bl_types
  use bl_constants_module
  use bl_error_module
  use network

contains

  subroutine burner(dens, temp, Xin, dt, time, Xout)

    implicit none

    real(kind=dp_t), intent(in) :: dens, temp, Xin(nspec), dt, time
    real(kind=dp_t), intent(out) :: Xout(nspec)

    integer :: n
    real(kind=dp_t) :: dX

    Xout(:) = Xin(:)

  end subroutine burner

end module burner_module

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib